2. 海洋石油勘探国家工程实验室, 陕西西安 710049;
3. 酒泉卫星发射中心, 甘肃酒泉 732750
2. National Engineering Laboratory for Offshore Oil Exploration, Xi'an 710049, China;
3. Jiuquan Satellite Launch Center, Jiuquan 732750, China
对地震信号进行稀疏表示可以有效揭示地震信号的本质特征, 有利于更为清晰和直观地认识和描述地震信号。在数字信号处理中, 常常将信号在给定的函数(或矢量)集上进行投影表示, 把这种在变换域上表达原始信号的方式称为信号表示或者信号变换。信号的稀疏表示就是在变换域上寻找尽量少的基函数重构逼近原始信号, 从而获得原始信号简洁而有效的表达形式。目前常用的变换方法主要有傅里叶变换、离散余弦变换(Discrete Cosine Transform, DCT)、小波变换、曲波(Curvelet)变换等, 以及联合这些变换方法组成超完备字典的信号分解方法。而稀疏性是压缩感知的前提, 如果采集很少一部分数据并且希望从这些少量数据中解压出大量信息, 就需要保证:①这些采集到的少量数据包含了原始信号的全局信息; ②存在一种算法能够从这些少量的数据中还原出原先的信息。因此, 信号的稀疏表示与压缩感知息息相关。
以上基于特定数学变换的信号表示方法, 根据各自变换基的特性, 可对相应结构模式的信号进行有效表示。然而, 考虑到实际地震数据的复杂性, 如叠前炮集中的体波、面波、工业干扰、抽油机振动噪声、滑动扫描接收方式中的谐波噪声等, 其波形结构差异明显。另外, 叠后地震数据中稳定沉积地层及与陆相储层相关的地质体(如点坝、河道等)的地震响应信号波形结构也存在明显差异。这样, 在分析实际地震信号时, 仅选择单一的变换方法, 往往很难对地震信号实现稀疏分解。如果将以上多种信号表示方法联合起来组成超完备冗余表示字典, 就有可能实现对复杂信号的稀疏表示。基于近年来信号稀疏表示方面研究取得的进展, 相关学者发展了形态成分分析(Morphological Component Analysis, MCA)的信号分析方法[1-2], 就是利用超完备的复合波形字典对信号进行稀疏分解。这种理论认为复杂信号可以表示成若干信号成分的线性组合, 假设每种成分都能找到对其进行稀疏表示的波形字典, 而该字典对其它成分不能稀疏表示, 那么复杂信号就能够由这些波形字典联合组成的超完备冗余字典进行稀疏表示。另外, 如果能够根据地震数据自身的特点, 自适应地构造出适用于这种地震数据的变换基, 也可以实现对地震数据的稀疏分解。OLSHAUSEN等[3]从待分析信号自身出发, 采用机器学习的方法自适应地构造出超完备字典, 带来信号表示的稀疏化。RUBINSTEIN等[4]结合数学变换, 给出基于稀疏字典学习的信号稀疏表示方法, 该方法提高了自适应波形字典的普适性, 降低了其对训练样本数量的要求, 并且使自适应字典具有快速的前向和伴随算子, 从而使得自适应字典方法适用于高维字典的构造和高维信号的分析。因此, 将复合型和学习型超完备波形字典的信号表示方法用于地震信号处理和解释, 有利于更好地理解地震信号包含的信息, 是非常有意义的发展方向。
本文根据地震数据中波形形态特征的差异, 构造了两种基于稀疏表示的超完备字典, 根据其功能分为:体波与面波干扰分离以及强反射屏蔽效应的消除, 并且利用实际数据验证了方法的有效性。
1 基于超完备字典的地震信号表示理论本文根据地震信号有效波和干扰波的波形结构特点, 构造了用于压制面波、工业噪声、谐波噪声及抽油机噪声等干扰的一维炮集数据信噪分离模型。对于炮集数据, 每道信号可看作有效信号和噪声两种分量的线性叠加, 地震记录s可表示为:
| $ \mathit{\boldsymbol{s = }}{\mathit{\boldsymbol{s}}_1} + {\mathit{\boldsymbol{s}}_2} + \mathit{\boldsymbol{n}} $ | (1) |
式中:s1为有效信号分量; s2为噪声分量; n为随机噪声分量。
(1) 式引入了加性随机噪声n, 从而提高了信号模型的普适性, 使其不仅能够描述受到噪声干扰的地震数据, 也允许实际数据包含适当的测量误差或其它杂波干扰。假设随机噪声分量n服从零均值的高斯白噪分布, 其在每个位置的取值都是一个N(0, σ2)随机过程。
(1) 式表示在噪声干扰下的一种多源混合信号恢复重建和分离问题, 这里的两种源信号分别指有效波信号s1和噪声s2。MCA的目标便是分别提取出s1, s2两种成分。根据s1和s2具有不同的形态特征, MCA理论假设s1和s2可以分别由字典Φ1和Φ2有效地稀疏表示, 但是用Φ2稀疏表示s1和用Φ1稀疏表示s2时稀疏性差, 那么有效波信号和干扰噪声可以在Φ1和Φ2联合的超完备波形字典中得到稀疏表示, 从而通过求解如下的最优化问题实现两种信号分量的恢复和分离:
| $ \begin{array}{l} \mathop {{\rm{argmin}}}\limits_{\{ \alpha 1,\alpha 2\} } \left\| {{\alpha _1}} \right\|{\;_1} + \left\| {{\alpha _2}} \right\|{\;_1}\\ {\rm{s}}.{\rm{t}}.\left\| {\mathit{\boldsymbol{s - }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_1}{\mathit{\boldsymbol{\alpha }}_1}\mathit{\boldsymbol{ - }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_2}{\mathit{\boldsymbol{\alpha }}_2}} \right\|\;_2^2 \le {\rm{ \mathsf{ ε} }} \end{array} $ | (2) |
式中:α1为有效波信号在字典Φ1中的向量表示形式; α2为噪声在字典Φ2中的向量表示形式; ε为信号重建的误差门限, 使得该优化问题在面对噪声干扰及数据缺失时仍然能够取得合理的稀疏解。
BRUCE等[5]提出了分块坐标松弛算法(Block-Coordinate-Relaxation, BCR)对(2)式进行求解, 其核心算法为通过每次迭代轮流对系数向量α1和α2进行更新。将(2)式中的优化问题表示为:
| $ \begin{array}{l} f({\mathit{\boldsymbol{\alpha }}_{{\rm{all}}}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{s - }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_1}{\mathit{\boldsymbol{\alpha }}_1}\mathit{\boldsymbol{ - }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_2}{\mathit{\boldsymbol{\alpha }}_2}} \right\|\;_2^2 + \lambda (\left\| {{\mathit{\boldsymbol{\alpha }}_1}} \right\|{\;_1} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{\alpha }}_2}} \right\|{\;_1}) = f({\mathit{\boldsymbol{\alpha }}_1},{\mathit{\boldsymbol{\alpha }}_2}) \end{array} $ | (3) |
其中,
假设k次迭代后系数向量更新为α1k, α2k, 进行k+1次迭代时, 首先更新α1而保持α2k不变, 即:
| $ f({\mathit{\boldsymbol{\alpha }}_1},\mathit{\boldsymbol{\alpha }}_2^k) = \frac{1}{2}\left\| {\mathit{\boldsymbol{s - }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_1}{\mathit{\boldsymbol{\alpha }}_1}\mathit{\boldsymbol{ - }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_2}\mathit{\boldsymbol{\alpha }}_2^k} \right\|\;_2^2 + \lambda \left\| {{\mathit{\boldsymbol{\alpha }}_1}} \right\|{\;_1} $ | (4) |
则α1k+1具有硬阈值函数解:
| $ \mathit{\boldsymbol{\alpha }}_1^{k + 1} = {F_\lambda }[\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_1^{\rm{T}}(\mathit{\boldsymbol{s - }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_2}\mathit{\boldsymbol{\alpha }}_2^k)] $ | (5) |
式中:Fλ为硬阈值算子。类似地, 下一次迭代时, 先更新α2, 保持α1k+1不变, 求解α2k+1的硬阈值函数解:
| $ \mathit{\boldsymbol{\alpha }}_2^{k{\rm{ + }}1} = {F_\lambda }[\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_2^{\rm{T}}(\mathit{\boldsymbol{s - }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_1}\mathit{\boldsymbol{\alpha }}_1^{k + 1})] $ | (6) |
如此反复迭代, 求解两部分的稀疏系数, 最终从混合信号s中将具有不同波形特征的信号s1和s2分量进行分离。
另外, 根据叠后地震资料中稳定沉积地层及与储层相关的地质体(河道、点坝等)对应的地震响应波形结构特点, 我们同时构造了用于提取高精度陆相沉积地震响应信号的叠后地震资料模型。其数学表达形式与(1)式类似。
2 基于Q可调小波的面波干扰压制面波是陆上地震记录中一种具有明显的高振幅和频散特性的规则干扰波, 具有能量强、频率低、视速度低等特点, 在单炮地震记录中一般呈扇形分布, 遮盖了绝大部分的体波信息。合理地压制面波干扰, 有利于提高地震资料的信噪比。常用的压制面波干扰的方法包括带通滤波、高通滤波和多通道滤波以及基于视速度特征差异的F-K滤波、沿时间-偏移距的加窗F-K滤波等[6-9]。基于体波及面波在变换域分布区域的不同, 小波变换、Curvelet变换、Radon变换等方法也被成功用于面波压制。
考虑到实际地震资料非常复杂, 仅选择单一的变换方法往往很难对地震信号实现稀疏表示。如果把以上多种信号表示方法联合起来组成超完备冗余表示字典, 就有可能实现对复杂信号的稀疏表示。陈文超等[10]提出了联合连续小波变换与局部余弦变换的超完备字典, 成功实现了面波及体波分量的分离。TRAD等[11-12]提出将线性Radon变换和双曲Radon变换联合构成一种混合Radon变换, 该方法可对符合模型要求的地震资料取得较好的波场分离效果。YARHAM等[13]将两种波场分量的初步估计用作Curvelet域中待求解真实波场分量的加权, 在Curvelet域建立稀疏优化分离模型, 得到较好的波场分离结果。WANG等[14]提出将一维平稳小波变换与局部离散余弦变换联合构成超完备波形字典, 初步探索了两种波场分量的自适应稀疏分解途径。
本文根据体波和面波信号的波形特征差异, 利用Q可调小波函数[15], 分别构造了能够稀疏表示体波和面波分量的Q小波函数。针对低震荡体波选择低Q值小波变换作为其稀疏表示字典Φ1, 针对高震荡面波选择高Q值小波变换作为其稀疏表示字典Φ2, 建立基于双波形字典的优化问题, 如公式(2)所示, 求解该稀疏优化问题, 按照公式(5)和公式(6)反复迭代求解, 得到最终的信噪分离结果。
图 1a和图 1b分别为能稀疏表示体波和面波的Q小波函数。由图 1可以看出, Q比较小的母小波函数具有宽带特性, 能稀疏表示频带较宽的体波信号; 而Q较大的小波函数能稀疏表示频带较窄的面波干扰。图 2为我国东部某油田实际炮集记录分离结果。其中, 图 2a为原始炮集记录, 该记录受到了强面波噪声的干扰; 图 2b为采用本文方法分离出的体波信号, 可以看出, 面波能量得到明显压制, 中深部目的层段同相轴连续性较好; 图 2c为分离出的面波干扰噪声剖面, 可以看出, 面波噪声中几乎观察不到有效信号残余。也就是说, 本文方法在明显压制面波干扰的同时, 对有效信号有良好的保真性能。
|
图 1 不同Q值的小波母函数 a Q为1的小波母函数; b Q为4的小波母函数 |
|
图 2 我国东部某油田实际炮集记录分离结果 a 原始炮集记录; b 本文方法分离出的体波信号; c 本文方法分离出的面波干扰 |
河道砂体的预测与识别是陆相储层预测面临的主要问题之一。近年来, 地震资料处理、解释方法在砂体预测与识别方面取得了显著的效果。例如, 相干体算法能清晰地刻画河道、断层等地质体边界; 曲率属性可以帮助解释人员识别沉积地层形态特征, 研究表明, 曲率属性对于描述断层、裂缝走向及分布十分有效; 频谱成像方法利用地层调谐效应可以增强对薄砂体的识别与预测能力。
上述地震处理、解释方法应用成功的前提条件是储层结构相对简单, 而且需要高信噪比地震数据作为基础[16]。然而实际勘探开发所面临的目标复杂。以大庆油田的扶余油层为例, 扶余油层与上部青山口组泥岩存在较明显的波阻抗差, 因而扶余油层顶面T2反射同相轴表现为连续的强反射特征。在T2强反射背景影响下, 扶余油藏的F11油层组砂体地震响应识别困难, 储层预测成功率低[17]。另外, 煤层的强反射特征也是影响煤系储层识别的主要因素[18]。目前多子波分解技术、匹配追踪方法及EMD分解等方法被用于压制强地震反射背景, 提高了储层识别率。
本文根据河道砂体等与储层相关的地质体空间展布范围较小的特点, 考虑到强背景反射连续性好、稳定等特征, 构造了能分别稀疏表示河道砂体等与储层相关的地质体及稳定地层的地震响应的超完备波形字典Φ1与Φ2, 建立稀疏优化问题, 如公式(2)所示, 通过公式(5)、公式(6)迭代求解, 实现了强背景干扰分离, 突出了砂体地震响应信号。
图 3a为能稀疏表示稳定沉积地层地震响应的Curvelet变换原子, 图 3b为能稀疏表示河道砂体等与储层相关的地质体地震响应的二维平稳小波变换原子。由图 3可以看出, Curvelet变换原子适宜于具有各向异性的曲线状结构信号的多尺度和多方向分析; 二维平稳小波变换原子适宜点状结构信号的多尺度分析。
|
图 3 Curvelet变换(a)与二维平稳小波变换(b)原子对比 |
图 4是我国东部某油田实际叠后地震剖面, 其中黑、白线之间区域为目标层。可以看出, 目标层内砂体地震响应完全被上覆泥岩强反射所掩盖, 根据原始数据很难刻画储层位置。图 5a为图 4中沿白线所示时间位置的沿层切片; 图 5b为用本文方法分离出的河道砂体等与储层相关的地震响应信号沿层切片。对比图 5a和图 5b可以看出, 经过强屏蔽层地震响应的压制, 河道砂体等地震响应信号明显凸显, 特别是图中黑色椭圆所示位置, 河道空间形态清晰。图 5c为本文方法分离出的屏蔽层地震响应信号沿层切片, 可以看出屏蔽层的空间大致变化情况。
|
图 4 我国东部某油田实际地震剖面 |
|
图 5 某油田压制强反射背景后的砂体地震响应切片 a 图 4所示地震剖面沿白线时间切片; b 本文方法分离出的河道砂体等与储层相关的地震响应信号沿层切片; c 本文方法分离出的屏蔽层地震响应信号沿层切片 |
本文根据地震数据中特征波形形态差异, 构造能分别稀疏表示数据组成成分的超完备字典, 建立了基于超完备字典的地震信号表示方法, 实现基于Q可调小波的面波干扰压制和基于强反射背景压制的砂体地震响应信号提取技术。通过实际数据分析, 检验了方法的有效性。与目前常用的基于单变换信号表示方法相比, 基于超完备字典的方法, 充分考虑了地震资料的复杂性, 结合了不同变换对地震数据中波形特征差异明显的多种成分的稀疏表示性能, 明显增强了在变换域中解决实际问题的能力。
致谢: 感谢大庆油田勘探开发研究院的陈树民副院长、王建民副总工、王成和赵忠华专家, 中国科学院地质与地球物理研究所李幼铭研究员及海洋石油勘探国家工程实验室副主任高静怀教授对本文工作的指导与支持。| [1] | STARCK J L, ELAD M, DONOHO D L. Redundant multiscale transforms and their application for morphological component analysis[J]. Advances in Imaging and Electron Physics, 2004, 132(4): 287-348 |
| [2] | STARCK J L, ELAD M, DONOHO D L. Image decomposition via the combination of sparse representations and a variational approach[J]. IEEE Transactions on Image Processing, 2005, 14(10): 1570-1582 DOI:10.1109/TIP.2005.852206 |
| [3] | OLSHAUSEN B A, FIELD D J. Emergence of simple-cell receptive field properties by learning a sparse code for natural images[J]. Nature, 1996, 381(6583): 607-609 DOI:10.1038/381607a0 |
| [4] | RUBINSTEIN R, ZIBULEVSKY M, ELAD M. Double sparsity:learning sparse dictionaries for sparse signal approximation[J]. IEEE Transactions on Signal Processing, 2010, 58(3): 1553-1564 DOI:10.1109/TSP.2009.2036477 |
| [5] | BRUCE A G, SARDY S, TSENG P. Block coordinate relaxation methods for nonparamatric signal denoising[C]//Proceedings of International Society for Optical Engineering. USA: Society of Photo-optical Instrumentation Engineers, 1998: 75-86 |
| [6] | EMBREE P, BURG J P, BACKUS M M. Wide-band velocity filtering-the pie-slice process[J]. Geophysics, 1963, 28(6): 948-974 DOI:10.1190/1.1439310 |
| [7] | GELISLI K, KARSLI H. F-K filtering using the Hartley transform[J]. Journal of Seismic Exploration, 1988, 7(2): 101-108 |
| [8] | TREITEL S, SHANKS J L, FRASIER C W. Some aspects of fan filtering[J]. Geophysics, 1967, 32(5): 789-800 DOI:10.1190/1.1439889 |
| [9] |
刘法启, 张关泉. 小波变换与F-K算法在滤波中的应用[J].
石油地球物理勘探, 1996, 31(6): 782-791 LIU F Q, ZHANG G Q. Application of wavelet transform and F-K algorithm in filtering[J]. Oil Geophysical Prospecting, 1996, 31(6): 782-791 |
| [10] |
陈文超, 王伟, 高静怀. 基于地震信号波形形态差异的面波噪声稀疏优化分离方法[J].
地球物理学报, 2013, 56(8): 2771-2782 CHEN W C, WANG W, GAO J H. Sparsity optimized separation of Ground-roll noise based on morphological diversity of seismic waveform components[J]. Chinese Journal of Geophysics, 2013, 56(8): 2771-2782 DOI:10.6038/cjg20130825 |
| [11] | TRAD D, SACCHI M, ULRYCH T J. A hybrid linear-hyperbolic Radon transform[J]. Journal of Seismic Exploration, 2001, 9(4): 303-318 |
| [12] | TRAD D, ULRYCH T J, SACCHI M. Latest views of the sparse radon transform[J]. Geophysics, 2003, 68(2): 386-399 |
| [13] | YARHAM C, BOENIGER U, HERRMANN F J. Curvelet-based ground roll removal[J]. Expanded Abstracts of 76th Annual Internat SEG Mtg, 2006: 2777-2781 |
| [14] | WANG W, GAO J H, CHEN W C, et al. Data adaptive ground-roll attenuation via sparsity promotion[J]. Journal of Applied Geophysics, 2012, 83(1): 19-28 |
| [15] | SELESNICK I W. Wavelet transform with tunable Q-factor[J]. IEEE Transactions on Signal Processing, 2011, 59(8): 3560-3575 DOI:10.1109/TSP.2011.2143711 |
| [16] | DE MATOS M C, YENUGU M, ANGELO S M, et al. Integrated seismic texture segmentation and cluster analysis applied to channel delineation and chert reservoir characterization[J]. Geophysics, 2011, 76(5): 11-21 DOI:10.1190/geo2010-0150.1 |
| [17] | GUAN X W, XIE C L, FENG X Y, et al. Application of wavelet decomposition technology in thin sand identification under strong reflection[C]. SPG/SEG Beijing 2016 International Geophysical Conference Technical Program Expanded Abstracts. Beijing: SPG & SEG, 2016: 417-420 |
| [18] | ZHANG Z J, ZHANG J H, GONG X N, et al. The method research of stripping strong reflection of coal seams[J]. Expanded Abstracts of 77th EAGE Conference and Exhibition, 2015: 1-5 |