基于模態分解法的長江與兩湖水沙特征分析

黃瑩華 萬迪文 丁玉堂

摘要：為研究長江和洞庭湖、鄱陽湖的水沙變化，利用模態分解方法（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