胡 玮，冯雪东，石 伟，郭 伟，刘 伟
（乌加河地震台，内蒙古自治区 巴彦淖尔市 015323）
数字形变仪器测量会受到多种因素的干扰，包括自然、人为、地球物理事件、仪器故障等因素[1-2]。由于小波分析具有能够根据分析对象自动调整有关参数的“自适应性”和能够根据观测对象自动“调焦”的特性，而广泛应用于干扰识别及噪声去除等研究领域[3]。宋治平和倪友忠等[4-5]的研究表明，小波分析可应用在地震震相识别和前兆异常分析等领域；
张军等[6]应用小波分析方法，在流体学科中较好地抑制了随机噪声；
刘建明等[7]通过运用小波分析将形变观测数据分解为不同的频段，提高了识震能力。去除形变观测中的干扰信号，对提高观测质量、识别前兆异常有一定的作用。目前对定点形变仪器的小波分析及阈值去噪的研究较少，本文应用db小波分析方法，对乌加河地震台DSQ水管倾斜仪所受的干扰进行了研究，总结分析其干扰规律和特征，并通过软阈值去噪方法剔除干扰，取得较好效果。
乌加河地震台（以下简称乌加河台）是国家级综合观测台站，地处狼山前东西向大断裂带和河套断陷盆地北缘阴山纬向构造带的中西段与狼山弧型构造带的复合部位[8]。观测学科有：形变、测震、地磁、地电和重力等。乌加河台形变山洞于1972年人工开凿，山洞洞室进深253 m，台基属花岗岩，完整性好，平均覆盖厚度189 m，观测条件较好（图1）。乌加河台DSQ水管仪于2001年进行了数字化改造。
图1 乌加河台山洞洞室平面示意图Fig.1 The cavity plan of the digital deformation instrument at Wujiahe Seismic Station
2.1 小波变换基本原理
小波变换是一种信号的时频分析方法，具有多尺度（分辨率）分析的特点，能有效区分平稳信号中叠加的非平稳部分，对数字化前兆资料的干扰识别与去除以及对不同频率的信息识别功能较强[4]。
离散信号的小波变换可表示为：
式中，f（x）为离散的信号序列，Ψa，b（x）为小波函数，a为尺度（伸缩）因子，b为平移因子。
采用Mallat算法，将数字信号f（x），分解成不同频率成分[9]：
式中，j为分解阶数，Aj f（x）是信号f（x）的频率不超过2-j的部分，称近似部分，Dj f（x）分是信号f（x）的频率介于2-j与2-j+1之间的部分，称细节部分，Aj与Dj分别为二进制离散近似与二进制离散细节，Aj与Dj满足如下关系：
式中，H为尺度函数对应的低通滤波器，G为小波函数对应的高通滤波器，J为最佳尺度。
设信号f（x）的采样率为Fs，根据采样定理，f（x）的最大频率为奈奎斯特频率fNyquist=Fs/2，记为fN，f（x）频带可表示为［0，fN］，则Mallat塔式算法可图示为：
由图2可知，与小波分解阶数j对应的近似部分Aj f（x）和细节部分Dj f（x）的频带分别为［0，fN/2j］、［fN/2j，fN/2j-1］。每一阶的近似部分Aj f（x）相当于对原始数据进行低通滤波，过滤了D1f（x）、D2f（x）…Dj f（x）等细节部分，只保留原始序列中频带为［0，fN/2j］的信号主体部分。
图2 应用Mallat塔式算法对信号进行小波多尺度分解示意图Fig.2 Schematic diagram of wavelet multi-scale decomposition of signal using Mallat tower algorithm
2.2 小波基的选择
小波变换与Fourier变换的区别之一是小波变换的小波基函数具备多样性，选择适合的小波基函数是利用小波分析识别并去除数字前兆异常信号的前提条件[9]。小波基函数选择的一般原则为：小波基函数具备紧支性、对称性、正则性（光滑性）和恰当的消失矩阶数等[10-11]。根据尹继尧、王嘉琦和孙伶俐等[12-14]的研究，应用Daubechies（dbN）小波分析方法可对形变观测数据中的噪声和干扰进行有效地识别。选取dbN小波为小波基函数可兼顾正交小波的紧支性和正则性[15]。
为选取合适的小波基，本文采用控制变量法：将分解层数均设置为8（对应的细节周期范围为2～512 min，可分解出固体潮信息），阈值均使用固定阈值，阈值去噪方法均为软阈值去噪，通过计算同一含噪信号的不同小波尺度均方根误差，确定最佳分解和重构尺度。
选取含噪信号为2020年11月25日水管仪NS向的数据，表1为计算的均方根误差结果。
表1 不同小波尺度下的均方根误差
表1中RMSE为均方根误差，该指标可体现原始信号与去噪后信号之间的整体偏差信息，数值越小，表明去噪效果越好[16]。其计算公式为：
式f（i）为原始含噪信号，f′（i）为去噪后信号，n为数据长度。
由表1可以得出，随着小波尺度的增大，RMSE先减小，后增大。当小波尺度选取为5时，RMSE最小。
为验证db5为合适的小波基，选取2018年6月18日水管仪NS向数据进行小波变换，图3为应用db5的小波重建信号及其误差。
图3 2018年6月18日水管仪NS向小波重建信号及其误差Fig.3 The wavelet reconstruction signal and its error of NS direction of water-tube tiltmeter on June 18，2018
由图3可以看出，小波重构后的数据与原始数据有一致的曲线形态，且小波重建误差非常微小。
为验证选择db5为小波基的合理性，表2计算了2018年6月18日水管仪NS向数据1-8阶小波分解下的RMSE。从表中可以得出，当阶数是5时，RMSE最小。
表2 不同小波尺度下的均方根误差
综上，应用db5对乌加河台水管仪数据进行小波变换是合适的。
乌加河台DSQ水管仪记录精度较高，仪器稳定，固体潮形态清晰，且能够清晰记录到各种干扰。通过查看水管仪观测日志和数据跟踪分析平台可知，水管仪的干扰因素主要有：自然环境干扰、地球物理事件干扰和人为干扰等。由于水管仪北南向和东西向观测干扰类型及其特征相似，故本文只选取了北南向数据进行了分析。本文统计了2018—2020年水管仪北南向受到的各类干扰及其主要特征（表3）。
表3 水管仪北南向干扰特征（2018—2020年）
在各类干扰中，人为干扰由于受观测人员进洞的影响，其干扰幅度是最大的，干扰形态也是最为明显的，易于识别。人为干扰时段内的数据是错误数据，根据形变学科要求，预处理时应将其做缺数处理，予以直接剔除，不必使用小波分析方法。2018年8月27日至9月1日的降雨事件，在乌加河地区属罕见天气，累计降雨量103.8 mm，造成水管仪观测数据趋势性变化，2019年2月才基本恢复正常。一般而言，前兆观测资料所受干扰的周期远小于固体潮，即有用信号为低频，而干扰为高频，通过小波分析，可将其高频干扰部分逐步分解出来并通过阈值去噪进行剔除，因此降雨干扰用小波分析方法效果不佳。
3.1 地球物理事件干扰
乌加河台水管仪对全球6级以上地震均能清晰记录其同震响应，干扰形态主要为突跳。图3为水管仪记录到2020年6月26日5时5分新疆和田地区（35.73°N，82.33°E）6.4级地震，震源深度10 km，震中距2316 km，最大变化幅度为6.02 ms。
从原始数据中很难看出受地球物理事件干扰的频段信息，应用小波分析方法，对水管仪北南向原始数据进行8层分解，分解出了不同频段的细节信息，d1～d8为该地震的小波分解1～8层细节，周期范围为2～512 min，对应Mallat塔式算法中的D1f（x）—D8f（x）（图2），f（x）为水管仪原始数据序列，采样率为1 min-1，其频带为0～0.5 min-1。
由图4可知，该次地球物理事件的干扰主要体现在细节d1～d3的相对高频部分，对应周期为2～16 min；
在细节d4～d8的相对低频部分，地球物理事件干扰基本被过滤出去，主要为水管仪记录到的包括固体潮在内的低频信息。
图4 2020年6月25日地球物理事件干扰及其小波分析Fig.4 Geophysical event interference and its wavelet analysis on June 25，2020
3.2 气压干扰
气压大幅度突变会引起形变仪器观测曲线的畸变，畸变程度和气压变化幅度成正相关。2019年6月9日15时至21时，乌加河地区气压发生较大幅度变化，短时内变化幅度达3.2 hPa。从图5可以看出：气压变化对水管仪数据造成干扰，干扰形态主要为噪声大，和气压曲线同步性较好；
气压干扰主要体现在细节d1～d3，对应周期为2～16 min；
在细节d4～d8的相对低频部分，气压干扰基本被过滤出去，主要为水管仪记录到的包括固体潮在内的低频信息。
图5 2019年6月9日至10日气压干扰及其小波分析Fig.5 Air Pressure Interference and its wavelet analysis from June 9 to 10，2019
3.3 大风干扰
乌加河台位于内蒙古西部地区，每年春秋两季大风天气较多，是形变数字化观测的主要干扰之一，尤其对水管仪等长基线仪器的干扰明显。2019年10月3日乌加河地区大幅降温，伴有6～7级大风天气，导致水管仪观测曲线出现明显扰动，表现为观测曲线上叠加有毛刺干扰信号，数据噪声变大（图6）。
图6 2019年10月2日至4日大风干扰及其小波分析Fig.6 Wind Interference and its wavelet analysis from October 2 to 4，2019
从图6可知，大风干扰主要体现在细节d1～d5，对 应 周 期 为：2～64 min；
细 节d6～d8为水管仪记录到的包括固体潮在内的低频背景信息。
大风干扰与气压干扰在干扰频次、幅度和形态上均表现出相似性（表1），但利用小波分析可看出二者的不同：气压扰动持续时间较短，扰动频段为1/16～1/2 min-1，扰动主要体现在1/8～1/2 min-1范围内；
大风干扰持续时间更长，分布在更宽的频带（1/64～1/2 min-1）上，能量主要集中在1/16～1/2 min-1，在1/64～1/16 min-1频段上仍存在其相对低频部分，而在这一频段上，气压干扰对水管仪几乎没有产生影响。
小波去噪的方法主要包括：屏蔽去噪法、模极大值检测法以及阈值去噪法等，其中阈值去噪法为最常见的去噪手段[17]。由前面的分析可知，噪声主要为相对高频的细节部分，信号主体在低频的近似部分。阈值去噪先计算出有用信号与噪声的小波系数之间的阈值，再将小于阈值的小波系数置零，实现去噪目的[16]。
阈值降噪应用最广泛的是Donoho提出的软、硬阈值降噪法。硬阈值去噪是将待处理的信号中绝对值小于阈值的小波系数置零，软阈值去噪则是在硬阈值的基础上，将待处理信号绝对值大等于阈值的小波系数再减去阈值，使数据处理边界连续光滑[18]。因此，本文采用软阈值降噪法处理数据。
本次研究针对每一个小波细节选取一个固定阈值，以此避免有用信号的小波系数损失过重，而噪声的小波系数没有得到有效压制[19]。
固定阈值Thr的计算公式[11]为：
式中，σ为标准偏差，n为信号长度。由上述公式，可计算出各种干扰类型的阈值（表4）。
表4 DSQ水管仪各类干扰的阈值
根据表4中的阈值，对水管仪地球物理事件、气压干扰和大风干扰数据进行去噪处理，在去除高频干扰的同时，较好地保留了原始数据中固体等有用信号，明显提高了数据的信噪比，去噪效果良好（图7）。在接下来的工作中，作者拟应用小波阈值去噪方法提取前震异常信号，开展深入研究。
图7 地球物理事件、气压干扰和大风干扰原始信号及小波阈值去噪信号对比Fig.7 Comparison of original signal and wavelet threshold denoising signal of Geophysical Event Interference,Air Pressure Interference and Wind Interference
通过以上分析，可以得到以下结论：（1）乌加河台水管倾斜观测主要受人为干扰、降雨干扰、地球物理事件干扰、气压干扰和大风干扰的影响；
（2）利用小波分解不同频段的信号可知，地球物理事件干扰主要分布在细节d1～d3，对应的周期为2～16 min；
气压干扰主要分布在细节d1～d3，对应的周期为2～16 min；
大风干扰主要分布在细节d1～d5，对应的周期为2～64 min；
（3）应用db5对乌加河台水管仪数据进行小波分析和软阈值去噪是适用的，基本去除了各种类型对原始数据的噪声干扰，去噪效果明显，去噪后的数据曲线光滑度高，同时保留了原始固体潮信息；
（4）运用小波分析和阈值去噪方法可以较好去除干扰信号，能清楚的观测出固体潮等有用信息，对定点形变观测异常识别有一定的作用，可供形变台站参考。
