一种拖缆地震数据和少量OBS数据联合的波形反演方法

摘要

本发明涉及一种拖缆地震数据和少量OBS数据联合的波形反演方法，属于海洋地震勘探的全波形反演速度建模领域。本发明主要包括如下步骤：1)获取观测的拖缆和OBS地震记录数据，并对其进行预处理；2)采用高阶统计方法提取地震子波；3)分别计算拖缆和OBS地震记录对应的正传波场和反传波场；4)基于伴随状态法计算目标函数的梯度场；5)计算共轭梯度场和速度模型的修正量；6)迭代更新速度模型；7)判断是否满足迭代终止条件，如果不满足迭代终止条件，则返回步骤3，否则，输出结果。本方法解决了只利用拖缆地震记录进行波形反演，由于缺少低频地震信息，反演结果不准确的问题。

说明书

技术领域

本发明涉及一种拖缆地震数据和少量OBS数据联合的波形反演方法，属于地震勘探的全波形反演速度建模领域。

背景技术

目前速度建模的方法有很多种，比如走时层析、斜率层析、菲尼尔体层析、偏移速度分析、全波形反演等。在这些方法中，全波形反演是近年来的研究热点之一，也是最有可能得到高分辨率、高保真度的速度建模方法之一。

然而，全波形反演严重依赖于初始速度模型的准确程度和观测地震记录中的低频地震信息。对于地质结构复杂的区域，利用现有的手段为全波形反演建立一个足够准确的初始速度模型通常比较困难。因此，为全波形反演提供足够低频的地震数据就显得尤为重要。但是，由于常用的采集技术的限制，由水听器记录的拖缆地震数据中通常缺少低频成份(< 5 Hz)。利用拖缆地震数据进行全波形反演，其反演通常收敛于局部解，反演结果不可信。

为此，许多学者提出了多种方法来构建低频地震信息。一种方法是利用高频地震资料通过非线性插值计算低频地震资料。该方法虽然在一定程度上可以抑制周期跳跃现象的发生，然而，插值得到的低频地震资料不满足波动方程。Beat tone 全波形反演利用数学变换从频率相近的两个高频信号中得到低频地震资料，然而构建的低频地震资料的质量严重依赖于高频地震数据。拉普拉斯域全波形反演在缺少低频地震资料时，可以得到足够平滑的速度模型，然而该方法需要大偏移距的地震数据。包络反演也可以构建低频信息，然而数学变换严重地破坏了地震波的波形。

这些方法主要通过数学变换来得到低频的地震信息，从而满足全波形反演的需要。然而，这些人工合成的低频地震记录与真实的低频地震信息存在一定程度的差异，使得反演结果中可能存在虚假的速度结构，降低了反演结果的可靠性。

发明内容

针对现有技术存在的上述缺陷，本发明提出了一种拖缆地震数据和少量OBS数据联合的波形反演方法，利用少量OBS观测获取实际的地震低频分量，解决了只使用拖缆数据进行全波形反演时，因缺乏低频分量使反演结果不可靠的问题。

本发明是采用以下的技术方案实现的：本发明所述一种拖缆地震数据和少量OBS数据联合的波形反演方法，包括如下步骤：

1）构建拖缆地震数据和少量OBS数据联合的混合域全波形反演目标函数，构建的目标函数如下：

，（1）

式中，E是目标函数值，α是权重系数，n_r是每一炮水听器的个数，s是炮号，r是水听器号，t是时间，x_s是模型空间中炮点s的位置，n_o是每一炮OBS的个数，o指OBS号，d_cal(r,t|x_s)是预测的拖缆地震记录，d_rea(r,t|x_s)是观测的拖缆地震记录，d_cal(o,t|x_s)是预测的OBS地震记录，d_rea(o,t|x_s)是观测的OBS地震记录；本发明通过最小化该目标函数来得到最优的反演结果；

2）获取野外观测的拖缆和OBS地震记录；

3）对观测数据进行预处理；通过带通滤波、F-K滤波等数据处理手段从观测数据中去除声波方程无法正演模拟的面波、转换波、随机噪声等；

4）采用高阶统计方法，从观测资料中提取地震子波；

5）计算权重系数；在反演迭代的初期，主要恢复模型的背景速度场，所以OBS资料应该起主导作用，而在反演后期，主要提高模型的分辨率，所以拖缆数据应该占主要地位；因此，本发明采用下式计算权重系数：

，（2）

式中，k是当前的迭代次数，n_ite是反演总共的迭代次数；

6）计算预测的拖缆地震记录及其对应的正传波场；利用给定的初始速度模型，计算正传波场，正传波场的计算公式如下：

，（3）

式中，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，v(x)是模型空间x点处的速度，P_fr(x,t|x_s)是拖缆地震记录对应的时间域正传波场，s_fr(x,t|x_s)是从预处理后观测的拖缆地震记录中提取的子波；预测的拖缆地震记录的计算公式如下：

，（4）

式中，r表示水听器号，t表示时间，x_s表示模型空间中炮点位置，x表示模型空间位置，x_r是水听器在模型空间的位置，P_fr(x,t|x_s)表示时间域正传波场，R表示将检波点处的波场值提取出来作为该时刻检波点的地震记录，d_cal(r,t|x_s)表示预测的拖缆地震记录；

7）计算预测的OBS地震记录及其对应的正传波场；利用给定的初始速度模型，计算正传波场，正传波场的计算公式如下：

，（5）

式中，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，v(x)是模型空间x点处的速度，P_fo(x,t|x_s)是OBS地震记录对应的时间域正传波场，s_fo(x,t|x_s)是从预处理后观测的OBS地震记录中提取的子波；预测的OBS地震记录的计算公式如下：

，（6）

式中，o表示OBS号，t表示时间，x_s表示模型空间中炮点位置，x表示模型空间位置，x_o是OBS在模型空间的位置，P_fo(x,t|x_s)表示时间域正传波场，R表示将OBS点处的波场值提取出来作为该时刻OBS的地震记录，d_cal(o,t|x_s)表示预测的OBS地震记录；

8）将预测地震记录和观测地震记录带入公式1，计算目标函数值；

9）计算拖缆地震记录的残差和拖缆地震记录对应的反传波场；拖缆地震记录残差的计算公式如下：

，（7）

式中，s是炮号，r是水听器号，t是时间，x_s表示模型空间中炮点位置；拖缆地震记录对应的反传波场的计算公式如下：

，（8）

式中，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，v(x)是模型空间x点处的速度，P_br(x,t|x_s)是拖缆地震记录对应的时间域反传波场，s_br(x,t|x_r)是第s炮拖缆地震记录的残差；

10）计算OBS地震记录的残差和OBS地震记录对应的反传波场；OBS地震记录残差的计算公式如下：

，（9）

式中，s是炮号，o是OBS号，t是时间，x_s表示模型空间中炮点位置；OBS地震记录对应的反传波场的计算公式如下：

，（10）

式中，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，v(x)是模型空间x点处的速度，P_bo(x,t|x_s)是OBS地震记录对应的时间域反传波场，s_bo(x,t|x_r)是第s炮OBS地震记录的残差；

11）将时间域拖缆和OBS地震记录对应的正传和反传波场变换到频率域，计算公式如下：

，（11）

式中，T是采样时间长度，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，ω表示角频率；

12）计算目标函数的梯度；梯度的计算公式如下：

，（12）

式中，α 是权重系数，n_r是水听器的个数，n_o是OBS的个数，ω_r是使用拖缆数据时选择的角频率，ω_o是使用OBS数据时选择的角频率，s是炮号，x是模型空间的位置，x_s是炮点在模型空间的位置，P_fr(x,ω_r|x_s)是拖缆数据对应的正传波场，P_fo(x,ω_o|x_s)是OBS数据对应的正传波场，P_br(x,ω_r|x_s)是拖缆数据对应的反传波场，P_bo(x,ω_o|xs)是OBS数据对应的反传波场，上标*表示复共轭转置；

13）计算共轭梯度场；共轭梯度场的计算公式如下：

，（13）

式中，d_k-1是第k-1次的共轭梯度，β_k是使相邻两次共轭梯度正交的系数，其计算公式如下：

，(14)

其中，g_k是梯度向量，d_k-1是共轭梯度向量，符号<,>表示内积运算；

14）迭代更新速度模型；速度模型更新的计算公式如下：

，（15）

式中，v_k+1是第k次迭代更新后的速度，λ是迭代更新的步长，可以通过线性搜索的方法得到；

15）判断是否满足迭代终止条件；如果满足终止条件，则输出计算结果，否则，将更新后的速度模型作为新的初始速度模型，返回步骤5）。

进一步地，所述步骤4）中，提取的地震子波有两个，分别是从观测的拖缆地震记录提取的子波s_fr和从观测的OBS地震记录中提取的地震子波s_fo。

进一步地，公式3、5、8和10都采用时间2阶、空间12阶规则网格有限差分方法求解，同时采用完全匹配层吸收边界条件来压制模型边界处的反射波。

本发明的有益效果是：海洋拖缆地震数据增加少量OBS地震数据的联合地震波形反演，在反演初期，主要利用少量OBS地震数据来反演模型的背景速度场，而在反演后期，主要利用拖缆地震数据来反演模型的高频成份；克服了只利用拖缆数据进行反演，由于缺少低频地震信息，导致反演结果不准的问题；其次，相对于现有方法所利用的人工计算的低频地震记录， OBS观测的地震记录的低频信息，更真实、更可靠，采用OBS地震记录进行反演，反演结果的可靠性高。

附图说明

图1为本发明的流程图。

图2为本发明建立的真实速度模型图。

图3为本发明所用的初始速度模型图。

图4为第17炮拖缆地震记录图。

图5为第17炮OBS地震记录图。

图6为常规全波形反演方法只使用拖缆数据的反演结果图。

图7为常规全波形反演方法只使用OBS地震记录的反演结果图。

图8为本发明的反演结果图。

具体实施方式

为了使本发明的目的、技术方案更加清楚明白，下面以理论模型为例并结合附图，对本发明作进一步详细说明。

本发明的流程图，如图1所示，主要包括如下步骤：

1）获取观测地震记录；给定一个真实速度模型，采用高频的雷克子波作为震源，利用常密度声波方程正演模拟，将得到的地震记录作为观测的拖缆地震记录；采用低频的雷克子波作为震源，利用常密度声波方程正演模拟，将得到的地震记录作为观测的OBS地震记录；正演模拟的计算公式如下：

，（1）

式中，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，v(x)是模型空间x点处的速度，P(x,t|x_s)表示时间域波场，s(x,t|x_s)表示震源函数；

2）计算权重系数；在反演迭代的初期，主要恢复模型的背景速度场，所以OBS资料应该起主导作用，而在反演后期，主要提高模型的分辨率，所以拖缆数据应该占主要地位；本发明计算权重系数的公式如下：

，（2）

式中，k是当前的迭代次数，n_ite是反演总共的迭代次数；

3）计算正传波场；给定初始速度模型和高频雷克子波s_fr，利用声波方程正演模拟得到拖缆地震记录对应的时间域正传波场P_fr；正传波场P_fr的计算公式如下：

，（3）

式中，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，v(x)是模型空间x点处的速度，s_fr(x,t|x_s)是高频雷克子波；

利用给定的初始速度模型和低频雷克子波s_fo，利用声波方程正演模拟得到OBS地震记录对应的正传波场P_fo；正传波场P_fo的计算公式如下：

，（4）

式中，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，v(x)是模型空间x点处的速度，s_fo(x,t|x_s)是低频的雷克子波；

4）计算预测的地震记录；预测的拖缆地震记录的计算公式如下：

，（5）

式中，r表示水听器号，t表示时间，x_s表示模型空间中炮点位置，x表示模型空间位置，x_r是水听器在模型空间的位置，P_fr(x,t|x_s)表示拖缆地震记录对应的时间域正传波场，d_cal(r,t|x_s)表示预测的拖缆地震记录，R(x|x_r,x_s)表示数据选择函数，当x等于x_r时，R等于1，否则为0；

预测的OBS地震记录的计算公式如下：

，（6）

式中，o表示OBS号，t表示时间，x_s表示模型空间中炮点位置，x表示模型空间位置，x_o是OBS在模型空间的位置，P_fo(x,t|x_s)表示OBS地震记录对应的时间域正传波场，d_cal(o,t|x_s)表示预测的OBS地震记录，R(x|x_o,x_s)表示数据选择函数当x等于x_o时，R等于1，否则为0；

5）计算地震记录残差；将预测的拖缆地震记录和观测的拖缆地震记录作差，得到拖缆地震记录的残差，计算公式如下：

，（7）

式中，s是炮号，r是水听器号，t是时间，x_s表示模型空间中炮点位置，d_cal(r,t|x_s)是预测的拖缆地震记录，d_rea(r,t|x_s)是观测的拖缆地震记录，d_err(r,t|x_s)是拖缆地震记录的残差；

将预测的OBS地震记录和观测的OBS地震记录作差，得到OBS地震记录的残差，计算公式如下：

，（8）

式中，s是炮号，o是OBS号，t是时间，x_s表示模型空间中炮点位置，d_cal(o,t|x_s)是预测的OBS地震记录，d_rea(o,t|x_s)是观测的OBS地震记录，d_err(o,t|x_s)是OBS地震记录的残差；

6）计算目标函数值，目标函数值的计算公式如下：

，（9）

式中，E是目标函数值，α是权重系数，n_r是每一炮水听器的个数，s是炮号，r是水听器号，t是时间，x_s是模型空间中炮点s的位置，n_o是每一炮OBS的个数，o指OBS号，d_cal(r,t|x_s)是预测的拖缆地震记录，d_rea(r,t|x_s)是观测的拖缆地震记录，d_cal(o,t|x_s)是预测的OBS地震记录，d_rea(o,t|x_s)是观测的OBS地震记录；

7）计算反传波场；将地震记录残差作为震源函数，利用给定的初始速度模型，计算反传波场；拖缆地震记录对应的反传波场P_br的计算公式如下：

，（10）

式中，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，v(x)是模型空间x点处的速度，s_br(x,t|x_r)是第s炮拖缆地震记录的残差；

OBS地震记录对应的反传波场P_bo的计算公式如下：

，（11）

式中，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，v(x)是模型空间x点处的速度，s_bo(x,t|x_r)是第s炮OBS地震记录的残差；

8）将时间域拖缆和OBS数据对应的正传、反传波场变换到频率域，其计算公式如下：

，（12）

式中，T是采样时间长度，x表示模型空间位置，t表示时间，x_s表示模型空间中炮点位置，ω表示角频率，P(x,t|x_s)表示时间域波场，P(x,ω|x_s)表示频率域波场；

9）计算目标函数的梯度；利用步骤3）计算出的正传波场和步骤7）计算出的反传波场，来计算目标函数的梯度场，计算公式如下：

，（13）

式中，α是权重系数，n_r是水听器的个数，n_o是OBS的个数，ω_r是使用拖缆数据时选择的角频率，ω_o是使用OBS数据时选择的角频率，s是炮号，x是模型空间的位置，x_s是炮点在模型空间的位置，P_fr(x,ω_r|x_s)是拖缆数据对应的正传波场，P_fo(x,ω_o|x_s)是OBS数据对应的正传波场，P_br(x,ω_r|x_s)是拖缆数据对应的反传波场，P_bo(x,ω_o|x_s)是OBS数据对应的反传波场，上标*表示复共轭转置；

10）计算目标函数的共轭梯度；共轭梯度的计算公式如下：

，（14）

式中，d_k-1是第k-1次的共轭梯度，β_k是一个使相邻两次共轭梯度正交的系数，其计算公式如下：

，(15)

式中，g_k是梯度向量，d_k-1是共轭梯度向量，符号<,>表示内积运算；

11）计算步长；步长通过线性搜索方法计算，计算公式如下：

，（16）

式中，ε是试探步长，向量d^’_cal是当λ=ε时，更新速度模型后，合成的拖缆地震记录，向量d_cal是合成的拖缆地震记录，向量d_rea是观测的拖缆地震记录；试探步长的计算公式如下：

，（17）

式中，向量v表示模型空间中所有点的速度，d_k表示共轭梯度向量，max表示从向量中选择最大值，符号||表示绝对值运算，即当向量中元素值为负数时，取其相反数；

12）更新速度模型；速度模型更新的计算公式如下：

，（18）

式中，v_k+1是第k次迭代更新后的速度，λ_k是第k次迭代更新的步长；

13）判断是否满足迭代终止条件；如果满足终止条件，则输出计算结果，否则，将更新后的速度模型作为新的初始速度模型，返回步骤2）。

其中，公式1、3、4、10和11中的波场通过时间2阶、空间12阶精度的规则网格有限差分方法求解，采用完全匹配层边界条件来压制模型边界处的反射波。

其中，所述步骤3）和7）中，反演所用的初始速度模型，通过对真实速度模型平滑得到。

为了更好的说明本发明，下面列举一实施例。

为了进一步说明本方法的实现思路及实现过程并证明方法的有效性，用marmousi2模型进行测试，并和常规全波形反演的结果进行比较。

S1：建立一个真实速度模型，如图2所示。真实速度模型宽度为9216 m，深度为3048m。采用正方形网格离散，网格大小为24 × 24 m。

S2：观测系统：采用海上拖缆数据采集观测系统，从水平方向3360 m处开始以96 m的间隔布设了51个炮点。炮点从左向右依次激发。拖缆上共有140个水听器接收地震记录，道间距24 m，最小偏移距为24 m，最大偏移距为3360 m。11个海底观测站以864 m的间隔均匀布设在水平方向120 m 至 8760 m的海底。地震记录采样时间为3 s，时间采样间隔为2ms。

S3：由真实速度模型（详见图2）和震源函数为10 Hz的雷克子波，利用公式1，通过时间2阶、空间12阶精度的规则网有限差分方法，采用完全匹配层边界条件，正演模拟得到观测的拖缆地震记录。图4是第17炮观测的拖缆地震记录。

S4：由真实速度模型和震源函数为5 Hz的雷克子波，利用公式1，采用同样的正演模拟方法，得到观测的OBS地震记录。图5是第17炮观测的OBS地震记录。

S5：利用公式2，计算权重系数。

S6：由初始速度模型（详见图3）和震源函数为10 Hz的雷克子波，利用公式3和5，采用同样的正演模拟方法，得到预测的拖缆地震记录和拖缆地震记录对应的正传波场。

S7：由初始速度模型和震源函数为5 Hz的雷克子波，利用公式4和6，采用同样的正演模拟方法，得到预测的OBS地震记录和OBS地震记录对应的正传波场。

S8：由观测的拖缆地震记录和预测的拖缆地震记录，利用公式7，计算拖缆地震记录的残差。

S9：由初始速度模型和拖缆地震记录的残差，利用公式10，计算拖缆地震记录对应的反传波场。

S10：由观测的OBS地震记录和预测的OBS地震记录，利用公式8，计算OBS地震记录的残差。

S11：由步骤8和10计算的地震记录残差，利用公式9，计算目标函数值。

S12：由初始速度模型和OBS地震记录的残差，利用公式11，计算OBS地震记录对应的反传波场。

S13：由S5、S6、S8和S10计算的时间域正传和反传波场，利用公式12，计算频率域正传和反传波场，并根据计算的频率波场，利用公式13，计算目标函数的梯度。在本实施例中，共使用了12个频率，分别是拖缆地震记录的6个频率，即6、7.32、8.94、10.91、13.32、16.26Hz，OBS地震记录的6个频率，即2、2.40、2.88、3.46、4.15、4.98 Hz。将12个频率分为6组（详见表1），对于每一个频率组，都迭代15次。

S14：利用公式14计算共轭梯度场。

S15：利用公式16和17计算步长。

S16：利用公式18更新速度模型。判断更新后的速度模型是否满足迭代终止条件，如果满足，则输出结果；如果不满足，则返回步骤5。

图6是只使用拖缆数据，常规全波形反演的结果。和真实速度模型对比，可以发现反演结果只有浅层接近真实速度模型，而600 m以下的区域并没有正确的反演出来。图7是只利用OBS数据的反演结果。和真实速度模型对比，可以发展，反演结果的分辨率较低，反演结果不能满足需要。图8是本发明的反演结果，和图6对比，可以发现，模型600 m以下区域的速度得到了正确的反演，和图7对比，可以发现模型的分辨率得到了显著的提高。

以上所述仅为本发明的较佳实施例而己，并不以本发明为限制，凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等，均应包含在本发明的专利涵盖范围内。

表1 本实施例使用的频率组