魯中山區CSLE模型中植被和地形因子的優化方法

邢先雙, 董明明, 郭 靜, 孟 琳, 張 玉, 趙登良, 龐海威, 邊 振

(1.山東省水文中心, 山東 濟南 250002; 2.濰坊市水文中心, 山東 濰坊 266071; 3.濟南大學 水利與環境學院, 山東 濟南 250022)

土壤侵蝕對生態系統的影響日益引起人們的重視。長期以來，對土壤侵蝕產生泥沙及泥沙運移的定量研究一直是相關學者的研究重點[1]。遙感技術作為開展大范圍、長時間序列土壤侵蝕研究的重要技術手段，在中國土壤侵蝕研究中發揮了巨大作用[2]。

目前，中國水土流失遙感監測方法主要有綜合評判法和中國土壤流失方程(CSLE， Chinese Soil Loss Equation)模型法[3]。CSLE充分考慮了中國的土壤流失特征，對模型中的各個變量都進行了適合中國自然地理狀況的系統性研究[4]，因此，相比于美國通用土壤侵蝕模型(USLE, universal soil loss equation)[5]和修正后的美國通用土壤侵蝕模型(RUSLE, revised universal soil loss equation)[6]，CSLE更適用于中國的水土流失狀況的評估。近些年，隨著遙感技術、地理信息系統技術、全球定位系統技術的日漸成熟，ArcGIS 10.7平臺被廣泛應用于土壤侵蝕因子的監測，將土壤侵蝕模型與ArcGIS 10.7平臺結合的研究也越來越多。

利用多源遙感數據的優勢，有助于解決長時間序列內高分辨率植被覆蓋信息獲取、提高土壤侵蝕因子監測的準確性[7]，解決了傳統水土流失動態監測主要采用中低分辨率遙感數據，其精度難以滿足當前水土流失監測與治理需要的問題。本文在前人研究的基礎上，將無人機低空遙感數據融入模型計算中，對CSLE模型中的植被覆蓋與生物措施因子和地形地貌因子進行參數優化及精度檢驗，并使用優化后的模型與國家監測的結果進行擬合分析，驗證優化后模型的可行性，優化模型方法對同類或相似地區開展水土流失動態監測和水土保持工作有重要的意義。

1 研究區概況

郝峪小流域位于山東省淄博市博山區池上鎮西池村(118°05′29″—118°08′06″E，36°19′52″—36°21′23″N)，屬小清河支流淄河上游源頭區，總面積14.4 km2，在全國水土保持區劃中屬于北方土石山區(北方山地丘陵區)—泰沂及膠東山地丘陵區—魯中南低山丘陵土壤保持區。研究區海拔高度范圍在350～657 m之間，相對高差307 m，地勢東南高西北低，整體向西北傾斜。

研究區屬暖溫帶半濕潤大陸性季風氣候，多年平均降水量715.1 mm，最大降水量1 332.7 mm(1964年)，最小降水量331.7 mm(1989年)，年平均蒸發量為1 201.4 mm，多年平均氣溫12.7 ℃，最高氣溫42 ℃，最低氣溫-21.8 ℃；年平均日照時數2 607 h；主導風向以南、西風為主，平均風速3.2 m/s。研究區上游頂部為片麻巖及黑云二長混合花崗巖，左半流域為各種陰影混合巖，右半流域為黑云斜長石片麻巖，流域內無斷層、裂隙漏水現象。由此母質發育形成的土壤以褐土和棕壤為主。研究區植被覆蓋較好，喬木有刺槐、楊樹、側柏、花椒、板栗、桃樹、蘋果樹等，灌木黃荊、酸棗等，草本有黃草、白背草、刺猬皮等，植被覆蓋率達90%。

2 研究方法

2.1 數據的獲取與處理

(1) 遙感數據(2019年)。主要包括用于解譯土地利用的GF-1衛星影像(分辨率2 m)；用于計算歸一化植被指數(NDVI)的Landsat系列影像(分辨率30 m)；用于對土地利用解譯核查、植被覆蓋度精確計算的研究區內無人機航拍數據(優于0.1 m)。遙感影像坐標及投影標準：CGCS 2000國家大地坐標系，1985國家高程基準，投影方式為正軸等面積割圓錐投影(Albers投影)。

(2) 土地利用解譯與核查。土地利用解譯采用室內人工目視解譯與野外調查相結合的方式進行。基于所獲取的高分辨率(優于2 m)遙感數據，結合ArcGIS 10.7平臺對研究區土地利用類型進行目視解譯，結合航拍數據進行修改、核查，確保解譯精度。解譯研究區2019年土地利用類型總圖斑數為119個，野外調查圖斑數為15個，占總圖斑數的12.61%。土地利用圖斑的解譯屬性及邊界吻合正確率為95.97%。

(3) 水土流失數據。淄博郝峪小流域控制站2019年產流產沙資料與國家水土流失動態監測成果中研究區2019年土壤侵蝕模數與各因子柵格數據。

(4) DEM數據。在地理空間數據云網站(http:∥www.gscloud.cn/)獲取研究區30 m分辨率DEM數據，進行檢查、裁剪、重采樣等處理，作為坡度、坡長數據獲取以及水土保持措施解譯和土壤侵蝕綜合分析的基礎。在流域內選取兩處調查樣區(圖1)，利用搭載多光譜相機的無人機獲取了樣區的DSM數據。

圖1 郝峪小流域無人機拍攝位置示意圖

2.2 土壤侵蝕模型CSLE計算方法

在《2019年度水土流失動態監測技術指南》的指導下，利用已獲取的研究區基礎數據，基于ArcGIS 10.7平臺，引入中國土壤侵蝕模型CSLE對研究區土壤侵蝕進行測算，計算模型為：

A=R·K·L·S·B·E·T

(1)

式中：A表示土壤侵蝕模數〔t/(hm2·a)〕;R表示降雨侵蝕力因子〔(MJ·mm)/(hm2·h·a〕;K表示土壤可蝕性因子〔(t·h)/(MJ·mm)〕;LS表示坡長坡度因子,無量綱;B表示植被覆蓋于生物措施因子,無量綱;E表示水土保持工程措施因子,無量綱;T表示耕作措施因子,無量綱。

2.2.1 地形地貌因子LS地形地貌因子包括坡長因子和坡度因子，表征的是地形對坡面土壤侵蝕速率的控制程度。采用符素華等[8]提出的公式計算坡度、坡長因子，計算公式分別為：

坡長因子計算公式：

(2)

式中：λ為坡長(m);m為坡長指數,隨坡度而變。

(3)

坡度因子計算公式：

(4)

式中：S為坡度因子(無量綱);θ為坡度(°)。

2.2.2 植被覆蓋與生物措施因子B植被蓋率指植物群落總體或各個體的地上部分的垂直投影面積與樣方面積之比的百分數。一切形式的植被覆蓋，均可不同程度地抑制水土流失[9]。基于Landsat系列中分辨率多光譜影像(包括藍、綠、紅和近紅外4個波段)及無人機低空遙感數據，根據土地利用類型，結合24個半月降雨侵蝕力因子比例計算B因子。園地、林地和草地采用公式計算，其余土地利用類型直接查表1進行賦值。

(1) 園地、林地B因子計算公式：

(5)

(6)

式中：WRi為第i個半月降雨侵蝕力占全年侵蝕力比例,取值范圍為0～1; SLRi為第i個半月園地、林地和草地的土壤流失比例,無量綱,取值范圍為0～1。

(2) 灌木林地SLRi計算公式：

(7)

(8)

(9)

式中：NDVI為歸一化植被指數; NIR為近紅外波段的反射率;R為可見光紅波波段的反射率; FVC為基于NDVI計算的植被覆蓋度,取值范圍為0～1。

(3) 果園、有林地和其他林地SLRi計算公式：

SLRi=0.444 68×e-3.200 96×GD-

0.040 99×eFVC(1-GD)+0.025

(10)

2.2.3 其他因子(R，K，E，T因子) 為研究L，S，B因子優化后對模型計算結果的影響，本研究的R，K，E，T因子與國家水土流失動態監測成果中的柵格數據保持一致，利用研究區邊界裁剪后代入公式進行計算。

2.3 模型評價方法

為使CSLE模型更適合魯中山區實際情況，本研究從提高原始數據空間分辨率角度，利用無人機正射影像與傾斜攝影技術，對模型中的地形地貌因子和植被覆蓋與生物措施因子進行了參數優化，并進行優化精度驗證。采用的因子優化及驗證方法為： ①擬合優度統計分析。基于SPSS分析平臺，使用相關性分析工具，確定兩組數據擬合曲線，并得出擬合優度的統計量：確定系數R2。R2衡量的是回歸方程整體的擬合程度，表達的是變量之間的總體關系其取值范圍為：0～1，越接近1說明模型的效率越高[10]。本研究使用該方法對CSLE模型中的植被覆蓋度進行了擬合分析，優化了模型中的植被覆蓋與生物措施因子。 ②直方圖匹配規則。在本研究所選取數據的時間內，流域的地形地貌可以認為是不變的，因此可以使用直方圖匹配方法對CSLE模型中參數地形地貌因子(LS)優化結果進行評價。直方圖指的是LS因子分布中每一因子值與其出現頻數的統計關系，用橫坐標表示研究區的LS因子值，縱坐標表示頻數。

對于像地形地貌因子這樣客觀不變的參數，其因子值的統計分布差異是由于研究人員所選取的研究方法導致。相關研究表明，基于不同分辨率的地形數據所提取的LS因子分布差異較大，地形數據的空間分辨率越高，LS因子值越精確[11]。在此研究基礎上，通過無人機航攝獲取的高程數據提取的高分辨率地形地貌因子LShigh可對基于30 m分辨率DEM數據提取的LSmedium進行修正。

表1 非園地、林地、草地的B因子賦值

基于ArcGIS 10.7平臺，提取LShigh和LSmedium的因子值累計頻率，根據直方圖匹配規則，即從累計頻率曲線中指定任意因子值LSh，查看其對應的LShigh的累計頻率，然后查看LSmedium中對應頻率值的因子值LSm，當無法得到對應的頻率值時，使用差分法計算得到LSm。用LSh替換修正LSm，即可得到修正后的，LSmedium修正，并得到參數修正的擬合曲線。通過對比分析尺度轉換前后的LS因子值頻率曲線的直方圖相似度(HS)，間接驗證尺度轉換模型的精度。用兩組數據的累計頻率直方圖面積差相對指數表示直方圖相似度，其計算公式為：

(11)

式中：Shigh表示LShigh的累計頻率分布曲線面積；Smedium表示LSmedium的累計頻率分布曲線的面積,可以通過微積分算得。HS值域為[0,1],HS值越接近1,兩組數據的相似度越高。

3 結果與分析

3.1 土壤侵蝕模型因子的計算與優化結果

3.1.1 地形地貌因子 為提高山東省水土保持動態監測項目中的模型計算精度，本研究以具有完整水文資料的魯中山區淄博市郝峪小流域為研究區，利用ASTER GDEM數據(空間分辨率30 m)，結合ArcGIS 10.7平臺中的水文分析工具提取并計算得出基于中分辨率數據的研究區地形地貌因子并重采樣為10 m分辨率，記為LSmedium。利用無人機對研究區中兩個有代表性地形的區域進行傾斜攝影，提取DSM數據(空間分辨率優于0.1 m)(圖2)。通過點云數據濾波處理，獲取地面點后編輯得到的DSM消除了樹木和建筑對地形的影響，可作為高分辨率的DEM數據。相同方法計算LS因子并重采樣為10 m分辨率，得出高分辨率數據的地形地貌因子LShigh。

在無人機傾斜攝影區域選取100個像元，利用ArcGIS 10.7工具將整個研究區及選取的100個像元的坡度進行分級，結合耕地和園林草的坡度分級將坡度分為8組(表2)。從表2可以看出，研究區和所選的100個像元坡度都以 8°～15°和15°～25°為主，坡度在0°～2°和>35°時最少。同樣，將整個研究區和選取的100個像元的LS因子分成8組(表3)。從表3可以看出，研究區和所選的100個像元中LS因子值都以10～20，20～30為主，LS值在70～80時最少。由此可以說明所選取的100個像元值具有代表性，可進行尺度轉換。

表2 研究區及所選數據坡度的像元數目分布

圖2 地理空間數據云提取的DEM以及無人機遙感影像提取的DSM數據

表3 研究區及所選數據LS因子的像元數目分布

將提取LSmedium和LShigh值進行相關曲線繪制，建立中分辨率與高分辨率地形地貌因子的尺度轉換公式，使用該公式對整個項目區LS因子進行修正，得到修正的研究區LSmedium修正。

(12)

利用直方圖匹配規則，設置修正步長選取為0.1，繪制LSmedium和LSmedium修正曲線(圖3)。

尺度轉換前后，通過LS的分布頻率和累計分布頻率圖(圖3)分析表明，LSmedium經過尺度轉換后生成的LSmedium修正頻率曲線與LShigh比較接近，這表明兩組數據的統計分布非常相似，研究區在LS因子值小于30的區域，LS因子值呈現錯位分布，研究區大于30的區域，二者分布幾乎重合，這說明LS因子尺度轉換模型有一定效果。

圖3 研究區LS因子尺度轉換前后頻率及累計頻率分布

經計算，尺度轉換后的LS因子值頻率曲線相似度為0.91，這表明兩組因子值的相似度較高，尺度轉換后的LS因子可以用來計算研究區的土壤侵蝕。選取無人機傾斜攝影范圍內不同高程梯度共30個像元的LSmedium修正值與無人機傾斜攝影提取的DSM數據進行擬合分析，線性相關曲線R2為0.697 1。進一步運用統計分布和直方圖相似度進行精度驗證。統計分布是指統計分析LS因子轉換前后的平均值、標準差、最大值和最小值以及頻率和累計頻率曲線圖。LShigh和LSmedium以及LSmedium修正的頻率曲線圖(圖3)表明，LSmedium頻率曲線峰值較小，且峰值附近頻率高，數據較集中于12～25之間。LShigh和LSmedium修正頻率曲線特征相似，峰值較LSmedium有所增大，且數值分布更加均勻。統計結果(表4)表明，修正后的LSmedium修正平均值、標準差相較LSmedium與LShigh更加接近，LSmedium修正整體較中分辨率的LSmedium增大4.6%，標準差較LSmedium增大8.5%。以上頻率特征的差異主要是由于中分辨率DEM數據由于空間分辨率較低，局部高程差異被消除，高程變化趨于平滑，高程變幅也有所減小，導致坡度、坡長因子值和空間變幅整體小于實際數據。而無人機DSM數據由于空間分辨率有較大提高，能較精確的反映細微的地形差異，因此基于無人機的DSM計算的LS因子LShigh和修正后的LSmedium修正能夠更準確地反映實際地形地貌信息，其精度較LSmedium有所提高。

表4 尺度轉換前后LS因子值統計結果

3.1.2 植被覆蓋因子與生物措施因子 基于不同分辨率的遙感數據提取的植被覆蓋度會有較大的差異，相較于Landsat系列衛星影像，搭載多光譜相機的無人機平臺可以提取更為精確的植被覆蓋數據。但是后者的耗時較大，成本較高，所以本研究于植物充分生長且云量相對較少的9月，在研究區選取兩個樣地進行無人機影像采集，經處理計算后生成分辨率為0.1 m的無人機正射影像。

基于Landsat影像提取的24期植被覆蓋度柵格數據中有兩期與無人機正射影像的時相接近。使用柵格計算器對兩期植被覆蓋度柵格數據求均值并重采樣后得到分辨率為10 m的中分辨率覆蓋度柵格數據FVCmedium。使用Creat Fishnet工具創建與FVCmedium像元邊界完全重合的解譯網格(大小為10 m×10 m)。使用Zonal工具分區統計各10 m×10 m網格內植被覆蓋度，得到10 m分辨率高精度植被覆蓋度柵格數據 FVCUVA。基于采樣空間分布均勻和各土地利用類型柵格數占柵格總數比例一致的原則，選取1 000個采樣點，使用ArcGIS 10.7中Sample工具提取FVCUVA及對應的FVCmedium像元值，共采集1 000個數值對。各土地利用類型采樣點選取情況詳見表5。

表5 各土地利用類型柵格數據采樣比例

對兩組柵格數據進行擬合分析，不同函數關系下的相關系數(表6)表明，FVCUVA和 FVCmedium呈現較好的相關性，各類函數的R2范圍在0.565 2～0.686 8之間，選擇R2最大的冪函數作為植被覆蓋度的尺度轉換模型，其表達式為：

(13)

式中：FVChigh為FVCUVA與FVCmedium最優擬合方程計算后的優化后植被覆蓋度柵格數據,FVCmedium為同時相中分辨率影像計算的植被覆蓋度。

FVChigh和FVCmedium的冪函數關系分布散點圖如圖4所示。

表6 FVChigh和FVCmedium擬合分析

圖4 FVChigh和FVCmedium柵格數值分布

利用上文公式對基于中分辨率的FVC進行修正，得到研究區高分辨率植被覆蓋度FVChigh。在無人機拍攝區域有林地和果園根據覆蓋度梯度抽取100個像元進行修正后植被覆蓋度與正射影像覆蓋度FVCUVA擬合分析(圖5)，R2為0.718 7，表明優化后的植被覆蓋度FVChigh更接近實際植被覆蓋度。

3.2 基于多源遙感的水土流失監測結果與分析

利用ArcGIS 10.7平臺將郝峪小流域的土壤侵蝕數據從山東省土壤侵蝕數據提取出，再通過優化后的模型計算出郝峪小流域的土壤侵蝕強度，結合《2019年度水土流失動態監測技術指南》中對土壤侵蝕量的分級標準，利用重分類工具進行侵蝕強度分級(圖6)。

圖5 郝峪小流域FVChigh和FVCUVA擬合分析

圖6 郝峪小流域原侵蝕強度和優化后侵蝕強度等級分布

利用ArcGIS 10.7平臺的分區統計工具對國家監測和優化后的的土壤侵蝕強度進行面積統計，結果詳見表7。從表7可以看出，優化LS因子和B因子后的模型計算所得水土流失面積較國家監測的水土流失面積增加0.85 km2，且全部為輕度侵蝕。為更好地驗證優化后模型的可行性，隨機從國家監測的土壤侵蝕數據和優化模型后的土壤侵蝕數據中提取100組數據進行擬合分析(圖7)，R2為0.700 3，說明優化后的土壤侵蝕數據與國家監測的土壤侵蝕數據具有良好的相關性，同時又存在一定差異。由線性擬合方程可知，優化模型因子后，侵蝕模數變化范圍有所減小，優化模型因子后的土壤侵蝕模數計算結果，在較小值區域大于國家監測數據，即國家監測成果中侵蝕模數小于947.53 t/(hm2·a)的區域經優化計算后，數值有所增大；反之在侵蝕模數值較大區域，優化計算結果小于國家監測成果。

表7 國家監測數據與優化模型后土壤侵蝕強度面積分析

魯中山區中度及以上侵蝕主要發生在灌木林地和旱地。原始數據空間分辨率的提高，使植被覆蓋度較中分辨率數據計算結果有所增大，減小B因子值；而地形因子計算中能夠更好地識別梯田區域從而減小LS因子值。輕度和微度侵蝕主要發生在有林地，為植被覆蓋度較高區域，優化后的B因子有所增大，導致侵蝕模數和輕度侵蝕強度面積的增加。同時由于微度和輕度侵蝕面積(合計14.41 km2)占研究區總面積(14.44 km2)的99.79%，因此微度和輕度侵蝕間相互轉移的概率最大，且侵蝕模數在某一閾值范圍內，對應的侵蝕強度不會發生變化，由實際計算結果統計可知，中度及以上侵蝕區域內的土壤侵蝕模數有所變化，但侵蝕強度未發生變化。

經過現場調查發現，研究區近年來旅游業發展和生態建設效果良好，不同坡面、坡度與土地利用類型區水土流失狀況差異不大，與優化后變化特征一致，但仍缺乏定量資料，驗證優化后模型的計算精度。

圖7 優化模型后的數據與國家監測數據的相關分析

4 結 論

(1) 本研究基于多元遙感數據，通過無人機高分辨率影像與中分辨率數據的擬合分析與精度驗證，從空間分辨率的提高方面，提高了DEM和FVC解譯精度。利用直方圖匹配規則，實現了LS因子的尺度轉換，驗證結果表明，LS因子頻率曲線相似度為0.91，可以進行研究區的土壤侵蝕模型計算。對多源遙感數據植被覆蓋度提取的尺度轉換進行了探討，得出R2為0.686 8的冪函數轉換模型。

(2) 對CSLE模型中因子進行優化的基礎上計算了研究區土壤侵蝕強度，并與該區域國家監測土壤侵蝕數據進行了對比分析。結果表明，優化模型因子后的研究區侵蝕模數空間異質性有所減小，輕度侵蝕面積有所增加，中度及以上侵蝕面積未變化，但侵蝕模數總體減小，這與LS因子和B因子優化后數值變化特征一致，說明模型因子的優化對監測結果具有顯著影響，但對侵蝕模數計算精度的提高是否具有顯著作用，仍需從土壤侵蝕機理上結合進一步試驗研究加以驗證。

(3) 不同空間分辨率的遙感數據源對土壤侵蝕的各因子計算均會造成一定的影響，優化后的LS因子和B因子更加接近無人機DSM數據處理的LS值和正射影像下的植被覆蓋情況，即以上因子精度均有所提高，繼而有助于提高模型計算結果精度，使得基于優化后因子計算的土壤侵蝕模數更符合當地實際情況。但本研究目前僅對地形地貌因子以及植被覆蓋與生物措施因子進行了尺度轉換和精度提升，對其他因子如何利用高分辨率數據源進行優化，并提高土壤侵蝕模型計算精度仍需進一步探討。