CN109358379B - 修正总变分模型约束下基于泛函重构的地球物理反演方法 - Google Patents

Abstract

修正总变分模型约束下基于泛函重构的地球物理反演方法，包括如下步骤：(1)针对修正的总变分模型约束的正则化反演问题，构建反演目标函数；(2)对修正的总变分稳定泛函项进行泛函重构；(3)根据目标函数最小化方程，对模型本身迭代求解新模型的计算公式；(4)设置反演各项参数，重复迭代获得符合期望的模型结果，用于推断地质构造或定位地质异常体；本发明创新性提出泛函重构法，保留了稳定泛函项对模型的约束性质，大幅简化了高度非线性给反演所带来的数值计算难度；反演采用对模型本身迭代的求解方式，相较于对模型校正量迭代的求解方式，正则化直接在模型上施加了稳定泛函约束，可以获得更加合理的模型结果。

Description

修正总变分模型约束下基于泛函重构的地球物理反演方法

技术领域

本发明涉及地球物理勘探技术领域，尤其涉及修正总变分模型约束下基于泛函重构的地球物理反演方法。

背景技术

反演在地球物理领域中占有非常重要的地位。地球物理反演问题根据测量数据或所观测的地球物理场求解场源体，以推断地质结构或定位异常体。反演待求解的模型参数表示地质体的物理性质，可以为密度、磁导率、电导率、弹性、热导率、放射性，与此相对应的勘探方法为重力勘探、磁法勘探、电法勘探、地震勘探、地温法勘探、核法勘探。由于通常观测数据量小于待求解模型参数的数量，地球物理反演问题存在不稳定性和多解性。正则化方法是获得稳定解的重要方法之一。正则化方法在反演目标函数中引入稳定泛函项对模型的解空间进行约束，并通过正则化因子调节拟合差项和稳定泛函项对目标函数的贡献。常见的稳定泛函有最小模型、最大平滑、修正的总变分、最小支撑、最小梯度支撑等。其中，最大平滑稳定泛函在地球物理反演中应用最为广泛，这种模型约束的求解简单、稳定，缺点是对模型参数突变界面刻画模糊；最小梯度支撑稳定泛函是一种有利突出陡变边界的约束方式，缺点是刻画的界面过于陡峭，对模型参数缓慢变化的边界不适应。

从理论上来说，修正的总变分稳定泛函对模型的约束效果介于最大平滑稳定泛函和最小梯度支撑稳定泛函之间，对于模型参数缓慢变化或者具有陡变边界的地下结构，都具有较好的适应性。由于修正的总变分稳定泛函具有高度非线性，带来了许多数值计算困难。研究学者们对该类问题的高效求解法进行了深入研究，现有求解方法中应用较广泛的是交替方向法，该方法将原始优化问题分解成若干子优化问题，对不同的目标变量交替优化，克服了一定的数值计算困难。这种求解方式仍然比其他稳定泛函约束下的反演问题的求解复杂得多，因此，修正的总变分稳定泛函并没有在地球物理反演领域中得到广泛应用。基于上述考虑，针对修正的总变分模型约束的反演问题，需要一种更加高效简单的求解方法，使其充分地应用于地球物理反演中。

发明内容

鉴于上述技术问题，本发明的目的在于提供一种修正总变分模型约束下基于泛函重构的地球物理反演方法，对修正的总变分稳定泛函做泛函重构，大幅简化了反演的数值计算难度；能够更加高效简单地求解修正的总变分模型约束的反演问题。

为了达到上述目的，本发明的技术方案为：

修正总变分模型约束下基于泛函重构的地球物理反演方法，包括以下步骤：

步骤A，针对修正的总变分模型约束的正则化反演问题，构建反演目标函数，目标函数表达式如下：

式中，P^α(m)为目标函数，d为观测数据，m为模型参数，m在地球物理勘探问题中表示地质体的物理性质，包括密度、磁导率、电导率、弹性、热导率或放射性，与此相对应的勘探方法为重力勘探、磁法勘探、电法勘探、地震勘探、地温法勘探或核法勘探，F(m)为正演函数，W_d为数据权重矩阵，▽为计算模型梯度的算子，▽m为模型参数梯度，α为正则化因子，β为不等于零的小数，

为L2范数的平方，|| ||_L1为L1范数，以L2范数表示的函数项为拟合差项，以L1范数表示的函数项为修正的总变分稳定泛函项；

步骤B，对修正的总变分稳定泛函项进行泛函重构，经过重构后，反演目标函数可以表示为如下形式：

式中，W_βTV(m)为变权函数，其表达式如下

上式中，ε为与计算机数值精度有关的一个很小正数；

步骤C，根据目标函数最小化方程，对模型本身迭代求解新模型的计算公式，具体包括：

子步骤C1：令初始模型为m₀，在m₁＝m₀+δm处将正演函数F(m)做泰勒一阶展开，有F(m₁)＝F(m₀)+J₀(m₁-m₀)+o||(δm)²||，其中，δm为关于初始模型m₀的校正矢量，J₀为正演函数F(m)关于m₀的偏导数矩阵，o||(δm)²||为泰勒展开式的高阶余项；

子步骤C2：舍弃F(m₁)的泰勒一阶展开式中高阶余项，将其代入目标函数中，同时，将m₀代入稳定泛函项中变权函数W_βTV(m)中得到固定权系数W_βTV(m₀)，在m₁处，令目标函数P^α(m)关于m₁的梯度为零，得到m₁关于m₀的求解公式；

子步骤C3：重复子步骤C1至C2，得到模型参数的迭代求解公式如下：

式中，m_i+1为第i次迭代后的新模型向量，m_i为第i次迭代的模型向量，J_i为正演函数F(m)关于m_i的偏导数矩阵，T为对矩阵做转置运算，W_d为数据权重矩阵，G为对模型向量做差分运算的矩阵，W_βTVi为将m_i代入变权函数W_βTV(m)中计算得到的加权系数对角矩阵，d为观测数据向量，F(m_i)为将m_i代入正演函数中计算所得的正演结果。

步骤D，设置反演各项参数，重复迭代直到获得符合期望的模型结果，最终模型结果用于推断地质构造或定位地质异常体。

所述的步骤D具体包括：

子步骤D1：读入观测数据d，构造差分运算矩阵G，设置数据权重矩阵W_d、初始模型m₀、正则化因子的数值集{α_n}、变权函数W_βTV(m)中β和ε的数值、反演所要求达到的拟合差φ；

子步骤D2：计算当前模型m_i的正演结果F(m_i)、偏导数矩阵J_i、模型加权系数矩阵W_βTVi；

子步骤D3：将正则化因子的所有数值{α_n}，代入m_i+1的迭代公式中，计算出所有相应的新模型m，筛选出数据拟合差最小的m，作为本次迭代的最终模型m_i+1；

子步骤D4：判断新模型m_i+1的数据拟合差是否满足小于反演所要求达到的拟合差φ的条件，若不满足，迭代继续，再次运行子步骤D2至D3；若满足，迭代终止，输出最终模型参数。

从上述技术方案可以看出，本发明修正的总变分模型约束下基于泛函重构的参数反演方法，具有以下有益效果：

(1)对修正的总变分稳定泛函做泛函重构，大幅简化了反演的数值计算难度；

(2)反演采用对模型本身迭代的求解方式，相较于对模型校正量迭代的求解方式，正则化直接在模型上施加了稳定泛函约束，可以获得更为合理的解。

附图说明

图1为根据本发明实施例修正的总变分模型约束下基于泛函重构的参数反演方法的步骤图。

图2A为本实施例中对瞬变电磁仿真数据的反演最终模型与预设真实模型对比图。

图2B为本实施例中由反演最终模型计算的正演结果与仿真数据对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的图号。附图中未绘示或描述的实现方式，为所属技术领域中普通技术人员所知的形式。另外，虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向。因此，使用的方向用语是用来说明并非用来限制本发明的保护范围。

在本发明的一个示例性实施例中，提供了对一组瞬变电磁仿真数据进行修正的总变分模型约束下基于泛函重构的地球物理反演方法的演示。图1是根据本发明实施例修正的总变分模型约束下基于泛函重构的地球物理反演方法的步骤图。请参照图1，本实施例包括：

步骤A，针对修正的总变分模型约束的正则化反演问题，构建反演目标函数。目标函数表达式如下：

式中，P^α(m)为目标函数，d为观测数据，m为模型参数，F(m)为正演函数，W_d为数据权重矩阵，▽为计算模型梯度的算子，▽m为模型参数梯度，α为正则化因子，β为不等于零的小数，

为L2范数的平方，||||_L1为L1范数，以L2范数表示的函数项为拟合差项，以L1范数表示的函数项为修正的总变分稳定泛函项。

在本实施例中，基于层状大地模型进行瞬变电磁数据仿真，将仿真数据表示为d＝[d₁,d₂,…,d_L]，L为数据总点数，将待反演的层状大地电阻率模型参数表示为m＝[m₁,m₂,…,m_N]，m_n对应深度z处的电阻率，z_n-1＜z＜z_n，n＝1,2,…,N，z₀＝0，N为总层数，通常取值为20-100，数据权重矩阵W_d可以表示为W_d＝diag{1/σ₁,1/σ₂,…,1/σ_L}，σ_l为第l个数据的数据误差，F(m)为层状大地模型的瞬变电磁正演函数。

步骤B，对修正的总变分稳定泛函项进行泛函重构。

关于修正的总变分稳定泛函s(m)的重构方式如下：

式中，ε为与计算机数值精度有关的一个很小正数，W_βTV(m)为变权函数，其表达式如下：

经过重构后，反演目标函数可以表示为如下形式：

步骤C，根据目标函数最小化方程，对模型本身迭代求解新模型的迭代公式。具体包括：

子步骤C1：令初始模型为m₀，在m₁＝m₀+δm处将正演函数F(m)做泰勒一阶展开，有F(m₁)＝F(m₀)+J₀(m₁-m₀)+o||(δm)²||，其中，δm为关于初始模型m₀的校正矢量，J₀为正演函数F(m)关于m₀的偏导数矩阵，o||(δm)²||为泰勒展开式的高阶余项；

子步骤C2：舍弃F(m₁)的泰勒一阶展开式中高阶余项，将其代入目标函数中，同时，将m₀代入稳定泛函项中变权函数W_βTV(m)中得到固定权系数W_βTV(m₀)，在m₁处，令目标函数P^α(m)关于m₁的梯度为零，可以得到m₁关于m₀的求解公式；

子步骤C3：重复子步骤C1至C2，可以得到模型参数的迭代求解公式如下：

式中，m_i+1为第i次迭代后的新模型矢量，m_i为第i次迭代的模型矢量，J_i为正演函数F(m)关于m_i的偏导数矩阵，T为对矩阵做转置运算，W_d为数据权重矩阵，G为对模型矢量做差分运算的矩阵，W_βTVi为将m_i代入变权函数W_βTV(m)中计算得到的加权系数对角矩阵，d为观测数据向量，F(m_i)为将m_i代入正演函数中计算所得的正演结果。

步骤D，设置反演各项参数，重复迭代直到获得符合期望的模型结果。具体包括：

子步骤D1：读入观测数据d，构造差分运算矩阵G，设置数据权重矩阵W_d、初始模型m₀、正则化因子的数值集{α_n}、变权函数W_βTV(m)中β和ε的数值、反演所要求达到的拟合差φ；

在本实施例中，待反演的层状大地电阻率模型参数为m＝[m₁,m₂,…,m_N]，为实现差分运算，构造的差分运算矩阵G表达式具体如下：

子步骤D2：计算当前模型m_i的正演结果F(m_i)、偏导数矩阵J_i、模型加权系数矩阵W_βTVi；

子步骤D3：将正则化因子的所有数值{α_n}，代入m_i+1的迭代公式中，计算出所有相应的新模型m，筛选出数据拟合差最小的m，作为本次迭代的最终模型m_i+1；

子步骤D4：判断新模型m_i+1的数据拟合差是否满足小于反演所要求达到的拟合差φ的条件，若不满足，迭代继续，再次运行子步骤D2至D3；若满足，迭代终止，输出最终模型参数。

图2A为本实施例中对瞬变电磁仿真数据的反演最终模型与预设真实地质模型对比图，图中采用对数坐标，横轴为深度，纵轴为电阻率值，红线为预设真实模型，蓝线为反演最终模型。从图2A可以看出，本发明的反演最终模型符合修正的总变分稳定泛函对模型的约束性质，并且与真实模型非常接近，符合反演的准确性要求。图2B为本实施例中由反演最终模型计算的正演结果与仿真数据对比图，图中采用对数坐标，横轴为时间，纵轴为磁场响应幅度。从图2B可以看出，由反演最终模型计算的正演结果F(m_last)与预设仿真数据d基本完全重合，这说明最终数据拟合差非常小，符合反演的准确性要求，其中，m_last表示反演最终模型。

结合本发明的具体实施步骤，可以发现，本发明涉及的数值计算和迭代过程相对简单。

至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明修正的总变分模型约束下基于泛函重构的参数反演方法有了清楚的认识。

此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换。

综上所述，本发明修正的总变分模型约束下基于泛函重构的地球物理反演方法，通过泛函重构法，保留了稳定泛函项对模型的约束性质，大幅简化了高度非线性给反演所带来的数值计算难度；反演采用对模型本身迭代的求解方式，相较于对模型校正量迭代的求解方式，正则化直接在模型上施加了稳定泛函约束，可以获得更加合理的最终模型。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.修正总变分模型约束下基于泛函重构的地球物理反演方法，其特征在于，包括以下步骤：

步骤A，针对修正的总变分模型约束的正则化反演问题，构建反演目标函数，目标函数表达式如下：

式中，P^α(m)为目标函数，d为观测数据，m为模型参数，m在地球物理勘探问题中表示地质体的物理性质，包括密度、磁导率、电导率、弹性、热导率或放射性，与此相对应的勘探方法为重力勘探、磁法勘探、电法勘探、地震勘探、地温法勘探或核法勘探，F(m)为正演函数，W_d为数据权重矩阵，

为计算模型梯度的算子，

为模型参数梯度，α为正则化因子，β为不等于零的小数，

为L2范数的平方，||||_L1为L1范数，以L2范数表示的函数项为拟合差项，以L1范数表示的函数项为修正的总变分稳定泛函项；

步骤B，对修正的总变分稳定泛函项进行泛函重构，经过重构后，反演目标函数可以表示为如下形式：

式中，W_βTV(m)为变权函数，其表达式如下

上式中，ε为与计算机数值精度有关的一个很小正数；

步骤C，根据目标函数最小化方程，对模型本身迭代求解新模型的计算公式；

步骤D，设置反演各项参数，重复迭代直到获得符合期望的模型结果，最终模型结果用于推断地质构造或定位地质异常体。

2.根据权利要求1所述的修正总变分模型约束下基于泛函重构的地球物理反演方法，其特征在于，所述的步骤C具体包括：

子步骤C1：令初始模型为m₀，在m₁＝m₀+δm处将正演函数F(m)做泰勒一阶展开，有F(m₁)＝F(m₀)+J₀(m₁-m₀)+o||(δm)²||，其中，δm为关于初始模型m₀的校正矢量，J₀为正演函数F(m)关于m₀的偏导数矩阵，o||(δm)²||为泰勒展开式的高阶余项；

子步骤C2：舍弃F(m₁)的泰勒一阶展开式中高阶余项，将其代入目标函数中，同时，将m₀代入稳定泛函项中变权函数W_βTV(m)中得到固定权系数W_βTV(m₀)，在m₁处，令目标函数P^α(m)关于m₁的梯度为零，得到m₁关于m₀的求解公式；

子步骤C3：重复子步骤C1至C2，得到模型参数的迭代求解公式如下：

式中，m_i+1为第i次迭代后的新模型向量，m_i为第i次迭代的模型向量，J_i为正演函数F(m)关于m_i的偏导数矩阵，T为对矩阵做转置运算，W_d为数据权重矩阵，G为对模型向量做差分运算的矩阵，W_βTVi为将m_i代入变权函数W_βTV(m)中计算得到的加权系数对角矩阵，d为观测数据向量，F(m_i)为将m_i代入正演函数中计算所得的正演结果。

3.根据权利要求1所述的修正总变分模型约束下基于泛函重构的地球物理反演方法，其特征在于，所述的步骤D具体包括：

子步骤D1：读入观测数据d，构造差分运算矩阵G，设置数据权重矩阵W_d、初始模型m₀、正则化因子的数值集{α_n}、变权函数W_βTV(m)中β和ε的数值、反演所要求达到的拟合差φ；

子步骤D2：计算当前模型m_i的正演结果F(m_i)、偏导数矩阵J_i、模型加权系数矩阵W_βTVi；

子步骤D3：将正则化因子的所有数值{α_n}，代入m_i+1的迭代公式中，计算出所有相应的新模型m，筛选出数据拟合差最小的m，作为本次迭代的最终模型m_i+1；

子步骤D4：判断新模型m_i+1的数据拟合差是否满足小于反演所要求达到的拟合差φ的条件，若不满足，迭代继续，再次运行子步骤D2至D3；若满足，迭代终止，输出最终模型参数。

Images