CN116538995A - 一种改进的似大地水准面模型精化方法 - Google Patents

Abstract

本发明公开了一种改进的似大地水准面模型精化方法，包括以下步骤：计算地面空间重力异常格网Δg；计算重力似大地水准面进行GNSS/水准高与重力似大地水准面融合，求解参数估计向量的迭代解。本发明基于Molodensky理论，研究快速实用的重力似大地水准面计算方法；采用抗差估计法，进行GNSS/水准与重力似大地水准面融合，本发明为利用GNSS技术精密测定控制点正常高展现了巨大潜力和良好前景。除某些工程测量必须的高精度水准测量外，繁重的水准测量在测绘作业中将逐步被GNSS测量所取代，可大大提高测绘生产效率，并产生巨大的经济效益。

Description

一种改进的似大地水准面模型精化方法

技术领域

本发明涉及物理大地测量技术领域，具体涉及一种改进的似大地水准面模型精化方法。

背景技术

大地水准面是与平均海平面最为密合的重力等位面，是大地测量基础设施中描述包括海洋在内的地球表面地形起伏的理想参考面，即高程的起算面。GNSS测定的大地高结合高精度似大地水准面模型可以快速获得精密海拔高程，据此将实现传统基于水准测量的地面标石高程基准向现代基于GNSS测量的数字高程基准转变。因此，不断精化似大地水准面是大地测量学的一项长期战略性任务。

似大地水准面计算是以地球重力场位理论为基础，通过解算相应的大地测量边值问题来确定的。Stokes边值理论是假定地球所有质量全部包含在大地水准面以内，以大地水准面为边界面，通过数学中的边值问题的解算。由于Stokes边值问题的前提是大地水准面以外不能有质量存在，所以在应用该理论计算大地水准面以前，必须先将大地水准面的外部质量消去。Molodensky边值问题从确定地球真正(自然)形状出发，引入似地形表面的概念，以似地形表面为边界面，直接以地面重力异常为边值确定高程异常，不需要重力归算和调整地壳质量，避开了对地壳密度的假设，克服了Stokes方法要求对大地水准面外部的质量进行归算的缺陷，建立新的更严密的地球形状和外部重力场理论。该方法虽然不要求已知地形质量密度，但需要求解扩展到整个地球表面的Molodensky积分方程，涉及边值条件方程线性化、球近似、斜向导数处理等一系列复杂数学问题，故而难用我们熟知的数学解析方法来表达，必须通过引入已知的似地球表面和正常重力位，将自由边值问题转化为固定边值问题，采用线性化方法建立线性边值条件。

在大地测量边值问题理论的基础上，以GNSS/水准确定的高精度、低分辨率的几何似大地水准面作为控制，将重力方法确定的高分辨率、低精度的重力似大地水准面与之拟合。高程拟合常用方法为最小二乘法。在观测值服从正态分布情况下，最小二乘平差是最优线性无偏估计。最小二乘法有极好的配赋误差的能力，但它不能抗拒粗差或异值的干扰，对粗差的存在非常敏感，可能会使拟合参数估值带有较大偏差。

发明内容

本发明的目的在于针对现有技术存在的上述问题，本发明提供一种改进的似大地水准面模型精化方法。

本发明的上述目的通过以下技术方案实现：

一种改进的似大地水准面模型精化方法，包括以下步骤：

步骤1、计算地面空间重力异常格网Δg；

步骤2、计算重力似大地水准面

步骤3、基于以下公式进行GNSS/水准高与重力似大地水准面融合：

v＝Ax-l

式中，H₁～H_n为第1～n个GNSS/水准点的GNSS/水准高，为第1～n个GNSS/水准点的重力似大地水准面，l为残差高程异常，l₁-l_n为第1～n个GNSS/水准点的残差高程异常，x为模型参数向量，A为系数矩阵，v为误差向量。

如上所述步骤1包括以下步骤：

步骤1.1、计算离散重力点的地面空间重力异常；

步骤1.2、计算离散重力点的地形改正和布格改正，得到离散重力点的布格重力异常：

Δg_B＝Δg_P+TC+BG

其中，Δg_P为离散重力点的空间重力异常，Δg_B为离散重力点的布格重力异常，TC为地形改正，BG为布格改正。

步骤1.3、对离散重力点的布格重力异常进行格网化，得格网布格重力异常；

步骤1.4、计算格网布格重力异常的格网中心处地形改正平均值和布格改正平均值，对格网布格重力异常进行恢复，求得地面空间重力异常格网。

如上所述步骤2中重力似大地水准面基于以下公式：

其中，ζ₀为零阶高程异常，ζ₁为一阶项G₁对高程异常的贡献。

如上所述零阶高程异常ζ₀基于以下公式：

其中，S(ψ)为以球面角距ψ为变量的Stokes核函数，R为地球平均半径，γ为似地形面上的正常重力，σ为单位球面，Δg_EGM为模型重力异常，ζ_EGM为高程异常。

如上所述一阶项G₁基于以下公式：

其中，h为流动点的正常高，h_P为计算点的正常高，l₀为计算点与流动点之间半径为R的球面距离。

如上所述一阶项G₁对高程异常的贡献ζ₁基于以下公式：

如上所述步骤3还包括求解参数估计向量的迭代解：

其中，k表示迭代次数，为第k次迭代的参数估计向量，i∈{H,N}，/>和/>分别为第k+1次迭代和第k次迭代的GNSS水准高程高的权矩阵，/>和/>分别为第k+1次迭代和第k次迭代的重力似大地水准面的权矩阵，/>为第k次迭代的GNSS水准高的改正量，/>为第k次迭代的重力大地水准面高的改正量，/>为单位权方差，I为单位矩阵，C为常数。

本发明相对于现有技术，具有以下有益效果：

Molodensky边值理论引入似地球表面与似大地水准面的概念，并利用地面上的重力异常求解地面上及外空扰动位，从而可以研究地球表面的真实形状。传统Molodensky边值理论计算重力似大地水准面的方法计算繁琐，难以精密计算。本发明针对Molodesnky边值理论的似大地水准面计算模型适用性问题，构建适合厘米级似大地水准面精化的Molodensky解实用简便的算法。

利用重力似大地水准面模型辅助内插GNSS/水准所得的几何大地水准面，常用方法为最小二乘法。在观测值服从正态分布情况下，最小二乘平差是最优线性无偏估计。最小二乘法有极好的配赋误差的能力，但它不能抗拒粗差或异值的干扰，对粗差的存在非常敏感，可能会使拟合参数估值带有较大偏差。抗差估计法可以在粗差不可避免的情形下，选择估计方法使未知量估值尽可能减免粗差的影响，提高大规模整体平差成果精度。

通过似大地水准面精化，利用GNSS技术结合高精度高分辨率似大地水准面模型，已成为高程测量的一种新方式，取代低等级水准测量，满足1：10000、1：5000、1：1000甚至更大比例尺测图的迫切需要。高精度似大地水准面的确定，对加快测绘工程的建设具有特别重要的科学意义，将产生巨大的社会效益和经济效益。

附图说明

图1为本发明的流程示意图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

一种改进的似大地水准面模型精化方法，包括以下步骤：

步骤1、地面空间重力异常格网计算：

步骤1.1、计算离散重力点的地面空间重力异常；

步骤1.2、利用SRTM3数字高程模型计算离散重力点的地形改正和布格改正，得到离散重力点的布格重力异常：

Δg_B＝Δg_P+TC+BG (1)

其中，Δg_P为离散重力点的空间重力异常，Δg_B为离散重力点的布格重力异常，TC为地形改正，BG为布格改正。

步骤1.3、采用Shepard格网化方法对离散重力点的布格重力异常进行格网化，得格网布格重力异常；

步骤1.4、利用SRTM3数字高程模型计算格网布格重力异常的格网中心处地形改正平均值和布格改正平均值，对格网布格重力异常进行恢复，求得地面空间重力异常格网。

步骤2、采用Molodensky理论计算重力似大地水准面：

步骤2.1、计算零阶高程异常ζ₀

由地面空间重力异常格网Δg与全球重力位模型计算的模型重力异常Δg_EGM的差异，按Stokes积分公式计算零阶剩余高程异常。零阶剩余高程异常加全球重力位模型计算的高程异常ζ_EGM得到零阶高程异常ζ₀。

式中，S(ψ)为以球面角距ψ为变量的Stokes核函数，R为地球平均半径，γ为似地形面上的正常重力，σ为单位球面。

步骤2.2、计算一阶项G₁：

式中，h为流动点的正常高，h_P为计算点的正常高，l₀为计算点与流动点之间半径为R的球面距离。

步骤2.3、由Stokes公式计算一阶项G₁对高程异常的贡献ζ₁：

步骤2.4、计算重力似大地水准面

步骤3、GNSS/水准高与重力似大地水准面融合：

由GNSS/水准高以及重力似大地水准面之间的拟合模型可表达为：

v＝Ax-l(6)

式中，H₁～H_n为第1～n个GNSS/水准点的GNSS/水准高，为第1～n个GNSS/水准点的重力似大地水准面，l为残差高程异常，l₁-l_n为第1～n个GNSS/水准点的残差高程异常，x为模型参数向量，A为系数矩阵，取决于拟合模型，一般选取多面函数法，v为误差向量。

通过最小二乘平差法，求出参数估计向量

以及GNSS水准高的改正量和重力大地水准面高的改正量/>

其中，P_H为GNSS水准高程高的权矩阵，P_N为重力似大地水准面的权矩阵，I为单位矩阵。反距离加权法是非规则分布点变成规则分布点常用的网格化方法。该方法的基本思想是离所估算的网格点距离越近的离散点对该网格点的影响越大，越远的离散点影响越小。

常用的最小二乘法对含粗差的观测量相当敏感，个别粗差就会对参数的估值产生较大的影响。抗差估计正是针对最小二乘法抗误差干扰性差这样一种缺陷而提出来的，其目的在于构造某种估计方法，使其对于模型误差，特别是对粗差具有较强的抵抗能力。抗差最小二乘估计的参数解式可表成迭代解：

其中k表示迭代次数，为第k次迭代的参数估计向量，P_i ^k(i∈{H,N},分别表示GNSS水准高和重力似大地水准面)为第k次迭代的等价权矩阵:

其中

式中，k表示迭代次数，为第k次迭代的参数估计向量，i∈{H,N}，/>和/>分别为第k+1次迭代和第k次迭代的GNSS水准高程高的权矩阵，/>和/>分别为第k+1次迭代和第k次迭代的重力似大地水准面的权矩阵，/>为第k次迭代的GNSS水准高的改正量，/>为第k次迭代的重力大地水准面高的改正量，/>为单位权方差，C为常数，通常取C＝3。

本文所描述的技术内容仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体技术方案做各种各样的修改或补充或采用类似的方式替代。但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (7)

1.一种改进的似大地水准面模型精化方法，其特征在于，包括以下步骤：

步骤1、计算地面空间重力异常格网Δg；

步骤2、计算重力似大地水准面

步骤3、基于以下公式进行GNSS/水准高与重力似大地水准面融合：

v＝Ax-l

式中，H₁～H_n为第1～n个GNSS/水准点的GNSS/水准高，为第1～n个GNSS/水准点的重力似大地水准面，l为残差高程异常，l₁-l_n为第1～n个GNSS/水准点的残差高程异常，x为模型参数向量，A为系数矩阵，v为误差向量。

2.根据权利要求1所述一种改进的似大地水准面模型精化方法，其特征在于，所述步骤1包括以下步骤：

步骤1.1、计算离散重力点的地面空间重力异常；

步骤1.2、计算离散重力点的地形改正和布格改正，得到离散重力点的布格重力异常：

Δg_B＝Δg_P+TC+BG

其中，Δg_P为离散重力点的空间重力异常，Δg_B为离散重力点的布格重力异常，TC为地形改正，BG为布格改正。

步骤1.3、对离散重力点的布格重力异常进行格网化，得格网布格重力异常；

步骤1.4、计算格网布格重力异常的格网中心处地形改正平均值和布格改正平均值，对格网布格重力异常进行恢复，求得地面空间重力异常格网。

3.根据权利要求2所述一种改进的似大地水准面模型精化方法，其特征在于，所述步骤2中重力似大地水准面基于以下公式：

其中，ζ₀为零阶高程异常，ζ₁为一阶项G₁对高程异常的贡献。

4.根据权利要求3所述一种改进的似大地水准面模型精化方法，其特征在于，所述零阶高程异常ζ₀基于以下公式：

其中，S(ψ)为以球面角距ψ为变量的Stokes核函数，R为地球平均半径，γ为似地形面上的正常重力，σ为单位球面，Δg_EGM为模型重力异常，ζ_EGM为高程异常。

5.根据权利要求4所述一种改进的似大地水准面模型精化方法，其特征在于，所述一阶项G₁基于以下公式：

其中，h为流动点的正常高，h_P为计算点的正常高，l₀为计算点与流动点之间半径为R的球面距离。

6.根据权利要求5所述一种改进的似大地水准面模型精化方法，其特征在于，所述一阶项G₁对高程异常的贡献ζ₁基于以下公式：

7.根据权利要求6所述一种改进的似大地水准面模型精化方法，其特征在于，所述步骤3还包括求解参数估计向量的迭代解：

其中，k表示迭代次数，为第k次迭代的参数估计向量，i∈{H,N}，/>和/>分别为第k+1次迭代和第k次迭代的GNSS水准高程高的权矩阵，/>和/>分别为第k+1次迭代和第k次迭代的重力似大地水准面的权矩阵，/>为第k次迭代的GNSS水准高的改正量，/>为第k次迭代的重力大地水准面高的改正量，/>为单位权方差，I为单位矩阵，C为常数。