随着气候变化、人口激增、环境恶化和能源短缺等一系列全球性问题的持续发酵，土地利用/覆盖变化(land use/cover change，LUCC)的研究逐渐成为全球环境变化研究中的核心部分^[1]。其中，应用多种方式构建LUCC的动态变化和预测模型，一直是的LUCC研究的热点^[2]。近年来，随着遥感影像向着高分辨率的方向发展，使得土地利用变化与遥感影像在时空尺度上具有了更高的一致性，利用高时间、高空间分辨率的影像数据来研究多时空尺度的土地利用变化已成众多学者的选择^[3-4]。多尺度的遥感影像同样被应用于土地利用/覆盖变化的各种预测模型^[5]。在这些预测模型中，Logistic回归模型、灰色模型、Markov模型等数量变化预测模型以及元胞自助机、多智能体系统、CLUE-S等空间预测模拟模型，均得到了较为广泛的应用。然而目前，国内外学者更青睐于采用多模型结合的耦合模型方法来深入到对土地利用的预测模拟研究中，如CA-Markov模型，MAS-CA模型，Logistic-CA-Markov和WLC-CA-Markov模型等。

从已有的LUCC相关的研究来看，理论和方法均已趋于成熟，但仍然存在一些问题，如：①土地单元面积的选择大小不一，未能对土地预测模拟中产生的尺度效应问题有效考虑；②研究区选择多为热点区域，空间尺度上至国家、江河经济带、城市结合群，下至市、区、县不等，而一些与水文、气候条件响应强烈的小流域尺度区域往往缺乏重视。位于海河上游的清水河流域具有典型小流域尺度特征，而且在过去20 a经历了剧烈的土地结构变化。一方面，该流域生态环境比较脆弱，全国400 mm年等降水量线从流域穿过，整个流域都处在半湿润向半干旱过渡的生态环境脆弱区，流域内的植被条件和水文环境对于气候变化和人类活动都极为敏感^[6]；另一方面，由于该区城乡建设加快以及2022冬奥会的筹办，该区已计划建设76个场馆相关项目，其中滑雪旅游和竞技滑雪可供开发的面积多达300 km²，2018—2019年仅雪季期间游客涌入人次多达2.80×10⁶，在未来冬奥会期间该情形将更加严峻。上述种种因素会使得流域内土地格局变化十分显著，相关的水文、环境条件也随之改变。然而该区地广人稀，地方政府在面对如此大量繁多的土地资源时对于土地的建设管理不能一概而论，因而需要有相应的针对性地制定土地资源的管理和可持续利用措施。目前针对该流域的土地利用变化的相关研究集中在崇礼区，而非与水文条件密切相关的清水河流域，时间也较为滞后，特定地关注其流域空间单元尺度效应以及模型预测的研究尚不多见。基于此，本文以清水河流域为研究区，利用遥感影像解译的土地利用数据分析流域内土地利用时空变化特征；同时，考虑到社会数据的可比性、可获得性等因素，选择与清水河流域边界几乎吻合的崇礼区作为土地预测模拟的研究区，结合该区域土地数据以及栅格化的驱动因子数据共设计了6个空间尺度并从中优选出最佳尺度，在该尺度下通过耦合Logistic-CA-Markov模型进行流域内的土地利用变化预测模拟，以期对清水河流域土地利用/覆盖方面的研究做必要补充，同时为当地有关土地资源开发利用决策提供一定参考。

1 材料与方法 1.1 研究区概况

海河上游清水河流域介于114°47′—115°31′E，40°49′—41°16′N，位于河北省西北部，地处内蒙古高原与华北平原过渡地带。流域内气候属东亚大陆性季风气候中温带亚干旱区，冬季干冷，春季回温较快，大风日数较多。多年平均降水量约442 mm，降雨主要集中于6—9月；多年平均气温3.9 ℃，年内季节均温最高相差40 ℃^[7]。流域总面积约2 200 km²，东高西底，海拔高度约为780~2 200 m^[8]，地貌属坝上坝下过渡型山区，山峰陡峭，海拔多在1 500~2 000 m之间。受山区地貌影响，局部冰雹、暴雨灾害时有发生。在地理、气候、人为等多重综合因素影响下，清水河流域水土流失较为严重，生态环境极为脆弱，曾在1983年被国家列入永定河流域重点治理区。

清水河流域崇礼区人口在1990—2018年期间由1.20×10⁵增长至近1.30×10⁵，整个行政区总面积约2 326 km²，人口密度在51.6~55.9人/km², 可谓是“地广人稀的燕赵大地”。2017年当地5大行业产值结构中，农、林、牧产值比例分别大致为53%，21%，25%，渔业和农林渔牧服务业占比均不到0.01%。其中，蔬菜产值占5大行业产值的近42%。

1.2 数据源及其处理 1.2.1 数据源

研究选取的1990，2000，2009，2018年4期Landsat OLI/TM影像以及DEM数据均来源于地理空间数据云(http://www.gscloud.cn/)，分辨率均为30 m。4期影像成像时间均选在10月份避开集中降雨期，云量均低于2%且成像质量较好。鉴于清水河流域边界和张家口市崇礼区行政边界基本吻合，流域内的社会数据和地理数据来源于《张家口统计年鉴》和国家基础地理数据中心(http://www.ngcc.cn/ngcc)。

1.2.2 影像数据处理

在ENVI软件中完成影像的几何校正、辐射校正、裁剪镶嵌等预处理工作；分类方法采用监督分类中分类精度和kappa系数较高的最大似然法^[9]；土地利用分类根据中国科学院土地利用/覆盖分类体系中一级分类标准并参考流域内的实际情况分为5类：耕地，林地，草地，水域，建设用地；将分类后的影像进行小斑块处理，并根据Google Earth中的邻近时段相同地区的高分辨影像结合目视解译法做分类后处理。通过生成的随机样本点在Google Earth对比目视判读进行分类精度验证。最后得到的4期影像分类精度均在85%~90%，可以满足后续的分析要求。

1.2.3 驱动因子栅格化

土地利用/覆盖变化驱动主要受自然和人类活动两方面因素影响^[10]。根据驱动因子的数据可获取性、时空一致性、因子显著性以及可量化性等原则，结合清水河流域崇礼区内的实际情况，共选择了7个因子：坡度、坡向、高程、人口密度和距河流、道路、村落距离。以上数据采用的ArcMap软件中的Aspect，Slope，Kernel Density，Distance等工具栅格化所得。因所选区域较小，降雨、温度等气象数据在栅格内变化较均质，故暂不考虑。同时将土地利用的5个地类数据单独成层并通过Reclass工具进行二值化，每一图层对应的地类栅格值为1，非该地类均为0。

1.2.4 空间尺度设计

栅格数据单元格大小均以30 m为基准，在ArcMap中利用Resample工具，设计了30 m×30 m，60 m×60 m，100 m×100 m，150 m×150 m，200 m×200 m，250 m×250 m共计6种尺度的栅格数据，部分空间尺度区间不一致主要是综合考虑尺度间隔的均匀性以保证更细节的变化特征同时在更大尺度范围进行对结果进行比较。所有的栅格数据经过Raster to ASCII工具转换成文本数据，具体过程是在6种不同尺度下，将各个土地利用类型二值化数据分别与7个驱动因子数据结合对应生成5组，经过Dyna-CLUE软件全部转换成SPSS中待分析的共计30个文本数据。

1.3 研究方法 1.3.1 土地利用时空变化研究

本文采用土地动态度来定量分析流域内土地利用的变化剧烈程度，用土地转移矩阵来研究各阶段时期内土地的动态转移变化。单一土地利用动态度指某一区域内在一定时间范围内某种土地利用类型的数量变化情况，此指标可定量描述研究区内某种土地类型在研究时段内变化的速率。表达式如下：

$K=\frac{U_{b}-U_{a}}{U_{a}} \times \frac{1}{T} \times 100 \% $

(1)

式中：K为研究时段内某一土地利用类型动态度；U_a，U_b分别为研究初期和研究末期某一种土地利用类型的数量；T为研究时段长^[11]。

土地转移矩阵反映了某一区域在研究初期和研究末期各类土地类型面积互相转化的动态信息，不仅可以定量地表明不同土地利用类型之间的转化情况，还可以揭示不同土地利用类型间的转移速率。转移矩阵表达式参考乔伟峰等^[12]的研究，表达如下：

$S_{i j}=\left[\begin{array}{cccc}S_{11} & S_{12} & \cdots & S_{1 n} \\ S_{21} & S_{22} & \cdots & S_{2 n} \\ \vdots & \vdots & \vdots & \vdots \\ S_{n 1} & S_{n 2} & \cdots & S_{n n}\end{array}\right] $

(2)

式中：S为面积(km²)；n为转移前后土地利用类型数；i，j(i，j=1，2，3，…，n)分别表示转移前后的土地类型；S_ij表示转移前的i类土地转换成转移后j类土地类型的面积。

1.3.2 土地利用预测模拟

尺度效应是指当空间数据经聚合而改变其粒度或栅格单元大小时，分析结果也随之变化。反映土地利用类型及其面积变化的数据必然涉及尺度效应问题^[13]。在土地利用预测模拟中，选择合适的空间尺度能够使得模型预测精度更高，结论也更合理。Logistic回归模型是一种广义的线性回归分析模型，被广泛应用于土地利用分析中，能够很好地揭示土地变化和驱动因素之间的定量关系^[14]。本文采用Logistic回归分析方法得到受试者工作曲线(receiver operating characteristic，ROC)作为参考来优选出最佳空间尺度^[15]，以ROC曲线下的面积值作为精度指标来检验回归方程的解释效果^[16]，从而判定预测模型的模拟精度。当0.50＜ROC值＜0.70之间时，回归方程具有一般解释效果；当ROC值＞0.70时，方程具有良好解释效果，值越高代表解释力度越大。

元胞自助机(cellular automata，CA)模型是一种在空间上相互作用、时间上又具有因果关系的一种网格动力学模型，具有处理复杂空间系统的能力^[17]。Markov模型基于Markov链，因其较好的稳定性和无后效性，能够预测土地利用变化中各个时刻的变动过程^[18]。CA-Markov模型综合了Markov模型的时序预测和CA模型的空间分布模拟，广泛被应用于土地利用数量变化、空间格局变化的模拟中^[19]。

通过耦合Logistic-CA-Markov模型，首先对崇礼区2009年的参数文本数据进行二元Logistic回归分析并分析各个尺度下的ROC曲线从中优选出最佳尺度，在该尺度下，将回归分析得到的5种地类分布适宜性图通过IDRISI软件中的collection editor工具组合生成各类土地分布概率适应性图集并作为CA模型的转换规则，同时以2009年土地利用分类结果和Markov工具中生成的2000—2009年面积转移矩阵作为预测模拟初始参数，共同参与到CA-Markov模型进行崇礼区2018年的土地预测模拟；将生成的土地分类预测结果与2018土地实际分类与进行叠置分析，采用Kappa系数、混淆矩阵等指标验证模型预测模拟的精度，在模拟精度良好的情况下利以同样方法进行2027年的预测模拟。

2 结果与分析 2.1 土地利用时空变化分析

1990—2018年清水河流域4期遥感影像解译土地利用分类结果如图 1所示，具体数据如表 1所示。

从整个时期的土地组成结构来看，清水河流域土地类型总体上主要以林地、草地、耕地为主，三者面积之和占据流域总面积96.00%以上，水域和建设用地面积在2000年共计占比不到2.00%，在2018年占比达到3.96%；从流域的土地结构可以分析，其与整个研究区的地貌存在一定相关性，清水河流域最低海拔在800 m左右，群山峻岭密布，人类活动区域相对较小，适宜林地、草地的生长；从总的变化形势上看，1990—2018年清水河流域土地利用变化中，林地、草地、和水域是处于减少态势，耕地面积变化不明显，而建设用地面积呈显著上升趋势，近30 a来增加了近51.72 km²；分阶段来看，1990—2000年，耕地变更的面积最大，10年间增加了约44.74 km²，建设用地变化的面积最小，仅3.20 km²；2000—2009年期间仅有耕地面积在减少，其余四种土地类型面积均有不同程度的增加；2009—2018年间，流域内仅有耕地和建设用地面积处于增加态势，增加面积分别为20.24，18.02 km²；林地、草地和水域面积分别减少了12.72，20.76，4.84 km²。

单一土地动态度反映的是某一类型土地变化的剧烈程度。由表 1可知，在1990—2018年整个时期内，清水河流域内建设用地的变化速率最为剧烈，其中以2000—2009年的变化较为显著；林地面积无明显变化，其动态度最高仅为0.21%；分阶段比较，2000—2009年土地整体变化剧烈程度较显著于其余2个时期，且在每一阶段内，土地变化速率大小按地类依次排序为：建设用地＞水域＞耕地＞草地＞林地。综合各个地类单一动态度可以看出，近30 a来清水河流域内农业用地更迭趋势比较平稳，而建设用地的增加较为迅猛，当地社会对交通用地、工矿地和居民用地等建设用地开发需求较大。

2.2 土地利用转移矩阵分析

考虑到社会数据的可获取性，仅对2000—2018年3期的清水河流域土地转移情况进行分析。通过叠置分析得到不同阶段的土地利用转移矩阵如表 2所示。从整个阶段来看，清水河流域在2000—2018年期间，耕地、草地和建设用地的的转换较明显，约有22.30 km²的草地和29.63 km²的建设用地由耕地转换而来；水域转为其它类型土地的程度不够显著，最多仅有1.17 km²水域的转出；农业用地到建设用地的转入较为显著，转入面积约53.25 km²，从实际情况考虑是由于近年来当地人口增多和冬奥会场地筹建等因素对建设用地的需求增大。

分阶段来看，2000—2009年和2009—2018年各土地流转程度比较相似，均有大量耕地和草地互相转化，转化面积在40~50 km²不等；其中，2000—2009年约有5.37 km²的水域由农业用地转入，这与近年来当地洪水等自然灾害频发以及农业经济发展对用水需求增多等因素进而引起对水库、水坝等水利设施的新建扩建有关；建设用地也基本由农业用地转入而来，转入面积的大小按地类排序依次为：耕地＞草地＞林地。其中有一定程度的林草地转入耕地，该区2000年后虽然已经实行退耕还林还草政策，但少量开荒复垦现象仍然存在，同时不排除分类精度带来的影响。

2.3 尺度模拟结果

经SPSS中二元Logistic回归分析得到结果如图 2所示。通过对比各个地类的ROC值可以观察到，随着空间尺度的不同，各个地类的ROC值也随之产生不同变化，即驱动因子对各个地类的解释水平随着空间尺度的不同产生了“尺度效应”的特征。

由图 2可知，在空间尺度为100 m×100 m时，水域、耕地和建设用地ROC值出现转折趋势，其中，水域ROC值的转折尤为明显且达到了峰值，说明该尺度下水域与驱动因子的拟合方程的可信水平较高；此外，在100 m之后，草地和建设用地ROC值接近最大值且变化趋势近乎平稳，尺度的改变对模型回归拟合的影响趋于微小。以上结果表明，该尺度下的二元Logistic回归模型相比其余5种尺度其拟合效果更优，因此可以认为100 m×100 m是进行本次土地预测模拟的最佳空间尺度。由图 3可知，在选定的最佳空间尺度下，各类土地ROC值分别为耕地0.67，林地0.65，草地0.66，水域0.83，建设用地0.79。通过比较，耕地、林地、草地的ROC值均接近0.70，水域和建设用地均大于0.70，回归方程拟合效果较好，ROC值满足检验标准，回归模型的结果可以作为CA-Markov模型中的转换规则进行下一步的土地利用变化预测模拟。

2.4 预测结果与精度验证

在ENVI软件中对2018年土地实际分类和预测分类结果。进行混淆矩阵运算，结果如表 3所示，得到预测图与实际图中像元一致的个数占比为84.68%，kappa系数为0.78。一般而言，当kappa系数大于0.75时可以认为两种分类图像一致性较高，模型模拟精度满足预测要求，可以进行下一时期的模拟。将数据对应的年份均更新至2018年，以此作为模拟的初始参数进行2027年的土地预测模拟，结果如图 3所示。

表 4为清水河流域崇礼区土地利用的预测结果。由表 4可知，2018年预测模拟得到的崇礼区土地利用变化面积中除了建设用地的误差稍大以外，其余4类土地预测的误差均在1.50%~6.02%，误差在可接受范围内。结合当地实际情况分析，对2018年建设用地面积预测误差较大的主要原因是因为崇礼区于当年新建了2022年冬奥会的滑雪场地，从而导致预测值的误差较大。因此，排除该因素的干扰，该模型预测模拟的实际精度会更高。对2027年进行预测模拟得到的土地分类结果如表 4所示，预计未来9 a，崇礼区内的耕地、林地和建设用地均会增加14.00~17.00 km²不等，草地和水域面积将分别减少42.44，5.43 km²。

3 讨论与结论 3.1 讨论

本研究中，由于所选研究区域较小，加之影像分辨率的限制，经目视判读选取的分类训练样本中会无可避免地掺杂进混合像元，对土地的分类精度造成一定影响；清水河流域自2000年实施退耕还林还草政策后，林地开荒、草地复垦现象得到显著改善但并未完全杜绝，经各期土地数据叠置分析后检验可知，受分类精度的影响出现部分林草地逆向转化为耕地的不合理转入，同时2000—2010年耕地较为显著下降印证了该政策取得了一定成效。此外，占地面积较小的水域和建设用地均在分类后处理中由目视解译所得，影像时间的选取已考虑降雨期和冬季结冰的影响均选择在10月份，所得分类数据中，清水河流域农业用地占地比例与建设用地面积与钟良子等^[20]，王磊等^[21]所得清水河流域以及崇礼区的数据结果基本一致。经过Google Earth高精度实地影像检测研究区内未利用地甚少，故未对其进行分类，同时考虑流域内河流时常发生断流、降水量不稳定，如2009年崇礼区降水量仅为334 mm，而其前后两年却高达540 mm以上，即使同一时期不同年份的河流面积变化仍可能不同，本研究参考水域二级分类将水库坑塘、水域设施均也纳入其中。在可控的分类精度影下本文结果具有一定参考性。

本文出于对数据可获取性以及可量化性等原则的考虑，假定了降雨、蒸发等气象因素在单元栅格面上对该类小空间尺度区域土地流转的影响较为均质，未能对该类因素在驱动土地利用变化的影响程度上做深入研究。在未来的研究中会考虑选择更高分辨率的影像和高精度的分类算法，对驱动因素的考量也会更加科学。

经分析，同类土地利用预测模拟研究中研究区多集中在市及市以上区域，以行政区为对象的研究对人为、社会政策因素驱动土地格局发展上具有很好的解释水平。经参考流域面积相近且同为CA-Markov模型的模拟研究^[22]，其土地结构与动态变化与本文相似，均较为稳定；对于地域较大的研究区域笔者认为其预测结果是由CA-Markov模型中单个的元胞基于Markov链与Logsitic模型概率公式对于空间单元的独立计算组成，市级以上的预测模拟近似于市内各县域模拟结果的共同叠加，其结果包含多个稳定子区域模型的误差累计，而以小范围的县域为研究区其驱动因素对于Logistic模型的表达更敏感，小尺度区域空间异质性也相对较低减弱了结果的不确定性。鉴于清水河流域崇礼区内2018年10月许多建设项目并未完工，而Logisric-CA-Markov模型是基于前一时期的参数进行后一时期的模拟，其结合流域内崇礼区前一时期自然、人为政策驱动因素下土地的变化，对未来主要由冬奥会驱动下的土地格局变化仍然具有相对程度的刻画。在驱动因素方面选择了3个自然因素和4个人类活动因素，本研究主要关注7个驱动因素与各个地类在Logistic模型中不同尺度下ROC曲线的变化以做尺度效应分析，未对7个因素与各个地类之间的驱动关系进行定量描述；模拟的时间点是基于模型预测时段与前一阶段保持一致的考量。综合比较，本文中对崇礼区2027年的土地格局的预测模拟具有一定的参考意义。此外，经查证目前崇礼区政府针对建设土地严格管控一系列措施的出台，预计在2022年冬奥会步入尾声时，清水河流域内建设用地的增加将会进入饱和状态，如果考虑该影响，则模型计算出的2027年建设用地面积预测值可能偏高。

3.2 结论

(1) 近30 a来，清水河流域土地利用类型主要以耕地、林地和草地为主，三者面积共计占比在96%以上，其与清水河流域平均海拔较高、地貌多为陡峭山区、人类可活动区域较少等因素有关；近几年来建设用地的转入较为明显，主要是由于冬奥会的筹备新建了滑雪场以及相关场馆加大了对建设用地的开发。

(2) 对二元Logistic回归模型中的ROC曲线进行分析，可以认为适合清水河流域这类小尺度流域研究区进行预测模拟的最佳空间尺度为100 m×100 m，该尺度下构建的方程模型对土地利用变化与驱动因素之间的定量描述具有较高的解释水平。

(3) 至2027年，预计清水河流域内耕地、林地和建设用地面积将呈现增加态势，草地会大幅度减少。通过调查2018年崇礼区的实际情况后考虑，当地虽然已建成万龙、云顶、太舞等冬奥会滑雪场地，但还未达到规划用地规模的最高值，预计未来一段时间流域内的农业用地将在冬奥会的用地需求下持续转化为建设用地。