王傳宇，郭新宇，杜建軍

（1. 北京農(nóng)業(yè)信息技術(shù)研究中心，北京 100097；2. 數(shù)字植物北京市重點實驗室，北京 100097）

0　引　言

葉面積指數(shù)（leaf area index，LAI）影響作物光合輻射截獲量，水分的吸收，潛熱和感熱通量，地表生態(tài)系統(tǒng)與大氣的 CO2交換。在一些作物生長模型中，LAI也是重要輸入?yún)?shù)和條件變量。LAI是冠層輻射傳輸模型的決定因素，在遙感數(shù)據(jù)同化中，已被作為冠層反射與作物生長模型的關(guān)鍵聯(lián)系[1]。測量作物整個生育時期的LAI指數(shù)時間變化規(guī)律，有助于作物產(chǎn)量、生物量預(yù)測模型的校正。

LAI獲取方法主要分為直接法和間接法2類，所謂直接法就是對直接組成部分進行量測，例如通過測量采樣區(qū)域植株的葉片面積計算 LAI[2-3]，直接法需要大范圍地破壞性采樣，消耗大量人力物力，測量指標主觀依賴性強，因此實際中間接方法應(yīng)用的更為廣泛深入[4-6]，Baret等[7]指出，冠層內(nèi)的輻射傳播主要決定于LAI和葉傾角分布函數(shù)，基于此，發(fā)展出了由冠層孔隙度（gap fraction）計算LAI的方式。冠層空隙度的獲取方法分為2種：1）測量冠層下方多個角度的太陽直接輻射值，如商業(yè)化冠層分析儀AccuPAR（Decagon Inc）；2）采用圖像同時獲得植株與天空/地面背景的像素投影，計算植株背景像素占整幅圖像的比例。其中半球圖像（hemispherical photography）能夠記錄冠層內(nèi)部多個角度的隙度，被眾多學(xué)者研究和應(yīng)用[8-10]。半球圖像的投影模型較為復(fù)雜[11]，對于處于營養(yǎng)生長前期的禾本科作物玉米，測量誤差較大，不能適應(yīng)整個生育時期玉米冠層LAI動態(tài)變化監(jiān)測的需求。有研究表明，可以采用冠層頂部單角度圖像獲得冠層孔隙度進而計算 LAI[12]，簡化了圖像投影模型算法。開展整個生育時期內(nèi)玉米冠層頂部圖像LAI動態(tài)監(jiān)測，計算精度主要受限于圖像分割算法，玉米植株自身形態(tài)及外部光照環(huán)境的顯著變化是影響圖像分割和冠層孔隙度計算精度的主要原因。針對上述研究中存在的不足，本文將研究一種單角度數(shù)字圖像LAI連續(xù)獲取方法，為解決田間變化光照對綠色植株分割的影響，采用紅外圖像去除植株葉片高光反射和陰影，基于高斯分布模型自動計算圖像分割閾值，推導(dǎo)了冠層頂視圖像孔隙度計算LAI的方法，同時建立了植株圖像像素與植株干質(zhì)量和鮮質(zhì)量的模型，并在多個圖像序列上對算法效果進行分析驗證，以期為田間環(huán)境下冠層參數(shù)的自動連續(xù)監(jiān)測提供解決方案。

1　材料與方法

1.1　數(shù)字圖像獲取裝置與獲取方法

成像設(shè)備置于玉米冠層頂部，距離地面6 m，由桁架結(jié)構(gòu)支撐固定，圖像采集指令由遠程上位機發(fā)出，無線通訊網(wǎng)絡(luò)負責指令和圖像數(shù)據(jù)的傳輸。成像單元采用JAI工業(yè)相機（AD-080CL，分辨率 1024×768像素，傳感器尺寸1/3英寸CCD，幀速30幀/s），鏡頭為富士能（HV-8M）12 mm紅外定焦鏡頭，電子光圈（F2.4-F22），工業(yè)相機上安裝有紅外/可見光成像轉(zhuǎn)換裝置，可根據(jù)外部光照條件或遠程指令獲取彩色圖像及近紅外黑白圖像（如圖1）。供試植株種植在北京市農(nóng)林科學(xué)院試驗田內(nèi)，2016年 6月23日播種，種植品種為zhongdan 2（60年代品種）和nongda108（90年代品種），種植密度為45 000株/hm2，正常水肥管理。從田間玉米播種后開始，每隔1 h獲取被監(jiān)測玉米小區(qū)可見光及紅外圖像，形成玉米整個生育期內(nèi)圖像序列（如圖2）。

圖1　玉米冠層紅外圖像序列獲取裝置示意圖Fig.1　Schematic diagram of corn canopy infrared image sequence acquisition device

圖2　序列圖像中不同時間節(jié)點圖像樣例Fig.2　Sample images of different time nodes in image sequences

冠層圖像處理程序由Visual Studio 2010開發(fā)，使用了圖像處理開源庫OpenCV 2.3，程序運行在PC機上，CPU核心頻率3.4 GHz，內(nèi)存4 GB。

1.2　冠層頂視圖像計算LAI的機理模型

冠層孔隙度（P0）與葉面積指數(shù) LAI的關(guān)系遵循Beer-Lambert定律[13]

式中P0(θ)是觀察方向與天頂夾角為θ時測量的冠層孔隙度；G(θ)是投影函數(shù)，即單位葉片面積在θ方向上的投影面積，?(θ)是葉片的叢生指數(shù)，代表了葉片隨機均勻分布變異程度。1/cos(θ)是光線穿過冠層的路徑長度。?(θ)LAI稱為有效葉面積指數(shù)（LAIe）[14-15]，從冠層孔隙度 P0(θ)計算 LAIe需要已知 G(θ)，文獻[16]提出測量一系列角度的冠層孔隙度，通過最優(yōu)化方法同時估計LAIe和G(θ)，文獻[17]對傾斜點樣方法進行分析研究，發(fā)現(xiàn)當葉片的平均傾斜角度約等于56°時，G(θ)函數(shù)在不同的天頂角θ方向上收斂于0.5，即此時G函數(shù)對天頂角的變化不敏感。這個角度約等于葉傾角為球狀分布時，葉片平均傾斜角度值。Goudriaan[18]和Spitters[19]指出球形葉片分布模型適合于描述包含玉米（不同生育期）在內(nèi)的多數(shù)作物冠層。因此當從垂直冠層頂部獲取冠層圖像計算得到孔隙度P0(0)時，cos(θ)值為 1，G(θ)值為 0.5，由式（1）可推導(dǎo)出下式

式中 P0(0)為垂直冠層頂部方向獲取的冠層圖像的孔隙度，LAIe為有效葉面積指數(shù)（間接方法測量結(jié)果均為LAIe，包括AccuPAR等設(shè)備，本文之后出現(xiàn)的LAI值，如無特殊說明均為有效葉面積指數(shù)）。

1.3　玉米冠層彩色圖像分割

冠層孔隙度P0的計算精度決定了LAIe的計算精度。冠層孔隙度的公式為

式中 Ps是冠層圖像上非植株像素數(shù)量，Pt是圖像上總像素數(shù)量。在田間條件下，冠層圖像植株及背景的顏色與亮度極易受到變化光照條件影響，文獻[20]中指出，不同光照條件對冠層孔隙度計算精度的影響可達 20%。冠層圖像獲取時間跨越整個生長季，變化的光照條件導(dǎo)致圖像呈現(xiàn)亮度差異，當光線較強時，植株頂部葉片有高光鏡面反射，當光線較弱時，植株底部葉片區(qū)域由于受光量少，在圖像中呈現(xiàn)陰影，植株上的高光和陰影區(qū)域容易被錯分為背景像素。Jeon和Lati等[21-22]指出田間光照條件下，植株上的高光鏡面反射和陰影對圖像分割是極大的挑戰(zhàn)。前人通過尋找優(yōu)化的顏色空間[23-25]，增強植株綠色分量比重[26-27]，及基于機器學(xué)習（machine learning）的圖像分割算法[28-29]進行冠層圖像分割。

支持向量機（SVM）是基于統(tǒng)計學(xué)習理論的模式識別方法，其主要思想是將低維不可分類向量映射到一個高維空間，在高維空間建立最大間隔超平面使得該平面兩側(cè)距平面最近的2個平行超平面距離最大化，SVM在小樣本、非線性及高維模式識別中具有一定優(yōu)勢。被監(jiān)測小區(qū)為3個，一個生長季90 d，每日5：00到19：00每隔1 h獲取一組圖像，在整體圖像集中，每3 d隨機從當日的圖像集合中抽取3幅，構(gòu)成90幅圖像訓(xùn)練樣本集，其他圖像作為測試集合。在90幅樣本圖像上，人工選擇植株像素和背景像素，構(gòu)建了24維的分類特征向量（如表1），計算兩類像素點特征向量各個分量值，作為SVM訓(xùn)練樣本輸入。

使用訓(xùn)練后的SVM分類器對圖2b進行分割，分割結(jié)果如圖3所示：

從分割結(jié)果圖像上看，葉片上的高光和遮擋產(chǎn)生的陰影導(dǎo)致了大量的錯誤分割結(jié)果，對圖2a-2f采用人工標注獲得植株圖像分割的真值，采用下式計算SVM分割方法精度

表1　植株圖像像素分類特征描述Table 1　Classifying features description of plant image pixels

圖3　07-28彩色圖像SVM分割結(jié)果Fig.3　SVM segmentation result of color image(07-28)

式中p 為像素灰度值，A是SVM方法分割得到的前景像素集(p=255)或背景像素集(p=0)，B為人工方法取得的圖像前景像素集(p=255)或背景像素集(p=0)，m、n分別為圖像的行數(shù)和列數(shù)，i、j分別為對應(yīng)的坐標。A和B的一致性越高，Qseg的值就越大，表明分割的精度就越高。圖2a-2f的均Qseq值為0.81?？梢钥闯?，圖像中的高光和陰影部分有較多的被錯誤分為背景，直接使用可見光彩色圖像序列計算冠層孔隙度誤差較大。

1.4　紅外圖像的高斯分布閾值分割方法

健康綠色植物在近紅外波段的光譜特征是反射率高（45%～50%），透過率高（45%～50%），吸收率低（＜5%）。在可見光波段與近紅外波段之間，即大約0.76 μm附近，反射率急劇上升，形成“紅邊”現(xiàn)象。玉米冠層紅外圖像中，由于植株和土壤背景的反射率不同，植株像素亮度高于圖像背景，存在于可見光彩色圖像中的葉片陰影區(qū)域由于反射率增高，亮度得到一定程度提升，彩色圖像中的高光部分在紅外波段光照強度削減，抑制了葉片高光反射現(xiàn)象。圖4給出了圖2c和2i的灰度圖像和紅外圖像的直方圖。

可見光圖像中土壤背景像素與植株像素的分界并不明顯，導(dǎo)致圖像分割時植株像素有可能被錯誤分割。在紅外圖像中，背景像素與植株像素呈現(xiàn)雙峰分布，分界較為明顯。

圖4　灰度圖像及紅外圖像的直方圖（08-03圖像）Fig.4 Histogram of gray image and infrared image(image on 08-03)

計算冠層圖像序列冠層孔隙度問題轉(zhuǎn)化為尋找冠層紅外圖像分割閾值問題，該類型的圖像灰度分布使用OTSU等閾值尋找算法結(jié)果并不理想，注意到土壤背景像素處于灰度直方圖的較暗一側(cè)，形狀上符合正態(tài)分布（背景像素點灰度值屬于大量隨機樣本）

若圖像分割閾值滿足超過 97.5%的背景像素能夠被正確分割，則

則能大大增加圖像二值化閾值分割的精度，計算直方圖中背景像素數(shù)量分布的峰值（第一個灰度分布峰值），記為μ，TL是像素數(shù)量累計為μ累積量的0.05時的像素值，TR為 TL相對于 μ的鏡像值，即 TR=2μ-TL，TR即為圖像的分割閾值?；诟咚狗植嫉膱D像分割計算簡便，適于在大量圖像序列上自動實施。

使用上述方法對圖2g-2l進行閾值分割，效果如圖5所示，圖5g-5l為灰度像素的直方圖。圖5a-5f的Qseg均值為0.94，分割精度得到較大提升。

圖5　紅外圖像灰度像素高斯分布閾值分割法效果Fig.5　Results of thresholding segmentation of infrared image of gray pixel Gaussian distribution

2　結(jié)果與分析

為了驗證本文所述裝置與方法測量LAI的精度，在2個品種（zhongdan2,nongda108）冠層生長過程中從苗期開始每間隔10 d，使用AccuPAR設(shè)備（Decagon Inc）測量同一冠層相同位置的LAI，對二者進行回歸分析，回歸直線如圖6所示。

圖6　本文與AccuPAR對LAI測量精度的回歸分析Fig.6　Regression analysis of LAI contained byAccuPAR and proposed method

回歸模型的決定系數(shù)為 0.94，證明本文所述方法測量LAI精度與AccuPAR接近，測量精度可滿足長期在線自動測量。

對 2016年田間種植的 zhongdan2、nongda108 品種使用本文所述的田間冠層連續(xù)監(jiān)測裝置及圖像處理算法計算LAI， zhongdan2為60年代的玉米栽培品種，冠層LAI連續(xù)測量值如圖7a所示，苗期到拔節(jié)期緩慢增大，拔節(jié)期后迅速增長，抽雄吐絲期LAI逐漸穩(wěn)定并達到峰值，吐絲后隨著營養(yǎng)生長向生殖生長過渡，葉片枯萎變黃，營養(yǎng)物質(zhì)向籽粒轉(zhuǎn)移，LAI呈現(xiàn)下降趨勢。nongda108為近現(xiàn)代栽培品種，測量結(jié)果如圖7b所示，LAI測量值變化趨勢前期與品種zhongdan2相近，后期（8月18日-9月7日）LAI下降趨勢更為緩慢，LAI維持在一定范圍內(nèi)持續(xù)時間較長。從 2個品種冠層的外觀上看，品種zhongdan2冠層后期葉片枯萎變黃趨勢更為明顯，品種nongda108葉片的“持綠性”較好，在籽粒形成和成熟過程中能夠提供更多的光合同化物質(zhì)，這也是近現(xiàn)代栽培品種產(chǎn)量高于歷史年代品種的一個原因。

圖7　不同年代玉米品種葉面積指數(shù)趨勢對比Fig.7　Comparison of LAI between two varieties in different decades

冠層覆蓋度（coverage fraction）定義為1減冠層空隙度，與冠層植株的干質(zhì)量和鮮質(zhì)量存在一定數(shù)量關(guān)系，連續(xù)監(jiān)測獲得冠層覆蓋度指標，可用于冠層植株生物量的無損估計。在上述 2個品種玉米植株生長過程中，從苗期開始每間隔 5 d于均一冠層中不同于紅外圖像監(jiān)測位置處破壞性采集植株地上部分（5株），分解后稱量鮮質(zhì)量，然后放置烘箱110 ℃殺青15 min后80 ℃烘干至恒重稱取干質(zhì)量，建立冠層覆蓋度與冠層植株鮮質(zhì)量/干質(zhì)量的曲線關(guān)系，如圖8所示。

圖8　冠層覆蓋度與生物量的關(guān)系Fig.8　Relationship of coverage fraction and biomass

從冠層覆蓋度與鮮質(zhì)量的預(yù)測模型決定系數(shù)上看，品種 zhongdan2高于品種 nongda108（0.96＞0.89），從圖上看模型的誤差主要來源于生長后期，由于品種zhongdan2后期葉片枯萎變黃較明顯，冠層覆蓋度與鮮質(zhì)量同步降低，而品種nongda108由于生長后期葉片持綠性較好，冠層覆蓋度下降趨勢平緩，葉片的鮮質(zhì)量仍在持續(xù)增加。冠層覆蓋度與干質(zhì)量的預(yù)測模型從決定系數(shù)上看，誤差主要來源于生長后期冠層覆蓋度下降，植株特別是籽粒干質(zhì)量逐漸增加。冠層覆蓋度的變化趨勢，基本可反映植株干質(zhì)量、鮮質(zhì)量的變化，可用于植株生物量變化的預(yù)測。

3　結(jié)論與討論

圖像法獲取冠層結(jié)構(gòu)參數(shù)的精度和穩(wěn)定性主要受到成像時光照條件的影響，使用SVM分類器對冠層彩色圖像分割的精度評價指標(Qseg)為 0.81，使用本文所述紅外冠層圖像結(jié)合背景正態(tài)分布分割方法，圖像序列分割精度評價指標(Qseg)為0.94，精度大大提高。在2個不同年代品種的田間實測對比試驗中，紅外圖像序列方法計算的LAI動態(tài)變化趨勢能夠很好地反應(yīng)2種不同株型玉米植株后期持綠性不同所導(dǎo)致冠層結(jié)構(gòu)差異。與目前間接測量原理的商業(yè)化冠層結(jié)構(gòu)測量設(shè)備AccuPAR進行精度對比，二者測量數(shù)據(jù)具有高度相關(guān)性，決定系數(shù)為0.94。

本文所采用的頂視圖孔隙度法估算冠層結(jié)構(gòu)參數(shù)（有效葉面積指數(shù)）具有以下幾個優(yōu)勢：

1）能夠自動化完成整個生育時期內(nèi)的測量，目前間接原理手持設(shè)備（如 AccuPAR）對測量的時間和光輻射強度有特定要求，不具備自動化功能，2）冠層圖像數(shù)據(jù)量大，包含內(nèi)容信息豐富易于系統(tǒng)集成，如利用顏色紋理信息進行作物營養(yǎng)診斷、利用形態(tài)信息進行水分供給狀況識別等。圖像法有利于建立起完整豐富的玉米冠層監(jiān)測系統(tǒng)。

冠層植株的“葉片重疊”造成真實葉面積指數(shù)與有效葉面積指數(shù)存在一定的偏差，這種偏差可通過計算孔隙度統(tǒng)計分布規(guī)律即“叢生指數(shù)”描述，結(jié)合品種及密度實驗的“叢生指數(shù)”計算將在下一步研究開展。

[參考文獻]

[1] Weiss M, Troufleau D, Baret F, et al. Coupling canopy functioning and radiative transfer models for remote sensing data assimilation[J]. Agricultural & Forest Meteorology,2001, 108(2): 113－128.

[2] 李云梅. 植被輻射傳輸理論與應(yīng)用[M]. 南京: 南京師范大學(xué)出版社, 2005.

[3] Kvet J, Marshall J K. Assessment of leaf area and other assimilating plant surfaces [C]//Plant Photosynthetic Production: Manual of Methods. The Netheriands: The Hague, 1971: 517－574.

[4] Song G Z M, Doley D, Yates D, et al. Improving accuracy of canopy hemispherical photography by a constant threshold value derived from an unobscured overcast sky[J]. Canadian Journal of Forest Research, 2013, 44(1): 17－27.

[5] Leblanc S G, Fournier R A. Hemispherical photography simulations with an architectural model to assess retrieval of leaf area index[J]. Agricultural and Forest Meteorology, 2014,194(6): 64－76.

[6] Woodgate W, Jones S D, Suarez L, et al. Understanding the variability in ground-based methods for retrieving canopy openness, gap fraction, and leaf area index in diverse forest systems[J]. Agricultural and Forest Meteorology, 2015,205(3): 83－95.

[7] Baret F, De S B, Lopezlozano R, et al. GAI estimates of row crops from downward looking digital photos taken perpendicular to rows at 57.5° zenith angle: Theoretical considerations based on 3D architecture models and application to wheat crops[J]. Agricultural & Forest Meteorology, 2010, 150(11):1393－1401.

[8] Gonsamo A, Jmn W, Pellikka P. CIMES: A package of programs for determining canopy geometry and solar radiation regimes through hemispherical photographs[J].Computers & Electronics in Agriculture, 2011, 79(2): 207－215.

[9] Campbell G S. Derivation of an angle density function for canopies with ellipsoidal leaf angle distributions[J].Agricultural and Forest Meteorology, 1990, 49(90): 173－176.

[10] Lang A R G. An instrument for measuring canopy structure[J]. Remote Sensing Reviews, 1990, 5(1): 61－71.

[11] 王傳宇，郭新宇，溫維亮，等. 田間光照條件下應(yīng)用半球圖像解析玉米冠層結(jié)構(gòu)參數(shù)[J]. 農(nóng)業(yè)工程學(xué)報, 2016,32(4): 157－162.Wang Chuanyu, Guo Xinyu, Wen Weiliang, et al.Application of hemispherical photography on analysis of maize canopy structural parameters under natural light[J].Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(4): 157－162. (in Chinese with English abstract)

[12] Sadeghi-Tehran P, Virlet N, Sabermanesh K, et al.Multi-feature machine learning model for automatic segmentation of green fractional vegetation cover for high-throughput field phenotyping[J]. Plant Methods, 2017,13(1): 103－119.

[13] Nilson T. A theoretical analysis of the frequency of gaps in plant stands[J]. Agricultural Meteorology, 1971, 8(71): 25－38.

[14] Chen J M, Black T A. Measuring leaf area index of plant canopies with branch architecture[J]. Agricultural & Forest Meteorology, 1991, 57(1–3): 1－12.

[15] Chen J M. Optically-based methods for measuring seasonal variation of leaf area index in boreal conifer stands[J].Agricultural & Forest Meteorology, 1996, 80(2-4): 135－163.

[16] Miller J B. A formula for average foliage density[J].Australian Journal of Botany, 1967, 15(1): 141－144.

[17] Wilson J W. Inclined point quadrats[J]. New Phytologist,1960, 59(1): 1－7.

[18] Goudriaan J. The bare bones of leaf-angle distribution in radiation models for canopy photosynthesis and energy exchange[J]. Agricultural & Forest Meteorology, 1988, 43(2):155－169.

[19] Spitters C J T. Separating the diffuse and direct component of global radiation and its implications for modeling canopy photosynthesis Part II. Calculation of canopy photosynthesis[J]. Agricultural & Forest Meteorology, 1986,38(1): 231－242.

[20] 胡凝, 呂川根, 姚克敏, 等. 利用魚眼影像技術(shù)反演不同株型水稻的冠層結(jié)構(gòu)參數(shù)[J]. 作物學(xué)報, 2014, 40(8): 1443－1451.Hu Ning, Lǖ Chuangen, Yao Kemin, et al. Inversion of rice canopy construction parameters from the Hemispherical photograph[J]. Acta Gronomica Sinica, 2014, 40(8): 1443－1451. (in Chinese with English abstract)

[21] Jeon H Y, Tian L F, Zhu H. Robust crop and weed segmentation under uncontrolled outdoor illumination[J].Sensors, 2011, 11(6): 6270－6283.

[22] Lati R N, Filin S, Eizenberg H. Robust methods for measurement of leaf-cover area and biomass from image data[J]. Weed Science, 2011, 59(2): 276－284.

[23] Philipp I, Rath T. Improving plant discrimination in image processing by use of different colour space transformations[J].Computers & Electronics in Agriculture, 2002, 35(2): 1－15.

[24] Liu Y, Mu X, Wang H, et al. A novel method for extracting green fractional vegetation cover from digital images[J].Journal of Vegetation Science, 2012, 23(3): 406－418.

[25] Panneton B, Brouillard M. Colour representation methods for segmentation of vegetation in photographs[J]. Biosystems Engineering, 2009, 102(4): 365－378.

[26] Meyer G E, Neto J C. Verification of color vegetation indices for automated crop imaging applications[J]. Computers &Electronics in Agriculture, 2008, 63(2): 282－293.

[27] Burgos-Artizzu X P, Ribeiro A, Guijarro M, et al. Real-time image processing for crop/weed discrimination in maize fields[J]. Computers & Electronics in Agriculture, 2011,75(2): 337－346.

[28] Ruiz-Ruiz G, Gómez-Gil J, Navas-Gracia L M. Testing different color spaces based on hue for the environmentally adaptive segmentation algorithm (EASA)[J]. Computers &Electronics in Agriculture, 2009, 68(1): 88－96.

[29] Zheng L Y, Tian K. Segmentation of Green Vegetation from Crop Canopy Images Based on Fisher Linear Discriminant[J].Key Engineering Materials, 2012, 500: 487－491.