張智韜 魏廣飛 姚志華 譚丞軒 王新濤 韓 佳
(1.西北農林科技大學水利與建筑工程學院, 陜西楊凌 712100;2.西北農林科技大學旱區農業水土工程教育部重點實驗室, 陜西楊凌 712100)
土壤鹽漬化已成為干旱-半干旱地區最嚴重的土壤問題之一,嚴重制約了區域經濟健康發展[1]。快速、準確地獲取鹽漬土時空分布信息成為合理防治鹽漬化的必要前提。傳統的土壤含鹽量測定方法費時、費力、成本高,且無法全面獲取數據。近年來,隨著光譜分析技術的發展,遙感技術在農業中的應用越來越廣泛,尤其在土壤鹽漬化監測方面[2-7]。已有研究結果顯示,在干旱和半干旱地區,植被的生長狀況可以間接反映土壤鹽漬化程度[8-10]。ALLBED等[11]研究表明,土壤含鹽量超過植被生長閾值時,致使葉片構造發生變化,在可見光和近紅外波段,植被冠層光譜反射率對土壤含鹽量響應敏感。ZHANG等[12]利用MODIS衛星影像數據構建增強性植被指數和EVI,反演了黃河三角洲植被區土壤鹽分。陳紅艷等[13]基于Landsat8 OLI多光譜影像引入第七波段對植被指數進行改進,構建土壤含鹽量的支持向量機模型,獲得了良好的土壤鹽分空間分布反演效果。以上研究大多以植被指數和鹽分指數作為指示因子建立與實測土壤含鹽量的關系[14-17]。王飛等[18]通過引入波段組、植被指數變量組、優選變量組等6個變量組,進行多變量分析,通過算法篩選植被指數、土壤鹽度指數等鹽度敏感變量,得出最優變量組合,并利用多個機器學習算法預測綠洲土壤鹽分,研究發現,在涉及的變量個數較多、且與各土壤含鹽量相關性較低時,機器學習算法的挖掘能力得到了充分的體現。BP神經網絡、支持向量機和隨機森林是3種常用的機器學習算法, BP神經網路算法具有非線性擬合能力和自學習能力,支持向量機算法可以避開從歸納到演繹的傳統過程,隨機森林算法具有可處理非線性數據、實現簡單、訓練速度快、抗擬合能力強的特點。王明寬等[16]研究表明,遙感影像的反射率與土壤含鹽量并不是單純的線性關系,BP神經網絡模型能很好地模擬土壤含鹽量與光譜數據的關系。厲彥玲等[19]利用遙感數據建立統計分析模型(多元線性回歸、偏最小二乘回歸),用機器學習模型(BP神經網絡、支持向量機模型和隨機森林)反演鹽分,結果表明,3種機器學習模型反演效果均取得了較高的精度,且反演效果遠優于統計分析模型。
國內外學者對區域鹽分定量監測已取得較多的成果,但這些成果多偏向于衛星遙感和高光譜遙感,IVUSHKIN等[20]通過在無人機上裝載熱成像、高光譜、激光雷達3種傳感器發現,無人機遙感在農田鹽分監測方面有較大潛力。但基于無人機多光譜遙感監測植被覆蓋下的土壤含鹽量研究較少,同時缺乏無人機多光譜遙感對植被覆蓋下的土壤含鹽量反演模型。
本文以河套灌區的沙壕渠灌域為研究區域,基于無人機搭載多光譜相機獲取的遙感數據和同步測定的土壤鹽分數據,引入敏感波段組、光譜指數組和全變量組作為模型輸入變量,利用多元線性回歸、BP神經網絡、支持向量機和隨機森林算法,構建土壤鹽分的定量反演模型,評定不同模型輸入變量和不同模型方法的土壤含鹽量反演模型的精度,優選土壤含鹽量最佳反演模型,以期為干旱及半干旱地區無人機遙感土壤鹽分定量反演提供參考。

圖1 研究區位置示意圖Fig.1 Sketch of study area
沙壕渠灌域是內蒙古河套灌區西北部解放閘灌域內部的一個獨立單元,位于40°52′~41°00′N,107°05′~107°10′E,灌域形狀近似為一個狹長的倒三角形(圖1),南窄北寬,南北約為15 km,東西平均寬度約4 km,地面較為平整,地勢走向為南高北低,海拔為1 034~1 037 m,總控制面積為52.4 km2,當地屬干旱-半干旱氣候,冬季嚴寒少雪,夏季高溫少雨,年均降雨量140 mm,年均蒸發量2 000 mm,年平均氣溫7℃。灌區土壤類型以粉壤土、砂壤土和壤土為主。研究區無霜期為120~150 d,主要農作物包括小麥、玉米、葵花、白菜等。沙壕渠灌域南部鹽漬化程度相對較輕,作物種植以小麥、玉米為主,而北部鹽漬化程度較重,以種植耐鹽的向日葵為主。
2018年5月,對沙壕渠灌域的耕地進行實地調研,根據土壤的鹽漬化程度,在灌域確定4塊不同鹽漬化梯度的試驗區域,并依次進行編號,每塊研究區域為16 hm2,每塊區域均勻布設15個作物(玉米、西葫蘆和葵花)覆蓋下的地面數據采集點,共計60個采樣點,如圖2所示。

圖2 采樣點分布示意圖Fig.2 Sampling point distribution maps
1.3.1多光譜遙感影像獲取
所用遙感平臺為大疆創新公司生產的經緯M600型六旋翼無人機,其搭載的傳感器為美國Tetracam公司生產的Micro-MCA多光譜相機(簡稱MCA),具有質量輕、體積小及遠程觸發的特點,非常適合在中小型無人機上進行搭載及拍攝,MCA相機包括490 nm(藍光)、550 nm(綠光)、680 nm(紅光)、720 nm(紅邊)、800 nm(近紅外)、900 nm(近紅外)6個波長的光譜采集通道。2018年8月12—16日于河套灌區沙壕渠灌域進行試驗,天氣情況分別為:晴轉多云、晴、晴、晴、晴轉多云,影像獲取時間選為11:00—15:00之間(均在晴朗時拍攝),無人機飛行模式按照提前規劃的航線飛行,拍照時相機鏡頭與地面呈90°,拍照模式為等時間間隔,主航線上和主航線間圖像重疊率均設置為80%以上,飛行高度經多次試飛后選定為120 m(此時影像地面分辨率為6.5 cm,影像刈幅寬度為66.5 m),采集光譜前用標準白板進行標定,4塊區域分別得到341、320、390、344幅遙感圖像。所用無人機與多光譜相機如圖3所示,無人機飛行路線(以4號地為例)及現場作業如圖4所示。

圖3 M600型無人機及Micro-MCA多光譜相機Fig.3 M600 unmanned aerial vehicle and micro-MCA multispectral camera

圖4 無人機飛行路線及現場作業圖Fig.4 UAV flight path and application scence of UAV over cropland
1.3.2多光譜遙感影像處理
利用多光譜相機配套的pixel Wrench2軟件對獲取的多光譜正射影像進行提取、配準與合成,然后將獲得的TIF格式的遙感圖像和每幅圖像的GPS數據導入pix4Dmapper軟件中完成4塊研究區域的圖像拼接,把拼接完成的圖像輸入ENVI Classic中作進一步處理,為排除單一點反射率產生的隨機誤差,由地面數據采集點的GPS裁剪出以采集點為中心的200像素×200像素圖像,采用監督分類的最大似然法剔除土壤背景后,將此區域6個波長平均灰度除以標準白板的灰度作為該采集點對應6個波長的反射率。
每塊研究區域均勻布設15個地面數據采集點,為保證遙感影像與地面數據的時間一致性,在獲取多光譜遙感影像的同時,用對角線5點取樣法在作物根系附近采集深度為0~20 cm、20~40 cm的土樣,并用手持式GPS儀記錄每個采樣點的位置,取出的土樣放置在標記的鋁盒中,并帶回實驗室干燥,采用土水比1∶5法配置土壤溶液,經過濾后靜置8 h,采用電導率儀(雷磁DDS-307A型,上海佑科儀器分公司)測定土壤溶液電導率(EC1∶5,dS/cm),并通過經驗公式計算土壤含鹽量(SSC,%)[21]:(SSC=0.288 2EC1∶5+0.018 3)。該時期作物的根系主要活動層深度為0~40 cm[22],故將此深度的含鹽量作為該試驗的土壤含鹽量數據(此時,0~40 cm的含鹽量為0~20 cm、20~40 cm含鹽量的平均值)。
采用遙感影像提取的各類光譜指數構建特征空間進行土壤含鹽量反演和檢測是目前土壤鹽分遙感監測的主要方法之一。本研究選擇10個光譜指數,分別為歸一化植被指數(Normalized difference vegetation index,NDVI)、差值植被指數(Difference vegetation index,DVI)、增強性植被指數(Enhanced vegetation index,EVI)、增強性植被指數2(Enhanced vegetation index 2,EVI2)、三角形植被指數(Triangular vegetation index,TVI)、土壤調節植被指數(Soil-adjusted vegetation index,SRVI)、歸一化差異綠度植被指數(Normalized difference greenness vegetation index,NDGI)、歸一化鹽分指數(Normalized difference soil index,NDSI)、鹽分指數(Salinity index,SI-T)、簡單比值指數(Simple ratio index,SR),其計算公式如表1所示。

表1 光譜指數及相關計算公式Tab.1 Spectral indexes and related formula
注:G、R、NIR、NIR2分別為550、680、800、900 nm波長處的光譜反射率。L為蓋度背景調節因子,取0.5。
將60個研究樣本按照含鹽量由小到大排序分組,根據建模集與驗證集為2∶1的比例進行等間隔取樣。選取40個樣本用于建模,其余20個樣本用于驗證。
本研究將模型輸入變量分為3組:敏感波段變量組、光譜指數變量組和全變量組。采用多元線性回歸模型、支持向量機、反神經傳播網絡和隨機森林等4種回歸方法分別建立基于3種模型輸入變量的土壤鹽分反演模型,共計12個鹽分反演模型。本研究的MLR建模與分析使用IBM SPSS 23軟件完成,3種機器學習算法模型均在R3.4.0軟件中完成。這4種回歸方法如下:
(1)多元線性回歸
MLR是用于兩個或兩個以上自變量時的一種多元回歸分析,由多個自變量的最優組合共同來預測因變量,比只用一個自變量進行預測或估計更有效、更符合實際[16],是遙感鹽分反演最常用的回歸方法之一[25-26]。
(2)支持向量機
SVM是根據統計學理論,以結構風險最小化為原則為理論基礎從線性可分擴展到線性不可分的一種新的機器學習方法,在圖像識別與分類上應用較廣,近年來在回歸問題上也得到了一些應用[27-28]。本文設定核函數類型為“poly”,采用訓練集交叉驗證和網絡搜素法(Grid search)進行參數尋優,根據方差最小原則確定懲罰系數為C=10 000,γ=0.01。
(3)反神經傳播網絡
BPNN是一種按誤差逆傳播算法訓練的多層前饋神經網絡,也已應用到了鹽分反演問題中[29-31]。本文采用3層網絡拓撲結構構建BPNN模型,按照訓練結果誤差相對較小的原則,權值的衰減參數和隱層節點數通過多次反復試驗確定,最終確定衰減參數為0.001、隱層節點數為4,其他參數均使用默認值。
(4)隨機森林
RF算法是將Bagging算法與決策樹算法進行結合所得到的集成學習算法,近年來多位學者將其應用到遙感技術中[31-32],此算法主要參數為決策樹的個數,采用訓練集交叉驗證和網路搜索法(Grid search)進行決策樹個數尋優,經過多次訓練確定決策樹個數為4。
利用20個驗證樣本對各類方法構建的12個鹽分反演模型進行檢驗分析,采用驗證集將模型預測值和實測值進行擬合,對模型建模和驗證的精度采用決定系數(Coefficient of determination,R2)和均方根誤差(Root mean square error, RMSE)兩個指標來評價。R2越接近1,RMSE越小,說明模型精度效果越好。
將60個土樣測得的土壤含鹽量分為3個等級[21],即:非鹽漬化(<0.2%)、輕度鹽漬化(0.2%~0.5%)和重度鹽漬化(0.5%~1.0%),土壤含鹽量的總樣本、建模集、驗證集的分析結果如表2所示。河套灌區總樣品統計結果顯示,非鹽漬土、輕度鹽漬土、重度鹽漬土占比分別為16.7%、63.3%、20%,且變異系數為0.49,表明表層土壤含鹽量呈現中度變異性。

表2 土壤鹽分的描述性統計分析Tab.2 Descriptive statistical analysis on soil salinity
在干旱和半干旱地區,生長在鹽漬化環境中的植物,其根系接觸到土壤溶液中各種可溶性鹽離子,當陰陽離子超過一定含量時,就會對植物生長產生各種危害,即土壤鹽漬化會在一定程度上影響地表作物生長長勢,進而影響冠層的光譜特征[33]。分別對多光譜六波段反射率、不同光譜指數分別與土壤含鹽量進行相關性分析,根據皮爾遜(Person)相關系數界值表,n=60時,當|r|>0.250時,表示在0.05水平上顯著,|r|>0.325時,表示在0.01水平上顯著,|r|>0.408時,表示在0.001水平上顯著,二者與土壤鹽分的Pearson相關系數如表3、4所示。

表3 多光譜六波段反射率與SSC的相關系數(n=60)Tab.3 Correlation coefficient of reflectivity in six bands of multispectral and SSC(n=60)
注:*表示顯著性檢驗p<0.05,** 表示p<0.01,*** 表示p<0.001,下同。
從表3可以看出,除了紅邊波長(720 nm)處相關性僅為0.094外,其他5個波長與SSC均表現出一定的相關性。在紅光波長(680 nm)、近紅外波長(800 nm)和近紅外波長(900 nm)處顯著性檢驗p<0.001,且相關系數均達到了0.5以上。六波長與SSC相關性由大到小依次為B6、B5、B3、B2、B1、B4。
從表4可看出,選取的10種光譜指數與SSC均呈現出良好的相關性,相關性由大到小依次為SI-T、TVI、EVI2、SRVI、NDVI、EVI、NDGI、NDSI、SR、DVI。其中,TVI、EVI2和SI-T相關系數均超過了0.6, NDVI和SRVI相關系數都接近0.6,說明SI-T、TVI、EVI2、SRVI、NDVI在一定程度上可表征土壤含鹽量。
將模型輸入變量分為3組。通過2.2節中多光譜六波長與SSC的相關性分析,選擇與土壤鹽分極顯著且Person相關系數大于0.5的B3、B5和B6波段,將這3個敏感波段作為第1組模型輸入變量,記為敏感波段組;通過2.2節中10種光譜指數與SSC的相關性分析,選擇SI-T、TVI、EVI2、SRVI和NDVI共5個相關性較高的光譜指數作為第2組模型輸入變量,記為光譜指數組;為了盡量減少因重要變量缺失而引起的模型偏差,將上述3個敏感波長和5個敏感光譜指數作為第3組模型輸入變量,共計8個變量,記為全變量組。

表4 光譜指數與SSC的相關系數(n=60)Tab.4 Correlation coefficient of spectral indexes and SSC(n=60)
2.3.1基于3個變量組的多元線性回歸土壤鹽分反演模型
將敏感波段組、光譜指數和全變量組分別作為模型輸入變量,建立基于多元線性回歸的鹽分反演模型,并對模型進行精度驗證,基于3種變量組的多元線性回歸模型結果見表5。


表5 基于不同變量組的多元線性回歸模型Tab.5 MLR model based on different groups of variables
注:敏感波段變量組中:X1、X2、X3分別表示B3、B5和B6波長對應的反射率;光譜指數組中:X1、X2、X3、X4、X5分別表示SI-T、TVI、EVI2、SRVI、NDVI;全變量組中:X1、X2、X3、X4、X5、X6、X7、X8分別表示B3、B5、B6、SI-T、TVI、EVI2、SRVI、NDVI。
綜合對比MLR模型的3個變量輸入組建模和驗證參數結果,發現基于光譜指數組建立的MLR模型較為穩定,是3個變量組中最佳的多元線性回歸鹽分反演模型。
2.3.2基于3個模型輸入變量組的機器學習鹽分反演模型
將敏感波段、光譜指數、全變量組作為模型輸入變量,分別建立SVM、BPNN、RF等3種機器學習回歸模型并進行模型精度驗證,模型結果如表6所示。

表6 基于不同變量組的機器學習模型Tab.6 Machine learning models based on different variable groups

采用敏感波段組、光譜指數組和全變量組作為模型輸入變量組,分別使用MLR、SVM、BPNN、RF共4種回歸方法,共構建12個土壤鹽分反演模型。將12個反演模型的土壤鹽分預測值與驗證集土壤鹽分實測值進行比較,結果如圖5所示。


圖5 實測值與預測值的散點圖Fig.5 Scatter plots between measured and predicted values

本研究取得了良好的鹽分反演效果,但是由于不同作物對鹽分的敏感響應程度不同,從而造成獲取的光譜特征也不同,所以下一步可以考慮對不同的作物進行分開鹽分反演。此外,本研究的鹽分反演模型局限于小灌域尺度,對大區域甚至是全國尺度的統一鹽分定量鹽分反演模型亟待研究。
(1)對比了基于3個模型輸入變量組構建的土壤含鹽量反演模型,光譜指數組取得了最優的反演精度。

