基于數(shù)據(jù)驅(qū)動(dòng)的陜北風(fēng)沙區(qū)禿尾河流域徑流變化歸因

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

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

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

禿尾河發(fā)源于陜西省神木縣瑤鎮(zhèn)西北的公泊海子，全長(zhǎng)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(來(lái)源于地理空間數(shù)據(jù)云，https://www.gscloud.cn/)、2個(gè)氣象站的關(guān)鍵氣象要素?cái)?shù)據(jù)(來(lái)源于中國(guó)氣象數(shù)據(jù)網(wǎng)，http://data.cma.cn/)以及9個(gè)雨量站的日降水和2個(gè)水文站的日徑流數(shù)據(jù)。日潛在蒸散發(fā)量根據(jù)Penman-Monteith方程[11]計(jì)算。各站點(diǎn)數(shù)據(jù)根據(jù)泰森多邊形法處理得到流域面平均系列。

1 研究方法

1.1 趨勢(shì)分析及突變檢測(cè)

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

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

1.2 天然期徑流模擬方法

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

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

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

Q=a×P+b×E+c

(1)

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

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

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

Q=K1×Pm1+K2×P+C1

(2)

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

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)系，以及非汛期徑流與非汛期降水、前一個(gè)汛期雨量與前一個(gè)汛期徑流深之差的相關(guān)關(guān)系，見(jiàn)式(4)—(8)：

QX=K4×PXm2+K5×QFpast

(4)

QF=K7×PFm3×HAm4

(5)

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

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

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

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

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

(9)

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

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)過(guò)多次迭代直到殘差平方和最小時(shí)停止迭代，得到最優(yōu)回歸參數(shù)。

1.3 徑流變化的歸因分析

對(duì)比天然期、人類(lèi)活動(dòng)影響期的實(shí)測(cè)和模擬徑流系列，建立徑流變化的歸因分析方法如下[3-5]：

ΔQT=QHR-QB

(11)

ΔQH=QHR-QHN

(12)

ΔQC=ΔQT-ΔQH

(13)

(14)

(15)

式中 ΔQT——徑流變化總量；ΔQH——人類(lèi)活動(dòng)對(duì)徑流的影響量；ΔQC——?dú)夂蜃兓瘜?duì)徑流的影響量；QB——天然期實(shí)測(cè)徑流量；QHR——人類(lèi)活動(dòng)影響期的實(shí)測(cè)徑流量；QHN——人類(lèi)活動(dòng)影響期模擬的徑流量；ηH、ηC——為人類(lèi)活動(dòng)和氣候變化對(duì)徑流影響的貢獻(xiàn)比例。

2 結(jié)果及分析

2.1 趨勢(shì)分析與突變檢測(cè)

1961—2015年流域年尺度降水和潛在蒸散發(fā)呈不顯著下降趨勢(shì)(統(tǒng)計(jì)特征量Z分別為-0.19和-0.16)，相應(yīng)的年徑流深為顯著下降趨勢(shì)(Z=-8.23)，見(jiàn)圖2。4種方法檢測(cè)得到禿尾河高家川站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)位置接近，暫不考慮)。將整個(gè)徑流系列劃分為天然期(1961—1979年)、影響期I(1980—1996年)和影響期II(1997—2015年)3個(gè)階段。支童等[14]在對(duì)禿尾河流域的降水、潛在蒸散發(fā)和徑流系列進(jìn)行趨勢(shì)分析和突變檢測(cè)時(shí)，也得出同樣的結(jié)論。據(jù)此可判斷上述具有一定可靠性。

突變檢測(cè)方法可能突變年份評(píng)分滑動(dòng)t檢驗(yàn)19790.519960.5OC法19790.819960.2Pettitt法19831.0SNHT法19791.0

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

采用上述7種數(shù)據(jù)驅(qū)動(dòng)的統(tǒng)計(jì)方法對(duì)天然期年徑流系列進(jìn)行擬合，選用Nash-Sutcliffe效率系數(shù)(NSE)、Pearson相關(guān)系數(shù)(PCC)、均方根誤差(RMSE)和平均絕對(duì)誤差(MAE)進(jìn)行擬合精度評(píng)定；利用非線性回歸擬合得到各方法的最優(yōu)參數(shù)值見(jiàn)表2，擬合精度統(tǒng)計(jì)結(jié)果見(jiàn)表3。結(jié)果表明，改進(jìn)后3種數(shù)據(jù)驅(qū)動(dòng)方法的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í)測(cè)降水資料為數(shù)據(jù)驅(qū)動(dòng)，利用改進(jìn)趙文林法重構(gòu)2個(gè)影響期的徑流系列，見(jiàn)圖3，計(jì)算得到氣候變化和人類(lèi)活動(dòng)對(duì)徑流減少的貢獻(xiàn)見(jiàn)表4。在影響期Ⅰ，年徑流模擬值略高于實(shí)測(cè)值，表明該時(shí)期人類(lèi)活動(dòng)的影響不顯著(對(duì)徑流減少的貢獻(xiàn)為19.1%)，而氣候變化影響占主導(dǎo)地位(貢獻(xiàn)率為80.9%)。相較于影響期I，影響期Ⅱ的多年平均年徑流量比天然期減少幅度更大(-47.1%)；年徑流模擬值遠(yuǎn)高于實(shí)測(cè)值，表明該時(shí)期人類(lèi)活動(dòng)的影響較大，對(duì)徑流減少的貢獻(xiàn)提高至58.8%，占主導(dǎo)地位。

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

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

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

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

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

時(shí)期實(shí)測(cè)值/mm模擬值/mm減少比例/%人類(lèi)活動(dòng)的貢獻(xiàn)絕對(duì)值/mm比例/%氣候變化的貢獻(xiàn)絕對(duì)值/mm比例/%天然期125.10124.80影響期Ⅰ92.8899.0525.86.1719.126.0580.9影響期Ⅱ66.17100.8447.134.6758.824.2641.2整個(gè)影響期78.78100.0037.021.2245.825.1054.2

3 結(jié)論

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