一种变倾角三维磁场的快速正演方法

本发明属于地球物理勘探，具体涉及一种变倾角三维磁场的快速正演方法。

背景技术：

1、磁法勘探作为常用的地球物理勘探方法，主要用来寻找和勘探矿产(例如铁矿、铅锌矿、铜锦矿等)、进行地质填图、研究与油气有关的地质构造、研究大地构造以及进行军事侦察(例如未爆弹、潜艇和水雷等)等。反演是解释磁性目标的重要过程，其主要是为了利用获得的磁性目标磁测数据评估地下或水下未知磁性体的轮廓形态。

2、正演作为反演的基础，正演计算的效率和精度直接影响了反演计算的效率和效果。基于bttb(block-toeplitz toeplitz-block)矩阵的重磁正演算法具有高精度、高效率的优势，众多学者对其进行了相关研究，其中袁洋等提出的基于bttb矩阵的快速高精度三维磁场正演算法，只能计算地磁场方向和磁化强度方向固定的情况，正演计算受磁场方向的影响较大。对于小区域并且磁性体仅包含感应磁化而忽略剩余磁化影响的情况，磁法正演时一般采用工区中心的磁参数代表整个工区的磁参数，由于工区内磁参数变化小，正演计算能够达到较高的计算精度。但是，当工作区域较大时或存在强剩磁磁性体时，工区内的磁参数存在较大的变化，如果强行用固定的磁参数，则会导致正演计算结果存有较大的误差。

3、因此，亟需提出一种变倾角三维磁场的快速正演方法，实现对大区域磁场或强剩磁磁场的高精度快速正演。

技术实现思路

1、本发明为了解决上述技术问题，提出了一种变倾角三维磁场的快速正演方法，实现了对大区域磁场和强剩磁磁场的高精度快速正演，有利于准确获取地下空间三维磁场的真实情况。

2、为了实现上述目的，本发明采用如下技术方案：

3、一种变倾角三维磁场的快速正演方法，具体包括如下步骤：

4、步骤1，获取工区的工况信息，确定工区观测平面中所有观测点的位置，对应各观测点的位置将工区地下场源空间划分为多个长方体模型，确定地下场源空间中长方体模型总磁场强度矢量的正演计算公式；

5、步骤2，根据工区的位置和日期，基于国际地磁参考场模型计算各长方体模型中心点处和各观测点位置的磁参数矩阵

6、步骤3，基于工区地下场源空间中长方体模型总磁场强度矢量的正演计算公式，结合国际地磁参考场模型所计算的各长方体模型中心点处和各观测点位置的磁参数矩阵，利用快速算法计算各长方体模型用于确定总磁场强度矢量，得到工区地下场源空间的总磁化强度t。

7、优选地，所述步骤1中，获取工区的工况信息，确定工区内观测点的总数量p以及各观测点的位置，以地表作为观测平面构建三维空间坐标系，沿z方向将工区地下场源空间剖分为p层，每层内将地下场源空间沿x方向将等间隔划分m份、沿y方向等间隔划分n份，划分后地下场源空间共设置有n个长方体模型，长方体模型的总数n为m×n×p，水平方向上各层中各长方体模型的位置与地表观测点的位置一一对应，工区内观测点的总数量p为m×n个；

8、在三维空间坐标系中，通过在观测平面上划分网格获取各观测点的坐标，确定长方体模型总磁场强度矢量的正演计算公式为体积分公式，长方体模型总磁场强度矢量的正演计算公式如公式(1)所示：

9、

10、其中，

11、

12、

13、

14、

15、r＝[(ξ-x)+(η-y)+(ζ-z)]1/2        (1-5)

16、式中，δt为总磁场，(x,y,z)为观测平面上网格点的坐标；(ξ,η,ζ)为长方体模型内场源点的坐标；(x0,y0,z0)为长方体模型中心点的坐标，a为长方体模型沿x方向的延伸长度，b为长方体模型沿y方向的延伸长度，c为长方体模型沿z方向的延伸长度；m为长方体模型的总磁化强度；μ0为真空中的磁导率；i0为总磁化强度倾角；d0为总磁化强度偏角；i1为地磁场方向倾角；d1为地磁场方向偏角；为观测点磁化方向在x方向的单位投影，为观测点磁化方向在y方向的单位投影，为观测点磁化方向在z方向的单位投影；vxx为长方体模型的磁力位在x方向上的二阶导数，vyy为长方体模型的磁力位在y方向上的二阶导数，vzz为长方体模型的磁力位在z方向上的二阶导数，vxy为长方体模型的磁力位在xy方向上的二阶导数，vyx为长方体模型的磁力位在yx方向上的二阶导数，vxz为长方体模型的磁力位在xz方向上的二阶导数，vzx为长方体模型的磁力位在zx方向上的二阶导数，vyz为长方体模型的磁力位在yz方向上的二阶导数，vzy为长方体模型的磁力位在zy方向上的二阶导数；tx为磁场在x方向的分量，ty为磁场在y方向的分量，tz为磁场在z方向的分量；π为圆周率；mx为长方体模型磁化率在x方向的分量，my为长方体模型磁化率在y方向的分量，mz为长方体模型磁化率在z方向的分量；为长方体模型磁化方向在x方向的单位投影，为长方体模型磁化方向在y方向的单位投影，为长方体模型磁化方向在z方向的单位投影；r为长方体内场源点与观测点之间的距离。

17、优选地，所述步骤2中，通过将长方体模型中心点处的位置以及各观测点的位置输入国际地磁参考场模型中，利用国际地磁参考场模型计算得到地下场源空间中各长方体模型中心点处的磁参数矩阵以及各观测点位置的磁参数矩阵

18、优选地，所述步骤3中，当地下场源空间中长方体模型与观测点具有不同的倾角和偏角时，将长方体模型的总磁场强度矢量表示为：

19、

20、式中，为矩阵尺寸p×p的对角矩阵，包括对角矩阵对角矩阵和对角矩阵为矩阵尺寸n×n的二维矩阵，包括对角矩阵对角矩阵和对角矩阵m为磁化率模型矩阵，矩阵尺寸为n×1；磁化率正演系数矩阵vrs，矩阵尺寸为p×n，包括核矩阵vxx、核矩阵vxy、核矩阵vxz、核矩阵vyx、核矩阵vyy、核矩阵vyz、核矩阵vzx、核矩阵vzy和核矩阵vzz；

21、当地下场源空间中长方体模型与观测点具有相同的倾角和偏角时，结合单位矩阵id和im，将长方体模型的总磁场强度矢量表示为：

22、δt＝vm             (3)

23、其中，

24、

25、式中，qx、qy、qz均为与倾角相关的标量。

26、优选地，所述步骤3中，所述磁化率正演系数矩阵vrs中的各核矩阵均为bttb矩阵，具有大量重复的矩阵元素，基于快速算法计算长方体模型的总磁场强度矢量，利用三维矩阵替代磁化率正演系数矩阵vrs，其中r,s∈{x,y,z}，三维矩阵的尺寸为(2m-1)×(2n-1)×p，再分别利用尺寸为m×n×p的三维矩阵替代向量m、尺寸为的m×n×p三维矩阵替代二维矩阵尺寸为m×n×p的三维矩阵替代二维矩阵后，计算地下场源空间中各长方体模型的总磁场强度矢量，具体包括以下步骤：

27、步骤3.1，当地下场源空间中长方体模型与观测点具有不同的倾角和偏角时，将等同于通过将长方体模型总磁场强度矢量计算公式中的同阶矩阵中的元素对应相乘，得到：

28、

29、式中，mt为计算磁化率分量，t∈{x,y,z}，计算磁化率分量mt包括矩阵mx、矩阵my和矩阵mz；

30、步骤3.2，将公式(5)中的trs＝vrsmt等同于其中，r,s,t∈{x,y,z}，k为层序号，*为卷积计算运算符，基于快速算法，通过快速傅里叶正变换和逆变换进行卷积计算，快速确定得到：

31、

32、式中，f为快速傅里叶正变换，f-1为快速傅里叶逆变换；

33、步骤3.3，将公式(6)中的等同于结束对工区地下场源空间的变倾角三维磁场快速正演，确定工区地下场源空间的总磁场强度t。

34、本发明所带来的有益技术效果：

35、本发明提出了一种变倾角三维磁场的快速正演方法，解决了现有变倾角三维磁场正演方法计算效率低下，导致难以开展大规模数据快速计算的问题。本发明方法通过将地下场源空间划分为多个长方体模型，将各长方体模型和观测点单独设置磁参数，通过快速且准确的获取各长方体模型的总磁场强度，提高了变倾角三维磁场快速正演的计算效率和计算精度，为地下或水下未知磁性体的探测评估奠定了基础。