康靖羚, 丁文峰,2, 韓昊宇, 郭宜薇, 王一然

(1.長江水利委員會長江科學(xué)院, 武漢 430010; 2.水利部 山洪地質(zhì)災(zāi)害防治工程技術(shù)研究中心, 武漢 430010)

黃河是我國的第二大河，華夏民族在此繁衍生息，孕育出了古老燦爛的中華文明。長期以來，強烈的人為活動以及全球氣候變化等因素導(dǎo)致出現(xiàn)水少沙多、水沙異源等問題，水沙關(guān)系極不協(xié)調(diào)[1]，黃土高原已成為我國土壤侵蝕最嚴(yán)重的地區(qū)[2]。鑒于水土流失的嚴(yán)重性，自20世紀(jì)50年代以來，我國政府將黃河流域治理作為全國水土保持工作的重點，開展了一系列生態(tài)治理工程，合理開發(fā)利用土地資源，有效控制水土流失。20世紀(jì)末期，經(jīng)過治理的黃土高原小流域示范區(qū)水土流失治理度由46.1%增加至80.2%，綜合治理成效顯著[3]。2000年，全國陸續(xù)啟動了“退耕還林”等生態(tài)工程，植被覆蓋率迅速增加[4]，截至2015年，黃河中游的水土流失治理度已達60%[5]。潼關(guān)控制站實測年均輸沙量也從20世紀(jì)前半葉的15.92億t減少到2000—2018年的2.4億t，其中2015年實測沙量僅有0.55億t[1,6]。

河流輸沙量作為水土流失的集中反映，已成為全球水沙關(guān)系研究的重要內(nèi)容，也是間接研究流域土壤侵蝕產(chǎn)沙的主要手段。由于大理河地處黃河中游粗泥沙集中來源區(qū)[7]，作為大理河的一級支流，岔巴溝流域成為研究黃土丘陵溝壑區(qū)土壤侵蝕以及水沙關(guān)系的理想?yún)^(qū)域[8]，此地區(qū)長期的觀測數(shù)據(jù)也為研究提供了良好基礎(chǔ)，了解岔巴溝流域水沙關(guān)系及其變化趨勢和影響因素對治理黃土丘陵溝壑區(qū)乃至黃土高原的水土流失具有重要意義。對岔巴溝流域的相關(guān)研究從20世紀(jì)80年代就已經(jīng)開始，目前為止許多學(xué)者在此流域?qū)斏沉孔兓归_研究，并在流域面積與水沙關(guān)系的空間變異性[9-11]、不同空間尺度下的泥沙量或泥沙濃度[12-14]、各類流域產(chǎn)沙模型[15-19]以及泥沙動力學(xué)等[20-21]方面取得許多成果。但大多數(shù)研究側(cè)重于不同空間尺度下的水沙關(guān)系及模型建立，對不同時間尺度下的流域水沙關(guān)系研究較少。本文通過分析岔巴溝流域1970—2016年的徑流泥沙變化規(guī)律，探討年時間尺度以及次洪水時間尺度下流域水沙關(guān)系變化趨勢和影響因素，以期為流域水沙模型建立提供依據(jù)。

1 研究區(qū)概況

岔巴溝流域位于陜西省子洲縣，屬于黃土高原丘陵溝壑區(qū)第一副區(qū)，地處東經(jīng)109°47′，北緯37°31′。流域面積205 km2，溝壑密度1.05 km/km2，主溝比降7.57‰，海拔高度900～1 100 m[22]。屬半干旱區(qū)，為溫帶大陸性季風(fēng)氣候，干燥少雨，降雨年內(nèi)分配不均，大部分集中于6—9月，且多是高強度短歷時暴雨[23]。

岔巴溝流域土壤由黃土發(fā)育而成，結(jié)構(gòu)疏松，易受侵蝕[21]。早在1959年，岔巴溝流域的水土流失治理工作就已開始，70年代后，流域內(nèi)大量水庫、淤地壩建成，并逐漸發(fā)揮效益，至1977年末，所修建庫壩數(shù)量達444座，總庫容2 651萬m3。至1980年，岔巴溝流域治理面積已達30.7 km2，占總面積的15%[23]。截至1993年，水土保持治理程度已達到70%以上[19]。90年代后期，岔巴溝流域響應(yīng)國家號召，開展了大規(guī)模的水土保持工作，生態(tài)環(huán)境明顯好轉(zhuǎn)[24]。

2 數(shù)據(jù)來源及處理

2.1 數(shù)據(jù)來源

本文所用數(shù)據(jù)資料出自《黃河流域水文年鑒》，包括雨量站降水?dāng)?shù)據(jù)以及水文站徑流和泥沙數(shù)據(jù)(1970—2016年)。本次研究選擇4個測站，分別為和民墕站、劉家坬站、姬家鹼站和曹坪站，具體信息見表1。

表1 測站信息

2.2 數(shù)據(jù)處理

2.2.1 次洪水?dāng)?shù)據(jù)處理 一次完整的洪水事件是指從徑流量增加至基流以上，到僅基流構(gòu)成徑流量為止[25]，故采用數(shù)字濾波法[26-27]去除基流，得到新的洪峰流量和徑流量。研究變量包括次洪水輸沙量(ESY)，采用式(1)計算；洪峰流量(Qw，m3/s)，徑流過程的最大值；徑流深(H，mm)，即徑流量與流域面積之比；次降雨量(P，mm)，即一次降雨總量；降雨歷時(D，h)，次降雨持續(xù)時長；平均降雨強度(I，mm/h)，即降雨量與降雨歷時之比以及徑流系數(shù)(Rc)，即徑流深和降雨量之比。

ESY=∑SqΔt

(1)

式中:ESY為次洪水輸沙量(萬t);Δt為水文站水文觀測的時間間隔(s);q為Δt時間間隔內(nèi)的流量(m3/s);S為Δt時間間隔內(nèi)的徑流含沙量(t/m3)。

從所有篩選出的數(shù)據(jù)中，選取洪峰流量大于1 m3/s、有對應(yīng)降雨記錄的洪水過程進行分析，共計106場洪水事件。

2.2.2 年尺度洪水?dāng)?shù)據(jù)處理 為研究年時間尺度下降雨、徑流和泥沙變量之間的關(guān)系，采用年平均輸沙量(ESYa，萬t)、年最大洪峰流量(Qwa，m3/s)，即一年中全部洪水事件中的最大洪峰流量值、年平均徑流深(Ha，mm)、年平均降雨量(Pa，mm)、年平均徑流系數(shù)(Rca)以及年最大24 h雨強(I24 h，mm/h)等研究變量進行分析。采用Mann-Kendall檢驗法[28-30]對年尺度研究變量進行趨勢分析和突變點檢測。數(shù)據(jù)處理均采用Excel 2016以及SPSS 25，圖像繪制采用ArcGIS 10.2以及Origin 2019。

3 結(jié)果與分析

3.1 年尺度研究變量分析

表2顯示了岔巴溝流域降雨量、徑流量和輸沙量的年際變化，可以看出降雨量呈增加趨勢，徑流量與輸沙量呈減小趨勢，降雨量、徑流量和輸沙量的最大值分別大于平均值13%，23%和117%，最小值分別比平均值小13%，20%和80%，年際變化較大。

表2 岔巴溝流域年降雨量、徑流量和輸沙量統(tǒng)計

3.1.1 年降雨量趨勢分析 根據(jù)實測降雨資料，在ArcGIS 10.2中采用泰森多邊形法計算流域面降雨量，并繪制年降雨量變化過程線，見圖1。M-K檢驗曲線見圖2。

圖1 年降雨量變化過程

圖2 年降雨量M-K檢驗曲線

統(tǒng)計量Z=2.75，|Z|>2.32，表明年降雨量在99%置信水平上呈顯著增加趨勢。顯著水平臨界線內(nèi)出現(xiàn)的交點分別位于1972—1973年與1986—1987年，但1972—1973年無前期研究數(shù)據(jù)，無法確定在此時間段內(nèi)發(fā)生突變，所以將第2個交點作為突變點，即降雨量于1986年出現(xiàn)突變。

3.1.2 年徑流量變化分析 圖3即流域年徑流量變化過程線，M-K檢驗曲線見圖4，統(tǒng)計量Z=-2.18，|Z|>1.96，表明年徑流量在95%置信水平上呈顯著減少趨勢。顯著水平臨界線內(nèi)出現(xiàn)的兩個交點分別位于1991—1992年和2001—2002年，即年徑流量突變點為1991年、2001年。

圖3 年徑流量變化過程

圖4 年徑流量M-K檢驗曲線

3.1.3 年輸沙量變化分析 流域年輸沙量變化過程線見圖5，M-K檢驗曲線見圖6，統(tǒng)計量Z=-3.13，|Z|>2.32，表明年輸沙量在99%置信水平上呈顯著減少趨勢。顯著水平臨界線內(nèi)出現(xiàn)一個交點，位于1982—1983年，故1982年是年輸沙量的突變點。

圖5 年輸沙量變化過程

圖6 年輸沙量M-K檢驗曲線

年降雨量的突變不受人為因素控制，但年輸沙量在1982年發(fā)生突變的原因是70年代大量修建的淤地壩在這個時期產(chǎn)生效用，從而大幅減少水土流失。通過調(diào)查岔巴溝流域水土保持措施得知[31]，1991年后，岔巴溝流域再無新建淤地壩，同時岔巴溝流域Landsat MSS/TM遙感資料[32]顯示，近年來植被恢復(fù)措施也在大力推進，植被覆蓋率逐漸由中低覆蓋度向中高覆蓋度轉(zhuǎn)化，至2007年，中覆蓋度、中高覆蓋度及高覆蓋度分別為51.33%，26.95%和20.65%。這就是年徑流量產(chǎn)生突變在1991年后短暫上升后又下降的原因。通過M-K檢驗得知，年降雨量呈增加趨勢，而徑流量和輸沙量卻始終呈現(xiàn)減少趨勢，這也印證了淤地壩建設(shè)以及植被恢復(fù)等水土保持措施對控制流域水土流失所起的積極作用。

3.2 次洪水研究變量趨勢分析

采用年輸沙量突變點1982年將研究時段劃分為兩個時期，第一時期是1970—1981年，共60場洪水，考慮到剩余時期時間序列過長，將其劃分為第二時期和第三時期，第二時期是1982—2000年，共28場洪水，第三時期是2001—2016年，共18場洪水。3個時期研究變量的平均值見表3。

表3 研究變量平均值

由表3可知，輸沙量和徑流系數(shù)呈減小趨勢，且減小幅度較大，與1970—1981年相比，2001—2016年輸沙量和徑流系數(shù)分別減少了約67%和83%；洪峰流量、徑流深、降雨量、降雨歷時在1982—2000年降低，后在2001—2016年又增加，但洪峰流量仍低于1970—1981年均值，約減少13%，徑流深、降雨量、降雨歷時均大于1970—1981年均值，分別增加了約40%，69%和55%。

比較不同時間尺度下研究變量的變化趨勢，發(fā)現(xiàn)降雨量和輸沙量在年尺度和次洪水尺度下表現(xiàn)一致，但徑流深卻在次洪水尺度下呈增加趨勢，究其原因，雖近年來徑流總量減少，但洪水發(fā)生次數(shù)也在減少，2001—2016年徑流總量約占70年代徑流總量的3/4，而發(fā)生的洪水次數(shù)僅為70年代的1/3，其減小幅度遠(yuǎn)大于徑流總量的減小幅度。

3.3 不同時間尺度水沙關(guān)系回歸分析

3.3.1 次洪水時間尺度回歸分析 構(gòu)建岔巴溝流域各研究時期次洪水降雨、徑流和泥沙變量的皮爾遜相關(guān)系數(shù)矩陣，1970—1981年的計算結(jié)果詳見表4。

表4 1970-1981年次洪水研究變量的相關(guān)系數(shù)矩陣

1970—1981年ESY與Qw，H和P具有極顯著的相關(guān)性，與D呈顯著相關(guān)性，表明Qw，H，P以及D在此研究時期是ESY的主要影響因素。1982—2000年ESY與Qw，H和P具有極顯著的相關(guān)性，2001—2016年ESY僅與Qw和H具有極顯著的相關(guān)性，與P呈顯著相關(guān)性，選擇上述與輸沙量具有顯著相關(guān)性的研究變量進行回歸分析，回歸結(jié)果見表5。

表5 研究變量與次洪水輸沙量的回歸方程

ESY在3個時期均與Qw，H以及P具有顯著的相關(guān)性，但整體來看，ESY與Qw和H的相關(guān)關(guān)系最為顯著和穩(wěn)定，雖然決定系數(shù)有所浮動，但始終在p<0.01水平上呈顯著相關(guān)性。從而建立ESY與Qw和H的回歸方程，回歸模型的決定系數(shù)均大于0.9，說明ESY始終與Qw和H存在良好的冪函數(shù)關(guān)系。H和Qw可以很好的反映輸沙能力，這在許多學(xué)者的研究成果中已經(jīng)體現(xiàn)，但與Zhang[11]、Zheng等[15]分析單個因素與輸沙量呈線性關(guān)系或比例關(guān)系不同，本文研究得出洪峰流量和徑流深與輸沙量呈冪函數(shù)關(guān)系，這與肖學(xué)年等人[10]的分析一致。考慮到本文采用H作為影響因素分析，與P，D同為次降雨指標(biāo)，它們之間具有相關(guān)性，故未將P和D作為回歸因子進行分析。

回歸方程中，回歸系數(shù)代表研究變量對輸沙量的影響程度，回歸系數(shù)在1982—2000年最大，在2001—2016年最小，表明上述研究變量在3個時期對輸沙量的影響整體呈先增大后減小的趨勢；變量所對應(yīng)的指數(shù)越大，表明它對輸沙量的影響越大，前兩個時期徑流深對輸沙量的影響大于洪峰流量，然而2001—2016年洪峰流量的指數(shù)明顯增大，徑流深則無顯著變化，表明洪峰流量對輸沙量的影響更大。同時，回歸模型的決定系數(shù)逐時期下降，說明岔巴溝流域的水沙關(guān)系逐漸變的復(fù)雜，而在上述影響輸沙量的相關(guān)變量的研究中，并未將水土保持因子考慮在內(nèi)，尤其是隨著時間推移，淤地壩與植被恢復(fù)工程發(fā)揮的水土保持效益應(yīng)呈現(xiàn)相反態(tài)勢，這也使得2001—2016年回歸模型的決定系數(shù)低于前兩個時期。

根據(jù)Li等[33]研究，暴雨是造成黃土高原水土流失的主要原因，一場嚴(yán)重暴雨造成的土壤流失可占年土壤流失的60%～90%。采用極端降雨閾值[34]評估岔巴溝流域3個時期的降雨量，具體而言，就是將研究時段內(nèi)降雨量大于0.1 mm的降雨事件升序排列，將第90個百分位數(shù)值認(rèn)定為該站點研究時段內(nèi)極端降雨的閾值，當(dāng)次降水量超過此閾值時，便認(rèn)為此次降雨為極端降雨事件。極端降雨事件在3個研究時期的發(fā)生頻率分別為10%，10.7%和11.1%，相關(guān)特征值見表6。

表6 極端降雨事件特征值

對比不同時期相似的極端降雨條件下的徑流量和輸沙量，如降雨事件2009-07-19與1987-08-26，以及2006-08-29與1978-08-07，發(fā)現(xiàn)2001—2016年極端降雨事件所產(chǎn)生的徑流量和輸沙量明顯減少，表明水土保持措施對減少徑流泥沙起著重要作用，但輸沙量占總輸沙量的比例卻明顯增多，分別為7%，18%和43%。在發(fā)生頻率大致相同的情況下，輸沙量占總輸沙量的比例遠(yuǎn)大于前兩個時期，這表明極端降雨事件在2001—2016年內(nèi)對流域水沙關(guān)系的影響明顯變大，導(dǎo)致水沙關(guān)系發(fā)生了顯著變化，極大的增加了水沙關(guān)系預(yù)估的不確定性，這與Zhao等[34]研究所得結(jié)論一致。

3.3.2 年時間尺度回歸分析 構(gòu)建岔巴溝流域各研究時期年降雨、徑流和泥沙變量的皮爾遜相關(guān)系數(shù)矩陣，1970—1981年的計算結(jié)果詳見表7。

表7 1970-1981年研究變量的相關(guān)系數(shù)矩陣

1970—1981年輸沙量與H和Qwa具有極顯著的相關(guān)性，與Rca呈顯著相關(guān)性，表明H，Qwa和Rca在此研究時期是輸沙的主要影響因素。1982—2000年輸沙量與H呈極顯著的相關(guān)性，與Rca呈顯著相關(guān)性，2001—2016年輸沙量與Qwa呈極顯著的相關(guān)性，與I24 h呈顯著相關(guān)性。對3個時期的年尺度水沙關(guān)系分別建立回歸函數(shù)模型，結(jié)果見表8。

表8 研究變量與年平均輸沙量的回歸方程

考慮到Pa，Rca與H同為次降雨指標(biāo)，故未將Pa和Rca作為回歸因子進行分析?；貧w方程中決定系數(shù)均大于0.96，說明ESYa始終與Qwa和Ha存在良好的多元冪函數(shù)關(guān)系。與目前Zheng等[15]在年時間尺度上通過單因素分析得出的線性關(guān)系相比，研究變量為徑流深和最大洪峰流量時回歸方程的決定系數(shù)更大，回歸效果更好，同時回歸模型也與次洪水尺度下相同，均為冪函數(shù)。然而回歸系數(shù)整體呈降低趨勢，表明研究變量對輸沙量的影響越來越??；在前兩個時期中徑流深的指數(shù)大于最大洪峰流量的指數(shù)，表明徑流深對輸沙量的影響更大，而2001—2016年回歸關(guān)系發(fā)生改變，最大洪峰流量的指數(shù)遠(yuǎn)大于徑流深的指數(shù)，此時最大洪峰流量對輸沙量的影響更大。

上述回歸分析表明，不同時間尺度下水沙關(guān)系變化趨勢不同。在年尺度下回歸方程的決定系數(shù)呈增加態(tài)勢，而次洪水尺度下呈降低趨勢。這種差異主要是由于研究數(shù)據(jù)和變量的時間尺度不同而導(dǎo)致的，年尺度回歸分析中采用的研究事件除已選取的次洪水事件之外，還包含其他因洪峰流量小、降雨歷時短以及缺乏降雨資料等原因未定義為一次洪水事件的徑流過程。統(tǒng)計所選取的次洪水事件，發(fā)現(xiàn)3個研究時期中次洪水事件的總輸沙量占年平均輸沙量的比例分別為76%，71%和59%。由此可見，年尺度回歸分析同次洪水回歸分析相比，關(guān)于流域輸沙過程研究的全面性和普遍性逐時期遞增，所以年尺度回歸分析的決定系數(shù)逐漸增加并大于次洪水回歸分析的決定系數(shù)。然而年尺度回歸分析無法細(xì)致的表征洪水事件的相關(guān)特征變量，這從回歸系數(shù)數(shù)值較低且逐時期減小的趨勢也可以看出，表明年尺度分析對流域水沙關(guān)系的研究無法深入，故目前的研究多集中于次洪水尺度下的水沙關(guān)系分析。

4 結(jié) 論

(1) 岔巴溝流域1970—2016年降雨量呈顯著增加趨勢，年徑流量和年輸沙量呈顯著下降趨勢，年輸沙量、年降雨量和年徑流量分別于1982年、1986年、1991年、2001年發(fā)生突變，106場次洪水事件分析結(jié)果表明降雨量和輸沙量的變化趨勢與年尺度上表現(xiàn)一致，但徑流深呈增加趨勢，主要是洪水次數(shù)大幅減少所致。

(2) 在次洪水時間尺度上，輸沙量與洪峰流量和徑流深始終呈現(xiàn)顯著相關(guān)關(guān)系，其中洪峰流量對輸沙量的影響逐漸增大，徑流深則無顯著變化。

(3) 在年時間尺度上，輸沙量整體與最大洪峰流量以及徑流深關(guān)系顯著，前期徑流深對輸沙量的影響較大，后期最大洪峰流量對輸沙量的影響較大。