重力异常的小波域优化位变滤波分离方法
【专利摘要】本发明公开一种重力异常的小波域优化位变滤波分离方法，其首先根据实测异常数据的区域特征和数据处理的需要，对数据进行分块评价，得到小波域位变滤波的转换函数，然后计算局部等效系数，并设计小波域优化位变滤波器，最后对重力异常进行滤波处理，实现重力异常分离。本发明提出的小波域优化位变滤波方法，具有空间变化滤波能力，在不同空间位置实现不同的滤波器特性，从而能实现局部数据频谱与全局数据频谱存在较大差异的重力异常分离问题。
【专利说明】
重力异常的小波域优化位变滤波分离方法
技术领域
[0001] 本发明涉及一种小波域优化位变滤波方法，具体涉及一种用于重力异常的小波域 优化位变滤波分离方法。
【背景技术】
[0002] 在重力勘探中，将由于地下岩石、矿物密度分布不均匀所引起的重力变化，或地质 体与围岩密度的差异引起的重力变化，称为重力异常。实测重力异常是不同埋深、不同规 模、不同形态的所有密度不均匀地质体重力效应的叠加，若要由实测重力异常反演某个特 定地质体，必须先从叠加的重力异常中分离出单纯由该目标地质体引起的异常，再用分离 出的异常进行反演和解释。由于重力场具有固有的多解性，从叠加异常中分离出目标地质 体的单一异常是困难的。重力异常的有效分离是目前重力资料处理解释中的难题之一。
[0003] 重力异常分离的方法很多，现在常用的主要包括最小二乘平滑法、解析延拓法、趋 势分析法、匹配滤波法、补偿圆滑滤波法、正则化滤波法、维纳滤波法、小波变换法、非线性 滤波法、优化滤波法等。这些方法具有各自的特点与优势，但同时也存在各自的局限性。例 如，维纳滤波法中，为了准确地分离出目标地质体引起的异常，需要先估计出目标地质体的 频谱，再以此为依据设计维纳滤波器。针对这个问题，Gapta和Ramani采用径向谱分析方法， 估计目标地质体的频谱;Pawlowski和Hansen利用钻孔数据等先验地质信息构建模型，估计 目标地质体的频谱，此后，Pawlowski又提出利用格林等效层理论估计目标地质体频谱的方 法;郭良辉等在优选延拓的理论基础上，进一步研究提出基于维纳滤波和格林等效层理论 的优化滤波法，该方法分离重力异常不需要已知延拓高度，可对重力异常进行指定频段的 低通带通和高通滤波。这些经典的频率域分离方法都是根据全局数据频谱特征设计滤波器 进行滤波分离，对于局部数据频谱和全局数据频谱差异较大的重力异常分离，频率域分离 方法的效果有待改进。熊光楚等曾针对线性滤波器局限性问题，研究并提出用于重磁异常 解释的空间域位变滤波器。只是当时考虑的还是异常组合不太复杂的简单情况，还不适合 大数据量及复杂异常的分离情况。
【发明内容】
[0004] 有鉴于此，有必要研究一种能根据局部数据频谱特征进行针对性的位变滤波，且 适用于大数据量和复杂异常的分离方法。
[0005] -种重力异常的小波域优化位变滤波分离方法，其特征在于，包括如下步骤：
[0006] S1、根据实测异常数据的区域特征和数据处理的需要，将实测异常数据进行分块， 相邻分块间不重叠；
[0007] S2、对每块实测数据进行频率域优化滤波法的特征评价，设计从第j块异常数据中 分离出第1层数据的优化滤波转换函数Η^( Γ);
[0008] S3、以优化滤波转换函数H^(r)作为小波域滤波器的转换函数，求得从第j块异常 数据中分离出第1层数据的局部滤波等效系数向量 ν(~) = (ν(&1，ι^)，ν(&2，ι^)，...， ν(&1， bj),...v(ak,bj))；
[0009] S4、将每个分块异常数据计算得到的局部滤波等效系数向量V(W)根据其空间位 置进行组合，得到全体重力异常的小波域位变滤波等效系数向量v(b);
[0010] S5、对全体重力异常数据做小波变换，得到小波变换系数W(a，b);
[0011] S6、根据小波域位变尺度滤波等效系数向量v(b)和小波变换系数W(a，b)，从而计 算得到小波域优化位变滤波分离出的第1层的异常数据。
[0012] 优选的，所述步骤Sl中，在将实测异常数据进行分块之前，需要对数据做小波变换 时频分析，确定各频率成分在空间的分布规律。
[0013] 优选的，所述步骤S2包括以下子步骤：
[0014] S21、对第j块异常数据进行傅立叶变换，计算第j块异常数据的对数功率谱h(r);
[0015] S22、分段线性拟合第j块异常数据的对数功率谱；
[0016] S23、用格林等效层对数功率谱拟合第j块异常数据实测对数功率谱h(r);
[0017] S24、如果目标层与其他深度层场源信息互不相关，则从第j块异常数据中分离出 第1层数据的优化滤波转换函数.
[0018] 优选的，所述从第j块异常数据中分离出第1层数据的优化滤波转换函数如下：
[0019]
[0020] 其中，Pj(r)是j块异常数据的功率谱密度函数，P^(r)是第j块异常数据第1层异常 数据的功率谱密度函数。
[0021] 优选的，步骤S22还包括以下子步骤：
[0022] S221、分析第j块异常数据的对数功率谱，确定频段数和各频段的频率范围，每一 频段对应一个格林等效层；
[0023] S222、利用直线Pji(r)=Kjir + Bji拟合每一频段异常数据对数功率谱，其中，Kji， Bjl分别是拟合第j块异常数据的第1频段的直线的斜率和截距。
[0024]优选的，步骤53中的所述局部滤波等效系数向量￥(13」）=(>(31，13」），￥(32，13」），...， v(ai，bj)，…v(ak，bj))由以下方程组求得：
[0025:
[0026] 其中：f为给定小波函数的频率域表达式;ai、a2. . . .ak为尺度因子;b为位置 因子。
[0027] 优选的，步骤S6中所述第1层的异常数据由以下方程组求得：
[0028]
[0029] 其中阳为尺度因子;b为位置因子。
[0030] 本发明所述重力异常的小波域优化位变滤波分离方法，通过在小波域滤波方法的 基础上增加位变滤波功能，并采用优化滤波方法评价分块数据的局部频谱特性，使所述的 方法具有空间变化滤波能力，在不同空间位置表现出不同的滤波器特性，从而适用于局部 数据频谱与全局数据频谱存在较大差异的重力异常分离问题。且在局部数据频谱与全局数 据频谱差异较大时，该方法相对于频率域优化滤波法具有一定优势且与反演解释成果吻合 得更好。
【附图说明】
[0031] 图1为本发明所述的重力异常的小波域优化位变滤波分离方法的流程图；
[0032]图2为图1中步骤S2的子流程图；
[0033]图3为图2中步骤S22的子流程图；
[0034]图4为组合模型重力异常及其剖面图；
[0035]其中，图4(a)为组合模型重力异常，图4(b)为组合模型剖面图；
[0036]图5为组合模型重力异常频率域优化滤波分离结果；
[0037]其中，图5(a)为低通滤波分离结果，图5(b)为带通滤波分离结果，图5(c)为高通滤 波分离结果；
[0038]图6为组合模型重力异常小波域优化位变滤波分离结果；
[0039] 其中，图6(a)为位变低通滤波分离结果，图6(b)为位变带通滤波分离结果，图6(c) 为位变高通滤波分离结果；
[0040] 图7为实测数据剖面重力异常；
[0041] 图8为现有技术的实测剖面重力异常频率域优化滤波异常分离结果；
[0042] 其中，图8(a)为现有技术中区域异常的低通滤波结果，图8(b)为现有技术中局部 异常的带通滤波结果，图8(c)为现有技术中噪声和拟合误差的高通滤波结果；
[0043]图9为采用本发明所述重力异常的小波域优化位变滤波分离方法设计小波域优化 位变滤波器得到的实测剖面重力异常分离结果；
[0044]其中，图9(a)为本发明得到的区域异常的低通滤波结果，图9(b)为本发明得到的 局部异常的带通滤波结果，图9(c)为本发明得到的噪声和拟合误差的高通滤波结果；
[0045]图10为实测剖面重力异常反演解释成果图；
[0046]图11为重力异常数据小波变换时频分析图。
【具体实施方式】
[0047] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对 本发明进行进一步详细说明，应当理解，此处所述的具体实施案例仅仅用于解释本发明，并 不用于限定本发明。
[0048] 为实施对重力异常的小波域优化位变滤波分离，其首先需要进行小波域滤波;小 波域滤波是基于小波变换的滤波方法，简要阐述如下。
[0049] 信号S(X)的连续小波变换定义为：
[0050]
(Ο[0051] 其中，i(X)是母小波，^(X)是f(X)的复共辄，a是尺度因子，b是位置因子。
[0053]公式（1)可用卷积方式描述如下：
[0052] (2、
[0057] 其中，多(6>)是s(x；>的频谱的频谱。[0058] 简记为：
[0054] (3)
[0055]
[0056] (4)
[0059]
(5)
[0060] 为实现滤波功能，将公式(5)两边同乘以向量V，V= (Vi,V2,···，Vk)，其第i个分量Vi 是第i个尺度ai的函数vi = v(ai)(i = l，2,…，k)，Vi称为滤波器的等效系数。可得：
[0065] 这里定义A W)为小波域滤波器转换函数。[0066] 对公式(8)两边做傅立叶反变换，可得到：
[0061] (6)
[0062] (7)
[0063]
[0064] (8)
[0067]
(9)
[0068] 其中：
表示对函数做傅里叶反变换，Y (X)为原信 号S(X)经过小波域滤波后的信号。
[0069]公式(9)表示，信号s'（X)可按如下步骤得到:先对原信号S(X)做小波变换，然后对 小波变换系数W(ai，b)加权求和，权重系数即为小波域滤波器等效系数Vl。
[0070]若能给定小波域滤波器的转换函数由公式(7)可得到小波域滤波器等效 系数V1，即完成小波域滤波器的设计。上述工作是我们进一步实现小波域优化位变滤波的 基础。
[0071 ]接着需要构建小波域位变滤波器，具体构建过程如下：
[0073] (10)
[0072] 假设公式(6)中滤波器等效系数V1不仅与尺度ai有关，而且与空间位置b有关，即为 函数v( ai，b)。则公式(6)可写作：
[0077] 这里定义为小波域位变滤波器转换函数，该转换函数与位置b有关。[0078] 对公式（12)两边做傅立叶反变换，可得到：
[0074] (11)
[0075]
[0076] (12).
[0079] (13)
[0080]
[0081 ]将公式（11)展开，写成如下的方程组形式：
(14)
[0083 ]给定小波函数的频率域表达式f (价乃和小波域位变滤波器转换函数
求解方程组（14)，可得到小波域位变滤波器等效系 数v(ai，b)。由此可以看出，设计小波域位变滤波器的关键是如何确定其位变转换函数 //.,(W) O
[0084] 在实际重磁异常分离中，人为给定小波域位变滤波器的转换函数往往是困难的。 为此，本文根据实测异常数据的区域特征和数据处理的需要，对数据进行分块评价，分析异 常在不同区域的特征差异，以便实现基于差异特征的与位置相关的位变滤波。具体实现中， 我们采用频率域优化滤波法对每块局部数据进行频谱分析，得到小波域位变滤波的转换函 数并计算局部等效系数，进而设计适应全局数据中具有局部特征差异的小波域位变 滤波器。
[0085] 如图1所示，本发明所述重力异常的小波域优化位变滤波分离方法包括如下步骤：
[0086] S1、根据实测异常数据的区域特征和数据处理的需要，对数据做小波变换时频分 析，确定各频率成分在空间的分布规律，并将实测异常数据进行分块，相邻分块间不重叠。
[0087] 具体的，如图4所示，对于图4(a)所示的重力异常，小波变换时频分析结果如图11 所示，从图11可以看出，以x = 3500m为界，左右两边的频率成分空间分布规律具有明显差 异，故将重力异常数据分为x〈 = 3500和x>3500两块。
[0088] S2、对每块实测异常数据进行频率域优化滤波法的特征评价，设计从第j块异常数 据中分离出第1层数据的优化滤波转换函数H^(r)。
[0089] S3、以优化滤波转换函数Hji(r)作为小波域滤波器的转换函数，即^⑷=%#)， 代入方程组（14)，求解该方程组得到第j块异常数据中分离出第1层数据的局部滤波等效系 数向量v(bj) = (v(ai，bj)，v(a2，bj)，· · ·，v(ai，bj)，· · .v(ak，bj))〇
[0090] 具体的，步骤33中的所述局部滤波等效系数向量￥(13」）=(>(31，13」），￥(32，13」），...， v(ai，bj)，...vUi^bj))由以下方程组求得：
[0091]
[0092]
[0093] 其中：f (i/fy)为给定小波函数的频率域表达式;ai、a2--ak为尺度因子;b为位置 因子。
[0094] S4、将每个分块异常数据计算得到的局部滤波等效系数向量V(b〇根据其空间位 置进行组合，得到全体重力异常的小波域位变滤波等效系数向量v(b)。
[0095] 具体的，若重力异常数据分为k块，则等效系数向量= ，v(b2)，...，v (bj)，...v(bk))，其中，V(bj)是第j块数据的局部等效系数向量。所以，简单按照下标的顺 序，可由v(bj)组合得到v(b)。
[0096] S5、对全体重力异常数据做小波变换，得到小波变换系数W(a，b)。
[0097] S6、根据小波域位变尺度滤波等效系数向量V(b)和小波变换系数W(a，b)，从而计 算得到小波域优化位变滤波分离出的第1层的异常数据。
[0098] 具体的，步骤S6中所述第1层的异常数据由以下方程组求得：
[0099]
[0100] 其中:ai为尺度因子;b为位置因子。
[0101] 如图2所示，所述步骤S2包括以下子步骤：
[0102] S21、对第j块异常数据进行傅立叶变换，计算第j块异常数据的对数功率谱Pjr)。
[0103] S22、分段线性拟合第j块异常数据的对数功率谱。
[0104] S23、用格林等效层对数功率谱拟合第j块异常数据实测对数功率谱Pjr)。
[0105] S24、假设目标层与其他深度层场源信息互不相关，则可设计从第j块异常数据中 分离出第1层数据的优化滤波转换函数。具体为
[0106]
(15)
[0107] 其中，Pj(r)是j块异常数据的功率谱密度函数，P^(r)是第j块异常数据第1层异常 数据的功率谱密度函数。
[0108] 具体的，如图3所示，步骤S22还包括以下子步骤：
[0109] S221、分析第j块异常数据的对数功率谱，确定频段数和各频段的频率范围，每一 频段对应一个格林等效层；
[0110] S222、利用直线Pji(r)=Kjir+Bji拟合每一频段异常数据对数功率谱，其中，Kji，Bji 分别是拟合第j块异常数据的第1频段的直线的斜率和截距。
[0111]下面根据附图对本发明所述重力异常的小波域优化位变滤波分离的方法进行具 体阐述：
[0112] 图4为组合模型的重力异常，其中，图4(b)所示的模型长方体1~4顶深分别为30 米、140米、180米、490米，矩形截面分别为65 X 65、138 X 138、385 X 100、1650 X 165,异常密 度分别为0.38/〇113、0.258/〇113、0.38/〇113、0.258/〇113，模型对应的重力异常如图4( &)所示。
[0113] 如图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归算为浅部场源，而小波位变滤波则具有这样的能力。
[0114] 图6所示的是按照本发明所述的方法设计小波域优化位变滤波器，模型重力异常 分离结果。其中，位变低通、带通滤波分离结果为深部区域异常和浅部局部异常，位变高通 滤波分离结果为拟合误差，由优化滤波拟合误差、小波域位变滤波误差和分块频谱估计所 引入误差的叠加而成。从图6所示的小波优化位变滤波分离结果可以看出，仅需采用一次小 波域位变滤波，即可有效分离出四个长方体的异常。其中，分离出的区域异常包含长方体2、 4的异常，局部异常包含长方体1、3的异常，与给定的模型符合性明显改进。
[0115] 图7所示的是某一重力异常剖面实测资料，在地表测点间距0.25km，剖面长度为 380km〇
[0116] 图8所示的是实测重力异常的频率域优化滤波法分离结果。从图8可以看出，区域 异常形态简单，仅反映出该剖面异常具有两边高、中间低的大尺度(波长大于180km)特征； 相应地，局部异常形态过于复杂，包含了波长从6km到180km的所有异常，局部异常分离结果 不够理想。
[0117] 图9所示的是按照本发明所述的方法设计小波域优化位变滤波器，实测剖面重力 异常分离结果。从图9可以看出，小波域优化位变滤波分离结果中，区域异常信息丰富，不仅 反映出该剖面异常所具有的两边高、中间低的大尺度特征，还包含较多的中等尺度异常;相 应的，局部异常则分辨得更加清晰，仅含有小尺度的浅部场源异常。
[0118] 图10所示的是实测剖面重力异常密度反演结果图。对比图8和图10则会发现，频率 域优化滤波分离出的区域异常过于简单，与反演结果不能很好的对应，而局部异常中残留 有较多区域异常，既包含深度小于IOkm的浅部场源异常，又包含深度10~50km的深部场源 异常。对比图9和图10可以看出，小波域优化位变滤波分离出的区域异常与反演结果中深度 大于10 km的场源分布规律一致，局部异常与反演结果中深度小于10 km的场源分布规律一 致，小波域优化位变滤波的分离结果与反演解释结果吻合得相当好。
[0119] 综上所述，本发明所述重力异常的小波域优化位变滤波分离方法，通过在小波域 滤波方法的基础上增加位变滤波功能，并采用优化滤波方法评价分块数据的局部频谱特 性，使所述的方法具有空间变化滤波能力，在不同空间位置表现出不同的滤波器特性，从而 适用于局部数据频谱与全局数据频谱存在较大差异的重力异常分离问题。附图实验表明， 在局部数据频谱与全局数据频谱差异较大时，该方法相对于频率域优化滤波法具有一定优 势且与反演解释成果吻合得更好。
[0120] 以上所述仅为本发明的较佳实施例，并不用于限制本发明，凡在本发明的精神和 原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
【主权项】
1. 一种重力异常的小波域优化位变滤波分离方法，其特征在于，包括如下步骤： 51、 根据实测异常数据的区域特征和数据处理的需要，将实测异常数据进行分块，相邻 分块间不重叠； 52、 对每块实测数据进行频率域优化滤波法的特征评价，设计从第j块异常数据中分离 出第1层数据的优化滤波转换函数田i(r); 53、 W优化滤波转换函数&i(r)作为小波域滤波器的转换函数，求得从第j块异常数据 中分离出第1层数据的局部滤波等效系数向量v(bj) = (v(ai，bj)，v(a2，bj)，...，v(ai， bj),...v(ak,bj))； 54、 将每个分块异常数据计算得到的局部滤波等效系数向量v(bj)根据其空间位置进行 组合，得到全体重力异常的小波域位变滤波等效系数向量v(b); 55、 对全体重力异常数据做小波变换，得到小波变换系数W(a，b); 56、 根据小波域位变尺度滤波等效系数向量v(b)和小波变换系数W(a，b)，从而计算得 到小波域优化位变滤波分离出的第1层的异常数据。2. 根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，所述 步骤S1中，在将实测异常数据进行分块之前，需要对数据做小波变换时频分析，确定各频率 成分在空间的分布规律。3. 根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，所述 步骤S2包括W下子步骤： 521、 对第j块异常数据进行傅立叶变换，计算第j块异常数据的对数功率谱P^r); 522、 分段线性拟合第j块异常数据的对数功率谱； 523、 用格林等效层对数功率谱拟合第j块异常数据实测对数功率谱P^r); 524、 如果目标层与其他深度层场源信息互不相关，则从第j块异常数据中分离出第1层 数据的优化滤波转换函数。4. 根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，所述 从第j块异常数据中分离出第1层数据的优化滤波转换函数如下：其中，P^r)是j块异常数据的功率谱密度函数，Pw(r)是第j块异常数据第1层异常数据 的功率谱密度函数。5. 根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，步骤 S22还包括W下子步骤： 5221、 分析第j块异常数据的对数功率谱，确定频段数和各频段的频率范围，每一频段 对应一个格林等效层； 5222、 利用直线Pji(r)=Kjir+Bji拟合每一频段异常数据对数功率谱，其中，Kji，Bji分别 是拟合第j块异常数据的第1频段的直线的斜率和截距。6. 根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，步骤 S3中的所述局部滤波等效系数向量v(bj) = (v(ai，bj)，v(a2，bj)，. . .，v(ai，bj)，. . .v(ak, bj))由W下方程组求得：其中：(化的为给定小波函数的频率域表达式;ai、a2. . . .ak为尺度因子;b为位置因子。7.根据权利要求1所述的重力异常的小波域优化位变滤波分离方法，其特征在于，步 骤S6中所述第1层的异常数据由W下方程组求得：其中：曰1为尺度因子;b为位置因子。
【文档编号】G01V7/06GK105842745SQ201610132386
【公开日】2016年8月10日
【申请日】2016年3月9日
【发明人】刘彩云, 熊杰
【申请人】长江大学