林草和梯田措施對小流域降雨—徑流—輸沙過程的影響分析

白 娟，張 亦 弛，楊 勝 天 ，劉 曉 燕，甘 甫 平，閆 柏 琨，郭 藝

(1.中國自然資源航空物探遙感中心，北京 100083；2.北京師范大學地理科學學部，北京100875；3.北京師范大學水科學研究院，北京100875；4.水利部黃河水利委員會，河南 鄭州 450004)

0 引言

水土流失是當前世界面臨的嚴重環境問題之一[1-4]，如何有效控制水土流失的發生和發展是流域管理的重點問題。坡改梯和林草種植是坡面易侵蝕區常用的水土保持措施，不僅可削減自身產水產沙，還可攔截上方坡面的來水來沙，對坡面乃至流域尺度的水沙過程具有顯著影響[5-9]。但目前對流域尺度林草、梯田措施減水減沙作用的研究仍有待完善：1)野外實驗觀測方法有助于揭示坡面尺度林草、梯田措施減水減沙的過程和機理[4,10-13]，但局限于小尺度研究[14]，受下墊面空間異質性影響較大[15]；2)水保法僅考慮林草、梯田措施自身的減水減沙作用，忽略了其對上方坡面來水來沙的削減，難以代表流域全局范圍內的水沙過程；3)水文法可得到下墊面變化導致減水減沙的整體效果，但無法給出單項措施的減水減沙貢獻，且計算結果受建模期雨量數據的代表性影響[16]；4)水文模型法側重于模擬坡面措施自身的減水減沙作用，忽視了對林草、梯田措施在匯流過程中攔水攔沙作用的參數化表達，且計算單一措施的減水減沙效率時，較難區分其他措施的影響，計算不同措施耦合的減水減沙效率時，會產生重復計算。

黃河中游位于黃土高原，水土流失嚴重，直接影響黃河的生態安全和區域社會經濟可持續發展[17-22]。為控制該區水土流失，國家實施了大規模的退耕還林還草和梯田建設工程，定量評價林草、梯田的減水減沙作用符合流域水沙調控的發展要求，對黃河流域水土流失防治和生態保護具有重要意義。Luo等應用修正通用土壤流失方程(Modified Universal Soil Loss Equation，MUSLE)，耦合分布式場次暴雨徑流LCM模型，構建了分布式LCM-MUSLE坡面水沙聯動模型，并對MUSLE方程中的作物管理因子(C)、地形因子(LS)和保持措施因子(P)進行了修訂，實現其在黃土高原應用的本土化，提高了產沙模擬精度[23]，但該模型忽略了匯流過程中坡面措施減水減沙作用；Bai 等對該模型的匯流方法進行改進，將等流時面內林草、梯田控制區域視為最小減水減沙單元，構建“控制區域—等流時面—子流域”一體化坡面匯流系統，考慮林草和梯田措施在匯流過程中對徑流和輸沙的削減作用，詮釋了水保措施導致的下墊面變化對流域水沙的影響[24]。在此基礎上，本研究以黃河中游多沙粗沙區為研究對象，分別應用原分布式LCM-MUSLE坡面水沙聯動模型和匯流方法改進后模型進行場次“降雨—徑流—輸沙”過程模擬，定量評估林草、梯田及其耦合措施在流域坡面匯流過程中減水減沙的貢獻率，并討論其主要影響因素，為黃河中游林草植被和梯田建設對流域水沙影響的定量評估提供科學依據。

1 研究區與數據

1.1 研究區概況

黃河中游多沙粗沙區暴雨多、洪水大、含沙量高、泥沙顆粒粗，是黃土高原水土保持治理的重點區域[25-27]。為控制水土流失，該區開展了造林、種草、封禁、修梯田等坡面措施以及修淤地壩等溝道措施[16]，致使黃河中游下墊面特征發生顯著變化[25]。自2000年黃河來沙量明顯減少，2000-2019年潼關來沙量為2.45億t/a，2000-2009年和2010-2019年的下墊面產沙能力相比1919-1959年均值分別下降了57%和80%[28]。綜合考慮研究區代表性、已有數據積累以及林草、梯田面積和空間分布的巨大變化等因素，本文選取偏關河偏關水文站以上流域和清澗河子長水文站以上流域為研究區(圖1)。

圖1 研究區位置Fig.1 Location of the study area

偏關河和清澗河均為黃河一級支流，均位于黃土高原丘陵溝壑區第1副區。偏關河流域面積為2 089.36 km2，偏關水文站控制面積1 922.81 km2，海拔984～2 162 m，屬溫帶半干旱氣候，多年平均降水量為429 mm，5-9月降水量占全年降水量的80.9%，產流類型為超滲產流，多年平均徑流量為3 948萬m3，多年平均輸沙量為1 258萬t，洪水輸沙模數為6 523 t/(km2·a)，水土流失嚴重[29]。清澗河流域面積約4 080 km2，其子長水文站控制面積約889.65 km2，海拔1 047～1 585 m，屬暖溫帶半干旱氣候，流域多年平均降水量為482.4 mm，汛期6-9月降水量占年降水量的70%以上，產流類型為超滲產流，多年平均徑流量為13 700萬m3，多年平均輸沙量為3 300萬t，汛期輸沙量占90%以上[30]。

1.2 數據來源

研究數據包括：1)流域內雨量站觀測的逐小時降雨數據以及水文卡口站觀測的徑流量和含沙量數據，對其進行時間重采樣和空間插值[31]；2)流域泥沙粒徑數據(D50)，由水利部門公布的各水文站粒徑數據經反距離加權插值后獲取；3)ASTER GDEM 數據，用于提取地形指數、等流時線、河網和子流域等空間信息；4)土地利用類型和植被覆蓋度數據，基于Landsat影像獲取，前者參照《土地利用現狀分類(GB/T21010-2007)》標準，后者基于歸一化植被指數(NDVI)，采用像元二分模型提取得到，研究中分別用1978年、2010年的土地利用和植被覆蓋度數據代表1980s和2010s的土地利用和植被覆蓋度狀況；5)梯田數據，采用水利部黃河水利委員會(YRCC)基于2012年ZY-3衛星影像解譯的結果；6)土壤數據，包括土壤類型、土壤機械組成、土壤有機碳含量和土壤飽和含水量，土壤類型來源于中國1∶100萬土壤圖，土壤機械組成和有機碳含量等屬性通過查詢世界土壤數據庫HWSD (Harmonized World Soil Database version 1.1)獲取，并運用SPAW(Soil-Plant-Atmosphere- Water Field & Pond Hydrology)模型計算各土壤類型的飽和含水量。

2 研究方法

2.1 改進的LCM-MUSLE坡面水沙聯動模型

圖2 改進的LCM-MUSLE模型結構Fig.2 Framework of modified LCM-MUSLE model

(1)

(2)

Vui=ΔAVcontrol,i/ΔAi

(3)

(4)

ΔTi,j=Ti,j-Ti,j-1

(5)

(6)

Tui=ΔATcontrol,i/ΔAi

(7)

Tci=ΔAterrace,i×di

(8)

采用等流時線法計算各子流域坡面匯流后，利用馬斯京根方法計算河道匯流[24]，匯沙過程還需考慮子流域泥沙顆粒大小和流域水力學特性。通過等流時線法計算各子流域的產沙量后，根據匯沙模型[36]計算河道匯沙過程前子流域i的輸沙量SYi：

(9)

式中：Yi為等流時線法坡面匯沙后計算出的子流域產沙量(t)；B為匯沙系數；Ti為子流域到流域出口的匯流時間(h)；D50i為子流域泥沙粒徑中值(mm)。

2.2 模型精度評價

基于小時尺度徑流量、輸沙量觀測數據，運用納西效率系數(Nash)對模型模擬結果進行精度評價，并分析其在研究區的適用性，計算公式為：

(10)

2.3 林草、梯田措施減水減沙作用評估

林草、梯田措施的減水減沙作用是指在相同降雨條件下實施措施后流域較天然時期減少的產水產沙量[37]。將流域植被基本穩定、水保措施相對較少、人類活動影響較小的1980s作為天然時期，將流域植被恢復明顯、水保措施建設密集、人類活動影響顯著的2010s作為現狀年。根據研究區在不同時期有無林草、梯田措施，提出O1(原模型)、R1(林草措施)、R2(梯田措施)、R3(林草+梯田措施)4種模擬情景。由于1980s的梯田數據不完整，且數據質量較差，故對于1980s的場次降雨僅模擬O1和R1情景。

通過構建削減率(RR)和單位面積削減量(AR)兩個指標，定量分析林草、梯田在匯流過程中的減水減沙作用。以徑流為例，將場次小時尺度徑流模擬值加總，獲取場次徑流總量，然后統計場次降雨在4種情景下的模擬徑流總量，分別計算R1、R2和R3相對O1的徑流削減率(式(11))，以及R1和R2的單位面積徑流削減量(式(12))；同理可計算R1、R2和R3相對O1的輸沙削減率，以及R1和R2的單位面積輸沙削減量。

RRi=(y-yi)/y×100%

(11)

ARi=(y-yi)/Ai

(12)

式中：RRi、ARi分別為流域場次降雨在i(i=1,2,3)措施下的徑流削減率(%)和單位面積徑流削減量；y為O1情景的模擬徑流總量(m3)；yi為R1、R2 或 R3情景的模擬徑流總量(m3)；A為流域林草或梯田的面積(km2)。

3 結果分析

3.1 場次“降雨—徑流—輸沙”模擬與精度驗證

經篩選，1981-2010年兩流域均有8個場次“降雨—徑流—輸沙”觀測數據用于模型模擬和精度驗證，每場次降雨以其起始時間(年/日/時)進行編號。偏關站以上流域以1981/203/17、1983/215/22和1983/235/16共3個場次作為參數率定期，1988/199/13、1989/203/19、2006/195/05、2006/224/08和2010/263/20共5個場次作為驗證期；子長站以上流域以1986/187/14、1987/238/00和1988/218/14共3個場次作為參數率定期，1988/197/04、1988/237/18、2006/188/20、2006/237/02和2006/263/20共5個場次作為驗證期。

由兩流域徑流量(圖3)和輸沙量(圖4)模擬值與觀測值對比可知：總體上，對于1980s的各場次降雨，R1的峰值流量和徑流總量比O1均有明顯削減，對于2010s的各場次降雨，徑流量峰值排序為O1>R2>R1>R3，表明林草措施對峰值流量和徑流總量的削減作用大于梯田措施，綜合考慮兩種措施對峰值流量和徑流總量的削減作用最強。對于場次輸沙量而言，2010s O1與R1兩種情景下輸沙量的差異小于二者徑流量的差異，這是由于在MUSLE產沙模型中，Rs、C和P因子的設定已體現對產沙過程的削減作用。對于1980s的各場次降雨，R1的輸沙量峰值均低于O1的峰值；對于2010s的各場次降雨，輸沙量峰值排序為O1>R2>R1>R3，表明林草措施對輸沙量峰值削減幅度比梯田措施大，兩種措施耦合對輸沙量的削減作用最強。

圖3 偏關站以上流域和子長站以上流域徑流量模擬值與觀測值對比Fig.3 Comparison of observed and simulated runoff in upstream watershed of Pianguan station and upstream watershed of Zichang station

圖4 偏關站以上流域和子長站以上流域輸沙量模擬值與觀測值對比Fig.4 Comparison of observed and simulated sediment discharge in upstream watershed of Pianguan station and upstream watershed of Zichang station

由兩流域場次“降雨—徑流—輸沙”模擬精度(表1)可知，偏關站以上流域在1980s和2010s的場次徑流量平均Nash系數分別從O情景的0.46、-15.29 提高到R3情景的0.62、0.39，場次輸沙量平均Nash系數分別從O情景的0.08、-2.47提高到R情景的0.22、0.37；子長站以上流域在1980s和2010s的場次徑流量平均Nash系數分別從O情景的-0.86、-5.28提高到R3情景的0.11、0.43，場次輸沙量平均Nash系數分別從O情景的-1.40、-0.34提高到R3情景的-0.24、0.34，表明原模型在現狀年模擬精度不高，改進后模型適用于現狀年的徑流(輸沙)量模擬。

表1 兩流域場次徑流量/輸沙量平均Nash系數Table 1 Average Nash coefficients of runoff and sediment discharge in upstream watershed of Pianguan station and upstream watershed of Zichang station

3.2 林草、梯田措施減水減沙作用分析

統計1980-2010年各場次降雨在不同情景下的徑流(輸沙)量，計算林草措施、梯田措施及二者耦合的徑流(輸沙)削減率，以及林草、梯田措施的單位面積徑流(輸沙)削減量(圖5、圖6)。可以看出：

圖5 林草、梯田措施徑流削減作用Fig.5 Runoff reduction effects of vegetation and terraces

圖6 林草、梯田措施輸沙削減作用Fig.6 Sediment reduction effects of vegetation and terraces

對于偏關站以上流域而言，1980s場次降雨林草措施的平均徑流削減率為17.86%，平均單位面積徑流削減量為1 378.84 m3/km2；2010s場次降雨林草措施的平均徑流削減率RR1為48.08%，受植被覆蓋度增加影響，比1980s增加了30.22%，平均單位面積徑流削減量為4 063.59 m3/km2；梯田措施的平均徑流削減率RR2為26.64%(小于林草措施)，平均單位面積徑流削減量為15 827.36 m3/km2(遠大于林草措施)；林草梯田耦合措施的平均徑流削減率RR3為69.26%。1980s場次降雨林草措施的平均輸沙削減率為12.02%，平均單位面積輸沙削減量為341.19 t/km2；2010s場次降雨林草措施的平均輸沙削減率RR1為32.59%，受植被覆蓋度增加影響，比1980s增加了20.57%，平均單位面積輸沙削減量為663.61 t/km2；梯田措施的平均輸沙削減率RR2為24.50%(小于林草措施)，平均單位面積輸沙削減量為3 384.71 t/km2(遠大于林草措施)；林草梯田耦合措施的平均輸沙削減率RR3為54.18%。

對于子長站以上流域而言，1980s場次降雨林草措施的平均徑流削減率為21.69%，平均單位面積徑流削減量為3 867.98 m3/km2；2010s場次降雨林草措施的平均徑流削減率RR1為61.01%，受植被覆蓋度增加影響，比1980s增加了39.32%，平均單位面積徑流削減量為7 594.16 m3/km2；梯田措施的平均徑流削減率RR2為8.55%(小于林草措施)，平均單位面積徑流削減量為21 638.68 m3/km2(遠大于林草措施)；林草梯田耦合措施的平均徑流削減率RR3為65.54%。1980s場次降雨林草措施的平均輸沙削減率為19.46%，平均單位面積輸沙削減量為1 423.69 t/km2；2010s場次降雨林草措施的平均輸沙削減率RR1為43.71%，受植被覆蓋度增加影響，比1980s增加了24.25%，平均單位面積輸沙削減量為541.07 t/km2；梯田措施的平均輸沙削減率RR2為13.94%(小于林草措施)，單位面積輸沙削減量為1 525.85 t/km2(遠大于林草措施)；林草梯田耦合措施的平均輸沙削減率RR3為55.84%。

總體上，兩流域徑流(輸沙)削減率排序均為林草+梯田措施>林草措施>梯田措施，且林草和梯田措施的徑流(輸沙)削減率之和略大于二者耦合結果，可能是當梯田起到減水作用后，同一等流時面內植被的徑流削減率會降低。此外，本文假設偏關站以上流域均為一類梯田，子長站以上流域梯田均無田埂，可能會影響梯田措施減水作用的估算結果。

3.3 林草、梯田減水減沙影響因素分析

流域林草植被面積和植被覆蓋度變化直接影響流域場次徑流量、輸沙量，林草植被的空間分布也會對匯流過程中林草的減水減沙作用產生影響。為分析流域場次降雨徑流(輸沙)量變化與林草植被數量及其空間分布的關系，分別計算林草面積比、植被覆蓋度和林草減水面積比3個指標與各場次降雨的林草徑流削減率、輸沙削減率的Pearson相關系數。其中，林草面積比為流域內林草地面積占易侵蝕面積比例，林草減水面積比為流域內林草減水面積(即林草面積及其集水區面積之和)占易侵蝕面積比例。結果顯示：在α=0.05的顯著性水平下，林草徑流削減率與林草面積比、植被覆蓋度和林草減水面積比的相關系數分別為0.81、0.80和0.34，林草輸沙削減率與上述3個指標的相關系數分別為0.76、0.71和0.48，說明流域林草徑流(輸沙)削減率與植被數量明顯相關，與植被空間分布的相關性較弱。

流域梯田面積和梯田田埂完整度會直接影響流域場次降雨徑流(輸沙)量，梯田的空間分布(用梯田集水區與流域面積比間接表征)也會對匯流過程中梯田的減水減沙作用產生影響。為分析流域場次降雨徑流(輸沙)量變化與梯田集水量及梯田空間分布的關系，分別計算梯田比、田埂完整度、梯田集水區與流域面積比3個指標與各場次降雨的梯田徑流(輸沙)削減率的Pearson相關系數。其中，梯田比為某地區水平梯田面積占區內輕度以上水土流失面積比例，該指標比梯田面積與流域面積比更能反映梯田對流域主要產沙區的控制程度[37,38]；研究區梯田均為水平梯田，故田埂完整度可反映梯田集水量。結果顯示：在α=0.05的顯著性水平下，梯田徑流削減率與梯田比、田埂完整度、梯田集水區與流域面積比的相關系數分別為0.98、0.56和0.94，梯田輸沙削減率與上述3個指標的相關系數分別為0.98、0.44和0.97，說明流域梯田徑流(輸沙)削減率與梯田面積和梯田空間分布明顯相關，與梯田田埂完整度的相關性較弱。

4 結論

本文分別采用改進的分布式LCM-MUSLE坡面水沙聯動模型和原模型在黃河中游偏關站以上流域和子長站以上流域開展場次降雨水沙過程模擬，計算林草、梯田及二者耦合措施的減水減沙作用并討論其主要影響因素。結論如下：1)改進后的分布式LCM-MUSLE坡面水沙聯動模型更適用于現狀年的徑流(輸沙)量模擬，可反映下墊面顯著變化(坡改梯和植被恢復)對徑流(輸沙)的影響；2)單位面積梯田的減水減沙量遠大于林草，林草措施總的徑流(輸沙)削減率顯著大于梯田措施；3)流域尺度梯田措施下徑流(輸沙)削減率影響因素排序為梯田比>梯田集水區與流域面積比>田埂完整度，林草措施下影響因素排序為林草面積比>植被覆蓋度>林草減水面積比。

本文一定程度上剖析了近年來黃河中游來水來沙量銳減的原因，建立了徑流量、輸沙量對于梯田、林草變化的響應機制，但仍存在如下不足：文中田埂完整度的設定來源于已發表文獻，與田埂真實完整程度可能存在偏差，造成徑流、泥沙削減率略高于真實值；未考慮水庫、淤地壩等工程措施以及開礦、修路等土地利用開發對河道及河床周圍水文過程的影響。未來將開展下墊面變化對河道及河床周圍水文過程的影響研究，以期更全面地刻畫多沙粗沙區“降雨—徑流—輸沙”過程。