CN110007365A - 一种基于信号数据稀疏空间快速计算的联合反演方法 - Google Patents

Abstract

本发明公开了一种基于信号数据稀疏空间快速计算的联合反演方法，包括如下步骤：S1.建立重磁联合反演目标函数；S2.建立密度和磁化率的初始模型；S3.对密度、磁化率进行正演计算；S4.计算交叉梯度偏导数和拉格朗日算子；S5.计算获得稀疏信号数据空间迭代模型；S6.对获得的稀疏信号数据空间迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出该数据空间联合反演模型，否则执行S3；本发明计算速度快、占用内存小，为重磁数据处理与解释提供了一种行之有效的计算方法。

Description

一种基于信号数据稀疏空间快速计算的联合反演方法

技术领域

本发明涉及地球科学及物理学交叉领域中，重力数据磁数据联合反演的处理技术，具体地说，是涉及一种基于信号数据稀疏空间快速计算的联合反演方法。

背景技术

地球物理方法已经广泛应用于深部地质探测和矿产、水、石油、地热等资源勘查，由于反演方法本身存在的多解性，单一反演方法存在着诸多问题，进而联合反演成为综合地球物理方法数据处理的主流方法。目前，科学家已经提出了多种物性参数之间结构藕合的联合反演方法，国内也实现了二种或多种物性参数之间的联合反演方法。

由于地球物理数据采集中，需要开展大数据量的处理解释，进而满足生产需要。对于大数据的解释，往往基于模型空间下进行反演计算。但是当反演物性参数种类由单一变成多种数据时，涉及的计算过程中的数据量将巨大增多，进而对计算机的计算速度和内存占用都将成为考验。该问题已经对综合地球物理方法的发展形成了挑战。

发明内容

本发明的目的在于提供一种基于信号数据稀疏空间快速计算的联合反演方法，以解决现有的联合反演方法计算速度慢、占用内存大等问题，同时更贴近实际地质模型。

为实现上述目的，本发明采用的技术方案如下：

一种基于信号数据稀疏空间快速计算的联合反演方法，包括如下步骤：

S1.建立重磁联合反演目标函数，所述联合反演目标函数如式(1)：

约束条件τ＝0；

(1)式中：

其中，角标g，m分别为密度、磁化率有关参数及模型，为对应的目标函数，W_d为观测数据协方差矩阵，A为正演响应，m为密度或磁化率模型参数，d为观测数据，α为正则化参数，为基于信号数据稀疏先验加权矩阵，W_m为模型参数的协方差矩阵，m_apr为先验模型；T为矩阵的转置，为梯度，τ为密度模型及磁化率模型的交叉梯度矩阵；

S2.建立密度、磁化率初始模型：

在笛卡尔坐标系下，二维数据处理中，沿x，z两个坐标轴分别将初始模型空间划分成为N_x、N_z个小的矩形网格单元，间距为：Δx_i(i＝1，2，…N_x)，Δz_i(i＝1，2，…N_z)；三维情况下将初始模型空间划分成为N_x、N_y、N_z个小的矩形网格单元，间距为：Δx_i(i＝1，2，…N_x)，Δy(i＝1，2，…N_y)，Δz_i(i＝1，2，…N_z)；

S3.对密度、磁化率进行正演计算：

采用Pick(1973)及Fu(2015)提出的立方体网格剖分下的密度及磁化率正演公式，即密度正演响应矩阵A_g的元素的计算公式为：

为A_g中的元素，G为万有引力常数，x，y，z为观测点的坐标，为第j个小立方体的x方向上的边界坐标，为第j个小立方体的y方向上的边界坐标，为第j个小立方体的z方向上的边界坐标；

磁化率正演响应矩阵A_mag的元素的计算公式为：

其中μ₀真空磁导率，

k₁＝cos I₀ sin A₀ sin I+sin I₀ cos I sin A，k₂＝cos I₀ sin A₀ sin I+sinI₀ cos I cos A，

k₃＝cos I₀ cos A₀ cos I sin A+cos I₀ sin A₀ cos I cos A，

k₄＝cos I₀ cos A₀ cos I cos A，

k₅＝cos I₀ sin A₀ cos I sin A，

k₆＝-sin I₀ sin I；

I₀，A₀，I，A分别为x方向相对于总场及地磁场磁偏角和磁倾角， 分别为第j个小立方体的x，y，z方向上立方体边界值；

S4.计算交叉梯度偏导数和拉格朗日算子：

将重磁交叉梯度偏导数利用一阶泰勒进行展开，其交叉梯度项表示拉格朗日算子t(σ，κ)为：

式(7)中，σ，κ分别代表重力密度参数和磁化率参数，B^g、B^m分别为密度参数及磁化率参数所对应的交叉梯度偏导数，且：

该偏导矩阵均表示为：

M为模型元素的个数，在二维情况下，令ii＝(k-1)·Nx·Ny+(j-1)·Nx+i，其中Nx，Ny，Nz分别为x，y，z方向上的单元个数，根据交叉梯度的定义式可得交叉梯度偏导数B^g：

B(ii，ii+Nx)＝-B(ii，ii-Nx)

B(ii，ii+NxNy)＝-B(ii，ii-NxNy)

(10)

同理求出交叉梯度偏导数B^m；

S5.计算获得稀疏信号数据空间迭代模型：

基于光滑l⁰算法对信号数据稀疏空间进行计算，选择期望值为零的高斯函数近似l⁰范数的光滑函数逼近不连续的l⁰范数；设σ为描述连续函数与不连续l⁰范数的近似程度的参数；

令连续函数表示为：

则有：

然后，定义：

则：

使用σ的递减序列，从而使F_σ最大化；前部数值较大的每个的σ，F_σ的最大化算法的初始值为对应的F_σ的最大值；当σ逐渐减小时，每个σ对应的F_σ的初始值从接近实际F_σ的最大值开始；得出最小l⁰范数的解；

引入在稀疏信号空间中紧支撑聚焦反演：

进行加权设定m^ω＝W_eW_mm，d^ω＝W_dd，其中W_e由光滑l⁰算法获得，将(2)(3)式形式表示为：

式中P^α(m^ω，d^ω)代表加权后的密度或磁化率反演方法的目标函数，A^ω为加权后的正演响应，m^ω为加权后的密度或磁化率模型参数，为加权后的先验模型，d^ω为加权后的观测数据矩阵。

将单一方法的变形式转变为：

进而基于交叉梯度的联合反演的变分形式变为：

其中F为A(m)的Frechet偏导数矩阵；

(15)-(17)式使用共轭梯度法进行优化计算，设定：

令最速下降方向为：

式(18)-(20)中，代表加权后的密度或磁化率模型在迭代过程中正演计算响应与数据模型的差值，代表权后的密度或磁化率模型与对应的先验模型的差值，代表模型迭代过程中的最速下降方向。

进而共轭梯度最速下降方向为：

其中且当n＝0，

共轭梯度搜索方向为：

其中E为单位矩阵。

进而迭代模型可表示为：

其中

S6.对获得的稀疏信号数据空间迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出该数据空间联合反演模型，否则执行S3；

其中，拟合差通过式(27)计算获得：

其中，N为观测数据个数。

与现有技术相比，本发明具有以下有益效果：

本发明提供的基于信号数据稀疏空间快速计算的联合反演方法能够将野外采集到的观测数据通过反演计算转换成地下二维、三维密度和磁化率模型的结构信息，从不同物性方面反映更为真实的地下结构信息。该反演方法克服了传统模型空间联合反演方法中，计算速度慢、占用内存大等问题，将模型空间存储和计算灵敏度矩阵通过稀疏空间转换处理进行存储和计算，由于观测数据的个数远远小于模型参数的个数，因此在内存存储量方面，数据空间方法要明显优于模型空间方法，从而提高了计算速度和效率。该反演方法在保证计算精度的条件下，计算效率很高，适用于目前国内外十分主流的windows等操作系统，可以在普通个人电脑上推广使用。该反演方法能够为探索成矿区带不同深度尺度上地球物理特征与成矿之间的关系，为解决深部资源勘查提供有力的依据，指明了找矿、水、石油、地热等资源的方向。

附图说明

图1为本发明的流程图。

图2为验证本发明提供的方法而设计的二维理论模型及计算结果图。图2(i)为模型设计及正演响应计算结果。

图3为验证本发明提供的方法而设计的不同方法对比的三维理论模型及计算结果图。

具体实施方式

下面结合附图说明和实施例对本发明作进一步说明，本发明的方式包括但不仅限于以下实施例。

实施例

如图1所示，本发明公开的一种基于信号数据稀疏空间快速计算的联合反演方法，其特征在于，包括如下步骤：

S1.建立重磁联合反演目标函数，所述联合反演目标函数如式(1)：

约束条件τ＝0；

(1)式中：

其中，角标g，m分别为密度、磁化率有关参数及模型，为对应的目标函数，W_d为观测数据协方差矩阵，A为正演响应，m为密度或磁化率模型参数，d为观测数据，α为正则化参数，为基于信号数据稀疏先验加权矩阵，W_m为模型参数的协方差矩阵，m_apr为先验模型；T为矩阵的转置，为梯度，τ为密度模型及磁化率模型的交叉梯度矩阵；

S2.建立密度、磁化率初始模型：

在笛卡尔坐标系下，二维数据处理中，沿x，z两个坐标轴分别将初始模型空间划分成为N_x、N_z个小的矩形网格单元，间距为：Δx_i(i＝1，2，…N_x)，Δz_i(i＝1，2，…N_z)；三维情况下将初始模型空间划分成为N_x、N_y、N_z个小的矩形网格单元，间距为：Δx_i(i＝1，2，…N_x)，Δy(i＝1，2，…N_y)，Δz_i(i＝1，2，…N_z)；

S3.对密度、磁化率进行正演计算：

采用Pick(1973)及Fu(2015)提出的立方体网格剖分下的密度及磁化率正演公式，即密度正演响应矩阵A_g的元素的计算公式为：

为A_g中的元素，G为万有引力常数，x，y，z为观测点的坐标，为第j个小立方体的x方向上的边界坐标，为第j个小立方体的y方向上的边界坐标，为第j个小立方体的z方向上的边界坐标；

磁化率正演响应矩阵A_mag的元素的计算公式为：

其中μ₀真空磁导率，

k₁＝cos I₀ sin A₀ sin I+sin I₀ cos I sin A，k₂＝cos I₀ sin A₀ sin I+sinI₀ cos I cos A，

k₃＝cos I₀ cos A₀ cos I sin A+cos I₀ sin A₀ cos I cos A，

k₄＝cos I₀ cos A₀ cos I cos A，

k₅＝cos I₀ sin A₀ cos I sin A，

k₆＝-sin I₀ sin I；

I₀，A₀，I，A分别为x方向相对于总场及地磁场磁偏角和磁倾角， 分别为第j个小立方体的x，y，z方向上立方体边界值；

S4.计算交叉梯度偏导数和拉格朗日算子：

将重磁交叉梯度偏导数利用一阶泰勒进行展开，其交叉梯度项表示拉格朗日算子t(σ，κ)为：

式(7)中，σ，κ分别代表重力密度参数和磁化率参数，B^g、B^m分别为密度参数及磁化率参数所对应的交叉梯度偏导数，且：

该偏导矩阵均表示为：

M为模型元素的个数，在二维情况下，令ii＝(k-1)·Nx·Ny+(j-1)·Nx+i，其中Nx，Ny，Nz分别为x，y，z方向上的单元个数，根据交叉梯度的定义式可得交叉梯度偏导数B^g：

B(ii，ii+Nx)＝-B(ii，ii-Nx)

B(ii，ii+NxNy)＝-B(ii，ii-NxNy)

(10)

同理求出交叉梯度偏导数B^m；

S5.计算获得稀疏信号数据空间迭代模型：

基于光滑l⁰算法对信号数据稀疏空间进行计算，选择期望值为零的高斯函数近似l⁰范数的光滑函数逼近不连续的l⁰范数；设σ为描述连续函数与不连续l⁰范数的近似程度的参数；

令连续函数表示为：

则有：

然后，定义：

则：

使用σ的递减序列，从而使F_σ最大化；前部数值较大的每个的σ，F_σ的最大化算法的初始值为对应的F_σ的最大值；当σ逐渐减小时，每个σ对应的F_σ的初始值从接近实际F_σ的最大值开始；得出最小l⁰范数的解；

引入在稀疏信号空间中紧支撑聚焦反演：

进行加权设定m^ω＝W_eW_mm，d^ω＝W_dd，其中W_e由光滑l⁰算法获得，将(2)(3)式形式表示为：

式中P^α(m^ω，d^ω)代表加权后的密度或磁化率反演方法的目标函数，A^ω为加权后的正演响应，m^ω为加权后的密度或磁化率模型参数，为加权后的先验模型，d^ω为加权后的观测数据矩阵。

将单一方法的变形式转变为：

进而基于交叉梯度的联合反演的变分形式变为：

其中F为A(m)的Frechet偏导数矩阵；

(15)-(17)式使用共轭梯度法进行优化计算，设定：

令最速下降方向为：

式(18)-(20)中，代表加权后的密度或磁化率模型在迭代过程中正演计算响应与数据模型的差值，代表权后的密度或磁化率模型与对应的先验模型的差值，代表模型迭代过程中的最速下降方向。

进而共轭梯度最速下降方向为：

其中且当n＝0，

共轭梯度搜索方向为：

其中E为单位矩阵。

进而迭代模型可表示为：

其中

S6.对获得的稀疏信号数据空间迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出该数据空间联合反演模型，否则执行S3；

其中，拟合差通过式(27)计算获得：

其中，N为观测数据个数。

下面从不同角度对本发明提供的基于信号数据稀疏空间重磁联合反演方法进行验证。

图2展示了本方法联合反演结果同单一反演方法的对比结果。如图2(i)所示，设计一个物性结构一致的二维双棱柱模型，左侧异常体剩余密度为1g/cm³，磁化率为1×10-⁴SI异常体大小为30×60m。右侧异常体剩余密度为2g/cm³，磁化率为-1×10-⁴SI异常体大小为30×60m大小为两个异常体顶面埋深均为100m，埋藏于剩余密度及磁化率均假定为0。

图(a)至图(d)为反演结果的对比图。其中(a)(c)分别为未进行联合反演获得的密度、磁化率反演结果，(b)(d)分别为联合反演后获得的密度、磁化率反演结果通过图可以看出基于数据空间快速计算的联合反演方法，异常体的空间几何形态均有更好的恢复，且能较好的反映处所设定的异常体的位置。

图3展示了基于稀疏空间的联合反演方法同其它联合反演方法的对比结果。如图3(a)所示，设计一个物性结构一致的三维双棱柱模型，左右侧异常体剩余密度均为1g/cm³，磁化率均为1×10^-3 SI异常体大小为200×200×100m。两个异常体顶面埋深均为100m，埋藏于剩余密度及磁化率均假定为0。

图(b)至图(g)为采用不同方法进行联合反演结果的对比图，包括光滑反演，常规紧支撑反演及基于信号数据稀疏空间快速计算的联合反演。通过图可以看出基于信号数据稀疏空间快速计算的联合反演，异常体的空间几何形态均有更好的恢复，且能较好的反映处所设定的异常体的位置。重、磁联合反演结果，异常体上下界面清晰可见，这是由于交叉梯度结构相似性的约束。同时，基于信号数据稀疏空间快速计算的联合反演相对于其它方法反演结果恢复效果更好。可以看出，光滑空间反演结果较为光滑，异常体边缘相对模糊，而数据空间反演结果聚焦程度更高，异常体边界清晰，异常体几何形态和真实模型基本吻合。在反演计算时间、迭代次数方面，基于信号数据稀疏空间快速计算的联合反演在内存消耗方面，RMS更低，收敛次数更小。可以看出，地下不同物性模型结构不一致时，交叉梯度联合反演也可以得到准确的反演结果，同时无论在计算时间还是在内存消耗方面数据空间联合反演更具有优势。

本发明通过上述设计，有效地解决了解决现有的联合反演方法计算速度慢、占用内存大等问题，提高了计算速度和效率。为重磁数据处理与解释提供了一种行之有效的计算方法。

上述实施例仅为本发明的优选实施方式之一，不应当用于限制本发明的保护范围，但凡在本发明的主体设计思想和精神上做出的毫无实质意义的改动或润色，其所解决的技术问题仍然与本发明一致的，均应当包含在本发明的保护范围之内。

Claims (2)

1.一种基于信号数据稀疏空间快速计算的联合反演方法，其特征在于，包括如下步骤：

S1.建立重磁联合反演目标函数，所述联合反演目标函数如式(1)：

约束条件τ＝0；

(1)式中：

其中，角标g，m分别为密度、磁化率有关参数及模型，为对应的目标函数，W_d为观测数据协方差矩阵，A为正演响应，m为密度或磁化率模型参数，d为观测数据，α为正则化参数，为基于信号数据稀疏先验加权矩阵，W_m为模型参数的协方差矩阵，m_apr为先验模型；T为矩阵的转置，为梯度，τ为密度模型及磁化率模型的交叉梯度矩阵；

S2.建立密度、磁化率初始模型：

在笛卡尔坐标系下，二维数据处理中，沿x，z两个坐标轴分别将初始模型空间划分成为N_x、N_z个小的矩形网格单元，间距为：Δx_i(i＝1，2，…N_x)，Δz_i(i＝1，2，…N_z)；三维情况下将初始模型空间划分成为N_x、N_y、N_z个小的矩形网格单元，间距为：Δx_i(i＝1，2，…N_x)，Δy(i＝1，2，…N_y)，Δz_i(i＝1，2，…N_z)；

S3.对密度、磁化率进行正演计算：

采用Pick(1973)及Fu(2015)提出的立方体网格剖分下的密度及磁化率正演公式，即密度正演响应矩阵A_g的元素的计算公式为：

为A_g中的元素，G为万有引力常数，x，y，z为观测点的坐标，为第j个小立方体的x方向上的边界坐标，为第j个小立方体的y方向上的边界坐标，为第j个小立方体的z方向上的边界坐标；

磁化率正演响应矩阵A_mag的元素的计算公式为：

其中μ₀真空磁导率，

k₁＝cosI₀sinA₀sinI+sinI₀cosIsinA，k₂＝cosI₀sinA₀sinI+sinI₀cosIcosA，

k₃＝cosI₀cosA₀cosIsinA+cosI₀sinA₀cosIcosA，

k₄＝cosI₀cosA₀cosIcosA，

k₅＝cosI₀sinA₀cosIsinA，

k₆＝-sinI₀sinI；

I₀，A₀，I，A分别为x方向相对于总场及地磁场磁偏角和磁倾角， 分别为第j个小立方体的x，y，z方向上立方体边界值；

S4.计算交叉梯度偏导数和拉格朗日算子：

将重磁交叉梯度偏导数利用一阶泰勒进行展开，其交叉梯度项表示拉格朗日算子t(σ，κ)为：

式(7)中，σ，κ分别代表重力密度参数和磁化率参数，B^g、B^m分别为密度参数及磁化率参数所对应的交叉梯度偏导数，且：

该偏导矩阵均表示为：

M为模型元素的个数，在二维情况下，令ii＝(k-1)·Nx·Ny+(j-1)·Nx+i，其中Nx，Ny，Nz分别为x，y，z方向上的单元个数，根据交叉梯度的定义式可得交叉梯度偏导数B^g：

B(ii，ii+Nx)＝-B(ii，ii-Nx)

B(ii，ii+NxNy)＝-B(ii，ii-NxNy)

(10)

同理求出交叉梯度偏导数B^m；

S5.计算获得稀疏信号数据空间迭代模型：

基于光滑l⁰算法对信号数据稀疏空间进行计算，选择期望值为零的高斯函数近似l⁰范数的光滑函数逼近不连续的l⁰范数；设σ为描述连续函数与不连续l⁰范数的近似程度的参数；

令连续函数表示为：

则有：

然后，定义：

则：

使用σ的递减序列，从而使F_σ最大化；前部数值较大的每个的σ，F_σ的最大化算法的初始值为对应的F_σ的最大值；当σ逐渐减小时，每个σ对应的F_σ的初始值从接近实际F_σ的最大值开始；得出最小l⁰范数的解；

引入在稀疏信号空间中紧支撑聚焦反演：

进行加权设定m^ω＝W_eW_mm，d^ω＝W_dd，其中W_e由光滑l⁰算法获得，将(2)(3)式形式表示为：

式中P^α(m^ω，d^ω)代表加权后的密度或磁化率反演方法的目标函数，A^ω为加权后的正演响应，m^ω为加权后的密度或磁化率模型参数，为加权后的先验模型，d^ω为加权后的观测数据矩阵。

将单一方法的变形式转变为：

进而基于交叉梯度的联合反演的变分形式变为：

其中F为A(m)Frechet偏导数矩阵；

(15)-(17)式使用共轭梯度法进行优化计算，设定：

令最速下降方向为：

式(18)-(20)中，代表加权后的密度或磁化率模型在迭代过程中正演计算响应与数据模型的差值，代表权后的密度或磁化率模型与对应的先验模型的差值，代表模型迭代过程中的最速下降方向；

进而共轭梯度最速下降方向为：

其中且当n＝0，

共轭梯度搜索方向为：

其中E为单位矩阵；

进而迭代模型可表示为：

其中

S6.对获得的稀疏信号数据空间迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出该数据空间联合反演模型，否则执行S3；

其中，拟合差通过式(27)计算获得：

其中，N为观测数据个数。

2.根据权利要求1所述的一种基于信号数据稀疏空间快速计算的联合反演方法，其特征在于，所述步骤S4中，如果在三维情况下，先对t(i，j，k)元素求偏导值，然后进行运算。