基于时间序列红外图像的玉米叶面积指数连续监测

王传宇，郭新宇，杜建军

（1. 北京农业信息技术研究中心，北京 100097；2. 数字植物北京市重点实验室，北京 100097）

0　引　言

叶面积指数（leaf area index，LAI）影响作物光合辐射截获量，水分的吸收，潜热和感热通量，地表生态系统与大气的 CO2交换。在一些作物生长模型中，LAI也是重要输入参数和条件变量。LAI是冠层辐射传输模型的决定因素，在遥感数据同化中，已被作为冠层反射与作物生长模型的关键联系[1]。测量作物整个生育时期的LAI指数时间变化规律，有助于作物产量、生物量预测模型的校正。

LAI获取方法主要分为直接法和间接法2类，所谓直接法就是对直接组成部分进行量测，例如通过测量采样区域植株的叶片面积计算 LAI[2-3]，直接法需要大范围地破坏性采样，消耗大量人力物力，测量指标主观依赖性强，因此实际中间接方法应用的更为广泛深入[4-6]，Baret等[7]指出，冠层内的辐射传播主要决定于LAI和叶倾角分布函数，基于此，发展出了由冠层孔隙度（gap fraction）计算LAI的方式。冠层空隙度的获取方法分为2种：1）测量冠层下方多个角度的太阳直接辐射值，如商业化冠层分析仪AccuPAR（Decagon Inc）；2）采用图像同时获得植株与天空/地面背景的像素投影，计算植株背景像素占整幅图像的比例。其中半球图像（hemispherical photography）能够记录冠层内部多个角度的隙度，被众多学者研究和应用[8-10]。半球图像的投影模型较为复杂[11]，对于处于营养生长前期的禾本科作物玉米，测量误差较大，不能适应整个生育时期玉米冠层LAI动态变化监测的需求。有研究表明，可以采用冠层顶部单角度图像获得冠层孔隙度进而计算 LAI[12]，简化了图像投影模型算法。开展整个生育时期内玉米冠层顶部图像LAI动态监测，计算精度主要受限于图像分割算法，玉米植株自身形态及外部光照环境的显著变化是影响图像分割和冠层孔隙度计算精度的主要原因。针对上述研究中存在的不足，本文将研究一种单角度数字图像LAI连续获取方法，为解决田间变化光照对绿色植株分割的影响，采用红外图像去除植株叶片高光反射和阴影，基于高斯分布模型自动计算图像分割阈值，推导了冠层顶视图像孔隙度计算LAI的方法，同时建立了植株图像像素与植株干质量和鲜质量的模型，并在多个图像序列上对算法效果进行分析验证，以期为田间环境下冠层参数的自动连续监测提供解决方案。

1　材料与方法

1.1　数字图像获取装置与获取方法

成像设备置于玉米冠层顶部，距离地面6 m，由桁架结构支撑固定，图像采集指令由远程上位机发出，无线通讯网络负责指令和图像数据的传输。成像单元采用JAI工业相机（AD-080CL，分辨率 1024×768像素，传感器尺寸1/3英寸CCD，帧速30帧/s），镜头为富士能（HV-8M）12 mm红外定焦镜头，电子光圈（F2.4-F22），工业相机上安装有红外/可见光成像转换装置，可根据外部光照条件或远程指令获取彩色图像及近红外黑白图像（如图1）。供试植株种植在北京市农林科学院试验田内，2016年 6月23日播种，种植品种为zhongdan 2（60年代品种）和nongda108（90年代品种），种植密度为45 000株/hm2，正常水肥管理。从田间玉米播种后开始，每隔1 h获取被监测玉米小区可见光及红外图像，形成玉米整个生育期内图像序列（如图2）。

图1　玉米冠层红外图像序列获取装置示意图Fig.1　Schematic diagram of corn canopy infrared image sequence acquisition device

图2　序列图像中不同时间节点图像样例Fig.2　Sample images of different time nodes in image sequences

冠层图像处理程序由Visual Studio 2010开发，使用了图像处理开源库OpenCV 2.3，程序运行在PC机上，CPU核心频率3.4 GHz，内存4 GB。

1.2　冠层顶视图像计算LAI的机理模型

冠层孔隙度（P0）与叶面积指数 LAI的关系遵循Beer-Lambert定律[13]

式中P0(θ)是观察方向与天顶夹角为θ时测量的冠层孔隙度；G(θ)是投影函数，即单位叶片面积在θ方向上的投影面积，Ω(θ)是叶片的丛生指数，代表了叶片随机均匀分布变异程度。1/cos(θ)是光线穿过冠层的路径长度。Ω(θ)LAI称为有效叶面积指数（LAIe）[14-15]，从冠层孔隙度 P0(θ)计算 LAIe需要已知 G(θ)，文献[16]提出测量一系列角度的冠层孔隙度，通过最优化方法同时估计LAIe和G(θ)，文献[17]对倾斜点样方法进行分析研究，发现当叶片的平均倾斜角度约等于56°时，G(θ)函数在不同的天顶角θ方向上收敛于0.5，即此时G函数对天顶角的变化不敏感。这个角度约等于叶倾角为球状分布时，叶片平均倾斜角度值。Goudriaan[18]和Spitters[19]指出球形叶片分布模型适合于描述包含玉米（不同生育期）在内的多数作物冠层。因此当从垂直冠层顶部获取冠层图像计算得到孔隙度P0(0)时，cos(θ)值为 1，G(θ)值为 0.5，由式（1）可推导出下式

式中 P0(0)为垂直冠层顶部方向获取的冠层图像的孔隙度，LAIe为有效叶面积指数（间接方法测量结果均为LAIe，包括AccuPAR等设备，本文之后出现的LAI值，如无特殊说明均为有效叶面积指数）。

1.3　玉米冠层彩色图像分割

冠层孔隙度P0的计算精度决定了LAIe的计算精度。冠层孔隙度的公式为

式中 Ps是冠层图像上非植株像素数量，Pt是图像上总像素数量。在田间条件下，冠层图像植株及背景的颜色与亮度极易受到变化光照条件影响，文献[20]中指出，不同光照条件对冠层孔隙度计算精度的影响可达 20%。冠层图像获取时间跨越整个生长季，变化的光照条件导致图像呈现亮度差异，当光线较强时，植株顶部叶片有高光镜面反射，当光线较弱时，植株底部叶片区域由于受光量少，在图像中呈现阴影，植株上的高光和阴影区域容易被错分为背景像素。Jeon和Lati等[21-22]指出田间光照条件下，植株上的高光镜面反射和阴影对图像分割是极大的挑战。前人通过寻找优化的颜色空间[23-25]，增强植株绿色分量比重[26-27]，及基于机器学习（machine learning）的图像分割算法[28-29]进行冠层图像分割。

支持向量机（SVM）是基于统计学习理论的模式识别方法，其主要思想是将低维不可分类向量映射到一个高维空间，在高维空间建立最大间隔超平面使得该平面两侧距平面最近的2个平行超平面距离最大化，SVM在小样本、非线性及高维模式识别中具有一定优势。被监测小区为3个，一个生长季90 d，每日5：00到19：00每隔1 h获取一组图像，在整体图像集中，每3 d随机从当日的图像集合中抽取3幅，构成90幅图像训练样本集，其他图像作为测试集合。在90幅样本图像上，人工选择植株像素和背景像素，构建了24维的分类特征向量（如表1），计算两类像素点特征向量各个分量值，作为SVM训练样本输入。

使用训练后的SVM分类器对图2b进行分割，分割结果如图3所示：

从分割结果图像上看，叶片上的高光和遮挡产生的阴影导致了大量的错误分割结果，对图2a-2f采用人工标注获得植株图像分割的真值，采用下式计算SVM分割方法精度

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

图3　07-28彩色图像SVM分割结果Fig.3　SVM segmentation result of color image(07-28)

式中p 为像素灰度值，A是SVM方法分割得到的前景像素集(p=255)或背景像素集(p=0)，B为人工方法取得的图像前景像素集(p=255)或背景像素集(p=0)，m、n分别为图像的行数和列数，i、j分别为对应的坐标。A和B的一致性越高，Qseg的值就越大，表明分割的精度就越高。图2a-2f的均Qseq值为0.81。可以看出，图像中的高光和阴影部分有较多的被错误分为背景，直接使用可见光彩色图像序列计算冠层孔隙度误差较大。

1.4　红外图像的高斯分布阈值分割方法

健康绿色植物在近红外波段的光谱特征是反射率高（45%～50%），透过率高（45%～50%），吸收率低（＜5%）。在可见光波段与近红外波段之间，即大约0.76 µm附近，反射率急剧上升，形成“红边”现象。玉米冠层红外图像中，由于植株和土壤背景的反射率不同，植株像素亮度高于图像背景，存在于可见光彩色图像中的叶片阴影区域由于反射率增高，亮度得到一定程度提升，彩色图像中的高光部分在红外波段光照强度削减，抑制了叶片高光反射现象。图4给出了图2c和2i的灰度图像和红外图像的直方图。

可见光图像中土壤背景像素与植株像素的分界并不明显，导致图像分割时植株像素有可能被错误分割。在红外图像中，背景像素与植株像素呈现双峰分布，分界较为明显。

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

计算冠层图像序列冠层孔隙度问题转化为寻找冠层红外图像分割阈值问题，该类型的图像灰度分布使用OTSU等阈值寻找算法结果并不理想，注意到土壤背景像素处于灰度直方图的较暗一侧，形状上符合正态分布（背景像素点灰度值属于大量随机样本）

若图像分割阈值满足超过 97.5%的背景像素能够被正确分割，则

则能大大增加图像二值化阈值分割的精度，计算直方图中背景像素数量分布的峰值（第一个灰度分布峰值），记为μ，TL是像素数量累计为μ累积量的0.05时的像素值，TR为 TL相对于 μ的镜像值，即 TR=2μ-TL，TR即为图像的分割阈值。基于高斯分布的图像分割计算简便，适于在大量图像序列上自动实施。

使用上述方法对图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　结果与分析

为了验证本文所述装置与方法测量LAI的精度，在2个品种（zhongdan2,nongda108）冠层生长过程中从苗期开始每间隔10 d，使用AccuPAR设备（Decagon Inc）测量同一冠层相同位置的LAI，对二者进行回归分析，回归直线如图6所示。

图6　本文与AccuPAR对LAI测量精度的回归分析Fig.6　Regression analysis of LAI contained byAccuPAR and proposed method

回归模型的决定系数为 0.94，证明本文所述方法测量LAI精度与AccuPAR接近，测量精度可满足长期在线自动测量。

对 2016年田间种植的 zhongdan2、nongda108 品种使用本文所述的田间冠层连续监测装置及图像处理算法计算LAI， zhongdan2为60年代的玉米栽培品种，冠层LAI连续测量值如图7a所示，苗期到拔节期缓慢增大，拔节期后迅速增长，抽雄吐丝期LAI逐渐稳定并达到峰值，吐丝后随着营养生长向生殖生长过渡，叶片枯萎变黄，营养物质向籽粒转移，LAI呈现下降趋势。nongda108为近现代栽培品种，测量结果如图7b所示，LAI测量值变化趋势前期与品种zhongdan2相近，后期（8月18日-9月7日）LAI下降趋势更为缓慢，LAI维持在一定范围内持续时间较长。从 2个品种冠层的外观上看，品种zhongdan2冠层后期叶片枯萎变黄趋势更为明显，品种nongda108叶片的“持绿性”较好，在籽粒形成和成熟过程中能够提供更多的光合同化物质，这也是近现代栽培品种产量高于历史年代品种的一个原因。

图7　不同年代玉米品种叶面积指数趋势对比Fig.7　Comparison of LAI between two varieties in different decades

冠层覆盖度（coverage fraction）定义为1减冠层空隙度，与冠层植株的干质量和鲜质量存在一定数量关系，连续监测获得冠层覆盖度指标，可用于冠层植株生物量的无损估计。在上述 2个品种玉米植株生长过程中，从苗期开始每间隔 5 d于均一冠层中不同于红外图像监测位置处破坏性采集植株地上部分（5株），分解后称量鲜质量，然后放置烘箱110 ℃杀青15 min后80 ℃烘干至恒重称取干质量，建立冠层覆盖度与冠层植株鲜质量/干质量的曲线关系，如图8所示。

图8　冠层覆盖度与生物量的关系Fig.8　Relationship of coverage fraction and biomass

从冠层覆盖度与鲜质量的预测模型决定系数上看，品种 zhongdan2高于品种 nongda108（0.96＞0.89），从图上看模型的误差主要来源于生长后期，由于品种zhongdan2后期叶片枯萎变黄较明显，冠层覆盖度与鲜质量同步降低，而品种nongda108由于生长后期叶片持绿性较好，冠层覆盖度下降趋势平缓，叶片的鲜质量仍在持续增加。冠层覆盖度与干质量的预测模型从决定系数上看，误差主要来源于生长后期冠层覆盖度下降，植株特别是籽粒干质量逐渐增加。冠层覆盖度的变化趋势，基本可反映植株干质量、鲜质量的变化，可用于植株生物量变化的预测。

3　结论与讨论

图像法获取冠层结构参数的精度和稳定性主要受到成像时光照条件的影响，使用SVM分类器对冠层彩色图像分割的精度评价指标(Qseg)为 0.81，使用本文所述红外冠层图像结合背景正态分布分割方法，图像序列分割精度评价指标(Qseg)为0.94，精度大大提高。在2个不同年代品种的田间实测对比试验中，红外图像序列方法计算的LAI动态变化趋势能够很好地反应2种不同株型玉米植株后期持绿性不同所导致冠层结构差异。与目前间接测量原理的商业化冠层结构测量设备AccuPAR进行精度对比，二者测量数据具有高度相关性，决定系数为0.94。

本文所采用的顶视图孔隙度法估算冠层结构参数（有效叶面积指数）具有以下几个优势：

1）能够自动化完成整个生育时期内的测量，目前间接原理手持设备（如 AccuPAR）对测量的时间和光辐射强度有特定要求，不具备自动化功能，2）冠层图像数据量大，包含内容信息丰富易于系统集成，如利用颜色纹理信息进行作物营养诊断、利用形态信息进行水分供给状况识别等。图像法有利于建立起完整丰富的玉米冠层监测系统。

冠层植株的“叶片重叠”造成真实叶面积指数与有效叶面积指数存在一定的偏差，这种偏差可通过计算孔隙度统计分布规律即“丛生指数”描述，结合品种及密度实验的“丛生指数”计算将在下一步研究开展。

[参考文献]

[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] 李云梅. 植被辐射传输理论与应用[M]. 南京: 南京师范大学出版社, 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] 王传宇，郭新宇，温维亮，等. 田间光照条件下应用半球图像解析玉米冠层结构参数[J]. 农业工程学报, 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] 胡凝, 吕川根, 姚克敏, 等. 利用鱼眼影像技术反演不同株型水稻的冠层结构参数[J]. 作物学报, 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.