基于高光譜成像的油菜苗期溫度脅迫檢測方法

張小敏 張延寧 姜海益 王怡田 林洋洋 饒秀勤

(1.浙江大學生物系統工程與食品科學學院， 杭州 310058； 2.農業農村部農產品產地處理裝備重點實驗室， 杭州 310058；3.浙江大學數學科學學院， 杭州 310058)

0 引言

油菜是我國大宗油料作物之一，長江中下游地區為我國重要的冬油菜生產基地，常年種植面積和產量均占80%以上的份額[1]。在10月中下旬適宜移栽期內，氣溫下降，容易出現持續寡照天氣，低溫脅迫會造成油菜種子發芽時間延長，容易對幼苗造成不可逆的破壞，這是油菜苗不健壯的主要原因之一[2-4]。因此，研究溫度脅迫對油菜苗期生長的影響以及脅迫檢測方法非常必要。

溫度脅迫會誘發油菜苗外部形態、生化組分含量、酶活性[5-7]等發生相應變化。在眾多脅迫檢測方法中，實驗室生化檢測是目前最為客觀準確的方法，但存在采樣復雜、時效性差、專業性強、成本過高等問題，不利于規模化、標準化作業[8]。機器視覺技術能夠快速、穩定地采集高精度的作物圖像信息，可節省人力和降低人為判斷的主觀性[9]。但由于傳統圖像數據大多只能反映作物外觀特征，故構建的脅迫檢測模型也存在局限性[10-12]。高光譜成像技術能夠提供待測作物更為詳細和全面的信息，基于光譜信息既可以獲取待測作物的外形，還可以反映作物生化組分含量和內部結構的變化[13-16]。

文獻[17]采集了小麥感染鐮刀菌400～1 000 nm波長范圍的高光譜圖像，基于主成分分析提取光譜子范圍特征，基于光譜角度匹配提取全譜特征，各染病時間平均檢測準確率為67%，不同感染程度的檢測準確率差異較大。文獻[18]采集了菠菜葉片冷藏過程中400～1 000 nm波長范圍的高光譜圖像，研究了冷藏條件下光譜特性的演化規律。文獻[19]采用400～900 nm波長范圍的高光譜來監測冷藏室和室溫下不同時長的玉米葉片光譜反射率，研究發現，1 h以后冷藏箱內的葉片反射率比正常室溫下的葉片反射率大。文獻[20]采集了番茄葉片早疫病、晚疫病380～1 023 nm波長范圍的高光譜圖像，基于連續投影算法(Successive projections algorithm, SPA)提取特征波長，建立了極限學習機分類器，對病害類別的多重分類準確率達97.1%。文獻[21]采集了玉米葉片干旱脅迫過程中400～1 000 nm波長范圍的高光譜圖像，提出一種基于聚類、監督波長選擇以及光譜相似性度量的檢測方法，研究表明，與傳統圖像相比，高光譜圖像能滿足更廣泛的作物生物脅迫和非生物脅迫的檢測需求[22-23]，通過觀察作物生長期內光譜的演化規律，有助于判斷脅迫類別、揭示脅迫過程、提高脅迫檢測效果，從而分析得出作物不健壯的原因[24-26]。

綜上，高光譜成像技術為脅迫檢測提供了一種快速、準確的檢測手段。本文以油菜為研究對象，以高光譜成像為技術手段，進行油菜苗溫度脅迫實驗。采用SPA和連續小波變換-逐步判別分析法(Continuous wavelet transform-stepwise discriminant analysis, CWT-SDA)在溫度脅迫敏感波段處提取特征波長，并通過持續觀測油菜苗期光譜特征的演化規律，研究分析溫度脅迫對油菜葉片光譜特征的影響，并建立脅迫判別模型，以期實現正常幼苗與脅迫幼苗的準確區分。

1 材料和方法

1.1 實驗材料與實驗設計

實驗時間為2019年12月12日至2020年1月20日，實驗地點為浙江大學紫金港校區生物系統工程與食品科學學院實驗室。以油菜品種浙油50作為實驗材料，播種前曬種4 h。配置基質含水率70%，將基質裝盤后，用壓穴模板在每個穴孔上壓出深1.0～1.5 cm的圓形播種穴，把顆粒飽滿、顏色一致、大小均勻油菜種子播種于穴盤中，1穴播種1粒種子，采用72穴的穴盤共播種18盤。播種完畢后均等置入智能人工氣候箱中培養40 d。

培養條件：晝溫(25±1)℃，夜溫(15±1)℃，相對濕度(70±5)%，光周期白天16 h、黑夜8 h，光密度350 μmol/(m2·s)。

待幼苗長至兩葉一心后，進入溫度脅迫實驗階段。將幼苗的育苗溫度分別設為20℃(T1)、25℃(T2)和30℃(T3)，其中T2組為對照組。分別在第1天(d1)至第21天(d21)，每天從每組幼苗中摘取20株油菜苗葉片，分別采用Q285型高光譜成像系統(德國Cubert公司)和SOC710 SWIR型高光譜成像系統(美國SOC公司)采集450～998 nm和917～1 717 nm的高光譜圖像。以每株油菜苗所有冠層葉片的平均光譜作為冠層光譜，構建冠層光譜樣本集，樣本總數均為1 260個。其中，T1、T2、T3組各420個。

1.2 高光譜數據采集

由于在可見光波段到近紅外波段范圍內，葉片光譜的反射特性與葉片生化組分的含量以及內部結構特性有關，不同波段能反映獨特的理化性質。為了全面地獲取油菜葉片光譜信息，本實驗采用Q285型和SOC710 SWIR型兩臺高光譜成像系統分別采集油菜葉片450～998 nm和917～1 717 nm的高光譜圖像。采集示意圖如圖1所示，葉片放置于黑色背景的載物臺上，采集時關閉暗箱門，單個樣本的總采集時間在1 min內。

其中，Q285型高光譜成像系統的光譜范圍為450～998 nm，光譜分辨率為4 nm，鏡頭焦距23 mm，像素分辨率為1 000像素×1 000像素；SOC710 SWIR型高光譜成像系統的光譜范圍為917～1 717 nm，光譜分辨率為2.75 nm，鏡頭焦距35 mm，像素分辨率為640像素×568像素。

1.3 光譜預處理

1.3.1感興趣區域提取

由于高光譜圖像采集系統采集的是整個視場范圍內的圖像，而感興趣區域是葉片區域，因此需要對葉片和背景進行分割。從樣本集中隨機抽取10個樣本，考慮到采集時背景均勻且不變，分別選取高光譜圖像中背景區域和葉片區域100像素×100像素，并分別求樣本的背景區域和葉片區域的平均光譜，如圖2所示。可以看出10個樣本在同一波段內的光譜曲線相似，圖2a中在550 nm和750 nm附近波段背景和葉片光譜曲線反射率差異較大，而圖2b中1 100 nm附近波段背景和葉片光譜曲線反射率差異較大，故上述波長圖像均可用于背景分割。

由于背景光譜反射率較小時，與葉邊緣區分度較好，最終選擇了546 nm和1 121 nm波長圖像分別對450～998 nm波段和917～1 717 nm波段高光譜圖像進行分割，對隨機抽取樣本中1號樣本的分割結果如圖3所示。由于不同葉片在同一波段內的光譜特性是相似的，故其余樣本都選用上述波長圖像進行圖像二值化操作，提取波段范圍下的感興趣區域。

1.3.2異常樣本剔除

提取冠層葉片的平均光譜后，去除兩端噪聲波段，選擇498～950 nm共114個波長和939～1 681 nm 共267個波長進行分析，并依次進行Savitzky-Golay平滑、標準化和歸一化處理。

為避免植株個體生長差異、儀器誤差、測量環境變化、操作失誤等原因產生的異常樣本影響光譜分析結果，需要檢測異常樣本并根據實際情況予以保留或剔除。采用主成分-馬氏距離法，根據Dixon準則或Chauvenet準則剔除異常樣本。計算T1、T2、T3各組平均光譜與每個樣本光譜的馬氏距離，按馬氏距離從小到大順序排列；指定檢出水平α=0.05，剔除水平α=0.01。首先分析各組光譜數據是否有明顯異常，根據Chauvenet準則檢出與平均光譜偏差在95%置信區間外的光譜。T1組馬氏距離分布情況如圖4所示。檢出異常樣本為長至兩葉一心后第17天的樣本(圖4中虛線所示)，可能是這天的數據錄入存在問題。

隨后分析各天光譜數據，根據Dixon準則檢出與平均光譜偏差在95%置信區間外的光譜。長至兩葉一心后第1天各組馬氏距離分布情況如圖5所示，未檢出異常樣本。類似地，依次檢測其余19 d樣本，均未檢出異常樣本。因此，最終有20 d樣本用于特征提取和健壯苗識別研究，每天60個樣本，T1、T2、T3組各20個。

1.4 溫度脅迫敏感波段提取

對光譜響應進行初步分析發現，在觀測初期，葉片生長對冠層光譜的影響較大，故需要區分生長和溫度脅迫的主要響應波段。按照觀測時間將冠層光譜樣本集劃分為d1～d21族，每族分為T1、T2、T3組。

1.4.1基于光譜反射率的顯著性分析

首先基于光譜反射率分析溫度脅迫敏感波長，指定顯著性水平α=0.05，依次對di族λj波長反射率r(di,λj)進行正態檢驗和方差齊性檢驗；若3組均滿足正態或輕微偏態分布和方差齊性要求，采用單因素方差分析，否則采用Mann-Whitney-Wilcoxon(MWW)檢驗，檢驗r(di,λj)在3組之間是否存在顯著差異；若單因素方差分析結果存在顯著差異，采用Student-Newman-Keuls(SNK)法進行多重比較，檢驗3組之間是否兩兩存在顯著差異；若兩兩存在顯著差異，則將λj作為di族溫度脅迫敏感波長，并統計d1～d21族λj出現的頻數。

保留頻數占總族數百分比不低于75%，即頻數不低于16的波長作為溫度脅迫敏感波長，溫度脅迫敏感波長分布情況如圖6所示。因此，基于光譜反射率共提取2個溫度脅迫敏感波段：554～714 nm、1 559～1 567 nm。

1.4.2連續小波變換分析

連續小波變換(CWT)是一種自適應的時頻分析方法，可以進行多分辨率分析，適合相似性檢測和奇異性分析。其中，一維連續小波分析通過用縮放和平移操作的一系列小波函數與原始光譜信號卷積，對于小波函數ψ(t)，連續小波變換公式為

(1)

(2)

式中ψa,b(λ)——連續小波基函數

λ——時刻

a——尺度因子

b——平移因子

w(a,b)——小波變換系數

x(λ)——連續平方可積函數

a和b可連續取值，a一般取整數，有別于二進制小波變換的取2整數冪。

采用CWT算法將冠層光譜在不同尺度下轉換成一系列的小波系數，圖7展示了分解尺度為1、5、10、15、20的小波系數圖。低尺度的小波系數反映光譜圖上局部峰谷變化，高尺度的小波系數極值發生偏移，反映光譜曲線變化趨勢和相對變化程度。

本文選擇Marr小波母函數，采用1～50連續的小波分解尺度，基于小波系數分析溫度脅迫敏感波段。以d1族為例：

(1) 計算d1族光譜樣本1～50分解尺度下的小波系數w(si,λj)，其中si表示分解尺度，得到小波系數矩陣。采用單因素方差分析，檢驗w(si,λj)在d1族3組之間是否存在顯著差異，得到P值矩陣，如圖8所示，P值越大，灰度越高，差異越不顯著。

(2)若存在顯著差異，采用SNK法進行多重比較，檢驗d1族3組之間是否兩兩存在顯著差異；若兩兩存在顯著差異，則將λj作為d1族溫度脅迫敏感波長，將w(si,λj)作為d1族敏感小波系數。如圖9所示，淺藍色、橙色、黃色區域分別代表3次比較中有1、2、3次在0.05水平上顯著。尤其是可見光波段各分解尺度下的小波系數兩兩比較幾乎都存在顯著差異，說明短期內溫度脅迫主要影響可見光波段。

(3)類似地，依次對其余20族進行分析，統計λj出現的頻數。保留頻數不低于16的波長，作為溫度脅迫敏感波長。基于連續小波變換提取的溫度敏感波長分布情況如圖10所示，此方法提取到的溫度脅迫敏感波長僅在近紅外波段。

基于CWT算法共提取2個溫度脅迫敏感波段：939～967 nm、1 210～1 213 nm。

因此，兩種算法共提取4個溫度脅迫敏感波段：554～714 nm、939～967 nm、1 210～1 213 nm和1 559～1 567 nm。將以上4個溫度脅迫敏感波段，用于光譜特征提取算法研究。

1.5 光譜特征提取

1.5.1樣本集劃分

油菜苗每隔一定時間生長出一片新葉，觀測期間，油菜生長經歷了從二葉期至五葉期的不同生長期，不同生長期具體的溫度敏感特征可能不同。根據生長期將冠層光譜樣本集劃分為C1、C2、C3、C4族，分別對應二葉期至五葉期。樣本集具體劃分情況如表1所示。

表1 冠層光譜樣本集劃分

1.5.2基于SPA算法的特征波長提取

SPA算法是在光譜中應用廣泛的特征選擇算法，通過迭代選取共線性低的波長組合來降低特征維度。首先將樣本分為用于算法迭代過程和確定最優變量的兩個子集；計算初始波長在其余波長上的投影，加入投影向量最大的波長，重復以上步驟直到迭代結束，建立基于波長組合的預測模型，采用均方根誤差判斷模型優劣。

采用SPA進行特征波長提取，以C1族為實例，指定特征波長數為1～20，根據均方根誤差篩選最優特征波長。優選過程中，均方根誤差隨特征波長數變化趨勢如圖11所示，當波長數大于7時，均方根誤差隨波長數的增加變化較小且穩定在較低值。SPA優選特征波長如圖12所示，篩選出622、642、942、953、959、962、1 567 nm共7個特征波長。SPA一定程度上消除了各波長變量之間的共線性影響，提取的特征代表性較強。

類似地，提取C2～C4族特征波長，具體結果如表2所示。

表2 C1～C4特征波長

1.5.3基于CWT-SDA算法的小波特征提取

以C1族為例，進行CWT處理，算法見1.4.2節。考慮到低尺度的小波系數容易受到噪聲的影響，最終選擇尺度不低于5、波長分布在溫度敏感波段內的小波系數。同時為避免特征個數較多和特征間相關性較強，導致矩陣求逆的精度下降，所建判別函數不穩定的問題，本文采用SDA算法篩選小波特征，尋找合適的特征子集。

SDA算法是一種迭代篩選變量的方法，通過逐個加入變量，達到篩選有效變量的目的，同時得到判別系數和常數。逐步判別分析對判別變量附加信息的檢驗采用Λ統計量，公式為

(3)

式中W——組內離差矩陣

T——總離差矩陣

(4)

(5)

式中nα——各組樣本總數

類似地，提取C2～C4族的小波特征，具體結果如表3所示。

表3 C1～C4小波特征

1.6 建模分析方法

1.6.1Fisher線性判別

本文采用Fisher判別器作為分類模型進行三分類，即分為T1、T2、T3組。Fisher分類器是一種結構簡單的線性有監督分類模型，以類內散度最小而類間散度最大作為準則進行投影降維，其對噪聲有著良好的魯棒性。對于本文所用到的Fisher線性判別，樣本類內散度矩陣、樣本類間散度矩陣和最大化目標分別定義為

(6)

(7)

(8)

式中Ci、μi——第i類的樣本集合、均值

mi——第i類的樣本個數

μ——這3個類全部樣本的均值

ω——需要求解的投影超平面

可以通過特征值分解，取前d個對應的特征向量求解，從而將N維樣本映射到d維進行分類。

1.6.2模型評價

本文采用各樣本集(C1、C2、C3、C4族)總體分類準確率對識別模型進行評價，以C1族為例，將樣本劃分成5個相等的子集(每個子集72個樣本)，每次取4個子集作為訓練集，余下的子集作為測試集，采用無放回隨機抽樣的5折交叉驗證法進行模型訓練和測試，5次分類準確率的平均值作為模型總體分類準確率。其中，在識別健壯苗時，T2組(25℃)作為正常幼苗樣本，T1組(20℃)和T3組(30℃)分別作為輕度低溫脅迫和輕度高溫脅迫幼苗樣本。

2 結果與分析

2.1 分類基線模型

分別對基于SPA算法提取的特征波長和CWT-SDA算法提取的小波特征，使用Fisher判別模型，構建基線模型，初步評價特征波長區分輕度低溫、輕度高溫和正常溫度下生長的油菜苗的能力。

兩種模型的分類準確率如表4、5所示。

表4 SPA算法優選波長在Fisher判別模型下的溫度脅迫檢測結果

表5 CWT-SDA優選小波特征在Fisher判別模型下的溫度脅迫檢測結果

分析兩個模型的分類準確率表明：SPA-Fisher模型和CWT-SDA-Fisher模型的各生長期平均總體分類準確率分別為85.31%和86.64%，后者的模型準確率和穩定性略優；T1組和T3組分類準確率隨脅迫時間增長提高，說明溫度脅迫造成的累積光譜響應差異逐漸明顯；T1組分類準確率高于T3組，說明低溫脅迫響應更易被區分，尤其是C1族T1組的分類準確率顯著高于T3組，葉片生長初期對低溫更加敏感。

2.2 分類模型優化

2.2.1波段特征提取方法

由于SPA和CWT-SDA算法提取的光譜特征在4個溫度敏感波段均有分布(圖6和圖10)，其中554～714 nm存在多個波長相鄰分布情況。相鄰波長往往具有較強的相關性，易導致特征冗余。本文通過進一步提取波段特征，減少特征冗余，進而提升模型穩定性。

如圖13所示，長至兩葉一心后以不同溫度處理2 d后的554～714 nm平均光譜曲線為例，不同溫度處理組的冠層光譜554～714 nm光譜曲線形狀相似，表現出綠色作物特有的綠峰-紅谷-紅邊特征，但不同溫度處理下綠峰、紅邊藍移、紅移速度不同，綠峰峰值下降的程度也不同，導致了680 nm附近紅谷曲率不同。因此，對554～680 nm分段提取波段MA曲線面積和波段CN正切特征值tanθ兩種特征，以更好地描述各組光譜曲線的差異。

(1)提取554～714 nm波段范圍內小波系數曲線的過零點。如圖14b局部圖所示，分別設置：A，582 nm；B，630 nm；C，646 nm。并設置：M，554 nm；N，714 nm。進一步提取MA和CN波段特征。

(2)采用3階多項式擬合CN。如圖15所示，為避免邊界點噪聲，各向內取一個點，計算C′、N′處切線，C′、N′與切線交點O構成的夾角∠C′ON′為CN的特征角，∠C′ON′角度的正切值為CN的波段特征tanθ。

(3)采用3階多項式擬合MA，計算MA曲線面積，作為MA的波段特征。

2.2.2光譜特征演化規律

以3 d為間隔，將觀測日期劃分為D1～D7，則D1～D5為育苗期，D6、D7為移栽期。對每一個特征進行測試，衡量該特征和分類值之間的關系，綜合考慮顯著性分析結果與單特征分類結果，篩選出以下特征，進一步分析特征均值演化規律。

(1)554～714 nm波段特征

MA曲線面積變化趨勢如圖16a所示。從圖中可以看出，隨觀測期推移，曲線面積先下降后穩定，穩定后從大到小依次為T3、T1、T2。tanθ變化趨勢如圖16b所示。從圖中可以看出，隨觀測期推移，tanθ先上升后穩定，穩定后從大到小依次為T1、T2、T3，與紅谷附近曲率有關，反映了綠峰、紅邊藍移、紅移速度和綠峰峰值下降程度的差異。

(2)特征波長

1 213 nm光譜反射率變化趨勢如圖17a所示。從圖中可以看出：隨觀測期推移，1 213 nm反射率先下降后穩定，穩定后從大到小依次為T2、T3、T1。1 567 nm光譜反射率變化趨勢如圖17b所示。從圖中可以看出：隨觀測期推移，1 567 nm反射率呈現先下降后穩定，穩定后從大到小依次為T3、T2、T1。

(3)小波特征

w(9, 967)變化趨勢如圖18a所示，不同觀測期內，w(9, 967)分布區間在3組間不重疊，移栽期時從大到小依次為T2、T1、T3。w(13, 1 213)變化趨勢如圖18b所示，隨觀測期推移，w(13, 1 213)先下降后上升，移栽期時從大到小依次為T2、T3、T1。w(7, 1 567)變化趨勢如圖18c所示，隨觀測期推移，w(7, 1 567)先下降后穩定，移栽期時T1組明顯小于T2、T3組。

2.3 受溫度脅迫的健壯苗識別模型

基于以上7個特征，即554～714 nm波段MA曲線面積、tanθ，1 213 nm和1 567 nm處反射率，以及小波特征w(9, 967)、w(13, 1 213)和w(7, 1 567)，建立多特征融合的Fisher判別模型，采用無放回隨機抽樣的5折交叉驗證方法進行模型建立和驗證，5次分類準確率的平均值作為模型總體分類準確率，如表6所示。

表6 多特征融合在Fisher判別模型下的溫度脅迫檢測結果

結果表明：多特征融合的Fisher判別模型平均分類準確率為88.68%，優于SPA-Fisher和CWT-SDA-Fisher模型，在生長期C2(即三葉期)時達到最佳檢測準確率，為95.56%；多特征融合的Fisher判別模型顯著提高了各生長期的T1組分類準確率，各生長期均達到90%以上。同時發現554～714 nm波段特征和特征波長在育苗期的演化曲線呈現一定的規律性，曲線特征對于脅迫檢測的效果有一定提升。

3 結束語

利用高光譜成像技術對油菜幼苗進行溫度脅迫檢測，篩選出554～714 nm波段MA曲線面積、正切特征值tanθ、1 213 nm和1 567 nm處反射率以及小波特征w(9,967)、w(13,1 213)、w(7,1 567)共7個特征，建立了多特征融合的溫度脅迫Fisher判別模型。結果表明，該模型平均分類準確率為88.68%，優于SPA-Fisher和CWT-SDA-Fisherm模型的分類效果，多特征融合的Fisher判別模型顯著提高了各生長期的T1組分類準確率。同時發現，554～714 nm波段特征和特征波長在育苗期的演化曲線呈現一定的規律性，且低溫脅迫檢測效果較好、脅迫響應更易被區分。綜合考慮各模型不同生長期的檢測效果，發現三葉期是較適合進行溫度脅迫檢測的時期。