基于CNN網絡和遙感數據的土壤侵蝕評估方法

周海宏，雷 磊，姚普及，吳 健，王 良，劉 濤，景 龑，潘曉彤

(1.國網陜西省電力公司，陜西 西安 710048；2.國網陜西省電力公司電力科學研究院，陜西 西安 710100；3.陜西省東莊水利樞紐工程建設有限責任公司，陜西 咸陽 713208)

土壤侵蝕是一個復雜的4階段動態過程，包括土壤的剝離、破壞、遷移和隨后的沉積物沉積[1].侵蝕是一種自然發生的過程，它通過水和風的物理力量或耕作等與農業活動有關的力量，侵蝕農田表層土，從而影響所有地貌.土壤侵蝕既可以是一個漸進的過程，持續相對不被注意，也可能以驚人的速度發生，造成表層土壤的嚴重流失.土壤侵蝕會去除有機質和重要營養物質，阻礙植被生長，從而對整個生物多樣性產生負面影響[2-3].這種特殊現象被認為是對土壤肥力和生產力的最大威脅.它改變了土壤的物理、化學和生物特性，導致潛在農業生產力下降，特別是如今世界人口不斷增長的情況下會引起人們對糧食安全的擔憂.此外，土壤侵蝕的影響超出了肥沃土地的損失，因為它還可能導致溪流和河流的污染和沉積加劇，關閉水道，導致魚類和其他海洋物種減少.退化的土地也降低了蓄水能力，有時導致洪水加劇.

目前，國內外學者在土壤侵蝕測量方面采用了許多不同的研究方法.最初，用于估計土壤可蝕性風險的最常見模型是通用土壤流失方程(USLE)，該方程后來被修訂為修訂的通用土壤流失方程(RUSLE)[4].綜合或單獨使用許多創新技術，如衛星遙感[5]、地質統計學[6]、野外光譜學、機器學習[7]以及現場和實驗室土壤反射率聯合測量等等.在試圖調查土壤侵蝕狀況時，需要收集大量空間分布良好的樣本進行進一步分析，這通常是一個昂貴和耗時的過程.因此，將統計學和遙感相結合是一個有效的土壤侵蝕檢測方法.特別是地球觀測(EO)衛星技術的日趨成熟，為地質監測提供豐富的數據基礎.因此，有學者已經綜合利用衛星遙感數據、地理信息系統方法和RUSLE方法監測和估計土壤侵蝕率[8-9].然而，上述方法采用的USLE和RUSLE方程都是經驗模型，其參數包含不確定性.

1 土壤流失估計

RUSLE代表了氣候、土壤、地形和土地利用如何影響由雨滴沖擊和地表徑流引起的細溝和溝間土壤侵蝕.它已被廣泛應用于估算土壤侵蝕損失和評估土壤侵蝕風險.此外，RUSLE可以指出更好的開發和保護計劃，以控制不同土地覆蓋條件下的侵蝕，如農田、牧場和受干擾的林地.RUSLE方程是：

A=R×K×LS×C×P.

(1)

式中，A為年平均土壤流失量，單位為t ha-1ya-1年；R為降雨侵蝕因子，單位為MJ mm ha-1h-1ya-1；K為土壤可蝕性因子，單位為t h MJ-1mm-1；LS為地形因子(無量綱)，C為覆蓋管理因子(無量綱)，P為保護措施因子(無量綱).

1.1 降雨因子R

在USLE和RUSLE模型中，降雨侵蝕力R因子是降雨侵蝕力的一個指標，它是降雨引起侵蝕的潛在能力.R因子計算公式是：

(2)

其中，R為降雨侵蝕指標，Pi為月降雨量.

1.2 土壤因子K

土壤侵蝕K值計算公式是：

K=10-3×(160.8-2.31x1+0.38x2+2.26x3+1.31x4+14.67x5).

(3)

其中，K為土壤侵蝕指標，x1為細礫(1～3 mm)%，x2為細砂(0.05～0.25 mm)%，x3為粗粉礫(0.01～0.05 mm)%，x4為細粉礫(0.005～0.01 mm)%，x5為有機質(10 g/kg).

1.3 地形因子LS

本文采用數字高程模型(DEM)計算LS因子.LS因子在每個DEM像素中計算公式如下，

(4)

式中，Lij-in是網格單元(i,j)的斜率長度；Aij-in是坐標(i,j)的網格單元面積，單位為m2；D是網格單元尺寸，單位為m；m是通用土壤流失方程(USLE)L因子的長度指數；xij為(sinαij+cosαij).

此外，在RUSLE中，m根據細溝和層間侵蝕的比率β而變化，m和β計算公式如下:

(5)

式中，β隨坡度變化而變化.

β值由以下公式得出：

(6)

此外，坡度由以下公式計算：

(7)

其中θ是坡度，單位為度.

1.4 地表覆蓋管理因子C

在USLE/RUSLE模型中，覆蓋系數(C)是一個反映種植方式對土壤侵蝕率影響的指數，它以土地利用為基礎.本文以NDVI、土壤調節植被指數(SAVI)和改良土壤調節植被指數(MSAVI)(高分遙感影像不同波段組合確定[10]),3種植被指數作為神經網絡的輸入層，C因子為輸出層.

1.5 土壤保持措施因子P

遙感影像基礎上P值需要根據土壤的類型和不同流域管理治理等資料來確定.一般情況下無保持措施地區取1.0，梯田地區取0.15.

2 基于CNN的C因子評估過程

2.1 數據集

數據集為2016年至2021年5景高分辨率地理參考正射影像(分辨率為0.25～0.5 m).所有圖像的特性(例如空間分辨率、顏色分布和光照條件)略有不同，但卻是在7月底至9月初的植被生長季節進行記錄.

此外，根據國家測繪局1980年代測繪的1∶5萬地形圖，建立了空間分辨率為30 m的數字高程模型(digital elevation model，DEM)，利用研究區的比例尺對土地覆蓋等級進行了評價.

2.2 訓練過程

用于訓練CNN模型的數據包括NDVI、SAVI、MSAVI數據和訓練標簽.訓練標簽是將土地利用圖和遙感影像相疊加，提取各個土地利用類型內植被指數值，并根據不同土地利用類型的C值作為標簽(表1所示).

表1 不同土地利用類型的C值

CNN網絡結構如圖1所示.首先，網絡由卷積層1和ReLU激活構成，后面是最大層化層.對于每個最大層化層，生成的特征映射的大小將減半.在后續卷積層中，應用一系列具有ReLU激活的轉置卷積層，以恢復原始圖像大小.最后，一個1×1卷積層和一個像素級的softmax激活函數提供最終的C因子輸出.softmax函數將每個數據的激活重新調整為[0，1]間隔.

圖1 CNN網絡結構

文中利用均方誤差損失對CNN進行訓練，則預測值f(x)與樣本真實值c的損失函數計算如下，

(8)

其中n為樣本的個數.ci和fxi分別為第i個樣本的真實值及其對應的預測值.

3 案例分析

選取密云水庫流域數字高程圖及遙感圖像，并對其進行土壤侵蝕分析.圖2所示為該區域覆蓋范圍示意圖.接下來根據第1節相關知識對用降雨因子R、地形因子LS、地表覆蓋管理因子C、土壤因子K、土壤保持措施因子P進行計算.

圖2 密云水庫流域區域覆蓋范圍示意圖

首先，通過GIS對整個流域的R因子輸入圖進行樣條插值.圖3所示為案例區域R值的空間分布.可以看出，R值隨降水特征由西北向東南逐漸增大.該區域的最小和最大R值分別為0.351和 206.214 MJ mm ha-1h-1ya-1.其次，根據所研究區域地貌特征及公式(3)，K值范圍為0.117至 0.397 5 t h MJ-1mm-1.圖4所示為案例區域K值的空間分布.

圖3 R因子空間分布示意圖

圖4 K因子空間分布示意圖

接著利用第2節方法對C因子進行評估.根據前述可知，數據集為2016年至2021年5景高分辨率地理參考正射影像，將其轉化為NDVI、SAVI、MSAVI數據和訓練標簽，共包含 1 292 個訓練樣本.仿真環境使用TensorFlow完成模型搭建，訓練過程參數如表2所示.圖5所示為使用CNN網絡計算的C因子空間分布圖，值在0.004 1～0.108 9之間，且受人為干擾影響，在低洼處較高.此外，由于包含區域有限，假定該區域的P因子值為1.

表2 網絡訓練部分參數

圖5 C因子空間分布示意圖

圖6所示為30個實測采樣點與CNN評估C值之間的散點圖.可以看出C因子與CNN反演的相關系數r為0.929，標準差為0.048.雖然有個別采樣點與評估結果有偏差，但大部分采樣點能夠正確反映C因子結果.

圖6 利用CNN計算C因子

接著采用空間分辨率為30 m的地形圖和式(4)、(5)繪制了圖7所示的坡度小于或等于21%的坡度因子(LS)空間分布圖.可以看出，LS值與地表起伏有直接關系.該區域儲層的最高LS值計算為20.63.

圖7 LS因子空間分布示意圖 圖8 土壤侵蝕空間分布圖

在完成數據輸入程序并編制R、K、C、P和LS圖作為數據層后，在GIS環境下對其進行倍增，繪制出顯示研究區土壤侵蝕空間分布圖，以確定每個區域的侵蝕風險，如圖8所示.可以看出在研究區內，黃土丘陵區由于土壤可蝕性而具有較大的侵蝕風險.密云水庫流域一半以上屬于極低侵蝕風險等級，土壤流失量小于10 t ha-1ya-1.土壤流失量從流域東南向西北增加，西北部土壤流失量最大，超過 150 t ha-1ya-1.

4 結語

本文對基于遙感圖像的土壤侵蝕評估進行了研究，在此基礎上，創新性的提出了利用CNN模型求解RUSLE中C因子.通過仿真分析，結果表明本文所提方法能夠利用遙感圖像進行土壤侵蝕風險評估.