一种岩巷超前探测三维偏移成像方法
阅读：46发布：2021-06-11
专利汇可以提供一种岩巷超前探测三维偏移成像方法专利检索，专利查询，专利分析的服务。并且本 发明 提供了一种岩巷超前探测三维偏移成像方法，通过获取三分量 地震 观测数据，利用三维速度分析公式，得到三维速度谱，获得从 震源 发出的波经过每个绕射点到各 检波器 的波场绕射时间和绕射射线出射方向，利用绕射射线出射方向与主偏振方向之间的差异构造方向权函数；在均匀各向同性介质中的克希霍夫积分法偏移的过程中，加入方向权函数项，得到极化深度偏移成像结果。本发明采用三维速度扫描分析来估算速度，能够更加准确的估计隧道介质速度，并利用初至纵波的极化信息来确定反射界面的几何参数，两者结合能够准确的确定初始速度模型。在克希霍夫积分偏移中引入了方向权函数，能够有效的压制隧道中克希霍夫积分成像方法中产生的偏移假象。，下面是一种岩巷超前探测三维偏移成像方法专利的具体信息内容。
1.一种岩巷超前探测三维偏移成像方法，其特征在于：包含如下步骤：
S1、获取三分量地震观测数据，将地震观测数据中的直达波剔除，并分离纵波；
S2、利用三维速度分析公式，对预处理后的地震观测数据进行速度扫描分析，得到三维速度谱，以估计地下介质地震波速度，再利用初至纵波的极化信息确定反射界面的几何参数，利用三维速度谱和反射界面的几何参数确定速度模型；
S3、确定时窗计算长度W，利用协方差矩阵特征分析法，求取三分量地震记录的波场极化信息，得到对应绕射时刻的主偏振方向；
S4、根据步骤S2得到的速度模型，利用射线追踪的方法，计算获得：从震源发出的波经过每个绕射点到各检波器的波场绕射时间和绕射射线出射方向；
S5、利用绕射射线出射方向与主偏振方向之间的差异构造方向权函数；
S6、在均匀各向同性介质中的克希霍夫积分法偏移的过程中，加入方向权函数项，得到极化深度偏移成像结果。
2.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S2中：三维速度分析公式ESpec，具体形式如下：
其中，i表示道序号，j表示采样点在窗口内的自然排列顺序序号，
表示各道的反射时间，M表示速度扫描计算时窗长度，t0表
示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向单位向量的Y分量；
首先计算各道的反射时间ti，再依据公式计算与t0、v、m对应的速度谱值，对上述三个变量选定不同值，重复以上步骤即可以计算地震剖面的三维速度谱，通过三维速度谱，估计反射界面到炮点之间的介质的速度大小。
3.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S2中：利用初至纵波的极化信息确定反射界面几何参数的公式，具体形式如下：
其中，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的任一点；
在计算过程中给出三维速度模型空间的大小限制，利用空间平面与速度模型空间的各棱边直线的交点来估计反射界面的具体位置，从而构建初始速度模型。
4.根据权利要求3所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述三维速度模型空间包含反射界面内所有的目标体。
5.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S3中：
所述协方差矩阵，具体形式如下：
其中，W表示时窗计算长度，u(x,y,z)表示三分量地震矢量波场， 表
示方差， 表示协方差，
表示数学期望， 表示采样点数，N表示每道的总采样点数，ΔT表示时间采
样间隔；
协方差矩阵为一对称矩阵，求取协方差矩阵的特征值λ1、λ2、λ3，且限定λ1≥λ2≥λ3，及其对应的单位特征向量为p1、p2、p3，协方差矩阵的三个特征值分别对应极化椭圆的三个轴的长度，最大特征值λ1对应极化椭圆的主轴，相应的主特征向量p1表示质点椭圆偏振的主方向，其余同理；
该时刻波的偏振类型用偏振系数γ(t)来衡量：
若γ值趋于1，主特征值λ1会远大于其它次特征值，并表现为极化椭球趋于细长，表示质点的空间振动过程趋于线性极化，此时主特征向量所对应的空间方向即为质点振动所在的方向；若γ值趋于0，主特征值会和次特征值相当，极化椭球表现为趋于球形，表示质点在空间振动时的非线性极化越强，此时质点振动方向没有一个优势方向，通过窗口在地震记录上的滑动，求取整个剖面上不同采样点处，不同时刻的波的极化特征信息。
6.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S5中，所述方向权函数的具体形式如下：
其中，Θ的具体形式如下：
θ表示射线追踪计算的绕射射线在检波器处的出射方向和根据绕射时间对地震数据进行极化分析所得的纵波极化方向两方向之间的方向夹角。
7.根据权利要求1所述的一种岩巷超前探测三维偏移成像方法，其特征在于：所述步骤S6中：
均匀各向同性介质中的克希霍夫积分解的具体形式为：
其中，第一项由波场的垂直梯度决定，第二项与波场函数值的延迟位有关，叫做近场源项，它以 衰减，这两项通常在克希霍夫积分法偏移中被忽略，而只使用第三项远场源项，第三项可以用波场观测数据项和权因子项W来表征：
其中， 表示权因子，即在积分过程中，各道观测值对最终各成像点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S表示积分曲面；
加入了方向权函数的克希霍夫积分偏移公式，具体形式如下：
其中，κ表示调节因子，用于调节权因子在成像过程中的压制效果，K(p,d)表示方向权函数， 表示权因子，为积分过程中各道观测值对最终各成像点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S表示积分曲面，由地面和一个向地下空间无限延伸的部分球面组成，同时，在无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。
一种岩巷超前探测三维偏移成像方法
技术领域
背景技术
[0002] 随着我国交通、水利水电事业的快速发展，大规模的铁路、公路和水电工 程建设数量不断增加，需要修建大量各种类型的隧道，在隧道施工建设中，掌 子面前方常隐含有不利地质构造对象，在开挖揭露过程中容易诱发地质灾害导 致安全事故，给工程建设造成重大损失，为了准确确定隧道掘进前方的异常地 质构造，就要在隧道的超前地震预报数据处理过程中，得到精准的偏移成像结 果。
[0003] 目前，岩巷超前地震探测中较常用的偏移方法是克希霍夫积分偏移，将射 线追踪方法与Kirchhoff积分偏移结合起来，可以实现叠前Kirchhoff深度偏 移来进行隧道地质预报；通过分析隧道环境下的地震波传播的时距关系，可以 通过插值的方法使用不同精度的网格剖分，以提高运算速度。由于隧道施工过 程中地下空间地质构造复杂且容易出现不连续的情况，从物理地震学的角度来 看，观测到的波场信息可以认为是地下所有点产生的绕射在观测点处叠加的结 果，因此，现有技术也有将散射成像技术应用于隧道地震超前地质预报中，当 观测系统的横向展布不够时，按照散射原理偏移，偏移结果中会呈现往外扩展 的圆弧，可以通过方向扫描的方法来改进这一缺陷。
[0004] 现有的各种隧道地震超前预报技术，仍存在一些问题：首先，在隧道地震 超前探测中，采用直达波估计得到的速度只代表隧道侧壁围岩表层介质的速 度，并不能反映波场在反射界面到测线之间的介质中的地震波传播速度。传统 速度分析方法以水平层状介质模型假设为前提，不适用于隧道环境下的速度估 计；其次，由于隧道施工环境限制，地震观测测线多沿隧道轴线布置于地下空 间，当沿用地面地震勘探中的偏移方法对数据进行偏移处理时，会产生关于测 线对称的偏移假象。
发明内容
[0005] 根据上述阐述，本发明的目的在于提供一种岩巷超前探测三维偏移成像方 法，解决偏移成像结果中存在假象的问题。
[0006] 本发明主要原理是利用波场的极化信息，引入了方向权函数来压制克希霍 夫积分偏移的假象，其具体提供的技术方案为：
[0007] 一种岩巷超前探测三维偏移成像方法，包含如下步骤：
[0008] S1、获取三分量地震观测数据，将地震观测数据中的直达波剔除，并分离 纵波；
[0009] S2、利用三维速度分析公式，对预处理后的地震观测数据进行速度扫描分 析，得到三维速度谱，以估计地下介质地震波速度，再利用初至纵波的极化信 息确定反射界面的几何参数，利用三维速度谱和反射界面的几何参数确定速度 模型；
[0010] S3、确定时窗计算长度W，利用协方差矩阵特征分析法，求取三分量地震 记录的波场极化信息，得到对应绕射时刻的主偏振方向；
[0012] S5、利用绕射射线出射方向与主偏振方向之间的差异构造方向权函数；
[0013] S6、在均匀各向同性介质中的克希霍夫积分法偏移的过程中，加入方向权 函数项，得到极化深度偏移成像结果。
[0014] 上述技术方案中，所述步骤S2中：三维速度分析公式ESpec，具体形式如 下：
[0015]
[0016] 其中，i表示道序号，j表示采样点在窗口内的自然排列顺序序号， 表示各道的反射时间，M表示速度扫描计算时窗 长度，t0
表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向 单位向量的Y分量；
表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向 单位向量的Y分量；
[0017] 首先计算各道的反射时间ti，再依据公式计算与t0、v、m对应的速度谱值， 对上述三个变量选定不同值，重复以上步骤即可以计算地震剖面的三维速度 谱，通过三维速度谱，估计反射界面到炮点之间的介质的速度大小。
[0018] 上述技术方案中，所述步骤S2中：利用初至纵波的极化信息确定反射界 面几何参数的公式，具体形式如下：
[0019]
[0020] 其中，N表示假设已知的直线束Li的数量，(i＝1，2，3，…，N)，Pi(xi,yi,zi) 表示每一直线上的某一已知坐标的点，Q(x,y,z)表示空间内的任意点， vi＝(li,mi,ni)表示各直线的单位方向向量，I表示N阶单位矩阵；
[0021] 具体为利用极化特征，确定地下界面的法线方向及其位置，先对波的类型 做出判别，再对纵波进行极化分析，以获取波场在不同位置的传播方向，然后 通过波场在不同位置的传播方向进行反向延长，求取交点，通过解如上线性方 程，求得反射射线的反向交点Q的坐标，交点与炮点关于界面镜像对称，由 交点与炮点确定的直线的单位方向的向量v为界面法线方向，所述界面法线方 向及交点与炮点构成线段的中点M(xM,yM,zM)，确定反射界面的解析表达式； 其中，点(xp,yp,zp)为空间平面上不同于点M的任一点。
[0022]
[0023] 在计算过程中给出三维速度模型空间的大小限制，利用空间平面与速度模 型空间的各棱边直线的交点来估计反射界面的具体位置，从而构建初始速度模 型。
[0024] 上述技术方案中，所述三维速度模型空间包含反射界面内所有的目标体。
[0025] 上述技术方案中，所述步骤S3中：
[0026] 所述协方差矩阵，具体形式如下：
[0027]
[0028] 其中，W表示时窗计算长度，u(x,y,z)表示三分量地震矢量波场， 表示方差， 表示协
方差， 表示数学期望， 表示采 样点
数，N表示每道的总采样点数，ΔT表示时间采样间隔；
方差， 表示数学期望， 表示采 样点
数，N表示每道的总采样点数，ΔT表示时间采样间隔；
[0029] 协方差矩阵为一对称矩阵，求取协方差矩阵的特征值λ1、λ2、λ3，且限定 λ1≥λ2≥λ3，及其对应的单位特征向量为p1、p2、p3，协方差矩阵的三个特征值 分别对应极化椭圆的三个轴的长度，最大特征值λ1对应极化椭圆的主轴，相应 的主特征向量p1表示质点椭圆偏振的主方向，其余同理；
[0030] 该时刻波的偏振类型用偏振系数γ(t)来衡量：
[0031]
[0032] 若γ值趋于1，主特征值λ1会远大于其它次特征值，并表现为极化椭球 趋于细长，表示质点的空间振动过程趋于线性极化，此时主特征向量所对应的 空间方向即为质点振动所在的方向；若γ值趋于0，主特征值会和次特征值相 当，极化椭球表现为趋于球形，表示质点在空间振动时的非线性极化越强，此 时质点振动方向没有一个优势方向，通过窗口在地震记录上的滑动，求取整个 剖面上不同采样点处，不同时刻的波的极化特征信息。
[0033] 上述技术方案中，所述步骤S5中，所述方向权函数的具体形式如下：
[0034]
[0035] 其中，Θ的具体形式如下：
[0036]
[0037] θ表示射线追踪计算的绕射射线在检波器处的出射方向和根据绕射时间 对地震数据进行极化分析所得的纵波极化方向两方向之间的方向夹角。
[0038] 上述技术方案中，所述步骤S6中：
[0039] 均匀各向同性介质中的克希霍夫积分解的具体形式为：
[0040]
[0041] 其中，第一项由波场的垂直梯度决定，第二项与波场函数值的延迟位有关，叫 做近场源项，它以 衰减，这两项通常在克希霍夫积分法偏移中被忽略，而 只使用第三项远场源项，第三项可以用波场观测数据项和权因子项W来表征：
[0042]
[0043] 其中， 表示权因子，即在积分过程中，各道观测值对最终各成像点幅 值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S表示 积分曲面；
[0044] 加入了方向权函数的克希霍夫积分偏移公式，具体形式如下：
[0045]
[0046] 其中，κ表示调节因子，用于调节权因子在成像过程中的压制效果，K(p,d)表 示方向权函数， 表示权因子，为积分过程中各道观测值对最终各成像 点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S 表示积分曲面，由地面和一个向地下空间无限延伸的部分球面组成，同时，在 无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。
[0047] 本发明与现有的岩巷克希霍夫积分偏移成像技术相比，优点为：采用 三维速度扫描分析来估算速度，能够更加准确的估计隧道介质速度，并利 用初至纵波的极化信息来确定反射界面的几何参数，两者结合能够准确的 确定初始速度模型。在克希霍夫积分偏移中引入了方向权函数，能够有效 的压制隧道中克希霍夫积分成像方法中产生的偏移假象。附图说明
[0048] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施 例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述 中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付 出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0049] 图1为本发明实施例的流程图；
[0050] 图2为虚实反射界面的反射纵波传播方向与反射纵波极化方向的差异示 意图；
[0051] 图3为与测线方向成小角度斜交的竖直断层模型A的示意图；
[0052] 图4为模型A经波场分离后的三分量纵波地震记录示意图；
[0053] 图5为模型A纯纵波地震记录经速度扫描所得的三维速度谱示意图；
[0054] 图6为模型A三维速度扫描结果ZOX平面切片图，其中(a)～(f)依次 对应于界面法线方向与Y轴夹角范围45°～70°，角度间隔5°；
[0055] 图7为模型A极化分析结果图，其中(a)表示不同时刻的极化程度参数 值，(b)～(d)表示主特征向量沿X、Y、Z三轴的方向余弦值；
[0056] 图8为三维单倾斜界面地震观测系统
[0057] 图9为利用真速度模型A进行极化偏移在深度z＝25m处XOY截面的成像结 果；
[0058] 图10为利用速度分析结果构建的速度模型A，得到的X分量地震记录采 用不同方向权函数时的偏移成像结果；
[0059] 图11为利用速度分析结果构建的速度模型A，得到的Y分量地震记录采 用不同方向权函数时的偏移成像结果；
[0060] 图12为利用速度分析结果构建的速度模型A，得到的Z分量地震记录采 用不同方向权函数时的偏移成像结果；
[0061] 图13为图9、图10、图11中，切片深度z＝25m的偏移成像结果切片图；
[0062] 图14为与测线方向成大角度斜交的竖直断层模型B的示意图；
[0063] 图15为模型B经波场分离后的三分量纵波地震记录示意图；
[0064] 图16为模型B纯纵波地震记录经速度扫描所得的三维速度谱示意图；
[0065] 图17为模型B三维速度扫描结果ZOX平面切片图，其中(a)～(f)依 次对应于界面法线方向与Y轴夹角范围45°～70°，角度间隔5°；
[0066] 图18为模型B极化分析结果图，其中(a)表示不同时刻的极化程度参数 值，(b)～(d)表示主特征向量沿X、Y、Z三轴的方向余弦值；
[0067] 图19为利用真速度模型B进行非极化偏移的结果，方向权因子K＝1，图 (a)～(c)依次对应反射纵波X、Y、Z分量地震记录的成像结果；
[0068] 图20为利用真速度模型B进行极化偏移的结果，方向权因子K＝f5，指数 κ＝1，(a)～(c)依次对应反射纵波X、Y、Z分量地震记录的成像结果；
[0069] 图21为利用速度分析结果构建的速度模型B进行极化偏移的结果，方向 权因子K＝f5，指数κ＝1，(a)～(c)依次对应反射纵波X、Y、Z分量地震记 录的成像结果；
[0070] 图22为图19、图20、图21的偏移成像结果切片图，切片深度z＝25m。
具体实施方式
[0071] 下面将结合本发明的附图，对本发明的技术方案进行清楚、完整地描述， 显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基 于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获 得的所有其他实施例，都属于本发明保护的范围。
[0072] 根据图1所示，本发明提供的一种岩巷超前探测三维偏移成像方法，包含 如下步骤：
[0073] S1、获取三分量地震观测数据，将地震观测数据中的直达波剔除，并分离 纵波；
[0074] S2、三维速度扫描分析和反射界面几何参数估计，利用三维速度分析公式， 对预处理后的地震观测数据进行速度扫描分析，得到三维速度谱，以估计地下 介质地震波速度，再利用初至纵波的极化信息确定反射界面的几何参数，利用 三维速度谱和反射界面的几何参数确定速度模型；
[0075] 三维速度分析公式ESpec，具体形式如下：
[0076]
[0077] 其中，i表示道序号，j表示采样点在窗口内的自然排列顺序序号， 表示各道的反射时间，M表示速度扫描计算时窗 长度，t0
表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向 单位向量的Y分量；
表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向 单位向量的Y分量；
[0078] 首先计算各道的反射时间ti，再依据公式计算与t0、v、m对应的速度谱值， 对上述三个变量选定不同值，重复以上步骤即可以计算地震剖面的三维速度 谱，通过三维速度谱，估计反射界面到炮点之间的介质的速度大小。
[0079] 利用初至纵波的极化信息确定反射界面几何参数的公式，具体形式如下：
[0080]
[0081] 其中，N表示假设已知的直线束Li的数量，(i＝1，2，3，…，N)，Pi(xi,yi,zi)表 示每一直线上的某一已知坐标的点，Q(x,y,z)表示空间内的任意点， vi＝(li,mi,ni)表示各直线的单位方向向量，I表示N阶单位矩阵；
[0082] 具体为利用极化特征，确定地下界面的法线方向及其位置，先对波的类型 做出判别，再对纵波进行极化分析，以获取波场在不同位置的传播方向，然后 通过波场在不同位置的传播方向进行反向延长，求取交点，通过解如上线性方 程，求得反射射线的反向交点Q的坐标，交点与炮点关于界面镜像对称，由 交点与炮点确定的直线的单位方向的向量v为界面法线方向，所述界面法线方 向及交点与炮点构成线段的中点M(xM,yM,zM)，确定反射界面的解析表达式； 其中，点(xp,yp,zp)为空间平面上不同于点M的任一点。
[0083]
[0084] 若能在计算过程中给出三维速度模型空间(或偏移目标区域)的大小限制， 利用空间平面与速度模型空间的各棱边直线的交点来估计反射界面的具体位 置，从而构建初始速度模型，三维速度模型空间应当包含反射界面内所有的目 标体，例如野外工作时，假设地下目标体形状为一个正方体，大小为 100m*100m*100m，此时速度模型空间只要大于该目标体的大小即可，例如可 以设置为120m*120m*120m，总而言之，只要模型空间大小包括所有需要被 测量的目标体即可。
[0085] S3、求取地震记录的波场极化信息，得到对应绕射时刻的主偏振方向，即 确定时窗计算长度W，利用协方差矩阵特征分析法，求取三分量地震记录的波 场极化信息，得到对应绕射时刻的主偏振方向；
[0086] 所述协方差矩阵，具体形式如下：
[0087]
[0088] 其中，W表示时窗计算长度，u(x,y,z)表示三分量地震矢量波场， 表示方差， 表示协
方差， 表示数学期望， 表示采 样点
数，N表示每道的总采样点数，ΔT表示时间采样间隔；
方差， 表示数学期望， 表示采 样点
数，N表示每道的总采样点数，ΔT表示时间采样间隔；
[0089] 协方差矩阵为一对称矩阵，求取协方差矩阵的特征值λ1、λ2、λ3，且限定 λ1≥λ2≥λ3，及其对应的单位特征向量为p1、p2、p3，协方差矩阵的三个特征值 分别对应极化椭圆的三个轴的长度，最大特征值λ1对应极化椭圆的主轴，相应 的主特征向量p1表示质点椭圆偏振的主方向，其余同理；
[0090] 该时刻波的偏振类型用偏振系数γ(t)来衡量：
[0091]
[0092] 若γ值趋于1，主特征值λ1会远大于其它次特征值，并表现为极化椭球 趋于细长，表示质点的空间振动过程趋于线性极化，此时主特征向量所对应的 空间方向即为质点振动所在的方向；若γ值趋于0，主特征值会和次特征值相 当，极化椭球表现为趋于球形，表示质点在空间振动时的非线性极化越强，此 时质点振动方向没有一个优势方向，通过窗口在地震记录上的滑动，求取整个 剖面上不同采样点处，不同时刻的波的极化特征信息。
[0093] S4、根据步骤S2得到的速度模型，利用射线追踪的方法，计算获得：从 震源发出的波经过每个绕射点到各检波器的波场绕射时间和绕射射线出射方 向；
[0094] S5、利用绕射射线出射方向与主偏振方向之间的差异构造方向权函数；
[0095] S6、在均匀各向同性介质中的克希霍夫积分法偏移的过程中，加入方向权 函数项，得到极化深度偏移成像结果。
[0096] 岩巷超前探测三维极化偏移成像方法的基本原理为：已知均匀各向同性介 质中的标量波动方程如下：
[0097]
[0098] 其中，U(x,y,z,t)描述空间中P(xP,yP,zP)点处的波场函数，v(x,y,z)描述了波传 播空间中的速度分布情况，F＝F(x,y,z)表示求解空间中的已知震源函数。
[0099] 根据惠更斯-菲涅尔原理：某时刻波前面上的所有质点可以看作广义绕射 震源，从这些震源发出的二次子波在下一时刻的一系列位置产生干涉叠加，这 些二次子波干涉叠加的包络即代表了该时刻波前面的新位置。该原理描述了这 样一个事实：当不知道震源确切情况时，若已知某空间的某些点(已知波前面) 的振动情况，可求该空间内待求目标点的振动情况。
[0100] 当在地面进行地震勘探时，对于三维标量波动方程的克希霍夫积分解中的 积分曲面，可以认为其由地面和一个向地下空间无限延伸的部分球面组成，在 无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。定义坐标系为 XOY面与地面重合，Z轴垂直地面指向地下空间为其正方向的情况下，对地面 的外法线方向的导数与对z坐标的导数只相差一个符号，即 此 处θ表示r与地面外法线方向n之间的夹角。当积分曲面内不包含有震源时， 即取F＝0，因此均匀各向同性介质中的克希霍夫积分解的具体形式为：
[0101]
[0102] 其中，第一项由波场的垂直梯度决定，第二项与波场函数值的延迟位有关，叫 做近场源项，它以 衰减，这两项通常在克希霍夫积分法偏移中被忽略，而 只使用第三项远场源项，第三项可以用波场观测数据项和权因子项W来表征：
[0103]
[0104] 其中， 表示权因子，即在积分过程中，各道观测值对最终各成像点幅 值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S表示 积分曲面。
[0105] 根据图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的方向。
[0106] 对于给定的已知速度模型，将模型中的每一个网格点看作时绕射点，而对 于有一定延续长度(面积)的界面，可以将其看作是由许多绕射点组成的一个 整体，从而可以把研究对象统一起来，使用各点的绕射叠加结果来表示反射信 息。对每个网格点进行射线追踪可以找到从每一个网格点(绕射点)到检波点 的绕射射线单位方向向量d(Rijk,Gn)。显然，对于某一个网格点来说，网格点到 各个检波点的方向不尽相同。
[0107] 对于观测系统(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)也与真实的波场传播方向存在 差别。
[0108] 根据上述分析，可以在克希霍夫积分法偏移的过程中加入上述两个方向的 差别有关的权函数K(p,d)，加入了方向权函数的克希霍夫积分偏移公式，具体 形式如下：
[0109]
[0110] 其中，κ表示调节因子，用于调节权因子在成像过程中的压制效果，K(p,d)表 示方向权函数， 表示权因子，为积分过程中各道观测值对最终各成像 点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S 表示积分曲面，由地面和一个向地下空间无限延伸的部分球面组成，同时，在 无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。
[0111] 方向权函数的具体形式如下：
[0112]
[0113] 其中，Θ的具体形式如下：
[0114]
[0115] θ表示射线追踪计算的绕射射线在检波器处的出射方向和根据绕射时间 对地震数据进行极化分析所得的纵波极化方向两方向之间的方向夹角。
[0116] 实施例采用模型A和模型B，模型A为与测线方向成小角度斜交的竖直断 层模型；模型B为与测线方向成大角度斜交的竖直断层模型，两个速度模型的 大小均为x×y×z＝
100m×100m×50m，模型在三个轴方向上的离散网格间距均为 0.5m。
100m×100m×50m，模型在三个轴方向上的离散网格间距均为 0.5m。
[0117] 模型A与模型B参数设置如表1所示。
[0118] 表1、模型A和模型B参数设置
[0119]
[0120] 在图3中为模型A，即与测线方向成小角度斜交的竖直断层模型，断层两 侧为岩性不同的两种介质，断层界面与YOZ平面所成角度角约为26.5°，分界 面法向单位向量为[0121] 在图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子波。
[0122] 具体实施方式为，首先将三分量纵波地震记录进行波场分离，结果如图4 所示，图4中的同相轴代表被分离出来的纵波，图4(a)、图4(b)、图4(c) 依次对应X、Y、Z三分量地震记录。然后将图4所示的分离纵波后的地震记录 进行速度扫描得到三维速度谱，结果如图5所示，图5表示了介质的速度分布 情况。
[0123] 图6为模型A三维速度扫描结果ZOX平面切片图，图6(a)～图6(f) 依次对应于界面法线方向与Y轴夹角范围45°～70°，角度间隔5°；
[0124] 三维速度扫描的具体做法如下：
[0125] 参考图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的对称点，坐标为 SR、RG表示从震源S发出的某一射线在R点反射后到达检波点G点，由反射 定律可知，连接S*、G两点的连线与界面ABC的交点即为反射点R(xR,yR,zR)。
[0126] 为了对图8中所示几何关系进行定量表示，首先，假设观测测线SG与XOZ 平面的交点坐标为(x0,0,z0)，根据上述几何关系，有：
[0127]
[0128] 其次，假设炮点S的零偏移距自激自收旅行时间为t0，且炮点到界面的垂 直距离为H，则有下述关系：
[0129] 2H＝vt0
[0130] 最后，假设炮点关于界面的镜像对称点S*坐标为(x′,y′,z′)。则有向量 显然，向量 与界面的法线方向n平行，因此有如 下关系式：
[0131]
[0132] 根据上述关系，可以用法线方向余弦表示各检波点处的反射旅行时间：
[0133]
[0134] 即：
[0135]
[0136] 从上式可以看出，当测线沿Y轴方向时，测线上任意点的反射旅行时间可 以由炮点自激自收时间t0、炮检点距离(yS-yG)以及界面单位法线方向在Y轴 方向上的投影值m三个参数完全表示，从而能在速度扫描过程中，避免多参 数引起的不便，也便于成图显示和分析。对于时间采样间隔为Δt共N道的共 炮点道集地震记录U(x,t)，设速度扫描计算时窗长度为M个采样点，依据平均 振幅能量准则，在经过速度扫描过程后，可以得到三维速度谱ESpec，具体形式 如下：
[0137]
[0138] 其中，i表示道序号，j表示采样点在窗口内的自然排列顺序序号， 表示各道的反射时间，M表示速度扫描计算时窗 长度，t0
表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向 单位向量的Y分量；
表示炮点自激自收时间、v表示界面上覆介质速度、m表示界面法向 单位向量的Y分量；
[0139] 首先计算各道的反射时间ti，再依据公式计算与t0、v、m对应的速度谱值， 对上述三个变量选定不同值，重复以上步骤即可以计算地震剖面的三维速度 谱。
[0140] 利用协方差矩阵法进行极化分析，如图7所示，求得极化程度参数值和主 特征向量的方向余弦值，图7(a)为计算得到的极化程度，其显示了在同相 轴处极化程度参数γ值接近1，表示在该同相轴所在的时间宽度范围内波场高 度极化，而在只有噪声存在的地方，参数γ值接近0，波场极化程度低。根据 极化方向的方向余弦图7(b)、(c)、图7(d)可知，沿空间X轴方向的方向 余弦整体上大于沿Y方向的方向余弦，而沿Z轴方向的方向余弦接近于零值， 即通过模型1模拟得到的波场中的纵波成分主要集中在X分量上，而Z分量 上的振幅值几乎为零，通过图7(c)可知纵波波场的Y分量存在一渐变过程： 当检波器越远离炮点位置，纵波波场的Y分量由大变小再增加，且对应的方 向余弦在70道附近发生符号的改变。
[0141] 具体做法如下：
[0142]
[0143] 其中，W表示时窗计算长度，u(x,y,z)表示三分量地震矢量波场， 表示方差， 表示协方差，
表示数学期望， 表示采 样点数，N
表示每道的总采样点数，ΔT表示时间采样间隔；
表示数学期望， 表示采 样点数，N
表示每道的总采样点数，ΔT表示时间采样间隔；
[0144] 协方差矩阵为一对称矩阵，求取协方差矩阵的特征值λ1、λ2、λ3，且限定 λ1≥λ2≥λ3，及其对应的单位特征向量为p1、p2、p3，协方差矩阵的三个特征值 分别对应极化椭圆的三个轴的长度，最大特征值λ1对应极化椭圆的主轴，相应 的主特征向量p1表示质点椭圆偏振的主方向，其余同理；
[0145] 该时刻波的偏振类型用偏振系数γ(t)来衡量：
[0146]
[0147] 若γ值趋于1，主特征值λ1会远大于其它次特征值，并表现为极化椭球趋 于细长，表示质点的空间振动过程趋于线性极化，此时主特征向量所对应的空 间方向即为质点振动所在的方向；若γ值趋于0，主特征值会和次特征值相当， 极化椭球表现为趋于球形，表示质点在空间振动时的非线性极化越强，此时质 点振动方向没有一个优势方向，通过窗口在地震记录上的滑动，求取整个剖面 上不同采样点处，不同时刻的波的极化特征信息。
[0148] 然后通过地震极化信息求模型A的反射界面参数，如表2所示。
[0149] 表2模型A极化分析计算得到的反射界面参数
[0150]
[0151] 通过地震极化信息求取反射界面几何参数的具体做法如下：
[0152] 设已知直线束Li，(i＝1，2，3，…，N)，每一直线上都有一已知坐标的点 Pi，其坐标为(xi,yi,zi)，且各直线的单位方向向量为vi＝(li,mi,ni)也已知。设 Q(x,y,z)为空间内的任意点，其到某一直线的距离hi可以表示为：
[0153]
[0154] 其中i、j、k为沿X、Y、Z三个坐标轴正方向的单位向量，且有：
[0155] α＝(y-yi)ni-(z-zi)mi
[0156] β＝(z-zi)li-(x-xi)ni
[0157] γ＝(x-xi)mi-(y-yi)li
[0158] 为了求得直线束Li的公共交点，取点Q到所有直线的距离之和 为最小，即令： 且 且
[0159]
[0160] 若令l、m、n分别表示如下列向量：
[0161]
[0162] 又因为vi为单位向量，有 上式左侧系数矩阵可进一步 化为(其中I为N阶单位矩阵)：
[0163]
[0164] 因此，方程可化为：
[0165]
[0166] 其中，N表示假设已知的直线束Li的数量，(i＝1，2，3，…，N)，Pi(xi,yi,zi)表 示每一直线上的某一已知坐标的点，Q(x,y,z)表示空间内的任意点， vi＝(li,mi,ni)表示各直线的单位方向向量，I表示N阶单位矩阵；
[0167] 对纵波进行极化分析，获取波场在不同位置的传播方向后，通过波场在不 同位置的传播方向进行反向延长，求取交点，通过解如上线性方程，求得反射 射线的反向交点Q的坐标，交点与炮点关于界面镜像对称，由交点与炮点确 定的直线的单位方向的向量v为界面法线方向，所述界面法线方向及交点与炮 点构成线段的中点M(xM,yM,zM)，确定反射界面的解析表达式；
[0168] 通过野外采集的三分量地震数据空间坐标范围大小可以给出三维速度模 型空间(或偏移目标区域)的大小限制，利用空间平面与速度模型空间的各棱 边直线的交点来估计反射界面的具体位置，通过求取空间反射界面与给定模型 A空间棱边的交点，得到表3所示的四个交点坐标。
[0169] 表3模型A空间棱边的交点
[0170]
[0171] 通过表2与表3上述求得的反射界面几何参数与图5中三维速度谱结合， 即可得到经速度分析后的偏移速度模型。
[0172] 已知速度模型后，利用射线追踪方法可以得到旅行时与绕射波出射方向， 利用绕射波出射方向和主偏振方向之间的差异可以构造方向权函数，在克希霍 夫积分偏移中加入方向权函数即可完成偏移中的假象压制，结果如图9～图13 所示。为了验证偏移算法有效性，首先利用真速度模型进行克希霍夫叠前深度 偏移成像。图9为利用真速度模型A进行极化偏移在深度z＝25m处XOY截面 的成像结果，虚线表示真实介质分界面与z＝25m处的水平切面的交线，点划 线表示真实介质分界面关于观测测线对称的面与z＝25m处的水平切面的交线， 垂直的线表示地震观测测线，可以看到，成像区域聚集真实界面(虚线)同一 侧，而在关于观测测线对称处(点划线)的成像假象几乎被完全压制；
[0173] 然后利用偏移速度模型进行克希霍夫叠前深度偏移成像，图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时， 对于某一分量的偏移结果，若仍以相同比例尺进行显示，在前述模型情形下， 已看不到明显的偏移假象存在，真实成像目标区域也进一步被缩小。
[0174] 利用加入了方向权函数的克希霍夫积分进行极化偏移，具体做法如下：
[0175]
[0176] 其中，κ表示调节因子，用于调节权因子在成像过程中的压制效果，K(p,d)表 示方向权函数， 表示权因子，为积分过程中各道观测值对最终各成像 点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(xP,yP,zP)点处的波场函数，S 表示积分曲面，由地面和一个向地下空间无限延伸的部分球面组成，同时，在 无限部分球面上的积分对所求点的波场值的贡献可以忽略不计。
[0177] 地震勘探中接收到的地震波场看作是地下无数二次绕射点源发出的波在 接收点处的叠加结果。若已知地下的速度结构，则可以将偏移目标区域离散成 一定数量的离散点，通过射线追踪的方法获得从震源发出的波经过每个绕射点 到各检波器的波场绕射时间，进而根据各离散点的绕射时间从地震记录上选取 振幅值进行叠加，即可以得到该点的偏移成像结果。对所有离散点进行上述过 程，则可以完成整个目标区域的偏移成像工作。
[0178] 方向权函数的具体形式如下：
[0179]
[0180] 其中，Θ的具体形式如下：
[0181]
[0182] θ表示射线追踪计算的绕射射线在检波器处的出射方向和根据绕射时间 对地震数据进行极化分析所得的纵波极化方向两方向之间的方向夹角。
[0183] 根据图14所示，竖直断层模型B为与测线方向成大角度斜交的断层，断 层两侧为岩性不同的两种介质，断层界面与YOZ平面所成角度角约为63.43°， 用于接收的检波器数量为121道，采样总时长为200ms，其余参数与模型A 的参数一样，同时，模型B的成像方法和过程也与模型A的一致。
[0184] 具体实施方式为，首先将三分量纵波地震记录进行波场分离，结果如图 15所示，图15中的同相轴代表被分离出来的纵波，图15(a)、(b)、(c)依 次对应X、Y、Z三分量地震记录。然后将图15所示的分离纵波后的地震记录 进行速度扫描得到三维速度谱，结果如图16所示。其中，图17为模型B三维 速度扫描结果ZOX平面切片图，图17(a)～图17(f)依次对应于界面法线 方向与Y轴夹角范围45°～70°，角度间隔5°；
[0185] 然后利用协方差矩阵法进行极化分析，结果如图18所示，求得极化程度 参数值和主特征向量的方向余弦值。结合图15所示纵波地震记录及图18所示 极化分析结果进行分析可知，纵波信号极化程度大，且主要在XOY平面内极 化。
[0186] 然后利用地震的极化信息来求取模型B的反射界面参数，如表4所示。
[0187] 表4模型B极化分析计算得到的反射界面参数
[0188]
[0189] 通过求取空间反射界面与给定模型B空间棱边的交点，得到表5所示的四 个交点坐标。
[0190] 表5模型A空间棱边的交点
[0191]
[0192] 通过表4与表5上述求得的反射界面几何参数与图16中三维速度谱结合， 即可得到经速度分析后的偏移速度模型。
[0193] 已知偏移速度模型后，利用射线追踪方法可以得到旅行时与绕射波出射方 向，利用绕射波出射方向和主偏振方向之间的差异可以构造方向权函数，在克 希霍夫积分偏移中加入方向权函数即可完成偏移中的假象压制，具体做法与模 型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分量地震记录的成像结果；
[0194] 图22为图19、图20、图21的偏移成像结果切片图，切片深度z＝25m， 虚线为反射界面所在位置。从图22(a1)、(a2)、(a3)可知，利用真速度模 型进行偏移，偏移结果能准确表示反射界面所在位置，但是存在偏移假象，且 假象与正确的像能量相当，若不利用极化信息来对波场方向进行分辨，则无法 确定真正成像位置的方位；图22(b1)、(b2)、(b3)表明，利用极化偏移能 较好压制偏移假象，但是，在压制假象的同时，引入了噪声，对成像结果造成 了干扰；图22(c1)、(c2)、(c3)所示为利用速度分析结果构建的速度模型 进行偏移结果，图中也显示了极化偏移对偏移假象的压制效果，但是，由于速 度模型不够准确，造成成像位置的准确程度不如图22(a1)、(a2)、(a3)中 所示的结果。
[0195] 以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于 此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到 变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应 所述以权利要求的保护范围为准。
| 标题 | 发布/更新时间 | 阅读量 |
|---|---|---|
| 基于多波束轮询的传输方法及通讯设备 | 2020-05-11 | 189 |
| 辐射源 | 2020-05-11 | 497 |
| 数据处理方法及相关设备 | 2020-05-12 | 214 |
| 时间提前值的管理 | 2020-05-13 | 994 |
| 带式供料器 | 2020-05-11 | 848 |
| 一种通信方法和装置 | 2020-05-12 | 15 |
| 高速移动性通信系统和方法 | 2020-05-13 | 63 |
| 用于数据传输的方法、发送端设备和接收端设备 | 2020-05-11 | 870 |
| 用于蒸气产生装置的可感应加热烟弹 | 2020-05-11 | 752 |
| 用于空对地系统的切换改进 | 2020-05-13 | 653 |
高效检索全球专利专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台，是提供专利分析，专利查询，专利检索等数据服务功能的知识产权数据服务商。
我们的产品包含105个国家的1.26亿组数据，免费查、免费专利分析。
分析报告
专利汇分析报告产品可以对行业情报数据进行梳理分析，涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度，系统通过AI智能系统对图表进行解读，只需1分钟，一键生成行业专利分析报告。
QQ群二维码
意见反馈
意见反馈