任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法与流程

本发明涉及磁场梯度张量数值模拟方法
技术领域：
，具体的涉及一种任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法。
背景技术：
：磁法勘探是地球物理勘探中一种重要的技术手段，在区域地质调查、普查找矿工作中发挥着重要作用。该方法通过测量由岩石、矿石等磁性体的磁化率差异所引起的磁异常，对测量数据进行反演、解释，进而研究地质构造和矿产资源的分布规律。随着传感器技术的不断发展，磁场梯度张量测量技术日渐成熟。相比磁异常和磁场三分量，磁场梯度张量具有更高的分辨率。利用磁场梯度张量数据进行反演、解释，是磁法勘探的重要研究方向。磁性体磁场梯度张量的数值模拟是进行反演的基础。任意磁化率分布复杂磁性体磁场梯度张量快速、高精度数值模拟一直是个难点问题。针对磁场数值模拟，众多国内外学者进行了研究。数值模拟首先对研究区域进行剖分，然后根据剖分方式，采用某种方法计算磁场。文献(姚长利,郝天珧,管志宁,张聿文.重磁遗传算法三维反演中高速计算及有效存储方法技术.地球物理学报,2003.46(2):252-258.)采用结构化剖分方式，根据离散化后数学问题的特点，提出了“格架分离”技术和“格架等效计算方案”，较好解决了计算效率和计算精度问题，但对于大规模剖分情形，该文献所给出的数值模拟方法的计算效率仍然比较低；文献(Tontini,F.C.,L.Cocchi,C.Carmisciano.Rapid3-DforwardmodelofpotentialfieldswithapplicationtothePalinuroSeamountmagneticanomaly(southernTyrrhenianSea,Italy).JournalofGeophysicalResearch,2009.114.)采用结构化剖分方法，采用三维傅里叶变换，给出了任意密度或者磁化率分布情形下重磁数值模拟的波数域表达式，借助三维快速傅里叶变换算法，实现了快速数值模拟，该方法效率极高，但为克服截断效应，使用该方法前需要对剖分区域进行扩边，影响了数值模拟精度；文献(Wu,L.,Tian,G.High-precisionFourierforwardmodelingofpotentialfields.Geophysics,2014,79(5):G59-G68.)提出了一种重磁数值模拟的Gauss-FFT方法，该方法有效克服了传统傅里叶变换方法的截断效应问题，提高了数值模拟精度，但计算效率降低，且计算场源内部磁场时误差比较大。剖分方式和计算方法共同决定了数值模拟的效率和精度。数值模拟的效率和精度是一对矛盾体，目前已有的数值模拟方法存在的最大问题是不能同时保证计算效率和精度，无法满足大规模数据三维磁化率精细反演、人机交互建模和解释的需求；同时，目前的方法不能准确计算场源内部的磁场梯度张量，不能满足井地联合反演的需求。因此，寻找一种计算效率高、同时能保证计算精度、且能精确计算全空间磁场梯度张量的数值模拟方法具有重要的现实意义。技术实现要素：本发明的目的在于提供一种任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法，该发明解决了目前磁场梯度张量数值模拟方法计算精度低、计算时间长，不能准确计算场源内部磁场，无法满足大规模观测数据精细反演、人机交互建模解释和井地联合反演需求的技术问题。本发明提供了一种任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法，包括以下步骤：步骤S100：建立包含目标区域的棱柱体模型，将棱柱体模型剖分为多个棱柱体，对各棱柱体设置磁化率，并计算各棱柱体的磁化强度得到组合棱柱体模型；步骤S200：按公式(5)采用二维离散卷积法计算各层棱柱体的磁场梯度张量其中，表示第r层(r＝1,2,…,L)棱柱体在高度面z0产生的磁场梯度张量，my(ξp,ηq,ζr)为坐标为(ξp,ηq,ζr)的棱柱体在x方向的磁化强度分量、my(ξp,ηq,ζr)为坐标为(ξp,ηq,ζr)的棱柱体在y方向的磁化强度分量，mz(ξp,ηq,ζr)为坐标为(ξp,ηq,ζr)的棱柱体在z方向的磁化强度分量，hxz(xm-ξp,yn-ηq,z0-ζr)为磁化强度XZ分量的加权系数、hyz(xm-ξp,yn-ηq,z0-ζr)为磁化强度YZ分量的加权系数、hzz(xm-ξp,yn-ηq,z0-ζr)为磁化强度ZZ分量的加权系数，M为目标区域x方向棱柱体的剖分个数，N为目标区域y方向棱柱体的剖分个数；步骤S300：按公式(16)累加各层棱柱体的磁场梯度张量得到组合棱柱体模型的模拟磁场梯度张量其中，L表示目标区域z方向棱柱体剖分个数。进一步地，加权系数按公式(6)～(8)计算：其中，(xm,yn,z0)为观测点坐标，z0为常值；x1＝ξp-0.5Δx-xm，x2＝ξp+0.5Δx-xm，y1＝ηq-0.5Δy-yn，y2＝ηq+0.5Δy-yn，z1＝ζr-0.5Δz-z0，z2＝ζr+0.5Δz-z0，μijk＝(-1)i(-1)j(-1)k，i＝1,2，j＝1,2，k＝1,2。进一步地，棱柱体模型为规则棱柱型。进一步地，根据目标区域的磁化率分布将每个棱柱体的磁化率设置为常值，并将位于目标区域空气部分的棱柱体的磁化率设为零。进一步地，棱柱体磁场强度的计算方法包括以下步骤：步骤S110：根据地球主磁场模型IGRF，计算棱柱体中心点的地球主磁场X轴的分量Tx(ξp,ηq,ζr)、Y轴的分量Ty(ξp,ηq,ζr)、Z轴的分量Tz(ξp,ηq,ζr)；步骤S120：按式(1)(2)(3)计算磁化强度mx(ξp,ηq,ζr)＝χ(ξp,ηq,ζr)Tx(ξp,ηq,ζr)(1)my(ξp,ηq,ζr)＝χ(ξp,ηq,ζr)Ty(ξp,ηq,ζr)(2)mz(ξp,ηq,ζr)＝χ(ξp,ηq,ζr)Tz(ξp,ηq,ζr)(3)其中，(ξp,ηq,ζr)表示目标区域中编号为(p,q,r)的棱柱体几何中心坐标，χ(ξp,ηq,ζr)表示该棱柱体的磁化率值，Tx(ξp,ηq,ζr)表示(ξp,ηq,ζr)处地球主磁场的X轴上的分量、Ty(ξp,ηq,ζr)表示(ξp,ηq,ζr)处地球主磁场的Y轴上的分量、Tz(ξp,ηq,ζr)表示(ξp,ηq,ζr)处地球主磁场的Z轴上的分量，mx(ξp,ηq,ζr)表示(ξp,ηq,ζr)处磁化强度的X轴上的分量，my(ξp,ηq,ζr)处磁化强度的Y轴上的分量、mz(ξp,ηq,ζr)处磁化强度的Z轴上的分量。本发明的技术效果：1、本发明提供的任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法，通过建立包含所有目标区域的规则棱柱体模型，使得目标区域(包含起伏地形)能完全嵌入在该棱柱体模型中，能很容易地刻画任意磁化率分布复杂磁性体以及起伏地形，从而提高后续模型模拟的计算精度；而且所构建模型剖分方法简单、灵活，能提高模拟效率。2、本发明提供的任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法，能够实现磁场梯度张量的快速、高精度数值模拟，可以满足大规模数据三维磁化率精细反演、人机交互建模和解释的需求；3、本发明提供的任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法，用于处理大规模数值模拟时，不但计算效率和计算精度高，并且所需计算机内存小；4、本发明提供的任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法，不仅可以应用于磁性体感应磁化产生的磁场梯度张量数值模拟，还可以应用于磁性体剩余磁化产生的磁场梯度张量数值模拟；5、本发明提供的任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法，不但可以精确模拟场源外部的磁场梯度张量，还可以精确模拟场源内部的磁场梯度张量，有利于开展井地联合反演的研究。具体请参考根据本发明的食品泵轴提出的各种实施例的如下描述，将使得本发明的上述和其他方面显而易见。附图说明图1为本发明提供任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法流程示意图；图2为本发明提供任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法优选实施例流程图；图3为本发明优选实施例中所处理复杂磁性体的模型表示图；图4为本发明优选实施例中磁性体剖分后所得棱柱磁性体剖面图；图5为本发明优选实施例的磁场梯度张量计算值；图6为本发明优选实施例的磁场梯度张量理论值；图7为本发明优选实施例中x＝0m观测线对应理论值与采用本发明提供方法所得计算值对比示意图。图例说明：L：表示z方向剖分小棱柱体个数；M：表示x方向剖分小棱柱体个数；N：表示y方向剖分小棱柱体个数；χ：表示磁化率。具体实施方式构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。参见图1，本发明提供的任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法包括以下步骤：步骤S100：建立包含目标区域的棱柱体模型，将棱柱体模型剖分为多个棱柱体，对各棱柱体设置磁化率，并计算各棱柱体的磁化强度得到组合棱柱体模型；复杂磁性体模型表示：建立包含目标区域的规则棱柱体模型，使得目标区域(包含起伏地形)完全嵌入在该棱柱体模型中；将该棱柱体模型划分成许多小棱柱体，每个小棱柱体磁化率为常值，不同棱柱体磁化率取值不同，以此刻画任意磁化率分布复杂磁性体模型；将位于空气部分的小棱柱体的磁化率值设为零，以此刻画起伏地形。采用这种剖分形式结合后续所用步骤能在保证剖分过程简单快捷的情况下，保证计算结果的精确性。有利于提供计算效率。磁化强度计算：可以按现有方法进行计算，优选的包括以下步骤：步骤S110：根据公布的地球主磁场模型IGRF，计算剖分所得各棱柱体中心点的地球主磁场三分量Tx(ξp,ηq,ζr),Ty(ξp,ηq,ζr),Tz(ξp,ηq,ζr)；步骤S120：计算磁化强度mx(ξp,ηq,ζr)＝χ(ξp,ηq,ζr)Tx(ξp,ηq,ζr)(1)my(ξp,ηq,ζr)＝χ(ξp,ηq,ζr)Ty(ξp,ηq,ζr)(2)mz(ξp,ηq,ζr)＝χ(ξp,ηq,ζr)Tz(ξp,ηq,ζr)(3)式(1)(2)(3)中，(ξp,ηq,ζr)表示目标区域中编号为(p,q,r)的棱柱体几何中心坐标；χ(ξp,ηq,ζr)表示该棱柱体的磁化率值；Tx(ξp,ηq,ζr)、Ty(ξp,ηq,ζr)、Tz(ξp,ηq,ζr)分别表示(ξp,ηq,ζr)处地球主磁场的x、y、z分量；mx(ξp,ηq,ζr)、my(ξp,ηq,ζr)、mz(ξp,ηq,ζr)分别表示(ξp,ηq,ζr)处磁化强度的x、y、z分量。按此方法计算磁化强度，计算效率得到提高。为了提高棱柱体计算精度，本发明提供了组合棱柱体模型磁场梯度张量计算公式为其中，L表示z方向棱柱体剖分个数；M表示x方向棱柱体剖分个数；N表示y方向棱柱体剖分个数。本发明提出的公式(4)，为计算组合棱柱体模型产生的场的精确公式，因而本发明提供的方法具有较高的计算精度。同时配合前述剖分棱柱体组合模型，能实现对复杂形状磁性体的刻画准确，从而提高了场就计算准确。步骤S200，按公式(5)采用二维离散卷积法计算各层棱柱体的磁场梯度张量其中，表示第r层(r＝1,2,…,L)棱柱体组合模型在高度面z0产生的磁场梯度张量，my(ξp,ηq,ζr)、my(ξp,ηq,ζr)、mz(ξp,ηq,ζr)为坐标为(ξp,ηq,ζr)的棱柱体在x方向、y方向、z方向的磁化强度分量，hxz(xm-ξp,yn-ηq,z0-ζr)、hyz(xm-ξp,yn-ηq,z0-ζr)、hzz(xm-ξp,yn-ηq,z0-ζr)分别为对应磁化强度三分量的加权系数，M表示目标区域x方向棱柱体剖分个数，N表示目标区域y方向棱柱体剖分个数。优选的，加权系数计算：加权系数计算公式为其中，(xm,yn,z0)表示观测点坐标，z0为常值；其它符号含义如下x1＝ξp-0.5Δx-xm，x2＝ξp+0.5Δx-xm，y1＝ηq-0.5Δy-yn，y2＝ηq+0.5Δy-yn，z1＝ζr-0.5Δz-z0，z2＝ζr+0.5Δz-z0，μijk＝(-1)i(-1)j(-1)k，i＝1,2，j＝1,2，k＝1,2。采用二维离散卷积计算方法来计算任一层(相对z方向而言)棱柱体组合模型磁场梯度张量。此处采用二维离散卷积计算方法配合所提出公式(4)能有效提高该方法的整体处理效率。二维离散卷积计算方法现有技术中主要用于处理数字图像处理，首次用于处理磁性体模拟计算。二维离散卷积计算方法可以按现有方法步骤进行，优选的，式(5)中包含三个二维离散卷积，步骤S200中的二维离散卷积快速计算方法，包括以下步骤：(1)将加权系数h(x1-ξp,y1-ηq,z0-ζr)排列成矩阵t，记为式(9)中，矩阵元素ti,j与加权系数h(x1-ξp,y1-ηq,z0-ζr)存在关系ti,j＝h(x1-ξi+1,yj+1-ηq,z0-ζr)(10)(2)将矩阵t补零扩展成矩阵将矩阵分成四个块矩阵，记为将块矩阵互换位置，得到矩阵cext(3)将第r层对应加权系数h(x1-ξp,y1-ηq,z0-ζr)的磁化强度分量m(ξp,ηq,ζr)，其中p＝1,2,…,M，q＝1,2,…,N，排列成矩阵g，矩阵元素gi,j与密度值存在关系gi,j＝m(ξi,ηj,ζr)(14)将矩阵g补零扩展成矩阵gext式(15)中，0M×N表示M×N零矩阵；(4)计算式中，fft2()表示二维快速傅里叶变换；(5)计算式中，“.*”表示对应元素相乘运算；(6)计算式中，ifft2()表示二维快速傅里叶反变换；(7)提取矩阵fext的前M行前N列，构成矩阵f，当矩阵f中h分别取hxz，hyz，hzz，而m分别对应取mx，my，mz时，即可得到式(5)中三个二维离散卷积：的累加结果。步骤S300，将各层棱柱体组合模型磁场梯度张量进行累加，得到整个组合模型的磁场梯度张量，即本发明提供方法：(1)通过对模型剖分形状和方式进行改进，给出了精确的棱柱体组合模型磁场梯度张量计算公式(4)，从而在保证处理效率的同时提高了计算精度；(2)本发明提供方法通过采用计算效率较高的二维离散卷积算法，进一步提高了处理效率。基于上述基础，本发明提供方法实现了任意磁化率分布复杂磁性体磁场梯度张量数值模拟在计算效率和计算精度上的统一。本发明提供的方法还能用于处理其余分量，以上仅为磁场梯度张量分量Bzz的处理说明。参见图2，下面结合附图对本发明提供的方法作进一步说明。1、复杂磁性体模型表示：首先，建立包含所有目标区域的规则棱柱体模型，确定棱柱体在x，y，z方向的起始位置，使得目标区域(包含起伏地形)完全嵌入在该棱柱体模型中；其次，根据实际问题需求，将棱柱体划分成许多规则小棱柱体(如图3所示)，确定小棱柱体的几何尺寸Δx，Δy，Δz；最后，根据目标区域的磁化率分布，对每个小棱柱体磁化率进行赋值，位于空气部分的小棱柱体，其磁化率设为零；2、磁化强度计算首先，根据公布的地球主磁场模型IGRF，计算目标区域的主磁场三分量Tx(ξp,ηq,ζr),Ty(ξp,ηq,ζr),Tz(ξp,ηq,ζr)；其次，计算磁化强度mx(ξp,ηq,ζr)＝χ(ξp,ηq,ζr)Tx(ξp,ηq,ζr)(1)my(ξp,ηq,ζr)＝χ(ξp,ηq,ζr)Ty(ξp,ηq,ζr)(2)mz(ξp,ηq,ζr)＝χ(ξp,ηq,ζr)Tz(ξp,ηq,ζr)(3)其中，(ξp,ηq,ζr)表示编号为(p,q,r)的小棱柱体几何中心坐标；χ(ξp,ηq,ζr)表示该棱柱体的磁化率值；Tx(ξp,ηq,ζr)、Ty(ξp,ηq,ζr)、Tz(ξp,ηq,ζr)分别表示(ξp,ηq,ζr)处地球主磁场的x、y、z分量；mx(ξp,ηq,ζr)、my(ξp,ηq,ζr)、mz(ξp,ηq,ζr)分别表示(ξp,ηq,ζr)处磁化强度的x、y、z分量；3、加权系数计算：加权系数计算公式为其中，(xm,yn,z0)表示观测点坐标，z0为常值；其它符号含义如下x1＝ξp-0.5Δx-xm，x2＝ξp+0.5Δx-xm，y1＝ηq-0.5Δy-yn，y2＝ηq+0.5Δy-yn，z1＝ζr-0.5Δz-z0，z2＝ζr+0.5Δz-z0，μijk＝(-1)i(-1)j(-1)k，i＝1,2，j＝1,2，k＝1,24、组合棱柱体模型磁场梯度张量计算磁场梯度张量计算公式为其中，L表示z方向棱柱体剖分个数；M表示x方向棱柱体剖分个数；N表示y方向棱柱体剖分个数；首先，采用二维离散卷积快速计算方法来计算一层(相对z方向而言)棱柱体组合模型磁场梯度张量，其计算公式为：其中，表示第r层(r＝1,2,…,L)棱柱体组合模型在高度面z0产生的磁场梯度张量；其次，将各层棱柱体组合模型磁场梯度张量进行累加，得到整个组合模型的磁场梯度张量，即其中包含三个二维离散卷积，以第一个为例，步骤4中所述的二维离散卷积快速计算方法，包括以下步骤：(1)将加权系数hxz(x1-ξp,y1-ηq,z0-ζr)排列成矩阵t，记为式(9)中，矩阵元素ti,j与加权系数hxz(x1-ξp,y1-ηq,z0-ζr)存在关系ti,j＝h(x1-ξi+1,yj+1-ηq,z0-ζr)(10)(2)将矩阵t补零扩展成矩阵将矩阵分成四个块矩阵，记为将块矩阵互换位置，得到矩阵cext(3)将第r层磁化强度mx(ξp,ηq,ζr)(p＝1,2,…,M，q＝1,2,…,N)排列成矩阵g，矩阵元素gi,j与密度值存在关系gi,j＝m(ξi,ηj,ζr)(14)将矩阵g补零扩展成矩阵gext式(15)中，0M×N表示M×N零矩阵；(4)计算式中，fft2()表示二维快速傅里叶变换；(5)计算式中，“.*”表示对应元素相乘运算；(6)计算式中，ifft2()表示二维快速傅里叶反变换；(7)提取矩阵fext的前M行前N列，构成矩阵f，即为二维离散卷积计算结果。参见图4～7，以下结合实例对本发明提供方法所得结果进行分析说明：为了说明本发明所提出的方法用于计算任意磁化率分布复杂磁性体磁场梯度张量时的效率和精度，设计了如下模型算例(图4所示)：目标区域有一个棱柱型磁性体，目标区域范围为：x方向从-1000m到1000m，y方向从-1000m到1000m，z方向从0m到500m(z轴向下为正)；磁性体展布范围为：x方向从-500m到500m，y方向从-500m到500m，z方向从100m到400m；磁化率为0.01(SI)；目标区域地球主磁场为50000nT，磁偏角为0度，磁倾角为90度。将目标区域剖分成400×400×200个大小相同的小棱柱体，计算高度为150m平面(图4中虚线所示，该平面经过场源)上的磁场梯度张量，计算点个数为400×400。本发明提供方法采用Fortran语言编程实现，运行程序所用的个人台式机配置为：CPU为i7-2620，主频为2.7GHz，内存为32GB，四核八线程。运行所需时间约为7秒，即可得的模拟结果，由此可见本发明提供方法计算效率高。磁场梯度张量模拟算法计算值和理论值分别如图5～6所示，两者形态一致，说明本发明提供的模拟方法精度较高。图7所示为x＝0m观测线上计算值和理论值对比，两者吻合度较高。用理论值减去本发明提供方法所得计算值，并对所得差值进行统计，统计结果列于表1中。表1本发明提供方法所得计算值与磁场梯度张量理论值的误差统计表(单位：nT/m)最大值最小值均值均方值理论值1.9-1.91.7*10-160.4误差6.6*10-13-6.4*10-136.1*10-162.9*10-13可知本发明提供方法的模拟结果与理论值误差较小，精度高。且由于所处理实例中所处理的观测面通过场源内部，计算得结果即为场源内部的场，因而本发明提供的方法也能精确计算场源内部磁场梯度张量。本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。通过对附图，说明书和权利要求书的研究，在实施本发明时本领域技术人员可以理解和实现所公开的实施例的变形。在权利要求书中，术语“包括”不排除其他步骤或元素，而不定冠词“一个”或“一种”不排除多个。在彼此不同的从属权利要求中引用的某些措施的事实不意味着这些措施的组合不能被有利地使用。权利要求书中的任何参考标记不构成对本发明的范围的限制。当前第1页1 2 3