在利用地震波研究地球深部结构时,经常需要绘制震相在深度剖面下的射线路径,同时也需要绘制地球内部的主要界面。

震相的射线路径可以用 TauP 提供的 taup_path 命令计算得到,然后在极坐标系 -JP 中绘制。难点在于如何绘制几个主要界面。

本示例将展示如何绘制震中距为30度的 PcP 和 PKiKP 震相的射线路径,同时绘制地球内的410、660界面以及 CMB 和 ICB。最终的绘图效果如下:

震相射线路径

绘图脚本如下:

#!/bin/bash
J=Pa10c/15
R=-10/40/0/6371
PS=earth-discontinuities.ps

gmt set MAP_GRID_PEN_PRIMARY 1p

# 绘制地表
gmt psbasemap -J$J -R$R -Byg6371 -Bs -K > $PS
# 绘制 410 界面
gmt psbasemap -J$J -R$R -Byg6371+5961 -Bs -K -O >> $PS
# 绘制 660 界面
gmt psbasemap -J$J -R$R -Byg6371+5711 -Bs -K -O >> $PS
# 绘制 CMB
gmt psbasemap -J$J -R$R -Byg6371+3480 -Bs -K -O >> $PS
# 绘制 ICB
gmt psbasemap -J$J -R$R -Byg6371+1221 -Bs -K -O >> $PS

# 计算并绘制震相的射线路径
taup_path -mod prem -ph PcP -h 300 -deg 30 -o PcP.raypath
taup_path -mod prem -ph PKiKP -h 300 -deg 30 -o PKiKP.raypath
gmt psxy PcP.raypath.gmt -W1p,blue -J$J -R$R -K -O >> $PS
gmt psxy PKiKP.raypath.gmt -W1p,red -J$J -R$R -K -O >> $PS

# 绘制震源和台站位置
gmt psxy -J$J -R$R -S -Gblack -N -K -O >> $PS << EOF
0 6071 0.4c a
30 6471 0.4c i
EOF

# 添加标注
gmt pstext -J$J -R$R -F+f11p+a -N -K -O >> $PS << EOF
-6 6170 20 Surface
-6 3280 20 CMB
-3 1020 20 ICB
16 4100 0 @;blue;PcP@;;
37 1600 0 @;red;PKiKP@;;
EOF

gmt psxy -J$J -R$R -T -O >> $PS
rm gmt.* *.raypath.gmt

脚本中使用了五个 psbasemap 命令分别绘制五个界面。

  • 第一个 psbasemap 用于绘制地表, -Byg6371 表明要在 Y 方向(极坐标下即 R 方向)以6371为间隔绘制网格线,由于地球半径是6371,所以理论上会在R=0和 R=6371 两处绘制网格线。 -Bs 的作用是只绘制 R=6371 处的网格线,且只绘制网格线而不绘制刻度;
  • 第二个 pabasemap 用于绘制410界面,其与前一命令基本相同,唯一的区别是 -Byg6371+5961 中多了 +5961,其作用是定义网格线的起算点,即此时将在 R=5961(即410界面)处绘制网格线。
  • 其他同理