1. 引言
Hilbert Huang Transform (HHT)是一种全新信号处理方法，对于处理非线性、非平稳信号有清晰的物理意义，能够得到信号包含时间、频率、振幅的三维离散时频谱，提供了清晰的局部细节的时频特征(Huang et al., 1998)。HHT方法具有较好的客观性与自适应性，适合于描述具有非线性非平稳变化特征的信号，该方法自推出以来已经成功地应用在地球物理学领域、生物医学领域和设备故障诊断领域 [1][2][3][4][5]等。国内学者在研究HHT在地震资料分析方面，也发表了大量的文章。在地震信号与资料的分析和处理方面、提高地震资料分辨率、噪声剔除以及提取地震瞬时属性等方面都有较好的应用 [6]- [13]。周挚等基于HHT提取昆明、下关重力固体潮的地震前兆信息，设计了重力固体潮地震前兆分析的特征参数，计算得到的实际参数和理论参数的频率比值 [14]；李东生等在HHT方法在山丹地电场数据处理中的应用 [15]清楚地分离出大地电场的主要周期成分；王宁等在时频分析方法在形变数据中的应用研究中认为，HHT方法保存了信号本身的物理特性，同时也具有很好的时间和频率分辨率 [16]，然而HHT方法在地震前兆观测数据分析中尚处于起步阶段，还没有得到普遍应用。
HHT方法在延寿钻孔倾斜分析尚属首次尝试，利用HHT时频分析方法提取固体潮汐特征，针对形变资料中可能存在的冰裂、突跳、库容变化和大气环境等干扰情况，进行EMD分解，提取重构数据再进行频谱分析，探讨可能存在的干扰因素在时频域的优势频率。
2. HHT基本原理
HHT由经验模态分解(Empirical Mode Decomposition, EMD)和Hilbert谱分析两部分组成。在应用过程中，首先将原始信号分解成一组本征模态函数(IMF)，然后对本征模态函数进行Hilbert变换。
2.1. 经验模态分解
HHT处理非平稳信号的基本过程为：① 利用EMD方法，将给定信号分解为若干固有模态函数(以Intrinsic Mode Function或IMF表示)，这些IMF是满足一定条件的分量；② 对每一个IMF进行Hilbert变换，得到相应Hilbert谱，即将每个IMF表示在联合的时频域中；③ 汇总IMF的Hilbert谱，得到原始信号的Hilbert谱 [17]。
2.2. Hilbert变换
Hilbert变换是信号检测理论研究和通信工程应用中的一个重要工具。通过Hilbert变换，可以在不造成信息损失的前提下，将一个实信号构造成一个复信号(解析信号)，使研究实信号的瞬时包络、瞬时相位和瞬时频率成为可能。将Hilbert变换应用于IMF，使之生成解析信号，然后计算出瞬时频率，进而求出Hilbert时频谱(亦称Hilbert幅值谱)和Hilbert边际谱 [18]。
3. HHT方法在延寿钻孔倾斜数据中的时频特征分析
延寿台倾斜观测使用的是CZB-1型竖直摆倾斜仪，仪器核心部分是采用一个竖直方向悬挂的重力摆。工作原理是当地表产生倾斜变化时，摆的支撑架随着产生倾斜，而高精度的位移传感器把二者得相对位移换成电压的变化，再用记录仪器把电压信号记录下来。钻孔井深71.6 m，孔径152 mm，井口距传感器距离61 m，一定程度上避开了地表干扰。但人为干扰、仪器零漂、周围环境和天气状况等干扰因素是无法避免的。
3.1. 基于HHT方法提取固体潮信号
选取延寿短时及连续一个月的倾斜整点值数据，基于HHT方法对固体潮汐信号进行提取。计算结果表明，原始数据与理论固体潮数据分析得到的频谱信号一致，均为22.64 µHz和10.57 µHz，即周期12.27小时的半日波信息和26.27小时全日波信息(图1，图2，图3)。
对2008年5月延寿钻孔倾斜EW向整点值数据与理论固体潮数据进行分析，由EMD (图2)的分解可以看出数据被分解成7个IMF分量。频率高的IMF总是比频率低的IMF优先被筛分出来。然而，在这些本征模态函数中并非都含有地球物理信息，这就需要对这些筛分出来的IMF进行希尔伯特变换分析。图3是计算得到的Hilbert边际谱。通过分析看出imf1、imf2为固体潮汐信号，Hilbert时频谱可以很明显的反映出固体潮频率随时间变化的情况及非潮汐信号对固体潮汐的影响。边际谱分析可以很好地识别数
(a) 理论固体潮 (b) 钻孔倾斜东西分量
Figure 1. Theoretical earth tide and Hilbert spectrum of Yanshou borehole inclined minute value during 18-19 February 2008
图1. 2008年2月18~19日理论固体潮与延寿钻孔倾斜分钟值Hilbert谱
(a) 延寿钻孔倾斜EMD分解 (b) 固体潮EMD分解
(c) 延寿钻孔倾斜边际谱 (d) 理论固体潮边际谱
Figure 2. Analysis of EW integer value and theoretical earth tidal spectrum of Yanshou borehole tilt in May 2008
图2. 延寿钻孔倾斜2008年5月EW向整点值与理论固体潮谱分析
Figure 3. Contrast diagram of original curve, theoretical curve and reconfiguration curve of Yanshou borehole tilt
图3. 延寿钻孔倾斜原始曲线、理论曲线和重构曲线对比图
据中的频率信息，对2008年5月延寿钻孔倾斜EW向整点值数据分析得到的频主要信号分量为22.37 µHz和11.51 µHz，即周期12.41小时和24.13小时半日波和日波信息。根据EMD分解结果，我们可以确定imf1为半日波，imf2为日波，imf3、imf4、imf5、imf6为其他周期信息，imf7可能为仪器零漂造成的趋势项。为了更好的体现HHT对固体潮信息的提取，基于imf1和imf2分量对数据进行重构，可以看出重构数据(图4)基本上和理论固体潮信息非常吻合。
(a) 倾斜imf1-3边际谱 (b)气压imf1-3边际谱
(c) 倾斜imf4-8边际谱 (d)气压imf4-8边际谱
Figure 4. Yanshou borehole tilt and barometric pressure marginal spectrum during 7-8 November 2009
图4. 2009年11月7~8日延寿倾斜与气压边际谱
3.2. 气压干扰提取
2009年11月7~8日17个小时内气压出现剧烈变化，变幅可达22.2 hpa。对气压分钟值数据与倾斜EW向分钟值数据分别做EMD分解，imf分解为9个分量，1~3 imf分量为高频信息，4~8 imf分量为低频信息，重构4~8 imf分量。分别对1~3 imf分量和4~8 imf分量进行谱分析(图4)，倾斜4~8 imf分量重构数据保留了完整的半日波周期数据，气压4~8 imf分量重构数据很好地体现了全日波周期数据。倾斜1~3 imf分量重构边际谱干扰信息主要体现在68.51~265.4 μHz (周期63~243分钟)，气压1~3 imf分量重构边际谱干扰信息主要体现在477.4 μHz (周期在35分钟)左右，因此气压不是延寿钻孔倾斜的主要干扰源。
3.3. 冰裂干扰
延寿台西南方向为长寿湖水库，每到冬季和春季由于水库封冻与融化会存在冰裂现象。对2008年4月13~14日钻孔倾斜数据进行EMD分解，将IMF5-9分量进行数据重构并进行边际谱分析，IMF5-9分量半日波及全日波信息优势明显，但IMF5分量仍留存冰裂干扰频谱，IMF6-9分量重构数据更为理想(图5)。
3.4. 库容干扰分析
延寿台西南方向为长寿湖水库，水库距钻孔倾斜观测井100 m左右，库容最高可达2500万m3，根据邱泽华老师提出的三维荷载模型 [19]，已计算得出水库库容日常变化对延寿钻孔倾斜的干扰小于限定标准 [20](图6)。
对2008年~2009年水库库容和延寿钻孔倾斜观测数据进行频谱分析(图7)，水库库容变化边际谱优势频率为0.115~0.1365 µHz (85~101天)，倾斜边际谱优势频率为0.2765~0.4066 µHz (28~41天)。库容变化并不是延寿钻孔倾斜的主要影响因素，但在月周期作用的频率成分与延寿钻孔倾斜频谱有所交叉。
4. 结论与讨论
1) 时频分析时实际观测数据经常包含较大噪声和各种干扰，有用的信号往往偏弱，这样很难直接利用HHT方法提取出我们需要的前兆信息，因此，非常必要对数据进行预处理，将明确原因的干扰先行剔除，并且
Figure 5. EMD decomposition and refactoring data of borehole tilt during 13-14 April 2008
图5. 2008年4月13~14日钻孔倾斜EMD分解及重构数据
Figure 6. The curve diagram of Reservoir Capacity and Yanshou borehole tilt during 2008-2009
图6. 2008~2009年水库库容与延寿钻孔倾斜曲线图
(a) 2008年库容边际谱 (b) 2009年库容边际谱
(c) 2008年倾斜边际谱 (d) 2009年倾斜边际谱
Figure 7. The Marginal Spectrum of Reservoir Capacity and Yanshou borehole tilt during 2008-2009
图7. 2008~2009年水库库容与延寿钻孔倾斜边际谱
结合数字滤波等其它有效方法联合分析，先将观测数据进行滤波，对滤出的某些频段信号再进行EMD分解。
2) HHT方法能够很好地提取固体潮信息，延寿钻孔倾斜原始数据与理论固体潮数据分析得到的频谱信号一致，均为22.64 μHZ和10.57 μHz，即周期12.27小时的半日波信息和26.27小时全日波信息。气压优势频率主要分布在477.4 μHz (周期在35分钟)左右。水库库容变化边际谱优势频率主要体现为长周期：0.115~0.1365 µHz (85~101天)。
3) 固体潮特征具体较为明确的物理意义，HHT方法也可以很好剔除部分干扰，提高数据分析质量。
基金项目
中国地震局监测、预测、科研三结合课题项目资助，课题编号：CEA-JC/3JH-160802。