三维激光扫描技术作为快速获取场景表面三维信息的手段，正以其独有的优势广泛地应用于虚拟现实、逆向工程、历史遗迹保护等众多行业领域，具有十分广阔的应用前景。而数码影像作为丰富几何信息和表面纹理信息的数据载体，与点云数据是两种跨模态的异源数据^[1]，两种数据结合可以充分发挥各自的优势，实现复杂对象的精细化三维重建。

针对如何获取带有丰富纹理信息的全景三维点云的问题，文献[2]利用尺度不变特征变换(scale-invariant feature transform, SIFT)^[2]特征对影像序列进行匹配，然后利用从运动中恢复结构(structure from motion, SfM)算法恢复稀疏的三维点云^[3]，最后将由数码影像生成的彩色密集点云与激光点云进行配准，得到带有纹理信息的全景彩色点云^[4]。该方法自动化程度较高，可以充分利用多张影像的纹理信息，且精度高，但其效率低下^[5]。文献[6]还利用在立体像对匹配点与三维扫描点云进行人工选点配准的方法实现点云与数码影像的融合，其原理与文献[4]类似。但由于在影像生成的稀疏点云上选点较为困难，使得得到的初始外参数错误，而导致纹理映射错误。文献[7]提出了一种局部区域表面一致性约束的纹理映射方法，但其对象为曲面重建后的三维模型表面。文献[8]提出了一种全景影像序列与点云的整体配准方法，直接利用全景成像数学模型对全景影像序列进行运动恢复结构，得到彩色密集点云，将其与激光点云进行最邻近迭代(iterative closest point, ICP)配准，从而得到具有纹理信息的彩色点云。文献[9]提出了一种利用点云生成的全景强度影像与光学影像在互信息衡量标准的基础上对两者进行配准的方法。然而影响点云强度图像灰度值的因素很多，因此将点云全景强度影像像素点作为配准基元难度较大，而且点云全景强度影像与光学影像的成像原理不同，SIFT、加速稳健特征(speeded up robust features, SURF)等特征检测算子不易获取正确特征点，因此难以与数码影像进行配准。文献[10]利用互信息对车载激光点云和全景影像进行配准。另外，还有学者将多张CCD数码影像拼接成全景影像，并且根据中心投影成像原理将全景影像转化为全景球面点云，在全景影像与三维激光扫描点云上建立对应控制点，使用罗德里格矩阵进行配准^[11]，得到彩色全景点云。由于数码影像拼接需要进行影像匹配，因此影像分辨率越高，影像拼接效率越低，而且得到的全景影像易产生畸变，对后续纹理映射的精度产生较大影响。文献[12]将普通数码相机固定在高精度旋转平台上集合成旋转全景相机，通过解算每张照片的外方位元素并加上前方交会的方法获得物方点坐标，通过该方法能够方便地实现野外摄影测量，并且有较高的精度。文献[13]提出了一种基于光源归一的投影回归方法，根据鱼眼镜头的光学系统构造了一个虚拟投影半球面，通过将激光扫描仪与鱼眼镜头的光源进行归一，从而实现鱼眼图像与激光点云在虚拟投影半球面上的数据匹配。然而这种方法对鱼眼镜头要求较高，相机成本较大。

纹理映射问题的关键在于如何快速、精确地获取摄像机的位置姿态，即其外参数。国内外已有很多学者对外参数标定方法做了大量深入的研究。在近景摄影测量领域，外参数标定被称为空间后方交会问题，即利用一定数量的2D-3D点对求解相机在世界坐标系下的姿态和位置(外方位元素)。基于摄影测量原理的空间后方交会方法主要有罗德里格矩阵法^[14]以及单位四元数^[15]解法。此类方法以罗德里格矩阵元素或者单位四元数代替外方位元素，然后再进行线性化迭代求解，然而此方法需要提供合适的迭代初值。直接线性变换(direct linear transformation, DLT)解法用11个无线性无关的参数代替内外方位元素，再利用参数间数学关系求解相机内外方位元素^[16-17]。尽管DLT方法可以将内外方位元素全部解出，但稳定性不强。由于DLT解法无需提供任何初值，计算复杂度不高，故其常作为一种初值求取方法与其他方法相结合^[18]。另外，文献[19]采用共线方程的变形形式，将其用二元多项式近似表示，通过点云与影像对应控制点来求取多项式系数。这种方法所需控制点较多，且精度不高。

国外学者对外参数标定也做了大量研究。这一问题于1841年由一位德国数学家首次提出，被称为n点透视问题(perspective n point, PnP)，在计算机视觉中又称为位姿估计。针对此问题，文献[20]首次提出P3P算法，文献[21]提出的随机采样一致(random sample consensus, RANSAC)算法进一步提高了该算法的鲁棒性，然而其解不唯一的特点一直是该算法的主要问题^[22-24]；文献[25]提出一种统一n点透视(unified perspective n point, UPnP）算法，该算法不仅可以求得摄像机的姿态参数，同时也可以求出摄像机焦距，具有较高的鲁棒性，相较于其他迭代解法^[26-27]，如最小二乘、直接线性变换等，在算法效率上具有优势。

针对三维激光点云与影像融合问题，综合上述方法的优劣之处，本文在采用三维激光扫描仪获取点云数据的同时外置数码相机，进行联合标定，求取每张影像相对于世界坐标系的绝对位置姿态，从而为三维全景点云与多张影像建立几何对应关系，实现全景影像与三维激光点云的纹理映射，获取具有丰富纹理信息的彩色点云，为后续利用点云进行精细化建模提供生动的纹理信息。

对于利用三维激光扫描仪与外置数码相机获取真三维彩色点云的纹理映射问题，外置相机为固定焦距的高分辨率数码相机，即内参数可视为已知，并且影像畸变可以忽略不计，与三维激光扫描仪的初始相对姿态固定。本文首先阐述基于虚拟控制点的初始影像位姿估计的基本原理，然后将标定结果与P3P、DLT、UPnP解法进行对比，利用此方法进行多像位姿估计，继而对经过Mask匀光法^[28-29]进行匀光后的多张影像进行全景纹理映射，最后利用线性融合的方法进行影像接缝处色彩一致性处理。

1 影像位姿估计与匀光匀色算法

本文在假设摄像机内参数已知的情况下，首先获取初始影像的平面坐标系与三维世界坐标系下的2D坐标与3D坐标对应控制点，通过摄像机位姿估计计算初始影像的姿态参数，即平移量t与旋转矩阵R，其余影像的姿态参数则可以根据激光扫描仪与外置摄像机的空间旋转法则，利用初始影像姿态参数计算得到；然后利用每张影像的姿态参数实现三维激光点云的全景纹理映射。另外，影像之间存在光照、颜色等的差异，影像匀光匀色算法较好地处理了影像接缝。

1.1 摄像机模型

本文算法使用的几何数学模型如图 1所示。世界坐标系为右手坐标系，原点为O_w，坐标轴朝向分别用X_w、Y_w、Z_w表示，坐标轴之间两两垂直；相机坐标系为右手系，坐标轴原点O_c为摄影中心，三轴朝向分别为X_c、Y_c、Z_c，世界坐标系原点O_w在相机坐标系下的坐标即为待求的平移矩阵t：

$ \mathit{\boldsymbol{t}} = [\begin{array}{*{20}{c}} {{\mathit{X}_\mathit{s}}}&{{\mathit{Y}_\mathit{s}}}&{{\mathit{Z}_\mathit{s}}]} \end{array} $

(1)

任意物方控制点P_i(i = 1，2…n)在世界坐标系下的坐标为p_i^w，在摄像机坐标系下的坐标为p_i^c，用齐次矢量可以表示为：

$ \left\{ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{p}}_i^w = {{\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {X_i^w}&{Y_i^w} \end{array}}&{\begin{array}{*{20}{c}} {Z_i^w}&1 \end{array}} \end{array}} \right]}^{\rm{T}}}\left( {i = 1,2 \cdots n} \right)}\\ {\mathit{\boldsymbol{p}}_i^c = {{\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {X_i^c}&{Y_i^c} \end{array}}&{\begin{array}{*{20}{c}} {Z_i^c}&1 \end{array}} \end{array}} \right]}^{\rm{T}}}\left( {i = 1,2 \cdots n} \right)} \end{array}} \right. $

(2)

影像坐标系为笛卡尔直角坐标系，其原点在影像左上角O_o，坐标轴u、v分别平行于影像扫描行列，u₀、v₀为影像中心点坐标。P_i在图像坐标系中的投影为：

u_i、v_i均以像素为单位。由针孔相机透视模型的几何关系得：

$ {u_i} = \frac{{{x_i}}}{{{\rm{d}}x}} + {u_0},{v_i} = \frac{{{y_i}}}{{{\rm{d}}y}} + {v_0} $

(3)

其中，(x_i，y_i)为实际的影像物理坐标；dx、dy为每个像素在u、v方向的物理长度，写为矩阵形式为：

$ \left[ {\begin{array}{*{20}{c}} {{u_i}}\\ {{v_i}}\\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\frac{1}{{{\rm{d}}x}}}&0&{{u_0}}\\ 0&{\frac{1}{{{\rm{d}}\mathit{y}}}}&{{v_0}}\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_i}}\\ {{y_i}}\\ 1 \end{array}} \right] $

(4)

式（4）表达了影像物理坐标与像素坐标之间的相互转换关系。同时，根据相似三角形原理，可以由相机焦距f得到三维场景中任意点P_i的摄像机坐标与其影像物理坐标之间的转换关系：

$ {x_i} = f\frac{{X_i^c}}{{Z_i^c}},{y_i} = f\frac{{Y_i^c}}{{Z_i^c}} $

(5)

写成矩阵形式为：

$ Z_i^c\left[ {\begin{array}{*{20}{c}} {{x_i}}\\ {{y_i}}\\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} f&0&0\\ 0&f&0\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {X_i^c}\\ {Y_i^c}\\ {Z_i^c} \end{array}} \right] $

(6)

另外，三维点P_i的摄像机坐标与世界坐标之间有如下关系：

$ \left[ {\begin{array}{*{20}{c}} {X_i^c}\\ {Y_i^c}\\ {Z_i^c} \end{array}} \right] = [\begin{array}{*{20}{c}} \mathit{\boldsymbol{R}}&{\mathit{\boldsymbol{t}}]} \end{array}\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {X_i^w}\\ {Y_i^w}\\ {Z_i^w} \end{array}}\\ 1 \end{array}} \right] $

(7)

式中，R、t即为摄像机外参数标定中所要求取的未知量，包含了世界坐标系在摄像机坐标系中的旋转角度与平移量。综合上述几式，可以得到三维点坐标与其投影在影像上的像素坐标之间的转换关系：

$ Z_i^c\left[ {\begin{array}{*{20}{c}} {{u_i}}\\ {{v_i}}\\ 1 \end{array}} \right] = \mathit{\boldsymbol{K}}[\begin{array}{*{20}{c}} \mathit{\boldsymbol{R}}&{\mathit{\boldsymbol{t}}]} \end{array}\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {X_i^w}\\ {Y_i^w}\\ {Z_i^w} \end{array}}\\ 1 \end{array}} \right] = \mathit{\boldsymbol{P}}\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {X_i^w}\\ {Y_i^w}\\ {Z_i^w} \end{array}}\\ 1 \end{array}} \right] $

(8)

式中，K为摄像机内参数矩阵；P为摄像机矩阵，包含摄像机内参数矩阵与外参数矩阵，本文假设摄像机内参数矩阵已知，摄像机透视投影模型如图 1所示。

1.2 控制点主成分分析

本文首先对控制点对进行处理，将多个控制点转换为4个虚拟控制点的权重之和，求取控制点在摄像机坐标系下的坐标，继而恢复2D-3D之间的映射关系。因此，虚拟控制点的选取成为关键。由虚拟控制点求取摄像机位姿的算法原理如图 2所示。

若已知n对2D-3D控制点，其在世界坐标系中的坐标记为：

$ \mathit{\boldsymbol{p}}_i^w,i = 1,2 \cdots n $

(9)

在摄像机坐标系中的坐标记为：

$\mathit{\boldsymbol{p}}_i^c,i = 1,2 \cdots n$

(10)

同时，选取的4个虚拟控制点在摄像机坐标系与世界坐标系中的坐标分别记为：

$ \mathit{\boldsymbol{c}}_i^c,\mathit{\boldsymbol{c}}_i^w,i = 1,2,3,4 $

(11)

因此可以利用虚拟控制点的权重和来表示其对应坐标系中的控制点坐标，分别为：

$ \mathit{\boldsymbol{p}}_i^w = \sum\limits_{j = 1}^4 {{\alpha _{ij}}} c_i^w,\quad \mathit{\boldsymbol{p}}_i^c = \sum\limits_{j = 1}^4 {{\alpha _{ij}}} c_i^c $

(12)

式中，$\sum\limits_{j = 1}^4 {{\alpha _{ij}}} = 1$，${\alpha _{ij}}$为每个3D控制点用虚拟控制点表示的权重。首先利用主成分分析法(principal components analysis, PCA)确定4个虚拟控制点在世界坐标系下的坐标，即c₁^w。主成分分析是一种统计方法，通过正交变换将一组可能存在相关性的变量转换成一组不相关的变量，同时得到的这组不相关的变量能最大程度反映原变量所代表的信息。本文采用主成分分析法确定虚拟控制点，求解步骤如下：

1) 求取所有控制点$\mathit{\boldsymbol{p}}_i^w$的质心$\mathit{\boldsymbol{Pc}} = \frac{1}{n}\sum\limits_{i = 1}^n {\mathit{\boldsymbol{p}}_i^w} $作为第一个虚拟控制点$\mathit{\boldsymbol{c}}_i^w$。

2) 将控制点$\mathit{\boldsymbol{p}}_i^w$按列组成n行m列矩阵M，求取每一列的均值。

3) 将M的每一行进行零值化，即减去这一列的均值。

4) 求出M的协方差矩阵$\mathit{\boldsymbol{Q}} = \frac{1}{n}{\mathit{\boldsymbol{M}}^{\rm{T}}}\mathit{\boldsymbol{M}}$。

5) 用奇异值分解法(singular value decomposition, SVD)分解协方差矩阵$\mathit{\boldsymbol{Q}} = \mathit{\boldsymbol{U \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{V}}^{\rm{T}}}$。其中，$\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}$中对角线上的元素σ₁、σ₂、σ₃为矩阵Q的奇异值，组成V^T的列向量v₁、v₂、v₃分别为奇异值对应的右奇异向量。

6) 利用求得的奇异值与奇异向量求取其余虚拟控制点$\mathit{\boldsymbol{c}}_{i + 1}^w = \mathit{\boldsymbol{c}}_1^w + \sqrt[{}]{{{\sigma _i}}}{\mathit{\boldsymbol{v}}_i},{\rm{i}} = 1,2,3$。

在求得4个虚拟控制点$\mathit{\boldsymbol{c}}_i^w$后，即可求得每个控制点的权重系数α_ij。将式（12）展开得：

$\left\{ {\begin{array}{*{20}{c}} {X_i^w = {\alpha _{i1}}c_{1x}^w + {\alpha _{i2}}c_{2x}^w + {\alpha _{i3}}c_{3x}^w + {\alpha _{i4}}c_{4x}^w}\\ {Y_i^w = {\alpha _{i1}}c_{1y}^w + {\alpha _{i2}}c_{2y}^w + {\alpha _{i3}}c_{3y}^w + {\alpha _{i4}}c_{4y}^w}\\ {Z_i^w = {\alpha _{i1}}c_{1z}^w + {\alpha _{i2}}c_{2z}^w + {\alpha _{i3}}c_{3z}^w + {\alpha _{i4}}c_{4z}^w} \end{array}} \right.$

(13)

另有：

${\alpha _{i1}} + {\alpha _{i2}} + {\alpha _{i3}} + {\alpha _{i4}} = 1$

(14)

联立式（13）、（14）即可解得α_ij。式（13）中，c_ix^w、c_iy^w、c_iz^w分别表示每个虚拟控制点在世界坐标系中3个坐标轴上的坐标。

1.3 求解摄像机坐标系虚拟控制点

对于每一个控制点，将式（12）代入式（4）和式（6）中，可得到：

${s_i}\left[ {\begin{array}{*{20}{c}} {{u_i}}\\ {{v_i}}\\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_x}}&0&{{u_0}}\\ 0&{{f_y}}&{{v_0}}\\ 0&0&1 \end{array}} \right]\sum\limits_{j = 1}^4 {{\alpha _{ij}}} \left[ {\begin{array}{*{20}{c}} {c_{ix}^c}\\ {c_{iY}^c}\\ {c_{iZ}^c} \end{array}} \right]$

(15)

其中，${f_x} = \frac{f}{{{\rm{d}}x}}$、${f_y} = \frac{f}{{{\rm{d}}y}}\;$分别表示以u、v轴像素尺寸表示的像素焦距；c_iX^c、c_iY^c、c_iZ^c分别表示每个虚拟控制点在摄像机坐标系中3个坐标轴上的坐标；${s_i} = \sum\limits_{j = 1}^4 {{\alpha _{ij}}} c_{iZ}^c$。展开得：

$ \left\{ {\begin{array}{*{20}{l}} {\sum\limits_{j = 1}^4 {{\alpha _{ij}}} \left[ {{f_x}c_{ix}^c + c_{iz}^c\left( {{u_0} - {u_i}} \right)} \right] = 0}\\ {\sum\limits_{j = 1}^4 {{\alpha _{ij}}} \left[ {{f_y}c_{iY}^c + c_{iZ}^c\left( {{v_0} - {v_i}} \right)} \right] = 0} \end{array}} \right. $

(16)

写成矩阵形式后提取系数矩阵记为M_2n×12，4个虚拟控制点在摄像机坐标系中的坐标组成12行1列的未知数向量X_12×1，有方程：

${\mathit{\boldsymbol{M}}_{2n \times 12}}{\mathit{\boldsymbol{X}}_{12 \times 1}} = 0$

(17)

当n > 6时，方程组为超定方程组，有唯一解；当n < 6时，方程组为欠定方程组，有无数个解，M_2n×12的零空间即为方程组解的集合。由于直接对M_2n×12进行SVD分解的时间复杂度为O (n³），而对M_2n×12^TM_2n×12进行SVD分解的时间复杂度为O (n)，因此本文采用对矩阵M^TM_{12 × 12}进行SVD分解求取其右奇异向量V来求解式（17）。

当n ≤ 6时，矩阵M_2n×12的零空间记为M^TM_{12 × 12}的右奇异向量的线性组合：

${\mathit{\boldsymbol{X}}_{12 \times 1}} = \sum\limits_{j = 1}^N {{\delta _j}} {V_{j,}}j = 1,2 \cdots N$

(18)

式中，V_j为M^TM_{12 × 12}零奇异值所对应的奇异向量；N是M^TM_{12 × 12}零空间的维数，根据虚拟控制点个数知，N可取值1、2、3、4；δ_j为V_j的未知系数。由于摄像机坐标系与世界坐标系均为三维笛卡尔直角坐标系，则虚拟控制点之间的距离在两个坐标系中是相等的。据此，有如下约束关系：

$ {\left\| {\mathit{\boldsymbol{c}}_i^c - \mathit{\boldsymbol{c}}_j^c} \right\|^2} = {\left\| {\mathit{\boldsymbol{c}}_i^w - \mathit{\boldsymbol{c}}_j^w} \right\|^2},i,j = 1,2,3,4 $

(19)

将式（18）代入式（19）得：

${\left\| {{\delta _i}{\mathit{\boldsymbol{V}}_i} - {\delta _j}\mathit{\boldsymbol{V}}} \right\|_j}^2 = {\left\| {\mathit{\boldsymbol{c}}_i^w - \mathit{\boldsymbol{c}}_j^w} \right\|^2}$

(20)

则4个虚拟控制点可以列出C₄² = 6个方程。

当N = 1时，有X_{12 × 1} = δ₁V₁。只有一个未知数δ₁，则有封闭解：

${\delta _1} = \frac{{\sum\limits_{\{ i,j\} \in [1,4]} {\left\| {{\mathit{\boldsymbol{V}}_i} - {\mathit{\boldsymbol{V}}_j}} \right\|} \cdot \left\| {\mathit{\boldsymbol{c}}_i^w - \mathit{\boldsymbol{c}}_j^w} \right\|}}{{\sum\limits_{\{ i,j\} \in [1,4]} {{{\left\| {{\mathit{\boldsymbol{V}}_i} - {\mathit{\boldsymbol{V}}_j}} \right\|}^2}} }}$

(21)

当N = 2时，式（18）展开得：

${\mathit{\boldsymbol{X}}_{12 \times 1}} = {\delta _1}{\mathit{\boldsymbol{V}}_1} + {\delta _2}{\mathit{\boldsymbol{V}}_2}$

(22)

代入式（19）中，方程组中有3个未知量：δ₁₁ = δ₁²，δ₁₂ = δ₁δ²，δ₂₂ = δ₂²。因为有4个虚拟控制点，可以列出6个上述线性方程，即可组成超定方程组：

${\mathit{\boldsymbol{L}}_{6 \times 3}}{\mathit{\boldsymbol{\delta }}_{3 \times 1}} = {\mathit{\boldsymbol{\rho }}_{6 \times 1}}$

(23)

利用最小二乘法进行迭代计算即可求得方程未知数δ_{3 × 1}。

当N=3时，6个距离约束条件可以列出超定方程组：

${\mathit{\boldsymbol{L}}_{6 \times 6}}{\mathit{\boldsymbol{\delta }}_{6 \times 1}} = {\mathit{\boldsymbol{\rho }}_{6 \times 1}}$

(24)

直接对L_{6 × 6}求逆即可解得未知数。

当N=4时，可列出线性方程组：

${\mathit{\boldsymbol{L}}_{6 \times 10}}{\mathit{\boldsymbol{\delta }}_{10 \times 1}} = {\mathit{\boldsymbol{\rho }}_{6 \times 1}}$

(25)

然而此方程组的未知数个数大于方程个数，属于欠定方程组，展开后未知数向量δ有10个元素，然而实际上真正的未知数个数只有4个，即：

$\mathit{\boldsymbol{\delta }} = {\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\delta _1}}&{{\delta _2}} \end{array}}&{\begin{array}{*{20}{c}} {{\delta _3}}&{{\delta _4}} \end{array}} \end{array}} \right]^{\rm{T}}}$

(26)

因此，本文采用未知数子集近似约束的方法求解。即在δ的10个元素中，只选取δ₁₁、δ₂₂、δ₃₃、δ₄₄作为未知数进行求解，此时式（25）变为超定方程组，利用最小二乘法或奇异值分解均可进行求解。

由于精度较低，上述解法求得的未知数$\mathit{\boldsymbol{\delta }} = {\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\delta _1}}&{{\delta _2}} \end{array}}&{\begin{array}{*{20}{c}} {{\delta _3}}&{{\delta _4}} \end{array}} \end{array}} \right]^{\rm{T}}}$需要进行优化以提高精度。本文采用高斯牛顿法将上述解法求得的未知数值作为初始值，通过最小化摄像机坐标系与世界坐标系中虚拟控制点两两之间的距离差值进行迭代优化，求得最终的系数向量δ，可表示为：

$e(\mathit{\boldsymbol{\delta }}) = \sum\limits_{(i,j) \le {\rm{s}}{\rm{.t}}{\rm{.}}i < j} {\left( {{{\left\| {\mathit{\boldsymbol{c}}_i^c - \mathit{\boldsymbol{c}}_j^c} \right\|}^2} - {{\left\| {\mathit{\boldsymbol{c}}_i^w - \mathit{\boldsymbol{c}}_j^w} \right\|}^2}} \right)} $

(27)

将上文求得的δ作为初始值δ₀代入，对δ求导并变形得：

$e'\left( \mathit{\boldsymbol{\delta }} \right){\rm{\Delta }}\mathit{\boldsymbol{\delta }} = - e\left( {{\mathit{\boldsymbol{\delta }}_0}} \right) = \mathit{\boldsymbol{\rho }} - \mathit{\boldsymbol{L}}{\mathit{\boldsymbol{\delta }}_0}$

(28)

4个虚拟控制点可得到6个距离约束，即可得线性方程组：

${\mathit{\boldsymbol{A}}_{6 \times 4}}{\mathit{\boldsymbol{x}}_{4 \times 1}} = {\mathit{\boldsymbol{b}}_{6 \times 1}}$

(29)

随后对A_{6 × 4}进行正交三角分解(QR decomposition, QR)分解，得未知数：

最后进行迭代即可求得系数δ_i。

1.4 绝对定向求解摄像机姿态

在求得了每个控制点与虚拟控制点的权重系数α_ij与矩阵M_{2n × 12}零空间的线性组合系数$\mathit{\boldsymbol{\delta }} = {\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\delta _1}}&{{\delta _2}} \end{array}}&{\begin{array}{*{20}{c}} {{\delta _3}}&{{\delta _4}} \end{array}} \end{array}} \right]^{\rm{T}}}$后，即可恢复虚拟控制点在摄像机坐标系中的坐标c_j^c，结合式（12）即可求得三维控制点在摄像机坐标系中的坐标p_j^c。同时，三维控制点在世界坐标系下的坐标已知，则摄像机姿态可采用文献[30]的绝对定向法求出。

首先需分别计算三维控制点在摄像机坐标系与世界坐标系下的重心，即：

$\mathit{\boldsymbol{p}}_c^c = \frac{{\sum\limits_{i = 1}^n {\mathit{\boldsymbol{p}}_i^c} }}{n},\mathit{\boldsymbol{p}}_c^w = \frac{{\sum\limits_{i = 1}^n {\mathit{\boldsymbol{p}}_i^w} }}{n}$

(30)

然后将所有控制点重心化，即：

$\mathit{\boldsymbol{q}}_i^c = \mathit{\boldsymbol{p}}_i^c - \mathit{\boldsymbol{p}}_c^c,\mathit{\boldsymbol{q}}_i^w = \mathit{\boldsymbol{p}}_i^w - \mathit{\boldsymbol{p}}_c^w$

(31)

接着计算矩阵H：

$\mathit{\boldsymbol{H}} = \sum\limits_{i = 1}^n {\mathit{\boldsymbol{q}}_i^c} \mathit{\boldsymbol{q}}_i^{w{\rm{T}}}$

(32)

最后对矩阵H进行SVD分解得到$\mathit{\boldsymbol{H}} = \mathit{\boldsymbol{U \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{V}}^{\rm{T}}}$，则摄像机坐标系与世界坐标系之间的旋转矩阵R、平移向量t可根据下式求得：

$ \mathit{\boldsymbol{R}} = \mathit{\boldsymbol{V}}{\mathit{\boldsymbol{U}}^{\rm{T}}},\mathit{\boldsymbol{t}} = \mathit{\boldsymbol{p}}_c^c - \mathit{\boldsymbol{Rp}}_c^w $

(33)

1.5 多像外参数解算

多像外参数解算的目的在于利用初始影像的位置姿态求取其余影像的外参数。若其中一张影像的安置参数已由上述位姿估计方法求得，那么其余影像的姿态就可以根据影像数据采集的特点依次按照空间向量旋转法则直接获取。

如图 3所示，假定仪器装配没有误差，即摄像机坐标系原点与仪器坐标系(世界坐标系)原点的连线与轴重合, 摄像机旋转过程中只改变各个相机的朝向，旋转过程中摄像机坐标系与世界坐标系之间的平移量没有发生改变，即：

$ \left\{ {\begin{array}{*{20}{c}} {{X_{Sj}} = {X_s}}\\ {{Y_{Sj}} = {Y_s}}\\ {{Z_{Sj}} = {Z_s}} \end{array}} \right.\;\;\;\left( {j = 1,2 \cdots m} \right) $

(34)

式中，m为相机模型的数量。

外置数码相机旋转模型如图 4所示。假设经过位姿估计的影像为A，所求旋转矩阵为R_a，由A经过顺次旋转所得影像依次为B、C…，它们在世界坐标系下的旋转姿态依次为R_b、R_c…，易知R_b、R_c…均由R_a绕Z轴旋转所得。假设顺次旋转角度为β，则根据几何关系可以得到：

$ \left\{ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}_b} = {\mathit{\boldsymbol{R}}_Z}{\mathit{\boldsymbol{R}}_a}}\\ {{\mathit{\boldsymbol{R}}_c} = \mathit{\boldsymbol{R}}_Z^2{\mathit{\boldsymbol{R}}_a}}\\ \vdots \end{array}} \right. $

(35)

$ {\mathit{\boldsymbol{R}}_Z} = \left[ {\begin{array}{*{20}{c}} {{\rm{cos}}\beta }&{{\rm{sin}}\beta }&0\\ { - {\rm{sin}}\beta }&{{\rm{cos}}\beta }&0\\ 0&0&1 \end{array}} \right] $

(36)

式中，R_Z为绕Z轴旋转的旋转矩阵。其余影像的旋转姿态可按式（35）类推。

另外应该注意的是，在摄像机进行水平旋转一周获取多张影像时，在过摄像机主光轴的直线上会同时出现两张影像，若不对其相对于扫描仪中心的前后关系进行判别，则会得到错误的纹理信息。因此, 本文通过判断摄像机中心与点云之间的方向向量与摄像机主光轴的夹角来区别扫描仪前后两张影像。若夹角小于90°，则该三维点在扫描仪前方，直接通过摄像机透视投影模型获取其红绿蓝色彩(red green blue, RGB)信息；若夹角大于90°，则舍弃。

1.6 影像匀光匀色

数码相机拍摄的存在不均匀光照的影像可以采用如下数学模型描述。

$ I{\rm{'}}\left( {x,y} \right) = I\left( {x,y} \right) + B\left( {x,y} \right) $

(37)

式中，I′(x，y)代表亮度分布不均匀的原始数码影像；I(x，y)为人们希望获得的光照均匀的图像；B(x，y)表示背景影像。那么，可知Mask匀光处理方法首先对原始影像进行低通滤波获得近似的背景影像，然后将原始影像与背景影像进行相减运算，此过程可表示为：

${I_{{\rm{dst}}}} = {I_{{\rm{src}}}} - {I_{{\rm{blur}}}} + {I_{{\rm{offset}}}}$

(38)

式中，I_offset代表偏移量，决定了最终输出影像的平均亮度，同时保证相减后的影像像素灰度值分布在0~255之间。若希望输出影像的平均亮度值与输入影像近似相等，那么偏移量的值可取输入影像的亮度均值。

经过匀光处理后的各张数码影像辐射亮度值已经被统一，在后续纹理映射过程中，每个点云按照摄像机投影模型利用摄像机的内外参数解算三维点云投影在影像上的坐标并获取其RGB值，为点云赋上纹理信息。

实际拍摄状况下，相邻影像间有一定的重叠度(本文中为0.5)。根据透视投影原理，在重叠区域内同一个点云P_i就可以投影到相邻的多张影像上去。虽然经过匀光和辐射纠正后的各个数码影像有着相同的亮度变化趋势，但同名点仍难免会有一定的色彩差异，实际效果中表现为色彩接缝。

本文通过线性融合的方法消除相邻影像间的色彩差异。线性融合方法的基本思想为：根据该点在影像中的位置判断其色彩权重，然后对该点的色彩值进行加权处理。

假设任意点云在影像A和B上的投影为点P，点P在A与B上的像素坐标分别为(X_A，Y_A)、(X_B，Y_B)，对应像素值分别为p_A、p_B，经融合后的点云像素值为p_P，线性融合公式为：

$ {p_P} = {W_A}{p_A} + {W_B}{p_B},{W_B} = 1 - {W_A} $

(39)

式中，${W_A} = \frac{{{D_a}}}{{{D_a} + {D_b}}}$，为A对应的像素权重。当点P位于影像的右半部分时，首先计算P点分别到影像A与B右侧边界的距离S_a、S_b，当S_b > S_a时，有如下关系：

$ {D_a} = {p_{{\rm{width}}}} - {X_A},{D_b} = {X_B} $

(40)

式中，p_width为相应影像的宽度。S_b < S_a的情况可类推获得，不再阐述。实际操作中可以根据具体情况只在接缝附近一定范围内进行拉伸融合，以防止重影出现。

2 位姿估计与全景纹理映射实验 2.1 实验准备

本文使用Rigel激光扫描仪加外置数码相机两套设备获取三维点云与数码影像。首先在室内墙面上均匀布设回光反射标志，Rigel扫描仪通过联机自动控制分别进行点云数据采集和数码影像获取，得到不含纹理信息的室内三维点云和彩色纹理影像，后期通过人工识别回光反射标志的方式获取可靠的同名点。

本文采用Visual Studio 2015以及C++语言进行程序编写，使用的操作系统为Windows 8.1 64位专业版，CPU为8核Intel(R)Core(TM)i7-4790。

2.2 单像位姿估计及多像位姿解算

本文首先利用4组真实数据进行单张影像位姿估计实验，然后将实验结果应用到多像位姿解算中。采用的实验数据如表 1所示。

本文主要从数值稳定性(平均反投影误差)、运算时间、实际贴图效果这3个方面对本文方法、UPnP^[25]、DLT^[16-17]、Gao’s P3P^[23]进行对比。为了能够更加精确地对比各种方法的运算效率，本文对不同方法的运算时间迭代50次并记录运算时间，所得到的各个方法的实际运算时间如图 5所示。

为了验证本文位姿估计算法的精度，选取24对控制点，前12个点求解影像位姿，后12个点作为检核点，将其三维坐标代入式（8）中反求其在影像上的像素坐标，并与在影像上选取的对应点像素坐标作差得到(Δx_i，Δy_i)。检核点的平均反投影误差r的计算公式为：

$ r = \frac{{\sum\limits_{i = 1}^n {\sqrt {\Delta x_i^2 + \Delta y_i^2} } }}{n},n = 12 $

(41)

反投影误差对比直方图如图 6所示。

由于控制点均由工作人员进行选取，故各个控制点的系统误差以及分布状况大致相同，所包含的噪声也被限制在了一定范围内。从图 5中可以看出，DLT运算时间较本文方法与UPnP长，这主要是因为其需要进行迭代计算，同时需要提供初始值，初始值的差异决定了迭代次数；P3P方法由于其数学模型的限制结果必然为多解，对结果进行筛选的过程也必然会降低计算效率，P3P解法和DLT解法常与RANSAC或反向随机抽样一致(a contrario random sample consensus, AC-RANSAC)结合进行最优解的筛选^[27]。而本文方法与UPnP均大量使用矩阵分解来计算摄像机位置姿态，算法复杂度较低。同时从图 5中可以看出本文方法较UPnP在算法效率上稍胜一筹。另外，在位姿估计精度方面，由图 6可以看出DLT反投影误差较大，这主要是因为DLT解法需要提供初始值进行迭代，若初始值与真实值差距过大，则在过少的迭代次数内可能无法得到精确值；而P3P解法只选择3个控制点进行解算，这种方法最多有4个解，需要另外加入至少一个控制点进行筛选，且3个控制点易分布在同一平面上，导致无法求解；UPnP解法在时间效率上与本文方法相当，但精度却稍逊一筹，这主要是由于其解算过程中较本文方法多了一个未知量——焦距f。

利用本文位姿估计方法，分别对后母戊鼎(HMWD)、太和门(THM)、北京建筑大学小鸟巢(Nest)、北京建筑大学至善亭(Pavilion)进行实验验证，均能得到摄像机在世界坐标系中正确的位置姿态，进而对其进行纹理映射验证其结果的正确性，如图 7所示。

为了验证§1.5多像外方位元素解算方法以及位姿旋转方法的可行性，利用安置的数码相机进行全景拍摄，总共摄取10张影像，则每张影像间的旋转角为β = 360°/10 = 36°。通过人工识别分别在点云与对应的初始影像上选取控制点，利用§1所述方法进行单像位姿估计得到初始影像外参数，而后利用多像位姿解算得到其余影像的外参数。但是后续影像的标定精度依赖于首张影像的标定精度，所以本文采用光束法平差^[31]的方法对后续影像的位姿参数进行优化。由于相机在旋转拍摄过程中，每两张影像之间均存在重叠区域，则可以通过上述方法获取的初始影像姿态参数得到每个三维点在哪两张影像上可以观测到，并记录其在影像上的像素坐标，这样在点云中便能够得到“二度重叠点”。最后，将“二度重叠点”作为真值、相机内外参数作为待优化的参数进行整体光束法平差^[29]优化。为了检验优化后的精度，本文在后续影像与点云上分别选取了4个控制点作为检核点，仍然采用式（8）的方法计算其反投影误差。选取的4个控制点如表 2所示。

使用平差前与平差后的相机内外参数计算每个点的反投影误差，如表 3所示。

表 3中，∆u与∆V为平差前的影像横纵坐标反投影误差，∆u'与∆V'为平差后的横纵坐标反投影误差。由表 3中数据可知，随着相机旋转拍摄，相机姿态误差越来越大。另外，经过光束法平差后，后续影像的纹理映射精度能够得到有效的提高。但是，经过平差后仍然存在一定的像素误差，这可能是由于在人工选取控制点时存在偶然误差，并且本文采用的影像像素分辨率较大，为2 000万左右，增加了精确选取影像控制点的难度。

最后利用摄像机透视投影模型对点云进行纹理映射，得到全景点云纹理映射结果。图 8所示为纹理映射前后全景点云效果对比。

纹理映射前后点云内部视角效果如图 9所示。

2.3 匀光匀色处理

由于数码相机在室内旋转拍摄的过程中光照条件变化剧烈，相邻影像产生辐射差异，导致纹理融合的过程中产生明显接缝，图 8与图 9中红色框即说明此类情况。因此，本文选取了目前最流行的3种匀光方法——Mask匀光、Wallis滤波、直方图规定化进行实验对比，发现对于本文所使用的数码影像，采用Mask匀光法能够很好地增强多张影像辐射值的整体一致性，而Wallis滤波在增强影像的对比度的同时抑制了影像的噪声，产生了多余的纹理以及较强的分块效应；直方图规定化选择第一张影像作为基础影像，其余影像的直方图与基础影像的直方图进行匹配，但由于基础影像的色调偏红，导致其余影像的红色通道所占比重增强，容易丧失原始色调。

经过综合考虑，本文采用Mask匀光方法先对影像进行整体辐射矫正，然后用§1.6所述的线性融合方法进行接缝处理，得到的结果如图 10所示，与之前对比，可以明显地看出接缝已经消除。

3 结语

针对三维激光点云与数码影像的多源数据融合问题，本文采用了三维激光扫描仪外置数码相机的方法，使数码相机随扫描仪水平旋转一周拍摄多张影像，通过基于虚拟控制点的摄像机位姿估计方法计算初始影像的位置姿态，从而获得其余影像相对世界坐标系的姿态，完成全景三维点云与多张数码影像的多源数据的融合，得到了具有丰富纹理信息的全景彩色点云。本文算法在效率与数值稳定性方面均优于P3P、UPnP、DLT等解法，并且对大倾角拍摄的状况稳定性良好，在这一点上优于传统的摄影测量解法。