一种地震叠前时间偏移中的射线追踪方法

摘要

本发明提供了一种地震叠前时间偏移中的射线追踪方法，包括：步骤1，获取M个反射界面的模型参数以及N个炮检距参数；步骤2，计算来自K个界面第j个炮检距的反射波传播路径初值参数a

说明书

技术领域

本发明涉及地震数据处理领域，尤其涉及一种地震叠前时间偏移中的射线追踪方 法。

背景技术

地震叠前偏移方法目前已经成为地震资料成像的主要方法，它可以分为叠前时间 偏移和叠前深度偏移方法。尽管叠前深度偏移能够对地下复杂构造进行准确成像，但 是由于其偏移时采用的是层速度模型，其成像精度对层速度模型精度要求很高。在一 些复杂地区，由于难以获得准确的层速度模型，因而从一定程度上限制了该方法的应 用。叠前时间偏移方法通常利用均方根速度模型进行偏移成像，速度模型相对容易获 得，成像效果对速度模型的依赖程度没有深度偏移方法那样严重，但其成像效果相对 于叠后偏移成像而言，能够得到显著改善，因而叠前时间偏移目前在实际地震资料成 像处理中得到了较为广泛的应用。叠前时间偏移主要采用克希霍夫积分法，通常基于 均方根速度模型，即将成像点上方的介质等效为均匀介质，速度为均方根速度，这类 方法是基于直射线的偏移方法。为了提高叠前时间偏移成像精度，发展了基于弯曲射 线路径的叠前时间偏移方法，将成像点上方的介质用多层水平层状模型来描述，此时 地震波的传播路径变成了折线段。为了获得地震波在水平层状介质中的传播时间，需 要通过射线追踪方法追踪出地震波传播路径，进而可以计算出旅行时间，用于叠前时 间偏移。在基于弯曲射线路径的叠前时间偏移中，由于要对地下各成像点、各炮检距 进行射线追踪，因而计算量很大，需要研究和利用快速射线追踪方法。

随着地震解释和地震反演向叠前发展，叠前合成记录制作变得较为重要。基于水 平层状介质的合成记录可以应用于纵波和转换波AVO分析、AVO反演和叠前反演 等，叠前合成记录制作的主要方法之一亦是射线追踪方法。

水平层状介质的射线追踪通常采用试射法。对于给定的炮检距（炮点和检波点）， 该方法给出反射波在第一层的入射角，根据斯奈尔定律计算出传播路径和在地面出射 点位置；当出射点位于炮检连线的线段外侧时，表明入射角偏大，需要减小入射角； 当出射点位于炮检连线的线段上时，表明入射角偏小，需要增大入射角；当出射点和 检波点位置误差达到给定误差时，即可获得反射波路径和时间。该方法在进行试射时， 需要确定试射入射角的范围，并利用二分法来更新入射角。实际上，试射入射角的范 围越大，所需的计算量越大。由于要针对不同界面、不同炮检距进行射线追踪，通常 按照界面深度由浅到深、炮检距由小到大进行射线追踪，因此，在实际追踪某界面、 某炮检距的反射波路径时，试射的最小入射角可以设为同一界面、相邻前一个炮检距 的入射角，试射的最大入射角可以设为同一炮检距、相邻前一个界面的入射角，由此 设计的入射角试射范围很小，计算量可以显著减小。但是，通过计算研究发现，要获 得一定精度的反射波传播路径和传播时间，仍然需要多次迭代。

总之，叠前时间偏移和合成记录制作需要进行射线追踪，常用的射线追踪方法需 要多次迭代，计算量较大。

发明内容

针对上述问题，本发明的目的是提供一种针对叠前时间偏移和合成记录制作的， 且具有计算速度快、精度高、始终收敛的射线追踪方法。

为了达到上述目的，本发明实施例提供一种地震叠前时间偏移中的射线追踪方 法，包括：步骤1，获取M个反射界面的模型参数以及N个炮检距参数；步骤2，根 据K个界面的模型参数和炮检距参数，计算来自K个界面第j个炮检距的反射波传 播路径初值参数a_i⁽¹⁾及传播时间i=1,2,...,M，k=1,2,...,M，j=1,2,...,N；步骤 3，根据设定的迭代系数，利用迭代法计算来自K个界面第j个炮检距的新的反射波 传播路径初值参数a_i⁽²⁾，并计算来自K个界面第j个炮检距的新的反射波传播时间 步骤4，如果ε为设定的误差，则即为来自K个界面第j个炮 检距的反射波传播时间，否则，返回所述步骤3继续迭代；步骤5，如果j<N，则转 向步骤2-4继续计算，否则，来自K个界面反射波传播时间计算完成；步骤6，如果 K<M，则转向步骤2-5继续计算，否则，最终的反射波传播时间计算完成。

进一步地，在一实施例中，在所述步骤1中，获取的M个反射界面的模型参数 包括：各层层速度v_k以及各层厚度h_k，其中，k=1,2,...,M。

进一步地，在一实施例中，在所述步骤2中，首先依据下述公式计算第k=1个 界面的各个炮检距的反射波时间：其中，x_j为第 k=1个界面上的N个由小到大排列的炮检距参数，v_k为第k=1个界面上的层速度， h_k为第k=1个界面的层厚度。

进一步地，在一实施例中，在所述步骤2中，根据K个界面的模型参数和炮检 距参数，计算来自K个界面第j个炮检距的反射波传播路径初值参数a_i⁽¹⁾及传播时间 包括：如果j=1，则${a_i}^{\left(1\right)}=\frac{v_i{h}_i}{\underset{l=1}{\overset{M}{\Sigma }}{v}_l{h}_l}{x}_1,$否则，${a_i}^{\left(1\right)}={a_i}^{\left(1\right)}+\frac{v_i{h}_i}{\underset{l=1}{\overset{M}{\Sigma }}{v}_l{h}_l}(x_j-x_{j-1}),$i=1,2,...,M；以及，$t_{k,j}^{\left(1\right)}=\underset{i=1}{\overset{M}{\Sigma }}{({\left[a_i^{\left(1\right)}\right]}^2+4h_i^2)}^{1/2}/v_i.$

进一步地，在一实施例中，在所述步骤3中，如果步骤2中的迭代次数小于最大 迭代次数iter，则根据设定的迭代系数，利用迭代法计算来自K个界面第j个炮检距 的新的反射波传播路径初值参数a_i⁽²⁾，包括：

当v_i≥v_i+1时，i=1,...,M-1，迭代系数$c_i={({r_i}^2+{r_i}^2\frac{{\left[a_i^{\left(1\right)}\right]}^2}{h_i^2}-\frac{{\left[a_i^{\left(1\right)}\right]}^2}{h_{i+1}^2})}^{1/2},$${r_i}^2=\frac{v_i^2{h}_i^2}{v_{i+1}^2{h}_{i+1}^2};$

当v_i<v_i+1时，i=1,...,M-1，迭代系数$c_i={({r_i}^2+{r_i}^2\frac{{\left[a_{i+1}^{\left(1\right)}\right]}^2}{h_{i+1}^2}-\frac{{\left[a_{i+1}^{\left(1\right)}\right]}^2}{h_i^2})}^{1/2},$${r_i}^2=\frac{v_{i+1}^2{h}_{i+1}^2}{v_i^2{h}_i^2};$

迭代系数d_i=a_i⁽¹⁾+a_i+1⁽¹⁾；

则根据设定的迭代系数c_i和d_i，计算来自K个界面第j个炮检距的新的反射波传 播路径初值参数：

${a_i}^{\left(2\right)}=\frac{c_i}{1+c_i}{d}_i,$${a_{i+1}}^{\left(2\right)}=\frac{1}{1+c_i}{d}_i,i=1,...,M-1.$

进一步地，在一实施例中，在所述步骤3中，计算来自K个界面第j个炮检距的 新的反射波传播时间包括：

根据公式计算来自K个界面第j个炮检距的新的反射 波传播时间，其中，a_i⁽²⁾为新的反射波传播路径初值参数。

进一步地，在一实施例中，在所述步骤4中，如果迭代次数大于最大迭代次数iter， 则结束迭代，否则转向步骤3继续迭代。

本发明实施例的地震叠前时间偏移中的射线追踪方法，基于高精度的初值公式和 快速逼近真解的迭代公式，具有计算速度快、精度高和始终收敛等优点，可以用于基 于弯曲射线的纵波和转换波叠前时间偏移、纵波和转换波叠前合成记录制作。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现 有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅 是本发明的一些实施例，对于本领域技术人员来讲，在不付出创造性劳动性的前提下， 还可以根据这些附图获得其他的附图。

图1为本发明实施例的地震叠前时间偏移中的射线追踪方法的流程示意图；

图2为本发明的地震叠前时间偏移中的射线追踪方法的另一个实施例的流程示 意图；

图3为该具体实施例中的模型示意图；

图4为利用本发明的迭代后的初值公式计算的反射波时间与精确值的对比图；

图5为在相同精度条件下，本发明的创新方法与传统方法（试射法）平均迭代次 数对比图；

图6为本发明的方法与传统方法（试射法）的CPU时间对比图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完 整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。 基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的 所有其他实施例，都属于本发明保护的范围。

图1为本发明实施例的地震叠前时间偏移中的射线追踪方法的流程示意图。如图 所示，所述方法包括以下步骤：

步骤S101，获取M个反射界面的模型参数以及N个炮检距参数；

步骤S102，根据K个界面的模型参数和炮检距参数，计算来自K个界面第j个 炮检距的反射波传播路径初值参数a_i⁽¹⁾及传播时间i=1,2,...,M，k=1,2,...,M， j=1,2,...,N；

步骤S103，根据设定的迭代系数，利用迭代法计算来自K个界面第j个炮检距 的新的反射波传播路径初值参数a_i⁽²⁾，并计算来自K个界面第j个炮检距的新的反射 波传播时间

步骤S104，如果ε为设定的误差，则即为来自K个界面第j 个炮检距的反射波传播时间，否则，返回所述步骤S103继续迭代；

步骤S105，如果j<N，则转向步骤S102-步骤S104继续计算，否则，来自K个 界面反射波传播时间计算完成；

步骤S106，如果K<M，则转向步骤S102-步骤S105继续计算，否则，最终的反 射波传播时间计算完成。

在本实施例中，由于初值公式较为精确、迭代公式能够快速逼近真解，因而具有 计算速度快、精度高以及始终收敛的优点。

在实施的过程中，在本实施例的步骤S101中，获取的M个反射界面的模型参数 包括：各层层速度v_k以及各层厚度h_k，其中，k=1,2,...,M。这些参数可以从地震资 料中直接获取得到。同时，可以获取到N个由小到大排列的炮检距参数x_j， j=1,2,...,N。

在本实施例的步骤S102中，首先依据下述公式计算第k=1个界面的各个炮检距 的反射波时间：$t_{1,j}^{\left(2\right)}={(x_j^2+4h_1^2)}^{1/2}/v_1,j=\mathrm{1,2},...,N;$

其中，x_j为第k=1个界面上的N个由小到大排列的炮检距参数，v_k为第k=1个 界面上的层速度，h_k为第k=1个界面的层厚度。

这种初值公式的计算较为精确，提高了计算的精度。

在本实施例的步骤S102中，根据K个界面的模型参数和炮检距参数，计算来自 K个界面第j个炮检距的反射波传播路径初值参数a_i⁽¹⁾及传播时间包括：

如果j=1，则${a_i}^{\left(1\right)}=\frac{v_i{h}_i}{\underset{l=1}{\overset{M}{\Sigma }}{v}_l{h}_l}{x}_1,$否则，${a_i}^{\left(1\right)}={a_i}^{\left(1\right)}+\frac{v_i{h}_i}{\underset{l=1}{\overset{M}{\Sigma }}{v}_l{h}_l}(x_j-x_{j-1}),i=\mathrm{1,2},...,M;$ 以及，计算反射波传播时间$t_{k,j}^{\left(1\right)}=\underset{i=1}{\overset{M}{\Sigma }}{({\left[a_i^{\left(1\right)}\right]}^2+4h_i^2)}^{1/2}/v_i.$

在本实施例的步骤S103中，如果步骤S102中的迭代次数小于最大迭代次数iter， 则根据设定的迭代系数，利用迭代法计算来自K个界面第j个炮检距的新的反射波传 播路径初值参数a_i⁽²⁾，包括：

当v_i≥v_i+1时，i=1,...,M-1，迭代系数$c_i={({r_i}^2+{r_i}^2\frac{{\left[a_i^{\left(1\right)}\right]}^2}{h_i^2}-\frac{{\left[a_i^{\left(1\right)}\right]}^2}{h_{i+1}^2})}^{1/2},$${r_i}^2=\frac{v_i^2{h}_i^2}{v_{i+1}^2{h}_{i+1}^2};$

当v_i<v_i+1时，i=1,...,M-1，迭代系数$c_i={({r_i}^2+{r_i}^2\frac{{\left[a_{i+1}^{\left(1\right)}\right]}^2}{h_{i+1}^2}-\frac{{\left[a_{i+1}^{\left(1\right)}\right]}^2}{h_i^2})}^{1/2},$${r_i}^2=\frac{v_{i+1}^2{h}_{i+1}^2}{v_i^2{h}_i^2};$

迭代系数d_i=a_i⁽¹⁾+a_i+1⁽¹⁾；

则根据设定的迭代系数c_i和d_i，计算来自K个界面第j个炮检距的新的反射波传 播路径初值参数：

${a_i}^{\left(2\right)}=\frac{c_i}{1+c_i}{d}_i,$${a_{i+1}}^{\left(2\right)}=\frac{1}{1+c_i}{d}_i,i=1,...,M-1.$

在本实施例的步骤S103中，计算来自K个界面第j个炮检距的新的反射波传播 时间包括：

根据公式计算来自K个界面第j个炮检距的新的反射 波传播时间，其中，a_i⁽²⁾为新的反射波传播路径初值参数。

在本实施例的步骤S104中，如果迭代次数大于最大迭代次数iter，则结束迭代， 否则转向步骤S103继续迭代。

为了使得本发明的目的、技术方案和优点更加清楚明白，下面结合一个具体的实 施例，对本发明进行进一步的详细说明。在此，本发明的示意性实施方式及其说明用 于解释本发明，但并构成对本发明不当的限定。

在本例中，提出了一种地震叠前时间偏移中的射线追踪方法，如图2所示，具体 包括：

步骤S10，给定有M个反射界面的模型参数和N个炮检距参数：M个反射界面 的模型参数包括：各层层速度v_k以及各层厚度h_k，其中，k=1,2,...,M；N个由小到 大排列的炮检距参数x_j，j=1,2,...,N。

步骤S20，计算来自第k=1个界面各个炮检距的反射波时间：j=1,2,...,N，其中，x_j为第k=1个界面上的N个由小到大排列的炮检距参数，v_k为 第k=1个界面上的层速度，h_k为第k=1个界面的层厚度；

步骤S30，k=k+1，j=1。

步骤S40，计算来自k个界面第j个炮检距的反射波传播路径初值参数 a_i(i=1,2,...,M)及传播时间

如果j=1，则${a_i}^{\left(1\right)}=\frac{v_i{h}_i}{\underset{l=1}{\overset{M}{\Sigma }}{v}_l{h}_l}{x}_1,$否则，${a_i}^{\left(1\right)}={a_i}^{\left(1\right)}+\frac{v_i{h}_i}{\underset{l=1}{\overset{M}{\Sigma }}{v}_l{h}_l}(x_j-x_{j-1}),i=\mathrm{1,2},...,M;$以及，计算反射波传播时间$t_{k,j}^{\left(1\right)}=\underset{i=1}{\overset{M}{\Sigma }}{({\left[a_i^{\left(1\right)}\right]}^2+4h_i^2)}^{1/2}/v_i.$

步骤S50，iter=1；

步骤S60，利用公式计算新的a_i(i=1,2,...,M)；当v_i≥v_i+1时，i=1,...,M-1，迭代 系数$c_i={({r_i}^2+{r_i}^2\frac{{\left[a_i^{\left(1\right)}\right]}^2}{h_i^2}-\frac{{\left[a_i^{\left(1\right)}\right]}^2}{h_{i+1}^2})}^{1/2},$${r_i}^2=\frac{v_i^2{h}_i^2}{v_{i+1}^2{h}_{i+1}^2};$

当v_i<v_i+1时，i=1,...,M-1，迭代系数$c_i={({r_i}^2+{r_i}^2\frac{{\left[a_{i+1}^{\left(1\right)}\right]}^2}{h_{i+1}^2}-\frac{{\left[a_{i+1}^{\left(1\right)}\right]}^2}{h_i^2})}^{1/2},$${r_i}^2=\frac{v_{i+1}^2{h}_{i+1}^2}{v_i^2{h}_i^2};$

迭代系数d_i=a_i⁽¹⁾+a_i+1⁽¹⁾；

则根据设定的迭代系数c_i和d_i，计算来自K个界面第j个炮检距的新的反射波传 播路径初值参数：

${a_i}^{\left(2\right)}=\frac{c_i}{1+c_i}{d}_i,$${a_{i+1}}^{\left(2\right)}=\frac{1}{1+c_i}{d}_i,i=1,...,M-1.$

步骤S70，计算来自k个界面第j个炮检距的反射波传播时间根据公式 计算来自K个界面第j个炮检距的新的反射波传播时间， 其中，a_i⁽²⁾为新的反射波传播路径初值参数。

步骤S80，判断或者iter>MAXiter，如果成立，则转入步骤S90，即为来自k个界面第j个炮检距的反射波传播时间，输出该时间；如果迭代次数超过 了最大迭代次数，则结束迭代；否则转入步骤S100，iter=iter+1后转向 步骤S60，继续迭代；

步骤S110，判断是否j<N，如果成立则转入步骤S120，j=j+1后转向步骤 S70，继续计算；如果不成立，计算来自k个界面反射波传播时间完成，转向步骤S130；

步骤S130，判断是否k<M，如果成立则转入步骤S140，k=k+1后转向步骤 S40，继续计算；如果不成立，计算反射波传播时间完成。

具体实施例：

图3为该具体实施例中的模型示意图。各层速度由浅至深分别为1500、2500、 2000、2500、2800、3000、4000、5000、4000、5000m/s。观测系统的参数为：最小 炮检距为40m，道间距为40m，最大炮检距为4000m。

图4为利用本发明的迭代后的初值公式计算的反射波时间与精确值的对比，图4 （a）为第2个界面的反射波时间的初值与真值的对比图（两条曲线基本重合），图4 （b）为第10个界面的反射波时间的初值与真值的对比图（两条曲线前部分基本重合， 后部分稍有区分），可见二者很接近，差别较小，表明了初值公式精度很高。图5 为在相同精度条件下，本发明的创新方法与传统方法（试射法）平均迭代次数对比， 可见本发明的方法的迭代次数很少即可达到给定精度。图6为本发明的方法与传统方 法（试射法）的CPU时间对比，可见本发明的方法计算时间明显少于传统方法。

本发明实施例采取以上技术方案，由于初值公式较为精确、迭代公式能够快速逼 近真解，因而具有计算速度快、精度高、始终收敛的优点。

在另外一个实施例中，还提供了一种软件，该软件用于执行上述实施例及优选实 施方式中描述的技术方案。

在另外一个实施例中，还提供了一种存储介质，该存储介质中存储有上述软件， 该存储介质包括但不限于：光盘、软盘、硬盘、可擦写存储器等。

显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以 用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算 装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而， 可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于 此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或 者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例 不限制于任何特定的硬件和软件结合。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特 定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包 括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一 个方框或多个方框中指定的功能。这些计算机程序指令也可装载到计算机或其他可编 程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计 算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程 图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例 的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术 人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述， 本说明书内容不应理解为对本发明的限制。