熊超 王欣 王鑫杰 吳和喜

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

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

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

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

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

1 小波閾值降噪算法原理

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

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

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

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

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

1.1 小波軟閾值降噪

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

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

1.2 小波硬閾值降噪

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

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

1.3 改進的小波閾值降噪

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 并行化小波變換

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

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

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

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

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

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

3.2 并行化閾值降噪

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

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

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

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

4 性能測試及分析

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

4.1 GPU加速效果測試

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

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

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

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

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

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

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

圖6 不同數(shù)據(jù)體量下總時間加速比以及改進閾值降噪函數(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ù)體量較小時GPU 的數(shù)據(jù)傳輸時間相對較短，對總計算時間的影響較小，不會對128×128的數(shù)據(jù)體量造成顯著影響。但隨著數(shù)據(jù)體量的增加，傳輸時間逐漸增加，特別是在256×256的數(shù)據(jù)體量下，GPU中的數(shù)據(jù)傳輸所占時間顯著增加，這導(dǎo)致傳輸時間在總計算時間的占比增大，GPU的計算性能無法彌補這一時間損失，因此加速比有所下降。而當(dāng)數(shù)據(jù)體量超過256×256 時，GPU 優(yōu)越的計算性能開始得到顯現(xiàn)，盡管傳輸時間仍然在增加，但性能提升的幅度更大，傳輸時間對加速比的影響遠小于計算性能提升帶來的加速比提升。此外，對于大數(shù)據(jù)體量，CPU串行執(zhí)行的指令數(shù)量非常龐大，計算時間不斷增加，而GPU可以同時支持多達896個線程執(zhí)行，這使得GPU 在相同指令數(shù)量的情況下能夠更快地完成計算任務(wù)。

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

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

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

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

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

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

5 閾值降噪效果實例驗證

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

隨后計算降噪后信噪比，通過信噪比來表征不同方法的降噪效果，其中，降噪前的信噪比計算公式為：

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

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

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

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

6 結(jié)語

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

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