基于Matlab的P-Ⅲ型頻率曲線優化適線程序的開發與應用

張曉艷 席秋義

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

1 引言

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

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

2 P-Ⅲ型曲線

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

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

當α、β、a0參數確定后，該P-Ⅲ型密度函數即可隨之確定。這3個參數與總體的統計參數、Cv、Cs具有如下關系：

3 xp推求

3.1 傳統推求方法

水文計算中，通常需要求出指定頻率相應的隨機變量取值xp，p與xp存在如下關系：

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

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

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

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

3.2 Matlab中xp的推導

Matlab本身有大量的預定義函數，因此，實現時可不采用傳統方法。下面進行公式推導：

Gamma概率密度函數如式（6）所示：

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

通過比較式（1）和式（8）可以看出，P-Ⅲ型概率密度函數經參數值轉換后，與Gamma函數完全一致。

由于Matlab中提供了已知概率計算相應分位值的Gamma反函數gaminv，具體形式為：

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

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

由上述過程可見，該方法可避免傳統復雜的查表插值計算過程，能有效提高xp計算速度和精度。

4 海森頻率格紙的繪制

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

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

Matlab語言正態分布的反函數為：

假設以作為海森格紙的繪制參考橫坐標0點，則各個頻率對應的海森格紙橫坐標由下式確定：

采用Matlab編制程序，海森格紙是通過控制背景網格的橫坐標來實現的，語句如下：

5 經驗點據與理論頻率曲線繪制

5.1 經驗頻率的計算

其中，P為≥xm的經驗頻率；m為xm的序號，即≥xm的項數；n為系列的總項數。

5.2 理論頻率曲線

理論頻率曲線統計參數的估算方法主要有兩種，一類是參數估計法，包括矩法、概率權重矩法等；另一類是適線法，包括經驗適線法、（殘差平方和最小）最小二乘法等。

實際水文頻率計算常通過目估適線，最終確定頻率曲線的統計參數，這在很大程度上依賴于計算者的實際經驗。為此，可以殘差平方和最小為適線準則，采用粒子群算法（PSO）對P-Ⅲ型頻率曲線參數進行優化，初始參數和搜索區間由矩法估算。PSO優化算法是一種基于群智能的優化算法，它對目標函數要求低，不需要進行求導計算，全局優化能力強，求解速度快。

所建模型目標函數如下：

其中，xi為實測數據；xpi為計算數據，由式（10）獲得。

5.3 經驗點據與理論頻率曲線繪制

海森格紙坐標中經驗點據的繪制命令如下：

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

5.4 程序流程

該P-III型頻率適線軟件實現流程如圖1所示：

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

圖2 PSO統計參數尋優迭代過程

圖3 漢江某水庫72h產流預報誤差特征點頻率分布圖

6 應用舉例

漢江某水庫其72小時產流預報誤差經檢驗，認為其服從P-Ⅲ型分布，利用矩法計算其統 計 參 數 分 別 為=1.0538、Cv=7.6412、Cs/Cv=0.04。

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

PSO粒子尋優迭代適應度變化過程如圖2所示，從圖中可見，采用PSO在30代內就計算得到了最優統計參數，速度很快。

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

7 結語

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

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

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

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

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

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

[6]華東水利學院.水文學的概率統計基礎[M],北京:水利出版社,1981.

[7]DL/T 5084-1998,電力工程水文技術規程[S].

[8]SL 44-93,水利水電工程設計洪水計算規范[S].