被动源反射地震勘探研究进展

阮小敏，陈明春，刘振东，王志辉，徐志伍，陈 淼，张新港

1 中国地质科学院，北京100037

2 中国地质调查局中国地质科学院 地球深部探测中心，北京100037

3 中石化石油工程地球物理有限公司 南方分公司，成都610016

0 引 言

地震勘探是地下结构探测和研究的关键方法，其中主动源是人工激发震源，其激发时刻和位置是精确已知的，信噪比高，但施工复杂，勘探成本高.而另一类被动源探测方法，无需人工激发震源，而是通过固定位置的地震仪台站或检波器等，在一段时间内连续接收天然噪声信号，然后利用地震干涉从背景噪声记录中提取波场传播信号（Campillo and Paul, 2003; Shapiro et al., 2005; Wapenaar and Fokkema, 2006）.传统上，被动源地震是开展大尺度地壳和岩石圈结构成像、反演地球深部物性的主要手段（Bensen et al., 2007）.近年来，被动源地震方法以其绿色环保、采集成本低、数据重构灵活性高和低频信息丰富等优点，在勘探地震领域受到越来越多的关注（Chamarczuk et al., 2022; Cheraghi et al., 2017）.

体波成像和面波成像是利用被动源地震勘探的基本方法.总的来说，由于面波传播距离远、衰减慢、能量大，利用噪声干涉法提取面波更容易，因此能够反演地下横波速度结构的面波勘探方法获得了更广泛应用（Chen et al., 2022; Hu et al., 2022;Le Feuvre et al., 2015; Xia, 2014）.体波信号相对较弱，从背景噪声中提取体波具有一定挑战性，但是体波具有更高的频率，考虑体波成分可以提高地下结构成像的分辨率（Nakata et al., 2015）.从被动源数据中提取体波记录进行反射波地震勘探，由于其成本相对较低，且能克服工区无法放炮问题和易于实施等优势，近些年成为了研究热点，并在地质环境复杂、散射较为严重的工区获得了较好效果，比如金属矿集区的勘探（Chamarczuk et al., 2022;Cheraghi et al., 2015; 刘国峰等, 2021; Olivier et al.,2015）、CO封存点四维动态监测（Boullenger et al., 2015; Cheraghi et al., 2017）、 油 气 勘 探（Draganov et al., 2009, 2013)和煤田勘探（Ning et al., 2021）等领域.

反射地震勘探中的高品质成像离不开高品质的原始炮集，不同于主动源激发的高信噪比炮集数据，基于被动源连续噪声信号构建的虚拟炮集往往质量不高，在实际应用中面临诸多技术挑战.比如：地下实际震源数量有限且分布不均匀，使得干涉法计算得到的虚拟炮记录中出现虚假反射；地表检波器接收到的实际背景噪声场往往主要受控于面波能量，反射体波则属于弱信号；大量检波器长时间观测的海量数据导致的计算和存储瓶颈等.针对这些挑战，近年来人们对被动源反射波勘探成像开展方法研究并取得一些进展，主要集中在不同干涉算法如何更好地重构虚拟炮集，从面波主导的背景噪声记录中有效甄别和提取反射信号，以及被动源原始数据预处理、虚拟炮集多次波压制和直接偏移成像等方面.最后通过几个实例了解被动源反射地震勘探的应用现状，并对将来的研究前景进行简单展望.

1 地震干涉法重构虚拟炮集

地震干涉（seismic interferometry）是指通过对原始地震信号进行干涉合成新记录，新记录保留了原始信号中相关性高的目标响应、并且能反映原始信号中所不具有的一些重要特征，比如地震干涉被广泛用在背景噪声面波成像、震源定位、多次波成像等诸多方面（Shapiro et al., 2005; Snieder and Vrijlandt, 2005; Wapenaar et al., 2004）.通过对被动源噪声数据进行地震干涉计算，可得到类似主动源激发的炮集记录，称之为被动源的虚拟炮集或拟炮集，而虚拟炮集的质量直接决定着最终反射剖面的成像效果，因此地震干涉法对被动源反射地震勘探而言十分关键.

1.1 地震干涉法的提出与发展

地震干涉技术的基本原理是通过对两个检波点接收到的地震记录做一定的数学运算（互相关、反褶积、互相干等），提取这两个检波点间的格林函数，得到其中一个检波点为新的虚拟震源，另一个检波点接收的模拟地震记录（图1）.该方法由Claerbout最早提出（Claerbout,1968），他阐述了通过地下埋藏震源产生的地震记录的自相关来求取地表格林函数的过程.Cole（1995）则尝试通过计算实测被动源数据来合成反射记录.Rickett和Claerbout（1999）将相关法推广到多维模型，并且称之为日光成像方法.Schuster在2001年的EAGE会议上，正式将该方法命名为地震干涉技术，并对被动源地震数据进行虚拟炮记录合成、偏移成像等，且取得了较好的效果（Schuster, 2001; Schuster et al., 2004）.Wapenaar等（2004）、Wapenaar和 Fokkema（2006）用积分理论和互易定理进行了严格的理论推导，证明了地震干涉技术在不同的非衰减介质和不同的震源条件下均成立.

图1 地震干涉示意图（修改自Schuster, 2009）Fig.1 Seismic interferometry schematic diagram (modified from Schuster, 2009)

图2 封闭区域模型的地震干涉（Wapenaar and Fokkema,2006）.格 林函数(xA,xB,ω)可由观测点 x A、xB和沿着边界∂ D 的 震源 x积分的互相关获得Fig.2 Seismic interferometry model enclosed by boundary(Wapenaar and Fokkema, 2006) Green's function(xA,xB,ω)can be obtained by cross-correlating the observations at x AandxB and integrating along source coordinate x at∂D

实际应用时通常对式（1）进行近似：（1）∂

D

边界之外是均匀介质；（2）远场近似；（3）围绕的介质参数均匀变化（Wapenaar and Fokkema,2006）.经过近似后，所得的重构两个检波点之间的格林函数公式为：

当∂

D

是一个球面，并且半径非常大时，公式（2）是比较精确的，不过当∂

D

是有限的时候，这些近似就会造成振幅上的误差（Ramírez and Weglein, 2009），但是并不影响相位，所以公式（2）在误差可接受的范围内仍适用于地震干涉方程（Thorbecke and Draganov, 2011）.该方程也由Derode等（2003）通过物理推理的方式得到.

地震干涉法的提出深化了对地震波传播规律的认识，很多方法都借助于地震干涉技术而提出或完善，大大丰富了地震学的理论与应用研究.一是基于背景噪声的研究.Shapiro等（2005）通过互相关方法提取出美国加州地震台站的瑞利波群速度频散曲线，发现反演得到的速度结构模型与研究区的构造特性对应关系较好，开启了背景噪声无源成像的应用时代.随后，利用背景噪声成像研究岩石圈速度结构呈爆发式增长（Bensen et al., 2007; Ritzwoller, 2009）.二是基于大地震尾波干涉研究.Snieder（2004, 2006）在总结前人工作基础上提出了大地震尾波干涉理论，进一步完善了地震干涉法理论.他们通过对地震台站记录的大地震尾波进行互相关，用来推断介质参数随时间的变化情况，即可实现面向目标介质的连续动态监测.尾波干涉法在震源定位（Snieder and Vrijlandt, 2005）、地震活动（Peng and Ben-Zion, 2006）、火山监测（Brenguier et al., 2008）、冰川活动（Jonsdottir et al.,2013）、地震波速动态变化（Wegler and Sens-Schönfelder, 2007）等研究领域中都有着广泛的应用.

近些年人们又在常规的地震干涉法基础上发展出了地震超越干涉法（Beyond Seismic Interferometry），该方法突破了传统地震干涉法对于虚拟震源必须置于检波器位置的限制，可实现虚拟震源或检波器在地下介质的任意位置分布.Broggini和Snieder（2012）首先将地震超越干涉法的Marchenko成像引入地球物理领域，其本质是通过虚源点与地表之间的直达波记录和地表反射响应反演得到上下行格林函数.Broggini等（2012）、Wapenaar等（2013）将该方法扩展到二维和三维，能够在无需了解上覆介质构造的情况下，基于一个宏观的非精确模型就可对目标区域的介质进行精细成像.Marchenko成像方法在多次波压制（Meles et al., 2015）、数据重构（靳中原等，2017）等方面有诸多应用，该方法也可被扩展到被动源数据，直接面向目标成像，避免检索炮集中的层间多次波影响（Jin, 2020）.

1.2 地震干涉法在勘探中的应用

地震干涉法的提出不仅极大促进了天然地震学的研究，对勘探地震领域也有重大理论和实际意义.比如在主动源地震勘探中，复杂的近地表如风化带、推覆体、小断块等通常制约着地下构造的精确成像，常规的勘探技术手段已经很难满足实际生产的需要.而地震干涉法可以通过逆时重构将井下检波器变为虚拟震源从而避免复杂近地表的影响，间接地达到面向目标成像的目的，提高了地震资料的信噪比和分辨率.由于这一过程中，不需要震源或地震子波的相关信息，该方法又称作“虚源法”（图3）.2004年，Calvert首次提出该方法（Calvert et al., 2004），随后被陆续用于实际的VSP剖面（Bakulin and Calvert, 2006）、OBC数据（Gaiser and Vasconcelos, 2010）和 OBS数据（Carrière and Gerstoft, 2013）等勘探实践.

图3 虚源法试验（曹辉等，2012）.（a）地表激发、井中接收；（b）虚源构建后的观测系统，消除了近地表复杂条件对波场的透射影响Fig.3 Virtual source method experiments (Cao et al., 2012).(a) The sources are fired on the surface, and the receivers are located in the borehole; (b) Acquisition geometry after virtual source construction, removing the effects of the complicated near surface

对于因为环保、成本等因素限制而无法实施主动源地震激发的勘探区域，被动源方法可以较好地填补这些区域地下结构成像的空白.被动源勘探采集的背景噪声信号粗看或许杂乱无章，但地震干涉法能够变无序为有序，即能够从背景噪声中重构有效的地震波场并使之成像，进而推断地下地质构造或监测储层变化等.自从 Wapenaar 等（2004）、Wapenaar和Fokkema（2006）对被动源地震干涉法在3D非均匀介质下进行严格的理论推导以来，被动源地震干涉技术不断取得进步，目前利用被动源环境噪声的面波勘探应用比较广泛，比如在城市环境中开展的工程勘察及各种无损检测（Chen et al., 2022; Cheng et al., 2015; Hu et al., 2022; 汪利民等, 2022; Xi et al., 2021）.但面波勘探方法分辨率相对较低，无法满足更高精度的勘探需求，因此被动源体波勘探逐渐兴起，这几年基于地震干涉技术的反射波提取方向成为了地震干涉领域新的热点.因为基于提取的被动源反射体波信息，不仅能得到速度结构，而且能获得高频率的反射剖面（Chamarczuk et al., 2022; Cheraghi et al., 2015,2017; 刘国峰等, 2021; Olivier et al., 2015）.

1.3 被动源虚拟炮集重构

基于地震干涉法对被动源噪声数据实际计算虚拟炮集的时候，为了提高信噪比，通常采用分段时窗互相关＋叠加来处理记录到的连续噪声记录.比如对于记录长度为

T

的

N

个台站数据，固定一检波点不动作为虚拟炮点，与剩余

N

-1个检波点在选取的小时窗

t

内作互相关，共进行

T

/

t

次互相关计算，最后将所有时窗内的互相关结果进行叠加就得到一个完整的虚拟炮集，再依次滚动计算生成其他位置的虚拟炮集记录（方捷等，2022）.计算可知，重构所有虚拟炮集需要求解

T

×

N

×(

N

-1)/2的互相关函数，对于长时间观测（一周以上）的大量检波器（上百甚至上千台）的被动源数据，总的计算和存储量都非常巨大.完成所有虚拟炮集重构之后，再借鉴常规反射地震方法进行处理，即可获得被动源反射成像剖面.由此可见，基于地震干涉法高效重构高质量的虚拟炮集是被动源反射波成像的关键.

1.3.1 多种干涉算法重构虚拟炮集

对不同检波器记录的连续波形数据进行多种干涉运算，比如经典的互相关（Schuster et al., 2004）、反褶积（Snieder, 2006）和互相干（Prieto et al.,2009）三种方法，都可以提取这两个检波点间的格林函数，进而获得地震虚拟炮集记录.

假设在检波器

r

点处的频率域波动波场

u

(

r

,

s

,ω)，是非扰动格林函数

G

(

r

,

s

,ω)和扩散波格林函数

G

(

r

,

s

,ω)分别与

s

处的震源函数

W

(

s

,ω)褶积求和得到的：

式子（3）中的

W

(

s

,ω)是一复杂的频率函数，随震源函数

s

的变化而变化； ω为频率；

G

为背景波场；

G

为扩散波场.

u

、

G

、

G

满足声波、弹性波、可衰减的波，或者是包含高频扩散的非均匀波.如果忽略频率 ω对震源函数

W

的影响，令

G

=

G

+

G

，则（3）式可简化为：

因此，频率域中A、B两点的波场互相关可表示为（Vasconcelos and Snieder, 2008）：

而反褶积干涉成像消除震源特性的影响，其在频率域可表示为（Snieder, 2006)：

同样，互相干地震干涉也消除了震源项的影响，相当于空间自相关，频率域的互相干干涉可表示为（Prieto et al., 2009)：

从公式（5）可知，互相关方法需要知道震源功率谱，但该方法稳定性高.而公式（6）和（7）中已经证明，反褶积和互相干方法不需要知道震源功率谱，且消除了震源项的影响，提高了分辨率.

很多学者曾对三种经典干涉算法的效果进行研究对比，并分析各方法的优缺点.Liu等（2020）用模型数据进行虚拟炮集重构试验，证实反褶积和互相干方法获得的虚拟炮集记录信号更接近于尖脉冲信号，分辨率和信噪比也更高（图4）.Snieder等（2009）认为反褶积方法提取的反射波响应可能会在参考检波点引入虚假反射.Cao和Askari（2019）比较不同方法效果后认为互相干方法是开展二氧化碳储存时移监测的更好选择.方捷等（2022）对内蒙地区实际数据进行分析对比后，也建议采用互相干方法，因为较之反褶积方法，它包含的波形信息更加丰富，较之互相关方法，其面波同相轴延伸更短，分布范围较小，中间无强噪声分布（图5）.

图4 不同地震干涉方法计算的模拟数据虚拟炮集（修改自Liu et al., 2020）.（a）互相关；（b）反褶积；（c）互相干Fig.4 Virtual shot gathers using noise-free data as input; generated by (a) cross-correlation, (b) deconvolution, and (c) cross-coherence (modified from Liu et al., 2020)

图5 不同地震干涉方法计算的实际数据虚拟炮集（方捷等，2022）.（a）互相关；（b）反褶积；（c）互相干Fig.5 Virtual shot gathers using low S/N data as input; generated by (a) cross-correlation, (b) deconvolution, and (c) cross-coherence (Fang et al., 2022)

在二维或三维介质中，公式（9）是病态方程，可通过最小二乘反演来求解：

其中*表示转置和复共轭，

W

是对角加权矩阵，

I

是单位矩阵， ∈是保证最小二乘反演稳定性的一个参数.这个矩阵反演的过程等同于在时间域沿着上边界∂

D

自由表面多次波响应与下传波场的反褶积.该方法也被称为多维反褶积.不规则分布震源（图6a）的地震干涉数值计算示例表明，多维反褶积方法的结果（图6d）较互相关计算结果（图6c）改善明显.但值得注意的是，多维反褶积方法需要同时对所有检波器的记录进行求逆运算，相较于传统干涉方法，在稳定性上略有不足（Minato et al., 2009）

图6 不规则分布震源的地震干涉数值计算示例（修改自Wapenaar et al., 2008）.（a）横向均匀的水平层状模型，不规则分布震源依次发出瞬态信号；（b）单个地下震源的地震记录响应；（c）互相关重构炮集（黑色点折线）与源在x=0处的直接模拟炮集（红线）对比；（d）多维反褶积炮集（黑色点折线）Fig.6 Numerical simulation of seismic interferometry with irregular sources (modified from Wapenaar et al., 2008).(a) Configuration with a horizontally layered target below a homogeneous overburden and a free surface.The irregularly distributed sources below the target emit transient signals sequentially.(b) Response of one of the sources.(c) Interferometry from cross-correlation (red traces) compared with the directly modeled response of a source at x=0 (black dashed traces).(d) Interferometry from multidimensional deconvolution

1.3.2 影响虚拟炮集重构的其它因素

除了不同干涉方法对重构虚拟炮集的影响之外，噪声震源也是影响干涉重构效果的关键因素，很多学者探究了如震源类型、震源分布、震源长度、震源数量等对重建虚拟炮集记录和最终成像的影响（Thorbecke and Draganov, 2011; 王德利等, 2012;张盼等, 2015; 朱恒等, 2012）.

首先，地下震源数量对干涉成像影响较大，较多的地下震源数量对虚拟炮集成像更有利.图7是Thorbecke和Draganov（2011）进行的一个被动源单炮模拟数值试验，对于震源在地下空间500～4 000 m随机分布的地下五层模型（图7a），将被动源重构虚拟炮集与模拟的主动源激发单炮进行比较（图7b），其中8 000个随机噪声震源重构的单炮反射较为清楚（图7c），1 000个震源重构单炮成像次之（图7d），50个震源重构单炮效果最差（图7f），且有很多假象.

其次，非相干的震源一般来说有助于地震干涉反射信号的重构，因此震源性质和位置分布也影响着虚拟炮集的成像.图7d、7f、7g的震源在数量上虽然都相同，不过震源之间的非相干性依次降低，其中图7d采用的是随机噪声震源，非相干性最强.图7f和7g的所有地下震源均采用相同的Ricker子波，但图7f的激发时间和震源位置都是随机的，仍能部分保证震源之间的非相干性，因此反射信号还可以较好地被重构；而图7g的震源均分布在2 700 m深度且同时激发，非相干性最低，因此反射数据检索最差，几乎都是假象.可见，非相干的地下震源对于构建高质量的虚拟炮集反射波十分重要.

图7 不同震源参数模拟的地震干涉结果（修改自Thorbecke and Draganov, 2011）.（a）地下震源随机分布的10×4 km数值模型；（b）地表中点激发的主动源单炮；（c-e）分别用8 000、1 000、50个随机位置地下噪声震源激发构建的被动源虚拟炮集；（f）1 000个随机位置雷克子波震源激发的重构虚拟炮集；（g） 1 000个固定2 700 m位置深度的雷克子波震源同时激发的重构虚拟炮集Fig.7 SI results from modeling by different sources (modified from Thorbecke and Draganov, 2011).(a) The 10×4 km model with sources positioned at random locations, visible as black dots.(b) A directly modeled reference result for an actual source at x=0 m.SI results for noise signature sources for varying numbers of (c) 8 000, (d) 1 000, and (e) 50.(f) SI results from 1 000 sources using the Ricker wavelet in random locations.(g) SI results by 1 000 sources using the Ricker wavelet at z=2 700 m;the sources are started simultaneously

基于被动源单炮数值模拟试验可知，噪声源分布和噪声源类型对于虚拟炮集的反射波提取起着关键作用，而数值试验本身也可为实际被动源数据的采集部署提供一定依据.比如较多的地下震源数量对虚拟炮集成像更有利，但是实际野外采集数据时，震源的数量人为并不可控，在构造稳定的区域往往缺少足够多的地下噪声源，而地下构造相对活跃的区域（比如地下正在开采的矿区等）一般来说更适于被动源体波反射成像.此外，适当增加记录时长通常能增加震源数量，可有效改善最后重构单炮记录的效果.但是记录时长和重构单炮记录的效果并不是呈简单线性关系，更多的震源可能意味着来自地面的震源占据大多数，从而使被动源信号中面波占据主导，即更多的噪声源分布在反射波干涉的非稳相区，使本来能量就很弱的体波被淹没，增加后期面波和体波分离的难度.因此被动源勘探的野外采集部署时应注意采集时间和最终成像效果两者的平衡，目前被动源反射波勘探采集的记录时长一般为几天（刘国峰等，2021）到几十天（Chamarczuk et al., 2022）.

2 被动源数据中弱体波反射信号的提取

由于实际环境噪声场往往主要受控于面波能量，地表附近震源往往无论从数量还是强度上均占主导，导致面波在被动源信号中的能量强（特别是检波器通常在地表布置），相对而言，体波反射则属于弱信号.而想用地震干涉法重建体波，噪声源必须分布在适当的平稳区域，其中大部分需要埋藏在地球内部（Shapiro, 2019）.Vidal等（2014）的数值模拟结果直观地反映了这一点（图8）.如果对原始噪声记录不加甄别和提取，直接对所有记录（包括面波能量占优和体波能量占优）的噪声窗口进行相关求和的虚拟炮集将以面波为主导（图8a），体波反射能量则被淹没其中，而对体波能量进行有效提取并单独重构的虚拟炮集反射则十分清晰（图8b），除了频率较低之外，几乎可与主动源激发的单炮媲美（图8c）.

图8 体波提取前后的地震干涉模拟结果（修改自Vidal et al., 2014）.（a）利用所有记录重构的虚拟炮集；（b）提取体波后的重构虚拟炮集；（c）主动源激发单炮Fig.8 SI result for selection of body wave retrieval (modified from Vidal et al., 2014).(a) Retrieved virtual common-source panel with all recorded data.(b) Retrieved virtual common-source panel after the selection of body wave retrieval.(c) Directly modeled reflection response

因此如何从大量背景噪声连续记录中更好地甄别和提取弱体波信号，就成为被动源反射波勘探成像的核心问题.如果对已重建完成的虚拟单炮记录再进行处理，因为被动源体波与面波在空间位置和频带范围都存在重合，FK等传统滤波方法效果往往并不理想.若是从原始数据出发，在互相关之前就对时窗内的数据先开展面波和体波成分的甄别和分离，再进行相干和叠加，则更有助于提高体波反射信号的信噪比（Chamarczuk et al., 2022; Cheraghi et al., 2015; 刘国峰等, 2021）.目前的研究主要基于被动源数据面波和体波的多方面差异特征，包括信噪比、速度以及方位角等差异来识别和提取弱体波反射信号.

2.1 基于信噪比差异提取被动源体波

被动源野外勘探实际采集的时候，不同时段的噪声主要来源往往有所差异，如白天时段人文交通等面波噪声更为突出集中.人们通过分析野外实际被动源采集数据时发现，面波占优时段与体波占优时段的信号频谱特征往往差异较大.例如图9a、9b分别是面波和体波占优数据段的信号的频谱.面波占优时，频谱呈梳状分布，而体波占优时，频谱形态正常（刘国峰等，2021）.这种差别可通过频率域的信噪比量化表示，见式（11）：

图9 被动源信号不同频率成分的功率密度谱特征（刘国峰等，2021）.（a）面波能量占优；（b）体波能量占优Fig.9 Power density spectrum characteristics for different frequency components of the passive seismic signal (modified from Liu et al., 2021).(a) Energy dominated by surface wave; (b) Energy dominated by body wave

其中，SNR为所述噪声信号的信噪比，|

X

(

f

)|为频率值

f

对应的功率值，频率值

f

＜

f

＜

f

＜

f

＜

f

＜

f

，其中

f

和

f

为功率谱密度计算的最小和最大频率值，

f

和

f

为有效信号的频带范围，

f

和

f

则是包括有效信号和随机噪声的频带范围（图10）.

图10 频率域计算信噪比参考图（张军华等，2009）Fig.10 Spectrum for the signal-to-noise ratio calculation in the frequency domain (Zhang et al., 2009)

刘国峰等（2021）等利用该方法对内蒙古浅覆盖区矿区的数据进行了处理，首先计算获得不同时窗内的噪声信号的信噪比，并基于信噪比差异将噪声信号分选为面波信号和体波信号.如果计算确定某段时窗内属于面波能量为主的信噪比范围，就可以在合成虚拟炮集之前去除这段时窗内的记录，最终可得到体波能量占比更高的虚拟炮集.图11a是未进行面波体波分离处理的单炮，图11b是分离处理后的记录，相比于图11a有较明显反射波能量出现.

图11 频率域分离计算后的虚拟炮集记录（刘国峰等，2021）.（a）以面波为主的虚拟炮集记录；（b）分离面波后突出体波的虚拟炮集记录，红色箭头为反射波信号Fig.11 Virtual shot gathers after frequency domain separation calculation (Liu et al., 2021).(a) Shot gather with surface wave as the main energy; (b) Shot gather with body wave as the main energy after surface wave separation.Red arrows indicate the reflection wave

2.2 基于速度差异提取被动源体波

除了信噪比特征差异外，面波速度一般要明显低于体波，人们也可用速度差异来识别和提取体波.例如Vidal等（2014）采用倾斜叠加照明分析的方法来区分二维测线上的低速面波和高速体波，该研究设计了复杂的地下弹性模型，并基于地表震源和地下震源数据分别构建虚拟炮集，再用不同的射线参数对虚拟炮集记录进行线性扫描叠加.这类似于对数据体进行Tau-

p

变换，其中Tau为截距，

p

为射线参数，然后求取Taup域上的极值对应的

p

值，即可知原始数据窗口内最强能量对应的速度值.地表震源激发的虚拟单炮对应着高绝对

p

值（低速），地下震源激发的虚拟单炮对应着低

p

值（高速），分别对应着面波和体波能量（图12）.

图12 倾斜叠加照明分析示意图（修改自Vidal et al., 2014）.（a）弹性地下模型，黄色倒三角地表检波器，黑点代表地下震源；（b）浅源数据（灰色星）重构的虚拟炮集；（c）虚拟炮集图（b）对应的慢度；（d-e）深源数据（青色星）的虚拟炮集和慢度，绿线表示面波和体波慢度的区分界限；（f）地下深源数据倾斜叠加照明分析的慢度分布图Fig.12 Slant-stack illumination diagnosis (modified from Vidal et al., 2014).(a) Elastic subsurface model with yellow triangles as receivers and black dots as sources.(b) The virtual source retrieved from the shallow source is marked by the grey star.Red lines indicate slownesses.Green lines are pre-defined limits between body and surface wave slownesses in the virtual source function.(c) Slowness representation of the virtual source function from (b).(d) As in (b), but for the source represented by the cyan star in (a).(e) As in (c), but from (d).(f) Illumination diagnosis, consisting of the results for the correlated common source panels from sources in (a) (black dots), with the panel from (b) (grey star) discarded, and the panel from (d) (cyan star) included

Chamarczuk等（2019）把该方法拓展到三维数据，研究分析了地表震源S1和地下震源S2对于5个检波器X1 至 X5组成的十字排列的走时响应，S1和S2震源对于X4、X2、X5组成的inline方向测线时距曲线相同，也就意味着难以从inline方向加以区分能量来自地表还是地下震源.但是两个震源对X1、X2、X3组成的crossline方向测线到时却有很大差别，其中地表S1震源是线性时距曲线，而地下S2震源是双曲时距曲线（图13）.因此不单单扫描inline测线方向的速度，同时也扫描crossline方向的速度，就能更好区分震源是来自于地表面波还是地下体波的能量.该研究表明，实际勘探中若能同时部署crossline测线，将更有助于从速度差异特征上准确识别和提取弱体波信号（Chamarczuk et al., 2019）.

图13 三维观测系统更好区分面波和体波示意图（修改自Chamarczuk et al., 2019）.（a）五个检波器 X1 至 X5 记录的来自地表震源S1和地下震源 S2 的直达走时.两源同时激发，且假定介质传播速度恒定5 km/s.（b）由五个检波器形成的十字排列观测系统.（c）由三条检波线组成的三维观测系统Fig.13 The diagram of event discrimination from surface- or body-waves (modified from Chamarczuk et al., 2019).(a) Direct arrivals from sources S1 and S2 recorded at five receivers X1 to X5 forming a cross-shaped array.We assume a constant propagation velocity of 5 km/s.(b) Configuration of the cross-shaped array formed by the five receivers.(c) 3D survey formed by the three receiver lines

2.3 基于方位角差异被动源体波

被动源实际观测的时候，地下震源往往不是均匀分布的，其中某些特定方位的波动可能贡献了大部分能量，因此可根据方位角差异来区分面波和体波.人们利用波束成形（beamforming）方法（Aki and Richards,1980; Cheraghi et al., 2015; Girard and Shragge, 2019, 2020）可以求取被动源数据中不同频率成分的波的方向和速度，再据此去除干扰面波和分离出有效反射波.计算公式如下（Cheraghi et al., 2015）：

式中，

B

(

v

,θ)为计算出的波束（ θ为传播方向，

v

为传播的平均速度），

f

表示信号的频率，

N

是接收器的数量，∂

f

是频率带宽，

f

是噪声的中心频率，

s

表示第

i

个检波器接收信号的傅里叶变换，其检波器的坐标分别为

x

、

y

.例如图14是加拿大一个CO存储监测项目实测数据的计算波束功率，圆的方位角表示波的来源方向，圆的半径代表不同的波速，振幅的大小表示该点（每个点都对应着一个速度和方向）的归一化波束功率值.该例子中，工区附近240°和 330°方位角之间的发电厂是环境噪声的主要来源，而且速度能量明显大于3 000 m/s，这部分能量以直达波、折射波、反射波等体波能量为主，有助于后续的被动源反射体波信号提取（Cheraghi et al., 2015）.

图14 波束成形的方向性分析（修改自Cheraghi et al.,2017）.正北方向为0°，角度沿逆时针方向增加Fig.14 Directional Beam-forming analysis (modified from Cheraghi et al., 2017).The north had an azimuth of 0°,which increased to the west (i.e., counterclockwise)

人们还可利用波束分析方法开展更精细研究，比如Cheraghi等（2015）对加拿大Lalor矿区的被动源数据细分成不同频率和时段来研究波动的主要方向和速度.首先将数据按小时段划分，并对各时段数据分别进行1～3 Hz、10～13 Hz 和 23～25 Hz三个频段的波束分析.图15表示研究区一天之中每个时段对应的最大波束功率谱及其对应的方位角.可见，三个频段的最大方位能量主要集中在两个方向，低频段能量主要来自于300°方位角（东北方向8 km）的人类活动，高频段能量主要来自于约80°方位角（正西方向）的地下采矿活动，而10～13 Hz能量则来自于这两个方向的综合贡献.基于波束分析的方位角信息，并结合速度差异特征（大于3.2 km/s的主要为体波能量），可以有效识别和提取弱体波信号，再对虚拟炮集数据进行处理获得最终的反射叠加剖面（图16b），与主动源剖面（图16a）对比可知，尽管被动源剖面的分辨率和精度要低，但该矿区的三套主要反射I1、J2、K3还是清晰可见.

图15 不同小时段的最大波束功率及其视速度和方位角（修改自Cheraghi et al., 2015）Fig.15 Calculated beamformers for the maximum beam power of each hourly panel and its apparent velocity and azimuth (modified from Cheraghi et al., 2015)

图16 加拿大Lalor矿区R141线的倾角时差校正叠加剖面（修改自Cheraghi et al., 2015）.（a）主动源剖面；（b）被动源剖面Fig.16 DMO stacked section along line R141 in Lalor deposit, Canada (modified from Cheraghi et al., 2015).(a) Active-source 3D survey and (b) passive interferometry survey

3 被动源反射地震数据处理和成像

3.1 被动源原始数据预处理

不同于天然地震学研究领域常用的宽频带地震仪（如60 s、120 s长周期地震仪），被动源地震勘探领域目前主要采用成本更为低廉、布设更为便捷、间距更为密集的节点式地震仪（如5 Hz、10 Hz检波器），旨在记录更多高频的体波反射信号.因此在预处理阶段，比如带通滤波等处理，需根据检波器的自然频率以及研究目标的特征选择合适的频带范围.由于被动源数据往往信噪比较低，尤其是连续记录中的一些强振幅噪声会在虚拟炮集干涉计算过程中产生一些虚假振幅，因此需要在数据预处理阶段尽可能加以去除.压制强振幅噪声一般是在对原始数据开展去均值、去趋势、端点归零等基本处理之后.其中一种较简便方法是对数据进行归一化处理，如Oren和Nowack（2017）在时间域使用sign-bit信号归一化来降低噪声权重，但该方法无法保留原始数据集的动态范围；另一种常见策略是利用窗口进行数据选择，即在互相关之前识别并去除（或降低权重）强噪声的窗口（Panea et al.,2014），该方法需要对所有数据窗口进行判断和选择，工作量比较大.最近几年，Girard和 Shragge（2020）、方捷等（2022）等开发自动识别和处理流程来提高效率，总体的处理原则是保持均方根振幅变化的趋势.比如计算每个滑动时窗内的振幅平均值，若该窗内的样点值偏离均值较大，就用均值代替当前的样点值，或者直接定义振幅的一个阈值来将高振幅记录降低到阈值以下.

除了连续强振幅噪声，还有一些在时间域以孤立的尖峰或毛刺形式出现的噪声，Girard和Shragge（2019）和方捷等（2022）使用基于L1迭代加权的最小二乘方法（Claerbout, 2014）处理.它采用拟合的思想对原始数据进行处理，可对高振幅信号进行压制，而低振幅信号不受影响.从图17可见，经过处理后强振幅噪声得以压制，而数据的动态范围仍然能够保持.

图17 时间域强振幅控制示意图（方捷等，2022）Fig.17 Strong amplitude control in the time domain (Fang et al., 2022)

在时间域去除强振幅噪声基础之上，Girard 和Shragge（2019）又在频率域压制检波点的固定频率噪声（如机电噪声等）.该方法首先对数据进行傅里叶变换并计算信号的振幅谱

M

和相位 ∅，见式（13）和（14），其中，R和I分别表示信号傅里叶变换的实部和虚部.然后与时间域处理类似，定义一个阈值将高谱幅值降低到阈值以下得到振幅谱

M

.最后根据公式（15）使用

M

重建信号：

式中，

F

是傅里叶反变换算子.图18是频率域强振幅控制的示意图，可见经过频率域处理后强振幅噪声得以进一步压制.

图18 频率域强振幅控制示意图（修改自Girard and Shragge, 2019）.（a）频率域强振幅控制前（蓝色）和后（红色）的时域信号；（b）频率域强振幅控制前（蓝色）和后（红色）的部分频谱Fig.18 Strong amplitude control in the frequency domain(modified from Girard and Shragge, 2019).(a) Timedomain trace before (blue) and after (red) debursting.(b) Partial spectrum before (blue) and after (red)debursting

Girard和Shragge（2019）展示了预处理不同阶段的虚拟道集（图19），分别是基于原始数据、窗口筛选后数据，以及强振幅压制（包括时间和频率域）后的数据.从图19共偏移距排列的虚拟道集可见，完成预处理后的数据噪声压制明显，其中0.5 s左右的反射同相轴的信噪比相对原始数据得以大幅提高.

图19 互相关+叠加得到的以共偏移距（1.0 km）排列的虚拟道集（修改自Girard and Shragge, 2019）.（a）原始数据；（b）窗口选择数据；（c）强振幅压制数据Fig.19 Single extracted cross-correlation + stack trace at a virtual source-receiver offset of 1.0 km (modified from Girard and Shragge, 2019).(a) Raw data.(b) After data selection.(c) After time and frequency debursting

3.2 被动源虚拟炮集的数据处理

经过预处理和面波、体波分离后构建的虚拟单炮尽管反射依然较弱，但可借鉴主动源反射地震数据的处理思路，利用多次叠加技术增强反射能量.被动源数据处理的主要目的是尽可能提高炮集反射的信噪比，改善最终成像质量，一般包括静校正、振幅垂向补偿和横向均衡、随机和相干噪声压制、多次波压制、速度分析、切除、叠加和反褶积等步骤.

被动源数据的静校正一般采用高程校正为主，将起伏地表虚拟炮集校正为水平地表虚拟炮集.值得注意的是，如果直接在被动源原始数据上先应用静校正再相关合成虚拟炮集，会导致静校正量发生改变，需对反射记录非因果部分进行两倍炮点校正，因果部分进行两倍检波点校正才能叠加生成正确的虚拟炮集（郑钧等，2022）.一般来说，通过被动源数据预处理过程中的去均值、去趋势已较好完成单炮能量的重新均衡.在此基础上，还可采用球面扩散补偿和地表一致性振幅补偿进行校正.

在虚拟炮集静校正和振幅均衡基础上可更好地开展噪声压制工作，其中随机噪声往往是被动源数据中面波之外最严重的噪声，我们可通过原始数据筛选、去均值、去趋势、强振幅控制、频率滤波等处理进行较好压制.除了随机噪声，虚拟炮集记录中的错误旅行时、干涉噪声还有多次波的存在都会大大影响后续地震处理和成像的精度.针对这些问题，Hu等（2020）和Wang等（2021）提出了一种迭代去噪方法，将聚焦变换引入到三维L1范数稀疏约束反演Closed-loop SRME（SRME：自由表面多次波压制技术）算法框架，实现被动源虚拟炮集的相干噪声和多次波一体化压制.首先，选取估计的一次波构建聚焦变换的反算子，使聚焦变换具有多次波降阶的作用，一次波变换为聚焦点，高阶多次波变换为低阶多次波，利用多次波降阶的深层信息纠正，重构被动源虚炮集中错误的远偏移距走时信息.其次，由于聚焦变换能够将有效信号聚焦到聚焦点，使有效信号与噪声的分布区域产生明显的区别，通过多次迭代的方式构建聚焦域切除算子逐步压制噪声，最终获得具有正确远偏移距走时信息与高信噪比的一次波虚拟炮集.该方法克服了虚拟炮集中相干噪声的限制，大大提高了噪声型被动源数据处理与解释的准确度.

Wang等（2021）通过对模拟数据试验，验证了基于聚焦变换和SRME三维稀疏反演方法对于被动源干涉噪声和多次波压制的有效性（图20）.图20a是用来试验的一个五层模型，该模型内部存在一个断层和凹陷式地层.299个随机震源布设于地下1.2～1.5 km的位置（图20b），用来模拟天然背景噪声（图20c）.图20d为互相关生成的虚拟炮集，图20e是利用该迭代去噪方法后的虚拟炮集，图20f则展示用来对比的SRME处理后的主动源记录.对比可知，经过处理后的虚拟炮集走时信息和同相轴曲率获得了纠正（白色箭头），且干涉噪声和多次波去除效果明显（黑色箭头）.经过多种方法对不同噪声压制之后，还需对被动源炮集记录两侧的虚假反射以及低速面波残留能量进行切除和叠加来进一步提高信噪比.

图20 基于聚焦变换和SRME三维稀疏反演的被动源数据去噪试验（修改自Wang et al., 2021）.（a）数值试验速度模型；（b）随机震源采用的子波；（c）震源分布位置；（d）互相关生成的被动源虚拟炮集；（e）应用该方法处理的虚拟炮集；（f）SRME处理后的主动源炮集Fig.20 Numerical tests of noise removal via focal-denoising closed-loop SRME based on the 3D sparse inversion of passive data(modified from Wang et al., 2021).(a) Model for simulating noise-source data and (b) wavelets of the random sources.(c) Source location.Virtual shot gathers from the noise-source data using (d) cross-correlation and (e) the proposed method.(f) Active-source data using SRME

此外，速度分析对最终叠加或偏移成像也十分重要，一般要对速度谱进行多轮、加密拾取.由于相关合成虚拟炮集过程中，可能会产生不同频率的干扰波，在生成速度谱时利用带通滤波只保留频率能量最集中的部分，这样有利于速度谱能量团的拾取.图21是一个实际被动源资料的速度分析例子.速度谱中有明显的能量团，共中心点道集有相对应的反射同相轴，而局部叠加结果中也可见连续性的反射（刘国峰等，2021）.

图21 被动源反射信号速度分析和局部叠加（刘国峰等，2021）Fig.21 Reflection velocity analysis and stacked local passive data (Liu et al., 2021)

3.3 被动源数据直接偏移成像

目前被动源成像一般采用间接偏移的方法，即首先将被动源数据重构为主动源虚拟炮集，再进一步采用主动源地震的处理和偏移方法来成像.但大部分地震干涉法重构的炮集仍然面临着数据量大、信噪比低以及假设条件过多等情况.另外一种思路是Artman（2006）最早提出的被动源直接偏移成像方法，他采用噪声源数据本身作为边界条件，在频率域直接利用被动源传动波场数据采用单程波偏移，得到了地下地质体构造图像.直接偏移成像方法免去了重构虚拟炮集时产生的二次噪声，提高了成像的准确性，不过该方法对于速度模型的精度要求比较高.

被动源直接偏移方法本质上是利用了多次波信息进行偏移成像，人们又将该方法扩展到逆时偏移成像上，从而适应更复杂的速度模型和高陡地下构造.Zheng等（2015）将逆时偏移结合多次波成像的边界条件应用于被动源地震数据，直接利用数据本身正传反传并进行偏移成像.Cai 等（2019）提出了被动源逆时偏移能量归一化成像条件，并采用波场分离的方法较好消除逆时偏移的低频噪声，同时利用并行分段叠加的成像策略提高计算效率.图22的数值试验表明该方法能较好克服被动源震源分布不均匀导致的成像能量比例失调问题，改善了模型深部与边缘成像的质量.

图22 被动源地震数据直接偏移的数值试验（修改自Cai et al., 2019）.（a）Marmousi模型，震源分布在左下方的黑框红框；（b）随机截取的一段4 s的被动源记录；（c）基于全波场分离归一化的被动源数据直接逆时偏移成像Fig.22 Numerical tests for the direct migration of passive data (modified from Cai et al., 2019).(a) Marmousi model with the passive data generated from source locations marked by black and red boxes.(b) Example of a passive seismic record with a length of 4 s.(c) The direct migration image by full wavefields decomposition normalized reverse time migration

Girard和 Shragge（2020）利用加拿大Lalor 湖附近矿区的被动源实际数据开展三维直接偏移成像研究.图23展示了主、被动源成像结果在矿体中心位置的[

E

,

N

,

Z

] = [2.36, 2.54, 0.8] km提取的水平切片和两个垂直深度剖面.可见，在0.5 km和1.5 km深度范围内，被动源三维数据体揭示工区西部是NE倾向反射，而东部则发育西倾构造，与主动源数据体的解释基本一致.被动源直接成像因为无需重构炮集，因此在诸如随钻测井、矿区现场等处理时效性要求高，背景速度又比较确定的研究场景具有较大应用潜力.

图23 加拿大Lalor湖附近矿区主、被动源水平切片和垂直剖面（修改自Girard and Shragge, 2020）.（a）被动源数据的直接偏移成像；（b）主动源数据的偏移结果Fig.23 Comparison of the horizontal slice and vertical section of 3D active and passive surveys in Lalor deposit, Canada (modified from Girard and Shragge, 2020).(a) Passive direct migration stack and (b) the active-source migration volume

4 被动源反射波勘探应用实例

随着便携节点地震仪等采集设备的日益普及，以及被动源反射地震处理成像方法的发展，被动源反射波勘探逐渐从理论探索走向实际应用.目前人们在开展被动源反射勘探实际研究时，往往还进行重合测线的主动源采集，这样可通过主、被动源成像剖面的对比，进一步分析被动源反射方法的成像精度，以及探讨其在不同勘探领域的应用潜力.

4.1 CO2封存场地监测研究实例

Cheraghi等（2017）利用被动源地震数据对加拿大Saskatchewan地区的一个CO封存场地开展研究，分析被动源地下结构成像的精度，以及被动源方法用于时移地震监测的潜力.该项目沿两条垂直检波线共采集两个时间段（2014年6月的23天、2015年2月的13天）的被动源数据.图24展示了L1线23天被动源数据的剖面结果，其在0.8～1.5 s目的层范围内成像精度较高，不仅在一定程度上可与带通滤波后的主动源地震剖面相比较，甚至与注入井的合成记录的相关系数也达到0.5～0.65.特别是一些关键层位，比如下科罗拉多页岩（LC）（0.85 s，约 930 m）、Vanguard 组砂岩（VF）（1.1 s，约 1 100 m）和 Bakken 页岩（BK）（1.45～1.5 s，约2 100 m）等标志层的反射响应均很明显.不过受限于信噪比，被动源剖面无法识别更深部的前寒武地层（PC）.

Roach等（2015）研究表明CO地下储层的监测研究需要重复度很高的时移地震数据，即不同时段采集的地震振幅变化要控制在10%～20%，这样才能刻画CO储层性质随着时间推移的变化（预期波阻抗变化范围约为-8%～-17%）.尽管该项目的被动源成像对关键层位的结构信息刻画比较清楚（图24），但是其精度还无法满足时移地震监测的需求，因为不同时段的成像结果对比表明，23天采集数据获得的剖面与时隔8月后的13天剖面成像差异很大（图25），而这显然不可能是CO储层性质变化所致.由于两次被动源成像采用的检波器和处理流程都相同，推测成像的巨大差异可能由于不同的采集时长、不同噪声强度以及噪声功率谱随时间的变化（Mehta et al., 2008）等原因引起.该研究表明目前被动源反射还只能达到构造成像的精度，尚很难通过被动源的时移反射地震来研究地下储层的性质变化，但Cheraghi等（2017）仍认为他们的研究可以通过使用更长持续时间的噪声记录来加以改进，且被动源地震由于其低成本优势，如果能开展不间断的连续记录，也许能在地下结构的长期动态监测方面更有作为.

图24 加拿大Saskatchewan的CO2封存场所的被动源研究（修改自Cheraghi et al., 2017）.（a）被动源剖面（红色为注入井的合成记录）；（b）井的合成记录（子波从被动源数据提取）；（c）主动源数据剖面；（d）深时转换的P波速度测井记录.LC：下科罗拉多页岩；VF：Vanguard组砂岩；BK：Bakken页岩；PC：前寒武纪岩石Fig.24 Passive surveys of aquistore CO2 storage site Saskatchewan, Canada.(modified from Cheraghi et al., 2017).(a) Passive section (red curve is the log-based synthetic seismogram).(b) Synthetic seismogram generated from a wavelet extracted using the passive data convolved with borehole log-based reflectivity.(c) Active-source cube.(d) Depth-to-time converted P-wave velocity log.LC, VF, BK, and PC represent the Lower Colorado shale, Vanguard Formation, Bakken shale, and Precambrian,respectively

图25 CO2封存场所的L1和L2线不同时段的被动源地震成像（修改自Cheraghi et al., 2017）Fig.25 Passive seismic imaging for different time periods at the CO2 storage site along lines L1 and L2 (modified from Cheraghi et al., 2017)

4.2 金属矿区应用实例

金属矿床多位于结晶岩地区和构造发育带，地下地质情况复杂，且金属矿体形态复杂，非均质性极强，一直是反射地震勘探的难点，相较于油气勘探领域，金属矿的主动源地震技术进展较为缓慢（汤聪等，2022）.不过对于被动源反射波成像来说，金属矿集区特别是那些正在开采的矿区因为能够提供十分关键的地下噪声源，反而是比较适合被动源勘探应用的潜在领域.

Chamarczuk等（2022）最近在芬兰东部的一个结晶岩多金属（Cu-Co-Zn-Ni-Ag-Au）矿床Kylylahti开展了一项被动源三维地震成像工作.该研究在3.5 km× 3 km的区域内以200 m线距和50 m点距总共布置了994个垂直检波器，并持续记录30天的被动源噪声数据.数据采集期间日常采矿活动持续进行，尤其是每天100～800 m深处的雷管爆破提供了大量地下噪声源.图26是基于所采数据构建的三维地震虚拟炮集与主动源数据的对比（均显示其中6条检波线）.由于矿区地下构造异常复杂，即使在主动源数据的单炮中也较难识别有效反射信号，好在开展的是三维地震采集，相邻检波线的数据可用来辅助追踪有效反射.总体上，主、被动源单炮的信噪比略有差异外，主动源单炮上的反射信号在被动源虚拟炮集也有所指示（图26绿色箭头），但在大偏移距的被动源炮集上则存在着一些假象（图26黑色箭头）.

图26 主、被动源三维共炮集数据对比（修改自Chamarczuk et al., 2022）.（a）主动源炮集；（b）30天记录重构的被动源虚拟炮集；（c）10天记录重构的被动源虚拟炮.VSG：虚拟炮集Fig.26 Comparison of 3D common-source gathers using active and passive data (modified from Chamarczuk et al., 2022).(a) Activeshot gather, (b) VSGs obtained using 30 d of noise, and (c) VSGs obtained using 10 d of noise.VSG:Virture Shot Gathers

基于重构虚拟炮集开展后续的三维反射地震数据处理，获得了最终的叠后偏移剖面（图27）.对比发现，主、被动源的剖面具有很好一致性，并且与该区基于钻孔的地质解释成果较为吻合，比如图27中多处红色箭头所指的岩体边界均能看到较明显反射能量.特别要指出的是，被动源剖面在该矿区有效弥补了主动源地震由于地形限制无法采集而导致的浅部数据缺口部分（图27蓝框部分），其中箭头4所示的位置清楚揭示了岩体在浅部的走向和边界.该研究证明了被动源反射地震勘探方法在金属矿地区有较大的应用潜力，因为最终成像不仅与基于钻孔的地质解释吻合度较高，甚至还提供了主动源地震成像缺失的一些关键信息.

图27 芬兰Kylylahti结晶岩矿集区主、被动源叠后偏移剖面（修改自Chamarczuk et al., 2022）.（a） 主动源剖面；（b）主动源剖面叠合地质解释；（c）被动源剖面叠合地质解释KAL沉积（蓝色）；OUM岩体（绿色）；S/MS矿化半矿化岩体（红色）Fig.27 Comparison of post-stack migrated sections obtained from active and passive surveys in Kylylahti polymetallic mine, Finland (modified from Chamarczuk et al., 2022).(a) Active survey; (b) Active survey with the geological model background;(c) Passive survey with the geological model background.Geology in the background is color-coded as follows: Kaleva Sedimentary Belt (KAL): blue; semi-massive to massive sulfide (S/MS) mineralization: red; Outokumpu ultra-mafics (OUM):green

4.3 煤矿地区应用实例

除了金属矿勘探方面的应用之外，最近Ning等（2021）还将被动源反射成像技术应用在某煤矿地区的地下结构成像.该研究在600 m× 800 m 的区域使用了202个检波器，间距为50～100 m，连续观测记录27天，期间煤矿开采活动正常进行.值得说明的是，该研究只选取其中一天的数据以两种不同方式开展试验，一种是首先应用波形偏移处理和tomoDD 地震定位方法 （Miao et al., 2021; Zhang and Thurber, 2003）从连续波形记录中检测到 576个高信噪比的地震事件，然后只利用包含这些高信噪比地震事件的记录进行虚拟炮集重构；另一种则是将连续地震记录以30 min时窗分段再互相关进行虚拟炮集重构.该研究采用Chamarczuk等（2021）的方法对两种方法重构的虚拟炮集进行相似度定量分析，其中相似度定义的范围在 -1～1之间变化，值越大说明两者越相似.实际计算所得的两组虚拟炮集的相似度高于0.7，说明两组炮集总体非常相近，也从侧面反映该研究提取的体波信号可信度较高.

基于提取体波信号重构虚拟炮集，并进行后续处理最终获得研究区的柯希霍夫时间偏移剖面（图28）.两条不同检波线的时间偏移成果均有两个明显的反射层（黑色虚线表示），对应着区域的煤层构造.其中较浅的煤层（箭头“1”表示）连续性较差，较深的煤层（箭头“2”表示）连续性相对较好，反映了该区较浅的煤层开采程度较高的现状.另外工区自西向东，煤层呈现出东倾的趋势，奥陶系地层处于第2个煤层之下，与该区已有的地质认识和工作面采掘资料吻合良好.该项研究表明正在开采中的煤矿能产生丰富的地下噪声源，这非常有助于体波信号的提取，被动源地震反射成像方法可以用来辅助调查区域内的煤层构造赋存形态和采空区范围.实际工作中还可进一步将被动源反射地震剖面与面波成像的速度反演结果等进行联合解释以提高研究精度.

图28 不同检波线的偏移剖面（修改自Ning et al., 2021）.黑色虚线表明煤层成像位置Fig.28 Migration sections of receiver lines (modified from Ning et al., 2021).The black dashed lines show the coal seams imaged using ambient-noise data

5 结论与展望

本文主要介绍了被动源反射波勘探成像的一些关键技术，包括地震干涉法重构虚拟炮集、基于多种差异特征分离面波和体波以及被动源原始数据预处理、虚拟炮集多次波压制和直接偏移成像等方法.如何在面波能量所主导的大量噪声连续记录中高效率地甄别和提取弱体波信号，是被动源反射地震成像的核心问题.通常需要在干涉和叠加之前对原始数据不同时窗内的能量加以判断和选择，数据存储和计算量都十分巨大.为了提高效率和精度，人们已开发各种自动化流程加以处理，同时大数据处理中的图像识别、机器学习、人工智能等技术已经被用于并将进一步在被动源反射地震勘探成像中发挥更重要作用（Chamarczuk et al., 2019）.

随着勘探技术的不断发展，一些新方法也被应用于被动源反射波成像.比如将地震超越干涉法的Marchenko成像应用到被动源数据，该方法能够直接面向目标成像，避免检索炮集中自由表面多次波和内部多次波的影响（Jin, 2020）.同时，以全波形反演为代表的高分辨率反演技术也被引进到被动源成像之中，如利用随机干涉重构（尚旭佳等，2021）和多尺度震源编码反演（蔡中正，2021）等技术获得较为准确的地下介质速度模型，再用该模型对被动源数据利用逆时偏移方法进行直接成像.而Shiraishi和Watanabe（2022）最近将被动源声波逆时偏移推广到弹性波逆时偏移，并获得了P波和S波的被动源成像.此外，主、被动源数据的联合反演和成像技术也在不断发展，比如可以利用被动源数据丰富的低频信号来降低高速屏蔽层对地震信号的散射和屏蔽作用，通过主、被动源据的联合应用，使其在照明和频带上优势互补，能显著提高成像和反演质量（张盼等，2019）.

除了成像方面研究之外，被动源数据的采集也值得关注.噪声源的分布是影响干涉重构效果的关键因素，如有可能尽量选择那些地下噪声源丰富的区域开展被动源反射成像研究，比如正在开采的金属矿区和煤田等.同时被动源勘探进行数据采集部署时应特别注意人工噪声源的分布角度，条件允许时应尽量部署十字排列或者三维观测系统.相比单纯的二维采集，多出的另一个方位信息有助于更好区分记录到的能量主要来自地表还是地下（Chamarczuk et al., 2019）.而在特定目标区进行主、被动源的混合采集也颇有前景，比如对构造稳定、缺少足够多地下噪声源的地区进行被动源勘探时，补充少量主动源数据就可以有效控制和改善被动源的成像效果.同时在稀疏炮点的主动源勘探中，有效利用被动源的信息也能够在成像中增加更多的细节信息，提高成像质量（张盼等，2015）.

被动源反射地震勘探，最近在地下矿产勘查等结构成像领域已初步取得实效，而在碳封存等需要动态监测储层性质变化的勘探领域仍在不断探索之中.尽管目前被动源反射成像与主动源成像的效果尚存在明显差距，但是作为一种低成本的绿色环保方法，且被动源勘探的部署和应用方式非常灵活，既可作为主动源地震的补充，又能结合已有的被动源面波速度成果进行联合分析，因此该技术的发展一直持续受到人们关注.特别是近年来由于节点检波器、分布式光纤等技术的日益普及，连续接收被动源地震记录成本已大幅度降低，以及高性能计算的大规模兴起，使得海量数据的存储和计算成本快速下降，被动源勘探在成本方面的优势得以进一步凸显，也在一定程度上刺激着相关处理新技术和成像新方法的不断发展，人们应该可以期待被动源反射地震勘探在未来会有更多更好的实际应用.

致谢

感谢三位评审专家提出的修改完善建议，感谢编辑部的大力支持.