土地利用/覆被变化(land use and land cover change，LUCC)是地球表面发生改变的表现之一，是影响区域生境质量变化的关键因素^[1]，而土地利用结构预测也在建设用地扩张的背景下成为全球关注的热点问题^[2]。生境质量是指生态系统在一定的时间和空间中，能够提供个体和种群持续适宜的生存条件的能力^[3]；其作为一切生态系统服务的重要基础，更是提升人类福祉的关键环节^[4]。近年来，众多学者运用SolVES，InVEST，FLUS，CA-Markov，MaxEnt等^[5-7]模型集中对城市、流域、县域尺度进行生境质量评价和土地利用变化预测。早期学者们主要以动植物为对象构建评价指标体系研究生境质量，随着数理模型方法及各种软件的发展与应用，生境质量的时空动态变化受到关注^[8]，InVEST模型的生境质量模块(habitat quality)是目前使用最多的生态系统服务评估模型^[9]。刘智方^[10]、刘春芳^[11]、王耕等^[12]结合土地利用转移矩阵和景观格局指数分析福建省、榆中县黄土丘陵区、老铁山自然保护区的生境质量时空格局，发现生境质量与社会经济发展、人类干扰强度存在一定联系。CA-Markov模型结合了元胞自动机在空间上模拟复杂变化的能力和马尔科夫链在数量上长期预测的优势^[13-14]，能够有效提高土地利用时空格局变化的模拟精度^[15]，得到国内外学者一致认可。Praveen Subedi^[16]、Nourqolipour等^[17]利用CA-Markov模型预测2015年美国佛罗里达州Saddle Greek流域及2020年马来西亚Kuala Langat地区油棕种植园的土地利用变化趋势；同时褚琳^[18]、武丹等^[19]结合InVEST和CA-Markov模型研究并预测武汉市和宁夏回族自治区中部干旱区未来土地利用结构及生境质量的时空演变特征。

中共十八大报告明确提出建设生态文明和美丽中国，需要加强环境管理和保护，坚持严格准入、限制开发，构建生态安全保护屏障。但目前学者们对生物资源丰富、发挥涵养水源作用的重点生态功能区生境质量研究较少，且没有将生境质量结果与土地利用类型中林地的空间相关性进行单独探讨。因此，本文以湖南省重点生态功能区且林地居多的南岳区作为研究对象，在采用土地利用转移矩阵及InVEST模型分析土地利用变化和生境质量时空格局的基础上，引入双变量空间自相关方法探究生境质量结果与林地两者的空间分布关系，并运用CA-Markov模型预测未来土地利用结构和生境质量变化趋势，为该地区的土地利用管理和生态保护提供科学依据。

1 研究区概况与数据来源 1.1 研究区概况

南岳区地处湖南省衡阳市北部(112°33′44″—112°46′34″E，27°11′29″—27°20′5″N)，总面积179.27 km²，现设一乡(寿岳乡)一镇(南岳镇)一街(祝融街道)，总人口为7.03万。该区属于大陆性中亚热带季风湿润气候，年均气温17.5 ℃，年日照1 628.6 h，年平均降水量2 153.4 mm；中部是衡山山脉的主峰祝融峰，由中部向东南与西北两侧逐渐放缓变平，形成中部高、两翼低的倒V地形。生态环境优良，资源丰富；有银杏、长白松、玉龙蕨等国家一级保护植物21种；黄腹角稚等珍稀濒危动物28种；森林覆盖率高达76.2%。南岳区作为生态保护屏障地区，是湖南省“重点生态功能区”“国家级风景公园”“国家级自然保护区”，但近年来该地区占用林地资源发展建设问题较为严重。在生态文明背景下，探究并预测重点生态功能区的生境质量时空特征对促进资源环境与社会经济的协调发展具有重要借鉴和指导意义。

1.2 数据来源

研究采用衡阳市自然资源和规划局提供的2009，2012，2015，2018和2021年5期土地利用变更数据，参考土地利用现状分类标准(GB/T21010-2017)及研究区实际情况将所有地类归并为耕地、林地、园地、草地、水域、建设用地、其他土地7个一级地类(表 1)。DEM高程数据从地理空间数据云平台(http://www.gscloud.cn/)下载，分辨率为30 m。人口数据来源于《南岳区统计年鉴》，生境质量均值、空间自相关指数及土地利用预测数据由ArcGIS，InVEST，IDRISI，GeoDa等软件进行运算处理后得到。

2 研究方法 2.1 土地利用转移矩阵

利用ArcGIS中的相交工具，将2009，2012，2015，2018和2021年5期土地变更调查数据进行叠加分析并导出Excel表，借助数据透视表构建土地利用转移矩阵，分析南岳区不同时段各土地利用类型之间的转入、转出去向及面积大小。

2.2 InVEST模型

InVEST模型通过模拟不同土地覆被情景下生态系统服务的变化情况，将生态系统服务功能评价结果以栅格图的形式表现出来，其主要包括陆地、海洋和淡水生态系统评估3大模块^[20]。InVEST模型的Habitat Quality模块是以土地利用数据为基础，根据各胁迫因子数据、生境类型的生境适宜性以及对胁迫因子的敏感性计算出生境质量指数来表征研究区土地利用动态变化与生境质量状况的关系^[21]。其计算的具体公式为：

$ Q_{x j}=H_j\left[1-\left(\frac{D_{x j}^z}{D_{x j}^z+k^z}\right)\right] $

(1)

式中：Q_xj为土地利用类型j在栅格x中的生境质量；D_xj为土地利用类型j在栅格x中的总威胁水平；H_j表示土地利用类型j的生境适宜性。k为半饱和系数, 取Q_xj试运行1次后所得到的最大值的1/2；z为模型默认参数。

$ D_{x j}=\sum\limits_{r=1}^R \sum\limits_{y=1}^{Y_r}\left(\omega_r / \sum\limits_{r=1}^R \omega_r\right) r_y i_{r x y} \beta_x S_{j r} $

(2)

式中：r为威胁因子；y为r威胁的栅格数；Y_r为威胁因子r所占的1组栅格数；ω_r为威胁因子r的权重；r_y为栅格y的威胁因子值；i_rxy为栅格y的威胁因子r对生境栅格x的威胁度；β_x为栅格x的可达性水平；S_jr为生境类型j对胁迫因子r的敏感性；i_rxy由以下计算公式得到：

$ i_{r x y}=1-\left(\frac{d_{x y}}{d_{r \max }}\right) $

(3)

式中：d_xy为栅格与栅格之间的直线距离；d_rmax为胁迫因子的最大影响距离。

生境质量的取值范围为0~1；分值的高低反映了在威胁源的影响下生境质量的好坏，得分越高，表明其生境质量越好；反之，越差。依据InVEST模型手册中生态因子的划分标准并结合研究区实际情况，本文选取了水田、旱地、公路用地、城镇用地、村庄用地5种人类活动较为集中、容易对生境质量产生影响的土地利用类型作为生境威胁因子；并选取了受人类活动干扰较少的有林地作为生境。以模型推荐值为基础，参考包玉斌^[22]、张学儒^[23]、吴哲^[24]、王军等^[25]的研究以及专家意见进行赋值，确定各威胁因子及其最大影响距离、权重及衰退类型(表 2)、不同土地利用类型的生境适宜度及其对不同威胁因子的相对敏感程度(表 3)。通常情况下，越接近自然生态环境系统的生境适宜性越高，敏感性相对较低，用0~1表示。该值越大，就代表生境对威胁因子敏感性越强。

2.3 双变量空间自相关

空间自相关分析可反映变量在空间上的分布特征及集聚性，分为全局空间自相关和局部空间自相关。Anselin提出的双变量空间自相关(Bivariate Moran’s I)能够探索多个变量在空间分布上的相关性^[26]。

双变量全局空间自相关是验证相关变量在整个研究区域的空间分布聚集趋势^[27]，从土地利用类型和生境质量空间分布结果来看，它们之间具有高度的相关性，运用全局Moran’s I分析土地利用类型中比例最大且作为生境的林地与生境质量结果的空间关联，计算公式为：

$ I=\frac{\sum\limits_{i=1}^n \sum\limits_{j=1}^n W_{i j}\left(x_i-\bar{x}\right)\left(y_j-\bar{y}\right)}{S^2 \sum\limits_{i=1}^n \sum\limits_{j=1}^n W_{i j}} $

(4)

式中：n表示研究的空间单元个数；W_ij为空间单元i和j的空间权重矩阵；x_i，y_j分别是自变量、因变量在空间单元i和j的观测值；s为观测样本x_i的方差。Global Moran’s I取值范围为-1~1, 若指数大于0, 则表示空间单元存在正相关关系, 且越接近1, 聚集程度越高；若指数趋向于0, 则表示该空间单元是随机分布, 不具备空间相关性；若指数小于0, 则表示空间单元属性呈负相关, 且越接近-1, 离散程度越高。

双变量局部空间自相关是衡量研究区域内各变量单元与其邻近单元之间的相关性，局部Moran’s I计算公式为：

$ I_i=z_i \sum\limits_{j=1}^n w_{i j} z_j $

(5)

式中：I_i代表空间单元自变量与因变量的局部空间自相关指数；z_i，z_j都是方差标准化的单元观测值。根据Local Moran’s I指数形成4种聚类模式, 并由此组成基于z检验的LISA(local indications of spatial association)分布图来呈现局部区域单元的关联性^[28]。聚类模式分为高高(HH)聚集, 即林地与优等生境质量区；低低(LL)聚集, 即非林地与低生境质量区；高低(HL)聚集, 即林地与低生境质量区；低高(LH)聚集, 即非林地与优等生境质量区。

2.4 CA-Markov模型 2.4.1 元胞自动机模型

CA(元胞自动机)是由元胞及其状态、邻域和转换规则4部分组成的一个时间、空间、状态都呈现离散特征的网络动力学模型，它按照定义的转换规则进行更新，模拟复杂系统时空演化过程^[29]。其表达式为：

$ S_{(t+1)}=f\left[S_{(t)}, N\right] $

(6)

式中：S为元胞有限、离散的状态集合；t，t+1表示不同时刻；N为元胞的邻域；f为局部空间的元胞转换规则。

2.4.2 马尔柯夫模型

Markov链是基于随机过程理论预测事件发生概率常用的一种方法，侧重于土地利用变化数量上的预测，但无法表达空间分布^[30]。它是通过矩阵分析当前状态以及变化趋势，对未来的可能性进行模拟，计算公式为：

$ S_{t+1}^{\prime \prime}=P_{i j} \times S_t^{\prime \prime} $

(7)

式中：S_t″，S_t+1″分别为t, t+1时刻的系统状态；P_ij为状态转移概率。确定状态转移概率矩阵是该模型研究的关键步骤, 计算公式为：

$ \boldsymbol{p}_{i j}=\left[\begin{array}{cccc} P_{11} & P_{12} & \cdots & P_{1 n} \\ P_{21} & P_{22} & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots \\ P_{n 1} & P_{n 2} & \cdots & P_{n n} \end{array}\right] $

(8)

式中：n表示土地利用类型; P_ij表示土地利用类型i转变为土地利用类型j的转移概率, 矩阵中每行为土地利用类型向其他类型转出的概率, 每列为土地利用类型由其他类型转入的概率。

2.4.3 CA-Markov模型

CA-Markov模型将空间维度和时间维度分析优势结合起来，可以更加精准地模拟长期且复杂的土地时空变化^[31]。本文使用IDRISI软件中的CA-Markov模块来预测研究区未来土地利用格局。主要步骤为：①将2009，2012，2015，2018和2021年5期土地利用栅格数据转换成IDRISI软件可识别的.rst栅格数据格式，然后对其进行重分类处理；利用Markov模块得到土地利用转移概率矩阵、面积矩阵和一系列条件概率图像。②结合研究区高程、坡度等实际情况制定CA模型转换规则，将各地类的适宜性图像通过Collection editor制成图集。③设置起始时刻和CA循环次数，使用5×5的滤波器模拟土地利用格局。运用Crosstab计算工具对预测图与实际图进行精度检验，得到2012，2015，2018和2021年的kappa系数均大于0.91；说明模拟结果与实际情况较为吻合，可运用此模型对未来土地利用格局进行预测。

3 结果与分析 3.1 土地利用变化分析

从土地利用变化数量上看，2009—2021年，研究区林地和耕地减少的面积最多(表 4)。其中林地面积共减少73.88 hm²；耕地面积由2 329.82 hm²减少至2 256.12 hm²，减少率为3.16%。然而，因研究期内建成区的扩张，南岳重点生态功能区建设用地面积逐年增长，由1 241.15 hm²增加至1 438.10 hm²，共增加196.95 hm²，增长率为15.87%。从土地利用变化空间分布上看(图 1)，林地遍布南岳区的北部、南部和西部，约占研究区总面积的75%。其次是建设用地和耕地，主要位于研究区的东部和东南部，呈“集中—点带”式分布。2009—2021年，土地利用类型的空间变化主要表现在原建设用地附近的林地和耕地转为建设用地。

从土地利用转移矩阵分时段来看(表 5)，2009—2012年，有99.38 hm²土地转为建设用地，主要由40.31 hm²林地和28.29 hm²耕地转入，分别占新增建设用地面积的40.58%，28.47%；这一时期建设占用了较多的林地资源，主要与当时南岳推动旅游经济发展，完善基础设施建设有关。2012—2015年，建设用地面积仍然保持较快的增长趋势，净增长72.30 hm²，有29.46 hm²林地和22.08 hm²耕地转入为建设用地，分别占新增建设用地面积的40.75%和30.54%；这一阶段建设占用林地的面积得到减少，但建设占用耕地的情况没有得到大幅缓解；

这是由于耕地大多分布在建成区周围，基础条件较好，便于开发利用。2015—2018年，建设用地扩张速率得到明显降低，仅为2009—2012年新增建设用地总面积的25.25%。2018—2021年，在严格管控的条件下，建设用地面积仅增长20.08 hm²，年增长率进一步下降。从土地利用类型的整体变化趋势来看，南岳作为林地居多的重点生态功能区，2009—2021年，生境适宜度最高的林地变化最大，其生境质量也将发生相应变化。

3.2 生境质量时空变化

根据南岳区实际情况和相关研究成果，采用自然断点法将研究区的生境质量结果(Q)划为4个等级：低(0≤Q＜0.2)、中(0.2≤Q＜0.5)、高(0.5≤Q≤0.8)、优(0.8＜Q≤1)。各年份的生境质量等级土地面积、比例及空间分布情况分别如表 6和图 2所示。通过栅格计算得到2009，2012，2015，2018和2021年研究期间南岳区的生境质量均值分别为0.617 6，0.608 7，0.601 1，0.598 8和0.596 9，呈“先快后慢”的持续下降趋势，降幅分别为1.44%，1.25%，0.38%，0.32%。

从时间尺度上看，2009—2021年南岳区生境质量变化主要表现为低等和中等生境质量区面积缓慢增加，优等生境质量区面积逐渐减少2.35%。2009—2012年，优等生境质量区比重降低了1.3%，低、中、高等生境质量区分别增加了0.38%，0.37%，0.55%；其中，高等生境质量区面积增加最多，为98.01 hm²。主要是由于这一时期建设用地面积增长速率较快，林地面积显著减少，生境质量下降较多，优等生境质量大多转为高等生境质量。2012—2015年，研究区生境质量变化趋势与2009—2012年基本保持一致，但变化速率有所下降。2012年11月党的十八大报告指出要“加大自然生态系统和环境保护力度”之后，建设用地面积增长趋势得到减缓，优等生境质量区面积变化也随之减少；但即使对人类干扰导致的环境恶化进行了一系列生态修复工程措施，也需要一定的时间才能见效。所以低等生境质量区面积增幅最大，为0.44%。2015—2018年和2018—2021年，随着土地利用变化中建设占用林地资源等情况的减少，生境质量退化速率也进一步降低。

从空间格局上看，研究区整体生境质量处于较高水平，以高等和优等为主，约占南岳区总面积的68%，呈现出“西部高、东南低”的空间分布特征。其中，优等生境质量区主要位于山脉西侧的林地部分，该地区的原始次生林为南岳重点生态功能区的定位奠定了重要基础。高等生境质量区的土地利用类型仍然以林地为主，但是由于距耕地、建设用地等威胁源较近，这部分林地的生境质量偏低，主要分布在山脉东侧的南岳镇以及寿岳乡零散耕地周围。中等生境质量区面积最少，呈“带状”分布在低等生境质量区与高等生境质量区两者的耕地及林地中；低等生境质量区的土地利用类型以建设用地为主，该区域地势平缓、人类活动频繁，导致生境质量偏低。

从空间变化上看，随着建成区的扩张，低生境质量区逐渐向其附近的中生境质量区扩散，使得该生境类型范围增大；中生境质量区沿着高生境质量区向西北方向递进；而高生境质量区逐渐向其西侧的优生境质量区蔓延。生境质量的等级变化主要发生在各等级接壤地带，出现了上一等级生境质量向其下一等级生境质量转换的现象。根据南岳区土地利用类型和生境质量空间分布图，可发现生境质量分区与土地利用类型的空间分布有着相似之处，进一步探索他们的内在联系，有利于加强重点生态功能区的保护与管控。

3.3 林地与生境质量的空间相关性

高等和优等生境质量区分布较广，且大多为林地；而低等生境质量区以建设用地和耕地居多。将研究区地类分为林地和非林地两部分，采用GeoDa软件的queen contiguity空间矩阵，测算出2009—2021年5个时间段南岳区林地与生境质量的双变量全局空间自相关指数，分别为0.204，0.207，0.209，0.213和0.216，且在1%的显著性水平下其正态性假设检验统计量Z值均通过了检验。双变量全局空间自相关检验表明林地与生境质量具有显著的空间依赖效应，2009—2021年研究区林地与生境质量在空间分布上并非完全随机，具有明显的聚集性，且其正相关性逐渐增强。由于近年来人类活动导致部分林地被建设占用，生境质量的降低与林地的减少存在着一定联系，为了揭示研究区内空间自相关性的具体位置，选择局部Moran’s I分析出现聚集区域的特征(图 3)。

从时间尺度上看，2009—2012年，林地被建设占用，导致高高聚集区域减少，转变为不显著；2012—2015年，生态文明建设推行的土地整治、退耕还林等措施，使林地增加，生境质量提升，高高聚集区域稍有增加；2015—2018年和2018—2021年，由于各项环境保护政策的实行，土地利用类型与生境质量变化较慢，研究区各聚集类型的空间分布变化也较小；2021年较2009年，高高聚集和低低聚集区域均有略微增加，说明部分林地转为非林地、生境质量降低；同时生态修复工程得以见效，部分林地生境质量升高。从空间格局上看，2009—2021年，高高聚集是南岳区林地与生境质量等级在空间位置上的主要关系，分布于山脉西部海拔较高且难以开发利用的生态区；低低聚集分布的区域也相对较广，主要在生境质量较低的耕地、建设用地等地势平缓的土地利用类型上集中分布；高低聚集地区的土地类型为林地，但是由于城镇开发的辐射作用，导致其生境质量偏低；低高聚集零散分布在优等生境质量区的园地、草地等其他地类中；高高聚集与低低聚集过渡带存在大量的聚集不显著地区，此处生境质量与林地的存在无直接关系。

3.4 土地利用结构与生境质量预测

根据空间自相关分析，发现重点生态功能区土地利用类型中的林地影响着生境质量的高低。通过验证模拟精度后，以2021年土地利用数据为基础，把生态红线范围的区域设置为限制因素；运用CA-Markov模型对研究区LUCC进行预测。并根据预测结果，得到2025—2030年南岳区土地利用变化率(表 7)及生境质量空间分布图(图 4—5)，以便引导政府未来管控。

从时间和数量上看，预测显示2025年研究区耕地、林地、园地、草地、水域、建设用地和其他土地的面积分别为2 242.08，13 397.73，263.61，51.54，193.26，1 500.86，278.82 hm²。其中，建设用地将从2021年的1 438.10 hm²增加至2025年的1 500.86 hm²；增加62.76 hm²，增长率为4.36%；而其他地类中，林地减少的面积最多，耕地减少的比重最大。低、高等生境质量区将呈现“圈层”扩散趋势，优等生境质量区向高等生境质量区转变、中等生境质量区向低等生境质量区转变，低、高等生境质量区面积将进一步增加。预测显示2030年研究区建设用地面积将持续增加到1 563.53 hm²，从转变类型上看，未来建设用地面积增加的主要来源为原建设用地周围的林地和耕地；在不考虑各项政策因素的影响下，随着建设用地面积的增加，未来南岳区生境质量将整体下降。

从空间格局上看，研究区西部将依然以林地为主，耕地在整个区域内呈破碎状分布，而草地、园地、水域和其他土地将继续以零散状态分布在整个区域内。结合南岳区“一山一核一轴一区”的发展格局及规划，未来土地利用变化将主要发生在地势较低的东南方；建设用地以“点状结合”的形式向北部和南部的耕地以及西部的林地扩散。研究区东南部生境质量受土地利用类型变化的影响将有所降低，中等生境质量区将继续向西侧高等生境质量区蔓延。但在重点生态功能区的定位背景下，各级生境质量变化速率将进一步降低。

4 讨论与结论 4.1 结论

本文以南岳重点生态功能区为例，基于2009，2012，2015，2018和2021年5期土地利用现状数据，借助ArcGIS软件、InVEST模型、双变量空间自相关和CA-Markov模型分析研究区2009—2021年土地利用结构和生境质量时空特征，探讨生境质量高低与林地之间的聚集关系，并预测该区域2025和2030年土地利用结构和生境质量的时空变化趋势，得出以下结论：

(1) 研究期内南岳区各地类间的转移以林地和耕地转为建设用地为主，林地减少的面积最多，耕地下降比重最大；2009—2021年建设用地面积增加196.95 hm²，增长率为15.87%，且45.37%，36.25%分别由林地和耕地转入。

(2) 研究区平均生境质量指数均超过0.59，优等生境质量区所占比例最大，分布较为集中；2009—2021年整体生境水平有所下降，生境质量指数由0.617 6下降至0.596 9；主要表现为优等生境质量区减少，低等生境质量区增加。

(3) 南岳区林地与生境质量的全局Moran’s I均大于0.2，具有显著的正相关性，且林地对生境质量的正向影响逐渐增强；其中，高高聚集为主要表现类型，但低低聚集区域的增加表明林地的减少会导致生境质量的降低。

(4) 2025和2030年研究区土地利用类型变化将主要发生在东南方地势平坦的城镇中心，低等生境质量区向周围中等生境质量区扩增；建设用地面积持续增加，林地和耕地面积继续减少，整体生境质量进一步降低。

4.2 讨论

本文研究发现在过去和将来，南岳作为旅游景区和重点生态功能区，生态环境保护都格外重要。但建设用地增加、林地减少导致的生境质量降低已成为难以避免的现实问题。从2009—2021年的土地利用转移矩阵来看，该区域土地利用变化以林地、耕地和建设用地3种土地类型为主；采用InVEST模型可将其生境质量评价结果直观地展示出来，为行政部门分区管控提供科学参考；探究林地与生境质量的空间关系，并预测研究区未来土地利用结构和生境质量特征，可进一步引导优化土地利用空间布局。

由于模型存在一定的局限性，生境质量中的威胁源不仅会受到人类社会经济活动的影响，还会受到自然环境中不确定性灾害的干扰。自然、社会、经济等诸多因素的存在，威胁因子和生境敏感性等参数的设置通常会受主观判断的影响，并且多种威胁的集体影响可能远远大于各种威胁水平之和，从而导致试验结果与实际状况存在偏差。对于未来土地利用结构及生境质量的预测，也会因政策及规划的不同而产生误差；但是作为湖南省国土空间规划确定的重点生态功能区，在追求经济效益的同时，更应该注重生态效益；适当控制建设用地扩张，加强土地资源管理，实现南岳区经济发展与生态保护“双赢”目标。