结合RANSAC算法与传统相似变换模型的稳定性分析

刘忠贺，李宗春，郭迎钢，何 华，赵文斌

(1.信息工程大学，河南 郑州 450001；2.中国科学院上海高等研究院，上海 201210)

加速器装置运营期间，磁铁等关键元件会受到外界诸多因素影响发生位移，如不能及时发现并做出调整，则装置运转性能会受到影响，严重时会导致粒子运行轨道发生偏转[1-2]。元件是否发生变形只能依靠控制网对其进行监测和调整，所以定期对控制网进行复测并对网点进行稳定性分析是一项非常重要的工作[3]。

德国测量学者Pelzer提出的平均间隙法[4]利用统计检验的方法对控制网几何图形作整体检验，判断控制网在两期观测间是否发生显著性变化，依次寻找变形点并将其剔除，直到检验通过[5-6]。基于平差后的两期坐标之间存在一定的平移、旋转、尺度信息，可采用传统相似变换(Traditional Similarity Transformation, TST)法计算两期控制网之间的变换参数，进而求出变形点的位移信息，但易受到位移较大点的影响，计算的模型参数不准确，稳定性分析结果不可靠。基于此，陈永奇[7]提出了迭代加权相似变换法(Iterative Weighted Similarity Transformation, IWST)，以所有点位移的一次范数最小为目标函数[8]，采用变权迭代的方式来计算点的位移量，在TST法基础上，增加对变形点的抗差性，提高稳定点判别的正确率。Duchnowski Robert[9]将基于中值的稳健估计理论应用到控制网的稳定性分析中，在一定程度上抵抗了粗差的影响。Krzysztof Nowel[10]提出观测值差异分析变形的稳健估计方法，是通过消除观测数据中存在的粗差来提高稳定性分析的可靠性。A.R.Amiri-Simkooei等人[11]同时将两期观测数据进行平差解算，将每一个点当作变形点，其变形量设为待求参数进行求解，依次剔除位移最大的点，然后进行整体检验，直到检验通过。该方法计算繁琐，类似于单点检验方法，并且有时会出现最大变形点判断错误的情况。郭迎钢[12]在加速器隧道控制网稳定性分析过程中加入了网点的空间灵敏度椭球，并分析了控制网的可监测性，结合整体稳定性检验和平均间隙法对隧道控制网进行了稳定性分析，但在求解网点变形量时未顾及较大变形点对变换参数的影响。

TST模型在求解相似变换参数时易受到较大变形点的影响，不具有抗差性，若能通过一定的方法将变形较大的点剔除，利用剩余较为稳定的点求解变换参数，则可以得到更加准确的网点变形信息。随机抽样一致性(Random Sample Consensus，RANSAC)算法[13]可用于数据的提纯，减少或避免异常点对模型的干扰[14-17]。文中结合RANSAC算法与TST模型，使该模型在应用时具备抗差性，并通过试验进行验证。

1 函数模型

1.1 TST模型

两期控制网之间可通过公共点求解转换参数，七参数坐标转换模型可表示为：

(1)

式中:(X，Y，Z)与(x，y，z)为同一点在两期观测下的坐标；两期观测之间绕X，Y，Z轴的旋转角分别为ω，φ，κ；(Δx，Δy，Δz)为平移量；m为尺度。两期观测间旋转矩阵R为：

(2)

将式(1)线性化，得：

(3)

式中：(vX,vY,vZ)为残差；(X0，Y0，Z0)为第二期坐标系下的坐标(x，y，z)根据两期之间转换参数初值转换至第一期坐标系下坐标。

TST模型可表示为：

(4)

式中：H表示系数矩阵；t=[dmdωdφdκdΔxdΔydΔz]T；d为位移向量。

1.2 IWST模型

由于TST模型在求解变换参数时容易受到变形点的影响，所以在其基础上，陈永奇提出了IWST模型[7]。在求解变换参数时，以网点位移向量d的一次范数最小为准则，根据位移的大小不断调整权值，通过迭代求解最终的位移。目标函数可表示为：

∑|di|=min(i=1,2,…,3n).

(5)

初始计算时，取权阵P(k=1)=I为单位阵，通过TST模型计算出位移向量d(k=1)。

对于第k+1次迭代，第i点权阵为：

(6)

为了避免权值无限大，分母上加了一个微小量ε。迭代过程如以下形式：

(7)

其中

S(k)=I-H(HTP(k)H)HTP(k).

(8)

式中：ΔX为两期坐标差；Q11，Q22分别为第一期和第二期平差后的协因数阵。

当前后两次迭代后的位移量之差小于某一阈值μ时，停止迭代，即:

|d(k+1)-d(k)|<μ.

(9)

最后，根据最终位移量，利用T检验法对网点的稳定性进行判断。

1.3 文中方法

1.3.1RANSAC算法思想

RANSAC算法是一种通过迭代计算过程来寻求样本中满足某一正确模型对应的最大内点集，并利用该集合中样本重新估算模型的过程[18]。基本过程：

1)给定一个数据集，从中随机抽取t个样本并计算待求参数模型初值；

2)设置阈值ε，利用求出的参数模型初值验证集合中剩余样本的模型误差是否超限。阈值范围内的样本归为内点；否则，归为外点；

3)统计内点集中样本数量，再次随机抽取t个样本，重复上述步骤；

4)经过r次迭代次数后，选取迭代过程中的最大内点集，判断该集合样本数量是否达到k值，若满足条件，利用该集合计算参数模型并筛选内点集，计算最终的模型；否则，算法失败。

上述过程涉及到的参数根据实际情况设定，代表含义如下：

t：计算模型参数所需最少样本数量；

ε：决定一致集样本数量多少；

k：表征正确模型对应一致集的最少样本数；

r：与外点所占样本集比例有关，决定模型的可靠程度。

1.3.2 稳定性分析流程

最小二乘估计是一种无偏估计，在无系统误差、粗差的情况下，可得到最优化参数估计结果，且计算过程简便。然而，该估计方法对粗差非常敏感，当观测数据中存在较大变形点时，变形量极易发生偏移、转移或分配到其他点的坐标信息中，对结果造成影响。为了得到准确的相似变换参数，必须将控制网中的变形点予以剔除，使其不参与变换参数的计算，通过稳定的点来估计正确的转换模型，在此基础上分析网点的稳定性，求解其绝对位移量。前面已经介绍了利用RANSAC算法提高坐标系之间转换参数精度的方法，将其应用到控制网的稳定性分析中来，对TST方法进行改进。基本过程如下：

1)设隧道控制网两期观测误差方程为：

(10)

采用重心基准，分别求解两期观测坐标X1,X2；单位权中误差σ1,σ2；协因数阵Q11,Q22。

2)经过平差后的两期坐标之间存在一定的平移、旋转及尺度缩放关系，通过公共点转换的方法对各参数进行求解，计算过程中采用RANSAC算法筛选相对稳定点，具体流程见图1。

筛选流程主要分为两步，首先通过迭代筛选出最大内点集，然后利用该内点集重新估算参数模型，迭代筛选最终的稳定点集。

3)利用筛选出的公共点求解变换参数，将所有点代入该参数模型，计算网点的位移量。

4)通过T检验法对每一个点进行稳定性判断，构造统计量。

(11)

(12)

式中：f1,f2为两期观测自由度。

选择一定显著性水平α，对网点稳定性进行判断。

Ti≤Fα(wi,f1+f2).

(13)

满足式(3)的点即为稳定点，反之，则为不稳定点。

TST的改进之处在于第2步求解转换参数过程中引入了RANSAC算法，使得TST模型在应用时具有稳健性，通过筛选并剔除某些变形点，进而提高转换参数的精度，使得两期坐标间的转换关系更加准确可靠。

2 试验与分析

2.1 模拟试验

为了验证文中方法的可行性，以粒子加速器隧道控制网为例，设计了模拟试验。

2.1.1 数据

在4 m×2.5 m×4 m范围内模拟了一组数据，

图1 相对稳定点筛选流程

包含S1～S6共6个测站以及A1～A16共16个控制点，为了使测量条件更加贴近实际，点位分布均匀且测站完全在控制点所包含的空间内部，网型及试验场景如图2所示。

图2 试验场景

理论数据见表1，并在各点的坐标分量上加入方向、大小不同的位移量，作为第二期理论数据。

根据表1，对A1,A4,A7,A9,A11,A13,A15共7个点的某一坐标分量加入了位移，考虑到5 m范围内，点位测量精度不超过0.05 mm，所以加入的位移量最小为0.05 mm，最大为3 mm，变形点占比43.75%。两期数据的理论坐标设计完成后，按照测角精度±(15+6×10-6·D)μm，测距精度为±0.5·Dμm/m的标称精度向数据添加随机误差，作为两期对比试验数据。

2.1.2 试验与分析

分别采用TST模型、IWST模型及文中方法对两期模拟数据进行处理，分析控制网的变形及稳定情况。计算过程中，取显著性水平α=0.05，在对各点进行T检验时，F0.05(3,420)=2.626 1，3种方法分析结果见表2。

表1 网点及测站设计坐标 mm

由表2可知，在16个点组成的控制网中，TST模型受到变形点的影响，最终稳定性分析效果较差，准确率达到50%。结合RANSAC算法与TST模型，最终结果接近实际情况，分析变形点结果准确率接近90%，误判了A5点。对于IWST模型，误判了两个点，其分析结果劣于文中方法。为了检验3种方法求解变换参数的准确性，对比其效果优劣，将3种方法计算的变形点位移同理论值进行对比，见表3。

表2 不同模型网点稳定性分析结果对比

表3 变形点位移量与理论值对比 mm

从表3可以发现，从整体效果来看，文中方法所计算的变形点位移量最接近真实变形情况。对于所有变形点，文中方法计算的A1点位移与理论值相差最大，达到0.025 mm，分析其原因，可能是在模拟数据时，添加的随机误差最大不超过0.05 mm，而在该点上随机误差大小接近0.025 mm，文中方法虽然是一种稳健估计方法，但作用对象是平差后的点坐标，对于平差后坐标本身所带有的误差并没有消除作用。然而，剩余变形点变形量均与理论变形值相差甚小。

2.2 实例分析

2.2.1 数据来源

文中采用2017年与2019年的SSRF复测数据，对两期控制网同名点进行稳定性分析。采用重心基准，分别对两期观测数据进行平差，结果见表4。

表4 平差结果

表4中自由度表示多余观测个数。经查询，两期测量数据中同名点数量共65个，平差后的同名点分布如图3所示。

图3 同名点分布图

从图3可以看出，同名点为储存环部分网点，以这些点为基础，采用不同方法分析其稳定性。

2.2.2 试验与分析

利用TST模型、IWST模型及文中方法对65个同名点进行处理。TST模型和IWST模型直接对所有点进行处理，文中方法处理的基本过程为：①筛选同名点中部分相对稳定点，将距离阈值设为0.3 mm，最终得到35个相对稳定点；②利用这些点计算变换参数；③计算剩余点在此变换参数下的位移量。计算时，取显著性水平α=0.05，3种方法分析结果见表5。

表5 不同方法稳定性分析结果

从表5可以发现，3种方法分析结果略有不同。3种方法均得出BD512点位变形点，且位移量基本相等，在3 mm左右；相比IWST模型，TST模型及文中方法分析结果中，均多出一个变形点，且不是同一点，其位移量较小。3种方案处理后的65个同名点变形量如图4所示。

图4 不同模型所得变形量

从图4直观地发现，3种方案处理结果均发现第42号点(BD512)发生了变形，且十分明显，而其余点变形相对较小。除BD512点外，1号点(BD303)变形最大，在2.5 mm左右，向后依次减小，趋于平稳，3种方法所得结果一致，说明控制网各网点之间相对稳定。为了进一步验证3种方法处理结果的正确性，利用SA软件对两期数据进行处理，将两期数据中的同名点进行最佳拟合变换，得到的变形量见图5。

图5 最佳拟合变换

从图5中可以看出，第42号点发生了较大变形，与以上3种方法处理结果一致。综合几种方法，该点可能发生变形，需要到实地勘测变形情况，确定其是否真正发生变形，对该点坐标信息进行更新并分析变形原因。

模拟试验虽然初步验证文中方法的有效性，但由于条件受限，加上TST模型及文中方法分析结果中所多出的变形点，对于实例分析部分，并未对相应变形点进行实地考察，不能确定点位是否真正发生了变形，只能从数据分析结果来判断，是文中的不足之处。

3 结束语

以粒子加速器隧道控制网为例，通过模拟试验与实例分析，验证文中方法在隧道控制网稳定性分析中有效性(亦可推广至其他场景)，并得出以下3点结论。

1)结合RANSAC算法与TST模型，使得该模型在计算相似变换参数时具备稳健性，能够避免较大变形点对变换参数的影响；

2)同TST模型及IWST模型相比，文中方法稳定点判别准确率最高，且计算得到的位移量与实际情况更加相符；

3)文中方法在TST模型基础上提高了坐标转换参数的精度和可靠性，在判断控制网稳定性方面更具可行性。