obspy地震走时软件计算直达P波的理论到时（附带程序）

   **写在前面：**TAUP地震走时软件是计算地震波传播路径和走时的开源软件，可以用于解决地震波传播路径和速度模型的相关问题。

其中中文的参考网址：Tapu_SeisMan

软件安装(Python 版本)Obspy

ObsPy TauP
Pyrocko Manual另一个python版本的地学软件
一小时介绍ObsPy全英视频

pip install obspy

主要用到这个包：obspy.taup - Ray theoretical travel times and paths

初始化以及定义模型

from obspy.taup import TauPyModel
model = TauPyModel(model="iasp91")

obspy中提供的一维速度模型有

简单的近源是速度模型，感觉不同模型之间没啥区别，当然你也可以自己定义模型。我没有找到什么很简单的方法，甚至无法在python上实现，只能用JAVA上的taup上面计算。下面是我所有的笔记，做了3天，实在没做出来，只能用笨办法了。
使用Pyrocko Manual插件（python）计算P波到时（没做出来版本，中途改用taup了）

实际应用

如果你没有什么更改原始模型的需求，那么只是想简单的计算一下到时，可以用我这个function
(免责声明：我也是刚接触python所以用的真的不是很精，如果有更好的方法，希望各位未来院士可以分享给我)
这个只能计算直达波P波

from obspy.taup import TauPyModel

def  get_p_and_s_arrival_times(source_depth_in_km, distance):
    """
    **Theoretical p&s wave arrival time**
    reference page:https://docs.obspy.org/packages/obspy.taup.html
    Parameters：
    - model: 地震模型，如 "ak135" 或 "prem"
    - source_depth_in_km: event deepth（单位：千米）
    - distance_in_degree: epicenter（unit：degree）--use epicenter(km)/111
    Return：
    - P wave arrival time---travel_time_p,travel_time_s
    - S wave arrival time-travel_time_s
    """
    distance_in_degree = distance / 111

    # 创建 TauPyModel 实例
    model = TauPyModel(model="herrin")

    # 获取 P 波到时
    try:
        arrival_p = model.get_travel_times(source_depth_in_km=source_depth_in_km,
                                           distance_in_degree=distance_in_degree,
                                           phase_list=["p"])
    except ValueError:
        arrival_p = model.get_travel_times(source_depth_in_km=source_depth_in_km,
                                           distance_in_degree=distance_in_degree,
                                           phase_list=["P"])

    travel_time_p = round(arrival_p[0].time, 3)

    # 获取 S 波到时
    try:
        arrival_s = model.get_travel_times(source_depth_in_km=source_depth_in_km,
                                           distance_in_degree=distance_in_degree,
                                           phase_list=["s"])
    except ValueError:
        arrival_s = model.get_travel_times(source_depth_in_km=source_depth_in_km,
                                           distance_in_degree=distance_in_degree,
                                           phase_list=["S"])

    travel_time_s = round(arrival_s[0].time, 3)

    return travel_time_p, travel_time_s