基于遙感RMMF模型的喀斯特地區水土保持功能評價

宮 晨,吳文瑾,段怡如,劉海江,何金軍,孫 聰,蔣 倩

1 中國科學院空天信息創新研究院,北京 100094

2 三亞中科遙感研究所,三亞 572029

3 中國環境監測總站,國家環境保護環境監測質量控制重點實驗室,北京 100012

4 鄂爾多斯市林業和草原科學研究所,鄂爾多斯 017010

5 北京航天長城衛星導航科技有限公司,北京 100032

水土流失是我國生態保護的重點關注問題,具有面積大、分布廣、治理難等特點。在我國的黃土高原、東北黑土區、長江經濟帶、桂黔滇等地區,水土流失問題持續突出[1]。水利部2018年發布的水土流失動態監測結果顯示,全國水土流失面積達273.69萬km2,占全國國土面積的28.6%,中度及以上水土流失面積為105.44萬km2,占總面積的38.5%[2]。為治理水土流失問題,我國劃定了水土保持重點生態功能區,并通過生態補償制度來獎勵各地在生態保護方面所做的貢獻,促進生態資源的高效利用,對提升區域生態環境狀況、推動可持續發展起到巨大作用。作為生態系統的重要能力指標之一,對水土保持功能進行有效的評價對生態效益補償工作的開展具有重要意義。

目前,國內外學者針對水土保持評價指標及模型已經開展了豐富的研究工作,通過對已調研水土保持評價方法進行歸納總結,可分為以下三類：

(1)土壤模型法

土壤模型法主要分為經驗統計模型和物理模型兩種。經驗模型通過對觀測數據的分析來尋求數據之間的特征關系,因其計算方便、結構簡潔,成為流域侵蝕研究的廣泛采用方法,能適用于觀測數據相對粗略的情況。物理模型是在明確侵蝕產沙機制前提下,利用沉積物質量守恒和流量動能等式,建立侵蝕模型,需要大量的觀測數據來支持。國際上流行的經驗統計模型包括RUSLE (Revised Universal Soil Loss Equation) 模型[3]、WEPP (USDA-Water Erosion Prediction Project) 模型[4]、PRISM (Probabilistic Symbolic Model Checker) 模型[5]等。其中最常見的USLE模型[6]是通過獲取年降水量、土壤侵蝕度、土地覆被和地形信息來估算土壤流失的年平均量,在世界范圍內得到了極為廣泛的應用。物理模型包括ANSWERS(Areal Nonpoint Source Watershed Environmental Response Simulation)模型[7]、CLEAMS(Chemicals, Loading, and Erosion from Agricultural Management Systems)模型[8]、MMF(Morgan-Morgan-Finney)模型[9]等。典型的LISEM (Limburg Soil Erosion Model) 模型[10]綜合考慮了土壤侵蝕的各個環節,將流域在空間離散化為一系列大小相等的柵格單元,對降雨侵蝕過程以等時間間隔分割,按照時間步長分時段模擬侵蝕過程。

(2)指標分析法

基于指標分析的生態功能評價方法通過因子加權記分綜合反映生態功能,該方法將得到的評估指標及權重排序構建研究區的生態功能評價模型,通過各項指標加權平均的綜合值,進行分級評價。常用的方法包括頻度分析法、專家咨詢法、層次分析法、灰色聚類分析法等。其中最為常見的層次分析法 (Analytic Hierarchy Process, AHP)[11]是將與決策有關的元素分解成目標、準則、方案等層次,在此基礎之上進行定性和定量分析的決策方法。

(3)回歸模型法

回歸模型法通過建立生態功能和其決定因素的回歸模型對生態功能進行評價,實現簡單,實用性強。例如植被保土作用系數C[12],C=0.779+0.595logC',式中,C'表示植被覆蓋度,可以反映植被抑制土壤侵蝕的能力；植被覆蓋管理因子是根據植被覆蓋度的不同,基于線性回歸模型分級得到。

目前我國大范圍的生態系統水土保持功能評價較多采用因子權重法,但這類方法主觀性較強,評價結果的物理意義也不夠明確,合理性和公平性易受到質疑。而基于土壤模型的方法難以區分自然氣候波動引起的水土流失量改變和生態系統水土保持能力變化引起的水土流失量改變,且在大范圍應用中面臨參數難以獲取的問題。特別地,現有評價中常采用的USLE模型并不適用于陡坡條件[13],給喀斯特等特殊地貌的評價帶來困難。綜合這些考慮,選取了可用于小流域、陡坡地貌的半物理土壤侵蝕模型RMMF[14],通過引入一系列優選遙感因子及其與模型輸入物理量之間經驗公式建立了遙感RMMF模型(RS-RMMF),基于這一模型建立了3個不受氣象動態干擾的水土保持功能評價指標,從而在兼顧模型評價客觀性和因子獲取便捷性的基礎上,對傳統評價工作中難度較大的喀斯特功能區進行了生態系統水土保持功能評價與時序分析。

1 研究方法

1.1 RMMF模型

RMMF土壤侵蝕模型通過將濺流和徑流的土壤剝離率與徑流輸送能力進行比較,可估算年土壤侵蝕速率并進一步實現對土壤侵蝕量的估計[15—16]。模型主要分為兩個階段模擬土壤侵蝕：計算降雨能量和地表流量的水相階段和計算土壤剝離和遷移率的沉積物階段。土壤侵蝕的年度總值由土壤剝蝕的年總速率(包括飛濺和徑流造成的土壤侵蝕)與徑流輸送能力之間的最小值給出。如圖1所示。

E=min{(F+H),TC}

(1)

式中,E為年度侵蝕總量(Mg hm-2a-1)；F是飛濺造成的土壤顆粒分離(Mg hm-2a-1)；H為徑流造成的土壤顆粒分離 (Mg hm-2a-1)；TC是徑流輸送能力(Mg hm-2a-1)。其中生態系統主要通過降水截流、減少徑流沖蝕量和運輸量來起到保持水土的作用,對相關參量的估算可以有效評估生態系統的水土保持功能。

1.2 RS-RMMF模型及指標計算

由于水土流失量受到降水動態的影響很大,因此水土流失量的大小并不與水土保持功能的強弱絕對相關。為了在水土保持生態功能評價中排除降雨的影響,基于RS-RMMF模型構建了3個新的指標對生態系統的水土保持功能進行評價,分別為區域單位降水截留率、區域單位徑流沖蝕量、以及區域單位徑流運輸量[17]。這3個指標計算所涉及的參數除土壤屬性外均可通過遙感手段獲取,并且土壤屬性變化極為緩慢、在短時序計算可認為近似恒定,這就使得模型可完全基于遙感數據實現動態評價。

(1)區域單位降水截留率PI0

植被覆蓋對單位面積降雨量的永久截留的比例稱為降水截留率,主要考慮植被覆蓋度及冠層垂直結構動態對降水截流量的影響。

(2)

IFmax=0.935+0.498LAI-0.00575(LAI)2

(3)

式中,PI0為植被降雨截留率,單位為%；CC為植被覆蓋度,單位為%；IFmax為從LAI估算的最大林冠存儲力,單位為mm；LAI為葉面積指數。

植被覆蓋度CC使用增強植被指數EVI計算得到,具體公式如下：

(4)

式中,EVImin表示理想無植被地表的EVI值,而EVImax則表示理想植被全覆蓋地表的EVI值。

(2)區域單位徑流沖蝕量H0

徑流沖蝕部分主要考慮植被覆蓋度和基巖裸露率變化對徑流沖蝕作用的影響及地形對泥沙的阻攔沉積作用,由土粒分散率、基巖裸露率、植被覆蓋度等得到沖量。

H0=(∑DRi×i)(sinS)0.3(1-CC)(1-EBC)×10-3

(5)

式中：H0為單位徑流沖蝕量,代表每mm徑流深沖刷產生的泥沙重量,單位為Mg/hm2；DR為土粒分散率；i分別表示粘粒、粉粒、砂粒；S為坡度因子,單位為°；CC為植被覆蓋度,單位為%；EBC為基巖裸露率,單位為%。

根據張莉等[18]在室內沖刷槽試驗實測數據得到不同類型土壤土粒分散率,結合各類型土壤粘粒、粉粒、砂粒含量反解出各級粒徑的土粒分散率,確定粘粒分散率為0.7,粉粒分散率為1.2,砂粒分散率為0.9。

坡度因子S采用DEM數據計算,具體公式如下：

(6)

式中,θ為地形坡度。

基巖裸露率EBC采用Huang等[19]模仿歸一化植被指數提出的歸一化巖石指數和張曉侖等[19]基于NDRI像元二分模型提取基巖裸露率的波段運算方法,其計算公式為：

(7)

(8)

式中,EBC為基巖裸露率；SWIR表示短紅外波段；R表示紅光波段；NDRI為歸一化巖石指數；NDRImax為歸一化巖石指數的最大值；NDRImin為歸一化巖石指數最小值。

(3)區域單位徑流運輸量TC0

主要考慮土地覆蓋/利用變化對徑流運輸作用的影響。

TC0=(∑DRi×i)×C×P×sinS×10-2

(9)

式中,TC0為每mm徑流深產生的單位徑流運輸量,單位為Mg/hm2；DR為土粒分散率；i分別表示粘粒、粉粒、砂粒；S為坡度因子,單位為°；C和P分別為植被覆蓋與管理因子和水土保持措施因子。

植被覆蓋管理因子C主要是基于植被覆蓋度CC得到最終的空間分布圖,使用的為蔡崇法等[21]建立的植被覆蓋度與C值之間的回歸關系。

(10)

對于非喀斯特地區,參考白雷超等[22]及黃杰等人[23]的研究,水土保持措施因子的賦值如表1所示。

表1 非喀斯特地區不同土地利用類型P值Table 1 P value of different land use types in non-karst areas

喀斯特地區屬于巖溶環境為主的特殊生態系統,區內土層瘠薄,降雨強度大,陡坡耕種普遍,水土流失非常嚴重。根據蔡崇法等人[20]的研究對P因子進行賦值,確定不同土地類型的P值,如表2所示。

表2 喀斯特地區不同土地利用類型P值Table 2 P value of different land use types in karst areas

1.3 生態系統水土保持功能綜合評價計算

為對區域生態系統水土保持功能進行綜合評估,將降水截留率PI0、徑流沖蝕量H0和徑流運輸量TC0三項指標歸一化為0-100之間無量綱數據,按照加權打分綜合得到區域水土保持能力得分。其中降水截留率PI0為正向指標,截留率值越大,水土保持效果越好；徑流沖蝕量H0和徑流運輸量TC0為負向指標,值越小,水土保持效果越好。由于生態系統對土壤的保持包含沖蝕截留和運輸截留兩部分,我們對兩個指標進行了整合,得到土壤截留能力得分。最后基于降水截留(保水)和土壤截留(保土)2項能力得分獲得水土保持的總評分。具體計算過程如下：

(1)對于正向指標歸一化處理：

(11)

(2)對于負向指標歸一化處理：

(12)

(3)截留能力計算：

土壤截留能力得分 = max(H0,TC0)

(13)

降水截留能力得分 = PI0

(14)

(4)綜合打分評估：

生態系統水土保持能力 = 0.5×土壤截留能力得分+0.5×降水截留能力得分

(15)

2 研究區概況

本研究選擇《全國主體功能區規劃》[24]劃定的9個水土保持功能區之一的桂黔滇喀斯特石漠化防治生態功能區(圖2)作為喀斯特地區的評價示例,評估其2011年至2019年生態系統水土保持功能狀況及變化特征。該功能區以貴州高原為中心,包括廣西西北部、云南東部和四川、重慶、湖南的部分地區。屬于巖溶環境為主的特殊生態系統,生態脆弱性極高,土地石漠化極為嚴重,土壤一旦流失,生態恢復難度極大。區內土層瘠薄,降雨強度大,陡坡耕種普遍,水土流失非常嚴重。

圖2 桂黔滇喀斯特水土保持區的地理位置Fig.2 Location of karst functional areas in Guangxi, Guizhou and Yunnan

3 結果與分析

3.1 數據來源

本研究中RS-RMMF模型輸入參量主要是通過國內外權威土壤和遙感數據集進行計算,數據來源如表3所示,其中的經驗系數結合參考文獻得到。

表3 輸入參量的數據來源Table 3 Data source of each input parameter

3.2 水土保持功能評價結果

首先基于RS-RMMF模型計算了桂黔滇喀斯特水土保持區的3個生態系統水土保持功能評價指標,即年度平均單位降水截留率、徑流沖蝕量和徑流運輸量。圖3—圖5展示了該功能區在2011年和2019年3個指標的空間分布情況,圖6和圖7分別為功能區的EVI和LAI空間分布圖,可以看到功能區東部和中部植被覆蓋較好,具有較高的降水截留能力,而西北部地表裸露情況較多。相比2011年,2019年功能區西北部植被指數和葉面積指數增加,降水截留率有較為明顯的提升,說明生態系統保水能力增強,可由植被垂直覆蓋度提升引起。在保土方面,功能區整體存在較多的土壤流失,相比2011年,2019年單位徑流沖蝕下降明顯,而運輸量也呈小幅下降,說明同等水量沖刷下土壤流失減少,生態系統的保土能力提升,可由植被水平覆蓋度提升、基巖裸露率下降引起。但可看出區域的保土能力動態在空間上存在較大的異質性和不連續性。功能區全區平均降水截留率由2011年的3.94%增加到2019年的5.58%；單位徑流沖蝕量由2011年的15.36×10-4Mg/hm2減少到2019年的9.40×10-4Mg/hm2；單位徑流運輸量由2011年的7.33×10-5Mg/hm2下降到7.27×10-5Mg/hm2。整體看來生態系統降水截留和沖蝕阻力大幅提升,運輸阻力有輕微提升,說明生態系統在保水、保土方面能力都有所增強,與空間分析結論一致。

圖3 喀斯特水土保持區2011和2019年平均單位降水截留率PI0Fig.3 Average rainfall interception rate PI0 of karst functional areas in 2011 and 2019

圖4 喀斯特水土保持區2011和2019年平均單位徑流沖蝕量H0Fig.4 Average runoff erosion amount H0 of karst functional areas in 2011 and 2019

圖5 喀斯特水土保持區2011和2019年平均單位徑流運輸量TC0Fig.5 Average runoff transportation volume TC0 of karst functional areas in 2011 and 2019

圖6 喀斯特水土保持區2011和2019年增強型植被指數EVIFig.6 Enhanced vegetation index EVI of karst functional areas in 2011 and 2019

圖7 喀斯特水土保持區2011和2019年葉面積指數LAIFig.7 Leaf area index LAI of karst functional areas in 2011 and 2019

進一步依照公式(15)計算了功能區在2011—2019年間生態系統水土保持功能的綜合得分,并采用自然分段法將總評分劃分6個區間,進一步計算了各區間的面積比例。圖8展示了2011和2019兩個年份的總評分空間分布。可以看到,得分主要集中在45—55區間,得分較高的區域主要分布于東部,而中部區域在2019年得分大幅增加。其中2011年水土保持功能綜合得分為49.75,2019年綜合得分為50.58。相比2011年,2019年綜合得分低于51的區域面積比例減少17.48%,得分在51—54區間的面積比例增加18.04%。圖9展示了2011—2019年總評分的統計分布和平均值的時序變化情況。結果顯示,該區域水土保持功能得分呈波動變化,線性趨勢不顯著。其中2012年的區域平均得分最低,為49.58,2017年的區域得分最高為52.92。 但從各等級分布的面積比例來看,從2011—2019年,得分處于45—48區間的區域比例呈線性下降,而得分處于51—54區間的區域比例線性增加,說明區域生態系統水土保持功能評分整體向高等級移動,生態功能向好,而區域均值中并未體現這一現象,可能由極大/極小值區域的評分波動引起。

圖8 喀斯特水土保持區2011和2019年水土保持功能綜合得分Fig.8 Comprehensive score of soil and water conservation function of karst functional areas in 2011 and 2019

圖9 2011—2019年喀斯特地區水土保持功能綜合得分變化值Fig.9 The comprehensive score changes of soil and water conservation function in karst areas from 2011 to 2019

3.3 與RUSLE結果對比分析

為了對比RS-RMMF模型相較于傳統經驗模型的評價差異,本文同時基于修正的通用土壤流失方程RUSLE[25-26],計算了桂黔滇喀斯特水土保持區的土壤侵蝕量,該模型計算公式為：

RUSLE=R×K×LS×C×P

(16)

式中,R為降雨侵蝕力因子,K為土壤可蝕性因子,LS,C和P分別為坡長-坡度因子,植被覆蓋與管理因子和水土保持措施因子,具體計算過程和數據來源與RS-RMMF模型一致。降雨侵蝕力因子所用數據是Terra Climate 0.25°分辨率月降水數據,土壤可蝕性因子所用數據同樣來自中國1:100萬土壤圖以及所附的土壤屬性表。

經RUSLE模型計算得到2011和2019年的土壤侵蝕分布情況,如圖8所示。2011年平均侵蝕模數為215.74 t hm-2a-1,2019年為209.37 t hm-2a-1,侵蝕模數下降2.95%。根據國家水利部2008年發布的SL190-207《土壤侵蝕分類分級標準》[27]對研究區土壤侵蝕強度進行分級,2011年微度、輕度、中度、強烈、極強烈、劇烈侵蝕面積分別占總面積的40.51%,2.33%,3.60%,4.19%,9.44%和39.93%；2019年的侵蝕分級面積比例分別為42.29%,2.12%,3.47%,4.04%,9.14%和38.94%。2011—2019年間功能區微度侵蝕面積比例增加,其余侵蝕面積比例減少,區域土壤流失量整體有所下降,與由RS-RMMF評價得出的結論基本一致。

但對比圖8和圖10可以看到,兩者在空間相對強弱的分布上存在顯著差異。其中,RS-RMMF模型的評價得分空間分布明顯由植被屬性主導,而在RUSLE的計算結果中,位于西北部的區域土壤侵蝕模數較低,但該區域植被生長并不好,而是由裸巖較多、土壤覆蓋少所導致的流失量較小。此外,RS-RMMF模型的評價與降水無關,而對比降水數據可知,RUSLE的結果在空間分布上很大程度受降水主導,如圖11所示,東南部區域降雨侵蝕力明顯高于西北部區域,土壤侵蝕模數也表現出西北部強于東南部的空間分布規律。同時RUSLE模型結果具有破碎化的景觀特征,由于喀斯特地區地勢起伏較為不規則,坡度-坡長因子波動較大,在陡坡等地形變化顯著區域坡度-坡長因子誤差也相應較大,進一步說明RUSLE模型并不適用于陡坡條件[28]。由此可見雖然兩者在2011年和2019年的對比結論基本一致,但受限于模型的計算機理,RUSLE的評價結果無論在時間上還是空間上都不是區域生態系統生態功能的絕對表征,難以被用于客觀評價人類對生態系統的保護和管理效果,從而并不適于被作為生態效益補償的計算指標[29—31]。與此同時,RS-RMMF模型不僅從降水截留、徑流沖蝕、運輸等物理過程對生態系統的水土保持功能做了詳細分析,從而可為有關部門有針對性地管理區域生態系統、提升其生態功能提供科學支撐,同時還分別評價了生態系統的保水和保土能力,有利于進一步對這些能力的價值進行精準核算。

圖10 喀斯特水土保持區2011和2019年土壤侵蝕模數分布圖Fig.10 Soil erosion modulus of karst functional areas in 2011 and 2019

圖11 喀斯特水土保持區2011和2019年降雨侵蝕力因子RFig.11 Rainfall erosion factor of karst functional areas in 2011 and 2019

4 結論與討論

喀斯特地區受長期溶蝕作用,土壤流失受植被覆蓋、土壤性質、徑流速率等影響,該區域水土流失具有隱蔽性和復雜性的特點,適用于大范圍、緩坡條件的USLE等模型難以做到準確評價。因此針對小流域、陡坡的喀斯特地貌環境,本研究基于RMMF模型,通過引入一系列優選遙感因子,與模型輸入物理量之間經驗公式建立了遙感RMMF模型。為在評估中客觀反映生態系統的水土保持功能,排除降水動態影響,基于RS-RMMF模型構建了區域單位降水截留率、區域單位徑流沖蝕量以及區域單位徑流運輸量3項指標,并通過綜合打分方法定量化計算區域水土保持功能。由RS-RMMF模型計算結果得知,相比2011年,2019年喀斯特功能區的區域單位降水截留率PI0升高了1.94%,區域單位徑流沖蝕量H0下降了5.96×10-4Mg/hm2,區域單位徑流運輸量TC0下降了6.0×10-7Mg/hm2,水土保持功能綜合得分增加0.83。該區域在此期間降水截留能力得到提升,降水截流能力不均衡性增加,單位沖蝕量向低等級移動,水土保持功能整體上得到一定程度改善。

與RUSLE土壤侵蝕模型比較發現,雖然兩者在2011年和2019年的對比結論基本一致,但在空間相對強弱的分布上存在顯著差異。其中RS-RMMF模型的評價得分空間分布由植被屬性主導,且受降水等因素影響較小。受限于模型機理,RUSLE的評價結果難以做到區域生態系統生態功能的絕對表征,不能客觀評價人類對生態系統的保護和管理效果,但RS-RMMF模型更多考慮了侵蝕過程中產沙流沙的物理機制,對土壤侵蝕描述更加清晰,同時通過遙感化改造克服了土壤侵蝕數據難收集的不足,從而可為有關部門客觀核算生態系統的水土保持能力提供科學支撐,進一步推進生態補償機制在各地的實施應用,促進生態資源的高效利用,推動生態環境可持續發展。

本研究目前僅在計算單位徑流運輸量時考慮了土地利用變化,暫未詳細區分自然條件和人為活動對水土保持生態功能的影響。下一步可以擴展加入社會經濟要素,使水土保持評價模型更加具有綜合性、動態性和可調控性,更能有針對性地對區域水土流失治理提出合理的建設建議。