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

(1.國(guó)網(wǎng)陜西省電力公司，陜西 西安 710048；2.國(guó)網(wǎng)陜西省電力公司電力科學(xué)研究院，陜西 西安 710100；3.陜西省東莊水利樞紐工程建設(shè)有限責(zé)任公司，陜西 咸陽(yáng) 713208)

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

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

1 土壤流失估計(jì)

RUSLE代表了氣候、土壤、地形和土地利用如何影響由雨滴沖擊和地表徑流引起的細(xì)溝和溝間土壤侵蝕.它已被廣泛應(yīng)用于估算土壤侵蝕損失和評(píng)估土壤侵蝕風(fēng)險(xiǎn).此外，RUSLE可以指出更好的開(kāi)發(fā)和保護(hù)計(jì)劃，以控制不同土地覆蓋條件下的侵蝕，如農(nóng)田、牧場(chǎng)和受干擾的林地.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為地形因子(無(wú)量綱)，C為覆蓋管理因子(無(wú)量綱)，P為保護(hù)措施因子(無(wú)量綱).

1.1 降雨因子R

在USLE和RUSLE模型中，降雨侵蝕力R因子是降雨侵蝕力的一個(gè)指標(biāo)，它是降雨引起侵蝕的潛在能力.R因子計(jì)算公式是：

(2)

其中，R為降雨侵蝕指標(biāo)，Pi為月降雨量.

1.2 土壤因子K

土壤侵蝕K值計(jì)算公式是：

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

(3)

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

1.3 地形因子LS

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

(4)

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

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

(5)

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

β值由以下公式得出：

(6)

此外，坡度由以下公式計(jì)算：

(7)

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

1.4 地表覆蓋管理因子C

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

1.5 土壤保持措施因子P

遙感影像基礎(chǔ)上P值需要根據(jù)土壤的類型和不同流域管理治理等資料來(lái)確定.一般情況下無(wú)保持措施地區(qū)取1.0，梯田地區(qū)取0.15.

2 基于CNN的C因子評(píng)估過(guò)程

2.1 數(shù)據(jù)集

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

此外，根據(jù)國(guó)家測(cè)繪局1980年代測(cè)繪的1∶5萬(wàn)地形圖，建立了空間分辨率為30 m的數(shù)字高程模型(digital elevation model，DEM)，利用研究區(qū)的比例尺對(duì)土地覆蓋等級(jí)進(jìn)行了評(píng)價(jià).

2.2 訓(xùn)練過(guò)程

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

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

CNN網(wǎng)絡(luò)結(jié)構(gòu)如圖1所示.首先，網(wǎng)絡(luò)由卷積層1和ReLU激活構(gòu)成，后面是最大層化層.對(duì)于每個(gè)最大層化層，生成的特征映射的大小將減半.在后續(xù)卷積層中，應(yīng)用一系列具有ReLU激活的轉(zhuǎn)置卷積層，以恢復(fù)原始圖像大小.最后，一個(gè)1×1卷積層和一個(gè)像素級(jí)的softmax激活函數(shù)提供最終的C因子輸出.softmax函數(shù)將每個(gè)數(shù)據(jù)的激活重新調(diào)整為[0，1]間隔.

圖1 CNN網(wǎng)絡(luò)結(jié)構(gòu)

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

(8)

其中n為樣本的個(gè)數(shù).ci和fxi分別為第i個(gè)樣本的真實(shí)值及其對(duì)應(yīng)的預(yù)測(cè)值.

3 案例分析

選取密云水庫(kù)流域數(shù)字高程圖及遙感圖像，并對(duì)其進(jìn)行土壤侵蝕分析.圖2所示為該區(qū)域覆蓋范圍示意圖.接下來(lái)根據(jù)第1節(jié)相關(guān)知識(shí)對(duì)用降雨因子R、地形因子LS、地表覆蓋管理因子C、土壤因子K、土壤保持措施因子P進(jìn)行計(jì)算.

圖2 密云水庫(kù)流域區(qū)域覆蓋范圍示意圖

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

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

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

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

表2 網(wǎng)絡(luò)訓(xùn)練部分參數(shù)

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

圖6所示為30個(gè)實(shí)測(cè)采樣點(diǎn)與CNN評(píng)估C值之間的散點(diǎn)圖.可以看出C因子與CNN反演的相關(guān)系數(shù)r為0.929，標(biāo)準(zhǔn)差為0.048.雖然有個(gè)別采樣點(diǎn)與評(píng)估結(jié)果有偏差，但大部分采樣點(diǎn)能夠正確反映C因子結(jié)果.

圖6 利用CNN計(jì)算C因子

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

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

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

4 結(jié)語(yǔ)

本文對(duì)基于遙感圖像的土壤侵蝕評(píng)估進(jìn)行了研究，在此基礎(chǔ)上，創(chuàng)新性的提出了利用CNN模型求解RUSLE中C因子.通過(guò)仿真分析，結(jié)果表明本文所提方法能夠利用遙感圖像進(jìn)行土壤侵蝕風(fēng)險(xiǎn)評(píng)估.