FK滤波变换

f-k变换实质上是一种二维傅里叶变换。
在时间域上，对一道地震信号做傅里叶变换，可以得到在时间上不同频度(称为频率)的 波动组分的振幅和相位信息。同样地，在空间上，我们也可以对多道地震信号做类似傅里叶变 换的数值变换，得到在空间上不同频度(称为波数)的波动组分的振幅和相位信息。进行这两种变换后，便可以分析一个炮集记录在频率-波数域的能量分布情况。

• 面波相速度和频率、波数的关系
  

def FKDisp(dat, dt, dx, vmin, dv, vmax, fmin, df, fmax):
    Nt, Ntrace = dat.shape  # 读取单炮记录的时间采样数与道数
    
    # 振幅归一化
    for i in range(Ntrace):
        dat[:, i] /= np.max(np.abs(dat[:, i]))
    
    nk = 4 * (2 ** int(np.ceil(np.log2(Ntrace))))  # 计算波数维度
    nf = 4 * (2 ** int(np.ceil(np.log2(Nt))))     # 计算频率维度
    
    f = np.arange(nf) / nf / dt  # 频率序列
    k = np.arange(nk) / nk / dx  # 波数序列
    wave_num = k[::-1]          # 倒序排列
    
    S = np.fft.fft2(dat, s=[nf, nk])  # 二维傅里叶变换
    f_k = S.T  # 转置
    
    v = np.arange(vmin, vmax + dv, dv)  # 速度显示区间
    ff = np.arange(fmin, fmax + df, df)  # 频率显示区间
    
    nff = len(ff)  # 频率序列长度
    nv = len(v)    # 速度序列长度
    
    disper = np.zeros((nv, nff))  # 定义频散矩阵
    
    # 插值到频率-速度显示区间
    for i in range(nff):
        interp_func = interp2d(f, wave_num, abs(f_k), kind='linear')
        disper[:, i] = interp_func(ff[i], ff[i] / v).flatten()  # 使用flatten()将结果转换为一维数组  v=f/k  k=f/v
    
    disper[np.isnan(disper)] = 0  # 找出nan值并赋值为0
    
    mn = np.min(disper, axis=0)
    mx = np.max(disper, axis=0)
    
    for i in range(len(mn)):
        disper[:, i] = (disper[:, i] - mn[i]) / (mx[i] - mn[i])  # 归一化
    
    EN = np.abs(disper)
    
    EN = np.flipud(EN)
    return EN, ff, v
# 加载数据

FK滤波

daspy fk
FK滤波方法是使用二维快速傅里叶变换将数据转换至频率-波数（FK）域，在FK域乘上滤波窗后再反变换回时间-距离域。若该地区交通信号的视波速小于1000m/s，地震波的视波速大于1000m/s，则可以以 vmin=(800, 1200) 来滤除交通信号，第一个输出为地震波，第二个输出为交通信号， 800 和 1200 为尖灭的边界。