基于極大似然法的負(fù)偏水文序列參數(shù)估計方法研究

基于極大似然法的負(fù)偏水文序列參數(shù)估計方法研究

胡詩松a，陳進(jìn)a，b，尹正杰a,c

(長江科學(xué)院 a.流域水資源與生態(tài)環(huán)境科學(xué)湖北省重點實驗室；b.院長辦公室；

c.水資源綜合利用研究所，武漢430010)

摘要：針對皮爾遜Ⅲ型曲線參數(shù)估計方法中尺度參數(shù)β>0，無法適用于潮位、水位等負(fù)偏水文序列(β<0)參數(shù)估計的問題，從正偏序列與負(fù)偏序列的關(guān)系出發(fā)，提出了一種極大似然估計的新方法。該方法是將負(fù)偏序列先后轉(zhuǎn)換為正偏序列、伽馬分布序列，并利用負(fù)偏序列確定ξ初值，通過伽馬分布序列迭代計算得到α,β的極大似然估計值。以某潮位站44 a的最高潮位數(shù)據(jù)為算例，比較了負(fù)偏極大似然法、負(fù)偏矩法、正偏查表法3種估計方法對負(fù)偏序列擬合的精度，結(jié)果發(fā)現(xiàn)負(fù)偏極大似然估計法的3項擬合指標(biāo)：PPCC檢驗值為0.993，OLS檢驗值為0.010，KS檢驗值為0.052，均優(yōu)于負(fù)偏矩法、正偏查表法。為此，負(fù)偏極大似然法對負(fù)偏序列擬合的精度為最優(yōu)。

關(guān)鍵詞：水文頻率；負(fù)偏序列；極大似然估計；皮爾遜Ⅲ型曲線；MATLAB軟件

中圖分類號：TV12

文獻(xiàn)標(biāo)志碼：A

文章編號：1001-5485(2015)06-0116-04

Abstract:The commonly used probability density function of Pearson-III distribution is invalid for negative-skewness water level and tide level series (β<0). To solve this problem, a means of maximum likelihood estimation (MLE) based on the relationship between negative-skewness series and positive-skewness series is proposed. In this method negative-skewnessPearson-Ⅲ distribution is translated into positive-skewnessPearson-Ⅲ and gamma

收稿日期：2015-04-02；修回日期：2015-04-13

作者簡介：李凌霄(1991-)，男，湖北武漢人，碩士研究生, 主要從事信息技術(shù)開發(fā)與應(yīng)用方面的研究，(電話)18627744329(電子信箱)lukelee0329@gmail.com。

DOI:10.3969/j.issn.1001-5485.2015.06.021

1研究背景

在水文計算中常需推求100 a一遇、甚至1 000 a一遇的設(shè)計值，水文頻率曲線是合理解決這一問題的有效“工具”。根據(jù)《水利水電工程設(shè)計洪水計算規(guī)范》(SL44—2006)，我國水文頻率計算線型采用皮爾遜Ⅲ型曲線。目前常用的皮爾遜Ⅲ型曲線尺度參數(shù)β>0，適用于Cs>0的正偏水文序列[1]，但分析水庫水位、潮位等水文序列時經(jīng)常出現(xiàn)負(fù)偏的情況，此時如仍沿用上述曲線配線，解得的理論頻率曲線和經(jīng)驗點偏離較大。為解決該問題，很多學(xué)者進(jìn)行了相關(guān)研究。如蔡體錄[2]研究了負(fù)偏和正偏頻率的關(guān)系，提出了一種基于查表的“三點法”；張濤等[3]提出用多項式代替皮爾遜Ⅲ型曲線配線以提高擬合精度；李興拼等[4]給出了基于矩法的正負(fù)偏轉(zhuǎn)換公式；張家鳴等[5]提出了一種直接求解負(fù)偏皮爾遜Ⅲ型曲線的方法。

國家標(biāo)準(zhǔn)規(guī)定使用其他線型配線需要單獨論證，由于多項式法通用性不強(qiáng)，而用矩法直接估計偏態(tài)系數(shù)Cs抽樣誤差又太大。因此，針對這2方面的問題，本文從正偏負(fù)偏頻率曲線的轉(zhuǎn)換關(guān)系出發(fā)，提出了比矩法精度更高的負(fù)偏水文序列皮爾遜Ⅲ型曲線參數(shù)估計的極大似然估計方法(后文統(tǒng)一簡稱為負(fù)偏極大似然法)，并以某潮位站44 a的最高潮位數(shù)據(jù)為算例，比較負(fù)偏極大似然法、負(fù)偏矩法、正偏查表3種方法對負(fù)偏序列擬合的精度，同時，對本文提出的方法加以了分析和驗證。

2計算原理

負(fù)偏水文序列的皮爾遜Ⅲ型曲線計算主要有2種思路：①尋找正偏負(fù)偏序列的關(guān)系，把負(fù)偏序列計算轉(zhuǎn)化為正偏，然后查表或者用矩法等方法估計參數(shù)；②直接根據(jù)負(fù)偏皮爾遜Ⅲ型曲線密度函數(shù)求解，需先計算該密度函數(shù)的積分然后用矩法或極大似然等方法估計參數(shù)。

皮爾遜Ⅲ型曲線密度函數(shù)為

(1)

式中：α為形狀參數(shù)(α>0 )；β為尺度參數(shù)；ξ為位置參數(shù)；x是樣本值；Γ(α)為α的伽馬函數(shù)[6]。

當(dāng)β<0時式(1)為負(fù)偏序列對應(yīng)的皮爾遜Ⅲ型曲線密度函數(shù)。采用第1種思路計算時，如果用矩法估計參數(shù)，取100個樣本點時α的抽樣誤差為40%～126%之間[1]，序列樣本數(shù)一般不會超過100，所以實際抽樣誤差會更大；采用第2種思路計算時，如果用極大似然法估計參數(shù)，由于式(1)右邊有負(fù)數(shù)，不能同時對兩端取對數(shù)，直接進(jìn)行極大似然估計將十分困難。為了兼顧計算精度和簡易性，本文結(jié)合了2種思路，采用先把負(fù)偏序列轉(zhuǎn)換為正偏序列，再對正偏序列采用極大似然法對P-Ⅲ型曲線密度函數(shù)進(jìn)行估計參數(shù)的計算。

2.1負(fù)偏與正偏序列頻率曲線關(guān)系

(2)

(3)

《工程水文學(xué)》[1]中給出的正偏皮爾遜Ⅲ型曲線密度函數(shù)為

(4)

(5)

(6)

2.2參數(shù)α,β,ξ的極大似然估計

與矩法相比，極大似然估計不僅利用了樣本點的信息，還利用了樣本分布函數(shù)的信息，因此，相較矩法有一定的優(yōu)勢[7]。極大似然估計是用樣本在總體中出現(xiàn)幾率最大時的參數(shù)模擬樣本總體參數(shù)，對于樣本x1,x2,…,xn，設(shè)他們的總體密度函數(shù)為f(x,α)，則有

(7)

式(7)稱為似然函數(shù)，使似然函數(shù)最大的參數(shù)α為樣本的極大似然參數(shù)估計值。把式(4)代入式(7)就得到了正偏皮爾遜Ⅲ型曲線的似然函數(shù)，即

(8)

(9)

把估計出的α，1/β初值代入式(9)就可以解出ξ。具體的轉(zhuǎn)換及迭代求解步驟如下：

(1) 對原始負(fù)偏序列x1,x2,…,xn用矩法計算一個ξ的初值ξ1；

(3) 把α1,β1代入式(9)解出ξ2，注意β1與式(9)中的β是互為倒數(shù)的關(guān)系；

(4) 利用ξ2重復(fù)步驟②至③直到相鄰2次ξi-ξi-1足夠小時停止迭代，最后一次得到的ξi,αi-1,1/βi-1即為所需的極大似然估計值。

3實例應(yīng)用

圖1　最高潮位頻率曲線擬合結(jié)果對比 Fig.1　Fitting results of maximum tide level frequency curves

估計方法分布參數(shù)αβξCs正偏查表38.10451.6310.3690.324負(fù)偏矩法7.384-13.6370.782-0.736負(fù)偏極大似然法4.923-17.9400.835-0.901

從圖1中可以直觀地發(fā)現(xiàn)負(fù)偏極大似然法和正偏查表2種方法擬合效果較好，正偏查表(取Cs=3Cv)在頻率區(qū)間兩端處有所偏離；負(fù)偏矩法擬合效果最差，整個概率區(qū)間上偏離經(jīng)驗點距較大。

結(jié)合表1分析上述3種估計方法產(chǎn)生不同擬合效果的原因：對負(fù)偏水文序列采用正偏查表擬合在頻率區(qū)間兩端處偏離經(jīng)驗點距，是因為我們忽略了樣本序列的負(fù)偏特性，強(qiáng)制令Cs=3Cv>0，而Cs越大皮爾遜Ⅲ型曲線的上段越陡，下段越平緩[9]；正偏查表和負(fù)偏極大似然擬合效果差別不顯著是因為本例中Cs=-0.736較小，加上潮位序列本身變化范圍很小(0.75~1.34 m)；負(fù)偏矩法α=7.384，與負(fù)偏極大似然法α=4.923誤差達(dá)50%，而α是影響皮爾遜Ⅲ型曲線形狀的敏感因子，因此， 負(fù)偏矩法考慮負(fù)偏序列特性擬合效果依然不理想。

為了更加客觀地比較3種估計方法的擬合效果，下面分別對3種估計方法進(jìn)行PPCC(概率點據(jù)相關(guān)系數(shù))檢驗[10]、OLS(離差平方和)檢驗[11]、KS(柯爾莫哥洛夫-斯摩洛夫)檢驗[12]。這3種檢驗方法從不同角度描述擬合效果，PPCC檢驗反映了模擬值和樣本值的相關(guān)性，數(shù)值越大擬合效果越好；OLS檢驗反映了離差的累積情況，數(shù)值越小擬合效果越好；KS檢驗反映了擬合最差點的離差，數(shù)值越小擬合效果越好。3種檢驗方法的公式及3種估計方法的檢驗值見表2。

表2擬合優(yōu)度檢驗值

Table 2Test values of fitting efficiency

估計方法PPCC檢驗OLS檢驗KS檢驗Σni=1[(xi-xm)(yi-ym)]Σni=1[(xi-xm)2(yi-ym)2]Σni=1[(xi-yi)]2 max(xi-yi)正偏查表0.9630.0310.147負(fù)偏矩法0.9922.1470.355負(fù)偏極大似然0.9930.0100.052

注:xi,yi分別為同頻率對應(yīng)經(jīng)驗曲線值和擬合曲線值；xm,ym分別為xi,yi均值。

表2結(jié)果表明:負(fù)偏矩法除了在PPCC檢驗所得檢驗值優(yōu)于正偏查表法外，其他指標(biāo)均劣于正偏查表法；負(fù)偏極大似然法3種檢驗方法所得檢驗值都是最優(yōu)的，擬合度最優(yōu)；3種估計方法的PPCC檢驗所得指標(biāo)都比較好，說明最高潮位曲線采用皮爾遜Ⅲ型曲線作為分布密度函數(shù)是符合實際情況的；OLS和KS檢驗法所得檢驗指標(biāo)正偏查表分別是負(fù)偏極大似然的3倍和2.7倍，說明本文提出的負(fù)偏極大似然法在累積離差和擬合最差點離差方面擬合效果比正偏查表有較大提升。負(fù)偏矩法KS檢測指標(biāo)為0.355已經(jīng)超過了《水文情報預(yù)報規(guī)范》(GB/T22482—2008)潮位的誤差限，因此不能直接使用該方法的計算結(jié)果。

4結(jié)論

正偏查表法估計負(fù)偏序列在頻率區(qū)間兩端處擬合效果較差，無論正偏負(fù)偏序列用矩法估計參數(shù)Cs都是無效的。本文基于先把負(fù)偏序列轉(zhuǎn)換為正偏序列，再用極大似然法對轉(zhuǎn)化后的正偏序列進(jìn)行P-Ⅲ型曲線密度函數(shù)參數(shù)估計的計算思路，提出了一種新的極大似然迭代求解方法，并經(jīng)過了某潮位站44 a的最高潮位數(shù)據(jù)驗證。所提出的負(fù)偏極大似然估計法3項擬合指標(biāo)：PPCC檢驗值為0.993，OLS檢驗值為0.010，KS檢驗值為0.052，均優(yōu)于負(fù)偏矩法、正偏查表法。為此，負(fù)偏極大似然法對負(fù)偏序列擬合的精度為最優(yōu)。

建議采用皮爾遜Ⅲ型曲線估計水文系列時應(yīng)先判斷該水文系列是正偏還是負(fù)偏，如果是負(fù)偏系列，應(yīng)轉(zhuǎn)化為對應(yīng)的正偏系列后采用極大似然估計密度函數(shù)參數(shù)，這既可以提高頻率分布曲線估計的精度也便于使用計算機(jī)自動化計算。

參考文獻(xiàn)：

[1]詹道江，徐向陽，陳元芳.工程水文學(xué)[M]. 北京：中國水利水電出版社，2010. (ZHAN Dao-jiang, XU Xiang-yang, CHEN Yuan-fang. Engineering Hydrology[M]. Beijing: China Water Power Press, 2010. (in Chinese))

[2]蔡體錄. 負(fù)偏頻率曲線的計算[J]. 東海海洋,1983,1(4):1-7.(CAI Ti-lu. Calculation of Negative Deviation Distribution Frequency[J]. Donghai Marine Science, 1983, 1(4): 1-7. (in Chinese))

[3]張濤,王世勛,王祥三，等. 水位負(fù)偏分布頻率計算方法分析與研究[J].水文,2008,28(2):5-9. (ZHANG Tao, WANG Shi-xun, WANG Xiang-san,etal. Analysis of Negative Deviation Distribution Frequency for Water Stage Calculation[J]. Journal of China Hydrology, 2008,28(2):5-9.(in Chinese))

[4]李興拼,鄭江麗. 淺論P-Ⅲ型負(fù)偏頻率曲線計算方法[J]. 廣東水利水電, 2010, 9(9): 17-18. (LI Xing-pin, ZHENG Jiang-li. Calculation Method of Negative-skewness Series with Pearson-Ⅲ Frequency Curve[J]. Guangdong Water Resources and Hydropower, 2010, 9(9): 17-18. (in Chinese))

[5]張家鳴,陳曉宏,葉長青. Pearson-Ⅲ型頻率曲線對負(fù)偏水文序列的計算[J]. 水利學(xué)報,2012, 43(11):1296- 1301. (ZHANG Jia-ming, CHEN Xiao-hong, YE Chang-qing. Calculation of Negative-skewness Hydrological Series with Pearson-Ⅲ Frequency Curve[J]. Journal of Hydraulic Engineering, 2012, 43(11): 1296-1301. (in Chinese))

[6]STEDINGER J R, VOGEL R M, FOUFOULA-GEORGIOU E. Handbook of Hydrology. New York: McGraw-Hill, 1993.

[7]余泱悅,賀信. P-Ⅲ曲線的極大似然估計及應(yīng)用[J]. 人民長江, 2012, 43(21): 21-23. (YU Yang-yue, HE Xin. Maximum Likelihood Estimation of Pearson-Ⅲ Curve and Its Application[J]. Yangtze River, 2012, 43(21): 21-23. (in Chinese))

[8]PONCE V M. Engineering Hydrology[M]. New Jersey: Englewood Cliffs, 1989: 214-216.

[9]RAO A R, HAMED K H. Flood Frequency Analysis[M]. Boca Raton: CRC Press, 2000.

[10]涂新軍,陳曉宏. 基于PPCC檢驗法的枯水徑流概率分布選型研究[J]. 水電能源科學(xué),2006,24(1):76-79. (TU Xin-jun, CHEN Xiao-hong. Study on Probability Distribution of Low Flow Based on PPCC Method[J].Water Resources and Power, 2006,24(1):76-79. (in Chinese))

[11]成靜清. 非一致性年徑流序列頻率分析計算[D].陜西：西北農(nóng)林科技大學(xué),2010. (CHENG Jing-qing. Hydrological Frequency Calculation Principle of Inconsistent Annual Runoff Series[D]. Shaanxi: Northwest A& F University, 2010. (in Chinese))

[12]郭生練. 設(shè)計洪水研究進(jìn)展與評價[M]. 北京：中國水利水電出版社，2005. (GUO Sheng-lian. Advance and Assessment of Design Flood Hydrograph Methods[M]. Beijing: China Water Power Press, 2005. (in Chinese))

(編輯：陳紹選)

Maximum Likelihood Estimation of Negative-skewnessHydrological Series

HU Shi-song1,CHEN Jin1,2,YIN Zheng-jie1,3

(1.Key Lab of Basin Water Resource and Eco-environmental Science in Hubei Province,

Yangtze River Scientific Research Institute, Wuhan 430010, China; 2.Administration Office,

Yangtze River Scientific Research Institute, Wuhan430010, China; 3.Water Resources Department,

Yangtze River Scientific Research Institute, Wuhan430010, China)

distribution in sequence. The initial value ofξis decided by the negative-skewness distribution and the MLE ofαandβare calculated through the gamma distribution. In the iterative process the sampling error with evaluating coefficient of skewness will be avoided. Moreover, the 44-year record of maximum tide level at a tide station is used as calculation example. The fitting accuracy of negative-skewness MLE, negative-skewness moment method and positive-skewness table look-up are compared. The fitting results of negative-skewness MLE (test value of PPCC is 0.993, OLS 0.010, and KS 0.052) are superior to those of the other two methods. In conclusion, negative-skewness MLE has the optimum fitting accuracy for negative-skewness hydrological series.

Key words: hydrological frequency; negative-skewness series; MLE; Pearson-Ⅲ frequency curve; MATLAB software

歡 迎 訂 閱歡 迎 投 稿

2015,32(06):120-126