基于1990—2019年Landsat影像的干旱区绿洲土地利用变化与模拟

宋奇， 冯春晖， 马自强， 王楠， 纪文君， 彭杰

(1.塔里木大学农学学院，阿拉尔 843300； 2.北京大学地球与空间科学学院遥感与地理信息系统研究所，北京 100871； 3.浙江大学环境与资源学院，杭州 310058； 4.中国农业大学土地科学与技术学院，北京 100083)

0 引言

土地利用变化是影响生态环境的重要因素之一，已成为全球变化和可持续发展研究的重要内容[1-2]。土地利用变化改变了区域的物质循环、自然景观和各种生态过程，进而导致全球生态系统发生变化[3-5]，使得土地利用变化与模拟成为研究区域和全球生态系统变化的主要趋势和潮流[6-8]。为降低区域土地利用变化对自然和社会环境的影响，有必要对连续长时序的区域土地利用变化进行分析并对未来土地利用变化情况进行模拟，以期为土地利用变化的调控和管理提供参考。

国内外研究者对土地利用变化的研究主要集中在土地利用时空变化特征[9]、土地利用变化机制或驱动力[10]、土地利用变化的环境影响[11]等方面，少有涉及连续长时间序列的土地内部结构详细转移过程与模拟。从研究区域上看，多数研究集中在东部沿海经济发达的城市区域[12-13]，而着眼于西北干旱区经济较为落后的典型人工绿洲区域的研究甚少。同时，城市区域的土地利用变化研究所使用的数据源多为典型时间断面的单时相遥感影像和单一的社会统计数据[14]，不仅未能有效捕捉连续性信息和解译更精细土地利用分类结果，特别是造成一些关键突变点或拐点信息的遗漏，而且也不能全面反映研究区地形地貌、路网和气象情况。因此，基于连续长时间序列的多时相遥感影像和多数据源(地形地貌、路网、气象和社会统计数据等)对人工绿洲区域的土地利用内部结构转变过程及定量模拟机制研究仍需进一步加强。

阿拉尔垦区位于我国西北干旱地区，属于典型的人工绿洲区域，垦区内荒漠区域的自然土壤表层常覆有一层盐结皮以防止水分流失，保护土壤内部结构，过度的土地开发利用破坏了表层盐结皮，加速了水分流失，影响了土壤水分和大气循环系统。同时，阿拉尔垦区紧邻塔克拉玛干沙漠，自然环境恶劣，常出现大风、扬尘和沙暴等极端天气，自然状态下的盐结皮能抵抗10级风的侵蚀，而表层盐结皮破坏后3级风就能引起扬尘，严重影响空气质量[15]。此外，植被也是研究区的天然生态屏障[16]，但自然植被对生态系统的维护能力大于人工植被，而阿拉尔垦区内的植被以人工植被(农作物和经济林)为主，自然植被(林草地)为辅，且农作物多为1 a生植物，极易受到破坏。因此，有必要对垦区进行土地利用变化和模拟研究，为未来土地利用规划提供科学依据。

本研究以西北干旱区典型人工绿洲阿拉尔垦区为例，收集1990—2019年所有可用Landsat遥感影像、地形地貌、路网、气象和社会统计数据，将逐年各月份影像进行多时相合成，使用支持向量机(support vector machine,SVM)分类方法，辅以2019年高分二号遥感影像及野外勘察对多时相合成后影像进行解译分类，获得逐年土地利用分类图，从面积变化、类型转化、空间动态变化3方面进行土地利用变化分析，选用CA(cellular automaton)-Markov模型对阿拉尔垦区2010年土地利用情况进行模拟并与真实情况比较，验证模拟模型的可靠性后，模拟2050年和2080年的土地利用情景，探究长时序土地利用类型面积的突变情况和驱动因素，进一步说明了CA-Markov模拟模型的优势，为典型人工绿洲区域的生态环境综合管理和宏观的土地资源与社会经济统筹管理提供理论依据。

1 研究区概况

阿拉尔垦区位于新疆维吾尔自治区阿克苏地区，地处塔克拉玛干沙漠和天山山脉的交汇过渡地带，地理坐标为E80°30′～81°58′，N40°22′～40°57′(图1)，总面积4 105.92 km2。阿拉尔垦区地势相对平坦，气象多变，常有沙暴、扬尘和霜冻等灾害，全年降水少蒸发大； 生态系统属于荒漠—绿洲型，土壤表层常常附有一层盐结皮，以防土壤水分蒸发； 垦区种植结构单一，大面积种植棉花。

图1 研究区地理位置

2 材料与方法

2.1 遥感数据源及其预处理

从美国USGS网站中(https: //earthexplorer.usgs.gov/)获取1990—2019年间阿拉尔垦区(轨道号为146/32，空间分辨率为30 m)的所有Landsat遥感影像，并从中筛选出云覆盖度低于40 %的可用影像，共259景，各影像云覆盖度情况和每年影像数情况如图2所示。

(a) 云覆盖率 (b) 每年影像数

从图2可看出各影像云覆盖度情况多低于10 %，由于Landsat系列传感器不断更新，因此本文1990—1998年选用Landsat5，1999—2012年选用Landsat7，2013—2019年选用Landsat8影像，避免了多传感器合成所带来的误差，且每年影像数至少有3景。

为了进一步提高土地利用分类精确度，将所有影像逐一进行预处理，其中包括辐射定标、大气校正、配准、裁剪等操作，并提取各景影像的归一化植被指数(normalized difference vegetation index,NDVI)，计算公式为：

NDVI=(NIR-R)/(NIR+R)，

(1)

式中：NIR为Landsat近红外波段反射率;R为红光波段反射率。

将每年各月份影像进行NDVI最大值合成，得到逐年植被覆盖最全面的情景，从数据源方面提高分类精度，相关影像的合成流程如图3所示，对影像合成前后的局部放大对比中可以看出，NDVI最大值合成后的影像植被覆盖度更高，利于解译分类。

图3 多时相遥感影像合成

2.2 土地利用分类方法

参考《新疆生产建设兵团统计年鉴》(1990—2019年)中第一师各土地利用面积统计情况和已有的新疆土地利用分类体系[17-18]，结合阿拉尔垦区土地利用类型特点，将研究区土地利用类型分为6类。各土地利用类型的地表特征、影像特征和实地照片情况如表1所示。

表1 阿拉尔垦区分类体系及解译标志

(续表)

以2019年为例，根据阿拉尔垦区分类体系及解译标志情况建立各土地利用类型的分类样本(图4)，图上括号内为各类样本所占的像素个数，共计1 024 951个分类像元数。根据各分类样本，辅以收集的相关统计数据、高分影像和实地调查，采用监督分类中的SVM分类方法进行土地利用分类。选用全连接条件随机场(fully connected conditional random field, FC-CRF)[19]这种分类后处理方法对分类后影像进行分类后处理。

图4 每类样本点个数和分布

由于本文采用了多源遥感影像(Landsat8,Landsat7和Landsat5)进行分类，此处以2019年的影像代表Landsat8，2010年的影像代表Landsat7，1990年的影像代表Landsat5，进一步分析不同数据源的分类结果(图5)和分类精度(表2)。从中可以看出不同数据源的总体分类精度和Kappa系数均大于80 %和0.8，并且各用地类型的制图精度和用户精度均大于75 %。整体而言，此分类方法得到的分类后结果优异，可用于后续相关研究分析。

(a) Landsat8(2019年)(b) Landsat7(2010年)(c) Landsat5(1990年)

图5 不同数据源的分类结果

表2 不同数据源的精度评价

2.3 土地利用变化分析方法

从各土地利用类型的面积变化、类型转化和空间动态变化3方面对阿拉尔垦区土地利用变化展开分析。基于ArcGIS软件的空间分析功能，对各时间段起始年份和终止年份的分类后影像进行分析，并制作各土地利用类型变化的转移矩阵和空间动态分布图。

2.4 土地利用模拟数据的收集与处理

为提高土地利用模拟精度，参考相关文献[20-22]并结合阿拉尔垦区的土地利用现状，选取高程、坡度、坡向、距高速、铁路、公路和水系距离、气温、降水、人口、国内生产总值(gross domestic product,GDP)、全社会固定资产投资、第一产业、农业生产总值和棉花价格共15种影响因素进行土地利用模拟。

高程数据从中国科学院计算机网络信息中心国际科学数据镜像网站中获取，在ArcGIS软件中对高程数据进行重采样得到坡度和坡向信息，利用Near工具得到每个栅格到高速、铁路、公路和水系的距离； 气温和降水数据从中国气象科学数据共享网中获取，并进行反距离权重(inverse distance weight,IDW)空间插值得到气温和降水的空间分布情况； 人口、GDP、全社会固定资产投资、第一产业、农业生产总值和棉花价格从《新疆生产建设兵团统计年鉴》(1990—2019年)中获取，并利用多因子权重分配法将以阿拉尔垦区为基本统计单元的统计数据分布到栅格单元上。

2.5 基于CA-Markov模型的土地利用模拟

马尔科夫(Markov)模型侧重于对土地利用变化数量的预测，缺乏对各土地利用类型的空间性表达，而CA模型能表达空间系统的时空动态演化过程。将Markov模型和CA模型结合(CA-Markov模型)[23]进行土地利用类型空间动态变化的模拟和预测，其既能够对土地利用的数量变化进行表达，也能对空间动态变化进行演示。

利用IDRISI软件中的CA-Markov功能，首先利用1990—1999年10期土地利用数据和收集的15种影响因素数据，模拟阿拉尔垦区2010年的土地利用变化情况并与2010年真实土地利用情况进行比较。对模拟图像进行精度检验，以验证CA-Markov模型模拟的可靠性，使用1990—2019年30期土地利用数据和相关影响因素数据模拟出2050年和2080年的土地利用变化情况。

3 结果与分析

3.1 阿拉尔垦区土地利用变化分析

3.1.1 面积变化分析

对阿拉尔垦区1990—2019年间的最大值合成影像进行土地利用分类并统计各土地利用类型面积，得到垦区各年份土地利用变化情况及其面积变化情况(图6和图7)。分类后结果的分辨率为30 m×30 m，从图7(a)中可看出阿拉尔垦区耕地、园地、水体和建设用地面积呈增加趋势，经统计分别增加了1 147.19 km2，674.19 km2，107.58 km2和36.54 km2，林草地和未利用地面积呈减少趋势，经统计分别减少了33.67 km2和1 931.83 km2，各土地利用类型在1990—1995年间年均变化率最大(图7(b))，说明此期间垦区内土地利用变化最为剧烈。

(a) 各土地利用类型面积(b) 各土地利用类型的年均变化率

3.1.2 类型转化分析

1990—2019年阿拉尔垦区各土地利用类型的转化情况如表3所示。阿拉尔垦区耕地面积增加1 516.19 km2，减少368.99 km2，净变化量为1 147.2 km2，耕地面积的净增加主要是由未利用地转化而来； 林草地面积减少了39.63 km2，主要被转化成耕地； 园地由未利用地转化了1 160.44 km2； 水体面积增加主要由未利用地转化； 建设用地面积的增加主要是由耕地转化而来； 未利用地减少了1 934.48 km2，主要被开垦为园地和耕地。

表3 1990—2019年阿拉尔垦区各土地利用类型的变化转换矩阵

不同时间段各土地利用类型变化特征不同(图8)，1990—2000年，除未利用地外，各土地利用类型面积净变化量均增加； 2000—2010年，林草地、园地和未利用地面积净变化量均减少，分别为-14.37 km2，-272.63 km2和-536.23 km2； 2010—2019年，林草地、水体和未利用地面积净变化量为减少； 1990—2019年，耕地、园地和未利用地净变化量最明显，分别为1 147.19 km2，674.19 km2和-1 931.84 km2。

(a) 1990—2000年(b) 2000—2010年

(c) 2010—2019年(d) 1990—2019年

图8 阿拉尔垦区在不同时间段的土地利用变化特征

3.1.3 空间动态变化分析

阿拉尔垦区不同时间段土地利用空间动态变化情况如图9所示。1990—2000年，未利用地发生了大面积开垦，主要分布在阿拉尔垦区中部； 2000—2010年，耕地、林草地和园地3种之间的相互转化主要分布在塔里木河沿岸区域； 2010—2019年，未利用地的开垦面积相对较少，主要分布在垦区东北和东南部； 1990—2019年，塔里木河沿岸区域主要发生的是耕地、林草地和园地之间的转化，沿岸区域之外的地方主要发生的是未利用地的开垦。

(a) 1990—2000年(b) 2000—2010年

(c) 2010—2019年(d) 1990—2019年

3.2 阿拉尔垦区2010年土地利用模拟和精度检验

将阿拉尔垦区2010年实际情况和模拟后结果进行对比(图10)，可以看出大部分区域一致，只有一些地区存在细微差异，如中心区域的建设用地、东北和东南部的林草地。

(a) 2010年实际情况(b) 2010年模拟情况

图10 2010年阿拉尔垦区土地利用的实际及模拟情况对比

通过建立混淆矩阵，进一步对2010年模拟后结果进行精度评价(表4)。表中对角线位置加黑为各类别正确分类的像素个数。通过表4可以看出，模拟精度的偏差主要是受到2005年城市规划政策新建工业园区的影响，建设用地面积增加，导致园地和建设用地精度较低。但整体看来，阿拉尔垦区的总体分类精度为89.82 %，Kappa系数为0.85，表明CA-Markov模型能够较好地模拟阿拉尔垦区土地利用变化情况，可对2050年和2080年研究区展开情景模拟。

表4 阿拉尔垦区2010年土地利用模拟的精度评价

3.3 阿拉尔垦区2050年和2080年的土地利用模拟

对比研究区2050年和2080年的土地利用模拟情况(图11)和各土地利用类型的面积变化情况(表5)，到2080年垦区未利用地仅占2.95%，耕地和园地大幅增加，共占垦区的75.41%，且空间分布较为均匀，到2080年建设用地面积迅速扩张到207.76 km2。这是由于阿拉尔垦区还受到总人口、农业生产总值、棉花价格、GDP、全社会固定资产投资、第一产业、年均气温和年均降水量等驱动因子的影响，导致模拟结果出现很大偏差。建议提前制定合理土地利用政策及生态修复措施，以提高各土地利用类型的资源利用率，促进土地资源均衡利用，开展生态环境治理，实现区域生态环境可持续发展。

(a) 2050年(b) 2080年

图11 阿拉尔垦区2050年和2080年的土地利用模拟情况

表5 阿拉尔垦区2050年和2080年各土地利用类型的面积变化情况

4 讨论

4.1 阿拉尔垦区土地利用类型面积变化的突变情况

探讨连续长时间序列土地利用类型面积变化的突变情况有利于区域土地信息监测、环境保护及制定可持续的土地利用政策[24]。采用累积距平法[25]分析连续长时间序列中的面积突变情况，从图12中可以看出耕地、园地、建设用地和未利用地在2005年出现明显转折，多数研究的突变情况则是由自然洪灾[26]和极端天气[27]影响所致，而本研究区的突变情况主要是由人为的政策干扰所致，这是由于阿拉尔垦区的政策导向所致，2004年推出西部大开发政策，新疆棉花产业扩大，农产品价格上涨，从2005年开始不断开荒扩地，耕地、园地和建设用地面积急剧增加，未利用地面积急剧减少。

(a) 耕地 (b) 园地

(c) 建设用地(d) 未利用地

4.2 土地利用变化的驱动因素

对土地利用变化的驱动因素进行深入研究，有利于土地信息监测、合理利用并为土地资源可持续发展提供科学依据[28]。本研究是以西北典型的农垦区域为对象展开的土地利用变化分析，而多数研究则是以东部沿海经济较为发达的城市区域[29]为对象展开的土地利用变化分析，驱动土地利用变化的因素多为工业化[30]和城市化[31]带来的经济发展。本文利用高程、坡度等15种影响因素数据和土地利用数据，运用CA-Markov模型进行土地利用模拟，能够显示出各影响因素的驱动力大小，最终发现本文总人口、农业生产总值和棉花价格为主要驱动因素。这是由于随着人口增加，机械化作业水平提高，农业生产总值增加，棉花价格上涨，大量的土地资源被开发，导致研究区土地利用发生剧烈变化。

4.3 不同土地利用模拟模型的比较

本文以CA-Markov模型和多智能体ABM(agent-based model)模型[32]这2种最常用的土地利用模拟方法进行比较讨论，分析本文所用模型的优势。CA-Markov模型不仅能够通过简单的转化规则模拟出复杂的土地利用类型变化模型，为更高级的土地利用演化提供良好的过程分析环境，而且能与GIS和RS数据嵌合，提高GIS分析复杂自然现象和长时序动态建模的能力。ABM模型能够通过微观智能体及其与地理空间环境相互作用的方式进行模拟，从而反映土地利用变化过程中复杂的人地关系，主要呈现研究区的人为因素影响[33]，缺乏对自然环境的考量。因此，CA-Markov模型是通过局部简单转换规则改变元胞状态，从宏观模拟土地利用变化结果，而ABM模型是通过智能体Agent交互作用，从微观模拟土地利用变化结果。这2种模型各具优势，但李少英等[34]分析近30 a相关土地利用与模拟的研究提出，CA-Markov模型结合了CA和Markov模型的优势，能够更精确地模拟土地利用变化情况。

5 结论

本文将每年各月份影像进行多时相合成后，运用支持向量机分类方法获得逐年土地利用分类图，从面积变化、类型转化、空间动态变化3方面展开土地利用变化分析，采用CA-Markov模型模拟2050年和2080年土地利用变化情况，探究土地利用类型面积变化的突变情况、驱动因素和不同模拟模型的优势，得出以下结论：

1)1990—2019年，阿拉尔垦区耕地、园地、水体和建设用地的面积呈增加趋势，林草地和未利用地的面积呈减少趋势； 耕地和园地面积的增加主要是由塔里木河沿岸区域之外的未利用地转换而来。

2)CA-Markov模型的总体分类精度可达89.82%，Kappa系数可达0.85，能够较好地模拟阿拉尔垦区土地利用变化情况； 模拟结果显示，至2080年垦区东北和东南部的未利用地将被逐渐开垦，耕地、园地和建设用地的面积将大幅增加。

3)土地利用类型面积在2005年发生转折性变化，耕地、园地和建设用地的面积急剧增加，未利用地面积急剧减少，这是由于阿拉尔垦区的政策导向所致； 垦区土地利用变化的主要驱动因素为总人口、农业生产总值和棉花价格。