前兆观测干扰信号频谱特征分析

李继业 武晓军 刘长生 任建辉 高研孙鹏宇 孙强 秦丽岩 战玉明

1）黑龙江省地震局，哈尔滨市南岗区鸿翔路24号 150090

2）哈尔滨市防震减灾技术中心，哈尔滨 150021

0 引言

在数字信号中，对随机信号的处理有重要意义。随机信号又可分为平稳随机信号和非平稳随机信号，而现实中大多为非平稳随机信号，所谓非平稳随机信号，就是其统计特性是时间的函数。由于计算机技术的飞速发展与普及，非平稳随机信号分析与处理的应用前景广阔（Cohen，1989）。

过往的地震前兆研究表明，在观测数据中一般都包含着动因复杂的各种信息，可将其分为正常动态及异常动态，而异常动态中又有与地震活动有关的震兆异常和与地震活动无关的干扰信号。地震前兆观测数据由于不同程度地受到观测站点的场地、自然环境、人为活动等各种因素的制约与干扰，其干扰信号出现的频次往往比前兆异常信息要多得多。若干扰信号幅度大于正常观测背景值时，地震活动前兆异常信息就难以提取。因此，如何识别和排除干扰，是提高地震监测研究水平的一项非常重要的基础性工作（邱永平，2011）。

前兆观测数据可视为不同频率、周期成分的叠加。数字化观测中各种干扰信号与震兆异常信息混杂，有些干扰信号无论是形态上还是变幅上，都极易与震兆异常信息混淆，从而给数字化数据的使用和地震预测研究带来很大的困难。本文选用基于傅里叶变换的傅里叶谱、功率谱以及基于希尔伯特-黄变换的边际谱方法，将不同周期成分分离，形成不同频率随时间变化的多个时间系列，分析频谱特征，进而提取黑龙江数字化水位和竖直摆倾斜测项常见干扰类型，初步总结了水位高频干扰、周期干扰以及竖直摆倾斜的风扰、同震应变等常见干扰信号的频谱特征，以期探索出用谱分析方法提取前兆信息的有效途径。

1 方法简述

1.1 基于希尔伯特-黄变换的边际谱方法

希尔伯特-黄变换（Hilbert-Huang Transform，简称 HHT）由 Huang等（1998、1999、2003）提出，是一种运用于非线性、非平稳信号的处理方法，是对傅里叶变换基本信号和频率定义做了创造性改进，其不再认为信号是由一组固定频率和振幅的正弦信号的加权和，而是由固有模态函数（intrinsic mode function，简称IMF）信号所组成。希尔伯特-黄变换方法由经验模态分解（empiricalmode decomposition，简写为 EMD）（Rilling et al，2003、2007）方法和 Hilbert变换组成。EMD方法即Huang变换，它依据信号本身时间尺度特征，将信号分解为含有不同时间尺度且满足以下2个定义条件的一组IMF，每个IMF可以被认为是信号中固有的一个模态函数：①对于一列数据极值点和过零点数目必须相等或至多相差一点；②在任意点，由局部极大点和极小点构成的2条包络线的平均值为零。

EMD分解过程：①找出信号x（t）的所有极大点和极小点，将其用3次样条函数分别拟合为原数据序列的上、下包络线，上、下包络线的均值为平均包络线m1，将原数据序列减去m1可得到一个去掉低频的新数据序列h1。一般h1不是一个平稳数据序列，为此重复以上过程n次，使所得到的平均包络线趋于零，此时的h1n就是第一个IMF（c1），它表示信号数据中的最高频成分。②用x（t）减去c1得到一个去掉高频成分的新数据序列，重复步骤①，得到一系列cn和最后一个不可分解的序列rn，rn代表x（t）的均值或趋势项。那么原序列x（t）可表示为IMF分量和一个残余项的和

信号经分解后得到多个IMF分量，对每个IMF分量进行Hilbert变换，即可得到每个IMF分量的瞬时频率，综合所有IMF分量的瞬时频谱得到Hilbert谱。

获得IMF分量以后，对每一阶IMF作 Hilbert变换，设为y（t）

x（t）和y（t）共同组合为一解析信号z（t），采用极坐标形式表示

定义瞬时频率

由式（6）可看出，ω（t）是时间t的单值函数，即某一时间对应某一频率，每个IMF序列在每一点的频率唯一。对每一阶IMF作Hilbert变换，并求出相应解析函数的幅值谱和瞬时频率，从而原始信号可以表示为

上式的频率ωj（t）和幅值aj（t）都是时间的函数，可构成幅值、频率、时间的3维时频图即Hilbert幅值谱，简称Hilbert谱，记为

基于Hilbert谱公式（8），通过对时间的积分，可以定义边际谱h（ω）为

边际谱表达了每个频率在全局上的幅度，它代表了在统计意义上的全部累加幅度。

1.2 基于傅里叶变换的傅里叶谱和功率谱

快速傅里叶变换（Fast Fourier Transform，简称FFT），是离散傅里叶变换的快速算法，也可用于计算离散傅里叶变换的逆变换。它是将时间域的信号变换成频率域的信号，得到与原信号长度相同的时间序列（万永革等，2009），因为有些信号在时域上很难看出什么特征，但当变换到频域之后，就完全不同了。FFT变换还可以将一个信号的频谱提取出来，这常用于频谱分析。由于随机信号是一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能应用确定信号的频谱计算方法去分析随机信号的频谱。这样就可以用功率谱来描述随机信号的频域特性，即信号能量特征随频率的变化。为行文简洁，傅里叶变换的具体算法、公式不再赘述。

2 干扰频谱特征

过往的地下水微动态和地倾斜研究表明，井水位和钻孔倾斜对地壳应力-应变的响应频率较宽，可记录到地震波中的瑞雷波（10～20s）、潮汐波动（⅟3～⅟2 d）、气压波动（几小时～几天）、风扰波动（十几～几十分钟）。除此之外，井水位和地倾斜观测还受降雨效应（补给、荷载）、地下水开采和观测仪器故障等干扰因素的影响（张淑亮等，2005；薄万举，2010）。因此，排除干扰信号是前兆观测和数据分析的重要前提。由于篇幅关系，本文仅例举具有普遍性、代表性的水位高频、周期性干扰以及竖直摆倾斜的风扰、同震应变等干扰信号的频谱分析结果。

在实际观测中，甘南台水位存在明显的周期性干扰信号。通过对甘南台水位2013年2月17日0时1分～8时59分的540个分钟值干扰信号去趋势（图1（a））后进行频谱特征分析发现，傅里叶谱和功率谱周期存在2.2～10min多个频谱值（图1（c）、（d），表1），而边际谱周期存在9～18min多个频谱值（图1（b）、表1）。傅里叶谱和功率谱谱峰分布在整个频域区间，并有多条等间距的谱峰出现，而边际谱却不明显，说明影响甘南水位的信号应该是一个固定周期的信号。可见，傅里叶谱和功率谱识别水位周期干扰信号效果更好，体现了傅立叶变换在平稳信号处理中的优势。

图1 典型的水位周期性干扰频谱特征

图2 典型的水位高频干扰信号的频谱特征

表1 干扰信号的频谱特征

在实际观测中，孙吴台水位存在明显的漏电高频干扰信号。通过对孙吴台水位2009年10月1日0时1分～8时59分的540个分钟值干扰信号去趋势（图2（a））后进行频谱特征分析发现，傅里叶谱和功率谱周期存在2～20min多个频谱值（图2（c）、（d），表1），而边际谱周期存在13～20min多个频谱值（图2（b）、表1）。傅里叶谱和功率谱谱峰分布在整个频域区间，而边际谱却不明显，说明影响孙吴台水位的高频干扰信号亦是一个有固定周期的信号，属于平稳信号。可见，傅里叶谱和功率谱识别水位高频干扰信号效果更好，体现了傅立叶变换在平稳信号处理中的优势。

形变观测受多种因素影响，包括气象变化、人为影响、仪器故障等因素，使排除干扰成为识别异常必不可少的工作。在目前认识水平下，干扰形态与地震前兆形态极为相似，给干扰与前兆的判断带来很大困难。通河台竖直摆倾斜存在明显的季节性大风干扰，通过对通河台竖直摆倾斜2012年8月29日14时1分～23时59分NS、EW2个分量各600个分钟值干扰信号去趋势（图3a、4a）后进行频谱特征分析发现，傅里叶谱和功率谱对风扰信号频谱特征提取效果较好，而边际谱提取效果并不明显。通河台竖直摆倾斜NS测项傅里叶谱和功率谱周期存在 11～26min多个频谱值（图 3（c）、（d），表 1），而边际谱频谱周期存在 8～27min多个谱值（图3（b）、表1）。EW测项傅里叶谱和功率谱周期存在 3.7～18min多个频谱值（图 4（c）、（d），表 1），边际谱周期存在 4～26min多个频谱值（图 4（b）、表 1）。功率谱和傅里叶谱出现的频谱周期均小于27min，边际谱存在77min的独立周期。可见，傅里叶谱和功率谱识别倾斜风扰信号效果更好，体现了傅立叶变换在平稳信号处理中的优势。

图3 典型的竖直摆倾斜风扰信号的频谱特征（NS）

图4 典型的竖直摆倾斜风扰信号的频谱特征（EW）

在倾斜仪记录中，由于地震引起地表变形使观测曲线偏离震前的位置，这种突变称倾斜同震应变。当地震发生时，孕震体释放能量，震源区应变场发生调整，由于地震动导致的应力-应变传递，使得各地应变场也随之调整，在震源区外的仪器能够记录到这种直接调整，这就是为什么在上千千米范围内的仪器都能观测到同震应变的原因（陈德福，1993）。前兆观测仪器记录到的同震应变对于正常背景观测数据而言也可以看作是一种干扰信号。

2013年 4月 6日 12时 42分，印度尼西亚（138.5°E、3.5°S）发生M7.0强震，震源深度70km，在震中距5560km的通河台竖直摆倾斜仪记录到了明显的同震应变波。通过对通河竖直摆倾斜2013年4月6日12时1分～16时59分的NS、EW2个分量各300个分钟值同震应变信号去趋势（图5（a）、6（a））后进行频谱特征分析发现，傅里叶谱和功率谱识别同震应变的频谱特征较为一致，其中NS测项傅里叶谱和功率谱周期存在2.9～12min多个频谱值（图 5（c）、（d），表 1），而边际谱周期存在 2.8～14min多个频谱值（图 5（b）、表 1）。EW测项傅里叶谱和功率谱周期存在 2.1～9min多个频谱值（图6（c）、（d），表1），边际谱周期存在2.2～14min多个频谱值（图6（b），表1）。可见，边际谱较傅里叶谱和功率谱识别同震应变信号效果更好，体现了边际谱在非平稳信号处理中的优势。

通过对黑龙江数字化水位周期、高频干扰以及竖直摆倾斜等的风扰、同震应变干扰进行频谱分析发现，甘南台水位和孙吴台水位存在固定周期的干扰信号叠加在正常背景观测信号中，且分布在整个频域内；而通河台竖直摆倾斜风扰和同震应变干扰频谱周期主要集中在10min以下。傅里叶谱和功率谱的频谱特征基本一致，这与它们均基于傅立叶变换有关，且功率谱识别干扰信号效果更好。由于边际谱体现了某一段频率总的幅值累积特征，边际谱谱值较傅里叶谱、功率谱谱值要高出几个数量级；同时边际谱存在77min的独立周期，这可能是希尔伯特-黄变换的虚假波动现象。

图5 竖直摆倾斜同震应变的频谱特征（NS）

图6 竖直摆倾斜同震应变的频谱特征（EW）

3 讨论与结论

前兆观测及研究表明，与现有的仪器观测精度和干扰水平相比，多数的地震前兆信息属于弱信息，接近噪声水平。为了更好地、真实地、连续可靠地提取地震前兆信息，对干扰信号进行合理的提取和分析，是非常关键和重要的环节。

通过采用基于傅里叶变换的功率谱、傅里叶谱以及基于希尔伯特-黄变换的边际谱方法，对黑龙江数字化水位和竖直摆倾斜典型干扰数据进行分析，提取水位周期干扰、高频干扰以及竖直摆倾斜的风扰、同震应变等常见干扰信号的频谱特征。通过对水位、竖直摆倾斜典型干扰信号进行频谱分析，发现干扰信号优势频率在1.7×10-3×10-3Hz以上，周期一般在10min以下。其中，水位周期性干扰信号优势频率为1.72×10-3～7.61×10-3Hz，周期一般为2～10min，高频干扰信号优势频率为 1.89×10-3～7.01×10-3Hz，周期一般为 2～9min；竖直摆倾斜大风干扰信号优势频率为1.51×10-3～4.23×10-3Hz，周期一般为 4～11min，竖直摆倾斜同震应变干扰信号优势频率为2.08×10-3～7.51×10-3Hz，周期一般为2～12min。

频谱变换方法在前兆观测资料干扰信号提取中各有优势。傅里叶变换可以将任何信号分解为正弦信号的加权和，而每一个正弦信号对应着一个固定的频率和固定的幅值。希尔伯特-黄变换对于处理非线性、非平稳信号有清晰的物理意义，能够得到信号的振幅-时间-频率分布特征，能使信号分解具有唯一性，又能在时域和频域中同时具有良好的局部化性质。由于希尔伯特-黄变换突破傅里叶变换仅对平稳信号有效的不足之处，使其能够更好的应用于非线性、非平稳信号的处理。研究发现，功率谱和傅里叶谱方法对提取干扰和噪声等线性、平稳信号更加有效，而边际谱对提取同震应变等非线性、非平稳信号更加显著，丰富了频谱分析方法在前兆观测信息提取中的应用。