傅里叶变换在深州井资料处理中的应用

尹宏伟，梁丽环，韩文英，郭学增，张 蕾，信世民

(1. 河北省地震局深州地震台，河北 深州 053800； 2. 河北省地震局兴济地震台，河北 沧州 061721)

傅里叶变换在深州井资料处理中的应用

尹宏伟1，梁丽环1，韩文英1，郭学增2，张 蕾2，信世民2

(1. 河北省地震局深州地震台，河北 深州 053800； 2. 河北省地震局兴济地震台，河北 沧州 061721)

深州井水位日常观测中，经常会因为仪器供电不稳、雷击等原因出现一些幅度很大的脉冲、突跳及数据缺测，严重影响了观测资料的内在质量。本文采用数字信号处理技术[1]对深州井水位资料进行分析处理，首先使用分段线性插值法对缺测数据进行插值处理，然后用快速傅里叶变换对水位数据进行频谱分析，绘制振幅谱，分析水位数据中的不同周期成分，选取恰当的截止频率进行滤波处理，去除了干扰，获得了井水位的变化动态，提高了资料质量。最后应用滤波后的水位数据做一些震例分析。

新泽5井；资料干扰；傅里叶变换；滤波；震例分析

0 引言

快速傅里叶变换的频谱分析方法能够从各种复杂的数据序列中分离出不同频率的信息成分，计算各种频率成分振幅的大小，并通过将频率域中某些频率成分的振幅置零，运用傅里叶逆变换到时间域而达到滤波的目的。

傅里叶变换因其强大的信号处理功能，广泛应用于通讯、电子、石油物探、生物医学等社会的各个研究领域，近年来国内越来越多的同行将其应用到地震前兆数据的分析中。肖本夫等[2]运用快速傅里叶变换对2008年5月12日汶川地震的波形数据进行频谱分析，并设计FIR低通数字滤波器进行滤波，去除干扰，提取了有效数字地震信号。李瑞华等[3]用快速傅里叶变换的频谱分析和数字滤波技术对山东省聊古1井2年的气氡逐时值与日均值进行分析处理，寻找其中隐含的周期结构，从中获得一些新的认识。李敬等[4]用快速傅里叶变换对汕头台数字地震仪记录的数据信号进行频谱分析，并设计IIR低通滤波器滤除汽车的干扰，获取有用信息。 而使用各种数学分析的方法对观测资料进行排除干扰的处理也被许多地震工作者所广泛采用[5-8]。笔者借鉴前人的经验，运用快速傅里叶变换对深州井水位数据进行频谱分析，绘制振幅谱，找出干扰信息的频段，进行滤波，去除脉冲、突跳，数据缺测等干扰，取得了较好的效果。

1 傅里叶变换原理

傅里叶变换(Fourier变换)是一种线性的积分变换，因其基本思想是由法国数学家和物理学家约瑟夫·傅里叶最早提出的，故以其名命名。“任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加”。根据这一原理创立的傅里叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。在不同的研究领域，傅里叶变换具有不同的变体形式，如连续傅里叶变换和离散傅里叶变换。而快速傅里叶变换(FFT)是离散傅里叶变换的快速算法，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进而获得的。

一个模拟信号，经过采样之后，变成数字信号，对数字信号就可以做FFT变换了。为了方便运算，通常采样点数N取2的整数次方。假设采样频率为Fs，信号频率为F，采样点数为N。那么FFT之后的结果就是一个为N点的复数。每一个点对应着一个频率点，某一点n(n从1开始)表示的频率为：Fn=(n-1)×Fs/N；该点的模值除以N/2就是对应该频率下的信号幅度；该点的相位即是对应该频率下的信号的相位。第1个点表示直流分量(即0 Hz)，而最后一个点N的再下一个点(实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第1个点分做两半分，另一半移到最后)则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。以频率作为横坐标，以信号的幅度作为纵坐标绘出的图形就称为Fourier振幅谱。

由于FFT结果的对称性，通常只使用前半部分的结果，即小于采样频率一半的结果。但是当使用FFT进行滤波时，则必需考虑所有的频率成分。

傅里叶变换是数字信号处理领域一种很重要的算法。有些信号在时域上很难看出什么特征，但是如果变换到频域之后，就很容易看出特征了。傅里叶变换将原来难以处理的时域信号转换成了易于分析的频域信号(信号的频谱)，并可以利用一些工具对这些频域信号进行处理、加工，然后利用傅里叶逆变换将这些频域信号再转换成时域信号，从而达到去除干扰信息、提取有用信息的目的。

2 深州新泽5井概况

新泽5井是河北省地震局地下水位观测井网的一口深井，位于河北省中部深州市境内，地处冀中凹陷深泽低凸起南缘(或深泽—刘村构造带南缘)，滹沱河冲积扇下游，地势微向东倾，第四系松散沉积物厚450 m，在地质上属于应力易于集中的地区。观测井为石油探井，井深3 364 m，主要含水层为奥陶系灰岩、白云岩，次为石炭系底部白云质灰岩。含水层为承压含水层，地下水类型为岩溶裂隙承压水，底部水温83.9 ℃，水质矿化度高 (含盐量大于10 g/L)。该井地下水补给区在太行山麓及周围边缘，排泄区在冀中凹陷文安、河间一带。由于该井深度大，隔水顶板埋深深(2 835.7 m)，套管水泥固体止水措施效果佳，所以含水层封闭性较好，不受降雨直接渗入补给的影响[9]。新泽5井目前采用DSC-1A型数字观测仪进行观测，观测内容有静水位、中层水温、气压。

3 新泽5井水位正常变化动态

深州井水位呈逐步下降的多年趋势变化；受固体潮和大气压力的影响明显，总体上呈夏高冬低的周期性年变。2011年水位观测基本未受到干扰因素的影响，水位变化曲线清晰、连续、完整(图1)。

图1 深州井2011年水位曲线图

4 干扰资料的分析处理

受雷击、仪器故障及供电电压不稳的影响，深州井2005年水位资料大量断记，而且出现了很多大幅度的脉冲、突跳(图2)，大幅度的脉冲突跳严重压制了水位正常的变化动态，震兆信息完全淹没在噪声干扰之中。

图2 深州井2005年水位曲线图

本文采用以下方法对水位数据中的各种干扰逐步进行处理。

4.1缺测数据的插值及插值后数据的频谱分析

线性插值是一种简单的插值方法，假设已知坐标(x0,y0)与(x1,y1), y0到y1随x线性变化。要得到[x0，x1]区间内某一位置x所对应的y值。公式(y-y0)/(x-x0)=(y1-y0)/(x1-x0)是成立的，因为等号两边同为直线的斜率。由上式可求得y=(y1-y0)*(x-x0)/(x1-x0)+y0，x0、y0、x1、y1是已知的，x是需要插值的点。x(时间序列)取不同的数值，便可求得其相应的y值(对应时间序列中不同时刻的观测数据值)。分段线性插值是把整个需要插值的区间分成若干个小区间，然后在每个小区间内分别进行线性插值。MATLAB中的interpl 函数可实现分段线性插值的计算。

地震前兆数据库中缺测数据用“999999”代替，找到深州井2005年水位数据中“999999”的信号，使用分段线性插值法对水位数据进行插值处理。然后利用快速傅里叶变换对插值后的水位数据进行频谱分析，在0～0.04 min-1频段绘制深州井2005年原始水位数据振幅谱(图3)。

图3 深州井2005年水位数据幅频曲线图

分析水位数据的幅频图，其特征是从低频向高频振幅迅速衰减，其优势频率集中在非常窄的低频范围内，是水位的正常变化动态。而振幅较小的高频部分为噪声干扰，需要进行滤波处理。但是图2中的几次大幅度负脉冲显然不是水位的客观变化，而是观测仪器问题引起的脉冲，如果携带这些大幅度脉冲直接进行滤波，得出的结果会在这些脉冲相对应的时序图上出现异常变化，资料分析时就会把这些异常变化误认为是地震异常，所以在滤波前要先把这些大幅度的负脉冲剔除掉。

4.2 错误数据的剔除

因为深州井水位近几年变化趋势为逐年下降的变化动态，根据2006—2007年水位数据的情况，2006年水位埋深最低为1.593 m，2007年最低为2.220 m，考虑可能出现的水位异常变化，所以把2005年水位数据中埋深大于2.6 m的脉冲突跳视为错误数据应该是合理的，对错误数据进行剔除，然后使用分段线性插值法对剔除的数据进行插值处理，结果如图4a所示，水位数据中大幅度的负脉冲被剔除后，可以基本反映出水位的变化趋势，但是仍然存在大量的小幅突跳、脉冲等高频干扰，掩盖了水位变化中的一些短周期信息，需要进一步的滤波处理。

4.3 干扰数据的滤波及滤波后数据的频谱分析

采用快速傅里叶变换对水位数据做低通滤波处理，为了达到最佳滤波效果，做到既能滤除脉冲、突跳等高频干扰信息，又能保留有用的信息，笔者经过反复多次试算，最后选取0.000 22 min-1作为截止频率，大于0.000 22 min-1频段的振幅置零，小于等于0.000 22 min-1频段的振幅保持不变，对数据进行低通滤波，并对滤波后的水位数据进行频谱分析，在0～0.01 min-1频段绘制滤波后水位数据的振幅谱，结果如图4b、图5所示。图4b是经过插值、剔除错误数据、再插值、快速傅里叶变换滤波等一系列处理后的水位数据图，与图2比较，可以看到，数据中的脉冲、突跳及缺测等各种干扰被完全排除了。

t/(×105 min) a 剔除突跳点后的深州井2005年水位曲线图

t/(×105 min)b FFT滤波后的深州井2005年水位曲线图

图5 经过FFT滤波后的深州井2005年水位数据幅频图

图5是经过滤波后的水位数据幅频图，与图3进行比较，可以看出经过快速傅里叶变换滤波后，水位数据中大于0.000 22 min-1频率的高频干扰被完全去除了，滤波效果非常好。相对于其它滤波器，快速傅里叶变换滤波虽然运算速度稍微慢些，但它是最“精确”、最“彻底”、最“干净”的滤波。

5 震例分析

通过对深州井水位数据进行快速傅里叶变换滤波处理，结果表明：在一些地震发生前后，深州井水位出现了不同程度的异常变化(图4、图6、图7)。分析这些震例发现：①记录到的大部分为震源比较浅的远大震；②每个地震发生前，临震阶段的水位异常变化形态特征都比较相似，大致分为2类，一类是“下降—转折上升—发震”，如2005年7月印尼苏门答腊北7.3级、2005年8月河北蔚县3.6级、2005年11月江西九江、瑞昌间5.7级、2006年7月河北文安5.1级地震；一类是“下降—发震—转折上升”，如2005年12月坦桑尼亚7.0级、2008年5月四川汶川8.0级、2008年10月新疆克孜勒苏柯尔克6.8级地震。③河北蔚县地震及汶川地震后水位均出现了大幅度的上升阶变，且短时间内未恢复到震前状态。④2006年7月河北文安地震发生前70 d，水位出现了大幅度的振荡型异常。

t/(×105 min) 图6 深州井2006年水位数据经分析处理后的曲线图

t/(×105 min)

地震是地壳介质在构造力作用下发生破裂的一种表现形式，而地震发生前能量积累的时间和空间尺度是很大的，其孕育过程中能量积累导致的应力场变化可能传递至很远的距离[10]。近年来华北地区的地下水动态在强震前短临阶段的异常大多表现为 “下降—转折上升—发震”这样一种模式[11]。而以上讨论的震例大部分遵循了这样的模式，其具体解释可能是地震孕育过程中能量积累导致应力场变化传递至深州井，作用在其含水层，引起含水层的岩石骨架发生拉张和压缩的变化，导致水位的下降和上升。以上讨论的情况仅仅是可能性，至于这种变化是否为地震前兆异常仅凭深州井单一手段的资料还难以下结论，尚需其它手段的观测资料进行佐证。

深州井对坦桑尼亚地震和新疆克孜勒苏柯尔克地震反映出来的水位下降而后转折上升的变化可能是同震效应；河北蔚县地震及汶川地震后出现的水位上升阶变，应属于震后效应。关于地下水位对地震的响应机理许多专家都做出了评述，兰双双[12]认为地下水位对远震的响应主要是由于含水层介质受到地震波应力的作用，其响应形态主要以振荡型和阶变型为主，异常出现的时间较晚。对近震的响应主要是含水层介质受到区域构造应力和地震波应力共同作用的结果，震中距越小，含水层受到震源构造应力场的控制作用越大，其响应形态主要以阶变型、脉冲型和振荡型为主。杨竹转[13-14]认为地震波的作用使含水层介质受到挤压，空隙堵塞，引起含水层渗透系数的减小是导致水位上升的可能原因。地震波对井孔含水层系统的变化仅起到触发作用，其变化方式是由观测井局部的地质构造和水文地质条件决定的，同一口井水位总是以固定的方式对地震波作用做出响应，即总是上升或总是下降。深州井水位的映震特征恰好印证了上述观点。

深州井水位对河北蔚县地震和四川汶川地震的响应特征与陕西渭南井水位对汶川地震的响应特征[12]非常相似，首先出现脉冲向下的异常，而后水位回升，表现出明显的阶升现象。其具体解释可能是震时含水层介质发生了弹性形变，从而导致地下水位表现为脉冲型异常，因为水位数字化观测是分钟值采样，周期小于1 min的水位震荡是记录不到的，所以短时的水震波可能会缺失，反映在数字化水位图上只是脉冲向下的变化。而震后的水位阶升现象很可能是含水层介质在地震波应力或区域构造应力或是二者的共同作用下，受到挤压，空隙堵塞，导致含水层渗透系数减小的原因造成的。

深州井水位对不同地震有不同的响应特征，比如有些地震发生后水位出现大幅度的上升阶变，而有些地震发生后水位没有反映或反映不明显。映震特征的不同不仅与井孔自身条件有关，还可能与震源机制、震源深度、震级以及地震波传播途经的地质构造有关。一个井孔的映震能力还有可能因为地震的触发作用而发生改变。

6 讨论与结论

运用快速傅里叶变换的方法对深州井水位观测资料进行分析处理，取得了较好的效果。经过插值处理后，补齐了缺测的数据；经过快速傅里叶变换滤波后的水位数据，无论在时间域，还是在频率域，均清晰地显示出其脉冲、突跳等高频干扰被完全排除了，水位曲线变得光滑了，更直观地看清了水位的中长期变化趋势，一些短周期的水位变化也清晰地凸显出来，提高了水位观测资料的内在质量，为地震预测研究提供了更为可靠的分析依据。

在运用快速傅里叶变换进行滤波时，截止频率的选取非常关键，要经过反复多次试算, 以得到最佳截止频率，使其滤波后既能有效地去除干扰，又能最大程度地保留有用的信息。滤波后的数据曲线在首尾两头儿出现了一些畸变，这是由数据截断效应造成的。

本文所讨论的水位数据中的干扰主要是高频噪声，所以只采用快速傅里叶变换的低通滤波方式去除高频干扰从而识别出低频信息，但是水位数据中高频部分也可能含有震兆信息，比如临震前的水位突跳，而此方法在去除干扰的同时，也会把处于相同频段的震兆信息一并滤除，这对于使用水位资料进行地震分析研究是不利的。如何从同一频段的噪声干扰中提取震兆信息？笔者将会在今后的工作中继续深入研究。

[1]万永革.数字信号处理的MATLAB实现[M].北京：科学出版社，2012.

[2]肖本夫，万永革，祁玉萍.强干扰环境下有效数字地震信号的提取[J].华北地震科学，2011,29(1)：6-9.

[3]李瑞华，张海燕.频谱分析方法在气氡周期分析中的应用[J].华北地震科学，2008，26(2)：40-44.

[4]李敬，甘延锋，黄友明，等.数字地震记录中干扰波的排除[J].防灾技术高等专科学校学报，2004,6(3)：20-25.

[5]邓亮，刘春平，万飞，等. 从水位观测数据中排除非构造因素影响的状态空间分析方法[J].华北地震科学，2009，27(2)：17-21.

[6]戴勇，高立新，杨彦明，等. 基于小波变换方法的包头台形变分析[J].华南地震，2013，33(4): 39-46.

[7]关玉梅，王紫燕，王青平，等. 数学形态滤波法在GPS基线分析中的应用[J]. 华北地震科学，2013，31(3):41-46.

[8]鲁权，邢西淳，李西京，等. 泾阳台数字地磁信号的干扰分析及去噪处理[J].华南地震，2013，33(1):49-54.

[9]汪成民，李宣瑚，王铁城，等．中国地震地下水位动态观测网[M]．北京：地震出版社，l990：300-301．

[10]白玉柱，徐锡伟，徐杰. 断裂分段之间不同相互作用对断裂运动的影响[J]. 华北地震科学，2012，30(2):33-38.

[11]张素欣，盛艳蕊，单连君，等.地下水动态多尺度地震预测技术初探[J].华北地震科学，2011，29(4)：1-6.

[12]兰双双，迟宝明，姜纪沂，等.地下水位对近震和远震异常响应的比较—以汶川地震和苏门答腊地震为例[J].吉林大学学报，2011，41(1)：145-152.

[13]杨竹转. 地震引起的地下水位变化及其机理初步研究 [D].北京：中国地震局地质研究所，2004.

[14]杨竹转，邓志辉，赵云旭，等.云南思茅大寨井水位同震阶变的初步研究[J].地震学报，2005，27(5)：569-574.

Application of Fourier Transform on Data Processing of Shenzhou Well

YIN Hong-wei1, LIANG Li-huan1, HAN Wen-ying1, GUO Xue-zeng2,ZHANG Lei2, XIN Shi-min2

(1.Shenzhou Seismic Station, Hebei Shenzhou 053800, China;2. Xingji Seismic Station, Hebei Cangzhou 061721, China)

Daily observations of Shenzhou well often have some sharp pulse, jumps and data missing because of unstable power supply, lightning and other reasons, which is a serious impact on the intrinsic quality of observational data. This paper uses digital signal processing techniques to analyze and process the data of Shenzhou well water level. First, we use piecewise linear interpolation method to interpolate the missing data. Then we use fast Fourier transform spectral analysis on the data of well level, draw its amplitude spectrum, and analyze its different period ingredients. Then by selecting appropriate cut-off frequency for filtering, we successfully get rid of the interference and get well level changes dynamically and improve the data quality. At last, we use the filtered water level data on several earthquake case analyses.

Xinze No.5 well; data interference; Fourier transform; filtering; earthquake case analysis

10.3969/j.issn.1003-1375.2014.04.008

2014-08-18

河北省地震局地震科研基金项目“河北省井水位破年变异常分析”

尹宏伟(1971—)，男，工程师，主要从事地下流体观测工作. E-mail：yinhongwei1971@163.com.

P315.723

A

1003-1375(2014)04-0039-05