CN113433596B - 一种基于空间域的重力梯度动态测量滤波方法 - Google Patents

Abstract

本发明涉及一种基于空间域的重力梯度动态测量滤波方法，包括如下步骤：步骤1、利用卫星导航信息将由与时间对应重力梯度原始信息通过插值的方法投影到与空间位置对应的重力梯度原始数据；步骤2、直接在空间域设计滤波器并对与空间位置对应的重力梯度原始数据进行滤波，并直接输出该测量结果。本发明将传统基于时间域滤波的方法改进为基于空间域滤波，消除因搭载平台速度波动造成时空不一致引起的测量误差，同时确保处理后的重力梯度数据各个位置的空间分辨率均为设定值，不会出现局部空间分辨率波动。

Description

一种基于空间域的重力梯度动态测量滤波方法

技术领域

本发明属于重力梯度仪动态测量时的数字滤波方法，具体涉及一种基于空间域的重力梯度动态测量滤波方法。

背景技术

本发明涉及一种重力梯度仪动态测量时的数字滤波方法，能够严格的按照测量规范要求的空间分辨率对重力梯度仪原始测量数据进行低通滤波，提高重力梯度动态测量精度。

背景技术：

重力梯度定义为重力加速度矢量的空间梯度，即重力位的二阶导数，表征重力矢量的空间变化率。在地理坐标系中，重力矢量

可以分解为x、y、z三个方向上的三个分量中，每一分量沿平行于坐标轴方向均有一个梯度。因此，重力梯度张量共有3×3个分量,如式(1)所示。

数学上，重力梯度表示为：

式中：

Γ——地球外部任意空间位置重力梯度张量矩阵；

——当前位置重力加速度矢量；

——当前位置矢量；

Γ_ij(i,j＝x,y,z)——重力梯度张量各个分量，表示重力分量g_i在j方向上的空间变化率。

在实施重力梯度信息动态测量时，由于重力梯度仪原始输出信号存在大量高频噪声，故需要对仪器原始输出信号进行低通滤波以提高重力梯度测量精度。但重力梯度信号是重力信号的空间微分，因此重力梯度信号比重力信号具有更高的空间分辨率。传统的滤波方式是直接对仪器输出低通滤波，但仪器输出是与时间相关的测量信息，而重力梯度信号是与空间相关的信息，二者的因变量不同，传统的滤波方法容易引入新的误差，同时会造成滤波后重力梯度信号局部空间分辨率的不统一；传统数据处理方法是通过约束载体在测量过程中搭载平台速度平稳性来确保时空一致，但搭载平台受种种因素速度波动不可避免。因此需要提出一种新的滤波方法，降低上述问题对重力梯度动态测量的影响。

发明内容

本发明的目的在于克服现有技术的不足之处，提供一种基于空间域的重力梯度动态测量滤波方法。

本发明的上述目的通过如下技术方案来实现:

一种基于空间域的重力梯度动态测量滤波方法，其特征在于，包括如下步骤：

步骤1、利用卫星导航信息将由与时间对应重力梯度原始信息通过插值的方法投影到与空间位置对应的重力梯度原始数据；

步骤2、直接在空间域设计滤波器并对与空间位置对应的重力梯度原始数据进行滤波，并直接输出该测量结果。

进一步的：步骤1包括如下步骤：

步骤1.1、根据重力梯度动态测量同步获得的卫导信息，卫导信息为“时间t，纬度Lat，经度Lon”序列；

步骤1.2、将在时间域的重力梯度测量数据“时间t，重力梯度原始测量值Γ序列”投影到空间域，得到“纬度Lat，经度Lon，重力梯度原始测量值Γ序列”；具体为：

1.2.1、根据测线方向为南北方向或东西方向选取主位置信息L，选取方法为：若测线方向为南北方向，则定义主位置信息为纬度信息，即L＝Lat；若测线方向为东西方向，则定义主位置信息为经度信息，即L＝Lon；

1.2.2、获得“主位置信息L，重力梯度原始测量值Γ序列”，具体的：

先基于主位置信息L构建等差的主位置信息D，其中序列D的第一个元素为主位置信息L的最小值，序列D的最后一个元素为主位置信息L的最大值，序列D为等差序列，其序列的公差为 d，即：

d＝D_i+1-D_i…………………………………(2)

式中：

d——等差序列D的公差；

D_i——序列D的第i个元素。

公差d可根据实际情况选取，拟设置为0.0005°。

然后再根据实测不等间隔的主位置信息序列L，与之对应的重力梯度原始测量值序列Γ和构建的等间隔的主位置信息序列D，通过分段插值的方法获得与序列D对应的重力梯度原始测量值序列

定义插值函数H为：

H(L_j)＝Γ_j,j＝0,1,2,···…………………………………(3)

式中：

H——插值函数；

L_j——序列L的第j个元素。

Γ_j——序列Γ的第j个元素。

而H(x)在每个子区间[L_j,L_j+1],(j＝0,1,2,···)上是线性函数，即：

由此得到与序列D对应的重力梯度原始测量值序列

具体公式为：

——通过插值获得与序列D对应的重力梯度原始测量值序列。

1.2.3、计算步骤2.2.2中公差d对应的空间位置间隔ds，具体方法为：

首选计算测线中心点位置主位置信息的地球曲率半径，定义中心测线中心点的纬度、经度和高度分别为Latc、Lonc和hc；当主位置信息为纬度时，地球曲率半径计算公式为：

式中：

R——本条测线地球曲率半径；

R_e——地球长半轴长度；

Latc——测线中心位置的纬度；

e——地球椭球偏心率，也称为椭球度。

当主位置信息为经度时，地球曲率半径计算公式为：

然后由地球曲率半径R计算空间位置间隔ds，具体公式为：

ds＝R·d……………………………………(8)

式中：

ds——测线主位置信息的空间位置间隔。

进一步的；步骤2包括：

步骤2.1、在空间域设计数字滤波器

设计中选择有限长单位脉冲响应(FIR)滤波器，滤波器等效截止频率为

等效采样频率为

步骤2.2、应用滤波器对重构的等位置间隔重力梯度原始数据

进行滤波，具体公式为：

式中：

——滤波后的重力梯度信息；

b——步骤2.1中设计的FIR低通滤波器；

——重构的等空间间隔重力梯度原始数据；

最终重力梯度动态测量结果为等间隔的测线主位置信息D和安装空间域滤波后的重力梯度数据

本发明具有的优点和积极效果：

本发明将传统基于时间域滤波的方法改进为基于空间域滤波，消除因搭载平台速度波动造成时空不一致引起的测量误差，同时确保处理后的重力梯度数据各个位置的空间分辨率均为设定值，不会出现局部空间分辨率波动。

附图说明

图1为重力梯度张量分量示意图。

具体实施方式

以下结合附图并通过实施例对本发明的结构作进一步说明。需要说明的是本实施例是叙述性的，而不是限定性的。

一种基于空间域的重力梯度动态测量滤波方法，其发明点为，包括如下步骤：

步骤1、利用卫星导航信息将由与时间对应重力梯度原始信息通过插值的方法投影到与空间位置对应的重力梯度原始数据。包括如下具体步骤：

步骤1.1、根据重力梯度动态测量同步获得的卫导信息，卫导信息为“时间t，纬度Lat，经度Lon”序列；

步骤1.2、将在时间域的重力梯度测量数据“时间t，重力梯度原始测量值Γ序列”投影到空间域，得到“纬度Lat，经度Lon，重力梯度原始测量值Γ序列”；具体为：

1.2.1、根据测线方向为南北方向或东西方向选取主位置信息L，选取方法为：若测线方向为南北方向，则定义主位置信息为纬度信息，即L＝Lat；若测线方向为东西方向，则定义主位置信息为经度信息，即L＝Lon；

1.2.2、获得“主位置信息L，重力梯度原始测量值Γ序列”，具体的：

其中主位置信息L是不等间隔的数据，无法直接将该数据进行数字滤波，因此构建等差的主位置信息D，其中序列D的第一个元素为主位置信息L的最小值，序列D的最后一个元素为主位置信息L的最大值，序列D为等差序列，其序列的公差为d，即：

d＝D_i+1-D_i…………………………………(2)

式中：

d——等差序列D的公差；

D_i——序列D的第i个元素。

公差d可根据实际情况选取，拟设置为0.0005°。

然后再根据实测不等间隔的主位置信息序列L，与之对应的重力梯度原始测量值序列Γ和构建的等间隔的主位置信息序列D，通过插值的方法获得与序列D对应的重力梯度原始测量值序列

插值的方法主要有Lagrange插值、Newton插值、Hermite插值、分段插值和三次样条插值等，可根据实际测量数据选取最优插值方法。本专利中选用分段插值的方法，定义插值函数H为：

H(L_j)＝Γ_j,j＝0,1,2,···…………………………………(3)

式中：

H——插值函数；

L_j——序列L的第j个元素。

Γ_j——序列Γ的第j个元素。

而H(x)在每个子区间[L_j,L_j+1],(j＝0,1,2,···)上是线性函数，即：

由此得到与序列D对应的重力梯度原始测量值序列

具体公式为：

——通过插值获得与序列D对应的重力梯度原始测量值序列。

1.2.3、计算步骤2.2.2中公差d对应的空间位置间隔ds，具体方法为：

首选计算测线中心点位置主位置信息的地球曲率半径，定义中心测线中心点的纬度、经度和高度分别为Latc、Lonc和hc；当主位置信息为纬度时，地球曲率半径计算公式为：

式中：

R——本条测线地球曲率半径；

R_e——地球长半轴长度；

Latc——测线中心位置的纬度；

e——地球椭球偏心率，也称为椭球度。

当主位置信息为经度时，地球曲率半径计算公式为：

然后由地球曲率半径R计算空间位置间隔ds，具体公式为：

ds＝R·d……………………………………(8)

式中：

ds——测线主位置信息的空间位置间隔。

步骤2、直接在空间域设计滤波器并对与空间位置对应的重力梯度原始数据进行滤波，并直接输出该测量结果。具体为：

步骤2.1、在空间域设计数字滤波器

为保证滤波器的线性相移，在重力梯度空间域滤波的设计中选择有限长单位脉冲响应 (FIR)滤波器，目前该原理数字滤波器的设计方法较为成熟，而本方法中对滤波器的设计方法不限，只需保证设计的滤波器等效截止频率为

等效采样频率为

步骤2.2、应用滤波器对重构的等位置间隔重力梯度原始数据

进行滤波，具体公式为：

式中：

——滤波后的重力梯度信息；

b——步骤2.1中设计的FIR低通滤波器；

——重构的等空间间隔重力梯度原始数据；

最终重力梯度动态测量结果为等间隔的测线主位置信息D和安装空间域滤波后的重力梯度数据

尽管为说明目的公开了本发明的实施例和附图，但是本领域的技术人员可以理解：在不脱离本发明及所附权利要求的精神范围内，各种替换、变化和修改都是可以的，因此，本发明的范围不局限于实施例和附图所公开的内容。

Claims (1)

1.一种基于空间域的重力梯度动态测量滤波方法，其特征在于，包括如下步骤：

步骤1、利用卫星导航信息将由与时间对应重力梯度原始信息通过插值的方法投影到与空间位置对应的重力梯度原始数据；

步骤2、直接在空间域设计滤波器并对与空间位置对应的重力梯度原始数据进行滤波，并直接输出该测量结果；

步骤1包括如下步骤：

步骤1.1、根据重力梯度动态测量同步获得的卫导信息，卫导信息为“时间t，纬度Lat，经度Lon”序列；

步骤1.2、将在时间域的重力梯度测量数据“时间t，重力梯度原始测量值Γ序列”投影到空间域，得到“纬度Lat，经度Lon，重力梯度原始测量值Γ序列”；具体为：

1.2.1、根据测线方向为南北方向或东西方向选取主位置信息L，选取方法为：若测线方向为南北方向，则定义主位置信息为纬度信息，即L＝Lat；若测线方向为东西方向，则定义主位置信息为经度信息，即L＝Lon；

1.2.2、获得“主位置信息L，重力梯度原始测量值Γ序列”，具体的：

先基于主位置信息L构建等差的主位置信息D，其中序列D的第一个元素为主位置信息L的最小值，序列D的最后一个元素为主位置信息L的最大值，序列D为等差序列，其序列的公差为d，即：

d＝D_i+1-D_i………………………………… (2)

式中：

d——等差序列D的公差；

D_i——序列D的第i个元素；

公差d可根据实际情况选取，拟设置为0.0005°；

然后再根据实测不等间隔的主位置信息序列L，与之对应的重力梯度原始测量值序列Γ和构建的等间隔的主位置信息序列D，通过分段插值的方法获得与序列D对应的重力梯度原始测量值序列

定义插值函数H为：

H(L_j)＝Γ_j,j＝0,1,2,…………………………………… (3)

式中：

H——插值函数；

L_j——序列L的第j个元素；

Γ_j——序列Γ的第j个元素；

而H(x)在每个子区间[L_j,L_j+1],(j＝0,1,2,…)上是线性函数，即：

由此得到与序列D对应的重力梯度原始测量值序列

具体公式为：

——通过插值获得与序列D对应的重力梯度原始测量值序列；

1.2.3、计算步骤2.2.2中公差d对应的空间位置间隔ds，具体方法为：

首选计算测线中心点位置主位置信息的地球曲率半径，定义中心测线中心点的纬度、经度和高度分别为Latc、Lonc和hc；当主位置信息为纬度时，地球曲率半径计算公式为：

式中：

R——本条测线地球曲率半径；

R_e——地球长半轴长度；

Latc——测线中心位置的纬度；

e——地球椭球偏心率，也称为椭球度；

当主位置信息为经度时，地球曲率半径计算公式为：

然后由地球曲率半径R计算空间位置间隔ds，具体公式为：

ds＝R·d…………………………………… (8)

式中：

ds——测线主位置信息的空间位置间隔；

步骤2包括：

步骤2.1、在空间域设计数字滤波器

设计中选择有限长单位脉冲响应(FIR)滤波器，滤波器等效截止频率为

等效采样频率为

步骤2.2、应用滤波器对重构的等位置间隔重力梯度原始数据

进行滤波，具体公式为：

式中：

——滤波后的重力梯度信息；

b——步骤2.1中设计的FIR低通滤波器；

——重构的等空间间隔重力梯度原始数据；

最终重力梯度动态测量结果为等间隔的测线主位置信息D和安装空间域滤波后的重力梯度数据

Images