0 引言

微地震事件通常具有持续时间短、能量微弱、信噪比(signal-to-noise ratio, SNR)低的特点(赵博雄等，2014)，有时微地震的有效信号甚至能被外部噪声所淹没，而微地震事件信号的信噪比极大地影响了震源定位的精度.因此，必须首先对监测到的微地震数据进行滤波去噪，提高微地震数据的信噪比，从而提高震源定位的精度.目前，已经有很多滤波方法应用于微地震数据去噪中，如卡尔曼滤波(Kalman Filter，KF)(Baziw and Weir-Jones, 2002；Baziw et al., 2004；邓小英等，2014)、自适应滤波(宋维琪等，2013)、小波分析(王正蕾，2009；Du and Zhang, 2010；宋广东等，2011)，其中有些方法是基于模型的滤波技术，如KF滤波已经成功应用在基于机理模型的微地震信号处理中.

2002年，Baziw和Weir-Jones (2002)提出了一种基于微地震信号机理模型的KF滤波算法.通过对微地震信号进行机理分析，利用衰减震荡周期过程模拟微地震子波，将微地震记录的背景噪声当作高斯马尔科夫过程进行处理和建模，在一系列的假设和线性化之下得到了微地震信号的机理模型.两年后，Baziw等(2004)改进了当前检测微地震P波初至的卡尔曼滤波器的设计，进一步提高了自动识别P波初至的能力.

由于机理建模是以分析大量的信号数据为基础的，使得机理模型在实际应用中很难获得.为了克服这个难题，本文采用衰减震荡周期信号模拟微地震子波，通过辨识方法，如Steiglitz-McBride迭代方法，建立微地震信号的ARMA模型.为了在微地震信号处理过程中应用KF滤波，将ARMA模型进一步地转换为适用于卡尔曼滤波算法的状态空间模型(邓自立，2005).在此基础上，设计了适合于微地震信号去噪的KF滤波算法.在微地震信号辨识建模过程中，消除了机理模型在线性化处理过程中所存在的一些假设，提高了模型的精度，从而改善了KF滤波效果.

通过对微地震合成信号和地面监测数据进行滤波处理，基于辨识模型的KF滤波算法能够有效抑制微地震监测数据中的随机噪声，显著提高信噪比，从而验证了该辨识模型和滤波算法的准确性和可行性.

1 微地震信号辨识建模 1.1 ARMA模型

自回归滑动平均模型(Auto-Regressive and Moving Average Model，ARMA)(邓自立，2005)是近代时序分析中描述离散平稳随机信号的最常用模型之一，可看做白噪声激励一个线性系统产生的时间序列模型.

离散系统的输入输出模型可用差分方程的形式表示为

(1)

将式写成更一般的形式，即

(2)

式中，z^-1为单位滞后算子，A(z^-1)=1+a₁z^-1+a₂z^-2+…+a_pz^-p, B(z^-1)=b₁z^-1+b₂z^-2+…+b_qz^-q, p和q分别是自回归部分和滑动平均部分的阶数，a_i(i=1, 2, …, p)和b_j(j=1, 2, …, q)分别是自回归和滑动平均系数.

特殊地，当ARMA模型中输入项的各个系数b_j都为0时，模型变为

(3)

式(3)称为p阶自回归模型，记为AR (p).

当ARMA模型中输出项各个系数a_i都为0时，模型变为

(4)

式(4)称为q阶滑动平均模型，记为AM (q).

参数辨识建模需要确定辨识模型的类型，模型的阶次p和q，进而估计出模型系数a_i(i=1, 2, 3…, p)和b_j(j=1, 2, 3…, q).模型确定后，需要对模型的参数进行估计.与信号的分析类似，模型参数的估计方法可分为时域和频域两大类，本论文采用时域方法进行参数估计.具体方法是在确定系统为AR模型或者是ARMA模型之后，从给定的系统响应信息，或者是输入输出序列的相关性中估计出模型参数的具体值.

1.2 微地震子波辨识建模

本文采用指数衰减周期震荡信号合成微地震子波(王维波等，2012)，如式(5)所示，与实际监测到的微地震信号波形比较相似.公式为

(5)

式中，A₀为信号的初始振幅，h为阻尼因子，w为角频率，t₀为信号的初始相位.通过对阻尼因子和角频率等参数选取合适的数值，可得如图 1所示的微地震子波.

通过对微地震合成子波的样本数据进行相关性分析(王少水等，2009)，结果表明，样本数据的自相关函数_k和偏自相关函数_kk都呈现拖尾特性，如图 2所示.

拖尾是指，当样本数据点数k增大时，样本数据的自相关函数_k以负指数的速度衰减，其曲线好像是拖着一条尾巴.截尾是指，当样本数据的点数k等于p时，样本数据的偏自相关函数_kk不为0，而当k>p时，样本偏自相关函数_kk等于0，好像是曲线被剪掉了尾巴一样.

由相关分析法的分析结果可知，样本数据的自相关函数_k和偏自相关函数_kk都呈现拖尾特性，又由表 1给出的相关分析法结果与辨识模型类别的对应关系，可以确定微地震合成子波的辨识模型的类别为ARMA模型.此外，相关分析法可以判定辨识模型的类别，但无法判定辨识模型的具体阶次.

为了确定辨识模型的阶次，本文通过奇异值分解法(Singular Value Decomposition，SVD)判定ARMA模型的阶次(王少水等，2010；王树青等，2012).奇异值分解法是应用比较广泛的完全正交分解法，通过构造特殊矩阵并进行SVD分解，将归一化处理后的奇异值按从大到小的顺序进行排序，画出归一化奇异值曲线(图 3).通过观察和分析奇异值曲线，找出曲线趋向于水平或者突然下降的点，该点所对应的奇异值的个数即为所求的模型阶次.该方法简单实用，判定结果也较为准确.

通过观察图 3所示的奇异值曲线，在第2个奇异值处曲线突然下降并趋向于水平，表明SVD法判定的模型阶次为2阶.在随后的建模仿真实验中，当辨识模型的阶次为2阶时，辨识模型的响应曲线与实验数据曲线高度吻合，所以最终选择辨识模型的阶次为2阶，即为ARMA (2, 1)，其输入输出表达式如式(6)所示，辨识结果如图 4所示.公式为

(6)

其中

采用该模拟信号合成微地震信号时，在有效信号A(t)上叠加均值为0，标准差为σ的高斯白噪声N(t)~N(0, σ²).用一个表达式描述微地震监测站点接收的微地震数据，即

(7)

从信号能量角度定义微地震信噪比R_SN为

(8)

即有效信号的能量与噪声能量之比.微地震记录的信噪比越高，微地震记录就越接近微地震信号，就越能有效地识别出P波、S波，从而提高震源定位的精度.

1.3 微地震子波状态空间模型

为了验证辨识模型的准确性和方法的可行性，将辨识得来的ARMA模型进一步转换为状态空间模型，公式为

以此状态空间模型合成微地震子波，并叠加一定的高斯白噪声(图 5).

2 Kalman滤波 2.1 Kalman滤波算法

基本KF滤波(Kalman，1960)只适用于线性离散系统，要求精确已知系统的数学模型和噪声的统计特性.假设系统中所有的噪声信号都服从高斯分布，并且相互独立.

设线性离散系统的状态空间模型(王静波，2012)为

(9)

(10)

其中，式(9)为系统的状态方程，式(10)为系统的观测方程.x_k和z_k分别为系统在k时刻的状态变量和观测变量，u_k为系统在k时刻受到的控制量.A为状态转移矩阵，B为控制矩阵，H为观测矩阵.假设系统噪声w_k和观测噪声v_k是互不相关的零均值高斯白噪声，并且它们是相互独立的.公式为

式中，Q和R分别为系统噪声w_k和观测噪声v_k的协方差.

KF滤波在确保最小均方误差的前提下，以迭代递推的形式对信号状态进行最优估计.滤波算法由时间更新和测量更新两个过程组成，时间更新是对当前时刻的系统状态和先验误差协方差进行预测，测量更新根据系统当前时刻的观测值修正时间更新过程中的系统预测值.滤波算法(王静波，2012)简述如下：

1) 系统的先验状态估计

(11)

为系统在k时刻的先验状态估计，为系统在k-1时刻的状态估计，u_k为系统在k时刻的控制作用.

2) 先验误差协方差估计

(12)

P_k^-为系统在k时刻的先验误差协方差估计，P_k-1为系统在k-1时刻的误差协方差，上标T表示矩阵的转置.

3) 计算卡尔曼滤波增益

(13)

K_k为卡尔曼滤波算法在k时刻的滤波增益，上标^-1表示矩阵的逆.

4) 更新系统状态

(14)

5) 更新误差协方差

(15)

单位矩阵I的维数与辨识模型的系数矩阵A保持一致.

系统初始条件为

(16)

(17)

Kalman滤波算法的实现步骤如下：

(1) 初始化系统状态及协方差矩阵P₀；

(2) 计算系统的先验状态估计和先验误差协方差P_k^-，完成时间更新过程；

(3) 计算卡尔曼滤波增益K_k、系统状态、误差协方差P_k，完成状态更新过程；

(4) 判断滤波算法是否结束，如果没有结束则返回步骤2，直至算法结束为止.

2.2 合成微地震信号滤波处理

为验证辨识模型的准确性和KF滤波算法的应用效果，对叠加不同程度高斯白噪声后信噪比为10%、30%、50%、100%的微地震合成信号进行KF滤波.结果表明，微地震信号中的高斯白噪声均得到明显抑制，效果显著.由于篇幅限制，本文只列出信噪比100%的微地震合成信号的滤波结果(图 6).

3 实际数据处理

为了进一步验证辨识建模方法的适用性和辨识模型的准确性，将该辨识建模方法和滤波算法应用于实际采集到的微地震监测信号中，笔者从采集到的高信噪比的微地震数据中提取出实际微地震信号，并对提取得到的实际微地震信号采用上述建模方法进行辨识建模，可得实际微地震信号的辨识模型为ARMA (15, 14)，其输入输出表达式如式(18)所示，实际微地震信号的辨识结果如图 7所示.

(18)

其中，

通过对实际采集到的微地震监测信号进行辨识建模，将得到的ARMA模型进一步转化为状态空间模型.在此基础上，设计了适于实际微地震监测数据去噪的卡尔曼滤波算法.为验证基于实际微地震信号辨识模型的卡尔曼滤波算法在实际微地震监测数据中的应用效果，笔者对某压裂气井的地面微地震监测信号(图 8)进行处理，滤波效果如图 9所示.

通过对实际采集到的微地震地面监测数据进行卡尔曼滤波处理，KF滤波对实际微地震数据中的随机噪声滤波效果明显，能够显著提高微地震记录的信噪比，从而验证了该辨识模型的准确性和滤波算法的可行性.

4 结论 4.1

本文采用衰减震荡周期信号合成微地震子波，针对微地震典型信号建立了ARMA辨识模型，辨识建模过程中消除了机理模型中因假设和线性化处理而存在的模型误差.为了在实际采集到微地震监测数据中应用KF滤波，笔者将辨识得来的ARMA模型进一步转化为适于卡尔曼滤波的状态空间模型，并研究和设计了卡尔曼滤波的实现算法.根据理论模型和对微地震子波信号的辨识结果可知，本文通过参数辨识方法建立的微地震典型信号的辨识模型误差很小，精度较高.

4.2

通过对叠加不同程度高斯白噪声、信噪比分别为10%、30%、50%、100%的微地震合成信号进行KF滤波，噪声得到显著抑制，滤波效果较好.实际采集到的微地震监测数据滤波结果表明，KF滤波可以较好地抑制微地震监测数据中的随机噪声，显著提高微地震监测数据的信噪比，从而验证了该辨识模型的准确性和基于辨识模型的KF滤波算法的可行性.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！