基于遥感与系统动力学模型的土地利用/覆被变化研究

钱国英，张长勤，董斌，李鑫，李欣阳

(安徽农业大学理学院，合肥 230036)

1 引 言

系统动力学模型应用于土地利用/覆被变化(Land-Use and Land-Cover Change，LUCC)的研究一直倍受学者关注。当前，国内外有学者利用元胞机和系统动力学模型构建了土地利用变化动力学模型LUSD(Land Use Scenarios Dynamics model)；有学者利用土地利用/覆被变化变化动态模拟CLUE_S和Markov模型模拟的数据进行分析比较；有利用气候、经济的变化和人类的行为对LUCC模型进行研究[1-5]。这些研究不少尚处于概念验证阶段和多个系统的理论比较阶段，仅单独利用理论数据和实际数据进行“静态陈述”阶段，很难从模拟和实际结果中得到更多的可用信息。本文以安徽升金湖国家级自然保护区为研究对象，利用ArcGIS、ERDAS和STELLA软件构建土地利用变化系统动力学模型，在分析研究区土地利用变化的基础上，利用典型的景观指数分析了土地利用景观格局，并构建系统动力学模型，对研究区的土地利用做出仿真和预测，利用仿真后的数据进行模型分析和验证，最后，分析了研究区土地利用的变化趋势及特征，为研究区乃至湿地自然保护区的土地科学利用提供理论依据。

2 研究区概况

安徽升金湖国家自然保护区(以下简称升金湖保护区或保护区)位于长江中下游南岸安徽省池州境内，坐落于东至县与贵池区相交处，地理位置东经116°55′至东经117°15′，北纬30°15′至北纬30°30′。本文研究区的面积为33329hm2，属亚热带季风气候，平均年降雨量1600mm，平均气温16.14℃，土壤主要以潮土和水稻土为主，保护区人口约2万余人，现有的水面平均面积13300hm2，荒滩荒地为1000hm2，农田为3400hm2，林地为4000hm2，其他11629hm2。研究区属内陆湿地与水域生态系统类型自然保护区。

3 研究方法

研究人类社会与自然景观空间之间的相互关系是景观生态学研究的重要内容之一[6]。在景观生态学中景观是一个具有时间和空间异质性的区域[7]。针对研究区域地理位置和研究的对象，采用了1986年、2004年和2011年的TM影像作为分析数据，以及1∶1万的地形图和1∶1万的SPOT的遥感数据参考数据，以上3期TM影像数据都是6月～8月份时期的数据，符合分析在相同时期下的要求。本文用ERDAS 9.2分别将1986年、2004年、2011年的TM影像的5、4、3波段合成假彩色影像，利用手持GPS实地勘测，然后利用ERDAS软件将研究区域裁剪下来，对遥感影像进行监督分类，并对各土地利用类型进行面积计算。应用Fragstats3.3软件对土地利用分类后的结果进行景观格局分析，计算各土地利用景观指数，利用STELLA软件对土地利用变化进行仿真和预测。图1为研究技术路线图。

图1 研究技术路线

3.1 土地利用景观格局分析

3.1.1 土地利用变化分析

研究区域的形成可追溯到地质史上300万年以前，经喜马拉雅运动后，长江沿岸抬升形成河流，河流汇聚发育形成湖泊，湖泊的平均面积为13300hm2。由于研究的区域主要以农田、渔业和林业，所以根据当地的实际情况把研究区地物类型分为建设用地、林地、耕地、水域、未利用地、草地6种类型[8]，分别对每年的分类结果进行精度评价，结果见表1。从表所示的总体Kappa系数都大于85%，符合研究的要求[9]。在选择训练样本时，使每一地类的训练样本中像元的地物类型与研究区的实际地物类型相一致。另外，选择的训练样本要具有相当高的代表性，能够精确地代表所选地物。根据Kappa系数的结果和训练样地的选择可确定本文的分类结果是真实可用的。根据分类结果初步得出3个时期的各地类的面积及所占比例如表2、图2所示。从表2和图2中可以看出，研究的区域中耕地和水域占主导地位。从表2和图3中可以看出建设用地由原来1986年的5.84%变化成2004年2.87%，到2011年18.34%，草地、耕地、林地、水域和未利用地都有不同幅度的变化。从表2可以看出2004年草地、林地、未利用地都有不同的减少，水域有大面积的增加，从中可判断水域是导致建设用地、林地和未利用地减少的主要原因，查阅2004年同遥感影像相同时期的池州水利史志年鉴。结果表明，该地区发生了洪涝灾害，该研究时期建设用地的减少是由洪涝引起的。

表1 不同年份的Kappa系数

表2 研究区域土地利用类型面积及所占比例

图2 不同地物类型所占面积的变化

从表2和图2可知2004年耕地占了48.94%，所占的比例比其他任何的地类的比例都要多，不单是2004年的耕地比例最大，在其他年份的农田所占比例也比其他土地类型比例大，因为我国是主要的农业大国，以发展农业为主，这与国家的政策是分不开的。

近20多年来升金湖出现过多次的洪涝，其中1997年～1998年发生了一次比较大型的洪涝，2004年也受到水涝的影响。所以在图3中可以明显看出2004年的水域所占比例比其他年份水域比例都大，达到了41.53%。由于未利用地大都是滩涂地，所以当2004年研究区域受灾后，未利用地也随之下降。

研究区域的地理位置决定了水域和耕地占重要的位置。未利用在逐步被开发，林地在21世纪初期有所上升，这是受政府退耕还林政策的影响，其中建设用地发展最快。

3.1.2 土地利用景观格局分析

本文以定量化的指标进行景观结构的描述和分析。本文采用景观格局软件Fragstats3.3进行分析，依据研究区的情况和相关的研究选择了斑块数(Number of Patches，NP)、斑块密度(Patch Density，PD)、总面积(Total Class Area，CA)、景观形状指数(Landscape Shape Index，LSI)、最大斑块指数(Largest Patch Index，LPI)、香农多样性指数(Shannon’s Diversity Index，SHDI)、香农均匀度指数(Shannon’s Evenness Index，SHEI)、边界密度(Edge Density，ED)、破碎化指数(Fractal Dimension Index，FDI)等进行分析[10-11]。

本文利用Fragstats3.3对升金湖保护区土地利用景观指数进行计算，结果如表3所示。从表3中可知1986年的破碎化程度比2011年要低，但比2004年要高。其中NP的变大可以改变物种间相互作用和协同共生的稳定性，其中PD的效果和NP的效果是一致的，NP变大则PD也是会变大它们成正比关系；从表3中可知LPI在2004年最大，其中LPI的大小决定了景观中的优势类型等。其值的大小决定着景观中内部的丰度等生态特征；ED的大小可以从侧面反应景观丰富度的生态因子，当ED的数值变大则其中的斑块数目会比较多，它反应的景观丰富度就会比较高。LSI是指形状指数的生态因子，但LSI的数值越大可以说明斑块的边缘发展越没有规则，但与其他类型的连通性会比较好。SHDI和SHEI所反应的效果比较接近，其中比较2004年和1986年的SHDI，2011年的SHDI数值、LPI和PDI是最大的，这表明2011年升金湖保护区的景观系统中，土地利用比之前的丰富，破碎化程度也比之前的要高。

表3 研究区主要景观指数变化状况

3.2 土地利用景观系统动力学模型仿真

国内外的系统动力学模型一直倍受关注，有在土地利用结构上优化的研究，有在天敌-稻纵卷叶螟系统上的研究[12-13]，这些研究在土地结构上的优化给决策部门提供了一定的建议和评价。本文则是利用STELLA软件通过现有的土地景观仿真2004年在未受到洪涝灾害的影响会是怎样的土地景观格局，分析该区受灾害的原因，并对以后的景观格局做一定的预测和建议。

3.2.1 土地利用变化系统动力学模型构建

本文利用STELLA软件进行建模，将6个地类视为该区域土地利用系统中的6个库(stock)，它们之间的相互转化视为流，影响状态变量变化的因素即为流速变量，(即本研究中的状态转移概率)，1986年～2011年的转移概率如表4所示，各土地利用类型间的相互转化构成系统内部的因果反馈关系。

表4 1986年～2011年土地利用状态转移概率矩阵(%)

为了定量地分析各土地类型间(库)的相互转换关系，利用ArcGIS中的栅格计算功能计算出1986年～2011年的转移面积，利用转移面积计算出这两个年代中6个地类的转移概率。

3.2.2 模型验证

系统动力学模型主要通过理论校验和历史仿真校验，其中理论校验主要是研究所建立的模型是否符合实际要求，模型变量的关系是否有意义以及所列的方程是否符合能表达变量。历史仿真校正选定某一时刻，将仿真得到的结果与实际结果进行对照，检查这两者是否吻合，验证模型是否能有效地代表实际结果[13-14]。

秦钟等[15]在基于系统动力学的土地利用变化研究中对STELLA的模拟数据和现实数据进行比较，结果显示STELLA模拟软件对土地利用类型转化的预测具有较高精度，模拟结果比较可信，并且把预测的结果和CA-Markov模型预测结果比较，相对误差也不到0.9%，说明STELLA软件预测的结果是确切可信。李志等[16]在基于系统动力学城市边缘化土地利用变化模拟与预测中也对系统动力模型进行了模型的验证，结果表明模型的仿真结果最大的误差为5%，符合模型的有效性。

本文以1986年升金湖地区土地利用格局为初始状态，以1986年升金湖地区土地利用格局为初始值，向经过验证的STELLA模型输入数据，进行系统动力学仿真和预测。仿真的结果见表5，将表5中2011年的仿真结果和由遥感观测的结果(表6)对比从中可知仿真结果和实际值的之间相差很少，变化幅度在6%以下，符合历史仿真检验范围，所以本文的系统动力学模型符合理论和历史的检验。

3.2.3 模型稳定性评价

系统动力学模型的稳定性的判定是通过模型的仿真结果来判断的，根据表6将2011年的仿真结果与实际结果比较可以得知，变化幅度在6%以下，符合模型精度范围要求，即所建的模型比较稳定[17]。另外，本文将表5中每年的各地类的面积相加，结果与研究区的实际面积相符。利用表6和表5的数据进行模型稳定性检验，结果表明，该模型的稳定性较高。

4 结果与分析

根据土地利用转移概率(表4)，利用系统动力学模拟软件STELLA进行仿真，结果见图3，运行时间段为1986年～2011年，利用Euler法求解微分方程。因为状态变量的变化每经过25即为一步，设置模拟步长为0.01年，故当n=1时，预测年份为2011年。仿真的年份是2004年，则状态变量变化要经过0.68步会与仿真的年份2004年是一样的。模型运行的结果见表5，将表5中2011年的仿真结果与表2中2011年的实际值进行比较得知，其中建设用地和林地的变化是最大的，而农田、草地、水域、未利用地变化较小。从中得知STELLA预测对农田、林地、水域和未利用地是最有效的。从总的结果来看，模型的模拟结果与实际比较相符，可以做进一步的分析预测和仿真。

图3 系统动力学模型

模型的动力方程式如下：

建设用地：js(t)=js(t-dt)+(f1+e1+d1+c1+b1-a5-a4-a3-a2-a1)*dt

草地：cd(t)=cd(t-dt)+(a1+f2+e2+d2+c2-b1-b5-b4-b3-b2)*dt

耕地：nt(t)=nt(t-dt)+(a2+f3+e3+d3+b2-c1-c5-c4-c3-c2)*dt

林地：sl(t)=sl(t-dt)+(a3+f4+e4+b3+c3-d1-d5-d4-d2-d3)*dt

水域：sy(t)=sy(t-dt)+(a4+f5+b4+c4+d4-e1-e5-e2-e3-e4)*dt

未利用地：wd(t)=wd(t-dt)+(a5+e5+d5+c5+b5-f1-f5-f4-f3-f2)*dt

(1)

上述方程中的参数如a1，a2，……表示建设用地转换为其他地类面积，即是模型中的“流”。建设用地转换成草地，耕地，林地水域，未利用地的面积分别用a1，a2，a3，a4，a5表示，t表示的是时间。利用表4中的转移概率直接就可以求出它们之间的转换面积。以建设用地转移其他用地为例：

表5 系统动力学(STELLA)模型模拟结果

a5=js*pwd1，a4=js*psy1，a3=js*psl1，a2=js*pnt1，a1=js*pcd1

(2)

上式中的pwd1，psy1，psl1，pnt1，pcd1分别表示在一定时期内由建设用地转移为未利用地，水域，耕地，林地，草地的概率。

本文运用遥感分类得到的土地利用类型结果利用STELLA进行仿真和预测，由于研究的区域经常会受到不确定因素的影响，如1998年和2003年～2004年研究区域受到了很大的洪涝灾害，而遥感影像有时效性，所以遥感影像在特殊的情况下不能准确的表达这一时期的土地利用状况，从而会给决策部门在土地利用决策中带来不准确的判断，为了避免出现上述的情况，利用系统仿真是非常必要的，另外利用仿真还可以推断出土地利用类型受到灾害的情况。

本文利用系统动力学模型对2004年的升金湖国家自然保护区的土地利用情况进行了仿真(表5)。将2004年的仿真结果和实际土地利用结果进行比较(表6)。从表中可知2004年受洪涝灾害最严重的是未利用地，其次是建设用地、草地、林地、水域，其中耕地受到的影响最小。从表6中可以看出研究区域2004年有77%的建设用地受到灾害影响，而草地、林地、耕地也在不同程度上受到影响，其中在2004年洪涝中升金湖的湖面面积增加47%，所增加的面积为4441.07hm2。在表6中把遥感分类的结果和仿真的结果进行比较可以很容易得出受灾害最严重的土地类型。因此将土地利用变化遥感分类结果和STELLA系统动力预测结果对比能够得出更多的可用信息。利用STELLA对每一年土地利用类型进行预测，从表5中知，可本文预测了2036年的土地利用情况，到2036年建设用地比2011年的建设用地上升了6.45%，逐渐趋于平稳状态，草地下降7.25%，林地上升了30.53%，这符合国家退耕还林的政策，水域趋于平稳下降，而未利用则得到快速地开发利用，到2036年未利用地下降了35.84%，这些预测的数据对以后土地利用类型的发展方向起到了参照和正确的指引作用。

表6 2004年和2011年实际土地利用面积与仿真结果比较

5 结束语

本文在分析土地利用景观格局过程中发现，在一定的时期它无法正确地描述土地变化之间的转换关系。本文利用同一时期两个年份能够真实反映当年土地利用情况的遥感数据，计算出它们内部的转换概率；再利用这些土地类型内部的转换关系建立系统转换模型；利用系统动力学软件STELLA将土地类型的面积和转换率对应输入系统动力学软件的各个“库”和“流”中，运行软件将不同时期的土地类型进行仿真和预测。本文利用STELLA对2004年升金湖的土地利用系统进行了仿真，结果表明STELLA模拟具有很好的开放性和灵活性，能够有效地仿真土地利用变化的格局，在土地利用规划方面具有良好的应用前景，利用STELLA的预测作用对升金湖的发展趋势做一个指导作用。

利用对土地利用变化的遥感观测数据、景观生态学的分析数据和STELLA的模拟预测数据对有关的土地利用变化进行研究还是不够的，应该综合社会、经济、文化等多种影响要素，构建土地利用驱动机制定量分析模型，通过分析土地利用的时空变化以及与导致这种变化的驱动力之间的因果反馈关系，才能充分揭示人地系统相互作用的内在机制，并以此来仿真和预测未来发展趋势。由于实际系统的复杂性，不可能每个因素都能够考虑在内，系统模型只能不断地完善，STELLA软件与GIS的融合应用还有待进一步的深入研究。希望系统动力学能够有效地应用在土地利用规划和其他各行各业，让它能够有效地仿真和预测，为人类做出正确的决策提供依据。

参考文献：

[1] 何春阳，史培军，陈晋，等.基于系统动力学模型和元胞自动机模型的土地利用情景模型研究[J].中国科学(D辑：地球科学)，2005，35 (5)：464-473.

[2] 邓祥征，刘纪远，战金艳，等.区域土地利用变化的多情景分析—以内蒙古自治区太仆寺旗为例[J].地球信息科学，2004，6(1)：81-87.

[3] 唐华俊，吴文斌，杨鹏，等.土地利用/土地覆被变化(LUCC)模型研究进展[J].地理学报，2009，64(4)：456-468.

[4] DARWIN R F.A global framework for integrated land use/cover modeling[C].GCTE-LUCC Open Science Conference，1998，14(18)：1-18.

[5] ZHOU Z.Landscape changes in a rural area in China[J].Landscape and Urban Planning，2000，47(1)：33-38.

[6] FORMAN RTT，OODRON M.Landscape ecology[M].NewYork：New York Willey，1986.

[7] PICKETT S T A，CADANASSO M L.景观生态学：生态系统的空间异质性[J].生态学杂志,1999，18(6)B：71-74.

[8] 冯春，陈建平.土地利用/土地覆盖研究中遥感图像分类精度的提高方法[J].地理与地理信息科学，2003，19(3)：26-28.

[9] 董斌，冯仲科，钱国英.基于遥感的山东黄河流域土地利用动态及有林地转移[J].福建农林大学学报(自然科学版)，2010，39(6)：646-650.

[10] 傅伯杰，陈利顶，马克明，等.景观生态学原理及应用[M].北京：科学出版社,2001.

[11] 孔维财，王让会，吴明辉.墨玉县绿洲人工林景观格局分析及其生态效益[J].遥感技术与应用，2011，26(1)：89-95.

[12] 秦钟，章家恩，骆世明，等.基于STELLA模型的天敌-稻纵卷叶螟系统的情景分析[J].中山大学学报，2012，51(3)：90-97.

[13] 江景波，华楠.城市土地利用总体规划——方法、模型、应用[M].上海同济大学出版社，1997：80-88.

[14] 张增峰，黄克龙.城乡结合部土地利用问题及对策探讨[J].南京师大学报(自然科学版)，2001，24(2)：110-118.

[15] 秦钟，章家恩，骆世明，等.基于系统动力学的土地利用变化研究[J].华南农业大学学报，2009，30(1)：89-93.

[16] 李志，周生路，卢长林，等.基于系统动力学城市边缘区土地利用变化模拟与预测——以南京市江宁区为例[J].土壤，2010，42(2)：314-318.

[17] 罗昌，贾素玲，王惠文.基于系统动力学的供应链稳定性研究[J].系统仿真学报，2008，20(14)：3815-3819.