阿拉尔垦区近30年耕地变化及其驱动因子分析

宋奇，冯春晖，高琪，王明玥，吴家林，彭杰

(塔里木大学植物科学学院，阿拉尔 843300)

0 引言

自1990年来中国人口数量不断增长，出现了土地资源匮乏、粮食供应不足、水资源减少、环境恶化等问题。中国作为发展中的农业大国，耕地问题不仅限制着土地的合理利用，更与农业发展和粮食安全息息相关，早有“民以食为天，食以农为源，农以土为本”[1]的说法，可知耕地在粮食安全、社会经济发展、农业可持续发展等方面占据着重要地位，是人类生存发展的物质基础[2]。随着自然因素和人为因素驱动，中国耕地开发重心向西部转移[3]。新疆维吾尔自治区土地辽阔，后备耕地资源丰富，近年来人口数量增加、机械化作业水平提高，大量后备耕地资源被开垦，在一定程度上解决了人口增长所带来的粮食安全问题，但也引发了诸如人地供需、生态需水与农业用水等矛盾。探明人口增长背景下耕地变化趋势及其驱动因子，对保障粮食安全、促进经济发展、缓解人地矛盾、实现环境保护及耕地可持续利用等方面均具有重要意义[4]。

南疆阿拉尔垦区土壤盐渍化严重，需要大量水资源淋洗土壤盐分以满足作物生长条件，而垦区水源主要来自天山山脉和昆仑山山脉的冰雪资源，淡水资源有限，耕地面积的不断增加进一步加剧了研究区生态用水不足、农业用水过度及环境恶化等问题。对研究区耕地变化信息进行监测，并分析影响耕地面积变化的驱动因子，有利于耕地合理利用与可持续发展。传统耕地调查基于实地勘探，不仅耗费大量人力、物力和财力，且不适合大面积研究，而遥感技术因其覆盖范围广、信息更新快、获取便捷、成本低、利用价值高等优势，已成为耕地变化信息监测的重要手段。目前遥感技术已在多个方面得到广泛应用，如李景刚等[5]利用遥感技术对中国近20 a间13省耕地变化信息进行监测，得出耕地面积变化主要由人类活动影响的结论；杨桂山[6]利用遥感技术监测了长江三角洲耕地近50 a的变化情况，研究表明耕地变化受到人为因素和自然因素驱动；Jia等[7]基于遥感技术分析了土地利用/覆被变化情况以及耕地面积增长对自然环境造成的影响；Wang等[8]基于遥感技术分析了土地利用/覆被变化、生态环境和人类活动之间的关系；禹丝思等[9]对比不同分类算法后证明分类后处理能够有效提高分类精度，并利用遥感技术对超大城市空间进行动态监测。以上研究为人们提供了良好的方法借鉴，但迄今针对南疆干旱区以灌溉农业为主、水土矛盾突出、种植结构以棉花和特色果林为主的绿洲农业区域的相关研究甚少。

本文针对西北地区土壤盐渍化、耕地结构复杂、种植作物多样的特点，以南疆阿拉尔垦区为典型研究区，基于最佳分类算法进行遥感解译，并分析研究区土地利用/覆被变化情况，结合统计数据和气象数据从人为因素和自然因素两方面探讨耕地面积与驱动因子以及各驱动因子之间的相互影响关系，探明驱动因子对耕地面积的直接和间接影响程度，以期为阿拉尔垦区耕地信息监测、合理利用和土地可持续发展提供科学依据，同时为西北类似干旱地区相关研究提供参考。

1 研究区概况及数据源

1.1 研究区概况

阿拉尔垦区(E80° 30′～81° 58′，N40° 22′～40° 57′，图1)位于中国新疆维吾尔自治区南疆地区，隶属新疆生产建设兵团第一师。北起天山南部，南抵塔克拉玛干沙漠北部，西至柯坪县，东到沙雅县，东西长281 km，南北宽180 km，总面积4 197.58 km2，拥有胜利、上游和多浪3大水库，水资源库存5.18×108m3[10]。紧邻阿克苏河和塔里木河，地势沿河道两侧有所抬升，西北高东南低，呈绿洲带状分布。研究区属于暖温带极端大陆性干旱荒漠气候，是典型绿洲农业地区，日照率5 869%、年均日照时数2 556.3～2 991.8 h，年均蒸发量1 876.6～2 558.9 mm，年均降水量40.1～82.5 mm，光、温、水和动植物资源丰富。农业用水以滴灌为主，有着得天独厚的种植优势，是新疆主要的优质棉产区之一[11]。

图1 研究区地理位置图Fig.1 Geographic position of study area

1.2 遥感数据

数据下载于美国地质调查局网站[12]，选取轨道号为146/32、空间分辨率为30 m、年份以整5 a为时间间隔、月份选取植物生长季的7—9月、云量低于10%的阿拉尔垦区遥感影像，共7景(图2)，分别为1990年8月2日的Landsat5、1994年8月13日的Landsat5、2000年8月5日的Landsat 7、2006年9月15日的Landsat5、2010年7月8日的Landsat5、2015年7月6日的Landsat8和2019年7月17日的Landsat8影像，其中1995年和2005年7—9月份所有阿拉尔垦区影像的云量均大于10%，故选取年份相近的1994年和2006年影像代替。在ENVI软件中对所有影像进行辐射定标、大气校正、几何纠正、裁剪、图像增强等预处理。

(a)1990年 (b)1994年 (c)2000年 (d)2006年

1.3 驱动因子数据获取及其分析方法

探明耕地面积与相关驱动因子的影响机制以及各因子间相互作用强度是耕地变化监测中急需解决的问题[13]。耕地面积变化的驱动因子主要分为人为驱动因子和自然驱动因子[14]，本文收集了《新疆生产建设兵团统计年鉴》(1990—2019年)第一师的总人口、非农业人口、国内生产总值(gross domestic product,GDP)、全社会固定资产投资、第一产业、农业生产总值和棉花价格数据代表人为驱动因子；收集《地面气象记录月报表》中1990—2019年阿拉尔垦区(区站号51730)的年均气温和年均降水量代表自然驱动因子，借助SPSS软件对收集的数据与耕地面积进行通径分析[15]，根据所得结果分析耕地面积与驱动因子、各驱动因子之间的相互影响关系。

2 研究方法

2.1 分类体系

为了区分不同地物类型，以期为阿拉尔垦区后备耕地资源定量分析提供可靠的数据支撑。本文结合阿拉尔垦区土地利用/覆被现状、《新疆生产建设兵团统计年鉴》中各师(局)土地利用情况以及新疆已有案例中的分类体系[16-17]，综合分析、对比，建立适合于阿拉尔垦区土地利用/覆被分类体系，包含耕地、林草地、园地等八类用地，详细信息如表1所示。表中，影像为2019年7月17日Landsat8 OLI B4(R),B3(G),B2(B)合成影像。以2019年为例，根据该分类体系建立相应的分类样本，每类样本点个数和分布情况如图3所示，各样本点的分布遵循均匀、全方位覆盖原则，各样本所占的像元数分别为：耕地73 940个、林草地31 591个、园地35 138个、水体31 759个、建设用地7 488个、沙地48 816个、盐碱地47 711个、其他26 503个，共302 946个分类像元，由此样本进行用地分类。

表1 阿拉尔垦区分类体系及解译标志Tab.1 Classification system and interpretation signs of Alar reclamation area

图3 每类样本点个数和分布图Fig.3 Number and distribution of sample points for each category

2.2 遥感影像分类方法

遥感影像解译通常是从影像的空间特征、色调特征以及不同植被生长差异特征着手，结合解译者对影像的认知进行地类识别[18]。遥感影像分类方法众多，不同学者在不同研究区得出的分类精度差异较大。基于本文的实际情况并参考已有研究案例中运用的解译方法[19]，经过综合分析，最终选用波谱角 (spectral angle mapper,SAM)、神经网络 (artificial neural net,ANN)、最小距离(minimum distance classification,MDC)、最大似然 (maximum likelihood classification,MLC)和支持向量机 (support vector machine,SVM)5种分类方法。全连接条件随机场(fully connected conditional random field,FC—CRF)[20]是一种具有独特优势的分类后处理方法，能够细化边界、恢复局部结构从而提高分类精度。因此本文采用FC—CRF处理方法对5种分类方法进行分类后处理，从中筛选出最佳算法来解译遥感影像。

2.3 耕地变化分析方法

1)面积变化和类型转化。耕地面积变化指从起始年份到终止年份耕地面积的数量和百分比变化；耕地类型转化指从起始年份到终止年份耕地与非耕地之间的面积转换量和空间分布特征。借助GIS的统计和空间分析功能对7期分类后影像进行耕地变化分析。

2)空间动态变化。对分类后影像进行空间分析，制作耕地空间动态变化图，用以直观展现阿拉尔垦区耕地转入和转出情况。

3 结果与分析

3.1 分类结果与精度评价

使用SAM—CRF，ANN—CRF，MDC—CRF，MLC—CRF和SVM—CRF 5种分类算法进行分类处理，以总体精度(overall accuracy,OA)和Kappa系数作为评价指标[21]，选取验证样本对分类结果进行精度验证，筛选出精度最高的算法，各分类算法的精度情况和分类结果见表2和图4。SAM—CRF算法的OA为0.65，Kappa系数为0.58，均为5种算法中最低；ANN—CRF,MDC—CRF和MLC-CRF这3种方法相较于SAM—CRF算法的精度有所提升，分类结果较好，这3种方法的总体精度为0.82～0.83，Kappa系数为0.77～0.86，SVM～CRF算法的精度最高，OA达到0.95，Kappa系数达到0.94。为验证研究区的实际精度以及FC—CRF分类后处理的有效性，选取了6个局部区域进行实地考察，对SVM和SVM—CRF分类算法结果的局部区域进行比较分析(图5)。通过对图5中局部放大区域进行实地调查，并结合Google Earth影像进行核实，证实了FC—CRF分类后处理可以将SVM分类后出现的错误斑块进行剔除并对边缘进行细化，从而提升分类结果的完整性和分类精度。

表2 2019年SAM—CRF,ANN—CRF,MDC—CRF,MLC—CRF和SVM—CRF分类结果的精度对比Tab.2 Accuracy comparison of SAM—CRF,ANN—CRF,MDC—CRF, MLC—CRF and SVM—CRF classification results in 2019

(a)SAM—CRF (b)ANN—CRF (c)MDC—CRF

(a)SVM (b)SVM—CRF

在整个研究区中随机选取119 075个测试样本，其中耕地的测试样本最多，以此保证耕地精度评价的可靠性，并建立混淆矩阵对SVM—CRF算法结果进行精度评价(表3)。可以看出各地类用户精度(user accuracy,UA)均在90%以上，其中耕地的UA值最高，为98.41%。塔里木河沿岸的河漫滩湿地容易被错分为水体，导致水体和其他两类地物类型的UA值相对较低，分别为92.81%和91.93%；制图精度(producer accuracy,PA)中也是耕地的分类精度最高，为99.77%。塔里木河沿岸的部分沙地容易被漏分，导致沙地的PA值最低，为89.49%，除沙地外，剩余地类的PA值均在90%以上。由于影像空间分辨率有限，部分地类所占像素点太少，难以分辨，容易出现错分和漏分现象，但经过SVM—CRF算法这种分类后处理避免了部分分类失误。通过总体对比分析，SVM—CRF算法表现出最好的分类效果。

表3 阿拉尔垦区用地类型精度评价①Tab.3 Evaluation of land type accuracy in Alar reclamation area

①对角线位置黑体为各类别正确分类的像素个数;纵向黑体为各类别用户精度(UA);横向黑体为各类别制图精度(PA)。

3.2 耕地变化特征

3.2.1 耕地面积变化

1990—2019年阿拉尔垦区土地利用/覆被情况、各类用地所占百分比和各时期年均变化率如图6—图8所示。30 a间耕地、林草地、园地、水体和建设用地不断增加，而沙地、盐碱地和其他不断减少。结果表明阿拉尔垦区耕地面积从233.81 km2(5.57%)增长到963.76 km2(22.96%)。不同时间段耕地面积变化情况不同，按所选影像的时间间隔将其分为6个阶段：

图8 1990—2019年阿拉尔垦区各土地利用类型的年均变化率Fig.8 Average annual change rate of land use types in Alar reclamation area from 1990 to 2019

1)1990—1994年。耕地面积的年均增长率达到6个时段中的最大值(14.95%)，从223.81 km2(5.57%)增加到373.58 km2(8.9%)。

2)1994—2000年。耕地面积的增长速率大幅降低，耕地面积从373.58 km2(8.9%)增加到516.72 km2(12.31%)，年均增长率为6.39%。

3)2000—2006年。耕地面积从516.72 km2(12.31%)增加到653.14 km2(15.56%)，共增加了136.42 km2(26.4%)。

图7 1990—2019年阿拉尔垦区各用地类型所占比例Fig.7 The proportion of each land type in Alar reclamation area from 1990 to 2019

(a)1990年 (b)1994年 (c)2000年 (d)2006年

4)2006—2010年。此时段的年均增长速率与上一时段相近，耕地面积从653.14 km2(15.56%)增加到765.64 km2(18.24%)，共增加了112.5 km2(17.22%)。

5)2010—2015年。耕地面积的年均增长率下滑为6个时段中的最小值(1.94%)，耕地面积从765.64 km2(18.24%)增加到839.94 km2(20.01%)。

6)2015—2019年。增长速率再次回升，耕地面积从839.94 km2(20.01%)增加到963.76 km2(22.96%)，年均增长率为3.69%。

3.2.2 耕地类型转化

1990—2019年，阿拉尔垦区耕地面积转入了854.2 km2，转出了124.25 km2。耕地与非耕地面积间转换频繁，如表4所示。各类用地转入成为耕地的面积由大到小排序为：其他(345.74 km2)>盐碱地(312.99 km2)>林草地(94.91 km2)>沙地(54.97 km2)>园地(32.95 km2)>水体(11.37 km2)>建设用地(1.27 km2)；耕地转出成为各类用地的面积由大到小排列为：林草地(52.49 km2)>园地(47.67 km2)>建设用地(10.59 km2)>其他(9.17 km2)>盐碱地(3.13 km2)>水体(1.08 km2)>沙地(0.12 km2)，结果表明30 a间大量未利用土地被开垦成为耕地，而原有耕地主要被转换成了林草地和园地。

表4 1990—2019年阿拉尔垦区各地类变化转换矩阵Tab.4 Transformation matrix in all parts of Alar reclamation area from 1990 to 2019 (km2)

3.2.3 耕地空间动态变化

1990—2019年间6个时段耕地空间动态变化情况如图9所示。

(a)[1990,1994)年 (b)[1994,2000)年 (c)[2000,2006)年

1)1990—1994年是耕地增长速率最快的时期，耕地转入远大于转出，主要分布在塔里木河沿岸地区。

2)1994—2000年，耕地持续向外扩张，耕地转入区域主要分布在东南部，其中由盐碱地转化成为耕地的面积最多。

3)2000—2006年，塔里木河沿岸的东北方向开始出现大量耕地转出情况，垦区西北部出现耕地转入迹象。

4)2006—2010年，垦区西北部大量耕地转入同时塔里木河沿岸地区出现连片耕地转出现象。

5)2010—2015年，塔里木河沿岸地区的耕地与非耕地间转化频繁，垦区西北部和塔里木河流域中段有大量耕地转出。

6)2015—2019年，耕地转化量为6个时期最大，分布范围最广，其中耕地转入主要来自东南部荒地开垦，转出主要发生在塔里木河流域东部。

3.3 耕地变化的驱动因子分析

3.3.1 驱动因子选择

根据阿拉尔垦区实际情况、数据可得性和代表性，经综合分析，以Y—耕地面积(km2)为因变量；人为驱动因子选取:X1—总人口(万人)、X2—非农业人口(万人)、X3—GDP(亿元)、X4—全社会固定资产投资(亿元)、X5—第一产业(亿元)、X6—农业生产总值(亿元)和X7—棉花价格(元/kg)；自然驱动因子选取：X8—年均气温(℃)和X9—年均降水量(mm)，共9个驱动因子作为自变量。

各因子单位不同，数值差距大，不能直接比较分析。为保证结果的可靠性和正确性，将各组数据进行标准化处理使各因子间存在可比性。应用无量纲标准化处理[22]，所得结果如表5所示，每组数据均值为0，方差为1，消除每组数据间的量纲关系。

表5 阿拉尔垦区1990—2019年耕地面积变化与相关指标标准化后数值Tab.5 Change of cultivated land area in Alar reclamation area from 1990 to 2019 and the value after standardization of relevant indexes

3.3.2 驱动因子分析

为探明耕地面积变化的主要驱动因子、作用路径和强度，不仅要分析各因子和耕地面积之间的相关性，还要考虑各因子之间的相互关系。进一步分析标准化后数据得到各变量间相关系数矩阵，如表6所示，表中黑体数值为相关性系数最高的值。两因子之间的数值越大表明相关性越高，结果表明各驱动因子和耕地面积以及各个因子间存在不同程度的相关性，其中非农业人口与第一产业、GDP与第一产业、GPD与全社会固定资产投资的相关系数最高，达到0.99，应对数据进一步分析，筛选出主要的驱动因子。

表6 耕地变化驱动力变量相关系数矩阵Tab.6 Variable correlation coefficient matrix of the driving force of farmland change

将因变量和自变量进行通径分析，得到以下表达式：

Y=0.57X1-1.51X3+ 1.12X4+ 0.63X6+ 0.26X7(R2=0.99，P<0.01)，

(1)

式中：Y为通径分析结果；X1为总人口(万人)；X3为GDP(亿元)；X4为全社会固定资产投资(亿元)；X6为农业生产总值(亿元)；X7为棉花价格(元·kg-1)；R2为决定系数；P为置信度。

可知X1—总人口、X3—GDP、X4—全社会固定资产投资、X6—农业生产总值和X7—棉花价格是阿拉尔垦区耕地面积变化的主要驱动因子，同时都属于人为驱动因子，表明人类活动是主要驱动因素。

基于5个主要驱动因子，进一步得到主要驱动因子和耕地面积之间的通径系数表(表7)。通过直接通径系数绝对值的比较得出5个主要驱动因子对耕地面积变化的直接影响的大小依次为：GDP(-1.51)>全社会固定资产投资(1.12)>农业生产总值(0.63)>总人口(0.57)>棉花价格(0.26)，其中GDP对耕地面积变化的直接影响最大，棉花价格最小，除GDP对耕地面积为负向作用外，其他4个因子对耕地面积为正向作用，表明GDP的增长会使得耕地面积减少。在间接影响中，全社会固定资产投资通过GDP对耕地面积变化的间接影响最大，间接通经系数为-1.51，说明GDP是耕地转出的主要驱动因子。除了各因子通过GDP对耕地面积的间接通径系数为负值外，其他因子的间接通径系数均为正值，对耕地面积有正向作用，推动耕地面积增加。各因素通过总人口、全社会固定资产投资、农业生产总值和棉花价格间接影响耕地面积增加，其中GDP通过全社会固定资产投资对耕地面积变化的间接影响最大，间接通经系数为1.12；全社会固定资产投资通过棉花价格对耕地面积变化的间接影响最小，间接通经系数为0.12。总的来说，各因子对耕地面积表现为正向作用，造成耕地面积转入大于转出，这也是阿拉尔垦区耕地面积在近30 a间持续增加的主要原因。

表7 驱动因子对耕地面积的通径系数①Tab.7 The path coefficient of driving factor to cultivated land area

4 结论

1)对SAM—CRF，ANN—CRF，MDC—CRF，MLC—CRF和SVM—CRF 5种分类算法进行精度比较，结果表明SVM—CRF算法的分类结果最佳(总体精度OA为0.95，Kappa系数为0.94)，SVM—CRF算法在原本分类精度较高的SVM算法基础上，尽可能地避免了错分和漏分现象，从而提高了分类精度，本文提出的SVM—CRF分类方法可为后续相关研究提供技术支持。

2)阿拉尔垦区耕地面积在近30 a间呈持续增加趋势，耕地面积从233.81 km2增加到963.76 km2，净变化量为729.97 km2(312.21%)。耕地与非耕地之间转换频繁，其中耕地的增加区域主要分布在垦区西北及东南部，主要以开垦未利用地为主；耕地减少的区域主要分布在塔里木河沿岸地区，耕地主要被转化成为园地和林草地。

3)阿拉尔垦区近30 a耕地面积变化的主要驱动因子为总人口、GDP、全社会固定资产投资、农业生产总值和棉花价格，除GDP对耕地面积变化有负向作用外，其他均为正向作用，5个主要驱动因子对研究区耕地面积整体表现为正向作用，这也是垦区近30 a耕地面积不断增加的主要原因。

本文是以7个典型时间断面进行30 a间的耕地变化分析，在各时间断面之间的耕地变化情况没能体现出来。因此，今后的相关研究将采用30 a间连续时间序列遥感数据进行耕地变化分析，从而能更全面反映研究区耕地的实际变化情况。