基于ttcrpy的跨孔CT高斯牛顿算法及python代码分享（2）

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的跨孔CT高斯牛顿算法及python代码分享（2）

一、ttcrpy正演

建立速度模型，利用python中ttcrpy进行正演。

1、模型一

速度模型如下：
在这里插入图片描述
正演结果如下：

1.2 模型二

速度模型如下：
在这里插入图片描述
ttcrpy正演结果如下：

二、ttcrpy反演

拉普拉斯算子（Laplace Operator）是n维欧几里德空间中的一个二阶微分算子，定义为梯度（▽f）的散度（▽·f）。拉普拉斯算子也可以推广为定义在黎曼流形上的椭圆型算子，称为拉普拉斯-贝尔特拉米算子。

作为一个二阶微分算子，拉普拉斯算子把C函数映射到C函数，对于k≥2时成立。算子Δ :C? →C?，或更一般地，定义了一个算子Δ :C(Ω) →C(Ω)，对于任何开集Ω时成立。

在反演中加入拉普拉斯算子可以使反演解更加收敛，得到的结果接近真实模型，结果更光滑。

1、模型一反演结果

1.1 普通高斯牛顿法反演结果
在这里插入图片描述
1.2 拉普拉斯算子约束的高斯牛顿算法反演结果

2、模型二反演结果

2.1 普通高斯牛顿法反演结果
在这里插入图片描述
2.2 拉普拉斯约束的高斯牛顿反演结果

3、收敛性对比

将两种反演结果的百分误差绘制成如下图：
在这里插入图片描述
由图可见，带约束的高斯牛顿法反演收敛性更好，反演效果更好。

三、python代码分享

1、ttcrpy正演代码

# -*- 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((4, 4))

v[8:12,4:8] = v1
v[18:22,4:8] = 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('G03_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 Mon Sep 12 13:48:33 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('G03_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
    zuo = A.T*A + 0.618*C
    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环境中运行。