江河水沙变化突变性与周期性分析方法及比较

刘 茜,王延贵

(1.北京交通大学土木建筑工程学院,北京 100044;2.国际泥沙研究培训中心,北京 100048)

江河水沙变化突变性与周期性分析方法及比较

刘 茜1,王延贵2

(1.北京交通大学土木建筑工程学院,北京 100044;2.国际泥沙研究培训中心,北京 100048)

为了更有效地分析江河水沙序列的突变性和周期性变化,在总结现有研究方法的基础上,以长江大通站和宜昌站的水沙变化为例,就江河水沙变化的突变性和周期性分析方法进行了分析与对比。不同方法的原理、计算结果及适用性的对比分析结果表明:有序聚类法计算简便,结果精确,最适用于突变性分析;小波分析法发展成熟,可用于一般的水沙序列周期性的分析计算。

江河水沙;突变性;周期性;长江;大通站;宜昌站

径流量和输沙量是江河水沙最重要的两个特征值,随着流域气候变化、人类活动频繁加剧等因素的影响,江河径流量和输沙量将不断变化。在长时段系列年中,江河水沙的长期变化特性主要表现为趋势性、突变性和周期性三方面,目前江河水沙变化研究也主要从上述三方面开展。鉴于江河水沙态势变化的重要性,相关的研究与研究方法较多,其中文献[1]已对趋势性分析方法做了详细的总结与比较。

常用于江河水沙变化突变的分析方法主要包括Mann-Kendall突变分析法(MK法)、有序聚类分析法、均值差异T检验法、Pettitt突变点检验法、Fisher最优分割法、累计相关曲线法等;周期性的研究方法主要包括Morlet小波分析法、最大熵谱法、周期图法、Lemper-Ziv复杂度方法。在突变性研究方面,丁文峰等[2]用均值差异T检验法对嘉陵江流域的径流量和输沙量序列进行了突变分析;王金花等[3]用MK法以内蒙古皇甫川为例进行了突变点计算,认为MK法计算水沙突变点具有确定性;陈远中等[4]对有序聚类法进行了改进,使其更适用于提取时间序列的转折点(突变点);王国庆等[5]运用斯波曼秩次相关检验法检验了黄河中游无定河流域水文序列的趋势性,并利用有序聚类分析法推估了人类活动对序列的显著影响干扰点;丁明军等[6]利用Pettitt突变点检验法,分析了1961—2007年鄱阳湖周边地区气温序列的变化特征。在周期性分析方面,刘锋等[7-9]利用小波分析法分析了黄河利津站、渭河流域、洞庭湖区域等水沙变化的周期性特点;邵骏等[10]应用最大熵谱估计法分析了岷江上游紫坪铺水文站1937—2004年实测径流量的周期变化规律。

综上所述,针对江河水沙态势变化的突变性和周期性,都曾采用多种方法开展研究。不同方法的原理与计算过程存在一定的差异,相应都有各自的特点和局限性,但相关的综合研究和对比分析仍然较少。本文结合《中国河流泥沙公报》[11]提供的资料,基于长江大通站和宜昌站水沙变化特点,就各种方法进行分析和比较。

1 江河水沙变化突变性分析方法

水沙序列各阶段之间的显著变化即为突变,分析突变年份对于查找水沙突变的原因有重要作用。以下采用几种常用的方法计算突变点,并利用洪水年份、水库建设及水土保持等方面的实际资料对计算结果进行佐证,判断突变点的合理性。

1.1 MK法

MK法主要是通过统计量子、方差σr和标准化变量M来实现突变计算的,计算公式为

式中:P为径流系列所有对偶观测值(Ri,Rj,i＜j)中Ri＜Rj出现的次数;N为系列长度。根据式(1),将水沙序列按时间顺序计算的M值组成曲线C1,按时间反序列计算所得M值乘以-1构成曲线C2,两条曲线交点即为突变点。若交点位于给定置信度(95%)水平线之间,则突变点在统计意义上是显著的。

图1和图2为MK法计算绘制的大通站和宜昌站水沙过程突变分析结果,可以看出,大通站径流量突变点出现在1955年、1980年、1987年、2003年附近。1954年、2002年长江干流发生洪水,MK法计算得到的突变点1955年和2003年与实际突变年份相差1年,1980和1987年未表现出突变。输沙量突变点为1994年,主要是1989年实施的长江上游水土保持重点防治工程(以下简称“长治”工程)和1988—1996年间长江流域修建的16座大型水库,导致大通站输沙量在20世纪90年代初发生了较显著的变化。

图1 大通站水沙MK法突变分析结果

图2 宜昌站水沙MK法突变分析结果

长江上游分别于1954年、1969年、1998年、2005年发生洪水,MK法计算的宜昌站径流突变点出现在1958年、1960年、1971年、2001年,除1960年外,虽略有偏差(3或4年),但都与实际情况接近。输沙量的突变点出现在2001年左右,这与2002年、2003年分别修建了紫坪铺水库和三峡水库吻合。

1.2 有序聚类分析法

采用有序聚类分析法推估水文序列的可能显著干扰点子0,实质上就是推求最优分割点,使同类之间的离差平方和最小,而类与类之间的离差平方和相对较大。最优点分割计算公式为

利用有序聚类分析法计算大通站和宜昌站的水沙过程总离差平方和,结果如图3和图4所示。可以看出,大通站径流量在1956年和2003年出现突变点,与实际突变点1954、2002年相差较小。输沙量在1991年处为一级突变点,二级突变点为1968年、1984年、2002年,其中一级突变点1991年与1989年实施的“长治”工程有一定关系;二级突变点1968年和2002年分别与1967年、2003年修建的丹江口水库、三峡水库相关,二级突变点1984年为无效的突变点。

图3 大通站有序聚类分析法突变分析结果

图4 宜昌站有序聚类分析法突变分析结果

从图4可以看出,宜昌站径流量变化的突变年份主要是在1956年、1968年、2000年和2005年,均与实际洪水年份较近。输沙量变化在2001年附近为一级突变点,二级突变点有1968年、1985年、1991年,其中2001年与紫坪铺水库和三峡水库的修建相关,1968年与当年洪水引起的大水大沙相关,1991年与1989年实施的“长治”工程以及1993—1996年间修建的6座大型水库有关,1985年为无效的计算点。

1.3 均值差异T检验法

定义样本长度为n的序列的突变指数AI和统计量t分别为

式中:X1、S1分别为基准年前的平均值和标准差; X2、S2分别为基准年后的平均值和标准差;m1、m2分别为基准年前后两段序列的样本长度;S2P为联合样本方差。统计量t服从自由度为m1+m2-2的T分布,当给出一定的显著水平α=0.05,如t＞tα,则在α的显著水平上,基准年两侧的均值有明显的差异,即在基准年处出现了突变。

采用T检验法对大通站和宜昌站1950—2010年径流量序列和输沙量序列作突变分析,结果如图5所示。大通站径流在1955年发生了一次突变,而输沙量在1956年后为持续明显减少的状态,不能根据此方法分析出明确的突变点,若取统计量的极值点为突变点,则1992年为突变点。宜昌站径流在1968年、1993年、1999年、2005年发生了突变;输沙量在1964年后统计量t均大于临界值,为持续明显减小的状态,若取统计量的最大值为突变点,则2001年为突变点,如还考虑其他极大值点,则有1968年、1993年也为突变点。

图5 均值差异T检验法突变分析结果

1.4 Pettitt突变点检验法

Pettitt突变点检验法基于非参数方法检测水文时间序列的突变点,计算较简便,可以明确突变的时间。对于具有n个样本量的时间序列(x1,x2,…, x1),构成一秩序列

若t0时刻满足kt0=max Sk(k=2,3,…,n),则t0点处为突变点,其显著性水平可以通过下式计算:

方娟娟等[12]运用Pettitt法对大通站1946—2009年径流进行分析,得出大通站年均流量无显著突变,洪季流量突变年份为1958年,枯季流量突变年份为1987年、1955年(二级突变点)。

1.5 Fisher最优分割法

Fisher最优分割法的分类依据是样本的总离差平方和最小,分割的原则是使得各类内部样本之间差异最小,而各类之间的差异最大。即对任意指定的分类数K,总将n个样本分为k类,且使各类直径(类内部样本之间的差异程度)的总离差平方和达到最小,即

D(ia-1+1,ia)

式中:直径D为样本离差平方和;P(n,k)为分类的目标函数;i为类的分割点。

万新宁等[13]应用该方法对1950—2000年宜昌和大通站的水沙量进行了分析,结果表明宜昌站径流量突变发生在1968年和1975年,其中1975年无突变特征,为无效的突变点;输沙量的突变点为1961年和1968年,但1961年没有发生突变,为无效的突变点。大通站径流量的突变点为1955年和1986年,其中1986年突变不成立;输沙量突变点为1966年和1980年,其中1980年无突变特征。

1.6 突变性分析方法综合比较

表1为各方法计算结果汇总。与实际情况相比,MK法计算的突变点数量不一致且偏差较大;有序聚类法在输沙量计算时有一个无效点,但突变年份偏差不超过2年,计算结果精确;均值差异T检验法计算径流量突变点数与实际情况较一致(只相差1个突变点),但无法计算出输沙量的有效突变点; Pettitt突变点检验法计算的大通站径流量突变点缺少2002年;Fisher分割法计算结果与实际突变年份相差较大,大通站计算结果有无效的突变点,宜昌站的计算结果也不准确。

从方法原理上看,MK法计算的突变点明显,但当交叉点在置信区间外时,突变点在统计意义上是不显著的;另当两条曲线在置信区间内多次相交时,判断的突变点过多。有序聚类分析法在水沙序列没有显著变化的条件下,能求得的突变点明显且较精确;当序列有显著的趋势性变化时,运用二级突变点的概念进行再分割,得出的二级突变点也可作为参考依据。均值差异T检验法通过统计量来判断突变点,计算结果较为精确,但计算序列持续明显的变化时,得不到有效的突变点。Pettitt法可以明确地计算出一个突变点,但计算的突变点数与实际情况相比往往要偏少。Fisher最优分割法能够根据定义的目标函数确定分几期较优,依据是总离差平方和最小,但其计算过程较为复杂,且精度较低。

综上所述,各突变性分析方法都有不足之处,但综合分析,有序聚类分析法计算简便,且能得到较为精确的结果,最适用于对江河水沙序列进行突变性分析计算。

2 江河水沙变化周期性分析方法

水沙量由于各种自然或人为因素的影响扰动而随时间上下波动,这个波动在某一时间尺度内有一定的重复性,即为周期性。周期性分析有利于水沙变化趋势的预测。目前用于周期性分析的方法有Morlet小波分析法、最大熵谱法、周期图法、Lemper-Ziv复杂度方法等。

表1 江河水沙突变性分析方法突变年份计算结果

2.1 Morlet小波分析法

小波分析的基本思想是用一簇小波函数系来表示或逼近某一信号或函数。对于给定的小波函数Ψ(t),水文时间序列f(t)∈L2(t)(定义在实轴上、可测的平方可积函数空间)的连续的小波变换为

式中:a为尺度因子,反映小波的周期长度;b为时间因子,反映时间上的平移;Wf(a,b)称为小波变换系数,能同时反映频域参数a和时域参数b的特性; ¯Ψ(t)为小波函数Ψ(t)的复共轭函数。以b为横坐标,a为纵坐标即可绘制关于Wf(a,b)的二维等值线图,即小波变换系数时频分布图。将时间域上关于a的所有小波系数的平方值在b域上积分,就可得到小波方差,即

小波方差随尺度a的变化过程,称为小波方差图,它能反映信号波动的能量随尺度a的分布。因此,小波方差图可用来确定信号中不同尺度扰动的相对强度和存在的主要时间尺度,即主周期,峰值越大,说明其周期震荡越强。

由图6可以看出大通站径流量和输沙量在不同时段上各时间尺度的强弱分布,其中径流量存在4个较为明显的峰值,依次对应着22 a、15 a、11 a和7 a的时间尺度;输沙量存在3个较为明显的周期,依次为22 a、12 a和5 a。

图6 大通站小波方差突变分析结果

2.2 最大熵谱法

熵谱是以熵的概念为基础进行的谱估计,其主要思想是,在观察时间之内的估计值等于观察值;在观察时间之外的取值不做任何假定,即保持最随机、最不确定性,也就是使得熵为最大,从而得到一种新的非线性谱估计法,即最大熵谱法(maximum entropymethod)[10]。在已知m+1个自相关函数值的条件下,最大功率谱S(T)满足下列关系:

式中:σ2m为方差;Δt为离散序列的时间间隔,在等间隔序列中一般取Δt=1;i为虚数单位;T为周期,bmj为m阶自回归系数,m称为自回归的截止阶。具体计算时自回归系数bmj采用预报误差过滤系数的递推计算法求得。

王盼成等[14-15]采用该方法计算了1923—2000年大通站的水沙序列,流量在频率为0.070和0.150时出现极大值,对应的周期为16.0 a和6.6 a,可认为年均流量具有16 a、7 a的周期性变化;输沙量极值点出现的频率为0.0635、0.1368、0.3174,对应的周期为15.8 a、7.4 a、3.2 a。

2.3 周期图法

周期图法[16]是将观测到的有限个样本序列(x1,x2,…,xn)利用FFT算法做傅立叶变换,按下式进行功率谱估计:

式中:PT(k)为(x1,x2,…,xn)的周期;XN(k)为(x1, x2,…,xn)的傅立叶变换;N为系列长度。在PT(k)出现峰值处,可能存在周期。

王盼成等[14-15]用周期图法计算了大通站1923—2000年水沙通量的周期。径流量变化周期为15.6 a和7 a;输沙量存在的周期为8 a、16 a和3.2 a。

2.4 Lemper-Ziv复杂度方法

Lemper-Ziv复杂度方法[17]是一种动态非线性时间序列分析方法,用于度量给定的有限长度符号序列的复杂性。现在普遍认为复杂度的物理意义在于复杂度反映了一个时间序列随着序列长度的增加出现新模式的速率。运用Lemper-Ziv复杂度法计算水沙变化的复杂度,如复杂度越大,说明水沙在窗口长度时期内随时间出现的新变化越多,发生新变化的速率越快,表明这一时期的水沙变化是无序而复杂的;反之,复杂度越小,则说明发生新变化的速率越慢,水沙变化是规则的,周期性越强。

文献[18]在Lemper-Ziv复杂度方法的基础上进行周期分析,认为嘉陵江北碚站年径流复杂度序列在20世纪50年代末至70年代末存在4 a左右的准周期,70年代以后的阶段则无显著周期。

2.5 周期性分析方法比较

表2为各周期性分析方法对大通站径流量和输沙量周期的计算结果。可以看出,大通站径流量的计算结果较为一致,都含有16 a和7 a的周期性,但小波分析法还包括了时间尺度较长的22 a,更为全面;而输沙量的周期性计算结果略有偏差,与计算的时间序列不一致有关。

表2 大通站径流量和输沙量周期的计算结果a

Morlet小波分析法发展成熟,具有明确的水文意义,通过小波变换系数,可直接得出任意长度的水沙序列的变化周期,且结果准确全面,目前应用最为广泛。最大熵谱法分辨率高,任何对周期的微小偏离都会使谱值迅速下降,且适应序列时间较短的情况;但与截止阶m的取值关系甚大,截止阶的确定存在一定的误差,可能导致结果产生较大偏差。周期图法不需要先估计自相关函数,适用于长时间序列的情况;但谱估计不满足一致性估计条件,方差不会随着样本长度的增大而趋于0,即序列长度一定时,要保证足够高的谱分辨率,谱估计的方差会很大,谱的正确性会很差。Lemper-Ziv复杂度方法物理意义明确,但计算过程还需利用功率谱估计复杂度的周期,过程烦琐,目前尚未得到推广应用。由上述各周期性分析方法的特点,水沙序列的周期分析时宜采用小波分析法,其计算过程清晰明了,且适用条件限制少,能得到较为精确的结果,适用于分析计算一般水文序列的周期性。

2.6 关于水沙周期性的讨论

目前,有不少学者对江河的水沙序列进行了周期性分析,但对于水沙的周期性仍有许多不同的看法。江河径流量的变化受降雨量的影响较大,降雨的周期变化会引起江河径流量的丰枯变化,因而径流量也会呈现周期性变化。天然条件下输沙量与径流量存在正向关系,径流量的周期导致输沙量的变化也存在一定的周期性,并与径流量周期基本对应;若流域人类活动频繁,如水库修建、水土保持、引水引沙等,会导致下垫面条件发生变化,引起水沙量大幅减少,将难以直接看出明显的周期性,不能直接对水沙序列进行周期性分析,需对输沙量序列进行趋势分量分离,分离后的序列(即线性回归后得到的残差序列)再进行周期性分析。但当水沙量具有明显增加或减少的趋势时,就没有利用水沙周期性预测变化趋势的必要性了。

3 结 论

a.MK法、有序聚类分析法、均值差异T检验法、Pettitt突变点检验法、Fisher最优分割法等突变性分析方法中,有序聚类法计算简便,且能得到较为精确的结果,最适用于对水沙序列进行突变性分析。

b.Morlet小波分析法、最大熵谱法、周期图法和Lemper-Ziv复杂度方法等周期性分析方法中,小波分析法计算简洁,适用性好,且结果精确,可用于一般的水沙序列周期性的分析计算。

[1]王延贵,刘茜,史红玲.江河水沙变化趋势分析方法与比较[J].中国水利水电科学研究院学报,2014,12 (1):41-47.(WANG Yangui,LIU Xi,SHI Hongling. Analyticalmethods and their comparisons for water and sediment variation trends[J].Journal of China Institute of Water Resources and Hydropower Research,2014,12 (1):41-47.(in Chinese))

[2]丁文峰,张平仓,任洪玉.近50年来嘉陵江流域径流泥沙演变规律及驱动因素定量分析[J].长江科学院院报,2008,25(3):23-27.(DING Wenfeng,ZHANG Pingcang,REN Hongyu.Quantitative analysis on evolution characteristics and driving factors of annual runoff and sediment transportation changes for Jialing River[J]. Journal of Yangtze River Scientific Research Institute, 2008,25(3):23-27.(in Chinese))

[3]王金花,张荣刚,张攀,等.两种水沙系列突变点算法的对比分析:以内蒙古皇甫川为例[J].中国水土保持, 2009(12):43-44.(WANG Jinhua,ZHANG Ronggang, ZHANG Pan,et al.Comparison on catastrophe point calculation of two water and sediment series[J].Soil and Water Conservation in China,2009(12):43-44.(in Chinese))

[4]陈远中,陆宝宏,张育德,等.改进的有序聚类分析法提取时间序列转折点[J].水文,2011,31(1):41-44. (CHEN Yuanzhong,LU Baohong,ZHANG Yude,et al. Improvement of sequential cluster analysis method for extracting turning point of time series[J].Journal of China Hydrology,2011,31(1):41-44.(in Chinese))

[5]王国庆,贾西安,陈江南,等.人类活动对水文序列的显著影响干扰点分析:以黄河中游无定河流域为例[J].西北水资源与水工程,2001,12(3):13-15.(WANG Guoqing,JIA Xi’an,CHEN Jiangnan,et al.Analysis onthe transition point of hydrological series impacted by human activities:a case study of Wudinghe basin in the middle reache of the Yellow River[J].Northwest Water Resources&Water Engineering,2001,12(3):13-15.(in Chinese))

[6]丁明军,郑林,杨续超.1961—2007年鄱阳湖周边地区气温变化趋势分析[J].中国农业气象,2010,31(4): 517-521.(DING Mingjun,Zhen Lin,YANG Xuchao. Trend analysis of long-term temperature time series in the area around Poyang Lake from 1961 to 2007[J].Chinese Journal of Agrometeorology,2010,31(4):517-521.(in Chinese))

[7]刘锋,陈沈良,彭俊,等.近60年黄河入海水沙多尺度变化及其对河口的影响[J].地理学报,2011,66(3): 313-323.(LIU Feng,CHEN Shenliang,PENG Jun,et al. Multi-scale variability of flow discharge and sediment load of Yellow River to sea and its impacts on the estuary during the past 60 years[J].Acta Geographica Sinica, 2011,66(3):313-323.(in Chinese))

[8]胡安焱,刘燕,郭生练,等.渭河流域水沙多年变化及趋势分析[J].人民黄河,2007,29(2):39-41.(HU Anyan,LIU Yan,GUO Shenglian,et al.Trend analysis of annual runoff and sediment transportation changes in Wei River Basin[J].Yellow River,2007,29(2):39-41.(in Chinese))

[9]李正最,谢悦波,徐冬梅,等.洞庭湖水沙变化分析及影响初探[J].水文,2011,31(1):45-53.(LI Zhengzui, XIE Yuebo,XU Dongmei,et al.Runoff-sediment variation and its effect on the Dongting Lake[J].Journal of China Hydrology,2011,31(1):45-53.(in Chinese))

[10]邵骏,袁鹏,李秀峰,等.基于最大熵谱估计的水文周期分析[J].中国农村水利水电,2008(1):30-33.(SHAO Jun,YUAN Peng,LI Xiufeng,et al.Application of maximum entropy method in the analysis of hydrologic period[J].China Rural Water and Hydropower,2008 (1):30-33.(in Chinese))

[11]中华人民共和国水利部.中国河流泥沙公报[G].北京:中国水利水电出版社,2010.

[12]方娟娟,李义天,孙昭华,等.长江大通站径流量变化特征分析[J].水电能源科学,2011,29(5):9-12.(FANG Juanjuan,LIYitian,SUN Zhaohua,et al.Analysis of runoff change characteristics at Datong Station of Yangtze River [J].Water Resources and Power,2011,29(5):9-12.(in Chinese))

[13]万新宁,王九发,何青,等.长江中下游水沙通量变化规律[J].泥沙研究,2003(4):29-35.(WAN Xinning, WANG Jiufa,HEQing,etal.Water and sediment fluxes in the middle and lower Yangtze River[J].Journal of Sediment Research,2003(4):29-35.(in Chinese))

[14]王盼成,贺松林.长江大通站水沙过程的基本特征:Ⅰ.径流过程分析[J].华东师范大学学报:自然科学版, 2004(2):72-80.(WANG Pancheng,HE Songlin.The basic character on process of runoff and sediment discharge at Datong Station of the Changjiang River:Ⅰ. analysis on runoff process[J].Journal of East China Normal University:Natural Science,2004(2):72-80.(in Chinese))

[15]贺松林,王盼成.长江大通站水沙过程的基本特征:Ⅱ.输沙过程分析[J].华东师范大学学报:自然科学版, 2004(2):81-86.(HE Songlin,WANG Pancheng.The basic character on process of runoff and sediment discharge at Datong Station of the Changjiang River:Ⅱ. analysis on sediment transport process[J].Journal of East China Normal University:Natural Science,2004(2):81-86.(in Chinese))

[16]余训锋,马大玮,魏琳.改进周期图法功率谱估计中的窗函数仿真分析[J].计算机仿真,2008,25(3):111-114.(YU Xunfeng,MA Dawei,WEI Lin.Simulation analysis of window function in power spectrum estimation based on modified periodogram[J].Computer Simulation,2008,25(3):111-114.(in Chinese))

[17]侯威,封国林,高新全,等.基于复杂度分析冰芯和石笋代用资料时间序列的研究[J].物理学报,2005,54 (5):2441-2447.(HOU Wei,FENG Guolin,GAO Xinquan,et al.Investigation on the time series of ice core and stalagmite based on the analysis of complexity[J]. Acta Physice Scientia,2005,54(5):2441-2447.(in Chinese))

[18]全球江河水沙变化与河流演变响应[R].北京:国际泥沙研究培训中心,2010.

Comparison of analyticalmethods of runoff and sediment load mutation and periodical variation

LIU Xi1,WANG Yangui2(1.School ofCivil Engineering,Beijing Jiaotong University,Beijing 100044,China;2.International Research and Training Center ofErosion and Sediment,Beijing 100048,China)

To analyze the mutation and periodical variation of runoff and sediment load more effectively,the existing methods are summarized and then a detailed comparison is carried on.The experiment data used comes from the Datong Station and Yichang Station of Yangtze River.By comparing the principle,result and applicability of each method,it is concluded that the sequential clustermethod ismost suitable formutation analysiswith both simple calculation and accurate result,while the wavelet analysismethod is highly developed and recommended for periodical variation analysis.

river runoff and sediment;mutation;periodical variation;Yangtze River;Datong Station;Yichang Station

TV14

A

1006-7647(2015)02-0017-07

10.3880/j.issn.1006 7647.2015.02.004

2013-12-16 编辑:熊水斌)

中国水利水电科学研究院研究专项(沙集1334)

刘茜(1990—),女,湖南邵阳人,硕士研究生,主要从事防灾减灾及防护工程研究。E-mail:liuxi08@foxmail.com