基于植被指數及多光譜紋理特征的降香黃檀葉片全鐵含量預測*

陳珠琳 王雪峰

(中國林業科學研究院資源信息研究所 北京 100091)

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

過量的金屬元素富集在土壤中，給植物的光合作用帶來不良影響。某些情況下，植被吸收的金屬離子會取代葉綠素分子中間的鎂原子，減小了聚光并導致光合作用崩潰。重金屬污染使植被光譜特性發生明顯改變，因此植物的生理敏感現象能夠對金屬污染做出快速的生理生態反應(屈永華等,2015)。例如，隨著銅鋅含量的增加，植物在可見光波段的光譜反射率上升，近紅外區間的光譜反射率下降，紅邊出現“藍移”(王慧等,2017;陳思寧等,2007)。但由于植物種類不同等原因，重金屬脅迫導致的光譜反射率特征變化并沒有相對一致的結論(Horleretal.,1980;Cleversetal.,2004)。

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

1 材料與方法

1.1 試驗地概況

研究區設立在海南省文昌市龍樓鎮的島東林場昌灑作業區(19°43′58″—19°44′58″N，110°57′34″—111°01′54″E)，該林場占地約1.9萬hm2，屬沿海平原地帶，熱帶海洋性氣候，年平均雨量為1 808.8 mm，土壤類型為濱海沙壤土，試驗地土壤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 試驗設計

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

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

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

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

(1)

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

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

1.4 自變量提取

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

圖1 攝影方法Fig.1 Screen method

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

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

1.5 變量篩選方法

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

1.6 PSO-BPNN算法

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

1.7 模型驗證

本研究按照2∶1的比例隨機抽取擬合數據和驗證數據。隨機抽取65個作為擬合數據，35個作為驗證數據。其中模型的擬合優度評價指標為：

(2)

(3)

獨立樣本檢驗的統計指標為：

2 結果與分析

2.1 不同鐵水平下降香黃檀植株生長量變化

整個試驗周期下降香黃檀在樹高、冠幅、地莖生長量的變化，結果如圖2所示。與對照組相比，施加鐵肥不同程度地促進了植株的生長。從CK到F2變化過程中，樹高、冠幅、地莖的增長量是持續上升的。但在F3梯度下，樹高增長速度與F2相比明顯降低，冠幅的增長量則小于F2，但地莖的增長較為明顯，增長量是F2的2倍以上。試驗結果說明適當添加鐵肥有助于提高生物量，但過量后樹高、冠幅生長減緩，地莖增大。

圖2 不同鐵肥梯度下降香黃檀樹高、冠幅、地徑生長量Fig.2 Growth of height,crown width and stem of Dalbergia odorifera with different iron fertilizer gradientsCK:對照組；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的變化趨勢

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

表2為圖2中各點對應的值。B波段光譜反射率在CK～F2梯度持續下降，但在F3梯度出現上升；G波段則相反，即在CK～F2梯度持續上升，在F3梯度下降；R波段反射率則在CK～F1梯度下降，F1～F2梯度上升，F2～F3基本無變化；RE波段則連續上升。這是由于在合理的鐵濃度下，葉片葉綠素含量升高，光合作用加強，對B、R、RE波段的吸收加強，G波段的反射率增強；而當全鐵出現過量時，植物葉綠體遭到破壞，葉綠素含量降低，則出現與之相反的結果。近紅外波段反射率受細胞結構的影響，隨著施鐵量的增加，植株的生長更加旺盛，細胞分裂與分化能力增強，葉片增大，細胞間空隙增大，所以近紅外波段反射率持續增加。

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

2.2 TIC與圖像參數的關系

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

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

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

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

① CK:對照組；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 不同波段下全鐵含量與紋理特征均值的顯著性檢驗①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 不同波段下全鐵含量與紋理特征方差的顯著性檢驗Tab.4 Significance test between total iron content and variance texture factors in different bands

2.3 變量相關性分析及篩選

通過公式及圖4A可以看出(0～30分別對應表1中的植被指數，順序為RVIB、RVIG、RVIR、RVIRE、DVIB…，以此類推)，前15個VIs之間的相關性較大，即變量之間具有較強的線性相關性。而后16個VIs之間的Pearson系數多數集中在0.7以下，變量之間的相關性較小。同理對與TIC在0.01、0.05水平上相關的15個TFMV，以及與TIC在0.01水平上相關的10個TFV分別進行相關性分析，結果如圖4(b)-(c)所示。圖4B中15個TFMV之間的相關性較大，圖像呈現出的顏色較深，大部分區域相關性在0.8以上；而圖4C中對角線附近的顏色較深，即說明部分相鄰參數之間的相關性較大。

圖4 參數相關性分析Fig.4 Correlation analysis of parameters

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

表5 變量篩選結果①Tab.5 Texture parameter selection results

①VIs代表植被指數；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神經網絡模型構建

對每一組試驗構建PSO-BP神經網絡模型，過程如下：

2) 模型優化。在BP神經網絡中內嵌一個種群規模為20，迭代次數為100的PSO算法。為提高算法的全局收斂能力，在PSO算法進化方程中加入慣性權重因子，設置初始慣性權重為0.9，結束慣性權重為0.1，使算法在迭代初期保持較強的全局搜索能力，在迭代后期能進行更精確的局部開發。

2.5 全鐵含量預測結果分析

表6 不同試驗下PSO-BP神經網絡模型預測結果①Tab.6 Prediction results of PSO-BP neural network model in different tests

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

圖5 CA與GA方法實測預測值Fig.5 Measured and predicted value of CA and GA method

2.6 不同特征組合對模型預測結果的影響

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

表7 不同特征組合下的預測結果Tab.7 Prediction results of different feature combination

3 結論

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