基于三維地震數(shù)據(jù)的并行希爾伯特黃變換算法研究

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

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

0 引 言

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

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

1 研究背景

1.1 HHT算法原理

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

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

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

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

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

于是可得

瞬時振幅

瞬時相位

瞬時頻率

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

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

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

2.1 GHHT的基本步驟

GHHT（GPU based Hilbert－Huang transform）并行算法流程圖如圖2所示，申請保存原始數(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)行希爾伯特頻譜分析，包括瞬時振幅譜、瞬時相位譜、瞬時頻率譜；

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

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

2.2 GHHT的加速方案

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

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

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

3 實驗結(jié)果與分析

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

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

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

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

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

3.2 GPU資源分析

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

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

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

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

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

圖7 處理前后的波形

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

圖8 處理前后的切片

4 結(jié)束語

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