利用拉东变换分离上下行波的代码

下面是一个基于Python语言的拉东变换分离上行波和下行波的代码示例：

import numpy as np
from scipy import fft, ifft

def rda_filter(seis_data, dx, dy, f_up, f_down):
    """
    利用拉东变换分离上行波和下行波
    :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)
    alpha = np.sqrt(kx**2 + ky**2)

    # 计算上下行波的拉东变换系数
    up_alpha = np.where(np.logical_and(alpha < 2*np.pi*f_up, alpha > 0), alpha, 0)
    down_alpha = np.where(np.logical_and(alpha > 2*np.pi*f_down, alpha <= np.max(alpha)), alpha, 0)

    # 将拉东变换系数矩阵转换为拉东变换矩阵
    up_matrix = np.exp(-1j * up_alpha * dx)
    down_matrix = np.exp(-1j * down_alpha * dx)

    # 对频率域数据进行拉东变换
    up_data = np.real(ifft.ifft2(up_matrix * freq_data))
    down_data = np.real(ifft.ifft2(down_matrix * freq_data))

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

    return up_data, down_data

使用示例：

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 = rda_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()

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