局部大地水准面确定方法研究

北京中铁第五勘察设计院集团有限公司 崔立鲁

局部大地水准面确定方法研究

北京中铁第五勘察设计院集团有限公司 崔立鲁

频域输入输出法和最小二乘配置法是融合多种重力场信息精化局部大地水准面的两种方法。本文，笔者主要介绍了这两种融合方法的基本原理，推导了相关的实用计算公式，并在理论上比较了两者之间的异同。以融合GPS水准和重力异常数据精化局部大地水准面为例，通过实测数据计算分析了两种方法在精化局部重力场方面的精度情况。

一、背景

大地水准面或似大地水准面是获取地理空间信息的高程基准面。高精度局部大地水准面将为测绘学、地球物理、地球动力学及海洋学等相关学科的发展和应用提供重要的基础地球空间信息。我国似大地水准面的精化工作虽然取得了阶段性的进展，但是与欧美国家相比还有较大的差距，同时也无法完全满足我国目前测绘生产的需要。

二、计算方法

1.频域输入输出法。根据输入输出端口的不同，频域输入输出法可以分为多种模型，常见的由单输入单输出、双输入单输出、多输入单输出和多输入多输出等模型。本文，笔者以双输入单输出模型为例，讨论频域输入输出法的原理和计算模型。双输入单输出模式（Double-input Single-output System）指输入信号为两种不同类型的信号，而输出信号为一种类型的信号。双输入单输出系统（有噪声）模型如图1所示。

在考虑噪声的情况下，图1中，n1，n2分别为两种输入信号的噪声，h′，Δg′分别表示带有误差的输入信号，即大地水准面高和重力异常。

式（1）、（2）、（3）中，F表示傅立叶变换，F{y}=Y，F{b1}+B1，F{h′}=H′，F{b2}=B2，F{Δg′}=ΔG， F{e}=E。

将（3）式改写为：

式（5）中，E*表示为E的共轭量。其中Pee表示功率谱密度（即PSD），计算公式为：

式（6）中，Cxy和Pxy分别表示信号x，y之间的协方差和功率谱密度。

在Pee=m in准则下，可求得系数B1，B2。

若已知相关变量的功率谱密度，则联合（3）、（7）、（8）式即可计算该系统的输出信号。

由公式（7）、（8）式可知，如需计算输出信号必须知道输入信号的功率谱密度，本文计算功率谱密度的公式如下。

式（9）、（10）中，f1（i, j）和f2（k＋i,l＋j）表示2个不同格网观测值，x， y分别表示格网数据的横纵方向和分别为观测值f1（i, j）和f2（k＋i,l＋j）的数学期望，M和N分别为数据点沿着方向上的个数，k和l分别为沿着x和y方向上的距离值。根据式（9）、（10）可计算出相应的功率谱密度。

2.最小二乘配置法。在物理大地测量学中，把用最小二乘预估来确定重力场的方法称为最小二乘配置法。在重力场逼近中，最小二乘配置法的最大优点是可以同时利用多种不同类型的重力观测值来确定地球外部重力场，运用该方法的关键是要获得一个符合实际情况的关于观测值的协方差函数。

根据Moritz H.的理论推导，最小二乘预估公式为：

式（11）中，s表示推估点的值，l为观测值的数据，Csl为推估值与观测值的互相关协方差，Cll为观测值之间的协方差。

在重力数据融合中，当观测值为大地水准面高和重力异常，推估值为任意值时，其相应公式为：

式（14）中，Cyh和CyΔg分别表示推估值与大地水准面高和重力异常的互协方差函数，Chh，ChΔg，CΔgh和CΔgΔg分别表示大地水准面高与重力异常的自协方差函数及其互协方差函数，Dhh、DΔgΔg分别表示大地水准面高与重力异常的噪声自协方差函数，其中ChΔg=ChΔgT。 N、 Δg′分别表示推估信号大地水准面高和重力异常，而Cee（N）和Cee（Δg′）分别表示大地水准面高和重力异常推估信号噪声协方差矩阵。

由最小二乘配置法原理可知，获得与观测值相关的协方差函数是该方法计算的关键。本文主要采用格网数据，其协方差计算公式如下。

式（15）、（16）中，f*是观测值f与该地区观测值的平均值之差，C（S）横和 C（S）纵分别表示横向和纵向协方差，分别按下式计算。

由（15）和（17）式得到离散的协方差值，再由多项式函数拟合出经验协方差函数p（x）：

式（18）中， p（x）为函数值；（p1，p2，··，pn，pn＋1）为函数的参数值，即需拟合的参数；x 为函数自变量。

3.两种方法的比较。从两种方法的基本原理，可以看出频域输入输出法与最小二乘配置法具有一定的相似性，但也存在着一些区别。下面从四个不同的角度对两种方法进行比较。

（1）计算原理不同。频域输入输出法是基于频域的数据处理方法，其基本原理是根据不同数据的频谱特性（即功率谱密度）确定其在计算结果中所占的比重；而最小二乘配置则是一种空域的数据处理方法，该方法主要是根据不同数据的数学特性（即协方差）来确定相应的比重。两者的区别在于数学运算法则在频域和空域内是不同的，其中空域的矩阵求逆到了频域就变成了矩阵相除，这样就极大地减少了计算所需要的软硬件资源，并节省了大量的时间，提高了计算效率。

（2）计算准则不同。频域输入输出法和最小二乘配置法的计算准则分别是 和 。根据频域输入输出法和最小二乘配置法原理，可知 和 在频域和空域内是相互对应的，分别表示的是输出信号噪声的功率谱密度和推估值的误差协方差。因此，两种方法的计算准则是一致的。

（3）数据处理类型不同。地球重力场实际上是一个各向异性场。由于理论上的假设，最小二乘配置法只能计算各向同性的观测数据，在计算前需要对数据进行中心化，而频域输入输出法不存在那样的假设，因此可以计算各向异性和各向同性的数据。而宁津生等分别对各向异性和各向同性数据恢复局部重力场，结果表明：利用各向异性的功率谱密度函数作为先验信息所确定的局部重力场精度要优于由各向同性功率谱密度函数所得到的精度。

（4）两种方法在计算前都需要知道观测数据的先验信息，即功率谱密度和协方差。在计算过程中，将根据不同观测数据之间的相互关系和自身的特性来确定各种类型数据在计算中所占的比重。而功率谱密度和协方差可以通过快速傅立叶变换（FFT）实现两者之间的转换。

三、数值计算与分析

实测数据为我国某地区53×101格网重力异常数据和65个GPS水准数据，其纬度范围为22°26′17.5″N- 22°55′00″N，经度为113°41′00″E- 114°39′55″E，格网的横、纵间隔分别为32.5″和35″。格网重力异常的标准差为1.600毫伽，GPS水准数据的标准差为0.010 m。另有29个大致均匀分布的GPS水准点作为外部检核点。

表1 融合结果误差

由表1可知，实测数据的计算结果，融合后的大地水准面内外符合精度均达到厘米级，满足目前大地水准面的精度要求。从实测数据计算结果可以看出：IOST的大地水准面高内符合精度为±0.020 2 m，而LSC的内符合精度为±0.064 4 m。IOST和LSC大地水准面高的外符合精度分别为±0.037 6 m和±0.042 1 m。融合结果为重力异常的计算精度也是如此。由此可知，频域输入输出法的融合精度要略优于最小二乘配置法。两种方法的计算精度均达到精化局部大地水准面的要求。