基于數據驅動的陜北風沙區禿尾河流域徑流變化歸因

余煌浩，李彬權

(河海大學水文水資源學院，江蘇 南京 210098)

受氣候變化和高強度人類活動影響，中國北方河流江河徑流量總體上呈現顯著性減少趨勢，在黃河流域，人類活動是徑流變化的主要驅動因素[1]；同時在空間分布上，人類活動的影響貢獻自上游往下游逐漸增大[2]。在黃河中游黃土高原地區，大部分支流流域徑流均表現為銳減趨勢，人類活動影響為主導因素。鮑振鑫等[3]基于VIC模型對窟野河流域徑流進行歸因分析發現，影響期與天然期相比，人類活動對徑流減少的貢獻為61%～76%。Li等[4]運用數理統計方法對無定河和蘆河流域徑流歸因發現，人類活動對兩流域徑流減少的貢獻率分別為85%～90%和83%～86%。在渭河流域，Gao等[5]采用雙累積曲線法計算得到人類活動對徑流減少的貢獻率為83%。然而在下墊面植被條件較好的延河流域，Wu等[6]基于SWAT模型計算得到氣候變化的貢獻率為56%，為徑流減少的主要影響因素。薛帆等[7]基于Budyko假設和分形理論對北洛河流域徑流進行歸因識別，研究得出流域上游和下游徑流減少主要受人類活動影響，而中游主要受氣候變化影響。目前，徑流變化歸因分析基本是建立在水文模擬基礎上的，很大程度上受模擬精度的限制。然而，黃土高原地區水文模型精度普遍不高[8-9]，特別是風沙區等特殊地貌條件下的水文模型適用性不強。本文從數據驅動途徑改進風沙區典型流域(以陜北禿尾河為例)徑流模擬方法，選擇適用性較強的徑流模擬方法對徑流變化進行歸因分析。

禿尾河發源于陜西省神木縣瑤鎮西北的公泊海子，全長140 km，流域面積3 294 km2(圖1)，屬于干旱半干旱風沙區，多年平均面雨量為404 mm，雨量在年內分配不均，汛期(6—9月)占比達75%[10]。研究數據包括90 m分辨率的數字高程數據DEM(來源于地理空間數據云，https://www.gscloud.cn/)、2個氣象站的關鍵氣象要素數據(來源于中國氣象數據網，http://data.cma.cn/)以及9個雨量站的日降水和2個水文站的日徑流數據。日潛在蒸散發量根據Penman-Monteith方程[11]計算。各站點數據根據泰森多邊形法處理得到流域面平均系列。

1 研究方法

1.1 趨勢分析及突變檢測

采用Mann-Kendall檢驗法(M-K法)[12]對降水、潛在蒸散發和徑流系列進行趨勢分析。根據M-K檢驗統計值Z大于0(或小于0)，確定系列呈上升(或下降)趨勢，顯著性水平α=0.1。

突變檢測選擇滑動t檢驗、有序聚類法(OC法)、Pettitt法、標準正態均一性檢驗(SNHT法)這4種方法[4,12]，并采用水文變異綜合診斷的專家評分法[13]確定最終的突變點，具體步驟為：①每種方法評分總分為1，對處在可信度較低區間的點給予較低的評分，對在可信度較高區間的點給予較高的評分；②若檢驗方法對樣本檢驗結果所處區間無要求，則對檢驗的變異點給予相同的評分；③對沒有可信度的檢測結果評分為0，然后采取加權平均和歸一化處理計算各可能變異點的綜合權重，綜合權重最大的點即為最可能突變點。

1.2 天然期徑流模擬方法

采用降水蒸發因子法、成因分割法、趙文林法、張經之法這4種方法[4,10]對天然期徑流系列進行擬合，同時引入潛在蒸散發因子改進后3種方法，使其更為合理。

1.2.1降水蒸發因子法

直接建立徑流量與降水量、潛在蒸散發量之間的統計關系，見式(1)：

Q=a×P+b×E+c

(1)

式中Q、P、E——年徑流深、年降水量、年潛在蒸散發量；a、b、c——天然期徑流系列的回歸參數。

1.2.2改進成因分割法

按徑流的不同形成過程，成因分割法將徑流分解為地下徑流(基流)和地表徑流(表流)，分別建立降水與地下徑流、地表徑流的相關關系，見式(2)、(3)：

Q=K1×Pm1+K2×P+C1

(2)

在式(2)中加入潛在蒸散發因子，則可改寫為：

Q=K1×Pm1+K2×P+K3×E+C1

(3)

式中Q、P、E——年徑流深、年降水量、年潛在蒸散發量；K1、K2、K3、m1、C1——天然期徑流系列的回歸參數。

1.2.3改進趙文林法

趙文林法將徑流系列劃分為汛期徑流和非汛期徑流，分別建立汛期徑流與汛期降水、上一期非汛期徑流的相關關系，以及非汛期徑流與非汛期降水、前一個汛期雨量與前一個汛期徑流深之差的相關關系，見式(4)—(8)：

QX=K4×PXm2+K5×QFpast

(4)

QF=K7×PFm3×HAm4

(5)

引入潛在蒸散發因子，則可改寫為

QX=K4×PXm2+K5×QFpast+K6×EX+C2

(6)

QF=K7×PFm3×HAm4+K8×EF+C3

(7)

Q=QX+QF

(8)

式中 QX、QF——汛期徑流深和非汛期徑流深；PX、PF——汛期降水量和非汛期降水量；EX、EF——汛期潛在蒸散發量和非汛期潛在蒸散發量；QFpast——前一時段非汛期徑流深；HA——前一個汛期降水量與前一個汛期徑流深之差；K4、K5、K6、K7、K8、m2、m3、m4、C2、C3——天然期徑流系列的回歸參數。

1.2.4改進張經之法

張經之法根據汛期、非汛期降水量及汛期最大1日降水量，建立年徑流量，見式(9)、(10)：

Q=K9(PX×fm5+PFm6)+C4

(9)

引入潛在蒸散發因子，則可改寫為：

Q=K9(PX×fm5+PFm6)+K10×E+C4

(10)

式中Q、E——年徑流深和年潛在蒸散發量；PX、PF——汛期降水量和非汛期降水量；f——汛期最大一日降水量與汛期降水量的比值；K9、K10、m5、m6、C4——天然期徑流系列的回歸參數。

最優回歸參數根據非線性回歸擬合得到。輸入模型表達式和參數初始值，經過多次迭代直到殘差平方和最小時停止迭代，得到最優回歸參數。

1.3 徑流變化的歸因分析

對比天然期、人類活動影響期的實測和模擬徑流系列，建立徑流變化的歸因分析方法如下[3-5]：

ΔQT=QHR-QB

(11)

ΔQH=QHR-QHN

(12)

ΔQC=ΔQT-ΔQH

(13)

(14)

(15)

式中 ΔQT——徑流變化總量；ΔQH——人類活動對徑流的影響量；ΔQC——氣候變化對徑流的影響量；QB——天然期實測徑流量；QHR——人類活動影響期的實測徑流量；QHN——人類活動影響期模擬的徑流量；ηH、ηC——為人類活動和氣候變化對徑流影響的貢獻比例。

2 結果及分析

2.1 趨勢分析與突變檢測

1961—2015年流域年尺度降水和潛在蒸散發呈不顯著下降趨勢(統計特征量Z分別為-0.19和-0.16)，相應的年徑流深為顯著下降趨勢(Z=-8.23)，見圖2。4種方法檢測得到禿尾河高家川站1961—2015年年徑流系列的突變點主要有1979、1983、1996年(表1)，它們的綜合權重分別為0.575、0.250、0.175，據此可確定主突變點為1979年、次突變點為1996年(1983年與主突變點位置接近，暫不考慮)。將整個徑流系列劃分為天然期(1961—1979年)、影響期I(1980—1996年)和影響期II(1997—2015年)3個階段。支童等[14]在對禿尾河流域的降水、潛在蒸散發和徑流系列進行趨勢分析和突變檢測時，也得出同樣的結論。據此可判斷上述具有一定可靠性。

突變檢測方法可能突變年份評分滑動t檢驗19790.519960.5OC法19790.819960.2Pettitt法19831.0SNHT法19791.0

2.2 天然期徑流系列擬合結果

采用上述7種數據驅動的統計方法對天然期年徑流系列進行擬合，選用Nash-Sutcliffe效率系數(NSE)、Pearson相關系數(PCC)、均方根誤差(RMSE)和平均絕對誤差(MAE)進行擬合精度評定；利用非線性回歸擬合得到各方法的最優參數值見表2，擬合精度統計結果見表3。結果表明，改進后3種數據驅動方法的NSE、PCC、RMSE和MAE 4種精度指標較原方法均有大幅度提升；總體上，改進趙文林法的擬合精度最高(NSE=0.77、PCC=0.88、RMSE=9.10、MAE=7.08)，可用于后續徑流變化歸因研究。

2.3 徑流減少的歸因分析

以實測降水資料為數據驅動，利用改進趙文林法重構2個影響期的徑流系列，見圖3，計算得到氣候變化和人類活動對徑流減少的貢獻見表4。在影響期Ⅰ，年徑流模擬值略高于實測值，表明該時期人類活動的影響不顯著(對徑流減少的貢獻為19.1%)，而氣候變化影響占主導地位(貢獻率為80.9%)。相較于影響期I，影響期Ⅱ的多年平均年徑流量比天然期減少幅度更大(-47.1%)；年徑流模擬值遠高于實測值，表明該時期人類活動的影響較大，對徑流減少的貢獻提高至58.8%，占主導地位。

從1980—2015年逐年氣候變化和人類活動影響貢獻柱狀圖來看，影響期Ⅰ中多數年份氣候變化貢獻率都大于人類活動；影響期Ⅱ中多數年份的人類活動貢獻率均大于氣候變化的貢獻，且呈增長趨勢，表明人類活動對徑流減少的影響程度呈增長趨勢。對比天然期與整個影響期徑流變化，結果表明多年平均年徑流減少37%，氣候變化和人類活動影響的貢獻率分別為54.2%、45.8%。

禿尾河流域在影響期Ⅰ中徑流減少的主要驅動因素為氣候變化，表現為多年平均降水量(339 mm)遠低于天然期的數值(420 mm)；但自1999年黃土高原退耕還林還草政策實施后，禿尾河流域植被覆蓋面積顯著上升，從而導致影響期Ⅱ中徑流減少的主要驅動因素轉變為人類活動[14]。

表2 7種數據驅動方法的參數取值

表3 天然期年徑流擬合精度統計

注：表格中括號內容為未改進前擬合精度。

時期實測值/mm模擬值/mm減少比例/%人類活動的貢獻絕對值/mm比例/%氣候變化的貢獻絕對值/mm比例/%天然期125.10124.80影響期Ⅰ92.8899.0525.86.1719.126.0580.9影響期Ⅱ66.17100.8447.134.6758.824.2641.2整個影響期78.78100.0037.021.2245.825.1054.2

3 結論

利用M-K趨勢分析方法發現，禿尾河流域1961—2015年年尺度降水、潛在蒸散發系列呈不顯著下降趨勢，而流域出口水文站高家川站年徑流量呈顯著減少趨勢；利用水文變異綜合診斷方法得到年徑流系列的突變點為1979年和1996年。引入潛在蒸散發因子，改進成因分割法、趙文林法和張經之法，增強了數據驅動方法在陜北風沙區禿尾河流域徑流模擬中的適用性，其中改進趙文林法精度最高(NSE=0.77)。影響期(1980—2015年)與天然期(1961—1979年)相比，多年平均年徑流量減少比例為37%，其中氣候變化的影響貢獻占主導地位(54.2%)；但在影響期II(1997—2015年)，人類活動的影響更大，對徑流減少的貢獻比例為58.8%。