小波分析方法在钻孔应变数字化资料处理中的应用

倪友忠，吴珊珊，教聪聪，方 韬，叶 青

(上海佘山基准地震台，上海市 201602)

0 引言

利用小波变换方法可对观测资料不同频率段内的信息进行分解,从而提取到观测资料不同频率段的信息，为形变数字化资料的分析处理提供一种新的行之有效的分析方法。文章应用小波变换理论，对佘山台历史记录到的钻孔应变资料，结合宋治平等[1]、张燕等[2]提到的方法和邱泽华[3]文中通过频率响应和小波分解来研究气象变化影响中采用的计算方法，求解得到小波分解细节部分d1-d4的频率段信息以及气象和应变的影响系数，分析及识别观测资料中存在的各种干扰现象。比如，系统故障、气象、地震(同震效应)干扰等。只有通过识别及排除这些干扰因素，才可以进一步对资料进行计算分析，从中识别和提取地震趋势与短临异常。在如何识别和提取地震趋势与短临异常的方法上，李杰等[4-5]在小波分解方法分析形变观测资料的正常背景变化特征和分析跨断层形变异常等文中，都有具体提及。张燕等[6]在大姚6级双震前地震前兆特征一文中，也具体阐述了在地震前震中附近的形变台站，观测资料在小波分解第三层都记录到了相同频段的异常信号。

1 小波分析的基本原理

小波分析是傅里叶方法的发展和延伸，分析信号频率随时间不断变化。其主要作用是能在信号的低频段具有较高的频率分辨率，在高频信号段具有较高的时间分辨率。能在观测信号中，探测到正常信号中夹带的瞬态间的反常信息，并将其反常的成分展示出来。观测台站记录到的形变资料f(x)组成的序列并不是平稳的，而是随时间不断变化的，其中夹杂着许多高频干扰。Fourier变换并不能提取f(x)中的奇异性和突变的信息。小波变换理论的出现，对于解决和展示信号的突变和反常成分，具有很强的实用性。可以利用离散小波变换DWT(Discrete Wavelet Transform)对观测资料的不同频率段的信号进行识别并加以分析。

运用公式：

(1)

计算中，采用a=2k。随着k的增加，信号从最高频向低频分解。当k=0时，信号为采样频率；k=1时将频率二等分，依此类推。对于数字信号f(x)，可以近似地表示为：

(2)

(3)

Mallat的重构算法为：

(4)

2 资料处理结果分析

2.1 小波变换对噪声的识别和排除

数字化前兆资料的特征是信息量大且干扰源多，所以在分析和研究数字化前兆资料时，必须先对干扰源进行有效的识别并加以排除。小波变换理论具备较强的干扰识别和排查能力，所以近年来在数字化前兆领域，出现了不少关于小波变换方面的研究。文章在分析处理钻孔应变前兆信号时，初步采用Morlet小波进行细节分析。

利用小波分析基本原理中上述公式(3)和(4)计算佘山台钻孔应变面应变Sa数据，在小波分解细节部分尺度为1-4阶等高频率段的信息。在图1中分别列出5张小图，分别是Sa面应变数据、d1-d4细节1阶至4阶高频信号数据。d1反映的是其中第1张小图Sa面应变，年变特征清晰，曲线趋势稳定。图中有2处台阶突跳，均为仪器电源故障导致，从图中很难识别出其他干扰信号。利用小波变换方法分解出不同频率域的信息中，发现在2017年10月30日后的Sa数据的d1-d3频率变大，特别是d1阶的频率和原来相比有很大幅度的变化，表明干扰源主要来自高频。d1-d3主要反映的还是高频信号段，这种影响一直持续到2018年12月12日仍未消失。在d4阶时，高频干扰信号有所减弱，固体潮信号段开始显示出来，表明小波分析可以清晰地分解出不同频率域的信息，对数据资料处理有较强的识别能力。

图1 钻孔应变面应变分钟变化及小波d1-d4细节Fig.1 Minute strain variation of borehole strain surface and details of wavelet d1-d4

放大2017年10月30日前后5 d的观测数据(见第24页图2)发现，佘山台四元件钻孔应变仪2017年10月30日电源故障导致断电，观测系统故障，数据造成较大幅度的漂移。仪器维修恢复后，元件3的原始观测曲线出现毛刺增多，曲线变粗。所以，在利用小波变化分析时，d1-d3阶频率段即2～16 min周期段的带通信号存在明显的高频异常现象。之后联系厂家了解情况，得知频繁的断电故障及之前使用UPS不间断电源供电，对钻孔元件3的传感器电容元件造成不可修复的损伤，导致元件3的传感器输出的观测信号曲线毛刺增多，曲线变粗，观测精度有所下降。

2.2 小波变换对气象变化影响的识别和排除

气象变化对钻孔应变观测的影响主要包括气压、气温及强降雨变化等的影响。文章利用小波变换，主要识别和排除气压对钻孔应变观测的影响。用小波方法分解气压和面应变两条观测曲线的观测时间序列，可以对比分析两条曲线d1-d8不同频率上气压和观测应变的变化情况及其对应关系。气象变化在已知研究领域中表明，对水平剪应变几乎没有影响，只对面应变产生影响。所以，只提取钻孔面应变某时段观测的时间序列，利用小波变换方法来加以对比，识别和排除气压对面应变的干扰。

第24页图3的小波分解将Sa面应变信号与气压信号分解至第9层，主要考察离散小波变换各层的细节系数，可以看成原始信号高频部分不同周期段的带通滤波信号。钻孔应变观测的采样间隔是1 min，第1层细节(d1)的信号实际上相当于周期为2～4 min的带通信号，第2层细节(d2)相当于周期为4～8 min的带通信号，依此类推，一直到第九层(d9)。图中，在面应变Sa第1层至第3层细节(d1-d3)中，可以一一对应地找到气压变化对面应变的高频干扰信号。第4层至第9层的细节(d4-d9)中，就显得不明显。第25页图4进行小波分解后，再分别对面应变和气压进行回归分析，得到各频带的气压影响系数。利用得到的一系列随周期变化的气压影响系数(kp)，绘制曲线图4b。

图3a、3b为面应变及d1-d8小波分解第1至第8阶细节部分展示。d1-d4主要是高频信号段展示，d5-d8为低频信号段展示。图3c、3d为气压及d1-d8小波分解第1至第8阶细节部分展示，d1-d4主要是高频信号段展示，d5-d8为低频信号段展示。在d1-d4高频信号区域，可以清晰地找到气压在面应变高频段的干扰信号，d1-d2的细节部分信号最为清晰，d5细节部分固体潮汐信号开始显现，干扰信号逐渐被排除，直至d8几乎看不到气压的干扰信号。由此可知，气压变化对钻孔应变观测的影响主要集中在高频信号段。为了更清晰准确地了解佘山台钻孔应变和气压的对应关系，利用小波分解及回归分析，求得应变气压的影响系数(kp)及曲线。

图2 钻孔应变四元件分钟值曲线Fig.2 Minute value curve of borehole strain four elements

图3 面应变、气压分钟值小波变化d1-d9细节部分展示Fig.3 Details of strain and pressure minute wavelet transform d1-d9

图4a为2016年1月1日至2017年9月30日面应变Sa及应变气压影响系数变化图，可见佘山台面应变的变化趋势大致与气压趋势同步，有明显年变规律；图4b的影响系数显示，佘山地区气压对四元件钻孔面应变Sa的影响。由此可见，应变气压影响系数周期在1 min，数值接近0.3 ns/hPa的时候，开始逐渐上升；周期10 min时，数值达到峰值1.5 ns/hPa；周期大于10 min后，应变气压影响系数反转开始下降；周期达到50 min时，应变气压影响系数趋向于0。该应变气压影响系数曲线同样清晰表明，在佘山台钻孔应变观测中，气压只在短周期段对钻孔面应变有影响，最大影响频率为10 min，对长周期数据无影响。

图4 面应变、气压分钟值及应变气压影响系数图Fig.4 Variation of strain, pressure minute value and strain pressure influence coefficient

图5为气压差分数据(气压变化影响kp×DP曲线)与佘山台面应变差分数据进行对比。从图5清晰找出气压对钻孔应变的影响，表明气压对佘山台面应变存在影响，影响周期主要集中在0～10 min的频率域范围内，小波分解细节d1-d3周期段。

图5 面应变分钟值差分、气压变化影响(kp×DP)曲线Fig.5 The strain minute difference and the effect of pressure change (kp×DP) curve

2.3 小波变换对地震(同震效应)的识别和排除

钻孔应变观测可以清晰地记录到远震、大震或近震、中小地震的同震效应(即地震面波)。但这种同震效应幅度大，夹杂在观测曲线中。对于前兆分析来说，也属于一种干扰，这种干扰在趋势异常分析时，应予以识别和排除。小波分解方法是排除此种干扰最直接有效的方法之一(见图6)。

图6 小波变换识别面应变同震效应d1-d4细节部分Fig.6 Wavelet transform identifies details of surface strain coseismic effect d1-d4

2018年9月28日18:02:44印度尼西亚发生7.4级地震，上海市佘山台钻孔应变四个元件测项，都清晰地记录到了同震效应。图6展示了钻孔应变仪元件1测项记录到该地震的同震效应,利用小波变换方法提取了第1阶至第4阶细节的高频信号。说明小波方法对于同震效应的识别、提取效果较好。

2.4 小波变换对钻孔应变前兆资料趋势与短临异常的识别

选取2013至2018年钻孔应变面应变Sa的日均值进行小波分析。第26页图7a、7b、7c分别是面应变Sa未剔除干扰的日均值变化图、Sa剔除干扰后的日均值变化图、Sa小波分解d3阶细节部分高频信号。Sa未剔除干扰日均值变化图中的干扰包含地震同震效应、仪器故障缺数、气象干扰等。在日均值图中表现为脉冲突跳，变化幅度比较大。利用matlab差分函数进行差分计算后，剔除部分干扰，展现出图7b面应变Sa的日均值长期变化趋势。从Sa的日均值长趋势变化图中，了解到佘山地区的面应变长期受压向下。图7c是Sa的日均值进行小波分解后第3阶细节部分的高频信息图。从图中可以看出，在2014年3月3日琉球M6.7地震前，Sa的日均值有明显鼓包现象，小波分解第三阶细节部分的高频信号幅度超出平时的2倍以上，连续且密集，持续了近半年左右，震后该高频脉冲信号幅度明显降低；2015年11月14日东海海域M7.2地震发生前，Sa的日均值鼓包不太明显，但Sa的日均值小波分解d3细节部分出现了持续4个月左右的密集，幅度超过平静日2倍以上的高频脉冲信号。2016年1至6月，Sa的日均值小波分解d3细节部分也曾陆续出现幅度较大的高频脉冲信号，但由于没有出现Sa的日均值鼓包，且高频信号不连续，所以该处的d3细节部分反映出来的高频信号不是地震前兆异常。

图7 面应变日均值变化图、小波分解第3阶细节部分Fig.7 Daily mean strain variation diagram, third-order detail of wavelet decomposition

3 结论与讨论

综上所述，得出如下结论。

(1) 中小波变换在数字化前兆资料分析中可以展现d1-d4阶细节部分的高频信号，能清晰地识别出地震(同震效应)、气象干扰(气压等)、仪器故障等干扰，具有较强的干扰识别与排除功能。

(2) 小波变换方法对钻孔应变长趋势的识别、提取和短临趋势异常也有较强的识别能力。图7在琉球群岛M6.7和东海海域M7.2地震前，在d3阶细节部分都明显地提取到了幅度较大，周期大概持续5～6个月的高频信号；且Sa的日均值有鼓包现象，地震发生在鼓包回落恢复的后期，震后该现象明显消失，故判断很可能为地震前兆异常。

(3) 2018年7月后，Sa的日均值第3阶细节部分开始出现明显的高频脉冲，Sa的日均值也有明显鼓包，这一异常现象是否为海域地震前兆异常，值得进一步关注研究。