2017, Vol. 32
2. 黑龙江省普通高校科技创新团队"断层变形、封闭性及流体运移", 黑龙江大庆 163318
3. 大庆钻探工程公司地球物理勘探一公司, 黑龙江大庆 163357
2. Heilongjiang Province College Science and Technology Innovation Team'Faults Deformation, Sealing Ability and Fluid Migration', Heilongjiang Daqing 163318, China
3. NO.1 Geophysical Exploration Company Daqing Drilling & Exploration Engineering Corporation, Heilongjiang Daqing 163357, China
逆时偏移方法因其能对复杂构造进行精确成像而得到关注，并被逐渐应用于实际地震数据处理中(Baysal et al., 1983；Whitmore，1983).当采用互相关成像条件时，逆时偏移的成像需要对震源波场和检波点波场进行同时刻互相关计算.震源波场可通过波场正演模拟获得，即沿时间方向正向传播，而检波点波场则通过对检波点信息的逆时传播来获取，因此，若要对上述两种波场进行任意时刻的互相关计算，则至少需要重建或保存一个时间方向上的全部波场信息.为此，众多学者针对逆时偏移波场的边界条件、存储策略以及重建方法展开了研究探讨.Symes(2007)提出了一种优化checkpoint方法，通过周期性设置checkpoint，保存部分中间时刻的波场信息，并以中间波场为节点，重建所需时刻的波场，该方法可有效降低波场存储量，但实现过程较为复杂，且需引入较多的额外计算.Dussaud等(2008)分析了当前几种逆时偏移存储策略，指出checkpoint方法最为节省存储空间，但会增加额外计算.Clapp(2008)提出了一种波场存储策略，结合有限差分方法，保存空间差分阶数一半层数的内部边界波场，即可实现震源波场的完全重建，该方法适用于任何吸收边界条件，且在保证波场重建精度的基础上，有效降低了波场存储量.Clapp(2009)提出了一种随机边界条件，通过在速度场外围引入随机速度层，可有效散射边界层的入射波，由于波场能量没有损失，因此仅需最后两个时刻的震源波场信息即可实现波场重建，大幅降低了波场存储量，但引入的随机噪声会降低逆时偏移结果的信噪比(李博等，2010；刘红伟等，2010).Feng和Wang(2012)提出了保存边界处单层波场的存储策略，降低了波场存储量，但由于震源波场正向传播和重建时采用了不同的差分格式，波场重建时会引入误差.王保利等(2012)提出了有效边界存储策略，借鉴checkpoint和缓存思想，降低了波场存储，且不会引入计算误差，但该方法会增加额外的计算量.结合GPU加速计算，可显著减少额外的计算耗时对算法计算量的影响，但并没有针对引入的缓冲区空间，给出参数选取的最优策略.Yang等(2014)将Clapp(2008)提出的波场存储策略应用于GPU(Graphic Processing Unit)平台进行逆时偏移计算，在获得较高计算效率的同时，也缓解了GPU端的波场存储压力，取得了较好的应用效果.综上所述，虽然逆时偏移的波场重建策略取得了长足的进步，但每种方法均存在不同程度的应用限制或没有提供相应方法的最优参数选取方案.针对该方向的研究，仍有很大的发展空间和重要的现实意义.
基于二阶常密度声波波动方程，本文提出了一种空间变阶数有限差分方法，在计算内部有效波场及上边界波场时采用高阶差分算子，计算左、右和下边界附近的波场时采用变阶数差分算子，差分阶数由高阶逐级降低至二阶，同时将这种变阶数差分应用于震源波场的正向传播和重建过程.存储策略方面，上边界保存差分阶数一半的边界层波场，左、右和下边界保存单层波场，从而保证了正传和重建震源波场的精度.为进一步降低波场存储量，借鉴缓冲区思想，利于对算法进行GPU加速，同时，给出了缓冲区空间的最优化策略，推导得出令波场存储量最小化的参数选取.通过对二维盐丘模型进行波场重建及逆时偏移测试，验证了本文方法的有效性.
1 理论方法 1.1 变阶数有限差分二阶常密度声波波动方程为
|
(1) |
式中，v(x)为速度场，U(x, t; xs)为波场，f(x, t; xs)为震源函数.其中，x为空间坐标，对应二维空间即为x=(x, z)；xs为震源坐标；t为时间.
常规的有限差分法离散波动方程表达式为
|
(2) |
式中，Pm, nk为离散形式下的波场，上标k代表时间采样点，下标m、n分别代表水平方向和垂直方向的空间离散采样点，Δt表示时间步长，Δx、Δz分别为空间水平方向和垂直方向网格间距，Cl为有限差分系数，L表示空间差分阶数的一半.
采用变阶数有限差分方法对(1)式离散可得到：
|
(3) |
与(2)式的不同在于，式(3)中的M、N均为空间差分阶数的一半，但不为恒定值，而是随空间网格点位置的变化而变化，D和E分别为变差分阶数情况下水平方向和垂直方向的有限差分系数，当L=M=N时，则有C=D=E.详细的空间差分阶数的变化分布如图 1所示.假设给定的初始波场范围如图 1中黑色实线框所示，该范围内的波场即为有效波场.为压制波场的边界反射，需引入一定厚度的边界层，如图 1中黑色虚线框与黑色实线框之间的区域，该区域的波场即为边界波场，通常可以采用吸收边界条件来实现对边界反射的吸收衰减(Liu and Sen, 2010；王维红等，2013).若常规计算区域的差分阶数为10阶，如图 1中灰色短实线所示，则变差分计算区域及差分阶数变化如图 1中彩色短实线表示.在有效波场内部及上边界区域，采用常规10阶有限差分格式，以保证波场的计算精度；在左、右及下边界区域，差分阶数沿边界外法线方向逐渐降低，直至有效波场边界处内侧时，差分阶数降为2阶，边界层内差分阶数的分布则延续有效波场边界处的差分格式，保持不变.
|
图 1 变阶数有限差分 Figure 1 Sketch of variable order finite difference |
依据网格点位置，震源波场可划分为两部分：有效波场(如图 1中黑色实线框内区域)和边界波场(如图 1中黑色虚线框和黑色实线框之间的区域).在有效波场范围内，需准确计算波场信息以满足成像精度；在边界波场范围内，进行指定边界条件的计算，以消除波场的边界反射的影响，本文采用的是PML(Perfectly Matched Layer)边界条件(Berenger，1994).逆时偏移方法中，成像区域仅限于有效波场内部，对有效波场范围内震源波场的精确重建是实现高精度成像的前提.
假设有效波场水平方向网格点数为Nx，垂直方向网格点数为Nz，边界层数为B，空间最大差分阶数为2L.震源波场正向传播过程中，保存所有时刻中有效波场范围内的上边界L层(以10阶差分为例，仅保存上边界5层波场)，左、右和下边界的单层波场，如图 2中黑色方块所示位置，以及最后两个时刻的全部有效波场.
|
图 2 波场存储策略 Figure 2 Sketch of wavefield storage strategy |
震源波场重建时，首先以最后两个时刻的有效波场作为初始条件，载入预先保存的对应时刻的有效边界波场作为边界条件.对公式(3)中Pm, nk+1和Pm, nk-1项进行互换，即可得到波动方程逆时传播的有限差分格式为：
|
(4) |
应用式(4)进行波场重建，所采用的差分阶数的分布如图 2中灰色及彩色实线所示(与图 1对应位置处的差分阶数分布相同)，由此即可完成有效波场范围内震源波场的逆时重建.
为进一步说明本文方法的有效性，以有效波场右边界为例，结合式(4)，用k时刻波场Pm, nk推算k-1时刻波场Pm, nk-1时，有效波场右边界处水平方向空间差分计算方式如图 3所示(假设最高差分阶数为10阶)，以有效波场边界处保存的单层波场(图 3中黑色方块)边界条件和计算得出的有效波场的内部波场值(图 3中紫色方块)为已知数据，空间差分阶数由外向内逐层提高，即可满足式(4)中的空间差分计算.因此，在震源波场反向传播时，只需保存有效波场上边界处内侧的5层波场值与左、右和下边界处的单层波场值，根据图 2所示的差分阶数分布，即可保证重建所得的有效波场(图 3中)与对应时刻的初始震源波场相同.
|
图 3 水平方向有效波场区域内变阶数有限差分计算 Figure 3 Sketch of variable order finite difference calculation within the valid wavefield region in horizontal direction |
目前，GPU加速技术在地震数据处理领域已取得了广泛应用(刘国峰等，2009；石颖等，2010；郭雪豹等，2016；Shi and Wang, 2016；柯璇和石颖，2017).GPU并行加速计算技术的引入极大的提升了地震波场数值模拟的计算效率，一定程度上缓解了逆时偏移方法的计算压力，降低了逆时偏移应用的计算成本.但受GPU显存容量的限制，有效边界波场存储策略所需的存储量仍然不可忽视，因此有必要进一步改进存储策略，以降低GPU显存设备的存储压力.
前文所述内容主要利用变阶数有限差分法和有效边界存储策略，在空间域降低存储量.然而在时间方向上，仍需要保存每一时刻的波场信息.为了压缩时间上的波场存储量，本文借鉴checkpoint思想，通过设置存储缓冲区，在时间方向上进一步降低波场存储需求.具体实现过程为：
(1) 首先在时间方向设置N个checkpoint，在进行震源波场正向延拓过程中，保存相应checkpoint处包含边界区域的震源波场信息(如图 2中黑色虚线框内所有网格点处的波场).另外设置缓冲区，如图 4所示，用于保存两个相邻的checkpoint之间的n个有效边界震源波场(如图 2中黑色方块所示位置的波场).当震源波场传播至如图 4中第N-1个checkpoint后，将第N-1至第N个checkpoint间的有效边界震源波场保存至缓冲区中.
|
图 4 Checkpoint+缓冲区存储策略 Figure 4 Sketch of checkpoint and buffer storage strategy |
(2) 当震源波场正向延拓至最大时刻后，根据第N个checkpoint中的波场，采用1.2节所述方法，利用(4)式计算逆时传播的波场，同时载入对应时刻的缓冲区中的有效边界震源波场，即可逆时重建出第N-1至第N个checkpoint间的有效震源波场.
(3) 再以第N-2个checkpoint中的波场信息为初始条件，进行震源波场的正向延拓，将第N-2至第N-1个checkpoint间的有效边界震源波场覆盖至缓冲区中，然后根据第N-1个checkpoint中的波场和缓冲区中的有效边界震源波场，利用(4)式即可逆时重建出第N-2至第N-1个checkpoint间的有效震源波场.检波点波场的反传计算与上述过程同时进行，并与对应时刻震源波场进行互相关计算，直至完成成像运算.
基于以上方法，避免了对所有时刻的有效边界震源波场的保存，可进一步降低波场存储量，但参数N和n的选取，直接影响checkpoint波场和缓冲区的存储占用量.因此，本文给出一种最优化参数选取方法.
定义时间方向采样点数为tmax，假设速度模型的水平和垂直网格点数分别为Nx和Nz，外边界层数为B，如图 2所示，则在checkpoint中需保存的一个含边界的震源波场切片的网格点数为(Nx+2B)(Nz+2B).由于本文采用时间方向2阶差分格式，每个checkpoint中需保存两个时刻的波场切片才能作为震源波场正向延拓的初始条件.为方便表示，以网格点数为波场存储量计量单位，则每个checkpoint中所需的波场存储量为2(Nx+2B)(Nz+2B)，而保存所有checkpoint波场需分配的存储空间为2N(Nx+2B)(Nz+2B).如图 2所示，可知保存在缓冲区中的一个有效边界震源波场切片需要存储的网格点数为(L+1)Nx+2(Nz－L-1)，不难推算，缓冲区所需占用的存储则为n[(L+1)Nx+2(Nz－L-1)].
根据前文所述，可知N和n满足关系：
|
(5) |
则波场存储量占用可由下式得到：
|
(6) |
由(6)式可知，波场存储量是一个关于N和n的函数.利用(5)式可设g(N, n)=tmax－Nn，将其作为限制条件，通过拉格朗日乘子法构建目标函数为
|
(7) |
将(6)式和g(N, n)=tmax－Nn代入(7)式可得：
|
(8) |
(8) 式为关于N、n和λ三个自变量的函数，分别对三个变量求偏导，并令其为0，可得：
|
(9) |
求解上式可得：
|
(10) |
根据(10)式选取N和n的值，即可实现对checkpoint和缓冲区的最优化分配.在计算过程中，N和n需取整数.因此，定义〈〉为四舍五入的取整运算，则当n=
采用二维盐丘模型测试波场重建方法，速度模型如图 5所示，水平方向网格点数为649，垂直方向网格点数为150，水平和垂直方向上的空间网格间距均为10 m，时间采样率为0.8 ms，震源为主频25 Hz的雷克子波，在地表 3.5 km处激发.
|
图 5 二维盐丘速度模型 Figure 5 2D salt velocity model |
当震源波场传播至0.48 s时，波场快照如图 6a所示，震源波场正传至最大时刻4.8 s后，采用1.2节所述方法，对震源波场进行逆时重建，逆时传播至0.48 s时的震源波场如图 6b所示.对正传波场和逆时重建波场进行误差分析，结果如图 6c所示，误差振幅与有效波场振幅相差6个数量级，可忽略不计，由此可知本文方法可实现波场的高精度重建.
|
图 6 波场重建结果 (a)0.48 s时正传震源波场；(b)0.48 s时重建震源波场；(c)0.48 s时重建震源波场误差. Figure 6 Results of wavefield reconstruction (a)Forward source wavefiled at 0.48 s; (b)Reconstructed source wavefield at 0.48 s; (c)The error of the reconstructed source wavefield at 0.48 s. |
本节对二维盐丘模型进行逆时偏移测试及存储策略分析，旨在验证变阶数差分方法的理论适用性和波场存储策略的有效性.偏移所采用的平滑速度场，如图 7所示，模型参数如2.1节所述.采用地表激发地表接收的观测方式，共450炮激发，炮间距为10 m，均匀分布在1.0~5.5 km范围内，如图 7中绿色实线所示；每炮649道接收，道间距10 m，检波器位置固定且均匀分布在0~6.49 km范围内，如图 7中蓝色虚线所示.时间采样率为0.8 ms，检波点记录时间采样点数为6000，时长为4.8 s.经Laplacian去噪后的逆时偏移剖面如图 8所示.
|
图 7 偏移平滑速度场 Figure 7 Smoothed velocity for reverse time migration |
|
图 8 逆时偏移结果 Figure 8 Result of reverse time migration |
本次测试中，采用1.2节所述有效边界存储策略和1.3节所述最优化有效边界存储策略，所需的波场存储量如表 1所示.由表中数据可知，在有效边界存储策略的基础上，通过本文的优化方法，存储量可降低约75%.
|
表 1 波场存储量对比 Table 1 Comparison of wavefield storage |
为定量分析本文方法对不同尺寸的模型应用逆时偏移方法时的波场存储量，假设模型水平方向和垂直方向网格点数相同(Nx=Nz)，外边界层数为50(B=50)，空间最大差分阶数为10阶(L=5)，时间方向采样点数为6000(tmax=6000).根据1.2和1.3节所述内容，可得出波场存储量随模型尺寸由100×100逐渐升至10100×10100的变化关系，如图 9所示.当模型网格点数为10010×10010时，最优化有效边界存储策略所需波场存储量约为320 Mb，仍小于目前绝大多数GPU设备的显存容量，由此可知本文方法具有较广的适用范围.
|
图 9 不同方法的波场存储量随模型尺寸变化 Figure 9 The variation of wavefield storage with model size for different methods |
地震资料逆时偏移方法对存储量的需求是限制其实用化的主要因素之一，同样的问题也存在于最小二乘逆时偏移、全波形反演等以逆时偏移为迭代核心的成像方法.本文提出了一种变阶数有限差分方法，上边界和内部波场均采用高阶差分以保证波场的计算精度，左、右及下边界处采用变阶差分，从而仅需保存上边界处的半差分阶数层波场和左、右、下边界处的单层边界波场即可实现波场逆时重建.在保证波场重建精度的同时，有效降低了波场存储量.同时，借鉴checkpoint和缓冲区思想，在时间方向上进一步压缩了有效边界存储策略的波场存储需求，并通过拉格朗日乘子法，给出了该存储策略中参数的最优化选取方案.
3.2二维盐丘模型的波场重建测试表明，本文所提出的变阶数差分方法，能够精确的重建波场，误差的量级在计算机系统误差范围内.同时，所提出的最优化有效边界存储策略，显著减少了波场存储量，降低了方法本身对硬件的需求，提升了逆时偏移方法的实用性.另外，本文中关于波场存储策略中参数最优化的推导方法，亦可应用于逆时偏移中的其他波场存储策略，具有广泛的适用性.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [] | Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration[J]. Geophysics, 48(11): 1514–1524. DOI:10.1190/1.1441434 |
| [] | Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics, 114(2): 185–200. DOI:10.1006/jcph.1994.1159 |
| [] | Clapp R G. 2008. Reverse time migration:Saving the boundaries[C].//Stanford Exploration Project, 136-144. |
| [] | Clapp R G. 2009. Reverse time migration with random boundaries[C].//79th Ann. Int. Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2809-2913, doi:10.1190/1.3255432. |
| [] | Dussaud E, Symes W W, Williamson P, et al. 2008. Computational strategies for reverse-time migration[C].//78th Ann. Int. Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2267-2271, doi:10.1190/1.3059336. |
| [] | Feng B, Wang Y H. 2012. Reverse time migration with source wavefield reconstruction strategy[J]. Journal of Geophysics and Engineering, 9(1): 69–74. DOI:10.1088/1742-2132/9/1/008 |
| [] | GUO Xue-Bao, LIU Hong, SHI Ying. 2016. Time domain full waveform inversion based on frequency attenuation[J]. Chinese Journal of Geophysics , 59(10): 3777–3787. DOI:10.6038/cig20161022 |
| [] | KE Xuan, SHI Ying. 2017. Forward simulation and reverse time migration imaging based on one step wavefield extrapolation[J]. Chinese Journal of Geophysics , 60(11): 4468–4479. DOI:10.6038/cjg20171131 |
| [] | LI Bo, LIU Hong-Wei, Liu Guo-Feng, et al. 2010. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU[J]. Chinese Journal of Geophysics , 53(12): 2938–2943. DOI:10.3969/j.issn.0001-5733.2010.12.017 |
| [] | LIU Guo-Feng, LIU Hong, LI Bo, et al. 2009. Method of prestack time migration of seismic data of mountainous regions and its GPU implementation[J]. Chinese Journal of Geophysics , 52(12): 3101–3108. DOI:10.3969/j.issn.0001-5733.2009.12.019 |
| [] | LIU Hong-Wei, LIU Hong, ZOU Zhen, et al. 2010. The problems of denoise and storage in seismic reverse time migration[J]. Chinese Journal of Geophysics , 53(9): 2171–2180. DOI:10.3969/j.issn.0001-5733.2010.09.017 |
| [] | Liu Y, Sen M K. 2010. A hybrid scheme for absorbing edge reflections in numerical modeling of wave propagation[J]. Geophysics, 75(2): A1–A6. DOI:10.1190/1.3295447 |
| [] | SHI Ying, LIU Hong, ZOU Zhen. 2010. Surface-related multiples prediction based on wave equation and adaptive subtraction investigation[J]. Chinese Journal of Geophysics , 53(7): 1716–1724. DOI:10.3969/j.issn.0001-5733.2010.07.023 |
| [] | Shi Y, Wang Y H. 2016. Reverse time migration of 3D vertical seismic profile data[J]. Geophysics, 81(1): S31–S38. DOI:10.1190/GEO2015-0277.1 |
| [] | Symes W W. 2007. Reverse time migration with optimal checkpointing[J]. Geophysics, 72(5): SM213–SM221. DOI:10.1190/1.2742686 |
| [] | WANG Bao-Li, GAO Jing-Huai, CHEN Wen-Chao, et al. 2012. Efficient boundary storage strategies for seismic reverse time migration[J]. Chinese Journal of Geophysics , 55(7): 2412–2421. DOI:10.6038/j.issn.0001-5733.2012.07.025 |
| [] | WANG Wei-Hong, KE Xuan, PEI Jiang-Yun. 2013. Application investigation of perfectly matched layer absorbing boundary condition[J]. Progress in Geophysics , 28(5): 2508–2514. DOI:10.6038/pg20130529 |
| [] | Whitmore N D. 1983. Iterative depth migration by backward time propagation[C].//53th Ann. Int. Mtg., Soc. Expi. Geophys. Expanded Abstracts, 382-385, doi:10.1190/1.1893867. |
| [] | Yang P L, Gao J H, Wang B L. 2014. RTM using effective boundary saving:A staggered grid GPU implementation[J]. Computers & Geosciences, 68: 64–72. DOI:10.1016/j.cageo.2014.04.004 |
| [] | 郭雪豹, 刘洪, 石颖. 2016. 基于频域衰减的时域全波形反演[J]. 地球物理学报, 59(10): 3777–3787. DOI:10.6038/cjg20161022 |
| [] | 柯璇, 石颖. 2017. 基于一步法波场延拓的正演模拟和逆时偏移成像[J]. 地球物理学报, 60(11): 4468–4479. DOI:10.6038/cjg20171131 |
| [] | 李博, 刘红伟, 刘国峰, 等. 2010. 地震叠前逆时偏移算法的CPU/GPU实施对策[J]. 地球物理学报, 53(12): 2938–2943. DOI:10.3969/j.issn.0001-5733.2010.12.017 |
| [] | 刘国峰, 刘洪, 李博, 等. 2009. 山地地震资料叠前时间偏移方法及其GPU实现[J]. 地球物理学报, 52(12): 3101–3108. DOI:10.3969/j.issn.0001-5733.2009.12.019 |
| [] | 刘红伟, 刘洪, 邹振, 等. 2010. 地震叠前逆时偏移中的去噪与存储[J]. 地球物理学报, 53(9): 2171–2180. DOI:10.3969/j.issn.0001-5733.2010.09.017 |
| [] | 石颖, 刘洪, 邹振. 2010. 基于波动方程表面多次波预测与自适应相减方法研究[J]. 地球物理学报, 53(7): 1716–1724. DOI:10.3969/j.issn.0001-5733.2010.07.023 |
| [] | 王保利, 高静怀, 陈文超, 等. 2012. 地震叠前逆时偏移的有效边界存储策略[J]. 地球物理学报, 55(7): 2412–2421. DOI:10.6038/j.issn.0001-5733.2012.07.025 |
| [] | 王维红, 柯璇, 裴江云. 2013. 完全匹配层吸收边界条件应用研究[J]. 地球物理学进展, 28(5): 2508–2514. DOI:10.6038/pg20130529 |