ABSTRACT

The south of Guizhou province suffers serious and frequent flood disasters. In order to provide an efficient and reliable theoretical basis for design of the flood control project in this area, the Partial L-moments are applied at the Bamao, Gaoche, Huishui, Guiyang, Libo, Caopingtou and Baben stations for flood frequency analysis. By estimating the parameters of Generalized Extreme Value (GEV) distribution and matching Partial L-Moments to annual maximum flow series of 7 hydrological stations, the design flood is calculated and the flood frequency curve is fitted. The cumulative squares error is used to evaluate the fitting ability of Partial L-moments and L-moments. The results show that as censored level F0 value increases, the relative deviation of the design value is smaller, except the hydrological stations of Gaoche and Caopingtou. Partial L-Moments can describe the data better in flood analysis and improve the estimation precision of design flood. Partial L-moments is a reasonable and effective method of flood frequency analysis in the south of Guizhou province.

Keywords:Flood Frequency Analysis, Generalized Extreme Value Distribution, Parameters Estimation, Partial L-Moments

部分线性矩在黔南部地区洪水频率分析中的应用

王俊珍

贵州省大坝安全监测中心，贵州 贵阳

收稿日期：2018年7月2日；录用日期：2018年7月13日；发布日期：2018年7月23日

摘 要

黔南部地区洪水灾害严重且频繁，研究部分线性矩在该区洪水频率分析中的应用，为该区的水利防洪工程设计建设提供可靠的理论依据。本文在介绍部分线性矩原理的基础上，选取黔南部八茂、高车、惠水、贵阳、荔波、草坪头及把本7个水文站的洪峰流量资料，选用广义极值分布，将部分线性矩应用于黔南部地区洪水序列参数估计中，计算设计洪水值，并拟合洪水频率曲线，利用累积相对偏差平方和评价所得洪水频率曲线对洪峰序列的拟合效果，并与普通线性矩法的拟合成果进行比较分析。利用部分线性矩和普通线性矩对黔南部地区7个水文站洪水进行的分析表明：除了高车和草坪头站外，随着删失值F₀的增大，设计值的相对偏差值越小，说明部分线性矩法对高尾部经验点据拟合效果更好，提高了大重现期设计值估算精度。因此，是一种合理有效的洪水频率分析参数估计方法，可用于黔南部地区的洪水频率计算。

关键词 :洪水频率分析，广义极值分布，参数估计，部分线性矩

Copyright © 2018 by author and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

洪水因暴雨形成，且陡坡流急，植被不足，每遇暴雨，即形成暴涨暴落的陡峻洪峰，易使两岸村镇、城市及人民生命财产毁于一旦 [1] 。黔南部地区洪水灾害严重且频繁，修建水库和堤防等水利工程是必要的防洪措施。而推求设计洪水是防洪规划和建设的依据 [2] [3] 。在水利工程中，往往需要推求大重现期的洪水设计值，而设计值可靠与否，会直接影响到工程的安危 [4] [5] 。目前，常用的参数估计方法有矩法、适线法、权函数法、概率权重矩法、线性矩法等，这些方法的估计精度及拟合效果优劣各异 [6] [7] 。矩法有较大的估计偏差，尤其是对Cs的估计；概率权重矩法具有良好的不偏性，但对于Cs的抽样误差较大，尤其是短序列；权函数法提高了Cs的估计精度，但是权函数属于单参数估计，不能解决均值，特别是Cv的估计精度问题；国内通常采用适线法，其中适线法可灵活调整参数，但主观任意性较大，结果因人而异；这些方法在一定程度上改善了估计精度，但对于大洪水段设计值拟合效果稍差 [8] 。因此，20世纪90年代，Bhattarai [9] [10] 建议使用部分线性矩(Partial L-Moments)来拟合分布，对于长序列洪水，该法以截取频率分布拟合大流量值，效果良好，使外延的大重现期设计洪水值精度得到提高，且借助于计算机，容易实现计算过程。

本文选用黔南部地区7个水文站洪峰流量序列为例，以部分线性矩法和线性矩进行洪水频率曲线拟合，进而评价所得洪水频率曲线对序列拟合效果及设计值误差，以期为该区水利建设提供计算依据。

2. 部分线性矩

部分线性矩是用于估计删失样本的参数估计方法。给定一个排序样本 $x_1\le x_2\le \cdots \le x_n$ ，对于低删失，前4阶部分线性矩分别为

${{\lambda }^{\prime }}_1=\frac{1}{{}^{n}C{}_1}{\displaystyle \sum _{i=1}^n{x}_{\left(i\right)}^{*}}$ (1)

${{\lambda }^{\prime }}_2=\frac{1}{2}\frac{1}{{}^{n}C{}_2}{\displaystyle \sum _{i=1}^n\left({}^{i-1}C{}_1-{}^{n-i}C{}_1\right)x_{\left(i\right)}^{*}}$ (2)

${{\lambda }^{\prime }}_3\text{}\text{}\text{}\text{}\text{}\text{}=\frac{1}{3}\frac{1}{{}^{n}C{}_3}{\displaystyle \sum _{i=1}^n\left({}^{i-1}C{}_2-2{}^{i-1}C{}_1{}^{n-i}C{}_1+{}^{n-i}C{}_2\right)x_{\left(i\right)}^{*}}$ (3)

${{\lambda }^{\prime }}_4=\frac{1}{4}\frac{1}{{}^{n}C{}_4}{\displaystyle \sum _{i=1}^n\left({}^{i-1}C{}_3-3{}^{i-1}C{}_2{}^{n-i}C{}_1+3{}^{i-1}C{}_1{}^{n-i}C{}_2-{}^{n-i}C{}_3\right)x_{\left(i\right)}^{*}}$ (4)

其中，

$x_{\left(i\right)}^{*}=\{\begin{array}{l}0x_{\left(i\right)}\le x_0\\ x_{\left(i\right)}{x}_{\left(i\right)}>x_0\end{array}$ (5)

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

式中： $x_0$ 为低删失门限值。

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}$ (7)

其逆函数形式为

$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}$ (8)

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

Wang (1990) [11] 提出了GEV分布部分概率权重矩，当 $k\ne 0$ 时，其表达式为

${{\beta }^{\prime }}_r=\left(\xi +\frac{\alpha }{k}\right)\frac{1}{r+1}\left(1-F_0^{r+1}\right)-\frac{\alpha }{k}\frac{\Gamma \left(1+k\right)}{{\left(1+r\right)}^{1+k}}P\left(1+k,-\left(1+r\right)\log F_0\right)$ (9)

式中： $F_0=F\left(x_0\right)$ ， $x_0$ 为低删失门限值， $P(\cdot )$ 为不完全gamma函数，即

$P\left(1+k,-\left(1+r\right)\log F_0\right)={\displaystyle {\int }_0^{-\left(1+r\right)\log F_0}\frac{x^k{\text{e}}^{-x}}{\Gamma \left(1+k\right)}\text{d}x}$ (10)

把 $r=0,\text{1},\text{2}$ 代入式(9)，可得

${{\beta }^{\prime }}_0=\left(\xi +\frac{\alpha }{k}\right)\left(1-F_0\right)-\frac{\alpha }{k}\Gamma \left(1+k\right)P\left(1+k,-\log F_0\right)$ (11)

$\frac{2{{\beta }^{\prime }}_1}{1-F_0^2}-\frac{{{\beta }^{\prime }}_0}{1-F_0}=-\frac{\alpha }{k}\left[\frac{P\left(1+k,-2\log F_0\right)}{2^k\left(1-F_0^2\right)}-\frac{P\left(1+k,-\log F_0\right)}{1-F_0}\right]$ (12)

$\frac{\frac{2{{\beta }^{\prime }}_1}{1-F_0^2}-\frac{{{\beta }^{\prime }}_0}{1-F_0}}{\frac{3{{\beta }^{\prime }}_2}{1-F_0^3}-\frac{{{\beta }^{\prime }}_0}{1-F_0}}=\frac{\frac{P\left(1+k,-2\log F_0\right)}{2^k\left(1-F_0^2\right)}-\frac{P\left(1+k,-\log F_0\right)}{1-F_0}}{\frac{P\left(1+k,-3\log F_0\right)}{3^k\left(1-F_0^3\right)}-\frac{P\left(1+k,-\log F_0\right)}{1-F_0}}$ (13)

令上式左边为

$z=\frac{\frac{2{{\beta }^{\prime }}_1}{1-F_0^2}-\frac{{{\beta }^{\prime }}_0}{1-F_0}}{\frac{3{{\beta }^{\prime }}_2}{1-F_0^3}-\frac{{{\beta }^{\prime }}_0}{1-F_0}}$ (14)

则有

$z=\frac{\frac{P\left(1+k,-2\log F_0\right)}{2^k\left(1-F_0^2\right)}-\frac{P\left(1+k,-\log F_0\right)}{1-F_0}}{\frac{P\left(1+k,-3\log F_0\right)}{3^k\left(1-F_0^3\right)}-\frac{P\left(1+k,-\log F_0\right)}{1-F_0}}$ (15)

当 $F_0\ne 0$ 时，给定k的取值范围 $-0.5\le k\le 0.5$ ，分别令 $F_0=0.1~0.5$ ，根据式(15)分别计算各组取值对应的z值，按式(16)拟合曲线，求得曲线拟合系数如表1所示。

$k=a_0+a_{1}z+a_2{z}^2+a_3{z}^3+a_4{z}^4$ (16)

Hosking [12] (1990)给出线性矩与概率权重矩的前4阶关系为

${\lambda }_1={\beta }_0$ (17)

${\lambda }_2=2{\beta }_1-{\beta }_0$ (18)

${\lambda }_3=6{\beta }_2-6{\beta }_1+{\beta }_0$ (19)

${\lambda }_4=20{\beta }_3-30{\beta }_2+12{\beta }_1-{\beta }_0$ (20)

式中： ${\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3$ 为概率权重矩， ${\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4$ 为前4阶线性矩，而部分概率权重矩与部分线性矩有类似的线性关系，式(14)中 ${{\beta }^{\prime }}_0,{{\beta }^{\prime }}_1,{{\beta }^{\prime }}_2$ 由 ${{\lambda }^{\prime }}_1,{{\lambda }^{\prime }}_2,{{\lambda }^{\prime }}_3$ 代替，并根据式(17)~式(20)的关系可得

$z=\frac{\frac{{{\lambda }^{\prime }}_2+{{\lambda }^{\prime }}_1}{1-F_0^2}-\frac{{{\lambda }^{\prime }}_1}{1-F_0}}{\frac{\frac{1}{2}\left({{\lambda }^{\prime }}_3+3{{\lambda }^{\prime }}_2+2{{\lambda }^{\prime }}_1\right)}{1-F_0^3}-\frac{{{\lambda }^{\prime }}_1}{1-F_0}}$ (21)

根据式(1)～式(3)，由样本计算 ${{\widehat{\lambda }}^{\prime }}_1,{{\widehat{\lambda }}^{\prime }}_2,{{\widehat{\lambda }}^{\prime }}_3$ ，并代入式(21)，可得z的估计量 $\widehat{z}$ 为

$\widehat{z}=\frac{\frac{{{\widehat{\lambda }}^{\prime }}_2+{{\widehat{\lambda }}^{\prime }}_1}{1-F_0^2}-\frac{{{\widehat{\lambda }}^{\prime }}_1}{1-F_0}}{\frac{\frac{1}{2}\left({{\widehat{\lambda }}^{\prime }}_3+3{{\widehat{\lambda }}^{\prime }}_2+2{{\widehat{\lambda }}^{\prime }}_1\right)}{1-F_0^3}-\frac{{{\widehat{\lambda }}^{\prime }}_1}{1-F_0}}$ (22)

表1. 不同 $F_0$ 下式(16)的拟合系数

上式计算出 $\widehat{z}$ 值，由式(16)计算参数k的估计值 $\widehat{k}$ ，其系数由表1查算。式(11)和式(12)可得参数 $\alpha ,\xi$ 的估计值 $\widehat{\alpha },\widehat{\beta }$ 分别为

$\widehat{\alpha }=-\frac{\widehat{k}}{\Gamma \left(1+\widehat{k}\right)}\frac{\frac{{{\widehat{\lambda }}^{\prime }}_2+{{\widehat{\lambda }}^{\prime }}_1}{1-F_0^2}-\frac{{{\widehat{\lambda }}^{\prime }}_1}{1-F_0}}{\frac{P\left(1+k,-2\log F_0\right)}{2^k\left(1-F_0^2\right)}-\frac{P\left(1+k,-\log F_0\right)}{1-F_0}}$ (23)

$\widehat{\beta }=\frac{{{\widehat{\lambda }}^{\prime }}_1}{1-F_0}+\frac{\widehat{\alpha }}{\widehat{k}}\left[\Gamma \left(1+\widehat{k}\right)\frac{P\left(1+k,-\log F_0\right)}{1-F_0}-1\right]$ (24)

当 $F_0=0.0$ 时，部分线性矩转化为普通线性矩。Hosking [12] (1990)推出了GEV分布下普通线性矩为

${\lambda }_1=\xi +\alpha \left[1-\Gamma \left(1+k\right)\right]/k$ (25)

${\lambda }_2=\alpha \left(1-2^{-k}\right)\Gamma \left(1+k\right)/k$ (26)

${\tau }_3=2\left(1-3^{-k}\right)/\left(1-2^{-k}\right)-3$ (27)

${\tau }_4=\left[5\left(1-4^{-k}\right)-10\left(1-3^{-k}\right)+6\left(1-2^{-k}\right)\right]/\left(1-2^{-k}\right)$ (28)

当 $-0.5<{\tau }_3<0.5$ 时，三个参数的估计量计算公式分别为

$\widehat{k}=7.8590C+2.9554C^2$ (29)

式中： $C=2/3\left(3+{\tau }_3\right)-\ln 2/\ln 3$ 。

$\widehat{\alpha }={\lambda }_2\widehat{k}/\left[\left(1-2^{-\widehat{k}}\right)\Gamma \left(1+\widehat{k}\right)\right]$ (30)

$\widehat{\xi }={\lambda }_1-\widehat{\alpha }\left[1-\Gamma \left(1+\widehat{k}\right)\right]/\widehat{k}$ (31)

式中： ${\lambda }_1,{\lambda }_2$ 可由样本进行计算。

4. 实例应用

本文选取黔南部7个水文测站的年最大流量资料，经还原处理，资料满足一致性要求，研究GEV分布部分线性矩法应用于洪水序列拟合的效果，资料长度见表2。

4.1. 绘制频率曲线

计算不同F0下各站年最大洪峰流量序列GEV分布参数估计值，由式(8)计算洪水设计值，并绘制理论频率曲线拟合图，如图1所示。

由图1看出，7个测站普通线性矩法(F₀ = 0.0)对低尾部经验点据拟合结果较好，随着F₀值的增大，洪水频率曲线对低尾部经验点据拟合结果变差，高尾部经验点据拟合结果得到显著改善，这表明部分线性矩法对高尾部洪水值拟合效果好。

4.2. 拟合效果分析

采用累积相对偏差平方和δ对上述拟合结果进行定量分析。式(32)为P = 50%~98%时，对应实测值与设计值累积偏差平方和的计算公式。

$\delta ={\displaystyle \sum _{i=i|P=50\%}^{i|P=98\%}{\left(\frac{x_i-{\widehat{x}}_i}{x_i}\right)}^2}$ (32)

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

由表2可得出，与普通线性矩(F₀ = 0.0)相比较，部分线性矩的设计值偏小。此外，除了高车和草坪头站，

表2. 不同F₀下各站资料长度和设计值误差比较

图1. 黔南部地区测站洪峰流量频率曲线图

其它5站随着F₀的增大，高尾部设计值计算偏差越小，即频率曲线与高尾部(P > 50%)经验点据越来越接近，提高了设计值估算精度，这与图1拟合曲线结果相一致。而高车和草坪头站，在F₀ ≤ 0.4时，高尾部经验点据拟合效果与其它6站变化趋势相同，但F₀ = 0.5时，设计值估算偏差变大，说明并非F₀越大越好，即并非截掉小洪水值越多越好，而是选择合适的F₀值，可提高设计值估算精度。因此，通过增大F₀值的方式来拟合高尾部洪水值是可行的，提高了该区设计值估算精度。

5. 结论

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

文章引用

王俊珍. 部分线性矩在黔南部地区洪水频率分析中的应用
The Application of Partial L-Moments for Flood Frequency Analysis in the South of Guizhou[J]. 水资源研究, 2018, 07(04): 418-425. https://doi.org/10.12677/JWRR.2018.74047

参考文献