CN114325829B - 一种基于双差思想的全波形反演方法 - Google Patents

Abstract

本发明涉及一种基于双差思想的全波形反演方法，该方法包括以下步骤：1)对原始地震数据进行预处理；2)建立初始速度模型：3)利用双差波形反演方法进行速度迭代更新，获得最终的速度模型。与现有技术相比，本发明能够在光滑模型下减弱地震子波误差对波形反演的影响，获得较高精度与分辨率的速度模型。

Description

一种基于双差思想的全波形反演方法

技术领域

本发明涉及勘探地震学中的速度建模领域，尤其是涉及一种基于双差思想的全波形反演方法。

背景技术

近地表的地下速度模型对于地震数据处理(静校正、偏移成像等)具有关键影响，然而由于地表、地下介质复杂，地下速度建模是一个极具挑战的问题，地球物理学发展了一系列基于反演的地下速度模型构建方法，包含了走时层析与波形反演方法，传统走时层析方法对地震波进行高频近似，利用走时信息，可高效的提供地下速度模型，但是分辨率低，在复杂速度环境下应用效果较差，随着勘探技术的提高，逐渐从利用走时信息走向利用波形信息来刻画地下速度结构，上世纪80年代，Tarantola等率先提出全波形反演理论框架(Full Waveform Inversion，简称FWI)(Tarantola等，1984)，根据全波形反演的定义(PannG S等，1986)，利用合成数据匹配观测数据指导更新速度模型，建立L2范数下观测地震记录与模拟地震记录误差的目标函数，利用优化算法降低目标函数来更新初始模型，因其利用波形信息进行反演，地震运动学与动力学保真，具有较高的精度与分辨率，成为近年来研究的热点，也广泛应用于各种复杂构造的速度建模。

然而全波形反演也存在一定局限性，一方面，它依赖于初始模型的建立(GauthierO等，1986；Mora P，1987)，另一方面如果引入系统误差，如地震子波选取不当，也会影响全波形反演的结果。

Lee，Kim等提出了不依赖子波的全波形反演方法，然而计算量大，需计算Jacobian矩阵(Lee K H等，2003)。Choi等在频率域提出利用褶积法和反褶积法，先对数据与参考道振幅相除再建立目标函数来消除子波影响(Choi Y等，2005)。随后也成功应用于时间域全波形反演(Choi Y等，2011)。敖瑞德等提出基于包络的不依赖子波方法，通过包络对数的目标函数建立全波形反演(敖瑞德等，2015)。除了直接消除子波影响，另一种思路在于尝试估计地震子波进行反演。刘立彬等介绍了一种针对时间域波形反演的子波估计方法，对一给定的试探地震子波，通过频率域对比观测数据来反演预测子波，可较好的吻合真实子波，在全波形反演中取得了不错的效果(刘立彬，2020)。但是与Pratt等的实验相似，需要观测数据良好的信噪比并建立在较准确的初始模型上。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于双差思想的全波形反演方法。

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

一种基于双差思想的全波形反演方法，该方法包括以下步骤：

1)对原始地震数据进行预处理；

2)建立初始速度模型：

3)利用双差波形反演方法进行速度迭代更新，获得最终的速度模型。

所述的步骤1)中，预处理包括去噪、道均衡和去多次波。

所述的步骤2)中，根据先验信息，采用层析成像建立初始速度模型。

所述的步骤3)具体包括以下步骤：

31)基于二维常密度声波方程，进行数值正演模拟计算，获得模拟地震记录d_cal；

32)根据真实地震记录d_obs与模拟地震记录d_cal的残差，计算每一炮的核函数，即伴随震源；

33)基于二维常密度声波方程，利用伴随震源进行波场反传，得到反传波场；

34)利用伴随状态法，计算基于双差波形反演目标函数的梯度；

35)更新速度模型；

36)开始下一次迭代，重复步骤31-35)，直到获得最终的速度反演结果。

所述的步骤31)中，二维常密度声波方程的表达式为：

其中，u(x，z，t)表示声波波场，v(x，z)为速度场，表示质点在(x，z)位置的纵波速度，s(t)表示震源函数，(x_s，z_s)表示震源所处位置，x表示横向位置，z表示深度位置，t表示时间，δ(·)为0/1函数，当自变量取值0时，其值为1，当自变量取值不为0时，则其值为0。

所述的步骤32)中，第k道伴随震源

的计算式为：

其中，Δd_k为第k道模拟地震记录与观测地震记录的残差，

分别为第k道模拟地震记录与观测地震记录，N_r为检波器数目。

所述的步骤34)中，目标函数E的表达式为：

其中，

分别为第k+1道模拟地震记录与观测地震记录，||·||²为L2范数度量的平方。

所述的目标函数E的梯度具体为：

其中，u代表模型下地震的正传波场，λ为以第k道伴随震源

作为震源按声波方程进行反传的反传波场。

所述的步骤35)中，采用LBFGS法更新速度模型。

所述的步骤35)具体包括以下步骤：

351)构造Hessian逆矩阵的正定近似矩阵H_p+1，满足关系式：

H_p+1y_p＝s_p

其中，

为p+1轮迭代速度模型下目标函数的梯度与上一轮迭代速度模型下目标函数的梯度的残差，s_p＝v_p+1-v_p，为速度模型残差。

352)利用有限数目的梯度残差和模型残差信息近似构造矩阵H_p：

其中，

I为单位矩阵，

为对Hessian逆矩阵的初始近似矩阵，r为常数，当迭代次数p大于r时，在计算出该次迭代的模型后，采用最新的一对向量{s_p，y_p}替换掉

中最旧的一对；

353)利用拟Newton法的迭代公式，进行速度模型更新，则有：

其中，v_p为第p轮迭代的速度模型，v_p+1为迭代完成下一次迭代的速度模型，H_p为Hessian逆矩阵的正定近似矩阵，

为第p轮迭代速度模型下目标函数的梯度。

与现有技术相比，本发明具有以下优点：

本发明基于双差思想，能够实现在光滑背景模型下减弱地震子波的误差对波形反演结果的影响，与全波形反演中子波误差项不同在于双差波形反演的子波误差项进行了一次道间作差，使得反演结果对比全波形反演更准确和稳定，获得较高精度与分辨率的速度模型。

附图说明

图1为本发明的方法流程图。

图2为实施例1中的真实速度模型。

图3为实施例1中的初始速度模型。

图4为实施例1中本发明提出的双差波形反演结果。

图5为实施例1中的全波形反演结果。

图6为实施例2中的真实速度模型。

图7为实施例2中的初始速度模型。

图8为实施例2中本发明提出的双差波形反演结果。

图9为实施例2中的全波形反演结果。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

如图1所示，本发明提出一种基于双差思想的全波形反演方法，该方法包括以下步骤：

1)对原始地震数据进行预处理，包括去噪、道均衡和去多次波；

2)建立初始速度模型；

3)利用双差波形反演方法进行速度迭代更新，获得最终的速度模型，具体包括以下步骤：

31)基于二维常密度声波方程，进行数值正演模拟计算，获得模拟地震记录u(x，t)，二维常密度声波方程为：

其中，u(x，z，t)表示声波波场；v(x，z)为速度场，表示质点在(x，z)位置的纵波速度；s(t)表示震源函数，(x_s，z_s)表示震源所处位置；

32)利用真实地震记录d_obs与模拟地震记录d_cal，计算每一炮的核函数，即伴随震源；

伴随震源的计算式，

第k道伴随震源

表示为：

式中，Δd_k为第k道模拟地震记录与观测地震记录的残差，

分别为第k道模拟地震记录与观测地震记录，N_r为检波器数目；

33)基于二维常密度声波方程，利用伴随震源进行波场反传；

34)利用伴随状态法，计算基于双差波形反演目标函数的梯度，目标函数E为：

其中，Δd_k表示第k道模拟地震记录与观测地震记录的残差。

分别为第k道模拟地震记录与观测地震记录，N_r为检波器数目，||·||²为L2范数度量的平方；

目标函数的梯度表达式为：

其中u代表模型下地震的正传波场，λ通过计算得知是以

(该道的模拟地震记录和观测地震记录的残差与下一道的模拟地震记录和观测地震记录的残差的差作为第该道与下一道的伴随震源)作为震源按声波方程进行反传的波场；

35)利用LBFGS法更新速度模型，具体包括以下步骤：

351)拟Newton法构造Hessian逆矩阵的正定近似矩阵H_p+1，满足关系式：

H_p+1y_p＝s_p

其中，

为p+1轮迭代速度模型下目标函数的梯度与上一轮迭代速度模型下目标函数的梯度的残差，s_p＝v_p+1-v_p，为速度模型残差。

352)利用有限数目的梯度残差和模型残差信息近似构造矩阵H_p：

其中，

I为单位矩阵，

为对Hessian逆矩阵的初始近似矩阵，单位矩阵I是一种最简单的初始矩阵选取方式，r＝3～20(本例中选取8)，当迭代次数p大于r时，在计算出该次迭代的模型后，采用最新的一对向量{s_p，y_p}替换掉

中最旧的一对；

353)利用拟Newton法的迭代公式，进行速度模型更新，则有：

其中，v_p为第p轮迭代的速度模型，v_p+1为迭代完成下一次迭代的速度模型，H_p为Hessian逆矩阵的正定近似矩阵，

为第p轮迭代速度模型下目标函数的梯度；

36)开始下一次迭代，重复步骤31-35)，直到获得最终的速度反演结果。

本发明基于双差思想实现在光滑背景模型下减弱地震子波的误差对波形反演结果的影响，得到较高精度与分辨率的速度模型，具体原理为：

定义地震记录为一个算子(格林函数)作用域地震子波：

u＝Gf

其值，u为地震记录，G代表格林函数，f代表震源函数。

假设用于波形反演的格林函数与地震子波为：

G′＝G₀+ΔG，f′＝f₀+Δf

式中，G₀为真实速度模型的格林函数，f₀为真实(观测记录)地震子波，G′为波形反演中初始速度模型的格林函数，f′为用于波形反演的地震子波，ΔG为初始模型与真实速度模型格林函数差，Δf为始模型与真实速度模型地震子波差。

在全波形反演中，目标函数(忽略高阶小项)可写作：

其中，ΔGf₀为由于速度模型差异导致的误差项，在全波形反演中指导更新速度模型，当速度模型等于真实速度模型，ΔG＝0，目标函数为零，但是，当存在子波误差项时，引入了G₀Δf，即当初始模型等于真实模型的情况下，目标函数仍不为零，会错误的更正模型，同时子波误差项的存在也让全波形反演的速度模型更新偏离实际，造成反演结果出错，不稳定。

由此，在本发明的双差波形反演中，第k、k+1地震道的目标函数(忽略高阶小项)可写作：

式中，G^k、G^k+1表示k、k+1道的格林函数，其中(ΔG^k+1-ΔG^k)f₀为由于速度模型差异导致的误差项，指导更新速度模型，

为子波误差项，与现有的全波形反演中子波误差项不同在于，本发明的双差波形反演的子波误差项进行了一次道间作差，一般模型下，由于

，该项不等于零，会遇到同全波形反演一样的问题，但是当真实速度模型较简单平滑，

子波误差项接近于零，在这种情况下，双差波形反演能一定程度上克服子波误差项对反演带来的影响，主要会根据速度误差项(ΔG^k+1-ΔG^k)f₀更新速度模型，使得反演结果对比全波形反演更准确和稳定。

实施例1：

本实施例以二维Marmousi理论模型作为真实模型(如图2所示)，该模型共有298×160个网格，网格间距为20米×20米，速度最大最小值分别为1500米/秒和5500米/秒，地下200米以上为1500米/秒的均匀常速水层，在该模型上进行声波正演模拟，共正演42炮。炮点均匀分布于地表，炮间距为140米，第一炮在水平位置40米处，检波点以20米间隔均匀分布于地表，初始速度模型为光滑后速度模型(如图3所示)，震源函数采用主频10Hz，能量振幅为5的雷克子波，观测时间共4s，时间采样间隔2ms，在应用本发明提出的双差波形反演方法(DDWI)的同时，基于伴随状态法的全波形反演方法(FWI)也同时被应用。

利用DDWI与FWI进行波形反演的结果分别如图4和5所示，双差波形反演与全波形反演均能有效进行速度反演，甚至前者深部的反演分辨率更高，可见Marmousi模型的典型特征和细节特征，尤其是中浅部，可见，对于复杂的Marmousi模型，双差波形反演与全波形反演结果相仿，都能有效准确地进行速度反演。

实施例2

本实施例以简单高速异常球体模型作为真实速度模型(如图6所示)，该模型共有101×101个网格，网格间距为5米×5米，背景速度为3000米/秒，在中心(250米，250米)位置设置一半径50米的异常球体，速度为3300米/秒，地下25米以上、地下475米以下为3000米/秒的均匀常速层，在该模型上进行声波正演模拟，共正演25炮，炮点均匀分布于地表，炮间距为20米，第一炮在水平位置10米处，检波点以5米间隔均匀分布于底部，初始速度模型为3000米/秒背景速度，无高速异常球体(如图7所示)，震源函数采用主频30Hz，能量振幅为5的雷克子波，在模拟地震记录中采用主频50Hz，能量振幅为1的雷克子波，观测时间共0.3s，时间采样间隔0.5ms，在应用本发明提出的双差波形反演方法(DDWI)的同时，基于伴随状态法的全波形反演方法(FWI)也同时被应用。

利用DDWI与FWI进行波形反演的结果分别如图8和9所示，全波形反演结果效果很差，反演结果不稳定，且不可见中心异常球体特征，偏离实际，相比之下，虽然双差波形反演结果也些许变差，中心异常球体轮廓尚存但是形状变大。但是总体上依旧符合真实速度模型特征。

因而双差波形反演在光滑背景模型下确实比全波形反演对子波误差不敏感，反演表现比传统全波形反演要好，即双差波形反演在子波误差存在的情况下，能一定程度上克服子波误差的影响。

Claims (6)

1.一种基于双差思想的全波形反演方法，其特征在于，该方法包括以下步骤：

1)对原始地震数据进行预处理；

2)建立初始速度模型：

3)利用双差波形反演方法进行速度迭代更新，获得最终的速度模型；

所述的步骤3)具体包括以下步骤：

31)基于二维常密度声波方程，进行数值正演模拟计算，获得模拟地震记录d_cal；

32)根据真实地震记录d_obs与模拟地震记录d_cal的残差，计算每一炮的核函数，即伴随震源；

33)基于二维常密度声波方程，利用伴随震源进行波场反传，得到反传波场；

34)利用伴随状态法，计算基于双差波形反演目标函数的梯度；

35)更新速度模型；

36)开始下一次迭代，重复步骤31-35)，直到获得最终的速度反演结果；

所述的步骤31)中，二维常密度声波方程的表达式为：

其中，u(x,z,t)表示声波波场，v(x,z)为速度场，表示质点在(x,z)位置的纵波速度，s(t)表示震源函数，(x_s,z_s)表示震源所处位置，x表示横向位置，z表示深度位置，t表示时间，δ(·)为0/1函数，当自变量取值0时，其值为1，当自变量取值不为0时，则其值为0；

所述的步骤32)中，第k道伴随震源

的计算式为：

其中，Δd_k为第k道模拟地震记录与观测地震记录的残差，

分别为第k道模拟地震记录与观测地震记录，N_r为检波器数目；

所述的步骤34)中，目标函数E的表达式为：

其中，

分别为第k+1道模拟地震记录与观测地震记录,||.||²为L2范数度量的平方。

2.根据权利要求1所述的一种基于双差思想的全波形反演方法，其特征在于，所述的步骤1)中，预处理包括去噪和道均衡。

3.根据权利要求1所述的一种基于双差思想的全波形反演方法，其特征在于，所述的步骤2)中，根据先验信息，采用层析成像建立初始速度模型。

4.根据权利要求1所述的一种基于双差思想的全波形反演方法，其特征在于，所述的目标函数E的梯度具体为：

其中，φ代表模型下地震的正传波场，λ为以第k道伴随震源

作为震源按声波方程进行反传的反传波场。

5.根据权利要求1所述的一种基于双差思想的全波形反演方法，其特征在于，所述的步骤35)中，采用LBFGS法更新速度模型。

6.根据权利要求5所述的一种基于双差思想的全波形反演方法，其特征在于，所述的步骤35)具体包括以下步骤：

351)构造Hessian逆矩阵的正定近似矩阵H_p+1，满足关系式：

H_p+1y_p＝s_p

其中，

为p+1轮迭代速度模型下目标函数的梯度与上一轮迭代速度模型下目标函数的梯度的残差，s_p＝v_p+1-v_p，为速度模型残差；

352)利用有限数目的梯度残差和模型残差信息近似构造矩阵H_p：

其中，

I为单位矩阵，

为对Hessian逆矩阵的初始近似矩阵，r为常数，当迭代次数p大于r时，在计算出该次迭代的模型后，采用最新的一对向量{s_p,y_p}替换掉

中最旧的一对；

353)利用拟Newton法的迭代公式，进行速度模型更新，则有：

其中，v_p为第p轮迭代的速度模型，v_p+1为迭代完成下一次迭代的速度模型，H_p为Hessian逆矩阵的正定近似矩阵，

为第p轮迭代速度模型下目标函数的梯度。

Images