一种剩余重力异常场分离方法与流程

本发明涉及油气藏生产、蒸汽腔发育和溶洞查找地球物理领域，具体地涉及微重力监测方案中剩余重力异常场分离技术。

背景技术：

微重力监测的前提是地下密度的变化，如油气藏开采、地下水迁移、水合物、矿产资源的开采和蒸汽腔发育监测等物理过程，在不同时期对测点进行时移微重力观测可以观察到重力值的变化。在油气藏资源开发过程中，时移微重力可以监测到地下注气注水动态变化、剩余油分布和低物性夹层分布情况，给油气生产提供指导性建议。重力异常分离在微重力数据处理环节中起着承上启下的重要作用，既能将野外重力观测数据转换为与地下密度体直接相关的参数，还能为后续反演工作建立精准的初始密度模型提供准确先验信息。

早期的趋势分析法、延拓法、圆周法和最小二乘法分离出的区域场分离不彻底，包含大量的残留信息，必须依靠经验和地质信息才能得到合理的解释。前人在异常分离方法方面上做了大量研究，熵滤波法、小波变换法和非线性滤波法等方法进一步提高了异常分离效果，优选延拓法和优化滤波方法结合格林等效层理论和维纳滤波在分离中削弱了长波信号的衰减，被证实具有更好的分离效果。然而，目前主流的异常分离方法依旧存在分离不彻底、平面位置发生偏移和深度上异常混叠等单种或多种问题，对后续的地质解释工作带来了极大的困扰。

技术实现要素：

为了解决上述之一问题，本发明提供了一种剩余重力异常场分离的方法，包括如下步骤：

步骤1接收原始观测重力异常场gg作为输入信号；

步骤2对该输入信号进行各项校正得到布格重力异常场gb；

步骤3对所述布格重力异常场gb进行扩边处理得到第0层区域重力场g⁰r；

步骤4对所述第0层区域重力场g⁰r逐次进行双树复小波变换；

步骤5根据所述双树复小波变换的分解层数j估算阈值函数λj；

步骤6对所述阈值函数用于估计第j层小波系数

步骤7取所述第j层估计小波系数重构低频信号场即区域异常场g^jr；

步骤8将所述第j-1层区域异常场g^j-1r减去所述区域异常场g^jr即剩余异常场g^jl.

优选地，步骤1中所述所接收的原始观测重力异常场gg是由精密的重力测量仪器如cg6、cg5或burris重力仪等在工区严格按照重力测量规范所测量得到的观测重力场。

优选地，步骤2中所述布格重力异常需要消除与时间和空间位置相关因素对重力场的影响，与时间相关的校正项具体包括格值校正、固体潮校正和零漂校正，与空间位置相关的校正项具体包括中间层校正、高度校正、地形校正和正常场校正,校正后的布格重力异常反映的是地下剩余密度体和深部基底起伏在地表引起的重力场。

优选地，步骤3中所述布格重力异常信号的行列点数需要满足条件为2ⁿ；当点数不够时，需要对数据进行扩边处理，直到行列点数达到2的幂次方.

优选地，步骤4中所述布格重力异常信号为一维曲线剖面时，使用一维双树复小波变换方法处理，当布格重力异常为二维平面数据时，按照矩阵排列并将信号输入两树中进行变换。

优选地，步骤5中所述阈值是随分解层数变换的函数，随着分解层数的增加，高频小波系数幅值急剧减小，低频小波系数幅值逐步增大。因此，阈值函数能适应性地调整以应对不同频率小波系数变化，确保分离结果准确性的同时保留更多细节信息。

优选地，步骤6中所述阈值函数采用软阈值法，使得小波系数变化更为平滑，具体表达式为：

其中，为小波系数，为估计小波系数，λj为阈值函数，i＝1,2,3，代表水平、垂直和对角三个方向，j为小波变换分解层数，j＝1,2,3…。

优选地，步骤7中所述分解层数j为从1开始的正整数，在确认总分解层数j之后，从第一层开始逐层进行双树复小波分解、阈值处理和重构依次处理，直至获得第j层极其以上所有层的区域剩余重力场和局部剩余重力场。

优选地，步骤8中所述取所述第j-1层区域剩余重力场g^j-1r减去所述第j层区域剩余重力场g^jr获得第j层局部剩余重力场g^jl,表述为：

优选地，步骤4中的两树高通滤波器和低通滤波器延迟需要满足相差半个采样周期，即g0＝h0(n-0.5)，使得所有采样点信息得到保留，在重构过程中具有完美重构性。

本发明涉及一种剩余重力场分离方法，在该方法中，在接收一个经过各项校正和数值转换后的布格重力异常数据gb，将所述布格重力异常场gb进行扩边处理得到第0层区域重力异常场g⁰r，取所述第0层区域重力异常场g⁰r从第0层开始逐层进行双树复小波变换，根据分解层数j估算阈值函数λ(j)并对第j层小波系数进行阈值处理得到估计小波系数，取所述估计小波系数重构低频信号场即第j层区域异常场g^jr，取所述第j-1层区域异常场g^j-1r减去所述第j层区域异常场g^jr即得到第j层局部异常场g^jl，如前所述依次得到第1～j层区域异常场和局部异常场。

本发明方法的依据是：

1)微重力数据具有以下特性：1.数据精度极高，对数据采集和数据处理具有很高的要求；2.剩余重力异常尺度较小，处理过程中轻微的偏移会在结果中造成巨大的畸变，对方法反映敏感。因此，现有的针对普通重力勘探的异常分离方法不适用于高精度和小尺度的微重力数据处。

2)小波变换类方法具有数学放大镜的特点，能多尺度地提取目标的细微特征，适用于有着高精度和小尺度需求的微重力剩余场分离，其各层分解的多个小波系数也为后期阈值设置提供更多可操作性。

3)小波分解过程中随着分解层数的增加，不同频率小波系数幅值变化趋势不同，高频小波系数幅值在不同尺度都有分量体现，在高层分解中幅值逐渐变小，因此单一的硬阈值方法无法适应分解特性。此外，硬阈值法由于缺乏连贯性，常会出现吉布斯效应和振铃效应。

4)现有的软阈值法提出使用尺度参数，进一步提高了不同频率信号分离能力，但在细节信息保留上仍存在问题，基于此本发明提出了一种新软阈值法。

本发明具有以下优点：

1)本发明中涉及的双树复小波变换计算效率更高，并且对多维信号进行滤波提供了更好的方向选择性。双树复小波变换在两树上分别进行相隔一个采样间隔的滤波器对信号进行滤波，不仅避免了信号能量的损失，具备完美重构性，计算效率也得到了大幅提升。

2)本发明涉及的双树复小波变换在剩余重力异常场提取方面具有更好的平面稳定性。其分离出的剩余重力场不会产生异常位置偏移，并且异常形态更准确清晰。由于双树复小波变换分解出的高频子带在多个方向的特性，其具备更好的方向选择性和更弱的深度混叠性。

3)本发明涉及的阈值约束方面，新软阈值法能在最大程度上对小波系数进行约束，使得小波变换分离出的剩余重力场削弱了随机噪声并突出异常值，提高分离准确性。

4)本发明涉及的剩余重力场分离方法对于微重力场分离具有极其重要的意义，其分离结果具有更好的平面稳定性、方向选择性和更弱的深度混叠性，分离结果有利于识别异常体，为后期圈划异常范围、地质体判识和微重力反演初始密度模型约束等提供关键信息。

附图说明

通过结合附图对本发明的实施例进行描述，可以更好地理解本发明，在附图中：

图1为本发明一个实施例中的剩余重力异常场分离方法的流程图；

图2为本发明设计的通用总重力场组成示意图，剩余重力异常分离目的是将目标体产生的重力异常从总异常场中分离出来，其中需要尽可能地去除近源体产生的重力场；

图3为本发明一个实施例中剩余重力异常场分离示意图，每一层局部剩余重力异常场均在前一层区域剩余重力异常场基础上进行分离；

图4为本发明一个实施例中新软阈值法数据处理示意图，具体显示了阈值处置形式；

图5为本发明一个实施例中所述通用模型的俯视图和透视图；

图6中a-b显示了本发明实施例中描述的方法分离出的b类长方体产生的局部剩余重力场与离散小波变换分离的局部剩余重力场对比图。

图7中a、b为本发明实施例中描述的方法分离出的bc两类长方体产生的局部剩余重力场与离散小波变换分离的局部剩余重力场对比图；

图8中a为本发明实施例中描述的方法分离出的bc两类长方体产生的剩余重力场剖面曲线峰值细节与离散小波变换分离结果峰值细节对比图，b为本发明实施例中描述的方法分离出的bc两类长方体产生的剩余重力场剖面曲线平稳细节与离散小波变换分离结果平稳细节对比图；

图9中a、b为离散小波变换和双树复小波变换分离出的某油田sagd区油层剩余重力异常图。

具体实施方式：

为了使本领域技术人员更好地理解本发明方案，下面将结合附图和说明对本发明技术方案进行清除、完整地描述。

本文所述实施例仅为发明一部分实施例，而非全部实施例。

本领域技术人员在无创造性劳动前提下所获得的其他实施例均属于本发明保护范围。

本文中使用的近似性语言可用于定量表述，表明在不改变基本功能的情况下可允许数量有一定的变动。因此，在使用“大约”、“左右”和“将近”等语言所修正的数值在正负百分之十(10％)的范围内变化，比如，“大约100”表示的可以使90至100之间的任何数值。此外，在“大约第一数值到第二数值”的表述中，“大约”同时修正第一数值和第二数值两个数值。在某些情况下，近似性语言可能与测量仪器的精度有关。

除有定义外，本文中所使用的的技术和科学术语具有与本发明所述领域技术人员普遍理解的相同含义。本文所述“一”或“一个”不表示数量的限定，而是表示存在一个的相关项目。

本发明的实施例涉及一种剩余重力异常场分离方法，其可精准地分离出区域剩余重力场和局部剩余重力场，该方法一般用于重力场数据处理流程中。

参见图2，本发明的实施例中，布格重力异常场gb由目标体及其周围剩余密度体产生的重力异常叠加而成。经发明人研究发现，布格重力异常场gb可视为区域剩余重力场gr和局部剩余重力异常场gl叠加组成，采取适当方法能将剩余场分离出来。

本发明的创新之处在于采用了一种新的软阈值法，小波变换随着分解层数增加，高频小波系数和低频小波系数幅值幅值变化趋势不同，根据分解层数使用尺度函数和指数函数设置阈值能适应小波分解特性，保留了更多的细节信息。

本发明的另一创新之处在于利用双树复小波多尺度和对细微目标特征具有高敏感度的特点，应用于有着高精度和高分辨率需求的微重力剩余场分离，并且双树复小波变换在每层行变换或者列变换过程中都会产生四个频率子带，对多方向、多频段小波系数进行阈值设置具备更高的敏感性和可操作性。

本文涉及双树复小波变换的知识，为了方便理解，以下对所述双树复小波变换方法进行介绍。

本发明中涉及的双树复小波变换采用了二叉树结构的双树离散小波的并行变换，一树生成变换的实部，一树生成变换的虚部。在分解过程中，两树都进行高通和低通滤波，如果滤波器的延迟是一个采样间隔，就能保证一树中采样的值恰好是另一树采样过程中丢失的数值。

所述双树复小波变换将两树的不同层次之间交替使用奇偶滤波器来保持对称性，巧妙地实现了两树对应滤波器间的半采样周期延迟。

所述输入信号为一维信号时，一维双树复小波变换利用两个实小波ψh(t)和ψg(t)形成一个复小波ψ(t)＝ψh(t)+jψg(t)，并且构成希尔伯特变换，滤波器的延迟为半个采样周期g0＝h0(n-0.5)。

所述输入信号为二维信号时，二维双树复小波变换通过两树分别作用于数据的行和列，一般在实际应用中会采用四个二维离散小波来实现二维双树复小波变换。

本发明涉及的双树复小波变换方法良好的平移稳定性，能完美地重构并且具有有限的冗余和高效的可计算性，并且在±15°、±45°和±75°六个方向具有良好的方向选择性。复小波变换域中的小尺度系数能准确地确定奇异点的位置和边缘，因而被广泛地运用到识别图像奇异性和边缘检测等领域，具有高精度和小尺度特性的微重力检测技术适用于双树复小波变换。

本发明涉及的双树复小波变换采用二叉树结构在两个树都进行高通和低通滤波，在两树分别设置低通滤波器h0和g0和高通滤波器h1和g1，当两树的滤波器延迟刚好是一个延迟时，就能保证一树采样损失的信号能量在另一树采样过程补充回来。

本发明所述一种剩余重力异常场分离方法不仅适用于二维重力场信号，也适用于一维重力场信号。

本发明所述一种剩余重力异常场分离方法不仅适用于地球物理领域中重力勘探方法，也适用于磁法勘探方法、图像去噪和图像特征提取等领域。

图1为一种剩余重力异常场分离方法的流程图。图3为图1所示的流程图中数据处理具体方式。图4为图3所示中处理步骤12的具体示意图。

所述所接收的原始观测重力异常场gg是由精密的重力测量仪器如cg6、cg5或burris重力仪等在工区严格按照重力测量规范所测量得到的观测重力场。

所述观测重力场受到与时间和空间位置相关的因素影响，因此在做进一步分析和研究前，需要计算各校正项消除干扰。所述观测重力场与时间相关校正项包括但不限于格值校正、固体潮校正和零漂校正，与空间位置相关校正项包括但不限于中间层校正、高度校正、地形校正和正常场校正。

所述经过各项校正后的布格重力异常场gb反映地下剩余密度体和深部基底起伏在地表引起的重力变化。为消除边界影响以及小波变换要求待处理数据长度需满足2ⁿ特性，所述布格重力异常场信号需添加数据长度进行扩边处理。

所述扩边处理后布格重力异常场gb可视为第0层区域重力场g⁰r，此时所述第0层区域重力场包括但不限于行和列的数据长度应满足2ⁿ要求。

所述区域异常场信号输入后，一树上的低通滤波器和高通滤波器h0和h1对信号作行卷积和系数为2的列下抽样，再使用另一树的低通滤波器和高通滤波器g0和g1对信号作列卷积和系数为2的行抽样。同样地，对另一树的输入信号使用低通滤波器和高通滤波器g0和g1对信号作行卷积和系数为2的列下抽样。再使用另一树上的低通滤波器和高通滤波器h0和h1作列卷积和系数为2的行抽样。每一次行变换或列变换都可以分别得到1个低频子带和3个高频子带，每对子带的和或者差构成两个低频系数和±15°、±45°、±75°系数的6个高频子带；

本发明涉及的双树复小波方法将输入的二维重力信号进行分解，各层分解产生多个低频小波系数和高频小波系数。

参见图4，所述输入小波系数要经过阈值处理，阈值λj的定义为：

其中，j为小波变换分解层数，j＝1,2,3…。i代表水平、垂直和对角三个方向，σj，i为第j层小波系数标准方差，其定义为：

本发明引入尺度函数α，使用了指数形式，将其定义为：

其中，j为总分解层数，lj为第j层信号数据大小，β的定义为：

β＝(2j+3j)/4j

本发明涉及的阈值函数中调节因子β与分解层数j和总分解层数j有关，随着分解层数j的增加而增大。

本发明涉及的阈值函数中尺度函数α与调节因子β、第j层信号数据大小、分解层数j和总分解层数j有关，采用了指数函数形式和以10为底的对数函数，自适应性地调节各层尺度函数大小。

通过小波系数与阈值λj对比，若则估计小波系数等于若则估计小波系数等于0。表达式为：

所述第j层估计小波系数用于重构第j层低频信号场，即第j层区域重力异常场g^jr。

需要说明的是，分解层数j为从1开始的正整数。在确认总分解层数j之后，从第一层开始逐层进行双树复小波分解、阈值处理和重构依次处理，直至获得第j层极其以上所有层的区域剩余重力场和局部剩余重力场。

取所述第j-1层区域剩余重力场g^j-1r减去所述第j层区域剩余重力场g^jr获得第j层局部剩余重力场g^jl。表述为：

本发明实施例中，布格重力异常场由各层区域剩余重力场和局部剩余重力场组成，后一层局部剩余重力异常场可由前一层区域剩余重力异常场减去后一层区域剩余重力异常场获得，如第三层局部剩余重力异常场g³l由第二层区域剩余剩余重力异常场g²r减去后一层区域剩余重力异常场g³r获得。

本发明实施例中，当分解层数为0时，即布格重力异常场尚未获得双树复小波变换分解处理，因此可视为第0层区域剩余重力异常场。

为了突出本发明方法的优势特点，以下结合具体场景实施例进行详细说明

参见图5，本发明一个实施例中所述通用模型的俯视图和透视图，其由不同大小、埋深和剩余密度的11个长方体组成，长方体分为a、b和c三类不同埋深的地质异常体，分别代表了深、中和浅层地质体异常。

本发明涉及的通用模型具备典型代表性，多方面考察分离方法特性，适用于对比分离效果。

图6中，a-b显示了本发明实施例中描述的方法分离出的b类长方体产生的局部剩余重力场与离散小波变换分离的局部剩余重力场对比图，离散小波方法分离结果虽然依稀能根据线框分辨出异常位置，但包含太多杂乱的干扰信息，如b1、和b5的两侧，b4的上下侧以及b3的左下角，并所有b类异常体位置都存在偏移现象，甚至还残留着较多的c类异常体的信号能量。本发明实施例所描述方法分离结果模型边界清晰干净，不受边界效应影响，没有虚假异常；异常清晰，异常边界分辨率高，与真实位置吻合良好；无其他深度异常体残留信息，分离更彻底。

图7中，a-b为本发明实施例中描述的方法分离出的bc两类长方体产生的局部剩余重力场与离散小波变换分离的局部剩余重力场对比图，离散小波变换分离结果中异常位置与异常体真实位置依旧存在明显偏移，异常边界模糊；本发明实施例所描述方法分离结果异常位置准确，基本无偏移。但归一化后的局部异常值较小，导致异常不明显。对于重要检验指标——b2和c2的异常分离，离散小波变换方法无法准确地将这两个异常体准确地分离开，不仅包含杂乱的虚假异常，而且b2和c2两个异常位置发生了偏移，而本发明实施例所描述方法能达到很好的分离效果。

图8中，a为本发明实施例中描述的方法分离出的bc两类长方体产生的剩余重力场剖面曲线峰值细节与离散小波变换分离结果峰值细节对比图。b1异常体处本发明实施例所描述方法曲线峰值顶点与参考线峰值顶点投影位置基本一致，离散小波变换曲线出现向右偏移现象；b2异常体处本发明实施例所描述方法曲线峰值顶点与参考线峰值顶点投影位置一致，离散小波变换曲线形态两侧不对称，且向右偏移；b3异常体处本发明实施例所描述方法曲线峰值顶点与参考线峰值顶点投影位置基本一致。

图8中，b为本发明实施例中描述的方法分离出的bc两类长方体产生的剩余重力场剖面曲线平稳细节与离散小波变换分离结果平稳细节对比图。在没有异常体的曲线段，参考线应当趋近于零。在b1～b2曲线段，离散小波变换曲线出现“振荡”，本发明实施例所描述方法曲线轻微波动后很快逼近参考线。在b2～b3曲线段，受异常体c2的影响，离散小波曲线残留了浅层异常体信号能量，形成槽状异常。本发明实施例所描述方法曲线更趋近于参考线。本发明实施例所描述方法具有更好的平面稳定性。

图9中a、b为离散小波变换和本发明实施例所描述方法分离出的某油田sagd区油层剩余重力异常图，其中黑色边框为油藏边界，红色线条为水平生产井，彩色点为观察井。

参见图9中b，在本发明实施例所描述方法分离结果中，剩余重力异常值较高的区域含剩余油较多，蒸汽腔尚未形成，对应的井温也较低，如well-2、well-7、well-11和well-12观察井；剩余重力异常值低值区的蒸汽腔发育良好，含剩余油较少，井温较高，如well-3、well-4和well-6观察井。其余观察井的井温也与异常值大小有显著的负相关关系。

参见图9中a为离散小波变换分离结果中部分区域存在异常形态变化、异常区域位置偏移和异常相对值变化等现象。如低温井well-2、well-8和高温井well-6与剩余异常对应，但高温井well-3、well-4和低温井well-11、well-12等井温数据与剩余异常无法对应。离散小波变换方法本身特性造成的分离结果异常形态形变、异常位置偏移和方向细节信息丢失等缺点无法在实际资料处理中获得准确可靠的结果。