高级检索

ISSN1001-3806CN51-1125/TN 网站地图

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

一种基于LS-WTLS的稳健平面拟合方法

欧江霞 邓雄文 蔡茂欣 邱敏

引用本文:
Citation:

一种基于LS-WTLS的稳健平面拟合方法

    作者简介: 欧江霞(1989-),男,硕士,注册测绘师,主要从事大地测量数据处理方法的研究。E-mail:oujiangxia666@163.com.
  • 基金项目:

    国家自然科学基金资助项目 41974214

  • 中图分类号: P207

A robust methods of fitting plane based on LS-WTLS

  • CLC number: P207

  • 摘要: 为了解决平面数据点位精度差异性及平面模型常数项解算精度较低的问题,提出了一种基于最小二乘-加权总体最小二乘(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%,表现出了良好的抗粗差干扰能力。该研究验证了新方法在平面拟合领域的优越性和可靠性。
  • 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
    下载: 导出CSV

    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
    下载: 导出CSV

    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
    下载: 导出CSV
  • [1]

    WU Ch, YUAN Y B, ZHANG M Y. Plane target positioning based on reflection intensity and k-means clustering method.Laser Technology, 2015, 39(3): 341-343(in Chinese).
    [2]

    OU J X, LIU W Ch, CAI M X. Study on influence factors of scanning quality for leica scanstation C10 terrestrial laser scanner [J]. Geospatial Information, 2017, 15(6): 103-106(in Chinese).
    [3]

    BURKHARD S, YARON A F. On the multivariate total least-squares approach to empirical coordinate transformations [J]. Journal of Geodesy, 2008, 82(7): 373-383.
    [4]

    LU T D, ZHOU Sh J. An iterative algorithm for total least squares sstimation [J]. Geomatics an Information Science of Wuhan University, 2010, 35(11):1351-1354(in Chinese).
    [5]

    GUAN Y L, LIU Sh T, ZHOU Sh J. Robust plane fitting of point clouds based on TLS.Journal of Geodesy and Geodynamics, 2011, 31(5):80-83(in Chinese).
    [6]

    GONG X Q. A roubust weighted total least squares method [J].Acta Geodaetica et Cartographica Sinica, 2014, 47(10):462-466(in Chinese).
    [7]

    BURKHARD S, YARON W. On weighted total least-squars adjustment for linear regression[J]. Journal of Geodesy, 2008, 82(6): 415-421.
    [8]

    SHEN Y Z, LI B F, FELUS Y A. An iterative solution of weighted total least-squares adjustment [J]. Journal of Geodesy, 2010, 85(4):229-238.
    [9]

    LI M F, OU J X, TAN D. Study on fixed weight methods in plane fitting of point clouds based on weighted total least squares [J]. Journal of Geodesy and Geodynamics, 2015, 35(6): 428-432(in Chinese).
    [10]

    OU J X, LI M F, WANG Y M. Plane fitting of point clouds based on robust weighted total least squares [J]. Journal of Geodesy and Geodynamics, 2014, 34(6): 160-163(in Chinese).
    [11]

    OU J X, LIU W Ch. Fitting of sphere point clouds by weighted total least squares based on IGGⅢ Scheme.Laser Technology, 2017, 41(5): 749-753(in Chinese).
    [12]

    GONG X Q, LI Zh L. A roubust mixed LS-TLS based on IGGⅡ scheme [J]. Geomatics an Information Science of Wuhan University, 2014, 39(4): 462-466(in Chinese).
    [13]

    WANG X Zh, TAO B Z, QIU W N. Advanced surveying adjustment [M]. Beijing: Survey and Mapping Press, 2006:73-89(in Chinese).
    [14]

    WANG H Y, ZHANG Y. Data processing of laser rangefinder based on robust estimation [J].Journal of Atmospheric and Environmental Optics, 2007, 2(3): 228-231(in Chinese).
    [15]

    WU J F, YANG Y X. Robust estimation for correlated GPS baseline vector network [J].Acta Geodaetica et Cartographica Sinica, 2001, 30(3): 247-251(in Chinese).
    [16]

    QIU W N, TAO B Z, YAO Y B. et al. The theory and method surveying data processing[M]. Wuhan: Wuhan University Press, 2008:54-64(in Chinese).
  • [1] 欧江霞刘伟诚 . 基于IGGⅢ方案的加权总体最小二乘点云球面拟合. 激光技术, 2017, 41(5): 749-753. doi: 10.7510/jgjs.issn.1001-3806.2017.05.026
    [2] 苍桂华岳建平 . 基于加权总体最小二乘法的点云平面拟合. 激光技术, 2014, 38(3): 307-310. doi: 10.7510/jgjs.issn.1001-3806.2014.03.005
    [3] 鲍鸿曾海涛白玉磊胡忠向志聪周延周申作春 . 基于概率密度最小二乘拟合的叶片后缘轮廓. 激光技术, 2016, 40(4): 555-559. doi: 10.7510/jgjs.issn.1001-3806.2016.04.021
    [4] 肖兴维马国鹭曾国英陆野 . 正交视觉与倾角仪组合空间位姿测量方法研究. 激光技术, 2020, 44(3): 278-282. doi: 10.7510/jgjs.issn.1001-3806.2020.03.002
    [5] 刘顺涛骆华芬陈雪梅徐静 . 结构光测量系统的标定方法综述. 激光技术, 2015, 39(2): 252-258. doi: 10.7510/jgjs.issn.1001-3806.2015.02.023
    [6] 胡林亭史德民李佩军任成才 . 激光监测系统测量精度的检测方法. 激光技术, 2008, 32(6): 670-672.
    [7] 郑春艳杨若夫刘艺吴健 . 液晶光栅相控阵偏转波前测量方法研究. 激光技术, 2012, 36(5): 645-648. doi: 10.3969/j.issn.1001-3806.2012.05.018
    [8] 柳静李明詹高伟肖武华韦庆玥 . 三坐标激光测量技术规范中参量确定的方法. 激光技术, 2015, 39(1): 140-144. doi: 10.7510/jgjs.issn.1001-3806.2015.01.028
    [9] 黄东杨凌辉罗文张晓日史慎东黄喆王姣 . 基于视觉/惯导的掘进机实时位姿测量方法研究. 激光技术, 2017, 41(1): 19-23. doi: 10.7510/jgjs.issn.1001-3806.2017.01.005
    [10] 盛旭波杨晖李然马生郑刚ZIVKOVIC Vladimir . 基于激光在线测量滚筒内颗粒流崩塌角新方法. 激光技术, 2016, 40(3): 344-348. doi: 10.7510/jgjs.issn.1001-3806.2016.03.009
    [11] 李儒颂马红梅叶文江 . 基于迈克尔逊干涉液晶双折射率的测量方法设计. 激光技术, 2016, 40(4): 487-490. doi: 10.7510/jgjs.issn.1001-3806.2016.04.007
    [12] 孙卿杨凌辉 . 基于测距传感器的WMPS非接触测量方法. 激光技术, 2016, 40(5): 670-675. doi: 10.7510/jgjs.issn.1001-3806.2016.05.011
    [13] 徐达何凯平熊伟李华高源 . 线激光散斑检测弹幕武器炮口振动测量方法. 激光技术, 2017, 41(6): 876-880. doi: 10.7510/jgjs.issn.1001-3806.2017.06.022
    [14] 杨志远卢荣军王生春 . 一种纵向共振光声池谐振频率测量方法. 激光技术, 2019, 43(3): 387-391. doi: 10.7510/jgjs.issn.1001-3806.2019.03.018
    [15] 王以忠安忠猛黄喆张晓日相健张凯伦 . 基于空间矢量约束的煤矿掘进机组合测量方法. 激光技术, 2019, 43(6): 804-808. doi: 10.7510/jgjs.issn.1001-3806.2019.06.014
    [16] 马晓春董俊良梁芳马宁 . 一种基于菲涅耳原理的光纤盐度测量方法. 激光技术, 2010, 34(3): 313-315. doi: 10.3969/j.issn.1001-3806.2010.03.008
    [17] 叶景峰胡志云张振荣王晟刘晶儒 . 一种提高羟基分子示踪速率测量精度的方法. 激光技术, 2012, 36(1): 64-66. doi: 10.3969/j.issn.1001-3806.2012.01.017
    [18] 杨宇郝晓剑武耀艳周汉昌 . 基于热电偶动态校准的非线性拟合方法研究. 激光技术, 2014, 38(2): 145-148. doi: 10.7510/jgjs.issn.1001-3806.2014.02.001
    [19] 谢苏隆钟鹰 . 泽尔尼克多项式拟合波面中采样点数目的研究. 激光技术, 2011, 35(2): 272-274,277. doi: 10.3969/j.issn.1001-3806.2011.02.035
    [20] 付华代巍 . 随机共振瓦斯微弱信号检测方法研究. 激光技术, 2016, 40(2): 213-218. doi: 10.7510/jgjs.issn.1001-3806.2016.02.013
  • 加载中
计量
  • 文章访问数:  180
  • HTML全文浏览量:  79
  • PDF下载量:  0
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-11-14
  • 录用日期:  2020-01-17
  • 刊出日期:  2020-11-25

一种基于LS-WTLS的稳健平面拟合方法

    作者简介: 欧江霞(1989-),男,硕士,注册测绘师,主要从事大地测量数据处理方法的研究。E-mail:oujiangxia666@163.com
  • 广州市地质调查院,广州 510440
基金项目:  国家自然科学基金资助项目 41974214

摘要: 为了解决平面数据点位精度差异性及平面模型常数项解算精度较低的问题,提出了一种基于最小二乘-加权总体最小二乘(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%,表现出了良好的抗粗差干扰能力。该研究验证了新方法在平面拟合领域的优越性和可靠性。

English Abstract

    • 平面拟合可用于建(构)筑物倾斜监测、边坡变形监测、隧道断面监测、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法求取常数项,有效抵抗点位拟合权值对常数项解算的影响。

    • 设平面方程为:

      $ {\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的列向量权阵P0A的行向量初始权阵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平面拟合方法。在误差项参量解算过程中,顾及平面数据各点精度存在差异,本文中依据稳健估计中各权函数特点及其适用范围,选用IGGⅢ方案对加权总体最小二乘平面拟合方法进行改进,其可在模型参量解算过程中,通过计算点与平面模型的相关关系,自适应地调整各点拟合权值,优化拟合权阵;对于常数项,先以三倍权因子中误差为阈值剔除观测数据中的粗差,再利用LS法进行解算。模拟平面数据及实测平面数据的拟合实例表明,该方法具备较好的可行性及优越性,利用该方法拟合得到的平面参量解可靠性更高。

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

参考文献 (16)

目录

    /

    返回文章
    返回