CN105572722A - 一种微地震双力偶震源的加载方法 - Google Patents
一种微地震双力偶震源的加载方法 Download PDFInfo
- Publication number
- CN105572722A CN105572722A CN201410542783.6A CN201410542783A CN105572722A CN 105572722 A CN105572722 A CN 105572722A CN 201410542783 A CN201410542783 A CN 201410542783A CN 105572722 A CN105572722 A CN 105572722A
- Authority
- CN
- China
- Prior art keywords
- partiald
- delta
- tau
- force
- lambda
- 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.)
- Granted
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
Abstract
本发明提供了一种微地震双力偶震源的加载方法,属于地震勘探基础应用领域。所述方法基于波动方程有限差分数值模拟实现微地震双力偶源的加载;所述方法包括:(1)确定断层面参数:走向Φ、倾角δ和滑动角λ;(2)根据所述断层面参数建立地震矩张量;(3)由所述地震矩张量的各个元素值获得双力偶源的加载方式;(4)通过集中力组合加载方式完成双力偶源的加载;(5)完成双力偶源波场模拟。
Description
技术领域
本发明属于地震勘探基础应用领域,具体涉及一种微地震双力偶震源的加载方法,通过模拟结果来认识微地震压裂波场特征和微地震震源形成机理,最终给微地震震源的反演问题提供技术指导。
背景技术
微地震监测技术通过对储层压裂引发的微地震震源成像可以显示压裂裂隙的延伸方向、高度、长度、不对称性等属性,有助于帮助油藏开发人员监测压裂施工效果,优化压裂施工设计、调整开发及注水井网部署,微地震监测技术是非常规油气资源特别是页岩气开发的关键技术之一。
对于微地震震源机理的研究经历了很长的时间,1923年日本一名学者首先提出了震源的单力偶力系,即在地震瞬间,震源处突然作用一个力偶,使断层两盘发生相对运动,扰动周围介质,辐射出地震波,用作用于震源处的一些集中力系来解释震源辐射地震波的特征。此后,日本另一位学者又提出双力偶力系,若在一个小的平面断层上发生一个突然的纯剪切错动,则会产生地震波辐射,这样的剪切错动震源产生的远场地震波与在震源处突然有一个双力偶的作用产生的地震波相同,由于地震学的震源理论与事实证明双力偶力系比较接近实际,因此现在比较常用的震源模型为双力偶力系点源模型。
为了监测水力压裂和弄清楚与地下介质有关的地震响应,通过加载双力偶源进行正演模拟生成合成数据来验证野外地震数据和在油藏范围内监测流体过程,运用模拟结果来分析不同双力偶源参数对应的波场特征,为后续的微地震震源机制解(运用实际资料反演地层压裂产状:震源等级,走向倾角(δ)和滑动角(λ))提供理论指导,因此实现基于波动方程有限差分数值模拟的微地震双力偶源的加载非常有意义。
发明内容
本发明的目的在于解决上述现有技术中存在的难题,提供一种微地震双力偶震源的加载方法,针对微地震特殊的震源机制,运用波动方程有限差分数值模拟技术,实现双力偶震源的加载,并且给出模拟结果,在此基础上,进行不同双力偶源参数的数值模拟,统计其对应的波场特征,为实际资料微地震震源机制解(震源等级,走向倾角(δ)和滑动角(λ))提供理论指导,震源参数的确定就可以确定裂缝的走向和方位,进一步可以确定压裂的有效体积,最终达到完成产量的评估作用。
本发明是通过以下技术方案实现的:
一种微地震双力偶震源的加载方法,基于波动方程有限差分数值模拟实现微地震双力偶源的加载;
所述方法包括:
(1)确定断层面参数:走向Φ、倾角δ和滑动角λ;
(2)根据所述断层面参数建立地震矩张量;
(3)由所述地震矩张量的各个元素值获得双力偶源的加载方式;
(4)通过集中力组合加载方式完成双力偶源的加载;
(5)完成双力偶源波场模拟。
所述步骤(1)中,走向Φ是由正北至断层线顺时针量取的角度(左脚在下盘,右脚在上盘),范围为0°至360°;倾角δ是断层面与地平面间的夹角,范围为0°至90°;滑动角λ是以断层走向为基础,向上为正,向下为负,范围为-180°至+180°。
所述步骤(2)是这样实现的:
地震矩张量为3×3形式的矩阵,共有9个元素,表示如下:
其中,Mij表示一个力偶,M11表示作用于1轴、平行于1轴的力偶;M12表示作用于2轴、平行于1轴的力偶;M13表示作用于3轴、平行于1轴的力偶;其他元素以此类推;
所述地震矩张量与断层面参数的表达式如下:
M11=-M0(sinδcosλsin2φ+sin2δsinλsin2φ)
M22=M0(sinδcosλsin2φ-sin2δsinλcos2φ)
M33=M0sin2δsinλ=-(M11+M22)
M13=-M0(cosδcosλcosφ+cos2δsinλsinφ)
M23=-M0(cosδcosλsinφ-cos2δsinλcosφ)
其中,M0表示地震力矩。
由于具有对称性(满足角动量守恒),因此只需这6个独立矩张量。
所述步骤(3)是这样实现的(针对二维情形-xoz平面内):
以水平方向和垂直方向间隔分别为dx和dz对速度介质进行网格状剖分,得到水平和垂直网格数为nx和nz;
选定震源点(nsx,nsz),则双力偶的四个力作用的网格分别为(nsx-1,nsz)、(nsx+1,nsz)、(nsx,nsz-1)和(nsx,nsz+1);
在所述步骤(2)得到的地震矩张量的各个元素值中,如果元素值不等于0,则表示这个作用力存在,如果元素值等于0,则表示此作用力不存在;根据各个元素的意义获取元素值不等于0的元素对应的位置和方向。
所述步骤(4)是这样实现的:
对步骤(3)得到的4个网格点加载作用力函数,具体如下:
从步骤(3)确定的4个网格点中找到元素值不等于0的网格点的位置,这些即为加载作用力的网格点位置,在每一个加载作用力的网格点位置处加载一个作用力,在作用力规定时间内,其值的大小都是所述作用力函数值。
所述步骤(5)是这样实现的:
采用基于波动方程有限差分数值模拟得到地震波场记录,具体步骤如下:
A)获得各向同性介质中的二维一阶速度-应力弹性波方程,即波动方程,如下所示;
式中:分别为质点振动速度在x和z方向上的分量;uz、uz分别为位移u在x和z方向上的分量;τxx和τzz为质点在x和z方向上的正应力;τxz为质点在xz平面内的剪切力;ρ为介质密度;
B)基于所述波动方程完成交错网格有限差分算法;
结合交错网格有限差分算法,得到所述波动方程的2N阶空间差分精度、二阶时间差分精度交错网格高阶有限差分格式,即
式中:表示沿x方向做向前差分;表示沿x方向做向后差分;Δx,Δz表示x,z方向的网格间距;Δt表示时间步长;其它表示形式以此类推。
所述作用力函数使用雷克子波函数。
与现有技术相比,本发明的有益效果是:本发明利用断层面的参数得到地震矩张量,进而得到双力偶源的加载方式,通过不同双力偶源的模拟,进行微地震震源机理的研究;通过对波场的模拟,进行微地震波场特征研究;为实际微地震震源机制解的反演问题提供技术指导,能更好地进行压裂有效裂缝的评估,效果明显。
附图说明
图1为地震矩张量各元素的示意图;
图2为断层面的走向、倾角及滑动角示意图;
图3(a)为垂直作用力
图3(b)为水平作用力
图3(c)为倾斜作用力;
图4为双力偶源加载示意图;
图5为双力偶震源离散网格加载方式;
图6a双力偶源模拟结果照明的Z分量
图6b为快照的Z分量;
图6c为离地面中心偏移距为500m的单道记录的Z分量;
图7本发明方法的步骤框图。
图8弹性波交错网格差分示意图
具体实施方式
下面结合附图对本发明作进一步详细描述:
(1)地震矩张量的建立
在远场和点源近似的情况下,在观测点x处质点位移可以简化为:
un(x,t)=Mkl*Gnk,l(1)
式中un(x,t)为质点位移函数、Mkl为地震矩张量、Gnk,l为格林函数。星号为两个时间函数的褶积运算,其中:
式中mkld∑具有力矩的量纲。
地震矩张量一般可以表示为3×3形式的矩阵,共有9个元素,各元素的意义如图1所示:
在各向同性介质中法向单位矢量为n、面积为S的断层上,发生的滑动(位错失量)D的位错源所对应的地震矩为:
Mkl={λD·nδkl+μ[Dlnk+Dknl]}S(4)
假设D矢量限定在沿断层面任意方向,并设断层面的走向、倾角和滑动角分别为φ、δ、λ(如图2所示),则有:
n=(-sinδsinφ,sinδcosφ,-cosδ)(5)
结合式(4)、(5)和(6)可以得到矩张量跟断层面参数的一般表达式,如下所示:
M11=-M0(sinδcosλsin2φ+sin2δsinλsin2φ)
M22=M0(sinδcosλsin2φ-sin2δsinλcos2φ)
M33=M0sin2δsinλ=-(M11+M22)
M13=-M0(cosδcosλcosφ+cos2δsinλsinφ)
M23=-M0(cosδcosλsinφ-cos2δsinλcosφ)(7)
由(7)式可知,断层面的参数决定震源地震矩张量的结果,因此知道地震矩张量也是可以确定断层面的参数,进而可以预测裂缝的方位。
(2)微地震双力偶源的模拟
1)集中力震源的加载实现
本发明主要是基于集中力震源的加载方式,结合断层面参数对应的地震矩张量进行双力偶源的实现。集中力源模拟是施加一个随时间变化的作用力。集中力源具有方向性,可以是垂直作用力也可以是水平作用力,还可以是倾斜作用力。集中力源在弹性介质中激发的弹性波场具有方向效应。集中力源的加载方式如图3(a)、图3(b)和图3(c)所示。
2)微地震双力偶源加载实现
运用上述集中力源的加载方式进行组合加载就可以得到(作用于xoz平面内)微地震中关于双偶力源的实现(作用方式见图4)。
如图7所示,详细实现步骤如下:
先以水平方向和垂直方向间隔分别为dx和dz(满足水平和垂直方向剖分后的网格数为整数即可)对介质(指的是速度介质;在剖分前不需要进行处理;需要纵波速度值)进行剖分成网格状,得到水平和垂直网格数为nx和nz,然后选定震源点(nsx,nsz)(根据模拟需要(地面或地下激发等)选择震源点位置,进而确定nsx、nsz的数值),双力偶(四个力)作用的网格分别为(nsx-1,nsz)、(nsx+1,nsz)、(nsx,nsz-1)和(nsx,nsz+1),最后按照这几个网格点加载(如图5所示,图5中给出了基于图4的作用方式完成双力偶源加载)作用力函数(可以使用的作用力函数有很多种,常用的就是雷克子波函数来表示作用力函数;具体加载方法就是:首先确定加载作用力的网格点位置,每一个位置处都要加载一个作用力(Fx表示水平方向,向右为正,向左为负;Fz表示垂直方向,向下为正,向上为负),在作用力规定时间内,其值大小都是雷克子波函数值)(见图5),进行数值模拟就能得到双力偶源的结果。
采用基于波动方程有限差分进行数值模拟,具体步骤如下:
A)获得各向同性介质中的二维一阶速度-应力弹性波方程,如下所示;
式中:分别为质点振动速度在x和z方向上的分量;uz、uz分别为位移u在x和z方向上的分量;τxx和τzz为质点在x和z方向上的正应力;τxz为质点在xz平面内的剪切力;ρ为介质密度;
B)基于上述波动方程完成交错网格有限差分算法;
各种波场分量和物性参数按照图8分布,再结合交错网格的差分思路,可以得到上述方程的2N阶空间差分精度、二阶时间差分精度交错网格高阶有限差分格式,即
式中:表示沿x方向做向前差分;表示沿x方向做向后差分;Δx,Δz表示x,z方向的网格间距;Δt表示时间步长;其它表示形式同理。
下面是本发明的一个实施例:
介质模型参数:模型网格数500*500,网格大小5m,主频40Hz,纵波速度为4000m/s,横波速度为2200m/s,密度为1000kg·m-3;断层面的走向、倾角和滑动角分别为270°、90°和270°,则对应的地震矩张量为M=[001;000;100],则对应图1中的M13和M31作用力存在,因此加载的位置和方向如图5所示,震源力函数为主频为40Hz的雷克子波。最后得到的结果如图6a-6c所示,按照此震源方式进行数值模拟就可以得到微地震记录,进行波场特征分析工作。
加载位置和方向的确定方法如下:
1)首先定下波场模拟时,震源点位置(即nsx和nsz的值,这个好理解);
2)再来确定双力偶源的加载位置,首先由上述的震源点位置的周围相邻位置确定双力偶源加载的大概位置范围(见图5示意图所示),具体位置和方向由矩张量确定;
3)由断层面参数计算得到矩张量的9个元素值,只要元素值不等于0,那么表示这个作用力存在,只要为0,就表示此作用力不存在;根据不等于0的元素对照图1所示的各元素示意图,来确定哪些力存在,在什么位置,方向如何(本专利中的实例可知,地震矩张量为M=[001;000;100],则对应图1中的M13和M31作用力存在,得到图4作用力示意图,再由图4确定图5加载方式。
上述技术方案只是本发明的一种实施方式,对于本领域内的技术人员而言,在本发明公开了应用方法和原理的基础上,很容易做出各种类型的改进或变形,而不仅限于本发明上述具体实施方式所描述的方法,因此前面描述的方式只是优选的,而并不具有限制性的意义。
Claims (7)
1.一种微地震双力偶震源的加载方法,其特征在于:所述方法基于波动方程有限差分数值模拟实现微地震双力偶源的加载;
所述方法包括:
(1)确定断层面参数:走向Φ、倾角δ和滑动角λ;
(2)根据所述断层面参数建立地震矩张量;
(3)由所述地震矩张量的各个元素值获得双力偶源的加载方式;
(4)通过集中力组合加载方式完成双力偶源的加载;
(5)完成双力偶源波场模拟。
2.根据权利要求1所述的微地震双力偶震源的加载方法,其特征在于:所述步骤(1)中,走向Φ是由正北至断层线顺时针量取的角度,范围为0°至360°;倾角δ是断层面与地平面间的夹角,范围为0°至90°;滑动角λ是以断层走向为基础,向上为正,向下为负,范围为-180°至+180°。
3.根据权利要求2所述的微地震双力偶震源的加载方法,其特征在于:所述步骤(2)是这样实现的:
地震矩张量为3×3形式的矩阵,共有9个元素,表示如下:
其中,Mij表示一个力偶,M11表示作用于1轴、平行于1轴的力偶;M12表示作用于2轴、平行于1轴的力偶;M13表示作用于3轴、平行于1轴的力偶;其他元素以此类推;
所述地震矩张量与断层面参数的表达式如下:
M11=-M0(sinδcosλsin2φ+sin2δsinλsin2φ)
M22=M0(sinδcosλsin2φ-sin2δsinλcos2φ)
M33=M0sin2δsinλ=-(M11+M22)
M13=-M0(cosδcosλcosφ+cos2δsinλsinφ)
M23=-M0(cosδcosλsinφ-cos2δsinλcosφ)
其中,M0表示地震力矩。
4.根据权利要求3所述的微地震双力偶震源的加载方法,其特征在于:所述步骤(3)是这样实现的:
以水平方向和垂直方向间隔分别为dx和dz对速度介质进行网格状剖分,得到水平和垂直网格数为nx和nz;
选定震源点(nsx,nsz),则双力偶的四个力作用的网格分别为(nsx-1,nsz)、(nsx+1,nsz)、(nsx,nsz-1)和(nsx,nsz+1);
在所述步骤(2)得到的地震矩张量的各个元素值中,如果元素值不等于0,则表示这个作用力存在,如果元素值等于0,则表示此作用力不存在;根据各个元素的意义获取元素值不等于0的元素对应的位置和方向。
5.根据权利要求4所述的微地震双力偶震源的加载方法,其特征在于:所述步骤(4)是这样实现的:
对步骤(3)得到的4个网格点加载作用力函数,具体如下:
从步骤(3)确定的4个网格点中找到元素值不等于0的网格点的位置,这些即为加载作用力的网格点位置,在每一个加载作用力的网格点位置处加载一个作用力,在作用力规定时间内,其值的大小都是所述作用力函数值。
6.根据权利要求5所述的微地震双力偶震源的加载方法,其特征在于:所述步骤(5)是这样实现的:
采用基于波动方程有限差分数值模拟得到地震波场记录,具体步骤如下:
A)获得各向同性介质中的二维一阶速度-应力弹性波方程,即波动方程,如下所示;
式中:分别为质点振动速度在x和z方向上的分量;ux、uz分别为位移u在x和z方向上的分量;τxx和τzz为质点在x和z方向上的正应力;τxz为质点在xz平面内的剪切力;ρ为介质密度;
B)基于所述波动方程完成交错网格有限差分算法;
结合交错网格有限差分算法,得到所述波动方程的2N阶空间差分精度、二阶时间差分精度交错网格高阶有限差分格式,即
式中:表示沿x方向做向前差分;表示沿x方向做向后差分;Δx,Δz表示x,z方向的网格间距;Δt表示时间步长;其它表示形式以此类推。
7.根据权利要求5所述的微地震双力偶震源的加载方法,其特征在于:作用力函数使用雷克子波函数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410542783.6A CN105572722B (zh) | 2014-10-14 | 2014-10-14 | 一种微地震双力偶震源的加载方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410542783.6A CN105572722B (zh) | 2014-10-14 | 2014-10-14 | 一种微地震双力偶震源的加载方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105572722A true CN105572722A (zh) | 2016-05-11 |
CN105572722B CN105572722B (zh) | 2017-10-20 |
Family
ID=55883063
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410542783.6A Active CN105572722B (zh) | 2014-10-14 | 2014-10-14 | 一种微地震双力偶震源的加载方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105572722B (zh) |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5183109A (en) * | 1991-10-18 | 1993-02-02 | Halliburton Company | Method for optimizing hydraulic fracture treatment of subsurface formations |
WO2010078577A2 (en) * | 2009-01-05 | 2010-07-08 | Services Petroliers Schlumberger | Processing time series data embedded in high noise |
GB2503903A (en) * | 2012-07-11 | 2014-01-15 | Schlumberger Holdings | Monitoring and characterising fractures in a subsurface |
CN103926620A (zh) * | 2014-05-09 | 2014-07-16 | 南京大学 | 基于阵列反褶积处理的水力压裂监测方法 |
-
2014
- 2014-10-14 CN CN201410542783.6A patent/CN105572722B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5183109A (en) * | 1991-10-18 | 1993-02-02 | Halliburton Company | Method for optimizing hydraulic fracture treatment of subsurface formations |
WO2010078577A2 (en) * | 2009-01-05 | 2010-07-08 | Services Petroliers Schlumberger | Processing time series data embedded in high noise |
GB2503903A (en) * | 2012-07-11 | 2014-01-15 | Schlumberger Holdings | Monitoring and characterising fractures in a subsurface |
CN103926620A (zh) * | 2014-05-09 | 2014-07-16 | 南京大学 | 基于阵列反褶积处理的水力压裂监测方法 |
Non-Patent Citations (1)
Title |
---|
JAMES T.RUTLEDGE ETC.: "《HYDRAAULIC STIMULATION OF NATURAL FRACTURES AS REVEALED BY INDUCED MICROEARTHQUAKES CARTHAGE COTTON VALLEY GAS FIEDL, EAST TEXAS》", 《地球物理学》 * |
Also Published As
Publication number | Publication date |
---|---|
CN105572722B (zh) | 2017-10-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Wang et al. | Parametric models for 3 D t o p o g r a p h i c a m p l i f i c a t i o n o f ground motions considering subsurface soils | |
Quintal et al. | Quasi‐static finite element modeling of seismic attenuation and dispersion due to wave‐induced fluid flow in poroelastic media | |
Masson et al. | Finite-difference modeling of Biot’s poroelastic equations across all frequencies | |
US20140129479A1 (en) | Method to aid in the exploration, mine design, evaluation and/or extraction of metalliferous mineral and/or diamond deposits | |
Ladopoulos | Elastodynamics for Non-linear Seismic Wave Motion in Real-Time Expert Seismology | |
Lee et al. | Dynamic analysis of a layered half-space subjected to moving line loads | |
Stabile et al. | A comprehensive approach for evaluating network performance in surface and borehole seismic monitoring | |
van Wees et al. | 3-D mechanical analysis of complex reservoirs: a novel mesh-free approach | |
Primofiore et al. | 3D numerical modelling for interpreting topographic effects in rocky hills for Seismic Microzonation: The case study of Arquata del Tronto hamlet | |
Mühlhaus et al. | Large amplitude folding in finely layered viscoelastic rock structures | |
Homma et al. | A physics-based Monte Carlo earthquake disaster simulation accounting for uncertainty in building structure parameters | |
Sidler et al. | A pseudo-spectral method for the simulation of poro-elastic seismic wave propagation in 2D polar coordinates using domain decomposition | |
Ronchin et al. | Solid modeling techniques to build 3D finite element models of volcanic systems: an example from the Rabaul Caldera system, Papua New Guinea | |
Burjánek et al. | Modeling the seismic response of unstable rock mass with deep compliant fractures | |
CN105676280A (zh) | 基于旋转交错网格的双相介质地质数据获取方法和装置 | |
Wu et al. | Two-dimensional finite-difference seismic modeling of an open fluid-filled fracture: Comparison of thin-layer and linear-slip models | |
Abul Khair et al. | The effect of present day in situ stresses and paleo-stresses on locating sweet spots in unconventional reservoirs, a case study from Moomba-Big Lake fields, Cooper Basin, South Australia | |
Li et al. | Elastic wave finite-difference simulation using discontinuous curvilinear grid with non-uniform time step: two-dimensional case | |
Sadovskii et al. | Analysis of Seismic Waves Excited in Near-Surface Soils by Means of the Electromagnetic Pulse Source ‘Yenisei’ | |
Anandarajah et al. | Calibration of dynamic analysis methods from field test data | |
Vasco et al. | On the use of adjoints in the inversion of observed quasi-static deformation | |
Corona-Romero et al. | Simulation of LP seismic signals modeling the fluid–rock dynamic interaction | |
Dujardin et al. | Simulation of the basin effects in the Po Plain during the Emilia-Romagna seismic sequence (2012) using empirical Green’s functions | |
CN105572722A (zh) | 一种微地震双力偶震源的加载方法 | |
Morozov et al. | Simulating the state of stress and strain in the epicentral zone of a large earthquake in Turkey (Izmit, 1999, M 7.4) |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |