CN102879801B - 一种基于射孔约束的EnKF微地震事件位置反演方法 - Google Patents

Abstract

本发明公开了一种基于射孔约束的EnKF微地震事件位置反演方法，该方法包括：(1)建立速度场；(2)选择第k个事件点；(3)通过求解程函方程得到第j个检波器接收到的该事件点的初至旅行时；(4)应用集合卡尔曼滤波方法将步骤(3)中计算得到的初至旅行时与第j个检波器实际接收到的初至旅行时进行拟合；(5)检查每一个检波器所接收到的该事件点的旅行时是否全部拟合，如果检查结果为否，则返回到步骤(3)，并使变量j增1，以对第j+1个检波器进行拟合，如果检查结果为是，则完成事件定位。

Description

一种基于射孔约束的EnKF微地震事件位置反演方法

技术领域

本发明属于石油天然气地震勘探领域，主要应用于石油地震勘探的微地震监测，具体地说，本发明涉及一种能够更精确地确定微地震事件位置的反演方法。

背景技术

在微地震压裂中，我们可以得到每一个检波器接收到的每一个事件的P波和S波的初至旅行时。通过三分量检波器，事件的方位角可以从矢量图中得到。P波和S波旅行时包含了从事件位置到检波器之间的距离信息，比如事件发生时间。由矢量图得到的方位角给出了事件位置的方位信息。这里我们用向量d_obs来代表测量数据向量，包括每一检波器接收到的每一事件的纵横波旅行时和方位角。

在微震监测中需要确定的参数包括每个事件的位置(x,y,z)和发生时间(t₀)。假设有Nr个检波器和N个事件点，则参数向量X的维数为4×N，如X＝[x_i,y_i,z_i,t_0i]^T，其中i＝1,2,…,N。如果P波和S波的旅行时以及方位角都是可用的，那么测量数据向量最大维数为3×Nr，如d＝[P_j,S_j,A_j]^T,其中j＝1,2,…,Nr。

以前的文献中，地震定位是基于位置参数和不同接收点的初至旅行时之间是线性关系的，其表达形式如下面的式1：

A^kΔX^k+1+S^k+1＝ΔT^k…………………………………………式1

其中T是剩余旅行时，即对每一个接收器计算得到的旅行时与观测到的旅行时之间的差异，可能包括P波或者S波旅行时之一，或者两者都包括。向量X包括所有事件的位置参数(x,y,z)和事件的起始时间t₀。矩阵A是敏感性矩阵，由第k次迭代预测得到的事件位置参数计算的旅行时得到。向量S包含由错误的速度结构引起的校正参数，向量S可以用不同的算法来处理得到，比如静态位置关系法(ST)，特定震源位置关系法(SSST)。当向量S是用静态位置关系法得到时，每一个检波器接收的信号s值与事件位置无关。当向量S是用特定震源位置关系法得到时，校正关系是接收器位置与事件位置间的一个方程。在两种相关定位方法中，都至少需要用迭代方法来求解X_K+1和S_K+1中的一个。每一次迭代，事件位置参数都用下面加入了校正关系S的式2进行更新：

X^k+1＝X^k+ΔX^k+1………………………………………………式2

用更新后的新的位置向量X，我们进行正演模型得到旅行时，重新计算T。S_k+1通过求解方程1得到。

上述的校正关系可以减少由不精确的速度模型反演得到的微地震事件位置的不确定性范围。然而，其也对更多参数进行优化，因此增加了反演问题的不适定性，例如，其使最优化问题趋于拥有多解性。因此上述方法可能会使计算得到的微地震事件位置远离真实的事件位置，即异常值。

发明内容

本发明的一方面提供了一种能够减少传统定位方法的不确定性，以更精确地确定微地震事件位置的基于射孔约束的EnKF微地震事件位置反演方法。

根据本发明的一方面，一种基于射孔约束的EnKF微地震事件位置反演方法包括下述步骤：(1)建立速度场；(2)选择第k个事件点；(3)通过求解程函方程得到第j个检波器(其中，j＝1，…，Nr，Nr表示检波器的总个数)接收到的该事件点的初至旅行时；(4)应用集合卡尔曼滤波方法将步骤(3)中计算得到的初至旅行时与第j个检波器实际接收到的初至旅行时进行拟合；(5)检查每一个检波器所接收到的该事件点的旅行时是否全部拟合，如果检查结果为否，则返回到步骤(3)，并使变量j增1，以对第j+1个检波器进行拟合，如果检测结果为是，则完成事件定位。

根据本发明的一方面，可通过计算声波测井和射孔时间数据来建立速度场。

根据本发明的一方面，步骤(4)可包括下述步骤：在贝叶斯框架下定义目标函数；应用集合卡尔曼滤波来进行事件定位，首先随机生成Ne个事件位置的实现；对第j个检波器接收到的微震事件初至旅行时进行扰动，生成Ne个观测数据向量i＝1...Ne；通过拟合第j个检波器接收到的事件的初至旅行时来更新事件位置向量X；使变量j增1，重复上述过程以对第j+1个检波器接收到的事件初至旅行时进行拟合。

根据本发明的一方面，Ne个事件位置的实现可为50-100个事件位置的实现。

根据本发明的一方面，目标函数可由下式表示：

$O\left(X\right)=\frac{1}{2}{(X-\stackrel{\‾}{X^p})}^T{C}_X^{-1}(X-\stackrel{\‾}{X^p})+\frac{1}{2}{(d_{\mathrm{obs}}-d_{\mathrm{calc}})}^T{C}_D^{-1}(d_{\mathrm{obs}}-d_{\mathrm{calc}})$

其中，O(X)表示为确定微地震事件位置所建立的目标函数，通过迭代使本目标函数的变化达到某一收敛条件ε_o(如ε_o＝10^-5)或l(如l＝10)。如l为迭代次数，X表示事件位置向量，表示上次迭代所更新的数据的平均值，第一次计算时为射孔位置参数，C_X是位置参数误差协方差矩阵 $C_X=\left[\begin{array}{cccc}{\mathrm{\δ}}_x^2& 0& 0& 0\\ 0& {\mathrm{\δ}}_y^2& 0& 0\\ 0& 0& {\mathrm{\δ}}_z^2& 0\\ 0& 0& 0& {\mathrm{\δ}}_t^2\end{array}\right],$ 其中是事件位置参数的标准偏差，其由裂缝大小计算得到，d_obs表示观测数据，d_calc表示正演计算得到的预测数据，C_D表示测量误差协方差矩阵，其中是观测数据的标准偏差，其中，k＝1,…,Nobs，Nobs为测量数据维数。

根据本发明的一方面，随机生成Ne个事件位置的实现的步骤可由下式表示：

$X_i^p=\stackrel{\‾}{X^p}+C_X^{1/2}Z,i=1,...,\mathrm{Ne},$

其中，向量Z是(0，1)内的一组服从高斯分布的随机数组成的随机偏移向量。表示上次迭代所更新的数据的平均值，在本方法中，由于是用射孔点约束来确定微地震事件位置，所以在第一次运算时，为射孔位置参数。C_X是位置参数误差协方差矩阵 $C_X=\left[\begin{array}{cccc}{\mathrm{\δ}}_x^2& 0& 0& 0\\ 0& {\mathrm{\δ}}_y^2& 0& 0\\ 0& 0& {\mathrm{\δ}}_z^2& 0\\ 0& 0& 0& {\mathrm{\δ}}_t^2\end{array}\right].$

根据本发明的一方面，拟合第j个检波器接收到的微震事件初至旅行时的步骤可包括：对Ne个实现中的每一个进行微震事件正演模拟的步骤，以得到Ne个预测的数据向量，该步骤由下式表示，

$d_i^p=d\left(X_i^p\right),i=1,...,\mathrm{Ne},$

其中，如上述，为生成的Ne个实现。

然后，对第j个检波器接收到的微震事件初至旅行时进行扰动，生成Ne个观测数据向量，该步骤由下式表示，

$d_{\mathrm{obs},i}^j=d_{\mathrm{obs}}^j+C_D^{1/2}Z,i=1,...,\mathrm{Ne},$

其中，C_D表示测量误差协方差矩阵，向量Z是(0，1)内的一组服从高斯分布的随机数组成的随机偏移向量。

根据本发明的一方面，更新事件位置向量X的步骤可由下式表示：

$X_i^u=X_i^p+C_{X,D}{(C_{D,D}+C_D)}^{-1}(d_{\mathrm{obs},i}^j-d_i^p),i=1,...,\mathrm{Ne},$

其中，表示更新后的事件位置向量，表示上一次更新得到的事件位置向量，C_X,D为事件位置参数与观测数据的协方差矩阵，C_D,D为观测数据与观测数据的协方差矩阵，C_D表示测量误差协方差矩阵。

根据本发明的一方面，协方差矩阵C_X,D和C_D,D可由下式表示：

$C_{X,D}=\frac{1}{N_e-1}\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}(X_i^p-\stackrel{\‾}{X^p}){(d_i^p-\stackrel{\‾}{d^p})}^T,$

$C_{D,D}=\frac{1}{N_e-1}\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}(d_i^p-\stackrel{\‾}{d^p}){(d_i^p-\stackrel{\‾}{d^p})}^T,$

其中，代表Ne个事件位置实现的平均值，代表Ne个预测数据实现的平均值。

根据本发明的一方面，和可分别由下式表示：

$\stackrel{\‾}{X^p}=\frac{1}{N_e}\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}\left(X_i^p\right),$

$\stackrel{\‾}{d^p}=\frac{1}{N_e}\underset{i=1}{\overset{N_e}{\mathrm{\Σ}}}\left(d_i^p\right).$

本发明应用集合卡尔曼滤波(Ensemble Kalman Filter，简称EnKF)方法来确定水力压裂时产生的微地震事件位置。对于微地震事件定位问题，集合卡尔曼滤波将吸收每个检波器接收的初至旅行时数据，并对其进行拟合，以使事件位置达到最优。因此，可以精确地确定微地震事件位置，可以更好地刻画水力压裂尺度和油藏改变体积。

本发明的方法与传统技术相比，最主要的优化之处在于：1)应用水力压裂点的射孔数据来调整微震事件位置参数，这可以极大地减少事件位置出现异常值的可能性。以使得最终得到的事件位置在允许的合理范围内；2)新技术应用有效的集合卡尔曼滤波方法作为优化方法，来解决微地震事件位置反演问题。

具体实施方式

下面将结对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

根据本发明的一个实施例，一种基于射孔约束的EnKF微地震事件位置反演方法包括下述步骤：(1)建立速度场；(2)选择第k个事件点；(3)通过求解程函方程得到第j个检波器(其中，j＝1，…，Nr，Nr表示检波器的总个数)接收到的该事件点的初至旅行时；(4)应用集合卡尔曼滤波方法将步骤(3)中计算得到的初至旅行时与第j个检波器实际接收到的初至旅行时进行拟合；(5)检查每一个检波器所接收到的该事件点的旅行时是否全部拟合，如果检测结果为否，则返回到步骤(3)，并使变量j增1，以对第j+1个检波器进行拟合，如果检测结果为是，则完成事件定位。

其中，速度场可通过计算声波测井和射孔时间数据得到，其具体建立方法和步骤对本领域技术人员来讲是公知的，因此在这里省略对速度场(速度模型)的详细描述。

根据本发明的一个实施例，应用集合卡尔曼滤波方法将计算得到的第j个检波器的初至旅行时与第j个检波器实际接收到的初至旅行时进行拟合的步骤可包括：在贝叶斯框架下定义目标函数；应用集合卡尔曼滤波来进行事件定位，首先需要随机生成Ne个(50-100个)事件位置的实现；为了拟合第j个检波器接收到的微震事件初至旅行时对其进行扰动，生成Ne个观测数据向量i＝1...Ne；通过拟合第j个检波器接收到的事件的初至旅行时来更新事件位置向量X；使变量j增1，重复上述过程以对第j+1个检波器接收到的事件初至旅行时进行拟合。

下面通过具体示例描述应用集合卡尔曼滤波方法将计算得到的第j个检波器的初至旅行时与第j个检波器实际接收到的初至旅行时进行拟合的步骤：

在贝叶斯框架下定义目标函数式3，

$O\left(X\right)=\frac{1}{2}{(X-\stackrel{\‾}{X^p})}^T{C}_X^{-1}(X-\stackrel{\‾}{X^p})+\frac{1}{2}{(d_{\mathrm{obs}}-d_{\mathrm{calc}})}^T{C}_D^{-1}(d_{\mathrm{obs}}-d_{\mathrm{calc}})$ ……………式3

其中，O(X)表示为确定微地震事件位置所建立的目标函数，通过迭代使本目标函数的变化达到某一收敛条件ε_o(如ε_o＝10^-5)或l(如l＝10)，例如l为迭代次数，X表示事件位置向量，表示上次迭代所更新的数据的平均值，第一次计算时为射孔位置参数，C_X是位置参数误差协方差矩阵 $C_X=\left[\begin{array}{cccc}{\mathrm{\δ}}_x^2& 0& 0& 0\\ 0& {\mathrm{\δ}}_y^2& 0& 0\\ 0& 0& {\mathrm{\δ}}_z^2& 0\\ 0& 0& 0& {\mathrm{\δ}}_t^2\end{array}\right],$ 其中是事件位置参数的标准偏差，其可由裂缝大小计算得到，d_obs表示观测数据，d_calc表示正演计算得到的预测数据，C_D表示测量误差协方差矩阵，其中是观测数据的标准偏差，其中,k＝1,…,Nobs，Nobs为测量数据维数。根据本发明的一个实施例，如果P波和S波的旅行时以及方位角都是可用的，那么Nobs最大为3。

式3右侧的第一项是与上一次迭代得到的位置信息的正则化项。第二项是观测数据(d_obs)与计算数据(d_calc)的不匹配性项。通过应用集合卡尔曼滤波方法可以有效的求取上述目标的最小值，来确定事件位置和事件初始时间，这是微地震领域的一个新的定位方法。

应用集合卡尔曼滤波来进行事件定位，具体地说，首先需要随机生成Ne个(大约50-100个)事件位置的实现，

…………………………………式4,

其中，向量Z是(0，1)内的一组服从高斯分布的随机数组成的随机偏移向量。表示上次迭代所更新的数据的平均值，在本方法中，由于是用射孔点约束来确定微地震事件位置，所以在第一次运算时，为射孔位置参数。C_X是位置参数误差协方差矩阵 $C_X=\left[\begin{array}{cccc}{\mathrm{\δ}}_x^2& 0& 0& 0\\ 0& {\mathrm{\δ}}_y^2& 0& 0\\ 0& 0& {\mathrm{\δ}}_z^2& 0\\ 0& 0& 0& {\mathrm{\δ}}_t^2\end{array}\right].$

对于Ne个实现中的每一个，均进行微震事件正演模拟，从而可以得到Ne个预测的数据向量，

………………………………………式5。

其中，如上所述，为生成的Ne个实现。

然后，为了拟合第j个检波器接收到的微震事件初至旅行时对其进行扰动，生成Ne个观测数据向量，

………………………………式6。

其中，C_D表示测量误差协方差矩阵，向量Z是(0，1)内的一组服从高斯分布的随机数组成的随机偏移向量。

随后，事件位置向量X可以通过拟合第j个检波器接收到的事件初至旅行时进行更新，

$X_i^u=X_i^p+C_{X,D}{(C_{D,D}+C_D)}^{-1}(d_{\mathrm{obs},i}^j-d_i^p),i=1,...,\mathrm{Ne}$ ……………式7，

其中，表示更新后的事件位置向量，表示上一次更新得到的事件位置向量，C_X,D为事件位置参数与观测数据的协方差矩阵，C_D,D为观测数据与观测数据的协方差矩阵，C_D表示测量误差协方差矩阵。

在式7中，协方差矩阵C_X,D和C_D,D定义如下：

$C_{X,D}=\frac{1}{N_e-1}\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}(X_i^p-\stackrel{\‾}{X^p}){(d_i^p-\stackrel{\‾}{d^p})}^T$ …………………………式8,

$C_{D,D}=\frac{1}{N_e-1}\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}(d_i^p-\stackrel{\‾}{d^p}){(d_i^p-\stackrel{\‾}{d^p})}^T$ …………………………式9,

其中代表Ne个事件位置实现的平均值，代表Ne个预测数据实现的平均值，其表达式如下：

…………………………………………式10,

…………………………………………式11,

重复上述过程(式5-式11)，并将可以对第j+1个检波器接收到的事件初至旅行时进行拟合。注意：其中在此将要更新为由先前接收器接收到的数据更新后的向量。

重复上述过程(式5-式11)，并由下式对进行更新，

…………………………………………………式12，

可以对第j+1个检波器接收到的事件初至旅行时进行拟合。注意：其中在此将要更新为由先前接收器(第j个接收器)更新后的向量。

下面将描述根据本发明的方法的一个具体示例，将本发明的方法应用到一个有6级水力压裂的水平井中。其速度模型通过计算声波测井和射孔时间数据得到，当射孔枪进行射孔时，12级检波器同时进行接收，得到较好质量的射孔事件初至旅行时数据。用射孔位置进行约束，用集合卡尔曼滤波方法对微震事件位置进行拟合。

虽然已经结合示例性实施例详细描述了本发明，但是本领域技术人员应当理解，在不脱离本发明的精神和范围的情况下，可以对这些实施例进行各种修改和改变，本发明的范围在权利要求书及其等同物中限定。

Claims (6)

1.一种基于射孔约束的EnKF微地震事件位置反演方法，包括下述步骤：

(1)建立速度场；

(2)选择第k个事件点；

(3)通过求解程函方程得到第j个检波器接收到的该事件点的初至旅行时其中，j＝1,...,Nr，Nr表示检波器的总个数；

(4)应用集合卡尔曼滤波方法将步骤(3)中计算得到的初至旅行时与第j个检波器实际接收到的初至旅行时进行拟合；

(5)检查每一个检波器所接收到的该事件点的旅行时是否全部拟合，如果检查结果为否，则返回到步骤(3)，并使变量j增1，以对第j+1个检波器进行拟合，如果检测结果为是，则完成事件定位，

其中，步骤(4)包括下述步骤：

在贝叶斯框架下定义目标函数；

应用集合卡尔曼滤波来进行事件定位，首先随机生成Ne个事件位置的实现，其中，随机生成Ne个事件位置的实现的步骤由i＝1,…,Ne表示，其中，向量Z是(0，1)内的一组服从高斯分布的随机数组成的随机偏移向量，表示上次迭代所更新的数据的平均值，在第一次运算时，为射孔位置参数，C_X是位置参数误差协方差矩阵 $C_X=\left[\begin{array}{cccc}{\mathrm{\δ}}_x^2& 0& 0& 0\\ 0& {\mathrm{\δ}}_y^2& 0& 0\\ 0& 0& {\mathrm{\δ}}_z^2& 0\\ 0& 0& 0& {\mathrm{\δ}}_t^2\end{array}\right];$

对第j个检波器接收到的微震事件初至旅行时进行扰动，生成Ne个观测数据向量其中，i＝1...Ne；

通过拟合第j个检波器接收到的事件的初至旅行时来更新事件位置向量X；

使变量j增1，重复上述过程以对第j+1个检波器接收到的事件初至旅行时进行拟合，

其中，所述目标函数由下式表示：

$O\left(X\right)=\frac{1}{2}{(X-\stackrel{\‾}{X^p})}^T{C}_X^{-1}(X-\stackrel{\‾}{X^p})+\frac{1}{2}{(d_{\mathrm{obs}}-d_{\mathrm{calc}})}^T{C}_D^{-1}(D_{\mathrm{obs}}-d_{\mathrm{calc}}),$

其中，O(X)表示为确定微地震事件位置所建立的目标函数，X表示事件位置向量，表示上次迭代所更新的数据的平均值，在第一次计算时，为射孔位置参数，C_X是位置参数误差协方差矩阵， $C_X=\left[\begin{array}{cccc}{\mathrm{\δ}}_x^2& 0& 0& 0\\ 0& {\mathrm{\δ}}_y^2& 0& 0\\ 0& 0& {\mathrm{\δ}}_z^2& 0\\ 0& 0& 0& {\mathrm{\δ}}_t^2\end{array}\right];$ 其中是事件位置参数的标准偏差，d_obs表示观测数据，d_calc表示正演计算得到的预测数据，C_D表示测量误差协方差矩阵，其中是观测数据的标准偏差，其中，k＝1,…,Nobs，Nobs为测量数据维数，

其中，对第j个检波器接收到的微震事件初至旅行时进行扰动，生成Ne个观测数据向量的步骤包括：

对Ne个实现中的每一个进行微震事件正演模拟，以得到Ne个预测的数据向量，该步骤由下式表示，

$d_i^p=d\left(X_i^p\right),i=1,...,\mathrm{Ne};$

对第j个检波器接收到的微震事件初至旅行时进行扰动，生成Ne个观测数据向量该步骤由下式表示，

$d_{\mathrm{obs},i}^j=d_{\mathrm{obs}}^j+C_D^{1/2}Z,i=1,...,\mathrm{Ne},$

其中，C_D表示测量误差协方差矩阵，向量Z是(0，1)内的一组服从高斯分布的随机数组成的随机偏移向量。

2.如权利要求1所述的方法，其中，通过计算声波测井和射孔时间数据来建立速度场。

3.如权利要求1所述的方法，其中，Ne个事件位置的实现为50-100个事件位置的实现。

4.如权利要求1所述的方法，其中，更新事件位置向量X的步骤由下式表示：

$X_i^u=X_i^p+C_{X,D}{(C_{D,D}+C_D)}^{-1}(d_{\mathrm{obs},i}^j-d_i^p),i=1,...,\mathrm{Ne},$

其中，表示更新后的事件位置向量，表示上一次更新得到的事件位置向量，C_X,D为事件位置参数与观测数据的协方差矩阵，C_D,D为观测数据与观测数据的协方差矩阵，C_D表示测量误差协方差矩阵，为正演得到的Ne个预测数据。

5.如权利要求4所述的方法，其中，协方差矩阵C_X,D和C_D,D由下式表示：

$C_{X,D}=\frac{1}{N_e-1}\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}(X_i^p-\stackrel{\‾}{X^p}){(d_i^p-\stackrel{\‾}{d^p})}^T,$

$C_{D,D}=\frac{1}{N_e-1}\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}(d_i^p-\stackrel{\‾}{d^p}){(d_i^p-\stackrel{\‾}{d^p})}^T,$

其中，代表Ne个事件位置实现的平均值，代表Ne个实现的预测数据的平均值。

6.如权利要求5所述的方法，其中，和分别由下式表示：

$\stackrel{\‾}{X^p}=\frac{1}{N_e}\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}\left(X_i^p\right),$

$\stackrel{\‾}{d^p}=\frac{1}{N_e}\underset{i=1}{\overset{N_e}{\mathrm{\Σ}}}\left(d_i^p\right).$