河南省豫西黃土丘陵區(qū)2015—2020年水土保持功能時(shí)空變化

杜 軍， 趙勝朝， 徐建昭， 郭 雷

（1.河南省科學(xué)院地理研究所，鄭州 450052； 2.河南省水土保持監(jiān)測(cè)總站，鄭州 450003）

習(xí)近平總書記在黃河流域生態(tài)保護(hù)和高質(zhì)量發(fā)展座談會(huì)上發(fā)表重要講話中指出，黃河流域是我國(guó)重要的生態(tài)屏障和經(jīng)濟(jì)地帶，是打贏脫貧攻堅(jiān)戰(zhàn)的重要區(qū)域，在我國(guó)經(jīng)濟(jì)社會(huì)發(fā)展和生態(tài)安全方面具有十分重要的地位，保護(hù)黃河是事關(guān)中華民族偉大復(fù)興的千秋大計(jì). 豫西黃土丘陵區(qū)是河南省沿黃區(qū)域生態(tài)環(huán)境比較脆弱的區(qū)域之一，也是河南省水土流失最嚴(yán)重的地區(qū)之一［1-2］，區(qū)域水土保持功能的改善與否，對(duì)于黃河流域生態(tài)保護(hù)與高質(zhì)量發(fā)展具有重要的影響［3-5］. 經(jīng)過近些年國(guó)土綠化、封山育林、坡改梯、淤地壩等水保措施［6-8］，區(qū)域生態(tài)環(huán)境有了較大改善，但對(duì)于區(qū)域水土保持功能空間變化以及近五年的時(shí)間變化缺乏相關(guān)研究，因此本研究在前人研究的基礎(chǔ)上，結(jié)合豫西黃土丘陵的實(shí)際情況，基于2015和2020年的遙感影像獲取研究區(qū)土地利用動(dòng)態(tài)變化數(shù)據(jù)，采用RUSLE 模型對(duì)豫西黃土丘陵地區(qū)2015—2020 年間區(qū)域水土保持功能的變化特征進(jìn)行評(píng)估，為豫西黃土丘陵區(qū)土壤侵蝕的治理提供一定的參考依據(jù).

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究范圍

豫西黃土丘陵區(qū)［9-10］（110°21′～113°30′、33°32′～35°17′）位于河南省西北部，行政區(qū)劃上包括鄭州市、洛陽(yáng)市、焦作市、三門峽市和濟(jì)源市的19 個(gè)縣（市），總面積2.71 萬(wàn)km2，約占河南省總面積的16.23%（圖1）. 屬于黃土高原向黃淮海平原過渡的交錯(cuò)地帶，水土流失較為嚴(yán)重，是河南省生態(tài)環(huán)境比較脆弱的區(qū)域之一. 區(qū)域內(nèi)地形起伏較大、地貌復(fù)雜、土地利用類型多樣，最低海拔85 m，最高海拔2392 m，地貌類型主要包括山地（中山和低山）、丘陵和平原等，土地利用類型主要以耕地、林地為主. 該區(qū)屬于溫帶半濕潤(rùn)半干旱氣候，多年平均氣溫8～15 ℃，多年平均降水約為574.0～634.4 mm.

圖1 豫西黃土丘陵區(qū)概況Fig.1 General situation of Loess Hilly Area in western Henan

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

本研究中利用2015年、2020年的遙感影像數(shù)據(jù)，對(duì)影像進(jìn)行融合、校正等預(yù)處理，利用監(jiān)督分類法并結(jié)合人工目視判讀對(duì)2期利用進(jìn)行解譯［11］，將研究區(qū)分為林地、草地、耕地、水域、建設(shè)用地和未利用地6種類型，解譯結(jié)果隨機(jī)采樣并進(jìn)行野外核查，總體精度大于95%.

2 研究方法

利用RUSLE模型［12-14］的土壤侵蝕評(píng)價(jià)方法，土壤侵蝕模型RUSLE考慮了降雨、土壤可蝕性、作物管理、坡度坡長(zhǎng)和水土保持措施五大因子，方程式如下：

式中：A為年平均土壤流失量，t/（km2·a）；R為降雨和徑流侵蝕因子，MJ·mm/（km2·h·a）；K為土壤可蝕性因子，此因子是指在標(biāo)準(zhǔn)狀態(tài)下（小區(qū)長(zhǎng)22.1 m，坡度9%，休閑地且順坡耕作）每年單位降雨侵蝕指標(biāo)所產(chǎn)生的侵蝕數(shù)量，單位為t·h/（MJ·mm）；L，S為地形因子，其中L為坡長(zhǎng)因子，S為坡度因子；C為作物管理因子；P為治理措施因子.

3 土壤侵蝕評(píng)價(jià)因子提取方法

3.1 降雨侵蝕力R因子

降水時(shí)雨滴擊濺、分離土壤顆粒以及徑流沖刷和轉(zhuǎn)運(yùn)是引起土壤侵蝕的主要?jiǎng)恿σ蛩? 美國(guó)在使用通用土壤流失方程預(yù)報(bào)土壤流失量時(shí)，一直用EI30作為經(jīng)典算法來計(jì)算降雨侵蝕力指標(biāo)，這種經(jīng)典算法最初是由Wischmeier 為確定降雨侵蝕力指標(biāo)［15］，對(duì)降雨量、降雨動(dòng)能、最大15 min 雨強(qiáng)、最大30 min 雨強(qiáng)等變量與土壤流失量之間進(jìn)行決定系數(shù)分析發(fā)現(xiàn)的，EI30 可以解釋土壤流失量變化的70.2%～89.2%. 經(jīng)典算法一般都需要降雨過程資料，必須用斷點(diǎn)降雨資料計(jì)算降雨動(dòng)能，工作量大且資料的獲取難度大，因此采用經(jīng)典算法計(jì)算降雨侵蝕力十分不易. 后來有學(xué)者逐步使用了簡(jiǎn)易算法，即采用年、月或日降雨資料進(jìn)行降雨侵蝕力的計(jì)算. 簡(jiǎn)易算法資料獲取更為容易，同時(shí)資料的整理及計(jì)算相對(duì)簡(jiǎn)單. 謝云等［16］對(duì)我國(guó)降雨情況進(jìn)行分析得出日降雨一般能夠替代次降雨計(jì)算降雨侵蝕力. 章文波和付金生［17］通過全國(guó)71 個(gè)氣象站點(diǎn)的降雨數(shù)據(jù)分析得出，在利用年、月、日降雨資料等5 種代表性雨量資料估算降雨侵蝕力的結(jié)果中，日雨量計(jì)算多年平均侵蝕力的精度最高，并建立了利用日雨量估算降雨侵蝕力的簡(jiǎn)易算法模型，提出了模型的參數(shù)估計(jì)方法. 模型在全國(guó)的決定系數(shù)平均為0.718，估算多年平均侵蝕力的相對(duì)誤差平均為4.2%.公式如下：

式中：Mi為第i個(gè)半月時(shí)段的侵蝕力值；α，β為模型參數(shù)；K為該半月時(shí)段內(nèi)的天數(shù)；Dj為半月時(shí)段內(nèi)第j天的日雨量，要求日雨量≥12 mm，否則以0計(jì)算，12 mm與中國(guó)侵蝕性降雨標(biāo)準(zhǔn)對(duì)應(yīng). 參數(shù)α和β反映了區(qū)域降雨特征，如下式：

式中：Pd12為日雨量≥12 mm的日平均雨量；Py12為日雨量≥12 mm的年平均雨量.

由于不是每次降雨都會(huì)產(chǎn)生侵蝕，少于一定數(shù)量的降雨經(jīng)常不會(huì)引起土壤侵蝕，因此在統(tǒng)計(jì)降雨資料時(shí)，存在一個(gè)侵蝕性降雨標(biāo)準(zhǔn). 美國(guó)侵蝕性降雨標(biāo)準(zhǔn)定為次降雨量大于6.5 mm，或一次降雨時(shí)間長(zhǎng)于12.7 min 或15 min. 中國(guó)在黃土高原的多年研究中，建立了相應(yīng)的侵蝕性降雨雨量標(biāo)準(zhǔn)為日降雨量大于12 mm，也是本文采用的標(biāo)準(zhǔn).

本研究利用章文波等的日降雨量估算方法計(jì)算降雨侵蝕力R因子. 降雨數(shù)據(jù)來自河南省2015—2020年的逐日降雨資料，并整理其中的侵蝕性降雨，按照以上公式進(jìn)行估算.

3.2 土壤可蝕性K因子

土壤可蝕性是土壤侵蝕預(yù)報(bào)中的一個(gè)重要方面，它是評(píng)價(jià)土壤是否易受侵蝕營(yíng)力破壞的性能，也是土壤對(duì)侵蝕營(yíng)力分離和搬運(yùn)作用的敏感性的體現(xiàn). 由于自然界存在的土壤類型多種多樣，在確定其可蝕性值大小時(shí)，不論人力物力還是時(shí)間上都很難現(xiàn)實(shí). 因此，國(guó)內(nèi)外的很多研究都致力于土壤可蝕性的估算.Wischmeier 和Smith［18］根據(jù)實(shí)測(cè)的土壤可蝕性因子與土壤性質(zhì)的相關(guān)性，建立了土壤可蝕性K值的諾謨方程，參數(shù)包括土壤質(zhì)地、土壤有機(jī)質(zhì)、土壤結(jié)構(gòu)和土壤通透性等. 但是這種方法需要較多參數(shù)，特別是土壤結(jié)構(gòu)級(jí)別和土壤滲透級(jí)別很難準(zhǔn)確地獲得，因此較為煩瑣. Williams 和Sharply［19］在侵蝕-生產(chǎn)力評(píng)價(jià)模型EPIC（Erosion productivity impact calculator）中，把土壤可蝕性因子K值的計(jì)算公式簡(jiǎn)化成只與土壤砂粒含量、粉粒含量、黏粒含量和土壤有機(jī)碳含量有關(guān)的公式，建立了一種相對(duì)簡(jiǎn)便的土壤可蝕性因子K值的計(jì)算方法. 我國(guó)學(xué)者也對(duì)土壤可蝕性K值進(jìn)行了大量研究，涉及研究區(qū)土壤K值的修正諾謨公式、修正EPIC模型以及只考慮幾何平均粒徑的計(jì)算方法等.

土壤數(shù)據(jù)（質(zhì)地、有機(jī)質(zhì)等）來源于黑河計(jì)劃數(shù)據(jù)管理中心（http：//westdc.westgis.ac.cn）的“基于世界土壤數(shù)據(jù)庫(kù)的中國(guó)土壤數(shù)據(jù)集（v1.1）”，從中截取了河南省地區(qū)數(shù)據(jù).

式中：SANd為砂粒含量，%；SILi為粉粒含量，%；CLAi為黏粒含量，%；C為有機(jī)質(zhì)含量，%；SN1=1-SAN/100.

3.3 地形因子LS因子

地形是自然地理要素中最基本的要素，對(duì)土壤侵蝕強(qiáng)度有重大影響. 各地形因子中，坡度、坡長(zhǎng)應(yīng)用最為廣泛，因?yàn)閺膭?dòng)力學(xué)角度分析，所有的地形因子的作用最終都取決于地形坡度. 因此在土壤侵蝕的研究中，使用坡度、坡長(zhǎng)相乘來計(jì)算地形因子LS值. 坡度和坡長(zhǎng)因子是降雨侵蝕作用中的加速因子，反映了地形對(duì)土壤侵蝕的影響，決定著徑流的運(yùn)動(dòng)狀態(tài)和方向.

數(shù)字高程模型（Digital Elevation Model，DEM）是描述地表每個(gè)單元的高程屬性和空間位置分布的有序集合. DEM數(shù)據(jù)的采集方式主要有地面測(cè)量、攝影測(cè)量、激光掃描、地形圖手扶跟蹤數(shù)字化、立體遙感、GPS、地形圖屏幕數(shù)字化、干涉雷達(dá)等. 其常用的表達(dá)形式有規(guī)則格網(wǎng)和不規(guī)則三角網(wǎng)兩種. 本文選取空間分辨率為30 m的DEM數(shù)據(jù).

3.3.1 坡度因子S坡度是影響土壤侵蝕的另一個(gè)重要參數(shù). 在土壤侵蝕量與坡度的關(guān)系上，一般建立的是冪函數(shù)或二次多項(xiàng)式關(guān)系. McCool等［20］通過分析大量的野外徑流試驗(yàn)小區(qū)資料和人工降雨資料，分別在0～9%和9%～25%坡度范圍內(nèi)建立了土壤侵蝕量與坡度的線性關(guān)系. Liu等［21］在McCool等的基礎(chǔ)上，結(jié)合中國(guó)黃土高原綏德、安塞和天水三個(gè)水土保持試驗(yàn)站的試驗(yàn)資料，建立了25%以上的陡坡坡度和土壤侵蝕量的關(guān)系式. 本文采用Liu等和McCool等結(jié)合的方法，計(jì)算坡度因子.

式中：θ為像元坡度，%.

3.3.2 坡長(zhǎng)因子L坡長(zhǎng)L因子是反映地形地貌特征對(duì)土壤侵蝕的影響，是水土保持、土壤侵蝕等研究中的基本因子. 坡長(zhǎng)能夠直接影響地面徑流的速度，進(jìn)而影響對(duì)地面土壤的侵蝕力. USLE模型中，坡長(zhǎng)L因子是指其他條件（降雨、土壤、耕作、水土保持措施、坡度等）相同的情況下，某一給定坡長(zhǎng)的坡面上，土壤流失量與標(biāo)準(zhǔn)徑流小區(qū)（坡長(zhǎng)22.13 m，坡度9%）典型坡面土壤流失量的比值介于0～1之間，反映地形地貌特征對(duì)土壤侵蝕的影響. Wischmeier和Smith［18］對(duì)坡長(zhǎng)指數(shù)的改進(jìn)結(jié)果在USLE模型中得到了廣泛應(yīng)用.

式中：L為坡長(zhǎng)因子；l為像元坡長(zhǎng)；m為坡長(zhǎng)指數(shù).m取值如下式：

傳統(tǒng)坡長(zhǎng)計(jì)算一般是在野外或地形圖上進(jìn)行，不僅費(fèi)時(shí)費(fèi)力，而且不能獲取研究區(qū)任意點(diǎn)的坡長(zhǎng). 地理信息系統(tǒng)、數(shù)字高程模型等技術(shù)的發(fā)展，使得坡長(zhǎng)計(jì)算的效率得到大幅提高，本文即是通過ArcGIS 軟件和DEM進(jìn)行坡長(zhǎng)提取. 由于GIS專業(yè)處理軟件中沒有直接求坡長(zhǎng)的功能，本文參考利用湯國(guó)安和楊昕［22］在水文分析功能和空間分析應(yīng)用時(shí)提出的一種快速近似計(jì)算坡長(zhǎng)的方法，在假定坡面水流方向與山脊線（分水嶺）垂直基礎(chǔ)上，計(jì)算每個(gè)點(diǎn)沿垂直方向到山脊線的實(shí)際匯流路徑長(zhǎng)度為坡長(zhǎng)的近似值.

計(jì)算步驟如下：

1）對(duì)DEM數(shù)據(jù)進(jìn)行填挖，再對(duì)填挖后的DEM數(shù)據(jù)提取負(fù)地形，獲得由原來溝道轉(zhuǎn)化成的山脊線.

2）求負(fù)地形的水流方向并計(jì)算流水累積量.

3）提取負(fù)地形的河網(wǎng). 使用條件函數(shù)工具con（accum>100，1），即像元流入量超過100 的所有像元均將成為河流網(wǎng)絡(luò)的一部分.

4）對(duì)山脊線進(jìn)行歐氏距離計(jì)算，得到近似坡長(zhǎng)l.

5）使用條件函數(shù)計(jì)算m值.

6）采用公式計(jì)算出坡長(zhǎng)因子L.

3.4 作物管理因子C

本研究區(qū)域的植被覆蓋度估算的數(shù)據(jù)源選用Landsat-8衛(wèi)星遙感影像，遙感植被指數(shù)選用歸一化差異植被指數(shù)（Normal Differential Vegetation Index，NDVI），方法選用像元二分模型. 植被是土壤侵蝕的主要抑制因素，C因子是這種抑制作用的定量化表達(dá)指標(biāo)之一.C因子指一定植被覆蓋條件耕作農(nóng)地上的土壤流失量與同等條件下實(shí)施清耕的連續(xù)休閑地的土壤流失量之比，其值大小變化于0～1 之間，C值越大，說明對(duì)應(yīng)的土地利用類型的土壤侵蝕越嚴(yán)重，其計(jì)算公式采用蔡崇法等［23］的研究結(jié)果.

3.5 治理措施因子P

USLE模型中，治理措施因子P定義為采用特定措施土地上的土壤流失量與順坡種植的土壤流失量的比值，介于0～1之間. 水土保持措施通過采取一定的耕作方式或改變小地形、攔蓄地表徑流、增加降雨入滲以充分利用光熱和水土資源，改善農(nóng)業(yè)生產(chǎn)條件達(dá)到保水保土保肥、減少土壤侵蝕和提高農(nóng)業(yè)生產(chǎn)的目的，主要包括耕作措施和工程措施兩大類. 據(jù)調(diào)查，流域內(nèi)的農(nóng)地主要為退耕還林后保留的基本農(nóng)田，大多已完成梯田改造；林地大部分是退耕還林工程造林，通常的侵蝕控制措施有梯田和淤地壩，未采取水土保持措施的土地利用類型P因子值為1. 具體參照劉寶元等的《土壤侵蝕預(yù)報(bào)模型》［24］并根據(jù)劃分的土地利用類型結(jié)合坡度信息對(duì)P因子進(jìn)行修正，確定修正后的P因子值，獲得P因子圖層.

4 結(jié)果與分析

4.1 土壤侵蝕強(qiáng)度等級(jí)變化

根據(jù)土壤侵蝕分類分級(jí)標(biāo)準(zhǔn)（SL 190—2007），豫西黃土丘陵地區(qū)的土壤侵蝕等級(jí)包括微度、輕度、中度、強(qiáng)度、極強(qiáng)度、劇烈等6 個(gè)級(jí)別［25-26］. 從表1 可以看出，2015 年到2020 年土壤侵蝕等級(jí)面積在豫西黃土丘陵地區(qū)發(fā)生變化，具體表現(xiàn)為：微度侵蝕中10 137.42 km2未變化，0.9 km2轉(zhuǎn)變?yōu)檩p度侵蝕；輕度侵蝕中6 566.85 km2未變化，100.08 km2轉(zhuǎn)變?yōu)槲⒍惹治g，1.34 km2轉(zhuǎn)變?yōu)橹卸惹治g，0.31 km2轉(zhuǎn)變?yōu)閺?qiáng)度侵蝕；中度侵蝕中3 580.29 km2未變化，7.22 km2轉(zhuǎn)變?yōu)槲⒍惹治g，69.59 km2轉(zhuǎn)變?yōu)檩p度侵蝕，0.9 km2轉(zhuǎn)變?yōu)閺?qiáng)度侵蝕，1.06 km2轉(zhuǎn)變?yōu)闃O強(qiáng)度侵蝕；強(qiáng)度侵蝕中2 469.92 km2未變化，0.69 km2轉(zhuǎn)變?yōu)槲⒍惹治g，35.33 km2轉(zhuǎn)變?yōu)檩p度侵蝕，17.39 km2轉(zhuǎn)變?yōu)橹卸惹治g，0.69 km2轉(zhuǎn)變?yōu)闃O強(qiáng)度侵蝕，0.46 km2轉(zhuǎn)變?yōu)閯×仪治g；極強(qiáng)度侵蝕中2 657.36 km2未變化，0.23 km2轉(zhuǎn)變?yōu)槲⒍惹治g，0.04 km2轉(zhuǎn)變?yōu)檩p度侵蝕，41.28 km2轉(zhuǎn)變?yōu)橹卸惹治g，4.59 km2轉(zhuǎn)變?yōu)閺?qiáng)度侵蝕，0.76 km2轉(zhuǎn)變?yōu)閯×仪治g；劇烈侵蝕中1 313.67 km2未變化，0.02 km2轉(zhuǎn)變?yōu)槲⒍惹治g，9.49 km2轉(zhuǎn)變?yōu)閺?qiáng)度侵蝕，2.51 km2轉(zhuǎn)變?yōu)闃O強(qiáng)度侵蝕.

表1 2015—2020年土壤侵蝕強(qiáng)度等級(jí)轉(zhuǎn)移矩陣Tab.1 Transfer matrix of soil erosion intensity grade from 2015 to 2020 單位：km2

4.2 土壤侵蝕時(shí)間變化

從表2可以看出，豫西黃土丘陵地區(qū)2015年和2020年從微度到劇烈的土壤侵蝕面積和占比的變化均呈整體減少趨勢(shì)，局部均呈先減少再增加再減少的雙峰波動(dòng)變化，峰值在強(qiáng)度和極強(qiáng)，土壤侵蝕面積和占比都在微度最大、劇烈最小. 從2015年到2020年，豫西黃土丘陵地區(qū)的土壤侵蝕面積和占比在不同侵蝕強(qiáng)度等級(jí)發(fā)生不同變化. 土壤侵蝕面積和占比分別為：微度分別增加105.62 km2、0.4%；輕度分別增加3.29 km2、0.02%；中度分別減少19.17 km2、0.07%；強(qiáng)度分別減少39.37 km2、0.14%；極強(qiáng)分別減少39.55 km2、0.15%；劇烈分別減少10.81 km2、0.04%. 由此可見，隨著時(shí)間推移，豫西黃土丘陵地區(qū)的強(qiáng)土壤侵蝕不斷減少，逐漸趨于輕微土壤侵蝕.

表2 2015年和2020年豫西黃土丘陵地區(qū)各土壤侵蝕強(qiáng)度的面積Tab.2 Area of soil erosion intensity in loess hilly area of western Henan from 2015 to 2020

4.3 土壤侵蝕空間變化

由圖2 可以看出，豫西黃土丘陵地區(qū)的大部分地區(qū)發(fā)生輕、微、中程度的土壤侵蝕，少部分地區(qū)發(fā)生強(qiáng)度以上土壤侵蝕. 多數(shù)輕微程度土壤侵蝕出現(xiàn)在研究區(qū)的東北部、中部以及零散分布在其他地區(qū)，而強(qiáng)度以上土壤侵蝕出現(xiàn)在研究區(qū)的北部和南部邊緣地帶以及東南部的部分地區(qū)，從行政區(qū)上看，中度侵蝕在盧氏縣、靈寶市最廣泛，強(qiáng)度侵蝕在盧氏縣、靈寶市、嵩縣最廣泛，極強(qiáng)度侵蝕在盧氏縣、嵩縣最廣泛，劇烈侵蝕在嵩縣最廣泛.

圖2 2020年豫西黃土丘陵地區(qū)的土壤侵蝕強(qiáng)度的空間分布圖Fig.2 Spatial distribution of soil erosion intensity in loess hilly area of western Henan in 2020

5 結(jié)論

本文利用RUSLE 模型對(duì)豫西黃土丘陵地區(qū)2015—2020年間區(qū)域水土保持功能的變化特征進(jìn)行研究，得到如下結(jié)論：

1）與2015年相比，區(qū)域的水土保持功能得到了顯著了提升，區(qū)域平均侵蝕模數(shù)總體屬于中度侵蝕量級(jí)，對(duì)全區(qū)域平均土壤侵蝕量下降率約為1.03%. 從不同侵蝕量級(jí)看，中度及以上面積都有所下降，總面積下降了約112 km2；微度面積增加明顯，增加了約106 km2.

2）2020年，豫西區(qū)域土壤侵蝕量級(jí)總體上以中度及以下為主，占比為76.09%，其中微度為37.93%，輕度24.69%，中度13.47%；強(qiáng)度及以上占比為23.91%，其中強(qiáng)度占9.2%，極強(qiáng)占9.85%，劇烈占4.87%. 土壤侵蝕防治重點(diǎn)區(qū)域主要分布于區(qū)域南部洛寧、汝陽(yáng)、嵩縣，西部的靈寶，北部的濟(jì)源、澠池等.