基于交互熵原理的水文分布參數(shù)估計研究

牛林森，宋松柏

(西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院，陜西 楊凌 712100)

基于交互熵原理的水文分布參數(shù)估計研究

牛林森，宋松柏

(西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院，陜西 楊凌 712100)

[摘要]【目的】 研究交互熵法進行水文分布參數(shù)估計的普適性。【方法】 應(yīng)用最小交互熵原理研究Gumbel分布參數(shù)估計，在此基礎(chǔ)上應(yīng)用蒙特卡洛試驗檢驗交互熵法統(tǒng)計性能，然后結(jié)合矩法和線性矩法等傳統(tǒng)參數(shù)估計方法，以陜西省關(guān)中地區(qū)周至、武功、蒲城、禮泉、白水、潼關(guān)6個水文站年降水序列為例，計算年降水量設(shè)計值并擬合實測值序列，利用累積相對偏差平方和評價理論年降水量頻率曲線對實測值序列的擬合效果。【結(jié)果】 蒙特卡洛試驗檢驗表明，交互熵法所求設(shè)計值的有效性指標估計量標準偏差(SE)和均方根誤差(RMSE)小于矩法和線性矩法,偏差指標控制在7%以內(nèi)；交互熵法估計周至、武功、蒲城、禮泉、白水、潼關(guān)6個水文站的累積相對偏差平方和分別為0.000 017 68，0.000 065 62，0.000 019 66，0.000 063 00，0.000 014 12和0.000 016 61，線性矩法估計上述6站的累積相對偏差平方和分別為0.000 087 62，0.000 093 55，0.000 086 52，0.000 101 39，0.000 065 15和0.000 069 05，矩法估計上述6站的累積相對偏差平方和分別為0.000 108 74，0.000 125 40，0.000 092 41，0.000 127 65，0.000 085 49和0.000 093 57。由此可知，交互熵法不僅具有較好的有效性與合理的不偏性，而且與實測序列的擬合效果也明顯優(yōu)于傳統(tǒng)方法?！窘Y(jié)論】 交互熵法是一種可行的水文分布參數(shù)估計方法,能有效提高Gumbel分布參數(shù)的估計精度。

[關(guān)鍵詞]水文分布；參數(shù)估計；交互熵法；Gumbel分布；分位數(shù)對約束

水文頻率分析是水文學(xué)研究的重要內(nèi)容，其通過研究和分析水文隨機現(xiàn)象，揭示其中蘊含的統(tǒng)計規(guī)律，并對未來可能的情勢做出預(yù)估，以滿足水利工程規(guī)劃、設(shè)計、管理以及水資源利用等工作的需要[1]。水文頻率分析包括分布線型選擇和參數(shù)估計[2]，其中水文分布參數(shù)估計方法研究是水文頻率分析的重要內(nèi)容之一?，F(xiàn)有的水文分布參數(shù)估計方法有矩法、極大似然法、概率權(quán)重矩法、線性矩法等[3]。除矩法之外，其他方法在求解參數(shù)時都會受到線型種類的影響，即對于不同分布線型需要分別推導(dǎo)參數(shù)求解表達式和(或)方程(組)，計算較為復(fù)雜。矩法雖然不受分布線型影響，但是精度較差(尤其是Cs值),一般用于初估參數(shù)值[4]。尋求參數(shù)求解簡便、精確度高的參數(shù)估計方法是水文頻率分析領(lǐng)域的重要課題。交互熵(Cross entropy)又稱相對熵(Relative entropy)、差異信息(Discrimination information)等，由Good[5]于1950年提出，表示兩個分布之間的概率距離,其是Shannon信息熵的一般化,能夠有效避免Shannon信息熵的缺點[6]。Kullback[7]于1959年提出了最小交互熵原理(Minimum cross entropy principle)，該原理借助了Jaynes的最大熵原理(Maximum entropy principle)，即所求分布即是不確定度最大的分布[8-9]，并采用類似貝葉斯法中先驗分布和后驗分布的概念，若先驗分布(亦稱參考分布)已知，且約束條件一定，則候選分布中能使交互熵函數(shù)值最小的分布即為所求分布,也稱后驗分布[10-14]。交互熵法廣泛應(yīng)用于結(jié)構(gòu)分析和光譜分析等領(lǐng)域[15-16]。1988年，Lind等[17]首次將交互熵法應(yīng)用于水文領(lǐng)域，其以分位數(shù)對為約束條件，依據(jù)最小交互熵原理，采用Gamma分布和Gumbel分布，以加拿大2條河流的年最大洪峰流量數(shù)據(jù)為例，運用交互熵法估計分布參數(shù)并繪制出擬合效果圖，然后將交互熵法所求熵值與傳統(tǒng)的矩法和極大似然法相對比，得出了交互熵法較優(yōu)的結(jié)論。美國學(xué)者Hosking于1990年提出的線性矩法是概率權(quán)重矩的線性組合，與概率權(quán)重矩法結(jié)果接近，并且具有良好的統(tǒng)計性能[18]。目前，國內(nèi)尚缺少將交互熵法應(yīng)用于水文領(lǐng)域的研究報道[19]。為此，本研究采用Gumbel分布作為參考分布，以矩法和線性矩法等傳統(tǒng)方法為對比，應(yīng)用蒙特卡洛試驗分析交互熵法的統(tǒng)計性能，并以陜西省關(guān)中地區(qū)6個水文站年降水序列為例，對交互熵法應(yīng)用于水文領(lǐng)域的普適性進行了分析，以期為陜西關(guān)中地區(qū)水利工程設(shè)計中設(shè)計洪水的準確計算提供參考。

1交互熵法原理

1.1交互熵函數(shù)

設(shè)隨機變量X的參考分布函數(shù)為P(x),密度函數(shù)為p(x),同時假設(shè)后驗分布為Q(x)，其對應(yīng)密度函數(shù)為q(x)，則交互熵函數(shù)D的形式為：

(1)

式中:D為交互熵函數(shù)，x為隨機變量X的實際觀測值。

1.2約束條件

選取分位數(shù)對為約束條件。將隨機變量X的觀測序列按升序排列，得到新序列S={xi},i=1，2，… ，r，其中xi∈R，r為觀測序列中觀測值的個數(shù)。假設(shè)x是一個可能發(fā)生的未知量，那么其落在X由xi分成的r+1個子區(qū)間[x0,x1)，[x1,x2),…,[xr,xr+1)內(nèi)的機率相等。根據(jù)樣本規(guī)則，對應(yīng)的分位數(shù)概率為i/(r+1)，假設(shè)Q(x|x1，x2，…，xr)是由序列S推導(dǎo)出的X的分布函數(shù),q(x|x1，x2，…，xr)為相應(yīng)密度函數(shù)。則有：

(2)

式中：r為X觀測序列中觀測值的個數(shù)。

1.3Kullback最小交互熵原理

給定一個參考分布函數(shù)P(x)和密度函數(shù)p(x)，在分位數(shù)x=xi上，P(x)值記為Pi。通過尋求后驗分布函數(shù)Q(x|x1，x2，…，xr)，使交互熵函數(shù)為：

(3)

令交互熵函數(shù)在滿足分位數(shù)對約束的條件下最小，則式(2)可寫為期望值的形式：

gi∶∫Ifi(x)q(x|x1,x2,…,xr)dx-

(4)

[logq(x|x1,x2,…,xr)-logp(x)]dx+

(5)

為使交互熵最小，根據(jù)變分法中的Euler-Lagrange方程，將式(5)僅對q(x)求偏導(dǎo)，并令其為0，有：

(6)

整理，有：

q(x|x1,x2,…,xr)=

(7)

其中：

(8)

當(dāng)x∈Ii時，有：

μ(x)=μi=exp(-1-λi)。

(9)

將式(7)代入式(4)的左邊，考慮x∈Ii時fi(x)=1，故有：

(10)

根據(jù)式(4)可知，式(10)等于0，故有：

(11)

式中：Pi和Pi+1分別為參考分布在x=xi和x=xi+1時的概率分布函數(shù)。因此，當(dāng)x∈Ii，有：

q(x|x1,x2,…,xr)=μip(x)。

(12)

將式(9)和式(10)代入式(2)中，有：

(13)

化簡整理，得：

D(q,p)=-log(r+1)-

(14)

令c=log(r+1)，c為非負數(shù)。再引入中間變量S(P)：

(15)

則有：

Dmin(q,p)=-c+S(P)/(r+1)。

(16)

式中：Dmin(q,p)可以稱為最小交互熵函數(shù)。

從式(16)可以看出，最小交互熵函數(shù)可以寫成參考分布的函數(shù)，并且S(P)值越小，最小交互熵函數(shù)值越小。在參考分布曲線類型選定的條件下，最小交互熵原理可解釋為最優(yōu)分布的最小交互熵函數(shù)值最小。

1.4參考分布

選擇Gumbel分布作為參考分布，其密度函數(shù)與分布函數(shù)分別為：

p(x)=αexp(-α(x-u)-e-α(x-u)),

-∞

(17)

P(x)=exp(-e-α(x-u))。

(18)

式中：α為尺度參數(shù)，μ為位置參數(shù)。

2蒙特卡洛試驗

本研究利用蒙特卡洛試驗進行交互熵法統(tǒng)計性能研究，并與傳統(tǒng)的參數(shù)估計方法，即矩法、線性矩法的估計結(jié)果進行比較。

2.1評價標準

若估計量的偏差越接近0，表明估計量的不偏性越好，在M次統(tǒng)計試驗中，其計算公式為：

(19)

本研究所選重現(xiàn)期有4個，分別是50，100，200和1 000年。由于本研究采用不及制累積概率，因此上述4個重現(xiàn)期對應(yīng)的頻率P分別為0.980，0.990，0.995，0.999。

估計量標準誤差(SE)越小，表明估計量的有效性越好，在M次統(tǒng)計試驗中，其計算公式為：

(20)

估計量均方根誤差(RMSE)越小，表明參數(shù)估計方法的有效性越好，其計算公式為：

(21)

RMSE表示的是估計值與真實值偏差平方的算術(shù)平均值再開方。α、μ2個參數(shù)值的誤差計算原理類似，此處不再贅述。

2.2試驗方案設(shè)計與試驗結(jié)果

為不失一般性，總體分布參數(shù)值取為α=1，u=0。樣本容量用n表示，取n=50，70，100，150，200，300，700，1 000共8組方案，各組的模擬次數(shù)為1 000[21-22]。按照上述方案進行蒙特卡洛模擬。因篇幅原因試驗結(jié)果不再全部列出，現(xiàn)選取具有代表性的n=50，70，100，200時的試驗結(jié)果列于表1。

表 1　不同樣本容量時蒙特卡洛模擬的不偏性和有效性計算結(jié)果

由表1可以看出：①隨著樣本容量的增加，3種方法對應(yīng)設(shè)計值的Bias、SE和RMSE都逐漸增加；②線性矩法的設(shè)計值偏差Bias最小，矩法次之，交互熵法最大。線性矩法設(shè)計值估計具有良好的不偏性，交互熵法的設(shè)計值偏差Bias控制在7%以內(nèi)，可以滿足精度要求；③交互熵法的設(shè)計值標準誤差SE最小，線性矩法次之，矩法最差。說明利用交互熵法估計設(shè)計值優(yōu)勢明顯；④交互熵法的均方根誤差RMSE最小，線性矩法次之，矩法最差。說明將交互熵法應(yīng)用于Gumbel分布參數(shù)估計時，其有效性優(yōu)于其余2種方法。

3實例應(yīng)用

本研究選取陜西省關(guān)中地區(qū)周至、武功、蒲城、禮泉、白水、潼關(guān)6個水文站年降水量資料(表2)，經(jīng)過“三性”審查，所有資料均滿足可靠性、代表性及一致性要求。利用這些資料研究交互熵法在Gumbel分布參數(shù)估計中的可行性，并結(jié)合矩法和線性矩法進行擬合效果評價。

表 2　陜西省關(guān)中地區(qū)6個水文站年降水量資料系列長度

3.1參數(shù)估計

以矩法所得參數(shù)作為初始值，利用Matlab 2010b中的優(yōu)化函數(shù)確定最小交互熵函數(shù)最小值和分布參數(shù)，并將結(jié)果與矩法、線性矩法計算結(jié)果進行對比，結(jié)果見表3。由表3可以看出：①矩法、線性矩法及交互熵法的S(P)和Dmin(q,p)(最小交互熵函數(shù)值)依次減??；②按照Kullback最小交互熵原理，所用到的3種推求Gumbel分布參數(shù)的方法中，矩法最差，線性矩法次之，交互熵法最優(yōu)。

表 3　陜西省關(guān)中地區(qū)6個水文站年降水量分布參數(shù)及熵值計算結(jié)果

3.2年降水量頻率曲線圖繪制

根據(jù)表3的參數(shù)值，推求年降水量設(shè)計值，并繪制年降雨量頻率曲線，結(jié)果如圖1所示。由圖1可以看出：選用Gumbel分布在3種方法下的擬合效果均表現(xiàn)良好；在實測數(shù)據(jù)的中低值部分，運用矩法和線性矩法求得的分布曲線擬合效果相差不大，比較周至、武功、蒲城、白水、潼關(guān)、禮泉6個水文站的擬合效果可以看出，運用交互熵法所求降水量頻率曲線明顯較前2種方法求得的曲線更加接近實測值；周至、潼關(guān)、武功、白水、禮泉5個站的擬合曲線表明：在實測數(shù)據(jù)的高值處，線性矩法優(yōu)于矩法；蒲城站2種方法擬合效果差異不明顯；在周至、武功、蒲城、禮泉、白水、潼關(guān)6個站點，用交互熵法計算得到的分布曲線在高值處的擬合效果均明顯優(yōu)于矩法和線性矩法。

總的來說，無論在傳統(tǒng)方法擬合效果較為接近的中低值部分，還是高值部分，交互熵法所求曲線的擬合效果最優(yōu)，線性矩法次之，矩法最差。

圖 1陜西省關(guān)中地區(qū)6個水文測站年降水量頻率曲線

Fig.1Precipitation frequency at 6 stations in Central Shaanxi

3.3擬合優(yōu)度評價

應(yīng)用累積相對偏差平方和(δ)分析上述參數(shù)估計方法的擬合效果。實測值與設(shè)計值累積偏差平方和(δ)的計算公式為:

(22)

不同參數(shù)估計方法在6個水文站點δ的計算結(jié)果如表4所示。

表 4　3種參數(shù)估計方法的累積相對偏差平方和的比較

從表4可以看出，在陜西關(guān)中的6個水文站，用交互熵法(CE)所求曲線的累計相對偏差平方和(δ)最小，其次是線性矩法，矩法的累積相對偏差平方和(δ)最大，說明交互熵法擬合效果最優(yōu)，這與圖1的擬合效果一致，并且與表3計算所得的最小交互熵函數(shù)值的表現(xiàn)規(guī)律相同，表明交互熵法是一種可行的參數(shù)估計方法。

4討論與結(jié)論

1)介紹了受分位數(shù)對約束的交互熵，以及基于Kullback最小交互熵原理求Gumbel分布參數(shù)的方法。利用蒙特卡洛試驗研究交互熵法的統(tǒng)計性能，并與矩法和線性矩法相對比，結(jié)果表明交互熵法具有良好的有效性，且其不偏性滿足精度要求，證明交互熵法是一種統(tǒng)計性能良好的參數(shù)估計方法。

2)選擇Gumbel分布曲線擬合陜西省關(guān)中地區(qū)6個水文站的年降水資料，所采用的3種方法的分布曲線擬合效果均可以達到要求，設(shè)計值的累積相對偏差平方和δ也很小。因此，Gumbel分布是一種能較好擬合關(guān)中地區(qū)年降水變量的分布曲線，可在以后的工作中推廣使用。

3)選定Gumbel分布曲線的前提下，以陜西省關(guān)中地區(qū)6個水文站為例，對交互熵法、矩法和線性矩法的擬合效果和累積相對偏差平方和進行對比分析。結(jié)果表明，交互熵法是一種簡單且有效的水文分布參數(shù)估計方法。

4)交互熵法與分布類型無關(guān)，并且原理簡單，計算簡便。本研究僅是一個交互熵法進行水文分布參數(shù)估計的探索性研究，對于我國水文計算普遍使用的P-Ⅲ分布，尚有待于進一步研究。

[參考文獻]

[1]李揚.水文頻率新型計算理論與應(yīng)用研究 [D].陜西楊凌：西北農(nóng)林科技大學(xué),2013.

Li Y.Research on new theory and application of hydrologic frequency analysis [D].Yangling,Shaanxi:Northwest A&F University，2013.(in Chinese)

[2]詹道江,徐向陽,陳元芳.工程水文學(xué) [M].北京:中國水利水電出版社,2011:141.

Zhan D J,Xu X Y,Chen Y F.Engineering hydrology [M].Beijing:China Water Power Press,2011:141.(in Chinese)

[3]Rao A R,Hamed K H.Flood frequency analysis [M].Florida:CRC Press LLC,2000:73-81.

[4]桑燕芳,王棟,吳吉春.水文頻率分析中參數(shù)估計 SAGA-ML 方法的研究 [J].水文,2009,29(5):23-29.

Sang Y F,Wang D,Wu J C.Research on SAGA-ML method for parameter optimition in hydrologic frequency analysis [J].Journal of China Hydrology,2009,29(5):23-29.(in Chinese)

[5]Good I J.Probability and the weighing of evidence [J].Biomet-rika,1950,38(3):170-171.

[6]Pandey M D.Extreme quantile estimation using order statistics with minimum cross-entropy principle [J].Probabilistic Engineering Mechanics,2001,16(1):31-42.

[7]Kullback S.Information theory and statistics [M].New York:Dover Publications,INC,1959.

[8]Jaynes E T.Information theory and statistical mechanics [J].Physical Review,1957,106(4):620-630.

[9]俞禮軍,嚴海,嚴寶杰.最大熵原理在交通流統(tǒng)計分布模型中的應(yīng)用 [J].交通運輸工程學(xué)報,2001,1(3):91-94.

Yu L J,Yan H,Yan B J.Maximum entropy method and its application in probability density function of traffic flow [J].Journal of Traffic and Transportation Engineering,2001,1(3):91-94.(in Chinese)

[10]Shore J,Johnson R.Properties of cross-entropy minimization [J].IEEE Transactions on Information Theory,1981,27(4):472-482.

[11]Deng J,Pandey M D.Using partial probability weighted moments and partial maximum entropy to estimate quantiles from censored samples [J].Probabilistic Engineering Mechanics,2009,24(3):407-417.

[12]Deng J,Pandey M D,Gu D.Extreme quantile estimation from censored sample using partial cross-entropy and fractional partial probability weighted moments [J].Structural Safety,2009,31(1):43-54.

[13]Deng J,Pandey M D.Cross entropy quantile function estimation from censored samples using partial probability weighted moments [J].Journal of Hydrology,2008,363(1):18-31.

[14]Pandey M D.Minimum cross-entropy method for extreme value estimation using peaks-over-threshold data [J].Structural Safety,2001,23(4):345-363.

[15]Shore J.Minimum cross-entropy spectral analysis [J].IEEE Transactions on Acoustics Speech & Signal Processing,1981,29(2):230-237.

[16]Lind N C,Solana V.Fractile constrained entropy estimation of distributions based on scarce data [J].Civil Engineering Systems,1990,7(2):87-93.

[17]Lind N C,Hong H P,Solana V.A cross entropy method for flood frequency analysis [J].Stochastic Hydrology and Hydraulics,1989,3(3):191-202.

[18]Hosking J R.M.L-moments:analysis and estimation of distributions using linear combinations of order statistics [J].Journal Royal Statistical Society,1990,52(2):105-124.

[19]牛林森，宋松柏.交互熵法在洪水頻率分布參數(shù)估計中的應(yīng)用 [J].水資源研究，2013,2(6):389-394.

Niu L S，Song S B.An application of cross entropy method to the parameters estimation in flood frequency analysis [J].Journal of Water Resources Research,2013,2(6):389-394.(in Chinese)

[20]Bhattarai K P.Partial L-moments for the analysis of censored flood samples [J].Hydrological Sciences Journal,2004,49(5):855-868.

[21]董雙林.Gumbel 分布的參數(shù)估計方法的統(tǒng)計分析 [J].水利學(xué)報,1989(11):35-42.

Dong S L.Statistical analysis of the methods in Gumbel distribution parameter estimation [J].Journal of Hydraulic Engineering,1989(11):35-42.(in Chinese)

[22]羅純,王筑娟.Gumbel 分布參數(shù)估計及在水位資料分析中應(yīng)用 [J].應(yīng)用概率統(tǒng)計,2005,21(2):169-175.

Luo C,Wang Z J.The estimates of the parameters of Gumbel distribution and their application to the analysis of the water level data [J].Chinese Journal of Applied Probability and Stasistics,2005,21(2):169-175.(in Chinese)

Application of cross entropy method in estimation of Gumbel distribution parameters

NIU Lin-sen，SONG Song-bai

(CollegeofWaterResourcesandArchitecturalEngineering,NorthwestA&FUniversity,Yangling,Shaanxi712100,China)

Abstract：【Objective】 This paper studied the universal applicability of the cross entropy method in estimation of hydrology frequency distribution parameters.【Method】 Kullback minimum cross entropy principle was used to estimate the parameters of Gumbel distribution.Then,Monte Carlo experiments were performed to verify its statistical performance.Observed precipitations at 6 hydrologic stations,Zhouzhi,Wugong,Pucheng,Liquan,Baishui and Tongguan in Central Shaanxi,were compared with predictions and the fitting results were evaluated using cumulative square error.【Result】 Monte Carlo experiments indicated that the SE and RMSE of cross entropy method were less than those of moments method and L-moments method.The bias of cross entropy method was less than 7%.The cumulative of square errors of cross entropy method for the 6 hydrologic stations of Zhouzhi,Wugong,Pucheng,Liquan,Baishui,and Tongguan were 0.000 017 68,0.000 065 62,0.000 019 66,0.000 063 00,0.000 014 12,and 0.000 016 61,those of L-moments method were 0.000 087 62,0.000 093 55,0.000 086 52,0.000 101 39,0.000 065 15,and 0.000 069 05,while those of moments method were 0.000 108 74,0.000 125 40,0.000 092 41,0.000 127 65,0.000 085 49,and 0.000 093 57,respectively.The cross entropy method not only possessed good effectiveness and reasonable non-biasedness,but also fitted better with the observed data series than traditional methods. 【Conclusion】 Cross entropy method is a feasible method for parameter estimation,and it can effectively improve the estimation precision of Gumbel distribution parameters.

Key words:hydrological distribution;parameter estimation;cross entropy method;Gumbel distribution;fractile-pair constraints

DOI：網(wǎng)絡(luò)出版時間:2016-05-0314:0510.13207/j.cnki.jnwafu.2016.06.028

[收稿日期]2014-11-03

[基金項目]國家自然科學(xué)基金項目(51179160，50879070，50579065)；高等學(xué)校博士學(xué)科點專項科研基金項目(20110204110017)

[作者簡介]牛林森(1990-)，女，河南駐馬店人，碩士，主要從事流域水文模擬及水文預(yù)報研究。E-mail:nls1990@163.com [通信作者]宋松柏(1965-)，男，陜西咸陽人，教授，博士，主要從事水文水資源研究。E-mail:ssb6533@nwsuaf.edu.cn

[中圖分類號]P333.9

[文獻標志碼]A

[文章編號]1671-9387(2016)06-0203-07

網(wǎng)絡(luò)出版地址:http://www.cnki.net/kcms/detail/61.1390.S.20160503.1405.056.html