CN109541677A - 一种保护激发频率的延时震源组合参数计算方法 - Google Patents

Abstract

本发明涉及一种保护激发频率的延时震源组合参数计算方法，该方法首先需要建立地质模型，在后期优化可能用到的坐标位置上布设震源，在以最浅目的层深度为半径组合中心为球心的球面上等角度布设检波器，完成所有震源的单炮正演；然后在给定组合个数和组合形态的基础上通过数值寻优算法或贪婪算法计算得到最佳组合参数，计算过程需要满足给定的主频约束条件。本发明克服了当前震源组合参数设计上的某些不足，寻找到了一种快速的组合参数定量化计算方法，同时具有良好地可操作性，可以为野外施工提供方案指导。

Description

一种保护激发频率的延时震源组合参数计算方法

技术领域

本发明涉及延时震源组合参数计算方法，具体涉及一种保护激发频率的延时震源组合参数计算方法，属于石油勘探领域。

背景技术

我国的陆地地震勘探中，炸药震源仍然是经常采用的能量源。炸药震源可视为点震源，在均匀介质中震源激发的能量向四周均匀扩散，波前面为一球面，非均匀介质存在波的额外反射，传播波场复杂，能量分布与介质结构有关。无论地下介质均匀与否，有一点是确定的，那就是地震勘探所接收到的反射能量仅仅只是原始震源激发能量的一小部分。炸药震源激发的能量很大一部分用于破碎周围岩石，这部分由岩石的力学性质决定无法人为控制，一小部分用于产生原始弹性震动能量，这部分能量中的又一小部分向地下传播到达地层产生反射波并最终被地表的检波器接收，形成地震记录。弹性能量的更大一部分直接传至地表或向两翼传播不产生反射，在降低能量利用效率的同时，会产生额外的杂乱反射叠加进有效波场，污染有效反射，降低记录信噪比。

因此，想要提升采集质量，就需要考虑改变弹性能量的空间配比，将更多的能量传导到地下地层。实践中改变能量空间分布的方法就是组合激发。

组合满足波场的线性叠加原理，各震源子波传播到空间某位置时具有不同的到达时间，时间同步则等相位叠加，能量加强，时间不同步则叠加能量减弱，相差1/2周期的简谐波叠加抵消为零。人为控制波场能量分配的原理就是设计合理的时间搭配关系，在需要的空间方向上加强能量，其它方向削弱能量。

组合能量的空间分布一般称为方向特性，组合频率的空间分布被称为频率特性。造成频率改变的原因本质上也是不同简谐波成分的叠加削弱与加强。总体而言，震源子波中的低频成分由于具有更长的周期，到达时间的一致性较容易得到满足，或者说对时间的不同步不敏感，与之相对应，具有短周期的高频成分则很容易被叠加削弱，从而导致整个频谱在组合后向低频端移动，也就是我们常说的组合具有降频效应。

目前生产中常用的震源组合方式为等深度的平面组合，这些施工方法存在诸多不足：

1、没有一个定量化的评价组合优劣的方法，组合好坏完全是试验多种组合后对比看效果；

2、没有考虑组合对频率的影响；

3、组合深度上的变化较少考虑；

4、延时激发较少采用。

发明内容

为了解决上述技术问题，本发明提供一种保护激发频率的延时震源组合参数计算方法，可以快速的将震源组合参数定量化，为野外施工提供设计方案指导。本发明的技术方案具体如下：

一种保护激发频率的延时震源组合参数计算方法，包括以下步骤：

步骤(1)、根据工区实际情况，建立三维地质模型；

步骤(2)、基于三维地质模型，确定待优化的组合中心，围绕该中心在三维地质模型的地质网格节点上布设一定数量的震源；

步骤(3)、获取最浅目的层深度信息

以步骤(2)确定的组合中心为球心，以最浅目的层深度为半径做球面，该球面与地质网格相交，在球面上以等倾角和水平角度间隔布设检波器，以获取最浅目的层深度信息；

步骤(4)、对步骤(2)中所有震源依次进行单程波动方程正演，记录下各检波器接收到的首波波形，存储待用；

步骤(5)、根据检波器接收排列长度和最浅勘探目标层的深度确定有效入射角α，根据组合激发图形中心点的深度和最浅勘探目标层的深度确定非有效入射角β；

步骤(6)、统计球面上投射至目的层锥角范围内的能量和其它方向能量的比值，以能量比值最大化为判别标准进行优化计算

根据步骤(5)得到的有效入射角α和非有效入射角β，将给定震源组合个数和平面组合形态作为已知量，统计球面上投射至目的层锥角范围内的能量和其它方向能量的比值，以能量比值最大化为判别标准，采用粒子群优化算法或贪婪算法计算延时震源组合参数，延时震源组合参数包括各个震源的三维空间坐标(x_i,y_i,z_i)、组合加权系数C_i和组合激发延迟时间△t_延迟i。

进一步地，步骤(2)中，震源布设的数目与范围取决于优化所需的震源坐标取值区间。

进一步地，步骤(5)中，有效入射角α按下式计算：

其中，L是检波器接收排列长度，H是最浅勘探目标层的深度；

非有效入射角β按下式计算：

其中，hs是组合激发图形中心点的深度(中心点的Z坐标)或者β＝90°-α。

进一步地，步骤(6)中，延时震源组合参数的计算按下式进行：

其中，θ为极坐标系下的倾角、为水平角，t为合成子波的传播时间，α为指向目的层的有效入射角度，β为指向目的层之外区域的非有效入射角度，T为某角度上从震源激发零时刻开始至合成子波波尾穿过检波器球面位置所用的时间，固定θ和求和的核函数f表达式为：

f_合(t)＝c₁f₁(t+Δt₁)+c₂f₂(t+Δt₂)+…+c_nf_n(t+Δt_n)，

，其中，C_i为各震源的加权系数，△t_i为各震源子波到达观测位置的时间差。

进一步地，步骤(6)中，△ti按下式计算：

Δt_i＝Δt_传播i+Δt_廷迟i，

其中，

式中，d_i为各震源到观测点的空间距离，v为介质速度，(x_i，y_i，z_i)为各震源点的空间坐标。

进一步地，步骤(6)中，延时震源组合参数的计算还需要进行频率约束，使得有效入射角度范围内所有角度上组合后子波主频不低于某一给定值，组合后某角度上的子波在频域表达式如下：

其主频位置可以由该角度下所有频率的谱能量最大值定位，再找出所有角度下的主频的最小值，让其满足大于给定频率即可，最小主频由下式定位：

与现有技术相比，本发明的有益效果具体如下：

(1)本发明可根据激发层及目的层地质信息设计合理的震源组合参数，提升有效信号能量的同时将其频率的降低控制在给定的范围内，从而在提高地震记录的信噪比的同时不损失分辨率，为后期的处理解释工作提供高质量的原始资料。

(2)本发明采用单程声波波动方程正演获取准确的空间地震波形，以最浅一套目的层的深度为半径，以组合中心为球心，统计该球面上投射至目的层锥角范围内的能量和其它方向能量的比值，以能量比值最大化为判别标准进行优化计算，同时满足该角度范围内主频降低不大于给定值。

(3)本发明考虑了介质非均质性、组合的空间位置关系、组合延迟时间因素、组合的频率效应，最重要的是给出了定量化的评价标准，填补了定量化震源组合参数设计技术的空白，同时也为野外施工提供了可靠的依据。

(4)本发明可根据简化的地质模型、激发子波主频率、观测系统排列长度、最浅目的层深度等信息快速计算出最优的震源组合方式，优化后的组合方式包含了震源在平面、深度、药量、激发时间上的变化。

附图说明

图1为实施例的震源子波空间叠加原理示意图；

图2为实施例的震源激发能量分区示意图；

图3为实施例的极坐标系下倾角和水平角的方向规定示意图，图3(a)为倾角，取值范围[-180°～180°]，图3(b)为水平角，取值范围[0°～360°]；

图4为实施例的波动方程正演计算时震源与检波器布设空间位置示意图；

图5为本实施例的工艺流程图；

图6为波动方程正演记录，抽取沿X轴的剖面(0°水平角)，-90°～90°倾角范围的检波器接收数据显示。图6(a)表示当震源位于组合中心位置时到各检波器距离相等，子波到达时间一致，图6(b)表示当震源偏离组合中心，与各检波器距离不等，子波到达时间波动起伏。

图7为组合后主频的空间分布示意图。图7(a)为三维显示，图7(b)为沿X轴切片二维显示。

图8为本实施例的9口震源组合的平面分布形态；

图9为本实施例的优化后的空间展布形态。

具体实施方式

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

如图1所示，由于满足线形叠加原理，各震源在观测位置的叠加子波表达式如式(1)所示：

f_合(t)＝c₁f₁(t+Δt₁)+c₂f₂(t+Δt₂)+…+c_nf_n(t+Δt_n) (I)

对其做傅里叶变换，得到频域表达式，如式(2)所示：

根据式(1)可以计算合成子波的能量，根据式(2)可以计算合成子波的频谱。

其中，C_i为各震源的加权系数(可近似正比于炸药量)，△t_i为各震源子波到达观测位置的时间差，其中包含因各震源空间位置差异导致的传播时间差和激发时间不同导致的延迟时间差两部分，如式(3)所示：

Δt_i＝Δt_传播i+Δt_延迟i (3)

激发延迟时间项采用生产中时间采样的整倍数(如1ms采样时，延迟时间取为-10ms、+3ms、+8ms等，负号代表提前激发，正号代表延后激发)，可以作为参数优化计算中的一个独立参数。

传播时间项与各震源到观测点位置的空间距离、介质速度结构有关，在匀速介质中传播时间可用式(4)表达：

式中d_i为各震源到观测点的空间距离，v为介质速度，(x_i，y_i，z_i)为各震源点的空间坐标，可见空间位置不同会引起传播时间不同，进而影响最终叠加合成的子波形态，空间坐标(x_i，y_i，z_i)也可以作为本发明组合参数优化计算中的三个独立参数。

非均匀介质中激发点到观测点的传播路径非直线无法用式(4)简单表述，可以采用波动方程正演记录的直达波初至时间。

可见，线形叠加原理下，影响最终合成子波能量和频率的只有加权系数和组合时间差两个因素，组合时差又与激发时间、震源空间位置有关，所以组合参数设计中待优化计算的参数一共包含x_i、y_i、z_i、△t_延迟i、c_i五个。

本实施例建立在如下的假设基础上：

(1)不考虑破碎区的影响，各震源的初始子波波形一致，为雷克子波；

(2)组合基距尺寸范围内介质均质；

(3)传播能量的衰减不考虑地层吸收因素的影响。

实施例1

本实施例涉及的计算过程均采用C++实现。

如图5所示，本实施例的保护激发频率的延时震源组合参数计算方法，包括如下步骤：

步骤(1)、根据工区实际情况，建立三维地质模型；

步骤(2)、基于三维地质模型，确定待优化的组合中心，围绕该中心在地质网格节点上布设一定数量的震源，震源布设的数目与范围取决于优化所需震源坐标取值区间，如设定组合中心坐标(0,0,0)，后期需要在X∈[-10m～10m]，Y∈[-10m～10m]，Z∈[-10m～10m]范围内寻优计算，地质网格间隔若为1m，则需要围绕组合中心布设11*11*11个震源；

如图8所示，本实施例以(0,0,0)作为组合中心坐标，设置9个震源，震源分布在围绕组合中心的X∈[-20m～20m]，Y∈[-20m～20m]，Z∈[-20m～20m]空间范围内，提取沿X轴剖面图，得到震源激发能量分区示意图如图2所示。

步骤(3)、获取最浅目的层深度信息

以步骤(2)确定的组合中心为球心，以最浅目的层深度为半径做球面，该球面与地质网格相交，在球面上以等倾角和水平角度(比如等1°)间隔布设检波器，以获取最浅目的层深度信息，如图4所示。

步骤(4)、对步骤(2)中所有震源依次进行单程波动方程正演，记录下各检波器接收到的首波波形，存储待用；

本实施例中，波动方程正演记录如图6所示，抽取沿X轴的剖面(0°水平角)，-90°～90°倾角范围的检波器接收数据显示。图6(a)表示当震源位于组合中心位置时到各检波器距离相等，子波到达时间一致，图6(b)表示当震源偏离组合中心，与各检波器距离不等，子波到达时间波动起伏。

图7为组合后主频的空间分布示意图，从图7可以看到，主频从垂直方向至水平方向逐渐降低。图7(a)为三维显示，图7(b)为沿X轴切片二维显示。

步骤(5)、根据检波器接收排列长度和最浅目的层深度信息确定有效入射角α，根据组合中心深度和最浅目的层深度确定非有效入射角β；其中，如图2所示，弧AB的张角范围(有效入射角)α按下式计算：

其中，L是检波器接收排列长度，H是最浅勘探目标层深度；

弧AD与BC的张角(非有效入射角)β大小相等，非有效入射角β按下式计算：

其中，hs是组合激发图形中心点的深度(中心点的Z坐标)。

由于震源深度相对于目的层埋深一般较浅，一般可简化为β＝90°-α。

步骤(6)、统计球面上投射至目的层锥角范围内的能量和其它方向能量的比值，以能量比值最大化为判别标准进行优化计算。

图2中，最大偏移距决定了地下地层的最大覆盖范围(EF区域)，只有在该范围内传播的能量才与有效反射相关，在剖面上该范围被OE与OF射线围限，射线与球面交于AB两点，二维情况下分别统计弧AB上接收的波形能量，统计弧AD与BC上接收的波形能量，二者相除，比值最大即为优化判别函数(即统计球面上投射至目的层锥角范围内的能量和其它方向能量的比值，以能量比值最大化为判别标准。)，三维情况下需要将各弧沿水平方向0°～360°累加，同样做能量比值，如式(5)所示：

其中，θ为极坐标系下的倾角、为水平角，t为合成子波的传播时间，α为指向目的层的有效入射角度，β为指向目的层之外区域的非有效入射角度，T为某角度上从震源激发零时刻开始至合成子波波尾穿过检波器球面位置所用的时间。有效入射角α和非有效入射角β根据步骤(5)得到，给定震源组合个数和平面组合形态作为已知量。

三个求和变量依次为以组合中心为坐标系原点的极坐标系下垂直倾角θ、水平方位角和传播时间t，倾角和水平角的正方向规定如附图3所示，传播时间的最大值T为各角度上子波波尾穿过球面的时间。

本实施例中，固定θ和为任意值(即在确定某一空间方向的情况下)，各震源在观测位置的叠加子波作为求和的核函数f表达式如下：

f_合(t)＝c₁f₁(t+Δt₁)+c₂f₂(t+Δt₂)+…+c₉f₉(t+Δt₉)

Δti＝Δt_传播i+Δt_延迟i

△t_传播i参照公式(4)，i指参与组合的第i个震源编号。

组合参数优化计算需要进行频率约束，保证有效入射角度范围内所有角度上组合后子波主频不低于某一给定值。本实施例组合后某角度上的子波在频域表达式如下：

频率约束满足有效入射角度范围内所有角度上组合后子波主频不低于某一给定值，某一角度下的子波主频位置可以由该角度下所有频率的谱能量最大值定位，再找出所有角度下的主频的最小值，让其满足大于给定频率即可，最小主频由式(6)定位：

如图8、图9所示，本实施例给定震源组合个数和平面组合形态作为已知量，以震源的三维空间坐标(x_i,y_i,z_i)、组合加权系数C_i、组合激发延迟时间△t_延迟i构成一组解向量，采用粒子群优化算法优化组合参数，目标函数为式(5)。

期间利用式(6)判断该解在有效入射角的范围内是否满足频率约束条件，若满足，记录下该解，并继续进行下一步操作，若不满足，放弃该解，继续进行下一步操作。

(7)重复步骤(6)达到设定的迭代次数或者计算精度，输出优化计算结果。

本实施例中，优化计算时采用30Hz雷克子波作为子波源，组合基距内介质纵波速度800m/s。本实施例中，组合个数和组合形态(比如圆形、十字交叉形、三角形等)不作为优化变量，直接给出，在特定组合数目特定平面组合形态的基础上优化各组元的空间坐标、加权系数、激发延迟时间参数。

所有优化变量的求解区间如表1，计算时各优化变量从解区间内随机生成一组任意值(或者由人为按照经验设定一组值)作为初值，按照步长变化寻优。

表1优化变量的求解区间

最终得到的优化结果如表2所示：

表2优化结果

实施例2

本实施例的保护激发频率的延时震源组合参数计算方法，步骤(6)中采用贪婪算法优化组合参数。

贪婪算法能保证找到全局最优解但计算速度慢，粒子群优化算法计算时间可控，但不能保证找到全局最优解，若给定逼近最优解的经验值作为初始输入，计算效果较好。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种保护激发频率的延时震源组合参数计算方法，其特征在于：包括以下步骤：

步骤(1)、根据工区实际情况，建立三维地质模型；

步骤(2)、基于三维地质模型，确定待优化的组合中心，围绕该中心在三维地质模型的地质网格节点上布设一定数量的震源；

步骤(3)、获取最浅目的层深度信息

以步骤(2)确定的组合中心为球心，以最浅目的层深度为半径做球面，该球面与地质网格相交，在球面上以等倾角和水平角度间隔布设检波器，以获取最浅目的层深度信息；

步骤(4)、对步骤(2)中所有震源依次进行单程波动方程正演，记录下各检波器接收到的首波波形，存储待用；

步骤(5)、根据检波器接收排列长度和最浅勘探目标层的深度确定有效入射角α，根据组合激发图形中心点的深度和最浅勘探目标层的深度确定非有效入射角β；

步骤(6)、统计球面上投射至目的层锥角范围内的能量和其它方向能量的比值，以能量比值最大化为判别标准进行优化计算

根据步骤(5)得到的有效入射角α和非有效入射角β，将给定震源组合个数和平面组合形态作为已知量，统计球面上投射至目的层锥角范围内的能量和其它方向能量的比值，以能量比值最大化为判别标准，采用粒子群优化算法或贪婪算法计算延时震源组合参数，延时震源组合参数包括各个震源的三维空间坐标(x_i,y_i,z_i)、组合加权系数C_i和组合激发延迟时间△t_延迟i。

2.根据权利要求1所述的延时震源组合参数计算方法，其特征在于：步骤(2)中，震源布设的数目与范围取决于优化所需的震源坐标取值区间。

3.根据权利要求1所述的延时震源组合参数计算方法，其特征在于：步骤(5)中，有效入射角α按下式计算：

其中，L是检波器接收排列长度，H是最浅勘探目标层的深度；

非有效入射角β按下式计算：

其中，hs是组合激发图形中心点的深度或者β＝90°-α。

4.根据权利要求1所述的延时震源组合参数计算方法，其特征在于：步骤(6)中，延时震源组合参数的计算按下式进行：

其中，θ为极坐标系下的倾角、为水平角，t为合成子波的传播时间，α为指向目的层的有效入射角度，β为指向目的层之外区域的非有效入射角度，T为某角度上从震源激发零时刻开始至合成子波波尾穿过检波器球面位置所用的时间，固定θ和求和的核函数f表达式为：

f_合(t)＝c₁f₁(t+Δt₁)+c₂f₂(t+Δt₂)+…+c_nf_n(t+Δt_n)，

其中，C_i为各震源的加权系数，△t_i为各震源子波到达观测位置的时间差。

5.根据权利要求4所述的延时震源组合参数计算方法，其特征在于：步骤(6)中，△ti按下式计算：

Δt_i＝Δt_传播i+Δt_延迟i，

其中，

式中，d_i为各震源到观测点的空间距离，v为介质速度，(x_i，y_i，z_i)为各震源点的空间坐标。

6.根据权利要求4所述的延时震源组合参数计算方法，其特征在于：步骤(6)中，延时震源组合参数的计算还需要进行频率约束，使得有效入射角度范围内所有角度上组合后子波主频不低于某一给定值，组合后某角度上的子波在频域表达式如下：

其主频位置可以由该角度下所有频率的谱能量最大值定位，再找出所有角度下的主频的最小值，让其满足大于给定频率即可，最小主频由下式定位：