CN102314674B - 一种地面激光雷达数据纹理影像配准方法 - Google Patents
一种地面激光雷达数据纹理影像配准方法 Download PDFInfo
- Publication number
- CN102314674B CN102314674B CN 201110250748 CN201110250748A CN102314674B CN 102314674 B CN102314674 B CN 102314674B CN 201110250748 CN201110250748 CN 201110250748 CN 201110250748 A CN201110250748 A CN 201110250748A CN 102314674 B CN102314674 B CN 102314674B
- Authority
- CN
- China
- Prior art keywords
- image
- partiald
- matrix
- laser radar
- unique point
- 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
Images
Landscapes
- Length Measuring Devices By Optical Means (AREA)
Abstract
本发明公开了一种地面激光雷达数据纹理影像配准方法,用于地面激光雷达点云和拍摄影像的融合,包括以下步骤:步骤一,在激光扫描影像形成的地面激光雷达点云上选取多个特征点A1~An,并且在拍摄影像上选取与多个特征点A1~An相对应的B1~Bn;步骤二,在被扫描的物方设定一个三维坐标原点,确定地面激光雷达点云上的各特征点A1~An的三维坐标(X,Y,Z),并且将拍摄影像的焦距设定为各特征点B1~Bn的第三维坐标(x,y,-f);重心化坐标,并计算缩放系数,以及计算正交矩阵R的配准参数;以上述参数为初始值,应用最小二乘迭代法解算配准参数的精确值。本发明在计算出上述参数的情况下,通过缩放、旋转和平移,就能够将拍摄影像融合入激光扫描影像中,形成彩色点云模型。
Description
技术领域
本方法主要应用在地面激光雷达点云和影像的融合,适用于任意拍摄角度的影像对于地面激光雷达点云的配准。
背景技术
激光扫描技术是上世纪九十年代中期开始出现的一项高新技术,它通过高速激光扫描测量的方法,大面积高分辨率地快速获取被测对象表面的三维坐标数据,广泛应用于模型重建、古建保护中,并逐渐成为三维城市数据模型获取的一种重要方法。如何获取激光雷达数据的纹理一直是该领域研究的重点。虽然很多地面激光扫描仪可以通过内置相机同步获取扫描点的纹理,但是其分辨率、摄影方式等不能满足应用的需要。目前,主要还是通过激光雷达数据与摄影影像的配准获取其纹理特征。这里的配准等同于摄影测量中的定向。对于3维与2维的配准,已经有很多研究成果,共线方程解法,基于共角条件的椎体解法,直接线性变换解法,基于罗德里格矩阵的直接接法,单位四元数法,奇异矩阵分解法等。地面激光雷达数据和其光学影像的配准存在着特殊性,首先,对于建筑物来说,对其内外、屋顶和地面分别进行摄本方法以点特征为配准基元,以罗德里格矩阵为旋转矩阵,采用分布解法进行配准。首先应用重心化特征点坐标计算缩放系数;然后应用重心化空间相似变换模型和罗德里格矩阵计算配准参数的初始值;最后应用共线方程最小二乘迭代计算配准参数的精确值。该方法需要的配准特征少,参数精度高,稳健性强,既避免了大角度不收敛问题,又避免了配准特征共面引起的参数不稳定的问题,适合任意影像拍摄角度和配准特征共面的激光雷达与纹理影像的配准。影时的,其角度是任意的,很多情况是大角度摄影。其次,某些扫描对象比如壁画,配准特征分布会出现共面现象。对于以上两个问题,现有的方法不能够完全解决。
发明内容
为了解决上述问题,本发明提供了一种地面激光雷达数据纹理影像配准法,用配准方于地面激光雷达点云和拍摄影像的融合,其特征在于,包括以下步骤:
步骤一,在激光扫描影像形成的地面激光雷达点云上选取多个特征点A1~An,并且在拍摄影像上选取与所述多个特征点A1~An相对应的B1~Bn,以形成A1/B1~An/Bn的n对特征点;
步骤二,在被扫描的物方设定一个三维坐标原点,根据该基准点确定地面激光雷达点云上的各特征点A1~An的三维坐标(X,Y,Z),以及确定拍摄影像上各特征点B1~Bn的二维坐标,并且将拍摄影像的焦距设定为各特征点B1~Bn的第三维坐标(x,y,-f);
步骤三,将各特征点A1~An的三维坐标(X,Y,Z)重心化为(Xm,Ym,Zm),将各特征点B1~Bn的三维坐标(x,y,-f)重心化为(xm,ym,0);
步骤四,根据步骤三中的重心化坐标,计算各特征点B1~Bn的三维坐标(x,y,-f)相对于各特征点A1~An的三维坐标(X,Y,Z)的缩放系数,并最终计算出拍摄影像相对于激光扫描影像的缩放系数λ;
步骤五,将拍摄影像相对于激光扫描的点云空间姿态设为罗德里格矩阵R,应用重心化空间相似变换模型和罗德里格矩阵计算配准参数;
步骤六,根据正交矩阵R和缩放系数λ计算出摄影像相对于激光扫描影像的平移矩阵参数(Xs,Ys,Zs);
步骤七,根据缩放系数λ、配准参数和平移矩阵参数,将地面激光雷达点云和拍摄影像融合。
优选的是,所述的地面激光雷达数据纹理影像配准方法中,在步骤四中,计算缩放系数λ的方式如下:
1)首先计算各特征点B1~Bn相对于各特征点A1~An的缩放系数λi每个特征点的缩放系数λi的计算公式为:
其中,(Xim,Yim,Zim)和(xim,yim,0)为特征点重心化坐标,i为标记特征点的标号,i从1~n,
2)随后计算影像缩放系数λ
影像缩放系数λ的计算公式为:
优选的是,在所述的地面激光雷达数据纹理影像配准方法中,所述步骤五中的重心化的空间相似变换模型为:
其中,(Xm,Ym,Zm)和(xm,ym,0)为重心化坐标,λ为缩放系数,R为标记拍摄影像相对于激光扫描点云空间姿态的设为罗德里格矩阵的正交矩阵。
优选的是,在所述的地面激光雷达数据纹理影像配准方法中,所述步骤五中的配准参数通过重心化空间相似变换模型和罗德里格矩阵计算后,得到初始值,再通过共线方程的最小二乘迭代,得到精确值。
优选的是,在所述的地面激光雷达数据纹理影像配准方法中,最小二乘迭代的终止条件为以配准参数标记的影像角元素改正数两次运算差值小于0.1秒(即为千分之一度)。
优选的是,在所述的地面激光雷达数据纹理影像配准方法中,所述步骤五中正交矩阵R被设定为
其中,a、b和c为配准参数。
优选的是,在所述的地面激光雷达数据纹理影像配准方法中,所述步骤五中正交矩阵R被设定为
其中,a、b和c为配准参数。
优选的是,在所述的地面激光雷达数据纹理影像配准方法中,其特征在于,所述共线方程为:
其中,a1、a2和a3分别依次代表R矩阵中的第一列从上向下的数值;b1、b2和b3分别依次代表R矩阵中的第二列从上向下的数值;c1、c2和c3分别依次代表R矩阵中的第三列从上向下的数值。
优选的是,在所述的地面激光雷达数据纹理影像配准方法中,其特征在于,对应特征点至少为4对。
通过本发明的方法,能够实现任意拍摄角度的影像对于地面激光雷达点云的配准。
附图说明
图1为物方坐标系和拍照示意图
具体实施方式
1.主要实施步骤
1.1确定配准参数初始值的原理及过程
本方法初始值的确定是通过重心化坐标化简空间相似变换模型,并应用反对称矩阵和罗德里格矩阵的性质对重心化的模型进行化简,直接得到计算配准参数的模型,并在有多余观测值的情况下进行最小二乘解算。
1.1.1重心化的空间相似变换模型
像空间坐标系相对于物方坐标系之间空间相似变换的一般模型为:
其中,(X,Y,Z)物方点的物方空间坐标,λ为坐标系间的缩放系数,R为两个坐标系间的设为罗德里格矩阵的旋转矩阵,(x,y,-f)为像空间坐标系坐标,(XS,YS,ZS)为摄站点的物方空间坐标。
坐标的重心化是摄影测量中经常采用的一种数据预处理的方法,用重心化坐标进行结算不仅可以减少坐标在计算中的有效位数,提高计算精度,还可以使法方程式的系数简化,个别项的数值为零,从而加快计算速度。将像点和物方点坐标进行重心化后,可以化简空间变换方程。此时,两空间坐标系的转换方程为式(1)变换为:
其中,(Xm,Ym,Zm)和(xm,ym,0)为重心化坐标。
1.1.2缩放系数λ的解算
对于传统摄影而言,整体地面模型相对于飞行高度可以近似看成平面,缩放系数应用摄影焦距和平均航高来计算。对于近景摄影,有的被摄物体的景深较大,每个点的缩放系数都不同,本文应用各个点缩放系数和的平均值来确定缩放系数。应用公式(2)和旋转矩阵为正交阵的特性,得到每个点的缩放系数的公式为:
其中,(Xim,Yim,Zim)和(xim,yim,0)为配准点重心化坐标。
则其影像λ为:
1.1.3基于罗德里格矩阵的配准参数的确定
设反对称矩阵为
R为正交矩阵,根据正交矩阵和反对称矩阵的运算性质有:
R=(I+S)(I-S)-1 (6)
把(5)式代入(6)式得到下式:
其中Δ=1+a2+b2+c2,此矩阵为罗德里格矩阵。传统的变换中,旋转矩阵在运算过程中,旋转矩阵由空间三个方向的角度函数构成,是非线性约束,经过线性化步骤后表达式特别复杂,计算效率也比较低。而组成罗德里格矩阵只需要进行四则运算,计算速度快,应用最广。
将(6)(7)式代入(2)式,得到如下方程式:
当有多余观测时,关于(a,b,c)的误差方程为:
V=AX-L (9)
其中令
求解误差方程式可得:
X=(ATA)-1ATL
根据(1)式求解平移向量:
此处f为相机焦距摄影时记录值。由于缩放系数和f存在误差,此时解算出的配准参数只能作为初始值。
1.2确定配准参数精确值的原理及过程
根据以上方法得到了配准参数的初始值,为共线方程的最小二乘迭代提供了先决条件,解决了大角度摄影带来的迭代计算不收敛的问题。本文最小二乘迭代模型不同于现有方法,应用罗德里格系数计算出影像外方位角元素的初始值,对角度进行迭代运算,而是推导了基于罗德里格参数的误差模型,这样更有利于方程解的稳定性,并且提高运算速度。共线方程公式为:
其中,a1、a2和a3分别依次代表R矩阵中的第一列从上向下的数值;b1、b2和b3分别依次代表R矩阵中的第二列从上向下的数值;c1、c2和c3分别依次代表R矩阵中的第三列从上向下的数值。
设基于共线方程的配准参数的误差方程为:
V=BX-L (11)
此时X=[ΔXS,ΔYS,ΔZS,Δa,Δb,Δc]
常数项为:L=[x-(x),y-(y)]T
根据罗德里格矩阵(6)式和共线方程的公式,推导出关于(a,b,c)误差方程系数项矩阵为:
其中Δ=1+a2+b2+c2
其它误差方程系数项与传统共线方程最小二乘迭代相同,此处不多做说明。以影像角元素改正数两次运算差值为0.1’为迭代的终止条件。为了满足初始值和精确值的解算,本方法最少需要4对配准特征。因为地面激光雷达分辨率一般在1cm左右,影像畸变引起的物方误差在分辨率范围内,所以本方法把内方位元素及畸变参数视为已知值,此处共线方程的最小二乘迭代在特征共面的情况下依然很稳定,避免了特征共面引起的配准参数解算不稳定的问题。
2.实例结果
2.1实验一
该实验应用一组大角度影像数据,分别应用本方法、直接线性变换方法和传统后方交会方法在精度、迭代次数进行比较,验证本方法在精度和速度方面的优越性。为了满足直接线性变换方法对特征点的要求,分别在影像和点云上选取8对同名特征点。实验结果如表1所示,本方法应用分步解法解决了大角度影像配准的问题,以罗德里格矩阵为旋转矩阵避免了大量的三角运算。本文方法较之直接线性变换方法和传统的后方交会方法在定向精度、迭代次数和计算速度上都有所提高。传统的后方交会方法在大角度的情况下,设其初始值为零,迭代不收敛。
表1本文方法与其它配准方法结果对照表
2.2实验二
实验二选择了某建筑内东西南北和屋顶5组数据,分别选用4对特征点进行配准,摄影大概方位如图1所示。相对于物方坐标系来说有大角度影像,有小角度影像,证明本方法的稳定性。5组数据通过应用本文方法都得到了符合精度的配准结果,结果如表2所示,并通过纹理映射及其可视化,证明了结果的正确性。本方法具有很好的稳定性,符合任意角度的激光雷达和纹理影像的配准。
表2应用室内不同方向的数据进行定向,证明该方法的稳健性
2.3实验三
实验三取用了两组控制点数据,应用本文方法和直接线性变换方法,验证配准特征共面的情况,其结果如表3所示。直接线性变换严禁所用的控制点布设在同一个平面内,否则会引起解得不稳定性,主要原因是误差方矩阵不是列满秩矩阵,列矩阵之间具有相关性。但是应用本文的方法对控制点的分布没有要求。第一组数据选择了8对分布在不同平面上的控制点,应用两种方法得到的配准参数经纹理映射可视化,结果正确。第二组数据从这8对点中去掉了两对点,此6对点趋近分布在同一个平面内,本文方法配准参数正确,直接线性变换方法虽然也收敛,但是引起了解的不稳定,可视化结果错误。实验证明,本方法对配准特征分布没有要求,稳定性强。
表3同名点对在同一个平面内和不在同一个平面内本方法和DLT方法的比较
尽管本发明的实施方案已公开如上,但其并不仅仅限于说明书和实施方式中所列运用,它完全可以被适用于各种适合本发明的领域,对于熟悉本领域的人员而言,可容易地实现另外的修改,因此在不背离权利要求及等同范围所限定的一般概念下,本发明并不限于特定的细节和这里示出与描述的图例。
Claims (4)
1.一种地面激光雷达数据纹理影像配准方法,用于地面激光雷达点云和拍摄影像的融合,其特征在于,包括以下步骤:
步骤一,在地面激光雷达形成的地面激光雷达点云上选取多个特征点A1~An,并且在拍摄影像上选取与所述多个特征点A1~An相对应的B1~Bn,以形成A1/B1~An/Bn的n对特征点;
步骤二,在被扫描的物方设定一个三维坐标原点,根据该坐标原点确定地面激光雷达点云上的各特征点A1~An的三维坐标(X,Y,Z),以及确定拍摄影像上各特征点B1~Bn的二维坐标,并且将拍摄影像的焦距设定为各特征点B1~Bn的第三维坐标,从而得到各特征点B1~Bn的三维坐标(x,y,-f);
步骤三,将各特征点A1~An的三维坐标(X,Y,Z)重心化为(Xm,Ym,Zm),将各特征点B1~Bn的三维坐标(x,y,-f)重心化为(xm,ym,0);
步骤四,根据步骤三中的重心化坐标,计算各特征点B1~Bn的三维坐标(x,y,-f)相对于各特征点A1~An的三维坐标(X,Y,Z)的缩放系数,并最终计算出拍摄影像相对于激光扫描影像的缩放系数λ;
步骤五,将拍摄影像相对于激光扫描的点云空间姿态设为罗德里格矩阵R,该矩阵为正交矩阵,应用重心化空间相似变换模型和罗德里格矩阵计算配准参数中的罗德里格矩阵参数a,b,c;
步骤六,再根据正交矩阵R和缩放系数λ计算出拍摄影像相对于激光扫描影像的配准参数中的平移参数Xs,Ys,Zs;
步骤七,以求出的配准参数作为初始值,依据共线方程的最小二乘迭代求解配准参数的精确值;
步骤八,根据配准参数精确值和共线方程,将地面激光雷达点云和拍摄影像融合;
其中,所述步骤五中的重心化的空间相似变换模型为:
其中,(Xm,Ym,Zm)和(xm,ym,0)为重心化坐标,λ为缩放系数,R为标记拍摄影像相对于激光扫描点云空间姿态的正交矩阵;
其中,所述步骤五中的配准参数通过重心化空间相似变换模型和罗德里格矩阵计算后,得到初始值,再通过共线方程的最小二乘迭代,得到精确值,迭代中罗德里格矩阵参数a,b,c法方程的系数为:
其中Δ=1+a2+b2+c2;
其中,所述步骤五中正交矩阵R被设定为:
其中,a、b和c为配准参数,也称作罗德里格矩阵参数;
其中,所述共线方程为:
其中,a1、a2和a3分别依次代表R矩阵中的第一列从上向下的数值;b1、b2和b3分别依次代表R矩阵中的第二列从上向下的数值;c1、c2和c3分别依次代表R矩阵中的第三列从上向下的数值。
2.如权利要求1所述的地面激光雷达数据纹理影像配准方法,其特征在于,在步骤四中,计算缩放系数λ的方式如下:
1)首先计算各特征点B1~Bn相对于各特征点A1~An的缩放系数λi
每个特征点的缩放系数λi的计算公式为:
其中,(Xim,Yim,Zim)和(xim,yim,0)为特征点重心化坐标,i为标记特征点的标号,i从1~n,
2)随后计算影像缩放系数λ
影像缩放系数λ的计算公式为:
3.如权利要求1所述的地面激光雷达数据纹理影像配准方法,其特征在于,最小二乘迭代的终止条件为以配准参数标记的影像角元素改正数两次运算差值小于0.1秒。
4.如权利要求3所述的地面激光雷达数据纹理影像配准方法,其特征在于,对应特征点至少为4对。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201110250748 CN102314674B (zh) | 2011-08-29 | 2011-08-29 | 一种地面激光雷达数据纹理影像配准方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201110250748 CN102314674B (zh) | 2011-08-29 | 2011-08-29 | 一种地面激光雷达数据纹理影像配准方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN102314674A CN102314674A (zh) | 2012-01-11 |
CN102314674B true CN102314674B (zh) | 2013-04-03 |
Family
ID=45427814
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN 201110250748 Active CN102314674B (zh) | 2011-08-29 | 2011-08-29 | 一种地面激光雷达数据纹理影像配准方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN102314674B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103017653B (zh) * | 2012-11-27 | 2015-06-24 | 武汉海达数云技术有限公司 | 一种球面全景影像与三维激光扫描点云的配准及测量方法 |
CN104964669B (zh) * | 2015-06-05 | 2017-07-07 | 北京建筑大学 | 类柱面文物对象正射影像生成方法 |
CN106910238A (zh) * | 2017-01-18 | 2017-06-30 | 北京建筑大学 | 基于大倾角近景影像的彩色纹理重建方法 |
CN107330929B (zh) * | 2017-06-08 | 2019-11-08 | 三峡大学 | 一种基于几何重心和质心距离比不变性的多尺度点云配准方法 |
CN110288636B (zh) * | 2019-05-05 | 2020-02-18 | 中国矿业大学 | 一种基于平面特征约束的LiDAR点云无初值配准方法 |
CN112578396B (zh) * | 2019-09-30 | 2022-04-19 | 上海禾赛科技有限公司 | 雷达间坐标变换方法及装置、计算机可读存储介质 |
CN110703245B (zh) * | 2019-10-15 | 2021-08-17 | 北京理工大学 | 基于同名点匹配与dem辅助的地基sar多角度图像配准方法 |
CN111949925B (zh) * | 2020-06-30 | 2023-08-29 | 中国资源卫星应用中心 | 基于罗德里格矩阵和最大凸包的影像相对定向方法及装置 |
CN114972446A (zh) * | 2021-12-23 | 2022-08-30 | 通用技术集团工程设计有限公司 | 一种基于罗德里格矩阵的光学影像与地面激光点云配准方法 |
Family Cites Families (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
GB0711531D0 (en) * | 2007-06-15 | 2007-07-25 | Qinetiq Ltd | Radar coordinate registration |
CN101533529B (zh) * | 2009-01-23 | 2011-11-30 | 北京建筑工程学院 | 基于深度图像的三维空间数据处理方法与装置 |
CN102062860A (zh) * | 2009-11-18 | 2011-05-18 | 中国科学院遥感应用研究所 | 基于单木位置与地表信息的地基激光雷达数据配准方法 |
CN101825442A (zh) * | 2010-04-30 | 2010-09-08 | 北京理工大学 | 一种基于移动平台的彩色激光点云成像系统 |
-
2011
- 2011-08-29 CN CN 201110250748 patent/CN102314674B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN102314674A (zh) | 2012-01-11 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102314674B (zh) | 一种地面激光雷达数据纹理影像配准方法 | |
CN109115186B (zh) | 一种针对车载移动测量系统的360°可量测全景影像生成方法 | |
US6922234B2 (en) | Method and apparatus for generating structural data from laser reflectance images | |
CN100557634C (zh) | 一种基于双一维靶标的摄像机标定方法 | |
CN109272574B (zh) | 基于投影变换的线阵旋转扫描相机成像模型构建方法和标定方法 | |
CN109727278B (zh) | 一种机载LiDAR点云数据与航空影像的自动配准方法 | |
WO2021004416A1 (zh) | 一种基于视觉信标建立信标地图的方法、装置 | |
Chatterjee et al. | Algorithms for coplanar camera calibration | |
CN103278138A (zh) | 一种复杂结构薄部件三维位置及姿态的测量方法 | |
CN109900205B (zh) | 一种高精度的单线激光器和光学相机的快速标定方法 | |
CN105300362A (zh) | 一种应用于rtk接收机的摄影测量方法 | |
CN103364012A (zh) | 一种带约束条件的多面阵航摄仪平台检校方法 | |
CN111754462A (zh) | 一种三维弯管的视觉检测方法及系统 | |
CN101794449A (zh) | 一种标定摄像机参数的方法及装置 | |
CN105118086A (zh) | 3d-aoi设备中的3d点云数据配准方法及系统 | |
CN112270698A (zh) | 基于最邻近曲面的非刚性几何配准方法 | |
CN105469386A (zh) | 一种确定立体相机高度与俯仰角的方法及装置 | |
CN113947638A (zh) | 鱼眼相机影像正射纠正方法 | |
CN101577004A (zh) | 一种极线矫正方法、装置和系统 | |
Gong et al. | A detailed study about digital surface model generation using high resolution satellite stereo imagery | |
CN104964669B (zh) | 类柱面文物对象正射影像生成方法 | |
CN111968182A (zh) | 一种双目相机非线性模型参数的标定方法 | |
CN107806861B (zh) | 一种基于本质矩阵分解的倾斜影像相对定向方法 | |
CN112819900B (zh) | 一种智能立体摄影内方位、相对定向和畸变系数标定方法 | |
Wu | Photogrammetry: 3-D from imagery |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant |