0 引 言

微地震正演模拟，目前一般是沿用地震勘探的有限差分方法(王守东，2003)和射线追踪(Leidenfrost et al.，1999；宋维琪等，2003)方法.随着微地震技术的应用的不断扩大，实际生产应用中发现了许多问题(Raab and Coyne，1997).其中最关键的问题是定位精度(宋维琪和杨晓东，2011)和实时监测的问题(Wolhart et al.，2006；De Meersman et al.，2009).

另外，微地震定位结果的最优解释，也是面临急需解决的问题(Pei et al，2009)，而这其中涉及微地震震源的反演，为了反演震源，走时反演已远不能满足问题的要求(张晓林等，2013)，快速的波场正演，是摆在我们面前的重要任务.考虑井中接收观测方式以及激发和接收边界条件的复杂性，有限差分方法正演结果的精度不能够保证.

鉴于以上分析，研究了微地震波场模拟的快速数值方法，把解微分方程的模式匹配理论(Liu et al.，1990；Fan et al.，2000；Li et al.，2005)引用到本研究中，将以上模型的正演问题，变成纵向采用解析方法横向采用有限元的组合方法.

1 基本方程和边界条件

对于井中观测的微地震监测问题，由于观测是沿垂直方向观测的，可以看作是二维问题，满足波动方程(王晨龙等，2013)

以及在XOZ平面上给出的边界条件.对于二维问题，不存在点震源.使用线震源来代替三维空间的点震源.线震源可以由点震源合成.波动方程(1)在柱坐标系中的表示为

当波场具有轴对称性质时，波函数与方位角无关，波动方程变为

将上式在时间方向进行傅里叶变换得：

其中F为某个位移分量的傅里叶变换.为了寻找方程的通解，可使用分离变量方法，令F=R(r)Z(Z)，则

式中，λ_m为本征值;R为相应的本征函数.

本征函数应满足如下边界条件：

(1)自然边界条件

具体为

2 正演模拟混合方法

如果设定地层的速度模型为水平层状模型，并且波场在径向变化以源为中心具有轴对称特点(Chew et al，1984).问题求解在垂直方向采用解析方法、径向采用数值方法的混合方法，公式为

对以上方程采用有限元法求得数值解.令： g _m=(g_m1，g_m2，…，g_mm)为径向上的一组完备的基函数，那么特征函数f_m(r)可以近似地表示为

经化简，方程(6)等价于如下矩阵形式的本征值方程:

其中：

式中，C _m为由N个本征向量 c _i组成的特征全矩阵，可以证明： C _m^T P _m C _m= I .

式(10)是一个广义特征值问题，通过求解广义特征值，可以得到各层的本征值和本征向量.考虑位移应力连续条件，方程(5)的解为

因此，由式(7)和(10)，第M层中 方程(2)的解可以写为

第M+1层中方程(4)的解可以写为

在Z=Zm的层界面上，位移和应力分量满足连续性条件，即，

对上式两边分别左乘 C _m^T g _m，并且从0到无穷对R积分，根据正交关系 C ^T BC = I，可以得到

其中， P _m+1，m= C_m ^T gC _m+1，h_m为第m层的厚度.将式(17)，(18)联立，可以解得

3 边界条件

研究问题是轴对称波场，除了速度模型具有轴对称特点外，震源能量辐射也应该满足轴对称的特点.微地震震源一般近似为能量均匀辐射的点震源，满足研究问题条件.对于微地震集中力源，等效垂直力源可以满足轴对称解条件，但对于水平方向的集中力源轴对称则满足条件，否则不满足条件.这里考虑满足轴对称源条件.

本文采用垂直方向的集中力源，即：

h(x，z)是一个随离开震源点距离呈指数衰减的函数，表达式为

式中：f是震源子波谐波频率，t₀微起始时刻，(x₀，z₀)为震源的位置，α>0为衰减系数. 3.1 边界条件与基函数

基函数采用有限元方法中讨论的HELMIT内插基函数，具体为

第一个节点有两个形函数，G1，G2，第二个节点有两个形函数，G3，G4，表达式为

式中r为径向距离.

定义于单元内部上述基函数，它满足下列条件：

(1)在节点i处，G_i=1，G′ _i=1，在其他节点处，G_i=0，G'_i=0，式中G_i为原函数基函数，G^′_i为原函数导数的基函数.

(2)能保证用它定义的未知量在相邻单元之间的原函数及其一阶导数连续.

(3)应包含任意线性项，使用它定义的单元唯一可满足边界条件.

(4)应满足下列等式：∑G_i=1，∑G^′_i=1∑=1.

3.2 径向截断边界条件处理

在地震波场有限差分数值正演方法中，对边界条件的处理大多采用吸收边界条件.考虑到有限单元的离散和组装特点，径向截断边界条件采无穷单元自然衰减边界条件.具体思路是：在径向有限范围内，离散成一系列的单元，对于截断范围之外，采用一个的较大的单元，根据自然边界条件，r→∞，U→0，保证波场自然衰减.

对于一个单元e，其两个节点是I，I+1单元的长度为L，采用局部坐标，记

已知节点位移，要求节点间任意一点的位移，利用插值方法来计算，公式为

同样得位移的导数为

组成的单元阵为

对于前述的 分别建立具体的单元阵，对各单元阵进行组装形成总矩阵，最终形成广义特征方程组，，解此广义特征方程组便可求出广义特征值和特征向量，把求得的特征值和前面求出的系数A′，A′代入以下方程便可求出一个谐波的位移，再对其进行傅里叶逆变换，最终求得时空域微地震位移数值解，公式为

4 模拟结果分析

为了验证方法的正确性，设计了水平层状模型且为了与一维射线追踪方法对比，设定横向速度不变.以源为轴对称中心，左右分别划分10个单元，每个单元横向大小为20 m.源放在检波器排列的上方埋深在1950 m处，坐标原点在检波器的位置.两种方法的走时对比结果如表 1所示.从结果看到，两种方法计算的走时是一致的.图 2为模拟的波形结果，图中第一个同相轴为P波，第二个同相轴为S波.随着观测点离开源位置的变化，逐渐出现了较弱的转换波.结果说明正演模拟方法是正确的，模拟结果是可靠的.

分析图 2各个结果图，垂直排列检波器记录微地震波的响应变化，随着离震源的不同距离，其变化特征明显不同，在震源附近，地震波接近垂直入射，没有转换波，PS波比较清楚，PS波的能量变化随着垂直排列检波器，由近到远能量逐渐增大.随着离开源的距离增大，出现转换波，这时P波的转换波和原始的S波干涉，使的地震记录同相轴零乱. 当距离很远时，只能看到S波和其转换波.当距离很近时，只有PS波.在源的上方P波顶能量较强，随着离开源的距离增大，能量逐渐衰减.在P波以后的区域波场相对复杂，在源的上方S波较单一，S波的极性为反极性.随着离开源距离的增大，S波和P波的转换波干涉，使的波形发生畸变，随着距离的依次增大，S波和转换的S波能量逐渐增大，之后又逐渐衰减，随着距离的再次增大，S波同相轴向上移动，出现折射S波.

5 结 论

5.1     对于微地震波动方程，在轴对称速度模型和震源均匀辐射情况下，在高频条件下，可以利用分离变量法把波动方程进行分解为垂直方向和径向方向两个方程.为了求得方程的解，利用模式匹配理论，把解微分方程的模式匹配理论引用到本研究中，研究形成了水平轴对称层状模型下，纵向采用解析方法横向采用有限元的混合方法.把计算范围截断边界条件，采取无穷大单元自然衰减边界条件的办法，使计算边界的响应变得较小.通过理论和实际论证和验证了研究方法的正确合理性，并和基于射线追踪的走时反演结果进行了对比分析，证明结果精度满足实际要求.

5.2     微地震波场的响应与空间距离关系是非线性的，本文在算法研究时，采用单元两节点线性插值计算，为了进一步提高结果的精度，建议再进行研究时，采用单元三节点二次插值方法.

致 谢 感谢审稿专家和编辑部的大力支持.