Kriging方法结合最小二乘配置在GPS高程拟合中的应用

第36卷第1期

2011年01月

测绘科学

Sc i ence o f Survey ing and M app i ng

V o l 36N o 1

Jan

作者简介:李军海(1985 ),男,湖北

黄冈人,硕士研究生,从事G PS 数据处理方面的研究。E m a i:l lij unha i_ch @126 co m 收稿日期:2010 04 13

基金项目:国家自然科学基金(40874012);武汉大学地球空间环境与大地测量教育部重点实验室开放研究基

金(07 05)

K riging 方法结合最小二乘配置在GPS 高程拟合中的应用

李军海

,文汉江 ,方爱平 ,刘焕玲

( 中国测绘科学研究院,北京 100830; 辽宁工程技术大学测绘与地理科学学院,辽宁阜新 123000) 摘 要 本文利用K r i g i ng 方法结合最小二乘配置将GPS 高程转换成正常高。研究了将K r i g i ng 方法中的变异函数用于计算最小二乘配置中的协方差的方法,并对一局部GP S 水准网的高程作了拟合计算。通过将最小二乘配置法与平面拟合模型和多面函数拟合模型等进行比较,其外符合精度从最大的 0 0277m 提高到 0 0162m 。 关键词 最小二乘配置;变异函数;协方差函数;高程拟合

中图分类号 P228 文献标识码 A 文章编号 1009 2307(2011)01 0099 03

1 引言

G PS 用于高程测量所测定的是测站相对于W GS 84椭球面的大地高,但是我国采用的是正常高系统,利用G PS 水准拟合方法可将大地高转换为正常高。常用的高程异常拟合方法有多项式平面拟合,多面函数法拟合,样条函数法拟合等[2]。多项式平面拟合法由于拟合函数始终只能是高程异常的趋势面,与高程异常的实际值比较必然会有一定差异[4]。多面函数法采用多曲面逼近原理,从理论上可以得到更切合实际的数学曲面,但是核函数中的光滑因子对G PS 水准拟合精度的影响十分明显,但是到目前为止,对光滑因子既无准确的定义,又无最优计算方法。只能根据经验反复试算得到[3]。近年来有些学者采用神经网络BP 算法转换GPS 高程。算法虽然具有简单、可塑等特点,但它是一种梯度下降搜索方法,还存在收敛速度慢、容易陷入局部最优、有时找不到全局最小值等缺点[9]。

最小二乘配置是基于统计线性最优无偏估计的方法,在推估过程中同时考虑非随机变量和随机变量,使得高程异常值的推估精度更高。在用最小二乘配置对高程异常拟合过程中关键是预先知道各信号间的协方差,严格地讲,协方差阵应该通过大量的观测资料经统计求得,在实际应用时,由于数据的有限性因而很难应用。K ri g ing 方法是法国著名学者M athe ron G 创立并发展起来的,它是一种基于随机过程的统计预测法,可对区域化变量求线性、无偏和最优内插估计值,具有平滑效应及估计方差最小的统计特征。对K rig i ng 方法中的变异函数进行研究发现变异函数实际上是协方差函数的另一种表达形式。因此,将变异函数作为最小二乘配置的协方差函数可有效地解决协方差的求定问题。

2 一些常用高程异常拟合方法

似大地水准面与椭球面之间的差距称为高程异常,用 表示。

大地高与高程异常关系为:

=H -h (1)

其中:H 为大地高,h 为正常高, 为高程异常。若能精确地求出GPS 点上的高程异常 就可以求出GPS 点的正常高。所以,高程异常的确定就成为实现GPS 高程转换的关键问题。

实际应用过程中对G PS 控制点进行几何水准联测,在数学上计算联测点的高程异常值,根据有限的高程联测点把整个测区的似大地水准面拟合为准面、多项式曲面或其他数学曲面,根据拟合的曲面内插区域内待求点的高程异常,从而将待求GPS 点测定的大地高转换为正常高。常用的方法有平面拟合法,多面函数法以及BP 神经网络法等。2 1 平面拟合法

本文主要是通过平面拟合法和多面函数法来进行比较。下面分别对这两种方法进行介绍。

平面拟合的数学模型为[2]:

i =a 0+a 1x i +a 2y i

(2)式中 i 表示高程异常;x i ,y i 为平面坐标。将式(2)转换成误差方程矩阵形式为:

v =BX - (3)

其中 =[ 1 2 n ] v =

[v 1 v 2 v n ]T X =[a 0 a 1 a 2]T B =

1x 1y 11x 2y 2

1x n y n

通过n(n 3)已知高程异常的水准重合点,用最小二乘法求得拟合系数a i ,就可利用式(2)求得未知点的高程异常。

2 2 多面函数拟合

多面函数拟合法的数学模型为[2]:

= k

i =1 i Q i (x,y,x i ,y i )

(4)

i 为待定系数;Q i (x,y,x i ,y i )为核函数,x,y 为待求点坐标,x i ,y i 为已知点坐标。k 为规则数学表面的总和。选择:

Q i (x,y,x i ,y i )=[(x -x i )2+(y -y i )2+ ]1/2

(5)

其中 为光滑因子,用来对核函数进行调整,本文选取参考点与待定点之间距离最大值作为光滑因子 的值。

将(5)式转写成误差方程式为:

v =Q - (6)

待定系数 可以根据已知点上的高程异常值按最小二乘法计算得到:

=(Q T Q )-1Q T (7)

测绘科学 第36卷

计算出核函数和待定系数后,就可以按照公式(4)计算待定点的高程异常。

3 最小二乘配置法原理

3 1 最小二乘配置

最小二乘配置法的函数模型为:

L =X +G Y + (8)

式中L 为观测向量, 为观测噪声, ~N (0 D ),Y 为待定的非随机参数,X 为观测信号,将未测点观测信号用X 表示。

由广义最小二乘原理得到观测方程为:

L =B Z +GY + (9)

其中:B =[I 0],Z =[X X ]T

已知先验信息E [X ]= X =0,var [X ]=D X ,E [X ]= X

=0,var [X ]=D X ,X 与X 的协方差为D XX =D T

X X 。

相应的误差方程为:

V =BZ ^

+GY ^

-L

(10)

其中:

Z ^

=[X ^

X ^

]T

根据最小二乘原理有:

V T P V +Z T P Z Z =m i n

(11)

其中:

P = 20D -1 ,P Z = 20D -1

Z ,D Z =

D X

D X X

D X X D X 于是结合(10)、(11)式构造极值函数: =V T P V +Z ^T P Z Z ^+2K T (B Z ^+GY ^-L -V )(12)将 分别对V,Z ^和Y ^求偏导:

V

=2V T P -2K T =0 V =P -1

K (13) Z ^

=2Z ^T P Z +2K T B =0 Z ^=-P -1Z B T

K (14)

Y ^

=2K T G =0 K T G =0(15)将(13)、(14)、(15)与(10)式联合计算得到估计值为:

Y ^=[G T (BD X B T +D )-1G ]-1G T (BD X B T +D )-1L

(16)

X ^=D X B T (BD X B T +D )-1

(L -G Y ^)(17)

X ^ =D X X B

T

(BD X B T +D )-1(L -G Y ^)(18)得到未测点上的完全信号为:

L ^ =G Y ^+X ^ (19)其中:G =1x 1-x 0y 1-

y 01x 2-x 0y 2-y 0

1x m -x 0y m -y 0

,m 为未知点个数。

(x 0,y 0):坐标平均值;(x i ,y i

):未测点的坐标值由(16)、(17)、(18)式可以看出,求解Y ^、X ^、X ^ 、的关键是要确定D Z 。严格来说,协方差阵各元素应该通过大量观测数据经统计得出。由于实际条件的限制,这不可能实现。因此本文采用变异函数来确定。本文直接将高程异常 作为无误差的观测值来进行处理。3 2 变异函数

K ri g ing 方法中的区域化变量z (x )在点x 和x +h 处的值z (x )与z (x +h )差的方差称为区域化变量z (x )的变异函数:

(h )=12E [(z (x )-z (x +h ))2]-1

2

E [z (x )]-E [z (x +h)]2

(20)

式中h 对的距离。由上式可知,二阶平稳或本

征假设条件下,随机变量z (x )对于任意的h 满足:

E [z (x +h )]=E [z (x )]

于是变异函数可以改写为:

(h )=1

2

E [(z (x )-(z (x +h)]2(21)

3 3 变异函数与协方差函数的关系

协方差函数D (h )是反映区域化变量z(x )与z (x +h )相关程度的量,通常是距离h 的减函数。变异函数r(h)是反映区域化变量z (x )与z (x +h )变异程度的量。

图1 协方差函数与变异函数可得变异函数与协方差函数的关系为:

(h)=D (0)-D (h )

(22)

实际应用中,可以根据已知高程点与距离构造函数

(h )的估计值 ^

(h )。依据不同的间隔

h 1,h 2, ,h m 计算各间隔上采样值的函数[1],然后采用球状模型拟合关于距离的变异函数:

^

(h m )=

12N (h m ) h m =|(x i ,y i )-(x j ,y j )|

(z (x i ,y i )-z (x j ,y j ))2(23)

其中N (h m )是研究区域距离为h m 的已知高程点对数,在实际工程中,由于进行水准拟合的点位分布比较散乱,因此,首先需要计算出任意两个数据点之间的距离h ij ,然后对距离进行分组h m

。本文采用的分组方法为:h m

=m ax h m -m in h m

N (h m )

选择球状函数作为理论模型:

(h)=C 0+C (3h 2a -h 3

2a 3

)(24)

式中C 0,C 和a 为未知数。

这样由公式(22)得到协方差函数的表达式为:

D (h)=C 0+C - (h)=C (1-(3h 2a -h 3

2a 3

))(25)

其中D (0)=C 0+C ,求出 ^

(h )后,对公式(24)进行最小二乘平差,将已知点代入函数模型组成观测方程,由观测方程可以得出该式的统一误差方程矩阵形式:

V =A - 0(h )

(26)其中 =(C 0,C,a)

在最小二乘准则下进行平差求解,我们可以求出参数向量 ,再将求得的参数C 、a 分别代入公式(25)中就可得到协方差函数表达式,在将所得到的协方差函数代入公式(16),(17),(18)就可得到各参数的估值,再将得到的估值代入公式(19)中即可求得未测点的值。

4 实例分析

本文GPS 高程拟合的算例中,数据来源为某区域平面图2 G PS 水准点

点位分布图

高程控制网。共有22个水准高程点,其坐标按国家G PS D 级网要求施测、正常高按四等水准测量施测。区域总面积近180k m 2,高程异常变化最大1 0863m 。从中选取15个控制点作为已知点进行拟合计算,其余7个点作为检核点。点位分布如图2

所示。100