CN116660980A - 一种基于改进程函方程的微地震定位方法 - Google Patents
一种基于改进程函方程的微地震定位方法 Download PDFInfo
- Publication number
- CN116660980A CN116660980A CN202310730999.4A CN202310730999A CN116660980A CN 116660980 A CN116660980 A CN 116660980A CN 202310730999 A CN202310730999 A CN 202310730999A CN 116660980 A CN116660980 A CN 116660980A
- Authority
- CN
- China
- Prior art keywords
- grid
- point
- grid point
- representing
- travel
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 48
- 238000012544 monitoring process Methods 0.000 claims abstract description 13
- 238000004364 calculation method Methods 0.000 claims description 16
- 230000001351 cycling effect Effects 0.000 claims description 3
- 238000001514 detection method Methods 0.000 claims description 3
- 230000004807 localization Effects 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 25
- 238000012360 testing method Methods 0.000 description 6
- 230000007547 defect Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000000977 initiatory effect Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/288—Event detection in seismic signals, e.g. microseismics
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Emergency Management (AREA)
- Business, Economics & Management (AREA)
- Acoustics & Sound (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明公开了一种基于改进程函方程的微地震定位方法，包括以下步骤：S1、将圈定的监测区域划分为二维离散网格，并在二维离散网格点上布设检波器；S2、依次将每个检波器作为新震源，采用球面波近似算法计算新震源附近网格点的走时；S3、根据步骤S2得到的新震源附近网格点的走时，采用多模板快速行进算法计算其余网格点的走时；S4、利用时间残差函数计算所有网格点走时的函数值，根据计算的函数值定位地震震源位置；本发明采用球面波近似算法计算新震源附近网格点的走时，同时采用多模板快速行进算法计算其余网格点的走时，降低了走时场在对角线方向以及新震源点附近的误差，提高了走时场的精确度，可以更加准确的定位微震震源位置。
Description
技术领域
本发明涉及微地震定位领域，具体涉及一种基于改进程函方程的微地震定位方法。
背景技术
对于非常规油气藏的勘探开发，微地震监测技术是定量分析水力压裂油藏改造过程中裂缝发育情况的有效办法。
微地震事件反演定位是微地震监测中的核心问题。根据不同的反演理论，大体上可以将微地震的反演问题分为走时反演和能量反演两种类型。最早出现的微地震反演定位方法借鉴了天然地震中的定位方法，如Geiger定位法、三角定位法等。此类方法属于走时反演类方法，其主要原理是根据微地震事件在地层中传播的空间、时间关系，将微地震信号从震源到检波器的理论走时与实际观测走时之差作为目标函数，而求解震源位置与发震时刻。随后出现了基于能量叠加及类偏移方法的微地震震源定位方法，如网格搜索法、逆时偏移法等，该类方法避免了微地震事件的初至拾取，寻找叠加后能量最大值点作为微地震震源。
传统走时反演类方法所使用的射线追踪方法在地质模型强非均匀时会存在焦散、阴影区等问题，导致射线路径难以确定。通过数值求解程函方程的方法可以准确获得非均匀介质中的射线路径。目前常用来解程函方程的方法是快速行进法，但该方法计算出的走时场在斜对角方向以及源点附近误差较大。
发明内容
针对现有技术中的上述不足，本发明提供了一种基于改进程函方程的微地震定位方法，解决传统射线追踪方法在地质强非均匀情况下存在阴影区的缺点，与快速行进算法解程函方程相比较，本发明提高了走时场的精确度，从而实现微地震信号的精确定位。
为了达到上述发明目的，本发明采用的技术方案为：
一种基于改进程函方程的微地震定位方法，包括以下步骤：
S1、将圈定的监测区域划分为二维离散网格，并在二维离散网格点上布设检波器；
S2、依次将每个检波器作为新震源，采用球面波近似算法计算新震源附近网格点的走时；
S3、根据步骤S2得到的新震源附近网格点的走时，采用多模板快速行进算法计算其余网格点的走时；
S4、利用时间残差函数计算所有网格点走时的函数值，根据计算的函数值定位地震震源位置。
进一步的，步骤S1具体包括：
根据微地震监测任务，圈定待监测的区域，建立二维直角坐标系，以微地震事件平面投影到原点的距离作为横坐标，以微地震事件的深度作为纵坐标，将圈定的监测区域划分为二维离散网格，每个网格的大小为d，并在二维离散网格点上以竖井方式布设检波器。
进一步的，步骤S2具体包括：
S21、将检波器(xk,zk)作为新震源；
S22、根据新震源到达网格点t1(x,z1)和网格点t2(x,z2)的走时，在两个网格点之间插入一个网格点t0(x,z0)，采用Lagrange插值法计算新震源到达网格点t0(x,z0)的走时，并计算新震源经过网格点t0(x,z0)到达网格点t(x+Δx,z2)的走时；
S23、最小化网格点t(x+Δx,z2)的走时，计算网格点t(x+Δx,z2)的走时对新震源到达网格点t0(x,z0)的深度z0的偏导逼近于0的点，得到新震源到达网格点t0(x,z0)的深度z0，进而求出当网格点t(x+Δx,z2)的走时最小时，网格点t0(x,z0)和网格点t(x+Δx,z2)的走时大小。
进一步的，步骤S22中新震源到网格点t0(x,z0)的走时的计算公式为：
其中，t0表示网格点t0(x,z0)的走时，t1、t2表示新震源分别到达这两个网格点的走时，z0表示从新震源到达网格点t0(x,z0)的深度，z1表示从新震源到达网格点t1(x,z1)的深度，z2表示从新震源到达网格点t2(x,z2)的深度。
进一步的，步骤S22中新震源到达网格点t(x+Δx,z2)的走时的计算公式为：
其中，t表示新震源附近网格点的走时，S表示单元格慢度，t0表示网格点t0(x,z0)的走时，z0表示从新震源到达网格点t0(x,z0)的深度，z2表示从新震源到达网格点t2(x,z2)的深度，Δx表示x坐标轴方向上的离散网格间距。
进一步的，步骤S23中最小化点t(x+Δx,z2)的走时的计算公式为：
其中，表示求网格点t(x+Δx,z2)的走时对新震源到达网格点t0(x,z0)的深度z0的偏导，W表示平均慢度，S表示单元格慢度，Δx表示x坐标轴方向上的离散网格间距，t0表示网格点t0(x,z0)的走时，z0表示从新震源到达网格点t0(x,z0)的深度，z1表示从新震源到达网格点t1(x,z1)的深度，z2表示从新震源到达点t2(x,z2)的深度。
进一步的，步骤S3具体包括：
S31、将步骤S2计算得到的新震源附近网格点标记为近点，将其余网格点标记为远点；
S32、多模板快速行进算法分别用模板1和模板2进行计算，将模板1用快速行进算法计算二维直角坐标系横轴或纵轴上，与近点距离为一个网格大小的网格点的走时，将模板2用快速行进算法计算二维直角坐标系对角线上，与近点距离为网格大小的网格点的走时，根据一阶差分公式和二阶差分公式进行计算，将计算得到的最小值作为网格点的走时，并将这些点标记为窄带点；
S33、在所有窄带点中寻找具有最小走时的网格点，将其标记为近点；
S34、寻找步骤S33中与近点相邻的远点标记为窄带点；
S35、用两个模板得出的一阶差分近似和二阶差分近似公式更新步骤S34中找到的与新增近点相邻的窄带点走时；
S36、循环步骤S33至步骤S35，直到窄带点为空。
进一步的，模板1的一阶差分近似公式如下：
其中，i表示沿x坐标轴方向上的离散网格点编号，j表示沿y坐标轴方向上的离散网格点编号，Δx表示x坐标轴方向上的离散网格间距，Δy表示y坐标轴方向上的离散网格间距，Ti,j表示网格点(i,j)的到时，T1、T2分别表示在网格点(i,j)处x坐标轴方向和y坐标轴方向上向前向后的最小时间，T1＝min(Ti-1,j,Ti+1,j)，T2＝min(Ti,j-1,Ti,j+1)，Si,j表示网格点(i,j)的慢度。
进一步的，模板1的二阶差分近似公式如下：
其中，Ti,j表示网格点(i,j)的到时，Δx表示表示x坐标轴方向上的离散网格间距，Δy表示y坐标轴方向上的离散网格间距，Si,j表示网格点(i,j)的慢度，T1、T2分别表示在网格点(i,j)处x坐标轴方向和y坐标轴方向上向前向后的最小时间，
进一步的，模板2的一阶差分近似公式如下：
其中，Ti,j表示网格点(i,j)的到时，Tv表示在网格点(i,j)处两条对角线方向上向前向后的最小时间，Si,j表示网格点(i,j)的慢度，d表示网格大小。
进一步的，模板2的二阶差分近似公式如下：
其中，Ti,j表示网格点(i,j)的到时，Tv表示在网格点(i,j)处两条对角线方向上向前向后的最小时间，Si,j表示网格点(i,j)的慢度，d表示网格大小。
进一步的，步骤S4具体包括：
寻找时间残差函数计算得到的最小函数值，找到的最小函数值所对应的网格点即为地震震源位置。
进一步的，步骤S4中时间残差函数具体为：
其中，fNs(x,z)表示点(x,z)的时间残差函数值，k表示检波器数量，Ns表示地震事件数量，Ti(x,z)表示第i个检波器计算得到的走时场中点(x,z)的时间，ti(j)表示第i个检波器检测到第j个微地震事件的到时，τ0(j)表示第j个微地震事件的发震时间。
本发明具有以下有益效果：
本发明采用球面波近似算法计算新震源附近网格点的走时，同时采用多模板快速行进算法计算其余网格点走时，降低了走时场在对角线方向以及新震源点附近的误差，提高了走时场的精确度，可以更加准确的定位微震震源位置。
附图说明
图1为本发明一种基于改进程函方程的微地震定位方法的流程示意图；
图2为球面波近似算法的示意图；
图3为多模板快速行进算法所用的模板1和模板2的示意图；
图4为测试阶段在T1模型下，震源位置为(51,51)的走时场误差示意图；
图5为测试阶段在T1模型下，震源位置为(1,1)的走时场误差示意图；
图6为测试阶段在T2模型下，震源位置为(51,51)的走时场误差示意图；
图7为测试阶段在T3模型下，震源位置为(51,51)的走时场误差示意图；
图8为测试阶段层状速度模型示意图；
图9为测试阶段的二维层状观测系统和定位效果对比示意图。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
如图1所示，一种基于改进程函方程的微地震定位方法，包括以下步骤S1-S4：
S1、将圈定的监测区域划分为二维离散网格，并在二维离散网格点上布设检波器。
在本发明的一个可选实施例中，本实施例中采用基于网格扩展的射线传播计算方法计算射线初至到时。
具体的，步骤S1包括：根据微地震监测任务，圈定待监测的区域，建立二维直角坐标系，以微地震事件平面投影到原点的距离作为横坐标，以微地震事件的深度作为纵坐标，将圈定的监测区域划分为二维离散网格，每个网格的大小为d，并在二维离散网格点上以竖井方式布设检波器，且布设检波器的个数为k。
S2、依次将每个检波器作为新震源，采用球面波近似算法计算新震源附近网格点的走时。
在本发明的一个可选实施例中，为解决快速行进算法计算走时在震源点附近误差较大的缺点，本实施例中根据波传播的互易原理，在同一模型中，直达波从震源传播至检波器的走时与该波从检波器传播至震源的走时是相等的，将检波器(xk,zk)作为新震源，采用球面波近似算法计算新震源点附近(xk-h:xk+h:zk-h:zk+h)网格点的走时t。且将地震波以球面波近似形式计算，更加符合地震波传播特性，该方法适用于任何复杂的各向同性介质，并且在震源点附近精确度高，但是计算效率低，所以本实施例仅采用该方法计算新震源点附近部分网格点的走时。
如图2所示，步骤S2包括以下步骤S21-S23：
S21、将检波器(xk,zk)作为新震源。
S22、根据新震源到达网格点t1(x,z1)和网格点t2(x,z2)的走时，在两个网格点之间插入一个网格点t0(x,z0)，采用Lagrange插值法计算新震源到达网格点t0(x,z0)的走时，并计算新震源经过网格点t0(x,z0)到达网格点t(x+Δx,z2)的走时；
在本发明的一个可选实施例中，本实施例中球面波近似算法提出一个网格点的走时，可以通过周围已知的网格点走时计算出来，所以本实施例用新震源到达网格点t1(x,z1)和网格点t2(x,z2)的走时计算新震源到达网格点t0(x,z0)的走时和新震源经过网格点t0(x,z0)到达网格点t(x+Δx,z2)的走时。
具体的，步骤S22中新震源到网格点t0(x,z0)的走时的计算公式为：
其中，t0表示网格点t0(x,z0)的走时，t1、t2表示新震源分别到达这两个网格点的走时，z0表示从新震源到达网格点t0(x,z0)的深度，z1表示从新震源到达网格点t1(x,z1)的深度，z2表示从新震源到达网格点t2(x,z2)的深度。
具体的，步骤S22中新震源到达网格点t(x+Δx,z2)的走时的计算公式为：
其中，t表示新震源附近网格点的走时，S表示单元格慢度，t0表示网格点t0(x,z0)的走时，z0表示从新震源到达网格点t0(x,z0)的深度，z2表示从新震源到达网格点t2(x,z2)的深度，Δx表示x坐标轴方向上的离散网格间距。
S23、最小化网格点t(x+Δx,z2)的走时，计算网格点t(x+Δx,z2)的走时对新震源到达网格点t0(x,z0)的深度z0的偏导逼近于0的点，得到新震源到达网格点t0(x,z0)的深度z0，进而求出当网格点t(x+Δx,z2)的走时最小时，网格点t0(x,z0)和网格点t(x+Δx,z2)的走时大小。
具体的，步骤S23中最小化点t(x+Δx,z2)的走时的计算公式为：
其中，表示求网格点t(x+Δx,z2)的走时对新震源到达网格点t0(x,z0)的深度z0的偏导，W表示平均慢度，仅与当前周围坐标和走时有关，脱离于虚拟震源的联系，S表示单元格慢度，Δx表示x坐标轴方向上的离散网格间距，t0表示网格点t0(x,z0)的走时，z0表示从新震源到达网格点t0(x,z0)的深度，z1表示从新震源到达网格点t1(x,z1)的深度，z2表示从新震源到达点t2(x,z2)的深度。
S3、根据步骤S2得到的新震源附近网格点的走时，采用多模板快速行进算法计算其余网格点的走时。
在本发明的一个可选实施例中，本实施例将步骤S2得到的新震源附近网格点的走时t作为初始区域，带入多模板快速行进算法计算其余网格点的走时，得到整个监测区域的网格点走时，而传统的多模板快速行进算法从震源一点开始计算其他点的走时，这样做会导致震源附近误差较大，而本实施例中提出的多模板快速行进算法是在传统的快速行进算法的基础上，结合了多模板和方向导数计算其余网格点的走时，这样计算既可以弥补多模板快速行进算法在震源附近精度不高的缺点，又可以提升计算效率。
如图3所示，步骤S3包括以下步骤S31-S36：
S31、将步骤S2计算得到的新震源附近网格点标记为近点，将其余网格点标记为远点。
S32、多模板快速行进算法分别用模板1和模板2进行计算，将模板1用快速行进算法计算二维直角坐标系横轴或纵轴上，与近点距离为一个网格大小的网格点的走时，将模板2用快速行进算法计算二维直角坐标系对角线上，与近点距离为网格大小的网格点的走时，根据一阶差分公式和二阶差分公式进行计算，将计算得到的最小值作为网格点的走时，并将这些点标记为窄带点；
在本发明的一个可选实施例中，由于快速行进算法和二阶快速扫描算法都忽略了对角线上网格点的信息，因此在对角线方向上容易产生较大的数值误差，所以，本实施例用模板1求覆盖的最近的邻点，用模板2求覆盖的对角线方向上的邻点。
具体的，模板1的一阶差分近似公式如下：
其中，i表示沿x坐标轴方向上的离散网格点编号，j表示沿y坐标轴方向上的离散网格点编号，Δx表示x坐标轴方向上的离散网格间距，Δy表示y坐标轴方向上的离散网格间距，Ti,j表示网格点(i,j)的到时，T1、T2分别表示在网格点(i,j)处x坐标轴方向和y坐标轴方向上向前向后的最小时间，T1＝min(Ti-1,j,Ti+1,j)，T2＝min(Ti,j-1,Ti,j+1)，Si,j表示网格点(i,j)的慢度。
具体的，模板1的二阶差分近似公式如下：
其中，Ti,j表示网格点(i,j)的到时，Δx表示表示x坐标轴方向上的离散网格间距，Δy表示y坐标轴方向上的离散网格间距，Si,j表示网格点(i,j)的慢度，T1、T2分别表示在网格点(i,j)处x坐标轴方向和y坐标轴方向上向前向后的最小时间，
具体的，模板2的一阶差分近似公式如下：
其中，Ti,j表示网格点(i,j)的到时，Tv表示在网格点(i,j)处两条对角线方向上向前向后的最小时间，Si,j表示网格点(i,j)的慢度，d表示网格大小。
具体的，模板2的二阶差分近似公式如下：
其中，Ti,j表示网格点(i,j)的到时，Tv表示在网格点(i,j)处两条对角线方向上向前向后的最小时间，Si,j表示网格点(i,j)的慢度，d表示网格大小。
S33、在所有窄带点中寻找具有最小走时的网格点，将其标记为近点。
S34、寻找步骤S33中与近点相邻的远点标记为窄带点。
S35、用两个模板得出的一阶差分近似和二阶差分近似公式更新步骤S34中找到的与新增近点相邻的窄带点走时。
S36、循环步骤S33至步骤S35，直到窄带点为空。
S4、利用时间残差函数计算所有网格点走时的函数值，根据计算的函数值定位地震震源位置。
在本发明的一个可选实施例中，本实施例中根据所有的检波器计算得到的走时场中网格点(x,z)的时间，并根据检波器检测到的事件所对应的到时和发震时间，利用时间残差函数计算这些网格点走时的函数值，寻找时间残差函数计算得到的最小函数值，找到的最小函数值所对应的网格点即为地震震源位置。
具体的，步骤S4中时间残差函数具体为：
其中，fNs(x,z)表示点(x,z)的时间残差函数值，k表示检波器数量，Ns表示地震事件数量，Ti(x,z)表示第i个检波器计算得到的走时场中点(x,z)的时间，ti(j)表示第i个检波器检测到第j个微地震事件的到时，τ0(j)表示第j个微地震事件的发震时间。
具体的，本实施例中采用三种速度模型来验证本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW的准确性，并与快速行进算法FMM和多模板快速行进算法MSFM进行比较，三种速度模型的真实走时场由表1给出的公式计算，其速度模型可由公式求得，表1如下所示：
表1三种速度模型真实走时场的计算公式
其中，x0和y0表示震源坐标，x和y表示网格区域内任意网格点，T1(x)表示在T1模型下的真实走时场，t2(x)表示在T2模型下的真实走时场，T3(x)表示在T3模型下的真实走时场。
本实施例为了测量计算得到的三种速度模型的走时场和真实走时场Ta(x)的解的误差，采用下述误差规范公式进行规范：
l∞＝mAx(|T-Ta|)
其中，n表示网格点的总数,T表示计算出的走时场，Ta表示真实走时场，L1、L2、L∞表示三种不同的误差规范。
本实施例中表2、表3通过对比本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW、快速行进算法FMM以及多模板快速行进算法MSFM这三种算法，在不同速度模型下震源位置不同时的误差，并结合图4、图5、图6、图7给出的在T1、T2、T3三种速度模型，震源在(51,51)和(1,1)位置时，上述三种算法的走时误差示意图，从中可以发现采用本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW计算得到的走时误差小于快速行进算法FMM和多模板快速行进算法MSFM的走时误差，而准确的走时是精确定位的前提，所以采用走时误差更小的算法可以使定位更加准确。
本实施例中表2如下所示：
表2从大小为101×101的网格的不同源点计算出的速度模型T1的误差规范
本实施例中表3如下所示：
表3从大小为101×101的网格中心计算出的速度模型T2、T3的误差规范
其中，图4表示在T1模型下，震源位置为(51,51)的走时场误差示意图，图4中(a)表示本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW的走时误差示意图，图4中(b)表示多模板快速行进算法MSFM的走时误差示意图，图4中(c)表示快速行进算法FMM的走时误差示意图；图5表示在T1模型下，震源位置为(1,1)的走时场误差示意图，图5中(a)表示本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW的走时误差示意图，图5中(b)表示多模板快速行进算法MSFM的走时误差示意图，图5中(c)表示快速行进算法FMM的走时误差示意图；图6表示在T2模型下，震源位置为(51,51)的走时场误差示意图，图6中(a)表示本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW的走时误差示意图，图6中(b)表示多模板快速行进算法MSFM的走时误差示意图，图6中(c)表示快速行进算法FMM的走时误差示意图；图7表示在T3模型下，震源位置为(51,51)的走时场误差示意图，图7中(a)表示本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW的走时误差示意图，图7中(b)表示多模板快速行进算法MSFM的走时误差示意图，图7中(c)表示快速行进算法FMM的走时误差示意图。
具体的，为了验证本发明层状模型的走时准确性，本实施例中建立如图8所示的层状速度模型，图8中倒三角形表示震源位置，黑色空心圆表示走时计算点，并设置震源位置为(400,321)，在观测系统中设置117个网格点检测走时，并将传统射线追踪方法计算的走时作为真实走时来做比较，得到三种方法的走时误差表4，表4如下所示：
表4从大小为800×600层状模型下117个网格点的走时误差
本实施例中通过分析表4可以得出本发明所用方法在层状模型下精确度依旧很高，能够满足定位要求，并且采用本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW得到的走时误差小于采用快速行进算法FMM和多模板快速行进算法MSFM得到的走时误差。
具体的，为验证本发明的定位准确性，本实施例中建立如图9(a)所示的二维层状观测系统，图9(a)中黑色倒三角形表示台站位置，黑色圆形表示地震事件真实位置，并利用12个台站对12个事件进行定位，使用传统射线追踪方法计算的走时作为真实走时，采用本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW以及快速行进算法FMM进行定位，定位效果如图9(b)所示，图9(b)表示两种算法的定位效果对比图，左三角形表示快速行进算法FMM定位位置，黑色叉号表示本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW定位位置，从中可以得到采用本发明结合的球面波近似算法和多模板快速行进算法MSFM_SW定位的误差会比采用快速行进算法FMM进行定位的误差更小，更接近真实位置。
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。
Claims (9)
1.一种基于改进程函方程的微地震定位方法，其特征在于，包括以下步骤：
S1、将圈定的监测区域划分为二维离散网格，并在二维离散网格点上布设检波器；
S2、依次将每个检波器作为新震源，采用球面波近似算法计算新震源附近网格点的走时；
S3、根据步骤S2得到的新震源附近网格点的走时，采用多模板快速行进算法计算其余网格点的走时；
S4、利用时间残差函数计算所有网格点走时的函数值，根据计算的函数值定位地震震源位置。
2.根据权利要求1所述的一种基于改进程函方程的微地震定位方法，其特征在于，步骤S1具体包括：
根据微地震监测任务，圈定待监测的区域，建立二维直角坐标系，以微地震事件平面投影到原点的距离作为横坐标，以微地震事件的深度作为纵坐标，将圈定的监测区域划分为二维离散网格，每个网格的大小为d，并在二维离散网格点上以竖井方式布设检波器。
3.根据权利要求1所述的一种基于改进程函方程的微地震定位方法，其特征在于，步骤S2具体包括：
S21、将检波器(xk,zk)作为新震源；
S22、根据新震源到达网格点t1(x,z1)和网格点t2(x,z2)的走时，在两个网格点之间插入一个网格点t0(x,z0)，采用Lagrange插值法计算新震源到达网格点t0(x,z0)的走时，并计算新震源经过网格点t0(x,z0)到达网格点t(x+Δx,z2)的走时；
S23、最小化网格点t(x+Δx,z2)的走时，计算网格点t(x+Δx,z2)的走时对新震源到达网格点t0(x,z0)的深度z0的偏导逼近于0的点，得到新震源到达网格点t0(x,z0)的深度z0，进而求出当网格点t(x+Δx,z2)的走时最小时，网格点t0(x,z0)和网格点t(x+Δx,z2)的走时大小。
4.根据权利要求3所述的一种基于改进程函方程的微地震定位方法，其特征在于，步骤S22中新震源到网格点t0(x,z0)的走时的计算公式为：
其中，t0表示网格点t0(x,z0)的走时，t1、t2表示新震源分别到达这两个网格点的走时，z0表示从新震源到达网格点t0(x,z0)的深度，z1表示从新震源到达网格点t1(x,z1)的深度，z2表示从新震源到达网格点t2(x,z2)的深度；
步骤S22中新震源到达网格点t(x+Δx,z2)的走时的计算公式为：
其中，t表示新震源附近网格点的走时，S表示单元格慢度，t0表示网格点t0(x,z0)的走时，z0表示从新震源到达网格点t0(x,z0)的深度，z2表示从新震源到达网格点t2(x,z2)的深度，Δx表示x坐标轴方向上的离散网格间距；
步骤S23中最小化点t(x+Δx,z2)的走时的计算公式为：
其中，表示求网格点t(x+Δx,z2)的走时对新震源到达网格点t0(x,z0)的深度z0的偏导，W表示平均慢度，S表示单元格慢度，Δx表示x坐标轴方向上的离散网格间距，t0表示网格点t0(x,z0)的走时，z0表示从新震源到达网格点t0(x,z0)的深度，z1表示从新震源到达网格点t1(x,z1)的深度，z2表示从新震源到达点t2(x,z2)的深度。
5.根据权利要求1所述的一种基于改进程函方程的微地震定位方法，其特征在于，步骤S3具体包括：
S31、将步骤S2计算得到的新震源附近网格点标记为近点，将其余网格点标记为远点；
S32、多模板快速行进算法分别用模板1和模板2进行计算，将模板1用快速行进算法计算二维直角坐标系横轴或纵轴上，与近点距离为一个网格大小的网格点的走时，将模板2用快速行进算法计算二维直角坐标系对角线上，与近点距离为网格大小的网格点的走时，根据一阶差分公式和二阶差分公式进行计算，将计算得到的最小值作为网格点的走时，并将这些点标记为窄带点；
S33、在所有窄带点中寻找具有最小走时的网格点，将其标记为近点；
S34、寻找步骤S33中与近点相邻的远点标记为窄带点；
S35、用两个模板得出的一阶差分近似和二阶差分近似公式更新步骤S34中找到的与新增近点相邻的窄带点走时；
S36、循环步骤S33至步骤S35，直到窄带点为空。
6.根据权利要求5所述的一种基于改进程函方程的微地震定位方法，其特征在于，模板1的一阶差分近似公式如下：
其中，i表示沿x坐标轴方向上的离散网格点编号，j表示沿y坐标轴方向上的离散网格点编号，Δx表示x坐标轴方向上的离散网格间距，Δy表示y坐标轴方向上的离散网格间距，Ti,j表示网格点(i,j)的到时，T1、T2分别表示在网格点(i,j)处x坐标轴方向和y坐标轴方向上向前向后的最小时间，T1＝min((Ti-1,j,Ti+1,j)，T2＝min(Ti,j-1,Ti,j+1)，Si,j表示网格点(i,j)的慢度；
模板1的二阶差分近似公式如下：
其中，Ti,j表示网格点(i,j)的到时，Δx表示表示x坐标轴方向上的离散网格间距，Δy表示y坐标轴方向上的离散网格间距，Si,j表示网格点(i,j)的慢度，T1、T2分别表示在网格点(i,j)处x坐标轴方向和y坐标轴方向上向前向后的最小时间，
7.根据权利要求5所述的一种基于改进程函方程的微地震定位方法，其特征在于，模板2的一阶差分近似公式如下：
其中，Ti,j表示网格点(i,j)的到时，Tv表示在网格点(i,j)处两条对角线方向上向前向后的最小时间，Si,j表示网格点(i,j)的慢度，d表示网格大小；
模板2的二阶差分近似公式如下：
其中，Ti,j表示网格点(i,j)的到时，Tv表示在网格点(i,j)处两条对角线方向上向前向后的最小时间，Si,j表示网格点(i,j)的慢度，d表示网格大小。
8.根据权利要求1所述的一种基于改进程函方程的微地震定位方法，其特征在于，步骤S4具体包括：
寻找时间残差函数计算得到的最小函数值，找到的最小函数值所对应的网格点即为地震震源位置。
9.根据权利要求1所述的一种基于改进程函方程的微地震定位方法，其特征在于，步骤S4中时间残差函数具体为：
其中，fNs(x,z)表示点(x,z)的时间残差函数值，k表示检波器数量，Ns表示地震事件数量，Ti(x,z)表示第i个检波器计算得到的走时场中点(x,z)的时间，ti(j)表示第i个检波器检测到第j个微地震事件的到时，τ0(j)表示第j个微地震事件的发震时间。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310730999.4A CN116660980A (zh) | 2023-06-19 | 2023-06-19 | 一种基于改进程函方程的微地震定位方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310730999.4A CN116660980A (zh) | 2023-06-19 | 2023-06-19 | 一种基于改进程函方程的微地震定位方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116660980A true CN116660980A (zh) | 2023-08-29 |
Family
ID=87727825
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310730999.4A Pending CN116660980A (zh) | 2023-06-19 | 2023-06-19 | 一种基于改进程函方程的微地震定位方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116660980A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117607957A (zh) * | 2024-01-24 | 2024-02-27 | 南方科技大学 | 基于等效慢度的快速推进法的地震波走时求解方法及系统 |
-
2023
- 2023-06-19 CN CN202310730999.4A patent/CN116660980A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117607957A (zh) * | 2024-01-24 | 2024-02-27 | 南方科技大学 | 基于等效慢度的快速推进法的地震波走时求解方法及系统 |
| CN117607957B (zh) * | 2024-01-24 | 2024-04-02 | 南方科技大学 | 基于等效慢度的快速推进法的地震波走时求解方法及系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106772577B (zh) | 基于微地震数据和spsa优化算法的震源反演方法 | |
| CN112328718B (zh) | 一种基于车辆动态轨迹跟踪的道路拓扑构建系统和方法 | |
| CN106885576B (zh) | 一种基于多点地形匹配定位的auv航迹偏差估计方法 | |
| Hu et al. | Acoustic emission source location and experimental verification for two-dimensional irregular complex structure | |
| CN103136393B (zh) | 一种基于网格划分的区域覆盖率计算方法 | |
| CN105241465B (zh) | 一种道路更新的方法 | |
| CN106597363A (zh) | 一种室内wlan环境下的行人定位方法 | |
| CN105093299A (zh) | 一种基于炮检距向量片技术优化观测系统的方法及装置 | |
| CN104807460A (zh) | 无人机室内定位方法及系统 | |
| CN106249295B (zh) | 一种井中微地震p、s波联合快速定位方法及系统 | |
| CN108984818A (zh) | 固定翼时间域航空电磁数据拟三维空间约束整体反演方法 | |
| CN105759311A (zh) | 一种近实时地震震源位置定位方法 | |
| CN105246153A (zh) | 一种室内指纹定位数据库高密度快速采集方法 | |
| CN106028446A (zh) | 室内停车场定位方法 | |
| CN104507097A (zh) | 一种基于WiFi位置指纹的半监督训练方法 | |
| CN104155694B (zh) | 一种反射转换横波共检波点叠加剖面的剩余静校正方法 | |
| CN108845358B (zh) | 断层及构造异常体识别方法及装置 | |
| CN116660980A (zh) | 一种基于改进程函方程的微地震定位方法 | |
| CN102706348B (zh) | 一种基于三角形的重力图快速匹配方法 | |
| CN113960532B (zh) | 一种基于假想源的二次定位计算的微地震定位方法 | |
| CN109085642B (zh) | 一种各向异性介质微地震事件定位方法 | |
| CN109212594B (zh) | 一种各向异性介质纵横波联合定位方法 | |
| CN109387872B (zh) | 表面多次波预测方法 | |
| CN113568041B (zh) | 时移地震三维拖缆采集数据的可重复性分析方法及系统 | |
| CN111189926B (zh) | 一种基于全域搜索辨识结构空洞位置的方法及系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |