南方丘陵山地带作为《全国主体功能区域》和“十二五”规划确定“两屏三带”为主体的生态安全战略格局的重要组成部分, 其主要功能是为我国华南和西南地区提供生态安全屏障, 这对于局部甚至全球区域的气候调节、植被修复、土壤保持及水源涵养均起到了至关重要的作用^[1-2]。生态系统服务是指各种类型的生态系统对人类福祉和效益的直接或间接贡献, 为人类提供了一系列的生态产品与生态功能, 而人类生存发展条件的好坏取决于生态系统服务的变化^[3-4]。然而, 随着当前社会经济的不断发展, 城市化进程逐渐加快, 城市扩张直接改变了土地利用模式, 进而导致诸多生态系统的结构发生改变^[5-6], 而人类对于自然生态系统提供的各种服务功能的依赖性逐渐增强^[7-8], 因此, 了解和评估城市化背景下的南方丘陵山地带生态系统服务价值并进一步分析其与城市化的时空关系是有效保护和提升该地区生态系统服务功能价值的前提。

量化评估生态系统服务价值有多种方法, 采用经济学的方法将其量化能够较直观的反映变化规律, Costanza等^[9]发表了全球生态系统服务与自然资本价值核算的相关研究内容, 推动了生态系统服务价值量化评估进展；谢高地等^[10]根据中国的实际情况制定了中国生态系统服务价值当量因子表, 为中国生态系统服务价值评估相关研究提供了参考。关于生态系统服务量化评估的研究从注重社会环境条件改变对生态系统服务产生的影响^[11-12]发展到运用模型并基于不同的尺度评估其服务功能价值, 例如, 在全球尺度^{[9, 13]}、区域尺度^{[4, 10]}、小流域尺度^[14-15]等尺度下进行空间变化分析。目前, 城市化与生态系统服务价值关系的研究越来越受到重视, 多数研究是将城市发展作为区域生态系统服务发生变动的驱动因素^[16], 并运用不同的统计模型探讨城市化因子对生态系统服务的影响^[17-19], 另外城市化对各种生态系统服务间权衡协同关系的变化所产生的影响^[20]以及这两个变量间的交互耦合作用^[21-22]也受到广泛关注, 这其中涉及到的生态系统包括城市^[23]、土壤^[24]、水体^[25]、农林^[26]以及文化生态系统^[27]等, 然而以往研究的地区多集聚在我国中东部城市群或经济相对发达的平原区域, 针对内陆局部山地区域的生态系统服务与城市化关系的研究还较少, 并且以往的研究尺度相对较大, 评估多基于行政单元的统计数据, 从而弱化了指标在区域内部空间上变化特征的探讨^[28-30]以及指标间的空间差异分析^[31], 因此本研究采用空间格网进行定量化评估, 突出研究指标的空间异质性特征, 便于进行空间关系的探讨。

南方丘陵山地带位于长江流域与珠江流域的分水岭及源头区, 贯穿广东、广西、江西及湖南等省市^[1], 粤北区域属于南方丘陵山地带的典型区域。地形复杂、资源丰富, 生态系统服务价值受地形的影响较大。因此, 本文以粤北为研究地区, 将该区域划分为5个地形梯度, 基于2km×2km的空间网格评估不同地形梯度上及不同年份的生态系统服务价值与城市化水平的时空关系。本研究对于促进粤北地区生态系统保护与经济建设的协调合理发展具有重要意义, 也为保障南方丘陵山地带的可持续发展提供理论依据。

1 研究区概况

粤北地区位于广东省北部, 地理位置在111°55′—114°81′E, 23°27′—25°57′N之间, 包括韶关与清远2个地级市, 其中包括2个市辖区、4个县级市、7个县以及3个少数民族自治县。全区总面积为37031 km², 占整个广东省面积的20.61%, 粤北整体地势呈现北部高南部低, 地势地貌复杂, 最高海拔为1900m, 最低为5m, 坡度范围在0—62°之间, 地形位指数在0.003—1.826之间(图 1)。该区域是典型的山地丘陵地带, 具有丰富的林业资源, 其森林覆盖面积达到了整个粤北地区的70%以上, 同时该地区也具有丰富的生物资源, 使得粤北地区逐渐形成了许多以森林为主体的动植物共生竞长的生态系统, 因此该区域也成为我国南方珍稀动植物的物种基因库。粤北全年四季如春, 气候温和, 雨量充沛, 年平均气温在19—23℃之间, 年平均降水量在2145—2648 mm之间。

2 研究方法 2.1 数据来源与处理

本研究利用1990、1995、2000、2005、2010年及2015年6期的Landsat TM/ETM遥感影像, 结合研究区域1∶25万地形图、土壤和植被专题图进行分析。根据土地利用覆盖变化(LUCC)分类体系与粤北土地利用现状, 结合ENVI软件通过对遥感影像图进行拼接、校正、裁剪与融合等相关处理, 将研究区的土地利用类型划分成耕地、林地、草地、水域、建筑用地及未利用地6类一级类型及13类二级用地类型, 对每期数据进行精度检验, 精度均达到85%以上, 满足研究对土地利用分类数据的精度要求。粤北相关社会经济数据均来源于广东省统计年鉴等。为便于研究生态系统服务价值与城市化数据的空间分布特征, 采用2km×2km网格进行空间采样, 利用Python编程语言对6期数据共计5.889万个网格样方内的数值进行计算。

2.2 地形梯度分级

地形位指数是综合海拔与坡度因子的地形因子, 其优势是可以全面综合地反映某一点的地形状况^[30]。本研究通过公式计算出粤北地区每个空间网格单元的地形位指数, 并基于ArcGIS 10.4中的自然断点工具将地形位指数计算结果划分成5个地形梯度等级(图 1)：0.003—0.368、0.368—0.639、0.639—0.890、0.890—1.154、1.154—1.826。后将地形梯度与ESV及UL的评估结果在GIS中进行空间叠置, 便于分析在地形梯度因子的影响下ESV与UL的动态规律及两者的关系变化。

地形位指数计算公式：

(1)

式中, T表示地形位指数, E与S分别代表某点的海拔(m)与坡度(°), E₀与S₀分别代表整个研究区的平均海拔(m)与平均坡度(°)。

2.3 土地利用变化图谱研究

图谱分析可以对不同时期的土地利用单元进行地图代数叠加实现图谱的融合^[32-33], 计算公式为：

(2)

式中, L代表在研究阶段内新的土地利用图谱单元属性图；A代表前一时期的土地利用图谱单元属性值；B代表后一时期的土地利用图谱单元属性值。

本研究将图谱的变化特征分为5种变化模式：前期变化型(土地利用只在1990—2000年发生变化)、中期变化型(土地利用只在2000—2010年发生变化)、后期变化型(土地利用只在2010—2015年发生变化)、反复变化型(1990与2015年土地类型一致且2010与2015年的土地类型一致；1990与2010年土地类型一致且2000与2015年土地类型一致)、连续变化型(1990、2000、2010及2015年中至少出现3种不同土地类型), 计算公式为：

(3)

式中, L代表在研究阶段新的土地利用图谱单元属性图；A′、B′、C′、D′分别代表1990年、2000年、2010年、2015年粤北土地利用图谱单元属性值^[34]。

2.4 城市化水平量化评估

本文借鉴已有的相关研究, 从人口城市化水平、经济城市化水平以及土地城市化水平三个层次衡量粤北地区的综合城市化水平^[28], 分别使用人口密度(Population Density, PD)、国内生产总值密度(Gross Domestic Product Density, GDPD)及建筑用地面积占比(Urban Land Proportion, ULP)三个指标来量化, 并依据自然断点法将评估结果分为5个等级, 即低水平、较低水平、中等水平、较高水平及高水平, 具体公式如下：

(4)

式中, UL是指城市化水平指数；PD′为标准化后的人口密度；GDPD′为标准化后的GDP密度；ULP′为标准化后的建筑用地面积占比。

2.5 生态系统服务价值评估

本研究采用价值当量法^[10], 以粤北地区2010年单位面积粮食产量为基准, 计算出粤北地区的生态系统服务当量修正系数为0.987。利用修正后的价值当量因子(表 1)计算每个空间网格的生态系统服务价值, 并使用自然断点法划分为5个等级：低值、较低值、中等值、较高值、高值。

生态系统服务价值评估公式：

(5)

式中, ESV表示所研究生态系统总服务价值；A_i表示第i种土地利用类型面积(hm²)；VC_i 表示第i种土地利用类型生态系统服务价值系数。

2.6 空间相关性研究

本文使用空间分析软件Geo Da中的双变量Moran′s I指数来探索粤北地区生态系统服务价值(ESV)与城市化水平(UL)的空间相关性。使用全局双变量研究整体研究区生态系统服务价值与城市化水平之间的空间关系, 而使用局部双变量研究整个地区内部不同空间单元的生态系统服务价值与城市化的空间关系^[35-36]。针对粤北地区生态系统服务价值与城市化水平的空间关系, 绘制了空间聚集图, 其聚集类型包括：不显著、高-高、高-低、低-低及低-高5种类型。

3 结果与分析 3.1 基于地形梯度的土地利用变化图谱分析

1990—2015年粤北地区土地利用变化图谱见图 2：(1)研究时段内各个土地利用类型均发生转移, 其中耕地与林地图谱单元转出率最高且在地形梯度1上最活跃, 其转移的图谱单元数量分别占发生转移的图谱单元总数量的41.3%与35.4%, 其次是草地转出率为20.4%且在梯度2、3上较活跃, 而水域与建筑用地同样在梯度1上最活跃, 但总转出率仅为2.9%。(2)研究时段, 建筑用地的转入率高达45.1%, 是占比最高的土地利用转入类型, 主要集中在地形梯度1的区域, 其次为林地, 转入率为28.3%且集中于梯度1、2区域, 草地主要由较高梯度的林地转移而来, 其转入率为17.5%, 而水域与耕地总转入率仅为9.1%且活跃于地形梯度1区域。(3)从图谱的5种转移模式来看, 其图谱单元数量占比由高到低依次为：中期变化(41.2%)、后期变化(33.2%)、前期变化(22.8%)、反复变化(2.1%)、连续变化(0.7%)。中后期的土地利用转移较为频繁, 对比其转移情况发现中后期主要的转移类型为耕地和林地转移为建筑用地, 其次是部分耕地转移为林地, 而这两种转移类型主要分布于1、2地形梯度区域；在前期主要是1、2梯度上的部分林地与水域被开发为耕地与草地；反复变化类型主要是2、3梯度上的耕、林、草之间重复转换, 且最终转移为林地或草地；连续变化型由于面积占比太小, 其变化具有偶然性。

3.2 基于地形梯度的城市化时空动态分析

1990—2015年粤北地区的城市化指数依次为0.018、0.018、0.019、0.038、0.037、0.043, 可以看出城市化水平整体呈上升的状态。5个地形梯度上的城市化水平见图 3, 1990—2015年期间尽管城市化程度不断加剧, 但粤北地区处于低城市化水平的区域面积占比较大, 且随地形梯度的升高, 低城市化水平区域面积随之增加；而较高及高城市化水平区域主要聚集于城市密集分布的地形梯度1上, 从1990年至2015年, 地形梯度1区域上的较高及高城市化水平区域面积占整体较高及高城市化水平区域面积比例分别为0.707、0.698、0.717、0.771、0.754、0.712, 可看出该区域面积在2000年后显著增加, 可见地形梯度1区域的城市化程度最高, 表明该区域适合人类开发与居住, 地形梯度对城市化的空间分布具有显著的影响。

3.3 基于地形梯度的生态系统服务价值时空动态分析

1990—2015年粤北地区生态系统服务价值呈轻微减少的趋势。各地形梯度上的生态系统服务价值随时间变化均呈现缓慢下降的趋势, 在研究时段内随着地形梯度的逐渐上升, 对应各梯度上的生态系统服务价值随时间变化的下降率逐渐减小。这主要是由于低梯度区域海拔低、坡度小, 城市、工业区与人口集中分布, 城市的快速发展以及较频繁的人类活动使得该区域自然资源的开发与利用程度较高, 进而导致该区域的生态系统服务价值下降率最高；而随着地形梯度的上升, 由于气候、坡度等环境条件的限制导致城市化程度降低, 使得原自然生态环境得以保留, 对于局部外来干扰具有较大的抵抗力与恢复力, 因此各梯度上的生态系统服务价值下降率也随之减少。

1990—2015年粤北地区生态系统服务价值空间分布情况见图 4。在研究期间, 粤北地区大面积区域的生态系统服务处于较高值与中等值, 其网格数量分别占总体网格数量的45%与35%左右, 将这些区域与地形梯度叠加发现其主要处于地形梯度3、4区域, 由于这里的地形与气候等环境条件可能限制了人类活动, 使得该区域丰富的自然资源与生态环境得以保护与维持, 因此粤北整体的生态系统服务处于良好的水平；然而, 低值区域与较低值区域集中于地形梯度1、2地区, 其网格数量在研究期间逐渐增多并且在2010年后增长率显著上升, 这说明低梯度区域的城市化进程加快, 可能导致了较大的生态系统服务价值损失, 因此该区域的生态系统服务处于较低水平；高值区域的面积占比最小, 而这类区域除了分布在梯度最高的区域外, 还分布在低梯度区域林地与水域集中的地区。

3.4 基于地形梯度的城市化与生态系统服务价值空间相关性分析

粤北地区生态系统服务价值与城市化水平全局空间相关性结果见图 5, 可以看出ESV与UL呈现较显著的负相关关系, 且两者的空间负相关效应趋于增强。数据的显著性检验结果见表 2, 由此可知ESV与UL在空间上的聚集效应也在增强。

粤北地区生态系统服务价值与城市化水平的局部空间关系见图 6, 本研究将其与粤北行政区划图进行空间叠加, 可知低-高型主要聚集在南部的清远市市辖区、中北部的韶关市市辖区及东北部的南雄市, 且这三个行政单元基本都位于地形梯度1区域, 其中两个市辖区是整个粤北城市与人口的核心聚集区域, 地势较低, 具有较高的城市化水平, 因此该区域提供较低的生态系统服务价值；而南雄市农业经济与旅游业发达, 综合开发程度较高, 导致该区域生态系统服务价值较低。高-低型的分布面积较广, 其主要聚集区为粤北北部的三个少数民族自治县及东部的始兴县, 该地区整体海拔较高, 森林资源丰富, 是少数民族村落聚集地, 城镇化率较低, 因此提供有较高的生态系统服务价值。另外, 高-高型与低-低型的分布面积较少, 高-高型聚集于粤北南部的清新县, 该区域首先具有丰富的森林资源与水力资源, 因而提供较高的生态服务价值；低-低型在研究区呈零散分布, 该区域多为农田与建筑用地, 资源类型单一, 生态系统服务价值较低。

4 讨论与结论 4.1 讨论

粤北地区属于丘陵山地带的典型区, 地形复杂, 因此本研究采用综合地形位指数反映每个空间网格评价单元的地形特征, 并基于空间单元网格评估了粤北地区的生态系统服务价值与城市化水平, 并对两者进行空间关系及聚类效应的研究。研究结果表明, 低梯度上的耕地主要是在前期(1990—2000年)进行扩张, 而到中后期(2000—2015年), 低梯度1, 2区域上的大面积耕地与林地则转化为建筑用地；在中后期阶段也有部分耕地转化为林地或草地。另外在研究期间, 粤北城市化水平随地形梯度的升高而降低, 且评估结果显示2000年之后的城市化水平显著上升, 这与图谱的研究结果相符, 而生态系统服务价值则是随梯度的升高而增加, 且粤北地区的生态系统服务价值整体处于较好水平。另外粤北的生态系统服务与城市化具有较强的负相关关系, 且在局部区域产生了显著的空间聚集效应。

粤北地区生态系统服务价值在空间上呈现出显著的分异规律, 其主要原因如下：自然环境条件的限制效应^[37]与城市化带来的胁迫效应^[38]。从自然环境条件方面来说, 由于粤北具有大面积的丘陵山地, 复杂的地形条件使得城市化与生态系统服务价值呈现显著的分异特征, 这主要是由于地形位指数高的区域海拔高, 坡度大, 再加上不适宜的气候条件等多种环境因素限制了人类的各种生产活动与开发建设^[37], 因此较高地形梯度上的自然生态系统受到了低程度的干扰与破坏, 这可能是山地的ESV一般要高于平原地区的原因^[39-40]。从城市化方面来说, 粤北的自然资源丰富, 导致人类对各种自然资源进行开发和利用, 并破坏局部自然环境进行城市建设与耕地扩张, 在此过程中也造成了一系列的环境污染问题, 比如水污染与大气污染等, 因此城市快速扩张与发展对自然生态环境产生了明显的胁迫效应^[38], 这类现象多聚集于地形梯度1、2区域, 而人口越密集, 胁迫效应越显著, 如清远市辖区与韶关市辖区等。在研究期间, 粤北经济也在快速发展, 然而粤北整体的ESV变化较稳定, 只有轻微下降, 这说明研究区原生生态环境本底条件较好, 城市化产生的胁迫作用只是对局部地区的ESV产生了显著的负面效应, 并没有对整个研究区的ESV产生太大的影响。

基于研究结论, 对于粤北地区的环境管理工作提出如下建议：(1)针对低梯度人口集中的区域(低ESV-高UL型与低ESV-低UL型), 建议稳定城市化建设发展的速度, 控制建筑用地面积的增长, 并且在城市建设的同时注重生态环境的保护, 加强绿色基础设施建设与生态恢复建设, 从而利于ESV水平的提升。(2)针对高梯度区域(高ESV-低UL型), 由于该区域受到的干扰程度低, 生态环境保持良好的状态, 因此造成的破坏容易得到恢复, 所以这部分区域应该着重进行保护, 尤其是林地与草地自然生态系统, 避免未来的过度开发与利用。(3)清新县(高ESV-高UL型)的经济发展模式可以为粤北地区未来的建设发展提供方向, 充分利用好粤北地区自然地理环境优势发展生态旅游等绿色经济产业, 将自然因素与社会经济因素协调合理化, 将绿色产业打造成该地区的优势产业。

本研究旨在探讨ESV与城市化的空间格局关系, 因此通过空间网格进行指标的量化评估, 但评估时的单位网格面积与土地利用数据的精度有密切关系, 本文1km²的土地分类栅格数据在目前的研究中已经相对粗略, 因此在未来的工作中要做到更细致的评估则需要提高基础数据的分辨精度^[41]。另外, 本文只是探讨了粤北ESV变化与城市化的关系, 在广阔的空间中ESV的变化不仅只受到城市化影响, 而且城市化也只是属于人类对自然干扰的一种方式, 除此之外也不可避免众多自然因子的参与, 因此ESV的时空变化机制及其驱动因素仍需要进行更多的研究来完善。本研究可为粤北地区相关环境决策部门和工作人员提供理论依据与相关参考。

4.2 结论

本文以粤北地区为例, 从地形梯度的角度出发, 首先利用图谱将1990—2015年粤北地区的土地利用转移动态可视化, 并基于空间分析评价了整个研究区的城市化水平及生态系统服务价值的时空变化特征, 并对两者的时空关系进行了研究, 明确了发生聚类效应的分布区, 为当地相关决策制定提供了较准确的指导。根据研究结果得出如下结论：

(1) 粤北地区土地利用转移动态在地形梯度1最为活跃, 而在研究期间主要的土地利用转移类型是耕地、林地转为建筑用地, 且中后期建筑用地的转入率明显上升。

(2) 粤北地区在研究期间城市化水平升高, 且城市化进程在2000年之后显著加快, 地形梯度1的城市化程度最高, 而随着地形梯度上升, 城市化水平逐渐下降。

(3) 粤北地区整体的ESV在研究期间略微下降, 呈现出较稳定的发展态势, 而ESV在各地形梯度上随时间变化同样呈现出下降的态势, 且在地形梯度1区域下降率最高。ESV的空间分布因地形因素呈现出显著的分异特征, 较高ESV分布于高梯度、林地草地等自然植被分布的山地区域, 而低ESV主要分布于低梯度、耕地与建筑用地等类型分布的平原区域。

(4) 空间相关性研究结果表明生态系统服务价值与城市化呈现较强的负相关性, 且在研究期间负相关性明显增强, 同时两者存在非常显著的空间聚集效应, 其中低-高型集中分布于粤北两个市辖区, 高-低型集中分布于人口相对较少, 城市化水平较低的少数民族自治县。叠加地形梯度与土地利用类型分布图可知低-高型分布在耕地与建筑用地聚集的地形梯度1、2区域, 该地区综合开发程度较高, 已严重影响到自然生态环境；高-低型地区海拔较高, 林地资源丰富, 人口相对分散, 应当继续保护, 禁止不合理的开发工程, 结合当地实际情况进行科学规划, 实现经济建设与自然环境平衡发展。