SSA 奇异普分析

原创已于 2022-11-10 14:34:50 修改 · 5.7k 阅读

16 ·

CC 4.0 BY-SA版权

文章标签：

#线性代数 #矩阵

于 2019-04-26 10:20:01 首次发布

Geophysics 专栏收录该内容

3 篇文章

订阅专栏

本文探讨了基于奇异谱分析(SSA)的地震数据降噪与插值技术，详细介绍了SSA方法的原理及其与f-x反褶积的联系。通过对地震数据进行傅立叶变换和矩阵低秩近似，SSA能有效抑制噪声，恢复信号。算法步骤包括傅立叶变换、低秩近似计算、系数平均及反变换。

Trickett S . F-xy Cadzow noise suppression[J]. Seg Technical Program Expanded Abstracts, 2008, 27(1):2586.
Trickett S R . Fx eigenimage noise suppression[J]. Seg Technical Program Expanded Abstracts, 2002, 21(1):2478.

简介

基于奇异普分析的地震噪声压制和插值方法也属于一大类方法，还有别的名字，如SSA，MSSA（多通道，三维），Cadzow方法（最原始的文献），本文将采样SSA，即奇异普分析这个名字。SSA的推导过程与 $f - x$ 反褶积类似，不同的地方在于将线性同相轴下的自回归性质用矩阵的低秩性质来表示，这里面涉及到一个向量到矩阵的变换，是通过Hankel矩阵实现的。另外在推导 $f - x$ 自回归滤波器时，我们也提到其中有个构造出来范德蒙的矩阵是低秩的。还有类外一类低秩方法，是通过块变换形成的矩阵的低秩性来做的，并不是通过线性同相轴推导的，因此不介绍的。

方法

之前我们在 $f - x$ 反褶积中推导过如下公式：
$u^(ω,kΔx)=v^(ω)∑j=1pe−ipjkΔxω\hat u(\omega,k\Delta x)=\hat v(\omega)\sum_{j=1}^pe^{-ip_jk\Delta x \omega}$
令 $tk=u^(ω,kΔx)t_k=\hat u(\omega,k\Delta x)$ ，构造如下矩阵：
$A=(t1t2t3⋯tN−n+1t2t3t4⋯tN−n+2t3t4t5⋯tN−n+3⋯tntn+1tn+2⋯tN)A=\left(\begin{matrix}t_1 & t_2 & t_3 &\cdots & t_{N-n+1}\\ t_2 & t_3 & t_4 &\cdots & t_{N-n+2}\\ t_3 & t_4 & t_5 & \cdots & t_{N-n+3} \\ \cdots \\ t_n & t_{n+1} & t_{n+2} & \cdots & t_{N} \\ \end{matrix}\right)$
此时， $A$ 的秩至多为 $p$ （Stephenson，1988）。一般取 $n = N /2$ 使得 $A$ 为方阵。

算法

对数据进行傅立叶变换形成矩阵 $A$
计算 $A$ 的低秩近似
在 $A$ 的反对角上进行平均得到得到新的傅立叶变换系数
傅立叶反变换得到去噪后的数据

 for k = ilow:ihigh; # 处理不同频率
 tmp  = DATA_FX_tmp(k,:)';
 for ic = 1:Ncol;
  for ir = 1:Nrow;
   M(ir,ic) = tmp(ir+ic-1); #形成Hankel阵
  end;
 end
 [U,S,V] = svd(M); #SVD分解
 SS(k,:) = diag(S); 
 for p=1:min([Nrow,Ncol]);
     if p~=P
         S(p,p)=0;
     end
 end;
 SSf(k,:) = diag(S);  # 阈值处理
 Mout = U*S*V';
Count = zeros(ntraces,1);
tmp = zeros(ntraces,1);
 for ic = 1:Ncol;
  for ir = 1:Nrow;
   Count(ir+ic-1,1) = Count(ir+ic-1,1)+1; # 反对角求和
   tmp(ir+ic-1,1)  = tmp(ir+ic-1,1) + Mout(ir,ic);
  end;
 end
 tmp = tmp./Count;
  DATA_FX_f(k,:) = tmp';
 end;

证明

若要证明 $A$ 的秩为 $p$ ，需证明任意一行均可以用其它任意 $p$ 行线性表示，取一个例子，假设第 $p + 1$ 行可以用前 $p$ 行线性表示，即：
$(tp+1tp+2tp+3⋯tN−n+p+1)=(a1a2a3⋯ap)(t1t2t3⋯tN−n+1t2t3t4⋯tN−n+2t3t4t5⋯tN−n+3⋯tptp+1tp+2⋯tN)\left(\begin{matrix} t_{p+1} & t_{p+2}& t_{p+3}&\cdots & t_{N-n+p+1}\end{matrix}\right) \\=\left(\begin{matrix} a_{1} & a_{2}& a_{3}&\cdots & a_{p}\end{matrix}\right) \left(\begin{matrix}t_1 & t_2 & t_3 &\cdots & t_{N-n+1}\\ t_2 & t_3 & t_4 &\cdots & t_{N-n+2}\\ t_3 & t_4 & t_5 & \cdots & t_{N-n+3} \\ \cdots \\ t_p & t_{p+1} & t_{p+2} & \cdots & t_{N} \\ \end{matrix}\right)$
将每一列分解出来，得到：
$tp+i=(a1a2a3⋯ap)(titi+1⋯ti+p−1)t_{p+i}=\left(\begin{matrix} a_{1} & a_{2}& a_{3}&\cdots & a_{p}\end{matrix}\right) \left(\begin{matrix} t_{i} \\ t_{i+1} \\ \cdots \\ t_{i+p-1} \\ \end{matrix}\right)$
这是一个回归关系，根据在 $f - x$ 论文中的推导，上式是成立的，证明完毕。
该证明也可以通过直接把线性组合的形式写出来。