CN107505651A - 地震初至波和反射波联合斜率层析成像方法 - Google Patents

Abstract

本发明涉及一种地震初至波和反射波联合斜率层析成像方法，属于油气勘探地震资料处理技术领域。本发明解决了现有的起伏地形速度模型的构建方法存在的在不同地震资料和方法建立的浅部和中深部速度模型存在差异，甚至会出现变化趋势不同的情况。本发明同时使用反射波斜率与初至波斜率，形成地震初至波和反射波联合斜率层析成像方法，同时建立由浅到深部的统一完整速度模型，适用于地形起伏的情况，即使是在近地表速度横向变化的情况，也能得到精确的宏观速度模型；联合斜率层析成像同时使用初至波和反射波信息，扩大了射线的角度覆盖并增加射线密度，获得比偏移速度分析更好的反演结果，对复杂情况的适用性更强。

Description

地震初至波和反射波联合斜率层析成像方法

技术领域

本发明涉及一种地震初至波和反射波联合斜率层析成像方法，属于油气勘探地震资料处理技术领域。

背景技术

斜率层析成像是1998年提出的应用反射波的走时和斜率信息获取宏观速度模型的层析成像方法，并在2000年应用于海上地震反射波数据，取得了良好的效果。同年，三维并行斜率层析成像得以实现，并于2002年应用到海上三维地震数据中。之后斜率层析成像的数据拾取工作在多种数据域中进行了实验，时间域或深度偏移域，叠后时间域或叠前时间偏移域等。斜率层析成像扩展到转换波分析、透射波层析、各向异性介质中波的传播等方面，形成了包括PP/PS各向异性斜率层析成像在内的多种方法。除以上研究成果外，斜率层析成像在约束方面也有一定进展，先后有学者提出了基于反射角的结构性平滑约束和倾斜约束。此外，在模型试验中，斜率层析成像还被用于为全波形反演提供初始速度模型。最近几年，斜率层析成像又逐步衍生出了基于三角网格的斜率层析成像和伴随状态的斜率层析成像。从斜率层析成像方法十多年的发展进程中可以看出，斜率层析成像方法研究主要在于海上反射波或透射波方面，面向陆上地震勘探的方法研究较少，目前还没有针对复杂地形的斜率层析成像方法。

另一方面，随着我国陆上油气勘探重点从东部转移到西部，近地表问题成为地震勘探中重点和难点。西部地区地形起伏，近地表结构复杂，速度变化大，容易引起反射同相轴扭曲、畸变，严重影响地震成像。此外，对于复杂的地表情况，传统的地表一致性静校正不能完全适用。于是，出现了基于起伏地形的地震成像方法相关研究。在该类算法中，深度成像需要建立由浅到深的速度模型。一种常用的做法是用反射波建立深部模型，用初至波建立浅层模型，再把深部和浅部模型融合在一起。但是应用不同地震资料和方法建立的浅部和深部速度模型难免存在差异，甚至会出现变化趋势不同的情况。

发明内容

针对现有技术存在的上述缺陷，本发明提出了一种地震初至波和反射波联合斜率层析成像方法,同时使用反射波斜率与初至波斜率，共同建立由浅到深部的完整统一的速度模型，以用于陆上地震勘探资料的叠前深度偏移成像。

本发明是采用以下的技术方案实现的：本发明所述的地震初至波和反射波联合斜率层析成像方法，包括如下步骤：

步骤一：地震数据的预处理：预处理工作包括去噪和增益，去噪目的是压制随机噪声、提高信噪比、去除相干噪声，而增益目的是平衡时间轴方向上能量分布；

步骤二：斜率和走时数据拾取与质量控制：

S1:斜率拾取是在共炮点道集和共检波点道集中实现的，采用对地震道包络进行倾斜叠加的方式，获得最大叠加能量对应的斜率；

S2:此后，通过人机交互方式进行质量控制，最终获得：

初至波数据：

和反射波数据：

其中，

分别为激发点S的横纵坐标以及初至波斜率；

分别为接收点R的横纵坐标以及初至波斜率；

双程初至走时为

N₁为初至波数据个数；

分别为激发点S的横纵坐标以及反射波斜率；

分别为接收点R的横纵坐标以及反射波斜率；

双程反射走时为

N₂为反射波数据个数；

步骤三：初始化速度模型：

S3:速度模型参数m_v的初始值采用随深度线性增大的函数v＝(v₀+k·z)km/s；其中，

v₀为地表速度；

k为速度垂直梯度；

z表示深度；

或采用已获得的先验速度作为初始速度；

S4:与此同时，基于激发点和接收点位置构建地形；

步骤四：射线段参数的初始化及其优化反演：

S5:设定初至波射线段参数表示为：

反射波射线段参数表示为：

其中，

和分别为初至波回折点的位置和反射点的位置；

和分别为初至波回折点c到激发点S、接收点R的单程走时和反射点C到激发点S、接收点R的单程走时，为反射射线在点C朝着激发点S和接收点R的出射角；

S6:取单程走时等于双程走时的一半，回折点c和反射点C的横坐标近似为炮点与检波点的中心位置，纵坐标选择是通过的给定一系列深度值进行正演，选取正演结果最符合实际数据的深度值作为初至值，射线从反射点出发的初始角度由反射点与炮点、检波点的位置决定；

S7:初始化后，基于初始速度，采用阻尼SVD法进一步优化射线段参数；

步骤五：速度模型和射线段参数的联合反演：

S8:对模型参数进行归一化，得到

将联合层析成像的目标函数设为下式：

其中，||.||为二范数运算，φ为权重，取值范围在0到1之间，C_d为数据协方差矩阵，d_{fst_cal}表示正演初至数据，d_{ref_cal}表示正演反射数据，λ为阻尼系数，L为速度扰动的拉普拉斯运算符，n_{v_prior}为先验的速度模型，添加该项用于抑制反演过程中的多解性问题；

S9:采用阻尼最小二乘方法原理，将以上非线性问题线性化，以迭代的形式进行求解；并设n_{v_prior}为上次反演得到的结果，可得第k次迭代的求解的方程组如下：

其中，

表示第k次迭代初至波数据空间对模型空间的敏感度矩阵；

表示第k次迭代反射波数据空间对模型空间的敏感度矩阵；

为第k次迭代待求的未知量；

为第k次迭代的正演初至波数据和拾取初至波数据之差；

为第k次迭代的正演反射波数据和拾取反射波数据之差；

S10:以上公式中的敏感度矩阵可由傍轴射线近似公式计算得到；

S11:构成得线性方程组采用LSQR算法求解；

S12:求得以上方程组后可进一步计算出模型更新量；

步骤六：检查反演结果：反演之后要及时检查反演结果是否收敛、合理，若反演结果合适，则反演结束，否则反演继续。

其中，步骤六的反演过程采用了多尺度策略，以达到较好的、需要的模型分辨率和反演精度；此外，反演后对反演所用数据进行筛选，对一些误差较大的数据进行剔除或重新赋值处理。

本发明的有益效果是：采用本发明所述的地震初至波和反射波联合斜率层析成像方法，(1)针对陆上地震勘探设计，有效地建立由浅到深的速度模型，适用于地形起伏的情况，即使是在近地表速度横向变化的情况，也能得到精确的宏观速度模型；(2)联合斜率层析成像同时使用初至波和反射波信息，扩大了射线的角度覆盖并增加射线密度，获得比基于反射波的偏移速度分析更好的反演结果，对复杂情况的适用性更强。

附图说明

图1是本发明的结构示意图。

图2(a)预处理前单炮记录AGC显示图。

图2(b)预处理后单炮记录AGC显示图。

图3初始速度模型图。

图4(a)联合斜率层析的反演速度图。

图4(b)联合斜率层析的深度偏移剖面图。

图5(a)偏移速度分析局部结果的反演速度图。

图5(b)偏移速度分析局部结果的深度偏移剖面。

图5(c)偏移速度分析局部结果的CIG道集。

图6(a)联合斜率层析局部结果的反演速度图。

图6(b)联合斜率层析局部结果的深度偏移剖面。

图6(c)联合斜率层析局部结果的CIG道集。

具体实施方式

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

如图1所示，本发明所述的地震初至波和反射波联合斜率层析成像方法，包括如下步骤：

步骤一：地震数据的预处理：预处理工作包括去噪和增益，去噪目的是压制随机噪声、提高信噪比、去除相干噪声，而增益目的是平衡时间轴方向上能量分布；

步骤二：斜率和走时数据拾取与质量控制：

S1:斜率拾取是在共炮点道集和共检波点道集中实现的，采用对地震道包络进行倾斜叠加的方式，获得最大叠加能量对应的斜率；

S2:此后，通过人机交互方式进行质量控制，最终获得：

初至波数据：

和反射波数据：

其中，

分别为激发点S的横纵坐标以及初至波斜率；

分别为接收点R的横纵坐标以及初至波斜率；

双程初至走时为

N₁为初至波数据个数；

分别为激发点S的横纵坐标以及反射波斜率；

分别为接收点R的横纵坐标以及反射波斜率；

双程反射走时为

N₂为反射波数据个数；

步骤三：初始化速度模型：

S3:速度模型参数m_v的初始值采用随深度线性增大的函数v＝(v₀+k·z)km/s；

其中，

v₀为地表速度；

k为速度垂直梯度；

z表示深度；

或采用已获得的先验速度作为初始速度；

S4:与此同时，基于激发点和接收点位置构建地形；

步骤四：射线段参数的初始化及其优化反演：

S5:设定初至波射线段参数表示为：

反射波射线段参数表示为：

其中，

和分别为初至波回折点的位置和反射点的位置；

和分别为初至波回折点c到激发点S、接收点R的单程走时和反射点C到激发点S、接收点R的单程走时，为反射射线在点C朝着激发点S和接收点R的出射角；

S6:取单程走时等于双程走时的一半，回折点c和反射点C的横坐标近似为炮点与检波点的中心位置，纵坐标选择是通过的给定一系列深度值进行正演，选取正演结果最符合实际数据的深度值作为初至值，射线从反射点出发的初始角度由反射点与炮点、检波点的位置决定；

S7:初始化后，基于初始速度，采用阻尼SVD法进一步优化射线段参数；

步骤五：速度模型和射线段参数的联合反演：

S8:对模型参数进行归一化，得到

将联合层析成像的目标函数设为下式：

其中，||.||为二范数运算，φ为权重，取值范围在0到1之间，C_d为数据协方差矩阵，d_{fst_cal}表示正演初至数据，d_{ref_cal}表示正演反射数据，λ为阻尼系数，L为速度扰动的拉普拉斯运算符，n_{v_prior}为先验的速度模型，添加该项用于抑制反演过程中的多解性问题；

S9:采用阻尼最小二乘方法原理，将以上非线性问题线性化，以迭代的形式进行求解；并设n_{v_prior}为上次反演得到的结果，可得第k次迭代的求解的方程组如下：

其中，

表示第k次迭代初至波数据空间对模型空间的敏感度矩阵；

表示第k次迭代反射波数据空间对模型空间的敏感度矩阵；

为第k次迭代待求的未知量；

为第k次迭代的正演初至波数据和拾取初至波数据之差；

为第k次迭代的正演反射波数据和拾取反射波数据之差；

S10:以上公式中的敏感度矩阵可由傍轴射线近似公式计算得到；

S11:构成得线性方程组采用LSQR算法求解；

S12:求得以上方程组后可进一步计算出模型更新量；

步骤六：检查反演结果：反演之后要及时检查反演结果是否收敛、合理，若反演结果合适，则反演结束，否则反演继续。

步骤六的反演过程采用了多尺度策略，以达到较好的模型分辨率和反演精度；此外，反演后用于对反演所用数据进行筛选，对一些误差较大的数据进行剔除或重新赋值处理。

该发明的处理流程如图1所示。首先是对来自起伏地区的陆上地震数据进行基本的预处理包括滤波、增益等，之后提取和整理得到反演所需的初至波走时斜率信息和反射波走时斜率信息，初始化速度模型和射线段参数，之后应用SVD算法对射线段参数进行优化，之后进行射线段和速度参数的联合反演，得到速度模型，之后选择是否剖分速度模型，若是则剖分速度网格，进而进行联合反演，若否则输出最终结果。

以中国西北某山前带的实际数据反演结果为例，说明实际应用效果。该工区范围长约20km，宽约1.5km，最大高程差约0.3km，地势大致为西低东高，工区的左侧部分为戈壁地貌，工区的右侧部分为山前带。实际资料原为三维地震数据，现抽取其中的一条测线进行处理。选取的资料共66炮，为单边接收方式，炮间距约为200m，每炮有120接收道，道间距约为50m。受地表情况的影响，炮点、检波点分布不均匀。图2(a)显示的是AGC(自动增益控制)后的某一单炮数据。由图2(a)可知，初至波波形清晰，便于初至数据拾取，最终得到7734组初至走时斜率数据。另一方面，地震记录中，噪声能量强，需要作去噪预处理。预处理后信噪比有所提高，该反射轴更为清晰，但资料依然存在相当强烈的噪声干扰，尤其是大于3s的部分，如图2(b)所示，所以在设置质量控制参数时需要小心。拾取的数据通过质量控制，剔除了一些异常值，最终选取了16119组反射数据作为层析反演的输入数据。

应用初至波反射波联合斜率层析成像方法解决地下速度建模问题，反演采用的初始速度为(2.5+0.5z)km/s，如图3所示，共有13×26个离散速度节点，节点的横向间距和纵向间距均为1km。初步反演之后，对离散速度模型进行剖分，最终速度节点的横向间距和纵向间距为0.25km。在联合斜率层析中，初至权重设为0.01，反射权重设为0.99，最终得到的反演速度如图4(a)所示。图4(b)是联合斜率层析的深度偏移剖面，剖面中反射同相轴清晰可见。

为对比反演效果，本文选取横向4km～9km范围内的部分剖面进行分析，并应用偏移速度分析得到的结果以进行比较。图5(a)显示的是偏移速度分析得到的速度，图6(a)显示的是联合斜率层析得到的速度，对比结果，联合斜率层析成像的反演结果存在较明显的浅层速度横向变化。

图5(b)为基于偏移速度分析结果的深度偏移剖面，图5(c)为其在x＝6.5km处共成像道集。图6(b)与图6(c)分别为联合斜率层析的深度偏移剖面及其在x＝6.5km处的共成像道集。对比图5(b)、图6(b)中显示的偏移剖面，对比黑圈中的部分，联合斜率层析成像在标志同相轴处的连续性更好一些，在其之下的同相轴清晰度也明显好于偏移速度分析的结果。对比图5(c)、图6(c)中显示的共成像道集，黑圈中的部分与标志同相轴对应，对比可以发现，基于联合斜率层析的CIG同相轴得到拉平，并且清晰、连续，而基于偏移速度分析的CIG同相轴虽然也被拉平，但其连续性明显差一些。

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

Claims (4)

1.一种地震初至波和反射波联合斜率层析成像方法，其特征在于，包括如下步骤：

步骤一：地震数据的预处理；

步骤二：斜率和走时数据拾取与质量控制：

S1:斜率拾取是在共炮点道集和共检波点道集中实现的，采用对地震道包络进行倾斜叠加的方式，获得最大叠加能量对应的斜率；

S2:此后，通过人机交互方式进行质量控制，最终获得：

初至波数据：

和反射波数据：

其中，

分别为激发点S的横纵坐标以及初至波斜率；

分别为接收点R的横纵坐标以及初至波斜率；

双程初至走时为

N₁为初至波数据个数；

分别为激发点S的横纵坐标以及反射波斜率；

分别为接收点R的横纵坐标以及反射波斜率；

双程反射走时为

N₂为反射波数据个数；

步骤三：初始化速度模型：

S3:速度模型参数m_v的初始值采用随深度线性增大的函数v＝(v₀+k·z)km/s；

其中，

v₀为地表速度；

k为速度垂直梯度；

z表示深度；

或采用已获得的先验速度作为初始速度；

S4:与此同时，基于激发点和接收点位置构建地形；

步骤四：射线段参数的初始化及其优化反演：

S5:设定初至波射线段参数表示为：

反射波射线段参数表示为：

其中，

和分别为初至波回折点的位置和反射点的位置；

和分别为初至波回折点c到激发点S、接收点R的单程走时和反射点C到激发点S、接收点R的单程走时，为反射射线在点C朝着激发点S和接收点R的出射角；

S6:取单程走时等于双程走时的一半，回折点c和反射点C的横坐标近似为炮点与检波点的中心位置，纵坐标选择是通过的给定一系列深度值进行正演，选取正演结果最符合实际数据的深度值作为初至值，射线从反射点出发的初始角度由反射点与炮点、检波点的位置决定；

S7:初始化后，基于初始速度，采用阻尼SVD法进一步优化射线段参数；

步骤五：速度模型和射线段参数的联合反演：

S8:对模型参数进行归一化，得到

将联合层析成像的目标函数设为下式：

其中，||.||为二范数运算，φ为权重，取值范围在0到1之间，C_d为数据协方差矩阵，d_{fst_cal}表示正演初至数据，d_{ref_cal}表示正演反射数据，λ为阻尼系数，L为速度扰动的拉普拉斯运算符，n_{v_prior}为先验的速度模型，添加该项用于抑制反演过程中的多解性问题；

S9:采用阻尼最小二乘方法原理，将以上非线性问题线性化，以迭代的形式进行求解；并设n_{v_prior}为上次反演得到的结果，可得第k次迭代的求解的方程组如下：

其中，

表示第k次迭代初至波数据空间对模型空间的敏感度矩阵；

表示第k次迭代反射波数据空间对模型空间的敏感度矩阵；

为第k次迭代待求的未知量；

为第k次迭代的正演初至波数据和拾取初至波数据之差；

为第k次迭代的正演反射波数据和拾取反射波数据之差；

S10:以上公式中的敏感度矩阵可由傍轴射线近似公式计算得到；

S11:构成得线性方程组采用LSQR算法求解；

S12:求得以上方程组后可进一步计算出模型更新量；

步骤六：检查反演结果：反演之后要及时检查反演结果是否收敛、合理，若反演结果合适，则反演结束，否则反演继续。

2.根据权利要求1所述的地震初至波和反射波联合斜率层析成像方法，其特征在于，所述步骤二的斜率和走时数据拾取与质量控制，斜率拾取是在共炮点道集和共检波点道集中实现的，采用对地震道包络进行倾斜叠加的方式，获得最大叠加能量对应的斜率。

3.根据权利要求1所述的地震初至波和反射波联合斜率层析成像方法，其特征在于，所述步骤五的速度模型和射线段参数的联合反演，通过同时应用初至波斜率数据和反射波斜率数据，构建联合层析成像的目标函数，通过求解局部优化的反演方程组，更新模型参数。

4.根据权利要求1所述的地震初至波和反射波联合斜率层析成像方法，其特征在于，所述步骤六的反演过程采用了多尺度策略，以达到需要的模型分辨率和反演精度。