高程面转换

高程面转换（共6篇）

高程面转换 篇1

摘要：利用粒子群算法(PSO)对径向基核函数的参数进行优选,结合多面函数对某GPS控制网的高程异常资料进行拟合分析,并与二次曲面拟合的结果进行对比。结果表明优化后的径向基核函数用于多面函数可获取更为精确的高程异常。

关键词：GPS,高程异常,多面函数,RBF,二次曲面,PSO

0 引言

凭借全天候、高精度等优点,全球定位系统(GPS)技术已被广泛地应用于大型控制网的建立中。目前GPS所建立的控制网平面精度可达到亚毫米级甚至更高,然而GPS的高程信息在实际工程中至今未能充分利用。究其原因是GPS高程和工程所需的正常高不属于同一个高程系统,GPS高程是基于WGS-84参考椭球的大地高,而正常高是基于似大地水准面的高程。并且正常高主要是通过水准仪所建立的水平视线进行高差观测获取,这也就意味着正常高与重力铅垂线有关系,而GPS三维坐标主要与测站和卫星之间的几何关系相关[2]。因此,近些年来国内外学者提出了各种将GPS大地高转换至正常高的方法,如地球重力场模型[3~6],数值拟合[7,8]和机器算法[9~11],并最终试图将GPS测高代替传统的水准测量。由于重力场模型的推算需要大量的重力数据,而且计算过程比较复杂,故重力法在通常情况下难于普及。另外,机器算法如神经网络需要设置大量参数(如隐含层数,结点个数以及权重等),而且最常用的BP学习算法还存在局部最优等缺点[12]。近些年来,多面函数以其模型的多样化和灵活性在地磁场建模[13]、地壳形变监测[14]等领域被广泛应用,但用于GPS高程转换则并不理想,原因主要受限于核函数及其参数的选取。不少文献[7,8]已证实了在平坦地区利用二次曲面和样条函数进行GPS水准拟合可达到厘米级精度,而多面函数用于GPS水准拟合的效果又将如何?因此,本文提出利用粒子群算法(PSO)对径向基核函数进行优选,然后探讨经优化的多面函数模型在GPS水准拟合中的应用,并将最终结果与二次曲面拟合进行对比分析。

1 大地高和正常高系统

如图1所示,GPS高程是地面某一点到参考椭球(如WGS-84椭球)的铅垂距离,即大地高HGPS,而实际工程中常采用的高程数据是地面某点到似大地水准面得铅垂距离,即正常高Hnor,两者相差了一个高程异常ε[1]:

根据大地高和正常高的关系,我们可以将GPS高程转换大地高,其原理主要是利用测区中一定数量且分布均匀的公共点(同时包含GPS大地高和水准正常高的水准点),进而按(1)式计算公共点的高程异常,并利用一定的数学法则(条件允许还可结合重力资料)构造某种曲面来逼近似大地水准面,然后推算出测区中未进行水准联测的GPS点(未测点)的高程异常,最终反算未测点的正常高。

2 多面函数拟合模型

多面函数[15]是美国的Hardy于1971年首次提出的,其基本思想是任何数学表面和任何不规则的圆滑表面,总可以用一系列有规则的数学表面总和以任意的精度逼近。

设某点(x,y)处的高程异常ξ的表达式为:

式中ki为待定系数;xi,yi为选定的中心点,n为中心点的个数;Q(x,y,xi,yi)为x,y的核函数,核函数可以根据实际情况而定,本文主要采用径向基核函数

其中λ为核参数。

当公共点的个数大于中心点的个数时,可采用最小二乘法求解待定系数,即:

当多面函数的待定系数求得后,可按式(2)推算待定点的高程异常。

3 PSO优化核参数

PSO[16]是由Eberhart和Kennedy提出的一种模拟鸟群捕食行为的进化算法,该方法不需进行大量参数的调整,而主要通过粒子追随自己搜寻最优解和整个粒子群的最优解来进行优化。本文将利用PSO对径向基核函数的参数进行优选,其主要步骤如下。

(1)初始化粒子群。对粒子进行随机定位(编码)方案,并将每个粒子的二进制表达转化为核参数及所选择的特征子集。

(2)进行多面函数拟合。根据粒子适应度(均方根误差RMSE),确定粒子个体最优pbest和全局最优gbest[17]。

(3)更新粒子的速度

其中:vk是粒子当前速度,rand()为均匀分布在(0,1)区间的随机数,Ti为最大迭代次数,c1,c2分别为加速因子,wmax为初始惯性权重,wmin为最后惯性权重,t为当前迭代次数,w为惯性权重[18],以避免最大速度设置不当引起的算法搜索最优解的能力。

(4)更新粒子的位置,即:

其中:xk是粒子当前位置。

(5)反复迭代步骤(2)~(4),直至达到最大迭代次数。

4 实例分析

4.1 试验数据

本文采用文献[19]中的GPS控制网数据进行建模分析。测区范围约300km2,且地势比较平缓。该控制网共有44个控制点,其空间分布如图2所示,其中控制点间的高差最大为0.5062m。

4.2 计算分析

在高程异常拟合之间先进行核参数的优化。在实际计算过程中,核参数的范围设为[1-108],粒子的规模设为10,加速因子c1和c2均设为2,wmax设为0.8,wmin设为0.2,最大迭代次数设为200,并作为算法终止条件。在优选步骤第(2)步中,笔者选取12个共同点(图2中“ο”注记)进行多面函数的系数求解,然后外加8个点(图2中“Δ”注记)进行适应度的计算,并优选核参数,其历史全局最优适应度轨迹见图3。考虑PSO算法的结果具有一定随机性,因此笔者进行多次运算,最终确定λ为23114000。

确定核函数后,再次运用多面函数对所有的高程异常进行拟合,其结果与二次曲面拟合模型进行了对比,结果见图4。另外,笔者采用内符合精度指标μ1和外符合精度指标μ2[8]对高程转换的精度进行了评定,其结果见表1。

5 结论

本文初步探讨了基于径向基核的多面函数在GPS高程转换中的应用,并检验了粒子群算法(PSO)用于核参数优选的效果。从文中结果可以看出,在核参数合理选取的情况下,基于径向基核的多面函数用于平坦地区GPS高程的转换精度可控制在1cm以内,其拟合精度要比二次曲面拟合模型更为理想。另外,笔者发现粒子群算法(PSO)对于初始参数的设置比较敏感,因此如何更合理地设置这些参数还有待进一步研究。

高程面转换 篇2

模糊神经网络在GPS高程转换巾的应用

介绍了T-S模糊神经网络的.基本原理以及如何确定GPS高程转换的模糊神经网络模型,并采用该模型对实测数据进行了计算分析.结果表明,模糊神经网络能够对小区域GPS高程做出比较准确定的拟合,从而能够为GPS高程转换提供一种较好的方法,能够满足实际工程需要.

作 者：方毅 雒卫民 FANG Yi LU Weimin  作者单位：上海勘测设计研究院,上海,34 刊 名：地理空间信息 英文刊名：GEOSPATIAL INFORMATION 年，卷(期)： 7(4) 分类号：P228.42 关键词：模糊神经网络   GPS高程   高程转换  

GPS高程转换的新方法研究 篇3

关键词：EIGEN-GL04C地球重力场模型,局部地形改正,移去-恢复,GPS高程转换

0 引言

为了克服水准测量在山区或丘陵地区实施困难的缺点, GPS大地高转化为正常高以代替水准测量将是一个有效的途径。GPS高程转换的关键在于精确求解高程异常, 现有的方法主要是GPS/水准几何拟合法和重力法, 几何拟合法包括很多种, 如曲线拟合法、曲面拟合法、多面函数拟合法、BP神经网络算法、遗传算法、蚁群算法等, 其本质都是用一定的数学曲面来拟合似大地水准面, 这种方法在地势平坦地区可达到四等水准测量的精度要求, 但在地形起伏较大地区精度难以保证。重力法是指用局部或全球重力场模型直接求解高程异常或大地水准面差距, 由于重力数据的缺乏以及重力场模型精度不高使得该方法一直得不到实际工程应用。

针对上述两种方法的不足, 本文采用“移去-恢复”的思路, 提出了一种综合利用EIGEN-GL04C地球重力场模型、局部地形改正以及几何拟合法进行GPS高程转换的方法。由于360阶次的地球重力场模型的分辨率大致为110km, 只能反映大地水准面的中、长波变化, 短波部分主要是地形的影响, 通过移去EIGEN-GL04C模型和局部地形改正所计算的高程异常, 剩余的模型误差 (为系统误差) 从理论上说主要是低阶项, 然后采用一定的几何拟合法对似大地水准面进行拟合, 其效果明显优于纯几何法的高程拟合, 最后再恢复所移去的高程异常。试验证实, 该方法在地形起伏较大地区比以往方法更为有效。

1 EIGEN-GL04C地球重力场模型的原理及精度分析

2003年以前, 高阶重力场模型EGM96被认为是全球最先进的重力场模型, 但由于没有采用卫星跟踪卫星技术以及卫星重力梯度技术, 并且缺乏覆盖全球的高精度地面重力测量数据, 使得该模型求解高程异常或大地水准面高的精度停留在分米级 (0.5～1m) 水平。随着高低卫星跟踪卫星CHAMP和低低卫星跟踪卫星GARCE的发射成功, 重力场模型的精度有了较大的提高。EIGEN-GL04C是GFZ (德国地学中心) 于2006年3月31号新推出的360阶全球重力场模型, 该模型不仅使用了GRACE卫星重力探测数据, 还使用了地球动力学卫星LAGEOS重力探测数据以及地面重力探测数据, 从而使该模型的精度提高了近一个数量级, 成为目前公布于民用的精度最高的模型之一。

EIGEN-GL04C模型精度较EGM96模型有了较大的提高, 但其解算似大地水准面的原理仍与EGM96模型类似, 都是将一个逼近地球质体外部引力位在无穷远处收敛到零值的调和函数展开成在理论上收敛的整阶次球谐函数的360阶次的级数。由完全展开到360阶次的球谐函数所求得的扰动位T可表示为:

$\begin{array}{l}{\rm T}(r,\theta ,\lambda )=\frac{G{\rm M}}{r}{\displaystyle \sum _{n=2}^{360}{\displaystyle \sum _{m=0}^n\left(\frac{a}{r}\right)}}{}^n\cdot (1)\\ [\overline{C}{}_{nm}^{*}\text{c}\text{o}\text{s}m\lambda +\overline{S}{}_{nm}\text{s}\text{i}\text{n}m\lambda ]\overline{{\rm P}}{}_{nm}(\cos \theta )\end{array}$

再根据Bruns公式可求得地面上任意一点的高程异常:

$\begin{array}{l}\zeta (r,\theta ,\lambda )=\frac{G{\rm M}}{r\gamma {}_0}{\displaystyle \sum _{n=2}^{360}{\displaystyle \sum _{m=0}^n\left(\frac{a}{r}\right)}}{}^n\cdot (2)\\ [\overline{C}{}_{nm}^{*}\text{c}\text{o}\text{s}m\lambda +\overline{S}{}_{nm}\text{s}\text{i}\text{n}m\lambda ]\overline{{\rm P}}{}_{nm}(\cos \theta )\end{array}$

式 (2) 中 (r, θ, λ) 是以地球质心为坐标原点、Z轴与地球自转轴重合的地心坐标系中的球坐标, $\overline{C}{}_{nm}^{*}‚\overline{S}{}_{nm}$为完全规格化的扰动位系数, 等于地球引力位与正常引力位球谐函数展开式中相应系数之差, γ0为该点的正常重力值, $\overline{{\rm P}}$nm (cosθ) 为完全规格化的缔合Legendre函数, GM= 3.98600415×1014m3s-2为地心引力常数, a=6.37813646×106m为参考椭球长半径。

以太中银铁路某隧道GPS控制网 (15个点) 为例, 分别运用EGM96与EIGEN-GL04C模型所解算的高程异常与实测高程异常比较如图1所示。虚线表示EGM96模型所解算的高程异常与实测值之差, 实线表示EIGEN-GL04C模型所解算的高程异常与实测值之差。通过比较可得, 利用EGM96模型所解算的高程异常与实测值之差要大于0.6m, 而EIGEN-GL04C模型所解算的高程异常与实测值之差则在0.1m以内, 证明该模型的精度要明显优于EGM96模型。

2 地形改正公式比较

地形改正指的是地面点附近地形起伏对该点大地水准面 (似大地水准面) 的影响, 可以看作是重力场的噪声。目前地形改正公式主要有三种, 一是根据地形起伏所引起的扰动位和Bruns公式来推导的, 直接采用牛顿引力公式进行积分计算, 如文献[1]公式如下 (推导略) :

$\begin{array}{l}{\rm T}{}_c=G\rho {\displaystyle {\iint }_{\sigma }{\displaystyle {\int }_{{\rm H}{}_r}^{?{\rm H}}\frac{1}{r}}}dzd\sigma (3)\\ \delta \zeta =\frac{{\rm T}{}_c}{\gamma }=\frac{G\rho }{\gamma }{\displaystyle {\iint }_{\sigma }\frac{h-h{}_r}{r{}_0}}dxdy-\\ \frac{G\rho }{6\gamma }{\displaystyle {\iint }_{\sigma }\frac{(h-h{}_r){}^3}{r{}_0^3}}dxdy(4)\end{array}$

式 (4) 中, Tc为扰动位, G为万有引力常数, r0为积分主体到地面点的平面距离;ρ为质量密度;h为流动点高程, hr为参考面高程;γ为正常重力值。上述公式在中央区域会引起积分奇异, 很多学者在此基础上进行改进, 推导出了能避免中央区域积分奇异的公式, 如文献[2]的分段积分方法和文献[3]的连续积分方法等, 这些方法所求解的结果与公式 (4) 基本相同。

二是将地形改正分为直接影响和间接影响。移去大地水准面外部地形质量的影响称为地形的直接影响, 而恢复地形质量带来的影响称为间接影响。直接影响一般应用Stokes公式来求解, 但需要地面点的实测重力数据, 对实际工程一般难以获得。间接影响采用Helmert凝聚改正法, 可得公式为:

$\begin{array}{l}\delta \zeta =\frac{\pi G\rho }{\gamma }h{}_p^2-\frac{G\rho }{6\gamma }{\displaystyle {\iint }_{\sigma }\frac{h{}^3-h{}_p^3}{r{}_3}}dxdy+\\ \frac{3G\rho }{40\gamma }{\displaystyle {\iint }_{\sigma }\frac{h{}^5-h{}_p^5}{r{}^5}}dxdy+\cdots (5)\end{array}$

式 (5) 中, hp为待求点的高程。实际工程应用时, 由于缺乏重力数据, 可暂时将直接影响忽略掉, 直接采用式 (5) 作为地形改正公式。

三是将地形对大地水准面 (似大地水准面) 的总影响 (直接影响和间接影响) 统一用一个公式来表示, 如下[5]:

$\delta \zeta =-\frac{2\pi G\rho }{\gamma }h{}_p^2(6)$

式 (6) 表明, 地形对大地水准面 (似大地水准面) 的总影响与计算点高程的平方成正比, 与式 (4) 和式 (5) 的方法相比, 该方法不需要积分, 更加简洁, 更加便于计算。

对照式 (4) 、式 (5) 和式 (6) , 式 (4) 考虑到参考面高程的取值, 但几乎没有考虑到待求点高程对其本身高程异常的影响, 而只跟其附近地形起伏有关, 所以说式 (4) 没有完全把地形的影响考虑进去;式 (6) 不需对参考面高程进行取值, 但是只考虑到待求点高程对其本身高程异常的影响, 而没有考虑到附近地形起伏对它的影响;式 (5) 不需对参考面高程进行取值, 并且指明了该点地形改正与其高程的平方项成正比, 同时也考虑了地形起伏对地形改正的影响。所以从理论上来说, 在进行地形改正时式 (5) 要优于式 (4) 和式 (6) 。

式 (4) 和式 (5) 的积分运算不可能严格完成, 实际计算时只能利用局部地区DEM数据, 根据测区情况将测区进行适当的网格划分, 求出各格网点地形起伏对地面点的影响值, 然后叠加, 进而算得达到一定精度需要的局部地形改正值。网格划分越密, DEM分辨率越高, 局部地形改正越准确, 若测区内没有DEM数据, 也可利用GPS点位上的大地高来内插所有格网的概略大地高, 这样必然存在内插误差, 但也可以适当地进行地形改正, 理论上来说总比不进行地形改正要好。

3 “移去-恢复”原理

根据物理大地测量学的理论, 高程异常可表示为[6]:

$\zeta =\zeta {}_{G{\rm M}}+\zeta {}_{RES}+\zeta {}_{{\rm T}}(7)$

其中ζGM是由重力场模型所计算的中长波项, ζRES、ζT分别为地形的直接影响和间接影响, 统称为地形改正δζ, 也叫做短波项。通过EIGEN-GL04C模型采用式 (2) 可计算出ζGM, 利用式 (4) 或式 (5) 可计算出局部地形改正δζ, 然后将ζGM, δζ移去可得到残差大地水准面ζc=ζGM-δζ, 残差大地水准面是由于模型误差、公式误差或系统误差而引起的, 具有一定的随机性, 可通过数学函数来拟合。实际工程中大多采用二次曲面函数来拟合, 拟合函数如下所示:

$\begin{array}{l}{\zeta }^{\prime }{}_c=\zeta {}_c+a{}_0+a{}_{1}x+a{}_{2}y+a{}_{3}xy+\\ a{}_{4}x{}^2+a{}_{5}y{}^2+\cdots (8)\end{array}$

若已知GPS水准点多于6个, 需要采用最小二乘原理求解。最后将求解的ζ′c加上之前移去的ζGM, δζ即可恢复为该点的高程异常。

4 实例分析

以太中银铁路某隧道洞外GPS控制网为例, 该山区属于低山区地貌, 地形起伏较大, 黄土冲沟发育, 多呈“V”形。隧道全长7.631km, 在隧道进出口以及斜井处共布设15个GPS点, 并按三等水准测量的要求对各点进行了联测, GPS网图如图2所示。

首先不加入EIGEN-GL04C模型及地形改正, 以测区内均匀分布的7个点 (进1、GPS202、GPS201、1X1、2X1、GPS7847、GPS7848) 作为已知高程异常的点, 其它8个点作为检核点, 直接采用二次曲面函数进行拟合, 称为方法一。然后只加入EIGEN-GL04C模型后重新利用二次曲面函数进行拟合, 称为方法二。最后分别利用地形改正公式 (4) 、 (5) 、 (6) , 对该测区进行1km×1km的网格化, 利用DEM数据内插出各格网点的高程, 从而采用“移去-恢复”方法进行GPS高程转换, 分别称为方法三、方法四、方法五。将五种方法所解算的高程异常与实测高程异常之差列于表1, 并分别计算各种方法的内符合精度和外符合精度, 公式如下:

$\mu =\pm \sqrt{\frac{[VV]}{n-1}}(9)$

式 (9) 中V表示已知点或检核点的拟合高程异常与实测值的差值, n表示已知点或检核点个数。通过表1的比较表明:

(1) 方法二较方法一的内外符合精度稍有提高, 但效果不明显, 这也说明EIGEN-GL04C模型解算的只是高程异常的中长波项, 并不能完全反映短波的变化。

(2) 方法三的数据结果表明, 利用式 (4) 所解算的地形改正把拟合数据带坏了, 外符合精度高达29cm, 说明式 (4) 不适合于该算例的地形改正。

(3) 方法四的数据结果表明, 利用式 (5) 所解算的地形改正很好地将高程异常的短波项反映出来, 虽然内符合精度有4mm, 但是外符合精度大大提高, 说明利用式 (5) 的“移去-恢复”方法适合于该算例的地形改正。

(4) 方法五的数据结果表明, 利用式 (6) 所解算的地形改正也把拟合数据带坏了, 说明该公式没有反映出地形起伏对高程异常的影响。

5 结论

(1) 综合利用EIGEN-GL04C地球重力场模型、局部地形改正以及几何法的“移去-拟合-恢复”方法是目前精度较高而实用的方法, 只要拟合区域有DEM数据和少数GPS与水准公共点, 无需实测重力数据, 就可应用于实际工程, 尤其是在地形起伏较大地区, 成为较理想的GPS高程转换途径。

(2) EIGEN-GL04C地球重力场模型能将我国大地水准面的中长波项很好地反映出来, 精度较以往模型有很大的提高, 部分地区可达到厘米级精度, 从而为GPS高程转换提供了一定的基础, 也为局部或区域性大地水准面的精化提供了一个参考模型。

(3) 无论从理论上还是本次实验验证上来说, 式 (5) 都比式 (4) 和式 (6) 更适合于地形改正。若有该地区的实测重力数据, 可同时算出地形对大地水准面的直接影响, 与式 (5) 累加一起作为地形改正公式, 更有利于提高GPS高程转换的精度。

参考文献

[1]徐绍铨, 李振洪, 吴云孙.GPS高程拟合系统的研究[J].武汉测绘科技大学学报, 1999, 24 (4) :338.

[2]熊永良, 路伯祥.高山区GPS网正常高求解方法[J].西南交通大学学报, 1997, 32 (2) :154~155.

[3]张同刚, 岑敏仪, 冯义从等.地形起伏GPS工程控制网高程异常的影响[J].铁道学报, 2005, 27 (2) :80.

[4]罗志才, 陈永奇, 宁津生.地形对确定高精度局部大地水准面的影响[J].武汉大学学报 (信息科学版) , 2003, 28 (3) :340~341.

[5]Wichiencharoen C.The Indirect Effects on the Computation of GeoidUndulations[R].Columbus:The Ohio State University, 1982.

高程面转换 篇4

利用GPS测量的方法代替常规的工程水准测量, 是目前GPS测量研究的一个热点。许多研究例子表明在较为平坦的地区和较小的作业范围, 采取拟合逼近的方法, GPS水准的结果可以达到常规工程水准的精度要求。但是在地势起伏的山区, 或者测量范围比较大, GPS水准测量的精度还不能满足工程水准的精度要求, 这是目前制约GPS测量在高程测量中应用的瓶颈。解决的方法是增加重力测量的数据, 在现有的全球重力场模型的基础上, 精化局部 (似) 大地水准面。但是增加重力测量无疑要加大测量的外业工作量, 而且为了确定局部精确 (似) 大地水准面模型, 所要进行重力测量的区域一般来说要大于工程测量的区域, 所以通过增加重力测量的方法在实际测量工程中的应用可行性不大。不过作为一个地方政府或者国家的基础测量建设, 这种通过联合重力和GPS水准的方法来精密确定似大地水准面的方法, 还越来越受到人们的重视。美国推出的GEOID96, 就是利用这一方法的典范。另一方面, 随着测量资料的丰富, 包括全球重力测量数据, 卫星测高数据等, 全球重力场模型的精度越来越高。EGM96模型重力场模型就是这样一个综合利用现有全球测量数据所计算出来的高精度全球重力场模型。按数据的来源划分, 求解重力异常的方法可以分为两大类, 即重力测量法和几何测量法。重力测量法就是在野外进行重力测量, 再根据斯托克司边值理论或者莫洛金斯基边值理论求解以确定重力异常;几何测量法就是用GPS确定点的大地高, 再进行水准联测确定点的正常高, 两者相减就得点的高程异常。

2 引入重力场模型改正的“移去-恢复”法

高程异常ξ可以分解为ξGM、ξ0G、ξT三部分, 即:

ξ=ξGM+ξ0G+ξT (1)

其中ξGM表示长波部分, 可以通过重力场模型计算得到;ξ0G表示中波部分, 可以通过求解重力异常的边值问题得到;ξT表示短波部分, 可通过求解地形改正得到。作者根据它们与实际的似大地水准面的逼近程度, 把ξGM部分描述为平滑大地水准面, 把ξGM+ξT部分描述为亚平滑似大地水准面, ξGM+ξ0G+ξT描述为详细似大地水准面。

2.1 “移去-恢复”法的原理

地球重力场模型是指重力位的球谐函数级数展开的系数, 简称位系数, 它是利用最新卫星跟踪数据、地面重力异常数据、卫星测高等重力场信息计算得到的。根据给定重力场模型的位系数 (Snm, Cnm) , 可用下式计算各个位置的高程异常:

undefined

式中:ρ, ψ, λ为计算点的地心向径、地心纬度和经度;GM为引力常数与地球质量的乘积;γ为计算点的正常重力值;a为参考椭球的长半轴;undefinednm, undefinednm为完全规格化位系数;undefinednm位完全规格化Lagrandre函数;N为地球重力长模型展开的最高阶数。

若给定一组位系数Snm, Cnm和其参考椭球的基本参数就意味着确定了一个的重力场模型, 也就可以求得相对应模型值。

一般来说, 利用全球重力场模型求解高程异常, 其绝对精度在米级, 因而难以直接用于生产应用。但重力场模型包含较准确的中长波信息, 可用于GPS高程转换中以改善转换的精度。所以, 笔者将GPS点的高程异常分为两部分求解, 即:

ξ=ξGM+ξC (3)

式中:ξGM为重力场模型求得的高程异常:ξC为实际高程异常与由模型求得的高程异常的差值。

2.2 “移去-恢复”算法改进过程

通过若干个己知大地高和正常高的GPS点, 则可以用“移去-恢复”法来求得其它未知点的高程异常, 最终得出未知点的正常高。其实现大体分以下三步:

(1) 移去:设有m个GPS水准联测点, 则可求得此m个点的高程异常ξk=hk-Hk (k=1, 2, …, m) 在这些点上用地球重力场模型, 根据式 (3) 计算出近似高程异常ξundefined, 最后得出剩余重力异常ξundefined=ξk-ξundefined。

(2) 拟合:以m个点的剩余重力异常ξC作为己知数据, 用常规拟合方法或者四参数模型、五参数模型计算出拟合模型的拟合系数, 再内插出未知点的剩余高程异常ξundefined。

ξundefined=x0+x1cosφicosλi+x2cosφ2sinλi+x3sinφi (4)

ξundefined=x0+x1cosφicosλi+x2cosφisinλi+x3sinφi+x4sin2φi (5)

式中xi为模型参数;φi为地心纬度;λi为地心经度。

(3) 恢复:在未知点上, 由公式 (4-2) 计算出近似高程异常, ξundefined, 再加土未知点的剩余重力常ξundefined, 得未知点的最终重力异常值:ξi=ξundefined+ξundefined。从而求未知点上的正常高:Hi=hi+ξi。

3 算例

为了更深层地理解“移去-恢复”法的方法原理及其在GPS高程转换过程中的可行性和有效性, 在这里我们用试验数据来进行算例分析。

用全球重力场模型和拟合的方法、多项式拟合的方法, 分别对试验数据运用“移去-恢复”法进行GPS高程转换。转换前各种模型计算的GPS水准联测点的点位分布图见图1, 各个模型的高程异常的残差比较见图2, 利用纯重力场模型和“移去-恢复法”对比的精度统计见表1。

4 结论

高程面转换 篇5

建立工程控制网坐标系，首先需考虑长度变形不超过±2.5cm/km。对于大部分高原地区，采用国家大地平面坐标系的投影归算面(如IAG-75参考椭球)不能满足这一要求[1,2]。为了减少投影变形，便于工程使用，要求由转换的控制点坐标直接反算的边长值与实测值之间的系统性差异较小，一般是以最重要的某一施工高程面作为坐标系统的最佳抵偿投影面[3]。因此，需要进行坐标投影高程面的转换计算[4,5]。

本文在研究高程面转换计算原理的基础上，提出一种基于TGO软件实现坐标投影高程面的转换计算方法。本文方法可以实现各种参考椭球、任意坐标投影高程面间、单点或批量投影高程面的坐标转换计算，这为快速、准确地进行坐标投影高程面转换计算提供了新工具，在实际工程中具有应用价值。

1 高程面转换计算原理

工程测量控制网在进行投影时，如果要得到某一指定高程面上的坐标，可以有两类方法，直接法和间接法。直接法包括椭球变形法、椭球平移法和椭球膨胀法;间接法则是进行相似变换。目前通常采用直接法中的椭球膨胀法[4]。而在椭球变换时，需要选定一个基准点，这个基准点可以是实际观测的点，也可以是使用多个点归算得到的一个等效的虚拟点。

本文对椭球膨胀法原理加以介绍，示意图如图1 所示。其中，E0为原基础椭球面，P0为地面上的基准点，E0沿P0的法线方向膨胀Δh到所定义的参考面Fh，Δh为图1中P1到P2点的距离，即在E0椭球下沿P0的法线方向穿过参考面Fh的距离。膨胀前后，椭球中心保持不动，椭球扁率α保持不变，椭球长半轴变化，即有:

椭球膨胀方法中计算长半轴膨胀大小的方法有多种，武汉大学测绘学院研制的CosaGPS V5.20提供了三种方法[5]。本文仅介绍一种常用的方法，如下:

其中，Δa为椭球膨胀后长半轴变化量;Δh为椭球面上基准点沿法线方向的膨胀量;dα为椭球扁率变化量。

由于椭球面具有各向异性，所以椭球膨胀后原法线P0→P1在E1椭球下不一定再与椭球面E1垂直。这时有:

其中，Rm、R′m分别为基础椭球和膨胀椭球的平均曲率半径，所对应的大地纬度分别为Bm、B′m;M、N分别为基础椭球的子午圈曲率半径和卯酉圈曲率半径;a、a1分别为基础椭球和膨胀椭球的长半轴;e为椭圆的第一偏心率;B′、B分别为膨胀椭球和基础椭球的大地纬度;dB为椭球膨胀后大地纬度的变化量;W为中间变量。这样就生成新椭球E1，在这个新椭球下，可以重新计算各点新的大地坐标(B,L,H)E1。

2 应用TGO实现高程面转换计算方法

以CGCS2000坐标系统为例，介绍基于美国天宝公司研制的随机软件Trimble Geomatics Office软件[6](简称TGO)实现坐标投影高程面转换计算方法的具体流程，概略流程图如图2所示。

具体实现步骤如下:

(1)椭球系统的建立。启动TGO实用程序“Coordinate System Manager”，进入坐标系统管理器。首先选择“椭球”标签，在空白处点击鼠标右键，选择“添加新椭球”。在弹出的椭球属性对话框中输入CGCS2000坐标系统相应椭球参数，即建立了基础椭球系统(如CGCS2000)。继续“添加新椭球”，采用椭球膨胀法，建立长半轴经过膨胀后的新椭球系统(如CGCS2000+1000，表示长半轴膨胀了1000m，扁率保持不变)。

(2)基准转换的建立。选择“坐标转换”标签，在左侧空白处点击鼠标右键，选择“添加新的基准转换参数→Molodensky”。在弹出的基准转换属性对话框中输入相关参数信息，椭球选择原基础椭球(如CGCS2000)。继续“添加新的基准转换参数→Molodensky”，输入相关参数信息，选择新椭球系统(如CGCS2000+1000)。

(3)坐标系统的建立。选择“坐标系统”标签，在左侧空白处点击鼠标右键，选择“添加新的坐标系统组”。在弹出的“坐标系统组参数”对话框中输入名称，如CGCS2000。选中建立好的坐标系统组，点击右键选中“添加新的坐标系统→横轴墨卡托投影”，弹出“投影带参数”对话框，输入名称，如117(表示中央子午线为117°)。基准名称选择“CGCS2000”，基准方法选择“Molodensky”，“下一步”进入“大地水准面模型”对话框，方法选择“无大地水准面模型”，“下一步”进入“投影”对话框，输入相应的参数。“下一步”进入“移位网格”对话框，选择移位网格名“无”，“完成”即建立了第I投影高程面区，投影面大地高为HⅠ。采用同样的方法即可建立第Ⅱ投影高程面区，取名为117ex(保持中央子午线不变，只变化高程面)，投影面大地高为HⅡ。然后，保存退出坐标系统管理器。

(4)高程面转换计算。打开“Trimble Geomatics Office”软件，新建一个工程项目，如“CGCS2000高程面转换计算”，在“项目属性”对话框中“坐标系统”下改变“坐标系统设置”，选中已配置的坐标系统组CGCS2000，投影带选第I高程面区，如117;在“单位和格式”标签中将坐标小数改为4位小数(便于与CosaGPS处理结果比较)。然后，选择“文件→导入”，选择“自定义”标签，将自定义格式的基于CGCS2000坐标系统下，投影面大地高为HⅠ的坐标文件(*.xyH)导入TGO里。如果遇到格式不对的点，会提示出错信息，修改后再继续。然后，选择“文件→导出”，选择“自定义”标签下的“选项”按钮，改变“导出数据坐标系统”为第Ⅱ投影高程面区，如CGCS2000,117ex，确定后即可获得第II高程面区下的基于CGCS2000坐标系统下的平面坐标。

3 算例验证

为了验证应用TGO软件实现坐标投影高程面转换计算方法的正确性和可靠性，本文算例选取某实测工程控制网某坐标系统(如BJ54和CGCS2000)下，中央子午线为L0I(如117°)，投影面大地高为HI(如0m)的97个测量控制点的平面坐标(x,y)数据资料。若每个测量控制点的高程为Hi(数值可以任意设置，假设Hi均为0m,i=1,2，…，97)，这样即可构成点位坐标文件(*.xyH)，记为第I投影高程面区;假定第Ⅱ投影高程面区是与第I投影高程面区具有相同的坐标系统，中央子午线为L0I，但投影面大地高为HII(如1000m)。为了减少投影变形，现须将这97个控制点的平面坐标(x,y)I从第I投影高程面区通过高程面的转换计算变换至第II投影高程面区下平面坐标(x,y)Ⅱ。

根据本文提出的高程面转换计算方法，我们将这97个测量控制点的平面坐标(x,y)从第I投影高程面区通过高程面的转换计算变换至第II投影高程面区。为了验证本文方法的有效性和可靠性，通过多次实验测试(不同参考椭球、不同投影高程面以及单点或批量转换计算)，并与CosaGPS软件“坐标换带与高程面转换计算”工具处理的计算结果进行了作差比较与分析。限于文章篇幅，仅将BJ54和CGCS2000坐标系统下中央子午线为L0I=117°，投影高程面大地高从HI=0m变换到投影面大地高为HⅡ=1000m的TGO计算结果与CosaGPS高程面转换计算结果的较差进行了数值统计与分析，比较结果如图3和图4所示。

从图3、图4可以反应出，应用TGO实现基于BJ54系、CGCS2000坐标系统的坐标投影高程面转换计算的结果与商用CosaGPS软件处理的结果在平面x方向最大差值分别为0.5mm和0.1mm，在平面y方向最大差值均为0.1mm。考虑到投影方式之间的微小差异或计算误差对高程面的转换计算带来的影响，并忽略对此的影响，在工程测量应用中我们可以认为，采用TGO软件实现坐标投影高程面的转换计算与CosaGPS软件处理的结果是一致的。同时，应用TGO软件对现有的其他坐标系统，如WGS-84、WRS-80、西安80，甚至用户自定义坐标系统所对应的参考椭球、任意投影高程面间转换、单点或批量的高程面转换计算做了多次实验测试与验证，所得结果是一致的。因此，应用TGO软件实现坐标投影高程面的转换计算方法是行之有效的，计算结果是正确可靠的。

4 结束语

本文研究了应用TGO实现坐标投影高程面的转换计算方法。以CGCS2000系统为例，介绍了本文方法的具体计算流程。结合某工程控制网控制点的平面坐标数据资料，对BJ54、CGCS2000以及现有其他系统如WGS-84、WRS-80、西安80，甚至用户自定义坐标系统所对应的参考椭球、任意高程面间转换计算、单点或批量的高程面转换计算做了多次实验测试，计算结果表明本文方法的计算结果与商用CosaGPS处理结果是一致的，从而验证了本文提出的方法是行之有效的，计算结果是正确可靠的。本文方法简单、实用，非常适用于工程测量中高程面的转换计算。同时，如何研究应用TGO软件实现单一的坐标投影换带计算方法以及同时改变中央子午线经度与投影高程面的高程，实现坐标换带与改变投影高程面两项内容的综合计算，将在另文讨论。

摘要：本文提出了一种基于TGO实现坐标投影高程面的转换计算方法。以CGCS2000为例,介绍了本文方法的具体计算流程,对某工程控制网的数据资料进行实验测试,结果表明本文方法的计算结果与CosaGPS处理结果是一致的,从而验证了本方法是行之有效的。

关键词：高程面转换,CGCS2000,坐标投影,TGO

参考文献

[1]赵建三,封良泉,陈宗成.具有抵偿面的任意带高斯投影法的应用[J].工程勘察,2009,(7):59～62.

[2]邱云峰,倪津.不同投影归算面间的坐标换算[J].测绘通报,2001,(9):12～13.

[3]李祖锋.限制边长投影变形最佳抵偿投影面的确定[J].工程勘察,2010,(2):75～78.

[4]田雪冬,郭际明,郭麒麟等.GNSS定位技术在水利水电工程中的应用[M].武汉:长江出版社,2009.

[5]郭际明,罗年学.GPS工程测量网通用平差软件包CosaGPSV5.20使用说明书[R].2010.

高程面转换 篇6

关键词：MATLAB拟合工具箱,GNSS,高程转换,精度分析

0 引言

工程建设所采用的高程一般为正常高, 目前常利用GNSS技术实现快速正常高测量, 其关键是利用一定的拟合方法或区域似大地水准面模型将GNSS测量获得的大地高转换为工程所用的正常高[1~5]。选择有效的拟合方法或者利用高精度的区域似大地水准面模型, 可以快速精确地获取地面点的正常高, 为工程项目节省了宝贵的时间, 提高了工作效率[6,7]。这两种方法的实现途径一般有两种, 利用商用的GNSS高程转换软件或自己编写GNSS高程转换软件[8,9]。若采用商用软件则要受商用软件功能的限制, 且需要购买, 而若自己编写程序, 则需要付出大量的时间, 且需要验证算法是否正确, 工作量比较大。MATLAB是一套高性能数值计算和可视化的软件, 它集数值分析、矩阵运算、信号处理和图形显示于一体, 在系统建模和仿真、科学和工程绘图及应用等方面有着广泛的应用[10~14], 其附带的拟合工具箱Surface Fitting Tool能够提供多种常用的拟合模型, 是科研和工程技术人员简捷的数据分析工具。本文主要探讨利用MATLAB附带的Surface Fitting Tool拟合工具箱实现GNSS高程转换的方法, 并分析其所能达到的精度。

1 GNSS高程转换原理

地面任意GNSS点P的大地高与正常高具有如下关系[1],

其中, H为大地高, 可利用GNSS方法测定, I为正常高, 可利用常规水准方法测定, ξ为该点对应的高程异常。GNSS高程转换的关键在于高程异常的求定, 在局部区域内, 高程异常ξ可用一组线性无关的函数表示[1]:

通过利用公共已知点数据 (GNSS/水准点) , 建立模型函数f (x, y) , 然后求出区域待求点的高程异常, 实现高程转换。

2 MATLAB拟合工具箱简介及使用

本文所使用的MATLAB拟合工具箱是MATLAB7.11.0 (R2010b) 版本附带的“Surface Fitting Tool”, 更高版本的MATLAB所附带的拟合工具箱具有相同的基本功能和操作界面[10~14]。Surface Fitting Tool适用于具有拟合要求的数据分析工作, 用户能够通过操作图形用户界面和调用函数这两种方式以实现数据处理[10~14], 该工具箱的交互界面简洁友好, 方便操作, 其运行界面见图1。Surface Fitting Tool包含了众多拟合函数, 能满足大部分的数据分析工作, 其常用的拟合函数见表1, 同时, 该工具箱还提供开放的公式编辑接口, 用户可以根据自身的要求而自定义一组拟合函数, 进行相应的拟合操作, 整个过程不需要任何编程工作。

Surface Fitting Tool的一般使用流程如下: (1) 按照格式准备已知数据和检验数据; (2) 输入相应的已知数据; (3) 选择拟合模型并设置参数; (4) 点击拟合命令进行拟合, 并得出结果; (5) 查看各种视图及统计信息, 检验拟合效果, 若模型符合要求, 则进行步骤 (6) , 否则返回步骤 (3) ; (6) 输入相应的检验数据; (7) 查看各种视图及统计信息, 检验拟合效果, 若模型符合要求, 则输出结果, 否则返回步骤 (3) 。关于Surface Fitting Tool更详细的使用说明可查阅MATLAB软件相关的帮助文档。

3 算例及精度分析

本文采用某一GNSS控制网数据进行实验, 该控制网位于华东地区, 共有54个GNSS/水准点, 控制网东西跨度约44km, 南北跨度约32km, 面积约1408km2, 整个测区地形起伏不大。为了研究工具箱中拟合函数的效果, 需要对GNSS/水准点进行分组。在测区内均匀地选择其中15个GNSS/水准点作为检验点, 剩余的39个作为备用点并根据不同的实验方案选取为已知点, 点位分布见图2。本文设计了两个实验, 其主要的实验步骤如下:首先将所有控制点的高程异常求出;然后按照格式准备好已知数据和检验数据;根据选点方案输入已知数据, 建立模型后, 利用检验数据作拟合模型的外部检验;最后对检验数据的RMS值进行整理, 得到相关的统计结果。

为了研究已知点数量对拟合效果的影响, 首先进行实验1, 实验方案如下:将上述的15个GNSS/水准检验点作为检验数据, 然后在剩余的39个备用点中选取15个作为已知点, 并建立拟合模型, 得到模型后, 利用检验点进行检验, 从而得到模型外部精度的RMS值, 接着在15个已知点的基础上继续等数目地增加已知点个数 (每次均匀地增加6个) , 重新建立模型并进行检验, 直至所有的备用点都被选取为止。需要注意的是, 由于工具箱中部分拟合函数不能进行外推, 且外推的精度较低, 所以第一轮选择的15个已知点需要尽量地分布在测区四周, 较好地涵盖测区范围。各方案的点位分布如图3所示, 得到RMS数据后绘制得图4。

从图4可以看出, 随着已知点个数的增加, 最邻近拟合模型和Biharmonic (v4) 样条拟合模型的外部精度有了明显的提高, 最邻近拟合模型的RMS值由13.7cm降到最优的7.5cm, Biharmonic (v4) 样条拟合模型的RMS值由6.6cm降到最优的3.5cm;但是, 当已知点个数达到一定数目时, 其精度逐渐趋向平稳, 变化不大。线性拟合模型、三次样条拟合模型、多项式拟合模型和线性/二次局部加权平滑拟合模型的精度变化不大。除了最邻近拟合模型和Biharmonic (v4) 样条拟合模型外, 其余模型均随着已知点个数的增加, 其外部精度的RMS值在3.0cm上下波动。在已知点个数为27时, 最优的RMS值为二次曲面拟合的2.6cm。对于每种选点方案, 通过比较可知, 在多项式拟合模型中, 二次曲面拟合模型达到的精度要优于平面拟合模型和三次曲面拟合模型, 另外, 窗宽为50%的线性/二次局部加权平滑拟合模型优于窗宽为25%时的情形。

(注:25%、50%为窗宽, 即是模型核函数的平滑因子, 数值越大, 拟合效果越平滑)

同样采用上述的15个GNSS/水准点作为检验数据, 然后通过不同的已知点选取方案, 研究工具箱中各种拟合模型对点位分布的敏感程度。已知点的个数为27, 共有5种选点方案, 其对应的点位分布见图5, 由此得到外部精度的RMS数据见图6。

从图6可以看出, 最邻近拟合模型和Biharmonic (v4) 样条拟合模型对点位的分布十分敏感, 随着点位的改变, 其外部精度的RMS值变化较大。最邻近拟合模型的RMS值最小为8.2cm, 最大为9.8cm, 两者相差达1.6cm;Biharmonic (v4) 样条拟合模型的RMS值最小为3.5cm, 最大为5.1cm, 两者相差同样为1.6cm。线性拟合模型和三次样条拟合模型的RMS值的变化幅度较小, 分别为0.7cm和0.8cm。相比之下, 线性/二次局部加权平滑拟合模型的稳定性都有较好的表现, 其RMS值变化幅度始终保持在0.5cm内。对于多项式拟合模型, 包括平面拟合模型、二次曲面拟合模型和三次曲面拟合模型, 其稳定性最优, RMS的变化幅度都在0.4cm内。结合图5和图6, 对5种选点方案比较可知, 方案1的已知点分布最为均匀, 故大部分拟合模型在该方案中均具有较优异的表现。因此, 在GNSS高程转换工作中, 需要均匀地选取已知点, 并且最好使用多项式拟合模型, 从而得到高精度的转换结果。

(注:25%、50%为窗宽, 即是模型核函数的平滑因子, 数值越大, 拟合效果越平滑)

4 结束语