用高分辨率谱分解方法研究三河-平谷大地震区的深部细结构

李 倩 酆少英 秦晶晶 田一鸣

(中国地震局地球物理勘探中心，郑州 450002)

0 引言

深地震反射技术是探测岩石圈精细结构及地球深部构造的有效技术手段之一(邓世广，2017)。然而深地震反射数据的探测深度较大，测线经过的区域构造复杂(如断裂带发育、 地形起伏变化大)，地震数据采集过程中的激发和接收条件均较差，最终导致深地震反射数据频率较低、 频带范围窄、 信噪比和分辨率低下及反射信号能量弱等(杨文采等，1999； Rossetal.，2004； 杨卓欣等，2006； 刘金凯等，2010； Kumaretal.，2011)。此外，由于地壳深部存在各种形态的复杂地质体，地震波在长距离传播中会发生能量吸收衰减等，导致深地震反射剖面中下地壳的反射波组经常存在反射能量较弱、 连续性较差、 呈交织状或带状分布的问题，这给地震资料的正确解释以及地质构造的认识带来了困难(Lietal.，1997； 徐明才等，1999； 王海燕等，2007； 李文辉等，2012)。

谱分解方法通常使用离散傅里叶变换，将地震数据从时域转换到时频域，利用不同频率的数据对不同地质情况的敏感性来描述不同的地质异常，是一种广泛应用于地震数据处理与解释的方法，常被用于进行储层预测、 定位复杂断层、 特殊地质特征分析以及油气识别等。近年来，随着深地震反射技术的发展，谱分解方法也逐渐成为一种提取深地震反射数据中弱信号的途径，并取得了一定的成果(Zhangetal.，2003； 邓攻等，2015； 曲中党等，2015； 邓世广，2017)。Rappin等(1991)分别使用加窗傅里叶变换和小波变换对ECORS的深反射地震数据进行了时频分析，以研究信号的频率成分，认为时频分析方法可以给出地壳内波的频率特性和下地壳结构的更多信息。Vasudevan等(2001)联合使用Wigner-Ville分布(WVD)技术和经验模态分解(EMD)技术研究了加拿大Alberta南部跨哈德逊造山带和阿尔伯塔基底横断面的深反射地震剖面，认为时频分析方法能够揭示在任何给定时间处随瞬时频率变化的能量横向变化规律； 同时也能揭示在任何瞬时频率处随时间变化的能量横向变化规律。邓攻等(2015)将S变换和Gabor变换应用于深地震反射数据中，经对比认为S变换比传统谱分解方法更具有优势，能够自动调节分辨率，有效提取深部弱信号，提高深部弱反射信号的分辨率和信噪比，并能够刻画低频细节特征，获得更好的成像效果。曲中党等(2015)使用基于S变换的软阈值滤波方法对深地震反射废炮的炮集数据进行处理，认为该方法能够有效压制频带范围内的混频干扰，提高反射波组的信噪比，增强反射波组的连续性，提高深部莫霍面反射波组的分辨率，为后续深地震反射剖面的精细处理奠定了良好基础。邓世广(2017)研究了基于OpenMP技术实现匹配追踪时频分析方法的并行计算方法，并将其应用于班公湖-怒江缝合带的深地震反射数据中，获得了深地震反射数据目标频段内的子波重构剖面与时频剖面。同时，其综合运用子波重构剖面、 时频剖面及纯波剖面对班公湖-怒江缝合带的深部反射结构特征进行了分析，并重点分析了上、 下地壳的反射特征，莫霍面内部反射特征以及亮点构造。尽管深地震反射数据中频率随时间的变化特征已经逐渐得到应用，但子波相位的变化特征尚未被利用，主要原因是深地震反射数据中相位的响应特征不明显，且难以提取。

随着高分辨率时频分析方法的出现，尝试利用新的时频分析方法对深地震反射数据进行时频分布特征分析及相位提取具有重要意义。稀疏约束反演谱分解(SCISD)是基于反演的谱分解方法，能够产生高分辨率的时频谱，其以地震信号时频谱的L1范数作为稀疏约束条件。该方法比基于变换的谱分解方法(如短时傅里叶变换、 连续小波变换及S变换等)具有更高的时频分辨率、 更好的时频聚焦性及抗噪性，能够提取出相应的时变子波相位信息，且不受时窗形状和长度的影响(韩利，2013； Lietal.，2016； 李倩等，2017)。稀疏约束反演谱分解方法已在石油地球物理勘探储层低频异常检测、 储层识别、 碳酸盐岩孔洞型储层识别等方面取得了一定成果。一方面，其能够提供高时间分辨率的结果，利于薄储层预测； 另一方面，能够提供子波相位变化特征以进行气藏检测及碳酸盐岩孔洞识别。然而，目前其尚未应用于深地震反射数据分析中。本文将稀疏约束谱分解方法应用于三河-平谷8.0级大地震区深地震反射数据中，分析深地震反射数据的时频、 相位特征，在目标频率范围内重构信号，提取深部低频有效弱信号，并将时频剖面和重构剖面相结合，研究三河-平谷8.0级大地震区深部反射细节特征。

1 基于非平稳褶积模型的反演谱分解算法

Bonar等(2010)提出了一种基于反演的复谱分解技术，该方法将传统的褶积模型拓展为一个多子波褶积模型，利用反褶积方法，使用不同频率零相位的雷克子波将地震记录在特定时刻分解成不同的频率分量，即

(1)

其中，s(t)是地震记录，k是与子波主频线性相关的参数，I表示参与计算的反射系数向量或子波的总个数，wk(t)是主频与角标k相对应的零相位雷克子波，rk(t)是具有地震记录中频率与相位信息的反射系数序列，n(t)指随机噪声。

将式(1)的褶积形式改写为乘积形式：

(2)

即

(3)

其中，W是复雷克子波库，r是包含了所有伪反射系数序列的矩阵(Portniaguineetal.，2004； Bonaretal.，2010； Hanetal.，2012； 韩利，2013； Lietal.，2016； 李倩等，2017)。

式(3)中矩阵r的元素个数远大于地震记录的元素个数，该方程是欠定的。本文采用L1范数进行正则化，将欠定问题转化成适定的，从而得到目标函数：

(4)

其中，‖r‖1是L1范数约束项，决定反射系数最终分布的形状。λ(λ>0)是正则化算子，即稀疏度参数，其控制解的稳定性。当使用稀疏约束条件(在本文中特指L1范数)求解方程(3)时，反演谱的分解方法称为稀疏约束反演谱分解方法。

对式(4)采用快速迭代软阈值算法(Becketal.，2009)求解，可以通过下式得到式(4)的解析解：

(5)

其中，α≥max eig(WHW)，即α的值必须大于等于WHW的最大特征值。yi是一个与ri和ri-1都有关的变量，用来提高计算效率。r的模即为地震记录的时频能量谱，而r的反正切为地震记录的时频相位谱。

2 深地震反射数据概况

三河-平谷8.0级大地震区位于北京东部华北沉降带与燕山隆起区相交会的平原区。在地质构造单元上位于华北裂陷盆地东北端与NW向张家口-蓬莱断裂带的复合、 交会部位，是新构造运动强烈的地区(刘保金等，2009)。该区破坏性地震发生较多，地震活动频繁，灾害程度较高(赵金仁等，2004)，区域地质构造见图1。

图1 研究区的地质构造和深地震反射剖面位置图

本文使用中国地震局地球物理勘探中心于2005年在北京地区布设的长约100km的深地震反射测线获得的反射地震数据开展研究。该数据采用爆炸震源激发，采样间隔为2ms，横向剖面桩号范围为6～99.6km(与刘保金等(2009)的工作一致)，测线穿过了三河-平谷8.0级大地震区以及北京地区的主要断裂构造(深地震反射测线位置见图1)。该区的深地震反射数据具有深反射数据的普遍特征，如深部有效信号频率较低、 中下地壳的地震波组能量较弱等。图2 中是三河-平谷大地震区深地震反射数据总的振幅谱，可以看出该实际数据的主要能量集中于10～35Hz频率范围，主频约为15Hz。本文将针对此数据进行时频分析与剖面重构研究。

图2 三河-平谷大地震区深地震反射数据的振幅谱

3 应用结果及分析

3.1 单道地震数据重构测试

地震信号在反演谱域是稀疏的，可以用高能量的稀疏成分代表几乎整个地震信号的能量，以实现地震信号的重构。为了测试稀疏约束反演谱分解方法的重构性，首先截取三河-平谷大地震区深地震反射数据中某道9.5～11.5s的深部数据单道进行测试(图3)。图3a 是原始单道记录，其地震数据波形特征明显，振幅强弱随时间变化，强振幅值分别出现在10.05s、 10.6s及10.95s处； 从图3b 所示的稀疏约束反演谱分解的时频能量谱中可以看出，该道地震数据频率范围主要集中在10～27Hz，主频约为16Hz，信号的频率随时间的变化特征较为明显，时频能量谱的能量聚焦性好，时间分辨率和频率分辨率同样较高，强能量团与图3a 中的强振幅值一一对应，能够有效识别出地震信号中的异常幅值； 图3c 是稀疏约束反演谱分解的时频相位谱，可以看出地震信号的相位信息丰富，相位特征变化明显，某些微小信号由于振幅能量较弱在时频能量谱上被忽略，在相位谱上也显现出相位大小及正负变化，使相位谱具有了实际意义，有利于地震信号的解释； 图3d 是原始单道(蓝色)和重构数据(红色)的对比，可以看出，通过稀疏约束反演谱分解方法可以较为准确地重构地震信号，信号的幅值没有得到太大的衰减。由此可见，稀疏约束反演谱分解方法在深地震反射数据中的应用具有一定的可行性。

图3 三河-平谷大地震区深地震反射单道数据的时频分析与重构

3.2 深地震反射剖面的重构与时频分析

将稀疏约束反演谱分解应用于整个深地震反射剖面(图4)，通过筛选合适的频率范围，得到重构剖面(图5)及时频相位谱(图6)。基本流程如下： 首先，对整个深地震反射数据进行稀疏约束反演谱分解，得到二维叠加剖面的时频分布特征； 其次，根据地震信号频谱分布范围，可知有效反射信号的能量主要集中在10～35Hz，因此选择10～35Hz的子波对剖面进行重构，得到10～35Hz重构剖面(图5)，并选取主频为15Hz的分频剖面(图6)进行地壳结构特征分析。

图4 三河-平谷大地震区的深地震反射时间叠加剖面

图5 三河-平谷大地震区的深地震反射数据重构剖面(10～35Hz)

图6 三河-平谷大地震区的深地震反射数据时频能量谱(15Hz)

对比分析原始叠加剖面(图4)、 重构剖面(图5)及15Hz时频剖面(图6)，尤其是在重构剖面(图5)及15Hz时频剖面(图6)上可以看出，反射同相轴能量在上地壳浅部(双程旅行时5s以浅)较强，且分段连续，同相轴中断、 位错明显，显示存在断裂，断裂密布，特征明显； 双程旅行时5.8～7.5s之间存在一条明显的反射条带，以近水平的叠层状反射为主，是上、 下地壳之间的分界面； 下地壳反射能量较弱，仅存在少量倾斜反射，可能由于下地壳为塑性状态，难以产生有效反射； 壳幔过渡带(莫霍过渡带)位于双程旅行时10～11.3s之间，在剖面上自双程旅行时约10s时开始出现，纵向上持续1～1.3s，表现为叠层状、 纵向上有一定延续时间的强反射条带； 壳幔过渡带之下的上地幔顶部表现为弱小、 凌乱的不连续反射，仅在剖面NW端约11.6s处出现了弧状强反射同相轴，延续至剖面桩号约75km处。对比图4—6可以看出，10～35Hz重构剖面及15Hz时频剖面比原始剖面具有更高的信噪比和分辨率，提高了同相轴的连续性，原叠加剖面上被噪声湮灭的低频细节特征也能够被清晰地刻画； 且其将整个地震剖面中值得关注的某一频段信号进行抽离并单独显示，有助于排除其他干扰分量、 深入分析目标频段的信号特征、 提高对地震剖面的整体认识。此外，壳幔过渡带反射波组的横向变化、 纵向的分布范围及内部的反射波组结构特征也比原始地震剖面上的更加清晰。

3.3 剖面的浅部构造及断层识别

振幅谱和相位谱共同组成了一个地震信号的完整频谱特征。然而，目前的地震资料解释工作一般使用振幅谱，由于相位谱复杂难懂，经常被忽视。反射地震信号的相位信息一般可用于确定地下地质构造的边界情况，但由于其分辨率不足且较为复杂，较少有学者利用信号的相位谱识别地下构造及断层信息。稀疏约束反演谱分解方法能够得到较高分辨率的相位谱，为我们利用相位谱进行构造及断层的判别提供了帮助，且可以与地震反射时间剖面上的反射波组特征相互验证，确定构造及断裂的相关信息。

为了验证时频相位谱是否能够有助于识别构造和断裂，选取三河-平谷大地震区深地震反射数据中浅部前4s构造及断层丰富的部分(横向剖面桩号范围为11～46km，截取时间范围为0～4s，如图7a 所示)，针对此部分数据进行稀疏约束反演谱分解，提取对应的时频相位谱剖面进行分析。

从图7a(红色实线为根据波组特征解释的断层)中可以看出，剖面反射能量较强，反射层位丰富，存在几组自西向东倾斜的密集反射层，具有典型的沉积盆地特征，断裂特征较为明显，然而仍不够清晰明了，容易漏掉有效信息，且无法判断上断点的位置。对该剖面数据进行稀疏约束反演谱分解处理，结合图2 所示的深地震反射数据的频率范围，分别提取主频16Hz、 小于主频10Hz和大于主频26Hz的时频相位谱，如图7b—d所示。与图7a 中的原始地震剖面进行比较可以看出，断层以及沉积盆地特征在时频相位谱中更加明显，与地震剖面相比也更容易识别。时频相位谱的相位变化明显，构造及断裂特征的时频能量分布在相位谱图上较为清晰，信号细节部分的变化也得到了体现，剖面横向及纵向上的相位特征变化清晰明了，断裂向上、 向下的延伸情况一目了然，有助于确定上断点的位置，也有利于判断构造边界。对比不同频率的时频相位谱(图7b—d)可以看出，10Hz的时频相位谱由于频率较低，缺失了一些高频信息，浅部有些位置的细节变化无法得到体现； 16Hz的时频相位谱由于在主频附近，信息较为全面且丰富，能量分布稍显均衡，相位变化的细节特征体现明显，构造和断裂信息较为丰富，需要仔细判断； 26Hz的时频相位谱上冗余信息较少，主要反射层特征明显，能够很容易分辨出构造边界、 断裂上断点位置及断距。对该实际地震数据的稀疏约束反演谱分解结果进行分析认为，该方法提供的时频相位谱简单易懂，为地震资料解释人员从相位的角度分析地下的地质构造信息提供了有效帮助，可对不同频率的时频相位谱与地震反射剖面进行综合分析，有助于识别地质构造及断裂。

3.4 剖面深部莫霍过渡带的特征识别

在地球物理学方面，学者们认为莫霍面通常是一个速度突变面，具有一定厚度，是由一组高速和低速薄层组成的速度梯级变化带，为地壳和地幔的分界面。莫霍面内部由一系列近水平、 相互结合平整的反射叠层组成，具有壳幔过渡带的性质(曾融生等，1985； 王椿镛等，1993； 张先康等，1996)。

在原始叠加剖面、 重构剖面及时频剖面(图4—6)中能够看到在下地壳及上地幔之间存在能量较强的反射波组(位于双程旅行时10～11.3s之间)，特别是在重构剖面以及时频剖面上得到了充分的体现，认为此处强反射波组是莫霍过渡带的显示。为了进一步分析研究莫霍过渡带内部的反射特征，利用10～35Hz的子波重构剖面及不同频率的单频相位谱突出刻画了9.5～11.5s范围内的反射信息(图8)。从本次处理得到的重构剖面及不同频率时频相位谱图中可以看出，莫霍过渡带表现为由彼此近平行、 相互结合平整的反射段组成的叠层反射结构，特别是在重构剖面和时频相位谱图上得到了突出的体现。该过渡带在重构剖面上约TWT 10s时开始出现，纵向上的持续时间约为1～1.2s，底界对应于莫霍面的位置。采用地壳平均P波速度为6.2km/s(赵金仁等，1999，2004； 张先康等，2002； 刘保金等，2011)计算莫霍过渡带的厚度及深度，根据双程旅行时可计算得到莫霍面的深度为34～35.5km，厚3.1～3.7km。重构剖面所揭示的莫霍过渡带整体较为平缓，横向上可连续追踪，莫霍过渡带的反射能量在剖面桩号约36km附近明显变弱。从时频相位谱图上可以看出，在莫霍过渡带内部，反射震相相位横向上变化明显，说明壳幔过渡带存在着强烈的横向不均匀性。剖面两端反射同相轴波组较为密集，多呈近水平状，并伴有错断、 相交以及弱能量等现象； 在剖面桩号约36km附近相位变化剧烈、 细碎且杂乱，使剖面NW和SE部分的反射特征被有效区分； 在剖面北西侧，反射波组呈近水平展布，在剖面东南侧，反射波组明显向NW倾斜，且与北西侧的反射波组存在一定的落差，暗示了该处可能有一条切穿莫霍面的深大断裂。此外，对比不同频率时频相位谱可以发现，低频的10Hz和15Hz时频相位谱能量较强，所含信息较为丰富，横向上相位变化明显，莫霍过渡带内部的变化显而易见； 而在25Hz的时频相位谱上只留下了一些能量较强的高频信息，缺失了相关的低频信息，细节变化无法得到体现。

图8 三河-平谷大地震区深地震反射莫霍面范围的局部剖面、 重构剖面及不同频率的时频相位谱剖面

4 结论

本文介绍了一种目前在油气储层预测、 碳酸盐岩孔洞识别等方面取得了一定成果的稀疏约束反演谱分解方法，并将其应用于深地震反射数据中，对三河-平谷8.0级大地震区深地震反射数据的单道及叠加剖面进行了时频特征分析以及数据重构，发现重构剖面及时频剖面比原始剖面具有更高的信噪比和分辨率，提高了剖面同相轴的连续性，不仅能标定深地震反射数据中弱信号不同时刻的频率分量，还能清晰刻画原叠加剖面上被噪声湮灭的低频细节特征； 高分辨率的时频相位谱具有丰富的相位信息，相位特征变化明显，可借助时频相位谱进行断层识别，有助于判断构造边界； 子波重构技术及高分辨率的单频相位谱有助于刻画出莫霍过渡带内部的反射结构特征。以上结论表明，稀疏约束反演谱分解方法适合用于深地震反射数据的谱分解。

对原始反射地震剖面进行时频特征分析，并通过地震信号的可重构性对深地震反射数据进行有目的的重构，能够提取深部有效弱信号； 对时频剖面和重构剖面的联合对比分析有助于深入研究深部精细构造，建立对深部构造的直观认识。该方法以新的视角分析深地震反射剖面，有助于取得新的地质认识，具有十分重要的研究意义，可以在深地震反射数据解释中进行推广应用。

致谢中国地震局地球物理勘探中心为本研究提供了数据支持； 刘保金研究员在地震资料解释方面给予了指导。在此一并表示感谢！