馬 蘇,劉軍會,康玉麟,劉 洋,鄭穎娟,周甲男
中國環(huán)境科學(xué)研究院環(huán)境信息研究所,北京 100012
土壤風(fēng)蝕是干旱與半干旱地區(qū)廣泛存在的嚴(yán)重環(huán)境問題之一[1]. 土壤風(fēng)蝕導(dǎo)致的土地沙化和土地退化降低了土壤的生產(chǎn)力,削弱了土壤保持服務(wù)功能[2-3],威脅著區(qū)域經(jīng)濟(jì)發(fā)展、糧食安全和人類福祉[4-5]. 根據(jù)2015年發(fā)布的《中國荒漠化和沙化狀況公報》,截至2014年,我國仍有荒漠化土地361.16×104km2和沙化土地172.12×104km2. 為遏制土地沙化與環(huán)境惡化的趨勢,國家相繼實施了一系列重大的生態(tài)保護(hù)與建設(shè)工程,如“三北”防護(hù)林、黃土高原中上游退耕還林還草、京津風(fēng)沙源治理等,使得我國北方地區(qū)土壤風(fēng)蝕得到了有效緩解[6-7],但是生態(tài)系統(tǒng)的恢復(fù)是漫長復(fù)雜的過程,及時監(jiān)測、評估工程區(qū)生態(tài)狀況及功能變化尤為重要[8]. 而現(xiàn)有研究多側(cè)重于植被覆蓋狀況的監(jiān)測,對區(qū)域防風(fēng)固沙功能的動態(tài)變化及其空間差異的研究相對較少.
防風(fēng)固沙是生態(tài)系統(tǒng)通過其結(jié)構(gòu)與過程減少由于風(fēng)蝕所導(dǎo)致的土壤侵蝕的作用[9],這種保持土壤、抑制風(fēng)蝕過程的功能即為防風(fēng)固沙功能[10-11]. 風(fēng)蝕模型是評估防風(fēng)固沙功能的主要技術(shù)手段,并已經(jīng)發(fā)展演化出通用風(fēng)蝕方程(WEQ)[12]、德克薩斯侵蝕分析模型(TEAM)[13]、Bocharov模型[14]、風(fēng)蝕預(yù)報系統(tǒng)(WEPS)[15]和修正的風(fēng)蝕模型(RWEQ)[16-17]等多種模型. 其中,RWEQ是目前應(yīng)用最廣泛的土壤風(fēng)蝕模型[18-23],該模型充分考慮氣候變化、植被覆蓋狀況、土壤質(zhì)地、地表粗糙度等因素,通過參數(shù)的修正和公式調(diào)整可以應(yīng)用到我國北方土壤風(fēng)蝕狀況的評估[24-32].
鄂爾多斯市地處我國北方防沙屏障帶的南側(cè),境內(nèi)分布有庫布齊沙漠和毛烏素沙地,土壤風(fēng)蝕嚴(yán)重,生態(tài)環(huán)境脆弱,是京津地區(qū)主要風(fēng)沙源之一,其防風(fēng)固沙能力不僅影響當(dāng)?shù)厣鷳B(tài)屏障建設(shè),而且關(guān)系到環(huán)京津地區(qū)的環(huán)境安全. 因此該研究基于長時間序列衛(wèi)星遙感影像以及氣象、植被、土壤等數(shù)據(jù)資料,采用修正的風(fēng)蝕模型(RWEQ)與GIS技術(shù),研究2000?2018年鄂爾多斯市防風(fēng)固沙功能的時空變化規(guī)律,并從氣候變化和人類活動角度探討驅(qū)動影響機理,以期為鄂爾多斯市科學(xué)實施生態(tài)修復(fù)和筑牢生態(tài)安全屏障提供科學(xué)參考.
鄂 爾 多 斯 市 位于106°42′40″E~111°27′20″E、37°35′24″N~40°51′40″N之 間,總 面 積 約8.68×104km2(見圖1). 地處黃河“幾字彎”腹地,三面黃河環(huán)抱,平均海拔1 000~1 500 m,北部為黃河沖積平原,東部為丘陵溝壑水土流失區(qū)和砒砂巖裸露區(qū),西北部為庫布齊沙漠,東南部為毛烏素沙地,沙化土地約占全市總面積的48%. 氣候類型為溫帶大陸性干旱半干旱氣候,年降水量170~350 mm,年蒸發(fā)量高達(dá)2 000~3 000 mm;多大風(fēng),年均風(fēng)速為3.6 m/s,最大風(fēng)速可達(dá)22 m/s. 植被以典型草原、荒漠草原和草原化荒漠為主. 土壤類型主要包括栗鈣土、棕鈣土和風(fēng)沙土等.
圖 1 鄂爾多斯市的地理位置示意Fig.1 The location of Ordos
遙感數(shù)據(jù):歸一化植被指數(shù)(NDVI)為MOD13Q1產(chǎn)品,該數(shù)據(jù)已經(jīng)經(jīng)過幾何精糾正、輻射校正、大氣校正等預(yù)處理. 利用MRT工具進(jìn)行投影和格式轉(zhuǎn)換批處理,然后采用最大值合成法(maximum value composite)獲得NDVI數(shù)據(jù),并基于像元二分法生產(chǎn)2000?2018年的植被覆蓋度.
土地覆被數(shù)據(jù):使用2000年的Landsat-TM遙感影像數(shù)據(jù)和2018年的Landsat8遙感影像數(shù)據(jù),結(jié)合高分遙感影像,采用CART決策樹分類分別得到2000年和2018年的土地覆被圖.
氣象數(shù)據(jù):來源于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http://cdc.cma.gov.cn)以及鄂爾多斯市氣象局8個氣象站點的數(shù)據(jù). 數(shù)據(jù)內(nèi)容包括鄂爾多斯市各氣象站點的編號、經(jīng)緯度和海拔,以及每個氣象站點在相應(yīng)分析時間尺度內(nèi)的降水量、風(fēng)速、氣溫、日照時數(shù)等;采用Kriging對長時間序列離散點的氣象數(shù)據(jù)進(jìn)行插值處理,以得到氣象要素的柵格圖像.
土壤數(shù)據(jù):來源于中國科學(xué)院南京土壤研究所提供的1:1 000 000土壤圖及所附的土壤屬性表,包括有機質(zhì)、砂粒和黏粒含量等屬性.
其他數(shù)據(jù):雪蓋數(shù)據(jù)采用中國科學(xué)院旱區(qū)寒區(qū)科學(xué)數(shù)據(jù)中心(http://westdc.westgis.ac.cn)的中國雪深長時間序列數(shù)據(jù)集計算得到. 數(shù)字高程(DEM)來自中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心(http://www.resdc.cn),空間分辨率為90 m.
所有數(shù)據(jù)在ArcGIS中進(jìn)行運算時投影坐標(biāo)系統(tǒng)均為Krasovsky_1940_Albers,數(shù)據(jù)均重采樣為250 m的空間精度.
RWEQ模型是一種以較高時空分辨率對區(qū)域土壤風(fēng)蝕狀況進(jìn)行長時間序列估算,從而有效預(yù)測防風(fēng)固沙量的模型[33-36]. 因此,該研究采用RWEQ模型測算防風(fēng)固沙量,進(jìn)而評估鄂爾多斯市防風(fēng)固沙功能及其動態(tài)變化,用單位面積防風(fēng)固沙量的多少代表防風(fēng)固沙服務(wù)功能的強弱. 防風(fēng)固沙量(SRA)為潛在土壤風(fēng)蝕量( SLp) 與實際土壤風(fēng)蝕量(SLr)的差值,具體計算公式如下:
式中:SLp為潛在風(fēng)蝕量,kg/m2;Z為下風(fēng)向距離;sr為潛在關(guān)鍵地塊長度,m;Qpmax為潛在風(fēng)力的最大輸沙能力,kg/m;sp為實際關(guān)鍵地塊長度,m;WF為氣象因子,kg/m;EF和SCF分別為土壤可蝕性因子和土壤結(jié)皮因子;K和C分別為土壤粗糙因子與植被因子;SLr為實際風(fēng)蝕量,kg/m2;Qrmax為實際風(fēng)力的最大輸沙能力,kg/m;SRA表示單位面積防風(fēng)固沙量,kg/m2.
a) 氣象因子(WF). 氣象因子表征了在考慮降雨、溫度、日照及雪蓋等因素下,風(fēng)力對土壤顆粒的搬運能力.
式中:g為重力加速度,m/s2,取值9.8 m/s2;ρ為空氣密度,kg/m3,該文取1.205 kg/m3;SW為土壤濕度因子;SD為雪蓋因子(無積雪覆蓋天數(shù)與研究總天數(shù)之比,雪蓋深度<25.4 mm為無積雪覆蓋),結(jié)合鄂爾多斯市的地理位置以及雪蓋數(shù)據(jù)的處理結(jié)果,該文取1;Ui為2 m處的風(fēng)速,m/s;Ut為2 m處的臨界風(fēng)速,筆者所在課題組于2021年8?9月在伊金霍洛旗選取22個土壤剖面,采用TDR100時域反射儀測度土壤水分含量,發(fā)現(xiàn)土壤含水量均在1%左右,而劉玉章等[37]的風(fēng)洞試驗研究表明,起沙風(fēng)速一般隨沙中含水量的增加而增加,在含水量為0.8%左右時,起沙風(fēng)速為4.5~5.0 m/s,所以鄂爾多斯市的臨界風(fēng)速取5 m/s;Nd為觀測風(fēng)速大于臨界風(fēng)速的天數(shù),d;N為一次試驗中觀察的總次數(shù),該研究中進(jìn)行了1年累計,所以設(shè)定N為365.
b) 土壤可蝕性因子(EF)與土壤結(jié)皮因子(SCF).土壤可蝕性因子與土壤結(jié)皮因子的計算公式分別如式(9)(10)所示.
式中:sa為土壤粗砂含量(0.2~2 mm),%;si為土壤粉砂含量,%;cl為土壤黏粒含量,%;[OM]為土壤有機質(zhì)含量,%;[CaCO3]為碳酸鈣含量,%.
c) 植被覆蓋因子(C). 植被覆蓋度是影響風(fēng)蝕的關(guān)鍵因子,其覆蓋程度直接影響近地表風(fēng)速及土壤粗糙度.
式中:SC為植被覆蓋度,%;NDVImin為裸地對應(yīng)NDVI值,取累計頻率為2%的值;NDVImax為NDVI最大值,取累計頻率為98%的值.
d) 地表粗糙因子(K). 地表粗糙因子的計算公式如式(13)所示.
式中,slope為地形坡度,根據(jù)DEM利用ArcGIS軟件中的Slope工具實現(xiàn).
該方法在柵格尺度上對防風(fēng)固沙量的變化趨勢進(jìn)行模擬,能夠反映研究時段內(nèi)防風(fēng)固沙量變化幅度的空間分布,其計算公式:
式中:SLOPE為回歸方程的斜率;m為數(shù)據(jù)集時間段的年數(shù);Fj為某像元第j年的防風(fēng)固沙功能,t/km2.若SLOPE為正,表示隨著時間變化,防風(fēng)固沙功能呈上升趨勢;反之,當(dāng)SLOPE值為負(fù)數(shù)時,表示隨著時間變化,防風(fēng)固沙功能呈現(xiàn)下降趨勢;SLOPE的數(shù)值大小反映防風(fēng)固沙功能變化的程度.
Pearson相關(guān)系數(shù)能夠表征2個變量之間的相關(guān)程度,該研究以像元為計算單元,分析降水對防風(fēng)固沙功能的影響,其計算公式為
式中:x為像元的防風(fēng)固沙功能,t/km2;y為降水量,mm.
為了更好地表達(dá)防風(fēng)固沙功能的空間差異,依據(jù)研究區(qū)防風(fēng)固沙功能值域及自然斷裂法,設(shè)置防風(fēng)固沙功能在0~0.7×104t/km2之間為低值區(qū),在0.7×104~2.9×104t/km2之間為較低區(qū),在2.9×104~5.0×104t/km2之間為一般區(qū),在5.0×104~7.8×104t/km2之間為較高區(qū),在7.8×104t/km2以上為高值區(qū).
評估結(jié)果表明,2018年鄂爾多斯市防風(fēng)固沙總量為95.07×108t,單位面積防風(fēng)固沙量為10.95×104t/km2. 空間上,防風(fēng)固沙功能存在明顯地域差異(見圖2),高值區(qū)域主要分布在鄂爾多斯市東部和北部黃河沿岸,包括準(zhǔn)格爾旗、達(dá)拉特旗東部和北部、東勝區(qū)、康巴什區(qū)、伊金霍洛旗東部等地區(qū),面積占比為25.14%;較高區(qū)域零散分布在鄂爾多斯市中部和南部,包括鄂托克旗西北部、鄂托克前旗大部、烏審旗南部、達(dá)拉特旗南部等區(qū)域,面積占比為33.03%;一般區(qū)域集中分布在杭錦旗、鄂托克旗和烏審旗的部分區(qū)域,面積占比為15.58%;較低區(qū)和低值區(qū)主要分布在庫布齊沙漠和毛烏素沙地,面積占比分別為15.27%、10.97%.
圖 2 2018年鄂爾多斯市防風(fēng)固沙功能的空間分布Fig.2 Spatial pattern of wind-breaking and sand-fixing function in the Ordos in 2018
從時間變化(見圖3)上看,2000?2018年,鄂爾多斯市防風(fēng)固沙功能呈顯著增加趨勢. 單位面積防風(fēng)固沙量由1.62×104t/km2增至10.95×104t/ km2,以平均每年26.96%的變化率波動增加. 從不同階段分析,2000?2009年,單位面積防風(fēng)固沙量變化相對平緩,年均變化率為31.68%,期間防風(fēng)固沙量增加了1.19×104t/km2;2009?2018年增長幅度較大,年均變化率為63.03%,期間單位面積防風(fēng)固沙量增加了8.14×104t/km2.
圖 3 2000?2018年鄂爾多斯市防風(fēng)固沙功能的變化Fig.3 Annual changes of wind-breaking and sandfixing function in Ordos from 2000 to 2018
從空間變化(見圖4)上看,2000?2018年鄂爾多斯市有54.89%的區(qū)域防風(fēng)固沙功能增加. 其中,防風(fēng)固沙功能明顯增加區(qū)主要分布在準(zhǔn)格爾旗和烏審旗的東南部等區(qū)域,其面積占鄂爾多斯市總面積的18.41%;防風(fēng)固沙功能一般增加區(qū)主要分布在杭錦旗的庫布齊沙漠北部、烏審旗的毛烏素沙地、鄂托克前旗、達(dá)拉特旗、伊金霍洛旗、東勝區(qū)、康巴什區(qū)等區(qū)域,面積約占鄂爾多斯市總面積的36.48%. 此外,相比2000年,2018年鄂爾多斯市有15.33%的區(qū)域防風(fēng)固沙功能下降,主要分布在杭錦旗北部的庫布齊沙漠和鄂托克旗西部. 其余29.76%的防風(fēng)固沙功能沒有發(fā)生明顯變化.
圖 4 鄂爾多斯市2000—2018年防風(fēng)固沙功能的空間分布Fig.4 The spatial distribution pattern of wind-breaking and sand-fixation function in Ordos from 2000 to 2018
4.3.1降雨量對防風(fēng)固沙功能的影響
利用R語言編程進(jìn)行降雨量和防風(fēng)固沙功能的顯著性檢驗,結(jié)果表明年降雨量與防風(fēng)固沙功能存在顯著相關(guān)(P<0.01,r=0.69). 2000?2018年鄂爾多斯市的降水量在200.58~430.96 mm之間,多年平均降雨量為317.3 mm,呈現(xiàn)波動增加的趨勢,且與防風(fēng)固沙能力變化趨勢呈現(xiàn)較好的一致性. 從防風(fēng)固沙能力變化較大的年份分析,2012年降雨量較2011年大幅增加,2012年防風(fēng)固沙功能也顯著增加;而2017年降雨量下降,防風(fēng)固沙功能也小幅下降(見圖5).
圖 5 2000—2018年鄂爾多斯市防風(fēng)固沙功能和降雨量的變化Fig.5 Annual changes of rainfall and WBSF function in Ordos from 2000 to 2018
圖 6 鄂爾多斯市降雨量與防風(fēng)固沙功能相關(guān)系數(shù)的空間分布Fig.6 The Spatial pattern of correlation coefficients between rainfall and WBSF function in Ordos
空間上,78.59%的區(qū)域相關(guān)系數(shù)為0.4~0.9(見圖6),主要分布在降雨量≥200 mm區(qū)域,該區(qū)域為防風(fēng)固沙功能一般區(qū)及以上,而相關(guān)系數(shù)為負(fù)值的占比較小,主要分布在庫布齊沙漠,該區(qū)域為防風(fēng)固沙功能低值區(qū). 說明該區(qū)域降雨量對防風(fēng)固沙功能的變化產(chǎn)生了積極影響.
4.3.2土地利用對防風(fēng)固沙功能的影響
根據(jù)鄂爾多斯市土地利用類型轉(zhuǎn)移矩陣(見表1),2000?2018年,鄂爾多斯土地利用類型變化的主要特點是城市擴(kuò)張、林地和草地增加以及沙地、耕地與濕地資源減少. 林地、草地、建設(shè)用地表現(xiàn)出不同程度的增加,其增幅分別為11.03%、1.45%、163.18%,其中林地和草地面積的增加主要來源于沙地,沙地整體呈現(xiàn)出凈轉(zhuǎn)出的特征,建設(shè)用地的增加主要來源于草地. 耕地和濕地總面積減少,其減幅分別為2.09%和6.71%,其中耕地主要轉(zhuǎn)換為建設(shè)用地以及退耕還草.
表 1 鄂爾多斯市2000?2018年土地利用類型轉(zhuǎn)移矩陣Table 1 Transferring matrix of land use type in Ordos from 2000 to 2018 km2
通過對2000?2018年土地利用類型變化和防風(fēng)固沙功能變化的疊加(見圖7),可統(tǒng)計分析不同土地利用類型之間的轉(zhuǎn)化對防風(fēng)固沙功能的影響. 結(jié)果表明,鄂爾多斯市土地利用類型發(fā)生變化的面積占研究區(qū)總面積的25.62%,而在土地利用類型發(fā)生變化的區(qū)域,其防風(fēng)固沙量的變化占研究區(qū)總防風(fēng)固沙量變化的29.86%,說明人類活動主導(dǎo)下的土地利用類型變化對防風(fēng)固沙功能的影響不容忽視. 影響防風(fēng)固沙功能的用地類型主要在林地、草地與耕地和沙地之間進(jìn)行. 鄂爾多斯市近20年實施了大面積植樹造林、京津風(fēng)沙源治理工程以及禁牧、休牧等措施,2018年沙地面積較2000年縮減了1 378.56 km2,減幅為6.51%,沙地主要轉(zhuǎn)變?yōu)椴莸?、林地和耕地,這部分區(qū)域防風(fēng)固沙功能增加了35.00%. 草地和林地增加區(qū)域防風(fēng)固沙功能分別增加了20.38%和8.78%. 而耕地轉(zhuǎn)為林地,防風(fēng)固沙功能降低了2.47%,耕地向林地轉(zhuǎn)換大多是人為造林的結(jié)果,新造林地在生長初期植被覆蓋度較低,導(dǎo)致短期內(nèi)地表植被的保護(hù)作用下降,固沙功能輕度降低. 綜合來說,與生態(tài)環(huán)境質(zhì)量提高相關(guān)的土地利用類型變化對區(qū)域防風(fēng)固沙功能起到了一定的增強作用.
圖 7 2000—2018年鄂爾多斯市土地利用類型和防風(fēng)固沙功能的空間變化Fig.7 The spatial change of land use and WBSF function in Ordos from 2000 to 2018
該研究評估了2000?2018年鄂爾多斯市防風(fēng)固沙功能變化及長時間序列的空間差異,發(fā)現(xiàn)鄂爾多斯市防風(fēng)固沙功能明顯提升,這與張照營[38]和張燕婷[39]得到的內(nèi)蒙古自治區(qū)防沙屏障帶和北方防沙帶的防風(fēng)固沙能力明顯增強的結(jié)論均一致. 該研究植被覆蓋度的變化與邵艷瑩等[40-42]的研究結(jié)果一致. 該研究計算得出的防風(fēng)固沙能力稍高于江凌等[43]估算的2000?2010年內(nèi)蒙古自治區(qū)的年均防風(fēng)固沙能力(4 880 t/km2),原因在于評估時段和區(qū)域不同. 另外,筆者采用小時數(shù)據(jù)作為氣象因子計算的數(shù)據(jù)來源,時間分辨率更高.
RWEQ模型是以物理過程理論為基礎(chǔ)的經(jīng)驗?zāi)P蚚44],其計算參數(shù)具有較強的地域性,具體應(yīng)用時需進(jìn)行必要的參數(shù)修正. 但遲文峰等[45]利用內(nèi)蒙古高原137Cs示蹤技術(shù)監(jiān)測成果驗證RWEQ模型反演結(jié)果,發(fā)現(xiàn)二者具有較好的擬合性(R2=0.83,P<0.01),從側(cè)面說明RWEQ模型在我國北方區(qū)域具有可行性. 筆者所在課題組通過在伊金霍洛旗實地采樣,測量22個土壤剖面的土壤含水量來修正當(dāng)?shù)仄鹕筹L(fēng)速,其他參數(shù)在綜合考慮與借鑒前人研究成果的基礎(chǔ)上,采用了我國北方沙化地區(qū)修正后的關(guān)鍵參數(shù)和計算公式,并且為滿足模型對數(shù)據(jù)空間和時間分辨率的要求,對輸入?yún)?shù)進(jìn)行了插值、擬合、漸進(jìn)等優(yōu)化,提高了模型的精確度. 在今后的研究中,可進(jìn)一步對模型參數(shù)如土地粗糙度以及土壤結(jié)皮和可蝕性因子的土壤粒徑轉(zhuǎn)換等問題進(jìn)行本地化調(diào)整,構(gòu)建具有更廣泛適用性的土壤風(fēng)蝕模型,以便更準(zhǔn)確地反映區(qū)域防風(fēng)固沙功能的變化狀況,進(jìn)而引導(dǎo)防風(fēng)固沙生態(tài)治理工程有針對性地投入.
a) 2018年鄂爾多斯市防風(fēng)固沙總量為95.07×108t,單位面積防風(fēng)固沙量為10.95×104t/km2;防風(fēng)固沙功能明顯存在區(qū)域差異,整體呈現(xiàn)由東向西遞減的趨勢,高值區(qū)面積占比為25.14%,低值區(qū)面積占比為10.97%.
b) 2000?2018年,鄂爾多斯市防風(fēng)固沙功能以平均每年26.96%的變化率波動增加. 54.89%的區(qū)域防風(fēng)固沙功能增強,主要分布在準(zhǔn)格爾旗和烏審旗的東南部;29.76%的區(qū)域防風(fēng)固沙功能穩(wěn)定,主要位于杭錦旗、鄂托克旗、伊金霍洛旗的部分區(qū)域;而防風(fēng)固沙功能降低區(qū)域為15.33%,主要位于杭錦旗北部的庫布齊沙漠和鄂托克旗西部地區(qū).
c) 在驅(qū)動因素中,鄂爾多斯防風(fēng)固沙功能變化與降雨量呈顯著相關(guān)(P<0.01,R2=0.69),多年平均降雨量為317.3 mm,呈現(xiàn)波動增加的趨勢,且與防風(fēng)固沙能力的變化趨勢呈現(xiàn)較好的一致性,78.59%的區(qū)域二者相關(guān)系數(shù)為0.4~0.9. 土地利用類型變化中,沙地面積減幅為6.51%,其防風(fēng)固沙功能增加了35.00%;草地和林地增加區(qū)域防風(fēng)固沙功能分別增加了20.38%和8.78%,退耕還林還草、沙地改善等有利于生態(tài)環(huán)境質(zhì)量提高的土地利用類型轉(zhuǎn)化,對改善區(qū)域防風(fēng)固沙功能有積極作用.
d) 生態(tài)保護(hù)和恢復(fù)工程建設(shè)應(yīng)該綜合考慮氣候變化和人類活動的復(fù)合影響,發(fā)揮氣候變化帶來的正面效應(yīng),結(jié)合鄂爾多斯市的自然地理環(huán)境特點,注重生態(tài)工程布局的區(qū)域適宜性,對生態(tài)系統(tǒng)進(jìn)行科學(xué)的保護(hù)與修復(fù),持續(xù)發(fā)揮北方防沙帶生態(tài)安全屏障作用.