山地岩溶区生态系统服务时空演变特征分析<br/>——以贵州省为例

山地岩溶区生态系统服务时空演变特征分析
——以贵州省为例

2023-10-07徐梓津张雪松陈明曼

生态环境学报 2023年7期

徐梓津，张雪松，陈明曼

1.华中师范大学城市与环境科学学院，湖北武汉 430079；2.贵州财经大学研究生院，贵州贵阳 550025；3.贵州财经大学管理科学与工程学院，贵州贵阳 550025；

生态系统服务是指人类直接或间接从生态系统中得到的惠益（Daily et al.，2009），人类社会的可持续发展在很大程度上取决于生态系统服务的可持续供给（WRI，2007）。在经济高速发展的背景下，人类活动显著改变了生态系统宏观结构并使其服务功能呈现复杂的时空变化特征（Clericia et al.，2019）。联合国开展的千年生态系统评估（Millennium Ecosystem Assessment，MA）结果指出：全球约60%的生态系统服务功能退化（Millenniu Ecosyste Assessment，2005）且极大地威胁了区域乃至全球的生态安全，人类社会可持续发展面临严峻挑战。在生态系统极为脆弱的山地岩溶区该问题尤为突出。山地的地表形态起伏变化剧烈，区域内海拔落差较大。而岩溶地貌是地下水与地表水对可溶性岩石溶蚀与沉淀以及重力崩塌等作用形成的特殊地貌。岩溶区生态极为脆弱、地表破碎度高、地域分异明显，生态系统相比其他区域更易遭到破坏（史莎娜等，2022）。深入了解山地岩溶区生态系统服务时空变化特征，对正确认识山地岩溶区生态系统服务发展规律以及提升区域生态系统服务具有重要意义（Bennett et al.，2009；Tian et al.，2016）。

近年来，学者们开始关注岩溶地貌区生态系统服务。岩溶区生态系统服务研究的热点集中在中小尺度且以小流域、县域行政区划为主（张明阳，2011；王月容等，2013；徐建宁，2016；徐烨，2019）；从生态系统类型来看，岩溶区森林生态系统（李阳兵等，2005；林清山等，2010）和流域生态系统是主要的研究对象。从生态系统服务类型来看，生态系统服务评估类型以涵养水源、水土保持、碳储存和产品供给为主（吴孔运等，2008；韦钧培等，2022）；在评估方法上，运用价值量评价法和物质量评价法较多（林子雁等，2020）。经过梳理发现，山地岩溶区生态系统服务时空变化大多关注水平空间的时序变化，缺乏对纵向空间生态系统服务变化的研究。山地岩溶区代表了具有浓缩的环境梯度、高度异质化的生境（梁红柱等，2022）。山地岩溶区的垂直带是表征岩溶生态景观垂直变化的经典地学模型，也是地学分异变化研究的重点内容之一（Kou et al.，2020）。所以，对山地岩溶区生态系统服务在纵向空间的变化开展研究显得尤为重要。基于此，本文选取典型山地岩溶区贵州省作为研究区，针对山地岩溶区生态系统存在的关键问题，同时考虑研究数据的可量化性和权威性，基于InVEST模型定量评估2000－2020年产水量服务、土壤保持服务和生境质量，揭示其在水平空间和纵向空间的演变特征。以期加深对山地岩溶区生态系统服务的认识和理解，为山地岩溶区制定相应的生态管理决策提供参考，提高典型区域生态系统效益。

1 数据与方法

1.1 研究区概况

贵州省（103°36′－109°35′E、24°37′－29°13′N）位于中国西南部，东毗湖南、南邻广西、西连云南、北接四川和重庆（见图1）。全省国土面积约1.76×106km2。占全国国土面积的1.8%（Zhang et al.，2021）。境内地势西高东低，平均海拔1100 m。贵州岩溶地貌发育典型，岩溶地貌面积约占全省国土面积的70%。气候属亚热带湿润季风气候，年均气温10－18 ℃，年降水量1000－1500 mm（贵州省统计局，2021）。贵州省共有9个市、州，88个县。2020年地区生产总值17.8×103亿元，比上年增长4.5%，高于当年全国的GDP增速2.3%（贵州省统计局，2021）。2020年末全省常住人口3.86×103万人，比上年末增加10万人。开展对山地岩溶区生态系统服务的研究，贵州省具备代表性和典型性。

图1 研究区位置示意图Figure 1 Location diagram of study area

1.2 数据来源与处理

结合研究区实际面积大小兼顾软件处理的可操作性和高效性，将研究所涉及的数据均进行网格化处理，分辨率统一调整为1 km×1 km。同时，将论文所涉及到的空间数据地理坐标统一调整为GCS_WGS_1984坐标，投影坐标统一调整为Mercator投影。本文中水平空间指地表平面空间，纵向空间指海拔高于地表的空间。

1.2.1 土地利用数据

研究区土地利用类型数据采用中国科学院资源环境科学数据中心数据注册与出版系统的解译数据。该数据解译精度达到94.3%（徐新良等，2018）可以满足本文的研究需要。基于土地利用类型原始数据的土地利用分类，结合研究区实际和研究需要，主要根据土地资源及其利用属性，将研究区土地利用类型划分为7个类别，包括耕地、林地、草地、湿地、水域、建设用地和裸地。

1.2.2 土壤数据

土壤相关数据来源于综合世界土壤数据库（Harmonized World Soil Database，HWSD），该土壤数据因其属性全面被广泛使用于各类相关研究中，可以满足论文研究的需要。

1.2.3 气象数据

降水数据来源于中国逐月降水量数据集（Yang et al.，2021）。潜在蒸散发数据（Peng et al.，2019）来源于中国1 km逐月潜在蒸散发数据集（1990－2020），该数据集是基于中国1 km逐月均温、最低温、最高温数据集（彭守璋，2020），采用Hargreaves潜在蒸散发计算式得到（Ding et al.，2020；Ding et al.，2021），数据品质能够满足论文研究的需要。

1.2.4 DEM数据

高程数据来源于地理空间数据云网站（http://www.gscloud.cn/）提供的由ASTER GDEM V3数据加工生成的数据高程模型数据（Digital Elevation Model，DEM），空间分辨率为30 m×30 m。该高程模型数据被广泛应用于地理、生态、环境等领域的研究，数据精度能够满足论文研究的需要。

1.2.5 其他数据

降雨侵蚀力R、土壤可侵蚀性因子K均由经验公式计算（Wischmeier et al.，1965；Williams et al.，1997），主要威胁因子、威胁源因子权重和不同土地利用类型对威胁因子的敏感性由参考文献确定（Anselin，1995；钟莉娜等，2017；Zhang et al.，2020a，2020b；朱增云等，2020），其他输入数据按InVEST模型说明书默认或推荐值输入。

1.3 研究方法

1.3.1 产水量

本研究利用InVEST模型中的Water Yield模块对研究区生态系统服务中的产水量进行评估。具体计算公式如下：

式中：

Wij——指i栅格单元上土地利用类型j的水源涵养量；

Aij——j类土地利用类型i栅格单元上的年均实际蒸散量；

Pi——i栅格单元上的年均降雨量。

式中：

Pi——i栅格单元的潜在蒸散量；

w——i栅格单元上的可利用水量和年降水量的比值，为非物理参数，无量纲。

Kij——土地利用类型j的i栅格单元上植物蒸散系数；

Ei——i栅格单元上参考蒸散量，是基于当地参考植被的蒸散量。

1.3.2 土壤保持

土壤保持是生态系统服务功能的一个重要方面。在岩溶区，土壤保持可以反映区域生态系统服务功能的强弱以及区域石漠化程度。本研究中土壤保持服务功能采用InVEST模型中土壤侵蚀模块计算，计算公式如下：

式中：

D——年均土壤保持量；

R——降雨侵蚀因子；

K——土壤可侵蚀因子；

L——坡长坡度因子；

C——植被覆盖和管理因子；

Ps——保土措施因子。

1.3.3 生境质量

本文用生境质量指数表征生境质量。InVEST模型中的Habitat Quality模块通过分析区域土地利用类型及其对生物多样性威胁程度计算得到区域生境质量指数。生境质量指数计算包括4个函数：（1）每种威胁因子的相对影响；（2）每种土地利用类型对每种威胁因子的相对敏感性；（3）土地利用类型与威胁源之间的距离；（4）土地受到保护的程度（包玉斌等，2015）。

式中：

Qij——土地利用j中栅格单元i的生境质量指数；

Hj——生境类型j的生境适宜度，取值范围0－1；

k——半饱和常数，一般为生境退化度最大值的1/2；

z——归一化常量，设置为2.5（Gao et al.，2017）。

Dij、R、br和yr——生境退化指数、威胁因子个数、威胁因子r的权重、威胁因子的栅格数；

friy——生境与威胁源之间的距离及影响；

β——因公共政策减轻威胁对生境质量影响的因素（即法律保护程度，受法律保护的区域为0，其余区域为1）；

gjr——生境类型j对威胁因子r的敏感度；

diy——栅格i与栅格y的直线距离；

dr,max——威胁源r的最大威胁距离。

1.3.4 海拔等级划分

研究区最低海拔148 m，最高海拔2885 m，海拔落差高达2737 m，地形起伏变化剧烈。为了便于数据统计分析，也为了使图表展示更加清晰，本文将研究区海拔分为27个梯度。每个海拔梯度平均间隔约100 m，详细分级情况见表1，表中的海拔梯度范围包含小值，不包含大值，比如海拔梯度1级指海拔在148－300 m之间，包含148 m但不包含300 m的范围。

表1 研究区海拔分级及各海拔等级面积占比Table 1 Altitude classification and area proportion of each altitude level in the study area

2 结果分析

2.1 研究区土地利用类型宏观结构变化

由于土地利用逐年变化差异不大，因此本文用2000、2005、2010、2015、2020年共5期数据分析西部大开发之后，典型山地岩溶区贵州省20年的生态系统结构及服务变化。2000－2020年，研究区土地里利用类型中以林地为主，其次为耕地和草地，面积最小的土地利用类型为湿地（见表2）。其中建设用地、湿地、水域和林地面积增加，变化最为明显的是建设用地类型，变化率为371%，主要在2010年以后快速增加；湿地类型面积在2000－2020年间也显著增加，变化率为196%。其余土地利用类型面积均表现出减少的趋势，其中草地和裸地面积减少较多，变化率分别为-9.72%、-7.5%。耕地的面积变化率为-0.201%。

表2 研究区2000－2020年土地利用类型面积变化Table2 Land use type area change in the study area from 2000 to 2020 hm2

2.2 水平尺度演变特征

2.2.1 产水量变化特征

研究区产水服务虽然受到各种复杂因素共同影响但表现出与区域内降水量的分布情况高度一致。总体上看，研究区2000－2020年产水量变化呈现出减少-增加-减少的波动趋势。研究期内不同土地利用类型的产水量差异显著，除建设用地产水量增加外，其余土地利用类型的产水量均减少，其中林地的产水量减少最多为25.1×109m3。各土地利用类型单位面积产水量显示，裸地单位面积产水量最高，其次是建设用地单位面积产水量（见表3）。这是因为典型的裸露岩溶地貌占研究区总面积61.9%（贵州省人民政府，2021），区域内大部分裸地为裸岩石质地，该类地表多为不透水面，降水渗入率低，易形成径流；同时，裸地与建设用地植被覆盖度较低，水量蒸散较少，因此表现出裸地及建设用地单位面积产水量较高。而水域的单位面积产水量最低。这是因为虽然水体下渗程度小于耕地、林地、草地等土地利用类型，但由于其主要蒸散形式是水面蒸散，水面蒸散速度与强度远高于植被蒸散（赵晓松等，2013），高蒸散量导致水域表现出单位面积产水量最低。

表3 2000－2020年贵州省各土地利用类型产水量及单位面积产水量Table 3 Water production and unit area water production of various land use types in Guizhou province from 2000 to 2020

在2000－2020年间，研究区产水服务表现出明显的空间异质性。总体而言，产水量服务在研究区西部区域明显偏少（见图2）。该区域位于贵州省西部，森林覆盖率高而降水量低。森林由于较高密度的植物覆盖产生较高的蒸散发，较低的降雨量水平和较高的蒸散发水平，使该区域产水量服务表现出较低的水平。由于产水量受区域降水量的直接影响，在2000－2020年间产水量高值区并不固定。2000年产水量高值区主要集中在研究区南部，2005年产水量高值区主要集中在研究区北部边界、西南边界及东南边界处；2010年产水量高值区主要集中在研究区东北部边界处和西南部区域；2015年产水量高值区主要集中在研究区南部和东北边界区域；2020年产水量高值区主要集中在研究区中部小范围区域和北部区域。

图2 2000－2020年研究区产水功能时空分布Figure 2 Temporal and spatial distribution of water production function in the study area from 2000 to 2020

2.2.2 土壤保持变化特征

研究区耕地、林地和草地的土壤保持量变化趋势一致。2000－2020年研究区不同土地利用类型中林地的土壤保持量最高（见表4），因为森林的林冠层能截留降水、根系能固持土壤，同时，研究区林地的人类活动强度较低，两者对土壤均能起到保护作用。表现出林地的土壤实际侵蚀量较小，土壤保持量则较高。水域、湿地、裸地和建设用地的土壤保持量极少。研究区内由于岩溶地貌广泛分布，裸地多为裸岩石质地，对土壤没有固持能力。建设用地为研究区的人口聚集地，植被覆盖度低而人类活动强度高，不利于土壤保持。在2000－2020年间，研究区林地、耕地和草地的土壤保持量变化趋势一致，均呈现出降低-增加-降低的变化趋势，在2015年达到高峰。水域由于几乎没有植被和土壤覆盖，其对土壤保持服务的作用极为有限，所以在分析时一般不予考虑水域的土壤保持服务。

表4 研究期内贵州省不同土地利用类型土壤保持量及单位面积土壤保持量Table 4 Soil conservation of different land use types and soil conservation per unit area in Guizhou Province during the study period

研究区土壤保持量空间分布图中，每个像元值表示对应位置土壤保持的能力大小，该像元值并不是泥沙拦蓄的实际数量，而是表征当前土地覆被等条件下，理论上可避免的土壤流失量。总体上看，研究区土壤保持服务的空间分布具有明显的空间异质性，整个研究期内土壤保持量的最大值存在明显差异，但不同年份之间存在相似的空间分布格局。从研究区土壤保持服务时空分布图中可以看出，土壤保持量的低值区分布在中部、中北部人类活动强度高的区域。另外的低值区呈现出线条状，经与研究区海拔、流域数据对比，该线条状区域为海拔较低的峡谷径流区（见图3）。

图3 2000－2020年研究区土壤保持功能时空分布Figure 3 Temporal and spatial distribution of soil conservation function in the study area from 2000 to 2020

2.2.3 生境质量变化特征

本文采用生境质量指数表征生境质量的高低。2000－2020年研究区林地、耕地、草地、湿地与水域的生境质量指数较高（见图4），其中林地的生境质量指数最高，且表现出继续增高的趋势；湿地的生境指数也表现出继续升高的趋势；耕地、草地的生境质量指数表现出略微下降的趋势。而建设用地的生境质量指数最低，其次是裸地，两者的生境质量指数明显低于其他地类，且表现出继续下降的趋势。

图4 不同土地利用类型的生境质量指数Figure 4 Habitat quality index for different land use types

2000－2020年研究区生境指数变化平稳，2000－2015年未有较大波动，2015－2020年表现出显著提升（见图5）。生境指数的高值和低值呈现出相似的空间分布，空间异质性明显。通过研究区生境质量服务时空分布图可以看出，2000、2005、2010、2015和2020年低值区均出现在研究区中部，该区域为贵州省省会城市所在地，人类活动强度较高，人类活动干扰并影响区域内的生境质量。研究区生境指数在2015－2020年间显著提升，这与退耕还林、退耕还草、设立生态保护区等生态保护政策严格实施密切相关。

图5 2000－2020年研究区生境质量时空分布Figure 5 Temporal and spatial distribution of habitat quality in the study area from 2000 to 2020

2.3 不同海拔等级下生态系统服务功能演变特征

研究区不同海拔梯度分级的面积占比显示，海拔梯度4－12级共9个海拔梯度的面积各占全域面积的比例均大于5%，且面积合计占研究区全域面积的73.9%，显示出研究区约四分之三的陆地生态系统及资源分布在500－1400 m的海拔范围内。其中海拔梯度7级面积占比最高达到10.4%。

2.3.1 产水量变化特征

2000－2020年研究区单位面积产水量在不同海拔梯度的变化趋势一致，随海拔升高呈现出降低-升高-降低-再升高的波动变化态势，见图6。在此期间，不同海拔梯度单位面积产水量最高点均在海拔等级为1的区域，即海拔300 m以下的区域。2000－2020年研究区不同海拔梯度单位面积产水量最低值除2020年外均分布在海拔等级为25的区域，即海拔在2600－2700 m的范围内；2020年有所不同，单位面积产水量最低值分布在海拔等级为26的区域，即海拔在2700－2800 m的范围内。在海拔梯度等级2－11的区域，研究区单位面积产水量变化较为平稳；在海拔等级24－27的高海拔区域，研究区单位面积产水量变化波动较为剧烈。考虑上述变化与不同海拔梯度的降水量变化、植被覆盖变化密切相关。

图6 2000－2020年研究区不同海拔等级单位面积产水量Figure 6 Water yield per unit area at different altitudes in the study area from 2000 to 2020

2.3.2 土壤保持变化特征

2000－2020年研究区单位面积土壤保持量在不同海拔梯度的变化趋势一致，随海拔升高呈现出波动式下降的变化态势。在此期间，研究区不同海拔梯度单位面积土壤保持量最高点均在海拔等级为2的区域。研究区不同海拔梯度单位面积土壤保持量最低值除均分布在海拔等级为27的区域（见图7）。虽然在不同海拔级别研究区单位面积土壤保持量呈现出波动式变化，但总体变化趋势为随海拔升高，单位面积土壤保持量降低。

图7 研究期贵州省不同海拔等级单位面积土壤保持量Figure 7 Soil conservation per unit area at different altitudes in Guizhou Province during the study period

2.3.3 生境质量变化特征

2000－2020年研究区不同海拔梯度的生境质量指数变化趋势较为一致，随海拔升高均呈现出浅U型的分布特征（见图8）。其中生境质量最低值均出现在海拔分级16的区域；最高值均出现在海拔分级27的区域。此变化特点与其他区域生境质量变化特点差异较大，这是由研究区山地地貌特点决定的。研究区地形起伏剧烈，高低落差大，生境指数低值区距威胁因子较近，受到干扰较大，表现出生境质量较差；高海拔和低海拔区距威胁因子距离增加，干扰强度减弱，因此生境质量指数也逐渐增加。

图8 2000－2020年研究区不同海拔等级生境质量Figure 8 Habitat quality at different altitudes in the study area from 2000 to 2020

3 讨论

从水平空间及纵向空间对山地岩溶区开展长时序的生态系统服务演变研究不仅能够加深对特殊区域的生态系统发展演变规律的了解，对区域可持续发展也具有重要意义。山地是浓缩的环境梯度，探索山地生态系统服务在纵向空间的变化，能够帮助我们更全面的了解山地生态系统发展变化规律。本文从水平空间和纵向空间对山地岩溶区贵州省生态系统服务变化开展研究，更为立体的了解山地岩溶区生态系统服务变化特点。

从水平空间看，不同土地利用类型的产水量差异显著。林地由于植被覆盖度高，能够截留降水，所以表现出单位面积产水量水平较低。该特点使得林地地表径流减少，水土流失减弱（熊康宁等，2012）。而裸地与建设用地在自然或人为因素作用下形成地表不透水层，表现出单位面积产水量较高（Li et al.，2021）。这导致大量地表径流从裸露地表向土壤覆盖区冲刷，加重山地岩溶区水土流失，形成恶性循环。研究区土壤保持服务的空间分布具有明显的空间异质性，但不同年份之间存在相似的空间分布格局。2000－2020年间研究区不同土地利用类型中林地的土壤保持量最高，因为森林的根系对土壤具有固持能力，能对土壤起到保护作用。研究区裸地和建设用地的土壤保持量极少。研究区裸地多为裸岩石质地，与建设用地一样对土壤没有固持能力。因此，研究区在生态保护和石漠化防治过程中应重点关注岩溶裸露区域的生态修复，增加该类地表的植被覆盖度；城镇建设用地则需要提升绿化率，同时可推进“海绵城市”建设（高江波等，2022），目的是降低不透水地表的面积，减缓该类地表的水土流失。

从纵向空间看，产水量、土壤保持和生境质量在不同海拔等级具有显著的空间异质性。这与高江波等（2019）的研究结果一致。研究区单位面积产水量与土壤保持随海拔升高而波动式降低，由于气温和降水在不同海拔区域会有显著差别，从而改变了不同海拔区域植被种类和数量，进而影响产水量及土壤保持。研究区生境质量随海拔上升呈现出浅U型的分布特征，这是由于研究区人类活动主要位于海拔1100－1500 m之间，正好处于“浅U型”的底部即生境质量的低值区，该海拔区域距威胁因子较近，景观破碎化严重且易受影响干扰，表现出生境质量较差。而在低海拔和高海拔区域，随着人类活动强度降低且距威胁因子距离增加，生境质量受到的干扰强度减弱，因此生境质量指数也逐渐增加。所以，研究区应依据不同垂直带的生态系统服务分布特点制定差异化的生态保护政策，使生态系统服务功能在不同垂直带均发挥最大效益。

本文采用较为成熟的InVEST模型评估研究区生态系统服务水平，虽然从水平尺度和纵向尺度对山地岩溶区生态系统服务长时序变化开展研究。但实际观测数据支撑有限，评估结果存在一定的误差。未来条件允许时将进一步增加野外观测工作，获取并补充实测数据以支撑研究结果。