生态系统服务是指通过生态系统的结构、过程和功能直接或间接得到的生命支持产品和服务^[1], 对人类健康与生存以及区域和全球生态安全至关重要^[2-4], 是生态学和地理学研究的前沿和热点^[5]。土地利用/覆被变化(Land Use/Cover Change, LUCC)作为人类活动行为与地球陆表自然生态系统最直接的表征形式^[6-7], 对生态系统服务价值(Ecosystem Service Values, ESV)起决定性作用, 其评估结果是高效、合理配置竞争性环境资源的基础^[8], 也是制定生态补偿^[9]、可持续发展规划政策^[10-12]以及生态系统服务权衡^[13-15]的重要前提。

目前, 国内外对于ESV定量评估的研究日趋成熟, Costanza等^[1]于1997年率先在全球尺度上对生态系统服务功能进行划分与价值评估, 成为后来研究ESV的重要基础。谢高地等^[16-17]在Costanza评价模型基础上改进并制定了“中国陆地生态系统单位面积生态服务价值当量表”, 被广泛应用在我国的森林^[18]、草原^{[8, 16]}、湖泊^[19-20]、城市圈^{[11, 21-22]}等区域ESV评估中。白杨等^[18]利用市场价值法、影子工程法和生产成本法等定量评价了海河流域森林生态系统服务功能的经济价值；李晋昌等^[8]以修正后的价值系数为基础, 对若尔盖高原ESV进行了评估；姚小薇等^[21]采用价值当量法核算了武汉城市圈ESV, 并运用空间计量方法探讨不同城镇化水平对ESV空间分异的影响；Zhang等^[23]采用双变量Moran等指数分析了武汉市ESV和城市化之间的空间相关性, 结果显示两者存在空间负相关关系和空间溢出效应。相关研究表明, 社会经济活动引起的LUCC是导致生态系统服务供给变化的主要因素, 不同土地利用格局会产生相应的生态过程, 从而对生态系统服务造成影响^{[5, 24]}。近年来随着研究的逐步深入以及3S技术的发展, 越来越多的学者也开始注重土地利用对生态服务价值时空分布格局及其影响因素的探讨^{[5, 11, 20]}, 但依然缺乏基于空间统计与分析来探索土地利用和ESV在空间上的聚集规律、关联模式等分布特征的定量研究^{[20, 25]}。

空间统计分析可定量揭示地理变量的空间关联性问题, 是研究区域土地利用变化与ESV空间格局特征的重要工具, 也是地理学领域研究的重要方向^{[19, 26-28]}。海南岛作为一个独立的地理单元, 生态系统丰富多样, 但随着国际旅游岛建设的不断推进, 土地利用变化剧烈, 与生态系统之间的矛盾也逐步凸显。位于海南岛东北部的海口、文昌2市是省会经济圈建设和国家航天发展战略所在地, 承担着区域经济和社会可持续发展重任。本文以海南岛东北部为研究对象, 运用空间自相关理论分析土地利用与ESV空间格局的分布特征, 旨在解决：(1)采用空间统计方法探索研究区土地利用空间自相关关系；(2)分析海南岛东北部地区ESV的空间分布变化；(3)通过双变量Moran′s Ⅰ探讨土地利用与ESV之间的空间相关性, 分析ESV对土地利用的空间依赖性。研究结果可为区域土地利用规划管理、生态系统服务价值的保育与恢复提供实际指导。

1 研究区概况

研究区位于海南岛东北部(110°07′—111°03′E, 19°20′—20°10′N), 包含海口、文昌2市全境, 是“海澄文”一体化经济圈的重要组成部分, 下辖41个镇、农场, 面积4748.27 km², 约占全省土地总面积的14%(图 1)。区域北濒琼州海峡, 东面南海, 西南接澄迈、定安和琼海3市县；地形西南高东北低, 以锅盖岭、马鞍岭为隆起高点, 大部以滨海平原地貌为主, 海拔在0—388 m；年均气温24℃, 年降水量1900 mm, 平均相对湿度84%, 属于热带海洋性季风气候。区域内有海南岛最长的河流南渡江从海口市中部穿过而入海, 水资源丰富, 且有海南东寨港、清澜港两大红树林自然保护区。2016年, 区域常住人口278万, 地区生产总值1332亿元, 分别占海南省的31%和36%, 是海南省政治、经济、交通、科技、文化中心和我国文昌航天发射场所在地。

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

研究采用的土地利用/覆被数据来源于2016年Landsat 8 OLI遥感影像(中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)下载)解译数据, 并结合海南省森林资源二类调查矢量数据进行修正, 数据精度达到95%以上。依据全国土地利用调查分类体系和研究区土地利用整体特征, 将土地利用类型分为7类, 分别为林地、耕地、草地、建设用地、湿地、水体及未利用地。借助ArcGIS 10.3 Create Fishnet工具将研究区划分为1 km×1 km正方形网格单元, 共计4967个, 计算每个网格单元中各土地利用类型的面积。

以研究区单元格各土地利用类型的面积百分比作为空间探索变量, 采用经验贝叶斯(Empirical Bayes)平滑方法^[29-30], 对面积百分比数据进行标准化处理, 并对标准化后的数据进行空间自相关分析。

2.2 研究方法 2.2.1 土地利用程度指数

土地利用程度在一定程度上反映着人为活动的干扰程度^[25]。根据庄大方和刘纪远^[31]提出的数量化土地利用程度分析方法, 将土地利用程度按照土地自然综合体在社会因素影响下自然平衡保持状态分为4级, 并分级赋予指数, 计算公式为：

(1)

式中, L为研究区域土地利用程度综合指数；n为土地利用类型的数量, 这里取7；A_i为第i类土地利用类型的面积；A_T为研究区域总面积；P_i为不同类型的土地利用程度参数, 其中：未利用地赋值1, 林地、草地、湿地和水体赋值2, 耕地赋值3, 建设用地赋值4^{[11, 31]}。

2.2.2 生态系统服务价值评估

生态系统服务价值包含供给服务、调节服务、支持服务和文化服务4大类型所提供的价值^[2]。目前, ESV评估方法大致可分为“直接市场价值法”和“条件价值法(假想市场法)”两种, 但国内相关研究仍以直接市场价值法为主^{[20, 32]}, 特别适用于区域及全球大尺度ESV评估^[33-35]；而后者在具体理论和方法上尚未形成突破^[32]。因此, 本研究依据谢高地等^[16-17]的价值当量因子换算方法, 并结合研究区实际情况进行修正^[26]。根据《海南统计年鉴》, 以2000—2015年海口、文昌两市的平均实际粮食产量和粮食平均价格, 计算出研究区ESV的当量为1307.34元, 从而确定了不同土地利用类型单位面积的生态系统服务价值系数(表 1), 其中建设用地的价值系数为0。

对每个网格分别计算生态系统服务功能价值和价值强度(按面积消除行政区边缘不完整网格的价值高低差异), 并进行空间特征分析^{[1, 36]}。

每个网格的生态系统服务功能价值：

(2)

每个网格的生态系统服务功能价值强度：

(3)

研究区生态系统服务功能总价值：

(4)

式中, A_ij为第j个网格i种土地利用类型的分布面积(hm²)；C_vi为第i种地类生态价值系数(元hm^-2 a^-1)；S为每个网格的面积；n为土地利用类型；m为网格个数。

2.2.3 空间自相关分析

空间自相关分析是用于衡量空间变量的分布是否具有集聚性, 包含全局空间自相关和局部空间自相关两方面^[21]。为了揭示多个变量之间的空间相关性, Anselin^[37]在此基础上提出双变量空间自相关, 揭示空间单元属性值与邻近空间上其他属性值的相关性^[22]。本文空间自相关分析采用GeoDa 1.6.7软件完成。

全局空间自相关能够反应整个研究区域内各个地域单元与邻近地域单元之间的相似性。全局Moran′s Ⅰ是被广泛应用的全局自相关统计量, 计算公式为：

(5)

局部空间自相关指标(Local Indicators of Spatial Association, LISA)常采用局部Moran′s Ⅰ统计量进行度量, 用以准确地把握局部空间要素的聚集性和分异特征, 在z检验的基础上(P < 0.05)绘制LISA分布图^[30], 计算公式为：

(6)

式中, n是空间单元数量, x_i和x_j分别表示单元i和单元j的观测值, (x_i- x)是第i个空间单元上的观测值与平均值的偏差, w_ij是基于空间k邻接关系建立的空间权重矩阵, 方差

2.2.4 半变异函数分析

半变异函数分析是反映区域化变量的空间相关性, 进行空间变量结构分析和最优模拟的重要工具, 计算公式及详细论述参看文献^[38]。变异函数理论模型的最优选择用决定系数R²来决定, 并综合考虑RSS、块金值C₀和变程A₀。空间变量经log变换后, 其分布符合正态分布规律, 满足半变异函数分析的前提条件, 借助地统计学软件GS+ 9.0应用最小二乘法以正北方向呈0°、45°、90°和135°4个角度为典型方向, 分别对实测半方差进行球面模型、指数模型、线性模型和高斯模型4种不同模型的变异函数理论模型拟合, 估计不同距离的半变异值, 进而采用ArcGIS对变量进行普通克里格法(Ordinary Kriging)插值分析^{[36, 39]}。

3 结果与分析 3.1 土地利用格局 3.1.1 全局空间自相关

从表 2检验结果显示, 海南岛东北部各土地利用类型全局Moran′s Ⅰ值均大于0, P值均小于0.001, 说明研究区内各土地利用类型整体上呈显著的正向空间自相关关系, 具有非常明显的聚集性, 在空间分布上并非完全随机(表 2)。其中建设用地和林地表现最为突出, 空间聚集性最强, 湿地和未利用地的空间自相关性相对较弱。

3.1.2 局部空间自相关

土地利用LISA分布图(P < 0.05)直观反映了7种土地利用类型在空间聚集与分异的位置分布特征(图 2)。其中, 林地主要聚集在研究区西南部锅盖岭和马鞍岭等较高海拔地区(HH聚集, 高-高聚集), 植被覆盖程度高, 面积百分比均在50%以上；而林地LL聚集(低-低聚集)主要分布在海口市建成区, 由于城市建设扩张造成林地占比很低。耕地HH聚集主要出现在文昌市北部的铺前、罗豆、锦山3镇, 在城市建成区、海岸带等不适宜耕作或农业发展的区域表现为LL聚集。草地HH型“组团”出现在城乡结合部或库塘周边, 这些区域很适合热带杂草生长或传播, 而在西南部则因乔木林地的高度覆盖造成纯草地分布很少。建设用地HH聚集明显分布于海口、文昌城区以及重要乡镇。湿地和水体则“组团”聚集分布在东寨港、清澜港、海岸带、南渡江及库塘等区域。而未利用地主要是城市待开发建设区域所形成的人工堆掘地或裸露地表, HH聚集集中分布在海口城区周边、文昌北部木兰湾及东部月亮湾等区域。

3.1.3 土地利用程度空间自相关

以土地利用程度综合指数为观测变量, 从图 3可以看出, Moran散点主要分布在第一象限(HH)和第三象限(LL), 第二象限(LH)和第四象限(HL)散点分布相对较少, 全局Moran′s Ⅰ指数为0.5687, 这说明海南岛东北部土地利用程度指数具有很强的空间正相关性。

从LISA分布图来看(图 4), 土地利用程度HH聚集主要分布在海口、文昌城区及交通便利的重要乡镇, 而在文昌市北部的铺前、罗豆、锦山三地因有大面积集约经营的耕地, 因此土地利用程度也较高。LL聚集区主要出现在未开发沿海地带以及东寨港、清澜港和库塘区域, 可能与自然保护区限制、利用困难或交通不便等因素有关, 造成土地利用程度较低。而HL聚集和LH聚集类型则在研究区内呈零星分布。

3.2 生态系统服务价值评估 3.2.1 各土地利用类型生态系统服务价值分析

从表 3可以看出, 2016年海南岛东北部生态系统服务总价值为132.09亿元, 占2016年研究区GDP的9.92%。其中林地价值最高, 占67.09%, 是研究区ESV的贡献主体；除建设用地外, 未利用地价值最低, 仅占0.13%。从ESV强度来看, 湿地和水体最高。每个网格单元(1 km²)的平均ESV为265.44万元, 最小值为0.02万元, 最大值达716.03万元, 各网格单元价值差异较大。

3.2.2 生态系统服务价值空间自相关分析

为避免研究区边缘网格单元面积大小不同造成ESV相差较大, 采用ESV强度进行空间自相关分析。结果显示, 海南岛东北部ESV强度的全局空间自相关Moran′s Ⅰ指数为0.5041, 同样表现出很强的空间正向自相关。从LISA分布图看(图 5), ESV强度较高的区域主要分布于东北部海岸带、东寨港和清澜港红树林、以及大型水库等, 在这些湿地和水体的主要分布地带, 表现为显著的HH聚集, 也反映出湿地和水体的单位ESV很高。LL聚集“组团”出现在海口城区和南渡江近入海口两侧, 人口分布密集, 建设用地面积比例很高, 导致ESV很弱。其次为文昌北部的耕地集中分布区和重要乡镇分布点, 同样因城镇用地和耕地集中分布也造成ESV较低。而在林地集中分布的研究区西南部区域也出现有零星的HH聚集, 但因差异不大, 表现并不十分显著。

3.2.3 半变异函数分析

研究区ESV强度 V_ES经变异函数理论模型拟合, 分析认为采用指数模型的拟合效果最好, 其块金值C₀=0.0001, 基台值C₀+C=0.1222, 决定系数R²=0.724, 残差RSS=7.809×10^-4, 变程A₀为2.6 km, 表明V_ES在1 km的网格间距内具有很强的空间自相关性(图 6)。因此, 运用ArcGIS空间分析模块中的普通克里格法进行插值生成 _ES的连续空间分布图(图 7)。

图 7所示, 海南岛东北部ESV强度呈现以海口城区、文昌北部农耕区为主要低值核心, 东寨港、清澜港及海岸带为高值分布地带的总体分布格局, 模型拟合的结果与ESV强度LISA分布格局(图 5)高度一致。高值区域明显聚集在湿地、海岸带等重要生态保护区, 生态系统服务功能很强；而低值区则容易出现在交通干线连接起来的人类活动密集区, 这些区域往往是城市建成区或耕地分布的重点区域, 在带来较好经济与社会效益的同时, 也容易导致ESV损失。由此可见, ESV强度在空间分布上与土地利用格局、区位因素和经济发展水平有很大关联。

3.3 双变量空间自相关关系

为进一步探索海南岛东北部土地利用程度与ESV之间的空间关系, 采用双变量空间自相关检验结果显示, Moran′s Ⅰ指数为-0.3937, 表明土地利用程度与ESV之间存在显著的空间负相关关系, 即随着土地利用程度的增强, ESV总体上呈现出下降趋势。从土地利用程度与ESV双变量LISA空间分布图来看(图 8), HL型主要分布在海口城区、文昌北部以及重要乡镇等, LH型分布在东寨港和清澜港红树林自然保护区、海岸带、以及大型水库等区域。从双变量LISA显著性水平上看(图 9), 在研究区城镇建设区、东寨港和清澜港红树林分布区、库塘、海岸带等土地利用强度和ESV呈HL型或LH型分布的区域表现为极显著相关(P < 0.01), 而在其周边部分区域表现为显著相关(P < 0.05)。

4 讨论 4.1 土地利用空间自相关格局

LUCC作为全球及区域气候与生态环境变化的主要根源, 受到自然、社会、经济等诸多因素影响^{[30, 40]}。Verburg等^[41]研究认为, 自然生态环境因子在土地利用变化的空间分布上起主导作用, 而社会经济因子在决定土地利用变化的数量特征方面为主。本文研究结果表明, 在研究区西南部的自然生态环境更适合林木生长的山地地区, 林地呈明显聚集状态, 水体则聚集在地势平缓东北部区域。受区域发展和交通等因素的影响, 在文昌北部的铺前、罗豆、锦山一带, 耕地呈“组团”聚集分布；而城镇建设用地则明显聚集在交通干线便利的人口密集区, 反映了区域土地利用的空间自相关格局和分布特征。事实上, LUCC是“人类-环境”综合作用的动态复杂系统^[42], 驱动因子和驱动机制多样、复杂^[43], 也造成土地利用呈明显的区域差异性或聚集性, 且具有很强的空间相关性。城镇用地扩张、人口聚集使得土地利用程度的高值在海口城区及周边地区的空间自相关性不断加强, 耕地的“组团”分布也形成文昌北部区域的土地利用程度较高；同时, 得益于近年来海南省海岸带生态修复和湿地保护政策, 使得东部海岸带及重要湿地周边的土地利用程度呈显著LL聚集。

4.2 生态系统服务价值空间分布特征

本文采用基于修正系数价值当量因子换算的直接市场方法进行ESV评估, 研究结果表明, 2016年海南岛东北部生态系统服务功能总价值为132.09亿元, 其中林地构成ESV的主体。与其他研究相比, 喻露露等^[26]对2012年海口市海岸带ESV的评估结果为39.24亿元, 从面积换算来看两者的评估结果比较接近；高玲等^[44]采用市场价格法等方法评估2008年海口市ESV达106.02亿元, 评估结果有较大偏差。这是由于不同评估方法在评估模型、参数设定等方面的不同, 造成结果差异较大^{[8, 25]}, 但这并不影响区域内ESV空间上的比较, 研究结论依旧具有可靠性。在空间分布上, ESV强度表现出明显的自相关分布特征, 高值区(HH聚集)位于湿地和水体的主要分布地带, 这可能跟生态系统内丰富生物多样性和环境组分之间物质、能量、信息相互作用有关；低值区(LL聚集)位于海口城区、文昌北部耕地分布区以及重要乡镇, 该区域受城市土地利用开发等人为活动影响, 打乱了原本系统内平衡的生态过程, 在带来较好经济与社会效益的同时, 也容易导致ESV减弱甚至退化。

4.3 生态系统服务价值对土地利用变化的响应

地理要素的空间交互作用中普遍存在空间溢出效应, 它是指由基于位置的邻近性所产生的空间外部性, 即一个单位从其邻居带来收益或成本^{[23, 45]}。大量研究表明, 土地利用程度对ESV呈负相关^{[20, 25]}, 且存在明显的空间自相关性和空间溢出效应。总体而言, 土地利用程度对ESV有负面的空间溢出效应, 也就是说, 一个地区的土地利用程度的提高可能会导致周边地区的ESV的退化。土地利用程度较高的如建设用地、耕地等是ESV的低值分布区, 而土地利用程度较低的如湿地、水域则为ESV的高值分布区。因此, 在高ESV附近的区域, 应避免耕地活动和城市建设开发产生的溢出效应对自然生态系统造成的负面影响。相反, 城市自然边缘应该致力于生态友好的土地用途(例如郊野公园、绿地系统)^[46]。根据陈利顶等^[47]提出的“源-汇”景观理论, 本文研究认为, 控制建设用地、耕地等“汇”(负服务)景观类型, 保护与恢复湿地、水域、林地等“源”(正服务)景观类型, 预防“汇”景观类型向“源”景观类型扩展或渗透, 搭建生态系统内部良性循环, 积极引导土地利用向ESV保值和增值方向发展, 是优化与提高区域ESV的有效途径。

4.4 LUCC与ESV空间关系对区域景观规划的启示

LUCC与ESV是一个相辅相成的综合系统, 开发建设还是生态保护？这是区域生态规划者和管理者必须面对的问题, 如何在不破坏生态可持续性的前提下分配建设活动, 应当是生态系统服务权衡的重要任务^{[23, 48-49]}。

根据土地利用程度与ESV双变量空间关系的四种聚类模式认为, 研究区内具有低土地利用程度和高ESV的区域(LH型)是生态优势区域, 可以被认为是区域性的“生态核心区”, 比如东寨港和清澜港红树林自然保护区、海岸带以及大型水库等, 应该采取最严格的保护措施, 禁止或限制城市扩张^{[23, 26, 50]}。在高土地利用程度和低ESV的区域(HL型), 比如海口城区和重要乡镇, 应当营造绿色空间或减少建筑强度来缓解人类对生态系统的压力, 引导土地利用向高ESV发展。在高土地利用程度和高ESV地区(HH型), 应鼓励生态友好的土地利用类型(例如城市森林、草地等), 而限制高污染土地利用类型。而在低土地利用程度和低ESV地区(LL型), 则应考虑更多的生态修复工程, 也可考虑在这些地区进行一些建设项目, 增强土地利用效率。此外, 还应规划和实施有效的综合生态修复工程, 根据生态系统退化程度的优先次序有计划地实施。

5 结论与展望

(1) 海南岛东北部各土地利用类型的全局Moran′s Ⅰ值均大于0.4, P值也均小于0.001, 具有显著的空间自相关性, 其中建设用地和林地的空间聚集性最强；受人为干扰和环境作用的综合影响, 土地利用程度有显著差异与空间聚集特征, 不同土地利用类型所表现的空间聚集或异常的区域及范围也明显不同。

(2) 研究区ESV有显著的空间正相关性, 空间集聚程度较高；高值区聚集于海岸带、东寨港和清澜港红树林、以及大型水库等, 低值区集中在海口城区、文昌北部农耕区域。

(3) 研究区土地利用程度与ESV呈空间负相关关系, 即土地利用程度越高, ESV越低, 且具有空间溢出效应。根据空间聚类模式提出区域保护措施, 为城市景观规划与修复提供实际指导。

本文运用空间自相关理论探讨了土地利用与ESV空间格局的分布特征与自相关关系, 根据研究结论建议继续实施海岸带生态修复与保护工程、红树林保护管理等政策, 优化建成区土地利用结构, 提升生态核心的服务价值, 对实现区域可持续发展和维护生态安全具有重要的实践意义。本研究也存在不足之处。首先, 生态学和地理学的特征与现象往往伴随着复杂的尺度效应, 土地利用与ESV的空间自相关也表现出明显的尺度依赖特征, 如何更科学的确定空间自相关的尺度还有待于深入研究。其次, 生态系统服务具有内部关联性, 但社会经济因子与生态系统之间联系不应被忽视, 尤其是不同土地利用变化驱动情景下的生态系统过程与服务的关系, 如何揭示其背后的生态过程和响应机理应当是今后重点研究的方向之一。

致谢: 感谢海南大学宋希强、王华锋和清华大学杨军老师对本文写作的帮助。