, 2. 中国石油塔里木油田分公司勘探开发研究院, 新疆库尔勒, 841000;
3. 中石化石油工程地球物理公司南方分公司, 四川成都, 610041
2. Research Institute of Exploration & Development, PetroChina Tarim Oilfield Company, Korla 841000, China;
3. Nangfang Branch of Sinopec Geophysical Corporation, Chengdu 610041, China
地震偏移成像旨在对地下构造的反射信号重新归位, 并根据地震资料刻画地下构造[1-6]。随着勘探要求的提升, 地震成像的目标从地层构造描述向地层属性描述转变[7-10]。将近年来快速发展的最小二乘逆时偏移算法与反演思想结合, 可用于精确描述地层属性。YAO等[11]提出了基于矩阵描述的最小二乘逆时偏移算法, 并指出该方法在消除低频噪声的同时更好地聚焦成像能量。郭振波等[12]提出了最小平方逆时偏移真振幅成像方法, 并验证了该方法在真振幅成像方面的明显优势。
目前, 针对最小二乘逆时偏移算法的研究主要集中于改善算法成像效果、提升算法适用性和计算效率等方面。ZHANG等[13]指出, 由于地下介质是一种变密度的粘弹性介质, 采用常规的声波方程模拟地震波场会导致振幅与实际情况匹配不佳, 因此提出了一种新的目标函数, 降低了振幅的不匹配对最优化算法的影响, 提升了算法的稳定性和数据的适应性。ZHANG等[14]提出了一种不依赖子波的最小二乘逆时偏移策略, 可有效降低由子波不匹配引起的噪声干扰。李庆洋等[15]提出了去均值归一化互相关最小二乘逆时偏移算法, 该方法修改了目标泛函, 利用去均值归一化算法, 降低了算法对子波能量的要求, 提升了算法的稳定性和可靠性。刘学建等[16]提出了表面多次波最小二乘逆时偏移算法, 并将多次波作为有效信号应用于成像算法, 增加了成像范围, 虽然初次迭代时产生了串扰, 但随着迭代次数的增加, 串扰得以消除。
最小二乘逆时偏移算法的迭代流程计算量大, 计算效率的提升对推动最小二乘逆时偏移算法发展至关重要。多炮数据同时计算是提升最小二乘逆时偏移算法效率的有效途径之一。BERKHOUT[17]定义了“面炮”偏移概念, 先对炮集数据合成叠加, 再进行偏移计算, 可有效降低偏移计算量, 该思路目前广泛应用于最小二乘逆时偏移算法; 平面波静态编码[18]、自适应奇异谱分析[19]和频率选择编码[20]等多种方法提升了最小二乘逆时偏移计算效率, 也有效压制了由炮集合成计算产生的串扰; 李闯等[21]推导了平面波最小二乘逆时偏移算法, 大幅提升了算法的执行效率, 改善了最小逆时偏移成像质量, 分析了多种编码策略, 总结了多种编码策略的优势。加快迭代误差下降速度的方法, 可降低计算量, 在迭代终止条件不变时, 可减少迭代计算次数, 以提升计算效率。LIU等[22]针对ZHANG等[13]提出的求取迭代步长参数的问题, 给出了解析步长(analytical step length, ASL)公式, 提升了迭代算法的误差下降效率; 预条件和规则化方法的应用, 有效加速了迭代算法的收敛, 并带来了计算效率的提升, 提高了深部成像分辨率和保幅性, 对于不规则的地震数据, 该方法有更好的适应性[23, 24]。GPU加速技术的发展, 从硬件方面提升了地震数据处理方法的计算效率:李博[25]、刘红伟[26]和SHI等[27]实现了基于GPU加速的逆时偏移算法; 石颖等[28, 29]将GPU加速技术应用于多次波的预测和衰减算法中; 郭雪豹等[30]采用GPU加速技术实现了基于频率衰减的全波形反演方法, 这为实现最小二乘逆时偏移算法的高性能计算带来新的契机。随着计算机技术的发展, 计算设备也逐渐升级, 如集群设备中, 单个高性能计算节点通常具备多核中央处理器(central processing unit, CPU)和多GPU, 但由于最小二乘逆时偏移算法相对复杂, 需频繁更新迭代参数, 从而导致了最小二乘逆时偏移算法仍缺乏一个相对完整的多GPU加速解决方案。
本文提出了一种多线程多GPU并行加速的最小二乘逆时偏移算法, 在GPU加速的最小二乘逆时偏移算法的基础上, 利用CPU的多核架构, 创建多线程协同操作, 调度多GPU进行并行加速计算和迭代参数的更新, 降低数据传输延迟, 大幅提升了计算效率。本文对炮集数据分块切割, 在GPU端实行粗粒度并行计算, 速度提升接近线性。Marmousi2截断模型和Marmousi模型的测试验证了该方法的有效性。
1 时空域最小二乘逆时偏移原理常密度声波方程表示如下:
| $ \frac{1}{{{v^2}\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}p\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}p\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = f\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},t} \right) $ | (1) |
式中:v(x)为速度场; p(x, t, xs)为波场; f(xs, t)为震源函数; t为时间; xs为震源位置。
假设速度场v(x)由背景速度场vb(x)和扰动速度场vs(x)叠加而成, 即v(x)=vb(x)+vs(x)。由波场叠加原理可知对应的波场也由背景波场pb(x, t, xs)和扰动波场ps(x, t, xs)叠加而成, 即p(x, t, xs)=pb(x, t, xs)+ps(x, t, xs)。令m(x)= [2vs(x)]/[vb(x)], 由波恩近似可得:
| $ \left\{ \begin{array}{l} \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \\ \;\;\;\;\;\;\;\;\;f\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},t} \right)\\ \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \\ \;\;\;\;\;\;\;\;\;m\left( \mathit{\boldsymbol{x}} \right)\frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}}\\ {d_{{\rm{cal}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{g}}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = {p_{\rm{s}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{g}}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \end{array} \right. $ | (2) |
式中:dcal(xg, t, xs)为模拟数据, xg为检波点位置。公式(2)的具体推导过程见附录A。为了简化表达, 公式(2)可采用矩阵向量形式表示为d=Lm, L为波恩正演算子, d为所有炮的模拟数据dcal(xg, t, xs)所构成的数据向量, m是m(x)的向量表达形式。
传统的偏移方法认为LT是波恩正演算子的伴随算子, 其表达式如下:
| $ \left\{ \begin{array}{l} \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = f\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},t} \right)\\ \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}{p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \\ \;\;\;\;\;\;\;\;\;{d_{{\rm{cal}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{g}}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)\\ m\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \sum\limits_t {\frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}}{p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \end{array} \right. $ | (3) |
式中:pr(x, t, xs)为根据模拟数据所得的检波点波场, 将所有炮的成像结果进行叠加计算, 即可获得最终的成像结果。为简化表达, 公式(3)也可表示为矩阵向量的形式:m=LTd。
根据最小二乘逆时偏移算法获得的扰动模型m来建立最小化目标函数:
| $ f\left( \mathit{\boldsymbol{m}} \right) = \frac{1}{2}{\left\| {\mathit{\boldsymbol{Lm}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|^2} $ |
式中:dobs为观测数据的矢量表示形式。
本文采用CLAERBOUT[31]提出的共轭梯度法, 由迭代算法获得扰动模型m的更新, 具体迭代流程如下:
| $ \mathit{\boldsymbol{r}} \leftarrow \mathit{\boldsymbol{LM}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $ | (4) |
| $ \begin{array}{l} {\rm{iterate}}\\ \left\{ {} \right.\\ \Delta \mathit{\boldsymbol{m}} \leftarrow {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{r}} \end{array} $ | (5) |
| $ \Delta \mathit{\boldsymbol{r}} \leftarrow \mathit{\boldsymbol{L}}\Delta \mathit{\boldsymbol{m}} $ | (6) |
| $ \begin{array}{l} \alpha = \left\{ {\left[ {\left( {{\mathit{\boldsymbol{S}}_k} \cdot {\mathit{\boldsymbol{S}}_k}} \right) \cdot \left( {\mathit{\boldsymbol{r}} \cdot \Delta \mathit{\boldsymbol{r}}} \right)} \right] - \left[ {\left( {{\mathit{\boldsymbol{S}}_k} \cdot \Delta \mathit{\boldsymbol{r}}} \right) \cdot } \right.} \right.\\ \;\;\;\;\left. {\left. {\left( {{\mathit{\boldsymbol{S}}_k} \cdot \Delta \mathit{\boldsymbol{r}}} \right)} \right]} \right\}/\left\{ {\left[ {\left( {\Delta \mathit{\boldsymbol{r}} \cdot \Delta \mathit{\boldsymbol{r}}} \right) \cdot \left( {{\mathit{\boldsymbol{S}}_k} \cdot {\mathit{\boldsymbol{S}}_k}} \right)} \right] - } \right.\\ \;\;\;\;\left. {\left[ {\left( {{\mathit{\boldsymbol{S}}_k} \cdot \Delta \mathit{\boldsymbol{r}}} \right) \cdot \left( {{\mathit{\boldsymbol{S}}_k} \cdot \Delta \mathit{\boldsymbol{r}}} \right)} \right]} \right\} \end{array} $ | (7) |
| $ \begin{array}{l} \beta = \left\{ {\left[ {\left( {{\mathit{\boldsymbol{S}}_k} \cdot \mathit{\boldsymbol{r}}} \right) \cdot \left( {\Delta \mathit{\boldsymbol{r}} \cdot \Delta \mathit{\boldsymbol{r}}} \right)} \right] - \left[ {\left( {\mathit{\boldsymbol{r}} \cdot \Delta \mathit{\boldsymbol{r}}} \right) \cdot } \right.} \right.\\ \;\;\;\left. {\left. {\left( {{\mathit{\boldsymbol{S}}_k} \cdot \Delta \mathit{\boldsymbol{r}}} \right)} \right]} \right\}/\left\{ {\left[ {\left( {\Delta \mathit{\boldsymbol{r}} \cdot \Delta \mathit{\boldsymbol{r}}} \right) \cdot \left( {{\mathit{\boldsymbol{S}}_k} \cdot {S_k}} \right)} \right] - } \right.\\ \;\;\;\left. {\left[ {\left( {{\mathit{\boldsymbol{S}}_k} \cdot \Delta \mathit{\boldsymbol{r}}} \right) \cdot \left( {{\mathit{\boldsymbol{S}}_k} \cdot \Delta \mathit{\boldsymbol{r}}} \right)} \right]} \right\} \end{array} $ | (8) |
| $ {\mathit{\boldsymbol{s}}_{k + 1}} = \alpha \Delta \mathit{\boldsymbol{m}} + \beta {\mathit{\boldsymbol{s}}_k} $ | (9) |
| $ {\mathit{\boldsymbol{S}}_{k + 1}} = \alpha \Delta \mathit{\boldsymbol{r}} + \beta {\mathit{\boldsymbol{S}}_k} $ | (10) |
| $ \mathit{\boldsymbol{m}} \leftarrow \Delta \mathit{\boldsymbol{m}} + {\mathit{\boldsymbol{s}}_{k + 1}} $ | (11) |
| $ \begin{array}{l} \mathit{\boldsymbol{r}} \leftarrow \Delta \mathit{\boldsymbol{r}} - {\mathit{\boldsymbol{S}}_{k + 1}}\\ \} \end{array} $ | (12) |
式中:r为数据残差; m为扰动模型, 也是迭代计算需要求取的目标解; Δm为模型域梯度; Δr为数据域共轭梯度; sk, Sk分别为第k次迭代中模型域和数据域的更新量; α, β分别为修正sk, Sk的步长参数, 〈·〉代表向量的点积运算。为提高计算效率, 在迭代过程中, 本文参照以下公式对梯度进行归一化补偿照明[32]:
| $ \Delta \mathit{\boldsymbol{m}} = \frac{{\Delta \mathit{\boldsymbol{m}}}}{{\sqrt {\sum\nolimits_{s = 1}^S {\int_0^{{t_{\max }}} {p_{\rm{b}}^2{\rm{d}}t} + {\gamma ^2}} } }} $ | (13) |
式中:γ为稳定性系数; s为当前炮数; S为总炮数; tmax为时间方向最大采样点数。
2 多GPU加速优化方法GPU加速技术能够有效提高并行算法的计算效率[33], 已在地震数据处理领域取得了较为广泛的应用[25-30]。关于GPU加速技术的基本流程不再赘述。本文采用基于CUDA平台的多GPU加速技术, 将GPU存储器优化和多线程多GPU并行加速方法应用于最小二乘逆时偏移算法, 利用GPU内部共享存储器和寄存器等高速存储器, 降低数据访问延迟, 提高计算效率, 结合CPU的多核架构, 分配多CPU线程协同调用多GPU进行加速计算。
2.1 存储器优化GPU包含多种存储器, 相较于全局存储器, 共享存储器和寄存器具有更高的数据传输带宽, 即更高的读写效率。因此, 采用共享存储器和寄存器作为数据存储器协助计算, 可获得更高的计算效率。
GPU加速计算时, 将数据的网格点划分为若干个Block(线程块), 各个Block的GPU线程可执行一对一的网格点数值计算。共享存储器是各个Block的内部存储器, 仅限于同一Block内的GPU线程访问, 寄存器则为每个GPU线程的私有存储器, 因此, 需合理分配存储器, 才能实现数据访问时的提速。
以(1)式的离散表达式为例(不考虑震源项):
| $ \begin{array}{l} P_{m,n}^{k + 1} = 2P_{m,n}^k - P_{m,n}^{k - 1} + \frac{{{v^2}\Delta {t^2}}}{{\Delta {x^2}}}\sum\nolimits_{l = - L}^L {C_l^LP_{m + l,n}^k} + \\ \;\;\;\;\;\;\;\;\;\;\frac{{{v^2}\Delta {t^2}}}{{\Delta {z^2}}}\sum\nolimits_{l = - L}^L {C_l^LP_{m,n + l}^k} \end{array} $ | (14) |
式中:上标k代表时间; 下标m, n分别代表水平方向和垂直方向; Pm, nk为离散形式下k时刻, 水平方向第m个网格点, 垂直方向第n个网格点的压力场; Δt表示时间步长; Δx, Δz分别为水平方向和垂直方向网格间距; v代表速度场; ClL为有限差分系数; L为空间差分阶数的一半。
假设压力场P水平方向网格点数为Nx, 垂直方向网格点数为Nz, 设某一个Block水平方向网格点数为Dx, 垂直方向网格点数为Dz。将压力场分割为若干个Block, 如图 1左半部分中黑色实线框所示, 每个Block中分配的共享存储器的网格尺寸为(Dx+2L)(Dz+2L)。由于不同Block中的共享存储器无法进行数据互访, 因此, 为满足(14)式中空间差分的计算条件, 根据图 1右半部分的方法分配Block中的共享存储器, 用于保存当前时刻的波场值Pm, nk。该共享存储器除保存Block中各GPU线程所对应坐标处的Pm, nk值(图 1中蓝色方块), 还根据差分阶数向外延伸并保存水平和垂直方向厚度为2L层的Pm, nk值。
|
图 1 Block划分及共享存储器分配示意 |
将全局存储器中对应坐标处的Pm, nk数据分别传输至各Block的共享存储器, 全局存储器中各个网格点处的速度数据、差分系数、时间步长和空间步长等参数传输至每个GPU线程的寄存器, 利用共享存储器和寄存器中的数据, 即可完成(14)式波动方程的离散计算。附录B为CUDA-C语言编程实现的GPU端存储器优化伪代码(基于空间8阶差分)。
2.2 多线程多GPU目前主流的计算设备中CPU端具备多核架构, 可同时触发多线程作业。同一计算设备中配备多个GPU设备即可支持多GPU并行运算。本文提出了多线程多GPU并行加速最小二乘逆时偏移算法, 根据CPU端多线程机制, 使每个CPU线程负责一个GPU的作业管理和数据传输, 将计算任务和数据分块, 传输至GPU端, 并以作业发送的方式, 实现多GPU并行运算, 具体情况如图 2所示。
|
图 2 多线程多GPU分配示意 |
最小二乘逆时偏移算法需进行迭代计算, 计算量随迭代次数线性增加, 本文采取的多GPU的并行策略为迭代计算时, 先对炮集数据进行分块, 再分派给各个GPU, 彼此独立地进行波场模拟计算。该策略既避免了多GPU间波场数据的实时交换, 又降低了频繁的数据传输引起的计算等待, 使得多GPU并行计算结果逼近线性加速效果。
共轭梯度法迭代时, 根据(7)式和(8)式, 对数据域维度(炮数×对应炮的道数×时间采样点数)的多个向量进行线性运算可求取参数α和β。编程实现时, 如采用BLAS库运行向量的线性运算, 需在GPU端的显存和CPU端的内存之间频繁地进行数据传输, 并且因等待数据传输而降低算法执行的效率。
向量点积运算满足(15)式~(17)式所示的性质, 式中A和B分别为两个向量, 对应的向量元素分别为a1, a2, …, an, an+1, an+2, …, a2n和b1, b2, …, bn, bn+1, bn+2, …, b2n, 该性质有利于多GPU并行计算的实现, 因此可以先分组计算, 再对各组结果求和获得最终结果。本文采用CUDA提供的线性代数程序库CUBLAS进行向量运算, 大部分运算均在GPU端执行, 减少了CPU端内存与GPU端显存之间的数据传输频率。
| $ \mathit{\boldsymbol{A}} = \left[ {{a_1},{a_2}, \cdots ,{a_n},{a_{n + 1}},{a_{n + 2}}, \cdots ,{a_{2n}}} \right] $ | (15) |
| $ \mathit{\boldsymbol{B}} = \left[ {{b_1},{b_2}, \cdots ,{b_n},{b_{n + 1}},{b_{n + 2}}, \cdots ,{b_{2n}}} \right] $ | (16) |
| $ \begin{array}{*{20}{c}} {\left\langle {\mathit{\boldsymbol{A}} \cdot \mathit{\boldsymbol{B}}} \right\rangle = \sum\nolimits_{i = 1}^{2n} {{a_i} \cdot {b_i}} = \sum\nolimits_{j = 1}^n {{a_j} \cdot {b_j}} + }\\ {\sum\nolimits_{k = n + 1}^{2n} {{a_k} \cdot {b_k}} } \end{array} $ | (17) |
本文采用Pthread接口进行CPU端的线程开发, 主线程负责数据同步和状态监控, 利用主线程调用函数“Pthread_create()”并创建多个并行子线程后, 各子线程负责对应数据块的读取, 随后调用对应的GPU进行计算。具体实现流程如图 3所示。
|
图 3 多线程多GPU最小二乘逆时偏移算法流程 |
主要步骤如下:
第1步, 主线程根据线程个数进行数据统计和分块;
第2步, 启动多线程, 各线程根据任务分配情况读取数据, 并将数据传输至GPU端;
第3步, 各线程调用GPU, 根据(5)式计算梯度;
第4步, 设置多线程阻塞函数“Pthread_join()”, 待所有线程执行完毕后调用CUBLAS库, 对各GPU端的梯度求和并同步;
第5步, 各线程调用GPU, 根据(6)式计算共轭梯度, 调用CUBLAS库在各GPU端完成共轭梯度法中所需的向量点积计算;
第6步, 设置多线程阻塞函数“Pthread join()”, 待所有线程执行完毕后, 对各GPU端计算向量点积结果, 并对各GPU所得点积结果进行求和同步, 然后根据(7)式和(8)式计算参数α和β;
第7步, 各线程调用GPU, 并根据(9)式~(12)式更新迭代结果及数据残差;
第8步, 判断是否满足迭代终止条件, 如果满足迭代终止条件, 则结束多线程, 输出数据; 否则重新迭代。
附录C为C语言编程实现的CPU端多线程作业触发的伪代码。
3 模型测试 3.1 最小二乘逆时偏移模型测试 3.1.1 模型测试1本文利用Marmousi2截断模型进行测试, 参数如下:纵、横向网格点数分别为296和600, 网格间距为15 m, 雷克子波主频为16 Hz, 时间采样间隔为1.0 ms, 采样点数为6 000, 设计20炮震源地表激发, 激发点均匀分布于水平方向1 485~7 185 m, 炮间距300 m, 每炮由199个检波器接收, 检波器均匀对称地分布于激发点两侧, 最大偏移距为1 485 m, 检波器间距15 m。
图 4a为准确速度模型, 即正演模型, 图 4b为背景速度模型, 即偏移算法采用的模型, 也是利用准确速度模型平滑所得的模型, 图 4c为扰动模型, 可利用准确速度模型和背景速度模型计算得到, 该扰动模型可视为最小二乘逆时偏移的理论解。图 5a为常规逆时偏移的成像结果; 图 5b为拉普拉斯去噪后的常规逆时偏移成像结果; 图 5c为最小二乘逆时偏移(50次迭代后)成像结果。对比图 5a, 图 5b和图 5c可以看出, 相较于常规逆时偏移, 最小二乘逆时偏移能够获得分辨率更高的成像结果, 能量更收敛, 照明范围也明显更广阔, 且振幅与理论值相近, 所得结果具有较为明确的物理意义。
|
图 4 模型参数 a 准确速度模型; b 背景速度模型; c 扰动模型 |
|
图 5 逆时偏移处理后得到的成像结果 a 常规逆时偏移; b 拉普拉斯去噪后的常规逆时偏移; c 最小二乘逆时偏移(50次迭代后) |
抽取水平方向3 km处不同迭代次数的最小二乘逆时偏移单道数据进行对比, 结果如图 6所示, 随着迭代次数的增加, 振幅和相位匹配逐渐逼近理论值。
|
图 6 水平方向3 km处单道不同迭代次数的最小二乘逆时偏移结果对比 a 迭代1次; b 迭代10次; c 迭代50次 |
图 7所示为数据残差的下降曲线, 可以看出, 本文方法可对残差数据进行有效更新, 数据残差随迭代次数增加而降低。
|
图 7 数据残差下降曲线 |
为进一步验证本文方法的适用性, 我们基于Marmousi模型进行测试, 参数如下:纵、横向网格点数分别为384和122, 网格间距为15 m, 雷克子波主频为16 Hz, 时间采样间隔为1.0 ms, 采样点数为3 001, 设计20炮震源地表激发, 激发点均匀分布于水平方向1 680~4 080 m, 炮间距120 m, 每炮由199个检波器接收, 检波器均匀对称地分布于激发点两侧, 最大偏移距为1 485 m, 检波器间距15 m。
图 8a为准确速度模型, 即正演模型, 图 8b为背景速度模型, 即偏移算法采用的模型, 也是准确速度模型平滑所得的模型, 图 8c为扰动模型m(x)。图 9a, 图 9b和图 9c分别为最小二乘逆时偏移1次、10次、50次迭代后的结果, 可以看出, 随着迭代次数的增加, 成像效果显著提升。
|
图 8 模型参数 a 准确速度模型; b 背景速度模型; c 扰动模型 |
|
图 9 不同迭代次数的最小二乘逆时偏移成像结果 a 迭代1次; b 迭代10次; c 迭代50次 |
如图 10和图 11所示, 将水平方向2.88 km处(图 8, 图 9中白色虚线处)和深度方向0.45 km处(图 8, 图 9中红色虚线处)理论数据和不同迭代次数的最小二乘逆时偏移结果进行对比, 可以看出, 随着迭代次数的增加, 本文方法所得结果的振幅和相位匹配逐渐逼近理论值, 也证明了本文方法对不同参数模型具普适性。
|
图 10 深度方向0.45 km处单道迭代次数分别为1次(a), 10次(b)和50次(c)的最小二乘逆时偏移结果对比 |
|
图 11 水平方向2.88 km处单道迭代次数分别为1次(a), 10次(b)和50次(c)次的最小二乘逆时偏移结果对比 |
本文分别采用常规GPU加速方法(调用全局存储器)和存储器优化方法(调用共享存储器和寄存器)对模型测试1中的数据进行1次迭代的最小二乘逆时偏移测试, 耗时情况如图 12所示, 常规GPU加速方法耗时约73.7 s, 存储器优化方法耗时约61.2 s。将本文提出的存储器优化方法应用于最小二乘逆时偏移算法的波场模拟后, 计算效率约提升17%。
|
图 12 常规GPU加速方法和存储器优化方法耗时对比 |
分别使用1~4个GPU进行1次迭代最小二乘逆时偏移计算。表 1所示为使用不同个数GPU时, 分配给各GPU的炮集数, 以1个GPU计算耗时为参考值, 计算各个GPU理论上的耗时情况, 当GPU个数为3时, 分配给GPU2设备的炮集数为7, 则理论上, GPU2的耗时应该为GPU2计算20个炮集数据所用时间的7/20倍, 以此类推。模型测试1中, 参与计算的GPU个数与耗时情况如图 13所示, 黑色实线为实际耗时情况, 可看出计算耗时随GPU个数增加而降低, 多GPU加速效果较为明显; 红色点划线为根据表 1中的分配情况, 多GPU并行算法达到理想化完全线性加速时理论上的耗时情况; 蓝色虚线所示为各个GPU计算每个炮集数据的平均耗时。从图 13中可看出, 本文方法的执行效率较高, 耗时情况接近线性加速, 但随着GPU个数的增加, 算法的实际执行效率有所降低。这是由于多GPU并行计算时, 需等待最慢的线程完成计算任务, 还需要执行梯度、更新步长等参数的同步计算, 因此无法完全实现线性加速, 从而导致随着GPU个数的增加, 平均每个炮集数据的计算耗时略有增加。
|
表 1 多GPU任务分配情况 |
|
图 13 多GPU实际、理论和平均计算耗时对比(基于模型测试1) |
模型测试2的模型尺寸和数据量均小于模型测试1的数据规模, 参与计算的GPU个数与耗时情况如图 14所示, 我们发现不同个数的GPU参与计算时, 实际耗时与理论耗时之间的差距会随着GPU个数的增加而增大, 模型测试2中二者差距的增幅大于模型测试1中的对应参数; 另外, 模型测试2中, 平均耗时随GPU个数的增加而提升的幅度大于模型测试1中的情况。由此说明, 随着模型尺寸和计算数据量的增加, 多线程多GPU并行计算方法的加速效率逐渐逼近线性加速, 主要原因在于随着计算量的提升, 并行算法的加速作用在整个计算过程中相对增加了, 而线程等待和数据同步等降低并行效率的计算相对减少了。
|
图 14 多GPU实际、理论和平均计算耗时对比(基于模型测试2) |
将GPU加速技术应用于二维时域最小二乘逆时偏移算法, 在CPU端采用多线程模式, 调用GPU, 实现了多GPU的并行加速计算, 迭代计算时, 调用了CUBLAS库函数协助计算, 在GPU端实现了数据的更新和同步计算, 减少了CPU与GPU间的数据传输, 降低了延迟, 明显提升了多线程多GPU最小二乘逆时偏移算法的计算效率。本文还采用了访问速度更快的共享存储器和寄存器, 进一步提升了GPU算法的执行效率。
本文提出的多线程多GPU并行加速最小二乘逆时偏移的思路也可应用于三维数据中, 但由于三维空间数据量大, 在进行多GPU加速计算时, 建议采用多GPU并行加速策略, 即对三维空间数据进行分块切割, 然后交由各GPU并行运算。多线程CPU支持数据的并行传输, 可继续充分发挥并行加速方法的优势。此外, 多线程多GPU并行加速方法可普遍应用于地震资料的迭代类成像算法以及全波形反演算法, 下一步的研究重点应集中于本文方法在实际地震资料处理中的应用。
附录A 波恩近似推导常密度波动方程:
| $ \frac{1}{{{v^2}\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}p\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}p\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = f\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},t} \right) $ | (A1) |
式中:x为空间坐标; xs为震源位置; t为时间; p为压力场; f为震源函数。
假设速度场v(x)由背景速度场vb(x)和扰动速度场vs(x)叠加而成:
| $ v\left( \mathit{\boldsymbol{x}} \right) = {v_{\rm{b}}}\left( \mathit{\boldsymbol{x}} \right) + {v_{\rm{s}}}\left( \mathit{\boldsymbol{x}} \right) $ | (A2) |
由波场叠加原理可知对应的波场由背景波场和扰动波场叠加构成:
| $ p\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = {p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) + {p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $ | (A3) |
则方程(A1)可转换为如下形式:
| $ \frac{1}{{{{\left[ {{v_{\rm{b}}}\left( \mathit{\boldsymbol{x}} \right) + {v_{\rm{s}}}\left( \mathit{\boldsymbol{x}} \right)} \right]}^2}}}\frac{{{\partial ^2}\left[ {{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) + {p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right]}}{{\partial {t^2}}} - {\nabla ^2}\left[ {{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) + {p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right] = f\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},t} \right) $ | (A4) |
对速度项进行二阶泰勒级数展开可得:
| $ \frac{1}{{{{\left[ {{v_{\rm{b}}}\left( \mathit{\boldsymbol{x}} \right) + {v_{\rm{s}}}\left( \mathit{\boldsymbol{x}} \right)} \right]}^2}}} \approx \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}} - \frac{{2{v_{\rm{s}}}\left( \mathit{\boldsymbol{x}} \right)}}{{{v_{\rm{b}}}\left( \mathit{\boldsymbol{x}} \right)}} $ |
将其代入式(A4)可得:
| $ \begin{array}{l} \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} + \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - \frac{{2{v_{\rm{s}}}\left( \mathit{\boldsymbol{x}} \right)}}{{v_{\rm{b}}^3\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - \frac{{2{v_{\rm{s}}}\left( \mathit{\boldsymbol{x}} \right)}}{{v_{\rm{b}}^3\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - \\ \;\;\;\;\;\;\;\;\;\;\;\;{\nabla ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) - {\nabla ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = f\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},t} \right) \end{array} $ | (A5) |
忽略高阶扰动波场项
| $ \left\{ \begin{array}{l} \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = f\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},t} \right)\\ \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \frac{{2{v_{\rm{s}}}\left( \mathit{\boldsymbol{x}} \right)}}{{v_{\rm{b}}^3\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} \end{array} \right. $ | (A6) |
令m(x)=
| $ \left\{ \begin{array}{l} \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = f\left( {{\mathit{\boldsymbol{x}}_{\rm{s}}},t} \right)\\ \frac{1}{{v_{\rm{b}}^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} - {\nabla ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = m\left( x \right)\frac{1}{{v_{\rm{b}}^3\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{p_{\rm{b}}}\left( {\mathit{\boldsymbol{x}},t,{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} \end{array} \right. $ |
#define BDIMX 16 //定义Block水平方向网格尺寸
#define BDIMY 16 //定义Block垂直方向网格尺寸
__global__ void step_forward_shared1(float *p0_d, float *p1_d, int nx, int ny, float dx, float dy, float dt, float *vb_d, float *coe_d) /*定义波动方程离散函数*/
{
/*定义全局存储器数据水平方向坐标*/
int tx = blockIdx.x*blockDim.x+threadIdx.x;
/*定义全局存储器数据垂直方向坐标*/
int ty = blockIdx.y*blockDim.y+threadIdx.y;
/*定义共享存储器数据水平方向坐标*/
int ix = threadIdx.x+4;
/*定义共享存储器数据垂直方向坐标*/
int iy = threadIdx.y+4;
__shared__ float s_data[BDIMY+2*4][BDIMX+2*4]; /*定义共享存储器*/
/*将数据由全局存储器传输至共享存储器*/
s_data[iy][ix]=p1_d[ty*nx+tx];
if(threadIdx.y＜4)
{
s_data[threadIdx.y][ix]=p1_d[(ty-4)*nx+tx];
s_data[threadIdx.y+BDIMY+4][ix]=p1_d[(ty+BDIMY)*nx+tx];
}
if(threadIdx.x＜4)
{
s_data[iy][threadIdx.x]=p1_d[ty*nx+tx-4];
s_data[iy][threadIdx.x+BDIMX+4]=p1_d[ty*nx+tx+BDIMX];
}
__syncthreads();
/*在寄存器中定义差分系数*/
float coe0, coe1, coe2, coe3, coe4;
/*将差分系数数据由全局存储器传输至寄存器*/
coe0=coe_d[0]; coe1=coe_d[1]; coe2=coe_d[2]; coe3=coe_d[3]; coe4=coe_d[4];
float v; /*在寄存器中定义速度数据*/
float tempx=0, tempz=0;
if(tx＞4&&tx＜nx-4&&ty＞4&&ty＜ny-4)
{
v=vb_d[ty*nx+tx]; /*将速度数据数据由全局存储器传输至寄存器*/
/*波动方程差分计算*/
tempx=coe0*s_data[iy][ix];
tempz=coe0*s_data[iy][ix];
tempx=tempx+coe1*(s_data[iy][ix+1]+s_data[iy][ix-1])+
coe2*(s_data[iy][ix+2]+s_data[iy][ix-2])+
coe3*(s_data[iy][ix+3]+s_data[iy][ix-3])+
coe4*(s_data[iy][ix+4]+s_data[iy][ix-4]);
tempz=tempz+coe1*(s_data[iy+1][ix]+s_data[iy-1][ix])+
coe2*(s_data[iy+2][ix]+s_data[iy-2][ix])+
coe3*(s_data[iy+3][ix]+s_data[iy-3][ix])+
coe4*(s_data[iy+4][ix]+s_data[iy-4][ix]);
p0_d[ty*nx+tx]=2.0 * s_data[iy][ix]
-p0_d[ty*nx+tx]
+v*v*dt*dt * tempx/dx/dx
+v*v*dt*dt * tempz/dy/dy;
}
}
附录C CPU端多线程作业触发伪代码void creat_thread(char *string, int nGPU, int GPUbeg, pthread_t *thread, void* f(void *), struct TGPUplan *arg) /*定义线程触发函数*/
{
int thread_status=0;
for(int iGPU=0;iGPU＜nGPU; iGPU++)
{
if((thread_status = pthread_create( & thread[iGPU], NULL, f, arg+iGPU))!= 0)
{printf("%sThread%dCreat Failed！\n", string, iGPU); }
else
{printf("%sThread%dCreat Successful！\n", string, iGPU); }
}
}
void joint_thread(char *string, int nGPU, pthread_t *thread) /*定义线程阻塞函数*/
{
for(int iGPU=0;iGPU＜nGPU; iGPU++)
{
if(thread[iGPU]!= 0)
{
pthread_join(thread[iGPU], NULL);
printf("%sThread %dRuns Over!\n", string, iGPU);
}
}
}
void* Function_name(void* arg) /*定义线程功能函数*/
{
TGPUplan*plan =(TGPUplan *) arg; /*传递线程功能函数的参数结构体*/
……
}
main()
{
……
int thread_status=0;
pthread_t thread[nGPU]; /*传递多线程*/
memset( & thread, 0, sizeof(thread)); /*初始化多线程*/
TGPUplan plan[nGPU]; /*创建多线程功能函数结构体*/
……
/*触发多线程函数*/
creat_thread("Function information", nGPU, GPUbeg, thread, Function_name, plan);
/*阻塞等待多线程函数*/
joint_thread("Function information", nGPU, thread);
……
}
| [1] |
SCHNEIDER W A. Integral Formulation for migration in two and three dimensions[J]. Geophysics, 1978, 43(1): 49-76. |
| [2] |
CLAERBOUT J. Toward a unified theory of reflector mapping[J]. Geophysics, 1971, 36(3): 467-481. DOI:10.1190/1.1440185 |
| [3] |
GAZDAG J. Wave equation migration with the phase-shift method[J]. Geophysics, 1978, 43(7): 1342-1351. DOI:10.1190/1.1440899 |
| [4] |
HEMON C. Equations d'onde et modeles[J]. Geophysical Prospecting, 1978, 26(4): 790-821. DOI:10.1111/gpr.1978.26.issue-4 |
| [5] |
BAYSAL E, KOSLOFF D D, SHERWOOD J W C. Reverse time migration[J]. Geophysics, 1983, 48(11): 1514-1524. DOI:10.1190/1.1441434 |
| [6] |
MCMECHAN G A. Migration by extrapolation of time-dependent boundary values[J]. Geophysical Prospecting, 1983, 31(3): 413-420. DOI:10.1111/gpr.1983.31.issue-3 |
| [7] |
TARANTOLA A. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 1984, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
| [8] |
LAMBARÉ G, VIRIEUX J, MADARIAGA R, et al. Iterative asymptotic inversion in the acoustic approximation[J]. Geophysics, 1992, 57(9): 1138-1154. DOI:10.1190/1.1443328 |
| [9] |
NEMETH T, WU C, SCHUSTER G T. Least-squares migration of incomplete reflection data[J]. Geophysics, 1999, 64(1): 208-221. DOI:10.1190/1.1444517 |
| [10] |
郭书娟, 马方正, 段心标, 等. 最小二乘逆时偏移成像方法的实现与应用研究[J]. 石油物探, 2015, 54(3): 301-308. GUO S J, MA F Z, DUAN X B, et al. Research of least-squares reverse-time migration imaging method and its application[J]. Geophysical Prospecting for Petroleum, 2015, 54(3): 301-308. DOI:10.3969/j.issn.1000-1441.2015.03.008 |
| [11] |
YAO G, JAKUBOWICZ H. Least-Squares Reverse-Time Migration[J]. Expanded Abstracts of 82nd Annual Internat SEG Mtg, 2012, 3406-3410. |
| [12] |
郭振波, 李振春. 最小平方逆时偏移真振幅成像[J]. 石油地球物理勘探, 2014, 49(1): 113-120. GUO Z B, LI Z C. True-amplitude imaging based on least-squares reverse time migration[J]. Oil Geophysical Prospecting, 2014, 49(1): 113-120. |
| [13] |
ZHANG Y, DUAN L, XIE Y. A stable and practical implementation of least-squares reverse time migration[J]. Geophysics, 2015, 80(1): V23-V31. |
| [14] |
ZHANG Q, ZHOU H, CHEN H, et al. Least-squares reverse time migration with and without source wavelet estimation[J]. Journal of Applied Geophysics, 2016, 134(1): 1-10. |
| [15] |
李庆洋, 黄建平, 李振春, 等. 去均值归一化互相关最小二乘逆时偏移及其应用[J]. 地球物理学报, 2016, 59(8): 3006-3015. LI Q Y, HUANG J P, LI Z C, et al. Mean-residual normalized cross-correlation least-squares reverse time migration and its application[J]. Chinese Journal of Geophysics, 2016, 59(8): 3006-3015. |
| [16] |
刘学建, 刘伊克. 表面多次波最小二乘逆时偏移成像[J]. 地球物理学报, 2016, 59(9): 3354-3365. LIU X J, LIU Y K. Least-squares reverse-time migration of surface-related multiples[J]. Chinese Journal of Geophysics, 2016, 59(9): 3354-3365. |
| [17] |
BERKHOUT A J. Areal shot record technology[J]. Journal of Seismic Exploration, 1992, 1(3): 251-264. |
| [18] |
黄建平, 李闯, 李庆洋, 等. 一种基于平面波静态编码的最小二乘逆时偏移方法[J]. 地球物理学报, 2015, 58(6): 2046-2056. HUANG J P, LI C, LI Q Y, et al. Least-squares reverse time migration with static plane-wave encoding[J]. Chinese Journal of Geophysics, 2015, 58(6): 2046-2056. |
| [19] |
LI C, HUANG J P, LI Z C, et al. Preconditioned prestack plane-wave least squares reverse time migration with singular spectrum constraint[J]. Applied Geophysics, 2017, 14(1): 73-86. |
| [20] |
DAI W, HUANG Y, SCHUSTER G T. Least-squares reverse time migration of marine data with frequency-selection encoding[J]. Expanded Abstracts of 83rd Annual Internat SEG Mtg, 2013, 3231-3236. |
| [21] |
李闯, 黄建平, 李振春, 等. 平面波最小二乘逆时偏移编码策略分析[J]. 石油物探, 2015, 54(5): 592-601. LI C, HUANG J P, LI Z C, et al. Analysis on encoding strategies of plane-wave least-square reverse time migration[J]. Geophysical Prospecting for Petroleum, 2015, 54(5): 592-601. DOI:10.3969/j.issn.1000-1441.2015.05.012 |
| [22] |
LIU Y S, TENG J W, XU T, et al. An efficient step-length formula for correlative least-squares reverse time migration[J]. Geophysics, 2016, 81(4): S221-S238. DOI:10.1190/geo2015-0529.1 |
| [23] |
李闯, 黄建平, 李振春, 等. 预条件最小二乘逆时偏移方法[J]. 石油地球物理勘探, 2016, 51(3): 513-520. LI C, HUANG J P, LI Z C, et al. Preconditioned least-squares reverse time migration[J]. Oil Geophysical Prospecting, 2016, 51(3): 513-520. |
| [24] |
WU D, YAO G, CAO J, et al. Least-squares RTM with L1 norm regularization[J]. Journal of Geophysics and Engineering, 2016, 13(5): 666-673. DOI:10.1088/1742-2132/13/5/666 |
| [25] |
李博, 刘红伟, 刘国峰, 等. 地震叠前逆时偏移算法的CPU/GPU实施对策[J]. 地球物理学报, 2010, 53(12): 2938-2943. LI B, LIU H W, LIU G F, et al. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU[J]. Chinese Journal of Geophysics, 2010, 53(12): 2938-2943. DOI:10.3969/j.issn.0001-5733.2010.12.017 |
| [26] |
刘红伟, 李博, 刘洪, 等. 地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报, 2010, 53(7): 1725-1733. LIU H W, LI B, LIU H, et al. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J]. Chinese Journal of Geophysics, 2010, 53(7): 1725-1733. DOI:10.3969/j.issn.0001-5733.2010.07.024 |
| [27] |
SHI Y, WANG Y H. Reverse time migration of 3D vertical seismic profile data[J]. Geophysics, 2016, 81(1): S31-S38. |
| [28] |
石颖, 刘洪, 邹振. 基于波动方程表面多次波预测与自适应相减方法研究[J]. 地球物理学报, 2010, 53(7): 1716-1724. SHI Y, LIU H, ZOU Z. Surface-related multiples prediction based on wave equation and adaptive subtraction investigation[J]. Chinese Journal of Geophysics, 2010, 53(7): 1716-1724. DOI:10.3969/j.issn.0001-5733.2010.07.023 |
| [29] |
石颖, 王维红, 李莹, 等. 基于波动方程三维表面多次波预测方法研究[J]. 地球物理学报, 2013, 56(6): 2023-2032. SHI Y, WANG W H, LI Y, et al. 3D surface-related multiple prediction approach investigation based on wave equation[J]. Chinese Journal of Geophysics, 2013, 56(6): 2023-2032. |
| [30] |
郭雪豹, 刘洪, 石颖. 基于频域衰减的时域全波形反演[J]. 地球物理学报, 2016, 59(10): 3777-3787. GUO X B, LIU H, SHI Y. Time domain full waveform inversion based on frequency attenuation[J]. Chinese Journal of Geophysics, 2016, 59(10): 3777-3787. DOI:10.6038/cjg20161022 |
| [31] |
CLAERBOUT J. Earth soundings analysis:processing versus inversion[M]. Boston: Blackwell Scientific Publication, 2004: 131-153.
|
| [32] |
YANG P L, GAO J H, WANG B L. Agraphics processing unit implementation of time-domain full-waveform inversion[J]. Geophysics, 2015, 80(3): F31-F39. DOI:10.1190/geo2014-0283.1 |
| [33] |
MICIKEVICIUS P.3D finite difference computation on GPUs using CUDA[C]//Proceedings of 2nd Workshop on General Purpose Processing on Graphics Processing Units.New York: ACM, 2008: 79-84 https://ci.nii.ac.jp/naid/10030070005
|