一种重磁自约束三维反演与联合解释方法与流程

本发明属于地质勘探技术领域，具体涉及一种重磁自约束三维反演与联合解释方法。

背景技术：

隐伏的深部地质构造对于矿产资源勘查、石油天然气勘查以及大型工程建设等都具有重要的意义。重磁勘探数据具有一定的深部信息，根据重磁勘探数据获取地下三维密度和磁性，从而对地质构造进行解释和划分，对于三维地质建模和透明化、预测油气资源和深部矿体的位置及赋存空间、大型工程建设选址等都极为重要。

重磁反演就是根据观测数据获取地下空间的密度和磁性。其中，物性反演将地下空间离散为若干个网格单元，只求解各单元相应的密度或磁性，它易于模拟复杂的地质体，成为三维反演的重要方向。目前，重磁三维物性反演主要存在两个方面的问题：其一，反演可归结为离散不适定问题，多解性和不稳定性严重；其二，随着观测数据量的增大，与之相应的计算量呈几何级增长，计算成本巨大。此外，在得到地下密度和磁性后，如何进行快速的联合解释，也是急需解决的问题之一。

围绕上述问题，国内外学者做了大量相关研究。在减少解的非唯一性与提高解的稳定性方面，主要基于tikhonov正则化求解不适定问题，并在反演中尽可能地利用先验信息对解模型施加约束，以及利用位场本身包含的潜在信息对解模型进行自约束反演等，但tikhonov正则化涉及到的计算量大，且需显示地存储大型稠密核矩阵，计算成本大。在节省计算成本方面，主要有并行计算、核矩阵压缩与重构、快速正演计算等方法。在密度与磁性联合解释方面，通常切割为许多剖面进行解释，这种做法固然有效，但毋庸置疑的是相应的工作量会很大。

技术实现要素：

为了解决现有技术存在的上述问题，本发明目的在于提供一种重磁自约束三维反演与联合解释方法。

本发明所采用的技术方案为：一种重磁自约束三维反演与联合解释方法，包括以下步骤：

g1.首先对目标地区的地下空间进行离散化处理剖切为多个网格单元；

g2.分别计算每个网格单元对于重力和磁法数据的mic值，将相应的mic值和深度加权系数作为自约束函数，分别形成重力和磁法反演的自约束核矩阵；

g3.通过正则化方法分别求解重力和磁法反演问题，得到重力和磁法的反演结果；

g4.根据重力和磁法反演结果，计算每个网格单元的密度和磁性；

g5.再计算沿深度方向上的密度和磁性之间的mic值，并根据mic值的平面分布规律进行分析和地质解释。

进一步的，所述步骤g2中因将目标地区的地下空间分为多个网格单元，故反演模型如下：

am＝d

式中：a为核矩阵；m为待求解的向量，即密度或磁性；d为观测数据向量。

反演的不适定性表现在：核矩阵是奇异和病态的；观测数据不可避免地存在误差等。使得上式直接解不唯一且不稳定，这就需要特殊的求解方法，即正则化方法。目前，求解不适定问题的正则化方法可分为直接正则化(如：tikhonov)与迭代正则化(如：lsqr)。

进一步的，所述步骤g2中的自约束函数表达式如下：

d＝dmicddepth

式中：所述dmic＝diag{1/mic1,1/mic2,λ,1/micm}，而micj为对应的网格单元的mic值；所述ddepth＝diag{1/(z1)^β/2,1/(z2)^β/2,λ,1/(zm)^β/2}，zj即为对应的网格单元的中心埋深。

其中，深度加权函数的表达式如下：

ddepth＝diag{1/(z1)^β/2,1/(z2)^β/2,λ,1/(zm)^β/2}

式中：所述ddepth为关于深度的自约束函数；而zj为对应单元的中心埋深。

进一步的，所述步骤g2中的反演模型如下：

a^*m^*＝d

式中：所述a^*＝ad^-1；所述m^*＝dm。

进一步的，所述步骤g3中采用lsqr正则算法进行计算。

目前，重磁三维反演绝大多数都基于tikhonov正则化。对于大数据量的重磁三维反演而言，tikhonov正则化会耗费大量的计算成本，主要表现在：每反演计算一次，都需计算矩阵的转置、多次矩阵的乘积和矩阵求逆等；每反演迭代一次，都需通过多次计算选择最佳正则化参数，而每个正则化参数的计算又相当于一次反演计算；需显式地存储大型稠密核矩阵。

lsqr法为基于krylov子空间投影法，是求解大型线性方程组问题的一种有效途径，同时它又具有天然的正则化性质，迭代次数即为正则化参数，具有对计算机内存要求低、迭代收敛快和求解结果稳定等优点。hilber、pascal矩阵以及地震层析成像数值实验结果表明，lsqr法较tikhonov正则化稳定可靠。因此，采用lsqr法求解不适定问题，其求解思路是首先把任意系数矩阵方程化为系数矩阵为方阵的方程，然后利用lanczos方法，求解方程的最小二乘解。利用lsqr法求解重磁反问题，主要有以下两个方面的特点：

1)反演时只需根据l曲线的拐点，返回拐点所对应迭代次数的求解结果即可，避免了正则化参数选择而额外耗费的计算成本；

2)lsqr算法仅涉及到矩阵与向量的乘积，大型稠密核矩阵可不用显式的表示出来，这对于大数据量的重磁三维反演是有利的；

进一步的，所述步骤g2-g4的具体步骤如下：

首先计算核矩阵元素a(i,j)和micj；

然后计算约束核矩阵a^*，根据a^*(i,j)＝a(i,j)·lj计算得到a^*，其中

基于lsqr算法求解不适定a^*m^*＝d，根据根据l曲线的拐点，返回对应迭代次数的解；

网格单元循环，根据式反求模型量，输出求解结果m，反演终止。

这里引入了mic系数进行自约束，因为反演具有多解性和不稳定性，在反演中尽可能地利用先验信息(包括地质、钻孔以及地震资料等)对解模型施加约束，可使得解更加符合实际情况，但在有些情况下，先验信息并不充足。paoletti等指出利用场源本身包含的潜在信息对解模型进行自约束，可提高反演结果的可靠性。事实上，地下空间离散化的网格单元与观测数据存在某种相关性，相关性越强，表明该单元对观测数据的贡献可能性越大，因此可利用相关系数值对网格单元进行加权约束。

由mic定义可知，某个网格单元与观测数据之间的mic值越大，表示该网格单元与观测数据的相关性越强，也就是说该网格单元对重磁异常的贡献可能性越大，反演时，该网格单元的权重相应的也就越大，从而利用mic进行自约束。写为

dmic＝diag{1/mic1,1/mic2,λ,1/micm}

考虑到位场“趋肤效应”，令自约束函数

d＝dmicddepth

式中：ddepth＝diag{1/(z1)^β/2,1/(z2)^β/2,λ,1/(zm)^β/2}，zj为对应单元中心埋深。

由于自约束矩阵d为对角矩阵，有d^-1d＝i(i为单位矩阵)，则am＝d式可写为：

ad^-1dm＝d

令a^*＝ad^-1，m^*＝dm，有：

a^*m^*＝d

上式中的矩阵a^*不同于核矩阵a，它包含了自约束信息，即自约束核矩阵，无疑使得信息更加丰富，有助于提高反演结果的可靠性。

此时，可利用迭代正则化方法lsqr法求解a^*m^*＝d，并根据l曲线的拐点返回求解结果m^*，再根据式m＝d^-1m^*反求模型量，舍弃了传统的tikhonov正则化，节省了大量计算成本。

进一步的，所述步骤g5的具体步骤如下：

(5.1)按照重磁反演的每个水平网格单元，取其沿深度方向上的密度和磁性，计算mic值；

(5.2)最后利用水平网格单元的平面坐标和mic值成图，划分地质体和断裂构造。

重磁三维反演结果提供的信息量相比二维反演大很多，在密度与磁性联合解释方面，按照传统的做法切割为多剖面进行解释固然奏效，但毋庸置疑的是相应的工作量也会很大。如能有效的提取和挖掘两种属性数据之间潜在的信息，可提供更加丰富的解译标识，近而提高联合解释的效率。

一般而言，主要岩石及地层的密度与磁性表现为较强的相关性，如沉积岩(低密度、低磁性)、花岗岩(低密度、低磁性)、基性超基性岩(高密度、高磁性)等；而断裂带则由于岩石比较破碎往往表现为低密度，其磁性则没有明显变化或因矿化蚀变而增高，也就是说断裂带的密度与磁性相关性弱或不相关。因此可以利用这种相关关系作为划分地质体和断裂的标识之一。

具体做法为：按照重磁反演的每个水平网格单元，取其沿深度方向上的密度和磁性，计算两者之间的mic值，然后根据水平网格单元的平面坐标和mic值成图，来划分地质体和断裂构造。当地质体没有遭受明显错断时，mic值较大且均一；而断裂带标识有：线性延伸的条带状、串珠状低mic值异常，具有环状的特征低mic值异常。

进一步的，所述mic值的计算式如下：

式中，xy<b(n)指网格分割细度小于b(n)；其中m(d)x，y为特征矩阵。

进一步的，所述m(d)x，y的表达式如下：

式中，i^*(d,x,y)是给定x列、y行情况下的最大互信息。

本发明的有益效果为：

(1)本发明引入大数据分析挖掘中的mic，反演前计算网格单元与观测数据之间的mic值，并将其作为自约束条件以提升反演结果的可靠性；

(2)本发明还将约束模型替换为约束核矩阵，利用lsqr算法求解反问题以节省计算成本。

(3)通过计算密度和磁性之间的mic值来挖掘两种属性数据之间的规律，以达到快速联合解释的目的。

附图说明

图1是本发明的步骤示意图；

图2是本发明实施例的合成模型结构图；

图3是本发明实施例中利用传统算法的反演结果图；

图4是本发明实施例中利用实施例2中的反演方法的反演结果图；

图5是本发明实施例中某矿区重力三维反演结果示意图；

图6是本发明实施例中的密度磁性的mic值平面图及其解释图；其中a:mic值平面图，图中灰色表示mic值大，黑色表示mic值小；b：断裂构造解译图；

图7是本发明反演的步骤示意图。

具体实施方式

下面结合附图及具体实施例对本发明做进一步阐释。

实施例：

本实施例提供一种重磁自约束三维反演与联合解释方法，如图1所示，具体步骤如下：

对目标区域的地下空间进行离散化处理，剖分为多个网格单元；分别计算每个网格单元对于重力和磁法数据的mic值，将相应的mic值和深度加权系数作为约束条件，分别形成重力和磁法反演的自约束核矩阵；通过lsqr方法分别求解重力和磁法反演问题，得出重力和磁法反演结果；根据重力和磁法反演结果，计算每个网格单元的密度和磁性；计算沿深度方向上的密度和磁性之间的mic值，根据mic值的平面分布规律进行分析和地质解释。

本实施例中以提升反演结果的可靠性、节省反演计算成本和快速进行联合解释为目的，引入大数据分析挖掘中的mic(最大信息数)，将其作为自约束条件以提升反演结果的可靠性，将约束模型替换为约束核矩阵，利用lsqr算法求解反问题以节省计算成本，通过计算密度和磁性之间的mic值来挖掘两种属性数据之间的规律，以达到快速联合解释的目的。最后，通过与传统做法的计算成本、合成模型反演结果的比较，以及实际应用效果来展示本发明专利的良好性能。

如图1所示，具体的反演步骤如下：

(1)将地下空间离散化为m个网格单元；

(2)计算核矩阵元素a(i,j)和micj；

(3)计算约束核矩阵a^*，其中a^*(i,j)＝a(i,j)·lj

(4)基于lsqr算法求解不适定a^*m^*＝d，根据l曲线的拐点，返回对应迭代次数的解；

(5)网格单元循环，根据式反求模型量，输出求解结果m。

而具体的根据反演结果联合解释方法如下：

(1)按照重磁反演的每个水平网格单元，取其沿深度方向上的密度和磁性，计算mic值；

(2)最后利用水平网格单元的平面坐标和mic值成图，划分地质体和断裂构造。

根据上述方法进行效果演算，其中先比较传统的tikhonov正则化与lsqr算法之间的区别。

如图7所示，用一个合成模型来展示本专利中的自约束反演算法效果，它由两个异常体组成，如图3所示，其左侧异常体剩余磁化强度为0.8a/m，右侧异常体剩余磁化强度为1a/m。假设总磁化强度方向的倾角(45°)和相对x轴偏角(50°)与地磁场一致。将数据加入10％高斯噪声，然后利用传统反演算法和本专利中的自约束反演算法进行反演计算。

由图4和图5可见，两种反演算法的求解结果差异较大，主要表现在以下两个方面：其一，求解的磁化强度范围差异大，传统反演算法求解范围在-0.1～0.4a/m之间(图4)，而自约束反演算法求解范围在-0.1～1a/m之间(图5)，后者与合成模型的磁化强度范围更加接近；其二，成像效果上存在差异，传统反演算法结果比较模糊，尤其在深部很难识别出由两个异常体组成，反演异常底部为发散的(图4)，自约束反演算法结果则明显与合成模型更加接近(图5)，根据反演结果能够轻易识别出真实异常体的位置。总之，本专利中的自约束反演算法求解结果的可靠性明显优于传统反演算法。

tikhonov正则化与分块矩阵lsqr方法计算成本对比

其中：(n×m)表示矩阵阶数，n为数据量，m为模型量；k为lsqr法的迭代次数，一般而言k<<min(m，n)；存储量不包含矩阵求逆所需空间；o表示时间或空间复杂度；—为不涉及。

上表中给出了传统基于tikhonov正则化与本专利反演算法所涉及的计算量和存储空间。其中，本实施例反演算法所需空间复杂度仅为o(2m+2n)，时间复杂度仅为o(2kmn)。由于一般情况下，lsqr算法的迭代次数k远小于数据量n和模型量m，因此，相比tikhonov正则化，本实施例中的反演算法节省了大量存储空间和计算时间。

假设测区网格数据为100×100的规模，模型为100×100×50的三维网格，lsqr算法迭代次数为1000次。则本专利反演算法的计算量约为传统tikhonov正则化的万分之一。若数据以双精度存储，则传统tikhonov正则化需约2tb内存储量，本专利反演算法仅需约100mb内存储量，是前者的二十分之一。可见，如此规模数据的反演计算任务，利用本专利反演算法在普通计算机上就可完成。

然后将本实施例中的方法具体应用在实际矿藏勘探中，具体应用如下：

以某矿区1:2.5万重力数据为例，来展示本专利对找隐伏铅锌矿体的效果。铅锌锑矿是某地区规模最大的多金属矿床，矿体赋存于下地层中，主要受近南北向张(扭)性高角度正断层控制，控矿断裂总体西倾。研究表明，铅锌矿(化)体虽然表现为高密度，但由于矿体赋存于断裂破碎带中，高密度的矿(化)体异常往往被低密度的断裂破碎带异常所淹没，相应的地表观测重力异常值往往表现为负。因此，通过重力三维反演来揭示近南北向的低密度异常位置，可达到间接找矿的目的。共2100个观测点数据，14万多个网格单元，在普通笔记本电脑(内存8gb，双核，处理器：inteli5-6300hq，主频2.30ghz)上计算，共耗时18min，拟合均方误差为0.133mgal，满足要求。

如图6所示，由三维密度成像结果可见，已知的ⅴ号矿体表现为西倾的低密度异常，且该低密度异常在走向上也与ⅴ号矿体一致，说明走向近南北向、倾向向西的低密度异常具有找矿指示作用。在ⅴ号矿体东南侧存在类似的低密度异常，结合其它资料认为该区具有找矿潜力。为此，布设了相应的钻孔进行验证，首钻在深度183～189m见到铅锌工业矿体，后续一系列钻孔验证效果良好，矿石平均品位达11％以上，为区内新发现的ⅹⅴ号矿体，储量达30万吨，实现社会经济效益达10亿元人民币。此外，在ⅴ号矿体和ⅹⅴ号矿体之间具有一个近椭圆型的低密度异常，其它资料表明该低密度异常为低导电性，推测为中酸性岩体，ⅴ号矿体和ⅹⅴ号矿体均位于岩体附近，这与该区岩浆热液成矿作用一致。上述实际数据反演效果说明，本实施例中的三维自约束反演算法具有良好的可靠性、准确性和应用前景。

以某矿集区1：5万重磁数据为例，通过本实施例三维自约束反演算法在得到地下密度和磁性后，按照每个水平网格单元，取其沿深度方向上的密度和磁性，计算两者之间mic值，然后绘制成灰度图，如图7所示，来划分地质体或地质构造。

图7中，黑色实线为根据地质资料确定的错那洞片麻岩穹隆边界(上拆离带)，实线为根据地质资料确定的断了带。由图7可见，错那洞片麻岩穹隆的mic值总体较大且均一，围绕穹隆存在明显的具有环状和线性延伸的条带状、串珠状低mic值异常。其中，环绕穹隆的环状低mic值异常与地质资料确定的穹隆边界位置基本一致；其它局部的呈线性延伸或串珠状分布低mic值异常与已知断裂位置基本一致；此外，还推断了一些隐伏断裂，这些推断断裂的局部位置已得到其他资料的印证。证明本实施例中利用密度和磁性之间的mic来识别达到快速划分地质体和断裂构造是可行和可靠的。

本发明不局限于上述可选的实施方式，任何人在本发明的启示下都可得出其他各种形式的产品。上述具体实施方式不应理解成对本发明的保护范围的限制，本发明的保护范围应当以权利要求书中界定的为准，并且说明书可以用于解释权利要求书。