余煌浩，李彬權(quán)

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

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

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

1 研究方法

1.1 趨勢分析及突變檢測

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

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

1.2 天然期徑流模擬方法

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

1.2.1降水蒸發(fā)因子法

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

Q=a×P+b×E+c

(1)

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

1.2.2改進(jìn)成因分割法

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

Q=K1×Pm1+K2×P+C1

(2)

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

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

(3)

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

1.2.3改進(jìn)趙文林法

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

QX=K4×PXm2+K5×QFpast

(4)

QF=K7×PFm3×HAm4

(5)

引入潛在蒸散發(fā)因子，則可改寫為

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——汛期潛在蒸散發(fā)量和非汛期潛在蒸散發(fā)量；QFpast——前一時段非汛期徑流深；HA——前一個汛期降水量與前一個汛期徑流深之差；K4、K5、K6、K7、K8、m2、m3、m4、C2、C3——天然期徑流系列的回歸參數(shù)。

1.2.4改進(jìn)張經(jīng)之法

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

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

(9)

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

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

(10)

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

最優(yōu)回歸參數(shù)根據(jù)非線性回歸擬合得到。輸入模型表達(dá)式和參數(shù)初始值，經(jīng)過多次迭代直到殘差平方和最小時停止迭代，得到最優(yōu)回歸參數(shù)。

1.3 徑流變化的歸因分析

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

ΔQT=QHR-QB

(11)

ΔQH=QHR-QHN

(12)

ΔQC=ΔQT-ΔQH

(13)

(14)

(15)

式中 ΔQT——徑流變化總量；ΔQH——人類活動對徑流的影響量；ΔQC——?dú)夂蜃兓瘜搅鞯挠绊懥浚籕B——天然期實測徑流量；QHR——人類活動影響期的實測徑流量；QHN——人類活動影響期模擬的徑流量；ηH、ηC——為人類活動和氣候變化對徑流影響的貢獻(xiàn)比例。

2 結(jié)果及分析

2.1 趨勢分析與突變檢測

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

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

2.2 天然期徑流系列擬合結(jié)果

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

2.3 徑流減少的歸因分析

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

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

禿尾河流域在影響期Ⅰ中徑流減少的主要驅(qū)動因素為氣候變化，表現(xiàn)為多年平均降水量(339 mm)遠(yuǎn)低于天然期的數(shù)值(420 mm)；但自1999年黃土高原退耕還林還草政策實施后，禿尾河流域植被覆蓋面積顯著上升，從而導(dǎo)致影響期Ⅱ中徑流減少的主要驅(qū)動因素轉(zhuǎn)變?yōu)槿祟惢顒覽14]。

表2 7種數(shù)據(jù)驅(qū)動方法的參數(shù)取值

表3 天然期年徑流擬合精度統(tǒng)計

注：表格中括號內(nèi)容為未改進(jìn)前擬合精度。

時期實測值/mm模擬值/mm減少比例/%人類活動的貢獻(xiàn)絕對值/mm比例/%氣候變化的貢獻(xiàn)絕對值/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 結(jié)論

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