黃土坡面細溝溝頭溯源侵蝕的量化研究

覃　超，何　超，鄭粉莉※，韓林峰，曾創(chuàng)爍

（1. 西北農(nóng)林科技大學水土保持研究所高原土壤侵蝕與旱地農(nóng)業(yè)國家重點實驗室，楊凌 712100；　2. 水利部黃土高原水土流失過程與控制重點實驗室，鄭州 450003；　3. 美國密西西比大學國家水科學計算與工程中心，牛津 38677；4. 重慶交通大學河海學院，重慶400074；5. 美國密西西比大學工程學院，牛津 38677）

0　引　言

細溝侵蝕是坡面侵蝕的重要方式，坡面一旦形成細溝，坡面流即由細溝間水流匯集成集中的細溝水流，坡面徑流流速和流深均不同程度增加，導致坡面土壤侵蝕加劇。在坡面細溝和細溝間侵蝕蝕區(qū)，細溝侵蝕一般占坡面侵蝕量的70%以上[1]。細溝發(fā)育包括溝頭溯源、溝底下切和溝壁擴張中的一個或多個子過程，而細溝發(fā)育不同階段的主導侵蝕方式不同[2-3]。溝頭溯源侵蝕作為一個高度紊亂的三維侵蝕體系[4]，是細溝侵蝕的開始，在坡面侵蝕產(chǎn)沙中占有重要地位[5]。在溝頭溯源至臨界坡長前，溝頭溯源侵蝕量占坡面細溝侵蝕量的比例可達60%[5-6]。

基于野外調(diào)查、定位觀測和航片解譯等手段，許多學者對野外坡面細溝溝頭在一場或多場降雨、匯流過程中的溯源過程、溝長增加速率、產(chǎn)沙貢獻、影響因素及溝頭的潛在發(fā)生位置和提取方法等進行了研究，取得了一批前瞻性成果[7-9]。Alonso等[10]在研究緩坡細溝溝頭溯源侵蝕時指出，溝頭溯源過程主要決定于土壤可蝕性，尤其是土壤中的黏粒含量，而上方匯流量對溯源過程的影響相對較小。Wirtz等[8]認為溝頭溯源侵蝕是使細溝水流含沙濃度大于徑流挾沙能力的原因之一，也是坡面侵蝕產(chǎn)沙的主要來源。Jia等[4]和Zhu等[11]分別用數(shù)值模擬和物理模型模擬的方法研究了流量、跌坎高度和土壤顆粒組成及土壤分層對溝頭溯源侵蝕的影響，豐富了溝頭溯源侵蝕的研究方法并取得了與物理實驗相近的模擬結(jié)果。為更好剖析溝頭溯源侵蝕機理，前人在室內(nèi)建造實體坡面模型和溝頭模型，以小尺度的水槽試驗對影響溝頭溯源過程的各個因素（上方匯流量和含沙濃度、跌坎高度、近地表水文條件、土壤粒級組成等）及溝頭的發(fā)生位置、間距等進行研究，闡明溝頭溯源侵蝕的過程機理[12-19]。他們的研究結(jié)果表明，溝頭溯源侵蝕距離與時間成線性函數(shù)增加趨勢，且細溝長度與土壤粒級分布、上方匯流流量及含沙濃度、土壤固結(jié)時間和跌坎高度有關(guān)[20]。由于溝頭形態(tài)及其附近水流流態(tài)的復雜性[4,21]和溯源侵蝕速率的多變性[5,12,22]，現(xiàn)有的侵蝕預報模型很難準確預測溝頭附近區(qū)域的土壤侵蝕量，因而現(xiàn)有模型在模擬溝頭溯源侵蝕時常常做了簡化和假設(shè)[5]。

現(xiàn)有溝頭溯源侵蝕研究多側(cè)重于溝頭形成及其溯源特征，而對溝頭溯源侵蝕過程中細溝溝槽底部進一步發(fā)育的二級溝頭的侵蝕機理及其對產(chǎn)沙的影響較少涉及[5,12-14,22-23]。黃土高原作為世界上水土流失最嚴重的地區(qū)之一，其主要土壤——黃綿土具有質(zhì)地輕、疏松綿軟、團粒結(jié)構(gòu)和抗蝕性差等特點[24]，該地區(qū)地形破碎而陡峻，暴雨集中，導致黃土坡面的土壤侵蝕往往由夏、秋季幾場侵蝕性暴雨造成。暴雨匯流使陡坡坡面集中水流能量被細溝溝頭溯源侵蝕消耗后又迅速增加，極易在一級溝頭下方區(qū)域生成新的二級溝頭，形成新的溝頭溯源侵蝕，具有細溝一旦產(chǎn)生就迅速發(fā)育的特點[25]。此外，由于一次暴雨過程中，坡面細溝發(fā)育速度快，其溝槽形態(tài)具有寬度窄和深度大的特征，且細溝水流含沙濃度大，導致動床條件下實時動態(tài)準確監(jiān)測溝頭形態(tài)和跌坎上、下方的水動力學特征相當困難[26-27]，阻礙了溝頭溯源侵蝕過程數(shù)據(jù)的準確獲取，也是造成細溝溯源侵蝕量化研究薄弱的主要原因之一。為此，本研究基于立體攝影測量技術(shù)獲取的高精度地面高程模型（DEM），采用人工模擬徑流沖刷的方法，量化細溝溝頭溯源侵蝕過程及其與坡面侵蝕產(chǎn)沙的關(guān)系，探討溝頭溯源侵蝕速率、溝頭跌坎高度和一級溝頭下方發(fā)生的二級溝頭特征對坡面侵蝕產(chǎn)沙的影響，以期為坡面土壤侵蝕過程預報模型的建立提供依據(jù)。

1　試驗材料與方法

1.1　試驗材料與設(shè)備

供試土壤為黃土高原丘陵溝壑區(qū)安塞縣的耕層黃綿土，其中有機質(zhì)含量為5.9 g/kg，黏粒、粉粒與砂粒含量分別為13.6%、58.1%和28.3%。供試土槽共有4個，長200 cm，寬30 cm，深50 cm，一端可升降，坡度調(diào)節(jié)范圍是 0～25°，每次隨機選擇其中的一個土槽進行試驗以減小系統(tǒng)誤差；土槽上方2.5 m處架有恒定水頭的供水裝置，通過調(diào)節(jié)水閥開度控制流量大小，流量調(diào)節(jié)范圍為0～10 L/min；在坡長70和120 cm處的土槽正上方（1.5 m高）平行架有兩臺能手動對焦的數(shù)碼照相機（Canon EOS 5D Mark II），基于立體攝影測量技術(shù)，獲取試驗過程中試驗土壤表面的高精度DEM。

1.2　試驗設(shè)計

根據(jù)黃土高原常見的短歷時、高強度侵蝕性降雨標準（10.5～234.8 mm/h）[28-29]，設(shè)計匯水流量為1、2、3、4 L/min（相當于在10 m匯水坡長、0.3 m坡寬、徑流系數(shù)為0.8的條件下分別發(fā)生25、50、75、100 mm/h降雨強度的侵蝕性暴雨）。根據(jù)已有研究，當?shù)乇砥露却笥?5°后，坡面細溝侵蝕發(fā)育強烈，再者在實施退耕還林工程近20年的今天，陜北無定河等流域還存在大量坡度大于15°的坡耕地，而25°是黃土高原退耕的上限；所以此研究設(shè)計試驗坡度為15°和20°。根據(jù)野外實測資料（黃土坡面耕層深度約15～20 cm），設(shè)計黃綿土填土厚度為20 cm；為模擬黃土坡面連續(xù)細溝形成前溝頭的溯源侵蝕過程，在坡長170 cm處建造高5 cm的雛形溝頭[12-14,25,30]。每個試驗處理重復 2次，本文所述結(jié)果為兩次重復試驗的平均值。整個試驗過程中兩臺相機均由紅外遙控器控制，每隔30～60 s同時拍照一次。

1.3　試驗步驟

1）裝填試驗土槽。試驗土槽深50 cm，其中0～30 cm裝填細沙，30～50 cm裝填黃綿土，裝土容重為 1.15 g/cm3，每5 cm填裝一層；在裝上層土之前，先抓毛下層土壤表面，以減少土壤分層現(xiàn)象。其中，在裝表層5 cm黃綿土時，首先需在坡長170 cm處放置一個細溝溝頭模型并固定其位置（模型長、寬、高分別為30、5和5 cm），然后分兩層在模型上方裝填黃綿土，第一層裝填3 cm，第二層裝填2 cm，且裝填前需將黃綿土與Ca(OH)2粉末混勻（Ca(OH)2的質(zhì)量百分數(shù)為5%）。添加Ca(OH)2粉末的目的是使表層土壤在預降雨時更好地產(chǎn)生結(jié)皮，從而確保一級溝頭上方坡面的土壤剝蝕速率接近于0[31]。建造雛形溝頭的目的是模擬黃土坡面細溝溝頭的溯源侵蝕過程。另外，用制作的這個細溝溝頭模型可保證每次試驗中人為建筑的細溝溝頭完全相同，保證不同試驗處理下墊面的一致性和研究結(jié)果的可比性。在完成試驗土槽裝土后，小心將0～30 cm坡長處的表層10 cm土壤取出，混以水泥（水泥和土的質(zhì)量比為5:2），然后回填坡面，其目的是使試驗段的土壤顆粒與過渡段的混合物顆粒緊密結(jié)合，減少上方匯流對坡頂過渡段的非正常沖刷[12,29-30]（圖1）。

2）為保證試驗前期土壤條件的一致性，使混有水泥的土壤表面充分凝固并使土壤表層形成結(jié)皮，正式降雨的前一天將土槽坡度調(diào)至3°，采用30 mm/h降雨強度進行預降雨至坡面產(chǎn)流為止。預降雨結(jié)束后，將試驗土槽調(diào)平（0°），靜置10 h[12-13]。

3）正式試驗開始前在土槽上方同一高度架設(shè)兩臺數(shù)碼相機，設(shè)置照片的拍攝規(guī)格為RAW，分辨率設(shè)置為相機的最大分辨率（2720×4080），調(diào)節(jié)照相機的方向并使其拍攝角度與坡面始終保持垂直；調(diào)節(jié)相機場景模式至“M手動”，并分別調(diào)節(jié)光圈（f/2.8）、ISO感光度（250）、快門（1/20），然后對焦，待圖像清晰后將對焦模式設(shè)置成手動，并使兩臺相機在整個試驗過程中的焦距保持24 mm不變（圖1）。

圖1　立體攝影測量相機同時拍攝的一組照片和溝頭橫剖面示意圖Fig.1　Example of paired photos obtained from photogrammetry during experiment and sketches of headcut cross section

4）在試驗土槽四周設(shè)置4個固定標靶和3個活動標靶，且使標靶與試驗土槽土壤表面保持平行，并保證任意 3個標靶不在同一條直線上（圖 1），試驗前用鋼尺（1 mm精度）分別測量每個標靶中心點之間的距離，以土槽右岸最下方的標靶中心點為坐標原點（0,0,0），分別計算其它坐標中心點的相對坐標，以保證后期照片的拼接精度。

5）調(diào)整土槽坡度，率定上方匯流量，當率定流量與設(shè)計目標流量的相對誤差小于5%時即可開始正式試驗。

6）試驗開始后即連續(xù)接取徑流泥沙樣并實時稱量徑流桶+徑流泥沙樣質(zhì)量，為使徑流泥沙樣的接樣時間與拍照時間間隔一致，設(shè)計每個徑流泥沙樣的接樣時間為30 s；在拍照間隙，以2 min為時間步長，30～50 cm為測距，分別在溝頭上方和下方循環(huán)測量細溝水流流速以及相對應(yīng)的流寬和流深。其中流速測量采用染色劑示蹤法，每個流速數(shù)據(jù)以 3次重復測量的表層流速的平均值乘以0.75為最終結(jié)果[15]；流深和流寬的測量采用普通鋼尺（1 mm精度），同樣以3次重復測量的平均值為最終結(jié)果。

7）沖刷結(jié)束后將徑流泥沙樣靜置12 h，去除徑流泥沙樣的上層清液并轉(zhuǎn)移至鋁制飯盒，在烘箱內(nèi)（105 ℃）烘干至恒重，然后計算徑流量和侵蝕量。

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

根據(jù)Bryan等和Slattery等[27,32]對A型溝頭（平整坡面上發(fā)育的溝頭）和B型溝頭（A型溝頭溯源過后下方溝槽內(nèi)發(fā)育的溝頭）的定義，本研究將由雛形跌坎發(fā)育而來的初始溝頭定義為一級溝頭，在一級溝頭溯源侵蝕過后產(chǎn)生的溝槽內(nèi)發(fā)育的溝頭定義為二級溝頭（圖1）。根據(jù)Bennet[22]對最大沖刷深度（maximum depth of scour）的定義，本研究將溝頭裂點（knickpoint）到水涮窩最低點的垂直距離定義為溝頭跌坎高度（圖1）。

試驗結(jié)束后將拍攝照片導入 Agisoft Photoscan Professional 1.2.4軟件，在完成照片標靶設(shè)置、圖像拼接、校正和坐標系設(shè)置后，生成高密度點云數(shù)據(jù)后導出（.txt格式）；然后將處理后的數(shù)據(jù)導入ArcGIS 10.1軟件，生成漁網(wǎng)并進行空間校正，然后生成高精度DEM（精度為2 mm×2 mm），通過三維分析和表面分析模塊，獲取坡面細溝形態(tài)特征。

為驗證擬合方程的有效性，選取模型有效性系數(shù) R2和納什系數(shù)ENS[33]作為評價指標：

式中Oi為試驗觀測值；Yi為模型預測值；O為試驗觀測值的平均值；Y為模型預測值的平均值；n為樣本個數(shù)。

R2表示觀測值與預測值之間關(guān)系的緊密程度，ENS表示觀測值和預測值在1:1線附近的分布情況。R2和ENS越接近1，表明模型的預測效果越好；R2和ENS趨向于0則表明模型預測值與實際觀測值之間有較大誤差。通常情況下，當R2大于0.6且ENS大于0.5時可認為模型預測達到可接受的精度[34]。

2　結(jié)果與討論

2.1　坡面產(chǎn)沙過程及細溝形態(tài)特征

不同試驗處理下的坡面產(chǎn)流率均小于設(shè)計流量，偏小幅度介于6.5%～28.0%，說明模擬試驗過程中有少量入滲發(fā)生，符合超滲產(chǎn)流規(guī)律，較好地模擬了裸露黃土坡面上發(fā)育細溝的產(chǎn)流過程（表1）。不同試驗處理下的坡面產(chǎn)沙率隨流量和坡度的增加而增大（表1，圖2）。在坡度不變的情況下，流量每增加1 L/min，坡面產(chǎn)沙率增加155～301 g/min，增加幅度介于0.59～5.34倍；流量每增加一倍，坡面產(chǎn)沙率增加155～529 g/min，增加倍數(shù)介于1.86～5.34倍。在流量不變的情況下，坡度增加5°，坡面產(chǎn)沙率增加 26～101 g/min，增加幅度介于 14.0%～89.7%。

表1　不同試驗處理下的坡面產(chǎn)流率和產(chǎn)沙率Table 1　Sloping runoff rate and sediment delivery rate at different treatments

坡面產(chǎn)沙率隨時間總體上呈上升的變化趨勢，但在不同的流量和坡度處理下，其上升的變化過程有所差異（圖2）。當流量小于或等于2 L/min時，坡面產(chǎn)沙率在試驗初期增加較快，隨著試驗的進行，產(chǎn)沙率緩慢上升，并伴有小幅波動。當流量大于2 L/min時，產(chǎn)沙率始終保持上升趨勢，且隨著流量和坡度的增加，波動上升幅度增大。產(chǎn)沙率隨時間變化波動的原因可歸納為以下幾點：1）當流量大于或等于2 L/min時（3 L/min流量和15°坡度處理除外），細溝溝槽內(nèi)均有二級溝頭產(chǎn)生，因此，在一級溝頭溯源侵蝕速率基本不變且溝頭跌坎高度僅緩慢增加的情況下，二級溝頭的生成增加了坡面產(chǎn)沙率。本結(jié)果證實了前人[5,27,32]的觀點，即細溝溝頭的形成是細溝侵蝕初期坡面產(chǎn)沙的主要來源；2）在相同試驗歷時下，大坡度和大流量的試驗處理擁有較長的細溝長度、較復雜的細溝溝槽形態(tài)，以及處于動態(tài)變化的以侵蝕為主的過程或以沉積為主的過程，造成坡面產(chǎn)沙率波動變化，這與Wells等[12-13]在緩坡（＜3°坡面）上的研究結(jié)果類似，他們發(fā)現(xiàn)在溝頭水涮窩下方存在明顯的泥沙沉積區(qū)；3）隨著細溝長度的增加，產(chǎn)沙率同時受溝頭溯源侵蝕速率、溝頭跌坎高度、侵蝕-沉積動態(tài)變化和二級溝頭產(chǎn)生、溯源的影響，因此波動更明顯。

圖2　不同試驗處理下產(chǎn)沙率隨時間的變化Fig.2　Time series of sediment delivery rate under different treatments

不同試驗處理下細溝的寬度、深度及寬深比分別介于3.0～10.0 cm、2.8～20.0 cm和0.5～3.0之間。試驗所得細溝寬度、深度符合前人有關(guān)黃土坡面細溝寬度（3～30 cm）和深度（3～20 cm）的野外調(diào)查和室內(nèi)模擬結(jié)果[1,3,35]。然而，在坡度為20°，流量為3、4 L/min和坡度為15°，流量為4 L/min的試驗處理下，由于一級溝頭下方的細溝溝槽內(nèi)產(chǎn)生了二級溝頭，而二級溝頭的進一步下切，使部分坡長處的細溝深度已達20 cm（溝底已下切至細沙層），由此導致部分坡段細溝寬深比略小于前人調(diào)查和模擬的結(jié)果（0.75～8.75）。

2.2　溝頭溯源侵蝕過程

為研究細溝長度隨時間、流量（Q）和坡度（S）的變化規(guī)律，且盡量減少坡頂過渡段對試驗結(jié)果的影響，繪制了不同試驗處理下細溝溯源長度達到100 cm（選取坡長70～170 cm為試驗段）的溝長變化過程（圖1）。結(jié)果表明，溝頭溯源長度達到100 cm所需的時間隨流量和坡度的增加逐漸縮短（圖3）。在15°坡度和流量為1、2、3、4 L/min時，溝頭溯源長度達到100 cm所需時間分別為38.6、19.6、18.5和15.9 min；在坡度為20°的4種流量下溝頭溯源長度達到100 cm所需時間分別為29.1、16.9、14.9和13.6 min。由此可見，在相同坡度下，當流量大于1 L/min時，溝頭溯源長度達到100 cm所需時間較流量為1 L/min的試驗處理縮短12 min以上；當流量為2、3、4 L/min時，不同試驗處理間溝頭溯源長度達到100 cm所需時間僅相差3.3～3.7 min，且數(shù)據(jù)點的分布在試驗前期十分接近，因此可以認為2 L/min是試驗條件下使溝頭溯源侵蝕速率明顯增加的一個臨界值；坡度對細溝長度隨時間的變化影響十分明顯，且這一影響隨流量的增加逐漸減弱。圖 3還表明，細溝長度隨時間基本呈線性增加趨勢，溝頭溯源侵蝕速率基本保持不變。

圖3　不同試驗處理下溝頭溯源距離達到100 cm的變化過程Fig.3　Time series of headcut retreat till 100 cm under different treatments

參考前人[13-14]的研究結(jié)果，即在一定的流量和坡度條件下，細溝溝頭溯源侵蝕速率隨時間推移基本保持不變，本研究選用線性函數(shù)L=aT來描述細溝長度隨時間的變化規(guī)律?；诖耍槍γ恳粋€試驗處理，擬合了細溝長度隨時間變化的方程（表2）。為擬合一個通用的細溝長度隨時間變化的方程，需要確定系數(shù)a的通用表達式。Zhu等[18]在研究細溝侵蝕與流量和坡度的關(guān)系時指出，坡度×流量的交互作用是預測溝頭溯源侵蝕量的良好因子。本研究通過分析a與單寬流量（q），s，qs和qs2等流量和坡度組合的關(guān)系得知，系數(shù)a受q和s的共同影響，且這一影響隨流量和坡度的增加呈減弱趨勢，因此兩者關(guān)系可由對數(shù)函數(shù)表達。隨機選取12場次的試驗數(shù)據(jù)（1、2、3、6、7、8場次及其重復）進行參數(shù)擬合，發(fā)現(xiàn)a可由下式表達：

式中q為單寬流量，L/min·m；s為坡度，m/m。

式中L為細溝長度，cm；T為試驗歷時，min。

表2　細溝長度隨時間變化方程的擬合及相關(guān)參數(shù)Table 2　Equations of rill length verses time and related parameters

為進一步驗證式（4）預測細溝長度隨時間變化的效果，選取剩余兩組且與建模場次相互獨立的試驗數(shù)據(jù)（4、5場次，用陰影表示）對式（4）進行驗證（表 2）。結(jié)果表明，場次4和5的方程模擬值與實際觀測值均有較好的相關(guān)關(guān)系，平均相對誤差分別為15.7%和 4.0%，最終相對誤差（即細溝溯源長度達到100 cm時的相對誤差）分別為 1.6%和?2.8%，決定性系數(shù) R2分別為 0.995和0.999，納什有效性系數(shù)ENS分別達到0.955和0.995，說明公式（4）符合模型有效性驗證的基本要求（圖 4），對預測不同流量和坡度條件下黃土坡面細溝長度隨時間的變化過程有一定的參考價值。

圖4　細溝長度的預測值和實際觀測值Fig.4　Observed and simulated values of rill length

2.3　溯源侵蝕速率、溝頭跌坎高度和二級溝頭對坡面產(chǎn)沙率的影響

有的研究指出[36-37]，以溝頭溯源侵蝕為主的坡面侵蝕過程其坡面產(chǎn)沙率主要受溝頭溯源侵蝕速率和溝頭跌坎高度影響。Bennett等和Wells等[5,12]在研究溝頭溯源侵蝕過程中溝頭形態(tài)時也指出，溝頭下方水涮窩深度和溝頭溯源速率受上方來水來沙和坡度的共同制約，并顯著影響產(chǎn)沙過程，產(chǎn)沙量、水涮窩深度和溝頭溯源侵蝕速率均隨流量和坡度的增加而增大。基于此，分別統(tǒng)計了不同試驗處理下的平均溝頭跌坎高度和平均溯源侵蝕速率（表3），并將兩者與產(chǎn)沙率進行回歸分析。結(jié)果表明，產(chǎn)沙率與溝頭跌坎高度和溯源侵蝕速率呈極顯著的正相關(guān)關(guān)系（P＜0.001），相關(guān)系數(shù)分別為0.727和0.902。然而，在大流量和大坡度的試驗處理組合下，溝頭跌坎高度并不隨流量和坡度的增加而嚴格增大。特別是當流量達到4 L/min時，15°坡度下的溝頭跌坎高度明顯大于20°坡度下的跌坎高度，而產(chǎn)沙率卻小于坡度為20°時的產(chǎn)沙率。分析原因，發(fā)現(xiàn)在流量和坡度較大時，一級溝頭下方發(fā)育二級溝頭的概率明顯增加，而二級溝頭數(shù)與產(chǎn)沙率也存在顯著相關(guān)關(guān)系（R=0.675，P=0.002）。Slattery等[27]和 Bryan等[32]也認為，二級溝頭的產(chǎn)生是細溝溝槽下切侵蝕的重要組成部分，在一級溝頭尚未溯源至上方匯水面積與坡度平方乘積AS2臨界值的情況下[38]，溝頭溯源侵蝕和溝底下切侵蝕同時主導坡面侵蝕過程，坡面產(chǎn)沙量達到峰值。以上結(jié)果還側(cè)面證實了Gordon等[31,39]的研究結(jié)論，溝頭跌坎高度的增加促使集中水流進入溝頭的入射角度增加，從而增大了固定渦流區(qū)域（captive eddy region）的徑流侵蝕力，使溝頭溯源速率和坡面產(chǎn)沙率增加。

綜上所述，以溝頭溯源侵蝕為主過程的坡面產(chǎn)沙率受溝頭溯源侵蝕速率、溝頭跌坎高度和溝頭下方溝槽內(nèi)發(fā)育的二級溝頭數(shù)及其相關(guān)特征影響，因此，在布設(shè)坡面水土保持措施時需要對溝頭上方的來水進行攔截，以防止溝頭溯源侵蝕并降低溝槽內(nèi)集中水流能量，從而阻止二級溝頭的形成及下切侵蝕的進一步發(fā)生，從而達到防治土壤侵蝕的目的。

表3　不同試驗處理下一級溝頭的溝頭跌坎高度、溯源侵蝕速率和一級溝頭下方發(fā)育的二級溝頭數(shù)Table 3　Height of initial headcut, headcut advancing rate and number of secondary headcuts under different treatments

2.4　以溝頭溯源侵蝕為主的坡面產(chǎn)沙率估算方程擬合

通過產(chǎn)沙率與溝頭跌坎高度、一級溝頭溯源侵蝕速率、二級溝頭數(shù)的回歸分析得知，產(chǎn)沙率隨溝頭跌坎高度或溝頭溯源侵蝕速率的增加呈冪函數(shù)（y=aXb）增加，而產(chǎn)沙率和二級溝頭數(shù)的關(guān)系，由于 0的存在，則根據(jù)情況可用一次函數(shù)或指數(shù)函數(shù)來表達。因此，初步擬定了以下幾種基本形式的方程：

式中H為一級溝頭跌坎高度，cm；V為一級溝頭溯源侵蝕速率，cm/min；N為一級溝頭下方溝槽內(nèi)發(fā)育的二級溝頭數(shù)；a，b，c，d，m，n為待定系數(shù)。

利用Origin9.0軟件的非線性（Nonlinear Curve Fit）擬合功能，隨機選取12組數(shù)據(jù)用于方程擬合（表3，非陰影表示），分別用上述初擬的關(guān)系式進行多元非線性回歸分析，發(fā)現(xiàn)方程（5）的決定性系數(shù) R2最大，達到0.932。該方程中的H和V分別符合冪函數(shù)關(guān)系，反映了溝頭溯源侵蝕過程中的重力侵蝕特征和徑流能量特征，有一定的物理意義[10-11]；為使方程左右兩邊量綱平衡，首先對產(chǎn)沙率（SL）、一級溝頭跌坎高度（H）和一級溝頭溯源侵蝕速率（V）進行無量綱處理：

式中SL為產(chǎn)沙率，g/min；SL′為無量綱的產(chǎn)沙率；ρ為細溝水流密度，取1.0 g/cm3；V0為使溝頭恰好發(fā)生溯源侵蝕的臨界流速，根據(jù)覃超等[40]的數(shù)據(jù)計算，取105.8 cm/min；A為土槽的有效試驗面積，取3 000 cm2；H0為溝頭最大可下切深度，與土層厚度相等，取20 cm；H′為無量綱的一級溝頭跌坎高度；V′為無量綱的一級溝頭溯源侵蝕速率。

將式（8）、（9）、（10）分別代入方程（5），得到其最終表達式

為驗證式（11）的模擬結(jié)果，選取剩余 4場次數(shù)據(jù)（表3）進行模型有效性驗證。結(jié)果表明，R2和ENS分別為0.892和0.799，符合模型有效性驗證的基本要求。由式（11）中的待定系數(shù)a、b和m、n的相對大小關(guān)系可知，溝頭溯源侵蝕速率對產(chǎn)沙率的影響大于溝頭跌坎高度對產(chǎn)沙率的影響，且兩者通過影響一級溝頭下方溝槽內(nèi)發(fā)育的二級溝頭數(shù)，共同決定溝頭溯源侵蝕過程中的坡面產(chǎn)沙率。Wirtz等[8]指出目前細溝侵蝕預報模型多用徑流功率和徑流剪切力代表細溝水流能量來預測細溝侵蝕量，對在坡面侵蝕過程中是否有溝頭產(chǎn)生及其溯源侵蝕的產(chǎn)沙貢獻尚未考慮。因此，上式在一定程度上有助于細化細溝侵蝕預報方程，并為建立基于物理過程的土壤侵蝕預報模型提供借鑒。

3　結(jié)　論

1）坡面產(chǎn)流率、產(chǎn)沙率、細溝溯源侵蝕速率隨流量和坡度的增加而增加，而溝頭跌坎高度和二級溝頭數(shù)卻不隨流量和坡度的增加而嚴格增大。流量每增加1 L/min，產(chǎn)沙率增加 0.59～5.34倍；坡度增加 5°，產(chǎn)沙率增加41.7%～122.1%。坡面產(chǎn)沙率隨時間總體呈上升的變化趨勢，但隨流量和坡度的增大，其上升速率和波動幅度均增大。

2）溝頭溯源長度達到100 cm所需時間隨流量和坡度的增加逐漸縮短，2 L/min是試驗條件下使溝頭溯源侵蝕速率明顯增加的一個流量臨界值，坡度對細溝長度隨時間變化的影響隨流量的增加逐漸減弱。

3）細溝長度隨時間的增加呈一次函數(shù)增加，且增加速率受流量Q和坡度S乘積的影響；以溝頭溯源侵蝕為主的坡面侵蝕過程受溝頭跌坎高度、一級溝頭溯源侵蝕速率、一級溝頭下方溝槽內(nèi)發(fā)育的二級溝頭數(shù)和土壤內(nèi)在性質(zhì)的共同制約，坡面產(chǎn)沙率與上述因子的關(guān)系可由一個多元非線性回歸方程表達。

4）在布設(shè)坡面水土保持措施時，首先需要對一級溝頭上方的坡面來水進行攔截以防止溯源侵蝕，此外還需考慮如何在充分發(fā)育的溝槽內(nèi)布設(shè)溝底防沖措施，以控制集中水流的流量和流速，降低二級溝頭的發(fā)生概率，從而達到防治坡面侵蝕的目的。

[參考文獻]

[1] 鄭粉莉，唐克麗，周佩華. 坡耕地細溝侵蝕影響因素的研究[J]. 土壤學報，1989，26(2)：109－116.Zheng Fenli, Tang Keli, Zhou Peihua. A study on factors affecting rill erosion on cultivated slope land[J]. Acta Pedologica Sinica, 1989, 26(2): 109－116. (in Chinese with English abstract)

[2] Bingner R L, Wells R R, Momm H G, et al. Ephemeral gully channel width and erosion simulation technology[J]. Natural Hazards, 2016, 80(3): 1949－1966.

[3] Qin Chao, Zheng Fenli, Xu Ximeng, et al. A laboratory study on rill network development and morphological characteristics on loessial hillslope[J]. Journal of Soils and Sediments, 2017. https://doi.org/10.1007/s11368-017-1878-y

[4] Jia Y, Kitamura T, Wang S. Numerical simulation of head-cut with a two-layered bed[J]. International Sediment Research, 2005, 20(3): 185－193.

[5] Bennett S J, Alonso C V, Prasad S N, et al. Experiments on headcut growth and migration in concentrated flows typical of upland areas[J]. Water Resources Research, 2000, 36(7):1911－1922.

[6] 韓鵬，倪晉仁，李天宏. 細溝發(fā)育過程中的溯源侵蝕與溝壁崩塌[J]. 應(yīng)用基礎(chǔ)與工程科學學報，2002，10(2)：115－125.Han Peng, Ni Jinren, Li Tianhong. Headcut and bank landslip in rill evolution[J]. Journal of Basic Science and Engineering,2002, 10(2): 115－125. (in Chinese with English abstract)

[7] Cox S E, Booth D T, Likins J C. Headcut erosion in Wyoming’s Sweetwater Subbasin[J]. Environmental Management, 2016, 57(2): 450－462.

[8] Wirtz S, Seeger K M, Ries J B. Field experiments for understanding and quantification of rill erosion processes[J].Catena, 2012, 91: 21－34.

[9] Knapen A, Poesen J. Soil erosion resistance effects on rill and gully initiation points and dimensions[J]. Earth Surface Processes and Landforms, 2010, 35: 217－228.

[10] Alonso C V, Bennett S J. Predicting head cut erosion and migration in concentrated flows typical of upland areas[J].Water Resources Research, 2002, 38(12): 31－39.

[11] Zhu Y H, Visser P J, Vrijling J K. Soil headcut erosion:process and mathematical modeling[M]//Kusuda T,Yamanishi H, Spearman J, et al. Sediment and Ecohydraulics:Intercoh 2005. Amsterdam: Elsvier, 2008: 125－136.

[12] Wells R R, Alonso C V, Bennett S J. Morphodynamics of headcut development and soil erosion in upland concentrated flows[J]. Soil Science Society of America Journal, 2009,73(2): 521－530.

[13] Wells R R, Bennett S J, Alonso C V. Modulation of headcut soil erosion in rills due to upstream sediment loads[J]. Water Resources Research, 2010, 46(12): 1－16.

[14] Wells R R, Bennett S J, Alonso C V. Effect of soil texture,tailwater height, and pore-water pressure on the morphodynamics of migrating headcuts in upland concentrated flows[J]. Earth Surface Processes and Landforms, 2009, 34(14): 1867－1877.

[15] Shen H, Zheng F, Wen L, et al. Impacts of rainfall intensity and slope gradient on rill erosion processes at loessial hillslope[J]. Soil and Tillage Research, 2016, 155: 429－436.

[16] 和繼軍，孫莉英，李君蘭，等. 緩坡面細溝發(fā)育過程及水沙關(guān)系的室內(nèi)實驗研究[J]. 農(nóng)業(yè)工程學報，2012，28(10)：138－144.He Jijun, Sun Liying, Li Junlan, et al. Experimental study on rill evolution process and runoff-sediment relationship for gentle slope[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012,28(10): 138－144. (in Chinese with English abstract)

[17] 王龍生，蔡強國，蔡崇法，等. 黃土坡面細溝形態(tài)變化及其與流速之間的關(guān)系[J]. 農(nóng)業(yè)工程學報，2014，30(11)：110－117.Wang Longsheng, Cai Qiangguo, Cai Chongfa, et al.Morphological changes of rill on loess slope and its relationship with flow velocity[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(11): 110－117. (in Chinese with English abstract)

[18] Zhu J C, Gantzer C J, Peyton R L, et al. Simulated small-channel bed scour and head cut erosion rates compared[J]. Soil Science Society of America Journal, 1995,59: 211－218.

[19] Wells R R, Momm H G, Bennett S J, et al. A measurement method for rill and ephemeral gully erosion assessments[J].Soil Science Society of America Journal, 2016, 80(1): 203.

[20] Babazadeh H, Ashourian M, Shafai-Bajestan M. Experimental study of headcut erosion in cohesive soils under different consolidation types and hydraulic parameters[J]. Environmental Earth Sciences, 2017, 76: 438.

[21] Nan L, Chen A, Xiong D, et al. The effects of film flow on headcut erosion of a gully in the dry-hot valley of Jinsha River, China[J]. Advances in Earth and Environment Science,2013, 189: 1013－1021.

[22] Bennett Sean J. Effect of slope on the growth and migration of headcuts in rills[J]. Geomorphology, 1999, 30(3): 273－290.

[23] Robinson K M, Bennett S J, Casalí J, et al. Processes of headcut growth and migration in rills and gullies[J].International Journal of Sediment Research, 2000, 15(1): 69－82.

[24] 黃自立. 陜北地區(qū)黃綿土分類的研究[J]. 土壤學報，1987，24(3)：266－271.Huang Zili. Discussion of classification of yellow cultivated loessial soils in Shaanxi[J]. Acta Pedologica Sinica, 1987,24(3): 266－271. (in Chinese with English abstract)

[25] He Jijun, Li Xiaojuan, Jia Lijuan, et al. Experimental study of rill evolution processes and relationships between runoff and erosion on clay loam and loess[J]. Soil Science Society of American Journal, 2014, 78: 1716－1725.

[26] R?mkens M J M, Prasad S N, Whisler F D. Process Studies in Hillslope Hydrology[M]//Anderson M. G., Burt T. P.Surface Sealing and Infiltration. New York: John Wiley,1990.

[27] Slattery M C, Bryan R B. Hydraulic conditions for rill incision under simulated rainfall: A laboratory experiment[J].Earth Surface Processes and Landforms, 1992, 17: 127－146.

[28] 周佩華，王占禮. 黃土高原土壤侵蝕暴雨標準[J]. 水土保持通報，1987，7(1)：38－44.Zhou Peihua, Wang Zhanli. Soil erosion storm rainfall standard in the Loess Plateau[J]. Bulletin of Soil and Water Conservation, 1987, 7(1): 38－44. (in Chinese with English abstract)

[29] Qin Chao, Zheng Fenli, Wells Robert R, et al. A laboratory study of channel sidewall expansion in upland concentrated flows[J]. Soil and Tillage Research, 2018, 178: 22－31.

[30] 覃超，鄭粉莉，徐錫蒙，等. 基于立體攝影技術(shù)的細溝與細溝水流參數(shù)測量及特征[J]. 農(nóng)業(yè)機械學報，2016，47(11)：150－156.Qin Chao, Zheng Fenli, Xu Ximeng, et al. Measurement and characteristics of rill geometry and flow parameter based on photogrammetry[J]. Transactions of the Chinese Society of Agricultural Machinery, 2016, 47(11): 150－156. (in Chinese with English abstract)

[31] Gordon L M, Bennett S J, Wells R R, et al. Effect of soil stratification on the development and migration of headcuts in upland concentrated flows[J]. Water Resources Research,2007, 43(7): 1－12.

[32] Brayan R B, Poesen J. Laboratory experiments on the influence of slope length on runoff, percolation and rill development[J]. Earth Surface Processes and Landforms,1989, 14(7): 211－231.

[33] Nash J E, Sutcliffe J V. River flow forecasting through conceptual models part I-A discussion of principles[J].Journal of Hydrology, 1970, 10(3): 282－290.

[34] Santhi C, Arnold J G, Williams J R, et al. Application of a watershed model to evaluate management effects on point and nonpoint source pollution[J]. Transactions of ASAE,2001, 44(6): 1559－1570.

[35] 和繼軍，呂燁，宮輝力，等. 細溝侵蝕特征及其產(chǎn)流產(chǎn)沙過程試驗研究[J]. 水利學報，2013，44(4)：398－405.He Jijun, Lv Ye, Gong Huili, et al. Experimental study on rill erosion characteristics and its runoff and sediment yield process[J]. Journal of Hydraulic Engineering, 2013, 44(4):398－405. (in Chinese with English abstract)

[36] Bennett S J, Casalí J. Effect of initial step height on headcut development in upland concentrated flows[J]. Water Resources Research, 2001, 37(5): 1475－1484.

[37] 張寶軍，熊東紅，楊丹，等. 跌水高度對元謀干熱河谷沖溝溝頭侵蝕產(chǎn)沙特征的影響初探[J]. 土壤學報，2017，54(1)：48－59.Zhang Baojun, Xiong Donghong, Yang Dan, et al.Preliminary study on effects of headcut height on soil erosion and sediment yield at gully heads in Yuanmou dry-hot valley region[J]. Acta Pedologica Sinica, 2017, 54(1): 48－59. (in Chinese with English abstract)

[38] Torri D, Poesen J. A review of topographic threshold conditions for gully head development in different environments[J]. Earth-Science Reviews, 2014, 130: 73－85.

[39] Gordon L M, Bennett S J, Bingner R L, et al. Simulating ephemeral gully erosion in AnnAGNPS[J]. Transactions of the ASABE, 2007, 50(3): 857－866.

[40] 覃超，吳紅艷，鄭粉莉，等. 黃土坡面細溝侵蝕及水動力學參數(shù)的時空變化特征[J]. 農(nóng)業(yè)機械學報，2016，47(8)：146－154，207.Qin Chao, Wu Hongyan, Zheng Fenli, et al. Temporal and spatial variation characteristics of rill erosion and hydrodynamic parameters on loessial hillslope[J].Transactions of the Chinese Society of Agricultural Machinery, 2016, 47(8): 146－154, 207. (in Chinese with English abstract)