阿拉爾墾區近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]。遙感影像分類方法眾多，不同學者在不同研究區得出的分類精度差異較大?；诒疚牡膶嶋H情況并參考已有研究案例中運用的解譯方法[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)；自然驅動因子選?。篨8—年均氣溫(℃)和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?？偟膩碚f，各因子對耕地面積表現為正向作用，造成耕地面積轉入大于轉出，這也是阿拉爾墾區耕地面積在近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間連續時間序列遙感數據進行耕地變化分析，從而能更全面反映研究區耕地的實際變化情況。