受扫描范围、物体遮挡等限制，激光扫描仪难以通过单次扫描获得完整模型，而是需要将不同视角的点云数据统一到同一坐标系中。通常在2台仪器的公共扫描范围内放置至少3个标靶球，先定位其球心坐标，再通过刚体变换实现配准。然而粗差会使球心定位出现明显偏差，因此在含有异常值的点云数据中，准确估计球心坐标是点云配准的关键^[1]。

在点云数据只含有随机误差的情况下进行标靶球定位时，一般采用最小二乘法(LS)计算其几何中心^[2]。然而系数矩阵中也含有存在误差的观测值，需要使用总体最小二乘法(TLS)建立EIV模型对球心进行定位。鲁铁定等^[3]使用奇异值分解方法得到TLS估计值；陈玮娴等^[4]基于点云入射角对坐标进行定权；陶武勇等^[5]使用双非线性TLS法进行2次线性化处理获得参数解。但上述TLS法均需要进行迭代求解，每次迭代均存在大量的矩阵运算，算法复杂度较高。

由于点云数据中难免含有粗差，因此需要将稳健估计理论应用到标靶球定位中。Zhou等^[6]提出基于M估计的稳健球心估计方法，能够确定含有粗差的点云数据中心坐标；王乐洋等^[7]提出适用于球心定位的稳健总体最小二乘法(RTLS)，在迭代过程中，需根据3倍中误差准则剔除粗差。然而，M估计和RTLS估计均依赖于迭代初值，容易收敛至错误的估计值，崩溃率有限，且TLS法需要迭代求解，半径估计值有偏。

针对上述问题，潘国荣等^[8]使用基于忠实距离的代数拟合方法，采用Pratt方法^[9]求解参数。结果表明，该方法基本等效于几何拟合法，且在计算上更具优势。在此基础上，本文进一步顾及算法的稳健性，融合加权Hyper方法^[10]和截断最小二乘法(LTS)实现快速稳健的球心定位。与2种常用的稳健球心定位算法相比，本文算法的崩溃率和计算效率都有显著提高。

1 加权Hyper方法

加权Hyper方法具有以下优点：1)具有与几何拟合相同的估计精度；2)参数估计值满足无偏性；3)计算简单，仅需计算五阶矩阵的特征值便能定位球心。常见的球面拟合准则分为几何准则和代数准则，几何准则必须进行迭代求解，而代数准则基于球的代数方程，无需迭代。代数准则的目标函数为：

$ \begin{gathered} \sum\limits_{i=1}^n f^2\left(x_i, y_i, z_i\right)= \\ \sum\limits_{i=1}^n\left[A\left(x_i^2+y_i^2+z_i^2\right)+\right. \\ \left.B x_i+C y_i+D z_i+E\right]^2=\min \end{gathered} $

(1)

式中，(x_i, y_i, z_i)为球面离散点观测值；A、B、C、D、E均为待估参数，待估参数与球心c和半径R的转换可表达为：

$ \boldsymbol{c}=\left[-\frac{B}{2 A}, -\frac{C}{2 A}, -\frac{D}{2 A}\right]^{\mathrm{T}} $

(2)

$ R=\sqrt{B^2+C^2+D^2-4 A E} / 2 A $

(3)

文献[8]中的最小化忠实距离平方和为$\sum\limits_{i=1}^n f^2 $，基本等效于几何方法。本文进一步顾及观测值精度，将目标函数定义为：

$ \sum\limits_{i=1}^n \frac{p_i f^2\left(x_i, y_i, z_i\right)}{4 R^2}=\min $

(4)

式中，p_i为点i的权值，本文借鉴入射角定权方法^[4]，设置p_i=cosθ_i。

令$k_i=x_i^2+y_i^2+z_i^2 $，Pratt^[9]和Taubin^[11]分别将式(2)转换为以下带约束的优化问题：

$ \boldsymbol{u}^{\mathrm{T}}\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{X}\right) \boldsymbol{u}=\min \quad \text { s. t. } \boldsymbol{u}^{\mathrm{T}} \boldsymbol{N}_P \boldsymbol{u}=1 $

(5)

$ \boldsymbol{u}^{\mathrm{T}}\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{X}\right) \boldsymbol{u}=\min \quad \text { s. t. } \boldsymbol{u}^{\mathrm{T}} \boldsymbol{N}_T \boldsymbol{u}=1 $

(6)

式中，$ \boldsymbol{u}=\left[\begin{array}{lllll} A & B & C & D & E \end{array}\right]^{\mathrm{T}}, \boldsymbol{X}=$ $\left[\begin{array}{ccccc} k_1 & x_1 & y_1 & z_1 & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ k_n & x_n & y_n & z_n & 1 \end{array}\right], \boldsymbol{P}=\left[\begin{array}{ccc} p_1 & 0 & 0 \\ \vdots & \ddots & \vdots \\ 0 & 0 & p_n \end{array}\right], $$\boldsymbol{N}_P=\left[\begin{array}{ccccc} 0 & 0 & 0 & 0 & -2 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ -2 & 0 & 0 & 0 & 0 \end{array}\right], $$\boldsymbol{N}_T=\left[\begin{array}{ccccc} 4 \bar{k} & 2 \bar{x} & 2 \bar{y} & 2 \bar{z} & 0 \\ 2 \bar{x} & 1 & 0 & 0 & 0 \\ 2 \bar{y} & 0 & 1 & 0 & 0 \\ 2 \bar{z} & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{array}\right] $。文献[10]对上述2种方法进行误差分析后发现，2种方案的半径估计值都有偏。进一步中和2种方法的固定偏差，提出无偏的加权Hyper方法，其目标函数为：

$ \boldsymbol{u}^{\mathrm{T}}\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{X}\right) \boldsymbol{u}=\min \quad \text { s. t. } \boldsymbol{u}^{\mathrm{T}}\left(\boldsymbol{N}_P-2 \boldsymbol{N}_T\right) \boldsymbol{u}=1 $

(7)

使用奇异值分解法分解矩阵$\left(\boldsymbol{N}_P-2 \boldsymbol{N}_T\right)^{-1}\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{X}\right) $，其中最小非负特征值对应的特征向量为u的估计值，具体证明见文献[12]。

2 截断最小二乘法LTS

尽管加权Hyper方法具备很多优点，但是容易因粗差而崩溃。为消除粗差影响，根据LTS的估计准则引入截断参数h，抽样得到子集并检验其目标函数。LTS方法的线性回归函数模型可表示为：

$ \boldsymbol{y}=\boldsymbol{A x}+\boldsymbol{e} $

(8)

式中，y、A分别为观测向量和系数矩阵，x为待估参数向量，e为独立等精度的误差向量。LS为最佳线性无偏估计，对粗差极其敏感，其准则为$\sum\limits_{i=1}^n e_i^2=\min $。为此，文献[13]提出LTS方法，选择一个由h个较小残差观测值组成的子集，并最小化子集中的误差平方和，其目标函数为：

$ \sum\limits_{i=1}^h e_i^2=\min \left(e_1^2<e_2^2<\cdots<e_n^2\right) $

(9)

LTS方法通过抽样检验求解参数，在数据被污染一半的情况下仍能抵抗粗差的干扰，具有渐近正态性和高崩溃率。

3 快速稳健的标靶球定位方法

本文结合Hyper方法和LTS方法提出一种新的快速稳健的标靶球定位方法。根据粗差比例κ设置子集中观测值的数量为$h=\lceil(1-\kappa) n\rceil $(符号$\lceil\rceil $表示向上取整)。若数据污染率未知，则设置κ=50%。为加快运算速度，随机选择m=4个点确定球的参数。根据概率论，需保证选择的4个点都不是粗差，抽样的次数l为：

$ l=\frac{\log (1-p)}{\log \left(1-(1-\kappa)^m\right)} $

(10)

式中，p为选到不受粗差污染的子集概率，实验中设置p=99%。重复以下3个步骤l次，选择目标函数最小时的估计值作为输出值：

1) 随机选择不在一个平面上的m=4个点，固定参数A=1，利用球面上的4点直接确定球面参数u。进一步计算点i(i=1, 2, …, n)的加权忠实距离平方r_i²：

$ r_i^2=\frac{p_i f^2\left(x_i, y_i, z_i\right)}{4 R^2} $

(11)

2) 将所有数据点的r_i²进行升序排列($r_1^2<r_2^2<\cdots<r_n^2 $)，保留h个加权忠实距离平方最小的点形成一个子集。

3) 使用加权Hyper方法计算步骤2)中的子集，得到u的估计值。将加权忠实距离平方按照升序排列，将前h个最小的加权忠实距离平方和作为目标函数。

步骤1)是为了减少抽样次数。根据式(10)可知，若直接设置抽样数量m=h，则抽样次数l会迅速增加；步骤2)的主要目的是保证步骤3)中有更多的正常观测值参与平差计算，增强有效性；步骤3)中的加权Hyper方法不同于几何拟合算法，无需迭代即可定位球心。该方法的目标函数为$\sum\limits_{i=1}^h r_i^2=\min $，多余观测数为(h-m)，最后评定精度为：

$ \hat{\sigma}=\sqrt{\sum\limits_{i=1}^h r_i^2 /(h-m)} $

(12)

4 实验与分析

使用RIEGL VZ-1000激光扫描仪扫描R=7.25 cm的标靶球，获得835个点云数据，经重心化处理后的球心为$ \tilde{\boldsymbol{c}}$=[0.031 8 0.026 1 -0.022 1]^T。在此基础上分别添加平面分布和密集分布的2种粗差，图 1为20%污染率的点云分布图。使用3种稳健估计方法求解参数：本文快速稳健球心定位算法(方案1)、基于粗差探测法的RTLS法^[7](方案2)和M估计法^[6](方案3)，其中方案2和方案3的最大迭代次数为200次。设置9组实验，粗差比例从5%依次增加至45%，每组实验重复N=100次。图 2和图 3分别为不同粗差比例下3种方案的单位权中误差$\hat{\sigma} $和平均球心偏差$\mathrm{D}(\hat{c})=\frac{1}{N} \sum\limits_{i=1}^N\left\|\hat{c}_i-\tilde{c}\right\| $。

由图 2、3可见，当污染率小于20%时，各方案均能较为准确地定位球心，表明3种算法都具有一定的稳健性。其中方案1和方案2的差异并不显著，说明本文基于忠实距离准则的算法等效于基于几何距离准则的TLS方法；当污染率增加到25%~30%时，方案2的2个指标大于其他2种方案；当污染率大于30%时，方案3的2个指标也明显增大，说明粗差探测法和M估计法的崩溃率有限。本文算法的2个指标波动较小，实验中算法崩溃率高达45%。

进一步比较3种算法的计算效率，图 4和图 5分别为3个方案的运行次数和平均运行时间t的对数lgt。由图 4可见，方案2和方案3的迭代次数较多，而方案1属于抽样检验算法，无需设置初值，运行次数最少。由图 5可见，3种算法的运行时间量级分别为10^-2 s、10¹ s和10⁰ s，本文算法运行时间远小于其他2种方案。这是因为在迭代过程中，RTLS需要进行复杂的大型矩阵乘法和求逆运算，点云数量较多，耗时较长；M估计法在迭代过程中需要重新线性化和定权，增加了参数解算时间；本文算法通过Hyper方法直接得到估计值，过程中仅进行五阶矩阵的运算，运算量显著减少，因此运行时间的量级明显低于其他2种算法。

选择WHU-TLS数据集^[14]中校园场景的球形建筑数据测试本文算法的性能。该数据有23 440个点，点云分布和球面拟合结果如图 6所示，其中下方环状点云和上方聚集点云在球心定位中属于粗差。从拟合效果上可以看出，本文算法能够有效区分粗差和正常观测值，具有较强的稳健性。实验抽样次数为19，总共耗时0.14 s，说明本文算法具有极高的运算效率，且效率并没有随点云数量的增加而快速下降。

5 结语

本文结合LTS方法和加权Hyper方法提出一种快速稳健的球心定位算法。该算法减少了传统算法复杂矩阵的迭代计算，能够有效抵御粗差干扰。与2种传统算法相比，本文算法的崩溃率更高且解算时间更短，适合处理带有粗差的标靶球点云数据。