基于多時相SAR的IEM和Oh模型地表土壤含水率反演與融合

黃 對,劉 九 夫，張 建 云,王 文,崔 巍

(1.南京水利科學研究院水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210029；2.河海大學水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210098)

0 引言

合成孔徑雷達(SAR)空間分辨率高、不受云雨影響，其后向散射系數與地表粗糙度、土壤含水率、植被等因素密切相關，故其L、C、X波段影像適用于大范圍地表土壤含水率反演[1，2]。基于SAR影像的土壤含水率反演實際是尋找后向散射系數與土壤含水率關系的過程。相關研究采用經驗模型[3-6]、理論模型[7-15]和半經驗模型[16-19]進行裸露地表土壤含水率反演，一些學者對理論模型進行了改進，如Baghdadi等[20]對地表粗糙度進行定標以獲得更真實的積分方程模型(Intergral Equation Model,IEM)模擬數據，后續有學者開展了基于粗糙度定標的IEM[21]、AIEM[22]土壤含水率反演。對于植被覆蓋地表，則需利用水云模型或MIMICS模型去除植被的影響，然后采用裸露地表反演模型進行反演[23-27]。

受模型條件、反演方法與雷達成像特征限制，目前研究主要有基于單期影像[6，16]、兩期影像[17，23]或多期影像[18，26]單一模型以及基于單期影像多模型[20，22，27，28]的反演方法，反演過程常需依賴實測土壤含水率數據。不同模型的反演結果必然存在差異，而只利用單期或兩期影像的反演也具有一定局限性，由于影像特性(極化方式、入射角)不同，即使基于同一模型其反演細節也會有差異，由此導致不同時相反演結果的精度差異與時空連續性問題，有必要探索基于多時相SAR影像的土壤含水率多模型/多反演方法及反演精度評估，并對時空連續性較差的數據進行融合研究，以獲取多時相高精度土壤含水率數據。鑒于此，本文首先利用水云模型分析植被因素的影響，分別基于地表粗糙度定標的IEM、Oh模型，在不利用實測土壤含水率的條件下建立多時相SAR影像的土壤含水率反演方法，利用實測數據分析反演精度，進一步對兩種模型的反演結果進行加權平均與算術平均融合，探討融合方法的適用性，為基于多時相SAR的多模型/方法土壤含水率反演以及多時相高精度土壤含水率數據獲取提供參考。

1 研究區域及數據

1.1 研究區域

研究區為美國亞利桑那州Walnut Gulch流域(31°40′～31°47′N，109°50′～110°10′W)，屬半干旱區，面積約150 km2。流域年均氣溫17.7 ℃，年均降雨量330 mm，海拔1 200～1 900 m，地表覆蓋為稀疏植被，主要為草地與灌叢[29]。該流域是美國農業研究中心(USDA-ARS，http://www.tucson.ars.ag.gov/dap/)和NASA陸地水文研究的長期研究區域，可提供豐富的觀測資料，為基于主動微波遙感反演土壤含水率提供了先驗知識與驗證資料。

1.2 研究數據

(1)土壤含水率數據。研究區共有19個土壤含水率觀測站點(圖1)，每20 min記錄0～5 cm的土壤體積含水率等參數，已經過鹽分等誤差校正，采樣時間為2004年7月1日-9月30日[30]。

圖1 研究區與地表覆蓋Fig.1 Study area and land cover

(2)SAR影像。采用5期ENVISAT-ASAR(C波段)交叉極化模式成像數據Level 1級產品，為雙極化、地距數字影像[31](表1)，下載網址為http://earth.esa.int。利用NEST軟件對影像進行預處理：以5×5像元窗口，采用GAMMA MAP濾波進行影像去噪；進行地理編碼與輻射定標，將影像DN值轉化為雷達后向散射系數值σ0(m2/m2)；結合DEM(SRTM 1sec grid數據，通過NEST軟件獲取)進行亮度訂正[32]，糾正地形導致的輻射與幾何畸變，并將σ0(m2/m2)轉化成分貝值σ0(dB)，同時將雷達影像重采樣為30 m，并利用研究區邊界矢量裁切影像。

表1 研究區合成孔徑雷達影像信息Table 1 Information of synthetic aperture radar images in the study area

(3)植被含水量數據。植被含水量(mveg)數據由歸一化差分紅外指數(NDII)和實測植被含水量數據擬合得到(式(1))[33]，NDII由TM影像第4、5波段光譜反射率計算得到(式(2))。研究區包含0805、0818、0824共 3期植被含水量數據，0714期、0720期數據缺失，以0729期的數據代替。

mveg=0.557692308×NDII+0.13923

(1)

NDII=(B4-B5)/(B4+B5)

(2)

2 研究方法

2.1 水云模型

(3)

(4)

γ(θ)2=exp(-2Bmvegsecθ)

(5)

(6)

2.2 IEM模型

基于IEM模型的土壤含水率反演可分為經驗[35]和半經驗方法[18]，二者均基于模型模擬數據的散射特征分析構建反演公式，區別在于構建時是否利用實測土壤含水率[23]。

(7)

(8)

式中：p為VH極化方式；A、B、C為對應入射角與極化方式下的系數。

(9)

2.3 Oh模型

(10)

(11)

本文基于Lambert′s law假設入射角與單位面積散射量的關系遵循余弦定律[37，38]，通過對鄰近期VH極化影像歸一化近似獲得VH影像，計算公式為：

(12)

2.4 模型精度評價與融合方法

利用均方根誤差(RMSE)、偏差(Bias)評價土壤含水率反演精度，計算公式為：

(13)

(14)

式中：Obsi、Simi分別為i樣本的觀測值和模型反演值；N為樣本數。

(15)

式中：i表示期數，本研究選取0720期、0805期、0824期反演結果進行融合。

3 結果與分析

3.1 去除植被影響的效果分析

計算各期影像各地類去除植被影響前后的后向散射系數差值(圖2)和土壤含水率反演結果差值(表2)，可以看出，各地類后向散射系數差值多在0.07 dB以下，土壤含水率差值多在0.002 cm3/cm3以下，說明研究區植被影響不大。結合植被覆蓋看，研究時段內研究區NDVI<0.23的面積平均占比約98%，以7月19日為例，沿河樹林與灌叢、灌叢與稀疏草地以及灌叢與草地3種混合地類的NDVI均值分別為0.192、0.191、0.186，雖然沿河樹林與灌叢NDVI最大值為0.24，但該地類像元數量少，且NDVI高值占比少(NDVI>0.2的像元占比約20%)。考慮研究區植被覆蓋度整體較低，導致去除植被影響后反演精度并未顯著提升。

圖2 各地類去除植被前后的后向散射系數差值Fig.2 Difference between backscatter coefficients after and before removing vegetation for various land cover

表2 各地類去除植被前后的土壤含水率反演結果差值Table 2 Difference between soil moisture inversion results after and before removing vegetation for various land cover cm3/cm3

3.2 0～5 cm土壤含水率反演結果

基于去除植被影響后反演的土壤含水率(圖3)進行分析，可以看出，IEM模型的反演值主要在0.1 cm3/cm3以下，時空分布較一致，其中0714期與0818期的土壤含水率較低；Oh模型的反演值多為0.1～0.2 cm3/cm3，0720期、0805期和0824期土壤含水率的空間分布總體相似，其中0805期的土壤含水率略高于其他兩期，0714期與0818期的空間分布不連貫。從反演方法看，基于IEM模型反演的0714期、0818期地表粗糙度反演公式有別于其他3期，基于Oh模型反演的0714期、0818期VH極化數據由相鄰極化數據進行歸一化處理近似獲得，且土壤含水率的反演公式也有別于其他3期，由此導致即使利用單一反演模型也存在反演方法差異；而基于同一模型同一反演方法的土壤含水率時空分布一致性強(0720期、0805期、0824期)，反演精度差異小。

圖3 基于IEM和Oh模型反演的土壤含水率分布Fig.3 Soil moisture retrieved from ASAR images using IEM and Oh model

由同一模型相鄰兩期的土壤含水率差值(圖4)可知，0824-0818期差值主要表現為土壤含水率減少，基于IEM、Oh模型反演的土壤含水率減少區域占比分別約為62%、45%，增加區域占比分別約為37%、35%；基于IEM、Oh模型的0818-0805期差值減少區域占比均約為40%，增加區域占比分別為59%、37%；0805-0720期土壤含水率增加區域占比分別約為44%、61%，0720-0714期土壤含水率減少區域占比分別約為63%、56%。可見IEM和Oh模型反演的相鄰兩期土壤含水率的主要變化一致。

圖4 基于IEM和Oh模型反演的土壤含水率差值Fig.4 Difference of soil moisture (Δmv) between adjacent dates retrieved by IEM and Oh model

從空間分布上看，Oh模型反演異常值主要位于西南和東北區域，IEM模型反演異常值主要位于西南區域。結合流域植被看，土壤含水率異常值區域與流域植被覆蓋度相對較高區域并不對應；結合地形因素看，IEM、Oh模型反演異常值區域對應影像后向散射系數極亮值區，這主要源于地形(坡度>15°)導致的成像高亮，雖然對影像進行了地形校正，但在海拔高且坡度大的區域效果有限。

3.3 實測數據驗證與融合分析

利用與影像獲取時間一致的19個站點的實測土壤含水率數據進行驗證(0714期實測數據缺失)(表3)，可以看出，IEM模型反演的土壤含水率RMSE和Bias值均低于Oh模型，R值均高于Oh模型，IEM模型反演的0805期土壤含水率存在一定程度的低估，而Oh模型反演的土壤含水率則普遍存在高估。

表3 基于IEM、Oh模型反演的土壤含水率精度對比Table 3 Comparison of accuracies of soil moisture retrieved by IEM and Oh model cm3/cm3

由于0714期缺實測驗證數據，0818期IEM模型反演結果精度高但Oh模型反演精度低，因此只對其他3期結果進行融合。由融合后的統計結果(表4)可知，0720期、0805期和0824期土壤含水率RMSE相比融合前減少了0.003～0.065 cm3/cm3；融合后Bias相對IEM增加、相對Oh減少；融合后R值均有不同程度的提高。其中，加權平均法的RMSE、Bias均小于算術平均法，R值均大于算術平均法，總體上結合實測數據的加權平均法優于算術平均法。以均值融合結果為例(圖5)，融合后的土壤含水率在不同時間上的空間連續性較好。受影像期數、間隔時間、特征限制，基于單個模型不同反演方法獲取土壤含水率的能力依然有限(如0805期)，融合方法可視為獲取更高精度土壤含水率的一種選擇。

表4 融合后的土壤含水率精度對比Table 4 Comparison of accuracies of retrieved soil moisture after fusion cm3/cm3

圖5 融合后的土壤含水率分布Fig.5 Distribution of retrieved soil moisture after fusion

4 結論

本文基于水云模型分析植被后向散射系數的影響，并分別利用IEM和 Oh模型建立適合不同時相SAR影像的多反演方法，進一步基于兩種模型反演結果建立土壤含水率融合方法，結合實測數據分析融合前后的土壤含水率精度，可為獲取高精度土壤含水率的多模型/方法反演、融合研究提供參考。主要結論如下：1)研究期內Walnut Gulch流域NDVI<0.23的面積占比達98%，植被覆蓋度整體較低，各地類植被去除前后的土壤含水率反演差值多期均值小于0.002 cm3/cm3，去除植被影響對反演結果的作用不明顯。2)基于IEM、Oh模型反演的0720期、0805期、0824期土壤含水率的時空分布較一致，0714期、0818期由于地表粗糙度和土壤含水率反演公式不同，導致其時空分布不同；基于同模型同一反演方法的土壤含水率時空分布的一致性強，反演精度差異小，Oh模型在影像數據未滿足其設定條件下的反演效果較差。3)實測數據驗證表明，基于IEM模型反演的土壤含水率的RMSE與Bias低于Oh模型，其中 0805期的土壤含水率存在低估，而Oh模型普遍存在高估；均值融合、加權融合方法能一定程度克服使用單模型的缺點，融合后土壤含水率在不同時間上的空間連續性較好，融合后的土壤含水率RMSE減小0.003～0.065 cm3/cm3，與實測站點的相關系數均有不同程度提高。本文所用融合方法可視為獲取多時相高精度土壤含水率的有效方式。

本文研究區為植被覆蓋度低的半干旱區，植被對反演結果的影響并不明顯；基于IEM模型的土壤含水率經驗反演方法和地表粗糙度等參數的反演方式在其他區域的適用性有待驗證；滿足Oh模型條件的最佳極化影像方法有待改進；SAR影像在地形陡峭地區的成像異常問題也有待解決。Sentinel SAR作為 ENVISAT-ASAR后續影像，時空分辨率高，下一步將利用該數據在國內高植被覆蓋區深入開展多時相SAR的多模型反演方法與融合方法的適用性研究。