CN110244351A - 一种不同约束地球物理反问题的统一构造反演方法 - Google Patents

Abstract

一种不同约束地球物理反问题的统一构造反演方法，包括如下步骤：A、在先验地质信息的基础上，选择稳定泛函，构建正则化反演目标函数；B、将稳定泛函构造为核函数和权函数的组合；C、构造统一形式的目标函数并求解其最小化问题；D、设置反演各项参数，重复迭代求解新模型，从而推断地质构造或定位地质异常体；本发明对不同约束地球物理反问题提供了一种统一构造反演方法，其中，统一构造形式使不同稳定泛函的异同更加明确，为反演中选择最优稳定泛函提供了理论依据；统一构造反演方法使不同约束地球物理反问题在统一框架下处理，有利于在实际数据处理中高效灵活地测试不同约束的反演效果，以获得最符合真实地电结构的结果。

Description

一种不同约束地球物理反问题的统一构造反演方法

技术领域

本发明涉及地球物理技术领域，尤其涉及一种不同约束地球物理反问题的统一构造反演方法。

背景技术

地球物理反问题是根据测量数据或所观测的地球物理场求解场源体，以推断地质结构或定位异常体。反问题待求解的模型参数表示地质体的物理性质，根据不同类型的地球物理勘探方法，其可以为密度、磁导率、电阻率率、弹性、热导率、放射性等。绝大部分地球物理反问题都是不适定的，为了克服不适定性并获得稳定解，正则化方法发挥着重要作用。在正则化理论的框架下，为了获得符合真实地电结构的反演结果，先验地质信息通过稳定泛函的形式包含在反演过程中。那么，这样的地球物理反问题有约束的，采用不同的稳定泛函进行模型约束，便构成了不同约束的地球物理反问题。目前，先验地质信息是反演中选择稳定泛函的前提条件。在此基础上，有些稳定泛函的性质是相似的，同时，不同约束地球反问题缺乏在统一架构下的理论对比与结果对比的依据，因此，如何选择最优稳定泛函以及在最优稳定泛函基础上获得最佳反演结果，仍然是待解决的难题。

发明内容

鉴于上述技术问题，本发明的目的在于提供了一种不同约束地球物理反问题的统一构造反演方法，以获得最优稳定泛函及最准确反映真实地电结构的结果。

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

一种不同约束地球物理反问题的统一构造反演方法，包括以下步骤：

步骤A，在先验地质信息的基础上，选择稳定泛函，构建正则化反演目标函数,目标函数表达式如下：

式中，P^α(m)为目标函数，为L2范数的平方和，d为观测数据，m为模型参数，m在地球物理反演问题中表示地质体的物理性质，包括密度、磁导率、电阻率、弹性、热导率或放射性，与此相对应的勘探方法为重力勘探、磁法勘探、电法勘探、地震勘探、地温法勘探或核法勘探，F(m)为正演函数，W_d为数据权重矩阵，α为正则化因子，s(m)为根据先验地质信息选择的约束模型构造的稳定泛函。

将不同稳定泛函统一采用二范数表示，具体形式如下：

(1)最小模型泛函s_min(m)

(2)一阶导数最光滑泛函s_maxG(m)

(3)二阶导数最光滑泛函s_maxL(m)

(4)修正的总变分泛函s_βTV(m)

(5)最小支撑约束泛函s_βMS(m)

(6)最小梯度支撑泛函s_βMGS(m)

其中，m为模型参数，为梯度运算符，为拉普拉斯运算符，β为模型聚焦调节因子。

步骤B，将稳定泛函构造为核函数和权函数的组合，不同稳定泛函的统一构造形式如下：

式中，w_s(m)为稳定泛函的权函数，k(m)为稳定泛函的核心函数。

将以上统一构造形式应用于步骤A中所介绍的几种稳定泛函，有：

(1)最小模型泛函s_min(m)的核函数和权函数

k(m)＝m

w_s(m)＝1

(2)一阶导数最光滑泛函s_maxG(m)的核函数和权函数

w_s(m)＝1

(3)二阶导数最光滑泛函s_maxG(m)的核函数和权函数

w_s(m)＝1

(4)修正的总变分泛函s_βTV(m)的核函数和权函数

(5)最小支撑泛函s_βMS(m)的核函数和权函数

k(m)＝m

(6)最小梯度支撑泛函s_βMGS(m)的核函数和权函数

其中，ε是与计算机数值精度有关的一个很小正数。

步骤C，构造统一形式的目标函数并求解其最小化问题，统一形式的目标函数如下：

以上目标函数最小化问题的求解过程具体包括：

子步骤C1：在统一构造的目标函数中，令稳定泛函的核函数为k(m)＝m；

子步骤C2：令初始模型为m₀，在m₁＝m₀+δm处，其中δm为关于初始模型m₀的校正矢量，将正演函数F(m₁)做泰勒一阶展开，舍弃高阶余项后，将泰勒一阶展开式F(m₁)＝F(m₀)+J₀(m₁-m₀)代入目标函数中，其中J₀是正演函数在m₀处的偏导数矩阵。同时，将m₀代入权函数w_s(m)中计算得到的权系数对角矩阵W_s0，在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_si为将m_i代入权函数w_s(m)中计算得到的加权系数对角矩阵，d为观测数据向量，F(m_i)为将m_i代入正演函数中计算所得的正演结果。

子步骤C4：针对稳定泛函的核函数是的情况，在子步骤C3的基础上，模型参数的迭代求解公式如下：

式中，G为梯度算子的矩阵表达式，对于m＝[m₁,m₂,…,m_N]，梯度矩阵G如下：

子步骤C5：针对稳定泛函的核函数是的情况，在子步骤C4的基础上，模型参数的迭代求解公式如下：

L＝G²

式中，L为拉普拉斯运算符的矩阵表达式。

步骤D，设置反演各项参数，重复迭代求解新模型，直到收敛，输出最终模型参数，然后图形化显示最终模型参数，观察到不同位置地质体的物理参数大小，从而推断地质构造或定位地质异常体。

步骤D具体包括：

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

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

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

子步骤D4：判断新模型m_i+1的数据拟合差是否满足小于反演所要求达到的拟合差φ的条件，若不满足，迭代继续，再次运行子步骤D2至D3；若满足，迭代终止，输出最终模型参数，然后图形化显示最终模型参数，可以观察到不同位置地质体的物理参数大小，从而推断地质构造或定位地质异常体。

(三)有益效果

从上述技术方案可以看出，本发明不同约束地球物理反问题的统一构造反演方法，具有以下有益效果：

(1)统一构造形式使得不同稳定泛函的异同更加明确，为反演中选择最优稳定泛函提供了理论依据；

(2)对不同稳定泛函进行统一构造，使不同约束地球物理反演问题可以在统一框架下进行求解；

(3)统一构造求解法有利于在实际数据处理中灵活地测试不同约束的反演效果，以获得最符合真实地电结构的结果。

附图说明

图1为根据本发明实施例不同约束地球物理反问题的统一构造反演方法的步骤图。

图2A为本实施例中基于最小模型泛函、最小支撑泛函的反演结果与真实模型对比图。

图2B为本实施例中基于一阶导数最光滑泛函、二阶导数最光滑泛函的反演结果与真实模型对比图。

图2C为本实施例中基于修正的总变分泛函、最小梯度支撑泛函的反演结果与真实模型对比图。

具体实施方式

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

在本发明的一个示例性实施例中，提供了对一组大地电磁仿真数据进行不同约束地球物理反问题的统一构造反演方法的演示。图1是根据本发明实施例不同约束地球物理反问题的统一构造反演方法的步骤图。请参照图1，本实施例包括：

步骤A，在先验地质信息的基础上，选择稳定泛函，构建正则化反演目标函数。目标函数表达式如下：

式中，P^α(m)为目标函数，为L2范数的平方和，d为观测数据，m为模型参数，m在地球物理反演问题中表示地质体的物理性质，包括密度、磁导率、电阻率、弹性、热导率或放射性，与此相对应的勘探方法为重力勘探、磁法勘探、电法勘探、地震勘探、地温法勘探或核法勘探，F(m)为正演函数，W_d为数据权重矩阵，α为正则化因子，s(m)为根据先验地质信息选择的约束模型构造的稳定泛函。

将不同稳定泛函统一采用二范数表示，具体形式如下：

(1)最小模型泛函s_min(m)

(2)一阶导数最光滑泛函s_maxG(m)

(3)二阶导数最光滑泛函s_maxL(m)

(4)修正的总变分泛函s_βTV(m)

(5)最小支撑约束泛函s_βMS(m)

(6)最小梯度支撑泛函s_βMGS(m)

其中，m为模型参数，为梯度运算符，为拉普拉斯运算符，β为模型聚焦调节因子。

在本实施例中，基于层状大地模型进行大地电磁数据仿真，将仿真数据表示为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为总层数，数据权重矩阵W_d可以表示为W_d＝diag{1/σ₁,1/σ₂,…,1/σ_L}，σ_l为第l个数据的数据误差，F(m)为层状大地模型的大地电磁正演函数，这里只采用视电阻率数据作为反演数据。为了检验本发明方法的有效性，我们将以上列举的6种不同稳定泛函约束的反演问题都依次进行求解并展示反演结果。

步骤B，将稳定泛函构造为核函数和权函数的组合。不同稳定泛函的统一构造形式如下：

式中，w_s(m)为稳定泛函的权函数，k(m)为稳定泛函的核心函数。

将以上统一构造形式应用于步骤A中所介绍的几种稳定泛函，有：

(1)最小模型泛函s_min(m)的核函数和权函数

k(m)＝m

w_s(m)＝1

(2)一阶导数最光滑泛函s_maxG(m)的核函数和权函数

w_s(m)＝1

(3)二阶导数最光滑泛函s_maxG(m)的核函数和权函数

w_s(m)＝1

(4)修正的总变分泛函s_βTV(m)的核函数和权函数

(5)最小支撑泛函s_βMS(m)的核函数和权函数

k(m)＝m

(6)最小梯度支撑泛函s_βMGS(m)的核函数和权函数

其中，ε是与计算机数值精度有关的一个很小正数。

在本实施例中，为了更好地展示6种稳定泛函约束的反演结果，根据稳定泛函的重构形式，将具有相似性的稳定泛函分组在一起以对比它们的反演结果，因此，基于最小模型泛函、最小支撑泛函的反演结果为一组，基于一阶导数最光滑泛函、二阶导数最光滑泛函的反演结果为一组，基于修正的总变分泛函、最小梯度支撑泛函的反演结果为一组。从这里可以看出，稳定泛函的统一构造形式从理论上明确了不同稳定泛函异同点。

步骤C，构造统一形式的目标函数并求解其最小化问题。统一形式的目标函数如下：

以上目标函数最小化问题的求解过程具体包括：

子步骤C1：在统一构造的目标函数中，令稳定泛函的核函数为k(m)＝m；

子步骤C2：令初始模型为m₀，在m₁＝m₀+δm处，其中δm为关于初始模型m₀的校正矢量，将正演函数F(m₁)做泰勒一阶展开，舍弃高阶余项后，将泰勒一阶展开式F(m₁)＝F(m₀)+J₀(m₁-m₀)代入目标函数中，其中J₀是正演函数在m₀处的偏导数矩阵。同时，将m₀代入权函数w_s(m)中计算得到的权系数对角矩阵W_s0，在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_si为将m_i代入权函数w_s(m)中计算得到的加权系数对角矩阵，d为观测数据向量，F(m_i)为将m_i代入正演函数中计算所得的正演结果。

子步骤C4：针对稳定泛函的核函数是的情况，在子步骤C3的基础上，模型参数的迭代求解公式如下：

式中，G为梯度算子▽的矩阵表达式，对于m＝[m₁,m₂,…,m_N]，梯度矩阵G如下：

子步骤C5：针对稳定泛函的核函数是的情况，在子步骤C4的基础上，模型参数的迭代求解公式如下：

L＝G²

式中，L为拉普拉斯运算符的矩阵表达式。

步骤D，设置反演各项参数，重复迭代求解新模型，直到收敛。最终模型结果用于推断地质构造或定位地质异常体。具体包括：

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

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

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

子步骤D4：判断新模型m_i+1的数据拟合差是否满足小于反演所要求达到的拟合差φ的条件，若不满足，迭代继续，再次运行子步骤D2至D3；若满足，迭代终止，输出最终模型参数，然后图形化显示最终模型参数，可以观察到不同位置地质体的物理参数大小，从而推断地质构造或定位地质异常体。

图2A-2C为本实施例中，针对大地电磁仿真数据，基于6种不同稳定泛函约束的反演结果与预设真实模型对比图，图中横轴为电阻率值并采用对数坐标，纵轴为深度且为实数坐标。其中，图2A为基于最小模型泛函、最小支撑泛函的反演结果与真实模型对比图，在图中依次标识为Minref、MSref、True Model。图2B为基于一阶导数最光滑泛函、二阶导数最光滑泛函的反演结果与真实模型对比图，在图中依次标识为MaxG、MaxL、True Model。图2C为基于修正的总变分泛函、最小梯度支撑泛函的反演结果与真实模型对比图，在图中依次标识为MTV、MGS、True Model。收敛要求的拟合差对这6组反演设为一致，均在10次迭代内达到收敛。从图2可见，同一分组的反演结果确实具有相似性，反映了本发明对稳定泛函统一构造的优势，即这种构造形式从理论上就突出了不同稳定泛函的异同；同时，反演结果也验证了本发明的统一构造求解法能够求解不同约束反演问题，并且这种统一构造求解法有利于高效灵活地测试不同约束的反演效果。综合考虑理论依据、反演过程、最终效果等各方面因素，可以获得最优稳定泛函及最准确反映地电结构的结果。

至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明不同约束地球物理反问题的统一构造反演方法有了清楚的认识。

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

综上所述，本发明对不同约束地球物理反问题提供了一种统一构造反演方法，其中，统一构造形式使得不同稳定泛函的异同更加明确，为反演中选择最优稳定泛函提供了理论依据；统一构造反演方法使不同约束地球物理反问题在统一框架下处理，有利于在实际数据处理中高效灵活地测试不同约束的反演效果，以获得最符合真实地电结构的结果。

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

Claims (5)

1.一种不同约束地球物理反问题的统一构造反演方法，其特征在于，包括以下步骤：

步骤A，在先验地质信息的基础上，选择稳定泛函，构建正则化反演目标函数,目标函数表达式如下：

式中，P^α(m)为目标函数，为L2范数的平方和，d为观测数据，m为模型参数，m在地球物理反演问题中表示地质体的物理性质，包括密度、磁导率、电阻率、弹性、热导率或放射性，与此相对应的勘探方法为重力勘探、磁法勘探、电法勘探、地震勘探、地温法勘探或核法勘探，F(m)为正演函数，W_d为数据权重矩阵，α为正则化因子，s(m)为根据先验地质信息选择的约束模型构造的稳定泛函；

步骤B，将稳定泛函构造为核函数和权函数的组合，不同稳定泛函的统一构造形式如下：

式中，w_s(m)为稳定泛函的权函数，k(m)为稳定泛函的核心函数；

步骤C，构造统一形式的目标函数并求解其最小化问题，统一形式的目标函数如下：

步骤D，设置反演各项参数，重复迭代求解新模型，直到收敛，输出最终模型参数，然后图形化显示最终模型参数，观察到不同位置地质体的物理参数大小，从而推断地质构造或定位地质异常体。

2.根据权利要求1所述的一种不同约束地球物理反问题的统一构造反演方法，其特征在于，所述步骤A具体为：

将不同稳定泛函统一采用二范数表示，具体形式如下：

(1)最小模型泛函s_min(m)

(2)一阶导数最光滑泛函s_maxG(m)

(3)二阶导数最光滑泛函s_maxL(m)

(4)修正的总变分泛函s_βTV(m)

(5)最小支撑约束泛函s_βMS(m)

(6)最小梯度支撑泛函s_βMGS(m)

其中，m为模型参数，▽为梯度运算符，▽²为拉普拉斯运算符，β为模型聚焦调节因子。

3.根据权利要求1所述的一种不同约束地球物理反问题的统一构造反演方法，其特征在于，所述步骤B中稳定泛函的统一构造形式应用于步骤A的稳定泛函，有：

(1)最小模型泛函s_min(m)的核函数和权函数

k(m)＝m

w_s(m)＝1

(2)一阶导数最光滑泛函s_maxG(m)的核函数和权函数

k(m)＝▽m

w_s(m)＝1

(3)二阶导数最光滑泛函s_maxG(m)的核函数和权函数

k(m)＝▽²m

w_s(m)＝1

(4)修正的总变分泛函s_βTV(m)的核函数和权函数

k(m)＝▽m

(5)最小支撑泛函s_βMS(m)的核函数和权函数

k(m)＝m

(6)最小梯度支撑泛函s_βMGS(m)的核函数和权函数

k(m)＝▽m

其中，ε是与计算机数值精度有关的一个很小正数。

4.根据权利要求1所述的一种不同约束地球物理反问题的统一构造反演方法，其特征在于，所述步骤C的目标函数最小化求解过程具体包括：

子步骤C1：在统一构造的目标函数中，令稳定泛函的核函数为k(m)＝m；

子步骤C2：令初始模型为m₀，在m₁＝m₀+δm处，其中δm为关于初始模型m₀的校正矢量，将正演函数F(m₁)做泰勒一阶展开，舍弃高阶余项后，将泰勒一阶展开式F(m₁)＝F(m₀)+J₀(m₁-m₀)代入目标函数中，其中J₀是正演函数在m₀处的偏导数矩阵；同时，将m₀代入权函数w_s(m)中计算得到的权系数对角矩阵W_s0，在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_si为将m_i代入权函数w_s(m)中计算得到的加权系数对角矩阵，d为观测数据向量，F(m_i)为将m_i代入正演函数中计算所得的正演结果；

子步骤C4：针对稳定泛函的核函数是k(m)＝▽m的情况，在子步骤C3的基础上，模型参数的迭代求解公式如下：

式中，G为梯度算子▽的矩阵表达式，对于m＝[m₁,m₂,…,m_N]，梯度矩阵G如下：

子步骤C5：针对稳定泛函的核函数是k(m)＝▽²m的情况，在子步骤C4的基础上，模型参数的迭代求解公式如下：

L＝G²

式中，L为拉普拉斯运算符▽²的矩阵表达式。

5.根据权利要求1所述的一种不同约束地球物理反问题的统一构造反演方法，其特征在于，所述步骤D具体包括：

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

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

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

子步骤D4：判断新模型m_i+1的数据拟合差是否满足小于反演所要求达到的拟合差φ的条件，若不满足，迭代继续，再次运行子步骤D2至D3；若满足，迭代终止，输出最终模型参数，然后图形化显示最终模型参数，可以观察到不同位置地质体的物理参数大小，从而推断地质构造或定位地质异常体。