李慶昕

(遼寧省沈陽水文局，遼寧 沈陽 110094)

水文學者針對P-Ⅲ型分布函數的求解做了大量的研究和分析，并取得了一定的成果，如李世才等[1]通過利用伽馬分布為數和函數可實現Kp值的轉化求解，并基于此給出了伽馬分布通用算法[2]的截斷誤差和分析的表達式。吳明官等[3]提出了一種計算速度大于麥克勞林法、分布積分和變量代換法的新的切比雪夫快速算法。劉俊哲等[4]利用龍貝格積分法計算快速、穩(wěn)定性好和便于操作的優(yōu)點，并結合不完全伽馬分布函數的分布積分特征提出了一種計算速度高于自適應辛普森法和切比雪夫多項式法的計算方法；任柏幟等[5]基于分布函數展開式的基本形式提出了連分式展開式和P-Ⅲ型函數的通用計算法；劉仕平等[6]通過選取誤差作為輸入參數提出了伽馬分布函數的通用表達式。對于適應性不長的參數利用該方法其計算結果較易溢出。王文川等[7]通過對水文頻率曲線參數進行尋優(yōu)計算和結果擬合估計驗證了該方法的適用性和可靠性。

對皮爾遜Ⅲ型分布最優(yōu)的參數值通常采用建立目標函數值，并進行參數尋優(yōu)方法進行求解。橫坐標離差和是優(yōu)化算法的適線準則，對最優(yōu)的參數值可通過尋優(yōu)計算進行求解。本研究對P-Ⅲ型曲線的參數分別利用自適應誤差分析法、龍貝格、辛普森以及梯形算法進行尋優(yōu)計算，并進行分析研究，以期為提高水文頻率曲線參數尋優(yōu)計算運算速度和計算精度提供一定的理論支持。

1 基本理論

1.1 理論頻率

利用下式表征皮爾遜-Ⅲ型分布特征。

(1)

為便于計算需進行格式統一和積分換元計算，分別引入t=β(x-a0)，dx=dt/β，若x=xp，則滿足t=β(xp-a0)=u，代入上述計算公式可有：

(2)

1.2 經驗頻率

對于一般洪水可利用下述公式進行洪水頻率求解。

(3)

式中，m、n—由大到小對水文變量排序并按自然數排序變出的序號和樣本容量。

對于特大洪水可利用下述公式進行經驗頻率的計算。

(4)

其余量(n-l)在實測系列中的經驗頻率可利用下式進行計算。

(5)

1.3 適線準則

對水文序列值所對應的頻率可利用上述伽馬函數的數值積分法進行求解并以此可構建目標函數，表達式為：

(6)

式中，n、Wm—樣本容量和權重系數；pm、pi—樣本經驗頻率值和樣本m的理論頻率值；k—冪級此，本研究取k值取2。

2 數值積分法

自適應誤差積分法、梯形法、龍貝格法以及辛普森法等是數值積分法的主要方法。

2.1 梯形算法

(7)

利用量算對式(7)中xi=iΔx;i=0，1，2，…，n進行求解，且式(7)成為復合梯形公式。

2.2 辛普森算法

通過將積分區(qū)間劃分為若干偶數個子區(qū)間，并引入Δx=[b-a]/2n作為子區(qū)間步長計算步長，點y0，y1，y2在曲線上對應的坐標為x0，x1，x2，由此，利用曲線上三個相鄰的作拋物線并對曲線進行擬合求解，將經過x0，x1，x2三個點形成的曲線段利用拋物線進行替代，利用牛頓-萊布尼茨公式對拋物線進行計算求解，其面積積分求解公式如下。

(8)

式中，S0—[a，b]區(qū)間兩端點的縱坐標之和；S0、S1—奇數項縱坐標和偶數項縱坐標之和；n—偶數。

當n為偶數時可對[a，b]區(qū)間在函數f1(x)和f2(x)所夾面積利用辛普森公式進行求解，其計算公式如下。

(9)

利用辛普森法是對弧線利用拋物線進行替代進而能夠更好的接近曲線的彎曲程度并減低積分誤差。

2.3 Romberg算法

(10)

其計算步驟和相應公式分別如下所示：

(2)對曲線梯形值進行求解，其表達式為：

2.4 自適應誤差積分法

如果α≠2條件成立，對實際步長可利用泰勒級數展開法進行近似求解，計算公式如下:

hp=C0+C1(t-α)+C2(t-α)

(11)

式中，α、t—參數和自變量；C0、C1、C2—待定參數，當α取值不同時其計算方法不同。

若α>2，則C0、C1、C2計算公式分別如下：

若α<2，則C0、C1、C2計算公式分別如下：

利用步長計算公式并分別代入C0、C1、C2值可對步長進行求解。

若α=2則可利用基本步長公式對函數曲線進行面積積分計算，并選取

(12)

當運行計算時若α>>180則數據溢出并停止計算，利用標準正態(tài)分布近似表征伽馬函數分布的狀態(tài)，即呈N(0，1)=φ(u)，可利用下述公式代表標準正態(tài)分布。

(13)

3 尋優(yōu)建模

SSO優(yōu)化算法的計算步驟：設定參數可進行搜索尋優(yōu)的空間維度為n，雌性和雄性蜘蛛的個體數分別為Nf和Nm，則整個蜘蛛群體內蜘蛛個體數目總量為N，且Nf和Nm滿足下式：

Nf=floor[(0.9-rand×0.25)N];
N=Nm+Nf

(14)

式中，rand—0～1之間的隨機數字；floor—實數和整數之間的映射關系。

設定S蜘蛛種群中有N個蜘蛛個體，初始雌性蜘蛛隨機樣本F={f1，f2，…，fNf}，初始雄性蜘蛛隨機樣本F={m1，m2，…，mNm}，通過下述計算公式設定交配范圍r。蜘蛛群體S={s1=f1，s2=f2，…，sNf=fNf，sNf+1=m1，sNf+2=m2，…，sN=mNm}。

(15)

群體內各蜘蛛的自身重量的計算如下：

(16)

式中，J(si)—通過計算所得到蜘蛛個體的目標函數適應值；bests—max(J(si)worsts=minJ(si))。

蜘蛛群體中根據蜘蛛個體的重量進行交配概率的確定，蜘蛛個體的重量越大則其繁衍后代的概率越大，并采用輪盤賭法對教派概率psi進行確定，計算公式如下：

(17)

通過上述計算過程，可分別對不同的參數組合下的參數進行優(yōu)化判定，分別對龍貝格算法、紫志英積分法、梯形算法以及辛普森法的各個序列值進行理論頻率值的計算；若滿足收斂性則計算終止，否則重新返回并計算，直至滿足要求。

4 數據檢驗

4.1 理想數據檢驗

本研究結合劉光文等相關研究成果，依據理想數據系列相關方法和理論進行統計試驗分析，其中判別標準為x0.01%與x0.1%之間的相對誤差低于12%～16.8%為基本依據。對曲線的參數分別利用上述方法進行數值積分求解，并將已知序列與理想樣本進行對比分析，計算對比結果見表1。

表1 理想數據參數計算結果

表2 在不同水文序列中各計算方法的計算結果值

4.2 討論分析

P-Ⅲ型曲線利用梯形積分算法的x0.1%決定平均誤差值為-16.62%，其中在合格區(qū)間范圍內的組數為15組，17組數據的誤差值最大為-31.68%；而x0.01%的誤差平均值為-17.36%，在合格誤差范圍內的有12組，其中誤差最大值為34.62%。數值設計值整體處于較低水平，誤差評價值處于控制水平以外，由此表明，利用梯形法求解值存在較大誤差，即擬合結果不符合設計精度要求。利用辛普森法的擬合結果顯示，x0.1%和x0.01%水平的誤差平均值分別為17.45%和12.26%，處于誤差合格控制范圍內的組數分別為20組和16組，二者的最大誤差值分別為-16.62%和-23.05%；利用辛普森算法的擬合結果值整體處于較低水平，且各擬合結果誤差處于合格范圍，相對于梯形算法利用辛普森法的擬合結果基本能夠符合設計有關精度要求。利用龍貝格法的擬合結果顯示，x0.1%和x0.01%水平的誤差平均值分別為13.18%和11.61%，處于誤差合格控制范圍內的組數分別為18組和15組，二者的最大誤差值分別為-16.64%和-17.49%；相對理論假象值設計值整體較低且誤差范圍處于合理控制范圍內，由此表明利用龍貝格法求解P-Ⅲ型積分曲線整體可以符合相關設計精度。利用自適應誤差積分法的擬合結果顯示，x0.1%和x0.01%水平的誤差平均值分別為-8.92%和-12.18%，處于誤差合格控制范圍內的組數分別為20組和20組，二者的最大誤差值分別為-12.66%和-10.44%%；利用自適應誤差積分法求解的P-Ⅲ型積分曲線低于理想假象值，由此表明利用該方法進行擬合的結果表現出良好的精度與可靠性并符合有關設計精度要求[11]。

5 實測數據

利用相關文獻[12]和資料中對水文序列的計算方法和理論可對曲線參數進行擬合估計，計算結果見表2，表中洪水序列1和2分別代表一般洪水和特大洪水序列。

由表2參數估計結果可知，利用文中所述4中方法的參數估計結果大致相同，且根據曲線變化趨勢可以發(fā)現，原始序列利用文中所述四類方法均能較好的進行擬合，由此表明，目標函數適線法利用橫坐標離差平方和具有較好的適用性和可靠性。另一方面結合離差平方和統計結果可知，自適應誤差積分法的離差平方和在一般洪水序列和特大洪水序列擬合中值較小，而采用其他三種方法的離差平方和相差較大并與理想數據統計值存在一定的差異；龍貝格和辛普森算法的擬合效果低于梯形算法。

由表2統計結果運算時間可知，相對于其他算法，自適應誤差分析法運算速度明顯較快，且積分步長與精度和積分區(qū)間存在密切的關系；在相同步長時梯形法僅僅是對積分區(qū)間求和而不表現出辛普森法的擬合特征，由于龍貝格法的運算量和計算時間隨迭代過程的減少而加快。綜上所述，利用自適應誤差積分法的離差平方和最小且運算速度較快，其擬合結果表現出較好的精度和可靠性。自適應誤差積分法在適線效果和離差平方和方面均優(yōu)于其他數值積分法。

6 結論

為進一步提高數值積分法的運算速度和計算精度，本文對水文頻率參數分別用4種數值積分法進行求解結算，并通過一般洪水、特大洪水序列的實測值與理想數據檢驗分析，探討了各數值積分法的適用性和可靠性，得出的主要結論如下：

(1)分割點的大小是決定辛普森和梯形算法的數值積分精度的關鍵因素；自適應誤差積分法可對擬合誤差精度結合實際需要進行調整，通過泰勒展開式對步長進行表示，可有效節(jié)約運算時間，并以此提高了運算的計算精度。

(2)對數值積分的步長，可利用誤差值的大小進行求解，并達到擬合結果的最優(yōu)；而利用龍貝格和辛普森法二者的擬合精度基本相同，均能符合相關設計要求，設計值擬合結果整體保持一致；利用梯形法的擬合結果表現出明顯的端距誤差由此表明該方法相對于其他方法其擬合精度較差。