CN103698809A - 一种无加速比瓶颈的克希霍夫叠前时间偏移并行方法 - Google Patents

Abstract

本发明是石油地球物理勘探地震数据处理采用的适合大规模、无加速瓶颈的克希霍夫叠前时间偏移并行方法，将输出空间细分为线号、炮检距等四维，根据最小炮检距、最大炮检距、炮检距增量参数计算炮检距个数，把所有参与计算的节点进行分组，每组内的节点数目均匀，同一组内的节点在节点名的顺序连续，每一组节点一次计算一个炮检距组分的叠前时间偏移。本发明实现地震数据最小量网络传递，彻底解决了克希霍夫叠前时间偏移的网络传递密集型问题，可以使计算硬件规模扩展，使叠前时间偏移地计算周期缩短。

Description

一种无加速比瓶颈的克希霍夫叠前时间偏移并行方法

技术领域

本发明涉及石油地球物理勘探地震数据的处理技术，是一种适合大规模、无加速瓶颈的克希霍夫叠前时间偏移并行方法。

背景技术

克希霍夫叠前时间偏移技术是地球物理勘探地震数据处理过程中最重要的一种时间域成像方法，该方法不但具有成像质量高、对速度模型要求低、不受野外数据观测方式的限制等优点，另外它还能输出基于炮检距的共成像点道集，这类道集又是进一步速度分析重要的信息。

在水平地表情况下，Kirchhoff积分法偏移公式为：

$I\left(\mathrm{\ξ}\right)={\∫}_{{\mathrm{\Ω}}_{\mathrm{\ξ}}}W(\mathrm{\ξ},m,h)D[t=t_D(\mathrm{\ξ},m,h),m,h]\mathrm{dmdh}---\left(1\right)$

其中成像点ξ=(x_ξ,y_ξ,z_ξ)，对于时间偏移为ξ＝(x_ξ,y_ξ,τ_ξ)，I(ξ)是成像结果，D[t,m,h]是野外观测的地震数据，m是共中心点，h是半炮检距，Ω_ξ是偏移孔径。叠前偏移过程就是对一系列观测数据进行加权求和，W(ξ,m,h)是加权因子。t_D(ξ,m,h)是由炮点到成像点再到接收点的旅行时间。

克希霍夫叠前时间偏移方法的数学表达虽然较为抽象，但其物理实现过程可简单描述为：地震道从等待偏移数据空间到偏移结果数据空间多对多的映射，映射关系为双平方根旅行时计算公式，在不考虑时变孔径的条件下，克希霍夫叠前时间偏移就是把每一道等待偏移地震数据（炮检距为h）映射到偏移结果数据集中的一个椭圆柱体中，圆柱体的轴心是当前等待偏移地震道按大地坐标在偏移结果空间中的投影位置。若偏移孔径采用圆形且孔径值为300，则一个等待偏移地震采样点将要映射到282600（3.14*300*300）个偏移结果地震采样点上。

在水平地表情况下，克希霍夫积分法偏移用公式为前面所述公式⑴。⑴式可以进行深度偏移，也可以完成时间偏移，当完成时间偏移式时，旅行时用图1图示，⑵式是计算公式。

$t=t_s+t_r=\sqrt{\frac{t_0^2}{4}+\frac{h_s^2}{v_{\mathrm{rms}}^2}}+\sqrt{\frac{t_0^2}{4}+\frac{h_r^2}{v_{\mathrm{rms}}^2}}---\left(2\right)$

公式中变量参看图1。

结合公式⑴和⑵可以看出，完成一个地震样点的叠前时间偏移要包括旅行时计算、振幅加权、积分求和等计算，另外还有反假频滤波没有在公式中表示出来。一般地，一个地震道包含几千个地震样点、一个工区包含几千万至几亿个地震道，所以一个工区叠前时间偏移的计算量是极其庞大的，传统计算机的计算能力无法满足叠前时间偏移计算量的要求。直到计算机集群集成技术的发展和普及后，地球物理勘探才开始研究并行克希霍夫叠前时间偏移方法，现在已形成了以下两类算法：

一是输入道并行算法，这种方法以观测面上包含的面元总数为单位把整个偏移结果数据空间平均分配到参与计算的各个节点的内存中，待偏移数据空间中的每一道地震数据都要传送到各个计算节点上，所有等待偏移的地震数据都偏移完成后，再把各个节点内存中的偏移结果合并成完整的偏移结果。主要缺点是网络间数据的传输量和节点个数成正比，对硬件的可扩展性受到了限制。

二是输出道并行算法，这类方法一次完成输出空间沿观测线方向一条线（Inline）的偏移，因为一条Inline线的数据量不算太大，因此可以在每一个参与计算节点的内存中为整条线申请空间。把等待偏移数据空间中对当前Inline线有贡献的所有数据按参与计算的节点个数进行平均分割，一个节点完成其中的一份数据的偏移计算任务。偏移结束后把各个节点内存中的偏移结果做求和操作并输出，当前Inline线偏移结束，再转入下一条输出线偏移。主要缺点是地震数据的读取量、网络传输数量和输出空间inline线的个数以及参与并行计算的节点个数成一定的比例关系，随着偏移工区增大、节点个数增大，地震数据的网络流量快速增长。

目前实际生产中叠前时间偏移会按照两种方式运行，一种是目标线偏移，目的是优化偏移速度和偏移参数试验等，另一种是体偏移，其目的是得到最终的叠前时间偏移结果。目标线偏移计算量较小采用输出道并行算法效率较高，体偏移计算量较大采用输入道并行算法效率较高。

不管是输入道并行算法还是输出道并行算法，随着计算节点个数增加地震数据网络流量将随之增大。当计算机集群带宽不能满足叠前时间偏移程序地震数据的网络流量要求时，叠前时间偏移程序就会遇到加速比瓶颈，这时再增加计算的节点个数，偏移程序的运行时间反而会增加、运行效率降低。上述两种方法都存在加速比瓶颈，国际权威机构统计的叠前时间偏移程序加速比瓶颈会在节点个数90附近出现。

发明内容

本发明的目的是提供一种硬件可扩展、无加速比瓶颈的克希霍夫叠前时间偏移并行方法。

本发明通过以下具体步骤实现：

1）野外采集得到的地震数据，进行静校正、去噪、反褶积后对该数据进行预处理；

步骤1）所述的预处理过程是：

（1）收集每一道数据的炮点位置X坐标、炮点位置Y坐标、检波点位置X坐标、检波点位置Y坐标、沿测线方向（Inline）线号，垂直测线方向（Xline）线号信息；

（2）把等待偏移地震数据按照炮检距（offset）、测线方向线号、垂直测线方向线号按照三级关键字分选，将每一个炮检距数据保存为一个独立的地震数据文件；

（3）建立（1）中的道头字信息和（2）中地震数据文件相互对应关系。

2）根据最小炮检距、最大炮检距、炮检距增量参数，按下式计算炮检距个数：炮检距个数=（最大炮检距-最小炮检距）/炮检距增量；

3）把所有参与计算的节点进行分组，每组内的节点数目均匀，每组内的节点个数不小于4，同一组内的节点在节点名的顺序连续；

4）根据炮检距个数和节点个数进行并行计算分配，每一组节点一次计算一个炮检距组分的叠前时间偏移；

步骤4）所述的任务分配按以下三种情况进行：

（1）当计算节点的分组个数等于叠前时间偏移的炮检距个数，每一组节点计算一个炮检距叠前时间偏移，当某一组节点完成计算后，其组内节点自动加入到当前进度最慢的一组协同计算；

（2）当计算节点的分组个数小于叠前时间偏移的炮检距个数，每一组节点计算一个炮检距叠前时间偏移，当某一组节点完成计算后，申请下一个炮检距的计算；

（3）当计算节点的分组个数大于叠前时间偏移输出的炮检距个数，每一组节计算一个炮检距叠前时间偏移，剩余节点组依次分配加入其它组协同计算。

5）将步骤3）分组的第一个节点确定为组长（master）节点，组内其它节点作为组员（worker）节点；

6）根据偏移孔径确定当前成像区间所需要的待偏移地震数距的范围，具体方法为：

待偏移地震数距的最小沿测线方向线号=当前成像区间的最小沿测线方向线号-（垂直测线方向的最大孔径/线间距）；

待偏移地震数距的最大沿测线方向线号=当前成像区间的最大沿测线方向线号+（垂直测线方向的最大孔径/线间距）；

7）每一组的组长节点逐道读取本炮检距待偏移地震数据并按需分发到本组所有组员组长和组员节点；

8）每一个节点独立完成各自分配到的克希霍夫叠前时间偏移计算任务；

9）重复步骤7）~8）完成一个炮检距的克希霍夫叠前时间偏移；

10）每一组的组长节点回收该组组员节点的叠前时间偏移结果并进行叠加完成该组炮检距数据的克希霍夫叠前时间偏移。

本发明有以下特点：

⑴地震数据网络流量小，解决了传统方法存在的加速比瓶颈问题，可以使并行规模无限制扩大。

⑵由于引入了共炮检距数据集的概念，从数据流方面考虑叠前时间偏移和常用的叠后时间偏移的概念比较接近，使叠前时间偏移的并行任务分配方式变得非常简单。

⑶将平均任务分配模式变为“能者多劳”式任务分配模式，更适合于异步通讯机制，从而解决了慢节点或故障节点对整体计算性能的影响。

本发明预先把地震数据按炮检距进行分组，使各组地震数据之间没有依赖关系，因此减少了传统并行克希霍夫叠前时间偏移方法地震数据在网络上大量的冗余传递，从而解决了网络传递密集型问题，可以使计算硬件规模扩展，使叠前时间偏移地计算周期缩短。

附图说明

图1水平地表旅行时计算图；

图2本方明并行任务分配图；

图3本发明的脉冲响应图；

图4输入道并行方法的脉冲响应图；

图5本发明和输入道并行方法地震数据网络流量对比图；

图6目标线偏移、体偏移运行效率对比测试结果图，

左上：目标线偏移运行时间图、右上：目标线偏移加速比曲线图；

左下：体偏移运行时间图、右下：体偏移加速比曲线图。

具体实施方式

以下结合附图具体说明本发明。

本发明的一种无加速比瓶颈的克希霍夫叠前时间偏移并行方法，从算法本身上大大减少了地震数据在计算节点之间的传输总量，彻底解决了困扰叠前时间偏移并行算法的网络密集型问题。

本发明通过以下具体步骤实现：

1）野外采集得到的地震数据，进行静校正、去噪、反褶积后对该数据进行预处理；

步骤1）所述的预处理过程是：

（1）收集每一道数据的炮点位置X坐标、炮点位置Y坐标、检波点位置X坐标、检波点位置Y坐标、沿测线方向（Inline）线号，垂直测线方向（Xline）线号信息；

（2）把等待偏移地震数据按照炮检距（offset）、测线方向线号、垂直测线方向线号按照三级关键字分选，将每一个炮检距数据保存为一个独立的地震数据文件；

（3）建立（1）中的道头字信息和（2）中地震数据文件相互对应关系。

2）根据最小炮检距、最大炮检距、炮检距增量参数，图2所示，按下式计算炮检距个数：

炮检距个数=（最大炮检距-最小炮检距）/炮检距增量

落在某一个炮检距区间的所有等待偏移的数据仅仅对该炮检距的偏移结果数据有贡献，不会贡献到其它偏移结果炮检距上。据此把等待偏移的数据和偏移结果数据按照炮检距区间进行分组，数据组之间没有任何依赖关系。

3）把所有参与计算的计算机节点名编辑到一个ascii文件中，每行一个节点名，形成节点名文件。对该文件中包含的所有节点进行分组。分组原则是：每组内的节点数目均匀；每组内的节点个数不小于4；同一组内的节点在节点名文件中的顺序应该连续。

4）根据步骤2）的炮检距个数和步骤3）的节点个数进行并行计算分配，每一组节点一次计算一个炮检距组分的叠前时间偏移；

步骤4）所述的任务分配按以下三种情况进行：

（1）当计算节点的分组个数等于叠前时间偏移的炮检距个数，每一组节点计算一个炮检距叠前时间偏移，当某一组节点完成计算后，其组内节点自动加入到当前进度最慢的一组协同计算；

（2）当计算节点的分组个数小于叠前时间偏移的炮检距个数，每一组节点计算一个炮检距叠前时间偏移，当某一组节点完成计算后，申请下一个炮检距的计算；

（3）当计算节点的分组个数大于叠前时间偏移输出的炮检距个数，每一组节计算一个炮检距叠前时间偏移，剩余节点组依次分配加入其它组协同计算。

5）根据步骤3）的计算节点分组情况，把每一个节点组内的第一个节点确定为组长（master）节点，组内其它节点作为组员（worker）节点；

6）根据偏移孔径确定当前成像区间所需要的待偏移地震数距的范围，具体方法为：

待偏移地震数距的最小沿测线方向线号=当前成像区间的最小沿测线方向线号-（垂直测线方向的最大孔径/线间距）；

待偏移地震数距的最大沿测线方向线号=当前成像区间的最大沿测线方向线号+（垂直测线方向的最大孔径/线间距）；

7）每一组的组长节点逐道读取本炮检距待偏移地震数据并按需分发到本组所有组员组长和组员节点；

8）每一个节点独立完成各自分配到的克希霍夫叠前时间偏移计算任务；

9）重复步骤7）~8）完成一个炮检距的克希霍夫叠前时间偏移；

10）每一组的组长节点回收该组组员节点的叠前时间偏移结果并进行叠加完成该组炮检距数据的克希霍夫叠前时间偏移。

本发明原理如下：

设：待偏移数据空间

$N=\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{j=1}{\overset{J_i}{\mathrm{\Σ}}}{N}_{\mathrm{ij}}---\left(3\right)$

其中：N----待偏移数据总大小

I----待偏移数据Inline线数

J_i----第i条Inline线的cmp个数

N_ij----第i条线第j个cmp道集内地震数据大小

设：偏移结果数据空间

$M=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}m---\left(4\right)$

其中：M----偏移结果数据总大小

K----偏移结果Inline线数

L----每条Inline线的crp个数

m----每个crp面元内地震大小

设：节点个数为I，总计算量为P

则，一个节点上要完成的计算量

P_i＝P/I    ⑸

设：传统的输入道并行方法的网络流量为T

则T=M*I+N    ⑹

从⑸式可以看出，单个节点上的计算量和节点个数成反比，随着节点个数的增加，单个节点上的计算量将减少，整体计算速度将提高。

从⑹式可以看出，网络流量和节点个数成正比，随着节点个数的增加，网络流量将正比例增加，数据传送的开销将要增加。当结点个数多于某个数值时，集群的网络资源将不能满足并行程序数据传输量的需求，使得节点的计算单元不能及时得到数据处于等待状态，从而降低整体计算速度。

综合分析⑸式和⑹式，输入道并行方法虽然在工业界发挥着重要的作用，当计算节点增加到一定个数时，增加的网络网络流量严重超越了集群网络资源传送能力，这就限制了输入道并行方法只能在中小规模的并行环境中运行。同样输出道方法也存在相同的问题，这里不再论述。

本发明的一种无加速比瓶颈的克希霍夫叠前时间偏移并行方法把输出空间细分为四维空间，四个坐标轴分别指的是：Inline线号、Xline线号、成像时间、炮检距。从⑷式可以看出，Inline线数和Xline线数是地震工区的成像范围，由用户编码提供，在地面上的投影是规则的矩形网格，同时每个CRP道集中包含的地震道数也是由用户填写成像炮检距区间和炮检距增量参数组确定。

设：MinOffset---最小炮检距

MaxOffset---最大炮检距

IncOffset---炮检距增量

则每个CRP道集中的地震道数

m=（MaxOffset-MinOffset）/IncOffset+1⑺

其中第i道的炮检距为

Offset_i=(i-0.5)*IncOffset+MinOffset    ⑻

把四维空间中相同道号或相同炮检距的所有数据抽取为一组，整个输出空间将变成一系列共道号或共炮检距的三维数据集，如图2下部所示。

本发明的一种无加速比瓶颈的克希霍夫叠前时间偏移并行方法以图2中的共道号或共炮检距数据集为单位进行并行任务划分，由一个或几个节点共同完成一个共炮检距数据集的三维叠前时间偏移的计算工作。

由于其它炮检距的待偏移数据不会给该炮检距的输出结果产生贡献，所以只需要在本组节点上传送该组炮检距的待偏移数据。

设：每个节点组内的节点个数为n，n由实现步骤3）确定

则本发明方法下地震数据的网络流量为：

T＝N+M*n    ⑼

比较⑹式和⑼式，一般来说待偏移数据N和偏移结果数据M的大小基本相当，⑹式中节点个数I远大于⑼式中n，所以⑼式表示的网络流量要远远小于⑹式的网络流量，随着节点个数增加，⑹式的网络流量成线性关系递增，⑼式网络流量保持不变。

本发明为了验证效果，选用一个理论数据进行了单道脉冲偏移对比试验，图3是本发明的脉冲响应，图4是输入道并行方法的脉冲响应，两附图中的脉冲响应完全相同，可以肯定本发明方法是完全正确的。

本发明为了验证效果，选用某一实际三维工区，工区待偏移数据大小174GB，在dn8集群上进行不同并行算法地震数据网络流量测试。dn8集群包含128个节点，本次测试使用了60个节点，前30个节点运行本发明的无加速比瓶颈的叠前时间偏移并行算法，后30个节点运行输入道叠前时间偏移并行算法。两个测试运行相同规模的体偏移作业，且同时发送。

图5展示的是本发明方法和输入道并行方法网络流量测试对比效果，其中每一个小图代表一个单独节点的资源使用情况，纵坐标为使用率，变化范围为0%~100%，横坐标为运行时间。黄色线条代表网络流量、红色线条代表CPU使用效率、绿色线条代表内存的使用情况。

从图中可以看出前面节点代表的本发明方法的地震数据网络流量远小于下部节点代表的输入道并行方法。所以本发明成功地解决了克希霍夫叠前时间偏移的网络密集型问题，仅剩下计算密集型，计算密集型通过增加计算节点就能解决。

本发明为了验证效果的硬件环境：测试集群有128个计算节点，每个节点有4个核。

测试数据：测试的地震数据是某三维工区，工区范围：Inline6~Inline1066，Xline139~Xline1479，目标线测试范围：Inline45~Inline265~Step55，即从45线起每隔55条线进行偏移总共偏移5条线。体偏移范围：Inline100~Inline200~Step1，总共偏移101条线。

图6是分别用本发明方法和输入道平行算法测试的目标线偏移、体偏移的运行时间和加速比。图6上部是目标线偏移的运行时间和加速比，运行时间计时单位为秒，蓝色代表输入道并行算法，红色代表本发明方法。图6下部是体偏移的运行时间和加速比，运行时间计时单位为小时，同样蓝色代表输入道并行算法，红色代表本发明方法。以16个节点为最小测试作业，按倍数关系依次测试16、32、64、128个节点的运行时间，并计算相应的加速比。

从图6上部的目标线偏移测试结果可以看出本发明方法随着节点个数的增加，运行时间基本按倍数关系减少，加速比曲线非常接近于理论值。输入道并行方法随着节点个数的增加，程序的运行时间在32节点处达到最小，再增加节点个数时，程序运行时间将不再减少，加速比曲线在32节点处达到最高值，继续增加节点时，该曲线严重偏离理论值。

从图6下部的体偏移测试结果可以看出，本发明方法和输入道并行方法当节点个数的增加，运行时间都在减少，很明显每一组相同节点个数作业的对比测试，本发明方法运行时间都小于输入道并行方法，但是不同组间的减少关系不是很明确。但从右边的加速比曲线可以看出，本发明的方法更接近于理论加速比。

Claims (3)

1.一种无加速比瓶颈的克希霍夫叠前时间偏移并行方法，特点是通过以下具体步骤实现：

1）野外采集得到的地震数据，进行静校正、去噪、反褶积后对该数据进行预处理；

2）根据最小炮检距、最大炮检距、炮检距增量参数，按下式计算炮检距个数：炮检距个数=（最大炮检距-最小炮检距）/炮检距增量；

3）把所有参与计算的节点进行分组，每组内的节点数目均匀，每组内的节点个数不小于4，同一组内的节点在节点名的顺序连续；

4）根据炮检距个数和节点个数进行并行计算分配，每一组节点一次计算一个炮检距组分的叠前时间偏移；

5）将步骤3）分组的第一个节点确定为组长节点，组内其它节点作为组员节点；

6）根据偏移孔径确定当前成像区间所需要的待偏移地震数距的范围，具体方法为：

待偏移地震数距的最小沿测线方向线号=当前成像区间的最小沿测线方向线号-（垂直测线方向的最大孔径/线间距）；

待偏移地震数距的最大沿测线方向线号=当前成像区间的最大沿测线方向线号+（垂直测线方向的最大孔径/线间距）；

7）每一组的组长节点逐道读取本炮检距待偏移地震数据并按需分发到本组所有组员组长和组员节点；

8）每一个节点独立完成各自分配到的克希霍夫叠前时间偏移计算任务；

9）重复步骤7）~8）完成一个炮检距的克希霍夫叠前时间偏移；

10）每一组的组长节点回收该组组员节点的叠前时间偏移结果并进行叠加完成该组炮检距数据的克希霍夫叠前时间偏移。

2.根据权利要求1的方法，特点是步骤1）所述的预处理过程是：

（1）收集每一道数据的炮点位置X坐标、炮点位置Y坐标、检波点位置X坐标、检波点位置Y坐标、沿测线方向线号，垂直测线方向线号信息；

（2）把等待偏移地震数据按照炮检距、测线方向线号、垂直测线方向线号按照三级关键字分选，将每一个炮检距数据保存为一个独立的地震数据文件；

（3）将以上步骤（1）中信息和以上步骤（2）中地震数据文件相互对应。

3.根据权利要求1的方法，特点是步骤4）所述的任务分配按以下三种情况进行：

（1）当计算节点的分组个数等于叠前时间偏移的炮检距个数，每一组节点计算一个炮检距叠前时间偏移，当某一组节点完成计算后，其组内节点自动加入到当前进度最慢的一组协同计算；

（2）当计算节点的分组个数小于叠前时间偏移的炮检距个数，每一组节点计算一个炮检距叠前时间偏移，当某一组节点完成计算后，申请下一个炮检距的计算；

（3）当计算节点的分组个数大于叠前时间偏移输出的炮检距个数，每一组节计算一个炮检距叠前时间偏移，剩余节点组依次分配加入其它组协同计算。

Images