基于ttcrpy(三维射线追踪)的跨孔CT高斯牛顿算法及python代码分享（1）
ttcrpy是加拿大学者伯纳德·吉鲁（Bernard Giroux）于2021年发布的开源python库，详见（https://github.com/groupeLIAMG），参考文献（Giroux B. 2021. ttcrpy: A Python package for traveltime computation and raytracing.
SoftwareX, vol. 16, 100834. doi:10.1016/j.softx.2021.100834
）。
ttcrpy库包含了三种射线追踪方法：快速扫描算法（FSM）、最短路径法（SPM）、动节点最短路径法（DSPM）。包含其二维与三维的实现。
ttcrpy库中给出了2D矩形网格和三角形网格、3D正六面体与四面体网格等网格剖分形式，对于非规则网格，要利用python中的vtk库和pygmsh库生成。
本博文借助ttcrpy中射线追踪算法，实现跨孔CT的高斯牛顿反演算法。
文章目录
一、ttcrpy正演
建立速度模型，高速体（4000m/s），低速体（2000m/s）,模型图如下：
利用ttcrpy对此模型进行射线追踪，得到一组初至走时数据：
射线路径：
二、ttcrpy反演
反演结果如下：
三、python代码分享
1、正演代码
# -*- coding: utf-8 -*-
"""
Created on Sat Sep 10 20:15:20 2022
@author: Shang'xiang
Stay Hungry Stay Foolish
"""
'''
此程序用于走时CT的正演
'''
import ttcrpy.rgrid as rg
import numpy as np
import matplotlib.pyplot as plt
# 创建网格
x = np.arange(0,31.0)
y = np.arange(0,21.0)
# 创建速度模型
v = 2000*np.ones((x.size-1,y.size-1))
v1 = 4000*np.ones((10, 4))
v[12:22,12:16] = v1
# 画模型图
fig, ax = plt.subplots()
cs = plt.pcolor(v,cmap='jet',edgecolor='w')
plt.colorbar()
plt.show()
# 给定发射点和接收点的坐标
nsrc = 30
nrcv = 30
srcs = np.ones(shape=(nsrc,2))
rcvs = 19*np.ones(shape=(nrcv,2))
for i in range(nsrc):
xsrc = 0.5 + 1*i
xrcv = 0.5 + 1*i
srcs[i][0] = xsrc
rcvs[i][0] = xrcv
# 离散网格
grid = rg.Grid2d(x, y, cell_slowness=True, method = 'SPM')
# 速度转换为慢度
slowness = 1./v
# tt_all = np.empty(0)
data_all = np.empty([900,5])
rays_all = list()
n = 0
for s in srcs:
src = np.array([s])
tt, rays = grid.raytrace(src, rcvs, slowness, return_rays=True)
data_all[30*n:30*(n+1),0:2] = np.repeat(src,30,axis = 0)
data_all[30*n:30*(n+1),2:4] = rcvs
data_all[30*n:30*(n+1),4] = tt
rays_all = rays_all + rays
n = n + 1
# 绘制走时图
plt.figure(2)
plt.plot(data_all[:,4], 'r-o')
plt.show()
# 保存走时
np.savetxt('G_FSM.txt', data_all)
plt.figure(3)
# 绘制射线路径图
for r in rays_all:
plt.plot(r[:,1],r[:,0],'b-',linewidth=0.25)
ax = plt.gca()
miloc = plt.MultipleLocator(1)
ax.xaxis.set_minor_locator(miloc)
plt.grid(linestyle = ':', color = 'black',which = 'both')
plt.xlim(0,20)
plt.ylim(0,31)
plt.show()
2、反演代码
# -*- coding: utf-8 -*-
"""
Created on Thu Aug 4 18:43:38 2022
@author: Shang'xiang
Stay Hungry Stay Foolish
"""
'''
此程序用于走时CT的GN反演
'''
# 初始化，读取数据
import ttcrpy.rgrid as rg
import numpy as np
from scipy.sparse.linalg import cg
import matplotlib.pyplot as plt
data = np.loadtxt('G_FSM.txt')
srcs = data[:,0:2]
rcvs = data[:,2:4]
tobs = data[:,4]
# 给定初始模型
# 创建网格
x = np.arange(0,31.0)
y = np.arange(0,21.0)
# 创建速度模型
v = 2000*np.ones((x.size-1,y.size-1))
# 离散网格
grid = rg.Grid2d(x, y, cell_slowness=True, method = 'SPM')
# 速度转换为慢度
slowness = 1./v
# =============================================================================
# # tcal,LL = grid.raytrace(srcs, rcvs, slowness, compute_L = True)
#
# # A = LL.todense()
#
# # I = np.eye(600)
#
# # zuo = A.T*A + 0.618*I
# # deltat = tobs - tcal
# # you = A.T.dot(deltat)
# # you = you.T
#
# # deltam,info = cg(zuo,you)
# =============================================================================
I = np.eye(600)
# =============================================================================
# # CT正演
# =============================================================================
def ctforward(slowness):
tcal,LL = grid.raytrace(srcs, rcvs, slowness, compute_L = True)
A = LL.todense()
zuo = A.T*A + 0.618*I
deltat = tobs - tcal
you = A.T.dot(deltat)
you = you.T
deltam,info = cg(zuo,you)
return deltam
# =============================================================================
# # 第一次计算
# # 计算模型修改量
deltam = ctforward(slowness)
deltam = deltam.reshape(30,20)
#
# # 修改慢度
slowness = slowness + deltam
# =============================================================================
# 开始循环
for i in range(10):
deltam = ctforward(slowness)
deltam = deltam.reshape(30,20)
slowness = slowness + deltam
v = 1./slowness
# 画模型图
fig, ax = plt.subplots()
cs = plt.pcolor(v,cmap='jet',edgecolor='w')
plt.colorbar()
plt.show()
注：此代码均在python spyder环境下运行。