CN112649848A - 利用波动方程求解地震波阻抗的方法和装置 - Google Patents

Abstract

本申请公开了利用波动方程求解地震波阻抗的方法和装置。该方法包括：根据波动方程建立目标函数φ(m)，其中，m为由傅里叶级数展开形式表示的波阻抗模型的傅里叶系数构成的矩阵；基于所述目标函数确定矩阵m中各个元素的值；将矩阵m中各个元素的值代入所述波阻抗模型。根据本申请，可以实现超深层少井或未钻进至目标层以及强非均质性储层高精度宽带波阻抗反演及应用，使得有井段与常规反演效果在波动方程反演的基础上对于局部边缘特征刻画分辨率更高，使得无井段在保持局部特征的同时提高深层的全局连续性，可为特深层断溶型储层的刻画及早期储量评估提供数据支撑。

Description

利用波动方程求解地震波阻抗的方法和装置

技术领域

本发明属于地震勘探领域，更具体地，涉及一种利用波动方程求解地震波阻抗的方法和一种利用波动方程求解地震波阻抗的装置。

背景技术

波阻抗信息作为联系地质和地球物理的一座桥梁，在地震解释中有重要作用。地震反演得到波阻抗，无论是叠前波阻抗反演还是叠后波阻抗反演，反演方法一般以基于褶积模型的方法为主。常用算法大体可分为以反褶积法为基础的直接反演(如递归反演，道积分等)、以模型为基础的反演和约束反演等，但褶积模型在平层假设的前提下，对于横向分辨率的提高存在局限性；其次常规方法缺乏超深层无井及少井情况下储层特征描述能力。对于波形反演类的方法在反演弹性参数时的主要问题在于二维及以上波动方程受多参数解耦及反演策略的限制，同时反演多个参数才求解多参数的稳定性和计算效率收到限制，函数的非凸性容易使求解过程陷入局部极小。

现有反演波阻抗的方法多是基于褶积模型，反演结果依赖于钻井记录构建的初始波阻抗模型。在数据采集不完整、少井或者钻井深度不到目标层的地区初始波阻抗模型不易构建，难以得到稳定可靠的波阻抗模型，不利于深部储层的解释与评价。

基于波动方程的波阻抗反演方法包括三大类：

(1)间接的从波动方程反演得到波阻抗参数,此方法意味着可能要同时反演多个参数，利用反演得到的其他参数间接得到波阻抗参数。例如，从波动方程中同时反演速度和密度参数，利用二者的乘积得到波阻抗模型；或者同时反演拉梅常数和密度参数，再计算波阻抗模型，Shi等(2007)同时反演了这三个参数，并将这些参数转化成泊松比等参数来寻找某地区数据的含气砂层；或者将方程转化成包含波阻抗参数和密度参数的方程，通过同时反演波阻抗参数和密度参数的方式得到波阻抗模型，Tatantola(1986)给出了同时反演P波阻抗、S波阻抗和密度三个参数的理论公式；等等。以上几种方法在反演过程中都涉及同时反演两个或两个以上参数的问题，考虑到多参数同时反演时目标函数对不同参数的敏感性不同、参数之间的耦合作用等因素的影响，多参数同时反演的难度较大，不易同时反演得到满足要求的多个参数，由这些参数转化得到的波阻抗参数的准确性可能得不到保证。

(2)有一种可以直接得到波阻抗参数的方法是反演波阻抗参数时假设速度参数已知，Plessix和Li(2013)利用谱分形的波形阻抗反演方法得到了声阻抗模型，该方法要求背景速度满足可以准确模拟反射的动力学特征的假设，即背景速度要求包含准确的低波数信息。

(3)另外一种可以避开多参数同时反演、直接反演得到波阻抗参数的方法是将原含有多个物性参数的方程转换为只含有一个物性参数(波阻抗)的波动方程，Santosa和Schwetlick(1982)将一维剪切波波动方程转换成只用剪切波阻抗表示的波动方程，从这个转换后的方程中可以直接反演阻抗参数。Santosa和Schwetlick(1982)用特征法反演了阻抗参数，李勤学等(1994)基于此方程应用逐段折叠方法得到了阻抗参数，张丽琴等(2004)利用同伦方法直接反演了波阻抗参数。

除以上几种反演方法之外，全波形反演方法也可以用来从此方程中反演波阻抗参数，目前还未看到相关研究。用全波形反演方法反演波阻抗参数时，可以选用时间域全波形反演方法(Tarantola，1984；Tarantola，1986；Gauthier等，1986；

2011)或者频率域全波形反演方法(Pratt&Worthington，1990；Pratt等，1998；Pratt，1999；Ravaut等，2004；Wang&Rao，2006；Wang&Rao，2009)。

相比基于褶积模型的波阻抗反演，由于波动方程更准确地描述了地震波的传播理论，所以基于波动方程的波阻抗反演可得到精度更高的波阻抗模型。但现有的各种基于波动方程的反演方法仍存在各自的缺陷，特别地，对于二维或三维数据体，利用现有基于一维波动方程反演得到的二维或三维波阻抗剖面可能会出现横向不连续。

发明内容

有鉴于此，本申请提出了一种能解决多维波阻抗剖面横向不连续问题的波阻抗反演方案。

根据本发明的一个方面，提供了一种利用波动方程求解地震波阻抗的方法，所述方法包括：根据波动方程建立目标函数φ(m)：

其中，m为由傅里叶级数展开形式表示的波阻抗模型的傅里叶系数构成的矩阵，

是坐标(pΔx,qΔy)处零偏移距的观测地震数据，

是所述波阻抗模型上坐标(pΔx,qΔy)处的合成地震数据，Δx、Δy和Δt分别表示在三个坐标维度下的采样步长，P、Q和H对应表示在三个坐标维度下的最大采样点数；

基于所述目标函数确定矩阵m中各个元素的值；

将矩阵m中各个元素的值代入所述波阻抗模型。

在一种可能的实施方式中，由傅里叶级数展开形式表示的波阻抗模型如下所示：

其中，I(x,y,t)是波阻抗参数，{x,y,t}分别是在三个坐标维度下的取值，{L,M,N}分别是三个坐标维度上谐波项的最大数，{Δk,Δω}是空间域的基本波数，(a_lmn,b_lmn,...,h_lmn)构成矩阵m。

在一种可能的实施方式中，基于所述目标函数确定矩阵m中各个元素的值，包括：

S21，计算目标函数φ(m)对每个傅里叶系数的一阶偏导

coef_i为矩阵m中第i个元素，i＝1,2,...,R，R为矩阵m中的元素个数；

S22，基于下式更新矩阵m_k，m_k为开始第k次迭代时矩阵m的当前取值：

m_k+1＝m_k+αΔm_k，其中，Δm_k为负梯度矩阵，Δm_k,i为负梯度矩阵Δm_k中的第i个元素，

α为步长，m₁为已给出的初始傅里叶系数矩阵；

S23，判断是否满足迭代终止条件，如果满足，则令矩阵m＝m_k+1，并结束迭代；如果不满足，则令k＝k+1，并回到步骤S21，进入下次迭代。

在一种可能的实施方式中，根据下式计算一阶偏导

其中，[Δu^pq]是所述观测地震数据与所述合成地震数据的数据残差，

是波场对傅里叶系数coef_i的Frechet导数。

在一种可能的实施方式中，所述迭代终止条件包括：

目标函数的值小于预设阈值；

或者迭代次数达到最大迭代次数。

根据本申请的另一方面，提供了一种利用波动方程求解地震波阻抗的装置，所述装置包括：

目标函数建立单元，用于根据波动方程建立目标函数φ(m)：

其中，m为由傅里叶级数展开形式表示的波阻抗模型的傅里叶系数构成的矩阵，

是坐标(pΔx,qΔy)处零偏移距的观测地震数据，

是所述波阻抗模型上坐标(pΔx,qΔy)处的合成地震数据，Δx、Δy和Δt分别表示在三个坐标维度下的采样步长，P、Q和H对应表示在三个坐标维度下的最大采样点数；

傅里叶系数确定单元，用于基于所述目标函数确定矩阵m中各个元素的值；

系数代入单元件，用于将矩阵m中各个元素的值代入所述波阻抗模型。

在一种可能的实施方式中，由傅里叶级数展开形式表示的波阻抗模型如下所示：

其中，I(x,y,t)是波阻抗参数，{x,y,t}分别是在三个坐标维度下的取值，{L,M,N}分别是三个坐标维度上谐波项的最大数，{Δk,Δω}是空间域的基本波数，(a_lmn,b_lmn,...,h_lmn)构成矩阵m。

在一种可能的实施方式中，在傅里叶系数确定单元中，基于所述目标函数确定矩阵m中各个元素的值，包括：

S21，计算目标函数φ(m)对每个傅里叶系数的一阶偏导

coef_i为矩阵m中第i个元素，i＝1,2,...,R，R为矩阵m中的元素个数；

S22，基于下式更新矩阵m_k，m_k为开始第k次迭代时矩阵m的当前取值：

m_k+1＝m_k+αΔm_k，其中，Δm_k为负梯度矩阵，Δm_k,i为负梯度矩阵Δm_k中的第i个元素，

α为步长，m₁为已给出的初始傅里叶系数矩阵；

S23，判断是否满足迭代终止条件，如果满足，则令矩阵m＝m_k+1，并结束迭代；如果不满足，则令k＝k+1，并回到步骤S21，进入下次迭代。

在一种可能的实施方式中，根据下式计算一阶偏导

其中，[Δu^pq]是所述观测地震数据与所述合成地震数据的数据残差，

是波场对傅里叶系数coef_i的Frechet导数。

在一种可能的实施方式中，所述迭代终止条件包括：

目标函数的值小于预设阈值；

或者迭代次数达到最大迭代次数。

根据本申请，通过计算波阻抗级数展开的傅里叶系数可有效提高反演结果全空间的横向分辨率；同时针对超深层少井或未钻进至目标层以及强非均质性储层。常规阻抗反演方法无法得到可靠反演结果，而本申请在推导时间-深度域波动方程正演和全空间波形反演单参数反演波阻抗公式基础上，开展了地震子波估计方法、宽带波阻抗的全空间同步反演方法、及阻抗参数波形反演方法研究。基于本申请，可以实现超深层少井或未钻进至目标层以及强非均质性储层高精度宽带波阻抗反演及应用，使得有井段与常规反演效果在波动方程反演的基础上对于局部边缘特征刻画分辨率更高，使得无井段在保持局部特征的同时提高深层的全局连续性，可为特深层断溶型储层的刻画及早期储量评估提供数据支撑。

附图说明

通过结合附图对本申请示例性实施方式进行更详细的描述，本申请的上述以及其它目的、特征和优势将变得更加明显，其中，在本申请示例性实施方式中，相同的参考标号通常代表相同部件。

图1示出根据本申请的一个实施例的一种利用波动方程求解地震波阻抗的方法的流程图。

图2示出根据本申请的一个实施例的一种利用波动方程求解地震波阻抗的装置的结构框图。

图3示出根据本申请的一个应用示例的示意图。

图4示出根据地震数据估算的地震子波的示意图。

图5A是根据测井数据得到的波阻抗模型；图5B是根据图5A的波阻抗模型得到的由傅里叶系数转换得到的波阻抗模型。

图6示出对地震数据进行均衡后得到的观测地震记录。

图7示出根据本申请得到的全空间波阻抗模型。

图8A示出观测地震记录，图8B示出利用根据本申请得到的波阻抗模型做波长模拟得到的模拟波形，图8C为某固定位置处的观测地震记录和根据本申请得到的合成地震记录的波形对比。

具体实施方式

下面将参照附图更详细地描述本申请的优选实施方式。虽然附图中显示了本申请的优选实施方式，然而应该理解，可以以各种形式实现本申请而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本申请更加透彻和完整，并且能够将本申请的范围完整地传达给本领域的技术人员。

Santosa&Schwetlick(1982)将一维剪切波波动方程转换成只用剪切波阻抗表示的波动方程，采用同样的方式，发明人可以将一维声波方程转换为只用声波阻抗表示的方程。得到由密度和体积模量参数表示的一维波动方程：

式中z是深度，κ(z)＝ρ(z)v²(z)代表体积模量，ρ(z)是体密度,v(z)是速度，u(z,t)是位移，t是旅行时。

令

其中τ是一个以时间为单位的新“深度”坐标，该变换将深度域的物理模型从深度域变换到时间域，则转换后的方程为：

式中Z(τ)是波阻抗，可以看出变换后的方程中只包含波阻抗一个物性参数。可用有限差分方法对上述转换后的方程(2)进行数值求解(Hall&Wang,2009；Levander,1988)。

从方程(2)可以看出，变换后的波动方程是一维方程，直接从该方程中反演波阻抗参数，只能得到一道的波阻抗模型。对于二维或三维数据体，利用一维波动方程逐道反演得到的二维或三维波阻抗剖面可能会出现横向不连续。因此，发明人以波动方程为基础，但不直接反演波阻抗参数，而是通过将波阻抗参数做傅里叶级数展开后反演傅里叶系数，最后用反演得到的傅里叶系数重构二维或三维波阻抗模型。

请参见图1。图1示出根据本申请的一个实施例的一种利用波动方程求解地震波阻抗的方法的流程图。如图所示，所述方法包括下列步骤。

步骤1，根据波动方程建立目标函数φ(m)：

其中，m为由傅里叶级数展开形式表示的波阻抗模型的傅里叶系数构成的矩阵，

是坐标(pΔx,qΔy)处零偏移距的观测地震数据，

是所述波阻抗模型上坐标(pΔx,qΔy)处的合成地震数据，Δx、Δy和Δt分别表示在三个坐标维度下的采样步长，P、Q和H对应表示在三个坐标维度下的最大采样点数。

在一种可能的实施方式中，用傅里叶级数展开形式表示的波阻抗模型如下所示：

其中，I(x,y,t)是波阻抗参数，{x,y,t}分别是在三个坐标维度下的取值，{L,M,N}分别是三个坐标维度上谐波项的最大数，{Δk,Δω}是空间域的基本波数，(a_lmn,b_lmn,...,h_lmn)构成矩阵m。

本领域技术人员可以理解的是，上述波阻抗模型不仅可用于表示三维波阻抗模型，也可用于表示二维模型，此时{x,y,t}中的一个参数可恒定取0，{L,M,N}中的对应参数同样恒定取0。类似地，上述波阻抗模型也可用于表示一维波阻抗模型。

步骤2，基于所述目标函数确定矩阵m中各个元素的值。

在一种可能的实施方式中，基于所述目标函数确定矩阵m中各个元素的值，具体可以包括下列步骤S21到步骤S23。

S21，计算目标函数φ(m)对每个傅里叶系数的一阶偏导

coef_i为矩阵m中第i个元素，i＝1,2,...,R，R为矩阵m中的元素个数。本领域技术人员可以理解的是，如果是三维波阻抗模型，则有R＝8；如果是二维波阻抗模型，则有R＝4。

具体地，可以根据下式计算一阶偏导

其中，[Δu^pq]是所述观测地震数据与所述合成地震数据的数据残差，

是波场对傅里叶系数coef_i的Frechet导数。

S22，基于下式更新矩阵m_k，m_k为开始第k次迭代时矩阵m的当前取值：

m_k+1＝m_k+αΔm_k，其中，Δm_k为负梯度矩阵，Δm_k,i为负梯度矩阵Δm_k中的第i个元素，

α为步长，步长α可以用线性化搜索方法计算(例如Nocedal&Wright,2006)得到，每次迭代后可根据当前结果重新搜索步长α，m₁为已给出的初始傅里叶系数矩阵。

S23，判断是否满足迭代终止条件，如果满足，则令矩阵m＝m_k+1，并结束迭代；如果不满足，则令k＝k+1，并回到步骤S21，进入下次迭代。

在一种可能的实施方式中，所述迭代终止条件包括：

目标函数的值小于预设阈值；

或者迭代次数达到最大迭代次数。

即只要满足上述条件中的任意一者，则可结束迭代。

步骤3，将矩阵m中各个元素的值代入所述波阻抗模型。

上述实施例在推导时间-深度域波动方程正演和全空间波形反演单参数反演波阻抗公式基础上，开展了地震子波估计方法、宽带波阻抗的全空间同步反演方法、及阻抗参数波形反演方法研究。基于上述实施例，可以实现超深层少井或未钻进至目标层以及强非均质性储层高精度宽带波阻抗反演及应用，使得有井段与常规反演效果在波动方程反演的基础上对于局部边缘特征刻画分辨率更高，使得无井段在保持局部特征的同时提高深层的全局连续性，可为特深层断溶型储层的刻画及早期储量评估提供数据支撑。

请参加图2。图2示出根据本申请的一个实施例的一种利用波动方程求解地震波阻抗的装置的结构框图。如图2所示，所述装置包括：目标函数建立单元202、傅里叶系数确定单元204和系数代入单元206。

目标函数建立单元202用于根据波动方程建立目标函数φ(m)：

其中，m为由傅里叶级数展开形式表示的波阻抗模型的傅里叶系数构成的矩阵，

是坐标(pΔx,qΔy)处零偏移距的观测地震数据，

是所述波阻抗模型上坐标(pΔx,qΔy)处的合成地震数据，Δx、Δy和Δt分别表示在三个坐标维度下的采样步长，P、Q和H对应表示在三个坐标维度下的最大采样点数；

傅里叶系数确定单元204用于基于所述目标函数确定矩阵m中各个元素的值；

系数代入单元件206用于将矩阵m中各个元素的值代入所述波阻抗模型。

在一种可能的实施方式中，由傅里叶级数展开形式表示的波阻抗模型如下所示：

其中，I(x,y,t)是波阻抗参数，{x,y,t}分别是在三个坐标维度下的取值，{L,M,N}分别是三个坐标维度上谐波项的最大数，{Δk,Δω}是空间域的基本波数，(a_lmn,b_lmn,...,h_lmn)构成矩阵m。

在一种可能的实施方式中，在傅里叶系数确定单元204中，基于所述目标函数确定矩阵m中各个元素的值，包括：

S21，计算目标函数φ(m)对每个傅里叶系数的一阶偏导

coef_i为矩阵m中第i个元素，i＝1,2,...,R，R为矩阵m中的元素个数；

S22，基于下式更新矩阵m_k，m_k为开始第k次迭代时矩阵m的当前取值：

m_k+1＝m_k+αΔm_k，其中，Δm_k为负梯度矩阵，Δm_k,i为负梯度矩阵Δm_k中的第i个元素，

α为步长，m₁为已给出的初始傅里叶系数矩阵；

S23，判断是否满足迭代终止条件，如果满足，则令矩阵m＝m_k+1，并结束迭代；如果不满足，则令k＝k+1，并回到步骤S21，进入下次迭代。

在一种可能的实施方式中，根据下式计算一阶偏导

其中，[Δu^pq]是所述观测地震数据与所述合成地震数据的数据残差，

是波场对傅里叶系数coef_i的Frechet导数。

在一种可能的实施方式中，所述迭代终止条件包括：

目标函数的值小于预设阈值；

或者迭代次数达到最大迭代次数。。

应用示例

下文以三维数据中的一条测线为例，对本申请的效果进行示例性说明。

图3示出根据本申请的一个应用示例的示意图。如图3所示，可先将地震数据301作为输入，对其进行傅里叶变换，利用功率谱和能量谱等信息估算构成广义子波的两个关键参数，再根据广义子波公式得到估算的地震子波302(如图4所示)。在此过程中可直接从地震数据中估算子波，而无需测井信息。

还可根据地震数据301得到初始波阻抗模型303(如图5A所示)。可基于地震子波302和初始波阻抗模型303，得到初始合成地震数据304。然后可对比初始合成地震数据304对地震数据301做均衡处理，得到观测地震数据305(如图6所示)，使其的振幅量级与初始合成地震数据304相近。

此外，可基于波阻抗模型303进行傅里叶级数展开，得到初始傅里叶系数矩阵306。图5B示出了基于初始傅里叶系数矩阵306构建的初始傅里叶级数波阻抗模型。然后根据本申请所公开的技术方案，得到由傅里叶级数展开形式表示的波阻抗模型307(如图7所示)。

利用波阻抗模型307和地震子波302做波长反演，得到合成地震数据308(如图8B所示)。图8A为观测地震记录。图8C为某固定位置处的观测地震数据305和根据本申请得到的合成地震数据308的波形对比，即反演结果评价309，可以看出两者差异很小，充分验证了本申请的有效性。

本申请可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本申请的各个方面的计算机可读程序指令。

计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是――但不限于――电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波(例如，通过光纤电缆的光脉冲)、或者通过电线传输的电信号。

这里参照根据本申请实施例的方法、装置(系统)和计算机程序产品的流程图和/或框图描述了本申请的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。

以上已经描述了本申请的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。

Claims (10)

1.一种利用波动方程求解地震波阻抗的方法，其特征在于，所述方法包括：

根据波动方程建立目标函数φ(m)：

其中，m为由傅里叶级数展开形式表示的波阻抗模型的傅里叶系数构成的矩阵，

是坐标(pΔx,qΔy)处零偏移距的观测地震数据，

是所述波阻抗模型上坐标(pΔx,qΔy)处的合成地震数据，Δx、Δy和Δt分别表示在三个坐标维度下的采样步长，P、Q和H对应表示在三个坐标维度下的最大采样点数；

基于所述目标函数确定矩阵m中各个元素的值；

将矩阵m中各个元素的值代入所述波阻抗模型。

2.根据权利要求1所述的方法，其特征在于，由傅里叶级数展开形式表示的波阻抗模型如下所示：

其中，I(x,y,t)是波阻抗参数，{x,y,t}分别是在三个坐标维度下的取值，{L,M,N}分别是三个坐标维度上谐波项的最大数，{Δk,Δω}是空间域的基本波数，(a_lmn,b_lmn,...,h_lmn)构成矩阵m。

3.根据权利要求1所述的方法，其特征在于，基于所述目标函数确定矩阵m中各个元素的值，包括：

S21，计算目标函数φ(m)对每个傅里叶系数的一阶偏导

coef_i为矩阵m中第i个元素，i＝1,2,...,R，R为矩阵m中的元素个数；

S22，基于下式更新矩阵m_k，m_k为开始第k次迭代时矩阵m的当前取值：

m_k+1＝m_k+αΔm_k，其中，Δm_k为负梯度矩阵，Δm_k,i为负梯度矩阵Δm_k中的第i个元素，

α为步长，m₁为已给出的初始傅里叶系数矩阵；

S23，判断是否满足迭代终止条件，如果满足，则令矩阵m＝m_k+1，并结束迭代；如果不满足，则令k＝k+1，并回到步骤S21，进入下次迭代。

4.根据权利要求3所述的方法，其特征在于，根据下式计算一阶偏导

其中，[Δu^pq]是所述观测地震数据与所述合成地震数据的数据残差，

是波场对傅里叶系数coef_i的Frechet导数。

5.根据权利要求3所述的方法，其特征在于，所述迭代终止条件包括：

目标函数的值小于预设阈值；

或者迭代次数达到最大迭代次数。

6.一种利用波动方程求解地震波阻抗的装置，其特征在于，所述装置包括：

目标函数建立单元，用于根据波动方程建立目标函数φ(m)：

其中，m为由傅里叶级数展开形式表示的波阻抗模型的傅里叶系数构成的矩阵，

是坐标(pΔx,qΔy)处零偏移距的观测地震数据，

是所述波阻抗模型上坐标(pΔx,qΔy)处的合成地震数据，Δx、Δy和Δt分别表示在三个坐标维度下的采样步长，P、Q和H对应表示在三个坐标维度下的最大采样点数；

傅里叶系数确定单元，用于基于所述目标函数确定矩阵m中各个元素的值；

系数代入单元，用于将矩阵m中各个元素的值代入所述波阻抗模型。

7.根据权利要求6所述的装置，其特征在于，由傅里叶级数展开形式表示的波阻抗模型如下所示：

其中，I(x,y,t)是波阻抗参数，{x,y,t}分别是在三个坐标维度下的取值，{L,M,N}分别是三个坐标维度上谐波项的最大数，{Δk,Δω}是空间域的基本波数，(a_lmn,b_lmn,...,h_lmn)构成矩阵m。

8.根据权利要求6所述的装置，其特征在于，在傅里叶系数确定单元中，基于所述目标函数确定矩阵m中各个元素的值，包括：

S21，计算目标函数φ(m)对每个傅里叶系数的一阶偏导

coef_i为矩阵m中第i个元素，i＝1,2,...,R，R为矩阵m中的元素个数；

S22，基于下式更新矩阵m_k，m_k为开始第k次迭代时矩阵m的当前取值：

m_k+1＝m_k+αΔm_k，其中，Δm_k为负梯度矩阵，Δm_k,i为负梯度矩阵Δm_k中的第i个元素，

α为步长，m₁为已给出的初始傅里叶系数矩阵；

S23，判断是否满足迭代终止条件，如果满足，则令矩阵m＝m_k+1，并结束迭代；如果不满足，则令k＝k+1，并回到步骤S21，进入下次迭代。

9.根据权利要求8所述的装置，其特征在于，根据下式计算一阶偏导

其中，[Δu^pq]是所述观测地震数据与所述合成地震数据的数据残差，

是波场对傅里叶系数coef_i的Frechet导数。

10.根据权利要求8所述的装置，其特征在于，所述迭代终止条件包括：

目标函数的值小于预设阈值；

或者迭代次数达到最大迭代次数。

Images