基于MATLAB的年最高潮水位極值Ⅰ型分布頻率分析計算

許弘喆, 程素珍

(1.蘭州大學, 甘肅 蘭州 730107; 2.山東省水利科學研究院， 山東 濟南 250013)

極值Ⅰ型、皮爾遜Ⅲ型分布曲線是進行最高潮水位頻率分析的主要理論頻率曲線[1]，目前利用MATLAB開發(fā)環(huán)境進行皮爾遜Ⅲ分布曲線計算較多[2]，如孫玉芳[3]的基于MATLAB開發(fā)環(huán)境的皮爾遜Ⅲ型頻率曲線分析軟件開發(fā)，李揚[4]的極值指數(shù)估計在太原站月降水頻率分析中的應用研究等，均利用MATLAB軟件進行皮爾遜Ⅲ型頻率曲線分析[5-11]。但有關(guān)利用極值Ⅰ型分布進行最高潮水位頻率分析的應用較少，且多側(cè)重于序列特征參數(shù)估計和計算方法對比研究，基本沒有涉獵基于極值Ⅰ型分布適線法在工程中的應用研究，現(xiàn)有相關(guān)規(guī)范均為參數(shù)查表、公式計算相應頻率下的水文特征設計值的計算方法。為推廣普及極值Ⅰ型分布在最高潮水位頻率分析中的應用，開發(fā)一種應用軟件是十分必要的。本文結(jié)合威海石島海陽觀測站年最高潮水位頻率分析，開發(fā)了一款基于MATLAB軟件應用環(huán)境下，極值Ⅰ型分布適線法在最高潮水位頻率分析計算中的應用，并成功在工程中應用，取得了理想效果。

1 極值Ⅰ型分布理論

極值分布是指在概率論中極大值(或者極小值)的概率分布，耿貝爾在1941年首次將極值Ⅰ型分布應用于洪水頻率分析工作，其密度函數(shù)符合降雨、徑流、潮水位等水文參數(shù)的分布規(guī)律。因此廣泛應用于水文特征參數(shù)的頻率分析計算。

假設某觀測站歷年水位特征觀測值為h1，h2，…h(huán)n，觀測值呈極值Ⅰ型分布，則極大值h的分布函數(shù)為：

F(hp)=f(h

(1)

水位序列均值和均方差采用矩法進行估算，見公式(2)、公式(3)。

(2)

(3)

式中：hi為序列第i年的年特征水位；n為年最高水位序列項數(shù)。

在水工設計中，通常用超越概率來表征設計頻率，即

p=p(h≥hp)=1-F(hp)=1-e-e-α(hp-β)

(4)

對式(4)兩邊取2次對數(shù)有：

(5)

把α、β化簡代入式(5)得：

(6)

令：

λpn=-(0.45+0.78×ln(-ln(1-p))) ，式(6)轉(zhuǎn)化為：

(7)

2 MATLAB實現(xiàn)最高潮水位頻率分析的原理和方法[12-16]

2.1 MATLAB實現(xiàn)參數(shù)估計

首先建立觀測數(shù)據(jù)樣本矩陣，x1、x2……，可以單獨對樣本進行分析，也可以通過cat(1,x1,x2,……)語句進行樣本矩陣合成分析。

矩法實現(xiàn)參數(shù)的統(tǒng)計語句是：

fori=1∶n

u=u+X(i)*X(i);

end

最大似然法通過mle語句實現(xiàn)參數(shù)的統(tǒng)計，即：

p=mle('ev',x);

σ=p(2)

2.2 MATLAB實現(xiàn)經(jīng)驗頻率計算

經(jīng)驗頻率是驗證各項計算成果的重要依據(jù)，根據(jù)經(jīng)驗頻率公式求出相應水位的經(jīng)驗頻率值，繪制經(jīng)驗頻率點，并與相應的計算成果相對比，進行相應的水位值和參數(shù)的合理性分析。經(jīng)驗頻率計算方法為：在n項連序水位系列中，按大小順序排位的第m項水位的經(jīng)驗頻率Pm，采用下列數(shù)學期望公式計算：

(8)

式中：m為水位連序系列中的序位；Pm為第m項水位的經(jīng)驗頻率。

MATLAB軟件的實現(xiàn)方式是：

n=length(X)；X=sort(X)；X=fliplr(X)；

Pm=[[1∶n]/(n+1)]

2.3 MATLAB實現(xiàn)極值Ⅰ型分布理論計算

2.4 MATLAB能同步實現(xiàn)適線法擬合

適線法是在一定的適線準則下，求解與經(jīng)驗點據(jù)擬合最優(yōu)的頻率曲線的統(tǒng)計方法，可避免因觀測序列短而產(chǎn)生計算結(jié)果可靠性差等問題。采用矩法或最大似然法對觀測數(shù)據(jù)序列進行統(tǒng)計，求出一組參數(shù)作為初值計算理論曲線，根據(jù)理論頻率曲線與經(jīng)驗頻率點據(jù)的配合情況，通過經(jīng)驗判斷調(diào)整適配參數(shù)，選定一條與經(jīng)驗點據(jù)擬合良好的頻率曲線。

適線法計算過程：(1) 統(tǒng)計水位連續(xù)觀測年最大值，建立水位觀測序列，并依數(shù)值大小排列成遞減序列，根據(jù)式(8)按序號確定逐項的經(jīng)驗頻率，將逐項數(shù)值和頻率點繪在黑森機率格紙上；(2) 統(tǒng)計計算的水位序列的均值、均方差一組統(tǒng)計參數(shù)，該組參數(shù)作為適線法初值；(3) 根據(jù)極值Ⅰ型分布理論，采用初值計算理論頻率曲線；根據(jù)理論頻率曲線和經(jīng)驗頻率點的擬合程度，保持均值不變，通過經(jīng)驗判斷調(diào)整均方差，一般采用0.1倍的數(shù)值增減，重新擬合頻率曲線，并把每次擬合的頻率曲線連同經(jīng)驗點據(jù)繪制在同一張機率格紙上，進行比較分析，選定一條與經(jīng)驗點據(jù)擬合良好的頻率曲線，這條線即為經(jīng)驗頻率曲線，通過該線查找相應設計頻率的潮水位。

通過初估、適線和綜合對比分析，可以得到比較合理的、能夠滿足工程設計要求的水位經(jīng)驗頻率曲線。通過選定的經(jīng)驗適線，根據(jù)極值Ⅰ型理論計算相應設計頻率的特征參數(shù)，并且可借助經(jīng)驗頻率曲線的延長，推求相應稀遇頻率的特征參數(shù)設計值。

適線法計算能夠通過MATLAB軟件中的Evinv、Evcdf函數(shù)實現(xiàn)，能同步完成數(shù)據(jù)的統(tǒng)計、分析、繪圖和曲線擬合，并采用Norminv函數(shù)生成黑森機率格紙，成果顯示在同一張機率格紙，便于評價和分析。

2.5 MATLAB軟件實現(xiàn)黑森機率格紙生成[13]

最高潮水位頻率計算中采用的黑森機率格紙是一種特殊的坐標系統(tǒng)，其縱坐標為均勻分格的常規(guī)數(shù)學坐標，橫坐標與頻率值的標準正態(tài)分布分位數(shù)有關(guān)。由于標準正態(tài)分布分位數(shù)在P=50%處為零，而黑森機率格紙在P=0.01%時的橫坐標值為零。它是實現(xiàn)最高潮水位頻率分析成果的重要工具，MATLAB軟件可通過調(diào)用Norminv函數(shù)實現(xiàn)，具體的程序語言為：

PL=[0.01 0.05 0.5 1 2 5 10 20 30 40 50 60 70 80 90 95 98 99 99.9 99.99];

xx=norminv(PL/100,0,1);

yy=[1∶0.25∶2.5];

axis([min(xx),max(xx),min(yy),max(yy)])

set(gca,‘xtick’,xx,‘ytick’,yy)

set(gca,‘xticklabel’,PL,‘yticklabel’,yy)

grid on

2.6 圖表成果生成[14-15]

圖、表是計算分析結(jié)果的主要表現(xiàn)形式，清晰、美觀、準確的圖表格式便于設計人員的利用，MATLAB軟件具有強大的圖表生成功能，可以根據(jù)客戶的要求進行定制，下面是一個簡單的表格生成程序語言：

mab1=cat(1,chengguo1,chengguo2,………);

data=mab1;

f=figure(1);

colnames={‘0.1%’,‘0.2%’,‘0.5%’,‘1%’,

‘2%’, ………};

rnames={‘First’,‘Second’,‘Third’,‘f’,‘t’};

t=uitable(f, ‘Data’, data, ‘ColumnName’,colnames, ...‘RowName’, rnames...‘Position’, [20 80 500 200]);

程序運行完成后，計算結(jié)果立即以圖、表的形成呈現(xiàn)，滿足設計人員的分析和評價，但是MATLAB生成的圖、表為圖片格式，不便編輯，為此MATLAB設計了與其他軟件的連接功能，便于格式的轉(zhuǎn)換。

MATLAB提供使其能與Excel互動操作的Excel-link宏。Excel-link 使得數(shù)據(jù)在MATLAB與Excel之間隨意交換，以及在Excel下調(diào)用MATLAB的函數(shù)，它將MATLAB強大的數(shù)值計算功能、數(shù)據(jù)可視化功能與Excel的數(shù)據(jù)Sheet功能結(jié)合在起，因此可以把MATLAB生成的計算數(shù)據(jù)導入Excel進行編輯和應用。

3 工程實例

(1) 運行MATLAB程序，調(diào)入編好的*.m程序文件。

(2) 把原始資料輸入*.m程序的原始參數(shù)矩陣，根據(jù)資料類型，修改預先定義的圖表和坐標名稱。

(3) 根據(jù)工程需要，選定序列特征參數(shù)的估算方法。

(4) 運行*.m程序，即生成**工程水位累計頻率曲線和和**工程成果匯總表。

(5) 根據(jù)經(jīng)驗頻率曲線選擇原則，對生成的各條頻率適線和經(jīng)驗頻率點據(jù)相比較，若成果符合要求，計算完畢；否則，調(diào)整適配參數(shù)，重新運行*.m程序，直到獲得滿意的成果，計算完畢。成果格式見圖1和表1。

圖1 年最高潮水位頻率曲線

(6) 把計算成果圖片導入Excel，進行編輯，表1為經(jīng)Excel編輯的表格格式。

表1 工程設計水位計算成果表

4 結(jié) 語

(1) 利用MATLAB軟件編程進行極值Ⅰ型頻率分析，能夠把資料匯總、經(jīng)驗頻率點據(jù)、適線參數(shù)估算、多條適線擬合計算和成果生成同時完成，具有程序簡單、高效以及計算精度高的特點，便于在實際工作中推廣應用。

(2) MATLAB簡單易學，編制的程序結(jié)構(gòu)簡單，功能齊全，修改方便，可以使工程設計人員從繁重的資料整理、復雜運算中解脫出來，大大提高水利工作者的工作效率。