牛林森，宋松柏

(西北农林科技大学 水利与建筑工程学院，陕西 杨凌 712100)

基于交互熵原理的水文分布参数估计研究

牛林森，宋松柏

(西北农林科技大学 水利与建筑工程学院，陕西 杨凌 712100)

[摘要]【目的】 研究交互熵法进行水文分布参数估计的普适性。【方法】 应用最小交互熵原理研究Gumbel分布参数估计，在此基础上应用蒙特卡洛试验检验交互熵法统计性能，然后结合矩法和线性矩法等传统参数估计方法，以陕西省关中地区周至、武功、蒲城、礼泉、白水、潼关6个水文站年降水序列为例，计算年降水量设计值并拟合实测值序列，利用累积相对偏差平方和评价理论年降水量频率曲线对实测值序列的拟合效果。【结果】 蒙特卡洛试验检验表明，交互熵法所求设计值的有效性指标估计量标准偏差(SE)和均方根误差(RMSE)小于矩法和线性矩法,偏差指标控制在7%以内；交互熵法估计周至、武功、蒲城、礼泉、白水、潼关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。由此可知，交互熵法不仅具有较好的有效性与合理的不偏性，而且与实测序列的拟合效果也明显优于传统方法。【结论】 交互熵法是一种可行的水文分布参数估计方法,能有效提高Gumbel分布参数的估计精度。

[关键词]水文分布；参数估计；交互熵法；Gumbel分布；分位数对约束

水文频率分析是水文学研究的重要内容，其通过研究和分析水文随机现象，揭示其中蕴含的统计规律，并对未来可能的情势做出预估，以满足水利工程规划、设计、管理以及水资源利用等工作的需要[1]。水文频率分析包括分布线型选择和参数估计[2]，其中水文分布参数估计方法研究是水文频率分析的重要内容之一。现有的水文分布参数估计方法有矩法、极大似然法、概率权重矩法、线性矩法等[3]。除矩法之外，其他方法在求解参数时都会受到线型种类的影响，即对于不同分布线型需要分别推导参数求解表达式和(或)方程(组)，计算较为复杂。矩法虽然不受分布线型影响，但是精度较差(尤其是Cs值),一般用于初估参数值[4]。寻求参数求解简便、精确度高的参数估计方法是水文频率分析领域的重要课题。交互熵(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]，并采用类似贝叶斯法中先验分布和后验分布的概念，若先验分布(亦称参考分布)已知，且约束条件一定，则候选分布中能使交互熵函数值最小的分布即为所求分布,也称后验分布[10-14]。交互熵法广泛应用于结构分析和光谱分析等领域[15-16]。1988年，Lind等[17]首次将交互熵法应用于水文领域，其以分位数对为约束条件，依据最小交互熵原理，采用Gamma分布和Gumbel分布，以加拿大2条河流的年最大洪峰流量数据为例，运用交互熵法估计分布参数并绘制出拟合效果图，然后将交互熵法所求熵值与传统的矩法和极大似然法相对比，得出了交互熵法较优的结论。美国学者Hosking于1990年提出的线性矩法是概率权重矩的线性组合，与概率权重矩法结果接近，并且具有良好的统计性能[18]。目前，国内尚缺少将交互熵法应用于水文领域的研究报道[19]。为此，本研究采用Gumbel分布作为参考分布，以矩法和线性矩法等传统方法为对比，应用蒙特卡洛试验分析交互熵法的统计性能，并以陕西省关中地区6个水文站年降水序列为例，对交互熵法应用于水文领域的普适性进行了分析，以期为陕西关中地区水利工程设计中设计洪水的准确计算提供参考。

1交互熵法原理

1.1交互熵函数

设随机变量X的参考分布函数为P(x),密度函数为p(x),同时假设后验分布为Q(x)，其对应密度函数为q(x)，则交互熵函数D的形式为：

(1)

式中:D为交互熵函数，x为随机变量X的实际观测值。

1.2约束条件

选取分位数对为约束条件。将随机变量X的观测序列按升序排列，得到新序列S={xi},i=1，2，… ，r，其中xi∈R，r为观测序列中观测值的个数。假设x是一个可能发生的未知量，那么其落在X由xi分成的r+1个子区间[x0,x1)，[x1,x2),…,[xr,xr+1)内的机率相等。根据样本规则，对应的分位数概率为i/(r+1)，假设Q(x|x1，x2，…，xr)是由序列S推导出的X的分布函数,q(x|x1，x2，…，xr)为相应密度函数。则有：

(2)

式中：r为X观测序列中观测值的个数。

1.3Kullback最小交互熵原理

给定一个参考分布函数P(x)和密度函数p(x)，在分位数x=xi上，P(x)值记为Pi。通过寻求后验分布函数Q(x|x1，x2，…，xr)，使交互熵函数为：

(3)

令交互熵函数在满足分位数对约束的条件下最小，则式(2)可写为期望值的形式：

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

(4)

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

(5)

为使交互熵最小，根据变分法中的Euler-Lagrange方程，将式(5)仅对q(x)求偏导，并令其为0，有：

(6)

整理，有：

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

(7)

其中：

(8)

当x∈Ii时，有：

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

(9)

将式(7)代入式(4)的左边，考虑x∈Ii时fi(x)=1，故有：

(10)

根据式(4)可知，式(10)等于0，故有：

(11)

式中：Pi和Pi+1分别为参考分布在x=xi和x=xi+1时的概率分布函数。因此，当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为非负数。再引入中间变量S(P)：

(15)

则有：

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

(16)

式中：Dmin(q,p)可以称为最小交互熵函数。

从式(16)可以看出，最小交互熵函数可以写成参考分布的函数，并且S(P)值越小，最小交互熵函数值越小。在参考分布曲线类型选定的条件下，最小交互熵原理可解释为最优分布的最小交互熵函数值最小。

1.4参考分布

选择Gumbel分布作为参考分布，其密度函数与分布函数分别为：

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

-∞

(17)

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

(18)

式中：α为尺度参数，μ为位置参数。

2蒙特卡洛试验

本研究利用蒙特卡洛试验进行交互熵法统计性能研究，并与传统的参数估计方法，即矩法、线性矩法的估计结果进行比较。

2.1评价标准

若估计量的偏差越接近0，表明估计量的不偏性越好，在M次统计试验中，其计算公式为：

(19)

本研究所选重现期有4个，分别是50，100，200和1 000年。由于本研究采用不及制累积概率，因此上述4个重现期对应的频率P分别为0.980，0.990，0.995，0.999。

估计量标准误差(SE)越小，表明估计量的有效性越好，在M次统计试验中，其计算公式为：

(20)

估计量均方根误差(RMSE)越小，表明参数估计方法的有效性越好，其计算公式为：

(21)

RMSE表示的是估计值与真实值偏差平方的算术平均值再开方。α、μ2个参数值的误差计算原理类似，此处不再赘述。

2.2试验方案设计与试验结果

为不失一般性，总体分布参数值取为α=1，u=0。样本容量用n表示，取n=50，70，100，150，200，300，700，1 000共8组方案，各组的模拟次数为1 000[21-22]。按照上述方案进行蒙特卡洛模拟。因篇幅原因试验结果不再全部列出，现选取具有代表性的n=50，70，100，200时的试验结果列于表1。

表 1　不同样本容量时蒙特卡洛模拟的不偏性和有效性计算结果

由表1可以看出：①随着样本容量的增加，3种方法对应设计值的Bias、SE和RMSE都逐渐增加；②线性矩法的设计值偏差Bias最小，矩法次之，交互熵法最大。线性矩法设计值估计具有良好的不偏性，交互熵法的设计值偏差Bias控制在7%以内，可以满足精度要求；③交互熵法的设计值标准误差SE最小，线性矩法次之，矩法最差。说明利用交互熵法估计设计值优势明显；④交互熵法的均方根误差RMSE最小，线性矩法次之，矩法最差。说明将交互熵法应用于Gumbel分布参数估计时，其有效性优于其余2种方法。

3实例应用

本研究选取陕西省关中地区周至、武功、蒲城、礼泉、白水、潼关6个水文站年降水量资料(表2)，经过“三性”审查，所有资料均满足可靠性、代表性及一致性要求。利用这些资料研究交互熵法在Gumbel分布参数估计中的可行性，并结合矩法和线性矩法进行拟合效果评价。

表 2　陕西省关中地区6个水文站年降水量资料系列长度

3.1参数估计

以矩法所得参数作为初始值，利用Matlab 2010b中的优化函数确定最小交互熵函数最小值和分布参数，并将结果与矩法、线性矩法计算结果进行对比，结果见表3。由表3可以看出：①矩法、线性矩法及交互熵法的S(P)和Dmin(q,p)(最小交互熵函数值)依次减小；②按照Kullback最小交互熵原理，所用到的3种推求Gumbel分布参数的方法中，矩法最差，线性矩法次之，交互熵法最优。

表 3　陕西省关中地区6个水文站年降水量分布参数及熵值计算结果

3.2年降水量频率曲线图绘制

根据表3的参数值，推求年降水量设计值，并绘制年降雨量频率曲线，结果如图1所示。由图1可以看出：选用Gumbel分布在3种方法下的拟合效果均表现良好；在实测数据的中低值部分，运用矩法和线性矩法求得的分布曲线拟合效果相差不大，比较周至、武功、蒲城、白水、潼关、礼泉6个水文站的拟合效果可以看出，运用交互熵法所求降水量频率曲线明显较前2种方法求得的曲线更加接近实测值；周至、潼关、武功、白水、礼泉5个站的拟合曲线表明：在实测数据的高值处，线性矩法优于矩法；蒲城站2种方法拟合效果差异不明显；在周至、武功、蒲城、礼泉、白水、潼关6个站点，用交互熵法计算得到的分布曲线在高值处的拟合效果均明显优于矩法和线性矩法。

总的来说，无论在传统方法拟合效果较为接近的中低值部分，还是高值部分，交互熵法所求曲线的拟合效果最优，线性矩法次之，矩法最差。

图 1陕西省关中地区6个水文测站年降水量频率曲线

Fig.1Precipitation frequency at 6 stations in Central Shaanxi

3.3拟合优度评价

应用累积相对偏差平方和(δ)分析上述参数估计方法的拟合效果。实测值与设计值累积偏差平方和(δ)的计算公式为:

(22)

不同参数估计方法在6个水文站点δ的计算结果如表4所示。

表 4　3种参数估计方法的累积相对偏差平方和的比较

从表4可以看出，在陕西关中的6个水文站，用交互熵法(CE)所求曲线的累计相对偏差平方和(δ)最小，其次是线性矩法，矩法的累积相对偏差平方和(δ)最大，说明交互熵法拟合效果最优，这与图1的拟合效果一致，并且与表3计算所得的最小交互熵函数值的表现规律相同，表明交互熵法是一种可行的参数估计方法。

4讨论与结论

1)介绍了受分位数对约束的交互熵，以及基于Kullback最小交互熵原理求Gumbel分布参数的方法。利用蒙特卡洛试验研究交互熵法的统计性能，并与矩法和线性矩法相对比，结果表明交互熵法具有良好的有效性，且其不偏性满足精度要求，证明交互熵法是一种统计性能良好的参数估计方法。

2)选择Gumbel分布曲线拟合陕西省关中地区6个水文站的年降水资料，所采用的3种方法的分布曲线拟合效果均可以达到要求，设计值的累积相对偏差平方和δ也很小。因此，Gumbel分布是一种能较好拟合关中地区年降水变量的分布曲线，可在以后的工作中推广使用。

3)选定Gumbel分布曲线的前提下，以陕西省关中地区6个水文站为例，对交互熵法、矩法和线性矩法的拟合效果和累积相对偏差平方和进行对比分析。结果表明，交互熵法是一种简单且有效的水文分布参数估计方法。

4)交互熵法与分布类型无关，并且原理简单，计算简便。本研究仅是一个交互熵法进行水文分布参数估计的探索性研究，对于我国水文计算普遍使用的P-Ⅲ分布，尚有待于进一步研究。

[参考文献]

[1]李扬.水文频率新型计算理论与应用研究 [D].陕西杨凌：西北农林科技大学,2013.

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

[2]詹道江,徐向阳,陈元芳.工程水文学 [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]桑燕芳,王栋,吴吉春.水文频率分析中参数估计 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]俞礼军,严海,严宝杰.最大熵原理在交通流统计分布模型中的应用 [J].交通运输工程学报,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]牛林森，宋松柏.交互熵法在洪水频率分布参数估计中的应用 [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 分布的参数估计方法的统计分析 [J].水利学报,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 分布参数估计及在水位资料分析中应用 [J].应用概率统计,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：网络出版时间:2016-05-0314:0510.13207/j.cnki.jnwafu.2016.06.028

[收稿日期]2014-11-03

[基金项目]国家自然科学基金项目(51179160，50879070，50579065)；高等学校博士学科点专项科研基金项目(20110204110017)

[作者简介]牛林森(1990-)，女，河南驻马店人，硕士，主要从事流域水文模拟及水文预报研究。E-mail:nls1990@163.com [通信作者]宋松柏(1965-)，男，陕西咸阳人，教授，博士，主要从事水文水资源研究。E-mail:ssb6533@nwsuaf.edu.cn

[中图分类号]P333.9

[文献标志码]A

[文章编号]1671-9387(2016)06-0203-07

网络出版地址:http://www.cnki.net/kcms/detail/61.1390.S.20160503.1405.056.html