雅鲁藏布江中游悬移质输沙率的周期性与不均匀性分析

杨超弈，陈 旻，张陵蕾，李 嘉，张耀文

（1.四川大学水力学与山区河流开发保护国家重点实验室，成都610065；2.四川城市职业学院建筑工程系，成都610110）

0 引言

雅鲁藏布江（以下简称“雅江”）是我国最大的高原河流，也是世界上海拔最高的大河之一，蕴藏9 000 多万kW 的水能资源，仅次于长江［1］，中游地区更是处于集西藏政治、经济、文化于一体的中心地带。然而雅江流域中游地区水土流失严重，流失面积约占流域面积的21.21%［2］，这是造成雅江内泥沙含量较高的主要原因，已有研究表明流域内年输沙量比黄河源多6.59%［3］，给下游水资源利用和生态保护带来了不确定性，因此在雅江中游开展泥沙输运的研究事关重大。

河流中悬移质泥沙含量通常占泥沙总量的90%［4］，其输移过程与土地覆盖、流域水文密切相关，通常被用作衡量流域内土地退化及水土流失的重要指标［5］。此外悬移质也是河流中生源物质的主要输运载体之一，是维系河流生态过程的重要因子，因此关于悬移质变化的研究一直广受水利工程、生态修复等领域的关注。目前，文安邦［6］等人发现了雅江在20世纪60-90年代悬移质输沙量呈现先减小后增大规律，张凡［7］等人则探究了悬移质输沙量在空间的变化规律，但受限地理环境等因素影响，完整悬移质输沙率资料难以获取，在流域并未开展更多深入的研究。而输沙率多时间尺度下的丰枯变化及年内分配方式将在不同程度上影响着水资源的开发利用和水库的调度运行方式，量化气候变化和人类活动对输沙率趋势变化的影响程度可为水利设施发展提供参考，因此在其他水电开发程度较大的流域上已有相关的研究开展，如郭文献［8］等采用小波分析探究长江中游输沙率在多时间尺度下的变化规律，郭彦［9］等人采用集中度和集中期的方法研究黄河内蒙古河段输沙率年内分配的不均匀性，孙倩［10］等采用累积量斜率变化率比较法量化了黄河中游降雨、人类活动对输沙率趋势变化影响的贡献程度。

目前，党中央关于制定“十四五”规划和二〇三五远景目标的建议中已经明确指出将在雅江进行水电开发，然而因为雅江流域内生态环境脆弱，为寻求合理的开发方式和运行方案，阐明雅江在天然状态下的输沙过程作为生态环境保护的参考基准十分必要。2000年后西藏迎来跨越式发展的“黄金时期”［11］，人均GDP 在20年间翻了近三番，相对而言在此之前雅江流域人类活动强度相对较小，尤其是河流中的水力连通关系基本保持为天然状态，因此这一时期内的河流过程具有较为显著的基线参考价值。

综上所述，本文基于雅江奴各沙、羊村、奴下水文站的1956-2000年悬移质输沙率（以下简称“输沙率”）、流量逐月数据，通过建立水沙关系曲线插补缺测数据，采用Morlet小波分析方法探究输沙率序列多时间尺度下丰枯变化规律，采用集中度和集中期方法分析输沙率在年内分配的不均匀性，采用累积距平法和累积量斜率变化率比较法分析输沙率变化趋势和影响因素对趋势变化的贡献程度，以期为未来流域内开展水土保持工作、水资源开发管理和河流水生态保护提供参考。

1 研究区域与方法

1.1 区域概况

雅江位于西藏自治区境内，地处北纬28°05′～31°16′，东经82°50′～96°40′，发源于西藏西南部喜马拉雅山北麓的杰马央宗冰川，在西藏境内全长2 057 km，流域面积24.048 万km2，分为上、中、下游三段。其中，中游段为里孜至米林县派镇段，该段地处温暖半干旱高原宽谷区段，该地区因地势平坦，土壤肥力较高，为西藏农业重点发展地区。中游段河长约1 340 km，均匀分布着奴各沙、羊村、奴下水文监测站。雅江中游河段水文站分布图见图1。

图1 雅鲁藏布江中游水文站分布图Fig.1 Distribution of the hydrological stations in the middle reaches of Yarlung Zangbo River basin.

1.2 数据来源

输沙率、流量数据均来源于《藏滇国际河流水文资料》，其测量方法严格遵循水文部门的标准准则。但受地理位置、技术条件、水文站搬迁等因素影响，在部分年、月存在缺测，缺测率约30%。研究［12］发现由流量与输沙率拟合的水沙关系曲线多呈QS=aQb幂函数形式，其中Qs为悬移质泥沙输沙率，kg/s；Q为流量，m3/s；a为系数、b为指数；刘家富［13］等人用水沙关系曲线对羊村站输沙率资料进行修正。为后文研究需要，本文采用水沙关系曲线将2000年以前缺测输沙率进行插补。

1.3 研究方法

1.3.1 小波分析

小波分析是一种研究非平稳随机时间序列的有效方法，而泥沙序列一般是非平稳的离散等间距序列，具有周期性［14］，因此，小波分析方法适用于泥沙序列的分析研究［15］。

对于给定的能量有限信号f(t) ∈L2(R)，其连续小波变换为：

式中：Wf(a，b)为小波变换系数；f(t)为一个信号或平方可积函数；a为伸缩尺度；b为平移参数。

Morlet 小波函数是常用的复值小波，其特点是具有良好的时、频域局部性［16］，因此在本研究中采用Morlet小波。其小波函数为：

式中：ω0为常数；i为虚数。

将时间域上关于时间尺度a的所有小波系数的平方进行积分，得到小波方差：

小波方差随时间尺度的变化过程称为小波方差图，能反映时间序列的波动能量随尺度的分布情况，通过小波方差图中的峰值对应的时间尺度，可以确定水文序列的主要周期，横坐标表示时间尺度，纵坐标是小波方差，小波方差越大，对应的时间尺度下水文序列波动越大，该尺度即为序列变化的主周期。

1.3.2 集中度与集中期法

集中度反映了输沙率在年内的集中程度，集中期反映了一年中最大输沙率出现的时间［17，18］。计算公式如下：

式中：RCDyear是输沙率集中度；RCPyear是输沙率集中期；Rx，Ry分别是12 个月的分量之和所构成的水平、垂直分量；ri是第i月的输沙率；θi是第i月输沙率的矢量角度。RCDyear越大表示输沙率分配越不均匀，各月份包含及代表的角度值见表1。

表1 RCPyear各月份包含和代表值Tab.1 The contains value and representative value of monthly RCPyear

1.3.3 灰色关联分析

灰色关联分析可用来判断不同序列之间的联系的紧密性［19］，在本文中将用来分析降雨量、径流量对含沙量的影响。

灰色关联度计算公式：

式中：γ0i为灰色关联度系数；ξ0i(k)为点关联系数。

把公式计算的各个水文要素进行排序，关联度系数越大表征自变量与因变量的相关性越密切。

1.3.4 累积距平法

累积距平法是一种判断水文序列离散数据点趋势变化的方法［20］，若累积距平曲线上升（下降），表示水文序列逐年增加（减小），曲线上由上升趋势变为下降趋势的点看作突变点，对应横坐标为突变年份，具体公式见下：

式中：LPi为第i年的距平累积值；Ri为第i年的输沙率；-R为输沙率的多年平均。

1.3.5 累积量斜率变化率比较法

该方法由王随继等［21］人提出，在量化自然因素和人类活动对输沙率趋势变化影响的贡献程度研究上广泛应用，以年份为自变量，以累积输沙率及累积降雨量为因变量，自然因素中常考虑仅有降雨量影响输沙率的趋势变化，基本公式见下：

式中：RS为累积输沙率斜率变化率，kg/（s·a）；KSb、KSa为累积输沙率-年份线性关系式在突变点前后两个时期的斜率，kg/（s·a）；RP为累积降雨量斜率变化率，mm/a；KPb、KPa分别为累积降雨量-年份线性关系式在突变点前后两个时期的斜率，mm/a；CP、CH分别为降雨量、人类活动对输沙率趋势变化的贡献率。

2 结果与分析

2.1 数据插补可靠性分析

根据奴各沙、羊村、奴下站多年月均流量和输沙率资料，点绘出流量与输沙率的关系图见图2。

图2 流量与输沙率关系图Fig.2 The relationship of flow and sediment transport rates

从空间分布来说，从奴各沙站往奴下站，系数a值先增大后减小，指数b值先减小后增大。研究表明，a、b具有实际的意义，a反映的是流域产沙特性，与水土保持情况及人类活动有关［22，23］。自上往下，3 个站点控制段林地面积依次为1.73 万、0.87万、3.07 万hm2；20年间由林地转化为建筑用地面积依次为300、1 650、150 hm2［24-26］，林地面积由小及大为：羊村＜奴各沙＜奴下，转化的面积由大及小为：羊村＞奴各沙＞奴下；因此对应系数a的大小：羊村＞奴各沙＞奴下。

指数b则反映了河流对于泥沙的搬运能力［22，23］，搬运能力通常与流速有关。虽然由上至下从奴各沙到羊村，受支流汇入和降雨影响，流量沿程增加，但由于羊村地处宽谷地带最宽处（宽约3 km），平均流速仅0.95 m/s，而奴下平均流速为1.05 m/s，奴各沙平均流速为1.44 m/s，即流速由大到小为：奴各沙＞奴下＞羊村，因此对应于奴各沙站的输沙能力最强，奴下站次之，羊村站相对最弱。

从水沙关系曲线的数据点拟合状况来看，羊村站表现为高流量处点比较分散，原因可能跟河床条件有关，奴各沙站、奴下站处河床由沙卵石、卵石组成，河床条件均较稳定，因此曲线拟合系数R2均大于0.94；而羊村站河床为沙夹卵石，河床不大稳定，冲淤变化较大。因此相较奴各沙站、奴下站，羊村站在高流量和高输沙率区间表现出相对离散特性。

综上所述，水沙关系曲线中系数与指数大小关系与实际情况相符，且拟合R2值均大于0.85，因此通过水沙关系曲线对缺测数据的插补具有合理性。

2.2 输沙率序列周期性分析

采用Morlet 小波变换对输沙率序列分析，结果见图3，可以看出3个水文站输沙率序列存在多时间尺度的变化周期。奴各沙站输沙率序列分别存在3～7、8～12、13～24、25～32 a 四类时间尺度的周期变化规律，除8～12 a 外的三类时间尺度的周期变化明显，在整个时间域内分布均匀；而8～12 a 的时间尺度仅在1956-1974年和1990-2000年间分布。同理羊村站输沙率序列分别存在3～7、8～12、13～23、25～32 a 四类时间尺度的周期变化规律，除8～12 a 外的三类时间尺度的周期变化明显且具有全域性；8～12 a 的时间尺度仅在1956-1970年和1995-2000年间分布。奴下站输沙率序列分别存在4～7、8～12、13～23、25～32 a 四类时间尺度的周期变化规律，且均具有全域性。

图3 输沙率小波等值线图Fig.3 Wavelet coefficients of sediment transport rates

3 个水文站输沙率小波方差曲线见图4。可以看出存在四处明显的峰值，第一大峰值对应的时间尺度为19 a，第二大峰值对应的时间尺度是26 a，第三大峰值对应的时间尺度是10 a，第四大峰值对应的时间尺度是6 a，对应输沙率第一主周期为19 a，第二主周期为26 a，第三主周期为10 a，第四主周期为6 a。

图4 输沙率小波方差图Fig.4 Wavelet variance of sediment transport rates

根据小波方差检验的结果，绘制出了控制段流域内输沙率变化的主周期小波变换系数图，结果见图5。以奴各沙站为例，可以看出输沙率在19 a的时间尺度下经历了3个低～高～低的交替循环，在26 a 的时间尺度下经历两个低～高～低的交替循环，在10 和6 a 的时间尺度下则经历更多的交替循环，羊村站和奴下站同理。在2000年末，3个水文站的19、26 a的特征时间尺度的小波系数为正值，因此流域在2000年后一段时间内存在一个多沙期，而这与Li［27］等人得出的“雅江中游2000年后的一段时间内将经历多径流期”相吻合。

图5 输沙率小波系数实部过程图Fig.5 Wavelet coefficients at different scales of sediment transport rates

综上所述可知，雅鲁藏布江流域中游的输沙率主要存在四个时间尺度的变化周期，第一和第二主周期分别为19 和26 a，属于较长的变化周期，而第三和第四主周期分别为10 第6 a，为雅江输沙率变化的主要短周期。

输沙率出现周期变化可能受到太阳活动和气候事件的影响。太阳活动主要是通过辐射向地球传递热能，由于垂向和纬向上热能收支的差异，造成了不同区域海表温度和气压的不同，因此出现了大气环流现象和各类气候事件［28］，而这些气候事件将通过影响大气环流的方式来影响水汽的输送，最终对各地降雨造成影响。降水主要从两方面影响输沙率：一方面通过对土壤的侵蚀作用从陆地携沙进入河流中，另一方面通过改变径流进而影响泥沙输移的变化，因此太阳活动和气候事件可以看作是输沙变化的遥相关因子，太阳活动与输沙关系见图6。对于雅江流域而言，已有研究［29-31］表明降雨量对ENSO 事件、北极涛动有一定响应，具体表现为：ENSO 事件通过加强（抑制）沃克环流，削弱（增强）南亚夏季风，最终减少（增加）了从孟加拉湾输送到西藏的水汽；类似地，当北极涛动的作用加强时，此时东亚大槽强度、西伯利亚高压、高原季风作用减弱，利于西南暖湿气流的加强，因此雅江降水出现6～10 的周期振荡［32］，这可能是受ENSO 的4～10 a 周期、北极涛动11 a 周期和太阳黑子活动平均11.07年周期的影响［33-35］。同样与之相符的是雅江的径流过程也存在4～6 和8～10 a 的周期变化［36］。雅江输沙率6 和10 a的变化周期基本与ENSO 和太阳活动的变化周期一致，因此可以认为全球性的气候事件在一定程度上影响着雅江中游输沙率的周期变化。

图6 太阳活动与输沙关系示意图Fig.6 The relationship of solar activity and sediment transport

2.3 输沙率序列年内分配不均匀性分析

采用集中度和集中期法对输沙率序列分析，结果见表2。根据年输沙率对应的角度，可以看出奴各沙站年内最大输沙率主要集中出现在8月，集中度对应的数值基本呈现先减小后增大的趋势，说明其最大输沙率对应的月份保持不变，仅在日期上先向前移动、后向后移动；羊村站和奴下站输沙率集中期在7月、8月交替出现，1961-1970、1971-1980年期间，羊村站与奴下站输沙率集中期均出现在7月。奴下站还在1991-2000年期间集中期出现在7月，其余时期集中期均发生在8月。从表1 可知，8月对应的角度为195°～225°，而奴各沙站、羊村站和奴下站发生在7月的集中期对应的数值与195°非常接近，说明集中期主要发生在7月末，接近8月。这与张小侠［37，38］等发现“6-8月的雨季降水占年总降水量的51%～80%，最大流量一般发生在7月或8月”相吻合。

表2 输沙率集中度与集中期Tab.2 Concentration frequency and concentration period of sediment transport rates

从奴各沙站到奴下站，集中期对应的数值在减小，说明从上游到下游，最大输沙率在日尺度上有向前出现的趋势，但在季节尺度上无明显向前或向后偏移的趋势；同时集中度也呈下降趋势，输沙率的集中程度沿程逐渐减弱，表明输沙率在年内分配上逐渐表现出均匀化。为探究输沙率在沿水流方向下降的原因，采用集中度计算方法对奴各沙、羊村、奴下站的降雨量和径流量序列进行计算，结果见表3，可以看出降雨量沿顺水流方向集中度降低，而径流量没有明显变化。降雨集中度下降的原因与水汽运动有关，雅江夏季降雨的水汽主要来源于印度洋的暖湿气流，受西南季风影响，水汽输送轨迹溯布拉马普特拉河而上后，经雅江下游向北移动，最后再在雅江大峡弯折往西北［39］，因此在雅江中游研究区域内，水汽依次到达奴下、羊村、奴各沙水文站，奴下的雨季开始的时间依次早于羊村、奴各沙水文站，据研究［40］表明雅江雨季的开始是下游早于上游约1～2个月，而雅江雨季的结束是下游晚于上游，这就造成降雨时期自上而下变长，降雨时期越长，降雨量在年内分配相对更均匀，因此沿河流方向，集中度大小为：奴各沙＞羊村＞奴下。而三个站点径流集中度没有明显变化，原因可能跟流域集水面积有关，奴各沙、羊村、奴下水文站的集水面积分别为10.637 8 万、15.319 1 万、18.984 3 万km2，集水面积越大，使得水体滞留时间变长，由降雨形成径流过程的调节作用就更强。虽然降雨的集中度在逐渐减小，但流域自上而下集水面积在增大，导致3个站点径流量的集中度并未出现明显变化［41，42］。因此可以推测输沙率在顺水流方向下降跟降雨有关。水力侵蚀是雅江泥沙的主要来源方式［6］，降雨、径流作为流域产沙的重要水力因子，也在一定程度上影响着河流中泥沙含量，通过计算含沙量和降雨量、径流量的灰色关联度，结果见表3，同一个站点内，降雨量与含沙量的关联度系数均大于径流量与含沙量的关联度系数，说明降雨对含沙量的影响更大。降雨量在年内分配得不均匀，导致河流含沙量的分配不均匀，进而影响输沙率的不均匀性。因此，可以初步判定雅江中游输沙率的集中度沿河流方向下降是因为降雨量集中度沿河下降。

表3 降雨量、径流量集中度和灰色关联度系数Tab.3 Concentration degree of precipitation、runoff and grey incidence coefficient

在长江八大支流之一的赣江流域上，水沙集中度基本吻合，最大输沙率出现在6月，与其雨季4-6月发生时间基本一致。受到水库修建运行影响，从20世纪60年代到70年代输沙率的集中度降低了16%［43］，但最大输沙率发生的时间在月尺度上基本上没变化，可见水库的调度运行对于集中度具有明显的改变作用，而集中期基本受气候影响。雅江干流在2000年前无水库工程的修建，因此输沙率的集中度、集中期主要受气候控制，而流域内输沙率的集中度、集中期受水库调节作用机制有待进一步研究。

2.4 降雨、人类活动对输沙率的影响

根据式（8）计算出1956-2000年输沙率累积距平的年际变化，结果见图7。奴各沙、羊村、奴下站的输沙率累积距平值均在1964年前呈现增大趋势，1964年后呈现减小趋势，即发生突变的年份为1964年。据此，可以将雅江中游流域输沙率的变化分为两个时期：1956-1964年为基准期，输沙率变化主要受降雨影响；1965-2000年为变化期，该时期输沙率的变化受到降雨和人类活动的影响。

图7 输沙率累积距平图Fig.7 Cumulative anomaly of sediment transport rates

各站点拟合关系情况见表4，从表4 中可以看出拟合的R2均大于0.9，且通过了显著性水平为0.01的Pearson显著性检验。奴各沙输沙率在基准期、变化期的斜率分别为571.9、262.3 kg/（s·a），变化期与基准期相比，斜率减小了309.6 kg/（s·a），减小率为54.14%，而降水量在基准期、变化期的斜率分别为45.33、35.85 mm/a，变化期与基准期相比斜率减小了9.48 kg/（s·a），减小率为20.91%，因此可以得到奴各沙站降水量的减小对输沙率减小的贡献率为38.62%，人类活动对输沙率减小的贡献率为61.38%。羊村站、奴下站同理结果见表5。可以看出3 个站点人类活动对于输沙率减小的贡献率均大于降雨变化的影响。雅江中游是农业重点发展的地区，该区农田面积占自治区农田总面积的60%以上［44］，据记载［45］，60年代耕地面积为21.25 万hm2，到20世纪90年代已达23.65万hm2，在30年间增长了2.4万hm2，为农业灌溉的需要，灌溉设施建设在60年代以后快速发展，60年代仅25 条千hm2、4 条万hm2水渠、4 座水塘在建，到80年代对应设施建设分别达66 条、10 条、16 座，数量上分别是60年代的2.64 倍、2.5 倍、4 倍。耕地面积的增加将会引发灌溉设施的兴建，这在一定程度上消耗水资源，间接影响输沙率的减小。而奴下站中人类活动的影响对输沙率减小的贡献率超过100%，主要是因为降雨的变化趋势是增加，而输沙率的变化趋势为减小，因此由降雨增加引起的输沙率减小的贡献率为-38.15%，而正是因为人类活动的作用，导致输沙率减小，因此人类活动对输沙率减小的贡献率为138.15%，超过100%。

表4 输沙率-年份、降雨量-年份拟合关系Tab.4 The fitting relationship sediment transport rate-year，precipitation-year

表5 降雨和人类活动对输沙率影响的分析结果Tab.5 Analysis of the effects of precipitation and human activities on sediment transport rates

3 结论

本文首先建立了雅鲁藏布江奴各沙、羊村、奴下水文站输沙率-流量的幂函数关系，对缺测数据进行插补，而后采用小波分析、集中度和集中期法、累积距平法、累积量斜率变化率比较法对雅江中游1956-2000年共45年的悬移质输沙率过程进行周期性分析、不均匀性分析、量化降雨和人类活动对输沙率趋势变化分析。结果表明：

（1）雅鲁藏布江中游输沙序列存在4 类时间尺度的变化特性，河段内主周期变化基本一致，第一主周期为19 a，第二主周期为26 a，第三主周期为10 a，第四主周期为6 a，输沙率的短周期变化与太阳活动、气候事件周期变化存在响应关系。

（2）雅江中游输沙率整体表现出较强的年内分配不均匀性，奴各沙站最大输沙率主要发生在8月，羊村站和奴下站的输沙率主要发生在7、8月。集中度在年代间变化不大，但呈现出从上往下游逐渐均匀化的规律，从奴各沙到奴下依次为0.83、0.77、0.74，这主要由于降水在局部区域之间差异影响含沙量的不均匀性造成。

（3）雅江中游输沙率在1956-2000年呈现先增长后减小的趋势，在1964年突变，由降雨影响导致输沙率减小的贡献率分别为：38.62%、41.04%、-38.15%，由人类活动影响导致输沙率减小的贡献率分别为：61.38%、58.96%、138.15%。□