CN115128700A - 一种基于gramian约束的重磁三维联合反演方法及系统 - Google Patents

Abstract

本发明涉及联合反演领域，提供一种基于gramian约束的重磁三维联合反演方法及系统，包括：S1：将初始密度模型作为当前密度模型，将初始磁化率模型作为当前磁化率模型；S2：构建联合反演目标函数，计算联合约束项、重力、磁法正则化项和各项对应正则化参数；S3：获得更新后的密度模型、更新后的磁化率模型和数据拟合差；S4：若数据拟合差小于目标拟合差，则输出当前的联合反演目标函数；否则返回步骤S2。本发明约束密度与磁化率模型物性参数或其梯度为线性关系；使用LSQR方法求解简化后的联合反演目标方程，减少了反演中的计算量；提高了联合反演的稳定性；提高了联合反演结果的可靠程度。

Description

一种基于gramian约束的重磁三维联合反演方法及系统

技术领域

本发明涉及联合反演领域，尤其涉及一种基于gramian约束的重磁三维联合反演方法及系统。

背景技术

重力和磁法是地球物理勘探中的重要方法，被广泛应用于金属矿产勘探、油气资源开发和地质构造研究等领域。由于观测数据的有限性与噪声的影响，重磁数据的三维反演往往存在较强非唯一性和不稳定性，严重影响地质解释的可靠程度。联合反演作为抑制非唯一性的有效手段，能够得到更可靠的反演模型，提高地质解释的精度。联合反演的优势主要由两方面决定，一是不同的地球物理方法探测不同的物性；二是不同数据包含的噪声存在差异；这些都可以在联合反演中互补。

目前的联合反演方法根据耦合方式可分为两大类，第一类为岩石物理耦合方法，该方法假设地下介质的不同物性之间满足已知的数学关系，这种岩石物理关系可以来自于理论关系、经验公式、岩样统计等方法。岩石物理耦合通过已知的数学关系将两种或多种物性联合起来，该方法可以在反演模型之间施加一个强力约束，但也过于依赖数学关系的正确性，当数学关系与实际地质情况不相符时，容易得到与真实情况相差较大的反演结果，且实际地质情况往往是复杂的，介质不同物性之间的关系难以用一个数学公式表示，因此，该方法有一定的局限性。

第二类为结构耦合方法，该方法假设地下介质表现的不同物性模型在结构上是相互联系的。交叉梯度约束是结构耦合中最典型的方法，该算子定义为两种模型参数的空间梯度的叉乘，当两种模型在任意空间位置的梯度方向相同或相反时，交叉梯度算子为0，随后基于交叉梯度约束的联合反演在地球物理中得到了广泛的应用。

基于Gramian约束的联合反演是一种更加广义的联合反演方法，Gramian算子是一个可以将两种或多种物性模型向量联系起来的行列式，将其加入目标函数中进行最小化可实现不同地球物理数据的联合反演。由于联合反演需要同时更新两个模型参数，所以计算所消耗的时间与内存是非常巨大的，因此开发高效可靠的联合反演算法是有必要的。

上述内容仅用于辅助理解本发明的技术方案，并不代表承认上述内容是现有技术。

发明内容

为解决上述技术问题，本发明提供一种基于gramian约束的重磁三维联合反演方法，包括：

S1：构建初始密度模型和初始磁化率模型，将初始密度模型作为当前密度模型，将初始磁化率模型作为当前磁化率模型；

S2：计算获得重力正则化参数和磁法正则化参数；引入联合约束项，通过当前密度模型、当前磁化率模型、联合约束项、重力正则化参数和磁法正则化参数构建联合反演目标函数；

S3：对联合反演目标函数进行优化更新，获得更新后的密度模型、更新后的磁化率模型和数据拟合差；

S4：若数据拟合差小于目标拟合差，则输出当前的联合反演目标函数，结束联合反演；否则将更新后的密度模型作为当前密度模型，将更新后的磁化率模型作为当前磁化率模型，返回步骤S2。

优选的，步骤S2具体为：

S21：计算获得重力加权矩阵W_m1，计算公式为：

W_m1＝diag(A₁ ^TA₁)^1/2

式中A₁为重力数据的线性正演算子；

通过重力加权矩阵计算获得重力正则化参数a₁，计算公式为：

式中，b₁和c₁为重力系数参数，x为随机列向量；

S22：计算获得磁法加权矩阵W_m2，计算公式为：

W_m2＝diag(A₂ ^TA₂)^1/2

式中A₂为磁法数据的线性正演算子；

通过磁法加权矩阵计算获得磁法正则化参数，计算公式为：

式中n为迭代次数，b₂和c₂为磁法系数参数，x为随机列向量；

S23：引入的联合约束项φ_G(m^n-1)表示为：

式中m^n-1为联合反演模型，

为当前密度模型，

为当前磁化率模型；G为Gramian算子；

S24：联合反演目标函数φ(m^n-1)的表达式为：

φ(m^n-1)＝φ_d(m^n-1)+a₁φ_m(m^n-1)+a₂φ_G(m^n-1)

式中φ_d(m^n-1)为数据拟合项，φ_m(m^n-1)为深度加权项。

优选的，步骤S3具体为：

S31：计算Gramian函数值对联合反演目标函数的梯度；

S32：通过高斯-牛顿法最小化联合反演目标函数，获得联合反演目标函数更新方向；

S33：通过联合反演目标函数更新方向计算获得更新后的联合反演目标函数，通过更新后的联合反演目标函数计算获得更新后的密度模型、更新后的磁化率模型和数据拟合差。

优选的，Gramian函数值对联合反演目标函数的梯度

的表达式为：

式中L_x,L_y,L_z分别是x、y、z方向的差分算子。

优选的，联合反演目标函数更新方向δm的表达式为：

式中W_d为数据加权矩阵，W_d＝diag(AA^T)^-1/2；W_m为深度加权矩阵，W_m＝diag(A^TA)^1/2；m为联合反演模型；α和γ分别为重力正则化参数和磁法正则化参数；G为Gramian算子；

A₁为重力数据的线性正演算子，A₂为磁法数据的线性正演算子。

优选的，更新后的联合反演目标函数mⁿ表示为：

式中n为迭代次数，δm为联合反演目标函数更新方向，I_d为最小化数据拟合差的步长，I_G为最小化Gramian算子的步长；

更新后的密度模型和更新后的磁化率模型的计算公式为：

式中

为更新后的密度模型，

为更新后的磁化率模型。

优选的，数据拟合差misfit的计算公式为：

式中W_d为数据加权矩阵，d_pre为联合预测数据，d_obs为联合观测数据。

一种基于gramian约束的重磁三维联合反演系统，包括：

模型构建模块，用于构建初始密度模型和初始磁化率模型，将初始密度模型作为当前密度模型，将初始磁化率模型作为当前磁化率模型；

联合反演目标函数构建模块，用于计算获得重力正则化参数和磁法正则化参数；引入联合约束项，通过当前密度模型、当前磁化率模型、联合约束项、重力正则化参数和磁法正则化参数构建联合反演目标函数；

模型更新模块，用于对联合反演目标函数进行优化更新，获得更新后的密度模型、更新后的磁化率模型和数据拟合差；

模型输出模块，若数据拟合差小于目标拟合差，则输出当前的联合反演目标函数，结束联合反演；否则将更新后的密度模型作为当前密度模型，将更新后的磁化率模型作为当前磁化率模型，返回联合反演目标函数构建模块。

本发明具有以下有益效果：

1、采取物性耦合或梯度耦合策略，约束密度与磁化率模型物性参数或其梯度为线性关系；

2、使用LSQR方法求解简化后的联合反演目标方程，减少了反演中的计算量；

3、演中使用了新的步长计算策略，提高了联合反演的稳定性；

4、有效降低了重磁三维联合反演中的非唯一性问题，提高了联合反演结果的可靠程度。

附图说明

图1为本发明实施例方法流程图；

图2为理论模型示意图；

图3为理论模型反演结果；

图4为理论模型反演结果物性或梯度交会图；

本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。

具体实施方式

应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

参照图1，Gramian算子是一个可以将两种或多种物性模型向量联系起来的行列式，将其加入目标函数中进行最小化可实现不同地球物理数据的联合反演；通过一定的转换，Gramian算子可以约束不同模型的物性、梯度或其他关系，也可以基于一定的岩石物理信息进行联合反演；该算法是一种广义的联合反演方法，为解决联合反演所需的巨大计算量问题，以及优化联合反演流程，提高联合反演效率，本发明提供一种基于gramian约束的重磁三维联合反演方法，包括：

S1：构建初始密度模型和初始磁化率模型，将初始密度模型作为当前密度模型，将初始磁化率模型作为当前磁化率模型；

S2：计算获得重力正则化参数和磁法正则化参数；引入联合约束项，通过当前密度模型、当前磁化率模型、联合约束项、重力正则化参数和磁法正则化参数构建联合反演目标函数；

S3：对联合反演目标函数进行优化更新，获得更新后的密度模型、更新后的磁化率模型和数据拟合差；

S4：若数据拟合差小于目标拟合差，则输出当前的联合反演目标函数，结束联合反演；否则将更新后的密度模型作为当前密度模型，将更新后的磁化率模型作为当前磁化率模型，返回步骤S2。

本实施例中，步骤S2具体为：

S21：计算获得重力加权矩阵W_m1，计算公式为：

W_m1＝diag(A₁ ^TA₁)^1/2

式中A₁为重力数据的线性正演算子；

通过重力加权矩阵计算获得重力正则化参数a₁，计算公式为：

式中，b₁和c₁为重力系数参数，x为随机列向量；

S22：计算获得磁法加权矩阵W_m2，计算公式为：

W_m2＝diag(A₂ ^TA₂)^1/2

式中A₂为磁法数据的线性正演算子；

通过磁法加权矩阵计算获得磁法正则化参数，计算公式为：

式中n为迭代次数，b₂和c₂为磁法系数参数，x为随机列向量；

S23：引入的联合约束项φ_G(m^n-1)表示为：

式中m^n-1为联合反演模型，

为当前密度模型，

为当前磁化率模型；G为Gramian算子；

Gramian算子的表达式为：

式中R是转换算子，可以将模型参数转换到其他空间，当模型参数满足以下关系时，Gramian算子等于0：

选择T为单位矩阵或梯度算子矩阵来约束模型参数或其梯度，约束反演模型物性参数或其梯度符合线性关系；

S24：联合反演目标函数φ(m^n-1)的表达式为：

φ(m^n-1)＝φ_d(m^n-1)+a₁φ_m(m^n-1)+a₂φ_G(m^n-1)

式中φ_d(m^n-1)为数据拟合项，φ_m(m^n-1)为深度加权项。

具体的，数据拟合项的表达式为：

其中，

d＝[d₁；d₂]，d₁,d₂分别为重力观测数据和磁法观测数据；

W_d为数据加权矩阵，表达式为：

W_d＝diag(AA^T)^-1/2

深度加权项φ_m(m^n-1)的表达式为：

其中，m_ref为先验模型，反演中被设置为0向量；W_m为深度加权矩阵，表达式为：

W_m＝diag(A^TA)^1/2

本实施例中，步骤S3具体为：

S31：计算Gramian函数值对联合反演目标函数的梯度；

S32：通过高斯-牛顿法最小化联合反演目标函数，获得联合反演目标函数更新方向；

S33：通过联合反演目标函数更新方向计算获得更新后的联合反演目标函数，通过更新后的联合反演目标函数计算获得更新后的密度模型、更新后的磁化率模型和数据拟合差。

本实施例中，Gramian函数值对联合反演目标函数的梯度

的表达式为：

式中L_x,L_y,L_z分别是x、y、z方向的差分算子。

具体的，假设在将三维离散模型转为向量时排列顺序为x、y、z，则L_x表示为：

其中d_x是每个单元x方向的长度，要的注意的是，L_x算子在第kN_x行，k∈R，N_x为x方向单元个数，需使用前一个网格进行差分，同理可得L_y与L_z的表达式。

本实施例中，使用LSQR求解最小二乘问题得到模型更新方向，根据数据拟合情况与联合约束效果得到更新步长，随后更新模型；

使用高斯-牛顿法最小化目标函数即求解以下方程：

其中δm＝[δm₁；δm₂]是模型更新方向。数据拟合项对模型的一阶和二阶导数分别为：

深度加权项对模型的一阶和二阶导数分别为：

对于联合约束项有：

其中D_G是Gramian算子对模型的一阶导数：

将式(13-18)带入式(12)得到：

[(W_dA)^TW_dA+αW_m ^TW_m+γD_G ^TD_G]δm＝-[(W_dA)^TW_d(Am^n-1-d)+αW_m ^TW_m(m^n-1-m_ref)+γD_G ^TG]]

对上式经过数学变换得到联合反演目标函数更新方向；

联合反演目标函数更新方向δm的表达式为：

式中W_d为数据加权矩阵，W_d＝diag(AA^T)^-1/2；W_m为深度加权矩阵，W_m＝diag(A^TA)^1/2；m为联合反演模型；α和γ分别为重力正则化参数和磁法正则化参数；G为Gramian算子；

A₁为重力数据的线性正演算子，A₂为磁法数据的线性正演算子。

本实施例中，更新后的联合反演目标函数mⁿ表示为：

式中n为迭代次数，δm为联合反演目标函数更新方向，I_d为最小化数据拟合差的步长，I_G为最小化Gramian算子的步长；

更新后的密度模型和更新后的磁化率模型的计算公式为：

式中

为更新后的密度模型，

为更新后的磁化率模型。

本实施例中，计算归一化数据拟合差，如果其小于目标拟合差或者达到高斯-牛顿最大迭代次数，则退出循环，输出反演结果，否则回到计算正则化参数这一步进行循环迭代。

为了评价反演模型对数据的拟合程度，数据拟合差misfit的计算公式为：

式中W_d为数据加权矩阵，d_pre为联合预测数据，d_obs为联合观测数据。

本发明提供一种基于gramian约束的重磁三维联合反演系统，包括：

模型构建模块，用于构建初始密度模型和初始磁化率模型，将初始密度模型作为当前密度模型，将初始磁化率模型作为当前磁化率模型；

联合反演目标函数构建模块，用于计算获得重力正则化参数和磁法正则化参数；引入联合约束项，通过当前密度模型、当前磁化率模型、联合约束项、重力正则化参数和磁法正则化参数构建联合反演目标函数；

模型更新模块，用于对联合反演目标函数进行优化更新，获得更新后的密度模型、更新后的磁化率模型和数据拟合差；

模型输出模块，若数据拟合差小于目标拟合差，则输出当前的联合反演目标函数，结束联合反演；否则将更新后的密度模型作为当前密度模型，将更新后的磁化率模型作为当前磁化率模型，返回联合反演目标函数构建模块。

实验验证：

图2展示了理论模型密度(a)与磁化率(b)的三维示意图，反演区域由围岩和四个块状异常体组成，异常体大小皆为0.9km×0.9km×0.6km，围岩的剩余密度为0g/cm³，剩余磁化率为0SI。异常体1剩余密度为1g/cm³，剩余磁化率为0.05SI，中心埋深0.7km；异常体2剩余密度为1g/cm³，剩余磁化率为0.05SI，中心埋深0.7km；异常体3剩余密度为1g/cm³，剩余磁化率为0.05SI，中心埋深1km；异常体4剩余密度为0.6g/cm³，剩余磁化率为0.03SI，中心埋深0.7km。该模型物性分布规律满足χ＝0.05ρ。假设地磁场倾角为90°，偏角为0°，数据采集点布设在地面，采集范围0km≤X≤4km，0km≤Y≤4km，采集间隔为0.1km，共计41×41＝1681个数据点。正演采用合成的重力Gz与磁法TMI数据，并添加2％高斯随机噪声。

根据理论进行联合反演，目标区域大小为4km×4km×2km，将其剖分为80×80×40共256000个单元，每个单元大小为0.05km×0.05km×0.05km，反演中初始密度和磁化率模型均设置为0，正则化衰减系数为3，联合约束项衰减系数为1，高斯牛顿迭代最大次数设置为30次，目标归一化数据拟合差为0.02，求解法方程的LSQR算法最大迭代次数为400，终止误差为0.05。

图3展示了理论模型反演结果在Y＝2.6km切片，(a-b)真实密度和磁化率模型；(c-d)单独反演；(e-f)物性耦合联合反演；(g-h)梯度耦合联合反演的密度与磁化率分布；该切片穿过异常体2和4。可以看到反演的剩余密度在0.7g/cm³左右，剩余磁化率在0.03SI左右，对比单独反演与联合反演的结果，联合反演后的密度模型的数值范围与形态上变化不大，磁化率模型形态发生了变化，其在垂向上的延伸减小，更加符合真实模型。

图4展示了反演物性和梯度交会图，单独反演、物性耦合、梯度耦合联合反演的(a-c)物性交会图；(d-f)X方向梯度交会图；(g-i)Y方向梯度交会图；(j-l)Z方向梯度交会图；图中黑点代表反演结果分布情况，绿色虚线表示真实模型分布情况。可以清楚地看出，通过联合反演后，物性及物性梯度能够更紧密地贴合在绿色虚线上，其分布规律也更加符合真实模型分布情况。

需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。

上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。词语第一、第二、以及第三等的使用不表示任何顺序，可将这些词语解释为标识。

以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (8)

1.一种基于gramian约束的重磁三维联合反演方法，其特征在于，包括：

S1：构建初始密度模型和初始磁化率模型，将初始密度模型作为当前密度模型，将初始磁化率模型作为当前磁化率模型；

S2：计算获得重力正则化参数和磁法正则化参数；引入联合约束项，通过当前密度模型、当前磁化率模型、联合约束项、重力正则化参数和磁法正则化参数构建联合反演目标函数；

S3：对联合反演目标函数进行优化更新，获得更新后的密度模型、更新后的磁化率模型和数据拟合差；

S4：若数据拟合差小于目标拟合差，则输出当前的联合反演目标函数，结束联合反演；否则将更新后的密度模型作为当前密度模型，将更新后的磁化率模型作为当前磁化率模型，返回步骤S2。

2.根据权利要求1所述的基于gramian约束的重磁三维联合反演方法，其特征在于，步骤S2具体为：

S21：计算获得重力加权矩阵W_m1，计算公式为：

W_m1＝diag(A₁ ^TA₁)^1/2

式中A₁为重力数据的线性正演算子；

通过重力加权矩阵计算获得重力正则化参数a₁，计算公式为：

式中，b₁和c₁为重力系数参数，x为随机列向量；

S22：计算获得磁法加权矩阵W_m2，计算公式为：

W_m2＝diag(A₂ ^TA₂)^1/2

式中A₂为磁法数据的线性正演算子；

通过磁法加权矩阵计算获得磁法正则化参数，计算公式为：

式中n为迭代次数，b₂和c₂为磁法系数参数，x为随机列向量；

S23：引入的联合约束项φ_G(m^n-1)表示为：

式中m^n-1为联合反演模型，

为当前密度模型，

为当前磁化率模型；G为Gramian算子；

S24：联合反演目标函数φ(m^n-1)的表达式为：

φ(m^n-1)＝φ_d(m^n-1)+a₁φ_m(m^n-1)+a₂φ_G(m^n-1)

式中φ_d(m^n-1)为数据拟合项，φ_m(m^n-1)为深度加权项。

3.根据权利要求1所述的基于gramian约束的重磁三维联合反演方法，其特征在于，步骤S3具体为：

S31：计算Gramian函数值对联合反演目标函数的梯度；

S32：通过高斯-牛顿法最小化联合反演目标函数，获得联合反演目标函数更新方向；

S33：通过联合反演目标函数更新方向计算获得更新后的联合反演目标函数，通过更新后的联合反演目标函数计算获得更新后的密度模型、更新后的磁化率模型和数据拟合差。

4.根据权利要求3所述的基于gramian约束的重磁三维联合反演方法，其特征在于，Gramian函数值对联合反演目标函数的梯度

的表达式为：

式中L_x,L_y,L_z分别是x、y、z方向的差分算子。

5.根据权利要求3所述的基于gramian约束的重磁三维联合反演方法，其特征在于，联合反演目标函数更新方向δm的表达式为：

式中W_d为数据加权矩阵，W_d＝diag(AA^T)^-1/2；W_m为深度加权矩阵，W_m＝diag(A^TA)^1/2；m为联合反演模型；α和γ分别为重力正则化参数和磁法正则化参数；G为Gramian算子；

A₁为重力数据的线性正演算子，A₂为磁法数据的线性正演算子。

6.根据权利要求3所述的基于gramian约束的重磁三维联合反演方法，其特征在于，更新后的联合反演目标函数mⁿ表示为：

式中n为迭代次数，δm为联合反演目标函数更新方向，I_d为最小化数据拟合差的步长，I_G为最小化Gramian算子的步长；

更新后的密度模型和更新后的磁化率模型的计算公式为：

式中

为更新后的密度模型，

为更新后的磁化率模型。

7.根据权利要求3所述的基于gramian约束的重磁三维联合反演方法，其特征在于，数据拟合差misfit的计算公式为：

式中W_d为数据加权矩阵，d_pre为联合预测数据，d_obs为联合观测数据。

8.一种基于gramian约束的重磁三维联合反演系统，其特征在于，包括：

模型构建模块，用于构建初始密度模型和初始磁化率模型，将初始密度模型作为当前密度模型，将初始磁化率模型作为当前磁化率模型；

联合反演目标函数构建模块，用于计算获得重力正则化参数和磁法正则化参数；引入联合约束项，通过当前密度模型、当前磁化率模型、联合约束项、重力正则化参数和磁法正则化参数构建联合反演目标函数；

模型更新模块，用于对联合反演目标函数进行优化更新，获得更新后的密度模型、更新后的磁化率模型和数据拟合差；

模型输出模块，若数据拟合差小于目标拟合差，则输出当前的联合反演目标函数，结束联合反演；否则将更新后的密度模型作为当前密度模型，将更新后的磁化率模型作为当前磁化率模型，返回联合反演目标函数构建模块。

Images