張 全 林柏櫟 楊 勃 彭 博 張 偉 涂 然

(①西南石油大學(xué)計算機(jī)科學(xué)學(xué)院， 四川成都 610500； ②電子科技大學(xué)信息與通信工程學(xué)院， 四川成都 611731； ③國家電網(wǎng)重慶市電力公司電力科學(xué)研究院，重慶 404100)

0 引言

隨著油氣勘探的不斷深入，地震勘探也面臨越來越大的挑戰(zhàn)，勘探目標(biāo)逐漸從簡單構(gòu)造轉(zhuǎn)向更復(fù)雜的地質(zhì)構(gòu)造。在構(gòu)造復(fù)雜的低信噪比地區(qū)，其發(fā)育的噪聲干擾往往會弱化或掩蓋有效波信息，導(dǎo)致地震資料無法準(zhǔn)確反映地下真實(shí)地質(zhì)狀況[1-2]。地震資料中的噪聲主要是面波、多次波及隨機(jī)噪聲，其中多次波的存在不僅嚴(yán)重干擾一次波的成像質(zhì)量，降低成像結(jié)果信噪比，且在地震剖面上易產(chǎn)生構(gòu)造假象，掩蓋中深部有利構(gòu)造，影響地震資料的精細(xì)解釋，所以在偏移前需壓制和去除多次波[2-3]。拋物線Radon變換[4]是現(xiàn)今壓制多次波的常用方法。

Radon變換[5]由奧地利數(shù)學(xué)家拉東(Radon)率先提出，后被廣泛地應(yīng)用于醫(yī)學(xué)、地球物理、核磁共振、天文和分子生物等領(lǐng)域。1971年，Claerbout等[6]將其引入地震數(shù)據(jù)處理中，奠定了Radon變換在該領(lǐng)域的應(yīng)用基礎(chǔ)。Thorson等[7]將Radon建模為一個稀疏逆問題，在Radon域?qū)崿F(xiàn)了高分辨率數(shù)據(jù)處理。隨后Scales等[8]利用迭代加權(quán)最小二乘算法求解了稀疏逆問題。Sacchi等[9]提出的頻域稀疏Radon變換方法被廣泛地應(yīng)用于地震信號處理。Trad等[10]通過混合Radon變換成功地壓制了噪聲。劉喜武等[11]基于最小二乘反演法并導(dǎo)出稀疏約束條件下的共軛梯度算法，實(shí)現(xiàn)了高分辨率Radon變換。Schonewille等[12]證明了Radon變換在時域的稀疏性高于頻域，故認(rèn)為Radon變換在時域具有更高分辨率。Lu[13]提出一種基于迭代收縮的快速稀疏Radon變換方法，提高了Radon模型的稀疏性。

在實(shí)際地震數(shù)據(jù)處理中，Radon變換存在分辨率低和計算效率不高兩大缺陷[14]。Radon模型分辨率的提高可減少空間假頻和假象的產(chǎn)生而取得更佳處理效果，但伴隨更大計算量。況且，地震勘探數(shù)據(jù)量的日益龐大，使拋物線Radon變換在地震數(shù)據(jù)處理上耗費(fèi)更長時間。

隨著GPU-CPU異構(gòu)結(jié)構(gòu)的出現(xiàn)，高性能計算技術(shù)取得快速且巨大的發(fā)展，并逐漸應(yīng)用于地震勘探領(lǐng)域。張軍華等[15]全面調(diào)研了地球物理勘探領(lǐng)域?qū)Ω咝阅苡嬎愕男枨蠛蛻?yīng)用現(xiàn)狀。近年來，更有眾多學(xué)者將CPU-GPU異構(gòu)并行方法應(yīng)用于疊前數(shù)據(jù)處理[16-20]。本文基于CPU-GPU異構(gòu)平臺，對拋物線Radon變換算法做并行化分析和實(shí)現(xiàn)。

1 算法原理及串行流程

1.1 算法原理

多次波的壓制通常是在做動校正后進(jìn)行的。經(jīng)合適的動校正后，一次波同相軸水平； 而校正不足時多次波同相軸沿拋物線軌跡分布，則可沿拋物線做Radon變換，并在Radon域壓制多次波。本文采用Lu[13]提出的拋物線Radon變換法?；旌嫌驋佄锞€Radon變換的數(shù)學(xué)模型表示為

d′=F-1[LF(m)]

(1)

式中：m為Radon域數(shù)據(jù)；d′為地震數(shù)據(jù)；L為頻域Radon算子； F(·)為傅里葉變換； F-1(·)為傅里葉逆變換。為了使Radon域數(shù)據(jù)盡量稀疏且Radon反變換后的數(shù)據(jù)d′與地震數(shù)據(jù)的誤差盡量小，所以最小化如下目標(biāo)方程

(2)

式中λ為正則化系數(shù)。

結(jié)合迭代收縮閾值算法[21]，對目標(biāo)方程最優(yōu)解的求解變?yōu)榈嬎鉳k，即有

mk=Tα{mk-1+2tF-1[A(F(d)-LF(mk-1))]}

(3)

式中：k表示迭代次數(shù)；mk表示迭代第k次的Radon域數(shù)據(jù)；t為迭代步長；A=(LTL)-1LT為算子L的廣義逆。

迭代收縮閾值算法是在每一步迭代過程中采用

(4)

(5)

式中： ave表示二維均值濾波； max表示取最大值； |m|表示對向量m每個元素取絕對值。且有

(6)

式中α為標(biāo)量。

1.2 串行算法流程

根據(jù)Lu[13]提出的一種基于迭代收縮的快速稀疏Radon變換方法的算法原理，算法實(shí)現(xiàn)流程如圖1所示，核心計算部分分為如下四步。

圖1 拋物線Radon變換串行算法流程

(1)首先從硬盤讀取一個地震道集數(shù)據(jù)到內(nèi)存，通過傅里葉變換將地震數(shù)據(jù)從時域轉(zhuǎn)換到頻域； 然后計算頻域Radon算子L。由于此時L的計算只與地震道集的道數(shù)和采樣點(diǎn)數(shù)有關(guān)，所以當(dāng)前道集與前一道集的道數(shù)和采樣點(diǎn)數(shù)相同時，不需再計算L，能直接用前一道集的計算結(jié)果，即可轉(zhuǎn)到步驟(3)； 不相同時需重新計算L，順接到步驟(2)。

(2)在頻域內(nèi)構(gòu)造Radon變換算子L，并計算其廣義逆A=(LTL)-1LT。

(3)基于迭代收縮閾值算法更新Radon域模型，根據(jù)步驟(2)計算結(jié)果，得到初始Radon域模型m0=F-1[AF(d)]，且k=k+1，利用式(3)得到第k次迭代的Radon域模型mk。當(dāng)k≤K時，不斷更新mk。該步驟為本算法的計算熱點(diǎn)。

(4)設(shè)定一個曲率q，將曲率大于q的Radon域數(shù)據(jù)作為多次波估計結(jié)果，用原始地震數(shù)據(jù)減去多次波估計結(jié)果，就得到一次波。所有道集處理完畢則結(jié)束，否則轉(zhuǎn)到步驟(1)。

2 并行拋物線Radon變換算法

2.1 GPU的并行實(shí)現(xiàn)

首先在單個GPU上實(shí)現(xiàn)拋物線Radon變換算法的并行化，其流程如圖2所示。

由于CPU與GPU不能相互訪問各自存儲器，所以通過PCIE總線用拷貝方式實(shí)現(xiàn)數(shù)據(jù)在內(nèi)存與顯存間的傳輸。首先在CPU上從硬盤讀取地震數(shù)據(jù)到內(nèi)存，然后將地震數(shù)據(jù)和所需參數(shù)拷貝到GPU，在GPU端實(shí)施拋物線Radon變換，最后將計算結(jié)果拷貝到CPU，由CPU將結(jié)果從內(nèi)存寫回硬盤。在GPU上實(shí)施拋物線Radon變換分為四個步驟。

(1)由于算法在頻域內(nèi)實(shí)現(xiàn)，首先要將時域地震數(shù)據(jù)轉(zhuǎn)換到頻域，本文采用CUDA提供的快速傅里葉變換(cuFFT)實(shí)現(xiàn)。在后述步驟中，地震數(shù)據(jù)和Radon域數(shù)據(jù)在計算過程中仍需在頻域與時域間進(jìn)行傅里葉正、逆變換，一律使用CUDA提供的cuFFT庫進(jìn)行處理。cuFFT庫的使用一方面能給算法帶來性能提升，另一方面可降低程序設(shè)計復(fù)雜度，提高開發(fā)效率。值得注意的是，GPU上利用cuFFT庫計算傅里葉變換的結(jié)果與CPU上的計算結(jié)果存在一定差異，但這種差異是由于兩種不同的處理器架構(gòu)在機(jī)器指令上的差異導(dǎo)致的，該差異與真值的比值在10-6以下，通常可忽略。

圖2 拋物線Radon變換GPU并行算法流程

(2)首先構(gòu)造頻率域算子L，再計算矩陣LTL和(LTL)-1，最后計算A=(LTL)-1LT。通過調(diào)用CUDA提供的cuBLAS庫、cuSOLVER庫做線性代數(shù)運(yùn)算。

(3)基于迭代收縮閾值算法更新Radon域模型。首先使用cuFFT庫和線性代數(shù)庫計算Radon初始估計模型m0=F-1[AF(d)]； 再通過式(3)進(jìn)行收縮閾值操作，這時需對Radon域數(shù)據(jù)(nq×np，nq為曲率個數(shù)，np為采樣點(diǎn)數(shù))進(jìn)行均值濾波處理，均值濾波模板尺寸為n×m。在做均值濾波前，需先擴(kuò)充Radon域數(shù)據(jù)矩陣，擴(kuò)充后大小為(nq+2n)×(np+2m)。矩陣擴(kuò)充后進(jìn)行均值濾波，分別用兩個核函數(shù)實(shí)現(xiàn)，一個核函數(shù)按模板m在列方向上累加求平均值，另一個核函數(shù)按模板n在行方向上累加求平均值，最終得到均值濾波結(jié)果； 在迭代更新Radon域模型時計算A[F(d)-LF(mk-1)]，使用cuFFT庫計算傅里葉變換，使用線性代數(shù)庫計算矩陣乘法； 拋物線Radon變換算法中Radon域模型更新為迭代過程，在迭代開始前將所需數(shù)據(jù)從CPU一次拷貝到GPU，避免了迭代時CPU與GPU之間數(shù)據(jù)的傳輸，使得迭代過程可全速進(jìn)行。

(4)對步驟(3)所得Radon域數(shù)據(jù)m中曲率大于q的數(shù)據(jù)作為多次波估計結(jié)果F-1[LF(m)]，通過cuFFT庫和CUDA提供的線性代數(shù)庫進(jìn)行計算。最后通過d-F-1[LF(m)]得到有效波。

2.2 CPU-GPU異構(gòu)平臺的協(xié)同并行實(shí)現(xiàn)

對于上述GPU的并行算法，在執(zhí)行過程中需要主機(jī)的一個線程完成GPU核函數(shù)的調(diào)用，主機(jī)的其他線程都處于空閑狀態(tài)，勢必造成計算資源的浪費(fèi)。為了充分利用多核CPU和多GPU異構(gòu)計算平臺的計算資源，本文進(jìn)一步提出一種協(xié)同計算方案，實(shí)現(xiàn)拋物線Radon變換算法。

該方案中，對于支持m個CPU線程和n個GPU的異構(gòu)計算平臺，有具體的并行優(yōu)化流程(圖3)。

(1)主機(jī)端主線程獲取CPU支持的線程數(shù)m，并對N個地震道集進(jìn)行任務(wù)的劃分，創(chuàng)建一個有N個地震道集的任務(wù)隊(duì)列。

(2)每個子線程從任務(wù)隊(duì)列取一個任務(wù)執(zhí)行，即每個子線程負(fù)責(zé)一個地震道集的拋物線Radon變換的計算，m個子線程并行執(zhí)行m個任務(wù)。

(3)每個子線程領(lǐng)取一個地震道集任務(wù)后，要將該地震道集數(shù)據(jù)從硬盤讀取到內(nèi)存； 然后申請CPU計算資源，若子線程成功獲取CPU計算資源，則對拋物線Radon變換算法進(jìn)行預(yù)處理； 否則，子線程被阻塞，直到計算機(jī)操作系統(tǒng)將其調(diào)度。

(4)子線程申請GPU計算資源，若成功獲取GPU計算資源，就可調(diào)用GPU核函數(shù)完成并行拋物線Radon變換算法； 否則，調(diào)用CPU計算資源完成串行拋物線Radon變換算法。

(5)每個子線程在計算結(jié)束后將地震數(shù)據(jù)從內(nèi)存寫回硬盤。最后判斷任務(wù)隊(duì)列是否為空，當(dāng)任務(wù)隊(duì)列為空，則任務(wù)結(jié)束； 否則，每個線程繼續(xù)獲取任務(wù)，重復(fù)上面的操作。

本實(shí)驗(yàn)中使用的CPU-GPU異構(gòu)平臺由12個CPU核心和兩個GPU組成。該平臺支持CPU超線程技術(shù)，每個核心支持雙線程，所以該異構(gòu)計算平臺支持24個CPU線程。在該實(shí)驗(yàn)異構(gòu)平臺上任務(wù)分配與執(zhí)行情況如圖4所示：輸入N個地震道集，主線程創(chuàng)建23個子線程，然后進(jìn)行任務(wù)分配并參與計算，從而實(shí)現(xiàn)24個線程并行執(zhí)行不同道集任務(wù)。

圖3 CPU-GPU異構(gòu)協(xié)同并行流程

圖4 異構(gòu)平臺任務(wù)分配與執(zhí)行

由于該實(shí)驗(yàn)平臺有兩塊GPU(GPU0和GPU1)，則由兩個CPU線程分別調(diào)用GPU0和GPU1對地震道集進(jìn)行計算。因?yàn)镚PU計算速度比CPU快，所以在一個CPU線程執(zhí)行周期，GPU計算了多個地震道集。最終，在CPU-GPU異構(gòu)計算平臺上，計算結(jié)果分別由CPU和GPU的計算結(jié)果構(gòu)成。

3 數(shù)據(jù)應(yīng)用

3.1 模擬數(shù)據(jù)

本文使用SeismicLab工具包中合成的CMP多次波模擬數(shù)據(jù)。該動校正后的CMP數(shù)據(jù)(圖5a)包括49道，每道有1001個采樣點(diǎn)。對比SeismicLab工具包中拋物線Radon變換壓制多次波的結(jié)果(圖5b)與本文方法的去噪結(jié)果(圖5c)，可見后者對多次波有更好壓制效果。

圖5 動校正后CMP道集(a)及SeismicLab拋物線Radon變換法(b)與本文算法(c)壓制多次波效果對比

3.2 實(shí)際數(shù)據(jù)

針對CPU-GPU異構(gòu)平臺拋物線Radon變換算法并行優(yōu)化的性能進(jìn)行測試。所有測試均在表1的軟硬件環(huán)境下進(jìn)行，所使用的實(shí)驗(yàn)平臺有兩個CPU處理器12個核心(支持24線程)和兩塊NVIDIA Tesla K20c。

表1 軟硬件測試環(huán)境

實(shí)測數(shù)據(jù)取自珠江口盆地白云深水區(qū)。該區(qū)地形崎嶇，地層傾角較大，橫向巖性變化劇烈，地質(zhì)構(gòu)造復(fù)雜多變，多次波極為發(fā)育，是影響該區(qū)地震資料成像效果的主要原因[22]。實(shí)際測線范圍是： Inline2200～2300，Xline1500～1700，且每個地震道集有51道，采樣點(diǎn)數(shù)為376，采樣間隔為4ms。

先測試并考察算法精度，通過比較CPU串行算法與單GPU并行算法對多次波的壓制結(jié)果，驗(yàn)證GPU并行結(jié)果的正確性。由于CPU-GPU并行結(jié)果是由CPU結(jié)果和GPU結(jié)果組成的，所以只需驗(yàn)證GPU結(jié)果即驗(yàn)證了CPU-GPU結(jié)果的正確性。選取該區(qū)某一地震道集，從其動校正后的原始道集數(shù)據(jù)(圖6a)中可看到明顯存在的多次波，在經(jīng)CPU串行算法(圖6b)、單GPU并行算法(圖6c)的處理結(jié)果中，多次波已被壓制。由于乘加運(yùn)算在GPU上由一條指令完成，而CPU需多條指令，從而導(dǎo)致計算結(jié)果有微小差異，但這種差異在高性能計算領(lǐng)域?qū)僬，F(xiàn)象，可忽略不計[23]。

以CPU處理結(jié)果為參照，將CPU運(yùn)行結(jié)果減去GPU運(yùn)行結(jié)果得到兩者的計算差異(圖6d)，其范圍在10-8量級，顯然該差異可忽略不計。對比多次波壓制前(圖7a)與在CPU-GPU異構(gòu)平臺上多次波壓制后(圖7b)的疊加剖面，可見后者的多次波被衰減(紅色箭頭)，層次變得清晰，使得一次波同相軸更突出。

圖6 CPU與GPU平臺上多次波去除效果

圖7 CPU-GPU異構(gòu)計算平臺多次波去除前(a)、后(b)疊加剖面

在保證結(jié)果正確的前提下，再對CPU-GPU異構(gòu)平臺并行加速的效果進(jìn)行測試。不計數(shù)據(jù)從硬盤讀寫的時間，分別測得CPU串行、單核CPU-單GPU、12核CPU和12核CPU-雙GPU的計算耗時，其結(jié)果如表2所示。并根據(jù)本文加速比公式

(7)

分別計算不同平臺加速比。式中TR和TC對應(yīng)異構(gòu)平臺運(yùn)行時間、串行運(yùn)行時間。

圖8為三種異構(gòu)平臺的加速比。分析表2和圖8可知：當(dāng)?shù)卣饠?shù)據(jù)范圍是Inline2200～2205、Xline1500～1510(即地震道集數(shù)為6×11)時，硬件資源未充分利用，加速比不穩(wěn)定； 當(dāng)?shù)卣饠?shù)據(jù)范圍是Inline2200～2300、Xline1500～1700(地震道集數(shù)為101×201)時，硬件資源得到充分利用，且在密集計算時CPU利用率達(dá)到100%，GPU利用率最高達(dá)88%，加速比趨于穩(wěn)定。此時單核CPU—單GPU異構(gòu)平臺的加速比達(dá)到13以上；12核CPU的平臺實(shí)際加速比達(dá)到12以上；12核CPU—雙GPU的加速比也超過30。由于Windows操作系統(tǒng)不是實(shí)時的，會受到其他程序影響，加速比會小幅波動。

圖8 不同異構(gòu)計算平臺加速比

表2 不同平臺上的計算耗時

4 結(jié)束語

拋物線Radon變換算法可壓制和去除多次波，提高地震資料的信噪比。但對于數(shù)據(jù)規(guī)模巨大的地震道集，進(jìn)行拋物線Radon變換需較長運(yùn)行時間。為此，本文基于CPU-GPU異構(gòu)平臺對拋物線Radon變換算法做并行化分析，通過CUDA優(yōu)化技術(shù)在GPU平臺實(shí)施并行化，且利用多核CPU多線程技術(shù)，提出一種通用的CPU-GPU異構(gòu)并行方案，充分地利用了CPU和GPU的計算資源，在保證道集質(zhì)量的前提下，大幅度提高了該算法的計算效率，加速比達(dá)30。

成都晶石石油軟件有限公司在課題研究中提供了GEOSCOPE地質(zhì)放大鏡軟件，在此深表感謝！