基于多线程多GPU并行加速的最小二乘逆时偏移算法

柯 璇,石 颖,张 伟,张 振,何 伟

(1.东北石油大学地球科学学院,黑龙江大庆 163318;2.中国石油塔里木油田分公司勘探开发研究院,新疆库尔勒,841000;3.中石化石油工程地球物理公司南方分公司,四川成都,610041)

地震偏移成像旨在对地下构造的反射信号重新归位,并根据地震资料刻画地下构造[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 时空域最小二乘逆时偏移原理

常密度声波方程表示如下:

(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)],由波恩近似可得:

(2)

式中:dcal(xg,t,xs)为模拟数据,xg为检波点位置。公式(2)的具体推导过程见附录A。 为了简化表达,公式(2)可采用矩阵向量形式表示为d=Lm,L为波恩正演算子,d为所有炮的模拟数据dcal(xg,t,xs)所构成的数据向量,m是m(x)的向量表达形式。

传统的偏移方法认为LT是波恩正演算子的伴随算子,其表达式如下:

(3)

式中:pr(x,t,xs)为根据模拟数据所得的检波点波场,将所有炮的成像结果进行叠加计算,即可获得最终的成像结果。为简化表达,公式(3)也可表示为矩阵向量的形式:m=LTd。

根据最小二乘逆时偏移算法获得的扰动模型m来建立最小化目标函数:

式中:dobs为观测数据的矢量表示形式。

本文采用CLAERBOUT[31]提出的共轭梯度法,由迭代算法获得扰动模型m的更新,具体迭代流程如下:

式中:r为数据残差;m为扰动模型,也是迭代计算需要求取的目标解;Δm为模型域梯度;Δr为数据域共轭梯度;sk,Sk分别为第k次迭代中模型域和数据域的更新量;α,β分别为修正sk,Sk的步长参数,〈·〉代表向量的点积运算。为提高计算效率,在迭代过程中,本文参照以下公式对梯度进行归一化补偿照明[32]:

(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)式的离散表达式为例(不考虑震源项):

(14)

图1 Block划分及共享存储器分配示意

2.2 多线程多GPU

目前主流的计算设备中CPU端具备多核架构,可同时触发多线程作业。同一计算设备中配备多个GPU设备即可支持多GPU并行运算。本文提出了多线程多GPU并行加速最小二乘逆时偏移算法,根据CPU端多线程机制,使每个CPU线程负责一个GPU的作业管理和数据传输,将计算任务和数据分块,传输至GPU端,并以作业发送的方式,实现多GPU并行运算,具体情况如图2所示。

最小二乘逆时偏移算法需进行迭代计算,计算量随迭代次数线性增加,本文采取的多GPU的并行策略为迭代计算时,先对炮集数据进行分块,再分派给各个GPU,彼此独立地进行波场模拟计算。该策略既避免了多GPU间波场数据的实时交换,又降低了频繁的数据传输引起的计算等待,使得多GPU并行计算结果逼近线性加速效果。

共轭梯度法迭代时,根据(7)式和(8)式,对数据域维度(炮数×对应炮的道数×时间采样点数)的多个向量进行线性运算可求取参数α和β。编程实现时,如采用BLAS库运行向量的线性运算,需在GPU端的显存和CPU端的内存之间频繁地进行数据传输,并且因等待数据传输而降低算法执行的效率。

图2 多线程多GPU分配示意

向量点积运算满足(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端显存之间的数据传输频率。

本文采用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,网格间距为15m,雷克子波主频为16Hz,时间采样间隔为1.0ms,采样点数为6000,设计20炮震源地表激发,激发点均匀分布于水平方向1485～7185m,炮间距300m,每炮由199个检波器接收,检波器均匀对称地分布于激发点两侧,最大偏移距为1485m,检波器间距15m。

图4a为准确速度模型,即正演模型,图4b为背景速度模型,即偏移算法采用的模型,也是利用准确速度模型平滑所得的模型,图4c为扰动模型,可利用准确速度模型和背景速度模型计算得到,该扰动模型可视为最小二乘逆时偏移的理论解。图5a为常规逆时偏移的成像结果;图5b为拉普拉斯去噪后的常规逆时偏移成像结果;图5c为最小二乘逆时偏移(50次迭代后)成像结果。

图4 模型参数a 准确速度模型; b 背景速度模型; c 扰动模型

图5 逆时偏移处理后得到的成像结果a 常规逆时偏移; b 拉普拉斯去噪后的常规逆时偏移; c 最小二乘逆时偏移(50次迭代后)

对比图5a,图5b和图5c可以看出,相较于常规逆时偏移,最小二乘逆时偏移能够获得分辨率更高的成像结果,能量更收敛,照明范围也明显更广阔,且振幅与理论值相近,所得结果具有较为明确的物理意义。

抽取水平方向3km处不同迭代次数的最小二乘逆时偏移单道数据进行对比,结果如图6所示,随着迭代次数的增加,振幅和相位匹配逐渐逼近理论值。

图6 水平方向3km处单道不同迭代次数的最小二乘逆时偏移结果对比a 迭代1次; b 迭代10次; c 迭代50次

图7所示为数据残差的下降曲线,可以看出,本文方法可对残差数据进行有效更新,数据残差随迭代次数增加而降低。

图7 数据残差下降曲线

3.1.2 模型测试2

为进一步验证本文方法的适用性,我们基于Marmousi模型进行测试,参数如下:纵、横向网格点数分别为384和122,网格间距为15m,雷克子波主频为16Hz,时间采样间隔为1.0ms,采样点数为3001,设计20炮震源地表激发,激发点均匀分布于水平方向1680～4080m,炮间距120m,每炮由199个检波器接收,检波器均匀对称地分布于激发点两侧,最大偏移距为1485m,检波器间距15m。

图8a为准确速度模型,即正演模型,图8b为背景速度模型,即偏移算法采用的模型,也是准确速度模型平滑所得的模型,图8c为扰动模型m(x)。

图8 模型参数a 准确速度模型; b 背景速度模型; c 扰动模型

图9 不同迭代次数的最小二乘逆时偏移成像结果a 迭代1次; b 迭代10次; c 迭代50次

图9a,图9b和图9c分别为最小二乘逆时偏移1次、10次、50次迭代后的结果,可以看出,随着迭代次数的增加,成像效果显著提升。

如图10和图11所示,将水平方向2.88km处(图8,图9中白色虚线处)和深度方向0.45km处(图8,图9中红色虚线处)理论数据和不同迭代次数的最小二乘逆时偏移结果进行对比,可以看出,随着迭代次数的增加,本文方法所得结果的振幅和相位匹配逐渐逼近理论值,也证明了本文方法对不同参数模型具普适性。

图10 深度方向0.45km处单道迭代次数分别为1次(a),10次(b)和50次(c)的最小二乘逆时偏移结果对比

图11 水平方向2.88km处单道迭代次数分别为1次(a),10次(b)和50次(c)次的最小二乘逆时偏移结果对比

3.2 最小二乘逆时偏移的多GPU储存器优化方法测试

本文分别采用常规GPU加速方法(调用全局存储器)和存储器优化方法(调用共享存储器和寄存器)对模型测试1中的数据进行1次迭代的最小二乘逆时偏移测试,耗时情况如图12所示,常规GPU加速方法耗时约73.7s,存储器优化方法耗时约61.2s。将本文提出的存储器优化方法应用于最小二乘逆时偏移算法的波场模拟后,计算效率约提升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)

4 结论与认识

将GPU加速技术应用于二维时域最小二乘逆时偏移算法,在CPU端采用多线程模式,调用GPU,实现了多GPU的并行加速计算,迭代计算时,调用了CUBLAS库函数协助计算,在GPU端实现了数据的更新和同步计算,减少了CPU与GPU间的数据传输,降低了延迟,明显提升了多线程多GPU最小二乘逆时偏移算法的计算效率。本文还采用了访问速度更快的共享存储器和寄存器,进一步提升了GPU算法的执行效率。

本文提出的多线程多GPU并行加速最小二乘逆时偏移的思路也可应用于三维数据中,但由于三维空间数据量大,在进行多GPU加速计算时,建议采用多GPU并行加速策略,即对三维空间数据进行分块切割,然后交由各GPU并行运算。多线程CPU支持数据的并行传输,可继续充分发挥并行加速方法的优势。此外,多线程多GPU并行加速方法可普遍应用于地震资料的迭代类成像算法以及全波形反演算法,下一步的研究重点应集中于本文方法在实际地震资料处理中的应用。