ACTA Scientiarum Naturalium Universitatis Pekinensis
起伏地形下地下间断面被动源叠前逆时偏移成像
梁作奎 盖增喜
北京大学地球与空间科学学院地球物理学系, 北京 100871; † 通信作者, E-mail: zge@pku.edu.cn
摘要 为了在起伏自由边界条件下应用逆时偏移成像, 提出一种边界元与有限差分混合的方法来计算地震波场的逆时传播。首先, 使用边界元法反传地表台阵数据, 可以得到地下一水平界面上的波场位移数据; 然后,将此水平面作为有限差分方法的边界, 并反传此水平面处的波形数据, 获得此水平面下空间的地震波场。为了应用成像条件, 利用地震波的性质, 将地震波场分离为矢量 P 波和 S 波成分, 对分离出的 P 波和 S 波分量,应用互相关成像条件可获得地下界面。通过对理论数据的测试, 验证了方法的有效性。关键词 起伏地形; 边界元法; 有限差分法; 被动源; 逆时偏移成像中图分类号 P315
CCP叠加结果。通常使用的平滑窗口处理和不同地震事件结果叠加方法虽然可以获得较平滑的结果,但结果的分辨率和准确性受到影响。
为了更深入地了解地下结构, 研究人员布设了密集的观测台阵, 以便获取更好的地下结构信息。由于存在上述问题, 接收函数方法不能充分利用密集台阵数据的优势, 需要提出新的方法来获得更精确的地下间断面形态。勘探地震学中使用的基于波动方程的逆时偏移成像方法, 能够充分利用密集台阵的波形数据, 获得精确的地下界面。其原理是对震源正演波场和接收点反传得到的波场应用成像条件, 获得地下界面。但是, 这种方法需要知道精确的震源时间函数，以便进行波场正演, 这在天然地震(被动源)研究中是很难满足的。
Shang 等[7]提出一种被动源叠前逆时偏移成像方法, 利用密集台阵, 在台阵下方存在复杂结构时仍能准确地获得间断面。这种方法使用有限差分法反传台阵波形数据, 获得地下空间波场快照, 然后从波场快照中分离出 P 分量和 S 分量, 并对这两个分量应用互相关类成像条件来获得地下界面。有限差分法是最常用的一种正演模拟方法。这种方法将波动方程中波场函数的空间导数和时间导数用相应空间和时间的差分代替[911], 与有限元法和边界元法相比, 计算速度非常快。有限差分法一般用于边界为水平或垂直的情况, 很难用于不规则边界条件。如果使用有限差分法, Shang 等[7]提出的方法在地表为水平状态或起伏可忽略时是可行的,但在地表有较大起伏时, 显然不能直接应用。然而, 在断裂带等研究热点地区, 往往具有较大的地形起伏, 限制了这种方法的应用范围。尽管目前发展了一些可以处理复杂边界的有限差分方法[12], 但往往比较复杂, 很难应用于任意起伏地形的情况。边界元法是一种用于层状介质地震波模拟的有效手段, 可以很容易地处理边界不规则的情况[1314]。与有限元法和谱元法等方法相比, 边界元法在边界处具有非常高的精度, 并很容易处理震源为无穷远处的情况[1516]。Ge 等[1516]提出并改进了使用全局反透射矩阵传播算子的边界元法, 具有较高的准确性, 即使在多层模型下也具有较高效率。为了在起伏地形情况下应用被动源逆时偏移成像方法, 本文结合边界元方法和有限差分方法的优点, 提出一种边界元与有限差分混合的方法进行波场重建, 以期获得台阵下方空间的地震波场。这一 方法可以在保留有限差分方法计算灵活、快速优点的基础上, 考虑不规则起伏地形的作用。根据其偏振特性, 可以将重建后的地震波场分离, 得到矢量P 波分量和 S 波分量, 然后应用成像条件, 获得地下间断面形态。
1 方法原理
当地震波入射到地球内部间断面时, 会发生透射和转换。当入射波为 P 波时, 在间断面上会产生透射的 P 波 Pp 和转换波 Ps, 这两种震相都会被地表的观测台站记录。因此, 可以利用这两种震相的到时差获取地球内部间断面的深度信息。接收函数方法可以将接收函数的到时直接转换为间断面的深度, 在地球内部间断面成像方面得到广泛的应用。但是, 接收函数方法假设 Pp 和 Ps 震相激发点的深度在同一水平面上。当地下间断面起伏较大, 甚至出现倾角较大的断层时, 这一假设不再适用。被动源叠前逆时偏移成像方法完全基于波动方程, 不对转换点深度做任何假设, 成像过程主要包括 3 个步骤: 波场重建、波场分离和应用成像条件。
1.1 波场重建
被动源逆时偏移成像, 就是利用地表台站的观测数据重建地下区域的波场, 其实质就是求解如下边值问题:
u ( x , t)
2 ( x , t ) 0, t2 (1) u ( x , t) u ( x , t ), 0其中, u0(x, t)为地表台站记录的三分量地震波场。研究者们通常采用有限差分法求解这一问题[78,17]。受有限差分网格剖分和差分格式的影响, 通常要求接收台站位于同一水平界面上[17]。然而, 实际情况是地球上存在各种起伏地形, 因此限制了这一方法的应用。由于边界元方法能够很好地描述地表和内部间断面的起伏特征, 本文提出一种边界元与有限差分相结合的方法来进行波场重建。已知起伏地表上的地震波场 u0(x1, x3(x1), t), 起
x3 x x1) [8,15]伏地形通过 = 3( 描述。地震波场的定理如下:
u ( x ,z 1) [ H ( x , x' )() u x' G ( x , x' ) t ( x' )]d ( x'), (2) 其中, G ( x , x')为位移场的格林函数, H ( x , x' ) 为应
即可作为逆时偏移成像方法的输入。本文采用水准量反褶积方法, 具体流程如图 1 所示。
2 数值计算
本文使用边界元法正演理论数据。震源时间函数是中心频率为0.4 Hz的 Ricker 子波, 平面P波入射。使用的模型均为两层介质的二维模型, 上下层分别具有相同的 P 波速度、S 波速度和密度: 上层的 P 波和 S 波速度分别为 5.8 和 3.2 km/s, 密度为2.6 g/cm3; 下层的 P 波和 S波速度分别为 8.1 和 4.5 km/s, 密度为 3.4 g/cm3。模型长度均为 300 km, 地形取自一个垂直于天山山脉走向的剖面, 地表均匀 分布着水平间隔为 0.5 km 的 601 个台站。
图 2 显示具有不同间断面的 3 种模型, 其中曲线表示地表起伏, 虚线表示地下 2 km 处的水平面,黑色折线表示间断面。间断面深度在 40~50 km 之间, 3种模型在水平方向 100 km附近存在断层。模型Ⅰ : 间断面倾角为 90°; 模型Ⅱ : 间断面倾角约63.4°; 模型Ⅲ : 间断面倾角约 45°。模型Ⅱ和模型Ⅲ在断层与水平间断面的连接处做了平滑处理。
2.1 波场重建
通过正演得到地表台阵记录后, 在进行波场延拓之前, 将台阵记录对视震源时间函数做反褶积,以去除震源特征影响。图 3显示使用模型Ⅱ且地震波入射角为 20°时正演得到的台阵记录 Z 分量和 X分量以及台阵记录反褶积后的Z 分量和 X分量。
从图 3(a)和(b)可以看出: P 波和 S 波波形包含旁瓣; 由于地形影响, 波形比较复杂; Z 分量和 X 分量都存在明显的 P 波, 而 S 波在 Z 分量不明显, 原因是入射角较小, 转换波能量主要分布在 X 方向。从图 3(c)和(d)可以看出, 反褶积后 P 波和 S 波波形的旁瓣已消除。
利用式(3)对起伏地形上的地震记录进行波场重建, 可以得到地下给定的水平面处(图2, 深度为2 km)的波场位移数据。图4表示正演得到的水平面处地震记录反褶积后的Z分量和X分量以及波场延拓后得到的水平面上波形的Z分量和 X 分量。通过对比可以看出, 重建后的波场(图 3)与理论模拟的波场(图 4)的两个分量的透射 P 波和转换Ps 波在波形和到时上都非常一致, 表明此水平面上
波形得到准确的恢复。需要说明的是, 正演模型没有地表作为边界, 也没有地形影响, 所以没有地表反射波, 波形也比较简单。
然后, 以水平面为边界, 使用时域有限差分法进行波场延拓, 即得到水平面下空间的波场。
2.2 波场分离与成像
对每一时刻的波场进行波场分离, 可以得到每个波场的 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 波