CN102778690A

CN102778690A - 一种基于混合基离散傅氏变换的波动方程叠前偏移性能优化方法

Info

Publication number: CN102778690A
Application number: CN2011101238653A
Authority: CN
Inventors: 杨祥森; 赵改善; 亢永敢; 刘百红
Original assignee: China Petroleum and Chemical Corp; Sinopec Geophysical Research Institute
Current assignee: China Petroleum and Chemical Corp; Sinopec Geophysical Research Institute
Priority date: 2011-05-13
Filing date: 2011-05-13
Publication date: 2012-11-14
Anticipated expiration: 2031-05-13
Also published as: CN102778690B

Abstract

本发明提供了一种基于混合基DFT的波动方程叠前偏移性能优化方法，属于地球科学中的油气勘探地球物理新技术的应用研究领域。所述方法首先通过分析基于混合基离散傅氏变换函数的变换时间与变换长度的关系得出相对最优变换长度及其计算时间，形成资源文件；然后根据用户给定的参考扩边数和单炮成像范围，用查找顺序表和资源文件的方式求取最优变换长度；最后通过调整波场延拓的扩边范围，使二维傅氏变换时间相对最小，从而提高频率-波数域波动方程叠前偏移的计算效率。本发明减低了二维付氏变换时间，提高了单炮偏移的计算效率，且不降低成像质量。虽然本发明是针对SSF延拓算子提出的，但完全可以扩展到所有频率—波数域单程波叠前偏移算法中。

Description

一种基于混合基离散傅氏变换的波动方程叠前偏移性能优化方法

技术领域

本发明属于地球科学中的油气勘探地球物理新技术的应用研究领域，具体涉及一种基于混合基离散傅氏变换（DFT）的波动方程叠前偏移性能优化方法，能在确保偏移成像质量的情况下，有效提高波动方程叠前偏移的串行或并行计算效率，缩短勘探地震数据的处理周期。

背景技术

近几年来，随着油气勘探开发工作的深入研究和计算机软硬件技术的高速发展，基于波动方程的叠前深度偏移越来越受到油气勘探开发业界的注目。由于波动方程叠前深度偏移方法能细致的描述波场在介质中的传播特征、适用于强横向变化介质的成像，且其成像效果优于Kirchhoff积分法，但其巨大的计算量和数据量也是业界必须面临的现实问题。

随着高性能计算机在石油勘探行业的应用，基于波动方程的偏移技术研究则越来越关注于并行计算，以提高计算性能。目前，偏移算法并行计算方面的研究主要集中于并行域的划分、任务动态分配、多文件并行I/O、编译参数优化和基于GPU环境的并行设计等。这些研究都是针对计算机的软、硬件系统和提高并行性的。但是，要提高波动方程叠前偏移的计算效率，首要任务是最大化优化串行算法。

目前，波动方程叠前偏移方法包括PS、PSPI、SSF、FD、FFD、GSP等。这些算子都是为了在确保成像质量前提下提高计算速度，或者在确保计算速度的前提下提高成像质量，是针对不同的地质介质，采用高频近似理论推导出来的。除了FD波场延拓算子是在时间—空间域计算外，其他算子都是在频率—波数域或频率—波数、时间—空间混合域中进行偏移计算的。另外，为解决计算速度，一些学者还使用多种近似计算、简化方法对原有算法进行改进，以提高计算速度，但同时也降低了成像质量。

对于频率—波数域波场延拓算子，波场延拓计算必须要进行二维正、逆傅氏变换，并且二维傅氏变换时间是偏移计算的计算“热点”或重要环节。因此，减少二维傅氏变换时间可以提高叠前偏移的计算效率。目前，研究较多的傅氏变换是基4FFT、基8FFT、分裂基FFT、素因子FFT、维诺格兰傅氏变换和算术傅氏变换等，这些算法都是为减少蝶形乘法运算，加快变换速度，但都增加了计算复杂度或存在其他缺陷，因此使用最多的傅氏变换为基2快速傅氏变换（基2FFT）。但是，当变换长度不满足2的幂次时，基2FFT需要补0，这就增加了计算量和内存需求。同时基2FFT无法对指定频率点进行DFT，且存在频谱泄露和栅栏效应，而混合基FFT刚好能解决基2FFT的这些缺陷。在中等和大型转换方面，Intel®MKL库中所有混合基DFT函数都已进行了高度优化，为多处理器系统提供了优秀的可伸缩性，其性能比FFTW（快速计算离散傅里叶变换的标准C语言程序集）更出色。因此它在频谱保真和计算速度上具有明显的优势，已逐渐得到广泛的应用。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种基于DFT的波动方程叠前偏移性能优化方法，以SSF波场延拓算子为例，根据Intel®MKL库中基于混合基DFT函数的变换时间与变换长度的关系，提高二维傅氏变换和整个偏移成像作业的计算效率。

本发明是通过以下技术方案实现的：

一种基于混合基离散傅氏变换的波动方程叠前偏移性能优化方法，所述方法首先通过分析基于混合基离散傅氏变换（DFT）函数的变换时间与变换长度的关系得出相对最优变换长度及其计算时间，形成资源文件；然后根据用户给定的参考扩边数和单炮成像范围，用查找顺序表和资源文件的方式求取最优变换长度；最后通过调整波场延拓的扩边范围，使二维傅氏变换时间相对最小，从而提高频率-波数域波动方程叠前偏移的计算效率。

所述方法包括以下步骤：

（1）读炮集，获得偏移参数和偏移数据；

（2）计算当前炮的成像范围（NX、NY），其中，NX、NY分别为单炮覆盖的CDP数和线数；

（3）扩边，根据参考扩边范围（NXP、NYP）求得偏移计算范围（NXX、NYY），其中，NXX、NYY分别为扩边后波场延拓的参考计算范围；

（4）判断与前一炮计算范围是否相等，如果相等，则转入步骤（10），如果不相等，则转入步骤（5）；

（5）向前或向后查顺序表，查找与此炮对应的最优变换长度；

（6）判断是否在顺序表中找到对应的最优变换长度，如果没有，则转入步骤（7），如果查到，则转入步骤（8）；

（7）从资源文件中查找最优变换长度值，并将找到的最优变换长度值插入到顺序表中；

（8）求得新的计算范围（NNX、NNY），即优化后的偏移计算范围；

（9）反算出实际的扩边数（NXP’、NYP’），其中，NXP’、NYP’分别为实际扩边CDP数和线数；

（10）偏移成像，分别用新的计算范围（NNX、NNY）和实际扩边范围（NXP’、NYP’）进行波场延拓和边界衰减。

所述顺序表为顺序数组或链表，每插入一组值就进行排序，插入的值包括成像范围和偏移计算范围。

所述资源文件是通过分析基于DFT函数的变换时间与变换长度的关系得出的相对最小变换时间及其对应的变换长度。

与现有技术相比，本发明的有益效果是：

（1）本发明主要针对Intel ® MKL库中DFT特点，通过调整SSF单程波叠前深度偏移的扩边范围，减少二维傅氏变换时间，提高单炮偏移的计算效率；

（2）虽然本发明是针对SSF延拓算子提出的，但完全可以扩展到所有使用二维傅氏变换的频率—波数域单程波叠前偏移算法中；

（3）由于用户给定的扩边数是由试验测试得到，一般能达到成像要求，而实际的扩边数往往大于等于参数扩边数，所以，实际扩边数并不会降低成像质量，甚至还会提高成像质量，同时降低了二维付氏变换时间，提高了单炮偏移的计算效率。

附图说明

图1 是本发明二维付氏变换长度与DFT变换时间关系图（ny=1024,nx=8~4096），图中，第二维长度定为1024，第一维长度从8递增到4096。其中横轴为第一维的变换长度，纵轴为变换时间；圆实心点为所有变换长度的计算时间，方实心点是变换长度为2的n次幂（n≥1）的计算时间，圆空心点为计算时间相对最小的变换长度及其计算时间。这些最优变换长度值及其计算时间将写入资源文件，供偏移计算时查表使用。

图2 是图1的局部放大图（ny=1024,nx=1998~2007），图中第二维长度定为1024，第一维长度从1998递增到2007。图中所示实心圆点（nx=2000）是以nx=1998为参考变换长度、搜索半径为10得到的变换时间相对最小的变换长度。

图3 是本发明中波场延拓范围示意图，图中，NX、NY分别为单炮覆盖的CDP数（共深度点）和线数，中间矩形区域为单炮偏移的成像范围；NXP、NYP为偏移计算的扩边CDP数和线数，周边区域为扩充区域，也称为波场衰减区域，但周边四个角落矩形区域与上下左右四个边界矩形区域的衰减系数不同；NXX、NYY为扩边后波场延拓的参考计算范围。

图4 是本发明数据采集过程中炮号与单炮覆盖范围关系示意图，每个工区中变观察系统采集的炮数一般较少；对于变观察系统采集的炮集，其覆盖范围一般呈逐渐变宽或变窄的趋势。

图5 是本发明方法的实施步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

如图1所示，一种基于混合基DFT的波动方程叠前偏移性能优化方法，所述方法先设定一个顺序表（顺序表为顺序数组或链表，每插入一组值（包括成像范围和偏移计算范围）要进行排序，顺序表的设置是为了提高查表速度），计算当前炮的成像范围NX、NY，并与前一炮的成像范围比较（对于第一炮，直接将其成像范围和偏移计算范围插入顺序表中，不用比较）；若与前一炮成像范围相等，则直接使用前一炮的偏移计算区域值和实际扩边数，否则在顺序表中向前或向后查找与之对应的最优变换长度；若在顺序表中未找到最优变换长度，则在资源文件中查找最优变换长度，并将找到的值插入到顺序表中。

如图5所示，所述方法包括以下步骤：

（1）读炮集，炮集包括道头和道数据，本步骤是获得偏移参数和偏移数据；

（3）扩边，根据参考扩边范围（NXP、NYP）求得偏移计算范围（NXX、NYY），其中，NXX、NYY分别为扩边后波场延拓的参考计算范围（也就是扩边后的CDP数和线数）；

（9）反算出实际扩边数（NXP’、NYP’），其中，NXP’、NYP’分别为实际扩边CDP数和线数；

（10）偏移成像，分别用新的计算范围（NNX、NNY）和实际扩边范围（NXP’、NYP’）进行波场延拓和边界衰减。（偏移成像即用偏移成像算子进行波场延拓与成像，本文采用SSF算子，具体的成像方法为本领域技术人员公知）。

所述资源文件是通过分析基于混合基DFT函数的变换时间与变换长度的关系得出的，具体分析如下：相邻的变换长度值的付氏变换时间差异很大，且在一个相对小的数值范围内总存在一个变换长度值，其对应的傅氏变换计算时间最短，如图2中的圆实心点是变换长度为1998~2007区间内变换时间最小的点。因此，本发明经过实际测试，找出相对最小的变换时间及其变换长度值（如图1中的圆空心点）作为资源文件，供后续的性能优化工作使用。

为减少边界效应，在波场延拓过程中一般都是对成像区域进行扩边（扩边范围如图3中的(周边矩形区域），每延拓一步就对扩充区域进行波场衰减。扩充区域必须足够大，否则将影响成像质量。为了提高二维付氏变换的计算效率，根据用户给定的参考扩边数（如图3中的NXP、NYP（NXX=NX+2*NXP、NYY=NY+2* NYP,可以从图3中看出））和成像区域（如图3中的中心矩形区域，即NX×NY）大小，用二次查表法求取最优二维傅氏变换的变换长度，并反算出实际扩边数（即NXP’、NYP’）；然后用实际扩边数和成像范围进行波场延拓和边界衰减。由于用户给定的扩边数是由试验测试得到，一般能满足成像精度要求，而实际的扩边往往大于等于参数扩边数。所以，实际扩边数并不会降低成像质量，甚至还会提高成像质量，同时降低了二维付氏变换时间，提高了单炮偏移的计算效率。

在勘探地震数据实际采集过程中，由于每个工区中变观察系统（变观察系统是相对于规则观测系统的地震数据采集观测系统,属地震勘探专有词汇，本领域技术人员公知）采集的炮数较少，且单炮覆盖范围一般呈逐渐变宽或变窄的趋势，如图3所示，所以本发明在求取最优变换长度值时采用二次查表法（即一查顺序表，二是查资源文件），将大大减少最优变换长度值检索次数，加快最优变换长度值求取速度。

由于用户（指开发的波动方程叠前偏移程序的使用者）给定的扩边数是由试验测试得到，一般能满足成像精度要求，本方法中求取的实际扩边数大于或等于参考扩边数，所以优化后的偏移技术不会降低成像质量，甚至还会提高成像质量，同时降低了二维付氏变换时间，提高了单炮偏移的计算效率。

本发明主要是根据NXX和NYY的值，用二次查表法找到最优变换长度NNX、NYY，然后反算出实际扩边数NXP’、NYP’。用NNX、NYY，NXP’、NYP’进行波场延拓和波场边界衰减计算。其中查顺序表和查资源文件即为二次查表过程，它是鉴于勘探地震数据的实际采集技术中单炮采集覆盖范围与炮号的关系设计的，可以减少最优变换长度值的查找次数，加快搜索速度。

下面通过两个表来说明本发明的效果：采用本发明的方法，表1、表2分别为SEG/EAEG盐丘模型数据优化前后的六项测试时间统计表，其中计算时间的单位是秒。表1、表2中，Sno为炮集的炮号，ntr为炮集的道数，sht为当前炮的计算时间，ipt为当前进程的总计算时间，nx、ny分别为扩边前的CDP点数和线数，即成像范围，nxx、nyy分别为扩边后的CDP点数和线数（即二维傅氏变换的第一、二维的变换长度），DFTt、PSt、PERt、IMSt、TAPt、CRSt分别为二维傅氏变换（包括正反变换）、相移、扰动项、成像与叠加、边界处理的计算时间和频率域波场数据拷贝时间。对比两表可以看出，各炮的计算效率都有不同程度的提高，其中599炮和603炮的计算效率提高最多，最高达2.8倍。

表1

表2

本发明主要针对Intel ® MKL库中混合基DFT特点，通过调整SSF单程波叠前深度偏移的扩边范围，减少二维傅氏变换时间，提高单炮偏移的计算效率。虽然本发明是针对SSF延拓算子提出的，但是对于三维地震成像，空间域变换到波数域都要用到二维傅氏变换，而本发明是针对傅氏变换做的性能优化，因此本发明的方法完全可以扩展到所有频率—波数域单程波叠前偏移算法中。

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

Claims

1.一种基于混合基离散傅氏变换的波动方程叠前偏移性能优化方法，其特征在于：所述方法首先通过分析基于混合基离散傅氏变换函数的变换时间与变换长度的关系得出相对最优变换长度及其计算时间，形成资源文件；然后根据用户给定的参考扩边数和单炮成像范围，用查找顺序表和资源文件的方式求取最优变换长度；最后通过调整波场延拓的扩边范围，使二维傅氏变换时间相对最小，从而提高频率-波数域波动方程叠前偏移的计算效率。

2.根据权利要求1所述的基于混合基离散傅氏变换的波动方程叠前偏移性能优化方法，其特征在于：所述方法包括以下步骤：

（1）读炮集，获得偏移参数和偏移数据；

3.根据权利要求1或2所述的基于混合基离散傅氏变换的波动方程叠前偏移性能优化方法，其特征在于：所述顺序表为顺序数组或链表，每插入一组值就进行排序，插入的值包括成像范围和偏移计算范围。

4.根据权利要求1或2所述的基于混合基离散傅氏变换的波动方程叠前偏移性能优化方法，其特征在于：所述资源文件是通过分析基于混合基离散傅氏变换函数的变换时间与变换长度的关系得出的相对最小变换时间及其变换长度。