单目相机 svd 从图像恢复3维位置_IMU辅助下的单目视觉坐标传递

22ee3fcdb63d900f733111a660b0e059.png

随着社会的进步和发展,基础设施建设日趋完善,各种管线管道的建设也日趋复杂。为了更方便地对管线设施进行维护,在其故障进行检测时及时确定故障位置,并建立数据库的需求迫在眉睫。数据库建立过程中,有效地测量显得尤其重要[1]

随着测绘技术的不断发展,可以应用于管线测量的方法很多,如全站仪、GNSS单点定位技术、RTK技术等[2-7]。全站仪测量可提供高精度坐标,但作业效率低;GNSS伪距单点定位精度在米级,不能满足精度要求;而GNSS精密单点定位技术需要较长时间固定模糊度,也不是一种高效的测量手段。

实时动态(real-time kinematic,RTK)载波相位差分技术是一种基于载波差分的GNSS技术[8],其实时性好,作业效率高,定位精度也比较高,在管线测量任务中被广泛应用。但是在GNSS信号条件不佳或地形不满足设站要求的情况下,RTK无法提供可靠的待测点坐标。这种情况下,可以利用激光测距仪测距并进行距离交会得到点位坐标,然而激光测距精度无法保证,并且单次测量仅能获取一个待测点坐标。因此,本文提出一种惯性测量单元(inertial measurement unit,IMU)辅助下的单目视觉坐标传递方法,在满足测量精度要求的同时,能够有效提高管线测量任务的作业效率。

1 算法框图及理论基础 

1.1 算法框图

算法流程如图 1所示,其具体过程为:①当相机内参和IMU相机间的外参未知时,需要对其进行标定,测得相机的内参矩阵和畸变参数及相机和IMU间的外参;②在两个已知点上分别架设相机,并朝向待测点方向拍照,尽量保证两帧图像有足够大的共视区域;③在两帧图像中提取特征并匹配;④使用对极几何的方法,利用匹配点对恢复两帧图像间的相对旋转与平移,但由于尺度未知,需要利用已知坐标点间的距离恢复真实尺度;⑤利用两帧图像间的旋转与平移三角化待测点,得到待测点在相机坐标系下的坐标;⑥利用IMU测得的重力矢量及相机光心间的矢量,可以计算出相机坐标系和真实地理坐标系间的相对旋转;⑦最后将待测点在相机坐标系下的坐标转到真实地理坐标系下。

4543b6232919d394ef77b984606bd7da.png
图 1 算法流程

1.2 对极几何恢复相机间相对旋转和平移

本文算法中,帧间旋转与平移的精度直接决定了系统精度。由于视场中没有控制点,只能利用两帧图像的匹配像点来恢复相对旋转和平移,此时只能通过2D-2D的方式,即利用对极几何方法[9-10]计算帧间旋转与平移。

对极几何约束如图 2所示,两帧图像I1I2对应的相机光心为O1O2,设空间点P,理想情况下其与两相机光心连线分别交照片于p1p2,这2个像点即为该空间点在照片上的投影点。此时,O1O2P三点构成的平面被称为极平面,而极平面与2个像平面之间的交线L1L2为极线。O1O2连线被称为基线,基线与两帧图像平面分别交于e1e2,它们被称为极点。

dfd7643055b9297d50dd47a28e5457ac.png
图 2 对极几何约束

对极约束可以表示为

cfc943f14f3523bfd4b9b9b700836c3e.png(1)

式中,x1x2为空间点的相机归一化平面坐标;p1p2为空间点的像素坐标;本质(essential)矩阵E=t(×)R,表达的是匹配成功的点的相机归一化坐标间的关系;基础(fundamental)矩阵F=K-Tt(×)RK-1,表达的则是匹配成功的点的像素坐标之间的关系[9-10]

考虑E矩阵和F矩阵仅相差一个已知的内参矩阵K,并且E的形式更为简单,因此通常情况下,利用E矩阵来恢复位姿。虽然E只有5个自由度,但是其内在性质是非线性的,为了线性求解时方便,只考虑单目视觉的尺度等价性,减少了一自由度,可利用八点法进行求解。

设相机归一化坐标系下的点坐标为

83513ed1aa2d60a73e0cabf40096f9d8.png

根据对极约束可得

d01572b3c587c934e4d46afab35c5083.png(2)

E展开成矢量形式为

62e953a5c5686d034c02dee3659886a4.png

则式(2)可以整理为

ddc41492e397dcf0efb63cfa679343ef.png

每有一对匹配点可以列一个等式,至少有8对点时,可以得到一个线性方程组,进而求解E矩阵。之后利用SVD分解便可以得到相机间相对位姿Rt

上述方法求得的t尺度未知,还需要利用已知点间的距离恢复真实尺度。设首帧图像的相机坐标系为参考坐标系,则

e680b49ac63220c0aa425d784ef66fb7.png

由于RTK测得的坐标(北向X、东向Y、高程h)不符合右手坐标系,为后续计算方便,本文定义地理坐标系为北向X、东向Y、垂直向下Z(-h), 则相机光心在地理坐标下表示为

c1e7f77c3336f3a672c4c9eb37dcd6c7.png(3)

其中h表示已知点的高程加上相机架设高度,则尺度因子为

aff8ac46f76b4e4386af4b3f8ebab895.png

最后得到真实尺度下的相机位姿

ec249ef8c723340d5aca518ab9b61a0b.png(4)

两相机光心在具有真实尺度的相机坐标系下的坐标分别为

f646de6598f8ef336090e4c5c337d326.png(5)

1.3 三角化

得到两帧图像的相对位姿后,需要三角化出待测点在相机坐标系下的三维坐标。设第一帧图像的相机坐标系为全局坐标系,则其旋转平移矩阵分别为

d3a13302c43e5421823c5863ce207d91.png

式中,I为单位阵;0为零向量。根据式(4),第二帧图像的旋转平移分别为

3135d78006e52eb13af5b8c3a16be52b.png

设相机系下待测点为PC,其齐次坐标为PC=[X Y Z 1]T,该点在两帧图像上像素坐标的齐次形式分别为p1=[u1v1 1]vp2=[u2v2 1]T,在齐次坐标下的投影关系可表示为

db1c2bb39540c98f549f3dc9cfd23ef7.png

对于此两帧图像,该关系可表示为

541b98e44ffeb5583d685382400684f6.png(6)

线性求解式(6)可得PC

1.4 计算相机坐标系与地理坐标系间的旋转

静基座解析初始对准[11]利用重力矢量和地球自转矢量分别在IMU坐标系和当地导航坐标系下的投影,可以计算两个坐标系间的相对旋转。

受其启发,本算法中利用重力矢量和相机光心间的矢量分别在地理坐标系和相机坐标系下的投影求解这两个坐标系间的旋转。相机光心间的矢量根据式(3)和式(5)可作如下定义

95b0a4152061727538c2d0f2d1b018fe.png

式中,pWpC分别为相机光心矢量在地理系和相机系下的投影。

对于重力矢量,IMU加速度计输出gb在相机坐标系下的投影为

a7d7edce336998be70d06209dc9f44bd.png

式中,RbC为IMU坐标系相对于相机坐标系的旋转矩阵,即标定得到的IMU相机间的外参。而重力矢量在地理坐标系下的投影为gW=[0 0-g]T

为了使线性计算简单,设

1850dbd7e773d102ac807f62051547ae.png

则有如下关系

4ccb9e3e4ed15dda6fca0ec8df02e6ad.png(7)

求解式(7)可得

74c1c5f99170bb76ab116ae7c57b6754.png

对线性解得的RWC正交化后转置可得

12733d5f8b6addeb233ecf7a7b7a26b9.png(8)

求得相机坐标系相对于地理坐标系的旋转矩阵,则地理坐标系下的待测点的坐标可表示为

bf1e32557389361d90c620c4a9d0820b.png(9)

2 试验结果分析

试验事先用RTK测得3个点的准确坐标,三点之间的两两距离大致为10m,坐标见表 1。设点1、2为已知点,点3为待测点。

表 1 RTK点精确坐标 

571eb76fb41f14536d4e5e30d01c5320.png

试验所用相机为小觅相机深度版,其集成了双目相机和IMU。本文试验只用其左目图像和IMU数据。试验中,在已知点上方距地面1m处架设相机,拍摄图像的同时采集IMU数据。

2.1 不同特征匹配方法效果对比

两帧图像间的特征匹配效果优劣直接影响本质矩阵的计算,并最终影响到后面三角化的精度,考虑KLT光流[12]跟踪只适用于小视差,本文尝试采用基于特征点的方式进行特征匹配,对比了ORB[13]、SIFT[14]和人工选点的特征匹配情况。

进行ORB和SIFT特征提取时每帧图像均提取2000个特征点,均采用暴力匹配的方式进行特征匹配,并利用随机采样一致性[15](random sample consensus,RANSAC)方法筛选匹配点;而对于人工选点,要求分布均匀的提取两张图像的特征大概15个。最后效果如图 3—图 5所示。

17148b6ae84fd11d0c680a54e4bbe98a.png
图 3 ORB特征匹配
e98ab79df2d88220efbb6f811fa4b913.png
图 4 SIFT特征匹配
062c3209c04227d09d3c838123b0454a.png
图 5 人工选取匹配特征点

由于此场景两帧图像的视差过大,特征匹配算法表现不尽如人意。ORB误匹配非常多;SIFT由于其描述子的计算方式更加复杂,描述子中包含了更多的信息,其匹配效果相对于ORB要好一些,但是特征分布不均匀,匹配到的几乎全是远处的点;而人工选点特征点数量虽少,但匹配精度很高,也可以人为控制使匹配点的分布较为均匀。因此,本文试验均采用人工选点的方式进行特征匹配。

2.2 试验结果精度分析

试验共进行了10次,定位结果见表 2,计算坐标各方向误差如图 6所示,平面和垂向误差如图 7所示。

表 2 定位结果 

975b09d926ab5174b8f3ad76dd150f10.png

bfe6decd30686aa6238a5543b3656c28.png
图 6 坐标误差
61ad661c22d6f1c299a91ca89b263699.png
图 7 平面和高程误差

综合表 2、图 6和图 7可知,北向和东向坐标误差最大值在0.15m左右,平面坐标误差几乎都小于0.2m,平面坐标误差的平均值为0.124m。相对而言,高程误差较平面坐标误差大一些,但也几乎在0.3m以下,高程误差平均为0.2m。计算得到的平面坐标分布情况如图 8所示。

43c614b89ee9e69c8becc44c083b20bb.png
图 8 点位平面坐标误差分布

图 8中,四处分布的方块表示利用本文算法计算得到的平面坐标的误差,外侧的虚线是以中心点为圆心、0.2m为半径的圆,内测实线圆则表示这10次试验的平面坐标误差的平均值,半径为0.124m。通过图 8可以更直观地看出这10组试验的定位效果,定位精度较高。

3 结语

本文针对管线测量中GNSS RTK由于信号遮挡、反射或测站不易架设等原因无法提供被测点精密坐标的场景,提出了一种IMU辅助下的单目视觉坐标传递方法。该方法通过在两个已知的RTK点上架设相机拍摄照片的同时采集IMU数据,计算后可得待测点在地理坐标系下的坐标,其平面坐标的平均误差为0.12m,高程的平均误差为0.2m,满足测量要求。

本文提出的方法相比于传统的全站仪或激光测距仪交会的方法,设备轻便、操作简单,且一次采集能够在事后获得共视区域内任何待测点的坐标,效率高,复测率低。故该方法具有独特的应用优势。

引文格式:

郭若南, 蒋郡祥, 牛小骥, 等. IMU辅助下的单目视觉坐标传递[J]. 测绘通报,2020(3):7-11. DOI: 10.13474/j.cnki.11-2246.2020.0068. 666309e7a4266a0d2d476aeb4cbf22cc.png

GUO Ruonan, JIANG Junxiang, NIU Xiaoji, et al. Monocular vision coordinate transfer method with IMU aiding[J]. Bulletin of Surveying and Mapping, 2020(3): 7-11. DOI: 10.13474/j.cnki.11-2246.2020.0068.666309e7a4266a0d2d476aeb4cbf22cc.png

———————&    End    &————————

微信投稿邮箱:imiaoh@163.com

欢迎加入《测绘通报》作者QQ群:  435838860

进群请备注:姓名+单位

  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
根据左右相机中成像坐标求解空间坐标的过程称为立体视觉三角化。具体步骤如下: 1. 标定相机,得到相机内参和外参矩阵。 2. 对左右两个相机拍摄到的同一场景进行特征点匹配,得到左右图像中的对应点。 3. 根据匹配点对,计算左右相机之间的基础矩阵或本质矩阵。 4. 对于一个目标点,假设其在左右相机中的成像坐标分别为(x_L, y_L)和(x_R, y_R),则可以通过以下公式计算其在空间中的坐标(X, Y, Z): X, Y, Z = triangulate(x_L, y_L, x_R, y_R, P_L, P_R) 其中,P_L和P_R分别为左右相机的投影矩阵,triangulate()函数为三角化函数,可以使用OpenCV等库中的函数实现。 需要注意的是,在进行三角化计算时,需要进行坐标系变换,将左右相机坐标系统一起来。 以下是Python代码示例: ```python import cv2 import numpy as np # 左右相机内参矩阵和投影矩阵 K_L = np.array([[fx, 0, cx], [0, fy, cy], [0, 0, 1]]) K_R = np.array([[fx, 0, cx], [0, fy, cy], [0, 0, 1]]) R = np.array([[r11, r12, r13], [r21, r22, r23], [r31, r32, r33]]) t = np.array([tx, ty, tz]) P_L = np.dot(K_L, np.hstack((np.identity(3), np.zeros((3, 1))))) P_R = np.dot(K_R, np.hstack((R, t.reshape(3, 1)))) # 左右相机中的匹配点对 pts_L = np.array([[x1_L, y1_L], [x2_L, y2_L], ... [xn_L, yn_L]]) pts_R = np.array([[x1_R, y1_R], [x2_R, y2_R], ... [xn_R, yn_R]]) # 计算基础矩阵或本质矩阵 F, mask = cv2.findFundamentalMat(pts_L, pts_R, cv2.FM_RANSAC) E = np.dot(np.dot(K_R.T, F), K_L) # E = np.dot(np.dot(R, np.diag([1, 1, 0])), R.T) # 如果已知R,则可直接计算本质矩阵 # 三角化计算 def triangulate(x_L, y_L, x_R, y_R, P_L, P_R): A = np.zeros((4, 4)) A[0, :] = x_L * P_L[2, :] - P_L[0, :] A[1, :] = y_L * P_L[2, :] - P_L[1, :] A[2, :] = x_R * P_R[2, :] - P_R[0, :] A[3, :] = y_R * P_R[2, :] - P_R[1, :] _, _, V = np.linalg.svd(A) X = V[-1, :3] / V[-1, 3] return X X = triangulate(x_L, y_L, x_R, y_R, P_L, P_R) ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值