f-x反褶积

f-x反褶积（链接）
Canales, 1984, Random noise reduction, 54.th. Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, pp. 525-527

简介

$f-x$ 反卷积是一种地震数据去噪方法，其基本假设是具有线性同向轴的地震数据在$f-x$域中的每个频率切片都可以进行自回归表示。$f-x$反卷积方法是很多算法的基础，因此理论推导写的比较详细。之所以称为反卷积，可以理解为求卷积系数，也可以理解为求降噪数据。（一般反卷积是指求数据，而不是系数）。

理论

一般理论

平面波满足如下微分方程：
$$p(x,t)\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}=0$$
假设$p(x,t)=p$为常数，则有
$$p\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}=0$$
该方程具有通解：
$$u(x,t)=v(t_0+t-px)$$
上式只适用于描述含有单独斜率的线性同相轴，对于含有多线性同向轴的二维地震数据，在$t-x$域可以标示成如下形式：
$$u(t,x)=\sum _{i=1}^{n}v(t_i+t-p_{i}x)$$
其中$n$为同向轴个数，$t_i$为初至，$p_i$为斜率,$v$为地震子波。在两侧分别对时间进行傅里叶变换，根据傅里叶变换的性质：
$$F(f(x-a))=e^{-ia\omega }F(f(x))$$
得到（注意，上标$i$表示虚数单位，下标$i$表示索引）：
$$\hat{u}(\omega ,x)=\sum _{i=1}^n\hat{v}(\omega )e^{i(t_i-p_{i}x)\omega }$$

简化情况

对于简单的情况，同向轴只有一条，并且经过原点，上式可以简化为：
$$\hat{u}(\omega ,x)=\hat{v}(\omega )e^{-ipx\omega }$$
在$x$方向上进行离散化，$x=k\Delta x$，得到：
$$\hat{u}(\omega ,k\Delta x)=\hat{v}(\omega )e^{-ipk\Delta x\omega }$$
令$U_k=\hat{u}(\omega ,k\Delta x)$，$K=e^{-ip\Delta x\omega }$则：
$$U_k=U_{k-1}K$$
则频率切片的自回归形式为：
$$U_k-K{U}_{k-1}=0$$
将所有关系式写到一起，形成一个线性方程组，然后通过最小二乘方法可以得到对所有$i$最优的$K$（注意$p$最开始是未知的，$K$也是未知的，所以需要通过最小二乘来求）。得到$K$之后，为了实现去噪，虽然自回归滤波器的长度为2，但我们可以利用更多的相邻到来实现叠加去噪，即：
$$U_i=\frac{1}{m}\sum _{k=i-m}^{k=i-1}{K}^{i-k}{U}_k$$

反推一般情形

此部分推导相对复杂一些，之前一直没整明白，这次写文档算是整明白了。假设数据中含有多个线性同向轴（并且都通过原点$t_i=0$），则频率域的关系式为：
$$\hat{u}(\omega ,k\Delta x)=\hat{v}(\omega )\sum _{i=1}^n{e}^{-i{p}_{i}k\Delta x\omega }$$
令$U_k=\hat{u}(\omega ,k\Delta x)$，$K_i=e^{-i{p}_i\Delta x\omega }$，$V=\hat{v}(\omega )$则：
$$U_k=V(\sum _{i=1}^n{K}_i^k)$$
需要证明在$U$存在自回归性质，即存在$a_0,a_1,\cdots ,a_{n-1}$，使得下面方程对任意$j$成立：
$$\sum _{k=j}^{j+n-1}{a}_{k-j}{U}_k=0$$
将$U_k$带入可以得到：
$$\sum _{k=j}^{j+n-1}{a}_{k-j}\sum _{i=1}^n{K}_i^k=0$$
由于后面的求和下标与$k$无关，所以可以把$a_{k-j}$放进去，得到：
$$\sum _{k=j}^{j+n-1}\sum _{i=1}^n{a}_{k-j}{K}_i^k=0$$
两组求和下标之间也没有关系，因此求和可以交换顺序：
$$\sum _{i=1}^n\sum _{k=j}^{j+n-1}{a}_{k-j}{K}_i^k=0$$
到这是关键的一步，我们假设对于第一个求和而言，每一项都为零，则求和也为零。得到：
$$\sum _{k=j}^{j+n-1}{a}_{k-j}{K}_i^k=0$$
两边再除以$K_i^j$得到：
$$\sum _{k=0}^{n-1}{a}_k{K}_i^k=0$$
我们发现对于不同的$i$，所需满足的方程是一样的。并且这个方程是有解的：可以得到一个具有$n$个未知数的方程，对于所有的$i$，共有$n$个方程，因此可以解出来$a_k$。
此时可知$U$是线性自回归的，并且回归系数与$p_i$有关。解出来回归系数之后，可以通过某一道与相邻道之间的线性回归关系来进行去噪。

同$n=1$的情况类似，在含有噪声的情况下，我们通常会选择更长的自回归滤波器，以达到更好的去噪效果。

$n=2$时的特例

$n=2$时的方程组为：
$$\begin{array}{rl}{a}_0+a_1{K}_1+K_1^2=& 0\\ a_0+a_1{K}_2+K_2^2=& 0\end{array}$$
可以得到：
$$\begin{array}{rl}{a}_0& =K_1{K}_2\\ a_1& =-K_1-K_2\end{array}$$

初至、通解

对于加入初至的情况，直接令$K_j=e^{i(t_j-p_j\Delta x)\omega }$即可。

一般情况也存在通解，只不过表达式比较复杂，并且并没有在算法中使用，所以就不推导了。

算法

基本算法是：

1. 将地震数据变换到$f-x$域；
2. 提取某一频率切片数据进行自回归表示（需要设定长度），求取自回归系数；
3. 利用求得的自回归系数对频率切片进行自回归；
4. 对每一频率切片重复2-3；
5. 将数据变回到$t-x$域。

如果数据非线性很强，需要分窗口处理。

代码

代码源自SeismicLab工具包，fx_decon程序，只摘录其核心部分。

  aux_in  = DATA_FX(k,:)'; %提取一个频率切片
  [aux_out_f,aux_out_b] = ar_modeling(aux_in,lf,mu); %自回归函数，lf为回归滤波器长度，mu应该是最二乘添加的约束项权重。

ar_modeling的核心代码为：

% backward ar-modeling 此部分为向后ar回归，同时可以再次进行向前ar回归，将两次结果作平均作为最后结果。

   y  = x(1:nx-lf,1);    %输出
   C  = x(2:nx-lf+1,1);  
   R  = x(nx-lf+1:nx,1);
   M = hankel(C,R);      %输入，M*ab=y结构如下
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   x_2          x_3 ...         x_lf+1      a_1     x_1
   x_3          x_4 ...         x_lf+2  *   a_2  =  x_2
   ...  
   x_nx-lf+1    x_nx-lf+2 ...   x_nx        a_lf    x_nx-lf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   B = M'*M;  beta = B(1,1)*mu/100;     %两边同时乘以M'，M'*M*ab=M'*y
   ab = (B + beta*eye(lf))\M'*y;        %ab=(M'M+mu*I)^-1*M'y 加入权重
   yb = M*ab;                           %重新做滤波得到去噪后数据

总结

在分析多线性同向轴的自回归模型是所用理论虽不困难，但也需要一些技巧，通过分别将含有$K_j^i$（对任意$i$）的项组合然后分别另其为零即可推导出自回归性质。这篇短文给出了所有的推导细节，作为之后的理论基础，在Spitz $f-x$插值（据作者所知，去噪称为$f-x$反褶（卷）积，插值称为$f-x$预测滤波，因为包含一步预测）、MSSA、降秩方法中都用到类似的推导。