CN102466818A - 一种利用井间地震数据对各向异性介质成像的方法 - Google Patents

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地震记录的纵波初至时间T_P反演模型纵波的垂直速度V_P(0°)， 

利用零偏VSP地震记录的横波初至时间T_s反演模型横波的垂直速度V_S(0°)， 

利用密度测井数据获得模型密度ρ， 

利用井间地震纵波初至时间数据获得纵波的水平速度V_P(90°)和P波45度速度V_P(45°)， 

然后以下公式(1)和(2)获得各向异性系数ε和δ建立的层状地质模型： 

3)将层状地质模型进行网格剖分； 

步骤3)的网格剖分是每个网格包含的模型参数有纵波的垂直速度V_P(0°)，横波的垂直速度V_S(0°)，各向异性系数ε和δ，介质密度ρ。 

4)用下式计算每个网格上随射线传播方向ψ变化的射线速度V(ψ)： 

V^-2(ψ)≈a₁+a₂cos²(ψ)-a₃cos⁴(ψ)    (3) 

其中： 

α₀＝V_P(0°) 

ε和δ为各向异性系数 

V_S[ψ(45°)]和ψ(45°)由公式(5)至(8)求出 

5)计算震源点所在列网格节点上的旅行时间； 

步骤5)首先确定震源点所在的网格结点(i_s，j_s)，再利用纵波垂直速度V_P(0°)和公式(9)计算震源点所在列的网格结点上的旅行时T_s(i_s，j)，j＝1，2，3Λj_s-1，j_s+1，ΛN，N为网格模型的行数。 

6)计算震源点所在列右边网格节点上的旅行时； 

步骤6)所述的计算震源点所在列右边网格节点上的旅行时是将步骤5)计算出的震源点所在列节点上的旅行时间T_s(i_s，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地震记录的纵波初至时间T_P反演模型纵波的垂直速度V_P(0°)， 

利用零偏VSP地震记录的横波初至时间T_s反演模型横波的垂直速度V_S(0°)， 

利用密度测井数据获得模型密度ρ， 

利用井间地震纵波初至时间数据获得纵波的水平速度V_P(90°)和纵波45度速度V_P(45°)， 

然后以下公式(1)和(2)获得各向异性系数ε和δ建立的层状地质模型： 

3)将层状地质模型进行网格剖分； 

步骤3)的网格剖分是每个网格包含的模型参数有纵波的垂直速度V_P(0°)，横波的垂直速度V_S(0°)，各向异性系数ε和δ，介质密度ρ。 

4)用下式计算每个网格上随射线传播方向ψ变化的射线速度V(ψ)： 

V^-2(ψ)≈a₁+a₂cos²(ψ)-a₃cos⁴(ψ)    (3) 

其中： 

α₀＝V_P(0°) 

ε和δ为各向异性系数 

V_S[ψ(45°)]和ψ(45°)由公式(5)至(8)求出 

5)计算震源点所在列网格节点上的旅行时间； 

步骤5)首先确定震源点所在的网格结点(i_s，j_s)，再利用纵波垂直速度V_P(0°)和公式(9)计算震源点所在列的网格结点上的旅行时T_s(i_s，j)，j＝1，2，3Λj_s-1，j_s+1，ΛN，N为网格模型的行数。 

6)计算震源点所在列右边网格节点上的旅行时； 

步骤6)所述的计算震源点所在列右边网格节点上的旅行时是将步骤5)计算出的震源点所在列节点上的旅行时间T_s(i_s，j)作为初始时间，计算震源 点所在列右边网格节点上的旅行时，计算顺序按列进行，逐列向右递推，直到计算区域的右边界。 

在每一列的计算过程中，先由上到下，如图3中的模式1所示，设点A纵坐标为z₁，B纵坐标为z₂，射线与AB交点的纵坐标为z₀，利用(10)式，则可以计算出待求点O上的旅行时T₁： 

$T_1=t_0+\frac{\sqrt{{(z_2-z_0)}^2+\mathrm{\Δ}{x}^2}}{V_P\left(\mathrm{\φ}\right)}---\left(10\right)$

其中， 

变量z₀可以利用二分法求解非线性方程(11)式得到： 

f(z₀)＝f₁(z₀)+f₂(z₀)+f₃(z₀)    (11) 

其中， 

$f_1\left(z_0\right)=\frac{{\mathrm{Wz}}_0}{\sqrt{|(z_0^2-z_1^2)t_2^2+(z_2^2-z_0^2)t_1^2|}}$

$f_2\left(z_0\right)=-\frac{(z_2-z_0)}{V_P\left(\mathrm{\φ}\right)\sqrt{{(z_2-z_0)}^2+\mathrm{\Δ}{x}^2}}$

$f_3\left(z_0\right)=\frac{V_P\left(\mathrm{\φ}\right)\·\mathrm{\Δx}\·[a_2\cos \left(\mathrm{\φ}\right)\sin \left(\mathrm{\φ}\right)-2a_3{\cos }^3\left(\mathrm{\φ}\right)\sin \left(\mathrm{\φ}\right)]}{\sqrt{\mathrm{\Δ}{x}^2+{(z_2-z_0)}^2}}$

$W=\frac{t_2^2-t_1^2}{z_2^2-z_1^2}$

φ为射线方向， 

V_P(φ)是射线方向为φ时的射线速度， 

a₁、a₂、a₃为由公式(4)求出的参数； 

如果相邻网格的速度大于当前计算网格的速度，计算沿界面传播的滑行波到达O点的旅行时T₂，取T₁和T₂中的最小值作为O点的旅行时T_s(i，j)。按照图3模式1向下逐行递推，直到模型底边界。再由下向上逐行进行计算，如图3模式2所示，将C点和B点的旅行时间分别当作t₁和t₂，向上逐 步递推，直到顶边界，得到相应网格节点上的旅行时间T_s(i，j)。 

7)计算震源点所在列左边网格节点上的旅行时； 

步骤7)同步骤6)，按照图3模式6由上到下计算，直到模型底边界，再按照图3模式5由下到上计算，直到顶边界，并按列逐步向左，计算震源点所在列左边网格节点上的旅行时间T_s(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)相同，逐步计算各个接收点到网格模型每个网格节点的旅行时间T_R(i，j)。 

12)形成接收点走时文件； 

步骤12)是将步骤11)计算得到的网格模型每个网格节点上的旅行时写到接收点直时文件中，顺序为：接收点号、网格节点横坐标、网格节点纵坐标、旅行时间，图4为计算得到的炮点和接收点的时间场。 

13)波场分离，将井间地震地震记录分离成不同类型的波场； 

步骤13)是利用频率波数滤波或中值滤波的方法，对井间地震记录进行波场分离，分离后的波场包括上行反射纵波，上行转换横波，下行反射纵波、下行转换横波，本发明利用到的波场为波场分离后的上行反射纵波和下行反射纵波，如图5所示。 

14)利用积分法完成各向异性介质成像。 

步骤14)所述的利用积分法成像是：对于任一道井间地震记录，对于任一网格节点r(i，j)，从炮点走时文件对应坐标位置读出该节点的炮点旅行时间T_s(i，j)，从接收点走时文件对应位坐标置读出该节点的接收点旅行时间T_R(i，j)，利用考虑近了场项的积分公式(12)进行积分法成像： 

$U(r,T_s)=\frac{1}{2\mathrm{\πC}}\∫\mathrm{dz}\cos \left(\mathrm{\θ}\right){\∫}_0^{\∞}\mathrm{dT}\·F\left(T\right)\·U(r_0,T+T_s+T_R)---\left(12\right)$

其中， 

T＝t-T_g，U(r₀，t)为在接收点r₀上记录的波场，θ为出射线与接收排列法向矢量的夹角。上行反射纵波和下行反射纵波是分别单独成像后再进行合并，对所有的地震道对进行成像之后，形成井间地震成像剖面，实现利用井间地震反射纵数据对各向异性介质成像的目的，图6为利用本发明方法得到的井间地震成像剖面。 

Claims (7)

1.一种利用井间地震反射波数据对各向异性介质成像的方法，特点是通过以下具体步骤：

1)采集井间地震数据；

2)利用采集的数据反演两口井之间的地质模型，包括：

利用零偏VSP地震记录的纵波初至时间T_P反演模型纵波的垂直速度V_P(0°)，

利用零偏VSP地震记录的横波初至时间T_s反演模型横波的垂直速度V_S(0°)，

利用密度测井数据获得模型密度ρ，

利用井间地震纵波初至时间数据获得纵波的水平速度V_P(90°)和P波45度速度V_P(45°)，

然后以下公式(1)和(2)获得各向异性系数ε和δ建立的层状地质模型：

3)将层状地质模型进行网格剖分；

4)用下式计算每个网格上随射线传播方向ψ变化的射线速度V(ψ)：

V^-2(ψ)≈a₁+a₂cos²(ψ)-a₃cos⁴(ψ)    (3)

其中：

α₀＝V_P(0°)

ε和δ为各向异性系数

V_g[ψ(45°)]和ψ(45°)由公式(5)至(8)求出

5)计算震源点所在列网格节点上的旅行时间；

6)计算震源点所在列右边网格节点上的旅行时；

7)计算震源点所在列左边网格节点上的旅行时；

8)由上到下按行计算炮点到每个网格结点上的旅行时；

9)由下到上按行计算炮点每个网格结点上的旅行时；

10)形成炮点走时文件；

11)计算接收点到每个网格节点的旅行时；

12)形成接收点走时文件；

13)波场分离，将井间地震地震记录分离成不同类型的波场；

14)利用积分法完成各向异性介质成像。

2.根据权利要求1所述的方法，特点是步骤1)所述的采集数据包括采集两口井的零偏VSP地震记录，多分量井间地震记录和密度测井数据。

3.根据权利要求1所述的方法，特点是步骤3)的网格剖分是每个网格包含的模型参数有纵波的垂直速度V_P(0°)，横波的垂直速度V_S(0°)，各向异性系数ε和δ，介质密度ρ。

4.根据权利要求1所述的方法，特点是步骤5)首先确定震源点所在的网格结点(i_s，j_s)，再利用纵波垂直速度V_P(0°)和公式(9)计算震源点所在列的网格结点上的旅行时T_s(i_s，j)，j＝1，2，3Λj_s-1，j_s+1，ΛN，N为网格模型的行数。

5.根据权利要求1所述的方法，特点是步骤6)所述的计算震源点所在列右边网格节点上的旅行时是将步骤5)计算出的震源点所在列节点上的旅行时间T_s(i_s，j)作为初始时间，计算震源点所在列右边网格节点上的旅行时，计算顺序按列进行，逐列向右递推，直到计算区域的右边界。

6.根据权利要求1所述的方法，特点是步骤10)所述的形成炮点走时文件是将步骤6)至步骤9)计算得到的网格模型每个网格节点上的旅行时间按照炮号、网格节点横坐标、网格节点纵坐标、旅行时间排列记录到炮点走时文件中。

7.根据权利要求1所述的方法，特点是步骤13)所述的波场分离是利用频率波数滤波或中值滤波的方法，对井间地震记录进行波场分离，分离后的波场包括上行反射纵波，上行转换横波，下行反射纵波、下行转换横波，本发明利用到的波场为上行反射纵波和下行反射纵波。

Images