基于區域特征相似度的微波土壤水分反演結果可信度評價

陳魯皖， 韓 玲， 秦小寶， 王文娟

(長安大學 地質工程與測繪學院， 陜西 西安 710064)

1 研究背景

土壤水分是農業、林業、氣象學、水文學和生態學等領域研究的重要參數[1]，利用微波遙感反演土壤水分是一種很有效的手段。利用SAR反演土壤水分的原理是土壤表面的后向散射系數與土壤含水量、地表粗糙度、地表物質介電特性、土壤物理特性、植被特性以及雷達系統參數密切相關[2]。目前通常利用野外地表實測數據結合同期SAR數據來擬合或訓練土壤水分反演經驗方程(關系)，以此得到更高精度的土壤水分反演結果[3-5]。反演結果的精度評價一般是基于實測數據的，即選取部分采樣點，將采樣點的土壤水分實測數據和反演數據進行比較，通常選擇均方根誤差(RMSE)和決定系數(R2)作為精度指標[6-9]。

野外實測布設樣區時，受多種因素的影響(衛星過頂時間、測量人員數目和測量經費等)，往往樣方的設置范圍和采樣點的數量都比較有限。根據主動微波探測土壤水分的物理機制，影響地表后向散射系數的各因子在不同地區數值不同，即反演區域內存在各影響因子的異質性。所以有限范圍的樣區中實測數據擬合或訓練得到的土壤水分反演經驗方程(關系)，在較大范圍的區域進行反演，不同地區的反演結果精度必然不同。目前常用的基于實測數據的土壤水分反演精度評價方法實際上只能反映實測數據所在樣區的反演精度，不能反映所有反演區域的反演精度。

本文提出了一種基于區域特征相似度的微波土壤水分反演結果可信度評價方法，該方法可以有效反映所有反演區域的土壤水分反演精度。這里的可信度，是反演區域內土壤水分反演結果的可信程度，表示了反演結果與土壤水分真實值之間的接近程度。本文方法的原理是，反演區域中不同地區與樣區在土壤濕度分布、地表粗糙度、地表物質介電特性、土壤物理特性、植被特性和雷達系統參數等方面具有不同的相似性，而土壤水分反演經驗方程(關系)又是基于樣區的實測數據得到的，因此可以通過不同反演地區與樣區的相似程度來反映土壤水分反演結果的可信度。通過分析TM影像、DEM數據、HWSD土壤質地數據集和EnvisatASAR數據，結合實際地表特征，選擇土壤濕度、地表溫度、NDVI、STI(SoilTextureIndex，土壤質地指數)、地形指數、雷達入射角和LandsatTM影像的4個波段(b3、b4、b5、b7)共10個影響因子，通過主成分分析得到前3個主成分，將前3個主成分合成為RGB影像，再基于該RGB影像使用分水嶺算法進行分割，構建各分割區域的特征向量，計算分割后各區域與樣本區域的馬氏距離，基于特征相似度數據集計算反映土壤水分反演結果的可信度。

2 試驗數據

本文所使用的試驗數據來源于中國科學院寒區旱區科學數據中心所提供的“黑河綜合遙感聯合試驗”。

2.1 SAR數據來源

本文選用歐空局的Envisat-1衛星上ASAR傳感器獲取的SAR影像作為土壤含水量反演的數據源。采用了2008年7月11日的ASAR數據，該影像的入射波段為C波段(f=5.331GHZ)，入射角范圍31.0°～36.3°，經度范圍99°28′E～100°43′E，緯度范圍38°42′N～39°48′N，地面分辨率為12.5m×12.5m，工作模式為AlternatingPolarization，極化方式為VV和VH兩種。

2.2 地面數據來源

為驗證本文所提出的方法的有效性，特地選取了具有多個采樣區的甘肅省黑河中游臨澤地區作為研究區。該地區有4個樣區，分別為360m×360m的樣方B、C、D、E，樣點間距為60m，B和C樣方地表類型是帶稀疏雜草的鹽堿地，D樣方地表類型是苜蓿，E樣方地表類型是大麥地。與ASAR影像同期的共187組地面樣方觀測數據，觀測數據包括地表土壤含水量、均方根高度和相關長度。為了對比驗證區域特征相似度是否能夠有效反映反演結果的可信度，分別基于B樣區和E樣區中的部分實測數據構建土壤水分反演經驗方程，其余的B、C、D、E樣區中的實測數據為驗證數據。

2.3 其他數據來源

本文選用了Landsat5衛星的TM影像，空間分辨率為30m×30m；還選用了中國科學院寒區旱區科學數據中心提供的“黑河流域HWSD土壤質地數據集”和“黑河流域ASTERGDEM數據集”。

3 研究方法

土壤水分反演結果可信度評價方法的流程如圖1所示。

圖1 研究方法流程圖

具體步驟：(1)利用多元數據獲取10個與地表土壤后向散射系數密切相關的影響因子；(2)將各影響因子的影像進行波段合成；(3)利用主成分分析法[10](PCA)對步驟(2)的結果圖像處理提取出主成分，將前3個主成分合成為RGB影像；(4)利用分水嶺算法[11]對步驟(3)的結果影像做分割，獲得過分割影像；(5)構建各分割區域和樣區的特征向量，分別計算不同分割區域與各樣本區域間的特征相似度；(6)對步驟(5)的結果進行處理，得到反演結果的可信度。

3.1 獲取影響因子

基于SAR反演土壤水分的物理機制，本文選擇了土壤濕度、地表溫度、NDVI、STI、地形指數、雷達入射角這6種影響因子。

TM影像4個波段(b3、b4、b5、b7)中，通過構建各波段之間二維光譜特征空間[12]，發現TM紅外、近紅外、短波紅外波段的光譜特征與土壤水分存在復雜而緊密的聯系，因而可以用來監測地表土壤濕度。

NDVI、地表溫度和土壤濕度可由TM反演計算得到；STI為本文提出的一種反映土壤質地的新指數，由“黑河流域HWSD土壤質地數據集” 中的土壤砂土、黏土、淤泥含量計算得到；地形指數由“黑河流域ASTERGDEM數據集”通過提取地形的坡度計算得到；雷達入射角由輻射定標后的SAR影像中得到。3.1.1 地表溫度和土壤濕度 利用TM影像和單窗算法可以反演地表溫度[13-14]。統計分析后發現土壤水分含量與其溫度之間存在相關性較好的線性關系，利用線性方程(1)可直接反演土壤濕度的分布。

Mv=cT+d

(1)

式中：Mv為土壤水分，cm3/cm3； T為地表土壤溫度，℃； c、d為方程的經驗系數。

用通過反演得到的地表土壤溫度來估算土壤濕度，可以反映出土壤濕度的分布規律。

3.1.2 土壤質地指數(SoilTextureIndex,STI) 土壤質地是根據土壤的顆粒組成劃分的土壤類型，通常使用砂土含量、黏土含量和淤泥含量3種指標來量化土壤質地。為綜合考慮這3種指標，減少參與后續主成分分析的因子數量，本文提出了土壤質地指數STI。

本文使用AIEM模型建立基于不同土壤砂土含量sv和黏土含量cv的后向散射系數模擬數據集，其中sv和cv的取值參照“黑河流域HWSD土壤質地數據集”中提供的黑河地區的sv和cv的數值。AIEM模型的地表參數和雷達系統參數為固定值，做如下設定：雷達入射角給定為30.0°，入射頻率5.331GHz，土壤水分Mv=30%，均方根高度S=1.2cm，相關長度l=12cm(粗糙度參數通過統計分析“黑河綜合遙感聯合試驗”、美國NASA與農業部的合作試驗“Washita92”可得)。

(2)

相關系數R2=0.999819，整理式(2)可得：

(3)

若設：

STI=3sv+cv

(4)

則可得后向散射系數與STI的線性方程：

(5)

式中：a、b為方程的經驗系數。

3.1.3 地形指數 地形指數是指反演區域地球表面積與投影面積之比，可通過坡度的余弦得到[15]。

Z=1/cos(S)

(6)

式中： Z為地形指數； S為坡度，rad。

3.1.4 雷達入射角 對分布目標而言，雷達回波在近距離較強，隨著向遠距離移動將逐漸減弱。雷達入射角可對微波遙感的反演結果產生較大影響[16]。而本文選取的EnvisatASAR影像，沿單幅圖像距離向入射角大約有5°的變化。因此，雷達入射角也是影響土壤后向散射系數的一個因素[2]。

本文通過使用NEST(NestESASARToolbox)軟件處理EnvisatASAR數據，可得到雷達入射角影像。

3.2 基于主成分分析的分水嶺算法影像分割

3.2.1 主成分分析 主成分分析法又稱K-L(Karhunen-Loeve)變換，可實現在盡可能不丟失信息的同時，用幾個綜合性分量代表多波段的原圖像[10]。一般情況下，前3個主成分(PC1、PC2、PC3)包含所有波段中90%以上的方差信息。

在對10個后向散射系數影響因子的圖層進行波段合成后，利用ENVI5.1軟件中的PrincipalComponents模塊進行主成分分析。

3.2.2 基于分水嶺算法的影像分割 分水嶺算法是一種建立在數學形態學基礎之上算法。傳統分水嶺算法只能用于灰度圖像的分割，無法分割彩色圖像。本文把主成分分析后得到的前3個主成分合成為RGB彩色影像，選擇基于形態學梯度的分水嶺彩色圖像分割方法，首先計算彩色圖像的梯度圖像，然后結合多尺度形態學和巴特沃斯低通濾波對彩色梯度進行修正得到新的彩色梯度圖，最后對修正后的梯度圖像進行分水嶺變換得到分割結果[11-17]。

3.3 計算區域特征相似度

在分割后得到的區域和樣區之間進行特征相似度的判斷，使用數值來有效地表示各分割區域和樣區間地表情況的相似程度，這便是區域特征相似度的計算問題。特征相似度大多可以表示成向量的形式，本文構建特征向量時選用了參與主成分分析的6個因子(土壤濕度、地表溫度、NDVI、STI、地形指數、雷達入射角)作為特征向量的分量。

通常采用距離法來表示特征的相似度量，常用的有歐氏距離、馬氏距離[18-19]。對后向散射系數有影響的因素中，土壤濕度、地形指數和NDVI值這幾個因子相對于其他因子對后向散射系數的影響更大[20]。馬氏距離主要用于特征向量的各分量具有相關性或各分量的權重不等的情況。

因此本文采用馬氏距離來計算區域特征相似度，馬氏距離[19]表示為：

(7)

式中：x和y為兩個準備計算的區域特征向量;E為影像的特征向量協方差矩陣;d(x,y)為向量x和y之間的馬氏距離。

3.4 計算反演結果可信度

反演結果可信度基于區域特征相似度進行計算，反映了土壤水分反演結果的可信程度，可通過下式計算：

(8)

式中：Re為反演結果可信度;d為馬氏距離，min為1/d中的最小值，max為1/d中的最大值。

3.5 土壤水分反演

根據文獻[3-5]，給定入射角時，裸土地表的后向散射系數與地表粗糙度、土壤含水量之間的關系可以表示為[4]：

(9)

Zs=S3/l2

(10)

4 試驗對比

4.1 獲取影響因子的試驗

本文使用Landsat5的TM影像，反演得到NDVI、地表溫度和土壤濕度影像。圖2(a)和2(b)分別是反演區域的地表溫度和土壤濕度影像。圖2(a)中紅色區域地表溫度最高，其土地覆蓋類型為裸土；其次是地表溫度較高的黃色區域，其土地覆蓋類型為稀疏草地；再次是溫度較低的綠色區域，土地覆蓋類型為小麥地、大麥地、玉米地和濕地。土壤濕度圖也有類似的分布情況。

根據“黑河流域HWSD土壤質地數據集”得到的研究區土壤砂土含量、黏土含量，采用公式(4)計算土壤質地指數STI。地形指數由“黑河流域ASTERGDEM數據集”得到；入射角由輻射定標后的SAR影像中得到。

4.2 基于主成分分析的分水嶺分割試驗

對10個影響因子進行主成分分析，得到表1。

從表1可以看出，前3個主成分的樣本方差累計貢獻率已經達到96.65%，把前3個主成分合成為RGB彩色圖。使用分水嶺算法對RGB影像進行分割，圖3為分割后的結果圖。從圖3中可以看出，裸地、草地、農田和水體都得到了較好的同質分割。

表1 主成分因子貢獻率

4.3 反演結果可信度計算

選擇臨澤地區的B樣區，計算每個分割區域特征向量與B樣區特征向量之間的馬氏距離，然后根據公式(9)計算得到反演結果可信度。為驗證反演結果可信度的有效性，對比基于不同樣區的實測數據得到的反演結果，又選擇E樣區，同樣計算每個分割區域特征向量與E樣區特征向量之間的馬氏距離和反演結果可信度。表2為部分分割區域與B、E樣區之間的馬氏距離和反演結果可信度，表2中各分割區域的地物類型參考“黑河生態水文遙感試驗：黑河流域土地利用覆被數據集(2011年7月)”。

表2 部分分割區域與樣本區域之間的特征相似度

B樣方地表類型是鹽堿草地，表2中與B樣方相似程度較高的區域為F4、F43、F5與F11。E樣方地表類型是大麥地，表2中與E樣方較為相似的區域為F24、F40、F22、F48。

各分割區域中，相似度較好分割區域的土地覆蓋類型基本與B樣方一致，這從一定程度上證明了本文方法的有效性。

為評估本文提出方法的有效性，采用包括臨澤地區4個樣方共187個樣點的地面觀測數據參與驗證試驗。選擇B樣方中的部分地面觀測數據作為訓練數據結合同期EnvisatASAR觀測影像，得到土壤水分反演經驗方程。

將B樣區訓練數據包括土壤含水量、均方根高度和相關長度等數據，代入到公式(9)、(10)中，可得到公式(11)、(12)，兩式聯立消參求解可得反演區域的土壤含水量。

0.96lnZslnMv+5.31

(11)

0.22lnZslnMv-12.86

(12)

圖4(a)為式(11)、(12)的土壤水分反演結果，圖4(b)為反演結果可信度。

同理，選擇E樣方中的部分地面觀測數據作為訓練數據，可反演得到土壤水分。圖5(a)為土壤水分反演結果，圖5(b)為反演結果可信度。

4.4 對比分析

首先將B樣方的其他土壤水分實測數據以及C、D、E樣方的土壤水分實測數據作為驗證數據，與基于樣區B的土壤水分反演結果中各采樣點的土壤水分反演數據進行對比；同樣地對基于樣區E的土壤水分反演結果進行對比，表3為對比的結果。

圖2 部分試驗影像及TM影像

圖3 RGB合成圖分割結果

訓練樣區驗證樣區R2分割區域特征相似度ReB0.7956F40.01321.0000BC0.7024F170.46130.3924D0.6121F221.06070.2939E0.4142F242.10760.1881B0.4211F42.10760.3694EC0.3875F172.55890.2812D0.6201F221.04830.4728E0.7442F240.02371.0000

圖4 土壤水分反演結果(基于樣區B)

圖5 土壤水分反演結果(基于樣區E)

對比的反演結果精度表明，反演結果可信度與R2正相關，隨著Re的增大，R2也隨之增大。由此可以證明，本文提出的基于區域特征相似度的微波土壤水分反演結果可信度評價方法可以有效地反映整個反演區域的土壤水分反演結果的精度。

5 結 論

為解決土壤水分反演結果的精度評價問題，本文提出了一種基于區域相似度的土壤水分反演結果可信度評價的方法。

(1)基于主成分分析的分水嶺算法影像分割的結果較好，說明本文選取的10個參與主成分分析的影響因子較為合適；

(2)本文提出的土壤質地指數(STI)，有效地減少了參與主成分分析的影響因子數，能夠綜合反映土壤砂土含量、黏土含量和淤泥含量3個土壤質地指標；

(3)通過計算分割區域的特征向量與樣區特征向量的馬氏距離，得到區域特征相似度，根據區域特征相似度計算反演結果可信度，試驗結果證明，本文提出的方法是有效可行的。

致謝：感謝中國科學院寒區旱區科學數據中心(http://westdc.westgis.ac.cn)所提供的“黑河綜合遙感聯合試驗”EnvisatASAR數據和同步野外觀測數據。

[1]ENGMANET，CHAUHANN.Statusofmicrowavesoilmoisturemeasurementswithremotesensing[J].RemoteSensingofEnvironment, 1995,51(1):189-198.

[2] 趙英時.遙感應用分析原理與方法[M].北京：科學出版社，2013.

[3] 余 凡，趙英時.合成孔徑雷達反演裸露地表土壤水分的新方法[J].武漢大學學報(信息科學版),2010，35(3)：317-321.

[4] 陳 晶，賈 毅，余 凡．雙極化雷達反演裸露地表土壤水分[J]．農業工程學報，2013，29(10)：109-115+298.

[5] 甄珮珮. 基于粗糙度參數的風沙灘地區土壤水分微波遙感反演模型研究[D]. 西安：長安大學,2016.

[6] 張桂欣，郝振純，祝善友，等.AMSR2 缺失數據重建及其土壤濕度反演精度評價[J]. 農業工程學報，2016,32(20)：137-143.

[7] 王 東.基于SAR圖像的植被覆蓋下土壤含水量反演方法研究[D].成都：電子科技大學，2014.

[8] 王軍戰，張友靜，鮑艷松，等. 基于ASAR雙極化雷達數據的半經驗模型反演土壤濕度[J].地理與地理信息科學，2009,25(2)：5-9+2.

[9] 孔金玲，甄珮珮,李菁菁.基于新的組合粗糙度參數的土壤水分微波遙感反演[J]. 地理與地理信息科學，2016,32(3)：34-38.

[10] 鄧書斌.ENVI遙感圖像處理方法[M].北京：科學出版社，2014.

[11] 徐天芝. 基于形態學梯度的分水嶺彩色圖像分割[J]. 計算機工程與應用，2016,52(11):200-203+208.

[12] 劉 英,吳立新,馬保東.基于TM/ETM+光譜特征空間的土壤濕度遙感監測[J].中國礦業大學學報，2013.42(2)：296-301.

[13] 覃志豪，ZHANGMinghua,KARNIELIA.用陸地衛星TM6數據演算地表溫度的單窗算法[J].地理學報，2001,56(4)：456-466.

[14] 喬平林,張繼賢,燕 琴，等.利用TM6進行土壤水分的監測研究[J].測繪通報，2003(7)：14-15+18.

[15] 洪奕豐,嚴恩萍,林 輝.浙東地貌形態與土地覆被格局關系的研究[J].中南林業科技大學學報，2012,32(3)：63-69.

[16] 王建社.SAR圖像入射角效應校正與分割方法研究[D]. 合肥：合肥工業大學，2011.

[17] 張聰聰. 基于改進分水嶺算法和NCut算法的彩色圖像分割研究[D]. 廣州：華南理工大學，2016.

[18] 莊越挺,潘云鶴,吳 飛.網上多媒體信息分析與檢索[M].北京：清華大學出版社，2002．

[19]XIANGShiming,NIEFeiping.LearningaMahalanobisdistancemetricfordataclusteringandclassification[J].PatternRecognition,2008,41 (12):3600-3612．

[20]DOBSONMC，ULABYFT，HALLIKAINENMT，etal.Microwavedielectricbehaviorofwetsoil[J].IEEETransactionsonGeoscienceandRemoteSensing,1985,23(1):35-46.