黃 對(duì),劉 九 夫，張 建 云,王 文,崔 巍

(1.南京水利科學(xué)研究院水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029；2.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098)

0 引言

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

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

1 研究區(qū)域及數(shù)據(jù)

1.1 研究區(qū)域

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

1.2 研究數(shù)據(jù)

(1)土壤含水率數(shù)據(jù)。研究區(qū)共有19個(gè)土壤含水率觀(guān)測(cè)站點(diǎn)(圖1)，每20 min記錄0～5 cm的土壤體積含水率等參數(shù)，已經(jīng)過(guò)鹽分等誤差校正，采樣時(shí)間為2004年7月1日-9月30日[30]。

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

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

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

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

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模型的土壤含水率反演可分為經(jīng)驗(yàn)[35]和半經(jīng)驗(yàn)方法[18]，二者均基于模型模擬數(shù)據(jù)的散射特征分析構(gòu)建反演公式，區(qū)別在于構(gòu)建時(shí)是否利用實(shí)測(cè)土壤含水率[23]。

(7)

(8)

式中：p為VH極化方式；A、B、C為對(duì)應(yīng)入射角與極化方式下的系數(shù)。

(9)

2.3 Oh模型

(10)

(11)

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

(12)

2.4 模型精度評(píng)價(jià)與融合方法

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

(13)

(14)

式中：Obsi、Simi分別為i樣本的觀(guān)測(cè)值和模型反演值；N為樣本數(shù)。

(15)

式中：i表示期數(shù)，本研究選取0720期、0805期、0824期反演結(jié)果進(jìn)行融合。

3 結(jié)果與分析

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

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

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

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

3.2 0～5 cm土壤含水率反演結(jié)果

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

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

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

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

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

3.3 實(shí)測(cè)數(shù)據(jù)驗(yàn)證與融合分析

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

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

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

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

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

4 結(jié)論

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

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