一种提取水平层状介质各向异性参数的方法与流程

本发明涉及微地震监测数据处理技术，是一种用井中微地震自身的信息提取水平层状介质各向异性参数的计算方法。

背景技术：

微地震压裂监测技术现在已成为致密储层油气田开发中一项常用的监测技术，其最主要的应用是水力压裂裂缝分布的监测，在优化压裂方案和开发井网部署方面起到重要作用。微地震监测有地面监测和井中监测两种方式。地面监测就是在监测目标(如压裂井)周围地面布置若干接收点或测线进行微地震监测。井中监测是在监测目标区域周围临近的井中布置接收检波器排列进行微地震监测。

微地震监测的首要目标是对微地震事件进行定位，而定位精度取决于模型的准确。井中微地震一般采用水平层状介质模型。由于观测井离压裂井较近(水平距离约为200米～600米)，并且受压裂对地层破坏的影响，井中微地震更易受各向异性因素影响。有些地区受条件限制(无地面三维地震资料、无正交偶极子声波测井)不能提供各向异性参数。

因此，本领域亟需一种提取水平层状介质各向异性参数的方法。

技术实现要素：

为了解决现有技术中存在的问题，本发明提出了一种利用井中微地震信息提取水平层状介质各向异性参数的方法。该方法只利用井中微地震自身的信息提取水平层状介质各向异性参数，特别适合没有三维地震和正交偶极子声波测井的地区。

根据本发明的一个方面，提供一种提取水平层状介质各向异性参数的方法，该方法包括：利用旅行时残差、射线路径长度向量、群慢度与各层中各向异性参数构造超定方程；求解超定方程获得各向异性参数。

进一步地，该方法在求解超定方程获得各向异性参数后获得修正量；判断是否满足模型的修正量充分小，不满足条件则更新模型，在初始模型各个参数加上所获得的修正量，然后再重新构造方程和求解；满足条件则输出模型。

进一步地，所述超定方程为：

式中为旅行时残差，δτⁱ为发震时刻误差；为射线路径长度向量；为p波或sv波的群慢度，即为或矢量；分别为垂向速度的修正向量、ε的修正向量和δ的修正向量。

进一步地，在水平层状介质中p波和sv波的精确相速度公式为：

式中p表示纵波，sv表示质点振动的垂直平面内的横波分量，α0为纵波垂向速度，β0为横波垂向速度。

进一步地，在水平层状介质中p波和sv波的群速度计算公式为：

θg＝θ+δθ(5)

式中，θg为群角。v为p波和sv波的相速度。

进一步地，通过给定井中检波器获得到达时利用震源位置、发震时刻和初始模型来计算理论到达时两者旅行时残差为：

进一步地，通过lu分解、svd分解、共轭梯度算法或lsqr方法求解超定方程。

本发明方法只利用井中微地震自身的信息提取水平层状介质各向异性参数，特别适合没有三维地震和正交偶极子声波测井的地区。

附图说明

通过结合附图对本公开示例性实施方式进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开示例性实施方式中，相同的参考标号通常代表相同部件。

图1显示了相速度和群速度之间的关系。

图2显示了井中微地震观测的示例图。

图3显示了根据本发明实施例的利用井中微地震信息提取水平层状介质各向异性参数的方法流程图。

附图仅用于示例性说明，不能理解为对本专利的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的；相同或相似的标号对应相同或相似的部件；附图中描述位置关系的用语仅用于示例性说明，不能理解为对本专利的限制。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施方式。虽然附图中显示了本公开的优选实施方式，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

下面参照图1和图2描述本发明方法的原理。首先，利用thomsen参数表示相速度和群速度。水平层状介质各向异性参数一般使用传统的thomsen参数化方式。p波和sv波的精确相速度公式为：

式中p表示纵波(正号)，sv表示质点振动的垂直平面内的横波分量(负号)。α0为纵波垂向速度，β0为横波垂向速度。θ为入射角，p波和sv波都用一个符号表示，下标表明了何种波。

在水平层状介质中p波和sv波的群速度计算公式为：

θg＝θ+δθ(5)

式中，θg为群角。v为p波和sv波的相速度。

接下来，求取各向异性参数。

在井中观测微地震时，已知的信息包括微地震事件的初至信息、射孔或导爆索的具体位置和记录到的初至信息。利用这些信息可以计算出thomsen参数(是一种各向异性参数表示法，指利用ε,δ参数描述各项异性)。

微地震事件i(i＝1,…,n)到接收点k(k＝1,…,m)的旅行时可以用射线理论表示为：

其中，u为群慢度场。对于微地震监测来讲，实际观测的旅行时为其中τⁱ为发震时刻。

给定井中检波器得到的到达时(p波、sv波)，利用震源(射孔、导爆索、微地震事件)位置，发震时刻和初始模型来计算理论到达时初始模型中的信息有p波和sv波的垂向速度(一般可利用声波测井建立)和ε,δ，如果没有地面地震信息或vsp信息，初始模型信息一般为零)。

两者旅行时的残差为

方程(8)为非线性方程，利用泰勒展开来线性化该方程，则残差可以用thomsen参数的扰动表示为：

式中δτⁱ为发震时刻误差；为射线路径长度向量，可根据射线追踪获得；为p波或sv波的群慢度(即)；为或矢量；由方程(3)求出；分别为垂向速度的修正向量、ε的修正向量和δ的修正向量。向量元素的个数等于层的个数。

求解方程(9)组成的方程即可获得所需的各向异性参数(实际就是要求解修正向量)，计算方法可选lu分解、svd分解、共轭梯度算法、lsqr方法。

下面以lsqr方法为例，描述求解过程。

利用方程(9)可构建稀疏线性方程组ax＝b。

lsqr方法具有计算量小的优点，并且能很容易地利用矩阵的稀疏性简化计算，因此适合求解大型稀疏矩阵。

lsqr原理如下：

(1)初始化

(2)双对角化：

βi+1ui+1＝avi-αiui

αi+1vi+1＝a^tui+1-βi+1vi

(3)正交变换：

(4)解矢量：

wi+1＝vi+1-wiθi+1/ρi

上述迭代公式中，αi,βi为对角化矩阵a的标量序列，ui,vi为矢量序列，λ为阻尼因子，其余为中间变量。利用lsqr方法求解方程(9)是本领域技术人员容易实现的，再次不再赘述。

为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。

可以利用声波测井建立水平层状模型，模型信息可包括p波和sv波的垂向速度，各层ε,δ初始化为零；接下来利用观测井三分量检波器记录的数据，拾取微地震事件初至，例如纵波初至和横波初至；然后构建超定方程，利用数值计算方法求解即可得到各向异性参数。

参照图3，描述根据本发明一个实施例的方法流程。首先，输入数据，所述数据可来自观测井三分量检波器记录的数据，例如三分量数据的旅行时；接下来通过射线追踪获得射线路径长度向量；构造超定方程(9)，例如方程(9)里面包括了射线路径的长度，以及各层中各向异性参数与群慢度的关系，换言之，方程(9)就是旅行时残差与各层参数之间的关系；通过lsqr方法求解方程(9)获得各向异性参数，即求解修正向量等；接下来判断是否满足模型的修正量充分小，该过程可根据本领域技术人员的经验判断；不满足条件则更新模型，例如在初始模型各个参数加上这些修正量(加法运算)，然后再重新构造方程和求解；满足条件则输出模型，即输出每层中p波和sv波的垂向速度和ε,δ。

下面结合具体实施方式对本发明作进一步的说明。其中，附图仅用于示例性说明，表示的仅是示意图，而非实物图，不能理解为对本专利的限制；为了更好地说明本发明的实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。