坡面回归插值

简介

Spitz插值是最经典的插值方法，它利用了上一篇博客中提到的具有线性同相轴的地震数据可以在频率–空间域进行线性回归的理论。（回顾：线性回归滤波器的长度等于同相轴（不同斜率）的数目加一，实际应用中会采用更长的自回归滤波器，可以使算法更稳定。自回归的好处在于其只用到了地震数据到本身，而不需要知道斜率信息）

Spitz插值的思路是，首先利用低频无假频部分计算出自回归系数（回归），根据道加密后相邻道之间的相移为原来一半，可以推导出加密后数据的自回归系数，然后在利用自回归关系从已知数据推导出未知数据（预测）。回顾和预测都通过最小二乘实现，因此算法比较稳定。缺点是：1，只能用于线性同相轴；2，只能用于规则网格。

理论

以下符号与Spitz论文中一致。

包含$L$个同相轴和$N$道的地震数据在频率域可以表示成：

$g_k(f)=\sum_{j=1}^La_j(f)z_j^{k-1}(f),k=1,...,N$

其中$g_k(f)$表示第$k$道数据，$a_j(f)$为第$j$个同相轴对应的子波频谱（第一道数据），$z_j(f)=exp(2\pi i f p_j)$对应于第$j$个同相轴相邻道之间的相移$p_j$（对应于斜率）。根据上一篇博客介绍的Canales方法，向前–向后单步预测滤波可以写成如下形式：
$g_k(f)=\sum_{j=1}^L P_j(f)g_{k-j}(f),k=L+1,...,N$
$g_k^*(f)=\sum_{j=1}^L P_j(f)g_{k+j}^*(f),k=1,...,N-L$
其中$P_j(f)$为预测滤波器的系数。（向前–向后滤波的共轭关系可以从附录中得到）

Spitz方法在原来的地震道中两道中间插入新的地震道，因此插值后的地震数据满足：
$g_k$
可以再一次写出向前–向后单步预测滤波形式：
$g$
$g$

上述方程为关于$g$的线性方程组，系数为$P$将已知值和未知值分别放到一块，整理得到：
$A(P$

如果$P$已知的话，以上表达式包含$2(2N-L-1)$个方程（向前向后滤波等式的个数）和$N-1$个未知数（文献中写成$2N-1$，不知道是写错了，还是包括了所有道），可以用最小平方方法求解（等价于一个从已知到未知的滤波器）：

$(A^+A)^{-1}A^+B$

其中$A^+$为$A$的复共轭转置。

虽然实际上$P$未知，但是可以用$P$来推导。注意输出（包括已知点）和输入剖面的相移关系为$z$（注：$z$，插值后相邻道之间的相移为原来的一半）。

由$P$与$z$的关系可以得到：

$P$

可以看到，只利用了无假频的低频部分。

以上推导由一些细节没有给出，待续。

附录：自回归的另外一种推导

Spitz用另外一种方法进行了Canales文章的推导。
附推导过程
对于地震数据$x^T(f)=[x_1(f),...,x_N(f)]$，假设$x(f)=g(f)+n(f)$。定义z变换$S_j^T(f)=[1,z_j(f),...,z_j^{N-1}(f)]$，其中，$z_j(f)=exp(2\pi ifp_j)$，则：
$g(f)=\sum_{j=1}^La_j(f)S_j(f)$
其中$a_j(f)$为第$j$个同向轴对应的子波的傅里叶变换。分开写为：
$\begin{aligned} g_1&=&\sum_{j=1}^La_j\cdot 1 \\ g_2&=&\sum_{j=1}^La_j\cdot z_j \\ &\cdots& \\ g_N&=&\sum_{j=1}^La_j\cdot z_j^{N-1} \end{aligned}$
写成矩阵的形式为：
$\left[ \begin{matrix} g_1 \\ g_2 \\ \vdots \\ g_n \end{matrix} \right]= \left[ \begin{matrix} 1 & 1& \cdots & 1 \\ z_1 & z_2& \cdots & z_L \\ \vdots \\ z_1^{N-1} & z_2^{N-1}& \cdots & z_L^{N-1} \end{matrix} \right] \times \left[ \begin{matrix} a_1 \\ a_2 \\ \vdots \\ a_n \end{matrix} \right] \tag{3}$
令其中的矩阵为$S$，即$S=[S_1, \cdots, S_L]$。将(3)记为矩阵形式：
$x=Sa+n$
形为$S$的矩阵称为范德蒙矩阵，$S_j$之间是线性无关的，则$S$的秩为$L$（$N>=L$）。$S$的任意一行都可以表示为其它$L$行的线性组合。例如，第$L+1$行可以表示为前$L$行的线性组合。写成如下形式：
$(z_1^L,\cdots,z_L^L)=(P_L,\cdots,P_1)\times \left[ \begin{matrix} 1 & 1& \cdots & 1 \\ z_1 & z_2& \cdots & z_L \\ \vdots \\ z_1^{L-1} & z_2^{L-1}& \cdots & z_L^{L-1} \end{matrix} \right] \tag{4}$
将上式右侧移到左侧，整理成关于$z_j$的方程，发现对任意$z_j$形式都一样，为：
$Z_L=z^L-P_1z^{L-1}-\cdots-P_L=0$
又因为$z_j$各不相同，上述方程有$L$个解（零点），所以$z_j$为上述方程的$K$个零点，即可以写成如下形式：
$Z_L=(z-z_1)(z-z_2)\cdots(z-z_L)$
对比上述两式可以得到
$\begin{aligned} P_1&=&z_1+\cdots+z_L \\ P_2&=&-(z_1z_2+\cdots+z_{L-1}z_L)\\ &\cdots& \\ P_L&=&(-1)^{L+1}z_1z_2\cdots z_L \end{aligned}$
上式即为我们之前没有推出来的通解。

代码解释：

核心代码：

PF = prediction_filter(x1,npf,pre1); #计算预测滤波器
  [y] = interpolate_freq(x2,PF,pre2); # 插值

自回归代码：

for j=1:ns-np  
    C(j,:)=VEC(j+np:-1:j);
 end  
 A = [C(:,2:np+1);conj(fliplr(C(:,1:np)))];  #向前向后（共轭）预测的输入
 B = [C(:,1);conj(C(:,np+1))];# 向前向后预测的输出
 R = A'*A; 
 g = A'*B;
 mu = (pre/100.)*trace(R)/np;
 PF =  (R+mu*eye(np))\g; #最小平方求解

插值代码：

TMPF1 = [fliplr(PF.'),-1];  #向前预测
 W1 = convmtx(TMPF1,ny-np);
 TMPF2 = conj(fliplr(TMPF1)); # 向后预测
 W2 = convmtx(TMPF2,ny-np);
 WT = [W1;W2];
 A = WT(:,2:2:ny); #未知数据
 B = -1*WT(:,1:2:ny)*x.'; #已知数据，这里可以看到，是用已知数据预测未知数据
 R= A'*A;
 g = A'*B;
 mu = (pre/100.)*trace(R)/(nx-1);
 y1 =  (R+mu*eye(nx-1))\g; #最小平方反演
 y = zeros(1,ny);
 y(1:2:ny)=x;
 y(2:2:ny)=y1.';