CN110780346A - 一种隧道超前探测复杂地震波场的分离方法 - Google Patents
一种隧道超前探测复杂地震波场的分离方法 Download PDFInfo
- Publication number
- CN110780346A CN110780346A CN201911144423.XA CN201911144423A CN110780346A CN 110780346 A CN110780346 A CN 110780346A CN 201911144423 A CN201911144423 A CN 201911144423A CN 110780346 A CN110780346 A CN 110780346A
- Authority
- CN
- China
- Prior art keywords
- wave
- time
- velocity
- frequency
- waves
- 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
- 238000000926 separation method Methods 0.000 title claims abstract description 31
- 238000001514 detection method Methods 0.000 title claims abstract description 28
- 230000010287 polarization Effects 0.000 claims abstract description 31
- 238000000034 method Methods 0.000 claims description 33
- 238000004458 analytical method Methods 0.000 claims description 28
- 238000001914 filtration Methods 0.000 claims description 21
- 238000001228 spectrum Methods 0.000 claims description 20
- 230000009466 transformation Effects 0.000 claims description 16
- 230000006870 function Effects 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 4
- 230000008878 coupling Effects 0.000 claims description 3
- 238000010168 coupling process Methods 0.000 claims description 3
- 238000005859 coupling reaction Methods 0.000 claims description 3
- 230000007547 defect Effects 0.000 claims description 3
- 230000008571 general function Effects 0.000 claims description 3
- 230000001902 propagating effect Effects 0.000 claims description 3
- 239000004576 sand Substances 0.000 claims description 3
- 238000012795 verification Methods 0.000 claims description 3
- 238000004519 manufacturing process Methods 0.000 claims description 2
- 230000001629 suppression Effects 0.000 claims description 2
- 230000001052 transient effect Effects 0.000 claims description 2
- 238000010276 construction Methods 0.000 abstract description 5
- 230000002194 synthesizing effect Effects 0.000 abstract description 3
- 238000010586 diagram Methods 0.000 description 5
- 238000013461 design Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000005284 excitation Effects 0.000 description 3
- 239000002245 particle Substances 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 229910052704 radon Inorganic materials 0.000 description 2
- SYUHGPGVQRZVTB-UHFFFAOYSA-N radon atom Chemical compound [Rn] SYUHGPGVQRZVTB-UHFFFAOYSA-N 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 210000005056 cell body Anatomy 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
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/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种隧道超前探测复杂地震波场的分离方法，包括获得炮点地震道集数据；求得瞬时方位角，根据方位角的特征，识别直达波与绕射横波；求得纵波速度与横波速度；根据纵波速度进行倾角的估算，综合利用偏移距与介质的速度，计算转换横波。本发明与隧道超前探测波场分离的现有技术相比，优点有：1、考虑了隧道施工环境的复杂地质条件和特殊的地震观测系统对地震波场有巨大影响，提高了隧道超前探测波场分离的适用性；2、综合地震波的时频特征、极化特征与视速度特征等地震波场多种特征识别有效波，提高了纵横波分离精度。
Description
技术领域
本发明属于地质物探技术领域，具体涉及利用S变换、F-K、极化分析方法综合实现地震波场分离的一种隧道超前探测复杂地震波场的分离方法。
背景技术
随着我国交通、水利水电事业的快速发展，大规模的铁路、公路和水电工程建设数量不断增加，需要修建大量各种类型的隧道，在隧道施工建设中，及时查明不良地质体变得尤为重要。
随着近几年多分量地震勘探技术的发展，波场分离技术也变得多样化、复杂化。常用的有极化滤波、Radon变换、奇异值分解、F-K滤波、中值滤波以及基于散度和旋度的滤波等。其中在极化滤波方法的发展里，利用单一的时间域或频率域的方法来识别有效反射波十分困难，一些基于时频分析的极化分析方法应运而生。Soma在2002年将小波变换与极化分析结合应用到多分量地震中，在时频域中构建多分量的协方差矩阵；Kulesh、Diallo等在之后几年相继提出通过时频谱的实部和虚部确定瞬时极化参数，再设计特定的滤波器进行转换波分离；Pinnegar在2006年提出利用S变换提取三分量地震信号极化属性的方法，该方法可以直接计算每一个采样点处的极化属性。
目前，关于纵横波波场分离的研究有一定的基础，但针对隧道超前探测的有效波波场分离的研究很少，极化滤波、Radon变换作为隧道超前探测的波场分离中最常用到方法，主要是根据极化特征与视速度特征进行波的分离，但具有一定的局限性。本发明基于几种常用的波场分离方法，并结合隧道超前探测的特殊观测系统，以及隧道的全空间波场特征，发展出了一种提高纵横波分离精度与适应度的新技术。
发明内容
根据上述阐述，本发明的目的在于提供一种隧道超前探测复杂地震波场的分离方法，是针对纵横波分离方法，对观测系统、施工条件适应性较强，且精度较高。
本发明提供的技术方案：
一种隧道超前探测复杂地震波场的分离方法，包括如下步骤：
S1、在隧道中，设计线性地震观测系统，获得至少一炮的炮点地震道集数据；
S2、读入地震记录并作S变换，求得瞬时方位角，根据方位角的特征，识别直达波与绕射横波；
S3、在单道时频记录中根据识别的时间逐一进行波形拾取，根据速度分析的原理，分别求得纵波速度与横波速度；
S4、对多道时频数据作傅里叶变换，得到时间、频率、波数域数据，通过改变时间节点，作频率波数域二维视速度滤波，得到压制了正视速度波的记录；
S5、在滤波后的数据中，当波形明显时，根据波的时频与速度特征，对反射纵波与转换横波直接拾取，则可直接完成隧道超前探测的波场分离；否则，根据波的极化特征与时频特征，初步拾取反射纵波，并根据纵波速度进行倾角的估算，综合利用偏移距与介质的速度，计算转换横波。
上述技术方案中，所述S5中，计算转换横波的旅行时，在地震记录中拾取转换横波，若该时间区间内未存在可分辨的特征或波形，则反射纵波的拾取存在问题，再次拾取并重复S2—S4步骤直至完成波场分离。
上述技术方案中，所述S2中，瞬时方位角求取为：
地震数据经傅里叶变换后，时频域数据是瞬时数据，将经S变换后的三分量地震数据进行极化计算，则瞬时极化主轴R(t,f)与瞬时极化次轴r(t,f)表示为：
瞬时极化主轴：
瞬时极化次轴：
其中，定义一个泛函数
表示取复数信号z的实部；u(t)是三分量地震数据；
为关于时频的实函数，当时，泛函数Φ(z)取得极值；
则瞬时方位角可表示如下：
识别直达波与绕射横波依据：随着波场的传播，有效波的方位角特征因受多个因素干扰而不确定，直达波和在隧道掌子面处激发的绕射波在传播过程中方位角是保持不变的，作为识别隧道周围直达波与绕射横波波场的依据。
上述技术方案中，所述S3中，采用叠加速度谱的方式进行速度扫描，通过设置规定的速度与时间间隔进行扫描，从记录中准确地拾取反射信号，得到正常时差，求取速度参数。
速度分析：隧道周围的波场，存在纵波与横波两种典型波形，在求取纵波速度时，可直接提取直达波进行速度分析，横波的波形通过直达波在掌子面处激发的绕射横波进行速度分析，直达波在地震波Y分量、绕射波横波在X分量上能量突出，
直达波的速度与时间存在以下关系：
绕射横波的速度与时间存在以下关系：
其中，震源到掌子面的垂直距离为L，偏移距为x，道间距为Δx，直达波传播到隧道掌子面R处的时间记为T0。
当vp、vs确定，将对应tp、ts时刻振幅值进行叠加，完成所有道的计算后，即可得到速度谱估计速度。
上述技术方案中，所述规定的速度是采用叠加能量或相似系数准则之一，作为速度分析的准则制作速度谱，选择的速度使得校正后时距曲线刚好拉平，叠后能量最强；或者相似系数准则中的使得各道反射信息最为相似。
上述技术方案中，所述S4中，将傅里叶变换后的时间、频率、波数域数据固定为时间与频率，通过对空间各处波动情况进行傅里叶变换，得到时间频率波数与数据，同时观察时间、频率、波数特征；查看任意时刻的频率波数域数据，即使某时刻存在波形耦合，其频率、视速度方向也可能不同，选择能够克服F-K滤波器易造成假频缺点的滤波器，在每个时间点，进行二维视速度滤波，随着时间的增长，某时刻的频率波数域数据受直达波影响变小。
上述技术方案中，所述S5中，拾取反射纵波的主要步骤如下：
Q1、设计半时间窗口长度Δt、半频率窗口长度Δf，选择初始道号i及反射纵波大致出现的时间ti与频率fi；
Q2、在第i道时频数据中，于时间[ti-Δt,ti+Δt]，频率[fi-Δf,fi+Δf]区域内搜索极值点，将其对应的时间Ti作为反射纵波在第i道中对应的时间；
Q3、在第i±1道时频数据中，于时间[Ti-Δt,Ti+Δt]，频率[fi-Δf,fi+Δf]区域内搜索极值点，将其对应的时间Ti±1作为反射纵波在第i±1道中对应的时间；
Q4、重复以上步骤，直至搜索完所有道数据，得到反射纵波在地震记录中的时间节点
Q5、设计拾取半波形的半时间窗口长度ΔT，在滤波后的地震记录中，设计滤波器拾取时间窗口[Ti-ΔT,Ti+ΔT]内的信号，得到反射纵波；
在滤波前直接拾取直达波，绕射波在滤波后拾取，拾取方法均采用以上步骤；
上述技术方案中，所述S5中，所述半时间窗口长度Δt应在至少包含半波长范围内，选择最小的，避免在波形交叉时，拾取干扰波；所述波形的半时间窗口长度ΔT，应至少大于一个波长范围。
上述技术方案中，所述Q5中，反射纵波经过压制正视速度波后，根据方位角特征与视速度特征拾取反射纵波，针对拾取的反射纵波的波形，计算倾角谱，具体实现步骤如下：
Q5.1、拾取Y分量的反射纵波，得到含有单一波形的记录；
Q5.2、设置初始角度α0、扫描倾角间隔Δα和倾角最大值αm；
Q5.3、设置扫描时间窗口长度T；
Q5.4、当界面倾角为α时，反射波旅行时为t2，则其关系表示如下：
根据反射波旅行时t2与速度分析得到的纵波速度vp，计算每一个扫描角度αi对应的波的传播时间ti，将地震记录中该时刻对应的振幅值进行叠加，记录在Eα(ti,αi)中；
Q5.5、扫描过全部倾角后，Eα即为扫描的界面倾角谱，据此得到倾角最佳值，以及对应的自激自收时间t0。
上述技术方案中，根据倾角分析估算地层倾角后，计算同一反射界面对应激发的转换横波的时间点，若计算出的时间区域内存在同相轴，且该区域内的时频特征符合转换横波的特点，则验证反射纵波的拾取正确性，反之，如果该时间区域内未存在明显同相轴，应当重新根据波形特征进行拾取，重复之前的步骤进行验证。
本发明与隧道超前探测波场分离的现有技术相比，优点有：1、考虑了隧道施工环境的复杂地质条件和特殊的地震观测系统对地震波场有巨大影响，提高了隧道超前探测波场分离的适用性；2、综合地震波的时频特征、极化特征与视速度特征等地震波场多种特征识别有效波，提高了纵横波分离精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1是本发明优选实施例的流程图；
图2为质点振动示意图；
图3为绕射横波传播路径示意图；
图4为倾斜断层模型的反射波传播路径示意图；
图5为模型的三维示意图；
图6为模型的观测系统布置图；
图7为模拟得到的部分地震合成记录图；
图8为不同偏移距的方位角绝对值示意图；
图9(a)和(b)分别为模拟数据的纵波和横波速度谱图；
图10为模拟地震数据的不同时刻的频率波数域数据图；
图11为模拟地震数据第一道时频数据滤波前(a)、滤波后(b)对比图；
图12为模拟地震记录中的反射纵波图；
图13为模型界面倾角的扫描结果图；
图14(a)、(b)和(c)分别是模拟地震数据的X、Y、Z分量的分离结果图。
具体实施方式
下面将结合本发明的附图，对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
根据图1所示本发明的基于S变换、F-K、极化分析方法综合实现隧道超前探测复杂地震波场分离的优选实施例的流程图。
本发明是综合地震波的时频特征、极化特征与视速度特征等地震波场多种特征进行反射波分离的方法，主要包括如下步骤：
S1、在隧道中，设计线性地震观测系统，采集获得一炮或多炮共炮点地震道集数据；
S2、读入地震记录并作S变换，求得瞬时方位角，根据方位角的特征识别直达波与绕射横波；
S变换：在S变换中，基本小波是由高斯函数与简谐波的乘积构成的。S变换的时窗函数可表示如下：
设h(t)为输入信号，则S变换可表示为：
令
u(t,f)＝h(t)e-2iπft
则S变换可以写作：
S(t,f)＝F-1[F[u(t,f)]·F[v(t,f)]]，F表示傅里叶变换，
若由时频域结果经S反变换得到原始信号h(t)，则首先作关于时间t的压缩求和，之后再进行傅里叶反变换。
S变换与S反变换的离散形式可表示为：
S正变换
S反变换
求取瞬时方位角：当地震数据经S变换后，其时频域数据已是瞬时数据，将经S变换后的三分量地震数据进行极化计算，则瞬时极化主轴R(t,f)与瞬时极化次轴r(t,f)可表示为：
其中，定义一个泛函数 
表示取复数信号z的实部；
u(t)是三分量地震数据；为关于时频的实函数，当
则瞬时方位角可表示如下：
根据图2所示的质点振动示意图，以纵波为例，r方向为射线方向，也为质点振动方向，其坐标为(x,y,z)，则方位角γ＝arctan(y/x)。
识别直达波与绕射横波：随着波场的传播，虽然有效波的方位角特征不定，受多个因素干扰，但将炮点置于测线上时，直达波在传播过程中方位角是保持不变的，且在隧道掌子面处激发的绕射波，其方位角也是近似不变的，这可作为识别隧道周围波场的依据。
S3、在单道时频记录中根据识别的时间逐一进行波形拾取，根据速度分析的原理，分别求得纵波速度与横波速度；
本发明采用叠加速度谱的方式进行速度扫描，通过设置一定的速度与时间间隔进行扫描，以叠加能量或相似系数等作为速度分析的准则制作速度谱。从记录中准确地拾取反射信号，得到正常时差，求取速度参数。
速度分析：隧道周围的波场存在纵波与横波两种典型波形，在求取纵波速度时，可直接提取直达波进行速度分析，而横波的波形则可通过直达波在掌子面处激发的绕射横波进行速度分析，且直达波在地震波Y分量、绕射波横波在X分量上能量突出。
直达波的速度与时间存在以下关系：
绕射横波的速度与时间存在以下关系：
根据图3所示的绕射横波传播路径，检波点G点为接收点，震源到掌子面的垂直距离为L，偏移距为x，道间距为Δx，直达波传播到隧道掌子面R处的时间记为T0。
当vp、vs确定，将对应tp、ts时刻振幅值进行叠加，完成所有道的计算后即可得到速度谱估计速度。
具体实现步骤如下：
①拾取X分量的绕射横波与Y分量的直达波，得到含有单一波形的记录；
②设置初始速度v0、扫描速度间隔Δv和速度最大值vm；
③设置扫描时间窗口长度T；
④当界面倾角为α时，反射波旅行时为t2，则其关系表示如下：
根据直达波和绕射波的tp、ts计算公式，计算每一个扫描速度vi对应的波的传播时间ti，将地震记录中ti时刻对应的振幅值进行叠加，记录在Ev(ti,vi)中；
⑤扫描过全部速度后，Ev即为对应波形的速度谱，据此得到波形的最佳速度。
S4、对多道时频数据作傅里叶变换，得到时间、频率、波数域数据，通过改变时间节点，作频率波数域二维视速度滤波，得到压制了正视速度波的记录；
进一步的，具体实施过程：将S变换后的时频域数据固定时间与频率，通过对空间各处波动情况进行傅里叶变换，得到时间频率波数与数据，可以同时观察时间、频率、波数特征。另外，查看任意时刻的频率波数域数据，即使某时刻存在波形耦合，其频率、视速度方向也可能不同，设计合适的滤波器，在每个时间点进行二维视速度滤波，在随着时间的增长，某时刻的频率波数域数据受直达波影响变小，克服直接使用f-k滤波的缺点。
S5、在滤波后的数据中，当波形明显时，可根据波的时频与速度特征对反射纵波与转换横波直接拾取，则可直接完成隧道超前探测的波场分离；否则，根据波的极化特征与时频特征初步拾取反射纵波，并根据纵波速度进行倾角的估算，综合利用偏移距与介质的速度，计算转换横波的旅行时，在地震记录中拾取转换横波，若该时间区间内未存在可分辨的特征或波形，则反射纵波的拾取存在问题，可再次拾取并重复S2～S4步骤直至完成波场分离。
(1)拾取反射纵波的主要步骤如下：
①设计半时间窗口长度Δt、半频率窗口长度Δf，选择初始道号i及反射纵波大致出现的时间ti与频率fi；
②在第i道时频数据中，于时间[ti-Δt,ti+Δt]，频率[fi-Δf,fi+Δf]区域内搜索极值点，将其对应的时间Ti作为反射纵波在第i道中对应的时间；
③在第i±1道时频数据中，于时间[Ti-Δt,Ti+Δt]，频率[fi-Δf,fi+Δf]区域内搜索极值点，将其对应的时间Ti±1作为反射纵波在第i±1道中对应的时间；
④重复以上步骤，直至搜索完所有道数据，得到反射纵波在地震记录中的时间节点；
⑤设计拾取波形的半时间窗口长度ΔT，在滤波后的地震记录中，设计滤波器拾取时间窗口[Ti-ΔT,Ti+ΔT]内的信号，得到反射纵波。
在拾取直达波时，可在滤波前直接拾取，而绕射波可在滤波后拾取，拾取方法均采用以上步骤。需要指出的是，取的半时间窗口长度Δt与ΔT是不同的，Δt应在至少包含半波长范围内尽量小，避免波形交叉时，拾取到干扰波，而ΔT应至少大于一个波长范围。
由于转换横波的能量相比反射纵波的能量很弱，在提取时存在一定的问题，但反射纵波经过压制正视速度波后根据方位角特征与视速度特征可拾取，针对拾取的波形可计算倾角谱，具体实现步骤如下：
A、拾取Y分量的反射纵波，得到含有单一波形的记录；
B、设置初始角度α0、扫描倾角间隔Δα和倾角最大值αm；
C、设置扫描时间窗口长度T；
D、当界面倾角为α时，反射波旅行时为t2，则其关系可表示如下：
根据图4所示，根据反射波旅行时t2与速度分析得到的纵波速度vp，计算每一个扫描角度αi对应的波的传播时间ti，将地震记录中该时刻对应的振幅值进行叠加，记录在Eα(ti,αi)中；
E、扫描过全部倾角后，Eα即为扫描的界面倾角谱，据此得到倾角最佳值，以及对应的自激自收时间t0。
根据以上倾角分析估算地层倾角后，计算同一反射界面对应激发的转换横波的时间点，若计算出的时间区域内存在同相轴，且该区域内的时频特征符合转换横波的特点，则可验证反射纵波的拾取正确性，反之，如果该时间区域内未存在明显同相轴，则存在两种可能性，一种是反射界面与测线垂直，使得波场传播到界面时，并未激发转换横波，另一种则是反射纵波拾取有误，应当重新根据波形特征进行拾取，重复之前的步骤进行验证。
为了进一步清晰地理解本发明，我们选择一个含有断层与隧道的模型进行详细说明。
设计的模型空间大小为100m×150m×50m，即x轴方向100m，y轴方向150m，z轴方向60m，隧道大小为10m×60m×10m，隧道中心位于(50，30，25)处，网格间距Δx＝Δy＝Δz＝0.5m，为含倾角与方位角断层模型，根据图5所示，模型具体物理参数见表1；我们采用交错网格有限差分正演模拟获得合成地震记录，采用线性观测系统，统一采用单炮激发、多道接收的方式，利用三分量检波器同时接收，即炮点和检波器在隧道壁内，埋深2米，炮点位于模型的(57,40,25)处，测线位于x＝57m,z＝25m处，具体设计根据图6；震源为雷克子波，主频为600Hz，道间距为0.5m，时间采样率为0.05ms，采样时长为100ms。
表1模型的物性参数
选择在x＝57m,y＝28～57.5m,z＝25m处接收到的地震记录进行试验，地震记录根据图7所示。
更改不同的道号，观察方位角的变换趋势，取方位角的绝对值根据图8所示，可大致分析波形特征。
根据方位角特征，拾取直达波与绕射波信息，进行速度分析。设置速度间隔为2m/s，得到的速度谱根据图9所示，纵波速度为3600m/s，横波速度为2326m/s。
根据图10所示，在单道的时频记录中，虽可以根据时间与频率的差异识别波形，但无法区分不同波形的视速度差异，故在每个时频点处作空间上的傅里叶变换，得到三维时间频率波数域数据，
通过不同的时间切片可以大致估计波形存在的时长，若波形明显，则可以推测其在地震记录中首次出现的时间，该时间为此波形在最远偏移距，即最靠近掌子面处检波器接收的记录。
图10列举典型时间点对应的频率波数图，由图可以看出，当t＝15ms时，于X分量可分辨明显具有负视速度的波形，Y分量可见能量较弱的负视速度波，而Z分量未见明显波形。根据表现出的时间和视速度特征与在三分量检波器接收到的能量分配上，X分量为绕射横波，Y分量为绕射纵波。绕射波持续出现一段时间后，逐渐消失，当t＝22.5～29.5ms时，Y分量显示出明显的负视速度波，且根据视速度大小推测，该波形为纵波，且Z分量可见微弱波形，在t＝32.5ms时，Y分量的能量逐渐减弱，X分量存在某一负视速度波，该波波速近似为横波速度，Z分量未见明显波形。
在三维时间频率波数域的每个时间点处设计合适的滤波器压制正视速度波，通过傅里叶反变换到时频域后可得到滤波后的多道时频记录，以第一道数据为例，滤波前后对比根据图11所示。在压制了正视速度波的记录中可见明显反射波波形，拾取反射纵波根据图12所示。
基于反射波的拾取结果，对前方异常体的倾角进行估算，在条件允许的情况下，可尽可能缩小扫描间隔与扫描区间，将结果表示地更加精确。本次试验设置倾角扫描间隔为0.05°，扫描区间为[0,90°]，倾角扫描结果根据图13所示，倾角大小为69.75°。
在得到倾角估计结果后，可通过时距曲线公式计算反射纵波与转换横波的时间节点，并设计合适滤波器进行拾取。得到的分离结果根据图14所示。
由图14可以看出，由于模型的界面设置了一定的方位角与倾角，有效波在三个分量上均有明显同相轴，且反射纵波在Y分量相对明显，转换横波在X分量更加明显，经过与未含有隧道模型的地震记录对比，拾取的有效波波形与极化特征等均符合要求，故可验证拾取的正确性。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。
Claims (10)
1.一种隧道超前探测复杂地震波场的分离方法，其特征在于：包括如下步骤：
S1、在隧道中，设计线性地震观测系统，获得至少一炮的炮点地震道集数据；
S2、读入地震记录并作S变换，求得瞬时方位角，根据方位角的特征，识别直达波与绕射横波；
S3、在单道时频记录中根据识别的时间逐一进行波形拾取，根据速度分析的原理，分别求得纵波速度与横波速度；
S4、对多道时频数据作傅里叶变换，得到时间、频率、波数域数据，通过改变时间节点，作频率波数域二维视速度滤波，得到压制了正视速度波的记录；
S5、在滤波后的数据中，当波形明显时，根据波的时频与速度特征，对反射纵波与转换横波直接拾取，则可直接完成隧道超前探测的波场分离；否则，根据波的极化特征与时频特征，初步拾取反射纵波，并根据纵波速度进行倾角的估算，综合利用偏移距与介质的速度，计算转换横波。
2.根据权利要求1所述的一种隧道超前探测复杂地震波场的分离方法，其特征在于：所述S5中，计算转换横波的旅行时，在地震记录中拾取转换横波，若该时间区间内未存在可分辨的特征或波形，则反射纵波的拾取存在问题，再次拾取并重复S2—S4步骤直至完成波场分离。
3.根据权利要求1所述的一种隧道超前探测复杂地震波场的分离方法，其特征在于：所述S2中，瞬时方位角求取为：
地震数据经傅里叶变换后，时频域数据是瞬时数据，将经S变换后的三分量地震数据进行极化计算，则瞬时极化主轴R(t,f)与瞬时极化次轴r(t,f)表示为：
瞬时极化主轴：
瞬时极化次轴：
其中，定义一个泛函数
表示取复数信号z的实部；u(t)是三分量地震数据；
为关于时频的实函数，当
时，泛函数Φ(z)取得极值；
则瞬时方位角可表示如下：
识别直达波与绕射横波依据：随着波场的传播，有效波的方位角特征因受多个因素干扰而不确定，直达波和在隧道掌子面处激发的绕射波在传播过程中方位角是保持不变的，作为识别隧道周围直达波与绕射横波波场的依据。
4.根据权利要求1所述的一种隧道超前探测复杂地震波场的分离方法，其特征在于：所述S3中，采用叠加速度谱的方式进行速度扫描，通过设置规定的速度与时间间隔进行扫描，作为速度分析的准则制作速度谱，从记录中准确地拾取反射信号，得到正常时差，求取速度参数；
速度分析：隧道周围的波场，存在纵波与横波两种典型波形，在求取纵波速度时，可直接提取直达波进行速度分析，横波的波形通过直达波在掌子面处激发的绕射横波进行速度分析，直达波在地震波Y分量、绕射波横波在X分量上能量突出，
直达波的速度与时间存在以下关系：
绕射横波的速度与时间存在以下关系：
其中，震源到掌子面的垂直距离为L，偏移距为x，道间距为Δx，直达波传播到隧道掌子面R处的时间记为T0；
当vp、vs确定，将对应tp、ts时刻振幅值进行叠加，完成所有道的计算后，即可得到速度谱估计速度。
5.根据权利要求1所述的一种隧道超前探测复杂地震波场的分离方法，其特征在于：所述规定的速度是采用叠加能量或相似系数准则之一，作为速度分析的准则制作速度谱，选择的速度使得校正后时距曲线刚好拉平，叠后能量最强；或者相似系数准则中的使得各道反射信息最为相似。
6.根据权利要求1所述的一种隧道超前探测复杂地震波场的分离方法，其特征在于：所述S4中，将傅里叶变换后的时间、频率、波数域数据固定为时间与频率，通过对空间各处波动情况进行傅里叶变换，得到时间频率波数与数据，同时观察时间、频率、波数特征；查看任意时刻的频率波数域数据，即使某时刻存在波形耦合，其频率、视速度方向也可能不同，选择能够克服F-K滤波器易造成假频缺点的滤波器，在每个时间点，进行二维视速度滤波，随着时间的增长，某时刻的频率波数域数据受直达波影响变小。
7.根据权利要求1所述的一种隧道超前探测复杂地震波场的分离方法，其特征在于：所述S5中，拾取反射纵波的主要步骤如下：
Q1、设计半时间窗口长度Δt、半频率窗口长度Δf，选择初始道号i及反射纵波大致出现的时间ti与频率fi；
Q2、在第i道时频数据中，于时间[ti-Δt,ti+Δt]，频率[fi-Δf,fi+Δf]区域内搜索极值点，将其对应的时间Ti作为反射纵波在第i道中对应的时间；
Q3、在第i±1道时频数据中，于时间[Ti-Δt,Ti+Δt]，频率[fi-Δf,fi+Δf]区域内搜索极值点，将其对应的时间Ti±1作为反射纵波在第i±1道中对应的时间；
Q4、重复以上步骤，直至搜索完所有道数据，得到反射纵波在地震记录中的时间节点
Q5、设计拾取半波形的半时间窗口长度ΔT，在滤波后的地震记录中，设计滤波器拾取时间窗口[Ti-ΔT,Ti+ΔT]内的信号，得到反射纵波；
在滤波前直接拾取直达波，绕射波在滤波后拾取，拾取方法均采用以上步骤。
8.根据权利要求1所述的一种隧道超前探测复杂地震波场的分离方法，其特征在于：所述S5中，所述半时间窗口长度Δt应在至少包含半波长范围内，选择最小的，避免在波形交叉时，拾取干扰波；所述波形的半时间窗口长度ΔT，应至少大于一个波长范围。
9.根据权利要求1所述的一种隧道超前探测复杂地震波场的分离方法，其特征在于：所述Q5中，反射纵波经过压制正视速度波后，根据方位角特征与视速度特征拾取反射纵波，针对拾取的反射纵波的波形，计算倾角谱，具体实现步骤如下：
Q5.1、拾取Y分量的反射纵波，得到含有单一波形的记录；
Q5.2、设置初始角度α0、扫描倾角间隔Δα和倾角最大值αm；
Q5.3、设置扫描时间窗口长度T；
Q5.4、当界面倾角为α时，反射波旅行时为t2，则其关系表示如下：
根据反射波旅行时t2与速度分析得到的纵波速度vp，计算每一个扫描角度αi对应的波的传播时间ti，将地震记录中该时刻对应的振幅值进行叠加，记录在Eα(ti,αi)中；
Q5.5、扫描过全部倾角后，Eα即为扫描的界面倾角谱，据此得到倾角最佳值，以及对应的自激自收时间t0。
10.根据权利要求1所述的一种隧道超前探测复杂地震波场的分离方法，其特征在于：根据倾角分析估算地层倾角后，计算同一反射界面对应激发的转换横波的时间点，若计算出的时间区域内存在同相轴，且该区域内的时频特征符合转换横波的特点，则验证反射纵波的拾取正确性，反之，如果该时间区域内未存在明显同相轴，应当重新根据波形特征进行拾取，重复之前的步骤进行验证。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911144423.XA CN110780346A (zh) | 2019-11-20 | 2019-11-20 | 一种隧道超前探测复杂地震波场的分离方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911144423.XA CN110780346A (zh) | 2019-11-20 | 2019-11-20 | 一种隧道超前探测复杂地震波场的分离方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN110780346A true CN110780346A (zh) | 2020-02-11 |
Family
ID=69392138
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911144423.XA Pending CN110780346A (zh) | 2019-11-20 | 2019-11-20 | 一种隧道超前探测复杂地震波场的分离方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110780346A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113640878A (zh) * | 2021-08-12 | 2021-11-12 | 西南石油大学 | 利用虚拟震源扫描构建方位角-视速度雷达图的方法 |
| CN115184992A (zh) * | 2022-08-18 | 2022-10-14 | 成都理工大学 | 基于自动拾取震相的远震震源快速定位方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6564150B2 (en) * | 1999-11-08 | 2003-05-13 | Board Of Regents For The University Of Texas | System and method for orienting seismic energy receivers to yield discriminated vertical shear waves |
| CN101246218A (zh) * | 2007-02-15 | 2008-08-20 | 中国石油化工股份有限公司 | 一种三分量vsp波场分离方法 |
| CN102103215A (zh) * | 2009-12-16 | 2011-06-22 | 中国石油天然气集团公司 | 三维高密度地震勘探记录叠前面波压制方法 |
| CN103645497A (zh) * | 2013-11-08 | 2014-03-19 | 中国石油大学(北京) | 一种基于出射角的多分量波场分离方法 |
| US20160320509A1 (en) * | 2015-04-30 | 2016-11-03 | Saudi Arabian Oil Company | Suppressing near-surface scattered surface waves |
| CN106249288A (zh) * | 2016-08-02 | 2016-12-21 | 中国地质大学（北京） | 基于Shearlet域的极化滤波面波压制方法 |
| CN107144880A (zh) * | 2017-05-12 | 2017-09-08 | 招商局重庆交通科研设计院有限公司 | 一种地震波波场分离方法 |
-
2019
- 2019-11-20 CN CN201911144423.XA patent/CN110780346A/zh active Pending
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6564150B2 (en) * | 1999-11-08 | 2003-05-13 | Board Of Regents For The University Of Texas | System and method for orienting seismic energy receivers to yield discriminated vertical shear waves |
| CN101246218A (zh) * | 2007-02-15 | 2008-08-20 | 中国石油化工股份有限公司 | 一种三分量vsp波场分离方法 |
| CN102103215A (zh) * | 2009-12-16 | 2011-06-22 | 中国石油天然气集团公司 | 三维高密度地震勘探记录叠前面波压制方法 |
| CN103645497A (zh) * | 2013-11-08 | 2014-03-19 | 中国石油大学(北京) | 一种基于出射角的多分量波场分离方法 |
| US20160320509A1 (en) * | 2015-04-30 | 2016-11-03 | Saudi Arabian Oil Company | Suppressing near-surface scattered surface waves |
| CN106249288A (zh) * | 2016-08-02 | 2016-12-21 | 中国地质大学（北京） | 基于Shearlet域的极化滤波面波压制方法 |
| CN107144880A (zh) * | 2017-05-12 | 2017-09-08 | 招商局重庆交通科研设计院有限公司 | 一种地震波波场分离方法 |
Non-Patent Citations (1)
| Title |
|---|
| 曾胜强 等: "隧道反射地震超前预报波场分离方法的数值模拟研究", 《地球物理学进展》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113640878A (zh) * | 2021-08-12 | 2021-11-12 | 西南石油大学 | 利用虚拟震源扫描构建方位角-视速度雷达图的方法 |
| CN113640878B (zh) * | 2021-08-12 | 2024-03-29 | 西南石油大学 | 利用虚拟震源扫描构建方位角-视速度雷达图的方法 |
| CN115184992A (zh) * | 2022-08-18 | 2022-10-14 | 成都理工大学 | 基于自动拾取震相的远震震源快速定位方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN102879817B (zh) | 基于地面地震数据获取地下裂缝信息的控制方法 | |
| CN103675910B (zh) | 一种水陆检波器地震数据标定因子反演方法 | |
| CA2783055A1 (en) | Method and device for wave fields separation in seismic data | |
| CN111025386B (zh) | 一种无分离假象的纵横波分离方法 | |
| CN116520419B (zh) | 一种热流体裂缝通道识别方法 | |
| CN112327362B (zh) | 速度域的海底多次波预测与追踪衰减方法 | |
| CN112305591B (zh) | 隧道超前地质预报方法、计算机可读存储介质 | |
| CN105510975B (zh) | 提高地震数据信噪比的方法及装置 | |
| CN110389377B (zh) | 基于波形互相关系数相乘的微震偏移成像定位方法 | |
| CN102073064A (zh) | 一种利用相位信息提高速度谱分辨率的方法 | |
| CN104570116A (zh) | 基于地质标志层的时差分析校正方法 | |
| CN110780346A (zh) | 一种隧道超前探测复杂地震波场的分离方法 | |
| CN114415234B (zh) | 基于主动源面波频散和h/v确定浅地表横波速度的方法 | |
| CN110780341B (zh) | 一种各向异性地震成像方法 | |
| CN114265118A (zh) | 一种用于海上随钻声波测井时差提取方法、装置及系统 | |
| CN107479091B (zh) | 一种提取逆时偏移角道集的方法 | |
| CN108693560B (zh) | 一种基于互相关道的散射波成像方法及系统 | |
| CN110579799B (zh) | 一种等旅行时间间隔的地震采集观测方法及系统 | |
| CN108919345B (zh) | 一种海底电缆陆检噪声的衰减方法 | |
| CN111352153A (zh) | 一种基于瞬时相位互相关加权的微地震干涉定位方法 | |
| CN108802822B (zh) | 方位各向异性介质中的保幅直接叠前时间偏移方法及装置 | |
| CN112526611A (zh) | 表层地震波品质因子的提取方法及装置 | |
| CN114114407A (zh) | 一种地震波探测用表面波和直达横波压制处理方法 | |
| Zhou et al. | Surface diffraction noise attenuation for marine seismic data processing with mathematical morphological filtering | |
| CN112327361B (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 | ||
| WD01 | Invention patent application deemed withdrawn after publication | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20200211 |