CN109269430B - 基于深度提取模型的多株立木胸径被动测量方法 - Google Patents

Abstract

本发明公开了一种基于深度提取模型的多株立木胸径被动测量方法，步骤一：对手机相机进行标定，获取相机内部参数和图像分辨率；步骤二：建立深度提取模型：

步骤三：采集并处理待测立木图像，计算摄影测量坐标系下立木深度值；步骤四：测量单幅图像中多株立木胸径值：

本发明的基于单目视觉和深度提取模型的多株立木胸径被动测量方法，能够适用于视场角、焦距、图像分辨率等参数不同的相机，提高测距精度和测量效率，为森林资源调查和林业资源管理提供支持。

Description

基于深度提取模型的多株立木胸径被动测量方法

技术领域

本发明涉及地面近景摄影测量领域，尤其是一种单目视觉***下基于单幅图像信息的多株立木胸径被动测量方法。

背景技术

胸高直径，亦称胸径(厘米)，是指树干距离地面根茎1.3米处树的直径^[1]。立木胸径是森林资源调查管理中一个重要的测量因子^[2-3]。传统的胸径测量方法测量值通常被视为胸径标准值^[4-5]，但方法耗时、耗力且效率低。随着传感技术、机器视觉及激光测距技术等的发展，非接触式测绘方法应运而生。非接触式胸径测量方法主要分为主动测量和被动测量两类^[6]。使用激光测距是主动测量的主流方法之一^[7-9]，利用激光雷达测距准确率高且数据处理速度快，但激光雷达测距具有***硬件集成成本高，移动端设备集成难度大等缺点^[10-11]。超声波技术依据声波传播技术及时间差计算立木距离并提取胸径，但超声波技术数据采集实时性差^[12-13]。机器视觉技术作为被动测量的主要途径，通过图像像素信息及相机成像原理进行距离估算并获取目标物尺寸，具有图像信息丰富、成本低等优点^[14-16]。机器视觉测量主要包括单目视觉和双目视觉测量两类^[17-19]。早期的图像信息提取方法主要是双目立体视觉原理和相机运动信息，需要多幅图像完成图像深度信息的提取^[20-22]。与双目视觉测量相比，单目测量图像采集不需要严格的硬件条件且设备集成便捷，因此近年来对单目视觉信息提取方法的研究逐渐深入^[23]。

采用单目视觉***测量图像中多株立木胸径需要先确定其深度距离信息。可以通过相机标定获取相机内外参数，结合投影模型求解出图像坐标系与世界坐标系之间的转换关系，从而计算目标物深度^[24-26]。此方法需要多次采集不同方位的标靶图像，并记录每个点在世界坐标系和图像坐标系中的对应坐标，标定对于测量精度影响较大。另外，由于移动端设备中集成了许多传感器，携带方便，因此许多学者基于移动端平台和单目视觉技术开发了一些立木因子测量平台，周克瑜等^[27]在测量胸径时基于三角函数原理，根据已知尺寸标定物进行对比计算，在一定程度上实现了胸径因子的快速测量，但是操作复杂，精度不高。黄晓东等^[16]利用数码相机获取待测样地上下方向的任意两张图片，并测量图像中任意一段物体长度用于推算摄影基线，计算真实空间比例关系，从而测量图像上样木胸径。该方法实现了图像中多株立木胸径的自动测量，但需要拍摄多张图片，且需预先测量图像中任意物体长度，该长度的精度对其它样木胸径测量值影响较大。

申请号为201710849961.3的发明申请，公开了一种改进的适用于智能移动端相机的相机标定模型及畸变矫正模型(以下简称：改进的带有非线性畸变项的标定模型)，该方法可以帮助矫正标定板图片，获取更高精度的相机内、外参数，不足的是，该方法没有扩展到对待测图像的非线性畸变校正及目标物的测量中。

参考文献：

[1]Library W P.Diameter at Breast Height[J].2016.

[2]史洁青，冯仲科，刘金成.基于无人机遥感影像的高精度森林资源调查***设计与实验[J].农业工程学报，2017，33(11)：82-90.

[3]Kunstler G，Falster D，Coomes DA，et al.Plant functional traits haveglobally consistent effects on competition[J].Nature，2016，529(7585)：204-207.

[4]Marsden C，Nouvellon Y，M’Bou A T，et al.Two independent estimationsof stand-level root respiration on clonal Eucalyptus stands in Congo：upscaling of direct measurements on roots versus the trenched-plot technique[J].New Phytologist，2008，177(3)：676-687.

[5]Saremi H，Kumar L，Turner R，et al.Impact of local slope and aspectassessed from LiDAR records on tree diameter in radiata pine(Pinus radiata，D.Don)plantations[J].Annals of Forest Science，2014，71(7)：771-780.

[6]刘文萍，仲亭玉，宋以宁.基于无人机图像分析的树木胸径预测[J].农业工程学报，2017，33(21)：99-104.

[7]Lin F，Dong X，Chen B M，et al.A Robust Real-Time Embedded VisionSystem on an Unmanned Rotorcraft for Ground Target Following[J].IEEE Trans onIndustrial Electronics，2012，59(2)：1038-1049.

[8]Kankare V，Puttonen E，Holopainen M，et a1.The effect of TLS pointcloud sampling on tree detection and diameter measurement accuracy[J].RemoteSensing Letters，2016，7(5)：495-502.

[9]孙俊灵，孙光民，马鹏阁，等.基于对称小波降噪及非对称高斯拟合的激光目标定位[J].中国激光，2017，44(6)：178-185.

[10]Liu C，Xing Y，Duanmu J，et al.Evaluating Different Methods forEstimating Diameter at Breast Height from Terrestrial Laser Scanning[J].Remote Sensing，2018，10(4)：513.

[11]Yang Xu，Zhang Yong，Xu Lu，et al.Increasing the range accuracy ofthree-dimensional ghost imaging ladar using optimum slicing number method[J].Chinese Physics B，2015，24(12)：319-324.

[12]Ahn S，Choi J，Doh N L，et al.A practical approach for EKF-SLAM inan indoor environment：fusing ultrasonic sensors and stereo camera[J].Autonomous Robots，2008，24(3)：315-335.

[13]Kumar S，Furuhashi H.Long-range measurement system usingultrasonic range sensor with high-power transmitter array in air[J].Ultrasonics，2017，74：186-195.

[14]石杰，李银伢，戚国庆，等.不完全测量下基于机器视觉的被动跟踪算法[J].华中科技大学学报，2017，45(6)：33-37.

[15]徐诚，黄大庆，孔繁锵.一种小型无人机无源目标定位及精度分析[J].仪器仪表学报，2015,36(5)：1115-1122.

[16]黄晓东，冯仲科.基于数码相机的样木胸径获取方法[J].农业机械学报，2015，46(9)：266-272.

[17]Zhang D.Monocular vision for pose estimation in space based oncone projection[J].Optical Engineering，2017，56(10)：1.

[18]Szeliski R.Computer vision：algorithms and applications[J].2010，21(8)：2601-2605.

[19]Sun W，Chen L，Hu B，et al.Binocular vision-based positiondetermination algorithm and system[C]//Proceedings of the 2012InternationalConference on Computer Distributed Control and Intelligent EnvironmentalMonitoring.Piscataway：IEEE Computer Society，2012：170-173.

[20]Ming A，Wu T，Ma J，et al.Monocular Depth-Ordering Reasoning withOcclusion Edge Detection and Couple Layers Inference[J].IEEE IntelligentSystems，2016，31(2)：54-65.

[21]Alexander E，Guo Q，Koppal S，et al.Focal Flow：Velocity and Depthfrom Differential Defocus Through Motion[J].International Journal of ComputerVision，2017(2)：1-22.

[22]Royden C S，Parsons D，Travatello J.The effect of monocular depthcues on the detection of moving objects by moving observers[J].VisionResearch，2016，124：7-14.

[23]李可宏，姜灵敏，龚永义.2维至3维图像/视频转换的图像深度提取方法综述[J].中国图象图形学报，2014，19(10)：1393-1406.

[24]Fathi H，Brilakis I.A Multi-Step Explicit Stereo CameraCalibration Approach to Improve Euclidean Accuracy of Large-Scale 3DReconstruction[J].Journal of Computing in Civil Engineering，2016，30(1)：04014120-04014120.

[25]Poulingirard A S，Thibault S，Laurendeau D.Influence of cameracalibration conditions on the accuracy of 3D reconstruction.[J].OpticsExpress，2016，24(3)：2678.

[26]黄小云，高峰，徐国艳，等.基于单幅立式标靶图像的单目深度信息提取[J].北京航空航天大学学报，2015，41(4)：649-655.

[27]周克瑜，汪云珍，李记，等.基于Android平台的测树***研究与实现[J].南京林业大学学报(自然科学版)，2016，40(4)：95-100.

[28]贺付亮，郭永彩，高潮，陈静.基于视觉显著性和脉冲耦合神经网络的成熟桑葚图像分割[J].农业工程学报，2017，33(6)：148-155.

[29]Harris C J.A combined corner and edge detector[J].Proc AlveyVision Conf，1988，1988(3)：147-151.

[30]Shi J，Tomasi.Good features to track[C]//Computer Vision andPattern Recognition.1994.Proceedings CVPR′94.1994IEEE Computer SocietyConference on.IEEE，2002：593-600.

[31]Geiger A，Moosmann F，Omer Car，et al.Automatic camera and rangesensor calibration using a single shot[C]//IEEE International Conference onRobotics and Automation.IEEE，2012：3936-3943.

[32]Zhang Z.Flexible Camera Calibration by Viewing a Plane fromUnknown Orientations[C]//The Proceedings of the Seventh IEEE InternationalConference on Computer Vision.IEEE.1999：666-673vol.1.

发明内容

本发明的目的是提供一种基于单目视觉和深度提取模型的多株立木胸径被动测量方法，能够适用于视场角、焦距、图像分辨率等参数不同的相机，提高测距精度和测量效率，为森林资源调查和林业资源管理提供支持。

为实现上述目的，本发明采用如下技术方案：

一种基于深度提取模型的多株立木胸径被动测量方法，包括如下步骤：

步骤一：对手机相机进行标定，获取相机内部参数和图像分辨率

采用张正友标定法，并引入改进的带有非线性畸变项的标定模型对相机内部参数进行校正

首先设定像平面上每个像素的物理尺寸大小为dx*dy(单位：mm)，图像坐标系(x，y)原点在像素坐标系(u，v)中的坐标为(u₀，v₀)，(x，y)是实际图像中像点归一化的坐标，图像中任意像素在两个坐标系中满足如下关系：

相机坐标系中任一点P_c(X_c，Y_c，Z_c)投影到图像坐标系上为(x_c，y_c，f)，图像坐标系平面与光轴z轴垂直，与原点距离为f，根据相似三角形原理可以得出：

引入所述改进的带有非线性畸变项的标定模型，包括由于镜头形状缺陷造成的径向畸变和由于光学***存在不同程度的偏心造成的切向畸变，径向畸变数学模型为：

其中r²＝x²+y²，(x’，y’)为矫正后不含畸变项的理想线性相机坐标系的归一化坐标值，径向畸变值与图像点在图像中的位置有关，图像边缘处的径向畸变值较大，

切向畸变数学模型为：

其中包含k₁、k₂、k₃、p₁、p₂共5个非线性畸变系数，由公式(3)、(4)得畸变矫正函数模型如下：

从世界坐标变换到相机坐标转换存在如下关系：

P_c＝R·(P_W-C)＝R·P_W+T (6)

结合式(1)～(6)，用齐次坐标与矩阵形式可表示为：

M_int、M_ext分别是相机标定内、外参数矩阵，其中相机内部参数包括图像中心点像素值u₀、v₀，f_x、f_y为x轴和y轴上的归一化焦距、通过Java结合OpenCV实现手机相机标定，获取手机相机所述的内部参数、相机镜头畸变参数和图像分辨率v_max、u_max；

步骤二：建立深度提取模型

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

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

α＝F(v，β)＝a·v+b (8)

其中参数a、b均与相机型号和相机旋转角度有关，

当α取最小值α＝α_min＝90-θ-β时，θ为相机垂直视场角的一半，即被拍摄物体投影到图片最底端时，v＝v_max(v_max为相机CMOS或CCD图像传感器列坐标有效像素数)，代入式(8)可得：

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

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

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

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

90-β+θ＝b (11)

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

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

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

步骤三：采集并处理待测立木图像，计算摄影测量坐标系下立木深度值

在立木图像采集步骤中，还包括对待测立木图像的非线性畸变校正和预处理，即：

通过手机相机进行图像采集，建立投影几何模型，其中f为相机焦距，θ为相机垂直视场角的一半，h为相机拍照高度，β为相机沿相机坐标系ox轴的旋转角，相机顺时针旋转β值为正，逆时针为负，β值通过相机内部重力传感器获取，α为目标物成像角度，γ为目标立木所在平面的坡度；

结合步骤一相机标定获取的相机镜头畸变参数，对图像存在的径向畸变和切向畸变误差进行非线性畸变矫正；将矫正后的理想线性归一化坐标值(x，y)代入公式(1)，求算出矫正后图像各点像素坐标值，通过双线性内插的方法对矫正后像素值进行插值处理从而得到矫正后图像；

将矫正后的图像RGB空间转换到Lab颜色空间，对于每个颜色通道L、a、b，计算各个像素点与整幅图像的平均色差并取平方，然后将这三个通道的值相加作为该像素的显著性值得到三个通道的均值图像，同时，将传统的基于频率调谐的视觉显著性中对图像的高斯滤波处理改为双边滤波处理，弥补该算法进行立木显著性表达时树干纹理所带来的影响，取三个通道的均值图像和滤波图像的欧氏距离并求和，得到图像的频率调谐视觉显著图，另外，将HSV颜色空间中S分量融合到频率调谐视觉显著图中，以增强图像显著性，对显著图进一步进行二值化处理后利用形态学膨胀和腐蚀组合运算达到连接邻近物体和平滑边界的作用，最后对立木树干轮廓进行边缘检测，得到目标立木轮廓提取结果，进而计算目标立木与地面接触的边缘的几何中心点纵坐标像素值v；

根据图像中立木树干的几何特征构建立木胸径测量坐标***

立木胸径测量坐标***包括5组坐标系，分别是像平面坐标系、像素坐标系、像空间坐标系、摄影测量坐标系、物方空间坐标系，立木胸径测量坐标***均属于右手坐标系，其中像平面坐标系、像素坐标系、像空间坐标系为图像中所有待测立木共同使用的坐标系，此外，本文还根据待测立木特征为每株立木建立其特属的摄影测量坐标系与物方空间坐标系；

以图像平面的左上角或左下角为原点建立像素坐标系o’-uv；像平面坐标系原点位于图像左上角o’，水平向右为u轴，垂直向下为v轴，各坐标轴均以像素为单位；以图像平面与光轴的交点o为原点建立像平面坐标系o-xy，水平向右为x轴，垂直向下为y轴，其单位为物理尺寸，原点o一般位于图像中心处，o在以像素为单位的图像坐标系中的坐标为(u₀，v₀)；像平面坐标系和像素坐标系属于同一个平面上原点不同的两个坐标系；以相机中心点为坐标原点建立像空间坐标系S-X_cY_cZ_c，图像平面与光轴Z_c轴垂直，和投影中心距离为f；将像空间坐标系点映射到投影平面上的过程为投影变换；另外，本文还为图像中N株待测量立木定义N组摄影测量坐标系D_i-U_iV_iW_i和物方空间坐标系T_i-X_iY_iZ_i(i≤N)，第i株立木对应的摄影测量坐标系原点D_i位于该立木底部几何中心点，V_i轴竖直向上垂直于水平地面，U_i轴沿立木成像面方向水平向右平行于像平面，摄影测量坐标系沿立木树干方向平移1.3m后，旋转一定的角度可得到该株立木的物方空间坐标系，物方空间坐标系Y_i轴沿立木树干方向竖直向上；

根据公式(14)深度提取模型计算出摄影测量坐标系下的立木深度值；

步骤四：测量单幅图像中多株立木胸径值

在针孔相机模型中，根据步骤4中立木胸径测量坐标***建立规则，研究成像点在不同坐标系中的转换关系，建立立木胸径测量模型，

根据立木胸径的定义和立木图像信息，确定第i株立木物方空间坐标系到摄影测量坐标系的旋转矩阵

和平移矩阵

从而计算该立木物方空间坐标系原点在其对应的摄影测量坐标系中的坐标值(U_i0，V_i0，W_i0)，

其中，平移矩阵

针对自然界中立木可能存在一定程度的倾斜角且树干为圆柱体的特征，定义立木胸径测量物方空间坐标系与摄影测量坐标系之间的旋转关系为以V_i轴为主轴的

***，即以V_i轴为主轴旋转0°，然后绕U_i轴旋转

角，最后绕W_i轴旋转Ψ角，因此旋转矩阵

可表示为：

若待测立木中存在倾斜立木，为了方便测量，结合立木树干的几何属性，在进行图片采集时可令倾斜立木树干深度值保持不变，将

角与Ψ角合并为同一个角，即令

摄影测量坐标系经平移和旋转到像空间坐标系是一种刚体运动，摄影测量坐标系分别沿U，V，W轴方向平移t_U，t_V，t_W距离(单位：毫米)后，从摄影测量坐标系旋转到像空间坐标系，其旋转角元素为(κ，β，ω)，旋转过程为以V轴为主轴旋转κ角度后，绕新的U轴旋转β角度，最后绕W轴旋转ω角度，第i株立木的摄影测量坐标系中的任意点(U_i，Y_i，W_i)在像空间坐标系下坐标(X_i，Y_i，Z_i)为：

式中，T表示由摄影测量坐标系到像空间坐标系的平移矩阵，R表示旋转矩阵：

[t_U，t_V，t_W]为：

根据摄影测量坐标***建立规则，摄影测量坐标系绕V轴旋转角度κ等于180°，且通过移动端相机进行图片采集时绕W轴旋转角度ω通常约等于0°，即进行图片采集时，手机竖直置于相机三脚架上，因此，结合公式(15)和公式(17)求算物方空间坐标系原点在像空间坐标系中的坐标值(X_i0，Y_i0，Z_i0)：

像空间坐标系中目标立木深度Z值为：

Z＝-(1300·cosψ-h+D·tanγ)sinβ+D·cosβ (21)

立木胸径测量处在像空间坐标系中Y值为：

Y＝(1300·cosψ-h+D·tanγ)cosβ+D·sinβ (22)

在立木胸径测量立体成像***中，点A’、B’为沿树干方向距地面1.3m处树干在像平面上的投影点，其在像平面上的坐标值为A’(x₁，y₁，f)、B’(x₂，y₂，f)，像点A’、B’在自然场景中的同名物点为A、B，AB所在的直线与图像平面平行，

令点A在相机坐标系下的坐标为(X，Y，Z)，点B的坐标为(X+T_x，Y，Z)，线段CD过AB的中点垂直于立木树干，CD的距离即为立木胸径真实值；在针孔相机成像模型中，基于步骤三中相机镜头组畸变的矫正，图像中像点、相机光心和物点3点共线，则有：

结合公式(1)、(20)～(23)，像点A’、B’的纵坐标像素值v_DBH可表示为：

Y值相同且深度值Z相等的两点A’、B’的水平视差d可表示为：

d＝x₂-x₁＝f·T_x/Z

＝(u₂-u₁)·d_x (25)

其中，u₁、u₂分别是像点A’、B’的横坐标像素值，在已知相机内部参数的情况下，

结合像空间坐标系中待测立木深度值Z，可计算点AB的距离T_x：

因此，该株立木的胸径可表示为：

DBH＝T_x·cosψ (27)。

与现有技术相比本发明的有益效果是：由于采用上述技术方案，

1、与传统立木胸径测量方法相比，该方法节约了劳动力提高了立木因子测量效率；

2、与其它非接触式立木胸径测量方法相比，该方法设备操作简单且便于携带，同时基于智能手机的立木因子测量设备便宜，普适性较高；

3、该方法不需要在测量场景中设定已知尺寸的标定物，测量方法更灵活，且不需要大场景标定场地，避免了数据拟合引起的误差；

4、深度提取模型具有设备通用性，并将相机旋转角引入模型，对于不同型号的相机，仅需第一次通过相机标定获取相机内部参数后，即可计算单幅图片上任意像点的深度；

5、现有的基于图像信息的立木胸径自动化测量方法多是单株立木胸径测量，该方法基于单目视觉单幅图像实现了多株立木胸径的自动化测量，提高了测量效率。

6、经验证，使用该方法测量深度信息在距离为2～10m时平均相对误差为2.32％，胸径小于20cm的树木，绝对误差小于0.237cm，胸径大于或等于20cm的树木，相对误差小于1.272％，满足国家森林资源连续清查中对胸径测量精度的要求。

附图说明

图1是本发明的测距方法流程图；

图2是新型标靶图；

图3是角点检测算法实现流程图；

图4是相机上视角高于水平线拍摄几何模型图；

图5是相机上视角低于水平线拍摄几何模型图；

图6是相机拍摄投影几何模型；

图7是待测立木原始图像；

图8是待测立木频率调谐视觉显著性图；

图9是待测立木树干主轮廓图；

图10是立木胸径测量坐标***；

图11是立木投影几何模型；

图12是三种型号设备物体纵坐标像素值与成像角度之间的关系示意图；

图13是不同相机旋转角度物体纵坐标像素值与实际成像角度之间的关系。

具体实施方式

为了使本发明的技术方案更加清晰，以下结合附图1至13，对本发明进行详细说明。应当理解的是，本说明书中描述的具体实施方式仅仅是为了解释本发明，并不是为了限定本发明的保护范围。

本发明是一种基于深度提取模型的多株立木胸径被动测量方法，包括如下步骤：

一、对手机相机进行标定，获取相机内部参数和图像分辨率。所述标定采用张正友标定法，并引入改进的带有非线性畸变项的标定模型对相机内部参数进行校正。

首先设定像平面上每个像素的物理尺寸大小为dx*dy(单位：mm)，图像坐标系(x，y)原点在像素坐标系(u，v)中的坐标为(u₀，v₀)，(x，y)是实际图像中像点归一化的坐标，图像中任意像素在两个坐标系中满足如下关系：

相机坐标系中任一点P_c(X_c，Y_c，Z_c)投影到图像坐标系上为(x_c，y_c，f)，图像坐标系平面与光轴z轴垂直，与原点距离为f，根据相似三角形原理可以得出：

引入所述改进的带有非线性畸变项的标定模型，包括由于镜头形状缺陷造成的径向畸变和由于光学***存在不同程度的偏心造成的切向畸变，径向畸变数学模型为：

其中r²＝x²+y²，(x’，y’)为矫正后不含畸变项的理想线性相机坐标系的归一化坐标值，径向畸变值与图像点在图像中的位置有关，图像边缘处的径向畸变值较大，

切向畸变模型数学模型为：

其中包含k₁、k₂、k₃、p₁、p₂共5个非线性畸变系数，由公式(3)、(4)得畸变矫正函数模型如下：

从世界坐标变换到相机坐标转换存在如下关系：

P_c＝R·(P_W-C)＝R·P_W+T (6)

结合式(1)～(6)，用齐次坐标与矩阵形式可表示为：

M_int、M_ext分别是相机标定内、外参数矩阵，其中相机内部参数包括图像中心点像素值u₀、v₀，f_x、f_y为x轴和y轴上的归一化焦距、通过Java结合OpenCV实现手机相机标定，获取手机相机所述的内部参数、相机镜头畸变参数和图像分辨率v_max、u_max；

二、通过对新型标靶图像的采集，建立相机深度提取模型。现有的标靶为长宽相等的黑白棋盘格标靶。本发明的新型靶标与现有靶标的区别在于，设定标靶距离相机最近的一排方格大小为d*dmm，随后每排方格的宽度为固定，长度后一排比前一排增加值为

下式中x_i为第i个角点到相机的实际距离，y_i为各方格的长度，则相邻方格长度的差值Δd_i为：

设各方格的计算长度与实际距离之间的关系为f(x)，根据公式(28)可得：

经Pearson相关性分析，长度与实际距离之间呈极显著线性相关关系(p＜0.01)，相关系数r等于0.975，通过最小二乘法可以求算出f(x)的导数f’(x)，

因此，当该标靶距离相机最近的一排方格大小为d*d mm(d的范围取30～60mm时测量精度最高)时，随后每排宽度固定、长度增加值

为d*f’(x)mm，新型标靶如图2所示，

拍摄水平地面上物体时存在透视变换现象使得Harris和Shi-Tomasi等常见的角点检测算法鲁棒性较差，且当相机沿相机坐标系ox轴逆时针旋转角度较大时也会检测失败，因此本发明结合Andreas Geiger等提出的基于生长的棋盘格角点检测方法和OpenCV提供的cornerSubPix()函数进行亚像素级角点位置检测，该算法鲁棒性高，对畸变程度较大的图片提取效果较好，

角点检测算法的实现流程如图3所示，本发明上述新型靶标的亚像素级角点提取步骤为：

1)根据图像中各像素点与模板的相似度参数在图像上寻找角点，定位标靶角点位置；

首先定义两种不同的角点模板，一种用于和坐标轴平行的角点，另一种用于旋转45°的角点，每个模板由4个滤波核{A，B，C，E}组成，用于后面和图像进行卷积操作；然后利用这两个角点模板来计算每个拐点与角点的相似度：

其中

表示卷积核X(X＝A，B，C，E)和模板i(i＝1，2)在某个像素点的卷积响应，

和

表示模板i的两种可能的拐点的相似度，计算图像中每个像素点的相似度可以得到角点相似图；利用非极大值抑制算法对角点像素图进行处理来获得候选点；然后用梯度统计的方法在一个局域的n x n邻域内验证这些候选点，先对局域灰度图进行sobel滤波，然后计算加权方向直方图(32bins)，用meanshift算法找到其中的两个主要的模态γ₁和γ₂；根据边缘的方向，对于期望的梯度强度

构造一个模板T。

(*表示互相关操作符)和角点相似度的乘积作为角点分值，然后用阈值进行判断就得到初始角点。

2)对角点的位置和方向进行亚像素级精细提取；

运用OpenCV中的cornerSubPix()函数进行亚像素级角点定位，将角点定位到子像素，从而取得亚像素级别的角点检测效果；为细化边缘方向向量，根据图像梯度值最小化其标准离差率：

其中

是相邻像素集，与模块i的梯度值m_i＝[cos(γ_i)sin(γ_i)]^T相匹配。(求算方案依据文献Geiger A，Moosmann F，Car

etal.Automatic camera and range sensor calibration using a single shot[C]//Robotics and Automation(ICRA)，2012IEEE International Conference on.IEEE，2012：3936-3943.)

3)最后是标记角点并输出其亚像素级坐标，根据能量函数生长并重建棋盘格，标记角点，输出亚像素级角点坐标；

根据文献“Geiger A，Moosmann F，Car

et al.Automatic camera and rangesensor calibration using a single shot[C]//Robotics and Automation(ICRA)，2012IEEE International Conference on.IEEE，2012：3936-3943.”提供的方法优化能量函数，重建棋盘格并标记角点，能量生长函数公式为：

E(x，y)＝E_corners(y)+E_struct(x，y) (32)

其中，E_corners是当前棋盘角点总数的负值，E_struct是两个相邻角点和预测角点的匹配程度；通过OpenCV输出角点像素值。

使用SPSS 22对物体成像角度、纵坐标像素值进行线性相关分析，输出Pearson相关系数，经验证，不同型号的设备和相机旋转角度下，物体纵坐标像素值与实际成像角度呈极显著负相关关系(p＜0.01)，另外，本发明还对不同的设备型号和相机旋转角度下物体纵坐标像素值与成像角度之间线性函数的斜率差异进行显著性检验，结果表明，不同设备型号和相机旋转角度下物体纵坐标像素值与成像角度之间线性函数的斜率差异极显著(p＜0.01)，说明不同型号的设备和相机旋转角度，其深度提取模型有所不同，

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

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

α＝F(v，β)＝a·v+b (8)

其中参数a、b均与相机型号和相机旋转角度有关，

当α取最小值α＝α_min＝90-θ-β时，θ为相机垂直视场角的一半，即被拍摄物体投影到图片最底端时，v＝v_max(v_max为相机CMOS或CCD图像传感器列坐标有效像素数)，代入式(8)可得：

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

当α_min+-2θ＞90°，即θ＞β时，此时相机上视角高于水平线，相机拍摄投影几何模型如图4，地平面无限远处，α无限接近于90°，此时v无限趋近于v₀-tanβ*f_y，f_y为像素单位下相机的焦距，β为负值即相机逆时针旋转时亦同理，因此，代入式(8)可得：

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

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

90-β+θ＝b (11)

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

公式(12)中L_CMOs为相机CMOS或CCD图像传感器的边长，结合式(9)～(12)，手机相机深度提取模型F(α，β)为：

公式(13)中δ为相机非线性畸变项误差，结合手机相机拍摄高度h，根据三角函数原理计算目标点深度值D为：

三、通过对待测目标物的图像采集，获取目标点像素值u、v。通过手机相机进行图像采集，建立投影几何模型如图6，其中f为相机焦距，θ为相机垂直视场角的一半，h为相机拍照高度，β为相机沿相机坐标系ox轴的旋转角，相机顺时针旋转β值为正，逆时针为负，β值通过相机内部重力传感器获取，α为目标物成像角度；结合第一步相机标定获取的相机镜头畸变参数，对图像存在的径向畸变和切向畸变误差进行非线性畸变矫正；将矫正后的理想线性归一化坐标值(x，y)代入公式(1)，求算出矫正后图像各点像素坐标值，通过双线性内插的方法对矫正后像素值进行插值处理从而得到矫正后图像；

图7为自然场景中立木树干相机拍摄原始图片灰度化结果，自然环境下多株立木的视觉显著性表达及轮廓检测，受复杂自然背景、树干纹理等因素的影响，直接使用感兴趣区域提取的方法，效果不佳。针对上述问题，本文以自然环境中的待测立木为研究对象，采用一种改进的频率调谐的视觉显著性轮廓检测算法，该算法包括颜色空间选取、频率调谐的视觉显著性表达以及立木树干轮廓提取三个部分。

根据Lab颜色模型中L、a、b三个分量的几何距离存在差异，且3个分量之间独立性高的特点，本文将图像RGB空间转换到Lab颜色空间，在此空间引入对自然环境中立木轮廓目标区域的频率调谐视觉显著性描述，使之成为识别目标的特征之一；采用Lab颜色空间作为图像特征，对于每个颜色通道L、a、b，计算各个像素点与整幅图像的平均色差并取平方，然后将这三个通道的值相加作为该像素的显著性值得到三个通道的均值图像；同时，由于依靠传统频率调谐显著性算法受树干纹理影响较大，本文将频率调谐的视觉显著性中对图像的高斯滤波处理改为双边滤波处理，弥补该算法进行立木显著性表达时树干纹理所带来的影响。取三个通道的均值图像和滤波图像的欧氏距离并求和，得到图像的频率调谐视觉显著图；另外，本文还将HSV颜色空间中S分量融合到频率调谐视觉显著图中，可以增强图像显著性，如图8所示。

对显著图进行二值化处理后利用形态学膨胀和腐蚀组合运算达到连接邻近物体和平滑边界的作用，最后对立木树干轮廓进行边缘检测，得到目标立木轮廓提取结果，图9是为了突出树干对轮廓进行填充的结果。

本发明定义了5组坐标系，分别是像平面坐标系、像素坐标系、像空间坐标系、摄影测量坐标系、物方空间坐标系。立木胸径测量坐标***均属于右手坐标系。在进行实际测量中，需要确定各坐标系之间的相对关系，各坐标系间的几何关系示意图如图10所示。

以图像平面的左上角或左下角为原点建立像素坐标系o’-uv。假设像平面坐标系原点位于图像左上角o’，水平向右为u轴，垂直向下为v轴，各坐标轴均以像素为单位；以图像平面与光轴的交点o为原点建立像平面坐标系o-xy，水平向右为x轴，垂直向下为y轴，其单位为物理尺寸，原点o一般位于图像中心处，o在以像素为单位的图像坐标系中的坐标为(u₀，v₀)；像平面坐标系和像素坐标系属于同一个平面上原点不同的两个坐标系。以相机中心点为坐标原点建立像空间坐标系S-X_cY_cZ_c，图像平面与光轴Z_c轴垂直，和投影中心距离为f。将像空间坐标系点映射到投影平面上的过程为投影变换；另外，本文还为图像中N株待测量立木定义N组摄影测量坐标系D_i-U_iV_iW_i和物方空间坐标系T_i-X_iY_iZ_i(i≤N)，第i株立木对应的摄影测量坐标系原点D_i位于该立木底部几何中心点，V_i轴竖直向上垂直于水平地面，U_i轴沿立木成像面方向水平向右平行于像平面，摄影测量坐标系沿立木树干方向平移1.3m后，旋转一定的角度可得到该株立木的物方空间坐标系，物方空间坐标系Y_i轴沿立木树干方向竖直向上。

四、测量单幅图像中多株立木胸径值。

在针孔相机模型中，根据步骤三中立木胸径测量坐标***建立规则，研究成像点在不同坐标系中的转换关系，建立立木胸径测量模型，

根据立木胸径的定义和立木图像信息，确定第i株立木物方空间坐标系到摄影测量坐标系的旋转矩阵

和平移矩阵

从而计算该立木物方空间坐标系原点在其对应的摄影测量坐标系中的坐标值(U_i0，V_i0，W_i0)，

其中，平移矩阵

针对自然界中立木可能存在一定程度的倾斜角且树干为圆柱体的特征，定义立木胸径测量物方空间坐标系与摄影测量坐标系之间的旋转关系为以V_i轴为主轴的

***，即以V_i轴为主轴旋转0°，然后绕U_i轴旋转

角，最后绕W_i轴旋转Ψ角，因此旋转矩阵

可表示为：

若待测立木中存在倾斜立木，为了方便测量，结合立木树干的几何属性，在进行图片采集时可令倾斜立木树干深度值保持不变，将

角与Ψ角合并为同一个角，即令

结合式(15)和(16)可得：

基于步骤三中提取的立木树干轮廓，可通过提取图像中立木轮廓的最小外接矩形获得旋转角Ψ。根据自然环境下的立木特点，本文采用优化的最小外接矩形直接计算方法提取最小外接矩形，该算法首先按照最小外接矩形直接计算方法获取立木树干轮廓区域外接矩形，此时该矩形倾斜角为0°；然后以外接矩形中心为中心点旋转该外接矩形，直至外接矩形面积取最小值时记录旋转角，从而找到最优矩形姿态，该旋转角即为旋转角Ψ。

摄影测量坐标系经平移和旋转到像空间坐标系是一种刚体运动，摄影测量坐标系分别沿U，V，W轴方向平移t_U，t_V，t_W距离(单位：毫米)后，从摄影测量坐标系旋转到像空间坐标系，其旋转角元素为(κ，β，ω)，旋转过程为以V轴为主轴旋转κ角度后，绕新的U轴旋转β角度，最后绕W轴旋转ω角度，第i株立木的摄影测量坐标系中的任意点(U_i，V_i，W_i)在像空间坐标系下坐标(X_i，Y_i，Z_i)为：

式中，T表示由摄影测量坐标系到像空间坐标系的平移矩阵，R表示旋转矩阵：

[t_U，t_V，t_W]为：

根据摄影测量坐标***建立规则，摄影测量坐标系绕V轴旋转角度κ等于180°，且通过移动端相机进行图片采集时绕W轴旋转角度ω通常约等于0°，即进行图片采集时，手机竖直置于相机三脚架上，因此，结合公式(15)和公式(17)求算物方空间坐标系原点在像空间坐标系中的坐标值(X_i0，Y_i0，Z_i0)：

像空间坐标系中目标立木深度Z值为：

Z＝-(1300·cosψ-h+D·tanγ)sinβ+D·cosβ (21)

立木胸径测量处在像空间坐标系中Y值为：

Y＝(1300·cosψ-h+D·tanγ)cosβ+D·sinβ (22)

如图11所示，在立木胸径测量立体成像***中，点A’、B’为沿树干方向距地面1.3m处树干在像平面上的投影点，其在像平面上的坐标值为A’(x₁，y₁，f)、B’(x₂，y₂，f)，像点A’、B’在自然场景中的同名物点为A、B，AB所在的直线与图像平面平行，

令点A在相机坐标系下的坐标为(X，Y，Z)，点B的坐标为(X+T_x，Y，Z)，线段CD过AB的中点垂直于立木树干，CD的距离即为立木胸径真实值；在针孔相机成像模型中，基于步骤三中相机镜头组畸变的矫正，图像中像点、相机光心和物点3点共线，则有：

结合公式(1)、(20)～(23)，像点A’、B’的纵坐标像素值v_DBH可表示为：

Y值相同且深度值Z相等的两点A’、B’的水平视差d可表示为：

d＝x₂-x₁＝f·T_x/Z

＝(u₂-u₁)·d_x (25)

其中，u₁、u₂分别是像点A’、B’的横坐标像素值，在已知相机内部参数的情况下，

结合像空间坐标系中待测立木深度值Z，可计算点AB的距离T_x：

因此，该株立木的胸径可表示为：

DBH＝T_x·cosψ (27)。

实施例1

下面以小米3(MI 3)手机为例，具体说明本发明的基于单目视觉的优化的深度提取和被动测距方法。

一、对手机相机进行标定，获取相机内部参数和图像分辨率

使用行列数为8*9尺寸大小为20*20的棋盘格标定板作为相机标定的实验材料，通过小米3手机相机采集20张不同角度的标定板图片，利用OpenCV根据上述改进的带有非线性畸变项的相机标定模型对小米3(MI 3)手机相机进行标定，

首先使用fin()函数读取标定板图片，并通过.cols、.rows获取第一张图片的图像分辨率；然后通过find4QuadComerSubpix()函数提取标定板图片中亚像素级角点，并用drawChessboardCorners()函数标记角点；调用calibrateCamera()函数对相机进行标定，将得到的摄像机内外参数用于对空间的三维点进行重新投影计算，得到新的投影点，计算新的投影点和旧的投影点之间的误差；输出并保存相机内参矩阵和畸变参数，

标定所得相机内部参数为：f_x＝3486.5637，u₀＝1569.0383，f_y＝3497.4652，v₀＝2107.9899，图像分辨率为3120×4208，相机镜头畸变参数为：[0.0981，-0.1678，0.0003，-0.0025，0.0975]，

二、通过对新型标靶图像的采集，建立相机深度提取模型

本发明使用45*45mm的传统棋盘格标定板作为标靶设计的初始实验材料，为计算相邻方格长度的差值，本发明设计了6组实验，提取方格大小为45*45mm的传统棋盘格角点值，并求算出相邻角点间单位像素在世界坐标系下代表的实际物理距离，为保证角点间纵坐标像素差值大致相等，各方格的长度yi的值如表1所示，

表1各方格的计算宽度

Table 1 Computing width of each grid

经Pearson相关性分析，长度与实际距离之间呈极显著线性相关关系(p＜0.01)，相关系数r等于0.975，通过最小二乘法可以求算出f(x)的导数f’(x)＝0.262，因此，当该标靶距离相机最近的一排方格大小为45*45mm时，随后每排宽度固定、宽度增加值Δd为11.79mm，

通过具体实施步骤中的角点提取算法提取该新型标靶的角点，

本发明分别选取小米、华为、iPhone三种不同型号的智能手机作为图像采集设备，相机旋转角度β＝{-10°，0°，10°，20°，30°}。使用所述角点检测算法采集数据，并对其关系进行函数拟合，图12为β＝10°时三种不同型号的智能手机纵坐标像素值与物体成像角度之间的关系，图13为不同相机旋转角度下纵坐标像素值与物体成像角度之间的关系，

不同型号的相机设备和相机旋转角度，随着纵坐标像素值的增加，物体成像角度呈递减趋势，且设备型号和相机旋转角度的不同使得像素值与成像角度之间呈现不同的线性函数关系，使用SPSS 22对物体成像角度、纵坐标像素值进行线性相关分析，输出Pearson相关系数r如表2所示。

表2物体纵坐标像素值与成像角度相关系数

Table 2 Pearson correlation coefficient of image ordinate pixelvalues and actual imaging angles

注：**表示极显著(p＜0.01)。

Note：**represents very significant correlation(p＜0.01).

经验证，同型号的设备和相机旋转角度下，物体纵坐标像素值与实际成像角度呈极显著负相关关系(p＜0.01)，相关系数r大于0.99。另外，本发明还对不同的设备型号和相机旋转角度下物体纵坐标像素值与成像角度之间线性函数的斜率差异进行显著性检验。结果表明，不同设备型号和相机旋转角度下物体纵坐标像素值与成像角度之间线性函数的斜率差异极显著(p＜0.01)，说明不同型号的设备和相机旋转角度，其深度提取模型有所不同。

根据具体实施例所述深度提取模型，将小米3手机相机内部参数代入模型，得到该设备的具体深度提取模型为：

将相机旋转角度β代入公式(34)，根据三角函数原理计算目标点深度值D(单位：mm)：

三、通过对待测目标物的图像采集，获取目标点像素值u、v。

使用小米手机3(MI 3)相机作为图片采集设备，通过相机三脚架进行图片采集，并测量相机到地面的高度h等于1300mm，相机旋转角β等于0°，

对图像存在的径向畸变和切向畸变误差进行非线性畸变校正；

根据第一步相机标定获取的相机镜头畸变参数：[0.0981，-0.1678，0.0003，-0.0025，0.0975]，根据公式(5)计算矫正后理想线性归一化坐标值：

结合公式(1)和(2)计算矫正后图像各点像素坐标值，通过双线性插值处理得到矫正后图像；

本发明以自然场景中的立木为例提取其树杆轮廓，测量胸径。首先对矫正后的立木图象进行视觉显著性表达，增强图像显著性，对显著图进行二值化处理后利用形态学膨胀和腐蚀组合运算达到连接邻近物体和平滑边界的作用，最后对立木树干轮廓进行边缘检测，得到目标立木轮廓提取结果，进而获取立木最底端几何中心点像素值；

定义立木胸径测量坐标***，将相机拍照高度h、旋转角度β和立木最底端几何中心点纵坐标像素值v代入公式(35)可以计算出各立木在摄影测量坐标系下的深度值D(单位：mm)实验结果如表3所示；

四、测量图像中多株立木胸径值

根据具体实施例中立木胸径测量方法，计算图像中每株立木胸径值，实验结果如表3所示：

表3立木胸径测量数据

Table 3 Data of trees DBH measurement

Claims (1)

1.一种基于深度提取模型的多株立木胸径被动测量方法，其特征在于包括如下步骤：

步骤一：对手机相机进行标定，获取相机内部参数和图像分辨率

采用张正友标定法，并引入改进的带有非线性畸变项的标定模型对相机内部参数进行校正

首先设定像平面上每个像素的物理尺寸大小为dx*dy，单位：mm，图像坐标系(x，y)原点在像素坐标系(u，v)中的坐标为(u₀，v₀)，(x,y)是实际图像中像点归一化的坐标，图像中任意像素在两个坐标系中满足如下关系：

相机坐标系中任一点P_c(X_c,Y_c,Z_c)投影到图像坐标系上为(x_c,y_c,f)，图像坐标系平面与光轴z轴垂直，与原点距离为f，根据相似三角形原理可以得出：

引入所述改进的带有非线性畸变项的标定模型，包括由于镜头形状缺陷造成的径向畸变和由于光学***存在不同程度的偏心造成的切向畸变，径向畸变数学模型为：

其中r²＝x²+y²，(x’,y’)为矫正后不含畸变项的理想线性相机坐标系的归一化坐标值，径向畸变值与图像点在图像中的位置有关，图像边缘处的径向畸变值较大，

切向畸变数学模型为：

其中包含k₁、k₂、k₃、p₁、p₂共5个非线性畸变系数，由公式(3)、(4)得畸变矫正函数模型如下：

从世界坐标变换到相机坐标转换存在如下关系：

P_c＝R·(P_W-C)＝R·P_W+T (6)

公式(6)中T表示由摄影测量坐标系到像空间坐标系的平移矩阵，R表示旋转矩阵，结合式(1)～(6)，用齐次坐标与矩阵形式可表示为：

M_int、M_ext分别是相机标定内、外参数矩阵，其中相机内部参数包括图像中心点像素值u₀、v₀，f_x、f_y为x轴和y轴上的归一化焦距、通过Java结合OpenCV实现手机相机标定，获取手机相机所述的内部参数、相机镜头畸变参数和图像分辨率v_max、u_max；

步骤二：建立深度提取模型

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

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

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

其中参数a、b均与相机型号和相机旋转角度有关，

当α取最小值α＝α_min＝90-θ-β时，θ为相机垂直视场角的一半，即被拍摄物体投影到图片最底端时，v＝v_max，v_max为相机CMOS或CCD图像传感器列坐标有效像素数，代入式(8)可得：

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

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

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

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

90-β+θ＝b (11)

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

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

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

步骤三：采集并处理待测立木图像，计算摄影测量坐标系下立木深度值

在立木图像采集步骤中，还包括对待测立木图像的非线性畸变校正和预处理，即：

通过手机相机进行图像采集，建立投影几何模型，其中f为相机焦距，θ为相机垂直视场角的一半，h为相机拍照高度，β为相机沿相机坐标系ox轴的旋转角，相机顺时针旋转β值为正，逆时针为负，β值通过相机内部重力传感器获取，α为目标物成像角度，γ为目标立木所在平面的坡度；

结合步骤一相机标定获取的相机镜头畸变参数，对图像存在的径向畸变和切向畸变误差进行非线性畸变矫正；将矫正后的理想线性归一化坐标值(x,y)代入公式(1)，求算出矫正后图像各点像素坐标值，通过双线性内插的方法对矫正后像素值进行插值处理从而得到矫正后图像；

将矫正后的图像RGB空间转换到Lab颜色空间，对于每个颜色通道L、a、b，计算各个像素点与整幅图像的平均色差并取平方，然后将这三个通道的值相加作为该像素的显著性值得到三个通道的均值图像，同时，将传统的基于频率调谐的视觉显著性中对图像的高斯滤波处理改为双边滤波处理，弥补该算法进行立木显著性表达时树干纹理所带来的影响，取三个通道的均值图像和滤波图像的欧氏距离并求和，得到图像的频率调谐视觉显著图，另外，将HSV颜色空间中S分量融合到频率调谐视觉显著图中，以增强图像显著性，对显著图进一步进行二值化处理后利用形态学膨胀和腐蚀组合运算达到连接邻近物体和平滑边界的作用，最后对立木树干轮廓进行边缘检测，得到目标立木轮廓提取结果，进而计算目标立木与地面接触的边缘的几何中心点纵坐标像素值v；

根据图像中立木树干的几何特征构建立木胸径测量坐标***

立木胸径测量坐标***包括5组坐标系，分别是像平面坐标系、像素坐标系、像空间坐标系、摄影测量坐标系、物方空间坐标系，立木胸径测量坐标***均属于右手坐标系，其中像平面坐标系、像素坐标系、像空间坐标系为图像中所有待测立木共同使用的坐标系，此外，根据待测立木特征为每株立木建立其特属的摄影测量坐标系与物方空间坐标系；

以图像平面的左上角或左下角为原点建立像素坐标系o’-uv；像平面坐标系原点位于图像左上角o’，水平向右为u轴，垂直向下为v轴，各坐标轴均以像素为单位；以图像平面与光轴的交点o为原点建立像平面坐标系o-xy，水平向右为x轴，垂直向下为y轴，其单位为物理尺寸，原点o位于图像中心处，o在以像素为单位的图像坐标系中的坐标为(u₀,v₀)；像平面坐标系和像素坐标系属于同一个平面上原点不同的两个坐标系；以相机中心点为坐标原点建立像空间坐标系S-X_cY_cZ_c，图像平面与光轴Z_c轴垂直，和投影中心的距离为f；将像空间坐标系点映射到投影平面上的过程为投影变换；另外，为图像中N株待测量立木定义N组摄影测量坐标系D_i-U_iV_iW_i和物方空间坐标系T_i-X_iY_iZ_i(i≤N)，第i株立木对应的摄影测量坐标系原点D_i位于该立木底部几何中心点，V_i轴竖直向上垂直于水平地面，U_i轴沿立木成像面方向水平向右平行于像平面，摄影测量坐标系沿立木树干方向平移1.3m后，旋转一定的角度可得到该株立木的物方空间坐标系，物方空间坐标系Y_i轴沿立木树干方向竖直向上；

根据公式(14)深度提取模型计算出摄影测量坐标系下的立木深度值；

步骤四：测量单幅图像中多株立木胸径值

在针孔相机模型中，根据步骤4中立木胸径测量坐标***建立规则，研究成像点在不同坐标系中的转换关系，建立立木胸径测量模型，

根据立木胸径的定义和立木图像信息，确定第i株立木物方空间坐标系到摄影测量坐标系的旋转矩阵

和平移矩阵

从而计算该立木物方空间坐标系原点在其对应的摄影测量坐标系中的坐标值(U_i0,V_i0,W_i0)，

其中，平移矩阵

定义立木胸径测量物方空间坐标系与摄影测量坐标系之间的旋转关系为以V_i轴为主轴的

***，即以V_i轴为主轴旋转0°，然后绕U_i轴旋转

角，最后绕W_i轴旋转Ψ角，因此旋转矩阵

可表示为：

若待测立木中存在倾斜立木，为了方便测量，结合立木树干的几何属性，在进行图片采集时可令倾斜立木树干深度值保持不变，将

角与Ψ角合并为同一个角，令

摄影测量坐标系经平移和旋转到像空间坐标系是一种刚体运动，摄影测量坐标系分别沿U,V,W轴方向平移t_U，t_V，t_W距离后，单位：毫米，从摄影测量坐标系旋转到像空间坐标系，其旋转角元素为(κ，β，ω)，旋转过程为以V轴为主轴旋转κ角度后，绕新的U轴旋转β角度，最后绕W轴旋转ω角度，第i株立木的摄影测量坐标系中的任意点(U_i,V_i,W_i)在像空间坐标系下坐标(X_i,Y_i,Z_i)为：

式中，T表示由摄影测量坐标系到像空间坐标系的平移矩阵，R表示旋转矩阵：

[t_U,t_V,t_W]为：

根据摄影测量坐标***建立规则，摄影测量坐标系绕V轴旋转角度κ等于180°，且通过移动端相机进行图片采集时绕W轴旋转角度ω等于0°，即进行图片采集时，手机竖直置于相机三脚架上，因此，结合公式(15)和公式(17)求算物方空间坐标系原点在像空间坐标系中的坐标值(X_i0,Y_i0,Z_i0)：

像空间坐标系中目标立木深度Z值为：

Z＝-(1300·cosψ-h+D·tanγ)sinβ+D·cosβ (21)

立木胸径测量处在像空间坐标系中Y值为：

Y＝(1300·cosψ-h+D·tanγ)cosβ+D·sinβ (22)

在立木胸径测量立体成像***中，点A’、B’为沿树干方向距地面1.3m处树干在像平面上的投影点，其在像平面上的坐标值为A’(x₁,y₁,f)、B’(x₂,y₂,f)，像点A’、B’在自然场景中的同名物点为A、B，AB所在的直线与图像平面平行，

令点A在相机坐标系下的坐标为(X,Y,Z)，点B的坐标为(X+T_x,Y,Z)，线段CD过AB的中点垂直于立木树干，CD的距离即为立木胸径真实值；在针孔相机成像模型中，基于步骤三中相机镜头组畸变的矫正，图像中像点、相机光心和物点3点共线，则有：

结合公式(1)、(20)～(23)，像点A’、B’的纵坐标像素值v_DBH可表示为：

Y值相同且深度值Z相等的两点A’、B’的水平视差d可表示为：

其中，u₁、u₂分别是像点A’、B’的横坐标像素值，在已知相机内部参数的情况下，结合像空间坐标系中待测立木深度值Z，可计算点A、B之间的距离T_x：

因此，该株立木的胸径可表示为：

DBH＝T_x·cosψ (27)。

Images