基于CA-Markov和InVEST模型的土地利用格局变化对生境的影响研究——以北京浅山区为例

邹天娇

倪 畅

郑 曦*

经济发展和城市扩张伴随着大规模土地利用方式、强度及格局的改变[1]。这会影响生境斑块间的物质循环、能量流动和信息传递，使生境分布、面积和功能等都发生变化[2]。评价变化多以生境质量作为主要指标，即生态系统可以给物种的个体种群、群落及人类持续生存和繁殖所需条件的能力[3]。

近年来针对生境变化的研究，有学者选取过小城镇[4]、河流[5]、自然保护区[6]等进行评估。随着3S技术和生态模型发展，方法趋于空间化、定量化和精细化，使用较多的包括生境适宜性模型HIS[7]、IDRISI软件中生物多样性评价模块[8]、SolVES模型[9]、InVEST模型[10]、ARIES[11]和MIMES[12]等。

其中InVEST模型被认为是目前最成熟、应用最多的生态系统服务评估模型[13]，它是自然资本项目(the Nature Capital Project)的成果，使用成本低、评估准确度高、空间分析能力强[14-15]。其中生境质量模型(Habitat Quality Model)更是广泛应用于研究基于生境胁迫而评估生境质量[16]。Terrado等发现InVEST模型中，生境质量模块与生物多样性显著相关[17]；Leh等用InVEST模型评估了非洲西部地区的生境质量[18]；包玉斌等分析了黄河湿地自然保护区的土地利用变化对生境的影响[19]；巩杰等研究过甘肃白龙江流域生境质量时空分异[20]。

针对生境适应性管理，需要对未来的土地利用变化进行预测。目前最具有地理学代表意义的是以元胞自动机(Cellular Automata，CA)预测土地利用变化的模型[21-23]，GIS的发展更推动了其进步[24]。在国外，Batty等提出城市发展动态模型，并模拟了美国纽约州西部城市布法罗的城市扩张过程[25]；Nourqolipour等以马来西亚油棕种植园为研究对象，利用模型模拟该地区2020年土地利用演变趋势[26]等。我国关于CA模型研究热潮始于20世纪90年代末期，后有孙贤斌等应用CA-Markov模型，研究挠力河流域不同时段土地利用对湿地的干扰[27]；赵建军等运用3期土地利用数据，较好地预测出了向海保护区2010、2018年的土地覆被类型等[28]。

InVEST模型和CA-Markov模型虽然各自都有较久的使用历史，但将二者结合的研究相对较少。本文拟结合InVEST模型和CA-Markov模型，先对1999和2017年北京浅山区生境质量、生境退化度等进行分析，并采用Logistic回归模型进行景观格局变化驱动力研究，模拟预测出浅山区2035年的生境情况，以期为相关部门进行区域生态规划和出台土地管理政策提供科学依据。

1 研究区概况与数据来源

1.1 研究区概况

北 京 位 于 华 北 平 原 西 北 部(115.7～117.4°E、39.4～41.6°N)，面积16 410.54km2。地势西北高耸，东南低缓。本次所研究的浅山区位于平原城区与西北山脉的过渡地带(图1)，面积约为3 917.84km2。研究区主要以海拔和地形作为遴选依据，区内海拔高度100～300m，以丘陵和台地地貌为主。气候上，浅山地区沿西北山麓形成长条形暖温带，沿太行山、燕山山脉的迎风坡形成一条弧形多雨带；风向受地形影响变化多样。因位于生态交错带内，生态系统多样化，物种种类和数量都胜于邻近区域，并且还是地表水源涵养、地下水补充、农作物生产和防风固沙的重要区域。

1.2 数据来源

1.2.1 土地覆盖类型分类

对1999、2017年的遥感影像进行处理，所用影像来源于地理空间数据云平台(http://www.gscloud.cn/)，选择云量少、月份接近的夏季数据，1999年采用Landsat7 ETM SLC-on影像，2017年采用Landsat8 OLI_TIRS影像，空间分辨率均为30m。

图2 研究思路框架

在ENVI5.0中对这2期遥感影像进行预处理与影像拼接，使用研究范围边界裁剪，得到研究区遥感影像数据。根据土地资源用途、利用方式和覆盖特征，参考《土地利用现状分类》(GB/T 21010—2017)和相关研究文献[29]，结合研究区实际情况，将土地分成5种类型，分别为水域、建设用地、林地、耕地和未利用地。采用监督分类与人工目视解译结合的方法，对研究区1999和2017年2期景观类别和格局信息进行解译提取。

1.2.2 驱动因子

结合实际情况选取与景观格局变化关联较大的指标。驱动力因子选取高程、坡度、距主要道路距离、距居民点距离和距河流水库距离5个因子。DEM数据来源于地理空间数据云平台(http://www.gscloud.cn/)，经投影转换、裁切得到研究区高程数据，并采用ArcGIS的Slope功能提取获得坡度。交通道路、居民点和河流水库分布图均来自地图，地图数据是从国家基础地理信息系统网站(http://nfgis.nsdi.gov.cn/)下载的1:400万中国地图和中国行政区划图，进行对应北京市环境图层的提取。采用GIS中Buffer进行缓冲分析和赋值，获得交通道路、居民点和河流水库的区位缓冲图。

1.2.3 威胁因子

参考前人对北京土地利用格局变化和对区域尺度生境质量的评估研究[29]，将城镇用地、农村居民点、主要交通干道、耕地和矿区定义为生境的威胁源(表1)。城镇用地、耕地的空间分布图来源于对土地覆被类型分布图重分类，交通道路和居民点的空间分布图均来源于地图，并确定了解译出的5种土地利用类型对不同威胁源的相对敏感程度(表2)。

2 研究模型与研究方法

2.1 研究框架与技术路线

提取研究区1999、2017年的景观类型信息，采用CA-Markov模型模拟浅山区2035年自然增长情景下的景观格局情况；结合InVEST模型，评价2期景观格局分布下的生境质量和变化，以及2035年模拟情景下的生境分布情况，总框架如图2所示。

表1 威胁源的影响范围及权重

表2 生境适宜度及其对不同威胁源的相对敏感程度

2.2 研究用到的模型

2.2.1 马尔柯夫模型

马尔柯夫(Markov)预测法来源于数学家马尔柯夫的随机过程理论，是一种特殊的随机运动过程[30]。该模型通过计算各种不同状态的初始概率和状态见转移概率关系，确定时间变化过程中的景观格局演变趋势，最终预测出景观格局的变化。

景观类型之间转化无向性和双向性强，类型和数量不断随机改变，难以用数学方式精确描述，符合了马尔柯夫研究前提，关键点是确定好每种景观之间转化的概率。景观格局变化过程中，每种状态间转化会形成一个转移概率矩阵，公式如下：

式中，n为景观类型数目；P x y为景观x类型变为y类型的概率。Pxy需要同时满足0≤Pxy≤1，且所有Pxy求和为1。

2.2.2 元胞自动机模型

元胞自动机(Cellular Automata，CA)模型，属于不连续的时空动力学模型[31-33]。每个位于元胞空间的细胞，其特定状态都是有限个的，会按所定义的局部规则同步更新。相互作用之下就产生了动态演化系统。CA模型的表达式如下：

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

2.2.3 CA模型与Markov模型的耦合

Markov模型和CA模型都有其局限性：Markov模型在预测土地利用变化时，只能预测到数量变化，而不能预测空间分布变化。CA模型主攻复杂系统空间演变，但没有Markov模型长期预测优势。二者结合则可以综合考虑到各种因素影响和土地利用的历史趋势，从而更精确地模拟土地利用变化的时空演变过程。

在栅格图中可以把一个个像元看作一个个元胞，元胞状态表示景观类型。在GIS中用转换面积矩阵和条件概率图像运算，探寻元胞状态的转移，据此模拟景观格局变化[34]，步骤如下。

1)以叠置分析获取景观类型转移概率矩阵、转移面积矩阵和其他各种条件概率图像(这些图像来源于转移概率矩阵，表示每个像元在下一时间点上被某一景观类型覆盖的概率)。

2)构造CA滤波器。分析邻居离元胞距离的远近，建立有显著空间意义的权重因子，使其作用于元胞，确定元胞状态改变情况。本研究采用5×5滤波器，即表示认定一个元胞周围的5×5个元胞组成的矩形空间会对该元胞状态改变产生显著影响。

3)确定起始时刻和CA循环次数。以研究区域1999—2017年景观类型转移矩阵为基础，将2017年作为预测的起始时刻，每年迭代一次，CA模型迭代次数定为18，在自然增长情景下模拟北京市浅山区2035年的景观格局。

为保证模拟结果的准确，需要验证模型，本研究中用Kappa指数检验了模型在模拟格局变化中的精度，公式如下：

式中，P0代表模拟正确的比例；Pc代表随机情景下模型模拟正确的比例；Pp代表理想分类情景下模拟正确的比例。

再运用IDRISI软件中的CROSSTAB模块，输入研究区2017年实际景观类型图与预测景观类型图，以对比检验模型精度。结果2017年Kappa系数结果为0.798 5，代表模拟效果良好，可使用验证后的CA-Markov预测2035年的景观类型。

2.2.4 InVEST模型

在InVEST模型中，Biodiversity模块用栖息地质量的优劣来表征生物多样性的弹性、恢复性、广度和深度的变化。其计算结果生境质量和生境稀缺性能用来反映生物多样性。

运行InVEST模型，按模型自带参数在Biodiversity模块里设置对应数据，模型会将景观类型信息与生物多样性威胁因子相结合计算，最终输出生境质量指数和生境退化度图层。

Habitatindex代表生境质量指数，公式如下：

式中，Qxj代表土地利用与土地覆盖j中栅格x的生境质量；Hj代表土地利用与土地覆盖j的生境适合性；k为半饱和常数，当k=0.5时，k值与D值相等；Dxj为景观类型j栅格x的生境胁迫水平，公式如下：

栅格y中胁迫因子r对栅格x中生境的胁迫作用为irxy：

式中，dxy是栅格x与y之间的线性距离；drmax是威胁因子r的最大作用距离，是威胁因子的权重；βx是栅格x的可达性水平，1表示极容易达到；Sjr是景观类型j对威胁因子r的敏感性，该值越接近1表示越敏感。

InVEST模型中假设认为，在一个生态系统中，地类对于威胁因子的敏感性程度越高，则该威胁因子对地类退化程度影响也就越大。生境退化程度由如下公式计算：

生境退化指数＝敏感性分布图层×威胁强度分布图层×权重值

3 结果与分析

3.1 景观格局的时空变化

由图3看出，1999—2017年，北京市浅山区的景观空间格局主要变化表现在以下几方面。1)在土地类型分布方面，农田和建设用地向浅山区内部扩张，在“凤凰岭-阳台山-鹫峰-小西山”“天龙潭-白虎涧”和“泉水河-牤牛河”地区最为明显，在密云水库南部、平谷顺义两区山前地带也多有碎片态建设用地出现。而多处林地退缩突出，整体而言景观破碎化趋势加强。2)在各类土地数量方面，建设用地和未利用地比例增长明显，分别增长1.04%(40.55km2)和4.12%(161.61km2)，反映出浅山地带不合理开发所导致的土地浪费现象较为严重。在5种研究的土地利用类型中，总量以耕地下降为最多，下降了140.65km2。消失的耕地大多转变为了建设用地，也有少量转为林地(表3、4)。

3.2 生境质量变化

生境质量高低代表的是生物多样性高低，为便于判断2个时间点的生境质量变化，根据输出的生境质量取值，并用ArcMap中的自然断裂法(Natural breaks)划分等级。

由图4看出，浅山区内生境质量较低的是农田和建设用地区域，生境质量较高的是林地和水域区域。生境质量变化与土地利用变化类型趋势基本吻合，2017年较1999年，浅山区生境质量分布更碎片化，生境质量较差区域向质量较好区域更深楔入。还有许多整块生境良好区域被低质量区域脉状冲散，这些区域大多为沟域地带。密云水库北岸区域质量变差最为明显。

3.3 生境退化度变化

生境退化程度分值的高低反映了该地类在当前保护程度下，受到人为威胁因子影响程度的大小，也就表示了其潜在的生境破坏与生境质量下降的可能性大小。平原与山区交界处的生境退化度十分突出，几乎连成一线。生境退化度分值最高的是农田及建设用地区域附近(图5)，表明其周边生境具有较高潜在退化的可能性，1999年主要是密云、门头沟区域高程在75～100m的地带退化度较高。2017年较1999年，密云水库、平谷区、房山区附近区域生境退化度明显增强。

3.4 景观生态格局预测与生境质量模拟

预测自然情景下浅山区2035年景观格局(图6)，发现2017—2035年的土地利用覆被变化基本会延续之前趋势，农田和建设用地进一步向浅山腹地延展，侵占部分有林地和疏林地，但速度较之前有明显减缓。之后预测分析了2035年浅山区的生境质量和生境退化度分布(图7)。发现生境质量和退化度中等过渡地区减少，有两极分化趋势，且分布更趋于碎片化，将不利于地区完整生态过程的实现和可持续发展。

3.5 生境稀缺性分析

分析了1999—2017、2017—2035两个时间段内的生境稀缺性(图8)。生境稀缺性得分是通过对其空间位置和影响程度的计算，找出破碎化程度较高，且土地利用类型变化较为频繁的地类斑块或部分栅格，得分值的大小表示应受重视和保护的不同程度。浅山区密云水库附近生境稀缺度较高，说明土地利用变化较频繁且明显，需要加强保护。其他得分值较高的区域大多分布在浅山近平原处的过渡地带，并向浅山内部延伸，表明人为干扰在未来一段时间内仍将是影响浅山地带生态环境的最主要因素，需要防微杜渐，分别因地制宜分析原因，制定保护性发展策略，实现土地精明利用和生态高效。

4 讨论

4.1 研究方法局限性

本研究方法的探究实践，尚存在一些局限性。

1)遥感数据主要来源于Landsat8 OLI夏季数据，遥感图像时相差别及30m空间分辨率对提取结果有影响，导致景观类型解译存在一定误差；解译选取的数据准确度和数量有待提高。

2)基于CA-Markov模型预测景观格局空间变化时，因为采用多地类模拟分析，各地类转换概率的确定需要以几期变化信息为主导，相关参数设定要参考地形、降雨和地区经济发展政策等。故已知数据年份选取宜增多，使模拟参数设定更精确。

3)InVEST模型更多倾向于植被多样性，认为生境质量与生物多样性有很强的正相关关系，但现实中却并非总是如此[35]。且InVEST模型中所有威胁因子均为叠加的，有时多种威胁的集体影响可能远大于个体之和的威胁水平，但模型未能考虑。

4)数据获取限制使得本研究的驱动因子指标体系有待完善。

4.2 研究结论与管理策略

4.2.1 研究结论

北京城市不断扩张会对浅山景观生态系统带来一定的破坏，在追求社会和经济发展的同时不能忽视保护生态环境。应限制盲目泛滥的建设，加强自然地区保育修复。研究发现1999—2017年，浅山区景观空间趋于碎片化，一定程度上影响了生态涵养功能的实现。农田和建设用地向浅山内部扩张，逐渐与林地交织分布。当前山后区整体生境质量优于山前过渡地带。预测2017—2035年仍将沿袭这一趋势，但变化速度有所减缓。

就目前来看，最值得提起注意的是密云水库北部广大地区，在过去有大量林地被人工侵入，生境质量和退化度不容乐观，而且稀缺性很高，即未来恶化可能性较大。还有一些线型的沟域地带中，多处建设用地沿沟谷向山内延展，使原来完整的生境“龟裂”，主要是因为人类开发区楔入原自然背景，过去不合理的经营模式可能造成了人地不和谐因素的深化。但另一方面，研究也发现有一些区域有生境质量好转的现象，推测应该属于生态型保育规划、退耕退工还林还草策略的成果。

4.2.2 对应管理策略

1)整体上：动态格局，持续精明发展。

从土地利用变化中提取规律，维护浅山生态安全和景观连续完整性，可更多采用在山水基质上“嵌入斑块”式发展，而非大面积地“冲撞、侵略”原自然生境。通过自然情景下的土地利用变化情况预测及生态效益评估，明确土地不同的规划目标和开发强度阈值，对于开发土地集约高效利用。

2)当前生境质量不佳区域：分类究因，修复防护并进。

表3 1999—2017年北京市浅山区土地利用类型面积及比例

表4 1999—2017年土地利用转移矩阵

对于当前生境质量表现为较差的区域，需对其深入调研，探究政策、开发和保护活动造成的影响，针对生物和非生物环境修正管理策略，控制不合理开发行为和不得当保护措施。探究其过去和未来的生境质量，如在过去或未来背景下有较好的情况，宜着重进行保育和修复，主要措施包括构建游憩廊道、绿道连接、对裸地增植补植、改造低质低效的林地为混合林等。对于过去和未来都情况不佳的地区，工作重点在于如何避免其不良态势的扩大，主要措施以建设防护隔离带等为主。

3)退化或稀缺度较高区域：规避风险，保障核心安全。

区内存在一些当前生境质量尚可，但有着较高退化度或稀缺度的区域，说明这些地区受到人为威胁因子影响程度大或土地变化频繁，未来恶化风险也较大，应该严加防范，尽力避免“侵蚀”扩大化。具体措施包括杜绝威胁区域安全的行为、设生物围栏，远期升级为核心区，严格保护原有生境和物种等。

注：文中图片均由作者绘制。

图8 1999—2017年(8-1)、2017—2035年(8-2)生境稀缺度分布

[20] 巩杰，马学成，张玲玲，等.基于InVEST模型的甘肃白龙江流域生境质量时空分异[J].水土保持研究，2018，25(3)：191-196.

[21] Deadman P, Brown R D, Gimblett H R. Modelling Rural Residential Settlement Patterns with Cellular Automata[J].Journal of Environmental Management, 1993, 37(2): 147-160.

[22] Batty M, Xie Y. From cells to cities[J].Environment and Planning B: Planning and Design, 1994, 21(7): S31-S48.

[23] White R, Engelen G. Cellular automata as the basis of integrated dynamic regional modelling[J].Environment and Planning B: Planning and Design, 1997, 24(2): 235-246.

[24] Lett C, Silber C, Dube P, et al. Forest dynamics: A spatial gap model simulated on a cellular automata network[J].Canadian Journal of Remote Sensing, 1999, 25(4): 403-411.

[25] Batty M, Xie Y, Sun Z. Modeling urban dynamics through GIS-based cellular automata[J].Computers Environment & Urban Systems, 1999, 23: 205-233.

[26] N o u r q o l i p o u r R, A b d u l R a s h i d B M S, Balasundram S K, et al. A GIS-based model to analyze the spatial and temporal development of oil palm land use in Kuala Langat District, Malaysia[J].Environmental Earth Sciences, 2015, 73: 1687-1700.

[27] 孙贤斌，刘红玉，李玉凤，等.基于CA-Markov模型土地利用对景观格局影响辨识[J].生态与农村环境学报，2009，25(1)：1-7；31.

[28] 赵建军，张洪岩，乔志和，等.基于CA-Markov模型的向海湿地土地覆被变化动态模拟研究[J].自然资源学报，2009，24(12)：2178-2186.

[29] 陈妍，乔飞，江磊.基于InVEST模型的土地利用格局变化对区域尺度生境质量的影响研究：以北京为例[J].北京大学学报：自然科学版，2016，52(3)：553-562.

[30] Hulst R V. On the Dynamics of Vegetation: Markov Chains as Models of Succession[J].Vegetatio, 1979, 40(1): 3-14.

[31] 侯西勇，常斌，于信芳.基于CA-Markov的河西走廊土地利用变化研究[J].农业工程学报，2004(5)：286-291.

[32] 韩玲玲，何政伟，唐菊兴，等.基于CA的城市增长与土地增值动态模拟方法探讨[J].地理与地理信息科学，2003(2)：32-35.

[33] 蔺卿，罗格平，陈曦.LUCC驱动力模型研究综述[J].地理科学进展，2005，24(5)：79-87.

[34] 杨国清，吴志峰，祝国瑞.广州地区土地利用景观格局变化研究[J].农业工程学报，2006，22(5)：218-221.

[35] 吴健生，曹祺文，石淑芹，等.基于土地利用变化的京津冀生境质量时空演变[J].应用生态学报，2015，26(11)：3457-3466.