CN102053269A - 一种对地震资料中速度分析方法 - Google Patents

Abstract

本发明属于地震勘探和开发领域，尤其涉及一种针对井间地震以及近地表的速度分析方法。本发明针对复杂地质条件地区，近地表横向速度变化大等问题，利用地面地震资料的初至走时，通过基于最短路径的SIRT层析成像技术，建立近地表速度模型，为后续的地震资料处理提供速度资料，从而使得开发的层析成像具有实用性。

Description

一种对地震资料中速度分析方法

技术领域

本发明专利属于地震勘探和开发领域，尤其涉及一种针对井间地震以及近地表的速度分析方法。

背景技术

在地震勘探领域里，对于复杂地质条件地区，如黄土塬等，由于其近地表速度横向变化大，通常无法通过微测井或小折射来获得，因此，如何直接从地震资料中求取复杂的近地表速度结构是近年来的研究热点，层析成像技术是其中之一。

层析成像技术依据其所利用地震资料中的不同参数，如走时或振幅等，通常分为基于射线追踪的走时层析成像和基于波动方程的衍射层析成像两大类。而衍射层析成像由于诸多困难，发展缓慢，目前尚处在理论研究阶段。

走时层析成像又可细分为初至走时层析成像和反射走时层析成像。反射层析成像主要应用于地下速度和反射层深度的反演，以及叠前或叠后偏移的速度分析之中。前者由于速度和深度之间的耦合关系，以及反射波到达时间及其层位难于拾取等，制约了反射波层析成像的应用范围。后者则是利用经过叠前或叠后CRP道集中同相轴未被拉平的剩余时差，经过层析成像来修正用于偏移的速度模型。这种构建速度模型的方法，目前正广泛应用于叠前深度或时间偏移中。

初至走时层析成像发展快速，方法较为简单直观，稳定性较好，且有很多的实际应用，如应用于井间地震以及近地表的速度分析之中等，为其它地震资料处理或提取其它物性参数，提供精确的速度模型。

发明内容

本发明为了解决现有技术中对地震资料速度模型无法切实完成资料准确性要求，而研发了一种对地震资料中速度分析方法。本发明针对复杂地质条件地区，近地表横向速度变化大等问题，利用地面地震资料的初至走时，通过基于最短路径的SIRT层析成像技术，建立近地表速度模型，为后续的地震资料处理提供速度资料，从而使得开发的层析成像具有实用性。

本发明为了解决上述目的，所采用的技术方案为：

本发明方法首先根据建立的观测系统，提取地震资料，并拾取初至走时；后将最短路径树射线追踪技术与SIRT平均化反演算法相结合，形成基于最短路径的SIRT层析成像技术，建立地震资料速度模型，为后续的地震资料处理提供正确的速度资料。

具体所述方法包括如下步骤，

①.建立地面地震观测系统；

②.地震资料预处理步骤

主要包括互相关、数据选排、道编辑等。

③.拾取初至走时步骤；

④.建立速度模型；

⑤.对速度模型进行最短路径树射线追踪；

i.对模型进行空间离散：用矩形来剖分模型，每个矩形内的速度为常数。设X方向的网格数为XN，Z方向的网格数为ZN，则网格总数为XN*ZN；

ii.对矩形进行角度离散：设每个矩形的边界上有M(M≥2)个节点，总节点数为N＝XN*(ZN+1)*M+ZN*(XN+1)*M；

iii.对各节点进行连接：相邻两节点可以相连仅当它们不在同一边界上，则每一个节点的最多连接数为6*M；

iv.计算全局走时及其射线路径：计算从震源s出发到模型中所有节点i的最小走时，即路径树；然后根据费马原理，即最小走时值，寻找射线路径；从震源s到模型中任一个节点i的最小走时应满足Bellman方程：

$t_i=\underset{j\≠i}{\min }(t_j+d_{\mathrm{ij}}/v),i,j\∈N---\left(1\right)$

式中，t_i表示从震源s到节点i的最小走时，t_j表示相邻节点j的最小走时，d_ij表示i、j节点间的距离或长度，

表示i、j节点间的速度平均值；

式(1)应满足的初始条件：t_s＝0，即震源所在的节点处走时为零；

⑥.根据残差进行SIRT层析成像；

SIRT层析成像公式如下：

${{\mathrm{\δS}}_j}^{\left(k\right)}=\frac{1}{M_j}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left\{\right[t_i-\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{S_j}^{\left(k\right)}]a_{\mathrm{ij}}/\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{a_{\mathrm{ij}}}^2\}---\left(2\right)$

式中，M_j表示穿过第j个网格的射线数，n表示射线总数，m表示模型离散化后的网格数；

⑦.对求取的慢度修正量进行平滑处理；

对上述每次求取的慢度修正量δS_j ^(k)进行光滑，其平滑公式如下：

${\mathrm{\δS}}_{i,j}^{\left(k\right)}=({\mathrm{\δS}}_{i-1,j-1}^{\left(k\right)}+{\mathrm{\δS}}_{i-1,j}^{\left(k\right)}+{\mathrm{\δS}}_{i-1,j+1}^{\left(k\right)}+{\mathrm{\δS}}_{i,j-1}^{\left(k\right)}+{\mathrm{\δS}}_{i,j}^{\left(k\right)}+{\mathrm{\δS}}_{i,j+1}^{\left(k\right)}+{\mathrm{\δS}}_{i+1,j-1}^{\left(k\right)}+{\mathrm{\δS}}_{i+1,j}^{\left(k\right)}+{\mathrm{\δS}}_{i+1,j+1}^{\left(k\right)})/9$

式中，i和j分别表示模型网格划分的第i行和第j列；

⑧.输出成像结果。

在具体的应用中，为了加快计算速度，并突出短射线的作用，优选下述平均化算法；

${{\mathrm{\δS}}_j}^{\left(k\right)}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{{\mathrm{\δt}}_i}^{\left(k\right)}}{N_i^4{L}_i}/\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{N_i}^{-4}---\left(3\right)$

式中， ${{\mathrm{\δt}}_i}^{\left(k\right)}=t_i-\underset{j=1}{\overset{m}{\mathrm{\Σ}}}\frac{L_i}{N_i}\·{S_j}^{\left(k\right)}$

$L_i=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}\left|a_{\mathrm{ij}}\right|$

分别表示走时残差δt_i ^(k)和第i条射线的总长度L_i，N_i表示第i条射线穿经模型的网格数。

本发明优点在于：(1)所追踪的走时为全局最小，且射线追踪的精度较高。(2)反演算法可靠稳定、收敛性好，效率高。(3)不仅加快了计算速度，更重要的是突出了短射线的作用。基于最短路径的SIRT初至走时层析成像，就是针对复杂地区近地表速度横向变化大等问题而提出和开发的。经过大量的数学模型试验，其结果表明：该层析成像算法高效可靠，稳定性好。实际资料处理结果表明：该方法能够用于地面地震资料的层析成像，以建立近地表速度模型，在精度和分辨率两方面都得到极大的改善，为后续的地震资料处理提供较为精确的速度资料。

附图说明

图1为本发明方法的流程框图；

图2为采用Marmousi速度模型及其基于最短路径的初至走时层析成像结果图；

图3为本发明的实施例：实际资料反演所用的初始速度模型及其层析成像结果；

图4为本发明的实施例：利用反演得到的近地表速度进行的炮点和检波点静校正图；

图5为利用现有技术对同样的实际资料进行的炮点和检波点静校正图上述各幅附图将结合具体实施方式加以说明

具体实施方式

根据发明内容所述的基于最短路径的初至走时层析成像方法原理，其技术实现的流程框图如图1所示。同时依次进行了数学模型和实际资料的基于最短路径的初至走时层析成像处理。

图2是数学模型的层析成像结果

利用Marmousi模型的浅部来进行数值试验，模型范围是4550m×2800m，如图2(a)所示为截取的Marmousi速度模型。观测系统：炮点检波点位于地面上，检波点距12.5米，炮点距12.5米。总炮数364个，检波器364个。首先利用二维声波有限差分模拟来生成其合成波场，模拟参数为：子波主频35Hz，采样率2ms，采样长度4000ms。

其次在获得的合成记录上进行初至拾取。最后，利用拾取的走时进行层析成像，其层析反演结果如图2(b)所示为基于最短路径的初至走时层析成像结果，与Marmousi速度模型比较，两者基本一致，但对于模型右侧中部的低速层，反演结果较差。

图3是本发明对实际资料的层析成像结果

利用某油田的一条高密度二维线资料来进行试处理，剖面长25km，共拾取走时1208909个，剖分网格为25米*25米，共1001*101个网格，反演所用的初始速度模型如图3(a)，速度值在深度上是渐变的，范围为1200-4200m/s。第三次层析反演结果如图3(b)所示，基于最短路径的初至走时层析成像结果。图3中黑线表示炮点和检波点位置，炮点与检波点最大高差可达323米，黑线上部区域没有意义，仅是为了进行层析成像而添加的。

结果表明：横向速度变化大，尤其在1000至24500m之间有一高速体，速度高达约4200m/s。而在0至1000m之间，速度较低，大约可分为四层，其速度分别约为：2000m/s，2500m/s，3200m/s和4200m/s。

同时，以此为基础进行了层析静校正，其结果如图4所示。其中图4a为层析反演计算得到的炮点静校正量曲线，纵坐标为静校正量(ms)。图4b为层析反演计算得到的检波点静校正量曲线，纵坐标为静校正量(ms)。

并与利用现有技术进行的静校正结果，如图5所示。其中，图5a为现有技术计算得到的炮点静校正量曲线，纵坐标为静校正量(ms)；图5b为现有技术计算得到的检波点静校正量曲线，纵坐标为静校正量(ms)。图4和图5进行比较，明显本发明的静校正结果较优。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的结构，因此前面描述的方式只是优选地，而并不具有限制性的意义。

Claims (3)

1.一种对地震资料中速度分析方法，其特征在于，所述方法首先根据建立的观测系统，提取地震资料，并拾取初至走时；然后将最短路径树射线追踪技术与SIRT平均化反演算法相结合，形成基于最短路径的SIRT层析成像技术，建立地震资料速度模型，为后续的地震资料处理提供正确的速度资料。

2.根据权利要求1所述的一种对地震资料中速度分析方法，其特征在于，所述方法包括如下步骤，

①.建立地面地震观测系统；

②.地震资料预处理步骤：包括互相关、数据选排、道编辑等；

③.拾取初至走时步骤；

④.建立初始速度模型；

⑤.对速度模型进行最短路径树射线追踪；

i.对模型进行空间离散：用矩形来剖分模型，每个矩形内的速度为常数。设X方向的网格数为XN，Z方向的网格数为ZN，则网格总数为XN*ZN；

ii.对矩形进行角度离散：设每个矩形的边界上有M(M≥2)个节点，总节点数为N＝XN*(ZN+1)*M+ZN*(XN+1)*M；

iii.对各节点进行连接：相邻两节点可以相连仅当它们不在同一边界上，则每一个节点的最多连接数为6*M；

iv.计算全局走时及其射线路径：计算从震源s出发到模型中所有节点i的最小走时，即路径树；然后根据费马原理，即最小走时值，寻找射线路径；从震源s到模型中任一个节点i的最小走时应满足Bellman方程：

$t_i=\underset{j\≠i}{\min }(t_j+d_{\mathrm{ij}}/v),i,j\∈N---\left(1\right)$

式中，t_i表示从震源s到节点i的最小走时，t_j表示相邻节点j的最小走时，d_ij表示i、j节点间的距离或长度，

表示i、j节点间的速度平均值；

式(1)应满足的初始条件：t_s＝0，即震源所在的节点处走时为零；

⑥.根据残差进行SIRT层析成像；

SIRT层析成像公式如下：

${{\mathrm{\δS}}_j}^{\left(k\right)}=\frac{1}{M_j}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left\{\right[t_i-\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{S_j}^{\left(k\right)}]a_{\mathrm{ij}}/\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{a_{\mathrm{ij}}}^2\}---\left(2\right)$

式中，M_j表示穿过第j个网格的射线数，n表示射线总数，m表示模型离散化后的网格数；

⑦.对求取的慢度修正量进行平滑处理；

对上述每次求取的慢度修正量δS_j ^(k)进行光滑，其平滑公式如下：

${\mathrm{\δS}}_{i,j}^{\left(k\right)}=({\mathrm{\δS}}_{i-1,j-1}^{\left(k\right)}+{\mathrm{\δS}}_{i-1,j}^{\left(k\right)}+{\mathrm{\δS}}_{i-1,j+1}^{\left(k\right)}+{\mathrm{\δS}}_{i,j-1}^{\left(k\right)}+{\mathrm{\δS}}_{i,j}^{\left(k\right)}+{\mathrm{\δS}}_{i,j+1}^{\left(k\right)}+{\mathrm{\δS}}_{i+1,j-1}^{\left(k\right)}+{\mathrm{\δS}}_{i+1,j}^{\left(k\right)}+{\mathrm{\δS}}_{i+1,j+1}^{\left(k\right)})/9$

式中，i和j分别表示模型网格划分的第i行和第j列；

⑧.如果满足停机准则，则输出层析成像结果；停机准则设定为最大迭代次数或走时残差小于设定的阀值。

3.根据权利要求2所述的一种对地震资料中速度分析方法，其特征在于，

在所述的步骤6中，为了加快计算速度，并突出短射线的作用，优选下述平均化算法；

${{\mathrm{\δS}}_j}^{\left(k\right)}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{{\mathrm{\δt}}_i}^{\left(k\right)}}{N_i^4{L}_i}/\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{N_i}^{-4}---\left(3\right)$

式中， ${{\mathrm{\δt}}_i}^{\left(k\right)}=t_i-\underset{j=1}{\overset{m}{\mathrm{\Σ}}}\frac{L_i}{N_i}\·{S_j}^{\left(k\right)}$

$L_i=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}\left|a_{\mathrm{ij}}\right|$

分别表示走时残差δt_i ^(k)和第i条射线的总长度L_i，N_i表示第i条射线穿经模型的网格数。

Images