零基础学习正演的数值模拟（含代码）

摘要: 本贴从零开始学习正演的数值模拟方法. 包括相应的偏微分基础、声波方程、雷克子波、均匀速度场的模拟、一般速度场的模拟.

1. 偏微分基础

本小节仅涉及高等数学相关知识, 与领域无关.

1.1 导数

引例: 物体从一维坐标的原点开始移动, 在 $t$ 时刻, 它在坐标轴的位置由函数 $s(t)$ 确定, 则速度为位置变化量与时间的比值:
$$v(t)=\frac{\mathrm{d}s(t)}{\mathrm{d}t}=\underset{\Delta t\to 0}{\lim }\frac{s(t+\Delta t)-s(t)}{\Delta t}\text{\hspace{1em}(1)}$$
加速度为速度变化量与时间的比值:
$$a(t)=\frac{\mathrm{d}v(t)}{\mathrm{d}t}=\underset{\Delta t\to 0}{\lim }\frac{v(t)-v(t-\Delta t)}{\Delta t}=\underset{\Delta t\to 0}{\lim }\frac{s(t+\Delta t)-2s(t)+s(t-\Delta t)}{\Delta t^2}\text{\hspace{1em}(2)}$$

推广 1: 给定一个单变量函数
$$y=f(x)\text{\hspace{1em}(3)}$$
其一阶导数记为
$$y^{\prime }=\frac{\mathrm{d}f(x)}{\mathrm{d}x}\text{\hspace{1em}(4)}$$
二阶导数记为
$$y^{\prime \prime }=\frac{{\mathrm{d}}^{2}f(x)}{\mathrm{d}{x}^2}\text{\hspace{1em}(5)}$$

1.2 偏导

给定一个二变量函数
$$z=f(x,y)\text{\hspace{1em}(6)}$$
其针对 $x$ 偏导的为
$$\frac{\partial z}{\partial x}=\underset{\Delta x\to 0}{\lim }\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}\text{\hspace{1em}(7)}$$
即 $x$ 发生了变化, 而 $y$ 并没变化.
进一步, 二阶偏导为
$$\begin{array}{ll}\frac{{\partial }^{2}z}{\partial x^2}& ={\lim }_{\Delta x\to 0}\frac{\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}-\frac{f(x,y)-f(x-\Delta x,y)}{\Delta x}}{\Delta x}\\ & ={\lim }_{\Delta x\to 0}\frac{f(x+\Delta x,y)-2f(x,y)+f(x-\Delta x,y)}{\Delta x^2}\end{array}\text{\hspace{1em}(8)}$$

另外有 (这两个式子在本贴里面不用, 写着让大家复习):
$$\frac{{\partial }^{2}z}{\partial x\partial y}=\underset{\Delta x\to 0,\Delta y\to 0}{\lim }\frac{f(x+\Delta x,y+\Delta y)-f(x,y+\Delta y)-f(x+\Delta x,y)+f(x,y)}{\Delta x\Delta y}\text{\hspace{1em}(9)}$$
$$\frac{{\partial }^{2}z}{\partial y\partial x}=\frac{{\partial }^{2}z}{\partial x\partial y}\text{\hspace{1em}(10)}$$
在进行数值模拟的时候, 不可能取 $\Delta x\to 0$, 因此 (8) 式简化为
$$\frac{{\partial }^{2}z}{\partial x^2}\approx \frac{f(x+\Delta x,y)-2f(x,y)+f(x-\Delta x,y)}{\Delta x^2}\text{\hspace{1em}(11)}$$

注 1: 为统一起见, 即使一元函数, 以后也常使用 $\partial$ 而不是 $\mathrm{d}$.

1.3 泰勒级数

显然, $\Delta x$ 越接近 $0$ 误差就越小, 但在实际地震数据采集中, 需要的设备就越多. 例如, $\Delta =5\mathrm{m}$, 即每隔 5m 放置一个检波器, 就要花费 $\Delta =20\mathrm{m}$ 时 4 倍的成本, 需要计算的网格点数也增加到后者的 16 倍 (2 维剖面) 或 64 倍 (3 维数据体).
从 (11) 式, 我们无法知道: 这个约等于究竟涉及多大的误差?
为了确切知道误差在哪个量级之内, 我们需要引入更高级的数学工具: 泰勒级数.
当函数 $f(x)$ 在 $x_0$ 处存在直到 $n$ 阶的导数, 则
$$f(x)=f(x_0)+\sum _{i=1}^n\frac{f^{(i)}(x_0)(x-x_0{)}^i}{i!}+o((x-x_0{)}^n)\text{\hspace{1em}(12)}$$
其中 $f^{(i)}(x_0)/i!$ 称为泰勒展开式的系数.
直观的解释: 如果函数 $f(x)$ 不是线性的, 则它的变化量不仅与斜率有关, 而且与斜率的变化率也有关. 更多内容就只有自己去找数学书啃了.
为与我们前面的内容相符, 将 (12) 式 $x_0$ 记作 $x$, $x$ 记作 $x+\Delta x$, 得到
$$f(x+\Delta x)=f(x)+\sum _{i=1}^n\frac{f^{(i)}(x)\Delta x^i}{i!}+o(\Delta x^n)\text{\hspace{1em}(13)}$$

为方便学计算机的同学理解, 来讲几个特例. 数学学院的同学忽略.
例 1. 验证二次函数
$$f(x)=a{x}^2+b\text{\hspace{1em}(14)}$$
$f(x+\Delta x)=(a{x}^2+b)+2ax\Delta x+2a\Delta x^2/2=a{x}^2+2ax\Delta x+a\Delta x^2+b=a(x+\Delta x{)}^2+b$
验证结束.
注意: 泰勒级数现实的意义不在于这种一直连续可导的函数, 而是在某些区域可导的函数.

1.4 2 阶精度

我们不知道 $n$ 的具体取值, 就会假设它比较大. 在数值计算的时候, 会忽略 (13) 式的 $o(\Delta x^n)$ 甚至前面某些累加式. 也就是说, 假设存在 10 阶导数, 但我们仅取 $n=2$ 的话, 误差就被控制在 $o(\Delta x^2)$, 这时候称为 2 阶精度. 仅取 $n=4$ 的话, 误差就被控制在 $o(\Delta x^4)$, 这时候称为 4 阶精度. 显然, 4 阶精度比 2 阶精度省略的值更小, 因此精度更高. 我们可以根据自己的需求, 设置相应的精度. 在地震数据数值模拟中, $n+2$ 阶比 $n$ 阶的误差大概低 1 个数量级. 正演模拟在传播过程中误差会不断累积, 严重的情况下出现“频散” (不应该存在的波, 见图 3), 所以我们还是倾向于多计算一些.
特别地, 仅考虑 2 阶精度的时候
$$f(x+\Delta x)=f(x)+\frac{\partial f}{\partial x}\Delta x+\frac{{\partial }^{2}f}{\partial x^2}\frac{\Delta x^2}{2}+o(\Delta x^2)\text{\hspace{1em}(15)}$$
将 (15) 式移项得到
$$f(x+\Delta x)-f(x)=\frac{\partial f}{\partial x}\Delta x+\frac{{\partial }^{2}f}{\partial x^2}\frac{\Delta x^2}{2}+o(\Delta x^2)\text{\hspace{1em}(16)}$$
同理得到
$$f(x)-f(x-\Delta x)=\frac{\partial f}{\partial x}\Delta x-\frac{{\partial }^{2}f}{\partial x^2}\frac{\Delta x^2}{2}+o(\Delta x^2)\text{\hspace{1em}(17)}$$
这里 $o(\Delta x^2)$ 的符号可正可负.
(16) 式减去 (17) 式, 然后将两边同时除以 $\Delta x^2$ 可得
$$\frac{{\partial }^{2}f}{\partial x^2}=\frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\Delta x^2}+O(\Delta x^2)\text{\hspace{1em}(18)}$$
这里的高阶无穷小除了相应的变量后, 成为同阶无穷小.
注 2: 奇数阶精度不用计算, 例如 3 阶与 2 阶精度的表达式是一样的, 仅把 $O(\Delta x^2)$ 替换为 $O(\Delta x^3)$ 即可.

1.5 4 阶精度

考虑到 4 阶导数的时候
$$f(x+\Delta x)=f(x)+\frac{\partial f}{\partial x}\Delta x+\frac{{\partial }^{2}f}{\partial x^2}\frac{\Delta x^2}{2}+\frac{{\partial }^{3}f}{\partial x^3}\frac{\Delta x^3}{6}+\frac{{\partial }^{4}f}{\partial x^4}\frac{\Delta x^4}{24}+o(\Delta x^4)\text{\hspace{1em}(19)}$$
$$f(x+\Delta x)-f(x)=\frac{\partial f}{\partial x}\Delta x+\frac{{\partial }^{2}f}{\partial x^2}\frac{\Delta x^2}{2}+\frac{{\partial }^{3}f}{\partial x^3}\frac{\Delta x^3}{6}+\frac{{\partial }^{4}f}{\partial x^4}\frac{\Delta x^4}{24}+o(\Delta x^4)\text{\hspace{1em}(20)}$$
$$f(x)-f(x-\Delta x)=\frac{\partial f}{\partial x}\Delta x-\frac{{\partial }^{2}f}{\partial x^2}\frac{\Delta x^2}{2}+\frac{{\partial }^{3}f}{\partial x^3}\frac{\Delta x^3}{6}-\frac{{\partial }^{4}f}{\partial x^4}\frac{\Delta x^4}{24}+o(\Delta x^4)\text{\hspace{1em}(21)}$$
(20) 式减去 (21) 式, 然后将两边同时除以 $\Delta x^2$ 可得
$$\frac{{\partial }^{2}f}{\partial x^2}=\frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\Delta x^2}-\frac{{\partial }^{4}f}{\partial x^4}\frac{\Delta x^2}{12}+o(\Delta x^2)\text{\hspace{1em}(22)}$$
将 $\Delta x$ 换为 $2\Delta x$, 同理可得
$$\frac{{\partial }^{2}f}{\partial x^2}=\frac{f(x+2\Delta x)-2f(x)+f(x-2\Delta x)}{4\Delta x^2}-\frac{{\partial }^{4}f}{\partial x^4}\frac{16\Delta x^2}{12}+o(\Delta x^2)\text{\hspace{1em}(23)}$$
(22) 式乘以 16 减去 (23) 式, 可得
$$\frac{{\partial }^{2}f}{\partial x^2}=\frac{16}{15}\frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\Delta x^2}-\frac{1}{15}\frac{f(x+2\Delta x)-2f(x)+f(x-2\Delta x)}{4\Delta x^2}+O(\Delta x^4)\text{\hspace{1em}(24)}$$
从这里可以看到, 通过引入 $2\Delta x$, 可以消去 4 阶偏导. 这是增加精度的核心技巧.

1.6 $2n$ 阶精度

通过前面的“核心技巧”, 将 (24) 式进一步推广, 可得 $2n$ 阶精度的表达式
$$\frac{{\partial }^{2}f}{\partial x^2}=\sum _{i=1}^n(-1{)}^{i+1}{c}_i\frac{f(x+i\Delta x)-2f(x)+f(x-i\Delta x)}{i^2\Delta x^2}+O(\Delta x^{2n})\text{\hspace{1em}(25)}$$

其中:

• 系数 $c_1,\dots ,c_n$ 没有给出理论值. 在实际工作中, 由于 $\Delta x$ 比较大, 所以要专门计算系数, 它们与差分格式有关. 也是这个方向的重要研究内容. 表 1 给出了中间差分格式的系数.
• 误差为 $O(\Delta x^{2n})$, 即 $n$ 越大误差越小, 计算量也越大 (不知道是否可以用 GPU 计算, 速度就会增快很多).
• $n$ 越大就涉及更远的点, 如果实际应用中的数据并没有对应那么高阶可导的函数, 效果不一定有那么好. 不过越远的点所对应的系数越小, 影响也没那么大.

2. 波动方程

2.1 弦振动 (横波) 方程

参见全波形反演的深度学习方法: 第 2 章 正演, 根据牛顿第二定律
$$F=ma\text{\hspace{1em}(26)}$$
弦振动方程为
$$\frac{{\partial }^{2}u(x,t)}{\partial t^2}=c^2\frac{{\partial }^{2}u(x,t)}{\partial x^2}+f(x,t)\text{\hspace{1em}(27)}$$
其中 $c^2=T/\rho$, $f(x,t)=F(x,t)/\rho$, 左式的物理意义是瞬时加速度 $a$, 右式第一项的物理意义是 单位质量所受的力 $F$, $c$ 的物理意义是速度.

进一步忽略重力 $F(x,t)$ 的作用, 可以推出一维齐次波动方程的解:
$$\frac{{\partial }^{2}u(x,t)}{\partial x^2}=\frac{1}{c^2}\frac{{\partial }^{2}u(x,t)}{\partial t^2}\text{\hspace{1em}(28)}$$

2.2 声波 (纵波) 方程

声波仅有纵波. 考虑二维的情况, 它满足
$$\frac{1}{v^2}\frac{{\partial }^{2}U}{\partial t^2}=\frac{{\partial }^{2}U}{\partial x^2}+\frac{{\partial }^{2}U}{\partial z^2}\text{\hspace{1em}(29)}$$
其中 $U$ 指压力. 注意该式是物理定律, 不是从其它式子推导过来的.

图 1 矩阵网格剖分

为了便于数值模拟, 将平面进行离散化, 仅考虑某些网格交叉点, 质量、压力等仅存在于这些点 (称为质点, 不知是否专业). 这样, 我们只考察第 $i$ 行第 $j$ 列的质点在时间 $k$ 的压力
$$U_{i,j}^k\text{\hspace{1em}(30)}$$
对于 2 阶精度, 将 (18) 式按照变量名改造后代入 (28) 式可得
$$\frac{1}{v^2}\frac{U_{i,j}^{k+1}-2U_{i,j}^k+U_{i,j}^{k-1}}{\Delta t^2}=\frac{U_{i+1,j}^k-2U_{i,j}^k+U_{i-1,j}^k}{\Delta x^2}+\frac{U_{i,j+1}^k-2U_{i,j}^k+U_{i,j-1}^k}{\Delta y^2}\text{\hspace{1em}(31)}$$
其中 $k+1$ 表示下一个时间点, $i+1$ 表示下一个质点.
对于 4 阶精度, 将 (24) 式按照变量名改造后代入 (28) 式可得
$$\begin{array}{ll}\frac{1}{v^2}\frac{U_{i,j}^{k+1}-2U_{i,j}^k+U_{i,j}^{k-1}}{\Delta t^2}=& c_1\frac{U_{i+1,j}^k-2U_{i,j}^k+U_{i-1,j}^k}{\Delta x^2}+c_1\frac{U_{i,j+1}^k-2U_{i,j}^k+U_{i,j-1}^k}{\Delta y^2}\\ & +c_2\frac{U_{i+2,j}^k-2U_{i,j}^k+U_{i-2,j}^k}{\Delta x^2}+c_2\frac{U_{i,j+2}^k-2U_{i,j}^k+U_{i,j-2}^k}{\Delta y^2}\\ & +O(\Delta x^4)\end{array}\text{\hspace{1em}(32)}$$
继续推广可获得 $2n$ 阶精度的方程
$$\frac{1}{v^2}\frac{U_{i,j}^{k+1}-2U_{i,j}^k+U_{i,j}^{k-1}}{\Delta t^2}=\sum _{m=1}^n\left(c_m\frac{U_{i+m,j}^k-2U_{i,j}^k+U_{i-m,j}^k}{\Delta x^2}+c_m\frac{U_{i,j+m}^k-2U_{i,j}^k+U_{i,j-m}^k}{\Delta y^2}\right)+O(\Delta x^{2n})\text{\hspace{1em}(32)}$$

利用式 (31) - (33), 根据 $k$ 与 $k-1$ 时刻各网格点的声压, 可以计算出 $k+1$ 时刻各网格点的声压. 这样我们就可以正演啦.

3. 代码与结果

3.1 雷克子波

振源需要产生一个波, 它持续一小段时间.

tic
clc
close all
clear all

end_time = 0.5; % Total simulation time, in seconds
delta_t = 0.0005; % Time step, in seconds
f0 = 30; % The wave frequency, in 10~40 Hz

ricker_wave = zeros(end_time / delta_t + 1);
i = 1;
for time = 0: delta_t: end_time
    ricker_wave(i) = 5.76 * f0^2 * (1 - 16 * (0.6 * f0 * time - 1)^2) * exp(-8 * (0.6 * f0 * time-1)^2);
    i = i + 1;
end

plot(ricker_wave);

图 1. 雷克子波.

根据代码生成的子波如图 1 所示. 用它可以模拟振源仅波动一次的情况, 从开始小的振幅, 到最大振幅, 又来一个小振幅, 就形成了一个完整的波. 当 time 比较大的时候, 波动几乎为 0.

3.2 声波 2 阶精度

根据 (31) 式可以写出第一个声波方程所对应的模拟过程.

% acousti_wave.m
% Forward simulation of acoustic wave.
tic
clc
close all
clear all

end_time = 0.5; % Total simulation time, in seconds
delta_t = 0.0005; % Time step, in seconds
delta_x = 6; % Space step in the X direction, in meters
delta_z = 6; % Space step in the Z direction, in meters

cnx = 301; % Number of grids in X direction
cnz = 301; % Number of grids in Z direction

v = 1500; % The velocity of the wave, in meters/second
sx = (cnx + 1)/2; % The X position of the wave source
sz = (cnz + 1)/2; % The Z position of the wave source

f0 = 30; % The wave frequency, in 10~40 Hz

% Initialization
u_now = zeros(cnx, cnz); % The pressure at the current moment, it is a matrix for all points
u_prev = zeros(cnx, cnz); % The pressure at the previous moment
u_next = zeros(cnx, cnz); % The pressure at the next moment

% Simulate
for time = 0: delta_t: end_time
  for i = 3: cnx - 2
    for j = 3: cnz - 2
      % Implement Eq. (31)
      part1 = (-2 * u_now(i, j) + u_now(i + 1, j) + u_now(i - 1, j)) / delta_x^2;
      part2 = (-2 * u_now(i, j) + u_now(i, j + 1) + u_now(i, j - 1)) / delta_z^2;
      u_next(i, j) = 2 * u_now(i, j) - u_prev(i, j) + v^2 * (part1 + part2) * delta_t^2;
    end
  end
  % The wave at the source
  u_next(sx, sz) = 5.76 * f0^2 * (1 - 16 * (0.6 * f0 * time - 1)^2) * exp(-8 * (0.6 * f0 * time-1)^2);
  % Update all points
  u_prev = u_now;
  u_now = u_next;
end

% Paint
surf(u_now)
shading interp;
view(2); %view(90,90)
colormap(gray);

toc

图 2. 波场快照.

图 2 给出了 0.5 秒时的波场快照. 和池塘里丢一颗石头相似, 振源来自于 (151, 151), 因此 0.5 秒时波传播到外面 黑、白、黑三个圈依次对应于图 1 的振幅负、正、负.

3.3 声波 4 阶精度时的频散

修改步长、网格点参数.

% acousti_wave.m
% Forward simulation of acoustic wave.
tic
clc
close all
clear all

end_time = 0.5; % Total simulation time, in seconds
delta_t = 0.0005; % Time step, in seconds
delta_x = 12; % Space step in the X direction, in meters
delta_z = 12; % Space step in the Z direction, in meters

cnx = 151; % Number of grids in X direction
cnz = 151; % Number of grids in Z direction

v = 1500; % The velocity of the wave, in meters/second
sx = (cnx + 1)/2; % The X position of the wave source
sz = (cnz + 1)/2; % The Z position of the wave source

f0 = 30; % The wave frequency, in 10~40 Hz

c_1 = 9/8; % 16/15, 1
c_2 = -1/24; % -1/15, 0
% Initialization
u_now = zeros(cnx, cnz); % The pressure at the current moment, it is a matrix for all points
u_prev = zeros(cnx, cnz); % The pressure at the previous moment
u_next = zeros(cnx, cnz); % The pressure at the next moment

position_x = 70;
position_z = 70;
temp_wave = zeros(end_time / delta_t + 1);

k = 1;
% Simulate
for time = 0: delta_t: end_time
  for i = 3: cnx - 2
    for j = 3: cnz - 2
      % Implement Eq. (31)
      part1_1 = (-2 * u_now(i, j) + u_now(i + 1, j) + u_now(i - 1, j)) / delta_x^2;
      part1_2 = (-2 * u_now(i, j) + u_now(i, j + 1) + u_now(i, j - 1)) / delta_z^2;
      part2_1 = (-2 * u_now(i, j) + u_now(i + 2, j) + u_now(i - 2, j)) / delta_x^2;
      part2_2 = (-2 * u_now(i, j) + u_now(i, j + 2) + u_now(i, j - 2)) / delta_z^2;
      
      parts = c_1 * (part1_1 + part1_2) + c_2 * (part2_1 + part2_2);
      u_next(i, j) = 2 * u_now(i, j) - u_prev(i, j) + v^2 * parts * delta_t^2;
    end
  end
  % The wave at the source
  u_next(sx, sz) = 5.76 * f0^2 * (1 - 16 * (0.6 * f0 * time - 1)^2) * exp(-8 * (0.6 * f0 * time-1)^2);
  % Update all points
  u_prev = u_now;
  u_now = u_next;
  
  temp_wave(k) = u_next(position_x, position_z);
  k = k + 1;  
end

plot(temp_wave);

toc

图 3. 采用系数 $(1, 0)$, 即只保留 2 阶精度的波形. 图 4. 采用系数 $(16/15, -1/15)$, 即理论保留 4 阶精度的波形. 图 5. 采用系数 $(9/8, -1/24)$, 即优化系数的波形.

比较图 3至图 5 可以发现, 使用更好的系数, 可以一定程度压制多余的子波, 即频散. 但要接近图 1 (振源) 的效果, 还需要其它的方法.

3.4 运用于一般的波场数据

图 2 这种波场快照仅仅是一个圆圈. 对于实际的数据, 不同点的速度是不一样的, 代码会复杂很多吧?
No no no !!!
只需要把前面代码里的代码 v 替换成一个二维数组

v = zeros(cnx, cnz);

然后从数据文件里面将它读入. 其它的代码不需要改, (29) 式的声波方程已经把所有东西都考虑了.
有没有再次被物理雷到?

图 6. 波场快照随着时间的改变, 以及地震数据的采集.

图 6 由张星移同学给出. 振源来自于中间顶部, 这符合我们在地球表面放炮的设定. 如左图所示, 波场快照开始的时候是一个半圆, 然后到深度 200 左右的时候遇到地层速度变化 (可能从泥土层到达花岗石层), 出现在反射、干涉等. 右图则记录了检波器 (均匀布置在地面上) 所获得的时序信号.

4. 小结

正演模拟看起来高大上, 但入门也没那么困难. 当然, 你要做深入了, 还是高大上!