0 引言

在地震资料处理中，分辨率一直是评判处理结果的关键指标之一。特别是在开发地震中，分辨率显得尤为重要。如高分辨率处理是薄储层、低幅度构造、小断层等小尺度地质体勘探的关键技术。

常用的提高分辨率方法主要有反褶积、Q补偿、谱白化、小波变换和广义S变换等。

赵玉敏^[1]认为传统反褶积技术无法较大幅度提高地震资料分辨率，同时反褶积方法是针对单道处理，容易破坏同相轴的横向连续性，且在有噪声时处理效果不稳定，因而提出一种信噪比约束下的提高分辨率方法。李红彩等^[2]先应用地表一致性稳健反褶积提高叠前地震数据分辨率，再对该叠后数据应用短时窗傅里叶变换，综合应用叠前、叠后两种提高分辨率方法，为墩塘地区岩性油气藏或隐蔽油气藏勘探提供物探支撑。魏忠宇^[3]提出一种径向道域反褶积方法，并应用于实际地震资料的高分辨率处理。

余连勇等^[4]基于地下介质的水平层状假设，根据叠前地震记录和VSP记录观测系统的特点，给出叠前地震记录反Q滤波和VSP记录(零井源距和非零井源距)反Q滤波的统一算法和公式，并通过模型验证了方法的有效性。董相杰等^[5]根据零井源距VSP下行初至波最大振幅随时间的变化趋势，将地层划分为强衰减和弱衰减两类，提出VSP强衰减地层反Q滤波。张固澜等^[6]基于Futterman提出的地震波振幅衰减和相速度频散表达式，从井震匹配角度出发，阐明了反Q滤波相位补偿的必要性；并且提出一种自适应增益限的反Q滤波振幅补偿方法^[7]。程志国等^[8]采用改进的谱比法技术系列从VSP下行波场准确提取Q值，利用井控约束建立三维Q场，最后使用稳定的反Q滤波对叠前数据进行处理，补偿后的剖面分辨率得到提高。陈树民等^[9]通过引入描述黏性吸收的等效Q值参数，发展了补偿介质吸收的叠前时间偏移方法，即基于地表观测数据建立等效Q值模型，应用补偿因子的光滑性阈值控制介质黏性吸收补偿的稳定性，提高了成像分辨率。孙明等^[10]针对冀东南堡4号中浅层河道砂体识别，采用VSP井驱动Q补偿、高保真分频去噪、分频剩余静校正、井约束反褶积及优势频带反褶积等技术手段，拓宽了频带、提高了信噪比及分辨率。赵秋芳等^[11]基于S变换的谱比法，探讨了获得理想Q估算值的最小含气储层厚度应不小于地震子波的主波长。周衍等^[12]运用稳定化反Q滤波方法对塔河三维地震数据进行提高分辨率处理，同步实施振幅补偿和子波相位校正，展示该方法能提高对碳酸盐岩缝洞储集体特征描述的精度。

王季^[13]利用Hilbert-Huang变换对地震信号进行时频分解，再通过白化滤波器对其Hilbert谱进行谱白化，增强了地震信号时域和频域的分辨率，使剖面更连续和清晰。颜中辉等^[14]对地震记录进行固有模态分解(EMD)，得到不同尺度的IMF分量；再利用常规谱白化方法对各分量根据瞬时频率进行合理的振幅均衡；将均衡后的IMF分量做反变换重构地震记录，从而得到高分辨地震数据。

路鹏飞等^[15]根据薄互层地质特点和地震响应特征，提出一种利用小波变换提高地震资料分辨率的方法，提高了老爷庙油田地震资料分辨率。黄捍东等^[16]基于广义S变换研究提高地震资料分辨率的方法。杨子鹏等^[17]将广义S变换理论与压缩感知理论相结合，提出一种新的叠后地震资料处理方法，并给出补偿计算公式。本文参考了该补偿计算公式，并做了相应改进。

通过上述调研发现，反褶积方法很容易破坏同相轴的横向连续性，且有噪声存在时，该方法很不稳定^[1]。反Q滤波法可补偿地震波振幅衰减，校正地震波的相位特征，从而改善地震资料同相轴的连续性，增强弱反射波的能量，提高地震资料分辨率，但现有的反Q滤波法因其补偿公式的缺陷，补偿后的地震记录会出现吉布斯效应及背景高频噪声被抬升等现象，导致地震资料信噪比降低^[18]。谱白化方法可提高地震资料分辨率，但在提频同时往往会放大噪声。应用小波变换可提高地震资料分辨率，且计算速度非常快。但该类方法易对资料中的低频弱振幅信号的特征造成损害，而这些信号所蕴含的信息往往与地质体的含油气特征息息相关^[19]。

在小波分析理论基础上，近年来相继出现了Ridgelet、Curvelet、Contourlet、Shearlet等变换的多尺度分析方法，其中Shearlet变换以优越的多分辨率、多尺度和多方向性特征，发展迅速，已广泛应用于信号和图像处理领域，并取得了良好的效果。同时，Shearlet变换在地学领域也有广泛应用，如数据重建^[20]、面波压制^[21]、随机噪声去除^[22-26]、河道检测^[27]，被认为是目前对地震信号最优表达方式之一。

已有的提高地震分辨率方法在提高地震信号分辨率的同时往往会放大噪声，降低地震资料的信噪比，需要噪声去除辅助工作。本文根据已有提高地震资料分辨率方法的优缺点，通过Shearlet变换将地震资料转换到Shearlet域，在Shearlet域将有效信号与随机噪声分开；同时，在Shearlet域将有效信号分解成多个尺度、多个分辨率和多个方向性信号。结合Shearlet变换这两个特点，首先舍弃Shearlet域随机噪声的系数，且仅对优势尺度(频带)和代表地层方向有效范围的Shearlet域系数做提频处理。这样既提高了地震资料分辨率，又保持了信噪比。合成数据和实际资料的处理结果表明，本文方法可有效提高叠后地震资料分辨率。

1 Shearlet变换提高地震资料分辨率 1.1 二维Shearlet变换原理

在二维空间，信号f的Shearlet变换为

$ {\rm{S}}{{\rm{H}}_\mathit{\Phi }}(f) = \left\langle {f,{\mathit{\Phi }_{a,s,\mathit{\boldsymbol{p}}}}} \right\rangle $

(1)

式中：a、s、p分别为对应尺度、方向、位置的参数；Φ_{a, s, p}为Shearlet母函数，且满足

$ \varPhi_{a, s, \boldsymbol{p}}(x)=\left|\operatorname{det} \boldsymbol{M}_{a, s}\right|^{\frac{-1}{2}} \varPhi\left(\boldsymbol{M}_{a, s}^{-1} x-\boldsymbol{p}\right) $

(2)

式中：a∈R⁺，s∈R，p∈R²，且R⁺、R、R²分别表示正实数集、实数集及二维实向量集；x为自变量；中间矩阵

$ \boldsymbol{M}_{a, s}=\left(\begin{array}{ll} a & \sqrt{a} s \\ 0 & \sqrt{a} \end{array}\right) $

(3)

式中M_{a, s}=S_sA_a，其中S_s为剪切矩阵，A_a为各向异性扩散矩阵，且有

$ \boldsymbol{S}_{s}=\left(\begin{array}{ll} 1 & s \\ 0 & 1 \end{array}\right) \quad \boldsymbol{A}_{a}=\left(\begin{array}{lc} a & 0 \\ 0 & \sqrt{a} \end{array}\right) $

(4)

Shearlet反变换为

$ f=\int_{\boldsymbol{R}^{2}} \int_{-\infty}^{\infty} \int_{0}^{\infty}\left\langle f, \varPhi_{a, s, \boldsymbol{p}}\right\rangle \varPhi_{a, s, \boldsymbol{p}} \frac{\mathrm{d} a}{a^{3}} \mathrm{~d} s \mathrm{~d} \boldsymbol{p} $

(5)

连续Shearlet变换可离散化。由于∀a>0，b∈A，可得对应的离散化尺度矩阵A和剪切矩阵B

$ \boldsymbol{A}=\left(\begin{array}{cc} a & 0 \\ 0 & \sqrt{a} \end{array}\right) \quad \boldsymbol{B}=\left(\begin{array}{ll} 1 & b \\ 0 & 1 \end{array}\right) $

(6)

当a=4，b=1时，$\boldsymbol{A}_{0}=\left(\begin{array}{ll} 4 & 0 \\ 0 & 2 \end{array}\right), \boldsymbol{B}_{0}=\left(\begin{array}{ll} 1 & 1 \\ 0 & 1 \end{array}\right)$

令a=2^－j，s=－l，则有

$ \varPhi_{j, l, \boldsymbol{k}}=\left|\operatorname{det} \boldsymbol{A}_{0}\right|^{\frac{j}{2}} \varPhi\left(\boldsymbol{B}_{0}^{l} \boldsymbol{A}_{0}^{j} \boldsymbol{x}-\boldsymbol{k}\right) $

(7)

式中：j, l∈Z，k∈Z²，离散化Shearlet变换为

$ {\rm{S}}{{\rm{H}}_\mathit{\Phi }}(f) = \left\langle {f,{\mathit{\Phi }_{j,l,\mathit{\boldsymbol{k}}}}} \right\rangle $

(8)

其反变换为

$ f = \left\langle {f,{\mathit{\Phi }_{j,l,\mathit{\boldsymbol{k}}}}} \right\rangle {\mathit{\Phi }_{j,l,\mathit{\boldsymbol{k}}}} $

(9)

1.2 利用Shearlet变换提高地震资料分辨率

地震数据s(t, x, y)经Shearlet变换，可表示为

$ {c_{j,l,\mathit{\boldsymbol{k}}}} = \left\langle {s,{\mathit{\Phi }_{j,l,\mathit{\boldsymbol{k}}}}} \right\rangle $

(10)

式中c_{j, l, k}是Shearlet变换后尺度为j、方向为l、位置为k的Shearlet域系数。

地震数据可从Shearlet域系数反变换得到

$ s = \sum\limits_{j,l,\mathit{\boldsymbol{k}}} {\left\langle {s,{\mathit{\Phi }_{j,l,\mathit{\boldsymbol{k}}}}} \right\rangle } {\mathit{\Phi }_{j,l,\mathit{\boldsymbol{k}}}} = \sum\limits_{j,l,\mathit{\boldsymbol{k}}} {{c_{j,l,\mathit{\boldsymbol{k}}}}} {\mathit{\Phi }_{j,l,\mathit{\boldsymbol{k}}}} $

(11)

对M×N个数据点的地震数据进行Shearlet变换，其变换后各尺度方向上的Shearlet系数都是M×N，Shearlet系数与原始地震数据在位置上具有一一对应的关系。因此，Shearlet变换对尺度方向表征更细腻，变换系数的物理意义更明确。

从式(11)可看到，对Shearlet域系数进行相应的处理，可得到新的地震数据。

地震波在地下传播过程中，若没有地层吸收和地震波频率及能量等的衰减，深、浅层反射波具有相同频宽的振幅谱，相位谱仅差一线性相位。如果把地震记录分成不同的频率(尺度)和不同的方向，所对应时间及空间的能量分布关系具有相似性，即对所有频率(尺度)而言，深层反射的能量与同一频率浅层反射能量之比应相同，所不同的只是各频率的绝对能量大小不一样。但由于地层的吸收及地震波的衰减，造成各频率能量对时间及空间的分布不同，在Shearlet域给表征地层的域系数乘以补偿因子，就起到了对地层吸收及地震波衰减的补偿作用。通过对Shearlet系数进行相应处理，可实现提高地震资料分辨率的目的。具体方法如下：

(1) 对地震数据进行Shearlet变换，得到Shearlet域系数。

(2) 因随机噪声本身无方向性且服从高斯分布，故随机噪声分解后，其Shearlet域系数值低且也是随机分布、无方向性，这在Shearlet域系数剖面中很易识别，应将该部分的域系数置为零。

(3) 分离完噪声与有效信号后，需对新的优势频带范围内的Shearlet域系数c′_{j, l, k}做能量补偿和提频处理。在Shearlet域系数中，最后一个Shearlet域系数数据代表原始信号的低频信息。为了提高中高频信号，不需对最后一个Shearlet域系数数据进行处理。地震数据的中低频信息能量相对较强，而中高频(尤其是高频)信息的能量相对较弱，但这些中高频信号含有更丰富的地层信息。通常，地震数据变换到Shearlet域，在Shearlet域中反映厚度较大地层的中低频信息域系数相对较强，反映厚度较薄地层信息的域系数较弱，需在Shearlet域对该部分信号进行补偿，使补偿后其能量与低、中频信号能量尽可能接近，从而拓宽地震信号的频宽，抬升中高频信号，进而提高地震资料的分辨率。

在Shearlet域补偿计算公式为

$ {c''_{j,l,\mathit{\boldsymbol{k}}}} = {c'_{j,l,\mathit{\boldsymbol{k}}}} \times \left( {\frac{{{\rm{MM}}}}{{(1 - \alpha )\left| {{{c'}_{j,l,\mathit{\boldsymbol{k}}}}} \right| + \alpha \times {\rm{MM}}}}} \right) $

(12)

式中：MM=max(|c′_{j, l, k}|)；α为提高分辨率因子，其取值视实情而定，通常的取值范围是0，1。α的取值主要考虑地震资料的信噪比和待分辨地层的厚度，若地震资料信噪比低，则α相对要大些；若地震资料信噪比较高，且要区分相对较薄地层，则α相对要小些。

(4) 对补偿后Shearlet域系数c″_{j, l, k}进行反变换，得到提高分辨率后的地震数据

$ {s_0} = \sum\limits_{j,l,\mathit{\boldsymbol{k}}} {{{c''}_{j,l,\mathit{\boldsymbol{k}}}}} {\mathit{\Phi }_{j,l,\mathit{\boldsymbol{k}}}} $

(13)

应用上述方法，既提高了地震资料分辨率，又保持了处理后信号的信噪比。

2 理论模型试验

为了验证此方法的可行性，设计了理论模型。测试重点是基于Shearlet变换提频方法的高频成分恢复能力。

构建如图 1所示的理论模型，横向范围是0~1000m，深度范围是0~300m。在模型100m处嵌入一套速度为2100m/s的砂岩体，砂岩体厚度范围是10~21m。在模型200m处嵌入两个叠置的速度为2500m/s的砂岩透镜体1和砂岩透镜体2。

图 2分别用主频为30Hz和45Hz的零相位Ricker子波进行正演，得到子波合成记录。通过对图 2的正演剖面做Shearlet变换提频处理以验证该方法的有效性。

将正演剖面进行Shearlet变换，其中Shearlet变换的尺度数为4，变换后形成49个域系数剖面。在该剖面中选取能在一定程度上反映砂体细节信息的Shearlet域系数，并用本文方法进行补偿(图 3)，可见在Shearlet域经本文方法补偿后(图 3b)的薄层信息体现得更清晰。

将补偿后的Shearlet域系数反变换得到提频后剖面(图 4a)。与图 2a相比，可见第一套砂体最薄处砂体边界完全分开：速度为2100m/s的砂体在第30~第52道之间地震同相轴已完全分开；速度为2500m/s的砂体在第32~第40道之间、第68~第87道之间地震同相轴也已完全分开；透镜体1和透镜体2的边界特征也与原模型匹配较好，两个透镜体细节信息也很好地展示出来。Shearlet变换提高分辨率方法效果(图 4a)与直接用45Hz主频子波正演的剖面(图 2b)基本一致。

另外，将此方法与小波变换提频进行对比，可见图 4b的第41~第51道之间速度为2500m/s的砂体地震同相轴未分开。应用小波变换后，两个透镜体附近的噪声干扰有所增强(图 4b)。

在理论模型中加入10%的随机噪声(图 5a)，验证该方法对随机噪声的稳定性。本文方法仍能分离出所有砂体的反射地震轴(图 5b)，并减少了噪声影响。而小波变换提频后噪声成分有所增强(图 6)。

为了比较Shearlet变换中不同分解尺度对提高地震资料分辨率的影响，分别以图 2的数据为试验对象，依次测试了尺度数为3、4、5的情形。图 7a和图 7c分别是尺度数为3和尺度数为5的处理结果，可见尺度为4是较合适(图 7b)。

从该模型试验可知，利用Shearlet变换提高分辨率处理方法能有效补偿有效地震波能量，有效提高薄砂体识别能力，提高地震资料分辨率。

3 实际资料应用

M区块前期已有地震数据无法满足现阶段的开发地质需求，因此将本文方法应用到该区块实际叠后地震数据，通过提频处理后再进行谱分解等属性提取及后续综合地质研究。

针对该区原始地震数据(图 8a)，应用本文方法提高分辨率后(图 8b)，地层波组特征更清晰，在组段内部，地质信息变得丰富，薄层及砂层组、单砂体呈现得更明显。从分辨率扩展前、后的频谱(图 9)可以看到，提频前的频带是8~70Hz，主频约35Hz，提频后的频带是8~90Hz，主频提高到约52Hz，提频后频带明显拓宽。

在图 10所示井旁合成地震记录中，地震子波是主频为45Hz的零相位Ricker子波。将井与原始地震数据及提高分辨率后地震数据分别进行比较，可见本文方法较好地提高了分辨率，且与井中信息吻合得较好。

通过对比、分析提高分辨率前、后地震数据的谱分解属性剖面(图 11)，可以看出本文方法显著地提高了三维地震资料的空间分辨率，属性剖面上河道(箭头处)特征更易于识别。

4 结束语

本文根据Shearlet变换的特点、地震资料中有效信号及优势频带与随机噪声在Shearlet域的差异，同时基于地层吸收及地震波衰减补偿原理，在Shearlet域对地层吸收及地震波衰减进行补偿。理论模型及实际地震数据的应用表明，利用Shearlet变换方法可以提高地震资料分辨率，并且此方法能保持较高信噪比。从实际资料的处理结果可看出，提高分辨率后剖面上地层和砂体特征变得更清晰，复波被有效分离，高频成分明显加强，低频部分也得到保持，在谱分解属性剖面上河道更清晰易辨，有利于后续开发地震及地质工作的开展。