基于雷達(dá)遙感的不同深度土壤含鹽量反演模型

張智韜 杜瑜燕 勞聰聰 楊 寧 周永財(cái) 楊亞龍

(1.西北農(nóng)林科技大學(xué)水利與建筑工程學(xué)院， 陜西楊凌 712100；2.西北農(nóng)林科技大學(xué)旱區(qū)農(nóng)業(yè)水土工程教育部重點(diǎn)實(shí)驗(yàn)室， 陜西楊凌 712100)

0 引言

土壤鹽漬化是灌溉農(nóng)業(yè)生產(chǎn)中的主要環(huán)境問(wèn)題，不僅影響灌區(qū)土地資源的經(jīng)濟(jì)效益，而且制約當(dāng)?shù)剞r(nóng)業(yè)的發(fā)展[1-2]，因此快速獲取土壤含鹽量信息對(duì)指導(dǎo)農(nóng)業(yè)生產(chǎn)具有重要意義。遙感技術(shù)具有快速、無(wú)創(chuàng)、宏觀的優(yōu)勢(shì)，為監(jiān)測(cè)土壤鹽漬化提供了技術(shù)手段[3-4]。

20世紀(jì)70年代，研究者開(kāi)始利用光學(xué)遙感進(jìn)行土壤鹽漬化的遙感反演[5-7]。EL HARTI等[8]用TM與OLI融合后，在鹽分指數(shù)(SI)的基礎(chǔ)上加入藍(lán)波段，構(gòu)建OLI-SR指數(shù)，提高了塔德拉灌區(qū)的土壤含鹽量反演精度。DOUAOUI等[9]利用Spot2的XS數(shù)據(jù)對(duì)阿爾及利亞下謝利夫盆地進(jìn)行了鹽度制圖，實(shí)施了針對(duì)謝利夫盆地地區(qū)土壤鹽分的時(shí)空監(jiān)測(cè)。陳紅艷等[10]利用Landsat8 OLI多光譜影像，通過(guò)加入短波紅外波段對(duì)植被指數(shù)進(jìn)行改進(jìn)，構(gòu)建了土壤含鹽量的支持向量機(jī)模型，獲得了較好的土壤鹽分空間分布的反演效果。以上研究都取得了較好的監(jiān)測(cè)效果，但主要集中于土壤表層的監(jiān)測(cè)，在對(duì)土壤含鹽量的反演中存在一定的局限性。

雷達(dá)遙感具有全天時(shí)、全天候和一定的穿透能力[11]，彌補(bǔ)了光學(xué)遙感的不足，利用雷達(dá)影像進(jìn)行土壤含鹽量監(jiān)測(cè)具有一定的優(yōu)勢(shì)[12]。土壤介電常數(shù)是微波輻射理論最重要、最基本的參數(shù)。土壤鹽分的積累會(huì)導(dǎo)致土壤介電常數(shù)虛部發(fā)生變化，進(jìn)而引起雷達(dá)后向散射系數(shù)的變化。已有研究表明，土壤介電常數(shù)是影響雷達(dá)影像后向散射系數(shù)的主要因素[13-14]。熊文成等[15]利用IEM后向散射理論分析發(fā)現(xiàn)，干旱期和多雨期后向散射系數(shù)差與介電常數(shù)虛部差值之間具有良好的線性關(guān)系，認(rèn)為雷達(dá)后向散射系數(shù)與土壤含鹽量具有良好的相關(guān)性，這為反演干旱區(qū)土壤含鹽量提供了一定的依據(jù)。何祺勝[16]探討了利用Radarsat雷達(dá)數(shù)據(jù)提取干旱區(qū)土壤鹽漬化信息的方法。依力亞斯江·努爾麥麥提等[17-19]以全極化合成孔徑雷達(dá)數(shù)據(jù)與地面實(shí)測(cè)0～10 cm深度范圍內(nèi)土壤含鹽量為數(shù)據(jù)源，進(jìn)行了渭干河-庫(kù)車(chē)河三角洲綠洲土壤含鹽量信息的提取分類(lèi)研究。郭曉靜[20]對(duì)比分析了不同極化組合的Radarsat-2雷達(dá)后向散射系數(shù)與河套灌區(qū)表層0～10 cm深度土壤含鹽量的關(guān)系，得出了提取鹽漬化土壤信息的最佳極化組合方式。劉全明等[21-22]利用Radarsat-2四極化雷達(dá)后向散射系數(shù)建立了內(nèi)蒙古河套灌區(qū)表層0～10 cm深度土壤含鹽量反演的人工智能模型。NURMEMET等[23]以新疆克里亞河地區(qū)為試驗(yàn)區(qū)，以表層0～20 cm為采樣深度，研究PALSAR合成孔徑雷達(dá)數(shù)據(jù)對(duì)土壤鹽漬化的監(jiān)測(cè)效果，認(rèn)為支持向量機(jī)分類(lèi)器是基于像素的分類(lèi)方法的最佳模型。以上研究雖然都建立了較好的監(jiān)測(cè)模型，但未分析雷達(dá)遙感對(duì)不同深度土壤含鹽量的監(jiān)測(cè)效果。

鑒于此，本文以河套灌區(qū)沙壕渠灌域裸土期不同深度的鹽漬化土壤為研究對(duì)象，分析Sentinel-1雷達(dá)衛(wèi)星后向散射系數(shù)和多種組合指數(shù)與不同深度土壤含鹽量之間的關(guān)系，構(gòu)建基于不同深度的土壤含鹽量反演模型，為雷達(dá)遙感監(jiān)測(cè)不同深度土壤含鹽量提供理論依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

試驗(yàn)在內(nèi)蒙古河套灌區(qū)沙壕渠灌域進(jìn)行(圖1)。沙壕渠總控制灌溉面積占整個(gè)解放閘灌區(qū)面積的3.7%，南北長(zhǎng)約15 km，東西約4 km，灌區(qū)地勢(shì)南高北低，平均海拔在1 034～1 037 m之間。沙壕渠處于中溫帶，受東南暖濕季風(fēng)影響較弱，西北干冷季風(fēng)影響較強(qiáng)，常年寒冷干燥，平均氣溫在7℃左右，年平均降水量140 mm，年平均蒸發(fā)量2 000 mm[24-25]。研究區(qū)主要以引黃河水漫灌為主，還輔以引潛水灌溉，長(zhǎng)期不恰當(dāng)?shù)墓喔确绞綄?dǎo)致灌區(qū)出現(xiàn)次生鹽漬化，研究區(qū)南部鹽漬化程度相對(duì)較輕，主要種植小麥、玉米；北部鹽漬化程度較重，主要種植耐鹽的向日葵。

圖1 研究區(qū)及采樣點(diǎn)分布Fig.1 Distribution of study area and sampling points

1.2 土壤采樣點(diǎn)設(shè)置與樣本處理

2019年4月12—21日，開(kāi)展河套灌區(qū)沙壕渠灌域的野外試驗(yàn)。選取沙壕渠灌域的大片裸露地區(qū)，采用五點(diǎn)法取0～10 cm深度的土層土樣與10～20 cm深度的土層土樣，共設(shè)置100個(gè)采樣點(diǎn)，土樣采集并編號(hào)后帶回實(shí)驗(yàn)室。每個(gè)采樣點(diǎn)通過(guò)手持GPS記錄采樣點(diǎn)位置信息及周?chē)h(huán)境信息，取樣期間天氣干燥，無(wú)降水。

將野外收集的土樣經(jīng)自然風(fēng)干后，配置土水質(zhì)量比為1∶5的土壤溶液[24]，經(jīng)攪拌、靜置、沉淀、過(guò)濾后，采用電導(dǎo)率儀(DDS-307A型，上海佑科儀器公司)測(cè)定土壤溶液電導(dǎo)率，對(duì)每個(gè)采樣點(diǎn)的5個(gè)土樣電導(dǎo)率取平均值作為該樣本樣點(diǎn)處的電導(dǎo)率，并通過(guò)經(jīng)驗(yàn)公式計(jì)算土壤含鹽量(SSC)[26]。將100個(gè)土壤樣本含鹽量按照鹽堿土分類(lèi)原則分為4個(gè)等級(jí)：非鹽漬化(0～0.2%)、輕度鹽漬化(0.2%～0.5%)、重度鹽漬化(0.5%～1.0%)和鹽漬化(>1.0%)，土壤含鹽量的樣本分析結(jié)果如表1所示。其中土壤深度為0～10 cm時(shí)非鹽漬化土、輕度鹽漬化土、重度鹽漬化土和鹽漬化土占比分別為65%、24%、8%和3%，且變異系數(shù)為1.072；土壤深度為10～20 cm時(shí)非鹽漬化土、輕度鹽漬化土、重度鹽漬化土和鹽漬化土占比分別為64%、27%、6%和3%，變異系數(shù)為0.966。

1.3 數(shù)據(jù)獲取與預(yù)處理

Sentinel-1衛(wèi)星具有雙極化、短重訪周期的特點(diǎn)[27]，可提供河套灌區(qū)區(qū)域連續(xù)的圖像(白天、夜晚和各種天氣)。試驗(yàn)所用雷達(dá)衛(wèi)星影像是2019年4月12日13:46過(guò)境的Sentinel-1，成像時(shí)間與地面采樣時(shí)間對(duì)應(yīng)，SAR影像覆蓋范圍為25 km×25 km，其地面分辨率為5 m×20 m，工作模式為干涉寬幅模式(Interferometric wide swath， IW)，數(shù)據(jù)級(jí)別為L(zhǎng)evel-1，GRDH，包括垂直極化模式VV(Vertical vertical)和交叉極化模式VH(Vertical horizontal)。采用歐空局開(kāi)發(fā)的SNAP軟件對(duì)GRDH數(shù)據(jù)進(jìn)行預(yù)處理，主要包括輻射校正、輻射地形校準(zhǔn)、應(yīng)用散斑過(guò)濾器和幾何地形校正等。去噪過(guò)程即相干波抑制采用Refined Lee濾波器對(duì)原始圖像進(jìn)行濾波，消除了大部分斑點(diǎn)，其中地形校正采用航天飛機(jī)雷達(dá)地形測(cè)繪使命(Shuttle radar topography mission，SRTM)生成的數(shù)字高程模型。

1.4 指數(shù)建立與灰度關(guān)聯(lián)

1.4.1指數(shù)建立

雷達(dá)后向散射系數(shù)反映地物特征值，經(jīng)過(guò)SNAP軟件處理和校正，得到雷達(dá)后向散射系數(shù)影像，把野外采集的點(diǎn)位信息的經(jīng)緯度輸入到雷達(dá)影像中，提取對(duì)應(yīng)點(diǎn)的后向散射系數(shù)，本文提取出VV極化與VH極化的雷達(dá)后向散射系數(shù)。有研究指出[20]，對(duì)于單極化雷達(dá)影像數(shù)據(jù)，如提取的土壤信息量相對(duì)較少，研究的結(jié)果會(huì)受到一定的影響，因此可以將雷達(dá)影像的極化方式進(jìn)行組合，從而提高雷達(dá)影像后向散射系數(shù)與土壤含鹽量的相關(guān)性[28-29]。對(duì)兩組雷達(dá)后向散射系數(shù)進(jìn)行極化組合變換，以達(dá)到預(yù)測(cè)土壤含鹽量的效果[30-31]。構(gòu)建的指數(shù)如表2所示，表中VV表示垂直極化的雷達(dá)后向散射系數(shù)，VH表示水平極化的雷達(dá)后向散射系數(shù)。

表2 極化組合指數(shù)Tab.2 Polarization combination index

1.4.2灰度關(guān)聯(lián)分析

灰色系統(tǒng)理論是一種多因素統(tǒng)計(jì)方法[32]，其目的是通過(guò)一定的方法原理確定系統(tǒng)中各因素的主要關(guān)系，用灰度關(guān)聯(lián)度刻畫(huà)因素間關(guān)系的強(qiáng)弱和次序，找出影響最大的因素[33]。主要計(jì)算步驟如下：設(shè)參考序列為x0={x0(t),t=1,2,…,n}，比較序列為xi={xi(t),t=1,2,…,n}，則x0和xi之間的灰度關(guān)聯(lián)度(Gray correlation degree，GCD)計(jì)算公式為

(1)

其中

(2)

式中ρ——介于0～1之間的分辨系數(shù)，取0.5

由文獻(xiàn)[34]可知，灰度關(guān)聯(lián)分析中，灰度關(guān)聯(lián)度會(huì)受到樣本數(shù)量的影響，本文中兩種深度的土壤含鹽量樣本數(shù)量各為100個(gè)，利用兩組后向散射系數(shù)進(jìn)行組合變換后得到20個(gè)指數(shù)。故采用0～10 cm深度含鹽量數(shù)據(jù)與10～20 cm深度含鹽量數(shù)據(jù)與構(gòu)建的20個(gè)指數(shù)樣本計(jì)算灰度關(guān)聯(lián)度，采用灰度關(guān)聯(lián)分析中常用的“均值化”數(shù)據(jù)預(yù)處理方法，消減組合指數(shù)數(shù)據(jù)與土壤含鹽量數(shù)據(jù)之間因量綱不一致而造成的分析誤差，進(jìn)而篩選出與土壤鹽分相關(guān)關(guān)系較好的指數(shù)。

1.5 鹽分反演模型構(gòu)建

1.5.1建模集與驗(yàn)證集劃分

將0～10 cm和10～20 cm按土壤含鹽量梯度進(jìn)行樣本劃分，每隔2個(gè)樣本取出一個(gè)作為驗(yàn)證樣本集，剩余土壤樣本作為建模樣本集。共選取67個(gè)樣本作為建模數(shù)據(jù)源，33個(gè)樣本作為驗(yàn)證數(shù)據(jù)源。各深度的土壤總樣本、建模集樣本和驗(yàn)證集樣本分布如圖2所示。0～10 cm土壤建模集與驗(yàn)證集樣本分布趨勢(shì)與總樣本分布趨勢(shì)一致；10～20 cm深度時(shí)總樣本、建模集樣本和驗(yàn)證集樣本分布也一致，說(shuō)明建模集與驗(yàn)證集樣本劃分合理。

圖2 樣本分布統(tǒng)計(jì)Fig.2 Sample distribution statistics

1.5.2模型構(gòu)建方法

將實(shí)測(cè)的0～10 cm與10～20 cm土壤含鹽量樣本數(shù)據(jù)各自按含鹽量劃分，隨機(jī)選取2/3作為模型輸入數(shù)據(jù)源，剩余的1/3作為模型驗(yàn)證數(shù)據(jù)源。并分別采用PLSR、QR和SVM 3種回歸方法構(gòu)建基于不同深度的土壤含鹽量反演模型。PLSR建模與分析使用The UnscramblerX 10.4軟件完成，QR與SVM兩種機(jī)器學(xué)習(xí)算法均在R3.4.0軟件中完成。使用R3.5.1軟件中的e1071包進(jìn)行SVM分析，設(shè)定核函數(shù)類(lèi)型為RBF，gamma為0.07，cost為100。

研究方法如圖3所示，雷達(dá)影像數(shù)據(jù)經(jīng)過(guò)輻射校正等步驟處理后，結(jié)合地面實(shí)際采樣點(diǎn)的經(jīng)緯度信息，提取對(duì)應(yīng)采樣點(diǎn)的后向散射系數(shù)并進(jìn)行極化組合；利用灰度關(guān)聯(lián)篩選最優(yōu)組合及部分地面實(shí)際采樣點(diǎn)土壤含鹽量數(shù)據(jù)作為模型的輸入數(shù)據(jù)源，剩余土壤含鹽量數(shù)據(jù)作為驗(yàn)證，得到最優(yōu)的估算模型。

圖3 研究方法流程圖Fig.3 Flow chart of research method

1.5.3模型精度評(píng)價(jià)

模型的建模和驗(yàn)證精度評(píng)估采用相關(guān)系數(shù)r、決定系數(shù)R2、均方根誤差RMSE。R2越大、RMSE越小說(shuō)明模型效果越好。

2 結(jié)果與分析

2.1 后向散射系數(shù)與土壤含鹽量的相關(guān)分析

根據(jù)地面采樣點(diǎn)的經(jīng)緯度，提取經(jīng)過(guò)SNAP軟件預(yù)處理后的VV和VH兩種極化方式下的雷達(dá)后向散射系數(shù)。將提取的后向散射系數(shù)分別與0～10 cm和10～20 cm兩個(gè)深度的土壤含鹽量建立相關(guān)關(guān)系，雷達(dá)后向散射系數(shù)與兩個(gè)深度的土壤含鹽量的相關(guān)關(guān)系如表3所示。

表3 不同極化方式的后向散射系數(shù)與土壤含鹽量的相關(guān)系數(shù)Tab.3 Correlation coefficient between backscattering coefficient of different polarization methods and soil salinity

從表3可以看出，雷達(dá)影像中，VH與VV極化方式下的后向散射系數(shù)與土壤含鹽量之間的相關(guān)關(guān)系不是特別好。其中10～20 cm深度的土壤含鹽量與VV極化模式下的后向散射系數(shù)之間的相關(guān)性最好，相關(guān)系數(shù)為0.264，達(dá)到了0.01顯著水平；0～10 cm深度的土壤含鹽量與VV極化模式下的后向散射系數(shù)之間的相關(guān)性較好，相關(guān)系數(shù)為0.235，達(dá)到了0.05顯著水平；而VH極化模式下的后向散射系數(shù)與0～10 cm深度和10～20 cm深度的土壤含鹽量之間的相關(guān)關(guān)系均較差。

圖4 基于不同深度土壤含鹽量的PLSR模型Fig.4 PLSR model based on soil salt content at different depths

2.2 基于灰度分析的極化組合指數(shù)選擇

在采樣點(diǎn)不同深度的條件下，利用灰色系統(tǒng)對(duì)構(gòu)建的所有極化組合指數(shù)與土壤含鹽量進(jìn)行灰色關(guān)聯(lián)分析，各極化組合指數(shù)與不同深度的土壤含鹽量(SSC)之間的灰色關(guān)聯(lián)度見(jiàn)表4。

從表4可以看出，同一深度，不同的極化組合指數(shù)與土壤含鹽量SSC的關(guān)聯(lián)度不相同；不同深度，不同的極化組合指數(shù)與土壤含鹽量SSC的關(guān)聯(lián)度也不相同?；叶汝P(guān)聯(lián)分析是一種多因素統(tǒng)計(jì)分析的方法，可通過(guò)比較各因素間的關(guān)聯(lián)度而得出系統(tǒng)中因素間的主次關(guān)系，找出影響最大的因素。從表4挑選出關(guān)聯(lián)度排序前10的極化組合指數(shù)進(jìn)行建模分析。

表4 不同深度土壤含鹽量與各指數(shù)的灰色關(guān)聯(lián)度Tab.4 Gray correlation statistics of soil salinity at different depths and each index

2.3 土壤含鹽量模型構(gòu)建

2.3.1基于PLSR的土壤含鹽量模型

綜合比較在0～10 cm和10～20 cm深度時(shí)的實(shí)測(cè)土壤樣本與預(yù)測(cè)樣本的擬合曲線，0～10 cm深度時(shí)的PLSR擬合曲線更接近于y=x函數(shù)，所以0～10 cm深度時(shí)的PLSR模型擬合精度高于10～20 cm。

2.3.2基于QR的土壤含鹽量模型

對(duì)于不同深度選定的極化組合指數(shù)建立分位數(shù)回歸模型，通過(guò)尋優(yōu)比較選取τ=0.6為最佳分位點(diǎn)。以不同深度下的雷達(dá)后向散射系數(shù)及多個(gè)極化組合指數(shù)為自變量，以土壤含鹽量為因變量，運(yùn)用QR模型中的0.6分位點(diǎn)進(jìn)行不同深度的土壤含鹽量估算，結(jié)果如圖5所示。

圖5 基于不同深度土壤含鹽量的QR模型Fig.5 QR model based on soil salt content at different depths

綜合比較在0～10 cm和10～20 cm深度時(shí)的實(shí)測(cè)土壤樣本與預(yù)測(cè)樣本的擬合曲線，0～10 cm時(shí)QR模型的擬合曲線更接近y=x函數(shù)分布，0～10 cm深度時(shí)的QR模型擬合精度高于10～20 cm。

2.3.3基于SVM的土壤含鹽量模型

以不同深度下的雷達(dá)后向散射系數(shù)和多個(gè)極化組合指數(shù)為自變量，以土壤含鹽量為因變量，進(jìn)行支持向量機(jī)建模分析。為了克服“離散值多”和“過(guò)學(xué)習(xí)”等問(wèn)題，尋優(yōu)選取參數(shù)gamma為0.007。SVM建模和驗(yàn)證結(jié)果如圖6所示。

圖6 基于不同深度土壤含鹽量的SVM模型Fig.6 SVM model based on soil salt content at different depths

綜合比較0～10 cm和10～20 cm深度時(shí)的實(shí)測(cè)土壤樣本與預(yù)測(cè)樣本的擬合曲線證明，0～10 cm深度時(shí)SVM模型的擬合曲線更接近y=x函數(shù)分布，0～10 cm深度時(shí)的SVM模型擬合精度高于10～20 cm。

2.4 模型綜合評(píng)價(jià)

基于不同采樣深度構(gòu)建PLSR模型、QR模型和SVM模型估算土壤含鹽量，對(duì)土壤含鹽量預(yù)測(cè)值與實(shí)測(cè)值進(jìn)行殘差分析，結(jié)果如圖7所示。

圖7 基于不同深度的土壤含鹽量估算模型的殘差分布Fig.7 Residual distribution of soil salt content estimation model based on different depths

由圖7可看出，在所有預(yù)測(cè)土壤含鹽量的模型中，0～10 cm深度條件下的PLSR模型、QR模型和10～20 cm深度條件下的PLSR模型殘差分布基本一致，殘差范圍跨度大，其中0～10 cm深度條件下的QR模型殘差范圍最廣，說(shuō)明預(yù)測(cè)效果最差。0～10 cm深度條件下的SVM模型殘差最集中并遵循標(biāo)準(zhǔn)正態(tài)分布，范圍最窄，預(yù)測(cè)效果最好。這與模型預(yù)測(cè)值的擬合曲線結(jié)論一致。

3 結(jié)論

(2)采用3種模型對(duì)不同深度土壤含鹽量進(jìn)行反演時(shí)，0～10 cm深度反演精度均高于10～20 cm深度的反演精度，其中0～10 cm深度時(shí)支持向量機(jī)模型是最優(yōu)模型。

(3)0～10 cm深度條件下的支持向量機(jī)模型適用于沙壕渠裸土期的土壤含鹽量反演。但是地表粗糙度對(duì)雷達(dá)后向散射系數(shù)有一定的影響，會(huì)影響雷達(dá)遙感反演土壤含鹽量的預(yù)測(cè)精度，這也是今后應(yīng)研究并解決的問(wèn)題。