CN103926620A - 基于阵列反褶积处理的水力压裂监测方法 - Google Patents

Abstract

本发明应用于页岩气、致密砂岩气、地热等地下资源的开发，提供一种监测由水力压裂导致的岩石破裂的方法，包括：用阵列反褶积处理，得到微地震震源机制相关散射系数和震源到各检波器的格林函数；将获得的格林函数延拓或者偏移，得到发震时间段的地下监测区域伪波场；通过该波场最大能量的变化拾取微地震震源空间位置；根据震源位置和震源机制相关散射系数搜索震源机制的矩张量，可用于分析震源机制。本发明可获得压裂产生的微地震的位置和震源机制。

Description

基于阵列反褶积处理的水力压裂监测方法

技术领域

本发明涉及地震勘探领域，应用于石油天然气和地热资源的勘探开发。具体地讲，本发明用于监测由水力压裂岩石形成的微地震，测量岩石破裂的空间位置和破裂机制等参数。 

背景技术

近二十多年来，随着人类对能源需求的增长和科学技术的不断发展，水力压裂方法已经广泛的应用于页岩气、致密砂岩气、地热、等地下资源的勘探开发。水力压裂是指将高压流体(水、二氧化碳或其他混合流体)压入富含资源的地下岩层，使的岩石破裂扩大岩石原有的裂隙系统或形成新的裂隙系统，以便资源的析出、释放。这些岩石破裂即定义为微地震。相对天然地震，这种由压裂导致岩石破裂的能量和岩石的位移都很小，但是这些能量仍然能被地震检波器监测到。微地震发生的时间空间顺序，能量大小，震源机制等信息是估计裂隙系统改造，压裂施工等技术环节的重要参考指标。 

地震学的方法技术是实现微地震监测的唯一有效手段，基本原理是：岩石破裂发生的地震波被布置在地表或监测井中的地震检波器观测到；结合监测区域的地震波速度结构，通过分析记录到的波形特征，得到微地震的震源位置和震源机制。与传统勘探地震学相比，微地震是非主动源，其发震时刻与位置都未知。虽然天然地震有非常成熟的方法求解震源位置和震源机制，但是，由于微地震监测的高精度要求，天然地震学的方法无法直接用于微地震监测。 

微地震监测有两个最基本的目的：一是确定微地震发生的空间位置，又称为微地震定位，二是震源的的破裂机制(张裂，走滑，还是闭合)。目前，微地震定位的技术方法主要有两个，有各自的优点也有一些缺陷。一种方法是拾取微地震首波在各检波器的到时，根据给定的速度模型，反演震源位置。这种方法非常依赖到时拾取的精度，通常是人工拾取，不便实现自动化定位。另外一种方法是将检波器记录的地震波向监测区域延拓或者偏移，拾取能量最大位置作为震源。这种方法对爆炸机制的源的定位非常有效，而微地震震源能大多是双力耦(double cuple)或非线性的，导致这种方法有一定的局限性。 

本发明提供一种方法，包括用阵列反褶积(array deconvolution)测量震源到检波器的格林函数(Green function)和震源机制的散射系数(Radiation parameter)的，并用这些测量结果计算微地震的震源位置和震源机制。 

发明内容

本发明提供的阵列反褶积处理技术是用来测量震源到检波器的格林函数和震源机制的散射系数，其中包括：选定包含微地震事件的时窗；截取各检波器对应时窗的记录；计算各记录的功率谱；功率谱归一后求的台阵平均功率谱估计；计算地震子波反褶积滤波因子；用各检波器的地震记录褶积反褶积滤波因子；从而估计格林函数和散射因子。本发明的核心技术是阵列反褶积处理提取格林函数和散射因子，其理论基础如下。 

如图1所示，微地震(02)发生在r_S，其矩张量为M，其P波首波被位于r_R检波器(01)接收到，则波形为： 

$u_k^P=\left(\frac{M_0}{{4\mathrm{\π\ρ\α}}^{3}r}{\mathrm{\γ}}_k\mathrm{\Σ}{\mathrm{\γ}}_l{\mathrm{\γ}}_j{m}_{\mathrm{lj}}\right)g(r_s,t_0;r_R,t)*w\left(t\right)+n\left(t\right);$

(l＝N，E，Z；j＝N，E，Z；k＝N，E，Z)其中，M₀是微地震震级，ρ是岩石密度，α是P波波速，r是地震波从震源传到检波器传播的距离，γ是角度因子(γ_N＝sinicosφ；γ_E＝sinisinφ；γ_Z＝cosi)，g(r_S，t₀；r_R，t)是从震源到检波器的格林函数，w(t)是地震子波，*是褶积运算符号，n(t)是随机噪声。对 第i个检波器接收到第n个微地震信号的垂直分量Z，波形可简写为： 

$u_{\mathrm{Zin}}^P=C_n{S}_{\mathrm{in}}{g}_{\mathrm{in}}\left(t\right)*w_n\left(t\right)+n\left(t\right),$ 其中， $C_n=\frac{M_0}{{4\mathrm{\π\ρ\α}}^3}$ 是场地、震级相关的因子， 是散射因子，格林函数 $g_{\mathrm{in}}\left(t\right)=\left\{\begin{array}{cc}1& t=T_{\mathrm{in}}\\ 0& t\≠T_{\mathrm{in}}\end{array}\right.$ 近似为脉冲函数，T_in是从震源到检波器的地震波走时。的功率谱为： 

$U_{\mathrm{Zin}}^P{U_{\mathrm{Zin}}^P}^{*}={\left(C_n{S}_{\mathrm{in}}\right)}^2\·\left[G_{\mathrm{in}}\left(t\right){G_{\mathrm{in}}}^{*}\left(t\right)\right]\·\left[W_n\left(t\right){W_n}^{*}\left(t\right)\right]+{\mathrm{\σ}}^2,$ 函数U、G、W分别是函数u、g、w的傅立叶变换，上角标*表示共轭，σ²是噪声功率。对观测台阵(I个检波器)，我们通过对检波器记录功率谱归一化后求平均，获得对微地震子波功率谱的平均估计：N_n是归一因子。采用温纳反褶积、最小二乘反褶积或其他反褶积方法计算[参考文献yilmaz，2001])，将估计出微地震子波归一化后的功率谱演算出地震子波w(t)的反褶积滤波器w′(t)，其中w(t)*w′(t)＝N_nδ(0)。推导得到其中序列g_in′(t)中g_in′(m)的绝对值满足：进一步可以导出，g_in(t)＝g_in′(t)/g_in′(m)，C_nN_nS_in＝g_in′(m)。求得g_in(t)和C_nN_nS_in分别为从震源到检波器的格林函数和散射因子参数。用类似的方法，将检波器波形记录的自相关函数归一化后求平均，也可以求得地震子波的反褶积滤波因子w′(t)。 

用从震源到检波器的格林函数，就可以定位震源位置，这是微地震监测中最重要的一个指标。受天然地震学方法研究大地震断层破裂过程成像的技术和勘探地震学地震波偏移成像技术的启发，本发明的微地震定位方法包括：把阵列反褶积方法得到的各检波器的格林函数作为地表波场(即地表伪地震波场，g_in(t))；以监测区域可靠的地震波速度结构为基础，将该地表伪地震波场延拓或者偏移；获得微地震震源附近的地震波场F_n(r，τ)，通过F_n(r，τ)振幅的最大值确定发震时刻和震源位置。下面以kirechoff偏移为例，说明定位原理。 

Kirchhoff积分：在点震源且仅考虑远场效应的情况下，曲面Ω上波场F_Ω(r，t)通过以下积分，可以得到震源S的波场其中t_TR是r_S到r_Ω的走时。我们用g_in(t+T_in)代替δ(t-t_TR)*F′(r_Ω，t)，得到离散时空下，微地震在震源波场表达式为： 即在震源附近的波场，震源位置发震时刻时的波场振幅最大，其他非震源位置的最大振幅都比震源的最大振幅小，这是本发明确定发震时刻和震源位置的理论根据。 

微地震的震源机制是微地震监测另外一个重要的指标，可确定应力方向，破裂机制，推测破裂位置的大致断面方向，是确定地下裂隙系统的重要参数。下面介绍本发明提供的利用阵列反褶积方法测量的g_in′(m)来求解震源机制的原理方法： 

理论上且归一化后震源的矩张量有六个独立变量m_NN，m_NE，m_NZ，m_EE，m_EZ，m_ZZ，设定这六个量的数值范围是[-1，1]。我们通过这六个量的不同数值的组合，计算理论计算的散射系数，然后比较理论计算与实际观测是否一致性，最终确定六个变量的最佳值。设理论计算的散射因子为： 

$\begin{array}{c}{\mathrm{RA}}_{\mathrm{in}}=\frac{1}{r_i}{\mathrm{\γ}}_Z\·\mathrm{\Σ}{\mathrm{\γ}}_i{\mathrm{\γ}}_j{m}_{\mathrm{ij}}\\ =\frac{1}{r_{\mathrm{in}}}\cos i_i\·\begin{array}{}\end{array}\left\{\begin{array}{c}{\sin }^2{i}_i[{\cos }^2{\mathrm{\φ}}_i{m}_{\mathrm{NN}}+{\sin }^2{\mathrm{\φ}}_i{m}_{\mathrm{EE}}+\sin \left({2\mathrm{\φ}}_i\right)m_{\mathrm{NE}}]\\ +{\cos }^2{i}_i{m}_{\mathrm{ZZ}}+\sin \left({2i}_i\right)\cos {\mathrm{\φ}}_i{m}_{\mathrm{NZ}}+\sin {\mathrm{\φ}}_{\mathrm{in}}{m}_{\mathrm{EZ}}\end{array}\right.\end{array}$

定义用理论计算值与观测值的比值表示二者的一致性特征数值： 

COH_in(m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ)＝RA_in(m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ)/g_in′(m) 

当m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ正好是震源矩张量的解时，对每个检波器的一致性特征的数值应当相同，都等于COH_in(m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ)＝1/(C_nN_n)，这组共I个检波器的一致性特征数值的方差最小，即m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ为微地震震源机制最佳解。 

附图说明

附图1，微地震震源、P波射线路径、检波器之间几何关系。 

附图2，基于阵列反褶积技术的水力压裂动态监测步骤；白底框为数据模块，灰色底为处理模块。 

附图3，阵列反褶积事实步骤。 

附图4，计算震源矩张量最佳解步骤。 

具体实施方式

现在对本发明进行清楚、详细的描述。虽然，所描述的实施例并非全部的步骤，但是，根据一下描述的实施例，本领域普通技术人员可以在没有做出创造性劳动的情况下完成实微地震监测工作。施步骤如附图2、3、4，其中相同的编号始终表示相同的数据模块或处理模块。本发明的具体实施例，包括以下步骤： 

步骤(I)，对原始地震记录(101)进行常规处理(102)：用包括去零漂、去直流、去线性趋势，和滤波等常规数字信号处理技术处理波形数据，以提高地震数据的信噪比；然后，用STA/LTA短长时窗平均能量比的方法[参考文献Freiberger,W.F.(1963)]、波形互相关方法[参考文献van Trees，H.L.(1968)]、或其他方法检测本时窗是否含有微地震；如果在该时窗内没有微地震事件，则到下一段地震记录重复步骤(I)；如果有微地震事件，则按照统一时窗T(t_n,t_n+1)截取包含微地震波形数据的I个记录(103)。特别说明：如果微地震监测采用的是井中监测系统，则还要采用偏振分析方法[参考文献Kanasewich，E.R.(1973)]三分量 地震记录转换到P-SV-SH坐标系，然后将3D定位问题转化为更简单的2D定位问题，再按照步骤(1)开始实施。 

步骤(II)，对微地震波形数据(103)进行阵列反褶积处理(104)，可按下面操作实施： 

a)计算时窗T(t_n,t_n+1)波形数据的功率谱(201)，计算方法包括直接用的傅立叶变换和该变换的共轭相乘，也可以采用最小二乘，最大熵等方法来估计功率谱； 

b)利用归一化的(201)计算地震子波的平均估计(202)，其中  是归一化因子，地震子波的平均估计的表达式是  $N_n^2\left[W_n\left(t\right){W_n}^{*}\left(t\right)\right]\≈\underset{i=1}{\overset{1}{\mathrm{\Σ}}}\frac{N_n^2\left(U_{\mathrm{Zin}}^P{U_{\mathrm{Zin}}^P}^{*}\right)}{I{\left(C_n{S}_{\mathrm{in}}\right)}^2};$

c)用地震子波的平均估计(202)来计算反褶积滤波因子(203)w′(t)；计算方法可以采用最小二乘反褶积、维纳反褶积或其他反褶积方法，可参考公开文献中详细阐述[参考文献Yilmaz，2001]； 

d)用反褶积滤波因子(203)褶积微地震波形数据(103)得到未归一的格林函数g_in′(t)(204)：  ${g_{\mathrm{in}}}^{\′}\left(t\right)=u_{\mathrm{Zin}}^P*w^{\′}\left(t\right);$

e)求未归一的格林函数(204)的归一化因子作为散射因子(105)：|g_in′(m)|＝Max(g_in′-(t)|)； 

f)求从震源r_Sn到r_Ri第i个检波器的格林函数(106)：g_in(t)＝g_in′(t)/g_in′(m)； 

步骤(III)：推导T(t_n,t_n+1)监测区域的伪地震波场：已经有非常多公开发表的论著[参考文献Biondi B.，(2006)；Fichtner A.，(2011)；Baysal E，et al.，(1983)]详细描述了地震波场的偏移和延拓理论原理和实施步骤，例如，kirchhoff偏移，有限差分偏移，逆时偏移等等。本发明就以kirchhoff偏移为例说明如何计算监测区域伪地震波场： 

g)根据监测区域的P波速度结构，对监测区域建立网格化的P波速度结构(107)V_P(x，y，z)，其中，P波速度结构还可包含P波的各向异性快波方向θ(x，y，z)和各向异性强度ρ(x，y，z)的等信息； 

h)计算各网格点r(x_P，y_P，z_P)到各检波器r_Ri的距离(03)r_P，Ri＝|r_Ri-r(x_P，y_P，z_P)|，地震射线倾角(05) $i_{P,\mathrm{Ri}}=\arccos \left(\frac{(x_P,y_P)-(x_{\mathrm{Ri}},y_{\mathrm{Ri}})}{r_{\mathrm{Ri}}-r(x_P,y_P,z_P)}\right)$ 和方位角(04) ${\mathrm{\φ}}_{P,\mathrm{Ri}}=\arccos \left(\frac{x_P-x_{\mathrm{Ri}}}{(x_P,y_p)-(x_{\mathrm{Ri}},y_{\mathrm{Ri}})}\right);$

i)计算各网格点r(x_P，y_P，z_P)到各检波器r_Ri的地震波走时T_Ri，P，计算方法包括根据Snell定律做射线追踪、有限差分、或FMM等方法，这些计算走时的方法也有大量公开发表[参考文献Yilmaz，(2001)；Sethian,(1999ab)；Zhao，(2005)]； 

j)假设网格点r(x_P，y_P，z_P)是微地震震源，根据kirchhoff积分原理和发明内容[0007]段到[0012] 段阐述的技术思路，时间段T(t_n,t_n+1)该网格点的伪地震波场振幅为： 

$F(r_P,t)=\frac{\mathrm{\Δx}\·\mathrm{\Δy}}{4\mathrm{\π}}\mathrm{\Σ}\frac{{\cos i}_{P,\mathrm{Ri}}}{{\mathrm{vr}}_{P,\mathrm{Ri}}}{g}_{\mathrm{in}}(t+T_{\mathrm{Ri},P})$

k)对监测区域每一个网格都按照步骤j)实施，就得到了监测区域的伪地震波场(110)F(r，t)； 

步骤(IX)，对监测区域伪地震波场(110)进行振幅时空变化特征的分析(111)；微地震的发震位置在发震时刻的振幅最大，即比较每个时刻每个网格点F(r，t)的振幅拾取最大值Max(|F(r，t)|)，对应的时刻和空间位置即震源发震时刻(112)t_Sn和震源位置(113)r_Sn。 

步骤(X)，采用搜索扫描的方法计算地震矩张量的最佳解(114)： 

L)利用震源空间位置(113)r_Sn和检波器的空间位置r_Ri计算从震源到检波器地震射线几何特征，传播距离、倾角和方位角： 

r_in＝|r_Ri-r_Sn| 

$i_{\mathrm{in}}=\arccos \left(\frac{(x_n,y_n)-(x_{\mathrm{Ri}},y_{\mathrm{Ri}})}{r_{\mathrm{Sn}}-r_{\mathrm{Ri}}}\right)$

${\mathrm{\φ}}_{\mathrm{in}}=\arccos \left(\frac{x_n-x_{\mathrm{Ri}}}{(x_n,y_n)-(x_{\mathrm{Ri}},y_{\mathrm{Ri}})}\right)$

m)假设震源矩张量M的六个独立变量分别为m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ，则对第i个检波器其理论散射系数为： 

$\begin{array}{c}{\mathrm{RA}}_{\mathrm{in}}(m_{\mathrm{NN}}^{\′},m_{\mathrm{NE}}^{\′},m_{\mathrm{NZ}}^{\′},m_{\mathrm{EE}}^{\′},m_{\mathrm{EZ}}^{\′},m_{\mathrm{ZZ}}^{\′})\\ =\frac{1}{r_{\mathrm{in}}}\cos i_{\mathrm{in}}\·\begin{array}{}\end{array}\left\{\begin{array}{c}{\sin }^2{i}_{\mathrm{in}}[{\cos }^2{\mathrm{\φ}}_{\mathrm{in}}{m}_{\mathrm{NN}}^{\′}+{\sin }^2{\mathrm{\φ}}_{\mathrm{in}}{m}_{\mathrm{EE}}^{\′}+\sin \left({2\mathrm{\φ}}_{\mathrm{in}}\right)m_{\mathrm{NE}}^{\′}]\\ +{\cos }^2{i}_{\mathrm{in}}{m}_{\mathrm{ZZ}}^{\′}+\sin \left({2i}_{\mathrm{in}}\right)\cos {\mathrm{\φ}}_{\mathrm{in}}{m}_{\mathrm{NZ}}^{\′}+\sin {\mathrm{\φ}}_{\mathrm{in}}{m}_{\mathrm{EZ}}^{\′}\end{array}\right.\end{array}$

一致性特征数值为： 

COH_in(m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ＝RA_in(m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ)/g_in′(m) 

对格点m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ的I个检波器的一致性特征数值方差为： 

$\begin{array}{c}{\mathrm{\σ}}_n^2(m_{\mathrm{NN}}^{\′},m_{\mathrm{NE}}^{\′},m_{\mathrm{NZ}}^{\′},m_{\mathrm{EE}}^{\′},m_{\mathrm{EZ}}^{\′},m_{\mathrm{ZZ}}^{\′})\\ =\underset{b=1}{\overset{I}{\mathrm{\Σ}}}{\left[\begin{array}{c}{\mathrm{COH}}_{\mathrm{bn}}(m_{\mathrm{NN}}^{\′},m_{\mathrm{NE}}^{\′},m_{\mathrm{NZ}}^{\′},m_{\mathrm{EE}}^{\′},m_{\mathrm{EZ}}^{\′},m_{\mathrm{ZZ}}^{\′})\\ -\frac{1}{I}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{\mathrm{COH}}_{\mathrm{in}}(m_{\mathrm{NN}}^{\′},m_{\mathrm{NE}}^{\′},m_{\mathrm{NZ}}^{\′},m_{\mathrm{EE}}^{\′},m_{\mathrm{ZZ}}^{\′})\end{array}\right]}^2\end{array}$

我们将以m_NN，m_NE，m_NZ，m_EE，m_EZ，m_ZZ张成六维网格化空间(-1≤m_ij≤1)，例如以0.1为网格间隔，对每一个格点计算该格点的I个检波器的一致性特征数值的方差，形成一个六维的一致性特征 数值方差的数据体(302)； 

n)搜索一致性特征数值方差的数据体(302)中方差最小的值，对应的m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ即地震矩张量的最佳解(115)。 

步骤(XI)，选取下一个时窗T(t_n+1，t_n+2)，重复进行上述步骤(I)至步骤(XI)操作。 

参考文献： 

Vallina，A.U.(1999)，Principles ofSeismology[M]，Cambridge University Press，ISBN9780521624787 

Ozdogan Yilmaz，Stephen M.Doherty,(2001)，Seismic data analysis：processing，inversion，and interpretation of seismic data[M]，Society of Exploration Geophysicists，ISBN1560800941，978156080094] 

Freiberger,W.F.(1963).An approximation method in signal detection[J]，Quart.J.Appl.Math.20，373-378. 

van Trees，H.L.(1968).Detection，Estimation and Modulation Theory[J]，John Wiley and Sons，Inc.，New York，249pp. 

Kanasewich，E.R.(1973)，Time Sequence Analysis in Geophysics[M]，Univ.Alberta Press. 

Biondi B.，(2006)，3D seismic imaging[M].Society of Exploration Geophysicists. 

Fichtner A.，(2011)，Full seismic waveform modelling and inversion[M].Berlin/Heidelberg：Springer. 

Baysal E，Kosloff D D，Sherwood J W C.，(1983)，Reverse time migration[J].Geophysics，48(11)：1514-1524. 

Sethian，J.A.，(1999a).Level Set Methods and Fast Marching Methods：EvolvingInterfaces in Computational Geometry,Fluid Mechanics，Computer Vision，and Materials Science[M]，Cambridge University Press. 

Sethian，J.A.&Popovici，A.M.，(1999b).Three dimensional traveltimes computation using the Fast Marching Method[J]，Geophysics，64，516-523. 

Zhao，H.K.，(2005)，A fast sweeping method for eikonal equations[J]，Mathematics of Computation，74，603—628。

Claims (8)

1.一种基于阵列反褶积处理的水力压裂监测方法，包括以下步骤： 

(I)对监测台网的数字地震记录进行常规的处理包括滤波、去零漂等，获得高信噪比记录，再应用STA/LTA或波形互相关方法等方法，截取含有微地震事件波形的记录； 

(II)采阵列反褶积方法处理从各检波器截取的含有微地震波形的记录，得到微地震震源到各检波器的散射因子和归一化后的格林函数； 

(III)以步骤(II)测量得到的归一化后的格林函数作为伪地震波场，采用地震波场的偏移或延拓等方法，计算得到微地震发生时间段震源附近的地震波场； 

(IX)根据计算得到的微地震发生时间段震源附近的地震波场，搜索最大振幅，拾取发震时刻和震源位置； 

(X)从地震矩张量解的理论取值范围中逐个计算理论的散射因子，并与步骤(II)中用阵列反褶积处理得到的散射因子进行比较，通过理论计算值与实际观测值之间的一致性，扫描搜索微地震矩张量的最佳解； 

(XI)对下一段微地震记录重复进行步骤(I)至步骤(XI)的操作。 

2.根据权利要求1所述的一种基于阵列反褶积处理的水力压裂监测方法，其中，步骤(I)包括以下特征：用包括去零漂、去直流、去线性趋势、滤波等常规数字信号处理技术处理波形数据，以提高地震数据的信噪比； 

用STA/LTA短长时窗平均能量比的方法、波形互相关方法、或其他方法检测记录中是否含有微地震事件；如果在该时窗内不存在微地震事件，则到下一段地震记录重复步骤(I)； 

如果存在微地震事件，以统一的时窗T(t_n,t_n+1)，截取完整微地震波形，得到I个微地震波形数据； 

特别说明，如果微地震监测仅采用的是井中监测系统，则还要采用偏振分析方法三分量地震记录转换到P-SV-SH坐标系，即将3D定位问题转化为更简单的在地震波传播平面内的2D定位问题，再按照步骤(1)开始实施。 

3.根据权利要求1所述的一种基于阵列反褶积处理的水力压裂监测方法，其中，步骤(II)包括以下特征：利用各个检波器波形记录的功率谱或者自相关函数来估计微地震地震子波的归一化功率谱或者归一化自相关函数；采用最小二乘反褶积、维纳反褶积、或其他反褶积方法，计算反滤波因子；将反滤波因子褶积各检波器记录；该褶积结果的归一化因子即散射系数，归一化后的结果即从震源到检波器的格林函数。 

4.根据专利要求1所述方法的步骤(II)或专利要求3所述特征，其中，以各检波器波形功率谱来计算散射系数和格林函数的包括以下步骤： 

a)采用法包括直接用波形的傅立叶变换和该变换的共轭相乘、最小二乘法或其他方法来估计时窗T(t_n,t_n+1)波形数据的功率谱

b)将归一化的功率谱求平均，可计算得到地震子波的平均估计其中  是归一化因子，地震子波的平均估计的表达式是 

c)采用最小二乘反褶积、维纳反褶积或其他反褶积方法，用地震子波的平均估计来计算反褶积滤波因子w′(t)； 

d)用反褶积滤波因子w′(t)褶积微地震波形数据得到未归一的格林函数g_in′(t)； 

e)求未归一的格林函数g_in′(t)的归一化因子：|g_in′(m)|＝Max(|g_in′(t)|)； 

f)求从震源r_Sn到r_Ri第i个检波器的格林函数：g_in(t)＝g_in′(t)/g_in′(m)。 

5.根据权利要求1所述的一种基于阵列反褶积处理的水力压裂监测方法，其中，步骤(III)包括以下特征： 

将权利要求1所述方法中步骤(II)计算得到的从震源到各检波器的格林函数作为伪地震波场； 

采用地震波场的偏移和延拓方法，包括Kirchhoff偏移，有限差分偏移，逆时偏移等偏移方法，将伪地震波场偏移到监测区域，获得微地震发生时间段在震源附近的地震波场。 

6.根据权利要求1所述方法的步骤(III)和权利要求5所述方法的，其中，以Kirchhoff偏移方法为例，包括以下特征： 

g)以均匀或非均匀网格建立对监测区域P波速度结构，或者建立包含各向异性特征的P波速的结构； 

h)计算各网格点r(x_P，y_P，z_P)到各检波器r_Ri的距离r_P，Ri，地震射线倾角i_P，Ri和方位角φ_P，Ri； 

i)采用包括Snell定律做射线追踪，用有限差分、或FMM等方法，计算各网格点r(x_P，y_P，z_P)到各检波器r_Ri的地震波走时T_Ri，P； 

j)假设网格点r(x_P，y_P，z_P)是微地震震源，在时间段T(t_n,t_n+1)内，该网格点的伪地震波场振幅变化为 

k)对监测区域每一个网格都按照步骤j)实施，就得到了监测区域的伪地震波场F(r，t)。

7.根据权利要求1所述的一种基于阵列反褶积处理的水力压裂监测方法，其中，步骤(IX)包括以下特 征： 

将权利要求1所述步骤(III)或权利要求4或权利要求5所述方法获得的监测区域的伪地震波场F(r，t)，进行包括逐个网格点比较等方法获得最大值Max(|F(r，t)|)，该最大值对应的时刻和空间位置即震源发震时刻t_Sn和震源位置r_Sn。 

8.根据权利要求1所述的一种基于阵列反褶积处理的水力压裂监测方法，其中，步骤(X)包括以下特征：将以m_NN，m_NE，m_NZ，m_EE，m_EZ，m_ZZ张成六维网格化空间，包括等间距或非等间距，m_ij的取值区间可以选择[-1，1]或者其他取值范围； 

对每一个格点，例如m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ，计算第i个检波器其理论散射系数，包括公式 或者其他近似表达式，其中r_in、i_in、φ_in分别是地震波从按照步骤(IX)从拾取的震源位置到第i个检波器的传播距离、地震射线倾角、和方位角； 

以检波器其理论散射系数与对应检波器观测值g_in′(m)比值，或者其他表达式作为表示理论计算值与实际观测值一致性的函数； 

对每一个格点，例如m′_NN，m′_NE，m′_NZ，m′_EE，m′_EZ，m′_ZZ，计算所有检波器的一致性特征数值，得到这个格点的一组一致性特征数值； 

对每一个格点，计算该格点那一组一致性特征数值的的方差或者其他统计特征，并用该统计特征的数值构成一个六维的表示理论计算值与实际观测值一致的数据体； 

根据统计特征，例如方差，搜索分析一致数据体中最小方差，对应的格点坐标即微地震矩张量的最佳解。