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

基于多元非線性空間建模的拉薩河流域溝蝕發(fā)生風(fēng)險(xiǎn)探測

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

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

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

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

0 引 言

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

1.2 數(shù)據(jù)來源

1.2.1 侵蝕溝調(diào)查與解譯

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

1.2.2 因子數(shù)據(jù)來源

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

圖2 地面調(diào)查、無人機(jī)影像與衛(wèi)星影像的侵蝕溝特征對比

表1 影響因子數(shù)據(jù)源

1.3 數(shù)據(jù)處理

1.3.1 變量離散化

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

1.3.2 多重共線性檢驗(yàn)

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

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

1.4 溝蝕發(fā)生風(fēng)險(xiǎn)探測

1.4.1 地理探測器

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

在ArcGIS柵格計(jì)算器中采用綜合指數(shù)法[3]將地理探測器作為一種二元統(tǒng)計(jì)模型用于評估溝蝕發(fā)生風(fēng)險(xiǎn)，見式（1）。通過自然間斷點(diǎn)法將結(jié)果劃分為非常低、低、中、高、非常高5個(gè)風(fēng)險(xiǎn)等級。

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

1.4.2 最優(yōu)尺度回歸

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

回歸系數(shù)可以表征各因子重要性，量化評分可以表征各分級/分類重要性?；赟PSS 20構(gòu)建最優(yōu)尺度回歸模型后，使用ArcGIS柵格計(jì)算器工具應(yīng)用該模型計(jì)算溝蝕風(fēng)險(xiǎn)并分級，見式（2）。

式中2為最優(yōu)尺度回歸得到的各柵格位置上的溝蝕發(fā)生風(fēng)險(xiǎn)，a為最優(yōu)尺度回歸中因子的系數(shù)，b為因子在各柵格位置上的量化評分，表示因子序號。

1.4.3 兩種方法的組合

在分別獨(dú)立使用兩種方法后，再將二者組合，分別用因子探測器的值與最優(yōu)尺度回歸因子各分級/分類的量化評分（地探×最優(yōu)）、最優(yōu)尺度回歸的系數(shù)與風(fēng)險(xiǎn)探測器因子各分級/分類的發(fā)生頻率（最優(yōu)×地探）來探測溝蝕發(fā)生的空間異質(zhì)性（圖3）。

式中3和4分比為上述兩種方法得到的各柵格位置上的溝蝕發(fā)生風(fēng)險(xiǎn)。

1.4.4 模型可靠性評價(jià)

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

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

2 結(jié)果與分析

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

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

表3 影響因子間多重共線性檢驗(yàn)

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

圖4 各因子最優(yōu)尺度回歸的系數(shù)和地理探測器q值

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

影響因子各分級/分類的重要性主要包括最優(yōu)尺度回歸的量化評分和地理探測器的發(fā)生頻率（圖5）?？梢钥闯?，各因子和溝蝕發(fā)生之間的關(guān)系基本上都是非線性的，因此量化因子各分級/分類重要性是必要的。雖然在數(shù)值上不同，但兩種方法得到的因子分級/分類重要性變化趨勢相對一致。

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

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

部分因子分級/分類的最優(yōu)尺度回歸量化評分和風(fēng)險(xiǎn)探測器發(fā)生頻率略有差異，可能是由于地理探測器本質(zhì)上是統(tǒng)計(jì)，各因子間是獨(dú)立的；而最優(yōu)尺度回歸各因子的影響要相互平衡，會給一些因子分級/分類打低分以平衡其他因子在該范圍的影響。所以，地理探測器對單個(gè)因子各分級/分類的探測效果更好，最優(yōu)尺度對整體風(fēng)險(xiǎn)探測更好。

2.3 溝蝕發(fā)生風(fēng)險(xiǎn)

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

圖6 4種方法的溝蝕發(fā)生風(fēng)險(xiǎn)圖

表4 4種方法對應(yīng)的溝蝕發(fā)生風(fēng)險(xiǎn)等級占比

由于海拔對溝蝕發(fā)生風(fēng)險(xiǎn)的空間分異最為重要（圖4），因此本研究單獨(dú)統(tǒng)計(jì)了最優(yōu)尺度回歸得到的溝蝕發(fā)生風(fēng)險(xiǎn)沿海拔梯度的分異特征（圖7）。91%的非常高風(fēng)險(xiǎn)區(qū)集中在海拔的1～3等級，即3 523～4 599 m；79%的高風(fēng)險(xiǎn)區(qū)集中在海拔的2～5等級，即3 978～5 059 m。72%的非常低風(fēng)險(xiǎn)區(qū)和71%的低風(fēng)險(xiǎn)區(qū)集中在海拔的5～7等級，即4 842～5 572 m?？傊?，高風(fēng)險(xiǎn)主要集中在5 059 m以下，低風(fēng)險(xiǎn)主要集中在5 059 m以上。

圖7 基于最優(yōu)尺度回歸的各海拔等級溝蝕風(fēng)險(xiǎn)等級柵格計(jì)數(shù)

2.4 模型可靠性

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

注：AUC表示曲線下面積。

3 討 論

3.1 侵蝕溝樣點(diǎn)選取

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

3.2 因子離散化程度對結(jié)果的影響

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

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

3.3 方法可靠性及應(yīng)用潛力

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

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

注：圖a黑框?yàn)閳Db-d范圍。

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

3.4 影響因素及防控對策

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

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

4 結(jié) 論

本文基于最優(yōu)尺度回歸、地理探測器和二者結(jié)合的方法，選取地形、下墊面、氣象、人類活動等共15個(gè)影響溝蝕的因子，對拉薩河流域溝蝕發(fā)生風(fēng)險(xiǎn)進(jìn)行了探測，得出以下結(jié)論：

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

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

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

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

[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ìn)展[J]. 土壤學(xué)報(bào)，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] 王小丹，鐘祥浩，范建容. 西藏水土流失敏感性評價(jià)及其空間分異規(guī)律[J]. 地理學(xué)報(bào)，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] 趙春敬，焦菊英，稅軍鋒，等. 西藏中南部侵蝕溝形態(tài)無人機(jī)航測與傳統(tǒng)地面測量的對比分析[J]. 水土保持通報(bào)，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] 王文娟，鄧榮鑫，張樹文. 東北典型黑土區(qū)溝蝕發(fā)生風(fēng)險(xiǎn)評價(jià)研究[J]. 自然資源學(xué)報(bào)，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]. 地理學(xué)報(bào)，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] 陳玉蘭，焦菊英，田紅衛(wèi)，等. 黃土高原歸一化植被指數(shù)與自然環(huán)境因子的空間關(guān)聯(lián)性：基于地理探測器[J]. 生態(tài)學(xué)報(bào)，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] 趙宇輝，張建軍，于洋，等. 晉西黃土區(qū)蔡家川小流域切溝的空間分布及形態(tài)特征[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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] 李建軍，焦菊英，曹雪，等. 柴達(dá)木盆地沙丘移動的空間分異及對形態(tài)參數(shù)的響應(yīng)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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] 楊麗娟，王春梅，張春妹，等. 基于遙感影像研究極端暴雨條件下新成切溝發(fā)生發(fā)展規(guī)律[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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]. 農(nóng)業(yè)工程學(xué)報(bào)，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] 謝健，劉景時(shí)，杜明遠(yuǎn)，等. 拉薩河流域高山水熱分布觀測結(jié)果分析[J]. 地理科學(xué)進(jìn)展，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

李建軍，陳玉蘭，焦菊英，等. 基于多元非線性空間建模的拉薩河流域溝蝕發(fā)生風(fēng)險(xiǎn)探測[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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

第二次青藏高原綜合科學(xué)考察研究（2019QZKK0603）；中國科學(xué)院戰(zhàn)略性先導(dǎo)科技專項(xiàng)（XDA20040202）

李建軍，博士生，研究方向?yàn)橥寥狼治g。Email：lijianjun@nwafu.edu.cn

焦菊英，博士，研究員，研究方向?yàn)榱饔蚯治g產(chǎn)沙、土壤侵蝕與植被關(guān)系及水土保持效益評價(jià)。Email：jyjiao@ms.iswc.ac.cn