基于延迟叠加算法的超声波束聚焦合成:揭秘DAS技术‌

Delay-and-SumDelay-and-Sum

基于延迟叠加算法的超声波束聚焦合成:揭秘DAS技术

本文将从超声成像的基本原理出发,系统介绍延迟叠加(Delay-and-Sum,简称 DAS)算法在超声波束形成(Beamforming)中的应用。文章包含数学推导、示意图与 Python 代码示例,帮助你直观理解 DAS 技术及其实现。

目录

  1. 引言
  2. 超声成像与束形成基础
  3. 延迟叠加(DAS)算法原理

    1. 几何原理与时延计算
    2. DAS 公式推导
  4. DAS 算法详细实现

    1. 线性阵列几何示意图
    2. 模拟点散射体回波信号
    3. DAS 时延对齐与叠加
  5. Python 代码示例与可视化

    1. 绘制阵列与焦点示意图
    2. 生成模拟回波并进行 DAS 波束形成
    3. 结果可视化
  6. 性能与优化要点
  7. 总结与延伸阅读

引言

超声成像在医学诊断、无损检测等领域被广泛应用,其核心在于如何从阵列换能器(Transducer Array)接收的原始回波信号中重建图像。波束形成(Beamforming)是将多个接收通道按照预先设计的时延(或相位)与加权方式进行组合,从而聚焦在某一空间点,提高信噪比和分辨率的方法。

延迟叠加(DAS)作为最经典、最直观的波束形成算法,其核心思路是:

  1. 对于每一个感兴趣的空间点(通常称为“像素”或“体素”),计算从这个点到阵列上每个元件(element)的距离所对应的声波传播时延;
  2. 将各通道的接收信号按照计算出的时延进行对齐;
  3. 对齐后的信号在时域上做简单加和,得到聚焦在该点的接收幅度。

本文将详细展示 DAS 算法的数学推导及 Python 实现,配合示意图帮助你更好地理解。


超声成像与束形成基础

  1. 超声成像流程

    • 发射阶段(Transmission):阵列的若干或全部换能元件发射聚焦波或游走波,激励超声脉冲进入组织。
    • 回波接收(Reception):声波遇到组织中密度变化会发生反射,反射波返回阵列,各通道以一定采样频率记录回波波形。
    • 波束形成(Beamforming):对多个通道的回波信号做时延补偿与叠加,从而将能量集中于某个方向或空间点,以提高对该点回波的灵敏度。
    • 成像重建:对感兴趣区域的各像素点分别做波束形成,得到对应的回波幅度,进而形成二维或三维图像。
  2. 阵列几何与参数

    • 线性阵列(Linear Array)平面阵列(Phased Array)圆弧阵列(Curvilinear Array) 等阵列结构,各自需要针对阵列元件位置计算时延。
    • 典型参数:

      • 元件数目 $N$。
      • 元件间距 $d$(通常为半波长或更小)。
      • 声速 $c$(例如软组织中约 $1540\~\mathrm{m/s}$)。
      • 采样频率 $f\_s$(例如 $20$–$40\~\mathrm{MHz}$)。
  3. 聚焦与分辨率

    • 接收聚焦(Receive Focus):只在接收端做延迟补偿,将接收信号聚焦于某点。
    • 发射聚焦(Transmit Focus):在发射阶段就对各换能元件施加不同的发射延迟,使发射波在某点聚焦。
    • 动态聚焦(Dynamic Focusing):随着回波时间增加,聚焦深度变化时,不断更新接收延迟。

延迟叠加(DAS)算法原理

几何原理与时延计算

以下以线性阵列、对焦在 2D 平面上一点为例说明:

  1. 线性阵列几何

    • 令第 $n$ 个元件的位置为 $x\_n$(以 $x$ 轴坐标表示),阵列位于 $z=0$。
    • 目标聚焦点坐标为 $(x\_f, z\_f)$,其中 $z\_f > 0$ 表示深度方向。
  2. 传播距离与时延

    • 声波从聚焦点反射到第 $n$ 个元件所需距离:

      $$ d_n = \sqrt{(x_n - x_f)^2 + z_f^2}. $$

    • 在速度 $c$ 的介质中,时延 $\tau\_n = \frac{d\_n}{c}$。
    • 若发射时不做发射聚焦,忽略发射时延,仅做接收延迟对齐,则各通道接收信号需要补偿的时延正比于 $d\_n$。
  3. 示意图

    线性阵列与焦点示意线性阵列与焦点示意

    图:线性阵列(横坐标 $x$ 轴上若干元件),焦点在 $(x\_f,z\_f)$。虚线表示波从聚焦点到各元件的传播路径,长度相差对应时延差。

DAS 公式推导

  1. 假设

    • 各通道采样得到离散时间信号 $s\_n[k]$,采样时间间隔为 $\Delta t = 1/f\_s$。
    • 目标像素点对应实际连续时刻 $t\_f = \frac{\sqrt{(x\_n - x\_f)^2 + z\_f^2}}{c}$。
    • 离散化时延为 $\ell\_n = \frac{\tau\_n}{\Delta t}$,可分为整数与小数部分:$\ell\_n = m\_n + \alpha\_n$,其中 $m\_n = \lfloor \ell\_n \rfloor$,$\alpha\_n = \ell\_n - m\_n$。
  2. 时延补偿(时域插值)

    • 对于第 $n$ 通道的采样信号 $s\_n[k]$,为了达到精确对齐,可用线性插值(或更高阶插值)计算延迟后对应时刻信号:

      $$ \tilde{s}_n[k] = (1 - \alpha_n) \, s_n[k - m_n] \;+\; \alpha_n \, s_n[k - m_n - 1]. $$

    • 若只采用整数延迟(或采样率足够高),则 $\alpha\_n \approx 0$,直接用:

      $$ \tilde{s}_n[k] = s_n[k - m_n]. $$

  3. 叠加与加权

    • 最简单的 DAS 即对齐后直接求和:

      $$ s_\text{DAS}[k] \;=\; \sum_{n=1}^N \tilde{s}_n[k]. $$

    • 实际中可给每个通道加权(例如距离补偿或 apodization 权重 $w\_n$):

      $$ s_\text{DAS}[k] \;=\; \sum_{n=1}^N w_n \, \tilde{s}_n[k]. $$

      常用的 apodization 权重如汉宁窗、黑曼窗等,以降低旁瓣。


DAS 算法详细实现

下面从示意图、模拟数据与代码层面逐步演示 DAS 算法。

线性阵列几何示意图

为了便于理解,我们绘制线性阵列元件位置和聚焦点的几何关系。如 Python 可视化所示:

Linear Array Geometry and Focal PointLinear Array Geometry and Focal Point

**图:**线性阵列在 $z=0$ 放置 $N=16$ 个元件(蓝色叉),焦点指定在深度 $z\_f=30\~\mathrm{mm}$,横向位置为阵列中心(红色点)。虚线表示从焦点到各元件的传播路径。
  • 横轴表示阵列横向位置(单位 mm)。
  • 纵轴表示深度(单位 mm,向下为正向)。
  • 从几何可见:阵列中心到焦点距离最短,两侧元件距离更长,对应更大的接收时延。

模拟点散射体回波信号

为直观演示 DAS 在点散射体(Point Scatterer)场景下的作用,我们用简单的正弦波模拟回波:

  1. 点散射体假设

    • 假定焦点位置处有一个等强度点散射体,发射脉冲到达焦点并被完全反射,形成入射与反射。
    • 可以简化成:所有通道都在同一发射时刻接收到对应于自身到焦点距离的时延回波。
  2. 回波信号模型

    • 每个通道接收到的波形:

      $$ s_n(t) \;=\; A \sin\bigl(2\pi f_c \, ( t - \tau_n )\bigr) \cdot u(t - \tau_n), $$

      其中 $f\_c$ 为中心频率(MHz)、$A$ 为幅度,$u(\cdot)$ 为阶跃函数表明信号仅在 $t \ge \tau\_n$ 时存在。

    • 离散采样得到 $s\_n[k] = s\_n(k,\Delta t)$。
  3. 示例参数

    • 中心频率 $f\_c = 2\~\mathrm{MHz}$。
    • 采样频率 $f\_s = 40\~\mathrm{MHz}$,即 $\Delta t = 0.025\~\mu s$。
    • 声速 $c = 1540\~\mathrm{m/s} = 1.54\~\mathrm{mm}/\mu s$。
    • 阵列元素数 $N = 16$,间距 $d=0.5\~\mathrm{mm}$。
    • 焦深 $z\_f = 30\~\mathrm{mm}$,焦点横向位于阵列中心。

DAS 时延对齐与叠加

  1. 计算每个元件的时延

    • 对第 $n$ 个元件,其位置 $(x\_n,0)$ 到焦点 $(x\_f,z\_f)$ 的距离:

      $$ d_n = \sqrt{(x_n - x_f)^2 + z_f^2}. $$

    • 对应时延 $\tau\_n = d\_n / c$(单位 $\mu s$)。
  2. 对齐

    • 对接收到的离散信号 $s\_n[k]$,计算离散时延 $\ell\_n = \tau\_n / \Delta t$,取整可先做粗对齐,如果需要更高精度可进行线性插值。
    • 例如:$m\_n = \lfloor \ell\_n \rfloor$,以 $s\_n[k - m\_n]$ 作为对齐结果。
  3. 叠加

    • 取所有通道在同一离散时刻 $k$ 上对齐后的样点,直接相加:

      $$ s_\text{DAS}[k] = \sum_{n=1}^N s_n[k - m_n]. $$

    • 对于固定 $k\_f$(对应焦点回波到达时间的离散索引),DAS 输出会在该时刻出现幅度最大的 “叠加峰”。

Python 代码示例与可视化

下面通过一段简单的 Python 代码,演示如何:

  1. 绘制线性阵列与焦点几何示意。
  2. 模拟点散射体回波信号。
  3. 基于 DAS 进行时延对齐 & 叠加。
  4. 可视化对齐前后信号与最终波束形成输出。

**提示:**以下代码在已安装 numpymatplotlib 的环境下可直接运行,展示两幅图:

  1. 阵列与焦点示意图。
  2. 多通道回波信号 & DAS 叠加波形。

绘制阵列与焦点示意图 & 模拟回波与 DAS 结果

import numpy as np
import matplotlib.pyplot as plt

# 阵列与信号参数
num_elements = 16          # 元件数量
element_spacing = 0.5      # 元件间距(mm)
focal_depth = 30           # 焦点深度(mm)
sound_speed = 1540         # 声速 (m/s)
c_mm_per_us = sound_speed * 1e-3 / 1e6   # 转换为 mm/μs
fs = 40.0                  # 采样频率 (MHz)
dt = 1.0 / fs              # 采样间隔 (μs)
f0 = 2.0                   # 中心频率 (MHz)

# 阵列元件位置 (mm)
element_positions = np.arange(num_elements) * element_spacing
focal_x = np.mean(element_positions)        # 焦点横坐标 (mm)
focal_z = focal_depth                       # 焦点深度 (mm)

# 时域采样轴
t_max = 40.0  # μs
time = np.arange(0, t_max, dt)  # 离散时间

# 模拟每个元件接收的回波信号(点散射体)
signals = []
delays_us = []
for x in element_positions:
    # 计算该通道到焦点距离及时延
    dist = np.sqrt((x - focal_x)**2 + focal_z**2)
    tau = dist / c_mm_per_us       # 时延 μs
    delays_us.append(tau)
    # 模拟简单正弦回波(t >= tau 时才有信号),幅度为1
    s = np.sin(2 * np.pi * f0 * (time - tau)) * (time >= tau)
    signals.append(s)

signals = np.array(signals)
delays_us = np.array(delays_us)

# DAS 对齐:整数时延补偿
delay_samples = np.round(delays_us / dt).astype(int)
aligned_signals = np.zeros_like(signals)
for i in range(num_elements):
    aligned_signals[i, delay_samples[i]:] = signals[i, :-delay_samples[i]]

# 叠加
beamformed = np.sum(aligned_signals, axis=0)

# 可视化部分
plt.figure(figsize=(12, 8))

# 绘制阵列几何示意图
plt.subplot(2, 1, 1)
plt.scatter(element_positions, np.zeros_like(element_positions), color='blue', label='Array Elements')
plt.scatter(focal_x, focal_z, color='red', label='Focal Point')
for x in element_positions:
    plt.plot([x, focal_x], [0, focal_z], color='gray', linestyle='--')
plt.title('Line Array Geometry and Focal Point')
plt.xlabel('Lateral Position (mm)')
plt.ylabel('Depth (mm)')
plt.gca().invert_yaxis()  # 深度向下
plt.grid(True)
plt.legend()

# 绘制模拟回波(示例几路通道)与 DAS 叠加结果
plt.subplot(2, 1, 2)
# 仅展示每隔 4 个通道的信号,便于观察
for i in range(0, num_elements, 4):
    plt.plot(time, signals[i], label=f'Raw Signal Element {i+1}')
plt.plot(time, beamformed, color='purple', linewidth=2, label='Beamformed (DAS)')
plt.title('Received Signals and DAS Beamformed Output')
plt.xlabel('Time (μs)')
plt.ylabel('Amplitude')
plt.xlim(0, t_max)
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

代码说明

  1. 阵列几何与时延计算

    dist = np.sqrt((x - focal_x)**2 + focal_z**2)
    tau = dist / c_mm_per_us
    • 先在平面中以 mm 为单位计算距离,再除以声速(mm/μs)得到回波时延(μs)。
  2. 生成点散射体回波

    s = np.sin(2 * np.pi * f0 * (time - tau)) * (time >= tau)
    • 采用简单的正弦信号模拟中心频率 $f\_0$ 的回波脉冲,实际系统可使用窗函数调制波包。
    • (time >= tau) 实现“在 $t < \tau$ 时无信号”(零填充)。
  3. DAS 对齐

    delay_samples = np.round(delays_us / dt).astype(int)
    aligned_signals[i, delay_samples[i]:] = signals[i, :-delay_samples[i]]
    • 将连续时延 $\tau$ 转为离散采样点数 $\ell = \tau/dt$,近似取整为整数延迟 $m = \lfloor \ell + 0.5 \rfloor$。
    • 整数对齐简单易行,但若需更高精度可插值。
  4. 叠加与可视化

    • 将对齐后的所有通道信号在时域上直接相加,形成 beamformed
    • 在第二幅图中,将若干通道的原始信号(尖峰位置不同)与叠加结果(峰值一致聚焦)放在同一子图,突出 DAS 聚焦效果。

结果可视化

运行上述代码后,你将看到两幅关键图像:

  1. 线性阵列与焦点示意图

    • 蓝色叉代表阵列上均匀分布的 16 个换能元件;
    • 红色叉代表聚焦点(深度 30 mm);
    • 虚线从各元件到焦点,直观说明不同元件回波时延不同。
  2. 多通道回波与 DAS 叠加输出

    • 上半图展示几个示例通道(如元素 1、5、9、13)的模拟回波信号,明显看到每路信号的到达时间不同;
    • 下半图(紫色曲线)为 DAS 对齐后加和的输出,在某一时刻出现峰值,说明成功聚焦到点散射体。

性能与优化要点

  1. 插值精度

    • 直接用整数时延对齐(附近点取值)简单,但会有量化误差;
    • 更精准的做法是线性插值或更高阶插值,对时延进行亚采样点对齐:

      $$ \tilde{s}_n[k] = (1-\alpha) s_n[k - m] + \alpha \, s_n[k - m -1],\quad \alpha \in [0,1]. $$

    • 插值虽能提升分辨率,但计算量增大。
  2. 加权策略(Apodization)

    • 为了抑制旁瓣,可以给每个换能元件一个加权系数 $w\_n$,如汉宁窗、黑曼窗:

      $$ s_\text{DAS}[k] = \sum_{n=1}^N w_n \, \tilde{s}_n[k]. $$

    • 通常 $w\_n$ 关于阵列中心对称,可以降低非焦点方向的能量。
  3. 动态聚焦

    • 当对不同深度进行成像时,焦点深度不断变化,每个深度都需要重新计算时延并叠加;
    • 实时成像时,需要针对每个像素点(或像素列)循环做 DAS,计算量大,可使用 GPU 加速或 FPGA 硬件实现。
  4. 多发多收与合成孔径

    • 不同聚焦位置往往需要多次发射(Tx)与接收(Rx),可合成多个 Tx-Rx 事件得到更复杂的波束合成。
    • 合成孔径(Synthetic Aperture)方式会在信噪比和分辨率上更出色,但更耗时。
  5. 并行加速

    • 在 CPU 上逐点做 DAS 速度较慢,可使用 GPU 或 SIMD 指令并行化:

      • 每个像素对应的多个通道时延计算、信号对齐与加权都可并行;
      • 多深度或多方向的计算也易并行分配。

总结与延伸阅读

  • DAS(Delay-and-Sum) 是经典、直观且易实现的超声波束聚焦算法,通过对各通道回波信号进行时延补偿后相加,实现空间聚焦。
  • 从几何原理到公式推导,再到 Python 代码可视化,本文详尽展示了 DAS 在点散射体场景下的原理与效果。
  • 实际超声成像中,需要动态聚焦、加权(Apodization)、插值对齐与多发多收策略等手段,以提升分辨率和旁瓣抑制。

延伸阅读建议:

  1. Jensen, J.A., “Field: A Program for Simulating Ultrasound Systems”, Medical & Biological Engineering & Computing, 1996.
  2. Boukerroui, D., Yessad, A.C., et al. “Ultrasound Beamforming: An Overview of Basic Concepts and State-of-the-Art in Fast Algorithms”, IEEE Access, 2020.
  3. Szabo, T.L., “Diagnostic Ultrasound Imaging: Inside Out”, 2nd Edition, Academic Press, 2013.
  4. 李庆等,《超声成像与成像技术》,科学出版社,2018。

评论已关闭

推荐阅读

DDPG 模型解析,附Pytorch完整代码
2024年11月24日
DQN 模型解析,附Pytorch完整代码
2024年11月24日
AIGC实战——Transformer模型
2024年12月01日
Socket TCP 和 UDP 编程基础(Python)
2024年11月30日
python , tcp , udp
如何使用 ChatGPT 进行学术润色?你需要这些指令
2024年12月01日
AI
最新 Python 调用 OpenAi 详细教程实现问答、图像合成、图像理解、语音合成、语音识别(详细教程)
2024年11月24日
ChatGPT 和 DALL·E 2 配合生成故事绘本
2024年12月01日
omegaconf,一个超强的 Python 库!
2024年11月24日
【视觉AIGC识别】误差特征、人脸伪造检测、其他类型假图检测
2024年12月01日
[超级详细]如何在深度学习训练模型过程中使用 GPU 加速
2024年11月29日
Python 物理引擎pymunk最完整教程
2024年11月27日
MediaPipe 人体姿态与手指关键点检测教程
2024年11月27日
深入了解 Taipy:Python 打造 Web 应用的全面教程
2024年11月26日
基于Transformer的时间序列预测模型
2024年11月25日
Python在金融大数据分析中的AI应用(股价分析、量化交易)实战
2024年11月25日
AIGC Gradio系列学习教程之Components
2024年12月01日
Python3 `asyncio` — 异步 I/O,事件循环和并发工具
2024年11月30日
llama-factory SFT系列教程:大模型在自定义数据集 LoRA 训练与部署
2024年12月01日
Python 多线程和多进程用法
2024年11月24日
Python socket详解,全网最全教程
2024年11月27日