激光雷达和相机均是实现快速实景复刻的测量设备. 激光雷达基于脉冲或相位式原理测量周围环境到测站的距离，并结合角度码盘生成三维点云，优点是数据精度高、抗干扰性强，缺点是缺乏色彩和纹理信息，不利于解读复杂场景 ^[1]. 全景图像由多台相机同步采集的图像拼接而成，具有360°环境下丰富的视觉信息，被广泛应用于目标检测、场景理解、三维重建等领域，但缺乏空间三维信息^[2]. 地面三维扫描仪惯以将激光雷达和相机同轴测量，并将点云与全景影像配准融合得到富有视觉特征的彩色点云，被应用于目标识别、三维建模、智能测绘等技术领域^[3-4]. 全景图像与点云融合的关键是构建配准方法，即建立数据间的映射关系，主要分为以下4个类别.

1）基于2D标靶. 建立激光雷达与相机外参数模型，利用特征点线解算模型参数，标靶制作容易且标定精度高，但标定过程复杂且自动化程度低. Zhang等^[5]利用棋盘格点在面上的几何约束完成二维激光雷达与相机外参数标定，但复杂环境的棋盘格点云提取精度受限. Geiger等^[6]在同一场景中放置多个棋盘格，提出基于生长点检测法实现了亚像素级测量，但标定结果受靶面位姿的影响. 镂空棋盘、ArUco标签、三角孔等2D标靶的特征提取精度均未摆脱点云密度的限制^[7-9].

2）基于3D标靶. 利用3D标靶标定激光雷达与相机的外参数能够解决2D标靶标定传感器外参数丢失Z轴约束的问题. 以特征点线为基元，标定结果精度高且稳定，但标定物制作困难，特征提取需要人工干预. Dong等^[10]制作非共面V形标定板，提出最小化点到平面距离法来获取外参数. 该方法对标定物的制作精度要求高，标定工作量和难度较大. Fan等^[11]设计球和空心标靶结合的标定物，利用空心标定板边缘线特征计算转换模型初始值，利用球心迭代优化参数，实现特征自动提取与配对，但该方法过于复杂，计算效率较低. 球面、格网、对折板等3D标靶的标定结果精度高^[12-14]，但都有特定的应用场景，通用性差.

3）基于无标靶. 利用自然场景的点线信息或灰度信息完成激光雷达与相机外参数标定，分为基于特征信息和基于互信息最大化2类标定方法. Bai等^[15]利用激光点云与图像中平行线上的点求解旋转矩阵，根据点在线上的约束得到平移向量，该方法要求应用场景中有平行线特征. Pandey等^[16]利用点云回波反射率与图像灰度值相关系数，使二者间联合直方图分散最小获得外参标定，在全景相机和激光系统中分别进行验证，该方法标定精度易受环境光照影响. 范光宇等^[17]提出基于灰度相似性的配准方法，依据全景影像与深度图在水平、垂直方向细分区域的灰度相似性完成配准. 该方法配准精度与分割密度相关，且随着分割密度增加计算量剧增，冗余计算过多且数据配准有限. 将车道线、天际线、楼角等作为基元的方法^[18-19]虽然灵活、简单，但标定精度普遍较低且稳定性差，依赖于特定环境，特征配对易出错.

4）基于深度学习. 这类算法通过设计目标函数评价点云与图像配准效果，利用优化方法更新外参数，逐步迭代改善匹配效果，取最优评价作为标定结果. Schneider等^[20]将深度学习算法用于激光雷达与相机外参标定问题，构建RegNet神经网络分别提取相机图像和激光雷达点云，利用回归计算完成标定. Yuan等^[21]提出基于RGGNet的激光雷达与相机外参标定方法，利用深度信息生成模型学习隐式容差模型，同时考虑标定误差和误差范围内容差，标定效果较好. 基于深度学习的标定还有CMRNet、CalibRCNN ^[22]. 该类方法多为无/半监督学习下的深度学习框架，自动化程度高，但泛化性有待提高.

上述4类配准方法均须建立相机与激光雷达的外参标定模型，利用点线或灰度特征解算模型参数，从而构建数据的映射关系，标定过程复杂，模型参数较多，不便于设备参数定期检校^[22]. 本研究提出直接构建扫描仪三维点云与图像间像素级映射关系的方法；将具有旋转不变性的球形标靶作为特征基元，分别利用色彩信息和几何信息从全景图像和点云中高精度提取基元；在空间球坐标系下构建基元三角形，利用最小角距差法完成基元配对，并建立点云与全景影像的初始映射关系，基于改进Levenberg-Marquardt（L-M)算法和自由形变法(free-form deformation，FFD)的混合算法（iLM-FFD）优化映射关系.

如图1所示为点云与全景图像映射关系的算法框架，其中B₁B₂B₃为基元三角形. 将三维点云与全景图像作为数据输入，通过1）在图像中提取特征基元，2）在点云中提取特征基元，3）基元配对与建立初始映射，4）混合算法优化映射关系4个主要步骤，实现在无需传感器外参数标定的情况下建立点云与全景图像间像素级的映射关系. 全景图像经预处理后，检测圆形特征感兴趣区域（region of interest，ROI），利用一维最大熵的改进Zernike矩亚像素边缘提取算法提取标靶中心，作为图像中的基元. 获取点云强度图，检测候选圆球ROI，利用球面点的三维信息从点云中提取可信靶球中心，作为点云中的基元. 利用最小角距差法将点云与图像的基元配对，并建立点云与全景图像的初始映射关系. 运用混合算法优化畸变导致的映射偏差，提升数据质量.

球形标靶具有旋转不变性，中心拟合精度高，颜色一般比周围环境反差大，为此通过检测圆形特征ROI和提取亚像素边缘来从全景图像中自动提取标靶中心. 全景图像信息单元较多全局处理效果欠佳，为此将图像分割成200张等尺寸子图，相邻图像间重叠100 pixel，检测每张子图的ROI. 每张子图像完成预处理、二值化、连通域处理等步骤，获得多个灰度区域；建立第$ i $个灰度区域的外接矩形包围框（宽$ {w_i} $、高$ {h_i} $、面积$ {S_i} $），以$ | {w_i} - {h_i} | \geqslant 2 $和$ {S_1} > {S_i} > {S_2} $为约束条件判断候选区域，$ {S_1}、{S_2} $由焦距和标靶距离估算；圆与外接矩形面积理论比值约为0.78，靶球所在像素的灰度值为1，因此当矩形框中像素比率$ r > 0.75 $时，将该矩形作为ROI.

(1)$ r = \frac{{{N_{\mathrm{g}}}}}{{{N_0}}}. $

式中：$ {N_{\mathrm{g}}} $为包围框内灰度值为1的像素数量，$ {N_0} $为包围框内所有像素数量. 提取ROI中的标靶边缘，获取图像中的特征基元. Zernike矩边缘定位精度较高、耗时少，具有良好的抗噪性能，但用以判定边缘的阶跃阈值$ {k_{\mathrm{t}}} $难以确定^[23]. 一维最大熵法是获取最大灰度信息熵的算法，能够自适应选取合理的边缘灰度^[24]. 本研究提出基于一维最大熵的改进Zernike矩亚像素检测算法，自适应提供$ {k_{\mathrm{t}}} $以判断标靶边缘. Teague^[25]根据Zernike正交实值多项式定义图像$ f(x,y) $的$ n $阶$ m $次Zernike矩为

(2)$ {Z_{nm}} = \frac{{n+1}}{{\text{π }}}\iint\limits_{{x^2}+{y^2} \leqslant 1} {f(x,y)V_{nm}^ * {{{\mathrm{d}}x{\mathrm{d}}y}}}. $

式中：$ V_{nm}^ * $为Zernike多项式$ {V_{nm}} $的共轭. 依据Zernike矩的旋转不变性实现亚像素边缘定位. 如图2所示为理想边缘模型. 阴影部分旋转$ \theta $前后的Zernike矩只改变相角，仍保持幅值不变，即$ {Z'_{nm}} = {Z_{nm}}{{\mathrm{e}}^{ -{\mathrm{ j}}m\theta }} $. 直线$ L $为理想边缘，$ L $两侧的灰度值分别为$ h $和$ h+k $，$ k $为阶跃灰度值，$ l $为原点到边缘的垂直距离，$ \theta $为$ l $和$ x $轴的夹角. 确定边缘需要通过$ {Z_{00}}、{Z_{11}}、{Z_{20}} $3个不同阶次的矩进行计算，由旋转前后Zernike的幅值不变性得到$ {Z'_{00}} = {Z_{00}} $、$ {Z'_{11}} = {Z_{11}}{{\mathrm{e}}^{{\mathrm{j}}\theta }} $、$ {Z'_{20}} = {Z_{20}} $，则旋转后各阶Zernike矩为

(3)$ \left. {\begin{array}{*{20}{l}} {{{Z'}_{00}} = h{{\text{π}}} +{{k{\text{π}} } / 2} - k\arcsin \;l - kl\sqrt {1 - {l^2}} ,} \\ {{{Z'}_{11}} = {{2k{{(1 - {l^2})}^{{3}/{2}}}} / 3},} \\ {{{Z'}_{20}} = {{2kl{{(1 - {l^2})}^{{3}/{2}}}} / 3}.} \end{array}} \right\} $

(4)$ \left. {\begin{array}{*{20}{l}} {\theta = \arctan \left(\dfrac{{{{\mathrm{Im}}}\,\, {Z_{11}}}}{{{{\mathrm{Re}}}\,\, {Z_{11}}}}\right),} \\ {\begin{array}{*{20}{l}} {d = \dfrac{{{Z_{20}}}}{{{{Z'}_{11}}}},}&{k = \dfrac{{3{{Z'}_{11}}}}{{2{{(1 - {d^2})}^{3/2}}}},} \end{array}} \\ {h = {{({Z_{00}} - \dfrac{{k{\text{π}} }}{2}+k\arcsin d+kd\sqrt {1 - {d^2}}) }}/{{\text{π}} }.} \end{array}} \right\} $

式中：$ {{\mathrm{Re}}}\,\, {Z_{11}} $和$ {{\mathrm{Im}}}\,\, {Z_{11}} $分别为$ {Z_{11}} $的实部和虚部. 满足$ d \leqslant {d_{\mathrm{t}}} $、$ k \geqslant {k_{\mathrm{t}}} $的像素点判定为亚像素边缘，$ {d_{\mathrm{t}}} $为小于1 Pixel的距离阈值，取$ {d_{\mathrm{t}}} $=0.5 pixel；$ {k_{\mathrm{t}}} $常通过人工选择，过小会检测出虚假边缘，过大会丢失边缘信息，须反复调试. 一维最大熵法可自动选取$ {k_{\mathrm{t}}} $，提高检测效率. 假设初始阈值为$ T $，$ {p_0}, {p_1},\cdots,{p_n} $为图像中各灰度等级的占比. 前景灰度A的占比为$ {p_{\mathrm{A}}} = {p_0}+{p_1}+ \cdots +{p_T} $，背景灰度B的占比为$ {p_{\mathrm{B}}} = {p_{T+1}}+ {p_{T+2}}+\cdots + {p_n} $，且$p_{\mathrm{A}}+p_{\mathrm{B}}=1 $. 遍历整个灰度层级，得到使$ H(T) $最大的$ T $，令$ {k_{\mathrm{t}}} = T $；

(5)$ H(T) = - \sum\limits_{i = 0}^T {\frac{{{p_i}}}{{{p_{\mathrm{A}}}}}} \lg \left(\frac{{{p_i}}}{{{p_{\mathrm{A}}}}}\right) - \sum\limits_{i = T+1}^n {\frac{{{p_i}}}{{{p_{\mathrm{B}}}}}} \lg \left(\frac{{{p_i}}}{{{p_{\mathrm{B}}}}}\right). $

采用$ N $×$ N $的矩阵模版卷积，积分圆区域半径为$ {N \mathord{\left/ {\vphantom {N 2}} \right. } 2} $，亚像素边缘坐标为

(6)$ \left[ {\begin{array}{*{20}{c}} {{x_{\mathrm{s}}}} \\ {{y_{\mathrm{s}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{x_0}} \\ {{y_0}} \end{array}} \right]+\frac{N}{2}d\left[ {\begin{array}{*{20}{c}} {\cos \theta } \\ {\sin \theta } \end{array}} \right]. $

利用边缘点拟合圆，剔除残差较大的点，得到最佳圆心坐标$ ({x_i},{y_i}) $. 依据球面成像模型将$ ({x_i},{y_i}) $转换为水平角、高度角$ \left\{ {\left( {{A_{{\mathrm{t}}i}},{H_{{\mathrm{t}}i}}} \right)} \right\}_{i = 1}^n $，作为图像中的基元，$ n $为图像中基元的数量.

利用球面点的三维信息从点云中高精度提取标靶中心. 以站心为原点将点云投影成强度图，将投影到同个像素的点云定义为点云群. 从强度图中获取候选ROI，利用ROI中的点云群拟合球面，得到拟合球半径$ {r_0} $与球心$ ({x_0},{y_0},{z_0}) $. 已知标靶中心偏移误差约为1~2 mm，将$ 4.8\;{\mathrm {cm}} < {r_0} < 5.2\;{\mathrm{ cm}} $作为约束条件. 若点云到球心的距离$ d(O,{p_i}) $与$ {r_0} $之差小于$ 2\sigma $ ($ \sigma $为设备三维点位误差，由仪器厂商标定)，该点可能是标靶上的点，进而将圆球度$ C > $90%的候选ROI作为可信标靶（此处90%为经验阈值）.

(7)$ C = \frac{{{M_1}}}{{{M_0}}} \times 100{\text{%}} . $

(8)$ \left. {\begin{array}{*{20}{l}} {{d_i} = \left| {d(O,{p_i}) - {r_0}} \right|,\,\, i = 1,2,\cdots,{M_0};} \\ {{e_{\mathrm{r}}} = \sqrt {{{\displaystyle\sum\limits_{i = 0}^N {{d_i}{d_i}} }}/{{{M_0}}}} .} \end{array}} \right\} $

式中：$ {M_1} $为$ \left| {d(O,{p_k}) - {r_0}} \right| < 2\sigma $的点数，$ {M_0} $为拟合总点数. 利用可信标靶的点云拟合球体，计算拟合残差$ {d_i} $和残差中误差$ {e_{\mathrm{r}}} $，剔除$ {d_i} > 2{e_{\mathrm{r}}} $的点，重复拟合球体，将最终拟合的球心坐标转换为方位角$ \left\{ {\left( {{A_{{\mathrm{p}}j}},{H_{{\mathrm{p}}j}}} \right)} \right\}_{j = 1}^m $作为点云中的基元，$ m $为点云中基元的数量.

准确配对点云与图像中的特征基元是建立数据映射关系的关键. 智能扫描仪的激光器中心与相机中心在结构上常同轴作业，杨必胜等^[4]研究发现二者数据贴合度较高. 如图3(a)所示为全景图像，圆点为靶球的位置，直线段表示图像空间球坐标系中靶球之间的距离；图3(b)为点云强度图，圆点为靶球位置，直线段表示点云空间球坐标系中靶球之间的距离. 由图可知，靶球基元在点云和图像空间球坐标系中的几何结构相似. 本研究提出最小角距差法以完成基元配对. 从$ \left\{ {\left( {{A_{{\mathrm{p}}j}},{H_{{\mathrm{p}}j}}} \right)} \right\}_{j = 1}^m $中选取3个基元组成球面基元三角形，计算点云基元间角距$ {G_{\mathrm{p}}} = $($ {g_{31}} $,$ {g_{12}} $,$ {g_{23}} $). 同理，从$ \left\{ {\left( {{A_{{\mathrm{t}}i}},{H_{{\mathrm{t}}i}}} \right)} \right\}_{i = 1}^n $中取3个基元，得图像基元间角距$ {G_{\mathrm{t}}} = $($ {g'_{31}} $,$ {g'_{12}} $,$ {g'_{23}} $). 将角距差平均值$ \Delta g $作为评价点云基元三角形与图像基元三角形相似度的依据，即将基元配对问题转化为最小角距差查找问题.

(9)$ \left. {\begin{array}{*{20}{l}} \begin{gathered} {g_{ij}} = \arccos\,\, (\cos {H_{{\mathrm{p}}i}}\cos {A_{{\mathrm{p}}i}}\cos {H_{{\mathrm{p}}j}}\cos {A_{{\mathrm{p}}j}}+ \\ \begin{array}{*{20}{c}} {}&{\cos {H_{{\mathrm{p}}i}}\sin {A_{{\mathrm{p}}i}}\cos {H_{{\mathrm{p}}j}}\sin {A_{{\mathrm{p}}j}}+\sin {H_{{\mathrm{p}}i}}\sin {H_{{\mathrm{p}}j}}),} \end{array} \\ \end{gathered} \\ \begin{gathered} {g'_{ij}} = \arccos\,\, (\cos {H_{{\mathrm{t}}i}}\cos {A_{{\mathrm{t}}i}}\cos {H_{{\mathrm{t}}j}}\cos {A_{{\mathrm{t}}j}}+ \\ \begin{array}{*{20}{c}} {}&{\cos {H_{{\mathrm{t}}i}}\sin {A_{{\mathrm{t}}i}}\cos {H_{{\mathrm{t}}j}}\sin {A_{{\mathrm{t}}j}}+\sin {H_{{\mathrm{t}}i}}\sin {H_{{\mathrm{t}}j}}).} \end{array} \\ \end{gathered} \end{array}} \right\} $

(10)$ \Delta g = {\mathrm{mean}}\Bigg(\sum\limits_{i = 1} {\sum\limits_{j = 1} {| {{g_{ij}} - {g'_{ij}}} |} } \Bigg) < \tau . $

式中：$ {g_{ij}} $、$ {g'_{ij}} $为第$ i $和$ j $个基元分别在点云空间球坐标系和图像空间球坐标系中的角距，$ \tau $为经验阈值，设$ \tau $=200". 遍历$ \left\{ {\left( {{A_{{\mathrm{t}}i}},{H_{{\mathrm{t}}i}}} \right)} \right\}_{i = 1}^n $中所有三角形组合，若有多组满足阈值的配对结果，将$ \Delta g $取最小时对应的三角形组合作为最终配对结果. 由于激光器与相机同轴作业，二者中心的距离相比于测量距离较小，在空间球面坐标系下可近似认为激光器与相机数据同心^[22]. 因此，点云空间方位角$ ({A_{\mathrm{p}}},{H_{\mathrm{p}}}) $与图像空间方位角$ ({A_{\mathrm{t}}},{H_{\mathrm{t}}}) $之间的初始映射关系可表示为

(11)$ {A_{\mathrm{p}}} = {A_{\mathrm{t}}}+{\mathrm{d}}A,\; {H_{\mathrm{p}}} = {H_{\mathrm{t}}}+{\mathrm{d}}H. $

式中：$ {\mathrm{d}}A,{\mathrm{d}}H $为映射参数，由基元对算得. 建立点云$ ({A_{\mathrm{e}}},{H_{\mathrm{e}}}) $到图像$ ({u_{\mathrm{e}}},{v_{\mathrm{e}}}) $的初映射关系，将图像RGB传递到点云，得到初始彩色点云，转换式如下，其中int为取整符号：

(12)$ \begin{split} {\mathrm{Con}}\;&(({A_{\mathrm{e}}},{H_{\mathrm{e}}}),({u_{\mathrm{e}}},{v_{\mathrm{e}}})) = \\& \left\{ {\begin{array}{*{20}{l}} {{v_{\mathrm{e}}} = {{\mathrm{int}}}\;((90 - ({H_{\mathrm{e}}} - {\mathrm{d}}H))/{{0.05}^ \circ });} \\ \begin{array}{*{20}{l}} {{u_{\mathrm{e}}} ={{\mathrm{int}}}\;(({A_{\mathrm{e}}} - {\mathrm{d}}A+{{180}^ \circ })/{{0.05}^ \circ }),} \\&{{\mathrm{d}}A < {0^ \circ },\; - {{180}^ \circ } < {A_{\mathrm{e}}} < {{180}^ \circ }+{\mathrm{d}}A;} \end{array} \\ {\begin{array}{*{20}{l}} {{u_{\mathrm{e}}} = {{\mathrm{int}}}\;(({A_{\mathrm{e}}} - {\mathrm{d}}A - {{180}^ \circ })/{{0.05}^ \circ }),} \\&{{\mathrm{d}}A < {0^ \circ },\;{{180}^ \circ } > {A_{\mathrm{e}}} \geqslant {{180}^ \circ }+{\mathrm{d}}A;} \end{array}} \\ {\begin{array}{*{20}{l}} {{u_{\mathrm{e}}} = {{\mathrm{int}}}\;(({A_{\mathrm{e}}} - {\mathrm{d}}A+{{180}^ \circ })/{{0.05}^ \circ }),} \\&{{\mathrm{d}}A > {0^ \circ },\; - {{180}^ \circ }+{\mathrm{d}}A \leqslant {A_{\mathrm{e}}} < {{180}^ \circ };} \end{array}} \\ {\begin{array}{*{20}{l}} {{u_{\mathrm{e}}} = {{\mathrm{int}}}\;(({A_{\mathrm{e}}} - {\mathrm{d}}A+{{540}^ \circ })/{{0.05}^ \circ }),} \\&{{\mathrm{d}}A > {0^ \circ },\; - {{180}^ \circ } < {A_{\mathrm{e}}} < - {{180}^ \circ }+{\mathrm{d}}A.} \end{array}} \end{array}} \right.\end{split} $

由于区域性图像畸变、错位造成部分区域的初映射关系存在一定偏差，达到像素级映射精度须进一步优化映射关系. 自由变形是局部特性变换模型，通过对变形图($ m \times n $)布设$ {d_{\mathrm{r}}} \times {d_{\mathrm{c}}} $的网格，以网格节点为控制点，每个控制点对应$ x $方向和$ y $方向2个控制参数$ {\varphi _x},{\varphi _y} $，使用变形图像每个像素点周围4×4个控制点的控制参数计算位置偏移^[26]. 在初映射关系的基础上，提出基于改进L-M与FFD的混合优化算法(iLM-FFD)，将优化映射关系的问题转换为以原始点云$ P(x,y,z) $强度图校准赋色后点云$ P'(x,y,z) $灰度图的问题. 以强度图为基准图$ {P_{{\mathrm{ref}}}} $，以赋色后点云灰度图为浮动图$ {P_{{\mathrm{def}}}} $，使用FFD作为变换模型对$ {P_{{\mathrm{def}}}} $逐像素变形，并与$ {P_{{\mathrm{ref}}}} $配准，并计算$ {P_{{\mathrm{ref}}}} $与变形后的$ {P_{{\mathrm{def}}}} $的相似度，解算使相似度最大的最优控制参数，同时基于改进的L-M算法加速参数优化过程，利用最优控制参数对$ {P_{{\mathrm{def}}}} $变形，进而优化映射关系.

归一化互相关系数(normalized cross-correlation , NCC)不受灰度值线性变化影响，适用于纹理、模糊图像，将该参数作为对变形后$ {P_{{\mathrm{def}}}} $与$ {P_{{\mathrm{ref}}}} $相似性的评价指标，以分块取平均的方式计算全景图像的归一化互相关系数:

(13)$ {{\mathrm{NCC}}} = \frac{1}{w}\displaystyle\sum\limits_{k = 1}^w {{\mathrm{NC}}{{\mathrm{C}}_k}\left(\frac{{\displaystyle\sum\limits_{i = 1}^m {\displaystyle\sum\limits_{j = 1}^n {({P_{{\mathrm{ref}}}}(x,y) - {{\overline P}_{{\mathrm{ref}}}})({P_{{\mathrm{def}}}}(x,y) - {{\overline P}_{{\mathrm{ref}}}})} } }}{{\sqrt {\displaystyle\sum\limits_{i = 1}^m {\displaystyle\sum\limits_{j = 1}^n {{{({P_{{\mathrm{ref}}}}(x,y) - {{\overline P}_{{\mathrm{ref}}}})}^2}} } } \sqrt {\displaystyle\sum\limits_{i = 1}^m {\displaystyle\sum\limits_{j = 1}^n {{{({P_{{\mathrm{def}}}}(x,y) - {{\overline P}_{{\mathrm{def}}}})}^2}} } } }}\right)} . $

式中，$ {\overline P_{{\mathrm{ref}}}}、{\overline P_{{\mathrm{def}}}} $为$ {P_{{\mathrm{ref}}}}、{P_{{\mathrm{def}}}} $块的灰度均值，$ m \times n $为像素总数，$ w $为图像分块数，$ {{\mathrm{NCC}}} \in [0,1] $. 由于L-M算法为求解目标的最小值，令模型的目标函数为$ F = 1 - {{\mathrm{NCC}}} $. 设置网格的尺寸为$ (r+3) $×$ (c+3) $，生成长度为$ 2(r+3) \times (c+3) = N $的一维数组由[−0.001,0.001]的随机数赋初始值. 基于FFD计算$ {{\mathrm{NCC}}} $，先计算像素点的网格坐标整数$ ({G_x},{G_y}) $和小数$ (u,v) $，利用B样条基函数$ F_{{\mathrm{BB}}}(s,l) $计算x方向和y方向权重系数$ {P_{x0}} $、$ {P_{x1}} $、$ {P_{x2}} $、$ {P_{x3}} $、$ {P_{y0}} $、$ {P_{y1}} $、$ {P_{y2}} $、$ {P_{y3}} $，再计算坐标偏移值并插值.

(14)$ \left. {\begin{array}{*{20}{c}} {x' = x+{\Delta} x = \displaystyle\sum\limits_{p = 0}^3 {\displaystyle\sum\limits_{q = 0}^3 {{P_{xq}}{P_{yp}}{\varphi _x}({G_x}+p,{G_y}+q)} } ,} \\ {y' = y+{\Delta} y = \displaystyle\sum\limits_{p = 0}^3 {\displaystyle\sum\limits_{q = 0}^3 {{P_{xq}}{P_{yp}}{\varphi _y}({G_x}+p,{G_y}+q)} } .} \end{array}} \right\} $

$ (x',y') $为对应基准图像的浮点坐标，须由双线性插值计算灰度值并赋值给$ {P_{{\mathrm{def}}}}(x,y) $.

(15)$ \left. {\begin{array}{*{20}{l}} {{x_1} = {\mathrm{floor}}\;(x'),\; {x_2} = {x_1}+1;} \\ {{y_1} = {\mathrm{floor}}\;(y'),\;\; {y_2} = {y_1}+1;} \\ \begin{gathered} {P_{{\mathrm{def}}}}(x,y) = ({x_2} - x')({y_2} - y'){P_{{\mathrm{ref}}}}({x_1},{y_1})+({x_1} - x'){\text{×}} \\ \begin{array}{*{20}{c}} {}&{({y_1} - y'){P_{{\mathrm{ref}}}}({x_2},{y_2}) - ({x_1} - x')({y_2} - y'){P_{{\mathrm{ref}}}}{\text{×}} } \end{array} \\ \end{gathered} \\ {\begin{array}{*{20}{c}} {}&{({x_2},{y_1}) - ({x_2} - x')({y_1} - y'){P_{{\mathrm{ref}}}}({x_1},{y_2}).} \end{array}} \end{array}} \right\} $

式中：floor为向下取整符号. 基于改进的L-M算法迭代优化控制参数. 差分法计算控制参数的梯度，$ \Delta d = 1 $.

(16)$ \begin{split} \frac{{\partial F}}{{\partial {X_i}}} = &\frac{{F({X_1},{X_2},\cdots,{X_i}+\Delta d,\cdots,{X_N})}}{{\Delta d}} - \\&{\frac{{F({X_1},{X_2},\cdots,{X_i},\cdots,{X_N})}}{{\Delta d}}.} \end{split} $

所有控制点的梯度组成一维梯度向量$ {{\mathbf{\nabla\varphi }}} = [{{\partial F}}/{{\partial {X_1}}}, {{\partial F}}/{{\partial {X_2}}},\cdots,{{\partial F}}/{{\partial {X_N}}}] $，也是改进L-M算法的雅克比矩阵$ {\boldsymbol{J}} $；计算当前控制参数的目标函数$ {f_v} = {{\boldsymbol{F}}_v}({X_1}, {X_2},\cdots,{X_N}) $和相关矩阵$ {\boldsymbol{h}} = - {({{\boldsymbol{J}}^{\text{T}}}{\boldsymbol{J}}+{\lambda _v}{\boldsymbol{I}})^{ - 1}}{{\boldsymbol{J}}^{\text{T}}}{f_v} $，$ {\boldsymbol{I}} $为的单位阵，$ {\lambda _v} $为控制步长. 将满足$ \left\| {{\boldsymbol{J}}_v^{\mathrm{T}}{{\boldsymbol{F}}_v}} \right\| \leqslant \varepsilon $解的$ {{\mathbf{\varphi }}_v}({X_1},{X_2},\cdots,{X_N}) $作为最优控制参数，$ \varepsilon $取很小的数. 根据计算结果，调整$ {\lambda _v} $，更新控制参数，计算目标函数值和步长因子$ \rho $.

(17)$ \left.\begin{array}{l}{{\boldsymbol{X}}}^{\prime }={{\boldsymbol{X}}}_{v}+{{\boldsymbol{h}}}^{{\mathrm{T}}}\text{，}{f}^{\prime }=F({{\boldsymbol{X}}}^{\prime }); \\\rho =\dfrac{{f}_{v}-{f}^{\prime }}{0.5{{\boldsymbol{h}}}^{{\mathrm{T}}}({\lambda }_{v}{\boldsymbol{h}}-{{\boldsymbol{J}}}^{{\mathrm{T}}}{f}_{v})}. \end{array}\right\} $

合理更新控制步长$ {\lambda _v} $有利于提高优化效率. Yamashita等^[27]证实$ {\lambda _v} = {\left\| {{{\boldsymbol{F}}_v}} \right\|^2} $时，在一定误差范围内L-M算法具有二次收敛性；Fan^[28]提出$ {\lambda _v} = {\mu _v}\left\| {{{\boldsymbol{F}}_v}} \right\| $，同样满足类似的收敛性，可逆转初始值远离解时收敛效果差的问题；Amini等^[29]提出$ {\lambda _v} = {{{\mu _v}\left\| {{{\boldsymbol{F}}_v}} \right\|} / {(1+\left\| {{{\boldsymbol{F}}_v}} \right\|)}} $，证明在解的一定误差范围内具有全局收敛和二次收敛性. 设控制步长为

(18)$ {\lambda _v} = \frac{{{\mu _v}{{\left\| {{{\boldsymbol{F}}_v}} \right\|}^2}}}{{1+{{\left\| {{{\boldsymbol{F}}_v}} \right\|}^2}}}. $

当结果远离解时，$ {\lambda _v} $接近标准L-M算法中的$ {\mu _v} $；当结果接近解时，$ {\lambda _v} $接近$ {\mu _v}{\left\| {{{\boldsymbol{F}}_v}} \right\|^2} $，与文献[27]类似，且$ {\left\| {{{\boldsymbol{F}}_v}} \right\|^2} $可加速迭代. 为了提高算法参数的更新效率，利用下降比率$ {r_v} $作为判断依据，每次计算的相似度与最差相似度比较：

(19)$ {r_v} = \frac{{\mathop {\max }\limits \left\{ {{{\left\| {{{\boldsymbol{F}}_j}} \right\|}^2}} \right\} - {{\left\| {{\boldsymbol{F}}({{\boldsymbol{X}}_v}+{\boldsymbol{h}})} \right\|}^2}}}{{{{\left\| {{{\boldsymbol{F}}_v}} \right\|}^2} - {{\left\| {{{\boldsymbol{F}}_v}+{\boldsymbol{Jh}}} \right\|}^2}}},\;j=1,\cdots, v. $

(20)$ \left. \begin{gathered} {{\boldsymbol{X}}_{v+1}} = \left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol{X}}_v}+{\boldsymbol{h}}},&{{r_v} \geqslant {p_0};} \\ {{{\boldsymbol{X}}_v}},&{{r_v} < {p_0}.} \end{array}} \right. \\ {\mu _{v+1}} = \left\{ {\begin{array}{*{20}{l}} {3{\mu _v}},&{{r_v} < {p_1};} \\ {{\mu _v}},&{{r_v} \in [{p_1},{p_2}];} \\ {\max \left\{ {{{\mu _v}}}/{3},{\mu _v}\eta \right\}} ,&{{r_v} > {p_2}.} \end{array}} \right. \\ \end{gathered} \right\} $

式中：$ \eta = 1 - {(2\rho - 1)^3} $，$ 0 < {p_0} \leqslant {p_1} \leqslant {p_2} < 1 $. 由式(20)更新参数. 进入下一步迭代计算. 得到最优控制参数$ {\mathbf{\varphi }({\boldsymbol{X}})} $，建立优化后的$ {P_{{\mathrm{def}}}} $与$ {P_{{\mathrm{ref}}}} $的映射关系$ {P_{{\mathrm{ref}}}} = {\mathrm{Con}}\;({P_{{\mathrm{def}}}},{\mathbf{\varphi }({\boldsymbol{X}})}) $. 以初始映射点云$ P'(x,y,z) $为中介，建立点云$ P(x,y,z) $与像素$ T(u,v) $间的映射关系，得到彩色点云：

(21)$ \left. \begin{gathered} P(x,y,z)\xleftarrow{{{P_{{\mathrm{ref}}}} = {\mathrm{Con}}\;({P_{{\mathrm{def}}}},{\mathbf{\varphi} ({\boldsymbol{X}})}})}P'(x,y,z), \\ P'(x,y,z)\xrightarrow{{{\mathrm{Con}}\;((A,H),(u,v))}}T(u,v). \\ \end{gathered} \right\} $

iLM-FFD能够让浮动图高精度、高效率地向基准图变形. 如图4(a)所示为纹理扭曲的浮动图和基准图经iLM-FFD处理得到变形结果图的过程. 可以看出，变形结果图与基准图无明显纹理差异. 由图4(b)可以看出，改进的L-M算法得到的归一化相关系数比梯度下降法得到的更接近标准值，收敛速度更快，相较于标准L-M算法计算效率提高约3倍，经过20次迭代达到标准值.

实验于2023年7月在华中地区某高校进行，利用自研地面三维扫描仪采集不同场景的点云和全景图像数据，包括操场、楼间、室内、树林、广场等多个常见工作环境，每站至少采集2组不同分辨率的点云数据，共20组，部分数据样例如图5所示. 在Windows 10中基于Visual C++编程平台实现并验证所建立点云与全景图像间的映射关系. 点云的高、中、低分辨率为3 mm@10 m、6 mm@10 m、12 mm@10 m，激光器的发散角为0.5 mrad，点云的点位误差为3 mm@100 m，全景图像分辨率为7200 pixel×3 600 pixel，观测范围为360°×300°，实验选用半径为5 cm±1 mm的靶球作为特征基元.

从点云和全景图像中准确提取基元是算法实施的第一步. 如图6所示为10组数据经过本研究所提算法处理后的基元识别和配对结果, 其中N_T为标靶数量. 实验结果表明：1)点云基元识别具有极高的准确率，标靶提取率达到96.7%. 点云中可信标靶均为真实靶球，原因是点云中信息丰富，可将强度、几何信息作为目标识别条件，大大增加了目标识别的准确性. 2)图像基元识别准确率高，但弱于点云识别，可以识别出76.7%的标靶. 由于标靶在图像识别中仅将RGB信息作为标靶判断依据，易受环境影响，有存在干扰目标的可能，影响了识别的准确性. 3)最小角距差法能够准确配对基元，被点云和图像同时识别的标靶均可以成功配对. 这是由于激光雷达与相机同轴测量，在空间球坐标中点云基元三角形与图像基元三角形稳定且具有较高的相似度，容易实现高成功率的基元配对.

为了定量分析映射结果的准确度，计算点云彩色图与点云强度图中对应标靶中心的像素偏移值，定义平均像素偏移为映射误差. 为了说明所建立的数据映射关系稳定性强，利用高分辨率数据建立点云与全景图像间的映射关系，改变测站位置并采集中、低分辨率数据验证映射关系，在3个场景中完成验证. 如图7所示为楼间场景下的映射实验结果，3个场景的高、中、低点云分辨率的映射误差分别为1.07、1.16、1.14 pixel，最大误差分别为1.48、1.62、1.86 pixel. 可以看出，1)高分辨率点云经初始映射得到的彩色点云质量良好，说明靶球中心提取精度高，映射算法合理，仅部分区域受图像畸变的影响出现较小偏移. 映射关系优化后得到的彩色点云无明显色彩偏移. 2)中、低分辨率点云依据映射关系能够直接得到高质量彩色点云，且不同点云分辨率的映射误差相差小于0.1 pixel. 结果说明所构建激光点云与全景图像间的映射关系稳定性强，不受测站和点云分辨率变化的干扰.

算法实现过程逐步提高点云与全景图像的映射精度. 如图8所示为不同算法处理下的映射误差对比. 20组对比数据相关误差（最大误差$ {e_1} $、最小误差$ {e_2} $、平均误差$ {e_{\mathrm{m}}} $、误差标准差$ {e_{\mathrm{r}}} $）的统计结果如表1所示. 可以看出，1)改进Zernike矩下的映射精度较传统Zernike下的初始映射精度有显著提高，平均映射误差减小了61.1%，因此依据一维最大熵得到的阶跃阈值更合理. 2)优化映射关系后，映射精度进一步提升，20组数据中最大误差为1.62 pixel，平均映射误差为0.96 pixel，误差标准差为0.12 pixel，说明映射结果精度较高，具有较好的精度一致性.

为了进一步说明iLM-FFD对点云映射关系的改善效果，在室内不同高度布设15个标靶，对比iLM-FFD处理前和处理后的标靶中心偏移O_F，结果如图9所示. 图中，H为高度角. iLM-FFD对由图像区域性畸变造成的色彩偏移具有修正作用，映射误差最多降低了53.6%，该算法处理无畸变区域的效果不明显.

为了对比iLM-FFD与直接线性变换法、灰度相似法以及基于变种ICP的处理结果和效率，利用这4种方法处理10组不同场景下的测量数据，某场景的处理结果如图10所示. 图中，光斑打在靶球边缘形成的黑圈是靶球在点云中的实际边界，白色圆为依据各算法的映射关系赋靶球色彩的点云，黑圈与白色圆的贴合度能够直观地表现算法的映射效果. 可以看出，在4种算法中，iLM-FFD得到的黑圈与白色圆贴合度更高，建立的映射关系更准确. 4种算法处理10组数据的精度对比与评价如表2所示. 表中，Acc为映射精度，T_T为耗时. 结果表明，1)传统标定方法多使用直接线性变换法解算模型参数，通过设定虚拟中心建立每张图像与点云之间的刚性配准数学模型，包括12个线性变换参数. 全景图像由多张不同位姿的图像拼接而成，标定过程烦琐，涉及的标定参数较多且不独立，模型未顾及如镜头畸变的非刚性误差的影响，映射精度为12.61 pixel. 2)灰度相似法无需建立传感器标定模型，依据全景影像与点云深度图在水平、垂直方向细分区域的灰度相似性完成数据配准. 配准精度与区域分割密度相关，理论上分割密度越大，数据配 准精度越高，但冗余计算量随分割密度而增加；此外，对于色彩重复度高的环境，灰度相似法出错的可能性较大. 实验中以2个像素为分割单元，映射精度为8.84 pixel. 3)变种ICP利用数据特征建立点云与图像的映射关系，将自序列全景影像虚拟生成多视立体密集匹配点云，继而使用变种ICP算法优化密集匹配点云与激光点云数据间3D-3D配准参数，间接优化全景图像与点云间的配准参数. 变种ICP的数据配准精度和自动化程度较高，映射精度为2.34 pixel，但由图像生成密集点云的过程较为复杂，数据密集匹配过程计算量较大，对数据处理终端的配置要求较高，数据映射关系受点云分辨率影响. 4) iLM-FFD实现在未标定传感器外参数的前提下建立点云与全景图像间像素级映射关系，先从点云和图像中自动提取和配对基元，再根据点云与全景图像在空间球坐标下相似性高的特点构建初始映射关系，最后逐像素优化映射关系. 相较于传统标定方法，iLM-FFD自动化程度和映射精度较高，初始映射参数仅2个，优化后得到稳定、准确的数据映射关系，点云与图像的映射精度为1.07 pixel，相比于灰度相似法提高了87.9%，相比于变种ICP算法提高了54.2%. 该结果能够用于生成真彩色点云，便于对地面三维扫描仪的数据质量进行定期检校.

（1）本研究从点云和全景图像中高精度提取靶球基元，利用最小角距差法将点云和全景图像中的基元准确配对. 该方法便于从数据中自动提取和配对基元，避免了出现错误基元对的问题. 相较于标准Zernike矩，基于一维最大熵的改进Zernike矩提取的基元中心的精度更高，初始映射精度提高了61.1%.

（2）混合优化算法对由于图像区域性畸变导致的初始映射偏移具有良好的修正作用. 相较于传统L-M算法，改进的L-M算法将优化效率提高约3倍，此外点云与全景图像的映射误差在初始映射关系的基础上最多降低了53.6%，有效提升了数据质量.

（3）本研究提出直接建立点云与图像映射关系的方法，无需对传感器外参数进行复杂标定即可建立点云与全景图像的映射关系. 相较于传统标定方法，所提方法自动化程度高，初始映射参数仅2个，优化后得到稳定且准确的数据映射关系，不受测站位置和点云分辨率变化的干扰. 在10组不同场景数据的处理结果中，所提方法数据映射精度为1.07 pixel，相比于灰度相似法提高了87.9%，相比于变种ICP算法提高了54.2%；所提方法的计算效率也较对比算法的高.

（4）相比于通过标定传感器间外参数建立数据转换模型，本研究实现了自动化提取和配对基元，建立和优化了映射关系，对设备数据质量的定期检校具有重要意义，对国产智能扫描仪的研发和传感器数据融合具有重要参考价值. 本研究所提算法对初始映射精度要求较高，同时须高精度提取标靶中心，标靶表面点云粗差数量会随着测量距离的增加而增加. 未来计划将抗差估计的方法引入标靶中心提取计算，致力于研究自适应的抗差估计算法拟合标靶，以保证特征点的提取精度.