田湘云，張 超，陳 棋，史小蓉，王 妍

(西南林業(yè)大學，云南 昆明 650224)

石漠化由袁道先院士在1991年提出，并于1995年用rock desertification進行表達，在2012年《中國石漠化公報》中將石漠化定義為：在熱帶、亞熱帶濕潤、半濕潤氣候條件和巖溶極度發(fā)育的自然條件下，植被受人為干擾破壞，引發(fā)土壤侵蝕及土地退化加劇，基巖裸露面積增大的情況[1]。石漠化的形成機理大致分為6類：巖溶致密堅硬、成土速率太慢、土壤沖刷嚴重、巖石與土壤之間存在“軟硬界面”黏著力差，土壤易流失；土壤及水分通過孔隙流失；人為因素[2]。云南省喀斯特地區(qū)與石漠化地區(qū)面積分別為110 875.7 km2和35 722.7 km2，居全國石漠化面積第2位[3]。由于地形、降水、土壤和植被覆蓋度對石漠化的發(fā)育有著密切關系：從環(huán)境結構來看，地形因子與土壤含水量、保水能力及需水量的聯(lián)系密不可分；降水后不同坡度對地表徑流作用力不同，影響土壤聚積，土壤侵蝕危害隨地形起伏度的增加而增加，從而導致石漠化的產(chǎn)生；降水量越大對土地石漠化濺蝕性越強；植被覆蓋度與土壤抗蝕能力程度呈正相關，植被覆蓋度高可以保持水土、涵養(yǎng)水源、減少地表徑流對石漠化造成影響，降低石漠化發(fā)育程度。水土流失造成區(qū)域土地退化、生產(chǎn)力降低，導致耕作土壤層被侵蝕，易發(fā)生旱澇，降低水資源保護工程效益，對石漠化地域的生產(chǎn)生活產(chǎn)生嚴重影響。地貌類型在大尺度上直接影響地表徑流及雨水沖刷力，對自然環(huán)境和水土流失有著關鍵性作用[4]。

地形起伏度是區(qū)域的海拔高程和地表切割程度的整體表征，各學者從多角度提出不同的定義及表達式[5-7]，筆者認為地形起伏度是宏觀描述整體地形變化情況，定量分析地貌形態(tài)、類型及演化的重要指標更為確切[8]。地形起伏度在各學術領域得到更多關注，如景觀格局分異特征[9]、地表形變[10]、地質災害[11-13]、貧困耦合關系[14-15]、人口與經(jīng)濟分布特征[11-18]、土壤水分[19]、土壤侵蝕水土流失強度[20-21]等方面有著廣泛應用。提取地形起伏度時具有尺度相關性[22-23]，最佳分析窗口是較為常用的方法之一，按照數(shù)據(jù)源及研究范圍不同，采用相應分析方法確定最佳分析窗口。地形起伏度是影響潛在水土流失的地形因子的重要指標，譚瑋頤等[24]通過對喀斯特高原山區(qū)的地形起伏度研究，利用最佳分析窗口對水土流失的地理分布格局和地形起伏度與水土流失強度的相互關系進行分析。劉新華等[25]通過地形起伏度對水土流失實用性進行研究，利用最佳分析窗口將土壤侵蝕圖、地形起伏度圖及地貌圖單獨比較，結果表明在地形因子中地形起伏度適用于水土流失評價。

1 材料與方法

1.1 研究區(qū)概況

云南省文山州是我國石漠化防治的重點地區(qū)，屬滇東南喀斯特較為復雜的地區(qū)。馬關縣、西疇縣及麻栗坡縣地處103°52′-105°18′E、22°42′-23°37′N，土地總面積6 494.7 km2，海拔高差2 569 m，境內(nèi)山巒起伏地形復雜；年均氣溫16.8 ℃，年均降水量1 236 mm，夏少炙熱，冬少酷寒。區(qū)內(nèi)石漠化現(xiàn)象嚴重，是南亞熱帶中低山河谷盆地的典型石漠化區(qū)[3]。

1.2 數(shù)據(jù)收集與預處理

研究采用的數(shù)據(jù)主要包括馬關縣、西疇縣及麻栗坡縣基礎地理信息矢量圖層(比例尺1∶25萬)、數(shù)字高程模型(DEM，空間分辨率30 m)及Sentinel-2遙感影像數(shù)據(jù)(2020年成像，空間分辨率10 m)。

根據(jù)上述基礎數(shù)據(jù)，提取坡度、海拔數(shù)據(jù)，得到地形起伏度。針對Sentinel-2影像數(shù)據(jù)的預處理主要包括大氣校正及輻射定標，利用 ENVI軟件進行鑲嵌、邊界裁剪處理。利用 ArcGIS軟件中的 Reclassify對空間分辨率為30 m的DEM數(shù)據(jù)進行重采樣，生成空間分辨率為10 m的DEM數(shù)據(jù)，使影像內(nèi)容與地物特征一致，便于對土地利用及植被覆蓋度提取?；谄露?、土地利用類型與植被覆蓋度等指標，依據(jù)《土壤侵蝕分級分類標準(SL190-2007修訂標準)》[26]，以決策樹的基本原理對水土流失程度類別進行分析及界定，確定馬關縣、西疇縣及麻栗坡縣水土流失強度空間分布特征。

1.3 地形起伏度計算

地形起伏度(relief degree of land surface，RDLS，公式中表示RDLS)，借鑒相關學者的研究方法，將我國的基準山體高視作500 m，將地形起伏度作為獨立的數(shù)值從而符合地質學意義[27]，它是在綜合相對高程、土地總面積及平地面積等地形要素的基礎上反映其地面起伏情況，計算公式如下

RDLSn=(Hnmax-Hnmin)×(1-MA/A)/500

(1)

式中：RDLSn為n像元領域分析地形起伏度；Hnmax為n像元領域分析最高高程值；Hnmin為n像元領域分析最低高程值；MA為研究區(qū)平地面積，即參考相關文獻和3個縣均為喀斯特地貌的條件下，將坡度小于5°的區(qū)域定義為平地；A為研究區(qū)土地總面積。

計算不同統(tǒng)計單元地形起伏度，運用領域分析法提取馬關縣、西疇縣及麻栗坡縣的地形起伏度，像元分析窗口為n×n(n=2,3,4，…，70，共69個數(shù)列)的矩形窗口，在 ArcGIS中利用焦點統(tǒng)計的RANGE計算各像元分析窗口最高與最低高程間差值，統(tǒng)計各像元窗口平均地形起伏度。

1.4 最佳統(tǒng)計單元分析

在繁雜的地形環(huán)境中，根據(jù)特殊的石漠化地區(qū)和數(shù)據(jù)類型選取適當?shù)泥徲蚍治龃翱谀軠蚀_地說明地形起伏度變化[28-29]。為確定馬關縣、西疇縣及麻栗坡縣的DEM最佳分析窗口，利用ArcGIS軟件中的Neighborhood Analysis工具，選擇矩形窗口，按順序抽取n×n(n=2,3,4，…，70)尺寸的像元網(wǎng)格海拔極值，同時計算各像元分析窗口的單位面積平均起伏度。采用均值變點法計算最佳像元分析窗口，過程如下。

1)依次將各像元分析窗口下的平均起伏度值除以對應窗口面積，獲取單位面積地形起伏度序列Tn

Tn=tn/sn(n=2，3，4，…，70)

(2)

式中：Tn為像元分析窗口單位面積起伏度/m；tn為像元分析窗口平均起伏度值/m；sn為像元分析窗口面積/m2；n為像元分析窗口領域值。

2)對數(shù)列Tn取對數(shù)ln(Tn)，得到數(shù)列Xn(n=2，3，4，…，70)，計算數(shù)列的離差平方和S

(3)

(4)

(5)

4)計算S與Si的差值ΔS

ΔS=S-Si(i=2，3，4，…，70)

(6)

1.5 水土流失強度分級

為確保水土流失強度等級劃分的準確性，借鑒舒天竹等[20]、王嬌等[30]的方法，以國家水利行業(yè)標準為基礎[26,31]，結合石漠化地區(qū)具體情況，將坡度、高程、植被覆蓋度和土地利用作為評價水土流失強度的因子(表1)。利用決策樹分析法，確定研究區(qū)水土流失強度的空間分布情況。

表1 水土流失強度分級指標

1.5.1 高程和坡度 基于DEM數(shù)據(jù)，借助ArcGIS中3D Analyst模塊獲取其高程及坡度信息。按照水土流失強度分級標準，將研究區(qū)坡度分為Ⅰ級(0°～5°)，Ⅱ級(5°～8°)，Ⅲ級(8°～15°)，Ⅳ級(15°～25°)，Ⅴ級(25°～35°)和Ⅵ級(>35°)共6個等級。

1.5.2 土地利用種類分類 借助ArcGIS和ENVI軟件對研究區(qū)進行決策樹分類，按照石漠化地區(qū)土地資源利用方式，經(jīng)由影像判讀、解譯等過程，將研究區(qū)土地利用種類分為水域、林地、草地、耕地、建設用地和裸地6種，對土地利用種類的水土流失強度依據(jù)《中國土壤侵蝕遙感監(jiān)測》[32]的標準進行分級。

1.5.3 植被覆蓋度 利用ENVI對植被覆蓋度進行提取，計算公式為

(7)

式中：NDVI為歸一化植被指數(shù)；NIR為近紅外波段；RED為紅光波段。

植被覆蓋度FVC(公式中表示為FVC)計算公式為：

(8)

式中：NDVIveg為完全植被覆蓋地表的NDVI；NDVIsoil為裸土區(qū)域的NDVI。根據(jù)水土流失強度的界定標準并結合研究區(qū)的實際情況，將植被覆蓋度分為<30%、30%～45%、45%～60%、60%～75%及>75%，共5類等級。

2 結果與分析

2.1 最佳分析窗口

運用領域分析法計算像元分析窗口面積與地形起伏度間關系(表2)。將表2結果中的像元窗口大小與平均地形起伏度用對數(shù)方程進行擬合，擬合模型為y=6.001 3+5.081x-0.024 3x2，擬合度R2=0.999 4，擬合效果較好(圖1)。在像元分析窗口(2×2)～(16×16)間，平均地形起伏度曲線隨像元分析窗口增大而大幅度增加；在像元分析窗口為(17×17)～(55×55)，平均地形起伏度曲線隨像元分析窗口的增大由急變緩；在像元分析窗口56×56之后，隨著像元分析窗口的增大，曲線變化逐步平緩。

表2 窗口面積與平均地形起伏度值對照

為找出不同窗口序列的離差平方和S與窗口序列分段后的統(tǒng)計量Si間差值最大時對應的點。運用均值變點對式(4)進行統(tǒng)計，計算得到S=55.01。根據(jù)式(5)、(6)計算Si與ΔS，計算結果見表3，ΔS與像元分析窗口的變化情況見圖2。

由表3和圖2可知，ΔS數(shù)列在窗口序列22時最大，即對應像元分析窗口為24×24，ΔS=40.014。此時，對于10 m×10 m空間分辨率的DEM數(shù)據(jù)，所獲取的地形起伏度最佳像元分析窗口為24×24，領域統(tǒng)計面積為0.057 6 km2。

表3 均值變點統(tǒng)計

2.2 地形起伏度與水土流失強度空間分布特征

研究區(qū)使用空間分辨率為10 m的DEM數(shù)據(jù)，研究區(qū)總面積為6 494.71 km2，根據(jù)上述統(tǒng)計的最佳分析窗口24×24的像元大小，根據(jù)公式(1)計算研究區(qū)RDLS值在0～2.41。由圖3可知，RDLS最大值2.41位于研究區(qū)的東北部、西南部及南部地區(qū)，分別分布于馬關縣的南部與西部，西疇縣的東北部，麻栗坡縣的東北部、中部與西南部，此區(qū)域多為林地。RDLS最小值0位于研究區(qū)南部及中部地區(qū)，分別分布于馬關縣東部、西疇縣中部和麻栗坡縣西部，此區(qū)域多為耕地。從圖4可知，RDLS在0.15～0.25的地區(qū)所占比例最多，為48.95%，當RDLS為0.35時，面積累計頻率將近96.14%，當RDLS為0.45時面積累計頻率高達99.29%。

根據(jù)研究區(qū)水土流失強度等級及面積統(tǒng)計表(表4)與水土流失強度分布圖(圖5)可知，現(xiàn)階段研究區(qū)水土流失強度概況為：微度等級面積為1 297.39 km2，占研究區(qū)總面積的19.98%，地理分布高度集中于馬關縣中部，西疇縣南部及麻栗坡縣西部地區(qū)；輕度等級面積為999.96 km2，占研究區(qū)總面積的15.40%，地理分布主要集中在馬關縣中部及西部，西疇縣南部及西南部，麻栗坡縣西部及南部地區(qū)；中度等級面積為1 208.00 km2，占研究區(qū)總面積的18.60%，地理分布相對均勻；強烈等級面積為1 226.27 km2，占研究區(qū)總面積的18.88%，地理分布較均勻；極強烈等級面積1 743.53 km2，占研究區(qū)總面積的26.85%，地理分布較均勻；劇烈等級面積19.56 km2，占研究區(qū)總面積的0.30%，地理分布高度集中于研究區(qū)馬關縣西部、南部及東南部，西疇縣東北部，麻栗坡縣東北部、中部及西南部地區(qū)。

表4 研究區(qū)水土流失強度基本特征

2.3 地形起伏度對水土流失的影響

地形起伏度是評價水土流失強度極其重要的宏觀地形因子，由圖3和圖5可知，水土流失強度的分布規(guī)律與RDLS存在高度相關性，水土流失強度高的地區(qū)RDLS在0.25～2.41，整體分布于馬關縣的西部及東南部，西疇縣的東北部，麻栗坡的縣東北部、中部及西南部地區(qū)。水土流失強度較低地區(qū)，RDLS在0～0.25時，整體趨勢隨著RDLS的增長水土流失強度逐漸增強。當RDLS在0～0.05時，水土流失強度上升幅度較大；當RDLS在0.05～0.15時，水土流失強度上升緩慢；當RDLS在0.15～0.25時，水土流失強度上升幅度較大；當RDLS大于0.25時，水土流失強度愈加趨向于平緩。為探討地形起伏度對水土流失強度分布情況的影響，借助于ArcGIS軟件將水土流失強度數(shù)據(jù)轉換為10×10柵格大小的GRID類型，并統(tǒng)一研究區(qū)RDLS與水土流失強度GRID數(shù)據(jù)的行列數(shù)及坐標系，利用區(qū)域分析中面積制表工具對2類數(shù)據(jù)進行統(tǒng)計分析(表5)。

從表5和圖6中可知，微度水土流失等級的面積在0～0.15時，所占面積隨RDLS的增長呈急劇增加趨勢，在0.05～0.15處到達最大值，之后面積隨RDLS的增長呈下降趨勢；輕度、中度、強烈和極強烈水土流失等級面積在0～0.25時，面積隨RDLS的增長呈急劇增加趨勢，并在0.15～0.25處到達最大值，之后面積隨RDLS的增長呈下降趨勢；劇烈水土流失等級的面積在0～0.15時，面積隨RDLS的增長呈逐步上升趨勢，并在0.05～0.15處到達最大值，之后面積隨RDLS的增長呈下降趨勢，且在0.35～2.41處面積趨近于0。

表5 各地形起伏度下水土流失強度分布

在研究區(qū)RDLS中，0～0.05、0.05～0.15、0.15～0.25區(qū)域面積呈上升趨勢，占研究區(qū)總面積的0.41%、29.29%和49.41%，且在RDLS中0.15～0.25處出現(xiàn)區(qū)域面積峰值；RDLS中0.05～0.15、0.15～0.25、0.25～0.35等級總計區(qū)域面積超過研究區(qū)總面積的95%，說明研究區(qū)水土流失強度區(qū)域高度集中于RDLS的0.05～0.15、0.15～0.25、0.25～0.35等級；在RDLS中0～0.05、0.45～0.55、0.55～0.65、0.65～2.41等級區(qū)域面積僅占研究區(qū)總面積的0.41%、0.56%、0.11%和0.05%，而劇烈水土流失強度在RDLS中0.35～2.41范圍內(nèi)面積趨近于0，水土流失強度區(qū)域占比最小，是由于地形起伏度越高的區(qū)域，植被覆蓋度越低，巖石裸露率越高，降水沖刷力越強，已幾近無土壤能流失，因而水土流失強度極小。RDLS能夠反映石漠化地區(qū)地形起伏度的基本情況，可以準確反映水土流失的強度特征，是水土流失程度分布的重要影響要素。

3 結論與討論

3.1 結論

基于30 m空間分辨率的DEM數(shù)據(jù)通過重采樣得到10 m空間分辨率后提取地形起伏度，得到最佳像元分析窗口大小為24×24，統(tǒng)計面積為0.057 6 km2。研究區(qū)境內(nèi)的地形起伏度0～2.41，而起伏度在0.15～0.25等級的區(qū)域總面積占比最大(49.41%)，0.05～0.15等級占比次之(29.29%)，RDLS能較好反映喀斯特地區(qū)的地形地貌特征。

現(xiàn)階段水土流失強度主要為極強烈等級，占研究區(qū)總面積26.85%，其次為微度(19.98%)、強烈(18.88%)、中度(18.60%)、輕度(15.40%)和劇烈(0.30%)。微度侵蝕連片分布，并在馬關縣的中部及東部地區(qū)、西疇縣的北部地區(qū)和麻栗坡縣的西南部地區(qū)有集中分布的特征；輕度侵蝕主要呈點狀零散分布于馬關縣中部及東部地區(qū)、西疇縣及麻栗坡縣的西南部地區(qū)；中度侵蝕與強烈侵蝕連片分布；極強烈侵蝕連片分布，但馬關縣及西疇縣的中部地區(qū)分布較少，麻栗坡縣的西南部地區(qū)分布較少；劇烈侵蝕集中分布于馬關縣的西部及東南部地區(qū)，西疇縣東北部地區(qū)，麻栗坡縣的東北部、中部及西南部地區(qū)。根據(jù)各水土流失強度等級的面積變化，充分說明研究區(qū)水土流失為極敏感，總體應加強水土流失的綜合防治。

研究區(qū)地形起伏度與水土流失強度的相關性具有顯著特征，在RDLS中0.05～0.35等級范圍是水土流失敏感區(qū)域，在0.15～0.25等級范圍水土流失強度面積出現(xiàn)峰值，因而該等級范圍對水土流失強度最為顯著，之后隨RDLS的增加逐漸降低，此區(qū)域應作為文山州南部及東南部地區(qū)水土保持重點生態(tài)功能區(qū)水土流失綜合防治核心區(qū)域。

3.2 討論

石漠化作為喀斯特地區(qū)的生態(tài)景觀之一，其水土流失強度與石漠化的形成及發(fā)展呈正相關關系，而RDLS影響水土流失的強度及速率。關于RDLS和水土流失強度的自然背景及表現(xiàn)特征方面國內(nèi)外學者已開展大量相關研究，文山州作為典型喀斯特地區(qū)，是國家石漠化重點防治區(qū)域，特殊地質地貌使地下漏失隱蔽且復雜，若要準確分析土壤及水分通過孔隙的流失量，還需對裂隙中土壤的成分及時間進行判別，關于喀斯特石漠化水土地下滲漏有待進一步研究。本研究僅以地形起伏度因子對水土流失強度進行剖析，利用領域法分析法計算得到最佳像元分析窗口，分析RDLS對水土流失強度的影響。研究結果表明，隨著RDLS的變化，石漠化地區(qū)水土流失強度的空間分布特征明顯，這與舒天竹等[20]和譚瑋頤等[24]的結論一致，且擬合模型效果更好，擬合度R2=0.999 4。

RDLS是影響水土流失的主要因素之一，可作為喀斯特地區(qū)水土流失強度分析評價的重要指標，水土流失亦受降水、土壤、植被覆蓋度及人為活動的綜合影響，如對水土流失強度進行更深入研究，需全面考慮地貌因子、氣候因子及人為因素?，F(xiàn)階段研究區(qū)水土流失強度占比最大為極強烈等級，文山州應繼續(xù)加強森林植被的修復與保護，因地制宜考慮各地形區(qū)域的獨特性，采取具有針對性的治理措施能夠優(yōu)化空間資源配置，如建立植被防護帶，利用根系網(wǎng)固土壤；坡耕地實施坡改梯，能有效減沙降低土壤流失，陡坡耕地退耕還林還草[33]；采用節(jié)水農(nóng)業(yè)技術，不僅可以緩解水資源緊缺，還可保土保肥；加強法制宣傳，建立環(huán)境保護意識。實現(xiàn)多方位針對性遏制水土的流失和石漠化的蔓延，為喀斯特地區(qū)水土流失防治提供科學參考。