根据研究区土地利用特点及用地需求,结合《宁武县土地利用总体规划》、《静乐县土地利用总体规划》,通过不断的调试Markov模型的转移概率以及转变适宜性图集,预测不同情景下2030年研究区的土地利用类型的数量及空间分布变化。

情景1:自然发展情景。根据2010—2017年土地利用转移概率矩阵以及适宜分布概率,以2017年土地利用现状数据为底图,基于CA-Markov模型预测2030年在自然发展情景下各土地利用类型的面积及空间分布。

情景2:生态保护情景。研究区为黄土丘陵区重要的退耕还林区域,也是华北平原生态屏障的核心区。为了进一步巩固退耕还林的成果,未来的土地利用要进一步加强对生态用地的保护。在该情景的设置中,林地、草地、水域向建设用地和未利用地的转移概率降低100%,耕地向林地的转换概率增加30%; 此外要严禁生境退化,将区域自然保护区、国家森林公园、高山湖泊群、地质公园等作为约束条件,限制生态用地的任意转换。

情景3:耕地保护情景。耕地是粮食安全的载体,该区域为重要的农业开发县,该情景是为保护耕地设置。一方面,基本农田区范围内的耕地严禁转出,将其作为约束条件; 另一方面,修正转移概率,将耕地向建设用地和未利用地的转移概率分别降低90%和100%,向林地、草地、水域的转移概率降低30%。

Logistic回归模型是解释自变量和因变量的关系和强度的空间模型,目前已被广泛地应用在土地利用变化的研究中。本文选择ROC曲线法检验空间Logistic模型的预测结果,ROC曲线下面积介于0.5～1.0,表明模型的可信度较高。其基本拟合方程^[22]为:

LOG[(P_i)/(1-P_i)]=β₀+β₁X₁+β₁X₁+β₂X₂+…+β_nX_n(1)

式中:P_i为每一个栅格可能出现某种地类i的概率; X为备选驱动因子。

CA-Markov模型的结合能够全面考虑自然和人文因素的影响,实现土地利用变化的动态模拟,使之既能模拟复杂系统空间变化的能力,也能发挥土地利用需求长期预测的优势。本文的元胞空间为汾河源头区域,元胞的状态为6种地类。借助Logistic回归诊断每一个栅格单元可能出现某种地类的概率,从而选择对土地利用变化影响较大的因子以及它们之间的定量关系与相互作用关系。

模型需要输入相关参数,数据处理如下:(1)运用Markov模型分析6种土地利用类型转换矩阵和转移状态概率。(2)驱动因子选取。考虑到指标的空间均衡性以及数据的可获取性,在提出相关性显著指标的基础上,本文选取海拔、坡度、坡向、土壤类型、距城镇的距离、距农村居民点的距离、距河流水系的距离、距道路的距离、年降水量、人口密度和人均GDP因子。(3)分布概率适宜图:以6种土地利用变化数据为自变量,11个空间化的驱动因子为因变量,运用IDRISI 17.0软件制作6种土地利用类型的适宜分布概率。(4)地类转换规则设置。其中,自然发展情景不设置转化规则,按照默认5×5邻域转化; 生态保护情景仅允许生境正向演替; 耕地保护情景禁止耕地的转出,但是允许生境的逆向演替。

生境质量指在生态系统为个体、种群、群落等提供生存繁衍所需条件的潜力。生境质量优劣与生物多样性的丰富程度呈正相关。InVEST模型假定生境质量是连续的,且在一定程度上生境质量的空间分布与生物多样性的空间分布基本具有一致性。因此本文运用InVEST模型中的生境质量模块评估,从外界威胁因子(表1)和生境敏感性来评价生境质量的优劣^[23]。

Dxj=∑Rr=1∑yry=1(Wr/∑Rr=1Wr)ryirxyβxSjr(2)

式中:D_xj为土地覆被类型j栅格x的生境退化度; R为威胁因子数量; y_r为r威胁栅格图上的一组栅格; W_r为威胁因子权重; r_y为栅格y的威胁因子值; β_x为威胁因子可达性; S_jr为生境类型j对威胁因子r的敏感性。i_rxy为生境类型栅格x与y的距离函数(包括指数衰减函数和线性衰减函数),计算公式为:

i_rxy=1-((d_xy)/(d_{r max}))(3)

i_rxy=exp((-(2.99)/(d_{r max}))d_xy)(4)

式中:d_xy为栅格像元x与y间的距离; d_{r max}为威胁r的最大威胁距离。

在InVEST模型中用Q_xj表示生境质量指数:

Q_xj=H_j[1-(D^z_xj)/(D^z_xj+k^z)](5)

式中:H_j为土地覆被类型j的生境适宜性; k和z为比例因子,其中k为半饱和常数,z为归一化常量。

土地利用变化类型对生境贡献率指某类土地利用发生变化使生境质量发生改变的比率。以2017年的生境质量指数为基期年,分析不同情景下生境质量指数的变化和土地利用变化对生态质量服务的贡献率,其表达式^[24]为

Q_LH=(Q_H,t+1-Q_H,t)×L/S(6)

式中:Q_LH为土地利用变化的生境贡献率; Q_H,t+1,Q_H,t分别为某种土地利用类型在变化初期和末期土地利用类型的生境质量指数; L为该变化类型的面积(km²); S为土地利用总面积(km²)。

表1 威胁因子属性

在剔除相关性指标和以及考虑到指标的可获取性的基础上,本文从自然条件和社会经济等角度选取11个驱动因子,分别为海拔、坡度、坡向、土壤类型、距城镇的距离、距农村居民点的距离、距河流水系的距离、距道路的距离、年降水量、人口密度和人均GDP因子。利用Logistic模型,分别计算了影响6种土地利用类型空间分布的11个驱动因子的回归系数。然后选择ROC曲线法检验Logistic模型的精度。ROC值介于0～1.0,若ROC值在0.5～1,则说明模型具有很好的预测能力; 相反如果ROC<0.5,说明构建的Logistic模型不具有预测功能。从表2中可以看出,6种土地利用类型的回归中的ROC值都在0.85以上,表明本文所构建的回归模型具有可信性,11个驱动因子具有很好的解释能力。耕地的分布概率与海拔、坡度、人口密度呈较强的负相关关系; 随着坡度和海拔的增加,林地的分布面积也增大; 坡度与草地的分布概率呈较强的正相关关系,与水域和建设用地呈较强的负相关关系,皆与实际情况相符,表明模型可信度较高。

表2 汾河源头区域逻辑回归结果及ROC检验值

为了验证模型的模拟结果,基于IDRISI软件中的CA-Markov模块,以2010年为基期年,利用2000—2010年的土地利用转移面积和概率矩阵,通过Logistic回归模型获得2000—2010年土地利用变化的转移适宜性图集,从而模拟2017年的土地利用的空间分布情况,并将其与2017年真实的土地利用分布数据相比较,以验证CA-Markov模型模拟的可靠性。经检验,Kappa系数为0.835,精度较高,说明CA-Markov模型参数设置较合理,可用来模型未来的土地利用空间分布情况(图1)。