謝 媛， 嚴(yán)頌華,2*， 馬秦生， 陳能成

(1.武漢大學(xué)電子信息學(xué)院， 武漢 430072； 2.武漢大學(xué)宇航科學(xué)與技術(shù)研究院， 武漢 430072；3.武漢大學(xué)測繪遙感信息工程國家重點(diǎn)實(shí)驗(yàn)室， 武漢 430072)

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

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

針對這一問題，依據(jù)主動微波遙感與光學(xué)遙感數(shù)據(jù)相結(jié)合，利用光學(xué)遙感圖像中反映的植被層信息，校正土壤水分反演中的植被層散射貢獻(xiàn)。其創(chuàng)新性在于既發(fā)揮了光學(xué)遙感高分辨率特性以及反映地表參數(shù)的優(yōu)勢，又充分結(jié)合了主動微波遙感全方位、長時間序列的觀測特性，實(shí)現(xiàn)對土壤水分準(zhǔn)確、大面積反演。與此同時采用經(jīng)驗(yàn)?zāi)Ｐ蛠碜鳛橥寥罎穸确囱菟惴ǜ哂欣碚撆c試驗(yàn)依據(jù)。

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

1 總體技術(shù)路線

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

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

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

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

2 數(shù)據(jù)處理

2.1 SAR數(shù)據(jù)的預(yù)處理

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

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

2.2 光學(xué)數(shù)據(jù)處理

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

2.2.1 常規(guī)預(yù)處理

光學(xué)數(shù)據(jù)的常規(guī)預(yù)處理主要對影像進(jìn)行大氣校正，校正光譜分辨率。利用Sen2Cor處理庫對Sentinel-2 L1C級數(shù)據(jù)進(jìn)行大氣校正，得到L2A級數(shù)據(jù)，采取10 m分辨率進(jìn)行處理，導(dǎo)出ENVI文件。

2.2.2 裁剪水體

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

(1)

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

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

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

3 算法研究

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

3.1 植被影響的消除

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

3.1.1 植被影響消除步驟

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

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

(2)

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

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

(3)

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

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

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

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

水云模型定義垂直極化后向散射系數(shù)由植被后向散射系數(shù)和植被下地表后向散射系數(shù)組成?？赏ㄟ^光學(xué)圖像求解每個像元植被覆蓋程度，考慮植被覆蓋度VFC，將一個像元看作是一部分植被覆蓋區(qū)域和另一部分裸土區(qū)域。因此本文中修正的水云模型定義垂直極化后向散射系數(shù)σ由以下三個部分組成：植被覆蓋地區(qū)的植被散射系數(shù)和經(jīng)過衰減之后的植被覆蓋下純粹地表散射系數(shù)以及裸土區(qū)域的純粹地表散射系數(shù)，表現(xiàn)形式如圖4所示。

圖4 單個像素區(qū)域后向散射系數(shù)分布示意圖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和植被層后向散射系數(shù)求解公式[19]為

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

(5)

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

(6)

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

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

(7)

3.1.2 植被消除影響試驗(yàn)示例

為驗(yàn)證3.1.1所述植被影響的消除算法，對試驗(yàn)場地豹澥區(qū)域做了相關(guān)性研究。在豹澥區(qū)域布設(shè)的九個原位土壤濕度計實(shí)測站點(diǎn)情況如圖5所示。

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

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

3.1.3 植被消除影響實(shí)驗(yàn)分析

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

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

(8)

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

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

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

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

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

3.2 裸土下土壤濕度反演

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

3.2.1 求解介電常數(shù)和均方根高度

(1)Dubois模型形式為

(9)

(10)

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

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

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

(11)

(12)

式中：

(13)

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

(14)

(15)

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

3.2.2 根據(jù)介電常數(shù)反演土壤濕度

根據(jù)介電常數(shù)ε和土壤濕度Mv的經(jīng)驗(yàn)?zāi)Ｐ蛯ν寥罎穸冗M(jìn)行反演，其中經(jīng)驗(yàn)?zāi)Ｐ蜑?/p>

(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。

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

4 土壤濕度反演實(shí)驗(yàn)及驗(yàn)證分析

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

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

4.1 時間序列土壤濕度反演結(jié)果分析

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

其他站點(diǎn)情況如下，繪制其他站點(diǎn)Mv和Mvs分別與SM的相關(guān)曲線。

圖7所示為6個站點(diǎn)去除植被前后土壤濕度反演曲線與實(shí)測土壤濕度變化曲線的相關(guān)性，可以得出以下結(jié)論：

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

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

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

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

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

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

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

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

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

由圖9可以看出，對于植被覆蓋度越大的地區(qū)，去除植被前后土壤濕度與實(shí)測土壤濕度相關(guān)性差值越大，且兩者相關(guān)性達(dá)到了0.8以上，及證明試驗(yàn)結(jié)果符合理論。

5 結(jié)論

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

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