CN101630014A - 一种利用垂直地震剖面数据对各向异性介质成像的方法 - Google Patents

Abstract

本发明是利用垂直地震剖面数据对各向异性介质成像的方法，先建立初始的深度－速度模型，再利用Walkaway VSP或三维VSP初至时间反演模型的层速度和各向异性系数建立初始速度模型，将三维空间内任意法线方向的平面转换为法线平行于Z轴的平面，加载观测系统信息，确定从炮点到接收点射线所穿过的反射层，求取近似射线参数，确定反射点在反射界面上的坐标并记录相对应的接收点序号，确定反射点到接收点的旅行时间及射线在反射点的出射角，输出成像道集数据。本发明计算简化，运算效率高，提高了VSP成像结果的精度与分辨率，可方便地对Walkaway VSP及三维VSP纵波、横波或转换波数据进行偏移成像。

Description

一种利用垂直地震剖面数据对各向异性介质成像的方法

技术领域

本发明涉及地球物理勘探技术，是地震资料处理方法，具体是一种利用垂直地震剖面数据对各向异性介质成像的方法。

背景技术

偏移成像是垂直地震剖面(VSP)资料常规处理的一个重要步骤，是通过对采集的VSP地震数据进行重排的反演运算，使地震波能量归位到其空间的真实位置，获取地下真实构造图像。自上世纪80年代以来，已经陆续发展了多种VSP成像方法，例如，VSP-CDP转换方法(Wyatt K.D.和Wyatt S.B.，1981；Cassell和Millahn，1984；Dillon and Thomson，1984，Genmeng Chen，2000)；利用波恩反演的F-K偏移法(Clayton和Yedlin，1980)；克希霍夫积分法((Dillon，1985；Wiggins和Levander，1984，Keho，1984)；逆时偏移法(Chang和McMechan，1986，1987；Yu-Taik Chon，2003)等。

VSP-CDP转换被认为是最为有用的VSP成像方法之一，被广泛地应用于各种VSP数据的成像，但是VSP-CDP转换是基于接收到的地震信号来自地下单一反射点这样一个假设，对于地震信号来自于多个地下反射点的情况，该方法不能使数据正确成确。克希霍夫积分法为另一种常用的VSP成像方法，该方法能处理任意炮点和接收点排列观测采集的地震数据，能对复杂的地质模型进行成像，并且利用正演计算，该方法还可以使折射波和散射波归位成像。

地球上的大多数岩石是各向异性的，它们的物理特性随方向而变化。最简单的各向异性是横向各向同性，就是在水平方向上物性分布均匀，在垂直方向上物性变化。然而，现有的传统VSP数据成像方法，无论是VSP-CDP转换还是克希霍夫积分法，都是基于各向同性介质模型的成像方法。对于零偏VSP和非零偏VSP观测，由于射线传播角度分布在很小的范围内，采集到的数据各向异性特征不明显，并且不存在多炮数据叠加成像的问题，所以可以采用各向同性介质模型对零偏VSP和非零偏VSP数据进行成像。但是，对于Walkaway VSP(井下接收仪器固定于接收井某一深度范围，震源以一定的间距向远离(渐近)井口的方向移动的VSP观测方式)及三维VSP观测，不同井源距的激发点的射线传播方向变化很大，采集到的VSP数据各向异性特征明显；并且Walkaway VSP及三维VSP数据成像需要进行多炮数据叠加，而不同激发点上采集的数据的各向异性强度又是不同的，所以利用传统的基于各向同性介质模型的成像方法对Walkaway VSP及三维VSP数据进行成像会给成像结果带来一些问题：一是由各性异性给成像位置带来较大误差；二是各向异性程度不同的不同炮的数据进行不同相叠加造成的分辨率降低及明显的“接缝”。

发明内容

本发明提供一种运算效率高，提高精度与分辨率的利用VSP数据对各向异性介质成像的方法。

本发明通过以下技术方案实现，具体步骤是：

1)采集VSP叠前地震记录；

步骤1)所述的地震记录为x(t)，t＝1，…，N，1，…，N为采样点。

2)波场分离，将VSP叠前地震记录分离成不同类型的波场；

步骤2)所述的波场包括上行反射纵波，上行转换波，下行转换波。

3)用零井源距VSP初至时间求得VSP层速度，建立初始的速度-深度模型，再利用Walkaway VSP或三维VSP的多条初至时间-深度时距曲线反演垂直对称轴各向异性(VTI)模型的层速度和各向异性系数；

4)用每个反射界面的深度、每一层的纵波速度、纵波各向异性系数、横波速度、横波各向异性系数，反射界平面的与Z轴的夹角

反射界面的与X轴的夹角θ建立初始速度模型；根据坐标变换公式将三维空间内任意法线方向的平面转换为法线平行于Z轴的平面；

步骤4)所述的坐标变换公式为：

$\left\{\begin{array}{c}x=l_{1}X+l_{2}Y+l_{3}Z\\ y=m_{1}X+m_{2}Y+m_{3}Z\\ z=n_{1}X+n_{2}Y+n_{3}Z\end{array}\right.$

上式中：

l₃＝sinθsinγ，

m₃＝-sinθcosγ，

n₂＝cosθ。

其中，X，Y，Z分别为坐标转换前平面的三个坐标轴，x，y，z分别为坐标转换后平面的三个坐标轴。θ、γ、

分别为平面的法线与X、Y、Z三个坐标轴的夹角。

5)加载要偏移的地震数据并利用地震数据道头加载观测系统信息；

步骤5)所述的观测系统信息包括炮点X坐标，炮点Y坐标，炮点Z坐标，炮点个数，接收点X坐标，接收点Y坐标，接收点Z坐标及接收点个数。

6)对于每一个炮点、每一个接收点和模型的每一个反射层，根据炮点、接收点和当前反射层在模型内的位置，确定从炮点到接收点射线所穿过的反射层；

7)用射线参数P(θ)的函数 $F\left[P\left(\mathrm{\θ}\right)\right]=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k{A}_k^2{V}_k^{2}P\left(\mathrm{\θ}\right)}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left(\mathrm{\θ}\right)}}-L,$ 利用二分法求取近似射线参数P(θ)，

上式中，k为模型层序号，H_k为第k层的厚度，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，L为炮点到接收井的水平距离；

8)做当前炮点到各个接收点的射线追踪，利用公式 $x=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k{A}_k^2{V}_k^{2}P\left({\mathrm{\θ}}_1\right)}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定反射点在第k个反射界面上的坐标，并记录相应的接收点序号，上式中，k为模型层序号，H_k为第k层的厚度，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，P(θ₁)为第1层的射线参数；

9)根据给定的成像面元大小，把第k个反射界面剖分成一系列等间隔的点，将这些点当作界面上的反射点，计算出每一个反射点的坐标，查找邻近的通过射线追踪计算出的反射点，并记录相对应的接收点序号；

10)做从当前炮点到第k个反射界上等间隔的反射点的射线追踪；

步骤10)所述的射线追踪是利用公式 $t=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定炮点到反射点的旅行时间，利用公式 $\tan {\mathrm{\Φ}}_k=\frac{A_k^2{V}_{k}P\left({\mathrm{\θ}}_1\right)}{\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定射线传播到反射点的入射角，上式中，k为模型层序号，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，P(θ₁)为第1层的射线参数；

11)做从第k个反射界上等间隔的反射点到与之相对应的接收点的射线追踪；

步骤11)所述的射线追踪是根据公式 $t=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定反射点到接收点的旅行时间，根据公式 $\tan {\mathrm{\Φ}}_k=\frac{A_k^2{V}_{k}P\left({\mathrm{\θ}}_1\right)}{\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定射线在反射点的出射角，上式中，k为模型层序号，H_k为第k层的厚度，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，P(θ₁)为第1层的射线参数。

12)根据步骤10)和步骤11)计算得到炮点到反射点的旅行时间与反射点到接收点的旅行时间之和，将加载的地震数据相应接收点上相应时间上的幅值乘以与射线传播方向相关的系数置于相应的成像点上；

13)做从第k个反射界上等间隔的反射点到与井径上各个接收点的射线追踪，采用与步骤11)相同的方法确定反射点到接收点的旅行时间及射线在反射点的出射角，如果步骤13)计算得到的射线出射角与步骤11)计算得到的射线出射角之差小于给定的偏移孔径角度，则按照步骤12)抽取地震数据进行成像；

14)根据给定的成像网格大小，将成像范围内的数据选排成道集，为各道成像数据置道头信息，输出成像道集数据。

本发明可简化计算，提高运算效率，消除由各向异性程度不同的各炮数据对偏移成像带来的影响，提高VSP成像结果的精度与分辨率，还可以方便的对Walkaway VSP及三维VSP纵波、横波或转换波数据进行偏移成像。

附图说明

图1(a)为本发明Kirchhoff积分法偏移成像流程图；

图1(b)为本发明Kirchhoff积分法偏移成像流程图；

图2为应用本发明对Walkaway VSP反射P波模型数据进行成像的结果；

图3是应用本发明对Walkaway VSP反射S波模型数据进行成像的结果；

图4是应用本发明对三维VSP反射P波模型数据进行成像的结果(第48线)；

图5是应用本发明对三维VSP反射P波模型数据进行成像的结果(第54线)；

图6是应用本发明对三维VSP反射S波模型数据进行成像的结果(第48线)；

图7是应用本发明对三维VSP反射S波模型数据进行成像的结果(第54线)；

图8是将应用本发明对Walkaway VSP反射P波实际数据进行成像的结果嵌入地面地震剖面的对比图。

具体实施方式

以下结合附图详细说明本发明。

本发明的具体实施方式为：

1)正演模拟Walkaway VSP和三维VSP合成地震记录，采集实际WalkawayVSP地震记录。

2)对合成和实际的VSP地震记录进行波场分离，分离的波场包括上行反射纵波，上行转换波。

3)用零井源距VSP初至时间求得VSP层速度，为实际VSP数据偏移建立初始的速度-深度模型，再利用Walkaway VSP或三维VSP的多条初至时间-深度时距曲线反演VTI(垂直对称轴各向异性)模型的层速度和各向异性系数。

4)加载初始速度模型。初始速度模型包括每个反射界面的深度、每一层的纵波速度、纵波各向异性系数、横波速度、横波各向异性系数，反射界平面的与Z轴的夹角

反射界面的与X轴的夹角θ。根据坐标变换公式：

$\left\{\begin{array}{c}x=l_{1}X+l_{2}Y+l_{3}Z\\ y=m_{1}X+m_{2}Y+m_{3}Z\\ z=n_{1}X+n_{2}Y+n_{3}Z\end{array}\right.$

将三维空间内任意法线方向的平面转换为法线平行于Z轴的平面，上式中：

l₃＝sinθsinγ，

m₃＝-sinθcosγ，

n₂＝cosθ。

其中，X，Y，Z分别为坐标转换前平面的三个坐标轴，x，y，z分别为坐标转换后平面的三个坐标轴。θ、γ、分别为平面的法线与X、Y、Z三个坐标轴的夹角。

5)加载要偏移的地震数据和观测系统。加载的地震数据为波场分离后的上行反射纵波和上行转换波数据，从输入的地震数据道头中获取观测系统信息，包括炮点X坐标，炮点Y坐标，炮点Z坐标，炮点个数，接收点X坐标，接收点Y坐标，接收点Z坐标及接收点个数。

6)对于每一个炮点、每一个接收点和模型的每一个反射层，根据炮点、接收点和当前反射层在模型内的位置，确定从炮点到接收点射线所穿过的反射层。

7)构造关于射线参数P(θ)的函数 $F\left[P\left(\mathrm{\θ}\right)\right]=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k{A}_k^2{V}_k^{2}P\left(\mathrm{\θ}\right)}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left(\mathrm{\θ}\right)}}-L,$ 利用二分法求取近似射线参数P(θ)。上式中，k为模型层序号，H_k为第k层的厚度，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，L为炮点到接收井的水平距离。

8)做当前炮点到各个接收点的射线追踪，利用公式 $x=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k{A}_k^2{V}_k^{2}P\left({\mathrm{\θ}}_1\right)}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ (将公式插入在这里)确定反射点在第k个反射界面上的坐标，并记录相应的接收点序号。上式中，k为模型层序号，H_k为第k层的厚度，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，P(θ₁)为第1层的射线参数。

9)根据给定的成像面元大小，把第k个反射界面剖分成一系列等间隔的点，将这些点当作界面上的反射点，计算出每一个反射点的坐标，查找邻近的通过射线追踪计算出的反射点，并记录相对应的接收点序号。

10)做从当前炮点到第k个反射界上等间隔的反射点的射线追踪，利用公式 $t=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定炮点到反射点的旅行时间，利用公式 $\tan {\mathrm{\Φ}}_k=\frac{A_k^2{V}_{k}P\left({\mathrm{\θ}}_1\right)}{\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ (将公式插入在这里)确定射线传播到反射点的入射角。

上式中，k为模型层序号，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，P(θ₁)为第1层的射线参数

11)做从第k个反射界上等间隔的反射点到与之相对应的接收点的射线追踪，根据公式 $t=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ (将公式插入在这里)确定反射点到接收点的旅行时间，根据公式 $\tan {\mathrm{\Φ}}_k=\frac{A_k^2{V}_{k}P\left({\mathrm{\θ}}_1\right)}{\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定射线在反射点的出射角。上式中，k为模型层序号，H_k为第k层的厚度，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，P(θ₁)为第1层的射线参数

12)根据步骤10)和步骤11)计算得到炮点到反射点的旅行时间与反射点到接收点的旅行时间之和，将加载的地震数据相应接收点上相应时间上的幅值乘以与射线传播方向相关的系数置于相应的成像点上。

13)做从第k个反射界上等间隔的反射点到与井径上各个接收点的射线追踪，采用与步骤11)相同的方法确定反射点到接收点的旅行时间及射线在反射点的出射角，如果步骤13)计算得到的射线出射角与步骤11)计算得到的射线出射角之差小于给定的偏移孔径角度，则按照步骤12)抽取地震数据进行成像。

14)根据给定的成像网格大小，将成像范围内的数据选排成道集，为各道成像数据置道头信息，输出成像道集数据。如附图所示，图2为应用本发明对Walkaway VSP反射P波模型数据进行成像的结果，图3是应用本发明对Walkaway VSP反射S波模型数据进行成像的结果，图4是应用本发明对三维VSP反射P波模型数据进行成像的结果(第48线)，图5是应用本发明对三维VSP反射P波模型数据进行成像的结果(第54线)，图6是应用本发明对三维VSP反射S波模型数据进行成像的结果(第48线)，图7是应用本发明对三维VSP反射S波模型数据进行成像的结果(第54线)，图8是将应用本发明对Walkaway VSP反射P波实际数据进行成像的结果嵌入地面。

本发明利用基于椭圆各向异性介质模型的Kirchhoff积分法对WalkawayVSP及三维VSP数据进行偏移成像，可以消除由各向异性程度不同的各炮数据对偏移成像带来的影响，提高VSP成像结果的精度与分辨率；并且该方法通过坐标系的转换，将三维空间内任意法线方向的平面转换为法线平行于Z轴的平面，具有极高的计算效率，完全能够满足实际生产中Walkaway VSP及三维VSP数据偏移成像对计算效率的要求。

Claims (7)

1、一种利用垂直地震剖面数据对各向异性介质成像的方法，其特征在于通过以下具体步骤实现：

1)采集VSP叠前地震记录；

2)波场分离，将VSP叠前地震记录分离成不同类型的波场；

3)用零井源距VSP初至时间求得VSP层速度，建立初始的速度-深度模型，再利用Walkaway VSP或三维VSP的多条初至时间-深度时距曲线反演垂直对称轴各向异性(VTI)模型的层速度和各向异性系数；

4)用每个反射界面的深度、每一层的纵波速度、纵波各向异性系数、横波速度、横波各向异性系数，反射界平面的与Z轴的夹角

反射界面的与X轴的夹角θ建立初始速度模型；根据坐标变换公式将三维空间内任意法线方向的平面转换为法线平行于Z轴的平面；

5)加载要偏移的地震数据并利用地震数据道头加载观测系统信息；

6)对于每一个炮点、每一个接收点和模型的每一个反射层，根据炮点、接收点和当前反射层在模型内的位置，确定从炮点到接收点射线所穿过的反射层；

7)用射线参数P(θ)的函数 $F\left[P\left(\mathrm{\θ}\right)\right]=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k{A}_k^2{V}_k^{2}P\left(\mathrm{\θ}\right)}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left(\mathrm{\θ}\right)}}-L,$ 利用二分法求取近似射线参数P(θ)，

上式中，k为模型层序号，H_k为第k层的厚度，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，L为炮点到接收井的水平距离；

8)做当前炮点到各个接收点的射线追踪，利用公式 $x=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k{A}_k^2{V}_k^{2}P\left({\mathrm{\θ}}_1\right)}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定反射点在第k个反射界面上的坐标，并记录相应的接收点序号，上式中，k为模型层序号，H_k为第k层的厚度，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，P(θ₁)为第1层的射线参数；

9)根据给定的成像面元大小，把第k个反射界面剖分成一系列等间隔的点，将这些点当作界面上的反射点，计算出每一个反射点的坐标，查找邻近的通过射线追踪计算出的反射点，并记录相对应的接收点序号；

10)做从当前炮点到第k个反射界上等间隔的反射点的射线追踪；

11)做从第k个反射界上等间隔的反射点到与之相对应的接收点的射线追踪；

12)根据步骤10)和步骤11)计算得到炮点到反射点的旅行时间与反射点到接收点的旅行时间之和，将加载的地震数据相应接收点上相应时间上的幅值乘以与射线传播方向相关的系数置于相应的成像点上；

13)做从第k个反射界上等间隔的反射点到与井径上各个接收点的射线追踪，采用与步骤11)相同的方法确定反射点到接收点的旅行时间及射线在反射点的出射角，如果步骤13)计算得到的射线出射角与步骤11)计算得到的射线出射角之差小于给定的偏移孔径角度，则按照步骤12)抽取地震数据进行成像；

14)根据给定的成像网格大小，将成像范围内的数据选排成道集，为各道成像数据置道头信息，输出成像道集数据。

2、根据权利要求1所述的利用垂直地震剖面数据对各向异性介质成像的方法，其特征在于步骤1)所述的地震记录为x(t)，t＝1，…，N，1，…，N为采样点。

3、根据权利要求1所述的利用垂直地震剖面数据对各向异性介质成像的方法，其特征在于步骤2)所述的波场包括上行反射纵波，上行转换波，下行转换波。

4、根据权利要求1所述的利用垂直地震剖面数据对各向异性介质成像的方法，其特征在于步骤4)所述的坐标变换公式为：

$\left\{\begin{array}{c}x=l_{1}X+l_{2}Y+l_{3}Z\\ y=m_{1}X+m_{2}Y+m_{3}Z\\ z=n_{1}X+n_{2}Y+n_{3}Z\end{array}\right.$

上式中：

l₃＝sinθsinγ，

m₃＝-sinθcosγ，

n₂＝cosθ。

其中，X，Y，Z分别为坐标转换前平面的三个坐标轴，x，y，z分别为坐标转换后平面的三个坐标轴。θ、γ、

分别为平面的法线与X、Y、Z三个坐标轴的夹角。

5、根据权利要求1所述的利用垂直地震剖面数据对各向异性介质成像的方法，其特征在于步骤5)所述的观测系统信息包括炮点X坐标，炮点Y坐标，炮点Z坐标，炮点个数，接收点X坐标，接收点Y坐标，接收点Z坐标及接收点个数。

6、根据权利要求1所述的利用垂直地震剖面数据对各向异性介质成像的方法，其特征在于步骤10)所述的射线追踪是利用公式 $t=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定炮点到反射点的旅行时间，利用公式 $\tan {\mathrm{\Φ}}_k=\frac{A_k^2{V}_{k}P\left({\mathrm{\θ}}_1\right)}{\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定射线传播到反射点的入射角，上式中，k为模型层序号，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，P(θ₁)为第1层的射线参数。

7、根据权利要求1所述的利用垂直地震剖面数据对各向异性介质成像的方法，其特征在于步骤11)所述的射线追踪是根据公式 $t=\underset{k=k_o}{\overset{k_r}{\mathrm{\Σ}}}\frac{H_k}{V_k\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定反射点到接收点的旅行时间，根据公式tan ${\mathrm{\Φ}}_k=\frac{A_k^2{V}_{k}P\left({\mathrm{\θ}}_1\right)}{\sqrt{1-A_k^2{V}_k^2{P}^2\left({\mathrm{\θ}}_1\right)}}$ 确定射线在反射点的出射角，上式中，k为模型层序号，H_k为第k层的厚度，V_k为第k层的垂向层速度，A_k为第k层的速度各向异性系数，P(θ₁)为第1层的射线参数。

Images