CN106405664A - 一种磁异常化极方法 - Google Patents

Abstract

本发明公开了一种磁异常化极方法。包括以下步骤：对采集到的磁场数据进行校正，得到磁异常数据；将磁异常数据形成规则的磁异常数据t(x,y)，记为矩阵T(1:m,1:n)；获得研究区的地磁场方向，并整理成与矩阵T(1:m,1:n)等大的矩阵I(1:m,1:n)；对矩阵T(1:m,1:n)进行镶边处理，转变成新矩阵T’(1‑m’:m+m’,1‑n’:n+n’)；设计频率域化极因子，并做反傅里叶变换，得到化极结果；计算化极结果的垂直梯度和总梯度，并以最大互相关系数确定剩磁方向；根据地磁场方向和剩磁方向，并对矩阵T’和I在频率域内进行变倾角化极处理，计算研究区磁异常的变倾角化极结果。通过上述方式，本发明能够提高磁法勘探的准确度。

Description

一种磁异常化极方法

技术领域

本发明涉及地球物理勘测技术领域，具体为磁法勘探领域，特别是涉及一种磁异常化极方法。

背景技术

在磁法勘探中，磁测数据的处理和解释(如化极、导数换算、分量换算、正演和反演等)通常需要已知磁场方向信息。磁场是一个矢量场，而利用磁力仪测得的磁场强度是标量，没有磁场方向。磁场方向是感应磁化方向和剩余磁化方向的矢量和。利用国际地磁参考场获得感应磁场方向，即地磁场方向；而剩余磁化方向无法直接获取。若磁源地质体剩磁不明显，那么磁场方向与感应磁场方向基本一致；否则，较强的剩磁会使磁场方向与感磁方向出现很大偏差，使得地下磁性体产生斜磁化现象。此时，若直接利用地磁场方向作为磁场方向进行磁异常化极，以及在此基础上进行的一系列磁异常数据处理与转换、正演和反演，其结果都无法反映真实的地质情况。

目前，除采集岩石标本并测定剩磁参数外，国内外学者提出了一些估计剩磁方向的方法。唐俊德(1986)提出利用磁异常三分量确定磁化方向的圆曲线法；Roest和Pilkington(1993)提出基于磁场总梯度模和磁源重力异常水平梯度互相关的估计方法；Medeiros和Silva(1995)提出利用等效源磁矩反演估计磁化方向；甘西(2001)利用磁场的垂直分量确定磁性体磁化方向；Phillips(2005)给出估计磁化方向的Helbig(1963)积分法的直接算法和间接算法；Dannemiller和Li(2006)提出基于化极磁异常垂直梯度与总梯度模互相关的估计方法；Gerovska等(2009)提出基于化极磁异常与总模量磁异常互相关的估计方法。

但是，磁异常化极的准确性还与研究区的纬度变化有关。若研究区的纬度范围较大，地磁倾角变化也会较大，斜磁化的程度会随纬度不同而变化。若采用单一的磁倾角作化极处理显然是不合理的，其结果将无法反映真实的地质情况，降低磁法勘探的准确度。

发明内容

本发明主要解决的技术问题是提供一种磁异常化极方法，能够利用地磁场方向和剩磁方向对实测磁异常进行变倾角化极，提高磁法勘探的准确度。

为解决上述技术问题，本发明采用的一个技术方案是：提供一种磁异常化极方法，包括以下步骤：采集研究区的磁场数据，对采集到的磁场数据进行校正，得到磁异常数据；将磁异常数据网格化，形成规则的磁异常数据t(x,y)，其中x,y分别表示磁异常数据的平面网格点坐标，记为矩阵T(1:m,1:n)，大小为m×n，m表示行数，n表示点数；获得研究区的地磁场方向，并将地磁场方向整理成与矩阵T(1:m,1:n)等大的矩阵，记为I(1:m,1:n)；对矩阵T(1:m,1:n)进行镶边处理，将矩阵T(1:m,1:n)转变成行数和列数均为2的最小整数幂的新矩阵T’(1-m’:m+m’,1-n’:n+n’)；设计频率域化极因子，并对频率域化极因子做反傅里叶变换，得到化极结果；计算化极结果的垂直梯度和总梯度，并计算其互相关系数，以最大互相关系数确定剩磁方向；根据地磁场方向和剩磁方向，并根据频率域化极因子对矩阵T’和I在频率域内进行变倾角化极处理，计算研究区磁异常的变倾角化极结果。

其中，将磁异常数据网格化的步骤包括：利用kriging插值方法将磁异常数据网格化。

其中，将地磁场方向整理成与矩阵T(1:m,1:n)等大的矩阵，记为I(1:m,1:n)的步骤包括：

获取地磁倾角范围为a-b；

设置矩阵

利用矩阵A的伪逆矩阵计算磁场方向；

P＝pinv(A)×V，pinv表示伪逆矩阵P为计算磁场方向时的参数矩阵，其使用方法见下式；

设置网格坐标阵[x,y]＝meshgrid(1:n,1:m)；

计算矩阵I＝P(1)+P(2)×y+P(3)×x。

其中，对矩阵T(1:m,1:n)进行镶边处理的步骤包括：

采用以下公式进行镶边：

其中，T位于T’的中心位置，即T(i,j)＝T’(i,j)，(1≤i≤m,1≤j≤n)。

其中，频率域化极因子为：

式中，ω_x和ω_y是x、y方向的圆波数， 是地磁方向的单位向量，L、M、N是地磁场的方向余弦：cosIcosD,cosIsinD,sinI，I和D分别表示地磁场倾角和地磁场偏角；是剩磁方向的单位向量，l,m,n分别是剩磁方向的方向余弦：cosicosd,cosisind,sini，i和d分别表示剩磁倾角和剩磁偏角。

其中，对频率域化极因子做反傅里叶变换的步骤包括：

采用以下的公式做反傅里叶变换：

其中，RTP表示化极结果。

其中，根据频率域化极因子对矩阵T’、I在频率域内进行变倾角化极处理的步骤之前包括：

从矩阵I中提取平均地磁倾角为平均地磁方向；

计算出研究区中的实际地磁方向与平均地磁方向的差

其中，利用泰勒公式进行变倾角化极，采用的公式是：

式中，RTP_diff是变倾角化极结果，RTP_mean是使用研究区平均地磁倾角和地磁偏角得到的化极结果。

本发明的有益效果是：区别于现有技术的情况，本发明提供一种磁异常化极方法，包括以下步骤：采集研究区的磁场数据，对采集到的磁场数据进行校正，得到磁异常数据；将磁异常数据网格化，形成规则的磁异常数据t(x,y)，其中x,y分别表示磁异常数据的平面网格点坐标，记为矩阵T(1:m,1:n)，大小为m×n，m表示行数，n表示点数；获得研究区的地磁场方向，并将地磁场方向整理成与矩阵T(1:m,1:n)等大的矩阵，记为I(1:m,1:n)；对矩阵T(1:m,1:n)进行镶边处理，将矩阵T(1:m,1:n)转变成行数和列数均为2的最小整数幂的新矩阵T’(1-m’:m+m’,1-n’:n+n’)；设计频率域化极因子，并对频率域化极因子做反傅里叶变换，得到化极结果；计算化极结果的垂直梯度和总梯度，并计算其互相关系数，以最大互相关系数确定剩磁方向；根据地磁场方向和剩磁方向，并根据频率域化极因子对矩阵T’和I在频率域内进行变倾角化极处理，计算研究区磁异常的变倾角化极结果。因此，本发明能够利用地磁场方向和剩磁方向对实测磁异常进行变倾角化极，提高磁法勘探的准确度。

附图说明

图1是本发明实施例提供的一种磁异常化极方法的流程图；

图2为南充-旺苍地区的磁异常资料；

图3为南充-旺苍地区的磁异常垂直梯度与总梯度的互相关系数；

图4为南充-旺苍地区考虑剩磁的变倾角化极磁异常；

图5为石柱-大竹地区的磁异常资料；

图6为石柱-大竹地区的磁异常垂直梯度与总梯度的互相关系数；

图7为石柱-大竹地区考虑剩磁的变倾角化极磁异常。

具体实施方式

请参阅图1，图1是本发明实施例提供的一种磁异常化极方法的流程图。如图1所示，本实施例的方法包括以下步骤：

步骤S1：采集研究区的磁场数据，对采集到的磁场数据进行校正，得到磁异常数据。

本步骤中，具体是通过磁力仪在野外采集研究区的磁场数据，并将获得的磁场数据导入计算机中，以进行各种校正，以得到异常数据。

步骤S2：将磁异常数据网格化，形成规则的磁异常数据t(x,y)，其中x,y分别表示磁异常数据的平面网格点坐标，记为矩阵T(1:m,1:n)，大小为m×n，m表示行数，n表示点数。

本步骤中，具体是通过surfer软件中的kriging插值方法将磁异常数据网格化。

步骤S3：获得研究区的地磁场方向，并将地磁场方向整理成与矩阵T(1:m,1:n)等大的矩阵，记为I(1:m,1:n)。

本步骤中，具体是利用国际地磁参考场，获得研究区的地磁场方向。本发明选用的国际地磁参考场(international geomagnetic reference field,IGRF)优选是第11代国际地磁参考场(IGRF-11)，由国际地磁学与高空物理学联合会(IAGA)于2009年12月发布。

其中，本步骤的将地磁场方向整理成与矩阵T(1:m,1:n)等大的矩阵，记为I(1:m,1:n)的步骤具体包括如下：首先获取地磁倾角范围为a-b，然后设置矩阵进而利用矩阵A的伪逆矩阵计算磁场方向；

P＝pinv(A)×V，pinv表示伪逆矩阵P为计算磁场方向时的参数矩阵，进一步设置网格坐标阵[x,y]＝meshgrid(1:n,1:m)，最后计算矩阵I＝P(1)+P(2)×y+P(3)×x。

步骤S4：对矩阵T(1:m,1:n)进行镶边处理，将矩阵T(1:m,1:n)转变成行数和列数均为2的最小整数幂的新矩阵T’(1-m’:m+m’,1-n’:n+n’)。

本步骤中，优选以“余弦衰减置零扩边法”给矩阵T(1:m,1:n)镶边。具体采用以下公式进行镶边：

其中，T位于T’的中心位置，即T(i,j)＝T’(i,j)，(1≤i≤m,1≤j≤n)。

本步骤中，利用“余弦衰减置零扩边法”给数据进行镶边，很好的压制了边界效应，提高了数据处理的精度。

步骤S5：设计频率域化极因子，并对频率域化极因子做反傅里叶变换，得到化极结果。

本步骤中，频率域化极因子具体为：

式中，ω_x和ω_y是x、y方向的圆波数， 是地磁方向的单位向量，L、M、N是地磁场的方向余弦：cosIcosD,cosIsinD,sinI，I和D分别表示地磁场倾角和地磁场偏角；是剩磁方向的单位向量，l,m,n分别是剩磁方向的方向余弦：cosicosd,cosisind,sini，i和d分别表示剩磁倾角和剩磁偏角。

其中，具体是采用以下的公式做反傅里叶变换：

其中，RTP表示化极结果。

本步骤中，通过频率域化极因子进行运算，频率域内的运算时间短，效率高，在计算机上更易编程实现，降低了方法的适用门槛。

步骤S6：计算化极结果的垂直梯度和总梯度，并计算其互相关系数，以最大互相关系数确定剩磁方向。

本步骤中，利用了磁异常的总梯度是垂直梯度的包络这一物理性质，二者在垂直磁化时具有对称性，当二者的互相关系数最大时，磁异常最接近垂直磁化状态；也就是说，此时的选用的剩磁方向能够使得化极结果最接近垂直磁化，因此该方向就是磁源体的真实剩磁方向。根据步骤S5的频率域化极因子中进行变倾角化极，求取化极结果的垂直梯度和总梯度，并计算其互相关系数，以最大互相关系数来剩磁方向。

步骤S7：根据地磁场方向和剩磁方向，并根据频率域化极因子对矩阵T’和I在频率域内进行变倾角化极处理，计算研究区磁异常的变倾角化极结果。

其中，在本步骤之前，首先从矩阵I中提取平均地磁倾角为平均地磁方向，然后计算出研究区中的实际地磁方向与平均地磁方向的差

本步骤具体利用泰勒公式进行变倾角化极，采用的公式是：

式中，RTP_diff是变倾角化极结果，RTP_mean是使用研究区平均地磁倾角和地磁偏角得到的化极结果。

由于本步骤利用泰勒公式，实现了变倾角化极，避免了用单一地磁倾角进行化极时的失真问题，从数学上解决了纬度跨度较大时研究区磁异常化极的问题，因此提高了磁法勘探的准确度。

本实施例的方法适用于大、中、小比例尺的磁异常的处理。

本发明还分别以南充-旺苍地区和石柱-大竹地区的磁异常资料为例，结合前文所述的方法，给予进一步说明。

首先，以南充-旺苍地区的磁异常资料为例，按照前文所述的方法，对磁异常资料进行处理。首先请参照图2-4。

首先利用磁力仪在野外采集资料，获得南充-旺苍地区的磁场数据资料，然后将获得的磁场资料导入计算机，对数据进行各种校正，得到南充-旺苍地区(南北210km×东西169km)的磁异常资料。然后利用surfer软件中的kriging插值方法将上述磁异常资料按照1km间距进行网格化，形成规则的磁异常资料，记为矩阵T(1:210,1:169)，详见图2。

进一步的，利用国际地磁参考场(IGRF)，获得该地区的地磁场方向。本地区磁异常数据采集自1971年，9月1日，纬度范围约为30.53°N～32.26°N，根据IGRF模型，地磁倾角范围约为44～46.5，地磁偏角约为2°。地磁倾角按照以下算法，整理成与磁异常矩阵T等大的矩阵I。

①设置矩阵

②利用矩阵A的伪逆矩阵计算磁场方向；

③P＝pinv(A)×V，pinv表示伪逆矩阵，P为计算磁场方向时的参数矩阵；

④设置网格坐标阵[x,y]＝meshgrid(1:169,1:210)；

⑤I＝P(1)+P(2)×y+P(3)×x。

进一步的，从矩阵I中提取平均地磁倾角计算出研究区中的实际地磁方向与平均地磁方向的差

进一步的，以“余弦衰减置零扩边法”给矩阵T(1:210,1:169)镶边，将矩阵T(1:210,1:169)转变成行数和列数均为2的最小整数幂的新矩阵T’(-22:233,-43:212)。具体镶边公式如前文所述，在此不再赘述。

设计的频率域化极因子，其中频率域化极因子的公式如前文所述，在此不再赘述。进一步对对频率域化极因子做反傅里叶变换，得到化极结果。具体反傅里叶变换的公式如前文所述，在此不再赘述。

进一步的，在频率域内计算化极结果的垂直梯度和总梯度，并求出其互相关系数C。C是剩磁方向的函数。C最大时对应的剩磁方向就是本地区磁源体的剩磁方向(i,d)。图3是互相关系数的结果，横坐标是偏角，纵坐标是倾角，只显示了最大互相关系数附近的图像。该地区的剩磁方向为剩磁倾角3.3°，剩磁偏角-0.5°。

进一步的，如图4，利用该地区剩磁方向和地磁场方向，通过矩阵T’、I在频率域内进行变倾角化极处理。计算该地区磁异常的变倾角化极结果。该地区地磁场倾角，即地磁场方向为45°左右，而剩磁倾角，即剩磁方向仅为3.3°，二者相差约14倍；对比图2与图4可以发现，磁异常幅值由-240～+420nT变为-300～+800nT，二者相差较大；正磁异常中心位置南充附近，在考虑剩磁影响后，磁异常峰值位于南部-仪陇之间，北偏达77.5km，且是唯一的等值线闭合处。原始磁异常(图2)与化极磁异常(图4)的巨大差异说明剩余磁化在该地区是很强的。

以下以石柱-大竹地区的磁异常资料为例，按照前文所述的方法，对磁异常资料进行处理。首先请参照图5-7。

首先利用磁力仪在野外采集资料，获得石柱-大竹地区的磁场数据资料，然后将获得的磁场资料导入计算机，对数据进行各种校正，得到石柱-大竹地区(南北168km×东西150km)的磁异常资料。然后利用surfer软件中的kriging插值方法将上述磁异常资料按照1km间距进行网格化，形成规则的磁异常资料，记为矩阵T(1:168,1:150)，详见图5。

进一步的，利用国际地磁参考场(IGRF)，获得该地区的地磁场方向。本地区磁异常数据采集自1971年，9月1日，纬度范围约为29.69°N～30.95°N，根据IGRF模型，地磁倾角范围约为43.5～45，地磁偏角约为2°。地磁倾角按照以下算法，整理成与磁异常矩阵T等大的矩阵I。

①设置矩阵

②利用矩阵A的伪逆矩阵计算磁场方向；

③P＝pinv(A)×V，pinv表示伪逆矩阵，P为计算磁场方向时的参数矩阵；

④设置网格坐标阵[x,y]＝meshgrid(1:150,1:168)；

⑤I＝P(1)+P(2)×y+P(3)×x。

进一步的，从矩阵I中提取平均地磁倾角计算出研究区中的实际地磁方向与平均地磁方向的差

进一步的，以“余弦衰减置零扩边法”给矩阵T(1:168,1:150)镶边，将矩阵T(1:168,1:150)转变成行数和列数均为2的最小整数幂的新矩阵T’(-43:212,-52:203)。具体镶边公式如前文所述，在此不再赘述。

设计的频率域化极因子，其中频率域化极因子的公式如前文所述，在此不再赘述。进一步对对频率域化极因子做反傅里叶变换，得到化极结果。具体反傅里叶变换的公式如前文所述，在此不再赘述。

进一步的，在频率域内计算化极结果的垂直梯度和总梯度，并求出其互相关系数C。C是剩磁方向的函数。C最大时对应的剩磁方向就是本地区磁源体的剩磁方向(i,d)。图6是互相关系数的结果，横坐标是偏角，纵坐标是倾角，只显示了最大互相关系数附近的图像。该地区的剩磁方向为剩磁倾角21°，剩磁偏角-41°。

进一步的，如图7，利用该地区剩磁方向和地磁场方向，通过矩阵T’、I在频率域内进行变倾角化极处理。计算该地区磁异常的变倾角化极结果。该地区地磁场倾角为44°左右，而剩磁倾角为21°，二者相差约2倍；对比图5与图7可以发现，磁异常幅值由-180～+200nT变为-320～+280nT，二者相差较小；正磁异常中心位置丰都-石柱之间，在考虑剩磁影响后，磁异常峰值位于垫江，北偏约25km。原始磁异常(图5)与化极磁异常(图7)差异不大，说明剩余磁化在该地区较弱。

综上所述，本发明能够利用地磁场方向和剩磁方向对实测磁异常进行变倾角化极，提高磁法勘探的准确度。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (8)

1.一种磁异常化极方法，其特征在于，包括以下步骤：

采集研究区的磁场数据，对采集到的所述磁场数据进行校正，得到磁异常数据；

将所述磁异常数据网格化，形成规则的磁异常数据t(x,y)，其中x,y分别表示磁异常数据的平面网格点坐标，记为矩阵T(1:m,1:n)，大小为m×n，m表示行数，n表示点数；

获得所述研究区的地磁场方向，并将所述地磁场方向整理成与矩阵T(1:m,1:n)等大的矩阵，记为I(1:m,1:n)；

对所述矩阵T(1:m,1:n)进行镶边处理，将矩阵T(1:m,1:n)转变成行数和列数均为2的最小整数幂的新矩阵T’(1-m’:m+m’,1-n’:n+n’)；

设计频率域化极因子，并对所述频率域化极因子做反傅里叶变换，得到化极结果；

计算化极结果的垂直梯度和总梯度，并计算其互相关系数，以最大互相关系数确定剩磁方向；

根据所述地磁场方向和所述剩磁方向，并根据所述频率域化极因子对所述矩阵T’和I在频率域内进行变倾角化极处理，计算所述研究区磁异常的变倾角化极结果。

2.根据权利要求1中所述的磁异常化极方法，其特征在于，所述将所述磁异常数据网格化的步骤包括：

利用kriging插值方法将所述磁异常数据网格化。

3.根据权利要求1中所述的磁异常化极方法，其特征在于，所述将所述地磁场方向整理成与矩阵T(1:m,1:n)等大的矩阵，记为I(1:m,1:n)的步骤包括：

获取地磁倾角范围为a-b；

设置矩阵

利用矩阵A的伪逆矩阵计算磁场方向；

P＝pinv(A)×V，pinv表示伪逆矩阵P为计算磁场方向时的参数矩阵，其使用方法见下式；

设置网格坐标阵[x,y]＝meshgrid(1:n,1:m)；

计算矩阵I＝P(1)+P(2)×y+P(3)×x。

4.根据权利要求1中所述的磁异常化极方法，其特征在于，对所述矩阵T(1:m,1:n)进行镶边处理的步骤包括：

采用以下公式进行镶边：

$T^{\′}(i,j)=T(i,1)\×cos\frac{(1-j)}{2n^{\′}}\mathrm{\π},(1\≤i\≤m,1-n^{\′}\≤j\≤0)$

$T^{\′}(i,j)=T(i,n)\×cos\frac{(j-n)}{2n^{\′}}\mathrm{\π},(1\≤i\≤m,1\≤j\≤n+n^{\′})$

$T^{\′}(i,j)=T(1,j)\×cos\frac{(1-i)}{2m^{\′}}\mathrm{\π},(1-m^{\′}\≤i\≤0,1-n^{\′}\≤j\≤n+n^{\′})$

$T^{\′}(i,j)=T(m,j)\×cos\frac{(i-m)}{2m^{\′}}\mathrm{\π},(m+1\≤i\≤m+m^{\′},1-n^{\′}\≤j\≤n+n^{\′})$

其中，T位于T’的中心位置，即T(i,j)＝T’(i,j)，(1≤i≤m,1≤j≤n)。

5.根据权利要求1中所述的磁异常化极方法，其特征在于，所述频率域化极因子为：

$\mathrm{\Ψ}=\frac{{\mathrm{\ω}}_x^2+{\mathrm{\ω}}_y^2}{(\stackrel{\→}{K}\·\hat{B})(\stackrel{\→}{K}\·\hat{M})}=\frac{{\mathrm{\ω}}_x^2+{\mathrm{\ω}}_y^2}{\left(i\right({\mathrm{L\ω}}_x+{\mathrm{M\ω}}_y)+N\sqrt{{\mathrm{\ω}}_x^2+{\mathrm{\ω}}_y^2})\left(i\right({\mathrm{l\ω}}_x+{\mathrm{m\ω}}_y)+n\sqrt{{\mathrm{\ω}}_x^2+{\mathrm{\ω}}_y^2})}$

式中，ω_x和ω_y是x、y方向的圆波数， 是所述地磁方向的单位向量，L、M、N是地磁场的方向余弦：cos I cos D,cos I sin D,sin I，I和D分别表示地磁场倾角和地磁场偏角；是剩磁方向的单位向量，l,m,n分别是剩磁方向的方向余弦：cos i cos d,cos i sin d,sin i，i和d分别表示剩磁倾角和剩磁偏角。

6.根据权利要求1-5任一项所述的磁异常化极方法，其特征在于，所述对所述频率域化极因子做反傅里叶变换的步骤包括：

采用以下的公式做反傅里叶变换：

其中，RTP表示化极结果。

7.根据权利要求1中所述的磁异常化极方法，其特征在于，根据所述频率域化极因子对矩阵T’、I在频率域内进行变倾角化极处理的步骤之前包括：

从所述矩阵I中提取平均地磁倾角为平均地磁方向；

计算出研究区中的实际地磁方向与平均地磁方向的差

8.根据权利要求7中所述的磁异常化极方法，其特征在于，利用泰勒公式进行变倾角化极，采用的公式是：

${\mathrm{RTP}}_{diff}={\mathrm{RTP}}_{mean}+\mathrm{\Δ}I\frac{\∂RTP}{\∂I}+0.5{\mathrm{\ΔI}}^2\frac{{\∂}^{2}RTP}{\∂I^2}+0.16666{\mathrm{\ΔI}}^3\frac{{\∂}^{3}RTP}{\∂I^3};$

式中，RTP_diff是变倾角化极结果，RTP_mean是使用研究区平均地磁倾角和地磁偏角得到的化极结果。