留言板

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

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

GPS高程转换的总体最小二乘拟合推估模型

王乐洋 吴飞 吴良才

王乐洋, 吴飞, 吴良才. GPS高程转换的总体最小二乘拟合推估模型[J]. 武汉大学学报 ● 信息科学版, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
引用本文: 王乐洋, 吴飞, 吴良才. GPS高程转换的总体最小二乘拟合推估模型[J]. 武汉大学学报 ● 信息科学版, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
WANG Leyang, WU Fei, WU Liangcai. Total Least Squares Fitting Estimation Model for GPS Height Transformation[J]. Geomatics and Information Science of Wuhan University, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
Citation: WANG Leyang, WU Fei, WU Liangcai. Total Least Squares Fitting Estimation Model for GPS Height Transformation[J]. Geomatics and Information Science of Wuhan University, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421

GPS高程转换的总体最小二乘拟合推估模型

doi: 10.13203/j.whugis20140421
基金项目: 

国家自然科学基金 41204003

江西省教育厅科技项目 GJJ150595

江西省教育厅科技项目 KJLD12077

江西省教育厅科技项目 KJLD14049

流域生态与地理环境监测国家测绘地理信息局重点实验室项目 WE2015005

对地观测技术国家测绘地理信息局重点实验室项目 K201502

测绘地理信息公益性行业科研专项 201512026

东华理工大学博士科研启动基金 DHBK201113

国家重点研发计划 2016YFB0501405

详细信息
    作者简介:

    王乐洋, 博士, 副教授, 主要研究方向为大地测量反演及大地测量数据处理。wleyang@163.com

  • 中图分类号: P228.41

Total Least Squares Fitting Estimation Model for GPS Height Transformation

Funds: 

The National Natural Science Foundation of China 41204003

Science and Technology Project of the Education Department of Jiangxi Province GJJ150595

Science and Technology Project of the Education Department of Jiangxi Province KJLD12077

Science and Technology Project of the Education Department of Jiangxi Province KJLD14049

Key Laboratory of Watershed Ecology and Geographical Environment Monitoring, NASG WE2015005

Key Laboratory of Mapping from Space, NASG K201502

National Department Public Benefit Research Foundation (Surveying, Mapping and Geoinformation) 201512026

Scientific Research Foundation of ECIT DHBK201113

National Key Research and Development Program 2016YFB0501405

More Information
    Author Bio:

    WANG Leyang, PhD, associate professor, specializes in geodetic inversion and geodetic data processing.wleyang@163.com

图(1) / 表(4)
计量
  • 文章访问数:  1070
  • HTML全文浏览量:  33
  • PDF下载量:  418
  • 被引次数: 0
出版历程
  • 收稿日期:  2014-09-09
  • 刊出日期:  2016-09-05

GPS高程转换的总体最小二乘拟合推估模型

doi: 10.13203/j.whugis20140421
    基金项目:

    国家自然科学基金 41204003

    江西省教育厅科技项目 GJJ150595

    江西省教育厅科技项目 KJLD12077

    江西省教育厅科技项目 KJLD14049

    流域生态与地理环境监测国家测绘地理信息局重点实验室项目 WE2015005

    对地观测技术国家测绘地理信息局重点实验室项目 K201502

    测绘地理信息公益性行业科研专项 201512026

    东华理工大学博士科研启动基金 DHBK201113

    国家重点研发计划 2016YFB0501405

    作者简介:

    王乐洋, 博士, 副教授, 主要研究方向为大地测量反演及大地测量数据处理。wleyang@163.com

  • 中图分类号: P228.41

摘要: 在现有的关于GPS高程转换的总体最小二乘方法研究中,通常是将高程异常转换参数的计算与待求点高程异常的计算分两步进行处理,并且只考虑由已知高程异常点的平面坐标组成的系数矩阵的误差,忽略了高程异常待求点的坐标误差。针对以上问题,本文提出了GPS高程转换的总体最小二乘拟合推估模型,将计算高程异常转换参数和待求点高程异常联合处理,且考虑到所有点的点位误差,最后采用拟合推估法进行求解。实验结果表明,本文方法能够有效地提高高程转换的精度。

English Abstract

王乐洋, 吴飞, 吴良才. GPS高程转换的总体最小二乘拟合推估模型[J]. 武汉大学学报 ● 信息科学版, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
引用本文: 王乐洋, 吴飞, 吴良才. GPS高程转换的总体最小二乘拟合推估模型[J]. 武汉大学学报 ● 信息科学版, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
WANG Leyang, WU Fei, WU Liangcai. Total Least Squares Fitting Estimation Model for GPS Height Transformation[J]. Geomatics and Information Science of Wuhan University, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
Citation: WANG Leyang, WU Fei, WU Liangcai. Total Least Squares Fitting Estimation Model for GPS Height Transformation[J]. Geomatics and Information Science of Wuhan University, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
  • 由于GPS可以在同一时间获得点位的三维坐标, 测量工作的效率得到了很大的提高, 但是GPS获得的是基于WGS84参考椭球的大地高, 在实际工程应用中, 我国高程系统普遍采用正常高, 如果大地高能够准确地转换为正常高, 现在的测量工作所需要的平面控制测量和高程控制测量就能够同时进行。

    对于大多数工程应用, 由于测区范围不大, 高程异常变化平缓, 多项式拟合便可满足精度要求[1], 可通过建立GPS高程异常值与平面坐标之间的多项式关系求解转换参数。在实际参数估计问题中, 对于大部分线性模型来说, 一般由于数据采样大小、模型化及测量等原因引起了系数矩阵的误差, 故其系数阵并非常数阵而是由其他方法计算得到的具有一定误差的变量矩阵[2], 通常可以采用总体最小二乘来进行解算。目前已有相关学者采用总体最小二乘方法进行GPS高程转换参数的估计:文献[3]考虑不同曲面拟合模型中的水平和竖直坐标中含有误差的情形, 分别基于二次多项式曲面、三次多项式曲面、移动二次多项式曲面, 利用最小二乘估计、总体最小二乘和加权总体最小二乘算法进行高程转换; 文献[4]引入基于加权总体最小二乘的二次多项式GPS高程转换, 并采用了一种简化的定权方法; 文献[5]考虑到二次多项式GPS高程转换模型系数矩阵的特点, 采用混合总体最小二乘的方法进行GPS高程转换, 在合理考虑系数矩阵和观测向量误差的情况下, 利用其提出的基于总体最小二乘平差的粗差探测方法进行粗差剔除后, 对GPS高程转换参数进行求解; 文献[6]考虑系数矩阵和观测数据的随机误差, 讨论基于总体最小二乘法建立的高斯-马尔科夫模型, 同时考虑由于粗差的存在, 观测数据精度不相等的情况, 提出稳健总体最小二乘的算法来计算转换参数。

    但是, 在现有的关于总体最小二乘方法在GPS高程转换中应用的研究中, 通常都是将转换参数的计算与待求点高程异常计算独立进行, 并且在对高程异常待求点进行高程转换时忽略了待求点本身的坐标误差。因此, 本文在文献[7]的基础上提出GPS高程转换的总体最小二乘拟合推估模型, 该模型考虑到所有点的点位误差, 且将计算转换参数和计算待求点高程异常联合处理, 最后采用拟合推估法进行求解, 并通过模拟算例对该方法的可行性进行了验证。

    • 基于二次多项式曲面的GPS高程转换模型为:

      (1)

      式中, ξi为坐标xi, yi的高程异常; λ=[a0a1a2a3a4a5]T为待估参数。式(1)的矩阵形式可写为:

      (2)

      其中, Ai=[1xiyixi2xiyiyi2]。

      由于控制点的坐标值及高程异常不可避免地受观测误差的影响, 对观测向量及系数矩阵引入误差向量可得:

      (3)

      其中, EAi=[0ExiEyiExi2ExiyiEyi2]。

      联合n个高程异常已知点求解高程异常转换参数, 则:

      (4)

      式中,

      假设有m个高程异常待求点, 计算它们的高程异常:

      (5)

      式中, 系数矩阵B=由高程异常待求点的坐标组成; 为对应的系数矩阵误差阵; g为高程异常待求点的转换值。从而建立GPS高程转换联合模型:

      (6)

      假设控制点的点位坐标误差相关, 且服从零均值正态分布, 则联合模型的各项误差的随机模型为:

      (7)

      式中, vec(·)为矩阵向量化算子; Q是协因数阵。

    • 拟合推估模型的标准形式为[7]:

      (8)

      式中, yey分别为n×1观测向量和误差向量; βk×1的参数向量; C是它的设计矩阵且列满秩; y0n0×1未知随机信号向量; ey0C0分别为它的误差向量和设计矩阵。通常要求n>k, 但并不要求n0 > k。模型(8)的求解包含平差和预报两个部分, 平差是指采用观测向量y估计参数β, 预报则是指根据信号y0与参数β的函数关系及其与观测向量y的随机关系对其做出合理的推估。若eyey0服从正态分布:

      (9)

      根据最小二乘准则:

      (10)

      求解的最优线性无偏参数估计值和预报值分别为:

      (11)
      (12)
    • GPS高程转换联合模型是非线性模型, 本文将根据高斯-牛顿法[7]和转换参数的数值特点, 分别将式(6)转化为线性模型, 并采用拟合推估法进行求解。

    • 采用高斯-牛顿法将联合模型线性化后进行迭代计算, 设第j次迭代后的估值为λ(j), 则参数λ表示为:

      (13)

      式中, δλ(j)为第j+1次迭代的参数改正数。再将式(13)代回式(6)可得:

      (14)

      尽管EAδλ(j)EBδλ(j)是二阶的小量, 但是完全忽略二阶以上的小量, 可能导致迭代发散或收敛错误。采用EAEB的第j次迭代估值EA(j)EB(j)代替EAEB计算系数矩阵可避免该缺陷。记A(j)=A-EA(j), B(j)=B-EB(j), 则(14)可变为:

      (15)

      将上式改写成拟合推估模型的标准形式:

      (16)

      式中, InIm分别表示n维和m维单位阵; 为Kronecker积算子, 推导利用了关系式vecABC=(CTA)vec(B)。对应推估模型的标准形式(8), 可得各个变量为:

      (17)

      顾及式(17) eyey0及式(7), 根据协方差传播律推出Qyy, Qy0y分别为:

      (18)
      (19)

      根据上述中的无偏参数和预报值的公式可以推出:

      (20)
      (21)
      (22)
      (23)

      其中, 式(20)、(21)式用于估计高程异常转换参数, 式(22)、(23)用于计算待求点的高程异常值。由于预报和估计是相对独立的过程, 因此只需要按照式(20)、(21)计算转换参数, 迭代终止后, 再按式(22)、(23)计算最后的推估值。

    • 求解两种模型需要已知QAAQBA, 定义变量分别为n个已知高程异常点和m个未知高程异常点的平面坐标误差。其对应的协方差阵为Q=

      由系数矩阵的表达式, 结合误差传播定律可以得到:

      (24)

      并可推出

      n个点系数矩阵误差可以写成:

      (25)

      分别令: 。其中, 可得到EA=J×Ez1, EB=K×Ez2

      根据误差传播定律, 可以得到:

      (26)

      式中, 为已知高程异常点和待求高程异常点的平面坐标之间的的协方差阵, 在实际生产中, 一般可以通过对所有点进行测边网平差或者GPS基线向量解算来获得, 在本文的算例中, 通过对模拟的测边网进行平差计算得到。

    • 选取14个点, 将其中的A、B两点作为起始坐标点, 模拟距离观测值, 建立一个测边网, 计算其他12个点的坐标平差值及其协方差阵, 选取2、6、8、9号点作为检验点, 各点的坐标真值如表 1所示。图 1为这14个点的空间分布, 可以看出所有数据点的空间分布较为均匀。令(b0, b1, b2, b3, b4, b5)=(-11, 0.05, 0.1, 0.01, 0.009, -0.02)为二次多项式曲面高程异常转换参数的真值[6], 并根据坐标真值及参数真值计算高程异常的真值。

      图  1  测边网点的空间分布

      Figure 1.  Spatial Distribution of Trilateration Net's Point

      表 1  各点坐标、高程异常/m

      Table 1.  Coordinates and Height Anomaly of Points/m

      点号 X Y ξ
      A 100 100 -
      B 120 120 -
      1 120 140 -87.8
      3 160 120 149.8
      4 170 200 -187.5
      5 150 190 -225
      7 130 220 -524.1
      10 110 180 -336.3
      11 170 160 35.3
      12 200 170 144
      2 140 160 -102.4
      6 190 140 220.9
      8 180 100 294
      9 200 220 -151

      反算两点间的边长作为距离的真值S0, 并对S0添加服从正态分布N(0, σ2)的随机误差作为距离的观测值, 其中, σ=a+b×S, 本文中计算取常数a=5 mm, b=0.05 mm/km。联合所有距离误差方程, 权取距离方差的倒数, 计算得到这12个点的坐标平差值及其协方差阵。再给高程异常的真值添加服从正态分布N0, σξ2随机误差, 作为高程异常的观测值, 其中σξ2为测边网平差过程中计算得到的单位权方差, 得到的平差坐标值及含有误差的高程异常值如表 2所示。

      表 2  各点平差坐标、加入误差的高程异常/m

      Table 2.  Adjusted Coordinates and Add the Abnormal Height of Error/m

      点号 X Y ξ
      A 100 100 -
      B 120 120 -
      1 119.997 9 140.000 6 -87.798 0
      3 159.989 1 120.001 7 149.795 9
      4 170.011 9 199.980 9 -187.495 4
      5 150.013 1 189.987 9 -224.997 6
      7 130.016 7 219.995 4 -524.100 7
      10 110.008 0 180.000 5 -336.299 9
      11 170.006 2 159.987 1 35.300 2
      12 200.013 5 169.973 9 144.003 0
      2 140.003 4 159.996 6 -102.396 2
      6 190.003 3 139.983 0 220.899 8
      8 179.962 3 100.003 3 293.999 2
      9 200.005 6 219.983 8 -151.001 1

      本文中采用的方案如表 3所示, 其中方案5所采用的定权方法, 是根据协因数传播定律求得协因数之间的相互关系, 并考虑到协因数会随点的不同而变化, 为了方便起见, 对定权的方法做了一定的简化[4]; 方案4则引用了文献[8]中采用的加权总体最小二乘迭代计算的方法, 该方法直接采用系数矩阵的协方差阵进行计算。

      表 3  各方案所用方法列表

      Table 3.  List of Methods Used in Different Schemes

      方案 方法
      方案1 本文提出方法
      方案2 最小二乘法
      方案3 总体最小二乘法
      方案4 加权总体最小二乘法[8]
      方案5 加权总体最小二乘法[4]

      各个方案的计算结果如表 4所示, 从计算结果的比较中可以看出, 采用传统最小二乘方法和总体最小二乘方法的计算结果差别并不大; 两种加权总体最小二乘法解算精度较好, 方案5采用的是简化的定权方法, 而方案4根据协因数阵进行计算, 方案4的计算方法相对方案5而言更加合理, 且方案4的转换效果也优于方案5。但是, 不论采用何种加权总体最小二乘方法, 仍然存在着目前关于总体最小二乘方法在GPS高程转换的研究中普遍存在的不足, 即将转换参数的计算与待求点高程异常计算分两步进行处理, 并且只考虑由高程异常已知点的平面坐标组成的系数矩阵的误差, 忽略了高程异常待求点的坐标误差。方案1将计算转换参数和计算检验点高程异常联合处理, 且考虑到所有点的点位误差, 从理论上来说该方案要更为严谨与合理, 从计算的结果可以看出, 该方案相对其他4种方案, 转换精度有所提高。

      表 4  检验点转换的结果及与真值的残差范数/m

      Table 4.  Checkpoint's Result of the Transformation and Norm of Residual Between the Estimated Value and True Value/m

      方案 点号
      2 6 8 9 ‖Δ
      真值 -102.4 220.9 294 -151 -
      方案1 -102.394 2 220.901 8 293.896 6 -151.149 4 0.181 7
      方案2 -102.413 5 220.914 5 293.929 7 -151.186 5 0.200 3
      方案3 -102.413 5 220.914 5 293.929 7 -151.186 5 0.200 3
      方案4 -102.405 1 220.893 1 293.900 7 -151.154 6 0.183 9
      方案5 -102.413 5 220.914 5 293.929 7 -151.186 5 0.199 3
    • 本文考虑到由高程异常待求点组成的系数矩阵中的误差, 并且将计算高程异常转换参数与计算待求点的高程异常联合处理, 提出了GPS高程转换拟合推估模型, 结合模拟数据进行转换参数的计算以及解算结果精度的评定。模拟算例的计算结果表明, 本文方法能够得到更高精度的高程转换结果。

参考文献 (15)

目录

    /

    返回文章
    返回