阮寧君，陳 俊，于文茂，余厚全，2，羅明璋，2，謝 凱，2

（1.長(zhǎng)江大學(xué) 電子信息學(xué)院油氣信息處理與識(shí)別研究所，湖北 荊州434023；2.油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，湖北 荊州434023）

0 引 言

在過(guò)去50年里，國(guó)內(nèi)油氣田被大規(guī)模的開(kāi)發(fā)，現(xiàn)在未被發(fā)現(xiàn)的油氣藏大都儲(chǔ)量較小、埋藏深、孔滲條件復(fù)雜且難于發(fā)現(xiàn)，那么為了更清楚的 “看見(jiàn)”那些存儲(chǔ)復(fù)雜儲(chǔ)量小的油氣藏，就必須對(duì)采集到的信號(hào)做出更準(zhǔn)確更精細(xì)的處理，以獲得更詳盡地層的信息。此前，在地震信號(hào)的時(shí)頻分析［1－2］領(lǐng)域里最常用的信號(hào)處理方法是傅立葉分析［3］和小波分析［4－5］。由于地震信號(hào)屬于非線性和非平穩(wěn)信號(hào)［6］，這就暴露了傅立葉分析和小波分析的許多自身難以克服的缺陷。希爾伯特［7－10］黃變換是基于信號(hào)局部特征的，并且是自適應(yīng)的，它獨(dú)特的處理方法不但使得它能識(shí)別出信號(hào)的不同震蕩模式，而且大大的提高了地震信號(hào)的分辨率。因此，黃變換在地震信號(hào)的分析和處理中具有非常重要的意義。

隨著處理內(nèi)容和算法復(fù)雜度的不斷增加，這為計(jì)算的速度帶來(lái)了巨大挑戰(zhàn)。并行處理技術(shù)日益引起石油地球物理界的廣泛關(guān)注，如何快速高效地并行處理大規(guī)模地震數(shù)據(jù)這一問(wèn)題已成為亟待解決的重大課題之一。當(dāng)前值得關(guān)注的是GPU［11－13］計(jì)算能力的飛速發(fā)展及其應(yīng)用的擴(kuò)展，特別是基于GPU并行處理技術(shù)的形成與發(fā)展，為地震數(shù)據(jù)的并行處理提供了平臺(tái)。任何有C語(yǔ)言基礎(chǔ)的用戶都很容易地開(kāi)發(fā) CUDA（compute unified device architecture）［14－15］的應(yīng)用程序?；贕PU的并行地震數(shù)據(jù)處理不僅大大的提高了性價(jià)比，而且為海量地震數(shù)據(jù)的處理提供了可能。

1 研究背景

1.1 HHT算法原理

HHT算法是由經(jīng)驗(yàn)?zāi)B(tài)分解和希爾伯特變換兩步組成：

（1）經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）

希爾伯特變換（hilbert huang transform，HHT）核心是對(duì)信號(hào)進(jìn)行 經(jīng)驗(yàn)?zāi)B(tài)分解 ，獲得有限數(shù)目的固有模態(tài)函數(shù)（IMF）。為保證IMF是單分量函數(shù) ，它必須滿足下列條件：① 極值點(diǎn)的數(shù)目和零交叉點(diǎn)的數(shù)目相同或者至多相差一個(gè)；②所有瞬時(shí)平均包絡(luò)序列M（t）的平方和與原序列Y（t）的平方和的比值小于0.1。

（2）希爾伯特頻譜分析（HSA）

對(duì)于任意一個(gè)時(shí)間序列X（t），都能得到它的Hilbert變換結(jié)果Y（t），即

于是可得

瞬時(shí)振幅

瞬時(shí)相位

瞬時(shí)頻率

1.2 CUDA的體系結(jié)構(gòu)

統(tǒng)一計(jì)算設(shè)備架構(gòu)（compute unified device architecture，CUDA）是NVIDIA公司提出的一種新的處理和管理GPU計(jì)算的硬件和軟件架構(gòu)，其是以C編程為基礎(chǔ)，可以直接用C語(yǔ)言進(jìn)行應(yīng)用程序的開(kāi)發(fā)。在CUDA的架構(gòu)下可以分為host端和device端，如圖1所示，device端是在GPU上執(zhí)行的部份，Host端則是在CPU上執(zhí)行的部份。在Device上執(zhí)行的程序稱為"kernel"。一般host端的程序會(huì)將數(shù)據(jù)準(zhǔn)備好后，復(fù)制到顯存中，再由GPU執(zhí)行device端程序完成對(duì)數(shù)據(jù)的處理，處理完后再由host端程序?qū)⒔Y(jié)果從顯存中取回到內(nèi)存中。

2 基于GPU的希爾伯特黃變換算法

2.1 GHHT的基本步驟

GHHT（GPU based Hilbert－Huang transform）并行算法流程圖如圖2所示，申請(qǐng)保存原始數(shù)據(jù)和中間變量的設(shè)備內(nèi)存。

圖1 CUDA執(zhí)行模型

圖2 GHHT算法流程

（1）將原始數(shù)據(jù)和中間變量的值拷貝到設(shè)備內(nèi)存中；

（2）執(zhí)行極大值與極小值的判斷、三次樣條插值、平均包絡(luò)；

（3）IMF分量終止條件的判斷；

（4）判斷條件是否滿足，是否退出循環(huán)；

（5）將分解后的數(shù)據(jù)進(jìn)行希爾伯特頻譜分析，包括瞬時(shí)振幅譜、瞬時(shí)相位譜、瞬時(shí)頻率譜；

（6）將結(jié)果從設(shè)備內(nèi)存中拷貝的主機(jī)內(nèi)存中。

在本次試驗(yàn)中，每一道數(shù)據(jù)有一個(gè)線程負(fù)責(zé)處理，但是雖說(shuō)是并行計(jì)算，并沒(méi)有滿足硬件設(shè)備的一些要求，且沒(méi)有發(fā)揮最大的運(yùn)算能力，硬件資源沒(méi)有充分的利用，所以加速比并不理想，GPU速度和CPU速度相當(dāng)。

2.2 GHHT的加速方案

（1）存儲(chǔ)器優(yōu)化：①在利用CUDA并行加速時(shí)。既要充分的利用硬件資源（如共享內(nèi)存、常量存儲(chǔ)器、寄存器、紋理緩存），也要合理的訪問(wèn)達(dá)到最高效率。其存儲(chǔ)器的存儲(chǔ)模型如圖4所示，其中包括多種不同類型的存儲(chǔ)器。每個(gè)thread都有自己的一定大小的register和local memory的空間。同一個(gè)block中的每個(gè)thread則有共享的一份share memory。此外，所有的thread（包括不同block的thread）都共享一份 global memory、constant memory、和texture memory。不同的grid則有各自的global memory、constant memory 和 texture memory。② 同 一 block 中 的thread可以通過(guò)共享內(nèi)存進(jìn)行通訊。不同的block可以通過(guò)全局存儲(chǔ)器進(jìn)行通訊。在不同的條件下合理搭配發(fā)揮存儲(chǔ)器的最大效率。例如共享內(nèi)存訪問(wèn)速度快，但是必須對(duì)齊訪問(wèn)和防止bank沖突，才能達(dá)到最佳效率，適用于多次重復(fù)訪問(wèn)的數(shù)據(jù)。常量存儲(chǔ)器擁有的一定大小的緩存的只讀存儲(chǔ)器，適用于多次隨機(jī)訪問(wèn)的常量。寄存器訪問(wèn)速度快，適用于存放中間變量，但數(shù)量相當(dāng)有限。每個(gè)線程用自己的可讀寫的寄存器于本地內(nèi)存。在同一塊中的線程可以通過(guò)讀寫共享內(nèi)存通訊，所有的線程可以讀寫全局的內(nèi)存（常量?jī)?nèi)存、紋理緩存、顯存）通信。能合理的組織數(shù)據(jù)的存儲(chǔ)至關(guān)重要。③本文將不變的常數(shù)如地震數(shù)據(jù)采樣點(diǎn)的個(gè)數(shù)、道數(shù)、切片數(shù)放入常量存儲(chǔ)器中，將存取頻率最大的地震數(shù)據(jù)放入共享存儲(chǔ)器中。為了充分的利用寄存器，本文將一些中間變量放入寄存器中如在三次樣條時(shí)會(huì)產(chǎn)生許多中間變量。

（2）線程模式優(yōu)化

在CUDA架構(gòu)下，GPU執(zhí)行時(shí)的最小單位是thread（線程）。Grid、block和thread的關(guān)系如圖3所示，多個(gè)thread可以組成一個(gè)block（塊）。而多個(gè)block組成一個(gè)grid。每一個(gè)block包含有限的thread數(shù)目。讓每個(gè)block處理一道數(shù)據(jù)，這樣就可以充分的利用共享內(nèi)存如圖4所示。把讀取頻繁的數(shù)據(jù)和一些中間變量放入共享內(nèi)存中，這樣就可以減少一些中間變量的存取時(shí)間。IMF判斷條件主要的計(jì)算量在于求過(guò)零點(diǎn)的個(gè)數(shù)和平均包絡(luò)與原始數(shù)據(jù)能量的比值。由每個(gè)block處理一道數(shù)據(jù)，線程之間的協(xié)調(diào)和通信至關(guān)重要。我們借助共享內(nèi)存將每個(gè)線程 “產(chǎn)生”的數(shù)據(jù)求和，在求和時(shí)采用樹(shù)狀加法如圖5所示，可以避免bank沖突達(dá)到最大效率。

3 實(shí)驗(yàn)結(jié)果與分析

利用GPU實(shí)現(xiàn)了對(duì)經(jīng)驗(yàn)?zāi)B(tài)分解的加速，為了驗(yàn)證本文算法的有效性，實(shí)驗(yàn)環(huán)境配置如下：

本次實(shí)驗(yàn)采用的設(shè)備GPU是GT240顯卡，有12個(gè)多核處理器、96個(gè)核心、512M內(nèi)存、15384個(gè)寄存器、65536b常量?jī)?nèi)存，時(shí)鐘頻率為1.51GHz。該顯卡著色器的理論計(jì)算能力為36GFLOP/s。對(duì)照實(shí)驗(yàn)采用的CPU為AMD Athlon II X2 245雙核 ，主頻2.89GHz。

3.1 GPU與CPU結(jié)果對(duì)比

在表1中，根據(jù)不同地震數(shù)據(jù)處理的時(shí)間對(duì)比，我們可以看出GPU的并行處理優(yōu)勢(shì)。與傳統(tǒng)CPU的串行處理相比，GPU并行處理的速度有了4倍左右的提升，性價(jià)比也大大提高。

表1 GPU與CPU的結(jié)果對(duì)比表

3.2 GPU資源分析

如圖6所示我們用軟件分析應(yīng)用程序測(cè)試GPU的資源占有率。在目前的GPU架構(gòu)中，所有的線程被組織成多個(gè)的warp塊，一個(gè)warp里面有32個(gè)threads，每16個(gè)線程為一個(gè)half－warp。而流多處理器是以half－warp為單位執(zhí)行的，因此每個(gè)流多處理器所分配的warp塊的數(shù)量直接影響著GPU資源的利用率。而每個(gè)流多處理器所分配的warp塊的數(shù)量又與每個(gè)塊（block）中的線程的大小密切相關(guān)。

圖6中橫坐標(biāo)表示每個(gè)block中分配線程的數(shù)目，縱坐標(biāo)表示每個(gè)流多處理器中Warp的數(shù)目。在本次試驗(yàn)中每個(gè)block分配256個(gè)線程，如圖6中三角符號(hào)的橫坐標(biāo)所示，其縱坐標(biāo)表示W(wǎng)arp的數(shù)目為32個(gè)?？芍Y源的配置已達(dá)到最佳配置狀態(tài)，說(shuō)明GPU資源已經(jīng)得到了充分的利用。

圖6 每個(gè)流多處理器當(dāng)前活動(dòng)的Warp塊數(shù)目

3.3 地震數(shù)據(jù)處理前后的結(jié)果對(duì)比

從微觀上看，處理前的波形與處理后的波形最明顯的差別在于穿越零點(diǎn)的個(gè)數(shù)，處理后的波形穿越零點(diǎn)的個(gè)數(shù)明顯增多，也就是縱向分辨率提高了，如圖7所示。

圖7 處理前后的波形

如圖8所示，左邊是處理前的數(shù)據(jù)，右邊是處理后的數(shù)據(jù)。處理前切片的斷層很模糊不容易被發(fā)現(xiàn)，處理后切片斷層很清晰易被發(fā)現(xiàn)，由此可見(jiàn)處理后地震數(shù)據(jù)的分辨率明顯提高。

圖8 處理前后的切片

4 結(jié)束語(yǔ)

本文在研究地震數(shù)據(jù)的并行處理時(shí)，采用了基于GPU的希爾伯特黃變換算法，實(shí)驗(yàn)結(jié)果表明經(jīng)過(guò)希爾伯特黃變換處理后的地震數(shù)據(jù)與處理前相比分辨率明顯提高，采用的GPU并行處理的速度與CPU串行處理的速度相比，有了大幅度的提升，大大的提高了性價(jià)比。我們?cè)谑煜ち私饬艘痪S希爾伯特黃變換的基礎(chǔ)上，今后會(huì)朝二維乃至三維方面發(fā)展。從二維到三維思想大致類似，這是黃變換的前瞻性理論，若能很好的處理，必將極大的推動(dòng)黃變換的發(fā)展，二維乃至三維的經(jīng)驗(yàn)?zāi)B(tài)分解的并行算法，有利于進(jìn)一步的發(fā)揮GPU的優(yōu)勢(shì)，從而進(jìn)一步的加速數(shù)據(jù)的處理。

［1］Ryzhak I S.Controlled multifunctional processing ofcausal signals on the basis of nonlinear（quadratic and cubic）time－frequency distributions［J］.Journal of Communications Technology and Electronics，2006，51（8）：1－5.

［2］Kostka P S，Tkacz E J.Feature selection based on time－frequency analysis in SVM classifier with rules extraction stage［C］.Munich，Germany：IFMBE Proceedings World Congress on Medical Physics and Biomedical Engineering，2009：1－5.

［3］Fethi Smach，CedricJean Paul Gauthier，et al.Generalized fourier descriptors with applications to objects recognition in SVM context ［J］.Journal of Mathematical Imaging and Vision，2008，30（1）：3－5.

［4］CHUANG Zhitao，LIU Youming.Wavelets with differential relation ［J］.Acta Mathematica Sinica－english Series，2011，27（5）：1011－1022.

［5］LI Jianping，TANG Yuanyan，YAN Zhonghong，et al.Uniform analytic construction of wavelet analysis filters based on sine and cosine trigonometric functions ［J］.Applied Mathematics and Mechanics，2001，22（5）：1－5.

［6］Seng Kah Phooi，Ang L M.New radial basis function neural network training for nonlinear and nonstationary signals ［G］.LNCS 4456：Computational Intelligence and Security，2007：1－5.

［7］Khalid Koufany.Application of Hilbert’s projective metric on symmetric cones ［J］.Acta Mathematica Sinica，2006，22（5）：1－5.

［8］Taras Banakh，Igor Zarichnyy.Topological groups and convex sets homeomorphic to non－separable Hilbert spaces ［J］.Central European Journal of Mathematics，2008，6（1）：1－5.

［9］ZHU Zheng，WANG Yilin，CAI Ping.Real－time implementation of vector Hilbert－Huang transform ［C］.International Workshop on Intelligent Networks and Intelligent Systems，2010：213－216.

［10］Pulung Waskito，Shinobu Miwa，Yasue Mitsukura，et al.Parallelizing Hilbert－Huang transform on a GPU ［C］.International Conference on Natural Computation，2010：184－190.

［11］PANG Waiman，QIN Jing，LU Yuqiang，et al.Accelerating simultaneous algebraic reconstruction technique with motion compensation using CUDA－enabled GPU ［J］.International Journal of Computer Assisted Radiology and Surgery，2011，6（2）：187－199.

［12］Elena Martinova.Realistic skin rendering on GPU ［J］.Intelligent Computer Graphics，2009，240：1－18.

［13］XIE K，WU P，YANG S.GPU and CPU cooperation parallel visualisation for large seismic data ［J］.Electronics Letters，2010，46（17）：1196－1197.

［14］NVIDIA CUDA C Programming Guide Version 4.0［S］.2011.

［15］CHE Shuai.A performance study of general－purpose applications on graphics processors using CUDA ［J］.Journal parallel and Distributed Computing，2008：1－5.