基于Python的优化适线法在水文频率分析中的应用

雷庆文，高培强，李建林

(1.河北工程大学水利水电学院，河北省邯郸市太极路19号 056038；2.河南理工大学资源环境学院，河南省焦作市世纪路2001号 454000)

水文水利计算是为水利、土木等工程建设提供设计参考数据的必要工作。其中，水文频率计算是利用数理统计方法，通过实测水文资料，假定其分布符合某种现有的概率分布模型，然后通过参数估计得到该模型中的未知参数，从而计算出所需重现期下的水文数据。可想而知，模型参数估计的精度直接关乎到工程建设的安全和投资的合理性。

目前，我国最为常用的是皮尔逊Ⅲ型分布，其具有3个参数。水文统计中常用来估计参数的方法有适线法、权函数法、熵估计法、概率权重矩法和线性矩法等。《水利水电工程水文计算规范》[1]规定：“水文频率分布曲线统计参数可先用矩法初估，再用适线法调整确定”。目估适线法依然是现行水文行业使用较多的适线方法。然而目估适线法存在着很大的随机性，计算结果因人而异，主要原因在于缺乏一个定量分析的优度指标。为此，在1993年后，水文工作者开始推荐使用计算机优化适线法，取得了较好的应用效果。基于一个优度指标(常采用离差平方和或离差绝对值)，通过智能优化算法，寻找使得优度指标最小的参数即为参数估计值。刘力等[2]利用粒子群(PSO)优化适线法对长江宜昌水文站1974～2003年共30年的年最大流量数据进行水文频率分析；姜铁兵[3]等提出用遗传算法实现水文频率计算参数优化的新方法,给出了遗传算法的具体寻优步骤以及软件设计的新思想。

文中采用物理模型中常用的参数拟合方法Levenberg-Marquardt法来计算非线性函数的最小二乘拟合，并基于Python语言给出完整的程序实现过程和模型拟合评价。

1 P-Ⅲ型分布的相关理论

1.1 分布密度曲线

P-Ⅲ型分布的概率密度函数为

(1)

式中：x为水文随机变量；α，β，x0分别为P-Ⅲ型分布的参数。其中α>0，x>x0。

各数字特征值(均值E(x)、方差D(x)、离势系数Cv、偏态系数Cs)的计算如下：

(2)

(3)

(4)

利用水文统计中常用的矩法估计，该分布的3个参数由式(2)-式(4)计算：

(5)

(6)

(7)

结合水文变量的物理性质，从理论上讲，水文统计中应用P-Ⅲ型分布还需保证x0>0,即Cs>2Cv，而Levenberg-Marquardt法讨论的是全局最优解。在求解关于固定Cs/Cv倍比系数以及x0>0等约束优化问题时，Python提供的解决方案是在优化器接口中限定参数优化的范围来求解局部最优解。

1.2 适线法的基本原理

1.2.1 经验频率计算

1.2.2 优化适线原理

现在面临的问题就是利用概率密度函数来求解频率为pm时的上α分位数。可以列出这样一个方程：

(8)

(9)

那么方程的解为：

(10)

至此，问题便转变为利用LSM(最小二乘法)来计算gaminv这个非线性函数的拟合问题。

Levenberg-Marquardt算法是非线性最小二乘拟合的常用方法，在选择了合理的参数迭代初始值后通常能给出很好的拟合结果。

2 基于Python的优化适线法

2.1 原理

Python是一种解释型脚本语言，其被广泛应用于人工智能、科学计算和统计分析等领域。文中将利用Numpy，Scipy库进行数值分析与计算，并利用Matplotlib库绘制海森概率格纸，可视化适线结果。

(2)目标函数minδ：

(11)

Python的科学计算库(Scipy)提供了非线性拟合的Levenberg-Marquardt最小二乘法计算模块。

2.2 步骤

基于Python的优化适线程序的实现过程包括4个环节。

(1)利用Numpy库中array方法将实测的水文数据转化为一个数组对象。

(2)通过Numpy中的sort与flip函数将(1)中转化的数组降序排列，并计算经验频率值。

(3)定义含待优化参数的拟合函数：函数中的形参分别对应于自变量p(频率值)，以及待优化的未知参数(P-Ⅲ型分布的3参数)。以下是该函数定义代码。

def optimize_(p,α,β,x0):

model=scipy.stats.gamma(α,1/β)

return model.ppf(1-p)+x0

(4)调用Scipy的智能优化模块optimize的curve_fit接口，设定方法为Levenberg-Marquardt算法。该方法选择初始迭代值后，会在迭代过程中不断降低均方损失值，最终输出最优化参数。

3 实例分析

以某水文站1949～2000年间47a的年径流量序列进行水文频率分析，构建求解优化参数问题的拟合模型，并利用Python语言实现优化问题的求解，获得最优估计参数。

3.1 优化适线法

将实测水文序列输入给Python的优化器计算接口，优化器回传优化拟合后P-Ⅲ型分布的3个参数值：α=6.656；β=0.017；x0=6.958。确定概率分布的3个参数后，计算出各频率下的理论值，并将计算结果绘制于海森概率格纸上，优化适线结果见图1。

图1 P-Ⅲ型分布曲线拟合对比Fig.1 Comparison of the curve fitting of P-Ⅲ distribution

3.2 矩法估计参数

矩法估计是一种古老而直观的方法，其构造统计量的原理与方法简单，不同总体分布都可使用。

(12)

(13)

(14)

(15)

通过样本观测值(实测水文序列)，利用式(12)-式(15)计算出相关参数，并代入式(5)-式(7)同样可以估计出3个参数(见表1)，将结果绘于图1概率格纸上。

表1 P-Ⅲ分布参数估计对比Tab.1 Comparison of the estimated parameters of P-Ⅲdistribution

由图1和表1对2种不同方法计算得到的参数估计值进行对比分析，不难发现矩法估计的曲线一直处在优化适线曲线的下方，两者Cv较为接近，Cs相差较大。这进一步验证了矩法估计所得的频率曲线总是系统偏小，其中Cs偏小尤为明显的说法[4]。如此一来就会使得计算设计值偏小，从而增加了工程的风险性。

3.3 拟合优度检验

以均方根误差(RMSE)和AIC准则为评价模型拟合精度的指标，均方根误差与AIC越小，说明模型效果越好。评价指标的计算结果见表2。

表2 参数估计精度评价Tab.2 Accuracy evaluation of parameter estimation

模型的残差分布状况也是判断模型拟合效果的重要标准，依其能判断原序列是否服从相应类型的概率分布，以检验模型的分布类型及其参数的选择是否合理。模型残差通常服从正态分布，可对残差序列进行正态检验，从而判断模型的效果。

Filliben相关系数能定量分析参数估计的合理性，其中，r1,r2,r3，…，rn为标准化服从正态分布的残差，ri为残差升序统计量，理论残差计算公式为Mi=φ-1(i/(n+1))。Filliben相关系数为

(16)

Filliben相关系数越大则模型的拟合优度越好，其计算结果如图2所示。从图2可以看出残差符合正态分布。

图2 正态检验QQ图Fig.2 QQ diagram of normal test

无论从表2中各指标值，还是从图2的正态分布检验结果来看，优化适线法对P-Ⅲ型分布的拟合效果均好于矩法，因此，该水文序列符合用Levenberg-Marquardt方法推求P-Ⅲ型概率分布假设。

以上分析表明，利用计算非线性函数最小二乘拟合的Levenberg-Marquardt方法来最优化P-Ⅲ型分布参数是合理的，计算结果可以满足水文频率计算的要求。

3.4 带约束的参数优化问题

(1)利用Levenberg-Marquardt最小二乘法求解的是全局最优拟合问题。但是，实际水文频率计算中要确保设计值不出现小于0的不合实际值，则需约束参数x0>0，因此，需要在此条件下寻找最优参数估计值。优化器curve_fit接口可指定参数bounds中的变量x0的范围为(0，+∞)，此情况下曲线优化拟合采用trf方法。

(2)实际计算中，一些地区给出了当地水文频率计算的Cs/Cv等值线图，因此，常需要结合经验数据与实际情况固定Cs/Cv倍比系数k。此时，通过数学计算，将参数x0进行变换。因Cs/Cv=k，由式(5)-式(7)得

k的取值应参照当地的水利计算手册,如河北省年降水量参数等值线图的编制及分析[8]中指出Cs/Cv分区图范围为2.0～3.5，便可将bounds中k的值设定在此区间，实现此区间内的最优化计算。

结合某水文站1949～2000年间47a的平均流量序列进行水文频率分析。限定Cs/Cv的范围在[2.5,4.5]区间，在此约束条件下计算使优度指标最小的参数，进而推求设计值。其结果对比如表3、图3所示。

图3 有约束条件优化对比Fig.3 Comparison of constrained conditions

表3 有约束条件指标对比Tab.3 Comparison of index under constrained conditions

单从优度指标分析，全局优化的拟合度肯定比局部优化的拟合度要好。但是，P-Ⅲ型分布模型仅是个统计模型，并不具备物理意义，所以，实际水利计算中有时须考虑当地的经验数据和实际情况。虽然这样会使得拟合效果变差，但却与实际更为相符，因此，究竟采用全局最优解还是局部最优解也是工程设计中需要考虑的，设计者可根据工程风险等来选择使用局部最优还是全局最优的设计值。

4 结论

(1)采用非线性函数曲线拟合的思想，可将P-Ⅲ型分布参数估计问题转化为一个求解Γ分布逆函数的非线性拟合问题，并借助于物理模型参数拟合中常用的Levenberg-Marquardt最小二乘法来拟合该非线性函数。

(2)通过与矩法估计获得参数进行比较，无论从拟合图像还是AIC和RMSE指标来看，优化适线法都要明显优于矩法估计。

(3)通过残差正态检验QQ图以及Filliben相关系数均表明利用Levenberg-Marquardt最小二乘法的优化适线法估计参数是合理性的。

综上，文中提出的基于Python实现P-Ⅲ型分布的优化适线法避免了传统目估适线法的盲目性和不确定性。相较于遗传算法和粒子群算法的优化适线方法，虽然其本质都是求解使得优度指标最小的参数组合，但本文中采用的方法更为简单方便。通过将问题转变为利用Python的优化拟合工具解决非线性函数的拟合问题，避免了繁琐的编程与参数调整过程，更加适合工程师在工程设计中使用。

在文中提出的优化方法中还存在不足，如未考虑具有历史洪水的计算问题，只对矩法估计做了比较。虽然存在一些不足，但文中的思路和方法在水文频率计算中具有较好的借鉴和应用意义。