CN109633762A - 基于正弦函数的相关性约束条件联合反演重磁数据的方法 - Google Patents

Abstract

本发明涉及一种基于正弦函数的相关性约束条件联合反演重磁数据的方法。提出了用密度参数的平方和磁性参数的平方代替密度参数和磁性参数的方法，正弦函数与两个参数向量的叉乘有关，且表达式和求导都更加简单直观，能有效提高联合反演的稳定性。与现有技术相比本发明将正弦函数的相关性约束应用于重磁数据联合反演方法中，在利用共轭梯度算法求解目标函数最优解的过程中，能有效避免奇点的出现，提高了联合反演稳定性。在计算新的相关性约束条件时，用密度参数的平方和磁性参数的平方代替密度参数和磁性参数的方法，使本发明在全局内地质体密度参数和磁性参数不总是同号变化导致的局部相关但全局不相关情况下也适用，故本发明的适用性更广。

Description

基于正弦函数的相关性约束条件联合反演重磁数据的方法

技术领域

本发明涉及一种重磁数据联合反演方法，主要是基于统计学中常规的皮尔森相关性概念，提出一种基于正弦函数的相关性约束条件，并用密度参数的平方和磁性参数的平方代替密度参数和磁性参数进行结构约束，提高了反演的稳定性和适用性。

背景技术：

重磁数据三维物性反演技术蓬勃发展，反演结果对矿产开发、灾害防治等具有重要意义。但是，Tikhonov等人(1977.Methods for solving ill-posed problems(p.12).John Wiley and Sons,Inc.)指出重磁等位场数据反演存在固有的多解性问题，为了降低多解性，可以根据先验信息引入约束条件，比如Li&Oldenburg(19966.3-D inversion ofmagnetic data.Geophysics,61(2):394-408)使用的最大光滑约束、Last&Kubik(1983.Compact gravity inversion.Geophysics,48(6):713-721)使用的最紧致约束等。近年来，通过引入不同种类的地球物理数据联合约束地质模型以降低多解性成为一种比较流行的手段，比较容易实现的是不同种类的地球物理数据对应同一种物性参数，比如Geng等人(2014.3D inversion of airborne gravity-gradiometry data usingcokriging.Geophysics,79(4):G37-G47)利用重力和重力梯度数据联合反演密度异常。另一种情况是，不同种类的地球物理数据对应的物性参数具有一定的相关性，Gardner等人(1974.Formation velocity and density-the diagnostic basics for stratigraphictraps.Geophysics,39:770-780)指出地震数据对应的速度参数和重力数据对应的密度参数具有一定的经验或统计相关，通过这种经验公式达到联合反演的目的。最常见的情况是，不同种类的地球物理数据对应不同的物性参数，而物性参数之间又没有一定的经验数学相关，以上提到的两种方法不再适用。这时，要考虑不同的地球物理数据异常是由同一个地质体产生的，也就是说不同物性参数异常体具有同一地质结构，基于这种理论发展出了两类比较流行的联合反演方法：第一类是Gallardo和Meju(2003.Characterization ofheterogeneous near-surface materials by joint 2D inversion of dc resistivityand seismic data.Geophysical Research Letters,30(13))提出的基于交叉梯度约束的联合反演方法，第二类是Oldenburg and Li(1999.Estimating depth of investigationin dc resistivity and IP surveys.Geophysics,64(2):403-416)以及Peter等人(2012.Joint inversion of seismic travel times and gravity data onunstructured grids with application to mineral exploration.Geophysics,77(1):K1-K15)提出的基于相关性约束的联合反演方法。本发明基于第二类联合反演方法进行了研究和改进，我们注意到常规的基于相关性约束条件本质上与两个参数向量组成的余弦函数有关，表达式中存在分母，这使得反演有出现奇点的风险。另外，常规的基于相关性约束条件是计算密度参数和磁性参数的相关系数，这导致常规的相关性约束方法在有些情况下不适用，例如地质体密度参数和磁性参数在全局内不总是同号导致的局部相关但全局不相关情况，这种全局不相关性限制了常规的基于相关性约束条件的适用性。

发明内容：

本发明的目的是针对上述现有技术的不足，提出一种基于正弦函数的相关性约束条件联合反演重磁数据的方法.

发明思想：正弦函数与两个参数向量的叉乘有关，这种与正弦函数有关的约束项能够避免奇点的出现，且表达式和求导都更加简单直观，能有效提高联合反演的稳定性。同时，选择由密度参数的平方和磁性参数的平方代替密度参数和磁性参数计算新的相关性约束，提高了本发明的适用性。

本发明的目的是通过以下技术方案实现的：

基于正弦函数的相关性约束条件联合反演重磁数据的方法是通过MATLAB平台实现的，包括以下步骤：

a、输入作为观测数据的重力数据g和磁数据T；

b、将正弦函数的相关性约束条件应用于重磁联合反演目标函数中，确定目标函数为：

φ₁＝φ_g+β₁φ_ρ+λ₁φ_JN＝(A₁ρ-g)^T(A₁ρ-g)+β₁(W_ρρ)^T(W_ρρ)+λ₁(|S||R|sinθ)²

φ₂＝φ_T+β₂φ_m+λ₂φ_JN＝(A₂m-T)^T(A₂m-T)+β₂(W_mm)^T(W_mm)+λ₂(|S||R|sinθ)²

其中，g为重力异常数据，ρ为恢复的密度模型，A₁为重力数据正演算子，W_ρ为密度模型的深度加权矩阵，T为磁异常数据，m为恢复的磁化强度模型，A₂为磁数据的正演算子，W_m为磁化强度模型的深度加权矩阵β₁、β₂为正则化因子，λ₁和λ₂为联合反演项平衡因子，其中 和分别为向量X和Y的均值，θ为向量S和R的夹角，其中X＝m²和Y＝ρ²。这个新的目标函数是本发明的核心公式。

c、设置正则化因子β₁、β₂和联合反演项平衡因子λ₁、λ₂的值；

d、设定初始密度参数ρ₀和磁性参数m₀；

e、在加权密度参数域和加权磁性参数域求解，计算ρ_w0＝W_ρρ₀，m_w0＝W_mm₀，A_1W＝A₁W_ρ，A_2W＝A₂W_m；

f、计算目标函数φ₁对加权密度ρ_W0的偏导数：

f_ρ0＝(A_1W ^TA_1W+β₁I)ρ_w0-A_1W ^Tg；计算沿ρ_W0的初始搜索方向：d_ρ0＝-f_ρ0；

计算对应的初始搜索步长：

g、计算目标函数φ₂对加权磁化强度m_W0的偏导数：

f_m0＝(A_2W ^TA_2W+β₁I)m_w0-A_2W ^TT；

计算沿m_W0的初始搜索方向：d_m0＝-f_m0；

计算对应的初始搜索步长：

h、设置最大迭代次数N_max，并记录迭代反演次数k，设置初始k＝1；

i、针对目标函数φ₁，更新加权参数域密度模型ρ_wk＝ρ_wk-1+t_ρk-1d_ρk-1；

j、计算目标函数φ₁对加权密度ρ_wk的一阶偏导数：

f_ρk＝(A_1W ^TA_1W+β₁I)ρ_wk-A_1W ^Tg+λ₁W_ρ ^-1[diag(|S|²R-(SgR)S)]ρ；

对目标函数求偏导是实现反演的关键问题，这个一阶偏导公式是基于本发明提出的目标函数φ₁求得的，为本发明首次计算和使用。

k、根据步骤j获得的偏导数，能够得到沿ρ_wk的搜索方向：

l、与步骤k中的搜索方向所对应的搜索步长为：

m、针对目标函数φ₂，更新加权参数域磁化强度模型m_wk＝m_wk-1+t_mk-1d_mk-1；

n、计算目标函数φ₂对加权磁化强度m_wk的一阶偏导数：

f_mk＝(A_2W ^TA_2W+β₁I)m_wk-A_2W ^TT+λ₂W_m ^-1[diag(|R|²S-(SgR)R)]m；

这个一阶偏导公式是基于本发明提出的目标函数φ₂求得的，也为本发明首次计算和使用。

o、根据步骤n获得的偏导数，能够得到沿m_wk的搜索方向：

p、与步骤o中的搜索方向所对应的搜索步长为：

q、判断当前反演次数小于最大迭代次数是否成立，即k＜N_max是否成立，成立则迭代次数k＝k+1，并重新回到步骤i；否则停止迭代过程，并输出ρ_wk和m_wk；

r、计算ρ＝W_ρ ^-1ρ_wk、m＝W_m ^-1m_wk，最终输出密度参数ρ和磁性参数m。

有益效果：与现有技术相比(1)本发明提出一种基于正弦函数的相关性约束条件的方法，并将这种基于正弦函数的相关性约束应用于重磁数据联合反演方法中，在利用共轭梯度算法求解目标函数最优解的过程中，这种约束方法能够有效避免奇点的出现，提高了联合反演的稳定性。(2)本发明在计算新的相关性约束条件时，提出用密度参数的平方和磁性参数的平方代替密度参数和磁性参数的方法，这使本发明在全局内地质体密度参数和磁性参数不总是同号变化导致的局部相关但全局不相关情况下也适用，故本发明的适用性更广。

附图说明

图1基于正弦函数的相关性约束条件联合反演重磁数据的方法流程图.

图2理论模型及其正演计算的地面观测数据.

(a)理论密度模型；(b)理论磁化强度模型；(c)理论密度模型正演计算的重力异常；(d)理论磁化强度模型正演计算的磁异常.

图3共轭梯度算法分别反演重力异常数据和磁异常数据的结果(a)和(c)是恢复的密度模型的水平切片和垂直切片；(c)和(d)是恢复的磁化强度的水平切片和垂直就切片.(水平切片h＝400m；垂直切片y＝1000m)

图4共轭梯度算法联合反演重力异常数据和磁异常数据的结果.

(a)和(c)是恢复的密度模型的水平切片和垂直切片；(c)和(d)是恢复的磁化强度的水平切片和垂直就切片.(水平切片h＝400m；垂直切片y＝1000m)

具体实施方式

下面结合附图和实例对本发明做进一步的详细的说明。

基于正弦函数的相关性约束条件联合反演重磁数据的方法，包括以下步骤：

a、输入观测到的重力数据g和磁数据T，重力数据g和磁数据T由理论模型给出，理论上的密度模型和磁化强度模型如图2(a)和1(b)，它们由密度模型和磁化强度理论模型正演计算产生，其具体参数为密度模型由5个300m×300m×300m的长方体组合，其中，1号、2号、3号和5号异常体顶部埋深都为200m，4号异常顶面埋深300m，1号、2号和5号异常体的密度差为1g/cm³，3号和4号异常体的密度差是-1g/cm³。磁化强度模型的几何形状、尺寸和位置与图2(a)中的密度模型相同，1号、4号、5号异常体磁化强度异常为1A/m，2号、3号异常体的磁化强度异常为-1A/m。这组密度模型和磁化强度模型的特点是：局部相关，全局几乎不相关，但它们的平方是全局相关的。而本发明提出的利用密度平方和磁化强度平方求取的相关性约束条件的技术适用于这种情况，而常规的直接利用密度和磁化强度求取相关性约束条件的技术在这种情况下不在适用。由理论密度模型产生的地表重力异常如图2(c)所示，假设磁化倾角I＝90°和磁化偏角D＝0°，磁化强度模型产生的地面总磁异常图2(d)所示。图2(c)和图2(d)中的数据采样间距都为100m，分别由665个数据组成，它们作为观测数据的重力数据g和磁数据T进入反演中。

b、将由本发明提出的基于正弦函数的相关性约束条件应用于重磁联合反演目标函数中，确定本发明的目标函数为：

φ₁＝φ_g+β₁φ_ρ+λ₁φ_JN＝(A₁ρ-g)^T(A₁ρ-g)+β₁(W_ρρ)^T(W_ρρ)+λ₁(|S||R|sinθ)² (1)

φ₂＝φ_T+β₂φ_m+λ₂φ_JN＝(A₂m-T)^T(A₂m-T)+β₂(W_mm)^T(W_mm)+λ₂(|S||R|sinθ)² (2)

其中，g为重力异常数据，ρ为恢复的密度模型，A₁为重力数据正演算子，W_ρ为密度模型的深度加权矩阵，T为磁异常数据，m为恢复的磁化强度模型，A₂为磁数据的正演算子，W_m为磁化强度模型的深度加权矩阵β₁、β₂为正则化因子，λ₁和λ₂为联合反演项平衡因子， 和分别为向量X和Y的均值，θ为向量S和R的夹角，X＝m²和Y＝ρ²代替m和ρ。这个新的目标函数是本发明的核心公式。

c、设置正则化因子β₁、β₂和联合反演项平衡因子λ₁、λ₂的值。确定β₁＝2×10⁵、β₂＝1.8×10²，联合反演项平衡因子λ₁＝7×10^-4、λ₂＝9；

d、设定初始密度参数ρ₀＝0和磁性参数m₀＝0；

e、在加权密度参数域和加权磁性参数域求解，计算ρ_w0＝W_ρρ₀，m_w0＝W_mm₀，A_1W＝A₁W_ρ，A_2W＝A₂W_m，其中密度模型的深度加权矩阵W_ρ和磁化强度模型的深度加权矩阵W_m能够有效克服反演结果的趋肤效应。

f、计算目标函数φ₁对加权密度ρ_W0的偏导数：

f_ρ0＝(A_1W ^TA_1W+β₁I)ρ_w0-A_1W ^Tg；计算沿ρ_W0的初始搜索方向：d_ρ0＝-f_ρ0；

计算对应的初始搜索步长：

g、计算目标函数φ₂对加权磁化强度m_W0的偏导数：

f_m0＝(A_2W ^TA_2W+β₁I)m_w0-A_2W ^TT；

计算沿m_W0的初始搜索方向：d_m0＝-f_m0；

计算对应的初始搜索步长：

h、设置最大迭代次数N_max＝50，并记录迭代反演次数k，设置初始k＝1；

i、针对目标函数φ₁，更新加权参数域密度模型ρ_wk＝ρ_wk-1+t_ρk-1d_ρk-1；

j、计算目标函数φ₁对加权密度ρ_wk的一阶偏导数：

f_ρk＝(A_1W ^TA_1W+β₁I)ρ_wk-A_1W ^Tg+λ₁W_ρ ^-1[diag(|S|²R-(SgR)S)]ρ；

这个一阶偏导公式是基于本发明提出的目标函数φ₁求得的，为本发明首次推导计算和使用。

k、根据步骤j获得的偏导数，能够得到沿ρ_wk的搜索方向：

l、与步骤k中的搜索方向所对应的搜索步长为：

m、针对目标函数φ₂，更新加权参数域磁化强度模型m_wk＝m_wk-1+t_mk-1d_mk-1；

n、计算目标函数φ₂对加权磁化强度m_wk的一阶偏导数：

f_mk＝(A_2W ^TA_2W+β₁I)m_wk-A_2W ^TT+λ₂W_m ^-1[diag(|R|²S-(SgR)R)]m；

这个一阶偏导公式是基于本发明提出的目标函数φ₂求得的，也为本发明首次推导计算和使用。

o、根据步骤n获得的偏导数，能够得到沿m_wk的搜索方向：

p、与步骤o中的搜索方向所对应的搜索步长为：

q、判断k＜N_max是否成立，成立则迭代次数k＝k+1，并重新回到步骤i；否则停止迭代过程，并输出ρ_wk和m_wk；

r、计算ρ＝W_ρ ^-1ρ_wk、m＝W_m ^-1m_wk，最终输出密度参数ρ和磁性参数m，成图结果切片如图3所示。

为了清楚表述本发明的效果和作用，在具体实施步骤的第c步中，将联合反演项平衡因子λ₁、λ₂设为0，这时联合反演项不起作用，目标函数(1)和(2)将退化为单独反演的目标函数，其他反演流程不变，将得到单独反演的密度模型和磁化强度模型，如图2，图2(a)-(d)分别展示了单独反演的密度和磁化强度分布的切片。从密度分布切片可以看出，1号和2号(或3号和4号)异常的中间边界很难区分。从恢复的磁化分布切片来看，2号和3号(或4号和5号)异常的中间边界也难以区分。总之，两个相邻且距离相对接近的同性异常体在单独反演中很难区分，无论埋深是否相同。本发明提出的基于正弦函数的相关性约束条件作为约束项在联合反演中发挥作用，将得到基于联合反演恢复的密度模型和磁化强度模型，如图3，图3(a)-(d)分别为恢复密度和磁化分布的切片。联合反演与单独反演对比，可以看出，应用本发明提出的基于正弦的相关性约束条件，密度模型中1号和2号异常以及3号和4号异常之间的界限能够被清晰的识别，磁化强度模型中2号和3号异常以及4号和5号异常之间的界限也变得更清晰。特别是密度模型和磁化模型中的4号异常埋藏较深，有4号异常体引起的观测异常较小，但联合反演仍可将其显现出来。可见，基于本发明的相关性约束联合反演具有明显的效果和优势。

Claims (1)

1.一种基于正弦函数的相关性约束条件联合反演重磁数据的方法，包括以下步骤：

a、输入观测到的重力数据g和磁数据T；

b、根据正弦函数的相关性约束条件确定联合反演重、磁数据的目标函数为：

φ₁＝φ_g+β₁φ_ρ+λ₁φ_JN＝(A₁ρ-g)^T(A₁ρ-g)+β₁(W_ρρ)^T(W_ρρ)+λ₁(|S||R|sinθ)²

φ₂＝φ_T+β₂φ_m+λ₂φ_JN＝(A₂m-T)^T(A₂m-T)+β₂(W_mm)^T(W_mm)+λ₂(|S||R|sinθ)²

其中，g为重力异常数据，ρ为恢复的密度模型，A₁为重力数据正演算子，W_ρ为密度模型的深度加权矩阵，T为磁异常数据，m为恢复的磁化强度模型，A₂为磁数据的正演算子，W_m为磁化强度模型的深度加权矩阵，β₁、β₂为正则化因子，λ₁和λ₂为联合反演项平衡因子，其中 和分别为向量X和Y的均值，θ为向量S与R的夹角，其中X＝m²和Y＝ρ²；

c、设置正则化因子β₁、β₂和联合反演项平衡因子λ₁、λ₂的值；

d、设定初始密度参数ρ₀和磁性参数m₀；

e、在加权密度参数域和加权磁性参数域求解，计算ρ_w0＝W_ρρ₀，m_w0＝W_mm₀，A_1W＝A₁W_ρ，A_2W＝A₂W_m；

f、计算目标函数φ₁对加权密度ρ_W0的偏导数：

f_ρ0＝(A_1W ^TA_1W+β₁I)ρ_W0-A_1W ^Tg，

计算沿ρ_W0的初始搜索方向：d_ρ0＝-f_ρ0，

计算对应的初始搜索步长：

g、计算目标函数φ₂对加权磁化强度m_W0的偏导数：

f_m0＝(A_2W ^TA_2W+β₁I)m_W0-A_2W ^TT，

计算沿m_W0的初始搜索方向：d_m0＝-f_m0，

计算对应的初始搜索步长：

h、设置最大迭代次数N_max，并记录迭代反演次数k，设置初始k＝1；

i、针对目标函数φ₁，更新加权参数域密度模型ρ_wk＝ρ_wk-1+t_ρk-1d_ρk-1；

j、计算目标函数φ₁对加权密度ρ_wk的一阶偏导数：

f_ρk＝(A_1W ^TA_1W+β₁I)ρ_wk-A_1W ^Tg+λ₁W_ρ ^-1[diag(|S|²R-(SgR)S)]ρ；

k、根据步骤j获得的偏导数，得到沿ρ_wk的搜索方向：

l、与步骤k中的搜索方向所对应的搜索步长为：

m、针对目标函数φ₂，更新加权参数域磁化强度模型m_wk＝m_wk-1+t_mk-1d_mk-1；

n、计算目标函数φ₂对加权磁化强度m_wk的一阶偏导数：

f_mk＝(A_2W ^TA_2W+β₁I)m_wk-A_2W ^TT+λ₂W_m ^-1[diag(|R|²S-(SgR)R)]m；

o、根据步骤n获得的偏导数，能够得到沿m_wk的搜索方向：

p、与步骤o中的搜索方向所对应的搜索步长为：

q、判断k＜N_max是否成立，成立则迭代次数k＝k+1，并重新回到步骤i；否则停止迭代过程，并输出ρ_wk和m_wk；

r、计算ρ＝W_ρ ^-1ρ_wk、m＝W_m ^-1m_wk，并输出密度参数ρ和磁性参数m。