CN102466818B - 一种利用井间地震数据对各向异性介质成像的方法 - Google Patents
一种利用井间地震数据对各向异性介质成像的方法 Download PDFInfo
- Publication number
- CN102466818B CN102466818B CN201010543171.0A CN201010543171A CN102466818B CN 102466818 B CN102466818 B CN 102466818B CN 201010543171 A CN201010543171 A CN 201010543171A CN 102466818 B CN102466818 B CN 102466818B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- travel time
- seismic
- wave
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 39
- 238000000034 method Methods 0.000 title claims abstract description 28
- 238000000926 separation method Methods 0.000 claims abstract description 13
- 238000004364 calculation method Methods 0.000 claims description 15
- 230000010354 integration Effects 0.000 claims description 9
- 238000001914 filtration Methods 0.000 claims description 6
- 239000004576 sand Substances 0.000 claims description 3
- 230000015572 biosynthetic process Effects 0.000 abstract 2
- 230000000704 physical effect Effects 0.000 description 3
- 230000005284 excitation Effects 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及物探技术种利用井间地震反射数据对各向异性介质成像的方法。利用采集的数据反演两口井之间的地质模型获得各向异性系数和建立的层状地质模型,将层状地质模型进行网格剖分,计算每个网格上随射线传播方向变化的射线速度,逐一计算震源点所在列网和邻近格节点上的旅行时间,形成炮点和接收点走时文件,波场分离,将井间地震地震记录分离成不同类型的波场,利用积分法完成各向异性介质成像。本发明保证了成像方法具有很好的模型适应性,并有很高的成像精度。
Description
技术领域
本发明涉及地球物理勘探技术,是地震资料处理方法,具体是一种利用井间地震反射数据对各向异性介质成像的方法。
背景技术
井间地震是在一口井中激发,在另外一口井或多口井中接收的地震勘探方法,由于是在井中激发和接收,激发能量传播距离短,传播路径避开低速带,观测接近探测目标,因此井间地震数据具有很高的频率和信噪比。随着油气资源需求与供给矛盾的日益严峻以及井中地震采集设备性能的提高,井间地震以其高精度和高分辨率的显著特点,正在成为解决构造精细成像、储层描述、油藏动态监测、了解剩余油分布等问题的一种关键技术。
井间地震成像处理包括初至波层析成像和反射波成像,常用的井间地震反射成像方法主要有:(1)XSP-CDP转换算法(Lazaratos,Khalil,1993),该方法利用声波测井或层析成像提供的速度模型,采用射线追踪的算法,对每一个共炮点道集或共接收点道集单独直接进行归位成像,速度模型的精度直接影响成像效果。(2)共中心深度点(CMD)叠加法(stewart,1993),该方法利用共中心深度点道内的各反射波旅行时相等的特点抽成CMD道集,然后计算道集内的各道所对应的反射点位置进行叠加成像。这种方法使用了井间的多道集叠加,提高了抗噪能力。该方法是在水平反射界面的假设下提出来的,只适用于反射界面倾角很小,地层速度变化不是很剧烈的情况。(3)共侧向点(CLP)叠加法(Smalley,1992,这种方法是从常速介质共炮点道集满足双曲线时距方程的前提出发,推导出水平层状均匀介质、倾斜层均匀介质的时距方程,通过水平和垂直的动校正和叠加速度分析进行成像,得到井间地震反射波时间成像剖面。这种方法的成像结果为时间剖面,并难以适应横向速度剧烈变化的情况。吴律公开了适应倾斜层的井间地震共深度点(DLCDP)成像算法,严又生等公开了适应非均匀介质的井间地震反射波的VSP-CDP成像方法,还有学者利用波动方程逆时偏移的方法对井间地震反射数据进行成像研究。
地球上的大多数岩石是各向异性的,它们的物理特性随方向而变化。最简单的各向异性是横向各向同性,就是在水平方向上物性分布均匀,在垂直方向上物性变化。然而,现有的井间地震反射数据成像方法,大都是基于各向同性介质模型的成像方法。但是,由于井间地震特殊的观测方式,通常接收来自各个方向的射线,观测到的井间地震数据各向异性特征非常明显,并且井间地震数据的频率远远高于地面地震,要求更准确的成像时间和成像方法,所以利用常规的基于各向同性模型的成像方法对井间地震反射数据进行成像会给成像结果较大误差,不能满足井间地震高分辨率的需要。
发明内容
本发明目的是提供一种模型适应性好,提高精度与分辨率的利用井间地震反射波数据对各向异性介质成像的方法。
本发明通过以下技术方案实现,具体步骤是:
1)采集井间地震数据;
步骤1)所述的采集数据包括采集两口井的零偏VSP地震记录,多分量井间地震记录和密度测井数据。
2)利用采集的数据反演两口井之间的地质模型,包括:
利用零偏VSP地震记录的纵波初至时间TP反演模型纵波的垂直速度VP(0°),
利用零偏VSP地震记录的横波初至时间Ts反演模型横波的垂直速度VS(0°),
利用密度测井数据获得模型密度ρ,
利用井间地震纵波初至时间数据获得纵波的水平速度VP(90°)和P波45度速度VP(45°),
然后以下公式(1)和(2)获得各向异性系数ε和δ建立的层状地质模型:
3)将层状地质模型进行网格剖分;
步骤3)的网格剖分是每个网格包含的模型参数有纵波的垂直速度VP(0°),横波的垂直速度VS(0°),各向异性系数ε和δ,介质密度ρ。
4)用下式计算每个网格上随射线传播方向ψ变化的射线速度V(ψ):
V-2(ψ)≈a1+a2cos2(ψ)-a3cos4(ψ)(3)
其中:
α0=VP(0°)
ε和δ为各向异性系数
Vg[ψ(45°)]和ψ(45°)由公式(5)至(8)求出
5)计算震源点所在列网格节点上的旅行时间;
步骤5)首先确定震源点所在的网格结点(is,js),再利用纵波垂直速度VP(0°)和公式(9)计算震源点所在列的网格结点上的旅行时Ts(is,j),j=1,2,3Λjs-1,js+1,ΛN,N为网格模型的行数。
6)计算震源点所在列右边网格节点上的旅行时;
步骤6)所述的计算震源点所在列右边网格节点上的旅行时是将步骤5)计算出的震源点所在列节点上的旅行时间Ts(is,j)作为初始时间,计算震源点所在列右边网格节点上的旅行时,计算顺序按列进行,逐列向右递推,直到计算区域的右边界。
7)计算震源点所在列左边网格节点上的旅行时;
8)由上到下按行计算炮点到每个网格结点上的旅行时;
9)由下到上按行计算炮点每个网格结点上的旅行时;
10)形成炮点走时文件;
步骤10)所述的形成炮点走时文件是将步骤6)至步骤9)计算得到的网格模型每个网格节点上的旅行时间按照炮号、网格节点横坐标、网格节点纵坐标、旅行时间排列记录到炮点走时文件中。
11)计算接收点到每个网格节点的旅行时;
12)形成接收点走时文件;
13)波场分离,将井间地震地震记录分离成不同类型的波场;
步骤13)所述的波场分离是利用频率波数滤波或中值滤波的方法,对井间地震记录进行波场分离,分离后的波场包括上行反射纵波,上行转换横波,下行反射纵波、下行转换横波,本发明利用到的波场为上行反射纵波和下行反射纵波。
14)利用积分法完成各向异性介质成像。
本发明采用的模型是与实际地质模型更为接近的VTI介质模型,旅行时计算采用的是基于VTI网格模型的球面波近似初至计算公式,成像则是采用了考虑近场项的积分公式,保证了成像方法具有很好的模型适应性,并有很高的成像精度。
附图说明
图1为VTI介质模型及井间地震观测系统示意图;
表1为VTI模型参数表;
图2(a)为利用准纵波波动方程正演模拟的井间地震记录(水平分量);
图2(b)为利用准纵波波动方程正演模拟的井间地震记录(垂直分量);
图3为旅行时计算模式示意图;
图4为VTI模型炮点(左)和第45个接收点(右的)时间场;
图5(a)为波场分离后的上行反射P波(垂直分量);
图5(b)为波场分离后的下行反射P波(垂直分量);
图6为利用各向异性时间场得到的井间地震成像剖面(垂直分量)。
具体实施方式
以下结合附图详细说明本发明。
本发明的具体实施方式为:
1)采集井间地震数据;
步骤1)所述的采集数据包括采集两口井的零偏VSP地震记录,多分量井间地震记录和密度测井数据。根据图1给出的VTI介质模型及井间地震观测系统及表1所示的模型参数,利用准纵波波动方程正演模拟合成井间地震记录,如图2所示。
表1VTI模型参数表
层序号 | vp(0)(m/s) | vs(0)(m/s) | ε | δ | ρ(g/cm3) |
1 | 2745 | 1508 | 0.103 | -0.018 | 2030 |
2 | 3048 | 1490 | 0.255 | -0.050 | 2420 |
3 | 3377 | 1490 | 0.200 | -0.075 | 2420 |
4 | 3794 | 2074 | 0.189 | 0.204 | 2560 |
5 | 4231 | 2539 | 0.200 | 0.100 | 2370 |
5 | 4231 | 2539 | 0.200 | 0.100 | 2370 |
2)利用采集的数据反演两口井之间的地质模型,包括:
利用零偏VSP地震记录的纵波初至时间TP反演模型纵波的垂直速度VP(0°),
利用零偏VSP地震记录的横波初至时间Ts反演模型横波的垂直速度VS(0°),
利用密度测井数据获得模型密度ρ,
利用井间地震纵波初至时间数据获得纵波的水平速度VP(90°)和纵波45度速度VP(45°),
然后以下公式(1)和(2)获得各向异性系数ε和δ建立的层状地质模型:
3)将层状地质模型进行网格剖分;
步骤3)的网格剖分是每个网格包含的模型参数有纵波的垂直速度VP(0°),横波的垂直速度VS(0°),各向异性系数ε和δ,介质密度ρ。
4)用下式计算每个网格上随射线传播方向ψ变化的射线速度V(ψ):
V-2(ψ)≈a1+a2cos2(ψ)-a3cos4(ψ)(3)
其中:
α0=VP(0°)
ε和δ为各向异性系数
Vg[ψ(45°)]和ψ(45°)由公式(5)至(8)求出
5)计算震源点所在列网格节点上的旅行时间;
步骤5)首先确定震源点所在的网格结点(is,js),再利用纵波垂直速度VP(0°)和公式(9)计算震源点所在列的网格结点上的旅行时Ts(is,j),j=1,2,3Λjs-1,js+1,ΛN,N为网格模型的行数。
6)计算震源点所在列右边网格节点上的旅行时;
步骤6)所述的计算震源点所在列右边网格节点上的旅行时是将步骤5)计算出的震源点所在列节点上的旅行时间Ts(is,j)作为初始时间,计算震源点所在列右边网格节点上的旅行时,计算顺序按列进行,逐列向右递推,直到计算区域的右边界。
在每一列的计算过程中,先由上到下,如图3中的模式1所示,设点A纵坐标为z1,B纵坐标为z2,射线与AB交点的纵坐标为z0,利用(10)式,则可以计算出待求点O上的旅行时T1:
其中,变量z0可以利用二分法求解非线性方程(11)式得到:
f(z0)=f1(z0)+f2(z0)+f3(z0)(11)
其中,
φ为射线方向,
VP(φ)是射线方向为φ时的射线速度,
a1、a2、a3为由公式(4)求出的参数;
如果相邻网格的速度大于当前计算网格的速度,计算沿界面传播的滑行波到达O点的旅行时T2,取T1和T2中的最小值作为O点的旅行时Ts(i,j)。按照图3模式1向下逐行递推,直到模型底边界。再由下向上逐行进行计算,如图3模式2所示,将C点和B点的旅行时间分别当作t1和t2,向上逐步递推,直到顶边界,得到相应网格节点上的旅行时间Ts(i,j)。
7)计算震源点所在列左边网格节点上的旅行时;
步骤7)同步骤6),按照图3模式6由上到下计算,直到模型底边界,再按照图3模式5由下到上计算,直到顶边界,并按列逐步向左,计算震源点所在列左边网格节点上的旅行时间Ts(i,j)。
8)由上到下按行计算炮点到每个网格结点上的旅行时;
步骤8)将步骤6)和步骤7)计算和网格节点上的旅行时间作为计算初始时间,在整个模型区域内,如步骤6)所述,先按照图3模式8,先由上到下,再由左到右逐步计算各个网格节点时间,将步骤8)计算的旅行时间与步骤6)和步骤7)计算的旅行时间比较,保留最小的旅行时;类似地,按照图3模式7,先由上到下,再由右到左逐步计算各个网格节点时间,将步骤8)计算的旅行时间与步骤6)和步骤7)计算的旅行时间比较,保留最小的旅行时。
9)由下到上按行计算炮点每个网格结点上的旅行时;
步骤9)与步骤8)类似,按照图3模式3和图3模式4逐行计算每个网格结点上的旅行时。
10)形成炮点走时文件;
步骤10)是将步骤6)至步骤9)计算得到的网格模型每个网格节点上的旅行时间写到炮点走时文件中,顺序为:炮号、网格节点横坐标、网格节点纵坐标、旅行时间。
11)计算接收点到每个网格节点的旅行时;
步骤11)与步骤6)至步骤9)相同,逐步计算各个接收点到网格模型每个网格节点的旅行时间TR(i,j)。
12)形成接收点走时文件;
步骤12)是将步骤11)计算得到的网格模型每个网格节点上的旅行时写到接收点直时文件中,顺序为:接收点号、网格节点横坐标、网格节点纵坐标、旅行时间,图4为计算得到的炮点和接收点的时间场。
13)波场分离,将井间地震地震记录分离成不同类型的波场;
步骤13)是利用频率波数滤波或中值滤波的方法,对井间地震记录进行波场分离,分离后的波场包括上行反射纵波,上行转换横波,下行反射纵波、下行转换横波,本发明利用到的波场为波场分离后的上行反射纵波和下行反射纵波,如图5所示。
14)利用积分法完成各向异性介质成像。
步骤14)所述的利用积分法成像是:对于任一道井间地震记录,对于任一网格节点r(i,j),从炮点走时文件对应坐标位置读出该节点的炮点旅行时间Ts(i,j),从接收点走时文件对应位坐标置读出该节点的接收点旅行时间TR(i,j),利用考虑近了场项的积分公式(12)进行积分法成像:
其中,T=t-Tg,U(r0,t)为在接收点r0上记录的波场,θ为出射线与接收排列法向矢量的夹角。上行反射纵波和下行反射纵波是分别单独成像后再进行合并,对所有的地震道对进行成像之后,形成井间地震成像剖面,实现利用井间地震反射纵数据对各向异性介质成像的目的,图6为利用本发明方法得到的井间地震成像剖面。
Claims (4)
1.一种利用井间地震反射波数据对各向异性介质成像的方法,其特征是,通过以下具体步骤:
1)采集井间地震数据;
2)利用采集的数据反演两口井之间的地质模型,包括:
利用零偏VSP地震记录的纵波初至时间TP反演模型纵波的垂直速度VP(0°),
利用零偏VSP地震记录的横波初至时间Ts反演模型横波的垂直速度VS(0°),
利用密度测井数据获得模型密度ρ,
利用井间地震纵波初至时间数据获得纵波的水平速度VP(90°)和P波45度速度VP(45°),然后以公式(1)和(2)获得各向异性系数ε和δ建立的层状地质模型:
3)将层状地质模型进行网格剖分;
步骤3)的网格剖分是每个网格包含的模型参数有纵波的垂直速度VP(0°),横波的垂直速度VS(0°),各向异性系数ε和δ,模型密度ρ
4)用下式计算每个网格上随射线传播方向ψ变化的射线速度V(ψ):V-2(ψ)≈a1+a2cos2(ψ)-a3cos4(ψ)(3)
其中:
α0=vp(0°)
ε和δ为各向异性系数
Vg[ψ(45°)]和ψ(45°)由公式(5)至(8)求出
5)计算震源点所在列网格节点上的旅行时;
步骤5)首先确定震源点所在的网格结点(is,js),再利用纵波垂直速度VP(0°)和公式(9)计算震源点所在列的网格结点上的旅行时Ts(is,j),j=1,2,3…js-1,js+1,…N,N为网格模型的行数;
6)计算震源点所在列右边网格节点上的旅行时;
7)计算震源点所在列左边网格节点上的旅行时;
8)由上到下按行计算炮点到每个网格结点上的旅行时;
9)由下到上按行计算炮点到每个网格结点上的旅行时;
10)形成炮点走时文件;
11)计算接收点到每个网格节点的旅行时;
12)形成接收点走时文件;
13)波场分离,将井间地震地震记录分离成不同类型的波场;所述的波场分离是利用频率波数滤波或中值滤波的方法,对井间地震记录进行波场分离,分离后的波场包括上行反射纵波、上行转换横波、下行反射纵波、下行反射横波,利用到的波场为上行反射纵波和下行反射纵波;14)利用积分法完成各向异性介质成像;
步骤14)所述的利用积分法成像是:对于任一道井间地震记录,对于任一网格节点r(i,j),从炮点走时文件对应坐标位置读出该节点的炮点旅行时Ts(i,j),从接收点走时文件对应坐标位置读出该节点的接收点旅行时TR(i,j),利用考虑近场项的积分公式(12)进行积分法成像:
其中, T=t-TR,U(r0,t)为在接收点r0上记录的波场,θ为出射线与接收排列法向矢量的夹角;上行反射纵波和下行反射纵波是分别单独成像后再进行合并,对所有的地震道对进行成像之后,形成井间地震成像剖面,实现利用井间地震反射波数据对各向异性介质成像。
2.根据权利要求1所述的方法,其特征是,步骤1)所述的采集数据包括采集两口井的零偏VSP地震记录,多分量井间地震记录和密度测井数据。
3.根据权利要求1所述的方法,其特征是,步骤6)所述的计算震源点所在列右边网格节点上的旅行时是将步骤5)计算出的震源点所在列节点上的旅行时Ts(is,j)作为初始时间,计算震源点所在列右边网格节点上的旅行时,计算顺序按列进行,逐列向右递推,直到计算区域的右边界。
4.根据权利要求1所述的方法,其特征是,步骤10)所述的形成炮点走时文件是将步骤6)至步骤9)计算得到的网格模型每个网格节点上的旅行时按照炮号、网格节点横坐标、网格节点纵坐标、旅行时排列记录到炮点走时文件中。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201010543171.0A CN102466818B (zh) | 2010-11-11 | 2010-11-11 | 一种利用井间地震数据对各向异性介质成像的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201010543171.0A CN102466818B (zh) | 2010-11-11 | 2010-11-11 | 一种利用井间地震数据对各向异性介质成像的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN102466818A CN102466818A (zh) | 2012-05-23 |
CN102466818B true CN102466818B (zh) | 2015-11-25 |
Family
ID=46070708
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201010543171.0A Active CN102466818B (zh) | 2010-11-11 | 2010-11-11 | 一种利用井间地震数据对各向异性介质成像的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN102466818B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103713312B (zh) * | 2012-10-09 | 2016-12-21 | 中国石油化工股份有限公司 | 一种虚源地震观测系统的设计方法 |
CN104484573B (zh) | 2014-12-30 | 2017-09-15 | 中国石油天然气股份有限公司 | 一种确定地层刚性系数的方法 |
CN109725345B (zh) * | 2018-11-15 | 2020-08-11 | 中国石油天然气集团有限公司 | 一种初至波正演模拟方法及装置 |
CN116148925B (zh) * | 2023-04-18 | 2023-06-30 | 山东省科学院海洋仪器仪表研究所 | 一种vti介质球面纵波反射系数解析方法 |
-
2010
- 2010-11-11 CN CN201010543171.0A patent/CN102466818B/zh active Active
Non-Patent Citations (1)
Title |
---|
井间地震交错网格高阶差分数值模拟及逆时偏移成像研究;张文波;《长安大学2005年度博士学位论文》;20061106;63-68,95 * |
Also Published As
Publication number | Publication date |
---|---|
CN102466818A (zh) | 2012-05-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104570125B (zh) | 一种利用井数据提高成像速度模型精度的方法 | |
CN107526101B (zh) | 一种获取地震反射波的采集和处理方法 | |
Regone et al. | Geologic model building in SEAM Phase II—Land seismic challenges | |
CN104133245B (zh) | 一种地震资料的静校正方法及系统 | |
CN101354444B (zh) | 一种确定地层岩性和孔隙流体的方法 | |
CN106094032B (zh) | 一种构建地层速度模型的方法 | |
CN102053263B (zh) | 调查表层结构的方法 | |
CN103645503B (zh) | 一种三维时间域照明分析及振幅补偿方法 | |
CN102841379B (zh) | 一种基于共散射点道集的叠前时间偏移与速度分析方法 | |
CN102466816A (zh) | 一种叠前地震数据地层弹性常数参数反演的方法 | |
CN105629325B (zh) | 前陆盆地冲积扇精细刻画与预测方法 | |
CN103513277B (zh) | 一种地震地层裂隙裂缝密度反演方法及系统 | |
CN103592680B (zh) | 一种基于正反演的测井数据和深度域地震剖面合成方法 | |
WO2012139082A1 (en) | Event selection in the image domain | |
CN103149588B (zh) | 一种利用井震标定计算vti各向异性参数的方法及系统 | |
CN104360388A (zh) | 一种三维地震观测系统评价方法 | |
CN102565852B (zh) | 针对储层含油气性检测的角度域叠前偏移数据处理方法 | |
CN114924315A (zh) | 一种多态式地震勘探方法和系统 | |
CN102466818B (zh) | 一种利用井间地震数据对各向异性介质成像的方法 | |
CN102877828A (zh) | 一种三维多井联合井地ct成像方法 | |
CN104422955B (zh) | 一种利用旅行时变化量进行各向异性参数提取的方法 | |
CN102053269A (zh) | 一种对地震资料中速度分析方法 | |
CN102798888A (zh) | 一种利用非零井源距数据计算纵横波速度比的方法 | |
CN104199088A (zh) | 一种提取入射角道集的方法及系统 | |
CN101937101B (zh) | 一种鉴定能否实施时移地震的方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant |