CN105842745A - 重力异常的小波域优化位变滤波分离方法 - Google Patents

Abstract

本发明公开一种重力异常的小波域优化位变滤波分离方法，其首先根据实测异常数据的区域特征和数据处理的需要，对数据进行分块评价，得到小波域位变滤波的转换函数，然后计算局部等效系数，并设计小波域优化位变滤波器，最后对重力异常进行滤波处理，实现重力异常分离。本发明提出的小波域优化位变滤波方法，具有空间变化滤波能力，在不同空间位置实现不同的滤波器特性，从而能实现局部数据频谱与全局数据频谱存在较大差异的重力异常分离问题。

Description

重力异常的小波域优化位变滤波分离方法

技术领域

本发明涉及一种小波域优化位变滤波方法，具体涉及一种用于重力异常的小波域优化位变滤波分离方法。

背景技术

在重力勘探中，将由于地下岩石、矿物密度分布不均匀所引起的重力变化，或地质体与围岩密度的差异引起的重力变化，称为重力异常。实测重力异常是不同埋深、不同规模、不同形态的所有密度不均匀地质体重力效应的叠加，若要由实测重力异常反演某个特定地质体，必须先从叠加的重力异常中分离出单纯由该目标地质体引起的异常，再用分离出的异常进行反演和解释。由于重力场具有固有的多解性，从叠加异常中分离出目标地质体的单一异常是困难的。重力异常的有效分离是目前重力资料处理解释中的难题之一。

重力异常分离的方法很多，现在常用的主要包括最小二乘平滑法、解析延拓法、趋势分析法、匹配滤波法、补偿圆滑滤波法、正则化滤波法、维纳滤波法、小波变换法、非线性滤波法、优化滤波法等。这些方法具有各自的特点与优势，但同时也存在各自的局限性。例如，维纳滤波法中，为了准确地分离出目标地质体引起的异常，需要先估计出目标地质体的频谱，再以此为依据设计维纳滤波器。针对这个问题，Gapta和Ramani采用径向谱分析方法，估计目标地质体的频谱；Pawlowski和Hansen利用钻孔数据等先验地质信息构建模型，估计目标地质体的频谱，此后，Pawlowski又提出利用格林等效层理论估计目标地质体频谱的方法；郭良辉等在优选延拓的理论基础上，进一步研究提出基于维纳滤波和格林等效层理论的优化滤波法，该方法分离重力异常不需要已知延拓高度，可对重力异常进行指定频段的低通带通和高通滤波。这些经典的频率域分离方法都是根据全局数据频谱特征设计滤波器进行滤波分离，对于局部数据频谱和全局数据频谱差异较大的重力异常分离，频率域分离方法的效果有待改进。熊光楚等曾针对线性滤波器局限性问题，研究并提出用于重磁异常解释的空间域位变滤波器。只是当时考虑的还是异常组合不太复杂的简单情况，还不适合大数据量及复杂异常的分离情况。

发明内容

有鉴于此，有必要研究一种能根据局部数据频谱特征进行针对性的位变滤波，且适用于大数据量和复杂异常的分离方法。

一种重力异常的小波域优化位变滤波分离方法，其特征在于，包括如下步骤：

S1、根据实测异常数据的区域特征和数据处理的需要，将实测异常数据进行分块，相邻分块间不重叠；

S2、对每块实测数据进行频率域优化滤波法的特征评价，设计从第j块异常数据中分离出第l层数据的优化滤波转换函数H_jl(r)；

S3、以优化滤波转换函数H_jl(r)作为小波域滤波器的转换函数，求得从第j块异常数据中分离出第l层数据的局部滤波等效系数向量v(b_j)＝(v(a₁,b_j),v(a₂,b_j),...,v(a_i,b_j),...v(a_k,b_j))；

S4、将每个分块异常数据计算得到的局部滤波等效系数向量v(b_j)根据其空间位置进行组合，得到全体重力异常的小波域位变滤波等效系数向量v(b)；

S5、对全体重力异常数据做小波变换，得到小波变换系数W(a,b)；

S6、根据小波域位变尺度滤波等效系数向量v(b)和小波变换系数W(a,b)，从而计算得到小波域优化位变滤波分离出的第l层的异常数据。

优选的，所述步骤S1中，在将实测异常数据进行分块之前，需要对数据做小波变换时频分析，确定各频率成分在空间的分布规律。

优选的，所述步骤S2包括以下子步骤：

S21、对第j块异常数据进行傅立叶变换，计算第j块异常数据的对数功率谱P_j(r)；

S22、分段线性拟合第j块异常数据的对数功率谱；

S23、用格林等效层对数功率谱拟合第j块异常数据实测对数功率谱P_j(r)；

S24、如果目标层与其他深度层场源信息互不相关，则从第j块异常数据中分离出第l层数据的优化滤波转换函数.

优选的，所述从第j块异常数据中分离出第l层数据的优化滤波转换函数如下:

$H_{jl}\left(r\right)=\frac{P_{jl}\left(r\right)}{P_j\left(r\right)}$

其中，P_j(r)是j块异常数据的功率谱密度函数，P_jl(r)是第j块异常数据第l层异常数据的功率谱密度函数。

优选的，步骤S22还包括以下子步骤：

S221、分析第j块异常数据的对数功率谱，确定频段数和各频段的频率范围，每一频段对应一个格林等效层；

S222、利用直线P_jl(r)＝K_jlr＋B_jl拟合每一频段异常数据对数功率谱，其中，K_jl,B_jl分别是拟合第j块异常数据的第l频段的直线的斜率和截距。

优选的，步骤S3中的所述局部滤波等效系数向量v(b_j)＝(v(a₁,b_j),v(a₂,b_j),...,v(a_i,b_j),...v(a_k,b_j))由以下方程组求得：

$\begin{array}{c}{\hat{H}}_{b_j}\left(\mathrm{\ω}\right)=H_{jl}\left(r\right)\\ \left[\begin{array}{cccc}{\hat{g}}^{*}\left(a_1{\mathrm{\ω}}_1\right)& {\hat{g}}^{*}\left(a_2{\mathrm{\ω}}_1\right)& ...& {\hat{g}}^{*}\left(a_k{\mathrm{\ω}}_1\right)\\ {\hat{g}}^{*}\left(a_1{\mathrm{\ω}}_2\right)& {\hat{g}}^{*}\left(a_2{\mathrm{\ω}}_2\right)& ...& {\hat{g}}^{*}\left(a_k{\mathrm{\ω}}_2\right)\\ ...& ...& ...& ...\\ {\hat{g}}^{*}\left(a_1{\mathrm{\ω}}_n\right)& {\hat{g}}^{*}\left(a_2{\mathrm{\ω}}_n\right)& ...& {\hat{g}}^{*}\left(a_k{\mathrm{\ω}}_n\right)\end{array}\right]\left[\begin{array}{c}v(a_1,b)\\ v(a_2,b)\\ ...\\ v(a_k,b)\end{array}\right]=\left[\begin{array}{c}{\hat{H}}_b\left({\mathrm{\ω}}_1\right)\\ {\hat{H}}_b\left({\mathrm{\ω}}_2\right)\\ ...\\ {\hat{H}}_b\left({\mathrm{\ω}}_n\right)\end{array}\right]\end{array}$

其中：为给定小波函数的频率域表达式；a₁、a₂....a_k为尺度因子；b为位置因子。

优选的，步骤S6中所述第l层的异常数据由以下方程组求得：

$s^{\′}\left(x\right)=\underset{i=1}{\overset{k}{\Σ}}v(a_i,b)W(a_i,b)$

其中：a_i为尺度因子；b为位置因子。

本发明所述重力异常的小波域优化位变滤波分离方法，通过在小波域滤波方法的基础上增加位变滤波功能，并采用优化滤波方法评价分块数据的局部频谱特性，使所述的方法具有空间变化滤波能力，在不同空间位置表现出不同的滤波器特性，从而适用于局部数据频谱与全局数据频谱存在较大差异的重力异常分离问题。且在局部数据频谱与全局数据频谱差异较大时，该方法相对于频率域优化滤波法具有一定优势且与反演解释成果吻合得更好。

附图说明

图1为本发明所述的重力异常的小波域优化位变滤波分离方法的流程图；

图2为图1中步骤S2的子流程图；

图3为图2中步骤S22的子流程图；

图4为组合模型重力异常及其剖面图；

其中，图4(a)为组合模型重力异常，图4(b)为组合模型剖面图；

图5为组合模型重力异常频率域优化滤波分离结果；

其中，图5(a)为低通滤波分离结果，图5(b)为带通滤波分离结果，图5(c)为高通滤波分离结果；

图6为组合模型重力异常小波域优化位变滤波分离结果；

其中，图6(a)为位变低通滤波分离结果，图6(b)为位变带通滤波分离结果，图6(c)为位变高通滤波分离结果；

图7为实测数据剖面重力异常；

图8为现有技术的实测剖面重力异常频率域优化滤波异常分离结果；

其中，图8(a)为现有技术中区域异常的低通滤波结果，图8(b)为现有技术中局部异常的带通滤波结果，图8(c)为现有技术中噪声和拟合误差的高通滤波结果；

图9为采用本发明所述重力异常的小波域优化位变滤波分离方法设计小波域优化位变滤波器得到的实测剖面重力异常分离结果；

其中，图9(a)为本发明得到的区域异常的低通滤波结果，图9(b)为本发明得到的局部异常的带通滤波结果，图9(c)为本发明得到的噪声和拟合误差的高通滤波结果；

图10为实测剖面重力异常反演解释成果图；

图11为重力异常数据小波变换时频分析图。

具体实施方式

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

为实施对重力异常的小波域优化位变滤波分离，其首先需要进行小波域滤波；小波域滤波是基于小波变换的滤波方法，简要阐述如下。

信号s(x)的连续小波变换定义为：

$W_{g|s}(a,b)=\frac{1}{\sqrt{a}}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}s\left(x\right){\tilde{g}}^{*}\left(\frac{x-b}{a}\right)dx---\left(1\right)$

其中，是母小波，是的复共轭，a是尺度因子，b是位置因子。

令

公式(1)可用卷积方式描述如下：

$W_{g|s}(a,b)={\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}s\left(x\right)g_a(b-x)dx=s\left(b\right)*g_a\left(b\right)---\left(3\right)$

对公式(3)两边做傅立叶变换，可得：

$\hat{W}(a,\mathrm{\ω})=\sqrt{a}{\hat{\tilde{g}}}^{*}(-a\mathrm{\ω})\hat{s}\left(\mathrm{\ω}\right)---\left(4\right)$

其中，是s(x)的频谱，是的频谱。

简记为：

$\hat{W}(a,\mathrm{\ω})={\hat{g}}^{*}\left(a\mathrm{\ω}\right)\hat{s}\left(\mathrm{\ω}\right)---\left(5\right)$

为实现滤波功能，将公式(5)两边同乘以向量v，v＝(v₁,v_2,…,v_k)，其第i个分量v_i是第i个尺度a_i的函数v_i＝v(a_i)(i＝1,2,…,k)，v_i称为滤波器的等效系数。可得：

$\underset{i=1}{\overset{k}{\Σ}}{v}_i{\hat{g}}^{*}\left(a_i\mathrm{\ω}\right)\hat{s}\left(\mathrm{\ω}\right)=\underset{i=1}{\overset{k}{\Σ}}{v}_i\hat{W}(a_i,\mathrm{\ω})---\left(6\right)$

令

则(6)式简写为：

$\hat{H}\left(\mathrm{\ω}\right)\hat{s}\left(\mathrm{\ω}\right)=\underset{i=1}{\overset{k}{\Σ}}{v}_i\hat{W}(a_i,\mathrm{\ω})---\left(8\right)$

这里定义为小波域滤波器转换函数。

对公式(8)两边做傅立叶反变换，可得到：

$s^{\′}\left(x\right)=\underset{i=1}{\overset{k}{\Σ}}{v}_{i}W(a_i,b)---\left(9\right)$

其中，F^-1(·)表示对函数做傅里叶反变换，s′(x)为原信号s(x)经过小波域滤波后的信号。

公式(9)表示，信号s'(x)可按如下步骤得到：先对原信号s(x)做小波变换，然后对小波变换系数W(a_i,b)加权求和，权重系数即为小波域滤波器等效系数v_i。

若能给定小波域滤波器的转换函数则由公式(7)可得到小波域滤波器等效系数v_i，即完成小波域滤波器的设计。上述工作是我们进一步实现小波域优化位变滤波的基础。

接着需要构建小波域位变滤波器，具体构建过程如下：

假设公式(6)中滤波器等效系数v_i不仅与尺度a_i有关，而且与空间位置b有关，即为函数v(a_i,b)。则公式(6)可写作：

$\underset{i=1}{\overset{k}{\Σ}}v(a_i,b){\hat{g}}^{*}\left(a_i\mathrm{\ω}\right)\hat{s}\left(\mathrm{\ω}\right)=\underset{i=1}{\overset{k}{\Σ}}v(a_i,b)\hat{W}(a_i,\mathrm{\ω})---\left(10\right)$

令

公式(10)简写为：

${\hat{H}}_b\left(\mathrm{\ω}\right)\hat{s}\left(\mathrm{\ω}\right)=\underset{i=1}{\overset{k}{\∑}}v(a_i,b)\hat{W}(a_i,\mathrm{\ω})$

这里定义为小波域位变滤波器转换函数，该转换函数与位置b有关。

对公式(12)两边做傅立叶反变换，可得到：

$s^{\′}\left(x\right)=\underset{i=1}{\overset{k}{\Σ}}v(a_i,b)W(a_i,b)---\left(13\right)$

其中，

将公式(11)展开，写成如下的方程组形式：

$\left[\begin{array}{cccc}{\hat{g}}^{*}\left(a_1{\mathrm{\ω}}_1\right)& {\hat{g}}^{*}\left(a_2{\mathrm{\ω}}_1\right)& ...& {\hat{g}}^{*}\left(a_k{\mathrm{\ω}}_1\right)\\ {\hat{g}}^{*}\left(a_1{\mathrm{\ω}}_2\right)& {\hat{g}}^{*}\left(a_2{\mathrm{\ω}}_2\right)& ...& {\hat{g}}^{*}\left(a_k{\mathrm{\ω}}_2\right)\\ ...& ...& ...& ...\\ {\hat{g}}^{*}\left(a_1{\mathrm{\ω}}_n\right)& {\hat{g}}^{*}\left(a_2{\mathrm{\ω}}_n\right)& ...& {\hat{g}}^{*}\left(a_k{\mathrm{\ω}}_n\right)\end{array}\right]\left[\begin{array}{c}v(a_1,b)\\ v(a_2,b)\\ ...\\ v(a_k,b)\end{array}\right]=\left[\begin{array}{c}{\hat{H}}_b\left({\mathrm{\ω}}_1\right)\\ {\hat{H}}_b\left({\mathrm{\ω}}_2\right)\\ ...\\ {\hat{H}}_b\left({\mathrm{\ω}}_n\right)\end{array}\right]---\left(14\right)$

给定小波函数的频率域表达式和小波域位变滤波器转换函数求解方程组(14)，可得到小波域位变滤波器等效系数v(a_i,b)。由此可以看出，设计小波域位变滤波器的关键是如何确定其位变转换函数

在实际重磁异常分离中，人为给定小波域位变滤波器的转换函数往往是困难的。为此，本文根据实测异常数据的区域特征和数据处理的需要，对数据进行分块评价，分析异常在不同区域的特征差异，以便实现基于差异特征的与位置相关的位变滤波。具体实现中，我们采用频率域优化滤波法对每块局部数据进行频谱分析，得到小波域位变滤波的转换函数并计算局部等效系数，进而设计适应全局数据中具有局部特征差异的小波域位变滤波器。

如图1所示，本发明所述重力异常的小波域优化位变滤波分离方法包括如下步骤：

S1、根据实测异常数据的区域特征和数据处理的需要，对数据做小波变换时频分析，确定各频率成分在空间的分布规律，并将实测异常数据进行分块，相邻分块间不重叠。

具体的，如图4所示，对于图4(a)所示的重力异常，小波变换时频分析结果如图11所示，从图11可以看出，以x＝3500m为界，左右两边的频率成分空间分布规律具有明显差异，故将重力异常数据分为x<＝3500和x>3500两块。

S2、对每块实测异常数据进行频率域优化滤波法的特征评价，设计从第j块异常数据中分离出第l层数据的优化滤波转换函数H_jl(r)。

S3、以优化滤波转换函数H_jl(r)作为小波域滤波器的转换函数，即代入方程组(14)，求解该方程组得到第j块异常数据中分离出第l层数据的局部滤波等效系数向量v(b_j)＝(v(a₁,b_j),v(a₂,b_j),...,v(a_i,b_j),...v(a_k,b_j))。

具体的，步骤S3中的所述局部滤波等效系数向量v(b_j)＝(v(a₁,b_j),v(a₂,b_j),...,v(a_i,b_j),...v(a_k,b_j))由以下方程组求得：

${\hat{H}}_{b_j}\left(\mathrm{\&omega;}\right)=H_{\mathrm{jl}}\left(r\right)$

$\left[\begin{array}{cccc}{\hat{g}}^{*}\left(a_1{\mathrm{\&omega;}}_1\right)& {\hat{g}}^{*}\left(a_2{\mathrm{\&omega;}}_1\right)& ...& {\hat{g}}^{*}\left(a_k{\mathrm{\&omega;}}_1\right)\\ {\hat{g}}^{*}\left(a_1{\mathrm{\&omega;}}_2\right)& {\hat{g}}^{*}\left(a_2{\mathrm{\&omega;}}_2\right)& ...& {\hat{g}}^{*}\left(a_k{\mathrm{\&omega;}}_2\right)\\ ...& ...& ...& ...\\ {\hat{g}}^{*}\left(a_1{\mathrm{\&omega;}}_n\right)& {\hat{g}}^{*}\left(a_2{\mathrm{\&omega;}}_n\right)& ...& {\hat{g}}^{*}\left(a_k{\mathrm{\&omega;}}_n\right)\end{array}\right]\left[\begin{array}{c}v(a_1,b)\\ v(a_2,b)\\ ...\\ v(a_k,b)\end{array}\right]=\left[\begin{array}{c}{\hat{H}}_b\left({\mathrm{\&omega;}}_1\right)\\ {\hat{H}}_b\left({\mathrm{\&omega;}}_2\right)\\ ...\\ {\hat{H}}_b\left({\mathrm{\&omega;}}_n\right)\end{array}\right]$

其中：为给定小波函数的频率域表达式；a₁、a₂....a_k为尺度因子；b为位置因子。

S4、将每个分块异常数据计算得到的局部滤波等效系数向量v(b_j)根据其空间位置进行组合，得到全体重力异常的小波域位变滤波等效系数向量v(b)。

具体的，若重力异常数据分为k块，则等效系数向量v(b)＝(v(b₁),v(b₂),...,v(b_j),...v(b_k))，其中，v(b_j)是第j块数据的局部等效系数向量。所以，简单按照下标的顺序，可由v(b_j)组合得到v(b)。

S5、对全体重力异常数据做小波变换，得到小波变换系数W(a,b)。

S6、根据小波域位变尺度滤波等效系数向量v(b)和小波变换系数W(a,b)，从而计算得到小波域优化位变滤波分离出的第l层的异常数据。

具体的，步骤S6中所述第l层的异常数据由以下方程组求得：

$s^{\&prime;}\left(x\right)=\underset{i=1}{\overset{k}{\&Sigma;}}v(a_i,b)W(a_i,b)$

其中：a_i为尺度因子；b为位置因子。

如图2所示，所述步骤S2包括以下子步骤：

S21、对第j块异常数据进行傅立叶变换，计算第j块异常数据的对数功率谱P_j(r)。

S22、分段线性拟合第j块异常数据的对数功率谱。

S23、用格林等效层对数功率谱拟合第j块异常数据实测对数功率谱P_j(r)。

S24、假设目标层与其他深度层场源信息互不相关，则可设计从第j块异常数据中分离出第l层数据的优化滤波转换函数。具体为

$H_{ji}\left(r\right)=\frac{P_{jl}\left(r\right)}{P_j\left(r\right)}---\left(15\right)$

其中，P_j(r)是j块异常数据的功率谱密度函数，P_jl(r)是第j块异常数据第l层异常数据的功率谱密度函数。

具体的，如图3所示，步骤S22还包括以下子步骤：

S221、分析第j块异常数据的对数功率谱，确定频段数和各频段的频率范围，每一频段对应一个格林等效层；

S222、利用直线P_jl(r)＝K_jlr+B_jl拟合每一频段异常数据对数功率谱，其中，K_jl,B_jl分别是拟合第j块异常数据的第l频段的直线的斜率和截距。

下面根据附图对本发明所述重力异常的小波域优化位变滤波分离的方法进行具体阐述：

图4为组合模型的重力异常，其中，图4(b)所示的模型长方体1～4顶深分别为30米、140米、180米、490米，矩形截面分别为65×65、138×138、385×100、1650×165，异常密度分别为0.3g/cm3、0.25g/cm3、0.3g/cm3、0.25g/cm3，模型对应的重力异常如图4(a)所示。

如图5所示的是组合模型重力异常的频率域优化滤波法分离结果。其中，低通滤波分离结果为深部区域异常，带通滤波分离结果为浅部局部异常，高通滤波分离结果为优化滤波中采用格林等效层近似实测数据所产生的拟合误差。从图5(a)可以看出，分离出的区域异常主要是长方体2、3、4的异常，同时含有少量的长方体1的异常；从图5(b)可以看出，分离出的局部异常主要是长方体1的异常，同时含有极少量的长方体3的异常。值得注意的是，在图5(a)所示的区域异常中，长方体3、4的异常叠加在一起，未能有效分离开。如果单从深浅场源分离的角度看，即将浅部场源1、3从相对较深的场源2、4中分离出来，这个效果也不好。因为比较而言，局部深的场源2(相对于场源1)与局部浅的场源3(相对于场源4)相比，场源2的频谱反而更高，常规的整体频谱分离方式，不可能将场源2归算为深部场源，以及将场源3归算为浅部场源，而小波位变滤波则具有这样的能力。

图6所示的是按照本发明所述的方法设计小波域优化位变滤波器，模型重力异常分离结果。其中，位变低通、带通滤波分离结果为深部区域异常和浅部局部异常，位变高通滤波分离结果为拟合误差，由优化滤波拟合误差、小波域位变滤波误差和分块频谱估计所引入误差的叠加而成。从图6所示的小波优化位变滤波分离结果可以看出，仅需采用一次小波域位变滤波，即可有效分离出四个长方体的异常。其中，分离出的区域异常包含长方体2、4的异常，局部异常包含长方体1、3的异常，与给定的模型符合性明显改进。

图7所示的是某一重力异常剖面实测资料，在地表测点间距0.25km，剖面长度为380km。

图8所示的是实测重力异常的频率域优化滤波法分离结果。从图8可以看出，区域异常形态简单，仅反映出该剖面异常具有两边高、中间低的大尺度(波长大于180km)特征；相应地，局部异常形态过于复杂，包含了波长从6km到180km的所有异常，局部异常分离结果不够理想。

图9所示的是按照本发明所述的方法设计小波域优化位变滤波器，实测剖面重力异常分离结果。从图9可以看出，小波域优化位变滤波分离结果中，区域异常信息丰富，不仅反映出该剖面异常所具有的两边高、中间低的大尺度特征，还包含较多的中等尺度异常；相应的，局部异常则分辨得更加清晰，仅含有小尺度的浅部场源异常。

图10所示的是实测剖面重力异常密度反演结果图。对比图8和图10则会发现，频率域优化滤波分离出的区域异常过于简单，与反演结果不能很好的对应，而局部异常中残留有较多区域异常，既包含深度小于10km的浅部场源异常，又包含深度10～50km的深部场源异常。对比图9和图10可以看出，小波域优化位变滤波分离出的区域异常与反演结果中深度大于10km的场源分布规律一致，局部异常与反演结果中深度小于10km的场源分布规律一致，小波域优化位变滤波的分离结果与反演解释结果吻合得相当好。

综上所述，本发明所述重力异常的小波域优化位变滤波分离方法，通过在小波域滤波方法的基础上增加位变滤波功能，并采用优化滤波方法评价分块数据的局部频谱特性，使所述的方法具有空间变化滤波能力，在不同空间位置表现出不同的滤波器特性，从而适用于局部数据频谱与全局数据频谱存在较大差异的重力异常分离问题。附图实验表明，在局部数据频谱与全局数据频谱差异较大时，该方法相对于频率域优化滤波法具有一定优势且与反演解释成果吻合得更好。

以上所述仅为本发明的较佳实施例，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种重力异常的小波域优化位变滤波分离方法，其特征在于，包括如下步骤：

S1、根据实测异常数据的区域特征和数据处理的需要，将实测异常数据进行分块，相邻分块间不重叠；

S2、对每块实测数据进行频率域优化滤波法的特征评价，设计从第j块异常数据中分离出第l层数据的优化滤波转换函数H_jl(r)；

S3、以优化滤波转换函数H_jl(r)作为小波域滤波器的转换函数，求得从第j块异常数据中分离出第l层数据的局部滤波等效系数向量v(b_j)＝(v(a₁,b_j),v(a₂,b_j),...,v(a_i,b_j),...v(a_k,b_j))；

S4、将每个分块异常数据计算得到的局部滤波等效系数向量v(b_j)根据其空间位置进行组合，得到全体重力异常的小波域位变滤波等效系数向量v(b)；

S5、对全体重力异常数据做小波变换，得到小波变换系数W(a,b)；

S6、根据小波域位变尺度滤波等效系数向量v(b)和小波变换系数W(a,b)，从而计算得到小波域优化位变滤波分离出的第l层的异常数据。

2.根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，所述步骤S1中，在将实测异常数据进行分块之前，需要对数据做小波变换时频分析，确定各频率成分在空间的分布规律。

3.根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，所述步骤S2包括以下子步骤：

S21、对第j块异常数据进行傅立叶变换，计算第j块异常数据的对数功率谱P_j(r)；

S22、分段线性拟合第j块异常数据的对数功率谱；

S23、用格林等效层对数功率谱拟合第j块异常数据实测对数功率谱P_j(r)；

S24、如果目标层与其他深度层场源信息互不相关，则从第j块异常数据中分离出第l层数据的优化滤波转换函数。

4.根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，所述从第j块异常数据中分离出第l层数据的优化滤波转换函数如下:

其中，P_j(r)是j块异常数据的功率谱密度函数，P_jl(r)是第j块异常数据第l层异常数据的功率谱密度函数。

5.根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，步骤S22还包括以下子步骤：

S221、分析第j块异常数据的对数功率谱，确定频段数和各频段的频率范围，每一频段对应一个格林等效层；

S222、利用直线P_jl(r)＝K_jlr+B_jl拟合每一频段异常数据对数功率谱，其中，K_jl,B_jl分别是拟合第j块异常数据的第l频段的直线的斜率和截距。

6.根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，步骤S3中的所述局部滤波等效系数向量v(b_j)＝(v(a₁,b_j),v(a₂,b_j),...,v(a_i,b_j),...v(a_k,b_j))由以下方程组求得：

其中：为给定小波函数的频率域表达式；a₁、a₂....a_k为尺度因子；b为位置因子。

7.根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特 征在于，步骤S6中所述第l层的异常数据由以下方程组求得：

其中：a_i为尺度因子；b为位置因子。