基于高光譜影像分解的土壤含水量反演技術

吳 見，劉民士，李偉濤

（滁州學院 地理信息與旅游學院，安徽 滁州239000）

土壤是農業生產的基礎，土壤的水分含量、質地、有機質含量［1－2］等性質會影響作物生長，其中，土壤水分不僅是植物生長和發育的必要條件，也是精準農業中重要的信息。因此，土壤含水量監測一直是學者們研究的重點［3－4］。傳統的土壤含水量測定方法，如張力計法、烘干法等均以點測量為基礎，準確性較高，但周期長、工作量大。遙感技術主要以可見光—近紅外、熱紅外和微波波段對土壤含水量進行反演，熱紅外波段反映的是地表很薄的一層信息，但也可利用熱慣量和輻射平衡估測表層以下的土壤水分，高光譜遙感具有精細的光譜分辨率，能根據可見光和近紅外波段的光譜信息，分析影響土壤含水量變化的光譜特征［5－6］，建立光譜特征與 土壤含水量的 關系模型［7－8］，進而定量反演土壤含水量。

目前，土壤含水量遙感監測研究大多是針對裸土進行的，對于衛星遙感來說，實際情況往往是植被與土壤同時存在于像元中，土壤含水量的監測不可避免地會受到植被光譜信息的影響。實際研究發現土壤含水量變化引起的光譜差異與植被含水量變化引起的光譜差異并不相同［9］，因此剔除植被光譜信息對土壤水分的干擾，是監測土壤含水量的有效途徑。本研究對高光譜影像進行光譜分解，剔除植被光譜的干擾，同時對土壤貢獻的光譜進一步處理，建立土壤含水量反演模型。

1 研究區概況及影像數據選取

1.1 研究區概況

懷柔區地處北京東北部，總面積達2 128.7km2，是華北平原、內蒙古高原以及燕山山脈的過渡地帶。該區北部環山，南部是草原，地形包括深山、淺山、平原3種類型，從南向北綿延128km，且山地面積占整個懷柔區的88.7%，宜林山場森林面積達41%。地勢呈北高南低，海拔在34～1 661m，變化較大。該區氣候冬天受西伯利亞冷空氣影響，寒冷干燥，夏天受海洋氣團影響，溫熱濕潤，春秋季節短且干旱多風。年平均氣溫在6～12℃，常年平均降水量在470～850mm。土壤主要包括草墊土、褐土、棕壤、風砂土等類型。

1.2 數據獲取

本研究選取北京懷柔部分地區2001年5月19日的EO-1Hyperion高光譜數據，該數據共有242波段，其中1～70波段是可見光近紅外波段（VNIR），71～242波段是短波紅外波段（SWIR），光譜分辨率為10nm，空間分辨率為30m。

首先刪除2個重復，20個受水汽影響嚴重及44個未定標波段，剩余176個波段；然后對剩余波段進行處理，包括壞線修復、條紋去除以及smile效應去除；最后，對處理后的圖像進行檢驗，繼續刪除質量差的波段7個，剩余169個波段。被刪除的具體波段為1～7，58～78，121～129，166～180，185～186，224～242。利用FLAASH軟件對剩余的169個波段進行大氣糾正，得到反射率圖像。大氣糾正后，采用1∶50 000地形圖對影像進行幾何糾正，總誤差是0.35個像元。

2 基于植被光譜信息剔除的土壤含水量反演

2.1 基于光譜匹配的分解算法

目前，全受限的線性光譜混合分解模型是最為常用的方法，該模型中像元的光譜表達成各基本組分光譜及其面積比率的線性組合［10］，即：

式中：Rb——b波段的光譜信息；N——基本組分數目；fi——基本組分i占所在像元的面積比例；Ri，b——基本組分i在b波段的光譜信息；eb——b波段的誤差。應用該模型的關鍵是盡可能減少模型中每個像元的誤差，即使均方根誤差RMSE最小化［10］：

式中：M——波段數。解該模型最常用的方法是最小二乘法，即求得式（2）RMSE最小的解。

將影像或光譜庫中選取的基本組分光譜（r1，r2，…，rn）根據已設定的比例（f1，f2，…，fn）組合，獲取若干測試光譜S1＝f1r1＋f2r2＋…＋fnrn，然后將待分解的像元光譜，即目標光譜S2與測試光譜S1進行匹配［10］，把匹配效果最佳的測試光譜的基本組分比例作為影像分解結果。因此，可將混合像元分解視為帶約束的非線性最優問題［10］：

式中：G——目標函數，描述目標與測試光譜差異的匹配模型。由最小二乘法原理可定義G為測試與目標光譜間的歐式距離［10］：

因S1b－S2b＝eb，根據式（2）可知GLS＝ ■M ×RMSE，說明最小二乘法也屬于光譜匹配的非線性最優化算法，但GLS采用歐式距離描述光譜匹配程度時，易受局部奇異值等因素影響。基于此，本研究引用楊偉等［10］提出的基于相關系數匹配的分解算法，設目標函數GSCM為：

式中：SCM——測試與目標光譜的相關系數。

采用相關系數SCM描述目標與測試光譜之間的匹配程度，以SCM最大值的測試光譜作為混合像元分解的結果。

2.2 土壤光譜信息處理

高光譜混合像元分解后，得到土壤和植被比例，以及裸土和純植被基本組分的光譜信息，因此可進一步剔除植被光譜的干擾，計算各像元中土壤貢獻的反射率信息。獲取土壤貢獻的反射率信息后，為突出土壤光譜的吸收特征，本研究采用光譜包絡線去除［11］和一階微分［12］兩種方法對光譜曲線做進一步處理。

2.2.1 反射率光譜信息 由于不同像元對應的土壤含水量不同，在不同含水量條件下裸土基本組分光譜是有差異的，因此不能采用土壤比例信息與裸土基本組分光譜直接計算不同比例土壤貢獻的反射率信息。由于研究區內植被類型較為單一，因此設定研究區域純植被基本組分光譜不變，從原始影像混合光譜中將植被光譜剔除，進而獲取土壤貢獻的反射率光譜信息，由像元二分模型［13］：

可知，剔除植被光譜后的土壤光譜信息為：

式中：Rsb，Rv——土壤和植被在波段b貢獻的光譜信息；fs，fv——土壤和植被在像元中的比例。

2.2.2 光譜一階微分 一階微分法是常用的光譜處理技術，不僅能夠識別重疊光譜，增強曲線坡度上的光譜細微變化，使光譜吸收峰參數便于提取，而且可以有效消除系統誤差，降低大氣吸收、散射、輻射等噪聲的影響［12］。一階微分公式為［12］：

式中：λb——波段b；FDRλb——波段b和b＋1之間的一階微分值；Rλb，Rλb＋1——波段b和b＋1的光譜反射率；Δλ——波段b和b＋1的波長之差。

2.2.3 包絡線去除法 包絡線去除法能夠增強光譜曲線的反射、吸收特征，且反射率歸一化在0～1.0之間，是一種有效突出目標吸收特征的光譜分析法，同時能夠將光譜吸收特征歸一化在同一背景上，方便不同光譜曲線之間的比較［11］。“包絡線”即以直線逐點連接光譜曲線上的峰值點，且峰值點上的折線外角大于180°。光譜去包絡是將原始曲線上的值除以包絡線上對應的光譜值，計算公式為［11］：

式中：λb——波段b；Rcb——波段b的包絡線去除值；Rb——波段b的反射率值；Rstart，Rend——光譜曲線起始、末端點的反射率值；λstart，λend——光譜曲線起始、末端點的波長；K——光譜曲線起始、末端波段間的斜率。

2.3 土壤含水量與光譜信息相關性分析

2.3.1 土壤含水量與反射率光譜相關性分析 從圖1可知，反射率光譜與土壤含水量在450～1 050nm的相關系數為正值，相關系數的峰值分別出現在波段456～465，496～505，786～795，1 004～1 013nm處，相關系數值分別為0.47，0.48，0.37，0.29；在1 050～2 380nm的相關系數為負值，相關系數的峰值分別出現在波段1 134～1 143，1 324～1 333，1 467～1 476，2 012～2 021，2 302～2 311nm處，相關系數值分別為－0.42，－0.38，－0.51，－0.61，－0.56。為避免“過度擬合”現象及共線性問題，并非全部選取相關系數值最大的波段，而是根據光譜曲線特征劃分區域，選取各區域中具有代表性或相關系數絕對值最大的波段作為敏感波段，即496～505，786～795，1 134～1 143，1 467～1 476，2 012～2 021nm處。

圖1 土壤含水量與反射率光譜相關系數

2.3.2 土壤含水量與一階微分光譜相關性分析 一階微分光譜與土壤含水量的相關系數在部分波段比反射率光譜明顯增加，部分波段相對減小，且一階微分光譜與土壤含水量之間相關性時正時負，不呈單一相關（如圖2所示）。相關系數的正峰值分別出現在波段1 004～1 013，1 174～1 183，1 194～1 203，1 467～1 476，1 497～1 506，1 787～1 796，1 797～1 806nm處，相關系數值分別為0.37，0.52，0.39，0.34，0.23，0.76，0.69；相關系數的負峰值分別出現在波段1 064～1 073，1 264～1 273，1 324～1 333，1 527～1 536，1 952～1 961，1 972～1 981nm處，相關系數值分別為－0.25，－0.48，－0.40，－0.25，－0.58，－0.47。選取各區域中具有代表性或相關系數絕對值最大的波段作為敏感波段，即1 004～1 013，1 174～1 183，1 324～1 333，1 787～1 796，1 952～1 961nm處。

2.3.3 土壤含水量與包絡線去除光譜相關性分析從圖3可以看出，包絡線去除光譜與土壤含水量在1 560～1 805nm為負相關，其余波長范圍內為正相關。相關系數的峰值分別出現在波段656～665，736～745，1 014～1 023，1 274～1 283，1 467～1 476，1 517～1 526，2 062～2 071，2 242～2 251nm處，相關系數值分別為0.66，0.67，0.66，0.62，0.26，0.26，0.54，0.49。選取各區域中具有代表性或相關系數絕對值最大的波段作為敏感波段，即656～665，1 014～1 023，1 517～1 526，2 062～2 071，2 242～2 251nm波段。

圖2 土壤含水量與一階微分光譜相關系數

圖3 土壤含水量與包絡線去除光譜相關系數

2.4 土壤含水量反演模型建立

選取反射率、一階微分、包絡線去除光譜與土壤含水量敏感的波段進行逐步回歸分析。通過比較，建立土壤含水量反演多項式（表1—3），選擇R2值最大的反演模型，作為土壤光譜與土壤含水量的關系模型。

表1 反射率光譜與土壤含水量逐步回歸分析模型

表2 一階微分光譜與土壤含水量逐步回歸分析模型

表3 包絡線去除光譜與土壤含水量逐步回歸分析模型

由表1—3可知，預測土壤含水量的最佳模型是以波段X661，X1019和X2067的土壤包絡線去除光譜為自變量建立的回歸方程：

其預測R2值為0.85。采用實測含水量數據與式（12）得到的預測值進行擬合，檢驗結果的R2值為0.80（如圖4所示）。

圖4 土壤含水量預測與真實值比較

若不對原始影像進行分解，直接采用未剔除植被光譜的像元反射率、一階微分、包絡線去除光譜建立土壤含水量估測模型，得到的最佳模型是以波段X541，X979和X1632的一階微分光譜為自變量建立的回歸方程：

其預測R2值為0.36。采用實測含水量數據與式（13）得到的估測值進行擬合，檢驗結果的R2值僅為0.06。因此，剔除植被光譜干擾能夠更有效地反演土壤含水量。

3 結論

（1）將混合像元分解問題歸結為一個基于光譜匹配的非線性最優化問題，并針對最小二乘法的不足，引用了一種基于相關系數匹配的混合像元分解技術。高光譜數據本身提供了豐富的光譜信息，可以預見在高光譜數據中利用基于相關系數匹配的混合像元分解算法能夠得到更為精確的結果。

（2）剔除植被光譜干擾后，土壤反射率、一階微分、包絡線去除光譜的部分波段能夠敏感地反映土壤含水量的變化，最佳模型的估測R2值為0.85；若直接采用原始影像光譜建立估測模型，估測R2值僅為0.36。

（3）通過高光譜影像分解剔除植被光譜干擾估測土壤含水量的方法是可行的，可為今后遙感估測土壤含水量的研究提供參考。

［1］ 周萍，王潤生，閻柏琨，等.高光譜遙感土壤有機質信息提取研究［J］.地理科學進展，2008，27（5）：27-34.

［2］ 陳榮毅，張元明，潘伯榮，等.古爾班通古特沙漠土壤養分空間分異與干擾的關系［J］.中國沙漠，2007，27（2）：257-265.

［3］ 宋韜，鮑一丹，何勇.利用光譜數據快速檢測土壤含水量的方法研究［J］.光譜與光譜學分析，2009，29（3）：675-677.

［4］ Bedidi A，Cervelle B，Madeira J，et al.Moisture effects on visible spectral characteristics of lateritic soils［J］.Soil Science，1992（153）：129-141.

［5］ 孫建英，李民贊，唐寧，等.東北黑土的光譜特性及其與土壤參數的相關性分析［J］.光譜學與光譜分析，2007，27（8）：1502-1505.

［6］ 劉煥軍，張柏，宋開山，等.黑土土壤水分光譜響應特征與模型［J］.中國科學院研究生院學報，2008，25（4）：503-509.

［7］ 王昌佐，王紀華，王錦地，等.裸土表層含水量高光譜遙感的最佳波段選擇［J］.遙感信息，2003（4）：33-36.

［8］ 王靜，何挺，李玉環.基于高光譜遙感技術的土地質量信息挖掘研究［J］.遙感學報，2005，9（4）：438-445.

［9］ 吳代暉，范聞捷，崔要奎，等.高光譜遙感監測土壤含水量研究進展［J］.光譜學與光譜分析，2010，30（11）：3067-3071.

［10］ 楊偉，陳晉，松下文經，等.基于相關系數匹配的混合像元分解算法［J］.遙感學報，2008，12（3）：454-461.

［11］ 丁麗霞，王志輝，葛宏立.基于包絡線法的不同樹種葉片高光譜特征分析［J］.浙江林學院學報，2010，27（6）：809-814.

［12］ 尹業彪，李霞，趙釗，等.沙質土壤含水率高光譜預測模型建立及分析［J］.遙感技術與應用，2011，26（3）：355-359.

［13］ 牛寶茹，劉俊蓉，王政偉.干旱半干旱地區植被覆蓋度遙感信息提取研究［J］.武漢大學學報：信息科學版，2005，30（1）：27-30.