基于自然恢复的采煤沉陷区植物群落演替模拟系统开发

米家鑫，张绍良，杨永均，侯湖平，丁忠义

(1.中国矿业大学 公共管理学院，江苏 徐州 221116；2.中国矿业大学 环境与测绘学院，江苏 徐州 221116；3.中国矿业大学 教育部矿区生态修复重点实验室,江苏 徐州 221116)

煤炭开发战略西移后，干旱半干旱区高强度地下开采引起了大面积的地表沉陷、水体破坏，严重胁迫着脆弱的生态系统，而大范围的人工修复往往成本高、效率低、效果不明显，因此一些学者提出了自然恢复的策略。尽管一些研究发现采煤沉陷区的植被的确具有自我恢复的能力，但是植被自然恢复过程时间长、效果和格局难以预测，这种不确定性极大地增加了采煤沉陷区植被自然恢复的风险，因此需要开发自然恢复效果评价技术，模拟矿区受损植物群落的自然演替过程成为了重要途径。

植物群落的演替过程十分复杂，涉及到了多个生态水文过程，根据模拟模型的原理主要分为经验模型与机理模型。经验模型主要依赖于植物群落类型、与其生长演替相关要素间的相关关系，在将其拟合为定量的经验函数后进行模拟。例如物种分布模型、Holdridge生命地带模型、BP神经网络模型等。但由于经验模型具有较强的区域性特征导致适用范围较小，在植物群落演替模拟中研究较少。另一方面，机理模型则通过模拟植物群落演替过程中发生的生态水文变化过程，最终输出植物群落的模拟结果。其中应用较为广泛的包括林窗模型和全球植被动态模型。林窗模型(Forest Gap Model)最早由BOTKIN与1972年提出，主要通过描述树木的生长、死亡过程来模拟森林的演替过程，最终模型将会输出模拟过程中不同树木的数量、树木的高度、森林的密度等参数的动态变化。全球植被动态模型(DGVM)，则主要从物质交换和能量流动的角度模拟植物群落的生长演替过程，通过将气象数据、CO浓度数据和土壤数据等输入植被动力模块以及生物地球物理-生物地球化学模块后，输出植物群落的类型和状态。总体来看，已有的植物群落演替模拟系统中主要依靠降雨、光照等气候数据模拟大尺度区域内的物质、能量流动来模拟植物群落的自然演替过程，而地形这类在小范围区域内存在空间差异的因素往往被忽视，没有考虑地形差异对各类生态水文过程的影响。因此这些模拟系统大多难以反映模拟单元内部植物群落在自然恢复过程中的空间分布及实现其可视化，更无法反映因采煤沉陷区地形变化而改变的植物群落演替过程及植被格局。

采煤沉陷区复杂地形条件下的水文生态过程是植被自然恢复模拟的难点。沉陷裂缝出现后改变了原有地形中坡度、坡向、高程等因素，从而影响着地表径流、水分蒸发、水土流失等多个水文过程，进一步影响土壤养分积累、种子萌发、群落发展等过程。因此小尺度下生态水文过程的变化是影响采煤沉陷区植被自然恢复的主要因素。而元胞自动机因其能够通过简单的局部转换规则来模拟十分复杂的空间结构，常被用于生态水文过程模拟中，例如坡面水蚀模型、河道汇流模拟模型、城市内涝模型等。一些学者也将其用于植被格局的模拟中，例如有学者基于元胞自动机模型模拟了高寒草甸退化过程中“秃斑”现象中出现的连通效应；有的结合逻辑斯第增长模型模拟了不同资源条件形成的植被空间格局；还有的构建了林火蔓延快速模型，模拟了火灾发生时树木的燃烧传递现象。然而，已有的基于元胞自动机的植被格局模拟研究大多使用植物群落与单一因素或单个生态水文过程之间的相关关系，作为演替的限制或驱动条件，却没有考虑复杂地形条件下地表径流、水分蒸发、水土流失等多个水文过程相互作用关系以及这一过程中植被的变化，从而难以完整的模拟采煤沉陷区的植被自然恢复过程及最终形成的植被格局。

笔者基于植被自然恢复规律和元胞自动机模型开发了植物群落演替模拟系统，并以山西省大同市云冈区的煤峪口矿采煤沉陷区为例，模拟了植被自然恢复下的群落演替过程，以及演替形成的植被格局，可以为矿区植被自然恢复策略的实施提供案例参考和决策工具。

1 研究区概况

大同市云冈区位于山西省北部，地理位置坐标为东经112.88°～113.19°，北纬40.17°～40.42°，如图1所示。全区地处黄土高原，属大陆性季风气候，全年平均气温6.4 ℃，年平均降水量394.6 mm，主要集中在6—9月份。笔者选择的沉陷裂缝区域位于大同煤业煤峪口矿沉陷裂缝区内，总面积约为0.13 km，区域内主要植被包括沙棘、狗尾草等灌木和草本植物。煤峪口矿年生产能力250万t，自2012年至今，煤峪口煤矿主要对14-2号煤层进行开采，煤层平均厚度为2.11 m，采深在310～376 m，采煤方法以长壁大冒顶机械化综合采煤为主。根据概率积分法的沉陷预计模型，煤峪口煤矿在2006—2018年间出现的下沉影响面积达到8.65 km(下沉大于10 mm)，沉陷达到2 000 mm以上的区域占到8.21%。

图1 研究区位置Fig.1 Location of study area

2 植物群落演替模拟系统设计

2.1 模拟系统总体框架

基于干旱半干旱区气候、水文、植被特征，本研究中的植物群落演替模拟系统由地表形变模块、降雨径流模块和植被生长模块组成，如图2所示。在干旱半干旱矿区，当地表形变发生后，首先以地表沉陷和地表裂缝为代表的地表形变通过根系拉伸等直接引起植被的死亡。此外，地表形变也导致了区域内高程、坡度和坡向的变化，并将对后续的水文过程产生影响。在所有生态过程和要素中，降雨既是地表植被水分的主要来源同时也是植物群落演替的主要限制因素。当降雨到达地表后，土壤入渗随即发生将地表水转化为土壤水；当单位时间的降雨量大于土壤入渗速率时一部分降雨转化为地表径流，并在重力作用下从高向低流出；高程决定了径流的流向，坡度则影响着入渗速率和径流速率，坡向则影响了土壤的蒸发速率。降雨结束后土壤蒸发随即发生，土壤水分逐渐降低至田间持水率。当植被立地条件满足不同类型的植被生长所必须的条件后，植被种子便有概率萌发。而在生长季结束之后，草本植被和一部分灌木植被死亡，一部分灌木植被和乔木植被则凋落，植被残枝和凋落物进入土壤转换为土壤养分；凋落的种子也将进入土壤种子库，从而为次年植被提供养分和种子。这些过程将会不断循环，最终促成了区域植物群落格局的形成。

图2 植物群落演替模拟系统框架Fig.2 Framework of simulation model for spontaneoussuccession of plant communities

本次研究构建的模拟系统主要基于元胞自动机(Cellular Automata，CA)模型，它是一种时间和空间上离散的动力学系统，主要由元胞空间、元胞邻居、元胞状态及转化规则组成。

2.2 模拟系统元胞单元设计

在本研究中，为了兼容地形数据特选择四边形网格的划分方式作为元胞空间(图3(a))。考虑到植物群落的划分尺度，将每个元胞大小设定为现实世界中1 m的方块，下文中植被覆盖度、等状态变量均以1 m大小元胞为最小空间单位。

图3 四边形网格元胞与摩尔型邻居模型Fig.3 Quadrilateral grid cell and Mooreneighbor model

2.3 模拟系统元胞邻居选择

元胞邻居定了元胞之间的传递关系及影响范围。考虑到现实情况下地表径流和植被生长过程中相邻斑块间的作用方式，本次研究选择摩尔型作为此模型的元胞邻居，如图3(b)所示。每个元胞受到与其相邻的8个邻居元胞的作用，在不同的模块中遵循不同的规则变化。

2.4 模拟系统元胞状态参数选取

元胞状态表示元胞在模型中的非空属性集合。每个元胞都有特定的状态值，并随着模型的演化时刻发生改变。本次研究中每个元胞的状态主要包括地形参数、入渗参数、径流参数、蒸发参数、土壤参数、植被参数，以及属于模拟系统的全局参数和降雨参数，见表1。

表1 模拟系统的元胞状态参数

2.5 模拟系统元胞规则设计

转换规则是元胞自动机的核心，其实质是根据元胞当前时刻的状态和邻居状态来决定下一个时刻该元胞的状态的动力学函数。在本次研究中，主要包含降雨规则、入渗规则、径流规则、蒸发规则、养分积累及流失规则、种子库积累规则、植被生长规则。

2.5.1 根系损伤规则

地表形变发生时不仅引起了地表高程、坡度、坡向的变化，同时也导致了部分植被的死亡。考虑到不同类型植被的生长方式及对外界扰动的恢复力差异，设置草本、灌木和乔木植被对沉陷的忍耐阈值分别为1.5，1.0和0.5 m，即当沉陷深度大于1.5 m时，区域内的草本植被即死亡，而灌木和乔木当沉陷超过1.0 m和0.5 m时死亡。在模拟模型中，将地表形变设置为模拟的第1天发生，死亡的植被将作为养分被土壤吸收，养分吸收规律与植被自然死亡时一致。

2.5.2 降雨规则

在真实的降雨场景中，单位时间落入地表的雨水呈现随机分布规律。将雨滴横截面作为单位面积并投影到其落入的地表，将秒作为单位时间并记录每一秒内落入地表的雨水数量，每一秒内落入雨水的地表是随机的，而单位时间内含有降雨落入地表的面积与地表总面积的比值即降雨落入概率。考虑到模型的整体运行速度，将运行的步长即单位时间设定为分钟，空间单位为1 m的元胞。在本次研究中，以云冈区过去50 a的降雨情况作为参考，分别确定不同季节的降雨量。由于云冈区降水70%集中在春夏季，且模拟系统仅将春夏季作为植被的生长期，因此仅将在春夏季时考虑降雨情况。其中设置春季为每年度的第59～151天，夏季为第151～243天，其余时间为秋冬季。具体参数见表2。

表2 生长季降雨参数设置

为了反映区域真实的降雨情况，通过区域的季度降雨总量、降雨天数、降雨时间、降雨强度、雨水落入概率以计算降雨强度，其具体公式为

=

(1)

(2)

式中，为每一步长内元胞获得的雨量;为雨强;为雨滴落入概率;为累计降雨量；为降雨天数；为降雨时间。

2.5.3 土壤入渗规则

土壤入渗是指地表水在重力势、基质势等作用下运移、存储变为土壤水的过程。对于植被覆盖的土壤来说，由于植被覆盖能有效减缓地表水的流动，从而大大增加其稳定入渗速率；坡度的增加则会减少来自重力的分量，从而引起入渗速率的减小，如图4所示。而对于无植被覆盖的土壤来说，由于其地表更容易受到降雨径流的侵蚀从而引起土壤质地的恶化，导致其土壤入渗速率的降低。因此针对裸地土壤增加了一个养分的调节参数，当土壤养分低于土壤初始条件时，其土壤入渗率将逐渐下降，以模拟现实情况下土壤侵蚀增加并逐渐沙化的现象。

图4 降水在坡面的重力加速度分量Fig.4 Gravity acceleration component of precipitation on slope

=(191+0326e-0215)cos

(3)

(4)

式中，，分别为植被覆盖和无植被覆盖区域的土壤入渗速率；为入渗时间；为土壤养分含量；为土壤初始养分含量。

为了保证模型运行的连贯性，在此模型中仅考虑表土土壤水分的变化。入渗结束后，表土土壤的水分增加量即等于累计入渗量，而表土为地表20 cm以内的土壤。当土壤水分达到土壤的饱和含水量时，土壤入渗过程将会暂停，多余的地表水会以径流的形式流向周边区域。通过对研究区域土壤的采样化验，测得植被覆盖和无植被覆盖土壤的饱和含水率分别0.35与0.30。为了便于降雨量、入渗量、含水量等不同量纲变量之间的转化，笔者定义的含水量均统一为单位体积内土壤水分的比例。由于降雨形成的地表水以及地表土壤的投影面积均为1 m，所以土壤含水量等于累计入渗量的深度与表土深度(20 cm)的比值。

2.5.4 径流规则

当地表水大于土壤入渗速率时开始出现径流。对于单个元胞来说，地表水输入来源包括降雨和来自其他元胞的径流，输出方式则包括土壤入渗及流出至其他元胞的径流。本文中降雨量的单位为mm，其含义为单位时间内落入单位面积内水量的深度。但对于地表径流，由于坡度的存在使得地表径流实际通过的面积将大于元胞的投影面积，如图5所示。因此，从降雨中获取的径流量需要进行转化，根据勾股定理可得，坡度为的坡面实际单位面积为1/cos，则径流公式为

′=

(5)

=′cos

(6)

式中，′为单位元胞通过降雨获得的理论径流量；为考虑地表坡度后单位元胞上的实际径流量。

图5 不同坡度元胞的径流状态Fig.5 Runoff status of cells with different slopes

当径流在一个元胞(,)内形成后，如果元胞(,)并非与其相邻的8个邻居元胞中的高程最小元胞，那么径流将在重力的作用下从该元胞流向相对于此元胞的最低点元胞(+1,+1)中。由于重力方向是固定的，所有每个元胞形成径流的流出方向只有1个，径流将始终流入高程相对较低的元胞，如图6所示。

图6 非相对最低点元胞径流的流出方向Fig.6 Outflow direction of the runoff from non relativelowest cellular

当元胞(,)的高程小于周围所有邻居元胞时，此元胞则成为径流的汇聚点，邻居元胞形成的径流都将流入元胞(,)中，且元胞(,)形成的径流则不会流出。但是，当元胞(,)内汇集的径流超出其所能汇集的最大值时，元胞(,)的邻居元胞则将停止流出。元胞(,)的邻居元胞内的径流由于无法再流出，这些邻居元胞也将逐渐汇集径流，从而出现径流由低向高增加的“溢出”现象，如图7所示(图中为汇聚点元胞的高程；为汇集点元胞周边的最高点元胞的高程；Drop为汇聚点元胞与最高点元胞之间的落差)。

图7 相对最低点元胞径流的流出方向Fig.7 Outflow direction of the runoff from relative lowest cellular

在径流过程中，元胞(,)的流出量与其径流量以及径流速率有关，而流入量则等于其邻居元胞的流向元胞(,)的流出量总和。当元胞内的径流量小于元胞的径流速率时，流出量等于此时的径流量当元胞内的径流量大于元胞的径流速率时，流出量则等于元胞的径流速率。其中，元胞(,)的径流速率等于单位时间内通过元胞的最大径流量，这一速率通过曼宁公式计算得到。曼宁公式是一个估测液体在开放管道(即明渠流)或非满管流中平均速度的经验公式，具体公式为

(7)

式中,为转换系数，国际标准制为1；为曼宁粗糙系数，通过查询曼宁系数表确定；为水力半径，计算方式是过水截面积与湿周长的比值；为水力坡线或是线性扬程损失的斜率。

在本次研究中，将每个元胞作为径流流动的基本空间单元，则为通过元胞的的水深，为了便于降雨量、径流量之间的水量交互与计算，对径流量与降雨量均统一使用毫米作为单位，因此在计算最大径流速率的公式中将径流量直接用于水深计算；经过查询曼妮系数表后，规定裸土的粗糙系数为0.02，植被覆盖土壤的粗糙系数为0.095；斜率则等于坡度的正切值tan。

2.5.5 蒸发规则

当土壤含水率高于土壤的田间持水率时则土壤水分开始蒸发，这一过程主要集中在无降雨期间。土壤的田间持水率主要受到覆盖类型的影响，存在植被覆盖的土壤其保水持水能力较强，而裸土则相对较弱。本研究中，测得植被覆盖和无植被覆盖土壤的田间持水率分别为0.22与0.15。蒸发速率在不同季节内有地形位置的不同而有所差异。由于不同坡向的差异，阳坡区域的土壤由于光照时间较长而温度较高，土壤的蒸发速率也较高。地表的植被覆盖则会一定程度的降低地表土壤的温度，从而降低土壤的蒸发速率。参考云冈区过去50 a的气象数据，并将蒸发速率转化为日蒸发量后，确定研究区不同季节、不同坡位及不同植被覆盖土壤的蒸发速率见表3。

表3 不同坡向的角度范围及蒸发速率

2.5.6 土壤养分积累及流失规则

土壤养分的积累主要来源于植被凋落物和残枝分解。当植被在秋冬季凋落枝叶时，其凋落物将被土壤中的微生物分解并转换为土壤养分；当植被死亡时，其残留的根茎也将被土壤吸收转化土壤养分其中，灌木和乔木植被包含植被死亡和凋落过程，而草本植被只通过死亡过程转化养分。此外，由于灌木植被和乔木植被的冠幅较大，其凋落和死亡后的枝叶通常也将落入周边区域，因此乔木和灌木植被的邻居元胞也将在其凋落和死亡时获得养分。本文使用土壤有机质含量反映土壤养分情况，基于研究区土壤养分水平情况发育规律，规定不同植被凋落或死亡后土壤养分的增加情况见表4。为了避免局部区域养分过高不符合现实情况，同时设置土壤养分的最大值为40 g/kg。

表4 植被死亡或凋落时土壤养分增加量

另一方面，径流对土壤的冲刷裹挟作用导致了土壤养分的改变。径流通过土壤时，一部分土壤以及土壤中的有机质、氮、磷、钾等养分会在径流的裹挟作用下融入到径流中，从而造成土壤养分的流失。径流的强度越大，其造成的养分流失则越明显；降雨的强度越大，土壤养分的速率越快。在这一过程中，地表植被则会起到减弱径流、稳固水土的作用，因此植被覆盖土壤的养分流失速率将明显低于裸土的流失速率。以往的研究针对径流引起的土壤养分流失积累大量预测模型，通过参考之前的研究，确定了本研究中不同类型土壤的养分流失速率(表5)。此外，分别设置5，30 g/kg作为土壤养分的最低值和最高值，为避免一些区域土壤养分持续流出或流入而造成养分过高或过低的情况。

表5 土壤养分流失速率

2.5.7 种子库积累规则

植被在生长过程中都会产生一定数量的种子，这些种子大多会凋落在植被周围并进入土壤种子库中；相似的，周边区域植被的种子也可能掉入其相邻区域，从而满足植被生长所需的土壤种子库条件。此外，风力运输等偶然事件同样会为某一区域的土壤提供种子。在本次研究中，规定植被死亡后其种子数量增加1；而一天内单个元胞因为偶然事件获得种子的概率为0.000 4，偶然事件发生后其种子数量增加1。对于来自相邻植被获得种子，则统一使用逻辑判断，即当周边存在相同类型的植被时，则满足植被生长的种子库条件。

2.5.8 植被生长规则

不同植被生长所需的土壤水分、养分、种子库以及种子萌发概率存在一定差异。植物群落的发展基本遵循草本—灌木—乔木的演替顺序，且灌木植被所需的土壤养分条件一般高于草本植被，但在水分条件上无明显差异，而乔木所需的土壤养分和水分条件均较高。此外为了模拟植物群落的“聚集”现象，同时考虑邻居元胞植被覆盖的影响。基于以上规则，在参考以往黄土丘陵区植被自然恢复的研究案例基础上，确定了不同类型植被所需的土壤养分、水分和种子库条件见表6。

表6 不同类型植被种子萌发条件

而对于不同类型植被的萌发概率来说，在一定范围内水分条件和养分条件越充足，其种子萌发的概率越高；邻居元胞中存在灌木和乔木的数量越多，种子萌发的概率同样越高。因此通过土壤实际的水分和养分与种子萌发所需的水分和养分的比值作为调节参数，以反映不同水分和养分的土壤在种子萌发条件上的差异。

2.6 模拟区域环境条件设置

此次研究中使用研究区已获取的无人机影像作为模拟区域的起始地形条件，并将其植被格局作为模拟结果的验证数据。研究区位于云冈矿区内煤峪口矿开采范围内，开采时间在2006—2008年。无人机影像采集时间为2019-08-22 T 09:00—12:00，距沉陷裂缝形成时间约为11 a，影像获取时区域内植被生长茂盛。无人机飞行高度为100 m，云台角度-90°，旁向重叠率65%，航向重叠率70%。获取无人机影像后，使用Agisoft PhotoScan软件进行影像预处理，并通过数字表面模型去除植被高度后获得数字地形模型。然后将数字地形模型影像作为模拟区域的高程数据，并在ArcMap中通过栅格分析功能转化为区域的坡度及坡向数据，如图8所示。模拟区域的初始植被条件设置为无植被，而验证数据则根据通过数字正射影像确定模拟区域实际的植被类型，如图8所示。为了与植物群落演替模拟系统中的元胞空间匹配，在ArcMap中将高程、坡度、坡向及植被类型数据重采样为1 m×1 m的栅格。考虑到模拟区域内存在少量人为建造的道路，在ArcMap中通过矢量化构建了掩膜图层并在模拟结束后添加至模拟结果中。

图8 沉陷裂缝区地形条件Fig.8 Terrains of subsidence area

2.7 模拟结果精度评价方法

由于植物群落自然演替的随机性，比较不同景观尺度下的植被格局进行精度评价较为合理。首先将模拟结果与验证数据划分为若干个大小的斑块(图9)，依次比较各个斑块内不同类型的植被覆盖度。当模拟结果与验证数据中对应斑块内植被覆盖度、草本覆盖度及灌木覆盖度相似，则说明此模拟系统能较好的模拟现实情况中的群落演替。考虑到不同景观尺度下植物群落格局的差异，本次研究中选择了10×10，20×20,50×50,80×80,100×100 m的5种斑块大小(P10，P20，P50，P80，P100)。

图9 植被覆盖类型及群落斑块划分情况Fig.9 Types and patches of plant communities

在本研究中，采用归一化均方误差反映模型的模拟精度。归一化均方误差是一种反映2张图像之间相似度的指数，其值域为0到正无穷，接近0则说明实际值与预测值之间差异越小，模型的模拟精度越高，具体公式为

(8)

当模拟结果中的植被覆盖度接近验证数据的植被覆盖度时，即停止模拟并将此时的植物群落格局记录并与验证数据进行比较。在实际情况中，植被覆盖度为50.79%，草本覆盖度为43.27%，灌木覆盖度为7.51%，而乔木覆盖度仅为0.01%。因此在进行模拟时仅将草本和灌木植被作为模拟结束和精度评价的标准，当模拟结果的草本植被覆盖度为35%～55%。且灌木覆盖度为7%～8%时即停止模拟，输出此时的植被格局并分别计算不同尺度下的。为了减少单次模拟带来的随机性，设置模拟次数为50次，取每次模拟结果与验证数据的NMSE平均值为最终值。

3 结果与讨论

3.1 模拟结果与系统精度

植物群落的模拟结果与实际分布对比如图10所示，其中图10(a)～(c)为随机选择的植被自然恢复模拟结果，图10(d)为实际的植被分布情况。从模拟结果与实际结果的植物群落整体来看，灌木植物群落均主要集中在沉陷裂缝区的低矮区域，且都表现出较为明显的群落聚集现象。但在模拟结果中，草本植被在整个区域内较为均匀分布，无明显的无植被覆盖区域。而实际情况中存在较为明显的无植被覆盖区域。实际调查发现，这些区域均为陡峭的边坡，经过长年的雨水冲刷以及成为裸露的岩石。这一差异主要由于此次模拟的初始土壤被设置为适宜草本植被的条件，但在实际情况下这些区域可能在自然恢复开始时便不适宜植被的生长，从而导致了模拟结果与实际情况在这些局部区域的差异。

图10 验证区域植物群落格局模拟结果Fig.10 Simulation results of vegetation pattern in stuyd area

从模拟结果与实际格局的来看，随着斑块尺度的增加模拟结果与实际情况的差异度随之降低，表明模拟精度也随之增加。总体模拟精度如图11所示，其中，灌木植被的空间分布格局差异度最高，5种斑块尺度下的在0.06～0.59；草本植被的总体较小，在0.05～0.14；而从总体的植被覆盖度来看，植被格局的模拟结果与实际情况差异最小，其在0.04～0.12。这主要来自于不同的植被覆盖度基数的差异，其中灌木覆盖度基数最小，因此在划分的板块中较小的差异都将被放大，尤其考虑到植被的自然演替本身便是充满随机性的过程，使得灌木植物群落格局的模拟在较小尺度精度较低。但从整体效果来看，目前设计的模拟系统能够还原现实世界中的植物群落格局。

图11 不同尺度下模拟结果和实际情况的差异度Fig.11 NMSE between simulated and actual vegetationpatterns at different scales

3.2 植物群落演替的时间变化模拟结果

沉陷裂缝区的植物群落在自然恢复初期以草本植被为主，在自然恢复约3 a后开始出现灌木植被。图12反映了沉陷裂缝区内不同类型植被覆盖度在自然恢复10 a间的动态变化，图13展示了植被自然恢复模拟过程中年度生长期最后一天的植被覆盖情况。在自然恢复初期，草本植物群落的发展速度较快，灌木植被由于对土壤养分和种子数量的要求较高，在自然恢复的前3 a间无灌木植被出现。但随着草本植被凋落死亡并转变为土壤养分，越来越多的区域内满足灌木的生长条件，因此在自然恢复5 a后灌木植被表现出明显的增长趋势。在本次研究中，群落演替模拟的停止标志为植被覆盖度介于35%～55%，灌木覆盖度介于7%～8%，50次模拟结果表明到达这一植被覆盖条件评价平均需要10.46 a，这与以往研究中的植被自然恢复时间相似。在以往的自然恢复实例中发现，植物群落的自然恢复基本遵循草本群落→草本-灌木群落→草本-灌木-乔木群落的演替顺序，其中灌木植被一般在自然恢复3～5 a后出现，而乔木植被的恢复需要十年乃至几十年。本次研究的50次模拟结果中，仅2次模拟结果在自然恢复10 a后开始出现乔木植被，这与实际的乔木植被自然恢复规律基本一致。

图12 模拟自然恢复过程中植被覆盖度的时序变化Fig.12 Dynamic of vegetation coverage during succession

图13 植被自然恢复过程中每一年度生长期末的植被覆盖情况Fig.13 Vegetation coverage at the last day of growth peroid in each year duing spontaneous succession

3.3 植物群落演替的空间分布模拟结果

沉陷裂缝区内低矮的沟壑区域由于水分与养分的汇集，成为了灌木植被的聚集区域。通过对比植被自然恢复模拟结果和沉陷裂缝区的地形条件可以发现，灌木植被主要集中分布在高程较低、坡度较小且坡向为北向的区域。由于本次研究中设置的初始土壤养分条件仅满足草本植被的生长，因此在植被自然恢复初期不存在灌木植被。但随着土壤获取草本植被凋落物而积累养分的过程中，地表径流通过裹挟区域内的土壤养分，在重力作用下不断向低矮区域汇集使其成为了养分和水分的富集区，因此高程较低的区域往往率先满足了灌木植被的生长条件。而坡度较小的区域则由于土壤入渗速率、径流速率较小，相较于坡度较大的土壤也够储存更多的水分和养分。坡向则主要影响着土壤的蒸发速率，这使得面向北方的阴坡更适宜植被的生长。因此不难发现，干旱、半干旱矿区沉陷裂缝的出现，在一定程度上反而促进了灌木植被的生长。沉陷裂缝一方面导致了部分区域成为了径流和养分汇集的富集区，同时也增加了周边区域的坡度，提高了径流速率而降低了土壤入渗速率，加剧了干旱半干旱矿区的水土流失现象，逐渐形成了植物群落集中于低矮沟壑区域的空间分布格局。

4 结 论

(1)基于开发的群落演替模拟系统模拟了大同市云冈矿区一处沉陷裂缝区的自然恢复后形成的植被格局，模拟结果与实际植被格局间的误差在0.04～0.59，其中灌木植被格局误差在0.06～0.59，草本植被格局误差在0.05～0.14，整体植被格局误差在0.04～0.12，模拟结果与实际的植被空间分布规律一致。

(2)以裸地为初始条件的50次模拟结果表明，研究区通过自然恢复使草本覆盖度达到35%～55%且灌木覆盖度达到7%～8%水平时，平均需要10.46 a，模拟结果与实际的植被动态变化规律一致。

(3)模拟结果中草本植被在自然恢复前期快速发展，灌木植被在恢复3 a后开始出现，且在5 a后增长加速，但乔木植被在50次模拟中仅有2次在自然恢复10 a后出现，模拟结果与实际的植物群落演替规律一致。