窟野河流域基于物理成因概念的次降雨徑流產(chǎn)沙統(tǒng)計(jì)模型研究

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

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

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

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

1 研究區(qū)及資料

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

該流域氣候?qū)俸疁貛Ц珊怠敫珊荡箨懶约撅L(fēng)氣候，降水多為暴雨形式，7—9 月降水占全年的68%。地質(zhì)地貌為沙質(zhì)（流域西部及偏北）、礫質(zhì)（流域北部及東北部）及黃土丘陵（南部，典型黃土丘陵溝壑地貌景觀）的組合。

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

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

圖1 窟野河流域概況

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

2.1 不同下墊面代表時(shí)期識(shí)別

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

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

圖2 窟野河流域場(chǎng)次洪水累積雨量-場(chǎng)次洪水累積沙量雙累積曲線(xiàn)

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

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

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

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

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

圖5 不同時(shí)期次洪水量與次洪沙量的關(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)沙是個(gè)復(fù)雜的物理過(guò)程，是流域下墊面包括地貌、土壤、植被、水保措施等以及降水因子綜合作用的結(jié)果，而流域侵蝕產(chǎn)沙會(huì)輸出徑流量和次洪含沙量，因此次洪平均含沙量和徑流深可以完整地表達(dá)流域產(chǎn)沙［18］，所以式（2）利用徑流深和影響次洪含沙量的因子的線(xiàn)性關(guān)系來(lái)模擬產(chǎn)沙是一種可行的研究方法。

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

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

2.4 模型構(gòu)建

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

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

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

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

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

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

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

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

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

3 模型驗(yàn)證

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

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

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

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

3.2 不同降雨等級(jí)模型模擬效果分析

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

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

表5 不同時(shí)期不同降雨等級(jí)模型模擬平均相對(duì)誤差絕對(duì)值

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

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

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

4 小結(jié)

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

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