基于MODIS和Landsat數(shù)據(jù)的湟水河流域土壤侵蝕時空變化研究

杜 梅, 趙健赟, 楊 靜, 丁圓圓, 劉文惠, 李國榮, 王祖順, 趙利江

(1.青海大學 地質(zhì)工程系, 青海 西寧 810016; 2.青海省剛察縣氣象站, 青海 剛察 812300;3.甘肅省祁連山生態(tài)環(huán)境研究中心, 甘肅 蘭州 730000; 4.青海省基礎(chǔ)測繪院, 青海 西寧 810001)

土壤侵蝕是土體在特定時間和空間條件下的遷移過程，是最重要的土地退化問題，它直接影響水、土資源的利用與保護，危害區(qū)域生態(tài)環(huán)境，妨礙區(qū)域社會經(jīng)濟和可持續(xù)發(fā)展，是區(qū)域水土保持、荒漠化防治以及生態(tài)保護領(lǐng)域的關(guān)鍵科學問題之一[1]。在全球范圍內(nèi)，因水力剝蝕、搬運作用形成的水蝕面積達1.09×108km2，因氣流沖擊、風力剝蝕形成的風蝕侵蝕面積近5.78×107km2，不考慮凍融侵蝕以及其他形式的侵蝕，全球土壤侵蝕面積已達1.64×108km2，而中國因各種因素造成不同形式的土壤侵蝕遍布全國、態(tài)勢嚴重[2]。湟水河流域作為黃河源區(qū)重要的生態(tài)屏障，地處高寒生態(tài)脆弱區(qū)[3]，隨著全球氣候及環(huán)境變化，流域內(nèi)土壤出現(xiàn)了不同程度的退化及侵蝕現(xiàn)象。因此，開展湟水河流域尺度的土壤侵蝕時空變化研究，對黃河上游地區(qū)水土流失治理具有重要的意義。

近幾十年來，國內(nèi)外學者在土壤侵蝕方面發(fā)展了通用土壤流失方程USLE(universal soil loss equation)、分布式土壤水力侵蝕預(yù)報模型WEPP(water erosion prediction project)等區(qū)域土壤侵蝕定量分析方法，在不同地區(qū)開展了土壤侵蝕調(diào)查、評估、預(yù)測等方面的研究[4-6]；同時，在土壤侵蝕動力機制與過程模擬，土壤侵蝕與物質(zhì)遷移，土壤侵蝕與氣候變化等方面取得了進展[7-9]，發(fā)現(xiàn)氣候變化可通過改變降水徑流、地表植被覆蓋度以及人類活動影響侵蝕過程，未來土壤侵蝕的變化可通過氣候變化模式與土壤侵蝕模型的耦合來進行預(yù)測[10-11]，而在激光雷達和高時空分辨率遙感等技術(shù)的支持下，能夠?qū)⒌乇硇畔⒕毣磉_，進一步增強土壤侵蝕的模擬[12]。另外，為適應(yīng)不同區(qū)域、不同尺度的土壤侵蝕預(yù)報，構(gòu)建了中國土壤流失方程CSLE(Chinese soil loss equation)[13]、AGNPS(agricultural non-point source)模型[14]和坡面尺度的水蝕預(yù)報模型[15]，等。修正的通用土壤流失方程RUSLE(revised universal soil loss equation)在不同地區(qū)土壤侵蝕研究中也有應(yīng)用，開展了諸如土壤侵蝕量估算、土壤侵蝕強度分析、土壤侵蝕時空演變及其空間分異特征等方面的研究[16-18]。但是，在土壤侵蝕模型模擬結(jié)果的檢驗、土壤侵蝕動力機制等方面，特別是應(yīng)用先進的對地觀測技術(shù)開展實地調(diào)查與驗證工作有待深入研究。本研究基于2000年和2018年的MODIS， Landsat，降雨、人口密度、經(jīng)濟等多源數(shù)據(jù)，利用修正的通用土壤流失方程RUSLE，分析湟水河流域土壤侵蝕的特點與時空變化規(guī)律，并基于無人機遙感技術(shù)對結(jié)果進行實地驗證與評價，以進一步明確湟水河流域近20 a來的土壤侵蝕時空變化特征及驅(qū)動因素，以期為黃河上游地區(qū)的水土保持與防治工作提供數(shù)據(jù)基礎(chǔ)與決策依據(jù)。

1 研究區(qū)概況

湟水河是黃河流域重要的水源地之一，也是黃河上游最大的一級支流，處于青藏高原與黃土高原的交錯地帶，位于北緯 36°05′—38°30′，東經(jīng)99°00′—103°43′，海拔1 545～5 232 m，年均氣溫5.2～6.8 ℃，年降水量350～550 mm，屬高海拔半干旱地區(qū)，地貌過渡性和復(fù)雜性明顯，流域水系呈樹葉狀分布，干流南北兩岸支流發(fā)育[19]。湟水河源區(qū)多沼澤濕地，流經(jīng)地區(qū)多峽谷、盆地相間分布，存在高寒草甸地貌，分布優(yōu)良牧場和良好草地；中游一帶平坦開闊，存在有較大面積的灌溉農(nóng)田和非灌溉農(nóng)田分布；中下游主要為石林山區(qū)、土石山區(qū)、黃土丘陵溝壑區(qū)；流域東南部多分布寒鈣土、新積土和冷鈣土。流域西部的山地以灰褐土和黑鈣土為主，流域中部的中高山地以黑鈣土、黃土和壤土為主。研究區(qū)總?cè)丝诩s為3.30×107人，農(nóng)田耕地超過3.00×106hm2，是青海省主要的人口居住區(qū)和產(chǎn)糧區(qū)[20-21]。

2 數(shù)據(jù)來源與研究方法

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

本研究主要使用MODIS，Landsat，數(shù)字高程數(shù)據(jù)DEM、降雨、人口密度和經(jīng)濟等數(shù)據(jù)(產(chǎn)品)。

MODIS數(shù)據(jù)來源于美國航空航天局數(shù)據(jù)中心(https:∥www.earthdata.nasa.gov/)，選用2000，2018年行列號為h26 v05和h25 v05的MODIS MOD13Q1產(chǎn)品數(shù)據(jù)，空間分辨率為250 m，采用最大值合成法獲得兩個年度的NDVI數(shù)據(jù)；

Landsat數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http:∥www.gscloud.cn/)，選取2000，2018年4月至10月軌道號為131—134，行編號為33—35，且云量小于5%的7幅影像，空間分辨率為30 m，對其進行輻射定標、幾何糾正、鑲嵌、分類等處理，獲得兩個年度的土地覆蓋數(shù)據(jù)；DEM數(shù)據(jù)來源地理空間數(shù)據(jù)云(http:∥www.gscloud.cn/)，空間分辨率為30 m，對該數(shù)據(jù)進行流域分析，獲得研究區(qū)流域邊界與河網(wǎng)信息； 2000，2018年研究區(qū)的降雨數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)(http:∥www.data.cma.cn)，人口密度和經(jīng)濟數(shù)據(jù)來源于中國資源環(huán)境科學與數(shù)據(jù)中心(http:∥www.resdc.cn)，對上述數(shù)據(jù)進行投影轉(zhuǎn)換，并利用流域邊界數(shù)據(jù)對其進行裁切，獲得研究區(qū)土壤侵蝕分析的數(shù)據(jù)集；

無人機遙感數(shù)據(jù)是利用大疆精靈4 RTK于2019年8月分別在研究區(qū)現(xiàn)場采集，設(shè)置航向重疊為80%，旁向重疊為70%，航高為120 m。將上述數(shù)據(jù)的空間分辨率統(tǒng)一重采樣至250 m，均使用 Albers Conical Equal Area投影。

2.2 研究方法

2.2.1 土壤侵蝕模型 本研究選用修正的通用土壤流失方程RUSLE計算研究區(qū)土壤侵蝕模數(shù)，公式如下[22]：

A=R×K×LS×C×P

(1)

式中：A為年均土壤侵蝕模數(shù)〔t/(hm2·a)〕;R為降雨侵蝕力因子〔(MJ·mm)/(hm2·h·a)〕;K為土壤可蝕性因子〔(t·h)/(MJ·mm)〕;LS為地形(坡長坡度)因子,無量綱;C為植被覆蓋和作物管理因子,無量綱;P為水土保持措施因子,無量綱。

2.2.2 降雨侵蝕力因子R值的確定 根據(jù)己有湟水河流域地區(qū)月度氣象數(shù)據(jù)資料和模型的適用性，本研究采用Wischmeier等[23]提出的R值計算方法，公式如下：

(2)

式中：R為降雨侵蝕力值;Pi為月均降雨量(mm);P為年均降雨量。

根據(jù)青海省2000年和2018年22個氣象站點的降雨數(shù)據(jù)，計算得到各個氣象站的R值，經(jīng)克里金空間插值計算得到研究區(qū)降雨侵蝕力因子R值分布，見圖1。

圖1 湟水河流域2000，2018年降水侵蝕力因子R值空間分布

2.2.3 土壤可蝕性因子K值的確定 土壤可蝕性因子K與土壤粉粒、黏粒等含量有關(guān)，其計算方法如下[24]：

(3)

其中：

(4)

式中：SAN為砂粒含量(%); SIL為粉粒含量(%); CLA為黏粒含量(%);C為土壤有機碳含量(%);Ka為土壤可蝕性因子。利用上述公式計算獲得研究區(qū)土壤可蝕性因子的空間分布(圖2)。

圖2 湟水河流域土壤可蝕性因子K值的空間分布

2.2.4 坡度坡長因子LS值的確定 基于DEM數(shù)據(jù)計算研究區(qū)坡度、坡長數(shù)據(jù)，進而通過分段函數(shù)計算坡度因子S值和坡長因子L值，進而計算坡度坡長因子LS，公式如下[2]：

(5)

(6)

(7)

LS=L×S

(8)

式中：λ為坡長(m);θ為傾斜角(°);S為百分比坡度(%)。利用上述公式計算獲得研究區(qū)坡度坡長因子LS值的空間分布(圖3)。

2.2.5 植被覆蓋與作物管理因子C值的確定 以獲得的研究區(qū)NDVI數(shù)據(jù)為基礎(chǔ)，經(jīng)處理得到植被覆蓋度后[22]，由蔡崇法等[25]提出的植被覆蓋度Fc與作物管理因子的轉(zhuǎn)換公式計算C值：

(9)

式中：C為植被覆蓋與作物管理因子;Fc為植被覆蓋度,其計算方法如下:

(10)

式中：Fc為年作物生長期的植被覆蓋度;N為像元NDVI值;Nsoil為裸地像元NDVI值;Nveg為純植被覆蓋像元NDVI值;Nveg和Nsoil分別為NDVI統(tǒng)計結(jié)果在累積概率為95%和5%時的取值。利用上述公式計算獲得研究區(qū)2000，2018年植被覆蓋與作物管理因子C值的空間分布(圖4)。

圖3 湟水河流域坡度坡長因子LS值的空間分布

圖4 湟水河流域2000，2018年植被覆蓋與作物管理因子C值空間分布

2.2.6 水土保持措施因子p值的確定 采用Lufafa等[6]提出的方法計算水土保持措施因子p值，公式如下：

P=0.2+0.03×θ

(11)

式中：θ為坡度,由此計算獲得研究區(qū)水土保持措施因子p值的空間分布(圖5)。

圖5 湟水河流域水土保持措施因子P值空間分布

2.2.7 土壤侵蝕強度的驗證 在土壤侵蝕典型區(qū)域選擇實地調(diào)查樣地進行模型模擬結(jié)果的驗證，無人機數(shù)據(jù)空間分辨率優(yōu)于0.05 m的影像數(shù)據(jù)，以250 m的格網(wǎng)為單元通過目視解譯法解譯調(diào)查樣地的土地利用[26]、水土保持措施，并計算植被覆蓋度和坡度，依據(jù)《水土保持遙感監(jiān)測技術(shù)規(guī)范(SL592-2012)》《土壤侵蝕分類分級標準(SL190-2007)》確定土壤侵蝕強度，計算總體精度[2]：

(12)

3 結(jié)果與分析

3.1 2000年土壤侵蝕模型計算分析

基于2000年土壤侵蝕模型因子計算結(jié)果，利用式(1)計算獲得研究區(qū)2000年土壤侵蝕模數(shù)，其最大值為15 706.50 t/(km2·a)，均值為477.81 t/(km2·a)。依據(jù)《土壤侵蝕分類分級標準》和祁連山地區(qū)土壤侵蝕的相關(guān)研究結(jié)果[27]，利用表1對計算獲得的土壤侵蝕模數(shù)結(jié)果進行等級劃分，結(jié)果見圖6。

表1 土壤侵蝕強度等級劃分

由圖6可知，2000年研究區(qū)整體土壤侵蝕強度不大，微度侵蝕主要分布在湟水河源區(qū)及山間河谷地區(qū)的耕地和高植被覆蓋度的草地、有林地；輕度、中度侵蝕主要分布在海拔較高、植被覆蓋少的山地、荒地。

圖6 2000年湟水河流域土壤侵蝕強度分布

3.2 2018年土壤侵蝕模型計算分析

基于2018年土壤侵蝕模型因子計算結(jié)果，利用式(1)計算獲得研究區(qū)2018年土壤侵蝕模數(shù)，其最大值為32 653.50 t/(km2·a)，均值為1 625.30 t/(km2·a)，利用表1對計算獲得的土壤侵蝕模數(shù)結(jié)果進行等級劃分，結(jié)果見圖7。

由圖7可知，2018年研究區(qū)主要以微度侵蝕為主，分布在湟水河源區(qū)、高山之間的河谷盆地地區(qū)，輕度侵蝕主要分布在山區(qū)及其周圍的耕地、坡耕地等區(qū)域，中度侵蝕主要分布在研究區(qū)東南部城鎮(zhèn)居民聚集地帶與河流灘地，強烈侵蝕和極強烈侵蝕零散分布于祁連山、達坂山等高山、禿嶺裸地地區(qū)。

圖7 2018年湟水河流域土壤侵蝕強度分布

3.3 土壤侵蝕模型的驗證分析

由以上獲得的研究區(qū)土壤侵蝕各個模型因子計算結(jié)果可知，降水侵蝕力因子R值呈東北低西南高的特點，土壤可蝕性因子K值較大區(qū)域分布在湟水河源區(qū)和中部河流兩側(cè)地區(qū)，而坡度坡長因子LS值在中部及中部偏北地區(qū)較高，植被覆蓋與作物管理因子C值東南高西北低，水土保持措施因子P值在湟水河源區(qū)普遍較低，在北部祁連山脈、中部達坂山和南部拉脊山一帶較高。因此，分別在研究區(qū)北部的門源、中部的大通和南部的湟中選擇總面積約3 km2的樣地，采集的無人機高分辨率影像數(shù)據(jù)(見圖8)，建立解譯標志。

圖8 湟水河流域土壤侵蝕強度無人機高分辨率驗證數(shù)據(jù)

依據(jù)解譯標志和土壤侵蝕強度驗證方法獲得49個格網(wǎng)的侵蝕強度，并依據(jù)公式(12)計算土壤侵蝕驗證總體精度P為91.8%，表明模型計算的研究區(qū)土壤侵蝕強度可靠。

3.4 土壤侵蝕強度時空變化分析

統(tǒng)計研究區(qū)2000年和2018年各個侵蝕強度等級的面積，繪制變化對比圖(圖9)，發(fā)現(xiàn)微度侵蝕面積比例顯著減少，中度、強烈、極強烈侵蝕面積顯著增加，而輕度、劇烈侵蝕的面積雖有增加，但幅度不大，表明研究區(qū)土壤侵蝕面積有向中度、強烈和極強烈變化的趨勢。對2000年和2018年研究區(qū)土壤侵蝕強度進行逐像元求差，得到0，1，2，4的運算結(jié)果，分別代表未發(fā)生變化、微度、中度和劇烈變化像元，空間分布如圖10所示。

圖9 2000—2018年湟水河流域各土壤侵蝕強度等級面積變化

圖10 2000—2018年湟水河流域土壤侵蝕強度空間變化

由圖10可知，在湟水河源區(qū)、祁連山與達坂山之間的支流大通河谷地等區(qū)域土壤侵蝕強度未發(fā)生明顯變化，在沿湟水河干流中下游的城鎮(zhèn)居民聚集地帶存在土壤侵蝕強度微度變化，面積約占28.12%；在祁連山、達坂山的高山裸地、禿嶺地區(qū)存在土壤侵蝕強度中度變化，面積約占8.91%，而土壤侵蝕強度劇烈變化區(qū)域則零散分布在高山裸地及城鎮(zhèn)地帶。

4 討論與結(jié)論

4.1 討 論

針對湟水河流域及其周邊地區(qū)的土壤侵蝕過程與機制問題，已有學者開展了一定程度的研究。如王雪璐[28]利用RUSLE模型估算青藏高原三江源區(qū)多年平均侵蝕量為3.10×109t/a；陳豪等[24]利用USLE模型對祁連山公園的土壤侵蝕進行了模擬計算，得到2015年和2019年的平均侵蝕模數(shù)為1 063 t/(km2·a)和1 068 t/(km2·a)；徐進等[29]對湟水河(西寧段)流域的土壤侵蝕危險性進行了綜合評價，發(fā)現(xiàn)土壤侵蝕危險性較高的地方主要集中在城鎮(zhèn)居民地聚集的市區(qū)及湟水河河谷地帶；魏素娟等[1]基于USLE模型對湟水流域土壤侵蝕進行了評估，得到湟水流域土壤侵蝕模數(shù)為523.8 t/(km2·a)，年侵蝕總量為8.37×107t/a。這些結(jié)果與本研究得到的湟水河流域土壤侵蝕模數(shù)和土壤侵蝕量基本一致，結(jié)果可靠。同時，本研究還利用低空無人機遙感技術(shù)對模型計算的土壤侵蝕強度結(jié)果進行了定量評價，具有一定的創(chuàng)新。

2000—2018年，研究區(qū)土壤侵蝕量增加，發(fā)生中、重度變化的區(qū)域主要集中在祁連山、達坂山的高山裸地、禿嶺地區(qū)，沿湟水河干流的城鎮(zhèn)居民地聚集地帶的侵蝕強度也存在加劇的趨勢。從本研究采用的RUSLE模型來看，土壤侵蝕過程是氣候、植被覆蓋等相對動態(tài)因子和地形、土壤等相對靜態(tài)因子共同驅(qū)動的結(jié)果。相對靜態(tài)因子主要控制土壤侵蝕的空間分布，不同的土壤類型孕育不同的植被，形成地域性特征；而相對動態(tài)因子控制土壤侵蝕的時空變化，即植被、氣候等因子對土壤侵蝕起降低或加劇作用。根據(jù)研究區(qū)氣象、人口、經(jīng)濟等數(shù)據(jù)，發(fā)現(xiàn)19 a來研究區(qū)平均降水量增加了113.12 mm，建設(shè)用地增加了5.25 km2，人口密度增加了30人/km2，而人口密度和建設(shè)用地的增加在一定程度上影響了地表植被覆蓋度，再加上氣候因子的持續(xù)性改變，導致相對動態(tài)因子逐步驅(qū)動研究區(qū)土壤侵蝕強度加劇。由此表明，研究區(qū)土壤侵蝕強度的增加，與氣候暖濕化、人類活動強度增加存在一定的關(guān)系。

2000—2018年，湟水河源區(qū)、大通河谷地等區(qū)域土壤侵蝕強度未發(fā)生明顯變化，在中下游城鎮(zhèn)居民聚集地帶發(fā)生微度變化，而高山裸地、禿嶺地區(qū)土壤侵蝕強度發(fā)生中度、強烈的變化，存在顯著的空間差異性。未發(fā)生明顯變化的區(qū)域人口密度低、人類活動以放牧為主，植被覆蓋穩(wěn)定，植被控制由于降雨量增加侵蝕強度變大的作用顯著；位于研究區(qū)中下游的城鎮(zhèn)居民區(qū)聚集有青海省60%的人口，但在生態(tài)工程建設(shè)的有效驅(qū)動下，增加了地表的植被覆蓋度和水土保持功能，在一定程度上阻滯了由于人類活動加劇和氣候變化帶來的侵蝕強度增大；而裸地、禿嶺地區(qū)由于缺乏有針對性的水保措施，土壤與植被相互影響和作用，出現(xiàn)侵蝕加劇與覆蓋度逐步降低的不良循環(huán)。

由于數(shù)據(jù)源時空分辨率的差異，導致模型計算結(jié)果的空間分辨率較低，且氣象站數(shù)據(jù)較少造成降雨侵蝕力因子的精度也不高，下一步將結(jié)合土壤侵蝕定位觀測試驗進行更加精細的研究，并開展湟水河流域土壤侵蝕驅(qū)動機制等方面的研究。

4.2 結(jié) 論

(1) 2000年研究區(qū)整體土壤侵蝕強度不大，以微度侵蝕為主，在祁連山、達坂山及拉脊山海拔較高、植被覆蓋較少的地區(qū)存在中度侵蝕；2018年研究區(qū)中度侵蝕主要分布在研究區(qū)東南部城鎮(zhèn)居民聚集地帶與河流灘地，強烈侵蝕和極強烈侵蝕主要分布在祁連山、達坂山等高山、禿嶺裸地地區(qū)。

(2) 2000—2018年，平均土壤侵蝕模數(shù)存在增加的趨勢，輕度侵蝕量占比顯著減少，強烈、極強烈和劇烈侵蝕量占比顯著增加，這與研究區(qū)氣候暖濕化、人類活動強度增加存在一定的關(guān)系。19 a間，湟水河源區(qū)、大通河谷地等區(qū)域土壤侵蝕強度未發(fā)生明顯變化，祁連山、達坂山的高山裸地、禿嶺地區(qū)的土壤侵蝕發(fā)生中、重度變化，沿湟水河干流的城鎮(zhèn)居民地聚集地帶的侵蝕強度也存在加劇的趨勢，存在較為明顯的空間差異性。這種差異性受區(qū)域人類活動特點、植被覆蓋情況、生態(tài)工程建設(shè)、自然地理特性等因素的影響。