🍍Python地震波逆问题解构算法复杂信号分析
🎯要点
🎯时域、时频域以及时间和频率相关联偏振特性分析三种算法 | 🎯时域波参数估计算法 | 🎯机器学习模型波形指纹分析算法 | 🎯色散曲线和频率相关波分析算法 | 🎯动态倾斜校正算法 | 🎯声学地震波像分析
🍪语言内容分比
🍇Python波场逆问题
在许多应用中,声波或电磁波被用来观察事物:水下声学、雷达、声纳、医学超声波、探地雷达、地震成像、全球地震学。在其中一些应用中,测量是被动的,我们记录调查对象发出的波,并试图推断其位置、大小等。在主动测量中,会产生波并记录对象的响应。
为了处理逆问题,我们必须首先了解正问题。在最基本的形式中,波动方程由下式给出
且
其中 表示波场, 是传播速度, 是源术语。我们假设源在空间和时间上都有紧支持并且是平方可积的。
为了定义(1)的唯一解,我们需要提供边界和初始条件。在上面讨论的应用中,通常考虑无界域 和柯西初始条件 。
在散射实验中,通常根据入射波场和散射波场重写波动方程 v=v_i+v_s 和散射势 :
在弱散射假设下,我们可以忽略u和v_s的相互作用,并将(2)中的替换为。
测量结果通常被视为限制为 的(分散)波场,其中 可以是流形或一组点。然后数据可以用表示。在主动实验中,通常的做法是考虑对源项集合的测量。源可以是点源,在这种情况下 且 。或者,事件场 可以由 给出和参数化。然后我们可以用 表示数据,其中和 。
基于这个基本设置,我们将讨论三个可能的逆问题:
逆源问题:从已知(且恒定) 的测量中恢复源。
逆散射:假设 已知,从多个(已知)源的散射场测量中恢复散射势 。
波形断层扫描:从多个(已知)源的总波场测量中恢复。
下面您将找到一些实践中出现的逆问题的典型示例。
地震定位:由源项 描述的地震由多个地震仪在位置 处记录。目标是恢复 以确定地震位置。
被动声纳:使用数组 记录从不明目标发出的声波,其中 表示数组的方向, 表示其宽度。目标是恢复源项 以确定源的起源和性质。
雷达成像:入射平面波,由方向 参数化,被发送到介质中,并由阵列记录其反射响应。目标是恢复散射势。
全波形反演:在勘探地震学中,目标是从表面总波场的测量中恢复。
超声断层扫描:目标是康复 𝑐 来自物体周围的源和接收器的总波场。
我们研究逆源问题的一个变体,其中 和 在上得到紧凑支持。常量 的前向运算符由下式给出
在半径为 的球体上进行测量。解决逆问题的一种流行技术是反向传播,它基于将前向算子的伴随应用于数据。在这种情况下,伴随运算符由下式给出:
我们看到 可以通过求解下式得到
使用时间反转格林函数并在 处求值,即 。
为了了解其原理,我们研究普通运算符 。在时域傅立叶域中,对于 ,运算符变为
和
于是,
对于 我们可以将其近似为 同样对于。这称为远场近似。引入单位球面,我们发现
此式:
有时称为点扩散函数。
这种过滤也可以通过乘以 在傅立叶域中实现。一个例子如下所示。
import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
from scipy.ndimage import laplace
def solve(q,c,dt,dx,T=1.0,L=1.0,n=1):
gamma = dt*c/dx
nt = int(T/dt + 1)
nx = int(2*L/dx + 2)
u = np.zeros((nx**n,nt))
for k in range(1,nt-1):
u[:,k+1] = A@u[:,k] + L@u[:,k] + B@u[:,k-1] + (dx**2)*C@q[:,k]
return u
def getMatrices(gamma,nx,n):
l = (gamma**2)*np.ones((3,nx))
l[1,:] = -2*(gamma**2)
l[1,0] = -gamma
l[2,0] = gamma
l[0,nx-2] = gamma
l[1,nx-1] = -gamma
if n == 1:
a = 2*np.ones(nx)
a[0] = 1
a[nx-1] = 1
b = -np.ones(nx)
b[0] = 0
b[nx-1] = 0
c = (gamma)**2*np.ones(nx)
c[0] = 0
c[nx-1] = 0
L = sp.diags(l,[-1, 0, 1],shape=(nx,nx))
else:
a = 2*np.ones((nx,nx))
a[0,:] = 1
a[nx-1,:] = 1
a[:,0] = 1
a[:,nx-1] = 1
a.resize(nx**2)
b = -np.ones((nx,nx))
b[0,:] = 0
b[nx-1,:] = 0
b[:,0] = 0
b[:,nx-1] = 0
b.resize(nx**2)
c = (gamma)**2*np.ones((nx,nx))
c[0,:] = 0
c[nx-1,:] = 0
c[:,0] = 0
c[:,nx-1] = 0
c.resize(nx**2)
L = sp.kron(sp.diags(l,[-1, 0, 1],shape=(nx,nx)),sp.eye(nx)) + sp.kron(sp.eye(nx),sp.diags(l,[-1, 0, 1],shape=(nx,nx)))
return A,B,C,L
L = 1.0
T = 1.0
dx = 1e-2
dt = .5e-2
nt = int(T/dt + 1)
nx = int(2*L/dx + 2)
c = 1
u = np.zeros((nx,nx))
u[nx//2 - 10:nx//2+10,nx//2 - 10:nx//2+10] = 1
q = np.zeros((nx*nx,nt))
q[:,1] = u.flatten()
w_forward = solve(q,c,dt,dx,T=T,L=L,n=2)
theta = np.linspace(0,2*np.pi,20,endpoint=False)
xs = 0.8*np.cos(theta)
ys = 0.8*np.sin(theta)
I = np.ravel_multi_index(np.array([[xs/dx + nx//2],[ys//dx + nx//2]],dtype=np.int), (nx,nx))
d = w_forward[I,:]
r = np.zeros((nx*nx,nt))
r[I,:] = d
r = np.flip(r,axis=1)
w_adjoint = solve(r,c,dt,dx,T=T,L=L,n=2)
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5), sharey=True)
plt.gray()
ax[0,0].plot(xs,ys,'r*')
ax[0,0].imshow(w_forward[:,2].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[0,0].set_title('t = 0')
ax[0,1].plot(xs,ys,'r*')
ax[0,1].imshow(w_forward[:,101].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[0,1].set_title('t = 0.5')
ax[0,2].plot(xs,ys,'r*')
ax[0,2].imshow(w_forward[:,200].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[0,2].set_title('t = 1')
ax[1,0].plot(xs,ys,'r*')
ax[1,0].imshow(w_adjoint[:,200].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[1,1].plot(xs,ys,'r*')
ax[1,1].imshow(w_adjoint[:,101].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[1,2].plot(xs,ys,'r*')
ax[1,2].imshow(w_adjoint[:,2].reshape((nx,nx)), extent=(-L,L,-L,L))
plt.show()
Last updated
Was this helpful?