(1)梁作奎 盖增喜†
北京大学地球与空间科学学院地球物理学系, 北京 100871; † 通信作者, E-mail: zge@pku.edu.cn
摘要 为了在起伏自由边界条件下应用逆时偏移成像, 提出一种边界元与有限差分混合的方法来计算地震波场的逆时传播。首先, 使用边界元法反传地表台阵数据, 可以得到地下一水平界面上的波场位移数据; 然后, 将此水平面作为有限差分方法的边界, 并反传此水平面处的波形数据, 获得此水平面下空间的地震波场。为了应用成像条件, 利用地震波的性质, 将地震波场分离为矢量P波和S波成分, 对分离出的P波和S波分量, 应用互相关成像条件可获得地下界面。通过对理论数据的测试, 验证了方法的有效性。
关键词 起伏地形; 边界元法; 有限差分法; 被动源; 逆时偏移成像
接收函数方法1979年由Langston[1]提出, 逐渐成为获取台站下方地壳以及上地幔速度间断面的常规手段。接收函数主要利用坐标旋转和反卷积的方法, 从远震 P 波中分离出速度间断面产生的 Ps 震相, 并将转换波震相反投影到地下转换点[2], 获得地下界面深度。其中, 共转换点(common conversion point, CCP)叠加方法[3]是在地壳和上地幔的研究中应用最广泛的接收函数偏移成像方法之一[3-6]。
由于接收函数方法假设地下某一区域为水平界面, 因此, 当地下界面起伏较为平缓或者为倾斜界面时可以获得较好的结果。但是, 当台站下方存在复杂结构(比如陡峭界面)时, CCP 叠加方法不能准确成像[7-8]。另一方面, 当地表起伏较大时, 台阵数据受其影响, 地震波场变得非常复杂, 最终会影响CCP 叠加结果。通常使用的平滑窗口处理和不同地震事件结果叠加方法虽然可以获得较平滑的结果, 但结果的分辨率和准确性受到影响。
为了更深入地了解地下结构, 研究人员布设了密集的观测台阵, 以便获取更好的地下结构信息。由于存在上述问题, 接收函数方法不能充分利用密集台阵数据的优势, 需要提出新的方法来获得更精确的地下间断面形态。勘探地震学中使用的基于波动方程的逆时偏移成像方法, 能够充分利用密集台阵的波形数据, 获得精确的地下界面。其原理是对震源正演波场和接收点反传得到的波场应用成像条件, 获得地下界面。但是, 这种方法需要知道精确的震源时间函数，以便进行波场正演, 这在天然地震(被动源)研究中是很难满足的。
Shang 等[7]提出一种被动源叠前逆时偏移成像方法, 利用密集台阵, 在台阵下方存在复杂结构时仍能准确地获得间断面。这种方法使用有限差分法反传台阵波形数据, 获得地下空间波场快照, 然后从波场快照中分离出P分量和S分量, 并对这两个分量应用互相关类成像条件来获得地下界面。
有限差分法是最常用的一种正演模拟方法。这种方法将波动方程中波场函数的空间导数和时间导数用相应空间和时间的差分代替[9-11], 与有限元法和边界元法相比, 计算速度非常快。有限差分法一般用于边界为水平或垂直的情况, 很难用于不规则边界条件。如果使用有限差分法, Shang等[7]提出的方法在地表为水平状态或起伏可忽略时是可行的, 但在地表有较大起伏时, 显然不能直接应用。然而, 在断裂带等研究热点地区, 往往具有较大的地形起伏, 限制了这种方法的应用范围。尽管目前发展了一些可以处理复杂边界的有限差分方法[12], 但往往比较复杂, 很难应用于任意起伏地形的情况。
边界元法是一种用于层状介质地震波模拟的有效手段, 可以很容易地处理边界不规则的情况[13-14]。与有限元法和谱元法等方法相比, 边界元法在边界处具有非常高的精度, 并很容易处理震源为无穷远处的情况[15-16]。Ge 等[15-16]提出并改进了使用全局反透射矩阵传播算子的边界元法, 具有较高的准确性, 即使在多层模型下也具有较高效率。
为了在起伏地形情况下应用被动源逆时偏移成像方法, 本文结合边界元方法和有限差分方法的优点, 提出一种边界元与有限差分混合的方法进行波场重建, 以期获得台阵下方空间的地震波场。这一方法可以在保留有限差分方法计算灵活、快速优点的基础上, 考虑不规则起伏地形的作用。根据其偏振特性, 可以将重建后的地震波场分离, 得到矢量P 波分量和 S 波分量, 然后应用成像条件, 获得地下间断面形态。
当地震波入射到地球内部间断面时, 会发生透射和转换。当入射波为 P 波时, 在间断面上会产生透射的 P 波 Pp 和转换波 Ps, 这两种震相都会被地表的观测台站记录。因此, 可以利用这两种震相的到时差获取地球内部间断面的深度信息。接收函数方法可以将接收函数的到时直接转换为间断面的深度, 在地球内部间断面成像方面得到广泛的应用。但是, 接收函数方法假设 Pp 和 Ps 震相激发点的深度在同一水平面上。当地下间断面起伏较大, 甚至出现倾角较大的断层时, 这一假设不再适用。被动源叠前逆时偏移成像方法完全基于波动方程, 不对转换点深度做任何假设, 成像过程主要包括 3 个步骤: 波场重建、波场分离和应用成像条件。
被动源逆时偏移成像, 就是利用地表台站的观测数据重建地下区域的波场, 其实质就是求解如下边值问题:
(1)其中, u0(x, t)为地表台站记录的三分量地震波场。研究者们通常采用有限差分法求解这一问题[7-8,17]。受有限差分网格剖分和差分格式的影响, 通常要求接收台站位于同一水平界面上[17]。然而, 实际情况是地球上存在各种起伏地形, 因此限制了这一方法的应用。由于边界元方法能够很好地描述地表和内部间断面的起伏特征, 本文提出一种边界元与有限差分相结合的方法来进行波场重建。
已知起伏地表上的地震波场 u0(x1, x3(x1), t), 起伏地形通过
=
描述。地震波场的定理[8,15] 如下:
(2)其中, 
为位移场的格林函数, 
为应力的格林函数。z1 代表地球内部某一水平界面。当
为自由表面时, 满足自由边界条件, 
=0, 所以
。 (3)
因此, 我们可以利用地表的观测数据, 重建地球内部虚拟界面上的波场, 并以此作为有限差分法的边界条件, 重建虚拟界面下方区域的波场。
与地震勘探中主动源逆时偏移成像方法不同, 被动震源逆时偏移成像方法利用转换点处透射 P 波和转换 Ps 波的相关性进行成像, 因此在波场重建以后, 需要从波场中分离出 P 波和 S 波成分, 并对地下任一点处的P波和S波应用互相关成像条件。
由 Helmholtz 定理, 任意矢量场
可以表示为一个无旋度的矢量场与一个无散度的矢量场之和:
, (4)其中, 
, 且
。那么,
(5)
。 (6)若已知空间的波场, 则有
(7)
对于P波, 在波数域有
, (8)k1~k3分别为x1~x3方向对应的波数。
, (9)
故
。 (10)对于二维模型, 式(9)和(10)简化为
, (11)
。 (12)对于S波, 在波数域有
, (13)
, (14)故
, (15)
。 (16)对于二维模型, 式(15)和(16)简化为
, (17)
。 (18)获得P波和S波的波场后, 利用间断面上转换点处的透射 P 波与转换波 Ps 波的相关性, 对波场中任一点应用如下互相关条件成像。考虑到 P 波与S 波的波场传播方向相近, 其偏振方向接近垂直, 所以本文采用的成像条件为
, (19)其中, T表示任一点处P波或S波波场的总持续时间, x表示地下空间的位置。式(19)的离散形式为
。 (20)
其中, npts 表示时间采样点数, 也等于波场快照个数; Dt 表示时间采样间隔, 如果为常数, 可省略。
需要注意的是, 在震源时间函数比较复杂的情况下, 若将台站波形数据直接用于逆时偏移成像, 结果往往是多道平行于间断面的条纹, 不直观[7]。可以借鉴接收函数的思想, 将台阵记录的径向分量和垂直分量都对震源时间函数做反褶积后的结果作为逆时偏移成像方法的输入。
对每一个地震事件, 按垂直分量 P 波到时对齐, 求均值后得到视震源时间函数, 将台阵记录的垂直分量和径向分量对视震源时间函数做反褶积, 即可作为逆时偏移成像方法的输入。本文采用水准量反褶积方法, 具体流程如图1所示。
本文使用边界元法正演理论数据。震源时间函数是中心频率为 0.4Hz 的 Ricker 子波, 平面 P 波入射。使用的模型均为两层介质的二维模型, 上下层分别具有相同的 P 波速度、S 波速度和密度: 上层的 P 波和 S 波速度分别为 5.8 和 3.2km/s, 密度为2.6 g/cm3; 下层的P波和S波速度分别为8.1和4.5 km/s, 密度为3.4 g/cm3。模型长度均为300 km, 地形取自一个垂直于天山山脉走向的剖面, 地表均匀分布着水平间隔为0.5 km的601个台站。
图1 被动源逆时偏移成像流程
Fig. 1 Flow chart of passive source reverse time migration
图 2 显示具有不同间断面的 3 种模型, 其中曲线表示地表起伏, 虚线表示地下 2km 处的水平面, 黑色折线表示间断面。间断面深度在 40~50km 之间, 3种模型在水平方向 100 km附近存在断层。模型Ⅰ: 间断面倾角为 90°; 模型Ⅱ: 间断面倾角约63.4°; 模型Ⅲ: 间断面倾角约 45°。模型Ⅱ和模型Ⅲ在断层与水平间断面的连接处做了平滑处理。
图2 具有不同间断面的3种模型
Fig. 2 Three forward models with different discontinuities
通过正演得到地表台阵记录后, 在进行波场延拓之前, 将台阵记录对视震源时间函数做反褶积, 以去除震源特征影响。图 3 显示使用模型Ⅱ且地震波入射角为 20°时正演得到的台阵记录 Z 分量和 X分量以及台阵记录反褶积后的Z分量和X分量。
从图 3(a)和(b)可以看出: P 波和 S 波波形包含旁瓣; 由于地形影响, 波形比较复杂; Z分量和X分量都存在明显的P波, 而S波在Z分量不明显, 原因是入射角较小, 转换波能量主要分布在X方向。从图3(c)和(d)可以看出, 反褶积后P波和S波波形的旁瓣已消除。
利用式(3)对起伏地形上的地震记录进行波场重建, 可以得到地下给定的水平面处(图 2, 深度为2km)的波场位移数据。图 4 表示正演得到的水平面处地震记录反褶积后的 Z分量和 X 分量以及波场延拓后得到的水平面上波形的 Z分量和 X分量。通过对比可以看出, 重建后的波场(图 3)与理论模拟的波场(图 4)的两个分量的透射 P 波和转换Ps 波在波形和到时上都非常一致, 表明此水平面上波形得到准确的恢复。需要说明的是, 正演模型没有地表作为边界, 也没有地形影响, 所以没有地表反射波, 波形也比较简单。
绿色为背景色, 黄色为正值, 蓝色为负值。下同
图3 台阵记录的Z分量(a)和X分量(b)以及反褶积后 Z分量(c)和X分量(d)
Fig. 3 Z component (a) and X component (b) of array record, Z component (c) and X component (d) after deconvolution of array record
图4 水平面处的波形: 正演得到的 Z分量(a)和 X分量(b)以及波场延拓得到的 Z分量(c)和 X分量(d)
Fig. 4 Waveforms on the horizontal plane: Z component (a) and X component (b) from forward modeling, Z component (c) and X component (d) from wavefiled continuation
然后, 以水平面为边界, 使用时域有限差分法进行波场延拓, 即得到水平面下空间的波场。
对每一时刻的波场进行波场分离, 可以得到每个波场的 P 波成分和 S 波成分, 并最终得到每一点随时间变化的 P 波波形和 S 波波形。图 5 显示两 个时刻(t1=32 s, t2=40 s)波场快照的 Z 分量和 X 分量、波场分离后P波成分的Z分量和X分量以及 S波成分的 Z 分量和 X 分量。通过比较分离前后波场中 P 波与 S 波的位置可以看出, 波场被很好地分离成 P 波成分和 S 波成分。
考虑到分离后的波场中 P 波和 S 波波形在转换点处也应具有相关性, 如果仅对t1和t2时刻分离出的 P 波和 S 波分量应用成像条件, 由于转换点处同时产生透射 P 波和转换 S 波, 两者在间断面上具有很强的相关性, 由式(20)计算得到的成像结果在转换点处应具有较大的绝对值, 因此可以根据成像结果中的幅度确定间断面位置。图6显示t1和t2时刻成像结果, 可以看出, 图像较大值(黄色)位于间断面附近, 且不同时刻处于不同位置。
对分离后的整个波场中每一点的 P 波和 S 波波形应用成像条件, 结果如图 7 所示。图像上部空白区域是由于有限差分法只计算 2km 深度以下空间的波场所致。图 7 中成像结果幅度较大处为间断面位置, 可见逆时偏移成像结果能够准确地显示间断面的位置。成像结果在间断面的断层处出现位相变化, 产生这一现象的原因如图8所示。
(a1)和(b1) 波场快照的Z分量; (a2)和(b2) 波场快照的X分量; (a3)和(b3) 分离后P波成分的Z分量;(a4)和(b4) 分离后P波成分的X分量; (a5)和(b5) 分离后S波成分的Z分量; (a6)和(b6) 分离后S波成分的X分量; (a1)~(a6) t1=32 s; (b1)~(b6) t2=40 s
图5 两个波场快照(t1, t2)及分离后的波场
Fig. 5 Two wavefield snapshots (t1, t2) and results after wavefield separation
红线表示间断面位置, 下同
图6 t1和t2时刻单个波场快照的成像结果
Fig. 6 Images from single wavefield snapshot at t1 and t2
图7 逆时偏移成像结果
Fig. 7 Result of reverse time migration imaging
图 8 中黑色实线表示间断面位置, 间断面下方为远震 P 波(黑色箭头线, 箭头表示传播方向), 间断面上方为透射 P 波(黑色箭头线)和转换 S 波(灰色箭头线)。当透射 P 波振动方向与传播方向一致时, 在位置 A, 转换 S 波振动方向向右(灰色虚线箭头); 在位置 B, 转换 S 波振动方向向左。由于成像条件为 P 波与 S 波振动的叉乘, 在 A 与 B 处叉乘结果正负相反, 表现为成像结果的位相变化。成像结果中位相变化处的像也具有较大幅度, 因此具有较大幅度的像共同给出间断面的位置。由于地震波垂直或以很小角度入射间断面时产生转换 Ps 波振幅很小, 因此应用成像条件后, 不能根据成像结果的幅度识别出这部分间断面, 但这并不意味着间断面不连续。另外, 根据位相变化产生的原理, 也可将出现位相变化作为间断面倾角变化的证据。
图8 成像位相变化原因
Fig. 8 Reason for phase change when imaging
CCP 叠加方法是接收函数中提高成像质量的重要手段。与 CCP 叠加方法类似, 为了提高成像质量, 可以叠加不同地震事件逆时偏移(reverse time migration, RTM)成像结果。图 9~11 显示分别使用图 2 所示模型Ⅰ、模型Ⅱ和模型Ⅲ, 地震波到达间断面前传播方向(竖直向上方向, 顺时针为正)不同时, CCP叠加方法和逆时偏移成像方法的结果。从图 9~11 可以看出, 对于 CCP 叠加方法, 当地震波从左下方入射时, 断层处成像结果具有断开特性, 当地震波从右下方入射时, 断层处成像结果为较平缓的斜面, 都不能获得准确的断层位置; 多个结果的叠加也不能获得准确的断层位置, 反而容易在断层处表现为双层特性, 给结果的解释带来困难。对于逆时偏移成像方法, 地震波从不同角度入射时, 都能够获得准确的间断面; 叠加后的结果也给出准确的间断面位置。
本文使用边界元与时域有限差分混合的方法, 有效地解决了存在起伏地形时的波场重建问题。通过对理论合成数据的测试, 表明在台阵密度足够时, 与 CCP 叠加方法相比, 从波场中分离出 P 波和 S波并对两者应用互相关成像条件, 被动源逆时偏移成像方法可以获得更准确的地下间断面位置。
(a1)~(a4) 地震波传播方向分别为30°, 20°, −20°和−30°时的CCP叠加方法结果; (a5) 叠加所有单个结果后的CCP叠加方法结果; (b1)~(b4) 地震波传播方向分别为30°, 20°, −20°和−30°时, 被动源逆时偏移成像方法结果; (b5) 叠加所有单个结果后的被动源逆时偏移成像方法结果
图9 CCP叠加方法和逆时偏移成像方法对模型Ⅰ的成像结果比较
Fig. 9 Comparison the imaging results of model Ⅰ between CCP stacking and RTM
由于无法获得准确震源波场, 被动源逆时偏移成像方法只能先进行波场重建, 再对从波场分离出的 P 波和 S 波应用互相关成像条件。根据本文的测试, 被动源逆时偏移成像方法可以获得较准确的间断面位置, 而 CCP 叠加方法能较准确地获得水平间断面的位置, 但无法准确地获得大倾角间断面的位置。为了展示本文提出的方法, 我们有意忽略了应用实际地震数据时会遇到的两个问题。
一是台站间距问题。为了准确地恢复波场, 台站间距应小于最大频率对应波长的一半。这通常需要根据台站间距, 对台阵记录进行带通滤波。在台阵密度足够时, 使用更高频率的数据可以提高成像分辨率, 但对于CCP叠加方法, 这不能解决在断层处成像失真的问题。对于被动源逆时偏移成像方法, 在使用更高频率数据提高成像分辨率的同时, 能够更精确地成像。
(a1)~(a4) 地震波传播方向分别为30°, 20°, −20°和−30°时的CCP叠加方法结果; (a5) 叠加所有单个结果后的CCP叠加方法结果; (b1)~(b4) 地震波传播方向分别为30°, 20°, −20°和−30°时, 被动源逆时偏移成像方法结果; (b5) 叠加所有单个结果后的被动源逆时偏移成像方法结果
图10 CCP叠加方法和逆时偏移成像方法对模型Ⅱ的成像结果比较
Fig. 10 Comparison the imaging results of Model Ⅱ between CCP stacking and RTM
二是获得较准确的速度模型的问题, 对构造比较复杂的区域尤为重要。通常, 可以参考其他研究者的结果, 在必要时也可先通过层析成像获得一个速度模型。另外, 在地震波入射方向垂直于间断面的位置, 转换波能量几乎为零, 应用互相关成像条件不能获得此处位置, 但可以通过叠加不同方向和不同震中距地震事件的结果更好地成像。
(a1)~(a4) 地震波传播方向分别为30°, 20°, −20°和−30°时的CCP叠加方法结果; (a5) 叠加所有单个结果后的CCP叠加方法结果; (b1)~(b4) 地震波传播方向分别为30°, 20°, −20°和−30°时, 被动源逆时偏移成像方法结果; (b5) 叠加所有单个结果后的被动源逆时偏移成像方法结果
图11 CCP叠加方法和逆时偏移成像方法对模型Ⅲ的成像结果比较
Fig. 11 Comparison of the imaging results of Model Ⅲ between CCP stacking and RTM
参考文献
[1] Langston C A. Structure under Mount Rainier, Wa-shington, inferred from Teleseismic body waves. Jour-nal of Geophysical Research, 1979, 84(B9): 4749- 4762
[2] 王兴臣, 丁志峰, 朱露培. 浅析接收函数共转换 点叠加成像方法. 地球物理学进展, 2015, 30(6): 2814-2819
[3] Zhu Lupei. Crustal structure across of the San An-dreas Fault, southern California from teleseismic converted waves. Earth and Planetary Science Letters, 2000, 179(1): 183-190
[4] 武岩, 丁志峰, 朱露培. 利用共转换点叠加方法研究华北地区地壳结构. 地球物理学报, 2011, 54(10): 2528-2537
[5] 张广成, 吴庆举, 潘佳铁, 等. 利用 H-K 叠加方法和 CCP 叠加方法研究中国东北地区地壳结构与泊松比. 地球物理学报, 2013, 56(12): 4084-4094
[6] Nábelek J, Hetényi G, Vergne J, et al. Underplating in the Himalaya-Tibet collision zone revealed by the Hi-CLIMB experiment. Science, 2009, 325: 1371-1374
[7] Shang X, Hoop M V D, Hilst R D V D. Beyond receiver functions: passive source reverse time mig-ration and inverse scattering of converted waves. Geophysical Research Letters, 2012, 39(15): L15308
[8] Shang Xuefeng. Inverse scattering: theory and app-lication to the imaging of the earth’s seismic discon-tinuities [D]. Boston: MIT, 2014
[9] Martin R, Komatitsch D. An unsplit convolutional per-fectly matched layer technique improved at grazing incidence for the viscoelastic wave equation. Geo-physical Journal International, 2009, 179(1): 333–344
[10] Martin R, Komatitsch D, Gedney S D. A variational formulation of a stabilized unsplit convolutional per-fectly matched layer for the isotropic or anisotropic seismic wave equation. Computer Modeling in En-gineering and Sciences, 2008, 37(3): 274-304
[11] 冯英杰, 杨长春, 吴萍. 地震波有限差分模拟综述. 地球物理学进展, 2007, 22(2): 487-491
[12] Li H, Zhang W, Zhang Z, et al. Elastic wave finite-difference simulation using discontinuous curvilinear grid with non-uniform time step: two-dimensional case. Geophysical Journal International, 2015, 202(1): 102-118
[13] Kawase H, Aki K. A study on the response of a soft basin for incident S, P, and Rayleigh waves with special reference to the long duration observed in Mexico City. Bulletin of the Seismological Society of America, 1989, 79(5): 1361-1382
[14] Hartmann F, Zotemantel R. The direct boundary ele-ment method in plate bending. International Journal for Numerical Methods in Engineering, 2010, 23(11): 2049-2069
[15] Ge Zengxi, Chen Xiaofei. Wave propagation in ir-regularly layered elastic models: a boundary element approach with a global reflection/transmission ma- trix propagator. Bulletin of Scismological Society of America, 2007, 97(3): 1025-1031
[16] Ge Zengxi, Chen Xiaofei. An efficient approach for simulating wave propagation with the boundary ele-ment method in multilayered media with irregular interfaces. Bulletin of Scismological Society of Ame-rica, 2008, 98(6): 3007-3016
[17] Yan Rui, Xie Xiaobi. An angle-domain imaging con-dition for elastic reverse time migration and its appli-cation to angle gather extraction. Geophysics, 2012, 77(5): S105-S115
Passive Source Reverse Time Migration Method to Image Discontinuity Beneath Irregular Free-Surface
LIANG Zuokui, GE Zengxi†
Deparment of Geophysics, School of Earth and Space Sciences, Peking University, Beijing 100871; † Corresponding author, E-mail: zge@pku.edu.cn
Abstract Because the most used finite difference method to calculate wave propagation cannot handle irregular boundaries easily, a hybrid boundary element and finite difference method is proposed. Firstly, the boundary element method is used to obtain the wavefileds along a horizontal boundary which is going to be the boundary condition of the finite difference method in the next stop. Then the finite difference method is applied to obtain the wavefields under the interface. The wavefields are separated to vector P-wave and S-wave constituents, and considering the coherence of P-wave and S-wave in a conversion point, the image condition of cross-correlation between P- and S-waves is applied to obtain discontinuity image underground. Synthetic tests show the effecti-veness of the proposed method.
Key words irregular free-surface; boundary element method; finite difference method; passive source; reverse time migration
中图分类号 P315
doi: 10.13209/j.0479-8023.2017.110
国家自然科学基金(41374045, 41774047)和中国地质调查局地质调查项目(DD20160082)资助
收稿日期: 2017-05-21;
修回日期: 2017-06-06;
网络出版日期: 2017-06-08