基于1990—2019年Landsat影像的干旱區(qū)綠洲土地利用變化與模擬

宋奇， 馮春暉， 馬自強， 王楠， 紀文君， 彭杰

(1.塔里木大學農學學院，阿拉爾 843300； 2.北京大學地球與空間科學學院遙感與地理信息系統(tǒng)研究所，北京 100871； 3.浙江大學環(huán)境與資源學院，杭州 310058； 4.中國農業(yè)大學土地科學與技術學院，北京 100083)

0 引言

土地利用變化是影響生態(tài)環(huán)境的重要因素之一，已成為全球變化和可持續(xù)發(fā)展研究的重要內容[1-2]。土地利用變化改變了區(qū)域的物質循環(huán)、自然景觀和各種生態(tài)過程，進而導致全球生態(tài)系統(tǒng)發(fā)生變化[3-5]，使得土地利用變化與模擬成為研究區(qū)域和全球生態(tài)系統(tǒng)變化的主要趨勢和潮流[6-8]。為降低區(qū)域土地利用變化對自然和社會環(huán)境的影響，有必要對連續(xù)長時序的區(qū)域土地利用變化進行分析并對未來土地利用變化情況進行模擬，以期為土地利用變化的調控和管理提供參考。

國內外研究者對土地利用變化的研究主要集中在土地利用時空變化特征[9]、土地利用變化機制或驅動力[10]、土地利用變化的環(huán)境影響[11]等方面，少有涉及連續(xù)長時間序列的土地內部結構詳細轉移過程與模擬。從研究區(qū)域上看，多數研究集中在東部沿海經濟發(fā)達的城市區(qū)域[12-13]，而著眼于西北干旱區(qū)經濟較為落后的典型人工綠洲區(qū)域的研究甚少。同時，城市區(qū)域的土地利用變化研究所使用的數據源多為典型時間斷面的單時相遙感影像和單一的社會統(tǒng)計數據[14]，不僅未能有效捕捉連續(xù)性信息和解譯更精細土地利用分類結果，特別是造成一些關鍵突變點或拐點信息的遺漏，而且也不能全面反映研究區(qū)地形地貌、路網和氣象情況。因此，基于連續(xù)長時間序列的多時相遙感影像和多數據源(地形地貌、路網、氣象和社會統(tǒng)計數據等)對人工綠洲區(qū)域的土地利用內部結構轉變過程及定量模擬機制研究仍需進一步加強。

阿拉爾墾區(qū)位于我國西北干旱地區(qū)，屬于典型的人工綠洲區(qū)域，墾區(qū)內荒漠區(qū)域的自然土壤表層常覆有一層鹽結皮以防止水分流失，保護土壤內部結構，過度的土地開發(fā)利用破壞了表層鹽結皮，加速了水分流失，影響了土壤水分和大氣循環(huán)系統(tǒng)。同時，阿拉爾墾區(qū)緊鄰塔克拉瑪干沙漠，自然環(huán)境惡劣，常出現(xiàn)大風、揚塵和沙暴等極端天氣，自然狀態(tài)下的鹽結皮能抵抗10級風的侵蝕，而表層鹽結皮破壞后3級風就能引起揚塵，嚴重影響空氣質量[15]。此外，植被也是研究區(qū)的天然生態(tài)屏障[16]，但自然植被對生態(tài)系統(tǒng)的維護能力大于人工植被，而阿拉爾墾區(qū)內的植被以人工植被(農作物和經濟林)為主，自然植被(林草地)為輔，且農作物多為1 a生植物，極易受到破壞。因此，有必要對墾區(qū)進行土地利用變化和模擬研究，為未來土地利用規(guī)劃提供科學依據。

本研究以西北干旱區(qū)典型人工綠洲阿拉爾墾區(qū)為例，收集1990—2019年所有可用Landsat遙感影像、地形地貌、路網、氣象和社會統(tǒng)計數據，將逐年各月份影像進行多時相合成，使用支持向量機(support vector machine,SVM)分類方法，輔以2019年高分二號遙感影像及野外勘察對多時相合成后影像進行解譯分類，獲得逐年土地利用分類圖，從面積變化、類型轉化、空間動態(tài)變化3方面進行土地利用變化分析，選用CA(cellular automaton)-Markov模型對阿拉爾墾區(qū)2010年土地利用情況進行模擬并與真實情況比較，驗證模擬模型的可靠性后，模擬2050年和2080年的土地利用情景，探究長時序土地利用類型面積的突變情況和驅動因素，進一步說明了CA-Markov模擬模型的優(yōu)勢，為典型人工綠洲區(qū)域的生態(tài)環(huán)境綜合管理和宏觀的土地資源與社會經濟統(tǒng)籌管理提供理論依據。

1 研究區(qū)概況

阿拉爾墾區(qū)位于新疆維吾爾自治區(qū)阿克蘇地區(qū)，地處塔克拉瑪干沙漠和天山山脈的交匯過渡地帶，地理坐標為E80°30′～81°58′，N40°22′～40°57′(圖1)，總面積4 105.92 km2。阿拉爾墾區(qū)地勢相對平坦，氣象多變，常有沙暴、揚塵和霜凍等災害，全年降水少蒸發(fā)大； 生態(tài)系統(tǒng)屬于荒漠—綠洲型，土壤表層常常附有一層鹽結皮，以防土壤水分蒸發(fā)； 墾區(qū)種植結構單一，大面積種植棉花。

圖1 研究區(qū)地理位置

2 材料與方法

2.1 遙感數據源及其預處理

從美國USGS網站中(https: //earthexplorer.usgs.gov/)獲取1990—2019年間阿拉爾墾區(qū)(軌道號為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 土地利用分類方法

參考《新疆生產建設兵團統(tǒng)計年鑒》(1990—2019年)中第一師各土地利用面積統(tǒng)計情況和已有的新疆土地利用分類體系[17-18]，結合阿拉爾墾區(qū)土地利用類型特點，將研究區(qū)土地利用類型分為6類。各土地利用類型的地表特征、影像特征和實地照片情況如表1所示。

表1 阿拉爾墾區(qū)分類體系及解譯標志

(續(xù)表)

以2019年為例，根據阿拉爾墾區(qū)分類體系及解譯標志情況建立各土地利用類型的分類樣本(圖4)，圖上括號內為各類樣本所占的像素個數，共計1 024 951個分類像元數。根據各分類樣本，輔以收集的相關統(tǒng)計數據、高分影像和實地調查，采用監(jiān)督分類中的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 %。整體而言，此分類方法得到的分類后結果優(yōu)異，可用于后續(xù)相關研究分析。

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

圖5 不同數據源的分類結果

表2 不同數據源的精度評價

2.3 土地利用變化分析方法

從各土地利用類型的面積變化、類型轉化和空間動態(tài)變化3方面對阿拉爾墾區(qū)土地利用變化展開分析。基于ArcGIS軟件的空間分析功能，對各時間段起始年份和終止年份的分類后影像進行分析，并制作各土地利用類型變化的轉移矩陣和空間動態(tài)分布圖。

2.4 土地利用模擬數據的收集與處理

為提高土地利用模擬精度，參考相關文獻[20-22]并結合阿拉爾墾區(qū)的土地利用現(xiàn)狀，選取高程、坡度、坡向、距高速、鐵路、公路和水系距離、氣溫、降水、人口、國內生產總值(gross domestic product,GDP)、全社會固定資產投資、第一產業(yè)、農業(yè)生產總值和棉花價格共15種影響因素進行土地利用模擬。

高程數據從中國科學院計算機網絡信息中心國際科學數據鏡像網站中獲取，在ArcGIS軟件中對高程數據進行重采樣得到坡度和坡向信息，利用Near工具得到每個柵格到高速、鐵路、公路和水系的距離； 氣溫和降水數據從中國氣象科學數據共享網中獲取，并進行反距離權重(inverse distance weight,IDW)空間插值得到氣溫和降水的空間分布情況； 人口、GDP、全社會固定資產投資、第一產業(yè)、農業(yè)生產總值和棉花價格從《新疆生產建設兵團統(tǒng)計年鑒》(1990—2019年)中獲取，并利用多因子權重分配法將以阿拉爾墾區(qū)為基本統(tǒng)計單元的統(tǒng)計數據分布到柵格單元上。

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

馬爾科夫(Markov)模型側重于對土地利用變化數量的預測，缺乏對各土地利用類型的空間性表達，而CA模型能表達空間系統(tǒng)的時空動態(tài)演化過程。將Markov模型和CA模型結合(CA-Markov模型)[23]進行土地利用類型空間動態(tài)變化的模擬和預測，其既能夠對土地利用的數量變化進行表達，也能對空間動態(tài)變化進行演示。

利用IDRISI軟件中的CA-Markov功能，首先利用1990—1999年10期土地利用數據和收集的15種影響因素數據，模擬阿拉爾墾區(qū)2010年的土地利用變化情況并與2010年真實土地利用情況進行比較。對模擬圖像進行精度檢驗，以驗證CA-Markov模型模擬的可靠性，使用1990—2019年30期土地利用數據和相關影響因素數據模擬出2050年和2080年的土地利用變化情況。

3 結果與分析

3.1 阿拉爾墾區(qū)土地利用變化分析

3.1.1 面積變化分析

對阿拉爾墾區(qū)1990—2019年間的最大值合成影像進行土地利用分類并統(tǒng)計各土地利用類型面積，得到墾區(qū)各年份土地利用變化情況及其面積變化情況(圖6和圖7)。分類后結果的分辨率為30 m×30 m，從圖7(a)中可看出阿拉爾墾區(qū)耕地、園地、水體和建設用地面積呈增加趨勢，經統(tǒng)計分別增加了1 147.19 km2，674.19 km2，107.58 km2和36.54 km2，林草地和未利用地面積呈減少趨勢，經統(tǒng)計分別減少了33.67 km2和1 931.83 km2，各土地利用類型在1990—1995年間年均變化率最大(圖7(b))，說明此期間墾區(qū)內土地利用變化最為劇烈。

(a) 各土地利用類型面積(b) 各土地利用類型的年均變化率

3.1.2 類型轉化分析

1990—2019年阿拉爾墾區(qū)各土地利用類型的轉化情況如表3所示。阿拉爾墾區(qū)耕地面積增加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年阿拉爾墾區(qū)各土地利用類型的變化轉換矩陣

不同時間段各土地利用類型變化特征不同(圖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 阿拉爾墾區(qū)在不同時間段的土地利用變化特征

3.1.3 空間動態(tài)變化分析

阿拉爾墾區(qū)不同時間段土地利用空間動態(tài)變化情況如圖9所示。1990—2000年，未利用地發(fā)生了大面積開墾，主要分布在阿拉爾墾區(qū)中部； 2000—2010年，耕地、林草地和園地3種之間的相互轉化主要分布在塔里木河沿岸區(qū)域； 2010—2019年，未利用地的開墾面積相對較少，主要分布在墾區(qū)東北和東南部； 1990—2019年，塔里木河沿岸區(qū)域主要發(fā)生的是耕地、林草地和園地之間的轉化，沿岸區(qū)域之外的地方主要發(fā)生的是未利用地的開墾。

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

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

3.2 阿拉爾墾區(qū)2010年土地利用模擬和精度檢驗

將阿拉爾墾區(qū)2010年實際情況和模擬后結果進行對比(圖10)，可以看出大部分區(qū)域一致，只有一些地區(qū)存在細微差異，如中心區(qū)域的建設用地、東北和東南部的林草地。

(a) 2010年實際情況(b) 2010年模擬情況

圖10 2010年阿拉爾墾區(qū)土地利用的實際及模擬情況對比

通過建立混淆矩陣，進一步對2010年模擬后結果進行精度評價(表4)。表中對角線位置加黑為各類別正確分類的像素個數。通過表4可以看出，模擬精度的偏差主要是受到2005年城市規(guī)劃政策新建工業(yè)園區(qū)的影響，建設用地面積增加，導致園地和建設用地精度較低。但整體看來，阿拉爾墾區(qū)的總體分類精度為89.82 %，Kappa系數為0.85，表明CA-Markov模型能夠較好地模擬阿拉爾墾區(qū)土地利用變化情況，可對2050年和2080年研究區(qū)展開情景模擬。

表4 阿拉爾墾區(qū)2010年土地利用模擬的精度評價

3.3 阿拉爾墾區(qū)2050年和2080年的土地利用模擬

對比研究區(qū)2050年和2080年的土地利用模擬情況(圖11)和各土地利用類型的面積變化情況(表5)，到2080年墾區(qū)未利用地僅占2.95%，耕地和園地大幅增加，共占墾區(qū)的75.41%，且空間分布較為均勻，到2080年建設用地面積迅速擴張到207.76 km2。這是由于阿拉爾墾區(qū)還受到總人口、農業(yè)生產總值、棉花價格、GDP、全社會固定資產投資、第一產業(yè)、年均氣溫和年均降水量等驅動因子的影響，導致模擬結果出現(xiàn)很大偏差。建議提前制定合理土地利用政策及生態(tài)修復措施，以提高各土地利用類型的資源利用率，促進土地資源均衡利用，開展生態(tài)環(huán)境治理，實現(xiàn)區(qū)域生態(tài)環(huán)境可持續(xù)發(fā)展。

(a) 2050年(b) 2080年

圖11 阿拉爾墾區(qū)2050年和2080年的土地利用模擬情況

表5 阿拉爾墾區(qū)2050年和2080年各土地利用類型的面積變化情況

4 討論

4.1 阿拉爾墾區(qū)土地利用類型面積變化的突變情況

探討連續(xù)長時間序列土地利用類型面積變化的突變情況有利于區(qū)域土地信息監(jiān)測、環(huán)境保護及制定可持續(xù)的土地利用政策[24]。采用累積距平法[25]分析連續(xù)長時間序列中的面積突變情況，從圖12中可以看出耕地、園地、建設用地和未利用地在2005年出現(xiàn)明顯轉折，多數研究的突變情況則是由自然洪災[26]和極端天氣[27]影響所致，而本研究區(qū)的突變情況主要是由人為的政策干擾所致，這是由于阿拉爾墾區(qū)的政策導向所致，2004年推出西部大開發(fā)政策，新疆棉花產業(yè)擴大，農產品價格上漲，從2005年開始不斷開荒擴地，耕地、園地和建設用地面積急劇增加，未利用地面積急劇減少。

(a) 耕地 (b) 園地

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

4.2 土地利用變化的驅動因素

對土地利用變化的驅動因素進行深入研究，有利于土地信息監(jiān)測、合理利用并為土地資源可持續(xù)發(fā)展提供科學依據[28]。本研究是以西北典型的農墾區(qū)域為對象展開的土地利用變化分析，而多數研究則是以東部沿海經濟較為發(fā)達的城市區(qū)域[29]為對象展開的土地利用變化分析，驅動土地利用變化的因素多為工業(yè)化[30]和城市化[31]帶來的經濟發(fā)展。本文利用高程、坡度等15種影響因素數據和土地利用數據，運用CA-Markov模型進行土地利用模擬，能夠顯示出各影響因素的驅動力大小，最終發(fā)現(xiàn)本文總人口、農業(yè)生產總值和棉花價格為主要驅動因素。這是由于隨著人口增加，機械化作業(yè)水平提高，農業(yè)生產總值增加，棉花價格上漲，大量的土地資源被開發(fā)，導致研究區(qū)土地利用發(fā)生劇烈變化。

4.3 不同土地利用模擬模型的比較

本文以CA-Markov模型和多智能體ABM(agent-based model)模型[32]這2種最常用的土地利用模擬方法進行比較討論，分析本文所用模型的優(yōu)勢。CA-Markov模型不僅能夠通過簡單的轉化規(guī)則模擬出復雜的土地利用類型變化模型，為更高級的土地利用演化提供良好的過程分析環(huán)境，而且能與GIS和RS數據嵌合，提高GIS分析復雜自然現(xiàn)象和長時序動態(tài)建模的能力。ABM模型能夠通過微觀智能體及其與地理空間環(huán)境相互作用的方式進行模擬，從而反映土地利用變化過程中復雜的人地關系，主要呈現(xiàn)研究區(qū)的人為因素影響[33]，缺乏對自然環(huán)境的考量。因此，CA-Markov模型是通過局部簡單轉換規(guī)則改變元胞狀態(tài)，從宏觀模擬土地利用變化結果，而ABM模型是通過智能體Agent交互作用，從微觀模擬土地利用變化結果。這2種模型各具優(yōu)勢，但李少英等[34]分析近30 a相關土地利用與模擬的研究提出，CA-Markov模型結合了CA和Markov模型的優(yōu)勢，能夠更精確地模擬土地利用變化情況。

5 結論

本文將每年各月份影像進行多時相合成后，運用支持向量機分類方法獲得逐年土地利用分類圖，從面積變化、類型轉化、空間動態(tài)變化3方面展開土地利用變化分析，采用CA-Markov模型模擬2050年和2080年土地利用變化情況，探究土地利用類型面積變化的突變情況、驅動因素和不同模擬模型的優(yōu)勢，得出以下結論：

1)1990—2019年，阿拉爾墾區(qū)耕地、園地、水體和建設用地的面積呈增加趨勢，林草地和未利用地的面積呈減少趨勢； 耕地和園地面積的增加主要是由塔里木河沿岸區(qū)域之外的未利用地轉換而來。

2)CA-Markov模型的總體分類精度可達89.82%，Kappa系數可達0.85，能夠較好地模擬阿拉爾墾區(qū)土地利用變化情況； 模擬結果顯示，至2080年墾區(qū)東北和東南部的未利用地將被逐漸開墾，耕地、園地和建設用地的面積將大幅增加。

3)土地利用類型面積在2005年發(fā)生轉折性變化，耕地、園地和建設用地的面積急劇增加，未利用地面積急劇減少，這是由于阿拉爾墾區(qū)的政策導向所致； 墾區(qū)土地利用變化的主要驅動因素為總人口、農業(yè)生產總值和棉花價格。