LSM与Born Modeling
LSM Born-Modeling
这篇文章会介绍基于有限差分方法进行波动方程模拟，同时使用Born近似的方式求解小扰动下的反射波场。基于这两种方法，进一步介绍RTM以及LSM方法的原理与代码实现。
1 波动方程与有限差分方法
对于最简单的声波方程情况：
\frac{1}{s^2} \nabla^2p + f = \frac{\partial^2 p}{\partial t^2}
有限差分方法会将这个方程中的空间微分用空间差分替代，在对时间二阶导使用同样的差分之后，就可以生成我们想要的迭代式。
空间微分近似表达式
公式在不同精度要求下对应着不同的空间差分格式，所以我们需要求得对应的待定系数。
- $f”(x)=\sum_{-n}^{n}c_{i}f(x+i\Delta x) + o(\Delta x ^{2n+1})$，写出对应的精度表达式。
- 代入f(x+i\Delta x)所对应的不同的泰勒展开式
- 得到系数矩阵的方程（下面展示的）c_0=-2 \sum_{i=1}^{n}c_i
\frac{1^2}{2!} & \frac{2^2}{2!} & \ldots & \frac{n^2}{2!}\\
\frac{1^4}{4!} & \frac{2^4}{4!} & \ldots & \frac{n^4}{4!}\\
\vdots & \vdots & \ddots & \vdots\\
\frac{1^{2n}}{2n!} & \frac{2^{2n}}{2n!} & \ldots & \frac{n^{2n}}{2n}
\end{pmatrix}
\begin{pmatrix}
c_1\\
c_2\\\vdots
\\ c_n
\end{pmatrix} =
\begin{pmatrix}
1\\0\\ \vdots\\0
\end{pmatrix}
将不同精度的空间差分系数与二阶差分\frac{\partial^2 p}{\partial t^2} = \frac{p(t+1)+p(t-1)-2p(t)}{\Delta t}，我们可以得到：
+(\frac{v\Delta t}{\Delta z})^2\sum_{i=1}^n c_i[p_t(x, z+i\Delta z)+p_t(x, -i\Delta z)] \\
+[(\frac{v\Delta t}{\Delta x})^2c_0+(\frac{v\Delta t}{\Delta z})^2c_0]p_t(x,z)+2p_t(x,z)-p_{t-1}(x,z) \\
=p_{t+1}(x,z)
边界条件与源
上面给出了如何进行时间迭代，但是一个完整的问题还需要提供边界条件与源的信息。对于地球物理问题，常见的有两种边界条件：
- 自由边界条件：这是对于地表情况的特殊处理。设定为\nabla u = p(x,0)=0。常用的处理方式是在上边界往外扩充模型，然后在进行计算的时候，赋予上面的模型与下面的数值相反。这样就相当于在边界处使用了单向泰勒表达式
- 吸收边界条件：这是模拟地下无限空间。防止地下边界处出现强反射现象。一个常用的方法是在边界外添加一个海绵层，在此海绵层中，p(t+1)=(1-\kappa)p(t)。这样就能一定程度上避免边界反射的问题。海绵层的常用设置\kappa=C\times(\frac{d_{i2b}}{D_{boundary}})^2。设置kappa随着边界逐渐变大。
而对于源的问题，只需要直接将源f(x,z,t)直接添加到上述公式，最终可以得到：
+(\frac{v\Delta t}{\Delta z})^2\sum_{i=1}^n c_i[p_t(x, z+i\Delta z)+p_t(x, -i\Delta z)] \\
+[(\frac{v\Delta t}{\Delta x})^2c_0+(\frac{v\Delta t}{\Delta z})^2c_0]p_t(x,z)\\
+(2-\kappa)p_t(x,z)-(1-\kappa)p_{t-1}(x,z)+s_{t}(x)\Delta t=p_{t+1}(x,z) \\
,where\ \kappa=(\frac{d_{i2b}}{D_{abc}})^2\times 3.0 \times 16 v_{min}\times \Delta t.\ s(x,t)\ is\ source
这里的\kappa,我选择用一个学长的设置
稳定性条件
–留空–
LSM 方法与Born Modeling
Born 近似
同样考虑在声波方程中，如果我的慢度场存在一个扰动 \delta s(x):
我们可以获得在进行近似（泰勒展开并忽略二阶极小）之后的表达式（注意p依然是满足波动方程的）
上面的表达式中，我们的原波动场p变成了扰动场q里面的源，二次反射(图里面的q(q))被忽略掉，变成了一个线性的关系（无论是对\delta s\ or\ p）。如果继续往下面计算，并使用格林函数的方法（或者线性矩阵）替换掉波动方程表达式，我们就可以得到\frac{\delta p}{\delta s}的表达式（Fréchet微分）。
LSM
在Fréchet微分的基础上，我们可以获得\Phi =(p-p_0)^2这个最小二乘问题下的梯度表达式，从而开始反演。但是，对于背景场的慢度信息，我们在四个地方需要使用：
- 1 获取正向波场p；这一个波场是从我们一开始设置的源开始的
- 2 获取反射率系数；这里的反射率代指从P场波动到q场源所乘上的系数
- 3 获取扰动波场q；
- 4 获取反传波场（以res为源）
如果我们每一次获取的梯度信息m=\frac{\delta \Phi}{\delta s}，都会更新这四个中的所有慢度信息，这样我们的问题就是一个非线性问题，对应着FWI方法。如果我们每一次都只画出m,不进行迭代，对应的就是RTM方法；而LSM方法是只更新反射率系数中对应的\delta s,不更新获取正传波场与扰动波场中的s(速度信息)，那么我们\frac{\delta \Phi}{\delta s}就是一个线性最小二乘问题，可以通过各种迭代方法进行。
Born Modelling
对于FWI等方法而言，只需要在理论推导中进行Born近似，但是在运算中，只需要计算正传播场p（这一步是123的综合，不需要进行近似）与反传波场（以geophones为源），然后计算双方的互相关即可得到梯度方向。但是对于LSM方法，我们只更新反射率，意味着我们的正传场p实际是不能改变的。这就要求我们记录下正传场，将之作为源，然后利用更改的反射率获取扰动场，这就是Born Modeling。Born Modeling代表利用Born近似，在正传场不变的情况下获取扰动场的方法。
代码实践
首先，我们需要整理基本的算法，这里我用python伪代码表示
import finite_diff, born_modelling
import correlation
# data
data,source = load('...')
# model
s,R0 = load('...')
p_field = finite_diff(sources, s)
for iter in range(iteration):
# get gradient
data_syn = born_modelling(p_field, s, R0)
data_res = data_syn-data
b_field = finite_diff(data_res, s)
m = correlation(p_field, b_field) # gradient or RTM images
R0 = R0 + alpha*m
output(R0)
Tricks
这个方法存在以下的一些问题：
- 当我们使用平滑的地下速度场计算波场时，直达波差异会比较大。由于表层的速度不精确，直达波的到时与振幅