CN111103620A - 一种岩巷超前探测三维偏移成像方法 - Google Patents

一种岩巷超前探测三维偏移成像方法 Download PDF

Info

Publication number
CN111103620A
CN111103620A CN201911143406.4A CN201911143406A CN111103620A CN 111103620 A CN111103620 A CN 111103620A CN 201911143406 A CN201911143406 A CN 201911143406A CN 111103620 A CN111103620 A CN 111103620A
Authority
CN
China
Prior art keywords
polarization
point
velocity
interface
wave
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
Application number
CN201911143406.4A
Other languages
English (en)
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.)
Individual
Original Assignee
Individual
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 Individual filed Critical Individual
Priority to CN201911143406.4A priority Critical patent/CN111103620A/zh
Publication of CN111103620A publication Critical patent/CN111103620A/zh
Pending legal-status Critical Current

Links

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/282—Application of seismic models, synthetic seismograms
    • 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/303—Analysis for determining velocity profiles or travel times
    • 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/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time

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

本发明提供了一种岩巷超前探测三维偏移成像方法，通过获取三分量地震观测数据，利用三维速度分析公式，得到三维速度谱，获得从震源发出的波经过每个绕射点到各检波器的波场绕射时间和绕射射线出射方向，利用绕射射线出射方向与主偏振方向之间的差异构造方向权函数；在均匀各向同性介质中的克希霍夫积分法偏移的过程中，加入方向权函数项，得到极化深度偏移成像结果。本发明采用三维速度扫描分析来估算速度，能够更加准确的估计隧道介质速度，并利用初至纵波的极化信息来确定反射界面的几何参数，两者结合能够准确的确定初始速度模型。在克希霍夫积分偏移中引入了方向权函数，能够有效的压制隧道中克希霍夫积分成像方法中产生的偏移假象。

Description

一种岩巷超前探测三维偏移成像方法
技术领域
本发明属于隧道地震地质预报技术领域，具体涉及利用三分量地震数据的 矢量特性以及不同类型波的极化特征，来推测反射界面的空间方位的一种岩巷 超前探测三维偏移成像方法。
背景技术
随着我国交通、水利水电事业的快速发展，大规模的铁路、公路和水电工 程建设数量不断增加，需要修建大量各种类型的隧道，在隧道施工建设中，掌 子面前方常隐含有不利地质构造对象，在开挖揭露过程中容易诱发地质灾害导 致安全事故，给工程建设造成重大损失，为了准确确定隧道掘进前方的异常地 质构造，就要在隧道的超前地震预报数据处理过程中，得到精准的偏移成像结 果。
目前，岩巷超前地震探测中较常用的偏移方法是克希霍夫积分偏移，将射 线追踪方法与Kirchhoff积分偏移结合起来，可以实现叠前Kirchhoff深度偏 移来进行隧道地质预报；通过分析隧道环境下的地震波传播的时距关系，可以 通过插值的方法使用不同精度的网格剖分，以提高运算速度。由于隧道施工过 程中地下空间地质构造复杂且容易出现不连续的情况，从物理地震学的角度来 看，观测到的波场信息可以认为是地下所有点产生的绕射在观测点处叠加的结 果，因此，现有技术也有将散射成像技术应用于隧道地震超前地质预报中，当 观测系统的横向展布不够时，按照散射原理偏移，偏移结果中会呈现往外扩展 的圆弧，可以通过方向扫描的方法来改进这一缺陷。
现有的各种隧道地震超前预报技术，仍存在一些问题：首先，在隧道地震 超前探测中，采用直达波估计得到的速度只代表隧道侧壁围岩表层介质的速 度，并不能反映波场在反射界面到测线之间的介质中的地震波传播速度。传统 速度分析方法以水平层状介质模型假设为前提，不适用于隧道环境下的速度估 计；其次，由于隧道施工环境限制，地震观测测线多沿隧道轴线布置于地下空 间，当沿用地面地震勘探中的偏移方法对数据进行偏移处理时，会产生关于测 线对称的偏移假象。
发明内容
根据上述阐述，本发明的目的在于提供一种岩巷超前探测三维偏移成像方 法，解决偏移成像结果中存在假象的问题。
本发明主要原理是利用波场的极化信息，引入了方向权函数来压制克希霍 夫积分偏移的假象，其具体提供的技术方案为：
一种岩巷超前探测三维偏移成像方法，包含如下步骤：
S1、获取三分量地震观测数据，将地震观测数据中的直达波剔除，并分离 纵波；
S2、利用三维速度分析公式，对预处理后的地震观测数据进行速度扫描分 析，得到三维速度谱，以估计地下介质地震波速度，再利用初至纵波的极化信 息确定反射界面的几何参数，利用三维速度谱和反射界面的几何参数确定速度 模型；
S3、确定时窗计算长度W，利用协方差矩阵特征分析法，求取三分量地震 记录的波场极化信息，得到对应绕射时刻的主偏振方向；
S4、根据步骤S2得到的速度模型，利用射线追踪的方法，计算获得：从 震源发出的波经过每个绕射点到各检波器的波场绕射时间和绕射射线出射方 向；
S5、利用绕射射线出射方向与主偏振方向之间的差异构造方向权函数；
S6、在均匀各向同性介质中的克希霍夫积分法偏移的过程中，加入方向权 函数项，得到极化深度偏移成像结果。
上述技术方案中，所述步骤S2中：三维速度分析公式ESpec，具体形式如 下：
Figure BDA0002281538920000031
其中，i表示道序号，j表示采样点在窗口内的自然排列顺序序号，
Figure BDA0002281538920000032
表示各道的反射时间，M表示速度扫描计算时窗 长度，t0表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向 单位向量的Y分量；
首先计算各道的反射时间ti，再依据公式计算与t0、v、m对应的速度谱值， 对上述三个变量选定不同值，重复以上步骤即可以计算地震剖面的三维速度 谱，通过三维速度谱，估计反射界面到炮点之间的介质的速度大小。
上述技术方案中，所述步骤S2中：利用初至纵波的极化信息确定反射界 面几何参数的公式，具体形式如下：
Figure BDA0002281538920000041
其中，N表示假设已知的直线束Li的数量，(i＝1，2，3，…，N)，Pi(xi,yi,zi) 表示每一直线上的某一已知坐标的点，Q(x,y,z)表示空间内的任意点， vi＝(li,mi,ni)表示各直线的单位方向向量，I表示N阶单位矩阵；
具体为利用极化特征，确定地下界面的法线方向及其位置，先对波的类型 做出判别，再对纵波进行极化分析，以获取波场在不同位置的传播方向，然后 通过波场在不同位置的传播方向进行反向延长，求取交点，通过解如上线性方 程，求得反射射线的反向交点Q的坐标，交点与炮点关于界面镜像对称，由 交点与炮点确定的直线的单位方向的向量v为界面法线方向，所述界面法线方 向及交点与炮点构成线段的中点M(xM,yM,zM)，确定反射界面的解析表达式； 其中，点(xp,yp,zp)为空间平面上不同于点M的任一点。
Figure BDA0002281538920000042
在计算过程中给出三维速度模型空间的大小限制，利用空间平面与速度模 型空间的各棱边直线的交点来估计反射界面的具体位置，从而构建初始速度模 型。
上述技术方案中，所述三维速度模型空间包含反射界面内所有的目标体。
上述技术方案中，所述步骤S3中：
所述协方差矩阵，具体形式如下：
Figure BDA0002281538920000051
其中，W表示时窗计算长度，u(x,y,z)表示三分量地震矢量波场，
Figure BDA0002281538920000052
表示方差，
Figure BDA0002281538920000053
表示协方差，
Figure BDA0002281538920000054
表示数学期望，
Figure BDA0002281538920000055
表示采 样点数，N表示每道的总采样点数，ΔT表示时间采样间隔；
协方差矩阵为一对称矩阵，求取协方差矩阵的特征值λ1、λ2、λ3，且限定 λ1≥λ2≥λ3，及其对应的单位特征向量为p1、p2、p3，协方差矩阵的三个特征值 分别对应极化椭圆的三个轴的长度，最大特征值λ1对应极化椭圆的主轴，相应 的主特征向量p1表示质点椭圆偏振的主方向，其余同理；
该时刻波的偏振类型用偏振系数γ(t)来衡量：
Figure BDA0002281538920000056
若γ值趋于1，主特征值λ1会远大于其它次特征值，并表现为极化椭球 趋于细长，表示质点的空间振动过程趋于线性极化，此时主特征向量所对应的 空间方向即为质点振动所在的方向；若γ值趋于0，主特征值会和次特征值相 当，极化椭球表现为趋于球形，表示质点在空间振动时的非线性极化越强，此 时质点振动方向没有一个优势方向，通过窗口在地震记录上的滑动，求取整个 剖面上不同采样点处，不同时刻的波的极化特征信息。
上述技术方案中，所述步骤S5中，所述方向权函数的具体形式如下：
Figure RE-GDA0002371817980000061
其中，Θ的具体形式如下：
Figure BDA0002281538920000062
θ表示射线追踪计算的绕射射线在检波器处的出射方向和根据绕射时间 对地震数据进行极化分析所得的纵波极化方向两方向之间的方向夹角。
上述技术方案中，所述步骤S6中：
均匀各向同性介质中的克希霍夫积分解的具体形式为：
Figure BDA0002281538920000063
其中，第一项由波场的垂直梯度决定，第二项与波场函数值的延迟位有关，叫 做近场源项，它以
Figure BDA0002281538920000064
衰减，这两项通常在克希霍夫积分法偏移中被忽略，而 只使用第三项远场源项，第三项可以用波场观测数据项和权因子项W来表征：
Figure BDA0002281538920000065
其中，
Figure BDA0002281538920000066
表示权因子，即在积分过程中，各道观测值对最终各成像点幅 值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S表示 积分曲面；
加入了方向权函数的克希霍夫积分偏移公式，具体形式如下：
Figure BDA0002281538920000067
其中，κ表示调节因子，用于调节权因子在成像过程中的压制效果，K(p,d)表 示方向权函数，
Figure BDA0002281538920000071
表示权因子，为积分过程中各道观测值对最终各成像 点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S 表示积分曲面，由地面和一个向地下空间无限延伸的部分球面组成，同时，在 无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。
本发明与现有的岩巷克希霍夫积分偏移成像技术相比，优点为：采用 三维速度扫描分析来估算速度，能够更加准确的估计隧道介质速度，并利 用初至纵波的极化信息来确定反射界面的几何参数，两者结合能够准确的 确定初始速度模型。在克希霍夫积分偏移中引入了方向权函数，能够有效 的压制隧道中克希霍夫积分成像方法中产生的偏移假象。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施 例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述 中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付 出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例的流程图；
图2为虚实反射界面的反射纵波传播方向与反射纵波极化方向的差异示 意图；
图3为与测线方向成小角度斜交的竖直断层模型A的示意图；
图4为模型A经波场分离后的三分量纵波地震记录示意图；
图5为模型A纯纵波地震记录经速度扫描所得的三维速度谱示意图；
图6为模型A三维速度扫描结果ZOX平面切片图，其中(a)～(f)依次 对应于界面法线方向与Y轴夹角范围45°～70°，角度间隔5°；
图7为模型A极化分析结果图，其中(a)表示不同时刻的极化程度参数 值，(b)～(d)表示主特征向量沿X、Y、Z三轴的方向余弦值；
图8为三维单倾斜界面地震观测系统
图9为利用真速度模型A进行极化偏移在深度z＝25m处XOY截面的成像结 果；
图10为利用速度分析结果构建的速度模型A，得到的X分量地震记录采 用不同方向权函数时的偏移成像结果；
图11为利用速度分析结果构建的速度模型A，得到的Y分量地震记录采 用不同方向权函数时的偏移成像结果；
图12为利用速度分析结果构建的速度模型A，得到的Z分量地震记录采 用不同方向权函数时的偏移成像结果；
图13为图9、图10、图11中，切片深度z＝25m的偏移成像结果切片图；
图14为与测线方向成大角度斜交的竖直断层模型B的示意图；
图15为模型B经波场分离后的三分量纵波地震记录示意图；
图16为模型B纯纵波地震记录经速度扫描所得的三维速度谱示意图；
图17为模型B三维速度扫描结果ZOX平面切片图，其中(a)～(f)依 次对应于界面法线方向与Y轴夹角范围45°～70°，角度间隔5°；
图18为模型B极化分析结果图，其中(a)表示不同时刻的极化程度参数 值，(b)～(d)表示主特征向量沿X、Y、Z三轴的方向余弦值；
图19为利用真速度模型B进行非极化偏移的结果，方向权因子K＝1，图 (a)～(c)依次对应反射纵波X、Y、Z分量地震记录的成像结果；
图20为利用真速度模型B进行极化偏移的结果，方向权因子K＝f5，指数 κ＝1，(a)～(c)依次对应反射纵波X、Y、Z分量地震记录的成像结果；
图21为利用速度分析结果构建的速度模型B进行极化偏移的结果，方向 权因子K＝f5，指数κ＝1，(a)～(c)依次对应反射纵波X、Y、Z分量地震记 录的成像结果；
图22为图19、图20、图21的偏移成像结果切片图，切片深度z＝25m。
具体实施方式
下面将结合本发明的附图，对本发明的技术方案进行清楚、完整地描述， 显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基 于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获 得的所有其他实施例，都属于本发明保护的范围。
根据图1所示，本发明提供的一种岩巷超前探测三维偏移成像方法，包含 如下步骤：
S1、获取三分量地震观测数据，将地震观测数据中的直达波剔除，并分离 纵波；
S2、三维速度扫描分析和反射界面几何参数估计，利用三维速度分析公式， 对预处理后的地震观测数据进行速度扫描分析，得到三维速度谱，以估计地下 介质地震波速度，再利用初至纵波的极化信息确定反射界面的几何参数，利用 三维速度谱和反射界面的几何参数确定速度模型；
三维速度分析公式ESpec，具体形式如下：
Figure BDA0002281538920000091
其中，i表示道序号，j表示采样点在窗口内的自然排列顺序序号，
Figure BDA0002281538920000101
表示各道的反射时间，M表示速度扫描计算时窗 长度，t0表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向 单位向量的Y分量；
首先计算各道的反射时间ti，再依据公式计算与t0、v、m对应的速度谱值， 对上述三个变量选定不同值，重复以上步骤即可以计算地震剖面的三维速度 谱，通过三维速度谱，估计反射界面到炮点之间的介质的速度大小。
利用初至纵波的极化信息确定反射界面几何参数的公式，具体形式如下：
Figure BDA0002281538920000102
其中，N表示假设已知的直线束Li的数量，(i＝1，2，3，…，N)，Pi(xi,yi,zi)表 示每一直线上的某一已知坐标的点，Q(x,y,z)表示空间内的任意点， vi＝(li,mi,ni)表示各直线的单位方向向量，I表示N阶单位矩阵；
具体为利用极化特征，确定地下界面的法线方向及其位置，先对波的类型 做出判别，再对纵波进行极化分析，以获取波场在不同位置的传播方向，然后 通过波场在不同位置的传播方向进行反向延长，求取交点，通过解如上线性方 程，求得反射射线的反向交点Q的坐标，交点与炮点关于界面镜像对称，由 交点与炮点确定的直线的单位方向的向量v为界面法线方向，所述界面法线方 向及交点与炮点构成线段的中点M(xM,yM,zM)，确定反射界面的解析表达式； 其中，点(xp,yp,zp)为空间平面上不同于点M的任一点。
Figure BDA0002281538920000111
若能在计算过程中给出三维速度模型空间(或偏移目标区域)的大小限制， 利用空间平面与速度模型空间的各棱边直线的交点来估计反射界面的具体位 置，从而构建初始速度模型，三维速度模型空间应当包含反射界面内所有的目 标体，例如野外工作时，假设地下目标体形状为一个正方体，大小为 100m*100m*100m，此时速度模型空间只要大于该目标体的大小即可，例如可 以设置为120m*120m*120m，总而言之，只要模型空间大小包括所有需要被 测量的目标体即可。
S3、求取地震记录的波场极化信息，得到对应绕射时刻的主偏振方向，即 确定时窗计算长度W，利用协方差矩阵特征分析法，求取三分量地震记录的波 场极化信息，得到对应绕射时刻的主偏振方向；
所述协方差矩阵，具体形式如下：
Figure BDA0002281538920000112
其中，W表示时窗计算长度，u(x,y,z)表示三分量地震矢量波场，
Figure BDA0002281538920000113
表示方差，
Figure BDA0002281538920000114
表示协方差，
Figure BDA0002281538920000115
表示数学期望，
Figure BDA0002281538920000116
表示采 样点数，N表示每道的总采样点数，ΔT表示时间采样间隔；
协方差矩阵为一对称矩阵，求取协方差矩阵的特征值λ1、λ2、λ3，且限定 λ1≥λ2≥λ3，及其对应的单位特征向量为p1、p2、p3，协方差矩阵的三个特征值 分别对应极化椭圆的三个轴的长度，最大特征值λ1对应极化椭圆的主轴，相应 的主特征向量p1表示质点椭圆偏振的主方向，其余同理；
该时刻波的偏振类型用偏振系数γ(t)来衡量：
Figure BDA0002281538920000121
若γ值趋于1，主特征值λ1会远大于其它次特征值，并表现为极化椭球 趋于细长，表示质点的空间振动过程趋于线性极化，此时主特征向量所对应的 空间方向即为质点振动所在的方向；若γ值趋于0，主特征值会和次特征值相 当，极化椭球表现为趋于球形，表示质点在空间振动时的非线性极化越强，此 时质点振动方向没有一个优势方向，通过窗口在地震记录上的滑动，求取整个 剖面上不同采样点处，不同时刻的波的极化特征信息。
S4、根据步骤S2得到的速度模型，利用射线追踪的方法，计算获得：从 震源发出的波经过每个绕射点到各检波器的波场绕射时间和绕射射线出射方 向；
S5、利用绕射射线出射方向与主偏振方向之间的差异构造方向权函数；
S6、在均匀各向同性介质中的克希霍夫积分法偏移的过程中，加入方向权 函数项，得到极化深度偏移成像结果。
岩巷超前探测三维极化偏移成像方法的基本原理为：已知均匀各向同性介 质中的标量波动方程如下：
Figure BDA0002281538920000131
其中，U(x,y,z,t)描述空间中P(xP,yP,zP)点处的波场函数，v(x,y,z)描述了波传播空间中的速度分布情况，F＝F(x,y,z)表示求解空间中的已知震源函数。
根据惠更斯-菲涅尔原理：某时刻波前面上的所有质点可以看作广义绕射 震源，从这些震源发出的二次子波在下一时刻的一系列位置产生干涉叠加，这 些二次子波干涉叠加的包络即代表了该时刻波前面的新位置。该原理描述了这 样一个事实：当不知道震源确切情况时，若已知某空间的某些点(已知波前面) 的振动情况，可求该空间内待求目标点的振动情况。
当在地面进行地震勘探时，对于三维标量波动方程的克希霍夫积分解中的 积分曲面，可以认为其由地面和一个向地下空间无限延伸的部分球面组成，在 无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。定义坐标系为 XOY面与地面重合，Z轴垂直地面指向地下空间为其正方向的情况下，对地面 的外法线方向的导数与对z坐标的导数只相差一个符号，即
Figure BDA0002281538920000132
此 处θ表示r与地面外法线方向n之间的夹角。当积分曲面内不包含有震源时， 即取F＝0，因此均匀各向同性介质中的克希霍夫积分解的具体形式为：
Figure BDA0002281538920000133
其中，第一项由波场的垂直梯度决定，第二项与波场函数值的延迟位有关，叫 做近场源项，它以
Figure BDA0002281538920000134
衰减，这两项通常在克希霍夫积分法偏移中被忽略，而 只使用第三项远场源项，第三项可以用波场观测数据项和权因子项W来表征：
Figure BDA0002281538920000135
其中，
Figure BDA0002281538920000136
表示权因子，即在积分过程中，各道观测值对最终各成像点幅 值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S表示 积分曲面。
根据图2所示，XOY面表示观测测线所在的水平面，并与地面平行，OR 表示垂直于地面的地质分界面，直线a表示是地震观测测线所在的直线，该直 线平行于隧道轴线，且与地面平行，虚线b所示为与地质界面关于测线所在竖 直面对称的假想界面，S表示炮点，G表示测线上任一检波点，检波点数量为 N。S*和S′分别表示炮点S关于真实分界面和假想分界面的对称点；带箭头 的实线c是根据射线定律所作的反射射线路径，同时也表示了波场的传播方 向；带箭头的实线d则表示与之对应的假想射线路径。双向箭头e表示反射纵 波在检波点处的极化方向。在炮点S(xS,yS,zS)激发的地震波，在真实地下界面 OR上发生反射，反射地震波场向着测线所在方向传播，当反射波场到达测线 上的检波点时，介质的震动状态被布置在该处的检波器记录下来。设点 Rijk(xR,yR,zR)为地下空间中的任意网格点(i、j、k为对地下空间进行网格化后 沿X、Y、Z三轴方向上的网格序号)，这些网格点均可作为绕射点存在，当其 处于跟炮检点位置有关的真实反射界面上时，表示其成为一反射点。在均匀介 质中，地震波沿直线方向传播，地震波的传播路径可以用射线路径S→R→G 来描述。第n个地震道Gn(xg,yg,zg)(检波点序号n＝1，2，…，N)点处检波器 接收到的震动信息相当于产生于炮点的振动信号滞后一个时间延迟后在G点 出现，延迟时间t可以通过地震记录上同相轴的位置来读取。显然，绕射射线 R→G的方向代表了R点的反射(或绕射)波动传播到检波器G的方向。
对于给定的已知速度模型，将模型中的每一个网格点看作时绕射点，而对 于有一定延续长度(面积)的界面，可以将其看作是由许多绕射点组成的一个 整体，从而可以把研究对象统一起来，使用各点的绕射叠加结果来表示反射信 息。对每个网格点进行射线追踪可以找到从每一个网格点(绕射点)到检波点 的绕射射线单位方向向量d(Rijk,Gn)。显然，对于某一个网格点来说，网格点到 各个检波点的方向不尽相同。
对于观测系统(S、G位置关系已知)及速度模型已知的情形，设由模型 中任一网格点Rijk到第n个地震道的射线追踪时间为Tn(S,Gn,Rijk)，从而，根据 该时间可以从地震记录上获取射线到各道的幅值信息。若某个网格点对应某一 真实反射点，则由射线追踪计算的对应于各个检波点的旅行时间Tn，可以从地 震记录上获取到一系列的振幅值un(S,Gn,Tn)。但是，在反射关系下，这N条绕 射射线中仅仅有唯一的真实反射射线R→Gk，该射线对应的振幅值uk(S,Gk,Tk) 为一峰值，其绕射射线的出射方向单位向量dk表示真实的纵波波场传播方向， 而其它射线对应的振幅值un(S,Gn,Tn)(其中n≠k)应该趋近于零值，它们相应 的绕射射线出射方向单位向量dn(其中n≠k)也与真实的波场传播方向存在 差别。
根据上述分析，可以在克希霍夫积分法偏移的过程中加入上述两个方向的 差别有关的权函数K(p,d)，加入了方向权函数的克希霍夫积分偏移公式，具体 形式如下：
Figure BDA0002281538920000151
其中，κ表示调节因子，用于调节权因子在成像过程中的压制效果，K(p,d)表 示方向权函数，
Figure BDA0002281538920000152
表示权因子，为积分过程中各道观测值对最终各成像 点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S 表示积分曲面，由地面和一个向地下空间无限延伸的部分球面组成，同时，在 无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。
方向权函数的具体形式如下：
Figure RE-GDA0002371817980000161
其中，Θ的具体形式如下：
Figure BDA0002281538920000162
θ表示射线追踪计算的绕射射线在检波器处的出射方向和根据绕射时间 对地震数据进行极化分析所得的纵波极化方向两方向之间的方向夹角。
实施例采用模型A和模型B，模型A为与测线方向成小角度斜交的竖直断 层模型；模型B为与测线方向成大角度斜交的竖直断层模型，两个速度模型的 大小均为x×y×z＝100m×100m×50m，模型在三个轴方向上的离散网格间距均为 0.5m。
模型A与模型B参数设置如表1所示。
表1、模型A和模型B参数设置
Figure BDA0002281538920000163
在图3中为模型A，即与测线方向成小角度斜交的竖直断层模型，断层两 侧为岩性不同的两种介质，断层界面与YOZ平面所成角度角约为26.5°，分界 面法向单位向量为
Figure BDA0002281538920000164
在图3中A部分表示低速弹性介质1；B部分表示高速弹性介质2，平行 于Y轴方向的直线D1表示地震探测测线方向。黑色小球体表示炮点S位置， 其坐标为(50m，5m，25m)，直线D1上的白色小球体表示部分检波器分布的 示意图，偏移距0.5m，道间距0.5m，共101道接收，检波器分布范围5.5～55.5m。 时间采样间隔Δt＝0.05ms，采样总时长为400ms，正演数值模拟采用有限差分 法进行，震源类型为爆炸震源，子波类型为主频300Hz的Ricker子波。
具体实施方式为，首先将三分量纵波地震记录进行波场分离，结果如图4 所示，图4中的同相轴代表被分离出来的纵波，图4(a)、图4(b)、图4(c) 依次对应X、Y、Z三分量地震记录。然后将图4所示的分离纵波后的地震记录 进行速度扫描得到三维速度谱，结果如图5所示，图5表示了介质的速度分布 情况。
图6为模型A三维速度扫描结果ZOX平面切片图，图6(a)～图6(f) 依次对应于界面法线方向与Y轴夹角范围45°～70°，角度间隔5°；
三维速度扫描的具体做法如下：
参考图8，其中X轴指向正北，Y轴指向正东方向，XOY面平为地面，Z 轴垂直地面向下为正。平面ABC表示地下某一倾斜界面，界面法线方向单位 向量为n＝(l,m,n)，其中l＝cos(n,X)，m＝cos(n,Y)＝cosθY，n＝cos(n,Z)。界面上 方介质速度为v。在隧道工程建设中时常遇到类似上图所示的地下地质构造面。 图中紫色直线表示布置在地下空间的沿Y轴方向延伸的直测线，测线穿过界面 ABC。炮点表示为S(xS,yS,zS)，测线上检波点表示为G(xG,yG,zG)。SS*垂直面 ABC于M点，且S*为炮点S关于空间界面ABC的对称点，坐标为
Figure BDA0002281538920000171
SR、RG表示从震源S发出的某一射线在R点反射后到达检波点G点，由反射 定律可知，连接S*、G两点的连线与界面ABC的交点即为反射点R(xR,yR,zR)。
为了对图8中所示几何关系进行定量表示，首先，假设观测测线SG与XOZ 平面的交点坐标为(x0,0,z0)，根据上述几何关系，有：
Figure BDA0002281538920000181
其次，假设炮点S的零偏移距自激自收旅行时间为t0，且炮点到界面的垂 直距离为H，则有下述关系：
2H＝vt0
最后，假设炮点关于界面的镜像对称点S*坐标为(x′,y′,z′)。则有向量
Figure BDA0002281538920000182
显然，向量
Figure BDA0002281538920000183
与界面的法线方向n平行，因此有如 下关系式：
Figure BDA0002281538920000184
根据上述关系，可以用法线方向余弦表示各检波点处的反射旅行时间：
Figure BDA0002281538920000185
即：
Figure BDA0002281538920000186
从上式可以看出，当测线沿Y轴方向时，测线上任意点的反射旅行时间可 以由炮点自激自收时间t0、炮检点距离(yS-yG)以及界面单位法线方向在Y轴 方向上的投影值m三个参数完全表示，从而能在速度扫描过程中，避免多参 数引起的不便，也便于成图显示和分析。对于时间采样间隔为Δt共N道的共 炮点道集地震记录U(x,t)，设速度扫描计算时窗长度为M个采样点，依据平均 振幅能量准则，在经过速度扫描过程后，可以得到三维速度谱ESpec，具体形式 如下：
Figure BDA0002281538920000191
其中，i表示道序号，j表示采样点在窗口内的自然排列顺序序号，
Figure BDA0002281538920000192
表示各道的反射时间，M表示速度扫描计算时窗 长度，t0表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向 单位向量的Y分量；
首先计算各道的反射时间ti，再依据公式计算与t0、v、m对应的速度谱值， 对上述三个变量选定不同值，重复以上步骤即可以计算地震剖面的三维速度 谱。
利用协方差矩阵法进行极化分析，如图7所示，求得极化程度参数值和主 特征向量的方向余弦值，图7(a)为计算得到的极化程度，其显示了在同相 轴处极化程度参数γ值接近1，表示在该同相轴所在的时间宽度范围内波场高 度极化，而在只有噪声存在的地方，参数γ值接近0，波场极化程度低。根据 极化方向的方向余弦图7(b)、(c)、图7(d)可知，沿空间X轴方向的方向 余弦整体上大于沿Y方向的方向余弦，而沿Z轴方向的方向余弦接近于零值， 即通过模型1模拟得到的波场中的纵波成分主要集中在X分量上，而Z分量 上的振幅值几乎为零，通过图7(c)可知纵波波场的Y分量存在一渐变过程： 当检波器越远离炮点位置，纵波波场的Y分量由大变小再增加，且对应的方 向余弦在70道附近发生符号的改变。
具体做法如下：
Figure BDA0002281538920000201
其中，W表示时窗计算长度，u(x,y,z)表示三分量地震矢量波场，
Figure BDA0002281538920000202
表示方差，
Figure BDA0002281538920000203
表示协方差，
Figure BDA0002281538920000204
表示数学期望，
Figure BDA0002281538920000205
表示采 样点数，N表示每道的总采样点数，ΔT表示时间采样间隔；
协方差矩阵为一对称矩阵，求取协方差矩阵的特征值λ1、λ2、λ3，且限定 λ1≥λ2≥λ3，及其对应的单位特征向量为p1、p2、p3，协方差矩阵的三个特征值 分别对应极化椭圆的三个轴的长度，最大特征值λ1对应极化椭圆的主轴，相应 的主特征向量p1表示质点椭圆偏振的主方向，其余同理；
该时刻波的偏振类型用偏振系数γ(t)来衡量：
Figure BDA0002281538920000206
若γ值趋于1，主特征值λ1会远大于其它次特征值，并表现为极化椭球趋 于细长，表示质点的空间振动过程趋于线性极化，此时主特征向量所对应的空 间方向即为质点振动所在的方向；若γ值趋于0，主特征值会和次特征值相当， 极化椭球表现为趋于球形，表示质点在空间振动时的非线性极化越强，此时质 点振动方向没有一个优势方向，通过窗口在地震记录上的滑动，求取整个剖面 上不同采样点处，不同时刻的波的极化特征信息。
然后通过地震极化信息求模型A的反射界面参数，如表2所示。
表2模型A极化分析计算得到的反射界面参数
Figure BDA0002281538920000211
通过地震极化信息求取反射界面几何参数的具体做法如下：
设已知直线束Li，(i＝1，2，3，…，N)，每一直线上都有一已知坐标的点 Pi，其坐标为(xi,yi,zi)，且各直线的单位方向向量为vi＝(li,mi,ni)也已知。设 Q(x,y,z)为空间内的任意点，其到某一直线的距离hi可以表示为：
Figure BDA0002281538920000212
其中i、j、k为沿X、Y、Z三个坐标轴正方向的单位向量，且有：
α＝(y-yi)ni-(z-zi)mi
β＝(z-zi)li-(x-xi)ni
γ＝(x-xi)mi-(y-yi)li
为了求得直线束Li的公共交点，取点Q到所有直线的距离之和
Figure BDA0002281538920000213
为最小，即令：
Figure BDA0002281538920000214
且
Figure BDA0002281538920000215
且
Figure BDA0002281538920000216
Figure BDA0002281538920000221
若令l、m、n分别表示如下列向量：
Figure BDA0002281538920000222
又因为vi为单位向量，有上式左侧系数矩阵可进一步 化为(其中I为N阶单位矩阵)：
Figure BDA0002281538920000224
因此，方程可化为：
Figure BDA0002281538920000225
其中，N表示假设已知的直线束Li的数量，(i＝1，2，3，…，N)，Pi(xi,yi,zi)表 示每一直线上的某一已知坐标的点，Q(x,y,z)表示空间内的任意点， vi＝(li,mi,ni)表示各直线的单位方向向量，I表示N阶单位矩阵；
对纵波进行极化分析，获取波场在不同位置的传播方向后，通过波场在不 同位置的传播方向进行反向延长，求取交点，通过解如上线性方程，求得反射 射线的反向交点Q的坐标，交点与炮点关于界面镜像对称，由交点与炮点确 定的直线的单位方向的向量v为界面法线方向，所述界面法线方向及交点与炮 点构成线段的中点M(xM,yM,zM)，确定反射界面的解析表达式；
通过野外采集的三分量地震数据空间坐标范围大小可以给出三维速度模 型空间(或偏移目标区域)的大小限制，利用空间平面与速度模型空间的各棱 边直线的交点来估计反射界面的具体位置，通过求取空间反射界面与给定模型 A空间棱边的交点，得到表3所示的四个交点坐标。
表3模型A空间棱边的交点
Figure BDA0002281538920000231
通过表2与表3上述求得的反射界面几何参数与图5中三维速度谱结合， 即可得到经速度分析后的偏移速度模型。
已知速度模型后，利用射线追踪方法可以得到旅行时与绕射波出射方向， 利用绕射波出射方向和主偏振方向之间的差异可以构造方向权函数，在克希霍 夫积分偏移中加入方向权函数即可完成偏移中的假象压制，结果如图9～图13 所示。为了验证偏移算法有效性，首先利用真速度模型进行克希霍夫叠前深度 偏移成像。图9为利用真速度模型A进行极化偏移在深度z＝25m处XOY截面 的成像结果，虚线表示真实介质分界面与z＝25m处的水平切面的交线，点划 线表示真实介质分界面关于观测测线对称的面与z＝25m处的水平切面的交线， 垂直的线表示地震观测测线，可以看到，成像区域聚集真实界面(虚线)同一侧，而在关于观测测线对称处(点划线)的成像假象几乎被完全压制；
然后利用偏移速度模型进行克希霍夫叠前深度偏移成像，图10～12的(a) 图为非极化偏移，K＝1方向权因子，图10～图12的(b)图方向权因子K＝cosθ， 指数κ＝2，图10～图12的(c)图中方向权因子K＝f5，指数κ＝1；图13(a1) —(a3)对应的为图10、图11、图12(a)图的偏移成像结果切片图，图13 (b1)—(b3)对应的为图10、图11、图12(b)图的偏移成像结果切片图 图13(c1)—(c3)对应的为图10、图11、图12(c)图的偏移成像结果切 片图，切片深度z＝25m；通过观察图10～12所示成像结果并对比同一地震分 量记录不同权函数的克希霍夫积分法叠前深度偏移成像结果可知：(1)当方向 权因子K＝1即非极化偏移时，与二维绕射扫描成像结果类似，介质分界面关 于地震观测测线对称处存在有假象；(2)利用射线追踪计算的波场传播方向单 位矢量与反射纵波主极化方向单位矢量的内积，即当方向权因子选为K＝cosθ， 指数κ＝2时，整体成像幅值会被削弱，假象被削弱的程度要大于真像区域， 即成像结果中的假象被抑制，但从图中可以看出，上述方向权因子并不能对假 象形成很好的选择性，假象仍然以较大幅值存在于成像结果中，这对后续对偏 移成像的结果进行地质解释造成了困难；(3)当方向权因子为K＝f5，κ＝1时， 对于某一分量的偏移结果，若仍以相同比例尺进行显示，在前述模型情形下， 已看不到明显的偏移假象存在，真实成像目标区域也进一步被缩小。
利用加入了方向权函数的克希霍夫积分进行极化偏移，具体做法如下：
Figure BDA0002281538920000251
其中，κ表示调节因子，用于调节权因子在成像过程中的压制效果，K(p,d)表 示方向权函数，
Figure BDA0002281538920000252
表示权因子，为积分过程中各道观测值对最终各成像 点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S 表示积分曲面，由地面和一个向地下空间无限延伸的部分球面组成，同时，在 无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。
地震勘探中接收到的地震波场看作是地下无数二次绕射点源发出的波在 接收点处的叠加结果。若已知地下的速度结构，则可以将偏移目标区域离散成 一定数量的离散点，通过射线追踪的方法获得从震源发出的波经过每个绕射点 到各检波器的波场绕射时间，进而根据各离散点的绕射时间从地震记录上选取 振幅值进行叠加，即可以得到该点的偏移成像结果。对所有离散点进行上述过 程，则可以完成整个目标区域的偏移成像工作。
方向权函数的具体形式如下：
Figure RE-GDA0002411727260000253
其中，Θ的具体形式如下：
Figure BDA0002281538920000254
θ表示射线追踪计算的绕射射线在检波器处的出射方向和根据绕射时间 对地震数据进行极化分析所得的纵波极化方向两方向之间的方向夹角。
根据图14所示，竖直断层模型B为与测线方向成大角度斜交的断层，断 层两侧为岩性不同的两种介质，断层界面与YOZ平面所成角度角约为63.43°， 用于接收的检波器数量为121道，采样总时长为200ms，其余参数与模型A 的参数一样，同时，模型B的成像方法和过程也与模型A的一致。
具体实施方式为，首先将三分量纵波地震记录进行波场分离，结果如图 15所示，图15中的同相轴代表被分离出来的纵波，图15(a)、(b)、(c)依 次对应X、Y、Z三分量地震记录。然后将图15所示的分离纵波后的地震记录 进行速度扫描得到三维速度谱，结果如图16所示。其中，图17为模型B三维 速度扫描结果ZOX平面切片图，图17(a)～图17(f)依次对应于界面法线 方向与Y轴夹角范围45°～70°，角度间隔5°；
然后利用协方差矩阵法进行极化分析，结果如图18所示，求得极化程度 参数值和主特征向量的方向余弦值。结合图15所示纵波地震记录及图18所示 极化分析结果进行分析可知，纵波信号极化程度大，且主要在XOY平面内极 化。
然后利用地震的极化信息来求取模型B的反射界面参数，如表4所示。
表4模型B极化分析计算得到的反射界面参数
Figure BDA0002281538920000261
通过求取空间反射界面与给定模型B空间棱边的交点，得到表5所示的四 个交点坐标。
表5模型A空间棱边的交点
Figure BDA0002281538920000271
通过表4与表5上述求得的反射界面几何参数与图16中三维速度谱结合， 即可得到经速度分析后的偏移速度模型。
已知偏移速度模型后，利用射线追踪方法可以得到旅行时与绕射波出射方 向，利用绕射波出射方向和主偏振方向之间的差异可以构造方向权函数，在克 希霍夫积分偏移中加入方向权函数即可完成偏移中的假象压制，具体做法与模 型A的成像过程相同，结果如图19～图22所示。图19为利用真速度模型B 进行非极化偏移的结果，方向权因子K＝1，图19(a)～(c)依次对应反射 纵波X、Y、Z分量地震记录的成像结果；图20为利用真速度模型B进行极化 偏移的结果，方向权因子K＝f5，指数κ＝1，图20(a)～(c)依次对应反射 纵波X、Y、Z分量地震记录的成像结果；图21为利用速度分析结果构建的速 度模型B进行极化偏移的结果，方向权因子K＝f5，指数κ＝1，图21(a)～ (c)依次对应反射纵波X、Y、Z分量地震记录的成像结果；
图22为图19、图20、图21的偏移成像结果切片图，切片深度z＝25m， 虚线为反射界面所在位置。从图22(a1)、(a2)、(a3)可知，利用真速度模 型进行偏移，偏移结果能准确表示反射界面所在位置，但是存在偏移假象，且 假象与正确的像能量相当，若不利用极化信息来对波场方向进行分辨，则无法 确定真正成像位置的方位；图22(b1)、(b2)、(b3)表明，利用极化偏移能 较好压制偏移假象，但是，在压制假象的同时，引入了噪声，对成像结果造成 了干扰；图22(c1)、(c2)、(c3)所示为利用速度分析结果构建的速度模型 进行偏移结果，图中也显示了极化偏移对偏移假象的压制效果，但是，由于速 度模型不够准确，造成成像位置的准确程度不如图22(a1)、(a2)、(a3)中 所示的结果。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于 此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到 变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应 所述以权利要求的保护范围为准。

Claims (7)

1.一种岩巷超前探测三维偏移成像方法，其特征在于：包含如下步骤：
S1、获取三分量地震观测数据，将地震观测数据中的直达波剔除，并分离纵波；
S2、利用三维速度分析公式，对预处理后的地震观测数据进行速度扫描分析，得到三维速度谱，以估计地下介质地震波速度，再利用初至纵波的极化信息确定反射界面的几何参数，利用三维速度谱和反射界面的几何参数确定速度模型；
S3、确定时窗计算长度W，利用协方差矩阵特征分析法，求取三分量地震记录的波场极化信息，得到对应绕射时刻的主偏振方向；
S4、根据步骤S2得到的速度模型，利用射线追踪的方法，计算获得：从震源发出的波经过每个绕射点到各检波器的波场绕射时间和绕射射线出射方向；
S5、利用绕射射线出射方向与主偏振方向之间的差异构造方向权函数；
S6、在均匀各向同性介质中的克希霍夫积分法偏移的过程中，加入方向权函数项，得到极化深度偏移成像结果。
2.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S2中：三维速度分析公式ESpec，具体形式如下：
Figure FDA0002281538910000011
其中，i表示道序号，j表示采样点在窗口内的自然排列顺序序号，
Figure FDA0002281538910000012
表示各道的反射时间，M表示速度扫描计算时窗长度，t0表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向单位向量的Y分量；
首先计算各道的反射时间ti，再依据公式计算与t0、v、m对应的速度谱值，对上述三个变量选定不同值，重复以上步骤即可以计算地震剖面的三维速度谱，通过三维速度谱，估计反射界面到炮点之间的介质的速度大小。
3.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S2中：利用初至纵波的极化信息确定反射界面几何参数的公式，具体形式如下：
Figure FDA0002281538910000021
其中，N表示假设已知的直线束Li的数量，(i＝1，2，3，…，N)，Pi(xi,yi,zi)表示每一直线上的某一已知坐标的点，Q(x,y,z)表示空间内的任意点，vi＝(li,mi,ni)表示各直线的单位方向向量，I表示N阶单位矩阵；
具体为利用极化特征，确定地下界面的法线方向及其位置，先对波的类型做出判别，再对纵波进行极化分析，以获取波场在不同位置的传播方向，然后通过波场在不同位置的传播方向进行反向延长，求取交点，通过解如上线性方程，求得反射射线的反向交点Q的坐标，交点与炮点关于界面镜像对称，由交点与炮点确定的直线的单位方向的向量v为界面法线方向，所述界面法线方向及交点与炮点构成线段的中点M(xM,yM,zM)，确定反射界面的解析表达式；其中，点(xp,yp,zp)为空间平面上不同于点M的任一点；
Figure FDA0002281538910000022
在计算过程中给出三维速度模型空间的大小限制，利用空间平面与速度模型空间的各棱边直线的交点来估计反射界面的具体位置，从而构建初始速度模型。
4.根据权利要求3所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述三维速度模型空间包含反射界面内所有的目标体。
5.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S3中：
所述协方差矩阵，具体形式如下：
Figure FDA0002281538910000031
其中，W表示时窗计算长度，u(x,y,z)表示三分量地震矢量波场，
Figure FDA0002281538910000032
表示方差，
Figure FDA0002281538910000033
表示协方差，
Figure FDA0002281538910000034
表示数学期望，
Figure FDA0002281538910000035
表示采样点数，N表示每道的总采样点数，ΔT表示时间采样间隔；
协方差矩阵为一对称矩阵，求取协方差矩阵的特征值λ1、λ2、λ3，且限定λ1≥λ2≥λ3，及其对应的单位特征向量为p1、p2、p3，协方差矩阵的三个特征值分别对应极化椭圆的三个轴的长度，最大特征值λ1对应极化椭圆的主轴，相应的主特征向量p1表示质点椭圆偏振的主方向，其余同理；
该时刻波的偏振类型用偏振系数γ(t)来衡量：
Figure FDA0002281538910000041
若γ值趋于1，主特征值λ1会远大于其它次特征值，并表现为极化椭球趋于细长，表示质点的空间振动过程趋于线性极化，此时主特征向量所对应的空间方向即为质点振动所在的方向；若γ值趋于0，主特征值会和次特征值相当，极化椭球表现为趋于球形，表示质点在空间振动时的非线性极化越强，此时质点振动方向没有一个优势方向，通过窗口在地震记录上的滑动，求取整个剖面上不同采样点处，不同时刻的波的极化特征信息。
6.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S5中，所述方向权函数的具体形式如下：
Figure RE-FDA0002371817970000042
其中，Θ的具体形式如下：
Figure RE-FDA0002371817970000043
θ表示射线追踪计算的绕射射线在检波器处的出射方向和根据绕射时间对地震数据进行极化分析所得的纵波极化方向两方向之间的方向夹角。
7.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S6中：
均匀各向同性介质中的克希霍夫积分解的具体形式为：
Figure FDA0002281538910000044
其中，第一项由波场的垂直梯度决定，第二项与波场函数值的延迟位有关，叫做近场源项，它以
Figure FDA0002281538910000051
衰减，这两项通常在克希霍夫积分法偏移中被忽略，而只使用第三项远场源项，第三项可以用波场观测数据项和权因子项W来表征：
Figure FDA0002281538910000052
其中，
Figure FDA0002281538910000053
表示权因子，即在积分过程中，各道观测值对最终各成像点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S表示积分曲面；
加入了方向权函数的克希霍夫积分偏移公式，具体形式如下：
Figure FDA0002281538910000054
其中，κ表示调节因子，用于调节权因子在成像过程中的压制效果，K(p,d)表示方向权函数，
Figure FDA0002281538910000055
表示权因子，为积分过程中各道观测值对最终各成像点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S表示积分曲面，由地面和一个向地下空间无限延伸的部分球面组成，同时，在无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。
CN201911143406.4A 2019-11-20 2019-11-20 一种岩巷超前探测三维偏移成像方法 Pending CN111103620A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911143406.4A CN111103620A (zh) 2019-11-20 2019-11-20 一种岩巷超前探测三维偏移成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911143406.4A CN111103620A (zh) 2019-11-20 2019-11-20 一种岩巷超前探测三维偏移成像方法

Publications (1)

Publication Number Publication Date
CN111103620A true CN111103620A (zh) 2020-05-05

Family

ID=70421503

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911143406.4A Pending CN111103620A (zh) 2019-11-20 2019-11-20 一种岩巷超前探测三维偏移成像方法

Country Status (1)

Country Link
CN (1) CN111103620A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112925010A (zh) * 2021-01-26 2021-06-08 云南航天工程物探检测股份有限公司 一种高精度相控阵弹性波隧道三维地质超前预报方法
CN113805232A (zh) * 2020-06-17 2021-12-17 中国石油化工股份有限公司 浅层地表的品质因子的估计方法、系统及存储介质
CN115631834A (zh) * 2022-12-23 2023-01-20 北京馐圣医学技术有限公司 一种血液标量波图像的通信方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8428880B2 (en) * 2007-12-14 2013-04-23 Shell Oil Company Method of processing data obtained from seismic prospecting
CN102495426B (zh) * 2011-12-02 2014-10-22 北京中科联衡科技有限公司 一种克希霍夫积分地震成像方法
CN105974472A (zh) * 2016-05-13 2016-09-28 中国矿业大学 一种基于反射信号的巷道超前探测速度建模方法
CN106772557A (zh) * 2016-11-29 2017-05-31 北京中矿大地地球探测工程技术有限公司 利用随掘信号探测煤矿掘进巷道各方向地质构造的方法
CN104166164B (zh) * 2014-08-08 2017-08-29 山东科技大学 煤巷掘进地质构造三分量多波反射三维地震超前探测方法
CN106443765B (zh) * 2016-08-30 2018-08-28 安徽惠洲地质安全研究院股份有限公司 基于多分量观测系统的城市工程地震探测综合成像方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8428880B2 (en) * 2007-12-14 2013-04-23 Shell Oil Company Method of processing data obtained from seismic prospecting
CN102495426B (zh) * 2011-12-02 2014-10-22 北京中科联衡科技有限公司 一种克希霍夫积分地震成像方法
CN104166164B (zh) * 2014-08-08 2017-08-29 山东科技大学 煤巷掘进地质构造三分量多波反射三维地震超前探测方法
CN105974472A (zh) * 2016-05-13 2016-09-28 中国矿业大学 一种基于反射信号的巷道超前探测速度建模方法
CN106443765B (zh) * 2016-08-30 2018-08-28 安徽惠洲地质安全研究院股份有限公司 基于多分量观测系统的城市工程地震探测综合成像方法
CN106772557A (zh) * 2016-11-29 2017-05-31 北京中矿大地地球探测工程技术有限公司 利用随掘信号探测煤矿掘进巷道各方向地质构造的方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JORG SCHLEICHER，MARTIN TYGEL，AND PETER HUBRAL: "3-D true-amplitude finite-offset migration", 《GEOPHYSICS》 *
刘盛东，章俊，李纯阳等: "矿井多波多分量地震方法与试验", 《煤炭学报》 *
王勃: "矿井地震全空间极化偏移成像技术研究", 《中国优秀博硕士学位论文全文数据库(博士) 基础科学辑》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113805232A (zh) * 2020-06-17 2021-12-17 中国石油化工股份有限公司 浅层地表的品质因子的估计方法、系统及存储介质
CN113805232B (zh) * 2020-06-17 2024-04-09 中国石油化工股份有限公司 浅层地表的品质因子的估计方法、系统及存储介质
CN112925010A (zh) * 2021-01-26 2021-06-08 云南航天工程物探检测股份有限公司 一种高精度相控阵弹性波隧道三维地质超前预报方法
CN115631834A (zh) * 2022-12-23 2023-01-20 北京馐圣医学技术有限公司 一种血液标量波图像的通信方法
CN115631834B (zh) * 2022-12-23 2023-05-05 北京馐圣医学技术有限公司 一种血液标量波图像的通信方法

Similar Documents

Publication Publication Date Title
Bishop et al. Tomographic determination of velocity and depth in laterally varying media
US7376539B2 (en) Method for simulating local prestack depth migrated seismic images
EP2673662B1 (en) Method of analyzing seismic data
US6826484B2 (en) 3D prestack time migration method
US8120991B2 (en) System and method for full azimuth angle domain imaging in reduced dimensional coordinate systems
Nishizawa et al. Seismic structure of rifting in the Okinawa Trough, an active backarc basin of the Ryukyu (Nansei-Shoto) island arc–trench system
US8289809B2 (en) Common reflection azimuth migration
US20100131205A1 (en) Method for identifying and analyzing faults/fractures using reflected and diffracted waves
Liu et al. A new 3D observation system designed for a seismic ahead prospecting method in tunneling
US8731838B2 (en) Fresnel zone fat ray tomography
CA2636250A1 (en) Traveltime calculation in three dimensional transversely isotropic (3d tti) media by the fast marching method
CN111103620A (zh) 一种岩巷超前探测三维偏移成像方法
Bleistein et al. Migration/inversion: think image point coordinates, process in acquisition surface coordinates
Share et al. Structural properties of the San Jacinto fault zone at Blackburn Saddle from seismic data of a dense linear array
Pica Fast and accurate finite-difference solutions of the 3D eikonal equation parametrized in celerity.
Chen et al. A compact program for 3D passive seismic source‐location imaging
Gong et al. Combined migration velocity model-building and its application in tunnel seismic prediction
CN109143362B (zh) 基于共散射角道集的散射波分离方法
US10379245B2 (en) Method and system for efficient extrapolation of a combined source-and-receiver wavefield
US20140297193A1 (en) Seismic methods and systems employing shallow shear-wave splitting analysis using receiver functions
US9791580B2 (en) Methods and systems to separate wavefields using pressure wavefield data
Rochlin et al. Imaging shallow-linear diffractors using 3D-converted wave data
Müller On the nature of scattering from isolated perturbations in elastic media and the consequences for processing of seismic data
Rochlin et al. Detection of subsurface lineaments using edge diffraction
Granger et al. Preliminary evaluation of azimuthal anisotropy over the Valhall Field using C-wave data

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: 20200505