近60年来潦河径流变化特征分析

Analysis of the Runoff Characteristics Variation of Liao River in Recent 60 Years

摘要: 在气候变化和人类活动的双重作用下，潦河径流过程发生了较为明显的改变。以近60年潦河万家埠水文站的逐月径流资料为基础，用不均匀系数、Mann-Kendall趋势检验、滑动t突变检验、功率谱周期分析等方法对径流量的年内分配、年际过程、年际突变和周期性等特征进行了分析。结果表明：万家埠站径流量的年内分配不均匀且呈单峰型；年径流量呈现显著上升趋势，并于1968年和1999年发生突变；年径流量存在21a的显著周期。

Abstract: Under the dual action of climate change and human activities, the runoff process of Liao River has significantly changed. Based on the monthly runoff data of Wanjiabu hydrological station in recent 60 years, the characteristics of annual distribution, interannual process, interannual abrupt change and periodicity of runoff were analyzed by using the methods of non-uniformity coefficient, Mann-Kendall trend test, sliding t abrupt change test and power spectrum period analysis. The results show the annual runoff distribution of Wanjiabu station not uniform with a unimodal pattern; the annual runoff has a significant upward trend and significant period of 21a, a sudden change is occurred in 1968 and 1999.

1. 引言

河川径流是重要的地表水资源，其形成和演化，对维系区域生态平衡和促进经济社会发展具有重要作用 [1]。近年来，鄱阳湖流域的赣江、修水、饶河和修水等河流径流序列均呈不显著的增加趋势 [2]，且降水是修河中上游径流变化的主要控制因素 [3]；随着修河流域年降水量、年平均气温的上升和柘林水库的蓄水运行，修河干流年内分配变得更为均匀且流域内河川径流的年内分配过程被大大改变了 [4]；其中潦河径流在气候变化以及人类活动的影响下也发生了变化，在1956~2008年间，万家埠站的径流量呈小幅度上升趋势 [5]，而年最大流量则呈现出了不显著的减少趋势 [6]，且年径流量于1968、1972和1976等年份发生了突变 [7]。故，潦河径流变化特征及其原因还需要进一步探索。

2. 研究区概况与数据

潦河发源于九岭山脉，为修水最大支流。潦河流域地处114˚87'E~115˚38'E，28˚38'N~29˚04'N (图1)，流域

Figure 1. Schematic diagram of Liao River system

图1. 潦河水系示意图

面积约为3548 km²，流域大致呈东西长，南北短的四边形；流域内地势西高东低，自西向东依次为中低山区、丘陵和冲积平原。流域位于中亚热带湿润季风气候区，气候温湿，四季分明，多年平均降水量1778 mm，年径流量42.8 × 10⁸ m³，4月~7月降水量约占全年总降水量的60% [4]。

本文用潦河下游万家埠水文站1956~2020年的逐月径流数据(2012年数据缺失)对潦河径流变化特征进行分析，使用数据均来自全国水雨情信息网。

3. 主要研究方法

3.1. Mann-Kendall检验法

Mann-Kendall检验法(简称M-K法)是目前较常用的水文变量趋势性检验方法，该方法优点是变量不需要服从一定的分布，适合于水文变量的趋势检验。Z为一个服从正态分布的统计量，当Z > 0表明原序列可能有上升趋势，Z < 0表示原序列可能有下降趋势；当|Z| ≥ 1.96和2.56时，分别表示通过了α = 0.05和α = 0.01的显著性水平检验 [8]。

3.2. 滑动t检验法

滑动t检验法是通过考察极端降水指数时间序列两组子序列的平均值差异是否超过了一定的显著性水平，来判断样本序列是否发生突变 [9]。此方法需要反复改变子序列长度进行试验比较，以提高计算结果的可靠性。给定显著性水平α，查t分布表可得到临界值t_α/2。若|t_i| < t_α/2，则认为基准点前后的两个序列均值无显著差异，否则认为在基准点时刻有突变发生。

3.3. 功率谱周期分析方法

功率谱周期方法是分析时间序列常用的方法，根据功率谱图的特点可以直观地揭示出离散数据序列的周期性并可以快速、准确地识别出径流序列的周期性，操作较为简单 [10]。其分析方法步骤如下：

给定离散时间序列 $Y\left(t\right)\left(t=\text{1},\text{2},\cdots ,n\right)$，并将径流量序列进行标准化，设序列Y(t)的均值和标准差分别为 $\bar{Y}$ 和S，则标准化后的新序列为：

$y^{\prime }=\left[Y\left(t\right)-\bar{Y}\right]/S\left(t=1,2,\cdots ,n\right)$ (1)

计算新序列的自相关函数：

$r\left(\tau \right)=\frac{1}{n-\tau }{\displaystyle {\sum }_{t=1}^{n-\tau }{y}^{\prime }\left(t\right)y^{\prime }\left(t+\tau \right)}\left(\tau =1,2,\cdots ,m\right)$ (2)

其中 $\tau$ 为自相关函数的时延，m为最大时延，一般取为 $\frac{n}{3}~\frac{n}{10}$，在本文中取 $m=\frac{n}{3}$。

通过以下公式计算不同时延的粗谱估计值：

${\hat{S}}_k=\{\begin{array}{l}\frac{1}{2m}\left[r\left(0\right)+r\left(m\right)\right]+\frac{1}{m}{\displaystyle {\sum }_{\tau =1}^{m-1}r\left(\tau \right)}\left(k=0\right)\\ \frac{1}{m}\left[r\left(0\right)+2{\displaystyle {\sum }_{\tau =1}^{m-1}r\left(\tau \right)\cos \frac{k\pi }{m}\tau }+r\left(m\right)\cos k\pi \right]\left(k=1,2,3,\cdots ,m-1\right)\\ \frac{1}{2m}\left[r\left(0\right)+{\left(-1\right)}^{m}r\left(m\right)\right]+\frac{1}{m}{\displaystyle {\sum }_{\tau =1}^{m-1}{\left(-1\right)}^{\tau }r\left(\tau \right)}\left(k=m\right)\end{array}$ (3)

式中：k为谐波波数。

粗谱估计值与真实谱存在一定的误差，故采用Hanning平滑系数进一步对 ${\hat{S}}_k$ 做平滑处理得到其精细谱值：

${\tilde{S}}_k=\{\begin{array}{l}0.5{\hat{S}}_0+0.5{\hat{S}}_1\left(k=0\right)\\ 0.25{\hat{S}}_{k-1}+0.5{\hat{S}}_k+0.25{\hat{S}}_{k+1}\left(1\le k<m\right)\\ 0.5{\hat{S}}_{k-1}+0.5{\hat{S}}_k\left(k=m\right)\end{array}$ (4)

以k为横坐标、精细谱 ${\tilde{S}}_k$ 为纵坐标绘制功率谱图，并根据图形的极大值确定序列得到周期T_k：

$T_k=\frac{2m}{k}$ (5)

同时比较标准谱与精细谱，进行周期的显著性检验。标准谱检验有以下两种情况：①如果式(2)中自相关函数r(1)为较大正值，则序列为“红噪声”，采用式(6)进行计算。

$S_{ok}=\bar{S}\left[\frac{1-r^2\left(1\right)}{1+r^2\left(1\right)+2r\left(1\right)\cos \frac{k\pi }{m}}\right]$ (6)

式中： $\bar{S}$ 为第m + 1个精细谱估计值的均值，即

$\bar{S}=\left[\frac{1}{2m}\left[\tilde{S}\left(0\right)+\tilde{S}\left(m\right)\right]+\frac{1}{m}{\displaystyle {\sum }_{k=1}^{m-1}\tilde{S}\left(k\right)}\right]$ (7)

②若r(1)接近零或为负值，序列则为“白噪声”，则采用式(8)进行计算。

${S^{\prime }}_{ok}=S_{ok}\left(\frac{{\chi }_{\alpha }^2}{\upsilon }\right)$ (8)

其中 $S_{ok}=\bar{S}$，自由度 $\upsilon$ 可通过式(9)计算。

$\upsilon =\frac{2n-\frac{m}{2}}{m}$ (9)

在给定显著性水平条件下，查表得到 ${\chi }_{\alpha }^2$ 值，计算 ${S^{\prime }}_{ok}$ ；若 ${\tilde{S}}_k>{S^{\prime }}_{ok}$，则表明波数k对应的周期是显著的。

4. 结果与讨论

4.1. 年内分配特征

4.1.1. 年内变化过程

万家埠站的径流量年内分配呈单峰型曲线(图2)，径流量集中在4~7月，该时期径流量为21.10 × 10⁸ m³，占年径流量的60%，其中6月径流量最大为6.57 × 10⁸ m³，12月径流量最小为1.04 × 10⁸ m³。近60年来径流量的年内分配出现了峰值后移现象，在1968年以前峰值位于5月，之后则位于6月，这可能是受江西省在1958~1966年间开展的水土保持工作的影响 [11]。

4.1.2. 年内不均匀系数变化

径流量年内分配不均匀系数的多年平均值为0.81 (图3)，这表明径流量年内分配极不均匀，其中最小值为0.42 (2001年)，最大值为1.22 (1958年)；将万家埠站的年径流数据按年代分为6个时段(表1)，其中1960 s不均匀系数值最大为0.88，1980 s最小为0.66，且自1960 s到1980 s不均匀系数值持续下降，随后又呈现出波动上升的趋势，但总体来说径流量年内分配不均匀系数呈现出微弱的下降趋势，其趋势线斜率为−0.029/10a，这说明径流量的年内分配不均匀性有所降低，这可能是受人类活动的影响。截至到2007年，潦河流域内已建了中型水库7座、小型水库229座，总库容达到了3.04 × 10⁸ m³，这对于潦河径流的年内分配具有一定的影响。

Table 1. Statistical characteristic values of interdecadal variation of runoff at Wanjiabu station

表1. 万家埠站径流量年代际变化统计特征值

Figure 2. Annual distribution of runoff of Liao River

图2.潦河径流量年内分配

Figure 3. Interannual variation of the Cv

图3. 不均匀系数年际变化

4.2. 年际变化分析

在1956~2020年，万家埠站的年径流量呈现出上升趋势(图4)，其变化率大致为1.36 × 10⁸ m³/(10a)。进一步利用Mann-Kendall检验法对径流量的变化趋势进行分析得出其Z值为2.14 (表2)，表明万家埠站的年径流量呈现出显著上升的趋势。在汛期(4~7月)径流量呈现出了不显著的上升趋势，其中径流量最大月(6月)呈现出微弱的上升趋势，7月则呈现出了显著的上升趋势；在非汛期(8月~次年3月)径流量则呈现非常显著的上升趋势，其中1、9和12月的上升趋势都非常显著。

Table 2. The results of Mannal-Kendall trend test

表2. Mannal-Kendall趋势检验结果

注：“*”和“**”分别表示序列变化的趋势显著和非常显著。

径流量的五年滑动平均值围绕年径流量趋势线呈现出“增–减”的波动(图4)。其中在1956~1973年、1978~1994年和2007~2013年间年径流量大致呈现出波动上升趋势，而在1973~1977年、1995~2006年和2014年后则大体呈现出下降趋势。径流量大致呈现出“丰–枯”交替变化的特征(表1，图5)，在1960 s、1980 s和2000 s径流量低于多年平均径流量35.03 × 10⁸ m³，表现为枯水期，其中2000 s的平均径流量最小值为28.21 × 10⁸ m³，低于多年平均径流量6.82 × 10⁸ m³。而在1970 s、1990 s和2010 s平均径流量则高于多年平均值，表现为丰水期，其中1990 s的平均径流量最大值为45.46 × 10⁸ m³，高出多年平均值10.43 × 10⁸ m³，此时鄱阳湖流域正处于偏湿时期 [12]。

Figure 4. Interannual variation process of runoff of Liao River

图4. 潦河径流年际变化过程

Figure 5. Annual cumulative anomaly of runoff

图5. 年径流量累计距平

4.3. 突变分析

利用滑动t检验，选取子序列 $n_{\text{1}}=n_{\text{2}}=\text{7}$，临界值为t_0.05 = ±2.18对万家埠站的径流可能的突变年份进行检验。1968和1999年的滑动t值分别为−2.25和2.81，达到了显著水平，而1977和1987年分别为1.35和−1.89，均没有达到显著水平(图6)。因此确定径流量的突变年份为1968和1999年，这与刘剑宇等的研究结果较为一致 [7]。1969~1999年间万家埠站的径流量发生了显著增加(图7)，相较于1956~1968年径流量增加了11.22 × 10⁸ m³，而在2000~2020年间平均径流量则下降至33.94 × 10⁸ m³，相较于1956~1968年径流量仍增加了6.32 × 10⁸ m³，这可能与鄱阳湖流域降雨量在1969年之前呈下降趋势，之后呈上升趋势有关 [13]。

4.4. 周期性分析

利用功率谱对万家埠站的径流量进行周期性分析，发现当年径流量在波数k为2、6和11时，精细谱取得极大值，表明径流与波数存在着对应的周期(图8)。根据式(5)计算得对应的周期分别为21a、7a和3.8a。由于万家埠站的年径流量自相关函数r(1)的值为0.22，故利用式(8)进行周期显著性检验。计算得出S_ok = 0.06和υ = 5.5，查表得 ${\chi }_{0.05}^2=11.07$，根据式(8)计算得出在α = 0.05水平下的 ${S^{\prime }}_{ok}=0.\text{11}$，将其与精细谱进行对比发现，当k = 2时即周期为21a时通过显著性检验，而7a和3.8a均没有通过显著性检验，故万家埠站的年径流量存在一个21a的显著周期和两个不显著周期。这与前人研究万家埠站最大径流量呈现22a的周期大致吻合 [6]。由于3.8a周期和修水流域实测降水量的3.6a周期相接近，故推测万家埠站的径流量变化可能是受到了修水流域降水量变化的影响 [14]，而21a的周期变化由于与鄱阳湖流域年降水变化周期22a相接近，故推测其可能是受到鄱阳湖流域降水量变化的影响 [15]。

Figure 6. Sliding t-test of annual runoff

图6. 年径流量的滑动t检验

Figure 7. Comparison of runoff before and after abrupt change

图7. 突变前后径流量对比

Figure 8. Annual runoff power spectrum

图8. 年径流量功率谱图

5. 结论

潦河径流量的年内分配为单峰型并出现了峰值后移的现象，在1968年以前峰值位于5月，之后则位于6月，且径流量年内分配的不均匀性呈现出微弱下降的年际变化趋势，这可能是由于人类活动的影响；近60年潦河年径流量呈现出显著的上升趋势，并大致呈现出“丰–枯”交替的变化特征，其中在1960 s、1980 s和2000 s表现为枯水期，在1970 s、1990 s和2010 s则表现为丰水期；潦河年径流量在1968年和1999年发生突变现象；潦河年径流量存在21a、7a和3.8a的变化周期，其中21a为显著周期，并推测该周期可能是受到鄱阳湖流域降水变化的影响。

基金项目

江西省教育厅科技项目(GJJ180903)。