2023 7.10~7.16 周报 (RTM研究与正演的Python复现) (8.3更新)

2023-11-10 17:51

本文主要是介绍2023 7.10~7.16 周报 (RTM研究与正演的Python复现) (8.3更新),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

8.3更新部分我将用蓝色字体表示.

0 上周回顾

上周简单阅读了论文《Deep-Learning Full-Waveform Inversion Using Seismic Migration Images》, 但是并没读完…因为这篇论文中提到一个技术吸引了注意力: RTM (Reverse-time migration)
于是计划下周去专门熟悉熟悉RTM的机制, 并且试着用Python复现这个操作.
另外, 在复现之前, 我们还需要用Python实现正演相关的操作, 因为RTM的过程也是一种逆向的正演.

1 本周计划

  1. Python复现简单的声波正演.
  2. 利用简单的声波正演手段实现RTM.
  3. Python复现包含复杂的边界吸收条件的弹性波正演.
  4. 试着看: 复杂的弹性波正演可否实现RTM.

2 完成情况

2.1 Python复现简单的声波正演

首先定义库和正演需要的基本参数.

采样间隔(s)x和z轴的间隔(m)震源频率(Hz)时间采样个数震源震动采样个数震源水平位置(m)
0.00110252000121150

然后速度模型是一个 201 × 301 201 \times 301 201×301图像, 通过上述参数不难得知, 这是一个约为 2 2 2km * 3 3 3km 地下区域.
在这里插入图片描述
然后在一共 2000 2000 2000个时域采样点中, 我们在前 121 121 121个采样点中设置了雷克子波样的震源. 即 2 2 2s中的前 0.121 0.121 0.121s内释放震源.
在这里插入图片描述
震源水平位置是 150 150 150, 也就是说在这片工区的中心位置释放震动.

from __future__ import division
import numpy as np
from bruges.filters.wavelets import ricker
import matplotlib.pyplot as plt
from scipy.signal import convolve
import skimage.filters
import matplotlib
import scipy.io
matplotlib.use('TkAgg')dt = 0.001                          # 采样的时间域间隔
dx = 10                             # 采样的x轴间隔
dz = 10                             # 采样的z轴间隔
f0 = 25
nt = 2000
tmax = nt * dt                      # 最大时间采样点 (实际是0.0999)
t_array = np.arange(0, nt)wav1 = ricker(duration = 0.120, dt = dt, f = f0)[0]# 定义震源波形 (雷克子波)
vmodel = scipy.io.loadmat("./vmodel1182.mat")["vmodel"]
vmodel = vmodel.T
nx, nz = vmodel.shape
isx = 150
isz = 0

然后我们假设一个二维波场平面, 这个平面在计算机中被若干个距离均匀的点表示, 即交错网络. 然后 u ( x , z , t ) u(x,z,t) u(x,z,t)作为这个平面中任意点的应力值, 如果你把这个平面放在三维空间中 (第三维为振幅), 那么 u ( x , z , t ) ≠ 0 u(x,z,t)≠0 u(x,z,t)=0就表现为平面中的突出与凹陷, 即振幅非零. 而 u ( x , z , t k ) u(x,z,t_k) u(x,z,tk)就意味着波场在变化的过程中第 t k t_k tk时刻的景象, 或者说, 时间从 0 0 0时刻一直进行到 t k t_k tk的时候突然被定格, u ( x , z , t k ) u(x,z,t_k) u(x,z,tk)就是定格时被拍摄下来的快照 (snapshot)
网格的横向与纵向间隔分别被表示为 Δ x \Delta x Δx Δ z \Delta z Δz, 那么如果 x i x_i xi可表征 i Δ x i \Delta x iΔx; z j z_j zj可表征 j Δ z j \Delta z jΔz; t k t_k tk可表征 k Δ t k \Delta t kΔt, 那么 u ( x i , z j , t k ) u(x_i, z_j, t_k) u(xi,zj,tk)也可表示为 u i , j k u_{i,j}^k ui,jk.
现在我们有个二维的波动方程
∂ 2 u ( x i , z j , t k ) ∂ z 2 + ∂ 2 u ( x i , z j , t k ) ∂ x 2 = 1 v 2 ∂ 2 u ( x i , z j , t k ) ∂ t 2 (1) \frac{\partial^{2} u(x_i, z_j, t_k)}{\partial z^{2}}+\frac{\partial^{2} u(x_i, z_j, t_k)}{\partial x^{2}}=\frac{1}{v^{2}} \frac{\partial^{2} u(x_i, z_j, t_k)}{\partial t^{2}} \tag{1} z22u(xi,zj,tk)+x22u(xi,zj,tk)=v21t22u(xi,zj,tk)(1)
然后通过的泰勒的二阶展开的项进行差值计算得到二阶中心差分.
f ( x ) − f ( x − h ) = ∂ f ( x ) ∂ x h − 1 2 ∂ f 2 ( x ) ∂ x 2 h 2 + o ( h 2 ) (2) f(x)-f(x-h)=\frac{\partial f(x)}{\partial x} h-\frac{1}{2} \frac{\partial f^{2}(x)}{\partial x^{2}} h^{2}+o\left(h^{2}\right) \tag{2} f(x)f(xh)=xf(x)h21x2f2(x)h2+o(h2)(2)
f ( x + h ) − f ( x ) = ∂ f ( x ) ∂ x h + 1 2 ∂ f 2 ( x ) ∂ x 2 h 2 + o ( h 2 ) (3) f(x+h)-f(x)=\frac{\partial f(x)}{\partial x} h+\frac{1}{2} \frac{\partial f^{2}(x)}{\partial x^{2}} h^{2}+o\left(h^{2}\right) \tag{3} f(x+h)f(x)=xf(x)h+21x2f2(x)h2+o(h2)(3)
这里差值是3式对2式求差, 这样就获得了一般的二阶中心差分形式, 这些形式表现在对 u u u Δ t \Delta t Δt, Δ x \Delta x Δx Δ z \Delta z Δz的导数可以表示为:
∂ 2 u i , j k ∂ t 2 = u i , j k + 1 − 2 u i , j k + u i , j k − 1 Δ t 2 + O ( Δ t ) (4) \frac{\partial^{2} u^{k}_{i,j}}{\partial t^{2}}=\frac{u^{k+1}_{i,j}-2 u^{k}_{i,j}+u^{k-1}_{i,j}}{\Delta t^{2}} + O(\Delta t)\tag{4} t22ui,jk=Δt2ui,jk+12ui,jk+ui,jk1+O(Δt)(4)
∂ 2 u i , j k ∂ x 2 = u i + 1 , j k − 2 u i , j k + u i − 1 , j k Δ x 2 + O ( Δ x ) (5) \frac{\partial^{2} u^{k}_{i,j}}{\partial x^{2}}=\frac{u^{k}_{i+1,j}-2 u^{k}_{i,j}+u^{k}_{i-1,j}}{\Delta x^{2}} + O(\Delta x)\tag{5} x22ui,jk=Δx2ui+1,jk2ui,jk+ui1,jk+O(Δx)(5)
∂ 2 u i , j k ∂ z 2 = u i , j + 1 k − 2 u i , j k + u i , j − 1 k Δ z 2 + O ( Δ z ) (6) \frac{\partial^{2} u^{k}_{i,j}}{\partial z^{2}}=\frac{u^{k}_{i,j+1}-2 u^{k}_{i,j}+u^{k}_{i,j-1}}{\Delta z^{2}} + O(\Delta z) \tag{6} z22ui,jk=Δz2ui,j+1k2ui,jk+ui,j1k+O(Δz)(6)
4, 5和6式带入到1式中可以得到
1 v 2 u i , j k + 1 − 2 u i , j k + u i , j k − 1 Δ t 2 = u i + 1 , j k − 2 u i , j k + u i − 1 , j k Δ x 2 + u i , j + 1 k − 2 u i , j k + u i , j − 1 k Δ z 2 (7) \frac{1}{v^2} \frac{u^{k+1}_{i,j}-2 u^{k}_{i,j}+u^{k-1}_{i,j}}{\Delta t^{2}} =\frac{u^{k}_{i+1,j}-2 u^{k}_{i,j}+u^{k}_{i-1,j}}{\Delta x^{2}} +\frac{u^{k}_{i,j+1}-2 u^{k}_{i,j}+u^{k}_{i,j-1}}{\Delta z^{2}} \tag{7} v21Δt2ui,jk+12ui,jk+ui,jk1=Δx2ui+1,jk2ui,jk+ui1,jk+Δz2ui,j+1k2ui,jk+ui,j1k(7)
稍微整理下, 我们希望知道下个时刻的波形结构, 因此我们只把这个式子中唯一的 u i , j k + 1 u_{i,j}^{k+1} ui,jk+1整理到左边, 其余全部往右侧移 (因为我们的交错网络每个点的分布是相互均匀的, Δ x = Δ z \Delta x = \Delta z Δx=Δz, 于是得用将 Δ z \Delta z Δz Δ x \Delta x Δx代替) . 得到:
u i , j k + 1 = r 2 ( u i + 1 , j k + u i , j + 1 k − 4 u i , j k + u i − 1 , j k + u i , j − 1 k ) + 2 u i , j k − u i , j k − 1 (8) u^{k+1}_{i,j} = r^2 \left( u^{k}_{i+1,j} + u^{k}_{i,j+1} - 4 u^{k}_{i,j} + u^{k}_{i-1,j} + u^{k}_{i,j-1}\right) + 2 u^{k}_{i,j} - u^{k-1}_{i,j} \tag{8} ui,jk+1=r2(ui+1,jk+ui,j+1k4ui,jk+ui1,jk+ui,j1k)+2ui,jkui,jk1(8)
这里 r = v Δ t Δ x r = \frac{v \Delta t}{\Delta x} r=ΔxvΔt.
同时, 8式括号里面的 i i i j j j是针对整个应力平面 u k + 1 u^{k+1} uk+1的更新一部分, 而这部分刚好可以通过一个卷积 D D D等价. 公式如下
u k + 1 = r 2 ( u k ∗ D ) + 2 u k − u k − 1 , D = ( 0 1 0 1 − 4 1 0 1 0 ) (9) u^{k+1} = r^2 \left( u^{k} * D\right) + 2 u^{k} - u^{k-1}, D= \begin{pmatrix}0& 1& 0\\1& -4& 1\\0& 1& 0\end{pmatrix}\tag{9} uk+1=r2(ukD)+2ukuk1,D= 010141010 (9)
这个公式里面, 我们把针对 i i i j j j的循环囊括了进来, 因此所有操作都是针对矩阵本身进行的. 这个过程可以通过下述代码实现

def solve_fd2d_withabc(v, w, vmodel, r, dt, dx, dz):'''Compute wave amplitude at the next k-th time stepwith boundary conditions:param v:           Snapshot of amplitude at last step (k-1):param w:           Snapshot of amplitude at previous step (k-2).:param vmodel:      Velocity model background:param r:           Mixing parameters:param dt:          Time sampling interval:param dx:          x-axis sampling interval:param dz:          z-axis sampling interval:return:'''D = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])u_out = 2 * v - 1 * w + (r ** 2) * convolve(v, D, mode='same')# bottomu_out[:, -1] = v[:, -2] + (r[:, -1] - 1) / (r[:, -1] + 1) * (u_out[:, -2] - v[:, -1])# rightu_out[-1, :] = v[-2, :] + (r[-1, :] - 1) / (r[-1, :] + 1) * (u_out[-2, :] - v[-1, :])# leftu_out[0, :] = v[1, :] + (r[0, :] - 1) / (r[0, :] + 1) * (u_out[1, :] - v[0, :])return u_out

后面三行是边界吸收代码 (低配版), 这个我在别处抄来的. 边界吸收算法可以很复杂, 这里算是超级简单的一种.
然后迭代代码:

u_abc = np.zeros((nx, nz, nt), dtype = float)
r = vmodel*dt/dx           
# 从第一个波场开始进行正演, 更新时间序列中每个节点的波场快照
for k in t_array:                                   # 一个时间采样点一个采样点地尝试if k >= 2:                                      # 需从步骤2开始, 因为步骤1和步骤0是初始条件v = u_abc[:, :, k - 1]w = u_abc[:, :, k - 2]u = solve_fd2d_withabc(v, w, vmodel, r, dt, dx, dz)if k < len(wav1):                           # 如果源处于活动状态, 则将其幅度添加到波场# 源是否活跃取决于雷克子波的影响时间u[isx, isz] = u[isx, isz] + wav1[k]u_abc[:, :, k] = u

完成迭代后, u_abc就存储不同时刻的波场快照了, 我们可以看到不同时刻的波场快照:
在这里插入图片描述
如果这个时候在图像的第一行设置一个系列检波器, 记录下每个检波器都在 2000 2000 2000个采样时刻下的波动情况, 那么我们就可以获得常规地震工区打炮后获得的地震数据图:

surface_record_raw = u_abc[:,1,:]

在这里插入图片描述
因为我们用的边界吸收条件是比较简易的, 因此这个地震数据存在很多不稳定的反射, 衍射和波形的 “尾流” (应该可以这么说吧…我感觉是挺像尾流的, 就是前向波的末端有抖动).

2.2 利用简单的声波正演手段实现RTM

什么是逆时偏移 (Reverse-time migration: RTM)?
我们这里使用的是叠前逆时偏移, 即对每个单炮记录进行逆时偏移, 然后对于各炮成像结果进行叠加.
单炮记录怎么进行逆时偏移?
它将采样的最后时刻的波场快照 u ( x , z , T ) u(x,z,T) u(x,z,T)消除前向波, 并以之作为起始平面, 按照时间反推, 并且令 u ( x , z = 0 , T ) u(x,z=0,T) u(x,z=0,T)作为震源进行正演材料得到 u ∗ ( x , z , 0 ) u^*(x,z,0) u(x,z,0)时刻的波场快照. 然后呢, 令 u ( x , z = 0 , T − 1 ) u(x,z=0,T-1) u(x,z=0,T1)作为震源进行正演材料得到 u ∗ ( x , z , 1 ) u^*(x,z,1) u(x,z,1)时刻的波场快照, 之后, 令 u ( x , z = 0 , T − 2 ) u(x,z=0,T-2) u(x,z=0,T2)作为震源进行正演材料得到 u ∗ ( x , z , 2 ) u^*(x,z,2) u(x,z,2)时刻的波场快照. 这就是一个Reverse-time的过程.

"令 u ( x , z = 0 , T ) u(x,z=0,T) u(x,z=0,T)作为震源进行正演材料"的含义有点两点:

  1. 震源不再是一个点, 而是横向一条线.
  2. 震源不再是长度为 t < T t<T t<T的雷克子波, 而是持续时间为 T T T的震源

这里, 可以简化理解一个点是: RTM的震源在 T T T时间内形成的二维波形图像其实就是 走时信息静默的 正演单炮山峰波形的上下翻转图.

这里放出代码来辅助理解 (说起来如此繁琐的东西, 用代码非常简单就表示了, 这就是代码的强大):
这里的逆时正演得到的波形图我们简称为上行波, 即 u ∗ u^* u, 为啥是上行波呢? 因为之前正演的波是向下传递的过程中接收器接受生成的, 而RTM震源是利用的接收器信息的逆过程.

vmodel_smooth = skimage.filters.gaussian(vmodel, 30)    # 平滑版的速度模型
r_smooth = vmodel_smooth*dt/dx                          # 得到平滑版的r
u_up = np.zeros((nx, nz, nt), dtype=float)
for k, tk in enumerate(t_array):if k >= 2:v = u_up[:, :, k - 1]w = u_up[:, :, k - 2]u = solve_fd2d_withabc(v, w, vmodel_smooth, r_smooth, dt, dx, dz)# 上行波是, 整个平面都提供震动源, 而非单一的点震源. 这个波是原本的反射波: [:, -1]->[:, -2]u[:, 0] = u[:, 0] + muted_gather[:, -(k - 1)]u_up[:, :, k] = u

强调!
这里进行上行波模拟的时候, 我基于的是模糊化后的速度模型来进行反向正演的! 这里是有原因的, 因为一般正演可以在现实世界中通过打炮检测来模拟出来 (毕竟震源可控, 只有单个炮点), 但是上行波的模拟很难实现 (震源众多, 波形难以控制). 因此这里采用了高斯平滑的速度模型来模拟人工猜测的低准确度的速度模型

这里的muted_gather是消除前向波, 也就是走时 (Traveltime) 信息的单炮地震图像.

什么是走时信息呢? 走时信息就是震源释放出来的最前面的波, 这个波没有任何反射内容, 就只是震源波的代表. 当然如果完全用前向波代替也不是很恰当, 因为向下传递的波也可以叫做前向波, 但是我们有个限定, 即这个走时信息的波是应当是沿着大地进行横向传播的前向波. 为何? 因为这个波形是每个接收器介绍到的第一个信息, 也就是左图中X[0, :]的每个点收到的第一次波动信息. 我们的单炮地震图像 (右图), 即共炮接受集合 (common shot gathers), 正是左图中X[0, :]的每个点在 T T T采样时间内接受的信号的绘制. 而走时信息就是这个图像的山峰的表面外层, 这个最清晰的部分
请添加图片描述要想消除这个信息不难, 因为由GIF图可知, 只要知道走时信息什么时候穿过来和震源有效波长时间, 然后把每个接收器的这段时间接受的信号设置为0 (静默) 就好了.
在这里插入图片描述

而走时信息之所以叫做"走时 (Traveltime) "信息, 其实就是波走 (Travel) 到一个位置的最短时间(time), 那么用路径除以时间就好.
代码如下:

# 前向波 (走时信息) 静默
muted_gather = surface_record_raw.copy()
x_array = np.arange(0, nx*dx, dx)                   # 基于采样进行标点, 得到地表处每个标点的数组 [0米, dx米, 2dx米, 3dx米, ... , (nx-1)dx米]
v0 = vmodel[:,0]                                    # v0表示地表一层的默认地层速度值 (长度为nx的数组)
# np.cumsum(dx/v0)表示在地表第一层区域内, 波从左到右传播的用时数组 [dx/v0秒, 2dx/v0秒, 3dx/v0秒 ...]
# ...[isx]自然是指的传到震源激发点的用时
# 相减后就得到中间出发向左右两边传递的波在不同位置的时间 [... -2dx/v0秒, -dx/v0秒, 0秒, dx/v0秒, 2dx/v0秒 ...]
traveltimes = abs(np.cumsum(dx/v0) - np.cumsum(dx/v0)[isx])
for traceno in range(len(x_array)):muted_gather[traceno, 0:int(traveltimes[traceno]/dt + len(wav1))] = 0   # 最后加上len(wav1)是雷克子波的时间, 因为雷克子波能量完全释放需要时间

请将这个代码放在上行波代码的前面. 最终得到静默走时信息的图像:
在这里插入图片描述
最后, RTM图像的获取需要上行波场和下行波的波场快照进行时间域上的互相关.
下行波的获取和常规正演一致, 只不过做RTM的时候我们对于完整速度模型是未知的, 因此这时正演的背景就是初始速度模型.
这种互相关是在时间域上互反的, 用代码来解释就是:

u_down = np.zeros((nx, nz, nt), dtype=float)
for k, tk in enumerate(t_array):if k >= 2:v = u_down[:, :, k - 1]w = u_down[:, :, k - 2]u = solve_fd2d_withabc(v, w, vmodel_smooth, r, dt, dx, dz)if k < len(wav1):u[isx, isz] = u[isx, isz] + wav1[k]u_down[:, :, k] = umigrated_image = np.zeros_like(u_up[:,:,0], dtype = float)
for i in range(nx):for j in range(nz):# 任何一个点在时间域上互相关: 下行波的t0 * 上行波的t_{end} + 下行波的t1 * 上行波的t_{end-1} + ... + 下行波的t_{end} * 上行波的t_{0}migrated_image[i, j] = np.sum(u_up[i, j, :] * u_down[i, j, ::-1])

可以发现, RTM图像的两个维度都是空间维度上的, 而非共源单炮集合 (common shot gathers) – 也就是单炮地震图像(山峰) – 是时间-空间维度意义上的. 这一点可以认为RTM在空间意义上更接近速度模型.
最终得到RTM图像如下所示 (包括速度模型对比):
在这里插入图片描述
因为我们的炮点位置在地表中点, 因此地下结构中的中间区域信息是比较准确的. 盐体中间部分竟可以只依靠正演得到如此准确的结果, 令人震惊, 上层的层次结构看似混乱, 但是也是有规可循的, 因为层边界似乎信号更集中一些.
我进一步将炮点的位置移动到75和225, 观测它们的RTM图像:
在这里插入图片描述
可以发现, 炮点正下方附近的信息预测都很不错, 轮廓细节与速度模型很接近, 但是距离炮点远点会因为衰减而预测不准确. 但是这些信息应当可以通过部分拼接来取长补短地进行进行融合! 于是我试着融合了一下:
在这里插入图片描述
左图是我直接把所有炮点图对应的RTM图像叠加取平均, 右图是我把炮点下方宽度为60的区域截取下来, 然后把所有RTM截取部分拼在一起 (重叠部分取平均).
效果非常amazing啊, 没想到仅仅是分析正演的波形就可以得到这样的效果.
但是声波信息相对简单, 而且地震波属于弹性波, 因此我们需要更真实的正演信息来进一步验证RTM的可行性.

2.3 Python复现包含复杂的边界吸收条件的弹性波正演.

于是我将一个拥有更高级的边界吸收能力的六阶中心差分的弹性波正演代码移植到了Python上.
源算法文件为.mat文件, 是师兄提供的.
这个算法相对来说就复杂很多了, 不方便介绍, 下面直接提供复现的代码. (一些注释来自@苗妮)

from __future__ import division
import numpy as np
from bruges.filters.wavelets import ricker
import matplotlib.pyplot as plt
from scipy.signal import convolve
import skimage.filters
import matplotlib
import scipy.io
import math
import cv2
matplotlib.use('TkAgg')vmodel = scipy.io.loadmat("./vmodel1182.mat")["vmodel"]
nz, nx = vmodel.shape
h = 10                                                                           # PML宽度
Nz, Nx = nz + 2 * h, nx + 2 * h                                                  # 套上PML层后的长度
sz, sx = 0, 150                                                                  # 震源位置
dx = 10
dz = 10
nt = 2000                                                                        # 时间采样个数
dt = 0.001
wave_source_ampl = 1                                                             # 震源振幅
nodr = 3                                                                         # 空间差分阶数的一半
f0 = 25                                                                          # 震源频率
sampling_time = dt * nt                                                          # 采样时长 = 采样间隔时间 x 采样个数##################################
# 计算矩阵C :进行差分计算时的系数 #
##################################
B = np.zeros([nodr, 1])
B[0][0] = 1
A = np.zeros([nodr, nodr])
for i in range(nodr):A[i, :] = np.power(np.arange(1, 2 * nodr, 2), 2 * i + 1)
C = np.dot(np.linalg.inv(A), B).reshape(-1)################################### 计算密度矩阵 #
##################################
rho = 1000 * np.ones([Nz, Nx])                                                   # 密度矩阵包含了套上的PML边
vmodel_pad = np.pad(vmodel, [h, h], 'edge')################################### 生成震源波 #
##################################
# 标准雷克子波
source_wav = ricker(duration = 0.08, dt=dt, f=f0)[0]
source_wave_duration = len(source_wav)
source_wav = np.concatenate((source_wav, 0 * np.ones(nt - len(source_wav))))################################### 预览 #
##################################
plt.plot(source_wav)
plt.show()
plt.imshow(vmodel)
plt.show()
plt.imshow(skimage.filters.gaussian(vmodel, smooth_value))
plt.show()################################### PML层的吸收系数计算 #
##################################
# The PML formula refers to the equations (2) and (3) of Marcinkovich and Olsen, 2003.
v_max = np.max(vmodel_pad)
dp_z = np.zeros([Nz, Nx])
dp_x = np.zeros([Nz, Nx])# 设置上下两层
dp0_z = 3 * v_max / dz * (8 / 15 - 3 * h / 100 + 1 * h**2 / 1500)
# 计算边缘的吸收因子
dp_z[:h, :] = np.dot(dp0_z * np.power(np.arange(h, 0, -1) / h, 2).reshape(-1, 1), np.ones([1, Nx]))
dp_z[(Nz - h):, :] = dp_z[h-1::-1, :]# 设置左右两层
dp0_x = 3 * v_max / dx * (8 / 15 - 3 * h / 100 + 1 * h**2 / 1500)
# 计算边缘的吸收因子
dp_x[:, :h] = np.dot(np.ones([Nz, 1]), dp0_x * np.power(np.arange(h, 0, -1) / h, 2).reshape(1, -1))
dp_x[:, (Nx - h):] = dp_x[:, (h-1)::-1]#################################### 依据广义胡克定律求的弹性系数 #
###################################
rho1 = rho.copy()
rho2 = rho.copy()
# Coeffi1 和 Coeffi2 是沿 x 轴和 z 轴方向的PML吸收因子的系数.
# 它们根据PML吸收因子(dp_x 和 dp_z)和时间步长 dt 计算得出. 这些系数在波场更新方程中用于考虑PML的吸收效果.
Coeffi1 = (2 - dt * dp_x) / (2 + dt * dp_x)
Coeffi2 = (2 - dt * dp_z) / (2 + dt * dp_z)
# Coeffi3 和 Coeffi4 是与密度(rho)和空间步长(dx 和 dz)相关的系数, 用于考虑波场更新方程中的空间导数项.
Coeffi3 = 1 / rho1 / dx * (2 * dt / (2 + dt * dp_x))
Coeffi4 = 1 / rho2 / dz * (2 * dt / (2 + dt * dp_z))
# Coeffi5 和 Coeffi6 是与密度(rho)和速度(vp)的平方以及空间步长(dx 和 dz)相关的系数, 用于考虑波场更新方程中的速度和应力项.
Coeffi5 = rho * np.power(vmodel_pad, 2) / dx * (2 * dt / (2 + dt * dp_x))
Coeffi6 = rho * np.power(vmodel_pad, 2) / dz * (2 * dt / (2 + dt * dp_z))###################################### 迭代前初始化 #
###################################### 设置外部空间: 所有波场值均为空, 防止越界
NZ = Nz + 2 * nodr
NX = Nx + 2 * nodr
# 套上PML层之后的有效索引区域
Znodes = np.arange(nodr, NZ - nodr, 1)
Xnodes = np.arange(nodr, NX - nodr, 1)
# 原图有效索引区域
znodes = np.arange(nodr + h, nodr + h + nz)
xnodes = np.arange(nodr + h, nodr + h + nx)
# 套上PML层和外部空间的震源位置
sz_pad = nodr + h + sz
sx_pad = nodr + h + sx
# 初始化应力矩阵和速度有向分量有关矩阵
Ut = np.zeros([NZ, NX])
Uz = np.zeros([NZ, NX])
Ux = np.zeros([NZ, NX])
Vz = np.zeros([NZ, NX])
Vx = np.zeros([NZ, NX])
U = -1 * np.ones([nz, nx, nt])
Psum = -1 * np.ones([Nz,Nx])def index_2dim(matrix, height_index_array, width_index_array):return matrix[height_index_array[:, np.newaxis], width_index_array]print("开始进行时间迭代...")
for cur_time in range(nt):Ux[sz_pad, sx_pad] = Ux[sz_pad, sx_pad] + wave_source_ampl * source_wav[cur_time] / 2# 震源点的外力补充情况Uz[sz_pad, sx_pad] = Uz[sz_pad, sx_pad] + wave_source_ampl * source_wav[cur_time] / 2Ut = Ux + Uz																 # Ut为Ux和Uz两个分力矩阵的结合U[:, :, cur_time] = index_2dim(Ut, znodes, xnodes)if (cur_time + 1) % 200 == 0:print("{}/{}".format(cur_time + 1, nt))Psum[:,:] = 0                                                                # 中间变量矩阵for i in range(1, nodr + 1):Psum = Psum + C[i-1] * (index_2dim(Ut, Znodes, Xnodes + i) - index_2dim(Ut, Znodes, Xnodes + 1 - i))Vx[nodr:NZ - nodr, nodr:NX - nodr] = Coeffi1 * index_2dim(Vx, Znodes, Xnodes) - Coeffi3 * PsumPsum[:, :] = 0for i in range(1, nodr + 1):Psum = Psum + C[i-1] * (index_2dim(Ut, Znodes + i, Xnodes) - index_2dim(Ut, Znodes + 1 - i, Xnodes))Vz[nodr:NZ - nodr, nodr:NX - nodr] = Coeffi2 * index_2dim(Vz, Znodes, Xnodes) - Coeffi4 * PsumPsum[:, :] = 0for i in range(1, nodr + 1):Psum = Psum + C[i-1] * (index_2dim(Vx, Znodes, Xnodes - 1 + i) - index_2dim(Vx, Znodes, Xnodes - i))Ux[nodr:NZ - nodr, nodr:NX - nodr] = Coeffi1 * index_2dim(Ux, Znodes, Xnodes) - Coeffi5 * PsumPsum[:, :] = 0for i in range(1, nodr + 1):Psum = Psum + C[i-1] * (index_2dim(Vz, Znodes - 1 + i, Xnodes) - index_2dim(Vz, Znodes - i, Xnodes))Uz[nodr:NZ - nodr, nodr:NX - nodr] = Coeffi2 * index_2dim(Uz, Znodes, Xnodes) - Coeffi6 * Psum
print("迭代完成!")common_shot_gather = U[1, :, :].T# 走时信息静默
muted_gather = common_shot_gather.copy()x_array = np.arange(0, nx*dx, dx
v0 = vmodel[0,:] 
traveltimes = abs(np.cumsum(dx/v0) - np.cumsum(dx/v0)[sx])
for traceno in range(len(x_array)):muted_gather[0:int(traveltimes[traceno]/dt + source_wave_duration), traceno] = 0# 对时间域进行下采样
# common_shot_gather = cv2.resize(common_shot_gather, dsize=(301, 400), interpolation=cv2.INTER_CUBIC) # 展示共炮集合图
fgr2, axs2 = plt.subplots(1,2, figsize = (16,6))
im1=axs2[0].imshow(common_shot_gather, cmap = plt.cm.seismic, vmin=-0.4, vmax=0.44, aspect='auto')
im2=axs2[1].imshow(muted_gather, cmap = plt.cm.seismic, aspect='auto')
axs2[0].set_title('Seismic wave structure (muting the forward wave)')
axs2[1].set_title('Seismic wave structure')
fgr2.colorbar(im1, ax = axs2[0])
fgr2.colorbar(im2, ax = axs2[1])
plt.show()

这个代码中, 我将走时静默的代码一并加入了进来, 左图为共源炮点图 (即单炮地震图像), 右图是它走时信息静默后的图像.

在这里插入图片描述

之前的单炮图像和走时信息静默图像, 可以发现, 采用更高级的边界吸收条件后, 波变"干净"了很多. 这是因为一些边界信息吸收不完全的地方和低阶差分对波拟合的误差都被很好地消除了.

2.4 复杂的弹性波正演可否实现RTM?

不知道为什么, 我尝试基于弹性波的正演过程还原RTM中上行波后进行互相关…但是效果并不这么好. 于是我有点怀疑代码的问题, 或者是我哪里理解有问题, 于是思考了几天.
最后我想出了一个思路, 不妨基于声波的二阶中心差分来构造上行波, 虽然上行波正演的过程基于的震源是弹性波的结构, 但是会不会有意想不到的效果呢? 当然, 推动我这个思路的缘由在于: 实际地质探测的过程中, 或者历史上的已有地震探测记录中, 可能根本没有上行波的地震资料, 因此没必要纠结上行波在正演的过程中是否足够真实, 毕竟RTM是一种手段, 并不需要多"像"真实情况, 因此这个波形可以脱离弹性波的范畴 (I guess…and I hope so).
结果如下:
在这里插入图片描述
效果非常amazing啊! 能得到完美的轮廓, 而且相比于声波波形的单炮RTM (炮点位置在中点) 那张图, 这里明显更"干净", 扰动少了很多. 于是, 照猫画虎, 我们试着模仿2.2中的RTM图像的叠加方法, 来看下弹性波RTM的叠加效果图:
在这里插入图片描述
效果还挺不错的. 相比于声波, 我们通过更加准确的弹性波, 更高阶的差分手段, 以及更高级的边界吸收条件为基础, 生成了更加 “干净” 的叠前RTM图像. 从结果上来说, RTM可以通过弹性波来实现, 而且效果非常好! (当然很大程度得以中心差分阶数和边界吸收条件).

3 存在的主要问题

简单总结这一周干的事情的意义:
目前来说, 我们可以在已知

  1. 正演的所有时间序列的波场信息 (弹性波) 地震观测数据 (CSG: Common-shot gathers)
  2. 初始速度模型

的情况下, 可以生成上面这种和速度模型高度接近的波形模拟图像.

  • 问题1: (应该好克服)
    我们需要提供一个初始速度模型, 这个倒是符合传统的FWI的观点, 但是有悖于DL-FWI的观念.
    但我们的思维也不见得如此死板, 如果说能以提供一个初始速度模型为代价, 来从而得到如此精确的速度模型波形模拟图像, 试想倘若用它来指导网络训练, 那么最终获得的速度模型将会多么精确?
    这个是可以谈的.

  • 问题2: (有点不确定)
    但是我最担忧的还是 " 正演的所有时间序列的波场信息 " 这个资料是否是工区现场能够提供的信息, 因为这个必须作为和上行波进行互相关才能得到RTM图像. 如果可以提供这个信息或者可以通过已知的单炮的山峰图逆向生成得到, 那么这个问题就可以推广到实际的预测. 否则RTM可能还是只是空中楼阁 .
    问题2已经得到了解决, 通过我的猜测并最终测试, 基本可以肯定并不需要现场的全时间域的波场快照就可以做出RTM图像. RTM需要上行波和下行波的波场快照进行互相关是没错的, 但这里的下行波波场快照实际没必要用原速度模型的正演资料进行, 准确说, 上下行波的正演都只基于初始的速度模型就可以做了.

在这里插入图片描述

最后再提一嘴RTM在空间意义上的优势:
叠前多炮数据本身的纵向维度是时间而非 速度模型 采用的深度, 这一点加大了端到端的神经网络拟合的难度 (要拟合不同时空域的关系), 而且也会被质疑解释性.
但是RTM图像本身和速度模型的维度完全一致, 并且边缘结构上存在大量相似信息. 若以之作为网络的输入的辅助资料, 将极大提高了网络拟合的效率和准确性, 降低了拟合边缘的难度 (边缘结构样貌的拟合是目前DL-FWI网络拟合的最大难点). 可解释性非常完美.

4 下一步工作

继续完成上周提到的论文阅读, 看看有关论文是如何利用RTM来完成DL-FWI的任务.
这一点的继续了解, 希望能解决问题2中的担忧.

这篇关于2023 7.10~7.16 周报 (RTM研究与正演的Python复现) (8.3更新)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/384313

相关文章

python: 多模块(.py)中全局变量的导入

文章目录 global关键字可变类型和不可变类型数据的内存地址单模块(单个py文件)的全局变量示例总结 多模块(多个py文件)的全局变量from x import x导入全局变量示例 import x导入全局变量示例 总结 global关键字 global 的作用范围是模块(.py)级别: 当你在一个模块(文件)中使用 global 声明变量时,这个变量只在该模块的全局命名空

poj3468(线段树成段更新模板题)

题意:包括两个操作:1、将[a.b]上的数字加上v;2、查询区间[a,b]上的和 下面的介绍是下解题思路: 首先介绍  lazy-tag思想:用一个变量记录每一个线段树节点的变化值,当这部分线段的一致性被破坏我们就将这个变化值传递给子区间,大大增加了线段树的效率。 比如现在需要对[a,b]区间值进行加c操作,那么就从根节点[1,n]开始调用update函数进行操作,如果刚好执行到一个子节点,

hdu1394(线段树点更新的应用)

题意:求一个序列经过一定的操作得到的序列的最小逆序数 这题会用到逆序数的一个性质,在0到n-1这些数字组成的乱序排列,将第一个数字A移到最后一位,得到的逆序数为res-a+(n-a-1) 知道上面的知识点后,可以用暴力来解 代码如下: #include<iostream>#include<algorithm>#include<cstring>#include<stack>#in

hdu1689(线段树成段更新)

两种操作:1、set区间[a,b]上数字为v;2、查询[ 1 , n ]上的sum 代码如下: #include<iostream>#include<algorithm>#include<cstring>#include<stack>#include<queue>#include<set>#include<map>#include<stdio.h>#include<stdl

【Python编程】Linux创建虚拟环境并配置与notebook相连接

1.创建 使用 venv 创建虚拟环境。例如,在当前目录下创建一个名为 myenv 的虚拟环境: python3 -m venv myenv 2.激活 激活虚拟环境使其成为当前终端会话的活动环境。运行: source myenv/bin/activate 3.与notebook连接 在虚拟环境中,使用 pip 安装 Jupyter 和 ipykernel: pip instal

【机器学习】高斯过程的基本概念和应用领域以及在python中的实例

引言 高斯过程(Gaussian Process,简称GP)是一种概率模型,用于描述一组随机变量的联合概率分布,其中任何一个有限维度的子集都具有高斯分布 文章目录 引言一、高斯过程1.1 基本定义1.1.1 随机过程1.1.2 高斯分布 1.2 高斯过程的特性1.2.1 联合高斯性1.2.2 均值函数1.2.3 协方差函数(或核函数) 1.3 核函数1.4 高斯过程回归(Gauss

hdu 1754 I Hate It(线段树,单点更新,区间最值)

题意是求一个线段中的最大数。 线段树的模板题,试用了一下交大的模板。效率有点略低。 代码: #include <stdio.h>#include <string.h>#define TREE_SIZE (1 << (20))//const int TREE_SIZE = 200000 + 10;int max(int a, int b){return a > b ? a :

AI行业应用(不定期更新)

ChatPDF 可以让你上传一个 PDF 文件,然后针对这个 PDF 进行小结和提问。你可以把各种各样你要研究的分析报告交给它,快速获取到想要知道的信息。https://www.chatpdf.com/

【学习笔记】 陈强-机器学习-Python-Ch15 人工神经网络(1)sklearn

系列文章目录 监督学习:参数方法 【学习笔记】 陈强-机器学习-Python-Ch4 线性回归 【学习笔记】 陈强-机器学习-Python-Ch5 逻辑回归 【课后题练习】 陈强-机器学习-Python-Ch5 逻辑回归(SAheart.csv) 【学习笔记】 陈强-机器学习-Python-Ch6 多项逻辑回归 【学习笔记 及 课后题练习】 陈强-机器学习-Python-Ch7 判别分析 【学

GIS图形库更新2024.8.4-9.9

更多精彩内容请访问 dt.sim3d.cn ,关注公众号【sky的数孪技术】,技术交流、源码下载请添加微信:digital_twin123 Cesium 本期发布了1.121 版本。重大新闻,Cesium被Bentley收购。 ✨ 功能和改进 默认启用 MSAA,采样 4 次。若要关闭 MSAA,则可以设置scene.msaaSamples = 1。但是通过比较,发现并没有多大改善。