第1章 地球物理资料滤波方法
随着现代仪器技术的不断发展，地球物理仪器观测的精度越来越高，比如高精度重磁测量方法已经能够探测地下规模较小的异常体，诸如小矿体、局部构造等，但是它对地表不均匀性、仪器噪声等外界干扰也表现得更加敏感。为了有效地提取目标异常，避免微弱异常的损失，人们已不再满足于传统的资料处理和解释方法，而是期望通过研究新的技术方法对地球物理资料进行精细处理、解释，以获取更丰富的地质信息。本章将在分析传统滤波方法的基础上，讨论小波域滤波方法、基于高阶统计量的滤波方法、Curvelet域滤波方法及基于L2范数的滤波（空间域）方法，力争避免或减弱传统滤波方法造成的“过圆滑”，实现数据的高保真处理。
1.1 小波变换与小波域滤波
小波分析方法是20世纪80年代以来发展起来的多尺度分析方法，“小波”的概念最早是由法国地球物理学家Morlet和Grossmann在70年代分析地震资料时提出的，其后经过Meyer、Mallat及Daubechies等多位数学家的大量工作，小波变换才有了系统的理论与计算方法，被广泛地研究并应用于图像与信号、地球物理、医学技术、航空航天技术、通信、计算机技术、故障监控等众多学科和相关领域，受到科研工作者的广泛关注。
1.1.1 小波变换原理
以时间t为变量的信号f的连续小波变换（continuous wavelet transform，CWT）定义为
（1.1.1）
设定
（1.1.2）
式（1.1.2）中：函数为小波函数（wavelet function），简称为小波（wavelet），它是由函数经过不同的时间尺度（t）伸缩和平移得到的；R表示实数域；是小波原型，称为母小波（mother wavelet）或基本小波（basic wavelet）；参数b表示时间平移，不同b值的小波沿时间轴移动到不同的位置；参数a表示时间轴的尺度伸缩，大的a值对应于小的尺度，相应的小波伸展较宽，小的a值对应的小波在时间轴上受到压缩；系数表示归一化因子，它的引入是为了使不同尺度的小波保持相等的能量。
综合式（1.1.1）和式（1.1.2），得到连续小波变换的简化定义式：
（1.1.3）
即信号关于的连续小波变换可以表达为与小波的内积，连续小波变换定量地表示了信号与小波函数系中每个小波的相关或接近程度。为方便表示，将(a，b)简记为Wf?(a，b)。
构造的母小波函数必须满足允许条件：
（1.1.4）
式中：为的傅里叶变换，如果是一个合格的窗函数，则是连续函数。
因此，允许条件意味着
（1.1.5）
其物理意义是为一个振幅衰减很快的“波”，“小波”的概念即因此而来。
令母小波是中心为、有效宽度为Dt的偶对称函数，经过伸缩平移后的小波的中心为，宽度为aDt，如图1.1.1（a）所示。如果把小波看成宽度随a位置随b变动的时域窗，那么连续小波变换可被看作连续变化的一组短时傅里叶变换的集合，这些短时傅里叶变换对不同的信号频率使用了不同宽度的窗函数。具体来说，高频使用窄时域窗，低频使用宽时域窗。这被称为小波变换的“变焦距”性质，即多尺度多分辨时频局部化特性。
图1.1.1 母小波与小波及其频率特性
另外，在频域中可以观察小波变换的性质，的傅里叶变换为
（1.1.6）
若母小波的傅里叶变换是中心频率为、宽度为Dω的带通函数，那么其傅里叶变换是中心为/a、宽度为Dω?/a的带通函数，如图1.1.1（b）所示。根据帕塞瓦尔恒等式，由式（1.1.3）可以得到
（1.1.7）
因此，连续小波变换给出了信号频谱在频域窗或内的局部信息。
设＞0，a为正实变量，可以把/a看作频率变量。的带宽与中心频率之比，即相对带宽，与尺度参数a或中心频率的位置无关，这就是“恒Q性质”。把看作频率变量后，“时间-尺度”平面等效于“时间-频率”平面。因此，连续小波变换的“时间-频率”定位能力和分辨率也能够用“时间-尺度”平面上的矩形窗口来描述，该窗口的范围是矩形窗口长为[即的有效宽度]，宽为[即的有效宽度]，面积为。面积大小与a无关，仅取决于的选择，因此，一旦母小波确定，窗口的面积也就随之确定。分析的矩形窗口宽度决定时间分辨率和时间定位能力。a越小（对应越高的频率），时间分辨率越高。因此，高频分析应采用窄的窗口。由于分析的矩形窗口面积恒定，当矩形窗口变窄时，其高度不可避免地增加，会降低频域分辨率和频率定位能力。图1.1.2就是从另一角度观察到的连续小波变换的变焦距性质，它与图1.1.1形成鲜明对照。
图1.1.2 小波变换的频域性质
分析的矩形窗口宽度随频率升高（尺度减小）而变窄
由连续小波变换重建原信号，其逆变换公式为
（1.1.8）
1.1.2 小波域滤波实现
实际采集的重磁异常中常含有较高频率的白噪声干扰，只有通过滤波处理，压制噪声干扰，才能有效地突出原信号中的地质信息。下面先介绍白噪声的几个特点。
（1）白噪声可以看作平稳的随机信号，记为，它在各测点处的值是一个随机量，取值大小与其他测点处的随机取值无关。因此，白噪声的随机性表明，不同的白噪声和之间是不相关的。
（2）白噪声可以看作能量无限且零均值的。白噪声在时间域中没有衰减性，因此它具有能量无限性。白噪声是随机变化的，因此有
（1.1.9）
（3）对于确定信号，白噪声的时间域特征是均匀密集的。
（4）对于有用信号，随机干扰具备较高的频率。
重磁勘探中的异常通常表现为中低频信号特征或一些比较平稳的信号，而噪声则表现为高频信号特征。对含有噪声的位场数据做小波分解后，噪声能量主要出现在小波分析的小尺度上，且主要集中于各个尺度分量的高频部分。在高频系数方面，随着尺度变大和分解层次增多，其幅值大约按2-1/2倍快速衰减。另外，噪声在不同尺度上的特征也是不相关的。
根据上面的分析，传统基于小波分析的滤波方法大致可以归结为两类。第一类是强制的切除法滤波，即把小波分解各尺度分量或某几个尺度的高频系数全部赋值为零，然后进行小波反变换至空间域。这类方法比较简单，重构后的信号也较平滑，但是这种“一刀切”的生硬处理方式，容易丢失原信号中的有用信息。第二类是通过阈值衰减方法压制噪声。这类方法是根据经验或某种法则设定阈值，对小波分解的系数进行阈值收缩，这符合噪声在高频部分均匀密集的特点。
在实际处理应用中，阈值方法是使用最为广泛的小波域滤波算法。常见的几种阈值计算方法如下。
（1）统一阈值法。对小波分解的各层采用单一的阈值：
（1.1.10）
式中：为阈值；为噪声强度；N为样点数。
（2）定义小波阈值为
（1.1.11）
式中：M为分解的层次；为各层所取的阈值；为噪声强度。
（3）定义小波阈值为
（1.1.12）
式中：为反映噪声的小波变换模在不同尺度j上的传播因子。
（4）将小波分解第一层的小波系数全部置为零，相当于将阈值t设为第一层小波系数模的最大值。
（5）非线性阈值法。硬阈值方法与软阈值方法是两种传统的阈值处理方法。硬阈值方法准则是保留大于阈值的系数，而将小于阈值的系数都赋值为零，具体表达式为
（1.1.13）
式中：Wf?(a，b)为原小波系数；为阈值处理后的小波系数。
软阈值方法准则是将小于阈值的系数都赋值为零，而将大于阈值的系数减去一个阈值大小的量，具体表达式为
（1.1.14）
由图1.1.3（a）、（b）可以看出：硬阈值方法是使大于阈值的系数不变，而小于阈值的系数取值为零，此方法存在间断点，会造成重构失真，滤波不彻底；软阈值方法虽然没有间断点，但是原系数减去了一个阈值大小的量，也会使重构失真。
横坐标为处理前的系数，纵坐标为处理后的系数
为了克服硬阈值方法和软阈值方法的不足，使信号变化尽可能平缓，可采用非线性阈值，如图1.1.3（c）所示。非线性阈值函数表达式为
（1.1.15）
式中：为在由1衰减到0的光滑的权函数。
1.2 高阶统计量滤波
早在20世纪50年代，一些学者就开始了高阶矩的研究。Rosenblatt等（1965）发表了关于双谱估计的文章，同年，Brillinger（1965）全面介绍了多谱理论。但是，直到80年代后期，这方面的研究才真正得到迅速发展与应用，迎来了高阶谱理论和应用研究的高潮。目前，其应用范围已涉及通信、生物医学、故障诊断、声呐等多个领域，近年来也已经逐步渗透到地球物理勘探领域。
所谓高阶统计量，通常是指高阶矩、高阶累积量及它们的谱——高阶矩谱和高阶累积量谱4种形式的统计量，此外，还有倒高阶累积量谱等。通过特征函数可以引出高阶统计量的定义，并推导它们的性质。本节主要讨论高阶累积量及高阶累积量谱（简称高阶谱），其他相关理论可参考Arivazhagan等（2006）和张贤达（1996，1995），此处不做详述。
1.2.1 高阶累积量及高阶累积量谱的概念
常见的k阶累积量，用来表示高阶累积量。
设{x(n)}为零均值的k阶平稳随机过程，其k阶累积量定义为
（1.2.1）
设{x(n)}为零均值的k阶平稳随机过程，且序列的k阶累积量是绝对可和的，即
（1.2.2）
则k阶累积量谱定义为k阶累积量的k-1维傅里叶变换，即
（1.2.3）
最常用的高阶累积量谱是三阶累积量谱（简称三阶谱）和四阶累积量谱（简称四阶谱），定义形式如下。
三阶累积量谱：
（1.2.4）
四阶累积量谱：
（1.2.5）
由于三阶累积量谱只有两个变量，四阶累积量谱只有三个变量，为了方便，特别称它们为双谱和三谱。
式（1.2.3）中的时，就是常见的功率谱：
（1.2.6）
高阶累积量具有如下几种性质。
性质1：相互独立的两随机序列的组合序列的累积量等于零。
设随机序列由相互独立的两随机序列与组成，则有
（1.2.7）
因此，对于一个由具有相同分布的相互独立的随机变量构成的随机序列，它的累积量为函数。
性质2：任何高斯过程的高阶累积量均等于零。
