被动源面波和体波成像在内蒙古浅覆盖区勘探应用

刘国峰， 刘语， 孟小红*， 刘澜波， 苏维俊， 王永志， 张致付

1 中国地质大学(北京)地球物理与信息技术学院， 北京 100083 2 美国康涅狄格大学地质系， 美国康涅狄格州 06269 3 吉林大学地球探测科学与技术学院， 长春 130061

0 引言

随着矿产勘探程度的不断深化，矿集区“透明化”立体探测成为国家深地资源勘探领域的重要研究方向.其对地球物理探测的重要需求是在3000 m深度范围内实现高精度、高分辨率成像.反射地震勘探具有深部高分辨探测能力，但是三维大规模应用于矿产勘探存在成本高、采集困难等问题.目前在矿集区立体探测中，一种典型方案是采用网状分布的二维反射地震，建立勘探区框架结构模型，或以解释的结构和物性成果为约束，辅助面积性的重、磁反演，完成矿集区三维立体结构或者物性填图.该方案利用了反射地震的勘探深度和分辨率优势，但当研究区结构横向变化大，地质结构复杂时，有限的地震剖面无法控制区域构造变化.因此，在矿集区实施经济、有效的三维反射地震勘探是目前立体探测的待突破瓶颈之一(吕庆田等，2019).

被动源地震是指不依赖主动源激发地震波，而是通过固定位置的地震仪台站/检波器等，在一段时间内连续接收天然噪声信号进行探测研究的技术方法.它是地震学进行大尺度地壳和岩石圈结构成像，反演地球深部物性的主要手段(Bostock et al.，2001;Moschetti et al.,2007;余大新等，2016).若将台阵更密集的摆放，可完成更精细尺度的探测应用，例如密集台阵和城市环境中进行的面波勘探(Wapenaar，2006;Wapenaar et al.，2010；王爽等，2018；Wang et al.,2020).相对于上述应用，从被动源中提取体波记录进行反射波地震勘探的应用，在近年成为了研究热点.

Claerbout早期曾指出：“通过地表两点放置的两个检波器接收的天然震动信号进行相关计算，可以构建以一点为检波器，另外一点为震源的主动源波场记录.”(Rickett and Claerbout，1999).同时提出该计算思路的还有Cole(Cole，1995)，不同之处在于，Cole是指提取面波，而Claerbout指的是提取体波.基于天然噪声提取体波的核心环节为相关计算(Interferometry)(Schuster，2001；Schuster et al.,2004;Schuster and Zhou，2006)，其由超声计算领域发展而来(Weaver and Lobkis，2001,2002)，主要内容为：相邻的两个检波器信号进行互相关计算等于两个检波器间激发接收格林函数与随机震动信号的褶积.

随着节点检波器的发展和普及(史子乐等，2019)，连续接收被动源地震记录成本大幅度降低，被动源反射地震探勘从理论走到了实际应用.自2007年以来，陆续有试验性应用文章发表，其普遍特点是在主动源地震采集时，进行重合测线的被动源采集，通过主动源和被动源的单炮及成像剖面对比，测试被动源反射波勘探在不同领域应用的可行性，这其中包括在油气勘探领域的应用(Draganov et al.,2009,2013)、二氧化碳存储四维动态监测的应用(Xu et al.,2012；Boullenger et al.,2015；Cheraghi et al.,2017)以及金属矿勘探的应用等(Cheraghi et al.,2015；Olivier et al.,2015；Roots et al.,2017).同时，正演计算的开展也在理论与实践应用之间搭建了论证纽带，分析了震源分布、震源长度、震源数量等对被动源成像质量的影响(Thorbecke and Draganov，2010；朱恒等，2012；张盼等，2015).同时，更精细化的被动源勘探地震数据处理方法也不断提出，以增强弱体波信号，进而提高反射波成像质量，这其中包括照明分析方法(Vidal et al.,2011)以及在原始道集数据进行的各种滤波处理(Girard and Shragge，2019a，b).

被动源地震勘探无需主动源激发，其采集成本不及主动源十分之一，是目前最具潜力的在活动矿区进行三维勘探的替代选择.但同时，被动源地震勘探也面临体波信号弱，现有可参考案例少等问题.基于此，我们在内蒙古集宁浅覆盖区进行了被动源地震勘探试验.本文推导了基于声波方程的被动源反射信息提取理论.在实际数据验证阶段，介绍了数据的采集、频率域面波和体波分离以及对拟炮集记录的处理等关键环节.面波反演和体波成像结果表明，被动源勘探在内蒙古浅覆盖区的应用具有广泛的前景.

1 被动源地震勘探拟炮集记录提取的理论基础

相关计算是被动源反射波勘探成像的计算基础，其可通过互易性定理，基于声波方程推导而来.

点源激发下的频率域声波方程可表达如下：

=-jωρδ(x-xA)，

(1)

其中，G(x,xA,ω)指震源在x处，检波点在xA,处的波传播格林函数的频率域表示，ρ(x)和v(x)分别是传播介质的密度和速度，ω是角频率，δ(x-xA)是脉冲震源.

根据互异性定理(Wapenaar,2006)，同属接收点的A,B两点间格林函数G(xA,xB,ω)为:

-(∂iG(xA,x，ω)*)G(xB,x，ω))nid2x,

(2)

式中∂D是任意封闭区间，ni=(n1,n2，n3)表示该区间的法向向量，R{}为取实部计算.公式(2)是A，B两点间格林函数的准确表达，其中包含了非均匀介质中传播的直达波、散射、反射以及多次反射等.但在实际应用中，往往对其进行简化后应用，其中包括假设应用介质为均匀介质、密度和速度为常量、∂D边界处介质光滑等(Wapenaar and Fokkema,2006),而公式(2)也随之简化为：

(3)

上述简化过程虽然会为波场带来异常振幅，但公式(3)中波场相位没有变化，因此，依然可以用来进行结构性成像计算.

被动源勘探中，震源为随机噪声，若其震源函数为N(x,ω)，则A，B两点的接收该震动的记录为：

u(xA,ω)=∮∂DG(xA,x，ω)N(x,ω)d2x,

u(xB,ω)=∮∂DG(xB,x，ω)N(x,ω)d2x,

(4)

若N(x,ω)为随机信号，在时间和空间上均互不相关，则满足：

〈N(x,ω)N(x′,ω)*〉=δ(x-x′)S(ω),

(5)

其中，S(ω)是随机噪声的功率谱.〈·〉代表空间的加权平均.

将公式(4)和(5)代入公式(3)中，可得：

(6)

式中，左端项是指接收点A，B之间的格林函数与噪声震源的褶积.右端项是检波器A，B记录数据的相关计算，其中，〈*〉代表不同互相关结果的叠加.

被动源反射波勘探中，依照公式(6)，固定位置A,将A处的震动记录与移动的B点震动记录分别进行固定窗口的互相关和叠加计算，可获得形如炮点位置在A处，接收点在变化的系列B处的炮集记录.因为A处并无实际的主动源激发，因此该记录称为拟炮集记录.按设计的观测系统变化A的位置以及与其互相关计算的移动B的范围，最终可获得勘探区的系列拟单炮数据集.图1表示了被动源体波成像计算的原理和步骤，其中图1a是被动源记录生成拟单炮记录的原理示意图，被动源在A处的记录与其经界面反射后的B处的记录褶积，可获得A处激发，B处接收的反射记录.图1b示意了A，B两点进行相干计算的过程.选定窗口长度，进行对应时间的相关和叠加，最终形成如图1c的拟单炮记录.将炮集记录按主动源反射地震流程处理，可获得反射成像剖面.

图1 被动源记录相干计算生成拟炮集记录示意图Fig.1 Schematic diagram of cross-coherent calculation of passive seismic records to generate virtual-shot records

上述应用随机噪声记录生成拟炮集记录计算的过程可以应用图2的模型，采用正演计算进一步说明.模型含三个不同速度和密度的层位，其中8000个地下震源随机分布.正演计算采用声波方程的有限差分计算完成(Thorbecke and Draganov，2011).地表记录随机噪声的检波点间距为50 m，记录长度为4 s，图3a为地表记录的随机噪声记录，图3b为震源在模型正中，所有检波点参与计算生成的拟炮集记录.若变换震源位置，可以获得按特点观测系统采集的炮集记录，进行常规处理可获得反射波成像结果.

图2 含有三个层位的模型，其中8000个震源随机分布Fig.2 A model with 3 layers, there are 8000 random sources in the model

图3 正演获得的随机噪声记录以及生成的拟炮集记录Fig.3 The random noise obtained by forward modeling and the generated virtual shot gather

在上述理论和模型验证的基础上，在内蒙古浅覆盖区对被动源成像方法进行了数据试验.

2 内蒙古浅覆盖区被动源勘探应用

2.1 被动源数据的采集

本文试验所用的被动源数据采集区域行政区划隶属于内蒙古集宁丰镇市官屯堡乡，如图4所示.该地区富含钼等多种金属矿产，属岩浆热液型矿床.在前期地质调查项目的支持下，该区进行了主动源二维地震采集工作.该次被动源试验选择与前期该区一条二维地震测线重合，可对成像结果进行互相验证、对比.

图4 数据采集区域位置图，图中红线为本次被动源采集的2 km测线Fig.4 Position of the passive seismic line. The red line is our 2 km length passive seismic line

被动源数据采集采用节点式单分量检波器，型号为SmartSolo的IGU-16.该节点检波器自然频率为5 Hz，1 ms采样下可连续采集25天.本次共摆放检波器100个，检波点间距是20 m.采集时设置采样间隔为1 ms，采集时长为连续的10天.为了能够避免地表风噪，将检波器埋于地下，埋藏深度约20 cm.

该测区无明显强活动震源.向北2.4 km为S24省道，其他方向交织村级公路.白天省道车流量相对较大，多有运煤卡车经过.采集时值秋季收割季节，有收割机在农田作业，附近有农民刨挖土豆.上述是地表可见的自然震源.附近矿场处在停工状态，无矿区活动.

图5是近村级公路500 m处检波器10月5日和6日两天的震动记录，总体白天震动强度要大于晚上，主要以公路交通和人文活动产生的噪声为主.

图5 某路边检波器记录的两天的被动源记录Fig.5 Two days of passive seismic records by a roadside receiver

2.2 基于频率域信噪比计算的面波和体波分离

被动源体波成像的理论计算中，假设天然震源在地下半空间均匀分布，但实际应用中，地表附近震源(以公路交通为主)无论从强度还是数量上均占优，导致面波在被动源信号中的能量强，而相对而言，体波则属于弱信号.将体波从面波中分离出来，有助于提高体波信号成像信噪比.在重建的拟单炮记录中，因为体波信号频率低，与面波在位置和频带范围都存在重合，滤波分离效果不佳.若从原始数据出发，在互相关前对时窗内的数据进行面波和体波成分的甄别和分离，再进行相干和叠加，则更为有效.

分析公共交通产生的面波特点，可见其具有与体波信号不同的频谱形态.图6a,b分别是面波和体波占优数据段的信号的频谱.面波占优时，频谱呈梳状分布，而体波占优时，频谱形态正常.这种差别可通过频率域的信噪比量化表示.面波占优时时窗数据的信噪比比体波占优时低，可基于此特点，对面波和体波数据进行甄选.本文采用张军华等(2009)频率域信噪比的计算公式，对参与互相关计算的时窗内信号信噪比进行计算，再根据选定的信噪比阈值，判断相关叠加结果分配到体波数据中还是面波数据中.公式(7)是频率域计算信噪比的公式，图7是其参数在功率密度谱上的表示.

图6 被动源信号不同频率成分的功率密度谱特征Fig.6 Power density spectrum characteristics of different frequency components of passive seismic signal

(7)

根据不同工区数据的统计和分析，确定甄选体波和面波的信噪比阈值ε，甄选和分流数据，分别相干计算得面波和体波为主的拟炮集记录，过程概括为公式(8)，

(8)

其中，Gk(r1,r2,t)是r1和r2两个检波器位置震动记录的相关叠加结果，t是记录时间，k表示波类型，当k为0时表示为面波，为1时表示为体波.n是叠加的总窗口数，ai(r1,r2,t)是某个时段的互相关结果.

图7 频率域计算信噪比参考图Fig.7 Spectrum for signal-to-noise ratio calculation in frequency domain

被动源数据量大，生成拟炮集过程的计算强度也很大.以本文采集为例，100个检波器，10天接收时间，1 ms采样接收的记录超过500 GB.公式(6)中的互相关在频率域计算比时间域具有更高的计算效率，辅以图形处理器GPU上的批量并行傅里叶变换，可实现更快的急速计算.本文介绍的基于频率域信噪比的面波和体波分离，可插入到上述傅里叶变换流程中进行，无需额外计算.

应用上述方法对本文采集的被动源数据进行拟单炮生成计算，计算时，互相关计算时窗长度选择为5 s，并通过上述的频率域信噪比方法进行面波和体波的分离，其中,各频率值为f1=0,fL=2,fH=15,f2=40,fc=60，信噪比阈值测试后选择为0.76,小于0.76的归于面波叠加，大于0.76的部分归于体波叠加.图8a是同一炮点计算获得的以面波为主的单炮，图8b是分流面波后的突出体波的记录，虽然其中依然有面波出现，但是相比于图8a而言，有明显的反射波能量出现.

图8 频率域分离计算后的拟炮集记录(a) 以面波为主的拟炮集记录; (b) 分离面波后突出体波的拟炮集记录，红色箭头为反射波信号.Fig.8 Virtual shot gathers after frequency domain separation calculation(a) Shot gather with surface wave as the main energy; (b) Shot gather with body wave as the main energy after surface wave separation， the red arrows point to reflection wave.

2.3 利用高频面波反演覆盖层厚度

被动源勘探处理过程能获得的辅助信息远少于主动源地震，因此，需要充分发掘不同信息成分的表征能力.在主动源地震勘探实践中，面波主要用来进行浅层结构探测，被动源重建的拟单炮记录面波，同样可以进行近地表结构探测应用.而在本文的勘探区域，面波反演结果可用来进行覆盖层厚度的估算.

面波反演的过程分为面波的识别，频散曲线的提取和横波速度反演三步.本文所用面波即为图8a所示计算的面波记录.图9是单炮面波记录和从中提取的频散曲线，从中可以看出，低阶面波频散曲线聚焦性、连续性都很好，易于识别，且可见相对连续的高阶面波.和主动源面波反演一样,一旦从被动源面波中提取了频散曲线,接下来就可对频散曲线进行反演以获取横波速度结构.本文采用遗传算法反演,遗传算法作为一种非线性全局优化方法,不要求苛刻的初始速度和深度模型,得到的最佳拟合模型稳定、可靠(赵东等，1995).反演中，仅应用了低阶面波，采用一维反演计算.对测线上所有单炮记录进行上述反演，综合成图，可得到图10所示的横波速度剖面.该剖面中，约100 m深度处可见明显的速度结构纵向变化界面，可判断其为覆盖层分界面，这与已知的地质符合良好(严昊伟等，2017).覆盖层厚度的准确计算对区域结构认识，以及面积重、磁的反演等都可提供约束.

图9 面波炮集记录及其频散图Fig.9 Shot gather with surface wave and its frequency-velocity spectrum

2.4 被动源反射波勘探成像

面波的反演是被动源勘探的传统应用领域，而基于被动源的反射波勘探则是目前的研究热点.经面波分离后的拟炮集记录虽然反射依然较弱，但可应用反射地震资料处理的多次叠加技术增强反射能量.在拟炮集计算过程中，设定炮间距和炮检距范围，则可形成测线反射勘探的观测系统.表1是本文采集数据生成的拟单炮记录的观测系统参数.为了提高反射波的信噪比，每个检波点位置都进行拟单炮记录的计算，图11是测线覆盖次数图，其中最大覆盖次数达到100.

图10 面波反演的横波速度模型Fig.10 Inverted shear wave velocity model with the surface wave

图11 拟炮集记录覆盖次数图Fig.11 Fold number of the virtual shot gathers

表1 被动源生成的拟单炮记录观测系统主要参数Table 1 The main geometry parameters of the virtual shot gathers from passive seismic

拟炮集数据处理采用主动源反射地震数据的处理流程，主要流程包括野外静校正，本文以高程校正为主，振幅垂向补偿和横向均衡，去噪，反褶积，速度分析及叠加处理.在去噪中，只进行了异常振幅消除，因为在共中心点道集中，同相轴的连续性差，因此没进行剩余静校正处理.考虑到速度对成像结果的影响，在与主动源对比时，只进行叠加剖面的对比.图12是处理中速度分析和局部叠加结果图，速度谱中有明显的能量团，相对应的共中心点道集有相对应的反射同相轴，而局部叠加结果中，可见连续性的反射.

图12 被动源反射信号速度分析和局部叠加Fig.12 Reflection velocity analysis and local stack for passive seismic data

经过上述处理流程处理，图13a是获得的被动源记录叠加剖面，而图13b是与之重合的主动源成像结果.被动源成像结果频率较低，所以以主动源记录上的明显反射为标准进行成像结果的对比.总体可见二者明显反射同相轴都具有较好的一致性.首先A处均无明显反射，且反射的起始位置一致.B和C处的同相轴形态和位置均较一致.D处被动源和主动源斜率稍有差别，但横向上的延展一致.E处被动源剖面上有与主动源对应的同相轴，但是斜率上稍有差别，被动源可能因为频率低的缘故，可见同相轴横向上的延展和连续.B和C相接处出现错段，可能解释为断层F，这种错段在被动源剖面上也可以清晰可见.

图13 被动源(a)和主动源(b)反射波叠加剖面对比Fig.13 Stack sections of (a) passive seismic and (b) active seismic for reflections

综合上述的分析对比，虽然该区无明显地下主动源记录，且分离的面波表明天然噪声中以面波为主，但经过去面波的拟炮集反射波成像结果依然与主动源记录有较好的一致性，这为被动源进一步应用提供了参考.但同时，被动源成像剖面的低频特性一定程度上限制了其单独应用，但在主动源二维框架的约束下，进行三维被动源勘探，联合主动源结果对主要结构的解释，或可是被动源地震在该区的一个应用方向.而对被动源数据进行拓频处理，可更易于对比和解释.

3 结论和建议

本文基于被动源地震方法，在内蒙古浅覆盖区进行了勘探试验.提出了一种频率域信噪比计算的分离面波和体波信息方法，分离的面波炮集用于覆盖层厚度反演.体波炮集经反射波常规地震资料处理实现多次覆盖叠加后，获得了与主动源记录一致性较好的反射地震叠加剖面.该试验表明被动源勘探方法在内蒙古浅覆盖区具有较好的应用效果.应用被动源进行反射波勘探不需要主动源激发，其采集成本是主动源方法的十分之一，可低成本的完成三维地震勘探，结合主动源二维框架的约束，可对研究区实施更大范围的三维勘探.但同时，被动源反射波勘探还存在诸如数据频率低、信噪比低等问题，需要更多的细致的、具有针对性的处理方法进行精细处理，特别是将体波信息从强面波干扰中分离出来的方法技术，以推动被动源反射地震勘探的生产应用.