2. 中国石油化工股份有限公司油田勘探开发事业部, 北京 100728
2. Oil field Exploration & Development Department Sinopec, Beijing 100728, China
如何从实测地震资料或环境噪声中, 提取地震响应(或格林函数), 一直是地球物理学家研究的焦点。近年来快速发展的地震干涉法[1-3]旨在解决此类问题, 其思想是, 通过两个地震记录的互相关运算, 来恢复该对记录所在接收点位置之间的地震响应, 即在其中一个接收点处激发(虚源), 而在另一个接收点处接收的地震响应。这一互相关过程被称为虚源构建, 是地震干涉法的本质, 该过程无需介质的速度模型, 也无需震源位置和激发时刻等信息[3]。
地震干涉法的起源由来已久, 它始于CLAERBOUT[4]提出的反射波与透射波响应之间的关系式; 他猜想该关系式能被推广应用于3D介质, 也讨论了其在日震学和油藏监测等领域的应用, 相关方法被命名为声波日光成像[5]方法, 即利用环境噪声对地下介质成像的方法。SCHUSTER等[6]利用稳相位方法证实了CLAERBOUT的猜想, 提出了“干涉成像方法”, 将该方法应用于VSP资料[7-9], 开创性地利用地震干涉进行地面地震、垂直地震、单井地震和井间地震采集基准面之间的转换[3]。WAPENAAR等[10-12]基于相关型和褶积型互换定理, 在数理上严格证明了地震干涉法相关型和褶积型控制方程, 奠定了其理论框架。SNIEDER等[13]利用天然地震尾波提取地下信息, 并根据稳相位方法, 详细解释了单散射介质中炮检基准面的重建过程[14]和虚源构建过程中产生非物理信号的原理[15], 极大地推动了地震干涉法的发展。CALVERT等[16]及BAKULIN等[17]受逆时声学的启发, 提出了虚源方法, 并获得美国专利授权, 该方法采用类似VSP的采集方式(如水平井或斜井), 通过构建虚源, 将地表炮点重建至井中接收点处, 避免了近地表上覆复杂介质的透射响应影响, 该方法主要应用于高陡构造成像和4D油藏监测[18-21]。自此, 地震干涉法进入蓬勃发展的阶段, 它既可应用于主动源地震资料, 也可应用于被动源地震资料, 期间出现了很多富有想象力的应用。
目前, 被动源地震干涉法主要应用于如海洋潮汐等50~100 s长周期环境噪声[22-24], 或如火山活动等0.2~2.0 Hz噪声[25-27], 或3~15 Hz的交通噪声[22, 28-30]条件下的面波恢复及其横波速度反演; 也可应用于反射体波恢复。DRAGANOV等[31-34]通过数值模型和实际资料展示了反射波恢复条件及注意事项, 王德利等[35]及张盼等[36]根据数值模型模拟的环境噪声资料, 采用互相关或反褶积算法恢复了反射波, 证明了利用背景噪声进行反射波恢复的可行性。
我们利用地震干涉法, 探索微地震地面监测资料的反射波恢复问题。针对水力压裂引起的地下岩石破裂问题, 如果仅取互相关的正延迟函数来生成虚源地震道, 将会丢失其负延迟中的反射波信息, 因此提出求和方法和相对位置方法来生成虚源地震道, 并利用数值模型和四川某工区采集的微地震地面监测资料对上述方法进行验证。
1 方法原理在任意的各向异性弹性无损介质中, 基于“归一化通量”假设条件恢复噪声地震干涉的弹性动力学格林函数可表示为[11]:
| $ \begin{array}{l} {\left\{G_{p, q}\left(x_{\mathrm{B}}, x_{\mathrm{A}}, t\right)+G_{p, q}\left(x_{\mathrm{B}}, x_{\mathrm{A}},-t\right)\right\} *} \\ {S_{\mathrm{N}}(t) \approx \frac{2}{\rho c_{\mathrm{P}}}\left\langle u_{p}\left(x_{\mathrm{A}},-t\right) * u_{q}\left(x_{\mathrm{B}}, t\right)\right\rangle} \end{array} $ | (1) |
式中:xA代表位置A; xB代表位置B; Gp, q(xA, xB, t)和Gp, q(xA, xB, －t)(p, q=1, 2, 3)分别表示在位置A处p方向激发(虚源)、位置B处接收的格林函数及其相应逆时格林函数的q分量; t表示时间; ρ, cP表示封闭面外任意选取的各向同性均匀声波介质的密度和P波速度, 可视作常数; SN(t)表示虚源子波, 它是噪声源的自相关函数; 〈·〉表示总体平均; *表示褶积。假设噪声震源不相关, 如果Nl(x, t)为作用方向l(l=1, 2, 3)的噪声震源, 那么在地表xA和xB位置处接收的p, q分量up(xA, t)和uq(xB, t)可表示为:
| $ \left\{ {\begin{array}{*{20}{l}} {\left\langle {{N_l}(x, - t)*{N_k}\left( {{x^\prime },t} \right)} \right\rangle = {\delta _{kl}}\delta \left( {x - {x^\prime }} \right){S_{\rm{N}}}(t)}\\ {{u_q}\left( {{x_{\rm{B}}},t} \right) = \oint\limits_{\partial D} {} {G_{q,k}}\left( {{x_{\rm{B}}},{x^\prime },t} \right)*{N_k}\left( {{x^\prime },t} \right){d^2}{x^\prime }}\\ {{u_p}\left( {{x_{\rm{A}}},t} \right) = \oint\limits_{\partial D} {} {G_{p,l}}\left( {{x_{\rm{A}}},x,t} \right)*{N_l}(x,t){d^2}x} \end{array}} \right. $ | (2) |
式中:δkl表示Kroneckerδ函数; δ(x－x′)表示狄拉克函数。
被动源地震干涉法虽然在理论上能够同时恢复包括面波和反射波在内的全格林函数, 但受震源数量及其分布等因素的制约, 实际是对面波和反射波分别进行恢复[11-12]。
本文将(1)式应用于微地震地面监测资料, 即利用在地表接收的包含自由表面多次波在内的透射波响应, 恢复地面地震镜像反射波。根据稳相位分析[1, 3, 14]结果可知, 当地下噪声震源是地表两个接收点A和B的稳相位震源时, 如果对这两个接收点的地震记录进行互相关, 即将在位置B记录的一阶自由表面多次波波至时间, 减去在接收点A的直达波波至时间, 那么可得到在位置A激发(虚源)、位置B接收的镜像反射波(图 1);以此类推, 对于高阶自由表面多次波, 被直达波所减, 可以得到次一阶多次波; 被次一阶多次波所减, 可以得到镜像反射波。这也解释了要求(2)式的“环境噪声记录的长度需足够长”的原因, 即能够记录到足够的稳相位噪声源地震波。环境噪声记录的具体长度[34]则要视该地区的微地震活动是否活跃而定, 可能为几个小时, 也可能为几周或几个月。
|
图 1 拟地面地震镜像反射波恢复的稳相位分析示意 |
在无耗损介质中, 格林函数对时间的偏导数是二阶的, 具有逆时不变性[2], 也就是说, 地震响应和它的逆时地震响应都是波动方程的解, 分别对应于(1)式右边波场互相关的正延迟和负延迟函数, 在满足“每个优势波长范围内至少有2个不相关震源, 规则分布于环绕接收点的一个封闭面上”[31]的条件下, 恢复的地震响应及其逆时地震响应是相同的。因此, 在生成虚源地震道时, 为节省计算时间, 通常只计算互相关的正延迟函数, 并将其作为恢复的地震响应, 上述方法称为常规方法。
上述方法的前提条件过于严格, 实际上往往不能得到满足。如在微地震地面监测中, 水力压裂引起的岩石破裂, 往往集中于压裂井段附近约200 m的水平范围内。如图 2所示, 震源位于地面测线左下角位置, 根据稳相位分析结果可知, 将接收点A作为虚源在另一个接收点B接收时, 根据(1)式, 需要将A处记录的波场在时间上反转(逆时), 然后与B处记录的波场褶积。互相关在数学上表示相位相减, 即B处地震记录的多次波波至时间减去A处地震记录的直达波波至时间, 由于从地下震源传播至地表B处的一阶鬼波旅行时间大于传播至地表A处的直达波时间, 因此重建的反射波(如图 2中的蓝色射线路径所示)位于互相关的正延迟函数中; 但当虚源位置不变而在位置C接收时, 因为C处记录的直达波波至时间小于A处记录的多次波波至时间, 故重建的反射波(如图 2中的红色射线路径所示)位于互相关的负延迟中。
|
图 2 稳相位分析示意 |
上述分析说明, 当震源分布位于测线左下角位置时, 如果利用常规方法来生成虚源地震道, 将造成互相关负延迟中恢复的镜像反射波丢失, 意味着虚源左侧地震道是噪声道, 使得覆盖次数减少一半。因此, 应同时计算互相关的正延迟和负延迟函数, 并将二者相加, 即将格林函数及其相应的逆时格林函数相加, 来获得虚源地震道U(xB, xA, t), 即:
| $ \begin{aligned} U\left(x_{\mathrm{B}}, x_{\mathrm{A}}, t\right)=&\left\{G\left(x_{\mathrm{B}}, x_{\mathrm{A}}, t\right)+G\left(x_{\mathrm{B}}, x_{\mathrm{A}},\right.\right.\\ &-t)\} * S_{\mathrm{N}}(t) \end{aligned} $ | (3) |
我们将上述方法称之为求和方法, 其优点是将丢失的恢复反射波重新利用, 覆盖次数增加了1倍; 缺点是计算量增加了1倍, 也增加了虚源道集的虚假同相轴, 降低了信噪比。
为了不降低信噪比, 也可根据虚源及其接收点位置与压裂井位置的相对关系(图 2), 来获得虚源地震道, 当接收点B位于地下震源和虚源位置A的右侧时, 取互相关的正延迟函数作为虚源地震道U(xB, xA, t), 即:
| $ U\left(x_{\mathrm{B}}, x_{\mathrm{A}}, t\right)=G\left(x_{\mathrm{B}}, x_{\mathrm{A}}, t\right) * S_{\mathrm{N}}(t) $ | (4a) |
当接收点C位于地下震源和虚源位置A之间时, 取互相关的负延迟函数作为虚源地震道U(xC, xA, t), 即:
| $ U\left( {{x_C},{x_{\rm{A}}},t} \right) = G\left( {{x_C},{x_{\rm{A}}}, - t} \right)*{S_{\rm{N}}}(t) $ | (4b) |
这种虚源地震道生成方法, 称为相对位置方法, 该方法既将丢掉的恢复反射波重新利用, 覆盖次数增加了1倍, 又不降低虚源地震道的信噪比。
2 数值模型验证采用THORBECKE等[34]设计的凹陷数值模型及5个噪声源子波(图 3), 比较常规方法、求和方法和相对位置方法的优劣, 所用的模型大小为10 000 m×4 100 m, 模型中黑点表示噪声源位置, 包含4个反射界面, 起伏界面为强反射面, 速度从上至下分别为:1 500, 2 000, 4 000, 3 000, 5 500 m/s。在地表等间距(20 m)分布着501个检波点, 接收来自地下被动源的透射波响应。噪声源子波(图 3b)为随机时间序列, 中心频率为35 Hz, 时间长度为5~90 s不等。利用声波有限差分模拟波场, 记录时间长度为120 s, 采样率为2 ms。
|
图 3 凹陷数值模型(a)及5个噪声源子波(b) |
首先, 测试理想情况下的反射波恢复能力。假设地下被动源共有1 000个, 随机分布在深度为500~4 000 m的模型空间内, 保证水平方向上每个优势波长范围内有2个被动源。采用常规方法生成的虚源(位于地表 5 000 m处)道集如图 4a所示, 利用主频为30 Hz的Ricker子波模拟的合成地震记录如图 4b所示, 二者相比可知, 400, 1 100, 2 100, 2 600 ms的4个界面的反射波均被恢复, 1 050 ms处的一阶鬼波和2 260m s的微曲多次波也被恢复, 此外, 恢复的直达波到时误差较大, 根据稳相位分析结果可知, 恢复直达波的稳相位震源区域位于地表[12], 而非地下。道集包含的噪声为线性或非线性噪声, 且信噪比低, 原因在于:在地表任意接收点接收的透射波场中, 除了包含直达波、一阶和二阶鬼波, 还包含面波、噪声震源下部的反射波、层间多次波、转换波以及随机噪声等, 根据(1)式, 对上述波场进行互相关, 不仅会产生反射波(直达波与一阶鬼波相关, 一阶鬼波与二阶鬼波相关)、鬼波(直达波与二阶鬼波相关)和层间多次波(直达波与层间多次波相关)等具有物理意义的同相轴, 而且还会产生很多非物理同相轴[15](直达波与反射波相关、反射波与鬼波相关等)。理论上, 如果噪声源数量及其分布满足被动源地震干涉的假设条件, 上述非物理同相轴在相关过程中将会相互抵消, 但实际上是不可能完全相互抵消的, 因此, 需要对实际资料做一些针对性的预处理, 如带通滤波和面波消除[32-33]等。由(2)式可知, 恢复反射波后的地震子波是噪声震源的自相关函数, 也是零相位子波, 其旁瓣可能有多个。
|
图 4 基于凹陷模型采用常规方法生成的虚源道集(a)与Ricker子波模拟的合成地震记录(b) |
其次, 将噪声震源(500个)置于凹陷模型的左下角位置(图 5a), 分别采用常规方法、求和方法和相对位置方法来生成虚源地震道, 结果如图 5所示。
|
图 5 左下部为500个噪声源的凹陷模型(a)和采用常规方法(b)、求和方法(c)和相对位置方法(d)生成的虚源道集(虚源位于地表第251道处) |
图 5b为采用常规方法生成的虚源道集, 正如上一节中分析的那样, 因为恢复的反射波位于互相关的负延迟中, 且被丢失了, 所以与噪声震源分布较为理想的虚源道集(图 4a)相比, 在虚源(位于地表第251道处)左侧的地震道上观察不到恢复的反射波, 仅是一些线性或非线性噪声; 而在虚源右侧的地震道上, 仅能够观察到恢复的部分反射波, 但恢复效果较差。浅部的假相干同相轴增多, 是因噪声震源分布于地表测线的左下角, 相对于地表接收点不对称所致[1]。
图 5c为采用求和方法生成的虚源道集, 可以看出在虚源两侧地震道上均恢复出了反射波, 将之前丢失的反射波重新加以利用, 可增加1倍的覆盖次数; 虚源左侧地震道上恢复的反射波较差, 与虚源右侧反射波不对称, 这是因为互相关函数不是关于时间的偶函数, 互相关过程不具备可交换性[37]; 互相关过程中引入了负延迟函数中的噪声(如右上倾的线性噪声), 降低了信噪比。
图 5d为采用相对位置方法生成的虚源道集, 虚源道集两侧恢复的反射波与采用求和方法恢复的反射波相同, 但却避免了噪声引入。采用相对位置方法未降低信噪比, 且显著改善偏移成像结果质量, 对微地震地面监测实际资料的反射波恢复效果良好。
3 实际资料试验利用四川某工区水力压裂采集的微地震地面监测资料, 根据(1)式构建反射波响应, 分别采用求和方法和相对位置方法生成虚源道集, 并对获得的虚源道集进行动校叠加, 测试两种方法的效果。
该工区共包括两段压裂, 深度为2 600~2 900 m。在地表布置了10条测线, 以压裂井为中心呈放射状等间隔分布; 每条测线123道, 每道放置12只垂向分量检波器(10 Hz低切滤波器)组合压制面波, 道间距为25 m; 共有1 230道。每段压裂开始前30分钟开始记录, 连续记录, 每30 s记录一个数据文件, 直至压裂停泵4小时后停止, 共获得2 018个有效数据文件。抽取第3条测线的监测资料, 并选取含有较强微地震事件的数据文件379个(共计5 910 s)。图 6为选取的第442个文件第25~30 s的原始地震记录及其频谱分析结果, 频谱成分复杂, 频带宽; 微地震事件的初至波清晰可辨, 能量强, 频率约为37 Hz, 有一定的延续性; 地震记录残留约16 Hz的高频面波, 位于靠近压裂井的第1~5、第20~30和第80~90地震道; 第115道附近存在着约50 Hz的工业电干扰以及其它持续时间短的随机噪声, 因而各地震道间振幅差异大, 且淹没了自由表面多次波。
|
图 6 四川某工区微地震地面监测原始地震记录(a)及其频谱分析结果(b) |
恢复反射波利用的主要是一阶及以上全程多次波, 但受传播距离长、介质的吸收衰减等因素影响, 反射波能量弱, 频率也低于微地震初至波, 且可能与噪声的频率部分折叠, 过度预处理将伤及有效波。为了不损害自由表面多次波, 我们选择保留大部分噪声, 虽然会一定程度上影响反射波恢复质量, 但可以在恢复反射波后再进行噪声消除。
由图 6b的频谱分析结果可知, 选取合适参数的带通滤波可以消除低频和甚高频随机干扰, 结果如图 7a所示。与图 6a的原始地震记录相比, 带通滤波后的地震记录基本保持原貌, 微地震初至的后续波稍有增强。为了能更好地进行地震道之间的互相关, 采用逐道均方根振幅归一化处理, 消除地震道间的振幅差异, 均衡各地震道间的振幅, 且保留每道的振幅相对关系, 由图 7b可见大于直达波斜率的多次波同相轴。
|
图 7 带通滤波后的地震记录(a)及振幅谱归一化后的地震记录(b) |
利用(1)式恢复反射波, 然后分别采用常规方法、求和方法和相对位置方法生成虚源地震道。常规方法生成的虚源道集如图 8所示:在近偏移距虚源道集(图 8a)上, 虚源位置右侧可见恢复的双曲线型反射波, 且有残留面波, 但形态不完整, 它是由原始道集(图 6a)上残留的直达面波所致; 在中偏移距虚源道集(图 8b)上, 虚源右侧地震道的浅层反射波清晰可辨, 但观察不到中深层的反射波, 它被强线性同相轴和面波所掩盖, 信噪比低, 而虚源左侧地震道则基本是噪声; 在远偏移距的虚源道集(图 8c)上, 整个剖面被噪声占据。
|
图 8 采用常规方法生成的虚源道集 a 虚源位于地表第5道处; b 虚源位于地表第63道处; c 虚源位于地表第120道处 |
图 9为采用求和方法生成的虚源道集, 与采用常规方法生成的虚源道集相比可知, 二者虚源右侧地震道上的反射波几乎一样, 但图 9虚源左侧地震道上的反射波却得以再现, 这些反射波来自采用常规方法丢失的互相关负延迟信号。
|
图 9 采用求和方法生成的虚源道集 a 虚源位于地表第5道处; b 虚源位于地表第63道处; c 虚源位于地表第120道处 |
图 10为采用相对位置方法生成的虚源道集, 将其与采用求和方法生成的虚源道集(图 9)相比可知, 二者几乎一样, 但前者信噪比更高。
|
图 10 采用相对位置方法生成的虚源道集 a 虚源位于地表第5道处; b 虚源位于地表第63道处; c 虚源位于地表第120道处 |
由于监测工区地下的反射面基本水平, 为了进一步观察不同虚源地震道生成方法的效果, 将虚源道集按共偏移距进行叠加, 结果如图 11所示。与图 8、图 9和图 10的虚源道集相比可知, 图 11中恢复的反射波得到了有效增强, 虚源处的面波也相应地得到了增强, 掩盖了其下恢复的反射波, 且削弱了线性噪声。
|
图 11 采用常规方法(a)、求和方法(b)和相对位置方法(c)得到的虚源道集叠加产生的共偏移距叠加剖面 |
对采用上述方法恢复的反射波虚源道集进行速度谱分析并建立偏移速度模型, 然后进行动校叠加, 其结果如图 12所示。相较于采用求和方法得到的叠加剖面(图 12a), 采用相对位置方法得到的叠加剖面(图 12b)成像效果明显得到了改善, 且与同位置的三维地震资料时间偏移剖面(图 12c)有很强的对比性。虽然采用相对位置方法得到的剖面信噪比低, 面波污染严重, 但两者反射波同相轴对应关系好, 主要表现在300 ms、900 ms、1 100 ms和2 700 ms处, 浅层反射波同相轴连续性更好, 时间大于2 700 ms处同相轴连续性差。
|
图 12 采用求和方法(a)和相对位置方法(b)得到的反射波动校叠加剖面与三维地震资料时间偏移剖面(c) |
为验证恢复结果的可靠性, 分别对与测线3相邻的测线2和测线4的地震资料采用求和方法处理, 得到的虚源道集动校正后再叠加, 生成的剖面如图 13所示。测线3的剖面与测线4的剖面高度相似, 但测线2与测线3的剖面略有不同(注意:恢复反射波前, 除了带通滤波和均方根振幅归一化处理, 该地震资料未进行其它预处理)。
|
图 13 采用求和方法对测线2(a)、测线3(b)和测线4(c)的地震资料处理得到的虚源道集动校正后的叠加剖面 |
研究了反射波恢复技术, 并提出相应的虚源地震道生成方法, 即求和方法和相对位置方法。凹陷模型和实际资料试验结果表明, 利用本文提出的虚源道集生成方法, 均能有效地将丢失的互相关负延迟信号重新加以利用, 增加1倍的覆盖次数, 提高偏移成像质量。此外, 相对位置方法避免了求和方法引入噪声的弱点, 具有很强的实用性, 简单易行。
由于采用了均方根振幅归一化处理, 却并未在恢复反射波或地震道互相关后, 立即进行反归一化处理, 因此, 恢复的反射波振幅既不精确也不保真, 目前仅可用于地下构造成像解释。在恢复反射波前应去除残留面波, 在恢复反射波后应去除线性噪声, 否则将严重污染虚源道集的偏移成像质量。由于进行了互相关处理, 故虚源道集的地震子波是零相位的, 为解决虚源道集的地震子波问题, 可在互相关过程中采用所谓的自反褶积, 也可在偏移成像之前先进行最小相位化处理, 再进行脉冲反褶积。
致谢: 衷心感谢中国石化石油物探技术研究院的微地震项目组为我们提供实际地面微地震监测资料。| [1] |
WAPENAAR K, DRAGANOV D, SNIEDER R, et al. Tutorial on seismic interferometry:Part 1-Basic principles and applications[J]. Geophysics, 2010, 75(5): A195-A209. |
| [2] |
WAPENAAR K, DRAGANOV D, SNIEDER R, et al. Tutorial on seismic interferometry:Part 2-Underlying theory and new advances[J]. Geophysics, 2010, 75(5): A211-A227. |
| [3] |
SCHUSTER G T. Seismic Interferometry[M]. Cambridge: Cambridge University Press, 2009: 1-29.
|
| [4] |
CLAERBOUT J F. Synthesis of a layered medium from its acoustic transmission response[J]. Geophysics, 1968, 33(2): 264-269. |
| [5] |
RICKETT J, CLAERBOUT J. Acoustic daylight imaging via spectral factorization:Helioseismology and reservoir monitoring[J]. The Leading Edge, 1999, 18(8): 957-960. |
| [6] |
SCHUSTER G T, YU J, SHENG J, et al. Seismic interferometric/daylight imaging[J]. Expanded Abstracts of 63rd EAGE Annual Conference, 2001, 345-349. |
| [7] |
SCHUSTER G T, YU J, SHENG J, et al. Interferometric/daylight seismic imaging[J]. Geophysical Journal of the Royal Astronomical Society, 2004, 157(2): 838-852. |
| [8] |
HE R, HORNBY B, SCHUSTER G. 3D wave-equation interferometric migration of VSP free-surface multiples[J]. Geophysics, 2007, 72(5): S195-S203. |
| [9] |
JIANG Z, SHENG J, YU J, et al. Migration methods for imaging different-order multiples[J]. Geophysical Prospecting, 2007, 55(1): 1-19. |
| [10] |
WAPENAAR K. Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross correlations[J]. Physical Review Letters, 2004, 93: 254301-1-254301-4. |
| [11] |
WAPENAAR K, FOKKEMA J. Green function representations for seismic interferometry[J]. Geophysics, 2006, 71(4): SI33-SI46. |
| [12] |
WAPENAAR K, JOOST N, ELMER R. Seismic interferometry by crosscorrelation and by multidimensional deconvolution:a systematic comparison[J]. Geophysical Journal International, 2011, 185(3): 1335-1364. |
| [13] |
SNIEDER R, GRET A, DOUMA H, et al. Coda wave interferometry for estimating nonlinear behavior in seismic velocity[J]. Science, 2002, 295(5563): 2253-2255. |
| [14] |
SNIEDER R. Extracting the Green's function from the correlation of coda waves:A derivation based on stationary phase[J]. Physics Review E, 2004, 69(4): 046610-1-046610-8. |
| [15] |
SNIEDER R, WAPENAAR K, LARNER K. Spurious multiples in seismic interferometry of primaries[J]. Geophysics, 2006, 71(4): SI111-SI124. |
| [16] |
CALVERT R W, BAKULIN A, JONES T C. Virtual sources, a new way to remove overburden problems[J]. Expanded Abstracts of 66th EAGE Annual Conference, 2004, 234-238. |
| [17] |
BAKULIN A, CALVERT R. The virtual source method:Theory and case study[J]. Geophysics, 2006, 71(4): SI139-SI150. |
| [18] |
BAKULIN A, MATEEVA A. Virtual source applications to imaging and reservoir monitoring[J]. The Leading Edge, 2007, 26(6): 732-740. |
| [19] |
BAKULIN A, CALVERT R. Virtual source method:overview of history and development[J]. Expanded Abstracts of 78th Annual Internat SEG Mtg, 2008, 2726-2730. |
| [20] |
陈国金, 曹辉, 吴永栓, 等. 虚源法地震技术及数值模型试验[J]. 地球物理学进展, 2013, 28(5): 2725-2733. CHEN G J, CAO H, WU Y S, et al. The seismic virtual source method and numerical experiments[J]. Progress in Geophysics, 2013, 28(5): 2725-2733. |
| [21] |
ALEXANDROV D. Correcting for non-repeatable source signatures in 4D seismic with buried receivers:virtual source case study from Saudi Arabia[J]. Expanded Abstracts of 85th Annual Internat SEG Mtg, 2015, 2726-2730. |
| [22] |
LIN F C, RITZWOLLER M H, SNIEDER R. Eikonal tomography:Surface wave tomography by phase front tracking across a regional broadband seismic array[J]. Geophysical Journal International, 2009, 177(3): 1091-1110. |
| [23] |
KOICHI H. Application of active and passive surface-wave methods to shallow to deep S-wave velocity estimation[J]. Expanded Abstracts of 87th Annual Internat SEG Mtg, 2017, 5202-5206. |
| [24] |
WANG J, WU G, CHEN X. Frequency-Bessel transform method for effective imaging of higher-mode Rayleigh dispersion curves from ambient seismic noise data[J]. Journal of Geophysical Research:Solid Earth, 2019, 124(4): 201-211. |
| [25] |
CAMPILLO M, PAUL A. Long-range correlations in the diffuse seismic coda[J]. Science, 2003, 299(5606): 547-549. |
| [26] |
SHAPIRO N M. High-resolution surface wave tomography from ambient seismic noise[J]. Science, 2005, 307(5715): 1615-1618. |
| [27] |
房立华, 吴建平, 吕作勇. 华北地区基于噪声的瑞利面波群速度层析成像[J]. 地球物理学报, 2009, 52(3): 663-671. FANG L H, WU J P, LV Z Y. Rayleigh wave group velocity tomography from ambient seismic noise in North China[J]. Chinese Journal of Geophysics, 2009, 52(3): 663-671. |
| [28] |
LIN F C, LI D, CLAYTON R W, et al. High-resolution 3D shallow crustal structure in Long Beach, California:Application of ambient noise tomography on a dense seismic array[J]. Geophysics, 2013, 78(4): Q45-Q56. |
| [29] |
CHANG J P, RIDDER S A, BIONDI B L. High-frequency Rayleigh-wave tomography using traffic noise from Long Beach, California[J]. Geophysics, 2016, 81(2): B1-B11. |
| [30] |
陈国金, 郭建, 张亚红, 等. 基于环境噪声的地震响应重建方法及应用[J]. 石油物探, 2017, 56(6): 798-803. CHEN G J, GUO J, ZHANG Y H, et al. Technique for seismic response retrieval from ambient noise and its application[J]. Geophysical Prospecting for Petroleum, 2017, 56(6): 798-803. |
| [31] |
DRAGANOV D, WAPENAAR K, THORBECKE J. Seismic interferometry:Reconstructing the earth's reflection response[J]. Geophysics, 2006, 71(4): SI61-SI70. |
| [32] |
DRAGANOV D, WAPENAAR K, MULDER W, et al. Retrieval of reflections from seismic background-noise measurements[J]. Geophysical Research Letters, 2007, 34(4): L04305-01-L04305-4. |
| [33] |
DRAGANOV D, CAMPMAN X, JAN T, et al. Reflection images from ambient seismic noise[J]. Geophysics, 2009, 74(5): A63-A67. |
| [34] |
THORBECKE J, DRAGANOV D. Finite-difference modeling experiments for seiemic interferometry[J]. Expanded Abstracts of 80th Annual Internat SEG Mtg, 2010, 4013-4017. |
| [35] |
王德利, 程浩, 朱恒, 等. 基于反褶积算法的地震干涉技术被动源成像[J]. 吉林大学学报(地球科学版), 2012, 42(6): 1920-1926. WANG D L, CHEN H, ZHU H, et al. Deconvolution-Based seismic interferometry passive source imaging[J]. Journal of Jilin University(Earth Science Edition), 2012, 42(6): 1920-1926. |
| [36] |
张盼, 韩立国, 周岩, 等. 基于被动源多窗谱方法的主动源地震低频数据重构[J]. 应用地球物理, 2015, 12(4): 585-597. ZHANG P, HAN L G, ZHOU Y, et al. Passive-source multitaper-spectral method based low-frequency data reconstruction for active seismic sources[J]. Applied Geophysics, 2015, 12(4): 585-597. |
| [37] |
YILMAZ O. Seismic data analysis[M]. Tulsa: SEG Publishing, 2001: 144-145.
|