生态系统服务价值(ecosystem services value，ESV)是生态系统服务功能的重要量化评估方式，能够直观反映生态系统服务功能效益的变化^[1-2]。随着我国城市人口数量不断增长以及工业化、城市化进程的加快，人类活动驱动下的土地利用变化正对生态系统平衡造成巨大压力，导致生态环境恶化及生态系统服务功能退化^[3]。以土地作为载体的生态系统结构在人类农业生产和各种经营活动下不断发生变化，致使ESV变化显著^[4]。研究ESV时空变化趋势，可以动态反映生态系统结构与功能的变化规律。

近年来，国内外众多学者对ESV定量评估的研究进行大量探索，并取得了丰硕成果。COSTANZA等^[5]在1997年率先提出全球ESV估算原理及方法，DAILY等^[6]、HIMES-CORNELL等^[7]、HOLMES等^[8]、TOLESSA等^[9]进一步对ESV的分类、价值化评估、核算方法以及土地利用变化对ESV的影响等多方面开展了深入研究。国内方面，谢高地等^[10]及欧阳志云等^[11]依据我国实际情况建立ESV当量表，同时确定经济价值核算方法并加以改进。这些研究引起国内学者的广泛借鉴，有效推动了ESV在国内的快速传播与发展。随着“3S”技术的发展，不少研究人员结合遥感信息数据与生态学模型，对城市区域的ESV评估等开展相关研究^[12-15]。已有研究对ESV的预测较少，所以通过空间统计与分析来探索城市扩张对未来土地利用空间格局变化的影响，对ESV在空间上的聚集特征、关联模式等进行预测模拟研究十分必要。

乌鲁木齐市位于天山山脉中段北麓，是西北典型的干旱区绿洲城市，随着城市人口的不断增长与建设用地的扩张，土地利用变化剧烈，城市生态系统之间的问题日渐突出，故开展土地利用变化与ESV之间的关系研究十分必要。笔者基于乌鲁木齐近20 a土地利用数据，运用CA-Markov模型模拟预测2025年土地利用变化并进行生态系统服务价值(ESV)评估，应用格网分析、空间自相关、Getis-Ord热点区分析等空间统计方法，进一步分析生态系统服务价值与土地利用空间格局的分布特征，更加直观展示出土地利用变化对ESV的影响，预测未来城市化进程中ESV的空间分布格局演变，可为干旱区绿洲城市生态建设与土地资源合理配置提供参考。

1 研究区概况

乌鲁木齐市位于新疆天山山脉中段北麓地区，坐落于准噶尔盆地南缘，北邻福海县，西接五家渠市与昌吉市，东靠阜康市，南与吐鲁番市高昌区、托克逊县、和静县接壤, 是丝绸之路经济带建设的核心城市。地势总体上表现为东南高、西北低，北部主要为低山河谷及冲积平原地域，南部为天山山区(图 1)，起伏较为悬殊，兼备山地城市和平原城市属性。属于中温带大陆性干旱气候区，昼夜温差较大，光照充足，寒暑变化剧烈。

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

基础数据选取1995和2005年Landsat-TM/ETM+和2015年Landsat-8 OLI遥感影像，条带号142-143，行编号29-30，选取6月中旬至9月下旬云量少、图像清晰的遥感影像，空间分辨率均为30 m，数据来源为中国科学院地理空间数据云平台。在ENVI软件支持下，对2005年遥感影像进行条带修复，对3期遥感影像进行图像拼接、裁剪等预处理。采用最大似然分类法与目视解译相结合的方法，参考《全国土地分类标准》和《土地利用现状分类标准》，并与乌鲁木齐的实际土地利用状况对比进行辅助纠正，将土地利用类型划分为耕地、林地、草地、水域、建设用地和未利用地6类，由此获得3期土地利用类型图(图 2)。基于混淆矩阵方法得到混淆矩阵精度评价表，对1995—2015年分类结果进行解译精度评价，1995、2005、2015年土地利用分类成果的总体精度分别为0.90、0.90、0.92，均达到分类标准精度，满足研究要求。

道路矢量数据由1995、2005、2015年遥感影像的第8个波段提取，通过ArcGIS软件空间分析功能中的欧式距离模块分析得到距道路(铁路、高速、公路)、城镇点、水域的距离；DEM高程数据由地理空间数据云(http://www.gscloud.cn)下载拼接获得，空间分辨率为30 m。

以1995和2005年土地利用数据为预测基础图像，利用Markov模块得到1995—2005年转移面积矩阵与概率矩阵，再采用CA-Markov模块预测2015年土地利用类型空间分布情况。随后对2015年预测结果与实际遥感解译土地利用图进行Kappa精度验证，精度高达0.962，说明模拟结果精准度较高，可用于未来土地利用变化预测。将2005和2015年土地利用数据作为预测基础图像，利用Markov模块得到2005—2015年转移面积矩阵与概率矩阵，以2015年为预测的起始时间，以10 a为时间间隔，预测2025年乌鲁木齐市土地利用变化。

2.2 研究方法 2.2.1 CA-Markov模型

Markov模型是一种随机模型，可根据事件的当前状况预测其未来各个时刻变动状况，其本质是对事件发生概率的预测^[16]。该模型主要依据初始时期到另一个不同时期的状态间转移概率与矩阵, 预测未来某时刻的土地利用变化趋势，公式为

$ {S_{t + 1}} = {S_t} \times {\mathit{\boldsymbol{P}}_{ij}}。$

(1)

式(1)中，S_t和S_t+1分别为t和t+1时刻土地利用系统的状态；P_ij为状态转移矩阵。

元胞自动机是一种时间、空间、状态都离散的局部格网动力学模型，具有模拟复杂系统时空演化过程的能力^[17]。模型具体公式为

$ {S_{t + 1}} = f\left( {{S_t}, n} \right)。$

(2)

式(2)中，S为元胞有限、离散的状态集合；n为元胞的邻域；f为局部空间的元胞转换规则。

CA-Markov模型综合了CA模型模拟复杂系统空间变化的能力与Markov模型定量长时间序列预测的优点^[18]，使用CA-Markov模型能更加准确地从时间和空间上模拟土地利用类型的变化情况。因此，以1995—2005及2005—2015年各土地利用类型之间的转化面积作为Markov状态转移概率矩阵元素，根据多准则评价模型选取高程、坡度、距道路距离、距水域距离作为驱动因子，并根据不同土地利用类型转换特点选取限制性因子，从而得到各土地利用类型空间分布概率适宜性图集，为下一步土地利用空间分布的模拟预测提供转移规则。

2.2.2 生态系统服务价值的估算方法

在谢高地等^[19]的研究成果基础上，采用改进型ESV方法的标准当量因子ESV价值量(D，取值为34.065万元·km^-2)，结合乌鲁木齐土地利用具体情况，将单个当量价值与其生态系统单位面积ESV当量表相乘，即可得到不同土地利用类型的单位面积生态服务价值系数(表 1)。

考虑该区域土地利用类型的种类，对计算指标进行调整，耕地与水田采用农田的指标计算，林地采用针叶林的指标计算，草地采用草原与灌草丛的指标平均值计算，水域采用水系和湿地的指标平均值计算，建设用地采用裸地的指标计算，未利用地采用荒漠的指标计算, V_ES的计算公式为

$ {V_{{\rm{ES}}}} = \sum\limits_{k = 1}^n {\left( {{A_k} \times {C_{{\rm{V}}, k}}} \right)} , $

(3)

$ {V_{{\rm{ES}}, f}} = \sum {\left( {{A_k} \times {C_{{\rm{V}}, jk}}} \right)} 。$

(4)

式(3)~(4)中，V_ES为生态系统服务价值, 万元；A_k为第k种土地利用类型面积，km²；C_V，k为第k类土地利用类型的生态系统价值系数，万元·km^-2；V_ES，f为研究区年度生态系统第f项服务功能总价值；C_V，jk为第k种土地利用类型第j项服务功能价值系数，万元·km^-2。

2.2.3 空间自相关

空间自相关分析能发现ESV分布的空间相关性是否具有集聚性，包括全局与局部空间自相关^[20]，全局空间自相关能反映整个区域的地域单元与相邻单元的总体聚集特征，但不能凸显局部区域单元的聚集状态^[21]。局部空间自相关(LISA)可以准确把握空间单元与邻近单元的聚集性与分异特征^[22]，局部Moran′s I值(I_i)计算公式为

$ {I_i} = \frac{{\left( {{x_i} - \bar x} \right)\sum\limits_{j = 1}^n {{\mathit{\boldsymbol{w}}_{ij}}} \left( {{x_j} - \bar x} \right)}}{{\frac{1}{n}\sum\limits_{i = 1}^n {{{\left( {{x_i} - \bar x} \right)}^2}} }}。$

(5)

式(5)中，n为空间单元总个数；x_i和x_j分别为空间单元i和j的ESV观测值；w_ij为基于要素i、j建立的空间权重矩阵。Moran′s I的取值范围为[-1, 1]，采用蒙特卡洛(Monte Carlo)模拟方法检验Moran′s I统计量的显著性，模拟次数为999，检验公式为

$ {Z_I} = \frac{{I - E(I)}}{{S(I)}}, $

(6)

$ S(I) = \sqrt {\mathit{Var}(I)} 。$

(7)

式(6)~(7)中，E(I)为观测变量自相关性的期望值；Var(I)与S(I)分别为方差和标准差；Z_I为标准差的倍数，在给定显著性水平下，Z值＞0时说明存在显著正相关，I值越接近1表明空间单元存在高-高或低-低值聚集，空间差异越小，性质越相似；Z值＜0时说明存在显著性负相关，I值越接近-1，表明空间单元存在高-低或低-高值空间聚集，空间差异越大；Z值等于0表示空间自相关性不显著^[23-24]。

2.2.4 冷/热点区Getis-Ord Gi^*演化分析

Getis-Ord Gi^*指数用于分析ESV变化的高/低空间聚集程度，即ESV动态变化的冷/热点区空间分布格局^[21]。Gi^*指数计算公式为

$ G_i^* = \frac{{\sum\limits_{j = 1}^n {{\mathit{\boldsymbol{w}}_{ij}}} {x_j} - \bar X\sum\limits_{j = 1}^n {{\mathit{\boldsymbol{w}}_{ij}}} }}{{\sqrt[{^{^{^{^{^{^{^{^{^s}}}}}}}}}]{{\left[ {n\sum\limits_{j = 1}^n {\mathit{\boldsymbol{w}}_{ij}^2} - {{\left( {\sum\limits_{j = 1}^n {{\mathit{\boldsymbol{w}}_{ij}}} } \right)}^2}} \right]/(n - 1)}}}}, $

(8)

$ \bar X = \frac{1}{n}\sum\limits_{j = 1}^n {{x_j}} , $

(9)

$ S = \sqrt {\frac{1}{n}\sum\limits_{j = 1}^n {x_j^2} - {{\bar X}^2}} 。$

(10)

式(8)~(10)中，G_i^*为输出统计Z得分；x_j为空间单元j的ESV指数变化量；w_ij为相邻空间单元i和j之间的空间权重。

3 结果与分析 3.1 生态系统服务价值变化分析

基于乌鲁木齐市1995—2025年土地利用分类结果与各土地利用类型的ESV计算指标(表 1)计算各年份不同土地利用类型的ESV变化(表 2)。结果显示，乌鲁木齐市1995、2005、2015、2025年的ESV总值分别为214.23亿、213.29亿、210.70亿和208.06亿元。从表 2、图 3可以看出，乌鲁木齐市土地利用变化引起ESV波动变化，其中草地是研究区ESV的最主要构成部分，其次是未利用地与林地，草地和林地对研究区的ESV有相当重要的生态意义。除建设用地维持稳定增长外，其他用地类型ESV均呈负增长趋势。1995—2025年乌鲁木齐市ESV总体上呈持续减少趋势，减少6.17亿元。各土地利用类型的ESV变化在不同时段存在一定差异，呈波动变化态势。其中2005—2015年为ESV总量急速下降阶段，此阶段对总量减少贡献最大的是耕地，其次是水域。主要原因在于该时段城市不断向外围扩张，大量耕地和水域转化为建设用地。推测到2025年前，研究区的林地和草地ESV值处于持续减少状态，但2015—2025年间下降速度有所减缓，表明近年来乌鲁木齐在城镇化扩展中生态保护意识逐渐加强。

3.2 生态系统服务价值总体分布格局

为了直观定量描述研究区ESV的空间分布特征，将乌鲁木齐市各期土地利用分类数据转入2 km×2 km的格网单元，结合ESV估算模型，分别计算每个格网单元内各土地利用类型的ESV值，通过ArcGIS软件的重分类工具，采用自然断点法将各格网单元ESV值从低到高分为低值、次低值、中值、次高值、高值5个等级(图 4)。如图 4所示，乌鲁木齐市ESV总体分布格局等级差异明显，建成区ESV的中值区被新增低值区不同程度包围，次高值区镶嵌分布于其他区域之中，高值区与次低值区分布相对集中，部分区域从高值区下降到中值区甚至低值区，且新增的低值区大多由原来的低值区向外扩展蔓延转变而成。

ESV高值区域主要分布于乌鲁木齐市南部伊连哈比尔尕山、东部博格达山北坡及乌拉泊湿地、萨尔乔克牧场、米东区林场附近等林地面积较多的地区，在城市周边低矮的山坡林地也有零星分布。低值区域主要分布在米东区北部古尔班通古特沙漠、达坂城区柴窝堡湖及盐湖等未利用地周边，乌鲁木齐中心四区(天山区、沙依巴克区、水磨沟区、头屯河区)周边区域因为受到建设用地不断向北部、南部和东部蔓延扩张及人口聚集的影响，成为低值区域。大量中值区转化成为低值区，ESV值显著下降。其原因在于近几年城市用地规划提出“南控、北扩、先两延、后东进”的发展战略，使得东北方向大量草地转化为建设用地。

乌鲁木齐县及沿柴窝堡—达坂城区形成的南北轴附近大部分区域的ESV值相对稳定，该地处于城市周边区域，受人为因素干扰及城市用地扩张影响程度较小；柴窝堡湖及沙依巴克区雅玛里克山、水磨沟山区、红光山附近区域ESV值在研究时段后期有所上升，这与乌鲁木齐市近年来大力实施“蓝天工程”“园林城市建设”“荒山绿化工程”等措施有极大的关系。2025年的ESV空间分布预测结果显示，次低值区呈增加趋势，随着城镇化的加速，新市区与头屯河区等受建成区域不断扩张的影响，其生态空间分布格局发生较大变化。由此看来，未来建设用地面积将不断增加，但城市扩张速度逐渐得到有效控制，各种生态用地类型之间的变化速率将趋于减缓，这得益于生态保护建设工程的实施、退耕还草还林政策及城市绿化等一系列生态恢复措施的推进，局部地区的生态环境将有所改善。

3.3 空间自相关与冷热点区空间分布格局 3.3.1 空间自相关

选择Moran′s I指数与散点图对研究区各年份的ESV进行空间自相关分析(表 3)，由OpenGeoDa软件绘制的Moran′s I散点图能够定性分辩区域各个格网单元与相邻格网单元ESV的相互关系(图 5)。研究区域4个年份Moran′s I指数的斜率均为正值，整体呈空间正相关，I值均大于0.77，在变化趋势上，Moran′s I指数由1995年的0.79下降到2005年的0.779，表明ESV空间自相关性与空间聚集现象在1995—2005年减弱，2015年上升到0.80，随后2025年下降到0.794，表明研究区ESV的空间分布具有较强的空间自相关性，且分布规律表现出相对一致性，大部分ESV散点分布于第1、3象限，呈现高值区聚集分布以及低值区聚集分布的现象，只有较少的异常值分布于第2、4象限。

然而，通过Moran′s I散点图无法判断各区域具体的空间位置信息以及相关性，因此需要使用Geoda软件研究区显著性水平(图 6)，旨在更为直观地观察ESV单元之间的空间聚集分布情况及具体的空间位置信息。从图 6可见，乌鲁木齐市绝大多数区域相关性不显著，相关性较强的地区主要为高-高值、低-低值聚集区域，这2种类型聚集分布与ESV高值、低值集区域高度重合，存在较强的空间正相关性。高值聚集区主要集中在乌鲁木齐县及沿柴窝堡-达坂城区形成的南北轴附近大部分区域，在空间范围上受人为因素干扰及城市用地扩张影响的程度较小。低值区主要聚集在乌鲁木齐市头屯河区、新市区和米东区附近区域(P < 0.01)，随着年份递增，达显著水平的低值聚集区范围明显扩大(表 3)，与其相邻的单元呈空间正相关性，说明这部分区域经济具有极强的带动性，主要原因是城市建成区从中心城区向北部、东部不断扩张，周围边缘区域内大量生态用地转化为新建区，原有的生态景观发生显著变化。

3.3.2 冷热点区空间分布格局

通过热点分析对乌鲁木齐市4个年份的ESV时空分布差异进行比较，以2 km×2 km格网单元的ESV变化量为对象，计算1995—2005、2005—2015、2015—2025年3个相邻时段各格网单元局域Gi^*统计量，进而识别ESV年际变化的热点区和冷点区分布情况。采用自然断点法将其划分为5类(高热区、次高热区、稳定区、次低冷区、低点区)，生成乌鲁木齐市冷热点空间分布格局(图 7)。

1995—2005年，乌鲁木齐市ESV高热区主要分布在米东区梧桐镇、乌鲁木齐县板房沟乡、水磨沟区红光山和沙依巴克区雅玛里克山附近、达坂城区乌拉泊及南山牧场周围地区。其主要原因为:一方面高热区加强植树造林，使林地面积有所增加；另一方面，受城市绿化及生态恢复措施影响，当地生态环境逐渐改善。而ESV低冷区则主要分布在新市区大部分区域、沙依巴克区西南部、米东区青格达湖乡及乌鲁木齐县周边零星区域。

2005—2015年，乌鲁木齐市ESV高热区及低冷区变化明显，ESV高热区主要分布在乌鲁木齐县甘沟、白杨沟等中部地区及达坂城区乌拉泊及南山牧场周围大部分地区，说明乌鲁木齐市荒山绿化、湿地修复等生态政策卓有成效。ESV低冷区则大幅增加，成片分布在建成区及其周边区域，主要是由于沿柴窝堡湖—盐湖附近地区不断占用草地发展旅游开发等项目。

2015—2025年预测结果显示，高热区与低冷区将频繁转换，稳定区将有所减少，呈现空间聚集效应。高热区主要分布在乌鲁木齐县南部地区，呈现由南部向北偏移趋势，这主要得益于退耕还林、湿地生态保护修复措施等政策的实施。ESV低冷区主要分布在头屯河区北部与米东区南部接壤的区域，在市辖区也有点状零星分布。未来对该区应严格控制建设用地的规模，逐步提升区域ESV值。

乌鲁木齐市ESV变化过程最为明显的是高值、低值区域的空间集聚现象，1995—2025年冷点区和热点区都呈现增加趋势，但值得关注的是热点区范围逐渐大于冷点区，稳定区趋于减少。ESV值变化显示最强烈的高热区主要位于乌鲁木齐县，近年来乌鲁木齐市大力推行湿地生态保护修复措施，加大退耕还林、还草实施力度，使乌鲁木齐生态环境状况明显好转，对乌鲁木齐市ESV增值起到重要作用。而冷点区域由老城区向新城区不断蔓延聚集分布，经济快速发展使中心城市不断占用周边区域的土地，是造成ESV损失的主要原因。

4 讨论与结论 4.1 讨论

土地利用变化深刻影响着自然生态系统结构和服务功能的演变^[25]，从而导致ESV变化。受近年来乌鲁木齐市大力实施“蓝天工程”“园林城市建设”“荒山绿化工程”等措施的影响，研究期间部分耕地与未利用地转移为林草地，对ESV增加贡献率较高，对于研究区东部天然牧场以及森林、河流和滩涂等ESV高的新增生态用地要大力加强保护力度，一旦被破坏，很难重建起原有的生态体系，同时影响生态屏障建设。随着该区域城镇化的不断推进以及土地利用开发强度的加剧，周围边缘地区的生态用地减少，造成ESV降低。从2025年变化趋势预测结果来看，建设用地尚处在不断扩张过程中，在国家提出“三区三线”空间规划体系的背景下，周围边缘区域将继续遭受不同程度的生态扰动。未来乌鲁木齐市土地利用规划应向存量建设用地开发模式转变，要重视建设用地的动态监测与预警，避免城市向边界的不断扩张，减少城市化对生态系统和生态本底的破坏。笔者基于空间统计与分析来探索未来土地利用变化对ESV空间格局的影响，并结合热点分析定量分析刻画了土地利用背景下的ESV未来时空变化，国内外众多学者在ESV时空变化作了大量研究，但对于CA-Markov预测结合热点分析来探究微观区域ESV的时空格局演变略有不足，该研究提供了一种新的研究思路与框架。

目前生态系统服务价值估算可分功能价值法和当量因子法2类。因为当量因子法具有可操作性、直观性，评估结果可比性高^[19]。与部分学者对于乌鲁木齐ESV的估算结果相比，笔者的计算结果存在一定差异性，其原因是对不同系统生态服务价值单价的计算方法不同，如阿依吐尔逊等^[26]将土地利用类型划分为7种，赋予其生态服务价值系数；李芝灵^[27]根据生态类型着重从5类生态功能出发，对乌鲁木齐市ESV进行评估；张浩等^[28]结合其他学者的相关研究成果，建立评价模型来计算乌鲁木齐市ESV。笔者直接采用谢高地等^[10]的评估方法与模型对乌鲁木齐市ESV进行估算，因此与其他学者的ESV计算结果存在差异，一定程度上丰富了ESV评估方法在干旱典型绿洲的应用。

由于生态系统的复杂性、动态性和非线性特征，当量因子法的估算结果存在一些局限^{[9, 29]}，预测模拟中土地利用类型分类的精度也会影响模拟结果。首先，笔者研究中未充分考虑当量因子法运用于不同区域时生态系统类型修正系数的差异，在一定程度上影响了ESV的评估精度。后续研究应采用植被净初级生产力、绿量指数等指标进行校正，获取本地化的当量因子和标准当量因子生态系统服务价值量，使评估结果更加符合现实。其次，该研究计算的ESV受生态服务当量与土地利用类型数据精度的影响，在今后的研究中需使用高分辨率的遥感影像获取土地利用数据，以提高计算精度。最后，该研究采用的土地利用数据仅分为6类，训练样本构建时仅获取研究区基础地理数据及部分自然影响因素，所选研究样本较少，约束条件不够全面，虽然对于小区域的土地利用变化和预测模拟来说具有一定的精度和可靠性，但在后续研究中需进一步细分土地利用类型，并加入社会经济因子来定量化模拟未来的土地利用情景。

4.2 结论

(1) 1995—2025年乌鲁木齐市ESV呈减少趋势，由214.23亿元减少到208.06亿元，净减少6.17亿元。ESV减少的主要原因在于耕地、林地、草地面积减少，其中林地和草地是研究区域ESV的最主要构成部分，未来要更加注重生态保护和生态恢复建设。

(2) 总体上乌鲁木齐ESV总体空间分布格局区域差异明显，ESV高值区域多年来主要分布在乌鲁木齐县林地面积较多的地区，ESV低值区域集中分布于乌鲁木齐中心四区(天山区、沙依巴克区、水磨沟区、头屯河区)周边区域，随着大力推进“园林城市建设”“荒山绿化工程”等政策的实施以及林草地生态系统的恢复，ESV在研究时段后期有所上升。

(3) 研究区ESV具有显著的空间正自相关性与空间聚集特征，高-高、低-低2种类型聚集分布分别与ESV高值、低值区域高度重合。1995—2025年，冷点区和热点区都表现为增加趋势，热点区域的增加与近年来乌鲁木齐大力推行生态保护及修复等多种措施有极大的关联性，冷点区的增加与城市建成区的增量扩张蔓延及土地开发等一系列人类活动有密切关系。