CN101697172A - 一种水文频率线型参数估计方法 - Google Patents

Abstract

本发明公开了一种水文频率线型参数估计方法，其将模拟退火遗传算法(SAGA)和极大似然法(ML)联合使用，建立了SAGA-ML法：即将似然函数相反数求解极小值的表达式作为目标函数，依据矩法估计参数取值范围作为约束条件，然后应用SAGA进行参数估计。与常规ML法思路有本质不同，SAGA-ML法通过遗传算法进行参数优化。通过蒙特卡罗试验，验证了SAGA-ML法在参数估计和不同频率设计值估计两个方面均具有很好的准确性；同时该方法不受线型类型、参数数目和约束条件的限制；可以避免应用常规ML法时出现似然方程无解等情况；且求解过程简便快捷，使ML法在理论上和实际应用中都成为有效的方法。

Description

一种水文频率线型参数估计方法

技术领域

本发明涉及一种水文频率分析方法，具体是一种水文频率线型参数估计方法。

背景技术

水文频率分析是研究某水文随机变量出现不同数值可能性的大小，为水利水电工程建设、水资源评价管理等提供具有概率意义的水文设计值。水文频率分析主要涉及两个问题，一是线型选择，另一个是参数估计。目前为止，国内外采用的线型达10多种(孙济良，秦大庸.水文频率分析通用模型研究[J].水利学报，1989，(4)：1-10；金光炎.水文频率分析述评[J].水科学进展，1999，10(3)：319-327)，但水文随机变量究竟服从何种概率分布还没有定论。本发明不对线型选择问题作讨论，而是以P-III型分布为例，探讨参数估计的新方法。

目前，参数估计有多种方法，主要是：①矩法(MOM)；②权函数法(WF)，包括单权函数法和双权函数法；③优化适线法(FIT)，包括最小二乘和最小一乘估计方法；④极大似然法(ML)等。此外还有近年来兴起的基于最大熵原理(POME)的参数估计法。不同方法的参数估计效果不同，有研究对各种方法的优劣进行对比研究(刘光文.皮尔逊III型分布参数估计[J].水文，1990，(4)：1-15；刘光文.皮尔逊III型分布参数估计(续完)[J].水文，1990，(5)：1-14；Singh V.P.Entropy-based parameter estimation inHydrology[M]，Kluwer Academic Publishers(Boston/London)，1998)，得到如下结论：POME与ML法效果相当，优于其他方法；相比于ML法，POME参数估计过程相对简便。综合分析可以看出，上述常用的参数估计方法主要存在以下问题：①除MOM外，几乎所有方法在求解参数时都受线型类型的影响，即对于不同分布线型需要分别推导参数求解表达式和(或)方程(组)。虽然应用MOM求解参数较为直观和方便，不因分布线型的不同而异，但其结果存在很大误差(特别是求解C_s时)，实际中一般用于初估参数值；②各种方法参数求解过程较为复杂困难。即使是相对较优的POME和ML法，一般情况下也需要解一非线性方程组，且方程中含有Digamma函数等非显式函数，造成求解困难；③各种方法参数求解的难易程度受参数个数影响很大；④虽然理论上讲，ML法的参数估计效果好且精度高，但由于常规做法是通过求偏导数估计参数，由此带来很多问题，实际应用中存在许多缺陷：一是求解过程困难；二是由于要通过求偏导数估计参数，一般认为C_s＞2时，无局部极值，似然方程无解，此时ML法无效；三是估计P-III型分布参数时，由于应用到三个一阶矩，使参数估计结果不灵敏(金光炎.水文频率计算中参数估计方法述评[J].安徽水利科技，2004，3：38-40)，等。因此使得ML法在实际应用中较少。

水文频率分析参数估计的实质就是参数优化问题。目前遗传算法(Holland J.H.Adaptation in Natural and Artificial System[M].Ann Arbor，Michigan：University ofMichigan Press，1975)是解决此类问题较优的方法。已有研究将遗传算法与适线法结合，用于水文频率分析的参数估计，其实质仍是基于最小二乘法或最小一乘法建立目标函数，因此仍不能克服适线法本身的缺陷。

而对于常规极大似然法，其主要分析思路如下：

水文随机变量服从的概率密度函数f(x，θ)确定之后，应用ML法估计参数θ，是以式(1)为基础进行求解。

$L\left(\mathrm{\θ}\right)=L(x_1,x_2,\·\·\·,x_n;\mathrm{\θ})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}f(x_i,\mathrm{\θ})---\left(1\right)$

其中L(θ)称为似然函数，x_i为随机变量值，θ为待估参数。

式(2)为极大似然估计式：

$L(x_1,x_2,\·\·\·,x_n;\stackrel{\‾}{\mathrm{\θ}})=\max \underset{i=1}{\overset{n}{\mathrm{\Π}}}f(x_i,\mathrm{\θ})---\left(2\right)$

满足式(2)的参数θ称为极大似然估计量，θ即为参数取值范围内的最优解。

因L(θ)和ln(L(θ))在同一θ处取得极值，实际中极大似然估计量θ常由式(3)估计：

$\frac{d}{\mathrm{d\θ}}\ln L\left(\mathrm{\θ}\right)=0---\left(3\right)$

式(4)为P-III型概率密度函数，α、β、a₀为三个参数。

$f\left(x\right)=\frac{{\mathrm{\β}}^{\mathrm{\α}}}{\mathrm{\Γ}\left(\mathrm{\α}\right)}{(x-a_0)}^{\mathrm{\α}-1}\·e^{-\mathrm{\β}(x-a_0)}---\left(4\right)$

应用ML法推求P-III型分布中的三个参数时，具体过程如下：

建立似然函数表达式：

$L\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}f(x_i,\mathrm{\θ})={\left[\frac{{\mathrm{\β}}^{\mathrm{\α}}}{\mathrm{\Γ}\left(\mathrm{\α}\right)}\right]}^n\×\underset{i=1}{\overset{n}{\mathrm{\Π}}}{(x_i-a_0)}^{\mathrm{\α}-1}\×\exp (-\mathrm{\β}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(x_i-a_0))---\left(5\right)$

式(5)两侧取对数得式(6)：

$\ln L\left(\mathrm{\θ}\right)=\mathrm{n\α}\ln \mathrm{\β}-n\ln \left(\mathrm{\Γ}\left(\mathrm{\α}\right)\right)+(\mathrm{\α}-1)\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\ln (x_i-a_0)-\mathrm{\β}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(x_i-a_0)---\left(6\right)$

根据式(3)建立参数估计方程组：

其中，

为Digamma函数。

根据式(7)可以求解出α、β、a₀三个参数。虽然该方法求解精度高，但由推导过程可以注意到：①由于a₀的取值范围为[0，X_min]，当在该范围内a₀无局部极大值时，则似然方程无解。一般认为C_s＞2时无解，此时只能在[0，X_min]内优选求解；②即使是C_s＜2的情况，由于式(7)的三个方程均为非线性方程，且第一个方程中含有Digamma函数，求解过程困难；③由于求解过程中使用三个一阶矩，使序列中的小值对参数的确定起主导作用，这与实际情况不符。因此，在P-III型分布参数估计时，很少用到ML法。

发明内容

本发明所要解决的技术问题是针对常规极大似然法存在的上述三个主要问题，通过将模拟退火遗传算法(SAGA)与极大似然法(ML)相结合，提供一种估计结果准确性高，稳定性较好，求解过程简便，应用范围广，适用性好的水文频率线型参数估计方法。

本发明所述的一种水文频率线型参数估计方法，包括以下步骤：

(1)依据水文物理成因机制，并结合流域和地区经验，选择合适的水文频率线型描述待分析的水文序列；

(2)根据极大似然法的基本原理，建立水文线型参数估计时对应的极大似然估计式；

(3)建立参数优化问题，其中，将极大似然估计式相反数求解极小值作为目标函数，然后依据矩法确定各参数相应的取值范围。通过建立水文线型参数估计的优化问题，有利于整体上对参数进行优化求解；

(4)应用模拟退火遗传算法求解步骤(3)所建立的参数优化问题，最终得到线型对应的参数值。引入遗传算法在参数优化估计方面的优势，可避免常规ML法出现似然方程无解等缺陷。

上述步骤2)具体过程是：

(5)首先确定水文随机变量服从的概率密度函数f(x，θ)；

(6)再建立似然函数L(θ)

$L\left(\mathrm{\θ}\right)=L(x_1,x_2,\·\·\·,x_n;\mathrm{\θ})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}f(x_i,\mathrm{\θ})---\left(1\right)$

其中，x_i为随机变量值，θ为待估参数；

则式(2)即为极大似然估计式：

$L(x_1,x_2,\·\·\·,x_n;\stackrel{\‾}{\mathrm{\θ}})=\max \underset{i=1}{\overset{n}{\mathrm{\Π}}}f(x_i,\mathrm{\θ})---\left(2\right)$

满足式(2)的参数θ称为极大似然估计量，θ即为参数取值范围内的最优解。上述步骤3)中参数优化的具体过程是：

(7)对式(1)取对数并取负，得到式(10)，此时式(10)的极小似然估计量和式(2)的极大似然估计量是一致的

$F=-\ln L\left(\mathrm{\θ}\right)=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\ln f(x_i,\mathrm{\θ})---\left(10\right)$

(8)求解式(10)的最小值优化问题，其目标函数为式(11)

$F=\min (-\ln L\left(\mathrm{\θ}\right))=\min (-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\ln f(x_i,\mathrm{\θ}))---\left(11\right)$

其参数的取值范围是：

a_i≤θ_i≤b_i  i＝1，2......n    (12)

式(11)～(12)中，x_i为随机变量值，θ_i为待估参数，a_i和b_i是由矩法确定的第i个变量的上下界，n为参数个数。

针对常规极大似然法存在的问题，本发明将模拟退火遗传算法和极大似然法结合，建立了一种参数估计新方法(SAGA-ML)。本发明相对现有技术具有以下有益效果：

(1)SAGA-ML法的参数估计值和不同频率设计值的估计结果准确性高，稳定性较好，与POME方法(Principle of Maximum Entropy)效果相当，优于其他方法。

(2)将模拟退火遗传算法与ML法结合，大大简化了传统ML法参数求解过程。

(3)与常规ML法通过求偏导数估计参数不同，SAGA-ML是应用遗传算法整体上进行参数寻优，因此避免了C_s＞2时似然方程无解的缺陷，使得ML法的应用范围扩大。

(4)此外，由于结合了遗传算法的参数估计思想，新方法不受线型分布、参数数目、约束条件等因素的影响，对不同的线型只需换作相应的目标函数(即极大似然估计式)即可，又具有很好的适用性。

(5)由于改进ML法利用模拟退火遗传算法进行参数优化，这与常规ML法通过求偏导数估计参数的思路有着本质的不同，因此它可以克服常规ML法存在的缺陷，使得ML法不单是理论上较优，同时使之更加简洁实用，在实际应用中也成为较优的方法。

具体实施方式

水文频率分析中参数估计的实质是一个参数优化问题。常规ML法是根据式(3)的一阶偏导数取零时寻找似然函数的极大值点，鉴于遗传算法在参数优化方面的独特优势，为克服和解决常规ML法存在的缺陷，与通过求偏导数估计参数的思路不同，此处将遗传算法引入求解似然函数极大值的优化问题，进而建立基于SAGA(Simulated AnnealingGenetic Algorithm)的ML法(SAGA-ML)。下面以P-III型分布的参数估计为例，对本方法做详细描述。

P-III型分布的特征值参数与分布参数之间的关系见式(8)：

$\left\{\begin{array}{c}\mathrm{\α}=4/{C_s}^2\\ \mathrm{\β}=2/C_s/{\mathrm{\σ}}_x\\ a_0=\stackrel{\‾}{x}-2{\mathrm{\σ}}_x/C_s\end{array}\right.---\left(8\right)$

将式(8)代入式(4)，可得到以特征值为参数的P-III型概率密度函数表达式：

$f\left(x\right)=\frac{{(2/C_s{\mathrm{\σ}}_x)}^{4/{C_s}^2}}{\mathrm{\Γ}(4/{C_s}^2)}{(x-\stackrel{\‾}{x}+2{\mathrm{\σ}}_x/C_s)}^{4/{C_s}^2-1}\×\exp (-2/C_s/{\mathrm{\σ}}_x\·(x-\stackrel{\‾}{x}+2{\mathrm{\σ}}_x/C_s))---\left(9\right)$

根据式(1)和式(9)建立似然函数，再取对数并取负，得到式(10)，此时式(10)的极小似然估计量和式(2)的极大似然估计量是一致的。

$F=-\ln L\left(\mathrm{\θ}\right)=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\ln f(x_i,\mathrm{\θ})---\left(10\right)$

求解式(10)最小值优化问题。其目标函数为式(11)，约束条件主要是参数的取值范围，可以用MOM首先初步估计参数，然后适当放大作为参数取值范围，最后用遗传算法进行优化求解，即可方便地确定相应特征参数值。为避免常规遗传算法的解早熟、局部寻优能力差等缺点，具体求解时采用模拟退火遗传算法(SAGA)。

$F=\min (-\ln L\left(\mathrm{\θ}\right))=\min (-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\ln f(x_i,\mathrm{\θ}))---\left(11\right)$

a_i≤θ_i≤b_i  i＝1，2......n    (12)

式(11)～(12)中，x_i为随机变量值，θ_i为待估参数，a_i和b_i是由矩法确定的第i个变量的上下界，n为参数个数。

上述参数估计方法称为基于模拟退火遗传算法的改进极大似然法(SAGA-ML法)。该方法结合了遗传算法求参的思想，因此，对于其他线型同样具有良好的适用性。对于不同类型的线型进行参数估计时，只需将式(9)换作相应的概率密度函数表达式，其余过程相同。

1.1Monte-Carlo方案设计

在选择算例时，考虑到实测水文时间序列对应的参数值以及不同频率设计值的真值未知，为便于对比不同方法参数估计效果的优劣，本发明选择模拟序列为算例，即生成不同特性的模拟序列之后再分别进行分析。

应用Monte-Carlo方法生成一系列服从P-III型分布的模拟水文序列，以比较SAGA-ML法与常用方法的优劣。生成模拟序列时，考虑了参数取值大小、不同参数组合以及序列长度等方面的影响因素，以便于综合分析和讨论。模拟序列长度取50年和100年两种情况，所有序列均值x都取100，变差系数C_υ分别取0.3、0.5、1.0三个不同值，变差系数和偏态系数比值(C_υ/C_s)取2和4两个值，共计12种模拟方案。每种方案的统计试验次数N为500次。生成模拟序列时，为消除初值的影响，每次模拟长度为550年，然后截取后50和100年序列作为研究对象。

然后分别应用矩法(MOM)、优化适线法(FIT)、权函数法(WF)、最大熵原理(POME)和本发明建立的SAGA-ML法，对各个模拟序列进行参数估计，之后再依据所估的参数值求得不同频率设计值的估计值。

为便于比较各种方法对应参数和不同频率设计值估计结果的优劣性，引入无偏性和稳定性两个判别依据。其中，选用均值作为无偏性判别的依据和标准。均值公式如下：

${\mathrm{EX}}_i=1/N\·\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{X}_{\mathrm{ij}}---\left(13\right)$

稳定性反映了该方法受随机抽样的影响程度，一般用相对均方误差描述。相对均方误差越小，表明方法的稳定性越好。公式如下：

${\mathrm{SX}}_i=\sqrt{(1/N\·\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(X_{\mathrm{ij}}-X_{i0})}^2)}/X_{i0}---\left(14\right)$

式(13～14)中，X_ij表示由第j个样本得到的第i个参数的估计值，X_i0表示第i个参数的设计值，N为统计试验次数。模拟序列参数估计和设计值估计结果分别见表1和表2。

1.2Monte-Carlo统计试验结果分析

1.2.1参数估计结果统计分析

对参数估计结果(表1)进行综合分析，可得到如下结论：

(1)各种方法对同一个方案参数估计结果的差异很大。整体上，不同方法的均值估计结果相同，但对C_υ和C_s两参数的估计结果差别较大，特别是结果稳定性的差异明显。这主要是因为C_υ和C_s两个参数涉及到高阶矩，求解过程存在误差累积现象。对C_υ估计时，MOM、WF、POME和SAGA-ML效果相近，FIT效果相对较差。对C_s的估计结果，POME和SAGA-ML两种方法效果最好，FIT次之，MOM和WF效果最差。

(2)保持C_υ和C_s两值不变，相比于50年模拟序列，当序列长度增长至100年之后，各种方法参数估计结果的准确性和稳定性都有所提高，尤其是稳定性有明显提高。这符合统计学规律，即随着统计样本数的增加，结果的均方差会逐渐减小，解的结果也逐渐趋于稳定。

(3)保持相同的模拟长度和C_υ值，随着C_s值的不断增大，各种方法的参数估计效果均相对变差。分析其原因：当C_s增大之后，序列中随机变量关于均值的不对称性增强，三阶中心矩的值明显变大，因此参数求解误差也会逐渐增大。

(4)保持模拟长度和C_υ/C_s一致，随着C_υ和C_s取值的不断增大，各种方法参数估计效果逐渐降低。造成参数估计效果逐渐降低的原因同(3)。

(5)由于SAGA-ML法是应用遗传算法进行参数优化，因此避免了C_s＞2时似然方程无解的情况。算例中对于C_s＞2的模拟序列都能很好的估计出相关参数，使得ML法的实际应用范围扩大。同时，该改进ML法也克服了常规ML法参数求解过程困难的缺陷。这样使得理论上较优的ML法在实际应用中更加适用。

(6)五种方法中，矩法(MOM)的参数估计效果最差，特别是对C_s估计的结果存在很大误差，当模拟序列的C_s值越大时，误差也越大。这主要是在求解C_s过程中应用三阶中心矩，存在误差累积造成的。实际中，该方法一般用于初步估计参数值。

1.2.2设计值估计结果统计分析

对设计值估计结果(表2)进行综合分析，可得到如下结论：

表1  模拟序列参数估计结果表

表2  模拟序列设计值估计结果表

(1)取相同的C_υ和C_s值，100年模拟序列对应的设计值估计结果优于50年模拟序列；序列长度和C_υ相同，随着C_s的逐渐增大，设计值的估计结果稳定性逐渐变差；序列长度和C_υ/C_s比值相同，随着C_υ和C_s取值的增大，设计值估计结果也逐渐变差。这三点结论与参数估计结果的结论一致，原因也相同。总体上，设计值估计结果和参数估计结果的优劣相互对应。

(2)同一方案，随着频率值的降低，对应的设计值估计结果逐渐变差，0.01％频率对应设计值的估计值精度和稳定性最差。这是由参数估计结果的误差引起的，参数估计结果的误差对小频率设计值的估计结果精度影响十分大。

(3)五种方法中，MOM对应的估计值总体上较真值偏小，其他四种方法对应结果的误差较小。其中POME和SAGA-ML两种方法的估计结果偏差接近，小于FIT和WF，且稳定性较好，因此这两种方法对应的结果最优。

1.3SAGA-ML方法分析

上述模拟统计试验显示：①SAGA-ML法的参数估计值和不同频率设计值的估计结果准确性高，稳定性较好，与POME效果相当，优于其他方法；②同时，将遗传算法与ML法结合，大大简化了传统ML法参数求解过程；③与常规ML法通过求偏导数估计参数不同，SAGA-ML是应用遗传算法整体上进行参数寻优，因此避免了C_s＞2时似然方程无解的缺陷，使得ML法的应用范围扩大；④此外，由于结合了遗传算法参数估计思想，新方法不受线型分布、参数数目、约束条件等因素的影响，对不同的线型只需换作相应的目标函数即可，又具有很好的适用性。

由于改进ML法利用模拟退火遗传算法进行参数优化，这与常规ML法通过求偏导数估计参数的思路有着本质的不同，因此它可以克服常规ML法存在的缺陷，使得ML法不单是理论上较优，同时使之更加简洁实用，在实际应用中也成为较优的方法，可根据具体问题选用之。

1.4结论

针对常用水文频率分析参数估计方法求解参数过程中的困难问题，本发明将模拟退火遗传算法与ML法相结合，建立了SAGA-ML法。通过模拟统计试验可以看出：与其他方法相比，该法的参数估计结果和不同频率设计值估计结果精度相对较高，稳定性较好，求解过程较为简单，且具有很好的适用性，因此在水文频率分析计算中，可以根据具体问题的需要，选择应用这种方法。但应该认识到，实测水文资料的审查和合理线型的选择仍是首要和关键，这样才能保证水文设计值的准确性和可靠性。

Claims (3)

1.一种水文频率线型参数估计方法，其特征在于包括以下步骤：

(1)依据水文物理成因机制，并结合流域和地区经验，选择合适的水文频率线型描述待分析的水文序列；

(2)根据极大似然法的基本原理，建立水文线型参数估计时对应的极大似然估计式；

(3)建立参数优化问题，其中，将极大似然估计式相反数求解极小值作为目标函数，然后依据矩法确定各参数相应的取值范围；

(4)应用模拟退火遗传算法求解步骤(3)所建立的参数优化问题，最终得到线型对应的参数值。

2.根据权利要求1所述的水文频率线型参数估计方法，其特征在于步骤2)具体过程是：

(5)首先确定水文随机变量服从的概率密度函数f(x，θ)；

(6)再建立似然函数

$L\left(\mathrm{\θ}\right)=L(x_1,x_2,...,x_n;\mathrm{\θ})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}f(x_i,\mathrm{\θ})---\left(1\right)$

其中L(θ)称为似然函数，x_i为随机变量值，θ为待估参数；

则式(2)即为极大似然估计式：

$L(x_1,x_2,...,x_n;\stackrel{\‾}{\mathrm{\θ}})=\max \underset{i=1}{\overset{n}{\mathrm{\Π}}}f(x_i,\mathrm{\θ})---\left(2\right)$

满足式(2)的参数θ称为极大似然估计量，θ即为参数取值范围内的最优解。

3.根据权利要求1或2所述的水文频率线型参数估计方法，其特征在于步骤3)中参数优化的具体过程是：

(7)对式(1)取对数并取负，得到式(10)，此时式(10)的极小似然估计量和式(2)的极大似然估计量是一致的

$F=-\ln L\left(\mathrm{\θ}\right)=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\ln f(x_i,\mathrm{\θ})---\left(10\right)$

(8)求解式(10)的最小值优化问题，其目标函数为式(11)

$F=\min (-\ln L\left(\mathrm{\θ}\right))=\min (-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\ln f(x_i,\mathrm{\θ}))---\left(11\right)$

其参数的取值范围是：

a_i≤θ_i≤b_i   i＝1，2......n    (12)

式(11)～(12)中，x_i为随机变量值，θ_i为待估参数，a_i和b_i是由矩法确定的第i个变量的上下界，n为参数个数。