CN102175213A - Gps水准高程传递方法 - Google Patents

Abstract

本发明提供了一种GPS水准高程传递方法，该方法包括：S1：由待测区域的地面重力观测数据，结合区域数字地面模型和全球重力场模型，基于莫洛金斯基理论构建区域重力似大地水准面模型；S2：利用GPS水准超定边值方法获得GPS大地高与水准高程的参考系统差异及其影响的高程异常改正数，并根据所述重力似大地水准面模型、所述高程异常改正数以及预置的GPS水准高程异常系统偏差构建实用似大地水准面模型；S3：根据所述重力观测数据和所述实用似大地水准面模型获得待测区域待定点的正常高。该方法可以实现GPS水准代替三角高程测量或者常规水准达到二等水准高差的传递精度。

Description

GPS水准高程传递方法

技术领域

本发明涉及大地测量应用技术领域，尤其涉及一种GPS水准高程传递方法。

背景技术

斯托克斯理论是物理大地测量的基本理论，它以大地水准面作为边界面，假设大地水准面外部没有质量，并以大地水准面上的重力测量基本微分方程为基础形成第三边值问题，用于求解扰动位和大地水准面高，从而确定地球的物理形状。斯托克斯理论下的大地水准面、正高系统等物理量具有严格的物理意义，但也存在如下基本缺点：1)其正高无法精确解算；2)解算是需要知道大地水准面上的实际重力，因此需要进行重力归算，对于地壳密度和地形质量需要估算；3)以大地水准面为边界面，属于自由边界问题，不易求解。

为克服上述问题，引入莫洛金斯基理论以及似大地水准面概念。莫洛金斯基边值问题以实际的地球表面为边界面，以地面重力异常为基本观测量解算边值问题，求解高程异常，从而确定似大地水准面的起伏，进而确定地球的形状。基于莫洛金斯基理论，利用正常重力的平均值代替实际重力在正高上的平均值作为“正常高”；由地面沿着正常重力线向下量取正常高所形成的曲面称为“似大地水准面”；似大地水准面到平均椭球面的距离称为“高程异常”。

工程建设和测图中需要使用的高程值是指“正常高”，即：测量点(地球外部空间任一点)到大地水准面铅垂线的长度，它是描述大地水准面起伏的唯一标准，具有物理意义。传统方法是利用水准测量一站一站地传递来获得所述“正常高”，这就要求前、后测站通视，并且地形起伏不能太大，自然人和仪器都需要到达现场。在我国西部山区道路桥隧建设中必须使用“正常高”，但由于西部山区高海拔高起伏的特殊地形，传统方法实施起来存在一定的难度。

因此，可以考虑利用GPS技术代替常规水准测量。然而，通过GPS技术观测所获得的“大地高”是相对于WGS-84椭球的高程值，即：测量点到该椭球面的距离，具有严格的数学意义，属于纯几何意义的高程系统。必须将通过GPS技术观测的“大地高”转换为相同点的“正常高”(GPS水准高程传递)，这就要求设法求得相对于我国似大地水准面的“高程异常”，即：区域似大地水准面模型。

目前，求解高程异常的方法主要有重力法(或称物理大地测量法)和几何拟合法等。重力法需要具有一定精度和分布的重力数据和地形数据；几何拟合法只需知道测区少数点的高程异常就能通过某种函数模型拟合出待定点的高程异常。国内相关领域采用的基于几何拟合法的高程异常求解方法包括：简单的GPS高程拟合，例如，线性内插法、平面拟合法、二次曲面拟合、多面函数法及神经网络法等；以及，顾及重力似大地水准面模型的强制拟合。上述方法运用于山区道路建设、大比例尺测图、跨海高程传递等可以达到四等水准高程传递精度，部分区域可达到三等水准高程传递精度，但仍然无法达到能够代替常规水准测量的二等水准传递精度。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是如何实现GPS水准代替三角高程测量或者常规水准达到二等水准高差的传递精度。

(二)技术方案

为解决上述技术问题，本发明的技术方案提供了一种GPS水准高程传递方法，包括：

S1：根据待测区域的地面重力观测数据，结合区域数字地面模型DTM和全球重力场模型EGM，基于莫洛金斯基理论构建区域重力似大地水准面模型；

S2：利用GPS水准超定边值方法获得GPS大地高与水准高程的参考系统差异及其影响的高程异常改正数，并根据所述重力似大地水准面模型、所述高程异常改正数以及预置的GPS水准高程异常系统偏差构建实用似大地水准面模型；

S3：根据所述重力观测数据和所述实用似大地水准面模型获得待测区域待定点的正常高。

优选地，步骤S1进一步包括：

S101：通过地面重力观测数据获得待测区域离散点的扰动重力；

S102：通过顾及地形影响的移去-恢复方法获得所述离散点的平均扰动重力；

S103：根据所述平均扰动重力和所述全球重力场模型EGM获得零阶剩余高程异常；

S104：根据所述零阶剩余高程异常和所述全球重力场模型获得莫洛金斯基零阶高程异常；

S105：对线性莫洛金斯基I/II项进行系统分解，并结合区域地形进行量化分析，获得高程异常地形改正；

S106：根据所述莫洛金斯基零阶高程异常和所述地面高程异常地形改正获得重力地面高程异常。

优选地，所述重力似大地水准面模型为：重力地面高程异常＝莫洛金斯基零阶高程异常+高程异常地形改正。

优选地，步骤S2中，所述GPS水准超定边值方法进一步包括：

S201：通过移去一恢复法，将“移去”并格网化的重力场元作为统计观测量，以Stokes公式或Hotine公式作为约束条件，建立GPS水准高程异常与格网化后的重力场元之间的条件方程；

S202：基于所述条件方程，通过条件平差法获得所述格网化的重力场元观测量的改正数的估值，并进一步获得格网重力场元平差值；

S203：根据重力场积分公式，由所述格网重力场元平差值计算获得地面高程异常改正数。

优选地，所述实用似大地水准面模型为：实用似大地水准面高＝重力地面高程异常+GPS水准高程异常系统偏差+地面高程异常改正数。

优选地，所述待定点的正常高＝大地高-实用似大地水准面高；其中，所述大地高通过所述重力观测数据获得。

(三)有益效果

与现有技术相比，本发明的技术方案具有以下优势：1)适用范围更广，特别是在我国西部高海拔、高起伏的山区；2)实现精度更高，传统的方法尽可达到三、四等水准高差精度，本方法可达到二等水准高差传递精度。

附图说明

图1是本发明的GPS水准高程传递方法的流程图；

图2是本发明的GPS水准高程传递方法中使用的GPS水准超定边值方法的流程图；

图3是本发明的GPS水准高程传递方法的一个实施例中待测区域的GPS水准和重力点分布示意图。

具体实施方式

本发明提出的GPS水准高程传递方法，结合附图和实施例说明如下。

如图1所示，本发明的GPS水准高程传递方法包括以下步骤S1-3：

S1：根据待测区域的地面重力观测数据，结合区域数字地面模型DTM和全球重力场模型EGM，基于莫洛金斯基理论构建区域重力似大地水准面模型；

本步骤的目的是实现基于莫洛金斯基Molodensky理论的重力似大地水准面模型解算。模型解算的数据源包括：待测区域离散点的重力观测数据(精度优于60毫伽)、全球重力场模型(EGM)、区域数字地面模型(DTM)。

具体地，如图2所示，本步骤的重力似大地水准面模型解算进一步包括：

S101：通过地面重力观测数据获得所述离散点的分别的扰动重力；

对于每一离散点，利用地面重力观测数据中的大地坐标和正常重力公式求得该点的正常重力值；则地面重力观测数据中的实际重力值和计算获得的正常重力值之差即为该点的扰动重力。

S102：通过顾及地形影响的移去-恢复方法获得所述离散点的平均扰动重力；

通过移去-恢复方法对步骤S101计算得到的扰动重力进行顾及地形影响的处理：先移去地形影响，后格网化处理，最后再恢复地形影响，从而得到格网化的平均扰动重力。顾及地形的移去-恢复的目的是降低离散点扰动重力格网化误差。移去-恢复方法属于现有技术，在此不做赘述。

S103：根据所述平均扰动重力和所述全球重力场模型EGM获得零阶剩余高程异常；

将步骤S102计算获得的平均扰动重力减去全球重力场模型EGM可得到剩余扰动重力；利用Hotine公式对该剩余扰动重力直接积分即可获得零阶剩余高程异常。

S104：根据所述零阶剩余高程异常和所述全球重力场模型获得莫洛金斯基零阶高程异常；

将全球重力场模型直接解算得到的模型高程异常加上通过步骤S103获得的零阶剩余高程异常即可得到莫洛金斯基零阶高程异常。即：莫洛金斯基零阶高程异常＝零阶剩余高程异常+模型高程异常。其算法包括数值积分、平面二维FFT算法、球面二维FFT算法以及球面一维FFT算法等。

S105：对线性Molodensky I/II项进行系统分解，并结合区域地形进行量化分析，获得莫洛金斯基零阶高程异常的地形改正；

地形影响及其改正算法的优劣是区域高精度似大地水准面模型精化的瓶颈。本发明的技术方案基于莫洛金斯基Molodenksy理论，对线性Molodensky I/II项进行了系统分解，并结合区域地形进行量化分析，引入地面起伏情况，对其一阶项的三个分项和其二阶项的四个分项做分别解算，获得莫洛金斯基零阶高程异常的地形改正。

具体地，似大地水准面高程异常的Molodensky级数解顾及二阶项可表示为：

ζ＝ζ₀+ζ₁+ζ₂    (1)

其中，高程异常的地形总影响ζ_rf可用下式计算：

ζ_rf＝ζ₁+ζ₂，(2)

ζ₁＝ζ₁₁+ζ₁₂₀+ζ₁₂₁，ζ₂＝ζ₂₁+ζ₂₂₀+ζ_23′+ζ₂₄    (3)

其中，

${\mathrm{\ζ}}_{11}=\frac{R}{4\mathrm{\π\γ}}\underset{\mathrm{\σ}}{\∫\∫}{G}_{11}S\left(\mathrm{\ψ}\right)\mathrm{d\σ},{\mathrm{\ζ}}_{120}=-\frac{\mathrm{\Δg}}{\mathrm{\γ}}h,{\mathrm{\ζ}}_{121}=-\frac{3h}{2R}{\mathrm{\ζ}}_0$

${\mathrm{\ζ}}_{21}=\frac{R}{4\mathrm{\π\γ}}\underset{\mathrm{\σ}}{\∫\∫}{G}_{21}S\left(\mathrm{\ψ}\right)\mathrm{d\σ},{\mathrm{\ζ}}_{220}=-\frac{G_1}{\mathrm{\γ}}h,{\mathrm{\ζ}}_{{23}^{\′}}=-\frac{R^2}{4\mathrm{\π}}\underset{\mathrm{\σ}}{\∫\∫}(\mathrm{\Δg}+\frac{3\mathrm{\γ}{\mathrm{\ζ}}_0}{2R})\frac{{(h-h_P)}^2}{l_0^3}\mathrm{d\σ}$

${\mathrm{\ζ}}_{24}=\frac{R}{4\mathrm{\π\γ}}\underset{\mathrm{\σ}}{\∫\∫}{G}_{24}S\left(\mathrm{\ψ}\right)\mathrm{d\σ},G_{11}=-\frac{{\mathrm{hR}}^2}{2\mathrm{\π}}\underset{\mathrm{\σ}}{\∫\∫}\frac{\mathrm{\Δg}-{\mathrm{\Δg}}_P}{l_0^3}\mathrm{d\σ}$

$G_{21}=-\frac{{\mathrm{hR}}^2}{2\mathrm{\π}}\underset{\mathrm{\σ}}{\∫\∫}\frac{G_1-G_{1P}}{l_0^3}\mathrm{d\σ},G_{24}=(\mathrm{\Δg}+\frac{{3\mathrm{\γ\ζ}}_0}{2R}){\mathrm{tg}}^2\mathrm{\β}$

(4)公式组(4)中，Δg为重力异常，即地面重力异常，将正常重力按照正常高由平均椭球体归算到似地球表面，在数值上与空间重力异常相等；γ为正常重力平均值，由正常重力和正常高求得即γ＝γ₀-0.1543h；h为正常高；R为地球平均半径；S(ψ)为斯托克斯函数，其中ψ为积分流动点和计算点相对于球心的夹角。

把我国地形分为平原、丘陵、一般山区、高山区以及高原平原和大倾角高山区，重力场特征按上述相应公式分析估计得出的我国大陆高程异常地形影响的一般规律如表1所示。

表1

在地形起伏大的山区，由于重力异常与地形的正相关关系得到加强，地形影响中的低频成分表现为正值；但高频成分，如二阶项，表现为正负相间。在山区，地形影响中的二阶项表现得比地形还粗糙的甚高频变化，这意味着它们的大小、起伏将随数字高程模型的分辨率不同发生较大变化，地形影响将变得非常复杂，此时，应根据地形特征选择合适的区域数字地面模型DTM的分辨率，才有可能满足似大地水准面的精度要求。例如，当地形倾角大约30度时，DTM分辨率不得高于10秒。

S106：根据所述莫洛金斯基零阶高程异常和所述地面高程异常的地形改正获得重力地面高程异常。

将步骤S104获得的莫洛金斯基零阶高程异常加上步骤S105获得的地面高程异常地形改正即为重力地面高程异常(重力似大地水准面高)。即：本发明技术方案构建的重力似大地水准面模型为：重力地面高程异常(重力似大地水准面高)＝莫洛金斯基零阶高程异常+高程异常地形改正。

S2：利用GPS水准超定边值方法获得表征GPS与水准的参考系统差异的偏差改正的高程异常改正数，并根据所述重力似大地水准面模型、所述高程异常改正数以及GPS水准高程异常系统偏差构建实用似大地水准面模型。

本步骤利用GPS水准超定边值方法实现GPS与水准的参考系统及其影响的偏差改正，从而有效降低不同基准差别以及系统偏差对高程传递的影响，特别是难以量化的参数影响；以及，有效抑制长波效应和远区影响。为有效地消除GPS水准实测高程异常与重力地面高程异常的不一致，需要计算两种高程异常之间的差异，并从差异中扣除系统偏差信息，获得的结果称为GPS水准残差高程异常(可按点值计算)。即：GPS水准残差高程异常＝GPS水准实测高程异常-重力地面高程异常-GPS水准高程异常系统偏差。此后，通过GPS水准与重力的联合平差(GPS水准重力融合)，得到实用的似大地水准面高模型。

具体地，本发明的技术方案提供的GPS水准超定边值方法(GPS/Leveling Over-determined BVP)表述如下：应用移去一恢复法，将“移去”并格网化的重力场元(重力异常、扰动重力等)作为统计观测量，并将Stokes公式或Hotine公式作为约束条件，建立GPS水准高程异常与格网化后的重力场元之间的条件方程；按条件平差法估计格网重力场元观测量的改正数，在此基础上按重力场积分公式由格网重力场元平差值计算获得地面高程异常改正数。

在本发明的一个实施例中，以Stokes积分作为超定边值问题的约束条件，并以格网重力异常作为统计观测量，则地面高程异常改正数的计算方法如下：在计算区域内设定一起算点，设区域性高程基准差别、水准高程从零点传递到起算点的系统偏差、以及其他因素引起大地高的系统偏差总和为a₀，在计算区域＜5°×5°时可认为基准差别引起大地高变化与纬度差ΔB和经度差ΔL成正比，认为从起算点到计算区域内任意GPS水准点j的水准高程传递系统差与高差大小Δh存在某种函数关系。则GPS水准点j的实测残差高程异常可表示为：

${\mathrm{\Δ\ζ}}_j={\mathrm{\Δ\ζ}}_0^j+a_0+a_1\mathrm{\ΔB}+a_2\mathrm{\ΔL}+a_3{D}^2\mathrm{sign}\left(\mathrm{\Δh}\right)\ln (1+{\mathrm{\Δh}}^2)---\left(5\right)$

公式(5)中，

为无基准差别和系统偏差时GPS水准点j的残差高程异常；a_k(k＝0，1，2，3)为待估参数；D为GPS水准点到起算点的水平距离；sign(·)为符号函数；a₃D²sign(Δh)ln(1+Δh²)为计算区域内水准高程传递系统偏差的特性函数。式(5)称为四参数模型，构造附有四个参数的条件方程，写成矩阵形式(6)：

$\underset{r\×n}{A}\underset{n\×1}{V}+\underset{r\×4}{C}\underset{4\×1}{a}-\underset{r\×1}{\mathrm{\Δ\ζ}}=0---\left(6\right)$

按附有未知参数的条件平差法，解r个条件方程，即可求得4个未知数a_k(k＝0，1，2，3)和残差重力异常改正数{v_i}，进而得到残差重力异常平差值

；再利用Stokes公式由

计算地面高程异常改正数估值，并进行精度评定。

此后，通过步骤S1获得的重力地面高程异常、步骤S2获得的地面高程异常改正数以及GPS水准高程异常系统偏差获得实用似大地水准面高。即：本发明技术方案构建的实用似大地水准面模型为：实用似大地水准面高＝重力地面高程异常+GPS水准高程异常系统偏差+地面高程异常改正数。

其中，GPS水准高程异常系统偏差与计算区域相关，例如，我国西部高大山区的GPS水准高程系统偏差可到20～30厘米。

S3：根据所述重力观测数据和所述实用似大地水准面模型获得正常高。

本发明技术方案的GPS水准高程传递方法的核心是公式：正常高＝大地高-高程异常。其中，“正常高”是工程建设中需要使用的高程值；“大地高”的获取方式参见步骤S1的记载；“高程异常”即为步骤S2中的“实用似大地水准面高”。输入待定点的GPS高程观测，利用本发明技术方案所构建的实用似大地水准面模型，即可得到待定点的水准高程，从而实现GPS水准代替三角高程测量或者常规水准实现高程控制。

下面将对本发明技术方案的GPS水准高程传递方法的传递精度进行检核。

由待测区域的数字地面模型中选取18个离散点参与计算(每个点都具有经纬度、水准和重力值)。其中，15个点参加拟合，3个点参加外部检核。所选取的15个点的内符合精度为0.0043米(标准差)，其中最大值为0.0108米，最小值为-0.0083米。检核点的选取原则为在GPS水准点较密的区域抽取，并尽量接近于线路中线。例如，图3所示为本发明的一个实施例中待测区域GPS水准和重力点分布示意图，检核点在图中使用圆形符号标识。

检核方法包括：输入检核点的GPS坐标(经度、纬度和大地高)，通过本发明技术方案所构建的实用似大地水准面模型，即可得到检核点的水准高；将计算得到的水准高与水准测量得到的水准进行比较。检核结果见表2。

表2

表2中的精度要求是根据二等水准每公里全中误差精度要求结合实际水准测量的闭合环数和闭合环周长反算而得。如表2所示，经校核，本发明技术方案的GPS水准高程传递方法的传递精度已经达到能够代替常规水准测量的二等水准传递精度。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (6)

1.一种GPS水准高程传递方法，其特征在于，所述方法包括：

S1：根据待测区域的地面重力观测数据，结合区域数字地面模型DTM和全球重力场模型EGM，基于莫洛金斯基理论构建区域重力似大地水准面模型；

S2：利用GPS水准超定边值方法获得GPS大地高与水准高程的参考系统差异及其影响的高程异常改正数，并根据所述重力似大地水准面模型、所述高程异常改正数以及预置的GPS水准高程异常系统偏差构建实用似大地水准面模型；

S3：根据所述重力观测数据和所述实用似大地水准面模型获得待测区域待定点的正常高。

2.如权利要求1所述的GPS水准高程传递方法，其特征在于，步骤S1进一步包括：

S101：通过地面重力观测数据获得待测区域离散点的扰动重力；

S102：通过顾及地形影响的移去-恢复方法获得所述离散点的平均扰动重力；

S103：根据所述平均扰动重力和所述全球重力场模型EGM获得零阶剩余高程异常；

S104：根据所述零阶剩余高程异常和所述全球重力场模型获得莫洛金斯基零阶高程异常；

S105：对线性莫洛金斯基I/II项进行系统分解，并结合区域地形进行量化分析，获得高程异常地形改正；

S106：根据所述莫洛金斯基零阶高程异常和所述地面高程异常地形改正获得重力地面高程异常。

3.如权利要求2所述的GPS水准高程传递方法，其特征在于，所述重力似大地水准面模型为：重力地面高程异常＝莫洛金斯基零阶高程异常+高程异常地形改正。

4.如权利要求1所述的GPS水准高程传递方法，其特征在于，步骤S2中，所述GPS水准超定边值方法进一步包括：

S201：通过移去一恢复法，将“移去”并格网化的重力场元作为统计观测量，以Stokes公式或Hotine公式作为约束条件，建立GPS水准高程异常与格网化后的重力场元之间的条件方程；

S202：基于所述条件方程，通过条件平差法获得所述格网化的重力场元观测量的改正数的估值，并进一步获得格网重力场元平差值；

S203：根据重力场积分公式，由所述格网重力场元平差值计算获得地面高程异常改正数。

5.如权利要求3或4所述的GPS水准高程传递方法，其特征在于，所述实用似大地水准面模型为：实用似大地水准面高＝重力地面高程异常+GPS水准高程异常系统偏差+地面高程异常改正数。

6.如权利要求5所述的GPS水准高程传递方法，其特征在于，所述待定点的正常高＝大地高-实用似大地水准面高；其中，所述大地高通过所述重力观测数据获得。

Images