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

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

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

0 引言

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

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

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

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

圖1 研究區(qū)位置Fig.1 Location of the study area

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

1.2 數(shù)據(jù)來(lái)源

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

2 研究方法

2.1 改進(jìn)的LCM-MUSLE坡面水沙聯(lián)動(dòng)模型

圖2 改進(jìn)的LCM-MUSLE模型結(jié)構(gòu)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)

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

(9)

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

2.2 模型精度評(píng)價(jià)

基于小時(shí)尺度徑流量、輸沙量觀測(cè)數(shù)據(jù)，運(yùn)用納西效率系數(shù)(Nash)對(duì)模型模擬結(jié)果進(jìn)行精度評(píng)價(jià)，并分析其在研究區(qū)的適用性，計(jì)算公式為：

(10)

2.3 林草、梯田措施減水減沙作用評(píng)估

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

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

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

(11)

ARi=(y-yi)/Ai

(12)

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

3 結(jié)果分析

3.1 場(chǎng)次“降雨—徑流—輸沙”模擬與精度驗(yàn)證

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

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

圖3 偏關(guān)站以上流域和子長(zhǎng)站以上流域徑流量模擬值與觀測(cè)值對(duì)比Fig.3 Comparison of observed and simulated runoff in upstream watershed of Pianguan station and upstream watershed of Zichang station

圖4 偏關(guān)站以上流域和子長(zhǎng)站以上流域輸沙量模擬值與觀測(cè)值對(duì)比Fig.4 Comparison of observed and simulated sediment discharge in upstream watershed of Pianguan station and upstream watershed of Zichang station

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

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

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

統(tǒng)計(jì)1980-2010年各場(chǎng)次降雨在不同情景下的徑流(輸沙)量，計(jì)算林草措施、梯田措施及二者耦合的徑流(輸沙)削減率，以及林草、梯田措施的單位面積徑流(輸沙)削減量(圖5、圖6)?？梢钥闯觯?/p>

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

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

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

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

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

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

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

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

4 結(jié)論

本文分別采用改進(jìn)的分布式LCM-MUSLE坡面水沙聯(lián)動(dòng)模型和原模型在黃河中游偏關(guān)站以上流域和子長(zhǎng)站以上流域開展場(chǎng)次降雨水沙過(guò)程模擬，計(jì)算林草、梯田及二者耦合措施的減水減沙作用并討論其主要影響因素。結(jié)論如下：1)改進(jìn)后的分布式LCM-MUSLE坡面水沙聯(lián)動(dòng)模型更適用于現(xiàn)狀年的徑流(輸沙)量模擬，可反映下墊面顯著變化(坡改梯和植被恢復(fù))對(duì)徑流(輸沙)的影響；2)單位面積梯田的減水減沙量遠(yuǎn)大于林草，林草措施總的徑流(輸沙)削減率顯著大于梯田措施；3)流域尺度梯田措施下徑流(輸沙)削減率影響因素排序?yàn)樘萏锉?梯田集水區(qū)與流域面積比>田埂完整度，林草措施下影響因素排序?yàn)榱植菝娣e比>植被覆蓋度>林草減水面積比。

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