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

許弘喆, 程素珍

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

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

1 極值Ⅰ型分布理論

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

假設某觀測站歷年水位特征觀測值為h1，h2，…hn，觀測值呈極值Ⅰ型分布，則極大值h的分布函數為：

F(hp)=f(h

(1)

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

(2)

(3)

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

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

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

(4)

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

(5)

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

(6)

令：

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

(7)

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

2.1 MATLAB實現參數估計

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

矩法實現參數的統計語句是：

fori=1∶n

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

end

最大似然法通過mle語句實現參數的統計，即：

p=mle('ev',x);

σ=p(2)

2.2 MATLAB實現經驗頻率計算

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

(8)

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

MATLAB軟件的實現方式是：

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

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

2.3 MATLAB實現極值Ⅰ型分布理論計算

2.4 MATLAB能同步實現適線法擬合

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

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

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

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

2.5 MATLAB軟件實現黑森機率格紙生成[13]

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

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]

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

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]);

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

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

3 工程實例

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

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

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

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

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

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

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

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

4 結 語

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

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