基于遥感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模型更多考虑了侵蚀过程中产沙流沙的物理机制,对土壤侵蚀描述更加清晰,同时通过遥感化改造克服了土壤侵蚀数据难收集的不足,从而可为有关部门客观核算生态系统的水土保持能力提供科学支撑,进一步推进生态补偿机制在各地的实施应用,促进生态资源的高效利用,推动生态环境可持续发展。

本研究目前仅在计算单位径流运输量时考虑了土地利用变化,暂未详细区分自然条件和人为活动对水土保持生态功能的影响。下一步可以扩展加入社会经济要素,使水土保持评价模型更加具有综合性、动态性和可调控性,更能有针对性地对区域水土流失治理提出合理的建设建议。