CN109285145A - 基于智能手机的多株立木高度测量方法 - Google Patents

Abstract

本发明提供一种基于智能手机的多株立木高度测量方法，立木高度测量时主要是利用智能终端的移动设备的图像采集功能和设备内置重力传感器的角度测量功能，将手机相机拍摄棋盘标定板进行标定和相机参数获取，利用非线性畸变模型对采集立木图片进行非线性畸变校正；利用图像分割处理，通过灰度化和二值化处理，识别标定物特征点；利用标定物特征点在图像中的像素位置，建立设标立木平面坐标系和设标立木高度计算模型，获取设标立木高度；分析图像中立木图像像素、物体成像角度和相机旋转角度，建立立木深度信息提取模型，获取待测立木深度信息；建立多株立木高度测量模型，将待测多株立木垂直最值点代入，进行多株立木高度测量。

Description

基于智能手机的多株立木高度测量方法

技术领域

本发明涉及的一种测量方法，具体涉及一种基于智能手机的多株立木高度测量方法。

背景技术

立木高度是评价立地质量与林木生长状况的重要依据。立木高度测量分为人工接触性测量和非接触性测量。目前，树高的测量仪器主要有布鲁莱斯测高器、电子经纬仪、全站仪等，但是这些测量仪器存在操作复杂、不易携带、成本较高等问题。目前国内外学者针对以上问题，基于近景摄影测量原理研究立木高度测量，开发出了一系列软、硬件系统和算法、模型，并且这些系统和算法、模型主要基于双目视觉测量原理，双目视觉测量精度受到相机性能、光照和基线长度的影响，算法程序复杂，设备昂贵，应用限制偏多。随着传感器性能大幅度提升，从单目图像中测量物体尺寸成为可能，王忠亮等基于单目机器视觉研究定位高压输电线路障碍物，但单目测量的方法和模型多数是针对摄像机等专业拍摄仪器和计算机处理平台。目前多株立木高度测量一般基于航空遥感影像进行模型构建，实现林分立木高度测量。近年来随着智能移动设备快速发展，许多学者基于智能移动设备和Android平台开发了测树软件，李亚东等研究基于Android手机传感器的实现林木单株树高测量的技术，周玉晨等设计与研究基于Android的角规测树及数据处理软件，周克俞等基于Android平台设计实现立木高度和胸径测量的系统。这些软件和方法在测量树高时基于三角函数原理，在一定程度上实现了树高因子的实时快速测量，但是操作复杂，鲁棒性差。

综上分析，目前基于智能移动设备的树高测量方法的研究多集中于单株立木高度测量，以图像中多株待测立木高度测量研究较少。

因此，需要对现有技术进行改进。

发明内容

本发明要解决的技术问题是提供一种高效的基于智能手机的多株立木高度测量方法。

为解决上述技术问题，本发明提供一种基于智能手机的多株立木高度测量方法，包括以下步骤；

1)、在设标立木上设置矩形标定物，执行步骤2)；

2)、图像采集，获得设标立木数字图像，执行步骤3)；

3)、拍摄参数获取，得到手机相机拍摄高度h和目标物成像角度α，执行步骤6)和步骤16)；

4)、智能手机标定，执行步骤5)；

5)、智能手机的相机参数获取，采用张正友标定法，获取手机相机所述的内部参数、外部参数、相机镜头畸变参数、图像分辨率vmax、umax和手机相机焦距f，

相机镜头畸变参数中k₁、k₂为径向畸变参数；P₁、P₂为切向畸变参数；S₁、S₂为薄棱镜畸变参数；

执行步骤6)；

6)、根据相机参数，对步骤2)采集的设标立木数字图像进行非线性畸变校正，得到校正图像Pi’(ui’,vi’)，执行步骤7)和步骤10)；

7)、对校正图像进行图像增强和去噪处理，得到增强去噪图像，执行步骤8)；

8)、对增强去噪图像进行轮廓提取和边缘检测处理，得到处理后图像，执行步骤9)；

9)、对处理后图像进行待测立木最值点识别，得到最高点t1和最低点t2，执行步骤14)、步骤17)和步骤18)；

10)、对校正图像进行图像区域分割，得到标定物图像，执行步骤11)；

11)、对标定物图像进行二值化处理和增强处理，得到二值图像，执行步骤12)；

12)、对二值图像进行标定物特征点获取，得到矩形的4个顶点的像在像素坐标系中的坐标值，执行步骤13)；

13)、灭点算法计算，得到设标立木所在平面即参考平面的方程，执行步骤14)；

14)、设标立木高度测量模型建立，执行步骤15)；

设最高点t1和最低点t2在相机坐标系中的坐标为(xt1,yt1,zt1)和(xd2,yd2,zd2)，在图像平面的坐标为t1’和t2’，在像素坐标系中的坐标为(ut1,vt1)和(ut2,vt2)；计算获取设标立木高度如式(17)所式：

15)、将最高点t1和最低点t2代入式(17)计算设标立木高度，执行步骤18)；

16)、立木深度信息提取模型建立，执行步骤17)；

设目标物成像角度α、纵坐标像素值v，得立木深度信息提取模型为：

β为相机内部重力传感器获取的倾斜角，fy为像素单位下相机的焦距，θ为相机垂直视场角的一半，V_max为图像垂直和水平方向的最大像素值，V₀为相机内部参数图像中心点像素值，δ为相机非线性畸变项误差；h为手机相机拍摄高度；

17)、设标立木深度信息，执行步骤18)；

将步骤9中获取待测立木的最值点中的最低点t2的纵坐标像素值，代入公式11中，得到的D作为立木深度；

18)、多株立木高度测量模型建立，执行步骤19)；

设立待测立木TreeA和待测立木TreeB，待测立木TreeA距离拍摄手机的前后距离距离较近，待测立木TreeB距离拍摄手机的前后距离距离较远；以设标立木作为待测立木TreeA或待测立木TreeB，另一株待测立木作为余下的待测立木TreeA或待测立木TreeB；

设待测立木TreeA水平方向在光轴形成虚拟立木A；同理待测立木TreeB和虚拟立木B，虚拟立木A、虚拟立木B和相机在同一平面；得多株立木高度测量模型：

Oa为立木TreeA树木位置根部，O_ABO_b为将立木TreeA轮廓最高点和立木TreeB轮廓最高点代入设标立木模型获取距离，O_AB为将图像中立木TreeA轮廓最高点投影到设标立木TreeB参考平面上的点；

Ob为立木TreeB树木位置根部，O_abO_a为将立木TreeA轮廓最高点和立木TreeB轮廓最高点代入设标立木模型获取距离，O_ab为将图像中立木TreeB轮廓最高点投影到设标立木TreeA参考平面上的点；

D_A和D_B为立木TreeA和立木TreeB的深度距离，hc为手机距离地面实际物理距离；O_BO_b为待测立木TreeB高度，O_AO_a为待测立木TreeA高度；

19)、将设标立木高度代入多株立木高度测量模型，得到待测立木高度。

作为对本发明基于智能手机的多株立木高度测量方法的改进，步骤13)包括：

建立空间坐标系(xw,yw,zw)和相机坐标系(xc,yc,zc)，zc为光轴方向，得到空间坐标系与相机坐标系重合；

T＝[0 0 0]^T (2)

手机相机采集立木图像，建立坐标系；图像坐标系(x,y)原点O位于图像中心，且在像素坐标系(u,v)；成像平面关于光轴zc的交点为O，相机焦距f＝OOc，图像坐标系和相机坐标系两者之间的关系为：

式(3)中：坐标点(u,v)为像素坐标系中的值；u0、v0为相机内部参数；dx、dy为单个像素点在图像坐标系中对应的x、y轴的物理尺寸；

将标定物所在平面w设为参考平面，在相机坐标系中的点P坐标为(xc,yc,zc)，在图像坐标系中P’(x,y)，两者之间满足如下关系：

式(4)中：K为相机内部参数矩阵；fx，fy为x和y方向上对应的焦距；

设标定物图像4个定点分别为P1、P2、P3、P4，在成像平面中形成图像坐标点分别为PI1、PI2、PI3、PI4，因手机相机采集图像产生透视畸变，PI1、PI2与PI3、PI4不平行且在延长线上相交记为点M1，同理可得PI1、PI4与PI2、PI3延长线上交点记为M2；基于小孔成像原理和光学成像原理，灭点和相机光心两点之间的连线平行于该灭点形成的空间平行线，即ocM1║P1P2，ocM2║P2P3；

在像素坐标系中PI1坐标点为(u1,v1)，同理PI2(u2,v2)、PI3(u3,v3)和PI4(u4,v4)，计算的灭点坐标为M1(um1,vm1)、M2(um2,vm2)，在相机坐标系中灭点坐标为M1(xm1,ym1,zm1)、M2(xm2,ym2,zm2)；通过归一化处理，在相机坐标系中计算两个灭点对应的空间连线ocM1和ocM2对应的单位方向向量为V1和V2，由式(5)可得标定参考平面在相机中的单位法向量Vp；

V_p＝V₁×V₂ (5)

根据矩形标定位的长L和宽H，利用最小二乘法原理，将式(6)、(7)联立，计算获取坐标点P1(x1,y2,z1)和其他3个坐标点P2、P3、P4的值；

根据单位法向量Vp和坐标点P1、P2、P3、P4联立，得定参考平面在相机坐标系中的平面方程：

Ax+By+Cz+D＝0 (8)

式(8)中：A、B、C、D为平面方程的特征系数，根据上述式(5)、(6)和(7)得出；

作为对本发明基于智能手机的多株立木高度测量方法的进一步改进，步骤16)包括：

当α取最小值α＝αmin＝90-θ-β时，θ为相机垂直视场角的一半，即被拍摄物体投影到图片最底端时，v＝vmax，代入式(17)可得：

90-β-θ＝a·v_max+b (18)

当αmin+2θ>90°时，v无限趋近于v0-tanβ*fy，fy为像素单位下相机的焦距，代入式(17)可得：

90＝a·(v₀-tanβ·f_y)+b (19)

当αmin+2θ<90°时，v＝0，代入式(17)可得：

90-β+θ＝b (20)

一半的相机垂直视场角θ的正切值等于相机CMOS或CCD图像传感器边长的一半除以相机焦距，故可以计算出θ：

公式(21)中L_CMOS为相机CMOS或CCD图像传感器的边长，结合式(18)～(21)，F(α,β)为：

公式(10)中δ为相机非线性畸变项误差，结合手机相机拍摄高度h，根据三角函数原理建立手机相机深度提取模型：

根据相机旋转角度β与一半的相机垂直视场角θ之间的大小关系，选择对应的深度提取模型，将步骤13)求算的相机内部参数图像中心点像素值v0，y轴上归一化的焦距fy以及图像分辨率vmax，和步骤16)求算的待测立木最低点纵坐标像素值v、相机旋转角度β和手机相机拍摄高度h代入所述深度提取模型，计算待测立木目标点深度值D。

本发明基于智能手机的多株立木高度测量方法的技术优势为：

本发明在分析现有立木高度测量方法的优点和不足的基础上，以快速、实时、自动获取立木高度为目标，提出一种基于智能移动设备的单目视觉系统的多株立木高度测量方法。首先将拍摄手机相机通过棋盘格角点检测方法和张正友标定进行非线性畸变校正，并获取相机参数，并用标定手机相机采集单幅的多株立木图像，结合图像处理技术进行立木图像特征点提取；然后通过灭点算法和图像处理技术计算设置标定物的立木(简称为设标立木)的高度，利用待测立木成像与纵坐标像素、相机拍照角度之间的映射关系建立多株立木深度信息提取模型，将图像中待测立木最低点代入深度信息提取模型得到待测立木深度信息距离；最后结合设标立木高度和获取的待测多株立木深度信息建立多株立木测量模型，将图像中多株立木垂直最值点代入模型，实现多株立木高度测量。

附图说明

下面结合附图对本发明的具体实施方式作进一步详细说明。

图1为本发明基于智能手机的多株立木高度测量方法的流程示意图；

图2为本发明标定物灰度化示意图；

图3为本发明标定物二值化示意图；

图4为本发明特征点识别示意图；

图5为本发明原图；

图6为本发明噪声处理示意图；

图7为本发明立木树冠处理示意图；

图8为本发明立木底端处理示意图；

图9为本发明边缘检测示意图；

图10为本发明灭点形成原理示意图；

图11为本发明待测立木和立木成像在光轴方向虚拟成像图；

图12为本发明待测立木和成像平面在二维世界关系图。

具体实施方式

下面结合具体实施例对本发明进行进一步描述，但本发明的保护范围并不仅限于此。

实施例1、基于智能手机的多株立木高度测量方法，如图1-12所示，

立木高度测量时主要是利用智能终端的移动设备的图像采集功能和设备内置重力传感器的角度测量功能，在进行立木高度高度测量时，首先需要确定测量点的合适位置，确保待测多株立木整株在视线范围内无遮挡，选取一株待测立木放置标定物(设为设标立木)，并用智能移动设备手机采集立木图片。将手机相机拍摄棋盘标定板进行标定和相机参数获取，利用非线性畸变模型对采集立木图片进行非线性畸变校正；对校正图像利用显著性分析图调整HSV空间特征分量和RGB颜色空间分布增加颜色的对比度和图像亮度，进行图像灰度化和二值化处理，锐化图像边缘，通过开、闭算法对图像进行孔洞填充，然后使用去噪算法和Canny算子进行图像轮廓提取和边缘检测，识别立木垂直最值点；利用图像分割处理，对标定物进行识别，通过灰度化和二值化处理，识别标定物特征点；利用标定物特征点在图像中的像素位置，通过灭点算法，建立设标立木平面坐标系和设标立木高度计算模型，代入设标立木垂直最值点，获取设标立木高度；分析图像中立木图像像素、物体成像角度和相机旋转角度，建立立木深度信息提取模型，代入立木最低点，获取待测立木深度信息；分析待测多株立木高度和立木深度信息之间关系，建立多株立木高度测量模型，将待测多株立木垂直最值点代入，进行多株立木高度测量。

包括以下步骤：

1)、设标立木设置，执行步骤2)；

制作矩形标定物，选取长为50cm和高为65cm的白色纸板，中间位置张贴黑色纸张长为45cm和高为60cm(大一些也可以)；将矩形标定物布置在相机视野范围内的任意一株待测立木平面内(简称为参考平面)。该步骤中矩形标定物只要挂在待测立木上即可。

2)、图像采集，获得设标立木数字图像，执行步骤3)；

利用智能手机获得待测立木的数字图像。(包括设标立木所在的所有待测立木，设标立木只是其中一株)。

4)、拍摄参数获取，执行步骤6)和步骤16)；

利用皮尺测量用三脚架放置手机相机进行拍摄图像时手机相机拍摄高度h，利用Android语言内嵌系统调用手机重力传感器获取智能手机镜头拍摄倾斜角度(目标物成像角度α)。

4)、智能手机标定，执行步骤5)；

智能手机拍摄棋盘标定板进行标定，用智能手机在不同方向对棋盘标定板进行拍摄获取20-30张图片；

5)、智能手机的相机参数获取，执行步骤6)；

采用张正友标定法，通过Java结合Opencv进行相机标定和相机参数获取，获取相机内参、外部参数、相机镜头畸变参数(k₁、k₂为径向畸变参数；P₁、P₂为切向畸变参数；S₁、S₂为薄棱镜畸变参数)、图像分辨率vmax、umax(图像垂直和水平方向的最大像素值)和手机相机焦距f。

6)、根据相机参数，对步骤2)采集的设标立木数字图像进行非线性畸变校正，得到校正图像，执行步骤7)和步骤10)；

立木图像特征点在世界坐标系中的坐标Pi＝(Xi,Yi,Zi)，其在图像中像素点的坐标表示为Pi，根据已知立木图像空间坐标系、相机内参M与目标物体相机坐标系的位置关系得到[R|t]，可以求得特征点在图像像素坐标系中的相应的图像点坐标值。根据光学成像原理和智能手机相机针孔模型可得：

式(6)中：M为相机内参矩阵；t为世界坐标系和相机坐标系之间的平移矩阵；R为世界坐标系和相机坐标系之间的旋转矩阵；αx＝f/dx，图像坐标系下u轴归一化焦距，αy＝f/dy，图像坐标系下v轴归一化焦距；dx，dy分别表示传感器u轴和v轴上单位像素的尺寸大小。

由于实际的工艺制作问题，制作出一个数学上理想透镜不大可能，并且从机械方面将透镜与成像仪不能完全平行安装，导致实际镜头的光学成像模型不是理想的小孔成像模型，因此空间图像点实际坐标值偏离小孔成像模型计算出来的坐标值。假设实际图像点坐标值为Pi’(ui’,vi’)，理论计算值为Pi(ui,vi)，两者之间的关系表示为：

式(7)中：Δu＝a_xx_nd，Δv＝a_yy_nd为非线性畸变值，u_i和v_i为理论图像坐标值，u_i’和v_i’为实际图像坐标值。

式(8)中：其中k₁、k₂为径向畸变参数；P₁、P₂为切向畸变参数；S₁、S₂为薄棱镜畸变参数。x_nd,y_nd为坐标点具有的非线性畸变值。r为相机内参(手机相机所述的内部参数)。

由公式(7)得到校正图像。

7)、对校正图像进行图像增强和去噪处理(现有技术)，得到增强去噪图像，执行步骤8)；

在手机相机采集的图像过程中，由于受到外界干扰，导致在采集的图像进行变换处理的过程中图像产生局部噪声，包括产生气泡、毛刺与凹陷等，这些噪声的产生对后续图像中立木轮廓识别提取产生影响。为了提高图像中立木垂直最值点识别的准确性，需对校正图像进行补偿处理。

将校正图像根据图像显著性分析法，利用分析图调整HSV空间特征分量和RGB颜色空间分布增加颜色的对比度和图像亮度。

图像进行增强处理之后，根据RGB色彩模式定义，按红色、绿色、蓝色3个颜色通道进行分离，分别得到单通道图像IR、IG、IB，将图像中像素点的3个单颜色通道灰度矩阵为[gR gG gB]，其中gR、gG、gB分别为R、G、B通道的灰度值。为了增强立木图像，通过高斯滤波的方式突出图像中低频部分，降低图像中高频部分，提高图像的信噪比，降低图像中的噪声。将去噪图像进行二值化处理，然后进行膨胀、腐蚀操作，去除孔、洞。得到增强去噪图像。

8)、对增强去噪图像进行轮廓提取和边缘检测处理，得到处理后图像，执行步骤9)；

图像部分处理后进行立木边缘检测和立木轮廓边缘点识别，本发明采用Canny算子进行检测运算。Canny算子相较于其他算子有良好的抗噪性，能够产生边缘梯度强度和方向两个信息。这样判断降低非边缘点判为边缘点的概率和将边缘点判为非边缘点的概率，即检测出的边缘点要尽可能在实际边缘的中心，并且单一边缘具有唯一响应和最大程度抑制对虚假边缘的响应。对边缘点定义获取的梯度值进行“非极大值抑制”，进行双阐值化和边缘连接处理。通过图像轮廓提取和边缘检测处理。

9)、对处理后图像进行待测立木最值点识别，执行步骤14)、步骤17)和步骤18)；

对处理后图像的图像轮廓进行检测，识别处理后图像中所有待测立木边缘的y轴方向最大值TmaxY和y轴方向最小值TminY(最大值TmaxY即为最高点t1，y轴方向最小值TminY即为最低点t2)，获取待测立木垂直最值点。

10)、对校正图像进行图像区域分割(现有技术)，得到标定物图像，执行步骤11)；

利用图像分割算法和边缘检测算法，进行标定物图像获取。

利用Grabcut算法方法，人为指定一个包含矩形标定物前景的矩形，随后用基于图切算法在图像中提取前景，OpenCV实现了这个算法，可以直接调用。

11)、对标定物图像进行二值化处理和增强处理，得到二值图像，执行步骤12)；

增强处理同步骤7。

12)、对二值图像进行标定物特征点获取，执行步骤13)；

利用hough变换原理计算矩形边界在图像坐标系中的直线方程。进而可以计算出矩形的4个顶点的像在像素坐标系中的坐标值。

13)、灭点算法计算，得到设标立木所在平面即参考平面的方程，执行步骤14)；

建立空间坐标系(xw,yw,zw)和相机坐标系(xc,yc,zc)，zc为光轴方向，得到空间坐标系与相机坐标系重合。根据相机标定得到相机内部参数、外部参数，其中R为旋转矩阵，T为平移矩阵。

T＝[0 0 0]^T (2)

手机相机采集立木图像，建立坐标系。图像坐标系(x,y)原点O位于图像中心，且在像素坐标系(u,v)。成像平面关于光轴zc的交点为O，相机焦距f＝OOc，图像坐标系和相机坐标系两者之间的关系为：

式(3)中：坐标点(u,v)为像素坐标系中的值；u0、v0为相机内部参数；dx、dy为单个像素点在图像坐标系中对应的x、y轴的物理尺寸。

将标定物所在平面w设为参考平面，在相机坐标系中的点P坐标为(xc,yc,zc)，在图像坐标系中P’(x,y)，两者之间满足如下关系：

式(4)中：K为相机内部参数矩阵；fx，fy为x和y方向上对应的焦距。

设标定物图像4个定点分别为P1、P2、P3、P4，在成像平面中形成图像坐标点分别为PI1、PI2、PI3、PI4，如图2所示。因手机相机采集图像产生透视畸变，PI1、PI2与PI3、PI4不平行且在延长线上相交记为点M1，同理可得PI1、PI4与PI2、PI3延长线上交点记为M2。基于小孔成像原理和光学成像原理，灭点和相机光心两点之间的连线平行于该灭点形成的空间平行线，即ocM1║P1P2，ocM2║P2P3。

在像素坐标系中PI1坐标点为(u1,v1)，同理PI2(u2,v2)、PI3(u3,v3)和PI4(u4,v4)，计算的灭点坐标为M1(um1,vm1)、M2(um2,vm2)，在相机坐标系中灭点坐标为M1(xm1,ym1,zm1)、M2(xm2,ym2,zm2)。通过归一化处理，在相机坐标系中计算两个灭点对应的空间连线ocM1和ocM2对应的单位方向向量为V1和V2。因ocM1║P1P2，ocM2║P2P3，在相机坐标系中标定物相邻两边的方向向量为V1和V2，由式(5)可得标定参考平面在相机中的单位法向量Vp。

V_p＝V₁×V₂ (5)

已知标定物P1、P2、P3、P4的长L和宽H，且在相机坐标系中坐标点P1和P2、P3、P4的关系如式(6)所示，且4个点在相机坐标系和像素坐标系中存在如在关系，如式(7)所示。利用最小二乘法原理，将式(6)、(7)联立，计算获取坐标点P1(x1,y2,z1)和其他3个坐标点P2、P3、P4的值。

上面已经计算出参考平面上其中一点P1在相机坐标系中的坐标，以及参考平面在相机坐标系中的单位法线向量如式(5)。因此，可以确定参考平面在相机坐标系中的平面方程(8)

Ax+By+Cz+D＝0 (8)

式(8)中：A、B、C、D为平面方程的特征系数，根据上述式(5)、(6)和(7)可计算得出A、B、C、D。

14)、设标立木高度测量模型建立，执行步骤15)；

获取设标立木轮廓最值点为最高点t1和最低点t2，最高点t1和最低点t2在相机坐标系中的坐标为(xt1,yt1,zt1)和(xd2,yd2,zd2)，在图像平面的坐标为t1’和t2’，在像素坐标系中的坐标为(ut1,vt1)和(ut2,vt2)。利用灭点在相机坐标系中的坐标，计算获取直线oct1和直线oct2在相机坐标系中的直线方程(即相机坐标系中点t1、t2对应的极线方程)，如式(16)所示：

将式(16)中的直线oct1和直线oct2方程与标定参考平面方程式(15)联立，得到设标立木轮廓垂直最值点t1、t2在相机坐标系中坐标值(xt1,yt1,zt1)和(xt2,yt2,zt2)，计算获取设标立木高度如式(17)所式：

15)、计算设标立木高度，执行步骤18)；

将设置标定的待测立木，根据步骤6中获取设标立木轮廓最高点和最低点的纵坐标像素，代入t1、t2中，得到设标立木高度。

16)、立木深度信息提取模型建立，执行步骤17)；

根据目标物成像角度α与纵坐标像素值v之间的线性关系设定抽象函数，建立含目标物成像角度α、纵坐标像素值v和相机旋转角β三个参数空间关系模型，即α＝F(v,β)，

不同型号的设备和相机旋转角度下，被拍摄物体纵坐标像素值与成像角度均呈极显著负线性相关关系，且该线性关系的斜率与截距有所不同，故设：

α＝F(v,β)＝a·v+b (17)

其中参数a、b均与相机型号和相机旋转角度有关(人为设置)；

当α取最小值α＝αmin＝90-θ-β时，θ为相机垂直视场角的一半，即被拍摄物体投影到图片最底端时，v＝vmax(为vmax图像垂直和水平方向的最大像素值)，代入式(17)可得：

90-β-θ＝a·v_max+b (18)

当αmin+2θ>90°，即θ>β时，此时相机上视角高于水平线，地平面无限远处，α无限接近于90°，此时v无限趋近于v0-tanβ*fy，fy为像素单位下相机的焦距，β为负值即相机逆时针旋转时亦同理，因此，代入式(17)可得：

90＝a·(v₀-tanβ·f_y)+b (19)

当αmin+2θ<90°，即θ<β时，此时相机上视角低于水平线，地平面无限远处目标物成像角度α取最大值，αmax＝αmin+2θ＝90-β+θ时，即被拍摄物体投影到图片最高点时，v＝0，代入式(17)可得：

90-β+θ＝b (20)

当θ＝β时拍摄不到相应的立木图像，不做考虑。

根据针孔相机构造原理，一半的相机垂直视场角θ的正切值等于相机CMOS或CCD图像传感器边长的一半除以相机焦距，故可以计算出θ：

公式(21)中L_CMOS为相机CMOS或CCD图像传感器的边长，结合式(18)～(21)，F(α,β)为：

公式(10)中δ为相机非线性畸变项误差，结合手机相机拍摄高度h，根据三角函数原理建立手机相机深度提取模型：

根据相机旋转角度β与一半的相机垂直视场角θ之间的大小关系，选择对应的深度提取模型，将步骤13)求算的相机内部参数图像中心点像素值v0，y轴上归一化的焦距fy以及图像分辨率vmax，和步骤16)求算的待测立木最低点纵坐标像素值v、相机旋转角度β和手机相机拍摄高度h代入所述深度提取模型，计算设标立木目标点深度值D。

17)、设标立木深度信息，执行步骤18)；

将步骤9中获取待测立木的最值点中的最低点的纵坐标像素值，代入公式11中。拍摄立木图像时，立木高度均大于拍摄时的相机高度，因此只取θ＞β时的情况即可，得到设标立木深度信息。

18)、多株立木高度测量模型建立，执行步骤19)；

设立待测立木TreeA和待测立木TreeB(待测立木TreeA和待测立木TreeB与拍摄手机的远近不同即可)，待测立木TreeA距离拍摄手机的前后距离距离较近，待测立木TreeB距离拍摄手机的前后距离距离较远；以设标立木作为待测立木TreeA或待测立木TreeB，另一株待测立木作为余下的待测立木TreeA或待测立木TreeB；

虚拟立木A、虚拟立木B和相机在同一平面。

在单目视觉多株立木高度测量中，基于地面平整且没有坡度的假设。测量人员通过手机相机进行图像信息采集，为获取整株立木图像需进行仰角拍摄使得成像面不与地面垂直，产生手机倾斜角β(倾斜角β由相机内部重力传感器获取)，手机距离地面实际物理距离为hc，平面γ为成像平面关于相机光心的对称平面(简称图像平面γ)。因此，在相机光学畸变校正的情况下，待测立木TreeA水平方向在光轴形成虚拟立木A。虚拟立木A和待测立木TreeA在成像面中形成像的大小完全一致，且两者之间呈平移关系，同理待测立木TreeB和虚拟立木B，如图11所示。

基于现有待测立木图像处理和特征点识别，结合小孔成像原理，将待测立木和立木成像平移到光轴方向形成虚拟像，此时虚拟的待测立木和成像平面在二维世界之间关系，如图12所示。根据光学成像原理可知，虚拟立木OAOa平移虚拟投影在相机成像面中成像为uAua，且最高点在相机成像面中坐标点为uA(xA,yA)。立木TreeB平移虚拟投影在相机成像面中成像为uBub，且最高点在相机成像面中坐标点为uB(xB,yB)，在像uBub中可以找到点uAB与坐标点uA完全重合。将点uA(xA,yA)代入设标立木高度测量模型可以计算高度O_ABO_b，同理在在像uAua在y轴延长方向中可以找到点uab与坐标点uB完全重合，将uB(xB,yB)代入立木高度测量模型可以计算高度O_abO_a。

设立木TreeB为设标立木，平移虚拟投影O_BO_b，已知设标立木实际高度HtB＝O_BO_b，可知所求待测多株立木距离手机相机的深度距离，即DA＝CmOa和DB＝CmOb。基于相似三角形原理，已知相机拍摄高度hc，可以计算得到立木TreeA虚拟投影高度O_AO_a(即立木TreeA测量高度)，如式(24)、(25)所示：

Oa为立木TreeA树木位置根部，O_ABO_b为将立木TreeA轮廓最高点和立木TreeB轮廓最高点代入设标立木模型获取距离，O_AB为将图像中立木TreeA轮廓最高点投影到设标立木TreeB参考平面上的点。

Ob为立木TreeB树木位置根部，O_abO_a为将立木TreeA轮廓最高点和立木TreeB轮廓最高点代入设标立木模型获取距离，O_ab为将图像中立木TreeB轮廓最高点投影到设标立木TreeA参考平面上的点。

D_A和D_B为立木TreeA和立木TreeB的深度距离，hc为手机距离地面实际物理距离；O_BO_b为待测立木TreeB高度，O_AO_a为待测立木TreeA高度；

设立木TreeA平移虚拟投影O_AO_a且为设标立木，已知设标立木实际高度HtA＝O_AO_a，可知所求多株立木距离手机相机的深度距离，即DA＝CmOa和DB＝CmOb。基于相似三角形原理，已知相机拍摄高度hc，可以计算得到立木TreeB虚拟投影高度O_BO_b(即立木TreeB测量高度)，如式(26)、(27)所示：

因此，结合式(24)～(27)，可得多株立木高度测量模型式(28)：

19)、立木最高点和最低点代入多株立木高度测量模型式，得到待测多株立木高度。

当设置的设标立木获取为远距离的立木时(根据立木距离拍摄手机的前后距离远近，人为判断)：

立木TreeB为设标立木时，根据设标立木模型公式，将立木TreeB的最高点和最低点代入公式，得到立木TreeB的高度。并且需要将待测立木中的立木(除了设标立木以外的立木)(举例：立木TreeA)的最高点和设标立木的最点像素代入设标立木模型公式，得到设标立木参考平面计算高度O_ABO_b。在已知立木TreeA和立木TreeB的深度距离为DA和DB，拍摄参数中相机高度hc，代入公式(25)，得到待测立木中的立木(立木TreeA)。

当设置的设标立木获取为近距离的立木时：

立木TreeA为设标立木时，根据设标立木模型公式，将立木TreeA的最高点和最低点代入公式，得到立木TreeA的高度。并且需要将待测立木中的立木(除了设标立木以外的立木)

(举例：立木TreeB)的最高点和设标立木的最点像素代入设标立木模型公式，得到设标立木参考平面计算高度O_abO_a。在已知立木TreeA和立木TreeB的深度距离为DA和DB，拍摄参数中相机高度hc，代入公式(27)，得到待测立木中的立木(立木TreeB)。

实验：

基于上述原理和算法流程，为验证方法的可行性及精度，本文测试手机为小米手机，型号为MI 2S，Android版本为4.1.2。本研究在Android系统平台上，利用Java语言结合C++语言，按照上述方法编写程序并调试后，使用小米手机相机作为图片采集设备在自然环境下进行精度验证。相机拍摄高度为测量拍摄所用三脚架高度，经相机标定计算可得小米3手机相机内部参数为：fx＝3486.5637，u0＝1569.0383，fy＝3497.4652，v0＝2107.98988，图像分辨率为3120×4208。

为了验证本文方法在真实场景应用中树高非接触性测量的精度和实时性，实验在校园林区内拍摄5幅图片，每张图像中包含2-3株立木，在已经识别的立木轮廓垂直最值点的基础上，验证立木高度测量精度。样本立木进行树高测量真值为待测多株立木通过长直杆和皮尺、布鲁莱斯测高器进行多次测量测得求平均获得。实验测量结果如表1所示。

本发明立木高度测量在距离为3～10m范围内平均相对误差为2.91％，树高测量精度为95.13％，相对误差为不超过5％，具有较高的测量精度和实时性，可以满足基于移动端的立木树高自动测量系统和森林资源连续清查要求。

表1待测立木树高测量精度

结论：

本发明提出一种基于智能手机的设标立木和深度信息的多株立木高度测量方法。该方法利用安卓手机在不同角度采集的棋盘格标定板，引用张正友标定法进行图像的非线性畸变校正。利用手机采集立木图像进行图像分割处理，获取的设标立木的标定物特征点在图像坐标系中的位置，根据灭点原理和相机内部参数，建立参考平面在相机坐标系中的平面方程，结合立木识别垂直最值点建立设标立木高度测量模型，并进行设标立木高度获取。本文根据相机成像原理选取了特殊成像点的实际成像角度和纵坐标像素值代入线性关系函数，通过相机标定获取相机内部参数，建立相应的立木深度信息提取模型，计算出待测立木的深度信息。然后结合投影几何模型和相机立体成像系统原理，分析设标立木与待测立木深度信息之间关系，计算待测多株立木垂直最值点虚拟投影像的物理高度，并根据勾股定理实现基于单目视觉系统的多株立木高度的非接触测量。为验证模型精度，本文真实场景中采集立木图像进行验证，树高测量相对误差小于5％。与同类文献相比，该方法还可以获取图像中多株待测立木高度。同时，该方法引入了立木深度信息提取模型，计算单幅图片上待测立木最低点到拍摄相机的深度信息，实现了多株立木高度的自动测量。

最后，还需要注意的是，以上列举的仅是本发明的若干个具体实施例。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (3)

1.基于智能手机的多株立木高度测量方法，其特征在于，包括以下步骤；

1)、在设标立木上设置矩形标定物，执行步骤2)；

2)、图像采集，获得设标立木数字图像，执行步骤3)；

3)、拍摄参数获取，得到手机相机拍摄高度h和目标物成像角度α，执行步骤6)和步骤16)；

4)、智能手机标定，执行步骤5)；

5)、智能手机的相机参数获取，采用张正友标定法，获取手机相机所述的内部参数、外部参数、相机镜头畸变参数、图像分辨率vmax、umax和手机相机焦距f，

相机镜头畸变参数中k₁、k₂为径向畸变参数；P₁、P₂为切向畸变参数；S₁、S₂为薄棱镜畸变参数；

执行步骤6)；

6)、根据相机参数，对步骤2)采集的设标立木数字图像进行非线性畸变校正，得到校正图像Pi’(ui’,vi’)，执行步骤7)和步骤10)；

7)、对校正图像进行图像增强和去噪处理，得到增强去噪图像，执行步骤8)；

8)、对增强去噪图像进行轮廓提取和边缘检测处理，得到处理后图像，执行步骤9)；

9)、对处理后图像进行待测立木最值点识别，得到最高点t1和最低点t2，执行步骤14)、步骤17)和步骤18)；

10)、对校正图像进行图像区域分割，得到标定物图像，执行步骤11)；

11)、对标定物图像进行二值化处理和增强处理，得到二值图像，执行步骤12)；

12)、对二值图像进行标定物特征点获取，得到矩形的4个顶点的像在像素坐标系中的坐标值，执行步骤13)；

13)、灭点算法计算，得到设标立木所在平面即参考平面的方程，执行步骤14)；

14)、设标立木高度测量模型建立，执行步骤15)；

设最高点t1和最低点t2在相机坐标系中的坐标为(xt1,yt1,zt1)和(xd2,yd2,zd2)，在图像平面的坐标为t1’和t2’，在像素坐标系中的坐标为(ut1,vt1)和(ut2,vt2)；计算获取设标立木高度如式(17)所式：

15)、将最高点t1和最低点t2代入式(17)计算设标立木高度，执行步骤18)；

16)、立木深度信息提取模型建立，执行步骤17)；

设目标物成像角度α、纵坐标像素值v，得立木深度信息提取模型为：

β为相机内部重力传感器获取的倾斜角，fy为像素单位下相机的焦距，θ为相机垂直视场角的一半，V_max为图像垂直和水平方向的最大像素值，V₀为相机内部参数图像中心点像素值，δ为相机非线性畸变项误差；h为手机相机拍摄高度；

17)、设标立木深度信息，执行步骤18)；

将步骤9中获取待测立木的最值点中的最低点t2的纵坐标像素值，代入公式11中，得到的D作为立木深度；

18)、多株立木高度测量模型建立，执行步骤19)；

设立待测立木TreeA和待测立木TreeB，待测立木TreeA距离拍摄手机的前后距离距离较近，待测立木TreeB距离拍摄手机的前后距离距离较远；以设标立木作为待测立木TreeA或待测立木TreeB，另一株待测立木作为余下的待测立木TreeA或待测立木TreeB；

设待测立木TreeA水平方向在光轴形成虚拟立木A；同理待测立木TreeB和虚拟立木B，虚拟立木A、虚拟立木B和相机在同一平面；得多株立木高度测量模型：

Oa为立木TreeA树木位置根部，O_ABO_b为将立木TreeA轮廓最高点和立木TreeB轮廓最高点代入设标立木模型获取距离，O_AB为将图像中立木TreeA轮廓最高点投影到设标立木TreeB参考平面上的点；

Ob为立木TreeB树木位置根部，O_abO_a为将立木TreeA轮廓最高点和立木TreeB轮廓最高点代入设标立木模型获取距离，O_ab为将图像中立木TreeB轮廓最高点投影到设标立木TreeA参考平面上的点；

D_A和D_B为立木TreeA和立木TreeB的深度距离，hc为手机距离地面实际物理距离；O_BO_b为待测立木TreeB高度，O_AO_a为待测立木TreeA高度；

19)、将设标立木高度代入多株立木高度测量模型，得到待测立木高度。

2.根据权利要求1所述的基于智能手机的多株立木高度测量方法，其特征在于，步骤13) 包括：

建立空间坐标系(xw,yw,zw)和相机坐标系(xc,yc,zc)，zc为光轴方向，得到空间坐标系与相机坐标系重合；

T＝[0 0 0]^T (2)

手机相机采集立木图像，建立坐标系；图像坐标系(x,y)原点O位于图像中心，且在像素坐标系(u,v)；成像平面关于光轴zc的交点为O，相机焦距f＝OOc，图像坐标系和相机坐标系两者之间的关系为：

式(3)中：坐标点(u,v)为像素坐标系中的值；u0、v0为相机内部参数；dx、dy为单个像素点在图像坐标系中对应的x、y轴的物理尺寸；

将标定物所在平面w设为参考平面，在相机坐标系中的点P坐标为(xc,yc,zc)，在图像坐标系中P’(x,y)，两者之间满足如下关系：

式(4)中：K为相机内部参数矩阵；fx，fy为x和y方向上对应的焦距；

设标定物图像4个定点分别为P1、P2、P3、P4，在成像平面中形成图像坐标点分别为PI1、PI2、PI3、PI4，因手机相机采集图像产生透视畸变，PI1、PI2与PI3、PI4不平行且在延长线上相交记为点M1，同理可得PI1、PI4与PI2、PI3延长线上交点记为M2；基于小孔成像原理和光学成像原理，灭点和相机光心两点之间的连线平行于该灭点形成的空间平行线，即ocM1║P1P2，ocM2║P2P3；

在像素坐标系中PI1坐标点为(u1,v1)，同理PI2(u2,v2)、PI3(u3,v3)和PI4(u4,v4)，计算的灭点坐标为M1(um1,vm1)、M2(um2,vm2)，在相机坐标系中灭点坐标为M1(xm1,ym1,zm1)、M2(xm2,ym2,zm2)；通过归一化处理，在相机坐标系中计算两个灭点对应的空间连线ocM1和ocM2对应的单位方向向量为V1和V2，由式(5)可得标定参考平面在相机中的单位法向量Vp；

V_p＝V₁×V₂ (5)

根据矩形标定位的长L和宽H，利用最小二乘法原理，将式(6)、(7)联立，计算获取坐标点P1(x1,y2,z1)和其他3个坐标点P2、P3、P4的值；

根据单位法向量Vp和坐标点P1、P2、P3、P4联立，得定参考平面在相机坐标系中的平面方程：

Ax+By+Cz+D＝0 (8)

式(8)中：A、B、C、D为平面方程的特征系数，根据上述式(5)、(6)和(7)得出。

3.根据权利要求2所述的基于智能手机的多株立木高度测量方法，其特征在于，步骤16)包括：

当α取最小值α＝αmin＝90-θ-β时，θ为相机垂直视场角的一半，即被拍摄物体投影到图片最底端时，v＝vmax，代入式(17)可得：

90-β-θ＝a·v_max+b (18)

当αmin+2θ>90°时，v无限趋近于v0-tanβ*fy，fy为像素单位下相机的焦距，代入式(17)可得：

90＝a·(v₀-tanβ·f_y)+b (19)

当αmin+2θ<90°时，v＝0，代入式(17)可得：

90-β+θ＝b (20)

一半的相机垂直视场角θ的正切值等于相机CMOS或CCD图像传感器边长的一半除以相机焦距，故可以计算出θ：

公式(21)中L_CMOS为相机CMOS或CCD图像传感器的边长，结合式(18)～(21)，F(α,β)为：

公式(10)中δ为相机非线性畸变项误差，结合手机相机拍摄高度h，根据三角函数原理建立手机相机深度提取模型：

根据相机旋转角度β与一半的相机垂直视场角θ之间的大小关系，选择对应的深度提取模型，将步骤13)求算的相机内部参数图像中心点像素值v0，y轴上归一化的焦距fy以及图像分辨率vmax，和步骤16)求算的待测立木最低点纵坐标像素值v、相机旋转角度β和手机相机拍摄高度h代入所述深度提取模型，计算待测立木目标点深度值D。