根据IPCC（政府间气候变化专门委员会）第5次评估报告，由人类活动引起的全球气候变暖毋庸置疑，随着全球变暖的持续发展，极端天气事件的强度和频率进一步增强，由此引发多种自然灾害，对社会经济生活和人民生命财产安全带来严重威胁。因此，准确预测不同重现期的气象水文要素极值是预防自然灾害、控制和降低灾害损失的重要基础性工作^[1–5]。

目前，常用的单变量气象水文频率分析模型包括皮尔逊Ⅲ分布、耿贝尔分布、广义极值分布等。相关应用研究有如：王银堂等^[6]构造了年径流量序列的水文频率分析模型，采用广义极大似然法估计模型参数；王芝兰等^[7]研究了广义极值分布理论在气象干旱检测中的适用性；张容焱等^[8]采用泊松–耿贝尔联合极值分布模型预测了沿海各气象站TC影响大风的多年一遇工程的设计风速，基于耿贝尔分布的北京夏季高温分布图研究^[9]。此外，干旱、洪水等多变量水文频率分析中的边缘分布模拟研究也属于单变量水文频率分析^[10–12]。以上研究用一些经典的分布函数如皮尔逊Ⅲ分布、耿贝尔分布、广义极值分布模拟年最大洪水量级、年最大洪峰、极端高水位、极端高温等气象水文要素的极值概率分布，常采用矩法估计、极大似然估计或权函数等传统参数估计法估计极值分布函数的未知参数，但这些参数估计方法对资料的长度要求很高，需要大量样本资料，对于中小样本可能导致误差增大、估计值不稳定等^[13–14]。故当数据样本稀少时如何估计极值水文频率分析模型的参数已经得到广泛关注^[15]，但针对小样本地区的极值水文频率分析模型应用研究尚不多见。

因此，以耿贝尔分布模型为例，提出一种新的参数估计方法—最大熵估计，用以解决极值水文频率分析模型在资料稀少地区模型参数难以估计的问题，并且通过实例研究和大量模拟试验证明最大熵估计的可靠性和鲁棒性。

1 基于最大熵估计的耿贝尔分布模型 1.1 耿贝尔分布理论基础

Fisher–TippetⅠ型分布，即耿贝尔分布常被用于水文统计中，耿贝尔分布函数和密度函数分别为：

$F\left( x \right) = {{\rm{e}}^{\left\{ { - {{\rm{e}}^{ - \alpha (x - \beta )}}} \right\}}}$

(1)

$f\left( x \right) = \alpha {{\rm{e}}^{ - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}}}}$

(2)

而 ${x_{\max }}$ 超过 $x$ 的概率为：

$ G\left( x \right) = P\left( {{x_{\max }} \ge x} \right) = 1 - F\left( x \right) = 1 - {{\rm{e}}^{ - {{\rm{e}}^{\alpha \left( {x - \beta } \right)}}}} $

(3)

式（1）～（3）中， $\alpha $ 和 $\,\beta $ 为未知参数。

1.2 参数估计

根据文献[16]，矩法和极大似然估计法常用来估计耿贝尔分布的未知参数。 $\alpha $ 和 $\,\beta $ 的计算表达式如下：

$\left\{\!\!\!\!\begin{array}{l} \alpha = \dfrac{{\text{π}} }{{\sqrt 6 }}\dfrac{1}{{{\sigma _x}}}, \\ \beta = E\left( x \right) - 0.45{\sigma _x} \\ \end{array}\right. $

(4)

根据矩法估计的原理^[13]，将 $E\left( x \right)$ 和 ${\sigma _x}$ 的矩估计量 $\overline x$ 和 ${s_x}$ 代入式（4）即可获得参数 $\alpha $ 和 $\,\beta $ 的矩法估计值， $\overline x$ 和 ${s_x}$ 的计算表达式如下：

$\left\{\!\!\!\!\begin{array}{l} \overline x = \dfrac{1}{n}\displaystyle \sum\limits_{i = 1}^n {{x_i}} ,\\ {s_x} = \sqrt {\dfrac{1}{{n - 1}}{{\left( {{x_i} - \overline x} \right)}^2}} \\ \end{array}\right. $

(5)

根据文献[13]，矩法估计的好坏取决于样本点的个数，样本量越多，参数估计越准确。

极大似然估计的原理如下：设 ${x_1}, {x_2}, \cdots , x{}_n$ 为相应于样本 ${X_1}, {X_2}, \cdots , X{}_n$ 的一个样本值，耿贝尔分布的似然函数为：

$\begin{gathered} L\!\left( {\alpha , \beta } \right)\! \!=\! \!L\left( {{x_1}, {x_2}, \cdots , {x_n};\alpha , \beta } \right) \! = \! \! \prod\limits_{i = 1}^n \!{\alpha {{\rm{e}}^{ - \alpha \left( {{x_i} \!- \!\beta } \right) \!-\! {{\rm{e}}^{ - \alpha \left( {{x_i}\! -\! \beta } \right)}}}}} \\ \end{gathered} $

(6)

若

$L\left( {{x_1}, {x_2}, \cdots , {x_n};\hat \alpha , \hat \beta } \right) \!=\! \max \left\{ {\prod\limits_{i = 1}^n {\alpha {{\rm{e}}^{ - \alpha \left( {{x_i} - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {{x_i} - \beta } \right)}}}}} } \right\}$

(7)

则称 $\hat \alpha \left( {{x_1}, {x_2}, \cdots , {x_n}} \right)$ 和 $\,\hat \beta \left( {{x_1}, {x_2}, \cdots , {x_n}} \right)$ 为 $\alpha $ 和 $\,\beta $ 的极大似然估计量，求解方程如下：

$\left\{ \begin{gathered} \frac{{{\rm{d}}\ln \;L\left( {\alpha , \beta } \right)}}{{{\rm{d}}\alpha }}=\frac{{\ln \left( {\displaystyle\prod\limits_{i = 1}^n {\alpha {{\rm{e}}^{ - \alpha \left( {{x_i} - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {{x_i} - \beta } \right)}}}}} } \right)}}{{{\rm{d}}\alpha }} = 0 ,\\ \frac{{{\rm{d}}\ln \;L\left( {\alpha , \beta } \right)}}{{{\rm{d}}\beta }} = \frac{{\ln \left( {\displaystyle\prod\limits_{i = 1}^n {\alpha {{\rm{e}}^{ - \alpha \left( {{x_i} - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {{x_i} - \beta } \right)}}}}} } \right)}}{{{\rm{d}}\beta }} = 0 \\ \end{gathered} \right.$

(8)

方程组（8）的求解同样需要大量样本数据，为此作者提出一种新的参数估计法，即最大熵估计法，其优点是不需要大量观测数据 ${x_i}\left( {i = 1, \;2,\; \cdots ,\; n} \right)$ ，只需要两个数据，即 ${x_i}$ 的最小值和最大值，数学原理如下：

对于连续信源的随机试验而言，耿贝尔分布熵的定义如下：

$\begin{aligned}[b] &\qquad H\left( {f\left( x \right)} \right) = - \int_c^d {f\left( x \right)\ln f\left( x \right){\rm{d}}x} =\\ & \! \! \! \!- \alpha \int_c^d {{{\rm{e}}^{ - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}}}} \cdot \left[ {\ln\; \alpha - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}}} \right]{\rm{d}}x} \end{aligned} $

(9)

式中， $c$ 和 $d$ 为气象水文要素 $x$ 的最小值和最大值，由于无法获取真实的最小值和最大值，在应用时一般取该气象水文要素观测样本的最小值和最大值。文献[17–18]指出，可以将最大熵原理作为求解模型参数的依据，因为其符合最大熵原理、熵增原理、第一原理、最大多重性原理和一致性要求。根据最大熵原理，当 $H\left( {f\left( x \right)} \right)$ 取得最大值时，得到的参数应是最优的。

基于以上分析，可建立如下的最优化模型：

$\begin{aligned}[b] \max H\left( {f\left( x \right)} \right) = - \alpha \displaystyle\int_c^d {{{\rm{e}}^{ - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}}}} \cdot } \\ \qquad\qquad\qquad\;\, \left\{ {\ln\; \alpha - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}}} \right\}{\rm{d}}x \end{aligned}$

(10)

根据多元函数极值理论，可得：

$\left\{\! \!\!\!\begin{array}{l} \dfrac{{\partial H\left( {f\left( x \right)} \right)}}{{\partial \alpha }} = \displaystyle\int_c^d {{{\rm{e}}^{ - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}}}} \cdot } \\ \qquad\left\{ {\ln\;\alpha - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}}} \right\}{\rm{d}}x + \\ \qquad\displaystyle\int_c^d \left\{ {\alpha {{\rm{e}}^{ - \alpha \left( {x \!-\! \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}}}} \cdot \left( {x\! -\! \beta } \right) \cdot \left[ {{{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}} \!- \!1} \right] \cdot } \right.\\ \qquad\left.\left[ {1 + \ln\; \alpha - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha - \alpha \left( {x - \beta } \right)}}} \right]+\right.\\ \qquad\left.{{\rm{e}}^{ - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x - \beta } \right)}}}} \right\}{\rm{d}}x = 0, \\ \dfrac{{\partial H\left( {f\left( x \right)} \right)}}{{\partial \beta }} \!=\! {\alpha ^2}\displaystyle\int_c^d {{{\rm{e}}^{ - \alpha \left( {x \!-\! \beta } \right) - {{\rm{e}}^{ - \alpha \left( {x \!-\! \beta } \right)}}}} \cdot \left[ {1 \!-\! {{\rm{e}}^{ - \alpha \left( {x \!-\! \beta } \right)}}} \right] \cdot } \\ \qquad\left[ {1 + \ln\; \alpha - \alpha \left( {x - \beta } \right) - {{\rm{e}}^{ - \alpha - \alpha \left( {x - \beta } \right)}}} \right]{\rm{d}}x = 0 \end{array} \right.$

(11)

根据式（11）可以看出， $\alpha $ 和 $\;\beta $ 的求解只需要提供随机变量 $x$ 的最小值 $c$ 和最大值 $d$ 。然而上述方程组的求解非常困难，可以利用遗传算法或MATLAB中的非线性优化函数进行寻优，本文利用MATLAB中的非线性优化函数对参数进行寻优。

1.3 拟合检验 1.3.1 拟合效果

将变量 $x$ 的实测资料从大到小排列，设排列后的数据序列为 ${x_1}, {x_2}, \cdots , {x_n}$ ，则经验频率^[16]为：

${H_i} = p\left( {X \le {x_i}} \right) = \frac{i}{{n + 1}}$

(12)

给定一系列的保证率 $p$ ，则对应的 ${x_p}$ 计算公式为^[17]：

${x_p} = \left( {\varphi {C_{\rm v}} + 1} \right)E\left( x \right)$

(13)

式中： $\varphi = - 0.779\;7\left\{ {0.577\;22 + \ln \left[ { - \ln \left( {1 - p} \right)} \right]} \right\}$ ；矩法估计中 $E\left( x \right)$ 和 ${C_{\rm v}}$ 为样本均值 $\overline x$ 和样本标准差 ${s_x}$ ，极大似然估计和最大熵估计中 $E\left( x \right)$ 和 ${C_{\rm v}}$ 的计算方法如下：

分别将极大似然法和最大熵估计法得到的参数代入式（4）即可求出 $E\left( x \right)$ 和 ${\sigma _x}$ ，再通过 ${C_{\rm v}} = \displaystyle\frac{{{\sigma _x}}}{{E\left( x \right)}}$ 即可求出 ${C_{\rm v}}$ 。

分别绘制理论频率曲线和经验频率曲线，可初步判断理论频率曲线与经验点据的拟合效果。经验频率和理论频率之间的均方误差（RMSE）计算表达式如下：

$RMS\!\!E = \sqrt {\frac{1}{n}{{\sum\limits_{i = 1}^n {\left( {{H_i} - {Q_i}} \right)^2} }}} $

(14)

式中， ${H_i}$ 和 ${Q_i}$ 分别为第 $i$ 个实测样本的经验频率和理论频率。

1.3.2 分布拟合检验

一个样本Kolmogorov–Smirnov检验（K–S检验）^[13]是用来检验样本是否来自正态分布、均匀分布和泊松分布的总体，是一种广泛使用的拟合优度检验方法。首先将矩法估计、极大似然估计和最大熵估计得到的参数和 ${x_1}, {x_2}, \cdots , {x_n}$ 分别代入式（3），得到3种参数估计方法的理论频率序列 $\left\{ {{f_1}\left( {i = 1, 2, \cdots , n} \right)} \right\}$ 、 $\left\{ {{f_2}} \right.\left( {i = } \right.$ $\left. {\left. {1,2, \cdots ,n} \right)} \right\}$ 和 $\left\{ {{f_3}\left( {i = 1, 2, \cdots , n} \right)} \right\}$ ，K–S检验用来判断这3个序列是否服从均匀分布。K–S双侧检验的原假设 ${H_0}$ 为对所有的 $x$ 值 $F\left( x \right) = F\left( {{x_0}} \right)$ ，则备择假设 ${H_1}$ 为至少存在一个 $x$ ，使得 $F\left( x \right) \ne F\left( {{x_0}} \right)$ 。其检验原理如下^[15]：

设 $S\!\!\left( x \right)$ 表示1组数据的经验分布，定义1组随机样本 ${x_1}, {x_2}, \cdots , {x_n}$ 的经验分布函数，用阶梯函数表示：

$S\!\!\left( x \right) = \frac{{{x_i} \le x}{\text{的个数}}}{n}$

(15)

式中， $S\!\!\left( x \right)$ 为对总体分布 $F\left( x \right)$ 的一个估计，检验统计量为：

$D = \mathop {\sup }\limits_x \left| {S\!\!\left( x \right) - {F_0}\left( x \right)} \right|$

(16)

在实际情形中，如果有 $n$ 个观测值，则可用式（17）的统计量代替 $D$ ：

$D \!=\! \mathop {\max }\limits_{1 \le i \le n} \left\{ {\max \left( {\left| {S\!\!\left( {{x_i}} \right)\! -\! {F_0}\left( {{x_i}} \right)} \right|, \left| {S\!\!\left( {{x_{i - 1}}} \right)\! -\! {F_0}\left( {{x_i}} \right)} \right|} \right)} \right\}$

(17)

当 $n \to \infty $ 时，大样本的渐进公式为：

$P ( {\sqrt n {D_n} < x} ) \to K\left( x \right)$

(18)

其分布函数的表达式为：

$K\left( x \right) = \left\{\!\!\!\! \begin{array}{l} 0 ,\;x < 0 {\text{；}}\\ \displaystyle \sum\limits_{j = - \infty }^\infty {{{\left( { - 1} \right)}^j}\exp \left( { - 2{j^2}{x^2}} \right),\;x > 0} \\ \end{array} \right.$

(19)

2 结果和分析 2.1 研究区概况及数据来源

本文选取黄河流域固原（经度106°16′，纬度36°）、平凉（经度106°40′，纬度35°33′）、定边（经度107°35′，纬度37°35′）和吴旗（经度108°10′，纬度36°55′）4个气象站点每天24 h累积降水量（数据来源于中国气象数据共享网），每个站点数据序列长度不一致（其中，固原站、定边站和吴旗站均为1957—2013年，平凉站为1951—2013年），从中筛选4个站点年最大日降水量序列。

2.2 参数估计

将固原、平凉、定边和吴旗4个气象站点年最大日降水量序列分别代入式（5）、（4）和（8），可得 $\alpha $ 和 $\,\beta $ 的矩法估计量和极大似然估计值；将4个站点年最大日降水量序列的最小值和最大值代入式（10），可得 $\alpha $ 和 $\,\beta $ 的最大熵估计值，如表1所示。

由表1可知，固原站、平凉站和定边站最大熵估计得到的参数值和矩法估计、极大似然估计结果比较接近，吴旗站最大熵估计得到的 $\,\beta $ 值和矩法估计、极大似然估计法的结果有点差异， $\alpha $ 值较接近。

2.3 拟合检验

图1为4个站点最大日降水量的经验点据和各种估计方法得到的耿贝尔分布理论频率曲线。根据式（3）、（12）和（14）可得4个站点理论频率和经验频率之间的均方误差及K–S拟合检验结果如表2所示。

由图1（a）可知，3种参数估计方法得到的分布曲线均与经验点据拟合较好，矩法估计和极大似然估计的拟合效果略优于最大熵估计，总体来说三者拟合效果较为接近。由图1（b）可知，极大似然估计得到的理论频率曲线与经验点据拟合最好，最大熵估计次之，矩法估计得到的理论频率曲线与经验频率曲线拟合最差，尤其当频率处于0到20%时，矩法估计得到的理论频率曲线与经验频率曲线有较大差异。

由图1（c）可知，矩法估计得到的理论频率曲线与经验点据拟合最好，极大似然估计次之，最大熵估计的拟合效果虽然不如矩法估计和极大似然估计，但是三者的拟合效果差别不大。由图1（d）可知，3种参数估计方法得到的分布曲线均与经验点据拟合较好，拟合效果几乎一致，总的来说，矩法估计和极大似然估计的拟合效果略优于最大熵估计。

由表2可知，3种参数估计方法得到的理论频率与经验频率之间的均方误差均比较小，其中平凉站最大熵估计的均方误差最小。由表2可知，矩法估计、极大似然估计和最大熵估计在4个站点的检验概率均小于0.05，故无足够证据推翻这3种方法得到的理论频率序列均服从均匀分布的假设，即3种方法在4个站点的拟合效果均比较满意。综合分析图1、表1和2可知，最大熵估计和矩法估计、极大似然估计的拟合效果接近，然而，最大熵估计只用了2个数据，即气象水文变量的最小值和最大值，而矩法估计和极大似然估计却用了一定数量的样本（均大于50）。

3 讨　论

虽然本文提出的基于最大熵估计的耿贝尔极值水文频率分析模型只需要气象水文变量的最小值和最大值即可获得与矩法估计和极大似然估计类似的拟合效果，但是当数据序列长度较小时，最大熵估计的拟合效果是否稳健？如果只能获得气象水文变量观测样本的最大最小值，而样本的最大最小值不一定能代表总体的最大最小值，尤其当样本数量很小时，样本的最大最小值和真实最大最小值之间的误差可能较大，最大熵估计的拟合效果是否仍然可以接受，仍需要评价最大熵估计在以上两种情形下的拟合效果。

3.1 采用不同数目样本

为比较矩法估计、极大似然估计和最大熵估计在不同样本长度时的拟合效果，样本长度分别为45、35、25、15、10和5，每种数目的样本随机抽取4次，共进行24次模拟实验。将矩法估计、极大似然估计和最大熵估计在这些样本长度下得到的参数值与其在样本数目为55时得到的参数进行比较，计算3种参数估计方法的相对误差，如表3和图2所示。其中最大熵估计用到的最大最小值为所取样本的最小最大值。由于每种数目样本随机取4次，每次所取样本的最大最小值可能会不同，将其与真实最小最大值进行比较（以样本为55时的最小最大值为真实的最小最大值），每种样本条件下选取最大相对误差，不同样本量条件下最小最大值与真实最小最大值之间的最大误差如表4所示。

由表3及图2可知，当样本长度大于25时，3种参数估计方法得到的参数相对误差及平均误差均比较接近。当样本长度小于15时，最大熵估计的平均相对误差明显低于矩法估计和极大似然估计。由图2可以进一步看出：当样本长度小于15时，最大熵估计对于样本长度变化并不敏感，最大平均相对误差不超过14%，极大似然估计和矩法估计则非常敏感；当只有5个样本时，极大似然估计的平均相对误差甚至达到166%。总的来说，当样本长度大于25时，3种参数估计的拟合效果几乎一致；当样本长度小于15时，最大熵估计表现出非常大的优越性，极大似然估计的拟合效果最差。

由表4可知，当样本长度小于15时，最小值与真实最小值之间的最大相对误差达到50.7%，最大值与真实最大值之间的最大相对误差达到33.8%（需要强调的是，最小值和最大值的最大相对误差并未同时达到），而最大熵估计得到的参数的最大平均相对误差低于14%。

3.2 最小最大值准确性的影响

为进一步评价最小最大值准确性对最大熵估计的影响，以样本为55时的最小最大值为真实的最小最大值，分别考虑当最小值与真实最小值之间的相对误差为25%、50%和75%时，最大值与真实最大值之间的相对误差为25%、50%和75%时，最小最大值同时与真实最小最大值之间的相对误差为25%、50%和75%时最大熵估计得到的参数估计值与实际值之间的相对误差，共9次模拟实验，如表5～7所示。

由表5可知，最大熵估计对最小值准确性的敏感性小，当最小值与真实最小值之间相对误差为75%时，参数值平均相对误差仅为28.9%。由表6可知，最大熵估计对最大值准确性较敏感，当最大值与真实最大值之间相对误差为75%时，参数值平均相对误差达46.4%，但是参数值准确率仍接近60%。由表7可知：当最小最大值同时与真实最小最大值之间相对误差为25%时，参数值平均误差仅为21.6%，参数准确率接近80%；当最小最大值同时与真实最小最大值之间相对误差为50%时，参数值准确率接近60%，基本合格；当最小最大值同时与真实最小最大值之间相对误差为75%时，平均误差仅为57%。由表4可知，矩法估计和极大似然估计得到的参数值最大平均误差分别高达113%和166%，此时最小最大值与真实最小最大值之间的最大误差仅为50.7%和33.8%，还未同时达到，进一步说明最大熵估计显著优于矩法估计和极大似然估计。

4 结　论

针对小样本观测资料条件下无法使用传统参数估计法的问题，以耿贝尔分布为例，提出一种新的参数估计方法，即最大熵估计，该方法只需要气象水变量的最小值和最大值。以黄河河流域4个站点的最大日降水量的水文频率分析为例，验证最大熵估计的效果，结果表明最大熵估计的拟合效果与传统参数估计方法几乎一样，而传统参数估计方法需要大量数据。为验证最大熵估计在小样本条件下的拟合效果，共进行33次模拟实验。实验表明最大熵估计具有如下潜力：1）当样本长度大于25时，3种参数估计方法的拟合效果几乎一致；当样本长度小于15时，最大熵估计表现出非常大的优越性，极大似然估计的拟合效果最差。2）最大熵估计对最小值准确性的敏感性小，对最大值准确性较敏感。

除了耿贝尔分布模型外，水文频率分析模型还有很多，如皮尔逊Ⅲ分布、广义伽马分布、对数伽马分布等，下一步作者将进一步研究这些模型在小样本观测资料条件下的参数估计方法。