基于Logistic回歸模型的陜西秦巴山區林地變化驅動力分析

鄧元杰，侯孟陽，張 曉，賈 磊，李園園，姚順波*，龔直文，劉廣全

(1.西北農林科技大學經濟管理學院，西北農林科技大學資源經濟與環境管理研究中心，陜西 楊凌 712100；2.中國水利水電科學研究院，北京 100038)

森林在維持陸地生態系統平衡、保障人類基本生存條件、減緩氣候變化和實現社會經濟可持續發展等方面起著重要的基礎性作用[1]。而林地作為森林資源的載體，不僅是林業可持續發展的基礎，同時也是維護森林生態服務和環境價值的前提[2]，因此辨析林地的時空變化特征，及其驅動力因素，對了解區域生態環境變化、合理利用資源和協調人地關系具有重要科學意義。

已有學者對林地變化的研究主要是描述研究區林地時空變化特征，并從研究區的自然社會等宏觀背景出發，定性分析可能引起林地變化的原因[2-3]。但這些研究無法定量揭示引起林地變化的驅動機制，導致決策者在制定林地保護政策時缺乏針對性與有效性，進而不利于林地資源保護及區域的可持續發展。有學者在林地時空變化特征描述的基礎上，綜合自然、人文、社會、經濟等因素，構建合適的計量模型探究引起林地變化的原因[4-6]，然而線性回歸模型容易忽略引起土地利用變化的驅動力(生物、自然和社會經濟)或制度(政策)在空間上存有的高度異質性，不能處理林地變化過程中的空間因素變量，導致難以客觀準確地揭示林地變化的原因。而Logistic回歸模型不僅能考慮引起土地變化驅動力因素的空間異質性，處理土地變化過程中的空間因素變量[7]，還可以在區域土地利用/覆被動態變化的基礎上定量分析其與社會、經濟、技術、政策、自然環境等驅動因子的相互關系，進而綜合不同時空尺度區域的土地利用/覆被變化過程和驅動力影響[8]。Shahbazian等[9]基于Logistic回歸建模，探究了城市空間擴張的驅動力因素，并進一步模擬了未來城市空間擴張的方向，而李云龍等[10]則聚焦于耕地，運用Logistic回歸模型剖析了影響耕地變化的驅動力因素。此外，還有學者從研究區土地利用變化的整體格局出發，對每一種土地利用類型變化均進行了Logistic回歸建模，識別出影響各土地利用類型變化的不同驅動因素[11]。

然而，目前較少有學者將Logistic模型運用在探究林地變化的驅動因素上，雖然有學者運用Logistic回歸模型對林地變化的驅動因素做了一定的解釋[12-13]，但他們在構建模型時忽略了林地增減變化與其所對應驅動因素的空間異質性，即在構建林地變化Logistic回歸模型時，將林地增加和減少均納入Logistic回歸模型中，而未將林地增加和減少進行單獨建模。使得在探究林地變化的驅動因素時，無法有效將林地增加和減少的驅動因素單獨分離，進而不利于對林地變化規律的準確把握?；诖?，本研究聚焦于陜西秦巴山區林地變化，在分析2000—2018年林地時空變化的基礎上，對林地增加和林地減少分別構建Logistic模型，探究影響陜西秦巴山區林地變化的驅動因素，以期為該區林地保護政策制定和林地資源可持續利用提供理論依據。

1 材料與方法

1.1 研究區概況

秦巴山區是秦嶺地區和大巴山地區的簡稱，位于我國中西部地區(105°30′～110°0′E,31°40′～34°40′N)，地跨陜西、甘肅、四川、重慶、湖北、河南5省1市，其主體位于陜西南部地區，北倚秦嶺，南靠巴山，漢江穿流而過，形成了兩山夾一川的地貌特征，地形起伏大，山地面積廣，盆地小且東西長，南北短，轄寶雞、西安、漢中、安康和商洛5市38個縣，該區土地面積約為86 995 km2，占陜西省總土地面積的42.31%。

1.2 數據來源與預處理

土地利用數據來源于地理國情監測云平臺(http://www.dsac.cn/)，共包括2000年和2018年兩期(圖1)，其空間分辨率為30 m。該數據集的生產是以2000年Landsat 7 ETM+和2018年Landsat 8 OLS_TIRS遙感影像為數據源，通過人工目視解譯生成，后期經過實地驗證，誤差修正后的土地利用以及類型綜合評價精度達到93%以上[14-15]。同時，利用Google Earth高分辨率衛星歷史影像與同期Landsat系列衛星影像對數據中存在的異常區域進行了再次驗證與修正，確保了土地利用數據的客觀性與準確性。土地利用數據分類系統根據“國家資源環境遙感宏觀調查與動態研究”項目中的分類標準[16]，將陜西秦巴地區土地利用類型重新分為耕地、林地、草地、水域、建設用地和未利用地6類，其中林地由有林地、灌木林、疏林地和其他林地4種林地類型構成。

本研究基于ArcGIS平臺，對兩期土地利用數據進行空間分析，將2018年相對于2000年林地增加柵格，即2000—2018年其他地類轉變為林地的柵格賦值為1，林地保持不變的柵格賦值為0，將其作為林地增加驅動力回歸模型的因變量(圖2A)；將2018年相對于2000年林地減少柵格，即2000—2018年林地轉變為其他地類的柵格賦值為1，林地保持不變的柵格賦值為0，將其作為林地減少驅動力回歸模型的因變量(圖2B)。

參考類似研究成果[12-13]，并基于秦巴山區實際特征和數據代表性、可獲取性，從自然、社會經濟、地理區位3個方面構建了林地變化驅動力體系(表1)。各具體指標數據來源及預處理方法如下：①自然驅動力，主要包括海拔、坡度兩個地形因子，年降雨量傾向率、年均氣溫傾向率兩個氣象因子以及土壤有機質含量。其中，地形因子數據來源于美國國家航空航天局在2020年最新發布的NASA DEM(https://lpdaac.usgs.gov);坡度數據為基于DEM數據，運用ArcGIS10.5軟件中的坡度函數計算獲得；氣象因子的傾向率用以表征氣象因子變化趨勢，正值表示趨向增加或提高，負值表示趨向減少或降低，傾向率的計算方法詳見文獻[17]，氣象因子數據來源于中國氣象數據網(https://data.cma.cn/)提供的秦巴山區及周邊區域12個氣象站點2000—2018年逐年氣候資料，計算得到降雨和氣溫的傾向率；土壤有機質含量提取自聯合國糧農組織(FAO)所構建的世界土壤數據庫HWSD V1.2(http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HWSD_Data/)。②社會經濟驅動力，主要包括陜西秦巴山區2000—2018年各縣總人口密度變化和年地均GDP變化。各縣總人口和GDP數據均來源于陜西省統計年鑒(2001年，2019年)。③地理區位驅動力，主要包括到道路的距離、到鐵路的距離、到城鎮的距離、到農村居民點的距離和到河流的距離。其中，道路、鐵路和河流的數據源于國家基礎地理信息中心公布的1∶100萬全國基礎地理數據庫(https://webmap.cn)，且在此基礎上結合OpenStreetMap和最新的Google Earth影像進行修正。農村居民點和城鎮數據則從秦巴山區2018年的土地利用數據中提取得到。為保證空間數據的一致性，以上所有數據空間分辨率均為30 m，空間坐標系均為WGS_1984_UTM_Zone_49N，所有數據均使用ArcGIS 10.5軟件實現了空間化表達。

表1 2000—2018年林地變化驅動力指標

1.3 Logistic回歸模型

二分類Logistic回歸模型可以對二分類因變量和自變量(或連續變量或混合變量)進行回歸建模，設第i個案例的自變量xi=(x1、x2……xn)，二分因變量Yi取值0或1(Yi=0表示主體事件不發生，Yi=1表示主體事件發生)。即林地增加時Y=1，否則為0；林地減少時Y=1，否則為0。該模型方法是基于抽樣數據，為各自變量產生回歸系數，并通過這些系數來討論模型中因變量與自變量的關系。設P為事件發生概率，取值范圍為0～1，則1-P為該事件不發生的概率，這種概率可以用Logistic函數計算，其表達式為[18]：

P=(Y=1|x1,x2……xn)=

(1)

通常對公式(1)進行Logit變換，以轉變為線性公式：

(2)

式中:x1、x2……xn為自變量；參數β0、β1……βn為待求的回歸系數；P/(1-P)為1組事件的發生比，發生比率expβ是衡量自變量對因變量的重要指標。expβ表示自變量每增加1個單位，事件發生比的變化倍數(expβ<1，發生比減少；expβ=1，發生比不變；expβ>1，發生比增加)[18]。Waldχ2用來評價每個自變量對因變量的解釋程度[18]。對Logistic回歸方程擬合度的檢驗選用Homsmer-Lemeshow 指標(HL)，當HL指標統計顯著表示模型擬合不好；相反，當HL指標統計不顯著表示模型擬合程度較好。除此之外，對Logistic回歸模型的解釋能力檢驗，通常采用受試者工作特征曲線(ROC曲線)進行驗證。該曲線常被用于度量 Logistic 回歸模型的吻合度, ROC曲線下面積的值為 0.5～1.0, 其值越高, 表明Logistic回歸分析模型的吻合度越好, 反之越差[19]。

2 結果與分析

2.1 2000—2018年秦巴山區林地變化特征

基于秦巴山區2000年和2018年兩期土地利用數據，借助ArcGIS 10.5軟件對其林地變化進行了空間分析，再運用Excel軟件對分析結果進行統計可得，2000—2018年秦巴山區林地面積由614.84萬hm2增長至618.88萬hm2，林地變化以轉入增加為主，18年間共轉入23.22萬hm2，轉出19.23萬hm2，凈增加39 951.72 hm2。由林地變化對地形因子的響應可得(圖3)，林地變化主要集中分布在海拔≥500～1 000 m、≥1 000～1 600 m，坡度≥15°～25°、≥25°范圍。其中，在海拔<500 m和坡度<6°范圍內，林地面積減少大于增加；而在海拔≥500 m和坡度≥15°范圍內，林地面積則為增加大于減少。

林地變化結構見圖4，陜西秦巴山區耕地轉為林地是林地面積增加最主要的原因，18年間共有21.83萬hm2耕地轉為林地，占林地轉入總面積的94%。此外也有少部分的草地(9 824.94 hm2)、未利用地(1 880.73 hm2)、建設用地(1 405.71 hm2)和水域(845.91 hm2)轉為林地。在林地轉出方面，林地轉為耕地是導致林地面積減少最主要的原因，18年間共有16.95萬hm2林地轉為耕地，占林地轉出總面積的88.17%。此外分別有11 303.91、4 829.04、3 553.20 和3 049.83 hm2的林地轉為草地、未利用地、水域和建設用地。從林地變化的空間分布來看，林地增加和減少在秦巴山區各處均有分布。但林地增加主要分布在秦巴山區南部和西南部(圖2A)；而林地減少主要位于秦巴山區東北部(圖2B)。

2.2 林地變化的Logistic回歸分析

在構建林地變化驅動力回歸模型之前，對各自變量進行共線性問題診斷，結果自變量的容忍度(tolerance)為0.50～0.96，方差膨脹因子(variance inflation factor, VIF)均小于 10，表明自變量間共線性不嚴重，均可以納入Logistic回歸模型。

2.2.1 林地增加的Logistic回歸分析

在林地增加的Logistic回歸模型中，坡度用4個虛擬變量分別代表坡度Ⅱ(≥2°～6°)、坡度Ⅲ(≥6°～15°)、坡度Ⅳ(≥15°～25°)和坡度Ⅴ(≥25°)，坡度Ⅰ(<2°)作為它們的參照對象。在HL檢驗中，Sig.值為0.355，大于 0.05，統計不顯著，即模型的擬合效果很好; 并由回歸結果可知，似然比卡方檢驗的觀測值為122.51(P<0.001)，且ROC曲線下的面積(AUC)大于0.7，該模型通過了各項檢驗，可用于林地增加的驅動力分析。

顯著性水平和Waldχ2統計量結果見表2，2000—2018年，秦巴山區林地增加較為重要的自變量依次為海拔、到城鎮的距離、年均氣溫傾向率、坡度、到農村居民點的距離、土壤有機質含量、到道路的距離和到鐵路的距離。其中，海拔、到城鎮的距離、年均溫傾向率、到道路的距離和到鐵路的距離對林地增加呈顯著負向影響。具體而言，海拔每升高1 m，林地面積增加的概率將降低0.19%；到城鎮、道路和鐵路的距離每增加1 m，林地面積增加的概率將分別降低0.002%、0.010%和0.001%；年均氣溫升高將不利于林地面積增加。坡度、到農村居民點的距離以及土壤有機質含量對林地增加呈顯著正向影響。具體而言，坡度Ⅱ(≥2°～6°)、坡度Ⅲ(≥6°～15°)、坡度Ⅳ(≥15°～25°)和坡度Ⅴ(≥25°)這4個坡度范圍對林地增加的影響均為顯著，其產生的正向影響均大于坡度Ⅰ(<2°)。且根據表3中Waldχ2統計量可得，坡度Ⅴ(≥25°)的Waldχ2統計量最大，表明其對林地面積增加的正向影響相比于其他坡度范圍更大。即以坡度Ⅰ(<2°)作為參照對象，林地面積增加概率在≥25°的坡度范圍內將會升高37.5%。到農村居民點的距離每增加1m，林地面積增加的概率將會升高0.004%；土壤有機質含量每增加1%，林地面積增加的概率將會升高4.64%。

表2 2000—2018年秦巴山區林地增加驅動力因素分析

2.2.2 林地減少的Logistic回歸分析

在林地減少的Logistic回歸模型中，坡度用4個虛擬變量分別代表坡度Ⅱ、坡度Ⅲ、坡度Ⅳ、坡度Ⅴ，坡度Ⅰ作為它們的參照對象。HL檢驗中，sig.值為0.148，大于 0.05，統計不顯著，即模型的擬合效果很好; 由回歸結果可知，似然比卡方檢驗的觀測值為103.45(P<0.001)，且ROC曲線下的面積(AUC)大于0.7，該模型通過了各項檢驗，可用于林地減少的驅動力分析。

根據顯著性水平和Waldχ2統計量(表3)可得，2000—2018年，秦巴山區林地減少較為重要的自變量依次為坡度、海拔、到城鎮的距離、到道路的距離、到農村居民點的距離和年均氣溫傾向率。其中，海拔、到城鎮的距離、到道路的距離及年均溫傾向率對林地面積減少呈顯著負向影響。具體而言，海拔每升高1m，林地面積減少概率將會降低0.12%；到城鎮、道路的距離每增加1m，林地面積減少的概率將會分別降低0.010%和0.002%；年均溫升高將有利于抑制林地面積減少。坡度和到農村居民點的距離對林地面積減少呈顯著正向影響。具體而言，坡度Ⅱ(≥2°～6°)、坡度Ⅲ(≥6°～15°)、坡度Ⅳ(≥15°～25°)和坡度(≥25°)這4個坡度范圍對林地面積減少的影響均為顯著的，其產生的正向影響均大于坡度Ⅰ(<2°)。且根據表3中Waldχ2統計量可得，坡度Ⅳ(≥15°～25°)的Waldχ2統計量最大，表明其對林地面積減少的正向影響也最大。即以坡度Ⅰ(<2°)作為參照對象，林地面積減少概率將會增加153.25%。到農村居民點的距離每增加1 m，林地面積減少的概率將會升高0.003%

表3 2000—2018年秦巴山區林地減少驅動力因素分析

3 討 論

1)2000—2018年，陜西秦巴山區林地面積增減共存，但總體呈增加態勢，且主要集中分布在中高海拔和坡度>15°范圍內，這與張碧桃等[20]研究結果一致。究其原因，陜西秦巴山區受退耕還林、長江天然林保護等生態修復工程的影響，坡耕地成為退耕還林的主要對象，再加之中國經濟的快速發展，山區農民進城務工趨勢明顯，耕地撂荒嚴重，農用地成為林地增長的主要來源[21]。與此同時，道路修建、居民點建設和非法采礦等人類社會經濟活動的開展是導致林地減少的重要原因[22]，但從總體來看，林地增加面積高于林地減少面積。此外，林地面積在低海拔和坡度<6°范圍內減少大于增加，這與城市擴張、交通設施建設等活動侵占了少量林地密切相關[23]。

2)Logistic回歸模型結果顯示，海拔和坡度是影響林地變化的主導驅動因素，這與前人的研究結論相同[12-13]，然而具體到秦巴山區，隨著海拔的增加，林地增加和減少的概率均呈逐步降低態勢，這一結果與本研究相吻合。林地增加概率隨海拔升高而降低，這與生態修復工程的開展多集中于開發更容易、成本較低的中低海拔地區有關，因為該區域水熱條件充沛，工程實施可使區域林地在短時間內得到增加，實現生態修復目標[13]。同時，社會經濟活動開發難度及其成本隨海拔的升高而增加，高海拔地區林地受人為活動擾動小，因此林地減少概率隨海拔上升而下降。林地增加概率隨坡度上升而增大，此結論與Xie等[24]的研究結論一致。中國自20世紀末在全國開展的以退耕還林為主的一系列生態修復工程是引起林地隨坡度上升而增加的重要原因[25]；相反，林地隨坡度增加而減少，這可能與道路修建、森林砍伐、居民點建設和非法采礦密切相關[24]。其他兩個自然因子氣溫傾向率和土壤有機質對林地變化也產生了較顯著的影響。其中，從年均氣溫傾向率對林地變化影響來看，2000—2018年，研究區內年均氣溫整體呈升高態勢。林地增加和減少的概率均隨著年均氣溫傾向率的增加而降低，雖然溫度適度升高有助于樹木生長[26]，并一定程度上阻緩林地退化，但秦巴山區的植被生長對氣溫變化敏感度較高，年均氣溫升高反而不利于植被的生長與林地面積的擴大[4,27]，因為溫度升高所誘導的干旱可能會導致研究區域樹木生長緩慢甚至死亡[28]。土壤有機質對林地增加的影響為正，對林地減少的影響卻不顯著，這是由于土壤有機質含量較高的林地土壤肥力綜合指數較大，秦巴山區人工林工程有效地改善了土壤養分，維持了土壤的肥力，提高了土壤持水的能力[29]，從而使得林地增加的可能性隨著土壤有機質含量的增加而升高。

3)地理區位因子也是影響秦巴山區林地變化的重要因素，到城鎮的距離以及到道路的距離，均對林地變化產生負向影響，表明林地增加和減少的概率隨著到道路和城鎮的距離增加而下降，隨著秦巴山區社會經濟和城鎮化的快速發展，一方面區域人口增長和經濟發展所產生的大量建設用地擴張需求導致林地被占用[23]；另一方面，陜西秦巴山區地方政府在生態文明建設等一系列生態保護政策的引導下，越來越重視生態環境的保護和治理，追求環境與經濟的協調發展[30]，通過在城市、道路周邊植樹造林，增加了城市周圍林地覆蓋率，改善了城市的生態環境。此外，模型結果顯示，林地減少概率隨著到農村居民點距離增加而增大，這可能與農戶作為退耕還林等生態修復政策的主要參與者，其退耕地離農村居民點較近有關[31]；同時，離農村居民點越近的林地，受人類活動干擾也越大，如道路建設、毀林開荒等[24]，因此距離農村居民點越近，林地增加可能性越小。

4)通過Logistic回歸分析，本研究雖較好地揭示了引起陜西秦巴山區林地轉化的驅動機理，但林地轉化的驅動機制極其復雜，與研究區的社會、經濟、人口、自然因素、政府政策與價值觀念等密切相關[8]，且受限于數據量化方法的可操作性，導致有些因素被遺漏或難以量化而未納入模型。首先，諸如政府政策、自然保護區影響程度等因素在特定區域與時間內對林地轉化驅動巨大，但難以精確合理地量化[10]，例如促使林地增加的退耕還林等生態修復政策；其次，受困于數據獲取，導致本研究所使用的社會經濟數據只能統計到區縣一級，如能使用更精確、細化的社會經濟數據，分析結果將更加精準。