拓十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

射频数据驱动的颈动脉超声分割与IMT精准测量

射频数据驱动的颈动脉超声分割与IMT精准测量

简介:基于MATLAB的射频数据颈动脉超声分割资料包,内容面向医学图像处理与超声信号分析方向的学生和入门研究者,系统覆盖RF数据理解、滤波去噪、特征提取、分割算法选择到评估优化的完整流程。资源共45个文件,以43个M脚本为主,另含1个FDA滤波器文件和1个Markdown说明文档,压缩包仅47KB,代码轻量、适合直接运行与对照阅读。已有163人学习浏览,对正在了解超声射频数据的读者具有一定参考价值。资料中除了Butterworth/FIR、自适应去噪等预处理方法,还介绍了希尔伯特变换包络提取、主动轮廓模型、阈值分割、区域生长、水平集以及结合SVM或随机森林的分割策略,并给出mylowpassfilter.m等可复用滤波脚本,帮助读者搭建完整的实验框架。此外利用Dice相似系数、Jaccard相似度等指标与手动标注对比,便于验证并提升分割精度。

1. 射频数据颈动脉超声分割:先把“原始信号”和“临床测量”之间的账算明白

做颈动脉超声分割的同行,绝大多数一开始都在B模式图像上画内中膜边界。B模式是超声设备把探头收到的射频信号经过波束合成、包络检波、对数压缩之后渲染出来的灰阶图,好看,但丢掉的原始信息远比我们以为的多。射频数据则是探头阵元直接采到的原始信号,幅度、相位、频率成分都在,理论上比B模式更适合区分内膜、中膜和外膜,尤其是管腔-内膜边界那种强反射界面。

给体检中心或者心内科做IMT自动测量时,医生要的不是一张漂亮的分割图,而是“这条血管壁厚度到底是多少毫米”。射频数据分割的价值就在于它有可能把IMT测量做得比手工标注更稳、更细。适合读这篇的人应该是医学影像算法工程师、超声设备厂商算法组,以及想摆脱繁琐手工标注的研究组。我会按一套可落地的技术路线来拆这件事:数据从哪来、预处理怎么做、模型怎么选、参数怎么调、坑在哪里。

2. 为什么拿RF数据做分割,而不是直接用B模式图像?

2.1 RF数据里到底多存了哪些信息

超声探头阵元采集到的射频信号,是组织对入射声波的反射回波叠加。频率范围通常落在2到12MHz,颈动脉表浅,常用7到10MHz线性探头。RF数据如果做了正交解调,会以IQ复数形式存下来,实部和虚部分别对应余弦和正弦分量的幅度;更原始的直接采集方式存实数采样点,每个采样点间隔由采样率决定,比如40MHz采样率下,相邻采样点对应约0.019mm的声传播距离(按声速1540m/s计算)。

B模式图像在生成时,会对RF包络做对数压缩,把动态范围从几十dB压到人眼可分辨的灰阶。这个过程保住了“回波强不强”,丢了相位和局部频率变化。而射频数据里同时保留了相位和幅度。相位对组织界面的微小位移极其敏感,0.01mm量级的位移都能体现为可测的相位旋转。对于颈动脉壁这种随心跳搏动的结构,相位信息意味着分割模型有机会学到“这个像素点是不是正处于搏动周期中的某一段”,这在B模式单帧图上是看不到的。

2.2 颈动脉内中膜边界在射频信号里的物理表现

临床上测IMT,标准位置是颈总动脉远壁,也就是远离探头的那侧血管壁。在射频信号里,这个区域有两个特征强烈的反射界面。

第一个是管腔血液与内膜的界面,血液是散射体,内膜是相对均质的组织层,这个界面临近声束方向时,回波极性从杂乱散射变成强正反射,幅度突变明显。第二个是中膜与外膜的界面,中膜是平滑肌,声阻抗跟外膜结缔组织差异不大,但射频信号的频谱在这个位置会出现规律性的高反射峰会与管腔侧弱得多的回波成为对照。两个界面之间的几何宽度就是内中膜厚度,正常值在0.6到1.0mm之间,在40MHz采样率下只有30到50个采样点,折算到B模式图像上往往只有3到5个像素。

这就是为什么射频分割在IMT任务上有结构性优势。B模式图像上的边界已经过点扩散函数模糊和对数压缩,边缘定位精度天生受限;而在RF域,边界是两个峰值之间的整数个采样点,物理位置可计算。许多医生手工测量IMT时反复放大图像取均值,本质上就是在跟这个定位精度较劲。

2.3 一条可落地的射频分割技术栈

整套流程可以写成一个流:线性阵列探头采集RF回波 → 按需做波束合成或直接用阵元数据 → 带通滤波去除频带外噪声 → 提取包络/相位特征 → 深度衰减补偿 → 送入分割网络 → 输出内膜和中膜外膜的逐像素掩膜 → 计算几何宽度得到IMT。

是否做波束合成是一个关键分叉。很多研究平台(比如Verasonics这类开放研究系统)能直接导出波束合成后的RF帧,每条线对应一个组织深度的一维信号。这类数据与B模式图像的几何关系简单,工程上最容易落地。如果真的拿到多阵元原始回波,则还要先做延时叠加,否则后续模型要多学一个波束形成的隐含关系。我的经验是,除非你有明确的相干成像算法要研究,否则直接用设备输出的波束合成后RF帧,也就是每条扫描线一个一维信号,能把工程复杂度降一半以上。

3. 把RF数据变成模型输入:预处理与数据集的完整做法

3.1 从采集设备导出RF数据的常见形式与读取方法

不同厂商导出的数据组织形式完全不同。常见有三种:int16实数RF采样、float32实数RF采样、float32复数IQ数据。有些设备把每帧组织成二维数组,形状是(采样点数, 扫描线数),有些则是把多条扫描线打包在一段连续内存里。拿到数据后第一件事不是训练,而是把数据格式、采样率、阵元中心频率、探头深度范围、线间距全部记录在数据集元数据里,缺了任何一项后面都会翻车。

我一般先把原始数据读成numpy数组,再做一次“冒烟检查”,用包络波形肉眼确认结构符合颈动脉解剖特征。

import numpy as np import scipy.signal as signal def load_rf_frame(filepath, n_lines, n_samples, dtype=np.float32): """ 读取单个RF帧。 n_lines: 扫描线数量,例如128条线 n_samples: 每条线深度方向采样点数,例如3117 注意:数据按(endian)保存,如果对不上请转换字节序 """ raw = np.fromfile(filepath, dtype=dtype) assert raw.size == n_lines * n_samples, \ f"文件长度 {raw.size} 与 {n_lines}x{n_samples} 不符" frame = raw.reshape(n_samples, n_lines, order='F') return frame # (深度采样点, 线数) rf = load_rf_frame("carotid_01_0001.bin", n_lines=128, n_samples=3117) print(rf.shape, rf.dtype, np.abs(rf).max())

这段代码的关键点是order='F'。许多设备按扫描线顺序写入数据,即先写完整第一条线的所有深度采样点,再写第二条线。如果文件的写入逻辑是列优先,而读取时用默认的C顺序reshape,波形会错位,训练出的模型必然学不到任何有效特征。字节序也是高频坑,int16数据经常是小端存储,在Windows和Linux上读出来符号位会反转,波形表现为剧烈振荡的噪声。

3.2 预处理管线:带通滤波、包络提取、归一化与纵深增益补偿

RF数据不能直接进分割网络。原始射频信号的频带中包含探头中心频率附近的组织回波,同时混有低频漂移和高频噪声。带通滤波的截止频率通常设置为中心频率的0.5到1.5倍。对于7MHz探头,我会用3.5MHz到10.5MHz的Butterworth带通滤波器,阶数4到6,避免过重的相位畸变。

def preprocess_rf(rf_frame, fs, fc, c=1540.0): """ rf_frame: (深度采样, 线数) fs: 采样率, 单位Hz fc: 探头中心频率, 单位Hz """ nyq = 0.5 * fs low = max(0.1, 0.5 * fc) / nyq high = min(0.95, 1.5 * fc) / nyq b, a = signal.butter(4, [low, high], btype='band', analog=False) filtered = signal.filtfilt(b, a, rf_frame, axis=0) # 零相位滤波 analytic = signal.hilbert(filtered, axis=0) # 解析信号 envelope = np.abs(analytic) # 包络 # 深度衰减补偿。接收回波随深度按指数衰减,近似补偿系数为 exp(2*alpha*depth*fc) depth = np.arange(rf_frame.shape[0]) / fs * c / 2.0 # 单位为m alpha = 0.5 # 经验衰减系数,单位 dB/(cm*MHz),可根据探头和频率调整 atten_db = 2.0 * alpha * depth * 1e2 * (fc / 1e6) tgc = np.exp(atten_db / 8.686) # dB转线性幅度 envelope_tgc = envelope * tgc[:, np.newaxis] # 对数压缩并线性归一化到[0, 1] log_env = np.log1p(envelope_tgc) out = (log_env - log_env.min()) / (log_env.max() - log_env.min() + 1e-6) return out.astype(np.float32)

逻辑说明分三层。第一层filtfilt做零相位带通滤波,保留了边界位置不漂移,这点对分割任务比最小相位滤波器更友好。第二层hilbert把实数RF信号变换为复数解析信号,取模得到包络,这就是B模式图像前身,但此时还没做对数坐标和后处理。第三层深度补偿是整个管线最容易错的地方:高频超声在组织中的衰减大约是每厘米每MHz衰减0.5dB,一个7MHz的探头在2cm深处已经损失了14dB回波强度,远端管壁的包络幅度会明显低于近端。如果不做补偿,模型很容易学成“只分割近端”的偏置。

3.3 标签对齐:如何把医生的B模式标注映射到RF坐标系

标注通常是医生在B模式超声图像上画的线,B模式图像的每个像素和RF帧的每个采样点之间不是简单的一一对应。线性阵列探头还好,像素的横向坐标对应扫描线索引,深度方向需要按像素分辨率反推采样点索引。但如果是凸阵探头或扇形扫描,B模式已经做了坐标变换,直接拿来用会整体偏移。

正确做法是记住设备导出的B模式图像规格,通常图像尺寸是W像素×H像素,对应的物理范围是横向宽度W_pitch×N_lines,深度范围D。通过线性映射把标注线的像素坐标换算成扫描线和深度采样点坐标。

def label_bmode_to_rf(points_px, bmode_h, rf_n_samples, depth_m, oversample=1.0): """ points_px: 标注点列表, 形状(N, 2), 第一列是深度像素, 第二列是横向像素 bmode_h: B模式图像深度方向像素数 rf_n_samples: RF帧深度采样点数 depth_m: RF帧对应的物理深度, 单位m """ sample_per_m = rf_n_samples / depth_m px_per_m = bmode_h / depth_m rf_coord = points_px[:, 0] / px_per_m * sample_per_m return rf_coord.astype(np.float32)

这段代码隐含一个前提:B模式图像深度方向和RF深度方向是一一对应的,只是分辨率不同。实际如果设备对RF数据做了纵向压缩或裁剪,单纯缩放就会错位。我一般用一个线模(wire phantom)去标定:把一根细线放在几个已知深度位置,扫描后分别在B模式图和RF包络上找它的位置,反推出两套坐标系的齐次映射。这个标定步骤看着繁琐,但能避免训练集里一半标签偏移一两个像素的问题,对于IMT这种只有几毫米的结构,一两个像素的误差就能让实验彻底失去意义。

4. 分割网络选型与训练参数:用数据规模定模型复杂度

4.1 用UNet改造射频多通道输入:通道怎么排、patch怎么裁

预处理后的包络只是一维信息,RF数据里能被模型利用的还有相位和瞬时频率。最稳妥的输入设计是把包络、归一化瞬时频率、深度位置索引堆叠成三通道图,前两个来自RF信号,第三个给模型一个显式的物理坐标提示。

import torch.nn as nn class RFSegUNet(nn.Module): def __init__(self, in_channels=3, out_classes=2): super().__init__() self.encoder = nn.Sequential( nn.Conv2d(in_channels, 32, kernel_size=3, padding=1), nn.BatchNorm2d(32), nn.ReLU(inplace=True) ) # 此处省略UNet中下采样、上采样等常规模块 # head输出两类:内膜边界和中外膜边界 def forward(self, x): x = self.encoder(x) return x

网络结构本身不特殊,特殊在输入。我需要强调的是三个通道的来源:edge_prob通道是包络经高斯拉普拉斯滤波后的响应,强化边界位置;inst_freq通道估计局部瞬时频率的偏移,反映组织衰减特性;depth通道直接填入深度归一化值。B模式图像里没有哪个通道能替代原始相位信息,所以这类三通道输入的提升主要在薄结构召回率上。

patch裁剪直接决定显存占用和边界上下文。颈动脉腔直径约5到7mm,加上前后壁组织,一个包含完整远壁的patch大约需要覆盖深度方向20到30mm。配合RF帧约0.019mm每采样点的轴向分辨率,也就是1000到1500个采样点,在横向上包含64到128条扫描线,输出尺寸为1024×128。这个尺寸对UNet已经是256倍下采样的上限,显存不足时优先沿深度方向切块,不要横向压缩。

4.2 损失函数与评价指标选择:IMT任务里DICE不是唯一标准

分割网络默认用DICE Loss,但对IMT这种厚度只有几个像素的带状结构,DICE对整体覆盖度敏感,对边缘位置不敏感。一个预测结果哪怕边界整体向外扩了一个像素,DICE可能只下降2%到3%,IMT测量误差却会放大到0.2mm以上,临床不可接受。

我推荐的组合是主损失DICE加上轮廓损失,轮廓损失可以简单用L1距离度量预测边界与标注边界之间的像素距离:

def boundary_aware_loss(pred_mask, true_label, true_boundary_dist): """ pred_mask: (B, 2, H, W) 网络输出经sigmoid的预测 true_label: (B, 2, H, W) one-hot标注 true_boundary_dist: (B, 1, H, W) 距离变换图 对接近真实边界的像素给予更高权重 """ weights = 1.0 + 3.0 * (true_boundary_dist < 5).float() ce = nn.functional.binary_cross_entropy(pred_mask, true_label, reduction='none') weighted_ce = (ce * weights).mean() dice = 1 - (2 * (pred_mask * true_label).sum(dim=(2, 3)) + 1) / \ ((pred_mask + true_label).sum(dim=(2, 3)) + 1) return weighted_ce + dice.mean()

这段代码里的距离变换图是预先用标注轮廓距离生成的,近边界五像素范围内权重提为4倍,其余位置权重为1。这样做会让网络优先把边界位置学准,而不是用一个粗略的覆盖去平摊损失。

评价指标在工科论文里当然要报DICE和IoU,但面向临床我还会同时报两项物理量:一是预测IMT与医生标注IMT的均值绝对误差(mean absolute error,MAE);二是Bland-Altman一致性界限(limits of agreement,LoA)。后者才能回答“这个分割结果能不能当测量工具用”,DICE在主管面前永远说不过这个。

4.3 训练超参数一页纸:学习率、步长、增强策略与显存控制

基于十来个病例数据起步的项目,初始条件下给出下面这组参数通常能跑通:批大小2,patch尺寸1024×128,AdamW优化器,初始学习率3e-4,权重衰减1e-5,使用余弦退火学习率调度器,训练轮次200。如果是自建小数据集,每轮做随机深度偏移、横向抖动、亮度扰动,注意不要做纵向拉伸,因为IMT的物理厚度不能随便缩放。

显存控制有两个诀窍。第一是不要一次调度整帧进GPU,沿深度方向随机切出1024点、在线方向取128线,上下文已经足够让网络看到完整的管腔前后壁。第二是梯度累积,batch size设为2实际心里计算梯度时累积到8再更新一步,模拟出等效batch size为8的效果,收敛稳定性好很多。RF数据量本来就比B模式大一个量级,预处理后可以直接存包络特征图,形状约每帧1024×128×3,一个病例上千帧占几个GB,切成patch存成h5或npy,训练时不临时现算,性能要稳得多。

5. 射频分割最常见的5个采坑:现象、原因、解决

5.1 训练DICE很高,测试IMT误差反而更大

我在这上面翻过一次车。模型分割结果和标注掩膜重叠率超过0.95,但算出来的IMT比医生手工值平均偏大0.3mm。原因在于DICE衡量的是区域覆盖度,但IMT测量的本质是两条边界之间的距离,这个物理量对边界整体内缩或者外扩极其敏感,跟覆盖度相关性弱。解决方法是改用带边界权重的损失函数,并且在验证集上以IMT绝对误差为早停标准,DICE只能作为次要观察指标。

5.2 近端管壁分割很好,远壁总是“漏底”

超声信号经过管腔血液散射后大幅衰减,即使加了固定TGC补偿,远壁回波仍可能比近壁低6dB以上。现象就是分割掩膜在中膜外膜边界上断断续续,尤其在血管较深或探头角度倾斜时更明显。解决思路有两个:一是把固定TGC补偿换成自适应增益,按每个采样点局部噪声底噪估计补偿增益;二是在训练时对远壁区域做双倍重采样,强制模型看到足够多的远壁阳性样本。后者操作更简单,通常立竿见影。

5.3 换一台设备或换一个探头型号,效果断崖式下跌

RF数据不像B模式图像那样有相对一致的视觉风格。不同设备的波束形成器、孔径切换逻辑、增益曲线差异极大,同一台机器换一个中心频率的探头,包络特征分布就变了。最直接的表现是训练时验证集DICE有0.93,拿到另一台机器采集的数据只剩0.75。解决路径分两步:预处理阶段做逐帧的幅值直方图匹配,把目标帧的包络分布对齐到训练集的聚合分布;训练阶段加入MixStyle或者通道级随机噪声增强,让网络不依赖某台设备的固定增益模式。跨中心的验证协议必须从项目第一天就设计进去,不要等到医院排队扫描了才发现不可用。

5.4 标注线和分割图高度重叠,但成对的内膜/外膜边界识别错了

颈动脉IMT远壁有两条边界,近管腔的是内膜,深一点的是中外膜。某些帧中外膜回波弱,B模式图像也看不太清,模型就把内膜同时当成了两条边界。DICE看起来还行,IMT值错了一半。我用的解决办法是在模型输出层之外加一个解剖先验检查器,预测两条边界的相对深度,要求外膜边界深度必须大于内膜边界,位置差不小于0.4mm、不大于2.0mm。违反这个约束的输出直接判为无效帧,交给系统重新做局部细化,而不是盲目上报告。

5.5 明明一帧一帧训练收敛良好,连续视频分割却像抽风一样跳

单帧分割的真实帧往往在相邻帧间剧烈变化,尤其收缩期管壁快速运动时,边界预测在相邻帧之间来回摆动,视频看起来像抖动。深层原因是训练数据里相邻帧来自不同心跳相位,模型学到的是“每个相位各自最可能的边界”,却没有学到边界在时间上的连续性。解决方案是训练阶段加入相邻帧一致性正则化,要求模型对同一条血管在时间上相近两帧的输出边界距离小于一个阈值;推理阶段再用指数移动平均对边界坐标做轻量平滑,平滑系数取0.6到0.8,在延迟和稳定性之间折中。

6. 验证与进阶:用多帧RF序列做时序分割,才算真正面向临床

6.1 多帧一致性评估与平滑

单帧指标漂亮还不够,临床测量通常在包含多个心动周期的视频片段上进行。我建议把所有病例按帧序列组织,额外计算三个指标:边界抖动幅度(相邻帧边界位置差分的标准差)、有效帧比例(通过解剖先验检查的帧占比)、以及整段视频IMT均值与医生在关键帧手工测量的差值。数值标准参考经验值:边界抖动小于0.1mm可以接受,有效帧比例低于80%说明预处理或模型稳定性有问题。

6.2 跨设备验证协议

只在自己科室数据上验证没有说服力。设计方案时最好让数据按设备和采集人员分层划分,至少保证验证集里包含一台训练时没见过的机器。报告呈现上列出每台设备的IMT MAE、LoA、DICE三个指标,比一个大一统的均值有意义得多。准备一份数据清单,逐项记录设备型号、探头中心频率、采样率、是否做TGC、是否做滤波,排查域偏移时这张表能救命。

6.3 用B模式知识辅助RF域数据集

RF数据标注成本高是落地卡点。一个便宜可靠的进阶路子是用同一患者的B模式图像做伪标签:另训练一个B模式分割模型,在大量历史B模式数据上学到稳健的IMT边界定位能力,然后把当前患者的B模式标注通过坐标映射转换到RF域,作为低置信度样本参与训练。我在实践中用这种半监督方式把小样本项目的数据规模扩了三倍,模型在换机场景下的稳定性提升明显。最终端的实用技巧是把“标注一致性检查”做成离线脚本,自动抓取两条边界厚度超出0.4到2.0mm物理范围的样本做复核。项目做久了最大的体会是,医生手工测量本身也有标准差,早期我以为模型应该无限逼近某个“真值”,后来学会用多帧均值和一致性界限去校准,反而报出的结论更可信。希望这一篇能帮你的射频分割项目少走一段弯路,落地顺畅。

本文还有配套的精品资源,点击获取

返回列表