CN114063160B - 一种地震速度反演方法及装置 - Google Patents

Abstract

本发明提供了一种地震速度反演方法及装置，属于油气物探工程领域。该方法包括：结合最大时移量、迭代总次数和组合叠加的个数得到第k次迭代的组合时移量；结合第k次迭代的合成极性编码超级炮和组合时移量及每个时移量对应的加权系数得到第k次迭代的组合编码超级炮

利用第k次迭代的组合时移量及每个时移量对应的加权系数、地震子波和极性编码矩阵得到第k次迭代的组合震源

并利用

和第k‑1次迭代的速度场得到第k次迭代的正演超级炮D_cal,k(v)；利用D_cal,k(v)和

得到第k次迭代的反演梯度场；结合第k次迭代的更新步长和反演梯度场以及第k‑1次迭代的速度场得到第k次迭代的速度场，循环迭代得到最终的速度场。本发明能增加反演过程的稳定性、提高反演精度。

Description

一种地震速度反演方法及装置

技术领域

本发明涉及一种地震速度反演方法及装置，属于油气物探工程领域。

背景技术

在当前的勘探领域中，地震速度场的反演对于地震数据的偏移成像和解释具有重要的意义，全波形反演能够充分利用地震数据中所包含的振幅、相位、走时、波形等信息，对地下模型参数进行精细刻画。

但是全波形反演存在计算量巨大、反演稳定性等诸多问题，由于采集成本、坏道、噪声、地形等诸多因素的影响，利用实际采集到的数据提取的子波往往不准确，并且初始速度模型与地下的真实模型差异较大，导致反演过程不稳定，给地震数据的速度反演造成了一定的困扰，而且巨大的计算量一直是限制FWI实用化的关键问题之一。

申请公布号为CN109541691A的发明专利申请文件中公开了一种地震速度反演方法，该方法对常规的维纳滤波全波形反演方法进行改进，对维纳滤波采用的参考道是一种加权的多震源地震道，不需要求取原始子波，能够减少子波估计的计算量，提高计算效率。虽然该方法能够提高地震速度反演的计算效率，但是在初始速度场与地下的真实模型相差较大的情况下，无法保证地震反演速度场的稳定性和准确性。

发明内容

本发明的目的在于提供一种地震速度反演方法及装置，用以解决在初始速度场与地下的真实模型相差较大的情况下，利用目前的地震速度反演方法获得的地震反演速度场的稳定性和准确性较差的问题。

为实现上述目的，本发明提供了一种地震速度反演方法，该方法包括以下步骤：

获取初始观测地震数据、初始速度场和极性编码矩阵；所述极性编码矩阵是一个K×N_s的矩阵，矩阵中的值为随机生成的正负1，K为迭代总次数，N_s为编码的总炮数；

根据所述初始观测地震数据、初始速度场和极性编码矩阵迭代求取速度场，直至达到迭代总次数或者生成的速度场满足设定要求，得到最终的速度场；

其中，第k次迭代的速度场通过以下步骤得到：

结合最大时移量、迭代总次数和组合叠加的个数得到第k次迭代的组合时移量，1≤k≤K；所述最大时移量利用所述初始观测地震数据的主频计算得到，所述组合时移量中包含的时移量个数等于组合叠加的个数；

利用所述极性编码矩阵对所述初始观测地震数据进行编码，得到第k次迭代的合成极性编码超级炮，并利用第k次迭代的组合时移量及组合时移量中每个时移量对应的加权系数和第k次迭代的合成极性编码超级炮，得到第k次迭代的组合编码超级炮；

利用第k次迭代的组合时移量及组合时移量中每个时移量对应的加权系数、主频与所述初始观测地震数据的主频相等的地震子波和极性编码矩阵，得到第k次迭代的组合震源，并利用第k次迭代的组合震源和第k-1次迭代的速度场得到第k次迭代的正演超级炮；其中，k＝1时，第k-1次迭代的速度场为所述初始速度场；

利用第k次迭代的正演超级炮与第k次迭代的组合编码超级炮做差反传，互相关得到第k次迭代的反演梯度场；

结合第k次迭代的更新步长、第k次迭代的反演梯度场和第k-1次迭代的速度场得到第k次迭代的速度场。

本发明还提供了一种地震速度反演装置，该装置包括处理器和存储器，所述处理器执行由所述存储器存储的计算机程序，以实现上述的地震速度反演方法。

本发明的有益效果是：本发明在地震速度反演过程中引入组合时移量，利用该组合时移量及组合时移量中每个时移量对应的加权系数分别对地震数据和子波进行时移组合叠加得到组合编码超级炮和正演超级炮，然后对正演超级炮与组合编码超级炮做差反传，互相关求取反演梯度场，进而迭代更新求取最终的反演速度场；由于引入了组合时移量，组合编码超级炮的同相轴明显变粗，相比合成极性编码超级炮，地震数据主频明显降低，从而就使得反演过程更加稳定；且由于计算地震反演速度场时使用的炮记录是多个经过时移的炮记录，与现有技术中仅使用某一时刻的炮记录相比，更加充分地利用了炮记录中的信息，而且，每个经过时移的炮记录还进行了加权操作，更能够突出炮记录中的重要信息，从而进一步增加反演过程的稳定性和反演精度；另外，本发明采用极性编码生成超级炮的形式进行反演，计算量大大减少，且能够加大临炮之间的相干性，较好的压制多震源间的串扰噪音，更精确的对复杂区域进行速度场的反演，提高反演精度。

进一步地，在上述方法及装置中，所述生成的速度场满足设定要求是指生成的速度场能使正演超级炮与组合编码超级炮的差小于设定阈值；在每次迭代后均判断本次迭代生成的速度场是否满足设定要求，若满足设定要求，则本次迭代的速度场即为最终的速度场；若达到迭代总次数时生成的速度场仍不满足设定要求，则从历次迭代得到的速度场中选取一个能使正演超级炮与组合编码超级炮的差最小的速度场作为最终的速度场。

进一步地，在上述方法及装置中，所述组合时移量的计算公式为：

式中，τ_i,k表示第k次迭代的组合时移量中的第i个时移量，i的取值范围为1～C_max，所有的τ_i,k共同构成组合时移量，C_max为组合叠加的个数，

为最大时移量，K为迭代总次数。

进一步地，在上述方法及装置中，所述组合时移量中每个时移量对应的加权系数的计算公式为：

λ_i＝2^[0.5*(C_max+1)-abs[0.5*(C_max+1)-i]-1]

式中，λ_i为第k次迭代的组合时移量中第i个时移量对应的加权系数，^为乘方运算，abs[]为取绝对值。

进一步地，在上述方法及装置中，所述第k次迭代的组合编码超级炮通过以下公式得到：

式中，

为第k次迭代的组合编码超级炮，d_n(t-τ_i,k)表示初始观测地震数据第n炮在t时刻经过τ_i,k时移量的炮记录，e_n,k表示极性编码矩阵第k次迭代第n炮的极性值，λ_i为第k次迭代的组合时移量中第i个时移量对应的加权系数，/>

为卷积运算。

进一步地，在上述方法及装置中，所述第k次迭代的组合震源的计算公式为：

式中，

为第k次迭代的组合震源，e_n,k表示极性编码矩阵第k次迭代第n炮的极性值，s_n(t-τ_i,k)表示第n炮主频与所述初始观测地震数据的主频相等的地震子波在t时刻经过τ_i,k时移量的子波，λ_i为第k次迭代的组合时移量中第i个时移量对应的加权系数，/>

为卷积运算。

进一步地，在上述方法及装置中，所述第k次迭代的正演超级炮通过以下公式得到：

式中，D_cal,k(v)为第k次迭代的正演超级炮，L为有限差分正演算子，v为第k-1次迭代的速度场。

进一步地，在上述方法及装置中，所述组合叠加的个数为3。

进一步地，在上述方法及装置中，所述最大时移量的计算公式为：

为最大时移量，f为初始观测地震数据的主频。

附图说明

图1是本发明方法实施例中的地震速度反演方法流程图；

图2是本发明方法实施例中的极性编码矩阵示意图；

图3是本发明方法实施例中使用的速度场模型图；

图4是本发明方法实施例中的初始速度场示意图；

图5是本发明方法实施例中的单震源观测地震数据示意图；

图6是本发明方法实施例中的合成极性编码超级炮示意图；

图7是本发明方法实施例中的组合编码超级炮示意图；

图8是本发明方法实施例中的主频为15hz的雷克子波示意图；

图9是本发明方法实施例中的组合震源示意图；

图10为图8雷克子波的频谱图；

图11为图9组合震源的频谱图；

图12是本发明方法实施例中常规方法得到的反演速度场；

图13是本发明方法实施例中本实施例方法得到的反演速度场；

图14是本发明装置实施例中的地震速度反演装置结构示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。

方法实施例：

如图1所示，本发明的地震速度反演方法包括以下步骤：

步骤1、获取初始观测地震数据、初始速度场和极性编码矩阵；

本实施例中的极性编码矩阵如图2所示，极性编码矩阵是一个K×N_s的矩阵，矩阵中的值为随机生成的正负1，K为迭代总次数，N_s为编码的总炮数；本实施例中的选用如图3所示的国际标准Marmous模型作为速度场模型，初始速度场如图4所示；在实际工程应用中，观测地震数据为单震源观测地震数据，本实施例中的初始观测地震数据为利用单震源打炮N_s次获取的N_s个单震源观测地震数据组成，某个单震源观测地震数据如图5所示。

步骤2、根据初始观测地震数据、初始速度场和极性编码矩阵迭代求取速度场，直至达到迭代总次数或者生成的速度场满足设定要求，得到最终的速度场；

其中，第k次迭代的速度场通过以下步骤得到：

(1)结合最大时移量、迭代总次数和组合叠加的个数得到第k次迭代的组合时移量，1≤k≤K；

其中，最大时移量利用初始观测地震数据的主频计算得到；具体地，利用傅里叶变换对图5所示的初始观测地震数据进行频谱分析，得到初始观测地震数据的主频f，则最大时移量

为：/>

本实施例中主频为15hz，则最大时移量为1/30＝0.033s，取时间间隔为0.001s，则最大时移量为1/30/0.001＝33个计算点。

其中，组合时移量中包含的时移量个数等于组合叠加的个数，组合叠加的个数C_max的取值可根据对地震速度反演结果的精度需求自行设置，本实施例中令C_max＝3就能取得理想的反演效果。

组合时移量的计算公式如下：

式中，τ_i,k表示第k次迭代的组合时移量中的第i个时移量，i的取值范围为1～C_max，所有的τ_i,k共同构成第k次迭代的组合时移量，C_max为组合叠加的个数，

为最大时移量，K为迭代总次数。

(2)利用极性编码矩阵对初始观测地震数据进行编码，得到第k次迭代的合成极性编码超级炮；

第k次迭代的合成极性编码超级炮

通过以下公式得到：

式中，d_n(t)表示初始观测地震数据第n炮在t时刻的炮记录，e_n,k表示极性编码矩阵第k次迭代第n炮的极性值，N_s为编码的总炮数。

第k次迭代的合成极性编码超级炮如图6所示，因为其将很多炮记录合成了一个超级炮，因此计算量大大减少。因为采用极性编码方式，加大了临炮之间的相干性，其可以较好的压制多震源间串扰噪音，更精确的对复杂区域进行速度场的反演，而且由于每次迭代采用的编码极性是随机的，因此其可以改善不同炮记录间的差异性，较好的压制串扰噪音。

(3)利用第k次迭代的组合时移量及组合时移量中每个时移量对应的加权系数和第k次迭代的合成极性编码超级炮，得到第k次迭代的组合编码超级炮；

其中，第k次迭代的组合时移量中第i个时移量对应的加权系数λ_i的计算公式为：

λ_i＝2^[0.5*(C_max+1)-abs[0.5*(C_max+1)-i]-1]

式中，λ_i为第k次迭代的的组合时移量中第i个时移量对应的加权系数，^为乘方运算，abs[]为取绝对值。

第k次迭代的组合编码超级炮

通过以下公式得到：

式中，

为第k次迭代的组合编码超级炮，d_n(t-τ_i,k)表示初始观测地震数据第n炮在t时刻经过τ_i,k时移量的炮记录，e_n,k表示极性编码矩阵第k次迭代第n炮的极性值，/>

为卷积运算。

第k次迭代的组合编码超级炮集如图7所示，因为其先将很多单震源观测地震数据合成为一个极性编码超级炮再进行时移组合，因此计算量大大减少，并且组合编码超级炮的同相轴明显变粗，相比合成极性编码超级炮，地震数据主频明显降低，从而就使得反演过程更加稳定。

(4)利用第k次迭代的组合时移量及组合时移量中每个时移量对应的加权系数、主频与初始观测地震数据的主频相等的雷克子波和极性编码矩阵，得到第k次迭代的组合震源；

本实施例中所用的地震子波为雷克子波，并且本实施例中初始观测地震数据的主频为15hz，因此这里结合第k次迭代的组合时移量、主频为15hz的雷克子波(如图8所示)和极性编码矩阵，得到第k次迭代的组合震源。

第k次迭代的组合震源的计算公式如下：

式中，

为第k次迭代的组合震源，e_n,k表示极性编码矩阵第k次迭代第n炮的极性值，s_n(t-τ_i,k)表示第n炮主频为15hz的雷克子波在t时刻经过τ_i,k时移量的子波。

本实施例中采用主频与初始观测地震数据的主频相等的雷克子波作为正演模拟用的地震子波，不需求取原始子波，一方面能减少子波估计的计算量提高计算效率，另一方面不依赖原始子波，能避免原始子波数据不准确对地震速度反演结果的干扰，提高地震速度反演准确性；作为其他实施方式，还可以利用现有的其他方式获得正演模拟所用的地震子波，例如直接从观测地震数据中提取地震子波。

本实施例中第k次迭代的组合震源如图9所示，发现组合震源的主频宽度明显变宽，通过频谱分析可以发现，组合震源的主频(图11)比未组合前雷克子波的主频(图10)向低频端移动，可以增加反演结果的稳定性。

(5)利用第k次迭代的组合震源和第k-1次迭代的速度场得到第k次迭代的正演超级炮；

具体地，利用第k次迭代的组合震源

对第k-1次迭代的速度场v_k-1利用有限差分进行正演生成第k次迭代的正演超级炮D_cal,k，如果k＝1，则v_k-1为初始速度场。

第k次迭代的正演超级炮通过以下公式得到：

式中，D_cal,k(v)为第k次迭代的正演超级炮，L为有限差分正演算子，v为第k-1次迭代的速度场。

(6)利用第k次迭代的正演超级炮与第k次迭代的组合编码超级炮做差反传，互相关得到第k次迭代的反演梯度场；

反演梯度场的过程即为通过不断迭代最小化一个目标函数的过程，这里的目标函数即为正演超级炮与组合编码超级炮的差，这两者做差的目标函数的具体公式如下：

式中，E(v)表示速度场v对应的数据残差二范数，

为第k次迭代的正演超级炮。

利用第k次迭代的正演超级炮与第k次迭代的组合编码超级炮做差反传，互相关得到反演梯度场，反演梯度场的具体公式如下：

式中，g_k为表示第k次迭代的反演梯度场，δ为求导算子，

为二阶偏导数，L为有限差分正演算子，/>

表示复共轭即共轭转置，<,>表示内积运算。

由于本实施例中采用超级炮反演方法，其计算效率相比普通单震源反演大大提高。

(7)结合第k次迭代的更新步长、第k次迭代的反演梯度场和第k-1次迭代的速度场得到第k次迭代的速度场。

其中，可以利用线性搜索方法或者抛物拟合方法等方法求取更新步长，使得第k次迭代的反演速度场对应的误差最小，本实施例中利用抛物插值方法求取第k次迭代的更新步长α_k，然后结合第k次迭代的反演梯度场和第k-1次迭代的速度场更新速度场，更新公式如下：

v_k＝v_k-1-α_kg_k

式中，v_k表示第k次迭代的反演速度场，g_k为第k次迭代的反演梯度场，v_k-1为第k-1次迭代的反演速度场，如果k＝1，则v_k-1为初始速度场。

其中，第k次迭代的更新步长α_k的计算过程如下：

首先，分别选取两个试探步长α₁、α₂，例如分别选取第k-1次迭代的速度场v_k-1的1％、0.5％作为α₁、α₂(这里的1％、0.5％可以根据实际需要调整为其他数值)；然后，分别将α₁、α₂作为α_k代入速度场更新公式v_k＝v_k-1-α_kg_k，求得α₁、α₂对应的更新后的速度场v_k，分别记为v_k1、v_k2，进而利用E(v)公式求得速度场v_k1、v_k2对应的数据残差二范数E₁、E₂，则第k次迭代的最优步长α_best(即更新步长α_k)可以利用下式获得：

式中，E₀表示第k次迭代的数据残差二范数E(v)。

本实施例中，在每次迭代后均判断此次迭代生成的速度场是否满足设定要求，即判断此次迭代生成的速度场是否能使正演超级炮与组合编码超级炮的差(以下简称速度场误差)小于设定阈值，若速度场误差小于设定阈值，则本次迭代的速度场即为最终的速度场；若达到迭代总次数时的速度场误差仍不小于设定阈值，则从历次迭代得到的速度场中选取一个误差最小的速度场作为最终的速度场。本实施例中的速度场误差是指目标函数E(v)的值。

利用本实施例方法进行反演的过程中，最大时移量逐渐变小，组合编码超级炮的主频逐步提高，实现了多尺度反演的目的，避免对地震数据直接分频滤波的截断效应，在对如图3所示的Marmousi模型进行反演的过程中，使用如图4所示的初始速度场进行反复迭代，应用本实施例方法得到第100次迭代的反演结果如图13所示，可以看出图13中对背景场及复杂构造层位反演更好，可以得到较为准确的反演速度场，消除了反演过程中不稳定的影响。如果不采用组合震源的方法(即采用常规方法)得到的反演结果如图12，可以看出反演的速度场相对于本实施例方法不太准确。

综上所述，本实施例对观测地震数据和地震子波采用极性编码矩阵，充分加大混叠炮之间的差异性，减少临炮之间的干扰，在混叠炮数据反演时较为好的压制串扰噪音，提高了多震源数据速度场反演时的精度。此外，本实施例首先对地震子波和观测地震数据进行时移组合叠加，再进行互相关梯度求取，与传统的多尺度反演方法相比，避免了对观测地震数据进行低通滤波的预处理步骤，其只需要一次逆时偏移的计算量，就可以完成梯度的求取，实现多尺度反演的方法，提高了地震波形反演速度场的精度与效率。在反演过程中，最大时移量逐渐变小，则组合编码超级炮主频逐步提高，实现了多尺度反演的目的，避免对地震数据直接分频滤波的截断效应，提高了地震波形反演速度场的精度。

装置实施例：

本实施例的地震速度反演装置，如图14所示，该装置包括处理器、存储器，存储器中存储有可在处理器上运行的计算机程序，所述处理器在执行所述计算机程序时实现上述方法实施例中的方法。

也就是说，以上方法实施例中的方法应理解为可由计算机程序指令实现地震速度反演方法的流程。可提供这些计算机程序指令到处理器，使得通过处理器执行这些指令产生用于实现上述方法流程所指定的功能。

本实施例所指的处理器是指微处理器MCU或可编程逻辑器件FPGA等的处理装置。

本实施例所指的存储器包括用于存储信息的物理装置，通常是将信息数字化后再以利用电、磁或者光学等方式的媒体加以存储。例如：利用电能方式存储信息的各式存储器，RAM、ROM等；利用磁能方式存储信息的的各式存储器，硬盘、软盘、磁带、磁芯存储器、磁泡存储器、U盘；利用光学方式存储信息的各式存储器，CD或DVD。当然，还有其他方式的存储器，例如量子存储器、石墨烯存储器等等。

通过上述存储器、处理器以及计算机程序构成的装置，在计算机中由处理器执行相应的程序指令来实现，处理器可以搭载各种操作系统，如windows操作系统、linux系统、android、iOS系统等。

Claims (10)

1.一种地震速度反演方法，其特征在于，该方法包括以下步骤：

获取初始观测地震数据、初始速度场和极性编码矩阵；所述极性编码矩阵是一个K×N_s的矩阵，矩阵中的值为随机生成的正负1，K为迭代总次数，N_s为编码的总炮数；

根据所述初始观测地震数据、初始速度场和极性编码矩阵迭代求取速度场，直至达到迭代总次数或者生成的速度场满足设定要求，得到最终的速度场；

其中，第k次迭代的速度场通过以下步骤得到：

结合最大时移量、迭代总次数和组合叠加的个数得到第k次迭代的组合时移量，1≤k≤K；所述最大时移量利用所述初始观测地震数据的主频计算得到，所述组合时移量中包含的时移量个数等于组合叠加的个数；

利用所述极性编码矩阵对所述初始观测地震数据进行编码，得到第k次迭代的合成极性编码超级炮，并利用第k次迭代的组合时移量及组合时移量中每个时移量对应的加权系数和第k次迭代的合成极性编码超级炮，得到第k次迭代的组合编码超级炮；

利用第k次迭代的组合时移量及组合时移量中每个时移量对应的加权系数、主频与所述初始观测地震数据的主频相等的地震子波和极性编码矩阵，得到第k次迭代的组合震源，并利用第k次迭代的组合震源和第k-1次迭代的速度场得到第k次迭代的正演超级炮；其中，k＝1时，第k-1次迭代的速度场为所述初始速度场；

利用第k次迭代的正演超级炮与第k次迭代的组合编码超级炮做差反传，互相关得到第k次迭代的反演梯度场；

结合第k次迭代的更新步长、第k次迭代的反演梯度场和第k-1次迭代的速度场得到第k次迭代的速度场。

2.根据权利要求1所述的地震速度反演方法，其特征在于，所述生成的速度场满足设定要求是指生成的速度场能使正演超级炮与组合编码超级炮的差小于设定阈值；在每次迭代后均判断本次迭代生成的速度场是否满足设定要求，若满足设定要求，则本次迭代的速度场即为最终的速度场；若达到迭代总次数时生成的速度场仍不满足设定要求，则从历次迭代得到的速度场中选取一个能使正演超级炮与组合编码超级炮的差最小的速度场作为最终的速度场。

3.根据权利要求1或2所述的地震速度反演方法，其特征在于，所述组合时移量的计算公式为：

式中，τ_i,k表示第k次迭代的组合时移量中的第i个时移量，i的取值范围为1～C_max，所有的τ_i,k共同构成组合时移量，C_max为组合叠加的个数，

为最大时移量，K为迭代总次数。

4.根据权利要求3所述的地震速度反演方法，其特征在于，所述组合时移量中每个时移量对应的加权系数的计算公式为：

λ_i＝2^[0.5*(C_max+1)-abs[0.5*(C_max+1)-i]-1]

式中，λ_i为第k次迭代的组合时移量中第i个时移量对应的加权系数，^为乘方运算，abs[]为取绝对值。

5.根据权利要求4所述的地震速度反演方法，其特征在于，所述第k次迭代的组合编码超级炮通过以下公式得到：

式中，

为第k次迭代的组合编码超级炮，d_n(t-τ_i,k)表示初始观测地震数据第n炮在t时刻经过τ_i,k时移量的炮记录，e_n,k表示极性编码矩阵第k次迭代第n炮的极性值，λ_i为第k次迭代的组合时移量中第i个时移量对应的加权系数，/>

为卷积运算。

6.根据权利要求4所述的地震速度反演方法，其特征在于，所述第k次迭代的组合震源的计算公式为：

式中，

为第k次迭代的组合震源，e_n,k表示极性编码矩阵第k次迭代第n炮的极性值，s_n(t-τ_i,k)表示第n炮主频与所述初始观测地震数据的主频相等的地震子波在t时刻经过τ_i,k时移量的子波，λ_i为第k次迭代的组合时移量中第i个时移量对应的加权系数，/>

为卷积运算。

7.根据权利要求6所述的地震速度反演方法，其特征在于，所述第k次迭代的正演超级炮通过以下公式得到：

式中，D_cal,k(v)为第k次迭代的正演超级炮，L为有限差分正演算子，v为第k-1次迭代的速度场。

8.根据权利要求3所述的地震速度反演方法，其特征在于，所述组合叠加的个数为3。

9.根据权利要求3所述的地震速度反演方法，其特征在于，所述最大时移量的计算公式为：

为最大时移量，f为初始观测地震数据的主频。

10.一种地震速度反演装置，其特征在于，该装置包括处理器和存储器，所述处理器执行由所述存储器存储的计算机程序，以实现如权利要求1-9任一项所述的地震速度反演方法。

Images