金沙江下游水系分布复杂并且存在明显的水沙异源现象^[1-3]，量化分析子区间水沙变化对流域出口水沙变化趋势的影响为当前研究难点^[4-6]。金沙江下游水沙分布的空间差异较大^[7-10]，下游流域的区间径流量占总流域的18.80%，而输沙量则占总流域的63.30%。金沙江下游流域内雅砻江区域为丰水少沙区，输沙模数为350 t/km²；攀枝花至华弹区间为平水丰沙区，输沙模数为2462 t/km²；华弹至屏山区间为多水多沙区，输沙模数为2180 t/km²。近年来随着金沙江流域内修建众多大型水库，下游来水来沙的空间分布发生明显变化^[1-3]。已有研究表明^[11-13]，对金沙江下游减沙存在决定性作用的大型水库有二滩水电站(雅砻江)、金安桥水电站(金沙江中游)、金沙江下游干流梯级水电站，上述水库拦沙导致1998年以后下游年均输沙量减小38.55%。

为明晰气候变化及水库拦沙对金沙江下游水沙变化的影响，众多学者运用Mann-Kendall(MK)趋势检验、累积距平、双累积曲线等方法，分析下游水沙变化趋势及水沙突变特性^[14-16]。研究发现，气候因素导致金沙江下游水沙在1983、1995、2003年发生较大变化，其中1983、1995年对应丰水期，2003年对应枯水期；水库拦沙因素导致金沙江下游水沙在1998、2013年左右发生突变，分别对应二滩、溪洛渡与向家坝水电站运行使用。在气候变化及水库拦沙共同作用下，2013年以后流域内98%以上的泥沙均被拦截至库区。以往研究阐述了金沙江下游水沙分布的总体变化趋势，但水库运用对下游水沙分布的影响尚不明晰^[17-19]。金沙江下游存在突出的水沙异源现象。近年来大型水库拦沙导致金沙江下游年出口输沙量锐减，流域水沙异源现象加剧。本文基于金沙江下游干支流水沙观测数据，运用MK趋势检验法、R/S分析法、双累积曲线法，分析水库运用对金沙江下游水沙时空变化趋势的影响，计算不同区间水沙变化对下游出口的减沙贡献度，掌握金沙江下游水沙变化规律，为开展精细化的水库泥沙调控、减少水库泥沙淤积、优化梯级水库淤积分布奠定基础。

1 研究区域及方法 1.1 研究区域

金沙江下游流域研究范围为攀枝花至屏山或向家坝水文站区域，区间集水面积为14.70万km²。1954—2021年金沙江下游干支流分布有22个水文站，包括8个干流水文站及14个支流水文站，流域水系及水文站空间分布如图 1所示。流域出口控制水文站为向家坝站，由于向家坝站2008年才建立，水文资料相对较少，所以长系列资料采用其上游30 km的屏山水文站(在向家坝水电站建成后屏山站改为水位站)作为补充。金沙江下游入口(三堆子水文站)主要水沙来源为金沙江中游及支流雅砻江，分别以攀枝花站、桐子林站为控制水文站，代表金沙江中游与雅砻江水沙量。对于金沙江下游干流沿程，由于2008年以后华弹水文站更替为白鹤滩站，因此2008年前后分别以华弹、白鹤滩坝址的观测资料，代表金沙江下游中间位置的水沙量，并将两者统称为白鹤滩站。

1.2 研究方法

根据22个干支流水文站实测流量、含沙量等资料^[20]，采用MK趋势检验法^[21]、R/S分析法^[22]、双累积曲线法^[23]，依次研究金沙江下游水沙变化趋势、水沙突变持久性及水沙关系变异特性。

1.2.1 MK趋势检验法

MK趋势检验法首先基于流域水沙时间序列构造秩序列，然后定义统计量分析流域水沙变化趋势。

$U F_k=\frac{S_k-E\left(S_k\right)}{\sqrt{\operatorname{var}\left(S_k\right)}}$

(1)

$S_k=\sum\limits_{i=1}^k \sum\limits_j^{i-1} a_{i, j} \quad(k=2, 3, \cdots, n)$

(2)

$a_{i, j}=\left\{\begin{array}{lll} 1 & y_i>y_j \\ 0 & y_i<y_j \end{array} \quad(1 \leqslant j \leqslant i)\right.$

(3)

式中，S_k为流域水沙的秩序列，y为流域水沙系列值，n为水沙序列长度，i、j、k均为序列值的时间脚标。UF_k为定义的统计量，E(S_k)为秩序列S_k的均值，var(S_k)为S_k的方差。当UF_k>0说明流域水沙变化存在增长趋势，反之则存在下降趋势。由于UF_k满足标准正态分布，在给定显著性水平α的条件下，说明序列变化趋势存在显著性。为分析流域水沙变化突变年份，将流域水沙系列进行逆序排列，计算逆序系列的统计量UB_k。分析UF_k与UB_k的时间变化曲线，两者交点即为流域水沙突变点对应年份。

1.2.2 R/S分析法

R/S分析法是一种基于长序列相关性时间序列趋势分析方法，能够从分形时间序列中区分出随机序列和非随机序列。运用R/S分析法分别计算流域水沙系列的极差、标准差，将两者的比值作为指标参数。

$(R / S)_n=\frac{1}{A} \sum\limits_{a=1}^A\left(R_{\mathrm{I} a} / S_{\mathrm{I} a}\right) $

(4)

$R_{\mathrm{I} a}=\max \left(X_{k, a}\right)-\min \left(X_{k, a}\right)$

(5)

$S_{\mathrm{I} a}=\sqrt{\frac{1}{n} \cdot \sum\limits_{k=1}^n\left(N_{k, a}-e_a\right)^2}$

(6)

式中，n为系列的划分长度，A为子区间数目，(R/S)_n为指标参数，R_Ia为子区间极差，max(X_{k, a})为子区间最大值，min(X_{k, a})为子区间最小值，S_Ia为子区间标准差，N_{k, a}为子区间内第k个序列值，e_a为子区间平均值。Hurst经验公式中指标参数R/S与区间长度n存在指数关系。

$R / S=(c n)^H$

(7)

式中，H为Hurst指数。若H=0.5，流域水沙系列为随机序列；0≤H＜0.5，流域水沙为反持久性序列，即序列前后时间段的变化趋势相反；0.5＜H≤1，流域水沙系列具有状态持续性；H越接近1，序列的持续性程度越强，反之越弱。

1.2.3 双累积曲线法

双累积曲线法首先计算流域水沙系列随时间的累积量变化：

$Q_k =\sum\limits_{t=1}^k Q_t(k=1, 2, 3, 4, \cdots, n)$

(8)

$S_k =\sum\limits_{t=1}^k S_t(k=1, 2, 3, 4, \cdots, n)$

(9)

式中，Q_k、S_k为流域水沙系列在第k年的水、沙累积量。双累积曲线法通过观测流域水、沙累积值的线性斜率变化，以此确定流域水沙关系的突变特性。由于大多数流域水沙关系满足线性或幂函数分布特性，因此流域水、沙累积量之间满足线性关系，而双累积曲线的拐点则对应流域水沙突变年份。

2 金沙江下游水沙时空分布特性 2.1 年际水沙过程

1954—2021年金沙江下游干流年水沙量变化过程如图 2所示。以向家坝站水、沙数据代表金沙江下游总水、沙量，流域年均径流量为1434.97亿m³，年均输沙量为2.11亿t。金沙江中游年均来水量568.68亿m³，雅砻江年均来水量590.05亿m³，两者与向家坝站年径流量占比为80.75%；中游年均来沙量4276.96万t，雅砻江年均来沙量3141.61万t，两者与向家坝站年输沙量占比为35.17%。分析金沙江下游出口水沙变化特性，年径流量过程不存在显著趋势性变化，而年输沙量则在1998年及2012年以后显著下降。对比各水文站输沙量变化过程可知，1998年以后雅砻江年输沙量由4352.77万t减小为1375.33万t，降低68.40%；2012年以后金沙江下游出口(向家坝站)年输沙量发生断崖式减少，而白鹤滩站的年输沙量不发生显著变化，说明流域出口年输沙量减少的主要原因为向家坝水电站运行。

2.2 年内水沙过程

1954—2019年金沙江下游干流各水文站的月均流量及输沙率的年内变化过程，如图 3所示。金沙江下游出口(向家坝站)各月流量均值为4371.38 m³/s，变化范围为1430.83~9708.89 m³/s；各月输沙率均值为6526.28 kg/s，变化范围为87.86~21644.74 kg/s。金沙江下游干流各水文站的年内水沙分布均集中于汛期(6—10月)，向家坝站汛期月均流量为7726.61 m³/s，月均输沙率为15016.50 kg/s；而非汛期月均流量仅为1974.78 m³/s，为汛期月均流量的25.56%，非汛期月均输沙率为461.83 kg/s，为汛期月均输沙率的3.08%。对比金沙江下游水、沙年内分布特征，金沙江下游干流沿程水沙分布在汛期的差异均大于非汛期，白鹤滩站在汛期的月均输沙率为5706.26 kg/s，仅为向家坝站月均输沙率的38.00%，说明白鹤滩—向家坝区间为金沙江下游的主要产沙区域。

2.3 水沙空间分布

分析1954—2021年金沙江下游年均水沙量的空间分布特性，结果如图 4所示。金沙江下游出口年径流量主要来源为中游及雅砻江，下游区间内(除雅砻江)已知支流的年均来水量为121.64亿m³，未控区间的来水量为154.59亿m³，占总径流量的10.77%。金沙江下游出口年输沙量的主要来源为金沙江中游、雅砻江、缺少输沙量观测的区间。下游区间中已知支流的年均来沙量为3517.80万t，占总输沙量的16.69%，输沙量未控区间的年均来沙量为9943.17万t，占总输沙量的48.12%。

3 金沙江下游水沙变化趋势 3.1 年径流量变化趋势

利用MK趋势检验法、R/S分析法，分析1954—2021年金沙江下游年径流量变化趋势，结果如图 5所示。金沙江下游不同区间的年径流量变化过程不存在显著变化，年均水量变化的Hurst指数均小于0.5，说明各区间年均水量变化均为非持久性趋势变化。金沙江下游出口(向家坝站)年径流量在2000年以前为下降趋势，2000—2012年为上升趋势，2012年以后又变为下降趋势，年径流突变时间为1987、1994、2005、2016年；中游(攀枝花站)来水量在1989年前后由下降趋势变为增加，突变时间为1986、2016年；雅砻江(桐子林站)来水量在1998年前基本为下降趋势，1998—2010年为上升趋势，2010年以后又变为下降趋势，径流突变时间为1966、1989、1994、2005年。分析金沙江中游与雅砻江对下游出口水量变化的影响，金沙江下游出口水量变化趋势与雅砻江基本一致，1987、2016年下游出口径流突变是中游水量变化的结果，而1994、2005年下游出口年水量突变由雅砻江水量变化引起。由于雅砻江与金沙江中、下游的年径流量变化趋势均在显著限以内，并且年径流量趋势变化与水库建设进程无对应关系，说明水库运用对金沙江下游年径流量变化不存在显著影响。

3.2 年输沙量变化趋势

研究1954—2021年金沙江下游年输沙量变化趋势，结果如图 6所示。金沙江下游出口年输沙量变化趋势在2013年以后超过显著限并且Hurst指数为0.16，表明金沙江下游年输沙量显著减小是由短期内人为活动导致。总变化趋势方面，金沙江下游出口沙量在1984年前为下降趋势，1984—2007年为上升趋势，2007年以后为下降趋势，尤其是2013年以后发生显著下降，2013年前后其年均输沙量由2.36亿t骤减为147.52万t；中游来沙量在1985年前为下降趋势，1985—2014年为增加趋势，2014年以后为下降趋势，输沙量突变时间为2017年左右；雅砻江来沙量在2003年前基本为增加趋势，2003年以后变为下降趋势，尤其2008年以后显著下降，输沙量突变时间为2007年。

分析下游入口来沙量变化对流域出口年输沙量变化趋势的影响：① 1984年以前，金沙江下游出口输沙量变化趋势与中游保持一致；② 1984—2007年，下游出口年输沙量总体为增加趋势，而2003年以后雅砻江输沙量变为减小趋势，说明2003年以后中游来沙量的增加趋势大于雅砻江的减沙作用；③ 2007—2013年，下游出口输沙量总体为下降趋势，这与雅砻江输沙量在2007年发生显著下降密切相关；④ 2013年以后，金沙江下游出口输沙量发生趋势性突变，并且在2013年以后显著下降。

3.3 水沙关系突变特性

为揭示水库运行对金沙江下游水沙关系的影响，运用双累积曲线法分析1954—2021年金沙江下游水沙关系的突变特性，结果如图 7所示。金沙江中游(攀枝花)、雅砻江(桐子林)的水沙累积曲线，分别在2010、1998年出现明显降低趋势，同时金沙江下游出口也存在下降趋势，但是下游出口水沙关系显著变化则发生于2013年。利用金沙江下游含沙量变化验证干、支流水沙突变特性，金沙江中游含沙量在2010年以后由0.89 kg/m³降低为0.16 kg/m³，雅砻江含沙量在1998年以后由0.79 kg/m³降低为0.22 kg/m³，而金沙江下游出口含沙量在1998、2010、2013年以后，由1.82 kg/m³依次降低为1.31、0.86、0.01 kg/m³。根据金沙江下游水沙变化过程分析下游出口含沙量的变化规律，1954—2021年金沙江下游年径流量不存在趋势性变化，而年输沙量则在1998、2012年显著减小^[20]，说明流域含沙量降低主要受输沙量变化影响。

结合金沙江下游水、沙变化趋势，分析金沙江下游水沙关系曲线发生明显降低趋势的主要原因。由于金沙江中游及雅砻江的来水量分别在2010、1998年以后为增加趋势，而同期内两者来沙量则为减小趋势，表明其水沙关系突变均由人为活动造成。在金沙江中游和雅砻江不发生水沙突变的条件下，金沙江下游出口径流量在2012年以后不存在明显趋势变化，而下游出口含沙量大幅下降，说明金沙江下游出口水沙关系突变也由人类活动影响导致。依据金沙江流域大型水库建设进程可知，雅砻江、金沙江中游、金沙江下游出口水沙关系突变，分别对应二滩水电站、金安桥水电站、溪洛渡水电站运行。

4 水库运用对金沙江下游水沙变化影响 4.1 干支流水沙关系变化规律

运用幂函数形式的水沙关系式^[24-26]，量化研究不同时期内金沙江下游不同区间的水沙变化规律，结果如图 8所示。

$W_{\mathrm{s}}=a \cdot Q^b $

(10)

式中，W_s为年输沙量(亿t)，Q为年径流量(亿m³)，a为系数，b为指数。由图 8可知，金沙江下游输沙量显著减小时，各水文站的水沙关系曲线均出现明显下凹，即相同径流条件下产沙量明显降低。1998—2009年雅砻江水沙关系突变，年均含沙量下降70.76%，白鹤滩站与向家坝站年均含沙量相比1998年以前分别下降16.71%、24.59%；2010—2012年金沙江中游水沙关系突变，其年均含沙量减少82.53%，白鹤滩站、向家坝站年均含沙量相比1998年以前分别下降46.68%、50.65%。

运用式(10)拟合不同时期的金沙江下游年际水沙关系，拟合参数见表 1。结合无量纲分析法及幂函数关系参数，分析金沙江下游水沙关系突变规律。

$\frac{W_{\mathrm{s}}^{\prime}-W_{\mathrm{s}}}{W_{\mathrm{s}}}=\frac{a^{\prime}}{a} \cdot Q^{b^{\prime}-b}-1$

(11)

式中，W_s、W′_s分别为水沙关系变异前后的年输沙量(亿t)，Q为年径流量(亿m³)，a′、a为水沙变异前后的水沙系数，b′、b为水沙指数。上式左侧为水沙关系突变后年输沙量的相对变化量，其比值越接近于0，表明水沙突变对输沙量的影响越小；等式右侧中a′/a与b′-b分别为水沙关系系数、指数的相对变化量，两者的值越接近于1和0，表明水沙关系变化程度越小。结合式(11)分析水沙关系突变对输沙量变化的影响，当水沙关系系数、指数同步增大或减小时，水沙关系突变对输沙量影响越大；反之，当水沙关系系数与指数的变化规律相反时，水库拦沙对输沙量的影响则受水沙关系限制。表 1中金沙江下游各水文站的水沙关系系数与指数变化始终保持负相关关系，即a′/a > 1时b′-b＜0，反之a′/a＜1时b′-b>0，由此表明水库拦沙条件下金沙江下游水沙关系存在自我调节作用^[27]。

4.2 干流含沙量沿程变化规律

通过分析金沙江下游沿程含沙量变化，表征下游干流水沙量沿程演变规律，结果如图 9所示。结合金沙江中游及雅砻江来水来沙条件，计算金沙江下游干流入口(三堆子站)的含沙量。

$S_0=\frac{W_{\mathrm{s} 1}+W_{\mathrm{s} 2}}{Q_1+Q_2}$

(12)

式中，S₀为下游入口含沙量(kg/m³)，W_s1、W_s2为金沙江中游、雅砻江年输沙量(亿t)，Q₁、Q₂为金沙江中游、雅砻江径流量(亿m³)。图 9中金沙江下游沿程相邻断面的含沙量均满足正相关关系，并且各时间段内三堆子至白鹤滩的含沙量变化率始终高于白鹤滩至向家坝段，说明三堆子至白鹤滩区间为下游来沙的主要区间。

运用线性函数拟合金沙江下游干流沿程含沙量关系，1998年以前三堆子至白鹤滩站的含沙量变化率为1.84，白鹤滩至向家坝站的含沙量变化率为1.21；1998—2009年雅砻江输沙率减小，下游沿程2个区间的含沙量变化率分别减小为1.75、1.02；2010—2012年金沙江中游水沙突变，下游沿程区间的含沙量变化率增加为2.02、1.15；2013年以后，泥沙绝大部分在4座梯级水库发生淤积，白鹤滩水文站和向家坝水文站含沙量都很小。

4.3 区间减沙对流域减沙贡献度

基于金沙江下游水沙关系及沿程水沙变化规律，计算不同时期内金沙江下游减沙率。首先将1998年和2010年以前作为金沙江中游及雅砻江水沙变化的基准期，并且根据1998年以前金沙江下游沿程含沙量关系，分别得到基准期内白鹤滩和向家坝站的含沙量表达式。

$S_1=1.84 \cdot \frac{a_1 \cdot Q_1^{b_1}+a_2 \cdot Q_2^{b_2}}{Q_1+Q_2} $

(13)

$S_2=1.21 \cdot S_1$

(14)

式中，S₁、S₂分别为白鹤滩、向家坝站在基准期的含沙量(kg/m³)，a₁、a₂为金沙江中游、雅砻江的水沙系数，b₁、b₂为金沙江中游、雅砻江的水沙指数，Q₁、Q₂为金沙江中游、雅砻江的年径流量(亿m³)。根据变化期与响应期的含沙量差异，计算不同时期内金沙江下游的总减沙率。

$\lambda=-\frac{W_{\mathrm{s}}-S \cdot Q}{S \cdot Q}$

(15)

式中，λ为金沙江下游(向家坝站)在响应期的总减沙率，Q、S为向家坝站在基准期的年径流量(亿m³)、含沙量(kg/m³)。运用式(15)计算金沙江中游、雅砻江、下游干流水库对流域的减沙作用，1998—2009年雅砻江水沙关系曲线发生明显降低趋势，对应流域总减沙率为35.33%；2010—2012年，在金沙江中游及雅砻江水沙变异条件下，下游总减沙率增加为51.44%；2013—2019年，金沙江下游干流陆续修建溪洛渡、向家坝水电站，下游总减沙率增至99.34%。为分析干支流水沙变化对下游减沙的影响，利用沙量守恒定律计算各区间对总减沙量的贡献度。

$\varphi(k)=\frac{W_{\rm{s}}^{\prime}(k)-W_{\rm{s}}(k)}{W_{\rm{s}}^{\prime}-W_{\rm{s}}}$

(16)

式中，φ(k)为k区间对金沙江下游流域的减沙贡献度，W′_s(k)为k区间的天然年输沙量(万t)，W_s(k)为区间k年输沙量观测值(万t)。运用式(16)计算不同区间对金沙江下游的减沙贡献度，结果见表 2。1998—2009年，金沙江下游出口年均减沙量为9749.70万t，相对1998年以前总减沙率为35.16%，其中由于二滩水电站的运用雅砻江年均减沙量为3206.01万t，对流域出口减沙贡献度为32.88%；由于水土保持措施和支流水库拦沙，三堆子—白鹤滩区间年均减沙量为4014.59万t，对流域减沙贡献度为41.18%；白鹤滩—向家坝区间年均减沙量为2529.10万t，对流域减沙贡献度为25.94%。2010—2012年金沙江下游出口年均减沙量1.05亿t，总减沙率为49.97%，其中由于水库拦沙作用金沙江中游年均减沙量为2253.37万t，对流域减沙贡献度为20.29%，雅砻江年均减沙量为2416.07万t，减沙贡献度为23.10%；水土保持措施和支流水库拦沙条件下，三堆子—白鹤滩区间年均减沙量为3207.90万t，对流域减沙贡献度为30.67%；白鹤滩—向家坝区间年均减沙量为2582.05万t，对流域减沙贡献度为24.69%。2013—2021年流域出口年均减沙量为24200.00万t，总减沙率为99.37%，其中三堆子以上区间受水库拦沙作用影响，金沙江中游年均减沙量为4747.21万t，对流域减沙贡献度为19.64%；雅砻江年均减沙量为3557.77万t，对流域减沙贡献度为14.72%；受水土保持措施和干支流水库拦沙影响，三堆子—白鹤滩区间年均减沙量为2984.05万t，对流域减沙贡献度为12.34%；在溪洛渡、向家坝水库拦沙条件下，白鹤滩—向家坝区间年均减沙量为12900万t，对流域减沙贡献度为53.30%。

5 结论

本文基于金沙江下游干支流水沙观测资料，利用MK检验、R/S分析法以及双累积曲线法，分析水库运用对流域水沙变化趋势的影响。通过分析干支流水沙关系和不同区间水库拦沙对流域出口减沙贡献度，得到如下结论：

1) 1954—2021年金沙江下游年径流量过程不存在显著趋势性变化，多年年均径流量为1434.97亿m³；年输沙量则在1998年和2012年以后显著下降，1954—1998年流域年均输沙量为2.55亿t，1999—2012年流域年均输沙量降低为1.73亿t，2013—2021年流域年均输沙量减小为0.015亿t。

2) 1954—2021年金沙江下游出口年输沙量大幅减小主要受大型水库运行的影响，雅砻江和金沙江中游、金沙江下游出口水沙关系分别在1998、2010、2012年发生突变，分别对应二滩水电站、中游梯级水库、下游梯级水电站运行，导致下游出口含沙量由1.82 kg/m³依次降低为1.31、0.86和0.01 kg/m³。

3) 二滩水电站运行后，1998—2009年金沙江下游出口总减沙率为35.16%，雅砻江的减沙贡献度为32.88%；金安桥等中游水电站运行后，2010—2012年下游出口总减沙率为51.43%，金沙江中游、雅砻江、三堆子至白鹤滩区间、白鹤滩至向家坝区间对下游减沙贡献度分别为21.54%、23.10%、30.67%、24.69%；金沙江下游梯级水库运行后，2013—2021年流域总减沙率为99.37%，不同区间对总减沙量的贡献度分别为19.64%、14.72%、12.34%、53.30%。