林 楠,劉海琪,楊佳佳,吳夢紅,劉翰霖
1. 吉林建筑大學(xué)測繪與勘查工程學(xué)院,吉林 長春 130118 2. 吉林大學(xué)地球科學(xué)學(xué)院,吉林 長春 130026 3. 東北大學(xué)資源與土木工程學(xué)院,遼寧 沈陽 110819 4. 中國地質(zhì)調(diào)查局沈陽地質(zhì)調(diào)查中心,遼寧 沈陽 110034
土壤養(yǎng)分含量快速估測對(duì)作物生長監(jiān)測和黑土科學(xué)管理具有重要的意義。 高光譜數(shù)據(jù)可以通過數(shù)千個(gè)緊密排列的波長記錄土壤信息。 眾多研究表明,土壤高光譜數(shù)據(jù)與土壤含水率[1]、養(yǎng)分含量[2-3]、重金屬含量[4]等土壤指標(biāo)存在密切關(guān)系,利用高光譜技術(shù)進(jìn)行土壤元素含量估測已取得了豐碩的研究成果。
目前,利用機(jī)器學(xué)習(xí)算法建立土壤高光譜特征和元素含量之間的反演模型是土壤元素含量估測的主要方法,其中支持向量機(jī)[5]、極限學(xué)習(xí)機(jī)[6]、隨機(jī)森林[7]等模型都得到了廣泛的應(yīng)用。 機(jī)器學(xué)習(xí)算法通常需要預(yù)先設(shè)定一些模型參數(shù)值,由于預(yù)定義的參數(shù)值很可能不包含全局最優(yōu)的參數(shù)值,導(dǎo)致機(jī)器學(xué)習(xí)模型達(dá)不到最佳效果。 為了克服機(jī)器學(xué)習(xí)模型在尋找最佳模型參數(shù)時(shí)所存在的問題,遺傳算法、粒子群算法等經(jīng)典優(yōu)化算法被用來優(yōu)化機(jī)器學(xué)習(xí)模型的內(nèi)部參數(shù)。 然而粒子群等經(jīng)典優(yōu)化算法對(duì)初始參數(shù)設(shè)置比較敏感,尋優(yōu)過程中容易陷入局部最優(yōu)解,導(dǎo)致算法后期的收斂速度變慢。 蝙蝠算法是一種新興的群體智能方法,在尋優(yōu)過程中,蝙蝠算法模仿蝙蝠聲波脈沖響度、頻率的自適應(yīng)調(diào)節(jié)過程,實(shí)現(xiàn)全局尋優(yōu)和局部尋優(yōu)過程的自由切換,從而使算法的全局搜索能力和局部搜索能力得到平衡,在模型參數(shù)尋優(yōu)計(jì)算方面表現(xiàn)優(yōu)異[8]。 本工作對(duì)黑土土壤有機(jī)質(zhì)、磷元素和鉀元素的光譜特征進(jìn)行詳細(xì)分析,選擇元素最佳光譜變換形式及相關(guān)性較高的敏感波段,將蝙蝠算法和Adaboost機(jī)器學(xué)習(xí)模型組合構(gòu)建土壤養(yǎng)分含量估測模型,利用蝙蝠算法求解Adaboost模型建模中的關(guān)鍵參數(shù),并比較模型參數(shù)優(yōu)化前后的估測精度,為黑土養(yǎng)分含量高光譜估測提供一種高效的新方法。
土壤樣品采集區(qū)位于黑龍江省訥河市,該區(qū)地處松遼平原的北端,大小興安嶺南緣,是黑土分布的典型區(qū)域。 區(qū)內(nèi)耕地面積廣闊,農(nóng)產(chǎn)品種類豐富,是甜菜和馬鈴薯的盛產(chǎn)地。 2019年4月在該區(qū)內(nèi)挑選典型黑土地塊,設(shè)置5公里大小的網(wǎng)格作為采樣單元,采樣點(diǎn)設(shè)計(jì)結(jié)合第二次全國土地調(diào)查成果合理分布,所布設(shè)的采樣點(diǎn)要能夠代表采樣單元的土壤性質(zhì),試驗(yàn)區(qū)內(nèi)共采集土壤樣本80個(gè)(圖1),采集樣品過程中確保采樣距離公路至少150 m,以每個(gè)采樣點(diǎn)位置為中心,在其周圍5 m×5 m范圍內(nèi)進(jìn)行樣品采集,每個(gè)采樣點(diǎn)位置共采集5份黑土樣本,取樣深度為表層土的15 cm以內(nèi),將樣本充分混合后裝入采樣袋。
圖1 研究區(qū)和樣品采集點(diǎn)地理位置Fig.1 Location of study area and sampling collection points
將樣品中的秸稈、砂礫等雜物剔除,風(fēng)干后研磨過篩,使得土壤的粒徑小于0.25 mm; 將樣品分成兩份,一份用于測定土壤元素含量,另一份用于室內(nèi)高光譜測量。 結(jié)合多目標(biāo)樣品分析質(zhì)量要求等技術(shù)規(guī)范,參照《土壤農(nóng)業(yè)化學(xué)分析方法》,選用重鉻酸鉀容量法測定土壤有機(jī)質(zhì)含量,選用X射線熒光光譜法進(jìn)行土壤磷元素和鉀元素含量測量[9]。 對(duì)元素含量的測量結(jié)果進(jìn)行統(tǒng)計(jì)分析(表1),將80份土壤樣品按元素含量由低到高分為20組,每組隨機(jī)抽取1個(gè)樣品放入驗(yàn)證集,共20個(gè)樣品作為驗(yàn)證集,其余60個(gè)樣品為訓(xùn)練集。
表1 土壤樣本養(yǎng)分元素含量基本信息Table 1 Statistics of nutrient contents in soil samples
在暗室中進(jìn)行土壤光譜測量,測量儀選用ASD FiledSpec4地物波譜儀,為提高光譜測量數(shù)據(jù)的精度,取5次光譜測量的均值作為土壤樣品的反射光譜數(shù)據(jù)。 由于光譜數(shù)據(jù)在350~399 nm以及2 451~2 500 nm范圍的噪聲較大,信噪比低,對(duì)土壤元素與反射率關(guān)系的分析會(huì)造成干擾,所以將400 nm前及2 450 nm后的數(shù)據(jù)剔除。 光譜儀采樣間隔為1 nm,即在400~2 450 nm范圍內(nèi)得到2 051個(gè)波段,由于光譜分辨率較高,波段數(shù)較多,相鄰波段之間可能存在信息重疊,更容易受到噪聲的影響,因此,將光譜數(shù)據(jù)進(jìn)行重采樣處理,設(shè)置采樣間隔為10 nm。 在去噪和重采樣處理的基礎(chǔ)上,將原始光譜反射率進(jìn)行一階微分、倒數(shù)對(duì)數(shù)等特征變換,不同的變換形式可以幫助準(zhǔn)確、快速地找到峰谷,并通過峰谷確定相應(yīng)的波長,從而確定敏感波段。
蝙蝠算法(Bat algorithm,BA)是一種啟發(fā)式搜索算法,它模擬蝙蝠使用聲納探測獵物和躲避障礙物,通過模擬蝙蝠飛行尋找獵物的過程來模擬優(yōu)化搜索過程,在計(jì)算過程中利用求解問題的適應(yīng)度值來選取蝙蝠的位置,利用優(yōu)勝劣汰的進(jìn)化過程來模擬較優(yōu)可行解代替較差可行解的迭代搜索過程[10]。 基于BA的基本原理,在算法各項(xiàng)參數(shù)被初始化后,從d維搜索空間中的一個(gè)隨機(jī)位置zl開始啟發(fā)式搜索。 以固定的頻率、不同的波長和音強(qiáng)搜索獵物,搜索過程中,蝙蝠根據(jù)接近獵物的距離自動(dòng)調(diào)整波長的大小。 經(jīng)全局搜索后更新每只蝙蝠的飛行速度和空間位置,并計(jì)算目標(biāo)函數(shù)的適應(yīng)度值,速度和空間位置更新公式如式(1)[11]
ft=fmin+(fmax-fmin)β
(1)
Adaboost算法基本思想是針對(duì)同一訓(xùn)練集訓(xùn)練能力一般的弱回歸器,通過疊加方法構(gòu)建一個(gè)訓(xùn)練能力很強(qiáng)的強(qiáng)回歸器。 其算法本身是通過改變樣本分布權(quán)重來實(shí)現(xiàn)的,根據(jù)每次訓(xùn)練集中每個(gè)樣本的預(yù)測精度以及上次的總體預(yù)測精度來計(jì)算每個(gè)弱回歸器的權(quán)重,同時(shí)更新每個(gè)樣本的分布權(quán)重,最后將每次訓(xùn)練得到的回歸器結(jié)果加權(quán)求和,作為強(qiáng)回歸器最后的輸出結(jié)果[12]。 建模過程中,弱學(xué)習(xí)器的最大迭代次數(shù)n和弱學(xué)習(xí)器的權(quán)重縮減系數(shù)v是兩個(gè)重要的參數(shù),迭代次數(shù)n設(shè)置過小會(huì)導(dǎo)致模型擬合不充分,n設(shè)置太大則會(huì)導(dǎo)致模型擬合過度,而較小的權(quán)重縮減系數(shù)v意味著需要更多的弱學(xué)習(xí)器的選代次數(shù),所以通常這兩個(gè)參數(shù)要一起進(jìn)行優(yōu)化調(diào)整[13]。
為了評(píng)價(jià)模型預(yù)測能力及模型的穩(wěn)定性,選擇決定系數(shù)(R2)和均方根誤差(RMSE)對(duì)模型的建模效果進(jìn)行評(píng)價(jià)
(2)
經(jīng)處理后的室內(nèi)土壤樣品光譜曲線如圖2(a)所示,從圖中可以看出,采集的黑土樣品原始反射率在0~0.9之間,各樣品光譜曲線波動(dòng)形狀相似,在可見光波段,隨著波長的增加,反射率逐漸增加,直到1 200 nm處反射率趨于穩(wěn)定。 近紅外區(qū)土壤光譜反射率總體高于可見光區(qū),兩個(gè)明顯的波谷分布在1 400和1 900 nm附近,主要是土壤中殘留的水分和空氣中水蒸氣的吸收造成的,而在2 200 nm處有輕微的凹陷,是受土壤中存在的粘土礦物的影響。 圖2(b)—(d)分別為原始反射率經(jīng)一階微分(R′)、倒數(shù)對(duì)數(shù)(lg1/R)、倒數(shù)對(duì)數(shù)一階微分[(lg1/R)′]變換后的光譜曲線,從變化結(jié)果可以看出,一階微分變換可以對(duì)原始光譜變化起到放大作用,經(jīng)變換后反射率在1 400,1 900和2 200 nm處的波動(dòng)更大。
圖2 反射率及其變換的土壤光譜曲線Fig.2 Soil spectral reflectance curves and its transformations
分別計(jì)算了土壤有機(jī)質(zhì)、磷、鉀含量與土壤反射率的相關(guān)系數(shù),并繪制相關(guān)系數(shù)曲線圖(圖3)。 從圖3(a)原始光譜的相關(guān)系數(shù)可以看出,土壤有機(jī)質(zhì)和磷元素含量與光譜反射率呈負(fù)相關(guān),而鉀元素含量則相反。 與原始光譜反射率相比,變換后的光譜數(shù)據(jù)與土壤有機(jī)質(zhì)、磷、鉀含量的相關(guān)性更高,其中,一階微分變換形式與土壤有機(jī)質(zhì)、磷、鉀含量的相關(guān)系數(shù)表現(xiàn)為正、負(fù)交叉,波峰和波谷較多,而且經(jīng)一階微分變換后,各元素的最高相關(guān)系數(shù)顯著提高。
圖3 土壤元素含量與光譜相關(guān)系數(shù)曲線Fig.3 Correlation coefficient curves of soil element content and spectral reflectance
選取相關(guān)系數(shù)大于0.4的敏感波段作為預(yù)測模型的樣本輸入數(shù)據(jù)(表2)。 由統(tǒng)計(jì)結(jié)果可以看出,經(jīng)過不同的光譜特征變換,土壤元素含量與光譜反射率的相關(guān)系數(shù)有所提高,其中,土壤有機(jī)質(zhì)對(duì)應(yīng)的最佳變換形式為一階微分,與光譜反射率數(shù)據(jù)的相關(guān)系數(shù)最高為0.796,波長范圍在1 370 mm附近; 土壤磷元素對(duì)應(yīng)的最佳變換形式為倒數(shù)對(duì)數(shù)的一階微分,相關(guān)系數(shù)最高為-0.688,波長范圍在690 mm附近; 土壤鉀元素對(duì)應(yīng)的最佳變換形式為倒數(shù)對(duì)數(shù)一階微分,相關(guān)系數(shù)最高為0.553,波長范圍在980 mm附近。
表2 最大相關(guān)系數(shù)和敏感波段Table 2 Maximum correlation coefficients and sensitive bands
將選取60個(gè)訓(xùn)練樣本最佳光譜變換形式的敏感波段和對(duì)應(yīng)元素含量值作為Adaboost模型的建模數(shù)據(jù),3種元素的建模波段數(shù)分別為85個(gè)、82個(gè)和51個(gè),選擇CART決策樹為Adaboost模型的弱回歸學(xué)習(xí)器,采用蝙蝠搜索算法對(duì)Adaboost模型的最大迭代次數(shù)n和弱學(xué)習(xí)器權(quán)重縮減系數(shù)v進(jìn)行尋優(yōu)。
基于BA的基本原理,首先需對(duì)模型各項(xiàng)參數(shù)進(jìn)行初始化,BA初始化參數(shù)較多,但除了迭代次數(shù)外其他參數(shù)敏感性不強(qiáng),可選擇默認(rèn)缺省參數(shù)進(jìn)行初始化。 BA算法的搜索空間是由最大迭代次數(shù)n和弱學(xué)習(xí)器權(quán)重縮減系數(shù)v為坐標(biāo)軸組成的二維空間,迭代搜索過程從搜索空間內(nèi)的L個(gè)隨機(jī)位置開始,在每次迭代過程中,利用每只蝙蝠在空間中位置的二維坐標(biāo)(n,v)作為Adaboost模型的初始化參數(shù),然后基于Adaboost模型對(duì)樣本數(shù)據(jù)訓(xùn)練建模,并計(jì)算預(yù)測結(jié)果。 根據(jù)計(jì)算結(jié)果選取R2值最大時(shí)對(duì)應(yīng)的位置作為蝙蝠當(dāng)前最優(yōu)位置,利用式(1)更新每只蝙蝠的空間位置。 設(shè)置種群大小L=20,脈沖頻率范圍fmin=0、fmax=1、脈沖音強(qiáng)范圍Amin=0,Amax=1,脈沖響度衰減系數(shù)α=0.9,脈沖頻度增加系數(shù)γ=0.9,作為缺省參數(shù),并計(jì)算了不同迭代次數(shù)對(duì)應(yīng)的R2值(圖4)。 從圖中可以看出,隨著迭代次數(shù)T的增加,R2值逐漸增大,3種元素中,鉀元素收斂速度最快,當(dāng)?shù)鷶?shù)達(dá)到10次的時(shí)候,R2值達(dá)到了最大,有機(jī)質(zhì)和磷元素分別是當(dāng)?shù)鷶?shù)達(dá)到12次和16次達(dá)到最大。 建模估測時(shí),選取R2最大時(shí)對(duì)應(yīng)的n和v的值作為Adaboost模型的建模參數(shù)。
圖4 BA-Adaboost模型R2變化曲線Fig.4 R2 variation curves of BA-Adaboost
表3 模型預(yù)測精度對(duì)比分析Table 3 Comparison of estimation accuracies
利用模型的預(yù)測值和實(shí)測值繪制預(yù)測結(jié)果擬合圖,進(jìn)一步對(duì)比分析3種元素模型擬合效果(圖5)。 從圖中可以看出,與Adaboost相比,優(yōu)化后的BA-Adaboost模型的預(yù)測值更緊密地分布在1∶1左右,數(shù)據(jù)擬合能力和穩(wěn)定性優(yōu)于Adaboost模型。 同時(shí)可以看出訓(xùn)練集的擬合效果比驗(yàn)證集更優(yōu)異,驗(yàn)證集樣本點(diǎn)分布相對(duì)分散,說明機(jī)器學(xué)習(xí)模型對(duì)訓(xùn)練樣本數(shù)據(jù)的可靠性要求較高。
以黑龍江省訥河市80個(gè)黑土樣品為研究對(duì)象,對(duì)黑土土壤有機(jī)質(zhì)、磷元素和鉀元素的光譜特征進(jìn)行分析,構(gòu)建BA-Adaboost模型對(duì)元素含量進(jìn)行估測,主要結(jié)論如下:
(1)光譜數(shù)據(jù)的特征變換處理可以增強(qiáng)光譜特征,一階微分變化有效地突出了光譜曲線的波峰和波谷,提高了光譜反射率與元素含量的相關(guān)系數(shù),有機(jī)質(zhì)、磷元素、鉀元素分別在1 370,690和980 mm處相關(guān)系數(shù)達(dá)到最大值。
圖5 元素含量實(shí)測值和預(yù)測值(a): OM含量Adaboost模型; (b): OM含量BA-Adaboost模型; (c): P含量Adaboost模型; (d): P含量BA-Adaboost模型; (e): K含量Adaboost模型; (f): K含量BA-Adaboost模型Fig.5 Measured and predicted values of element content(a): Adaboost of OM content; (b): BA-Adaboost of OM content; (c): Adaboost of P content; (d): BA-Adaboost of P content; (e): Adaboost of K content; (f): BA-Adaboost of K content
(2)蝙蝠算法在參數(shù)尋優(yōu)過程中能夠動(dòng)態(tài)控制全局尋優(yōu)和局部尋優(yōu)過程之間的自由切換,利用蝙蝠算法對(duì)Adaboost模型參數(shù)進(jìn)行優(yōu)化計(jì)算,避免了估測模型陷入局部極小值,與傳統(tǒng)參數(shù)尋優(yōu)方法相比,該方法可更快更準(zhǔn)確地找到全局最優(yōu)參數(shù)。
(3)將BA與Adaboost模型相結(jié)合,構(gòu)建BA-Adaboost土壤含量估測模型,該組合模型只需設(shè)置搜索空間,然后自動(dòng)搜索模型最優(yōu)參數(shù)值; 對(duì)比BA算法優(yōu)化前后的估測精度可以看出,優(yōu)化后BA-Adaboost模型R2增大、RMSE變小,估測精度明顯提高,說明BA-Adaboost模型在土壤元素含量高光譜估測中具有一定的適用性,擴(kuò)展了機(jī)器學(xué)習(xí)模型在土壤成分估測研究中的應(yīng)用。