0
点赞
收藏
分享

微信扫一扫

基于ttcrpy(三维射线追踪)的跨孔CT高斯牛顿算法及python代码分享(1)

夹胡碰 2022-09-18 阅读 695

基于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环境下运行。

举报

相关推荐

0 条评论