張曉艷 席秋義

(1.山西黃河河務(wù)局臨猗黃河河務(wù)局 山西 臨猗 044100；2.國網(wǎng)陜西省電力公司電力科學(xué)研究院 陜西 西安 710054)

1 引言

P-Ⅲ型頻率曲線廣泛應(yīng)用于水文、氣象頻率分析中。針對這種頻率曲線，國內(nèi)已開發(fā)了一些專用的分析計(jì)算軟件，如基于EXCEL和Visual Basic平臺開發(fā)的該類軟件[1～3]。這些軟件在計(jì)算某一頻率對應(yīng)的隨機(jī)變量值時，大多采用離均系數(shù)φ，在讀取離均系數(shù)表的基礎(chǔ)上，再通過插值計(jì)算得出，這樣做不僅編程比較繁瑣而且計(jì)算精度較低。此外，這些軟件對源代碼進(jìn)行了封裝，用戶僅能實(shí)現(xiàn)其規(guī)定的頻率計(jì)算基本功能，不能很好地適應(yīng)頻率分析的多種拓展需求。

Matlab是一種專業(yè)的計(jì)算機(jī)程序，用于工程科學(xué)的矩陣數(shù)學(xué)運(yùn)算。由于它提供了一個極其廣泛的預(yù)定義函數(shù)庫，眾多應(yīng)用顯示，在解決工程技術(shù)問題方面，Matlab比其它任何計(jì)算機(jī)語言（包括Fortran和C）都簡單高效[4]。為解決P-III型頻率曲線參數(shù)傳統(tǒng)計(jì)算方法工作量大，精度低的問題，本文探討了基于Matlab的P-III型頻率曲線適線軟件的實(shí)現(xiàn)方法，對某一頻率對應(yīng)的隨機(jī)變量值計(jì)算方法根據(jù)Matlab語言特點(diǎn)進(jìn)行改進(jìn)。

2 P-Ⅲ型曲線

皮爾遜Ⅲ型（簡稱P-Ⅲ型）分布P-Ⅲ型曲線是一條一端有限一端無限的不對稱單峰，正偏曲線，是規(guī)范中規(guī)定采用的一種線型，其概率密度函數(shù)為[5]：

其中，Γ（α）為伽瑪函數(shù)；α、β、a0分別為 P-Ⅲ型分布的形狀、尺度和位置參數(shù)，其中，α＞0，β＞0；a0稱為位置參數(shù)，是 x總體的最小值，即概率密度曲線起點(diǎn)的橫坐標(biāo)。

當(dāng)α、β、a0參數(shù)確定后，該P(yáng)-Ⅲ型密度函數(shù)即可隨之確定。這3個參數(shù)與總體的統(tǒng)計(jì)參數(shù)、Cv、Cs具有如下關(guān)系：

3 xp推求

3.1 傳統(tǒng)推求方法

水文計(jì)算中，通常需要求出指定頻率相應(yīng)的隨機(jī)變量取值xp，p與xp存在如下關(guān)系：

直接由上積分形式計(jì)算是非常繁雜的，實(shí)際做法是通過變量轉(zhuǎn)換，根據(jù)擬定的Cs值進(jìn)行積分，并將成果制成專用表格供查用，變換成下面的積分形式：

其中，被積函數(shù)只含有一個特定參數(shù)Cs，其它兩個參數(shù)、Cv都包含在Φ中，因而只要假定一個Cs值，便可從式（4）通過積分求出P與Φp之間的關(guān)系[6]。

在頻率計(jì)算時，由已知的Cs值，查Φ值表得出不同P的Φp值，然后利用已知的、Cv值，通過下式即可求出與各種p相應(yīng)的xp值。

其中，Φp為離均系數(shù)，其均值為零，均方差為1。

3.2 Matlab中xp的推導(dǎo)

Matlab本身有大量的預(yù)定義函數(shù)，因此，實(shí)現(xiàn)時可不采用傳統(tǒng)方法。下面進(jìn)行公式推導(dǎo)：

Gamma概率密度函數(shù)如式（6）所示：

代入式（1），則得到下式：

通過比較式（1）和式（8）可以看出，P-Ⅲ型概率密度函數(shù)經(jīng)參數(shù)值轉(zhuǎn)換后，與Gamma函數(shù)完全一致。

由于Matlab中提供了已知概率計(jì)算相應(yīng)分位值的Gamma反函數(shù)gaminv，具體形式為：

先說說我對“夕陽”的定位。50歲左右，竊以為不能算作“夕陽”，應(yīng)是“下午四五點(diǎn)鐘的太陽”。但是，以我長期所處的縣、鄉(xiāng)級為例，很多這個年齡段的人，都以不同形式“下崗”了，他們往往也自稱“夕陽”，但非常勉強(qiáng)。65歲之前，只是剛踏上“夕陽”的邊，雖然精力、體力比不上以往，但仍有自己的優(yōu)勢，還不失為人生的一個“黃金期”，仍可“大有作為”。

因此，可根據(jù)P-Ⅲ型分布函數(shù)的參數(shù)α和β，代入式（7）得到α'和β'，通過函數(shù)式（9）即可得到不同頻率p對應(yīng)的Gamma函數(shù)值x'p，P-Ⅲ型頻率p對應(yīng)的xp可通過式（10）得到：

由上述過程可見，該方法可避免傳統(tǒng)復(fù)雜的查表插值計(jì)算過程，能有效提高xp計(jì)算速度和精度。

4 海森頻率格紙的繪制

在普通方格紙上繪制P-Ⅲ頻率曲線，因頻率曲線的兩端特別陡峭，且受圖幅的限制，對于特小頻率或特大頻率，尤其是特大頻率的點(diǎn)子很難點(diǎn)在圖上。為此，水文計(jì)算中采用海森頻率格紙繪制P-Ⅲ頻率曲線。

正態(tài)頻率曲線在普通格紙上是一條規(guī)則的S形曲線，它在50%前后的曲線方向雖然相反，但形狀完全一樣。海森格紙其橫坐標(biāo)的分劃是按反標(biāo)準(zhǔn)正態(tài)頻率曲線拉成一條直線的原理計(jì)算出來的。這種頻率格紙的縱坐標(biāo)仍是普通分格，但橫坐標(biāo)的分格是不相等的，中間分隔較密，越往兩端分格越稀，其間距關(guān)于P=50%是對稱的。

Matlab語言正態(tài)分布的反函數(shù)為：

假設(shè)以作為海森格紙的繪制參考橫坐標(biāo)0點(diǎn)，則各個頻率對應(yīng)的海森格紙橫坐標(biāo)由下式確定：

采用Matlab編制程序，海森格紙是通過控制背景網(wǎng)格的橫坐標(biāo)來實(shí)現(xiàn)的，語句如下：

5 經(jīng)驗(yàn)點(diǎn)據(jù)與理論頻率曲線繪制

5.1 經(jīng)驗(yàn)頻率的計(jì)算

其中，P為≥xm的經(jīng)驗(yàn)頻率；m為xm的序號，即≥xm的項(xiàng)數(shù)；n為系列的總項(xiàng)數(shù)。

5.2 理論頻率曲線

理論頻率曲線統(tǒng)計(jì)參數(shù)的估算方法主要有兩種，一類是參數(shù)估計(jì)法，包括矩法、概率權(quán)重矩法等；另一類是適線法，包括經(jīng)驗(yàn)適線法、（殘差平方和最?。┳钚《朔ǖ取?/p>

實(shí)際水文頻率計(jì)算常通過目估適線，最終確定頻率曲線的統(tǒng)計(jì)參數(shù)，這在很大程度上依賴于計(jì)算者的實(shí)際經(jīng)驗(yàn)。為此，可以殘差平方和最小為適線準(zhǔn)則，采用粒子群算法（PSO）對P-Ⅲ型頻率曲線參數(shù)進(jìn)行優(yōu)化，初始參數(shù)和搜索區(qū)間由矩法估算。PSO優(yōu)化算法是一種基于群智能的優(yōu)化算法，它對目標(biāo)函數(shù)要求低，不需要進(jìn)行求導(dǎo)計(jì)算，全局優(yōu)化能力強(qiáng)，求解速度快。

所建模型目標(biāo)函數(shù)如下：

其中，xi為實(shí)測數(shù)據(jù)；xpi為計(jì)算數(shù)據(jù)，由式（10）獲得。

5.3 經(jīng)驗(yàn)點(diǎn)據(jù)與理論頻率曲線繪制

海森格紙坐標(biāo)中經(jīng)驗(yàn)點(diǎn)據(jù)的繪制命令如下：

理論頻率曲線的繪制命令如下：

5.4 程序流程

該P(yáng)-III型頻率適線軟件實(shí)現(xiàn)流程如圖1所示：

圖1 P-III型頻率曲線適線繪制程序流程

圖2 PSO統(tǒng)計(jì)參數(shù)尋優(yōu)迭代過程

圖3 漢江某水庫72h產(chǎn)流預(yù)報(bào)誤差特征點(diǎn)頻率分布圖

6 應(yīng)用舉例

漢江某水庫其72小時產(chǎn)流預(yù)報(bào)誤差經(jīng)檢驗(yàn)，認(rèn)為其服從P-Ⅲ型分布，利用矩法計(jì)算其統(tǒng) 計(jì) 參 數(shù) 分 別 為=1.0538、Cv=7.6412、Cs/Cv=0.04。

粒 子 搜 索 區(qū) 間 設(shè) 為 ：Cv：[Cv-0.2Cv，Cv+0.2Cv]；Cs/Cv：[Cs/Cv-0.2Cs/Cv，Cs/Cv+0.2Cs/Cv]；均值固定不變。粒子群體大小Swarmsize=20，iter=30。

PSO粒子尋優(yōu)迭代適應(yīng)度變化過程如圖2所示，從圖中可見，采用PSO在30代內(nèi)就計(jì)算得到了最優(yōu)統(tǒng)計(jì)參數(shù)，速度很快。

利用該程序繪制的該水庫72小時產(chǎn)流預(yù)報(bào)誤差分布曲線圖，如圖3所示，其顯示的海森機(jī)率格紙和頻率曲線形象逼真直觀，能滿足頻率分析計(jì)算的需要。

7 結(jié)語

本文利用Matlat強(qiáng)大的預(yù)定義函數(shù)庫功能、圖形顯示功能和運(yùn)算速度快、易于編程的特點(diǎn)，開發(fā)了P-Ⅲ型頻率曲線適線及繪制程序，該程序克服了傳統(tǒng)頻率計(jì)算中查表插值計(jì)算工作繁瑣、精度低的缺點(diǎn)，結(jié)合PSO群智能算法，很方便地實(shí)現(xiàn)對其統(tǒng)計(jì)參數(shù)的優(yōu)化。通過matlab本身預(yù)定義的函數(shù)直接計(jì)算頻率相應(yīng)的隨機(jī)變量值，編程量相對較小，可準(zhǔn)確體現(xiàn)水文、電力工程有關(guān)規(guī)范對水文頻率計(jì)算和繪制的標(biāo)準(zhǔn)化、規(guī)范化需求。本程序開發(fā)時，適線準(zhǔn)則僅采用了殘差平方和最小準(zhǔn)則，殘差絕對值和最小和相對殘差平方和最小準(zhǔn)則并沒有考慮，用戶可根據(jù)需要將目標(biāo)函數(shù)轉(zhuǎn)換成適合的類型。陜西水利

[1]耿鴻江.Excel在水文計(jì)算中的應(yīng)用[J],水資源研究.2004,25(3):48～49.

[2]劉煥彬.氣象參數(shù)極值理論頻率曲線的Excel實(shí)現(xiàn)[J],河南氣象.2006(2):70～72.

[3]趙培穎,金冶,張忠孝.Visual Basic在繪制PIII型頻率曲線中的應(yīng)用[J],水利規(guī)劃與設(shè)計(jì).2008(2):55～57.

[4]Stephen J.Chapman.MATLAB編程（第二版）[M],北京:科學(xué)出版社,2003.

[5]葉守澤.水文水利計(jì)算[M],北京:水利電力出版社,1992.

[6]華東水利學(xué)院.水文學(xué)的概率統(tǒng)計(jì)基礎(chǔ)[M],北京:水利出版社,1981.

[7]DL/T 5084-1998,電力工程水文技術(shù)規(guī)程[S].

[8]SL 44-93,水利水電工程設(shè)計(jì)洪水計(jì)算規(guī)范[S].