基于多源遙感協(xié)同反演的區(qū)域性土壤鹽漬化監(jiān)測(cè)

馮雪力 劉全明

(1.內(nèi)蒙古建筑職業(yè)技術(shù)學(xué)院，呼和浩特 010070； 2.內(nèi)蒙古農(nóng)業(yè)大學(xué)土木建筑工程學(xué)院，呼和浩特 010018)

0 引言

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

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

在前期研究基礎(chǔ)上，本文優(yōu)選高光譜土壤鹽分特征波段，融合微波后向散射特性，構(gòu)建基于多源遙感數(shù)據(jù)的土壤鹽分預(yù)測(cè)模型，為利用多源遙感數(shù)據(jù)協(xié)同反演鹽漬土?xí)r空分布提供科學(xué)基礎(chǔ)。為此，在高光譜鹽分特征波段優(yōu)選方法對(duì)比的基礎(chǔ)上，創(chuàng)建基于多源遙感數(shù)據(jù)的土壤鹽分經(jīng)驗(yàn)?zāi)Ｐ秃腿斯ぶ悄芊蔷€性預(yù)測(cè)模型，從而將雷達(dá)后向散射系數(shù)、光譜反射率及其變換值等參數(shù)轉(zhuǎn)換為土壤鹽分，以提高鹽漬化監(jiān)測(cè)的精度與廣適性。

1 研究區(qū)概況

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

2 材料與方法

2.1 數(shù)據(jù)獲取與處理

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

圖1 試驗(yàn)區(qū)雷達(dá)影像圖Fig.1 Radar image of experimental area

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

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

野外采樣點(diǎn)土壤反射率光譜數(shù)據(jù)的采集使用美國(guó)Analytical Spectral Device公司生產(chǎn)的Field Spec4便攜式光譜儀，可測(cè)量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進(jìn)行了3種不同形式的變換：對(duì)數(shù)(lgR)、一階導(dǎo)數(shù)(R′)、二階導(dǎo)數(shù)(R″)，并將光譜反射率及上述3種變換光譜作為本研究的輸入光譜數(shù)據(jù)，后期分別將上述4種輸入值與對(duì)應(yīng)采樣點(diǎn)鹽分值進(jìn)行相關(guān)性分析。

2.2 主成分回歸

主成分回歸是考察多個(gè)變量間相關(guān)性的一種多元統(tǒng)計(jì)方法，研究如何通過(guò)少數(shù)幾個(gè)主成分來(lái)揭示多個(gè)變量間的內(nèi)部結(jié)構(gòu)，即從原始變量中導(dǎo)出少數(shù)幾個(gè)主成分，使它們盡可能多地保留原始變量的信息，且彼此間互不相關(guān)。通常數(shù)學(xué)上的處理就是將原來(lái)多個(gè)指標(biāo)作線性組合，作為新的綜合指標(biāo)。

2.3 逐步回歸分析

逐步回歸分析是將變量逐個(gè)引入模型，每引入一個(gè)解釋變量后都要進(jìn)行F檢驗(yàn)，并對(duì)已經(jīng)選入的解釋變量逐個(gè)進(jìn)行t檢驗(yàn)，當(dāng)原來(lái)引入的解釋變量由于后面解釋變量的引入變得不再顯著時(shí)，則將其刪除。以確保每次引入新的變量之前回歸方程中只包含顯著性變量。這是一個(gè)反復(fù)的過(guò)程，直到既沒(méi)有顯著的解釋變量選入回歸方程，也沒(méi)有不顯著的解釋變量從回歸方程中剔除為止，以保證最后所得到的解釋變量集是最優(yōu)的。

2.4 偏最小二乘回歸

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

2.5 BP人工神經(jīng)網(wǎng)絡(luò)

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

3 特征波段選取建模與土壤鹽分預(yù)測(cè)

3.1 特征波段選取建模

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

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

利用PCR對(duì)光譜的一、二階導(dǎo)數(shù)中的特征波段進(jìn)行主成分分析，結(jié)果如圖3所示?？梢园l(fā)現(xiàn)，無(wú)論是一階還是二階導(dǎo)數(shù)變換的特征波段，都可以從中選取3特征波段來(lái)涵括，累積貢獻(xiàn)率分別為99.31%和99.87%；其中，第一主成分一階導(dǎo)數(shù)變換特征波段1 088～1 092 nm與二階導(dǎo)數(shù)變換特征波段618～622 nm的特征值和貢獻(xiàn)因子分別為4.07×10-7、1.51×10-7和72.19、74.08，皆具有較高代表性。

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

利用MSR對(duì)光譜的一、二階導(dǎo)數(shù)中的特征波段進(jìn)行逐步回歸分析，結(jié)果如圖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，研究表明，光譜一階導(dǎo)數(shù)的1 806～1 810 nm、2 138～2 142 nm 2個(gè)特征波段和二階導(dǎo)數(shù)的618～622 nm、1 802～1 806 nm、2 344～2 348 nm 3個(gè)特征波段更具有代表性。

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

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

圖5 反射率導(dǎo)數(shù)變換與土壤全鹽量的相關(guān)性及PLSR權(quán)重分析Fig.5 Analysis on correlation and PLSR weight between reflectance derivative transformation and soil salinity

PCR和MSR的擬合結(jié)果對(duì)比如表2所示。不難發(fā)現(xiàn)，二階導(dǎo)數(shù)變換的擬合程度高于一階導(dǎo)數(shù)變換，均方根誤差也都有所下降，說(shuō)明二階導(dǎo)數(shù)變換更適合于土壤全鹽量的擬合；綜合比較使用光譜二階導(dǎo)數(shù)變換特征波段的MSR經(jīng)驗(yàn)?zāi)Ｐ湍M全鹽量效果最優(yōu)，其決定系數(shù)(R2)和均方根誤差(RMSE)分別為0.343和9.414 g/kg。

表2 多模型擬合結(jié)果對(duì)比Tab.2 Results contrast of multi-model

圖6 土壤全鹽量實(shí)測(cè)值與模擬值Fig.6 Measured and simulated values of soil salinity

使用MSR選取的3個(gè)特征中心波段反射率、6個(gè)SAR后向散射系數(shù)及其組合值和地表組合粗糙度，建立土壤全鹽量與光譜二階導(dǎo)數(shù)的經(jīng)驗(yàn)?zāi)Ｐ?/p>

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

同時(shí)建立神經(jīng)網(wǎng)絡(luò)模型，由圖6可見(jiàn)，神經(jīng)網(wǎng)絡(luò)模型R2為0.890 8，遠(yuǎn)優(yōu)于MSR模型，體現(xiàn)了神經(jīng)網(wǎng)絡(luò)模型的優(yōu)越性。

3.2 土壤鹽分預(yù)測(cè)

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

表3 土壤鹽漬化分類(lèi)Tab.3 Classification of soil salinity

注：Ⅰ～Ⅴ級(jí)分別代表非鹽漬土、弱鹽漬土、中鹽漬土、強(qiáng)鹽漬土和鹽堿土。

圖7 研究區(qū)土壤全鹽量預(yù)測(cè)圖Fig.7 Prediction image of soil salinity in study area

表4 土壤鹽漬化預(yù)測(cè)分類(lèi)占比Tab.4 Proportion of all levels of soil salinity prediction classification

4 結(jié)論

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

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