FLUS模型是由中山大学黎夏教授团队开发的是用于模拟人类活动与自然影响下的土地利用变化以及未来土地利用情景的模型,具有很强的时空动态模拟能力,是在其前身——地理模拟与优化系统GeoSOS的基础上的发展与传承,该软件在处理地理空间模拟、土地空间布局优化等方面有着良好的效果。该模型首先采用神经网络算法(ANN)通过土地利用现状数据、驱动因子等来计算各类用地类型在研究范围内的适宜性概率,然后通过耦合系统动力学模型(SD)和元胞自动机(CA)模型以提高模型的适用性,可以减少误差传递的发生,并能有效处理非线性关系^[22-23]。主要的计算步骤如下:

基于神经网络的适宜性概率计算模块计算公式为:

sp(p,k,t)=∑_jω_j,k×sigmoid(net_j(p,t))

=∑_jω_j,k×1/(1+e^-net_j(p,t))(2)

式中:sp(p,k,t)为k类型用地在时间t、栅格p下的适宜性概率; ω_j,k是输出层与隐藏层之间的权重; sigmoid()是隐藏层到输出层的激励函数; net_j(p,t)表示第j个隐藏层栅格p在时间t上所接到的信号。

基于自适应惯性机制的元胞自动机模块:

FLUS模型中,最终的土地之间的转化概率不仅取决于上述神经网络输出的分布概率,还受到邻域密度、惯性系数、转换成本、以及土地之间的竞争所影响。自适应惯性系数是该机制的重要内容,根据各种土地类型与需求数量之间的差距决定各自的惯性系数并在迭代中自适应调整,使得各类用地的数量在模拟过程中向目标数量逐渐接近^[24],其公式为:

Interia^t_k={Interia^t-1_k,|D^t-2_k|≤|D^t-1_k|

Interia^t-1_k×(D^t-2_k)/(D^t-1_k),0>D^t--2_k>D^t-1_k

Interia^t-1_k×(D^t-1_k)/(D^t-2_k),D^t-1_k>D^t-2_k>0(3)

式中:Intreia^t_k表示第k种用地在代时间t上的惯性系数; D^t-1_k、D^t-2_k分表示上一次、上两次迭代时第k种土地利用类型的栅格数量与需求量之间的差距。

在FLUS模型的验证方面主要参考OA,ROC及Kappa等3个参数,值越接近1,模拟准确度越高。

土地利用类型发生转化与驱动因子有着密切的关系,因此在选取驱动因子时要尽量可能地选取与研究区域关联密切的因素。本文通过阅读相关文献再结合研究区的实际情况根据数据的可获取性和一致性、驱动因子的可量化性、地域分异性等原则^[25]共选取区位、自然、社会经济3个方面的因素,共有8个驱动因子,具体见表1。

为了检验选取的驱动因子对种土地利用类型变化相关性,进行Logistic回归分析(表2):

表1 万州区驱动因子选择及含义

表2 Logistic回归分析结果

由表2看出本文中各种地类的ROC检验值均大于0.7,说明选择的上述驱动因子能够描述研究区的土地利用变化情况。

表示土地利用类型之间型转换的难易程度,它是土地类型发生变化的重要因素。根据研究区的土地利用动态实际情况和阅读参考相关文献^[26],设定所有土地利用类型之间均可以相互转化。

在FLUS软件中,经过多次调试及验证,得到模拟精度较高时时的邻域因子参数(表3),参数范围为0～1,因子值越接近1表示该土地类型的扩张能力越弱,越不容易转化为其他类型土地^[27]。

表3 万州区土地利用变化邻域因子参数

参考相关学者^[4,18-19]的研究确定的生态服务价值系数并根据重庆市万州区的实际情况和咨询相关学者的建议,建立适合研究区的生态系统服务价值系数体系(表4)。其中耕地直接采用朱康文对重庆市两江新区的价值系数为6 203元/hm²; 园地的价值系数根据相关研究及专家意见定为7 143元/hm²; 万州区作为重要的生态涵养保护区,林地在生态系统中有更大的作用,因此林地价值系数在参考谢高地等^[18]的基础上增加20%为23 538元/hm²; 草地的生态价值系数参考张凤太等^[19]对重庆主城区和朱康文等^[4]对重庆两江新区的研究取平均值4 261.5元/hm²; Costanza等^[17]、谢高地等^[18]等将城镇和工矿用地的生态价值系数确定为0,而朱康文等^[4]在张凤太等^[19]等的基础上没hm²降低8%,根据专家意见结合研究区的实际情况本研究取对重庆主城区和两江新区研究的平均值为-5 586元/hm²; 万州区水系分布众多,水域系数根据专家意见在朱康文等^[4]的基础上增加5%为6 090元/hm²; 其他土地综合参考几位学者的研究确定为876.5元/hm²。

表4 生态系统服务价值系数元/hm2

生态系统服务价值计算公式^[28]为:

V=∑⁷_iA_i×C_v,i(4)

式中:V表示区域生态系统服务的评估总价值; A_i表示第i种土地利用类型的面积; C_v,i含义表示第i种土地利用类型的生态系统服务价值系数。

生态价值敏感度(C_s)表示生态系统服务总价值对价值系数的依赖程度。如果C_s>1表示该生态系统服务价值富于弹性; C_s<1表示生态系统服务价值缺乏弹性^[29]。计算公式:

C_s=|((V_g-V_f)/V_f)/((C_v,g-C_v,f)/C_v,f)|(5)

式中:V表示估算的生态系统总服务价值,单位为元; C_v表示生态服务价值系数; f和g分别表示在价值系数调整前后生态系统服务价值情况。

生态价值变化度(R_c)反映某类土地利用面积的增加或减少引起的生态系统总服务价值的变化状况。如果R_c>1时表明某类土地利用面积的变化对生态系统总服务价值影响明显,如果R_c<1时表明某类土地利用面积的变化对生态系统总服务价值影响不明显^[29]。计算公式:

R_c=|((V_g-V_f)/V_f)/((A_g-A_f)/A_f)|(6)

式中:A为某土地利用类型的面积。

根据研究区2009年和2017年两期的地类数据,以两期的现状矢量图为底图,利用ArcGIS软件首先对属性表中各地类字段合并在一起,采用Analysis Tools工具对合并以后的数据进行空间叠加分析,打开叠加后结果中的dbf格式文件导入到Excel表格中通过计算得出2009—2017年万州区的各地类的转移概率矩阵(见表6)。

表6 万州区土地利用类型转移概率矩阵

借助表6求得的矩阵,以2009年的土地利用类型数据为基础数据,运用马尔科夫链预测得出万州区2017年和2025年耕地、园地、林地等7种土地的数量。

其中在设定的自然发展情境下,主要根据历史发展趋势进行发展; 在生态保护情境下,是在自然发展的情境下设定园地、林地两种土地利用类型向建设用地转化的数量减少20%,水域向建设用地转化的数量减少25%,同时加入生态保护区的限制因素。具体数量见表7。

表7 万州区2017年、2025年土地需求量预测hm2

根据万州区生态红线的划定成果,将生态红线范围内的区域作为生态红线约束下限制发展的区域(下图中表示为0的区域),主要包括自然保护区、部分风景名胜区、重点生态功能区、水土流失敏感区、饮用水水源保护区、三峡水库消落区、森林公园以及生态公益林等,主要布局在西部的响水镇、甘宁镇,中北部的新田镇、茨竹乡、九池乡,北部的高梁镇、天城镇等和东部的长滩镇、龙驹镇、梨树乡以及地宝土家族乡。在生态保护情境下将该数据作为限制条件输入模型中。进行标准化处理过的生态红线保护区数据见图2。

图2 万州区生态红线保护区

本研究先将研究区2009年土地数据作为基期初始数据,将驱动因子结果进行标准化处理,归一化到[0,1],设置神经网络训练的采样比例为20,占研究区有效采样数的百分之二,神经网络的隐藏层数设置为12。运行软件得到2009年土地利用适宜性概率图集,其中均方根误差RMSE=0.23,说明数据训练结果可信。在FLUS软件的基于自适应惯性机制的元胞自动机模(CA)模块中导入2009年土地利用类型栅格图,加载适宜性概率图集,不加任何限制性因素,输入土地利用变化邻域因子参数,最后设置模型参数,运行软件得到2017年研究区自然发展模式下预测的现状栅格图并与2017年万州区实际现状图对比(图3)。之后借助软件的Precision Validation模块,加载2017年的真实与预测的结果进行对比验证,得出应用总体精度(OA)、FOM指数及Kappa指数。OA和Kappa值通常介于0～1,数值越接近1,模型的,仿真精度越高。当Kappa大于或等于0.75时,模型的仿真精度较高。当Kappa值在0.5～0.75时,模型的仿真精度一般,其中FOM指数大小受模拟年数影响,每模拟年数增加一年,FOM指数增加不大于0.01为标准水平。对比万州区2017年现状和预测图,两张图基本一样,验证精度的指数见表7。由精度验证结果OA值为0.906 2大于0.8且Kappa为0.843 7大于0.75,FOM值处于标准水平,说明FLUS模型可用于万州区土地利用变化模拟^[30]。然后根据表8求得的2025年不同情景下的土地需求数量,以2017年的土地利用现状为初始状态,模拟2025年两种情景模式下的土地利用类型空间分布状况。结果见图4。

根据模拟结果,在自然发展情景下,耕地、园地、草地和水域面积都有不同幅度的减少,其中耕地面积减少最多为1 288.02 hm²,园地、草地、水域分别减少301.53 hm²,144.38 hm²,231.17 hm²,建设用地、其他土地、林地的面积有所增加,其中建设用地面积增加幅度最大为1 919.97 hm²,其他土地和林地只有小幅度的增加。这是由于在自然的发展情境下,随着社会经济的发展和人口的增长,建设用地势必会不断的扩张来占用其他土地类型的面积,主要呈现为在研究区中部沿长江两岸地区的城镇区域建设用地向四周蔓延,主要是由于研究区地势起伏大,山地丘陵较多,建设用地的发展只能沿长江两岸一些地势缓和的区域蔓延发展。其余土地利用情况除其他土地数量增长以外,园地、林地、草地和水域都有小幅度的减少,是区域土地利用变化的特征现象。在生态保护情境下,耕地、园地、其他土地的面积有所减小,其中耕地面积减小幅度最大为1 268.13 hm²,园地和其他土地减少分别为64.83 hm²,246.18 hm²,其他土地类型的面积都有所增加,其中建设用地增加最多为1 093.41 hm²,林地、草地、和水域分别增加342.19 hm²,2.41 hm²,141.04 hm²。这是由于在生态保护情景下,加入了生态红线保护区限制了林地向其他用地类型的转化,增加了林地的稳定性。同时林地、园地、水域用地较自然情形面积有一定的增加,林地增加面积最大主要分布在东南和西北的密林地区,同时建设用地的发展得到了较好的限制,其他用地和草地依然有较小幅度的减少趋势,在该情景下的模拟结果与研究区作为主要生态涵养区和保护区的宏观战略相贴合。总体来说两种情境下的用地变化情况符合重庆地区的发展趋势,对研究区未来的土地利用规划有一定的参考价值。

图3 万州区2017年土地现状及模拟结果

图4 万州区土地利用变化情景模拟结果

通过计算(表8)2017年万州区的生态系统服务总价值约为41.06亿元,2025年设定的两种情景下的总价值有所减少,分别为40.83亿元、41.00亿元。两种情景中耕地和园地的生态服务价值又有所下降原因是城镇的发展、新增建设用地要占用一定的耕地和园地; 林地的生态服务价值有所增加在生态保护情景中增加的价值比自然保护情景多0.079亿元,说明林地在生态保护方面有着更重要的作用; 草地在自然增长情景呈现生态服务价值下降的趋势在生态保护情境下增加了0.006亿元,说明草地在生态保护系统中也具有重要的生态价值; 建设用地在两种情景中都表现出负向价值,在自然增长情境中的负向价值更大比生态保护情景多了0.046亿元,这是由于在生态保护情境下生态红线划定的生态红线保护区域在一定程度上限制了建设用地的发展,有利于经济的可持续发展; 水域在生态保护情境中生态服务价值多了0.009亿元,呈现出增加的趋势。生态红线划定下生态保护情景生态的总价值比自然保护情景下多了0.169亿元说明此种情景下生态红线的划定增加了生态系统的服务总价值,对研究衡量生态红线的划定价值意义重大。

表8 万州区生态系统服务价值亿元

根据表4中确定的价值系数,分别增加和减少5%来计算总价值,计算出敏感性值结果见表9。7种土地的值都小于1且林地的值最高其他土地的值最低,根据敏感性数值可以得出万州区的生态系统服务价值对价值系数的变化缺乏弹性。根据王亚娟等^[31]对宁夏彭阳县的研究得出的结论绝大多数地区的价值系数的变化对地区的生态系统服务的总价值影响较小,说明结果准确。

表9 万州区生价值敏感度

以2017年各种土地利用类型的面积为基础,通过增加5%和减少5%来重新计算各种土地地类的面积,并计算出生态系统服务总价值来分析生态系统服务价值对地类面积的敏感性(表 10)。由表格中的结算结果可以看出来只有林地的价值变化度接近于2.0,说明在万州区林地的面积变化对于生态系统服务价值的影响最为显著,其中园地、耕地和水域的价值变化度大于0.5说明这3种土地利用类型的面积变化对于生态系统服务价值有一定的影响,其余的土地类型的影响较小。所以政府部门在制定相关规划时可以考虑通过增加水域用地和林地的比例来提高区域的生态环境质量。万州区作为重庆乃至全国重要的生态保护和生态脆弱区,林地面积较大,在土地利用中占有重要的比例,而生态保护红线的划定可以有效地巩固和加强林地等增强生态系统服务价值地类的面积,对于保护生态功能和生态敏感区具有重要的意义。因此政府等决策部门应该加强对生态空间的保护,严格限制划定的生态红线内的区域的开发,构建良好的生态环境,实现社会经济的可持续、健康发展。

表 10 万州区生态价值变化度