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

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

0　引　言

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

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

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

1　試驗(yàn)材料與方法

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

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

1.2　試驗(yàn)設(shè)計(jì)

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

1.3　試驗(yàn)步驟

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

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

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

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

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

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

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

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

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）的定義，本研究將溝頭裂點(diǎn)（knickpoint）到水涮窩最低點(diǎn)的垂直距離定義為溝頭跌坎高度（圖1）。

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

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

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

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

2　結(jié)果與討論

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

不同試驗(yàn)處理下的坡面產(chǎn)流率均小于設(shè)計(jì)流量，偏小幅度介于6.5%～28.0%，說明模擬試驗(yàn)過程中有少量入滲發(fā)生，符合超滲產(chǎn)流規(guī)律，較好地模擬了裸露黃土坡面上發(fā)育細(xì)溝的產(chǎn)流過程（表1）。不同試驗(yàn)處理下的坡面產(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　不同試驗(yàn)處理下的坡面產(chǎn)流率和產(chǎn)沙率Table 1　Sloping runoff rate and sediment delivery rate at different treatments

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

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

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

2.2　溝頭溯源侵蝕過程

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

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

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

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

式中L為細(xì)溝長度，cm；T為試驗(yàn)歷時(shí)，min。

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

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

圖4　細(xì)溝長度的預(yù)測值和實(shí)際觀測值Fig.4　Observed and simulated values of rill length

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

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

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

表3　不同試驗(yàn)處理下一級溝頭的溝頭跌坎高度、溯源侵蝕速率和一級溝頭下方發(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ù)來表達(dá)。因此，初步擬定了以下幾種基本形式的方程：

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

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

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

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

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

3　結(jié)　論

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

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

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

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

[參考文獻(xiàn)]

[1] 鄭粉莉，唐克麗，周佩華. 坡耕地細(xì)溝侵蝕影響因素的研究[J]. 土壤學(xué)報(bào)，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] 韓鵬，倪晉仁，李天宏. 細(xì)溝發(fā)育過程中的溯源侵蝕與溝壁崩塌[J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，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] 和繼軍，孫莉英，李君蘭，等. 緩坡面細(xì)溝發(fā)育過程及水沙關(guān)系的室內(nèi)實(shí)驗(yàn)研究[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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] 王龍生，蔡強(qiáng)國，蔡崇法，等. 黃土坡面細(xì)溝形態(tài)變化及其與流速之間的關(guān)系[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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]. 土壤學(xué)報(bào)，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] 周佩華，王占禮. 黃土高原土壤侵蝕暴雨標(biāo)準(zhǔn)[J]. 水土保持通報(bào)，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ù)的細(xì)溝與細(xì)溝水流參數(shù)測量及特征[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，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] 和繼軍，呂燁，宮輝力，等. 細(xì)溝侵蝕特征及其產(chǎn)流產(chǎn)沙過程試驗(yàn)研究[J]. 水利學(xué)報(bào)，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]. 土壤學(xué)報(bào)，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] 覃超，吳紅艷，鄭粉莉，等. 黃土坡面細(xì)溝侵蝕及水動(dòng)力學(xué)參數(shù)的時(shí)空變化特征[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，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)