CN109407161B - 用于提取地球物理磁异常场边界的磁场刻痕分析方法 - Google Patents

Abstract

本发明公开了一种用于提取地球物理磁异常场边界的磁场刻痕分析方法，利用磁场与重力场之间的解析关系即泊松方程，推导出磁场表面二阶谱矩的三个分量，获得磁场表面刻痕分析的结果参数。该方法包括四个步骤：一是磁场分量转换；二是总磁场的模量异常ΔT总水平分量模表面局部面元二阶谱矩；三是总磁场的模量异常ΔT总水平分量模表面局部面元统计不变量；四是磁场表面提取的边界系数四个步骤。本发明的有益效果是：利用磁场刻痕分析方法对磁场进行边界提取，保留了重力场刻痕分析对弱边界具有更强的增强能力、分辨能力和噪音压制能力的同时，消除了由于伴生异常引起的虚假边界，可为地质构造的解释提供更有效的参考信息。

Description

用于提取地球物理磁异常场边界的磁场刻痕分析方法

技术领域

本发明涉及地球物理信息处理技术领域，具体地说，涉及一种用于提取地球物理磁异常场边界的磁场刻痕分析方法。它可以从地球物理磁异常表面提取边界信息，消除由于磁场伴生异常引起的虚假边界，从而可以为断裂划分、磁性体边界圈定及地质解释等提供更为准确、客观的参考信息。

背景技术

因为从重、磁异常上提取的边界可用于地质构造解释，如推断断裂构造、指示大地构造单元圈定地质体边界等，所以国内外学者对重、磁异常的边界提取进行了大量的研究并发展了多种方法技术。重、磁异常边界识别大致可划分为三个阶段：第一，早期利用简单导数方法提取边界阶段，如垂直导数法(Evjen,1936)、水平梯度法(cordell,1979)、二维解析信号(Nabighian,1972)及三维解析信号(Roest，1997)等；第二，对边界进行均衡提取阶段，如倾斜角方法(Miller和Sigh,1994)、倾斜角的水平总梯度(Verduzco，2004)以及Thetamap方法(Verduzco，2004)；第三，精细化提取边界阶段，如Cooper和Cowan(2006)利用水平与垂直导数的均方差(NSTD)来增强对弱异常边界的识别，马国庆(2012,2013)提出增强型均衡滤波器和增强型局部相位滤波器，孙艳云(2016)提出重力场刻痕分析方法可更精细地识别重力场表面的边界。

上述重、磁异常上边界提取方法各有优缺点。总体上，这些不足的地方可以归纳为以下几点：首先，提取的边界较宽；其次，提取的边界主要反映具有较大幅值异常的构造边界，而较小幅值异常的边界识别不清晰；最后，这些方法对噪音较敏感。针对上述不足，本人于2016年将信息学中表面识别技术引入重力场领域，提出了重力场刻痕分析方法。该方法具有更强的均衡能力、分辨能力及对噪音的压制能力。

由于与重力场相比，磁场与地下异常体之间存在更复杂的关系，所以如果将重力场刻痕分析方法直接用于磁场进行边界提取，就会出现许多与伴生异常相关的虚假边界，这些虚假边界将为后续的地质解释提供错误的信息。目前，大多数边界识别方法在对磁场进行边界提取时都存在这种问题，并且这种虚假边界随着边界识别方法对弱边界增强的能力增大而变多变强，该问题在业界尚未得到解决。

发明内容

本发明正是为了解决上述技术问题而设计的一种用于提取地球物理磁异常场边界的方法(磁场刻痕分析方法)，它利用磁场与重力场之间的解析关系(泊松方程)推导出磁场表面二阶谱矩的三个分量，获得磁场表面刻痕分析的结果参数，对弱边界具有更强的增强能力和分辨能力，同时可消除由于磁场伴生异常引起的虚假边界，可为后续的地质解释提供更有效的参考信息。

本发明解决其技术问题所采用的技术方案是：

所述用于提取地球物理磁异常场边界的磁场刻痕分析方法，其利用磁场与重力场之间的解析关系即泊松方程，推导出磁场表面二阶谱矩的三个分量，获得磁场表面刻痕分析的结果参数，对弱边界具有更强的增强能力和分辨能力，同时可消除由于磁场伴生异常引起的虚假边界，可为后续的地质解释提供更有效的参考图件，该方法具体步骤如下：

步骤一：磁场分量转换

磁场各分量之间的转换通常可以通过两种不同的计算方式来实现，即空间域和频率域两种方式。由于频率域运算比较简单，在频率域中进行磁场各分量转换是目前比较常用的方法。

由磁场与磁位的关系可以得到以下磁场各分量之间的关系式：

其中X_a、Y_a、Z_a分别为磁场异常水平和垂直分量；t₀为地磁场方向的单位矢量。

若设S_x(u,v,z)、S_y(u,v,z)、S_z(u,v,z)以及S_T(u,v,z)分别为H_ax(x,y,z)、H_ay(x,y,z)、Z_a(x,y,z)及其ΔT(x,y,z)的频谱。

利用频谱微分定理可得到上述磁场各分量导数在频率域内相应的换算关系式：

其中q_t0＝2π[i(L₀u+M₀u)+N₀(u²+v²)^1/2]，而L₀、M₀、N₀为地磁场单位矢量t₀的方向余弦；

实际磁测一般更易于精确测定总磁场的模量ΔT，由上述公式便可得到H_ax和H_ay分量。

步骤二：计算总磁场的模量异常ΔT的总水平分量模B，并在总水平分量模B上计算二阶谱矩和统计不变量；

B为总磁场的模量异常ΔT总水平分量模，计算如下：

在总水平分量模B表面，选取滑动窗口w₁用于计算B表面局部面元内的二阶谱矩和统计不变量，二阶谱矩的离散表达式为：

其中i＝1,2,...N_α,j＝1,2,...M_α,并且B(x_j,y_i)是滑动窗口w₁截取的B表面局部面元；N_α、M_α分别为滑动窗口w₁内x和y方向的点数；

B表面二阶谱矩的三个面元都有其各自的物理含义；m₂₀、m₀₂分别表示B表面x和y方向上的斜率的方差，m₁₁为B表面在x和y方向上斜率的协方差。

由于实验表明频率域计算表面谱矩的方式具有局限性，因而本发明在计算二阶谱矩时采用空间域的计算方式。

步骤三：总磁场的模量异常ΔT总水平分量模表面局部面元统计不变量

由ΔT总水平分量模表面局部面元二阶谱矩的三个分量，可计算获得ΔT总水平分量模表面局部面元的统计不变量：

值得注意的是表面二阶谱矩的三个分量随着坐标系的转动而改变，但是统计不变量是不随着坐标系的选取而变化的，因而称为表面的统计不变量。表面的统计不变量M₂是表面斜率的方差，△₂是表面方向效应的一种度量。

步骤四：磁场表面提取的边界系数

磁场的边界系数Λ可以对较小的磁异常的边界进行刻画，计算如下：

其中△为拉普拉斯算子，-sign(ΔB)是为了区分是由“谷”还是“脊”增强获得的边界。这是由“脊”增强而得的边界才是目标边界。磁场表面的边界系数Λ与

成正比，体现了表面形态的脊形化特征；与统计不变量M₂成反比，可增加对磁场表面弱边界的识别能力；

所述用于提取地球物理磁异常场边界的磁场刻痕分析方法，其滑动窗口w₁根据实际的数据网格密度确定，在网格密度[250点×250点]/[1°×1°]时，滑动窗口半径3-5个点。

值得注意的是，本公式是由作者在发明重力场刻痕分析方法时构建的结果参数，在磁场刻痕分析中还仍继续该结果参数的表达形式，但是公式中各变量的计算方式已经不同，具体见上述各步骤。

本发明提供一种用于提取地球物理磁异常场边界的方法(磁场刻痕分析方法)，下面结合实际例子说明其优点：

(1)消除由于伴生异常造成的虚假边界，更有效地从磁场表面提取边界，凸显出有用的构造信息

由于伴生异常而导致的虚假边界现象普遍存在，例如直接利用重力场刻痕分析方法从航磁异常ΔT及其化极异常提取的边界结果(图11和图12)、总水平梯度(图14)、thetamap(图15)和NSTD(图16)方法提取的结果中都存在虚假边界，例如椭圆A圈定位置。边界增强方法的增强能力越强，这种虚假边界越多，越影响后续的构造解释。由图11、图14-图16可知，重力刻痕分析方法的增强效果，总水平梯度方法增强的能力最弱。另外，与图11相比，在化极数据上提取的结果(图12)上，虚假边界受到少许的压制(如图11中箭头a指示位置)。这是由于化极可以压制伴生异常，但是不能消除伴生异常。

与上述方法的结果相比，本发明的方法(图13)极大地消除了由于伴生异常导致的虚假边界(例如图13中椭圆A圈定的区域)；有效的提取实际边界的同时更便于有效边界的构造解释，例如图13中矩形B所圈定的边界为连云港断裂，该断裂在图11中很明显，但是由于伴生异常导致的虚假边界的影响，该断裂在其他方法提取的结果很难清晰识别。

(2)具有更强的弱边界增强能力

本发明的方法源于重力场刻痕分析方法，但是克服了它只能用于重力场的弱点，从而可以直接从磁场上提取边界。因而该方法也具有很强的弱边界增强能力。例如，总水平梯度(图14)和Thetamap方法(图15)在研究区的东南海水覆盖区域，提取的边界信息都很微弱，不利于后期的构造解释。而本文发明的方法依然可以很清楚的识别出这些弱边界信息。

(3)提取的边界较窄并且连续，更便于地质解释

与总水平梯度(图14)和Thetamap(图15)方法提取的边界对比，本发明提取的边界比较细窄，从而增强了方法的水平空间分辨能力，不会让较弱的边界被较强的边界覆盖。Thetamap方法(图15)对弱边界的增强能力很强，但是该方法提取的边界能量不均一，表现出断断续续。与它对比可以发现，本发明提出的方法在提取边界信息时，不但增强了弱边界，而且提取的边界更加的连续。由于本发明的方法消除了伴生异常，而其他方法提取的边界均受到虚假异常的干扰，因而无法具体示例对比说明。

创新点：本发明针对重力场刻痕分析方法的局限性，提出磁场表面刻痕分析方法，利用泊松公式直接推导出磁场表面二阶谱矩的三个分量，获得磁场表面边界系数的直接表达形式，可消除由于伴生异常造成的虚假边界，从磁场表面提取更有效的边界。

本发明的有益效果是：利用磁场刻痕分析方法对磁场进行边界提取，保留了重力场刻痕分析对弱边界具有更强的增强能力、分辨能力和噪音压制能力的同时，消除了由于伴生异常引起的虚假边界，可为地质构造的解释提供更有效的参考图件。

附图说明

图1.为本发明流程框图。

图2.为日照-连云港地区总磁场的模量异常ΔT图。

图3.为由ΔT转换的磁场分量Hax。

图4.为由ΔT转换的磁场分量Hay。

图5.为ΔT总水平分量模。

图6.为二阶谱矩统计不变量M₂。

图7.为二阶谱矩统计不变量Δ₂。

图8.为二阶谱矩分量m₂₀。

图9.为二阶谱矩分量m₀₂。

图10.为二阶谱矩分量m₁₁。

图11.为ΔT上由重力场刻痕分析方法提取结果。

图12.为ΔT化极后由重力场刻痕分析方法提取结果。

图13.为ΔT上由磁场刻痕分析方法提取结果。

图14.为ΔT化极后由总水平梯度导数方法提取结果。

图15.为ΔT上由Thetamap方法提取结果。

图16.为ΔT化极后由NSTD方法提取结果。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

如图1所示，本发明一种用于提取地球物理磁异常场边界的磁场刻痕分析方法，利用磁场与重力场之间的解析关系即泊松方程，推导出磁场表面二阶谱矩的三个分量，获得磁场表面刻痕分析的结果参数，对弱边界具有更强的增强能力和分辨能力，同时可消除由于磁场伴生异常引起的虚假边界，可为后续的地质解释提供更有效的参考图件，该方法具体步骤如下：

步骤一：磁场分量转换

磁场各分量之间的转换通过频率域方式实现；

由磁场与磁位的关系可以得到以下磁场各分量之间的关系式：

其中X_a、Y_a、Z_a分别为磁场异常水平和垂直分量；t₀为地磁场方向的单位矢量；

设S_x(u,v,z)、S_y(u,v,z)、S_z(u,v,z)以及S_T(u,v,z)分别为H_ax(x,y,z)、H_ay(x,y,z)、Z_a(x,y,z)及其ΔT(x,y,z)的频谱；

利用频谱微分定理可得到上述场各分量导数在频率域内相应的换算关系式：

其中q_t0＝2π[i(L₀u+M₀u)+N₀(u²+v²)^1/2]，而L₀、M₀、N₀为地磁场单位矢量t₀的方向余弦；

用磁场测量仪器测定总磁场的模量ΔT，由上述公式便可得到H_ax和H_ay分量；

步骤二：总磁场的模量异常ΔT总水平分量模表面局部面元二阶谱矩

设B为总磁场的模量异常ΔT总水平分量模，计算如下：

在总水平分量模B表面，选取滑动窗口w₁用于计算B表面局部面元内的二阶谱矩和统计不变量，二阶谱矩的离散表达式为：

其中i＝1,2,...N_α,j＝1,2,...M_α,并且B(x_j,y_i)是滑动窗口w₁截取的B表面局部面元；

B表面二阶谱矩的三个元都有其各自的物理含义；m₂₀、m₀₂分别表示表面x和y方向上的斜率的方差，m₁₁为表面在x和y方向上斜率的协方差；

步骤三：总磁场的模量异常ΔT总水平分量模表面局部面元统计不变量

由ΔT总水平分量模表面局部面元二阶谱矩的三个分量，可计算获得ΔT总水平分量模表面局部面元的统计不变量：

表面的统计不变量M₂是表面斜率的方差，△₂是表面方向效应的一种度量；

步骤四：磁场表面提取的边界系数

磁场的边界系数Λ可以对较小的磁异常的边界进行刻画，计算如下：

其中△为拉普拉斯算子，-sign(ΔB)是为了区分是由“谷”还是“脊”增强获得的边界；由“脊”增强而得的边界才是目标边界；表面的边界系数Λ与

成正比，体现了表面形态的脊形化特征；与统计不变量M₂成反比，可增加对磁场表面弱边界的识别能力。

所述用于提取地球物理磁异常场边界的磁场刻痕分析方法，其滑动窗口w₁根据实际的数据网格密度确定，在网格密度[250点×250点]/[1°×1°]时，滑动窗口半径3-5个点。

实施例：

以日照-连云港地区总磁场的模量异常△T为例，说明本发明在实际工程应用中的具体操作步骤。本发明流程框图如图1所示，具体步骤如下：

步骤一：总磁场的模量异常△T数据的输入

选取日照-连云港地区的总磁场的模量异常△T数据作为输入数据，数据为“.grd”格式，网格尺寸为0.004°×0.004°，网格大小为366×496个网格点。日照-连云港的总磁场的模量异常△T示如图2。

步骤二：磁场分量转换与总水平分量模计算

将日照-连云港地区的总磁场模量异常△T数据进行分量转换。转换时选取的国际地磁参考场IGRF参数为：磁倾角I＝52.27°、磁偏角D＝-6.03°。根据公式3可得到磁场分量Hax和Hay，分别见图3和4。

根据获得的磁场水平分量Hax和Hay，利用公式4容易计算出总磁场的模量异常ΔT总水平分量模，如图5。

步骤三：表面二阶谱矩计算

步骤二获得的ΔT总水平分量模数据是本步骤计算的原始数据。由本发明的流程框图可知，在计算二阶谱矩前需使用滑动窗口来截取局部面元，利用公式5计算出该局部面元的二阶谱矩的三个分量，放置于该面元的中心点位置。逐点滑动该窗口并重复上述计算直至获得总水平分量模表面所有点的二阶谱矩三个量。在本示例中，选用半径为3个采样点的圆形滑动窗口，计算获得二阶谱矩三个分量m₂₀、m₀₂、m₁₁分别示如图8-10.m₂₀、m₀₂分别表示总水平分量模表面x和y方向上的斜率的方差，m₁₁为表面在x和y方向上斜率的协方差。

步骤四：统计不变量计算

将步骤三获得的总水平分量模表面二阶谱矩的三个分量输入到公式6中，便可得到表面的统计不变量M₂(图6)和Δ₂(图7)。统计不变量M₂实际上反映的是总水平分量模表面较大异常的边界；而统计不变量Δ₂反映的是航磁总水平分量模表面的各向异性。

步骤五：结果生成

利用公式7可计算获得反映磁异常边界的最终的结果参数，即边界系数。值得注意的是对总水平分量模表面求取散度是为了区别最终计算出的边界来自于总水平分量模表面的“脊”部还是“谷”部，来自“脊”部的边界才是我们真正需要的边界。由日照-连云港地区总磁场的模量异常Δ₂表面提取的边界示如图13。

本发明不局限于上述最佳实施方式，任何人在本发明的启示下得出的其他任何与本发明相同或相近似的产品，均落在本发明的保护范围之内。

Claims (1)

1.一种用于提取地球物理磁异常场边界的磁场刻痕分析方法，其特征在于：利用磁场与重力场之间的解析关系即泊松方程，推导出磁场表面二阶谱矩的三个分量，获得磁场表面刻痕分析的结果参数，对弱边界具有更强的增强能力和分辨能力，同时可消除由于磁场伴生异常引起的虚假边界，可为后续的地质解释提供更有效的参考图件，该方法具体步骤如下：

步骤一：磁场分量转换

磁场各分量之间的转换通过频率域方式实现；

由磁场与磁位的关系可以得到以下磁场各分量之间的关系式：

其中X_a、Y_a、Z_a分别为磁场异常水平和垂直分量；t₀为地磁场方向的单位矢量；

设S_x(u,v,z)、S_y(u,v,z)、S_z(u,v,z)以及S_T(u,v,z)分别为H_ax(x,y,z)、H_ay(x,y,z)、Z_a(x,y,z)及其ΔT(x,y,z)的频谱；

利用频谱微分定理可得到上述磁场各分量导数在频率域内相应的换算关系式：

其中q_t0＝2π[i(L₀u+M₀u)+N₀(u²+v²)^1/2]，而L₀、M₀、N₀为地磁场方向单位矢量t₀的方向余弦；

用磁场测量仪器测定总磁场的模量异常ΔT，由上述公式便可得到H_ax和H_ay分量；

步骤二：计算总磁场的模量异常ΔT的总水平分量模B，并在总水平分量模B上计算二阶谱矩和统计不变量；

B为总磁场的模量异常ΔT总水平分量模，计算如下：

在总水平分量模B表面，选取滑动窗口w₁用于计算B表面局部面元内的二阶谱矩和统计不变量，二阶谱矩的离散表达式为：

其中i＝1,2,...N_α,j＝1,2,...M_α,并且B(x_j,y_i)是滑动窗口w₁截取的B表面局部面元；N_α、M_α分别为滑动窗口w₁内x和y方向的点数；

B表面二阶谱矩的三个面元都有其各自的物理含义；m₂₀、m₀₂分别表示B表面x和y方向上的斜率的方差，m₁₁为B表面在x和y方向上斜率的协方差；

步骤三：总磁场的模量异常ΔT总水平分量模表面局部面元统计不变量

由ΔT总水平分量模表面局部面元二阶谱矩的三个分量，可计算获得ΔT总水平分量模表面局部面元的统计不变量：

表面的统计不变量M₂是表面斜率的方差，△₂是表面方向效应的一种度量；

步骤四：磁场表面提取的边界系数

磁场的边界系数Λ可以对较小的磁异常的边界进行刻画，计算如下：

其中△为拉普拉斯算子，-sign(△B)是为了区分是由“谷”还是“脊”增强获得的边界；由“脊”增强而得的边界才是目标边界；磁场表面的边界系数Λ与

成正比，体现了表面形态的脊形化特征；与统计不变量M₂成反比，可增加对磁场表面弱边界的识别能力；

滑动窗口w₁根据实际的数据网格密度确定，在网格密度[250点×250点]/[1°×1°]时，滑动窗口半径为3-5个点。

Images