20201120翻译_disba基于Python的面波正演模拟程序包

于 2020-12-16 20:52:56 发布
阅读量3.1k 收藏 24
点赞数 3
分类专栏： 地震学软件
原文链接：https://pypi.org/project/disba/
版权

disba 面波正演模拟程序包(Python)

项目描述

  disba是一个计算效率高的Python库，用于(正演)模拟面波频散，它实现了Computer Programs in Seismology (CPS) 程序代码的子集，这些代码是用Python用numba的just-in-time进行编译的。这样的实现减轻了其他基于CPS的库(例如pysurf96、srfpython和PyLayeredModel)通常需要的Fortran编译器的先决条件，这通常会导致进一步的安装故障排除，尤其是在Windows平台上。
  disba的目标是在不影响性能的前提下实现轻便化和可移植。例如，它与CPS用f2py编译的用于Rayleigh-wave的surf96程序的速度相似，但是随着层数的增加，(处理)Love-wave的速度要快得多。disba还实现了瑞雷波快速delta矩阵fast delta matrix算法，虽然慢得有点讽刺，但更稳健，且能较好处理了低速度带引起的相速度反转（的情况）。

特点

正演模拟：

  • 使用Dunkin矩阵 或 快速增量矩阵(fast delta matrix) 算法计算瑞雷波相/群速度频散曲线
  • 使用Thomson-Haskell方法计算勒夫波相/群速度频散曲线
  • 计算瑞雷波椭圆率

特征函数及敏感核：

  • 计算瑞雷波和勒夫波特征函数
  • 计算瑞雷波和勒夫波相/群速度；瑞雷波椭圆率相对于层厚、P波速度、S波速度及密度的敏感核

安装

推荐的安装disba及其所有依赖是通过Python包索引：

pip install disba[full] --user

另外，复制并提取disba包，然后在包所在位置运行：

pip install .[full] --user

用法

  • 1.下面的例子计算了瑞雷波/勒夫波前3阶(基阶、一阶、二阶)相速度频散曲线(PhaseDispersion)：
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#usage:The following example computes the Rayleigh- and Love- wave 
#phase velocity dispersion curves for the 3 first modes
#计算瑞雷波/勒夫波前三个阶次的相速度频散曲线(根据输入模型计算理论频散曲线)
import numpy as np
from disba import PhaseDispersion           #计算相速度频散
#from disba import GroupDispersion          #计算群速度频散
import matplotlib
from matplotlib import pyplot as plt

# Velocity model
# thickness, Vp, Vs, density
# km, km/s, km/s, g/cm3
#velocity_model = numpy.array([
#   [10.0, 7.00, 3.50, 2.00],
#   [10.0, 6.80, 3.40, 2.00],
#   [10.0, 7.00, 3.50, 2.00],
#   [10.0, 7.60, 3.80, 2.00],
#   [10.0, 8.40, 4.20, 2.00],
#   [10.0, 9.00, 4.50, 2.00],
#   [10.0, 9.40, 4.70, 2.00],
#   [10.0, 9.60, 4.80, 2.00],
#   [10.0, 9.50, 4.75, 2.00],
#])

##0.read velocity model读取速度模型
velocity_model=np.loadtxt('simple_mod.txt',skiprows=1)
#print('velocity_model=',velocity_model)

# Periods must be sorted starting with low periods
#t = numpy.logspace(0.0, 3.0, 100)
t = np.logspace(0.0, 3.0, 100)          #创建等比数列(周期)
#print('t=',t)

# Compute the 3 first Rayleigh- and Love- wave modal dispersion curves 计R/L频散曲线
# Fundamental mode corresponds to mode 0(基阶mode=0)
pd = PhaseDispersion(*velocity_model.T)
#print(pd)
#pd = GroupDispersion(*velocity_model.T)
cpr = [pd(t, mode=i, wave="rayleigh") for i in range(3)]
cpl = [pd(t, mode=i, wave="love") for i in range(3)]

##开始画图
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(12,6),sharey=False)

##1.Rayleigh Dispersion
ax1.set_xscale("log")
ax1.set_xlim(1,1000)
ax1.set_xlabel("Period[s]",fontsize=15)
ax1.set_ylim(3.2,4.8)
ax1.grid(axis="y")
ax1.set_ylabel("Phase Velocity [km/s]",fontsize=15)
ax1.plot(cpr[0][0],cpr[0][1],color="blue",linewidth=1,label="Fundamental")
ax1.plot(cpr[1][0],cpr[1][1],color="orange",linewidth=1,label="Mode 1")
ax1.plot(cpr[2][0],cpr[2][1],color="green",linewidth=1,label="Mode 2")
ax1.legend(loc="lower right",fontsize=10)
ax1.set_title("Rayleigh-wave",fontsize=20)
#print("cpr =",cpr)
#print("mode 0 cpr=",cpr[0])
#print("mode 1 cpr=",cpr[1])
#print("mode 2 cpr =",cpr[2])	
#print("mode 2 cpr period =",cpr[2][0])			#period
#print("mode 2 cpr velocity =",cpr[2][1])		#velocity
#print("mode 2 cpr mode =",cpr[2][2])			#mode
#print("mode 2 cpr wave =",cpr[2][3])			#wave
#print("mode 2 cpr type =",cpr[2].type)			#type
# pd returns a namedtuple (period, velocity, mode, wave, type)

##2.Love Dispersion
#print("cpl =",cpl)
ax2.set_xscale("log")
ax2.set_xlim(1,1000)
ax2.set_xlabel("Period[s]",fontsize=15)
ax2.set_ylim(3.2,4.8)
ax2.grid(axis="y")
ax2.set_ylabel("Love Velocity [km/s]",fontsize=15)
ax2.plot(cpl[0].period,cpl[0].velocity,color="blue",linewidth=1,label="Fundamental")
ax2.plot(cpl[1].period,cpl[1].velocity,color="orange",linewidth=1,label="Mode 1")
ax2.plot(cpl[2].period,cpl[2].velocity,color="green",linewidth=1,label="Mode 2")
ax2.legend(loc="lower right",fontsize=10)
ax2.set_title("Love-wave",fontsize=20)

##3.保存图片
#plt.subplot_tool()
#plt.show()
#plt.tight_layout()
plt.savefig('1syn_RL_phaseDisp.png',dpi=300)

R/L相速度频散曲线
同样的，GroupDispersion可以用于计算群速度(频散曲线)。

  • 2.disba的应用程序接口(API)在其所有类中保持一致，它们以同样的方式初始化及调用。因此，特征函数可以如下计算，特别注意resample函数(用于对输入模型进行重采样)的调用：
import numpy as np
from disba import EigenFunction
import matplotlib
from matplotlib import pyplot as plt

from disba._helpers import resample         #resample对模型进行重采样

# input velocity model
#velocity_model = np.array([
#   [10.0, 7.00, 3.50, 2.00],
#   [10.0, 6.80, 3.40, 2.00],
#   [10.0, 7.00, 3.50, 2.00],
#   [10.0, 7.60, 3.80, 2.00],
#   [10.0, 8.40, 4.20, 2.00],
#   [10.0, 9.00, 4.50, 2.00],
#   [10.0, 9.40, 4.70, 2.00],
#   [10.0, 9.60, 4.80, 2.00],
#   [10.0, 9.50, 4.75, 2.00],
#])
##0.读取速度模型
velocity_model=np.loadtxt('simple_mod.txt',skiprows=1)
#print("velocity_model1=",velocity_model)
velocity_model_thickness=velocity_model[0]
#print("velocity_model_thickness1=",velocity_model_thickness)
##0.1.对原始速度模型进行重采样
velocity_model=resample(*velocity_model.T,dz=0.5)
#print("velocity_model2=",velocity_model)
velocity_model_thickness=velocity_model[0]
#print("velocity_model_thickness2=",velocity_model_thickness)
velocity_model_vp=velocity_model[1]
#print("velocity_model_vp2=",velocity_model_vp)
velocity_model_vs=velocity_model[2]
#print("velocity_model_vs2=",velocity_model_vs)
velocity_model_rho=velocity_model[3]
#print("velocity_model_rho2=",velocity_model_rho)

#eigf = EigenFunction(*velocity_model.T)
#eigf = EigenFunction(velocity_model)
eigf = EigenFunction(velocity_model_thickness,velocity_model_vp,velocity_model_vs,velocity_model_rho)
eigr = eigf(20.0, mode=0, wave="rayleigh")      #基阶
#print(eigr)
eigl = eigf(20.0, mode=0, wave="love")          #基阶

# eigf returns a namedtuple
#  - (depth, ur, uz, tz, tr, period, mode) for Rayleigh-wave
#  - (depth, uu, tt, period, mode) for Love-wave

##开始画图
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(12,6),sharey=False)

##1.Rayleigh-wave eigenfunction
ax1.set_xlim(-2.0,2.0)
ax1.set_xlabel("Normalized eigenfunction",fontsize=15)
ax1.set_ylim(90,0)
ax1.set_ylabel("Depth(km)",fontsize=15)
ax1.grid(axis="y")
ax1.set_title("Rayleigh-wave",fontsize=20)

ax1.plot(eigr.ur,eigr.depth,color="blue",linewidth=1,label="UR")
ax1.plot(eigr.uz,eigr.depth,color="orange",linewidth=1,label="UZ")
ax1.plot(eigr.tz,eigr.depth,color="green",linewidth=1,label="TZ")
ax1.plot(eigr.tr,eigr.depth,color="red",linewidth=1,label="TR")
ax1.legend(loc="lower right",fontsize=10)

##2.Love-wave eigenfunction
ax2.set_xlim(-2.0,2.0)
ax2.set_xlabel("Normalized eigenfunction",fontsize=15)
ax2.set_ylim(90,0)
ax2.set_ylabel("Depth(km)",fontsize=15)
ax2.grid(axis="y")
ax2.set_title("Love-wave",fontsize=20)


ax2.plot(eigl.uu,eigr.depth,color="blue",linewidth=1,label="UU")
ax2.plot(eigl.tt,eigr.depth,color="orange",linewidth=1,label="TT")
ax2.legend(loc="lower right",fontsize=10)

#print("eigr.uz=",eigr.uz)
#print("eigr.depth=",eigr.depth)

#plt.show()
plt.savefig('2syn_RL_eig.png',dpi=300)

20s周期R/L特征函数

  • 3.相速度敏感核(PhaseSensitivity)(GroupSensitivity用于计算群速敏感核)：
import numpy as np
from disba import PhaseSensitivity
import matplotlib
from matplotlib import pyplot as plt

from disba._helpers import resample

##0.读取速度模型
#velocity_model = numpy.array([
#   [10.0, 7.00, 3.50, 2.00],
#   [10.0, 6.80, 3.40, 2.00],
#   [10.0, 7.00, 3.50, 2.00],
#   [10.0, 7.60, 3.80, 2.00],
#   [10.0, 8.40, 4.20, 2.00],
#   [10.0, 9.00, 4.50, 2.00],
#   [10.0, 9.40, 4.70, 2.00],
#   [10.0, 9.60, 4.80, 2.00],
#   [10.0, 9.50, 4.75, 2.00],
#])
#velocity_model=np.loadtxt('simple_mod.txt',skiprows=1)

velocity_model=np.loadtxt('simple_mod.txt',skiprows=1)
#print("velocity_model1=",velocity_model)
velocity_model_thickness=velocity_model[0]
#print("velocity_model_thickness1=",velocity_model_thickness)

velocity_model=resample(*velocity_model.T,dz=0.5)
#print("velocity_model2=",velocity_model)
velocity_model_thickness=velocity_model[0]
#print("velocity_model_thickness2=",velocity_model_thickness)
velocity_model_vp=velocity_model[1]
#print("velocity_model_vp2=",velocity_model_vp)
velocity_model_vs=velocity_model[2]
#print("velocity_model_vs2=",velocity_model_vs)
velocity_model_rho=velocity_model[3]
#print("velocity_model_rho2=",velocity_model_rho)

#ps = PhaseSensitivity(*velocity_model.T)
ps = PhaseSensitivity(velocity_model_thickness,velocity_model_vp,velocity_model_vs,velocity_model_rho)
parameters = ["thickness", "velocity_p", "velocity_s", "density"]
skr = [ps(20.0, mode=0, wave="rayleigh", parameter=parameter) for parameter in parameters]
#print("skr=",skr)
skl = [ps(20.0, mode=0, wave="love", parameter=parameter) for parameter in parameters]

# ps returns a namedtuple (depth, kernel, period, velocity, mode,wave, type, parameter)

##开始画图
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(12,6),sharey=False)

##1.Rayleigh-wave sensitivity kernels
ax1.set_xlim(-0.02,0.02)
ax1.set_xlabel("Sensitivity kernel",fontsize=15)
ax1.xaxis.get_major_formatter().set_powerlimits((0,1))
ax1.set_ylim(90,0)
ax1.set_ylabel("Depth[km]",fontsize=15)
ax1.grid(axis="y")
ax1.set_title("Rayleigh-wave",fontsize=20)
ax1.plot(skr[0].kernel,skr[0].depth,color="blue",linewidth=1,label=r'${\partial c}/{\partial d}$')
ax1.plot(skr[1].kernel,skr[1].depth,color="orange",linewidth=1,label=r'${\partial c}/{\partial \alpha}$')
ax1.plot(skr[2].kernel,skr[2].depth,color="green",linewidth=1,label=r'${\partial c}/{\partial \beta}$')
ax1.plot(skr[3].kernel,skr[3].depth,color="red",linewidth=1,label=r'${\partial c}/{\partial \rho}$')
ax1.legend(loc="lower right",fontsize=10)

##2.Love-wave sensitivity kernels
ax2.set_xlim(-0.02,0.02)
ax2.xaxis.get_major_formatter().set_powerlimits((0,1))
ax2.set_xlabel("Sensitivity kernel",fontsize=15)
ax2.set_ylim(90,0)
ax2.set_ylabel("Depth[km]",fontsize=15)
ax2.grid(axis="y")
ax2.set_title("Love-wave",fontsize=20)

ax2.plot(skl[0].kernel,skl[0].depth,color="blue",linewidth=1,label=r'${\partial c}/{\partial d}$')
ax2.plot(skl[1].kernel,skl[1].depth,color="orange",linewidth=1,label=r'${\partial c}/{\partial \alpha}$')
ax2.plot(skl[2].kernel,skl[2].depth,color="green",linewidth=1,label=r'${\partial c}/{\partial \beta}$')
ax2.plot(skl[3].kernel,skl[3].depth,color="red",linewidth=1,label=r'${\partial c}/{\partial \rho}$')
ax2.legend(loc="lower right",fontsize=10)

#plt.show()
plt.savefig('3syn_RL_phaseSensitivity.png',dpi=300)

20s周期R/L相速度敏感核

  • 4.椭圆率及椭圆率敏感核：
import numpy as np
from disba import Ellipticity, EllipticitySensitivity
import matplotlib
from matplotlib import pyplot as plt

from disba._helpers import resample

#velocity_model = numpy.array([
#   [10.0, 7.00, 3.50, 2.00],
#   [10.0, 6.80, 3.40, 2.00],
#   [10.0, 7.00, 3.50, 2.00],
#   [10.0, 7.60, 3.80, 2.00],
#   [10.0, 8.40, 4.20, 2.00],
#   [10.0, 9.00, 4.50, 2.00],
#   [10.0, 9.40, 4.70, 2.00],
#   [10.0, 9.60, 4.80, 2.00],
#   [10.0, 9.50, 4.75, 2.00],
#])
#velocity_model=np.loadtxt('simple_mod.txt',skiprows=1)

velocity_model=np.loadtxt('simple_mod.txt',skiprows=1)
#print("velocity_model1=",velocity_model)
velocity_model_thickness=velocity_model[0]
#print("velocity_model_thickness1=",velocity_model_thickness)

velocity_model=resample(*velocity_model.T,dz=0.5)
#print("velocity_model2=",velocity_model)
velocity_model_thickness=velocity_model[0]
#print("velocity_model_thickness2=",velocity_model_thickness)
velocity_model_vp=velocity_model[1]
#print("velocity_model_vp2=",velocity_model_vp)
velocity_model_vs=velocity_model[2]
#print("velocity_model_vs2=",velocity_model_vs)
velocity_model_rho=velocity_model[3]
#print("velocity_model_rho2=",velocity_model_rho)

t = np.logspace(0.0, 3.0, 100)

#ell = Ellipticity(*velocity_model.T)
ell = Ellipticity(velocity_model_thickness,velocity_model_vp,velocity_model_vs,velocity_model_rho)
rel = ell(t, mode=0)
#print("rel=",rel)

# ell returns a namedtuple (period, ellipticity, mode)
parameters = ["thickness", "velocity_p", "velocity_s", "density"]
#es = EllipticitySensitivity(*velocity_model.T)
es = EllipticitySensitivity(velocity_model_thickness,velocity_model_vp,velocity_model_vs,velocity_model_rho)
#print("es=",es)
ek = [es(20.0, mode=0, parameter=parameter) for parameter in parameters]
#print("ek=",ek)

# es returns a namedtuple (depth, kernel, period, velocity, mode,wave, type, parameter)

##开始画图
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(12,6),sharey=False)

##1.Ellipticity
ax1.set_xlim(1,1000)
ax1.set_xscale("log")
ax1.set_xlabel("Period[s]",fontsize=15)
#ax1.xaxis.get_major_formatter().set_powerlimits((0,0))
ax1.set_ylim(0.55,0.85)
ax1.set_ylabel("Ellipticity [H/V]",fontsize=15)
ax1.grid(axis="y")
ax1.set_title("Ellipticity",fontsize=20)
ax1.plot(rel.period,rel.ellipticity,linewidth=1,color="blue")

##2.
ax2.set_xlim(-0.02,0.02)
ax2.xaxis.get_major_formatter().set_powerlimits((0,1))
ax2.set_xlabel("Sensitivity kernel",fontsize=15)
ax2.set_ylim(90,0)
ax2.set_ylabel("Depth [km]",fontsize=15)
ax2.grid(axis="y")
ax2.set_title("Ellipticity sensitivity",fontsize=20)
ax2.plot(ek[0].kernel,ek[0].depth,color="blue",linewidth=1,label=r'${\partial \chi}/{\partial d}$')
ax2.plot(ek[1].kernel,ek[1].depth,color="orange",linewidth=1,label=r'${\partial \chi}/{\partial \alpha}$')
ax2.plot(ek[2].kernel,ek[2].depth,color="green",linewidth=1,label=r'${\partial \chi}/{\partial \beta}$')
ax2.plot(ek[3].kernel,ek[3].depth,color="red",linewidth=1,label=r'${\partial \chi}/{\partial \rho}$')

ax2.legend(loc="lower right",fontsize=10)

#plt.show()
plt.savefig('4syn_ellipticity_ellipticitySensKer.png',dpi=300)

不同周期瑞雷波椭圆率及20s周期瑞雷波椭圆率敏感核

ucaslilong
关注
  • 3
    点赞
  • 踩
  • 24
    收藏
    觉得还不错? 一键收藏
  • 8
    评论
专栏目录
评论 8
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值