皮尔逊Ⅲ型曲线参数优化估计方法对比研究

2019-01-16李庆昕

水利规划与设计 2018年12期

李庆昕

(辽宁省沈阳水文局，辽宁沈阳 110094)

水文学者针对P-Ⅲ型分布函数的求解做了大量的研究和分析，并取得了一定的成果，如李世才等[1]通过利用伽马分布为数和函数可实现Kp值的转化求解，并基于此给出了伽马分布通用算法[2]的截断误差和分析的表达式。吴明官等[3]提出了一种计算速度大于麦克劳林法、分布积分和变量代换法的新的切比雪夫快速算法。刘俊哲等[4]利用龙贝格积分法计算快速、稳定性好和便于操作的优点，并结合不完全伽马分布函数的分布积分特征提出了一种计算速度高于自适应辛普森法和切比雪夫多项式法的计算方法；任柏帜等[5]基于分布函数展开式的基本形式提出了连分式展开式和P-Ⅲ型函数的通用计算法；刘仕平等[6]通过选取误差作为输入参数提出了伽马分布函数的通用表达式。对于适应性不长的参数利用该方法其计算结果较易溢出。王文川等[7]通过对水文频率曲线参数进行寻优计算和结果拟合估计验证了该方法的适用性和可靠性。

对皮尔逊Ⅲ型分布最优的参数值通常采用建立目标函数值，并进行参数寻优方法进行求解。横坐标离差和是优化算法的适线准则，对最优的参数值可通过寻优计算进行求解。本研究对P-Ⅲ型曲线的参数分别利用自适应误差分析法、龙贝格、辛普森以及梯形算法进行寻优计算，并进行分析研究，以期为提高水文频率曲线参数寻优计算运算速度和计算精度提供一定的理论支持。

1 基本理论

1.1 理论频率

利用下式表征皮尔逊-Ⅲ型分布特征。

(1)

为便于计算需进行格式统一和积分换元计算，分别引入t=β(x-a0)，dx=dt/β，若x=xp，则满足t=β(xp-a0)=u，代入上述计算公式可有：

(2)

1.2 经验频率

对于一般洪水可利用下述公式进行洪水频率求解。

(3)

式中，m、n—由大到小对水文变量排序并按自然数排序变出的序号和样本容量。

对于特大洪水可利用下述公式进行经验频率的计算。

(4)

其余量(n-l)在实测系列中的经验频率可利用下式进行计算。

(5)

1.3 适线准则

对水文序列值所对应的频率可利用上述伽马函数的数值积分法进行求解并以此可构建目标函数，表达式为：

(6)

式中，n、Wm—样本容量和权重系数；pm、pi—样本经验频率值和样本m的理论频率值；k—幂级此，本研究取k值取2。

2 数值积分法

自适应误差积分法、梯形法、龙贝格法以及辛普森法等是数值积分法的主要方法。

2.1 梯形算法

(7)

利用量算对式(7)中xi=iΔx;i=0，1，2，…，n进行求解，且式(7)成为复合梯形公式。

2.2 辛普森算法

通过将积分区间划分为若干偶数个子区间，并引入Δx=[b-a]/2n作为子区间步长计算步长，点y0，y1，y2在曲线上对应的坐标为x0，x1，x2，由此，利用曲线上三个相邻的作抛物线并对曲线进行拟合求解，将经过x0，x1，x2三个点形成的曲线段利用抛物线进行替代，利用牛顿-莱布尼茨公式对抛物线进行计算求解，其面积积分求解公式如下。

(8)

式中，S0—[a，b]区间两端点的纵坐标之和；S0、S1—奇数项纵坐标和偶数项纵坐标之和；n—偶数。

当n为偶数时可对[a，b]区间在函数f1(x)和f2(x)所夹面积利用辛普森公式进行求解，其计算公式如下。

(9)

利用辛普森法是对弧线利用抛物线进行替代进而能够更好的接近曲线的弯曲程度并减低积分误差。

2.3 Romberg算法

(10)

其计算步骤和相应公式分别如下所示：

(2)对曲线梯形值进行求解，其表达式为：

2.4 自适应误差积分法

如果α≠2条件成立，对实际步长可利用泰勒级数展开法进行近似求解，计算公式如下:

hp=C0+C1(t-α)+C2(t-α)

(11)

式中，α、t—参数和自变量；C0、C1、C2—待定参数，当α取值不同时其计算方法不同。

若α>2，则C0、C1、C2计算公式分别如下：

若α<2，则C0、C1、C2计算公式分别如下：

利用步长计算公式并分别代入C0、C1、C2值可对步长进行求解。

若α=2则可利用基本步长公式对函数曲线进行面积积分计算，并选取

(12)

当运行计算时若α>>180则数据溢出并停止计算，利用标准正态分布近似表征伽马函数分布的状态，即呈N(0，1)=φ(u)，可利用下述公式代表标准正态分布。

(13)

3 寻优建模

SSO优化算法的计算步骤：设定参数可进行搜索寻优的空间维度为n，雌性和雄性蜘蛛的个体数分别为Nf和Nm，则整个蜘蛛群体内蜘蛛个体数目总量为N，且Nf和Nm满足下式：

Nf=floor[(0.9-rand×0.25)N];
N=Nm+Nf

(14)

式中，rand—0～1之间的随机数字；floor—实数和整数之间的映射关系。

设定S蜘蛛种群中有N个蜘蛛个体，初始雌性蜘蛛随机样本F={f1，f2，…，fNf}，初始雄性蜘蛛随机样本F={m1，m2，…，mNm}，通过下述计算公式设定交配范围r。蜘蛛群体S={s1=f1，s2=f2，…，sNf=fNf，sNf+1=m1，sNf+2=m2，…，sN=mNm}。

(15)

群体内各蜘蛛的自身重量的计算如下：

(16)

式中，J(si)—通过计算所得到蜘蛛个体的目标函数适应值；bests—max(J(si)worsts=minJ(si))。

蜘蛛群体中根据蜘蛛个体的重量进行交配概率的确定，蜘蛛个体的重量越大则其繁衍后代的概率越大，并采用轮盘赌法对教派概率psi进行确定，计算公式如下：

(17)

通过上述计算过程，可分别对不同的参数组合下的参数进行优化判定，分别对龙贝格算法、紫志英积分法、梯形算法以及辛普森法的各个序列值进行理论频率值的计算；若满足收敛性则计算终止，否则重新返回并计算，直至满足要求。

4 数据检验

4.1 理想数据检验

本研究结合刘光文等相关研究成果，依据理想数据系列相关方法和理论进行统计试验分析，其中判别标准为x0.01%与x0.1%之间的相对误差低于12%～16.8%为基本依据。对曲线的参数分别利用上述方法进行数值积分求解，并将已知序列与理想样本进行对比分析，计算对比结果见表1。

表1 理想数据参数计算结果

表2 在不同水文序列中各计算方法的计算结果值

4.2 讨论分析

P-Ⅲ型曲线利用梯形积分算法的x0.1%决定平均误差值为-16.62%，其中在合格区间范围内的组数为15组，17组数据的误差值最大为-31.68%；而x0.01%的误差平均值为-17.36%，在合格误差范围内的有12组，其中误差最大值为34.62%。数值设计值整体处于较低水平，误差评价值处于控制水平以外，由此表明，利用梯形法求解值存在较大误差，即拟合结果不符合设计精度要求。利用辛普森法的拟合结果显示，x0.1%和x0.01%水平的误差平均值分别为17.45%和12.26%，处于误差合格控制范围内的组数分别为20组和16组，二者的最大误差值分别为-16.62%和-23.05%；利用辛普森算法的拟合结果值整体处于较低水平，且各拟合结果误差处于合格范围，相对于梯形算法利用辛普森法的拟合结果基本能够符合设计有关精度要求。利用龙贝格法的拟合结果显示，x0.1%和x0.01%水平的误差平均值分别为13.18%和11.61%，处于误差合格控制范围内的组数分别为18组和15组，二者的最大误差值分别为-16.64%和-17.49%；相对理论假象值设计值整体较低且误差范围处于合理控制范围内，由此表明利用龙贝格法求解P-Ⅲ型积分曲线整体可以符合相关设计精度。利用自适应误差积分法的拟合结果显示，x0.1%和x0.01%水平的误差平均值分别为-8.92%和-12.18%，处于误差合格控制范围内的组数分别为20组和20组，二者的最大误差值分别为-12.66%和-10.44%%；利用自适应误差积分法求解的P-Ⅲ型积分曲线低于理想假象值，由此表明利用该方法进行拟合的结果表现出良好的精度与可靠性并符合有关设计精度要求[11]。

5 实测数据

利用相关文献[12]和资料中对水文序列的计算方法和理论可对曲线参数进行拟合估计，计算结果见表2，表中洪水序列1和2分别代表一般洪水和特大洪水序列。

由表2参数估计结果可知，利用文中所述4中方法的参数估计结果大致相同，且根据曲线变化趋势可以发现，原始序列利用文中所述四类方法均能较好的进行拟合，由此表明，目标函数适线法利用横坐标离差平方和具有较好的适用性和可靠性。另一方面结合离差平方和统计结果可知，自适应误差积分法的离差平方和在一般洪水序列和特大洪水序列拟合中值较小，而采用其他三种方法的离差平方和相差较大并与理想数据统计值存在一定的差异；龙贝格和辛普森算法的拟合效果低于梯形算法。

由表2统计结果运算时间可知，相对于其他算法，自适应误差分析法运算速度明显较快，且积分步长与精度和积分区间存在密切的关系；在相同步长时梯形法仅仅是对积分区间求和而不表现出辛普森法的拟合特征，由于龙贝格法的运算量和计算时间随迭代过程的减少而加快。综上所述，利用自适应误差积分法的离差平方和最小且运算速度较快，其拟合结果表现出较好的精度和可靠性。自适应误差积分法在适线效果和离差平方和方面均优于其他数值积分法。