波动方程的有限差分解（地震数据正演）

原创Xenny约 3966 字大约 17 分钟

波动方程的有限差分解（地震数据正演）

有限差分法

有限差分法是一种求解微分方程的数值方法，对于微分方程来说微分项很难处理，有限差分的思路便是用近似的方法处理微分项，例如对于一个一阶微分方程
$\begin{aligned} u^{'}(x) + c(x)u(x) &= f(x)\\ u(a) &= d \end{aligned}\tag{1}$
由导数定义有
$u^{'}(x) = \lim\limits_{\Delta x \rightarrow 0}\frac{u(x+\Delta x) - u(x)}{\Delta x} \approx \frac{u(x + \Delta x) - u(x)}{\Delta x} \tag{2}$
这里的几何意义便是使用割线斜率近似替代切线斜率，显然 $\Delta x$ 越小则越精确，在后续求解中，我们将 $\Delta x$ 记为步长。
图1. 微分近似表示
由此在求解时我们要对求解区域进行离散化，即将区间按步长分为一个个离散点。设对于区间 $[a,b]$ 共分为 $n$ 个区间，步长为 $h$ ，坐标为 ${(x_0, u_0), (x_1, u_1),\dots, (x_n,u_n)}$ ，则我们可以得到各点处导数值得近似表示
$u_i^{'} = \frac{1}{h}(u_{i+1}- u_i)\tag{3}$
带入原方程得
$(u_{i+1} - u_i)/h + c(x_i)u_i = f(x_i)$
其中 $x_i = a+ih, i=0,1,2,\dots$ ，则我们有 $n$ 个方程 $(i=0,1,\dots,n-1)$ ，包含 $n$ 个未知数 $(u_1,u_2,\dots,u_n)$ 。此时我们对方程变形得到
$u_{i+1} = (1-hc(x_i))u_i + hf(x_i)\tag{4}$
则上述方程组可写为矩阵的形式
$\begin{bmatrix} 1&&&\\ C_1&1&&\\ &\ddots&\ddots&\\ &&c_{n-1}&1 \end{bmatrix}\begin{pmatrix} u_1\\ u_2\\ \vdots\\ n_n \end{pmatrix} = \begin{pmatrix} hf(x_0)\\ hf(x_1)\\ \vdots\\ hf(x_{n-1}) \end{pmatrix} - \begin{pmatrix} C_0u_0\\ 0\\ \vdots\\ 0 \end{pmatrix}\tag{5}$
其中 $C_i = hc(x_i)-1$ ，上述矩阵方程记为 $\mathbf{A}\boldsymbol{u} = \mathbf{F}$ ，则解为 $\boldsymbol{u} = \mathbf{A}^{-1}\mathbf{F}$ （这里显然 $r(\mathbf{A}) = n$ 即可逆）。

差分公式

对于有限差分来说，其原理大致和上述相同，只是不是简单的使用割线斜率来表示切线斜率，而是使用泰勒级数对函数进行展开，从而得到各阶的近似公式也即“差分公式”。
由泰勒级数我们有函数 $u(x)$ 在点 $x$ 处的展开为（前提是该点邻域 $n$ 阶可导）
$u(x+h) = u(x) + hu^{'}(x) + \frac{h^2}{2}u^{''}(x) + \cdots + \frac{h^k}{k!}u^{(n)}(x) + R_n\tag{6}$
根据需要来截取不同的近似阶数。
对于一阶微分，则有
$u(x+h) = u(x) + hu^{'}(x) + \frac{h^2}{2}u^{''}(x+\xi)\tag{7}$
其中 $\xi \in (0,h)$ ，则变形即可得到一阶微分的向前差分公式：
$u^{'}_F(x) = \frac{u(x+h)-u(x)}{h} - \frac{h}{2}u^{''}(x+\xi)\tag{8}$
同时将 $h$ 替换为 $-h$ ，我们可以得到一阶微分的向后差分公式：
$u^{'}_B(x) = \frac{u(x) - u(x-h)}{h} + \frac{h}{2}u^{''}(x+\xi)\tag{9}$
联立式(8)与式(9)，可以得到一阶微分的中心差分公式：
$u^{'}_C(x) = \frac{u(x+h) - u(x-h)}{2h} + O(h^2)\tag{10}$
同理也可以得到其他高阶微分的差分公式，其实和高数里面利用泰勒求高阶导类似。

偏微分方程

不过在后面的波动方程中，它是一个偏微分方程，所以我们还需要了解偏微分方程中的有限差分是如何近似的。
偏微分的概念其实也就是提高了维度，在一维中我们是得到导数近似的差分公式后对原方程离散化处理求解，这里也是一样的，假设我们要求一个二元偏微分方程，此时则将求解区域 $x\in [a,b], y\in [c,d]$ 进行网格化，在每个格点使用微分差分公式来近似表达偏导进行求解。
例如对于方程组
$\begin{cases} \Delta u = f(x, y)\\ u(0,y) = c(y),\ u(a,y) = d(y)\\ u(x,0) = g(x),\ u(x,b) = h(x) \end{cases}\tag{11}$
其中 $x \in [0,a], y \in[0,b]$ ，则分别将 $x,y$ 分为 $m,n$ 份，记为 $x_i(i=0,1,\dots,m), y_j(j=0,1,\dots,n)$ 。
由二阶微分中心差分公式有
$\frac{\partial^2 u}{\partial x^2} = \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h_{x}^2}, \frac{\partial^2 u}{\partial y^2} = \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h_{y}^2}\tag{12}$
其中 $h_x=a/m,h_y=b/n$ 分别为 $x$ 方向和 $y$ 方向的步长。令 $p=1/h_x^2,q=1/h_y^2,r=-2(p+1),f(x_i,y_j) = f_{i,j}$ 和式(12)带入式(11)得
$\begin{cases} qu_{i,j-1} + pu_{i-1,j} + ru_{i,j} + pu_{i+1,j} + qu_{i,j+1} = f_{i,j}\\ u_{0,j} = c_j,\ u_{m,j} = d_j\\ u_{i,0} = g_i,\ u_{i,n} = h_i \end{cases}\tag{13}$
显然我们可以得到 $m\times n$ 个这样的方程组，我们将上述方程组转换为矩阵形式为
$\begin{bmatrix} \mathbf{U}&\mathbf{V}&\\ \mathbf{V}&\mathbf{U}&\mathbf{V}&\\ &\ddots&\ddots&\ddots\\ &&\mathbf{V}&\mathbf{U}&\mathbf{V}\\ &&&\mathbf{V}&\mathbf{U} \end{bmatrix}\begin{pmatrix} \vec{u}_{i,1}\\ \vec{u}_{i,2}\\ \vdots\\ \vec{u}_{i,n-2}\\ \vec{u}_{i,n-1} \end{pmatrix} = \begin{pmatrix} \vec{f}_{i,1}\\ \vec{f}_{i,2}\\ \vdots\\ \vec{f}_{i,n-2}\\ \vec{f}_{i,n-1}\\ \end{pmatrix} - q\begin{pmatrix} \vec{u}_{i,0}\\ 0\\ \vdots\\ 0\\ 0 \end{pmatrix} - q\begin{pmatrix} 0\\0\\\vdots\\0\\\vec{u}_{i,m} \end{pmatrix}\\\ \\- p\begin{pmatrix} \vec{u}_{0,1}\\ \vec{u}_{0,2}\\ \vdots\\ \vec{u}_{0,n-2}\\ \vec{u}_{0,n-1} \end{pmatrix}- p\begin{pmatrix} \vec{u}_{m,1}\\ \vec{u}_{m,2}\\ \vdots\\ \vec{u}_{m,n-2}\\ \vec{u}_{m,n-1} \end{pmatrix}\tag{14}$
其中
$\mathbf{U} = \begin{bmatrix} r&p&\\ p&r&p\\ &\ddots&\ddots&\ddots&\\ &&p&r&p\\ &&&p&r \end{bmatrix}_{(m-1)\times (m-1)},\mathbf{V} = \begin{bmatrix} q\\ &q\\ &&\ddots\\ &&&q\\ &&&&q \end{bmatrix}_{(m-1)\times (m-1)}\\ \vec{u}_{i,j} = (u_{1,j}, u_{2,j}, \dots, u_{m-1,j})^T, \vec{f}_{i,j} = (f_{1,j},f_{2,j},\dots,f_{m-1,j})^T\\ \vec{u}_{i,0} = (u_{1,0}, u_{2,0}, \dots, u_{m-1,0})^T, \vec{u}_{i,n} = (u_{1,n},u_{2,n},\dots,u_{m-1,n})^T\\ \vec{u}_{0,j} = (u_{0,j}, 0,\dots, 0)^T,\vec{u}_{m,j} = (0,0,\dots,0,u_{m,j})^T$
这里用的是知乎文章^[1]的例子，里面没给推导过程。这里记录一下推导过程
首先是固定 $j$ ，将 $j$ 看成常量转为矩阵形式有
$\begin{bmatrix} r&p&\\ p&r&p\\ &\ddots&\ddots&\ddots&\\ &&p&r&p\\ &&&p&r \end{bmatrix}\begin{pmatrix} u_{1,j}\\ u_{2,j}\\ \vdots\\ u_{m-2,j}\\ u_{m-1,j} \end{pmatrix} = \begin{pmatrix} f_{1,j}\\ f_{2,j}\\ \vdots\\ f_{m-2,j}\\ f_{m-1,j} \end{pmatrix} - q\begin{pmatrix} u_{1,j-1}\\ u_{2,j-1}\\ \vdots\\ u_{m-2,j-1}\\ u_{m-1,j-1} \end{pmatrix} - q\begin{pmatrix} u_{1,j+1}\\ u_{2,j+1}\\ \vdots\\ u_{m-2,j+1}\\ u_{m-1,j+1} \end{pmatrix}\\ -p\begin{pmatrix} u_{0,j}\\ 0\\ \vdots\\ 0\\0 \end{pmatrix}-p\begin{pmatrix} 0\\ 0\\ \vdots\\ 0\\ u_{m,j} \end{pmatrix}\tag{15}$
这样的矩阵方程有 $n-1$ 个，此时便可以将每个向量换成上面的记法，
$\mathbf{U}\vec{u}_{i,j} = \vec{f}_{i,j} - \mathbf{V}\vec{u}_{i,j-1} - \mathbf{V}\vec{u}_{i,j+1} - p\vec{u}_{0,j} - p\vec{u}_{m,j}\tag{16}$
此时便可以固定 $i$ ，将 $i$ 看成常量再转为一个矩阵方程
$\begin{bmatrix} \mathbf{U}&\mathbf{V}\\ &\mathbf{V}&\mathbf{U}&\mathbf{V}\\ &&\ddots&\ddots&\ddots\\ &&&\mathbf{V}&\mathbf{U}&\mathbf{V}\\ &&&&\mathbf{V}&\mathbf{U} \end{bmatrix}\begin{pmatrix} \vec{u}_{i,1}\\ \vec{u}_{i,2}\\ \vdots\\ \vec{u}_{i,n-2}\\ \vec{u}_{i,n-1} \end{pmatrix} = \begin{pmatrix} \vec{f}_{i,1}\\ \vec{f}_{i,2}\\ \vdots\\ \vec{f}_{i,n-2}\\ \vec{f}_{i,n-1}\\ \end{pmatrix} - p\begin{pmatrix} \vec{u}_{0,1}\\ \vec{f}_{0,2}\\ \vdots\\ \vec{f}_{0,n-2}\\ \vec{f}_{0,n-1}\\ \end{pmatrix} - p\begin{pmatrix} \vec{u}_{m,1}\\ \vec{f}_{m,2}\\ \vdots\\ \vec{f}_{m,n-2}\\ \vec{f}_{m,n-1}\\ \end{pmatrix}\\ - q\begin{pmatrix} \vec{u}_{i,0}\\ 0\\ \vdots\\ 0\\ 0 \end{pmatrix} - q\begin{pmatrix} 0\\0\\\vdots\\0\\\vec{u}_{i,m} \end{pmatrix}\tag{17}$
至此便得到了该二阶微分方程的解。

误差分析

对于的二维泊松方程，取 $f(x,y) = -2\pi^2\sin(\pi x)\sin(\pi y),c(y)=d(y)=g(x)=h(x)=0,(x,y)\in [0,1]\times [0,1]$ ，此时有对应的解析解 $u(x,y) = \sin(\pi x)\sin(\pi y)$ ，令 $m=n=100$ 时解析解和FDM得到的解对比如图2所示。

结果

代码

import numpy as np
import matplotlib.pyplot as plt
from numpy import *

M, N = 100, 100
a, b = 1, 1  
hx, hy = a/M, b/N
p, q = 1/hx**2, 1/hy**2
r = -2*(p + q)

U = zeros((M-1, M-1))
for i in range(M-1):
    U[i, i] = r
    if i < M-2: U[i, i+1] = p
    if i > 0:   U[i, i-1] = p 
V = diag([q]*(M-1))
Zero_mat = zeros((M-1, M-1))

A_blc = empty((N-1, N-1), dtype=object) # 矩阵A的分块形式
for i in range(N-1):
    for j in range(N-1):
        if i == j:
            A_blc[i, j] = U
        elif abs(i-j) == 1:
            A_blc[i, j] = V
        else:
            A_blc[i, j] = Zero_mat
            
A = vstack([hstack(A_i) for A_i in A_blc]) # 组装得到矩阵A

x_i = linspace(0, a, M+1)
y_i = linspace(0, b, N+1)
F = vstack([-2*pi**2 * sin(pi*x_i[1:M].reshape((M-1, 1))) * sin(pi*j) for j in y_i[1:N]])

u = dot(linalg.inv(A), F).reshape(M-1, N-1)
u_f = vstack([zeros((1,M+1)),  # 最后组装边界条件得到全域的解
            hstack([zeros((N-1,1)), u, zeros((N-1,1))]), 
            zeros((1,M+1))])

nx, ny = 101, 101
# 创建网格
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
X, Y = np.meshgrid(x, y)
# 计算函数值
u = np.sin(np.pi * X) * np.sin(np.pi * Y)

fig, axs = plt.subplots(nrows=1,ncols=2)

# 绘制电势图
axs[0].imshow(u_f, extent=[0, 1, 0, 1])
fig.colorbar(axs[0].images[-1], ax=axs[0], location='right')
axs[0].set_xlabel('X')
axs[0].set_ylabel('Y')
axs[0].set_title('FDM solution')

axs[1].imshow(u, extent=[0, 1, 0, 1])
fig.colorbar(axs[1].images[-1], ax=axs[1], location='right')
axs[1].set_xlabel('X')
axs[1].set_ylabel('Y')
axs[1].set_title('Analytic solution')

plt.tight_layout()
plt.show()

可以看出几乎没有什么误差，同时还可以来分析一下不同 $m,n$ 取值时的误差，这里选择使用L2距离计算误差，取 $m=n=s$ ，得到的误差变化如图3所示。

地震正演

雷克子波

雷克（Ricker）子波是地震子波的一种，是由震源激发、经地下传播再被地面捕获的地震波。在FWI中，通常用地震子波于反射系数的卷积表示地震数据。
其表达式为
$s(t) = e^{-(\pi f t)^2}(1-2(\pi f t)^2)$
图像
图4. 雷克子波
代码
```
import numpy as np
import matplotlib.pyplot as plt

def ricker_wavelet(t, f):
    pi = np.pi
    gamma = 3
    return np.exp(-(2 * pi * f / gamma) ** 2 * t ** 2) * np.cos(2*pi * f * t)

dt = 0.001  # 时间步长
t = np.arange(0, 0.2, dt)
f0 = 25  # 主频25Hz

ricker_wave = ricker_wavelet(t, f0)
plt.figure(figsize=(10, 4))
plt.plot(t, ricker_wave, label='Ricker wavelet')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.title('Ricker Wavelet')
plt.grid(True)
plt.legend()
plt.show()
```
这是雷克子波的一个周期，我们可以将输出范围调节至[-0.2, 0.2]来查看雷克子波的全貌。
图5. 雷克子波对称区间图像
在正演时我们将使用雷克子波作为初始震源，然后通过求解波动方程得到在指定的速度模型中波信号的变化。

波动方程

首先先看一下声波波动方程：
$\frac{1}{v^2(x, z)} \frac{\partial^2 p(x, z, t)}{\partial t^2} = \frac{\partial^2 p(x,z,t)}{\partial x^2} + \frac{\partial^2 p(x,z,t)}{\partial z^2}\tag{19}$
其中 $(x,z)$ 代表空间坐标， $p$ 代表波场， $v$ 表示速度场， $t$ 代表时间。而我们要做的便是求解该微分方程从速度 $v$ 反解得到 $p$ 。
由二阶微分中心差分公式有
$\begin{aligned} \frac{\partial^2 p}{\partial x^2} = \frac{p(x+1,z,t) - 2p(x,z,t) + p(x-1,z,t)}{h_x^2}\\ \frac{\partial^2 p}{\partial z^2} = \frac{p(x,z+1,t) - 2p(x,z,t) + p(x,z-1,t)}{h_z^2}\\ \frac{\partial^2 p}{\partial t^2} = \frac{p(x,z,t+1) - 2p(x,z,t) + p(x,z,t-1)}{h_t^2} \end{aligned}\tag{20}$
代入原方程得
$\frac{1}{v^2(x, z)} \frac{p(x,z,t+1) - 2p(x,z,t) + p(x,z,t-1)}{h_t^2}\tag{20} = \\\frac{p(x+1,z,t) - 2p(x,z,t) + p(x-1,z,t)}{h_x^2} + \frac{p(x,z+1,t) - 2p(x,z,t) + p(x,z-1,t)}{h_z^2}$
整理可得
$\begin{aligned} p(x,z,t+1) &= 2p(x,z,t) - p(x,z,t-1)\\ &+\frac{v^2(x, z)h_t^2}{h_z^2}(p(x,z+1,t) - 2p(x,z,t) + p(x,z-1,t))\\ &+\frac{v^2(x, z)h_t^2}{h_x^2}(p(x+1,z,t) - 2p(x,z,t) + p(x-1,z,t)) \end{aligned}$
这样我们便得到了在地震正演中如何任意位置上一时刻到下一时刻的波信号变化。

单炮数据正演

首先需要生成速度模型，这里我直接用了链接1^[2]的模拟速度模型(vmodel20.mat)。

然后我们编写python代码进行正演，然后我们可以得到如下结果

结果

代码

import numpy as np
import matplotlib.pylab as plt
import matplotlib.animation as animation
import scipy.io

data = scipy.io.loadmat(r'vmodel20.mat')

velocity_model = data['vmodel']

dx = dz = 10.0
dt = 0.001
nz, nx = velocity_model.shape
nt = 800

wf = np.zeros((nz, nx, nt))
source_x, source_z = nx // 2, 1
f0 = 25.0
source_time = np.zeros(nt)
t = np.arange(0, nt * dt, dt)

source_time = np.exp(-(np.pi * f0 * t) ** 2) * (1-2*(np.pi*f0*t)**2)

for t in range(2, nt-1):
    wf[source_z, source_x, t] += source_time[t]
    for i in range(1, nx - 1):
        for j in range(1, nz - 1):
            v = velocity_model[j, i]

            wf[j, i, t+1] = 2*wf[j, i, t] - wf[j, i, t-1] \
                            + v**2 * dt**2 / (dz**2) * (wf[j+1, i, t] + wf[j-1, i, t] - 2*wf[j, i, t]) +  \
                            + v**2 * dt**2 / (dx**2) * (wf[j, i+1, t] + wf[j, i-1, t] - 2*wf[j, i, t])

fig, ax = plt.subplots()

im = ax.imshow(wf[:, :, 0], extent=[0, nx, nz, 0], cmap='seismic', aspect='auto')
ax.set_xlabel('X (m)')
ax.set_ylabel('Z (m)')
ax.set_title('Wavefield Animation')
fig.colorbar(im, ax=ax, label='Amplitude')

# 初始化函数，设置初始状态
def init():
    im.set_data(wf[:, :, 0])
    return im,
# 更新函数，用于每一帧的更新
def update(frame):
    im.set_data(wf[:, :, frame])
    return im,
# 创建动画
ani = animation.FuncAnimation(fig, update, frames=range(1, nt),
                            init_func=init, blit=True, interval=20)
# 显示图形和动画
plt.show()

波场快照即 $t$ 时刻波的传播情况，这里我们用叠加动画的方式进行展示。可以看到震源处的波随着时间而传播开，在通过不同的介质时也获得了不同的速度得到不同的波形以及形成反射波，但是会发现波在碰到边界后也产生反射波，这是因为正常的波在地下传播是无限的，但我们人为截取了一个边界导致边界反射问题，这个问题将在后续单独的篇章中记录。

这里得到了波场快照之后，我们便可以设置检波点来模拟实际获取数据的效果，因为波场是我们使用波动方程推导的，而实际获取数据时需要在地面或地下设置具体的检波器来获取数据。检波器的设置很简单，获取固定位置的波场数据即可。

结果

代码

import numpy as np
import matplotlib.pylab as plt
import matplotlib.animation as animation
import scipy.io

data = scipy.io.loadmat(r'vmodel20.mat')

velocity_model = data['vmodel']

dx = dz = 10.0
dt = 0.001
nz, nx = velocity_model.shape
nt = 2000

wf = np.zeros((nz, nx, nt))
source_x, source_z = nx // 2, 1
f0 = 25.0
source_time = np.zeros(nt)
t = np.arange(0, nt * dt, dt)

source_time = np.exp(-(np.pi * f0 * t) ** 2) * (1-2*(np.pi*f0*t)**2)

rnums = nx
rx = np.linspace(0, nx-1, rnums, dtype=int)
rz = 1

seismic_records = np.zeros((rnums, nt))


for t in tqdm(range(2, nt-1)):
    wf[source_z, source_x, t] += source_time[t]
    for i in range(1, nx - 1):
        for j in range(1, nz - 1):
            v = velocity_model[j, i]
            wf[j, i, t+1] = 2*wf[j, i, t] - wf[j, i, t-1] \
                            + v**2 * dt**2 / (dz**2) * (wf[j+1, i, t] + wf[j-1, i, t] - 2*wf[j, i, t]) +  \
                            + v**2 * dt**2 / (dx**2) * (wf[j, i+1, t] + wf[j, i-1, t] - 2*wf[j, i, t])
            
    for i, _rx in enumerate(rx):
        seismic_records[i, t+1] = wf[rz, _rx, t+1]
fig, ax = plt.subplots(figsize=(6.2,8), dpi=120)

seismic_records = np.rot90(seismic_records, k=3)

im = ax.imshow(seismic_records, extent=[0, nx, nt, 0], cmap='seismic',aspect='auto', vmin=-0.4, vmax=0.44)
ax.set_xlabel('X (m)')
ax.set_ylabel('time')
ax.set_title('Wavefield Animation')
plt.show()

波动方程的有限差分解（地震数据正演）

# 波动方程的有限差分解（地震数据正演）

# 有限差分法

# 差分公式

# 偏微分方程

# 误差分析

# 地震正演

# 雷克子波

# 波动方程

# 单炮数据正演

# 参考

波动方程的有限差分解（地震数据正演）

有限差分法

差分公式

偏微分方程

误差分析

地震正演

雷克子波

波动方程

单炮数据正演

参考