赤峰打印机价格协会

论文推荐|赵爽:基于大比例尺航空影像共面约束条件的相机自检校方法

测绘学报2022-01-11 13:09:19

《测绘学报》

构建与学术的桥梁        拉近与权威的距离


基于大比例尺航空影像共面约束条件的相机自检校方法

赵爽1,2, 李学军1, 刘涛1, 谢剑薇1     

1. , 北京 101416; 
2. 北京跟踪与通信技术研究所, 北京 100094

收稿日期:2016-05-09; 修回日期:2015-11-15

基金项目:国防预研项目(513150701)

第一作者简介: 赵爽(1979-), 女, 博士生, 研究方向为计算机图形学、遥感影像处理。

E-mail:zhsh2002@163.com

摘要:提出了基于大比例尺航空影像共面约束的相机自检校方法,该方法使用所有立体像对同名点基于共面约束对相机的内方位元素及畸变系数进行解算。首先进行航空影像同名点匹配,构建立体像对;然后基于共面约束使用直接解法和迭代优化进行相对定向,解算相机位置与姿态;最后使用最小二乘优化方法解算相机内方位元素和畸变系数。对于高分辨率大尺寸航空影像,图像中心及边缘的畸变差异较大,为了进一步提高解算精度,对图像进行网格区域划分解算畸变。使用大比例尺航空影像进行解算能真实精确反映航空摄影测量时所获取图像的相机参数和畸变系数,避免检校环境与使用环境不同解算得到的相机畸变参数不能真实反映所获取影像的畸变问题;使用所有同名点解算,避免由于选择不同特征点或控制点对检校精度的影响;通过区域网格划分,进一步提高了解算精度。对检校结果进行了分析,该方法精度较高,与基于室外检校场的精度相当,能真实精确反映航空摄影测量时所获取图像的相机参数和畸变系数,提高了三维重建的精度。

关键词: 大比例尺航空影像     共面约束     相机自检校     网格区域划分    

Camera Self-calibration Method Based on Large Scale Aerial Images Coplanarity Constraint

ZHAO Shuang1,2, LI Xuejun1, LIU Tao1, XIE Jianwei1     

Abstract: In this paper an automated camera self-calibration method based on large scale aerial images coplanarity constraint is proposed. The method uses stereo images corresponding points to compute camera inner parameters and the distortion coefficients based on coplanarity constraint. Firstly, aerial images corresponding points are matched and stereo images are constructed. Then the relative orientation of the camera is resolved by using the direct method and the iterative optimization based on the coplanarity constraint. Finally, the least square optimization method is used to calculate the inner parameters and the distortion coefficients of the camera. For large image the distortion of image boundary is more severe than center, in order to compute parameters more precisely mesh partition in arbitrary area of image is proposed. Compared with classical techniques which use expensive equipment or calibration field, the proposed technique can reflect the actual inner parameters of the camera in the real environment. Using all the corresponding points avoid the influence of choosing different feature points or control points on the calibration accuracy.Experimental results indicate that it is effective and flexible.The accuracy of the method is equivalent to that of calibration field, which can accurately reflect the camera parameters and the distortion coefficients of the image acquired by aerial photogrammetry.

Key wordslarge scale aerial image     coplanarity constraint     camera self-calibration     mesh partition    

相机检校是摄影测量及计算机视觉至关重要的一个问题,相机检校是否精确,会影响到后续测图及三维重建的精度。影响相机检校结果的因素很多,包括相机自身因素和外部因素等,如相机制造工艺、相机焦距、调节光圈、更换镜头、相机晃动、更换电池、光照条件、检校靶标数量、靶标点提取精度、安装角、测量误差、相机模型选择等[1-4],因此在进行摄影测量时,在获得相机的初始参数后,还需进行相机检校,减少相机参数对测量精度的影响[5-6]。摄影测量领域的传统方法通常是在专用检校场(鉴定场)进行相机检校[7-8],但检校场建造维护代价很高,且要综合考虑气候条件、地形地貌等的影响。为了满足相机鉴定的迫切需求,使相机鉴定更加方便灵活,文献[9-12]使用靶标进行相机鉴定。基于靶标的相机鉴定结果精度很高,稳健性强,但需要制作特定靶标,而且需要精密仪器测量世界坐标。

最好、最适宜的相机检校方法取决于实际应用类型[13],立体相机根据两个相机之间固定相对关系建立约束进行标定[14],存在“空间平行线”等建筑物场景时,则可利用灭点进行相机检校[15]。为了真实精确反映航空摄影测量时所获取图像的相机参数和畸变系数,提高解算精度,进行精确三维建模。本文提出基于航空影像的相机自检校方法,联合所有航空影像同名点基于共面约束进行相机检校。对于高分辨率大尺寸航空影像,中心及边缘的畸变差异较大,为了提高畸变参数解算精度,提出分区域网格划分检校法。

1 相机检校

航空影像是现代信息化战场的重要信息来源,为了提高数据处理精度,保证数据的可信性,关键一步就是对相机进行严格的检校,求出相机的内方位元素及各项畸变系数。在数据处理中消除影像的畸变差,提高同名像点的空间解算精度,实现高精度三维重建。本文提出基于航空影像的相机自检校方法,具体步骤如下:首先进行航空影像同名点匹配,构建立体像对;然后基于共面约束使用直接解法和迭代优化进行相对定向,解算相机位置与姿态;最后对航空影像进行区域划分,使用最小二乘优化方法解算相机内方位元素和畸变系数。由于该方法联合所有影像同名点进行解算,为了避免大量同名点对存储和计算速度的影响,优化了矩阵解算策略,使方程解算与同名点数量无关。该方法操作简单,对图像无特殊要求,不需要制作精良靶标,也不需要使用特定检校场。使用大比例尺航空影像进行解算能真实精确反映摄影测量时所获取图像的相机参数和畸变系数,避免检校环境与使用环境不同解算得到的相机畸变参数不能真实反映所获取影像的畸变问题;使用所有同名点解算,避免由于选择不同特征点或控制点对检校精度的影响[16],同时,焦距、像主点和相机畸变系数可以独立解算,降低了参数之间的耦合度[12]

在进行相机检校前首先必须根据实际应用选择合理的相机畸变模型,模型的选择要比参数估计更重要,同时还需要选择合适的评价指标衡量选择的模型[17]。本文研究了小孔模型、经典的畸变模型及在xy方向使用不同的径向畸变系数的畸变模型,同时提出区域网格划分法进行畸变解算,并对检校结果进行了比较。同名点误差及共面约束误差对三维重建精度影响很大,因此本文选择同名点误差及共面误差衡量各模型的性能。

1.1 相机模型

被广泛使用的相机模型为小孔成像模型,如式所示

(1)

式中,s为空间点深度或尺度因子;(uv)为以像素为单位的图像坐标;(XwYwZw)为世界坐标系中的坐标;(u0v0)为像主点坐标;αu轴和v轴不垂直因子;(fufv)分别是相机以像素为单位的横坐标轴和纵坐标轴上的焦距;R为旋转矩阵;t为平移向量。

本文设u轴和v轴具有相同的尺度因子,即fu=fv,并假设成像平面的横、纵轴垂直,则最终的相机投影变换模型为

(2)

式中,f为等效焦距;K为内方位元素矩阵或投影变换矩阵;M为外参数矩阵,包含旋转矩阵R和平移向量t

由于制造工艺等多种因素,相机通常都会存在畸变,如式(3)所示

(3)

式中,(xy)为中心化的像点坐标,x=x0-cxy=y0-cy;(x0y0)为实际像点坐标;cxcy为光心坐标;r为像点到光心的距离,k1k2为径向畸变系数;p1p2为偏心畸变系数;s1s2为薄棱镜畸变系数。

对于高分辨率大尺寸航空影像,xy方向的径向畸变系数可能并不相同,因此在进行解算时,设xy方向具有不同的径向畸变系数kx1kx2ky1ky2,如式(4)所示

(4)

由于径向畸变的影响远大于其他畸变[18-19],因此本文只考虑径向畸变。

1.2 相机检校


1.2.1 同名点匹配

同名点匹配是进行相机检校的基础。输入图像为无序的遥感影像集,因此必须首先提取特征点,进行同名点匹配,进而确定各影像间的邻接关系,建立立体模型。从无序影像集构建立体模型的策略如下:首先提取所有影像的特征点;对于任意两幅影像,进行特征匹配,因为从不同视点拍摄的两幅影像存在缩放、旋转等仿射变换,基于RANSAC算法求解它们之间的仿射变换关系,以判断匹配的有效性:若两幅图像无重叠或误点太多,则无法得到一个有效的仿射变换,该像对被剔除;对于有重叠并匹配成功的两幅影像,通过RANSAC算法可以剔除少量误匹配点,最终获得一定数量的可靠同名像点,对应的两幅影像即构成了一个立体像对[20]。在立体像对模型内进行稠密同名点匹配,获取大量的观测值,可以为后续的最小二乘优化提供大量的观测方程。


1.2.2 相对定向

根据解析法相对定向,同名射线在空间对对相交,射线opo′p′及摄影基线B共面,如图 1所示,根据三矢量共面,它们的混合积为0,可得约束公式

图 1 共面约束示意图Fig. 1 Coplanarity constraint

图选项


(5)

式中,opo′p′为同名射线矢量;B为基线矢量。

首先采用直接解法相对定向获取各相机位置与姿态初值,然后使用最小二乘进行迭代优化,依据共面约束,建立优化函数如式(6)所示

(6)

式中,F=B·op×o′p′

具体过程如下:①选择系列图像中中间位置的图像作为主图,相邻的一幅图作为辅图,邻接图为相互有重叠的图,从其他未定位图中找一幅邻接点最多的邻接图作为候选图。通过提取主图与邻接图间的同名点,由相对定向过程,即可确定邻接图对应相机的相对位置和姿态;重复上述过程,即可确定全部相机的相对位置和姿态。②以位置与姿态参数为优化目标,基于共面约束建立优化函数,解算各相机位置与姿态参数。对于多幅航空影像,为了快速解算位置与姿态,迭代解算能够收敛,避免参数过多相互干扰,采用如下策略:对主、辅图定位后,主、辅图只绕Z轴旋转κ1κ2,不进行平移变换,其他图像进行平移调整(txityitzi),同时绕Z轴旋转κi,经过迭代调整后使同名像点重合。此次位置与姿态解算不加入相机内方位元素。主、辅图各对应1个位置与姿态参数,其他图像对应4个位置与姿态参数,对于M幅图像N对同名点,则可生成包含(4M-6)个未知数N个方程的线性约束方程组,如式(7)所示

(7)

式中

经过直接解算与迭代优化后可以解算出相机相对位置与姿态。


1.2.3 相机内方位元素解算

在相对定向解算的位置与姿态基础上加入相机内方位元素(焦距、像主点、畸变系数)进行联合解算,进一步优化位置与姿态参数,解算相机内方位元素。约束方程和优化目标不变,以相机相对于公共坐标系的旋转角参数、平移向量及相机内方位元素为未知量,作为最小二乘参数迭代优化,如式(8)所示

(8)

式中

在调整过程中,主图绕XY轴旋转ω1φ1,辅图绕3个轴旋转ω2φ2κ2,其他图像进行平移调整(txityitzi),同时绕3个轴旋转ωiφiκi。主、辅图对应5个位置与姿态参数,其他图像对应6个位置与姿态参数,共有7个相机参数(fcxcykx1kx2ky1ky2),对于M幅图像N对同名点,则可生成包含6M个未知数、N个方程的线性约束方程组。

经过最小二乘迭代优化,可以解算出相机的内方位元素和畸变系数及更精确的位置与姿态参数。


1.2.4 分区域畸变校正

对于高分辨率大尺寸航空影像,图像中心及边缘的畸变差异较大,为了准确地解算相机畸变系数,进行畸变校正,对图像进行了区域网格划分,分区域进行畸变校正。将航空影像划分为m×n矩形网格,mn分别为xy方向的网格数。联合所有影像同名点使用双线性插值和最小二乘迭代优化算法解算网格交点畸变量,然后根据网格交点处畸变量校正该网格区域内像点。通过联合计算所有航空影像并通过合理划分矩形网格,能够最真实地反映摄影时相机的内方位元素及畸变系数,对像点进行畸变校正。

图 2中,设网格(ij)交点pij的畸变量为δxpijδypij,则根据双线性插值,同名像点Qk处的畸变量为

图 2 区域网格划分Fig. 2 Mesh partition

图选项


(9)

式中,uv分别为Qk到直线Pi+1, jPi+1, j+1Pij+1Pi+1, j+1的归一化距离。然后进行最小二乘优化迭代,在位置与姿态参数、相机内方位元素的基础上加入所有网格交点处畸变量进行迭代优化。xy方向的网格数分别为mn,则网格交点畸变量参数个数为2×m×n,对于N对同名点,构成的线性约束方程组如式(10)所示

(10)

式中


2.2.5 优化策略

联合所有航空影像大规模同名点进行解算,对存储和计算要求较高。本文采用最小二乘方法进行迭代优化,方程AX=BATAX=ATB同解,为了节省存储空间,不保存系数矩阵A而直接构造ATA,方程右边直接构造ATB,对于式(8)最终所需的存储单元仅为6M×(6M+1)。常规方法需要的存储单元为2N×(6M+1),这样的处理使此方法的方程解算规模与引入的同名点数量无关,从而极大地减少了存储量,提高了计算速度。多次迭代直至||δX2||<ε(ε为某一小值),进而得到未知量的最优解。

2 试验2.1 试验数据

试验数据1采用威海天福山试验场拍摄的航空影像[21],试验场面积为2 km2,场区内均匀布设了183个精确量测的地面控制点,每个控制点都制作了明显的地面标志,控制点坐标精度优于0.01 m。用有人驾驶飞机对场区进行航拍,得到了5条航线、共115幅航空影像,原始影像数据量约为20 GB;图像分辨率为8956×6708像素。图像的航向重叠率约为60%,旁向重叠率大于20%,影像空间拍摄相机为哈苏H4D-60,镜头焦距为50 mm,像元大小为6 μm。无POS记录。图 3(a)为试验用的部分影像。

图 3 航空影像数据Fig. 3 Aerial images

图选项


试验数据2采用使用佳能5DⅡ相机对青岛王家村进行航空摄影拍摄的航空影像,全航区共拍摄了78(13×6航线)幅影像,镜头焦距为35 mm,像元大小为6.4 μm,图像大小为5616×3744像素。部分试验数据如图 3(b)所示。

2.2 试验结果

从数据1威海天福山试验场的22幅航空影像64个同名点文件中匹配14幅图像32个同名点文件进行了相机内方位元素及畸变联合解算,同名点对数量约为640 000对。参数初始值(区域网格划分时网格交点畸变量初值为0)及试验结果如表 1所示。畸变量(δxδy)为位于图片角上的像点在xy方向的径向畸变,用于评估畸变大小。共面误差(MSE)为左像同名点和基线构成平面与右像同名点的距离的均方误差。方法1和方法2为文献[18]基于室外检校场联立全部影像和4幅影像的检校结果;方法*-1(*=1, 2, 3)为本文方法单独解算焦距(f)和像主点(cxcy);方法*-2(*=1, 2, 3)为使用传统畸变模型解算相机内方位元素和畸变系数;方法*-3为xy方向使用不同的径向畸变参数解算相机内方位元素和畸变系数;方法*-4、*-5和*-6为采用区域网格划分进行畸变解算,网格数分别为16×12,32×24,64×48。

表 1 文献[18]和本文方法检校结果Tab. 1 Intrinsic parameters of camera obtained from refs. [18] and our method

方法焦距/像素像主点/像素畸变系数畸变量/像素共面误差
/像素

fcxcykx1(k1)kx2(k2)ky1ky2δxδyMSE
方法18 303.04 433.33 415.53.21e-11-1.94e-18--4.013.00-
方法28 313.24 431.93 417.23.96e-11-1.43e-18--0.720.54-
初始值18 308.34 477.53 353.50000---
方法1-18 332.24 433.83 413.3------0.441
方法1-28 336.94 433.73 413.63.03e-116.40e-19--7.055.280.435
方法1-38 343.94 433.03 413.73.83e-115.30e-192.80e-119.72e-197.696.130.435
方法1-4-------8.156.260.343
方法1-5-------7.896.170.341
方法1-6-------7.766.150.338
初始值28 303.04 433.33 415.50000---
方法2-18 331.84 433.43 413.5------0.438
方法2-28 337.44 433.73 413.53.05e-116.43e-19--7.095.310.435
方法2-38 344.64 433.13 413.73.87e-115.32e-192.82e-119.76e-197.766.160.435
方法2-4-------8.076.280.341
方法2-5-------7.796.190.338
方法2-6-------7.706.130.336
初始值38 333.344783 354.00000---
方法3-18 332.24 433.73 413.2------0.441
方法3-28 336.84 433.73 413.53.02e-116.40e-19--7.045.270.435
方法3-38 343.94 433.03 413.73.83e-115.31e-192.79e-119.73e-197.696.120.434
方法3-4-------8.186.300.368
方法3-5-------7.886.150.342
方法3-6-------7.856.170.337

表选项


表 1试验结果可以看出,本文方法可以有效求出相机内方位元素,且精度较高。像主点解算结果非常稳定,与基于室外检校场的解算精度相当。基于区域网格划分的畸变校正方法共面误差明显小于其他畸变校正方法。相机焦距的检校结果与初始值及文献[18]相差较大,经分析主要原因如下:一是由于相机镜头加工、装配过程中引入的误差,使实际值与标称值存在误差;二是在不同的环境下拍摄,摄影测量时外部条件变化等引入的误差[122-27];如不同的摄影距离下主距和像主点也会发生变化,不同的检校场、图像集检校结果也会有差异[27-28]

表 2为本文方法焦距和结果分析,可以看出该方法对优化算法初值不敏感,解算的结果非常稳定。表 3图 4为各图定向后同名点误差,使用区域网格划分的畸变解算算法精度明显高于其他方法,同名点误差显著减少,解算结果对初值不敏感,所有最小二乘迭代优化均小于7次收敛。

表 2 焦距、像主点检校结果分析Tab. 2 Analysis of the result of focal length and principal point

像素
参数均值标准差最小值中值最大值
cx4 433.40.344 433.04 433.74 433.8
cy3 413.70.683 413.23 413.63 415.5
f8 337.75.268 331.88 336.98 344.6

表选项


表 3 畸变校正后同名点误差Tab. 3 The error of corresponding points after calibration

像素
方法图像编号
90629063906490659088908990909091909291129113911491159116
方法1-10.3840.3880.4130.4650.5600.3890.3890.4010.3650.5680.4050.3610.4070.570
方法1-20.4030.3980.4190.4880.5760.4070.4160.4200.3770.5780.4170.3690.4090.571
方法1-30.4000.3960.4180.4870.5780.4070.4160.4190.3770.5780.4180.3720.4100.573
方法1-40.3220.3250.3540.4130.4980.3120.3080.3290.3160.5480.3620.2710.2490.444
方法1-50.3190.3240.3530.4110.4960.3100.3060.3270.3130.5430.3570.2670.2440.441
方法1-60.3160.3210.3500.4090.4940.3080.3040.3240.3100.5400.3540.2650.2430.439
方法2-10.4060.3990.4170.4810.5810.4140.4220.4230.3910.5850.4220.3830.4240.578
方法2-20.4030.3980.4190.4880.5760.4070.4160.4200.3770.5780.4170.3690.4090.571
方法2-30.4000.3960.4180.4870.5780.4070.4160.4190.3770.5780.4180.3720.4100.573
方法2-40.3190.3230.3530.4110.4950.3100.3060.3260.3140.5460.3560.2660.2440.439
方法2-50.3160.3220.3520.4100.4940.3090.3040.3230.3100.5410.3530.2630.2410.437
方法2-60.3130.3200.3500.4080.4920.3070.3020.3200.3070.5380.3500.2610.2390.435
方法3-10.4060.3990.4170.4810.5810.4140.4220.4230.3910.5850.4220.3830.4240.578
方法3-20.4030.3980.4190.4880.5760.4070.4160.4200.3770.5780.4170.3690.4090.571
方法3-30.4000.3960.4180.4870.5780.4070.4160.4190.3770.5780.4180.3720.4100.573
方法3-40.3440.3470.3720.4280.5090.3340.3320.3550.3320.5520.3930.3080.2910.478
方法3-50.3210.3250.3540.4130.4960.3120.3080.3290.3150.5440.3580.2670.2450.442
方法3-60.3140.3200.3500.4080.4920.3070.3020.3210.3080.5390.3510.2620.2400.436

表选项


图 4 同名点误差柱状图Fig. 4 The error of corresponding points after calibration

图选项


表 4-表 6图 5为畸变校正前后共面误差,采用本文方法进行畸变校正后共面误差显著减小,使用分区域畸变校正后,共面误差在0.4像素以内。

表 4 方法1畸变校正前后共面误差比较Tab. 4 Error of coplanarity constraint for first method

像素
方法方法1-1方法1-2方法1-3方法1-4方法1-5方法1-6
畸变校正前8.2568.2568.2568.2568.2568.256
畸变校正后0.4410.4350.4350.3430.3410.338

表选项


表 5 方法2畸变校正前后共面误差比较Tab. 5 Error of coplanarity constraint for second method

像素
方法方法2-1方法2-2方法2-3方法2-4方法2-5方法2-6
畸变校正前8.2358.2358.2358.2358.2358.235
畸变校正后0.4380.4350.4350.3410.3380.336

表选项


表 6 方法3畸变校正前后共面误差比较Tab. 6 Error of coplanarity constraint for third method

像素
方法方法3-1方法3-2方法3-3方法3-4方法3-5方法3-6
畸变校正前8.2478.2478.2478.2478.2478.247
畸变校正后0.4410.4350.4340.3680.3420.337

表选项


图 5 共面误差比较Fig. 5 Error of coplanarity constraint

图选项


为了进一步验证本文方法对普通数码相机的标定精度,使用数据2中10幅图像进行了试验,试验结果如表 7所示。从试验结果可以看出,该方法对普通数码相机标定精度也很高。

表 7 佳能5DⅡ拍摄影像试验结果Tab. 7 Intrinsic Parameters of Camera Canon 5DⅡObtained by Our Method

方法焦距/像素
像主点/像素
畸变系数
畸变量/像素
共面误差
/像素
fcxcykx1(k1)kx2(k2)ky1ky2δxδyMSE
初始值15 468.0
2 807.51 871.5









方法1-15 457.22 768.31 860.7






方法1-25 460.62 768.41 860.52.48E-09-1.66E-16--18.112.01.012
方法1-35 462.82 768.61 860.82.52E-09-1.75E-162.21E-09-1.41E-1615.912.40.891
方法1-4






16.312.80.612
方法1-5






16.812.70.548

表选项


3 结论

本文提出了基于大比例尺航空影像利用共面约束的相机自检校方法。与传统方法相比,该方法操作简单,精度较高,能真实反映航空摄影测量时相机参数。通过对高分辨率大尺寸航空影像进行区域划分解算相机畸变量,进一步提高了解算精度。通过实际图像运算结果可以看出,畸变校正后同名点匹配精度有了显著提高,共面误差在0.4像素以内,该方法已应用于大规模航空影像拼接,进行三维重建。

【引文格式】赵爽,李学军,刘涛,等。 基于大比例尺航空影像共面约束条件的相机自检校方法[J]. 测绘学报,2017,46(1):98-106. DOI: 10.11947/j.AGCS.2017.20160206

更多精彩内容:

审稿专家必看,如何写出一篇优质的评审报告!


两会 | 政协委员王名:建议使用遥感技术监测机动车尾气


行业动态 | 欧洲环境监测卫星“哨兵-2B”成功入轨


论文推荐|唐炉亮:网络空间中线要素的核密度估计方法


论文推荐|龙江平:联合干涉相位和相干性幅度的极化干涉SAR最优相干性估计


招聘|2017年中国地图出版集团应届毕业生招聘公告


招聘 | 黑龙江工程学院招聘40名教师


论文推荐|朱庆:顾及纹理特征的航空影像自适应密集匹配方法


活动 |“司南导航杯”书香测绘 有奖征文活动征稿启事


政务 | 关于推荐2017年全国优秀测绘工程奖的通知


两会上,测绘与导航领域院士专家的声音!


论文推荐|于英:影像连接点均衡化高精度自动提取


行业动态|北斗拓展民用市场 产业化应用全面提速


招聘 | 南京大学国际地球系统科学研究所博士后招聘


论文推荐|王密:高分四号静止轨道卫星高精度在轨几何定标


招生 | 南师大中美合作举办地理信息科学硕士学位教育项目招生简章


权威 | 专业 | 学术 | 前沿

微信投稿邮箱 | song_qi_fan@163.com



微信公众号中搜索「测绘学报」,关注我们,长按上图二维码,关注学术前沿动态。

欢迎加入《测绘学报》作者QQ群: 297834524


进群请备注:姓名+单位+稿件编号