2. 北京建筑工程学院 现代城市测绘国家测绘地理信息局重点实验室，北京 100044
2. Key Laboratory for Urban Geomatics of National Administration of Surveying, Mapping and Geoinformation, Beijing University of Civil Engineering and Architecture, Beijing 100044, China
1 引 言
激光扫描技术通过高速激光扫描测量的方法，大面积高分辨率地快速获取被测对象表面的三维坐标数据[1]，广泛应用于模型重建、古建保护中，并逐渐成为三维城市数据模型获取的一种重要方法[2]。如何获取准确的激光雷达数据纹理一直是该领域研究的重点。虽然很多地面激光扫描仪可以通过内置相机同步获取扫描点的纹理，但是其分辨率、摄影方式、精准度等不能满足应用的需要。目前，主要还是通过激光雷达数据与摄影影像的配准获取其纹理特征。这里的配准等同于摄影测量中的定向。传统近景摄影测量的影像定向应用全站仪等获取控制点的三维信息，本文是从地面激光雷达点云上获取控制点三维信息。配准同名点信息分别从点云和影像上手动选取。对于一般扫描对象而言，若扫描分辨率为1cm，按照影像与物方分辨率1∶10的一般要求，摄影影像分辨率为1mm。配准误差主要来源于控制点选点误差和地面激光雷达的仪器误差，若点云获取时的仪器误差为5mm，配准误差在5个像素内方能达到精度要求。对于大型古建筑激光雷达对象，为了获取完整的、高分辨率的纹理影像需要对其进行全方位、多角度的摄影，摄影常常是大角度；同时，由于控制点选点误差和激光点云的仪器误差，配准精度往往达不到要求。
关于三维与二维的配准，已经有很多研究成果：共线方程解法[3]、角锥体法[4]、直接线性变换解法[5, 6]、基于罗德里格矩阵的直接解法[7, 8, 9]与单位四元数法[10, 11, 12, 13]等。通过分析，共线方程解法与角锥体法需要有较好的初值，通过迭代完成计算[3]；直接线性变换解法不能应用于控制点分布或近似分布在一个平面的场合，需要6对以上同名点[6]；基于罗德里格矩阵的直接解法是通过3个控制点，用矩阵的加、减、乘、除计算出影像的外方位元素，其变换模型不够严密[7, 8]；单位四元数法虽然能够解决大旋角影像的定向问题，但其采用传统等权的最小二乘平差进行迭代求解[12]，对于本文的研究对象，不能完全满足配准精度要求。
从激光点云和影像中获取的控制点信息存在选点误差和仪器误差，被视为小粗差。传统等权最小二乘平差只是把粗差分散到其他的测量值中[14]，不能消除小粗差对配准结果的影响。选权迭代法是一种变权的最小二乘估计,它可以对粗差进行定位和去除，同时亦可对观测值进行降权处理，提高定位精度[15]。选权迭代主要有Lq迭代法、丹麦法、带权data snooping法等，其中，丹麦法比较适合于摄影测量平差[16]。
综合以上分析，本文提出了一种稳健的配准方法。应用重心化的空间相似变换模型和正交旋转矩阵与反对称矩阵的关系，推导出解算点云与影像配准角度参数的模型。然后，应用改进的丹麦法选权迭代小粗差降权的方法对初始参数进行精确配准。影像配准角度参数模型适合任意摄影角度影像的配准，应用所有控制点按最小二乘平差求解，提高配准参数初始值的准确度；改进的丹麦法选权迭代精配准可以使配准精度满足本文研究对象配准的要求，有很强的稳健性。 2 粗配准参数模型 2.1 重心化的空间相似变换模型
像空间坐标系相对于物方坐标系之间空间相似变换的一般模型为[3]
式中，(X,Y,Z)物方点的物方空间坐标；λ为坐标系间的缩放系数；R为两个坐标系间的旋转矩阵；f为相机焦距，为标称值；(x,y,-f)为物方点像点的像空间坐标系坐标；(XS,YS,ZS)为摄站点的物方空间坐标。重心化坐标可以减少坐标在计算中的有效位数，提高计算精度，还可以使法方程式的系数简化，加快计算速度[17]。将像点和物方点坐标进行重心化后，式(1)变换为
式中，(Xm,Ym,Zm)和(xm,ym,0)为重心化坐标。 2.2 缩放系数λ的解算
传统摄影测量缩放系数应用摄影焦距和平均航高来计算。近景影像景深较大，每个点的缩放系数都不同，所以分别计算每个点的缩放系数应用于后续的计算理论上更加严密。对式(2)两边取转置，再与式(2)相乘，由λ>0及R为正交旋转矩阵的特性，每对配准点的缩放系数计算公式为
2.3 粗配准参数模型的推导
设反对称矩阵为[7]
R为正交旋转矩阵，根据正交矩阵和反对称矩阵的运算性质有[7]
将式(4)、式(5)代入式(2)，推导出配准角度参数(a,b,c)的模型方程式
该模型简单，运算速度快，程序容易实现。对所有控制点按照式(6)列方程，每个控制点按照式(3)计算缩放系数，应用最小二乘平差求解。若式(6)简写为
AX-L=0
则其法方程解的表达式为
得到(a,b,c)参数后，按照式(5)计算旋转矩阵，应用旋转矩阵各个元素与角元素的关系，计算影像摄影时空间姿态角元素[3]。根据角元素和式(1)求解配准参数线元素。
式(6)数学模型稳健，可以解决任意角度影像与激光点云的配准，应用最小二乘，提高模型的严密性。由于选点误差和激光扫描仪本身的误差，此时解算出的配准参数只能作为初始值，需要进行精确配准。 3 改进丹麦法选权迭代的精确配准 3.1 选权迭代
本文配准控制点存在小粗差，传统的等权最小二乘平差精度会因此而降低[18]。选权迭代的平差从最小二乘法开始，每次平差后，应用一个与关于残差的权函数计算观测值的权值，进而参加下一次迭代。随着迭代，含有粗差观测值的权将越来越小，迭代终止时，相应的残差将直接指出残差的数值，而平差结果将不再受到粗差的影响[16]。该方法的条件是[14]
∑pivi2→min
式中，权公式为
i为迭代次数；v是改正数；f(v)是权函数。
选权迭代的权函数的选择对于整个平差模型的抗差性是十分关键的[19]。本文选用丹麦法标准化残差权函数进行选权迭代，并对该函数进行了改进。 3.2 改进的丹麦法权函数
文献[14]和[16]均指出一种适合于摄影测量平差的丹麦法权函数，公式如下
式中，v为观测值残差；p0为权因子；m0为量测值的标准偏差。
文献[20]指出，残差并不能使改正数在同一个概率水平。文献[21]指出，标准残差属于同方差母体的子样，应用标准化残差替换残差作为权函数的参数更为合理。标准化残差的一般公式为
式中，vj为标准化残差；vj为观测值残差；qjj为标准化因子，其值为Qvv的主对角元素，Qvv为余差协因数阵。当权矩阵是单位阵的时候，以式(10)计算标准化残差。但是，选权迭代在第1次迭代后，权矩阵不再是单位阵。文献[21]指出，此时的标准化因子要用QvvP矩阵的主对角外元素来计算，即
式中，gjk为QvvP矩阵中的相应元素；m为观测值个数。文献[21]指出，一个好的观测值当有较大残差的时候，会被赋予较小的权，其QvvP的行元素也随之增大，余差量值也可能增大，用式(11)计算标准残差更合理，在一定程度上起到抑制弃真的作用。
经过以上分析，改进的丹麦法权函数为
3.3 基于改进丹麦法的选权迭代
根据第2部分解算出的粗配准参数，进行基于共线方程的改进丹麦法选权迭代精确配准。迭代过程中，按照式(12)的要求，以第3次迭代为分界点，采用不同的权函数。当角元素改正数两次运算差小于0.1′时迭代终止。 4 试验结果与分析 4.1 试验数据和试验工具
本文以故宫建筑为研究对象，应用徕卡HDS6000型号扫描仪进行扫描，扫描分辨率为1cm,扫描距离约为6m,此扫描距离仪器误差约为±2mm。应用Canon5D照相机进行拍照，影像像素为4368×2912，影像像素大小约0.008mm，拍摄镜头为24～70变焦镜头。按照引言中的论述，此配准数据的精度应在±4个像素内，即±0.032mm。本文方法需经过粗配准和精配准两个过程，至少需要4对控制点。为了验证本文方法的优越性共进行了4组试验，应用C#编程语言和Visual Studio 2005开发平台进行验证。 4.2 大角度影像配准试验结果对比与分析
为了验证本文方法对大角度影像的优势，以故宫“褉赏亭”内4组大角度数据为试验对象，分别应用本文方法、直接线性变换方法在精度、迭代次数和时间方面进行比较。选取8对控制点进行配准，采用基于共线方程的等权最小二乘迭代。表 1为4组数据的配准结果。
| 配准数据 | 本文方法 | 直接线性变换方法 | ||||
| 配准参数解算结果 | 单位权中误差/mm | 迭代次数 | 配准参数解算结果 | 单位权中误差/mm | 迭代次数 | |
| 1 | XS=21.77069，YS=11.19925，ZS=4.88670，φ=1.64615，ω=0.12391，κ=-1.55774 | 0.022 | 5 | XS=21.80073，YS=11.19693，ZS=4.86654，φ=1.62685，ω=0.11994，κ=-1.53983 | 0.041 | 7 |
| 2 | XS=23.22429，YS=14.09691，ZS=4.98056，φ=2.66291，ω=-1.50452，κ=-0.47469 | 0.043 | 7 | XS=23.21396，YS=14.11326，ZS=4.10047，φ=2.64338，ω=-1.48643，κ=-0.48521 | 0.069 | 10 |
| 3 | XS=24.98769,YS=14.82660，ZS=4.94876，φ=-1.62121，ω=-0.05529，κ=1.58519 | 0.030 | 4 | XS=24.96386,YS=14.83775，ZS=4.97362，φ=-1.59874，ω=-0.04998，κ=1.57957 | 0.051 | 7 |
| 4 | XS=23.30224,YS=15.57575，ZS=4.91999，φ=-2.22736，ω=1.49753，κ=2.26139 | 0.040 | 6 | XS=23.28985,YS=15.58447，ZS=4.9167，φ=-2.22663，ω=1.48994，κ=2.25438 | 0.056 | 8 |
| 均值 | 0.0338 | 5.5 | 0.0543 | 8 | ||
从表 1中可以看出，本文方法和直接线性变换方法都可以解决影像大角度的问题，并且都不需要初始值，但从每组数据单位权中误差、迭代次数、多组数据单位权中误差均值和迭代次数均值上可以看出，本文方法在配准精度和配准效率上都优于直接线性变换变换方法。对于本文方法，从多组大角度影像配准精度上看，应用等权最小二乘迭代有的能够满足要求，有的不能满足精度要求，这与选点精度和仪器误差等因素有关。这个问题应用本文提出的改进丹麦法选权迭代来解决。 4.3 全方位影像配准试验结果与分析
为了验证本文方法对任意角度影像配准方法的稳健性，选用故宫神武门内东、西、南、北、天花5个方向的影像数据，如图 1所示。每个方向选择4张影像，分别选用4对和8对控制点，采用基于共线方程的等权最小二乘迭代进行配准。配准结果如表 2所示。从室内全方位影像配准结果中可以看出，本文方法适合于任意角度的影像配准，有很好的稳健性，通过纹理映射，可视化结果正确。对于控制点的数量，从每组影像的单位权中误差均值中可以看出，4对控制点的精度要高于8对控制点的精度，证明控制点越多，带入的误差会增大，精度会降低。从每组数据精度均值中还可以发现，相同数量控制点的配准精度大体上相同。
|
| 图 1 物方坐标系和摄影方位示意图 Fig. 1 Object space coordinates and the photograph direction |
| 摄影方向 | 4对控制点 | 8对控制点 | ||
| 本文方法配准参数角元素解算结果 | 单位权中误差及均值/mm | 本文方法配准参数角元素解算结果 | 单位权中误差及均值/mm | |
| 东 | φ=1.47532，ω=0.02094，κ=-1.55370 φ=7.76693，ω=0.01841，κ=4.64034 φ=1.45371，ω=0.02472，κ=-1.51356 φ=1.50337，ω=0.02764，κ=-1.54573 | 0.051 0.037 0.030 0.037 0.029 | φ=1.45555，ω=0.02508，κ=-1.55423 φ=1.46384，ω=0.02371，κ=-1.53281 φ=1.42555，ω=0.02649，κ=-1.57494 φ=1.48648，ω=0.02838，κ=-1.56236 | 0.056 0.041 0.043 0.031 0.035 |
| 南 | φ=-7.03995，ω=-1.60465，κ=-10.16003 φ=-0.62167，ω=-1.52710，κ=2.50316 φ=0.17249，ω=-1.51467，κ=-2.95113 φ=0.19108，ω=-1.51723，κ=2.91796 | 0.019 0.034 0.049 0.031 0.038 | φ=-6.70961，ω=-1.61086，κ=-9.83166 φ=-3.40037，ω=-1.60628，κ=-0.27820 φ=0.16868，ω=-1.52366，κ=-3.08495 φ=-0.17660，ω=-1.61241，κ=2.95358 | 0.053 0.041 0.042 0.030 0.039 |
| 西 | φ=-1.49146，ω=-0.09035，κ=1.59764 φ=-1.51732，ω=-0.08743，κ=1.54668 φ=-1.45339，ω=-0.09231，κ=1.56382 φ=-1.47532，ω=-0.08884，κ=1.61264 | 0.043 0.034 0.030 0.035 0.028 | φ=-1.52652，ω=0.01165，κ=1.58846 φ=-1.51473，ω=0.09335，κ=1.56483 φ=-1.48663，ω=0.01834，κ=1.59931 φ=-1.49446，ω=0.01203，κ=1.55482 | 0.046 0.042 0.035 0.039 0.048 |
| 北 | φ=0.29723，ω=1.71856，κ=-0.28612 φ=1.79108，ω=1.60512，κ=1.81018 φ=-2.87354，ω=1.65540，κ=2.92751 φ=6.52314，ω=1.61991，κ=-6.52912 | 0.033 0.033 0.031 0.044 0.025 | φ=0.31849，ω=1.42922，κ=-0.30505 φ=1.81235，ω=1.59422，κ=1.80446 φ=-2.94856，ω=1.64916，κ=2.96735 φ=6.67701，ω=1.60822，κ=-6.68155 | 0.038 0.046 0.059 0.026 0.061 |
| 天花 | φ=-3.16159，ω=-0.05157，κ=3.14547 φ=-3.13327，ω=0.06213，κ=-3.13072 φ=-3.12015，ω=-0.03513，κ=0.01523 φ=3.13569，ω=-0.01503，κ=0.015487 | 0.021 0.021 0.029 0.011 0.025 | φ=-3.18227，ω=-0.05949，κ=3.14283 φ=-3.09629，ω=0.06612，κ=-3.12995 φ=3.15384，ω=-0.03978，κ=0.01532 φ=3.13984，ω=-0.02161，κ=0.01526 | 0.032 0.032 0.033 0.028 0.037 |
为了验证本文提出的改进丹麦法选权迭代在提高配准精度方面的优势，应用一组大角度影像数据，在点云和影像上选取了12对控制点，并设置了8对检查点。分别应用传统的最小二乘迭代和本文提出的改进丹麦法选权迭代进行试验对比，试验结果如表 3所示。
表 2中，定向点拟合误差的公式为
Vx和Vy为定向点上的残差。检查点的公式与其相同。
从表 3试验结果可以看出，改进的丹麦法选权迭代虽然迭代次数要多于最小二乘迭代，但是配准结果要高于等权最小二乘迭代，配准精度从0.04mm提高到了0.023mm，满足了4个像素的配准要求。同时，从配准参数中误差和定向点的拟合误差和检查点拟合误差结果可以再次证明，改进丹麦法选权迭代的配准参数精度要高于最小二乘迭代，证明了选权迭代的优越性。经过大量的实际数据验证，应用改进丹麦法选权迭代精确配准，都能够使精度达到要求，实用性很强。
| 方法 | 单位权中误差/mm | 参数中误差mXS mYS mZS/mmmφmωmκ/(°) | 定向点拟合误差ξx/mmξy/mm | 检查点拟合误差ζx/mm ζy/mm | 迭代次数 |
| 最小二乘迭代 | 0.040 | 0.05145 0.14338 0.00569 0.14363 0.05092 0.05448 | 0.07343 0.01434 | 0.05447 0.08995 | 5 |
| 改进丹麦法选权迭代 | 0.023 | 0.02799 0.11254 0.00287 0.11264 0.04434 0.05020 | 0.04654 0.00815 | 0.02048 0.01836 | 35 |
本文方法通过降权处理，配准精度满足要求，这对于大型点云和全方位连续影像配准有很重要的现实意义。本组试验以两张连续影像为对象，分别应用现有的纹理贴图软件和本文方法进行影像配准和点云纹理映射可视化试验。图 2(a)和图 2(b)分别为应用现有软件和本文方法的连续影像点云纹理映射可视化结果。从图 2(a)的左侧纹理映射点云中可以明显看出有纹理接边现象，而应用本文的方法图 2(b)右侧纹理映射点云没有纹理接边现象。从模型映射局部放大图中，图 2(a)可以明显看出纹理接边，而图 2(b)的接边现象已经很小。通过大量的试验表明，连续影像在平行摄影，比例尺大致相同的基础上，对于激光点云可达到无缝纹理映射。
|
| 图 2 纹理映射可视化对比图 Fig. 2 Comparison of visual texture mapping |
本文提出了一种地面激光雷达纹理影像的稳健、高精度配准方法，解决了低分辨率点云和高分辨率影像的配准问题：
(1) 最少4对控制点即可完成配准，提高了效率；
(2) 把像主点和相机镜头的畸变视为标称值，避免了控制点在同一平面造成法方程病态问题；
(3) 通过先粗配准再精配准的策略提高精度，使大型点云和全方位连续影像达到了无缝纹理映射。
本文方法只满足点云与影像间的纹理映射，由于分辨率的差异，配准结果对点云相对应的三角网模型的纹理映射还没有达到完全无缝纹理映射，这个问题将在其他论文中予以阐述和解决，在此不做过多介绍。
| [1] | WANG Yanmin, ZHU Guang, GUO Ming et al. Digital Protection of Historic Building by Combining Terrestrial Lidar with Digital Image[C]//Optical 3D Measurement Techniques VⅢ: Ⅱ. Berne:[s.n.], 2007: 264-271. |
| [2] | VOSSELMAN G, KESSELS P, GORTE B G H. The Utilization of Airborne Laser Scanning for Three Dimensional Mapping[J]. International Journal of Applied Earth Observation and Geoinformation, 2005, 6(324): 177-186. |
| [3] | ZHANG Jianqing, PAN Li, WANG Shugen. Photogrammetry[M]. Wuhan: Wuhan University Press, 2003. (张剑清, 潘励, 王树根. 摄影测量学[M]. 武汉: 武汉大学出版社, 2003.) |
| [4] | LI Deren, ZHENG Zhaobao. Analytical Photogrammetry[M]. Beijing: Surveying and Mapping Press, 1992. (李德仁, 郑肇葆. 解析摄影测量学[M]. 北京: 测绘出版社, 1992.) |
| [5] | ZHANG Fan. High Fidelity Texture Reconstruction in 3D Modeling from Laser Scanning Data[D]. Wuhan: Wuhan University. (张帆. 激光扫描数据三维建模中高保真度纹理重建研究[D].武汉: 武汉大学, 2008.) |
| [6] | FENG Wenhao. Photogrammetry[M]. Wuhan: Wuhan University Press, 2002. (冯文灏. 近景摄影测量[M]. 武汉: 武汉大学出版社, 2002.) |
| [7] | YU Zhilu, YAO Jili, LU Changguang. The Direct-solution of Spatial-resection Based on Lodriguez Matrix[J]. Engineering of Surveying and Mapping, 2005, 14(2): 50-52. (于志路, 姚吉利, 吕长广. 罗德里格矩阵在空间后方交会直接解法中的应用[J]. 测绘工程, 2005, 14(2): 50-52.) |
| [8] | YAO Jili, SUN Yating, WANG Shuguang. The Direct Solution of Spatial Resection Based on Luodrigues Matrix[J]. Journal of Shandong University of Technology: Science and Technology, 2006, 20(2): 36-39. (姚吉利, 孙亚廷, 王树广. 基于罗德里格矩阵的数码像片直接定向的方法[J]. 山东理工大学学报: 自然科学版, 2006, 20(2): 36-39.) |
| [9] | YANG Fan, LI Guangyun, WANG Li, et al. A Method of Least Square Iterative Coordinate Transformation Based on Lodrigues Matrix[J]. Geotechnical Investigation & Surveying, 2010(9): 80-84. (杨凡, 李广云, 王力, 等. 一种基于罗德里格矩阵的最小二乘迭代坐标转换方法[J]. 工程勘察, 2010(9): 80-84.) |
| [10] | WANG Yong, JIANG Ting, JIANG Gangwu, et al. Space Resection of Single Image Based on the Description of Unit Quaternions[J]. Journal of Geomatics Science and Technology, 2007, 24(2): 133-135. (王勇, 姜挺, 江刚武, 等. 基于单位四元数描述的单像空间后方交会[J]. 测绘科学技术学报, 2007, 24(2): 133-135.) |
| [11] | LV Haiqing, ZOU Zhengrong, LUO Faming. Research on the Construction of Rotation Matrix in Close-range Photogrammetry[J]. Science of Surveying and Mapping, 2007, 32(3):15-17. (闾海庆, 邹峥嵘, 罗发明. 近景摄影测量中旋转矩阵构成方法的研究[J]. 测绘科学, 2007, 32(3): 15-17.) |
| [12] | JIANG Gangwu, JIANG Ting, WANG Yong, et al. Space Resection Independent of Initial Value Based on Unit Quaternions[J]. Acta Geodaetica et Cartographica Sinica, 2007, 36(2): 169-175. (江刚武, 姜挺, 王勇, 等. 基于单位四元数的无初值依赖空间后方交会[J]. 测绘学报, 2007, 36(2): 169-175.) |
| [13] | GUAN Yunlan, CHENG Xiaojun, ZHOU Shijian, et al. A Solution to Space Resection Based on Unit Quaternion[J]. Acta Geodaetica et Cartographica Sinica, 2008, 37(1): 30-35. (官云兰, 程效军, 周世健, 等. 基于单位四元数的空间后方交会解算[J]. 测绘学报, 2008, 37(1): 30-35.) |
| [14] | WANG Zhizuo.Photogrammetry Principle Sequel[M]. Wuhan: Wuhan University Press, 2007. (王之卓. 摄影测量原理续编[M]. 武汉: 武汉大学出版社, 2007.) |
| [15] | TANG Yanmei, WANG Jian, LIU Chao, et al. Research on Danish Robust Navigation Model Based on S/N Weighted[J]. GNSS World of China,2010(6): 12-15. (唐艳梅, 王坚, 刘超, 等. 基于信噪比加权的丹麦法抗差导航模型研究[J]. 全球定位系统,2010(6): 12-15.) |
| [16] | LI Deren. Gross Error Location by Means of the Iteration Method with Variable Weights[J]. Journal of Wuhan Technical University of Surveying and Mapping. 1984(1): 46-67. (李德仁. 利用选择权迭代法进行粗差定位[J]. 武汉测绘科技大学学报, 1984(1): 46-67.) |
| [17] | LIN Junjian, CANG Guihua. Photogrammetry[M]. Beijing: National Defence Industry Press, 2006. (林君建,苍桂华. 摄影测量学[M]. 北京:国防工业出版社,2006.) |
| [18] | WANG Renxiang. Functional Essentiality of Robust Adjustment with Weight Function[J]. Journal of People’s Liberation Army Institute of Surveying and Mapping, 1986(2): 39-47. (王任享.选权迭代剔除粗差的实质[J]. 解放军测绘学院学报, 1986(2): 39-47.) |
| [19] | LUO Rui, LIU Jianhua, FAN Dongjuan. Iterated Adjustment with Weight Function of SPOT Resection[J]. Journal of the PLA Institute of Surveying and Mapping.1998,15(3): 30-32. (罗睿, 刘建华, 范东娟. 引入选权迭代的SPOT图像后方交会[J]. 测绘学院学报, 1998, 15(1): 30-32.) |
| [20] | EARLS C J. A Precision Potential of Close-range Photogrammetry System[J]. The Photogrammetric Record, 1983, 62(1):14-23. |
| [21] | WANG Renxiang. Effects of Parameters of Weight Function for Blunder Detection by the Recursive Weighted Least Squares Method[J]. Journal of Wuhan Technical University of Surveying and Mapping, 1988, 13(4): 42-50. (王任享. 选权迭代定位粗差时权函数参数之功能[J]. 武汉测绘科技大学学报, 1988, 13(4): 42-50.) |