基于GPS高程拟合的总体最小二乘算法研究

基于GPS高程拟合的总体最小二乘算法研究
崔立鲁;杜艺婕;江雪梨;袁梼;杜石
【摘要】针对GPS测量噪声影响高程拟合精度的问题,本文详细阐述了最小二乘方法、总体最小二乘算法和加权总体最小二乘算法三种误差处理方法的基本原理及计算公式.根据常用的两种曲面拟合模型,通过对实测数据拟合结果分析GPS测量噪声对高程拟合精度的影响,并对比上述三种算法的结果.数值计算结果表明在顾及GPS 测量噪声的情况下,总体最小二乘算法能够很好地削弱其对高程拟合的影响,从而提高拟合精度.
【期刊名称】《北京测绘》
【年(卷),期】2019(033)002
【总页数】4页(P144-147)
【关键词】总体最小二乘;高程拟合;测量噪声;曲面拟合
【作者】崔立鲁;杜艺婕;江雪梨;袁梼;杜石
【作者单位】成都大学建筑与土木工程学院,四川成都610106;武汉大学测绘学院,湖北武汉430079;成都大学建筑与土木工程学院,四川成都610106;成都大学建筑与土木工程学院,四川成都610106;成都大学建筑与土木工程学院,四川成都610106;成都大学建筑与土木工程学院,四川成都610106
【正文语种】中文
【中图分类】P227
0 引言
利用全球定位系统(Global Positioning System, GPS)技术进行高程测量存在着高程转换问题,即将GPS测量得到的大地高转换为常用的水准高[1]。

常用的方法是建立某一点高程异常和平面坐标(x,y)之间的数学关系,然后根据已知点数据解算数学模型中的未知参数,建立适应该测区的高程转换模型[2-3]。

最小二乘方法是经典的测量数据处理算法,常常被用于求解高程转换模型的未知参数,但已知点的坐标数据存在着由外部环境或观测方法引起的测量噪声[4]。

而最小二乘方法将由已知点的平面坐标数据组成的系数矩阵定义为常数阵,没有考虑到GPS测量噪声对计算结果的影响,因此对参数拟合精度造成了一定的影响[5-6]。

鉴于测量噪声对参数拟合精度的影响,本文分别介绍了最小二乘方法(Least Squares, LS)、总体最小二乘方法(Total Least Squares, TLS)和加权总体最小二乘方法(Weighted Least Squares, WLS)数学模型,然后利用实测数据分别对基于二次多项式曲面和三次多项式曲面的拟合模型进行数值计算,最后将三种算法的计算结果进行比较和分析。

1 GPS高程转换数据模型
1.1 多项式曲面拟合模型
二次多项式曲面拟合方程为:
ζ=a0+a1x+a2y+a3x2+a4xy+a5y2
(1)
式中,ζ为某点的高程异常,(x,y)为该点的平面坐标,ai(i=0,…,5)为待求拟合参数。

将其转换为矩阵形式,即
ζ=Aα
(2)
其中,
,
α,
ζ
三次多项式曲面拟合模型为:
ζ=a0+a1x+a2y+a3x2+a4xy+a5y2+
a6x2y+a7xy2+a8x3+a9y3
(3)
其矩阵形式同理可得。

1.2 LS法
以二次多项式曲面拟合方程为例,根据式(2)和LS法基本原理,可得参数计算公式如下:
ζ
(4)
从式(4)可知,矩阵A是由高程点的平面坐标组成的,在这里把平面坐标值作为真值。

但是在实际工程中,平面坐标值是通过测量得到的,因此本身是存在着观测噪声的。

忽略此观测噪声必然会对拟合结果造成影响。

因此,LS法在理论上并不严密。

1.3 TLS法
TLS方法的数学模型为[7-9]:
ζ=(A-EA)α+eζ
(5)
式中,α为待求参数向量,ζ为高程异常向量,A为系数矩阵,eζ为高程异常的误差向量,EA为系数矩阵的误差矩阵。

求解式(5)中的待估参数需要通过迭代法实现,其计
算步骤如下:
步骤 1 将LS法计算得到的参数值作为计算初值:
ζ
(6)
步骤 2 计算迭代过程中的中间变量:
ζζ
(7)
步骤 3 计算新的待求参数值:
ζ
(8)
判断迭代终止的阈值条件为:ε。

如果满足阈值要求,则为最终结果,否则回到步骤2,重新计算直到满足要求为止。

1.4 WTLS法
考虑到各个已知点到待求点距离不同,其对高程转换参数计算的贡献也不同。

因此在定权时将权阵P简单处理为单位阵是不合理的,需要根据实际情况进行定权。

设测区内共有n个已知点,第i个点到待求点的距离为di,相应的权为,则权阵为[8-9]: Pζ
(9)
设P0=diag(1,1,1,1,1,1)和Px=In×n分别表示系数矩阵A的列向量的权阵和行向量之间的权阵,即两者对参数计算的贡献是等同的。

再假设,,Qζ分别为权阵的逆矩阵,则WLS法的求解步骤为:
步骤 1 计算迭代的初值:
,ζ
(10)
步骤 2 计算中间变量:
α(i))Qx)-1(ζ
(11)
(12)
步骤 3 计算新的参数值:
AT(Qζζ
(13)
判断迭代终止的阈值条件为:ε。

如果满足阈值要求,则为最终结果,否则回到步骤2,重新计算直到满足要求为止[10-14]。

2 数值计算与分析
本文采用某地实测高程数据,具体情况如表1所示。

利用前12个点作为已知数据,后6个作为待求点。

分别根据二次多项式曲面拟合模型和三次多项式曲面拟合模型对待估参数进行拟合求解,最后求出待求点的高程异常值,并将拟合结果与真值进行对比分析。

已知点和待求点点位分布如图1所示。

图1中对所有点的平面坐标进行了标准化处理。

表1 实测点位高程数据单位:m点号XYζ146.049189.6995-
3.9739240.857699.6299-3.5376338.484192.1108-3.7673717.003298.8406-
3.3261822.1838101.9830-3.2858919.0450106.2452-
3.05611046.5016105.1646-3.39981332.3969110.7557-
3.00361527.9913123.9794-2.42622343.8827143.2226-
1.94002623.8734138.8852-1.73962725.5653145.7765-
1.52241436.3441123.2001-
2.58761934.5588126.6291-
2.40682034.9654131.8336-2.21042130.5073129.2705-
2.23222238.5731134.9822-2.13792531.7270137.7843-1.8762
表2 二次多项式曲面拟合结果点号252120191422真实值/m-1.8762-2.2322-2.2104-2.4068-2.5876-2.1379LS估计值/m-1.9333-2.2467-2.2167-2.4116-2.5703-2.1553误差/mm-57.10-14.53-6.33-4.8017.33-17.42中误差
/mm26.26TLS估计值/m-1.9273-2.2315-2.2048-2.3950-2.552-2.1504误差/mm-51.100.735.5611.8335.17-12.48中误差/mm26.18WTLS估计值/m-
1.9333-
2.2467-2.2167-2.4116-2.5703-2.1553误差/mm-57.10-14.53-6.33-4.8017.33-17.42中误差/mm26.26
表3 三次多项式曲面拟合结果点号252120191422真实值/m-1.8762-2.2322-2.2104-2.4068-2.5876-2.1379LS估计值/m-1.9245-2.2345-2.1879-2.3851-2.5400-2.1185误差/mm-48.27-2.2622.4721.6747.6419.41中误差
/mm31.51TLS估计值/m-1.9306-2.2399-2.1940-2.3899-2.5433-2.1240误差/mm-54.36-7.6616.4316.8744.2913.89中误差/mm30.88WTLS估计值/m-
1.9245-
2.2345-2.1880-2.3852-2.5400-2.1185误差/mm-48.31-
2.3022.4221.6247.6019.36中误差/mm31.49
图1 实测数据点位分布
本文分别利用LS法、TLS法和WTLS法对二次多项式曲面拟合模型的待求参数进行求解,再根据得到的参数计算出待定点高程异常值,结果见表2。

由表2可知,TLS
法所求的高程异常与真实值之间的差值最大值为-51.0979 mm,最小值为0.7320 mm,平均值为-2.6420 mm,中误差为26.18 mm;LS法结果的差值最大值为-57.0995 mm,最小值为-6.3347 mm,平均值为-12.7545 mm,中误差为26.26 mm;WTLS法结果的差值最大值为-57.0995 mm,最小值为-6.3347 mm,平均值为-12.7545 mm,中误差为26.26 mm。

由表可知,LS法与真值之间的差值较
大,WTLS法的计算结果与LS法基本相当,而TLS法的计算结果更接近真值。

因此,TLS法在二次曲面拟合中具有较高的精度和准确性。

本文分别利用LS法和TLS法对三次多项式曲面拟合模型待求参数进行求解,再根据求解得到的参数计算待定点高程异常,结果见表3。

由表3可知,TLS法所求的高程异常与真实值之间的差值最大值为-54.3621 mm,最小值为-7.6591 mm,平均值为2.1717 mm,中误差为30.88 mm;LS法结果的差值最大值为-48.2698 mm,最小值为-2.2638 mm,平均值为6.3642 mm,中误差为31.51 mm;WTLS法结果的差值最大值为-48.3059 mm,最小值为-2.3005 mm,平均值为6.3288 mm,中误差为31.49 mm。

虽然最大值和最小值上,LS法和WTLS法的结果要高于TLS法,但是在平均值方面TLS法结果精度却远远高于其他两种方法。

再对比LS法和WTLS法的结果,会发现WTLS法结果的平均值要高于LS法的。

因此,从总体上来说TLS法在三次曲面拟合中具有较高的精度和准确性,WTLS法其次,LS法最后。

综上所述,利用二次多项式曲面拟合模型进行参数估计中,TLS法的解算精度是最高的,而LS法和WLS法几乎相同;而换成三次多项式曲面拟合模型以后,还是TLS法的精度最高,但是WTLS法虽然低于TLS法,但是略高于LS法。

本文实验结果表明已知点的平面坐标测量噪声的确对参数估计精度产生了影响,因此在进行高程转换计算时必须考虑这一影响,否则会导致计算结果精度变差。

总的来说,在处理相关测量误差方面,TLS法要优于WTLS法和LS法,并且三次多项式曲面拟合模型的拟合效果要明显高于二次多项式曲面拟合模型。

3 结束语
本文针对GPS测量噪声对高程转换参数精度的影响,基于二次多项式曲面拟合模型和三次多项式曲面拟合模型,分别运用LS法、TLS法和WLS法对实测数据进行高程转换参数的估计。

由数值计算结果可知,GPS平面坐标测量噪声会对转换参数拟合精度产生一定的影响,必须对这部分噪声进行处理。

相对于LS法和WTLS法,在削弱这种噪声影响方面TLS法具有很好的效果,而WTLS法又要优于LS法。

参考文献
【相关文献】
[1] 李征航,黄劲松.GPS测量与数据处理[M].湖北武汉:武汉大学出版社,2016.
[2] 徐绍铨,李振洪,吴云孙.GPS高程拟合系统的研究[J].武汉测绘科技大学学报,1999,24(4):336-340.
[3] 崔立鲁,何秀凤,罗志才,等.Galileo系统定位性能仿真模拟分析[J]测绘信息与工程,2007,32(4)10-12.
[4] 崔立鲁,罗志才,钟波,等.观测噪声对频域输入输出法数据融合的影响[J].大地测量与地球动力学,2009,29(1):79-82.
[5] 崔立鲁.观测噪声对重力场数据融合方法的影响[J].成都大学学报(自然科学版),2016,35(2):183-186.
[6] GOLUB H G, VAN LOAN F C. An Analysis of the Total Least-Squares Problem[J].SIAM JNumer Anal,1980(17):883-893.
[7] 丁海勇,孙景领.GPS高程转换的总体最小二乘方法研究[J].大地测量与地球动力
学,2013,33(3):52-55.
[8] 鲁铁定.总体最小二乘平差理论及其在测绘数据处理中的应用[D].湖北武汉:武汉大学,2010.
[9] MARKOVSKY I. The Element-wise Weighted Total Least-Squares Problem[J]. Comput Stat Data Anal,2006(50):181-209.
[10] SCHAFFRIN B, WIESER A. On Weighted Tatal Least-Squares Adjustment for Linear Regression[J].Journal of Geodesy,2008，82(7):415-421.
[11] 崔立鲁,林景泓,安家春,等.两种GPS精密星历拟合算法的性能比较[J].北京测
绘,2018,32(3):266-268.
[12] 崔立鲁,熊旭,杜石,等.GPS卫星数据仿真[J].北京测绘,2018,32(6):670-673.
[13] 崔立鲁,高先明,魏朋志.基于滑动式算法的精密星历拟合方法比较研究[J].成都大学学报(自然科学版),37(1):52-55,83.
[14] 崔立鲁,陈冬梅,汪晓龙,等.BDS三频伪距/载波线性组合周跳探测与修复[J].测绘工
程,2018,27(12):1-3,7.。