CN107728200B - 地面微地震压裂裂缝动态展布实时监测方法 - Google Patents

地面微地震压裂裂缝动态展布实时监测方法 Download PDF

Info

Publication number
CN107728200B
CN107728200B CN201710915453.0A CN201710915453A CN107728200B CN 107728200 B CN107728200 B CN 107728200B CN 201710915453 A CN201710915453 A CN 201710915453A CN 107728200 B CN107728200 B CN 107728200B
Authority
CN
China
Prior art keywords
microseism
seismic
data
micro
real
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.)
Active
Application number
CN201710915453.0A
Other languages
English (en)
Other versions
CN107728200A (zh
Inventor
尚新民
刁瑞
芮拥军
单联瑜
王常波
韩宏伟
王兴谋
李红梅
崔庆辉
李强
王东凯
胡晓婷
赵翠霞
揭景荣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Geophysical Research Institute of Sinopec Shengli Oilfield Co
Original Assignee
China Petroleum and Chemical Corp
Geophysical Research Institute of Sinopec Shengli Oilfield Co
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, Geophysical Research Institute of Sinopec Shengli Oilfield Co filed Critical China Petroleum and Chemical Corp
Priority to CN201710915453.0A priority Critical patent/CN107728200B/zh
Publication of CN107728200A publication Critical patent/CN107728200A/zh
Application granted granted Critical
Publication of CN107728200B publication Critical patent/CN107728200B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/288Event detection in seismic signals, e.g. microseismics
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/40Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
    • G01V1/44Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging using generators and receivers in the same well
    • G01V1/48Processing data
    • G01V1/50Analysing data
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/10Aspects of acoustic signal generation or detection
    • G01V2210/12Signal generation
    • G01V2210/123Passive source, e.g. microseismics
    • G01V2210/1234Hydrocarbon reservoir, e.g. spontaneous or induced fracturing
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/10Aspects of acoustic signal generation or detection
    • G01V2210/14Signal detection
    • G01V2210/142Receiver location
    • G01V2210/1425Land surface
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/10Aspects of acoustic signal generation or detection
    • G01V2210/14Signal detection
    • G01V2210/142Receiver location
    • G01V2210/1429Subsurface, e.g. in borehole or below weathering layer or mud line
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/10Aspects of acoustic signal generation or detection
    • G01V2210/16Survey configurations
    • G01V2210/161Vertical seismic profiling [VSP]
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/10Aspects of acoustic signal generation or detection
    • G01V2210/16Survey configurations
    • G01V2210/163Cross-well
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/64Geostructures, e.g. in 3D data cubes
    • G01V2210/646Fractures

Landscapes

  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Environmental & Geological Engineering (AREA)
  • Acoustics & Sound (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Business, Economics & Management (AREA)
  • Emergency Management (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供一种地面微地震压裂裂缝动态展布实时监测方法,包括:获取连续微地震数据,将压裂区域网格剖分并建立块状高精度速度模型;计算剖分网格点的旅行时间、方位角度和出射角度;采用谱分解方法进行分频段处理,计算三个频段数据的三分量波形协方差矩阵及线偏振特征;根据线偏振特征对分频段数据进行极化滤波;任选剖分网格点,对极化滤波后的微地震数据进行矢量偏移成像,获取该网格点的成像能量值;将梯度下降值大于一预先设定阈值的区域视为裂缝密集区,得到该时间段内的裂缝展布图像。该地面微地震压裂裂缝动态展布实时监测方法对地表微地震监测具有很好的效果,可以实现监测微裂缝的动态发展。

Description

地面微地震压裂裂缝动态展布实时监测方法
技术领域
本发明涉及油气田开发技术领域,特别是涉及到一种地面微地震压裂裂缝动态展布实时监测方法。
背景技术
非常规储层的压裂改造是提高非常规油气采收率的关键技术,微地震技术是监测非常规油气开采过程中储层压裂改造效果的主要技术,微地震监测技术能够实现裂缝的动态监测。对于储层压裂改造过程中的微地震监测而言,关键环节是准确识别微地震事件,并准确定位微地震事件的位置,从而进一步监测微地震裂缝发生、发展的动态过程。
申请号为CN94105769的中国专利申请液压致裂的地层裂缝监测系统中,提到类似的裂缝监测方法,它在液压井孔周围设置的至少三口监测井孔的套管地面端分别装设声波检测器采集地层压裂时的声发射信号,经信号处理得到各路信号接收时差后由计算机判定压裂点的平面坐标并绘制分布图,显示出地层裂缝的方位与形态。该发明具有节省费用、提高效率、数据可靠等诸多优点。但是该申请是在井中进行微地震监测,井中监测具有选井困难、监测费用高、施工要求高等不利因素,地面微地震监测具有一定的优势。
申请号为CN201010151905的中国专利申请油井压裂微地震地面—井中裂缝监测方法,提到类似的裂缝监测方法:利用弹性波场在地层介质中的传播特性,微震波被布置在井口周围的监测分站收到,根据各分站微震波的到达时差,会形成一系列的方程组,求解这一系列方程组,就可确定微震震源位置,进而确定出裂缝分布的方位、长度、高度及地应力方向等地层参数。井中的背景噪声小,通过井下的观测信号对地面的观测信号进行滤波处理,将压裂微震信号和地面干扰进行筛选区分,放置在井中的检波器是整套系统采集信号的一个基准值,用这个基准值来对地面拾震器所采集到的信号进行识别,滤除干扰信号,使得监测的结果更加的真实可靠。该申请采用地面—井中联合的方式进行裂缝监测,通过井下观测信号对地面观测信号进行滤波处理,这就要求必须具备井中和地面两种观测数据,对施工提出了更高的要求,更强调地面和井中数据的联合使用。
为了便于施工,同时提高裂缝监测精度,通常在地表布置微地震监测台站,微地震信号在近地表沉积层内衰减很快,造成了微地震信号的高频剧烈衰减,但是常规地面微地震监测台站装置一般采用10Hz模拟检波器,在自然频率以下振幅响应快速下降,不利于低频端有效信号的接收。在地表布置天然地震台站或具有自主知识产权的宽频检波器(已获实用新型专利授权,授权号:ZL 201520971428.0),以此来监测致密油气藏水力压裂是一种新的尝试,监测频率带宽一般在1Hz–200Hz左右,能够有效接收到微地震低频信号,弥补了常规自然频率10Hz检波器的不足。
地表布置天然地震台站或宽频检波器进行微地震监测,具有:观测系统设计有针对性、台站布设灵活、有效监测低频信号、安装简单、监测成本较低等优点。但是,在地面进行微地震监测时,易于受到地表噪音、人为干扰的影响,在地表很难记录到高信噪比的微地震信号,并且,压裂改造过程中产生的绝大部分微地震事件属于极其微弱的信号。如何从地表微地震数据中识别、提取微弱的微地震信号,并用来监测压裂改造过程中人工裂缝的动态,是地面微地震监测技术的关键。传统基于能量扫描的方法难以准确的识别微弱能量的微地震事件,无法实现微地震事件的准确定位,因此需要发展新的方法来识别微弱信号。为此我们发明了一种新的地面微地震压裂裂缝动态展布实时监测方法,解决了以上技术问题。
发明内容
本发明的目的是提供一种实现了微弱信号的有效识别和准确定位,达到了微地震监测压裂裂缝动态展布目的的地面微地震压裂裂缝动态展布实时监测方法。
本发明的目的可通过如下技术措施来实现:地面微地震压裂裂缝动态展布实时监测方法,该地面微地震压裂裂缝动态展布实时监测方法包括:步骤1,采用优选观测方式获取连续微地震数据,将压裂区域网格剖分并建立块状高精度全三维速度模型;步骤2,结合速度模型各向异性及震源机制解,计算剖分网格点的旅行时间、方位角度和出射角度;步骤3,采用谱分解方法进行微地震数据分频段自动化处理,计算三个频段数据的三分量波形协方差矩阵及线偏振特征;步骤4,根据线偏振特征分别对分频段数据进行极化滤波,然后将分频段数据进行合并;步骤5,任选剖分网格点,对极化滤波后的微地震数据进行矢量偏移成像,获取该网格点的成像能量值;步骤6,将梯度下降值大于一预先设定阈值的区域视为裂缝密集区,得到该时间段内的裂缝展布图像。
本发明的目的还可通过如下技术措施来实现:
在步骤1中,在压裂区域上方地表或浅井中,采用分散式布设台站仪器,以获取连续微地震数据。
在步骤1中,在储层水力压裂区域的上方地表或浅井中,采用均匀、对称、充分、连续的方式,布设天然地震台站仪器或宽频检波器,以获取致密储层水力压裂改造的x,y,z三分量微地震宽频数据。
采用高分辨率、宽频带的三分量天然地震台站仪器或宽频检波器进行观测方式布设时:一方面,以地面微地震观测系统最优化设计为指导思路,采用均匀、对称、充分、连续的最佳观测方式布设,并采取挖坑和埋置措施,降低地表噪音的干扰,所布设台站监测孔径要大于压裂段垂深的1.5倍以上;另一方面,在均匀、对称、充分、连续的基础上,根据该地区地质构造、速度模型与震源机制的压裂段能量传播特征,在分散式布设台站仪器时,保证微地震监测能量在起伏地表的均匀分布,以获取高品质的连续微地震数据。
在步骤1中,将储层水力压裂区域范围进行网格点剖分,并结合三维地面地震、井中垂直地震和测井资料,建立压裂区域的块状精细速度模型,并采用块状精细速度模型方法,实现全三维速度的逐渐变化,全三维速度模型包含:起伏真地表、复杂低速层、渐变降速带、平滑高速层顶界面、各项异性中深层速度模型。
在步骤2中,结合速度模型各向异性及震源机制解,基于射线高斯束追踪方法,采用速度逐渐变化的块状精细速度模型,计算每一个剖分网格点经过块状全三维速度模型后,到达地表或浅井中台站的旅行时间、方位角度和出射角度。
在步骤3中,选取一个台站的三分量微地震宽频数据,采用高精度谱分解方法,将宽频微地震数据自动化分解为三个频段的微地震数据,分别计算三个频段数据的三分量波形协方差矩阵,分别利用协方差矩阵的特征值计算三个频段微地震数据的线偏振特征。
在步骤3中,采用高精度谱分解方法,将宽频微地震数据分解为三个频段的微地震数据,三个频段的微地震数据分别是微地震事件可靠频段、微地震事件优势频段、微地震事件劣势频段,在高精度谱分解之后,计算不同频段数据的三分量波形协方差矩阵,利用协方差矩阵的特征值计算频段微地震数据的线偏振特征,利用线偏振特征进行极化滤波处理,增强微地震事件,提高微地震数据的信噪比。
在步骤3中,在高精度谱分解之后,计算不同频段数据的三分量波形协方差矩阵,计算一段时间内的分频段数据的协方差矩阵时间采样点总数为N,计算公式为:
其中:
L=(N-1)/2,L是时间窗口长度的一半;为微地震三分量的分频段数据,l,m=1,2,3;计算协方差矩阵的特征值(λ123)和相应的特征向量(V1,V2,V3);利用(λ123)计算该时刻内三分量微地震记录的线偏振特征;
采用全局偏振参数Mτ或线性偏振度MRL来表示,
其中:MQ决定MRL对不同偏振程度的敏感度,MQ≤1。
在步骤4中,根据线偏振特征对三个频段,即微地震事件可靠频段、微地震事件优势频段、微地震事件劣势频段的微地震数据进行极化滤波处理,然后利用谱分解原理将分频段数据合并。
在步骤4中,根据线偏振特征对三个频段的微地震数据进行极化滤波处理,极化滤波表示为:
A(i)取线偏振特征Mτ或MRL,或者其它表示线偏振特征的量,为微地震三分量的分频段数据,其中:i=1,2,3,...,N,为微地震数据的离散采样;k=1,2,3,...,M,为微地震台站序号;j=1,2,3,共计3个分量,即两个水平分量和一个垂直分量。
该地面微地震压裂裂缝动态展布实时监测方法还包括,在步骤4之后,重复步骤3和4,对所有台站的微地震宽频数据进行相同处理,获取所有微地震台站仪器的分频段极化滤波后微地震宽频数据。
在步骤5中,选取一个剖分网格点,根据该剖分网格点到不同台站仪器的旅行时间、方位角度和出射角度,对所有台站仪器的微地震宽频数据进行分频段极化滤波处理,并对分频段极化滤波后的微地震宽频数据进行矢量偏移成像,获取该剖分网格点的偏移成像能量值。
在步骤5中,选取剖分网格点,将宽频微地震数据分解为三个频段的微地震数据,对所有台站仪器的微地震宽频数据进行分频段极化滤波处理,并对分频段极化滤波后的微地震宽频数据进行矢量偏移成像,对P波进行偏移成像,或对S波进行偏移成像,或将P波和S波进行联合偏移成像。
在步骤5中,网格点x的能量值利用公式(6)表示:
其中:c=P or S,τ为不同台站所对应的旅行时间,L=(N-1)/2,L是时间窗口长度的一半,M为台站仪器个数,wc为权重值,为微地震数据的极化滤波结果,j=1,2,3,共计3个分量;令θ为出射角度,φ为方位角度,则射线矢量可记为:
Vr=(sinθ,cosθcosφ,cosθsinφ) (7)
权重值wc的计算公式为:
其中:Vr为射线矢量,V1为协方差矩阵的特征向量。
该地面微地震压裂裂缝动态展布实时监测方法还包括,在步骤5之后,重复步骤5,计算所有剖分网格点在某一时间段内的偏移成像能量值,并计算该时间段内的能量梯度下降值。
在步骤6中,在压裂区域预先设定阈值,当能量梯度下降值大于该阈值时,表明该区域范围内为压裂裂缝的密集区,从而得到该时间段内的压裂裂缝展布图像。
该地面微地震压裂裂缝动态展布实时监测方法还包括,在步骤6之后,沿采集时间段进行处理,对后续时间段的地面微地震宽频数据重复步骤3-6,计算不同时间段内压裂裂缝展布图像,监测压裂裂缝动态展布的发展情况。
本发明中的地面微地震压裂裂缝动态展布实时监测方法,采用天然地震台站或宽频检波器,采集到x,y,z三分量数据,根据微地震数据中波形相似特性,通过极化滤波方法、纵横波联合应用、偏移成像和梯度下降数学原理,对已有常规方法进行了改进和完善,实现了微弱信号的有效识别和准确定位,达到了微地震监测压裂裂缝动态展布的目的,具有重要的意义和良好的应用前景。本发明的目的在于克服地面微地震监测中出现的问题,高分辨率的天然地震台站仪器或宽频检波器能够有效监测低频信号,基于三分量微地震数据、谱分解分频段极化滤波方法、纵横波联合应用、偏移成像和梯度下降技术,即可使用P波或S波数据,可以有效获取极低信噪比的微地震信号,克服了地表监测微地震数据信噪比低的问题,有效实现了储层压裂中裂缝的动态监测。
高精度谱分解方法将宽频微地震数据自动化分解为三个频段的微地震数据,分别是微地震事件可靠频段、微地震事件优势频段、微地震事件劣势频段,三个频段的微地震数据中微地震有效信号的占比是不同的,在后续的计算过程中赋予不同的比例因子。在高精度谱分解之后,计算不同频段数据的三分量波形协方差矩阵,利用协方差矩阵的特征值计算频段微地震数据的线偏振特征,利用线偏振特征进行极化滤波处理,可大幅增强微地震事件,同时提高微地震数据的信噪比,有助于微地震事件的准确识别。
根据实施例方法,该技术在压裂施工之前建立剖分网格信息数据库,大大提高了计算效率,在监测过程中由计算机自行工作,可以实时监测压裂改造中裂缝的发生、发展,实时评价压裂效果。对地表微地震监测具有很好的效果,可以实现监测微裂缝的动态发展,该方法并适用于识别微地震信号,同时适用于地表、浅井、深井或井中监测的微地震三分量数据。
附图说明
图1为本发明的地面微地震压裂裂缝动态展布实时监测方法的一具体实施例的流程图;
图2为本发明的一具体实施例中三分量地面微地震宽频数据进行谱分解分频段极化滤波处理前、后的对比效果图;
图3为本发明的一具体实施例中地面微地震压裂裂缝动态展布监测效果图。
具体实施方式
为使本发明的上述和其他目的、特征和优点能更明显易懂,下文特举出较佳实施例,并配合附图所示,作详细说明如下。
如图1所示,图1为本发明的地面微地震压裂裂缝动态展布实时监测方法的流程图。
在步骤101,在储层水力压裂区域的上方地表或浅井中,采用均匀、对称、充分、连续的方式,布设天然地震台站仪器或宽频检波器,以获取致密储层水力压裂改造的x,y,z三分量微地震宽频数据;采用高分辨率、宽频带的三分量天然地震台站仪器或宽频检波器,在压裂井上方地表或浅井中,以地面微地震观测系统最优化设计为指导思路,采用均匀、对称、充分、连续的最佳观测方式布设,并采取适当、有效的挖坑和埋置措施,尽可能降低地表噪音的干扰,所布设台站监测孔径要大于压裂段垂深的1.5倍以上。在均匀、对称、充分、连续的观测布设基础上,根据该地区地质构造、速度模型与震源机制的压裂段能量传播特征,在分散式布设台站仪器时,保证微地震监测能量在起伏地表的均匀分布,以获取高品质的连续微地震数据。
在步骤102,将储层水力压裂区域范围进行网格点剖分,并结合三维地面地震、VSP地震(井中垂直地震)、井间地震和测井资料,建立压裂区域的块状精细全三维速度模型;并采用块状精细速度模型方法,实现全三维速度的逐渐变化,全三维速度模型包含:起伏真地表、复杂低速层、渐变降速带、平滑高速层顶界面、各项异性中深层速度模型等。
在步骤103,结合速度模型各向异性及震源机制解,基于射线高斯束追踪方法,采用速度逐渐变化的块状精细速度模型,计算每一个剖分网格点经过块状全三维速度模型后,到达地表台站的旅行时间、方位角度和出射角度;需要压裂区域内的三维高精度速度模型,因此需要结合三维地面地震、VSP地震(井中垂直地震)、井间地震和测井资料,并采用块状精细全三维速度模型方法,实现速度的逐渐变化,避免层状介质速度模型误差大的缺点。块状精细速度建模方法摒弃了传统层状介质的假设条件,以多边形立方体为基础,形成高精度速度模型,块状速度值具备了逐渐变化的特征,更加逼近于真实的地质构造和岩性特征。
在步骤104,选取一个台站的三分量微地震宽频数据,采用高精度谱分解方法将宽频微地震数据自动化分解为三个频段的微地震数据,分别计算三个频段数据的三分量波形协方差矩阵,分别利用协方差矩阵的特征值计算三个频段微地震数据的线偏振特征;采用高精度谱分解方法,将宽频微地震数据分解为三个频段的微地震数据,三个频段的微地震数据分别是微地震事件可靠频段、微地震事件优势频段、微地震事件劣势频段,在高精度谱分解之后,计算不同频段数据的三分量波形协方差矩阵,利用协方差矩阵的特征值计算频段微地震数据的线偏振特征,利用线偏振特征进行极化滤波处理,可大幅增强微地震事件,同时提高微地震数据的信噪比,有助于微地震事件的准确识别。
在步骤105,根据线偏振特征对三个频段的微地震数据进行极化滤波处理,然后利用谱分解原理将分频段数据合并。如图2所示,利用线偏振特征进行极化滤波处理后,微地震事件的能量大幅增强,微地震数据的信噪比得到提高,有助于微地震事件的准确识别和定位精度的提高。
在步骤106,重复步骤104和105,对所有台站的微地震宽频数据进行相同处理,获取所有微地震台站仪器的分频段极化滤波后微地震宽频数据。
在步骤107,选取一个剖分网格点,根据该剖分网格点到不同台站仪器的旅行时间、方位角度和出射角度,对所有台站仪器的微地震宽频数据进行分频段极化滤波处理,并对分频段极化滤波后的微地震宽频数据进行矢量偏移成像,获取该剖分网格点的偏移成像能量值;选取剖分网格点,将宽频微地震数据分解为三个频段的微地震数据,对所有台站仪器的微地震宽频数据进行分频段极化滤波处理,并对分频段极化滤波后的微地震宽频数据进行矢量偏移成像,既可以对P波进行偏移成像,也可以对S波进行偏移成像,还可以将P波和S波进行联合偏移成像。将分频段极化滤波与P/S波进行联合应用,一方面,大幅提高微地震事件识别准确度;另一方面,通过矢量偏移成像大幅提高微地震事件定位精度。
在步骤108,重复步骤107,计算所有剖分网格点在某一时间段内的偏移成像能量值,并利用最优化方法计算该时间段内的能量梯度下降值;通过该时间段内的能量梯度下降值,实现微地震事件定位精度,通过梯度下降代替能量值,消除了异常振幅噪音的干扰,提高了微地震数据定位精度。整个处理流程实现了完全机器流程作业,没有人工干预,大大提高数据处理能力,可以实现压裂裂缝动态展布实时监测、指导储层压裂作业。
在步骤109,在压裂区域预先设定阈值,当能量梯度下降值大于该阈值时,表明该区域范围内为压裂裂缝的密集区,从而得到该时间段内的压裂裂缝展布图像。
在步骤110,沿采集时间段进行处理,对后续时间段的地面微地震宽频数据重复步骤104-109,计算不同时间段内压裂裂缝展布图像,从而达到监测压裂裂缝动态展布的发展情况。本发明的微地震监测方法采用了三分量微地震数据、分频段极化滤波方法、纵横波联合应用、偏移成像和梯度下降技术,根据不同分量之间的相关性,同时利用P波和S波分量进行微地震信号的识别与定位,通过偏移成像和梯度下降方法消除噪音的影响,从而可以极大地提高微地震事件的识别能力,提高储层压裂裂缝的监测精度,达到裂缝动态展布监测的目的。如图3所示,通过该专利的应用实例,可实时监测压裂改造中裂缝的发生、发展及展布过程,可实时评价压裂效果,提供指导建议。
在应用本发明的一具体实施例中,采用如下地面微地震数据:时间采样点数为N,时间采样间隔为0.5ms、1ms或2ms,可根据情况进行设置;台站仪器个数为M,三分量微地震数据为:其中:i=1,2,3,...,N;k=1,2,3,...,M;j=1,2,3,共计3个分量,即两个水平分量和一个垂直分量。该实施例包括了以下步骤:
在步骤S001时,在储层水力压裂区域的上方地表或浅井中,采用均匀、对称、充分、连续的方式,布设天然地震台站仪器或宽频检波器,共计M个台站,以获取致密储层水力压裂改造的x,y,z三分量微地震宽频数据
作为优化,在步骤S001,采用高分辨率、宽频带的三分量天然地震台站仪器或宽频检波器,在压裂井上方地表或浅井中,以地面微地震观测系统最优化设计为指导思路,采用均匀、对称、充分、连续的最佳观测方式布设。为了降低地表人为活动噪音的干扰,采取适当、有效的挖坑和埋置措施,尽可能降低地表噪音的干扰,一般需要挖一个约1m3的坑,在坑底用混凝土做一个基底,将天然地震台站仪器安置于基底之上,然后用泡沫包围上,上面再用土掩埋,减少温度变化和空气流动对天然地震台站仪器的干扰。所布设台站监测孔径要大于压裂段垂深的1.5倍以上。在均匀、对称、充分、连续的基础上,根据该地区地质构造、速度模型与震源机制的压裂段能量传播特征,在分散式布设台站仪器时,保证微地震监测能量在起伏地表的均匀分布,以获取高品质的连续微地震数据。
在步骤S002,根据压裂工区情况,采用预定大小的网格尺度,将储层水力压裂区域范围进行网格点剖分。网格大小的选择根据定位精度要求要预先确定,定位精度要求越高,网格划分的越小。并结合三维地面地震、VSP地震(井中垂直地震)、井间地震和测井资料,建立压裂区域的块状精细全三维速度模型。
在步骤S003,基于射线高斯束追踪方法,计算每一个剖分剖分网格点经过块状全三维速度模型后,到达地表台站的旅行时间、方位角度和出射角度。
作为优化,在步骤S003,需要压裂区域内的全三维高精度速度模型,因此需要结合三维地面地震、VSP地震(井中垂直地震)、井间地震和测井资料,并采用块状精细全三维速度模型方法,实现速度的逐渐变化,避免层状介质速度模型误差大的缺点。块状精细速度建模方法摒弃了传统层状介质的假设条件,以多边形立方体为基础,形成高精度速度模型,块状速度值具备了逐渐变化的特征,更加逼近于真实的地质构造和岩性特征。
由于剖分网格点信息的计算量巨大,可以在压裂前计算出所有剖分网格点到所有台站仪器的P波和S波旅行时间、方位角度和出射角度,建立一个剖分网格点信息数据库,在后续计算中直接读取这个数据库,可以大大降低后续计算时间,提高计算效率。
在步骤S004,采用高精度谱分解方法,将宽频微地震数据自动化分解为三个频段的微地震数据,三个频段的微地震数据分别是微地震事件可靠频段、微地震事件优势频段、微地震事件劣势频段,并将三个频段数据分别赋予不同的占比权重。在高精度谱分解之后,计算不同频段数据的三分量波形协方差矩阵,计算一段时间内的分频段数据协方差矩阵时间采样点总数为N,计算公式为:
其中:
L=(N-1)/2,L是时间窗口长度的一半;msk l(i)、msk m(i)为微地震三分量的分频段数据(l,m=1,2,3)。计算协方差矩阵的特征值(λ123)和相应的特征向量(V1,V2,V3)。利用(λ123)可以计算该时刻内三分量微地震记录的线偏振特征。
作为优化,步骤S004中采用全局偏振参数Mτ或线性偏振度MRL来表示,
其中:MQ决定MRL对不同偏振程度的敏感度,一般情况下:MQ≤1。
在步骤S005,根据线偏振特征对三个频段的微地震数据进行极化滤波处理,极化滤波可以表示为:
A(i)取线偏振特征Mτ或MRL,或者其它表示线偏振特征的量,为微地震三分量的分频段数据,其中:i=1,2,3,...,N,为微地震数据的离散采样;k=1,2,3,...,M,为微地震台站序号;j=1,2,3,共计3个分量,即两个水平分量和一个垂直分量。
然后利于谱分解原理将分频段数据合并,利用协方差矩阵的特征值计算频段微地震数据的线偏振特征,利用线偏振特征进行极化滤波处理,可大幅增强微地震事件能量,同时提高微地震数据的信噪比,有助于微地震事件的准确识别。如图2所示,图2为本发明的一具体实施例中三分量地面微地震宽频数据进行谱分解分频段极化滤波处理前、后的对比效果图,利用线偏振特征对分频段数据进行极化滤波处理后,微地震事件的能量大幅增强,微地震数据的信噪比得到提高,有助于微地震事件的准确识别和定位精度的提高。
在步骤S006,重复步骤S004和S005,对所有台站的微地震宽频数据进行相同处理,获取所有微地震台站仪器的分频段极化滤波后微地震宽频数据
在步骤S007,任选一个剖分网格点x,根据该剖分网格点到不同台站仪器的旅行时间、方位角度和出射角度,对所有台站仪器的微地震宽频数据进行分频段极化滤波处理,并对分频段极化滤波后的微地震宽频数据进行矢量偏移成像,获取该剖分网格点的偏移成像能量值。
作为优化,网格点x的能量值可以利用公式(6)表示:
其中:c=P or S,τ为不同台站所对应的旅行时间,L=(N-1)/2,L是时间窗口长度的一半,M为台站仪器个数,wc为权重值,为微地震数据的极化滤波结果,j=1,2,3,共计3个分量;令θ为出射角度,φ为方位角度,则射线矢量可记为:
Vr=(sinθ,cosθcosφ,cosθsinφ) (7)
权重值wc的计算公式为:
其中:Vr为射线矢量,V1为协方差矩阵的特征向量。
选取剖分网格点,将宽频微地震数据分解为三个频段的微地震数据,对所有台站仪器的微地震宽频数据进行分频段极化滤波处理,并对分频段极化滤波后的微地震宽频数据进行矢量偏移成像,既可以对P波进行偏移成像,也可以对S波进行偏移成像,还可以将P波和S波进行联合偏移成像。将分频段极化滤波与P/S波进行联合应用,一方面,大幅提高微地震事件识别准确度;另一方面,通过矢量偏移成像大幅提高微地震事件定位精度。
在步骤S008,重复S007,计算所有剖分网格点在某一时间段内的偏移成像能量值,并利用最优化方法计算该时间段内的能量梯度下降值。
作为优化,通过该时间段内的能量梯度下降值,实现微地震事件定位精度,通过梯度下降代替能量值,消除了异常振幅噪音的干扰,提高了微地震数据定位精度。
在步骤S009,在压裂区域预先设定阈值,当能量梯度下降值大于该阈值时,表明该区域范围内为压裂裂缝的密集区,从而得到该时间段内的压裂裂缝展布图像。由于不同微地震发生的位置很接近,因此微地震数据也有很强的相似性,大量微地震信号叠加后将会大大增强微地震信号,提高信噪比,能够识别更多的微地震信号。
作为优选,可以将阈值定义为Mthres=f×TMc,其中f为一系数,根据储层地质情况决定,一般可取1<f<10,TMc为该时刻模型空间所有格点x上的能量Uc(x)的梯度下降值。
在步骤S010,沿采集时间段进行处理,对后续时间段的地面微地震宽频数据重复S004-S009,计算不同时间段内压裂裂缝展布图像,从而达到监测压裂裂缝动态展布的发展情况。如图3所示,地面微地震压裂裂缝动态展布监测效果图,可实时监测压裂改造中压裂裂缝的发生、发展及展布过程,可实时评价压裂效果,提供指导建议。
如上所述,根据实施例的微地震监测方法采用了谱分解自动化分频段极化滤波技术和矢量偏移成像技术,从而可以极大地提高压裂裂缝展布实时监测精度,达到监测裂缝动态发展的目的。
根据实施例方法,整个处理流程实现了完全机器流程作业,没有人工干预,大大提高数据处理能力,可以实现压裂裂缝动态展布实时监测、指导储层压裂作业。避免了传统定位技术中微地震事件难以识别的不足,对于信噪比较低的地表微地震监测具有广阔的应用前景,满足工业化生产的需求。
在进行微地震信号定位后,得到一系列的离散微地震震源点,针对离散震源点,可以利用解释方法进行裂缝解释等,实现微地震裂缝的监测,微地震事件的定位和解释方法不能替代本方案,无法实现裂缝动态展布的监测。本发明实现了微地震处理和解释方法的紧密结合,简化了微地震监测流程,提高了现场工作效率,更有利于实现裂缝动态展布监测和指导储层压裂优化。
虽然本发明已经示出并描述了实施例的示例,但是本领域技术人员应该理解的是,实施例不限于此,不能认定本发明的具体实施只局限于这些说明,在不脱离本发明原理和构思的前提下,在不脱离如权利要求的精神和范围的情况下,可以进行各种改变和修改,都应视为属于本发明的保护范围。

Claims (16)

1.地面微地震压裂裂缝动态展布实时监测方法,其特征在于,该地面微地震压裂裂缝动态展布实时监测方法包括:
步骤1,采用观测方式获取连续微地震数据,将压裂区域网格剖分并建立块状全三维速度模型;
步骤2,结合速度模型各向异性及震源机制解,计算剖分网格点的旅行时间、方位角度和出射角度;
步骤3,采用谱分解方法进行微地震数据分频段自动化处理,计算三个频段数据的三分量波形协方差矩阵及线偏振特征;
步骤4,根据线偏振特征分别对分频段数据进行极化滤波,然后将分频段数据进行合并;
步骤5,任选剖分网格点,对极化滤波后的微地震数据进行矢量偏移成像,获取该网格点的成像能量值;
重复步骤5,计算所有剖分网格点在某一时间段内的偏移成像能量值,并计算该时间段内的能量梯度下降值;
步骤6,将能量梯度下降值大于一预先设定阈值的区域视为裂缝密集区,得到该时间段内的裂缝展布图像。
2.根据权利要求1所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤1中,在压裂区域上方地表或浅井中,采用分散式布设台站仪器,以获取连续微地震数据。
3.根据权利要求2所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤1中,在储层水力压裂区域的上方地表或浅井中,采用均匀、对称、充分、连续的方式,布设天然地震台站仪器或宽频检波器,以获取致密储层水力压裂改造的x,y,z三分量微地震宽频数据。
4.根据权利要求3所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,采用高分辨率、宽频带的三分量天然地震台站仪器或宽频检波器进行观测方式布设时:一方面,以地面微地震观测系统最优化设计为指导思路,采用均匀、对称、充分、连续的最佳观测方式布设,并采取挖坑和埋置措施,降低地表噪音的干扰,所布设台站监测孔径要大于压裂段垂深的1.5 倍以上;另一方面,在均匀、对称、充分、连续的基础上,根据该区域地质构造、速度模型与震源机制的压裂段能量传播特征,在分散式布设台站仪器时,保证微地震监测能量在起伏地表的均匀分布,以获取高品质的连续微地震数据。
5.根据权利要求1所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤1中,将储层水力压裂区域范围进行网格点剖分,并结合三维地面地震、井中垂直地震和测井资料,建立压裂区域的块状全三维速度模型,并采用块状全三维速度模型方法,实现全三维速度的逐渐变化,全三维速度模型包含:起伏真地表、复杂低速层、渐变降速带、平滑高速层顶界面、各项异性中深层速度模型。
6.根据权利要求1所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤2中,结合速度模型各向异性及震源机制解,基于射线高斯束追踪方法,采用速度逐渐变化的块状全三维速度模型,计算每一个剖分网格点经过块状全三维速度模型后,到达地表或浅井中台站的旅行时间、方位角度和出射角度。
7.根据权利要求1所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤3中,选取一个台站的三分量微地震宽频数据,采用高精度谱分解方法,将宽频微地震数据自动化分解为三个频段的微地震数据,分别计算三个频段数据的三分量波形协方差矩阵,分别利用协方差矩阵的特征值计算三个频段微地震数据的线偏振特征。
8.根据权利要求7所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤3中,采用高精度谱分解方法,将宽频微地震数据分解为三个频段的微地震数据,三个频段的微地震数据分别是微地震事件可靠频段、微地震事件优势频段、微地震事件劣势频段,在高精度谱分解之后,计算不同频段数据的三分量波形协方差矩阵,利用协方差矩阵的特征值计算频段微地震数据的线偏振特征,利用线偏振特征进行极化滤波处理,增强微地震事件,提高微地震数据的信噪比。
9.根据权利要求7所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤3中,在高精度谱分解之后,计算不同频段数据的三分量波形协方差矩阵,计算一段时间内的分频段数据的协方差矩阵时间采样点总数为N,计算公式为:
其中:
L=(N-1)/2,L是时间窗口长度的一半;为微地震三分量的分频段数据,其中参数k为微地震台站序号,k=1,2,3,...,M,M为微地震台站总个数,l,m=1,2,3;计算协方差矩阵的特征值(λ123)和相应的特征向量(V1,V2,V3);利用(λ123)计算该时间内三分量微地震记录的线偏振特征;
采用全局偏振参数Mτ或线性偏振度MRL来表示,
其中:MQ决定MRL对不同偏振程度的敏感度,MQ≤1。
10.根据权利要求1所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤4中,根据线偏振特征对三个频段,即微地震事件可靠频段、微地震事件优势频段、微地震事件劣势频段的微地震数据进行极化滤波处理,然后利用谱分解原理将分频段数据合并。
11.根据权利要求10所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤4中,根据线偏振特征对三个频段的微地震数据进行极化滤波处理,极化滤波表示为:
A(i)取线偏振特征Mτ或MRL,其中Mτ为全局偏振参数,MRL为线性偏振度,或者其它表示线偏振特征的量,为微地震三分量的分频段数据,其中:i=1,2,3,...,N,为微地震数据的离散采样;k=1,2,3,...,M,为微地震台站序号;j=1,2,3,共计3个分量,即两个水平分量和一个垂直分量。
12.根据权利要求1所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,该地面微地震压裂裂缝动态展布实时监测方法还包括,在步骤4之后,重复步骤3和4,对所有台站的微地震宽频数据进行相同处理,获取所有微地震台站仪器的分频段极化滤波后微地震宽频数据。
13.根据权利要求12所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤5中,选取一个剖分网格点,根据该剖分网格点到不同台站仪器的旅行时间、方位角度和出射角度,对所有台站仪器的微地震宽频数据进行分频段极化滤波处理,并对分频段极化滤波后的微地震宽频数据进行矢量偏移成像,获取该剖分网格点的偏移成像能量值。
14.根据权利要求13所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤5中,选取剖分网格点,将宽频微地震数据分解为三个频段的微地震数据,对所有台站仪器的微地震宽频数据进行分频段极化滤波处理,并对分频段极化滤波后的微地震宽频数据进行矢量偏移成像,对P波进行偏移成像,或对S波进行偏移成像,或将P波和S波进行联合偏移成像。
15.根据权利要求1所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,在步骤6中,在压裂区域预先设定阈值,当能量梯度下降值大于该阈值时,表明该区域范围内为压裂裂缝的密集区,从而得到该时间段内的压裂裂缝展布图像。
16.根据权利要求1所述的地面微地震压裂裂缝动态展布实时监测方法,其特征在于,该地面微地震压裂裂缝动态展布实时监测方法还包括,在步骤6之后,沿采集时间段进行处理,对后续时间段的地面微地震宽频数据重复步骤3-6,计算不同时间段内压裂裂缝展布图像,监测压裂裂缝动态展布的发展情况。
CN201710915453.0A 2017-09-29 2017-09-29 地面微地震压裂裂缝动态展布实时监测方法 Active CN107728200B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710915453.0A CN107728200B (zh) 2017-09-29 2017-09-29 地面微地震压裂裂缝动态展布实时监测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710915453.0A CN107728200B (zh) 2017-09-29 2017-09-29 地面微地震压裂裂缝动态展布实时监测方法

Publications (2)

Publication Number Publication Date
CN107728200A CN107728200A (zh) 2018-02-23
CN107728200B true CN107728200B (zh) 2019-03-29

Family

ID=61209477

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710915453.0A Active CN107728200B (zh) 2017-09-29 2017-09-29 地面微地震压裂裂缝动态展布实时监测方法

Country Status (1)

Country Link
CN (1) CN107728200B (zh)

Families Citing this family (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108680950B (zh) * 2018-05-16 2019-07-26 吉林大学 一种基于自适应块匹配的沙漠地震信号位置检测方法
CN109100787A (zh) * 2018-06-29 2018-12-28 安徽万泰地球物理技术有限公司 一种等离子体脉冲谐振增透的微震监测评价方法
US11353612B2 (en) * 2019-03-11 2022-06-07 Saudi Arabian Oil Company Nonstationary maximum likelihood method to estimate dispersion spectra for full waveform sonic logging
CN110208857B (zh) * 2019-05-28 2021-01-05 哈尔滨工程大学 一种海底沉积层区域地震计信号振幅测量方法
CN110727028A (zh) * 2019-09-17 2020-01-24 河南理工大学 一种基于地面微地震监测的煤储层裂隙评价方法
CN110568487A (zh) * 2019-09-19 2019-12-13 中国科学技术大学 基于天然地震波形的活跃断层结构成像方法
CN112558145B (zh) * 2019-09-25 2024-07-09 中国石油化工股份有限公司 基于图像处理的微地震有效事件识别方法及系统
CN112647935B (zh) * 2019-10-12 2024-06-18 中国石油化工股份有限公司 压裂裂缝参数计算方法及系统
CN110927785A (zh) * 2019-12-11 2020-03-27 山东省煤田地质规划勘察研究院 一种矿区水力压裂裂缝微地震监测应用方法
CN113009578B (zh) * 2019-12-19 2022-11-08 新奥科技发展有限公司 生产井的井身轨迹获取方法及生产井的连通方法
CN111025392B (zh) * 2019-12-27 2020-07-24 中国矿业大学 一种利用微震信号的煤岩体压裂裂缝实时快速监测评价方法
CN111175816B (zh) * 2020-01-06 2022-04-15 中国石油化工股份有限公司 油藏改造实时构建微地震裂缝网络的方法及装置
CN111175815B (zh) * 2020-01-06 2022-04-15 中国石油化工股份有限公司 油藏改造微地震监测裂缝震源机制求解方法及系统
CN113050159B (zh) * 2021-03-23 2021-11-16 中国矿业大学 一种煤岩水力压裂裂缝微震定位及扩展机理监测方法
CN113253343B (zh) * 2021-05-12 2022-05-31 中油奥博(成都)科技有限公司 基于微地震监测技术识别地下储气库断层活动的方法
CN113945966B (zh) * 2021-05-25 2022-05-03 中国矿业大学(北京) 人工压裂裂缝网络构建方法及装置

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6473696B1 (en) * 2001-03-13 2002-10-29 Conoco Inc. Method and process for prediction of subsurface fluid and rock pressures in the earth
CN102193108B (zh) * 2010-03-19 2013-07-03 中国石油天然气集团公司 一种提高石油勘探资料处理信噪比的方法
US9784866B2 (en) * 2013-07-28 2017-10-10 Geokinetics Usa, Inc. Method and apparatus for enhanced monitoring of induced seismicity and vibration using linear low frequency and rotational sensors
CN105445801B (zh) * 2014-09-01 2018-06-01 中国石油化工股份有限公司 一种消除二维地震资料随机噪音的处理方法
CN106154334B (zh) * 2015-04-13 2018-02-16 中石化石油工程地球物理有限公司胜利分公司 基于网格搜索的井下微地震事件实时反演定位方法
CN105044769B (zh) * 2015-06-10 2017-11-10 中国石油集团川庆钻探工程有限公司地球物理勘探公司 提高地震信号的分辨率的方法
CN105954795A (zh) * 2016-04-25 2016-09-21 吉林大学 一种用于微地震定位的网格逐次剖分方法
RU2620785C1 (ru) * 2016-06-14 2017-05-29 Общество с ограниченной ответственностью "Макросейс" Способ определения местоположения очага микросейсмического события

Also Published As

Publication number Publication date
CN107728200A (zh) 2018-02-23

Similar Documents

Publication Publication Date Title
CN107728200B (zh) 地面微地震压裂裂缝动态展布实时监测方法
Parker et al. Active‐source seismic tomography at the Brady geothermal field, Nevada, with dense nodal and fiber‐optic seismic arrays
CN107526101B (zh) 一种获取地震反射波的采集和处理方法
CN108037526B (zh) 基于全波波场vsp/rvsp地震资料的逆时偏移方法
CN106094029B (zh) 利用偏移距矢量片地震数据预测储层裂缝的方法
CN104280775B (zh) 一种基于全波形矢量偏移叠加的微地震监测定位方法
Wang et al. Current developments on micro-seismic data processing
CN107490808B (zh) 一种高可靠性地震勘探观测系统的建立方法
CN102645670B (zh) 一种基于叠加响应分析的观测系统优化设计方法
CN110133715A (zh) 一种基于初至时差和波形叠加的微地震震源定位方法
CN103336297B (zh) 微破裂向量扫描方法
EP2818897B1 (en) Method and apparatus for estimating a parameter of a subsurface volume
CN109884710B (zh) 针对激发井深设计的微测井层析成像方法
CN108089226B (zh) 一种基于道间能量叠加的微地震事件自动识别方法
CN107728204A (zh) 基于叠前纵波各向异性的裂缝预测方法及系统
CN102053263A (zh) 调查表层结构的方法
CN110007340A (zh) 基于角度域直接包络反演的盐丘速度密度估计方法
CN107728205B (zh) 一种地层压力预测方法
CN109633745A (zh) 一种三维构造图的制图方法及装置
CA2961168A1 (en) Integrating vertical seismic profile data for microseismic anisotropy velocity analysis
Zeng et al. High‐resolution shallow structure at brady hot springs using ambient noise tomography (ANT) on a trenched distributed acoustic sensing (DAS) array
CN109469477B (zh) 一种人工裂缝延伸方向的预测方法和装置
CN112230274B (zh) 面向随钻导向的声波方程频率域逆时偏移快速成像方法
CN104155690B (zh) 基于椭球展开的三维地震数据叠加速度求取方法
CN106483570A (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
GR01 Patent grant
GR01 Patent grant