在地震资料分析中，地震信号的频率时常是随时间和距离变化而变化的，震源发出的地震信号在地层中穿越而被吸收衰减，大地对地震信号进行连续的滤波作用使得被检波器接收的信号时频信息受到严重的影响，高频信息被大量衰减，波形响应时延.这些对于识别有效信号和干扰信号都是非常不利的，然而通过时频分析，掌握其时频特性后，这样的变化对研究地下地质情况又是十分有帮助的(兰金涛，2011; 杨兆斌等，2011; 李玲利等，2012).通过这些变化能够反映出地下构造，岩性等方面的变化.

面波是地震波沿地表传播的低频、低速、高频散、强振幅波动，由于其低速和强振幅特点在炮集记录上通常呈扫帚状(夏洪瑞等，2013).这样在原始的炮集记录上，有效波往往被面波掩盖，导致在后续的处理中带来一定的干扰，对剖面的解释就会造成误判，从而难以得到正确的解释结果而前功尽弃.学者们通过研究面波的传播特性以便更好地去找到利用和压制它的方法(李媛媛，2004; 李晶，2006; 石颖和王维红，2012).面波与有效波之间相互混杂使得面波的压制变得十分困难，要想把面波压制的同时对有效波的损害达到最小是所有滤波方法都追求的终极目标(董烈乾等，2011; 聂鹏飞等，2012).

面对庞大的地震数据，借助于计算机来处理这些数据变得十分必要，然而如此海量的数据要实现高精度而快速的处理却是高性能计算机的一大挑战，如何能够获得快速而精准的答案是每个地球物理学者所探索的问题.GPU与生俱来的特性为它在地震勘探方面的应用提供了可能，学者纷纷尝试将算法移植到GPU上(李博等，2009; 李肯立等，2009).虽然GPU计算开始于近几年，但它迅猛的发展已展示其通用计算领域的成功，一些算法在适合GPU移植时动辄效率达几十倍的提高，更甚者达百倍之多(张广智和陈雷，2011; 钟勇和陈磊，2011; 陈召曦等，2012; 王保利等，2012; 刘守伟等，2013; 刘文卿等，2013). 1 GPGPU简介及CUDA 1.1 GPGPU简介

GPGPU全称General Purpose Graphics Processing Unit，即通用计算图形处理器，是一种可编程图形处理器.可编程显卡的出现为图形硬件在通用计算方面提供了可能，它按照所编写的代码来执行图形任务，输入的数据为图形的形状描述和属性参数，以向量和矩阵方式输入，经过一系列的向量和矩阵运算来完成图像变换(仇德元，2011).由于不同的显卡其硬件接口存在差异，必然催生出图形API(Application Programming Interface).OpenGL(Open Graphics Language)和Direct3D就是图形API的典形代表，图形API为开发人员提供了与底层硬件无关的接口(仇德元，2011).图形计算时所采用向量、矩阵的方式被开发者所关注，若能将GPU出色的实时图形处理能力用于计算其它领域的问题将会有什么结果？GPU高度并行化的架构和可编程的着色器使人们渐渐尝试将它用于科学计算，此时可编程的着色器和着色语言成了核心技术，将图形数据用所要计算的科学数据替代输入GPU，通过GPU计算后的数据作结果输出，而不是向图形输出设备转移.

从2006年以后的GPU基本都将通用计算考虑在内，NVIDIA公司作为其代表在同年发布了统一着色器模型，每个着色器都能完成原来只有特定着色器才能完成的任务，这使得GPU在计算时同时有许多着色器一起在工作，效率得到了很大提高，同时统一计算设备架构语言(Compute Unified Device Architecture，简称CUDA)发布.随后由苹果公司和其它几个公司发起的开放式计算语言(Open Computing Language，简称OpenCL)也发布，这两种语言在接下来的几年中成为GPGPU语言的代表，并被广泛使用.起初ATI公司也发布了Stream作为其GPU计算标准，后来转向OpenCL，除此之外，DirectCompute、BrookGPU等也是不同的设备厂商推出的计算标准(仇德元，2011). 1.2 CUDA简介

CUDA是由NVIDIA(英伟达)公司推出的通用并行计算架构，该架构使得GPU能够解决许多复杂的计算问题.它包含了CUDA的指令集架构以及GPU内部的并行计算引擎.开发人员目前可以使用C语言和Fortran语言来为CUDA架构编写程序，C语言和Fortran语言都是应用广泛的高级编程语言，因此用户不必再学习更多关于GPU的知识就能写出能够在GPU上执行的优秀代码，同时CUDA也支持当前流行的科学计算软件Matlab.CUDA的工具箱也包含了cuBLAS、cuFFT等数学库，使得人们使用GPU运算更加便捷.CUDA中支持OpenCL，但其依赖于支持CUDA的设备，即只有计算机中有支持CUDA的显卡才能进行GPU计算，而OpenCL则在某种意义上已经超出了GPU计算的范畴，它支持各种异构资源(仇德元，2011).CUDA也支持图形的互操作(即可以将处理的数据传递给显示设备，也可以对图形数据加工处理)，它可以跟OpenGL和DirectX实现无缝连接，使得科学计算可视化得到这种桌面级超级计算机的加速，同时计算的过程和结果可以可视化监控，实时修正，直观被呈现，方便人机交互.

CUDA的线程管理由许多线程(Thread)组成线程块(Block)，然后又由许多线程块组成线程格子(Grid)，在对数据做GPU计算之前我们首先要准备计算的数据，然后根据GPU的线程结构特点将数据载入GPU的显存位置，由各线程执行核函数的各副本对数据做处理并得到结果. 2 时间S变换

S变换是对信号加上一个宽度随频率变化而变化的高斯窗，认为在窗内，函数是一个平稳信号，然后对平稳信号做傅里叶变换，因此窗函数越宽可以分辨出频率越低的信号，反之则可以分辨较高频率的信号(Stockwell et al.，1996).其连续变换为

因为高斯函数的傅式变换还是高斯函数，因此由卷积定理，S(τ，f)是可以通过信号的傅里叶变换来求，即所谓的频率S变换(Stockwell et al.，1996)：

同时也可以定义信号h(t)加窗的函数集(Schimmel and Gallart，2005)：

连续时间反S变换：

为了对时间S变换的时频特性进行探讨和与频率S变换对比，分别对时间S变换和频率S变换的重构误差、滤波边缘差异进行对比分析，并运用时间S变换做数值试验和实际地震资料处理.

图 1a是设计的一个时间采样间隔为1 ms，50 Hz单位正弦信号.图 1b是分别使用时间S变换和频率S变换的结果来重构信号，用重构信号与原始信号相减，求出重构误差.计算结果显示频率S变换重构误差为零(蓝色星形线)，而时间S变换存在一个百分之五左右的误差(红色实线)，该误差来源于对时间S变换的理论做了一个近似 .

图 1

Fig. 1

图 1 S变换的重构误差 (a)原始信号；(b)时间S变换与频率S变换重构误差.Fig. 1 The reconstruction errors of S transform (a)The original signal；(b)The reconstruction errors of time S transform and frequency S transform.

为了能清楚的看出时间S变换和频率S变换在信号滤波时的局部抖动而设计一个单位正弦信号，其时间采样间隔为2 ms，时间延续长度为1024 ms，频率为50 Hz，现在分别用频率S变换和时间S变换去对100 ms到150 ms滤波，得到的结果如图 2a和图 2b所示.

图 2

Fig. 2

图 2 S变换的边缘效应 (a)频率S变换滤波后重构信号；(b)时间S变换滤波后重构信号.Fig. 2 Side effects of S transform (a)The reconstruction signal after frequency S transform filtering；(b)The reconstruction signal after time S transform filtering.

频率S变换需要用到信号的傅里叶变换结果，滤波后重构信号时是在频率上对时间累加，然后再反傅里叶变换回时域，在反变换过程中，由于傅里叶变换对平稳信号是全局的，而非平稳信号在重构时会被以平稳信号的方式重构，由于傅式变换的一些特点，因此在局部上会出现曲线的抖动，如图 2a所示.时间S变换是用了一种近似，在重构时会产生一定的误差，但是整个计算过程中不使用傅里叶变换，因此信号局部特征不会被全局影响，能很好的得到保留，比频率S变换要好很多，如图 2b所示.

设计一个信号，时间采样间隔为1 ms，时间长度512 ms，前256 ms为一个50 Hz的单位余弦信号，后面是250 Hz 的单位余弦信号，其中41 ms到104 ms上有一 300 Hz的信号如图 3所示.

图 3

Fig. 3

图 3 混合信号Fig. 3 The hybrid signal

对上述信号做S变换得到时频谱如图 4a所示，可以看出，时间上对应的各频率已经能够区分，从高频到低频都能清楚描绘，尤其是各频率信号的起始时间和终止时间，这对精确的时频分析和滤波是非常重要的.滤掉其中41 ms到104 ms上的300 Hz信号后时频谱如图 4b所示.

图 4

Fig. 4

图 4 信号的S变换时频谱 (a)混合信号的时间S变换时频谱;(b)滤波后的时间S变换时频谱.Fig. 4 S transform spectrum map of signal (a)Time S transform spectrum map of hybrid signal;(b)Time S transform spectrum map after filtering.

将滤波后的时频谱数据反变换到时间域即可得到滤掉干扰后的信号，从图 5的结果可以看到已经滤去了目标频率信号，只在起始时间和终止时间处有一很小的突变.

图 5

Fig. 5

图 5 滤波后的信号Fig. 5 The signal after filtering

上述数值实验表明，时间S变换与频率S变换都具有很好的时频特性，虽然时间S变换在振幅精度上不如频率S变换，但是误差具有一个范围，且在滤波时的局部信号保留上有绝对的优势.下面则通过一个实际的地震资料处理来展示时间S变换的滤波效果.

图 6为面波能量很强的原始单炮记录，通过第一道数据的时间S变换，见图 7a，可以看到数据的主频在30 Hz左右，随着时间持续，高低频信号都衰减强烈，而高频更甚.

图 6

Fig. 6

图 6 原始单炮带面波数据Fig. 6 The original seismic single shot data with surface wav

图 7

Fig. 7

图 7 实际资料的S变换时频谱 (a)第一道数据时频谱(无低频面波能量);(b)第四十道数据时频谱(有低频面波能量).Fig. 7 S transform spectrum map of real data (a)The time S transform spectrum maps of the first trace(without low frequency of surface wave);(b)The time S transform spectrum maps of the fortieth trace(with low frequency of surface wave).

与第一道数据的时间S变换时频谱对比，第四十道数据的时频谱中，见图 7b，在时间的终点出现了低频的强能量信号，对比面波的时间和频率特性可判断这就是面波在时频谱上的表现，若去掉该处的信息后反变换回时域即可滤掉对应的面波.依次对每道做这样的处理即可滤掉炮集上的面波.

从实际资料处理的结果来看，通过时间S变换处理，面波已经被很好的压制，有效信号特征突出，相比原始剖面，得到的结果是一个信噪比较高的剖面，如图 8a所示.而经过S变换所滤掉的面波如图 8b所示.可见，经过S变换分离面波是合理的，在分离掉面波的同时，有效信号得到极大保护.

图 8

Fig. 8

图 8 实际资料滤波效果 (a)滤波后的单炮数据;(b)滤掉的面波.Fig. 8 The filtering effect of real data (a)The seismic single shot data after filtering;(b)The surface wave of filtered out.

实验中使用的机器是ACER 4750G，操作系统为Windows 7，CUDA版本为5.0，CPU和GPU参数如表 1.

表 1(Table 1)

表 1 机器配置参数 Table 1 The computer configuration parameters 型号 内存 内存总线宽度 主频 计算 能力 CPU Intel i3-2330M 3G — 2.20 GHz — GPU GeForce GT 540M 2G 128-bit 1344 MHz 2.1

表 1 机器配置参数 Table 1 The computer configuration parameters

表 1中的数据表明CPU的主频和内存都比GPU要好很多，根据时间S变换计算的特点，下面分别使用全局内存和零拷贝内存来处理相同数据，通过处理不同数据量所消耗的时间与相同条件下在CPU上所消耗的时间相除，得到一个加速比，通过加速比来判断加速效果.

实验表明虽然CPU的内存和主频都高于GPU，但是其计算速度却不如GPU快.通过多次测时计算的平均加速比如图 9和图 10所示，图中可以看出，随数据量的增加，加速比都快速增加，从图 9中可以看到，加速比增加到一定值后会缓慢下降，其原因与GPU的硬件结构有关，数据量太小不利于线程缓存的隐藏，数据的转移消耗的时间所占比例过大，而GPU单个计算单元的计算能力没有CPU强.数据增加时对GPU的性能利用率上升，直到达到最大值，随着数据量的继续增大，GPU的性能已经充分利用，这样的情况下，GPU的计算加速比就会缓缓下降，如图 9所示.

图 9

Fig. 9

图 9 正变换加速比随采样点变化Fig. 9 The speed-up ratio of forward transform vary with sample points

图 10

Fig. 10

图 10 反变换加速比随采样点变化Fig. 10 The speed-up ratio of inverse transform vary with sample points

从图 9和图 10中可以看到，全局内存计算的加速比曲线在零拷贝内存加速比曲线下方，也就是使用零拷贝内存有更高的加速效果，其加速曲线变化过程的原因与全局内存的原因相似.就图 10而言，反变换中零拷贝内存的加速效果更加明显，其原因与GPU的数据传输机制有关，因此对于不同的算法，我们在GPU上应采取不同的存储方式来获得最大的计算效率.

我们对野外采集的原始炮集记录做面波压制，在相同条件下GPU使用全局内存加速达到21倍，使用零拷贝内存也达到67倍，相比CPU的计算效率有很大提升，由于滤波过程中的数据预处理使得加速比相比单纯的变换有所降低. 5 结 论

5.1    S变换可分为频率S变换和时间S变换.时间S变换是一种近似，其重构误差在信号最大振幅的百分之四到五左右，相比频率S变换而言，该误差在容许范围内.但是，时间S变换在对数据局部保真方面比频率S变换要好很多，因此时间S变换具备更强的优势.

5.2     小数据量时GPU加速效果不太突出，随着计算的数据量增加，加速效果明显，这是因为小数据量不利于线程缓存的隐藏，GPU计算是不合算的.

5.3     时间S变换适合移植到GPU上执行，其加速比增加十分明显，已经是CPU运算的上百倍.而使用全局内存的反变换相对零拷贝内存的加速比明显较低，但总的来说加速效果也是惊人的.在做数据处理时对不同的数据最好使用不同的存储对比，以便使用最快的存储方式.

5.4     数据量的持续增加，由于GPU的硬件条件，加速比是会有所降低，但总的来说，加速效果是突出的.

5.5     为了能发挥GPU的高性能计算，应当在计算时让GPU一直处于“忙碌”状态，对它的显存和线程都要充分利用.