匹配滤波和曲波去噪相结合的气枪主动源弱信号提取
2020-09-26谭俊卿杨润海向涯王彬姜金钟
谭俊卿 杨润海 向涯 王彬 姜金钟
摘要：提出由一维模板匹配滤波技术（MFT）和二维曲波（Curvelet）变换法穿插的数据处理方法，即先通过一维模板匹配滤波方法得到相关系数，将相关系数组成二维数据并用曲波变换法处理，最后将各道对应相关系数分别与模板信号褶积，得到高信噪比恢复信号。将该方法用于处理仿真数据和宾川主动源气枪信号，结果表明：在宾川主动源气枪信号的处理中，本方法较单一数据去噪方法恢复能力更好，可对气枪主动源的模拟与实际信号进行更好的数据提取与恢复，得到噪声干扰更少的地下介质波速变化，提高低信噪比数据的利用率与可分析性。
关键词：主动源气枪信号;模板匹配滤波技术;曲波变换;波速变化
中图分类号：P315.61文献标识码：A文章编号：1000-0666（2020）04-0701-10
0引言
地震波是照亮地球内部的一盏明灯，是少数能穿透整个地球的信号之一，携带了丰富的地下介质结构和物性信息，是研究地球内部结构的最有效工具（陈颙等，2007a，b，c）。地下介质波速变化可以通过高度相似的地震信号计算得到，为此，国内外已开展了许多主动源探测项目，其中，气枪震源具有绿色环保、信号低频、能量转换效率高等优点，在激发条件不变时信号重复性高、激发时刻精确、激发位置可控，为精细获取区域地下介质结构变化提供数据基础。在陈颙院士“地下明灯计划”的推动下，大容量气枪震源被移植到陆地水库中，进行地下介质变化的观测与研究（Chen et al，2008;王彬等，2012;王宝善等，2016;王伟涛等，2017;丘学林等，2007;向涯等，2017;杨微等，2016;林建民等，2010）。
地下介质在地震发生前都会有裂隙和孔隙压力的变化，这个变化在介质波速上有所体现，这种波速的细微变化能被主动源信号所捕获，有助于认识地震孕育的物理过程（黄亦磊等，2017，Simmons，1964;Yukutanke，1988）。地震记录中与地下介质变化相关的信息较微弱，如何从强背景噪声中识别与提取地震弱信号，是区域结构探测和地下介质物性结构变化研究的关键（刘必灯等，2011;徐逸鶴等，2016;翟秋实等，2016;栾奕等，2016;刘自凤等，2015;魏芸芸等，2016;姚佳琪等，2017;Wang et al，2010）。
运用叠加、滤波等方法压制噪声，是提高地震资料信噪比和增强弱信号提取能力的主要途径。信号叠加可以有效识别淹没在强背景噪声下的气枪信号，但长时间的叠加受人类不定期活动、固体潮变化、大气压等因素影响，叠加信号精度一般。滤波方法有很多种，其中频域滤波通过Fourier变换得到原始信号频谱，在分析非稳态的地震信号时具有一定局限，其滤波方式简单，但效果一般;小波变换弥补了Fourier变换的不足，在频域和时域都能较好地表达局部化特征，在表示“线奇异”和“面奇异”目标时，不能获得最优的非线性逼近;曲波变换是继小波变换、脊波变换等之后发展起来的一种数据稀疏表示方法，具有多分辨率、时频局部性、多方向性等特点，克服了小波变换的不足，在图像和指纹识别等领域应用广泛（Candes et al，2004;袁艳华等，2010;林春，王绪本，2009）。模板匹配滤波技术（以下简称匹配滤波）是一种基于互相关的信号识别方法，该方法在滑动窗互相关（Sliding-Window Cross-Correlation，简称SCC）检测技术的基础上发展而来，是低信噪比数据中识别、提取弱信号的一种方法，在医学、电子通信、图像识别等众多领域应用较广（李璐等，2017;Yang et al，2009;Abbott et al，2016;Hoover et al，2000;Dong et al，2008;Avadhanulu，Sreenivas，2013;Gibbons，Ringdal，2006;Shelly et al，2007;Peng，Zhao，2009）。
常用地震数据去噪方法是根据信号主要特征进行单一的数据去噪，效果往往不够理想，本文考虑将不同的去噪方法穿插结合，介绍了匹配滤波法和曲波变换法相结合（以下简称“匹配曲波结合法”）的主动源气枪弱信号提取过程，并进行仿真实验和实际数据的计算分析，得到了比单一方法计算更为精确的地下介质波速变化。
地震研究43卷第4期谭俊卿等：匹配滤波和曲波去噪相结合的气枪主动源弱信号提取1方法原理简介
1.1匹配滤波基本原理
匹配滤波是以已知信号波形为模板，与各道低信噪比信号进行互相关计算，对于某一模板信号，利用下式计算模板信号和待检测信号的相关系数（Cross-Correlation coefficient，简称CC系数）：CC=∑t1t0[X（t）-X]*[Y（t）-Y]∑t1t0[X（t）-X]2*∑t1t0[Y（t）-Y]2（1）式中：t0，t1是计算相关系数窗口的开始与结束时刻;X（t），Y（t）为t0～t1的模板信号和待检测信号。计算时每次在待检测信号上移动一个采样点，得到该分量的一道相关系数波形（李璐等，2017;Peng，Zhao，2009;Yao et al，2015）。
本文模板信号通过数据叠加得到，根据模板信号对原始信号进行匹配识别，计算得到相关系数，最后将模板信号与相关系数褶积得到恢复信号。相关系数对恢复信号的信噪比有着直接的影响，系数加窗可实现噪声的压制。Kaiser窗是一种最优化窗，其频带内能量主要集中在主瓣中，有着最好的旁瓣抑制性能，本文用该窗函数对相关系数中的干扰进行初步压制。相关系数处理的另一方法是曲波变换，曲波变换是一种良好的二维数据稀疏处理方法，将一维相关系数转为二维数据格式，进行曲波处理，可去除相关系数中的大部分干扰。基本处理流程如图1所示。1.2曲波变换基本原理
曲波变换和小波变换、脊波变换原理一样，都属于系数理论范畴，通过基函数与信号的内积来实现信号的稀疏表示。在二维连续空间R2中，r，θ为频域中的极坐标，ω为频域变量，x为空间位置变量。定义频域中平滑、非负、实值的半径窗W（r），角度窗V（t）满足（Candes，et al，2004;申阳，2011;柳慧谱等，2014;庄哲民等，2014;仝中飞等，2008）：∑∞j=-∞W2（2jr）=1r∈（34，32）
∑∞l=-∞V2（t-l）=1t∈（-12，12）（2）式中：j为半径;l为角度参数。
对所有的j≥1，定义Fourier频域中频域窗为：Uj（r，θ）=2-34jW（2-jr）V（2[j2]θ2π）（3）式中：[j2]表示j2的整数部分;Uj在频域极坐标中的支撑区域为一个“楔形窗”，如图2所示，在频域的第j级、第k角度，位置为k=（k1，k2）∈Z2的Curvelet变换系数定义为：C（i，j，k）=1（2π）2∫f^（ω）Uj（Rθω）ej
笛卡尔坐标中的局部频域窗定义为：Uj～，l=Wj～（ω）Vj（Sθlω）（5）式中：Wj（ω）=[2j+1-2j（ω）]
Vj（ω）=V（2[j2]gω1ω2）（j≥0）（6）式中：定义为2个一维低通窗口的内积。
离散Curvelet 变换定义为：C（j，l，k）=∫f～（ω）Uj～（STθlω）ej（S-Tθlb，ω）
=∫f～（Sθlω）Uj～（ω）ej（b，ω）dω（7）离散曲波变换算法有2种，本文使用的是基于Wrapping的快速曲波变换算法，其运算速度快、算法效率高。
地震信号沿时空域的不同角度具有强相关性，当曲波变换系数分布在曲波域的有限区域时，系数为强振幅，随机噪声在时空域中不具有相关性;当曲波变换系数分布于整个曲波域时，系数为弱振幅，利用信号和噪声在曲波系数的幅值差异，去除一定比例的与噪声干扰相对应的弱幅值系数，进行曲波重构，实现曲波去噪。
在曲波数据处理中，如何在控制保留系数的比例、尽可能在保留有效信息的前提下实现数据去噪，是曲波处理中的难点，根据曲波变换的特性，曲波变换域中约10%的曲波系数就几乎能完整地恢复出信号的各项信息。对低信噪比信号通过曲波变换得到的相关系数，这个比例甚至可以缩小到1%～5%，再经与模板信号的褶积计算，仍可较好地恢复出信号波形和对应的走时差变化。匹配曲波结合法变相使曲波变换的稀疏表达能力增强，可以用更少的、与信号相关的系数进行数据重构恢复，提高数据中噪声与信号的分离度与信号信噪比。
2仿真数据计算
2.1仿真数据
本文的研究方法需要数据的初至具有良好的时间一致性，信号道间具有高度重复性，不同信号道间尽量满足近等间隔采样，这样才能得到更准确的数据处理结果。主动源气枪格林函数在气枪源激发条件不变时具有高重复性，信号间的相关性可高于0.99，满足高重复性要求。根据等间隔采样要求，首先模拟3年每天一次的气枪主动源激发过程，二维数据格式如图3a所示，并在每天的模拟信号中添加波速相对变化，该变化包含了周期为一年的大变化和周期为一个月的小变化（图3c），每道信号采样率为100 Hz，采样时长为14 s;对仿真信号添加高斯白噪声，本文以-5 dB信噪比数据（图 3b，d）为例说明。
2.2匹配曲波结合法与其它方法的比较
对幅值归一化信号应用一维频域滤波、曲波去噪、匹配滤波及匹配曲波结合法进行对比说明。应用频域滤波对信号主要频段2～7 Hz进行滤波处理;匹配滤波只进行简单加窗处理;曲波去噪保存约5%的大幅值曲波域系数。
按照图1所示流程，对-5 dB模拟信号处理计算，匹配识别的相关系数在曲波处理前后如图4所示，图4a为匹配识别相关系数第1道，图4c为匹配二维相关系数，可见较明显的噪声干扰;图4b为曲波处理后相关系数第1道，图4d为曲波处理后的二维相关系数，可见噪声干扰减少，最后用图4d对应的相关系数与叠加模板信号褶积计算。对于高重复信号，利用模板信号计算的相关系数在理论上是左右对称的，本文只采用右边一半的相关系数进行计算，即可得到对应的重构信号。不同数据处理方法的结果如图5所示，频域滤波去噪后的二维数据噪点最多（图5a），曲波去噪中有小部分亮点出现（图5b），其它2种方法噪点不明显（图5c，d）。在不同方法恢复的单道与理论信号的对比中（图6），频域滤波处理后信号残留了较多干扰，曲波去噪数据波形与理论变化最接近;匹配滤波与匹配曲波结合法结果相当，但在部分位置出现了信號幅值和相位的异常变化，这两种方法在计算过程中都使用到了叠加参考模板信号参与的褶积计算，导致了重构信号的固定偏差，重构信号在应用互相关法计算信号波速相对变化时对波速变化的计算几乎无影响。
根据恢复信号计算得到精确的信号走时差变化是信号处理的最终目的。对于不同滤波方法，结果不同方法计算得出对应的波速变化曲线，如图7所示。分析可知：频域滤波的变化最杂乱，匹配滤波略优于频域滤波;曲波去噪总体年变化趋势与理论值相当，月变化细节信息丢失;匹配曲波结合法去噪数据的波速相对变化与理论值最接近，仅在信号波速相对变化幅值上出现细微差别。其中，曲波去噪结果对应的单道恢复波形中存在少量噪声干扰，说明此时曲波系数幅值的保留量已经处于主要信号和噪声临界点位置，而曲波去噪后波速变化的月变化信息被作为噪声干扰而去除，增大系数保留量则会增加干扰，一般的曲波系数保留方法难以实现对干扰信息的精确去除。综合对比认为，匹配曲波法可较好地克服匹配滤波去噪能力较弱和曲波去噪细节丢失的问题。
为了定量地对数据结果进行分析对比，制定评价参数如下：Eva=Corr·lg（SNR）Resid（8）式中：SNR为恢复信号的信噪比;Corr为计算的波速变化与理论值的相关系数，0 分析表1中不同方法的评价结果发现，匹配曲波结合法在处理信号的相对波速变化的恢复上性能最好，曲波去噪恢复数据的信噪比最高，但在波速相对变化中细节信息丢失较严重，匹配曲波结合法信噪比表现一般，与褶积重构时的固定偏 差对恢复波形的影响有关，此偏差使得重构恢复信号的相位和振幅与理论信号间存在差异。 2.3不同去噪方法的综合评价分析 计算不同方法在不同信噪比下恢复数据的综合评价值Eva，如图8所示。由图可见，频域滤波最低，曲波去噪次之，匹配曲波结合法最高。随着信噪比的增加，频域滤波、匹配滤波、曲波去噪法Eva值均近似线性增加，说明这3种方法受原始数据信噪比的影响较大。而匹配滤波和曲波变换结合法在仿真实验中可对-8 dB信噪比信号进行数据的恢复重构，且能较精确地计算出波速相对变化。 3实际数据处理 以云南省内距宾川主动源气枪发射台约50 km处的53284台记录数据为例进行计算说明，选取2012—2016年部分记录数据（图9），数据采样率为100 Hz。53284台距离激发台较远，信号较弱，噪声干扰多（图9a）。首先使用频域滤波初步处理并叠加得到模板信号（图9b）;计算模板信号与各道信号的相关度，去除约10%相关性较低、异常信号道;分别应用曲波去噪、匹配滤波、匹配曲波结合法等处理，结果如图10所示。由图10a的二维去噪结果可见，频域滤波效果最差，二维信号和单道信号噪声残留多（图10a-1，b-1）;曲波去噪后数据同相轴清晰度增加（图10a-2），单道信号干扰减少（图10b-2）;匹配滤波重构信号同相轴较清晰（图10a-3），观测对应单道信号可以发现，此时恢复信号与叠加模板信号相似度最大（图10b-3）;匹配曲波结合法结果中，同相轴连续性进一步增强（图10a-4）。在匹配滤波的褶积信号重构中，需用模板信號参与计算，而叠加得到的模板信号中难去除的干扰会被带入到重构信号中，干扰重构信号的波形变化，所以对图10a-4结果再次应用曲波去噪做适当处理，结果如图10a-5所示，此时单道信号噪声干扰明显减少（图10b-5）。对上述数据的处理结果分别计算其波速相对变化曲线，并与该地的水位变化曲线进行对比，如图11所示。由图11a可见，频域滤波和匹配滤波结果计算的波速变化较杂乱、干扰较多，波形变化趋势难以分辨;而曲波去噪优于前2种方法，但计算的波速变化率曲线中仍残留了较多噪声干扰（图11b）;匹配曲波结合法处理的数据计算的波速变化与再经曲波去噪处理的结果相接近，但后者的异常变化更少（图11b中绿线）。同时该地区的水位变化曲线与波速变化率曲线具有近似正相关关系，特别是2014年6月前波速和水位的总体变化趋势较接近。 4结论与讨论 本文以宾川主动源气枪信号为基础，介绍了一种匹配滤波和曲波变换相结合的、在强背景噪声下提取地下介质波速变化弱信号的方法，可提高强背景噪声下信号的信噪比和地下介质波速相对变化的计算精度。通过实验分析得出以下结论： （1）通过仿真实验发现，匹配曲波结合法在信号波速相对变化的恢复精度上远高于单独使用曲波去噪或匹配滤波方法，波速相对变化与理论变化接近程度高。 （2）以云南省内距宾川主动源气枪发射台约50 km处的53284台记录数据为例，处理实际低信噪比地震数据得到波速相对变化曲线，并与水位变化曲线对比，发现53284台的波速变化曲线与水位变化曲线存在近似正相关关系。 仿真数据是一天一次激发的等间隔采样数据，而实际数据记录受采集环境限制，为非等间隔采样数据。二维数据处理过程中，将非等间隔采样数据作为等间隔数据处理后又恢复到非等间隔状态过程中，在时间跨度比较大的地方无法准确估计出波速相对变化的差异，如空缺几个月数据的波速变化曲线前后值比较接近。低信噪比信号叠加得到的模板信号存在较多难以去除的干扰，这些干扰会影响匹配滤波的信号恢复重构，降低恢复信号的精度。 参考文献： 陈颙，王宝善，葛洪魁，等.2007a.建立地震发射台的建议[J].地球科学进展，22（5）：441-446. 陈颙，王宝善，姚华健.2017b.大陆地壳结构的气枪震源探测及其应用[J].中国科学：地球科学，47（10）：1153-1165. 陈颙，张先康，邱学林，等.2007c.陆地人工激发地震波的一种新方法[J].科学通报，52（11）：1317-1321. 黄亦磊，运乃丹，周仕勇.2017.台差法提高低信噪比下走时变化探测精度的研究[J].地震研究，40（3）：595-604. 李璐，王宝善，侯金欣.2017.模板匹配滤波技术在地震数据处理中的应用[J].中国地震，33（1）：14-22. 林春，王绪本.2009.第二代Curvelet变换在地震随机噪声衰减中的应用[J].计算机工程与应用，45（25）：222-224. 林建民，王宝善，葛洪魁，等.2010.大容量气枪震源子波激发特征分析[J].地球物理学报，53（2）：342-349. 刘必灯，李小军，周正华，等.2011.气枪震源地震动效应分析[J].地震学报，33（4）：539-544. 刘自凤，苏有锦，王宝善，等.2015.宾川主动源地震波走时变化分析方法研究[J].地震研究，38（4）：591-597. 柳慧谱，吕健，潘伟杰，等.2014.基于第二代curvelet变换的自适应阈值去噪研究[J].微型机与应用，（10）：76-79. 栾奕，杨宏峰，王宝善.2016.大容量气枪主动源波形资料处理（一）：云南宾川[J].中国地震，32（2）：305-318. 丘学林，陈颙，朱日祥，等.2007.大容量气枪震源特征及地震波传播的震相分析[J].地球物理学报，51（1）：206-212. 申阳.2011.基于Curvelet变换的图像去噪算法研究[D].哈尔滨：哈尔滨工程大学. 仝中飞，王德利，刘冰.2008.基于Curvelet变换阈值法的地震数据去噪方法[J].吉林大学学报（地球科学版），38（S1）：48-52. 王宝善，葛洪魁，王彬，等.2016.利用人工重复震源进行地下介质结构及其变化研究的探索和进展[J].中国地震，32（2）：168-179. 王彬，杨润海，王宝善，等.2012.地震波走时变化精确测量的实验研究[J].云南大学学报（自然科学版），34（S2）：168-179. 王伟涛，王宝善，蒋生淼，等.2017.利用气枪震源探测大陆浅部的地震学研究回顾与展望[J].地震研究，40（4）：514-5123. 魏芸芸，王海涛，苏金波，等.2016.新疆2次中强地震前气枪震源反射波震相走时异常变化初步研究[J].中国地震，32（2）：270-281. 向涯，王彬，杨润海，等.2017.气枪主动源与天然地震信号的传播特征对比研究[J].地震研究，40（4）：605-612. 徐逸鹤，王宝善，王伟涛.2016.利用固定台站分析长江激发气枪信号特征[J].中国地震，32（2）：282-294. 杨微，王宝善，刘政一，等.2016.不同激发环境下井中气枪震源特征研究[J].中国地震，32（2）：231-240. 姚佳琪，曹文忠，叶泵，等.2017.零相位化和非线性叠加对提取气枪震源地震波信号波形的影响研究[J].地震研究，40（4）：581-594. 袁艳华，刘伊克，王一博.2010.基于curvelet变换的角度控制去噪方法研究[C]// 中國地球物理2010——中国地球物理学会第二十六届年会、中国地震学会第十三次学术大会论文集. 翟秋实，姚华建，王宝善.2016.气枪震源资料反褶积方法及处理流程研究[J].中国地震，32（2）：295-304. 庄哲民，姚伟克，杨金耀，等.2014.基于自适应阈值的Curvelet医学超声图像去噪算法[J].中国医疗器械杂志，38（6）：398-401. Abbott B P，Abbott R，Abbott T D，et al.2016.Observation of gravitational waves from a binary black hole merger[J]，Phys Rev Lett，116（6）：061102. Avadhanulu J，Sreenivas T V.2013.Matched filtering approach to robust speech recognition[J].Journal of the Indian Institute of Science，79（3）：185. Candes E J，Demanet L，Donoho D L.2004.New tight frame of curvelets and optional representation of objects with piecewise singularities[J].Pure and Applied Mathematics，57（2）：219-266. Chen Y，Liu L B，Ge H K，et al.2008.Using an arigun array in a land reservoir as the seismic source for seismotectonic studies in northern China：experiments and preliminary results[J].Geophysical Prospecting，56（4）：601-612. Dong F，Pierpaoli E，Gunn J E，et al.2008.Optical cluster finding with an adaptive matched-filter technique：algorithm and comparison with simulations[J].The Astrophysical Journal，676（2）：868. Emmanuel C，Laurent D，David D，et al.2006.Fast discrete curvelet transforms[J].Multiscale Modeling & Simulation，5（3）：861-899. Gibbons S J，Ringdal F.2006.The detection of low magnitude seismic events using array-based waveform correlation[J].Geophys J Int，165（1）：149-166. Hoover A，Kouznetsova V，Goldbaum M.2000.Locating blood vessels in retinal images by piecewise threshold probing of a matched filter response[J].IEEE Transactions on Medical Imaging，19（3）：203-210. Peng Z G，Zhao P.2009.Migration of early aftershocks following the 2004 Parkfield Earthquake[J].Nature Geoscience，2（12）：877-881. Shelly D R，Beroza G C，Ide S.2007.Non-volcanic tremor and low-frequency earthquake swarms[J].Nature，446（7133）：305-307. Simmons G.1964.Velocity of shear waves in rocks to 10 Kilobars[J].Joural of Geophysical Research Atmospheres，69（6）：1123-1130. Wang B S，Yang W，Yuan S，et al.2010.An experimental study on the excitation of large volume airguns in a small volume body of water[J].Journal of Geophysis and Engineering，7（4）：388-394. Yang H F，Zhu L P，Chu R S.2009.Fault-plane determination of the 18 April 2008 Mount Carmel，Illinois，Earthquake by detectingand relocating aftershocks[J].Bull Seism Soc Am，99（6）：341-3420. Yao D D，Peng Z G，Meng X F.2015.Remotely triggered earthquakes in South-Central Tibet following the 2004 MW9.1 Sumatra and 2005 MW8.6 Nias Earthquakes[J].Geophys J Int，201（2）：543-551. Yukutanke H.1988.In Situ Measurements of elastic wave velocity in a Mine，and the effects of water and stress on their variation[J].Tectonophysics，149（1-2）：165-175. Combining MFT and Curvelet Transform Method to Extract Weak Signal in Active Source of Air Gun TAN Junqing1，YANG Runhai2，XIANG Ya3，WANG Bin2，JIANG Jinzhong2 （1.School of Earth Sciences，Yunnan University，Kunming 650091，Yunnan，China） （2.Yunnan Earthquake Agency，Kunming 650224，Yunnan，China） （3.Key Laboratory of Earthquake Geodesy，Institute of Seismology，China Earthquake Administration，Wuhan 430071，Hubei，China） Abstract In this paper，we introduced a data processing method interspersed with two method that based on one-dimensional template Matching Filtering Technology（MFT）and two-dimensional Curvelet Transform method.Firstly，the correlation coefficients are obtained by the one-dimensional MFT，then the correlation coefficients are composed into two-dimensional data and processed by Curvelet Transform method.Finally，the corresponding correlation coefficients are respectively folded with the template signal to obtain the recovery signal.Then，we applied this method to the processing of simulation signal and air gun signal of Binchuan active source.Experiments show that this method is more better than that processing by one single method，with better recovery ability in Binchuan Air Gun Experimental Base，which can get better recovery signal in simulation and actual recorded data signal processing，and get the wave velocity variation with less noise interference in underground media.It is beneficial to the analysis and research of the following seismologists and improve the utilization and analyzability of low SNR data. Keywords：weak air gun signal of active source;Matched Filtering Technology;Curvelet Transform denoising;velocity changes
