向友珍 李汪洋 臺 翔 安嘉琪 王 辛 陳俊英

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

0 引言

土壤鹽漬化是指在自然和人為作用下,鹽分在土壤表層不斷積累的現(xiàn)象,是造成土壤退化、農(nóng)業(yè)減產(chǎn)和生態(tài)環(huán)境惡化的重要因素[1]。河套灌區(qū)土壤鹽漬化問題嚴(yán)重,對鹽漬化土壤進(jìn)行治理和開發(fā)利用是農(nóng)業(yè)可持續(xù)發(fā)展與生態(tài)文明建設(shè)高質(zhì)量發(fā)展的有力保障[2]。實(shí)時、高效地監(jiān)測土壤鹽漬化狀況,是治理土壤鹽漬化的重要前提。

人工定點(diǎn)采樣法是土壤鹽漬化監(jiān)測的方式之一,該方法監(jiān)測周期長、工作量大,且監(jiān)測的精度與采樣點(diǎn)的數(shù)量、布設(shè)方式有很大關(guān)系,很難反映鹽漬化在空間上的實(shí)際分布。隨著科學(xué)技術(shù)的不斷發(fā)展,遙感被廣泛應(yīng)用于鹽漬化研究和監(jiān)測,為大面積時空分布動態(tài)化的鹽漬化監(jiān)測提供了新的途徑。然而廣泛可用的衛(wèi)星圖像無法提供高空間分辨率和數(shù)據(jù)采集時間的靈活性,在農(nóng)業(yè)應(yīng)用方面有很大的局限性[3]。

除衛(wèi)星遙感外,利用無人機(jī)平臺搭載的小型多光譜遙感設(shè)備對區(qū)域土壤鹽漬化監(jiān)測也是重要手段[4-5]。與衛(wèi)星圖像相比,基于無人機(jī)多光譜圖像的鹽漬化監(jiān)測具有空間分辨率高、光譜分辨率強(qiáng)、波段連續(xù)性強(qiáng)等優(yōu)點(diǎn),可以獲得多維、高精度的鹽漬化檢測信息,實(shí)現(xiàn)對土壤鹽漬化的動態(tài)監(jiān)測;同時,無人機(jī)傳感器攜帶多光譜相機(jī)平臺,具有實(shí)時、高分辨率、移動性和靈活性的優(yōu)勢,在高精度鹽漬化監(jiān)測方面具有優(yōu)勢[6]。劉旭輝等[7]利用無人機(jī)多光譜遙感影像建立的機(jī)器學(xué)習(xí)模型反演土壤含鹽量,得出不同季節(jié)的土壤含鹽量會有所不同。劉楠等[8]將局部區(qū)域的相關(guān)性分析和多維光譜變化特性采用二進(jìn)制編碼的形式來表征不同類型地物紋理特征,實(shí)現(xiàn)了區(qū)分和判別影像上不同紋理的目的,并證明該方法對提取地物紋理特性具有一定可行性。萬亮等[9]利用無人機(jī)多光譜融合植被指數(shù)和紋理特征提高了水稻含水量預(yù)測結(jié)果;周聰?shù)萚10]通過實(shí)踐證明,高分辨率遙感圖像的紋理特征可以作為估算植被生長參數(shù)的更有效指標(biāo)。陳鵬飛等[11]通過實(shí)踐證實(shí)了基于無人機(jī)多光譜圖像的剔除土壤背景和增加紋理特征處理可以提高棉花植株氮濃度的反演精度。但是基于無人機(jī)多光譜圖像剔除土壤背景和紋理特征是否可以提高土壤含鹽量精度方面研究較少。

本文以河套灌區(qū)沙壕渠灌域4塊不同土壤鹽漬化程度的典型鹽漬化試驗(yàn)地為試驗(yàn)研究區(qū),以灌區(qū)典型地物為研究對象,在灌區(qū)作物生長期內(nèi)采集試驗(yàn)區(qū)鹽漬土樣本并獲取土壤鹽漬化信息。同時利用無人機(jī)多光譜采集遙感數(shù)據(jù),對比研究基于鹽漬化土壤和耐鹽作物光譜的土壤含鹽量反演模型的性能,建立基于無人機(jī)多光譜遙感平臺的區(qū)域鹽漬化監(jiān)測體系,為區(qū)域鹽漬化研究提供可靠依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

本研究區(qū)域位于巴彥淖爾市沙壕渠灌域(40°52′～41°00′N,107°05′～107°10′E),隸屬于河套灌區(qū)解放閘灌域,如圖1所示。灌域因不合理的灌排方式,氣候、土質(zhì)和地貌等因素的綜合影響,土壤鹽漬化問題突出。區(qū)內(nèi)種植作物以向日葵、玉米等耐鹽糧油作物為主。沙壕渠灌域土壤類型為粉壤土、砂壤土和壤土。沙壕渠灌域面積約52.4 km2,南北跨度約15.1 km,東西橫跨4.2 km,為典型的溫帶大陸性氣候,多年平均氣溫3.7～7.6℃,多年平均降水量210～290 mm,多年平均年蒸發(fā)量2 100～3 080 mm,多年平均日照時長為3 000～3 200 h。

1.2 試驗(yàn)地布設(shè)

將試驗(yàn)地以土壤含鹽量為基礎(chǔ)從低到高依次分為1、2、3、4,每塊試驗(yàn)地面積為15.5 hm2左右。4塊試驗(yàn)地主要種植作物為向日葵和玉米。每塊試驗(yàn)地均勻布設(shè)30個土壤采樣點(diǎn)。采樣點(diǎn)布設(shè)見文獻(xiàn)[12]。

1.3 光譜圖像采集和處理

使用的無人機(jī)為深圳市大疆創(chuàng)新科技有限公司生產(chǎn)的M600型六旋翼無人機(jī),其最大上升速度5 m/s,最大下降速度3 m/s,最大飛行速度18 m/s,飛行承載質(zhì)量6 000 g,飛行高度2 500 m,單次飛行時間35～40 min。多光譜遙感相機(jī)采用美國Tetracam公司生產(chǎn)的6通道Micro-MCA多光譜相機(jī),包括藍(lán)光波段、綠光波段、紅光波段、紅邊波段、近紅外1波段、近紅外2波段共6個遙感波段,波長分別為490、550、680、720、800、900 nm。試驗(yàn)時間為2022年7月16—20日。每次試驗(yàn)均在11:00—14:00進(jìn)行,試驗(yàn)日晴朗無風(fēng),以確保充分的輻射強(qiáng)度,盡量減小植被陰影對光譜的影響。根據(jù)提前規(guī)劃好的航線,設(shè)置無人機(jī)飛行高度120 m,對應(yīng)多光譜相機(jī)分辨率為6.5 cm,相機(jī)拍攝速率為18～19幅/min,每次試驗(yàn)均設(shè)有白板進(jìn)行圖像標(biāo)定。

1.4 光譜指數(shù)計算

光譜指數(shù)是綜合考慮地物的各波段光譜特征,對不同波段反射率進(jìn)行數(shù)學(xué)變換組合,以增強(qiáng)地物特定的信息[13]。分別選擇8種植被指數(shù)和10種鹽分指數(shù),其計算公式如表1所示。

表1 光譜指數(shù)Tab.1 Spectral index

1.5 灰度共生矩陣法

灰度共生矩陣(GLCM)是一種基于統(tǒng)計數(shù)據(jù)的圖像紋理特征提取方法,由HARALICK等[26]提出。

GLCM的元素是不同組合出現(xiàn)的頻數(shù),為了計算圖像紋理的特征參數(shù),對GLCM進(jìn)行歸一化。已知i、j分別為(x,y)、(x+Δx,y+Δy)的像素值,(i,j)的個數(shù)為V(i,j),將歸一化GLCM記為P,則

(1)

(2)

式中Vi,j——矩陣元素N——矩陣元素個數(shù)

基于歸一化灰度共生矩陣,計算圖像紋理特征參數(shù),本研究采用8個紋理特征參數(shù),其計算公式和特性如表2所示。

表2 圖像紋理特征參數(shù)Tab.2 Texture feature parameters of images

1.6 Otsu算法

Otsu算法是一種計算簡單、自適應(yīng)強(qiáng)且已得到最廣泛使用的圖像閾值自動選取方法[29]。設(shè)某一灰度級對應(yīng)閾值為T,類間方差計算公式為

(3)

式中Wb——閾值T下背景占整幅圖像的比重

Wf——閾值T下前景占整幅圖像的比重

1.7 全子集篩選法

全子集篩選法是利用全子集回歸分析,對自變量不同的組合,用最小二乘法進(jìn)行建模分析,篩選最優(yōu)的變量組合。選擇最優(yōu)模型的評價標(biāo)準(zhǔn)為:①似然函數(shù)最大化。②模型未知參數(shù)最小化[30]。本文通過R編程語言進(jìn)行全子集篩選,利用決定系數(shù)R2和貝葉斯信息準(zhǔn)則(Bayesian information criterion,BIC)來評價篩選結(jié)果,對比分析R2越大、BIC越小的篩選結(jié)果,為最優(yōu)變量組合。

1.8 建模策略

建模策略分別為:未剔除土壤背景的光譜指數(shù)(策略1)、剔除土壤背景后的光譜指數(shù) (策略2)、未剔除土壤背景的光譜指數(shù)+圖像紋理特征(策略3)、 剔除土壤背景的光譜指數(shù)+圖像紋理特征(策略4)。

1.9 模型方法和評價指標(biāo)

極限學(xué)習(xí)機(jī)(Extreme learning machine, ELM)是一種基于最小二乘學(xué)習(xí)算法的隱層前饋網(wǎng)絡(luò)[31]。ELM收斂速度比傳統(tǒng)算法快,因?yàn)樗鼰o需迭代即可學(xué)習(xí),同時,隨機(jī)隱藏節(jié)點(diǎn)保證了全局逼近能力[32]。本文ELM模型采用R語言elmNNRcpp包構(gòu)建。支持向量機(jī)(Support vector machines,SVM)是用于監(jiān)督學(xué)習(xí)的強(qiáng)大計算工具[33]。本研究采用非線性支持向量機(jī)模型,利用R語言e1071包構(gòu)建。通過決定系數(shù)R2、均方根誤差(RMSE)、標(biāo)準(zhǔn)均方根誤差(NRMSE)指標(biāo)評價模型精度。R2越接近1,RMSE越接近0,說明模型效果越好。

本文利用R語言編程構(gòu)建SVM和ELM土壤含鹽量反演模型,其中實(shí)測土壤含鹽量為因變量,光譜變量為自變量。將深度0～20 cm平均土壤含鹽量按照比例2∶1隨機(jī)劃分建模集和驗(yàn)證集,通過調(diào)整參數(shù),獲得每個條件下最佳模型。建模流程如圖2所示。

圖2 建模流程圖Fig.2 Model building flowchart

2 結(jié)果與分析

2.1 土壤含鹽量

對試驗(yàn)地總計87個采樣點(diǎn)的土壤含鹽量進(jìn)行統(tǒng)計,將計算表層和深度10～20 cm實(shí)測土壤含鹽量的平均值作為深度10～20 cm的平均土壤含鹽量。將深度0～20 cm平均土壤含鹽量按照比例2∶1隨機(jī)劃分建模集和驗(yàn)證集,結(jié)果如表3所示。將各總集、建模集、驗(yàn)證集樣本點(diǎn)土壤含鹽量劃分為4個等級:非鹽土(D1,SSC(土壤含鹽量)小于等于0.2%)、輕度鹽漬化(D2,SSC為(0.2%、0.5%])、重度鹽漬化(D3,SSC為(0.5%、1%])和鹽土(D4,SSC大于1.0%)[34]。非鹽土、輕度鹽漬化、重度鹽漬化和鹽漬土占比分別為34.1%、55.7%、9.1%和1.1%。含鹽量變異系數(shù)均處于中等差異(變異系數(shù)CV反映樣點(diǎn)值的離散程度,CV<0.1為弱變異性;0.11為強(qiáng)變異性)[34]。

表3 土壤含鹽量特征統(tǒng)計分析Tab.3 Statistical analysis of soil salt content characteristics

對各鹽漬化等級下建模集、驗(yàn)證集和總集的含鹽量進(jìn)行統(tǒng)計分析,如圖3(圖中1.5IOR表示1.5倍的四分位距)所示。從圖3可以看出,D1、D2、D3等級的建模集、驗(yàn)證集和總集的含鹽量分布、值域和均值相近,確保建模集和驗(yàn)證集數(shù)據(jù)的代表性。

圖3 土壤含鹽量特征統(tǒng)計Fig.3 Statistical map of soil salinity characteristics

2.2 不同植被指數(shù)分類結(jié)果

基于近紅外和可見光波段構(gòu)建4個植被指數(shù),利用對土壤背景敏感的4個植被指數(shù)(NDVI、MSAVI、DVI、CRSI)進(jìn)行圖像分類,驗(yàn)證不同指數(shù)的圖像分類精度。為了對比4個植被指數(shù)圖像分類的精度,從87個采樣點(diǎn)中隨機(jī)選取40個樣本,進(jìn)行如圖4所示,4個植被指數(shù)的精度存在顯著差異,其中NDVI的圖像分類精度最高,40個樣本的總體精度均大于93.4%,均值為96.9%;Kappa系數(shù)在0.84～0.99之間,均值為0.92。MSAVI的分類精度略低于NDVI,總體精度和Kappa系數(shù)均值分別為95.8%和0.91;DVI和CRSI 2個植被指數(shù)的分類效果相對較差,總體精度均值小于91%,Kappa系數(shù)均值均小于0.88。基于以上分析,NDVI圖像分類結(jié)果表現(xiàn)最優(yōu)。

圖4 基于不同植被指數(shù)的圖像分類結(jié)果Fig.4 Evaluation based on image classification results of different vegetation indices

分類結(jié)果評價。分別對40個樣本的多光譜圖像進(jìn)行2.1節(jié)的操作,獲得分類結(jié)果,并在ENVI 5.3中對分類結(jié)果進(jìn)行評價,得到4個植被指數(shù)的總體精度和Kappa系數(shù)。

2.3 圖像紋理計算

灰度共生矩陣具有豐富的特征參數(shù),能從不同的角度對紋理進(jìn)行細(xì)致刻畫[35]。由鹽漬化土壤覆蓋植被冠層的光譜特征分析可知,植被冠層光譜中紅光波段與土壤含鹽量的相關(guān)性最高,對鹽漬化的變化有最顯著的響應(yīng),故選擇波段Band3計算多光譜圖像紋理。為了更好地體現(xiàn)植株間冠層圖像的紋理差異,紋理特征參數(shù)計算的窗口選擇尺度為11像素×11像素,尺寸為71.5 cm;每個特征參數(shù)都有4個不同方向的值,取其平均值作為方向無關(guān)的特征值。將采樣點(diǎn)的多光譜圖像輸入ENVI 5.3軟件,利用二階統(tǒng)計濾波工具計算2個波段灰度共生矩陣的8個特征參數(shù)。如圖5、6所示,分別為剔除土壤背景的植被冠層圖像(E1)和原始多光譜圖像(E2)處理的可見光示意圖,均值(MEA)、方差(VAR)、均勻性(HOM)、對比度(CON)、差異(DIS)、熵(ENT)、二階矩(SEC)和相關(guān)性(COR)特征參數(shù)灰度圖。從圖中可看出,灰度共生矩陣的特征參數(shù)表現(xiàn)出圖像豐富的紋理特征,剔除土壤背景前后的多光譜圖像紋理特征存在明顯差異,統(tǒng)計灰度圖的紋理特征參數(shù)數(shù)值,對紋理特征進(jìn)行進(jìn)一步分析。

圖5 剔除土壤背景多光譜圖像灰度共生矩陣特征參數(shù)灰度圖Fig.5 Gray scale of characteristic parameters of co-occurrence matrix of multispectral image with soil background eliminated

圖6 原始多光譜圖像灰度共生矩陣特征參數(shù)灰度圖Fig.6 Gray level co-occurrence matrix characteristic parameters gray level

2.4 圖像紋理特征

圖像紋理是用于識別圖像中感興趣的對象或區(qū)域的重要特征之一[36],為了對比剔除土壤背景以及不同鹽漬化等級下多光譜圖像紋理特征的差異,對E1和E2處理的多光譜圖像Band3的紋理特征參數(shù)進(jìn)行統(tǒng)計,計算2個處理下3個鹽漬化等級樣本點(diǎn)多光譜圖像的紋理特征參數(shù)的平均值,如表4所示。

表4 多光譜圖像紋理特征參數(shù)Tab.4 Multispectral image texture feature parameters

分別對比2個處理下3個鹽漬化等級的紋理特征參數(shù)平均值。整體上,相比剔除土壤背景前的E2,剔除土壤背景后的E1處理的紋理特征值中,SEC、CON、MEA特征值減小,DIS、VAR、ENT、COR、HOM特征值增大。其中CON、ENT、COR 3個參數(shù)的變化最顯著,CON減小表示圖像紋理溝紋變淺,圖像清晰度下降;ENT變大,表示圖像紋理增多;COR增大,表示圖像紋理變粗。

2.5 基于全子集回歸篩選光譜自變量

通過R語言全子集回歸模型算法,分別對E1和E2處理的 a (植被指數(shù))和b (植被指數(shù)和圖像紋理特征參數(shù)) 數(shù)據(jù)集進(jìn)行全子集變量篩選,篩選評價指標(biāo)選擇R2和 BIC 兩個指標(biāo)。綜合衡量R2最大和 BIC最小的變量組合,得到每個數(shù)據(jù)集中最優(yōu)的變量組,結(jié)果如表5所示。

表5 全子集篩選最優(yōu)變量組合統(tǒng)計Tab.5 Optimal variable combination statistics screened by full subset method

對比E1處理下2個變量種類的篩選結(jié)果,變量組a篩選得到的敏感變量組含有NDVI、SRVI、NDSI 3個植被指數(shù);變量組b的敏感變量組中包含NDVI、

SRVI、ARVI、TVI 4個植被指數(shù),篩選掉了紋理特征參數(shù)。對比E2處理下2個變量種類的篩選結(jié)果,變量組a篩選得到的敏感變量組包含DVI、EVI、SRVI 3個植被指數(shù);變量組b敏感變量組中包含RVI、NDSI兩個植被指數(shù)和ENT、COR、SEC 3個圖像紋理特征參數(shù)。

2.6 基于圖像紋理特征和植被指數(shù)的SSC反演模型

通過不同植被指數(shù)分類結(jié)果分析可以得出NDVI可以達(dá)到良好的圖像分類效果。因此,利用NDVI對多光譜圖像進(jìn)行分類,并對圖像分類得到的土壤像元進(jìn)行掩膜處理,剔除土壤像元,獲得純凈度更高的植被冠層圖像。分別將E1和E2輸入ENVI 5.3軟件提取6個波段的反射率。將多光譜可見光和近紅外1波段的反射率代入植被計算公式(表1),得到相應(yīng)的植被指數(shù),涉及的植被指數(shù)有NDVI、DVI、RVI、MSAVI、ARVI、EVI、CRSI。并基于全子集回歸變量模型篩選結(jié)果,構(gòu)建基于敏感變量組的SSC反演模型。

基于E1和E2處理下a、b兩種變量組全子集法篩選的變量結(jié)果,構(gòu)建基于敏感變量的SVM和ELM土壤含鹽量反演模型,各模型建模集和驗(yàn)證集的R2和NRMSE如表6所示。

表6 基于不同變量組的SVM、ELM反演模型精度Tab.6 SVM and ELM inversion model based on different variable groups

基于全子集篩選的2種敏感變量組,利用SVM和ELM兩種機(jī)器學(xué)習(xí)方法,構(gòu)建E1和E2處理共8個土壤含鹽量反演模型,驗(yàn)證集實(shí)測值和模擬值如圖7所示。

圖7 土壤含鹽量實(shí)測值和模擬值對比Fig.7 Comparison of measured and simulated soil salt contents

E1、E2處理下2個變量組的SVM和ELM機(jī)器學(xué)習(xí)土壤含鹽量反演模型中,ELM算法的E1-a、E1-b、E2-a、E2-b模型的R2分別為0.602、0.657、0.584、0.681,均大于對應(yīng)的SVM模型。ELM模型的標(biāo)準(zhǔn)均方根誤差分別為0.181、0.156、0.201、0.191,均小于SVM模型驗(yàn)證集的標(biāo)準(zhǔn)均方根誤差。綜上,ELM模型的表現(xiàn)更好。

3 討論

表7 剔除土壤背景前后紋理特征參數(shù)平均值Tab.7 Average value of image texture characteristic parameters before and after removing soil background

本研究僅基于無人機(jī)遙感多光譜反射率、圖像紋理特征與植被指數(shù)的土壤含鹽量反演,對于作物的耕種方式、作物不同生育期受土壤鹽漬化的影響,灌水量對土壤含鹽量的影響尚未考慮。本研究建立了特定情形下的土壤含鹽量反演模型,模型結(jié)果較好,但模型是否具有廣泛的適用性,不僅需要更加深入的研究,還需要綜合考慮更多因素,進(jìn)一步提高反演模型精度。

4 結(jié)論