FWI传统方法

LBFGS

1. 目标：解决无约束优化问题：$\min_{x\in\mathbb{R}^n} f(x)$

2. 牛顿法：

$$f(x_k+p)\approx f(x_k)+\nabla f(x_k)^Tp+\frac{1}{2}p^T\nabla^2f(x_k)p$$

最小化该二次模型，可以得到搜索方向$p_k=-[\nabla^2f(x_k)]^{-1}\nabla f(x_k)$，然后逐步更新$x_{k+1}=x_k+p_k$。该方法需要计算和存储Hessian矩阵和逆，计算复杂度$O(n^3)$ 3. *拟牛顿法：构造Hessian矩阵逆的近似$H_k\approx[\nabla^2f(x_k)]^{-1}$，理想情况下希望它满足拟牛顿条件：

$$H_{k+1}(x_{k+1}-x_k)\approx\nabla f(x_{k+1})-\nabla f(x_k)$$

4. BFGS:

$$H_{k+1}=(I-\rho_ks_ky_k^T)H_k(I-\rho_ks_ky_k^T)+\rho_ks_ks_k^T$$

其中$s_k=x_{k+1}-x_{k},y_k=\nabla f(x_{k+1})-\nabla f(x_k),\rho_k=\frac{1}{y_k^Ts_k}$，避免计算Hessian，但仍然需要存储稠密的$n\times n$的矩阵$H_k$

5. LBFGS: 目标是计算$p_k=-H_k\nabla f(x_k)$，而不显式的存储$H_k$。迭代流程如下:
   • 保存最近m次的迭代对$\{s_i,y_i\}$，初始逆矩阵近似$H_k^0$：通常为简单的标量矩阵，例如$H_K^0=\gamma_kI,\gamma_k=\frac{s_{k-1}^Ty*{k-1}}{y^T*{k-1}y\_{k-1}}$
   • 初始化$q=\nabla f(x_k)$
   • 对于$i=k-1, k-2, \cdots, k-m$，令$\alpha_i=\rho_is_i^Tq,q=q-\alpha_iy_i,r=H_k^0q$
   • 对于$i=k-m, k-m+1, \cdots, k-1$，令$\beta_i=\rho_iy_i^Tr,r=r+s_i(\alpha_i-\beta_i)$，最终得到的r就是我们想要的$H_k\nabla f(x_k)$，即$p_k=-r$该方法的收敛速度为超线性收敛，空间复杂度为O(mn)

Wasserstein-p距离

1. 定义：假设地面上$\mathcal{X} = \mathbb{R}^2$堆了一些石子，石子的分布我们用$\mu : \mathcal{X} \to \mathbb{R}$来表示，采取这样的表示方法对于地面上的任意一块面积$A \subseteq \mathcal{X}$，$\mu(A)$表示这块面积上放置了质量为多少的石子。同样的我们可以定义目标石子堆的分布$\nu$。定义一个输运方案$T: \mathcal{X} \to \mathcal{X}$把现有的石子堆变成目标石子堆。$T(A) = B$表示把原来放在A处的石子都运到B处放好，类似地可以定义反函数$T^{-1}(B) = A$。该输运方案成立需要满足

$$\nu(B) = \mu(T^{-1}(B)), \quad \forall B \subseteq \mathcal{X}$$

即任意位置的石子通过输运过后都刚好满足分布$\mu$的要求。这也可以写为$T\#\mu = \nu$。由此，两堆石子之间的距离可以被定义成把一堆石子挪动或另外一堆所需要的最小输运成本

$$W_p(\mu, \nu) = \left(\inf_{T: T\#\mu = \nu} \int_{\mathcal{X}} \|x - T(x)\|^p \mu(x)dx\right)^{1/p}$$

对于石子不可分的情况，我们有如下公式

$$W_2(\mu, \nu) = \left(\inf_{\gamma \in \Gamma(\mu, \nu)} \int_{X \times Y} d(x,y)^2 d\gamma(x,y) \right)^{1/2}$$

1. $\Gamma(\mu, \nu)$是所有联合分布的集合，其边缘分布分别为$\mu$和$\nu$, $d(x,y)$是x和y之间的距离度量, $\gamma$表示如何将分布$\mu$"传输"到分布$\nu$. 一维情况下有如下公式：

$$W_2^2(P,Q)=\int_0^1|F^{-1}(u)-G^{-1}(u)|^2du$$

其中$F$和$G$为相应的概率分布函数

2. 一维情形下我们有如下的闭式解：

$$W_p^p(P,Q)=\int_0^1|F^{-1}(u)-G^{-1}(u)|^pdu$$

F和G分别为与P，Q对应的分布函数。离散情况下，假设我们有从P中取出的N个样本$\{x_1,\cdots,x_n\},$从Q中取出的N个样本$\{y_1,\cdots,y_n\}$,我们有经验分布函数：

$$FN(x)=\frac{1}{N}\sum_{i=1}^N \mathbb{1}_{\{x_i<x\}}$$

相应的有经验分位数函数

$$FN^{-1}(u)=x_{(k)},\frac{k-1}{N}<u\le\frac{k}{N}$$

同理可得$G^{-1}(u)$。由此我们有

$$W_p^p(P_N,Q_N)=\int_0^1|F_N^{-1}(u)-G_N^{-1}(u)|^pdx=\sum_{k=1}^N\int_{(k-1)/N}^{k/N}|x_{(k)}-y_{(k)}|^p du=\frac{1}{N}\sum_{k=1}^N|x_{(k)}-y_{(k)}|^p$$

3. 二维情形可以采用如下两种方法: 1. Sliced Wasserstein Distance：将二维分布做投影到一维分布，然后对一维分布投影计算相应的W2距离，然后做积分。具体的数学公式为

$$SW_p(P,Q)=(\int_{\mathbb{S}^{d-1}}W_p^p(P_\theta,Q_\theta)d\theta)^{1/p}$$

其中$P_\theta$为$P$在方向$\theta$上的一维投影分布：

$$P_\theta=\text{Distribution of } \theta^TX,\; X\sim P$$

实践中我们采用随机采样方向来近似这个积分

$$SW_p(P,Q)=(\frac{1}{L}\sum_{l=1}^LW_p^p(P_\theta,Q_\theta)d\theta)^{1/p}$$

也即我们先随机生成L个随机方向$\theta_1,\cdots,\theta_L\in \mathcal{S}^{d-1}$, 然后将P中样本$\{x_1,\cdots,x_N\}\subset\mathbb{R}^d$，Q中样本$\{y_1,\cdots,y_N\}\subset \mathbb{R}^d$分别在这L个方向上做投影（也即做点乘），得到的投影值即为新的一维分布，这样就可以按照上面的一维情形求解了2. sinkhorn：引入熵正则化。原始的离散Wasserstein 距离需要解线性规划：

$$W_p^p(P,Q)=\min_{\gamma\in\Gamma(P,Q)}\sum_{i=1}^m\sum_{j=1}^n\gamma_{ij}C_{ij}$$

其中$C_{ij}=|x_i-y_j|^p,\sum_j\gamma_{ij}=p_i,\sum_i\gamma_{ij}=q_j,\gamma_{ij}\ge0$。sinkhorn引入熵正则化：

$$W_{p,\epsilon}^p(P,Q)=\min_{\gamma\in\Gamma(P,Q)}[(\sum_{i=1}^m\sum_{j=1}^n\gamma_{ij}C_{ij})-\epsilon H(\gamma)]$$

其中熵$H(\gamma)=-\sum_{i,j}\gamma_{ij}(\log\gamma_{ij}-1)$，该问题有显式解

$$\gamma=\text{diag}(u)K\text{diag}(v)$$

其中$K_{ij}=\exp(-C_{ij}/\epsilon)$，将该解代入约束可得

$$\sum_j u_i K_{ij}v_j =p_i\Rightarrow u_i=\frac{p_i}{\sum_j K_{ij}v_j}$$ $$\sum_i u_i K_{ij}v_j =q_j\Rightarrow v_j=\frac{q_j}{\sum_i K_{ij}u_i}$$

由此可以得到迭代公式$u\leftarrow \frac{p}{Kv},v\leftarrow \frac{q}{K^Tu}$

关于误差函数

波场是难以拟合好的，根本原因在于波场容易错开，具体来讲我们可以看下面的例子：

在这张图片中，我们可以看到橙色的波场是蓝色波场的一个偏移，在拟合的时候，梯度下降会优先将波做压缩而非平移，他并不知道平移的结果会使得误差函数降得更低。并且当偏移量足够大的时候，会发现MSE几乎不动，比如下面的图：

我们可以看到的是我们更想要橙色的波场，因为离蓝色的波场更近，但是计算橙色和绿色的波场的MSE都是差不多的，我们可以画出偏移量和MSE的关系图：

可以看到稍微大点的偏移量MSE几乎没啥改变。所以我们要想办法修改一下MSE。

如果波场都是正的，是不是我直接沿着x轴做cumsum就可以了？所以我尝试了一下取绝对值再然后再累积求和，具体的数学公式如下：假设我们有观测波场$u_o(t)$，合成波场$u_p(t)$（波场均为一维的），我们定义如下的距离：

$$d(u_o,u_p)=\int_0^T(\int_0^t|u_o(s)|-|u_p(s)|ds)^2dt$$

果然解决了这个问题：

画出来的误差图像如下：

这个操作能解决波场偏移的问题，波峰的高低可能可以再用什么trick处理一下？实际实验似乎有那么一点效果，但应该还能改进（用这个误差dps好像试出来是有戏的）。有什么见解可以评论一下。这里会发现t比较小的时候的波场占比会比较大，应该可以再做某些调整。

W2的做法似乎也是差不多这个想法，把波场当成分布函数，但是波场这里不一定都是正的，是不是得有某些操作处理一下，而且我也没法从波场做采样。另外，个人感觉这里不应该当成二维图像做sinkhorn，而应该直接对那个1000维的向量求误差，再把总共的 $5\times 70=350$个都加起来。

几种经典传统方法的结果对比（curvevelA版本）

省流

1. 如果有好的初值估计的话，无脑对逐通道使用W2，效果奇佳。
2. 这里的W2使用的是加上波场的最小值再做归一化，直接取绝对值效果不佳
3. 在没有好的初值的情况下，不管用什么损失函数，不管是否加tv正则化，都是寄（层数少的话用W2可能可以勉强做一做）。
4. 在异侧的情况不见得比在同侧的情况好求解。
5. 不推荐使用sinkhorn（不排除我参数没找好的情况）。
6. 我现在也不知道哪个算子简单了，难以描述，建议看下面的结果。个人感觉如果本来初值的表面就很接近真实解的话（也即面波本来的误差就非常小），打点在表面甚至更好做。尤其是逐通道的W2.
7. 不知道为什么积分loss失效了，完全跑不通（每次都是一开始相对误差就不下降，最后得到的图像也是完全不对的）（现在的结果我感觉是取绝对值的问题，换成$W^2$那种减去最小值就能降了，还是能得到形状的，效果倒是没有W2那么好，当然不排除是调参的原因）

震源和接收器都在左

无正则化，不添加噪声

1. 使用$L^2$ norm: 1. Adam优化器：1. 初始化为真解的模糊化，具体公式为

   $$\text{output}(x,y) = \frac{1}{\iint \frac{1}{v_{\text{true}}(u,v)} \cdot G(x-u,y-v) \, du \, dv}$$

其中$G$为高斯核函数：

$$G(x,y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2 + y^2}{2\sigma^2}}$$

对应的离散格式为

$$\text{output}[i,j] = \frac{1}{\sum_{m=-k}^{k} \sum_{n=-k}^{k} \frac{1}{v_{\text{true}}[i-m,j-n]} \cdot G[m,n]},\quad G[m,n] = \frac{1}{2\pi\sigma^2} e^{-\frac{m^2 + n^2}{2\sigma^2}}$$

此处$\sigma$越大，初始解就越模糊，离真实解就越远。取$\sigma=10$，初始相对误差为12.1%，迭代的结果如下（分别为结果，解的误差，波场的误差），相对误差为7.7%（注意，此时并不会出现所谓半收敛性的线性的现象），相对误差虽然一直在降，但是很明显整个图已经脱离数据集了。理论上应该还能降）

1. 取$\sigma=100$，初始相对误差为26.6%，迭代的结果如下，相对误差为26.7%

2. 初始化为常数值（这里我取的中间值3000），也会在迭代后期出现上面速度太小导致的警告，虽然强行运行也是能得到一个结果的

这时我们就能明显地看到半收敛性现象了（当上面的$\sigma$较大时，也会出现相对误差不降的现象）> deepwave的建议是一个完整波形的长度至少需要6个网格点，否则可能无法准确模拟波形。根据波长=速度/频率，速度太小时，波长也会相应的小，所包含的网格点也就少，所以会出现上面的问题。此时在物理上无法进行有效的模拟2. SGD：结果有点amazing，相对误差几乎不带降，也就是基本没做移动（这里我学习率设置为50），把学习率调到5000，loss开始降了，但相对误差几乎不变（应该是在局部极小点附近了）

1. lbfgs优化（原来torch也是有lbfgs优化器的）：我们继续重复上面的操作，结果也是很amazing，如果使用默认的使用 strong_wolfe进行线搜索，也会出现上面和上面SGD同样的情况。自己实现的lbfgs，使用原步长结果也是不太能动。

2. 使用W2 metric（以下用的Adam优化器）

   1. 使用sinkhorn：初始解依然设置为真解的模糊化，使用adam优化器，$\sigma=10$时出现了速度过小的情况，而且是在前几步就出现了这个问题，后面的结果也不收敛。

   

3. 把每个接收器和震源的pair对应的一维信号当作分布函数，直接求相应的分布函数，然后求逆，然后再直接计算误差

   1. 取$\sigma=10$，结果如下，相对误差为5.1%

   2. 取$\sigma=100$，结果如下，相对误差为22.1%

   3. 取初始值为常数值3000，结果如下

添加TV正则化

1. 定义：增加一个惩罚项：图像梯度的L2范数

$$TV(u)=\sum\_{i,j} \sqrt{(u[i+1, j] - u[i, j])^2 + (u[i, j+1] - u[i, j])^2}$$

该惩罚项能促进梯度域的稀疏性。鼓励大部分像素梯度为0（平滑区），允许少数像素梯度很大（边缘），即保边去噪。此时的loss_func我们定义为：

$$J(c)=MSE(d,u(c))+\lambda TV(c)$$

这里的d为观测到的波场，迭代过程中的速度场，u(c)为相应的模拟波长，$\lambda$为正则化项的系数（这里由于速度场的值本身比较大(1500~4500)，所以必须取一个比较小的值（比如0.001就太大了，相对误差会一直上升）。2. 使用L2 norm

1. Adam优化器：

2. 我们首先取$\lambda=1e-6$，$\sigma=10$，相对误差为7.3%，得到的图像如下。误差略微下降，并且图像看起来是会更接近数据集一点。

1. 接下来取$\lambda=1e-5$，$\sigma=10$，相对误差为8.4%，得到的图像如下。误差反而增加了，d但图像看起来是会更接近数据集一点。

3. 使用W2 norm

   1. 取$\lambda=1e-6$，使用sinkhorn: ，取$\sigma=10$，结果并没有得到改善：
   2. 取$\lambda=1e-5$
   3. 使用sinkhorn: ，取$\sigma=10$，结果并没有得到改善：

4. 逐通道使用W2norm

   1. 取$\sigma=10$，结果如下，相对误差为3.7%

   

不同正则化参数误差对比

mse：

Sinkhorn:

震源在左，接收器在右

无正则化，不添加噪声

1. 使用L2 norm
   1. Adam优化器：
2. 初始化为真解的模糊化，取$\sigma=10$，初始相对误差为12.1%，迭代的结果如下，相对误差为8.8%，比震源和接收器都在左侧高。本来以为应该是比都放在左侧简单的，但实验结果好像不是这样的，并且出现了半收敛性？

1. 取$\sigma=100$，还是寄

2. SGD：相对误差几乎不带降，也就是基本没做移动（这里我学习率设置为50），把学习率调到5000，loss开始降了，但相对误差几乎不变（应该是在局部极小点附近了）

3. lbfgs优化（原来torch也是有lbfgs优化器的）：我们继续重复上面的操作，如果使用默认的使用 strong_wolfe进行线搜索，也会出现上面和上面SGD同样的情况。自己实现的lbfgs，使用原步长结果也是不太能动。

4. 使用W2 metric（以下用的Adam优化器）

   1. 使用sinkhorn：初始解依然设置为真解的模糊化，使用adam优化器，$\sigma=10$时出现了速度过小的情况，而且是在前几步就出现了这个问题，后面的结果也不收敛。

取$\sigma=100$，结果寄：

1. 使用W2：（不如用MSE）
   1. 使用sinkhorn：（不管模糊化程度多少，这个都是寄）$\lambda=1e-5$得到的结果如下，相对误差大小为33.7%（寄）：
   

震源在地表，接收器在地底（source_depth=2, receiver_depth=68)

无正则化，不添加噪声

1. 使用MSE loss

   1. 初始化为真解的模糊化，取$\sigma=10$，初始相对误差为12.1%，迭代的结果如下，相对误差为8.7%，似乎也没有很简单

   

2. 使用W2 metric（以下用的Adam优化器）

   1. 使用sinkhorn：初始解依然设置为真解的模糊化，使用adam优化器，$\sigma=10$时出现了速度过小的情况，而且是在前几步就出现了这个问题，后面的结果也不收敛。
   2. 使用swd：1.$\sigma=10$，效果也还可以，相对误差大小为7.2%

添加TV正则化

1. 使用MSE loss

   1. 初始化为真解的模糊化，取$\sigma=10$，$\lambda=1e-6$，初始相对误差为12.1%，迭代的结果如下，相对误差为7.5%

   

2. 使用W2 metric（以下用的Adam优化器）

   1. 使用sinkhorn：初始解依然设置为真解的模糊化，使用adam优化器，$\sigma=10$时出现了速度过小的情况，而且是在前几步就出现了这个问题，后面的结果也不收敛。
   2. 使用swd：1.$\sigma=10$，效果也还可以，相对误差大小为7.2%

震源和接收器均在地表（source_depth=2, receiver_depth=2)

无正则化，不添加噪声

1. 使用MSE loss

   1. 初始化为真解的模糊化，取$\sigma=10$，初始相对误差为12.1%，迭代的结果如下，相对误差为7.0%

   

2. 使用W2 metric（以下用的Adam优化器）

   1. 使用sinkhorn：初始解依然设置为真解的模糊化，使用adam优化器，$\sigma=10$时出现了速度过小的情况，而且是在前几步就出现了这个问题，后面的结果也不收敛。

3. 逐通道使用W2norm

   1. 初始化为真解的模糊化，取$\sigma=10$，初始相对误差为12.1%，迭代的结果如下，相对误差为4.2%

   

添加TV正则化

1. 使用MSE loss

   1. 初始化为真解的模糊化，取$\sigma=10$，$\lambda=1e-6$，初始相对误差为12.1%，迭代的结果如下，相对误差为6.5%

   

2. 使用W2 metric（以下用的Adam优化器）

   1. 使用sinkhorn：初始解依然设置为真解的模糊化，使用adam优化器，$\sigma=10$时出现了速度过小的情况，而且是在前几步就出现了这个问题，后面的结果也不收敛。

3. 逐通道使用W2norm

   1. 初始化为真解的模糊化，取$\sigma=10$，$\lambda=1e-6$，初始相对误差为12.1%，迭代的结果如下，相对误差为3.2%

   

包络反演

地震数据的包络

1. 希尔伯特变换：

   $$H\left\{ f\left( t \right) \right\}=\frac{1}{\pi}\mathrm{P.V.} \int_{-\infty}^{+\infty}\frac{f\left( \tau \right)}{t-\tau}\text{d}\tau = f (\frac{1}{\pi t})$$

其中$\mathrm{P.V.}$代表柯西主值，定义如下（假设c处是奇点）：

$$\mathrm{P.V.} \int_a^b f(x)\, dx = \lim_{\varepsilon \to 0^+} \left[ \int_a^{c-\varepsilon} f(x) \, dx + \int_{c+\varepsilon}^b f(x) \, dx \right]$$

希尔伯特变换相当于把信号所有频率分量相位推迟90度，对于$f(t)=e^{i(wt+\phi)},H\{f(t)\}=e^{i(wt+\phi-\frac{\pi}{2})}$

2. 对于实信号$f(t)$解析信号 $\hat{f}\left( t \right)=f\left( t \right)+\text{i}H\left\{ f\left( t \right) \right\}$ 我们有如下定义：
3. 瞬时振幅（也称为信号的包络）：$E\text{}\left( t \right)=\sqrt{{{\left[ f\left( t \right) \right]}^{2}}+{{\left[ H\{f\left( t \right)\} \right]}^{2}}}$
4. 瞬时相位：$\Phi \left( t \right)=\arctan\frac{H\left\{ f\left( t \right) \right\}}{f\left( t \right)}$
5. $\hat{f}\left( t \right)=E\left( t \right){{\text{e}}^{^{\text{i}\mathit{\Phi}(t)}}}$