A robust methods of fitting plane based on LS-WTLS

平面拟合可用于建(构)筑物倾斜监测、边坡变形监测、隧道断面监测、3-D激光扫描点云滤波^[1-2]等工程领域，其精度会对变形分析、点云滤波和3-D建模效果产生重要影响。鉴此，有必要在传统拟合方法基础上，顾及误差数据比例及其大小对拟合结果产生的影响，构造最佳平面拟合方法，提高所求平面参量的准确度。

目前，常用平面拟合方法包括传统最小二乘(least squares，LS)法、总体最小二乘(total least squares，TLS)法^[3-5]、混合总体最小二乘(mixed TLS, MTLS)法^[6]等，以上方法均将各点当做独立等精度观测值处理，由于受观测环境、系统误差等因素影响，观测数据3-D坐标x, y, z 3个方向上均含有误差且各点点位精度均不相同，若采用上述方法进行平面拟合，所得参量解并非平面参量的最或然值。因此，国内外学者在加权总体最小二乘(weighted total least squares，WTLS)模型^[7-8]基础上，顾及各点点位精度差异性，结合基于方差-协方差膨胀的抗差估计，提出了加权总体最小二乘平面拟合算法，同时对抗差估计中的权因子函数进行了深入探讨，构造了多种定权准则及与其相对应的平面拟合方法^[9-11]，一定程度上提高了平面参量解的准确度。但对于各WTLS法，其已虽将平面模型中的常数项当做无需修正项对待，但实际拟合过程中，常数项因受点位拟合权值影响，准确度较低。

本文中提出了一种基于LS-WTLS混合解算平面参量的方法，该方法将平面参量划分为含误差与不含误差(常数项)两部分，综合利用WTLS法、LS法对这两部分参量进行解算，即顾及到各点点位精度，选用WTLS模型，在此基础上，以稳健估计中的IGGⅢ(由中国科学院测量与地球物理研究所(Institute of Geodesy & Geophysics, IGG)提出的一种抗差估计方法, IGGⅢ是第3种衍生方法)方案权因子为定权函数^[11-15]，在参量迭代解算过程中确定各点拟合权重，并自适应的对拟合权阵进行合理修正，提高误差项参量解算精度，与此同时，又以三倍权因子中误差为阈值剔除粗差数据，利用LS法求取常数项，有效抵抗点位拟合权值对常数项解算的影响。

设平面方程为：

式中，a, b, c为待求平面参量，定义a, b为误差项参量，c为常数项参量。下标(i)表示迭代。顾及观测向量误差与系数矩阵误差的平面方程变量误差模型为：

式中，Z为含有随机误差e矩阵的n×1维观测向量；A=[A₁, A₂]为系数矩阵，A₁为误差项系数矩阵，A₂为常数项系数矩阵；ξ=[ξ₁, ξ₂]^T为代估参量矩阵，ξ₁为误差参量矩阵，ξ₂为常数项矩阵：

式中，v表示计算值与实测值之差。在满足${{\left\| \left( \Delta {{\mathrm{A}}_{1}},\Delta \mathrm{Z} \right) \right\|}_{F}}=\text{min}$的约束条件下(${{\left\| \text{ }\!\!\cdot\!\!\text{ } \right\|}_{F}} $为F范数)，可求得ξ的值。

参照WTLS模型迭代解算方法，可将(2)式可改写成如下形式：

式中，A(_i)=A-E(_i)，δξ=ξ_(i+1)－ξ_i，根据(4)式构造拉格朗日目标函数：

式中，λ为n×1维拉格朗日乘数, I是单位矩阵，T表转置。通过对(5)式求导、求极值即可计算出λ及参量ξ_(i)的代估值。

顾及粗差对解算精度产生的影响，通过引入稳健估计IGGⅢ方案权因子自适应的修正拟合权阵以确保权阵准确性，又通过设定阈值剔除粗差数据以确保数据可靠性，提出了基于LS-WTLS混合解算平面参量的方法(简称LS-WTLS混合法)，其是在通过WTLS法得到代估参量${{\mathit{\boldsymbol{\hat \xi }}}_{(i)}} $基础上，提取${{\mathit{\boldsymbol{\hat \xi }}}_{(i)}} $中的含误差部分参量${{\mathit{\boldsymbol{\hat \xi }}}_{1(i)}}$，再利用LS法计算得到常数项${{\mathit{\boldsymbol{\hat \xi }}}_{2(i)}}$，完成所有参量解算。上标^表示实际值的最或然值(下同)。该方法具体解算步骤见下。

(1) 利用LS法计算平面参量估值$ {{\mathit{\boldsymbol{\hat \xi }}}_{(0)}} = {\left[ {{{\mathit{\boldsymbol{\hat \xi }}}_{1(0)}}, {{\mathit{\boldsymbol{\hat \xi }}}_{2(0)}}} \right]^{\rm{T}}}$，之后分别按照下式设定系数矩阵A的列向量权阵P₀及A的行向量初始权阵P_A(1)、观测向量初始权阵P_Z(1):

(2) 计算参量$ {{\mathit{\boldsymbol{\hat \xi }}}_{(1)}} = {\left[ {{{\mathit{\boldsymbol{\hat \xi }}}_{1(1)}}, {{\mathit{\boldsymbol{\hat \xi }}}_{2(1)}}} \right]^{\rm{T}}}$的迭代初始值，其中${{{\mathit{\boldsymbol{\hat \xi }}}_{1(1)}}}$为$ {{\mathit{\boldsymbol{\hat \xi }}}_{(1)}}$的含误差部分，且定义$ {{{\mathit{\boldsymbol{\hat \xi }}}_{2(1)}}}={{{\mathit{\boldsymbol{\hat \xi }}}_{2(0)}}}$：

式中，Q_Z, Q_A, Q₀分别为P_Z, P_A, P₀的广义逆。

(3) 计算$ {{{\mathit{\boldsymbol{\hat \mu }}}_{{\rm{(}}\mathit{i}{\rm{)}}}}}$，${{{\mathit{\boldsymbol{\hat \lambda }}}_{{\rm{(}}\mathit{i}{\rm{)}}}}} $，$ {{\mathit{\boldsymbol{\hat v}}}_{{\rm{(}}\mathit{i}{\rm{)}}}} $：

(4) 计算${{\mathit{\boldsymbol{\hat \xi }}}_{1(i + 1)}}$，即$ {{\mathit{\boldsymbol{\hat \xi }}}_{(i + 1)}}$的含误差部分：

(5)利用(10)式计算各残差V_(i):

(6) 根据V_(i)及IGGⅢ方案(见(11)式)计算权因子ω_(i)，令P_(i)=P_(i)ω_(i)，之后根据准则重新设定权阵P_A(i+1), P_Z(i+1):

式中，V_(i)为残差，σ为中误差，k₀=1.0~1.5，k₁=2.5~3.0，本文中取k₀=1.5, k₁=2.5。

(7) 计算权因子中误差σ_ω(i)，并剔除ω_(i)>3σ_ω(i)的点，同时根据下式计算$ {{\mathit{\boldsymbol{\hat \xi }}}_{2(i + 1)}}$，综合步骤(4)可得参量:

(8) 重复步骤(3)~步骤(7)，直到$\left| {{{\mathit{\boldsymbol{\hat \xi }}}_{(i + 1)}} - {{\mathit{\boldsymbol{\hat \xi }}}_{(i)}}} \right| < {\delta _0} $(δ₀为给定阈值，本文中取δ₀=10^－6)。

经迭代运算得到拟合平面的参量解之后，可利用(13)式计算单位权中误差$ {{\hat \sigma }_0}$及平面拟合精度$ {{\hat \sigma }_p}$进行精度评定：

式中，d_(i)为点至拟合平面的距离。

为对本文中所构建平面拟合方法的适用性及优越性进行验证，分别利用其对模拟平面数据、实测平面数据进行拟合。

设要拟合的平面方程为Z=3X+4Y+5，依次抽取4组观测数据，并分别按照0%, 5%, 10%, 20%的比例加入粗差(大小为3σ₀~5σ₀，σ₀为每组观测数据的中误差)，每个观测数据中，X, Y分别为[0, 1000], [0, 2000]区间内的随机整数，相应的Z为[5 11005]区间内的随机整数，同时在模拟数据中加入随机误差e，且随机误差e服从均值为0、标准差为σ²I的正态分布(σ=0.3，I为单位阵)。

分别利用LS法、TLS法(奇异值分解(singular value deeoniposition, SVD))、LS-TLS法、IGGⅢ-LS-TLS法^[5]、LS-WTLS混合法(本文中算法)对上述4组观测数据进行拟合实验，各方法计算得到的平面参量及精度评定指标如表 1所示。

表 1中，当粗差比例为0%时，各方法所得参量解非常接近实际参量值，单位拟合中误差与平面拟合精度较小，由此说明各方法的拟合效果均比较好，仅由于数据取位原因导致各方法所得参量值与实际参量值之间存在一定偏差。

当粗差比例由5%递增到20%时，各方法所得参量解与实际参量值的偏差越来越大，单位权中误差与平面拟合精度呈递增趋势，由此说明，拟合效果随粗差的增加而变差。其中，LS法受粗差影响，拟合效果较粗差比例为0%时明显下降；TLS法过多顾及到系数矩阵不含误差部分，整体拟合效果最差；LS-TLS法虽顾及平面参量的常数项不含误差，对两部分参量进行独立解算，但其未考虑到各观测值精度差异，因此拟合效果不稳定；IGGⅢ-LS-TLS法综合考虑了含误差参量与常数项的区别，对其进行了独立解算，同时根据各观测值的精度，确定其参量平面拟合的权重，因此获得了较好的拟合效果；LS-WTLS混合法在IGGⅢ-LS-TLS法的基础上，引入加权总体最小二乘模型，以拟合残差为依据、IGGⅢ方案为定权准则，在含误差参量迭代解算过程中，其可自适应地修正观测向量权阵及系数矩阵权阵，以获取最能合理反映各观测值精度的权值，同时又以三倍权因子中误差为阈值，剔除观测数据中的异常数据，再利用LS法计算得到常数项的值，最终，通过有限次的迭代计算，获得了最为可靠参量解，该方法较好抵抗了粗差干扰，各项拟合指标均优于其它算法，拟合精度最高。

表 2中^[16]，X为坝体温度实测值、Y为水位压力实测值、Z为大坝水平位移实测值(X, Y, Z无具体单位), 3个向量构成Z=aX+bY+c的平面关系。再次利用LS法、TLS法(奇异值分解(SVD))、LS-TLS法、IGGⅢ-LS-TLS法、LS-WTLS混合法(本文中算法)对表 2中的实测数据进行拟合实验，所得平面参量及精度评定指标如表 3所示。

由于实际平面数据的受测量仪器、外界环境、人为干扰等多种因素影响，因此数据纯度较低，其包含的偶然误差及粗差的比例存在不确定性。由表 3可看出，TLS法将常数项当作误差项处理，导致参与拟合的系数矩阵不准确，因此，拟合效果最不理想，其拟合结果不具参考性；LS法与LS-TLS法的计算结果、拟合精度评定指标$ ({{\hat \sigma }_0}, {{\hat \sigma }_{\rm{p}}})$较为接近，由此说明，当未顾及各观测值精度的差异时，LS-TLS法虽将两部分参量独立解算，但其仅仅保证了系数矩阵的相对正确性，未对最终参量的准确度产生积极影响；IGGⅢ-LS-TLS法与LS-WTLS混合法所得参量解及拟合精度评定指标同样较为接近，说明两种方法都比较可靠，拟合精度较高，但对比表中数据可发现，LS-WTLS混合法以点与拟合平面的相关关系为原则，结合IGGⅢ方法合理设定拟合权阵、剔除观测数据粗差，获得了最为精确、可信度最高的参量解，相对于LS法、TLS法、LS-TLS法、IGGⅢ-LS-TLS法，该方法的单位权中误差分别提高了53.6%, 195.0%, 47.5%, 5.1%，平面拟合精度分别提高了49.4%, 179.3%, 48.7%, 46.99%，进一步表明该方法的拟合精度最高，效果最好。

提出了一种对平面含误差与不含误差(常数项)两部分参量进行独立解算的LS-WTLS平面拟合方法。在误差项参量解算过程中，顾及平面数据各点精度存在差异，本文中依据稳健估计中各权函数特点及其适用范围，选用IGGⅢ方案对加权总体最小二乘平面拟合方法进行改进，其可在模型参量解算过程中，通过计算点与平面模型的相关关系，自适应地调整各点拟合权值，优化拟合权阵；对于常数项，先以三倍权因子中误差为阈值剔除观测数据中的粗差，再利用LS法进行解算。模拟平面数据及实测平面数据的拟合实例表明，该方法具备较好的可行性及优越性，利用该方法拟合得到的平面参量解可靠性更高。

基于LS-WTLS的稳健平面拟合方法较最小二乘法、总体最小二乘法、混合总体最小二乘法更为稳健，但在解算过程中，当数据量过大时，由于权值的自适应修正过程较为复杂，迭代计算较为繁琐，解算所需时间较多，如何提高解算效率值得进一步深入研究。