融合合成孔徑雷達和光學影像的植被覆蓋區土壤濕度反演

謝 媛， 嚴頌華,2*， 馬秦生， 陳能成

(1.武漢大學電子信息學院， 武漢 430072； 2.武漢大學宇航科學與技術研究院， 武漢 430072；3.武漢大學測繪遙感信息工程國家重點實驗室， 武漢 430072)

在農業生產中，土壤水分對生態環境和水平衡起著至關重要的作用，也是研究土地和大氣之間物理過程的重要變量[1-4]。然而由于土壤濕度在時間、空間上具有時空易變性[5]，對土壤水分進行大尺度、高精度且連續的時空監測成為各大學者爭相研究的熱點。合成孔徑雷達(synthetic aperture radar,SAR)在估計土壤參數方面展現了良好的優勢[6-11]。對于裸土下土壤濕度反演模型，到目前為止，已形成了多種反演算法[12]，大致分為經驗方法、物理方法以及半經驗半物理方法幾個主要算法。最常用的是IEM(integral equation model)、CIEM(calibrate integral equation model)、AIEM(advanced integral equation model)、Oh和Dubois模型。這些模型的共同之處在于建立SAR數據和土壤參數之間的函數。同時，結合光學圖像等多源數據，有助于提高SAR反演土壤濕度的精度。

目前的研究方法中，包括主被動微波遙感相結合共同改進土壤濕度數據測量方法[13]，以及基于被動遙感與光學遙感相結合的手段進行降尺度獲取更高分辨率的土壤水分數據[14]取得了一定的進展，但是這些方法大多依賴于額外植被指數產品，缺乏光學遙感數據本身在植被方面突出表現的考慮。實際上，在地表幾何結構的影響中，植被的遮擋使得SAR接收到的垂直極化后向散射系數不僅包含了土壤表層后向散射系數，還包含了植被層散射系數。而對于半干旱，或有農田遮擋區域消除植被的影響能更準確地獲取土壤表面水分情況。

針對這一問題，依據主動微波遙感與光學遙感數據相結合，利用光學遙感圖像中反映的植被層信息，校正土壤水分反演中的植被層散射貢獻。其創新性在于既發揮了光學遙感高分辨率特性以及反映地表參數的優勢，又充分結合了主動微波遙感全方位、長時間序列的觀測特性，實現對土壤水分準確、大面積反演。與此同時采用經驗模型來作為土壤濕度反演算法更具有理論與試驗依據。

首先給出了總體技術路線，其次是數據處理部分，然后進行算法研究，其中包括用于求取植被覆蓋區地表后向散射系數的植被影響消除算法以及用于得到土壤濕度反演結果的裸土下土壤濕度反演算法，之后是基于算法研究給出的試驗結果及分析，最后給出結論。

1 總體技術路線

使用Sentinel-1A衛星作為主要數據源，Sentinel-2多光譜數據作為次要數據源，對武漢市豹澥區域的土壤濕度進行反演研究。如圖1所示為總體技術路線。

圖1 技術路線流程圖Fig.1 Flow chart of technical route in this paper

首先對SAR數據以及光學影像進行處理，然后進行植被影響的消除。即從光學影像中提取各個區域植被覆蓋度以及各項植被指數。利用考慮了植被覆蓋度的修正水云模型，對每個像素點的純粹土壤表面后向散射系數進行求解。相比于傳統的水云模型去除植被影響辦法[15-18]，考慮植被覆蓋度因素，使得對于每個像元，植被覆蓋部分的影響可以達到定量分析，將垂直極化后向散射系數由原始的植被層后向散射系數與經過衰減之后的裸土后向散射系數2個組成部分變成了3部分的組合——植被覆蓋部分植被層散射系數、植被覆蓋部分經過衰減的純粹地表散射系數、非植被覆蓋部分純粹地表散射系數。最后得出純粹地表后向散射系數。

將純粹地表后向散射系數代入由Dubois模型和Oh模型結合的裸土土壤濕度反演模型。因此對校正了植被層散射貢獻的研究區域進行土壤濕度反演，從而得到時間序列上與實測土壤濕度更為接近、更為準確的土壤濕度反演趨勢。

2 數據處理

2.1 SAR數據的預處理

Sentinel-1是歐洲航天局哥白尼計劃中的地球觀測衛星，由兩顆衛星組成，載有C波段合成孔徑雷達，本文中下載的數據模式為IW SLC模式，空間分辨率為5 m×20 m。圖2所示為2020年3月3號的Sentinel-1A數據干涉寬幅模式的斜距單視復數產品圖。

對Sentinel-1號數據預處理過程如下：在ESA的專用軟件SNAP里面對Sentinel-1號數據進行輻射定標、多視、濾波以及幾何校正，經過預處理后從SAR影像中可以得到每個像元的兩個參數信息：垂直極化后向散射系數σ和入射角λ。

2.2 光學數據處理

Sentinel-2號是高分辨率多光譜成像衛星，攜帶一枚多光譜成像儀，可覆蓋13個光譜波段。對光學數據進行處理目的是對原始影像進行大氣校正，再利用歸一化水體指數(normalized difference water index, NDWI)得到研究區域非水體影像。因此，在常規預處理基礎上，添加對光學數據進行水體裁剪的處理部分。光學數據處理分為兩步：常規預處理、裁剪水體。

2.2.1 常規預處理

光學數據的常規預處理主要對影像進行大氣校正，校正光譜分辨率。利用Sen2Cor處理庫對Sentinel-2 L1C級數據進行大氣校正，得到L2A級數據，采取10 m分辨率進行處理，導出ENVI文件。

2.2.2 裁剪水體

由于對于土壤濕度的反演工作針對的是陸地部分，若不進行水體的去除，則在反演過程中，水體會被認為是具有植被覆蓋部分，因此，水體的存在會對反演土壤濕度產生干擾。利用光學數據進行水體的去除手段是利用歸一化水體指數NDWI進行決策樹的分類。首先在ENVI(environment for visualizing images)軟件中，利用Layer Stacking模塊進行紅、綠、藍和近紅外四個波段的融合，然后計算NDWI，計算公式為

(1)

式(1)中：Rgreen和Rnir為綠波段和近紅外波段的光譜反射率；σNDWI為歸一化水體指數。

根據光譜圖像NDWI，將其進行分類，即在光學影像中可以認為σNDWI<0的部分為非水體部分，σNDWI≥0的部分為水體部分，再進行決策樹的分類，得到非水體部分shp圖，對原始圖像進行剪裁，即可得到去除水體之后的圖像，也就是光學影像處理完最后結果。

剔除水體之后的影像即為下一步處理的基礎影像，其中包含了四個波段的反射率Rred、Rgreen、Rblue、Rnir。

3 算法研究

實驗中所涉及的算法主要分為兩個步驟：第一步旨在通過修正的水云模型，從垂直極化后向散射系數中消除植被對于后向散射系數的影響；第二步則通過Dubois模型和Oh模型結合，利用介電常數與土壤濕度之間經驗關系，結合第一步求解的純粹土壤表面散射系數反演植被區域土壤濕度。

3.1 植被影響的消除

經過對Sentinel-1A SAR圖像進行預處理得到垂直極化后向散射系數σ與入射角θ，以及對光學影像數據處理得到新的光學影像之后，可以利用3.1.1節進行純粹土壤表面散射系數的求取，3.1.2節是實驗結果示例，3.1.3節則利用實驗結果分析了植被對于后向散射系數的影響與植被覆蓋度的相關性。

3.1.1 植被影響消除步驟

將預處理之后的光學影像作為輸入，求解純粹土壤表面散射系數步驟如下：

(1)利用波段反射率數據求解歸一化植被指數(normailized differential vegetation index,NDVI)。作為植被影響因素的重要參數，求解公式為

(2)

式(2)中：Rred和Rnir分別表示紅波段和近紅外波段的光譜反射率；σNDVI代表歸一化植被指數。

(2)利用NDVI指數，根據植被覆蓋度的定義求出每個像素點的植被覆蓋程度σVFC，計算方式為

(3)

式(3)中：σVFC為每個像素點的植被覆蓋程度參數；σNDVI,s為沒有植被覆蓋的像素點純粹裸土表面NDVI；σNDVI,veg為完全由植被覆蓋的像素點NDVI。本文定義σNDVI,s是圖像中累積概率為5%的NDVI，而σNDVI,veg是累積概率為95%的NDVI。

對于2020年3月20日的武漢小部分區域植被覆蓋度如圖3所示。

圖3 2020年3月20日武漢小部分區域(包含豹澥)植被覆蓋度分布圖Fig.3 Distribution map of vegetation coverage in a small area of Wuhan (including Baoxie area) on March 20, 2020

(3)利用NDVI以及VFC求解植被層散射系數以及去除植被之后純粹地表散射系數。

水云模型定義垂直極化后向散射系數由植被后向散射系數和植被下地表后向散射系數組成。可通過光學圖像求解每個像元植被覆蓋程度，考慮植被覆蓋度VFC，將一個像元看作是一部分植被覆蓋區域和另一部分裸土區域。因此本文中修正的水云模型定義垂直極化后向散射系數σ由以下三個部分組成：植被覆蓋地區的植被散射系數和經過衰減之后的植被覆蓋下純粹地表散射系數以及裸土區域的純粹地表散射系數，表現形式如圖4所示。

圖4 單個像素區域后向散射系數分布示意圖Fig.4 A pixel grid back scattering coefficient distribution model of modified water cloud model

由上述分析以及模型示意圖，可以得出修正水云模型形式為

σ=(1-σVFC)σs+σVFCσveg+σVFCτ2σs

(4)

式(4)中：雙層衰減因子(透過率)τ2和植被層后向散射系數求解公式[19]為

τ2=exp[-σNDVI/cosθ]

(5)

σveg=σNDVIcosθ(1-τ2)

(6)

式(4)中第1項為模型圖給出的非植被覆蓋地區純粹地表散射系數；第2項為模型圖中植被層散射系數；第3項為經過植被雙層衰減后純粹地表散射系數。σs為需要求解的純粹土壤表面散射系數。

(4)最后，可以根據SAR雷達影像總體后向散射系數σ和入射角度θ，以及由光學影像參數得到植被覆蓋區域純粹地表散射系數表達式為

(7)

3.1.2 植被消除影響試驗示例

為驗證3.1.1所述植被影響的消除算法，對試驗場地豹澥區域做了相關性研究。在豹澥區域布設的九個原位土壤濕度計實測站點情況如圖5所示。

圖5 原位土壤濕度計布設情況Fig.5 Layout of in situ soil hygrometer

其中每個站點相對于原始的SAR圖像中的垂直極化后向散射系數，去除植被散射系數之后都有一定程度的衰減，將光學圖像中每個站點于2020年3月20日的植被覆蓋度與在時間序列上植被影響的極化值平均衰減量做相關性驗證，所有站點中由于站點3、6、7原位土壤濕度計電池沒電故不予考慮。實驗數據主要由2019年11月28日至2020年4月20日這一時間段，組成時間序列。以2020年3月3日為例，利用植被消除算法對每個站點垂直極化后向散射系數σ進行植被影響的消除，得到純粹地表后向散射系數σs，對于每個站點，坐標以及去除植被前后垂直極化值、入射角情況如下表1所示。

3.1.3 植被消除影響實驗分析

下面分析植被覆蓋度與時間序列各站點垂直極化后向散射系數平均衰減量相關性。其中時間序列垂直極化后向散射系數平均衰減量用Δσ表示，植被覆蓋度選取的是2020年3月20日光學影像得到的每個站點的VFC，各站點Δσ計算方法如式(8)所示，計算結果如表2所示。

表1 2020年3月3日植被影響消除算法處理結果Table 1 Processing results of vegetation impact elimination algorithm on March 3, 2020

(8)

式中：求和范圍j為每個站點在時間序列2019年11月28日到2020年04月20日上的12個周期；σj為第j個周期站點i的后向散射系數，即植被消除算法輸入；σj,s為第j個周期站點i的植被消除算法輸出。

繪制植被覆蓋度與垂直極化值平均衰減量Δσvv相關曲線如圖6所示。

表2 原位土壤濕度計站點植被覆蓋度與ΔσTable 2 Vegetation coverage of in-situ soil hygrometer sites and Δσ

圖6 植被覆蓋度與垂直極化平均衰減量相關性曲線Fig.6 Correlation curve between vegetation coverage and average attenuation of vertical polarization

由表2、圖6可知，站點植被覆蓋度與時間序列植被對垂直極化值平均衰減量有著非常強的相關性，皮爾森相關系數達到了0.83，因此，對于植被覆蓋大的區域，更有必要去考慮植被的影響并消除。

3.2 裸土下土壤濕度反演

對于裸土下土壤濕度的反演，利用Dubois模型和Oh模型作為介電常數和均方根高度的求解方法。其中，Dubois模型給出了垂直極化值σvv和水平極化值σhh分別關于介電常數和土壤均方根高度的表達形式。Oh模型則給出了σvv和σhh之間的關系表達式。因此通過結合這兩個模型，首先對介電常數和均方根高度進行求解，其次，利用介電常數和土壤濕度之間的經驗方程，對裸土下土壤濕度值進行反演。

3.2.1 求解介電常數和均方根高度

(1)Dubois模型形式為

(9)

(10)

式中：k為波數；λ為波長；θ為入射角；ε為介電常數；s為均方根高度。

為以Dubois模型為基底建立方程組，需要求解水平極化值。對于水平極化值，引入Oh模型，依據Oh模型中垂直極化值與水平極化值比值來求解。

(2)Oh模型的表達形式為

(11)

(12)

式中：

(13)

(3)建立方程組。結合Dubois模型[式(9)、式(10)]以及Oh模型[式(11)]，也就是垂直極化值與水平極化值的比值，可以得到關于介電常數和均方根高度的二元方程，即

(14)

(15)

根據式(14)、式(15)組成的方程組，代入后向散射系數和入射角，即可求解出均方根高度s和介電常數ε。

3.2.2 根據介電常數反演土壤濕度

根據介電常數ε和土壤濕度Mv的經驗模型對土壤濕度進行反演，其中經驗模型為

(1.993+0.002a+0.015b)+(38.086-

0.176a-0.633b)Mv+(10.720+

(16)

式(16)中：a為黏土含量,%；b為沙粒含量,%，a、b在不同頻段下取值不同，本文中取頻率為5.55 GHz時含量，即a=0.515 1,b=0.134 3。

由此，即給出了利用垂直極化后向散射系數和入射角信息求解裸土區域土壤濕度的反演算法。

4 土壤濕度反演實驗及驗證分析

為了驗證土壤濕度反演趨勢的正確性， 選取2019年11月28日至2020年4月20日這一時間段進行分析，總共12景數據，其中，有2月20號1幅影像元數據丟失。土壤濕度反演步驟如第1節中技術路線所述，具體步驟如第2、3節詳細說明。

與實測土壤濕度比較結果分為兩個部分：時間序列土壤濕度反演趨勢結果分析和空間序列土壤濕度反演結果分析。

4.1 時間序列土壤濕度反演結果分析

根據以上處理結果，對所選時間序列土壤濕度進行反演。進行單個站點時間序列方面的分析。表3以站點1為例，列舉時間序列上反演的土壤濕度信息，反演模型輸入部分包括垂直極化后向散射系數、植被影響消除后純粹土壤表面散射系數、入射角，代入模型可以得到植被影響去除前后土壤濕度反演值。

其他站點情況如下，繪制其他站點Mv和Mvs分別與SM的相關曲線。

圖7所示為6個站點去除植被前后土壤濕度反演曲線與實測土壤濕度變化曲線的相關性，可以得出以下結論：

表3 站點1時間序列土壤濕度反演結果以及實測信息Table 3 Soil moisture inversion results and measured information on time series of site 1

(1)土壤濕度反演趨勢與原位濕度計的實測值變化趨勢都有很高的相關性。大部分站點皮爾遜相關系數均達到了0.6左右，即強相關。對于站點2 相關系數達到了0.8以上，呈現極強相關趨勢。

(2)相比于植被影響消除之前，去除植被影響后反演的土壤濕度值與實測值變化趨勢更相關，即證明實驗結果對于利用SAR影像反演土壤濕度的反演結果有了很好的提升。也證明利用Sentinel-1 SAR影像和Sentinel-2多光譜數據結合來反演土壤濕度一定時間內變化趨勢相比單利用SAR微波后向散射系數反演具有更高可信度。

4.2 空間序列植被去除前后土壤濕度反演圖像結果分析

圖8所示為2020年3月3日豹澥區域植被去除前后土壤濕度反演結果。

圖8 植被去除前后土壤濕度反演結果對比Fig.8 Comparison chat of soil moisture inversion results before and after vegetation removal

由于豹澥區域植被覆蓋面積較少，選取附近紅色區域發現相比于去除植被影響之前總體后向散射系數反演的土壤濕度，去除植被影響后，對于土壤濕度黃土部分和植被覆蓋代表的濕部分邊界更為明顯，整塊區域有明顯的幾何形狀，而未去除前的土壤濕度的區域劃分則很模糊。

同時，對于植被覆蓋度越大地區，理論上植被對于土壤濕度反演影響則越大，因此，繪制每個站點P2與P1差值曲線，與站點植被覆蓋度做相關性分析，結果如圖9所示。

圖9 植被覆蓋度與土壤濕度反演提升效果相關曲線Fig.9 Correlation curve of VFC and soil moisture inversion lifting effect

由圖9可以看出，對于植被覆蓋度越大的地區，去除植被前后土壤濕度與實測土壤濕度相關性差值越大，且兩者相關性達到了0.8以上，及證明試驗結果符合理論。

5 結論

利用Sentinel-1 SAR影像數據與Sentinel-2光學影像數據對土壤濕度反演過程中的植被影響進行了研究。通過對SAR圖像進行預處理，提取所研究區域的垂直極化后向散射系數以及入射角信息，同時利用Sentinel-2號的光學影像求解歸一化植被指數NDVI以及植被覆蓋度VFC，得到植被層散射系數和垂直極化衰減量，代入修正水云模型求解出植被覆蓋區域純粹地表后向散射系數，最后利用裸土下土壤濕度求解算法，對研究區域進行土壤濕度的反演。證實所述方法對比直接利用SAR影像垂直極化后向散射系數反演土壤濕度，提升了與實測值的相關系數，對于利用SAR影像反演土壤濕度的研究領域有著更好的參考價值和實際意義。

方法不足之處在于：沒有對植被的類型進行更進一步的研究，對于不同植被覆蓋產生的影響只是做了統一性處理，其次，在分辨率方面，最終結果采用的是5 m×20 m的空間分辨率，因此對于更加精細的區域進行土壤濕度的反演還存在一定的研究空間。在該論述方案中，對于空間序列方面土壤粗糙度的影響缺乏探討，克服這些問題，將是下一步研究的重點。