基于太赫茲成像技術(shù)的大豆葉片水分含量測(cè)定

步正延，李臻峰，宋飛虎，李 斌，李 靜

(1.江南大學(xué) 機(jī)械工程學(xué)院，江蘇 無錫 214000； 2.江蘇省食品先進(jìn)制造裝備技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 無錫 214122； 3.國(guó)家農(nóng)業(yè)信息化工程技術(shù)研究中心，北京 100097)

葉片作為作物自身最重要的生理器官之一，其內(nèi)部的水分信息既可以直接反映出整個(gè)作物的水分情況，又可以間接了解到作物生理生化狀況和生長(zhǎng)發(fā)育過程[1-2]。農(nóng)學(xué)中傳統(tǒng)的葉片水分測(cè)量方法主要有烘干法、化學(xué)法、電容法等，但這些方法都具有破壞性，且非連續(xù)，又很浪費(fèi)時(shí)間[3-4]。如何快速、精準(zhǔn)地獲取作物葉片的水分信息是一個(gè)亟待解決的科研問題。

太赫茲(terahertz，THz)一般是指頻率范圍為0.1～10 THz(1 THz=1012Hz)的一段電磁波，位于微波和遠(yuǎn)紅外之間，是光子學(xué)向電子學(xué)的過渡區(qū)域，在電磁波波譜中占有很特殊的位置[5-6]。水分子在平衡位置附近處平動(dòng)和轉(zhuǎn)動(dòng)的弛豫時(shí)間處于皮秒、亞皮秒量級(jí)，在與太赫茲波相互作用時(shí)，會(huì)進(jìn)行選擇性吸收[7-9]。因此，在太赫茲波段研究水的特性譜，對(duì)測(cè)量葉片水分含量、研究葉片內(nèi)部組織特性等問題，具有重要意義。

目前，國(guó)內(nèi)外已有許多學(xué)者基于太赫茲技術(shù)開展了葉片檢測(cè)的探索研究[10-11]。左劍等[12]利用太赫茲時(shí)域光譜技術(shù)分別對(duì)全綠和泛黃的葉片進(jìn)行測(cè)量，發(fā)現(xiàn)太赫茲光譜在葉片不同位置處的吸收率和折射率與葉片中的水分分布呈現(xiàn)出某種特殊的規(guī)律。Gente等[13]建立了一種基于太赫茲時(shí)域光譜技術(shù)的葉片含水量預(yù)測(cè)模型，用這種方法測(cè)量得到的結(jié)果與葉片真實(shí)水分含量具有很好的一致性。不難看出，用太赫茲技術(shù)來檢測(cè)葉片含水量是可行的。然而，已有的研究主要都是圍繞著太赫茲光譜技術(shù)來展開葉片水分測(cè)量分析的，目前尚缺少?gòu)奶掌澇上窠嵌葘?duì)葉片進(jìn)行水分檢測(cè)的研究。

本文基于太赫茲成像技術(shù)，以大豆葉片作為研究對(duì)象，采用主成分分析法挑選出水分敏感波段，對(duì)特征頻率下的圖像進(jìn)行特征參數(shù)提取和分析，并結(jié)合不同算法建立水分預(yù)測(cè)模型，對(duì)比分析各模型性能，挑選出一種基于最優(yōu)參數(shù)的葉片含水率預(yù)測(cè)模型，為作物葉片水分含量檢測(cè)提供新方法和新思路。

1 材料與方法

1.1 試驗(yàn)樣本采集與稱量

本試驗(yàn)所用大豆葉片樣本均采摘于北京農(nóng)林科學(xué)院小湯山試驗(yàn)基地溫室，大豆品種為中黃13。選取32盆長(zhǎng)勢(shì)良好的開花期大豆植株，每一盆分別剪取上層、中層、下層各一片葉片，總計(jì)96片葉片。在剪取過程中，盡可能保證葉片大小、厚度相似、形狀完好。剪取完后，立即用密封袋將葉片封裝，裝入保濕箱，帶回實(shí)驗(yàn)室。測(cè)試環(huán)境室溫為22.5～24.5 ℃，相對(duì)濕度為48%～52%，測(cè)試過程在國(guó)家農(nóng)業(yè)信息化工程技術(shù)研究中心太赫茲實(shí)驗(yàn)室完成。

試驗(yàn)前，先對(duì)葉片進(jìn)行干燥處理，以得到不同含水量的葉片樣本。葉片在干燥過程中，若失水速度較快，會(huì)發(fā)生蜷縮現(xiàn)象，使得葉片表面不平整、出現(xiàn)褶皺，導(dǎo)致葉片圖像不易采集。為解決這一問題，本試驗(yàn)采用自然干燥處理，將葉片背面朝上置于室溫環(huán)境下，通過控制時(shí)間變量得到不同含水量的葉片。

將干燥后的葉片置于精度0.1 mg的電子天平(上海菁海儀器有限公司，型號(hào)為FA2004N)上稱量，重復(fù)3次，取平均值，記為m1。葉片稱量完成后，以高聚乙烯板(high polyethylene sheet，HDPE)作為背景板，將固定有葉片的背景板夾持在二維掃描平移臺(tái)上，等待系統(tǒng)掃描成像。

1.2 試驗(yàn)系統(tǒng)與太赫茲圖像獲取

試驗(yàn)所用的太赫茲光譜成像系統(tǒng)主要包括透射式太赫茲時(shí)域光譜系統(tǒng)(terahertz time-domain spectroscopy system，THz-TDS，德國(guó)門羅系統(tǒng)公司，型號(hào)為TERA K15)、二維掃描平移臺(tái)(德國(guó)門羅系統(tǒng)公司，型號(hào)為TERA Image)和PC機(jī)(日本富士通公司，型號(hào)為USD-D3164)。其中，透射式THz-TDS主要包含以下4部分：太赫茲源、太赫茲發(fā)射器、太赫茲探測(cè)器，以及相關(guān)的光路系統(tǒng)。數(shù)據(jù)采集與保存均在LabVIEW軟件中完成。

如圖1所示，由鈦寶石飛秒激光器發(fā)射出飛秒激光，經(jīng)分束鏡分為2束。一束為抽運(yùn)光，經(jīng)過可變延遲線，入射到光電導(dǎo)天線裝置砷化鎵(GaAs)晶體上，產(chǎn)生THz波脈沖，然后經(jīng)過2組離軸拋物面鏡，最終將攜帶葉片樣品信息的THz波聚焦到光電檢測(cè)轉(zhuǎn)置碲化鋅(ZnTe)晶體探測(cè)器上。另一束為探測(cè)光，經(jīng)多次反射后，由硅片將其反射到探測(cè)晶體上，使THz波與探測(cè)光最終共聚于探測(cè)晶體之上。通過電光取樣探測(cè)出太赫茲脈沖的整個(gè)時(shí)域波形信號(hào)，經(jīng)鎖相放大器傳至PC機(jī)中。

采集一次數(shù)據(jù)完畢后，再通過PC機(jī)控制二維平移臺(tái)移動(dòng)，采用雙程方式掃描整片葉片(100 pixel×108 pixel)，記錄數(shù)據(jù)。每片葉片的掃描時(shí)間約為20 min。試驗(yàn)中為消除高聚乙烯板作為背景板對(duì)測(cè)試結(jié)果的影響，根據(jù)式(1)進(jìn)行太赫茲光譜圖像的校正。

式(1)中：Rimg為原始的太赫茲光譜圖像，Rbgb為掃描標(biāo)準(zhǔn)高聚乙烯板得到的標(biāo)定太赫茲光譜圖像，R表示校正后的太赫茲光譜圖像。

在測(cè)試過程中，為避免空氣中水分等極性分子對(duì)太赫茲波的強(qiáng)烈吸收，將試驗(yàn)系統(tǒng)裝置密封，封，充滿氮?dú)?，保證測(cè)試環(huán)境相對(duì)濕度小于4%、溫度在22 ℃左右。

圖1 透射式THz-TDS系統(tǒng)裝置圖Fig.1 Diagram of transmission-type THz-TDS system

1.3 試驗(yàn)樣本含水率的測(cè)定

數(shù)據(jù)采集完成后，采用烘干法測(cè)定葉片含水率。將葉片放于干燥箱(廣州格蘭仕公司，型號(hào)為KWS1525X-F2R)中，設(shè)置烘干溫度為100 ℃，分別烘干1 h 和1.5 h后，記錄2次烘干后的質(zhì)量。若2次稱量值之差小于0.5 mg，則認(rèn)為葉片已烘干至恒重，記錄此時(shí)葉片質(zhì)量(即烘干1.5 h后的干重)為葉片干重，記為m2。根據(jù)式(2)計(jì)算葉片含水率(w)。

1.4 主成分分析

主成分分析(principal component analysis， PCA)是一種數(shù)據(jù)降維思想，它用少數(shù)的幾個(gè)綜合變量替代原來的大部分變量，排除眾多變量中互相重疊的信息，以達(dá)到數(shù)據(jù)降維的目的。將葉片太赫茲圖像看成是一個(gè)擁有多個(gè)層面、按一定順序疊合而成的圖像集，采用PCA對(duì)其進(jìn)行降維處理，壓縮圖像信息，從太赫茲原始波段中尋找出最能表征葉片水分信息的特征波段，去除其他冗余波段，簡(jiǎn)化分析模型。

1.5 回歸建模

為了挑選最優(yōu)建模參數(shù)，本試驗(yàn)嘗試對(duì)特征波段下的葉片圖像進(jìn)行分離，將其分為葉脈圖像與葉肉圖像，與葉片圖像一道求取各圖像的灰度均值。將各圖像灰度均值參數(shù)作為輸入、葉片含水率作為輸出，分別代入以下回歸模型：多元線性回歸(multiple linear regression， MLR)，反向傳播神經(jīng)網(wǎng)絡(luò)(back propagation-artificial neural network， BP-ANN)和最小二乘支持向量機(jī)(least squares-support vector machine， LS-SVM)，對(duì)各模型進(jìn)行對(duì)比分析，最終選取最優(yōu)模型。

圖像處理、數(shù)據(jù)處理，以及回歸建模與分析均在Matlab 2014a中完成。

2 結(jié)果與分析

2.1 大豆葉片太赫茲光譜曲線分析

由于葉片內(nèi)部不同部位組織結(jié)構(gòu)、表面紋理、厚度等因素的差異，太赫茲光譜在同一葉片不同部位所測(cè)量到的吸收系數(shù)與折射率差異較大[12]。試驗(yàn)選擇圖2所示的葉脈區(qū)域、葉肉區(qū)域和背景板區(qū)域這3個(gè)位置處5 pixel × 5 pixel的區(qū)域，計(jì)算其平均光譜，得到葉片的時(shí)域譜曲線和頻域譜曲線(圖3)。其中，背景板區(qū)域的光譜信息代表參考信號(hào)，葉片區(qū)域的光譜信息代表樣品信號(hào)。從圖3-A中可以看出，樣品信號(hào)的振幅強(qiáng)度要弱于參考信號(hào)，且在相位上也有延遲，而就樣品信號(hào)自身而言，葉脈信號(hào)相對(duì)于葉肉信號(hào)在幅值上表現(xiàn)出一定程度的衰減，在時(shí)間上也略滯后。這說明太赫茲脈沖信號(hào)透射葉片樣本時(shí)，葉片對(duì)其存在吸收和散射現(xiàn)象，且葉脈區(qū)域的吸收強(qiáng)度要明顯高于葉肉區(qū)域。

圖2 大豆葉片圖像Fig.2 Image of soybean leaf

在頻域譜中，隨著波段的變化，樣品在不同頻率下對(duì)太赫茲波的吸收迥異，即當(dāng)頻率發(fā)生變化時(shí)，葉片太赫茲圖像的成像對(duì)比度也發(fā)生變化，從視覺角度來分析，選用該種模式更利于后續(xù)的成像研究[14]。如圖3-B所示，在0.2～1.6 THz波段，3條譜線的幅值情況基本與圖3-A類似：參考信號(hào)>葉肉>葉脈。說明在這一波段范圍內(nèi)，三者間具有非常好的成像對(duì)比度，能較好地區(qū)分出葉脈與葉肉部分，且樣品區(qū)域與參考區(qū)域也具有較好的區(qū)分度，而1.6 THz以后的波段信號(hào)強(qiáng)度較弱，信噪比也較低，成像對(duì)比度較差；因此，只截取0.2～1.6 THz，共149個(gè)波段，作為光譜有效波段作進(jìn)一步研究。

2.2 太赫茲特征圖像的選取

針對(duì)每個(gè)大豆葉片樣本的太赫茲圖像，本試驗(yàn)選取了圖2所示的太赫茲圖像中的上、下、左、右4個(gè)10 pixel × 10 pixel的矩形感興趣區(qū)域的光譜平均值作為樣本原始光譜數(shù)據(jù)(96維×149維)。采用高光譜圖像分析中較為常見的分析方法——PCA[2]，對(duì)太赫茲原始光譜數(shù)據(jù)進(jìn)行降維處理，根據(jù)方差貢獻(xiàn)率來挑選出葉片樣本的太赫茲特征波段，用于后續(xù)試驗(yàn)處理。

圖3 太赫茲時(shí)域(A)、頻域(B)光譜圖Fig.3 THz time domain (A) and frequency domain (B) spectrogram

經(jīng)主成分分析，得到如圖4所示的96個(gè)樣本的前3幅主成分圖像PC1、PC2、PC3(累積貢獻(xiàn)率大于99%)下各波段的平均權(quán)重系數(shù)，可以看出，雖然PC1包含的信息量最大，但無明顯的特征峰，因此，綜合PC2和PC3，挑選主成分載荷絕對(duì)值大于0.1的這些點(diǎn)所對(duì)應(yīng)的頻率作為水分敏感特征波段[15]，分別為0.557、1.098、1.163、1.411 THz。

圖4 主成分載荷曲線Fig.4 Principal component load curve

A， 0.557 THz; B， 1.098 THz; C， 1.163 THz; D， 1.411 THz.圖5 特征頻率下的大豆葉片太赫茲圖像Fig.5 THz images of soybean leaf at characteristic frequencies

圖5為上述4個(gè)特征頻率下的大豆葉片太赫茲圖像。前3個(gè)頻率下的太赫茲圖像均具有較好的成像效果，圖中脈絡(luò)分布清晰，且能夠較為直觀地看出葉片中水分分布的情況，而圖5-D中的葉片圖像與前三者相比，成像效果較差，不僅整副圖像灰度偏低，較難區(qū)分出葉脈和葉肉區(qū)域，而且圖像中還出現(xiàn)了較多噪聲點(diǎn)，嚴(yán)重降低了圖像的質(zhì)量。因此，本試驗(yàn)剔除了1.411 THz，最終選用0.557、1.098、1.163 THz這3個(gè)頻率作為葉片水分的敏感特征波段。

2.3 特征圖像處理

以0.557 THz頻率下的太赫茲圖像為例，進(jìn)行處理與分析。在圖像采集過程中，由于系統(tǒng)自身等問題，會(huì)產(chǎn)生條紋偏差，使得圖像質(zhì)量下降，因此在對(duì)圖像進(jìn)行處理前，須先作預(yù)處理。本試驗(yàn)采用非局部均值濾波法對(duì)圖像進(jìn)行濾波處理，圖6-A所示為經(jīng)過濾波后的葉片圖像。采用二維最大信息炯算法對(duì)其進(jìn)行二值化加形態(tài)學(xué)腐蝕和灰度反轉(zhuǎn)處理，得到如圖6-B所示的葉片二值圖像。將二者進(jìn)行逐像點(diǎn)相乘運(yùn)算處理，得到如圖6-C所示的葉片目標(biāo)圖像。為了將圖中的葉脈與葉肉區(qū)域分離，對(duì)圖6-C作自適應(yīng)閾值分割及圖像灰度反轉(zhuǎn)處理，得到葉脈二值圖像(圖6-D)。再將圖6-B與圖6-D做逐像點(diǎn)相減運(yùn)算處理，得到葉肉二值圖像(圖6-E)。最后，分別對(duì)圖6-D和圖6-E做與圖6-C的逐像點(diǎn)相乘運(yùn)算處理，得到葉脈目標(biāo)圖像(圖6-F)和葉肉目標(biāo)圖像(圖6-G)。

2.4 特征參數(shù)提取

當(dāng)太赫茲脈沖信號(hào)透射葉脈與葉肉時(shí)，由于葉脈與葉肉自身含水量的差異，兩者在幅值上存在著較為明顯的區(qū)別(圖3)。因此，本試驗(yàn)考慮對(duì)葉片圖像中的葉脈區(qū)域與葉肉區(qū)域進(jìn)行分離，分別就0.557、1.098、1.163 THz這3個(gè)特征波段下的葉片目標(biāo)圖像、葉脈目標(biāo)圖像，葉肉目標(biāo)圖像，求取出各自的圖像灰度均值，分別記為l1、l2、l3，v1、v2、v3，m1、m2、m3(l代表葉片，v代表葉脈，m代表葉肉，下標(biāo)中的1、2、3分別指代0.557、1.098、1.163 THz特征頻率)。將上述求取的9個(gè)灰度均值特征參數(shù)分成3組：葉片特征模型組(l1、l2、l3)，記為G1；葉脈特征模型組(v1、v2、v3)，記為G2；葉肉特征模型組(m1、m2、m3)，記為G3。對(duì)各組分別進(jìn)行建模分析，尋找葉片最佳水分預(yù)測(cè)模型。

2.5 回歸模型建立與結(jié)果分析

將96份大豆葉片樣本，根據(jù)含水率差異，分成5組：0～20%、20%～40%、40%～60%、60%～80%、80%～100%。校正集由每個(gè)梯度類別中隨機(jī)選取3/4的樣本組成，共計(jì)72個(gè)，預(yù)測(cè)集則由每個(gè)類別中剩下的1/4樣本組成，共計(jì)24個(gè)。將G1、G2和G3中的各特征參數(shù)作為回歸模型的自變量、葉片含水率預(yù)測(cè)值作為因變量，分別建立基于MLR、BP-ANN、LS-SVM的水分預(yù)測(cè)模型。具體模型處理與分析如下。

采用Matlab軟件中的regress函數(shù)，對(duì)原始數(shù)據(jù)進(jìn)行處理分析，得到G1、G2和G3的MLR回歸方程：

A～G分別為濾波后圖像、葉片二值圖像、葉片目標(biāo)圖像、葉脈二值圖像、葉肉二值圖像、葉脈目標(biāo)圖像和葉肉目標(biāo)圖像。A-G were filtered image， leaf binary image， leaf target image， vein binary image， mesophyll binary image， vein target image and mesophyll target image， respectively.

y1=110.2948-0.1332l1-0.0486l2-0.2114l3;

(3)

y2=95.9628-0.1138v1-0.0237v2-0.2919v3;

(4)

y3=114.3278-0.1575m1-0.0206m2-0.2117m3.

(5)

式(3)～(5)中y1、y2、y3分別為G1、G2、G3所對(duì)應(yīng)的葉片含水率預(yù)測(cè)值。

采用Matlab軟件中的BP神經(jīng)網(wǎng)絡(luò)分類器，構(gòu)建3層BP神經(jīng)網(wǎng)絡(luò)，設(shè)置如下：網(wǎng)絡(luò)輸入層結(jié)點(diǎn)數(shù)為3，隱含層結(jié)點(diǎn)數(shù)為4，輸出層節(jié)點(diǎn)數(shù)為1，訓(xùn)練步長(zhǎng)為0.01，動(dòng)量因子為0.90，目標(biāo)誤差為0.001，訓(xùn)練次數(shù)為5 000，隱含層的傳遞函數(shù)為sigmoid型函數(shù)。BP-ANN模型中3組變量的實(shí)際訓(xùn)練次數(shù)分別達(dá)到3 232、3 765和3 523時(shí)，便已達(dá)到指定目標(biāo)誤差，完成收斂。

采用Matlab軟件LS-SVM工具箱中的tunelssvm、trainlssvm和simlssvm等函數(shù)，依次進(jìn)行參數(shù)尋優(yōu)、正則化定參、學(xué)習(xí)訓(xùn)練，以及模型預(yù)測(cè)，試驗(yàn)中，LS-SVM模型采用默認(rèn)的RBF核函數(shù)，當(dāng)3組輸入變量的gamma和sig2參數(shù)分別設(shè)置為8.236 5、2.013 8，8.985 6、12.432 1，35.028 8、6.509 8時(shí)，校正集模型決定系數(shù)(R2)最高，分別為0.955 6、0.857 6和0.967 8。

表1所示為上述9種預(yù)測(cè)模型性能的比較。從圖像特征模型組別來看，G3預(yù)測(cè)模型效果普遍較好，G1次之，G2效果較差，即基于葉肉特征模型組的水分預(yù)測(cè)效果要明顯優(yōu)于葉片組與葉脈脈組。這是因?yàn)楫?dāng)葉片處于失水狀態(tài)下，葉肉組織相比葉脈組織水分流失較快，導(dǎo)致在太赫茲圖像中葉脈部分的灰度要明顯低于葉肉部分，且隨著時(shí)間的變化，這種灰度的差異性也在發(fā)生變化。直接采用葉片圖像建模的話，會(huì)造成圖像中灰度變化的不一致性，將葉脈與葉肉兩者分開建模就避免了這種差異性。此外，由于葉片中80%～90%的區(qū)域均是葉肉部分，葉脈只占很少的一部分，因此葉肉圖像相比于葉脈圖像更能表征整片葉片的信息。綜合這兩方面，基于G3的預(yù)測(cè)模型可以較好地預(yù)測(cè)葉片含水率。

表1 不同模型的葉片含水率建模結(jié)果比較Table 1 Comparison of modeling results of leaf moisture content based on different models

分析G3模型組下3種算法的預(yù)測(cè)效果(表1)，LS-SVM模型的預(yù)測(cè)集的R2值最高，達(dá)到了0.963 2，對(duì)應(yīng)的均方根誤差(root mean square error， RMSE)值最低，為0.046 5，且該模型校正集的結(jié)果與預(yù)測(cè)集大體一致。由此可知， LS-SVM模型基本符合模型評(píng)判中R2高、RMSE低的特點(diǎn)，預(yù)測(cè)效果較好。相較于前者，另外2組模型的預(yù)測(cè)結(jié)果不僅R2值有所下降，且RMSE值也略有提高，模型預(yù)測(cè)效果一般。究其原因，可能是因?yàn)長(zhǎng)S-SVM算法通過最小二乘回歸方法，消除了自變量間多重線性關(guān)系的影響，簡(jiǎn)化了模型結(jié)構(gòu)，尤其是在小樣本回歸分析中，LS-SVM不僅具有良好的泛化能力，同時(shí)也避免了ANN中的局部最優(yōu)問題。因此，本試驗(yàn)選用基于G3的LS-SVM模型作為最佳模型用于預(yù)測(cè)葉片中的含水率，圖7為該模型的校正與預(yù)測(cè)結(jié)果。

圖7 基于G3的LS-SVM模型的校正與預(yù)測(cè)結(jié)果Fig.7 Calibrated and predicted results of LS-SVM model based on G3

3 結(jié)論與討論

利用大豆葉片的太赫茲圖像信息，通過PCA提取出光譜特征波段，處理特征波段圖像，求取圖像特征并進(jìn)行分組，分別通過MLR、BP-ANN和LS-SVM這3種算法進(jìn)行建模。對(duì)比分析各模型性能發(fā)現(xiàn)：1)葉肉特征模型組建模效果最好，葉片特征模型組次之，葉脈特征模型組效果最差：2)從模型決定系數(shù)和均方根誤差來看，基于G3的LS-SVM模型預(yù)測(cè)精度高且穩(wěn)定性好，明顯優(yōu)于其他模型：3)綜合各模型的預(yù)測(cè)結(jié)果可以看出，用太赫茲成像技術(shù)這一手段來檢測(cè)葉片中的含水率是可行的，且具有較好的預(yù)測(cè)效果。本研究同樣適用于其他作物葉片水分含量的測(cè)定，為作物葉片水分檢測(cè)提供了一種新的思路。