TauP - teoretické časy příchodů a paprsky

pro 1-D modely Země


In [1]:
%pylab inline
from __future__ import print_function
import matplotlib.pylab as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 12, 8
import warnings
warnings.filterwarnings("ignore")

Populating the interactive namespace from numpy and matplotlib
In [2]:
from obspy.taup import TauPyModel

Nejdřív si načteme rychlostní model, TauP má předdefinováno několik 1-D rychlostních modelů:

  • 1066a, see [GilbertDziewonski1975]
  • 1066b, see [GilbertDziewonski1975]
  • ak135, see [KennetEngdahlBuland1995]
  • ak135f, see [KennetEngdahlBuland1995], [MontagnerKennett1995], and http://rses.anu.edu.au/seismology/ak135/ak135f.html (not supported)
  • ak135f_no_mud, ak135f with ak135 used above the 120-km discontinuity; see the SPECFEM3D_GLOBE manual at https://geodynamics.org/cig/software/specfem3d_globe/
  • herrin, see [Herrin1968]
  • iasp91, see [KennetEngdahl1991]
  • jb, see [JeffreysBullen1940]
  • prem, see [Dziewonski1981]
  • pwdk, see [WeberDavis1990]
  • sp6, see [MorelliDziewonski1993]

Anebo si člověk může sestavit vlatní rychlostní model, nejdřív do formátů .tvel nebo .nd a použije funkci build_taup_model().

In [3]:
model = TauPyModel(model="iasp91")

K použití jsou podobné funkce jako u klasického TauP:

  • get_travel_times
  • get_ray_paths
  • get_pierce_points

U variace s koncovkou _geo se místo distance zadávají dvojice (lat,lon).

In [4]:
arrivals = model.get_travel_times(source_depth_in_km=30, distance_in_degree=90, phase_list=['ttp+'])
print(arrivals)

7 arrivals
	P phase arrival at 776.509 seconds
	PcP phase arrival at 777.557 seconds
	pP phase arrival at 786.160 seconds
	sP phase arrival at 789.862 seconds
	PKiKP phase arrival at 1072.436 seconds
	pPKiKP phase arrival at 1082.368 seconds
	sPKiKP phase arrival at 1086.009 seconds

Objekt Arrivals obsahuje jednotlivé objekty Arrival s atributy:

  • phase Phase that generated this arrival
  • distance Actual distance in degrees
  • time Travel time in seconds
  • purist_dist Purist angular distance (great circle) in radians
  • ray_param Ray parameter in seconds per radians
  • name Phase name
  • purist_name Phase name changed for true depths
  • source_depth Source depth in kilometers
  • incident_angle Angle (in degrees) at which the ray arrives at the receiver
  • takeoff_angle Angle (in degrees) at which the ray leaves the source
  • pierce Points pierced by ray (computed by model.get_pierce_points())
  • path Path taken by ray (computed by model.get_ray_paths())
In [5]:
first = arrivals[0]
print(first.name)
first.__dict__

P
Out[5]:
{'distance': 90,
 'incident_angle': 14.004672319343223,
 'name': u'P',
 'path': None,
 'phase': <obspy.taup.seismic_phase.SeismicPhase at 0xb82a860>,
 'pierce': None,
 'purist_dist': 1.5707963267948966,
 'purist_name': u'P',
 'ray_param': 265.82560300719064,
 'ray_param_index': 235,
 'receiver_depth': 0.0,
 'source_depth': 30,
 'takeoff_angle': 15.812559069077675,
 'time': 776.5089537581711}

Zkratky pro předdefinované skupiny fází:

  • ttp = {p, P, Pn, Pdiff, PKP, PKiKP, PKIKP}
  • ttp+ = ttp + {PcP, pP, pPdiff, pPKP, pPKIKP, pPKiKP, sP, sPdiff, sPKP, sPKIKP, sPKiKP}
  • tts = {s, S, Sn, Sdiff, SKS, SKIKS}
  • tts+ = tts + {sS, sSdiff, sSKS, sSKIKS, ScS, pS, pSdiff, pSKS, pSKIKS}
  • ttbasic = ttp+ + tts+ + {ScP, SKP, SKIKP, PKKP, PKIKKIKP, SKKP, SKIKKIKP, PP, PKPPKP, PKIKPPKIKP}
  • ttall = ttbasic + {SKiKP, PcS, PKS, PKIKS, PKKS, PKIKKIKS, SKKS, SKIKKIKS, SKSSKS, SKIKSSKIKS, SS, SP, PS}

Můžeme vytvářet i jiné fáze, např. s využitím rozhraní v 410 km a 660 km.

Příklady:

  • P410s ... P se obrací pod 410, při zpětném průchodu skrz 410 konverze na S
  • Pv660P ... P se odráží od 660
  • 4kmps ... povrchová vlna s rychlostí 4 km/s

Více viz oddíl Phase naming v https://docs.obspy.org/packages/obspy.taup.html

Vykreslování paprsků



Vykreslovací funkce:

  • plot_travel_times() - časy příchodů vln v závisloti na epicentrální vzdálenosti
  • plot_ray_paths() - paprsky pro jeden zdroj a plno přijímačů
  • arrivals.plot_times() - časy příchodů vln v závisloti na jedné epicentrální vzdálenosti
  • arrivals.plot_rays()
    - paprsky pro daný zdroj a přijímač, zobrazení je buďto kartézské nebo sférické
    - ještě předtím nutno použít model.get_ray_paths()
In [6]:
from obspy.taup import plot_travel_times, plot_ray_paths
In [7]:
plt.rcParams['figure.figsize'] = 12, 8
plot_travel_times(source_depth=30, phase_list=['P','Pdiff','PcP','PKP','PKIKP','PKiKP', 'PP'], npoints=200)
Out[7]:
<matplotlib.axes._subplots.AxesSubplot at 0xb88ec88>

Přiřazení původních názvů PKP větví k názvům fází v ObsPy:

PKPab = PKP
PKPbc = PKP
PKPdf = PKIKP
PKPcd = PKiKP

In [8]:
plot_travel_times(source_depth=30, phase_list=['PKP','PKIKP','PKiKP'], npoints=200, 
                  min_degrees=130, max_degrees=160)
Out[8]:
<matplotlib.axes._subplots.AxesSubplot at 0xd7ebc18>
In [9]:
arrivals = model.get_ray_paths(source_depth_in_km=30, distance_in_degree=147, phase_list=['PKP','PKIKP','PKiKP'])
arrivals.plot_rays(plot_type='spherical', legend=True)
Out[9]:
<matplotlib.axes._subplots.PolarAxesSubplot at 0xec0c1d0>

Kartézské zobrazení a rozhraní 410 km a 660 km

In [14]:
phases = ['P']
phases = ['P', 'p^410P']
phases = ['P', 'p^410P', 'Pv660p^410Pv660P']
phases = ['P', 'p^410P', 'Pv660p^410Pv660P', 'sP']
phases = ['P', 'p^410P', 'Pv660p^410Pv660P', 'sP', 'sPP']

arrivals = model.get_ray_paths(source_depth_in_km=500, distance_in_degree=20, phase_list=phases)
arrivals.plot_rays(plot_type='cartesian', legend='lower right')
Out[14]:
<matplotlib.axes._subplots.AxesSubplot at 0x10fe4978>

Trocha umění na závěr

In [15]:
plot_ray_paths(source_depth=30, phase_list=['SKKS'], npoints=100, legend=True, plot_all=False)
Out[15]:
<matplotlib.axes._subplots.PolarAxesSubplot at 0x11f09d68>
In [16]:
plot_ray_paths(source_depth=30, phase_list=['PKJKS'], npoints=100, legend=True)
Out[16]:
<matplotlib.axes._subplots.PolarAxesSubplot at 0x14620d30>
In [17]:
plot_ray_paths(source_depth=30, phase_list=['PcPPcPPcP'], npoints=60, legend=True, plot_all=False)
Out[17]:
<matplotlib.axes._subplots.PolarAxesSubplot at 0x23e56780>
In [18]:
plot_ray_paths(source_depth=30, phase_list=['Pdiff', 'PP'], npoints=100, legend=True, plot_all=True)
Out[18]:
<matplotlib.axes._subplots.PolarAxesSubplot at 0x1734b7b8>