CN103424777A - 一种提高地震成像分辨率的方法 - Google Patents

Abstract

一种提高地震成像分辨率的方法，应用于地震勘探中反射地震资料处理。该方法通过对偏移叠加数据体进行频率相关的幅值恢复和频散校正，整体提升成像的分辨率，并使得中、深层构造的成像分辨率达到与浅层接近的程度。该方法基于粘弹性理论进行幅值恢复和频散校正，但不需要事先提供非均匀Q值场。该方法引入等效Q值描述地震波的粘性吸收，通过建立单一的数值指标，利用参数扫描得到计算所需的非均匀等效Q值场。与非稳态反褶积等方法相比，所建议的方法在恢复高频信息的同时也完成频散校正；由于引入了稳定控制算法，该方法保证计算的稳定性。该方法能更好地指示中、深层构造的断裂和地层沉积样式，对油气、矿产资源勘探有重要应用价值。

Description

一种提高地震成像分辨率的方法

技术领域

本发明属于地震勘探中反射地震资料处理技术领域，涉及地震资料处理过程中的偏移成像技术范畴，是一种能提高地震成像分辨率的地震资料高分辨处理方法。

背景技术

在油气勘探中，对小断层、小断裂的识别，是认识油气疏导体系，进而识别有利储层的重要环节。因此，在地震勘探中反射地震资料处理过程中，提高地震成像的分辨率一直是一个重要的努力目标。

实际地球介质存在粘性吸收，地球介质的小尺度非均匀也会产生类似于粘性吸收的幅值衰减效应。这些客观存在导致地震波在传播过程中发生幅值的吸收衰减和随频率的传播速度变化；幅值的衰减对地震波的不同频率成份是不同的，频率越高，衰减就越强，这导致接收到的反射地震资料的有效频带随反射深度逐渐变窄；而不同频率成份以不同的速度传播，也导致了地震子波的频散，这一频散现象也是反射构造越深，频散越严重。由于常规偏移方法没有补偿粘性吸收导致的幅值衰减，也没有校正频散，因而偏移成像结果的分辨率较低，其特点是，构造越深，分辨率就越低。

为提高地震成像的分辨率，已发展了许多方法，可分为针对叠前资料和偏移叠加数据体的两类技术。针对叠前资料的方法从理论上讲更合理，但由于叠前资料信噪比较低，而三维叠前资料数据量巨大，大规模实际应用还存在一些困难。由于偏移叠加数据体的信噪比高，数据量较叠前资料减少了许多，针对这类数据的提高成像分辨率研究较活跃；已发展了谱白化反褶积，非稳态反褶积，反Q滤波和基于统计假设或测井资料的各类拓宽频带技术。谱白化以及各类拓频技术是通过引入地震资料以外的信息提升地震成像的分辨率，尽管可获得更高的视分辨率，但其可靠性尚待提高。

非稳态反褶积和反Q滤波是针对粘性吸收导致的分辨率降低而发展的提高分辨率方法，有较坚实的物理基础。但非稳态反褶积在估计空变的非稳态子波上存在较多困难，这一方法一般不能同时实现频散校正；由于很难获得非均匀的Q值，现行的反Q滤波方法更多地采用均匀Q值和层状Q值假定，这不能反应实际地球介质的非均匀Q值情况，也限制了反Q滤波的粘性补偿成效；稳定性和噪音放大是这两类方法实际应用遇到的另一类问题。

发明内容

本发明的目的是：提供一种提高地震成像分辨率的方法，它通过对偏移叠加数据体进行频率相关的幅值恢复和频散校正，整体提升成像的分辨率，并使得中、深层构造的成像分辨率达到与浅层接近的程度；这一方法基于粘弹性理论，通过使用非均匀Q进行粘性吸收补偿，来实现幅值恢复和频散校正，但不需要事先提供非均匀Q值场。

本发明采用的技术方案是：提高地震成像分辨率的方法，具体步骤包括：

（1）用多条拖缆或测线采集人工震源激发的反射地震信号，将地震信号记录到磁带上；

（2）从磁带上读取地震信号到计算机，对叠前地震资料做常规的噪音衰减和叠前时间偏移处理，得到三维偏移叠加数据体；

（3）在偏移叠加数据体的剖面图像上，在不同水平位置，考虑同相轴的清晰程度和连续性，确定一个最浅的清晰和连续的同相轴；剖面图像的水平坐标为平行于测线和垂直于测线的两个距离，深度坐标为旅行时，二维显示时取垂直于测线距离为常量得到平行测线方向切片的剖面图，取平行于测线距离为常量得到垂直测线方向切片的剖面图；根据给定的平行和垂直间距，在剖面图像的水平坐标面上确定一组水平样点，水平样点的水平坐标应在选定同相轴的高信噪比部位，且两个相邻样点的距离应在0.5到1.5倍的给定间距之间；在每个水平样点，从最浅的选定同相轴开始，通过选取清晰和连续的同相轴，确定一组旅行时由浅到深变化的采样点，定义最浅的采样点为浅部参考点，其它采样点为Q值拾取点；

（4）确定等效Q值的取值范围，按1/Q等间距选取系列Q值；定义三维偏移叠加数据体中相同水平坐标下的一组随旅行时变化的数据为一个成像道，用每个Q值，对每个选定的水平样点的成像道应用频散校正反Q滤波，在每个水平样点处得到一组随旅行时和Q值变化的二维道集；

（5）利用步骤4得到的二维道集和基于测井资料的合成地震记录，使用谱比的频率导数指标，得到全部浅部参考点处的等效Q值；

（6）用步骤4得到的二维道集和步骤5得到的浅部参考点的等效Q值，确定全部Q值拾取点处的等效Q值；

（7）根据步骤5和步骤6得到的所有浅部参考点和Q值拾取点处的等效Q值，用基于1/Q的光滑插值，得到三维偏移叠加数据体所包含的全部采样点上的等效Q值；

（8）用步骤7得到的非均匀等效Q值场，对三维偏移叠加数据体的每个成像道应用稳定的非均匀反Q滤波，得到更新的三维成像数据体；

（9）通过显示软件将三维成像数据体转换为地下反射构造的剖面图像，剖面图像将指示地下构造的微小断裂和精细的地层沉积样式，用于确定地下生、储油构造和识别油气储层。

所述的确定等效Q值的取值范围，按1/Q等间距选取系列Q值；定义三维偏移叠加数据体中相同水平坐标下的一组随旅行时变化的数据为一个成像道，用每个Q值，对每个选定的水平样点的成像道应用频散校正反Q滤波，在每个水平样点处得到一组随旅行时和Q值变化的二维道集是这样实现的：设等效Q值的最小、最大可能值为Q_min、Q_max，则选取的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

Q_i是无量纲参数，i和N是正整数；

用无量纲数组I(x,y,kΔτ)表示三维偏移叠加数据体，其中，(x,y)是水平坐标，单位米，k是正整数，代表第k个旅行时采样点，Δτ是偏移叠加数据体的双程旅行时时间采样，单位秒，则f_n(kΔτ)=I(x_n,y_n,kΔτ)代表一个成像道，其中(x_n,y_n)是第n个水平样点的水平坐标；设偏移叠加数据体的旅行时采样点数为M，而满足快速傅立叶算法的时间采样点数应为M₀，当M₀>M时在成像道的尾部添加零值，则傅立叶变换后的角频率采样为Δω=2π/(M₀Δτ)，单位弧度/秒；令偏移叠加数据体的有效频带下限是ω_min，而期望通过粘性吸收补偿达到的频带上限是ω_max，其单位是弧度/秒；通过求整得到正整数l₁和l₂，使得l₁Δω与ω_min和l₂Δω与ω_max的差的绝对值最小；对每个f_n(kΔτ)应用快速傅立叶变换算法，得到离散的频率域成像道为F_n(mΔω)，m是离散的频率样点序号；对每一成像道F_n(mΔω)，将存储(l₂+15)-l₁+1个非零复数值；用全部水平样点的平均F_n(mΔω)的随频率变化的模，来确定偏移叠加数据体的主频ω_p；

计算参数 $g=\frac{\mathrm{\π}}{M_0{Q}_i},g_0=\frac{g{\mathrm{\ω}}_p}{\mathrm{\Δ\ω}},a=\frac{2}{M_0{Q}_i}$ 和 $b=\frac{2\mathrm{\π}}{M_0}-\frac{2}{M_0{Q}_i}\ln \left(\frac{\mathrm{\Δ\ω}}{{\mathrm{\ω}}_p}\right),$ 其单位为弧度；再计算正整数

若m₀>15，则令m₀=15；对选定的水平样点循环，在每个水平样点(x_n,y_n)，对频率域成像道的高频端应用平滑处理，即

F_n(mΔω)=F_n(mΔω)exp[-0.06(m-l₂)²]，m=l₂+1,l₂+15

对频率循环，m=l₁,l₂+m₀，计算弧度c=m(b-aln(m))，进而计算实数d₁=exp(m·g)cos(c)和d₂=exp(m·g)sin(c)；在每一频率，对旅行时采样点循环，k=1,M，递推计算复数

P(kΔτ,mΔω)=P[(k-1)Δτ,mΔω](d₁+jd₂)

式中j是单位虚数，当k=1时P[(k-1)Δτ,mΔω]=F_n(mΔω)；

对每个k，计算正整数

若m_k>15，则令m_k=15；再对P(kΔτ,mΔω)的实部进行累加，得到新的时间序列

${\stackrel{\‾}{f}}_n(\mathrm{k\Δ\τ},Q_i)=\underset{m=l_1}{\overset{l_2+m_k}{\mathrm{\Σ}}}\mathrm{Re}\left\{P(\mathrm{k\Δ\τ},\mathrm{m\Δ\ω})\right\}$

式中Re代表取实部；完成对全部N个Q_i和选定的水平样点循环，就在每个水平样点处得到一组随旅行时和Q变化的二维道集

所述的利用步骤4得到的二维道集和基于测井资料的合成地震记录，使用谱比的频率导数指标，得到全部浅部参考点处的等效Q值是这样实现的：（1）根据已有的井资料对选定的浅部参考点分组，按与井点位置的直线距离分为与已有井资料的个数相同的组；对每个井和每组浅部参考点，选取该组中与井直线距离最近的浅部参考点为基准点，令ω_r=0.6(ω_p+ω_m)，其中ω_p和ω_m分别是偏移叠加数据体的主频和有效频带上限，单位是弧度/秒；用以ω_r为主频的雷克子波，利用测井资料得到以基准点的旅行时为中点的旅行时时窗内的合成地震记录，用基准点的对应时窗内的随Q变化的波形与合成地震记录对比，按照波形相似程度和合理性，确定基准点处的等效Q值；（2）在每个选定的水平样点，对浅部参考点处旅行时时窗内的随Q变化的波形做短时傅立叶变换，计算傅立叶变换的模，在有效频带内，求模的自然对数；模的自然对数由如下方法计算：令有效频带内，随频率变化的模D(ω)的平均值是实数D₀，当D(ω)≥0.22D₀，直接计算其对数ln[D(ω)]；当D(ω)<0.22D₀时，计算实数d=5D(ω)D₀，进而求得自然对数为

$\ln \left[\mathrm{D\&omega;}\right]=\ln \left(0.2D_0\right)+(d-1)(1-\frac{1}{2}(d-1)+\frac{1}{3}{(d-1)}^2)$

（3）对计算得到的每组对数做中值滤波处理；（4）对每组浅部参考点，以基准点处求得的等效Q值对应的一组对数值作为参考值，将其它浅部参考点的随Q变化的每组对数值分别与该参考值相减；对每个Q值对应的随频率变化的相减结果，在有效频带内确定计算频率导数的有效区间，要求区间内突变的峰值是正、负成对出现的；对相减结果做高斯平滑，然后用中心差分法计算其对频率的一阶导数，对得到的导数值做中值滤波，然后求导数的算术平均值，平均值的绝对值就是确定等效Q值的单一的频率导数指标；（5）以水平坐标为水平面，等效Q值为纵坐标，在每个水平样点处，在指定的横向窗口内，将频率导数指标做成随Q值变化的曲线，根据频率导数指标接近零和等效Q值的横向连续性，得到各个浅部参考点处的等效Q值。

所述的用步骤4得到的二维道集和步骤5得到的浅部参考点的等效Q值，确定全部Q值拾取点处的等效Q值是这样实现的：（1）在每个水平样点，根据步骤5得到的浅部参考点的等效Q值，从随旅行时和Q值变化的二维道集中拾取该等效Q值对应的浅部参考点处的短时波形，并按照步骤5中描述的方法求得中值滤波处理后的一组对数值，定义该组对数值作为参考值；（2）对各个Q值拾取点处随Q变化的波形，按照步骤5中描述的方法得到中值滤波处理后的对数值组；（3）把各个Q值拾取点的随Q变化的对数值组分别与浅部参考点的参考值相减；（4）对每个随频率变化的相减结果，在有效频带内确定计算频率导数的有效区间，要求区间内突变的峰值是正、负成对出现的；对相减结果做高斯平滑，然后用中心差分法计算其对频率的一阶导数，对得到的导数值做中值滤波，然后求导数的算术平均值，平均值的绝对值就是对应等效Q值的频率导数指标；（5）形成水平样点处不同时间深度上随Q值变化的频率导数指标图；该图做法如下：在各个Q值拾取点的旅行时处，在指定的纵向窗口内，将频率导数指标做成随Q值变化的曲线；（6）基于这一图形，在频率导数指标接近零的邻域，考虑等效Q值的横向连续性并保证等效Q值应随旅行时增加而增加，直接拾取等效Q值；（7）对每个水平样点采用相同步骤，得到全部Q值拾取点处的等效Q值。

所述的用步骤7得到的非均匀等效Q值场，对三维偏移叠加数据体的每个成像道应用稳定的非均匀反Q滤波，得到更新的三维成像数据体是这样实现的：根据地震资料的信噪比，确定幅值补偿的最大倍数A，可取为1000.0，该参数无量纲；设偏移叠加数据体的旅行时采样点数为M，而满足快速傅立叶算法的时间采样点数应为M₀；定义无量纲参数

式中各参数与步骤4中所述的相同；计算参数

其单位为弧度；计算正整数

和

其中Q₀是非均匀等效Q值场的最小值；若n₁≥n₃，则n₁=n₃和n₂=n₃，若n₂≥n₃，则n₂=n₃；设Δτ是偏移叠加数据体的双程旅行时时间采样，单位秒，则傅立叶变换后的角频率采样为Δω=2π/(M₀Δτ)，单位弧度/秒；令偏移叠加数据体的有效频带下限是ω_min，而期望通过粘性吸收补偿达到的频带上限是ω_max，其单位是弧度/秒；通过求整得到正整数l₁和l₂，使得l₁Δω与ω_min和l₂Δω与ω_max的差的绝对值最小；设ω_p是偏移叠加数据体的主频，单位是弧度/秒；

建立无量纲数组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₂

a_i=1.1A,n₂<i≤n₃

对三维偏移叠加数据体的全部成像道循环，对每个成像道应用快速傅立叶变换算法，得到离散的频率域成像道为F_n(mΔω)，其中整数n是成像道的序号，整数m=l₁,l₂+15是离散的频率样点序号；对频率域成像道的高频端应用平滑处理，即

F_n(mΔω)=F_n(mΔω)exp[-0.06(m-l₂)²]，m=l₂+1,l₂+15

计算和

其单位为弧度；计算无量纲参数 $c_1=\ln \left(\frac{\mathrm{\&Delta;\&omega;}}{{\mathrm{\&omega;}}_p}\right);$

对旅行时采样点循环，k=1,M，从非均匀等效Q值场中拾取该采样点处的等效Q值，用q₀对等效Q值的到数取整，得整数n_q；在每一旅行时采样点，对频率循环，m=l₁,l₂，计算弧度c=m·k(e-c₂(c₁+ln(m))和正整数n₀=m·k；在数组a_i中拾取第n₀个元素值a(n₀)，计算无量纲参数d₁=a(n₀)cos(c)和d₂=a(n₀)sin(c)，再计算复数

P(kΔτ,mΔω)=F_n(mΔω)(d₁+jd₂)

在上述频率循环中记录主频ω_p对应的无量纲补偿因子a(ω_p)和最大频率对应的补偿因子a_max，其中

据此计算正整数 $m_k=\mathrm{int}\left(\sqrt{\frac{1}{0.06}\ln (a_{\max }/a\left({\mathrm{\&omega;}}_p\right))+25}\right);$ 若m_k>15，则令m_k=15；

继续对频率循环，m=l₂+1,l₂+m_k，计算弧度c=m·k(e-c₂(c₁+ln(m))，计算无量纲参数d₁=a_maxcos(c)和d₂=a_maxsin(c)，再计算复数

P(kΔτ,mΔω)=F_n(mΔω)(d₁+jd₂)

对P(kΔτ,mΔω)的实部进行累加，得到新的时间序列

${\stackrel{\&OverBar;}{f}}_n\left(\mathrm{k\&Delta;\&tau;}\right)=\underset{m=l_1}{\overset{l_2+m_k}{\mathrm{\&Sigma;}}}\mathrm{Re}\left\{P(\mathrm{k\&Delta;\&tau;},\mathrm{m\&Delta;\&omega;})\right\}$

式中Re代表取实部，

即是粘性吸收补偿后的成像道；完成对全部成像道的循环，就得到更新的三维成像数据体。

本发明的提高地震成像分辨率的方法，能整体提升成像的分辨率，并使得中、深层构造的成像分辨率达到与浅层接近的程度。

本发明的提高地震成像分辨率的方法，通过对偏移叠加数据体进行频率相关的幅值恢复和频散校正，并保证幅值恢复的稳定性，来提高地震成像的分辨率。

本发明的提高地震成像分辨率的方法，通过引入等效Q值描述地震波的粘性吸收，利用参数扫描得到计算所需的非均匀等效Q值场。

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

本发明的核心有三点，一是从粘弹性波传播原理出发，给出等效Q值的定义和基于等效Q值的地震波幅值衰减与频散的解析表达式；二是基于一维深度偏移方法，给出基于非均匀Q值的稳定的频散校正和噪音压制反Q滤波方法；三是发展了确定非均匀等效Q值的方法。其具体实现原理如下：

1．等效Q值和地震波幅值衰减与频散

地震波在地球介质中的衰减可用基于Q值的粘弹性波理论描述，一般假设Q值不随地震波的频率变化，即使用衡Q理论。粘弹性波理论通过引入随频率变化的复相速度和实相速度，来描述地震波的幅值衰减与频散；衡Q理论分别定义复相速度和实相速度为

$\frac{1}{c\left(\mathrm{\&omega;}\right)}=\frac{1}{v_r\left(\mathrm{\&omega;}\right)}(1-\frac{j}{2Q})---\left(1a\right)$

$v_r\left(\mathrm{\&omega;}\right)=v_r\left({\mathrm{\&omega;}}_c\right)[1+\frac{1}{\mathrm{\&pi;Q}}\ln \left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_c}\right)---\left(1b\right)$

式中，j是单位虚数，ω是角频率，ω_c是实相速度变化的截止角频率，超过这个频率，实相速度就不随频率变化；c(ω)是介质的复相速度，v_r(ω)是介质的实相速度；Q就是品质因子，即Q值，它是一个无量纲参数，反应了地震波传播过程中原有能量与被吸收能量的比值。

实际地球介质满足1/Q<<1，若令ω_p是地震资料的主频，则由式（1b）可得：

$v_r\left(\mathrm{\&omega;}\right)=v_r\left({\mathrm{\&omega;}}_p\right)\frac{1+[1/\left(\mathrm{\&pi;Q}\right)]\ln (\mathrm{\&omega;}/{\mathrm{\&omega;}}_c)}{1+[1/(\mathrm{\&pi;Q}\left)\right]\ln ({\mathrm{\&omega;}}_p/{\mathrm{\&omega;}}_c)}\&ap;\frac{v_r\left({\mathrm{\&omega;}}_p\right)}{1-[1/\left(\mathrm{\&pi;Q}\right)]\ln (\mathrm{\&omega;}/{\mathrm{\&omega;}}_p)}---\left(2\right)$

将式（2）带入式（1a）得

$\frac{1}{c\left(\mathrm{\&omega;}\right)}=\frac{1}{v}(1-\frac{j}{2Q})[1-\frac{1}{\mathrm{\&pi;Q}}\ln \left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_0}\right)---\left(3\right)$

式中v是主频ω_p对应的实相速度。由于v可通过偏移速度建模方法得到，式（3）代表了可由地震资料求得的复相速度。

偏移叠加数据体可看作垂直入射和出射的反射地震信号，因此可用一维介质理论处理叠加数据体。若将非均匀介质分为若干均匀介质层，在频率域，基于粘弹性波理论，入射波场沿深度的传播可表示为：

$\begin{array}{c}\tilde{P}(\mathrm{\&omega;},z=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{z}_i)=S\left(\mathrm{\&omega;}\right)\exp (-j\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{z}_i\frac{\mathrm{\&omega;}}{c_i\left(\mathrm{\&omega;}\right)})\\ \&ap;S\left(\mathrm{\&omega;}\right)\exp \{-\mathrm{j\&omega;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&Delta;}{z}_i}{v_i}(1-\frac{1}{\mathrm{\&pi;}{Q}_i}\ln \left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_p}\right)-\frac{j}{2Q_i})\}\end{array}---\left(4\right)$

式中

是深度z处的频率域波场，Δz_i是各层介质的厚度，n是目的层以上的层状介质的层数，c_i(ω)是各层介质的复相速度，v_i为各层介质的实相速度，Q_i为各层介质的品质因子，S(ω)是震源子波的傅立叶变换。

若用垂直旅行时T表达深度，有Δτ_i=Δz_i/c_i，

则式（4）可表达为：

$\tilde{P}(\mathrm{\&omega;},T)=S\left(\mathrm{\&omega;}\right)\exp \{-\mathrm{j\&omega;T}(1-\frac{1}{\mathrm{\&pi;}{Q}_{\mathrm{eff}}}\ln \left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_p}\right))-\frac{\mathrm{\&omega;T}}{2Q_{\mathrm{eff}}}\}---\left(5\right)$

式中Q_eff是本发明引入的等效Q值，其表达式为：

$\frac{1}{Q_{\mathrm{eff}}}=\frac{1}{T}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&Delta;}{\mathrm{\&tau;}}_i}{Q_i}---\left(6\right)$

引入等效Q值并用垂直旅行时表达深度，使得式（5）的粘弹性波场传播表达式仅与单个参数Q_eff有关，这使得可用扫描方法获得等效Q值。式（5）中的

即是频散项，而

是幅值衰减项。

2．基于非均匀Q值的反Q滤波

可基于一维深度偏移方法处理偏移叠加数据体中的每个成像道。式（5）给出了垂直入射的下行波场，同理，可得反传的记录波场为(即深度延拓）

$\tilde{U}(\mathrm{\&omega;},T)=\tilde{f}\left(\mathrm{\&omega;}\right)\exp \{\mathrm{j\&omega;T}(1-\frac{1}{\mathrm{\&pi;}{Q}_{\mathrm{eff}}}\ln \left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_p}\right))+\frac{\mathrm{\&omega;T}}{2Q_{\mathrm{eff}}}\}---\left(7\right)$

式中

是单个成像道的时域信号的傅立叶变换；定义三维偏移叠加数据体中相同水平坐标下的一组随旅行时变化的数据为一个成像道。

一般情况下震源子波是未知的。既然反褶积可剔除子波的影响，可在成像中忽略震源子波项，即式（5）中的S(ω)。对式（5）和式（7）应用反褶积成像条件，可得粘性吸收补偿后的成像道为

$I\left(\mathrm{\&tau;}\right)=\underset{\mathrm{\&omega;}}{\mathrm{\&Sigma;}}\mathrm{Re}\left\{\tilde{f}\left(\mathrm{\&omega;}\right)\exp \left(\mathrm{j\&omega;\&tau;}\right)\exp (-\frac{\mathrm{j\&omega;\&tau;}}{\mathrm{\&pi;}{Q}_{\mathrm{eff}}}\ln \left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_p}\right))\exp \left(\frac{\mathrm{\&omega;\&tau;}}{2Q_{\mathrm{eff}}}\right)\right\}---\left(8\right)$

式中用双倍旅行时替代T，即τ=2T；Re代表取实部。式（8）与常规反Q滤波有些相识，不同之处是：1）Q_eff是随τ变化的，2）增加了一个频散校正项；此外，在式（8）的数值实现中，我们引入了针对幅值补偿项的稳定处理和压制由频率域折返效应导致的高频噪音的方法。

下面首先讨论保证幅值补偿稳定的方法。当τ较大或Q_eff较小时，式（8）中的补偿因子，即

将趋于无穷大，为此提出一个光滑性阈值控制的稳定性算法。具体实现如下：

令幅值补偿的最大倍数为A，设计一个光滑函数，使其在幅值小于等于A时，即当

时，与补偿因子完全一致；当

时，其幅值光滑过渡到1.1A，然后保持恒定。为避免重复计算这一函数，可事先计算相应的结果并将其存储在一个一维表中，每次按

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

其中q为拾取等效Q值时使用的1/Q间距的十分之一，M是成像道的时间采样点数，可按下式计算一维表中的元素：

a_i=exp(iΔ),i≤n₁   （9）

a_i=A(1-lnA-2.5(lnA)²)+A(1+5lnA)iΔ-2.5A(iΔ)²,n₁<i≤n₂   （10）

上两式中式（10）的函数保证了函数数值和其一阶导数连续。根据最高频率和最小等效Q值，可得到一维表的上限n₃；因此，在构建一维表时，若n₁≥n₃，则n₁=n₃和n₂=n₃，若n₂≥n₃，则n₂=n₃。

实际计算时，可对

取整；若该整数大于n₂，取λ=1.1A；否则，在一维表中拾取相邻两点的值，插值得到λ，用λ替代式（8）中的

进行计算。这样即避免了指数运算，又可保证补偿因子小于或等于1.1A，保证了幅值补偿的无条件稳定。

式（8）的计算采用离散的频率值，频率采样是由快速傅立叶算法决定的。设满足快速傅立叶算法的时间采样点数为M₀，令Δτ是成像道双程旅行时的时间采样，则频率采样为Δω=2π/(M₀Δτ)。

令偏移叠加数据体的有效频带的下限是ω_min，而期望通过粘性吸收补偿达到的频带上限是ω_max，可通过求整得到整数l₁和l₂，使得l₁Δω和l₂Δω与ω_min和ω_max最接近，则式（8）的离散表达式为

$I\left(\mathrm{k\&Delta;\&tau;}\right)=\underset{m=l_1}{\overset{l_2}{\mathrm{\&Sigma;}}}\mathrm{Re}\left\{\tilde{f}\left(\mathrm{m\&Delta;\&omega;}\right)\exp (\mathrm{jkm\&Delta;\&omega;}\&CenterDot;\mathrm{\&Delta;\&tau;})\exp (-\frac{\mathrm{j\&omega;m\&Delta;\&omega;}\&CenterDot;\mathrm{\&Delta;\&tau;}}{\mathrm{\&pi;}{Q}_{\mathrm{eff}}}\ln \left(\frac{\mathrm{m\&Delta;\&omega;}}{{\mathrm{\&omega;}}_p}\right))\exp \left(\frac{\mathrm{km\&Delta;\&omega;}\&CenterDot;\mathrm{\&Delta;\&tau;}}{2Q_{\mathrm{eff}}}\right)\right\}---\left(11\right)$

若不引入补偿因子，式（11）的右端项的模在高频截止点m=l₂处已平滑地接近零，其对应的傅立叶反变换结果，即I(kΔt)，将不产生频率域的折返效应。由于引入补偿因子，式（11）的右端项的模在高频截止点将远大于零，此时继续使用式（11）将产生由折返效应导致的高频噪音。为此，本发明在高频截止点增加光滑的过渡带，籍此衰减高频噪音。具体实现如下：

对成像道超出高频截止点的频率成份，引入衰减带如下：

$\tilde{f}=\left(\mathrm{m\&Delta;\&omega;}\right)=\tilde{f}\left(\mathrm{m\&Delta;\&omega;}\right)\exp [-0.06{(i-l_2)}^2],i=l_2+1,l_2+15---\left(12\right)$

在式（11）的频率域累加计算中，高频端衰减带的实际宽度是根据补偿因子的相对大小而变化的；令地震资料主频对应的补偿因子为a(ω_p)，最大频率点对应的补偿因子为a_max，据此计算整数m_k，

$m_k=\mathrm{int}\left(\sqrt{\frac{1}{0.06}\ln (a_{\max }/a\left({\mathrm{\&omega;}}_p\right))+25}\right);$ 若m_k>15，令m_k=15   （13）进而计算

$I_0\left(\mathrm{k\&Delta;\&tau;}\right)=\underset{m=l_2+1}{\overset{l_2+m_k}{\mathrm{\&Sigma;}}}\mathrm{Re}\left\{\tilde{f}\left(\mathrm{m\&Delta;\&omega;}\right)\exp (\mathrm{jkm\&Delta;\&omega;}\&CenterDot;\mathrm{\&Delta;\&tau;})\exp (-\frac{\mathrm{jkm\&Delta;\&omega;}\&CenterDot;\mathrm{\&Delta;\&tau;}}{\mathrm{\&pi;}{Q}_{\mathrm{eff}}}\ln \left(\frac{\mathrm{m\&Delta;\&omega;}}{{\mathrm{\&omega;}}_p}\right))\exp \left(\frac{k{l}_2\mathrm{\&Delta;\&omega;}\&CenterDot;\mathrm{\&Delta;\&tau;}}{2Q_{\mathrm{eff}}}\right)\right\}---\left(14\right)$

以I(kΔt)+I₀(kΔt)为更新的成像道，将有效衰减折返效应导致的高频噪音。

当引入阈值控制的稳定性算法时，式（11）和（14）中的稳定因子是由一维表中拾取得到的。

当Q_eff是常量时，式（11）和（14）的计算应先固定频率，然后对时间，即对k循环；这样对每一频率，可事先计算复数

$g_0=\exp (\mathrm{jm\&Delta;\&omega;}\&CenterDot;\mathrm{\&Delta;\&tau;})\exp (-\frac{\mathrm{jm\&Delta;\&omega;}\&CenterDot;\mathrm{\&Delta;\&tau;}}{\mathrm{\&pi;}{Q}_{\mathrm{eff}}}\ln \left(\frac{\mathrm{m\&Delta;\&omega;}}{{\mathrm{\&omega;}}_p}\right))\exp \left(\frac{\mathrm{m\&Delta;\&omega;}\&CenterDot;\mathrm{\&Delta;\&tau;}}{2Q_{\mathrm{eff}}}\right)$

然后重复利用该复数进行累加计算，即，式（11）和（14）中随k变化的复数可有下式递推得到

g_k+1=g_k·g₀，

这样可极大地减少式（11）和（14）的计算量。

3．等效Q值拾取

估计Q值是一个非常困难的问题，特别是在Q值非均匀的情况下。为此我们提出了等效Q值的概念，由于在每一样点的吸收补偿是由该点的等效Q值唯一决定的，使得可利用扫描方法直接求取等效Q值。

在利用透射信号，如VSP资料，估计Q值时，可利用频谱形状改变等信息决定Q值，因此发展了谱比法等估计算法。由于薄层调谐等现象存在，将谱比法等方法应用于反射地震资料则产生了新的困难；主要原因是此时主频的移动或频谱形状的变化受到粘性吸收和薄层调谐的共同影响，而很多时候薄层调谐的影响更大。为此，我们提出了用谱比的频率导数指标确定等效Q值的方法。

谱比的频率导数指标由如下方法计算：设粘性补偿后样点的时域波形的短时傅立叶变换的模为G(ω)，而参考点处，采用正确Q值补偿后，应用短时傅立叶变换得到的模为G₀(ω)，则谱比的频率导数指标为

$\mathrm{\&eta;}=\underset{\mathrm{\&omega;}}{\mathrm{\&Sigma;}}\frac{d}{\mathrm{d\&omega;}}(\ln G\left(\mathrm{\&omega;}\right)-\ln G_0\left(\mathrm{\&omega;}\right))---\left(15\right)$

与谱比法相比，式（15）的频率导数指标可剔除薄层调谐对谱比计算结果的影响；这是因为薄层调谐项在求频率导数后变为

算术平均即可剔除该项。为避免求导计算对噪音的放大，需先对随频率变化的ln(G(ω))等做中值滤波和平滑等处理；在求算术平均值前，也需对导数值做中值滤波处理。为避免薄层调谐产生的零值所导致的对数计算的不稳定，式(15)中对数的计算需采用稳定的对数算法。

若指标η的绝对值为零或最小，即表明此时的G(ω)是得到了正确补偿，由此可确定该样点正确的等效Q值。因此，η可作为从扫描计算中拾取准确等效Q值的指标。

本发明的有益效果：该方法能整体提升成像的分辨率，并使得中、深层构造的成像分辨率达到与浅层接近的程度，从而获得更高分辨率的地下构造图像，能更好地识别小断层和薄储层。该方法的应用不需要事先提供非均匀Q值场和其它参数。该方法所需计算量较少，适用于大规模的三维偏移叠加数据体。该方法对陆相薄互层油气储层的勘探和开发有重要应用价值。

附图说明

图1a是东部油田某区块典型的三维偏移叠加剖面图像，是垂直测线方向坐标2km处的平行测线方向切片的剖面图；图1b是东部油田某区块典型的三维偏移叠加剖面图像，是平行测线方向坐标2km处的垂直测线方向切片的剖面图。

图2是图1a剖面中4－8km部分上的浅部参考点和Q值拾取点，图中叉号标记代表浅部参考点和Q值拾取点。

图3a是得到的等效Q值场的等值线图，是垂直测线方向坐标2km处的平行测线方向切片的等值线图；图3b是得到的等效Q值场的等值线图，是平行测线方向坐标2km处的垂直测线方向切片的等值线图；图3a和图3b中等值线的数字代表1/Q值大小。

图4a是应用本发明处理后的三维偏移叠加剖面图像，是垂直测线方向坐标2km处的平行测线方向切片的剖面图；图4b是应用本发明处理后的三维偏移叠加剖面图像，是平行测线方向坐标2km处的垂直测线方向切片的剖面图。图4a和图4b的成像剖面切片的位置与图1a和图1b相同。

图5是应用本发明前后偏移叠加剖面图像中部（时窗1.2－1.5s）的频谱比较图。图中点线是应用本发明前的频谱，实线是应用本发明后的频谱；从图中可见，主频明显提高，最高频率提高了10Hz。

具体实施方式

实施例1：提高地震成像分辨率的方法，以东部陆上油田为例，具体为以下步骤：

（1）用多条拖缆或测线采集人工震源激发的反射地震信号，将地震信号记录到磁带上。具体采集参数是，12条测线同时记录地震信号，测线间距80m，每条测线有2801个检波器组，检波器组间距40m；沿垂直测线方向布炮线，炮点沿炮线的间距40m，炮线间的线间距40m，共采集5679炮，记录长度2.5s，时间采样1ms。

（2）从磁带上读取地震信号到计算机，对叠前地震资料做常规的噪音衰减和叠前时间偏移处理，得到三维偏移叠加数据体。三维偏移叠加数据体共含149160个成像道，在平行测线和垂直测线方向的道间距分别是20m和40m，覆盖面积120平方千米；用旅行时表达的深度为2.5s，旅行时采样1ms。图1a是偏移叠加数据体平行测线方向切片的剖面图，图1b是偏移叠加数据体垂直测线方向切片的剖面图。

（3）在偏移叠加数据体的剖面图像上，在不同水平位置，考虑同相轴的清晰程度和连续性，确定一个最浅的清晰和连续的同相轴；剖面图像的水平坐标为平行于测线和垂直于测线的两个距离，深度坐标为旅行时，二维显示时取垂直于测线距离为常量得到平行测线方向切片的剖面图，取平行于测线距离为常量得到垂直测线方向切片的剖面图；根据给定的平行和垂直间距，在剖面图像的水平坐标面上确定一组水平样点，水平样点的水平坐标应在选定同相轴的高信噪比部位，且两个相邻样点的距离应在0.5到1.5倍的给定间距之间；在每个水平样点，从最浅的选定同相轴开始，通过选取清晰和连续的同相轴，确定一组旅行时由浅到深变化的采样点，定义最浅的采样点为浅部参考点，其它采样点为Q值拾取点。图2给出了垂直测线方向坐标2km处的平行测线方向切片的局部剖面上的浅部参考点和Q值拾取点。

（4）确定等效Q值的取值范围，按1/Q等间距选取系列Q值；定义三维偏移叠加数据体中相同水平坐标下的一组随旅行时变化的数据为一个成像道，用每个Q值，对每个选定的水平样点的成像道应用频散校正反Q滤波，在每个水平样点处得到一组随旅行时和Q值变化的二维道集。具体是，令等效Q值的取值范围为50到300，共选取12个1/Q等间距的系列Q值。

设等效Q值的最小、最大可能值为Q_min、Q_max，则选取的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$

Q_i是无量纲参数，i和N是正整数；

用无量纲数组I(x,y,kΔτ)表示三维偏移叠加数据体，其中，(x,y)是水平坐标，单位米，k是正整数，代表第k个旅行时采样点，Δτ是偏移叠加数据体的双程旅行时时间采样，单位秒，则f_n(kΔτ)=I(x_n,y_n,kΔτ)代表一个成像道，其中(x_n,y_n)是第n个水平样点的水平坐标；设偏移叠加数据体的旅行时采样点数为M，而满足快速傅立叶算法的时间采样点数应为M₀，当M₀>M时在成像道的尾部添加零值，则傅立叶变换后的角频率采样为Δω=2π/(M₀Δτ)，单位弧度/秒；令偏移叠加数据体的有效频带下限是ω_min，而期望通过粘性吸收补偿达到的频带上限是ω_max，其单位是弧度/秒；通过求整得到正整数l₁和l₂，使得l₁Δω与ω_min和l₂Δω与ω_max的差的绝对值最小；对每个f_n(kΔτ)应用快速傅立叶变换算法，得到离散的频率域成像道为F_n(mΔω)，m是离散的频率样点序号；对每一成像道F_n(mΔω)，将存储(l₂+15)-l₁+1个非零复数值；用全部水平样点的平均F_n(mΔω)的随频率变化的模，来确定偏移叠加数据体的主频ω_p；

计算参数 $g=\frac{\mathrm{\&pi;}}{M_0{Q}_i},g_0=\frac{g{\mathrm{\&omega;}}_p}{\mathrm{\&Delta;\&omega;}},a=\frac{2}{M_0{Q}_i}$ 和 $b=\frac{2\mathrm{\&pi;}}{M_0}-\frac{2}{M_0{Q}_i}\ln \left(\frac{\mathrm{\&Delta;\&omega;}}{{\mathrm{\&omega;}}_p}\right),$ 其单位为弧度；再计算正整数

若m₀>15，则令m₀=15；对选定的水平样点循环，在每个水平样点(x_n,y_n)，对频率域成像道的高频端应用平滑处理，即

F_n(mΔω)=F_n(mΔω)exp[-0.06(m-l₂)²]，m=l₂+1,l₂+15

对频率循环，m=l₁,l₂+m₀，计算弧度c=m(b-aln(m))，进而计算实数d₁=exp(m·g)cos(c)和d₂=exp(m·g)sin(c)；在每一频率，对旅行时采样点循环，k=1,M，递推计算复数

P(kΔτ,mΔω)=P[(k-1)Δτ,mΔω](d₁+jd₂)

式中j是单位虚数，当k=1时P[(k-1)Δτ,mΔω]=F_n(mΔω)；

对每个k，计算正整数

若m_k>15，则令m_k=15；再对P(kΔτ,mΔω)的实部进行累加，得到新的时间序列

${\stackrel{\&OverBar;}{f}}_n(\mathrm{k\&Delta;\&tau;},Q_i)=\underset{m=l_1}{\overset{l_2+m_k}{\mathrm{\&Sigma;}}}\mathrm{Re}\left\{P(\mathrm{k\&Delta;\&tau;},\mathrm{m\&Delta;\&omega;})\right\}$

式中Re代表取实部；完成对全部N个Q_i和选定的水平样点循环，就在每个水平样点处得到一组随旅行时和Q变化的二维道集

（5）利用步骤4得到的二维道集和基于测井资料的合成地震记录，使用谱比的频率导数指标，得到全部浅部参考点处的等效Q值。

采用以下5个步骤：1）根据已有的井资料对选定的浅部参考点分组，按与井点位置的直线距离分为与已有井资料的个数相同的组；对每个井和每组浅部参考点，选取该组中与井直线距离最近的浅部参考点为基准点，令ω_r=0.6(ω_p+ω_m)，其中ω_p和ω_m分别是偏移叠加数据体的主频和有效频带上限，单位是弧度/秒；用以ω_r为主频的雷克子波，利用测井资料得到以基准点的旅行时为中点的旅行时时窗内的合成地震记录，用基准点的对应时窗内的随Q变化的波形与合成地震记录对比，按照波形相似程度和合理性，确定基准点处的等效Q值；2）在每个选定的水平样点，对浅部参考点处旅行时时窗内的随Q变化的波形做短时傅立叶变换，计算傅立叶变换的模，在有效频带内，求模的自然对数；模的自然对数由如下方法计算：令有效频带内，随频率变化的模D(ω)的平均值是实数D₀，当D(ω)≥0.22D₀，直接计算其对数ln[D(ω)]；当D(ω)<0.22D₀时，计算实数d=5D(ω)/D₀，进而求得自然对数为

$\ln \left[\mathrm{D\&omega;}\right]=\ln \left(0.2D_0\right)+(d-1)(1-\frac{1}{2}(d-1)+\frac{1}{3}{(d-1)}^2)$

3）对计算得到的每组对数做中值滤波处理；4）对每组浅部参考点，以基准点处求得的等效Q值对应的一组对数值作为参考值，将其它浅部参考点的随Q变化的每组对数值分别与该参考值相减；对每个Q值对应的随频率变化的相减结果，在有效频带内确定计算频率导数的有效区间，要求区间内突变的峰值是正、负成对出现的；对相减结果做高斯平滑，然后用中心差分法计算其对频率的一阶导数，对得到的导数值做中值滤波，然后求导数的算术平均值，平均值的绝对值就是确定等效Q值的单一的频率导数指标；5）以水平坐标为水平面，等效Q值为纵坐标，在每个水平样点处，在指定的横向窗口内，将频率导数指标做成随Q值变化的曲线，根据频率导数指标接近零和等效Q值的横向连续性，得到各个浅部参考点处的等效Q值。

（6）用步骤4得到的二维道集和步骤5得到的浅部参考点的等效Q值，确定全部Q值拾取点处的等效Q值。

采用以下7个步骤：1）在每个水平样点，根据步骤5得到的浅部参考点的等效Q值，从随旅行时和Q值变化的二维道集中拾取该等效Q值对应的浅部参考点处的短时波形，并按照步骤5中描述的方法求得中值滤波处理后的一组对数值，定义该组对数值作为参考值；2）对各个Q值拾取点处随Q变化的波形，按照步骤5中描述的方法得到中值滤波处理后的对数值组；3）把各个Q值拾取点的随Q变化的对数值组分别与浅部参考点的参考值相减；4）对每个随频率变化的相减结果，在有效频带内确定计算频率导数的有效区间，要求区间内突变的峰值是正、负成对出现的；对相减结果做高斯平滑，然后用中心差分法计算其对频率的一阶导数，对得到的导数值做中值滤波，然后求导数的算术平均值，平均值的绝对值就是对应等效Q值的频率导数指标；5）形成水平样点处不同时间深度上随Q值变化的频率导数指标图；该图做法如下：在各个Q值拾取点的旅行时处，在指定的纵向窗口内，将频率导数指标做成随Q值变化的曲线；6）基于这一图形，在频率导数指标接近零的邻域，考虑等效Q值的横向连续性并保证等效Q值应随旅行时增加而增加，直接拾取等效Q值；7）对每个水平样点采用相同步骤，得到全部Q值拾取点处的等效Q值。

（7）根据步骤5和步骤6得到的所有浅部参考点和Q值拾取点处的等效Q值，用基于1/Q的光滑插值，得到三维偏移叠加数据体所包含的全部采样点上的等效Q值。图3a是平行测线方向切片的等效Q值的等值线图，图3b是垂直测线方向切片的等效Q值的等值线图；图3a和图3b中数字代表1/Q值的大小。

（8）用步骤7得到的非均匀等效Q值场，对三维偏移叠加数据体的每个成像道应用稳定的非均匀反Q滤波，得到更新的三维成像数据体。

根据地震资料的信噪比，确定幅值补偿的最大倍数A，可取为1000.0，该参数无量纲；设偏移叠加数据体的旅行时采样点数为M，而满足快速傅立叶算法的时间采样点数应为M₀；定义无量纲参数

式中各参数与步骤4中所述的相同；计算参数 $\mathrm{\&Delta;}=\frac{\mathrm{\&pi;}}{M_0}{q}_0,$ 其单位为弧度；计算正整数 $n_1=\frac{\ln A}{\mathrm{\&Delta;}},n_2=\frac{\ln A+0.2}{\mathrm{\&Delta;}}$ 和 $n_3=\frac{\mathrm{\&pi;}{l}_2}{\mathrm{\&Delta;}{Q}_0},$ 其中Q₀是非均匀等效Q值场的最小值；若n₁≥n₃，则n₁=n₃和n₂=n₃，若n₂≥n₃，则n₂=n₃；设Δτ是偏移叠加数据体的双程旅行时时间采样，单位秒，则傅立叶变换后的角频率采样为Δω=2π/(M₀Δτ)，单位弧度/秒；令偏移叠加数据体的有效频带下限是ω_min，而期望通过粘性吸收补偿达到的频带上限是ω_max，其单位是弧度/秒；通过求整得到正整数l₁和l₂，使得l₁Δω与ω_min和l₂Δω与ω_max的差的绝对值最小；设ω_p是偏移叠加数据体的主频，单位是弧度/秒；

建立无量纲数组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₂

a_i=1.1A,n₂<i≤n₃

对三维偏移叠加数据体的全部成像道循环，对每个成像道应用快速傅立叶变换算法，得到离散的频率域成像道为F_n(mΔω)，其中整数n是成像道的序号，整数m=l₁,l₂+15是离散的频率样点序号；对频率域成像道的高频端应用平滑处理，即

F_n(mΔω)=F_n(mΔω)exp[-0.06(m-l₂)²]，m=l₂+1,l₂+15

计算

和

其单位为弧度；计算无量纲参数 $c_1=\ln \left(\frac{\mathrm{\&Delta;\&omega;}}{{\mathrm{\&omega;}}_p}\right);$

对旅行时采样点循环，k=1,M，从非均匀等效Q值场中拾取该采样点处的等效Q值，用q₀对等效Q值的到数取整，得整数n_q；在每一旅行时采样点，对频率循环，m=l₁,l₂，计算弧度c=m·k(e-c₂(c₁+ln(m))和正整数n₀=m·k；在数组a_i中拾取第n₀个元素值a(n₀)，计算无量纲参数d₁=a(n₀)cos(c)和d₂=a(n₀)sin(c)，再计算复数

P(kΔτ,mΔω)=F_n(mΔω)(d₁+jd₂)

在上述频率循环中记录主频ω_p对应的无量纲补偿因子a(ω_p)和最大频率对应的补偿因子a_max，其中

据此计算正整数 $m_k=\mathrm{int}\left(\sqrt{\frac{1}{0.06}\ln (a_{\max }/a\left({\mathrm{\&omega;}}_p\right))+25}\right);$ 若m_k>15，则令m_k=15；

继续对频率循环，m=l₂+1,l₂+m_k，计算弧度c=m·k(e-c₂(c₁+ln(m))，计算无量纲参数d₁=a_maxcos(c)和d₂=a_maxsin(c)，再计算复数

P(kΔτ,mΔω)=F_n(mΔω)(d₁+jd₂)

对P(kΔτ,mΔω)的实部进行累加，得到新的时间序列

${\stackrel{\&OverBar;}{f}}_n\left(\mathrm{k\&Delta;\&tau;}\right)=\underset{m=l_1}{\overset{l_2+m_k}{\mathrm{\&Sigma;}}}\mathrm{Re}\left\{P(\mathrm{k\&Delta;\&tau;},\mathrm{m\&Delta;\&omega;})\right\}$

式中Re代表取实部，

即是粘性吸收补偿后的成像道；完成对全部成像道的循环，就得到更新的三维成像数据体。

（9）通过显示软件将三维成像数据体转换为地下反射构造的剖面图像，剖面图像将指示地下构造的微小断裂和精细的地层沉积样式，用于确定地下生、储油构造和识别油气储层。图4a和图4b分别是平行测线方向和垂直测线方向切片的剖面图，对应的位置与图1a和图1b相同；对比图1a和图4a、图1b和图4b可见，应用本发明后，分辨率明显提高，叠合的同相轴被更好地分离，断点更加清晰。图5是应用本发明前后剖面图像中部（时窗1.2－1.5s）的频谱比较，可见主频明显提高，最高频率提高了10Hz。

Claims (5)

1.一种提高地震成像分辨率的方法，其特征在于：包括以下步骤：步骤A、用多条拖缆或测线采集人工震源激发的反射地震信号，将地震信号记录到磁带上；步骤B、从磁带上读取地震信号到计算机，对叠前地震资料做常规的噪音衰减和叠前时间偏移处理，得到三维偏移叠加数据体；步骤C、在偏移叠加数据体的剖面图像上，在不同水平位置，考虑同相轴的清晰程度和连续性，确定一个最浅的清晰和连续的同相轴；剖面图像的水平坐标为平行于测线和垂直于测线的两个距离，深度坐标为旅行时，二维显示时取垂直于测线距离为常量得到平行测线方向切片的剖面图，取平行于测线距离为常量得到垂直测线方向切片的剖面图；根据给定的平行和垂直间距，在剖面图像的水平坐标面上确定一组水平样点，水平样点的水平坐标应在选定同相轴的高信噪比部位，且两个相邻样点的距离应在0.5到1.5倍的给定间距之间；在每个水平样点，从最浅的选定同相轴开始，通过选取清晰和连续的同相轴，确定一组旅行时由浅到深变化的采样点，定义最浅的采样点为浅部参考点，其它采样点为Q值拾取点；步骤D、确定等效Q值的取值范围，按1/Q等间距选取系列Q值；定义三维偏移叠加数据体中相同水平坐标下的一组随旅行时变化的数据为一个成像道，用每个Q值，对每个选定的水平样点的成像道应用频散校正反Q滤波，在每个水平样点处得到一组随旅行时和Q值变化的二维道集；步骤E、利用步骤D得到的二维道集和基于测井资料的合成地震记录，使用谱比的频率导数指标，得到全部浅部参考点处的等效Q值；步骤F、用步骤D得到的二维道集和步骤E得到的浅部参考点的等效Q值，确定全部Q值拾取点处的等效Q值；步骤G、根据步骤E和步骤F得到的所有浅部参考点和Q值拾取点处的等效Q值，用基于1/Q的光滑插值，得到三维偏移叠加数据体所包含的全部采样点上的等效Q值；步骤H、用步骤G得到的非均匀等效Q值场，对三维偏移叠加数据体的每个成像道应用稳定的非均匀反Q滤波，得到更新的三维成像数据体；步骤I、通过显示软件将三维成像数据体转换为地下反射构造的剖面图像，剖面图像将指示地下构造的微小断裂和精细的地层沉积样式，用于确定地下生、储油构造和识别油气储层。

2.根据权利要求1所述的一种提高地震成像分辨率的方法，其特征在于：所述的确定等效Q值的取值范围，按1/Q等间距选取系列Q值；定义三维偏移叠加数据体中相同水平坐标下的一组随旅行时变化的数据为一个成像道，用每个Q值，对每个选定的水平样点的成像道应用频散校正反Q滤波，在每个水平样点处得到一组随旅行时和Q值变化的二维道集是这样实现的：设等效Q值的最小、最大可能值为Q_min、Q_max，则选取的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

Q_i是无量纲参数，i和N是正整数；

用无量纲数组I(x,y,kΔτ)表示三维偏移叠加数据体，其中，(x,y)是水平坐标，单位米，k是正整数，代表第k个旅行时采样点，Δτ是偏移叠加数据体的双程旅行时时间采样，单位秒，则f_n(kΔτ)=I(x_n,y_n,kΔτ)代表一个成像道，其中(x_n,y_n)是第n个水平样点的水平坐标；设偏移叠加数据体的旅行时采样点数为M，而满足快速傅立叶算法的时间采样点数应为M₀，当M₀>M时在成像道的尾部添加零值，则傅立叶变换后的角频率采样为Δω=2π/(M₀Δτ)，单位弧度/秒；令偏移叠加数据体的有效频带下限是ω_min，而期望通过粘性吸收补偿达到的频带上限是ω_max，其单位是弧度/秒；通过求整得到正整数l₁和l₂，使得l₁Δω与ω_min和l₂Δω与ω_max的差的绝对值最小；对每个f_n(kΔτ)应用快速傅立叶变换算法，得到离散的频率域成像道为F_n(mΔω)，m是离散的频率样点序号；对每一成像道F_n(mΔω)，将存储(l₂+15)-l₁+1个非零复数值；用全部水平样点的平均F_n(mΔω)的随频率变化的模，来确定偏移叠加数据体的主频ω_p；

计算参数 $g=\frac{\mathrm{\&pi;}}{M_0{Q}_i},g_0=\frac{g{\mathrm{\&omega;}}_p}{\mathrm{\&Delta;\&omega;}},a=\frac{2}{M_0{Q}_i}$ 和 $b=\frac{2\mathrm{\&pi;}}{M_0}-\frac{2}{M_0{Q}_i}\ln \left(\frac{\mathrm{\&Delta;\&omega;}}{{\mathrm{\&omega;}}_p}\right),$ 其单位为弧度；再计算正整数

若m₀>15，则令m₀=15；对选定的水平样点循环，在每个水平样点(x_n,y_n)，对频率域成像道的高频端应用平滑处理，即

F_n(mΔω)=F_n(mΔω)exp[-0.06(m-l₂)²]，m=l₂+1,l₂+15

对频率循环，m=l₁,l₂+m₀，计算弧度c=m(b-aln(m))，进而计算实数d₁=exp(m·g)cos(c)和d₂=exp(m·g)sin(c)；在每一频率，对旅行时采样点循环，k=1,M，递推计算复数

P(kΔτ,mΔω)=P[(k-1)Δτ,mΔω](d₁+jd₂)

式中j是单位虚数，当k=1时P[(k-1)Δτ,mΔω]=F_n(mΔω)；

对每个k，计算正整数

若m_k>15，则令m_k=15；再对P(kΔτ,mΔω)的实部进行累加，得到新的时间序列

${\stackrel{\&OverBar;}{f}}_n(\mathrm{k\&Delta;\&tau;},Q_i)=\underset{m=l_1}{\overset{l_2+m_k}{\mathrm{\&Sigma;}}}\mathrm{Re}\left\{P(\mathrm{k\&Delta;\&tau;},\mathrm{m\&Delta;\&omega;})\right\}$

式中Re代表取实部；完成对全部N个Q_i和选定的水平样点循环，就在每个水平样点处得到一组随旅行时和Q变化的二维道集

3.根据权利要求1所述的一种提高地震成像分辨率的方法，其特征在于：所述的利用步骤D得到的二维道集和基于测井资料的合成地震记录，使用谱比的频率导数指标，得到全部浅部参考点处的等效Q值是这样实现的：（1）根据已有的井资料对选定的浅部参考点分组，按与井点位置的直线距离分为与已有井资料的个数相同的组；对每个井和每组浅部参考点，选取该组中与井直线距离最近的浅部参考点为基准点，令ω_r=0.6(ω_p+ω_m)，其中ω_p和ω_m分别是偏移叠加数据体的主频和有效频带上限，单位是弧度/秒；用以ω_r为主频的雷克子波，利用测井资料得到以基准点的旅行时为中点的旅行时时窗内的合成地震记录，用基准点的对应时窗内的随Q变化的波形与合成地震记录对比，按照波形相似程度和合理性，确定基准点处的等效Q值；（2）在每个选定的水平样点，对浅部参考点处旅行时时窗内的随Q变化的波形做短时傅立叶变换，计算傅立叶变换的模，在有效频带内，求模的自然对数；模的自然对数由如下方法计算：令有效频带内，随频率变化的模D(ω)的平均值是实数D₀，当D(ω)≥0.22D₀，直接计算其对数ln[D(ω)]；当D(ω)<0.22D₀时，计算实数d=5D(ω)/D₀，进而求得自然对数为

$\ln \left[\mathrm{D\&omega;}\right]=\ln \left(0.2D_0\right)+(d-1)(1-\frac{1}{2}(d-1)+\frac{1}{3}{(d-1)}^2)$

（3）对计算得到的每组对数做中值滤波处理；（4）对每组浅部参考点，以基准点处求得的等效Q值对应的一组对数值作为参考值，将其它浅部参考点的随Q变化的每组对数值分别与该参考值相减；对每个Q值对应的随频率变化的相减结果，在有效频带内确定计算频率导数的有效区间，要求区间内突变的峰值是正、负成对出现的；对相减结果做高斯平滑，然后用中心差分法计算其对频率的一阶导数，对得到的导数值做中值滤波，然后求导数的算术平均值，平均值的绝对值就是确定等效Q值的单一的频率导数指标；（5）以水平坐标为水平面，等效Q值为纵坐标，在每个水平样点处，在指定的横向窗口内，将频率导数指标做成随Q值变化的曲线，根据频率导数指标接近零和等效Q值的横向连续性，得到各个浅部参考点处的等效Q值。

4.根据权利要求1所述的一种提高地震成像分辨率的方法，其特征在于：所述的用步骤D得到的二维道集和步骤E得到的浅部参考点的等效Q值，确定全部Q值拾取点处的等效Q值是这样实现的：（1）在每个水平样点，根据步骤E得到的浅部参考点的等效Q值，从随旅行时和Q值变化的二维道集中拾取该等效Q值对应的浅部参考点处的短时波形，并按照权利要求3中描述的方法求得中值滤波处理后的一组对数值，定义该组对数值作为参考值；（2）对各个Q值拾取点处随Q变化的波形，按照权利要求3中描述的方法得到中值滤波处理后的对数值组；（3）把各个Q值拾取点的随Q变化的对数值组分别与浅部参考点的参考值相减；（4）对每个随频率变化的相减结果，在有效频带内确定计算频率导数的有效区间，要求区间内突变的峰值是正、负成对出现的；对相减结果做高斯平滑，然后用中心差分法计算其对频率的一阶导数，对得到的导数值做中值滤波，然后求导数的算术平均值，平均值的绝对值就是对应等效Q值的频率导数指标；（5）形成水平样点处不同时间深度上随Q值变化的频率导数指标图；该图做法如下：在各个Q值拾取点的旅行时处，在指定的纵向窗口内，将频率导数指标做成随Q值变化的曲线；（6）基于这一图形，在频率导数指标接近零的邻域，考虑等效Q值的横向连续性并保证等效Q值应随旅行时增加而增加，直接拾取等效Q值；（7）对每个水平样点采用相同步骤，得到全部Q值拾取点处的等效Q值。

5.根据权利要求1所述的一种提高地震成像分辨率的方法，其特征在于：所述的用步骤G得到的非均匀等效Q值场，对三维偏移叠加数据体的每个成像道应用稳定的非均匀反Q滤波，得到更新的三维成像数据体是这样实现的：根据地震资料的信噪比，确定幅值补偿的最大倍数A，可取为1000.0，该参数无量纲；设偏移叠加数据体的旅行时采样点数为M，而满足快速傅立叶算法的时间采样点数应为M₀；定义无量纲参数

式中各参数与权利要求2中所述的相同；计算参数

其单位为弧度；计算正整数

和

其中Q₀是非均匀等效Q值场的最小值；若n₁≥n₃，则n₁=n₃和n₂=n₃，若n₂≥n₃，则n₂=n₃；设Δτ是偏移叠加数据体的双程旅行时时间采样，单位秒，则傅立叶变换后的角频率采样为Δω=2π/(M₀Δτ)，单位弧度/秒；令偏移叠加数据体的有效频带下限是ω_min，而期望通过粘性吸收补偿达到的频带上限是ω_max，其单位是弧度/秒；通过求整得到正整数l₁和l₂，使得l₁Δω与ω_min和l₂Δω与ω_max的差的绝对值最小；设ω_p是偏移叠加数据体的主频，单位是弧度/秒；

建立无量纲数组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₂

a_i=1.1A,n₂<i≤n₃

对三维偏移叠加数据体的全部成像道循环，对每个成像道应用快速傅立叶变换算法，得到离散的频率域成像道为F_n(mΔω)，其中整数n是成像道的序号，整数m=l₁,l₂+15是离散的频率样点序号；对频率域成像道的高频端应用平滑处理，即

F_n(mΔω)=F_n(mΔω)exp[-0.06(m-l₂)²]，m=l₂+1,l₂+15

计算

和

其单位为弧度；计算无量纲参数 $c_1=\ln \left(\frac{\mathrm{\&Delta;\&omega;}}{{\mathrm{\&omega;}}_p}\right);$

对旅行时采样点循环，k=1,M，从非均匀等效Q值场中拾取该采样点处的等效Q值，用q₀对等效Q值的到数取整，得整数n_q；在每一旅行时采样点，对频率循环，m=l₁,l₂，计算弧度c=m·k(e-c₂(c₁+ln(m))和正整数n₀=m·k；在数组a_i中拾取第n₀个元素值a(n₀)，计算无量纲参数d₁=a(n₀)cos(c)和d₂=a(n₀)sin(c)，再计算复数

P(kΔτ,mΔω)=F_n(mΔω)(d₁+jd₂)

在上述频率循环中记录主频ω_p对应的无量纲补偿因子a(ω_p)和最大频率对应的补偿因子a_max，其中据此计算正整数 $m_k=\mathrm{int}\left(\sqrt{\frac{1}{0.06}\ln (a_{\max }/a\left({\mathrm{\&omega;}}_p\right))+25}\right);$ 若m_k>15，则令m_k=15；

继续对频率循环，m=l₂+1,l₂+m_k，计算弧度c=m·k(e-c₂(c₁+ln(m))，计算无量纲参数d₁=a_maxcos(c)和d₂=a_maxsin(c)，再计算复数

P(kΔτ,mΔω)=F_n(mΔω)(d₁+jd₂)

对P(kΔτ,mΔω)的实部进行累加，得到新的时间序列

${\stackrel{\&OverBar;}{f}}_n\left(\mathrm{k\&Delta;\&tau;}\right)=\underset{m=l_1}{\overset{l_2+m_k}{\mathrm{\&Sigma;}}}\mathrm{Re}\left\{P(\mathrm{k\&Delta;\&tau;},\mathrm{m\&Delta;\&omega;})\right\}$

式中Re代表取实部，

即是粘性吸收补偿后的成像道；完成对全部成像道的循环，就得到更新的三维成像数据体。

Images