CN111366130A - 一种高程异常内插推估方法 - Google Patents

Abstract

本发明提供一种高程异常内插推估方法，包括：以离散高程异常待定点为中心，在似大地水准面格网数值模型中选取与该点的大地经度差和大地纬度差分别小于预定数目格网间距的格网结点的大地坐标和高程异常；将选取的拟合控制点和高程异常待定点的大地坐标进行高斯投影，获得点位的投影平面坐标；根据拟合控制点的投影平面坐标计算拟合区域的高程异常二次多项式拟合参数；根据区域拟合参数和高程异常待定点的平面坐标计算该点的高程异常。本发明顾及了似大地水准面的局部特征。并且插值精度高。此外，本发明选用二次多项式拟合法进行拟合，计算方法简洁，有利于程序实现。

Description

一种高程异常内插推估方法

技术领域

本发明涉及工程应用技术领域，尤其涉及一种在似大地水准面格网数值模型中的高程异常内插推估方法。

背景技术

在工程应用中，我国法定高程系统为正常高系统，正常高是以似大地水准面为基准定义的高程系统。而工程应用中，GPS(GPS,Global Positioning System)技术获得的高程信息是椭球面的大地高，因此，需要将椭球面的大地高转换为正常高才能直接应用于工程建设。高精度、高分辨率的似大地水准面数值模型，给出区域等经纬度格网的高程异常值，使用时需要对高程异常待定点进行高程异常内插计算，似大地水准面格网插值就是基于邻近高程异常点之间存在相关性，由邻近点的高程异常内插出待定点的高程异常，如果GPS大地高测定的精度是一定的，GPS高程转换的关键是高程异常内插的精度。如何选取有效的方法，在高程异常内插计算过程中尽量少损失精度是提高GPS高程转换精度的关键。

似大地水准面是一个连续的、光滑的、不规则的曲面，高程异常内插推估的方法很多，工程应用中常用的有反距离加权插值法、线性插值法、二次曲面拟合法和人工神经网络方法等。反距离加权插值的缺点是插值是在单个格网，曲面会发生变形，曲面是断开的，不连续、不光滑，这种差值方法不能顾及趋势性的影响；线性插值法的缺点插值计算是在单个格网中进行中，插值结果不是光滑的曲面，而是扭面，格网边缘存在明显的分割线，类似于折痕，而在结点处形成尖点，导致曲面不光滑；二次曲面拟合法是一种较常用的高程异常插值方法，但是似大地水准面格网数值模型一般是等经纬度的格网，在地球表面，相同的经度差在不同的纬度上距离并不相等，给插值计算带来一定的误差。人工神经网络方法理论研究较多，但缺点是计算复杂，编程实现困难，不利于实际应用。

发明内容

本发明的目的在于解决上述现有技术存在的缺陷，提供一种在似大地水准面格网数值模型中的高程异常内插推估方法。

一种高程异常内插推估方法，包括以下步骤：

(1)拟合控制点的选取，以离散高程异常待定点为中心，在似大地水准面格网数值模型中选取与该点的大地经度差和大地纬度差分别小于预定数目格网间距的格网结点的大地坐标和高程异常；

所述预定数目格网间距的格网结点以两个格网间距为条件选取获得；

(2)点位投影，将选取的拟合控制点和高程异常待定点的大地坐标进行高斯投影，获得点位的投影平面坐标；

(3)计算拟合参数，根据拟合控制点的投影平面坐标计算拟合区域的高程异常二次多项式拟合参数；

(4)计算待定点的高程异常，根据区域拟合参数和高程异常待定点的平面坐标计算该点的高程异常。

进一步地，如上所述的高程异常内插推估方法，所述步骤(1)包括：

设定离散高程异常待定点D的大地坐标为(B_D,L_D)，选取的拟合控制点大地坐标为(B_i,L_i)，似大地水准面数值模型经、纬度格网间距为Δβ，则点位选取依据如下算法：

按照以上拟合控制点选取原则，在似大地水准面格网中取出16个格网节点高程异常作为拟合控制点。

进一步地，如上所述的高程异常内插推估方法，所述步骤(2)包括：

选用高程异常待定点实际经度作为投影的中央子午线，并且以中央子午线和赤道的交点o作为坐标原点，以中央子午线的投影为纵坐标x轴，以赤道的投影为横坐标y轴,形成高斯平面直角坐标系。

进一步地，如上所述的高程异常内插推估方法，在点位投影过程中，地球椭球的长半轴：a＝6378137m；地球扁率：1/298.257 222 101，地球自转角速度ω＝7.292115×10^{- 5}rad/s；地心引力常数GM＝3.986004418×10¹⁴m³/s²；

点位投影就是根据椭球面上A点的大地坐标(B,L)，求该点在高斯投影平面上的直角坐标(x,y)，即

的坐标变换，精度精确至0.5mm时的公式为:

式中:X为子午线弧长；l为椭球面上P点的经度与中央子午线L₀的经度差：l＝L-L₀,P点在中央子午线之东,l为正，在西则为负；N为地球椭球的卯酉圈曲率半径,t和η为简记符号；式(2)中其它量的计算公式为：

t＝tanB

η²＝e′²cos²B

在式(2)中，设有子午线上两点O和H，O在赤道上，H的纬度为B，a为椭球长半轴，b为椭球短半轴,O和H间的子午线弧长X计算公式

式(3)中，参数的计算公式为：

进一步地，如上所述的高程异常内插推估方法，所述步骤(3)包括：

假定在高程异常待定点邻近区域的高程异常ζ与坐标(x,y)存在如下函数关系：

式中a₀,a₁,a₂,a₃,a₄,a₅为多项式系数,

利用拟合控制点的投影平面坐标和高程异常组成多项式拟合方程，16个拟合控制点组成16个方程的方程组，根据最小二乘原理求解6个多项式系数的值。

进一步地，如上所述的高程异常内插推估方法，所述步骤(4)包括：

D点高程异常根据经纬度由多项式函数可求得

有益效果：

(1)本申请在似大地水准面格网数值模型中选取与高程异常待估点的位置的大地经度差和大地纬度差分别小于预定数目的格网结点的大地坐标和高程异常，其中，所述预定数目的似大地水准面格网结点以两个格网间距为条件选取获得；使拟合面为光滑曲面。

(2)将高程异常待估点和格网结点的大地坐标投影到平面上，解决了相同的经度差在不同的纬度上距离不相等的问题，以点位间实际位置关系进行拟合，提高插值精度。

(3)按照选取的格网节点进行投影后的实际位置计算二次多项式拟合参数，计算方法简单，有利于编程实现。

附图说明

图1为本发明拟合控制点选取示意图；

图2为本发明平面直角坐标系示意图；

图3为本发明子午线弧长示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

似大地水准面格网插值就是基于邻近高程异常点之间存在相关性，由邻近点的高程异常内插推估待定点的高程异常，如何选取有效的方法，在高程异常内插计算过程中尽量少损失精度是提高GPS高程转换精度的关键。

本发明提供一种在似大地水准面格网数值模型中的高程异常内插推估方法，插值高程异常的选点投影拟合方法，主要步骤包括：(1)拟合控制点的选取，以离散高程异常待定点为中心，在似大地水准面格网数值模型中选取与该点的大地经度差和大地纬度差分别小于预定数目格网间距的格网结点的大地坐标和高程异常；(2)点位投影，将选取的拟合控制点和高程异常待定点的大地坐标进行高斯投影，获得点位的投影平面坐标，解决相同的经度差在不同的纬度上距离不相等的问题；(3)计算拟合参数，根据拟合控制点的投影平面坐标计算拟合区域的高程异常二次多项式拟合参数；(4)计算待定点的高程异常，再根据区域拟合参数和高程异常待定点的平面坐标计算该点的高程异常。具体步骤如下：

1、选取拟合控制点

以离散的高程异常待定点为中心，在似大地水准面格网数值模型中选取与该点的大地经度差和大地纬度差分别小于预定格网间距数目的格网结点的大地坐标和高程异常，本发明中，所述预定数目的似大地水准面格网结点以小于两个格网间距为条件选取获得；

如图1，设定离散高程异常待定点D的大地坐标为(B_D,L_D)，选取的拟合控制点大地坐标为(B_i,L_i)，似大地水准面数值模型经、纬度格网间距为Δβ，则点位选取依据如下算法：

按照以上拟合控制点选取原则，在似大地水准面格网中取出A1-A16共16个格网节点高程异常作为拟合控制点。

2、点位投影

似大地水准面模型以一定间距的经纬度格网数值给出，选出的拟合控制点包括点位的大地坐标和高程异常值。大地坐标是位于椭球面上的，相同的经度差在不同的纬度上点间距离是不相等的，将离散的高程异常待定点和选取的拟合控制点的大地坐标投影到平面上，才能得到各点位的实际位置关系，按点位实际位置进行拟合可提高拟合精度。

首先建立投影平面直角坐标系。如图2所示，为减少投影变形，选用高程异常待定点实际经度作为投影的中央子午线。中央子午线和赤道的投影都是直线，并且以中央子午线和赤道的交点o作为坐标原点，以中央子午线的投影为纵坐标x轴，以赤道的投影为横坐标y轴,形成高斯平面直角坐标系。

在点位投影过程中，地球椭球参数选用2000国家大地坐标系采用的地球椭球参数：长半轴：a＝6378137m；地球扁率：1/298.257 222 101(

a为椭球长半轴，b为椭球短半轴)地球自转角速度ω＝7.292115×10^-5rad/s；地心引力常数GM＝3.986004418×10¹⁴m³/s²。

点位投影就是根据椭球面上A点的大地坐标(B,L)，求该点在高斯投影平面上的直角坐标(x,y)，即

的坐标变换。精度精确至0.5mm时的公式为:

式中:X为子午线弧长；l为椭球面上P点的经度与中央子午线(L₀)的经度差：l＝L-L₀,P点在中央子午线之东,l为正，在西则为负,式(2)中其它量的计算公式为：

t＝tanB

η²＝e′²cos²B

在式(2)中，子午线弧长X的计算较复杂,如图3所示,设有子午线上两点O和H，O在赤道上，H的纬度为B，O和H间的子午线弧长X计算公式

式(3)中，参数的计算公式为：

3、计算拟合参数

似大地水准面是一个连续的、光滑的、不规则的曲面，计算拟合参数就是根据选取的拟合控制点建立数学函数逼近该区域的似大地水准面。本发明采用二次多项式法进行拟合计算。假定在高程异常待定点邻近区域的高程异常ζ与坐标(x,y)存在如下函数关系：

式中a₀,a₁,a₂,a₃,a₄,a₅为多项式系数。

利用拟合控制点的投影平面坐标和高程异常组成多项式拟合方程，16个拟合控制点组成16个方程的方程组，根据最小二乘原理求解6个多项式系数的值。

4、计算高程异常待定点的值

D点高程异常根据经纬度由多项式函数可求得。

这样获得的搞成一场顾及了似大地水准面起伏的趋势性且能有效减少拟合插值的模型误差，计算生成的高程异常具有叫高精度精度。D点的高程异常和在该点测量的GPS大地高可计算该点的正常高，满足工程建设应用。

相较于现有技术，本发明的在似大地水准面格网数值模型中插值高程异常的选点投影拟合法通过设置以离散的高程异常待定点为中心，在似大地水准面格网数值模型中选取与该点的大地经度差和大地纬度差分别小于预定数目格网间距的格网结点的大地坐标和高程异常，将选取的拟合控制点和高程异常待定点的大地坐标进行投影，获得各点的平面坐标，在利用拟合控制点的平面坐标和高程异常计算逼近区域内似大地水准面的二次多项式拟合参数，再根据拟合参数计算待定点的高程异常。本发明选取与插值点相近的格网节点，顾及了似大地水准面的局部特征。并且将高程异常待估点和格网结点的大地坐标投影到平面上，点位间实际位置关系进行拟合，从而提高了插值精度。此外，本发明选用二次多项式拟合法进行拟合，计算方法简洁，有利于程序实现。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (6)

1.一种高程异常内插推估方法，其特征在于，包括以下步骤：

(1)拟合控制点的选取，以离散高程异常待定点为中心，在似大地水准面格网数值模型中选取与该点的大地经度差和大地纬度差分别小于预定数目格网间距的格网结点的大地坐标和高程异常；

所述预定数目格网间距的格网结点以两个格网间距为条件选取获得；

(2)点位投影，将选取的拟合控制点和高程异常待定点的大地坐标进行高斯投影，获得点位的投影平面坐标；

(3)计算拟合参数，根据拟合控制点的投影平面坐标计算拟合区域的高程异常二次多项式拟合参数；

(4)计算待定点的高程异常，根据区域拟合参数和高程异常待定点的平面坐标计算该点的高程异常。

2.根据权利要求1所述的高程异常内插推估方法，其特征在于，所述步骤(1)包括：

设定离散高程异常待定点D的大地坐标为(B_D,L_D)，选取的拟合控制点大地坐标为(B_i,L_i)，似大地水准面数值模型经、纬度格网间距为Δβ，则点位选取依据如下算法：

按照以上拟合控制点选取原则，在似大地水准面格网中取出16个格网节点高程异常作为拟合控制点。

3.根据权利要求1所述的高程异常内插推估方法，其特征在于，所述步骤(2)包括：

选用高程异常待定点实际经度作为投影的中央子午线，并且以中央子午线和赤道的交点o作为坐标原点，以中央子午线的投影为纵坐标x轴，以赤道的投影为横坐标y轴,形成高斯平面直角坐标系。

4.根据权利要求3所述的高程异常内插推估方法，其特征在于，在点位投影过程中，地球椭球的长半轴：a＝6378137m；地球扁率：1/298.257 222 101，地球自转角速度ω＝7.292115×10^-5rad/s；地心引力常数GM＝3.986004418×10¹⁴m³/s²；

点位投影就是根据椭球面上A点的大地坐标(B,L)，求该点在高斯投影平面上的直角坐标(x,y)，即

的坐标变换，精度精确至0.5mm时的公式为:

式中:X为子午线弧长；l为椭球面上P点的经度与中央子午线L₀的经度差：l＝L-L₀,P点在中央子午线之东,l为正，在西则为负；N为地球椭球的卯酉圈曲率半径,t和η为简记符号,式(2)中其它量的计算公式为：

t＝tanB

η²＝e′²cos²B

在式(2)中，设有子午线上两点O和H，O在赤道上，H的纬度为B，a为椭球长半轴，b为椭球短半轴,O和H间的子午线弧长X计算公式

式(3)中，参数的计算公式为：

5.根据权利要求1所述的高程异常内插推估方法，其特征在于，所述步骤(3)包括：

假定在高程异常待定点邻近区域的高程异常ζ与坐标(x,y)存在如下函数关系：

式中a₀,a₁,a₂,a₃,a₄,a₅为多项式系数,

利用拟合控制点的投影平面坐标和高程异常组成多项式拟合方程，16个拟合控制点组成16个方程的方程组，根据最小二乘原理求解6个多项式系数的值。

6.根据权利要求1所述的高程异常内插推估方法，其特征在于，所述步骤(4)包括：

D点高程异常根据经纬度由多项式函数可求得

Images