Python地震波逆问题解构算法复杂信号分析

🎯要点

🎯时域、时频域以及时间和频率相关联偏振特性分析三种算法 | 🎯时域波参数估计算法 | 🎯机器学习模型波形指纹分析算法 | 🎯色散曲线和频率相关波分析算法 | 🎯动态倾斜校正算法 | 🎯声学地震波像分析

📜Python,MATLAB和C/C++物理信号滤波合成检测模拟用例

📜Python和MATLAB数字信号波形和模型模拟

📜Python量化噪声卷积信号和傅里叶时频分析

📜C++数字化声音信号处理和数控振荡合成

📜Python嵌入式动态用户调制解调响应式射频信号

📜Python嵌入式片上系统逻辑电路物理信号处理

📜Arduino生物波反馈和环境检测外套

📜Python线性代数傅里叶分析和动态系统模拟分析之一

📜Python线性代数数字图像和小波分析之二

📜Python数值和符号算法计算及3D视图物理数学波形方程

📜MATLAB_ESP32有限脉冲响应FIR无限脉冲响应IIR滤波器

🍪语言内容分比

🍇Python波场逆问题

在许多应用中，声波或电磁波被用来观察事物：水下声学、雷达、声纳、医学超声波、探地雷达、地震成像、全球地震学。在其中一些应用中，测量是被动的，我们记录调查对象发出的波，并试图推断其位置、大小等。在主动测量中，会产生波并记录对象的响应。

为了处理逆问题，我们必须首先了解正问题。在最基本的形式中，波动方程由下式给出
$$L[c]v(t,x)=q(t,x)(1)$$
且
$$L[c]=\left[\frac{1}{c(x{)}^2}\frac{{\partial }^2}{\partial t^2}-{\nabla }^2\right]$$
其中 $v:R\times R^n\to R$ 表示波场，$c:R^n\to R$ 是传播速度，$q:R\times R^n\to R$ 是源术语。我们假设源在空间和时间上都有紧支持并且是平方可积的。

为了定义(1)的唯一解，我们需要提供边界和初始条件。在上面讨论的应用中，通常考虑无界域 $\left(x\in R^n\right)$ 和柯西初始条件 $v(0,x)=0,\frac{\partial v}{\partial t}(0,x)=0$。

在散射实验中，通常根据入射波场和散射波场重写波动方程 $v=v_i+v_s$ 和散射势 $u(x)=c(x{)}^{-2}-c_0^{-2}$ :
$$L\left[c_0\right]v_i(t,x)=q(t,x),L\left[c_0\right]v_s(t,x)=-u(x)\frac{{\partial }^{2}v}{\partial t^2}(t,x)(2)$$
在弱散射假设下，我们可以忽略$u$和$v_s$的相互作用，并将(2)中的$v$替换为$v_i$。

测量结果通常被视为限制为 $[0,T]\times \Delta$ 的（分散）波场，其中 $\Delta \subset R^n$ 可以是流形或一组点。然后数据可以用$f(t,r)$表示。在主动实验中，通常的做法是考虑对源项集合的测量。源可以是点源，在这种情况下 $q(t,x)=w(t)\delta (x-s)$ 且 $s\in \Sigma$。或者，事件场 $v_i$ 可以由 $s\in \Sigma$ 给出和参数化。然后我们可以用 $f(t,s,r)$ 表示数据，其中 $s\in \Sigma 、r\in \Delta$ 和 $t\in [0,T]$。

基于这个基本设置，我们将讨论三个可能的逆问题：

• 逆源问题：从已知（且恒定）$c(x)\equiv c_0$ 的测量中恢复源$q$。
• 逆散射：假设 $c_0$ 已知，从多个（已知）源的散射场测量中恢复散射势 $u(x)$。
• 波形断层扫描：从多个（已知）源的总波场测量中恢复 $c(x)$。

下面您将找到一些实践中出现的逆问题的典型示例。

• 地震定位：由源项 $w(t)q(x)$ 描述的地震由多个地震仪在位置 Δ = { r k } k = 1 n 处记录 \Delta=\left\{r_k\right\}_{k=1}^n 处记录 Δ={rk​}k=1n​处记录。目标是恢复 $q$ 以确定地震位置。
• 被动声纳：使用数组 Δ = { x 0 + r p ∣ r ∈ [ − h , h ] } \Delta=\left\{x_0+r p \mid r \in[-h, h]\right\} Δ={x0​+rp∣r∈[−h,h]} 记录从不明目标发出的声波，其中 $p \in S ^2 $ 表示数组的方向，$h$ 表示其宽度。目标是恢复源项 $w(t)q(x)$ 以确定源的起源和性质。
• 雷达成像：入射平面波，由方向 $s\in \Sigma \subset S^2$​ 参数化，被发送到介质中，并由阵列记录其反射响应。目标是恢复散射势。
• 全波形反演：在勘探地震学中，目标是从表面总波场的测量中恢复$c$： Σ = Δ = { x ∣ n ⋅ x = 0 } \Sigma=\Delta=\{x \mid n \cdot x=0\} Σ=Δ={x∣n⋅x=0}。
• 超声断层扫描：目标是康复 𝑐 来自物体周围的源和接收器的总波场。

我们研究逆源问题的一个变体，其中 $q(t,x)=\delta (t)u(x)$ 和 $u$ 在 $\Omega \subset R^n$ 上得到紧凑支持。常量 $c$ 的前向运算符由下式给出
$$Ku(t,x)={\int }_{\Omega }u\left(x^{\prime }\right)g\left(t,x-x^{\prime }\right)d{x}^{\prime }$$
在半径为 $\rho$ 的球体上进行测量。解决逆问题的一种流行技术是反向传播，它基于将前向算子的伴随应用于数据。在这种情况下，伴随运算符由下式给出：
$$K^{\ast }f(x)={\int }_{\Delta }{\int }_0^{T}g\left(t^{\prime },x^{\prime }-x\right)f\left(t^{\prime },x^{\prime }\right)d{x}^{\prime }d{t}^{\prime }$$
我们看到 $p=K^{\ast }f$ 可以通过求解下式得到
$$L[c]w(t,x)={\int }_{\Delta }f\left(t,x^{\prime }\right)\delta \left(x-x^{\prime }\right)d{x}^{\prime }$$
使用时间反转格林函数并在 $t=0$ 处求值，即 $p(x)=w(0,x)$。

为了了解其原理，我们研究普通运算符 $K^{\ast }K$。在时域傅立叶域中，对于 $c=1$，运算符变为
$$\hat{f}(\omega ,x)={\int }_{\Omega }u\left(x^{\prime }\right)\frac{\exp \left(ı\omega \left|x-x^{\prime }\right|\right)}{\left|x-x^{\prime }\right|}d{x}^{\prime }$$
和
$$u(x)={\iint }_{\Delta }\hat{f}\left(\omega ,x^{\prime }\right)\frac{\exp \left(-ı\omega \left|x^{\prime }-x\right|\right)}{\left|x^{\prime }-x\right|}d{x}^{\prime }d\omega$$
于是，
$$K^{\ast }Ku(x)={\iint }_{\Delta }{\int }_{\Omega }f\left(x^{\prime }\right)\frac{\exp \left(ı\omega \left|x^{\prime \prime }-x^{\prime }\right|\right)}{\left|x^{\prime \prime }-x^{\prime }\right|}\frac{\exp \left(-ı\omega \left|x^{\prime \prime }-x\right|\right)}{\left|x^{\prime \prime }-x\right|}d{x}^{\prime }d{x}^{\prime \prime }d\omega$$
对于 $\left|x^{\prime \prime }\right|\gg |x|$ 我们可以将其近似为 $\left|x^{\prime \prime }-x\right|\approx \left|x^{\prime \prime }\right|+x\cdot x^{\prime \prime }/\left|x^{\prime \prime }\right|$ 同样对于 $\left|x^{\prime \prime }-x^{\prime }\right|$。这称为远场近似。引入${\xi }^{\prime \prime }=x^{\prime \prime }/\left|x^{\prime \prime }\right|=x^{\prime \prime }/\rho$单位球面，我们发现
$$K^{\ast }Ku(x)=\rho {\int }_{\Omega }u\left(x^{\prime }\right)\iint \exp \left(w{\xi }^{\prime \prime }\cdot \left(x^{\prime }-x\right)\right)d{\xi }^{\prime \prime }d\omega d{x}^{\prime }$$
此式：
$$k\left(x-x^{\prime }\right)=\iint \exp \left(u\omega {\xi }^{\prime \prime }\cdot \left(x^{\prime }-x\right)\right)d{\xi }^{\prime \prime }d\omega$$

有时称为点扩散函数。

这种过滤也可以通过乘以 $|\xi {|}^2$ 在傅立叶域中实现。一个例子如下所示。

import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
from scipy.ndimage import laplacedef solve(q,c,dt,dx,T=1.0,L=1.0,n=1):gamma = dt*c/dxnt = 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 udef getMatrices(gamma,nx,n):l = (gamma**2)*np.ones((3,nx))l[1,:] = -2*(gamma**2)l[1,0] = -gammal[2,0] = gammal[0,nx-2] = gammal[1,nx-1] = -gammaif n == 1:a = 2*np.ones(nx)a[0] = 1a[nx-1] = 1b = -np.ones(nx)b[0] = 0b[nx-1] = 0c = (gamma)**2*np.ones(nx)c[0] = 0c[nx-1] = 0L = sp.diags(l,[-1, 0, 1],shape=(nx,nx))else:a = 2*np.ones((nx,nx))a[0,:] = 1a[nx-1,:] = 1a[:,0] = 1a[:,nx-1] = 1a.resize(nx**2)b = -np.ones((nx,nx))b[0,:] = 0b[nx-1,:] = 0b[:,0] = 0b[:,nx-1] = 0b.resize(nx**2)c = (gamma)**2*np.ones((nx,nx))c[0,:] = 0c[nx-1,:] = 0c[:,0] = 0c[:,nx-1] = 0c.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 = 1u = 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()

👉参阅&更新：计算思维 | 亚图跨际