向友珍 李汪洋 臺(tái) 翔 安嘉琪 王 辛 陳俊英

(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ū)土壤鹽漬化問(wèn)題嚴(yán)重,對(duì)鹽漬化土壤進(jìn)行治理和開發(fā)利用是農(nóng)業(yè)可持續(xù)發(fā)展與生態(tài)文明建設(shè)高質(zhì)量發(fā)展的有力保障[2]。實(shí)時(shí)、高效地監(jiān)測(cè)土壤鹽漬化狀況,是治理土壤鹽漬化的重要前提。

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

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

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

1 材料與方法

1.1 研究區(qū)概況

本研究區(qū)域位于巴彥淖爾市沙壕渠灌域(40°52′～41°00′N,107°05′～107°10′E),隸屬于河套灌區(qū)解放閘灌域,如圖1所示。灌域因不合理的灌排方式,氣候、土質(zhì)和地貌等因素的綜合影響,土壤鹽漬化問(wèn)題突出。區(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,多年平均日照時(shí)長(zhǎng)為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個(gè)土壤采樣點(diǎn)。采樣點(diǎn)布設(shè)見文獻(xiàn)[12]。

1.3 光譜圖像采集和處理

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

1.4 光譜指數(shù)計(jì)算

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

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

1.5 灰度共生矩陣法

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

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

(1)

(2)

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

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

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

1.6 Otsu算法

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

(3)

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

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

1.7 全子集篩選法

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

1.8 建模策略

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

1.9 模型方法和評(píng)價(jià)指標(biāo)

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

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

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

2 結(jié)果與分析

2.1 土壤含鹽量

對(duì)試驗(yàn)地總計(jì)87個(gè)采樣點(diǎn)的土壤含鹽量進(jìn)行統(tǒng)計(jì),將計(jì)算表層和深度10～20 cm實(shí)測(cè)土壤含鹽量的平均值作為深度10～20 cm的平均土壤含鹽量。將深度0～20 cm平均土壤含鹽量按照比例2∶1隨機(jī)劃分建模集和驗(yàn)證集,結(jié)果如表3所示。將各總集、建模集、驗(yàn)證集樣本點(diǎn)土壤含鹽量劃分為4個(gè)等級(jí):非鹽土(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)計(jì)分析Tab.3 Statistical analysis of soil salt content characteristics

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

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

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

基于近紅外和可見光波段構(gòu)建4個(gè)植被指數(shù),利用對(duì)土壤背景敏感的4個(gè)植被指數(shù)(NDVI、MSAVI、DVI、CRSI)進(jìn)行圖像分類,驗(yàn)證不同指數(shù)的圖像分類精度。為了對(duì)比4個(gè)植被指數(shù)圖像分類的精度,從87個(gè)采樣點(diǎn)中隨機(jī)選取40個(gè)樣本,進(jìn)行如圖4所示,4個(gè)植被指數(shù)的精度存在顯著差異,其中NDVI的圖像分類精度最高,40個(gè)樣本的總體精度均大于93.4%,均值為96.9%;Kappa系數(shù)在0.84～0.99之間,均值為0.92。MSAVI的分類精度略低于NDVI,總體精度和Kappa系數(shù)均值分別為95.8%和0.91;DVI和CRSI 2個(gè)植被指數(shù)的分類效果相對(duì)較差,總體精度均值小于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é)果評(píng)價(jià)。分別對(duì)40個(gè)樣本的多光譜圖像進(jìn)行2.1節(jié)的操作,獲得分類結(jié)果,并在ENVI 5.3中對(duì)分類結(jié)果進(jìn)行評(píng)價(jià),得到4個(gè)植被指數(shù)的總體精度和Kappa系數(shù)。

2.3 圖像紋理計(jì)算

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

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

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

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

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

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

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

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

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

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

通過(guò)不同植被指數(shù)分類結(jié)果分析可以得出NDVI可以達(dá)到良好的圖像分類效果。因此,利用NDVI對(duì)多光譜圖像進(jìn)行分類,并對(duì)圖像分類得到的土壤像元進(jìn)行掩膜處理,剔除土壤像元,獲得純凈度更高的植被冠層圖像。分別將E1和E2輸入ENVI 5.3軟件提取6個(gè)波段的反射率。將多光譜可見光和近紅外1波段的反射率代入植被計(jì)算公式(表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個(gè)土壤含鹽量反演模型,驗(yàn)證集實(shí)測(cè)值和模擬值如圖7所示。

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

E1、E2處理下2個(gè)變量組的SVM和ELM機(jī)器學(xué)習(xí)土壤含鹽量反演模型中,ELM算法的E1-a、E1-b、E2-a、E2-b模型的R2分別為0.602、0.657、0.584、0.681,均大于對(duì)應(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ī)H基于無(wú)人機(jī)遙感多光譜反射率、圖像紋理特征與植被指數(shù)的土壤含鹽量反演,對(duì)于作物的耕種方式、作物不同生育期受土壤鹽漬化的影響,灌水量對(duì)土壤含鹽量的影響尚未考慮。本研究建立了特定情形下的土壤含鹽量反演模型,模型結(jié)果較好,但模型是否具有廣泛的適用性,不僅需要更加深入的研究,還需要綜合考慮更多因素,進(jìn)一步提高反演模型精度。

4 結(jié)論