森林资源数据采集是森林资源监测的基础, 提高数据采集的效率和质量是森林资源信息化的重要内容^[1-2]。在各类森林资源调查中, 立木的胸径测量是最重要的内容之一, 其精度直接影响森林资源蓄积量的消长等情况^[3]。传统的人工接触式测量方法^[4]劳动强度大、人力成本高、效率低。使用全站仪等精密测量设备存在操作复杂、不易携带、成本较高等问题。近景摄影测量方法^[5-6]虽然可以很好地解决这些问题, 但标定算法复杂, 且需要在计算机上操作完成, 步骤繁琐。目前国内外学者在立木测量方法上做了大量的研究并取得了诸多成果, 开发出了一系列软、硬件系统, 但这些系统及其算法、模型多数是针对计算机和数码相机, 无法在较为普适的移动终端设备上完成测量^[7-8]。近年来, 随着智能手机的快速发展, 利用智能手机来进行立木测量成为可能, 许多学者基于Android平台开发了测树软件^[9-10], 这些软件基于三角函数原理测量树高, 在测量胸径时需要根据已知尺寸的标定物进行对比计算, 一定程度上实现了树高、胸径因子的实时测量, 但是操作复杂, 精度不高。针对上述问题, 本研究以立木胸径(DBH)的快速、自动获取为目标, 提出一种借助普适率高的智能移动端设备(智能手机), 结合机器视觉技术、相机标定技术和三维重建技术的立木胸径自动测量方法, 以期提升森林资源外业调查的作业精度和效率, 同时可让非专业人士借助非专业装备完成立木胸径的测量。

1.   测量方法设计

1.1   胸径定义

立木胸径是森林资源调查中的一个重要因子, 指乔木主干离地表面1.3 m处的直径。胸径测量可为林业相关作业提供基础数据和科学依据, 如林分蓄积量测定、林木采伐、林区精准施药、林区地图绘制等。同时, 林区立木的胸径生长数据反映了林间周围的气候、光照、土地营养状态等环境条件。

1.2   胸径测量方法设计

本研究提出的立木胸径测量方法流程如图 1所示。首先通过智能手机相机采集待测立木(悬挂棋盘格标定板)图像, 根据智能手机相机特点建立一种改进的带有非线性畸变项的相机标定模型, 矫正图像中各像素的偏移, 标定高精度的相机内、外参数。利用相机内、外参数和二维图像坐标信息进行三维世界坐标重建, 计算单位像素对应实际物理尺寸; 同时, 采用Lab(CIE Lab color model)颜色模型和HSV(色调H, 饱和度S, 明度V)颜色空间模型作为图像特征, 分别提取图像的L(亮度), a(洋红色至绿色范围)和b(黄色至蓝色范围)3个分量及H分量, 利用各像素点(x, y)的Lab三分量的算数平均值与各像素经3 × 3算子卷积运算得到的下采样图之间的差异, 得到高频图。将色调H分量均衡化处理后进行特征融合, 得到融合图, 从而增强立木轮廓特征, 再将融合图二值化。利用膨胀和腐蚀组合运算进行去噪和平滑处理, 并根据立木主轮廓和干扰轮廓周长差距大的特点, 提取主轮廓, 完成立木轮廓的图像分割; 提取立木轮廓的最小外接矩形, 获取矩形树干高度方向的像素参数。由单位像素物理尺寸信息和树干高度方向的像素参数信息, 获取树高1.3 m位置处的胸径像素参数。结合树高1.3 m处胸径方向像素信息和单位像素物理尺寸信息, 完成胸径计算。

图  1  采用相机标定技术和机器视觉技术的胸径自动测量方法流程图

Figure  1.  Flow chart of automatic DBH measurement method based on camera calibration technology and machine vision technology

下载: 全尺寸图片 幻灯片

2.   立木胸径测量关键技术

2.1   智能手机相机标定模型及畸变矫正模型

现有的相机标定方法主要分为传统标定物标定法、基于主动视觉标定法和自标定法。为了适应移动端相机镜头组的特点, 使之能够实现畸变矫正并获取高精度相机参数, 本研究提出一种经改进的带有非线性畸变项的相机标定模型来实现相机标定。下式分别描述了畸变矫正模型(1), 忽略了高阶项的相机径向、切向畸变模型函数(2)和(3)及由式(1)~式(3)得到的畸变矫正模型函数(4)。

其中(x_u, y_u)为线性针孔模型计算出的理想图像坐标点, (x, y)是实际的图像坐标点, (δ_x, δ_y)是非线性畸变值, 它与坐标的位置有关。(δ_xr, δ_yr)为忽略了高价项的非线性径向畸变值, (δ_xd, δ_yd)为忽略了高价项的非线性切向畸变值。k₁, k₂, p₁, p₂代表 4个非线性畸变系数。

图像坐标系上的点(x, y)在像素坐标系(u, v)中的坐标为(u₀, v₀), 像面上每一个像素点在x轴, y轴方向上的物理尺寸为d_x, d_y。可知图像中任意一个像素在2个坐标系中满足如下关系:

结合式(5), 线性相机模型可化为齐次坐标与矩阵形式(6), 其中(x_w, y_w, z_w)表示三维世界坐标下的点坐标, f表示相机焦距, R表示旋转矩阵, T表示平移向量。

2.2   相机参数计算及优化

将棋盘格标定板固定于立木树干上, 通过智能手机获取标定图片(图 2), 进行相机内、外参数标定。根据ZHANG^[11]标定法中的相机标定解决方法。先利用模板平面与图像间的单应性关系作为约束条件获得封闭解, 在封闭解的基础上计算初始估值内、外参数矩阵; 然后根据最大似然估计理论给出非线性的最优化解, 并考虑透镜的径向畸变, 得到解析解和非线性解; 最后应用非线性最小二乘的Levenberg-Marquardt算法来迭代获得更精确的相机内、外参数。选取1组(12幅)标定模板进行标定实验的数据, 其中图像的各角点反投影误差和每幅模板图像的平均像素误差如图 3所示。

图  2  标定图片获取示意图

Figure  2.  Schematic drawing of calibration picture

下载: 全尺寸图片 幻灯片

图  3  标定模板图像的角点和平均像素误差

Figure  3.  Corner reprojection errors and average pixels errors of calibration images

图 3A的x和y坐标分别代表模板图像的2个方向，图中的各点表示模板图像中各角点在x和y方向上的偏移方向和偏移量。图 3B的横坐标表示模板图像的编号，纵坐标表示每幅图像中所有角点的平均像素误差

下载: 全尺寸图片 幻灯片

相机标定方法能获取相机的内、外参数, 而畸变矫正能解决智能手机相机镜头组存在畸变的问题, 对图片中的所有像素进行偏移纠正, 获得经畸变矫正的图像。

2.3   由图像坐标系重建三维世界坐标系

经过相机标定, 可计算图像上2个像素点在三维世界坐标系下的真实物理距离。算法过程为:根据智能手机的相机成像传感器尺寸L_CCD推出每个像素的物理尺寸L_pixels; 由2.2节获取的经畸变矫正后的手机相机参数, 计算摄像机焦距f_c的物理尺寸L_fc; 然后计算像素平面上待测目标物两点间的距离L_distance; 最后可根据二维图像信息计算对应三维世界中2个像素点之间的真实物理距离L_3D-real。其中P_sum表示待测目标两点之间的像素个数; R和T分别表示相机标定所得外参数中的旋转矩阵和平移向量。三维世界重建算法(图 4)可计算单位像素对应的三维世界坐标的物理尺寸。

图  4  三维世界坐标重建算法流程图

Figure  4.  Flow chart of 3D coordinate reconstruction algorithm

下载: 全尺寸图片 幻灯片

2.4   立木树干轮廓视觉显著性表达及色调分量融合

相比工业环境下的机器视觉图像分割, 自然环境下立木轮廓的视觉分割具有目标对象形态各异, 图像背景干扰因素多和光照不均匀等特点, 将现有的图像分割方法直接应用于自然环境进行感兴趣区域检测及图像分割的效果不佳。本研究以自然环境中的立木图像为对象, 提出一种结合频率调谐的视觉显著性^[12]的轮廓检测方法, 具体为:①颜色空间选取。利用Lab颜色模型中L, a和b 3个分量的几何距离差异, 区分图像中不同颜色区域; 由于3个分量之间独立性高, 在目标显著性的特征表达上效果较好。因此, 本研究引入对自然环境中立木轮廓目标区域的频率调谐视觉显著性描述, 使之成为识别目标的特征之一。Lab三分量可通过式(7)提取, 其中R, G和B分别代表该目标像素的红色, 绿色和蓝色值。

② 频率调谐的视觉显著性表达^[13-14]。以Lab颜色空间为图像特征, 计算各颜色通道(L, a, b)上每个像素点(x, y)与整幅图像的平均色差并取平方; 然后将这3个通道的值相加作为该像素的显著性值。采用3 × 3算子对图像进行卷积运算, 得到一次下采样图, 并构建高斯金字塔, 对图像进行多次高斯平滑处理, 最终得到高频图。

式(8)中:S(x, y)表示像素点(x, y)的频率调谐视觉显著性值, I_μ为图像特征的几何平均向量, I_ωhc为对原始图像的高斯模糊, 采用3 × 3算子; (L, a, b)各分量下标μ表示采集原始图像下各特征分量的算术平均值, 下标ωhc表示高斯滤波后图像每个像素点的特征分量。将此方法应用于目标立木轮廓提取的特征表达, 得到的视觉显著图(图 5B)能有效捕捉到各颜色空间分量的剧烈跳变点, 突出立木轮廓。③色调分量均衡化与特征融合。仅依靠Lab颜色模型构建立木图像的视觉显著图, 受光照影响较大。如果在视觉显著图的基础上对图像进行色调分量均衡化与特征融合处理, 可提取较高质量立木树干轮廓; 如将HSV颜色空间中H分量融合到目标立木树干特征表达中, 可以弥补在Lab空间上表达立木显著性时光照不均匀所带来的影响。提取图像H通道分量(图 5C), 通过对比度受限自适应直方图均衡化调整^[15], 增强图像中立木树干部分的颜色对比度, 捕获更多棕褐色系的立木树干与绿色系的背景之间的细节差异。将立木树干图像的均衡化H分量与频率调谐视觉显著图通过式(9)进行融合:

图  5  融入H分量的特征表达图及树干轮廓提取示意图

Figure  5.  Feature expression diagram integrating into H components and a sketch map of tree trunk extraction

下载: 全尺寸图片 幻灯片

式(9)中:I_fusion(x, y)表示每一像素点(x, y)的融合特征, H(x, y)表示像素点(x, y)的均衡化色调, S(x, y)该像素点(x, y)的频率调谐视觉显著性值。通过特征调制, 减少光照强度变化对于视觉显著性的影响, 从而增强目标立木树干轮廓。融入H分量的特征表达图(图 5D)。

2.5   噪声处理及树干轮廓检测

分割后的图像仍会呈现出立木轮廓不完整现象和残留部分噪声, 包括树干纹理造成的空隙、树叶遮挡后不完整的轮廓等。利用形态学膨胀和腐蚀组合运算, 通过开操作(先腐蚀后膨胀)分割出独立的图像元素, 消除细小物体, 达到在纤细处分离物体和平滑较大物体边界的效果。通过闭操作(先膨胀后腐蚀)连接图像中相邻的元素, 填充物体内的细小空洞, 达到连接邻近物体和平滑边界的效果。由于背景复杂, 上述去噪处理后图像仍然可能存在较大的干扰噪声块, 利用立木主轮廓和干扰轮廓的周长差距大的特点, 剔除周长较小的干扰项, 保留主轮廓, 完成立木轮廓的输出(图 6)。

图  6  树干轮廓提取示意图

Figure  6.  Sketch map of tree trunk outline extraction

下载: 全尺寸图片 幻灯片

在完成立木轮廓检测后, 可通过提取立木轮廓的最小外接矩形^[16]获得立木树干高度的属性, 在2.3节的基础上进行树干高度(指所获待测图像中树干部分的高度, 并非树高)参数提取, 并获取树干高度方向的像素值。由于目标立木在图像中存在形状不规则和位置、方向自由等特点, 使得难以准确定位立木位置及提取立木的边缘轮廓, 因此通过目标主轴方向估计求出立木轮廓最小外接矩形。根据主轴确定目标的初始矩形, 以初始矩形的几何中心为旋转中心对其进行旋转, 从而找到最优矩形姿态, 并平移优化旋转后的矩形(图 7)。

图  7  立木轮廓最小外接矩形提取及胸径位置定位和胸径计算

Figure  7.  Minimum outer rectangleof tree trunk outline extraction, DBH location and DBH calculation

下载: 全尺寸图片 幻灯片

2.6   胸径检测与计算

提取立木轮廓的最小外接矩形, 将矩形中树高方向定义为树干高度方向, 获取最小外接中树干高度方向的像素参数P_H, 结合经三维重建测得的单位像素对应的三维世界的物理尺寸信息L_3D-real, 即可计算树干高度:H=P_H × L_3D-real。胸径(DBH)的测量算法如下:由单位像素物理尺寸信息和树干高度方向的像素参数信息L_3D-real, 获取树高1.3 m位置处的胸径像素参数:P_DBH=1 300/L_3D-real, 结合树高1.3 m处胸径方向像素信息PDBH和单位像素物理尺寸信息L_3D-real, 可得胸径为:D_BH= P_DBH×L_3D-real。此处, P_H与P_DBH指在未经采样处理的原图像中提取的像素信息。

3.   实验验证

基于上述原理和算法流程, 本研究在Android系统平台上, 利用Java语言和C++语言等开发适用于Android智能手机的树干高度、胸径测量的测试系统, 用于本研究方法的精度验证。

3.1   树干高度测量精度验证

随机选定12株立木, 编号后用皮尺测量各树干高度; 利用本研究的方法测量图片中的立木高度以验证测量精度。使用经校准的实验手机对12株立木进行测量, 相对误差均小于4.65%(表 1)。

表  1  树干高度测量数据

Table  1.  Measurement data of the trunk height

样本号	树干高度真值/m	树干高度测量值/m	相对误差/%
1	2.700	2.715±0.148	3.89
2	2.800	2.665±0.007	4.82
3	3.600	3.620±0.099	3.95
4	3.900	4.000±0.014	2.57
5	4.200	4.170±0.226	3.81
6	4.400	4.225±0.049	3.98
7	4.400	4.605±0.007	4.65
8	4.900	5.085±0.064	3.78
9	5.200	5.225±0.163	2.21
10	5.500	5.340±0.057	2.91
11	5.800	5.670±0.028	2.25
12	6.100	6.330±0.028	3.77
说明：树干高度指所获取图像内的树干部分的高度，并非指立木因子中的树高

下载: 导出CSV 

| 显示表格

3.2   胸径测量精度验证

随机选定10株立木, 编号后用轮尺测量各立木1.3 m高度处的胸径作为真实胸径值。用实验手机对各样本多次测量, 计算平均值作为测量值。由表 2可知:10个样本的胸径测量结果相对误差均小于2.50%, 绝对误差均小于5 mm, 符合森林资源调查工作的精度要求。

表  2  胸径测量数据

Table  2.  Measurement data of the DBH

样本号	胸径真值/ mm	胸径测量值/mm	平均绝对误差/mm	相对误差/%
1	163.000	161.000±0.265	-2.000	1.23
2	148.000	145.900±0.400	-2.100	1.42
3	122.000	119.900±0.557	-2.100	1.72
4	152.000	149.467±0.404	-2.500	1.67
5	188.000	183.533±1.343	-4.500	2.37
6	132.000	129.233±1.242	-2.800	2.10
7	139.000	137.367±0.611	-1.600	1.17
8	203.000	200.500±0.200	-2.500	1.23
9	204.000	201.367±0.814	-2.600	1.29
10	228.000	225.933±0.551	-2.100	0.77

下载: 导出CSV 

| 显示表格

4.   结束语

本研究提出的结合相机标定技术、机器视觉技术及三维重建技术的立木胸径自动测量方法, 能够更加方便地实现林业智能测量, 有效提升森林资源外业调查的作业精度和效率。与现有的基于Android平台的测树软件相比, 该方法的优势主要体现在:①结合相机标定技术, 消除了相机镜头存在的畸变, 减少了系统误差; ②结合机器视觉技术, 将图像处理技术应用于林业资源调查工作, 在实现胸径快速测量的同时, 获取并保存了立木图像信息; ③相比三角函数方法, 本方法的测量效率和精度更高; 相比比例模型测量胸径的方法, 本方法使用图像增强, 感兴趣区域检测及相机标定和三维重建等技术, 使得鲁棒性和测量精度更高。本研究虽然较好地实现了立木胸径的快速测量, 但仅能在一定距离内获取简单的直接测树因子, 希望在下一步的研究中, 实现立木高度、森林蓄积量、生物量等数据的自动计算, 为“智慧林业”的数据采集提供支持。