从仿真动画到代码：我是如何用Python复现超声波束合成中的Grating Lobe和Side Lobe的

从仿真动画到代码：Python复现超声成像中的Grating Lobe与Side Lobe现象

超声成像技术作为现代医学诊断的重要工具，其核心原理是通过声波的发射与接收构建人体内部结构的图像。在这个过程中，波束合成的质量直接影响成像的清晰度和准确性。本文将带你用Python一步步实现超声成像中的两个关键现象——Grating Lobe（栅瓣）和Side Lobe（旁瓣）的仿真，通过代码直观展示这些现象的产生原理及其对图像质量的影响。

1. 超声成像基础与环境搭建

在开始编码之前，我们需要理解几个基本概念。超声探头由多个阵元(element)组成，通过控制各阵元的发射延时形成聚焦声束。接收时，则通过动态聚焦技术将各阵元接收到的信号合成。这一过程中，阵元间距、发射频率等参数会直接影响波束特性。

首先配置Python环境，我们将使用以下主要库：

import numpy as np import matplotlib.pyplot as plt from scipy.signal import gausspulse from scipy.ndimage import gaussian_filter

• numpy：用于数值计算和数组操作
• matplotlib：用于数据可视化
• scipy.signal：提供信号处理工具
• scipy.ndimage：用于图像处理

提示：建议使用Python 3.8+版本，并创建独立的虚拟环境以避免依赖冲突。

2. 超声场建模与波束形成

2.1 探头参数设置

我们先定义一个64阵元的线阵探头，设置基本参数：

# 探头参数 num_elements = 64 # 阵元数量 pitch = 0.3e-3 # 阵元间距(米) center_freq = 3e6 # 中心频率3MHz sampling_freq = 50e6 # 采样频率50MHz c = 1540 # 声速(m/s)

2.2 发射波束形成

实现偏转45度的发射波束形成：

def create_transmit_delays(num_elements, pitch, angle, focus_depth, c): """计算各阵元的发射延时""" positions = np.arange(num_elements) * pitch - (num_elements-1)*pitch/2 delays = np.sqrt(positions**2 + focus_depth**2 - 2*positions*focus_depth*np.sin(angle))/c return delays - delays.min() # 45度偏转，聚焦深度30mm angle = np.radians(45) focus_depth = 30e-3 tx_delays = create_transmit_delays(num_elements, pitch, angle, focus_depth, c)

2.3 声场计算

计算空间各点的声压分布：

def calculate_pressure_field(tx_delays, num_elements, pitch, freq, c, grid_size=(100,100)): """计算声场压力分布""" x = np.linspace(-20e-3, 20e-3, grid_size[0]) z = np.linspace(0, 40e-3, grid_size[1]) X, Z = np.meshgrid(x, z) pressure = np.zeros_like(X) for i in range(num_elements): elem_pos = i * pitch - (num_elements-1)*pitch/2 distance = np.sqrt((X - elem_pos)**2 + Z**2) arrival_time = distance / c - tx_delays[i] pressure += np.cos(2 * np.pi * freq * arrival_time) * (arrival_time >= 0) return x, z, pressure x, z, pressure = calculate_pressure_field(tx_delays, num_elements, pitch, center_freq, c)

3. Grating Lobe现象仿真与分析

3.1 可视化声场分布

将计算得到的声场可视化：

plt.figure(figsize=(10,8)) plt.imshow(20*np.log10(np.abs(pressure)+1e-6), extent=[x.min()*1e3, x.max()*1e3, z.max()*1e3, z.min()*1e3], cmap='gray', vmin=-40, vmax=0) plt.colorbar(label='声压(dB)') plt.xlabel('横向位置(mm)') plt.ylabel('深度(mm)') plt.title('45度偏转发射声场分布') plt.show()

3.2 Grating Lobe形成机制

Grating Lobe是由于阵元间距过大导致的伪影。当阵元间距d与波长λ满足以下关系时，会出现明显的Grating Lobe：

d/λ > 1/(1 + |sinθ|)

其中θ为偏转角度。我们可以通过调整参数观察Grating Lobe的变化：

3.3 参数影响实验

通过代码系统研究各参数影响：

def analyze_grating_lobe_conditions(): """分析不同参数下Grating Lobe表现""" conditions = [ {'pitch':0.2e-3, 'freq':5e6, 'angle':30}, {'pitch':0.3e-3, 'freq':3e6, 'angle':45} ] fig, axes = plt.subplots(1, 2, figsize=(15,6)) for ax, params in zip(axes, conditions): tx_delays = create_transmit_delays(num_elements, params['pitch'], np.radians(params['angle']), 30e-3, c) _, _, pressure = calculate_pressure_field(tx_delays, num_elements, params['pitch'], params['freq'], c) ax.imshow(20*np.log10(np.abs(pressure)+1e-6), cmap='gray', vmin=-40, vmax=0) ax.set_title(f"间距{params['pitch']*1e3:.1f}mm, {params['freq']/1e6}MHz, {params['angle']}°") plt.tight_layout() plt.show() analyze_grating_lobe_conditions()

4. Side Lobe现象与抑制技术

4.1 Side Lobe的产生

Side Lobe是主瓣周围出现的次级能量峰，主要由孔径有限导致的衍射引起。我们可以通过变迹(Apodization)技术来抑制Side Lobe。

实现汉宁窗变迹：

def apply_apodization(pressure_field, window_type='hanning'): """应用变迹窗函数""" if window_type == 'hanning': window = np.hanning(num_elements) elif window_type == 'uniform': window = np.ones(num_elements) apodized_pressure = np.zeros_like(pressure_field) for i in range(num_elements): apodized_pressure += pressure_field[:,:,i] * window[i] return apodized_pressure

4.2 变迹效果对比

比较不同变迹方式的效果：

# 计算无变迹声场 _, _, pressure_uniform = calculate_pressure_field(tx_delays, num_elements, pitch, center_freq, c) # 计算汉宁窗变迹声场 pressure_hanning = apply_apodization(pressure_uniform, 'hanning') # 绘制对比 plt.figure(figsize=(12,5)) plt.subplot(121) plt.imshow(20*np.log10(np.abs(pressure_uniform)+1e-6), cmap='gray', vmin=-60, vmax=0) plt.title('无变迹') plt.subplot(122) plt.imshow(20*np.log10(np.abs(pressure_hanning)+1e-6), cmap='gray', vmin=-60, vmax=0) plt.title('汉宁窗变迹') plt.show()

变迹技术对成像质量的影响：

1. 优点：

   • 显著降低Side Lobe幅度
   • 提高图像对比度
   • 减少伪影干扰

2. 代价：

   • 主瓣略微展宽
   • 轻微降低分辨率

5. 完整成像仿真流程

5.1 靶点模型建立

创建一个包含多个理想靶点的测试模型：

def create_target_model(grid_size=(100,100)): """创建靶点测试模型""" model = np.zeros(grid_size) # 添加靶点 model[30, 50] = 1 # 中心靶点 model[30, 70] = 1 # 右侧靶点 model[50, 50] = 0.8 # 下方靶点 return gaussian_filter(model, sigma=1) target_model = create_target_model()

5.2 完整成像仿真

实现从发射到接收的完整成像流程：

def simulate_imaging(target_model, tx_angle, tx_focus, rx_focus, apodization='hanning'): """完整成像仿真""" # 发射波束形成 tx_delays = create_transmit_delays(num_elements, pitch, np.radians(tx_angle), tx_focus, c) # 声场传播 _, _, tx_pressure = calculate_pressure_field(tx_delays, num_elements, pitch, center_freq, c) # 接收信号模拟 rx_signals = np.zeros((num_elements, *target_model.shape)) for i in range(num_elements): elem_pos = i * pitch - (num_elements-1)*pitch/2 # 简化的接收信号模型 rx_signals[i] = tx_pressure * target_model * np.exp(-0.02 * distance_from_element(elem_pos)) # 接收波束形成 rx_delays = create_transmit_delays(num_elements, pitch, np.radians(tx_angle), rx_focus, c) if apodization == 'hanning': weights = np.hanning(num_elements) else: weights = np.ones(num_elements) beamformed = np.zeros(target_model.shape) for i in range(num_elements): beamformed += rx_signals[i] * weights[i] * compensate_delay(rx_delays[i]) return beamformed

5.3 结果分析与优化

对比不同参数下的成像质量：

# 不同参数组合成像 params_combinations = [ {'tx_angle':0, 'tx_focus':30e-3, 'rx_focus':30e-3, 'apodization':'uniform'}, {'tx_angle':45, 'tx_focus':30e-3, 'rx_focus':30e-3, 'apodization':'uniform'}, {'tx_angle':45, 'tx_focus':30e-3, 'rx_focus':30e-3, 'apodization':'hanning'} ] results = [] for params in params_combinations: img = simulate_imaging(target_model, **params) results.append(20*np.log10(np.abs(img)+1e-6)) # 可视化比较 fig, axes = plt.subplots(1, 3, figsize=(18,6)) titles = ['0度无变迹', '45度无变迹', '45度汉宁窗'] for ax, img, title in zip(axes, results, titles): ax.imshow(img, cmap='gray', vmin=-40, vmax=0) ax.set_title(title) plt.show()

在实际项目中，我们发现当使用3MHz中心频率、0.3mm阵元间距时，45度偏转会出现明显的Grating Lobe伪影。通过调整阵元间距至0.2mm或提高频率至5MHz，可以显著改善这一问题。同时，汉宁窗变迹能有效抑制Side Lobe，虽然会使主瓣略微展宽，但在大多数临床应用中是可以接受的折衷方案。