技术领域
[0001] 本
发明涉及一种在石油勘探过程中使用的地震
数据处理方法，特别是关于一种曲波域中
地震波衰减补偿方法。
背景技术
[0002] 在石油勘探过程中，
地层对地震波的吸收造成了
震源子波时变，使资料的
分辨率降低，特别是深层
信号，由于地层的低通滤波作用，高频信号衰减严重。因此，地表记录到的不同深度反射信号及其频带，已经不能反映地下的真实情况，影响了后续地震资料的处理和解释工作，也影响了石油勘探结果的准确性。为了提高地震资料的分辨率，就必须对地震波衰减进行吸收补偿处理。
[0003] 目前，地震数据处理中采用的地震波衰减补偿方法主要有常规球面发散与吸收补偿方法（牟永光.地震资料数字处理方法.石油工业出版社.1981.3）、常规反Q滤波（Futterman W I.Dispersive body wave[J].Geophysics Res.1962.67：5279-5291）、常规谱白化和单道或多道反褶积（Robinson.E.A.Principles of digital Wiener filtering.Geophysical prospecting[J].1967.15(3)：311-333）等。常规球面发散与吸收补偿方法仅补偿地震波振幅随时间变化的衰减，而不能补偿随
频率变化的衰减。常规反Q（Q为品质因子，表示地层吸收的量）滤波尽管可以补偿由于地震波衰减引起的振幅减弱与速度频散，但要求品质因子Q已知，而品质因子Q值不易求准，且反Q滤波降低了数据的
信噪比。常规谱白化是在假设反射系数白噪的条件下，通过对分频数据进行自动增益来补偿振幅随时间和频率的吸收衰减，因此难以保持子波和振幅
能量的信息，而且会相应增强炮间的干扰能量。单道或多道反褶积是基于时窗统计内平均条件下的时频补偿，难以较好的逐点补偿振幅衰减和频率吸收。
[0004] 近年来，国内学者提出几种基于时频分析的地层吸收补偿方法，如基于短时
傅立叶变换地层吸收补偿方法（白桦，李鲲鹏.基于时频分析的地层吸收补偿[J].石油地球物理勘探.1999，34（6）：642-648）、基于小波包分解的地层吸收补偿方法（李鲲鹏，李衍达，张学工.基于小波包分解的地层吸收衰减[J].地球物理学报.2000.43(4)：542-549）和基于广义S变换的地层吸收补偿方法（刘喜武，年静波，刘洪.基于广义S变换的吸收衰减补偿方法[J].石油物探.2006.45（1）：9-14）等。这三种方法利用时频分析和地震波衰减特性求得地层吸收系数，可以一定程度上消除地震波衰减对地震记录的时变影响，改善地震资料的分辨率，但是上述几种方法对信号信噪比有一定要求，在补偿地震波衰减的同时也增强了噪声和随机干扰的能量。
[0005] 综上所述，现有的技术对地震波衰减补偿的处理均有不足之处，如常规球面发散与吸收衰减补偿方法仅能补偿随时间变化的能量衰减，而不能补偿随频率的吸收衰减；常规反Q滤波方法Q值不易求准，且处理的同时降低了地震信号的信噪比；常规谱白化方法难以保持振幅；单道或多道反褶积方法不能点对点补偿吸收衰减变化；基于时频分析的吸收补偿方法对地震数据信噪比有一定要求，而且在补偿地震波衰减的同时相应提高了噪声和随机干扰的能量。
发明内容
[0006] 针对上述问题，本发明的目的是提供一种受地震数据信噪比限制较小且只针对有效信号进行补偿，使得不同深度反射波的
波形基本一致的曲波域中地震波衰减补偿方法。
[0007] 为实现上述目的，本发明采取以下技术方案：一种曲波域中地震波衰减补偿方法，其包括以下步骤：
[0008] 1）收集地震信号并进行预处理，得到地震的时域信号x(t)，其中t为二维空间中时域信号x(t)的空间
位置变量；
[0009] 2）将时域信号x(t)做曲波变换，得到曲波域信号ψj,l(ωj,Tk)，并逐点补偿地震波衰减，其中，j是曲波变换的尺度因子且j＝1,2,...,N；l是曲波变换的
角度因子且l＝1,2,...,N，ωj是频带，Tk是反射时间；
[0010] 3）将曲波域信号ψj,l(ωj,Tk)进行平滑处理得到曲波系数Φj,l(ωj,Tk)，减小反射波相互重叠的影响；
[0011] 4）对曲波系数Φj,l(ωj,Tk)进行
阈值迭代处理，得到去除噪声后的信号Ψj,l(ωj,Tk)，减少随机噪声干扰；
[0012] 5）逐点递推求每一频段补偿角度内去除噪声后的信号Ψj,l(ωj,Tk)的深浅层能量比δj,ll(ωj,Tk)，使高低频信号具有相同的深浅层能量比；
[0013] 6）求取补偿角度内的所有频段补偿因子系数矩阵Dj,ll(ωj,Tk)，使得补偿角度系数范围内的各频段深浅层能量的比值相同；
[0014] 7)将所有频段补偿角度内的补偿因子系数矩阵Dj,ll(ωj,Tk)的倒数加权，得到曲波域加权后的信号Cj,l(ωj,Tk)，实现相应频段相应角度信号的分频定向补偿；
[0015] 8)对曲波域加权后的信号Cj,l(ωj,Tk)进行反曲波变换，得到分频定向补偿后的时域信号xx(t)。
[0016] 所述步骤1）中，时域信号x(t)通过以下方法获得：对叠前地震资料进行处理得到二维
叠加剖面数据，将所得二维地震数据按照
共中心点记录方式排列，得到地震的时域信号x(t)。
[0017] 所述步骤2）中，曲波域信号ψj,l(ωj,Tk)通过以下方法得到：
[0018] 在二维空间中，设时域信号x(t)的空间位置变量是t，频率域变量是ω，对应的极坐标是(r,θ)，半径窗是W(r)，角时窗是V(t)；半径窗W(r)在 范围内
支撑，角时窗V(t)在t∈(-1,1)范围内支撑，且满足容许性条件：
[0019]
[0020]
[0021] 对于每一个j≥1，利用傅立叶变换定义频率窗函数：
[0022]
[0023] 其中， 是j/2的整数部分，极角θ范围是0～2π；频率窗函数Uj(r,θ)的支撑区间为一楔形区域，受半径窗W(r)和角时窗V(t)的支撑区间所限制；定义 其中，频率窗函数Uj(r,θ)为Uj(ω)的极坐标形式，在j＝1时， 为时域信号x(t)在频率域对应的曲波变换“母”函数， 经旋转平移得到其它尺度（j＝2,3,...,N）和所有角度（l＝1,2,...,N）上的曲波基函数 其中平移参数k＝(k1,k2)∈Z2表示位置，k1表示纵坐标的位置，k2表示横坐标的位置，Z2表示数据在二维空间上的分布区间的位置；
[0024] 定义旋转角度序列 范围在0＜θl＜2π；平移参数k＝(k1,k2)∈Z2处的曲波变换为：
[0025]
[0026] 其中， 表示θ弧度的旋转， 是Rθ的逆，且 时域信号x(t)做曲波变换后，转换为曲波域信号ψj,l(ωj,Tk)：
[0027]
[0028] 所述步骤3）中，曲波系数Φj,l(ωj,Tk)通过以下方法得到：
[0029] 采用Parzen窗作为一维平滑窗函数：
[0030]
[0031] 曲波域信号ψj,l(ωj,Tk)与平滑窗函数d(τ)相褶积，得到曲波系数Φj,l(ωj,Tk)：
[0032]
[0033] 其中，T是数据时间长度，且T＝5s。
[0034] 所述步骤4）中，去除噪声后的信号Ψj,l(ωj,Tk)通过以下方法得到：
[0035] 利用曲波域阈值迭代衰减随机噪声影响的基本原理描述为求解如下的优化反演问题：
[0036]
[0037] 其中，Pε表示求解ε的优化反演问题，ε表示极小值，A＝CT是反曲波变换，x是曲波系数Φj,l(ωj,Tk)，y是所述步骤1）中含噪声的地震的时域信号x(t)，ε是一任意小量， 是最终输出结果；
[0038] 求解上述优化问题，反演一组满足经反曲波变换后与原始数据差的二范数小于某一噪声参量的l1最小的曲波系数Φj,l(ωj,Tk)，阈值迭代去噪声后的结果
[0039] Pε问题进一步优化为如下的优化反演问题：
[0040]
[0041] Pλ表示求解λ的优化反演问题；
[0042] 求解Pλ问题中采用迭代阈值法，其迭代序列：
[0043] x＝Tλ(x+AT(y-Ax))
[0044] Tλ(x):＝sgn(x)·max(0,|x|-|λ|)
[0045] 其中，Tλ是软阈值函数；
[0046] 阈值迭代去除噪声后，输出的结果 表示为去除噪声后的信号Ψj,l(ωj,Tk)。
[0047] 所述步骤5）中，深浅层能量比δj,ll(ωj,Tk)通过以下方法得到：
[0048] 去除噪声后的信号Ψj,l(ωj,Tk)在曲波域中的地震记录为Ψ(ω)：
[0049]
[0050] 其中，rk(k＝0,1......K)是反射系数序列，K是
地震道的道数，Tk是一个反射系数rk的反射时间， 是Q滤波的振幅响应，其中，Q为品质因子，表示地层吸收的量， Qeq(Tk)为反射时间为Tk处的等效Q值，A0(ω,0)是初
始时刻的等效振幅衰减，ω(ω)是曲波域中的地震子波；
[0051] 经过所述步骤3）平滑处理后，地震记录上相邻反射面反射波互不重叠，则每一频段特定补偿角度内的深浅层能量比δj,ll(ωj,Tk)：
[0052]
[0053] Ψj,ll(ωj,Tk)是每个频段补偿角度内的曲波系数矩阵，Ψj,ll(ωj,0)是每个频段补偿角度内的曲波系数矩阵的浅层系数值； 是不同频段补偿角度系数的振幅响应，ll是曲波变换补偿角度系数，曲波变换补偿角度系数ll的选取原则为在时间域的角度范围只要不小于地震数据同相轴的角度范围即可，此角度系数一般选为±45°。
[0054] 所述步骤6）中，补偿因子系数矩阵Dj,ll(ωj,Tk)通过以下方法得到：
[0055] 补偿因子系数矩阵Dj,ll(ωj,Tk)是高频信号（j＝3,4,5）特定补偿角度内的深浅层能量比δj,ll(ωj,Tk)同基准频带信号对应补偿角度内的深浅层能量比δjj,lll(ωjj,Tk)之比：
[0056]
[0057] 其中，j＝3,4,5，jj是基准频段的尺度值，lll是高频段曲波变换补偿角度系数ll对应的基准频段的对应角度系数值；不同尺度角度间的对应关系为 曲波变换补偿角度系数ll的值是对方程右边取整数。
[0058] 所述步骤7）中，曲波域加权后的信号Cj,l(ωj,Tk)通过以下方法得到：
[0059] 曲波域加权后的信号为Cj，ll(ωj,Tk)，补偿角度系数之外的信号为Ψj，l-ll(ωj,Tk)，曲波变换的角度因子l为曲波变换的尺度因子j上0～2π方向，经过分频定向补偿处理后的曲波域加权后的信号Cj,l(ωj,Tk)：
[0060]
[0061] Cj,l(ωj,Tk)＝Cj,ll(ωj,Tk)+Ψj,l-ll(ωj,Tk)。
[0062] 所述步骤8）中，分频定向补偿后的时域信号xx(t)通过以下方法得到：
[0063] 对曲波域加权后的信号Cj,l(ωj,Tk)进行反曲波变换，得到分频定向补偿处理后的时域信号xx(t)；
[0064] 由曲波变换原理得到反曲波变换的实现过程：
[0065]
[0066] 其中，M表示指标集，φj,l,k表示由指标(j,θ,k)确定的曲波族。
[0067] 本发明由于采取以上技术方案，其具有以下优点：1、由于地震波振幅的衰减是逐点时变的，本发明通过曲波变换将时域信号变换得到相应的频域信号，因此能在曲波域逐点补偿地震波的能量衰减。2、本发明的曲波变换能多层次划分频带，将高频成分进一步细分，因此具有比常用的时频变换（如傅立叶变换和
小波变换等）更高的时频分辨率，能够更好地补偿高频衰减，提高地震信号的分辨率和薄层解释的准确性。3、本发明采用阈值迭代法进行去噪，设定补偿角度提取有效信号进行处理，实现了对地震波衰减的定向补偿，即只定向补偿地震波有效信号的衰减，而未增强噪声和随机干扰的能量。4、本发明由于对曲波变换数据的信噪比没有太大要求，因此不需要在处理前进行去噪处理，不但可以减少有效信号预前去噪过程中的损失，而且还可以降低地震波衰减分频定向补偿处理的工作量。5、由于本
申请人已经拥有有Matlab（Matrix&Laboratory，矩阵实验室）和Linux(一种网络
操作系统)两种环境下的程序包，因此，地球物理数据的处理和解释人员都可以借助程序包进行本发明的地震波衰减分频定向补偿处理，具有较好的使用和推广价值。本发明可以广泛用于石油勘探过程中对各种地震波衰减进行吸收补偿处理的过程中。
附图说明
[0068] 图1是本发明浅层数据处理前（a）后(b)对比示意图
[0069] 图2是本发明中深层数据处理前（a）后(b)对比示意图
[0070] 图3是本发明深层数据处理前（a）后(b)对比示意图
[0071] 图4是本发明数据处理前（a）后(b)第50道地震信号的
频谱分析示意图具体实施方式
[0072] 下面结合附图和
实施例对本发明进行详细的描述。
[0073] 本发明曲波域中地震波衰减补偿方法包括以下步骤：
[0074] 1）收集地震信号并进行预处理，得到地震的时域信号x(t)
[0075] 从磁带上读取地震信号，并对叠前地震资料进行常规处理，得到二维叠加剖面数据，将所得二维地震数据按照共中心点记录方式排列，此时所得信号为地震的时域信号x(t)。如图1(a)、图2(a)、图3(a)所示，分别为地震的浅、中、深层的时域信号数据，其中，各图的横轴CDP表示采用110道数据叠加而成，纵轴t表示时间，单位为秒，图1(a)截取的是地震波衰减分频定向补偿处理前后时间轴上0.7s到1.8s处的浅层数据；图2(a)截取的是地震波衰减分频定向补偿处理前后时间轴上2.0s到3.1s处的中深层数据；图3(a)截取的是地震波衰减分频定向补偿处理前后时间轴上3.4s到4.5s处的深层数据。
[0076] 2）将时域信号x(t)做曲波变换（Curvelet变换），得到曲波域信号ψj,l(ωj,Tk)，并逐点补偿地震波衰减
[0077] 地震波在地层中传播时，振幅总体趋势呈指数规律衰减，且在高频和低频中具有不同的衰减趋势，即地震波振幅的衰减是逐点时变的。因此，将时域信号x(t)变换成曲波域信号ψj,l(ωj,Tk)才可以真正实现地震波衰减的逐点补偿，其中，j是曲波变换的尺度因子，即频率域变量为ω的频带数，范围是j＝1,2,...,N；l是曲波变换的角度因子，表示曲波变换的尺度因子j内的曲波变换角度数，范围是l＝1,2,...,N，ωj是不同频带，Tk是反射时间。
[0078] 在二维空中，设时域信号x(t)的空间位置变量是t，频率域变量是ω，对应的极坐标是(r,θ)，半径窗是W(r)，角时窗是V(t)。半径窗W(r)在 范围内支撑，角时窗V(t)在t∈(-1,1)范围内支撑，且满足容许性条件：
[0079]
[0080]
[0081] 对于每一个j≥1，利用傅立叶变换定义频率窗函数：
[0082]
[0083] 其中， 是j/2的整数部分，极角θ范围是0～2π。频率窗函数Uj(r,θ)的支撑区间为一楔形区域，受半径窗W(r)和角时窗V(t)的支撑区间所限制。定义 其中，频率窗函数Uj(r,θ)为Uj(ω)的极坐标形式，在j＝1时， 为时域信号x(t)在频率域对应的曲波变换“母”函数， 经旋转平移得到其它尺度（j＝2,3,...,N）、所有角度（l＝1,2
2,...,N）上的曲波基函数 其中平移参数k＝(k1,k2)∈Z表示位置，k1表示纵坐标的位置，k2表示横坐标的位置，Z2表示数据在二维空间上的分布区间的位置。
[0084] 定义旋转角度序列 范围在0＜θl＜2π。平移参数k＝(k1,k2)∈Z2处的曲波变换为：
[0085]
[0086] 其中， 表示θ弧度的旋转， 是Rθ的逆，且 时域信号x(t)做曲波变换后，转换为曲波域信号ψj,l(ωj,Tk)：
[0087]
[0088] 本实施例中，曲波变换的尺度因子j为5，曲波变换的角度因子l为16，即时域信号x(t)经曲波变换后，对应的曲波域信号ψj,l(ωj,Tk)有5个频带，第一个频带有1个方向，方向间隔为2π；第二个频带共有16个方向，方向间隔为2π/16；第三个频带共有32个方向，方向间隔为2π/32；第四个频带共有32个方向，方向间隔为2π/32；第五个频带共有64个方向，方向间隔为2π/64。本实施例中考虑第二频带。
[0089] 3）将曲波域信号ψj,l(ωj,Tk)进行平滑处理得到曲波系数Φj,l(ωj,Tk)，减小反射波相互重叠的影响
[0090] 实际地震记录中，相邻反射界面的反射波是有重叠的，重叠部分的地震信号会产生能量畸变和
相位失真。实际数据处理中需要对信号进行平滑，减小反射波相互重叠的影响。
[0091] 对信号的平滑需要通过加平滑窗的方式来解决，因此平滑窗的选取十分重要，本发明方法选择满足四阶累积量对称性质和归一化条件，且满足在所有频率内函数值不小于0的Parzen窗（ParzenWindow，又名核
密度估计，Kernel Density Estimation）来作为一维平滑窗函数对曲波域信号ψj,l(ωj,Tk)进行平滑，其中平滑窗函数Parzen窗的表达式为：
[0092]
[0093] 曲波域信号ψj,l(ωj,Tk)与平滑窗函数d(τ)相褶积，可以得到曲波系数Φj,l(ωj,Tk)：
[0094]
[0095] 其中，T是数据时间长度，且T＝5s。
[0096] 褶积后得到曲波系数Φj,l(ωj,Tk)比单个曲波域信号ψj,l(ωj,Tk)和平滑窗函数d(τ)都光滑。特别地，当平滑窗函数d(τ)为具有紧致集的光滑函数时，且曲波域信号ψj,l(ωj,Tk)为局部可积时，它们的褶积曲波系数Φj,l(ωj,Tk)也是光滑函数。利用这一性质，对原始曲波域信号ψj,l(ωj,Tk)进行平滑可以有效修正地震信号在不同反射界面的反射波相互重叠处的能量畸变和相位失真，使反射波相互重叠的影响降至最低。
[0097] 4）对曲波系数Φj,l(ωj,Tk)进行阈值迭代处理，得到去除噪声后的信号Ψj,l(ωj,Tk)，减少随机噪声干扰
[0098] 由步骤2）中的公式（3）可以导出曲波基函数 曲波基函数 的支撑区间在频率域满足长度≈宽度×宽度，这种“长条形”的支撑区间是方向性的一种体现。长条形的“
各向异性”基使得曲波域信号ψj,l(ωj,Tk)在时域信号x(t)与曲波基函数 方向大致相同时取得较大的曲波系数Φj,l(ωj,Tk)，时域信号x(t)与曲波基函数 方向垂直或近似垂直时取得较小的曲波系数Φj,l(ωj,Tk)。
[0099] 地震数据同相轴可以近似认为是曲线反射层，因此有效信号（即地震信号同相轴）曲波系数Φj,l(ωj,Tk)较大，而没有确定方向的噪声和随机干扰则系数较小。为了衰减噪声的影响，保留有效信号，针对曲波系数Φj,l(ωj,Tk)进行阈值迭代处理，即选取一个阈值λ，它决定了曲波系数Φj,l(ωj,Tk)经过
阀值迭代处理后保留的比例系数。略去较小系数，保留较大系数，即可以衰减噪声的影响，避免对噪声过多补偿而降低地震数据的信噪比。
[0100] 利用曲波域阈值迭代衰减随机噪声影响的基本原理可以描述为求解如下的优化反演问题：
[0101]
[0102] 其中，Pε表示求解ε的优化反演问题，ε表示极小值，A＝CT是反曲波变换，x是曲波系数Φj,l(ωj,Tk)，y是步骤1）中含噪声的地震的时域信号x(t)，ε是一任意小量， 是最终输出结果。
[0103] 求解上述优化问题，就是要反演一组满足经反曲波变换后与原始数据差的二范数小于某一噪声参量的l1最小的曲波系数Φj,l(ωj,Tk)，最终的输出结果即为阈值λ迭代去噪声后的结果
[0104] Pε问题可以进一步优化为如下的优化反演问题：
[0105]
[0106] Pλ表示求解λ的优化反演问题。
[0107] 求解Pλ问题中采用迭代阈值法，其迭代序列：
[0108] x＝Tλ(x+AT(y-Ax)) （10）
[0109] Tλ(x):＝sgn(x)·max(0,|x|-|λ|) （11）
[0110] 其中，Tλ是软阈值函数。
[0111] 计算过程中，阈值λ随着循环次数的增加而降低，不断迭代更新曲波系数Φj,l(ωj,Tk)，直到取得满意的结果为止。
[0112] 阈值λ迭代去除噪声后，输出的结果 表示为去除噪声后的信号Ψj,l(ωj,Tk)。
[0113] 5）逐点递推求每一频段补偿角度内去除噪声后的信号Ψj,l(ωj,Tk)的深浅层能量比δj,ll(ωj,Tk)，使高低频信号具有相同的深浅层能量比
[0114] 如果地震波在地层中传播时，没有经过衰减，则不同深度的反射波具有相同的振幅谱。对地震记录做时频变换，不同频段的信号应具有相似的能量-时间分布关系，即所有频段内深层反射波与浅层反射波能量之比相同。
[0115] 但是，由于地层的低通滤波作用，高频信号的衰减远远大于低频信号，使得不同频段内的地震信号具有不同的能量衰减趋势，如果对高频信号乘以一个补偿因子矩阵，使得高低频信号具有相同的深浅层能量比，即补偿了地震波的衰减。若要求取有效信号的补偿因子系数矩阵，则首先要求取曲波域所有频段特定补偿角度范围内的深浅层能量比。
[0116] 去除噪声后的信号Ψj,l(ωj,Tk)在曲波域中的地震记录为Ψ(ω)：
[0117]
[0118] 其中，rk(k＝0,1......K)是反射系数序列，K是地震道的道数，Tk是一个反射系数rk的反射时间， 是Q滤波的振幅响应，其中，Q为品质因子，表示地层吸收的量， Qeq(Tk)为反射时间为Tk处的等效Q值，A0(ω,0)是初
始时刻的等效振幅衰减，ω(ω)是曲波域中的地震子波。
[0119] 经过步骤3）平滑处理后，可以认为地震记录上相邻反射面反射波互不重叠，则每一频段特定补偿角度内的深浅层能量比δj,ll(ωj,Tk)：
[0120]
[0121] Ψj,ll(ωj,Tk)是每个频段补偿角度内的曲波系数矩阵，Ψj,ll(ωj,0)是每个频段补偿角度内的曲波系数矩阵的浅层系数值。 是不同频段补偿角度系数的振幅响应，ll是曲波变换补偿角度系数，可以根据实际情况设定，应用曲波变换补偿角度系数ll可以实现对地震波衰减的定向补偿，即在此曲波变换补偿角度系数ll范围内的曲波域地震信号会进行处理，而此曲波变换补偿角度系数ll范围外保持不变。曲波变换补偿角度系数ll的选取原则为在时间域的角度范围只要不小于地震数据同相轴的角度范围即可。对于叠后同相轴
水平或近水平的地震数据，此角度系数一般选为±45°。
[0122] 本实施例中，第一频段为一个矩阵，故无角度，第二、三、四频段的补偿角度系数为±45°，第五频段的系数为±60°，这样可以在有效补偿中高频信号吸收衰减的同时，尽量多的补偿高频信号的衰减。
[0123] 6）求取补偿角度内的所有频段补偿因子系数矩阵Dj,ll(ωj,Tk)，使得补偿角度系数范围内的各频段深浅层能量的比值相同
[0124] 由于大地具有低通滤波作用，低频成分衰减较慢，可以认为低频信号基本未衰减，故选取为基准频带。为使补偿角度系数范围内的各频段深浅层能量的比值相同，因此需要求出补偿因子系数矩阵Dj,ll(ωj,Tk)。
[0125] 补偿因子系数矩阵Dj,ll(ωj,Tk)是高频信号（j＝3,4,5）特定补偿角度内的深浅层能量比δj,ll(ωj,Tk)同基准频带信号对应补偿角度内的深浅层能量比δjj,lll(ωjj,Tk)之比：
[0126]
[0127] 其中，j＝3,4,5，jj是基准频段的尺度值，lll是高频段曲波变换补偿角度系数ll对应的基准频段的对应角度系数值。不同尺度角度间的对应关系为 曲波变换补偿角度系数ll的值是对方程右边取整数。
[0128] 7)将所有频段补偿角度内的补偿因子系数矩阵Dj,ll(ωj,Tk)的倒数加权，得到曲波域加权后的信号Cj,l(ωj,Tk)，以实现相应频段相应角度信号的分频定向补偿[0129] 由步骤6）中的公式（14）可知，补偿因子系数矩阵Dj,ll(ωj,Tk)是高频信号（j＝3,4,5）特定补偿角度内的深浅层能量比δj,ll(ωj,Tk)同基准频带信号对应补偿角度内的深浅层能量比δjj,lll(ωjj,Tk)之比，那么使各个频段相应补偿角度内每个频段的曲波系数矩阵Ψj,ll(ωj,Tk)倒数加权，使得各频段深浅层能量的比值相同，分频定向补偿功能即可实现。
[0130] 加权后的信号为Cj,ll(ωj,Tk)，补偿角度系数之外的信号为Ψj,l-ll(ωj,Tk)，曲波变换的角度因子l为曲波变换的尺度因子j上0～2π方向，故经过分频定向补偿处理后的曲波域加权后的信号Cj,l(ωj,Tk)：
[0131]
[0132] Cj,l(ωj,Tk)＝Cj,ll(ωj,Tk)+Ψj,l-ll(ωj,Tk（) 16）
[0133] 8)对曲波域加权后的信号Cj,l(ωj,Tk)进行反曲波变换，得到分频定向补偿后的时域信号xx(t)
[0134] 曲波域加权后的信号为Cj,l(ωj,Tk)，对其进行反曲波变换即可以得到分频定向补偿处理后的时域信号xx(t)。
[0135] 由曲波变换原理可以得到反曲波变换的实现过程：
[0136]
[0137] 其中，M表示指标集，φj,l,k表示由指标(j,θ,k)确定的曲波族。
[0138] 如图1所示，对比图1（a）和图1（b）：CDP在1-110区域内，时间轴上0.9s-1.1s之间，可以看出同相轴（地震记录上各道振动相位相同的极值的连线，即图上能量较强的黑色部分）连续性明显增强，且比处理前更加清晰，即表示处理后的浅层数据分辨率有了明显提升。对比图1（a）和图1（b）：CDP在1-110区域内，时间轴上1.3s-1.6s之间，可以看出处理后的数据除了同相轴连续性大大增强外，有效信号的能量也有少许增强（CDP在1-110区域内，时间轴上1.3s-1.4s之间，处理后同相轴的能量有一定加强）。处理前杂乱无章的信号经过处理后能量明显降低（如对比图1（a）和图1（b）：CDP在1-110区域内，时间轴1.1s-1.2s之间），说明处理后随机噪声得到了明显衰减。
[0139] 由图1（a）和图1（b）对比结果可知：曲波域地震波衰减分频定向补偿处理后浅层地震记录同相轴的分辨率有明显提升，连续性有很大增强，有效信号的能量有少许加强，随机噪声得到了明显衰减。
[0140] 如图2所示，对比图2（a）和图2（b）：CDP在1-110区域内，时间轴上2.4s-2.5s之间，可以看出处理后的数据除了同相轴连续性有较大增强外，一些原本模糊不清的同相轴也变得清晰起来，这说明处理后中深层数据的分辨率有了较大改善。处理后同相轴的能量相比处理前都得到了一定增强，即有效信号的能量衰减得到了较好的补偿（在CDP从1-110区域内，时间轴2.7s-2.8s之间，处理后同相轴的能量有较大加强），噪声和随机干扰也得到了较好的压制（在CDP从1-110区域内，时间轴2.7s-2.8s之间，杂乱无章的信号得到了较大的衰减）。
[0141] 由图2（a）和图2（b）对比结果可知：处理后的中深层数据同相轴分辨率有较大改善，一些原本模糊的同相轴也变的清晰，噪声和随机干扰得到了较好的压制，有效信号的能量衰减得到了较好的补偿。
[0142] 如图3所示，对比图3（a）和图3（b）：CDP在1-110区域内，时间轴在3.6s-3.8s之间，可以看出处理后的数据同相轴的连续性和分辨率有较大增强外，原本模糊不清的同相轴也清晰可辨，噪声和随机干扰的能量也得到了较好的衰减，这说明处理后深层数据的分辨率和信噪比有了极大提高。由于分辨率的提高，许多地质薄层也有了清楚的表示（如两图CDP在1-110区域内，时间轴上4.0s-4.1s），可以看出处理后，原先混在一起的地质薄层得以分离开来，这对于后续的地质解释和储层圈定工作有较大的意义。
[0143] 由图3（a）和图3（b）对比结果可知：处理后地震信号的分辨率和信噪比有极大的提高，模糊不清的同相轴变的清晰可辨，地质薄层也有了清楚的表示。
[0144] 如图4所示，对比图4（a）和图4（b）可知：频率相同时，处理后的数据振幅普遍高于补偿前，尤其是中高频（40Hz-60Hz）处，处理后地震信号的振幅值有了明显的提升，这说明使用本发明可以衰减严重的中高频信号得到有效补偿，频带展宽。
[0145] 如图1～4所示的结果可以看出：本发明可以使地震波衰减得到较好补偿，不同深度反射波的波形基本一致，中深层反射波的分辨率和信噪比有明显提高，高频加强，频带展宽。
[0146] 综上所述，可以将本发明归纳为以下步骤：
[0147] 1）收集地震信号并进行预处理，得到地震的时域信号x(t)；
[0148] 2）将地震的时域信号x(t)做曲波变换，得到曲波域信号ψj,l(ωj,Tk)，并逐点补偿地震波衰减；
[0149] 3）将曲波域信号ψj,l(ωj,Tk)进行平滑处理得到曲波系数Φj,l(ωj,Tk)，减小反射波相互重叠的影响；
[0150] 4）对曲波系数Φj,l(ωj,Tk)进行阈值迭代处理，得到去除噪声后的信号Ψj,l(ωj,Tk)，减少随机噪声干扰；
[0151] 5）逐点递推求每一频段补偿角度内去除噪声后的信号Ψj,l(ωj,Tk)的深浅层能量比δj,ll(ωj,Tk)，使高低频信号具有相同的深浅层能量比；
[0152] 6）求取补偿角度内的所有频段补偿因子系数矩阵Dj,ll(ωj,Tk)，使得补偿角度系数范围内的各频段深浅层能量的比值相同；
[0153] 7)将所有频段补偿角度内的补偿因子系数矩阵Dj,ll(ωj,Tk)的倒数加权，得到曲波域加权后的信号Cj,l(ωj,Tk)，实现相应频段相应角度信号的分频定向补偿；
[0154] 8)对曲波域加权后的信号Cj,l(ωj,Tk)进行反曲波变换，得到分频定向补偿后的时域信号xx(t)。
[0155] 由上述步骤可知，本发明只要确定曲波变换的尺度因子j、曲波变换的角度因子l，阈值λ和曲波变换补偿角度系数ll即可实现曲波域中地震波衰减补偿。
[0156] 上述各实施例仅用于说明本发明，凡是在本发明技术方案的
基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。
