基于CUDA的航空γ能谱数据小波降噪并行加速算法

熊超 王欣 王鑫杰 吴和喜

1（东华理工大学 核科学与工程学院 南昌 330013）

2（苏州大学 医学部放射医学与防护学院 苏州 215123）

航空γ能谱数据通常表现为γ场背景、矿致异常和噪声的耦合叠加，因测量时域批次性等因素影响，探测所得数据常含有以异常形态呈现的高值条带。此类噪声和假异常叠加会导致矿致异常信息被淹没，影响矿产资源勘探和辐射环境的评估的准确性［1］。熊超等［2］利用多重分形理论对航空γ 能谱数据进行功率谱分形分析，并通过计算得出的截止频率设计滤波器，从而实现对原始数据的逐层滤波。该方法在特定频率区间内能够识别已知矿点的矿致异常，同时能够排除由飞行时域批次性因素引起的虚假异常。熊超等［3］提出一种依据标准差变异系数来选择小波基函数和分解层数的方法，旨在对航空γ能谱数据进行分解和重构，并使用古铀量分布情况对提取的异常区域进行修正。该方法能够显著表征空间位置信息和元素含量信息，从而提取范围更小且精度更高的异常区域，与成矿环境相吻合且无假异常区域。赵思文等［4］提出一种基于奇异谱分析（Singular Spectrum Analysis，SSA）的γ 能谱去噪算法，通过SSA 对γ 能谱数据进行预处理以识别和消除噪声，从而提高能谱数据的质量。

与此同时，由于航空γ 能谱测量效率高且范围广造成实测原始数据体量巨大，航空γ 能谱数据后处理需要占用大量计算资源，包括计算机处理器和存储。近年来随着计算机硬件的飞速发展，利用图形处理单元（Graphics Processing Unit，GPU）进行数据并行处理和分析的方法已在诸多领域得到广泛应用。GPU具有大量计算核心，可同时并行处理多个任务从而加速数据处理进程。许明炬等［5］提出一种基于广义交叉验证（Generalized-Cross-Validation，GCV）准则和通用并行计算架构（Compute Unified Device Architecture，CUDA）技术的离散小波变换方法用于压制地震记录中的面波噪声，该方法可最大程度地保留有效信号，且自适应阈值处理降低了人?工数据处理对地震记录的主观偏差，提高处理结果的精确性和可靠性；Honzátko 等［6］提出了一种利用CUDA 加速的块匹配和三维块匹配（Block Matching and 3D filtering，BM3D）算法，使BM3D 方法更贴近实际应用；邱霁岩等［7］提出了一种基于GPU的离散小波变换算法，在一张2048×2048分辨率的图像中利用CUDA实现在GPU上的并行计算，达到最大106.34 倍的加速比；Xu 等［8］提出了一个通用的升降小波变换并行计算框架（GPCF-LWT），利用GPU 集群和CUDA 以解决在大规模测量数据和有限计算资源下实时执行动态时间规整（Dynamic Time Warping，DWT）的问题。该框架在计算效率、操作稳健性和任务通用性方面有显著改进；为解决离散小波变换过程耗时且不利于实际工程应用的问题，张金霜［9］提出利用GPU 平台的CUDA 加速技术对小波变换算法进行并行化改造，发现基于GPU/CUDA技术的并行小波Mallat算法比串行小波变换算法的执行速度最高提升了50余倍，且算法效率与计算量呈正向关系；巴振宁等［10］研究了一种基于CUDA 编程平台的CPU-GPU 异构并行方法用于模拟复杂场地下的近断层地震动，模拟结果显示，近断层地震动的集中性、破裂的方向性、速度脉冲和永久位移等特征以及真实地形对地震动的影响都得到了清晰体现。由此可见，CPU-GPU异构并行方法可以有效地提高谱元法模拟的计算效率，并适用于大尺度复杂场地的地震波场模拟。

本研究基于航空γ 能谱数据体量巨大的特点，拟采用基于GPU加速的小波阈值降噪算法对航空γ能谱数据进行并行后处理计算，优选适用于航空γ能谱数据处理的小波基函数对原始数据的高频系数进行调整以达到降噪效果，还原异常信息形态、降低噪声和条带干扰。针对小波变换和阈值降噪可并行程度高的特点，将并行化程度高的部分移植至GPU进行计算，搭建CPU-GPU 异构计算平台，提升处理速度，降低时间成本。

1 小波阈值降噪算法原理

在非平稳信号中，信号的频谱会随时间变化，传统的傅里叶变换（Fourier Transform，FFT）无法反映这种变化，并且FFT 不适用于处理有限持续信号。而短时傅里叶变换（Short-time Fourier Transform，STFT）为传统的傅里叶变换对非平稳过程的局限性提供了新的方案：变换局部化。STFT可以提供信号在时间和频率上的局部信息，但由于窗口固定化等原因，STFT 具有时间分辨率和频率分辨率矛盾、选取窗函数困难和信号长度有限等缺陷。面对FFT和STFT 的缺陷，小波变换（Wavelet Transform，WT）的提出克服了窗口大小不随频率变化等缺点，继承和发展了STFT 的变换局部化思想，将傅里叶变换无限长的三角函数基换成有限长度且会衰减的小波基，通过平移和拉伸小波基函数，自适应调整分析窗口大小和形状以适应信号的局部特征。小波变换［11］也因其优越的时频特性、低熵性、多分辨率特性、去相关性和基函数选择灵活性等特点，在图像降噪领域受到广泛的关注。

信号经过小波变换后，所得的小波系数包含信号中的重要信息。小波分解后，信号的小波系数通常较大，而噪声的小波系数相对较小，且噪声的小波系数值通常小于信号的小波系数（图1）。通过选取适当的阈值，可以将小波系数分为两类：大于阈值的小波系数被视为包含信号信息，应予以保留；而小于阈值的小波系数则被认为是由噪声引起，因此将其置零以达到去噪目的。小波阈值降噪的实质为抑制信号中无用部分、增强有用部分的过程。小波阈值去噪过程为：1）分解过程，即选定一种小波对信号进行n层小波分解；2）阈值处理过程，即对分解的各层系数进行阈值处理，获得估计小波系数；3）重构过程，根据去噪后的小波系数进行小波重构，获得去噪后的信号（图2）。

图1 小波降噪原理示意图 (a) 小波变换分解过程，(b) 小波变换重构过程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

其中，阈值和阈值函数的选择是小波阈值去噪的关键，直接影响着重构信号的质量。传统的阈值降噪有小波软阈值降噪和小波硬阈值降噪两种方法。

1.1 小波软阈值降噪

软阈值降噪函数如下：

式中：w是系数；λ是阈值；sgn()是符号函数，当小波系数的绝对值小于给定阈值时，令其为0，大于阈值时，系数为正则减去阈值，系数为负则加上阈值。

1.2 小波硬阈值降噪

硬阈值降噪函数如下：

式中：当小波系数的绝对值小于给定阈值时，令其为0，大于阈值时，系数值不变。

1.3 改进的小波阈值降噪

在对原始信号进行小波变换后，噪声通常存在于水平高频系数、垂直高频系数和对角高频系数中。为实现二维离散小波阈值降噪，需将这些高频系数与预设的阈值进行比较和处理。对高频分量进行适当处理后，便可通过重构过程以获取降噪后的重构数据［12］。

然而，当采用硬阈值函数时，由于其不连续性可能会导致重构信号出现局部震荡，表现出伪吉布斯现象。而当采用软阈值函数时，则会出现与真实小波函数系数之间的恒定偏差，从而导致重构后信号的精度下降。因此这里需要对阈值函数进行改进，参考曹栋等［13］提出的一种改进的小波阈值函数：

式中：μ是改进阈值函数的调整因子，且μ大于0。

2 小波阈值降噪算法的实现

在进行阈值降噪处理前需先确定降噪阈值。阈值的选取是小波阈值降噪的核心流程之一［14］，如果选择的阈值过大，会导致有用的信号被当作噪声滤除；阈值过小则容易导致噪声的滤除不够彻底。阈值选取的最优结果是刚好大于噪声的最大水平，这里采用曹栋等人提出的阈值选取方法：

式中：j是分解层数；W是高频分量。

由§1.3 可知，数据和滤波器序列传输至GPU 的全局内存中时，经过行列方向的一维离散小波重构，会产生近似分量（cLL）、垂直细节分量（cLH）、水平细节分量（cHL）和对角细节分量（cHH），后三者均为高频分量。cHH可以理解为是信号与高频滤波序列两次卷积的结果，更接近噪声，故在阈值求解中，令W等于cHH。通过该方法可以通过程序计算求解阈值，并将阈值用于后续的阈值降噪函数中。

在多尺度小波变换中，传统的阈值选择策略是全局阈值策略，即对多层分解后的所有高频分量采用相同的固定阈值。但不同分解层级的分量可能需要不同的阈值，因此全局阈值在这种情况下可能不适用。此外，全局阈值策略也无法对不同尺度的子信号进行差异化处理，可能导致某些尺度的子信号被过度压缩或保留过多细节。因此，此处采用局部自适应阈值选择策略，即分层阈值选择策略对每一层分别选取不同的阈值。该策略可以更好地满足不同尺度下信号降噪的需求。

基于上述分析及式（1）～（3）构建的阈值降噪算法流程图如图3所示，其中硬阈值降噪、软阈值降噪和改进后的阈值降噪均是通过读取每一层的细节信号值，而后将各细节值与阈值对比，并在满足条件后改变细节值。高频信号经过处理后再与cLL重构得到降噪后的结果。

图3 阈值降噪算法流程图Fig.3 Flowchart of threshold denoising algorithm

3 基于CUDA 的小波变换及降噪并行算法的实现

中央处理器（Central Processing Unit，CPU）注重灵活性和高效的任务切换。它使用了更多的晶体管来优化控制逻辑和缓存管理。GPU 则专注于并行效率，GPU 在算术逻辑单元上使用大量的晶体管，因此拥有比CPU更多的处理核心，这也是GPU能够实现大规模小波降噪并行计算的基础［15］。在小波变换算法实现层面，通过循环迭代，每次取与滤波器序列长度一致数量的航空γ能谱数据进行卷积运算生成低频和高频分量的系数。为解决大循环导致的时间复杂度增加问题，通过CUDA 中的多线程计算实现并行化小波变换，并定义在GPU 中运行的核函数，通过线程索引数据位置实现简化大循环，提高性能。而在阈值降噪方面，则是通过优化处理实现并行化，利用CPU-GPU 异构平台，将滤波器序列定义为常量以提高速度。并在GPU 中定义与原始数据数组大小一致的内存，通过CUDA 函数进行数据传输，并调用核函数进行计算，确定线程网格和线程块的大小以提高阈值降噪的计算速度。

简而言之，利用CUDA中的多线程计算和GPU的并行计算能力将降噪任务分解为独立的线程并在GPU 上执行以提高小波变换和阈值降噪的计算性能。这里需要注意的是，二维离散小波变换需先从行方向再从列方向进行一维离散小波变换，重构则相反。

3.1 并行化小波变换

每次循环取选取与滤波器序列长度一致数量的数据，随后与对应滤波器序列对应位置的数据相乘，累加后进行卷积运算，形成低频或者高频分量的一个系数，然后进行下一次循环。多个循环的叠加会使计算量骤增，加大了时间复杂度O（f（n）），但在变换过程中，每次大循环可以通过独立的单线程执行。根据这一思路可以并行化小波变换［16］，即并行化行方向与列方向一维小波变换，每个线程索引各自的数据与滤波器序列卷积，简化大循环。同理阈值降噪是让每一个高频系数与阈值作比对。在CPU 中计算这一过程时需通过一个大循环索引数据位置，得到高频系数。由于过程中需要对全部高频系数索引，故而在对每一层系数索引时，循环次数等于单个高频分量的行数乘以列数（宽width×高height），即矩阵阶数。这意味着循环次数和时间复杂度会随着航空γ能谱原始数据体量的增加呈几何式增长。而通过并行多个单线程执行索引与之比对，所有线程相互独立且没有数据关联及通信，在高效执行循环过程的同时可大幅降低时间复杂度。

在CUDA 中，使用多线程进行计算主要通过定义在GPU 中运行的函数，为实现CPU-GPU 异构平台的搭建，使用_global_限定符定义在CPU 中调用且在GPU 中运行的核函数。函数的正常运行需预先定义线程的索引，该索引表示线程在网格中所处的位置，网格维度是一维的，也可以是二维的，此处主要使用二维网格。对线程的索引为：

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

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

后来，在酒店里。丁小强提出想和杜一朵一起洗澡。杜一朵不干，说原来你的演出就是一起洗澡？太没有情调了。杜一朵就提出打牌，茶几上的确有一副纸牌。她说，我们来打牌，输了就罚酒。

其中，col 和row 表示信号数据的坐标，进行数据索引，即data［col］［row］，threadIdx 是线程的索引，blockIdx是线程块的索引，blockDim是线程块维度。具体的数据索引按照（col，row）进行索引访问，并通过（col，row）确定结果的位置，将结果传输到相应的位置。此时每个线程对数据的访问并无关联性，可以独立执行单个任务，即访问数据后与滤波器卷积。

3.2 并行化阈值降噪

并行化阈值降噪需要对运行在CPU 中的步骤进行优化处理，而在CPU-GPU异构平台上处理可以将速度提升到最优。并行化的小波阈值降噪流程图如图4 所示，在运行kernel 函数之前，需将主机端的数据传输到设备端的全局内存。因此，首先在GPU中定义与主机数据数组内存大小一致的内存，用于存储从主机端传输的数据。在CUDA 中，可以使用cudaMallocPitch 函数来定义这个数组，并且通过cudaMemcpy2D 函数将数据从主机端传输至设备端。当需要将设备端数据传回主机端时，同样可以使用cudaMemcpy2D 函数。为了提高计算速度，将滤波器序列定义为常量。对于常量的传输，可以使用cudaMemcpyToSymbol函数来实现。完成数据传输和结果矩阵内存的定义后，可以调用核函数进行计算。最后需要将计算得到的结果数据传输回主机端。

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

内核函数的调用需要确定线程网格和线程块的大小，为方便编程，CUDA 中使用dim3 类型内建变量threadIdx和blockIdx。dim3是用于表示三维线程块或网格大小的结构体，它包含三个unsigned int类型的成员x、y、z，分别表示线程块或网格在x、y、z方向上的大小。为实现线程能够索引全部的数据，需要足够多的线程，通过dim3结构体定义grid和block的大小：

dim3 block_size（THREAD_X， THREAD_Y）；//线程块block的大小

dim3 grid_size（（height） / THREAD， （width） /THREAD）；//网格grid大小

对于一些不可并行、数据传输和逻辑运算和控制等任务将其放在CPU上执行，而高度并行化的任务放在GPU上运行，实现CPU-GPU异构平台，提升计算速度。

4 性能测试及分析

测试平台CPU 采用i59300H CPU，4 核心8 线程，2.40 GHz主频，内存DDR43200 MHz，32 Gb，图形显示卡选取Nvidia GTX1650，该显卡采用12 nm工艺，拥有896个流处理单元，核心频率1485 MHz，最大支持线程数为1048576 个，显存容量4 GB，显存频率为8000 MHz。

4.1 GPU加速效果测试

通过测试发现阈值降噪中选用硬阈值、软阈值或者改进阈值函数所消耗的时间基本一致，在计算不同block尺寸在GPU中运行阈值降噪函数的运行时间时，此处采用改进阈值降噪函数进行测试，同时小波基函数选择为“db15”，小波分解层数为2 层。数据体量为5122时，不同block尺寸在函数中的计算时间，结果如表1所示。

表1 不同线程尺寸对于数据体量5122在改进阈值降噪函数的计算时间Table 1 Calculation time for improved threshold denoising function with different thread sizes for a data volume of 5122

改进阈值降噪函数计算时间的对数以及总计算时间对数随block尺寸的变化如图5所示。

图5 改进阈值降噪函数计算对数时间及总计算对数时间随block尺寸变化图Fig.5 Variation curves of logarithmic and total logarithmic calculation times of the improved threshold denoising function with block size

通过表1 和图5 可以发现，在数据体量为5122时，最佳的block 在642～1282之间。block 尺寸小于642时，尺寸越大，计算所消耗的时间越少，在82～322之间，计算时间基本持平并且block尺寸大于1282时消耗时间也基本持平，在大于2562时，时间有所回升。并且由表1可以发现，数据体量为5122时，GPU所计算的时间最好可以达到35 ms 左右，对于人类反应时间200 ms来说，远小于达到实时处理对人类反应时间的要求，这为后续航空γ 能谱数据实时处理提供了理论依据，也验证了董荦等［17］得出的block中thread数是32的倍数时速度更快的结论。

对于不同数据体量，GPU 加速的效果不同，这里分别使用1282、2562、5122、10242、20482和40962大小体量的数据，利用上述条件分别计算各函数在CPU中执行和在GPU中的执行时间和加速比（GPU执行时间除以CPU执行时间）。

对改进阈值降噪函数加速比对数据体量成图以及总时间加速比随数据体量的变化如图6所示。

图6 不同数据体量下总时间加速比以及改进阈值降噪函数加速比的变化Fig.6 Changes of total time acceleration ratio and acceleration ratio of improved threshold denoising function with data volume

由表2 和图6 可知，在256×256 的数据体量下，与之前的128×128 相比，加速比有所降低。这是由于在数据体量较小时GPU 的数据传输时间相对较短，对总计算时间的影响较小，不会对128×128的数据体量造成显著影响。但随着数据体量的增加，传输时间逐渐增加，特别是在256×256的数据体量下，GPU中的数据传输所占时间显著增加，这导致传输时间在总计算时间的占比增大，GPU的计算性能无法弥补这一时间损失，因此加速比有所下降。而当数据体量超过256×256 时，GPU 优越的计算性能开始得到显现，尽管传输时间仍然在增加，但性能提升的幅度更大，传输时间对加速比的影响远小于计算性能提升带来的加速比提升。此外，对于大数据体量，CPU串行执行的指令数量非常庞大，计算时间不断增加，而GPU可以同时支持多达896个线程执行，这使得GPU 在相同指令数量的情况下能够更快地完成计算任务。

表2 不同体量的数据在GPU和CPU中的加速时间Table 2 Acceleration times for different volumes of data on GPU and CPU

4.2 GPU加速处理航空γ能谱数据效果分析

不同小波基函数的滤波序列长度不一样，分解一层后的数据大小与上一层数据大小的关系为：

因此，不同长度的滤波序列下，下一层的数据大小不一，时间损耗的结果也会有差异，这种差异会对GPU加速比产生一定的影响。此处数据选用1312×101 大小的航空γ 能谱原始数据集，采用4 层分解。为测试不同小波基函数对航空γ 能谱数据的适用性，这里选取45种不同小波基函数分别对实测数据进行改进阈值降噪加速测试，不同小波基下改进阈值降噪函数的GPU加速比如表3所示。

表3 不同小波基函数下改进阈值降噪函数的GPU加速时间比Table 3 GPU acceleration time ratios for improved threshold denoising function with different wavelet basis functions

由表3可知，随着滤波序列长度的增长，各函数的加速比也在不断增加。其中，coif5小波基函数的加速比最好，高达350倍左右，且对改进阈值降噪函数达到了570 倍左右的加速比，而coif5 小波基函数也是所有小波中滤波器序列最长的小波，滤波器序列长度为30，说明不同滤波器序列对加速的效果也有较大的影响。同时，相比于总数据插值后的数据体量或者在核应急和辐射环境评估中低空无人机飞行的数据，此处所选用的1312×101 数据体量较小。数据实验表明，在GPU算力范围内，数据体量越大，加速比越高。因此，可以认为，将该算法应用于大区域航空γ能谱数据中加速比将更为可观。

5 阈值降噪效果实例验证

实例选取AGS-863型航空γ能谱仪在内蒙古锡林郭勒盟某矿区测量得到的试验飞行数据。为了更好地体现降噪效果，此处人工添加均值为0、标准差为1的白噪声，噪声分布如图7所示。

随后计算降噪后信噪比，通过信噪比来表征不同方法的降噪效果，其中，降噪前的信噪比计算公式为：

式中：xi，j是原始信号值；yi，j是染噪后的信号值，降噪后的信噪比计算公式为：

式中：zi，j是降噪后的信号值。通过公式计算，发现bior2.4 小波基函数降噪后的信噪比值最小，降噪后效果如图8 所示，bior3.1 小波基函数降噪后的信噪比值最大，其降噪效果如图9所示。

图8 bior2.4小波基函数的三种阈值降噪的降噪效果(a) 原始数据，(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小波基函数的三种阈值降噪的降噪效果(a) 原始数据，(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 小波基函数对测区西南方的噪声和条带有所消除，软阈值的消除效果最好，但同时丢失了部分异常，硬阈值和改进阈值在消除条带的同时，保留了一些基本地质特征，由图8 可知，三种阈值降噪对条带消除的效果均较好，但对于测区北部数据降噪处理的结果并不理想，基本与原始图像一致；图9 中，使用bior3.1 小波基函数，出现了过度降噪导致的图像失真现象，并且引入了行方向条带，硬阈值在测区西南方还引入了大面积原始数据中不存在的噪声。因此，这里选取降噪后信噪比在中间位置的小波基滤波器序列，其降噪效果如图10所示。

图10 不同小波基函数下各阈值降噪效果(a) 原始数据图，(b) 基于小波基函数bior3.7的改进阈值降噪效果，(c) 基于小波基函数coif1的软阈值降噪，(d) 基于小波基函数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 中，三种阈值函数在东南方主条带处的降噪效果均不佳，通过与实测区域地质图对比，发现此处是其他测量误差导致的假异常，且此处软阈值降噪出现失真，反而衍生出条带；其次，在西南方位置处的大部分条带消除效果很好，但是硬阈值相较于其他两种函数，图像更加粗糙；与此同时，在测区西北部三种函数对条带的消除均取得了显著的效果，但在特征信息的保留上，改进阈值效果更佳。综上所述，改进阈值降噪对于条带噪声的消除情况比其他两种函数好。

6 结语

针对航空γ能谱数据体量巨大而CPU执行数据后处理时计算效率低的问题，采用了基于GPU的二维离散小波阈值并行降噪技术对航空γ能谱数据进行处理。首先进行了GPU加速效果的测试，结果显示，不同的block 尺寸对GPU 计算时间产生了显著影响；其次测试了不同数据体量下的加速效果，结果发现，随着数据体量的增大，相较于CPU，GPU的加速比不断增加，但在数据体量为256×256 时加速比相对于128×128 和512×512 较低，这是由于数据传输所带来的影响；为了更充分利用GPU 性能，建议尽可能增大数据体量。在航空γ能谱数据处理测试中，不同小波基函数对加速比亦产生了显著影响。除了部分滤波器序列较短的小波基函数外，80%的小波基函数总时间加速比都达到了100 以上，91%的小波基函数总时间加速比达到90倍以上，特别是coif5 小波基函数的加速比达到350 倍，对于阈值降噪函数的加速比接近570 倍；最后对不同降噪函数的处理结果进行了对比，结果显示：所有函数都存在信噪比较低时降噪不足和信噪比较高时过度降噪的情况；使用硬阈值选取小波基函数coif5、软阈值选取小波基函数coif1 以及改进阈值选取小波基函数bior3.7 进行处理，都取得了显著的降噪效果。其中改进阈值的效果最佳，但在异常条带部分的消除效果相对较差，后续需要进一步研究改进。

作者贡献声明熊超负责研究的提出及设计、数据的收集和整理；王欣负责文章的起草和最终版本的修订；王鑫杰负责程序设计及实验数据的处理；吴和喜负责最终版本的修订和项目的监督及管理。