接收函数

对于远震P波记录,将R分量反卷积Z分量即得到P波接收函数。P波接收函数是台站下方介质对远震入射 P波的介质响应。

以最简单的单层地壳模型为例,接收函数中主要的震相是直达 P,Ps, PpSs, PsPs 等。 一个典型的接收函数如图所示:

接收函数反映了台站下方的介质信息,因而可以用于约束台站下方结构。接收函数中可用的信息包括但不限于:

  1. 震相走时差
  2. 震相相对振幅
  3. 波形

基本原理

原理

台站处的地震记录可以表示为震源项、路径项和仪器响应的卷积:

$$ U_R(t) = S(t)*E(t)*E_R(t)*I(t) $$ $$ U_Z(t) = S(t)*E(t)*E_Z(t)*I(t) $$

其中 $U_R(t)$ 和 $U_Z(t)$ 分别表示地震记录的R分量和Z分量。 $S(t)$ 为震源项,$E(t)$ 表示从 震源到台站下方的介质响应。 $E_R(t)$ 和 $E_Z(t)$ 代表台站下方的介质响应,$I(t)$ 代表仪器响应。

假设 $E_Z(T) \approx \delta(t)$,则可得到:

$$U_Z(t) \approx S(t) * I(t)$$

由于时间域的卷积等于频率域相乘,因而得到

$$\frac{U_R(w)}{U_Z(w)}=\frac{S(w)E_R(w)E(w)I(w)}{S(w)E(w)I(w)}= E_R(w) $$

即将远震P波观测数据的R分量反卷积Z分量可得到台站下方的介质响应。

基本假设

接收函数的基本假设是台站下方介质响应的Z分量近似为脉冲函数。这一近似成立的条件是P波近垂直入射。

当P波近垂直入射时,由于P波偏振方向平行于波传播方向,故而Z分量中P震相振幅很大,而Ps 震相偏振方向垂直于波传播方向,故而Z分量Ps震相振幅很小。入射角越小,Z分量上的Ps震相 震相越小,因而将Z分量近似为脉冲函数。故而只有当震中距较大时才能满足这一近似条件,故而接收 函数要求震中距大于30度。

而当P波垂直入射时,此时P波到S波的转换系数近似为0,因而R分量将看不到Ps震相,无法 用于接收函数的研究,故而P波入射角也不能太小。实际情况下震中距上限设为95度。

基本流程

  1. 数据筛选
    1. Mw>5.5
    2. 震中距: 30-95
  2. 以P波走时作为参考,截取-30 s到120 s的波形
  3. 滤波至 0.05-1.5 Hz
  4. 旋转至RTZ分量或LQT分量
  5. R分量反卷积Z分量得到接收函数
  6. 利用接收函数反演地下结构

基本参数

利用 TauP 可以计算接收函数的一些基本参数:

$ taup_time -mod prem -ph P,P24.4s -h 100 -deg 36
$ taup_time -mod prem -ph P,P24.4s -h 100 -deg 98
$ taup_pierce -mod prem -ph P,P24.4s -h 100 -deg 36
$ taup_pierce -mod prem -ph P,P24.4s -h 100 -deg 98
震中距 Ps-P走时(s) Ps-P入射角(度) Ps-P 离台站距离(度) Ps-P 射线参数 (s/deg)
36度 414.37-411.03=3.34 14.17 vs 26.35 0.06 vs 0.12 8.507 vs 8.511
98度 806.64-803.46=3.17 7.33 vs 13.36 0.03 vs 0.06 4.430 vs 4.431

一些总结:

  1. Ps-P走时差随震中距变化很小,因而在简单情况下可以对接收函数做直接迭加
  2. Ps, P在Moho的穿透点与台站的水平距离在13km 之内 => 不同震中距的接收函数采样了台站下方半径13km的圆锥
  3. Ps-P走时差、P和Ps入射角以及Ps-P与台站的距离对震源深度不敏感 => 不同深度的接收函数可以直接叠加

反卷积

计算接收函数的核心步骤是反卷积。反卷积有多种算法。

频率域除法

时间域的反卷积等效于频率域的除法。因而,接收函数理论上可用如下公式计算:

$$ R(w)=\frac{U_R(w)}{U_Z(w)} $$

但频率域除法在遇到高频噪声以及频谱中低值时很不稳定,类似于数学中的除0。

water level 法

water level法的做法是人为给谱域除法加上一个常数,使得分母上的谱不小于这一常数值。

基本公式:

$$ r(\omega) = (1+c) \frac{R(\omega)S^*(\omega)}{|S(\omega)|^2+c\sigma_0^2} e^{-\frac{\omega^2}{4\alpha^2}} $$

几点说明:

  1. 引入 $S^*(\omega)$ 的目的是为了使分母变成实数,将复数除法变成乘法
  2. $e^{-\frac{\omega^2}{4\alpha^2}}$ 对得到的结果做高斯低通滤波以消除高频噪声的影响,根据 $\alpha$ 的值可以计算出滤波器的截止频率
  3. 为了防止除零或除以小值导致结果不稳定,加入了 $c\sigma_0^2$ 项,c的选取需要 根据频率域除法的结果以及经验判断,通常取0.1-0.5
  4. $1+c$ 的目的在于补偿由于water level造成的振幅损失

时间域迭代反卷积

Ligorría & Ammon, 1999, BSSA

将接收函数看做是一堆脉冲的叠加:

$$R(t)=Z(t)*RF(t)=Z(t)*\sum_m \delta(t_m)=\sum_mZ(t)*\delta(t_m)=\sum_m Z(t_m)$$

由此可以得到,R(t)可以认为一堆Z(t)的叠加。因而时间域迭代反卷积的做法是:

  1. 将Z(t)与R(t)做相关,相关系数最大的时间点对应两个接收函数中最大的脉冲之间的时间差
  2. 计算这个时刻的脉冲所对应的振幅,得到第一个接收函数
  3. 将此接收函数卷积Z分量,然后将其从R(t)中减去
  4. 重复步骤1-3,直到最终得到的接收函数和Z(t)的卷积与Z(t)相似为止

反演

震相相对走时直接反演法

根据地震学基本公式 $dt = p dx + \eta dz$ 以及射线参数保持不变这一特性,得到

$$ T = p X + \int \eta dz $$

其中 T 是震相走时,X 是震中距,p 是射线参数,$\eta$ 是垂直慢度。

进一步得到:

$$ T = p X + \int_h^{hmax} \eta dz + \int_{hmax}^{H} \eta dz + \int_{H}^0 \eta dz $$

其中 h 是震源深度,hmax 是射线最深穿透深度,H 是 Moho 深度。

下面以36度为例,分别讨论上面四项对 Ps-P 走时差的贡献:

  1. $p_1 * X - p_0 * X = (8.507-8.511) * X = -0.004*X=-0.15 s$
  2. $\eta_1 - \eta_2 = \sqrt{\frac{1}{Vp^2}-p_1^2} - \sqrt{\frac{1}{Vp^2}-p_1^2}$ 由此可以估算出第二项以及第三项的贡献之和大约是0.1 s

由此可知,第一、二、三项的贡献很小,可以忽略,进而得到:

$$ t_{Ps-P} = H * (\sqrt{\frac{1}{V_s^2}-p^2} - \sqrt{\frac{1}{V_p^2}-p^2}) $$

即,接收函数中 Ps-P 走时差与地壳内平均 Vp、Vs、射线参数 p 和地壳厚度 H 有关。

由此,可根据 Ps-P 走时差反演得到 Vp、Vs 和 H。

此方法的缺点在于:

  1. Ps 震相微弱,难以拾取
  2. Moho深度与Vp/Vs之间存在trade-off,即存在多个 H 与 Vp/Vs 的组合计算得到的 Ps-P 走时差与观测相同。

h-k 法

Zhu & Kanamori, 2000, JGR

如果只使用 Ps-P 走时差,由于 Moho深度与 Vp/Vs 存在 trade-off,所以无法约束介质结构。

hk 方法额外引入了 PsPs 和 PpSs+PsPs 震相。对于每个Moho深度 H 和 k 值,可计算理论的三个 走时残差,然后将三个走时所对应的振幅按照如下公式

$$ s(H,k) = w_1 r(t_1) + w_2 r(t_2) + w_3 r(t_3)$$

加权起来作为该点的值,最终得到 h-k 图。

此法适用于同一个台站多个地震的多个接收函数,对于每个接收函数,均可计算出 $s(H,k)$,然后对所有 计算得到的 $s(H,k)$ 做叠加即可。

共转换点叠加法(CCP)

基本原理:认为接收函数中的每个波形的振幅正比于特定深度处界面的波阻抗,因而可以将接收函数从时间 域反投回深度域。

具体操作如下:

  1. 根据背景速度模型计算出Ps转换波射线路径
  2. 将接收函数的振幅做入射角校正,并将振幅值放在 Ps 射线路径上对应的位置
  3. 将整个区域分成固定大小的网格,并计算每个网格内振幅的均值

缺点:

  1. 假定了一维速度模型
  2. 只能处理Ps震相
  3. 地壳多次反射波会被成像在更深的地方,因而CCP结果中会出现假界面。这一界面深度在100-250 km左右,恰好与LAB深度差不多,因而P波接收函数不能用于成像LAB

S接收函数

S波接收函数的原理与P接收函数的原理类似,其优点在于Moho、LAB、410和660界面的转换波均出现在S波到时的前面,而地壳多次波的到时在S波之后,因而S波接收函数中不会出现 地壳多次波的假界面现象,适合研究LAB。

缺点:

  1. S波位于P波的尾波内,噪声较大
  2. S波比P波更低频,因而分辨率更低
  3. S波到P波的转换波仅在一个小震中距范围内存在

数据筛选:

  1. S: 55-85°
  2. SKS: >85°
  3. ScS: 50-75°

接收函数的应用

  1. 地壳厚度变化
  2. 地壳内平均波速比变化
  3. 地壳地幔内界面(Moho、LAB、410、660等)
  4. 反演S波速度

可能的研究方向

  1. 更精确的反卷积方法
  2. 利用接收函数更精确地得到地下成像
  3. 消除浅部结构的多次波造成的假象
  4. 计算接收函数相对于模型的偏导,用于线性反演中
  5. 接收函数与其他方法的联合反演

参考文献

  1. Zhu & Kanamori, 2000, JGR
  2. Ligorría & Ammon, 1999, BSSA
  3. Zhu, 2000, EPSL
  4. Yuan et. al., 2006, GJI