基于模态分解法的长江与两湖水沙特征分析

黄莹华 万迪文 丁玉堂

摘要：为研究长江和洞庭湖、鄱阳湖的水沙变化，利用模态分解方法（Welch能谱法及本征正交分解POD法），对长江和两湖中19个主要水文站1956～2018年的水沙序列变化特征进行分析。结果表明：长江和两湖的径流序列中，存在明显的8 a和20 a周期震荡。利用POD法分析得到的长江和鄱阳湖的径流变化率分别为2%～10%和10%～20%;而输沙变化率都接近90%。洞庭湖的水沙变化率要远远大于长江和鄱阳湖，其中荆江三口水沙变化率在140%以上。

关 键 词：

径流量; 输沙量; 模态分解; 长江; 洞庭湖; 鄱阳湖

中图法分类号： TV142

文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2021.06.006

0 引 言

流域系统是人类赖以生存的重要场所，具有调节区域气候、记录区域环境变化和维持区域生态平衡的功能[1]。流域系统中最为敏感的是水沙的变化，其间接地体现气候变化和人类活动对流域系统的影响，分析流域水沙的时空特征对河势稳定、河道演变、江湖关系等具有重要作用[2-3]。

长江流域的流域面积为180万km2，约占我国陆地面积的1/5。作为中华文明的发祥地，长江在中国的历史、文化和经济中具有举足轻重的地位。中国两个最大的淡水湖（洞庭湖和鄱阳湖）与长江直接相连，它们在维护长江流域水资源可持续性、水生生态系统健康及当地经济发展方面发挥着重要作用。

国内外对长江流域水沙的变化过程己有大量的研究。张强等[4]利用Mann-Kendall方法，分析了近40 a（1963～2004年）长江流域主要水文站点的流量及输沙率的变化特征。许全喜和童辉[3]利用长江干、支流28个主要水文站1956～2010年径流量和输沙量实测资料系列，较为系统地研究了近50 a来长江流域不同河段、不同时段的水沙变化特性。王海斌[2]采用Mann-Kendall检验方法和累积距平方法分析了长江流域干流6个水文站（1956～2009年）年径流量和年输沙量的变化情况。覃红燕[5]利用洞庭湖水沙1951～2009年序列，分析了洞庭湖水沙特征，湖南四水水沙演变规律及其成因。李微等[6]利用鄱阳湖入、出湖1956～2011年水沙系列数据，分析水沙变化特征及其原因。彭俊[1]对鄱阳湖入流五河水文站的年径流量和年输沙量序列（1950～2012）的变化进行了分析，并探讨了水沙变化的原因。邴建平[7]分析了长江中游干流和鄱阳湖的江湖水位关系、湖泊调蓄洪水能力等水系统要素的长历时时空演变特征及趋势。魏丽等[8]根据三峡水库蓄水前后水沙资料，对长江上游水沙特性进行了分析，结果表明上游水文站近期年径流量变化不大，年输沙量有减小的趋势。毛海涛等[9]基于2002～2016年期间三峡水库上、下游的11所水文站资料，分析了三峡水库的冲淤特性，探究了泥沙特性变化的原因。Zhang等[10]和Zhao等[11]认为宜昌站、汉口站和大通的年径流量没有显著趋势，但输沙量变化却很不相同。Zhang等[12]使用Mann-Kendall方法分析了鄱阳湖的入流与出流，发现入流在1993年前呈下降趋势，而此后呈上升趋势。Dai等[13]发现宜昌站、汉口站和大通的年径流量序列突变都发生在2000年左右。Wang等[14]分析了寸滩站和宜昌站日径流量序列，发现了周期为1周到1 a的震荡。

尽管已有许多研究讨论了长江流域径流量和输沙的变化趋势，但是采用的方法主要是定性描述、线性回归和Mann-Kendall法，基于这些方法得到结果主要是定性的，不能定量提取隐藏在时间序列中的特征模态。此外，这些论文中的水文站个数较少，不能全面反映长江流域水沙状况。因此，本文收集了长江流域21个水文监测站（长江5个，洞庭湖10个和鄱阳湖6个，如图1所示）在1956～2018年期间的年径流量和年输沙量序列，并利用模态分解方法（Welch能谱法及本征正交分解POD法），对水沙序列变化特征进行了分析。

1 研究区域及方法

1.1 研究区域

本次研究收集到长江和两湖共21个水文站的年径流量和年输沙量序列（1956～2018年），其中长江共5个，上游3个（朱沱、北碚、武隆），中游1个（宜昌），下游1个（大通）;洞庭湖共10个，其中四水4个（湘潭、桃江、桃源及石门），荆江三口5个（松滋口2个，太平口1个，藕池口2个），洞庭湖出口1个（城陵矶）;鄱阳湖6个，其中五河5个（外洲、万家埠、李家渡、梅港及虎山），鄱阳湖出口1个（湖口）。为下文叙述方便，设x（n）（n=1，2，3，…，N，N=63）为某一变量的时间序列（径流或输沙）。

1.2 数据分析方法

本文共使用了3种数据分析方法（M-K趋势检测法、能谱及POD方法）提取时间序列中的特征。M-K趋势检测法是统计检测方法，用来检测数据增减趋勢，并将无显著变化的数据用于能谱分析;能谱和POD法是模态分解方法，能谱是基于正交三角函数空间的模态分解，而POD法是基于信号序列自相关系数矩阵的特征向量空间的模态分解。

M-K方法为非参数性检测，且稳定性高，目前被广泛用于检测时间序列（如温度和降雨等）的趋势。M-K方法中的S变量定义为

S=∑Nn=2∑n-1j=1sign（xn-xj）（1）

式中：sign为符号函数;S是所有递增和递减趋势的总和，它反映了时间序列的总变化趋势。

当样本数量大于8时，S变量的分布近似正态分布，且其统计值us与σ2s如下：

μS=0，σ2S=NN-12N+15/18 （2）

式中：μs为正态分布的期望;σ2s为正态分布的方差。

由S变量估计的趋势置信度可以由变量Z估计，其计算值如下：

Z=S-1/σS，  S>00     ， S=0S+1/σS， S<0（3）

当Z>0时表示序列存在递增趋势，当Z<0时表示序列存在递减趋势。检测趋势的置信水平，当Z>Z1-p/2时，可以拒绝零假设，表示时间序列递增或递减的置信水平为p，其中Z1-p/2是标准正态离差，可以由标准正态分布表查得。当p=0.05时Z0.975=1.96，当p=0.01时Z0.995=2.58。

能谱计算最常用的是Welch法，通过添加窗函数、将数据分段并使各段之间有重叠，达到改善能谱曲线的光滑性及减小方差的目的。设Welch法分段计算的数据长度为M，每段数据间重合度为O1（0≤O1≤1），总计算段数L=（N-MO1）/（M-MO1），采用的窗函数为w（j）（1≤j≤M），采样频率为Fs，则Welch平均能谱公式为

P（f）=1FsMWLLl=1Mj=1xl（j）w（j）e-ifj2 （4）

式中：W=Mj=1w（j）2/M;Fs/M为能谱分辨率。

将Welch能谱无量纲化，得归一化能谱PSD（f）=10·lgP（f）/maxP（f），dB。

Proper Orthogonal Decomposition（簡称POD）为本征正交分解，与Karhunen-Loeve分解、主成分分析及奇异值分解等类似，主要用于提取序列中的含能模态。POD方法将信号序列投影到最优函数空间，使得信号在该最优函数空间投影的平均能量比其他函数空间（如Fourier基空间）都大，而该最优函数空间的解为信号序列的空间相关系数矩阵的特征向量空间，可通过特征分解进行求解，具体步骤如下。

将序列x（n）分解成若干相互重叠的小段，组成新矩阵A（X，tl），如下：

A（X，tl）=x（1）x（2）…x（Ns）x（2）x（3）…x（Ns+1）x（M）x（M+1）…x（N）（5）

式中：X表示每一个分段，每段长为M;tl表示分段的排列顺序，总共有Ns个分段，且Ns=N-M+1。A（X，tl）在特征向量空间的POD分解公式如下：

A（X，t）=Mk=1ak（t）Ψk（X）

ak（t）=[A（X，t），Ψk（X）]（6）

式中：Ψk（X）为A（X，t）的时间相关系数矩阵的特征向量;akt为信号序列A（X，t）在Ψk（X）上的投影系数。

Ψk（X）的计算公式如下：

Cs=AAT

[Φ，Λ]=eigen（Cs）

Ψ=（A·Φ）·Λ-1/2  （7）

式中：Cs表示A（X，t）的时间相关系数矩阵;eigen表示对Cs矩阵计算特征值对角阵Λ及特征向量矩阵Φ;对角阵Λ中的特征值λk（1≤k≤M）按从大到小排列，代表了时间序列在对应特征向量（第k阶模态）上投影能量的大小。

序列总能量、投影系数及特征值的相关关系为

E=Mk=1Lt=1ak（t）2=Mk=1λk（8）

式中：takt2是在Ψk（X）上投影能量的总和。

令第k阶模态的含能比例Ek=λk/E。

2 结果分析

2.1 M-K趋势分析

图2给出了年径流量和年输沙量序列的M-K趋势分析结果。对于长江上、中、下游的年径流量而言，不存在显著的变化趋势，Z值都位于95%的置信区间以内;但相比下游而言，长江上、中游的径流存在轻微的递减趋势。与年径流量不同，除湖口外，其他各位置的年输沙量存在急剧递减趋势，且置信度均高于99%。长江流域近半个世纪以来，由于气候变化及人类活动导致的泥沙输移量骤减情况甚为显著。

对于洞庭湖径流而言，仅湘、资、沅、澧四水不存在明显变化，其他位置均存在显著递减趋势，其中荆江三口和城陵矶的递减置信度均高于99%。对于洞庭湖的输沙而言，四水、荆江三口及城陵矶存在非常明显的递减趋势。

对于鄱阳湖的径流而言，五河（赣江、抚河、信江、饶河、修河）的变化趋势减小，而湖口的径流增加趋势较为显著。对于鄱阳湖的输沙而言，仅五河存在显著递减趋势，而湖口几乎不存在显著变化趋势。

2.2 能谱分析

长江和两湖的年径流量及输沙的能谱计算结果见图3，为清晰显示，已将横坐标的频率参数转化成了对应的周期。图3中仅给出了图2中不存在明显变化趋势的径流量和输沙量序列的能谱值，其中Welch计算采用的分段长度M=40 a，各分段间的重叠系数ol

=97.5%，窗函数采用hamming窗。图3给出了长江上、中、下游径流的能谱图，可见，长江上、中、下游径流中存在明显的8 a周期震荡。与长江类似，洞庭湖的入流四水及鄱阳湖的五河及湖口的年径流量中也存在较为明显的8 a周期震荡，但除此之外，它们还存在明显的20 a周期震荡。对于鄱阳湖湖口的年输沙量序列而言，其中存在明显的4.5 a及40 a周期震荡。

Mote等[15]认为，径流中的多年及几十年周期震荡特性与全球性气候变化相关，如南半球厄尔尼诺现象的周期为2～7 a，而太平洋气候震荡的周期为20～30 a。对比可知，长江和两湖的年径流量的变化周期正好位于全球性气候变化的周期内，也许全球性气候的周期变化引发了长江和两湖径流的周期变化。

2.3 POD分析

鉴于图3中径流量序列的最大周期为20 a，本节POD计算采用的分段序列长度设为20 a，各分段间隔为1 a，总段数为44。POD各阶模态的含能比例见图4，为清晰显示，仅给出前15阶模态，前3阶模态的具体含能百分数见表1。无论是年径流量序列还是年输沙量序列，其POD模态的含能比例都随着模态阶数的增大而减小，且绝大多数的能量集中在一阶模态。年径流量序列的一阶模态含能在90%以上，而其他模态含能都小于1%;年输沙量序列的一阶模态含能在84%以上，而其他数模态含能都小于2%;由一阶模态含能比例数值的高低推测可知，年径流量序列的有序性要高于年输沙量序列。

由于一阶模态包含了时间序列的绝大多数能量（大于85%），一阶模态可以优良近似原始序列，故分析一阶模态即可推得原始序列的趋势。归一化的一阶模态见图5。长江上、中游，荆江三口及城陵矶的年径流量存在较为明显的下降趋势;而长江下游、洞庭四水、鄱阳五河及湖口的年径流量却存在上升的趋势。对于年输沙量序列，除去湖口外，其他位置均存在明显的下降趋势。值得指出的是，图5中的一阶模态的斜率并不代表真实的变化率，必须乘以特征值才能得到真实值。

用最小二乘法拟合出图5中一阶模态的斜率，结合一阶模态的特征值得出年径流量及年输沙量序列的年变化率和总变化率（见表2）。联合图4及表2可知，荆江三口及城陵矶的年径流量急剧减少，年递减率分别为-19.7亿m3及-15.8亿m3，而洞庭四水年径流量却以每年0.7亿m3的速度增长;长江中上游径流的年递减率约为-5亿m3/a，但长江下游径流量以每年4.8亿m3的速度增长;鄱阳湖的五河和湖口都存在径流增加趋势，年递增率分别为2.1亿，5.3亿m3。

通过年变化率可以计算序列63 a的总变化率，即63 a的总变化值除以63 a的均值。由表2可知，长江的年徑流量序列在过去的63 a时间里变化很小，总变化率的绝对值都小于10%;鄱阳五河的径流变化也在10%左右，而湖口变化率约20%;洞庭四水总变化率很小，但城陵矶径流递减了30%，而荆江三口径流递减了近140%。

结合图4及表3可知，长江输沙量的递减率约为-600万t/a，远远高于洞庭湖和鄱阳湖的输沙递减率，但荆江三口的输沙递减率（约-380万t/a）相对较大。63 a来，长江上中下三游的输沙总变化率可达-100%，说明来沙条件发生了极为剧烈的变化;与之类似，鄱阳湖五河的输沙总变化超过-90%;洞庭湖的输沙总变化率更为显著，尤其是荆江三口，其总变化可达-230%。由此可知，自20世纪90年代起，受长江上游骨干水库群的逐步建成运行与拦沙、人工采砂、水土保持减沙、及流域降雨量时空分布改变的影响，导致了长江和两湖输沙量急剧减少。

3 结 论

利用统计检测和模态分解方法对长江和两湖水沙序列进行了特征分析，主要结论如下：

（1） 在1%的显著水平下，对于年径流量而言，除荆江三口与城陵矶外，长江、洞庭湖及鄱阳湖都不存在明显改变;对于年输沙量而言，除湖口外，长江与两湖的年输沙量总量呈显著下降趋势。

（2） 对长江和两湖中不存在明显增减趋势的年径流量与输沙量序列进行能谱分析可知，年径流量序列中存在明显的8 a和20 a周期震荡，而湖口的年输沙量序列中存在明显的4.5 a周期震荡，这些震荡应该与南半球厄尔尼诺现象和太平洋气候震荡相关。

（3） 由POD法计算结果可知，年径流量序列的一阶模态含能在90%以上，而年输沙量序列的一阶模态含能在84%以上。利用一阶模态及其特征值计算年径流量及年输沙量序列的变化特征可知，长江和鄱阳湖的径流变化率分别为2%～10%和10%～20%;而输沙变化率都接近90%。洞庭湖的水沙变化率要远远大于长江和鄱阳湖，其中荆江三口水沙变化率在140%以上。三峡工程对长江和荆江三口的减沙存在显著的影响。

参考文献：

[1] 彭俊.1950年以来鄱阳湖流域水沙变化规律及影响因素分析[J].长江流域资源与环境，2015，24（10）：1751-1761.

[2] 王海斌.长江流域干流水沙阶段性分析（1956～2009年）[J].地下水，2012（3）：150-152.

[3] 许全喜，童辉.近50年来长江水沙变化规律研究[J].水文，2012，32（5）：38-47，76.

[4] 张强，陈桂亚，姜彤，等.近40年来长江流域水沙变化趋势及可能影响因素探讨[J].长江流域资源与环境，2008（2）：257-263.

[5] 覃红燕.近50余年洞庭湖水文环境演变及其成因分析[D].长沙：湖南农业大学，2013.

[6] 李微，李昌彦，吴敦银，等.1956～2011年鄱阳湖水沙特征及其变化规律分析[J].长江流域资源与环境，2015，24（5）：832-838.

[7] 邴建平.长江—鄱阳湖江湖关系演变趋势与调控效应研究[D].武汉：武汉大学，2018.

[8] 魏丽，卢金友，刘长波.三峡水库蓄水后长江上游水沙变化分析[J].中国农村水利水电，2010（6）：1-4.

[9] 毛海涛，王正成，林荣，等.三峡水库蓄水后上下游河段水沙特性变化及影响因素分析[J].水资源与水工程学报，2019，30（5）：161-169.

[10] ZHANG Q，XU C，BECKER S，et al.Sediment and runoff changes in the Yangtze River basin during past 50 years[J].Journal of Hydrology，2006（331）：511-523.

[11] ZHAO G，MU X，SU B，et al.Analysis of streamflow and sediment flux changes in the Yangtze River basin[J].Water International，2012，37（5）：537-551.

[12] ZHANG Q，YE X，WERNER A D，et al.An investigation of enhanced recessions in Poyang Lake：Comparison of Yangtze River and local catchment impacts[J].Journal of Hydrology，2014（517）：425-434.

[13] DAI Z，LIU J T，XIANG Y.Human interference in the water discharge of the Changjiang（Yangtze River），China[J].Hydrological Sciences Journal，2015，60（10）：1770-1782.

[14] WANG G，JIANG T，BLENDER R，et al.Yangtze 1/f discharge variability and the interacting river-lake system[J].Journal of Hydrology，2008（351）：230-237.

[15] MOTE P W，PARSON E A，HAMLET A F，et al.Preparing for climatic change：the water，salmon，and forests of the Pacific Northwest[J].Climatic Change，2003（61）：45-88.

（編辑：江 文）

Analysis on runoff and sediment load in Changjiang River，Dongting Lake and

Poyang Lake based on modal decomposition method

HUANG Yinghua1，WAN Diwen2，DING Yutang3

（1.Jiepai Hydropower Hub Management Office of Jiangxi Provincial Port and Shipping Administration，Yingtan 335001，China; 2.Xiajiang Water Conservancy Project Management Bureau in Jiangxi Province，Nanchang 330046，China; 3.State Key Laboratory of Hydrology，Water Resources and Hydraulic Engineering，Nanjing Hydraulic Research Institute，Nanjing 210029，China）

Abstract：

In order to analyze trend of annual runoff and sediment load in the Changjiang River and its two river-connected lakes，Dongting Lake and Poyang Lake，we adopted modal decomposition method to extract variation characteristic of 19 major hydrological stations in the Changjiang River，Dongting Lake and Poyang Lake during 1956～2018.The results showed that predominant periods of 8 and 20 years were discovered in the runoff series.On the basis of the principal component（capturing more than roughly 90% total energy），total change ratios of runoff and sediment loads during the last 63 years were evaluated.For the runoff，absolute total change ratios were 2%～10% and 10%～20% for the Changjiang River and Poyang Lake respectively;while for the sediment load，they were around 90%.However，the change rate of flow and sediment in Dongting Lake were much larger，especially for the three outlets on the Jingjiang River reach（over 140%）.

Key words：

runoff;sediment load;modal decomposition;Changjiang River;Dongting Lake;Poyang Lake