皮爾遜Ⅲ型曲線參數(shù)優(yōu)化估計方法對比研究

李慶昕

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

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

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

1 基本理論

1.1 理論頻率

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

(1)

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

(2)

1.2 經(jīng)驗頻率

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

(3)

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

對于特大洪水可利用下述公式進(jìn)行經(jīng)驗頻率的計算。

(4)

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

(5)

1.3 適線準(zhǔn)則

對水文序列值所對應(yīng)的頻率可利用上述伽馬函數(shù)的數(shù)值積分法進(jìn)行求解并以此可構(gòu)建目標(biāo)函數(shù)，表達(dá)式為：

(6)

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

2 數(shù)值積分法

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

2.1 梯形算法

(7)

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

2.2 辛普森算法

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

(8)

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

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

(9)

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

2.3 Romberg算法

(10)

其計算步驟和相應(yīng)公式分別如下所示：

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

2.4 自適應(yīng)誤差積分法

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

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

(11)

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

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

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

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

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

(12)

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

(13)

3 尋優(yōu)建模

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

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

(14)

式中，rand—0～1之間的隨機(jī)數(shù)字；floor—實數(shù)和整數(shù)之間的映射關(guān)系。

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

(15)

群體內(nèi)各蜘蛛的自身重量的計算如下：

(16)

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

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

(17)

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

4 數(shù)據(jù)檢驗

4.1 理想數(shù)據(jù)檢驗

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

表1 理想數(shù)據(jù)參數(shù)計算結(jié)果

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

4.2 討論分析

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

5 實測數(shù)據(jù)

利用相關(guān)文獻(xiàn)[12]和資料中對水文序列的計算方法和理論可對曲線參數(shù)進(jìn)行擬合估計，計算結(jié)果見表2，表中洪水序列1和2分別代表一般洪水和特大洪水序列。

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

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

6 結(jié)論

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

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

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