陳珠琳 王雪峰

(中國(guó)林業(yè)科學(xué)研究院資源信息研究所 北京 100091)

降香黃檀(Dalbergiaodorifera)是海南省特有的極珍貴用材樹(shù)種，屬于國(guó)家二級(jí)保護(hù)植物(鄭堅(jiān)等,2016;苗靈鳳,2017),由于具有極高的用材價(jià)值，人們對(duì)其進(jìn)行肆意破壞，造成資源急劇減少(楊菁等,2015)。近些年來(lái)，海南省大力推廣種植珍貴樹(shù)種，但由于缺乏先進(jìn)的管理監(jiān)測(cè)技術(shù)，推廣的效果并不理想。因此，為提高降香黃檀的生長(zhǎng)質(zhì)量與結(jié)香速率，必須提出方便快捷的營(yíng)養(yǎng)診斷方法。傳統(tǒng)的營(yíng)養(yǎng)診斷法需要實(shí)驗(yàn)室化學(xué)測(cè)定，耗費(fèi)大量人力物力財(cái)力(賈良良等,2009)。近些年來(lái)，機(jī)器自主營(yíng)養(yǎng)診斷發(fā)展迅速。主要方法有：數(shù)字圖像處理、手持光譜儀和高光譜成像。從數(shù)據(jù)的豐富度來(lái)講，高光譜成像技術(shù)無(wú)疑是最佳的選擇，但不適用于小尺度監(jiān)測(cè)；手持光譜儀技術(shù)可以獲取400～2 500 nm之間的植被光譜反射率，但數(shù)據(jù)之間的冗余度大，易出現(xiàn)休斯現(xiàn)象；數(shù)字圖像處理技術(shù)通常指普通的可見(jiàn)光圖像，但由于波段范圍窄且?guī)捿^大而受限。相比之下，多光譜相機(jī)波段比普通相機(jī)豐富，也可以獲取各波段的紋理、形狀等信息，可用于地面和航空攝影多種尺度，為植物營(yíng)養(yǎng)診斷提供一種獲取數(shù)據(jù)的新方法。

過(guò)量的金屬元素富集在土壤中，給植物的光合作用帶來(lái)不良影響。某些情況下，植被吸收的金屬離子會(huì)取代葉綠素分子中間的鎂原子，減小了聚光并導(dǎo)致光合作用崩潰。重金屬污染使植被光譜特性發(fā)生明顯改變，因此植物的生理敏感現(xiàn)象能夠?qū)饘傥廴咀龀隹焖俚纳砩鷳B(tài)反應(yīng)(屈永華等,2015)。例如，隨著銅鋅含量的增加，植物在可見(jiàn)光波段的光譜反射率上升，近紅外區(qū)間的光譜反射率下降，紅邊出現(xiàn)“藍(lán)移”(王慧等,2017;陳思寧等,2007)。但由于植物種類不同等原因，重金屬脅迫導(dǎo)致的光譜反射率特征變化并沒(méi)有相對(duì)一致的結(jié)論(Horleretal.,1980;Cleversetal.,2004)。

鐵(Fe)是植物生長(zhǎng)發(fā)育必不可少的微量元素之一，直接參與機(jī)體的生理生化活動(dòng)(袁程等,2012)。盡管土壤中的鐵元素總量不低，但在堿性環(huán)境中，F(xiàn)e2+離子含量低，會(huì)造成植物缺鐵(周曉今等,2012)；而當(dāng)土壤呈現(xiàn)出弱酸性時(shí)，植物又會(huì)因“鐵毒”出現(xiàn)生理病變。目前，有關(guān)葉片金屬元素含量預(yù)測(cè)的研究存在以下不足：1) 大量的研究均將植物葉片或冠層光譜反射率作為自變量，忽視了對(duì)紋理特征的分析。2) 輸入變量的相關(guān)性對(duì)模型精度的影響較大，如何減小自變量之間的冗余也是建模的重點(diǎn)。3) 在評(píng)價(jià)效果相似的模型時(shí)，忽略了與植物營(yíng)養(yǎng)需求的最佳區(qū)間相結(jié)合，不利于實(shí)際應(yīng)用。鑒于此，本文以幼齡降香黃檀為試驗(yàn)對(duì)象，提取多光譜圖像特征，主要分析葉片全鐵含量與光譜反射率變化，以及植被指數(shù)、各波段圖像紋理特征與全鐵含量之間的關(guān)系，提出合理的葉片全鐵含量預(yù)測(cè)方法，為珍貴樹(shù)種重金屬營(yíng)養(yǎng)診斷提供參考。

1 材料與方法

1.1 試驗(yàn)地概況

研究區(qū)設(shè)立在海南省文昌市龍樓鎮(zhèn)的島東林場(chǎng)昌灑作業(yè)區(qū)(19°43′58″—19°44′58″N，110°57′34″—111°01′54″E)，該林場(chǎng)占地約1.9萬(wàn)hm2，屬沿海平原地帶，熱帶海洋性氣候，年平均雨量為1 808.8 mm，土壤類型為濱海沙壤土，試驗(yàn)地土壤pH 5.0～6.6，有效氮98.3～114.8 mg·kg-1，有效磷3.38～4.56 mg·kg-1，速效鉀69.9～78.2 mg·kg-1，有效鐵2.33～4.89 mg·kg-1。

1.2 試驗(yàn)設(shè)計(jì)

試驗(yàn)采用2年生的降香黃檀實(shí)生苗，平均苗高為1.07 m，2017年3月選取健壯且生長(zhǎng)一致的幼樹(shù)移植至島東林場(chǎng)。定植2周后施加等量的水溶性復(fù)合肥進(jìn)行培養(yǎng)，1個(gè)月后進(jìn)行脅迫處理。將試驗(yàn)地劃分為4個(gè)正方形小區(qū)，間隔為5 m，小區(qū)面積為25 m2，每個(gè)小區(qū)內(nèi)種植25株降香黃檀。單因素試驗(yàn)設(shè)計(jì)，鐵源為乙二胺二鄰羥苯基大乙酸鐵鈉(EDDHA-FeNa)，采用根部施埋法，梯度分別為0、10、15、20 g每株(對(duì)應(yīng)4個(gè)小區(qū))，記為CK(對(duì)照組)、F1(低鐵)、F2(中鐵)、F3(高鐵)。

1.3 樣本采集及多光譜圖像獲取

試驗(yàn)于2017年4月17日和7月2日進(jìn)行鐵脅迫試驗(yàn)，9月1日進(jìn)行取樣，由于植物體內(nèi)養(yǎng)分的運(yùn)輸機(jī)制導(dǎo)致不同高度的葉片之間養(yǎng)分含量有一定的差異，所以本研究取植株的上中下層，并兼顧內(nèi)外層，進(jìn)行后期圖像獲取和全鐵含量化驗(yàn)。

試驗(yàn)使用美國(guó)生產(chǎn)的MicaSense RedEdge3多光譜相機(jī)獲取葉片圖像。該相機(jī)有5個(gè)波段(B、G、R、RE和NIR)，其中心波段分別為475、560、668、717和840 nm；帶寬分別為20、20、10、40和10 nm。將葉片放入暗箱中，光源為對(duì)應(yīng)中心波段的單波段LED燈珠組成的燈排，如圖1所示。獲取葉片圖像前進(jìn)行白板和黑板校正，各通道的光譜反射率計(jì)算公式為：

(1)

式中：IC為校正后的反射率，O代表目標(biāo)物的顏色均值，B代表黑板的顏色均值，W代表白板的顏色均值。

獲取圖像后，將葉片按編號(hào)裝入信封并放入烘箱內(nèi)，105 ℃下烘干0.5 h進(jìn)行殺青，然后75 ℃烘干至恒質(zhì)量。將烘干后的葉片磨成粉狀，采用原子吸收分光光度法測(cè)定全鐵含量(簡(jiǎn)稱TIC)。

1.4 自變量提取

葉片圖像使用PhotoShop 5.0進(jìn)行處理，之后提取不同數(shù)據(jù)源下的圖像參數(shù)。其中包括：植被指數(shù)(以下簡(jiǎn)稱VI)，如表1所示；多光譜圖像的各通道紋理特征均值(texture feature mean value,TFMV)與紋理特征方差(texture feature variance,TFV)，包括均值、方差、同質(zhì)性、對(duì)比度、相異性、熵值、二階矩、相關(guān)性。

圖1 攝影方法Fig.1 Screen method

表1 植被指數(shù)類型及與降香黃檀葉片全氮含量相關(guān)性分析①Tab.1 Vegetation indexes and correlation analysis with total nitrogen content in sandalwood leaves

①：*代表在0.05水平上顯著相關(guān)；**代表在0.01水平上顯著相關(guān)。下同。* represent significant correlation at 0.05 level;** represent significant correlation at 0.01 level.The same below.

1.5 變量篩選方法

當(dāng)數(shù)學(xué)模型的輸入自變量多且不相互獨(dú)立時(shí)，會(huì)出現(xiàn)過(guò)擬合現(xiàn)象，從而導(dǎo)致建立的模型精度低、建模時(shí)間長(zhǎng)等問(wèn)題。因此，在建立模型之前，需要將冗余的變量去除，選擇最能反映輸入輸出關(guān)系的自變量參與建模。本研究使用4種變量篩選方法：相關(guān)性分析法(correlation analysis，CA)、主成分分析法(principal component analysis,PCA)、平均影響值法(mean impact value,MIV)和遺傳算法(genetic algorithm,GA)，并比較分析其特點(diǎn)(李哲等,2018;顏勝科等,2017)。

1.6 PSO-BPNN算法

BPNN(back propagation neural network)是目前應(yīng)用最廣的神經(jīng)網(wǎng)絡(luò)技術(shù)(賈方方等,2016)，具有自適應(yīng)性、自學(xué)習(xí)、自組織和非線性映射的功能，適合于模擬各種錯(cuò)綜復(fù)雜的非線性關(guān)系(朱云芳等,2017)，同時(shí)也存在收斂速度慢、易陷入局部收斂、網(wǎng)絡(luò)的學(xué)習(xí)和記憶能力不穩(wěn)定等缺點(diǎn)(馬廉潔等,2016)。粒子群優(yōu)化算法(particle swarm optimization,PSO)具有參數(shù)少、易實(shí)現(xiàn)、精度高、搜索速度快等優(yōu)點(diǎn)，對(duì)多峰問(wèn)題、非線性均有較好的全局搜索能力(劉耿耿等,2018)，結(jié)合BPNN可以解決其學(xué)習(xí)速度慢以及局部收斂性不足等缺點(diǎn)。將二者相結(jié)合，既具有適應(yīng)性強(qiáng)、泛化能力好的優(yōu)點(diǎn)，又能減小陷入過(guò)擬合的概率，加強(qiáng)預(yù)測(cè)能力。

1.7 模型驗(yàn)證

本研究按照2∶1的比例隨機(jī)抽取擬合數(shù)據(jù)和驗(yàn)證數(shù)據(jù)。隨機(jī)抽取65個(gè)作為擬合數(shù)據(jù)，35個(gè)作為驗(yàn)證數(shù)據(jù)。其中模型的擬合優(yōu)度評(píng)價(jià)指標(biāo)為：

(2)

(3)

獨(dú)立樣本檢驗(yàn)的統(tǒng)計(jì)指標(biāo)為：

2 結(jié)果與分析

2.1 不同鐵水平下降香黃檀植株生長(zhǎng)量變化

整個(gè)試驗(yàn)周期下降香黃檀在樹(shù)高、冠幅、地莖生長(zhǎng)量的變化，結(jié)果如圖2所示。與對(duì)照組相比，施加鐵肥不同程度地促進(jìn)了植株的生長(zhǎng)。從CK到F2變化過(guò)程中，樹(shù)高、冠幅、地莖的增長(zhǎng)量是持續(xù)上升的。但在F3梯度下，樹(shù)高增長(zhǎng)速度與F2相比明顯降低，冠幅的增長(zhǎng)量則小于F2，但地莖的增長(zhǎng)較為明顯，增長(zhǎng)量是F2的2倍以上。試驗(yàn)結(jié)果說(shuō)明適當(dāng)添加鐵肥有助于提高生物量，但過(guò)量后樹(shù)高、冠幅生長(zhǎng)減緩，地莖增大。

圖2 不同鐵肥梯度下降香黃檀樹(shù)高、冠幅、地徑生長(zhǎng)量Fig.2 Growth of height,crown width and stem of Dalbergia odorifera with different iron fertilizer gradientsCK:對(duì)照組；F1:低鐵梯度；F2:中鐵梯度；F3:高鐵梯度。CK:Control group;F1:Low level of iron fertilization gradient;F2:Middle level of fertilization gradient;F3:High level of fertilization gradient.

2.1 光譜反射率隨TIC的變化趨勢(shì)

將降香黃檀葉片5個(gè)波段的多光譜反射率連接成折線圖，如圖3所示。NIR波段差異明顯，隨著施肥量的增加，光譜反射率呈現(xiàn)上升趨勢(shì)。可見(jiàn)光波段(B、G、R、RE)變化量相對(duì)較小，其大小無(wú)法識(shí)別，圖3中將其進(jìn)行了放大，如灰色矩形框中所示，各波段中不同脅迫對(duì)應(yīng)的反射率排名不同，例如，與其他脅迫組相比，CK在B波段的反射率最高，在G、RE和NIR波段得到的反射率最低。

表2為圖2中各點(diǎn)對(duì)應(yīng)的值。B波段光譜反射率在CK～F2梯度持續(xù)下降，但在F3梯度出現(xiàn)上升；G波段則相反，即在CK～F2梯度持續(xù)上升，在F3梯度下降；R波段反射率則在CK～F1梯度下降，F(xiàn)1～F2梯度上升，F(xiàn)2～F3基本無(wú)變化；RE波段則連續(xù)上升。這是由于在合理的鐵濃度下，葉片葉綠素含量升高，光合作用加強(qiáng)，對(duì)B、R、RE波段的吸收加強(qiáng)，G波段的反射率增強(qiáng)；而當(dāng)全鐵出現(xiàn)過(guò)量時(shí)，植物葉綠體遭到破壞，葉綠素含量降低，則出現(xiàn)與之相反的結(jié)果。近紅外波段反射率受細(xì)胞結(jié)構(gòu)的影響，隨著施鐵量的增加，植株的生長(zhǎng)更加旺盛，細(xì)胞分裂與分化能力增強(qiáng)，葉片增大，細(xì)胞間空隙增大，所以近紅外波段反射率持續(xù)增加。

圖3 不同全鐵含量下降香黃檀葉片光譜反射特征Fig.3 Spectral reflectance of D.odorifera leaves in different TIC

2.2 TIC與圖像參數(shù)的關(guān)系

表1中的皮爾森系數(shù)為由不同波段計(jì)算得到的植被指數(shù)與TIC之間的Pearson系數(shù)，與其他波段相比，由RE和NIR波段計(jì)算得到的植被指數(shù)最能夠反映降香黃檀葉片TIC，除MCARIRE和MTVRE外，其他均與TIC在0.01水平上顯著相關(guān)；而G和NIR得到的植被指數(shù)與TIC的現(xiàn)象關(guān)系最不明顯，僅有4個(gè)指數(shù)與TIC在0.05水平上相關(guān)。R-NIR組合及B-NIR組合得到的植被指數(shù)均有7個(gè)與TIC在0.05水平上相關(guān)，但B-NIR組合在數(shù)值上略占優(yōu)勢(shì)。

本文對(duì)葉片TIC與各波段TFMV與TFV進(jìn)行顯著性檢驗(yàn)，結(jié)果見(jiàn)表3和表4。TFMV與TIC在0.01水平上相關(guān)的有4個(gè)，在0.05水平上相關(guān)的有11個(gè)；而TFV與TIC在0.01水平上相關(guān)的有10個(gè)，在0.05水平上相關(guān)的有12個(gè)，這說(shuō)明變量的離散程度更能反映出葉片TIC值。由表3可知，從波段角度可以得出，RE波段和NIR波段呈現(xiàn)出的相關(guān)性最好，二者均在方差和對(duì)比度值方面與TIC呈現(xiàn)0.01水平相關(guān)，在相異性方面與TIC呈現(xiàn)0.05水平相關(guān)。從特征參數(shù)角度分析則得出同樣的結(jié)論，即方差和對(duì)比度值相關(guān)性最強(qiáng)，其次為相異性。由表4可知：從波段角度得到，NIR波段呈現(xiàn)出的相關(guān)性最好，分別在均值、方差、對(duì)比度、相異性方面呈現(xiàn)出0.01水平相關(guān)。其次是R和RE波段，G波段中呈現(xiàn)相關(guān)的參數(shù)最多，但均分布在0.05水平。從特征角度分析得到，均值效果最佳，在5個(gè)波段中均與TIC具有相關(guān)性，其次為相異性、同質(zhì)性、方差。

綜合表1、表3和表4的結(jié)論可以得出，RE和NIR波段是指示降香黃檀葉片TIC的重要指示波段，在植被指數(shù)以及紋理特征方面均優(yōu)于其他波段。

表2 不同鐵脅迫下降香黃檀葉片各波段光譜反射率①Tab.2 Spectral reflectance of leaves of D.odorifera under different iron stresses

① CK:對(duì)照組；F1:低鐵梯度；F2:中鐵梯度；F3:高鐵梯度。CK:Control group;F1:Low level of iron fertilization gradient;F2:Middle level of fertilization gradient;F3:High level of fertilization gradient.

表3 不同波段下全鐵含量與紋理特征均值的顯著性檢驗(yàn)①Tab.3 Significance test between total iron content and mean texture factors in different bands

①NIR:Near IR 波段；RE:紅邊波段。下同。NIR:Near IR band;RE denotes red edge band.The same below.

表4 不同波段下全鐵含量與紋理特征方差的顯著性檢驗(yàn)Tab.4 Significance test between total iron content and variance texture factors in different bands

2.3 變量相關(guān)性分析及篩選

通過(guò)公式及圖4A可以看出(0～30分別對(duì)應(yīng)表1中的植被指數(shù)，順序?yàn)镽VIB、RVIG、RVIR、RVIRE、DVIB…，以此類推)，前15個(gè)VIs之間的相關(guān)性較大，即變量之間具有較強(qiáng)的線性相關(guān)性。而后16個(gè)VIs之間的Pearson系數(shù)多數(shù)集中在0.7以下，變量之間的相關(guān)性較小。同理對(duì)與TIC在0.01、0.05水平上相關(guān)的15個(gè)TFMV，以及與TIC在0.01水平上相關(guān)的10個(gè)TFV分別進(jìn)行相關(guān)性分析，結(jié)果如圖4(b)-(c)所示。圖4B中15個(gè)TFMV之間的相關(guān)性較大，圖像呈現(xiàn)出的顏色較深，大部分區(qū)域相關(guān)性在0.8以上；而圖4C中對(duì)角線附近的顏色較深，即說(shuō)明部分相鄰參數(shù)之間的相關(guān)性較大。

圖4 參數(shù)相關(guān)性分析Fig.4 Correlation analysis of parameters

為了減少數(shù)據(jù)之間的冗余，本文首選選擇出與TIC在0.05和0.01水平上相關(guān)的變量，再使用4種變量篩選方法進(jìn)行二次選擇。其中相關(guān)性分析法中，將VIs和紋理特征的相關(guān)性閾值分別設(shè)置為0.7和0.8(即比較相關(guān)性大于閾值的變量，選擇與TIC的Pearson系數(shù)最高的變量)；主成分分析保留前三主成分；MIV算法中選擇在原值的基礎(chǔ)上加減10%構(gòu)成新的訓(xùn)練樣本；GA算法中染色體長(zhǎng)度為30，種群大小為20，最大進(jìn)化代數(shù)設(shè)置為100。為避免不同自變量個(gè)數(shù)對(duì)最終預(yù)測(cè)結(jié)果的影響，本研究選擇每種方法篩選出的前8位變量，最終結(jié)果如表5所示。

表5 變量篩選結(jié)果①Tab.5 Texture parameter selection results

①VIs代表植被指數(shù)；TFMV代表紋理特征均值；TFV代表紋理特征方差。下同。VIs represents vegetation indexes;TFMV represents texture feature mean value;TFV represents texture feature variance value.The same below.

2.4 PSO-BP神經(jīng)網(wǎng)絡(luò)模型構(gòu)建

對(duì)每一組試驗(yàn)構(gòu)建PSO-BP神經(jīng)網(wǎng)絡(luò)模型，過(guò)程如下：

2) 模型優(yōu)化。在BP神經(jīng)網(wǎng)絡(luò)中內(nèi)嵌一個(gè)種群規(guī)模為20，迭代次數(shù)為100的PSO算法。為提高算法的全局收斂能力，在PSO算法進(jìn)化方程中加入慣性權(quán)重因子，設(shè)置初始慣性權(quán)重為0.9，結(jié)束慣性權(quán)重為0.1，使算法在迭代初期保持較強(qiáng)的全局搜索能力，在迭代后期能進(jìn)行更精確的局部開(kāi)發(fā)。

2.5 全鐵含量預(yù)測(cè)結(jié)果分析

表6 不同試驗(yàn)下PSO-BP神經(jīng)網(wǎng)絡(luò)模型預(yù)測(cè)結(jié)果①Tab.6 Prediction results of PSO-BP neural network model in different tests

①CA：相關(guān)性分析法；PCA：主成分分析法；MIV：平均影響值法；GA：遺傳算法；PSO：粒子群優(yōu)化算法。下同。CA：Correlation analysis;PCA：Principal component analysis;MIV：Mean impact value;GA：Genetic algorithm;PSO：Particle swarm optimization.The same below.

圖5 CA與GA方法實(shí)測(cè)預(yù)測(cè)值Fig.5 Measured and predicted value of CA and GA method

2.6 不同特征組合對(duì)模型預(yù)測(cè)結(jié)果的影響

由表1～3可知，植被指數(shù)對(duì)擬合效果的貢獻(xiàn)率比較大，所以該參數(shù)是必不可少的自變量。本研究以CA和GA為例，設(shè)計(jì)4組對(duì)比試驗(yàn)進(jìn)行分析，驗(yàn)證紋理特征方面，TFMV和TFV對(duì)模型的精度是否有較大程度的提高，結(jié)果如表7所示。僅使用VI建模并不能很好的表明降香黃檀葉片TIC，加入紋理特征可以提高預(yù)測(cè)精度，并且加入TFV的效果要優(yōu)于TFMV，并且加入2種紋理特征后，各指標(biāo)提高的幅度較小，這說(shuō)明紋理特征的離散程度更能表達(dá)全鐵對(duì)葉片的影響。

表7 不同特征組合下的預(yù)測(cè)結(jié)果Tab.7 Prediction results of different feature combination

3 結(jié)論

以海南省島東林場(chǎng)種植的2年生降香黃檀為研究對(duì)象，使用多光譜相機(jī)獲取葉片圖像，并建立了反演葉片全鐵含量的神經(jīng)網(wǎng)絡(luò)模型，得到以下結(jié)論：適當(dāng)添加鐵肥有助于提高降香黃檀生物量，但過(guò)量后出現(xiàn)樹(shù)高、冠幅生長(zhǎng)減緩，地莖增大的現(xiàn)象。隨著葉片含鐵量的增大，可見(jiàn)光波段(B、G、R和RE)的反射率變化不明顯，NIR波段反射率持續(xù)增大。但除NIR波段外，可見(jiàn)光的RE波段也是指示鐵脅迫的有效波段，通過(guò)上述2個(gè)波段計(jì)算得到的植被指數(shù)和紋理參數(shù)與全鐵含量的相關(guān)性最大。在變量篩選時(shí)，使用相關(guān)性分析法得到的結(jié)果作為自變量得到的預(yù)測(cè)模型在150～300 mg·kg-1區(qū)間更接近真實(shí)值，而其他方法在該區(qū)間得到的結(jié)果偏低，若用于施肥指導(dǎo)，則會(huì)增大鐵毒產(chǎn)生的概率。另外，僅使用植被指數(shù)得到的模型預(yù)測(cè)結(jié)果較差，加入紋理特征后很明顯的提高了擬合優(yōu)度及預(yù)測(cè)精度，而且紋理特征方差對(duì)模型精度的影響更大，即片紋理的離散程度可以較好的作為預(yù)測(cè)全鐵含量的輔助信息。