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

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

0 引言

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

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

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

1.1 研究區(qū)概況

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

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

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

1.2 數(shù)據(jù)來源

研究數(shù)據(jù)包括：1)流域內(nèi)雨量站觀測的逐小時降雨數(shù)據(jù)以及水文卡口站觀測的徑流量和含沙量數(shù)據(jù)，對其進(jìn)行時間重采樣和空間插值[31]；2)流域泥沙粒徑數(shù)據(jù)(D50)，由水利部門公布的各水文站粒徑數(shù)據(jù)經(jīng)反距離加權(quán)插值后獲取；3)ASTER GDEM 數(shù)據(jù)，用于提取地形指數(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ù)，采用水利部黃河水利委員會(YRCC)基于2012年ZY-3衛(wèi)星影像解譯的結(jié)果；6)土壤數(shù)據(jù)，包括土壤類型、土壤機(jī)械組成、土壤有機(jī)碳含量和土壤飽和含水量，土壤類型來源于中國1∶100萬土壤圖，土壤機(jī)械組成和有機(jī)碳含量等屬性通過查詢世界土壤數(shù)據(jù)庫HWSD (Harmonized World Soil Database version 1.1)獲取，并運(yùn)用SPAW(Soil-Plant-Atmosphere- Water Field & Pond Hydrology)模型計算各土壤類型的飽和含水量。

2 研究方法

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

圖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)

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

(9)

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

2.2 模型精度評價

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

(10)

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

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

通過構(gòu)建削減率(RR)和單位面積削減量(AR)兩個指標(biāo)，定量分析林草、梯田在匯流過程中的減水減沙作用。以徑流為例，將場次小時尺度徑流模擬值加總，獲取場次徑流總量，然后統(tǒng)計場次降雨在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 結(jié)果分析

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

經(jīng)篩選，1981-2010年兩流域均有8個場次“降雨—徑流—輸沙”觀測數(shù)據(jù)用于模型模擬和精度驗證，每場次降雨以其起始時間(年/日/時)進(jìn)行編號。偏關(guān)站以上流域以1981/203/17、1983/215/22和1983/235/16共3個場次作為參數(shù)率定期，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個場次作為參數(shù)率定期，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，表明林草措施對峰值流量和徑流總量的削減作用大于梯田措施，綜合考慮兩種措施對峰值流量和徑流總量的削減作用最強(qiáng)。對于場次輸沙量而言，2010s O1與R1兩種情景下輸沙量的差異小于二者徑流量的差異，這是由于在MUSLE產(chǎn)沙模型中，Rs、C和P因子的設(shè)定已體現(xiàn)對產(chǎn)沙過程的削減作用。對于1980s的各場次降雨，R1的輸沙量峰值均低于O1的峰值；對于2010s的各場次降雨，輸沙量峰值排序為O1>R2>R1>R3，表明林草措施對輸沙量峰值削減幅度比梯田措施大，兩種措施耦合對輸沙量的削減作用最強(qiáng)。

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

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

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

表1 兩流域場次徑流量/輸沙量平均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)計1980-2010年各場次降雨在不同情景下的徑流(輸沙)量，計算林草措施、梯田措施及二者耦合的徑流(輸沙)削減率，以及林草、梯田措施的單位面積徑流(輸沙)削減量(圖5、圖6)?？梢钥闯觯?/p>

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

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

對于偏關(guān)站以上流域而言，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(遠(yuǎn)大于林草措施)；林草梯田耦合措施的平均徑流削減率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(遠(yuǎn)大于林草措施)；林草梯田耦合措施的平均輸沙削減率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(遠(yuǎn)大于林草措施)；林草梯田耦合措施的平均徑流削減率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(遠(yuǎn)大于林草措施)；林草梯田耦合措施的平均輸沙削減率RR3為55.84%。

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

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

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

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

4 結(jié)論

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

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