阅读说明：本技术 一种航空重力测量数据的弱信息提取方法 (Weak information extraction method for aviation gravity measurement data ) 是由 王冠鑫 罗锋 周锡华 熊盛青 王林飞 于 2019-10-28 设计创作，主要内容包括：一种航空重力测量数据的弱信息提取方法,包括以下步骤：(S1)对DGPS数据和加速度数据进行改正；(S2)对改正后的DGPS数据和加速度数据进行卡尔曼滤波处理；(S3)对步骤(S2)中得到的结果进行卡尔曼平滑处理,得到重力异常数据。本发明一是将传统的频率域滤波转换成时间域滤波,避免了截止频率的选择及信号的丢失；二是从飞机运动的角度进行真实状态的估计,可以直接得到重力异常信号,简化了航空重力数据处理流程,提高了运算效率。三是给出了全新的航空重力异常提取的数学模型,避免了DGPS数据和加速度数据产生错位,即使针对起伏飞行数据,也可精确地提取出重力异常信号。(A weak information extraction method of aviation gravity measurement data comprises the following steps: (S1) correcting the DGPS data and the acceleration data; (S2) performing kalman filtering processing on the corrected DGPS data and acceleration data; (S3) kalman smoothing is performed on the result obtained in the step (S2) to obtain gravity anomaly data. The invention firstly converts the traditional frequency domain filtering into time domain filtering, thus avoiding the selection of cut-off frequency and the loss of signals; and secondly, the real state is estimated from the angle of the aircraft motion, so that the gravity anomaly signal can be directly obtained, the aviation gravity data processing flow is simplified, and the operation efficiency is improved. And thirdly, a brand-new mathematical model for extracting the aviation gravity anomaly is provided, the DGPS data and the acceleration data are prevented from being staggered, and even if the fluctuating flight data are targeted, the gravity anomaly signal can be accurately extracted.)

一种航空重力测量数据的弱信息提取方法

技术领域

本发明属于航空重力勘探技术领域，涉及一种航空重力测量数据的弱信息提取方法。

背景技术

重力勘探是一种十分重要的地球物理勘探手段，主要包括地面重力测量和航空重力测量，当遇到复杂山地、沼泽、海洋等人工无法通过的情况或需要进行快速扫面测量时，通常要采用后者对重力场数据进行采集。航空重力测量以飞机作为搭载平台，利用重力、定位传感器组合系统测量获得自由空间的航空重力原始测量数据；然后对所获得的原始测量数据进行解算，得到空间重力异常数据(即弱信息)；最终以此为基础，推演出有用的地球物理信息。实际的飞行测量过程中，由于飞机同时受到自身震动、起伏及转弯飞行、气流作用等多方面影响，无论是采用捷联式重力仪、矢量重力仪或是平台式重力仪进行飞行测量，其所获得的原始重力信号中均含有大量的噪声，噪信比高达几千甚至上万级。而对于航空重力测量结果，工程勘探工作要求重复线内符合精度达到0.8mGal，对于新研制的国产重力仪要求重复线内符合精度达到0.6mGal，这就要求数据解算技术达到非常高的水准。另一方面，重力仪测得的垂向加速度信息中包含了重力场信息、飞机的垂向加速度以及高频噪声。采用传统的频域滤波器(例如当前普遍采用的fir100s滤波)，一是无法去除原始测量数据中的包含的飞机垂向加速度，需要在解算过程中对其单独处理，二是无法确定准确的截止频率，只能大概提取原始测量数据的低频分量，而将高频分量中含有的有效信息一并滤除，导致后期数据解算与应用效果变差。因此，如何从原始测量数据中获取满足可用精度的重力异常数据，仍旧是制约航空重力勘探领域深入发展的一个关键技术难题。

发明内容

本发明的目的在于克服现有的技术缺陷，提供一种数据处理过程简洁、信息提取准确度与精度高、计算过程占用硬件资源少的航空重力测量数据弱信息提取方法。

为实现上述目的，本发明采用了如下技术方案：

一种航空重力测量数据的弱信息提取方法，包括以下步骤：

(S1)对DGPS数据和加速度数据进行改正；

(S2)对改正后的DGPS数据和加速度数据进行卡尔曼滤波处理；

(S3)对步骤(S2)中得到的结果进行卡尔曼平滑处理，得到重力异常数据。

步骤(S1)中所述对DGPS数据和加速度数据进行改正的方法，包括如下方法：偏心改正、厄特沃斯改正、正常场改正、零漂改正、水平加速度改正。

进一步，步骤(S2)中所述卡尔曼滤波处理的方法如下：

(S2.1)构建航空重力异常的数学模型：

其中，Δg为重力异常值，f_Σ为各项改正后垂向重力值，包括正常场改正、高度改正、水平加速度改正、厄缶改正、零漂改正、基点改正等，v_u为飞机的垂向加速度，q_Σ为各类噪声总和；

(S2.2)构建卡尔曼滤波状态方程：

式中，f₁'、f₂'、f₃'为两个水平加速度及垂向加速度的测量值；重力仪测得的垂向加速度误差可分为低频分量δf_T和高频分量δf，其中δf_T是由于低频扰动所造成的，例如飞机俯仰、横滚运动造成的角加速度；KF₁、KF₂是材料疲劳和温度变化造成的平台安装误差角的变化量；

(S2.3)构建卡尔曼滤波测量方程：

h'＝h+δh (3)

式中，h'为实际测量高度，h为真实高度，δh为实际测量高度的误差；

(S2.4)将步骤(S2.2)建立的状态方程和步骤(S2.3)建立的测量方程代入卡尔曼滤波公式(4)和(5)中进行求解：

X_k＝Φ_k,k-1X_k-1+B_k-1u_k-1+Γ_k-1W_k-1 (4)

Z_k＝H_kX_k+V_k (5)

式中，Φ_k,k-1、B_k-1为常数矩阵；u_k-1为控制项；Γ_k-1为系统噪声驱动阵；W_k为系统激励噪声序列，；V_k为量测噪声序列，H_k为量测阵；其中，卡尔曼滤波的输入值为f₁'、f₂'、f₃'、f_Σ、h'，输出值为h、v_u、KF₁、KF₂、Δg、

Δg为航空重力勘探需要提取的重力异常值(即弱信息)；

进一步，将所述步骤(S2.4)中采用的卡尔曼滤波器状态方程，由

X_k＝Φ_k,k-1X_k-1+B_k-1u_k-1+Γ_k-1W_k-1 (4)

替换为：

X_k＝Φ_k,k-1X_k-1+B_ku_k+ΓW (6)

式中，X_k为被估计状态；Φ_k,k-1、B_k-1、B_k为常数矩阵；u_k-1、u_k为控制项；Γ_k-1、Γ为系统噪声驱动阵；W_k、W为系统激励噪声序列。

本发明一种航空重力测量数据弱信息提取方法，具有以下有益效果：一是将传统的频率域滤波转换成时间域滤波。航空重力数据属于位场数据，重力异常信号频率分布于整个频带，传统的频域率滤波(如现阶段工程上常用的fir100s滤波)只是提取了信号的低频成分，对高频的部分进行了舍弃，造成了有用信息的丢失。而本技术方案中所选用的卡尔曼滤波则是通过飞行中重力传感器的受力及运动情况来对其真实状态进行估计，从而提取重力异常信号，避免了截止频率的选择及信号的丢失。二是简化了航空重力数据处理流程，提高了运算效率。采用传统的频率域滤波无法去除飞机运动产生的垂向加速度，需要对飞机的垂向加速度单独处理。而本技术方案所采用的卡尔曼滤波器从飞机运动的角度进行真实状态的估计，可以直接得到重力异常信号，无需其他处理。三是对通用的带控制项卡尔曼滤波方程(参见《卡尔曼滤波与组合导航原理(第3版)》，第49页，公式(2.2.35))进行了改进，本发明给出了航空重力异常提取的数学模型，根据实际特点将前一时刻的控制项修改为当前时刻的控制项，避免了DGPS数据和加速度数据产生错位。改进后的算法即使针对起伏飞行中获取的测量数据，也可精确地提取出重力异常信号。

附图说明

图1是实施例1中一种航空重力测量数据的弱信息提取方法的流程示意图；

图2是实施例2中采用传统fir100s滤波器的处理结果；

图3是实施例2中采用实施例1中所述方法的处理结果；

图4是实施例2中采用实施例1中所述带控制项公式改进前、后起伏飞行数据的处理结果。

具体实施方式

以下结合附图，进一步说明本发明一种航空重力测量数据的弱信息提取方法的具体实施方式。本发明一种航空重力测量数据的弱信息提取方法不限于以下实施例的描述。

实施例1：

本实施例给出一种航空重力测量数据的弱信息提取方法的原理及步骤。

如图1所示，航空重力测量数据的弱信息提取方法主要包括各项重力改正、卡尔曼滤波及卡尔曼平滑三个步骤。

(S1)各项重力改正：本步骤所采用的方法为航空重力领域现有的方法，例如《航空重力勘探理论方法及应用》(熊盛青等，第七章)中所述，此处不再赘述。

(S2)卡尔曼滤波：本步骤提出的方法的原理如下：

设被估计状态X_k在t_k时刻受到系统噪声序列W_k-1和不确定输入项u_k-1的驱动，则驱动机理的状态方程为：

X_k＝Φ_k,k-1X_k-1+B_k-1u_k-1+Γ_k-1W_k-1 (4)

式中，Φ_k,k-1、B_k-1为常数矩阵；u_k-1为控制项；Γ_k-1为系统噪声驱动阵；W_k为系统激励噪声序列，其方差为Q_k；V_k为量测噪声序列，其方差为R_k。且W_k和V_k满足：

其中假设Q_k为非负定，R_k为正定。Z_k和X_k满足线性关系，量测方程为：

Z_k＝H_kX_k+V_k (5)

式中，H_k为量测阵。

在航空重力实际测量中，根据牛顿第二定律对测量载体的受力及飞行状态进行分析，可构建航空重力异常的数学模型为：

其中，Δg为重力异常值，f_Σ为各项改正后垂向重力值(其中包括正常场改正、高度改正、水平加速度改正、厄缶改正、零漂改正)，

为飞机的垂向加速度，q_Σ为各类噪声总和。

从式(1)可以看出每一时刻的重力异常Δg需由该时刻的改正后垂向重力值f_Σ、垂向加速度

及噪声q_Σ求出，因此在航空重力测量中，仍根据式(4)中由前一时刻的控制项u_k-1和噪声W_k-1对当前时刻的状态进行驱动显然是不合适的。针对此问题，本子课题将对(4)式进行修正，即由当前时刻的控制项u_k和噪声W_k对当前状态进行驱动：

X_k＝Φ_k,k-1X_k-1+B_ku_k+Γ_kW_k (8)

实际上，重力异常近似于一个平稳均匀的随机过程，可以认为根据(1)所构建的状态方程的系统误差为一个恒定的状态，则有：

X_k＝Φ_k,k-1X_k-1+B_ku_k+ΓW (5)

其中k＝1,2，...，n。

将式(4)改成式(5)之后，在飞机平稳平飞时，数据解算结果略有改善，在飞机起伏飞行时，数据解算结果改善明显。

航空重力测量与地面重力测量的差别之一是后者在开展测区测量前进行标定，而前者只在实验室进行标定。而在实际测量中，航空重力传感器会受到温度、湿度、材料疲劳度等一系列外在因素的影响，可能会使实验室标定的参数存在偏差。因此，需要在式(1)的基础上增加额外的考虑因素。公式(9)是利用测量误差关系建立的等式：

式中，f₁'、f₂'、f₃'为两个水平加速度及垂向加速度的测量值；f_z1、f_z2、f_z3为两个水平加速度及垂向加速度的真实值；KF₃为重力仪比例因子的变化量；δf₁、δf₂为水平加速度测量误差；重力仪测得的垂向加速度误差可分为低频分量δf_T和高频分量δf，δf_T是由于低频扰动造成的(如飞机俯仰、横滚运动造成的角加速度)，这部分δf是由于振动、噪声及外部干扰等因素造成的。

为了便于解算，需要利用力学关系(10)对公式(9)中的非已知值f_z1、f_z2、f_z3进行替换：

f_z3＝f₃-(α₁+KF₁)f_z2+(α₂+KF₂)f_z1 (10)

式中，α₁、α₂分别为绕x、y轴倾斜的平台误差角；KF₁、KF₂是由于材料疲劳和温度变化造成的平台安装误差角的变化量，在一次飞行中基本不变，可当做常数处理。联立公式(1)、(9)、(10)，即可得到状态方程：

重力异常是一个随机过程，可以用成型滤波器：

对其进行建模，对式(2)进行补充，进而将重力异常Δg作为被估计状态。这里

为滤波状态向量，F_g、Γ_g为常数矩阵，q_g为产生强度为Q_g的白噪声。

取为二维，则有：

联合式(2)、(12)，可得状态方程：

离散处理后得到被估计状态

本实施例给出的卡尔曼滤波状态方程充分考虑了飞机的飞行状态，材料的疲劳度等因素的影响，而传统频率域滤波仅从信号频带的角度对数据进行处理。

(S3)卡尔曼平滑处理：由于卡尔曼滤波达到收敛需要一定的时间，在算法达到稳态之前，会出现小段未收敛的值，本实施例采用RTS平滑技术对该小段未收敛的值进行处理，如《卡尔曼滤波与组合导航原理(第3版)》，第189页所述。即利用已经收敛的被估计值对未收敛的被估计值进行重新估计。

本实施例给出的一种航空重力测量数据的弱信息提取方法在保留更多有用的信息的同时，能够更好的去除噪声干扰，能够应对更为复杂的飞行情况。

实施例2：

本实施例给出一种采用实施例1所述方法对某海域实际飞行获取的航空重力测量数据进行处理的具体结果。

如图2和图3所示，分别是传统fir100s滤波器的和实施例1所述卡尔曼滤波器处理结果。飞机对同一测线进行了4次测量，得到了4条重复线数据。通过对比4条重复线数据的内符合精度可以看出，卡尔曼滤波的处理结果较传统fir100s滤波器更为平滑，重复线内符合精度为0.471mGal,而传统fir100s处理结果的重复线内符合精度为0.719mGal，明显使用卡尔曼滤波更能满足实际勘探以及新研制国产重力仪的需求。由此可见，采用卡尔曼滤波相比于传统fir100s，在处理航空重力测量数据方面具有更优的效果。

如图4所示，分别是带控制项公式改进前和改进后的平飞及起伏飞行数据处理结果。通过对比可以看出，在将卡尔曼滤波公式修改后，在飞机平稳飞行阶段处理效果略有改善，在飞机起伏飞行阶段，处理效果改善明显。由此可见，采用式(5)相比于式(4)，在进行复杂飞行的情况下具有更优的效果。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。