基于CUDA的航空γ能譜數(shù)據(jù)小波降噪并行加速算法

熊超 王欣 王鑫杰 吳和喜

1（東華理工大學(xué) 核科學(xué)與工程學(xué)院 南昌 330013）

2（蘇州大學(xué) 醫(yī)學(xué)部放射醫(yī)學(xué)與防護(hù)學(xué)院 蘇州 215123）

航空γ能譜數(shù)據(jù)通常表現(xiàn)為γ場(chǎng)背景、礦致異常和噪聲的耦合疊加，因測(cè)量時(shí)域批次性等因素影響，探測(cè)所得數(shù)據(jù)常含有以異常形態(tài)呈現(xiàn)的高值條帶。此類(lèi)噪聲和假異常疊加會(huì)導(dǎo)致礦致異常信息被淹沒(méi)，影響礦產(chǎn)資源勘探和輻射環(huán)境的評(píng)估的準(zhǔn)確性［1］。熊超等［2］利用多重分形理論對(duì)航空γ 能譜數(shù)據(jù)進(jìn)行功率譜分形分析，并通過(guò)計(jì)算得出的截止頻率設(shè)計(jì)濾波器，從而實(shí)現(xiàn)對(duì)原始數(shù)據(jù)的逐層濾波。該方法在特定頻率區(qū)間內(nèi)能夠識(shí)別已知礦點(diǎn)的礦致異常，同時(shí)能夠排除由飛行時(shí)域批次性因素引起的虛假異常。熊超等［3］提出一種依據(jù)標(biāo)準(zhǔn)差變異系數(shù)來(lái)選擇小波基函數(shù)和分解層數(shù)的方法，旨在對(duì)航空γ能譜數(shù)據(jù)進(jìn)行分解和重構(gòu)，并使用古鈾量分布情況對(duì)提取的異常區(qū)域進(jìn)行修正。該方法能夠顯著表征空間位置信息和元素含量信息，從而提取范圍更小且精度更高的異常區(qū)域，與成礦環(huán)境相吻合且無(wú)假異常區(qū)域。趙思文等［4］提出一種基于奇異譜分析（Singular Spectrum Analysis，SSA）的γ 能譜去噪算法，通過(guò)SSA 對(duì)γ 能譜數(shù)據(jù)進(jìn)行預(yù)處理以識(shí)別和消除噪聲，從而提高能譜數(shù)據(jù)的質(zhì)量。

與此同時(shí)，由于航空γ 能譜測(cè)量效率高且范圍廣造成實(shí)測(cè)原始數(shù)據(jù)體量巨大，航空γ 能譜數(shù)據(jù)后處理需要占用大量計(jì)算資源，包括計(jì)算機(jī)處理器和存儲(chǔ)。近年來(lái)隨著計(jì)算機(jī)硬件的飛速發(fā)展，利用圖形處理單元（Graphics Processing Unit，GPU）進(jìn)行數(shù)據(jù)并行處理和分析的方法已在諸多領(lǐng)域得到廣泛應(yīng)用。GPU具有大量計(jì)算核心，可同時(shí)并行處理多個(gè)任務(wù)從而加速數(shù)據(jù)處理進(jìn)程。許明炬等［5］提出一種基于廣義交叉驗(yàn)證（Generalized-Cross-Validation，GCV）準(zhǔn)則和通用并行計(jì)算架構(gòu)（Compute Unified Device Architecture，CUDA）技術(shù)的離散小波變換方法用于壓制地震記錄中的面波噪聲，該方法可最大程度地保留有效信號(hào)，且自適應(yīng)閾值處理降低了人?工數(shù)據(jù)處理對(duì)地震記錄的主觀偏差，提高處理結(jié)果的精確性和可靠性；Honzátko 等［6］提出了一種利用CUDA 加速的塊匹配和三維塊匹配（Block Matching and 3D filtering，BM3D）算法，使BM3D 方法更貼近實(shí)際應(yīng)用；邱霽巖等［7］提出了一種基于GPU的離散小波變換算法，在一張2048×2048分辨率的圖像中利用CUDA實(shí)現(xiàn)在GPU上的并行計(jì)算，達(dá)到最大106.34 倍的加速比；Xu 等［8］提出了一個(gè)通用的升降小波變換并行計(jì)算框架（GPCF-LWT），利用GPU 集群和CUDA 以解決在大規(guī)模測(cè)量數(shù)據(jù)和有限計(jì)算資源下實(shí)時(shí)執(zhí)行動(dòng)態(tài)時(shí)間規(guī)整（Dynamic Time Warping，DWT）的問(wèn)題。該框架在計(jì)算效率、操作穩(wěn)健性和任務(wù)通用性方面有顯著改進(jìn)；為解決離散小波變換過(guò)程耗時(shí)且不利于實(shí)際工程應(yīng)用的問(wèn)題，張金霜［9］提出利用GPU 平臺(tái)的CUDA 加速技術(shù)對(duì)小波變換算法進(jìn)行并行化改造，發(fā)現(xiàn)基于GPU/CUDA技術(shù)的并行小波Mallat算法比串行小波變換算法的執(zhí)行速度最高提升了50余倍，且算法效率與計(jì)算量呈正向關(guān)系；巴振寧等［10］研究了一種基于CUDA 編程平臺(tái)的CPU-GPU 異構(gòu)并行方法用于模擬復(fù)雜場(chǎng)地下的近斷層地震動(dòng)，模擬結(jié)果顯示，近斷層地震動(dòng)的集中性、破裂的方向性、速度脈沖和永久位移等特征以及真實(shí)地形對(duì)地震動(dòng)的影響都得到了清晰體現(xiàn)。由此可見(jiàn)，CPU-GPU異構(gòu)并行方法可以有效地提高譜元法模擬的計(jì)算效率，并適用于大尺度復(fù)雜場(chǎng)地的地震波場(chǎng)模擬。

本研究基于航空γ 能譜數(shù)據(jù)體量巨大的特點(diǎn)，擬采用基于GPU加速的小波閾值降噪算法對(duì)航空γ能譜數(shù)據(jù)進(jìn)行并行后處理計(jì)算，優(yōu)選適用于航空γ能譜數(shù)據(jù)處理的小波基函數(shù)對(duì)原始數(shù)據(jù)的高頻系數(shù)進(jìn)行調(diào)整以達(dá)到降噪效果，還原異常信息形態(tài)、降低噪聲和條帶干擾。針對(duì)小波變換和閾值降噪可并行程度高的特點(diǎn)，將并行化程度高的部分移植至GPU進(jìn)行計(jì)算，搭建CPU-GPU 異構(gòu)計(jì)算平臺(tái)，提升處理速度，降低時(shí)間成本。

1 小波閾值降噪算法原理

在非平穩(wěn)信號(hào)中，信號(hào)的頻譜會(huì)隨時(shí)間變化，傳統(tǒng)的傅里葉變換（Fourier Transform，F(xiàn)FT）無(wú)法反映這種變化，并且FFT 不適用于處理有限持續(xù)信號(hào)。而短時(shí)傅里葉變換（Short-time Fourier Transform，STFT）為傳統(tǒng)的傅里葉變換對(duì)非平穩(wěn)過(guò)程的局限性提供了新的方案：變換局部化。STFT可以提供信號(hào)在時(shí)間和頻率上的局部信息，但由于窗口固定化等原因，STFT 具有時(shí)間分辨率和頻率分辨率矛盾、選取窗函數(shù)困難和信號(hào)長(zhǎng)度有限等缺陷。面對(duì)FFT和STFT 的缺陷，小波變換（Wavelet Transform，WT）的提出克服了窗口大小不隨頻率變化等缺點(diǎn)，繼承和發(fā)展了STFT 的變換局部化思想，將傅里葉變換無(wú)限長(zhǎng)的三角函數(shù)基換成有限長(zhǎng)度且會(huì)衰減的小波基，通過(guò)平移和拉伸小波基函數(shù)，自適應(yīng)調(diào)整分析窗口大小和形狀以適應(yīng)信號(hào)的局部特征。小波變換［11］也因其優(yōu)越的時(shí)頻特性、低熵性、多分辨率特性、去相關(guān)性和基函數(shù)選擇靈活性等特點(diǎn)，在圖像降噪領(lǐng)域受到廣泛的關(guān)注。

信號(hào)經(jīng)過(guò)小波變換后，所得的小波系數(shù)包含信號(hào)中的重要信息。小波分解后，信號(hào)的小波系數(shù)通常較大，而噪聲的小波系數(shù)相對(duì)較小，且噪聲的小波系數(shù)值通常小于信號(hào)的小波系數(shù)（圖1）。通過(guò)選取適當(dāng)?shù)拈撝?，可以將小波系?shù)分為兩類(lèi)：大于閾值的小波系數(shù)被視為包含信號(hào)信息，應(yīng)予以保留；而小于閾值的小波系數(shù)則被認(rèn)為是由噪聲引起，因此將其置零以達(dá)到去噪目的。小波閾值降噪的實(shí)質(zhì)為抑制信號(hào)中無(wú)用部分、增強(qiáng)有用部分的過(guò)程。小波閾值去噪過(guò)程為：1）分解過(guò)程，即選定一種小波對(duì)信號(hào)進(jìn)行n層小波分解；2）閾值處理過(guò)程，即對(duì)分解的各層系數(shù)進(jìn)行閾值處理，獲得估計(jì)小波系數(shù)；3）重構(gòu)過(guò)程，根據(jù)去噪后的小波系數(shù)進(jìn)行小波重構(gòu)，獲得去噪后的信號(hào)（圖2）。

圖1 小波降噪原理示意圖 (a) 小波變換分解過(guò)程，(b) 小波變換重構(gòu)過(guò)程Fig.1 Principle diagram of wavelet denoising(a) Decomposition process of wavelet transform, (b) Reconstruction process of wavelet transform

圖2 小波閾值降噪過(guò)程Fig.2 Wavelet threshold denoising process

其中，閾值和閾值函數(shù)的選擇是小波閾值去噪的關(guān)鍵，直接影響著重構(gòu)信號(hào)的質(zhì)量。傳統(tǒng)的閾值降噪有小波軟閾值降噪和小波硬閾值降噪兩種方法。

1.1 小波軟閾值降噪

軟閾值降噪函數(shù)如下：

式中：w是系數(shù)；λ是閾值；sgn()是符號(hào)函數(shù)，當(dāng)小波系數(shù)的絕對(duì)值小于給定閾值時(shí)，令其為0，大于閾值時(shí)，系數(shù)為正則減去閾值，系數(shù)為負(fù)則加上閾值。

1.2 小波硬閾值降噪

硬閾值降噪函數(shù)如下：

式中：當(dāng)小波系數(shù)的絕對(duì)值小于給定閾值時(shí)，令其為0，大于閾值時(shí)，系數(shù)值不變。

1.3 改進(jìn)的小波閾值降噪

在對(duì)原始信號(hào)進(jìn)行小波變換后，噪聲通常存在于水平高頻系數(shù)、垂直高頻系數(shù)和對(duì)角高頻系數(shù)中。為實(shí)現(xiàn)二維離散小波閾值降噪，需將這些高頻系數(shù)與預(yù)設(shè)的閾值進(jìn)行比較和處理。對(duì)高頻分量進(jìn)行適當(dāng)處理后，便可通過(guò)重構(gòu)過(guò)程以獲取降噪后的重構(gòu)數(shù)據(jù)［12］。

然而，當(dāng)采用硬閾值函數(shù)時(shí)，由于其不連續(xù)性可能會(huì)導(dǎo)致重構(gòu)信號(hào)出現(xiàn)局部震蕩，表現(xiàn)出偽吉布斯現(xiàn)象。而當(dāng)采用軟閾值函數(shù)時(shí)，則會(huì)出現(xiàn)與真實(shí)小波函數(shù)系數(shù)之間的恒定偏差，從而導(dǎo)致重構(gòu)后信號(hào)的精度下降。因此這里需要對(duì)閾值函數(shù)進(jìn)行改進(jìn)，參考曹棟等［13］提出的一種改進(jìn)的小波閾值函數(shù)：

式中：μ是改進(jìn)閾值函數(shù)的調(diào)整因子，且μ大于0。

2 小波閾值降噪算法的實(shí)現(xiàn)

在進(jìn)行閾值降噪處理前需先確定降噪閾值。閾值的選取是小波閾值降噪的核心流程之一［14］，如果選擇的閾值過(guò)大，會(huì)導(dǎo)致有用的信號(hào)被當(dāng)作噪聲濾除；閾值過(guò)小則容易導(dǎo)致噪聲的濾除不夠徹底。閾值選取的最優(yōu)結(jié)果是剛好大于噪聲的最大水平，這里采用曹棟等人提出的閾值選取方法：

式中：j是分解層數(shù)；W是高頻分量。

由§1.3 可知，數(shù)據(jù)和濾波器序列傳輸至GPU 的全局內(nèi)存中時(shí)，經(jīng)過(guò)行列方向的一維離散小波重構(gòu)，會(huì)產(chǎn)生近似分量（cLL）、垂直細(xì)節(jié)分量（cLH）、水平細(xì)節(jié)分量（cHL）和對(duì)角細(xì)節(jié)分量（cHH），后三者均為高頻分量。cHH可以理解為是信號(hào)與高頻濾波序列兩次卷積的結(jié)果，更接近噪聲，故在閾值求解中，令W等于cHH。通過(guò)該方法可以通過(guò)程序計(jì)算求解閾值，并將閾值用于后續(xù)的閾值降噪函數(shù)中。

在多尺度小波變換中，傳統(tǒng)的閾值選擇策略是全局閾值策略，即對(duì)多層分解后的所有高頻分量采用相同的固定閾值。但不同分解層級(jí)的分量可能需要不同的閾值，因此全局閾值在這種情況下可能不適用。此外，全局閾值策略也無(wú)法對(duì)不同尺度的子信號(hào)進(jìn)行差異化處理，可能導(dǎo)致某些尺度的子信號(hào)被過(guò)度壓縮或保留過(guò)多細(xì)節(jié)。因此，此處采用局部自適應(yīng)閾值選擇策略，即分層閾值選擇策略對(duì)每一層分別選取不同的閾值。該策略可以更好地滿足不同尺度下信號(hào)降噪的需求。

基于上述分析及式（1）～（3）構(gòu)建的閾值降噪算法流程圖如圖3所示，其中硬閾值降噪、軟閾值降噪和改進(jìn)后的閾值降噪均是通過(guò)讀取每一層的細(xì)節(jié)信號(hào)值，而后將各細(xì)節(jié)值與閾值對(duì)比，并在滿足條件后改變細(xì)節(jié)值。高頻信號(hào)經(jīng)過(guò)處理后再與cLL重構(gòu)得到降噪后的結(jié)果。

圖3 閾值降噪算法流程圖Fig.3 Flowchart of threshold denoising algorithm

3 基于CUDA 的小波變換及降噪并行算法的實(shí)現(xiàn)

中央處理器（Central Processing Unit，CPU）注重靈活性和高效的任務(wù)切換。它使用了更多的晶體管來(lái)優(yōu)化控制邏輯和緩存管理。GPU 則專(zhuān)注于并行效率，GPU 在算術(shù)邏輯單元上使用大量的晶體管，因此擁有比CPU更多的處理核心，這也是GPU能夠?qū)崿F(xiàn)大規(guī)模小波降噪并行計(jì)算的基礎(chǔ)［15］。在小波變換算法實(shí)現(xiàn)層面，通過(guò)循環(huán)迭代，每次取與濾波器序列長(zhǎng)度一致數(shù)量的航空γ能譜數(shù)據(jù)進(jìn)行卷積運(yùn)算生成低頻和高頻分量的系數(shù)。為解決大循環(huán)導(dǎo)致的時(shí)間復(fù)雜度增加問(wèn)題，通過(guò)CUDA 中的多線程計(jì)算實(shí)現(xiàn)并行化小波變換，并定義在GPU 中運(yùn)行的核函數(shù)，通過(guò)線程索引數(shù)據(jù)位置實(shí)現(xiàn)簡(jiǎn)化大循環(huán)，提高性能。而在閾值降噪方面，則是通過(guò)優(yōu)化處理實(shí)現(xiàn)并行化，利用CPU-GPU 異構(gòu)平臺(tái)，將濾波器序列定義為常量以提高速度。并在GPU 中定義與原始數(shù)據(jù)數(shù)組大小一致的內(nèi)存，通過(guò)CUDA 函數(shù)進(jìn)行數(shù)據(jù)傳輸，并調(diào)用核函數(shù)進(jìn)行計(jì)算，確定線程網(wǎng)格和線程塊的大小以提高閾值降噪的計(jì)算速度。

簡(jiǎn)而言之，利用CUDA中的多線程計(jì)算和GPU的并行計(jì)算能力將降噪任務(wù)分解為獨(dú)立的線程并在GPU 上執(zhí)行以提高小波變換和閾值降噪的計(jì)算性能。這里需要注意的是，二維離散小波變換需先從行方向再?gòu)牧蟹较蜻M(jìn)行一維離散小波變換，重構(gòu)則相反。

3.1 并行化小波變換

每次循環(huán)取選取與濾波器序列長(zhǎng)度一致數(shù)量的數(shù)據(jù)，隨后與對(duì)應(yīng)濾波器序列對(duì)應(yīng)位置的數(shù)據(jù)相乘，累加后進(jìn)行卷積運(yùn)算，形成低頻或者高頻分量的一個(gè)系數(shù)，然后進(jìn)行下一次循環(huán)。多個(gè)循環(huán)的疊加會(huì)使計(jì)算量驟增，加大了時(shí)間復(fù)雜度O（f（n）），但在變換過(guò)程中，每次大循環(huán)可以通過(guò)獨(dú)立的單線程執(zhí)行。根據(jù)這一思路可以并行化小波變換［16］，即并行化行方向與列方向一維小波變換，每個(gè)線程索引各自的數(shù)據(jù)與濾波器序列卷積，簡(jiǎn)化大循環(huán)。同理閾值降噪是讓每一個(gè)高頻系數(shù)與閾值作比對(duì)。在CPU 中計(jì)算這一過(guò)程時(shí)需通過(guò)一個(gè)大循環(huán)索引數(shù)據(jù)位置，得到高頻系數(shù)。由于過(guò)程中需要對(duì)全部高頻系數(shù)索引，故而在對(duì)每一層系數(shù)索引時(shí)，循環(huán)次數(shù)等于單個(gè)高頻分量的行數(shù)乘以列數(shù)（寬width×高h(yuǎn)eight），即矩陣階數(shù)。這意味著循環(huán)次數(shù)和時(shí)間復(fù)雜度會(huì)隨著航空γ能譜原始數(shù)據(jù)體量的增加呈幾何式增長(zhǎng)。而通過(guò)并行多個(gè)單線程執(zhí)行索引與之比對(duì)，所有線程相互獨(dú)立且沒(méi)有數(shù)據(jù)關(guān)聯(lián)及通信，在高效執(zhí)行循環(huán)過(guò)程的同時(shí)可大幅降低時(shí)間復(fù)雜度。

在CUDA 中，使用多線程進(jìn)行計(jì)算主要通過(guò)定義在GPU 中運(yùn)行的函數(shù)，為實(shí)現(xiàn)CPU-GPU 異構(gòu)平臺(tái)的搭建，使用_global_限定符定義在CPU 中調(diào)用且在GPU 中運(yùn)行的核函數(shù)。函數(shù)的正常運(yùn)行需預(yù)先定義線程的索引，該索引表示線程在網(wǎng)格中所處的位置，網(wǎng)格維度是一維的，也可以是二維的，此處主要使用二維網(wǎng)格。對(duì)線程的索引為：

int col = threadIdx.x + blockIdx.x * blockDim.x；

int row = threadIdx.y + blockIdx.y * blockDim.y

后來(lái)，在酒店里。丁小強(qiáng)提出想和杜一朵一起洗澡。杜一朵不干，說(shuō)原來(lái)你的演出就是一起洗澡？太沒(méi)有情調(diào)了。杜一朵就提出打牌，茶幾上的確有一副紙牌。她說(shuō)，我們來(lái)打牌，輸了就罰酒。

其中，col 和row 表示信號(hào)數(shù)據(jù)的坐標(biāo)，進(jìn)行數(shù)據(jù)索引，即data［col］［row］，threadIdx 是線程的索引，blockIdx是線程塊的索引，blockDim是線程塊維度。具體的數(shù)據(jù)索引按照（col，row）進(jìn)行索引訪問(wèn)，并通過(guò)（col，row）確定結(jié)果的位置，將結(jié)果傳輸?shù)较鄳?yīng)的位置。此時(shí)每個(gè)線程對(duì)數(shù)據(jù)的訪問(wèn)并無(wú)關(guān)聯(lián)性，可以獨(dú)立執(zhí)行單個(gè)任務(wù)，即訪問(wèn)數(shù)據(jù)后與濾波器卷積。

3.2 并行化閾值降噪

并行化閾值降噪需要對(duì)運(yùn)行在CPU 中的步驟進(jìn)行優(yōu)化處理，而在CPU-GPU異構(gòu)平臺(tái)上處理可以將速度提升到最優(yōu)。并行化的小波閾值降噪流程圖如圖4 所示，在運(yùn)行kernel 函數(shù)之前，需將主機(jī)端的數(shù)據(jù)傳輸?shù)皆O(shè)備端的全局內(nèi)存。因此，首先在GPU中定義與主機(jī)數(shù)據(jù)數(shù)組內(nèi)存大小一致的內(nèi)存，用于存儲(chǔ)從主機(jī)端傳輸?shù)臄?shù)據(jù)。在CUDA 中，可以使用cudaMallocPitch 函數(shù)來(lái)定義這個(gè)數(shù)組，并且通過(guò)cudaMemcpy2D 函數(shù)將數(shù)據(jù)從主機(jī)端傳輸至設(shè)備端。當(dāng)需要將設(shè)備端數(shù)據(jù)傳回主機(jī)端時(shí)，同樣可以使用cudaMemcpy2D 函數(shù)。為了提高計(jì)算速度，將濾波器序列定義為常量。對(duì)于常量的傳輸，可以使用cudaMemcpyToSymbol函數(shù)來(lái)實(shí)現(xiàn)。完成數(shù)據(jù)傳輸和結(jié)果矩陣內(nèi)存的定義后，可以調(diào)用核函數(shù)進(jìn)行計(jì)算。最后需要將計(jì)算得到的結(jié)果數(shù)據(jù)傳輸回主機(jī)端。

圖4 并行化閾值降噪流程圖Fig.4 Flowchart of parallel threshold denoising process

內(nèi)核函數(shù)的調(diào)用需要確定線程網(wǎng)格和線程塊的大小，為方便編程，CUDA 中使用dim3 類(lèi)型內(nèi)建變量threadIdx和blockIdx。dim3是用于表示三維線程塊或網(wǎng)格大小的結(jié)構(gòu)體，它包含三個(gè)unsigned int類(lèi)型的成員x、y、z，分別表示線程塊或網(wǎng)格在x、y、z方向上的大小。為實(shí)現(xiàn)線程能夠索引全部的數(shù)據(jù)，需要足夠多的線程，通過(guò)dim3結(jié)構(gòu)體定義grid和block的大?。?/p>

dim3 block_size（THREAD_X， THREAD_Y）；//線程塊block的大小

dim3 grid_size（（height） / THREAD， （width） /THREAD）；//網(wǎng)格grid大小

對(duì)于一些不可并行、數(shù)據(jù)傳輸和邏輯運(yùn)算和控制等任務(wù)將其放在CPU上執(zhí)行，而高度并行化的任務(wù)放在GPU上運(yùn)行，實(shí)現(xiàn)CPU-GPU異構(gòu)平臺(tái)，提升計(jì)算速度。

4 性能測(cè)試及分析

測(cè)試平臺(tái)CPU 采用i59300H CPU，4 核心8 線程，2.40 GHz主頻，內(nèi)存DDR43200 MHz，32 Gb，圖形顯示卡選取Nvidia GTX1650，該顯卡采用12 nm工藝，擁有896個(gè)流處理單元，核心頻率1485 MHz，最大支持線程數(shù)為1048576 個(gè)，顯存容量4 GB，顯存頻率為8000 MHz。

4.1 GPU加速效果測(cè)試

通過(guò)測(cè)試發(fā)現(xiàn)閾值降噪中選用硬閾值、軟閾值或者改進(jìn)閾值函數(shù)所消耗的時(shí)間基本一致，在計(jì)算不同block尺寸在GPU中運(yùn)行閾值降噪函數(shù)的運(yùn)行時(shí)間時(shí)，此處采用改進(jìn)閾值降噪函數(shù)進(jìn)行測(cè)試，同時(shí)小波基函數(shù)選擇為“db15”，小波分解層數(shù)為2 層。數(shù)據(jù)體量為5122時(shí)，不同block尺寸在函數(shù)中的計(jì)算時(shí)間，結(jié)果如表1所示。

表1 不同線程尺寸對(duì)于數(shù)據(jù)體量5122在改進(jìn)閾值降噪函數(shù)的計(jì)算時(shí)間Table 1 Calculation time for improved threshold denoising function with different thread sizes for a data volume of 5122

改進(jìn)閾值降噪函數(shù)計(jì)算時(shí)間的對(duì)數(shù)以及總計(jì)算時(shí)間對(duì)數(shù)隨block尺寸的變化如圖5所示。

圖5 改進(jìn)閾值降噪函數(shù)計(jì)算對(duì)數(shù)時(shí)間及總計(jì)算對(duì)數(shù)時(shí)間隨block尺寸變化圖Fig.5 Variation curves of logarithmic and total logarithmic calculation times of the improved threshold denoising function with block size

通過(guò)表1 和圖5 可以發(fā)現(xiàn)，在數(shù)據(jù)體量為5122時(shí)，最佳的block 在642～1282之間。block 尺寸小于642時(shí)，尺寸越大，計(jì)算所消耗的時(shí)間越少，在82～322之間，計(jì)算時(shí)間基本持平并且block尺寸大于1282時(shí)消耗時(shí)間也基本持平，在大于2562時(shí)，時(shí)間有所回升。并且由表1可以發(fā)現(xiàn)，數(shù)據(jù)體量為5122時(shí)，GPU所計(jì)算的時(shí)間最好可以達(dá)到35 ms 左右，對(duì)于人類(lèi)反應(yīng)時(shí)間200 ms來(lái)說(shuō)，遠(yuǎn)小于達(dá)到實(shí)時(shí)處理對(duì)人類(lèi)反應(yīng)時(shí)間的要求，這為后續(xù)航空γ 能譜數(shù)據(jù)實(shí)時(shí)處理提供了理論依據(jù)，也驗(yàn)證了董犖等［17］得出的block中thread數(shù)是32的倍數(shù)時(shí)速度更快的結(jié)論。

對(duì)于不同數(shù)據(jù)體量，GPU 加速的效果不同，這里分別使用1282、2562、5122、10242、20482和40962大小體量的數(shù)據(jù)，利用上述條件分別計(jì)算各函數(shù)在CPU中執(zhí)行和在GPU中的執(zhí)行時(shí)間和加速比（GPU執(zhí)行時(shí)間除以CPU執(zhí)行時(shí)間）。

對(duì)改進(jìn)閾值降噪函數(shù)加速比對(duì)數(shù)據(jù)體量成圖以及總時(shí)間加速比隨數(shù)據(jù)體量的變化如圖6所示。

圖6 不同數(shù)據(jù)體量下總時(shí)間加速比以及改進(jìn)閾值降噪函數(shù)加速比的變化Fig.6 Changes of total time acceleration ratio and acceleration ratio of improved threshold denoising function with data volume

由表2 和圖6 可知，在256×256 的數(shù)據(jù)體量下，與之前的128×128 相比，加速比有所降低。這是由于在數(shù)據(jù)體量較小時(shí)GPU 的數(shù)據(jù)傳輸時(shí)間相對(duì)較短，對(duì)總計(jì)算時(shí)間的影響較小，不會(huì)對(duì)128×128的數(shù)據(jù)體量造成顯著影響。但隨著數(shù)據(jù)體量的增加，傳輸時(shí)間逐漸增加，特別是在256×256的數(shù)據(jù)體量下，GPU中的數(shù)據(jù)傳輸所占時(shí)間顯著增加，這導(dǎo)致傳輸時(shí)間在總計(jì)算時(shí)間的占比增大，GPU的計(jì)算性能無(wú)法彌補(bǔ)這一時(shí)間損失，因此加速比有所下降。而當(dāng)數(shù)據(jù)體量超過(guò)256×256 時(shí)，GPU 優(yōu)越的計(jì)算性能開(kāi)始得到顯現(xiàn)，盡管傳輸時(shí)間仍然在增加，但性能提升的幅度更大，傳輸時(shí)間對(duì)加速比的影響遠(yuǎn)小于計(jì)算性能提升帶來(lái)的加速比提升。此外，對(duì)于大數(shù)據(jù)體量，CPU串行執(zhí)行的指令數(shù)量非常龐大，計(jì)算時(shí)間不斷增加，而GPU可以同時(shí)支持多達(dá)896個(gè)線程執(zhí)行，這使得GPU 在相同指令數(shù)量的情況下能夠更快地完成計(jì)算任務(wù)。

表2 不同體量的數(shù)據(jù)在GPU和CPU中的加速時(shí)間Table 2 Acceleration times for different volumes of data on GPU and CPU

4.2 GPU加速處理航空γ能譜數(shù)據(jù)效果分析

不同小波基函數(shù)的濾波序列長(zhǎng)度不一樣，分解一層后的數(shù)據(jù)大小與上一層數(shù)據(jù)大小的關(guān)系為：

因此，不同長(zhǎng)度的濾波序列下，下一層的數(shù)據(jù)大小不一，時(shí)間損耗的結(jié)果也會(huì)有差異，這種差異會(huì)對(duì)GPU加速比產(chǎn)生一定的影響。此處數(shù)據(jù)選用1312×101 大小的航空γ 能譜原始數(shù)據(jù)集，采用4 層分解。為測(cè)試不同小波基函數(shù)對(duì)航空γ 能譜數(shù)據(jù)的適用性，這里選取45種不同小波基函數(shù)分別對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行改進(jìn)閾值降噪加速測(cè)試，不同小波基下改進(jìn)閾值降噪函數(shù)的GPU加速比如表3所示。

表3 不同小波基函數(shù)下改進(jìn)閾值降噪函數(shù)的GPU加速時(shí)間比Table 3 GPU acceleration time ratios for improved threshold denoising function with different wavelet basis functions

由表3可知，隨著濾波序列長(zhǎng)度的增長(zhǎng)，各函數(shù)的加速比也在不斷增加。其中，coif5小波基函數(shù)的加速比最好，高達(dá)350倍左右，且對(duì)改進(jìn)閾值降噪函數(shù)達(dá)到了570 倍左右的加速比，而coif5 小波基函數(shù)也是所有小波中濾波器序列最長(zhǎng)的小波，濾波器序列長(zhǎng)度為30，說(shuō)明不同濾波器序列對(duì)加速的效果也有較大的影響。同時(shí)，相比于總數(shù)據(jù)插值后的數(shù)據(jù)體量或者在核應(yīng)急和輻射環(huán)境評(píng)估中低空無(wú)人機(jī)飛行的數(shù)據(jù)，此處所選用的1312×101 數(shù)據(jù)體量較小。數(shù)據(jù)實(shí)驗(yàn)表明，在GPU算力范圍內(nèi)，數(shù)據(jù)體量越大，加速比越高。因此，可以認(rèn)為，將該算法應(yīng)用于大區(qū)域航空γ能譜數(shù)據(jù)中加速比將更為可觀。

5 閾值降噪效果實(shí)例驗(yàn)證

實(shí)例選取AGS-863型航空γ能譜儀在內(nèi)蒙古錫林郭勒盟某礦區(qū)測(cè)量得到的試驗(yàn)飛行數(shù)據(jù)。為了更好地體現(xiàn)降噪效果，此處人工添加均值為0、標(biāo)準(zhǔn)差為1的白噪聲，噪聲分布如圖7所示。

隨后計(jì)算降噪后信噪比，通過(guò)信噪比來(lái)表征不同方法的降噪效果，其中，降噪前的信噪比計(jì)算公式為：

式中：xi，j是原始信號(hào)值；yi，j是染噪后的信號(hào)值，降噪后的信噪比計(jì)算公式為：

式中：zi，j是降噪后的信號(hào)值。通過(guò)公式計(jì)算，發(fā)現(xiàn)bior2.4 小波基函數(shù)降噪后的信噪比值最小，降噪后效果如圖8 所示，bior3.1 小波基函數(shù)降噪后的信噪比值最大，其降噪效果如圖9所示。

圖8 bior2.4小波基函數(shù)的三種閾值降噪的降噪效果(a) 原始數(shù)據(jù)，(b) 軟閾值，(c) 硬閾值，(d) 改進(jìn)閾值Fig.8 Denoising effects of three thresholding methods with the bior2.4 wavelet basis function(a) Original data, (b) Soft thresholding, (c) Hard thresholding,(d) Improved thresholding

圖9 bior3.1小波基函數(shù)的三種閾值降噪的降噪效果(a) 原始數(shù)據(jù)，(b) 軟閾值，(c) 硬閾值，(d) 改進(jìn)閾值Fig.9 Denoising effects of three thresholding methods with the bior3.1 wavelet basis function(a) Original data, (b) Soft thresholding, (c) Hard thresholding,(d) Improved thresholding

圖8 中，bior2.4 小波基函數(shù)對(duì)測(cè)區(qū)西南方的噪聲和條帶有所消除，軟閾值的消除效果最好，但同時(shí)丟失了部分異常，硬閾值和改進(jìn)閾值在消除條帶的同時(shí)，保留了一些基本地質(zhì)特征，由圖8 可知，三種閾值降噪對(duì)條帶消除的效果均較好，但對(duì)于測(cè)區(qū)北部數(shù)據(jù)降噪處理的結(jié)果并不理想，基本與原始圖像一致；圖9 中，使用bior3.1 小波基函數(shù)，出現(xiàn)了過(guò)度降噪導(dǎo)致的圖像失真現(xiàn)象，并且引入了行方向條帶，硬閾值在測(cè)區(qū)西南方還引入了大面積原始數(shù)據(jù)中不存在的噪聲。因此，這里選取降噪后信噪比在中間位置的小波基濾波器序列，其降噪效果如圖10所示。

圖10 不同小波基函數(shù)下各閾值降噪效果(a) 原始數(shù)據(jù)圖，(b) 基于小波基函數(shù)bior3.7的改進(jìn)閾值降噪效果，(c) 基于小波基函數(shù)coif1的軟閾值降噪，(d) 基于小波基函數(shù)coif5的硬閾值降噪Fig.10 Different denoising effects with various wavelet basis functions (a) Original data, (b) Improved threshold denoising based on bior3.7, (c) Soft threshold denoising based on coif1,(d) Hard threshold denoising based on coif5

圖10 中，三種閾值函數(shù)在東南方主條帶處的降噪效果均不佳，通過(guò)與實(shí)測(cè)區(qū)域地質(zhì)圖對(duì)比，發(fā)現(xiàn)此處是其他測(cè)量誤差導(dǎo)致的假異常，且此處軟閾值降噪出現(xiàn)失真，反而衍生出條帶；其次，在西南方位置處的大部分條帶消除效果很好，但是硬閾值相較于其他兩種函數(shù)，圖像更加粗糙；與此同時(shí)，在測(cè)區(qū)西北部三種函數(shù)對(duì)條帶的消除均取得了顯著的效果，但在特征信息的保留上，改進(jìn)閾值效果更佳。綜上所述，改進(jìn)閾值降噪對(duì)于條帶噪聲的消除情況比其他兩種函數(shù)好。

6 結(jié)語(yǔ)

針對(duì)航空γ能譜數(shù)據(jù)體量巨大而CPU執(zhí)行數(shù)據(jù)后處理時(shí)計(jì)算效率低的問(wèn)題，采用了基于GPU的二維離散小波閾值并行降噪技術(shù)對(duì)航空γ能譜數(shù)據(jù)進(jìn)行處理。首先進(jìn)行了GPU加速效果的測(cè)試，結(jié)果顯示，不同的block 尺寸對(duì)GPU 計(jì)算時(shí)間產(chǎn)生了顯著影響；其次測(cè)試了不同數(shù)據(jù)體量下的加速效果，結(jié)果發(fā)現(xiàn)，隨著數(shù)據(jù)體量的增大，相較于CPU，GPU的加速比不斷增加，但在數(shù)據(jù)體量為256×256 時(shí)加速比相對(duì)于128×128 和512×512 較低，這是由于數(shù)據(jù)傳輸所帶來(lái)的影響；為了更充分利用GPU 性能，建議盡可能增大數(shù)據(jù)體量。在航空γ能譜數(shù)據(jù)處理測(cè)試中，不同小波基函數(shù)對(duì)加速比亦產(chǎn)生了顯著影響。除了部分濾波器序列較短的小波基函數(shù)外，80%的小波基函數(shù)總時(shí)間加速比都達(dá)到了100 以上，91%的小波基函數(shù)總時(shí)間加速比達(dá)到90倍以上，特別是coif5 小波基函數(shù)的加速比達(dá)到350 倍，對(duì)于閾值降噪函數(shù)的加速比接近570 倍；最后對(duì)不同降噪函數(shù)的處理結(jié)果進(jìn)行了對(duì)比，結(jié)果顯示：所有函數(shù)都存在信噪比較低時(shí)降噪不足和信噪比較高時(shí)過(guò)度降噪的情況；使用硬閾值選取小波基函數(shù)coif5、軟閾值選取小波基函數(shù)coif1 以及改進(jìn)閾值選取小波基函數(shù)bior3.7 進(jìn)行處理，都取得了顯著的降噪效果。其中改進(jìn)閾值的效果最佳，但在異常條帶部分的消除效果相對(duì)較差，后續(xù)需要進(jìn)一步研究改進(jìn)。

作者貢獻(xiàn)聲明熊超負(fù)責(zé)研究的提出及設(shè)計(jì)、數(shù)據(jù)的收集和整理；王欣負(fù)責(zé)文章的起草和最終版本的修訂；王鑫杰負(fù)責(zé)程序設(shè)計(jì)及實(shí)驗(yàn)數(shù)據(jù)的處理；吳和喜負(fù)責(zé)最終版本的修訂和項(xiàng)目的監(jiān)督及管理。