利用f-k滤波分离上行波和下行波的代码

下面是一个基于Python语言的f-k滤波分离上行波和下行波的代码示例：

import numpy as np
from scipy import fft, ifft

def fk_filter(seis_data, dx, dy, f_up, f_down):
    """
    利用f-k滤波分离上行波和下行波
    :param seis_data: 地震数据，二维数组
    :param dx: X方向采样间隔
    :param dy: Y方向采样间隔
    :param f_up: 上行波频率范围，单位Hz
    :param f_down: 下行波频率范围，单位Hz
    :return: 上行波数据，下行波数据
    """
    # 快速傅里叶变换，得到频率域数据
    freq_data = fft.fft2(seis_data)

    # 获取数据长度和宽度，以及X和Y方向上的采样点数
    nx, ny = seis_data.shape[1], seis_data.shape[0]
    nf, nky = nx//2+1, ny//2+1

    # 创建频率-波数网格
    kx = np.arange(0, nx//2+1) * 2 * np.pi / (nx*dx)
    ky = np.arange(0, ny//2+1) * 2 * np.pi / (ny*dy)
    kx = np.concatenate((kx, kx[1:-1][::-1]))
    ky = np.concatenate((ky, ky[1:-1][::-1]))
    kx, ky = np.meshgrid(kx, ky)

    # 计算频率-波数网格上的频率
    freq = np.sqrt(kx**2 + ky**2)

    # 根据上下行波的频率范围选择相应的频率-波数域数据
    up_mask = np.logical_and(freq >= f_up-dx/2, freq <= f_up+dx/2)
    down_mask = np.logical_and(freq >= f_down-dx/2, freq <= f_down+dx/2)
    up_data = np.zeros((ny, nx), dtype=np.complex)
    down_data = np.zeros((ny, nx), dtype=np.complex)
    up_data[up_mask] = freq_data[up_mask]
    down_data[down_mask] = freq_data[down_mask]

    # 对选择的频率-波数域数据进行反f-k变换，得到上行波和下行波的时域数据
    up_wave = np.real(ifft.ifft2(up_data))
    down_wave = np.real(ifft.ifft2(down_data))

    # 对时域数据进行滤波，去除噪声和干扰信号，这里使用一个简单的高通滤波器
    up_wave = np.where(up_wave > 0, up_wave, 0)
    down_wave = np.where(down_wave > 0, down_wave, 0)

    return up_wave, down_wave

使用示例：

import matplotlib.pyplot as plt

# 生成一个模拟地震数据
t = np.linspace(0, 1, 100)
x = np.linspace(0, 1, 50)
y = np.linspace(0, 1, 50)
xx, yy = np.meshgrid(x, y)
seis_data = np.sin(2*np.pi*10*t) * np.exp(-((xx-0.5)**2 + (yy-0.5)**2)/0.1)

# 分离上行波和下行波
up_wave, down_wave = fk_filter(seis_data, 0.02, 0.02, 10, 20)

# 绘制结果
fig, axs = plt.subplots(2, 2, figsize=(8, 8))
axs[0][0].imshow(seis_data, cmap='seismic', aspect='auto')
axs[0][0].set_title('Original data')
axs[0][1].imshow(up_wave, cmap='seismic', aspect='auto')
axs[0][1].set_title('Up wave')
axs[1][0].imshow(down_wave, cmap='seismic', aspect='auto')
axs[1][0].set_title('Down wave')
axs[1][1].imshow(up_wave+down_wave, cmap='seismic', aspect='auto')
axs[1][1].set_title('Up wave + Down wave')
plt.tight_layout()
plt.show()

该示例中，我们生成了一个模拟地震数据，并利用f-k滤波将上行波和下行波分离出来。最终绘制了原始数据、上行波、下行波和上行波加下行波的图像。