-
设平面方程为:
$ {\mathit{\boldsymbol{Z}}_{(i)}} = a{\mathit{\boldsymbol{X}}_{(i)}} + b{\mathit{\boldsymbol{Y}}_{(i)}} + c, (i = 1, 2, \ldots , n) $
(1) 式中,a, b, c为待求平面参量,定义a, b为误差项参量,c为常数项参量。下标(i)表示迭代。顾及观测向量误差与系数矩阵误差的平面方程变量误差模型为:
$ \mathit{\boldsymbol{Z}}-\mathit{\boldsymbol{e}} = ({\mathit{\boldsymbol{A}}_1} + {\mathit{\boldsymbol{E}}_{{A_1}}}){\mathit{\boldsymbol{\xi }}_1} + {\mathit{\boldsymbol{A}}_2}{\mathit{\boldsymbol{\xi }}_2} $
(2) 式中,Z为含有随机误差e矩阵的n×1维观测向量;A=[A1, A2]为系数矩阵,A1为误差项系数矩阵,A2为常数项系数矩阵;ξ=[ξ1, ξ2]T为代估参量矩阵,ξ1为误差参量矩阵,ξ2为常数项矩阵:
$ \begin{matrix} \underset{n\times 1}{\mathop{\boldsymbol{Z}}}\,=\left[ \begin{matrix} {{z}_{1}} \\ {{z}_{2}} \\ \vdots \\ {{z}_{n}} \\ \end{matrix} \right],\underset{n\times 2}{\mathop{{{\boldsymbol{A}}_{1}}}}\,=\left[ \begin{matrix} {{x}_{1}} &{{y}_{1}} \\ {{x}_{2}} &{{y}_{2}} \\ \vdots &\vdots \\ {{x}_{n}} &{{y}_{n}} \\ \end{matrix} \right],\underset{n\times 1}{\mathop{{{\boldsymbol{A}}_{2}}}}\,=\left[ \begin{matrix} 1 \\ 1 \\ \vdots \\ 1 \\ \end{matrix} \right], \\ \underset{n\times 2}{\mathop{\boldsymbol{E}}}\,=\left[ \begin{matrix} {{v}_{{{x}_{1}}}} &{{v}_{{{y}_{1}}}} \\ {{v}_{{{x}_{2}}}} &{{v}_{{{y}_{2}}}} \\ \vdots &\vdots \\ {{v}_{{{x}_{n}}}} &{{v}_{{{y}_{n}}}} \\ \end{matrix} \right],\underset{n\times 3}{\mathop{\boldsymbol{e}}}\,=\left[ \begin{matrix} {{v}_{{{z}_{1}}}} \\ {{v}_{{{z}_{2}}}} \\ \vdots \\ {{v}_{{{z}_{n}}}} \\ \end{matrix} \right] \\ \end{matrix} $
(3) 式中,v表示计算值与实测值之差。在满足${{\left\| \left( \Delta {{\mathrm{A}}_{1}},\Delta \mathrm{Z} \right) \right\|}_{F}}=\text{min}$的约束条件下(${{\left\| \text{ }\!\!\cdot\!\!\text{ } \right\|}_{F}} $为F范数),可求得ξ的值。
-
参照WTLS模型迭代解算方法,可将(2)式可改写成如下形式:
$ \mathit{\boldsymbol{Z}} - \mathit{\boldsymbol{e}} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\xi }}_{(i)}} + {\mathit{\boldsymbol{A}}_{(i)}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\xi }} - \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{\xi }}_{(i)}} $
(4) 式中,A(i)=A-E(i),δξ=ξ(i+1)-ξi,根据(4)式构造拉格朗日目标函数:
$ \begin{array}{l} \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(\mathit{\boldsymbol{e}}, \mathit{\boldsymbol{E}}, \mathit{\boldsymbol{\lambda }}, \mathit{\boldsymbol{\zeta }}) = {\mathit{\boldsymbol{e}}^{\rm T}}{\mathit{\boldsymbol{Q}}_Z}^{ - 1}e + {\mathit{\boldsymbol{E}}^{\rm T}}\mathit{\boldsymbol{Q}}_A^{ - 1}\mathit{\boldsymbol{E}} + 2{\mathit{\boldsymbol{\lambda }}^{\rm{T}}}\mathit{\boldsymbol{Z}} - \\ \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\xi }}_{(i)}} - {\mathit{\boldsymbol{A}}_{(i)}}\delta \mathit{\boldsymbol{\xi }} - \mathit{\boldsymbol{e}} + \left( {{\mathit{\boldsymbol{\xi }}_{(i)}}^{\rm{T}} \otimes \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{e}}] \end{array} $
(5) 式中,λ为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的列向量权阵P0及A的行向量初始权阵PA(1)、观测向量初始权阵PZ(1):
$ \left\{ {\begin{array}{*{20}{l}} {\mathop {{\mathit{\boldsymbol{P}}_0}}\limits_{3 \times 3} = {\rm{diag}}\left[ {\begin{array}{*{20}{c}} 1&1&0 \end{array}} \right]}\\ {\mathop {{\mathit{\boldsymbol{P}}_{A(1)}}}\limits_{n \times n} = \mathop {{\mathit{\boldsymbol{P}}_{Z(1)}}}\limits_{n \times n} = {\rm{diag}}\left[ {\overbrace {\begin{array}{*{20}{c}} 1&1& \ldots &2 \end{array}}^n} \right]} \end{array}} \right. $
(6) (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)}}}$:
$ \left\{ {\begin{array}{*{20}{L}} {{{\mathit{\boldsymbol{\hat v}}}_{{\rm{(0)}}}} = 0}\\ {{{\mathit{\boldsymbol{\hat \xi }}}_{{\rm{(0)}}}} = {{({\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{Z}}{\rm{(1)}}}}\mathit{\boldsymbol{A}})}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{Z}}{\rm{(1)}}}}\mathit{\boldsymbol{Z}}}\\ {{{\mathit{\boldsymbol{\hat \mu }}}_{{\rm{(0)}}}} = {{\left\{ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{Z}}{\rm{(1)}}}} + \left[ {{{({{\mathit{\boldsymbol{\hat \xi }}}_{{\rm{(0)}}}})}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_0}{{\mathit{\boldsymbol{\hat \xi }}}_{{\rm{(0)}}}}} \right]{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{A}}{\rm{(1)}}}}} \right\}}^{ - 1}}}\\ {{{\mathit{\boldsymbol{\hat \xi }}}_{{\rm{(1)}}}} = \left[ {{{({{\mathit{\boldsymbol{\hat \mu }}}_{{\rm{(0)}}}}\mathit{\boldsymbol{A}})}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{{\mathit{\boldsymbol{\hat \mu }}}_{{\rm{(0)}}}}} \right]\mathit{\boldsymbol{Z}}} \end{array}} \right. $
(7) 式中,QZ, QA, Q0分别为PZ, PA, P0的广义逆。
(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{)}}}} $:
$ \left\{ {\begin{array}{*{20}{L}} {{{\mathit{\boldsymbol{\hat \mu }}}_{{\rm{(}}\mathit{i}{\rm{)}}}} = {{\{ {\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{Z}}(i)}} + [{{({{\mathit{\boldsymbol{\hat \xi }}}_{(i)}})}^T}{\mathit{\boldsymbol{Q}}_0}{{\mathit{\boldsymbol{\hat \xi }}}_{(i)}}]{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{A}}(i)}}\} }^{-1}}}\\ {{{\mathit{\boldsymbol{\hat \lambda }}}_{(i)}} = {{\mathit{\boldsymbol{\hat \mu }}}_{{\rm{(}}\mathit{i}{\rm{)}}}}\cdot\left( {\mathit{\boldsymbol{Z}}\mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat \xi }}}_{(i)}}} \right)}\\ {{{\mathit{\boldsymbol{\hat v}}}_{{\rm{(}}\mathit{i}{\rm{)}}}} = {{\left( {{{\mathit{\boldsymbol{\hat \lambda }}}_{(i)}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{A}}(i)}}{{\mathit{\boldsymbol{\hat \lambda }}}_{(i)}}} \end{array}} \right. $
(8) (4) 计算${{\mathit{\boldsymbol{\hat \xi }}}_{1(i + 1)}}$,即$ {{\mathit{\boldsymbol{\hat \xi }}}_{(i + 1)}}$的含误差部分:
$ {{\mathit{\boldsymbol{\hat \xi }}}_{1(i + 1)}} = {({A^{\rm{T}}}{{\mathit{\boldsymbol{\hat \mu }}}_{{\rm{(}}\mathit{i}{\rm{)}}}}A - {{\mathit{\boldsymbol{\hat v}}}_{{\rm{(}}\mathit{i}{\rm{)}}}}\cdot{\mathit{\boldsymbol{Q}}_0})^{ - 1}}\cdot({A^{\rm{T}}}{{\mathit{\boldsymbol{\hat \mu }}}_{{\rm{(}}\mathit{i}{\rm{)}}}}Z) $
(9) (5)利用(10)式计算各残差V(i):
$ {\mathit{\boldsymbol{V}}_{(i)}} = \frac{{\mathit{\boldsymbol{Z-}}\mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat \xi }}}_{(i + 1)}}}}{{\sqrt {1 + {{({{\mathit{\boldsymbol{\hat \xi }}}_{1(i + 1)}})}^{\rm{T}}}{{\mathit{\boldsymbol{\hat \xi }}}_{1(i + 1)}}} }} $
(10) (6) 根据V(i)及IGGⅢ方案(见(11)式)计算权因子ω(i),令P(i)=P(i)ω(i),之后根据准则重新设定权阵PA(i+1), PZ(i+1):
$ {\rm{ }}{\omega _{(i)}} = \left\{ \begin{array}{l} 1, (\left| {\frac{{{\mathit{\boldsymbol{V}}_{(i)}}}}{\sigma }} \right| < {k_0})\\ \frac{{{k_0}}}{{\left| {\frac{{{\mathit{\boldsymbol{V}}_{(i)}}}}{\sigma }} \right|}}(\frac{{{k_1} - \left| {\frac{{{\mathit{\boldsymbol{V}}_{(i)}}}}{\sigma }} \right|}}{{{k_1}{k_0}}}), ({k_0} \le \left| {\frac{{{\mathit{\boldsymbol{V}}_{(i)}}}}{\sigma }} \right| < {k_1})\\ 0, (\left| {\frac{{{\mathit{\boldsymbol{V}}_{(i)}}}}{\sigma }} \right| > {k_1}) \end{array} \right. $
(11) 式中,V(i)为残差,σ为中误差,k0=1.0~1.5,k1=2.5~3.0,本文中取k0=1.5, k1=2.5。
(7) 计算权因子中误差σω(i),并剔除ω(i)>3σω(i)的点,同时根据下式计算$ {{\mathit{\boldsymbol{\hat \xi }}}_{2(i + 1)}}$,综合步骤(4)可得参量:
$ {{\mathit{\boldsymbol{\hat \xi }}}_{2(i + 1)}} = - {({\mathit{\boldsymbol{A}}_2}^{\rm{T}}{\mathit{\boldsymbol{A}}_2})^{ - 1}}({\mathit{\boldsymbol{A}}_2}^{\rm{T}}{\mathit{\boldsymbol{A}}_2} - {\mathit{\boldsymbol{A}}_2}^{\rm{T}}\mathit{\boldsymbol{Z}}) $
(12) (8) 重复步骤(3)~步骤(7),直到$\left| {{{\mathit{\boldsymbol{\hat \xi }}}_{(i + 1)}} - {{\mathit{\boldsymbol{\hat \xi }}}_{(i)}}} \right| < {\delta _0} $(δ0为给定阈值,本文中取δ0=10-6)。
经迭代运算得到拟合平面的参量解之后,可利用(13)式计算单位权中误差$ {{\hat \sigma }_0}$及平面拟合精度$ {{\hat \sigma }_p}$进行精度评定:
$ \left\{ {\begin{array}{*{20}{l}} {{{\hat \sigma }_0} = \sqrt {\frac{{{\mathit{\boldsymbol{\lambda }}_{(\mathit{i})}}^{\rm{T}}(\mathit{\boldsymbol{Z - A}}{{\mathit{\boldsymbol{\hat \xi }}}_{(i + 1)}})}}{{n - 3}}} }\\ {{{\hat \sigma }_{\rm{p}}} = \sqrt {\frac{{\sum\limits_{i = 1}^n {{d_{(\mathit{i})}}^2} }}{n}} } \end{array}} \right. $
(13) 式中,d(i)为点至拟合平面的距离。
-
为对本文中所构建平面拟合方法的适用性及优越性进行验证,分别利用其对模拟平面数据、实测平面数据进行拟合。
-
设要拟合的平面方程为Z=3X+4Y+5,依次抽取4组观测数据,并分别按照0%, 5%, 10%, 20%的比例加入粗差(大小为3σ0~5σ0,σ0为每组观测数据的中误差),每个观测数据中,X, Y分别为[0, 1000], [0, 2000]区间内的随机整数,相应的Z为[5 11005]区间内的随机整数,同时在模拟数据中加入随机误差e,且随机误差e服从均值为0、标准差为σ2I的正态分布(σ=0.3,I为单位阵)。
分别利用LS法、TLS法(奇异值分解(singular value deeoniposition, SVD))、LS-TLS法、IGGⅢ-LS-TLS法[5]、LS-WTLS混合法(本文中算法)对上述4组观测数据进行拟合实验,各方法计算得到的平面参量及精度评定指标如表 1所示。
Table 1. Plane parameters and fitting precision of simulated data
proportion of σ Δ|a| Δ|b| Δ|c| $ {{{\hat \sigma }_0}}$ $ {{{\hat \sigma }_p}}$ 0% LS 0 0 0.0007 0.0152 0.0029 TLS 0 0 0.0009 0.0152 0.0029 LS-TLS 0 0 0.0007 0.0030 0.0029 IGGⅢ-LS-TLS 0 0 0.0007 0.0024 0.0029 LS-WTLS 0 0 0.0003 0.0023 0.0022 5% LS 0.0020 0 0.0047 3.5046 0.6770 TLS 0.0158 0.0064 15.2664 6.4890 1.2569 LS-TLS 0.0020 0 0.0047 3.5046 0.6770 IGGⅢ-LS-TLS 0.0003 0.0002 0.3013 0.2115 0.6843 LS-WTLS 0.0001 0 0.0028 0.0565 0.0556 10% LS 0.0036 0.0008 1.7493 4.1745 0.8067 TLS 0.0193 0.0081 19.1327 7.4930 1.4523 LS-TLS 0.0035 0.0008 1.7493 4.1745 0.8067 IGGⅢ-LS-TLS 0 0.0001 0.7808 0.1841 0.8314 LS-WTLS 0 0.0001 0.0718 0.0823 0.0810 20% LS 0.0113 0.0020 2.6627 11.0755 2.1427 TLS 0.1209 0.0527 124.3950 44.9742 8.8819 LS-TLS 0.0111 0.0019 2.5019 2.1755 2.1427 IGGⅢ-LS-TLS 0.0016 0.0006 4.6832 0.94298 2.2148 LS-WTLS 0.0014 0.0005 0.4423 0.51715 0.5080 表 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所示。
Table 2. Observed data of plane
X Y Z 1 11.2 36.0 -5.0 2 10.0 40.0 -6.8 3 8.5 35.0 -4.0 4 8.0 48.0 -5.2 5 9.4 53.0 -6.4 6 8.4 23.0 -6.0 7 3.1 19.0 -7.1 8 10.6 34.0 -6.1 9 4.7 24.0 -5.4 10 11.7 65..0 -7.7 11 9.4 44.0 -8.1 12 10.1 31.0 -9.3 13 11.6 29.0 -9.3 14 12.6 58.0 -5.1 15 10.9 37.0 -7.6 16 23.1 46.0 -9.6 17 23.1 50.0 -7.7 18 21.6 44.0 -9.3 19 23.1 56.0 -9.5 20 19.0 36.0 -5.4 21 26.8 58.0 -16.8 22 21.9 51.0 -9.9 Table 3. Plane parameters and fitting precision of observed data
${\hat a} $ ${\hat b} $ ${\hat c} $ ${{{\hat \sigma }_0}} $0 ${{{\hat \sigma }_{\rm{p}}}} $ LS -0.2710 0.0085 -4.2789 2.1536 1.9316 TLS -0.2558 0.2616 -15.9565 4.1357 3.6093 LS-TLS -0.3108 0.0215 -4.2807 2.0680 1.9219 IGGⅢ-LS-TLS -0.2111 0.0279 -5.8982 1.4735 1.8997 LS-WTLS -0.2384 0.0388 -5.8379 1.4018 1.2924 由于实际平面数据的受测量仪器、外界环境、人为干扰等多种因素影响,因此数据纯度较低,其包含的偶然误差及粗差的比例存在不确定性。由表 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的稳健平面拟合方法
A robust methods of fitting plane based on 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%,表现出了良好的抗粗差干扰能力。该研究验证了新方法在平面拟合领域的优越性和可靠性。Abstract: In order to solve the positional accuracy difference of plane data and the low calculation precision of plane model's constant term, a robust plane fitting method based on least squares-weighted total least square (LS-WTLS) was proposed. This method uses least squares model and robust estimation of IGGⅢ scheme to calculate the error parameters of plane model. Meanwhile, after rejecting the gross error data by setting the threshold, the constant term of plane model was calculated by using least square model. And based on this model, the accuracy of plane parameters was further improved. The new method shows favorable resistant to gross errors in experiments of fit the simulated plane data, meanwhile, the observed plane data fitting experiments show that compared with LS method, TLS method, LS-TLS method, IGGⅢ-LS-TLS method, the new method's mean square error of unit weight increased by 53.6%, 195.0%, 47.5%, and 5.1%, respectively, and its plane fitting accuracy increased by 53.6%, 195.0%, 47.5%, and, 5.1%, respectively. The results effectively verify this new method's superiority and reliability in the field of plane fitting.
-
Table 1. Plane parameters and fitting precision of simulated data
proportion of σ Δ|a| Δ|b| Δ|c| $ {{{\hat \sigma }_0}}$ $ {{{\hat \sigma }_p}}$ 0% LS 0 0 0.0007 0.0152 0.0029 TLS 0 0 0.0009 0.0152 0.0029 LS-TLS 0 0 0.0007 0.0030 0.0029 IGGⅢ-LS-TLS 0 0 0.0007 0.0024 0.0029 LS-WTLS 0 0 0.0003 0.0023 0.0022 5% LS 0.0020 0 0.0047 3.5046 0.6770 TLS 0.0158 0.0064 15.2664 6.4890 1.2569 LS-TLS 0.0020 0 0.0047 3.5046 0.6770 IGGⅢ-LS-TLS 0.0003 0.0002 0.3013 0.2115 0.6843 LS-WTLS 0.0001 0 0.0028 0.0565 0.0556 10% LS 0.0036 0.0008 1.7493 4.1745 0.8067 TLS 0.0193 0.0081 19.1327 7.4930 1.4523 LS-TLS 0.0035 0.0008 1.7493 4.1745 0.8067 IGGⅢ-LS-TLS 0 0.0001 0.7808 0.1841 0.8314 LS-WTLS 0 0.0001 0.0718 0.0823 0.0810 20% LS 0.0113 0.0020 2.6627 11.0755 2.1427 TLS 0.1209 0.0527 124.3950 44.9742 8.8819 LS-TLS 0.0111 0.0019 2.5019 2.1755 2.1427 IGGⅢ-LS-TLS 0.0016 0.0006 4.6832 0.94298 2.2148 LS-WTLS 0.0014 0.0005 0.4423 0.51715 0.5080 Table 2. Observed data of plane
X Y Z 1 11.2 36.0 -5.0 2 10.0 40.0 -6.8 3 8.5 35.0 -4.0 4 8.0 48.0 -5.2 5 9.4 53.0 -6.4 6 8.4 23.0 -6.0 7 3.1 19.0 -7.1 8 10.6 34.0 -6.1 9 4.7 24.0 -5.4 10 11.7 65..0 -7.7 11 9.4 44.0 -8.1 12 10.1 31.0 -9.3 13 11.6 29.0 -9.3 14 12.6 58.0 -5.1 15 10.9 37.0 -7.6 16 23.1 46.0 -9.6 17 23.1 50.0 -7.7 18 21.6 44.0 -9.3 19 23.1 56.0 -9.5 20 19.0 36.0 -5.4 21 26.8 58.0 -16.8 22 21.9 51.0 -9.9 Table 3. Plane parameters and fitting precision of observed data
${\hat a} $ ${\hat b} $ ${\hat c} $ ${{{\hat \sigma }_0}} $0 ${{{\hat \sigma }_{\rm{p}}}} $ LS -0.2710 0.0085 -4.2789 2.1536 1.9316 TLS -0.2558 0.2616 -15.9565 4.1357 3.6093 LS-TLS -0.3108 0.0215 -4.2807 2.0680 1.9219 IGGⅢ-LS-TLS -0.2111 0.0279 -5.8982 1.4735 1.8997 LS-WTLS -0.2384 0.0388 -5.8379 1.4018 1.2924 -