盛 艷，秦富倉，劉林甫

(內蒙古農業(yè)大學沙漠治理學院，內蒙古呼和浩特 010019)

耕地對于人類的生存和發(fā)展是必不可少的資源，也是一種重要的土地利用形式。耕地資源對我國的糧食安全問題、經濟發(fā)展持續(xù)性和生態(tài)平衡起重要作用[1-3]。影響耕地的數量與質量變化的因素是自然、社會、經濟、技術等各種要素綜合。由于人口的增加、經濟迅速發(fā)展和城市化進程的加快，對建設用地的需求量不斷增加，大量的耕地轉化為建設用地，導致耕地資源數量在不斷減少，隨之生態(tài)環(huán)境也在惡化，使得耕地保護的責任與壓力不斷加大[4-5]。當前的首要任務是保護耕地，耕地演變及空間格局變化的研究對耕地資源的合理持續(xù)利用、優(yōu)化農業(yè)生產結構和實現區(qū)域農業(yè)的可持續(xù)發(fā)展具有重要意義。

關于耕地動態(tài)變化及格局演變等方面的研究我國學者已經從宏觀尺度、中尺度以及在一些經濟發(fā)達區(qū)域進行了大量相關研究，在經濟相對欠發(fā)達的低山丘陵區(qū)及微觀尺度上研究較少[6-8]。多數學者從宏觀層面揭示耕地演變及空間格局分布狀況,從耕地多時空尺度和空間關聯(lián)性整合方面去揭示耕地空間格局演變規(guī)律研究甚少?？臻g自相關是用來度量某一要素在空間上是否具有一定規(guī)律的空間變量，檢驗要素屬性值與其臨近空間要素上的屬性值是否具有顯著關聯(lián)性[9-11]。將GIS 技術與空間自相關方法相互結合，能夠全面地反映空間單元的位置及與其相鄰空間單元之間的關聯(lián)性。筆者以赤峰市松山區(qū)作為研究區(qū)，運用空間自相關方法對研究區(qū)耕地變化空間集聚性和異質性進行研究，將耕地演變空間聚集與時間變化過程有機結合,繪制 LISA 聚類圖,動態(tài)揭示多時空尺度下的耕地空間聚集狀態(tài),更為直接地反映空間自相關變化差異性及趨勢[12-13]，有效揭示耕地時空演變的特征和規(guī)律，對研究區(qū)合理利用保護耕地具有重要的理論價值和實踐意義。

1 資料與方法

1.1 研究區(qū)概況赤峰市松山區(qū)位于赤峰市東南部，地理坐標為117°47′～119°39′E、42°01′～42°43′N。松山區(qū)地處松遼平原、內蒙古高原的中間地帶，屬北溫帶大陸性季風氣候，年平均氣溫為5.6 ℃，年無霜期為155 d，年降水量為377.60 mm；年蒸發(fā)量為1 852 mm；年日照時數為2 799 h。松山區(qū)氣候特征為四季分明，春季干旱多風且蒸發(fā)量大；夏季雨熱同期,降水較為集中；秋季較為短促,氣溫下降較快,初霜降來臨較早；冬季較為漫長且寒冷,日照比較充足。研究區(qū)熱量、水分條件的時空分布差異顯著,熱量分布是東南部好于西北部,降水量的分布則以西北部偏多、東南部偏少。松山區(qū)全區(qū)總面積5 618 km2，宜農地占總面積的22.8%，全區(qū)轄9鎮(zhèn)5鄉(xiāng)、7個街道辦事處，2019年，全區(qū)戶籍人口608 883人，全區(qū)實現地區(qū)生產總值270.5億元，增長5%，其中，第一產業(yè)產值48.7億元，增長2.1%。目前，松山區(qū)已成為國家商品糧基地之一,全區(qū)糧食作物播種面積13.01萬hm2。

1.2 數據來源研究區(qū)遙感數據主要來源于美國國家航空航天局(NASA)陸地衛(wèi)星數據,包括1990、2000、2010年的Landsat-5 TM影像以及2020年的Landsat-8 OLI_TIRS 影像。利用ERDAS軟件對遙感影像進行處理，為了使影像中的地物特征顯現更為清晰，通過人機交互解譯,同時結合研究區(qū)的實際狀況,基于ArcGIS 10.8軟件對6種土地利用類型(耕地、林地、草地、建設用地、水域、未利用地)進行數據提取，主要依據國家現行標準《土地利用現狀分類》(GB/T 21010—2017)進行研究區(qū)土地利用分類[14]。選取了6種典型地類的訓練樣區(qū)，運用最大似然法對遙感影像進行監(jiān)督分類。結合實地抽樣調查對分類結果進行精度驗證，總體精度達到90.26%，kappa系數為0.90，可以較好地滿足研究的需求，同時收集了1990—2020年研究區(qū)的歷史統(tǒng)計數據及相關文獻資料。

1.3 研究方法

1.3.1動態(tài)度。耕地的變化速度主要利用動態(tài)度來進行定量描述，主要反映研究區(qū)耕地的變化程度。通過對研究區(qū)耕地變化速度分析，可以了解到研究區(qū)耕地變化的劇烈程度。計算公式如下：

(1)

式中,K為T 時段內的耕地利用動態(tài)度；Ub、Ua分別為研究期末和研究期初耕地的面積。

1.3.2空間自相關分析?？臻g自相關屬于地學空間統(tǒng)計方法中的一種，能夠有效揭示某一要素屬性值與其臨近空間要素上的屬性值是否存在明顯的相關性，計算某空間單元的特征值與其臨近的單元空間自相關性程度，可以揭示空間單元分布現象的特征及空間單元相互之間的作用機理[15-16]。空間自相關分為全局空間自相關和局部空間自相關。全局空間自相關能夠揭示整個區(qū)域總體的空間關聯(lián)及差異的顯著性，檢驗整個區(qū)域中各空間單元的屬性值是否存在集聚特性。全局空間自相關在空間單元屬性值空間分布上出現不平穩(wěn)時就不能反映單元之間的空間關聯(lián)性[17]。局部空間自相關能夠揭示局部區(qū)域內空間各單元之間的相關性，反映各個空間單元的屬性值在空間分布上的異質性[18]。

Moran’sⅠ指數是空間自相關性的度量指標，Moran’sⅠ指數的值在[-1,1][19]，在一定顯著性水平下，當 Moran’sⅠ指數 > 0時，表示各空間單元在空間分布上存在正的空間自相關且呈現集聚狀態(tài)分布；當Moran’sⅠ指數<0時，表示各空間單元在空間分布上存在負的空間自相關且呈現離散狀態(tài)； 當Moran’sⅠ指數=0時，表示各空間單元在空間分布上不存在空間自相關且呈隨機分布狀態(tài)[20]。其計算公式如下：

(2)

局部空間自相關的程度可以用Moran 散點圖來反映，Moran 散點圖的第一象限為H-H類型區(qū)，局部空間自相關呈現正的，主要表現為高值的空間單元臨近也為高值空間單元的空間聯(lián)系形式。第二象限為L-H類型區(qū)，局部空間自相關表現為負的，主要表現為低值的空間單元臨近為高值的空間單元的空間聯(lián)系形式。第三象限為L-L類型區(qū)，局部空間自相關表現為正的，主要表現為低值的空間單元臨近也為同樣的低值的空間單元的空間聯(lián)系形式。第四象限為H-L類型區(qū)，局部空間自相關表現為負的，主要表現為高值的空間單元被低值的空間單元所包圍的空間聯(lián)系形式。用 Moran 散點圖與 LISA 顯著性水平相結合的方法，能夠繪制Moran 顯著性水平圖，具有顯著性水平的 LISA 區(qū)域可以顯現，局部空間的異質狀況就可以測度，同時診斷局部空間集聚的“熱點”與“冷點”。

2 結果與分析

2.1 松山區(qū)耕地時空演變特征

2.1.1耕地資源時間序列演變。由圖1和圖2可知，松山區(qū)耕地面積總體變化態(tài)勢是急劇增加—增加—減少的變化過程。1990—2000年松山區(qū)新增耕地凈增速度最快的時期，耕地開墾量與凈增量都達到了峰值，共增長了34 459.02 hm2，增長率為28.05%，其中水澆地增長率為54.85%，增長最快；其次是旱地增長率為18.36%。2000—2010年，耕地總量增加幅度有所減緩，共增加了23 911.47 hm2，增長率為 15.20%，其中，水澆地的增長速度最快，增長率為33.45%，而水田則減少了36.33%，旱地增加了6.50%。2010—2020年耕地呈現減少態(tài)勢，減少了7.13%。影響耕地演變的主要因素是國家政策，在20世紀90年代初期，由于畜牧產品的市場價格的嚴重下跌，國家對農產品市場價格的保護，開荒種地比經營畜牧業(yè)能夠帶來更高收益，導致出現了大面積的開荒。2000年之后，各地大量過度開墾開荒使得生態(tài)環(huán)境日益惡化，國家和各級政府都采取了積極的應對措施，在生態(tài)脆弱區(qū)和農牧交錯區(qū)實施退耕還林還草措施，使研究區(qū)耕地增加的趨勢有所減緩，隨著經濟的快速發(fā)展及城鎮(zhèn)化水平的不斷提高，2010—2020年耕地出現了減少的態(tài)勢。

圖1 不同時期土地利用類型

圖2 1990—2020年松山區(qū)耕地面積變化量和變化率

由圖3可知，1990—2000年，耕地增長率超過50%的鄉(xiāng)鎮(zhèn)有王府鎮(zhèn)、當鋪地滿族鄉(xiāng)、太平地鎮(zhèn)、初頭朗鎮(zhèn)和安慶鎮(zhèn)，耕地增長率超過20%的鄉(xiāng)鎮(zhèn)有大廟鎮(zhèn)、夏家店鄉(xiāng)、哈拉道口鎮(zhèn)、崗子鄉(xiāng)、老府鎮(zhèn)和上官地鎮(zhèn)，而耕地呈現減少的鄉(xiāng)鎮(zhèn)有大夫營子鄉(xiāng)、穆家營子鎮(zhèn)和城子鄉(xiāng)。據文獻記載,主要源于1995年以來松山區(qū)存在持續(xù)的墾荒現象，導致耕地在不斷增長。由于不當開發(fā)利用導致土地退化及沙化嚴重而出現耕地撂荒行為，農地生產力下降,生態(tài)環(huán)境隨之惡化。2000—2010年，崗子鄉(xiāng)、夏家店鄉(xiāng)和上官地鎮(zhèn)耕地依然在增長，而穆家營子鎮(zhèn)耕地減少29.43%，主要由于穆家營子鎮(zhèn)屬于城郊區(qū)，建設用地占用了大量的優(yōu)質耕地；同時王府鎮(zhèn)、大廟鎮(zhèn)、城子鄉(xiāng)和大夫營子鄉(xiāng)耕地也在減少。2010—2020年，大夫營子鄉(xiāng)、城子鄉(xiāng)和穆家營子鎮(zhèn)耕地減少速度較快，城子鄉(xiāng)和穆家營子鎮(zhèn)減少的耕地主要轉化為建設用地，大夫營子鄉(xiāng)耕地退化現象嚴重，部分耕地退耕還林還牧。

圖3 1990—2020年各鄉(xiāng)鎮(zhèn)耕地面積變化

2.1.2各鄉(xiāng)鎮(zhèn)耕地變化速度分析。利用4期的耕地數據，運用公式(1)計算了各鄉(xiāng)鎮(zhèn)的耕地動態(tài)度(圖4)。依據計算結果，將各個鄉(xiāng)鎮(zhèn)的耕地變化速度劃分為3個級別，具體如下：①耕地動態(tài)度處于[-4.79%,0)是耕地面積減少區(qū)，主要有城子鄉(xiāng)、大廟鎮(zhèn)、穆家營子鎮(zhèn)和大夫營子鄉(xiāng)。城子鄉(xiāng)和穆家營子鎮(zhèn)耕地減少主要是位于城郊區(qū)，由于建設用地的不斷擴張而占用了大量的耕地，耕地的后備資源嚴重不足而未能及時補充,從而導致耕地減少。大廟鎮(zhèn)和大夫營子鄉(xiāng)處于山地丘陵區(qū)，部分耕地為坡耕地，極易產生水土流失，土壤肥力下降，生產率低，被棄耕撂荒及退耕，導致耕地減少。②耕地動態(tài)度處于[0，3.00%)是耕地面積緩慢增長區(qū)，主要有老府鎮(zhèn)、哈拉道口鎮(zhèn)和王府鎮(zhèn)。為了保持耕地總量動態(tài)平衡，實施生態(tài)災毀補充耕地和“占一補一”的制度，這3個鄉(xiāng)鎮(zhèn)大量開墾后備耕地資源，所以耕地有小幅增加。③耕地利用動態(tài)度處于[3.00%，4.91%]是耕地面積快速增長區(qū)，主要有7個鄉(xiāng)鎮(zhèn)，主要分布于松山區(qū)的中東部，海拔較低及地勢平坦，耕地后備資源充足，主要以傳統(tǒng)農業(yè)為主，由于設施農業(yè)的快速發(fā)展，在經濟利益驅使下，大面積的開墾導致耕地數量大幅增加。

圖4 1990—2020年各鄉(xiāng)鎮(zhèn)耕地動態(tài)度

2.2 耕地變化的空間自相關分析

2.2.1耕地變化的全局空間自相關分析。通過GeoDa軟件，計算得到松山區(qū)1990—2000、2000—2010、2010—2020年3個時間段耕地面積變化的全局Moran’sⅠ指數,通過顯著性檢驗(α=0.05)，3個時間段耕地面積變化的全局Moran’sⅠ指數分別為0.146、0.176和0.188，表明松山區(qū)耕地面積在空間上整體變化呈現顯著的空間聚集特征，耕地變化顯著的區(qū)域其周邊區(qū)域的耕地變化也呈現相似特征。由于3個時間段的Moran’sⅠ指數持續(xù)上升，表明松山區(qū)耕地變化在空間上的集聚態(tài)勢不斷加強。

2.2.2耕地面積變化的局部空間自相關分析?；贕eoDa軟件計算松山區(qū)3個時間段各鄉(xiāng)鎮(zhèn)耕地變化的局部Moran’sⅠ指數，并進行顯著性檢驗。從圖5可以看出，1990—2000年松山區(qū)各鄉(xiāng)鎮(zhèn)耕地面積變化的局部Moran’sⅠ指數分布在[-1.302 9,2.871 7]；6個鄉(xiāng)鎮(zhèn)的耕地變化呈現顯著的空間集聚性，而8個鄉(xiāng)鎮(zhèn)耕地變化呈現顯著的空間異質性，與相近區(qū)域耕地變化存在明顯差異。2000—2010 年各鄉(xiāng)鎮(zhèn)耕地變化的局部Moran’sⅠ指數分布在[-1.675 6,3.583 7]，有8個鄉(xiāng)鎮(zhèn)的耕地變化呈現顯著的集聚性，其他鄉(xiāng)鎮(zhèn)與臨近區(qū)域的耕地變化存在明顯差異；2010—2020 年各鄉(xiāng)鎮(zhèn)耕地變化的局部Moran’sⅠ指數分布在[-1.432 7,3.167 1]，有7個鄉(xiāng)鎮(zhèn)的耕地變化呈現明顯的空間集聚性。2000—2020年各鄉(xiāng)鎮(zhèn)的局部Moran’sⅠ指數的正值所占比重在增加，表明各鄉(xiāng)鎮(zhèn)耕地變化在空間上集聚性在增強。各鄉(xiāng)鎮(zhèn)耕地演變格局呈現顯著的空間分異性，1990—2000年松山區(qū)東部各鄉(xiāng)鎮(zhèn)耕地變化具有較明顯的集聚性，2000—2020 年聚集性逐漸轉移到松山區(qū)的中南部地區(qū)。

圖5 不同時間段各鄉(xiāng)鎮(zhèn)耕地變化的局部Moran’s Ⅰ指數空間分布

為反映各個鄉(xiāng)鎮(zhèn)耕地演變的空間自相關程度，通過Moran散點圖，繪制了各鄉(xiāng)鎮(zhèn)耕地變化的LISA聚類圖(圖6)，通過LISA聚類圖，能直觀地揭示各個鄉(xiāng)鎮(zhèn)單元的集聚模式。計算Anselin Local Moran’sⅠ指數，通過顯著性檢驗，表明研究區(qū)1990—2000年耕地變化空間集聚性不明顯，而2000—2020年耕地變化在空間分布上集聚性較為顯著，表現為正的空間自相關性。

圖6 不同時間段各鄉(xiāng)鎮(zhèn)耕地變化的LISA聚類圖

(1)第一象限為高值聚集區(qū)(H-H類型區(qū))，在此類型區(qū)的各鄉(xiāng)鎮(zhèn)耕地變化的空間差異不顯著，各鄉(xiāng)鎮(zhèn)的耕地變化呈現正的空間自相關，表明各鄉(xiāng)鎮(zhèn)與周邊鄉(xiāng)鎮(zhèn)的耕地變化幅度都較大。1990—2000年安慶鎮(zhèn)、哈拉道口鎮(zhèn)、初頭朗鎮(zhèn)和太平地鎮(zhèn)耕地變化屬于高值聚集區(qū)；2000—2010年上官地鎮(zhèn)、太平地鎮(zhèn)、當鋪地滿族鄉(xiāng)、安慶鎮(zhèn)和夏家店鄉(xiāng)耕地變化屬于高值聚集區(qū)；2010—2020年當鋪地滿族鄉(xiāng)、哈拉道口鎮(zhèn)、安慶鎮(zhèn)、太平地鎮(zhèn)和夏家店鄉(xiāng)屬于高值聚集區(qū)。以上鄉(xiāng)鎮(zhèn)主要集中分布在松山區(qū)的東中部，總體趨勢表現為耕地增加較快地區(qū)從松山區(qū)東北部轉移到東中部地區(qū)。

(2)第二象限為低值離群點區(qū)(L-H類型區(qū))，在此類型區(qū)的各鄉(xiāng)鎮(zhèn)耕地變化空間差異顯著，各鄉(xiāng)鎮(zhèn)的耕地變化表現為負的空間自相關，各鄉(xiāng)鎮(zhèn)耕地變化比周邊鄉(xiāng)鎮(zhèn)耕地變化慢，出現了局部的異質“冷點”。1990—2000年夏家店鄉(xiāng)、大夫營子鄉(xiāng)、大廟鎮(zhèn)和老府鎮(zhèn)耕地變化屬于低值離群點區(qū)，這些鄉(xiāng)鎮(zhèn)主要分布于松山區(qū)的西北部。2000—2020 年耕地變化分布于低值離群點區(qū)的鄉(xiāng)鎮(zhèn)有大夫營子鄉(xiāng)、哈拉道口鎮(zhèn)和大廟鎮(zhèn)。

(3)第三象限為低值聚集區(qū)(L-L類型區(qū))，在此類型區(qū)的各鄉(xiāng)鎮(zhèn)耕地出現減少趨勢，分布于此類型區(qū)的鄉(xiāng)鎮(zhèn)與其周邊鄉(xiāng)鎮(zhèn)相比較，耕地面積變化相對緩慢，耕地變化的局部空間差異不顯著，表現為正的空間自相關性。1990—2000年穆家營子鎮(zhèn)和城子鄉(xiāng)分布于此類型區(qū)，2000—2020 年穆家營子鎮(zhèn)、初頭朗鎮(zhèn)、城子鄉(xiāng)和王府鎮(zhèn)分布于此類型區(qū)。低值聚集區(qū)逐漸由松山區(qū)南部向北部擴展。

(4)第四象限為高值離群點區(qū)(H-L類型區(qū))，在此類型區(qū)的各鄉(xiāng)鎮(zhèn)耕地面積變化局部空間差異性較為顯著，表現為負的空間自相關，各鄉(xiāng)鎮(zhèn)耕地變化比周圍鄉(xiāng)鎮(zhèn)耕地變化幅度大，形成局部的異質“熱點”。1990—2000 年當鋪地滿族鄉(xiāng)、崗子鄉(xiāng)、上官地鎮(zhèn)和王府鎮(zhèn)屬于高值離群點區(qū)；2000—2010 年崗子鄉(xiāng)、老府鎮(zhèn)和初頭朗鎮(zhèn)屬于高值離群點區(qū)；2010—2020 年屬于高值離群點區(qū)的鄉(xiāng)鎮(zhèn)有崗子鄉(xiāng)、上官地鎮(zhèn)和老府鎮(zhèn)，以上鄉(xiāng)鎮(zhèn)主要零星分散于松山區(qū)的北部和西部。

3 結論與討論

(1)松山區(qū)耕地面積經歷了急劇增加—增加—減少的態(tài)勢,城郊地區(qū)的耕地呈現減少趨勢，耕地演變的空間差異性顯著。依據耕地動態(tài)度對14個鄉(xiāng)鎮(zhèn)的耕地變化速度進行分級，分為耕地面積減少區(qū)、耕地面積緩慢增長區(qū)、耕地面積快速增長區(qū)，揭示耕地快速增長的動因為經濟利益的驅使。

(2)松山區(qū)耕地演變空間分布存在顯著正相關,耕地變化明顯的區(qū)域其相鄰區(qū)域也呈現相同特征，空間集聚態(tài)勢呈現“弱(1990—2000 年)—強(2000—2010 年)—強(2010—2020 年)”的態(tài)勢；耕地演變低值或周圍高值的鄉(xiāng)鎮(zhèn)表現出明顯集聚態(tài)勢,局部空間異質性有所增強；各鄉(xiāng)鎮(zhèn)的LISA聚類圖直接反映各個鄉(xiāng)鎮(zhèn)耕地面積變化的集聚模式，“L-H”集聚相對穩(wěn)定，部分“H-L”集聚逐漸向“L-L”集聚轉變；松山區(qū)耕地面積變化在空間分布上的集聚態(tài)勢在不斷加強，且積聚態(tài)勢由東部各鄉(xiāng)鎮(zhèn)轉移到中南部地區(qū)。

(3)耕地演變驅動機制較為錯綜復雜，既有自然因素的影響，也受人類活動干擾，需要進一步開展耕地驅動機制研究，尤其是不同空間尺度內的耕地演變，其驅動機制是否具有空間尺度效應，探索多空間尺度內耕地利用管理模式。