CN116309844A - 一种基于无人机单张航空图片的三维测量方法 - Google Patents

Abstract

本发明提供一种基于无人机单张航空图片的三维测量方法，涉及计算机视觉技术领域。本发明通过无人机相机获得待检测图像，构造小孔成像模型，将二维空间坐标转化至三维空间坐标，将已知的若干控制点的坐标系进行数据转换，进行相机标定获得相机内外参数，矫正待处理图像的畸变问题，并将待处理图像进行灰度化预处理，构造了一种端到端的单目三维物体测量模型，将处理后图像输入至三维物体测量模型中，优化了视频中三维测量技术，不限于目标物体的类型，大幅度提高了检测正确率。并且能够以无监督的方式从简单的视频中学习，容易迁移到其他数据集上。

Description

一种基于无人机单张航空图片的三维测量方法

技术领域

本发明涉及计算机视觉技术领域，尤其涉及一种基于无人机单张航空图片的三维测量方法。

背景技术

计算机视觉是指利用计算机来实现人类的相关视觉功能，计算机视觉技术作为一门新兴学科，它涉及到信号处理、人工智能、图像处理和模式识别等多个领域。近年来，随着人工智能技术以及数字图像处理技术的飞速发展，计算机视觉技术已经被广泛应用于各个行业如航空航天等。三维视觉测量的开创性工作是从六十年代中期开始的，美国麻省理工学院的Robert完成了三维景物分析开创性工作，把对二维图像推广到了三维图像，三维视觉技术由此诞生。在之后的二十年中，三维视觉测量迅速发展成为一门学科。D.Marr创立的立体视觉计算理论对三维视觉技术的发展产生了深远影响。三维视觉测量在军事、航空航天和目标识别等领域应用变得非常广泛。当前对于三维视觉测量方法的性能在不断提高，其理论和方法也处于不断发展和完善之中。

当代所生产的无人机具有成本费用低、灵活度高、容易操作、体积小等特点，目前被应用在了更多多元化的场景中。无人机航拍是极具创新性的成像方式，它的智能分析处理可以快速高效地提取所拍摄图像中的物体信息。无人机航空拍摄系统可用如数码相机、磁测仪或激光扫描仪等机载遥感设备，获取指定区域范围内的图像，再进行相关图像数据的处理。恢复一张图像中某一个物体的三维信息是无人机航空拍摄系统一个非常重要的研究方向，而物体的三维信息研究的一个非常重要的基础环节就是测量目标物体的长宽高。三维测量信息图像中的目标分割、目标行为分析与场景动态分析，对以无人机为背景的其他技术具有重要的辅助。

发明内容

针对现有技术的不足，本发明提供一种基于无人机单张航空图片的三维测量方法。本发明能够提高测量精度，解决了现有三维物体测量需要人工参与的问题。该方法不限制于检测物体的类型，不需要考虑物体类别的差别性。

一种基于无人机单张航空图片的三维测量方法，具体包括以下步骤：

步骤1：利用无人机的相机进行单张图像拍摄并记录无人机拍摄时的飞行高度s。

步骤2：构造光学二维空间小孔成像模型，建立物距与焦距的关系，来确定空间中目标三维物体信息与图像二维成像之间的点面关系。

步骤2.1：建立物距与像距和焦距的关系表达式

其中焦距设为f，物体到镜头光心称为物距设为u，成像平面到镜头光心称为像距设为v；

步骤2.2：若己知物距v，设x₁、y₁为物体在三维空间的长和宽，x₂、y₂分别为物体摄像头所拍摄的二维图像的长和宽，则求得物距与焦距的关系公式，即光学二维空间小孔成像模型为：

步骤3：将二维空间坐标转化至三维空间坐标，将已知的若干控制点的坐标系进行数据转换。

步骤3.1：建立像素坐标系：在图像上定义直角坐标系，每一像素的坐标点(a,b)表示其坐标值，a和b分别是这个像素在数组中的行数与列数；

步骤3.2：建立图像坐标系：X轴和Y轴分别平行于二维图像的两条垂直边，用(x,y)表示像素坐标值，设(a₀,b₀)为二维图像的中心，相机中感光器件的尺寸为x*y，像素坐标系与图像坐标系之间的关系表示为：

其中a-a₀和b-b₀分别表示某以像素点坐标数值与中心点坐标数值的X轴差值与Y轴差值；

将上式改写成矩阵形式得：

将偏移项带入进乘积项中，构造一个三维矩阵乘法形式，转化为齐次坐标形式为：

步骤3.3：建立相机坐标系：相机坐标系中心点设为O_x，图像坐标系中心点设为O_t，O_x与的O_t相交点的延长线为Z轴，即步骤3.2所述图像坐标系的延伸，Z轴的方向服从物理学中的右手规则，根据欧拉旋转定理知三维空间的旋转顺序不可逆性，X和Y轴在一平面内且平行，Z轴与X和Y轴所称平面不平行，三维空间中的点P在二维平面的投影点设为P'，再通过透视投影关系算出坐标，设点P在三维图像中的坐标为(x,y,z)，点P在二维平面投影点P'的坐标(x',y')用步骤2.3中物距与焦距的关系公式及相似三角形比例关系求出透视投影关系公式：

用齐次坐标与矩阵表示前述的透视投影关系公式为:

步骤3.4：建立世界坐标系；

所述世界坐标系是描述摄像头所在位置的一个基准坐标系，由X轴、Y轴、Z轴组成；摄像机坐标系与世界坐标系之间的关系在用旋转矩阵R与偏移向量T来描述，即物体只需要被旋转和平移；若在三维空间中存在一点P，其在世界坐标系下与像机坐标系下的齐次坐标分别是(x_s,y_s,z_s,1)与(x',y',z')，则存在如下关系：

建立像素坐标系与世界坐标系之间的变换关系：

式中，其中m及n是相机内部参数；

建立投影矩阵M为：

其中a₀、b₀为二维图像的中心点坐标，只与摄像机内部结构有关，为摄像机内部参数，是摄像头感光芯片安装时与理论中心的水平与垂直差值，

由摄像头相对于世界坐标系的方位决定，是相机相对于世界坐标系的平移变换关系，即摄像机外部参数；

步骤3.5：从图像坐标系转换到世界坐标系，己知摄像机的内外参数，获得投影矩阵M，对三维空间中的点P＝(x,y,z)，最终需要得到其在世界坐标系下的坐标P_s＝(x_s,y_s,z_s,1)，求出在图像坐标系点P的位置(x,y)，先从像素坐标系变换到相机坐标系，再由相机坐标系变换到世界坐标系，具体为：

其中，s为从世界坐标系转换到图像坐标系的尺度因子，t为平移矩阵。

步骤4：进行无人机相机标定，求出相机参数；

步骤4.1：根据步骤3.4中建立的像素坐标系与世界坐标系之间的变换关系，使用张氏标定法，每一个角点的像素坐标设为(a,b)，每一个角点在世界坐标系下的物理坐标(A,B,C)＝0来进行相机的标定，获得相机的内参矩阵和径向畸变参数。

步骤4.2：将2D空间中的图像点设为n＝[a,b]^T，3D空间中的图像点设为N＝[x,y,z]，齐次坐标分别为

三维空间坐标映射到二维图像坐标转换为/>

其中s为世界坐标系到图像坐标系的尺度因子，α和β分别为焦距与像素横纵比，γ为径向畸变参数，/>

为相机内参矩阵；

步骤4.3：求解相机内参，三维空间坐标到二维图像坐标的映射为：H＝A[r1 r2t]，其中r1和r2为两个相互正交的旋转向量，A＝s*M，解Homographic矩阵进行单位性转换，将平面的点映射到另一个二维平面的投影，定义旋转矩阵的第i列为r_i，棋盘格位于Z轴上，表示为：s*[a b 1]^T＝A[r₁r₂ r₃ t]*c*[x y 0 1]^T。Homography有8个自由度，r₁和r₂正交且模相等，得到约束式为：

式中，h₁、h₂分别表示矩阵H的第一列与第二列；

构造一个矩阵B，使得

且约束条件为：

式中，v代表图像中棋盘格的角点纵坐标，因棋盘格是长方形，所以只需要得知这个棋盘格的其中三个角点坐标就可以确定剩余一个角点坐标。v11、v12、v22分别表示棋盘格的左上方、右上方、右下方角点的纵坐标，b表示用奇异值分解最小奇异值对应的特征向量。

内参表达式如下：

其中，B₁₁、B₁₂、B₁₃、B₂₂、B₂₃、B₃₃公式表示矩阵B中不同位置元素，v₀、u₀表示中心点的实际位置坐标，λ为λ＝s^-1即尺度因子的逆矩阵，γ为径向畸变参数；

步骤4.4：求解相机外参：

r₁＝λM^-1h₁

r₂＝λM^-1h₁

r₃＝h₁*h₂

t＝λM^-1h₃

其中，r₁、r₂、r₃分别表示三个旋转向量，h₁、h₂、h₃分别表示矩阵H的第一列、第二列和第三列；

根据旋转矩阵的性质得：||r₁||＝||λM^-1h₁||＝1即

用h₁求解推导出r₁、r₂、r₃和t。

步骤5：解决图像畸变问题，对图像进行矫正；

步骤5.1：为解决因相机摄像头前方的透镜对光的传播路径的影响，对图像进行径向畸变矫正。

对于图像中某一点为点P，设x_j为发生径向畸变后的该点横坐标，y_j为发生径向畸变后的该点纵坐标；x为校正后的该点横坐标，y为校正后的该点纵坐标，归一化极坐标表示为[r,θ]，其中r为矢量长度，t是矢量的水平夹角，则径向畸变的矫正公式为：

x_j＝x*(1+k₁r²+k₃r⁶+k₂r⁴)

y_j＝y*(1+k₁r²+k₃r⁶+k₂r⁴)

其中，k₁、k₂、k₃是径向畸变参数；

步骤5.2：为解决因在机械组装中的误差产生的影响，对图像进行切向畸变矫正。

对于图像中某一点为点P，设x₀为发生横向畸变后的该点横坐标，y₀为发生横向畸变后的该点纵坐标；x为校正后的该点横坐标，y为校正后的该点纵坐标；用另外两个参数p₁、p₂进行纠正，则切向畸变的矫正公式为：

x₀＝x+2p₁y+p₂(r²+2x²)

y₀＝y+2p₂x+p₁(r²+2y²)

步骤5.3：联合步骤5.1和步骤5.2中的径向畸变、切向畸变矫正公式，对于相机坐标系中的点P[x,y,z]，寻找它对应点的像素坐标的过程时将三维空间点投影到归一化图像平面，该点的像素坐标为：u＝f_x*x+c_x，v＝f_y*y+c_y。其中，u为像素点的横坐标，v为像素点的纵坐标，c_x为图像像素中心点的横坐标，c_y为图像像素中心点的横坐标。

步骤6：将图像进行预处理；

步骤6.1：直方图是在灰度范围[0,...,k,...,225]区间的功能函数，列出原图像的灰度级j，统计出原图像中每个灰度级所拥有的像素个数n_j；

步骤6.2：计算出原图像的各灰度级的灰度分布概率

其中n为直方图中的像素总个数，再计算出原图像的各灰度累积分布概率/>

其中k表示灰度级，，求出均衡化后的灰度级i＝INT[(L-1)S(i)+0.5]，INT()表示取整函数、L表示总灰度级数目、S(i)表示在i级的离散灰度，将(L-1)*p(k)进行四舍五入计算。

步骤6.3：确定原图像的灰度级与均衡化后的灰度级的映射关系j→i，再计算变换后的图像直方图

步骤7：设计一种端到端的单目三维物体测量模型，将预处理后的航拍图像输入到单目三维物体测量模型中，完成航空图片的三维测量。

步骤7.1：建立单目三维物体测量模型，所述单目三维物体测量模型包括backbone、关键点检测分支和3D框回归分支；其中所述backbone层使用DLA-34结构来提取特征，将所获取的二维的RGB图像输入至backbone层中，通过DLA-34输出特征图，其长和宽分别为输入图像的1/4，通道数为256；在backbone中加入GN层，GN层将channel分组，在每个group中做LN，再做归一化处理；

所述二维RGB图像中，设W代表图像宽度，H代表图像高度；用(d,w,l,x,y,z,θ)表示目标三维图像，其中(d,w,l)的分别表示目标三维图像的高、宽、长，(x,y,z)表示目标三维图像在相机坐标系中的中心坐标点，θ表示拍摄图像时摄像头角度，拍摄图像时的飞行高度为s，相机内参设为K，

步骤7.2：将步骤7.1中backbone层输出的特征图输入关键点检测分支网络层，所述关键点检测分支网络层用于寻找目标的3D中心点在二维图像中的投影点位置；在关键点检测分支网络层寻找输入图像中的关键点并标记绿色，将目标三维图像的中心点在二维图像平面的映射点[x_cy_c]^T设为关键点，目标三维图像的中心点坐标设为[x y z]^T，用坐标[x_b,1 y_b,1]^T、[x_b,2 y_b,2]^T、[x_b,3 y_b,3]^T、[x_b,4 y_b,4]^T、[x_b,5 y_b,5]^T、[x_b,6 y_b,6]^T、[x_b,7 y_b,7]^T、[x_b,8 y_b,8]^T分别表示目标三维图像的8个顶点在图像中的位置，并用

确定高斯函数的标准差，其中/>

表示为表示3D框的最小长方形。

步骤7.3：将步骤7.1中backbone层输出的特征图像及步骤7.2中所得得到的目标3D中心点输入至3D框回归分支网络层中，设δ_z表示深度偏移，δ_xc和δ_yc为量化误差，δ_d、δ_w、δ_l分别表示高、宽、长的残差，通过深度偏移求得目标的深度值z＝μ_z+δ_zσ_z，μ_z、σ_z分别尺度因子和平移因子，物体的目标三维图像的中心点坐标表示为：

再求目标物体的长、宽、高表达式为：

在式中，

分别为公开的KITTI数据集中计算出的平均长度、宽度、高度，e^δl、e^δw、e^δd分别为长宽高的残差值。

步骤7.5：根据角度和目标在三维空间中的位置求得摄像头角度

α是预测的观测角，则拍摄该物体时，无人机距离目标物体的直线距离m为/>

通过将航向角矩阵R绕x轴旋转得到矩阵R_θ，根据物体的长宽高[l w d]^T和中心点坐标[x y z]^T得到目标三维物体的8个角点坐标矩阵表示为：

采用上述技术方案所产生的有益效果在于：

本发明提供一种基于无人机单张航空图片的三维测量方法，本方法本申请提供了一种根据视频中给定的图像进行分布式三维测量的方法，能够提高测量精度，解决了现有三维物体测量需要人工参与的问题。该方法不限制于检测物体的类型，不需要考虑物体类别的差别性。

附图说明

图1为本发明实施例提供的单张图像三维测量方法中方案总体步骤示意图；

图2为本发明实施例提供的单张图像三维测量方法中灰度级处理步骤示意图；

图3为本发明实施例提供的单张图像三维测量方法中将预处理后的图像输入三维物体测量模型步骤示意图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

一种基于无人机单张航空图片的三维测量方法，如图1所示，具体包括以下步骤：

步骤1：利用无人机的相机进行单张图像拍摄并记录无人机拍摄时的飞行高度s。

本实施例中利用安装相机无人机进行图像拍摄，且无人机由控制端操控。无人机在飞行过程中其摄像头将当前环境画面实时直播至控制端，控制端找到目标物体后，无人机利用其拍摄设备进行图像拍摄，不局限于拍摄设备，同时需要记录无人机拍摄目标物体时的飞行高度为s。

步骤2：构造光学二维空间小孔成像模型，建立物距与焦距的关系，来确定空间中目标三维物体信息与图像二维成像之间的点面关系。

步骤2.1：建立物距与像距和焦距的关系表达式

其中焦距设为f，物体到镜头光心称为物距设为u，成像平面到镜头光心称为像距设为v；

步骤2.2：若己知物距v，设x₁、y₁为物体在三维空间的长和宽，x₂、y₂分别为物体摄像头所拍摄的二维图像的长和宽，则求得物距与焦距的关系公式，即光学二维空间小孔成像模型为：

焦距的大小决定着视场角的大小，焦距数值与视场角大小成反比。

步骤3：将二维空间坐标转化至三维空间坐标，将已知的若干控制点的坐标系进行数据转换。摄像头拍到的图像的是二维的平面，最终想要将二维平面图像变为3D立体图像并进行三维测量就需要在不同参考坐标系进行相关转换，具体需从图像坐标系经过其他坐标系转换到世界坐标系。

步骤3.1：建立像素坐标系：在图像上定义直角坐标系，每一像素的坐标点(a,b)表示其坐标值，a和b分别是这个像素在数组中的行数与列数；

步骤3.2：建立图像坐标系：X轴和Y轴分别平行于二维图像的两条垂直边，用(x,y)表示像素坐标值，设(a₀,b₀)为二维图像的中心，相机中感光器件的尺寸为x*y，像素坐标系与图像坐标系之间的关系表示为：

其中a-a₀和b-b₀分别表示某以像素点坐标数值与中心点坐标数值的X轴差值与Y轴差值；

将上式改写成矩阵形式得：

将偏移项带入进乘积项中，构造一个三维矩阵乘法形式，转化为齐次坐标形式为：

步骤3.3：建立相机坐标系：相机坐标系中心点设为O_x，图像坐标系中心点设为O_t，O_x与的O_t相交点的延长线为Z轴，即步骤3.2所述图像坐标系的延伸，Z轴的方向服从物理学中的右手规则，根据欧拉旋转定理知三维空间的旋转顺序不可逆性，X和Y轴在一平面内且平行，Z轴与X和Y轴所称平面不平行，三维空间中的点P在二维平面的投影点设为P'，再通过透视投影关系算出坐标，设点P在三维图像中的坐标为(x,y,z)，点P在二维平面投影点P'的坐标(x',y')用步骤2.3中物距与焦距的关系公式及相似三角形比例关系求出透视投影关系公式：

上式描述了点P和点P'之间的空间关系；

用齐次坐标与矩阵表示前述的透视投影关系公式为:

步骤3.4：建立世界坐标系；

所述世界坐标系是描述摄像头所在位置的一个基准坐标系，由X轴、Y轴、Z轴组成；摄像机坐标系与世界坐标系之间的关系在用旋转矩阵R与偏移向量T来描述，即物体只需要被旋转和平移；若在三维空间中存在一点P，其在世界坐标系下与像机坐标系下的齐次坐标分别是(x_s,y_s,z_s,1)与(x',y',z')，则存在如下关系：

建立像素坐标系与世界坐标系之间的变换关系：

式中，其中m及n是相机内部参数，与具体的摄影场景无关。

建立投影矩阵M为：

其中a₀、b₀为二维图像的中心点坐标，只与摄像机内部结构有关，为摄像机内部参数，是摄像头感光芯片安装时与理论中心的水平与垂直差值，

由摄像头相对于世界坐标系的方位决定，是相机相对于世界坐标系的平移变换关系，即摄像机外部参数；

步骤3.5：从图像坐标系转换到世界坐标系，己知摄像机的内外参数，获得投影矩阵M，对三维空间中的点P＝(x,y,z)，最终需要得到其在世界坐标系下的坐标P_s＝(x_s,y_s,z_s,1)，求出在图像坐标系点P的位置(x,y)，先从像素坐标系变换到相机坐标系，再由相机坐标系变换到世界坐标系，具体为：

其中，s为从世界坐标系转换到图像坐标系的尺度因子，t为平移矩阵。

步骤4：进行无人机相机标定，求出相机参数；

为确定图像中的空间物体表面某点的三维几何位置与其在图像中对应点之间的相互关系，需要建立摄像机成像的几何模型，利用步骤2和步骤3得到的数值和公式，其中的几何模型参数就是摄像机参数。求解参数的过程就称之为相机标定。这个步骤非常重要，会影响到最终测量结果的准确性，主要是求拍摄设备的内参、外参以及畸变参数。求从世界坐标系换到图像坐标系的过程，也就是求最终的投影矩阵P的过程。

步骤4.1：根据步骤3.4中建立的像素坐标系与世界坐标系之间的变换关系，使用张氏标定法，每一个角点的像素坐标设为(a,b)，每一个角点在世界坐标系下的物理坐标(A,B,C)＝0来进行相机的标定，获得相机的内参矩阵和径向畸变参数。

步骤4.2：将2D空间中的图像点设为n＝[a,b]^T，3D空间中的图像点设为N＝[x,y,z]，齐次坐标分别为

三维空间坐标映射到二维图像坐标转换为/>

其中s为世界坐标系到图像坐标系的尺度因子，α和β分别为焦距与像素横纵比，γ为径向畸变参数，/>

为相机内参矩阵；

步骤4.3：求解相机内参，三维空间坐标到二维图像坐标的映射为：H＝A[r1 r2t]，其中r1和r2为两个相互正交的旋转向量，A＝s*M，解Homographic矩阵进行单位性转换，将平面的点映射到另一个二维平面的投影，定义旋转矩阵的第i列为r_i，棋盘格位于Z轴上，表示为：s*[a b 1]^T＝A[r₁r₂ r₃ t]*c*[x y 0 1]^T。Homography有8个自由度，r₁和r₂正交且模相等，得到约束式为：

式中，h₁、h₂分别表示矩阵H的第一列与第二列；

构造一个矩阵B，使得

且约束条件为：

式中，v代表图像中棋盘格的角点纵坐标，因棋盘格是长方形，所以只需要得知这个棋盘格的其中三个角点坐标就可以确定剩余一个角点坐标。v11、v12、v22分别表示棋盘格的左上方、右上方、右下方角点的纵坐标，b表示用奇异值分解最小奇异值对应的特征向量。

内参表达式如下：

/>

其中，B₁₁、B₁₂、B₁₃、B₂₂、B₂₃、B₃₃公式表示矩阵B中不同位置元素，v₀、u₀表示中心点的实际位置坐标，λ为λ＝s^-1即尺度因子的逆矩阵，γ为径向畸变参数；

步骤4.4：求解相机外参：

r₁＝λM^-1h₁

r₂＝λM^-1h₁

r₃＝h₁*h₂

t＝λM^-1h₃

其中，r₁、r₂、r₃分别表示三个旋转向量，h₁、h₂、h₃分别表示矩阵H的第一列、第二列和第三列；

根据旋转矩阵的性质得：||r₁||＝||λM^-1h₁||＝1即

λ由步骤4.3求出，矩阵投影M可由步骤3.4求出，继而可求得h₁，用h₁求解推导出r₁、r₂、r₃和t。

步骤5：解决图像畸变问题，对图像进行矫正；无人机在拍摄目标图像的过程中受到透镜及相机自身出厂组装影响从而导致所拍摄的图像出现畸变，为了获得最佳的成像效果，需要在相机镜头前加透镜，但透镜会对光线的传播路线产生影响。透镜种类不同，成像的畸变程度也不同，为解决此影响应进行径向畸变矫正。为解决因在机械组装中的误差产生的影响，例如透镜和成像平面不平行等，应进行切向畸变矫正。两种畸变共同之处在于它们都会导致空间点成像画面不在实际线性模型所描述的位置，此时需要对图像进行矫正。

步骤5.1：为解决因相机摄像头前方的透镜对光的传播路径的影响，对图像进行径向畸变矫正。

径向畸变主要分为桶形畸变和枕形畸变，穿过图像中心和光轴由交点的延长直线能保持形状不变。对于图像中某一点为点P，设x_j为发生径向畸变后的该点横坐标，y_j为发生径向畸变后的该点纵坐标；x为校正后的该点横坐标，y为校正后的该点纵坐标，归一化极坐标表示为[r,θ]，其中r为矢量长度，t是矢量的水平夹角。可看成是某一个坐标点沿着长度方向发生了δr，其中δr＝r'-r即矢量长度的变化变化也就是它距离远点的长度发生了变化，都是随着与中心之间的距离增加而增加，则径向畸变的矫正公式为：

x_j＝x*(1+k₁r²+k₃r⁶+k₂r⁴)

y_j＝y*(1+k₁r²+k₃r⁶+k₂r⁴)

其中，k₁、k₂、k₃是径向畸变参数；

在本式中，纠正畸变比较小的图像中心区域时k₁起主要作用，纠正畸变比较大图像边缘区域时k₂起主要作用。

步骤5.2：为解决因在机械组装中的误差产生的影响，对图像进行切向畸变矫正。

切向畸变不存在对称的性质。可以卡成是坐标点沿着切线的方向发生了变化，水平夹角变化了δθ。对于图像中某一点为点P，设x₀为发生横向畸变后的该点横坐标，y₀为发生横向畸变后的该点纵坐标；x为校正后的该点横坐标，y为校正后的该点纵坐标；用另外两个参数p₁、p₂进行纠正，则切向畸变的矫正公式为：

x₀＝x+2p₁y+p₂(r²+2x²)

y₀＝y+2p₂x+p₁(r²+2y²)

步骤5.3：联合步骤5.1和步骤5.2中的径向畸变、切向畸变矫正公式，对于相机坐标系中的点P[x,y,z]，寻找它对应点的像素坐标的过程时将三维空间点投影到归一化图像平面，该点的像素坐标为：u＝f_x*x+c_x，v＝f_y*y+c_y。其中，u为像素点的横坐标，v为像素点的纵坐标，c_x为图像像素中心点的横坐标，c_y为图像像素中心点的横坐标。

步骤6：将图像进行预处理，图像质量的好坏影响到最后的测量数据。因图像是彩色图像，我们若直接对其三个通道R、G、B进行依次处理将会非常耗时，所以应对图像进行灰度化。

将图像进行预处理。也就是想增强图像的质量或突出图像某个特性，为了得到更加清晰的图像，需要通过某技术对图像进行处理，使用对比度增强的方法来处理图像，即对图像输出的灰度级进行放大。通过调整图像灰度值的分布来丰富所获取图像的灰度层次，目的是改善图像的检测效果。本实施例采用直方图均衡化的算法，将原图像的直方图修正为均匀的直方图，然后按均衡直方图修正原图像。将一副图像的直方图分布通过累积分布函数变成近似均匀分布，从而增强图像的对比度。均衡化图像后图像的直方图是平直的，灰度级具有均匀的概率分布图像就更清晰。如图2所示，本实施例的方法如下所述。

步骤6.1：直方图是在灰度范围[0,...,k,...,225]区间的功能函数，列出原图像的灰度级j，统计出原图像中每个灰度级所拥有的像素个数n_j；

步骤6.2：计算出原图像的各灰度级的灰度分布概率

其中n为直方图中的像素总个数，再计算出原图像的各灰度累积分布概率/>

其中k表示灰度级，，求出均衡化后的灰度级i＝INT[(L-1)S(i)+0.5]，INT()表示取整函数、L表示总灰度级数目、S(i)表示在i级的离散灰度，将((L-1)*p(k))进行四舍五入计算。

步骤6.3：确定原图像的灰度级与均衡化后的灰度级的映射关系j→i，再计算变换后的图像直方图

步骤7：设计一种端到端的单目三维物体测量模型，如图3所示，将预处理后的航拍图像输入到单目三维物体测量模型中，完成航空图片的三维测量。

步骤7.1：建立单目三维物体测量模型，所述单目三维物体测量模型包括backbone、关键点检测分支和3D框回归分支；其中所述backbone层使用DLA-34结构来提取特征，以便对不同层之间的特征进行融合。其中，DLA-34增加了FPN结构以此提高小物体检测的性能；将所获取的二维的RGB图像输入至backbone层中，通过DLA-34输出特征图，其长和宽分别为输入图像的1/4，通道数为256；在backbone中加入GN层，GN层将channel分组，在每个group中做LN，再做归一化处理；以此减小网络性能对batch_size的敏感程度，还可以减少训练时间。

所述二维RGB图像中，设W代表图像宽度，H代表图像高度；用(d,w,l,x,y,z,θ)表示目标三维图像，其中(d,w,l)的分别表示目标三维图像的高、宽、长，(x,y,z)表示目标三维图像在相机坐标系中的中心坐标点，θ表示拍摄图像时摄像头角度，拍摄图像时的飞行高度为s，相机内参设为K，

其中f_x及f_y在步骤4.3已经求出，c_x及c_y在步骤5.3已经求出。/>

步骤7.2：将步骤7.1中backbone层输出的特征图输入关键点检测分支网络层，所述关键点检测分支网络层用于寻找目标的3D中心点在二维图像中的投影点位置；在关键点检测分支网络层寻找输入图像中的关键点并标记绿色，将目标三维图像的中心点在二维图像平面的映射点[x_cy_c]^T设为关键点，目标三维图像的中心点坐标设为[x y z]^T，用坐标[x_b,1y_b,1]^T、[x_b,2y_b,2]^T、[x_b,3y_b,3]^T、[x_b,4y_b,4]^T、[x_b,5y_b,5]^T、[x_b,6y_b,6]^T、[x_b,7y_b,7]^T、[x_b,8y_b,8]^T分别表示目标三维图像的8个顶点在图像中的位置，并用

确定高斯函数的标准差，其中/>

表示为表示3D框的最小长方形。

步骤7.3：将步骤7.1中backbone层输出的特征图像及步骤7.2中所得得到的目标3D中心点输入至3D框回归分支网络层中，设δ_z表示深度偏移，δ_xc和δ_yc为量化误差，δ_d、δ_w、δ_l分别表示高、宽、长的残差，通过深度偏移求得目标的深度值z＝μ_z+δ_zσ_z，μ_z、σ_z分别尺度因子和平移因子，都为预先设定好的值，物体的目标三维图像的中心点坐标表示为：

再求目标物体的长、宽、高表达式为：

(单位为米)

在式中，

分别为公开的KITTI数据集中计算出的平均长度、宽度、高度，e^δl、e^δw、e^δd分别为长宽高的残差值。

步骤7.5：根据角度和目标在三维空间中的位置求得摄像头角度

α是预测的观测角，则拍摄该物体时，无人机距离目标物体的直线距离m为/>

通过将航向角矩阵R绕x轴旋转得到矩阵R_θ，根据物体的长宽高[l w d]^T和中心点坐标[x y z]^T得到目标三维物体的8个角点坐标矩阵表示为：

以上描述仅为本公开的较佳实施例以及对所运用技术原理的说明。本领域技术人员应当理解，本公开的实施例中所涉及的发明范围，并不限于上述技术特征的特定组合而成的技术方案，同时也应涵盖在不脱离上述发明构思的情况下，由上述技术特征或其等同特征进行任意组合而形成的其它技术方案。例如上述特征与本公开的实施例中公开的(但不限于)具有类似功能的技术特征进行互相替换而形成的技术方案。

Claims (7)

1.一种基于无人机单张航空图片的三维测量方法，其特征在于，包括以下步骤：

步骤1：利用无人机的相机进行单张图像拍摄并记录无人机拍摄时的飞行高度s；

步骤2：构造光学二维空间小孔成像模型，建立物距与焦距的关系，来确定空间中目标三维物体信息与图像二维成像之间的点面关系；

步骤3：将二维空间坐标转化至三维空间坐标，将已知的若干控制点的坐标系进行数据转换；

步骤4：进行无人机相机标定，求出相机参数；

步骤5：解决图像畸变问题，对图像进行矫正；

步骤6：将图像进行预处理；

步骤7：设计一种端到端的单目三维物体测量模型，将预处理后的航拍图像输入到单目三维物体测量模型中，完成航空图片的三维测量。

2.根据权利要求1所述的一种基于无人机单张航空图片的三维测量方法，其特征在于，所述步骤2具体包括以下步骤：

步骤2.1：建立物距与像距和焦距的关系表达式

其中焦距设为f，物体到镜头光心称为物距设为u，成像平面到镜头光心称为像距设为v；

步骤2.2：若己知物距v，设x₁、y₁为物体在三维空间的长和宽，x₂、y₂分别为物体摄像头所拍摄的二维图像的长和宽，则求得物距与焦距的关系公式，即光学二维空间小孔成像模型为：

3.根据权利要求1所述的一种基于无人机单张航空图片的三维测量方法，其特征在于，所述步骤3具体包括以下步骤：

步骤3.1：建立像素坐标系：在图像上定义直角坐标系，每一像素的坐标点(a,b)表示其坐标值，a和b分别是这个像素在数组中的行数与列数；

步骤3.2：建立图像坐标系：X轴和Y轴分别平行于二维图像的两条垂直边，用(x,y)表示像素坐标值，设(a₀,b₀)为二维图像的中心，相机中感光器件的尺寸为x*y，像素坐标系与图像坐标系之间的关系表示为：

其中a-a₀和b-b₀分别表示某以像素点坐标数值与中心点坐标数值的X轴差值与Y轴差值；

将上式改写成矩阵形式得：

将偏移项带入进乘积项中，构造一个三维矩阵乘法形式，转化为齐次坐标形式为：

步骤3.3：建立相机坐标系：相机坐标系中心点设为O_x，图像坐标系中心点设为O_t，O_x与的O_t相交点的延长线为Z轴，即步骤3.2所述图像坐标系的延伸，Z轴的方向服从物理学中的右手规则，根据欧拉旋转定理知三维空间的旋转顺序不可逆性，X和Y轴在一平面内且平行，Z轴与X和Y轴所称平面不平行，三维空间中的点P在二维平面的投影点设为P'，再通过透视投影关系算出坐标，设点P在三维图像中的坐标为(x,y,z)，点P在二维平面投影点P'的坐标(x',y')用步骤2.3中物距与焦距的关系公式及相似三角形比例关系求出透视投影关系公式：

用齐次坐标与矩阵表示前述的透视投影关系公式为:

步骤3.4：建立世界坐标系；

所述世界坐标系是描述摄像头所在位置的一个基准坐标系，由X轴、Y轴、Z轴组成；摄像机坐标系与世界坐标系之间的关系在用旋转矩阵R与偏移向量T来描述，即物体只需要被旋转和平移；若在三维空间中存在一点P，其在世界坐标系下与像机坐标系下的齐次坐标分别是(x_s,y_s,z_s,1)与(x',y',z')，则存在如下关系：

建立像素坐标系与世界坐标系之间的变换关系：

式中，其中m及n是相机内部参数；

建立投影矩阵M为：

其中a₀、b₀为二维图像的中心点坐标，只与摄像机内部结构有关，为摄像机内部参数，是摄像头感光芯片安装时与理论中心的水平与垂直差值，

由摄像头相对于世界坐标系的方位决定，是相机相对于世界坐标系的平移变换关系，即摄像机外部参数；

步骤3.5：从图像坐标系转换到世界坐标系，己知摄像机的内外参数，获得投影矩阵M，对三维空间中的点P＝(x,y,z)，最终需要得到其在世界坐标系下的坐标P_s＝(x_s,y_s,z_s,1)，求出在图像坐标系点P的位置(x,y)，先从像素坐标系变换到相机坐标系，再由相机坐标系变换到世界坐标系，具体为：

其中，s为从世界坐标系转换到图像坐标系的尺度因子，t为平移矩阵。

4.根据权利要求1所述的一种基于无人机单张航空图片的三维测量方法，其特征在于，

步骤4.1：根据步骤3.4中建立的像素坐标系与世界坐标系之间的变换关系，使用张氏标定法，每一个角点的像素坐标设为(a,b)，每一个角点在世界坐标系下的物理坐标(A,B,C)＝0来进行相机的标定，获得相机的内参矩阵和径向畸变参数；

步骤4.2：将2D空间中的图像点设为n＝[a,b]^T，3D空间中的图像点设为N＝[x,y,z]，齐次坐标分别为

三维空间坐标映射到二维图像坐标转换为

其中s为世界坐标系到图像坐标系的尺度因子，α和β分别为焦距与像素横纵比，γ为径向畸变参数，/>

为相机内参矩阵；

步骤4.3：求解相机内参，三维空间坐标到二维图像坐标的映射为：H＝A[r1 r2 t]，其中r1和r2为两个相互正交的旋转向量，A＝s*M，解Homographic矩阵进行单位性转换，将平面的点映射到另一个二维平面的投影，定义旋转矩阵的第i列为r_i，棋盘格位于Z轴上，表示为：s*[a b 1]^T＝A[r₁r₂ r₃ t]*c*[x y 0 1]^T；Homography有8个自由度，r₁和r₂正交且模相等，得到约束式为：

式中，h₁、h₂分别表示矩阵H的第一列与第二列；

构造一个矩阵B，使得

且约束条件为：

式中，v代表图像中棋盘格的角点纵坐标，因棋盘格是长方形，所以只需要得知这个棋盘格的其中三个角点坐标就可以确定剩余一个角点坐标；v11、v12、v22分别表示棋盘格的左上方、右上方、右下方角点的纵坐标，b表示用奇异值分解最小奇异值对应的特征向量；

内参表达式如下：

其中，B₁₁、B₁₂、B₁₃、B₂₂、B₂₃、B₃₃公式表示矩阵B中不同位置元素，v₀、u₀表示中心点的实际位置坐标，λ为λ＝s^-1即尺度因子的逆矩阵，γ为径向畸变参数；

步骤4.4：求解相机外参：

r₁＝λM^-1h₁

r₂＝λM^-1h₁

r₃＝h₁*h₂

t＝λM^-1h₃

其中，r₁、r₂、r₃分别表示三个旋转向量，h₁、h₂、h₃分别表示矩阵H的第一列、第二列和第三列；

根据旋转矩阵的性质得：||r₁||＝||λM^-1h₁||＝1即

用h₁求解推导出r₁、r₂、r₃和t。

5.根据权利要求1所述的一种基于无人机单张航空图片的三维测量方法，其特征在于，所述步骤5具体包括以下步骤：

步骤5.1：为解决因相机摄像头前方的透镜对光的传播路径的影响，对图像进行径向畸变矫正；

对于图像中某一点为点P，设x_j为发生径向畸变后的该点横坐标，y_j为发生径向畸变后的该点纵坐标；x为校正后的该点横坐标，y为校正后的该点纵坐标，归一化极坐标表示为[r,θ]，其中r为矢量长度，t是矢量的水平夹角，则径向畸变的矫正公式为：

x_j＝x*(1+k₁r²+k₃r⁶+k₂r⁴)

y_j＝y*(1+k₁r²+k₃r⁶+k₂r⁴)

其中，k₁、k₂、k₃是径向畸变参数；

步骤5.2：为解决因在机械组装中的误差产生的影响，对图像进行切向畸变矫正；

对于图像中某一点为点P，设x₀为发生横向畸变后的该点横坐标，y₀为发生横向畸变后的该点纵坐标；x为校正后的该点横坐标，y为校正后的该点纵坐标；用另外两个参数p₁、p₂进行纠正，则切向畸变的矫正公式为：

x₀＝x+2p₁y+p₂(r²+2x²)

y₀＝y+2p₂x+p₁(r²+2y²)

步骤5.3：联合步骤5.1和步骤5.2中的径向畸变、切向畸变矫正公式，对于相机坐标系中的点P[x,y,z]，寻找它对应点的像素坐标的过程时将三维空间点投影到归一化图像平面，该点的像素坐标为：u＝f_x*x+c_x，v＝f_y*y+c_y；其中，u为像素点的横坐标，v为像素点的纵坐标，c_x为图像像素中心点的横坐标，c_y为图像像素中心点的横坐标。

6.根据权利要求1所述的一种基于无人机单张航空图片的三维测量方法，其特征在于，所述步骤6具体包括以下步骤：

步骤6.1：直方图是在灰度范围[0,...,k,...,225]区间的功能函数，列出原图像的灰度级j，统计出原图像中每个灰度级所拥有的像素个数n_j；

步骤6.2：计算出原图像的各灰度级的灰度分布概率

其中n为直方图中的像素总个数，再计算出原图像的各灰度累积分布概率/>

其中k表示灰度级，求出均衡化后的灰度级i＝INT[(L-1)S(i)+0.5]，其中INT()表示取整函数、L表示总灰度级数目、S(i)表示在i级的离散灰度，将(L-1)*p(k)进行四舍五入计算；

步骤6.3：确定原图像的灰度级与均衡化后的灰度级的映射关系j→i，再计算变换后的图像直方图

7.根据权利要求1所述的一种基于无人机单张航空图片的三维测量方法，其特征在于，所述步骤7具体包括以下步骤：

步骤7.1：建立单目三维物体测量模型，所述单目三维物体测量模型包括backbone、关键点检测分支和3D框回归分支；其中所述backbone层使用DLA-34结构来提取特征，将所获取的二维的RGB图像输入至backbone层中，通过DLA-34输出特征图，其长和宽分别为输入图像的1/4，通道数为256；在backbone中加入GN层，GN层将channel分组，在每个group中做LN，再做归一化处理；

所述二维RGB图像中，设W代表图像宽度，H代表图像高度；用(d,w,l,x,y,z,θ)表示目标三维图像，其中(d,w,l)的分别表示目标三维图像的高、宽、长，(x,y,z)表示目标三维图像在相机坐标系中的中心坐标点，θ表示拍摄图像时摄像头角度，拍摄图像时的飞行高度为s，相机内参设为K，

步骤7.2：将步骤7.1中backbone层输出的特征图输入关键点检测分支网络层，所述关键点检测分支网络层用于寻找目标的3D中心点在二维图像中的投影点位置；在关键点检测分支网络层寻找输入图像中的关键点并标记绿色，将目标三维图像的中心点在二维图像平面的映射点[x_c y_c]^T设为关键点，目标三维图像的中心点坐标设为[x y z]^T，用坐标[x_{b, 1}y_b,1]^T、[x_b,2y_b,2]^T、[x_b,3y_b,3]^T、[x_b,4y_b,4]^T、[x_b,5y_b,5]^T、[x_b,6y_b,6]^T、[x_b,7y_b,7]^T、[x_b,8y_b,8]^T分别表示目标三维图像的8个顶点在图像中的位置，并用

确定高斯函数的标准差，其中/>

表示为表示3D框的最小长方形；

步骤7.3：将步骤7.1中backbone层输出的特征图像及步骤7.2中所得得到的目标3D中心点输入至3D框回归分支网络层中，设δ_z表示深度偏移，δ_xc和δ_yc为量化误差，δ_d、δ_w、δ_l分别表示高、宽、长的残差，通过深度偏移求得目标的深度值z＝μ_z+δ_zσ_z，μ_z、σ_z分别尺度因子和平移因子，物体的目标三维图像的中心点坐标表示为：

再求目标物体的长、宽、高表达式为：

在式中，

分别为公开的KITTI数据集中计算出的平均长度、宽度、高度，e^δl、e^δw、e^δd分别为长宽高的残差值；

步骤7.5：根据角度和目标在三维空间中的位置求得摄像头角度

α是预测的观测角，则拍摄该物体时，无人机距离目标物体的直线距离m为/>

通过将航向角矩阵R绕x轴旋转得到矩阵R_θ，根据物体的长宽高[l w d]^T和中心点坐标[x y z]^T得到目标三维物体的8个角点坐标矩阵表示为：

Images