CN108227024A - 一种采用全张量磁梯度数据反演地下磁化率的方法及系统 - Google Patents

Abstract

本发明涉及一种采用全张量磁梯度数据反演地下磁化率的方法及系统，该方法包括：获取全张量磁梯度数据；根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型；基于改进的共轭梯度算法对所述稀疏正则化模型进行迭代求解，反演得到磁化率。还涉及一种系统，该系统包括：数据获取模块、模型建立模块、磁化率反演模块。通过本发明设置一些先验约束条件，建立稀疏正则化模型，这样可以考虑物理参数的有界性约束，大大提高了反演磁化率的精度。

Description

一种采用全张量磁梯度数据反演地下磁化率的方法及系统

技术领域

本发明属于磁法勘探和地球物理勘探领域，尤其涉及一种采用全张量磁梯度数据反演地下介质磁化率的方法及系统。

背景技术

磁法勘探地球物理勘探领域中的重要组成部分，在矿物、石油勘探中有着广泛的应用。例如，在矿物勘探中，磁方法可以用于决定地下介质体参数；在石油勘探中，磁方法可以用来映射地下的沉积特性和缺点，用来控制沉积盆地的沉积背景。

磁法勘探是通过观测和分析由岩石、矿石等磁性差异所引起的磁异常，进而研究地质构造和矿产资源或其他探测对象分布规律的一种地球物理方法。根据磁测数据反演出地下物理参数，例如磁化率和磁化强度，是磁法勘探中的主要科学问题。实际应用中，能够反演目标地理区域的磁化率，得到该目标地理区域的磁化率模型，通过该磁化率模型能够反映该目标地理区域内各个地理位置的磁化率分布，进而将磁化率符合预设要求的某地理位置选择出来，作为进一步的研究对象。相关技术中通常基于总磁场强度(TMI)数据反演求解磁化率，而磁梯度张量数据可以有效抑制甚至消除区域性磁场干扰和时变磁场影响，跟传统磁测方法相比具有较大优势。且传统磁化率反演技术只考虑统计正则化和基于变换的滤波方法来重构物理参数，没有考虑模型的系属性以及有界性约束，从而导致反演得到的磁化率不够精确。

发明内容

本发明所要解决的技术问题是：传统基于TMI数据的磁化率反演方法精度较低；现有磁化率反演技术中只是考虑统计正则化和基于变换的滤波方法来重构物理参数，但没有考虑关于物理参数的有界性约束，从而导致反演得到的磁化率不够精确。

为解决上面的技术问题，本发明提供了一种采用全张量磁梯度数据反演地下磁化率的方法，该方法包括：

S1，获取全张量磁梯度数据；

S2，根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型；

S3，基于改进的共轭梯度算法对所述稀疏正则化模型进行迭代求解，反演得到磁化率。

本发明的有益效果：通过上述的方法，设置一些先验约束条件，建立稀疏正则化模型，这样可以考虑物理参数的有界性约束，大大提高了反演磁化率的精度。

进一步地，所述步骤S2中建立所述稀疏正则化模型的具体公式为：

其中，m表示磁化率向量，L表示离散化紧算子，d表示所述全张量磁梯度数据，表示l₂范数，·||·||₁表示l₁范数，并且α，ξ＞0是指正则化参数。上述进一步地有益效果：通过增加了对Tikhonov正则化模型的稀疏约束，对地下介质进行稀疏假设，使得反演结果更加精确。

进一步地，步骤S3包括：

步骤S31，预设初始磁化率m₀，预设迭代停止条件阈值e＞0，预设正则化参数α和ξ，预设迭代次数k＝0，计算初始梯度值▽J^α,ξ(m_k)，并判断||▽J^α,ξ(m_k)||是否小于或等于预设迭代停止条件阈值e，若是，则执行步骤S36，若否，则执行步骤S32，g表示所述正则化模型的梯度；

步骤S32，判断预设迭代次数k是否大于0，当k＞0，则确定一个用于计算梯度时所需的步长β_k-1＝||▽J^α,ξ(m_k)||²/||▽J^α,ξ(m_k-1)||²，否则，β_k-1＝0；

步骤S33，根据所述β_k-1和所述初始梯度值▽J^α,ξ(m_k)，计算搜索方向h_k，其中所述h_k＝-▽J^α,ξ(m_k)+β_k-1h_k-1；

步骤S34，计算搜索步长τ_k，其中所述搜索步长 其中，τ_k表示所述迭代步长β₁和β₂均表示正实数，H表示所述海森矩阵，g表示所述梯度；

步骤S35，更新磁化率m_k，计算磁化率m_k+1＝m_k+τ_kh_k，令k＝k+1，重复执行步骤S1，直到所述||▽J^α,ξ(m_k)||小于预设迭代停止条件阈值e；

步骤S36，停止计算，并输出磁化率m_k。

进一步地，判断预设迭代次数k是否大于0，当k＞0，则计算β_k-1＝||P^1/2▽J^α,ξ(m_k)||²/||P^1/2▽J^α,ξ(m_k-1)||²，其中预收敛条件矩阵P；

则所述步骤S33为：根据所述β_k-1和所述初始梯度值▽J^α,ξ(m_k)，计算搜索方向h_k，其中所述h_k＝-▽J^α,ξ(m_k)+β_k-1h_k-1；

其中P为预收敛条件矩阵，这里选取L^TL+αI作为预条件矩阵。

进一步地，所述步骤S32中还包括：当时，将所述β_k替换为

其中β^hybrid表示共轭方向混合参数变量，β^HS表示HS型共轭方向参数变量，β^DY表示DY型共轭方向参数变量，u(·)表示单位阶跃函数。

上述进一步地有益效果：通过使用混合模型参数，可以增大共轭梯度算法的步长，提高了算法的精度。

进一步地，所述步骤S35还包括：

计算所述磁化率m_k+1＝P_Π(m_k+τ_kh_k)，其中，P_Π表示凸集投影算子。

本发明还涉及一种采用全张量磁梯度数据反演地下磁化率的系统，该系统包括：数据获取模块、模型建立模块、磁化率反演模块；

所述数据获取模块，用于获取全张量磁梯度数据：

所述模型建立模块，用于根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型；

所述磁化率反演模块，用于基于改进的共轭梯度算法对所述稀疏正则化模型进行迭代求解，反演得到磁化率。

本发明的有益效果：通过上述的系统，设置一些先验约束条件，建立Tikhonov正则化模型，这样可以考虑物理参数的有界性约束，大大提高了反演磁化率的精度。

进一步地，所述模型建立模块，用于根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型，其中建立所述稀疏正则化模型的具体公式为：

其中，m表示磁化率向量，L表示离散化紧算子，d表示所述全张量磁梯度数据，表示l₂范数，||·||₁表示l₁范数，并且α，ξ＞0是指正则化参数。

还涉及一种计算机设备，该计算机设备包括：处理器、存储器及存储在所述存储器上并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如上述任一项所述方法的步骤。

附图说明

图1为本发明的一种采用全张量磁梯度数据反演地下磁化率的方法的流程图；

图2为本发明的一种采用全张量磁梯度数据反演地下磁化率的系统的结构示意图。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

如图1所示，本发明实施例1提供的是一种采用全张量磁梯度数据反演地下磁化率的方法，该方法包括：

S1，获取全张量磁梯度数据；

S2，根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型；

S3，基于改进的共轭梯度算法对所述稀疏正则化模型进行迭代求解，反演得到磁化率。

通过上述的方法，设置一些先验约束条件，建立稀疏正则化模型，这样可以考虑物理参数的有界性约束，大大提高了反演磁化率的精度。

建立稀疏正则化模型，下面是对建立该模型的方式：

假设地下模型I被网格化为N个单元，并且在每个单元中磁化率是常数，以下线性代数方程：

Lm＝d, (1)

其中，L是M×N矩阵，L表示离散化紧算子；m是磁化率向量，长度为N；d是长度为M的全张量磁梯度数据。

由于离散方程(1)是病态的。在现实中，全张量磁梯度数据d可能带有噪声，例如，d＝d_true+n，其中，d_true表示真实的磁张量场，n表示噪声。解决方程(1)常常归结为正则化问题

其中，m表示磁化率向量，L表示离散化紧算子，d表示所述全张量磁梯度数据，表示l₂范数，||·||₁表示l₁范数，并且α，ξ＞0是指正则化参数。

关于约束项的梯度可以推导如下，

g_Ω(m)＝m (3)

最小方差的梯度可以计算得到如下形式，

g_ρ(m)＝L^TLm-L^Td (4)

地下介质异常体的非光滑信息f(m)＝||m||₁，为了计算方便，我们将||m||₁近似为其中l是向量m的长度。f(m)＝||m||₁的梯度可以计算得到如下形式，

因此，目标函数Jα(m)＝ρ²(Lm,d)+αΩ(m)+ξg_f(m)的梯度可以写成如下形式：

g(m)＝▽J^α,ξ(m)＝g_ρ(m)+αg_Ω(m)+ξg_f(m) (5)

有两种方法可以用来解决目标函数J^α,ξ(m)＝ρ²(Lm,d)+αΩ(m)+ξf(m)(2)的最小化问题：梯度下降算法和牛顿算法。牛顿算法是最快速收敛的，但是需要计算Hessian矩阵并需要大量的存储空间；对于3维的磁梯度张量这样的大规模计算问题，采用梯度下降算法是必要的。

基于改进的共轭梯度算法对所述Tikhonov正则化模型进行迭代求解，反演得到磁化率，其过程如下：

步骤S31，预设初始磁化率m₀，预设迭代停止条件阈值e＞0，预设正则化参数α和ξ，预设迭代次数k＝0，计算初始梯度值▽J^α,ξ(m_k)，并判断||▽J^α,ξ(m_k)||是否小于或等于预设迭代停止条件阈值e，若是，则执行步骤S36，若否，则执行步骤S32，g表示所述正则化模型的梯度；

步骤S32，判断预设迭代次数k是否大于0，当k＞0，则确定一个用于计算梯度时所需的步长β_k-1＝||▽J^α,ξ(m_k)||2/||▽J^α,ξ(m_k-1)||²，否则，β_k-1＝0；

步骤S33，根据所述β_k-1和所述初始梯度值▽J^α,ξ(m_k)，计算搜索方向h_k，其中所述h_k＝-▽J^α,ξ(m_k)+β_k-1h_k-1；

步骤S34，计算搜索步长τ_k，其中所述搜索步长 其中，τ_k表示所述迭代步长β₁和β₂均表示正实数，H表示所述海森矩阵，g表示所述梯度；

步骤S35，更新磁化率m_k，计算磁化率m_k+1＝m_k+τ_kh_k，令k＝k+1，重复执行步骤S1，直到所述||▽J^α,ξ(m_k)||小于预设迭代停止条件阈值e；

步骤S36，停止计算，并输出磁化率m_k。

在另一实施例中考虑预条件矩阵P，由此还可以推导出如下公式：

步骤S31，计算初始梯度值▽J^α,ξ(m_k)，并判断||▽J^α,ξ(m_k)||是否小于或等于预设迭代停止条件阈值e，若是，则执行步骤S36，若否，则执行步骤S32，其中所述初始梯度值▽J^α,ξ(m_k)为▽J^α,ξ(m_k)＝g_ρ(m_k)+αg_Ω(m_k)+ξg_f(m_k)，所述α和ξ分别为正则化参数，g表示所述正则化模型的梯度；

步骤S32，判断预设迭代次数k是否大于0，当k＞0，则计算β_k-1＝||P^1/2▽J^α,ξ(m_k)||²/||P^1/2▽J^α,ξ(m_k-1)||²，其中预收敛条件矩阵P，这里选取L^TL+α I作为预条件矩阵。；

步骤S33：根据所述β_k-1和所述初始梯度值▽J^α,ξ(m_k)，计算搜索方向h_k，其中所述h_k＝-P▽J^α,ξ(m_k)+β_k-1h_k-1，其中P为预收敛条件矩阵；

步骤S34，根据所述搜索方向h_k和非精确线搜索方法，计算搜索步长τ_k，其中所述搜索步长其中，τ_k表示所述迭代步长，k表示迭代次数，β₁和β₂均表示正实数，H表示所述海森矩阵，g表示所述梯度；

步骤S35，更换新的所述正则化参数α和ξ，计算磁化率m_k+1＝m_k+τ_kh_k，其中令k＝k+1，重复执行步骤S1，直到所述||▽J^α,ξ(m_k)||等于预设迭代停止条件阈值e；

步骤S36，停止计算，并输出磁化率m_k。

其中我们将上述算法称之为预条件共轭梯度算法，其中P是预收敛条件矩阵。

下面是一些具体的数据输入的例子：输入：给出初始磁化率向量m₀，预条件矩阵P，输入磁梯度张量数据d，设置α＝0.001，ξ＝0.1，ε＞0，并设置迭代索引k＝0；

Step1：计算初始梯度值▽J^α,ξ(m_k)＝g_ρ(m_k)+αg_Ω(m_k)+ξg_f(m_k)，如果||▽J^α(m_k)||≤e，转到Step6；

Step2：如果k＞0，计算β_k-1＝||P^1/2▽J^α,ξ(m_k)²/||P^1/2▽J^α,ξ(m_k-1)||²；

Step3：计算搜索方向

Step4：由计算τ_k；

Step5：计算m_k+1＝m_k+τ_kh_k，并设置k:＝k+1，转到step1。

Step6：计算结束，输出磁化率向量m_k

可以看到，将||m||₁近似为后，目标函数J^α,ξ(m)＝ρ²(Lm,d)+αΩ(m)+ξf(m)是连续且二阶可导的，并且▽J^α,ξ(m_k)满足Lipstichz条件。另外，由预条件共轭梯度算法产生的迭代点是有界的，例如，对于任意的m∈F＝{m:J^α,ξ(m)≤J^α,ξ(m₀)}，存在δ＜∞使得||m_k||＜δ。因此，可以证明进而可以得到由此可以得到上述算法是收敛的。

可选地，在另一实施例中由于上述共轭梯度算法的步长可能会很小，以至于引起欠定反演问题收敛所需的迭代次数增加。为了克服这个缺陷并加速共轭梯度算法的收敛速度，对此进行了一些改进。在用于共轭梯度算法的公式基于选择合适的参数βk的混合模型。标记和其中，上标HS和DY分别表示Hestenes和Stiefel公式以及Dai和Yuan公式。对于前面所述的预条件技术，这两个参数可以重新写成如下形式：

其中，P是预条件算子可以如前述方法进行选择，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型共轭方向参数变量，β^DY表示DY型共轭方向参数变量。

由此，我们的改进的共轭梯度算法公式可以基于如下公式：

选择上述公式的原因在于HS方法即使在简单的二维最小化问题上也会出现发散现象。但是HS方法具有共轭性：对于任何线性搜索都适用。DY方法对于普遍的线性搜索技术是全局收敛的。本实施例中有一个显著的优势，其可以避免小步长的出现；尤其对于病态问题，小步长可能会引起误差累计并导致迭代步数增加才能收敛。

可选地，在另一实施例中计算所述磁化率m_k+1＝P_Π(m_k+τ_kh_k)，其中，P_Π表示凸集投影算子。

如图2所示，在本发明实施例中还涉及一种采用全张量磁梯度数据反演地下磁化率的系统，该系统包括：数据获取模块、模型建立模块、磁化率反演模块；

所述数据获取模块，用于获取全张量磁梯度数据：

所述模型建立模块，用于根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型；

所述磁化率反演模块，用于基于改进的共轭梯度算法对所述稀疏正则化模型进行迭代求解，反演得到磁化率。

需要说明的是，本实施例中的系统适用于上述提及到的所有的方法实施例。

可选地，在另一实施例中所述模型建立模块，用于根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型，其中建立所述稀疏正则化模型的具体公式为：

其中，m表示磁化率向量，L表示离散化紧算子，d表示所述全张量磁梯度数据，表示l₂范数，||·||₁表示l₁范数，并且α，ξ＞0是指正则化参数。

本发明实施例中还涉及一种计算机设备，该计算机设备包括：处理器、存储器及存储在所述存储器上并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如权利要求1-6中任一项所述方法的步骤。

在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种采用全张量磁梯度数据反演地下磁化率的方法，其特征在于,该方法包括：

S1，获取全张量磁梯度数据；

S2，根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型；

S3，基于改进的共轭梯度算法对所述稀疏正则化模型进行迭代求解，反演得到磁化率。

2.根据权利要求1所述的方法，其特征在于，所述步骤S2中建立所述稀疏正则化模型的具体公式为：

其中，m表示磁化率向量，L表示离散化紧算子，d表示所述全张量磁梯度数据，表示l₂范数，||·||₁表示l₁范数，并且α，ξ＞0是指正则化参数。

3.根据权利要求2所述的方法，其特征在于，步骤S3包括：

步骤S31，预设初始磁化率m₀，预设迭代停止条件阈值e＞0，预设正则化参数α和ξ，预设迭代次数k＝0，计算初始梯度值▽J^α,ξ(m_k)，并判断||▽J^α,ξ(m_k)||是否小于或等于预设迭代停止条件阈值e，若是，则执行步骤S36，若否，则执行步骤S32，g表示所述正则化模型的梯度；

步骤S32，判断预设迭代次数k是否大于0，当k＞0，则确定一个用于计算梯度时所需的步长β_k-1＝||▽J^α,ξ(m_k)||²/||▽J^α,ξ(m_k-1)||²，否则，β_k-1＝0；

步骤S33，根据所述β_k-1和所述初始梯度值▽J^α,ξ(m_k)，计算搜索方向h_k，其中所述h_k＝-▽J^α,ξ(m_k)+β_k-1h_k-1；

步骤S34，计算搜索步长τ_k，其中所述搜索步长 其中，τ_k表示所述迭代步长β₁和β₂均表示正实数，H表示所述海森矩阵，g表示所述梯度；

步骤S35，更新磁化率m_k，计算磁化率m_k+1＝m_k+τ_kh_k，令k＝k+1，重复执行步骤S1，直到所述||▽J^α,ξ(m_k)||小于预设迭代停止条件阈值e；

步骤S36，停止计算，并输出磁化率m_k。

4.根据权利要求3所述的方法，其特征在于，所述步骤S32还包括：

判断预设迭代次数k是否大于0，当k＞0，则计算β_k-1＝||P^1/2▽J^α,ξ(m_k)||²/||P^1/2▽J^α,ξ(m_k-1)||²，其中预收敛条件矩阵P；

则所述步骤S33为：根据所述β_k-1和所述初始梯度值▽J^α,ξ(m_k)，计算搜索方向h_k，其中所述h_k＝-▽J^α,ξ(m_k)+β_k-1h_k-1；

其中P为预收敛条件矩阵，这里选取L^TL+αI作为预条件矩阵。

5.根据权利要求4所述的方法，其特征在于，所述步骤S32中还包括：当时，将所述β_k替换为

其中β^hybrid表示共轭方向混合参数变量，β^HS表示HS型共轭方向参数变量，β^DY表示DY型共轭方向参数变量，u(·)表示单位阶跃函数。

6.根据权利要求3-5任一所述的方法，其特征在于，所述步骤S35还包括：

计算所述磁化率m_k+1＝P_Π(m_k+τ_kh_k)，其中，P_Π表示凸集投影算子。

7.一种采用全张量磁梯度数据反演地下磁化率的系统，其特征在于，该系统包括：数据获取模块、模型建立模块、磁化率反演模块；

所述数据获取模块，用于获取全张量磁梯度数据：

所述模型建立模块，用于根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型；

所述磁化率反演模块，用于基于改进的共轭梯度算法对所述稀疏正则化模型进行迭代求解，反演得到磁化率。

8.根据权利要求7所述的系统，其特征在于，所述模型建立模块，用于根据预设的先验约束条件和所述全张量磁梯度数据，建立稀疏正则化模型，其中建立所述稀疏正则化模型的具体公式为：

其中，m表示磁化率向量，L表示离散化紧算子，d表示所述全张量磁梯度数据，表示l₂范数，||·||₁表示l₁范数，并且α，ξ＞0是指正则化参数。