基于多元非線性空間建模的拉薩河流域溝蝕發生風險探測

李建軍，陳玉蘭，焦菊英,，陳同德,3，陳一先,4，趙文婷，趙春敬,5，尚天赦，簡金世,，曹 雪

基于多元非線性空間建模的拉薩河流域溝蝕發生風險探測

李建軍1，陳玉蘭2，焦菊英1,2※，陳同德1,3，陳一先2,4，趙文婷1，趙春敬1,5，尚天赦1，簡金世1,2，曹 雪2

（1. 西北農林科技大學水土保持研究所，黃土高原土壤侵蝕與旱地農業國家重點實驗室，楊凌 712100；2. 中國科學院水利部水土保持研究所，楊凌 712100；3. 青海民族大學政治與公共管理學院，西寧 810007；4. 西北大學城市與環境學院，西安 710127；5. 黃河水利科學研究院，鄭州 450003）

青藏高原的生態環境面臨著氣候變暖和人類活動增加的雙重壓力，增加了土壤侵蝕風險。溝蝕是土壤侵蝕最為劇烈的表現形式，為調查當地溝蝕現狀和主控因素，該研究選擇拉薩河流域作為代表，通過野外調查和遙感解譯建立2 171個樣點，并首次基于最優尺度回歸、地理探測器和兩者的組合共4種方法對15個影響溝蝕的因子及其分級/分類的重要性和溝蝕發生風險進行了探測。結果發現：1）在最優尺度回歸中，因子系數前三位分別為海拔（0.442）、土壤類型（0.168）和歸一化植被指數（0.156）；在地理探測器中，海拔（0.263）、土壤類型（0.251）和人類足跡（0.174）排在前三位。2）最優尺度回歸和地理探測器的受試者工作特征曲線下面積（Area Under Curve，AUC）值分別為0.899和0.833，兩種組合方法AUC值分別為0.866和0.848，各方法探測效果均良好，都適用于空間建模。3）拉薩河流域有9.52%～13.97%的區域有著非常高的溝蝕風險，主要集中在拉薩河下游河谷兩岸和當雄盆地等相對低海拔地區。研究結果可為青藏高原生態安全屏障建設和水土保持工作提供參考。

溝蝕；遙感；空間建模；最優尺度回歸；地理探測器；青藏高原

0 引 言

青藏高原被譽為“世界屋脊”，高海拔帶來的高寒缺氧和巨大的地形起伏使得該地區適宜人類生存的范圍十分有限，土地資源緊缺[1-2]。同時，由于海拔高、氣溫低，巖石風化程度較低，成土過程緩慢，土層薄，土壤資源十分稀缺[3]。目前，該地區正面臨著氣候變化和人類活動增加的雙重壓力，生態環境非常脆弱，有限的土地資源面臨的土壤侵蝕和土地退化風險加劇[4]。

溝蝕是一個全球性的環境問題[5]，是土壤侵蝕的主要形式之一，對侵蝕總量的貢獻很大[6]。侵蝕溝一旦產生，之后的侵蝕過程會大大加速。溝蝕不僅會造成土壤資源的流失，還會破壞農田、草場等寶貴的土地資源，威脅糧食安全，其造成的破壞幾乎是不可逆的。現階段，青藏高原的土壤侵蝕研究主要集中在通用土壤流失方程的應用上[4,7-8]，對溝蝕的研究很少，溝蝕現狀和驅動因素還不明確[9-10]。因此，有必要識別溝蝕發生的高風險區和驅動因素，以便更有針對性地控制溝蝕[5]。

溝蝕是一個多因素作用下的復雜過程[5-6,11]，其影響因素包括地形、土壤、巖性、植被、降水、人類活動等[6]，探測溝蝕發生風險和量化影響因素是土壤侵蝕研究中的重要組成部分[5,12-13]。過去的溝蝕歸因分析中，地形因子在研究中的應用最為廣泛[14-15]，人類活動的影響多采用與道路距離來量化，缺乏綜合的指標；而且對于土壤、土地利用等類型變量，往往難以納入到傳統的線性回歸模型中。

研究中的常用方法主要為機器學習、二元統計模型和多準則決策方法[5]。這類研究大部分集中在伊朗，例如Amiri等[11]基于機器學習算法評估了伊朗Maharloo流域的溝蝕發生風險，發現土地利用是影響當地溝蝕的主要因素。此外，國內學者也采用機器學習方法評估了黃土高原小流域尺度[16]和東北黑土區中大流域[14]、小流域尺度[15]的溝蝕發生風險。但上述方法有著操作復雜的缺點，且部分方法基于主觀認識，可信度不足[17]。而地理探測器和最優尺度回歸在一定程度上可以消除上述的缺點，是研究多元非線性問題和類型變量的有力工具[18-19]。由于兩種方法操作簡單、結果可靠，被廣泛應用于自然科學、社會科學和二者的交叉領域[19-20]。但地理探測器和最優尺度回歸的現有應用多局限在對現象的歸因上[20-21]，缺乏在空間規律上的應用，而兩種方法在包括溝蝕風險在內的許多領域有著很大的應用潛力，因而有必要探索其在空間建模中的適應性，且目前青藏高原溝蝕發生風險的空間建模研究還是空白。為此，本研究基于地理探測器和最優尺度回歸等方法，量化影響溝蝕的各因子及其各分級/分類的重要性，進行空間建模，繪制拉薩河流域溝蝕發生風險圖，以期探明拉薩河流域溝蝕發生的空間異質性及驅動因子，為青藏高原國家生態安全屏障保護與建設和土壤侵蝕防控提供參考。

1 材料與方法

1.1 研究區概況

拉薩河流域（29°20′～31°15′ N，90°05′～93°20′ E）位于青藏高原南部（圖1），面積約32 000 km2，是西藏自治區的政治經濟中心。拉薩河是雅魯藏布江中段重要的支流之一，全長568 km。研究區地處高原半干旱溫帶季風氣候區，年平均降水量在396～726 mm之間，主要集中在6－9月，年平均氣溫在?8.4～9.1 ℃之間，年平均降雪量8.1 mm，融雪期一般在4－8月[22]，蒸發量大于2 000 mm[23]。地形以中山、高山、極高山為主。基巖包括白堊紀花崗巖、蛇綠巖和第三紀火山巖等[24]。地表覆被主要以草地和裸地為主，植被以高山草甸和灌叢草原為主[25]。土壤以砂土和砂壤土為主，礫石含量高[3]。土壤、土地資源十分稀缺，人類活動主要集中在拉薩河下游河谷階地和洪積扇上[1]。

注：青藏高原邊界來自青藏高原數據中心，由中國科學院地理科學與資源研究所張鐿鋰研究員繪制并提供。

1.2 數據來源

1.2.1 侵蝕溝調查與解譯

在2018－2019年，對拉薩河流域的溝蝕情況開展了野外考察，并通過大疆精靈Phantom 4 Pro無人機航拍了15個樣區的無人機影像。基于實地調查和無人機影像，建立了侵蝕溝在遙感影像上的識別特征（圖2）。之后，基于Google Earth軟件提供的亞米級高分辨率遙感影像解譯了2 171個樣點，其中溝蝕發生（有溝）樣點1 025個，未發生（無溝）樣點1 146個（圖1）。Google Earth提供的高分影像在流域不同區域多分布在2010—2020年間，在提取時選擇最新的影像。首先由ArcGIS生成1 000個隨機點，對其判別后發現，基本都是無溝樣點，之后又通過主動識別在流域手動添加了1 171個樣點，其中有溝樣點位置標在侵蝕溝橫斷面正中間。在分析中，有溝樣點賦值為1，無溝樣點賦值為0。

1.2.2 因子數據來源

溝蝕是一個受自然和人為因素綜合影響的復雜過程。因此本研究所用的因子數據主要包括地形、地表覆被、氣候、人類活動等數據（表1）。其中海拔數據來自先進陸地觀測衛星（Advanced Land Observing Satellite，ALOS）獲取的12.5 m分辨率數字高程模型（Digital Elevation Model，DEM），其余地形因子基于DEM數據計算得到。土壤數據來自世界土壤數據庫（Harmonized World Soil Database，HWSD），巖性數據來自土壤-土地數字化數據庫（Soil and Terrain Digital Database，SOTER），植被數據在Google Earth Engine平臺基于Landsat 8遙感衛星2020年度遙感影像通過最大值合成法獲得[11]。所有因子數據通過上述矢量點位置和ArcGIS多值提取至點工具獲取。

圖2 地面調查、無人機影像與衛星影像的侵蝕溝特征對比

表1 影響因子數據源

1.3 數據處理

1.3.1 變量離散化

多元線性回歸既不能納入類型變量，也無法表征復雜的非線性關系。而地理探測器和最優尺度回歸可用于類型變量分析，能夠更好地反映變量間的復雜關系，使得建模結果更準確更可靠[19-20]。因此，將連續數值型變量通過離散化的方式轉換為類型變量。除坡向為9類外，其余連續變量均采用自然間斷點法分為8個等級[21]（表2）。

1.3.2 多重共線性檢驗

建模時，因子間的多重共線性影響模型的可靠性，并干擾因子重要性排序。因此，本研究在數據離散化前后，基于SPSS 20軟件的多元線性回歸工具對所有因子間的多重共線性進行了檢驗，剔除了地表粗糙度指數、年平均氣溫等存在共線性的因子。通過容忍度和方差膨脹因子量化因子間共線性的程度。當容忍度小于0.2或方差膨脹因子大于5時，說明因子間存在共線性[11]。

表2 影響因子的分級/分類信息

1.4 溝蝕發生風險探測

1.4.1 地理探測器

地理探測器（Geodetector）的核心原理是通過空間異質性來探測因變量與自變量之間分布格局的一致性，包含因子探測器、交互探測器和風險探測器等[19]。因子探測器可以探測變量的空間分層異質性，以及探測某因子在多大程度上解釋變量的空間分異，用值度量，可表征各因子重要性；交互探測器原理是將兩個因子的分級/分類疊加，以疊加后的新分級/分類來計算值，以分析兩個因子的交互作用會增強還是削弱對變量的解釋力；風險探測器原理是計算各分級/分類內的屬性均值。本研究中屬性值為0（無溝）和1（有溝），因此風險探測器結果為有溝點占總樣點的比例，即溝蝕發生頻率（），可表征各分級/分類重要性。

在ArcGIS柵格計算器中采用綜合指數法[3]將地理探測器作為一種二元統計模型用于評估溝蝕發生風險，見式（1）。通過自然間斷點法將結果劃分為非常低、低、中、高、非常高5個風險等級。

式中1為地理探測器得到的各柵格位置上的溝蝕發生風險，q為因子探測器中因子的值，r為風險探測器得到的因子在各柵格位置上的溝蝕發生頻率，表示因子序號。

1.4.2 最優尺度回歸

最優尺度回歸，又稱分類回歸（Categorical Regression，CATREG），是通過為類別指定數值來量化分類數據，從而生成轉換后變量的最優線性回歸方程。其原理是基于模型效果最優化的原則，首先對原始變量進行變換，將各變量轉換為適當的、最佳的量化評分，然后使用量化評分代替原變量進行回歸分析。這樣不僅能將原本不能放入各種回歸模型的分類變量加入到線性回歸模型中，還能對連續變量離散后的類型量進行評分[21]。

回歸系數可以表征各因子重要性，量化評分可以表征各分級/分類重要性。基于SPSS 20構建最優尺度回歸模型后，使用ArcGIS柵格計算器工具應用該模型計算溝蝕風險并分級，見式（2）。

式中2為最優尺度回歸得到的各柵格位置上的溝蝕發生風險，a為最優尺度回歸中因子的系數，b為因子在各柵格位置上的量化評分，表示因子序號。

1.4.3 兩種方法的組合

在分別獨立使用兩種方法后，再將二者組合，分別用因子探測器的值與最優尺度回歸因子各分級/分類的量化評分（地探×最優）、最優尺度回歸的系數與風險探測器因子各分級/分類的發生頻率（最優×地探）來探測溝蝕發生的空間異質性（圖3）。

式中3和4分比為上述兩種方法得到的各柵格位置上的溝蝕發生風險。

1.4.4 模型可靠性評價

使用ArcGIS多值提取至點工具提取4種模型并獲得各矢量點溝蝕風險數值，并通過受試者工作特征（Receiver Operating Characteristic，ROC）曲線評價上述4種模型可靠性。ROC曲線由敏感度和1－特異度組成，曲線下的面積（Area Under Curve，AUC）被用來對比模型間可靠性的優劣[11]。AUC值越接近1，模型越可靠。AUC>0.7，意味著模型可接受；AUC>0.8，意味著模型良好；AUC>0.9，意味著模型優秀[26]。本研究的ROC曲線在Origin 2021b軟件中繪制。

注：ROC為受試者工作特征；q、r分別為地理探測器的因子探測器q值和風險探測器發生頻率；a、b分別為最優尺度回歸中因子系數和各分類量化評分；R1、R2、R3、R4為地理探測器、最優尺度回歸、地探q×最優a、最優b×地探r得到的各柵格位置上的溝蝕發生風險。下同。

2 結果與分析

2.1 影響因子共線性和重要性

各因子間多重共線性檢驗的容忍度均大于0.2，方差膨脹因子（Variance Inflation Factor，VIF）均小于5（表3），這表明不存在明顯的共線性，所有因子可同時應用于溝蝕發生風險的空間建模。

表3 影響因子間多重共線性檢驗

各因子的重要性主要包括最優尺度回歸的系數和因子探測器的值（圖4）。在最優尺度回歸中，因子系數前三位分別為海拔（0.442）、土壤類型（0.168）、歸一化植被指數（0.156），所有因子系數的值均小于0.001，達到了極顯著水平，通過了統計學顯著性檢驗。在地理探測器中，海拔、土壤類型、人類足跡占據了前三位，重要性（值）分別為0.263、0.251、0.174。除與河流距離<0.01外，其余因子值的顯著性檢驗值均小于0.001。所有因子值均達到了極顯著水平，說明各因子均有統計學意義。

圖4 各因子最優尺度回歸的系數和地理探測器q值

2.2 影響因子的分級或分類重要性

影響因子各分級/分類的重要性主要包括最優尺度回歸的量化評分和地理探測器的發生頻率（圖5）。可以看出，各因子和溝蝕發生之間的關系基本上都是非線性的，因此量化因子各分級/分類重要性是必要的。雖然在數值上不同，但兩種方法得到的因子分級/分類重要性變化趨勢相對一致。

圖5 不同影響因子各分級/分類的重要性

由圖5可知，海拔各等級的發生頻率和量化評分均隨海拔的升高而遞減。坡度、地形位置指數、地形起伏度、與河流距離均為第2、3等級對溝蝕風險的重要性最高；坡向中的東北和西最重要；地形濕度指數第7級最重要。人類足跡的第6、7等級，而不是第8等級溝蝕風險最高，這是由于人類足跡最高的地區為硬化地表，不易發生溝蝕。歸一化植被指數（NDVI）的第3等級，即0.2附近時，對溝蝕發生的重要性最高；由于植被最差的區域為冰川和水體，植被最差和最好的分類重要性均較低。各巖性中，砂巖、雜砂巖、長石砂巖/頁巖的重要性最高。各土壤類型中，簡育黑土、飽和疏松巖性土、黑色石灰薄層土和松軟薄層土的重要性較高。各土地利用中，灌木林地重要性最高，水體和不透水面最低。與居民點距離越近，溝蝕發生頻率越高。年平均降水量在第4、5等級，即高于500 mm左右后對溝蝕的重要性開始下降。

部分因子分級/分類的最優尺度回歸量化評分和風險探測器發生頻率略有差異，可能是由于地理探測器本質上是統計，各因子間是獨立的；而最優尺度回歸各因子的影響要相互平衡，會給一些因子分級/分類打低分以平衡其他因子在該范圍的影響。所以，地理探測器對單個因子各分級/分類的探測效果更好，最優尺度對整體風險探測更好。

2.3 溝蝕發生風險

從空間上來看，拉薩河谷下游河谷和當雄盆地被識別為高風險區，流域東北部和其他高海拔地區被識別為低風險區（圖6）。4種方法繪制的溝蝕發生風險的空間分異特征整體上相近。但最優尺度回歸、地探×最優兩種方法斑塊破碎化程度更高（圖6a、6c），另外兩種方法且斑塊更連續（圖6b、6d）。此外，地理探測器、最優×地探兩種方法高估了下游河道和河岸階地的風險。溝蝕發生的非常高風險等級占比在9.52%～13.97%之間，高風險等級占比在12.64%～20.99%之間。4種方法的風險等級占比略有不同。其中，最優尺度回歸和地探×最優的各等級占比基本一致，而地理探測器和最優×地探的各等級占比接近（表4）。

圖6 4種方法的溝蝕發生風險圖

表4 4種方法對應的溝蝕發生風險等級占比

由于海拔對溝蝕發生風險的空間分異最為重要（圖4），因此本研究單獨統計了最優尺度回歸得到的溝蝕發生風險沿海拔梯度的分異特征（圖7）。91%的非常高風險區集中在海拔的1～3等級，即3 523～4 599 m；79%的高風險區集中在海拔的2～5等級，即3 978～5 059 m。72%的非常低風險區和71%的低風險區集中在海拔的5～7等級，即4 842～5 572 m。總之，高風險主要集中在5 059 m以下，低風險主要集中在5 059 m以上。

圖7 基于最優尺度回歸的各海拔等級溝蝕風險等級柵格計數

2.4 模型可靠性

4種方法的ROC曲線和AUC值見圖8。其中，最優尺度回歸的表現最好，AUC值為0.899；地理探測器的AUC值為0.833。兩種組合方法中，地探×最優的表現相對更好，AUC值為0.866。4種方法都有良好的可靠性。此外，與最優尺度回歸結合可提高地理探測器的探測效果。

注：AUC表示曲線下面積。

3 討 論

3.1 侵蝕溝樣點選取

本研究的對象是遙感影像能夠識別的侵蝕溝，其規模上接近傳統意義的切溝和沖溝[27]。但當地侵蝕溝與黃土高原和東北黑土區的侵蝕溝形態、地貌、發育條件都有所區別，所以不能直接使用其他地區的分類系統，在提出成熟的分類系統之前，本文暫時都稱作侵蝕溝。本研究的侵蝕溝寬度多在5～50 m之間，而Google Earth提供的0.5 m分辨率高分影像能夠支持本研究侵蝕溝的判讀，其可靠性在以往研究中也得到了廣泛檢驗[28-29]。以往同類研究中，往往只關注該樣點是否為侵蝕溝[5,11,15]，而忽視了遙感影像與各因子數據來源不一致可能導致的輕微偏移。同時，10 m左右分辨率的影響因子數據在同類研究中表現最佳，分辨率更高則可能受數據偏移影響，分辨率更低則舍棄了一部分數據信息[30]。因此，將有溝點標繪在侵蝕溝橫斷面的正中間能夠代表該位置的真實情況。

3.2 因子離散化程度對結果的影響

分類的數量會在一定程度上影響建模結果[19,21]。隨著海拔分級數量增加，因子解釋力值逐漸增大，風險探測結果的細節也越多（圖9）。海拔分為8級時，最低海拔等級的溝蝕發生頻率最高（圖9a）。分為15級時，最低等級的發生頻率明顯低于第2～4等級（圖9b）。而分為30級時，最低等級海拔的溝蝕發生頻率僅0.35（圖9c）。但海拔30級時，后兩個等級沒有樣點（圖9c），建模結果不完整。因此為滿足所有因子各等級均有樣點，本研究將海拔等因子離散化為8級。今后的研究中可分別探索各因子值最高的分級/分類數量，以提高建模精度。對于樣點不足的可以考慮用鄰近等級的量化評分或發生頻率代替。

圖9 不同離散程度的各海拔等級的溝蝕發生頻率

3.3 方法可靠性及應用潛力

本研究首次將地理探測器和最優尺度回歸應用到空間建模上。最優尺度回歸的表現優于地理探測器和二者的組合，這表明它作為一種回歸模型要比統計模型有著更好的擬合效果。地理探測器高估了下游河道與兩岸階地溝蝕風險（圖10b、10d），可能是因為地理探測器是獨立統計各因子及其分級/分類的重要性[19]；而最優尺度回歸的各因子間能相互制衡，避免高估風險[18]（圖10c）。

溝蝕發生風險的研究中，常用的是機器學習、二元統計模型和多準則決策三類方法[5]。機器學習方法中，表現最佳的隨機森林等算法ROC曲線的AUC值大多能達到0.89以上[11]。二元統計模型中，AUC值大多能達到0.79以上[31]。多準則決策方法中，AUC值通常在0.85以上[17]。與以上3類方法相比，本研究4種方法的AUC值均達到了良好水平，最優尺度回歸的表現已經相對接近機器學習算法的擬合優度。機器學習是一個“黑箱”方法，而本文用到的方法可以定量評價各因子、及其各分級/分類的重要性，得到的模型比機器學習更直觀，更能揭示機理。此外，本研究的方法操作簡單，是兩種無代碼的建模方式，都適用于多元非線性建模。

注：圖a黑框為圖b-d范圍。

本研究中，海拔單因子能解釋因變量26.3%的變異性（<0.001）。而雙因子解釋力中，海拔與歸一化植被指數疊加的交互效應最高，能解釋因變量37.3%的變異性（<0.001）。最優尺度回歸的決定系數（2）為0.472（<0.001），表示最優尺度回歸模型所有自變量能解釋47.2%的因變量變異。溝蝕是一個多因素長時間驅動的復雜過程[5-6]。因此，對于二分類問題來說，通過了顯著性檢驗，且能解釋將近一半的因變量變異，說明本研究所選用的因子具有合理性[21]。

3.4 影響因素及防控對策

本研究中，海拔是影響溝蝕發生的最主要的因子，且海拔越低，溝蝕發生風險越高，這與伊朗等地[32]的結果相似。巨大的海拔落差是拉薩河流域的一個特點，因此海拔對溝蝕最重要的原因可能是其影響了徑流的分配，與高海拔地區相比，低海拔區徑流匯集范圍更大。另外，海拔還影響了人類活動范圍、土壤類型、土壤凍融和降水分布[1,33-34]。例如，謝健等[35]發現，在拉薩河流域西北的當雄縣，降水主要集中在距山谷地面高800 m以下的范圍。土壤類型在地理探測器和最優尺度回歸中的重要性都排在了第二位，說明土壤是影響溝蝕發生的重要因素。土壤是發生侵蝕的下墊面，是被侵蝕物質，其抗蝕性決定了侵蝕的難易程度[36]。研究發現簡育黑土、飽和疏松巖性土、黑色石灰薄層土和松軟薄層土的更易發生溝蝕。可能是由于這幾種土壤均較為疏松，且均分布于低海拔地區，發育程度較高[34]，導致更容易被侵蝕。此外，最優尺度回歸中NDVI的重要性排在海拔和土壤類型之后，這與黃土高原和東北黑土區的研究結果都相近。Zhao等[37]發現，在中國黃土高原，溝壑密度本身與地形和植被覆蓋顯著相關。而中國東北地區，溝蝕空間分異主要受地形起伏度、海拔等地形因子和 NDVI、年平均降水量等共同影響[12]。當NDVI在0.2左右時，對溝蝕發生的重要性最高，超過該閾值后，溝蝕發生的風險逐漸降低。這是由于植被覆蓋和根系固土能夠控制侵蝕。地理探測器中的人類足跡和與居民區距離的重要性較高。而本研究中溝蝕高風險區集中在拉薩河下游河谷和當雄盆地等人類活動強度較大的地區，說明人類活動對溝蝕發生有著重要影響。在兩個方法中，年平均降水量的重要性分別為第四和第六位。當年平均降水量高于500 mm后，對溝蝕發生的重要性開始下降，可能是由于500 mm以上的降水量促進植被生長，明顯抑制了溝蝕發生[9]。人類足跡和年平均降水量對溝蝕發生有重要影響，這意味著在人類活動增加和氣候變化的背景下青藏高原地區的溝蝕發生風險有著較大的變化潛力。

青藏高原可供人類利用的土壤土地資源十分稀缺[1]，而溝蝕是一個不可逆的資源破壞損失過程。因此，為實現青藏高原土壤土地資源的保護和合理利用，提出以下主要建議：1）拉薩河下游河谷是人類活動的主要區域，而且溝蝕發生風險最高，應重點關注，具體的防控措施包括減少表土擾動和促進植被恢復。2）與坡面上部相比，下游河谷兩岸下部坡面承受著更多的徑流，因此有著更高的溝蝕發生風險。可以在坡面水流路徑上恢復植被，并在已經形成的侵蝕溝底部布設引水、蓄水、排水的工程和植物措施[38]，以改變溝蝕匯水區的水文連通性，削減徑流量和徑流動能，降低侵蝕驅動力[9]。3）對于降水量較大、人類活動較少、且溝蝕發生風險較低的東北部地區，建議主要以減少和預防過度放牧帶來的表土擾動和植被破壞為主。

4 結 論

本文基于最優尺度回歸、地理探測器和二者結合的方法，選取地形、下墊面、氣象、人類活動等共15個影響溝蝕的因子，對拉薩河流域溝蝕發生風險進行了探測，得出以下結論：

1）在最優尺度回歸中，因子系數前三位分別為海拔（0.442）、土壤類型（0.168）、歸一化植被指數（0.156）；在地理探測器中，海拔（0.263）、土壤類型（0.251）、人類足跡（0.174）占據了因子重要性前三位；所有因子對溝蝕發生的影響均達到極顯著水平（<0.01）。

2）海拔3 523～5 059 m的區域、歸一化植被指數在0.2左右、土壤類型中的簡育黑土、飽和疏松巖性土、黑色石灰薄層土和松軟薄層土、以及人類活動的高等級區有著最高的溝蝕發生風險。

3）4種方法的擬合效果均較好，曲線下的面積（Area Under Curve，AUC）均大于0.8，且操作簡單，使用門檻低，適合廣泛應用于多元非線性的空間建模，其中最優尺度回歸的探測效果最佳（AUC=0.899）。

4）拉薩河流域有9.52%～13.97%的區域有著非常高的溝蝕風險，高風險主要集中在拉薩河下游河谷兩岸和當雄盆地等人類活動較多海拔較低地區，有必要對重要城鎮周邊采取改變水文連通性、減少表土擾動和適當植被恢復等措施。

[1] Chen T, Jiao J, Chen Y, et al. Distribution and land use characteristics of alluvial fans in the Lhasa River Basin, Tibet[J]. Journal of Geographical Sciences, 2021, 31(10): 1437-1452.

[2] Wang L, Liu H. Quantitative evaluation of Tibet's resource and environmental carrying capacity[J]. Journal of Mountain Science, 2019, 16(7): 1702-1714.

[3] Chen T, Jiao J, Zhang Z, et al. Soil quality evaluation of the alluvial fan in the Lhasa River Basin, Qinghai-Tibet Plateau[J]. Catena, 2022, 209: 105829.

[4] 陳同德，焦菊英，王顥霖，等. 青藏高原土壤侵蝕研究進展[J]. 土壤學報，2020，57(3)：547-564.

Chen Tongde, Jiao Juying, Wang Haolin, et al. Progress in research on soil erosion in Qinghai-Tibet Plateau[J]. Acta Pedologica Sinica, 2020, 57(3): 547-564. (in Chinese with English abstract)

[5] Vanmaercke M, Panagos P, Vanwalleghem T, et al. Measuring, modelling and managing gully erosion at large scales: A state of the art[J]. Earth-Science Reviews, 2021, 218: 103637.

[6] Poesen J, Nachtergaele J, Verstraeten G, et al. Gully erosion and environmental change: Importance and research needs[J]. Catena, 2003, 50(2): 91-133.

[7] He Q, Dai X A, Chen S. Assessing the effects of vegetation and precipitation on soil erosion in the Three-River Headwaters Region of the Qinghai-Tibet Plateau, China[J]. Journal of Arid Land, 2020, 12(5): 865-886.

[8] 王小丹，鐘祥浩，范建容. 西藏水土流失敏感性評價及其空間分異規律[J]. 地理學報，2004，59(2)：183-188.

Wang Xiaodan, Zhong Xianghao, Fan Jianrong. Assessment and spatial distribution of sensitivity of soil erosion in Tibet[J]. Acta Geographica Sinica, 2004, 59(2): 183-188. (in Chinese with English abstract)

[9] Li J, Zhao C, Chen T, et al. Gully erosion on alluvial fans can be mitigated by altering the hydrological connectivity between an alluvial fan and the contributing catchment: A study in the Lhasa River basin[J]. Land Degradation & Development, 2022, 33(8): 1170-1183.

[10] 趙春敬，焦菊英，稅軍鋒，等. 西藏中南部侵蝕溝形態無人機航測與傳統地面測量的對比分析[J]. 水土保持通報，2019，39(5)：120-127.

Zhao Chunjing, Jiao Juying, Shui Junfeng, et al. Comparative analysis on morphological characteristics of erosion gullies measured by an unmanned aerial vehicle and traditional ground survey in South Central Tibet[J]. Bulletin of Soil and Water Conservation, 2019, 39(5): 120-127. (in Chinese with English abstract)

[11] Amiri M, Pourghasemi H R, Ghanbarian G A, et al. Assessment of the importance of gully erosion effective factors using Boruta algorithm and its spatial modeling and mapping using three machine learning algorithms[J]. Geoderma, 2019, 340: 55-69.

[12] Zhou Y, Zhang B, Qin W, et al. Primary environmental factors controlling gully distribution at the local and regional scale: An example from Northeastern China[J]. International Soil and Water Conservation Research, 2020, 9(1): 58-68.

[13] Poesen J. Soil erosion in the Anthropocene: Research needs[J]. Earth Surface Processes and Landforms, 2018, 43(1): 64-84.

[14] 王文娟，鄧榮鑫，張樹文. 東北典型黑土區溝蝕發生風險評價研究[J]. 自然資源學報，2014，29(12)：2058-2067.

Wang Wenjuan, Deng Rongxin, Zhang Shuwen. Preliminary research on risk evaluation of gully erosion in typical black soil area of northeast China[J]. Journal of Natural Resources, 2014, 29(12): 2058-2067. (in Chinese with English abstract)

[15] Huang D, Su L, Zhou L, et al. Assessment of gully erosion susceptibility using different DEM-derived topographic factors in the black soil region of Northeast China[J/OL]. International Soil and Water Conservation Research, (2022-4-17)[2022-7-24]https://doi.org/10.1016/j.iswcr.2022.04.001.

[16] Jiang C, Fan W, Yu N, et al. Spatial modeling of gully head erosion on the Loess Plateau using a certainty factor and random forest model[J]. Science of The Total Environment, 2021, 783: 147040.

[17] Arabameri A, Pradhan B, Rezaei K, et al. Gully erosion susceptibility mapping using GIS-based multi-criteria decision analysis techniques[J]. Catena, 2019, 180: 282-297.

[18] van der Kooij A J, Meulman J J, Heiser W J. Local minima in categorical multiple regression[J]. Computational Statistics & Data Analysis, 2006, 50(2): 446-462.

[19] 王勁峰，徐成東. 地理探測器：原理與展望[J]. 地理學報，2017，72(1)：116-134.

Wang Jingfeng, Xu Chengdong. Geodetector: Principle and prospective[J]. Acta Geographica Sinica, 2017, 72(1): 116-134. (in Chinese with English abstract)

[20] Olivares B O, Calero J, Rey J C, et al. Correlation of banana productivity levels and soil morphological properties using regularized optimal scaling regression[J]. Catena, 2022, 208: 105718.

[21] 陳玉蘭，焦菊英，田紅衛，等. 黃土高原歸一化植被指數與自然環境因子的空間關聯性：基于地理探測器[J]. 生態學報，2022，42(9)：3569-3580.

Chen Yulan, Jiao Juying, Tian Hongwei, et al. Spatial correlation analysis between vegetation NDVI and natural environmental factors based on geographical detector on the Loess Plateau[J]. Acta Ecologica Sinica, 2022, 42(9): 3569-3580. (in Chinese with English abstract)

[22] 何曉紅，馬艷鮮. 拉薩降雪年際變化特征分析[J]. 西藏科技，2005(10)：54-56.

He Xiaohong, Ma Yanxian. Analysis of inter-annual variation of snowfall in Lhasa[J]. Tibet Science and Technology, 2005(10): 54-56. (in Chinese with English abstract)

[23] 張核真，卓瑪，向飛，等. 1981－2013年氣候因子變化對西藏拉薩河徑流的影響[J]. 冰川凍土，2015，37(5)：1304-1311.

Zhang Hezhen, Zhuoma, Xiang Fei, et al. Effect of climate factors on the runoff over Lhasa River basin during 1981-2013[J]. Journal of Glaciology and Geocryology, 2015, 37(5): 1304-1311. (in Chinese with English abstract)

[24] Sun J, Li S, Muhs D R, et al. Loess sedimentation in Tibet: provenance, processes, and link with Quaternary glaciations[J]. Quaternary Science Reviews, 2007, 26(17): 2265-2280.

[25] 林紅，焦菊英，陳同德，等. 西藏拉薩河流域中下游洪積扇植被的物種組成與多樣性特征[J]. 水土保持研究，2021，28(5)：67-75.

Lin Hong, Jiao Juying, Chen Tongde, et al. Species composition and diversity of vegetation of diluvial fan in the Lhasa River Basin of Tibet[J]. Research of Soil and Water Conservation, 2021, 28(5): 67-75. (in Chinese with English abstract)

[26] Alireza A, Wei C, Marco L, et al. Comparison of machine learning models for gully erosion susceptibility mapping[J]. Geoscience Frontiers, 2020, 11(5): 1609-1620.

[27] 趙宇輝，張建軍，于洋，等. 晉西黃土區蔡家川小流域切溝的空間分布及形態特征[J]. 農業工程學報，2022，38(4)：151-158.

Zhao Yuhui, Zhang Jianjun, Yu Yang, et al. Spatial distribution and characteristics of the gullies in Caijiachuan watershed in loess region of Western Shanxi Province, China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(4): 151-158. (in Chinese with English abstract)

[28] 李建軍，焦菊英，曹雪，等. 柴達木盆地沙丘移動的空間分異及對形態參數的響應[J]. 農業工程學報，2021，37(7)：309-314.

Li Jianjun, Jiao Juying, Cao Xue, et al. Spatial regionalization and response to morphological parameters of dune migration in the Qaidam Basin of China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(7): 309-314. (in Chinese with English abstract)

[29] 楊麗娟，王春梅，張春妹，等. 基于遙感影像研究極端暴雨條件下新成切溝發生發展規律[J]. 農業工程學報，2022，38(6)：96-104.

Yang Lijuan, Wang Chunmei, Zhang Chunmei, et al. Occurrence and development of newly formed gullies under extreme rainstorm conditions using remote sensing images[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(6): 96-104. (in Chinese with English abstract)

[30] Garosi Y, Sheklabadi M, Pourghasemi H R, et al. Comparison of differences in resolution and sources of controlling factors for gully erosion susceptibility mapping[J]. Geoderma, 2018, 330: 65-78.

[31] Zabihi M, Mirchooli F, Motevalli A, et al. Spatial modelling of gully erosion in Mazandaran Province, northern Iran[J]. Catena, 2018, 161: 1-13.

[32] Alireza A, Fatemeh R, Subodh C P, et al. Modelling of piping collapses and gully headcut landforms: Evaluating topographic variables from different types of DEM[J]. Geoscience Frontiers, 2021, 12(6): 135-152.

[33] 孫寶洋，吳志廣，李占斌，等. 凍融對土壤分離能力及侵蝕阻力的影響[J]. 農業工程學報，2020，36(11)：57-65.

Sun Baoyang, Wu Zhiguang, Li Zhanbin, et al. Effects of freeze-thaw on soil detachment capacity and erosion resistance[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(11): 57-65. (in Chinese with English abstract)

[34] Sklar L S, Riebe C S, Genetti J, et al. Downvalley fining of hillslope sediment in an alpine catchment: Implications for downstream fining of sediment flux in mountain rivers[J]. Earth Surface Processes and Landforms, 2020, 45(8): 1828-1845.

[35] 謝健，劉景時，杜明遠，等. 拉薩河流域高山水熱分布觀測結果分析[J]. 地理科學進展，2009，28(2)：223-230.

Xie Jian, Liu Jingshi, Du Mingyuan, et al. Analysis of the Observation Results of Temperature and Precipitation over an Alpine Mountain, the Lhasa River Basin[J]. Progress in Geography, 2009, 28(2): 223-230. (in Chinese with English abstract)

[36] Zhu T, Xu X, Zhu A, et al. Assessment of soil shear strength in gully systems: A case-study in the Wangjiagou (WJG) watershed on the Loess Plateau of China[J]. Land Degradation & Development, 2020, 31(17): 2648-2661.

[37] Zhao J, Vanmaercke M, Chen L, et al. Vegetation cover and topography rather than human disturbance control gully density and sediment production on the Chinese Loess Plateau[J]. Geomorphology, 2016, 274: 92-105.

[38] Molina A, Govers G, Putte A, et al. Assessing the reduction of the hydrological connectivity of gully systems through vegetation restoration: Field experiments and numerical modelling[J]. Hydrology and Earth System Sciences, 2009, 13(10): 1823-1836.

Detecting gully occurrence risks using multivariate nonlinear spatial modeling in the Lhasa River Basin of China

Li Jianjun1, Chen Yulan2, Jiao Juying1,2※, Chen Tongde1,3, Chen Yixian2,4, Zhao Wenting1, Zhao Chunjing1,5, Shang Tianshe1, Jian Jinshi1,2, Cao Xue2

(1.,,,712100,; 2.,,712100,; 3.,,810007,; 4.,,710127,; 5.,450003,)

The risk of soil erosion is ever increasing in the Tibetan Plateau, due to the ecological environment under climate warming and human activities. Among them, gully erosion has been the most severe type of soil erosion. Taking the Lhasa River Basin of China as a sample, this study aims to detect the gully occurrence risk using multivariate nonlinear spatial modeling. 2171 gully points were selected to determine the controlling factors of gully erosion using the field survey and remote sensing interpretation. Gully points were set as the value of 1, whereas, no gully point was the value of 0. The Categorical Regression (CATREG), Geodetector, and their combination were utilized to quantify and classify the importance of 15 factors for the risk map of gully occurrence. The 15 factors were the elevation, slope, aspect, topographic position index, topographic wetness index, topographic relief, Normalized Difference Vegetation Index (NDVI), lithology, soil type, permafrost thermal stability, land use, human footprint, distance to the residential area, distance to the river and mean annual precipitation, all of which passed the collinearity test before modeling. Continuous numerical variables were also converted to the type variables by means of discretization. Eight classes were divided by the Jenks Natural Breaks Classification, except for nine classes of aspect. All key factors were converted to the 12.5 m resolution raster datasets. The factor datasets of all sample points were extracted by multi-value extraction to point tool of ArcGIS 10.2 software for data analysis. The Geodetector was run in the Excel software. The factor and risk detectors were used as two parameters of the binary statistical model. The CATREG was performed on the SPSS 20 software. The fitting parameters were obtained, including the factor coefficient, and each class of factors. The improved model was executed in the Raster Calculator Tool in ArcGIS to obtain the map of gully occurrence risk. Finally, the goodness of fit of each was evaluated by the receiver operating characteristic (ROC) curve. The results showed that: 1) The top three factor coefficients in the CATREG were the elevation (0.442), soil type (0.168), and NDVI (0.156). By contrast, the elevation (0.263), soil type (0.251), and human footprint (0.174) were the top three among the Geodetector. 2) The regression performed the best (AUC=0.899), followed by the two combined methods (AUC=0.866/0.848). But, the Geodetector performed relatively poorly (AUC=0.833). All these methods were suitable for spatial modeling in this case. Only a relatively few influences were found in the classification number of factors on the modeling performance. Correspondingly, the classification number with the highestvalue of each factor in the Geodetector can be explored to improve the modeling accuracy in the future. In addition, the quantitative score or frequency of adjacent grade can be considered for the grade of sample deficiency. 3) About 9.52% to 13.97% of the Lhasa River Basin was at a very high risk of gully occurrence, mainly distributed in the lower Lhasa River valley and Damxung Basin. 91% of the very high risk areas were clustered in the elevation classes 1-3, whereas, 79% of the high risk areas were clustered in the elevation classes 2-5. Therefore, the hydrological connectivity can be altered to restore the vegetation for less topsoil disturbance in the low-elevation area around large towns. The high-elevation area with a low risk can be used to prevent overgrazing and vegetation destruction because of the high precipitation in the northeast of the basin. The finding can provide a strong reference to construct the ecological security barrier, as well as the soil and water conservation on the Tibetan Plateau.

gully erosion; remote sensing; spatial modeling; Categorical Regression; Geodetector; Tibetan Plateau

10.11975/j.issn.1002-6819.2022.17.008

S157.1

A

1002-6819(2022)-17-0073-10

李建軍，陳玉蘭，焦菊英，等. 基于多元非線性空間建模的拉薩河流域溝蝕發生風險探測[J]. 農業工程學報，2022，38(17)：73-82.doi：10.11975/j.issn.1002-6819.2022.17.008 http://www.tcsae.org

Li Jianjun, Chen Yulan, Jiao Juying, et al. Detecting gully occurrence risks using multivariate nonlinear spatial modeling in the Lhasa River Basin of China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(17): 73-82. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2022.17.008 http://www.tcsae.org

2022-07-24

2022-08-25

第二次青藏高原綜合科學考察研究（2019QZKK0603）；中國科學院戰略性先導科技專項（XDA20040202）

李建軍，博士生，研究方向為土壤侵蝕。Email：lijianjun@nwafu.edu.cn

焦菊英，博士，研究員，研究方向為流域侵蝕產沙、土壤侵蝕與植被關系及水土保持效益評價。Email：jyjiao@ms.iswc.ac.cn