阅读说明：本技术 一种基于克希霍夫积分解的地震槽波深度偏移的成像方法 (Imaging method for seismic channel wave depth migration based on kirchhoff product decomposition ) 是由 李志勇 于 2019-11-20 设计创作，主要内容包括：本发明提供了一种基于克希霍夫积分解的地震槽波深度偏移的成像方法,包括输入炮集槽波地震记录；对槽波地震信号进行包络计算；在所获取的槽波地震记录中进行子波反褶积运算,压缩槽波,使其变换为有正有负的波形；建立槽波的速度模型；运用波前快速推进法计算旅行时和射线路径；对槽波数据进行叠前克希霍夫偏移成像。本发明通过希尔伯特包络计算、同态法子波提取和整形反褶积各个算法实现了“压缩”和消除频散的槽波；适用于煤矿井下煤巷的地震槽波超前探测条件下的克希霍夫深度偏移成像,首次为槽波反射波超前探测提供了一种成像手段,使得地质和地球物理工程师能够直观地推断工作面前方复杂的地质情况,为煤巷的安全掘进提供了重要的技术保障。(The invention provides an imaging method of seismic channel wave depth migration based on kirchhoff product decomposition, which comprises the steps of inputting shot gather channel wave seismic records; carrying out envelope calculation on the channel wave seismic signals; performing wavelet deconvolution operation in the acquired groove wave seismic record, and compressing the groove wave to convert the groove wave into a positive waveform and a negative waveform; establishing a velocity model of the channel wave; calculating travel time and ray paths by using a wavefront fast propulsion method; and performing pre-stack kirchhoff migration imaging on the slot wave data. The invention realizes the 'compression' and the elimination of frequency dispersion groove waves through various algorithms of Hilbert envelope calculation, homomorphic wavelet extraction and shaping deconvolution; the method is suitable for kirchhoff depth migration imaging under the seismic channel wave advanced detection condition of the coal mine underground coal roadway, provides an imaging means for channel wave reflected wave advanced detection for the first time, enables geological and geophysical engineers to intuitively infer complex geological conditions in front of a working face, and provides important technical support for safe tunneling of the coal roadway.)

一种基于克希霍夫积分解的地震槽波深度偏移的成像方法

技术领域

本发明属于地质物探技术领域，具体涉及用于在巷道存在的情况下，实 现槽波超前探测的一种基于克希霍夫积分解的地震槽波深度偏移的成像方 法。

背景技术

自1972年8月，在美国芝加哥召开快速掘进与隧道工程会议后，隧道 施工中的地质条件超前预报日益受到重视，上世纪80年代以来，德国、英 国、澳大利亚、匈牙利、捷克斯洛伐克、苏联、美国等都先后设立专门机构， 从事槽波地震勘探(in-seam seismic，ISS)研究，促进了井下试验、理论研 究和仪器研制的蓬勃开展，逐步地把ISS投入了实际应用(Giese and Klose, 2003)。目前，隧道地震超前预报、地震超前探测技术在地质工程的很多方 面有着广泛的运用，但目前地震超前探测技术研究仍是以隧道工程为主，对 于煤矿巷道应用鲜有使用，其原因为：多数地震超前探测技术主要以体波为 主，在煤层中槽波能量非常强，体波能量相对较弱，无法利用体波去对煤巷 完成超前探测。如果利用槽波进行煤巷超前探测会比体波更加有效。然而如 何将已有的成熟隧道探测技术移植到煤矿井下，还需要大量的基础研究与试 验。

偏移方法可分为时间偏移、深度偏移和叠前偏移、叠后偏移等。叠后偏 移方法应用于水平叠加之后，需要两个理论假设：一、输入数据要满足水平 层状介质假设条件；二、其输入为零炮检距的自激自收记录，那么当地层结 构形态复杂时，这些假设条件就难以被满足，而叠前时间偏移技术避免了这 些理论假设，更适应实际复杂情况，叠前偏移方法有利于提高构造解释与成 图精度，确保叠前属性提取和叠后地震数据反演结果的真实可靠性。

随着地震勘探理论研究的不断深入，地震偏移成像技术进入空前发展的 新阶段，逆时偏移技术、叠前深度偏移技术已经成为地震偏移的研究主流， 随之也将会衍生出各种其他的地震偏移技术。

发明内容

根据上述阐述，本发明的目的在于提供一种基于克希霍夫积分解的地震 槽波深度偏移的成像方法，是针对纵横波分离方法，对观测系统、施工条件 适应性较强，且精度较高。

本发明提供的技术方案：

一种基于克希霍夫积分解的地震槽波深度偏移的成像方法，包括如下步 骤：

S1、输入炮集槽波地震记录，对三分量槽波地震记录进行预处理，完成 振幅恢复、振幅平衡与滤波处理；

S2、对槽波地震信号进行包络计算，将波形序列改变成能量序列，使包 络后地震波频率远远低于原始记录中主要能量的频率；

S3、在所获取的槽波地震记录中提取子波，进行子波反褶积运算，压缩 槽波，并使其变换为有正有负的波形；

S4、建立槽波的速度模型；

S5、运用波前快速推进法计算旅行时和射线路径；

S6、对槽波数据进行叠前克希霍夫偏移成像。

上述技术方案中，所述步骤S2中，采用Z分量记录进行数据处理和断 层定位。

上述技术方案中，所述步骤S1中，滤波处理的目的是滤除直达槽波， 压制高能量的直达波干扰，提高后续槽波深度偏移成像的效果。

上述技术方案中，所述步骤S2中，进行包络计算的目的是消除槽波频 散特征所引起的波列不稳定性及叠加后的相消干涉结果，包络计算方法如 下：

若x(t)为实际地震道信号，构造一个复信号c(t)

c(t)＝x(t)+iy(t)

则信号的包络：

式中y(t)为实信号x(t)的希尔伯特变换：

上述技术方案中，所述步骤S3中，采用同态法提取子波，提取方法为：

设一道地震记录为x(t)，对其进行傅里叶变换，并取对数，然后进行傅 里叶反变换，得：

x'(t)＝IFT[lnX(ω)]

x_i'(t)具有非常重要的性质：x_i'(t)的任一段记录，由于子波序列的复赛谱 近似不变，而反射系数序列的复赛谱是随机变化的，所以在时间域内对其进 行线性滤波，可得到近似子波的时间域时间序列：

w(t)＝L[x'(t)]

其中，L[·]表示线性滤波器，其滤波目的是保留子波，由此可以获得子 波

lnW′(ω)＝FT[w′(t)]

最后，通过傅里叶逆变换，得到估计的地震子波：

w(t)＝IFT(ln(W′(ω))。

上述技术方案中，所述步骤S3中，子波整形反褶积的过程如下：

设x(k)为要整形的槽波地震道，设计一个整形算子m(k)作用于x(k)，

z(k)＝m(k)*x(k)

可得整形反褶积的结果z(k)；

上式中“*”表示反褶积。

上述技术方案中，所述步骤S6中，叠前克希霍夫偏移的基本公式是：

其中，

表示权因子，即在积分过程中，各道观测值对最终各成 像点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(x_P,y_P,z_P)点处的波场函数， S表示积分曲面。因子

可以通过射线追踪计算来获得；

实现叠前克希霍夫偏移成像的算法主要步骤是：

Q1，将地下介质网格化，地下每个网格看作一个绕射点，即成像点；

Q2，基于步骤S4建立的槽波速度模型，通过波前快速推进法计算地下 网格剖分绕射点的时间场，得到绕射点到炮点以及绕射点到检波点的时间；

Q3，根据Q2中计算得到的绕射点的时间场，计算得到炮点到绕射点的 最短时间路径以及绕射点到检波点的最短时间路径，根据路径信息可以得到 克希霍夫偏移算法中的振幅权值；

Q4，循环计算地下所有绕射点的时间场，得到所有绕射点到炮点和检波 点的时间和传播路径，保存时间和路径信息；

Q5，根据绕射点的时间场得到绕射点的绕射时距曲线以及射线路径，用 克希霍夫积分偏移公式成像，叠加得到成像点的能量；

Q6，循环计算完成地下介质所有绕射点的成像，得到叠前克希霍夫深度 偏移成像结果。

本发明的优点有：通过希尔伯特包络计算、同态法子波提取和整形 反褶积各个算法实现了“压缩”和消除频散的槽波；适用于煤矿井下煤 巷的地震槽波超前探测条件下的克希霍夫深度偏移成像，首次为槽波反 射波超前探测提供了一种成像手段，使得地质和地球物理工程师能够直 观地推断工作面前方复杂的地质情况，为煤巷的安全掘进提供了重要的 技术保障。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实 施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面 描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲， 在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本发明优选实施例的流程图。

图2为实施例之模型1的三维显示图。

图3为模型1的XOY水平切平面(Z＝5m)。

图4为模型1模拟计算得到的三分量槽波地震记录。

图5为模型1提取的反射槽波及包络计算结果和反褶积计算结果。

图6为模型1二维克希霍夫深度偏移结果对比。

图7为实施例模型2三维显示图。30°断层速度模型三维显示图。

图8为模型2的XOY水平切平面图示(Z＝5m)。

图9为根据模型2数值计算得到的三分量槽波地震记录。

图10为模型2提取的反射槽波及包络计算结果和反褶积计算结果。

图11为模型2二维克希霍夫深度偏移结果。

具体实施方式

下面将结合本发明的附图，对本发明的技术方案进行清楚、完整地描述， 显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。 基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下 所获得的所有其他实施例，都属于本发明保护的范围。

本发明基本原理是：常规的地震偏移方法无法适用于反射槽波的偏移， 因为地震槽波是煤层中发育的一种由体波干涉产生的导波，具有能量强、传 播距离远、频散的特点，频散的存在使常规的波动方程偏移的结果变得不稳 定。

本发明的基本思想是将频散的槽波经过处理，变成没有频散的波形，然 后再进行槽波的叠前深度偏移。

消除频散的基本原理是，利用希尔伯特变化，对槽波地震信号进行包络 计算、再用同态法提取子波，并整形反褶积以实现频散槽波的压缩，得到无 频散且波形确定的槽波。

在此基础上，本发明提出的基于波动方程的克希霍夫积分解实现槽波的 叠前深度偏移的算法如下：

一种基于克希霍夫积分解的地震槽波深度偏移的成像方法，包括如下步 骤：

S1、输入炮集槽波地震记录，对三分量槽波地震记录进行预处理，完成 振幅恢复、振幅平衡与滤波处理；

S2、对槽波地震信号进行包络计算，将波形序列改变成能量序列，使包 络后地震波频率远远低于原始记录中主要能量的频率；

S3、在所获取的槽波地震记录中提取子波，进行子波反褶积运算，压缩 槽波，并使其变换为有正有负的波形；

S4、建立槽波的速度模型；

S5、运用波前快速推进法计算旅行时和射线路径；

S6、对槽波数据进行叠前克希霍夫偏移成像。

具体为：

(1)Z分量槽波数据

根据我们对地震槽波的研究，在超前探测中采集的三分量槽波地震记 录，其Z分量主要由Rayleigh型槽波构成，波形较为单一，用于数据处理和 成像相对简单，故一般采用Z分量记录进行偏移成像和断层定位；

(2)包络计算

对地震槽波信号进行包络计算的目的是，消除槽波频散特征所引起的波 列不稳定性。

包络计算的方法用希尔伯特变换将槽波信号相移90°，然后槽波信号作 为实部，相移后的信号作为虚部，构造一个复信号，复信号的振幅即是原地 震信号的包络。

具体过程如下。

设有一道地震槽波信号x(t)，做希尔伯特得到

构造复地震信号c(t)：

c(t)＝x(t)+iy(t) (2)

取函数c(t)的模得到包络e(t)，即：

(3)子波提取

采用同态法提取子波，同态法属于基于高阶的统计性子波提取方法，其 优点是对地震子波具有较少的限制条件，当时窗取得较小时，能克服求对数 谱的相位时遇到的多值问题。

一道地震记录x(k)可以表示为地震子波b(k)和反射系数序列r(k)的褶 积：

x(k)＝b(k)*r(k)

上式的傅里叶变换结果为：

X(ω)＝B(ω)R(ω)

取其对数可得：

ln X(ω)＝ln B(ω)+ln R(ω)

对上式，进行傅里叶反变换得：

IFT[lnX(ω)]＝x′(k)＝b′(k)+r′(k)

上式中所表示的地震信号，在某一截取的时窗中，由于子波序列的复赛 谱近似不变，而反射系数序列的复赛谱是随机变化的，所以在时间域内对其 进行线性滤波L[·]，可得到近似子波的时间域信号：

y_i'(k)＝L[x'(k)]

然后再次傅里叶变换，再次得到对数谱：

ln Y′(ω)＝FT[y_i'(k)]

最后，进行指数变换和傅里叶逆变换，得到估计的地震子波：

Y(ω)＝exp(lnY′(ω))

y(k)＝IFT(Y(ω))

(4)整形反褶积

设x(k)为经过包络计算后的地震数据，设计一个整形算子m(k)，通过 反褶积

z(k)＝m(k)*x(k)

获得整形后的地震数据道z(k)，其中m(k)可以通过下式

R_xxM＝R_yx

的计算获得，其中R_xx为x(k)自相关函数矩阵，R_yx为x(k)与先前提取的 地震子波y(k)的互相关函数矩阵，M为整形算子向量，与m(k)的关系是：

解矩阵方程：

可获得m(k)。

(5)克希霍夫深度偏移

均匀各向同性介质中的标量波动方程的Kirchhoff积分解的具体形式为：

式中，第一项由波场的垂直梯度决定，第二项与波场函数值的延迟位有关， 叫做近场源项，它以

衰减。在实际工作中，反射波的传播距离一般都足够 远，第一和第二项通常在克希霍夫积分法偏移中被忽略，而只使用第三项远 场源项。第三项可以用波场观测数据项和权因子项W来表征：

其中，

表示权因子，即在积分过程中，各道观测值对最终各成 像点幅值的贡献比例，U(x,y,z,t)表示空间中任一P(x_P,y_P,z_P)点处的波场函数， S表示积分曲面。因子可以通过射线追踪计算来获得。

实际数据的偏移成像过程如下：

实现叠前克希霍夫偏移成像的算法主要步骤是：

Q1，将地下介质网格化，地下每个网格看作一个绕射点，即成像点；

Q2，基于步骤S4建立的槽波速度模型，通过波前快速推进法计算地下 网格剖分绕射点的时间场，得到绕射点到炮点以及绕射点到检波点的时间；

Q3，根据Q2中计算得到的绕射点的时间场，计算得到炮点到绕射点的 最短时间路径以及绕射点到检波点的最短时间路径，根据路径信息可以得到 克希霍夫偏移算法中的振幅权值；

Q4，循环计算地下所有绕射点的时间场，得到所有绕射点到炮点和检波 点的时间和传播路径，保存时间和路径信息；

Q5，根据绕射点的时间场得到绕射点的绕射时距曲线以及射线路径，用 克希霍夫积分偏移公式成像，叠加得到成像点的能量；

Q6，循环计算完成地下介质所有绕射点的成像，得到叠前克希霍夫深度 偏移成像结果。

根据图2和图3所示的型1和模型2具体形状和测线位置，本发明以模 型1和模型2作为实施例，表1为模型1和模型2所采用的模型物性参数表。

表1为模型1和模型2所采用的模型物性参数

在上述两个实例的正演模拟中，采用交错网格有限差分正演模拟获得槽 波记录。

实施例1：90°竖直断层模型1大小为40m×80m×10m，X轴方向40m、 Y轴方向80m、Z轴方向10m；网格间距dx＝dy＝dz＝0.5m；时间采样间隔 0.05ms，记录时长为100ms；震源采用主频为500Hz的雷克子波；煤层大小 为40m×40m×2m，中心位置在(20,20,5)处；巷道大小为3m×30m×3m， 中心位置在(17.5,15,5)处；断层在巷道前方40m处，断层走向与巷道掘进 方向垂直；震源(20,5,5)位于距巷道壁1m深的煤层中，测线位于x＝20m， z＝5m，y＝0～29m，共30道检波器接收，道间距为0.5m。

根据图4所示，通过有限差分正演数值模拟得到的三分量槽波地震记录， Z分量主要是Rayleigh型槽波，波形较为单一，接收Z分量信号用于数据处 理较为简便。偏移实验采用Z分量记录进行数据处理、成像和断层定位。

从震源激发的地震槽波，直接传播到检波器的称为槽波直达波，继续传 播并经断层反射回来的为槽波反射波。

根据图5所示，三分量合成槽波记录可以发现：反射槽波与直达槽波的 视速度符号相反，即直达槽波具有正视速度，反射槽波具有负视速度。

由于模型1比较简单，可以很明显的分辨出反射槽波，直接截取出Z分 量地震记录中的反射槽波。

取三分量槽波地震记录的前25道作为分析对象。

根据图5所示，图5(a)显示了分离所得的反射槽波；图5(b)所示 为根据反射槽波计算得到的波形包络，图5(c)所示为反褶积计算结果。从 图5(c)可以看出，频散严重的槽波反射波经过包络、反褶积处理后，其频 散效应被大大压缩，反射波同相轴清晰可见，从图5(c)中可以看到两条主 要同相轴，它们之间的时间延迟差别较大，但是它们其实代表来自同一反射 界面的槽波反射信息。

为了验证偏移结果的准确性，采用正演所用煤层速度(v＝2300m/s)制 作均匀速度模型(二维模型大小为X×Y＝40×80m2)以及使用正演所用真 速度模型进行偏移计算，见图6(a)，偏移结果对比如图6所示，图(b)为 槽波偏移处理结果，可见槽波的偏移能够正确恢复断层的位置(虚线为断层 的位置)。

根据图7所示，实施例2：30°断层模型2大小为40m×80m×10m，X 轴方向40m、Y轴方向80m、Z轴方向10m；网格间距dx＝dy＝dz＝0.5m；时间 采样间隔0.01ms，记录时长为120ms；震源雷克子波频率为500Hz；煤层大 小为40m×40m×2m，中心位置在(20,20,5)处；巷道大小为3m×30m×3m， 中心位置在(10,15,5)处；断层在巷道前方60m处，断层走向与巷道掌子 面方向成30度夹角；震源(12.5,5,5)位于距巷道壁1m深的煤层中，测 线位于x＝12.5m，z＝5m，y＝0～29m，共150道接收，道间距0.2m。

图8所示为在深度z＝5m处的水平切片图(XOY面)，且断层位置、巷道 位置、测线、炮点均已在图中标示出。

图9展示了正演数值模拟计算得到的三分量合成槽波记录，从图9中可 以看出，Z分量波形较为单一，主要是Rayleigh型槽波，偏移实验采用Z 分量记录进行数据处理、成像和断层定位。

图10(a)展示了滤波分离出的反射槽波地震记录，及其包络计算结果 (图11(b))和反褶积结果(图11(c))。

图11展示了偏移计算的结果，阴影区表示围岩，斜线区域表示巷道所 处位置，其它位置表示煤层。从图中可以看出，偏移结果清晰地显示出了反 射界面的位置。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限 于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易 想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护 范围应所述以权利要求的保护范围为准。