CN102590862A - 补偿吸收衰减的叠前时间偏移方法 - Google Patents

Abstract

补偿吸收衰减的叠前时间偏移方法，应用于地震勘探中反射地震资料处理。该方法在偏移过程中补偿地球介质粘性、薄层散射导致的高频地震波幅值衰减，恢复被衰减的高频成份，使得中、深层构造的成像分辨率达到与浅层接近的程度，且可保证稳定性并避免噪音放大。该方法不对Q值的空间分布做层状或均匀的假设，通过引入描述地震波幅值吸收衰减的等效Q值，将补偿吸收衰减与叠前时间偏移有效地结合到一起。该方法能基于地面采集的反射地震资料确定等效Q值，规避了现行各类吸收衰减补偿方法在确定Q值上的困难。该方法能更好地指示中、深层构造的断裂和地层沉积样式，对油气、矿产资源勘探有重要应用价值。

Description

补偿吸收衰减的叠前时间偏移方法

技术领域

本发明属于地震勘探中反射地震资料处理技术领域，涉及地震资料处理过程中的叠前偏移成像技术范畴，是一种能在偏移计算中补偿地震波传播的吸收衰减的叠前时间偏移方法。

背景技术

地震勘探中反射地震资料处理过程中，叠前偏移成像是关键的环节，而叠前时间偏移是叠前偏移成像中的一种重要方法。叠前时间偏移方法可对一类断层较为复杂但速度横向变化不是很剧烈的地质构造较好成像。与叠前深度偏移方法相比，除具有较高的计算效率外，其主要的优点是只需使用叠加(均方根)速度；这样可简单地通过速度扫描等方式得到恰当的速度模型，回避了使用叠前深度偏移方法面临的一个主要困难：速度建模。因此，叠前时间偏移方法已成为地震勘探领域广泛应用的关键技术。

实际地球介质存在粘性吸收，地球介质的小尺度非均匀也产生类似于粘性吸收的幅值衰减效应。这些客观存在导致地震波在传播过程中发生幅值的吸收衰减；衰减对地震波的不同频率成份是不同的，频率越高，衰减的越强。因此，地表记录到的、从不同深度反射的地震信号其频带是不同的；这导致构造越深，常规偏移成像的分辨率就越低。在油气勘探中，对小断层、小断裂体系的识别，是认识油气疏导体系，进而识别有利储层的重要环节，因此石油工业界对提高成像分辨率的努力一直在进行着。

为提高地震成像的分辨率，已发展了许多方法，包括成像叠加剖面的谱白化反褶积、非稳态反褶积、基于统计假设或测井资料的各类拓宽频带技术、反Q滤波、粘弹性叠前深度偏移等。谱白化以及各类拓频技术是通过引入地震记录以外的信息提升地震方法的分辨率，尽管可获得更高的视分辨率，但其可靠性尚待提高；非稳态反褶积一般仅在叠加剖面上能得到较好效果；此外，各类拓频技术使用的前提是地震记录的子波是不变的，因此即使应用该类技术，也必须首先补偿地震波的吸收衰减，实现地震子波的一致性。

反Q滤波可同时应用于叠前地震资料和叠后的偏移叠加剖面。这一方法是从补偿地震波幅值的粘性吸收出发，具有坚实的物理基础。但就用于叠前地震资料的反Q滤波而言，它忽略了地震波传播路径对幅值衰减的影响，实际上仅在均匀Q值情况下是准确的；叠后资料的反Q滤波可处理层状Q值模型情况，但由于叠加过程已将不同传播路径，即已将存在不同程度吸收衰减的幅值相叠加，这一处理不能完全消除吸收衰减的影响。

粘弹性叠前深度偏移准确考虑了地震波的传播和粘性衰减过程，理论上是可以较好地补偿地震波的粘性吸收。但粘弹性深度偏移方法需要深度域层间Q值模型，这导致了和叠前深度偏移速度建模相同，甚至更大的困难，因此粘弹性深度偏移的实际应用还存在相当的困难。

就实际勘探中遇到的大量地质目标而言，尽管地表复杂、断裂发育，但地层速度的横向变化并不是很剧烈，叠前时间偏移方法是可以在这些地质目标的勘探中发挥作用的。但现行的叠前时间偏移方法不具有补偿地震波粘性吸收的能力。

在偏移成像过程中恢复地震波被衰减的高频成份是提高地震成像分辨率的关键。它可以真正地提高地震勘探方法对小断层、小断裂体系的识别能力。

发明内容

本发明的目的是：提供一种补偿吸收衰减的叠前时间偏移方法，它可以在偏移过程中补偿地震波传播过程中的高频幅值衰减，恢复被衰减的高频成份，使得中、深层构造的成像分辨率达到与浅层接近的程度；这一方法可处理实际存在的Q值非均匀情况，能基于地面采集的反射地震资料自主确定等效Q值，规避了现行各类衰减补偿方法在确定Q值上的困难。

本发明采用的技术方案是：补偿吸收衰减的叠前时间偏移方法，具体步骤包括：

(1)用单条拖缆或测线记录人工震源激发的反射地震信号，记录到磁带上；

(2)从磁带上读取地震信号，对叠前地震资料做常规的噪音衰减处理。针对部分共中心点，抽取共中心点道集，对抽取的道集作常规的动校正速度拾取，对所得到速度的倒数做横向平均，将得到的横向均匀的速度场作为初始偏移速度；

(3)将噪音衰减处理后的叠前地震资料按偏移距排序、分组，将不同组地震资料放置到集群计算机的不同计算节点上。在每个计算节点，依据初始偏移速度，用常规的叠前时间方法进行偏移计算，收集各计算节点的偏移结果，形成共反射点道集；

(4)对共反射点道集，利用初始偏移速度做反动校，再做动校正得到更新后的速度，对各个CDP点更新后的速度的倒数做空间平滑，所得到的平滑速度场即是最终的偏移速度场；

(5)依据偏移速度场，对噪音衰减处理后的叠前地震资料做常规叠前时间偏移计算。收集各计算节点的偏移结果，形成共反射点道集，对共反射点道集做剩余动校和拉伸切除，得到偏移叠加剖面。在经增益处理后的偏移叠加剖面上选定清晰、连续的同相轴。根据给定的采样间距，确定等间距分布的拾取Q值的CDP点，依据选定的同相轴，在这些CDP点上确定从浅到深的幅值对比点和浅层的参考点，并在各点处选定应用短时傅立叶变换的时窗。从浅到深的幅值对比点处的同相轴需展现由细变粗的变化趋势；

(6)对噪音衰减处理后的叠前地震资料应用时间域快速傅立叶变换，求半导数、施加频率窗衰减带、计算地震道能量、剔除奇异道；

(7)确定等效Q值的取值范围，按1/Q等间距选取系列Q值，用每个Q值，对步骤6处理后的频率域叠前地震资料，做稳定的压噪反Q滤波计算，得到新的时间域地震资料；

(8)依据偏移速度场，对每个Q值对应的反Q滤波后的地震资料做常规叠前时间偏移计算。收集各计算节点的偏移结果，形成共反射点道集，对共反射点道集做剩余动校，切除剩余动校后的拉伸部分。将需拾取Q值的CDP点的相邻共反射点道集叠加，形成各CDP点的超道集。定义小、中、大三个射线参数p₀、p₁、p₂，对全部需拾取Q值的CDP点的幅值对比点和参考点计算偏移距

其中i＝0，1，2，T为该点的单程旅行时，单位秒，v_rms为该点的偏移速度，单位米/秒，射线参数p_i的单位是秒/米，偏移距d_i的单位是米。在各点选定的时窗内，按偏移距范围[0，d₀]、[d₀，d₁]和[d₁，d₂]对超道集分别进行叠加，在每点处得到三个短时时间序列；

(9)重复步骤7和步骤8，收集各幅值对比点和参考点的对应小、中、大射线参数的三个短时时间序列，按1/Q值大小排列，在各点形成三个随1/Q值变化的小、中、大射线参数的短时Q道集；

(10)根据小、中、大射线参数的短时Q道集拾取等效Q值，得到成像区域的等效Q值场；

(11)在每个计算节点，对步骤6处理后的地震资料的全部地震道循环，应用稳定的变频带频率域算法计算各成像点的频率域波场实部；

(12)由各成像点频率域波场的实部，用变频带、保能量成像条件求得各成像点的偏移幅值；

(13)将偏移幅值按偏移距大小累加到存放偏移结果的数组中对应的偏移距上。收集各计算节点的偏移结果，形成全部CDP点的共反射点道集。对该道集做剩余动校正和拉伸切除，将不同偏移距的偏移结果叠加，形成偏移叠加剖面；

(14)通过显示软件将偏移叠加剖面数值转换为地下反射构造的剖面图像，剖面图像将指示中、深层构造的微小断裂体系和精细的地层沉积样式，用于确定地下生、储油构造和识别油气储层。

所述的对噪音衰减处理后的叠前地震资料应用时间域快速傅立叶变换，求半导数、施加频率窗衰减带、计算地震道能量、剔除奇异道是这样实现的：在每个计算节点，对全部地震道循环，应用快速傅立叶变换算法，对地震道的时间序列离散信号做傅立叶变换；做傅立叶变换时，根据时间采样点数，确定满足快速傅立叶算法的采样点数M，增加采样点时添加零值。设Δτ是地震资料的时间采样，单位秒，则角频率采样Δω＝2π/(MΔτ)，其单位是弧度/秒；令地震资料有效频带的下和上限分别是ω_min和ω_max，其单位是弧度/秒，通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近。

设

是频率域的地震道，采用拖缆记录时其单位是帕，采用测线记录时其单位是米/秒；对i＝l₁，l₂循环，求半导数后的地震道为

式中j是单位虚数，Δω是角频率采样，单位弧度/秒。

设频率窗低频端的衰减带含m_d个频率采样点，对i＝l₁，l₁+m_d，计算

在频率窗的高频端，对i=l₂+1，l₂+15，计算

对每一地震道，将存储(l₂+15)-l₁+1个复数值。

计算地震道的能量

将E₀记入该地震道的道头字中。

完成对所有地震道的循环，即得到计算半导数和施加频率窗衰减带后的频率域叠前地震资料。

地震资料是按偏移距排序分配到各个计算节点的。对每一组共偏移距地震道，设定单个地震道能量的变化间距，按E₀的大小将地震道分为若干组，一般取50～100组。计算每组中平均的地震道道数，若最大和最小能量组中的道数小于平均数的百分之一，则判定该组地震道均是奇异道，予以剔除；从最大和最小能量的两端剔除奇异道组不限于一次，连续进行，至到最大和最小能量组中的道数超过平均数的百分之一停止。

所述的确定等效Q值的取值范围，按1/Q等间距选取系列Q值，用每个Q值，对步骤6处理后的频率域叠前地震资料，做稳定的压噪反Q滤波计算，得到新的时间域地震资料是这样实现的：设等效Q值的最小、最大可能值为Q_min、Q_max，等效Q值是无量纲参数，则选取的N个系列Q值为

$\frac{1}{Q_i}=\frac{1}{Q_{\max }}+\frac{i-1}{N-1}(\frac{1}{Q_{\min }}-\frac{1}{Q_{\max }}),i=1,N$

令地震记录满足快速傅立叶算法的时间采样点数为M，Δτ是地震资料的时间采样，单位秒，则角频率采样Δω＝2π/(MΔτ)，其单位是弧度/秒；令地震资料有效频带的下和上限分别是ω_min和ω_max，其单位是弧度/秒，通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近；根据地震资料的信噪比，确定幅值补偿的最大倍数A，一般取为2000.0，A是无量纲数。

对每个Q_i，计算

n₁＝lnA/Δ，n₂＝(lnA+0.2)/Δ，n₃＝l₂M；若n₁≥n₃，则n₁＝n₃和n₂＝n₃，若n₂≥n₃，则n₂＝n₃。建立数组a_i(i＝l₁，n₃)，有

a_i＝exp(iΔ)，i≤n₁

a_i＝A(1-lnA-2.5(lnA)²)+A(1+5lnA)iΔ-2.5A(iΔ)²，n₁＜i≤n₂

a_i＝1.1A，n₂＜i≤n₃

对每一地震道，在对该道做反Q滤波前，事先计算地震道的模

对时间采样点循环，τ＝mΔτ(m＝1，M)，τ的单位是秒，计算

d₀＝cos(Δω·τ)，b₀＝sin(Δω·τ)

d₁＝cos(l₁Δω·τ)，b₁＝sin(l₁Δω·τ)

对频率循环，ω＝iΔω(i＝l₁，l₂)，递推计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

计算n＝m·i，在数组a_i中拾取第n个元素值a_n，得波场实部为

式中

和

分别是步骤6处理后地震资料的实部和虚部。

在上述频率循环中记录峰值频率点对应的补偿因子a_p和最大频率点对应的

据此计算整数 $k=\mathrm{int}\left(\sqrt{\frac{1}{0.06}\ln (a_{(m\·l_2)}/a_p)+25}\right);$ 若k＞15，则令k＝15。

继续对频率循环，ω＝iΔω(i＝l₂+1，l₂+k)，同上由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

进而计算

进行累加计算

从该地震道的道头字中拾取E₀，计算

P(τ)＝(E₀/E_m)P₀(τ)

完成全部时间采样点循环，得到的P(τ)就是新的时间域地震资料的一个地震道，其单位与原始地震道相同。在各个计算节点对所有地震道循环，就得到对应Q_i的新的时间域地震资料。

所述的根据小、中、大射线参数的短时Q道集拾取等效Q值，得到成像区域的等效Q值场是这样实现的：(1)对短时Q道集中的各道做傅立叶变换，得到频率信号；(2)计算傅立叶变换的模，在有效频带内，求模的自然对数。定义有效频带为

其中m_d为低频端的衰减带所含的频率采样点数，

由各短时Q道集对应的射线参数p_i，单位秒/米，对应幅值对比点的单程旅行时T，单位秒，偏移速度v_rms，单位米/秒，各道对应的Q_i，以及频带上限l₂Δω得到：

${\stackrel{\‾}{\mathrm{\ω}}}_{\max }=\min \{\frac{Q_i\ln A}{T}\sqrt{1-v_{\mathrm{rms}}^2{p}_i^2},l_2\mathrm{\Δ\ω}\}.$

为计算的稳定性，模的自然对数由如下方法计算：令D₀是有效频带内模D(ω)的平均值，当D(ω)≥0.22D₀，直接计算其对数ln[D(ω)]；当D(ω)＜0.22D₀时，令d＝5D(ω)/D₀，计算

$\ln \left[D\left(\mathrm{\ω}\right)\right]=\ln \left({0.2D}_0\right)+(d-1)(1-\frac{1}{2}(d-1)+\frac{1}{3}{(d-1)}^2)$

(3)对计算得到的对数做中值滤波处理；(4)收集同属于同一CDP点的短时Q道集的随频率变化的对数值，将各幅值对比点的各道的对数值与浅层的参考点处对应道的对数值相减。在有效频带内，对随频率变化的相减结果先做高斯平滑，然后用中心差分法计算其对频率的一阶导数。对得到的导数值做中值滤波，然后求导数的算术平均值，平均值的绝对值就是该道的频率导数值。在求算术平均值时，若有效频带高频端出现没有相反符号峰值对应的峰值，应将该峰值段剔除；(5)形成需拾取Q值的全部CDP点的频率导数道集图。该图做法如下：在各个幅值对比点的时间深度上，在指定的窗口内，将小射线参数Q道集各道的频率导数值做成随Q值变化的曲线，把中、大射线参数Q道集各道的频率导数值标在图上对应的位置上；(6)基于这一图形，在频率导数值接近零的邻域，考虑小、中、大射线参数结果间的接近程度，直接拾取等效Q值。拾取时要同时考虑等效Q值的横向连续性并保证等效Q值应随旅行时增加而增加；(7)对拾取的Q值，做有关1/Q的空间平滑，即得到成像区域的等效Q值场。

所述的在每个计算节点，对步骤6处理后的地震资料的全部地震道循环，应用稳定的变频带频率域算法计算各成像点频率域波场实部是这样实现的：令地震记录满足快速傅立叶算法的时间采样点数为M，Δτ是地震资料的时间采样，单位秒，则角频率采样Δω＝2π/(MΔτ)，其单位是弧度/秒；令地震资料有效频带的下和上限分别是ω_min和ω_max，其单位是弧度/秒，通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近；根据地震资料的信噪比，确定幅值补偿的最大倍数A，一般取为2000.0，A是无量纲数；设拾取等效Q值时采用的1/Q间距为 $\mathrm{\Δq}=\frac{1}{N-1}(\frac{1}{Q_{\min }}-\frac{1}{Q_{\max }}),$ Δq是无量纲数，计算

其中(Q_eff)_min是由Q值场确定的成像区域内最小的等效Q值；若n₁≥n₃，则n₁＝n₃和n₂＝n₃，若n₂≥n₃，则n₂＝n₃；建立数组a_i(i＝1，n₂)，有

a_i＝exp(iΔ)，i≤n₁

a_i＝A(1-lnA-2.5(lnA)²)+A(1+5lnA)iΔ-2.5A(iΔ)²，n₁＜i≤n₂

对每一地震道，在对该道成像前，事先计算地震道的模

采集是在二维剖面上进行的。令炮点水平坐标为x_s，单位米，检波点水平坐标为x_r，单位米，成像点的水平坐标和单程旅行时分别为x₀和T，其单位分别是米和秒，成像点处的偏移速度为v_rms，单位米/秒，成像点处的等效Q值为Q_eff；则炮点至成像点再反射到检波点的走时为

$\mathrm{\τ}=T(\sqrt{1+\frac{{(x_s-x_0)}^2}{{\left(v_{\mathrm{rms}}T\right)}^2}}+\sqrt{1+\frac{{(x_r-x_0)}^2}{{\left(v_{\mathrm{rms}}T\right)}^2}})$

成像仅需要频率域波场的实部，成像点(x₀，T)处的频率域波场实部由如下快速方法求得：

计算

d₀＝cos(Δω·τ)，b₀＝sin(Δω·τ)

d₁＝cos(l₁Δω·τ)，b₁＝sin(l₁Δω·τ)

$b_3=\frac{\mathrm{\Δ\ω}\·\mathrm{\τ}}{{2Q}_{\mathrm{eff}}\mathrm{\Δ}}$

$\stackrel{\‾}{l_2}=\mathrm{int}(n_2/b_3)$

若 $\stackrel{\‾}{l_2}>l_2,$ 则 $\stackrel{\‾}{l_2}=l_2.$

对频率循环，ω＝iΔωω的单位是弧度/秒，由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

对i·b₃取整，在数组a_i中拾取相邻两点的值，插值得到补偿因子λ_i，而波场实部则为

继续对频率循环，

此时补偿因子取为常量1.1A，继续由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

有

前面各式中Δω是频率采样，单位弧度/秒，τ是地震波走时，单位秒，

和

分别是步骤6处理后地震资料的实部和虚部。

在上述频率循环中记录峰值频率点对应的补偿因子λ_p和最大频率点i＝l₂对应的补偿因子

据此计算整数若m＞15，则令m＝15。

继续对频率循环，ω＝iΔω(i＝l₂+1，l₂+m)，继续由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

进而计算成像点的频率域波场实部为

在频率循环时同时求得波场的能量

所述的由各成像点频率域波场的实部，用变频带、保能量成像条件求得各成像点的偏移幅值是这样实现的：(1)令炮点水平坐标为x_s，检波点水平坐标为x_r，成像点坐标为(x₀，T)，成像点的偏移速度为v_rms；计算无量刚的成像权系数(2)对频率域波场的实部进行累加，计算

(3)从地震道的道头字中拾取E₀，再利用计算得到的波场能量E，最终得到成像点(x₀，T)的偏移幅值为I(x₀，T)ρ(E₀/E)。

本发明的补偿吸收衰减的叠前时间偏移方法，能恢复中、深层反射波中被衰减的高频成份，使得中、深层构造的成像分辨率达到与浅层接近的程度，并可保证补偿的稳定性且避免噪音放大。

本发明的补偿吸收衰减的叠前时间偏移方法，能基于地面采集的反射地震资料自主确定等效Q值，规避了现行各类衰减补偿方法在确定Q值上的困难。

本发明的具体实现原理如下：

本发明的核心有三点，一是从粘弹性单程波方程和稳相点原理出发，给出等效Q值的定义和基于等效Q值的地震波幅值吸收衰减的解析表达式；二是基于叠前深度偏移的成像条件，发展了稳定的、避免噪音放大的高频恢复成像算法；三是发展了利用地面采集的反射地震资料自主确定等效Q值的方法。其具体实现原理如下：

1.等效Q值与幅值衰减

首先基于粘弹性单程波方程和稳相点原理，给出等效Q值的定义和基于等效Q值的地震波幅值衰减的解析表达式。

采集是在二维剖面上进行的，将针对二维问题开展讨论。假设非均匀介质可近似为层状介质，不失一般性，令炮点或检波点放置在坐标原点。在波数-频率域，基于粘弹性单程波方程，炮点或检波点的波场深度延拓可表示为：

式中

是深度z处的波数-频率域波场，Δz_i是各层介质的厚度，n是目的层以上的层状介质的介质层数，ω是角频率，c_i(ω)是各层介质的复速度，k_x是水平波数，j是单位虚数，

是炮点或检波点的时域信号的傅立叶变换。对粘性吸收采用衡Q模型，并忽略实相速度随频率的微小变化，有

$\frac{\mathrm{\ω}}{c_i\left(\mathrm{\ω}\right)}=\frac{\mathrm{\ω}}{v_i}-j\frac{\mathrm{\ω}}{{2v}_i{Q}_i}---\left(2\right)$

式中v_i为各层介质的实相速度，Q_i为各层介质的品质因子。

用旅行时表示波场的延拓步长，即令Δz_i＝c_iΔτ_i，则深度z处的单程旅行时为

式(1)中的波数-频率域波场可进一步表示为：

式中p_x＝k_x/ω为射线参数。引入下式的复数开方近似式：

$\sqrt{1-\mathrm{jx}}\≈(1+x^2/8)-\mathrm{jx}/2---\left(4\right)$

式(3)中右端指数项中的相移量可近似为：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Δ\τ}}_i\sqrt{{(1-j\frac{1}{{2Q}_i})}^2-v_i^2{p}_x^2}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Δ\τ}}_i\sqrt{1-v_i^2{p}_x^2-\frac{1}{{4Q}_1^2}-\frac{j}{Q_i}}---\left(5\right)$

$\≈T\sqrt{1-v_{\mathrm{rms}}^2{p}_x^2}-j\frac{T}{{2Q}_{\mathrm{eff}}}\frac{1}{\sqrt{1-v_{\mathrm{Qrms}}^2{p}_x^2}}$

式中：

$v_{\mathrm{rms}}^2=\frac{1}{T}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{v}_i^2{\mathrm{\Δ\τ}}_i,$ $\frac{1}{Q_{\mathrm{eff}}}=\frac{1}{T}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\Δ\τ}}_i}{Q_i},$ $v_{\mathrm{Qrms}}^2=\frac{Q_{\mathrm{eff}}}{T}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{v_i^2{\mathrm{\Δ\τ}}_i}{Q_i}---\left(6\right)$

式(6)中第一个参数就是大家熟知的均方根速度，第二个参数是本发明引入的等效Q值，第三个参数是与Q值相关的均方根速度；为减少偏移计算中参数估计的工作量，一般假定第三个参数与均方根速度相等，即

将式(5)代入式(3)并做空间傅里叶反变换，可得空间-频率域波场为：

式中x是水平坐标，p_x是射线参数。式(7)是一个震荡积分，可利用稳相点原理求得渐进解为：

式(8)中

$\mathrm{\λ}\left(p_x\right)=\exp (-\frac{\mathrm{\ωT}}{{2Q}_{\mathrm{eff}}}\frac{1}{\sqrt{1-v_{\mathrm{rms}}^2{p}_x^2}})---\left(10\right)$

而

是的零点，可解得为：

$p_x^0=-\frac{x}{v_{\mathrm{rms}}\sqrt{x^2+T^2{v}_{\mathrm{rms}}^2}}---\left(11\right)$

将(11)式中的解代入式(8)，可得到地震波从坐标源点(0，0)到成像点(x，T)的走时项τ、几何扩散幅值项A和粘性吸收幅值项λ为：

$\mathrm{\τ}=T\sqrt{1+\frac{x^2}{{\left(v_{\mathrm{rms}}T\right)}^2}}---\left(12\right)$

$\mathrm{\λ}=\exp (-\frac{\mathrm{\ωT}}{{2Q}_{\mathrm{eff}}{c}_{\mathrm{rms}}}\sqrt{1+\frac{x^2}{{\left(v_{\mathrm{rms}}T\right)}^2}})=\exp (-\frac{\mathrm{\ω\τ}}{{2Q}_{\mathrm{eff}}})---\left(14\right)$

地震波传播的幅值因子为Aλ。当Q值为正无穷大时，λ＝1，此时式(8)退化为与常规叠前时间偏移完全相同。

尽管(12)、(13)和(14)式的走时和幅值是基于(1)式的层状介质假设推出的，但通过允许均方根速度v_rms和等效Q值参数Q_eff横向变化，式(12)、(13)和(14)的解可代表横向弱非均匀粘性介质中地震波的走时和幅值衰减。

2.成像算法

2.1成像条件

若将单个地震道看作是仅有一个接收道的单炮记录，则由(8)、(13)和(14)式的解析公式分别得到炮域偏移的下行波场和反传波场为：

和

式中

是震源子波的傅立叶变换，

是地震道的傅立叶变换，τ_s和τ_g以及A_s和A_g分别是(12)和(13)式表达的炮点至成像点和检波点至成像点的走时和几何扩散幅值。波场反传时幅值的衰减变为增加。

一般情况下震源子波是未知的。既然反褶积可剔除子波的影响，我们可在成像中忽略震源子波相关项的影响，即

应用反褶积成像条件，可得单道数据的成像结果(即脉冲响应)为

式中A_g/A_s是成像权系数，可利用(13)式的幅值公式得到。但(13)式的幅值是基于地震波在二维空间中传播得到的。本发明虽然针对二维测线，但由于震源仍是点源，地震波的几何扩散项需考虑地震波在三维空间中的传播，利用三维问题的幅值公式，可得成像权系数为 $\mathrm{\ρ}=\frac{{\left(v_{\mathrm{rms}}T\right)}^2+{(x_s-x_0)}^2}{{\left(v_{\mathrm{rms}}T\right)}^2+{(x_r-x_0)}^2}.$

若令地震资料有效频带的下和上限分别是l₁Δω和l₂Δω。式(15)可近似为：

式中Re代表取实部。

式(15)中的

是补偿地震波幅值衰减的补偿因子；忽略粘性，即当Q_eff趋于无穷大，该值为1，式(15)退化为与常规叠前时间偏移方法相同的形式。

当Q_eff不随时间深度和水平位置变化时，若假设

则式(15)可写为

式中g(τ_s+τ_g，τ_s+τ_g)就是对地震道f(t)做反Q滤波后的结果。式(17)表明，均匀Q值时，成像可通过对反Q滤波后的数据做常规叠前时间偏移来实现，即先做半导数，然后拾取τ_s+τ_r时刻的值，再乘权系数。因此。“反Q滤波+叠前时间偏移”是本发明补偿吸收衰减的偏移方法在均匀Q值情况下的特例。

2.2稳定性处理

当τ_s+τ_r较大或Q_eff较小时，式(15)存在稳定性问题，即补偿因子趋于无穷大；稳定性是各类补偿幅值衰减方法必须克服的问题。本发明为此提出了一个光滑性阈值控制的稳定性方法。具体实现如下：

令地震记录满足快速傅立叶算法的时间采样点数为M，Δτ是地震资料的时间采样，则角频率采样Δω＝2π/(MΔτ)(弧度/秒)；令地震资料有效频带的下和上限分别是l₁Δω和l₂Δω。根据地震资料的信噪比，确定成像过程中幅值补偿的最大倍数A。设计一个光滑函数，令在幅值小于等于A时，即当时，其与补偿因子

完全一致；当

时，其幅值光滑过渡到1.1A，然后保持恒定。

为避免重复设计这一函数，可事先计算相应的结果并将它们存在一个一维表中，每次按

的值到表中拾取对应数值。取表的间距为

其中Δq为拾取等效Q值时1/Q的间距。可按下式计算一维表中的元素：

a_i＝exp(iΔ)，i≤n₁                      (18)

a_i＝A(1-lnA-2.5(lnA)²)+A(1+5lnA)iΔ-2.5A(iΔ)²，n₁＜i≤n₂       (19)

上两式中

式(19)的函数保证了函数数值和一阶导数连续。当表的间距固定时，根据最高频率和最小等效Q值(Q_eff)_min，可得到拾取数值的最大下标为

因此，在构建表时，若n₁≥n₃，则n₁＝n₃和n₂＝n₃，若n₂≥n₃，则n₂＝n₃。

实际成像时，可根据计算得到的炮点至成像点再反射到检波点的走时τ和频率iΔω，对

取整；若该整数大于n₂，取λ＝1.1A；否则，在一维表中拾取相邻两点的值，插值得到λ_i，用λ_i做补偿因子。这样即避免了指数运算，又可保证补偿因子小于等于1.1A，保证了无条件稳定。

2.3噪音压制

尽管通过阈值控制可避免补偿因子过大导致的失稳，但较大的补偿因子在恢复幅值的同时也放大了噪音。本发明通过引入保能量和变频带成像条件，在恢复高频分量的同时避免了噪音放大。具体实现如下：

计算成像点处波场的能量

由式(20)可知，不进行粘性补偿时各成像点的能量就是地震道自身的能量，即是一常量；若对成像结果乘上因子(E₀/E)，则可实现能量不变。保能量的成像改变了频率成份的相对大小，因而实现了拓宽有效频带的目的。

成像是在频率域、按有限带宽完成的。为避免频率窗两端的折返效应导致的噪音，需对频率窗口的两端施加衰减带，使其平滑过渡到零。在偏移成像算法中，一般是通过对叠前资料做频率域衰减处理来实现的。但在补偿粘性衰减的成像中，随时间深度的增加，高频成份得到更大的补偿，频率窗高频端处原来接近零的值将变的远大于零，这导致高频噪音产生。为此本发明提出了一个随补偿因子大小变化的变频带成像方法，以此保证高频端平滑过渡到零。具体实现如下：

对叠前地震资料超出有效高频的部分，引入较宽的衰减带，计算

式中l₂Δω是最大有效频率，是求半导数后频率域的地震道。在成像过程的频率域累加计算中，高频端衰减带的厚度是根据补偿因子的相对大小而变化的，即衰减带所含点数为

$m=\mathrm{int}\left(\sqrt{\frac{1}{0.06}\ln ({\mathrm{\λ}}_{l_2}/{\mathrm{\λ}}_p)+25}\right)---\left(22\right)$

式中λ_p是峰值频率点处的补偿因子，

是最大频率点l₂处的补偿因子；式中用

来反应幅值补偿的大小。式(22)表明，最小的衰减带厚度为5个点，由式(21)可知，最大的衰减带厚度为15个点。由于处于衰减带中的频率成份仅是保证幅值平滑过渡到零，对处于衰减带中的频率成份，补偿因子将不随频率增加而增加，直接取为

2.4空道和坏道剔除

保能量的成像条件增加了除去地震道的能量。当存在空道时，这一除法运算将产生问题。为避免这一问题的发生，必须更严格地剔除空道和坏道。为此本发明也给出了更好地剔除空道和坏道的方法。具体实现如下：

地震资料是按偏移距排序分配到各个计算节点的。对每一组共偏移距地震道，设定单个地震道能量的变化间距，按E₀的大小将地震道分为若干组，一般应小于100。计算每组中平均的地震道道数，若最大和最小能量组中的道数小于平均数的百分之一，可认为该组地震道均是奇异道，可予以剔除；从最大和最小能量的两端剔除奇异道组不限于一次，可连续进行，至到最大和最小能量组中的道数超过平均数的百分之一停止。

2.5稳定的压噪反Q滤波

反Q滤波的幅值补偿算法是上述补偿粘性衰减的成像方法的特例，因此，可将保持稳定和压制噪音的成像方法直接应用于本发明的反Q滤波计算中。由于反Q滤波对应着均匀Q值，此时，补偿因子表的间距可直接取为其中Q_i是反Q滤波使用的Q值，这样拾取补偿因子时就避免了取整和插值。

3.等效Q值拾取

估计地球介质的Q值是一个非常困难的问题，特别是在Q值非均匀的境况下。本发明为此提出了等效Q值的概念，由于每一成像点的吸收补偿由该成像点的等效Q值唯一决定，使得可利用扫描方法直接拾取等效Q值，这极大地减少了非均匀Q值估计的困难。

在利用透射信号，如VSP资料，决定Q值时，可根据主频移动、频谱形状等信息决定Q值；但由于薄层调谐等现象存在，基于反射地震资料决定Q值则变的更加困难。主要原因是此时主频的移动以及频谱形状的变化受到粘性吸收和薄层调谐的共同影响，而很多时候薄层调谐的影响更大。本发明因此提出了一套完整的基于反射地震资料的等效Q值拾取方法。

在利用扫描方法直接拾取等效Q值的均匀Q值偏移计算中，如在2.1节中讨论的那样，我们利用了“反Q滤波+叠前时间偏移”来完成偏移计算，这一方法可高效地得到任一成像点随不同Q值变化的叠前成像道集。

在基于不同Q值的叠前成像道集拾取准确Q值方面，我们发展了基于谱比值的频率导数指标拾取等效Q值的方法。

不失一般性，下面将以均匀Q值情况，参考层是单个界面反射，而紧接着的反射事件是单个薄层反射为例讨论这一方法。就层状介质而言，由Snell定理可知，无论各层介质速度如何变化，射线参数p将保持不变。因此，在射线参数域，可将地震波入射和反射作为一平面波问题讨论。

假设平面波到参考层的双程走时为t₁，到达单个薄层的双程走时为t₂，波在薄层内传播的双程走时为Δτ，则地面记录的参考层的频率域反射响应为

式中r₁(p)是参考层的反射系数，是地震子波的频谱，p是射线参数。地面记录的薄层的频率域反射响应为

设反Q滤波使用的Q值为Q_i，则在成像道集中，参考层的频率域响应为

薄层的频率域响应为

对式(25)和(26)的频率域反射响应分别求模，然后求其对数并相减，最后对相减结果求频率导数，可得

$B(\mathrm{\ω},p)=\frac{t_2-t_1}{2}(\frac{1}{Q_i}-\frac{1}{Q})+\frac{\mathrm{\Δ\τ}}{2\tan (\mathrm{\ω\Δ\τ}/2)}---\left(27\right)$

式中第二项的零点在ω＝π/Δτ，3π/Δτ，5π/Δτ，在零点两端数值是反对称的；这样，沿频率轴对B(ω，p)求和，第二项的和将为零，若Q_i等于实际的Q值，则频率导数B(ω，p)的算术平均值将为零。因此，频率导数的平均值是否为零，可作为拾取Q值的重要指标；如式(27)所示，这一方法有效地克服了薄层调谐的影响。

式(26)在ω＝0，2π/Δτ，4π/Δτ处存在零点，为保证求对数的稳定性，我们发展了稳定的求对数算法。这一方法的核心是利用级数展开近似式：

$\ln x\≈(x-1)-\frac{1}{2}{(x-1)}^2+\frac{1}{3}{(x-1)}^3---\left(28\right)$

令D₀是有效频带内模D(ω)的平均值，当D(ω)≥0.22D₀，直接计算其对数ln[D(ω)]；当D(ω)＜0.22D₀时，令d＝5D(ω)/D₀，此时d＜1.1，可利用式(28)近似计算

$\ln \left[D\left(\mathrm{\ω}\right)\right]=\ln \left({0.2D}_0\right)+(d-1)(1-\frac{1}{2}(d-1)+\frac{1}{3}{(d-1)}^2)---\left(29\right)$

式(29)的稳定算法也避免了式(27)中分母为零的问题。

当有效频率范围不能包含式(27)中第二项在其零点两端的全部区域时，在计算频率导数的算术平均值时第二项的贡献就不是零。为此，在求算术平均值时可加以判断，若有效频带高频端出现没有相反符号峰值对应的峰值，则将该峰值段剔除。此外，式(27)中第二项对频率导数平均值的贡献随射线参数的变化改变较小；因此，当拾取到正确的Q值时，不同射线参数下的频率导数的平均值相差不大。这一分析表明，在根据频率导数拾取Q值时，即要导数值接近零，又要考虑不同射线参数间导数值的接近程度。

本发明的有益效果：该方法能在偏移过程中恢复地震波被衰减的高频成份，使得中、深层构造的成像分辨率达到与浅层接近的程度，从而获得更高分辨率的地下构造图像，实现更好地识别小断层和薄储层。该方法能基于地面采集的反射地震资料自主确定等效Q值，规避了现行各类吸收补偿方法在确定Q值上的困难。该方法对陆相薄互层油气储层的勘探和开发有重要应用价值。

附图说明

图1是南海某区块典型的共中心点道集，CDP点号为256。

图2是偏移速度场的等值线图。图中数字是偏移速度值。

图3是CDP256点邻近经增益处理后的偏移叠加剖面和选定的幅值对比点。图中箭头指示从浅到深的幅值对比点位置。

图4是CDP256点的4个不同时间深度上幅值对比点的小射线参数对应的短时Q道集(不包括参考点)，图中时间轴是各幅值对比点处选取的时窗，横轴是1/Q×10000。

图5是CDP256点的频率导数道集图。图中粗实线是由该图拾取的随时间深度变化的等效Q值，相邻细点线间即是频率导数值变化的窗口，细实线是小射线参数Q道集中各道的频率导数值随Q值的变化曲线，×号代表中射线参数Q道集各道的频率导数值随Q值的变化，·号代表大射线参数Q道集各道的频率导数值随Q值的变化。

图6是拾取得到的等效Q值场的等值线图。图中数字是1/Q值的大小。

图7是南海某区块应用补偿吸收衰减的叠前时间偏移方法得到的偏移成像结果。图中可见中、深层构造得到较清晰的成像，陡倾角构造清晰成像。

图8a是该区块应用常规叠前时间偏移方法得到的中深层局部成像结果，图8b是补偿吸收衰减的叠前时间偏移方法在相同局部区域的偏移成像结果。对比图8a和8b可知，补偿吸收衰减偏移方法使得同相轴变细，断点更清晰，成像分辨率明显提高。

图9a是该区块应用常规叠前时间偏移方法得到的深层局部成像结果，图9b是补偿吸收衰减的叠前时间偏移方法在相同局部区域的偏移成像结果。对比图9a和9b可知，补偿吸收衰减偏移方法使得深层的同相轴变的更细，断点也更清晰，成像分辨率明显提高，陡倾角构造得到更好成像。

具体实施方式

实施例1：补偿吸收衰减的叠前时间偏移方法，针对南海某区块为例，具体为以下步骤：

(1)用单条拖缆记录人工震源激发的反射地震信号，记录到磁带上。具体反射地震信号是，采用单条拖缆采集，每条拖缆930道，道间距3.125m；在拖缆的前端用气枪阵震源激发，气枪阵与第一个接收道的距离为130m，记录长度5s，时间采样1ms。共激发和记录666炮，炮点间距12.5m。

(2)从磁带上读取地震信号，对叠前地震资料做常规的噪音衰减处理。针对部分共中心点，抽取共中心点道集，对抽取的道集作常规的动校正速度拾取，对所得到速度的倒数做横向平均，将得到的横向均匀的速度场作为初始偏移速度。图1是典型的共中心点道集。

(3)将噪音衰减处理后的叠前地震资料按偏移距排序、分组，将不同组地震资料放置到集群计算机的不同计算节点上。在每个计算节点，依据初始偏移速度，用常规的叠前时间方法进行偏移计算，收集各计算节点的偏移结果，形成共反射点道集。具体是，将地震资料按偏移距大小分为35组，放置到集群计算机的35个节点上。

(4)对共反射点道集，利用初始偏移速度做反动校，再做动校正得到更新后的速度，对各个CDP点更新后的速度的倒数做空间平滑，所得到的平滑速度场即是最终的偏移速度场。图2是偏移速度场的等值线图。图中数字是偏移速度值。

(5)依据偏移速度场，对噪音衰减处理后的叠前地震资料做常规叠前时间偏移计算。收集各计算节点的偏移结果，形成共反射点道集，对共反射点道集做剩余动校和拉伸切除，得到偏移叠加剖面。在经增益处理后的偏移叠加剖面上选定清晰、连续的同相轴。根据给定的采样间距，确定等间距分布的拾取Q值的CDP点，依据选定的同相轴，在这些CDP点上确定从浅到深的幅值对比点和浅层的参考点，并在各点处选定应用短时傅立叶变换的时窗。从浅到深的幅值对比点处的同相轴需展现由细变粗的变化趋势。具体是，拾取Q值的CDP点的采样间距是10个CDP点。图3是CDP256点邻近经增益处理后的偏移叠加剖面和选定的幅值对比点；图中箭头指示从浅到深的幅值对比点位置。

(6)对噪音衰减处理后的叠前地震资料应用时间域快速傅立叶变换，求半导数、施加频率窗衰减带、计算地震道能量、剔除奇异道。

在每个计算节点，对全部地震道循环，应用快速傅立叶变换算法，对地震道的时间序列离散信号做傅立叶变换；做傅立叶变换时，根据时间采样点数，确定满足快速傅立叶算法的采样点数M，增加采样点时添加零值。设Δτ是地震资料的时间采样，单位秒，则角频率采样Δω＝2π/(MΔτ)，其单位是弧度/秒；令地震资料有效频带的下和上限分别是ω_min和ω_max，其单位是弧度/秒，通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近。

设

是频率域的地震道，采用拖缆记录时其单位是帕，采用测线记录时其单位是米/秒；对i＝l₁，l₂循环，求半导数后的地震道为

式中j是单位虚数，Δω是角频率采样，单位弧度/秒。

设频率窗低频端的衰减带含m_d个频率采样点，对i＝l₁，l₁+m_d，计算

在频率窗的高频端，对i＝l₂+1，l₂+15，计算

对每一地震道，将存储(l₂+15)-l₁+1个复数值。

计算地震道的能量

将E₀记入该地震道的道头字中。

完成对所有地震道的循环，即得到计算半导数和施加频率窗衰减带后的频率域叠前地震资料。

地震资料是按偏移距排序分配到各个计算节点的。对每一组共偏移距地震道，设定单个地震道能量的变化间距，按E₀的大小将地震道分为100组。计算每组中平均的地震道道数，若最大和最小能量组中的道数小于平均数的百分之一，则判定该组地震道均是奇异道，予以剔除；从最大和最小能量的两端剔除奇异道组不限于一次，连续进行，至到最大和最小能量组中的道数超过平均数的百分之一停止。

(7)确定等效Q值的取值范围，按1/Q等间距选取系列Q值，用每个Q值，对步骤6处理后的频率域叠前地震资料，做稳定的压噪反Q滤波计算，得到新的时间域地震资料。具体是，令等效Q值的取值范围为100到500，共选取36个系列Q值。

设等效Q值的最小、最大可能值为Q_min、Q_max，等效Q值是无量纲参数，则选取的N个系列Q值为

$\frac{1}{Q_i}=\frac{1}{Q_{\max }}+\frac{i-1}{N-1}(\frac{1}{Q_{\min }}-\frac{1}{Q_{\max }}),i=1,N$

令地震记录满足快速傅立叶算法的时间采样点数为M，Δτ是地震资料的时间采样，单位秒，则角频率采样Δω＝2π/(MΔτ)，其单位是弧度/秒；令地震资料有效频带的下和上限分别是ω_min和ω_max，其单位是弧度/秒，通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近；根据地震资料的信噪比，确定幅值补偿的最大倍数A，一般取为2000.0，A是无量纲数。

对每个Q_i，计算

n₁＝lnA/Δ，n₂＝(lnA+0.2)/Δ，n₃＝l₂M；若n₁≥n₃，则n₁＝n₃和n₂＝n₃，若n₂≥n₃，则n₂＝n₃。建立数组a_i(i＝l₁，n₃)，有

a_i＝exp(iΔ)，i≤n₁

a_i＝A(1-lnA-2.5(lnA)²)+A(1+5lnA)iΔ-2.5A(iΔ)²，n₁＜i≤n₂

a_i＝1.1A，n₂＜i≤n₂

对每一地震道，在对该道做反Q滤波前，事先计算地震道的模

对时间采样点循环，τ＝mΔτ(m＝1，M)，τ的单位是秒，计算

d₀＝cos(Δω·τ)，b₀＝sin(Δω·τ)

d₁＝cos(l₁Δω·τ)，b₁＝sin(l₁Δω·τ)

对频率循环，ω＝iΔω(i＝l₁，l₂)，递推计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

计算n＝m·i，在数组a_i中拾取第n个元素值a_n，得波场实部为

式中

和

分别是步骤6处理后地震资料的实部和虚部。

在上述频率循环中记录峰值频率点对应的补偿因子a_p和最大频率点对应的

据此计算整数 $k=\mathrm{int}\left(\sqrt{\frac{1}{0.06}\ln (a_{(m\·l_2)}/a_p)+25}\right);$ 若k＞15，则令k＝15。

继续对频率循环，ω＝iΔω(i＝l₂+1，l₂+k)，同上由递推法计算exp(j·iΔω·τ)的实部di和虚部bi为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

进而计算

进行累加计算

从该地震道的道头字中拾取E₀，计算

P(τ)＝(E₀/E_m)P₀(τ)

完成全部时间采样点循环，得到的P(τ)就是新的时间域地震资料的一个地震道，其单位与原始地震道相同。在各个计算节点对所有地震道循环，就得到对应Q_i的新的时间域地震资料。

(8)依据偏移速度场，对每个Q值对应的反Q滤波后的地震资料做常规叠前时间偏移计算。收集各计算节点的偏移结果，形成共反射点道集，对共反射点道集做剩余动校，切除剩余动校后的拉伸部分。将需拾取Q值的CDP点的相邻共反射点道集叠加，形成各CDP点的超道集。定义小、中、大三个射线参数p₀、p₁、p₂，对全部需拾取Q值的CDP点的幅值对比点和参考点计算偏移距

其中i＝0，1，2，T为该点的单程旅行时，单位秒，v_rms为该点的偏移速度，单位米/秒，射线参数p_i的单位是秒/米，偏移距d_i的单位是米。在各点选定的时窗内，按偏移距范围[0，d₀]、[d₀，d₁]和[d₁，d₂]对超道集分别进行叠加，在每点处得到三个短时时间序列。具体是，定义小、中、大三个射线参数为2.5e-4，3.5e-4，5.0e-4。

(9)重复步骤7和步骤8，收集各幅值对比点和参考点的对应小、中、大射线参数的三个短时时间序列，按1/Q值大小排列，在各点形成三个随1/Q值变化的小、中、大射线参数的短时Q道集。图4是CDP256点的4个不同时间深度上幅值对比点的小射线参数对应的短时Q道集，图中时间轴是各幅值对比点处选取的时窗。

(10)根据小、中、大射线参数的短时Q道集拾取等效Q值，得到成像区域的等效Q值场。图5是CDP256点的频率导数道集图。图6是拾取得到的等效Q值场的等值线图，图6中数字是1/Q值的大小。

采用以下7个分步骤：1)对短时Q道集中的各道做傅立叶变换，得到频率信号；2)计算傅立叶变换的模，在有效频带内，求模的自然对数。定义有效频带为

其中m_d为低频端的衰减带所含的频率采样点数，

由各短时Q道集对应的射线参数p_i，单位秒/米，对应幅值对比点的单程旅行时T，单位秒，偏移速度v_rms，单位米/秒，各道对应的Q_i，以及频带上限l₂Δω得到：

${\stackrel{\‾}{\mathrm{\ω}}}_{\max }=\min \{\frac{Q_i\ln A}{T}\sqrt{1-v_{\mathrm{rms}}^2{p}_i^2},l_2\mathrm{\Δ\ω}\}.$

为计算的稳定性，模的自然对数由如下方法计算：令D₀是有效频带内模D(ω)的平均值，当D(ω)≥0.22D₀，直接计算其对数ln[D(ω)]；当D(ω)＜0.22D₀时，令d＝5D(ω)/D₀，计算

$\ln \left[D\left(\mathrm{\ω}\right)\right]=\ln \left({0.2D}_0\right)+(d-1)(1-\frac{1}{2}(d-1)+\frac{1}{3}{(d-1)}^2)$

3)对计算得到的对数做中值滤波处理；4)收集同属于同一CDP点的短时Q道集的随频率变化的对数值，将各幅值对比点的各道的对数值与浅层的参考点处对应道的对数值相减。在有效频带内，对随频率变化的相减结果先做高斯平滑，然后用中心差分法计算其对频率的一阶导数。对得到的导数值做中值滤波，然后求导数的算术平均值，平均值的绝对值就是该道的频率导数值。在求算术平均值时，若有效频带高频端出现没有相反符号峰值对应的峰值，应将该峰值段剔除；5)形成需拾取Q值的全部CDP点的频率导数道集图。该图做法如下：在各个幅值对比点的时间深度上，在指定的窗口内，将小射线参数Q道集各道的频率导数值做成随Q值变化的曲线，把中、大射线参数Q道集各道的频率导数值标在图上对应的位置上；6)基于这一图形，在导数值接近零的邻域，考虑小、中、大射线参数结果间的接近程度，直接拾取等效Q值。拾取时要同时考虑等效Q值的横向连续性并保证等效Q值应随旅行时增加而增加；7)对拾取的Q值，做有关1/Q的空间平滑，即得到成像区域的等效Q值场。

(11)在每个计算节点，对步骤6处理后的地震资料的全部地震道循环，应用稳定的变频带频率域算法计算各成像点的频率域波场实部。

令地震记录满足快速傅立叶算法的时间采样点数为M，Δτ是地震资料的时间采样，单位秒，则角频率采样Δω＝2π/(MΔτ)，其单位是弧度/秒；令地震资料有效频带的下和上限分别是ω_min和ω_max，其单位是弧度/秒，通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近；根据地震资料的信噪比，确定幅值补偿的最大倍数A，一般取为2000.0，A是无量纲数；设拾取等效Q值时采用的1/Q间距为 $\mathrm{\Δq}=\frac{1}{N-1}(\frac{1}{Q_{\min }}-\frac{1}{Q_{\max }}),$ Δq是无量纲数，计算

其中(Q_eff)_min是由Q值场确定的成像区域内最小的等效Q值；若n₁≥n₃，则n₁＝n₃和n₂＝n₃，若n₂≥n₃，则n₂＝n₃；建立数组a_i(i＝1，n₂)，有

a_i＝exp(iΔ)，i≤n₁

a_i＝A(1-lnA-2.5(lnA)²)+A(1+5lnA)iΔ-2.5A(iΔ)²，n₁＜i≤n₂

对每一地震道，在对该道成像前，事先计算地震道的模

采集是在二维剖面上进行的。令炮点水平坐标为x_s，单位米，检波点水平坐标为x_r，单位米，成像点的水平坐标和单程旅行时分别为x₀和T，其单位分别是米和秒，成像点处的偏移速度为v_rms，单位米/秒，成像点处的等效Q值为Q_eff；则炮点至成像点再反射到检波点的走时为

$\mathrm{\τ}=T(\sqrt{1+\frac{{(x_s-x_0)}^2}{{\left(v_{\mathrm{rms}}T\right)}^2}}+\sqrt{1+\frac{{(x_r-x_0)}^2}{{\left(v_{\mathrm{rms}}T\right)}^2}})$

成像仅需要频率域波场的实部，成像点(x₀，T)处的频率域波场实部由如下快速方法求得：

计算

d₀＝cos(Δω·τ)，b₀＝sin(Δω·τ)

d₁＝cos(l₁Δω·τ)，b₁＝sin(l₁Δω·τ)

$b_3=\frac{\mathrm{\Δ\ω}\·\mathrm{\τ}}{{2Q}_{\mathrm{eff}}\mathrm{\Δ}}$

$\stackrel{\‾}{l_2}=\mathrm{int}(n_2/b_3)$

若 $\stackrel{\‾}{l_2}>l_2,$ 则 $\stackrel{\‾}{l_2}=l_2.$

对频率循环，ω＝iΔω

ω的单位是弧度/秒，由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

对i·b₃取整，在数组a_i中拾取相邻两点的值，插值得到补偿因子λ_i，而波场实部则为

继续对频率循环，

此时补偿因子取为常量1.1A，继续由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

有

前面各式中Δω是频率采样，单位弧度/秒，τ是地震波走时，单位秒，

和

分别是步骤6处理后地震资料的实部和虚部。

在上述频率循环中记录峰值频率点对应的补偿因子λ_p和最大频率点i＝l₂对应的补偿因子

据此计算整数

若m＞15，则令m＝15。

继续对频率循环，ω＝iΔω(i＝l₂+1，l₂+m)，继续由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

进而计算成像点的频率域波场实部为

在频率循环时同时求得波场的能量

(12)由各成像点频率域波场的实部，用变频带、保能量成像条件求得各成像点的偏移幅值。

采用以下3个分步骤：1)令炮点水平坐标为x_s，检波点水平坐标为x_r，成像点坐标为(x₀，T)，成像点的偏移速度为v_rms；计算无量刚的成像权系数

2)对频率域波场的实部进行累加，计算

3)从地震道的道头字中拾取E₀，再利用计算得到的波场能量E，最终得到成像点(x₀，T)的偏移幅值为I(x₀，T)ρ(E₀/E)。

(13)将偏移幅值按偏移距大小累加到存放偏移结果的数组中对应的偏移距上。收集各计算节点的偏移结果，形成全部CDP点的共反射点道集。对该道集做剩余动校正和拉伸切除，将不同偏移距的偏移结果叠加，形成偏移叠加剖面。

(14)通过显示软件将偏移叠加剖面数值转换为地下反射构造的剖面图像，剖面图像将指示中、深层构造的微小断裂体系和精细的地层沉积样式，用于确定地下生、储油构造和识别油气储层。图7是南海某区块应用补偿吸收衰减的叠前时间偏移方法得到的偏移成像结果。图7中可见中、深层构造得到较清晰的成像，陡倾角构造清晰成像。图8和图9是该区块同样地震资料应用常规叠前时间偏移方法和补偿吸收衰减的叠前时间偏移方法的偏移成像结果的局部对比图。由图8和图9可知，中、深层的同相轴变细，断点更清晰，成像分辨率明显提高，陡倾角构造得到更好成像。

Claims (6)

1.一种补偿吸收衰减的叠前时间偏移方法，其特征在于：包括以下步骤：步骤A、用单条拖缆或测线记录人工震源激发的反射地震信号，记录到磁带上；步骤B、从磁带上读取地震信号，对叠前地震资料做常规的噪音衰减处理；针对部分共中心点，抽取共中心点道集，对抽取的道集作常规的动校正速度拾取，对所得到速度的倒数做横向平均，将得到的横向均匀的速度场作为初始偏移速度；步骤C、将噪音衰减处理后的叠前地震资料按偏移距排序、分组，将不同组地震资料放置到集群计算机的不同计算节点上；在每个计算节点，依据初始偏移速度，用常规的叠前时间方法进行偏移计算，收集各计算节点的偏移结果，形成共反射点道集；步骤D、对共反射点道集，利用初始偏移速度做反动校，再做动校正得到更新后的速度，对各个CDP点更新后的速度的倒数做空间平滑，所得到的平滑速度场即是最终的偏移速度场；步骤E、依据偏移速度场，对噪音衰减处理后的叠前地震资料做常规叠前时间偏移计算；收集各计算节点的偏移结果，形成共反射点道集，对共反射点道集做剩余动校和拉伸切除，得到偏移叠加剖面；在经增益处理后的偏移叠加剖面上选定清晰、连续的同相轴；根据给定的采样间距，确定等间距分布的拾取Q值的CDP点，依据选定的同相轴，在这些CDP点上确定从浅到深的幅值对比点和浅层的参考点，并在各点处选定应用短时傅立叶变换的时窗；从浅到深的幅值对比点处的同相轴需展现由细变粗的变化趋势；步骤F、对噪音衰减处理后的叠前地震资料应用时间域快速傅立叶变换，求半导数、施加频率窗衰减带、计算地震道能量、剔除奇异道；步骤G、确定等效Q值的取值范围，按1/Q等间距选取系列Q值，用每个Q值，对步骤F处理后的频率域叠前地震资料，做稳定的压噪反Q滤波计算，得到新的时间域地震资料；步骤H、依据偏移速度场，对每个Q值对应的反Q滤波后的地震资料做常规叠前时间偏移计算；收集各计算节点的偏移结果，形成共反射点道集，对共反射点道集做剩余动校，切除剩余动校后的拉伸部分；将需拾取Q值的CDP点的相邻共反射点道集叠加，形成各CDP点的超道集；定义小、中、大三个射线参数p₀、p₁、p₂，对全部需拾取Q值的CDP点的幅值对比点和参考点计算偏移距

其中i＝0，1，2，T为该点的单程旅行时，单位秒，v_rms为该点的偏移速度，单位米/秒，射线参数p_i的单位是秒/米，偏移距d_i的单位是米；在各点选定的时窗内，按偏移距范围[0，d₀]、[d₀，d₁]和[d₁，d₂]对超道集分别进行叠加，在每点处得到三个短时时间序列；步骤I、重复步骤G和步骤H，收集各幅值对比点和参考点的对应小、中、大射线参数的三个短时时间序列，按1/Q值大小排列，在各点形成三个随1/Q值变化的小、中、大射线参数的短时Q道集；步骤J、根据小、中、大射线参数的短时Q道集拾取等效Q值，得到成像区域的等效Q值场；步骤K、在每个计算节点，对步骤F处理后的地震资料的全部地震道循环，应用稳定的变频带频率域算法计算各成像点的频率域波场实部；步骤L、由各成像点频率域波场的实部，用变频带、保能量成像条件求得各成像点的偏移幅值；步骤M、将偏移幅值按偏移距大小累加到存放偏移结果的数组中对应的偏移距上；收集各计算节点的偏移结果，形成全部CDP点的共反射点道集；对该道集做剩余动校正和拉伸切除，将不同偏移距的偏移结果叠加，形成偏移叠加剖面；步骤N、通过显示软件将偏移叠加剖面数值转换为地下反射构造的剖面图像，剖面图像将指示中、深层构造的微小断裂体系和精细的地层沉积样式，用于确定地下生、储油构造和识别油气储层。

2.根据权力要求1所述的一种补偿吸收衰减的叠前时间偏移方法，其特征在于：所述的对噪音衰减处理后的叠前地震资料应用时间域快速傅立叶变换，求半导数、施加频率窗衰减带、计算地震道能量、剔除奇异道是这样实现的：在每个计算节点，对全部地震道循环，应用快速傅立叶变换算法，对地震道的时间序列离散信号做傅立叶变换；做傅立叶变换时，根据时间采样点数，确定满足快速傅立叶算法的采样点数M，增加采样点时添加零值；设Δτ是地震资料的时间采样，单位秒，则角频率采样Δω＝2π/(MΔτ)，其单位是弧度/秒；令地震资料有效频带的下和上限分别是ω_min和ω_max，其单位是弧度/秒，通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近；设

是频率域的地震道，采用拖缆记录时其单位是帕，采用测线记录时其单位是米/秒；对i＝l₁，l₂循环，求半导数后的地震道为

式中j是单位虚数，Δω是角频率采样，单位弧度/秒；

设频率窗低频端的衰减带含m_d个频率采样点，对i＝l₁，l₁+m_d，计算

在频率窗的高频端，对i＝l₂+1，l₂+15，计算

对每一地震道，将存储(l₂+15)-l₁+1个复数值；

计算地震道的能量

将E₀记入该地震道的道头字中；

完成对所有地震道的循环，即得到计算半导数和施加频率窗衰减带后的频率域叠前地震资料；

地震资料是按偏移距排序分配到各个计算节点的；对每一组共偏移距地震道，设定单个地震道能量的变化间距，按E₀的大小将地震道分为若干组，一般取50～100组；计算每组中平均的地震道道数，若最大和最小能量组中的道数小于平均数的百分之一，则判定该组地震道均是奇异道，予以剔除；从最大和最小能量的两端剔除奇异道组不限于一次，连续进行，至到最大和最小能量组中的道数超过平均数的百分之一停止。

3.根据权力要求1所述的一种补偿吸收衰减的叠前时间偏移方法，其特征在于：所述的确定等效Q值的取值范围，按1/Q等间距选取系列Q值，用每个Q值，对步骤F处理后的频率域叠前地震资料，做稳定的压噪反Q滤波计算，得到新的时间域地震资料是这样实现的：设等效Q值的最小、最大可能值为Q_min、Q_max，等效Q值是无量纲参数，则选取的N个系列Q值为

$\frac{1}{Q_i}=\frac{1}{Q_{\max }}+\frac{i-1}{N-1}(\frac{1}{Q_{\min }}-\frac{1}{Q_{\max }}),i=1,N$

令地震记录满足快速傅立叶算法的时间采样点数为M，Δτ是地震资料的时间采样，单位秒，则角频率采样Δω＝2π/(MΔτ)，其单位是弧度/秒；令地震资料有效频带的下和上限分别是ω_min和ω_max，其单位是弧度/秒，通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近；根据地震资料的信噪比，确定幅值补偿的最大倍数A，一般取为2000.0，A是无量纲数；

对每个Q_i，计算

n₁＝lnA/Δ，n₂＝(lnA+0.2)/Δ，n₃＝l₂M；若n₁≥n₂，则n₁＝n₂和n₂＝n₃，若n₂≥n₃，则n₂＝n₃；建立数组a_i(i＝l₁，n₃)，有

a_i＝exp(iΔ)，i≤n₁

a_i＝A(1-lnA-2.5(lnA)²)+A(1+5lnA)iΔ-2.5A(iΔ)²，n₁＜i≤n₂

a_i＝1.1A，n₂＜i≤n₃

对每一地震道，在对该道做反Q滤波前，事先计算地震道的模

对时间采样点循环，τ＝mΔτ(m＝1，M)，τ的单位是秒，计算

d₀＝cos(Δω·τ)，b₀＝sin(Δω·τ)

d₁＝cos(l₁Δω·τ)，b₁＝sin(l₁Δω·τ)

对频率循环，ω＝iΔω(i＝l₁，l₂)，递推计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

计算n＝m·i，在数组a_i中拾取第n个元素值a_n，得波场实部为

式中

和

分别是步骤F处理后地震资料的实部和虚部；

在上述频率循环中记录峰值频率点对应的补偿因子a_p和最大频率点对应的据此计算整数 $k=\mathrm{int}\left(\sqrt{\frac{1}{0.06}\ln (a_{(m\·l_2)}/a_p)+25}\right);$ 若k＞15，则令k＝15；

继续对频率循环，ω＝iΔω(i＝l₂+1，l₂+k)，同上由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

进而计算

进行累加计算

从该地震道的道头字中拾取E₀，计算

P(τ)＝(E₀/E_m)P₀(τ)

完成全部时间采样点循环，得到的P(τ)就是新的时间域地震资料的一个地震道，其单位与原始地震道相同；在各个计算节点对所有地震道循环，就得到对应Q_i的新的时间域地震资料。

4.根据权力要求1所述的一种补偿吸收衰减的叠前时间偏移方法，其特征在于：所述的根据小、中、大射线参数的短时Q道集拾取等效Q值，得到成像区域的等效Q值场是这样实现的：(1)对短时Q道集中的各道做傅立叶变换，得到频率信号；(2)计算傅立叶变换的模，在有效频带内，求模的自然对数；定义有效频带为

其中m_d为低频端的衰减带所含的频率采样点数，

由各短时Q道集对应的射线参数p_i，单位秒/米，对应幅值对比点的单程旅行时T，单位秒，偏移速度v_rms，单位米/秒，各道对应的Q_i，以及频带上限l₂Δω得到： ${\stackrel{\‾}{\mathrm{\ω}}}_{\max }=\min \{\frac{Q_i\ln A}{T}\sqrt{1-v_{\mathrm{rms}}^2{p}_i^2},l_2\mathrm{\Δ\ω}\};$

为计算的稳定性，模的自然对数由如下方法计算：令D₀是有效频带内模D(ω)的平均值，当D(ω)≥0.22D₀，直接计算其对数ln[D(ω)]；当D(ω)＜0.22D₀时，令d＝5D(ω)/D₀，计算

$\ln \left[D\left(\mathrm{\ω}\right)\right]=\ln \left({0.2D}_0\right)+(d-1)(1-\frac{1}{2}(d-1)+\frac{1}{3}{(d-1)}^2)$

(3)对计算得到的对数做中值滤波处理；(4)收集同属于同一CDP点的短时Q道集的随频率变化的对数值，将各幅值对比点的各道的对数值与浅层的参考点处对应道的对数值相减；在有效频带内，对随频率变化的相减结果先做高斯平滑，然后用中心差分法计算其对频率的一阶导数；对得到的导数值做中值滤波，然后求导数的算术平均值，平均值的绝对值就是该道的频率导数值；在求算术平均值时，若有效频带高频端出现没有相反符号峰值对应的峰值，应将该峰值段剔除；(5)形成需拾取Q值的全部CDP点的频率导数道集图；该图做法如下：在各个幅值对比点的时间深度上，在指定的窗口内，将小射线参数Q道集各道的频率导数值做成随Q值变化的曲线，把中、大射线参数Q道集各道的频率导数值标在图上对应的位置上；(6)基于这一图形，在频率导数值接近零的邻域，考虑小、中、大射线参数结果间的接近程度，直接拾取等效Q值；拾取时要同时考虑等效Q值的横向连续性并保证等效Q值应随旅行时增加而增加；(7)对拾取的Q值，做有关1/Q的空间平滑，即得到成像区域的等效Q值场。

5.根据权力要求1所述的一种补偿吸收衰减的叠前时间偏移方法，其特征在于：所述的在每个计算节点，对步骤F处理后的地震资料的全部地震道循环，应用稳定的变频带频率域算法计算各成像点频率域波场实部是这样实现的：令地震记录满足快速傅立叶算法的时间采样点数为M，Δτ是地震资料的时间采样，单位秒，则角频率采样Δω＝2π/(MΔτ)，其单位是弧度/秒；令地震资料有效频带的下和上限分别是ω_min和ω_max，其单位是弧度/秒，通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近；根据地震资料的信噪比，确定幅值补偿的最大倍数A，一般取为2000.0，A是无量纲数；设拾取等效Q值时采用的1/Q间距为

Δq是无量纲数，计算

其中(Q_eff)_min是由Q值场确定的成像区域内最小的等效Q值；若n₁≥n₃，则n₁＝n₃和n₂＝n₃，若n₂≥n₃，则n₂＝n₃；建立数组a_i(i＝1，n₂)，有

a_i＝exp(iΔ)，i≤n₁

a_i＝A(1-lnA-2.5(lnA)²)+A(1+5lnA)iΔ-2.5A(iΔ)²，n₁＜i≤n₂

对每一地震道，在对该道成像前，事先计算地震道的模

采集是在二维剖面上进行的；令炮点水平坐标为x_s，单位米，检波点水平坐标为x_r，单位米，成像点的水平坐标和单程旅行时分别为x₀和T，其单位分别是米和秒，成像点处的偏移速度为v_rms，单位米/秒，成像点处的等效Q值为Q_eff；则炮点至成像点再反射到检波点的走时为

$\mathrm{\τ}=T(\sqrt{1+\frac{{(x_s-x_0)}^2}{{\left(v_{\mathrm{rms}}T\right)}^2}}+\sqrt{1+\frac{{(x_r-x_0)}^2}{{\left(v_{\mathrm{rms}}T\right)}^2}})$

成像仅需要频率域波场的实部，成像点(x₀，T)处的频率域波场实部由如下快速方法求得：

计算

d₀＝cos(Δω·τ)，b₀＝sin(Δω·τ)

d₁＝cos(l₁Δω·τ)，b₁＝sin(l₁Δω·τ)

$b_3=\frac{\mathrm{\Δ\ω}\·\mathrm{\τ}}{{2Q}_{\mathrm{eff}}\mathrm{\Δ}}$

$\stackrel{\‾}{l_2}=\mathrm{int}(n_2/b_3)$

若 $\stackrel{\‾}{l_2}>l_2,$ 则 $\stackrel{\‾}{l_2}=l_2;$

对频率循环，ω＝iΔωω的单位是弧度/秒，由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

对i·b₃取整，在数组a_i中拾取相邻两点的值，插值得到补偿因子λ_i，而波场实部则为

继续对频率循环，ω＝iΔω

此时补偿因子取为常量1.1A，继续由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

有

前面各式中Δω是频率采样，单位弧度/秒，τ是地震波走时，单位秒，

和

分别是步骤F处理后地震资料的实部和虚部；

在上述频率循环中记录峰值频率点对应的补偿因子λ_p和最大频率点i＝l₂对应的补偿因子据此计算整数

若m＞15，则令m＝15；

继续对频率循环，ω＝iΔω(i＝l₂+1，l₂+m)，继续由递推法计算exp(j·iΔω·τ)的实部d_i和虚部b_i为

d_i＝d_i-1d₀-b_i-1b₀，b_i＝b_i-1d₀+d_i-1b₀

进而计算成像点的频率域波场实部为

在频率循环时同时求得波场的能量

6.根据权力要求1所述的一种补偿吸收衰减的叠前时间偏移方法，其特征在于：所述的由各成像点频率域波场的实部，用变频带、保能量成像条件求得各成像点的偏移幅值是这样实现的：(1)令炮点水平坐标为x_s，检波点水平坐标为x_r，成像点坐标为(x₀，T)，成像点的偏移速度为v_rms；计算无量刚的成像权系数

(2)对频率域波场的实部进行累加，计算

(3)从地震道的道头字中拾取E₀，再利用计算得到的波场能量E，最终得到成像点(x₀，T)的偏移幅值为I(x₀，T)ρ(E₀/E)。

Images