蘇婷婷，吳志勇，何 海，孫昭敏

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

據(jù)研究，林草過牧濫伐、城鎮(zhèn)無序擴(kuò)張等不合理的人類活動是引起土壤侵蝕的主要原因，適度控制人類活動、提高植被覆蓋率是土壤侵蝕治理及水土流失防治的重要手段［1，2］。20 世紀(jì) 80 年代中期以來，由于黃河流域降雨等自然因素和水利、水土保持工程等人類活動的影響，黃河天然徑流量和來沙量急劇減少，水土流失得到改善，同時水沙關(guān)系發(fā)生重大變化，產(chǎn)生一系列新問題，如下游河槽萎縮，河道輸沙能力顯著降低，威脅流域的防洪安全，影響區(qū)域經(jīng)濟(jì)的可持續(xù)發(fā)展。治理水土流失事關(guān)經(jīng)濟(jì)社會可持續(xù)發(fā)展和中華民族長遠(yuǎn)福祉，因此，研究黃河流域侵蝕產(chǎn)沙是新時期黃河治理亟需解決的基礎(chǔ)問題。

現(xiàn)有研究主要關(guān)注的是黃河沙量銳減的原因及未來某一降雨情景下的多年平均來沙量預(yù)測［3-6］，對場次暴雨可能來沙量預(yù)測卻少有研究。侵蝕產(chǎn)沙模型是研究流域侵蝕產(chǎn)沙的重要工具。目前，流域產(chǎn)沙模型主要分為經(jīng)驗統(tǒng)計模型和物理成因模型。物理成因模型從大量觀測和定量描述入手，基于物理基礎(chǔ)建立數(shù)學(xué)模型模擬土壤侵蝕過程，但物理成因復(fù)雜，加上水沙動力過程的復(fù)雜性，模擬結(jié)果依然存在諸多不確定因素。與物理成因模型相比，經(jīng)驗統(tǒng)計模型結(jié)構(gòu)簡單、使用方便，仍然是目前流域侵蝕產(chǎn)沙研究的主要工具［5］。黃土丘陵溝壑區(qū)地貌復(fù)雜，真實下墊面參數(shù)獲取難度大［6］，眾多學(xué)者根據(jù)各研究區(qū)域的產(chǎn)輸沙特點，提取影響土壤侵蝕和流域產(chǎn)沙的主要影響因素，建立了不同的統(tǒng)計模型［4，7-17］。如，江忠善等［7］基于陜北、晉西、隴東南黃土丘陵區(qū)溝道小流域?qū)崪y資料，將洪水徑流總量、流域平均坡度、黃土中沙粒和粉粒含量、植被作用系數(shù)作為產(chǎn)沙量預(yù)報的因子；牟金澤等［8］在岔巴溝流域主要考慮洪水徑流模數(shù)、洪峰模數(shù)、主溝道平均比降及流域長度來構(gòu)造洪水產(chǎn)沙量綜合預(yù)報公式。本研究以窟野河流域為對象，利用流域1961—2016年的水文泥沙數(shù)據(jù)分析水沙關(guān)系變化，劃分不同下墊面代表時期，確定次含沙量的影響因子，分別構(gòu)建不同下墊面代表時期的流域次降雨徑流產(chǎn)沙經(jīng)驗?zāi)Ｐ停瑥牟煌a(chǎn)沙水平及不同降雨等級2 種角度分析模型模擬效果，為客觀認(rèn)識黃河流域侵蝕產(chǎn)沙特點提供科學(xué)支撐。

1 研究區(qū)及資料

窟野河是黃河中游水土流失嚴(yán)重的一條多沙粗沙支流，發(fā)源于內(nèi)蒙古南部沙漠地區(qū)，于陜西省神木縣匯入黃河。干流全長242 km，流域面積為8 706 km2。流域多年平均降水量為410 mm，多年平均徑流量為7.47 億 m3，年均輸沙量為1.36 億t，占黃河總輸沙量的6.9%。

該流域氣候?qū)俸疁貛Ц珊?、半干旱大陸性季風(fēng)氣候，降水多為暴雨形式，7—9 月降水占全年的68%。地質(zhì)地貌為沙質(zhì)（流域西部及偏北）、礫質(zhì)（流域北部及東北部）及黃土丘陵（南部，典型黃土丘陵溝壑地貌景觀）的組合。

流域內(nèi)共設(shè)有38 個雨量站和6 個水文站，出口控制站為溫家川站。水文資料采用窟野河流域1961—2016 年降水摘錄表、洪水要素摘錄表，包含時段降水量、流量、含沙量等數(shù)據(jù)，由黃河水利委員會水文局提供。降水量摘錄表只記錄汛期數(shù)據(jù)，主要集中在6—9 月。水沙資料主要源于流域出口控制站溫家川站實測徑流泥沙資料。流域概況見圖1。

以窟野河流域溫家川站1961—2016 年夏汛期（6—9 月）的場次洪水為研究對象，洪峰流量大于100 m3/s 的洪水全部入選。共計收集到97 場洪水資料，其中參與建模的1961—1979 年共23 場洪水，1980—1996 年共 20 場洪水，1997—2016 年共 30 場洪水，未參與建模的共24 場洪水，用于驗證模型模擬效果。

圖1 窟野河流域概況

2 次降雨產(chǎn)沙經(jīng)驗?zāi)Ｐ蜆?gòu)建

2.1 不同下墊面代表時期識別

采用降雨徑流泥沙雙累積曲線法分析流域水沙序列突變點，識別不同下墊面的代表時期。雙累積曲線法是利用同期累積降雨量與累積徑流量（累積輸沙量）的關(guān)系曲線斜率變化分析水沙序列的變化趨勢，斜率發(fā)生變化表明由于人類活動或環(huán)境因素導(dǎo)致下墊面產(chǎn)水產(chǎn)沙關(guān)系發(fā)生改變，轉(zhuǎn)折點為突變點，以突變點為準(zhǔn)劃分為不同的時期。

圖2 為場次洪水累積雨量-場次洪水累積沙量雙累積曲線，圖3 為年徑流量和年輸沙量的累積距平曲線，2 條曲線均在1979 年和1997 年出現(xiàn)轉(zhuǎn)折點，以此將徑流和輸沙序列劃分為3 個階段。將1961—1979 年作為天然時期，1997 年之后流域產(chǎn)沙環(huán)境大體穩(wěn)定，故1997—2016 年作為現(xiàn)狀下墊面代表時期。據(jù)研究，黃河流域產(chǎn)沙環(huán)境在20 世紀(jì)70 年代以前持續(xù)惡化，70 年代以后逐漸好轉(zhuǎn)，1997 年以后顯著改善，因此1997 年之后產(chǎn)沙量和產(chǎn)沙強度大幅降低［15］，此結(jié)論與本研究相符。

圖2 窟野河流域場次洪水累積雨量-場次洪水累積沙量雙累積曲線

圖3 1960—2016 年窟野河流域年徑流量和年輸沙量累積距平曲線

2.2 水沙關(guān)系分析

降雨產(chǎn)生徑流，徑流是泥沙的載體。將流域內(nèi)場次降雨面雨量與次洪沙量一一對應(yīng)，生成面雨量和次洪沙量關(guān)系散點圖（圖4）。通過對比1961—2016 年3 個不同時期面雨量和次洪沙量關(guān)系散點位置，可直觀反映流域的面雨量和次洪沙量關(guān)系的變化。若面雨量和次洪沙量關(guān)系散點位置偏上，在面雨量相同的前提下，則表明該時期的輸沙量較多，反之則較少。分析不同時期面雨量與次洪沙量關(guān)系（圖4）可知，隨時間推移，點群向右下偏移，呈現(xiàn)3 個帶狀分布，說明在現(xiàn)狀下墊面條件下，同樣量級的降雨，與天然時期相比產(chǎn)沙量明顯減少。以2012 年7月和2016 年8 月發(fā)生的2 場洪水為例，次洪沙量點據(jù)明顯向右下偏離歷史洪水形成的趨勢線。

圖4 不同時期面雨量與次洪沙量的關(guān)系

分析1961—2016 年3 個不同時期次洪水量與次洪沙量關(guān)系（圖5）可知，不同時期次洪水量與次洪沙量關(guān)系發(fā)生明顯變化，點群逐漸向X 軸方向偏移，近年來相同量級場次洪水的含沙量與天然時期相比顯著減小。

圖5 不同時期次洪水量與次洪沙量的關(guān)系

2.3 次洪含沙量影響因素分析

由次洪輸沙量的定義可知，徑流事件中一次暴雨徑流的水沙有如下關(guān)系：

式中，M為次洪產(chǎn)沙量（kg），C為次洪含沙量（kg/m3），Q為次洪徑流量（m3），等式兩邊同除以流域面積得：

式中，Ms為輸沙模數(shù)（t/km2）；Cs為次洪平均含沙量（kg/m3）；H為徑流深（mm）。

流域侵蝕產(chǎn)沙是個復(fù)雜的物理過程，是流域下墊面包括地貌、土壤、植被、水保措施等以及降水因子綜合作用的結(jié)果，而流域侵蝕產(chǎn)沙會輸出徑流量和次洪含沙量，因此次洪平均含沙量和徑流深可以完整地表達(dá)流域產(chǎn)沙［18］，所以式（2）利用徑流深和影響次洪含沙量的因子的線性關(guān)系來模擬產(chǎn)沙是一種可行的研究方法。

次洪含沙量與次洪水流挾沙能力和流域可被攜帶土壤分散量有關(guān)。由于黃河流域沙源豐富，故水流挾沙能力決定了次洪含沙量的大小。水動力學(xué)研究表明，水流挾沙能力與流量關(guān)系密切［19-22］。參與相關(guān)分析的因子有反映洪水徑流特征的平均流量、洪峰流量、徑流總歷時及反映徑流輸沙特征的最大含沙量。采用1997—2016 年30 場洪水參與相關(guān)分析，采用SPSS 軟件進(jìn)行相關(guān)分析，結(jié)果（表1）表明，次洪含沙量與洪峰流量、最大含沙量、徑流總歷時在0.01 水平上顯著相關(guān)，故選用這3 個要素指標(biāo)和徑流深模擬產(chǎn)沙。

表1 窟野河流域次洪平均含沙量與其影響因素的相關(guān)系數(shù)

2.4 模型構(gòu)建

基于公式（2），采用徑流深和與次洪含沙量相關(guān)的最大含沙量、洪峰流量、徑流總歷時的非線性函數(shù)形式模擬產(chǎn)沙，利用SPSS 軟件中的非線性回歸過程模塊確定輸沙模數(shù)與各因子之間的關(guān)系，建立流域輸沙模型。

由輸沙模數(shù)與影響因素相關(guān)關(guān)系（圖6）可知，輸沙模數(shù)與各影響因素間并非簡單的線性相關(guān)。輸沙模數(shù)與最大含沙量的對數(shù)形式擬合較好，與洪峰流量和徑流總歷時的多項式形式擬合較好?；谳斏衬?shù)影響因素相關(guān)關(guān)系及流域輸沙關(guān)系表達(dá)式假設(shè)輸沙模數(shù)與各因子間的模型表達(dá)式為：

式中，Ms為輸沙模數(shù)（t/km2）；H為徑流深（mm）；Cm為最大含沙量（kg/m3）；Qm為洪峰流量（m3/s）；T為徑流總歷時（min），a1、b1、b2、b3為模型參數(shù)。

圖6 窟野河流域輸沙模數(shù)與其影響因素的關(guān)系

建模具體步驟如下：①為所有未知參數(shù)指定一個初始值，將原方程按泰勒級數(shù)展開，并只取一階各項作為線性函數(shù)的逼近，其余各項均歸入誤差；②采用最小二乘法對模型中的參數(shù)進(jìn)行估計，用參數(shù)估計值替代初始值，得到1 個新的曲線函數(shù)；③將新得到的曲線函數(shù)展開，線性化，求出參數(shù)估計值；④如此反復(fù)，直至參數(shù)估計值收斂，得到最優(yōu)解。

基于現(xiàn)狀下墊面代表時期（1997—2016年）30場洪水的模擬結(jié)果見表2。由表2可知，Ms與H、Qm的擬合關(guān)系式的R2=0.960，均高于其余2個模型，可見該模型能解釋96.0%的變異，模型的擬合效果較好。最終確定模型結(jié)構(gòu)為Ms=109.174H+0.113HQm-10.225。

表2 基于現(xiàn)狀下墊面（1997—2016 年）代表時期的產(chǎn)沙模型擬合關(guān)系式

以同樣的方法在1961—1979 年（共挑選23 場洪水參與建模）、1980—1996 年（挑選20 場洪水參與建模）分別構(gòu)建模型，模型表達(dá)式見表3。

表3 不同時期產(chǎn)沙模型擬合關(guān)系式

3 模型驗證

3.1 不同產(chǎn)沙水平模擬效果分析

為了驗證不同產(chǎn)沙水平的模型模擬效果，從窟野河1961—2016 年的97 場長系列輸沙資料中選擇未參與模型構(gòu)建的24 場洪水，其中1961—1979 年11 場，1980—1996 年 7 場，1997—2016 年 6 場。以常規(guī)分類方法將產(chǎn)沙水平分為3類，分別為<500 t/km2、500～1 000 t/km2、>1 000 t/km2，驗證模型對不同產(chǎn)沙水平次降雨產(chǎn)沙量的模擬效果，模擬結(jié)果如圖7 所示。由圖7 可知，模型在不同的產(chǎn)沙水平的模擬誤差上表現(xiàn)出一定的變異特性，總體而言，平均相對誤差絕對值為30.6%，相對誤差最小達(dá)9.7%，隨著產(chǎn)沙水平的增大，模型模擬相對誤差逐漸減小。3個時期不同產(chǎn)沙水平的平均相對誤差絕對值見表4，模擬規(guī)律與長系列輸沙資料一致，具體表現(xiàn)為隨著產(chǎn)沙水平的增大，模型模擬相對誤差逐漸減小。上述規(guī)律說明小產(chǎn)沙事件與大產(chǎn)沙事件的產(chǎn)沙機制有所差異，由于黃土高原丘陵溝壑區(qū)高含沙水流的存在，使得不同產(chǎn)沙水平下流量含沙量之間的作用發(fā)生變化。綜上可知，模型對大產(chǎn)沙事件的模擬效果優(yōu)于小產(chǎn)沙事件。

圖7 1961—2016 年窟野河流域24 場不同水平產(chǎn)沙事件與產(chǎn)沙模型的相對誤差

表4 不同時期不同產(chǎn)沙水平模型模擬平均相對誤差絕對值

3.2 不同降雨等級模型模擬效果分析

為驗證所建立的次降雨產(chǎn)沙經(jīng)驗?zāi)Ｐ湍M不同等級降雨侵蝕產(chǎn)沙的效果，從窟野河1961—2016 年的97 場長系列輸沙資料中選擇未參與模型構(gòu)建的24 場洪水，其中 1961—1979 年 13 場，1980—1996 年6 場，1997—2016 年5 場。按不同降雨等級將洪水進(jìn)行分類分析，分類標(biāo)準(zhǔn)：中小雨，1 d（24 h）降雨量<25 mm；大雨，1 d（24 h）降雨量在25～50 mm；暴雨，1 d（24 h）降雨量在50～100 mm。用3 個不同時期的模型分別模擬對應(yīng)時期的場次洪水產(chǎn)沙量，將實測值與模型模擬值作對比，模擬相對誤差見圖8，平均相對誤差見表5?？傮w而言，模型模擬平均相對誤差為32.8%，相對誤差最小達(dá)0.03%，模型對大雨及暴雨的模擬效果略優(yōu)于中小雨。以19910720 次（1991年7 月20 日）洪水為例，本時期模型模擬輸沙模數(shù)為273 t/km2，1960—1979 年模型模擬輸沙模數(shù)為725.3 t/km2，同樣等級的降雨不同時期模型模擬值的差異除了模型本身誤差，可認(rèn)為是不同時期下墊面差異導(dǎo)致的。

圖8 1961—2016 年24 場不同降雨等級典型產(chǎn)沙事件模擬結(jié)果相對誤差

表5 不同時期不同降雨等級模型模擬平均相對誤差絕對值

3.3 不同時期模型參數(shù)差異分析

不同時期次降雨徑流產(chǎn)沙模型參數(shù)差異的原因：①分析輸沙模數(shù)與模型自變量（洪峰流量、徑流深）在1961—1996 年和1997—2016 年的變化趨勢可知（圖9），1961—1996 年洪峰流量、輸沙模數(shù)及徑流深總體變化不大，1997 年之后顯著減小，變化明顯；②分析不同時期的模型參數(shù)，發(fā)現(xiàn)1961—1979 年與1980—1996 年的模型參數(shù)相差無幾，1997—2016 年模型項的系數(shù)明顯增大，導(dǎo)致輸沙模數(shù)急劇減小。據(jù)研究［23］，窟野河流域 1979 年以前，年徑流量和輸沙量處于上升趨勢，該階段受人類活動影響較??；1980—1996 年，年徑流量和輸沙量呈微弱下降趨勢，主要是20 世紀(jì)80 年代開展了水土保持治理及80 年代末開始的煤礦開采等造成的影響；1997 年以后，年徑流量和輸沙量大幅下降，究其原因主要是20 世紀(jì)末大規(guī)模的煤礦開采對流域水文地質(zhì)造成一定的影響，進(jìn)而影響徑流和輸沙過程。

圖9 輸沙模數(shù)、洪峰流量及徑流深在不同時期的變化趨勢

4 小結(jié)

1）雙累積曲線、累積距平曲線均在1979 年和1997 年出現(xiàn)轉(zhuǎn)折點，以此將徑流和輸沙序列劃分為3 個階段。將 1961—1979 年作為天然時期，1997 年之后流域產(chǎn)沙環(huán)境大體穩(wěn)定，1997—2016 年為現(xiàn)狀下墊面代表時期。

2）研究基于Ms=CsH，采用非線性回歸法，構(gòu)建以徑流深與洪峰流量為變量的流域場次洪水產(chǎn)沙統(tǒng)計模型。模型R2均大于0.9，平均相對誤差絕對值在30%左右；通過模擬不同時期不同產(chǎn)沙水平的場次洪水可知，隨著產(chǎn)沙水平的增大，模型模擬相對誤差逐漸減小，模型對大產(chǎn)沙事件的模擬效果優(yōu)于小產(chǎn)沙事件；通過模擬不同等級降雨侵蝕產(chǎn)沙發(fā)現(xiàn)，模型對大雨及暴雨的模擬效果略優(yōu)于中小雨，為客觀認(rèn)識黃河流域侵蝕產(chǎn)沙特點提供科學(xué)支撐。