乡镇尺度大娄山区生态系统服务价值时空变化研究

姜栋栋， 马伟波， 邹凤丽， 李海东， 张龙江， 刘桂建

1.中国科学技术大学， 安徽 合肥 230026

2.生态环境部南京环境科学研究所， 江苏 南京 210042

3.生态环境部华东督察局， 江苏 南京 210019

4.武汉大学遥感信息工程学院， 湖北 武汉 430079

生态系统服务价值(Ecosystem Service Value，ESV)核算可定量表达生态系统获得的生命支持产品和服务[1]，是开展生态保护、自然资产核算和生态补偿的重要依据和基础，是区域生态安全的重要指标[2]，因此科学评估ESV可为区域发展政策的制定提供重要支撑[3]. 自从Costanza等[4]在1997年率先开展ESV核算研究以来，ESV引起普遍关注[5-9]. 以谢高地等[10]提出的中国ESV当量为代表，通过当量法测算ESV的相关研究被广泛探索，如ESV与土地利用响应[11-12]、ESV权衡与协同[13-14]、ESV驱动机制[15-23]等，研究成果为生态环境保护政策的制定提供了重要支撑作用[24-27]. 空间尺度是定量核算ESV的重要影响因素[25,28]，不同空间尺度下生态系统结构、过程和功能差异显著，较大尺度易忽视内部差异，而较小尺度则不能代表整体特征，因此不同空间尺度ESV响应特征不同，同时所对应的宏观政策调控也不同. 黄木易等[29]研究表明，乡镇是大别山区ESV空间分异的特征尺度，ESV在乡镇尺度信息量丰富.

在山区村镇建设中，地形因子对聚落分布、水热条件分配[30]、地质灾害[31]，尤其对生态系统的结构和功能[32-33]等有显著的约束和影响. 关于县域尺度ESV与地形梯度的研究指出，ESV在山区和平原有较大差异[34]，但对山区型乡镇尺度ESV空间分异特征分析以及地形效应还缺乏深入探索. 生态宜居是我国乡村振兴战略实施的主要目标之一，因此开展山区乡镇尺度的ESV及其地形梯度效应研究，分析不同乡镇ESV的功能差异和地形特征，解析山区型村镇建设生态约束条件，对山区型乡镇开展“三生”空间优化和生态服务功能提升具有重要意义.

为阐述影响村镇建设的生态约束条件，分析山区型乡镇尺度ESV地形效应特征，笔者基于2000年和2018年2期土地利用数据，通过ESV评估、地貌、地表粗糙度和空间统计等方法，分析了乡镇域ESV时空变化及其地形效应和空间关联特征，以期为大娄山区水源涵养与生物多样性保护功能区(简称“大娄山区”)[35]村镇建设以及开展生态安全格局优化和差异化功能调控政策的制定提供支撑.

1 材料与方法

1.1 研究区概况

大娄山区位于川、滇、黔、渝三省一市接壤地带(见图1)，地理位置为26°35′N～28°10′N、104°60′E～107°12′E，总面积3.29×104km2，属西南季风区，降雨主要集中在北部，年均降雨为 1 080 mm，大娄山区地势南高北低，海拔介于198～2 865 m之间. 大娄山区是赤水河与乌江水系、横江水系的分水岭以及重要的水源涵养区[30]，该地区中南部是典型的喀斯特山区，北部是丹霞地貌；受过度垦殖影响，局部地区水土流失严重. 大娄山区包括17个县市，涉及357个乡镇，其中乡镇行政区边界被完整包括的有293个. 2018年大娄山区内人口总计 1 368.9 万人，其中贵州省人口占比超过54%. 大娄山区中包含的四川省叙永县、古蔺县，贵州省大方县、赫章县、纳雍县以及云南镇雄县、威信县是我国集中连片特困区，经济发展水平较低. 大娄山区土地利用以林地为主，耕地和草地次之，2018年林地总面积为1.80×104km2.

图1 大娄山区水源涵养与生物多样性保护重要区地理位置

1.2 数据获取

该研究中高程数据采用SRTM (Shuttle Radar Topographic Mission)90 m的数字高程模型(Digital Elevation Data，DEM)，通过SRTM官网(http:srtm.csi.cgiar.orgsrtmdata)下载.

2000年和2018年2期土地利用数据通过中国科学院资源环境科学与数据中心获取，空间分辨率为30 m，整体分类精度超过85%. 将原水域二级分类中的滩涂和滩地划分为湿地类型，以满足ESV核算中对土地利用类型和生态系统的划分. 由表1可见，林地是大娄山区最主要的自然生态系统.

表1 2000年和2018年大娄山区土地利用情况

该研究采用MODIS传感器MOD173AH的植被净初级生产力指数(Net Primary Production，NPP)数据对区域性服务功能价值系数进行修正. NPP数据经镶嵌、投影转换、裁剪等预处理，单位为kgm2(以C计)，数据下载网址为https:e4ftl01.cr.usgs.govMOLTMOD17A3HGF.006.

粮食产量、播种面积、粮食价格等主要来源于《遵义市统计年鉴》《泸州市统计年鉴》《重庆市统计年鉴》和《中国农产品价格调查年鉴2018年》.

1.3 研究方法

该研究采用Costanza等[4]提出的当量因子法对大娄山区ESV进行估算，计算公式：

(1)

该研究通过地区农作物产值和面积校正地区ESV当量因子，以及通过NPP食物生产、原材料生产、气体调节、气候调节、净化环境、维持养分循环、维持生物多样性、提供美学景观等8种服务功能价值系数当量因子，无人力投资的自然生态系统输出价值为农田食物生产服务单价的17，标准当量修正公式：

(2)

式中：T为研究区粮食总价值，元，该研究价格水平参照《中国农产品价格调查年鉴2018年》；X为研究区粮食播种面积，hm2.

NPP条件因子计算公式：

λ=BB0

(3)

Ek=λ×E0k

(4)

式中：λ为生态系统服务当量的区域修正系数；B0为全国平均NPP数值，kgm2(以C计)；B为大娄山区NPP平均值，kgm2(以C计)；Ek为第k类生态系统经过区域修正后的生态系统服务当量；E0k为第k类生态系统全国生态系统服务当量平均值，k=1,2,…,8，分别代表食物生产、原材料生产、气体调节、气候调节、净化环境、维持养分循环、维持生物多样性、提供美学景观服务功能. 2000—2018年，遵义市、泸州市、毕节市和重庆市农作物种植面积为 1 543.65×104hm2，对应总产值为 2 848.83×108元，得到单位面积粮食产量的经济价值为 2 577.21元hm2. 2000—2018年全国NPP平均值为 4 327 kgm2(以C计)，大娄山区同期NPP平均值为 6 671 kgm2(以C计)，因此食物生产、原材料生产、气体调节、气候调节、净化环境、维持养分循环、维持生物多样性、提供美学景观功能修正系数均为1.54. 综上，得到大娄山区单位面积生态服务价值系数如表2所示.

表2 大娄山区单位面积生态服务价值系数

Table 2 ESV coefficient per unit area in Dalou Mountain area 元(hm2·a)

表2 大娄山区单位面积生态服务价值系数

一级类型二级类型耕地林地草地湿地裸地水域食物生产3 373.571 230.361 508.182 024.140.003 175.12供给服务原料生产1 587.562 817.922 222.591 984.450.00912.85水资源供给51.54953.57798.946 674.970.0021 365.07气体调节2 659.179 326.927 818.747 540.9279.383 056.06调节服务气候调节1 428.8127 901.3920 677.9914 288.050.009 088.79净化环境396.897 898.126 826.5114 288.05396.8922 027.41水文调节695.859 046.019 844.9462 445.8077.32263 493.95土壤保持2 654.537 370.826 185.305 953.3651.542 396.81支持服务维持养分循环476.27873.16714.40714.400.00277.82生物多样性515.9610 319.158 652.2131 235.2779.3810 120.70文化服务美学景观238.134 524.553 810.1518 772.9139.697 501.23

(5)

(6)

式中：m为R半径内DEM栅格的数量；zp为第p个栅格处高程，m. 鉴于研究区山地地形复杂的特点，通过3～100个像元半径单位的测试，最终选择R为30个像元尺度，即半径为 2 700 m的圆形卷积窗口，该尺度较好地兼顾了细分地貌与景观. 地貌的变化源于坡面的变化，根据研究区实际情况，改进了Dickson等[37]对地貌的划分方法，将大娄山区划分为山脊、陡坡、缓坡和峡谷底. 地貌划分综合了高程、slope (研究区绝对数值坡度)以及TPI所蕴含的地形特征(见表3).

表3 基于TPI的地貌划分

地表粗糙度定义为地表单元的曲面面积与其在水平面上的投影面积之比，粗糙度是反映地形起伏变化的宏观地形因子，研究区内所有坡度取正割计算后的平均值.

利用探索性空间数据分析大娄山区空间自相关系数，分析大娄山区ESV空间分布格局的空间集聚程度. 选择Moran′s I测度ESV全局空间自相关特征，Gi*测度ESV变化的聚集与分异特征，即“热点”与“冷点”分布格局[38]. 计算公式：

(7)

(8)

(9)

(10)

2 结果与讨论

2.1 ESV时空变化特征

2.1.1ESV时间尺度变化

由表4可见：2000—2018年大娄山区ESV整体呈稳定且微弱增长趋势，由 1 891.44×108元增至 1 895.73×108元，增幅为0.23%. 从生态系统服务功能增幅来看，水文调节ESV增幅(1.63%)最大，增加了3.44×108元；水资源供给ESV增幅(1.6%)次之；而食物生产ESV降幅为1.34%，ESV减少了0.88×108元. 从ESV数值变化来看，2000—2018年水文调节ESV增值最大，从211.27×108元增至214.72×108元，增加了3.45×108元；而食物生产ESV降幅最大. 从2018年各单项ESV占总ESV比例来看，气候调节ESV在2018年为592.46×108元，占大娄山区总ESV的31.25%；水资源供给占比(1.13%)最低，其ESV为21.36×108元. 2000—2018年大娄山区耕地面积减少了257.7 km2，导致食物生产和原材料生产服务功能降低，但林地、水域和湿地面积均增加，使得气候调节、水文调节和生物多样性的ESV均增加，尤其是湿地面积增加使水文调节ESV增加了3.44×108元.

表4 2000—2018年大娄山区ESV时序变化

2.1.2ESV空间分布变化

由图2(a)(b)可见，2018年大娄山区ESV较高的地区主要分布在大娄山区北部，与该地区以林地自然生态系统为主有关，主要涉及古蔺县、叙永县和赤水市. 大娄山北部地区分布有赤水桫椤国家级自然保护区、长江上游珍稀特有鱼类国家级自然保护区、赤水桫椤国家级自然保护区、四川画稿溪国家级自然保护区和遵义习水中亚热带常绿阔叶林国家级自然保护区，保护区面积达 1 450 km2，得益于保护区的建立，古蔺县、叙永县、赤水市、习水县等地区林地资源得到了较好的保护. 因此，该地区乡镇在净化环境、提供美学景观、气候调节和生物多样性和水资源供给等方面ESV整体高于大娄山区其他乡镇. ESV评估存在很多主观因素，不同评价方法在考虑视角、参数选择等方面的侧重，使得评估结果存在较大差异，但得出的变化趋势是一致的，不影响研究结果的相互比较. 王金凤等[39-40]对我国西南地区的ESV研究结果与笔者所得ESV分布趋势在整体上吻合.

图2 大娄山区ESV分布格局

将2018年ESV划分为5个区间，分别为0.3×106～3.0×106、3.0×106～4.5×106、4.5×106～5.5×106、5.5×106～6.5×106和6.5×106～8.8×106元km2〔见图2(b)〕，分布在5个区间内的乡镇个数分别为18、65、121、111和40个. 分布在ESV为6.5×106～8.8×106元km2区间的乡镇主要集中在2个地区，一处是以赤水市习水县为主体的川、黔交界处，包括元厚镇、东皇镇、黄荆乡、官渡镇、福宝镇、四面山镇等17个乡镇；另一处是汇川区、桐梓县和绥阳县相邻的大娄山东北部边缘地区，主要包括高坪镇、毛石镇、沙湾镇、茅石镇等12个乡镇. 分布在ESV为0.3×106～3.0×106元km2区间的乡镇分布比较分散，主要包括大方县顺德街道、洪山街道等大娄山西北部边缘乡镇.

将ESV增幅超过3%、降幅超过3%的地区分别划分为上升和下降地区，其余地区划分为稳定地区. 由图2(c)可见：2000—2018年大娄山区乡镇尺度ESV总体呈基本稳定态势，有309个乡镇ESV基本保持稳定，面积为2.8×104km2，占大娄山区总面积的85.74%；有28个乡镇ESV呈下降趋势，其中22个分布在贵州省，包括七星关区青场镇等11个地区，赤水市复兴镇等3个地区. 以云南省以勒镇、贵州省复兴镇和茅台镇为例，自2016年以来，当地均开展了旅游等经济开发建设活动，土地利用格局发生了较大变化，因而对其ESV造成显著影响. 而七星关区的清水铺镇位于省际毗邻地区，近些年追求经济发展，耕地、林地、草地和湿地土地利用转移频率较高，其中湿地转化为耕地开展农业发展是其土地利用格局变化的主要因素. 有20个乡镇ESV呈上升趋势，且乡镇分布较散，主要分布在贵州省(15个)，增幅最大的2个乡镇为重庆市的永新镇和云南省的赤水源镇，ESV增幅均超过10%. 云南省赤水源镇ESV上升的主要原因来自于赤水河下游经济发达地区对上游地区的补偿措施，促使其在保持较好的生态环境状况下开展以烟叶种植和高原经济作物为主的产业. 其他ESV上升乡镇在不同程度上保持了土地利用格局变化较小的特征.

2.1.3单项ESV的空间分布

由图3可见：净化环境ESV、提供美学景观ESV、气候调节ESV和生物多样性维护ESV虽然价值量差异较大(以气候调节ESV最高，提供美学景观ESV最低)，但分布态势相似；同时，食物生产ESV、维持养分循环ESV、原材料生产ESV分布格局相似，整体价值量水平较低，且与净化环境ESV、提供美学景观ESV、气候调节ESV和生物多样性ESV的分布态势相反. 土地利用格局是影响生态系统服务功能及其价值的重要因素，其中，净化环境ESV、提供美学景观ESV、气候调节ESV和生物多样性维护ESV较高的乡镇，其土地利用类型主要以林地为主；而食物生产ESV、维持养分循环ESV、原材料生产ESV较高的乡镇，其土地利用类型主要以耕地为主. 气体调节ESV和土壤保持ESV在空间上呈相似的分布格局.

图3 大娄山区2018年11种乡镇尺度ESV空间分布

由于大娄山区中叙永县等7个县是我国集中连片特困区，近年来为积极脱贫，在当地经济发展从低水平急剧增加的过程中，建设用地扩张迅猛，且土地利用在林地耕地和草地之间转换频率增加，ESV整体保值压力增大；同时，经实地调研，大娄山区内的贵州省水土保持功能整体较好，而云南省镇雄县和威信县的水土流失问题较为严重. 若这种发展趋势不被纠正，生态系统将无法支撑社会经济的可持续性发展，需结合生态系统特征和ESV格局制定公共服务政策[41].

2.2 ESV地形效应

2.2.1地形地貌

通过自然断点法将地表粗糙度值按[1.01～1.03)、[1.03～1.05)、[1.05～1.07)、[1.07～1.09)、[1.09～1.13)划分为5个等级，分别用Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ表示，5个地表粗糙度等级内分布的乡镇数分别为43、106、104、78和24个〔见图4(a)〕. 依据坡度和TPI结果，将大娄山区划分为4种地貌〔见图4(b)〕，结合表4可知，峡谷底、缓坡、陡坡和山脊区域的面积分别为 3 318.87、13 581.62、12 979.76 和 3 030.36 km2，以缓坡和陡坡面积占比最高，分别为41.27%和39.44%.

图4 大娄山区乡镇尺度地表粗糙度和地貌划分

2.2.2地貌效应

由表5可见：在大娄山区4种地貌类型中，以缓坡区域的ESV最高(2018年)，为811.69×108元；陡坡区域的ESV最低(2018年)，为158.65×108元. 2000—2018年缓坡区域ESV增加了17.15×108元，而陡坡区域ESV下降了15.92×108元，山脊和峡谷底ESV变化较小. 2018年，区域单位面积ESV由高到低依次为峡谷底(0.21×108元km2)>山脊(0.07×108元km2)>缓坡(0.06×108元km2)>陡坡(0.01×108元km2). 由于水域和湿地价值系数较高的土地类型主要分布在大娄山区峡谷底区域，所以峡谷底单位面积生态服务功能极高.

表5 大娄山区ESV地貌效应

2.2.3地表粗糙度效应

由表6可见：2018年大娄山区地表粗糙度以Ⅲ级ESV最高，为630.24×108元；以Ⅴ级区域的ESV最低，为141.11×108元. 2000—2018年地表粗糙度Ⅲ级区域ESV增加了5.68×108元，增幅为0.91%，而Ⅰ级区域ESV下降了3.28×108元，Ⅲ和Ⅱ级区域分别增加了5.68×108和2.47×108元，其中Ⅲ级区域增幅最大，达0.91%.

表6 大娄山区ESV地表粗糙度效应

去除大娄山区边缘面积低于10 km2的乡镇，绘制2018年乡镇地表粗糙度与ESV散点图以分析二者之间的线性关系. 由图5可见，大娄山区乡镇粗糙度与单位面积ESV呈线性相关，随着地表粗糙度的增加ESV呈增长趋势，二者之间的拟合优度R2为0.49.

图5 大娄山区乡镇尺度单位面积ESV与地表粗糙度散点图

大娄山区乡镇地表粗糙度与ESV拟合优度低于杨锁华等[34]以长江中游地区县域尺度地形起伏度与ESV的拟合优度(R2为0.53)，主要原因是长江中游地区范围超过50×104km2，其中山地平原等地理特征差异显著，因此地形起伏度指标设置合理；而笔者前期试验中地形起伏度对于乡镇尺度不够敏感，选择地表粗糙度则能在一定程度指示ESV的地形效应.

2.3 ESV空间统计

2.3.1空间自相关

由表7可见，2000年和2018年大娄山区ESV全局Moran′s I值分别为0.111和0.109，且P值均小于0.05，表明研究区有显著的空间自相关性和空间聚集效应. Moran′s I值为正值，但数值较小，表明研究区ESV的空间分布具有一定程度的正向自相关性，即大娄山区各乡镇ESV在空间上具有一定的聚集性，高值区趋于聚集，低值区趋于相邻. 2000—2018年全局Moran′s I值达到0.140，且P值小于0.05，说明在2000—2018年大娄山区各乡镇ESV的空间自相关性整体呈增强趋势，即2000—2018年大娄山区乡镇的社会经济发展导致ESV的空间分布聚集性不断增强.

表7 大娄山区ESV全局相关性

2.3.2ESV变化冷点与热点

由图6可见：2000年大娄山区乡镇尺度ESV分布有通过99%置信度的热点，主要分布在大方县的顺德街道和七星关区岔河镇等. 2018年，大娄山区乡镇尺度冷点热点呈一定程度的聚集，总体上有3个冷点集聚区和5个热点集聚区. 冷点集聚区主要分布在原2000年热点集聚区的周边地区，主要涉及六龙镇等16个乡镇，这些地区逐渐由不显著性冷点转换为显著性冷点；热点集聚区主要分布在后山镇、太平镇、榕右乡、柏林镇和双木镇为中心的地区，涉及34个乡镇. 2000—2018年，乡镇尺度ESV变化主要有5个空间热点聚集区和2个冷点聚集区，热点聚集区涉及36个乡镇，冷点聚集区涉及19个乡镇，且在大娄山区东北部以温水镇和高大坪乡为中心的热点聚集区呈一定的区域连通趋势，即大娄山区东北部地区的乡镇在2000—2018年ESV变化呈较强的空间聚集态势.

图6 2000—2018年大娄山区乡镇尺度ESV冷点和热点空间分布

3 结论

a) 2000—2018年大娄山区ESV整体呈稳定且微弱增长趋势，由 1 891.44×108元增至 1 895.73×108元；有28个乡镇ESV呈下降趋势，其中22个分布在贵州省.

b) 在大娄山区4种地貌类型中，以缓坡区域的ESV最高(2018年)，为811.69×108元；而区域单位面积ESV以峡谷底最高，为0.21×108元km2. 2000年和2018年大娄山区ESV全局Moran′s I值分别为0.111和0.109，且P值均小于0.05，表明研究区有显著的空间自相关性和空间聚集效应. 2000—2018年，乡镇尺度ESV变化主要有5个空间热点聚集区和2个冷点聚集区

c) 建议对大娄山区28个ESV下降的乡镇、大娄山区ESV冷点聚集区所涉及的乡镇、大娄山区Ⅲ级地表粗糙度和缓坡重叠区域应优先开展“三生”空间优化，实施“山水林田湖草”系统性生态保护修复工程，加强对湿地、水域和林地等高ESV的保护与恢复，维护和强化大娄山区水源涵养和生物多样性生态系统服务功能，引导土地利用向ESV保值或增值方向发展.

d) 该研究在ESV评估时仅分析了2000年和2018年2个时间节点，此期间逐年份ESV及其地形和空间特征的信息捕捉不足，后期应积极补充相关资料数据，开展时间序列分析，挖掘时间尺度上乡镇ESV特征信息，总结乡镇尺度ESV功能提升策略.