基于地理探測器的江蘇省沿江地區土壤侵蝕定量分析

吳鎮宇, 羅夢琦, 郭紅麗, 吳 芳, 黃昱楠, 戴玉婷, 劉 霞

(1.南京林業大學 南方現代林業協同創新中心 江蘇省水土保持與生態修復重點實驗室， 江蘇 南京210037; 2.江蘇省水文水資源勘測， 江蘇 南京 210029)

長江經濟帶發展是國家《促進中部地區崛起“十三五”規劃》中“三大戰略”之一。水土資源是制約長江經濟帶經濟社會發展的主要因素之一，根據水利部發布的2019年中國水土保持公報[1]，長江經濟帶水土流失面積3.95×105km2，占區域土地總面積19.18%。在“以生態優先、綠色發展為引領，依托長江黃金水道，推動長江上中下游地區協調發展和沿江地區高質量發展”總體要求下[2]，區域資源環境承載力、水土流失綜合整治以及生態環境保護任務更加緊迫。研究區域土壤侵蝕空間分布特征、識別土壤侵蝕的主要影響因素對指導區域水土流失綜合防治具有重要的意義。任樂萌、王慧琴[3-4]分別采用地理加權回歸模型、灰色關聯度判定、相關性分析等方法開展丹江口庫區、華北北部區域土壤侵蝕驅動因子識別及影響程度分析。何琪琳等[5]采用聚類分析與回歸分析進行東北黑土區坡耕地土壤侵蝕影響因素識別。傳統的回歸分析、相關系數分析等方法可以進行土壤侵蝕驅動力因子識別及影響評價，但難以定量分析多個變量間的耦合作用及對土壤侵蝕空間分異的影響[6-7]。地理探測器是一種可探測變量空間分異性及其驅動力的統計學方法，既能探測單個自變量對因變量的影響程度，也能判斷兩兩因子交互作用的強弱與類型，并進行風險區域識別[8]。目前，地理探測器在喀斯特地貌類型區[9-11]、陜北黃土丘陵溝壑區已有探索[8]，但在江蘇沿江地區未見研究與實踐。江蘇沿江地區是長江經濟帶的重要組成部分，區域經濟發達，開發利用程度高，人口密度大，人為活動強烈，經濟社會發展需求與水土資源承載能力的矛盾比較突出。為此，本文以江蘇沿江地區為研究區域，基于地理探測器對區域土壤侵蝕狀況與驅動力的多因子耦合影響進行定量歸因，為區域生態環境保護及社會經濟可持續發展提供理論依據。

1 材料與方法

1.1 研究區概況

江蘇沿江地區位于江蘇省南部(118°12′—122°0′E，30°45′—33°14′N)，包括南京市、鎮江市、常州市、無錫市、蘇州市、揚州市、泰州市和南通市8市59個縣市區，土地面積5.10×104km2，占江蘇省土地總面積的47.57%。研究區地勢平緩，主要地貌類型為平原、丘陵山區、河湖水域，高程一般在200 m以下，分布有老山山脈、寧鎮山脈、茅山山脈和宜溧山脈，主要河流水系有長江與太湖水系；屬亞熱帶季風氣候，多年平均降雨量1 147.19 mm；以棕壤土、潮土、水稻土、褐土等為主；地帶性植被為亞熱帶常綠闊葉林，主要樹種香樟(Cinnamomumbodinieri)、馬尾松(Pinusmassoniana)、杉木(Cunninghamialanceolata)、水杉(Metasequoiaglyptostroboides)、毛竹(Phyllostachysheterocycla)等。研究區人口密度983人/km2，區域生產總值7.30×1012元，占江蘇省的78.83%。

1.2 研究方法

1.2.1 數據源與處理 基于2 m分辨率GF-1與ZY-3遙感影像(2017年3—11月，2018年1—6月)，結合野外調查，通過人機交互解譯獲取土地利用數據。采用250 m空間分辨率MODIS植被指數產品計算標準植被指數(NDVI)與植被覆蓋度(FVC)。采用5 m分辨率數字高程模型(DEM)提取坡度、坡長。采用1981—2017年研究區188個站點逐日降雨量數據計算多年平均降雨量及降雨侵蝕力。土壤可蝕性因子數據來源于第一次全國水利普查成果[9]。

1.2.2 中國土壤流失方程 采用中國土壤流失方程(CSLE)計算土壤侵蝕模數，計算公式為：

A=RKLSBET

(1)

式中：A為土壤流失量〔t/(hm2·a)〕;R為降雨侵蝕力因子〔(MJ·mm)/(hm2·h·a)〕;K為土壤可蝕性因子〔(t·hm2·h)/(hm2·MJ·mm)〕;L為坡長因子(無量綱);S為坡度因子(無量綱);B為植被覆蓋與生物管理措施因子(無量綱);E為工程措施因子(無量綱)；T為耕作措施因子(無量綱)。模型各因子計算方法參照文獻《中國水力侵蝕抽樣調查》[10]。

1.2.3 地理探測器原理及應用 地理探測器法是通過探測要素空間異質性來揭示其驅動力的空間統計學方法，包括4個模塊[7]，采用解釋力(q值)度量自變量對因變量空間分異性的影響程度。本文以地貌形態類型、坡度、植被覆蓋度、多年平均降雨量、土地利用類型為自變量，土壤侵蝕模數為因變量，采用因子探測器、交互作用探測器對土壤侵蝕影響因素的定量歸因進行探索，采用風險探測器識別土壤侵蝕風險區域。

根據地理探測器要求，需要對土壤侵蝕影響因子數據進行離散化處理[7]。結合研究區特征，將地貌形態類型劃分為平原區、丘陵山區與水域3類；坡度分為<2°，2°～5°，5°～10°，10°～15°，15°～20°，20°～25°，25°～30°，30°～35°，>35°共9個等級；植被覆蓋度以0.1為步長，將區間范圍0～1按等距法分為10個等級；多年平均降雨量以50 mm為步長，將區間范圍993.04～1 452.55 mm按等距法分為10個等級；土地利用類型劃分水田、旱地、果園、茶園、林地、草地、城鎮村建設用地、生產建設活動擾動、交通運輸用地、水域及水利設施用地10個類型。通過對數據重采樣及坐標系轉換，在研究區范圍內采用行列數均為300的漁網均勻布設17 492個樣點，提取各影響因子基礎信息，作為地理探測器運行數據。各因子如附圖1—5所示。

2 結果與分析

2.1 區域水土流失狀況

2.1.1 水土流失空間分異性 研究區水土流失面積823.93 km2，占土地總面積1.62%，空間分異性明顯，主要分布在區域西南部，以南京、鎮江、常州、無錫為主(附圖6)。侵蝕強度以輕度為主，占水土流失面積的77.1%；中度及其以上侵蝕強度面積占比較小。

2.1.2 不同地貌形態類型下水土流失特征 研究區平原區占總面積的70.2%；丘陵山區占17.96%，主要分布于區域西南部；水域占11.84%。水土流失面積丘陵山區占91.73%；平原區占8.27%。

2.1.3 不同土地利用類型下水土流失特征 研究區水田占總面積的33.77%，水域及水利設施用地占27.53%，城鎮村建設用地占18.43%，林地與交通運輸用地分別占8.76%，7.06%，果園、茶園、旱地等僅占4.45%。水土流失主要分布于林地、生產建設活動擾動、果園等土地利用類型，分別占水土流失面積的51.1%，17.28%和13.33%，均以輕度侵蝕為主；旱地、茶園、草地、城鎮村建設用地與交通運輸用地的水土流失面積較少，共占18.29%；中度及其以上強度侵蝕集中分布于生產建設活動擾動。但從水土流失與土地利用類型面積占比看，生產建設活動擾動、果園、茶園水土流失占該土地利用類型面積的比例均在10%以上，說明這3類是區域內容易發生水土流失的土地利用類型。

2.2 土壤侵蝕影響因子顯著性研究

以地貌形態類型與土地利用類型為切入點，運行因子探測器對區域土壤侵蝕狀況進行分類分區定量歸因。考慮水域地貌形態類型無侵蝕、水田與水域及水利設施用地土地利用類型無侵蝕，本文僅圍繞平原區與丘陵山區，對8個土地利用類型開展分析。

2.2.1 不同地貌形態類型下土壤侵蝕影響因子顯著性 不同地貌形態類型下研究區土壤侵蝕的主導影響因子均為土地利用類型和坡度，但解釋力q值有差異(表1)。平原區以坡度為主導因子(0.055 1)，其次是土地利用類型(0.013 0)，由于平原區土壤侵蝕的空間分異性較小，各影響因子的解釋力均較小；丘陵山區則以土地利用類型為主導因子(0.112 6)，其次是坡度(0.036 8)。在丘陵山區，除坡度因子外，其余因子解釋力均高于平原區；可能因為丘陵山區地形較為復雜，導致植被覆蓋、土地利用的空間差異性增大，多重因素組合對土壤侵蝕產生影響，增加了土地利用、植被覆蓋、降雨量同土壤侵蝕模數空間變異的相似性。

表1 不同地貌形態類型影響因子顯著性

2.2.2 不同土地利用類型下土壤侵蝕影響因子顯著性 不同土地利用類型下研究區土壤侵蝕的主導影響因子及解釋力q值有顯著差異(表2)。茶園以植被覆蓋度為主導因子，除茶園外的其他類型土地均以坡度為主導因子。從不同影響因子來看，坡度對果園、生產建設活動擾動的影響最大，解釋力分別為0.607 1和0.590 4，其次為草地和城鎮村建設用地，茶園和林地解釋力在0.1以下；地貌形態類型對草地、生產建設活動擾動的影響最大，解釋力分別為0.131 2和0.109 2，其余土地利用類型解釋力均在0.1以下；降雨量對果園的影響最大，為0.167，其余土地利用類型解釋力均在0.1以下。林地各影響因子解釋力均在0.1以下，與各因子的空間相似性都比較低有關，可能是因為其層次結構比較復雜，受林下灌草和枯落物影響較大。

2.3 土壤侵蝕影響因子交互作用研究

2.3.1 不同地貌形態類型下土壤侵蝕影響因子交互作用 采用交互探測器可知，不同地貌形態類型下，土壤侵蝕影響因子兩兩交互作用均呈現非線性增強的趨勢，以坡度與土地利用類型的協同作用為首位，遠高于其他交互作用(表3)。由表3可知，平原區主導交互作用解釋力(0.700 6)明顯高于丘陵山區主導交互作用(0.484 9)，說明平原區更易受坡度和土地利用類型的影響，與平原區水土流失主要位于生產建設活動擾動有關。平原區土壤侵蝕的交互作用前3位以坡度為核心疊加其他因子，表明在地勢相對平坦的區域，坡度是土壤侵蝕的主導影響因子。丘陵山區土壤侵蝕的交互作用前3三位以土地利用類型為核心疊加其他因子，可能因為丘陵山區地形復雜，土地利用類型的空間差異性增大，與坡度、降雨量、植被覆蓋因子疊加增加了與土壤侵蝕空間變異的相似性。

表2 不同土地利用類型影響因子顯著性

表3 不同地貌形態類型土壤侵蝕影響因子交互作用及其解釋力(q)

2.3.2 不同土地利用類型下土壤侵蝕影響因子交互作用 不同土地利用類型下，土壤侵蝕影響因子交互作用呈現非線性增強和雙因子增強趨勢，排在前三位的交互作用詳見表4。由表4可知，城鎮村建設用地、果園、生產建設活動擾動、交通運輸用地、旱地、林地的交互作用排在第一位的均為坡度∩降雨量，城鎮村建設用地與果園解釋力最高，分別為0.959 4和0.766 8；茶園與草地的交互作用排在第一位的均為坡度∩植被覆蓋度，解釋力分別為0.919 1和0.725 8。除果園、茶園、林地與草地外，其他土地利用類型排在第二位的交互作用均為坡度∩地貌形態類型。果園、茶園排在前3位的交互作用解釋力較強，均在0.6以上；林地的交互作用解釋力均低于0.3，與林地各影響因子空間差異性較小有關。

2.4 土壤侵蝕高風險區域識別

采用風險探測器對研究區水土流失集中分布的果園、茶園、林地與生產建設活動擾動等土地利用類型進行高風險區域識別。果園、茶園、林地的土壤侵蝕高風險區域均位于坡度>20°，降雨量1 300～1 400 mm、植被高覆蓋度的丘陵山區。其原因如下： ①由于這些區域植被覆蓋度雖高，但林下蓋度較低，林分結構相對單一，林下灌草與枯落物少； ②由于果園、林地土壤侵蝕為坡度主導，茶園為坡度∩植被覆蓋度主導，丘陵山區坡度較大導致土壤侵蝕風險較高。生產建設活動擾動的土壤侵蝕高風險區域位于坡度25°～30°，降雨量1 300～1 400 mm的丘陵山區。

3 結 論

(1) 根據因子探測器結果，平原區土壤侵蝕以坡度為主導因子，但解釋力不高(0.055 1)；丘陵山區土壤侵蝕以土地利用類型為主導因子，解釋力為0.112 6。不同土地利用類型土壤侵蝕影響因子的層間差異顯著。茶園土壤侵蝕以植被覆蓋度為主導因子，解釋力高達0.867 1；除茶園外的其他類型土地的土壤侵蝕以坡度為主導因子，解釋力范圍為0.037 7～0.607 1。

(2) 兩兩因子交互作用均能增強對土壤侵蝕的解釋力。平原區與丘陵山區主導交互作用的因子為坡度疊加土地利用類型，解釋力分別為0.700 6和0.484 9，且顯著高于其他交互作用。不同土地利用類型的主導交互作用的因子主要為坡度疊加降雨量，對城鎮村建設用地與果園影響較大，解釋力分別為0.959 4和0.766 8；其次為坡度疊加植被覆蓋度，對茶園、草地影響較大，解釋力分別為0.919 1和0.725 8。

(3) 不同土地利用類型土壤侵蝕高風險區域均位于坡度>20°，降雨量1 300～1 400 mm，植被高覆蓋度的丘陵山區。

4 討 論

對土壤侵蝕影響因子進行定量歸因，有助于開展區域水土流失綜合防治工作。但土壤侵蝕是由多因子耦合影響的。目前研究土壤侵蝕影響因素主要采用多元回歸分析、相關分析、主成分分析、灰色關聯度等方法。上述方法對數據有限制，且僅能分析自變量對因變量的影響，難以反映其耦合影響和空間分異情況。本文使用的地理探測器對自變量與因變量無需假設與限定，可探測數據型數據與定性數據。該方法操作簡便，結果包括單因子影響、兩因子交互作用及風險區域分析。目前已在西部喀斯特地貌類型區、黃土丘陵溝壑區有過嘗試，在東部江蘇沿江地區進行探索十分必要。本文僅從區域層面分析土壤侵蝕影響因子，但江蘇沿江地區的平原沙土區地貌以及生產建設項目密集的特征，同時也是該區域土壤侵蝕關注的重點，可作為下一步的研究方向，結合地理探測器進行研究和探索。