CN111427096A

CN111427096A - 全张量重力梯度仪数据质量评价和滤波处理方法

Info

Publication number: CN111427096A
Application number: CN202010257866.6A
Authority: CN
Inventors: 袁园
Original assignee: Second Institute of Oceanography MNR
Current assignee: Second Institute of Oceanography MNR
Priority date: 2020-04-03
Filing date: 2020-04-03
Publication date: 2020-07-17

Abstract

本发明提供了一种全张量重力梯度仪数据质量评价和滤波处理方法。该方法包括：通过奇偶网格分析法分别对全张量重力梯度仪测量的各个梯度分量进行噪声估计；根据由噪声估计结果确定的各重力梯度分量权值，以及拉普拉斯方程的傅里叶级数解，构建加权联合反演滤波方程；利用最优线性反演方法求解加权联合反演滤波方程中的傅里叶级数系数；利用计算得到的傅里叶级数系数，正演计算得到滤波结果。本发明提供的全张量重力梯度仪数据质量评价和滤波处理方法能够实现在去除全张量重力梯度仪测量数据的高频噪声的同时保留有效的高频有效重力梯度信息，显著提高滤波后的成果数据的精度和分辨率。

Description

全张量重力梯度仪数据质量评价和滤波处理方法

技术领域

本发明涉及重力梯度数据处理技术领域，特别是涉及一种全张量重力梯度仪数据质量评价和滤波处理方法。

背景技术

全张量重力梯度测量相对传统的重力测量的主要优势在于能同时获得六个不同方向的重力梯度分量信息，从而从不同方向反映地下同一地质目标。

各个重力梯度分量满足引力位所满足的拉普拉斯方程约束条件，具有内部一致性，且具有高分辨率、高精度等优势，所包含的高频信号成分能更好的描述小的异常特征。

目前，全张量重力梯度数据的质量评价主要利用重力梯度数据满足的拉普拉斯特性去掉实际的重力梯度信号，对残余的噪声进行整体估计。该方法认为各个重力梯度张量受噪声的干扰水平相同。然而，实际情况下，不同方向的重力梯度数据受平台的影响并不相同，其噪声水平也不尽相同。

在全张量重力梯度数据滤波处理过程中，去除高频干扰噪声的同时保留有效的高频重力梯度信号，且保持各梯度分量之间的内部一致性将非常有挑战。传统的低通滤波方法通过假设各重力梯度分量真实异常与噪声拥有不同的能谱特性，使用不同截止频率的低通滤波器分别对各梯度张量单独进行滤波处理。然而，不同梯度分量原始数据的信噪比差别较大，噪声频谱特征存在差别，对应低通滤波的截止频率存在差异，经滤波处理后的各梯度分量将不再满足拉普拉斯约束条件，及不满足内部一致性。

随着重力梯度测量技术的发展，多种联合多个梯度张量分量的滤波方法被提出，可有效地保留各梯度分量之间的内部一致性。最简单的方法是根据引力位与重力梯度各分量之间的微分、积分关系分别对各梯度分量积分产生引力位，通过平均引力位来加强信噪比，从而微分获取滤波后的梯度张量分量。除此之外，联合多个重力梯度分量的反演滤波方法也用于相应的数据处理，但都是基于平面网格数据，对于机载沿山区起伏飞行和船载随波浪起伏测量的数据，需事先进行调平和网格化处理，从而引入数据后处理误差，如调平误差和网格差值误差等。

因此，为避免产生后处理误差影响，同时保证处理后的重力梯度数据保持内部一致性，研究基于三维空间不规则离散采样点测量数据的直接反演滤波方法至关重要。

发明内容

本发明要解决的技术问题是提供一种全张量重力梯度仪数据质量评价和滤波处理方法，能够实现在去除全张量重力梯度仪测量数据的高频噪声的同时保留有效的高频有效重力梯度信息，显著提高滤波后的成果数据的精度和分辨率。

为解决上述技术问题，本发明提供了一种全张量重力梯度仪数据质量评价和滤波处理方法，所述方法包括：通过奇偶网格分析法分别对全张量重力梯度仪测量的各个梯度分量进行噪声估计；根据由噪声估计结果确定的各重力梯度分量权值，以及拉普拉斯方程的傅里叶级数解，构建加权联合反演滤波方程；利用最优线性反演方法求解加权联合反演滤波方程中的傅里叶级数系数；利用计算得到的傅里叶级数系数，正演计算得到滤波结果。

在一些实施方式中，通过奇偶网格分析法分别对全张量重力梯度仪测量的各个梯度分量进行噪声估计，分析其噪声水平，包括：奇偶网格分析法利用奇数测线和偶数测线分别对数据进行网格，利用两个独立网格数据的和与差来实现对噪声的估计。

在一些实施方式中，利用两个独立网格数据的和与差来实现对噪声的估计，包括：根据如下公式对噪声的均方根误差进行估计：

n_d(i,j)＝[s_o(i,j)+n_o(i,j)]-[s_e(i,j)+n_e(i,j)]＝n_o(i,j)-n_e(i,j)

其中，n_d(i,j)为奇偶网格差，s_o(i,j)为奇数网格信号成分，s_e(i,j)为偶数网格信号成分，n_o(i,j)为奇数网格噪声成分，n_e(i,j)为偶数网格噪声成分。

在一些实施方式中，加权联合反演滤波方程具有以下形式：

Ac＝d

其中，A为方程系数，c为要拟合的系数,d为测量的重力梯度分量值。

在一些实施方式中，利用最优线性反演方法求解加权联合反演滤波方程中的傅里叶级数系数，还包括：对目标函数添加正则化项。

在一些实施方式中，正则化项具有如下形式：

其中，μ为正则化参数。

在一些实施方式中，利用最优线性反演方法求解加权联合反演滤波方程中的傅里叶级数系数，还包括：根据噪声估计结果，定义不同重力梯度分量的方向加权矩阵。

在一些实施方式中，方向加权矩阵具有如下形式：

W＝diag(1/σ₁,…,1/σ₆)

其中，σ_i为第i个数据的标准偏差。

采用这样的设计后，本发明至少具有以下优点：

本发明可以实现在去除全张量重力梯度仪测量数据的高频噪声的同时保留有效的高频有效重力梯度信息，显著提高滤波后的成果数据的精度和分辨率，为我国自主获取可用于地质解释的成果梯度数据的关键，并为国内自主研发重力梯度仪能尽早开展动态环境试验和测量数据处理提供技术支撑和方法积累。

附图说明

上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。

图1是全张量重力梯度仪数据质量评价和滤波处理方法的流程图；

图2为实测全张量重力梯度数据；

图3为全张量重力梯度数据反演滤波后结果。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

参见图1，利用全张量重力梯度仪数据质量评价和滤波处理方法主要包含以下步骤：

步骤1：采用奇偶网格分析法分别对全张量重力梯度仪测量的各个梯度分量进行噪声估计，分析测量过程中各个重力梯度分量的噪声水平；

奇偶网格分析法利用奇数测线和偶数测线分别对数据进行网格，利用两个独立网格数据的和与差来实现对噪声的估计。本方法可认为奇测线网格和偶数测线网格的数据中反映地质体的重力梯度信号成分相同，不同部分为噪声影响。

在网格数据中，每一个点(i,j)的数据中：奇数网格的信号成分为：s_o(i,j)；偶数网格的信号成分为：s_e(i,j)；其中，s_o(i,j)＝s_e(i,j)。奇数网格的噪声成分为：n_o(i,j)；偶数网格的噪声成分为：n_e(i,j)。奇数网格中每一个点(i,j)的数据为s_o(i,j)+n_o(i,j)；偶数网格中每一个点(i,j)的数据为s_e(i,j)+n_e(i,j)。

其中，奇偶网格差的表达式为：

n_d(i,j)＝[s_o(i,j)+n_o(i,j)]-[s_e(i,j)+n_e(i,j)]＝n_o(i,j)-n_e(i,j) (1)

重力梯度数据的网格差n_d(i,j)的均方根误差N_d为：

奇偶网格和的表达式为：

网格数据和中的每一个点(i,j)的数据为

S(i,j)＝[s_o(i,j)+n_o(i,j)+s_e(i,j)+n_e(i,j)]/2 (3)

网格数据和中的噪声n_c(i,j)为：

n_c(i,j)＝[n_o(i,j)]/2+[n_e(i,j)]/2 (4)

重力梯度数据的网格和中噪声的均方根误差N_c为：

联合奇偶网格差中噪声的均方根误差N_d与奇偶网格和中噪声的均方根误差N_c，我们可以得

N_c＝N_d/2 (6)

因此，奇偶网格差中的噪声水平为奇偶噪声和中的噪声水平的2倍。

步骤2：根据引力位V在自由空间满足的拉普拉斯方程求解各重力梯度分量在笛卡尔坐标系下的傅里叶级数解。

全张量重力梯度数据为引力位V在x，y，z三个正交方向上的二阶偏导数，构成了重力梯度张量矩阵：

其中

且

所以，矩阵T为对称矩阵。

由位场理论可知，引力位V在自由空间满足拉普拉斯方程：

根据方程(7)和(8)，可以知道T中只含有五个独立的梯度张量成分T_xx、T_xz、T_yx、T_yy、T_yz。全张量重力梯度仪提供了T_xx、T_xz、T_yx、T_yy、T_yz和T_zz六个张量成分。

在笛卡尔坐标系下求解引力位V和重力梯度张量分量T_αβ的解分别为：

其中，引力位V与重力梯度张量分量T_αβ的傅里叶级数解的系数关系见表1。

表1引力位V和重力梯度T_αβ之间的系数关系

表1中的各个离散采样点的波数通过勘探区域内x和y正交方向上的最大波长λ_xmax和λ_ymax决定。最大波长λ_xmax和λ_ymax决定x和y方向的最小波数

和

引力位V与重力梯度张量分量T_αβ的傅里叶级数解中的正、余弦基函数和指数衰减项都为已知量，与采样点的坐标x，y，z相关。未知量为系数a_nm、b_nm、c_nm、d_nm。

步骤3：联合公式(10)和表1，利用重力梯度分量的傅里叶级数解在三维空间离散测量点建立加权联合反演滤波方程：

Ac＝d (11)

方程中，A为方程系数，由傅里叶级数展开的基函数和指数部分组成，只与采样点的坐标x，y，z和波数相关；c为要拟合的系数，a_nm、b_nm、c_nm、d_nm构成的；d为测量的重力梯度分量值。

步骤4：采用最优线性反演方法求解傅里叶级数系数a_nm、b_nm、c_nm、d_nm，再通过正演计算求得的梯度值T_αβ来拟合实测重力梯度张量数据，将拟合结果作为滤波后的重力梯度异常。

在步骤4中，为保证反演计算的稳定性，对目标函数添加正则化项：

其中μ为正则化参数。使用L曲线和GCV准则估计最优μ值使目标函数

最小，实现系数对a_nm、b_nm、c_nm、d_nm的最优估计。

在最优线性反演过程中，通过步骤1中分析的噪声水平，定义不同重力梯度分量的方向加权矩阵，发挥不同张量分量各自的方向优势来提高反演滤波结果的可信度。数据方向加权矩阵为对角矩阵：

W＝diag(1/σ₁,…,1/σ₆) (13)

其中σ_i为第i个数据的标准偏差。

采用美国Bell Geospace公司提供的路易斯安娜州Vinton Dome的全张量重力梯度实测数据来证明本发明内容的适用性，实测数据和拉普拉斯滤波结果见图2和图3。

本发明针对Vinton Dome地区的全张量重力梯度数据的反演滤波，滤波结果更能准确的与各个重力梯度成分的地球物理意义相对应，更好的反映了各个重力梯度张量成分的异常形态,证明了基于拉普拉斯方程的傅里叶级数解的多个重力梯度分量联合反演滤波方法的适用性和有效性。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。

Claims

1.一种全张量重力梯度仪数据质量评价和滤波处理方法，其特征在于，包括：

通过奇偶网格分析法分别对全张量重力梯度仪测量的各个梯度分量进行噪声估计；

根据由噪声估计结果确定的各重力梯度分量权值，以及拉普拉斯方程的傅里叶级数解，构建加权联合反演滤波方程；

利用最优线性反演方法求解加权联合反演滤波方程中的傅里叶级数系数；

利用计算得到的傅里叶级数系数，正演计算得到滤波结果。

2.根据权利要求1所述的全张量重力梯度仪数据质量评价和滤波处理方法，其特征在于，通过奇偶网格分析法分别对全张量重力梯度仪测量的各个梯度分量进行噪声估计，分析其噪声水平，包括：

奇偶网格分析法利用奇数测线和偶数测线分别对数据进行网格，利用两个独立网格数据的和与差来实现对噪声的估计。

3.根据权利要求2所述的全张量重力梯度仪数据质量评价和滤波处理方法，其特征在于，利用两个独立网格数据的和与差来实现对噪声的估计，包括：

根据如下公式对噪声的均方根误差进行估计：

n_d(i,j)＝[s_o(i,j)+n_o(i,j)]-[s_e(i,j)+n_e(i,j)]＝n_o(i,j)-n_e(i,j)

4.根据权利要求1所述的全张量重力梯度仪数据质量评价和滤波处理方法，其特征在于，加权联合反演滤波方程具有以下形式：

Ac＝d

5.根据权利要求4所述的全张量重力梯度仪数据质量评价和滤波处理方法，其特征在于，利用最优线性反演方法求解加权联合反演滤波方程中的傅里叶级数系数，还包括：

对目标函数添加正则化项。

6.根据权利要求5所述的全张量重力梯度仪数据质量评价和滤波处理方法，其特征在于，正则化项具有如下形式：

其中，μ为正则化参数。

7.根据权利要求4所述的全张量重力梯度仪数据质量评价和滤波处理方法，其特征在于，利用最优线性反演方法求解加权联合反演滤波方程中的傅里叶级数系统，还包括：

根据噪声估计结果，定义不同重力梯度分量的方向加权矩阵。

8.根据权利要求7所述的全张量重力梯度仪数据质量评价和滤波处理方法，其特征在于，方向加权矩阵具有如下形式：

W＝diag(1/σ₁,…,1/σ₆)

其中，σ_i为第i个数据的标准偏差。