[发明专利]一种基于信号数据稀疏空间快速计算的联合反演方法

摘要

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

主权项

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为观测数据个数。