2020, Vol. 63
, CHEN ZhanGuo, LEI ChaoYang, ZHANG YaHong, ZHANG Jie, SUN ZhenTao, ZHANG WeiHong, WU YongShuan
随着分布式光纤声学传感(DAS)技术(Mateeva et al., 2012)的不断进步和完善，3D DAS-VSP采集(Mateeva et al., 2013)已成为石油工业界的研究热点之一.光缆能永久安装于油水井或压裂作业井中，抗电磁干扰、耐高温，实现全井段采集或4D监测，高效低廉.2012年，Shell在墨西哥湾Mars油田进行海底节点(OBN)采集的同时，采用QinetiQ公司的OptaSense问询器和已经安装的分布式温度传感系统(DTS)进行海上3D DAS-VSP的首次采集试验(Mateeva et al., 2013a)，并探索将之作为一种独立的地震工作方法和技术，应用于构造的高精度成像和4D油气藏动态/永久监测(Mateeva et al., 2014, 2015；Wu et al., 2017；Zwartjes et al., 2018).东方物探公司与成都电子科技大学于2019年联合发布成功研制了“u-DAS分布式光纤传感地震仪”，并完成20余次现场试验.
相比于地面地震，VSP是在井中接收地震波，由于地震波仅经过近地表一次，其分辨率较高，且具有很多优点(朱光明，1988；刘宝和，2008)，主要应用于井周地层标定和各向异性等研究.但受井中检波器串级数局限，VSP采集较为困难、成本高，反射波成像范围窄，也不能对井中最浅接收点上方区域成像；同时随着离开井距离的增加，叠加次数减少，成像精度逐渐变低，导致VSP成像同相轴难以追踪，不能有效地应用于构造成像和4D油藏监测.
为了显著扩大VSP成像范围，需要借助井中接收器记录到的多次波进行成像，以便准确连续地追踪成像同相轴、提取地震属性，提高地质解释可靠性.然而，由于多次波旅行路径特殊，至少碰触地下反射面一次，多次波直接成像尚有诸多困难.幸运的是，近十几年蓬勃发展的地震干涉法，能为多次波成像提供另一种选择，已逐步发展成为一种称之为“地震干涉成像”的方法技术，尤其适用于VSP多次波成像(Jiang et al., 2005；Yu and Schuster, 2006；He et al., 2007).已有研究结果(Wu et al., 2015；Zhan et al., 2015)表明：3D DAS-VSP地震干涉成像范围和分辨率明显优于3D DAS-VSP反射波成像，与3D OBN资料的反射波成像相当，显著好于3D VSP(井中检波器串)反射波成像.
地震干涉法的本质是构建虚震源(Schuster，2009)，即两个接收点记录的互相关，就像是在一个接收点处激发另一个接收点处接收的记录一样.虚震源构建过程类似于基于模型的基准面重建，但本质上不同(Schuster and Zhou, 2006)：它是实测数据驱动的，无需速度模型，也无需震源位置和激发时刻信息.
地震干涉法起源于Claerbout(1968)提出的方法，即通过自相关地震台站所记录的1D透射波响应，重建其反射波响应.猜想这一方法可推广应用于3D介质，并将此方法命名为“声波日光成像”(Rickett and Claerbout, 1999).Schuster(2009)和Snieder(2002)利用稳相位分析，证实了地震干涉的可行性，分别应用于VSP、地面地震、单井地震和井间地震之间的采集基准面重建(Schuster，2009)以及地震尾波的面波恢复及其横波速度反演(Snieder，2004)，极大地推动了地震干涉法的发展.Bakulin和Calvert(2006)基于逆时声学，从物理直觉上提出虚震源方法，主要应用于高陡构造(如盐丘侧翼)的成像和4D油藏监测(Bakulin and Calvert, 2008)，避免上覆复杂介质的透射影响.
自2002年始，Wapenaar发表了关于地震干涉理论的系列论文，从无损介质至衰减介质，根据相关型(Wapenaar，2006)或褶积型(Wapenaar，2010)互换定理，严格证明了格林函数恢复表示，统一了地震干涉的理论框架(Wapemaar，2010a)，既可应用于被动源(随机波场)地震，也可应用于主动源(确定性波场)地震(Wapenaar et al., 2010b, 2010c, 2011).被动源地震干涉利用海洋潮汐等长周期自然界环境噪声或交通等人类文明噪声，通过两个接收点记录相关并对记录时间段求和，恢复这两个接收点间的Rayleigh型面波格林函数，进而利用提取的频散曲线反演横波速度(Campillo and Paul, 2003；Shapiro et al., 2005；房立华等，2009；徐佩芬等，2009；陈国金等, 2013, 2017；Wang et al., 2019)，或恢复反射波用来成像等(Draganov et al., 2007, 2009；齐诚等，2007).主动源地震干涉则是基于实测数据，通过两个接收点记录相关并对震源求和，重建这两个接收点间的地滚波或反射波地震响应，主要应用于地震勘探中的地滚波消除(Forghani et al., 2010；Halliday et al., 2010；Chiffot et al., 2017)，或地质构造成像和4D油藏监测等(Bakulin and Calvert, 2006；Wu et al., 2015).
本文从Wapenaar和Fokkema(2006)提出的相关型格林函数恢复出发，利用VSP上下行波场分离技术，研发了VSP地震干涉成像方法技术，可应用于常规VSP或3D DAS-VSP采集的资料.通过楔形体和水平排列的不同尺度块体模型，进行了成像分辨率测试.利用新疆地区某个井的声波测井及零偏移距VSP资料，设计了一个含有不同深度不同尺度溶洞的水平层状模型，讨论了观测井段接收器道数和位置对干涉成像的影响.最后，将该方法应用于新疆地区某个井采集的3D VSP资料，并与地面地震和VSP反射波的成像结果比较.
1 方法原理VSP地震干涉成像方法是一种对多次波间接成像的方法，它基于近年来发展起来的地震干涉理论，将VSP地震资料重建为拟地面地震资料，然后偏移成像.而地震干涉法中有多种格林函数恢复算法可供选择，如互相关(Wapenaar and Fokkema, 2006)、反褶积(Vasconcelos，2008)和多维反褶积(Wapenaar，2011)等.它们各具特点(Wapenaar et al，2010a, 2010b)：互相关算法简单且稳健，但仅适用于无损介质，且恢复的地震响应不够精确；反褶积和多维反褶积算法不仅适用于有损介质，且恢复的地震响应准确，但反褶积算法不够稳健且仅适用于1D介质，而多维反褶积算法则需要进行反演，计算量巨大.
VSP采集通常是三分量观测，虽然可以利用弹性波地震干涉来恢复拟地面地震弹性动力学格林函数，但由于井中接收点数有限(通常是16个3C检波器串)以及为了简化后续的拟地面地震资料偏移成像问题，实际上往往只进行声波格林函数恢复，因此，本文采用Wapenaar和Fokkema(2006)提出的声波地震干涉方法.声波地震干涉方法是根据Rayleigh互换定理(De Hoop，1988；Fokkema and Van Den Berg，1993)和无损介质的逆时不变性原理(Bojarski，1983；Fink，1997)推导而来.如图 1所示，假设在自由表面∂D0和任何形状边界∂D1组成的封闭界面∂D的体V内介质是无损的，并假设封闭界面∂D是一个具有很大半径的球面以及∂D上和外部的介质是均匀的，那么，在高频近似条件下的频率域声波格林函数恢复关系式，可表示为
|
(1) |
|
图 1 声波地震干涉示意图.图中射线表示完全的地震响应，包括直达波和∂D1内部和外部不均匀性引起的多次散射以及自由表面∂D0的反射.引自Wapenaar和Fokkema(2006)的图 2 Fig. 1 Schematic diagram of acoustic seismic interferometry. The rays in the figure represent a complete seismic response, including direct waves and multiple scattered waves caused by ∂D1 internal and external heterogeneity, and reflections from the free surface ∂D0. From Figure 2 of Wapenaar and Fokkema (2006) |
式中，
方程(1)说明：在多个震源均匀分布于围绕两个接收点的封闭震源边界上的情况下，互相关两个接收点的记录并对所有震源求和，就能获得在其中一个接收点处激发另一个接收点处接收的地震响应及其逆时地震响应的叠加.那么，将方程(1)应用于VSP地震资料，可将VSP地表震源重建于井中接收点处，以获得在井中激发接收的拟单井地震资料，但这不是本文的目的.
根据炮-检互换原理，将VSP震源和接收点位置互换，得到逆VSP(IVSP)资料，然后进行虚震源构建，即利用方程(1)来获得在地表激发接收的拟地面地震资料(Schuster，2009).实际上，也可以将VSP共炮点道集转换为井中共接收点道集，这等价于IVSP共炮点道集，然后对同一个井中接收点记录的来自两个不同地表震源的地震波场互相关并对接收点位置积分(如图 2所示)，来获得在地表一个震源xa处激发(虚震源)另一个震源xb处接收的拟地面地震响应，即将井中接收点位置zr重建至地表的一个震源xb位置，那么，用实际的地震波场替代格林函数，并将方程(1)改写为时间域，即
|
(2) |
|
图 2 利用VSP资料重建拟地面地震资料示意.红色和黑色虚线表示重建的拟地面地震镜像反射波路径；蓝色实线表示自由表面鬼反射波地下反射点位置；绿色实线表示VSP反射点位置 Fig. 2 Schematic of pseudo-surface seismic data reconstruction using VSP data. The red and black dashed lines represent the reconstructed pseudo-surface seismic mirror reflection wave path; the solid blue line represents the position of the underground reflection point of the free surface ghost reflection wave; the solid green line represents the position of the VSP reflection point |
式中，G(xb, xa, t)表示在xa处激发(虚源)xb处接收的重建格林函数，对应于方程(2)右边互相关的正延迟部分；G(xb, xa, -t)则表示相应的逆时格林函数，对应于方程(2)右边互相关的负延迟部分；S(t)表示在xa处的虚震源时间函数；C表示比例因子；*表示褶积；u(zr, xb, t)表示在地表xb处震源激发的井中某个接收点zr处接收的波场，u(zr, xa, -t)表示在地表xa处震源激发的井中某个接收点zr处接收的逆时波场；L表示观测井中接收点zr组成的轨迹.
地震勘探中的有效波是反射波，因此，对于重建的拟地面地震反射波而言，根据稳相位分析，方程(2)右边干涉积分的主值来源于积分面上的平稳点(Snieder，2004; Snieder et al., 2006；Schuster，2009).据此，如图 2所示，当井中接收点zr是其两个地表震源xa和xb的稳相位位置，即地表震源xa至井中接收点zr的下行直达波路径与地表震源xb至井中接收点zr的下行鬼反射波路径中的最后一段重合时，方程(2)右边干涉积分对重建的拟地面地震反射波贡献最大；换言之，由于地震波场互相关运算在数学上表示旅行时间相减，即地表震源xb激发井中接收点zr接收的自由表面一阶下行鬼反射波至时间，减去地表震源xa激发井中接收点zr接收的下行直达波至时间，就得到方程(2)左边的地表xa激发(虚震源)xb接收的镜像反射波.显然，与VSP反射波成像相比，VSP干涉成像的成像范围(VSP拟地面资料成像轨迹与观测井之间的区域)要大得多，且可对井中最浅接收点上方区域成像.
但是，VSP通常是三分量观测，需要对VSP资料作偏振分析，将P波数据从中分离出来.其次，分离后的P波数据，不仅包含直达P波和自由表面多次P波，还包含上行反射等其他波，如果直接利用方程(2)来重建拟地面地震反射波，将会产生很多假信号(Snieder et al., 2006)，导致重建的拟地面地震反射波资料信噪比降低.假设井中接收波场中包含了直达波D、上行反射波Ru和一阶下行鬼反射波Md，那么，两个地表震源xa和xb引起的在井中接收点zr的地震波场可表示为
|
(3a) |
|
(3b) |
根据方程(2)，对这两个地震波场进行互相关，则有
|
(4) |
根据稳相位分析，即假设井中接收点zr是地表震源xa和xb的稳相位点，那么，来自震源xa的下行直达波Da和震源xb的下行鬼反射波Mdb互相关的项Mdb⊗Da，就是希望获得的镜像反射波；而其他互相关的项不具有物理意义，将产生假信号，称为“串扰项或交叉项”.
为了能够最大限度避免产生假信号，解决串扰项的有效方法是进行波场分离(Mehta et al., 2007)，将下行直达波和一阶的下行鬼反射波从实测的VSP资料中分离出来.本文采用简单的中值滤波方法，来分离VSP中的上下行波场，然后，通过移动时窗，从分离后的下行波场中分别截取下行直达波和一阶的下行鬼反射波.
总之，VSP地震干涉成像方法是一种基于地震干涉的对多次波成像的方法，包含三个基本步骤：VSP上下行波场分离，拟地面地震资料重建及其偏移成像.其中，拟地面地震资料重建是VSP地震干涉成像的关键环节，而VSP上下行波场分离的好坏则直接决定拟地面地震资料重建质量，同时也决定了VSP地震干涉成像的质量.关于拟地面地震资料的偏移成像，则可采用现有的任何地面地震资料时间域或深度域方法来实现，本文不再赘述.
2 成像分辨率测试空间分辨率是衡量成像质量的重要指标之一，通常利用楔形体和水平排列不同尺度块体的成像来测试其垂向和水平极限分辨尺度.为了能比较不同深度异常体的VSP干涉成像结果，测试模型的大小都是2000 m×3000 m，采用同样的VSP观测参数和地震波场模拟参数.VSP观测参数：炮点位于地表，共201个，炮间距10 m；在模型左侧的井中接收，共101个接收点，道间距10 m，观测井段1000~2000 m.VSP声波有限差分地震波场模拟参数：网格间距5 m，Ricker子波中心频率40 Hz，记录长度3 s，采样率1 ms.由于在VSP地震波场模拟过程中已经结合Poyting矢量实现了上下行波分离，因此，无需再进行上下行波场分离.
2.1 垂向分辨率设计的楔形体模型如图 3所示，其中楔形体分别位于观测井段上方、中间和下方的模型区域，称为楔形体模型1、楔形体模型2和楔形体模型3.
|
图 3 楔形体模型示意，楔形体分别位于500 m、1500 m和2500 m深度，震源位于地表，观测井段位于模型左侧1000~2000 m. (a)楔形体模型1; (b)楔形体模型2; (c)楔形体模型3 Fig. 3 The wedge models: wedges are located at depths of 500 m, 1500 m, and 2500 m, the source is at the surface, and the observation well section is 1000~2000 m to the left of the model |
图 4a表示楔形体模型1模拟的VSP下行波地震波场:可清晰看到强的直达波(500 ms)和楔形体上下界面弱的一阶(1000 ms)和二阶(1500 ms)鬼反射波；由于VSP观测井段位于楔形体下方，故接收不到来自楔形体的反射波.
|
图 4 楔形体模型1：(a)模拟的零偏移距VSP共炮点道集；(b)重建的虚震源(位于第1道)道集；(c)与(b)相对应的模拟地面共炮点道集 Fig. 4 The wedge model 1: (a) Synthetic zero offset VSP common shot point gather; (b) Reconstructed virtual seismic source (located in trace 1) gather; (c) Synthetic surface common shot point gather with (b) |
将模拟的VSP共炮点资料分选为共接收点道集，通过移动时窗，分别截取下行直达波和下行多次波，使得方程(2)右边的u(zr, xa, t)只有下行直达波，而u(zr, xb, t)只包含下行鬼反射波，然后利用方程(2)重建拟地面地震反射波资料，共获得201个虚震源道集，每一个道集有201道，与VSP地表震源数目一致.图 4b表示虚震源位置位于模型左上角的重建拟地面地震道集，与其相应的模拟地面共炮点道集(图 4c，地震波场模拟参数与VSP的相同)比较：两者在反射波及多次波的到达时间(500 ms和1000 ms)一致，但子波不同，重建虚震源子波持续时间较长.这是因为：假设VSP观测地表震源子波为w(t)，那么，根据方程(2)，虚震源子波可表示为S(t)=|w(t)|2，是地表震源子波自相关的功率谱，它是零相位子波.其次，重建虚震源道集(图 5b)上没有直达波，是因为参与互相关的两个地震道，其中一个仅含有VSP下行直达波而另一个仅包含下行鬼反射波的缘故.此外，图 4b中还可看到由于互相关自身产生的线性干扰.
|
图 5 楔形体模型的VSP地震干涉成像结果. (a)楔形体模型1; (b)楔形体模型2; (c)楔形体模型3 Fig. 5 VSP seismic interferometric imaging results of wedge models. (a) Wedge model 1; (b) Wedge model 2; (c) Wedge model 3 |
为了观察拟地面地震资料中的噪声对偏移成像的影响，未作任何资料预处理，仅进行速度谱分析建立偏移速度模型，然后偏移成像，3个楔形体模型的VSP地震干涉成像结果如图 5所示.与相应的楔形体模型(图 3)比较：两者吻合良好，且垂向分辨率可达约10 m，不随楔形体深度增加而变化，只是成像范围稍窄.对于地面地震资料偏移成像，具有很宽范围的反射波入射角度是很重要的，较小角度的反射波有利于垂向分辨率，而较大角度的反射波有利于水平分辨率(Bleistein et al., 2001)，这一条件通常可通过震源或接收器的孔径宽度来满足.而VSP地震干涉成像则是利用多次波来满足该条件的(Jiang，2006)，对于同一个井中接收器，能接收到其近偏移距震源小角度入射的下行鬼反射波，且随着楔形体深度增加，下行鬼反射波的入射角变小，提高了其垂向分辨率.
此外，在楔形体模型1的VSP地震干涉成像结果(图 5a)中可以隐约看到：深度1000 m处有一个能量弱的成像楔形体，但比500 m的楔形体长度要小，这是模拟的VSP下行波场中含有高阶的鬼反射波所致.因为干涉积分中下行鬼反射波与下行直达波间的互相关运算，相当于对高阶的鬼反射波进行降阶处理，将一阶的鬼反射波降阶为拟地面地震资料中的反射波，将二阶的鬼反射波降阶为一阶的鬼反射波等等.而在楔形体模型2和3的VSP地震干涉成像结果(图 5b和c)，则未能观察到该现象，因为来自深度为1500 m或2500 m楔形体的二阶的鬼反射波到达时间，大于模拟记录的时间长度(3 s).
2.2 横向分辨率块体模型如图 6所示，其中块体分别位于观测井段上方、中间和下方的模型区域，称为块体模型1、块体模型2和块体模型3.
|
图 6 块体模型示意，块体分别位于500 m、1500 m和2500 m深度，震源位于地表，观测井段位于模型左侧1000~2000 m. (a)块体模型1; (b)块体模型2; (c)块体模型3 Fig. 6 Block models: blocks are located at depths of 500 m, 1500 m, and 2500 m, the source is at the surface, and the observation well section is 1000~2000 m to the left of the model. (a) Block model 1; (b) Block model 2; (c) Block model 3. |
图 7a表示块体模型1模拟的VSP地震波场，由于模型中有一个水平反射面以及10个大小不一的小块体，因此，与楔形体模型1的波场(图 4a)相比，块体模型1模拟的VSP下行波场要复杂得多，不仅有下行直达波，还有水平反射面和块体顶界面的下行多次波，以及来自块体上部两个角的绕射波等.图 7b表示虚震源位置位于模型左上角的重建拟地面地震道集，与其相应的模拟地面共炮点道集(图 7c)比较：两者在反射波的到达时间(500 ms)一致，但子波不同，重建的虚震源子波持续时间较长；其次，重建的远离井位置的块体上界面反射波(1000 ms以下)很弱，则是因其旅行路径较长，未进行几何扩散校正的缘故.
|
图 7 块体模型1：(a)模拟的零偏移距VSP共炮点道集；(b)重建的虚震源(位于第1道)道集；(c)与(b)相对应的模拟地面共炮点道集 Fig. 7 Block model 1: (a) Synthetic zero offset VSP common shot point gather; (b) Reconstructed virtual seismic source (located in trace 1) gather; (c) Synthetic surface common shot point gather with (b) |
图 8分别表示3个块体模型的VSP地震干涉成像结果，由于成像块体尺寸较小，将之垂向放大后置于成像目标的上方或下方.与相应的块体模型(图 6)比较：两者吻合良好，能分辨的块体宽度可达20 m；但随着块体深度增加，其分辨率变差，这是因为：对于同一个井中接收器，能接收到远偏移距大角度入射的多次波，随着块体深度增加，其入射的下行鬼反射波角度变小，使得水平分辨率变差.
|
图 8 块体模型的VSP地震干涉成像结果.图中标有“放大区域”是块体成像的垂向放大显示. (a)块体模型1; (b)块体模型2; (c)块体模型3 Fig. 8 VSP seismic interferometric imaging results of block models. The "enlarged area" marked in the figure is a vertically enlarged display of the block imaging. (a) Block model 1; (b) Block model 2; (c) Block model 3 |
综上所述，无论异常体位于井中观测段的上方、中间和下部区域，其VSP地震干涉成像的垂向分辨率可达10 m，横向分辨率可达20 m，这得益于近偏移距震源的小角度多次波入射和远偏移距震源的较大角度多次波入射.因此，在靠近观测井的区域，其垂向分辨率较高，而在较远区域则水平分辨率较高；随着异常体深度增加，入射多次波的角度变小，垂向分辨率提高但降低水平分辨率.
3 主要采集参数影响分析VSP观测中的接收点数量及其观测井段位置直接影响VSP观测效率和施工成本，也是VSP地震干涉成像的主要影响因素.在中国新疆地区，油气储层常常位于地下深部，且是缝洞型的，常规VSP反射波成像难以奏效.为此，根据新疆地区某工区的声波测井和零偏移距VSP资料，设计了一个含有5个不同深度溶洞(速度为2400 m·s-1)的水平层状模型(如图 9a所示)，目的是为了检验VSP地震干涉成像对地下深部溶洞的分辨能力，以及VSP观测参数对地震干涉成像质量的影响.
|
图 9 (a) 含有5个不同深度溶洞的水平层状模型; (b)模拟的地面地震资料成像结果 Fig. 9 (a) A horizontal layered model containing karst caves with five different depths; (b) Migration results of synthetic surface seismic data |
采用与该工区VSP实际观测相同的地表炮点参数，即地表炮点130个，炮间距50 m，最小偏移距50 m；井中观测井段位于模型左部，深度2000至5500 m，比实际观测井段3500~3730 m范围大得多，道间距10 m，故共有351个井中接收点.VSP声波地震波场模拟参数为：中心频率为35 Hz的Ricker子波，记录长度8.4 s，采样率1.2 ms.注意：由于VSP地震干涉成像是对自由表面鬼反射波成像，因此，VSP采集的记录时间长度通常是常规采集的1.5倍.
此外，为了检验VSP地震干涉成像的效果，也模拟了在地表激发接收的地面地震资料，激发点与接收点重合，都是130个，间距50 m，与利用VSP资料重建的拟地面地震资料中虚震源及其接收点的位置一致.对于模拟的地面地震资料，其深度域叠前偏移结果如图 9b所示.
3.1 接收器道数的影响首先对模拟的351道VSP资料分别进行VSP地震干涉成像和反射波成像，结果如图 10所示.与模拟的地面地震资料成像剖面(图 9b)相比：VSP地震干涉成像无论是水平反射面，还是低速溶洞，两者基本相同；仅成像范围有一定差异，VSP干涉成像范围呈倒梯形，比地面地震的稍小.但是，VSP干涉成像的垂向分辨率却比地面地震的高，尤其在靠近观测井区域，而水平分辨率比地面地震的稍低，这是因观测方式导致的，即VSP近偏移距震源的多次波入射角度，较地面地震反射波入射角度要小.与VSP反射波成像结果(图 10b)相比，VSP地震干涉成像的成像范围要大得多，能对观测井中最浅接收点上方区域成像；而且，水平反射面的成像同向轴连续性更好，分辨率更高；更重要的是，位于深度3200 m的两个大小分别为50 m×100 m和100 m×100 m的低速溶洞成像清晰，而深度6300 m处的3个大小分别为50 m×100 m、100 m×100 m和200 m×100 m的低速溶洞也能成像，仅比3200 m深度处的溶洞稍差.因为对于VSP观测方式，只能接收井中接收点下方地层的反射波，但却能接收到接收点上方和下方地层的鬼反射波；对于同一个反射点，具有更多角度入射的鬼反射波，且其入射角度小于反射波.
|
图 10 351道接收的VSP干涉成像与VSP反射波成像结果比较. (a) VSP地震干涉成像; (b) VSP反射波成像 Fig. 10 Comparison of 351 received VSP interferometric imaging and VSP reflection imaging results. (a) VSP seismic interferometric imaging; (b) VSP reflection imaging |
其次，从模拟的351道VSP资料中抽取观测井段中间的24道，分别进行VSP地震干涉成像和反射波成像，其结果如图 11所示.与351道的VSP地震干涉成像(图 10a)相比：两者无本质上的区别，只是信噪比稍低一些，成像范围稍窄一些，因为根据方程(2)和稳相位分析，观测井段中间24道接收到的VSP下行多次波足以重建拟地面地震的镜像反射波，只是重建的截断误差稍大一些.但是，24道的VSP反射波成像(图 11b)却比351道的差了很多，因为覆盖次数减少很多以及反射波入射角度变窄.
|
图 11 24道接收的VSP干涉成像与VSP反射波成像结果比较. (a) VSP地震干涉成像; (b) VSP反射波成像 Fig. 11 Comparison of 24 received VSP interferometric imaging and VSP reflection imaging results. (a) VSP seismic interferometric imaging; (b) VSP reflection imaging |
同样，从模拟的351道VSP资料中分别抽取观测井段顶部和底部12道资料进行VSP地震干涉成像，其结果如图 12所示.可以看到：底部12道接收的成像范围较顶部的要宽很多，但信噪比相对较低；两者对水平界面还有溶洞的成像效果差别不大，因为对于同一个地表震源和反射面，如图 12中的射线所示，较浅井中接收器接收到的是较大角度较短旅行路径的鬼反射波，而较深接收点接收到的则是较小角度和旅行路径较长的鬼反射波.
|
图 12 不同观测位置VSP干涉成像结果比较，射线表示最远偏移距震源的多次波路径. (a)观测井段顶部12道; (b)观测井段底部12道 Fig. 12 Comparison of VSP interferometric imaging results at different observation positions. The ray represents the multiple path of the farthest offset from the source. (a) 12 channels at the top of the observation well section; (b) 12 channels at the bottom of the observation well section |
另外与24道的VSP地震干涉成像(图 11a)相比，不管是顶部还是底部12道，其成像信噪比进一步降低，成像分辨率也随之有一定程度的降低，尤其是深度6300 m两个靠得很近的50 m×100 m和100 m×100 m溶洞，难以清晰地识别，这是因为随着接收点数的减少，方程(2)干涉积分的截断误差增大的缘故.
上述讨论表明，与VSP反射波成像相比，VSP地震干涉成像在如下三方面更具优势：①利用12至24级检波器串的VSP常规采集，VSP地震干涉成像就足以获得可接受的成像结果；②成像范围显著扩大，尤其是井中最浅接收点上方区域，与相应地面地震的相当或稍小；③成像同相轴连续性更好，分辨率更高，尤其是对深部较小溶洞较为精确的成像能力.这些优势是在不增加采集成本和不降低采集效率的条件下取得的，能够极大地拓展VSP的应用领域和范围，如可应用于难以进行大规模地面地震采集的山地，或应用于储层动态监测.
4 VSP实际资料试验2016年，在新疆地区进行高精度3D地面地震资料采集的同时，也进行了3D VSP联采.本文抽取其中一条过井的西向2D VSP变偏移距资料，用于VSP地震干涉成像试验，并与VSP反射波成像以及相同位置的高精度3D地面地震偏移成像结果进行了对比分析.
如图 13所示，五角星表示观测井位置，观测井段为3500~3730 m，24道3C检波器串接收，道间距10 m，记录长度6 s，采样率1 ms.红色点线表示地表炮线，共有激发点138个，炮点距为30 m，最小偏移距30 m，最大偏移距5250 m，激发点从左至右从小至大进行编号.注意：在第17道与第18道、第72道与第73道以及第91道与第92道之间，存在三处距离不等的炮点缺失，这将导致重建的拟地面地震资料在缺口下方接收点的覆盖次数为零，而其CDP覆盖次数呈非线性，影响VSP地震干涉成像的质量.
|
图 13 VSP观测地面炮点分布(红色点线)和接收点覆盖次数及相应的重建拟地面地震资料CDP覆盖次数.五角星表示观测井位置，虚线表示地表炮点位置 Fig. 13 The distribution of VSP ground shots (red dotted lines) and the number of coverages of receiving points and the corresponding number of CDPs of reconstructed pseudo-surface seismic data. The five-pointed star indicates the position of the observation well, and red dotted line indicates surface shot point |
图 14表示不同偏移距单炮(近、中、远偏移距)的原始垂直分量记录，近偏移距存在一定的干扰，中远偏移距波场上下行P波清晰，能量较强，容易识别地震波组类型，总体上而言，采集质量较高.
|
图 14 不同偏移距的VSP原始单炮垂直分量记录. (a)近偏移距; (b)中偏移距; (c)远偏移距 Fig. 14 VSP original single shot vertical component records with different offsets. (a) Near offset; (b) Middle offset; (c) Far offset |
根据本文方法原理，对于实测的VSP地震资料，其地震干涉成像处理流程如图 15所示.其中，VSP资料预处理及上下行波场分离可采用任何已有的VSP处理模块，本文采用ProMAX商业软件中的相关VSP处理模块；而重建拟地面地震资料的偏移成像，则是根据美国Colorado大学波现象研究组开发的开源平地表时间域叠前偏移成像算法改编的.此外，利用方程(2)进行拟地面地震资料重建，无需速度模型，也无需井中接收点的具体位置信息，因此，无需进行井中接收点静校正.这是VSP地震干涉成像方法的特点之一，避免了井中接收点静校正因速度和位置不准确而引入的误差.
|
图 15 VSP地震干涉成像处理流程示意 Fig. 15 Workflow of VSP seismic interferometric imaging |
根据VSP地震干涉成像处理流程，在利用方程(2)进行拟地面地震资料重建之前，需要对VSP地震记录进行上下行波场分离.经过激发点时差校正、初至时间拾取、球面扩散补偿、三分量记录合成、地震波组识别和分类等处理后，选择简单的中值滤波方法，进行上下行P波波场分离，结果如图 16所示，下行P波、上行反射P波都较为单一，特征明显，地震波组层次清晰，主频约为25 Hz.值得一提的是，上行P波可用于VSP反射波成像，而下行P波则可用于VSP地震干涉成像，其成像结果可相互验证，进一步确保偏移成像的可靠性.
|
图 16 VSP资料中一个共炮点道集的上下行波场分离结果. (a)下行P波; (b)上行P波 Fig. 16 Wavefield separation of a common shot point gather in VSP data. (a) Down-going P wave fields; (b) Up-going P wave fields |
对于获得的下行P波，通过移动时窗，分别截取下行直达波和下行多次波，使得方程(2)右边的u(zr, xa, t)只有下行直达波，而u(zr, xb, t)只包含下行多次波，然后利用方程(2)重建拟地面地震反射波，共获得138个虚震源道集，每个虚震源道集有138道.注意：重建的拟地面地震资料中虚震源及其接收点位置与VSP的地表炮点位置完全重合.
图 17a表示虚震源位于第138道的道集，隐约可见具有双曲线时差的重建反射波同相轴，但剖面上有很多的线性或非线性强干扰，是由互相关自身以及VSP的层间或短程下行多次波等与下行直达波相关所致，屏蔽了部分有效反射波信号，使得信噪比很低.对于此类线性干扰，可采用f-k滤波器来消除，结果如图 17b所示.由图 17b可知，基本消除了线性干扰，再现了重建的反射波，主频约为24 Hz，比上下行波分离后的下行波主频稍低.其次，可观察到第72道与第73道之间有较大的时移，在第17道与第18道之间和第91道与第92道之间有不同的时移，是因地表炮线上的3处激发点缺失引起的.
|
图 17 重建的虚震源(位于138道)道集(a)及其线性干扰消除后的相同道集(b) Fig. 17 The reconstructed virtual source (located at 138 channels) seismic gather (a) and the same gather (b) after the elimination of linear noise |
对重建的拟地面地震资料进行速度谱分析，以建立速度模型，平滑后的偏移成像速度模型如图 18所示；然后偏移成像，其时间域的叠前偏移结果如图 19中的深灰色梯形区域所示，图中，底图表示相同测线的3D地面地震时间域叠前偏移结果，浅灰色锥形区域表示VSP反射波时间域叠前偏移结果.与VSP反射波成像相比，VSP地震干涉成像的成像范围显著扩大.与相同位置的地面地震成像相比，VSP地震干涉成像的分辨率较低，这主要是因未对重建的拟地面地震资料作子波反褶积处理造成的.在成像剖面浅中部(0~2000 ms)，信噪比较高，但由于受地表 3处激发点缺失影响，VSP地震干涉成像同相轴畸变严重，呈现陡倾角形态或空缺；而在成像剖面深部(大于2000 ms)，两者的主要同相轴基本一致.
|
图 18 根据速度谱分析建立的平滑速度模型 Fig. 18 Smooth velocity modeled by velocity spectrum analysis |
|
图 19 拟地面地震资料时间域叠前偏移成像结果，即VSP地震干涉成像结果(深灰色梯形区域).底图表示相同位置的3D地面地震资料时间域叠前偏移结果；浅灰色锥形区域表示VSP反射波时间域叠前偏移结果 Fig. 19 Time-domain prestack migration imaging results of pseudo-surface seismic data, that is, VSP seismic interferometric imaging results (dark gray trapezoidal area). The base map shows the time-domain prestack migration results of 3D surface seismic data at the same location; the light gray cone area shows the VSP time-domain prestack reflection migration results |
其次，将VSP地震干涉成像剖面与相同位置的地面地震成像结果进行比较，如图 20所示，分别表示成像剖面上部(0~1500 ms)和下部(1800~4500 ms)的对比.为了能更好地进行比较，将VSP地震干涉成像井旁道叠加，作为井旁道标定(重复显示的6个蓝色地震道)置于图中.可以看到：在成像剖面上部(图 20a)0~400 ms之间，由于井中检波器串接收不到来自远偏移距地表震源一阶下行鬼反射波，VSP地震干涉成像效果差；随深度增大，其成像效果也随之改善，甚至比地面地震的成像同相轴更加连续，信噪比也较高，如400~900 ms之间的同相轴；但在1100 ms处有一定的差异，是VSP地表激发点缺失引起的.在成像剖面下部(图 20b)，VSP地震干涉成像受地表激发点缺失的影响相对较小，且其成像同相轴大体上与地面地震的相吻合，但也有差异，主要表现在局部同相轴增多，可能是因为局部异常体作为次级震源，增加了入射波角度的多样性，从而改善了成像效果，也有可能是由分离后的下行P波中还含有层间下行多次波导致的，具体原因尚有待进一步试验和分析.
|
图 20 VSP地震干涉成像剖面与相同位置地面地震成像剖面的比较. (a)成像剖面上部(0~1500 ms)的对比; (b)成像剖面下部(1800~4500 ms)的对比. (a)和(b)图中左边的表示地面地震成像结果，中间重复显示的6个地震道表示VSP地震干涉成像井旁道的叠加，右边表示VSP地震干涉成像结果 Fig. 20 Comparison of VSP seismic interferometric imaging profile with ground seismic imaging profile at the same location. The left side of (a) and (b) shows the results of surface seismic imaging. The six seismic traces repeatedly shown in the middle represent the superposition of the well sideways of the VSP seismic interferometric imaging, and the right side shows the results of VSP seismic interferometric imaging. (a) Comparison of the upper part of the imaging section (0~1500 ms); (b) Comparison of the lower part of the imaging section (1800~4500 ms) |
最后，将VSP地震干涉成像剖面与相同位置VSP反射波成像结果进行比较，如图 21所示.图中间重复显示的6个VSP井旁道标定，是根据VSP反射波成像结果的井旁地震道叠加而成.可以看到：VSP地震干涉成像范围更大，与VSP反射波的成像同相轴能较好对应，且同相轴连续性更好，但主频略低.因为VSP反射波只经过近地表一次，而下行鬼反射波则需经过近地表三次，近地表的衰减吸收使得其频率降低.
|
图 21 VSP地震干涉成像剖面与相同位置VSP反射波成像剖面的比较.图中左边为VSP地震干涉成像结果，中间重复显示的6个地震道表示VSP反射波成像井旁道的叠加，右边是VSP反射波成像结果 Fig. 21 Comparison of VSP seismic interferometric imaging profile with VSP reflection imaging profile at the same location. The left side of the figure shows the results of VSP seismic interferometric imaging. The six seismic traces repeatedly shown in the middle represent the superposition of the well sideways of the VSP reflection imaging, and the right side shows the results of VSP reflection imaging |
本文研发了VSP地震干涉成像方法技术，通过数值模型测试了它的纵横向极限分辨率，并讨论了主要采集参数的影响.通过新疆地区某个VSP实际资料的地震干涉成像，建立了其成像处理流程，并与相应的VSP反射波成像和地面地震成像结果进行了对比分析，有如下结论：
① 地震干涉成像的垂向和水平极限分辨率分别可高达约10 m和20 m，且能分辨6500 m深度处50 m×100 m的小溶洞.这对于中国新疆地区深部缝洞型储层的成像具有重要意义.
② 利用12或24道井中检波器串的VSP观测，其地震干涉成像的范围比VSP反射波成像范围显著扩大，尤其是在井中最浅接收位置上方区域，成像同相轴在浅部更加连续，与地面地震资料的成像范围相当.这有利于成像同相轴的追踪对比、地震属性提取和地质解释，使得VSP技术更具适用性和实用性，可进一步应用于油气藏的动态监测或EOR过程中.
③ VSP地震干涉成像无需对井中接收点做静校正，也不需要知道震源位置和激发时刻信息，因为拟地面地震资料重建(互相关)过程本身能够自动消除接收点静校正误差，这对于3D DAS-VSP数据的多次波成像尤为重要，因为DAS技术的主要缺点之一是难以准确定位接收点位置.同时，这也意味着VSP干涉成像可应用于如地面微地震监测或环境噪声的被动源成像，因为如地下岩石破裂或断层滑动等的震源与地表接收器之间能够组成一个IVSP观测方式.
但是新方法在实际应用中仍面临着一些技术挑战，如受VSP观测系统的局限，井中接收器分布是1D的，只能重建2D的拟地面地震资料.因此，采用以井位置为中心的柱坐标方位角，抽取共方位角度道集，分别进行2D拟地面地震资料重建及其成像，然后通过插值形成3D偏移成像剖面，而且在重建拟地面地震资料之前，还需要校正3D的几何扩散影响；再如VSP地震干涉成像是基于无损耗介质的声波相关型互换方程的，未考虑介质衰减效应，会导致重建的振幅不精确；同时还可以考虑二阶的鬼反射波压制或利用，以及克服复杂介质的散焦问题等.随着这些问题的解决，势必会大大提升新方法的成像效果.
Bakulin A, Calvert R. 2006. The virtual source method:Theory and case study. Geophysics, 71(4): SI139-SI150. DOI:10.1190/1.2216190 |
Bakulin A, Calvert R. 2008. Virtual Source Method: overview of history and development.//78th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2726-2730.
|
Blestein N, Cohen J, et al. 2001. Mathematics of Multidimensional Seismic. Imaging, Migration and Inversion, Berlin, Springer Verlag. |
Bojarski N. 1983. Generalized reaction principles and reciprocity theorems for the wave equations, and the relationship between the time-advanced and time-retarded fields. Journal of the Acoustical Society of America, 74: 281-285. DOI:10.1121/1.389721 |
Claerbout J. 1968. Synthesis of a layered medium from its acoustic transmission response. Geophysics, 33(2): 264-269. DOI:10.1190/1.1439927 |
Campillo M, Paul A. 2003. Long-range correlations in the diffuse seismic coda. Science, 299: 547-549. DOI:10.1126/science.1078551 |
Chen G J, Cao H, Wu Y S, et al. 2013. The seismic virtual source method and numerical experiments. Progress in Geophysics (in Chinese), 28(5): 2725-2733. |
Chen G J, Guo J, Zhang Y H, et al. 2017. Technique for seismic response retrieval from ambient noise and its application. Geophysical Prospecting for Petroleum (in Chinese), 56(6): 798-803. |
Chiffot C, Prescott A, et al. 2017. Data-driven interferometry method to remove spatially aliased and non-linear surface waves.//87th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 4980-4984.
|
De Hoop A T. 1988. Time-domain reciprocity theorems for acoustic wave fields in fluids with relaxation. Journal of the Acoustical Society of America, 84: 1877-1882. DOI:10.1121/1.397152 |
Draganov D, et al. 2009. Reflection images from ambient seismic noise. Geophysics, 74(5): A63-A67. DOI:10.1190/1.3193529 |
Draganov D., Wapenaar K, et al. 2007. Retrieval of reflections from seismic background-noise measurements. Geophysical Research Letters, 34: L04305-4. |
Fink M. 1997. Time reversal acoustics. Phys. Today, 50: 30-40. DOI:10.1063/1.881845 |
Fokkema J T, et al. 1993. Seismic applications of acoustic reciprocity. Elsevier Science Publishing Company, Inc. |
Forghani F, Snieder R. 2010. Underestimation of body waves and feasibility of surface-wave reconstruction by seismic interferometry. The Leading Edge, 29: 790-794. DOI:10.1190/1.3462779 |
Fang L H, Wu J P, Lü Z Y. 2009. Rayleigh wave group velocity tomography from ambient seismic noise in North China. Chinese Journal of Geophysics (in Chinese), 52(3): 663-671. DOI:10.1002/cjg2.1388 |
He R Q, Hornby B, Schuster G T. 2007. 3D wave-equation interferometric migration of VSP free-surface multiples. Geophysics, 72(5): S195-S203. DOI:10.1190/1.2743375 |
Halliday D, Curtis A, et al. 2010. Interferometric ground-roll removal:Attenuation of scattered surface waves in single-sensor data. Geophysics, 75(2): SA15-SA25. |
Jiang Z Y, Yu J H, Schuster G T, et al. 2005. Migration of multiples. The Leading Edge, 24(3): 315-318. DOI:10.1190/1.1895318 |
Jiang Z.2006. Migration and attenuation of surface related and interbed multiples[Ph.D. dissertation], University of Utah, Salt Lake City, Utah.
|
Liu B H. 2008. Encyclopedia of Petroleum Exploration and Development in China (in Chinese). Beijing: Petroleum Industry Press.
|
Mateeva A, Mestayer J, Cox B, et al. 2012. Advances in distributed acoustic sensing (DAS) for VSP.//82nd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1903-1907.
|
Mateeva A, Mestayer J, Yang Z, et al. 2013. Dual-well 3D VSP in deepwater made possible by DAS.//83rd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 5062-5066
|
Mateeva A, Lopez J, Mestayer J, et al. 2013a. Distributed acoustic sensing for reservoir monitoring with VSP. The Leading Edge, 32(10): 1278-1283. DOI:10.1190/tle32101278.1 |
Mateeva A, Lopez J, Potters H, et al. 2014. Distributed acoustic sensing for reservoir monitoring with vertical seismic profiling. Geophysical Prospecting, 62(4): 679-692. DOI:10.1111/1365-2478.12116 |
Mateeva A, Hornman J C, et al. 2015. Frequent seismic monitoring for pro-active reservisor management.//85th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 4817-4821.
|
Mehta K A, Bakulin J, Sheiman R, et al. 2007. Improving the virtual source method by wavefield separation. Geophysics, 72(4): V7-V86. |
Qi C, Chen Q F, Chen Y. 2007. A new method for seismic imaging from ambient seismic noise. Progress in Geophysics (in Chinese), 22(3): 771-777. DOI:10.1016/S1872-5791(07)60044-X |
Rickett J, Claerbout J. 1999. Acoustic daylight imaging via spectral factorization:Helioseismology and reservoir monitoring. The Leading Edge, 18: 957-960. DOI:10.1190/1.1438420 |
Shapiro N, Campillo M, Stehly L. 2005. High-resolution surface-wave tomography from ambient seismic noise. Science, 307: 1615-1618. DOI:10.1126/science.1108339 |
Schuster G T, Zhou M. 2006. A theoretical overview of model-based and correlation-based redatuming methods. Geophysics, 71(4): SI103-SI110. DOI:10.1190/1.2208967 |
Schuster G T. 2009. Seismic Interferometry. Cambridge: Cambridge University Press.
|
Snieder R, et al. 2002. Coda wave interferometry for estimating nonlinear behavior in seismic velocity. Science, 295: 2253-2255. DOI:10.1126/science.1070015 |
Snieder R. 2004. Extracting the Green's function from the correlation of coda waves:A derivation based on stationary phase. Physical Review E, 69: 046610. DOI:10.1103/PhysRevE.69.046610 |
Snieder R, Wapenaar K, Larner K. 2006. Spurious multiples in seismic interferometry of primaries. Geophysics, 71(4): SI111-SI124. DOI:10.1190/1.2211507 |
Vasconcelos I, et al. 2008. Imaging internal multiples from subsalt VSP data-examples of target-oriented interferometry. Geophysics, 73(4): S157-S168. DOI:10.1190/1.2944168 |
Wapenaar K, Fokkema J. 2006. Green's function representations for seismic interferometry. Geophysics, 71(4): SI33-SI46. DOI:10.1190/1.2213955 |
Wapenaar K, Draganov D, Snieder R, et al. 2010. A representation for Green's function retrieval by multidimensional deconvolution. J. Acoust. Soc. Am, 128(6): EL366-EL371. DOI:10.1121/1.3509797 |
Wapenaar K, et al. 2010a. On seismic interferometry:the generalized optical theorem and the scattering matrix of a point scatterer. Geophysics, 75(3): SA27-SA35. DOI:10.1190/1.3374359 |
Wapenaar K, Draganov D, Snieder R, et al. 2010b. Tutorial on seismic interferometry:Part 1-Basic principles and applications. Geophysics, 75(5): 75A195-75A209. DOI:10.1190/1.3457445 |
Wapenaar K, Draganov D, Snieder R, et al. 2010c. Tutorial on seismic interferometry:Part 2-Underlying theory and new advances. Geophysics, 75(5): 75A211-75A227. DOI:10.1190/1.3463440 |
Wapenaar K, et al. 2011. seismic interferometry by crosscorrelation and by multidimensional deconvolution:a systematic comparison. Geophysical Journal International, 185: 1335-1364. DOI:10.1111/j.1365-246X.2011.05007.x |
Wang J N, Wu G X, Chen X F. 2019. Frequency-Bessel transform method for effective imaging of higher-mode Rayleigh dispersion curves from ambient seismic noise data. Journal of Geophysical Research:Solid Earth, 124(4): 3708-3723. DOI:10.1029/2018JB016595 |
Wu H, Wong W F, Yang Z H, et al. 2015. Dual-well 3D vertical seismic profile enabled by distributed acoustic sensing in deepwater Gulf of Mexico. Interpretation, 3(3): SW11-SW25. DOI:10.1190/INT-2014-0248.1 |
Wu X, Willis M E, Palacios W, et al. 2017. Compressional-and shear-wave studies of distributed acoustic sensing acquired vertical seismic profile data. The Leading Edge, 36(12): 987-993. DOI:10.1190/tle36120987.1 |
Xu P F, Li C J, Lin S Q, et al. 2009. Mapping collapsed columns in coal mines utilizing Microtremor Survey Methods. Chinese Journal of Geophysics (in Chinese), 52(7): 1923-1930. |
Yu J H, Schuster G T. 2006. Crosscorrelogram migration of inverse vertical seismic profile data. Geophysics, 71(1): S1-S11. |
Zhan G, Kommedal J, Nahm J. 2015. VSP field trials of distributed acoustic sensing in Trinidad and Gulf of Mexico.//85th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 5539-5543.
|
Zhu G M. 1988. Vertical Seismic Profiling (in Chinese). Beijing: Petroleum Industry Press.
|
Zwartjes P, Mateeva A, Chalenski D, et al. 2018. Frequent, multi-well, stand-alone 3D DAS VSP for low-cost reservoir monitoring indeepwater.//88th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 4948-4952.
|
陈国金, 曹辉, 吴永栓, 等. 2013. 虚源法地震技术及数值模型试验. 地球物理学进展, 28(5): 2725-2733. |
陈国金, 郭建, 张亚红, 等. 2017. 基于环境噪声的地震响应重建方法及应用. 石油物探, 56(6): 798-803. DOI:10.3969/j.issn.1000-1441.2017.06.004 |
房立华, 吴建平, 吕作勇. 2009. 华北地区基于噪声的瑞利面波群速度层析成像. 地球物理学报, 52(3): 663-671. |
刘宝和. 2008. 中国石油勘探开发百科全书-勘探卷. 北京: 石油工业出版社.
|
齐诚, 陈祺福, 陈顒. 2007. 利用背景噪声进行地震成像的新方法. 地球物理学进展, 22(3): 771-777. DOI:10.3969/j.issn.1004-2903.2007.03.017 |
徐佩芬, 李传金, 凌胜群, 等. 2009. 利用微动勘察方法探测煤矿陷落柱. 地球物理学报, 52(7): 1923-1930. DOI:10.3969/j.issn.0001-5733.2009.07.028 |
朱光明. 1988. 垂直地震剖面法. 北京: 石油工业出版社.
|