重力场模型辅助反距离加权似大地水准面格网插值方法

技术领域

本发明属于测绘科学与技术领域，尤其涉及一种重力场模型辅助反距离加权似大地水准面格网插值方法。

背景技术

目前，在工程应用中，我国以正常高作为法定高程系统，似大地水准面是我国的高程基准面，正常高是地面点位至似大地水准面的距离。全球导航卫星系统GNSS(GlobalNavigation Satellite System)技术获得的高程信息是相对于椭球面的大地高，大地高转换为正常高才能直接应用于工程建设。高程异常是似大地水准面至参考椭球面的距离，大地高转换为正常高的关键问题是需要获得GNSS观测点的高程异常，GNSS观测的大地高结合该点的高程异常即可获得该点位的正常高。高精度、高分辨率的似大地水准面模型，可以看作一种测定正常高的参考框架，给出任一点的高程异常。如果GNSS大地高测定的精度是一定的，GNSS高程转换的关键是高程异常的精度。

完成精化的区域似大地水准面为格网数值模型，模型按照一定经纬度格网给出高程异常值，GNSS观测一般为离散的观测点，并不位于似大地水准面格网数值模型的结点上，使用的时候需要利用似大地水准面格网数值模型进行高程异常内插计算。似大地水准面格网插值计算就是基于邻近点之间高程异常存在相关性，由邻近格网结点的高程异常内插出待定点的高程异常。高程异常内插推估的方法很多，有反距离加权插值、线性插值和人工神经网络插值等方法，在似大地水准面格网插值的方法中，反距离加权法是根据距离插值点最近的四个格网结点的高程异常，以距离的倒数为权进行的加权平均，该方法以思路明晰、计算简单且易于实现成为似大地水准面格网插值的一种常用方法，但是该方法的缺点是受似大地水准面的趋势性的影响较大，致使插值结果精度较低。

因此，如何提高反距离加权似大地水准面格网插值的精度是似大地水准面应用中需要解决的问题。

通过上述分析，现有技术存在的问题及缺陷为：现有技术受似大地水准面的趋势性的影响较大，致使插值结果精度较低。

解决以上问题及缺陷的难度为：由于地球密度和质量分布不均匀，似大地水准面为不规则曲面，无法用精确的数学模型来表达，任何拟合内插方法都有误差，致使插值结果和真实值相比存在一定的偏差。

解决以上问题及缺陷的意义为：提高似大地水准面格网插值的精度，获得较高精度的高程异常值，可将GNSS观测信息中的大地高转变为工程中应用的正常高，即海拔高程，对于水准测量难以实施的地区利用GNSS观测信息确定正常高，代替费时费力的水准测量，改变工程测量中高程测量的传统模式，具有重要意义。

发明内容

针对现有技术存在的问题，本发明提供了一种重力场模型辅助反距离加权似大地水准面格网插值方法。

本发明是这样实现的，一种重力场模型辅助反距离加权似大地水准面格网插值方法，所述重力场模型辅助反距离加权似大地水准面格网插值方法，包括：

步骤一，搜索距离插值点最近的四个格网结点，在似大地水准面格网数值模型中搜索距离插值点最近的四个格网结点，并根据似大地水准面格网数值模型给出四个格网结点的大地坐标和似大地水准面高程异常；大地坐标用于计算模型高程异常。

步骤二，计算格网结点与插值点的模型高程异常，利用参考重力场模型位系数，根据插值点和格网结点的大地坐标计算插值点和四个格网结点的模型高程异常；模型高程异常用于去掉插值区域的趋势性。

步骤三，计算格网点的残差高程异常，格网结点的高程异常减去参考重力场模型计算的模型高程异常，得到格网结点的残差高程异常；残差高程异常较似大地水准面高程异常起伏小，适于进行拟合内插。

步骤四，反距离加权法计算插值点的残差高程异常，依据四个格网结点的残差高程异常，按照距离倒数定权的方法计算插值点的残差高程异常；

步骤五，计算插值点的高程异常，插值点的残差高程异常加上参考重力场位系数模型计算的插值点的模型高程异常，即可得到插值点的高程异常结果。

进一步，所述步骤一中，搜索距离插值点最近的四个格网结点，具体过程为：

利用反距离加权方法进行似大地水准面格网插值，离散的高程异常待定点的高程异常与似大地水准面格网数值模型中距离插值点较近的格网点是相关的，利用相关的格网点的高程异常值推估插值点的高程异常；反距离加权方法进行似大地水准面格网插值首先应该搜索距离插值点最近的四个格网结点；

以离散的高程异常待定点S为中心，在似大地水准面格网数值模型中选取与该点的大地经度差和大地纬度差分别小于预定格网间距数目的格网结点，预定数目的似大地水准面格网结点以小于一个格网间距为条件选取获得，得到距离插值点最近的四个格网结点G1、G2、G3和G4的大地坐标和高程异常；

S为插值点，其大地坐标为(B

按照以上搜索原则，可在似大地水准面格网中搜索出G1-G4共4个格网结点。

进一步，所述步骤二中，计算格网点与插值点的模型高程异常具体过程为：

根据广义布隆斯公式，地球表面上任一点A的高程异常ζ为：

式中：扰动位T

进一步，所述计算格网点与插值点的模型高程异常中采用2000国家大地坐标系椭球作为参考椭球，根据参考椭球参数和大地坐标可以计算椭球面的正常位U

求地面任一点A的重力位W

式中：ρ为矢径；θ为极距；λ为地心经度；fM为地球引力常数，

利用参考重力场位系数模型，计算插值点S和距离插值点最近的四个格网点G1、G2、G3和G4的模型高程异常Mζ

进一步，所述步骤三中，计算格网点的残差高程异常具体过程为：

四个格网点G1、G2、G3和G4的似大地水准面高程异常ζ

Δζ

似大地水准面是包括高程异常起伏全波段的连续的、光滑的、不规则的曲面，似大地水准面起伏结构是长波占优，短波分量相对较小；如果首先由格网点高程异常ζ减去模型高程异常Mζ，获得起伏较小的的残差高程异常。

进一步，所述步骤四中，反距离加权法计算插值点的残差高程异常具体过程为：

依据四个格网结点的残差高程异常Δζ

式中,p

进一步，所述格网结点和插值点都是地球椭球面上的点，地球椭球面上两点间距离可用大地线的长度代替，两点间大地线长度是地球椭球面上两点间的最短线，大地线长度可根据椭球面上两点的大地经纬度进行计算，设椭球面上两点的大地坐标分别为P

由P

σ＝arccos[sin u

P

根据式(9)计算的大地线长度即为椭球面上两点之间最短距离，作为反距离加权法插值方法定权的依据。

进一步，所述步骤五中，计算插值点的高程异常具体过程为：

插值点S的残差高程异常加上参考重力场位系数模型计算的模型高程异常Mζ

ζ

利用起伏度较小的残差高程异常场进行反距离加权插值，提高拟合内插的精度。

本发明的另一目的在于提供一种接收用户输入程序存储介质，所存储的计算机程序使电子设备执行所述重力场模型辅助反距离加权似大地水准面格网插值方法，包括下列步骤：

步骤一，搜索距离插值点最近的四个格网结点，在似大地水准面格网数值模型中搜索距离插值点最近的四个格网结点，并根据似大地水准面格网数值模型给出四个格网结点的大地坐标和高程异常；

步骤二，计算格网结点与插值点的模型高程异常，利用参考重力场模型位系数，根据插值点和格网结点的大地坐标计算插值点和四个格网结点的模型高程异常；

步骤三，计算格网点的残差高程异常，格网结点的高程异常减去参考重力场模型计算的模型高程异常，得到格网结点的残差高程异常；

步骤四，反距离加权法计算插值点的残差高程异常，依据四个格网结点的残差高程异常，按照距离倒数定权的方法计算插值点的残差高程异常；

步骤五，计算插值点的高程异常，插值点的残差高程异常加上参考重力场位系数模型计算的插值点的模型高程异常，即可得到插值点的高程异常结果。

本发明的另一目的在于提供一种存储在计算机可读介质上的计算机程序产品，包括计算机可读程序，供于电子装置上执行时，提供用户输入接口以实施所述重力场模型辅助反距离加权似大地水准面格网插值方法。

结合上述的所有技术方案，本发明所具备的优点及积极效果为：本发明采用重力场模型辅助，移去似大地水准面格网数值模型中起伏性大，带有趋势性的长波部分，获得起伏性较小的残差高程异常进行内插计算，可以大幅度提高高程异常内插结果的精度。本发明利用椭球面上的大地线长度作为两点间距离进行定权，并给出根据椭球面上两点的大地坐标计算大地线长度的计算方法，计算较为简洁。

本发明利用残差高程异常格网数值，采用反距离加权法计算插值点的残差高程异常，不需要解求复杂的误差方程组，计算容易实现。距离插值点近的格网点给予较大的权值，使插值结果更加可靠。在根据大地线长度定权的过程中，计算权值的分母中引入距离平滑因子，避免当插值点与格网点距离很小的情况下，权值接近无穷大。插值计算过程中，采用距离倒数二次方定权，较常规的用距离倒数法定权计算结果精度有所提升。

附图说明

图1是本发明实施例提供的重力场模型辅助反距离加权似大地水准面格网插值方法流程图。

图2是本发明实施例提供的选取距离插值点最近的格网点示意图。

图3是本发明实施例提供的某区域高程异常等值线示意图。

图4是本发明实施例提供的区域残差高程异常等值线示意图。

图5是本发明实施例提供的仿真实验等值线示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

针对现有技术存在的问题，本发明提供了一种重力场模型辅助反距离加权似大地水准面格网插值方法，下面结合附图对本发明作详细的描述。

本发明提供的重力场模型辅助反距离加权似大地水准面格网插值方法业内的普通技术人员还可以采用其他的步骤实施，图1的本发明提供的重力场模型辅助反距离加权似大地水准面格网插值方法仅仅是一个具体实施例而已。

如图1所示，本发明实施例提供的重力场模型辅助反距离加权似大地水准面格网插值方法，包括：

S101：在似大地水准面格网数值模型中搜索距离插值点最近的四个格网结点，根据似大地水准面格网数值模型给出四个格网结点的大地坐标和似大地水准面高程异常；

S102：利用参考重力场模型位系数，根据插值点和搜索的格网结点的大地坐标计算格网结点与插值点的模型高程异常；

S103：格网点的似大地水准面高程异常减去模型高程异常即可获得四个格网结点的残差高程异常；

S104：根据格网结点残差高程异常采用反距离加权法计算插值点的残差高程异常；插值点的残差高程异常加上该点的模型高程异常即得到插值点的高程异常最终插值结果。

本发明实施例提供的S101中，搜索距离插值点最近的四个格网结点，具体过程为：

利用反距离加权方法进行似大地水准面格网插值，认为离散的高程异常待定点(插值点)的高程异常与似大地水准面格网数值模型中距离插值点较近的格网点的高程异常是相关的，利用相关的格网点的高程异常值推估插值点的高程异常。反距离加权方法进行似大地水准面格网插值首先应该搜索距离插值点最近的四个格网结点。

以离散的高程异常待定点S为中心，在似大地水准面格网数值模型中选取与该点的大地经度差和大地纬度差分别小于预定格网间距数目的格网结点，本发明中，所述预定数目的似大地水准面格网结点以小于一个格网间距为条件选取获得，得到距离插值点最近的四个格网结点G1、G2、G3和G4的大地坐标和高程异常，如图2所示。

如图2所示，S为插值点，其大地坐标为(B

按照以上搜索原则，可在似大地水准面格网中搜索出G1-G4共4个格网结点。

本发明实施例提供的S102中，计算格网点与插值点的模型高程异常具体过程为：

根据广义布隆斯(Bruns)公式，地球表面上任一点A的高程异常ζ为：

式中：扰动位T

计算中，采用2000国家大地坐标系椭球作为参考椭球，根据参考椭球参数和大地坐标可以计算椭球面的正常位U

求地面任一点A的重力位W

式中：ρ为矢径；θ为极距；λ为地心经度；fM为地球引力常数，

表1中所列为EGM96重力场模型的前5阶次的位系数

利用参考重力场位系数模型，计算插值点S和距离插值点最近的四个格网点G1、G2、G3和G4的模型高程异常Mζ

本发明实施例提供的S103中，计算格网点的残差高程异常具体过程为：

四个格网点G1、G2、G3和G4的高程异常ζ

Δζ

似大地水准面是包括高程异常起伏全波段的连续的、光滑的、不规则的曲面，似大地水准面起伏结构是长波占优(相对幅度大于90％)，短波分量相对较小。如果首先由格网点高程异常ζ减去模型高程异常Mζ，获得起伏较小的的残差高程异常(这个过程称为移去)。图3为某区域高程异常等值线，包含似大地水准面全波段的起伏；图4为该区域残差高程异常等值线，总体起伏较小，更突出高程异常的局部细节特征。

本发明实施例提供的S104中，反距离加权法计算插值点的残差高程异常具体过程为：

依据四个格网结点的残差高程异常Δζ

式中,p

格网结点和插值点都是地球椭球面上的点，地球椭球面上两点间距离可用大地线的长度代替，两点间大地线长度是地球椭球面上两点间的最短线，大地线长度可根据椭球面上两点的大地经纬度进行计算，设椭球面上两点的大地坐标分别为P

由P

σ＝arccos[sin u

P

根据式(9)计算的大地线长度即为椭球面上两点之间最短距离，作为反距离加权法插值方法定权的依据。

本发明实施例提供的S105中，计算插值点的高程异常具体过程为：

插值点S的残差高程异常加上参考重力场位系数模型计算的模型高程异常Mζ

ζ

该方法实质是利用起伏度较小的残差高程异常场进行反距离加权插值，来提高拟合内插的精度。

下面结合仿真实验对本发明的技术方案作进一步的描述。

利用EGM2008重力场模型计算某区域(东经100°-107°，北纬28°-32°)2.5′似大地水准面格网数值模型，计算时EGM2008重力场模型应用到最大阶次。该区域高程异常起伏较大，呈无规律变化，等值线如图5所示：

在该区域随机选取100离散点，这100个点也利用EGM2008重力场模型计算模型高程异常，作为高程异常“真值”，进行插值结果的比对。利用似大地水准面格网数值模型对100个离散点利用不同方法进行高程异常的拟合内插，结果与模型计算值进行比较。

反距离加权插值高程异常的公式为：按照距离倒数定权的方法进行加权计算，距离插值点S近的点给予较大权值；其具体表达式为：

式中,p

表2无辅助的反距离加权插值结果与高程异常“真值”之差的统计值/m

从统计结果数据来看，在似大地水准面插值计算中取权的二次方较直接取一次方精度更高。为验证重力场模型辅助插值的有效性，首先，采用EIGEN-CHAMP05S重力场模型辅助插值，辅助重力场模型应用到140阶次，重力场辅助插值结果与高程异常“真值”之差的统计值列入表3中：

表3重力场辅助插值结果与高程异常“真值”之差统计值/m

从表2和表3的统计结果数据来看，利用重力场模型辅助似大地水准面插值计算的高程异常结果比无模型辅助的插值结果精度有明显改善。但由于CHAMP卫星主要探测的地球重力场的长波部分，包含的中短波分量较少，EIGEN-CHAMP05S模型辅助高程异常内插精度提高幅度并不显著。再采用EGM_DIR_R2重力场模型辅助，EGM_DIR_R2重力场模型是采用GOCE卫星观测数据计算的模型，包含较多的重力场中短波分量。EGM_DIR_R2辅助重力场模型应用到140阶次，EGM_DIR_R2重力场辅助插值结果与高程异常“真值”之差的统计值列入表4中：

表4 EGM_DIR_R2重力场辅助插值结果与高程异常“真值”之差的统计值/m

从表2、表3和表4的统计结果数据来看，利用包含重力场中短波分量的重力场模型辅助似大地水准面插值计算的高程异常结果比无辅助的插值结果精度有非常明显改善。EGM_DIR_R2模型辅助反距离加权似大地水准面拟合内插较无模型辅助的插值结果精度提高40％多。

应当注意，本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的普通技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现，例如在诸如磁盘、CD或DVD-ROM的载体介质、诸如只读存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、诸如逻辑芯片、晶体管等的半导体、或者诸如现场可编程门阵列、可编程逻辑设备等的可编程硬件设备的硬件电路实现，也可以用由各种类型的处理器执行的软件实现，也可以由上述硬件电路和软件的结合例如固件来实现。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。