CN101609167B - 基于起伏地表的井间地震波动方程叠前深度偏移成像方法 - Google Patents

Abstract

本发明公开了一种基于起伏地表的井间地震波动方程叠前深度偏移成像方法，包括：A、预设基准面，将反射波场进行初始化后得到初始反射波场；B、对所述初始反射波场使用波动方程延拓算子计算出波场值；C、在某一延拓步长上，如果检测到炮点或检波器，则计算出炮点或者检波器对应的横向位置的波场值；D、将成像空间范围内所有的炮点或检波器重复步骤C，得到所有炮点或者检波器对应的横向位置的波场值；E、对步骤D得到的所有波场值都进行处理后进行叠加成像。本发明方法不仅能够适用于复杂地质构造和横向变速介质，而且具有成像精度高，井间地震反射波振幅保持好、计算效率较高的优点。

Description

基于起伏地表的井间地震波动方程叠前深度偏移成像方法

技术领域

本发明属于井间地震资料处理领域，尤其涉及一种基于起伏地表的井间地震波动方程叠前深度偏移成像方法。

背景技术

井间地震反射波资料处理主要由预处理、反射波场分离和反射波成像三部分组成。反射波成像则是利用走时场延拓成像和偏移成像、VSP-CDP成像等技术，将分离出的上、下行反射波转换成类似于地面地震资料的井间地震反射波剖面。在井间地震反射波成像中，VSP-CDP成图方法是通常采用的方法，VSP-CDP成图的优点是算法稳定、容易实现，但其缺点也很明显，主要原因是该方法基于水平层状常速介质假设，因而成像精度较低。井间地震的POSTMAP成像是在VSP-CDP成像结果基础上，再进一步使绕射波收敛归位，其成像质量稍有提高，但是难以适应速度场的强横向变化。

叠前深度偏移成像技术按方法可分为Kirchhoff积分法和波动方程延拓法。Kirchhoff积分法根据射线追踪等获得偏移成像所需的走时信息，可以对陡倾角反射层进行成像，但在复杂地质条件下，地震波的多路径问题难以解决，在地震波振幅保持成像方面难度也较大。

发明内容

本发明实施例的目的在于，提供一种基于起伏地表的井间地震波动方程叠前深度偏移成像方法，用以解决现有偏移成像方法中存在的难以适应速度场横向变化、成像精度偏低及井间地震反射波振幅保持不好的问题。

为实现上述目的，本发明实施例提供一种基于起伏地表的井间地震波动方程叠前深度偏移成像方法，该方法包括：

A、预设基准面，将反射波场进行初始化后得到初始反射波场；

B、对所述初始反射波场使用波动方程延拓算子计算出波场值；

C、在某一延拓步长上，如果检测到炮点或检波器，则计算出炮点或者检波器对应的横向位置的波场值；

D、将成像空间范围内所有的炮点或检波器重复步骤C，得到所有的炮点或检波器的波场值；

E、对步骤D得到的所有波场值都进行处理后进行叠加成像。

其中，所述反射波场可为频率-空间域下行反射波场和上行反射波场。

所述波动方程延拓算子可为：

$\frac{\∂\tilde{u}}{\∂z}=\±\frac{\mathrm{j\ω}}{v(x,y,z)}[Q+E]\tilde{u}$

其中

为频率-空间域反射波场；v(x，y，z)为介质速度；ω为频率；j为虚数单位；±符号分别对应频率-空间域下行反射波场和上行反射波场； $Q=1+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_x}{1+{\mathrm{\β}}_i{R}_x}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_y}{1+{\mathrm{\β}}_i{R}_y},$ 其中，α_i、β_i是连分式展开系数； $E=\sqrt{1-{\left(\frac{v{k}_x}{\mathrm{\ω}}\right)}^2-{\left(\frac{v{k}_y}{\mathrm{\ω}}\right)}^2}-[1-\frac{\mathrm{\α}{\left(\frac{v{k}_x}{\mathrm{\ω}}\right)}^2}{1-\mathrm{\β}{\left(\frac{v{k}_x}{\mathrm{\ω}}\right)}^2}-\frac{\mathrm{\α}{\left(\frac{v{k}_y}{\mathrm{\ω}}\right)}^2}{1-\mathrm{\β}{\left(\frac{{\mathrm{vk}}_y}{\mathrm{\ω}}\right)}^2}],$ k_x、k_y分别是频率-波数域x和y方向的波数。

步骤C所述计算出炮点或者检波器对应的横向位置的波场值具体可为：将步骤B计算出的所述波场值加上所述某一延拓步长位置处的检波器对应的地震道上得到的波场值或者炮点对应的地震道上得到的波场值。

步骤E所述对步骤D得到的所有波场值都进行处理后进行叠加成像可为：

E1、根据相关成像条件提取每个波场值的成像值，并将每个成像值放置到对应的横向位置的波场中，得到成像剖面图；

E2、对每个炮点或者检波器重复步骤A、B、C、D及E1直至得到炮点或者检波器的多个成像剖面图；

E3、对所述多个成像剖面图进行叠加成像。

所述延拓步长可为不大于检波器间距或者炮点间距的最小值。

步骤A之前还可包括：

A1、建立炮点和检波器的深度点。

本发明实施例具有以下有益效果：该方法既考虑了井间地震反射波资料的运动学特征(时间信息等)，又考虑了井间地震反射波资料的动力学特征(振幅信息等)，因而能够适用于复杂地质构造和横向变速介质，成像精度高，井间地震反射波振幅保持好，并且计算效率也较高。

附图说明

图1为本发明实施例提供的基于起伏地表的井间地震波动方程叠前深度偏移成像方法示意图；

图2为本发明实施例提供的上行反射波的炮域波动方程叠前深度偏移成像方法示意图；

图3为胜利油田某一地区地震地质模型；

图4为对图3中的地质模型使用带误差补偿的频率-空间域有限差分波场延拓算子对上行反射波场进行延拓的成像示意图。

具体实施方式

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

本发明的叠前深度偏移成像方法，将用于起伏地表的“逐步-累加”波场延拓法用到井间地震波动方程叠前深度偏移成像方法中，这样由于本发明的方法既考虑了井间地震反射波资料的运动学特征(时间信息等)，又考虑了井间地震反射波资料的动力学特征(振幅信息等)，因而能够适用于复杂地质构造和横向变速介质，成像精度高，井间地震反射波振幅保持好；并且本发明还对传统的傅里叶有限差分波场延拓算子进行改进得到带误差补偿的频率-空间域有限差分波场延拓算子，该带误差补偿的频率-空间域有限差分波场延拓算子相对于传统的算子具有计算效率较高的优点。

用于起伏地表的“逐步-累加”波场延拓法原理为：

将初始反射波场从一个设定的水平基准面向下延拓，在出现检波器(或炮点)的某一深度延拓步长上，将初始反射波场根据波动方程延拓算子计算出的反射波场值加上该延拓步长位置处检波器(或炮点)对应的地震道上的波场值得到该位置处的反射波场值，这样就可以对起伏地形上的记录进行叠前偏移。

图1为本发明实施例提供的基于起伏地表的井间地震波动方程叠前深度偏移成像方法示意图，该方法包括以下步骤：

S11、预设基准面，将反射波场进行初始化后得到初始反射波场；

步骤S11中，所述反射波场为频率-空间域下行反射波场和上行反射波场。并且所述基准面的设置与起伏地表的“逐步-累加”波场延拓法中基准面的设置相同或者相类似，本领域技术人员完全可以明白如何来预设基准面，在此不再赘述。

S12、对所述初始反射波场使用波动方程延拓算子计算出波场值；

步骤S12中，波动方程延拓算子可以采用传统的分步傅里叶算子、傅里叶有限差分波场延拓算子等来计算，当然为了提高计算效率，也可以采用经过本发明改进的带误差补偿的频率-空间域有限差分波场延拓算子来计算，该带误差补偿的频率-空间域有限差分波场延拓算子形式为：

$\frac{\∂\tilde{u}}{\∂z}=\±\frac{\mathrm{j\ω}}{v(x,y,z)}[Q+E]\tilde{u},$ 其中，

为频率-空间域反射波场；v(x，y，z)为介质速度；ω为频率；j为虚数单位；±符号分别对应频率-空间域下行反射波场和上行反射波场；

$Q=1+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_x}{1+{\mathrm{\β}}_i{R}_x}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_y}{1+{\mathrm{\β}}_i{R}_y},$ α_i、β_i是连分式展开系数；

$E=\sqrt{1-{\left(\frac{v{k}_x}{\mathrm{\ω}}\right)}^2-{\left(\frac{v{k}_y}{\mathrm{\ω}}\right)}^2}-[1-\frac{\mathrm{\α}{\left(\frac{v{k}_x}{\mathrm{\ω}}\right)}^2}{1-\mathrm{\β}{\left(\frac{v{k}_x}{\mathrm{\ω}}\right)}^2}-\frac{\mathrm{\α}{\left(\frac{v{k}_y}{\mathrm{\ω}}\right)}^2}{1-\mathrm{\β}{\left(\frac{{\mathrm{vk}}_y}{\mathrm{\ω}}\right)}^2}],$ k_x、k_y分别是频率-波数域x和y方向的波数。

S13、在某一延拓步长上，如果检测到炮点或检波器，则计算出炮点或者检波器对应的横向位置的波场值；

步骤S13中，为了使波场值能准确的对应到横向位置(该横向位置是指成像点对应的所在速度场中的地面水平位置)，所述延拓步长可为不大于检波器间距或者炮点间距的最小值，而如果所述延拓步长不能满足上述要求时，则对速度场进行插值即可。计算出炮点或者检波器对应的横向位置的波场值可为：将步骤S12计算出的所述波场值加上所述某一延拓步长位置处的检波器对应的地震道上得到的波场值或者炮点对应的地震道上得到的波场值即可。

S14、对成像空间范围内所有的炮点或检波器重复步骤S13，得到所有的炮点或检波器的波场值；

S15、对步骤S14得到的所有波场值都进行处理后进行叠加成像。

步骤S15中，所述对步骤S14得到的所有波场值都进行处理后进行叠加成像具体为：

S151、根据相关成像条件提取每个波场值的成像值，并将每个成像值放置到对应的横向位置的波场中，得到成像剖面图；

例如，在频率域，将震源波场和检波波场相乘，并把所有频率成分相加提取成像值(即不同频率成分加和后的振幅值)放到对应的横向位置(即成像点对应的所在速度场中的地面水平位置)即可得到成像剖面图。

S152、对每个炮点或者检波器重复步骤S11-S14、S151直至得到所有炮点或者检波器的多个成像剖面图；

S153、对所述多个成像剖面图进行叠加成像。

在上述成像方法过程中，步骤S11之前还可以包括建立炮点及检波器深度点，并且炮点及检波器对应的深度点位置可以从道头字的震源、检波器的坐标中读出深度值。

下面以上行反射波的炮域波动方程叠前深度偏移成像过程为例，对本发明成像过程进行详细说明：

对于一个炮点来说，在井中的垂直深度上有多个检波器用于接收该炮点激发产生的上行反射波场，并得到很多的地震道。每个检波器位置处的波场值是已知的，炮点处为震源波场，如图2所示，则上行反射波的炮域波动方程叠前深度偏移成像步骤为：

S21、建立每个炮点和检波器的深度点；

S22、初始化上行反射波场，得到初始上行反射波场；

S23、对所述初始上行反射波场利用带误差补偿的频率-空间域有限差分波场延拓算子进行延拓计算出波场值；

S24、在某一延拓步长上，如果检测到检波器，则同时根据所述S21建立的炮点和检波器的深度点来判断在该延拓步长的位置深度是否与选定的一个炮点和检波器的位置深度一致，如果是则进行步骤S25，否则使得S23延拓计算出的波场值即为检波器对应的横向位置的波场值后进行步骤S26；

步骤S24中，如果检测不到检波器，则继续进行下一个延拓步长，直到检测到检波器。

S25、将S23延拓计算出的波场值加上该延拓步长的位置深度处的检波器对应的地震道上得到的波场值，得到检波器对应的横向位置的波场值；

S26、将检波器对应的横向位置的波场值和震源波场相乘，并把所有频率成分的成像值相加(即相加不同频率成分的振幅值)，然后放到对应的横向位置即成像点对应的所在速度场中的地面水平位置；

S27、对井中的垂直深度上的多个检波器重复步骤S23-S26，直至所有的检波器重复完毕；

S28、得到某个炮点的成像剖面图；

S29、对下一个炮点重复步骤S21-S28直至所有的炮点重复完毕，从而得到多个成像剖面图；

S30、对所有的成像剖面图进行叠加成像。

这样，上行反射波的炮域波动方程叠前深度偏移成像过程完成。

同样，对于一个检波器来说，会对应很多个炮点，其波动方程叠前深度偏移成像过程与上述炮域的类似，在此不再赘述，在此仅强调一下，在检测到炮点时，该炮点的深度与延拓步长的位置深度及选定的一个检波器的位置深度一致时，该炮点的位置深度处的反射波场值为：延拓计算出的波场值加上该炮点的位置深度处从炮点对应的地震道上得到的波场值。

下面将推导带误差补偿的频率-空间域有限差分波场延拓算子过程进行详细描述：

以空间域的三维声波方程出发：

$\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2}+\frac{{\∂}^{2}u}{\∂z^2}=\frac{1}{v^2(x,y,z)}\frac{{\∂}^{2}u}{\∂t^2}---\left(1\right)$

其中，u(x，y，z)为空间域的反射波场，v(x，y，z)为介质的速度。

在频率-空间域中，式(1)可以表示为

$\frac{\∂\tilde{u}}{\∂z}=\±\frac{\mathrm{j\ω}}{v(x,y,z)}[1+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_x}{1+{\mathrm{\β}}_i{R}_x}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_y}{1+{\mathrm{\β}}_i{R}_y}]]\tilde{u}---\left(2\right)$

其中， $R_x=\frac{v^2(x,y,z)}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2},$ $R_y=\frac{v^2(x,y,z)}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂y^2},$

为频率-空间域反射波场，v(x，y，z)为介质速度，ω为频率，j为虚数单位，±符号分别对应频率-空间域下行反射波场和上行反射波场，α_i，β_i是连分式展开系数。

在频率-空间域中上述(2)式还可以表示为：

$\frac{\∂\tilde{u}}{\∂z}=\±\frac{\mathrm{j\ω}}{v(x,y,z)}\sqrt{1+R_x+R_y}\tilde{u}=\±\frac{\mathrm{j\ω}}{v(x,y,z)}P\tilde{u}---\left(3\right)$

其中， $P=\sqrt{1+R_x+R_y}$

由(2)式可以得出：

$\frac{\∂\tilde{u}}{\∂z}=\±\frac{\mathrm{j\ω}}{v(x,y,z)}Q\tilde{u}---\left(4\right)$

其中， $Q=1+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_x}{1+{\mathrm{\β}}_i{R}_x}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_y}{1+{\mathrm{\β}}_i{R}_y}$

则P和Q两者之间的差分算子误差E为

$E=P-Q=\sqrt{1+R_x+R_y}-(1+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_x}{1+{\mathrm{\β}}_i{R}_x}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_y}{1+{\mathrm{\β}}_i{R}_y})---\left(5\right)$

当(5)式中的n＝1时，在频率-波数域，则差分算子误差E可以表示为：

$E=\sqrt{1-{\left(\frac{{\mathrm{vk}}_x}{\mathrm{\ω}}\right)}^2-{\left(\frac{v{k}_y}{\mathrm{\ω}}\right)}^2}-[1-\frac{\mathrm{\α}{\left[\frac{{\mathrm{vk}}_x}{\mathrm{\ω}}\right]}^2}{1-\mathrm{\β}{\left(\frac{{\mathrm{vk}}_x}{\mathrm{\ω}}\right)}^2}-\frac{\mathrm{\α}{\left(\frac{{\mathrm{vk}}_y}{\mathrm{\ω}}\right)}^2}{1-\mathrm{\β}{\left(\frac{v{k}_y}{\mathrm{\ω}}\right)}^2}]---\left(6\right)$

其中，k_x、k_y是空间x和y方向的波数。

差分算子误差E可以按(6)式在一延拓步长或若干延拓步长上进行相移校正。

因此，带误差补偿的频率-空间域有限差分波场延拓算子可以表示为：

$\frac{\∂\tilde{u}}{\∂z}=\±\frac{\mathrm{j\ω}}{v(x,y,z)}[Q+E]\tilde{u}---\left(7\right)$

本发明实施例提供的方法具有的有益效果：由于该方法既考虑了井间地震反射波资料的运动学特征(时间信息等)，又考虑了井间地震反射波资料的动力学特征(振幅信息等)，因而能够适用于复杂地质构造和横向变速介质，成像精度高，井间地震反射波振幅保持好；并且本发明实施例还对传统的傅里叶有限差分波场延拓算子进行改进得到带误差补偿的频率-空间域有限差分波场延拓算子，该带误差补偿的频率-空间域有限差分波场延拓算子相对于传统的算子具有计算效率较高，从而最终能够提高成像精度的优点。

图3为胜利油田某一地区地震地质模型；图4为对图3中的地质模型使用带误差补偿的频率-空间域有限差分波场延拓算子对上行反射波场进行延拓的成像示意图。如图3和图4所示，使用本发明实施例提供的成像方法对胜利油田某一地区进行了成功的井间成像。从图3中可以看到两井之间存在两个较大断距的断层、多个透镜体、砂体尖灭、砂泥岩薄互层和倾斜产状的地层等。而图4中可以看到对图3中所示的两个较大的断层和砂泥岩薄互层以及倾斜产状的地层都有很好的成像，成像效果较好。

显然，上述实施例仅是本发明优选实施例而已，它并不限制本发明的保护范围，在本发明的保护范围内，所述领域技术人员还可以对本发明的方法做出各种改进及润饰，当然这些改进及润饰也视为本发明的保护范围。

Claims (4)

1.一种基于起伏地表的井间地震波动方程叠前深度偏移成像方法，其特征在于，该方法包括以下步骤：

A、预设基准面，将反射波场进行初始化后得到初始反射波场；

B、对所述初始反射波场使用带误差补偿的频率-空间域有限差分波场延拓算子计算波场值，所述带误差补偿的频率-空间域有限差分波场延拓算子具体为：

$\frac{\∂\tilde{u}}{\∂z}=\±\frac{\mathrm{j\ω}}{v(x,y,z)}[Q+E]\tilde{u},$

其中：

为频率-空间域反射波场；v(x，y，z)为介质速度；ω为频率；j为虚数单位；±符号分别对应频率-空间域下行反射波场和上行反射波场； $Q=1+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_x}{1+{\mathrm{\β}}_i{R}_x}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_i{R}_y}{1+{\mathrm{\β}}_i{R}_y},$ 其中， $R_x=\frac{{\mathrm{\ν}}^2(x,y,z)}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2},$ $R_y=\frac{{\mathrm{\ν}}^2(x,y,z)}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂y^2},$ α_i、β_i是连分式展开系数； $E=\sqrt{1-{\left(\frac{v{k}_x}{\mathrm{\ω}}\right)}^2-{\left(\frac{v{k}_y}{\mathrm{\ω}}\right)}^2}-[1-\frac{\mathrm{\α}{\left(\frac{v{k}_x}{\mathrm{\ω}}\right)}^2}{1-\mathrm{\β}{\left(\frac{v{k}_x}{\mathrm{\ω}}\right)}^2}-\frac{\mathrm{\α}{\left(\frac{v{k}_y}{\mathrm{\ω}}\right)}^2}{1-\mathrm{\β}{\left(\frac{v{k}_y}{\mathrm{\ω}}\right)}^2}],$ k_x、k_y分别是频率-波数域x和y方向的波数；

C、在某一延拓步长上，如果检测到炮点或检波器，则计算出炮点或者检波器对应的横向位置的波场值，所述计算出炮点或者检波器对应的横向位置的波场值具体为：将步骤B计算出的所述波场值加上所述某一延拓步长位置处的检波器对应的地震道上得到的波场值或者炮点对应的地震道上得到的波场值；

D、对成像空间范围内所有的炮点或检波器重复步骤C，得到所有炮点或者检波器对应的横向位置的波场值；

E、对步骤D得到的所有波场值都进行处理后进行叠加成像，具体为：

E1、根据相关成像条件提取每个波场值的成像值，并将每个成像值放置到对应的横向位置的波场中，得到成像剖面图；

E2、对每个炮点或者检波器重复步骤A、B、C、D及E1直至得到所有炮点或者检波器的多个成像剖面图；

E3、对所述多个成像剖面图进行叠加成像。

2.根据权利要求1所述的方法，其特征在于，所述反射波场为频率-空间域下行反射波场和上行反射波场。

3.根据权利要求1所述的方法，其特征在于，所述延拓步长为不大于检波器间距或者炮点间距的最小值。

4.根据权利要求1所述的方法，其特征在于，步骤A之前还包括：

A1、建立炮点和检波器的深度点。

Images