跳转至

反演

权重文件

在事件目录下,需要一个与之对应的权重文件 weight.dat。权重文件中指定了每个台站的5段波形数据在反演中的权重。

权重文件格式

权重文件的格式为:

station_name dist w1 w2 w3 w4 w5 tp ts
  • 第1列是台站名,需要与 SAC 文件名(后缀除外)一致
  • 第2列指定了该台站要使用哪个震中距的格林函数。gCAP 在反演时,会根据这里的震中距去找对应的格林函数文件 dist.grn.?,并将其与台站的实际观测数据进行对比。因而,理论上此处的 dist 应等于实际观测数据的震中距。但实际使用时为了避免每一个地震都要算一遍格林函数,因而 dist 只需要取格林函数库中与实际观测数据震中距最近的值即可。
  • w1w5 代表了 PnlZ、PnlR、Z、R、T 五段波形的权重,默认值为 1。可以根据各段波形的数据质量将权重设置为 0 到 1 中的任意值。
  • tp 是正数,则其表示 Pnl 波的到时
  • ts 是 S 波实际到时减理论到时

额外的说明:

  1. w2 值为 -1,就表示对应的台站是远震台,只使用 P (PnlZ) 和 SH (T) 波形段,此时若 ts 为正值则表示 S 波到时

生成权重文件

example/weight.pl 是一个用于生成权重文件的脚本。执行如下命令,则会生成权重文件 example/20080418093658/weight.dat

$ perl weight.pl 20080418093658

生成权重文件之前,可以手动观察一下数据质量。如果不想使用某数据,可以在相应的 Z 分量的头段 t9 上随意标一个到时。脚本在发现某数据的 t9 头段值不是 -12345 后,权重值会定为 0,从而在反演时不使用这个数据。

权重文件的内容是:

  IU_WCI      140 1 1 1 1 1  20.2   0.0
  NM_BLO      140 1 1 1 1 1  20.4   0.0
 NM_SIUC      145 1 1 1 1 1  21.2   0.0
  NM_SLM      210 1 1 1 1 1  29.0   0.0
  NM_FVM      235 1 1 1 1 1  31.8   0.0
  IU_WVT      260 1 1 1 1 1  35.0   0.0
 NM_PVMO      280 1 1 1 1 1  37.7   0.0
  IU_CCM      300 1 1 1 1 1  40.3   0.0
  NM_MPH      415 1 1 1 1 1  54.2   0.0

几点说明:

  1. 此权重文件为自动生成,实际情况中需要根据波形质量以及反演结果对其中的参数进行微调
  2. 示例数据中,每 5 km 震中距计算一个格林函数,因而权重文件中第2列都是5的倍数
  3. 第8列中 tp

    • 默认使用SAC头段 T1 的值(即用户手动标定的P波到时)
    • T1 未定义,则使用头段 T0 的值(即 TauP 标记的理论到时)
    • 若未安装 TauP 而导致头段 T0 未定义,则默认值为 0

时窗截取

CAP 会把完整的波形记录切割成两个时间窗口:Pnl 和面波窗口。时间窗口的定义方式如下:

  1. 若 SAC 文件 Z 分量的 t1t2 以及 T 分量的 t3t4 有定义,则将四个时间依次作为 Pnl 起点、终点和面波的起点、终点
  2. 若条件 1 不满足,且在反演时用 -V 选项给了视速度,则会用下面的公式计算 t1、t2、t3 和 t4:

    t1 = dist/vp - 0.2*m1, t2 = ts
    t3 = dist/vLove - 0.3*m2, t4 = dist/vRayleigh + 0.7*m2
    
  3. 若条件 1 和 2 均不满足,则使用下面的公式计算:

    t1 = tp - 0.2*m1,  t2 = ts
    t3 = ts - 0.3*m2,  t4 = t3+m2
    

说明:

  1. 以上几个公式的正确性尚待确认
  2. m1m2 决定了 Pnl 和面波时间窗口的最大长度,由 -T 选项控制
  3. 上述公式内的 dist 是震源位置到台站的直线距离,程序在实际计算时是等于 sqrt( distance^2 + depSqr),其中 distance 是震中距,而 depSqr 理应是震源深度的平方,但是程序实际是把 depSqr 固定为 25。尚不确定是否是BUG,但实际影响很小
  4. tpts 是用的格林函数文件中的值,而不是权重文件中的值

反演

执行反演

example/ 目录下执行命令:

$ cap.pl -H0.2 -P0.3 -S5/10/0 -T35/70 -D2/1/0.5 -C0.05/0.3/0.02/0.1 -W1 -X10 -Mmodel_15/5.0 20080418093658

即可得到反演结果。

结果解释

执行以上命令后,会看到如下输出:

20080418093658 15 1
Warning: flag=21 => the minimum 295.1/90.0/  1.7 is at boundary
inversion done

第二行的警告是因为 NM_MPH 台的位置很特殊,刚好在节面上,可以打开 sac,看看这个台的波形可以看到其 P 波非常弱。可以将权重文件中 NM_MPH 台的五个权重全部改为 0,警告就会消失。如果没有看到第三行 inversion done 这样的提示,说明并没有进行反演,多半是因为忘记生成了权重文件。

打开 example/20080418093658/model_15.ps,可以看到反演的结果。震源球右侧的三行文字分别表示:

  • 第一行给出了事件ID(20080418093658),反演使用的模型 model 以及震源深度 15 km
  • 机制解(FM)的参数是 295 90 2,即走向是 295 度,倾角是 90 度,滑移角是 2 度,说明这是一个走滑断层。震级为 5.20 Mw。残差为 1.881e-2

下面的每一排波形代表一个台站的拟合情况,波形图中红色的理论波形,黑色的是实际的波形。台站名下方的两个数字是震中距和偏移时间,如 IU_WCI 下面的 137.5/-1.89 表示震中距是 137.5 公里,波形向前移动了 1.89 秒,每一列震相下面的数字是拟合度和该震相的波形偏移。

有一点需要说明,也许你做出来的结果和我的有微小的差异。因为要写这个项目,这个事例,我做过多次,发现并不是每一次做出来的结果都一样,这一点,我尚不能解释,不过这种差异非常非常微小。

反演多个深度的结果

执行 example/inversion.pl 脚本即可反演多个深度的结果:

$ perl inversion.pl 20080418093658

执行 example/get_depth.pl 即可求出最佳深度和误差估计:

$ perl get_depth.pl 20080418093658