基于多源遙感協同反演的區域性土壤鹽漬化監測

馮雪力 劉全明

(1.內蒙古建筑職業技術學院，呼和浩特 010070； 2.內蒙古農業大學土木建筑工程學院，呼和浩特 010018)

0 引言

土壤鹽漬化是干旱、半干旱農業灌區面臨的主要土地退化問題。內蒙古河套灌區是西北地區典型的寒旱灌區，傳統灌溉多為含鹽黃河水漫灌，且受降水量少、蒸發量大、排水排鹽不暢、地勢平坦而地下水位居高等自然客觀條件限制，造成灌區鹽漬化和生態環境惡化，從而影響了農作物生長和糧食生產安全。因此，精準快速地獲取鹽漬土鹽分時空分布信息，對灌區的鹽漬土防治、中低產田改造與農業可持續發展具有重大的指導作用。傳統野外土壤定點調查鹽漬化方式不僅耗費時力、破壞性強，而且測點少，難以反映區域尺度的鹽漬化分布，而多傳感器、多波段、多時相的遙感技術為大面積實時動態監測鹽漬土狀況提供了可能。

國外從20 世紀70 年代[1-4]就開始研究鹽堿土遙感監測，國內從20世紀80 年代開始集中在利用可見光波段進行遙感鹽分反演。鹽漬化土壤光譜反射率和土壤參數間的轉換函數是復雜非線性的關系，故人工神經網絡逐漸被更多地用于獲取土壤鹽漬化參數[5-9]。吐爾遜·艾山[10]以新疆渭干河-庫車河三角洲綠洲作為研究區，利用實測光譜數據和地下水埋深、礦化度等因子構建了基于BP人工神經網絡(Back propagation artificial neural networks，BPANN)的土壤鹽分反演模型，模型精度為80.77%。王靜等[11]以松嫩平原西部長嶺縣為例，利用光譜導數變換，獲取了表征鹽漬土鹽分信息的最佳波段，所建立的3 層神經網絡鹽漬土鹽分模型預測結果優于傳統回歸模型。丁建麗等[12]研究發現，實測高光譜土壤、植被一階微分光譜變換對土壤鹽漬化響應度最敏感，基于實測綜合光譜指數的鹽漬化監測高光譜模型可以準確提取土壤鹽漬化信息，明顯優于傳統遙感方法中單純利用植被指數或者土壤鹽分指數的模型。姚遠等[13]以新疆渭干河-庫車河三角洲綠洲不同鹽漬化程度的土壤為研究對象，以電磁感應儀(EM38)測得的鹽漬土電導率數據和ASD光譜儀測得的鹽漬土高光譜數據為基礎數據源。王爽等[14]以新疆渭干河-庫車河三角洲綠洲區域范圍內不同鹽漬化程度的土壤高光譜反射率及其土壤含鹽量為基礎數據源，構建了地表光譜建模的最佳土壤鹽漬化監測模型。劉全明等[15]以受鹽漬化影響較嚴重的內蒙古河套灌區解放閘灌域為試驗區，構建了BPANN土壤含鹽量的定量反演模型。彭杰等[16]通過分析新疆、浙江、吉林3 個不同地區鹽漬化土壤的高光譜特征，研究了鹽漬化土壤高光譜特征的區域異質性，并構建高精度的跨區域土壤鹽分高光譜定量反演模型。彭杰等[17]以南疆地區水稻土為研究對象，通過分析土樣的高光譜數據和室內測定的鹽分與電導率數據，研究了耕作土壤含鹽量與電導率的關系，并比較了含鹽量和電導率與不同光譜指標的相關性以及二者高光譜反演的精度。研究者們使用高光譜和電磁感應數據進行土壤鹽分模擬，獲得了不少成果[18-22]。

在前期研究基礎上，本文優選高光譜土壤鹽分特征波段，融合微波后向散射特性，構建基于多源遙感數據的土壤鹽分預測模型，為利用多源遙感數據協同反演鹽漬土時空分布提供科學基礎。為此，在高光譜鹽分特征波段優選方法對比的基礎上，創建基于多源遙感數據的土壤鹽分經驗模型和人工智能非線性預測模型，從而將雷達后向散射系數、光譜反射率及其變換值等參數轉換為土壤鹽分，以提高鹽漬化監測的精度與廣適性。

1 研究區概況

內蒙古河套灌區地處內蒙古自治區的巴彥淖爾市，位于40°12′～41°20′N，106°10′～109°30′E。灌區地形平坦，西南高，東北低，海拔1 007～1 050 m，土壤以鹽漬化淺色草甸土和鹽土為主。夏季高溫干旱、冬季嚴寒少雪，年降水量130～215 mm，蒸發量高達2 100～2 300 mm，無霜期120～150 d，封凍期長，是典型的溫帶大陸性氣候。河套灌區東西長250 km，南北寬40～60 km，土地總面積18 890 km2，其中灌溉面積5 583 km2。河套灌區土壤養分含量地域分布極不均勻，降水量小、蒸發量大、土壤次生鹽漬化嚴重，快速了解灌區夏灌(4—6月)前土壤鹽分分布狀況對灌區灌溉制度和土壤改良有著深刻影響和重要意義。

2 材料與方法

2.1 數據獲取與處理

試驗區選在河套灌區內的解放閘灌域。為保證野外采樣時間與雷達成像時間對應一致，提前購置2016年4月6日RADARSAT-2精細四極化SLC(Single look complex)格式雷達影像一景，影像尺寸為25 km×25 km，其地面分辨率為8 m。考慮到土壤鹽漬化空間分布的不均勻性，在野外采樣前使用Google Earth對灌域的地形地貌進行分析，結合前期土壤鹽分統計資料并制定野外采樣路線，89個采樣點分布如圖1紅色五角星所示。野外通過手持GPS接收機測量采樣點在WGS84坐標系下的經緯度坐標，獲取采樣點地表粗糙度等數據，并同時在每個樣點采集3份土壤表層土，取平均值以消除取樣代表性誤差，把每個土壤樣本編號裝入鋁盒，采用干燥法測定土壤水溶性鹽分含量。同時使用帶有厘米格網的剖面板方法來測量地表粗糙度。

圖1 試驗區雷達影像圖Fig.1 Radar image of experimental area

試驗區所使用的RADARSAT-2影像數據是由C波段合成孔徑雷達系統獲取，其影像處理主要包括輻射定標、幾何校正、斜距轉地距和濾波等。本研究使用ENVI軟件的SAR Scape模塊對原始SLC影像進行多視、濾波、地理編碼和輻射定標等處理，在Google影像上選取地面控制點(Ground control point,GCP)對影像進行幾何配準后，得到標準四極化后向散射系數影像數據，表1所示為試驗區部分采樣點四極化后向散射系數與土樣全鹽量、地表組合粗糙度。表中SHH、SHV、SVH、SVV分別代表HH(水平極化方式發射，水平極化方式接收)、HV(水平極化方式發射，垂直極化方式接收)、VH(垂直極化方式發射，水平極化方式接收)、VV(垂直極化方式發射，垂直極化方式接收)4種極化方式的后向散射系數。

表1 雷達影像后向散射系數與土壤全鹽量原始數據Tab.1 Data of radar back-scattering coefficients and soil saltcontents

野外采樣點土壤反射率光譜數據的采集使用美國Analytical Spectral Device公司生產的Field Spec4便攜式光譜儀，可測量350～2 500 nm波段，其中，350～1 000 nm波段的光譜采樣間隔為1.4 nm，光譜分辨率為3 nm；1 000～2 500 nm波段的光譜采樣間隔為2 nm，光譜分辨率為10 nm。利用光譜儀自帶的ViewSpec Pro Version 6.0軟件將原始光譜反射率R進行了3種不同形式的變換：對數(lgR)、一階導數(R′)、二階導數(R″)，并將光譜反射率及上述3種變換光譜作為本研究的輸入光譜數據，后期分別將上述4種輸入值與對應采樣點鹽分值進行相關性分析。

2.2 主成分回歸

主成分回歸是考察多個變量間相關性的一種多元統計方法，研究如何通過少數幾個主成分來揭示多個變量間的內部結構，即從原始變量中導出少數幾個主成分，使它們盡可能多地保留原始變量的信息，且彼此間互不相關。通常數學上的處理就是將原來多個指標作線性組合，作為新的綜合指標。

2.3 逐步回歸分析

逐步回歸分析是將變量逐個引入模型，每引入一個解釋變量后都要進行F檢驗，并對已經選入的解釋變量逐個進行t檢驗，當原來引入的解釋變量由于后面解釋變量的引入變得不再顯著時，則將其刪除。以確保每次引入新的變量之前回歸方程中只包含顯著性變量。這是一個反復的過程，直到既沒有顯著的解釋變量選入回歸方程，也沒有不顯著的解釋變量從回歸方程中剔除為止，以保證最后所得到的解釋變量集是最優的。

2.4 偏最小二乘回歸

偏最小二乘回歸是一種多因變量對多自變量的回歸建模方法，可以較好地解決許多用普通多元回歸無法解決的問題。和PCR一樣，兩種方法都采用得分因子作為原始預測變量線性組合的依據，用于建立預測模型的得分因子之間必須線性無關。而PLSR與PCR的不同之處在于得分因子的提取方法不同，PCR產生的權重矩陣反映的是預測變量之間的協方差，PLSR產生的權重矩陣反映的是預測變量與響應變量之間的協方差。在建模當中，PLSR產生了權重矩陣，矩陣的列向量用于計算變量列向量的得分矩陣。通過不斷地迭代計算這些權重，使得響應與其相應的得分因子之間的協方差達到最大。PLSR是使用波段的權重來表示波段對于土壤實測值的敏感性，進而選取對鹽分敏感的特征波段，并進行土壤鹽分的回歸擬合。

2.5 BP人工神經網絡

BP人工神經網絡能夠辨識復雜系統輸入和輸出數據組間的非線性關系，不需任何假定或構建數學方程，可在一定程度上消除特異值的影響，具備強大的適應性、并行處理及抗干擾能力和突出的非線性擬合的特點，現已被各學科領域廣泛應用。本文經綜合優化試驗，確定選用3 層結構的BP模型，該模型由輸入層、隱含層和輸出層組成。考慮土壤含鹽量與高光譜光學特征波段及SLC影像的四極化后向散射系數、地表粗糙度有著顯著的響應關系，最終確定輸入層由11個神經元組成，包括4個光學特征波段反射率，6個后向散射系數及組合值SHH、SVV、SHV、SVH、SHH/SVV、SHV/SVH和地表組合粗糙度；隱含層神經元數的確定由試算法完成，本文經試算法確定隱含層神經元為10 個；輸出層神經元為1 個，對應為采樣點土樣的全鹽量。計算采用Matlab自帶神經網絡工具箱(nprtool)，采用原始數據的70%參與學習與訓練網絡模型，其余30%作為驗證。

3 特征波段選取建模與土壤鹽分預測

3.1 特征波段選取建模

圖2反映了土壤全鹽量與反射率的不同變換形式的相關關系。不難發現，相較于原始光譜和對數變換，光譜的一、二階導數可以獲得更好的相關性。在土壤全鹽量的相關性中，在一階導數變換的1 088～1 092 nm、1 806～1 810 nm、2 138～2 142 nm、2 289～2 295 nm 這4個波段表現最佳，相關系數分別為0.29、0.36、0.3、0.27，在二階導數變換的618～622 nm、1 802～1 806 nm、2 169～2 173 nm、2 344～2 348 nm這4個波段相關性最高，相關系數分別為0.37、0.28、0.39、0.27。

圖2 反射率的不同變換形式與土壤全鹽量的相關性分析Fig.2 Analysis on correlation of reflectance transforms and soil salinity contents

利用PCR對光譜的一、二階導數中的特征波段進行主成分分析，結果如圖3所示。可以發現，無論是一階還是二階導數變換的特征波段，都可以從中選取3特征波段來涵括，累積貢獻率分別為99.31%和99.87%；其中，第一主成分一階導數變換特征波段1 088～1 092 nm與二階導數變換特征波段618～622 nm的特征值和貢獻因子分別為4.07×10-7、1.51×10-7和72.19、74.08，皆具有較高代表性。

圖3 特征波段的主成分分析Fig.3 Principal component analysis of characteristic bands

利用MSR對光譜的一、二階導數中的特征波段進行逐步回歸分析，結果如圖4所示。其中圖4a中X1～X4分別表示1 088～1 092 nm、1 806～1 810 nm、2 138～2 142 nm和2 289～2 295 nm波段，圖4b中X1～X4分別表示618～622 nm、1 802～1 806 nm、2 169～2 173 nm、2 344～2 348 nm，研究表明，光譜一階導數的1 806～1 810 nm、2 138～2 142 nm 2個特征波段和二階導數的618～622 nm、1 802～1 806 nm、2 344～2 348 nm 3個特征波段更具有代表性。

圖4 特征波段的多元逐步回歸分析Fig.4 MSR analysis of characteristic bands

再使用PLSR選取特征波段并擬合土壤全鹽量，PLSR使用波段的權重來表示波段對于土壤實測值的敏感性。如圖5所示，PLSR篩選的波相較PCR以及MSR選取的波段延后，且存在極為敏感的波段，如一階導數變換中1 850～1 854 nm波段權重可達4，二階導數變換中1 853～1 857 nm與2 483～2 487 nm權重分別可達15和5.7，高權重波段與上文中分析的高相關性波段也均有延后現象。

圖5 反射率導數變換與土壤全鹽量的相關性及PLSR權重分析Fig.5 Analysis on correlation and PLSR weight between reflectance derivative transformation and soil salinity

PCR和MSR的擬合結果對比如表2所示。不難發現，二階導數變換的擬合程度高于一階導數變換，均方根誤差也都有所下降，說明二階導數變換更適合于土壤全鹽量的擬合；綜合比較使用光譜二階導數變換特征波段的MSR經驗模型模擬全鹽量效果最優，其決定系數(R2)和均方根誤差(RMSE)分別為0.343和9.414 g/kg。

表2 多模型擬合結果對比Tab.2 Results contrast of multi-model

圖6 土壤全鹽量實測值與模擬值Fig.6 Measured and simulated values of soil salinity

使用MSR選取的3個特征中心波段反射率、6個SAR后向散射系數及其組合值和地表組合粗糙度，建立土壤全鹽量與光譜二階導數的經驗模型

式中X618、X1 802、X2 346——618～622 nm、1 802～1 806 nm、2 344～2 348 nm波段的反射率二階導數算術平均值

同時建立神經網絡模型，由圖6可見，神經網絡模型R2為0.890 8，遠優于MSR模型，體現了神經網絡模型的優越性。

3.2 土壤鹽分預測

根據表3的河套灌區土壤鹽漬化分類標準，將采樣點所覆蓋的試驗區域分為5類鹽漬化類型，即非鹽漬土、輕度鹽漬土、中度鹽漬土、強鹽漬土和鹽堿土。圖7所示為MSR選取的光譜二階導數的3個特征波段和雷達參數建立的回歸方程模型模擬的全鹽量預測結果，先將MSR提取的3個特征光譜波段中心反射率二階導數和地表組合粗糙度用ArcGIS軟件地質統計模塊克立格插值生成5個柵格文件，同時提取生成6個SAR后向散射系數柵格文件，最后基于光譜二階導數的回歸方程，采用ENVI遙感圖像軟件進行11個柵格文件運算，其模擬結果的統計分析見表4，可見受到鹽漬化影響的面積占96.23%，鹽堿土占21.47%，應對重點區域采取科學的鹽堿土改良措施進行治理，以保證農業生產的可持續發展。

表3 土壤鹽漬化分類Tab.3 Classification of soil salinity

注：Ⅰ～Ⅴ級分別代表非鹽漬土、弱鹽漬土、中鹽漬土、強鹽漬土和鹽堿土。

圖7 研究區土壤全鹽量預測圖Fig.7 Prediction image of soil salinity in study area

表4 土壤鹽漬化預測分類占比Tab.4 Proportion of all levels of soil salinity prediction classification

4 結論

(1)通過不同形式的光譜變換處理可以使光譜原反射率與土壤鹽分含量的關系得以提升。分別利用PCR、MSR和PLSR選取特征光學波段，發現光譜的一、二階導數具有更好的相關性，尤其是二階導數變換的618～622 nm、1 802～1 806 nm、2 169～2 173 nm、2 344～2 348 nm這4個波段相關性較強；PLSR篩選的波段相較相關系數法偏后，其二階導數變換模型擬合度小于MSR。

(2)通過協同光譜特征波段反射率二階導數、雷達后向散射系數和地表組合粗糙度諸影響因子，對比土壤鹽分的PCR、MSR和PLSR模型，其中MSR模型的決定系數和均方根誤差分別為0.343和9.414 g/kg，優于另兩種模型，而BP人工神經網絡模型為最佳預測模型，其模型預測R2為0.890 8，模型的穩定性和預測精度較好。融合高光譜特征波段與SAR后向散射系數及地表組合粗糙度參數而構建的多源遙感數據源神經網絡模型具有明顯的優勢。