全球黑土區(qū)耕地動態(tài)變化及其驅(qū)動力Logistic 回歸分析

楊佳佳，呂駿超，金珊合，劉 洋，李文博

中國地質(zhì)調(diào)查局沈陽地質(zhì)調(diào)查中心，遼寧沈陽 110034

0 引言

土地利用/土地覆蓋變化（land use and land cover change，LUCC）是全球氣候變化、環(huán)境變化研究內(nèi)容的重點.地球系統(tǒng)科學(xué)的原理告訴我們，土地資源與生態(tài)環(huán)境是緊密相連的，土地資源數(shù)量與質(zhì)量的變化是生態(tài)環(huán)境變化的重要特征，生態(tài)環(huán)境信息也反映出土地資源的開發(fā)、利用狀況.只有掌握土地資源數(shù)量、質(zhì)量和可持續(xù)利用的綜合信息，才能為區(qū)域的持續(xù)發(fā)展提供有力保障［1］.另外，土地利用變化會對區(qū)域經(jīng)濟(jì)增長、收入水平以及自然資源生態(tài)環(huán)境產(chǎn)生多方面影響，同時，這些自然、社會、經(jīng)濟(jì)因素變化又會反饋于土地利用系統(tǒng)，造成土地利用變化.所以，對土地利用變化機(jī)制的深入研究有助于決策者在不同尺度水平上實施政策干預(yù)，保證土地利用的可持續(xù)［2］.

黑土土壤肥沃，生產(chǎn)力高，已成為一種寶貴的農(nóng)業(yè)經(jīng)濟(jì)資源.世界四大黑土區(qū)均為重要的農(nóng)牧業(yè)生產(chǎn)基地，保護(hù)黑土已經(jīng)是相關(guān)國家土地管理者的共識.目前，國內(nèi)外學(xué)者在黑土區(qū)開展了水土流失［3-5］、土地流轉(zhuǎn)［6］、土壤屬性［7-8］、土壤光譜［9］、有機(jī)碳變化［10-12］等諸多研究，但針對全球黑土區(qū)耕地利用變化的驅(qū)動力研究尚無.耕地變化模擬與預(yù)測對黑土耕地可持續(xù)發(fā)展及更好地理解黑土區(qū)土地利用系統(tǒng)機(jī)制有重要研究意義.

土地利用空間格局變化模擬研究一是預(yù)測土地利用變化將在什么地方發(fā)生，二是研究導(dǎo)致土地利用變化的驅(qū)動因素.土地利用這種空間變量的增減變化無法定量化描述，一般的線性回歸模型因變量必須為連續(xù)變量［13］.Logistic 回歸模型作為一種非概率性線性回歸，是針對二分類或多分類因變量進(jìn)行回歸分析時采用的分類統(tǒng)計方法，其中的自變量既可以是連續(xù)變量，也可以是分類變量.通過Logistic 回歸分析，可以計算自變量的相關(guān)權(quán)重，確定相關(guān)驅(qū)動力的影響作用，同時還可以解釋因變量的發(fā)生概率，能夠有效地解決以往土地利用變化驅(qū)動力研究中存在的局限性問題.該方法應(yīng)用于土地利用變化驅(qū)動力研究中在國內(nèi)外已有較多資料可尋［14-18］.部分學(xué)者將元胞自動機(jī)、馬爾科夫模型、神經(jīng)網(wǎng)絡(luò)等方法與Logistic 回歸模型結(jié)合來預(yù)測未來土地利用的空間格局［19-23］.本研究基于全球黑土區(qū)不同時期的土地利用數(shù)據(jù)、自然環(huán)境數(shù)據(jù)和社會經(jīng)濟(jì)數(shù)據(jù)，系統(tǒng)分析2005—2019 年間耕地變化的過程和特征，構(gòu)建全球黑土區(qū)耕地變化的Logistic 回歸模型，定量分析耕地變化的空間驅(qū)動因素，并對模型預(yù)測能力進(jìn)行檢驗，以期為理解全球黑土區(qū)耕地變化機(jī)制提供科學(xué)參考.

1 研究區(qū)概況及數(shù)據(jù)來源

1.1 全球黑土區(qū)耕地分布范圍

美國土壤分類系統(tǒng)（Soil Taxonomy，ST，1999 年）目前被聯(lián)合國糧農(nóng)組織（FAO）、聯(lián)合國教科文組織（UNESCO）以及全球45 個國家普遍采用.該分類系統(tǒng)將所有土壤分為12 個土綱，其中暗沃土綱（Mollisols）定義為草原生態(tài)土壤，以深厚、暗色表層為顯著特征，因其英文名稱來自拉丁語，這一肥沃的表層也稱為軟土層，是源于植物根系的有機(jī)質(zhì)長期加入的結(jié)果.暗沃土表層厚度不小于25 cm，該層有機(jī)質(zhì)含量大于或等于1%（C有機(jī)≥0.6%），鹽基飽和度（NH4OAc 法）大于或等于50%，潤態(tài)的明度和彩度小于或等于3［24］.其定義與中國土壤系統(tǒng)分類（Chinese Soil Taxonomy，CST，1999 年）中均腐土綱（Isohumosols）、世界土壤資源參比基礎(chǔ)（World Reference Base for Soil Resources，WRB，1998 年）中黑土（Phacozems）、黑鈣土（Chemozems）、栗鈣土（Kastanozems）的定義比較接近但又有區(qū)別.暗沃土被視為廣義的黑土，按照其定義，全球共分布四大片黑土區(qū)，包括北美洲密西西比河流域地區(qū)、中國松遼流域地區(qū)、烏克蘭與俄羅斯大平原以及阿根廷和烏拉圭的潘帕斯大草原地區(qū)［25］（圖1），其分布區(qū)內(nèi)耕地將作為本次全球黑土研究的明確范圍.

圖1 基于暗沃土綱（Mollisols）的全球黑土區(qū)耕地分布范圍Fig.1 Global cropland distribution in black soil（Mollisols）regions

1.2 全球黑土區(qū)耕地分布面積

以WGS84 為坐標(biāo)系，采用正弦投影（sinusoidal projection），依據(jù)2019 年MODIS 土地覆蓋類型三級數(shù)據(jù)500 m 產(chǎn)品（MCD12Q1）統(tǒng)計面積如下：全球黑土區(qū)耕地分布總面積為375.66×104km2，在六大洲皆有分布，但不均衡，其中歐洲分布147.32×104km2，占比39.22%；北美洲分布119.43×104km2，占比31.79%；亞洲分布74.83×104km2，占比19.92%；其他地區(qū)包括南美洲分布25.24×104km2，非洲分布7.44×104km2，大洋洲分布1.40×104km2，總計占比9.07%（圖2）.

圖2 全球各大洲黑土區(qū)分布占比Fig.2 Distribution percentage of black soil areas by continents

1.3 項目數(shù)據(jù)來源

研究采用的數(shù)據(jù)說明：MODIS 三級數(shù)據(jù)描述土地覆蓋類型全球500 m 產(chǎn)品（MCD12Q1），根據(jù)一年的Terra 和Aqua 衛(wèi)星觀測所得的數(shù)據(jù)進(jìn)行處理，信息提取采用監(jiān)督?jīng)Q策樹分類方法.選取國際地圈生物圈計劃（IGBP）作為全球植被分類方案，該土地覆被數(shù)據(jù)集中包含了17 個主要土地覆蓋類型，其中包括11 個自然植被類型，3 個土地開發(fā)和鑲嵌的地類和3 個非草木土地類型定義類.

驅(qū)動力因子來源：高程（DEM）數(shù)據(jù)來源于美國航空航天局（NASA）提供的SRTM3 數(shù)據(jù)，精度為90 m；世界人口密度來源于WorldPop 發(fā)布的地圖，是目前比較可靠的長時間序列數(shù)據(jù)；氣象因子數(shù)據(jù)來源于美國國家海洋和大氣管理局（NOAA）官方網(wǎng)站；道路、河流數(shù)據(jù)來源于Google 全球地圖.

2 研究方法

2.1 土地利用數(shù)據(jù)預(yù)處理

根據(jù)LUCC 分類體系（參考中科院分類標(biāo)準(zhǔn)），將研究區(qū)土地類型合并為6 個一級地類：耕地、草地、林地、水域、建設(shè)用地和未利用地.其中將IGBP 中常綠針葉林、常綠闊葉林、落葉針葉林、落葉闊葉林、混交林、郁閉灌叢、稀疏灌叢7 個地類合并為林地類；將多樹草原、稀樹草原、典型草原3 個地類合并為草地類；將永久濕地、農(nóng)田/植被（由農(nóng)田、森林、灌木叢和草原組成的土地，其中沒有一個組成部分占景觀的60%以上）、冰雪、裸地或低植被覆蓋地4 個地類合并為未利用地類；剩下3 種地類分別為水域（水體）、耕地、建設(shè)用地（城市和建筑）.

2.2 Logistic 回歸模型

Logistic 回歸模型是一種因變量與自變量之間為非線性關(guān)系的分類統(tǒng)計方法［13，26］.其中，因變量類型為二分類、多分類變量，自變量類型為連續(xù)和分類變量，最終以事件發(fā)生概率的形式提供結(jié)果.

以耕地為例，運用Logistic 逐步回歸方法對耕地分布的驅(qū)動因子進(jìn)行分析.此方法解決了一般線性回歸方法無法計算自變量為分類變量的問題，對于每個解釋變量生成單獨系數(shù)，并通過一定的權(quán)重運算法則被解釋為耕地類型的變化概率［13，27］.

式中：pi表示空間上每個圖斑可能轉(zhuǎn)為耕地的概率；α 為截距；xi表示相應(yīng)各候選解釋變量；β 為回歸系數(shù)；k 為解釋變量的種類.

發(fā)生比（odds）是發(fā)生頻數(shù)與不發(fā)生頻數(shù)之間的比，用來對各種自變量的Logistic 回歸系數(shù)進(jìn)行解釋，可以較直觀看出自變量對事件概率的作用，運用在本研究中可反映出各解釋變量對耕地影響程度的差異.

式中：odds 為發(fā)生比，即每個圖斑成為耕地的概率和成為其他土地類型的概率之比.由式（2）可知，每當(dāng)驅(qū)動因子增加一個單位時，發(fā)生比相應(yīng)變化exp（β）倍.exp（β）用來解釋驅(qū)動因子的作用，即為發(fā)生比率.

3 耕地變化驅(qū)動力研究

3.1 基于Logistic 模型的耕地變化驅(qū)動力分析

3.1.1 耕地數(shù)量動態(tài)變化

2005—2019 年間，全球黑土區(qū)耕地面積發(fā)生很大變化.耕地由2005 年的434.42×104km2銳減到2019年的375.66×104km2，共計減少58.77×104km2.其中2010—2015 年間耕地面積消減最大，共計減少58.19×104km2.從全球分布情況來看，亞洲黑土區(qū)耕地減少占比31.70%，北美黑土區(qū)耕地減少占比29.46%，南美黑土區(qū)耕地減少占比27.79%，其他黑土區(qū)耕地共計減少11.06%.從地類轉(zhuǎn)變情況來看，其中有76.04%的耕地轉(zhuǎn)變?yōu)槲蠢玫兀òㄓ谰脻竦亍⑥r(nóng)田/植被、冰雪、裸地或低植被覆蓋地），有19.95%的耕地轉(zhuǎn)變?yōu)椴莸兀?.94%轉(zhuǎn)變?yōu)榻ㄔO(shè)用地，1.90%轉(zhuǎn)變?yōu)榱值兀?.16%轉(zhuǎn)變?yōu)樗w.造成這種情況的原因有很多，比如耕地退化、土地災(zāi)害等自然條件變化因素，生態(tài)退耕、土地結(jié)構(gòu)調(diào)整、政策干預(yù)等人為干涉因素.

3.1.2 模型建立與抽樣過程

模型建立所需的因變量（耕地變化狀態(tài)）與自變量（影響因子）如表1 所示.將2005—2010、2010—2015、2015—2019 年3 個時期的耕地分布數(shù)據(jù)及各驅(qū)動因子用ArcGIS 軟件組織成匹配的空間數(shù)據(jù)，針對各時期內(nèi)的耕地變化情況分別建立Logistic 模型.為減弱數(shù)據(jù)空間自相關(guān)性，依據(jù)兩種因變量結(jié)果（耕地增加賦值為1、耕地減少賦值為0）在設(shè)定20 km 距離的條件下隨機(jī)抽取30%的點，再用這些點去提取相應(yīng)自變量的值，通過異常值去除、多重共線性診斷后導(dǎo)入SPSS 進(jìn)行分析.SPSS 提供了卡方值（wald）統(tǒng)計量作為回歸系數(shù)的估計，表示在模型中每個解釋變量的相對權(quán)重，用來評價每個解釋變量對事件預(yù)測的貢獻(xiàn)力［28］.需要說明的是：在Logistic 模型中設(shè)置坡度Ⅱ（2～6°）、坡度Ⅲ（6～15°）、坡度Ⅳ（15～25°）、坡度Ⅴ（＞25°）為啞變量，坡度Ⅰ（≤2°）作為它們的參照水平；設(shè)置坡向Ⅱ（陽坡）、坡向Ⅲ（半陽坡）、坡向Ⅳ（半陰坡）、坡向Ⅴ（陰坡）為啞變量，坡向Ⅰ（平坡）作為它們的參照水平；在對DEM、溫度、降水提取過程中采用了雙線性插值方法.

表1 耕地利用變化及影響因子變量Table 1 Cultivated land use changes and impact factor variables

3.1.3 結(jié)果分析

模型似然比檢驗對整體模型進(jìn)行有效性分析發(fā)現(xiàn)（表2）：各模型p 值均小于0.01，模型有效；對比各階段模型AIC 值（衡量統(tǒng)計模型擬合優(yōu)良性的一種標(biāo)準(zhǔn)），2015—2019 年階段模型構(gòu)建最好（AIC 值越低越好）.

表2 Logistic 回歸模型似然比檢驗結(jié)果Table 2 Results of logistic regression model likelihood ratio test

3 個階段模型中影響耕地變化的驅(qū)動因子類型變動較大，其中到最近道路距離、DEM、溫度分別出現(xiàn)2次，降水、坡度Ⅲ、到最近河流距離分別出現(xiàn)1 次.另外，在3 個階段模型中，人口密度變化這一變量都沒有進(jìn)入模型，可能的原因是該人口密度數(shù)據(jù)未能反映從事農(nóng)業(yè)人口變化的真實情況（表3）.

表3 耕地變化的logistic 回歸模型相關(guān)系數(shù)Table 3 Parameters of logistic regression model for cultivated land change

第1 階段（2005—2010 年）：本階段耕地變化重要的解釋變量按照重要程度（Wald 統(tǒng)計量）排列順序是到最近道路距離、降水、DEM、溫度、坡度Ⅲ.其中，到最近道路距離、DEM、溫度、降水會對耕地變化產(chǎn)生顯著的正向影響關(guān)系，坡度Ⅲ會對耕地變化產(chǎn)生顯著的負(fù)向影響關(guān)系.其中，到最近道路距離每增加1 km，其他地類轉(zhuǎn)為耕地的概率將增加0.99 倍（即1 km 的發(fā)生率是e0.011=1.011），分析原因可能是在這一階段，隨著農(nóng)業(yè)經(jīng)濟(jì)發(fā)展，道路建設(shè)投入增加成為耕地變化的主導(dǎo)原因.

第2 階段（2010—2015 年）：本階段耕地變化重要的解釋變量按照重要程度排列順序是到最近河流的距離、DEM，該兩項因子均會對耕地變化產(chǎn)生顯著的負(fù)向影響關(guān)系.其中，到最近河流距離每增加1 km，耕地轉(zhuǎn)為其他地類的概率將增加1.01 倍，分析原因可能是在這一階段，隨著區(qū)域社會、工業(yè)化與城鎮(zhèn)化的發(fā)展，河流附近土地利用的轉(zhuǎn)變成為耕地變化的主導(dǎo)原因.

第3 階段（2015—2019 年）：本階段重要的解釋變量按照重要程度排列順序是溫度、到最近道路的距離，該兩項因子均會對耕地變化產(chǎn)生顯著的正向影響關(guān)系.其中，溫度每增加1 ℃，其他地類轉(zhuǎn)為耕地的概率將增加0.84 倍，分析原因可能是在這一階段，區(qū)域社會的發(fā)展趨于穩(wěn)定，人為因素的干擾作用在減退，溫度的變化反而成為耕地變化的主導(dǎo)原因.

3.2 耕地轉(zhuǎn)入模型空間格局模擬

利用2005—2010 年的Logistic 耕地轉(zhuǎn)入模型，以2005 年為基期更新驅(qū)動因子，來檢驗在2005—2019年階段可能轉(zhuǎn)為耕地的準(zhǔn)確性，并將預(yù)測結(jié)果與2005—2019 年實際變化的耕地進(jìn)行空間比對，檢測模型的預(yù)測能力.根據(jù)2005—2010 年耕地轉(zhuǎn)入模型運行結(jié)果，耕地轉(zhuǎn)入概率P 的預(yù)測模型為：

其中，Y=ln［P/（1-P）］=-1.206＋0.011×距道路距離＋0.848×DEM＋0.024×溫度＋1.750×降水-1.311×坡度Ⅲ.

通過研究區(qū)耕地轉(zhuǎn)入概率模擬結(jié)果與2000—2019 年耕地變化分布的對比發(fā)現(xiàn)（圖3、4），2005—2019 年間呈片狀增加的耕地大部分落入耕地轉(zhuǎn)入預(yù)測概率高值區(qū)，其間在實際轉(zhuǎn)入的耕地中，有71.64%圖斑的預(yù)測概率值在0.5 以上，52.07%圖斑的預(yù)測概率值在0.6 以上.由于土地利用系統(tǒng)的復(fù)雜性，土地利用空間分布的預(yù)測準(zhǔn)確性很難達(dá)到85%以上［2，29］.可見，雖然2 個階段影響耕地變化的解釋變量不同，但從模擬的結(jié)果來看，該模型對轉(zhuǎn)入的耕地預(yù)測效果較好，說明在驅(qū)動過程保持不變或者近似的條件下，近期的耕地增加模型可以有效地預(yù)測未來一定時期內(nèi)耕地轉(zhuǎn)入空間分布概率.

圖3 2005—2019 年中國黑土區(qū)耕地轉(zhuǎn)入預(yù)測結(jié)果及概率分布Fig.3 Predicted probability distribution of cultivated land in black soil regions of China during 2005-2019

圖4 2005—2019 年北美黑土區(qū)耕地轉(zhuǎn)入預(yù)測結(jié)果及概率分布Fig.4 Predicted probability distribution of cultivated land in black soil regions of North America during 2005-2019

4 結(jié)論

1）全球黑土區(qū)3 個階段（2005—2010、2010—2015、2015—2019 年）耕地變化顯著，耕地變化特征先后出現(xiàn)小幅增加—大幅減少—小幅增加的趨勢.14 年間耕地共減少58.77×104km2，其中2010—2015 年間耕地面積消減58.19×104km2.從全球分布情況來看，亞洲黑土區(qū)耕地減少幅度最大，占比31.70%.從地類轉(zhuǎn)變情況來看，耕地轉(zhuǎn)變?yōu)槲蠢玫卣急茸畲?，約占76.04%.

2）全球黑土區(qū)耕地變化在各個階段有著不同的重要驅(qū)動因素.在2005—2010 年階段最主要的驅(qū)動因素是到最近道路距離、降水和DEM；2010—2015 年主要的驅(qū)動因素是到最近河流距離和DEM；2015—2019年主要的驅(qū)動因素是到最近道路距離和溫度.通過對耕地變化Logistic 回歸模型的檢驗，得出Logistic 回歸模型進(jìn)行耕地變化驅(qū)動力分析具有可構(gòu)建性、可操作性等優(yōu)點，能較好地揭示全球耕地變化過程的主要驅(qū)動力.同時，本研究選取的自然因子、社會因子指標(biāo)中還缺乏一些重要變量，比如土壤類型、地貌類型、國內(nèi)生產(chǎn)總值（GDP）等，在一定程度上影響了模型的精度.

3）全球黑土區(qū)耕地變化是一個相當(dāng)復(fù)雜的系統(tǒng)過程，與全球各地區(qū)的經(jīng)濟(jì)、社會、人口、技術(shù)、自然因素、政治結(jié)構(gòu)與價值觀念等直接相關(guān)，有些因素難以量化納入模型，因此為了提高耕地變化的預(yù)測概率，我們應(yīng)考慮建立驅(qū)動力的動態(tài)模型.這些驅(qū)動力應(yīng)包括政策的制定與調(diào)控、農(nóng)產(chǎn)品的價格變化等動態(tài)變化因素，這樣，代表變化可能性原因的驅(qū)動力在空間模型中就可以顯現(xiàn).因此，從這方面選取一些可以量化的指標(biāo)應(yīng)是下一步研究的重點內(nèi)容.