基于Python实现Marmousi速度模型的2D移动平均滤波方法咨询

一、2D移动平均滤波器的实现 #

2D移动平均滤波的核心是对每个网格点，计算其邻域窗口内所有速度值的算术均值并替换原网格值。以下是两种实用的实现方式：

1. 基于SciPy的快速实现 #

利用scipy.ndimage.uniform_filter可以高效完成滤波，这是地震数据处理中最常用的方式，无需手动遍历窗口：

import numpy as np
from scipy.ndimage import uniform_filter

# 假设marmousi_vel是加载后的速度模型二维数组（shape: [nz, nx]）
# 设置窗口大小，格式为(纵向窗口数, 横向窗口数)
window_size = (7, 11)
# 执行滤波，mode='reflect'用于避免边缘失真
smoothed_vel = uniform_filter(marmousi_vel, size=window_size, mode='reflect')

2. 手动实现（原理演示用） #

如果需要自定义边界逻辑或特殊窗口形状，可以手动遍历网格计算均值：

import numpy as np

def moving_average_2d(vel_model, z_win, x_win):
    nz, nx = vel_model.shape
    z_half = z_win // 2
    x_half = x_win // 2
    smoothed = np.zeros_like(vel_model, dtype=np.float32)
    
    for i in range(nz):
        z_start = max(0, i - z_half)
        z_end = min(nz, i + z_half + 1)
        for j in range(nx):
            x_start = max(0, j - x_half)
            x_end = min(nx, j + x_half + 1)
            smoothed[i, j] = np.mean(vel_model[z_start:z_end, x_start:x_end])
    return smoothed

# 使用示例
smoothed_vel = moving_average_2d(marmousi_vel, z_win=7, x_win=11)

二、适用于Marmousi模型的窗口大小选择 #

Marmousi模型的典型网格间距为25m（横向x）和25m（纵向z），窗口大小需平衡噪声抑制与地质构造保留：

• 纵向（z方向）: 选择5_{9个网格（125}225m）。浅部地层厚度较小（100_{200m），窗口取5}7个网格；深部地层厚度增大，可放宽至7~9个网格，避免过度模糊深部的速度渐变层。
• 横向（x方向）: 选择9_{17个网格（225}425m）。Marmousi模型的横向地质构造（如盐丘边缘、断层）尺度多在500m以上，该范围窗口既能滤除小尺度高频噪声，又不会模糊关键构造的边界。
• 特殊场景: 若高频噪声主要集中在浅部，可采用自适应窗口：浅部用小窗口，深部用大窗口，通过条件判断动态调整窗口大小。

注意：窗口大小建议取奇数，确保滤波后网格点与原位置对齐，避免偏移误差。

内容的提问来源于stack exchange，提问作者Anthony rebello