漢江流域安康段降雨徑流特征分析及預測

劉易文，李家科*，丁 強，郝改瑞

(1.西安理工大學省部共建西北旱區生態水利國家重點實驗室，陜西 西安 710048；2.陜西省環境監測中心站，陜西 西安 710054)

降雨被認為是水文循環的主要組成部分，是河川徑流的主要來源，徑流對水資源規劃、管理、開發、模型構建非常重要[1-3]。而識別降雨徑流特征變化以及精確模擬河川徑流是水資源開發利用、水文模型構建的基礎工作之一。由于降雨徑流受氣候、人為等諸多不確定因素的影響，常常伴隨著非線性變化和較強的隨機性，屬于高度的非線性系統[4-6]。傳統的徑流預測方法只能反映出線性時間序列和簡單的非線性時間序列，而實際徑流往往受多因素影響，徑流序列與各影響因素之間的非線性特征較強。因此，對復雜系統而言，僅憑傳統單一模型，難以精確模擬整個徑流序列。基于深度學習和機器學習的組合模型在很大程度上緩解了徑流模擬精度低的問題，成為現階段廣受關注的熱點之一[7-8]。組合模型更能體現時間序列的非線性特征，而且具有很強的非線性映射能力，適用于隨機性較強的水文過程模擬[9]。此外，基于降雨徑流演變規律和趨勢特征，組合模型模擬的徑流序列更為準確、可靠[10]。

近年來，隨著深度學習和機器學習的推廣，眾多學者利用人工神經網絡所衍生出來的組合模型進行水文序列模擬。王少麗等[11]利用BP神經網絡對年徑流系數進行了模擬，并與線性主成分回歸模擬進行比對，研究發現BP神經網絡的模擬結果更好，精度更高。梁川等[12]利用偏最小二乘回歸及神經網絡對高寒區域河川徑流進行預測，結果表明偏最小二乘回歸的模擬結果合理，但神經網絡的模擬精度優于偏最小二乘回歸。李秀峰等[13]將徑流序列進行小波降噪處理后，利用偏最小二乘回歸對徑流進行了預測，結果表明該方法精度高且穩定性較好，可作為徑流預測的有效方法。

總體看，徑流預測的研究熱點仍為基于深度學習、機器學習的組合模型。但少有在突變成分附近或者在徑流局部對模擬效果進行優化的組合模型。基于此，本文依據漢江流域安康水文站控制斷面以上1956—2017年的實測降雨、氣溫及徑流數據，采用Mann-Kendall突變檢驗法識別降雨徑流突變成分，并結合R/S分析法對降雨徑流進行趨勢分析，利用Mrolet小波提取降雨徑流變化主周期。選用能較好識別時間序列的偏最小二乘回歸(PLSR)及BP神經網絡——偏最小二乘回歸(BP-PLSR)對徑流進行模擬，根據已識別的突變成分，對比2種方法在峰值處和突變成分附近的模擬精度，以期為流域水資源管理、水資源利用提供依據以及支撐。

1 研究區概況

漢江屬于長江一級支流，全長1 532 km，發源地位于陜西省寧強縣冢山，是南水北調工程的重要水源。漢江流域涉及湖北、陜西、甘肅、四川、重慶5個省(市)，面積達到15.9萬km2。漢江流域整體地勢西北高，東南低，屬于亞熱帶季風氣候區，平均年降雨量達873 mm，降雨充沛，汛期集中在5—10月，汛期徑流量可占全年75%，年際變化大。安康水文站(109°E，33.67°N)位于陜西省安康市漢濱區(圖1)。

圖1 漢江流域安康段

2 材料和方法

2.1 數據來源

本文主要收集了漢江流域安康水文站控制斷面以上1956—2017年多年實測降雨、氣溫及徑流數據資料。其中降雨、氣溫數據來自于中國科學院資源環境科學數據中心(http://www.resdc.cn/)，徑流數據來自于長江流域水文年鑒。

2.2 突變識別及趨勢分析

2.2.1Mann-Kendall突變檢驗法

M-K突變檢驗是一種非參數方法，不需要測試樣本遵循特定的分布[14]。同時，也不會受到少數樣本數值的影響，適用于多種類型且計算較為簡單[15]。

2.2.2R/S分析法

在M-K突變檢驗的基礎上，采用R/S法分析漢江流域安康水文站控制斷面以上降雨徑流的趨勢特征。R/S分析法的優勢在于可以定量描述時間序列的趨勢性，該趨勢特征由Hurst指數表征，具有良好的穩定性[16-17]。

2.3 降雨徑流周期分析

Mrolet小波是基于傅里葉變換的基小波函數與高斯函數推導而來，很好地克服了傳統譜分析方法的缺點，在實際應用中比實數形式小波更具優勢。Mrolet小波能夠提高降雨徑流序列周期變化檢驗水平，如今被廣泛運用于信號處理、模式識別、數據分析等領域。

2.4 徑流預測

2.4.1偏最小二乘回歸(PLSR)

PLSR有機結合了多種主流分析方法，包括多元線性回歸分析、主成分分析、相關分析，該方法優于最小二乘回歸[18]，且適用于時間序列的模擬。首先構建自變量和因變量的數據表，分別為X和Y；對二者實施線性組合，得到PLSR第一成分t1與u1，進行X對t1的回歸以及Y對u1的回歸，當第一成分對應的回歸方程滿足精度要求時終止計算；若誤差較大，利用自變量數據表與第一成分回歸后的殘余信息進行第二次成分提取，重復上述步驟至誤差最小[19-21]。并根據所提取的成分信息，形成回歸方程。

2.4.2BP神經網絡-偏最小二乘回歸(BP-PLSR)

BP神經網絡是多層前饋神經網絡，通常由輸入、隱藏、輸出層構成基礎網絡結構。該方法的核心在于反向傳遞(Back Propagation)，即依靠誤差反向傳播不斷優化神經網絡的權重和閾值，以期減小誤差。尤其對于非線性較強的徑流序列，BP神經網絡能更好地識別局部誤差，避免結果過擬合。BP神經網絡構建如下。

a)對數據進行歸一化處理：

(1)

式中xmin、xmax——原始數據的最小值及最大值。

b)本例中選用了3層神經網絡結構，選用S型傳遞函數(Log-Sigmoid)作為激活函數，激活正向傳遞過程，通過反傳誤差函數不斷調節網絡的權值以及閾值，使誤差函數E最小。

(2)

yi=f(Sj)

(3)

(4)

(5)

(6)

式中Sj——正向傳遞函數；Wi,j——結點i與節點j之間的權值；Wj——第j個節點輸出值；f——節點激活函數；dj——第j個輸出結果；b——閾值；N——隱藏層節點數；m——輸入層節點數；n——輸出層節點數；a——常數(范圍取1～10)。

c)PLSR能較好地表征時間序列，但該方法在徑流預測時容易陷入局部最優，造成序列過擬合。考慮到BP神經網絡Levenberg-Marquardt算法(L-M算法)具有自適應性噪聲消除功能，借此可以較好地消除目標序列噪聲，防止PLSR的輸出結果過擬合。借助Matlab中的Neural Net工具箱，將降雨、徑流及氣溫值作為輸入序列，選擇Levenberg-Marquardt算法對原始數據進行處理，此時輸出值為降噪后的徑流序列。處理后的徑流序列并不等同于徑流模擬，將其標準化后作為PLSR的目標序列運行模型。標準化公式如下：

(7)

(8)

式中F0——Y的標準化矩陣；E0——X的標準化矩陣；E(y)——Y的均值；E(xi)——X的均值；Sy、Sxi——Y、X的均方差；n——樣本數。

3 結果與討論

3.1 降雨徑流過程分析

漢江流域安康水文站控制斷面以上1956—2017年多年平均降雨為815.61 mm，流量為587.25 m3/s；其中年均最大降雨量出現在2010年(1 231.9 mm)，年均最大徑流量出現在1983年(1 294.42 m3/s)，年均最小降雨量及徑流量均出現在1999年，分別為525.8 mm、247.75 m3/s，二者極差分別為706.1 mm、1 046.67 m3/s。多年降雨、徑流過程見圖2。

根據圖2a，1956—2017年降雨量出現了4個變化階段。經歷了增加—減小—增加—減小的變化過程。其中，1974—1983年及1999—2001年降雨的上升趨勢最為明顯，2個階段平均值分別為898.9、844.8 mm。通過62年整體變化趨勢，判斷降雨呈現出整體增長趨勢。圖2b中，62年徑流變化同樣存在4個階段，與降雨序列一致，呈現出相同的變化過程。90年代前，平均徑流量大于62年整體均值的年份明顯多于90年代后。說明90年代后，流量整體減少趨勢明顯。2011—2016年，累計減少255.3 m3/s。徑流序列整體呈現下降趨勢，降雨徑流變化過程見表1。

表1 降雨徑流變化特征

3.2 降雨徑流突變成分識別

漢江流域安康水文站控制斷面以上降雨徑流突變識別見圖3。結果顯示：M-K檢驗結果中降雨突變成分過多，存在偽突變點，結合累計距平法剔除偽突變點(圖4)。對比2種方法分析結果，得到降雨突變點出現在1973年、1984年、2002年。1973年為先減小后增加的轉折點，從1974年開始，降雨呈現持續增加現象，直到1984年，降雨增加17.4%。1984年突變前后降雨相對變化為30%。隨后降雨開始突變減小，直到第三次突變，降雨減小23.2%。

徑流序列突變點出現在1977年、1979年、1985年。分析發現1977—1979年徑流量呈現小幅增加趨勢，結合累計距平檢驗結果，認為1979年為偽突變點，固將其剔除。突變成分主要體現為先增加后減小，直到1985年，徑流突變減小37.9%。

結合了M-K突變檢驗與累計距平檢驗，發現累計距平檢驗結果與M-K檢驗結果呈包含關系，結合分析降雨徑流數據，對偽突變點進行剔除。根據M-K突變檢驗，發現降雨徑流2個序列的統計量UF與UB的交點均在顯著性水平線之間，所以突變檢驗通過顯著性分析。

降雨徑流突變檢驗結果與安康市近50年降雨量變化趨勢相吻合[22-23]。安康市60年代初與80年代初降雨量為極大時期，同時60年代至80年代為偏豐時期，而在80年代以后降雨偏少。其中1974年突變點經歷了降水的先減少后增大，從1984年突變點開始，降水從極大值突變減小。經檢驗，徑流突變點滯后于降雨突變，對比其他學者的徑流突變點分析結果[24]，由于徑流突變檢驗的時間段選擇與方法不同，安康徑流突變僅出現在1984年，與本文結果相差1年，且徑流突變均滯后于降雨突變，分析結果相似。

3.3 趨勢分析

選用R/S分析法，設置lgm為橫坐標，lg(R/S) 為縱坐標，并結合最小二乘法繪制散點圖。根據散點擬合直線，直線斜率即為Hurst指數。Hurst指數不同，對應的趨勢變化則不同，其不同取值對應的趨勢見表2。R/S計算結果見圖5。

表2 Hurst指數取值

a)降雨

a)降雨

a)降雨

a)降雨

a)降雨

b)徑流

對R/S結果進行分析可得：降雨Hurst指數接近0，說明與未來趨勢相反。降雨突變減小的趨勢減弱，降雨序列在研究時段內將呈現增長趨勢；徑流序列Hurst指數為0.025 1，H指數介于0～0.5之間，說明該序列未來趨勢將與過去相反，且H指數接近0，表明這種反持續狀態較為強烈(表3)。根據徑流量突變分析可以看出徑流整體呈現減少趨勢，且第二突變成分之后再未出現徑流量跳躍增加或跳躍減小，同時驗證了這種突變趨勢正在減弱，該結果與張玨等[16]的分析結果一致。

表3 Hurst計算結果

3.4 小波分析

利用Mrolet小波對降雨徑流進行周期分析，獲得降雨徑流的周期性變化特征(圖6)。

b)徑流

1956—2017年，可看出降雨量有6、10 a 2個周期。在10 a以下降雨量周期變化規律較為明顯，2個周期的波動控制著安康降雨在整個時間域內的變化特征，震蕩周期附近都存在著明顯的豐、枯交替變換，其中偏豐期、偏枯期交替突變的點出現在1958年、2013年。對于徑流周期變化，62年內安康徑流量存在周期14 a，此處震蕩最強，為主周期。5 a處震蕩較弱，為副周期。16 a以下徑流周期變化規律較為明顯，其中偏豐期、偏枯期交替突變的點出現在1979年、1983年。

3.5 徑流預測

首先對氣溫序列、降雨及徑流序列進行相關性檢驗，分析得到徑流與氣溫以及降雨與徑流之間均存在相關性(P<0.005)，通過了顯著性檢驗(表4)，本例中所選用的自變量和因變量之間存在較為明顯的相關性。自變量之間存在相關性會造成共線性，進而影響回歸結果，而偏最小二乘回歸會減小自變量共線性對結果造成的影響。

表4 相關性檢驗(新增降雨及氣溫之間的相關性檢驗)

基于1956—2017年62 a的實測年降雨及氣溫作為基礎輸入數據，分別利用PLSR及BP-PLSR 2種方法對徑流量進行預測。將1956—2004年設置為訓練期，2005—2017年設置為驗證期。預測精度指標選擇均方根誤差(Root-Mean-Square Error,RMSE)及納什效率系數(Nash-Sutcliffe Efficiency Coefficient ,NSE)，精度指標計算公式如下：

(9)

(10)

利用PLSR對徑流進行模擬。首先根據平均值與均方差對原始數據進行標準化處理，提取出2個成分，得到PLSR回歸方程：

y*=1207.9+0.9x1-83x2

(11)

式中y*——預測年徑流量；x1——實測年降雨量；x2——實測年氣溫值。

PLSR徑流模擬結果見圖7。在模型訓練期，RMSE為153.093，NSE為0.489；驗證期內，RMSE為70.275，NSE為0.259；在整個模擬時期內，PLSR模擬結果精度較低，RMSE為152.182，NSE為0.456，納什效率系數(NSE)低于0.5且R2為0.500 3。模擬精度低的主要原因為：62 a的時間周期內，徑流量出現突變，在徑流峰值前后出現了過擬合現象，且在突變點附近出現局部最優問題。1983年為徑流最大年份，徑流峰值處模擬值為944.58 m3/s，與實測值的誤差為349.84 m3/s；在突變點附近，1977、1979、1985年的模擬值分別為678.30、781.00、635.19 m3/s，與實測值的誤差分別為288.6、297.8、254.9 m3/s。

a)不同年份下徑流量實測、模擬值

b)訓練期實測-模擬

c)驗證期實測-模擬

d)實測-模擬

采用BP-PLSR對徑流進行模擬，模擬結果見圖8。所生成的回歸方程為：

y*=1333.7+0.6x1-78.7x2

(12)

式中y*——預測年徑流量；x1——實測年降雨量；x2——實測年氣溫值。

在訓練期RMSE為24.371，相比PLSR在驗證期的結果減少了84.4%。NSE為0.784，明顯大于PLSR在訓練期的效率系數；在驗證期，RMSE為12.343，比PLSR在驗證期的結果減少了82.5%。NSE為0.891，BP-PLSR模擬結果精度明顯優于PLSR。在整個模擬時期，RMSE為92.863，誤差下降了40%。NSE明顯增加，達到0.797，模擬結果良好。1983年峰值處，實測值為1 294.42 m3/s，模擬值為1 114.60 m3/s。相比PLSR，誤差下降了48.6%。在突變點附近，1977年、1979年、1985年模擬值分別為613.94、498.77、593.92 m3/s。相比PLSR，誤差分別降低了22.30%、94.76%、91.40%。

未對原始數據進行預處理，用PLSR進行模擬后，發現模擬結果精度較低，特別是在突變成分附近及峰值處的模擬結果精度低。對原始數據進行BP神經網絡預處理，可以解決數據在峰值處的過擬合問題，能較好地避免模型在突變成分附近陷入局部最優。BP-PLSR模擬結果與實測值最為接近，模擬結果可信度高。

BP-PLSR模擬效果優于傳統PLSR，模擬精度有所提高。但在徑流局部區域難免出現過擬合現象，只有在1977年處，2種方法的模擬結果相對誤差均大于50%(表5、圖9)。在訓練期，BP-PLSR對過擬合及局部最優問題控制較好，未出現明顯泛化誤差(表6)。由于PLSR容易在徑流峰值處及突變點附近出現明顯泛化誤差，而BP-PLSR可以優化上述問題，能明顯提高模擬精度。

表5 突變成分附近及峰值處的對比

a)不同年份下的徑流量實測、模擬值

b)訓練期實測-模擬

c)驗證期實測-模擬

d)實測-模擬

圖9 局部模擬結果

表6 不同模擬方法對應的徑流模擬精度指標

將突變成分附近及峰值附近實測徑流序列作為目標序列，選擇對應時段降雨及氣溫值進行模擬，將結果回代到全序列以判斷BP-PLSR的適用性(圖10、11)，得到的回歸方程為：

y*=1452.3+0.9x1-101.5x2

(13)

所選時段為1964—2002年，NSE系數為0.776，決定系數R2為0.809 4。

在突變點1977年、1979年、1984年處，模擬值與實測值誤差為-224.185、-15.570、-21.310 m3/s，1977年模擬值誤差較大。峰值處模擬與實測值誤差為179.77 m3/s。整體模擬效果較好。將結果回代到全序列中時，得到NSE值為0.798，決定系數R2為0.829 5。整體模擬效果理想，但是仍然存在個別突變點及峰值處模擬誤差較大的現象。通過2種模擬方式，即利用全序列體現局部效果及依靠局部回代全序列，都體現了較好的模擬結果，NSE系數均接近0.8，R2超過0.8，模擬結果良好，說明了該組合模擬適用性較好。

a)1956—2017年徑流量實測、模擬值

b)實測-模擬

a)1964—2002年徑流量實測、模擬值

b)實測-模擬

4 結論

a)62 a漢江流域安康水文站控制斷面以上多年平均降雨為815.61 mm，1974—1983年、1999—2001年降雨上升趨勢最為明顯，2個階段平均值分別為898.9、844.8 mm，降雨呈現出整體增長趨勢。62 a徑流均值為587.25 m3/s ，90年代前有19 a的徑流量大于整體平均值587.25 m3/s，而90年代后只有5 a的徑流大于587.25 m3/s。徑流序列多個階段均存在減少，徑流整體呈現下降趨勢。

b)降雨突變成分出現在1973年、1984年、2002年，徑流序列突變成分出現在1977年、1985年。選擇顯著性水平為5%，降雨徑流突變識別通過顯著性檢驗。

c)利用R/S分析，降雨徑流序列Hurst指數均介于0～0.5之間，說明未來將存在相反趨勢；H指數接近0，表明這種反持續狀態較為強烈。降雨徑流趨勢分析相關系數分別為0.941 5、0.894 3，趨勢分析效果良好。

d)PLSR的徑流模擬結果精度較低，RMSE為152.182、NSE為0.456。NSE值低于0.5，誤差較大，模擬結果受到了過擬合及局部最優影響。對原始數據進行BP處理后，能有效提高模擬精度，RMSE為92.863、NSE為0.797；RMSE降低40%且NSE大于0.5。在峰值處，模擬精度提升48.6%；在突變點附近，精度提高了69.5%。模擬結果可信度增加，有效解決和避免了過擬合及局部最優問題。

5 結語

基于漢江流域安康水文站1956—2017年62 a實測降雨、氣溫及徑流數據，利用M-K突變檢驗進行突變點識別，選用R/S分析法以及小波分析對安康水文站控制斷面以上降雨徑流序列進行趨勢及周期分析。采用PLSR、BP-PLSR對徑流進行模擬，并對兩者進行了對比。由于水文序列隨機性強，不確定因素較多，受氣候、下墊面以及人類活動的影響大，傳統方法難以較好模擬非線性序列，建議充分利用深度學習和機器學習，優化現有模擬方法，提高模擬精度。