基于微地震数据和SPSA优化算法的震源反演方法

摘要

本发明涉及了一种基于微地震数据和SPSA优化算法的震源反演方法，该方法包括以下步骤：根据测井资料及地下岩石性质资料分析建立水平层状速度模型；对微地震事件进行正演模拟，以二分法调整射线角度优选射线追踪路径，计算得到初至走时，然后计算相邻两检波器的走时之差；将走时之差加入随机扰动作为观测值；应用SPSA算法将迭代更新得到的微地震震源位置依照前述方法计算得到初至走时之差，通过拟合计算数据与观测值不断优化震源位置，检查计算数据与观测值之差，如果满足精度或者达到最大迭代次数则停止更新，输出结果，完成定位，否则继续进行下一次迭代。本发明通过快速高效的求解微地震事件源，提高了震源定位的准确性。

权利要求

1.一种基于微地震数据和SPSA优化算法的震源反演方法，包括以下步骤：步骤一：根据测井资料及地下岩石性质分析建立水平层状速度模型；

步骤二：对微地震事件进行正演模拟，计算得到初至走时，然后计算相邻两检波器的走时之差；

步骤三：将步骤二计算得到的走时之差加入在(0,1)内的一组服从高斯分布的随机数组成的随机扰动，以此作为观测值；

步骤四：应用SPSA算法，通过拟合计算数据与观测值，优化求解最优震源位置；

其特征在于：

步骤一包括以下步骤：

①根据测井资料和地下岩石性质分析对地层进行分类划分层位，地层共分为M个层位，由深到浅依次编号，介质分界面都是水平面，给定每一层距离地面深度；

②假定油层和上覆岩石均质，且孔隙介质中充满流体，则根据Raymer模型可以建立地层纵波(P波)速度场；

步骤二包括以下步骤：

(1)、基于Snell定律的射线追踪定位方法：假定试验区域内地层共分为M层，根据Snell定律可以计算出地震波从S点出发到达第i级检波器该射线在地层内传播角度、射线传播路径以及初至走时；

(2)、改进打靶法进行射线追踪时，首先给定射线初始位置的发射角度θ：地震波在层与层之间的传播遵循Snell定律，判断射线的终止位置是否为监测点的位置，如果射线终止位置不是监测点位置，则根据二分法不断调整射线角度，重新计算直至满足要求；

(3)、由于压裂过程持续时间比较长，震源发生时间不容易确定，因此，通过计算检波器之间的初至走时之差规避震源的发震时刻，而且该方法还能有效降低地层中的干扰程度，提高反演的稳健性；

步骤三包括以下步骤：

在(0,1)内生成一组服从高斯分布的随机数组成的随机扰动，把计算得到的初至走时之差加上随机扰动即可得到观测值；

步骤四包括以下步骤：

(1)、首先，设置初始参数，按照步骤二所述计算得到初至走时，进而得到初至走时之差，然后计算得到观测数据与计算数据之差；

(2)、更新迭代参数，同样按照步骤二所述计算得到初至走时、初至走时之差以及计算数据与观测数据之差，选定搜索迭代方向；

(3)、根据选取的迭代方向进行计算，得到更新后的微地震事件位置；

(4)、按照步骤二所述计算得到更新后的初至走时，进而得到更新后的初至走时之差以及计算数据与观测数据之差；

(5)、判断计算数据与观测数据之差是否满足精度，如果满足，停止更新；反之，选定新的搜索方向继续迭代计算。

说明书

基于微地震数据和SPSA优化算法的震源反演方法

技术领域

[0001] 本发明属于石油天然气地震勘探领域，具体地，涉及一种基于微地震数据和SPSA优化算法的震源反演方法。

背景技术

[0002] 在低渗透油气田生产过程中，压裂工艺是增产增注的一项重要措施，而压裂所产生裂缝(裂缝方位又与地应力有着密切的关系)和压裂规模，是井网部署的重要参考依据，因而深入地研究裂缝方位及形态，适时地进行井网调整，这一直是油田亟待解决的课题。

[0003] 压裂井人工微地震实时监测评价技术是建立在微地震监测技术基础上的一项油田生产动态监测技术。微地震监测是目前储层压裂中最精确、最及时、信息最丰富的监测手段。实时微地震成像可以及时指导压裂工程，适时调整压裂参数；对压裂的范围、裂缝发育的方向、大小进行追踪、定位，客观评价压裂工程的效果，对下一步的生产开发提供有效的指导，降低开发成本，更好的为后续油田管理提供依据。

[0004] 微地震监测的数据处理最核心的问题是震源定位，最常用的是基于走时的反演方法，如纵横波时差法或同型波时差法，这类方法是微地震震源反演中的常规方法，但是该方法需要求解线性或非线性方程组，在方程组的求解过程中可能存在超定、欠定或病态问题，此外该方法严重依赖检波器上所检测到的事件时间，定位精度不能满足要求。

[0005] 目前常用的反演方法也多种多样，包括一些梯度类算法，如牛顿法、共轭梯度法等，该类方法在优点是收敛较快，但是容易陷入局部最优解，而且对初值的选取依赖性大；还有一些随机类算法，如粒子群算法、遗传算法、模拟退火法等，该类算法可以很好地解决局部最优解的问题，但是该类算法求解的速度较慢。

发明内容

[0006] 为克服现有技术存在的缺陷，本发明提供一种基于微地震数据和SPSA优化算法的震源反演方法，能够快速高效实现震源定位。

[0007] 为实现上述目的，本发明采用下述方案：

[0008] 基于SPSA算法的微地震震源定位方法，包括以下步骤：

[0009] 步骤一：根据测井资料及地下岩石性质分析建立水平层状速度模型；

[0010] 步骤二：对微地震事件进行正演模拟，计算得到初至走时，然后计算相邻两检波器的走时之差；

[0011] 步骤三：将步骤二中计算得到的走时之差加入在(0，1)内的一组服从高斯分布的随机数组成的随机扰动，以此作为观测值；

[0012] 步骤四：应用SPSA算法，通过拟合计算数据与观测值，优化求解最优震源位置。

[0013] 在射线追踪定位方面，针对局部射线追踪法严重依赖于射线密度因素导致准确性差的缺点，本发明提出使用二分法改进了基于Snell定律的射线追踪的打靶法，首先是给定初始射线角度，然后根据二分法不断调整，经检验该方法不仅提高了收敛速度，而且精度也大幅度提升，提高了定位的速度和精确性。

[0014] 在目标函数的选取上，本发明使用同一震源在相邻检波器监测到的走时之差作为拟合参数，这样不仅规避了求解震源发生时刻的繁琐，而且降低了地震波在地层传播过程中收到的干扰程度，提高了抗噪能力。

[0015] 此外，本发明还提出了使用SPSA算法来反演震源位置，该方法把梯度类算法和随机类算法结合在一起考虑，利用随机类算法给定搜索方向及随机梯度，然后利用梯度优化求解，既能考虑全局最优，又能实现目标函数的快速求解，收敛速度较快，能够实现快速准确的震源定位，具有计算速度快、效率高的优点。

附图说明

[0016] 图1为基于微地震数据和SPSA优化算法的震源反演方法流程示意图；

[0017] 图2为基于Snell定律的射线追踪定位示意图；

[0018] 图3为SPSA算法优化反演流程图；

[0019] 图4为正演模拟得到的射线追踪路径图；

[0020] 图5为正演模拟得到的初至走时图示；

[0021] 图6为正演模拟得到的初至走时之差图示；

[0022] 图7为加入随机扰动后得到观测数据图示；

[0023] 图8为真实震源位置与反演得到震源位置图示。

具体实施方式

[0024] 如图1所示，基于微地震数据和SPSA优化算法的震源反演方法，包括如下步骤：

[0025] 步骤一：根据测井资料及地下岩石性质资料分析建立水平层状速度模型；

[0026] 具体方法如下：

[0027] ①根据测井资料和地下岩石性质分析对地层进行分类划分层位，地层共分为M个层位，由深到浅依次编号，介质分界面都是水平面，每一层距离地面深度为Layer_h＝[Layer_h1Layer_h2……Layer_hM]；

[0028] ②假定油层和上覆岩石均质，孔隙介质中充满流体，则根据Raymer模型可以建立地层速度场：

[0029] v＝(1-φ)2·vsolid+φ·vfluid

[0030] 其中，

[0031]

[0032] Kfliud＝Kw×Sw+Ko×So+Kg×Sg；

[0033]

[0034]

[0035]

[0036] ρfliud＝ρw×Sw+ρo×So+ρg×Sg；

[0037] 式中，v为速度(m/s)，φ为孔隙度，K为弹性模量(Pa)，μ为剪切模量(Pa)，ρ为密度(kg/cm3)，fi为固体中第i种组成所占的体积分数，Nsolid为固体中所包含的组分数，Sm为相m的饱和度；下标solid、fliud分别表示固体、流体，w、o、g分别表示水、油、气三相。

[0038] 由上式可以看出，Raymer模型综合考虑了地下固体和流体的基本特征，通过测井资料和地下岩石性质分析即可得到以上参数值，这样就可以计算得到地震波在每一层传播的速度。由于P波(纵波)具有传播速度快、易识别的优点，所以通常选择层界面的P波来进行震源定位，这样速度场从浅到深分别为Layer_v＝[Layer_v1Layer_v2……Layer_vM]。

[0039] 步骤二：对微地震事件进行正演模拟，计算得到初至走时，然后计算相邻两检波器的走时之差；

[0040] 具体方法如下：

[0041] ①基于Snell定律的射线追踪定义如下：

[0042] 试验区域内地层共分为M层，由浅到深每层中P波(纵波)的传播速度为Layer_v＝[Layer_v1Layer_v2……Layer_vM]，检波器级数为N，地震波从S点出发到达第i级检波器该射线在地层内传播角度、传播路径以及走时的计算公式如下：

[0043]

[0044]

[0045]

[0046] 其中，si是地震波从S点出发到达第i级检波器该射线在地层内的路径，ni是检波器编号(自上而下在所对应地层中的编号)，hij是射线在纵向投影高度，θij是地震波传播到第i级检波器过程中经过第j层地层时的角度(与垂直方向的夹角)，Layer_vj是地震波在第j层地层中传播的速度，Toibs是地震波从S点出发到达第i级检波器该射线在地层内的走时。

[0047] ②改进打靶法进行射线追踪时，首先给定射线初始位置的发射角度θ：

[0048] a＝0,b＝90

[0049] θ＝(a+b)/2(θ∈[a,b])

[0050] 其中，θ为射线初始位置的发射角度，a、b为初始位置发射角度的上下限。

[0051] 从震源点出发，地震波在层与层之间的传播遵循Snell定律，判断射线的终止位置是否为接收点的位置，如果射线终止位置不是接收点位置，若高于接收点位置，则有：

[0052] b＝(a+b)/2

[0053] θ＝(a+b)/2(θ∈[a,b])

[0054] 反之，则有：

[0055] a＝(a+b)/2

[0056] θ＝(a+b)/2(θ∈[a,b])

[0057] 根据二分法不断调整射线角度，重新计算直至满足要求。

[0058] ③由于压裂过程持续时间比较长，震源发生时间不容易确定，因此通过计算检波器之间的初至走时之差规避震源的发震时刻，而且还能有效降低地层中的干扰程度，计算如下：

[0059]

[0060] 其中， 为相邻两级检波器走时之差。

[0061] 步骤三：将步骤二中走时之差加入在(0,1)内的一组服从高斯分布的随机数组成的随机扰动作为观测值；

[0062] 具体方法如下：

[0063] 在(0,1)内生成一组服从高斯分布的随机数组成的随机扰动

[0064] Δt＝[Δt1Δt2……ΔtN-1]

[0065] 则观测值为

[0066]

[0067]

[0068] 其中，Δt为随机扰动向量，ΔTobs为观测数据向量。

[0069] 步骤四：应用SPSA算法，通过拟合计算数据与观测值，优化求解最优震源位置；

[0070] 具体方法如下：

[0071] ①随机给定微地震事件初始位置X0＝(x0,y0,z0)，按照步骤(2)所述计算得到初至走时 进而得到初至走时之差

[0072]

[0073] 得到计算数据与观测数据之差

[0074]

[0075] 给定k＝0(k为当前搜索次数)，maxIter＝100(最大搜索次数)。

[0076] 其中，ΔTcal0为初始位置的初至走时之差向量，||ΔT0||为此时的计算数据与观测数据之差。

[0077] ②搜索次数k加1，随机给定一个迭代方向α(在[-1,1]之间且不为0)及方向步长ΔX，则有

[0078] k＝k+1

[0079] X1＝X0+αΔX

[0080] 同样按照按照步骤二所述计算得到初至走时 进而得到初至走时之差

[0081]

[0082] 同样可得计算数据与观测数据之差

[0083]

[0084] 如果||ΔT0||＞||ΔT1||，则选取方向为 否则

[0085] 其中，ΔTcal1为计算更新得到的初至走时之差向量，||ΔT1||为此时的计算数据与观测数据之差，α为迭代方向，ΔX为迭代方向的步长。

[0086] ③计算更新微地震事件位置

[0087] Xupdate＝X0+αgStep

[0088] 其中，Xupdate为更新后的震源位置，Step为迭代步长。

[0089] ④按照按照步骤二所述计算得到初至走时 进而得到初至走时之差

[0090]

[0091] 得到计算数据与观测数据之差

[0092]

[0093] 其中，g(Xupdate)为迭代更新的目标函数，ΔTupdate为更新的初至走时之差向量，ΔTobs为观测数据向量。

[0094] ⑤判断g(Xupdate)是否满足精度，如果满足则输出Xbest＝Xupdate；否则X0＝Xupdate，返回继续执行③迭代更新微地震事件位置，直至g(Xupdate)不再下降，此时返回执行②，判断k是否大于maxIter，若是则输出Xbest＝Xupdate，否则重新选取迭代方向α，搜索次数k加1。此外，在迭代更新过程中如果出现g(Xupdate)满足精度要求则输出Xbest＝Xupdate，停止更新。

[0095] 实施例

[0096] 基于步骤一，根据测井资料及地下岩石性质分析建立水平层状速度模型：本次算例选长400米、宽400米、高400米的区域，地层共分为4个层位，由深到浅依次编号，介质分界面都是水平面，每一层距离地面深度为Layer_h＝[2100 2150 2300 2400]m，层界面的P波(纵波)速度从浅到深分别为Layer_v＝[1600 2000 2400 2800]m/s。

[0097] 基于步骤二，对给定微地震事件进行正演模拟，计算得到初至走时，然后计算相邻两检波器的走时之差：

[0098] 设置检波器参数如下表1所示：

[0099] 表1检波器参数

[0100]级数 1 2 3 4 5 6 7 8 9 10 11 12
X 0 0 0 0 0 0 0 0 0 0 0 0
Y 0 0 0 0 0 0 600 600 600 600 600 600
Z 2017 2032 2047 2062 2077 2092 2107 2122 2137 2152 2167 2182
级数 13 14 15 16 17 18 19 20 21 22 23 24
X 600 600 600 600 600 600 600 600 600 600 600 600
Y 0 0 0 0 0 0 600 600 600 600 600 600
Z 2197 2212 2227 2242 2257 2272 2287 2302 2317 2332 2347 2362

[0101] Z坐标为相对度地面的深度，第一级检波器深度是-2017米，第24级检波器深度是-2362米，每级检波器高差15米。

[0102] 设置真实震源事件位置如下表2所示：

[0103] 表2真实震源事件参数

[0104]

[0105] 根据改进的射线追踪打靶法，并基于二分法调整射线角度，正演模拟得到射线追踪路径见图4。

[0106] 在得到路径后，除以所在层的地震波传播速度，得到震源到检波器位置处的初至走时分布见图5，进而得到初至走时之差分布见图6。

[0107] 基于步骤三，将步骤二计算得到的走时之差加入在(0,1)内的一组服从高斯分布的随机数组成的随机扰动，以此作为观测值；

[0108] 在(0,1)内生成一组服从高斯分布的随机数组成的随机扰动Δt＝[Δt1Δt2……ΔtN-1]，则观测值为

[0109]

[0110] 加入随机扰动后得到观测数据分布见图7。

[0111] 基于步骤四，应用SPSA算法，通过拟合计算数据与观测值，优化求解最优震源位置：

[0112] 利用SPSA优化方法反演震源位置流程图见图3，按照图3所示步骤执行，反演得到震源参数如表3所示：

[0113] 表3反演得到的震源位置

[0114]

[0115] 真实震源位置与反演得到震源位置见图8。

[0116] 由表3和图8可以直观看出，由于在真实值的基础上加入了一定的扰动之后，真实震源位置与反演得到的震源位置误差最大为21.9219m，最小为4.0864m，经程序测试反演总共用时约为140s，反演的结果在合理范围内。