欢迎Star!
PyFMM 是一个基于 Fast Marching/Sweeping Method 求解程函方程
PyFMM is a C/Python package for solving eikonal equation using Fast Marching/Sweeping Method, with examples and annotations.
PyFMM 目前可在Linux、Mac系统上运行,在Windows系统上可借助WSL兼容。
我还制作了一个简易图形界面 PyFMM-GUI 计算二维走时场,初学者可更好的理解射线追踪,也可更方便、直观地看到不同速度场下射线的扭曲形态。
我主要使用 PyFMM 计算地震波从震源出发在复杂介质中传播形成的初至波走时场, 并使用梯度下降获得满足费马原理的射线路径,故代码中的一些术语偏专业性。 类似的原理也可用于其它方面,如计算点到曲线/面的距离,或光学、电磁学等。
-
Python语言的便携、可扩展性与C语言的计算高效特点结合。 C程序被编译链接成动态库 libfmm.so ,PyFMM 再基于Python的 ctypes 标准库实现对C库函数的调用。再基于第三方库 NumPy、 SciPy 等可很方便地完成对C程序结果的数据整合;
-
C代码采取模块化编写,各功能分在不同代码文件中,方便移植到其它程序;
-
支持二维和三维情况;2D and 3D
-
支持直角坐标系和球坐标系;Cartesian and Spherical Coordinate
-
中文注释及示例;
为方便使用,我建立了在线文档,包括简易安装、API的介绍以及使用示例。
要求系统上有gcc
,make
,openmp
等开发工具,支持pip
命令一键安装:
pip install -v pyfmm-kit
对于不同系统以及安装问题,详见文档安装说明。
更多使用示例详见在线文档。
import pyfmm
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
pyfmm.logger.myLogger.setLevel('ERROR')
# 定义网格
nx, ny, nz = 401, 1, 101
xarr = np.linspace(0, 200, nx)
yarr = np.array([0.0])
zarr = np.linspace(0, 50, nz)
# 定义1D速度
vel1d = np.array([
[0.0, 3.2],
[5.0, 5.8],
[15.0, 6.5],
[30.0, 6.8],
[35.0, 8.1],
[80.0, 8.2]
])
# 插值1d分层速度
# _idxs = np.searchsorted(vel1d[:,0], zarr)
# velocity = vel1d[_idxs, 1]
# OR
# 插值1d梯度速度
velocity = interpolate.interpn((vel1d[:,0],), vel1d[:,1], zarr)
# 慢度数组
slowness = np.empty((nx, ny, nz))
slowness[...] = 1.0/velocity[None,None,:]
# 定义震源位置
srcloc = [0.0, 0.0, 0.0]
# 计算时间场
TT = pyfmm.travel_time_source(
srcloc,
xarr, yarr, zarr, slowness)
#====================================================================
# 绘制走时场和射线
fig, ax1 = plt.subplots(1, 1)
cs = ax1.contour(xarr, zarr, TT[:, 0, :].T, levels=30, linewidths=0.5)
ax1.clabel(cs)
for x in np.arange(5, 200, 5):
# 射线追踪
rcvloc = [x, 0, 0]
travt, rays = pyfmm.raytracing(
TT, srcloc, rcvloc, xarr, yarr, zarr, 0.1)
ax1.plot(rays[:,0], rays[:,2], c='r', lw=0.8, ls='--')
ax1.set_aspect('equal')
ax1.set_xlim(0, 200)
ax1.set_ylim(0, 50)
ax1.yaxis.set_inverted(True)
代码是我在研二写的,如果遇到bug,欢迎联系我([email protected]),我会完善! 也欢迎提出建议和更多示例!
基于PyFMM的体波走时反演以及面波反演后续也会开源。