受到自然环境剧烈变化、水文过程循环以及人类活动的影响，流域的产水产沙会呈现阶段性的趋势变化，同时具备周期性的特点^[1-3]。宁夏别称“塞上江南”，在黄河与清水河交汇处形成一个得天独厚的灌溉区，即卫宁灌域。对小流域内部的水沙变化成因进行量化的深入研究，既有助于加深对小流域水沙关系认知，又有利于卫宁灌域水沙调控机制的发展与完善，为灌域的灌溉提供更有利的条件。近60年来，黄河上游的年径流量和年输沙量呈现波动下降趋势^[4]，泥沙量的减少与宁夏段大型水库等水利设施的运行有关^[5]，而流域年径流量与气候变化密切相关，可通过纵向对比降雨、泥沙以及年内分配进行分析^[6]，SWAT分布式水文模型也被用于宁夏地区流域的径流和泥沙预测^[7-8]。但目前大多数宁夏境内的水沙研究集中在中南部流域，关于卫宁灌域的水沙研究较少，且缺少多方法结合的时间尺度分析研究。

以往的区域水文分析工作中，原始值、距平值、距平百分率和标准化距平等指标应用较多。理论上，空间差异小的范围使用绝对指标和相对指标的分析结果差异不明显。因此笔者对水沙的累计距平值分析，径流结果可显示持续增加、持续减少或比较稳定，并相应判断丰水期、枯水期和平水期^[9]，泥沙结果的变化也可利用累计距平值分析说明。为了进一步揭示研究区水沙变化特征，使用Mann-Kendall趋势分析及突变检验(M-K检验)和小波方差分析探究水沙突变情况及周期性。M-K检验是非参数检验，其不要求样本分布，也能忽略少数异常值或缺失值的干扰^[10]；而小波变换则可以对序列周期性进行解释，并分析序列变化特征^{[1 1-13]}。国内外不少学者已经利用M-K检验对降雨、径流等进行特征研究。在南非的Limpopo流域^[14]以及印度的Delhi和Mumbai^[15]，均使用M-K检验分析降雨序列和径流序列的趋势以及突变点；豫南地区的气候变化也利用M-K检验和小波分析的方法^[16]；黄河上游的水沙趋势^[4]、长江大通站的径流变化特征^[17]亦是使用M-K检验和小波分析进行研究。

笔者选用位于宁夏中卫市的泉眼山水文站1989—2019年的水沙时间序列资料，利用累计距平法、Mann-Kendall趋势分析及突变检验和小波方差分析，从多时间尺度系统分析卫宁灌域水沙变化趋势及周期特征，以期得出较为准确的变化规律，并阐明特殊时间点水沙变化的原因，以期为当地不同研究尺度水文预测的需求提供数据支持。

卫宁灌域是清水河流域的子流域之一，位于宁夏回族自治区中卫市内，涉及沙坡头区及中宁县，处于黄河主河道与清水河干流交汇处，为耕地的灌溉提供优良丰厚的水资源。流域总面积为1 353.74 km²，泉眼山水文站位于黄河与清水河交汇处(E 105.545°、N 37.471°)，多年平均降雨量为233.55 mm(图 1)。

图 1(Fig. 1)

图 1 研究区及泉眼山水文站地理位置图 Fig. 1 Location of the study area and Quanyanshan hydrologic station

假设时间序列x有n个样本，构建秩序列^[16]：

式中：s_k为x时间序列中第i个数值大于第j个数值的数量累计值；UF₁=0；s_k为秩序列s_k的平均值；Var(s_k)为秩序列s_k的方差。把序列x倒置形成一个新序列x_n，x_n-1，…，x₁，对其进行与UF相同的运算则得到UB^[18]。

笔者采用连续小波方差分析流域水沙变化周期特征，所用小波函数为Morlet函数^[19-20]。对时间序列x进行一维连续小波变换：

式中：W_f(a, b)为f(t)在相平面(a, b)处的小波变换系数；a为尺度因子；b为时间因子；V_ar(a)为小波方差。在研究的时间序列中，可用V_ar(a)变化的峰值对应周期去反映起主要作用的周期^[21]。

基于小波变换所使用的的一元线性回归在未消除短期高频的噪音时，无论是月序列还是年序列的分析结果都未能达到较好的可信度^[22-23]，因此笔者先对年径流量与年输沙量序列进行标准化处理，以达到降噪的目的，随后进行数据延伸，减弱“边界效应”，最后再进行小波变换。

1989—2019年31 a间泉眼山水文站实测年径流量的年际变化情况与实测年均泥沙量的年际变化情况相似(图 2)。结合流域水沙距平值及累计距平值的变化趋势图(图 3)，年径流量的距平值曲线表面在1995和1996年的径流量比多年平均高，可以认为年径流量震荡上升，但年径流量的上升趋势并不明显，上升幅度非常有限(0.039亿m³/a)。1994、1995和1996年连续3 a的年径流量远高于平均值，令1996年往后的12 a间累计距平值均＞0，可认为1994、1995和1996这3 a的高年径流量是拉动研究段31 a的径流量趋势上升的主要原因，因丰水年(距平值>0的年份)的径流量占比优势不明显，故整体上升幅度不大。测站实际泥沙量变化趋势与径流量相似，在2012年前两者的变化距平值曲线是较为吻合的。与径流量不同，2011年后径流量逐渐上升，使总体的径流变化趋势为缓慢增加，而泥沙量在2011年后绝大部分时间呈现减少趋势，导致2011—2019年的泥沙量显著减少(P＜0.05)，最终令泥沙整体变化趋势为减少。

图 2(Fig. 2)

图 2 近31 a卫宁灌域水沙标准化结果对比 Fig. 2 Comparison of water and sediment normalized results in Weining irrigated watershed in recent 31 years

图 3(Fig. 3)

图 3 近31 a卫宁灌域水沙变化趋势 Fig. 3 Variation trend of water and sediment in Weining irrigated watershed in recent 31 years

如图 4所示，利用M-K突变检验法对泉眼山31 a的实测径流量进行突变分析，UF呈标准正态分布，显著性水平为0.05，观察UF曲线，在UF曲线和UB曲线相交前位于0.05置信水平临界线以内，未超出置信水平，UF曲线和UB曲线第1次相交后，不久又出现数次相交，分别发生在1990、1991、1992、1994和2001年，但结合距平值趋势分析，这些交点都并非突变点。在1989—1991年间UF曲线由正变负，说明在这期间年径流量都为增加趋势，而往后在1995年前却＜0，说明趋势慢慢趋向减弱。但无论增加还是减弱都没有超出置信水平，证明变化趋势不明显，这与距平值分析结果一致。而后从1996—2007年间。UF曲线均＞0，证明在此12 a间，径流量趋势为增加，且呈现较大幅度增长，但在2008—2016年间UF曲线却一直＜0，证明在这期间径流量又开始呈现下降趋势。2017—2019年UF曲线恢复＞0，证明径流量又开始呈增加趋势。UF曲线由始至终并未超过置信水平，证明研究31 a间径流量的变化并不明显，变化趋势较小，这与年距平值变化相一致。

图 4(Fig. 4)

图 4 近31 a卫宁灌域年径流量的Mann-Kendall曲线 Fig. 4 Mann-Kendall curve of annual runoff in Weining irrigated watershed in recent 31 years

如图 5所示，对研究区31 a的实测径流量进行M-K突变分析，UF为标准正态分布，显著性水平为0.05。UF曲线在1989—1993年间在0值上下波动，说明在这期间泥沙量的变化为先增后减，累计量近0，这与距平值和累计距平值曲线分布结果一致。1993年后、UF曲线与UB曲线相交即2004年前，UF曲线由负到正，且累计正值大于负值，证明这段时间泥沙量变化慢慢由减变增、由少变多。而在2005年后，UF曲线均位于0值以下，并在2011年超出0.05置信水平临界线，说明2005年后泥沙量呈现逐步下降的趋势，且在2011年后发生显著的下降变化，这与距平值分析结果相一致。UF曲线与UB曲线存在1个交点，发生在2005年前。结合距平值变化可知，2005年前距平值的线性回归趋势为正值，2005年后为负值，且累计距平值在2005年前波动增加，在2005年后不断下降，实测年泥沙量均值在2005年前后由0.039亿t/a变为0.011亿t/a，下降0.028亿t/a。对该序列进行突变检验，结果通过0.95的置信水平，证明近31 a研究区的泥沙量发生显著突变，从2005年后泥沙量开始显著减少^{[11, 24]}。

图 5(Fig. 5)

图 5 近31 a卫宁灌域年泥沙量的Mann-Kendall曲线 Fig. 5 Mann-Kendall curve of annual sediment in Weining irrigated watershed in recent 31 years

根据卫宁灌域泥沙量距平值趋势以及M-K突变检验结果可知，流域内在1989—2019年仅发生过1次泥沙量突变事件，突变时间点是2005年前。

卫宁灌域处于沙坡头区与中宁区交界处，流域内有黄河干道以及清水河出口，水资源丰富的同时，水土流失也较为严重。在20世纪未开展系统整治前，黄河经历长度约为19 km的沙坡头段后，水中泥沙含量增加约115万621.5 t/a^[25]，但在2000年，为响应国家西部大开发战略，黄河沙坡头水利枢纽工程计划(以下简称“枢纽工程”)在中卫市境内开展，2000年底临时工程开始施工，至2004年9月主体工程竣工，期间枢纽第一台机组发电时间为2004年3月26日^[26]。在水沙变化趋势分析中已经了解到，流域内年径流量在2001—2005年总体虽有上升，但是上升幅度较比2000年以前要低，而年泥沙量在与年径流量变化保持同步的基础上，于2005年前发生突变，这与枢纽工程的实际动工时间吻合。在2001年动工不久时，工程需要截流施工，水文站所监测到的年径流量与年泥沙量对应减少；而2002—2003年，在不需要截流施工的情况下，年径流量与年泥沙量在施工影响下有回复趋势但幅度较小；2004年9月工程竣工，年径流量对应大幅度减少，而年泥沙量也在这个时间段后发生突变。

由此可以认为造成流域内水沙变化的重要原因是枢纽工程。在水沙变化趋势分析中，径流量无明显突变的原因在于枢纽工程建设目的是为引水和水力发电，在径流本身年变化趋势就已经较弱的基础上强化年内调控^[27]，对年际变化影响不明显。但泥沙具有沉积的特点，在枢纽工程的沙坡头水库中大量泥沙沉积，在逐年减少水库库容的同时也减少流域内的流动泥沙^[26]，这使得枢纽工程竣工(2004年9月)后流域内部的泥沙发生明显的减少，并且逐年减少。

如图 6所示，在1989—2019年这31 a间，卫宁灌域的年径流量震荡更加频繁，而年泥沙量的震荡比径流量要弱。在这31 a的演化期间，径流量存在着6~15 a的主震荡周期，而泥沙量的震荡周期则集中在6 a。

图 6(Fig. 6)

图 6 卫宁灌域年径流量与泥沙量的小波变换系数实部时频分布与小波变换方差 Fig. 6 Time-frequency distribution and wavelet transform variance of real part of wavelet transform coefficient of annual runoff and sediment in Weining irrigated watershed

1) 年径流量的震荡周期在6 a的时间尺度下，年径流量变化较为稳定，极点出现比较集中，证明以6 a为尺度分析周期变化的话能够比较清晰地显示丰水年和枯水年的变化，具有全局性。根据小波变换系数实部时频分布图的结果(图 6a)显示，1995与2019年的波动周期最强，偏低径流波动周期信号年份为1996和2007年。卫宁灌域径流量在研究段的31 a间有着较为复杂的震荡周期变化，存在多重时间尺度的嵌套结构，根据小波变换方差图的结果显示，该区年径流量有着明显的6、9、15 a和准30 a的震荡周期信号。

2) 年泥沙量的震荡周期也是在6 a的时间尺度下，极点出现较为集中，而且极值变化较大，证明在研究时间段31 a内研究区的年泥沙量以6 a为周期变化明显，且不同周期泥沙量相差较大。在2005年后小波变换系数实部时频分布极点数量减少，2011年后变为0，说明在2005—2019年泥沙的周期变化逐渐减弱，并在2011年后不再出现周期变化。这与泥沙距平值、M-K检验的结果相一致。年泥沙量的周期震荡未出现嵌套结构，较为明显的仅有6 a 1个周期，1993年出现了最强震荡小波信号，而在1989和1996年出现2个偏低泥沙波动周期信号。其中，泥沙1996年的偏低信号是与径流量的偏低信号同步发生的。

根据水沙周期变化时间尺度图，以年径流量表现出的周期为参考值，设定固定时间尺度T(6、9、15 a)值，绘制水沙不同时间尺度下的小波系数实部变化过程曲线(图 7)，以表示水沙量的波动特性。

图 7(Fig. 7)

图 7 各时间尺度下卫宁灌域年径流量与年泥沙量小波系数实部变化过程 Fig. 7 Real parts of wavelet coefficients of annual runoff and annual sediment in Weining irrigated watershed at different time scales

在小波系数表征下，卫宁灌域年径流量与泥沙量序列各时间尺度小波系数实部变化曲线的过零点对应突变点，可见年径流量与泥沙量的小波系数实部变化过程高度相似。而根据小波系数实部时频分布分析结果，结合实际年径流量泥沙量的变化，在本节年径流量分析周期为6、9和15 a，年泥沙量分析周期为6 a。

1) 对年径流量首先分析较大的15 a时间尺度，在1989—1991年(枯)、1992—1997年(丰)、1998—2002年(枯)、2003—2007年(丰)、2008—2013年(枯)、2013—2019年(丰)共出现6次震荡，后续曲线向下且斜率绝对值变小，证明未来年径流量有由丰转枯，但是枯水期不会太严重或者太久；在9 a的时间尺度上，2003年前的震荡幅度较为强烈，而后震荡幅度逐渐减小，呈现缓慢弱化趋势；在6 a的时间尺度上，1995年前震荡幅度较弱，1996—2003年间震荡幅度较为强烈，2003年后震荡幅度再次减弱且趋于平缓，呈现弱化趋势。

根据M-K突变检验、小波系数实部变化过程和累计距平法分析出卫宁灌域在1989—2019年没有发生径流量突变事件。M-K检验结果显示，在2003年后UB统计量突破0.05置信水平临界线，证明从2003年开始研究区年径流量有了显著的变化，但UB与UF曲线并未有相交现象，说明2003年的变化同比整个研究时间段还未达到突变水平。而小波分析的结果与累计距平值的结果也证实2003年后的变化趋于稳定但同比整个31 a间未有突变水平。

2) 对研究区的31 a年泥沙量进行小波分析后，决定以唯一周期6 a进行研究。从年泥沙量的6 a的小波系数实部变化过程曲线中可以发现，1995年前震荡幅度较弱，1996—2005年间震荡幅度较为强烈，2005年后震荡幅度再次减弱且趋于平缓，2011年后震荡变幅呈现弱化趋势，这与M-K突变检验的分析结果相一致。且曲线末端为下降趋势，证明未来年泥沙量的变化仍以减少为主。

根据M-K突变检验、小波系数实部变化过程和累计距平法，分析出卫宁灌域在1989—2019年发生1次年泥沙量突变事件，突变时间在2005年前。在流域河床未发生显著变化的前提下^[28]，卫宁灌域在2005年泥沙量发生突变的因素在于人为干预的程度加大，在2000—2004年9月间，黄河沙坡头水利枢纽工程从开工至竣工运行，为泥沙量的减少奠定基础^[27]；以清水河流域为对象开展治沙减沙的工程^[29]，修建大量的淤地坝，拦沙量逐年上升，并在2005年达到极值^[30]，卫宁灌域也有着同样的泥沙变化，并在2005年前集中表现出突变情况。在2005年后，研究区的年泥沙量并未出现明显突变，仅表现为缓慢减少。

1) 1989—2012年研究区的年径流量与年泥沙量都发生相似的震荡，不同在于幅度大小，且泥沙变化发生稍微滞后于径流变化的发生时间点。但在2012—2019年的变化过程中，水沙趋势产生分化，年径流量微量增加，而年泥沙量显著降低，这与研究时间段的水沙总体变化相符。说明2012年后水沙的趋势变化决定1989—2019年31 a间水沙的总体趋势。

2) 在时间尺度上，径流量表现出6次“丰-枯”的变化过程，总体上呈现上升趋势，但非常微小，在周期中表现出“大-中-小”的嵌套结构；而泥沙量虽有不同时间尺度的周期变化，但结合分析6 a的变化最为明显，属于单周期变化。

3) 年径流量在31 a间未发生明显的突变，而年泥沙量则在2005年前发生较为显著的突变。突变原因与流域内部水利工程主要是黄河沙坡头水利枢纽工程有关，枢纽工程的完工令流域内部泥沙发生显著减少，放在年际变化上就表现为在2005年前突变为减少趋势。

4) 水沙的周期变化有着6 a的共同时间尺度，说明在6 a的尺度上泥沙量的变化是与径流量的变化有关的，但从长时间尺度上来讲，由于水土保持措施工程的实施与水利设施的运行，最终会导致泥沙的变幅比径流较大，人为因素占主导的情况下无法将泥沙量与径流量的变化直接联系。