高阶线性矩法在陕北地区洪水频率分析中的应用

Using Higher-Order L-Moments for Flood Frequency Analysis in Northern Shaanxi

摘要: 在介绍高阶线性矩原理的基础上，选取陕北交口河、张家山、赵石窖、绥德、刘家河、张村驿、林家村及神木8个水文站的洪峰流量资料进行广义极值分布高阶线性矩的参数估计，评价拟合效果和设计值的计算偏差，并与普通线性矩法拟合结果进行比较分析。利用高阶线性矩和普通线性矩对陕北8个水文站洪水频率进行的分析表明：随着阶数η的增大，设计值的相对偏差值越小，说明高阶线性矩法对洪峰序列的大洪水值拟合效果更好，提高了设计值估算精度。因此，高阶线性矩法是一种合理有效的洪水频率分析参数估计方法，可为大重现期设计洪水值的计算提供依据。

Abstract: In order to provide an efficient and reliable theoretical basis for design floods in the northern Shaanxi province, the higher-order L-Moments are applied in flood frequency analysis based on the principles of higher-order L-Moments. The annual maximum flood series of 8 hydrological stations at Jiaokou, Zhangjiashan, Zhaoshiyao, Sudie, Liujia, Zhangcunyi, Linjiacun and Shenmu Rivers are selected for case study. The parameters of Generalized Extreme Value (GEV) distribution and the design floods are estimated. The flood frequency curves are fitted, and the cumulative of squares error is regard as an indicator to evaluate the effect, and compared with the traditional Method of moments. The results show that higher-order L-Moments possess good statistical performance, which can describe the data series much better than lower-order L-Moments in flood analysis. Consequently, this method is reasonable and feasible, and would be provided the basis for the flood quantile calculation.

1. 引言

洪水是我国最频繁和严重的自然灾害之一。洪水因暴雨形成，具有突发性，破坏性和危险性，可以使村镇、城市及人民生命财产毁于一旦 [1] 。陕北地区洪水灾害严重且频繁，修建水库和堤防等水利工程是必要的防洪措施。而推求设计洪水是防洪规划和建设的依据 [2] [3] 。洪水频率分析法是推求设计洪水的主要途径，在分布函数给定的情况下，参数估计是设计洪水计算的重要工作 [4] 。国内外学者在此方面进行了大量研究，并相继提出了概率权重矩法、权函数法、线性矩法、极大似然法等参数估计方法。实际中，理论频率曲线与实测大洪水点据在大重现期段拟合难以获得满意的结果。为了解决这一问题，国外学者Wang [5] [6] ，Pandey [7] [8] ，Deng [9] [10] 等先后进行了大洪水段的计算研究，主要包括部分概率权重矩(Partial Probability Moments)法，部分线性矩(Partial Linear Moments)法，最小交互熵(Minimum Cross-Entropy)法等。这些方法认为小洪水值在设计洪水估算中具有滋扰行为，可应用截取概率分布进行相对大洪水值的推求。研究结果显示，在大洪水段取得较好的拟合效果，但计算复杂。国内工程设计采用适线法，通过调整参数或给定权重进行经验适线，这些方法适线结果存在一定的主观性和经验性 [11] [12] [13] [14] 。20世纪90年代，高阶线性矩法 [15] (higher-order L-Moments)开始引入洪水分布参数计算，通过提高线性矩阶数来改善大洪水段拟合效果，使外延的大重现期设计洪水值精度得到提高，且借助于计算机，容易实现计算过程。

本文选用陕北地区8个水文站洪峰流量序列为例，以高阶线性矩法进行洪水频率曲线拟合，进而评价所得洪水频率曲线对序列拟合效果及设计值误差。

2. 高阶线性矩

高阶线性矩是高阶概率权重矩 [16] 的线性组合。给定样本容量为n的样本， $x_{in}$ 表示n个样本点中第i个由小到大的排序值，变量服从 $F\left(x\right)=P\left(X\le x\right)$ 分布。Hosking(1990)提出 $x_{in}$ 的数学期望为

$E\left[X_{i,n}\right]=\frac{n!}{\left(i-1\right)!\left(n-i\right)!}{\displaystyle {\int }_0^{1}x\left(F\right)F^{i-1}{\left(1-F\right)}^{n-i}\text{d}F}$ (1)

高阶线性矩定义为

${\lambda }_1^{\eta }=E\left[X_{\left(\eta +1\right)\left(\eta +1\right)}\right]$ (2)

(3)

${\lambda }_3^{\eta }=\frac{1}{3}E\left[X_{\left(\eta +3\right)\left(\eta +3\right)}-2X_{\left(\eta +2\right)\left(\eta +3\right)}+X_{\left(\eta +1\right)\left(\eta +3\right)}\right]$ (4)

${\lambda }_4^{\eta }=\frac{1}{4}E\left[X_{\left(\eta +4\right)\left(\eta +4\right)}-3X_{\left(\eta +3\right)\left(\eta +4\right)}+3X_{\left(\eta +2\right)\left(\eta +4\right)}-X_{\left(\eta +1\right)\left(\eta +4\right)}\right]$ (5)

式中： ${\lambda }_1^{\eta }$ 表示样本容量为 $\eta +1$ 中最大变量的数学期望， ${\lambda }_2^{\eta }$ 表示样本容量为 $\eta +2$ 中最大变量的数学期望，表示样本容量为 $\eta +3$ 中最大变量的数学期望， ${\lambda }_4^{\eta }$ 表示样本容量为 $$ 中最大变量的数学期望。

当 $\eta =0$ 时，高阶线性矩转化为普通线性矩(Hosking, 1990)。随着 $\eta$ 增高，高阶线性矩对随机变量的较大值更为依赖。高阶线性矩的变差系数 ${\tau }_2^{\eta }$ ，偏态系数 ${\tau }_3^{\eta }$ 和峰态系数 ${\tau }_4^{\eta }$ 分别为

${\tau }_2^{\eta }=\frac{{\lambda }_2^{\eta }}{{\lambda }_1^{\eta }}$ (6)

(7)

${\tau }_4^{\eta }=\frac{{\lambda }_4^{\eta }}{{\lambda }_2^{\eta }}$ (8)

给定一个排序样本 $x_1\le x_2\le \cdots \le x_n$ ，高阶线性矩的估计量分别为

${\hat{\lambda }}_1^{\eta }=\frac{1}{{}^{n}C{}_{\eta +1}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{}^{i-1}C{}_{\eta }{x}_{\left(i\right)}}$ (9)

${\hat{\lambda }}_2^{\eta }=\frac{1}{2}\frac{1}{{}^{n}C{}_{\eta +2}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left({}^{i-1}C{}_{\eta +1}-{}^{i-1}C{}_{\eta }{}^{n-i}C{}_1\right)x_{\left(i\right)}}$ (10)

${\hat{\lambda }}_3^{\eta }=\frac{1}{3}\frac{1}{{}^{n}C{}_{\eta +3}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left({}^{i-1}C{}_{\eta +2}-2{}^{i-1}C{}_{\eta +1}{}^{n-i}C{}_1+{}^{i-1}C{}_{\eta }{}^{n-i}C{}_2\right)x_{\left(i\right)}}$ (11)

${\hat{\lambda }}_4^{\eta }=\frac{1}{4}\frac{1}{{}^{n}C{}_{\eta +4}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left({}^{i-1}C{}_{\eta +3}-3{}^{i-1}C{}_{\eta +2}{}^{n-i}C{}_1+3{}^{i-1}C{}_{\eta +1}{}^{n-i}C{}_2-{}^{i-1}C{}_{\eta }{}^{n-i}C{}_3\right)x_{\left(i\right)}}$ (12)

${}^{n}C{}_i=\frac{n!}{i!\left(n-i\right)!}$ (13)

3. 广义极值分布及其高阶线性矩

广义极值分布(Generalized extreme value distribution, GEV)的分布函数为

$F\left(x\right)=\{\begin{array}{l}\exp \left\{-{\left[1-\frac{k}{\alpha }\left(x-\xi \right)\right]}^{\frac{1}{k}}\right\};k\ne 0\\ \exp \left\{-\exp \left[-\frac{1}{\alpha }\left(x-\xi \right)\right]\right\};k=0\end{array}$ (14)

其逆函数形式为

$x\left(F\right)=\{\begin{array}{l}\xi +\frac{\alpha }{k}\left[1-{\left(-\ln F\right)}^k\right];k\ne 0\\ \xi -\alpha \ln \left(-\ln F\right);k=0\end{array}$ (15)

式中：k为形状参数， $\alpha$ 为尺度参数， $\xi$ 为位置参数。

当 $k\ne 0$ 时，即为Hosking (1985) [17] 给出的GEV分布概率权重矩(PWM)

${\beta }_r={\displaystyle {\int }_0^{1}x\left(F\right)F^r\text{d}F}$ (16)

把式(14)和式(15)代入式(16)，可推出

(17)

当时，式(14)为EVI型分布(Gumbel分布)。Greenwood (1979) [18] 给出了GEV分布的PWM

$\left(r+1\right){\beta }_r=\xi +\alpha \left[\epsilon +\ln \left(r+1\right)\right]$ (18)

式中， $\epsilon =0.5772156649\cdots$ 为欧拉常数。

将式(17)代入式(1)~式(5)，可推出 $k\ne 0$ 时，前4阶线性矩分别为

${\lambda }_1^{\eta }=\xi +\frac{\alpha }{k}\left[1-\Gamma \left(1+k\right){\left(\eta +1\right)}^{-k}\right]$ (19)

${\lambda }_2^{\eta }=\frac{\left(\eta +2\right)\alpha \Gamma \left(1+k\right)}{2!k}\left[-{\left(\eta +2\right)}^{-k}+{\left(\eta +1\right)}^{-k}\right]$ (20)

${\lambda }_3^{\eta }=\frac{\left(\eta +3\right)\alpha \Gamma \left(1+k\right)}{3!k}\left[-\left(\eta +4\right){\left(\eta +3\right)}^{-k}+2\left(\eta +3\right){\left(\eta +2\right)}^{-k}-\left(\eta +2\right){\left(\eta +1\right)}^{-k}\right]$ (21)

$\begin{array}{c}{\lambda }_4^{\eta }=\frac{\left(\eta +4\right)\alpha \Gamma \left(1+k\right)}{4!k}[-\left(\eta +6\right)\left(\eta +5\right){\left(\eta +4\right)}^{-k}+3\left(\eta +5\right)\left(\eta +4\right){\left(\eta +3\right)}^{-k}\\ -3\left(\eta +4\right)\left(\eta +3\right){\left(\eta +2\right)}^{-k}+\left(\eta +3\right)\left(\eta +2\right){\left(\eta +1\right)}^{-k}]\end{array}$ (22)

同样，将式(18)代入式(1)~式(5)，可推出k = 0时，前4阶线性矩分别为

${\lambda }_1^{\eta }=\xi +\alpha \left[\epsilon +\ln \left(\eta +1\right)\right]$ (23)

${\lambda }_2^{\eta }=\frac{\left(\eta +2\right)\alpha }{2!}\left[\ln \left(\eta +2\right)-\ln \left(\eta +1\right)\right]$ (24)

${\lambda }_3^{\eta }=\frac{\left(\eta +3\right)\alpha }{3!}\left[\left(\eta +4\right)\ln \left(\eta +3\right)-2\left(\eta +3\right)\ln \left(\eta +2\right)+\left(\eta +2\right)\ln \left(\eta +1\right)\right]$ (25)

$\begin{array}{c}{\lambda }_4^{\eta }=\frac{\left(\eta +4\right)\alpha }{4!}[\left(\eta +6\right)\left(\eta +5\right)\ln \left(\eta +4\right)-3\left(\eta +5\right)\left(\eta +4\right)\ln \left(\eta +3\right)\\ +3\left(\eta +4\right)\left(\eta +3\right)\ln \left(\eta +2\right)-\left(\eta +3\right)\left(\eta +2\right)\ln \left(\eta +1\right)]\end{array}$ (26)

给定一个样本，由式(21)，式(20)，式(24)，式(25)和式(7)可推出

$k=a_0+a_1\left[{\tau }_3^{\eta }\right]+a_2{\left[{\tau }_3^{\eta }\right]}^2+a_3{\left[{\tau }_3^{\eta }\right]}^3{a}_4{\left[{\tau }_3^{\eta }\right]}^4$ (27)

取k为 $-0.5\le k\le 0.5$ ，分别令 $\eta =0,1,2,3,4，5$ ，计算各阶取值对应的 ${\tau }_3^{\eta }$ 值，按照(27)拟合曲线，求得式(27)的系数，如表1所示。

应用式(10)和式(11)计算 ${\tau }_3^{\eta }$ ，根据表1系数，由式(27)计算参数k的估计值 $\hat{k}$ 。利用式(20)计算参数 $\alpha$ 的估计值 $\hat{\alpha }$ ( $\hat{k}\ne 0$ )为

$\hat{\alpha }=\frac{{\lambda }_2^{\eta }\times \hat{k}\times 2!}{\Gamma \left(1+\hat{k}\right)\left(\eta +2\right)\left(-{\left(\eta +2\right)}^{-\hat{k}}+{\left(1+\eta \right)}^{-\hat{k}}\right)}$ (28)

由式(19)和式(28)可推出参数 $\xi$ 的估计值 $$ ( $\hat{k}\ne 0$ )为

$\hat{\xi }={\hat{\lambda }}_1^{\eta }-\frac{\hat{\alpha }}{\hat{k}}\left(1-\Gamma \left(1+\hat{k}\right){\left(\eta +1\right)}^{-\hat{k}}\right)$ (29)

4. 实例应用

采用陕西省交口河、张家山、赵石窖、绥德、刘家河等8个水文测站的年最大洪峰资料，经分析，资料满足一致性要求，研究GEV型分布的高阶线性矩法进行洪水序列大洪水段的拟合，资料的基本概况如表2所示。

4.1. 绘制频率曲线

根据表2结果，由式(15)计算洪水设计值，利用期望公式(30)估计经验频率，其表达式为

$P=\frac{m}{n+1}\times 100$ (30)

式中，P为频率，m为系列样本中的第m项，n为系列样本中的总数。

根据以上计算成果，绘制各站不同阶线性矩的GEV分布理论频率曲线拟合图，详见图1。

根据拟合的频率曲线图，对各站不同阶线性矩法进行比较，分析理论频率曲线与经验点据的拟合情况。由图1可以看出，在P = 0~50%的频率区间，8个测站均随着线性矩阶数的增大，理论频率曲线和经验点据的拟合

Table 1. Coefficients , a 1 , a 2 , a 3 and a 4 for Eq. (27)

表1. 式(27)的拟合系数 $a_0$ ， $a_1$ ， $a_2$ ， $a_3$ 和 $a_4$

Table 2. Comparison of quantile errors using different orders L-Moments

表2. 不同阶线性矩法的设计值偏差比较

(a) 交口河(b) 张家山 (c) 赵石窖(d) 绥德 (e) 刘家河(f) 张村译 (f) 林家村(g) 神木

Figure 1. Fitting of the GEV distribution to annual maximum flows for different orders L-Moments

图1. 年最大洪峰序列各阶线性矩频率曲线拟合图

效果不显著；在P = 50%~99%的频率区间，8个测站均随着线性矩阶数的增大，理论频率曲线和经验点据的拟合效果得到显著改善，说明本方法计算大重现设计洪水值的偏差小。因此，从实测序列值拟合结果可知，采用高阶线性矩法进行拟合计算时，大洪峰段取得较好的拟合效果。

4.2. 拟合效果分析

应用累积相对偏差平方和 $\delta$ 分析上述不同阶线性矩法对经验点据的拟合效果。式(31)为P = 50%~98%时，对应实测值与设计值累积偏差平方和的计算公式。

$\delta ={\displaystyle \underset{i=i|P=50\%}{\overset{i|P=98\%}{\sum }}{\left(\frac{x_i-{\hat{x}}_i}{x_i}\right)}^2}$ (31)

式中： $x_i$ 为实测值； ${\hat{x}}_i$ 为设计值。计算结果如表2所示。

由表2看出，除了张村译站以外，其余7个站点的 $\delta$ 均随着线性矩阶数 $\eta$ 的升高而减小，即大洪水段实测值与设计值之间的偏差逐渐减小，说明大重现期设计值在阶数较高时与实测值更为接近，这与图1中的拟合结果相一致。而张村驿站，在 $\eta \le 4$ 时，大洪水段实测值拟合效果与其它7站变化趋势相同，但 $\eta =5$ 时，设计值估算偏差变大，说明并非 $$ 越大越好，即并非给与的权重越大越好，而是选择合适的 $$ 值，可提高设计值估算精度。所以，通过提高线性矩阶数的方式来拟合序列大洪水段是可行的，可为外延的大重现期设计值提供有利依据。

综合以上分析，在研究区大重现期设计值计算中，建议选用不超过4阶线性矩。

5. 结论

以陕北地区8个水文测站的年最大洪峰流量序列为例，以高阶线性矩法对洪峰序列进行参数估计，评价其拟合效果及设计值的偏差。结果表明，高阶线性矩法通过增大值的方式来拟合大洪水值是可行的，且随着 $\eta$ 的增大，大重现期设计值计算偏差越小，可为研究区防洪工程建设提供理论依据。但并非 $\eta$ 越大越好，应选择合适的 $\eta$ 值，可提高设计值估算精度。但高阶线性矩法还存在一定的问题，由于P-III型分布的逆函数无法表达为显式函数，推导P-III型分布参数与高阶线性矩法的关系式存在较大的困难。因此，在今后的研究中，应重点研究高阶线性矩法在P-III型分布参数估计中的应用研究。