一种基于多层级方法重磁快速反演方法

本发明属于地球物理及勘探技术领域，尤其涉及重磁勘探技术领域。

背景技术：

为满足大规模、多尺度、大数据量的快速精细三维重磁反演的需要，现有技术中，主要通过数据降维和提高设备计算性能的方式开展快速反演研究。

采用数据降维方法，可以减少核系数矩阵内存消耗，提高计算效率，即使用等效几何构架、小波变换、足印反演和傅里叶变换等技术，改变核系数矩阵的实现方法，在减少核系数矩阵的所需存储空间的同时，能极大地提高反演实施效率，但如小波压缩技术在一定程度上存在损失了信号精度的问题。

利用并行设备加速，在多核心cpu、集群计算机及gpu等平台上，进行大规模并行加速研究，以提高计算设备处理海量矩阵运算的能力，从而有效降低反演解释耗时，但存在硬件门槛过高的问题。

通过约束函数加快反演收敛速度，如通过深度加权函数，优化反演结果的空间分布；通过构造预条件矩阵，改善反演目标函数的收敛性；添加密度约束函数，优化反演结果的上下界。通过多尺度分析方法，如使用八叉树方法，去除某些不必要的系数，使核系数矩阵得以稀疏化，从而大大减少内存需求和加快反演实施；引入l1多尺度小波和全变差技术构造混合正则化反演目标函数，通过添加多尺度小波系数阈值选择函数，提取地质构造边界信息。然而，重力资料多尺度小波处理，对小波母函数的选择依赖比较大，高阶母小波有利于识别场源的细节信息，但可能会放大噪声，继而影响场源识别结果的准确性，多尺度小波系数阈值亦存在类似的问题。

基于以上描述，亟需一种快速反演方法，以解决现有的快速反演方法存在的地球物理多尺度快速反演过程中信号精度低，硬件门槛过高以及场源识别结果的准确性低的问题。

技术实现要素：

本发明目的在于提供一种基于多层级方法重磁快速反演方法,该方法可以提高地球物理多尺度快速反演过程中信号精度，增加场源识别结果准确度，且对硬件要求较低。

为解决上述技术问题，本发明的一种基于多层级方法重磁快速反演方法的具体技术方案如下：

一种基于多层级方法重磁快速反演方法，包括以下步骤：

s1、确定多尺度网格层数n，离散地球物理模型；

s2、确定重磁位场计算公式；

s3、设置虚拟线线观测网格，获得核系数矩阵gq,m；

s4、设置虚拟面面观测网格，获得核系数矩阵

s5、构造基于多尺度haar小波的插值算子和限制算子；

s6、通过插值算子和限制算子将大尺度观测数据向下采样，并在不同层间逐一实施预光滑反演；

s7、实施校正反演。

作为优选，步骤s1采用点元法使用平行的轴截平面将地球物理模型分割成多个直立六面体。

作为优选，步骤s2中重磁位场计算公式为矩阵形式：

d＝gm

式中，d为观测数据，均为长度nd的向量；m为长度为nm物性参数向量；g为相应的正演算子，为nd×nm大小的矩阵；nm和nd分别为剖分网格数和观测数据点个数；nm＝nx×ny×nz，nd＝nx×ny，其中nx、ny和nz分别为模型沿三轴向的剖分数。

作为优选，所述观测数据d为gz、gxx、gxy、gxz、gyy、gyz或gzz。

作为优选，步骤s3中gq,m为第q条测线及所对应的近地表第m列网格的核系数矩阵，gq,m与对应的物性单元mq,m的乘积gq,m*mq,m为：

式中，和f为傅里叶变换对；作为替代核系数矩阵而需要存储的部分为：

作为优选，步骤s4中为第t层观测点及所对应的第n层物性网格观测方案的核系数矩阵，与对应的物性单元的乘积为：

作为替代核系数矩阵而需要存储的部分

作为优选，相邻两层核系数矩阵的转换为：

gⁱ⁺¹＝rⁱgⁱpⁱ(14)

式中，r为限制算子，p为传递算子或插值算子。

作为优选，步骤s5中构造基于多尺度haar小波的插值算子和限制算子的方法为：

选择haar小波作为r^t和p，则：

式中，j,k＝1,2，w和w^t分别为haar小波矩阵正变换和逆变换，上标i表示预处理的层号，i＝0对应于最精细网格，i＝n为最粗网格。

作为优选，通过步骤s6中预光滑反演获得的近似结合数据泛函和模型泛函，引入深度加权函数，采用正则化方法得：

式中：

为小波域参考模型；

对角矩阵diag(wz)＝(z+z0)^-κ/2为深度加权矩阵,z为物性网格中心埋深,κ为场源相关的衰减系数,z0为与物性网格大小和观测高度相关的标量；

wm为粗糙度矩阵，常为1阶，2阶差分算子；

wd为稀疏对角矩阵：

式中，为观测数据的标准差，ε为极小量。

作为优选，步骤s7中校正反演：

有益效果：

本发明充分利用重磁位场正演核矩阵为特殊结构矩阵和多尺度haar小波的特点，将大尺度重磁正演核矩阵转换至小波域，实现基于toeplitz的小波多尺度域的快速正演，且由于在多层级反演中将下一层反演的结果作为参考模型，从而可以实现快速反演。除此之外，本发明所提供的技术方案还具有以下优点：

1.使用haar多尺度小波，保证不同层的核矩阵仍然具有特殊结构，构造基于多尺度haar小波的插值算子和限制算子。

2.通过插值算子和限制算子将大尺度观测数据向下采样，且将下一层反演的结果作为参考模型，不同层间逐一实现快速反演。

3.进行校正反演，解决向下采样的数据可能丢失精度的问题。

4.将大尺度、多面积和大数据量数据进行逐级分块处理，大大降低对硬件的要求，从而使传统计算机亦能处理大尺度、多面积和大数据量数据。

附图说明

图1为本实施例提供的基于多层级方法重磁快速反演优选方法流程图；

图2为本实施例提供的虚拟线线观测方案示意图；

图3为本实施例提供的虚拟面面观测方案示意图；

图4为本实施例提供的正演模型示意图；

图5为第一层快速反演结果结果；

图6为第二层快速反演结果结果；

图7为第三层快速反演结果结果。

具体实施方式

为了更好地了解本发明的方法及步骤，下面结合附图，对本发明一种基于多层级方法重磁快速反演方法做进一步详细的描述。

本发明基于多层级方法重磁快速反演方法理论及推导如下：

在笛卡尔坐标系下，存在一剩余质量为m，体积为v，剩余密度为ρ的三度体。采用点元法使用一系列平行的轴截平面将地球物理模型分割成大量的直立六面体。任意观测点p对地下空间内任意长方体[ξ1→ξ2,η1→η2,ζ1→ζ2]的全重力梯度张量无解析奇点正演计算公式，以重力gz为例：

式中，

μijk＝(-1)^i+j+k，

xi＝x-ξi，yj＝y-ηj，zk＝z-ζk。

对于所有观测点的异常响应而言，可根据叠加原理逐一计算地下空间内各直立六面体对相应观测点的异常响应。因而，可将重力及重力张量正演计算写成矩阵形式：

d＝gm(1)

式中，d为观测数据，可具体为gz、gxx、gxy、gxz、gyy、gyz和gzz等，均为长度nd的向量；m为长度为nm物性参数向量，即物性网格的物性参数，如密度值、磁化强度或磁化率值等；g为相应的正演算子，为nd×nm大小的矩阵；nm和nd分别为剖分网格数和观测数据点个数。nm＝nx×ny×nz，nd＝nx×ny，其中nx、ny和nz分别为模型沿三轴向的剖分数。

对于传统重磁勘探而言，于直角坐标系下，以重力场正演计算为例，当正演核函数计算公式确定时，核函数仅与观测点和网格单元角点两者空间相对位置有关。这里，引入三轴向网格剖分数标记<*,*,*>，其为关于nx,ny的函数。以沿笛卡尔坐标系三轴向剖分数分别为<l,m,n>和<p,q,t>分别描述第j个物性网格q和第i个观测点p，则：

式中，<1,1,1>为计算原点。

由式(2)可知，利用核函数内存的对称性，任一观测点于任一网格的核函数的计算，可通过换算至计算原点计算，这极大的简便了某些具有特定对称特性分量，但并不是所有分量均具有该特性。

如图2所示，设置虚拟线线观测方案。这里所谓的线线观测方案为单一勘探测线所对应地球物理模型中单一列物性网格的观测方案。从而，基于平移等效性的等效几何构架计算公式为：

式中，1≤l≤nx，1≤m≤nx，1≤p≤ny，1≤q≤nx，1≤n≤nz，t≥1。

对于常规勘探方法而言，观测数据d是二维的，即t＝1；对于等维反演而言，观测数据d是三维的，t≥1。采用公式(3)将极大地简化不同分量的正演计算。

式中，i'＝(q-1)nx+(t-1)nxny，j'＝(m-1)ny+(n-1)nxny。

在式(3)中，p和l分别以δp和δl(且δp＝δl)为步长沿x轴自1至nx开始移动，由式(4)可知：

根据式(5)，对于任意一q和m，都能在式(4)左侧矩阵中找到各元素相等的对角线，如q＝1且m＝1时，δp自1至nx沿x轴开始移动，即有：

按照相类似的方法，可以构建其它对角线，即其任意对角线的核系数幅值都一致。因而，式(4)左侧矩阵可以由2nx-1个元素表达，可将其重写为：

式中，gq,m对应着第q条测线及所对应的近地表第m列网格的核系数矩阵，a-p与al为式(4)下/上三角矩阵中对角线的第一个元素。

根据高等数学知识，将gq,m与对应的物性单元mq,m的乘积gq,m*mq,m可写为：

式中，和f为傅里叶变换对；作为替代核系数矩阵而需要存储的部分为：

这里将式(8)定义为函数t(gq,m,mq,m)，对于多个不同埋深物性网格列所对应的mq,m矩阵(即gq,m,n)与单一向量v的一一乘积之和而言，通过类似式(4)～(8)推导，可以得出：

∑n＝1gq,m,nv＝(∑n＝1gq,m,n)v(9)

上式主要应用于在位场反演中g^td的计算，可进一步极大地加快反演的实施。

如图3所示，设置虚拟面面观测方案。这里所谓的面面观测方案为由多条勘探测线组成的观测面与所对应的地球物理模型中的单层物性网格的观测方案，该方案在现实地球物理勘探中并不存在，仅为推导快速正演计算方法而提出的。

通过式(7)构建其它测线与各列物性网格的关系，并将其扩展至第t层观测点及所对应的第n层物性网格观测方案的核系数矩阵

同样可以采用式(8)类似的推导过程可导出的计算公式：

作为替代核系数矩阵而需要存储的部分

对于传统位场勘探而言，线线观测方案、面面观测方案和等维反演方案等三类观测方案，其中的核系数矩阵分别为gq,m、和因gq,m为toeplitz矩阵，故而gq,m和也为toeplitz矩阵。

针对在处理大尺度、多面积、大数据量的重磁梯度及其张量数据时，点元法难于计算或存储核矩阵这一瓶颈问题,借鉴基于多尺度金字塔和代数多重网格等方法多尺度分层思想，将大尺度、多面积和大数据量数据进行逐级分块处理，从而使传统计算机亦能处理大尺度、多面积和大数据量数据。则式(1)可以改写为：

gⁱmⁱ＝dⁱ,0≤i≤n(13)

式中，上标i表示预处理的层号，i＝0对应于最精细网格，即g⁰＝g，i＝n为最粗网格。

相邻两层核系数矩阵的转换，可以写为：

gⁱ⁺¹＝rⁱgⁱpⁱ(14)

式中，r为限制算子，p为传递算子，或插值算子，即由细网格至粗网格间核函数的传递，反之亦然。

当选择haar小波作为r^t和p时：

将式(13)改写到频率域，则：

式中，w和w^t分别为haar小波矩阵正变换和逆变换,和

这里，利用w1和w2将式(15)写成块矩阵乘法形式：

其中，和

对于j,k＝1,2，可以被证明为典型的toeplitz矩阵和bttb矩阵的组合,按虚拟线线观测方案将gⁱ展开为矩阵形式：

这里矩阵的toeplitz向量为：

进而,有：

其中，表示矩阵e和f的张量积，vec(e)为对矩阵e的向量化运算。

由于计算式(18)需要将使用原有核矩阵因而在内存消耗上并没有改进。为此，将式(18)进行改写获得的toeplitz向量

式中，为wj小波基。

为求解式(16)中的将其改写为：

由于是未知参数,为此，通过预光滑反演获得的近似则可将式(20)改写为：

对于最粗网格而言，为避免趋肤效应，尽可能的保存反演结果深度分辨率，为此，结合数据泛函和模型泛函，引入深度加权函数，采用正则化方法求解式(21)的最小化问题：

式中，为小波域参考模型。对角矩阵diag(wz)＝(z+z0)^-κ/2为深度加权矩阵,z为物性网格中心埋深,κ为场源相关的衰减系数,z0为与物性网格大小和观测高度相关的标量.wm为粗糙度矩阵，常为1阶，2阶差分算子,稀疏对角矩阵wd：

这里，为观测数据的标准差，ε为极小量，以避免对小幅值观测数据施加过大的权重。

为求解(16)中根据式(21),有：

由于在多层级反演过程中，通过插值算子和限制算子，实现观测数据的向下插值和向上插值，这可能导致数据精度的丢失，为避免因数据精度的丢失导致反演结果错误，引入校正反演：

为结合实施例对本发明进行说明，建立理论模型，如图1所示，本发明基于多层级方法重磁快速反演方法具体实施流程如下：

1)确定多层级方法层数，如n＝3，按nx＝ny可被2³整除划分地球物理模型；

2)确定重磁位场计算公式，这里以重力gz为例；

3)以虚拟线线观测方案，确定核矩阵需要使用的toeplitz向量t⁰；

4)于不同层构造toeplitz向量tⁱ；

5)于不同层构造多尺度小波算子

6)于小波域，实现不同块矩阵的快速正演

7)进入多层级反演，多层级反演流程如下：

第1层至第3层的快速反演结果如图5至图7所示。

可以理解，本领域技术人员知悉的，在不脱离本发明的精神和范围的情况下，可以对该方法进行各种改变或等效替换。另外，在本发明的教导下，可以对该方法进行修改以适应具体的情况而不会脱离本发明的精神和范围。因此，本发明不受此处所公开的具体实施例的限制，所有落入本申请的权利要求范围内的实施例都属于本发明所保护的范围内。