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

Python地震波逆问题解构算法复杂信号分析 | 亚图跨际亚图跨际 on Notion

🏈指点迷津 | Brief

🎯要点

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

🍪语言内容分比

🍇Python波场逆问题

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

为了处理逆问题，我们必须首先了解正问题。在最基本的形式中，波动方程由下式给出

$$L[c] v(t, x)=q(t, x)\qquad(1)$$

且

$$L[c]=\left[\frac{1}{c(x)^2} \frac{\partial^2}{\partial t^2}-\nabla^2\right]$$

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

我们研究逆源问题的一个变体，其中 $q(t, x)=\delta(t) u(x)$ 和 $u$ 在$\Omega \subset R ^n$上得到紧凑支持。常量$c$ 的前向运算符由下式给出

$$K u(t, x)=\int_{\Omega} u\left(x^{\prime}\right) g\left(t, x-x^{\prime}\right) d x^{\prime}$$

在半径为 $\rho$ 的球体上进行测量。解决逆问题的一种流行技术是反向传播，它基于将前向算子的伴随应用于数据。在这种情况下，伴随运算符由下式给出：

$$K^* 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^* 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^*K$。在时域傅立叶域中，对于 $c=1$，运算符变为

$$\widehat{f}(\omega, x)=\int_{\Omega} u\left(x^{\prime}\right) \frac{\exp \left(\imath \omega\left|x-x^{\prime}\right|\right)}{\left|x-x^{\prime}\right|} d x^{\prime}$$

和

$$u(x)=\iint_{\Delta} \widehat{f}\left(\omega, x^{\prime}\right) \frac{\exp \left(-\imath \omega\left|x^{\prime}-x\right|\right)}{\left|x^{\prime}-x\right|} d x^{\prime} d \omega$$

于是，

$$K^* K u(x)=\iint_{\Delta} \int_{\Omega} f\left(x^{\prime}\right) \frac{\exp \left(\imath \omega\left|x^{\prime \prime}-x^{\prime}\right|\right)}{\left|x^{\prime \prime}-x^{\prime}\right|} \frac{\exp \left(-\imath \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^* K u(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 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()