CN109269430B - 基于深度提取模型的多株立木胸径被动测量方法 - Google Patents
基于深度提取模型的多株立木胸径被动测量方法 Download PDFInfo
- Publication number
- CN109269430B CN109269430B CN201810913273.3A CN201810913273A CN109269430B CN 109269430 B CN109269430 B CN 109269430B CN 201810913273 A CN201810913273 A CN 201810913273A CN 109269430 B CN109269430 B CN 109269430B
- Authority
- CN
- China
- Prior art keywords
- image
- coordinate system
- camera
- value
- standing
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 210000000481 breast Anatomy 0.000 title claims abstract description 58
- 238000000605 extraction Methods 0.000 title claims abstract description 36
- 238000000691 measurement method Methods 0.000 title abstract description 9
- 238000005259 measurement Methods 0.000 claims abstract description 55
- 238000000034 method Methods 0.000 claims abstract description 49
- 238000012545 processing Methods 0.000 claims abstract description 19
- 210000000038 chest Anatomy 0.000 claims abstract description 7
- 238000003384 imaging method Methods 0.000 claims description 44
- 239000011159 matrix material Substances 0.000 claims description 30
- 238000012937 correction Methods 0.000 claims description 18
- 230000000007 visual effect Effects 0.000 claims description 16
- 238000004422 calculation algorithm Methods 0.000 claims description 15
- 230000003287 optical effect Effects 0.000 claims description 15
- 238000001914 filtration Methods 0.000 claims description 10
- 238000013519 translation Methods 0.000 claims description 10
- 230000000694 effects Effects 0.000 claims description 7
- 230000009466 transformation Effects 0.000 claims description 7
- 238000013178 mathematical model Methods 0.000 claims description 6
- 230000007547 defect Effects 0.000 claims description 5
- 238000006243 chemical reaction Methods 0.000 claims description 4
- 230000007797 corrosion Effects 0.000 claims description 4
- 238000005260 corrosion Methods 0.000 claims description 4
- 238000003708 edge detection Methods 0.000 claims description 4
- 238000013507 mapping Methods 0.000 claims description 4
- 230000000877 morphologic effect Effects 0.000 claims description 4
- 230000002146 bilateral effect Effects 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 3
- 238000006467 substitution reaction Methods 0.000 claims description 3
- 230000010339 dilation Effects 0.000 claims description 2
- 238000011835 investigation Methods 0.000 abstract description 3
- 238000001514 detection method Methods 0.000 description 9
- 238000005516 engineering process Methods 0.000 description 7
- 244000062793 Sorghum vulgare Species 0.000 description 6
- 235000019713 millet Nutrition 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000012886 linear function Methods 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 4
- 238000013459 approach Methods 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 230000010354 integration Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000012360 testing method Methods 0.000 description 3
- 238000012795 verification Methods 0.000 description 3
- 239000002023 wood Substances 0.000 description 3
- 238000010220 Pearson correlation analysis Methods 0.000 description 2
- 238000010219 correlation analysis Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 239000000284 extract Substances 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 240000000249 Morus alba Species 0.000 description 1
- 235000008708 Morus alba Nutrition 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000007935 neutral effect Effects 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
- 239000013598 vector Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/08—Measuring arrangements characterised by the use of optical techniques for measuring diameters
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/80—Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Length Measuring Devices By Optical Means (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
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)中的坐标为(u0,v0),(x,y)是实际图像中像点归一化的坐标,图像中任意像素在两个坐标系中满足如下关系:
相机坐标系中任一点Pc(Xc,Yc,Zc)投影到图像坐标系上为(xc,yc,f),图像坐标系平面与光轴z轴垂直,与原点距离为f,根据相似三角形原理可以得出:
引入所述改进的带有非线性畸变项的标定模型,包括由于镜头形状缺陷造成的径向畸变和由于光学***存在不同程度的偏心造成的切向畸变,径向畸变数学模型为:
其中r2=x2+y2,(x’,y’)为矫正后不含畸变项的理想线性相机坐标系的归一化坐标值,径向畸变值与图像点在图像中的位置有关,图像边缘处的径向畸变值较大,
切向畸变数学模型为:
其中包含k1、k2、k3、p1、p2共5个非线性畸变系数,由公式(3)、(4)得畸变矫正函数模型如下:
从世界坐标变换到相机坐标转换存在如下关系:
Pc=R·(PW-C)=R·PW+T (6)
结合式(1)~(6),用齐次坐标与矩阵形式可表示为:
Mint、Mext分别是相机标定内、外参数矩阵,其中相机内部参数包括图像中心点像素值u0、v0,fx、fy为x轴和y轴上的归一化焦距、通过Java结合OpenCV实现手机相机标定,获取手机相机所述的内部参数、相机镜头畸变参数和图像分辨率vmax、umax;
步骤二:建立深度提取模型
根据目标物成像角度α与纵坐标像素值v之间的线性关系设定抽象函数,建立含目标物成像角度α、纵坐标像素值v和相机旋转角β三个参数空间关系模型,即α=F(v,β),
不同型号的设备和相机旋转角度下,被拍摄物体纵坐标像素值与成像角度均呈极显著负线性相关关系,且该线性关系的斜率与截距有所不同,故设:
α=F(v,β)=a·v+b (8)
其中参数a、b均与相机型号和相机旋转角度有关,
当α取最小值α=αmin=90-θ-β时,θ为相机垂直视场角的一半,即被拍摄物体投影到图片最底端时,v=vmax(vmax为相机CMOS或CCD图像传感器列坐标有效像素数),代入式(8)可得:
90-β-θ=a·vmax+b (9)
当αmin+2θ>90°,即θ>β时,此时相机上视角高于水平线,地平面无限远处,α无限接近于90°,此时v无限趋近于v0-tanβ*fy,fy为像素单位下相机的焦距,β为负值即相机逆时针旋转时亦同理,因此,代入式(8)可得:
90=a·(v0-tanβ·fy)+b (10)
当αmin+2θ<90°,即θ<β时,此时相机上视角低于水平线,地平面无限远处目标物成像角度α取最大值,αmax=αmin+2θ=90-β+θ时,即被拍摄物体投影到图片最高点时,v=0,代入式(8)可得:
90-β+θ=b (11)
根据针孔相机构造原理,一半的相机垂直视场角θ的正切值等于相机CMOS或CCD图像传感器边长的一半除以相机焦距,故可以计算出θ:
公式(12)中LCMOS为相机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在以像素为单位的图像坐标系中的坐标为(u0,v0);像平面坐标系和像素坐标系属于同一个平面上原点不同的两个坐标系;以相机中心点为坐标原点建立像空间坐标系S-XcYcZc,图像平面与光轴Zc轴垂直,和投影中心距离为f;将像空间坐标系点映射到投影平面上的过程为投影变换;另外,本文还为图像中N株待测量立木定义N组摄影测量坐标系Di-UiViWi和物方空间坐标系Ti-XiYiZi(i≤N),第i株立木对应的摄影测量坐标系原点Di位于该立木底部几何中心点,Vi轴竖直向上垂直于水平地面,Ui轴沿立木成像面方向水平向右平行于像平面,摄影测量坐标系沿立木树干方向平移1.3m后,旋转一定的角度可得到该株立木的物方空间坐标系,物方空间坐标系Yi轴沿立木树干方向竖直向上;
根据公式(14)深度提取模型计算出摄影测量坐标系下的立木深度值;
步骤四:测量单幅图像中多株立木胸径值
在针孔相机模型中,根据步骤4中立木胸径测量坐标***建立规则,研究成像点在不同坐标系中的转换关系,建立立木胸径测量模型,
其中,平移矩阵针对自然界中立木可能存在一定程度的倾斜角且树干为圆柱体的特征,定义立木胸径测量物方空间坐标系与摄影测量坐标系之间的旋转关系为以Vi轴为主轴的***,即以Vi轴为主轴旋转0°,然后绕Ui轴旋转角,最后绕Wi轴旋转Ψ角,因此旋转矩阵可表示为:
摄影测量坐标系经平移和旋转到像空间坐标系是一种刚体运动,摄影测量坐标系分别沿U,V,W轴方向平移tU,tV,tW距离(单位:毫米)后,从摄影测量坐标系旋转到像空间坐标系,其旋转角元素为(κ,β,ω),旋转过程为以V轴为主轴旋转κ角度后,绕新的U轴旋转β角度,最后绕W轴旋转ω角度,第i株立木的摄影测量坐标系中的任意点(Ui,Yi,Wi)在像空间坐标系下坐标(Xi,Yi,Zi)为:
式中,T表示由摄影测量坐标系到像空间坐标系的平移矩阵,R表示旋转矩阵:
[tU,tV,tW]为:
根据摄影测量坐标***建立规则,摄影测量坐标系绕V轴旋转角度κ等于180°,且通过移动端相机进行图片采集时绕W轴旋转角度ω通常约等于0°,即进行图片采集时,手机竖直置于相机三脚架上,因此,结合公式(15)和公式(17)求算物方空间坐标系原点在像空间坐标系中的坐标值(Xi0,Yi0,Zi0):
像空间坐标系中目标立木深度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’(x1,y1,f)、B’(x2,y2,f),像点A’、B’在自然场景中的同名物点为A、B,AB所在的直线与图像平面平行,
令点A在相机坐标系下的坐标为(X,Y,Z),点B的坐标为(X+Tx,Y,Z),线段CD过AB的中点垂直于立木树干,CD的距离即为立木胸径真实值;在针孔相机成像模型中,基于步骤三中相机镜头组畸变的矫正,图像中像点、相机光心和物点3点共线,则有:
结合公式(1)、(20)~(23),像点A’、B’的纵坐标像素值vDBH可表示为:
Y值相同且深度值Z相等的两点A’、B’的水平视差d可表示为:
d=x2-x1=f·Tx/Z
=(u2-u1)·dx (25)
其中,u1、u2分别是像点A’、B’的横坐标像素值,在已知相机内部参数的情况下,
结合像空间坐标系中待测立木深度值Z,可计算点AB的距离Tx:
因此,该株立木的胸径可表示为:
DBH=Tx·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)中的坐标为(u0,v0),(x,y)是实际图像中像点归一化的坐标,图像中任意像素在两个坐标系中满足如下关系:
相机坐标系中任一点Pc(Xc,Yc,Zc)投影到图像坐标系上为(xc,yc,f),图像坐标系平面与光轴z轴垂直,与原点距离为f,根据相似三角形原理可以得出:
引入所述改进的带有非线性畸变项的标定模型,包括由于镜头形状缺陷造成的径向畸变和由于光学***存在不同程度的偏心造成的切向畸变,径向畸变数学模型为:
其中r2=x2+y2,(x’,y’)为矫正后不含畸变项的理想线性相机坐标系的归一化坐标值,径向畸变值与图像点在图像中的位置有关,图像边缘处的径向畸变值较大,
切向畸变模型数学模型为:
其中包含k1、k2、k3、p1、p2共5个非线性畸变系数,由公式(3)、(4)得畸变矫正函数模型如下:
从世界坐标变换到相机坐标转换存在如下关系:
Pc=R·(PW-C)=R·PW+T (6)
结合式(1)~(6),用齐次坐标与矩阵形式可表示为:
Mint、Mext分别是相机标定内、外参数矩阵,其中相机内部参数包括图像中心点像素值u0、v0,fx、fy为x轴和y轴上的归一化焦距、通过Java结合OpenCV实现手机相机标定,获取手机相机所述的内部参数、相机镜头畸变参数和图像分辨率vmax、umax;
二、通过对新型标靶图像的采集,建立相机深度提取模型。现有的标靶为长宽相等的黑白棋盘格标靶。本发明的新型靶标与现有靶标的区别在于,设定标靶距离相机最近的一排方格大小为d*dmm,随后每排方格的宽度为固定,长度后一排比前一排增加值为
下式中xi为第i个角点到相机的实际距离,yi为各方格的长度,则相邻方格长度的差值Δdi为:
设各方格的计算长度与实际距离之间的关系为f(x),根据公式(28)可得:
经Pearson相关性分析,长度与实际距离之间呈极显著线性相关关系(p<0.01),相关系数r等于0.975,通过最小二乘法可以求算出f(x)的导数f’(x),
拍摄水平地面上物体时存在透视变换现象使得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算法找到其中的两个主要的模态γ1和γ2;根据边缘的方向,对于期望的梯度强度构造一个模板T。(*表示互相关操作符)和角点相似度的乘积作为角点分值,然后用阈值进行判断就得到初始角点。
2)对角点的位置和方向进行亚像素级精细提取;
运用OpenCV中的cornerSubPix()函数进行亚像素级角点定位,将角点定位到子像素,从而取得亚像素级别的角点检测效果;为细化边缘方向向量,根据图像梯度值最小化其标准离差率:
其中是相邻像素集,与模块i的梯度值mi=[cos(γi)sin(γi)]T相匹配。(求算方案依据文献Geiger A,Moosmann F,Caretal.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,Caret 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)=Ecorners(y)+Estruct(x,y) (32)
其中,Ecorners是当前棋盘角点总数的负值,Estruct是两个相邻角点和预测角点的匹配程度;通过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=vmax(vmax为相机CMOS或CCD图像传感器列坐标有效像素数),代入式(8)可得:
90-β-θ=a·vmax+b (9)
当αmin+-2θ>90°,即θ>β时,此时相机上视角高于水平线,相机拍摄投影几何模型如图4,地平面无限远处,α无限接近于90°,此时v无限趋近于v0-tanβ*fy,fy为像素单位下相机的焦距,β为负值即相机逆时针旋转时亦同理,因此,代入式(8)可得:
90=a·(v0-tan β·fy)+b (10)
当αmin+2θ<90°,即θ<β时,此时相机上视角低于水平线,相机拍摄投影几何模型如图5,地平面无限远处目标物成像角度α取最大值,αmax=αmin+2θ=90-β+θ时,即被拍摄物体投影到图片最高点时,v=0,代入式(8)可得:
90-β+θ=b (11)
根据针孔相机构造原理,一半的相机垂直视场角θ的正切值等于相机CMOS或CCD图像传感器边长的一半除以相机焦距,故可以计算出θ:
公式(12)中LCMOs为相机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在以像素为单位的图像坐标系中的坐标为(u0,v0);像平面坐标系和像素坐标系属于同一个平面上原点不同的两个坐标系。以相机中心点为坐标原点建立像空间坐标系S-XcYcZc,图像平面与光轴Zc轴垂直,和投影中心距离为f。将像空间坐标系点映射到投影平面上的过程为投影变换;另外,本文还为图像中N株待测量立木定义N组摄影测量坐标系Di-UiViWi和物方空间坐标系Ti-XiYiZi(i≤N),第i株立木对应的摄影测量坐标系原点Di位于该立木底部几何中心点,Vi轴竖直向上垂直于水平地面,Ui轴沿立木成像面方向水平向右平行于像平面,摄影测量坐标系沿立木树干方向平移1.3m后,旋转一定的角度可得到该株立木的物方空间坐标系,物方空间坐标系Yi轴沿立木树干方向竖直向上。
四、测量单幅图像中多株立木胸径值。
在针孔相机模型中,根据步骤三中立木胸径测量坐标***建立规则,研究成像点在不同坐标系中的转换关系,建立立木胸径测量模型,
其中,平移矩阵针对自然界中立木可能存在一定程度的倾斜角且树干为圆柱体的特征,定义立木胸径测量物方空间坐标系与摄影测量坐标系之间的旋转关系为以Vi轴为主轴的***,即以Vi轴为主轴旋转0°,然后绕Ui轴旋转角,最后绕Wi轴旋转Ψ角,因此旋转矩阵可表示为:
基于步骤三中提取的立木树干轮廓,可通过提取图像中立木轮廓的最小外接矩形获得旋转角Ψ。根据自然环境下的立木特点,本文采用优化的最小外接矩形直接计算方法提取最小外接矩形,该算法首先按照最小外接矩形直接计算方法获取立木树干轮廓区域外接矩形,此时该矩形倾斜角为0°;然后以外接矩形中心为中心点旋转该外接矩形,直至外接矩形面积取最小值时记录旋转角,从而找到最优矩形姿态,该旋转角即为旋转角Ψ。
摄影测量坐标系经平移和旋转到像空间坐标系是一种刚体运动,摄影测量坐标系分别沿U,V,W轴方向平移tU,tV,tW距离(单位:毫米)后,从摄影测量坐标系旋转到像空间坐标系,其旋转角元素为(κ,β,ω),旋转过程为以V轴为主轴旋转κ角度后,绕新的U轴旋转β角度,最后绕W轴旋转ω角度,第i株立木的摄影测量坐标系中的任意点(Ui,Vi,Wi)在像空间坐标系下坐标(Xi,Yi,Zi)为:
式中,T表示由摄影测量坐标系到像空间坐标系的平移矩阵,R表示旋转矩阵:
[tU,tV,tW]为:
根据摄影测量坐标***建立规则,摄影测量坐标系绕V轴旋转角度κ等于180°,且通过移动端相机进行图片采集时绕W轴旋转角度ω通常约等于0°,即进行图片采集时,手机竖直置于相机三脚架上,因此,结合公式(15)和公式(17)求算物方空间坐标系原点在像空间坐标系中的坐标值(Xi0,Yi0,Zi0):
像空间坐标系中目标立木深度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’(x1,y1,f)、B’(x2,y2,f),像点A’、B’在自然场景中的同名物点为A、B,AB所在的直线与图像平面平行,
令点A在相机坐标系下的坐标为(X,Y,Z),点B的坐标为(X+Tx,Y,Z),线段CD过AB的中点垂直于立木树干,CD的距离即为立木胸径真实值;在针孔相机成像模型中,基于步骤三中相机镜头组畸变的矫正,图像中像点、相机光心和物点3点共线,则有:
结合公式(1)、(20)~(23),像点A’、B’的纵坐标像素值vDBH可表示为:
Y值相同且深度值Z相等的两点A’、B’的水平视差d可表示为:
d=x2-x1=f·Tx/Z
=(u2-u1)·dx (25)
其中,u1、u2分别是像点A’、B’的横坐标像素值,在已知相机内部参数的情况下,
结合像空间坐标系中待测立木深度值Z,可计算点AB的距离Tx:
因此,该株立木的胸径可表示为:
DBH=Tx·cosψ (27)。
实施例1
下面以小米3(MI 3)手机为例,具体说明本发明的基于单目视觉的优化的深度提取和被动测距方法。
一、对手机相机进行标定,获取相机内部参数和图像分辨率
使用行列数为8*9尺寸大小为20*20的棋盘格标定板作为相机标定的实验材料,通过小米3手机相机采集20张不同角度的标定板图片,利用OpenCV根据上述改进的带有非线性畸变项的相机标定模型对小米3(MI 3)手机相机进行标定,
首先使用fin()函数读取标定板图片,并通过.cols、.rows获取第一张图片的图像分辨率;然后通过find4QuadComerSubpix()函数提取标定板图片中亚像素级角点,并用drawChessboardCorners()函数标记角点;调用calibrateCamera()函数对相机进行标定,将得到的摄像机内外参数用于对空间的三维点进行重新投影计算,得到新的投影点,计算新的投影点和旧的投影点之间的误差;输出并保存相机内参矩阵和畸变参数,
标定所得相机内部参数为:fx=3486.5637,u0=1569.0383,fy=3497.4652,v0=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)中的坐标为(u0,v0),(x,y)是实际图像中像点归一化的坐标,图像中任意像素在两个坐标系中满足如下关系:
相机坐标系中任一点Pc(Xc,Yc,Zc)投影到图像坐标系上为(xc,yc,f),图像坐标系平面与光轴z轴垂直,与原点距离为f,根据相似三角形原理可以得出:
引入所述改进的带有非线性畸变项的标定模型,包括由于镜头形状缺陷造成的径向畸变和由于光学***存在不同程度的偏心造成的切向畸变,径向畸变数学模型为:
其中r2=x2+y2,(x’,y’)为矫正后不含畸变项的理想线性相机坐标系的归一化坐标值,径向畸变值与图像点在图像中的位置有关,图像边缘处的径向畸变值较大,
切向畸变数学模型为:
其中包含k1、k2、k3、p1、p2共5个非线性畸变系数,由公式(3)、(4)得畸变矫正函数模型如下:
从世界坐标变换到相机坐标转换存在如下关系:
Pc=R·(PW-C)=R·PW+T (6)
公式(6)中T表示由摄影测量坐标系到像空间坐标系的平移矩阵,R表示旋转矩阵,结合式(1)~(6),用齐次坐标与矩阵形式可表示为:
Mint、Mext分别是相机标定内、外参数矩阵,其中相机内部参数包括图像中心点像素值u0、v0,fx、fy为x轴和y轴上的归一化焦距、通过Java结合OpenCV实现手机相机标定,获取手机相机所述的内部参数、相机镜头畸变参数和图像分辨率vmax、umax;
步骤二:建立深度提取模型
根据目标物成像角度α与纵坐标像素值v之间的线性关系设定抽象函数,建立含目标物成像角度α、纵坐标像素值v和相机旋转角β三个参数空间关系模型,即α=F(v,β),
不同型号的设备和相机旋转角度下,被拍摄物体纵坐标像素值与成像角度均呈极显著负线性相关关系,且该线性关系的斜率与截距有所不同,故设:
α=F(v,β)=a·v+b (8)
其中参数a、b均与相机型号和相机旋转角度有关,
当α取最小值α=αmin=90-θ-β时,θ为相机垂直视场角的一半,即被拍摄物体投影到图片最底端时,v=vmax,vmax为相机CMOS或CCD图像传感器列坐标有效像素数,代入式(8)可得:
90-β-θ=a·vmax+b (9)
当αmin+2θ>90°,即θ>β时,此时相机上视角高于水平线,地平面无限远处,α无限接近于90°,此时v无限趋近于v0-tanβ*fy,fy为像素单位下相机的焦距,β为负值即相机逆时针旋转时亦同理,因此,代入式(8)可得:
90=a·(v0-tanβ·fy)+b (10)
当αmin+2θ<90°,即θ<β时,此时相机上视角低于水平线,地平面无限远处目标物成像角度α取最大值,αmax=αmin+2θ=90-β+θ时,即被拍摄物体投影到图片最高点时,v=0,代入式(8)可得:
90-β+θ=b (11)
根据针孔相机构造原理,一半的相机垂直视场角θ的正切值等于相机CMOS或CCD图像传感器边长的一半除以相机焦距,故可以计算出θ:
公式(12)中LCMOS为相机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在以像素为单位的图像坐标系中的坐标为(u0,v0);像平面坐标系和像素坐标系属于同一个平面上原点不同的两个坐标系;以相机中心点为坐标原点建立像空间坐标系S-XcYcZc,图像平面与光轴Zc轴垂直,和投影中心的距离为f;将像空间坐标系点映射到投影平面上的过程为投影变换;另外,为图像中N株待测量立木定义N组摄影测量坐标系Di-UiViWi和物方空间坐标系Ti-XiYiZi(i≤N),第i株立木对应的摄影测量坐标系原点Di位于该立木底部几何中心点,Vi轴竖直向上垂直于水平地面,Ui轴沿立木成像面方向水平向右平行于像平面,摄影测量坐标系沿立木树干方向平移1.3m后,旋转一定的角度可得到该株立木的物方空间坐标系,物方空间坐标系Yi轴沿立木树干方向竖直向上;
根据公式(14)深度提取模型计算出摄影测量坐标系下的立木深度值;
步骤四:测量单幅图像中多株立木胸径值
在针孔相机模型中,根据步骤4中立木胸径测量坐标***建立规则,研究成像点在不同坐标系中的转换关系,建立立木胸径测量模型,
摄影测量坐标系经平移和旋转到像空间坐标系是一种刚体运动,摄影测量坐标系分别沿U,V,W轴方向平移tU,tV,tW距离后,单位:毫米,从摄影测量坐标系旋转到像空间坐标系,其旋转角元素为(κ,β,ω),旋转过程为以V轴为主轴旋转κ角度后,绕新的U轴旋转β角度,最后绕W轴旋转ω角度,第i株立木的摄影测量坐标系中的任意点(Ui,Vi,Wi)在像空间坐标系下坐标(Xi,Yi,Zi)为:
式中,T表示由摄影测量坐标系到像空间坐标系的平移矩阵,R表示旋转矩阵:
[tU,tV,tW]为:
根据摄影测量坐标***建立规则,摄影测量坐标系绕V轴旋转角度κ等于180°,且通过移动端相机进行图片采集时绕W轴旋转角度ω等于0°,即进行图片采集时,手机竖直置于相机三脚架上,因此,结合公式(15)和公式(17)求算物方空间坐标系原点在像空间坐标系中的坐标值(Xi0,Yi0,Zi0):
像空间坐标系中目标立木深度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’(x1,y1,f)、B’(x2,y2,f),像点A’、B’在自然场景中的同名物点为A、B,AB所在的直线与图像平面平行,
令点A在相机坐标系下的坐标为(X,Y,Z),点B的坐标为(X+Tx,Y,Z),线段CD过AB的中点垂直于立木树干,CD的距离即为立木胸径真实值;在针孔相机成像模型中,基于步骤三中相机镜头组畸变的矫正,图像中像点、相机光心和物点3点共线,则有:
结合公式(1)、(20)~(23),像点A’、B’的纵坐标像素值vDBH可表示为:
Y值相同且深度值Z相等的两点A’、B’的水平视差d可表示为:
其中,u1、u2分别是像点A’、B’的横坐标像素值,在已知相机内部参数的情况下,结合像空间坐标系中待测立木深度值Z,可计算点A、B之间的距离Tx:
因此,该株立木的胸径可表示为:
DBH=Tx·cosψ (27)。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810913273.3A CN109269430B (zh) | 2018-08-12 | 2018-08-12 | 基于深度提取模型的多株立木胸径被动测量方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810913273.3A CN109269430B (zh) | 2018-08-12 | 2018-08-12 | 基于深度提取模型的多株立木胸径被动测量方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109269430A CN109269430A (zh) | 2019-01-25 |
CN109269430B true CN109269430B (zh) | 2020-10-09 |
Family
ID=65153657
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810913273.3A Active CN109269430B (zh) | 2018-08-12 | 2018-08-12 | 基于深度提取模型的多株立木胸径被动测量方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109269430B (zh) |
Families Citing this family (22)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109816779A (zh) * | 2019-01-30 | 2019-05-28 | 北京林业大学 | 一种使用智能手机重建人工林森林模型获取单木参数的方法 |
CN109919990B (zh) * | 2019-02-19 | 2021-03-16 | 北京工业大学 | 利用深度感知网络和视差遥感影像进行森林高度预测方法 |
CN110030928A (zh) * | 2019-04-11 | 2019-07-19 | 接楚添 | 基于计算机视觉的空间物体定位和测量的方法和*** |
CN110097516B (zh) * | 2019-04-25 | 2021-02-12 | 上海交通大学 | 内孔壁面图像畸变纠正方法、***及介质 |
CN110345921B (zh) * | 2019-06-12 | 2020-06-26 | 中国农业大学 | 立体视场视觉测量及垂轴像差和轴向像差校正方法及*** |
CN110378901B (zh) * | 2019-08-16 | 2023-08-08 | 扬州大学 | 一种基于深度相机的球形花木中心测量与修剪控制方法 |
CN110617772A (zh) * | 2019-10-09 | 2019-12-27 | 南京天创电子技术有限公司 | 一种非接触式线径测量装置及方法 |
CN111161305A (zh) * | 2019-12-18 | 2020-05-15 | 任子行网络技术股份有限公司 | 无人机智能识别追踪方法及*** |
CN113538477B (zh) * | 2020-04-14 | 2023-08-29 | 北京达佳互联信息技术有限公司 | 平面位姿的获取方法、装置、电子设备及存储介质 |
CN111504208A (zh) * | 2020-05-21 | 2020-08-07 | 南京鸿亦沄智能科技有限公司 | 一种基于计算机视觉的非接触式树木胸径测量方法及*** |
CN111612720B (zh) * | 2020-05-21 | 2023-11-07 | 烟台艾睿光电科技有限公司 | 一种广角红外图像优化方法、***及相关组件 |
CN112001654B (zh) * | 2020-08-31 | 2024-04-30 | 上海市园林科学规划研究院 | 一种城市搬迁地绿化树种地下生物量模型的构建方法 |
CN111784622B (zh) * | 2020-09-07 | 2021-01-26 | 成都纵横自动化技术股份有限公司 | 一种基于无人机单目倾斜的图像拼接方法及相关装置 |
CN112556589B (zh) * | 2020-12-03 | 2022-01-07 | 厦门大学 | 一种高度测量方法和高度缺陷判定方法 |
CN113484867B (zh) * | 2021-06-25 | 2023-10-20 | 山东航天电子技术研究所 | 一种基于成像声呐封闭空间下鱼群密度探测方法 |
CN113776483B (zh) * | 2021-08-13 | 2023-11-14 | 华电电力科学研究院有限公司 | 一种不圆度测量装置及测量方法 |
CN114279431A (zh) * | 2021-11-22 | 2022-04-05 | 北京林业大学 | 一种用于密郁闭度森林的林木位置测量方法及其装置 |
CN114529616B (zh) * | 2022-04-22 | 2022-07-26 | 武汉精视遥测科技有限公司 | 基于内壁刻度的广角镜头参数标定方法、***及计算机 |
CN116977403B (zh) * | 2023-09-20 | 2023-12-22 | 山东科技大学 | 一种基于双目视觉的薄膜生产幅宽检测和控制方法 |
CN117607063B (zh) * | 2024-01-24 | 2024-04-19 | 中国科学院地理科学与资源研究所 | 一种基于无人机的森林垂直结构参数测量***和方法 |
CN118031847B (zh) * | 2024-04-10 | 2024-06-04 | 甘肃祁连山国家级自然保护区管护中心大河口自然保护站(大熊猫祁连山国家公园甘肃省管理局张掖分局大河口保护站) | 一种周期式树木生长直径监测方法和*** |
CN118111356B (zh) * | 2024-04-28 | 2024-07-16 | 慧诺云谱(海南)科技有限公司 | 植物表型图像测量方法、装置、电子设备及存储介质 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105513078A (zh) * | 2015-12-15 | 2016-04-20 | 浙江农林大学 | 基于图像的立木信息采集方法和装置 |
CN105571561A (zh) * | 2015-12-15 | 2016-05-11 | 浙江农林大学 | 立木信息采集方法和装置 |
CN107084672A (zh) * | 2017-05-08 | 2017-08-22 | 北京林业大学 | 一种用于树木胸径测量的图像采集装置、***和方法 |
CN107607053A (zh) * | 2017-09-20 | 2018-01-19 | 浙江农林大学 | 一种基于机器视觉和三维重建技术的立木胸径测量方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US9198363B2 (en) * | 2012-12-12 | 2015-12-01 | The Boeing Company | Tree metrology system |
-
2018
- 2018-08-12 CN CN201810913273.3A patent/CN109269430B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105513078A (zh) * | 2015-12-15 | 2016-04-20 | 浙江农林大学 | 基于图像的立木信息采集方法和装置 |
CN105571561A (zh) * | 2015-12-15 | 2016-05-11 | 浙江农林大学 | 立木信息采集方法和装置 |
CN107084672A (zh) * | 2017-05-08 | 2017-08-22 | 北京林业大学 | 一种用于树木胸径测量的图像采集装置、***和方法 |
CN107607053A (zh) * | 2017-09-20 | 2018-01-19 | 浙江农林大学 | 一种基于机器视觉和三维重建技术的立木胸径测量方法 |
Non-Patent Citations (4)
Title |
---|
基于Android平台的测树***研究与实现;周克瑜 等;《南京林业大学学报( 自然科学版)》;20160731;第40卷(第4期);第95-100页 * |
基于单目视觉的近景摄影立木测量方法研究;周洋;《中国优秀硕士学位论文全文数据库 信息科技辑》;20180415;第I138-2507页 * |
基于数码相机的样木胸径获取方法;黄晓东 等;《农业机械学报》;20150930;第46卷(第9期);第266-272页 * |
基于智能手机图像分析的树木胸径测量研究;沈亚峰 等;《江苏林业科技》;20170228;第44卷(第1期);第28-33页 * |
Also Published As
Publication number | Publication date |
---|---|
CN109269430A (zh) | 2019-01-25 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109269430B (zh) | 基于深度提取模型的多株立木胸径被动测量方法 | |
CN109035320B (zh) | 基于单目视觉的深度提取方法 | |
CN109146980B (zh) | 基于单目视觉的优化的深度提取和被动测距方法 | |
CN112396664B (zh) | 一种单目摄像机与三维激光雷达联合标定及在线优化方法 | |
CN110118528B (zh) | 一种基于棋盘靶标的线结构光标定方法 | |
Wu et al. | Passive measurement method of tree diameter at breast height using a smartphone | |
WO2019105044A1 (zh) | 一种镜头畸变矫正和特征提取的方法及*** | |
CN109255818B (zh) | 一种新型标靶及其亚像素级角点的提取方法 | |
CN107977996B (zh) | 基于靶标标定定位模型的空间目标定位方法 | |
CN108171715B (zh) | 一种图像分割方法及装置 | |
CN112132908B (zh) | 一种基于智能检测技术的相机外参数标定方法及设备 | |
CN110889829A (zh) | 一种基于鱼眼镜头的单目测距方法 | |
Beekmans et al. | Cloud photogrammetry with dense stereo for fisheye cameras | |
CN109859137B (zh) | 一种广角相机非规则畸变全域校正方法 | |
CN109448043A (zh) | 平面约束下的立木高度提取方法 | |
CN106971408A (zh) | 一种基于时空转换思想的摄像机标定方法 | |
CN114396875B (zh) | 一种基于深度相机垂直拍摄的长方形包裹体积测量方法 | |
CN108320310B (zh) | 基于图像序列的空间目标三维姿态估计方法 | |
Xinmei et al. | Passive measurement method of tree height and crown diameter using a smartphone | |
CN115097421A (zh) | 一种相机-激光雷达外参标定装置及方法 | |
CN115201883A (zh) | 一种运动目标视频定位测速***及方法 | |
Abanay et al. | A calibration method of 2D LIDAR-Visual sensors embedded on an agricultural robot | |
CN116129037A (zh) | 视触觉传感器及其三维重建方法、***、设备及存储介质 | |
CN114998448A (zh) | 一种多约束双目鱼眼相机标定与空间点定位的方法 | |
CN107941241B (zh) | 一种用于航空摄影测量质量评价的分辨率板及其使用方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |