简介基于射频数据的颈动脉超声分割MATLAB方案面向医学图像处理研究者、超声成像相关专业学生及工程技术人员系统覆盖射频信号预处理、特征提取、分割算法选择与评估优化等环节。资源内包含43个.m源码文件实现Butterworth/FIR滤波、希尔伯特变换包络提取、峰值检测、主动轮廓模型等关键技术并涵盖区域生长、阈值分割、水平集以及支持向量机、随机森林等机器学习分割思路另有1个.fda滤波器定义用于低通滤波配置1个.md说明文档梳理实验结构与代码用法。压缩包共45个文件整体仅47KB轻量且目录清晰。目前已有163人学习通过该资源可以快速搭建射频数据处理流程理解分割算法在超声图像上的落地方法并参考其参数设置与Dice、Jaccard等评估指标适合希望借鉴完整MATLAB实现的研究者与开发者。1. 射频数据做颈动脉超声分割为什么值得从B超图退回原始信号颈动脉超声分割这个任务放在B超图上做远远没有听起来那么顺利。颈动脉前壁和后壁的内膜、中膜、外膜三层在屏幕上只差几个像素灰度值同一个患者换一台设备窗宽窗位一变灰度梯度整个换掉。我最初只是在显示图像上抠轮廓Dice稳定在75%怎么调模型都上不去。后来把原始射频数据RF data拿来做输入同样一组标注分割Dice爬到82%以上。射频数据是探头接收到的回波信号比B超图多保留了幅度细节和相位信息尤其是内膜-中膜那种弱回声边界RF包络的梯度变化比灰度图可靠得多。下面要拆解的是一条完整路径从RF数据怎么读取和预处理到用什么分割网络和损失再到评估、避坑以及最后如何从分割结果里量出颈动脉IMT。适合算法工程师、超声研究组和想用原始信号做产品的开发者。2. 射频数据到训练张量包络、IQ双通道与坐标对齐2.1 先分清手里的RF是什么三种数据形态差别很大超声设备能导出的数据至少有三层。最底层的原始RF是探头每个阵元收到的回波经放大和AD采样后的实数序列采样率通常是中心频率的3到4倍比如7.5MHz探头配40MHz采样。中间一层是IQ基带数据把RF信号和本振混频后低通得到一对I,Q复数值等效于每个采样点用两个实数表示幅度和相位。最上层是B超图也就是把IQ包络检测、对数压缩、扫描变换后得到的显示灰度图。做颈动脉分割时直接用B超图最省事但也最容易遇到机器相关性问题B超图为了“看起来平滑”把80dB的原始动态范围压到人眼能接受的20dB左右弱回声区域被压到和噪声底一个级别。颈动脉壁的内膜-中膜边界恰恰就是弱回声向中等回声的过渡区在几步压缩后灰度梯度变钝。而RF或IQ数据保留了这个梯度的原始量化值网络能从更宽的动态范围里学到边界。另一个常见情况是设备导出的“RF数据”实际已经是IQ。比如研究模式下一帧IQ大小可能是256条扫描线2048个深度点每条线对应一个复信号。这时不需要自己重做正交解调I和Q两个通道直接可以作为CNN输入。要做的只是把坐标和B超图对应起来。下面先按最常见的一维RF到包络的处理写。2.2 从RF帧到包络图预处理流程与常用参数把RF数据变成一张能进U-Net的单通道包络图核心步骤是去直流、带通滤波、包络检测、对数压缩和归一化。下面这段代码我通常放在数据加载模块里import numpy as np from scipy.signal import hilbert, butter, sosfiltfilt # rf_data: shape (num_scanlines, num_depth_samples) # 很多设备导出的是int16先转float再把按扫描线存储的帧转成二维矩阵 rf np.load(carotid_rf_001.npy).astype(np.float32) # 1) 去直流 rf rf - rf.mean(axis1, keepdimsTrue) # 2) 带通滤波只保留探头中心频率附近的频带 fs 40e6 # AD采样率 fc 7.5e6 # 探头中心频率 sos butter(4, [0.6 * fc, 1.2 * fc], btypebandpass, fsfs, outputsos) rf sosfiltfilt(sos, rf, axis1) # 3) 包络检测沿深度维做hilbert envelope np.abs(hilbert(rf)) # 4) 对数压缩 按帧归一化到[0,1] log_env 20.0 * np.log10(envelope 1e-6) log_env (log_env - log_env.min()) / (log_env.max() - log_env.min() 1e-8)逻辑说明去直流是消掉信号里的固定偏移否则会在深度方向形成一条亮带带通滤波用四阶Butterworth通带设在0.6到1.2倍中心频率把高于探头频响的噪声扔掉hilbert在深度维计算返回复解析信号模就是射频包络。对数压缩是模拟B超的做法但这里保留了更宽的动态范围归一化按帧独立算能减轻不同帧之间的平均强度差异。参数说明fs必须等于设备AD采样率通带上下限按探头带宽改如果探头发射脉冲带宽较小建议下限收到0.7倍中心频率上限收到1.1倍。sosfiltfilt对整帧做零相位滤波会让边界略微平滑但换来的是包络不容易出现相位失真。如果文件很大建议把这步做成离线预处理保存成npy训练时直接读包络图。2.3 IQ双通道输入保留相位信息的另一种思路用包络图会丢掉回波相位而相位信息在中膜和外膜分界上起着“方向性”的作用。如果手头是IQ基带数据我一般直接喂双通道iq np.load(carotid_iq_001.npy) # shape (scanlines, depth, 2) i iq[..., 0].astype(np.float32) q iq[..., 1].astype(np.float32) # 按帧归一化注意I/Q要共用比例 scale np.max(np.sqrt(i**2 q**2)) 1e-8 i / scale q / scale # 变成 (2, H, W) 的浮点张量 input_tensor np.stack([i, q], axis0) # (2, depth, scanlines)如果你只有原始实数RF又想得到IQ可以做一次正交解调但我不太建议自己在工程里重新造这个轮子。原因是正交解调需要发射脉冲的中心频率和采样率精确匹配稍微偏一点就会在图像里留下严重的条带伪影。常见做法是先用设备的IQ导出或者用上一节的hilbert包络作为单通道替代。逻辑说明共用scale而不是分别归一化I和Q是为了保持瞬时幅度sqrt(I^2Q^2)的相对关系否则会扭曲组织反射强度。参数说明如果IQ幅度分布非常不均匀可以取第99百分位数做scale避免个别强反射点把所有通道压暗。输入张量里H是深度维W是扫描线维神经网络不关心这两个维度的语义差异只要训练和推理保持一致即可。2.4 标注坐标映射线阵探头怎么从B超图映射回RF标注阶段多半还是在B超图上画轮廓然后要把轮廓坐标送到RF/IQ帧里。线阵探头的映射相对简单B超图的水平坐标和扫描线序号基本一一对应深度方向需要换算采样点。# 假设B超图上某个标记点depth_px对应显示图像的第depth_px行 mm_per_px 0.1 # 显示图像每个像素对应的毫米 depth_mm depth_px * mm_per_px sonic_speed 1540 # 软组织平均声速m/s fs 40e6 # 采样率 sample_idx int(depth_mm * 1e-3 * 2 * fs / sonic_speed)这个公式里的2来自超声的往返时间发射到接收走了一个来回。计算出的sample_idx就是RF矩阵深度维的标号。注意这个映射只在B超图未经扫描变换、且深度轴没有经过插值时才严格成立。很多设备显示图上会做横向插值和动态范围增强像素坐标和原始深度不是线性关系。遇到这种情况我会先在包络图上找两个强反射点血管前后壁回声用线性回归拟合显示深度到sample_idx的偏移再用于所有标注点。逻辑说明深度转换的顺序是从像素到毫米再从毫米到采样点中间不要跳过声速。凸阵探头则不能这么算凸阵B超图是扇形扫描变换后的直角坐标必须用极坐标逆变换把标注点转回“角度-深度”域再对应扫描线。这也是为什么我建议颈动脉分割优先选线阵临床测IMT本来就用线阵映射干净算法才能把边界误差控制在0.1mm以内。预处理之后下一步就是让模型直接消费这些张量。3. 分割网络与损失函数从U-Net开始相位信息怎么用3.1 为什么不是自己提特征而是让网络直接吃RF包络早年做超声斑块分割常见思路是算一堆纹理特征比如灰度共生矩阵、局部二值模式再喂给分类器。这套路的瓶颈不在分类器而在特征只能刻画局部小窗口。颈动脉壁是一个纵向连续、厚度只有0.6到1mm的结构需要上下文把前后壁连接起来。U-Net这类编码器-解码器结构天然能聚合从3像素到上百像素的感受野。RF或IQ输入相对B超图来说视觉上信噪比显得差但这种“差”是动态范围大带来的误解。网络里的卷积可以自动选择它需要的频段并不需要人先做一套复杂的特征工程。我的建议是直接用包络或IQ双通道输入先跑通单通道包络再换IQ双通道对比。如果IQ并没有带来提升也不要惊讶——很多数据集里B超图本身已经保留了足够的纹理信息RF的优势更多体现在跨设备泛化上。另外U-Net并不是唯一选择但对颈动脉分割来说它足够。Attention U-Net、U-Net会有收益但它们更适合边界极其模糊的病灶。RF数据里的血管壁层状边界相对锐利普通U-Net加上一个合理的损失函数就能跑出可信结果。先把这条基线做扎实再谈加模块。3.2 一个能跑的最小U-Net输入两通道IQ输出四分类下面这个U-Net为了快速试错通道数压得很低。它的输入是B,2,H,W输出是B,4,H,W对应背景、管腔、内膜中膜复合体、外膜四个类别。如果你只有包络图把in_channels改成1。import torch import torch.nn as nn import torch.nn.functional as F class DoubleConv(nn.Module): def __init__(self, in_ch, out_ch): super().__init__() self.conv nn.Sequential( nn.Conv2d(in_ch, out_ch, 3, padding1, biasFalse), nn.BatchNorm2d(out_ch), nn.ReLU(inplaceTrue), nn.Conv2d(out_ch, out_ch, 3, padding1, biasFalse), nn.BatchNorm2d(out_ch), nn.ReLU(inplaceTrue), ) def forward(self, x): return self.conv(x) class UNetRF(nn.Module): def __init__(self, in_channels2, n_classes4): super().__init__() self.enc1 DoubleConv(in_channels, 16) self.enc2 DoubleConv(16, 32) self.enc3 DoubleConv(32, 64) self.pool nn.MaxPool2d(2) self.up2 nn.ConvTranspose2d(64, 32, kernel_size2, stride2) self.dec2 DoubleConv(64, 32) self.up1 nn.ConvTranspose2d(32, 16, kernel_size2, stride2) self.dec1 DoubleConv(32, 16) self.out nn.Conv2d(16, n_classes, kernel_size1) def forward(self, x): e1 self.enc1(x) e2 self.enc2(self.pool(e1)) e3 self.enc3(self.pool(e2)) d2 self.dec2(torch.cat([self.up2(e3), e2], dim1)) d1 self.dec1(torch.cat([self.up1(d2), e1], dim1)) return self.out(d1)逻辑说明编码层做了两次2倍池化输入尺寸从H,W缩到H/4,W/4所以输入H和W都需要能被4整除。BatchNorm对RF这种尺度波动的输入很重要它把每层激活归一化减少设备间增益差异带来的分布偏差。代码中的通道数16、32、64是一个保险的起点显存不够时可以把16改成8需要更强表达则改成32。参数说明in_channels2对应IQ双通道如果你只使用包络图把输入stack成1,H,W并把in_channels改成1。n_classes按标签数量调整我这里用的是四分类0背景、1管腔、2内膜中膜复合体、3外膜。如果你暂时只想做二分类就把n_classes改成2自己同步改损失。3.3 损失函数软Dice与Cross Entropy混合颈动脉壁在整帧RF图像里只占很小一块面积直接用交叉熵训练网络很容易把所有像素都预测成背景。混合一个软Dice损失能强制优化前景区域的重叠率。下面是我常用的损失定义class SoftDiceLoss(nn.Module): def __init__(self, n_classes, smooth1e-6): super().__init__() self.n_classes n_classes self.smooth smooth def forward(self, logits, targets): probs F.softmax(logits, dim1) # (B,C,H,W) mask targets.squeeze(1).long() # (B,H,W) one_hot F.one_hot(mask, num_classesself.n_classes) # (B,H,W,C) one_hot one_hot.permute(0, 3, 1, 2).float() inter (probs * one_hot).sum(dim(2, 3)) union probs.sum(dim(2, 3)) one_hot.sum(dim(2, 3)) dice (2 * inter self.smooth) / (union self.smooth) return 1.0 - dice.mean()使用方式ce nn.CrossEntropyLoss() dice_loss SoftDiceLoss(n_classes4) for x, y in dataloader: logits model(x) # (B,4,H,W) y y.long() # (B,H,W) 或 (B,1,H,W) loss ce(logits, y.squeeze(1)) dice_loss(logits, y)逻辑说明软Dice对每个类别独立计算再取平均天然处理类别不平衡。外膜往往只占像素的5%到10%如果没有Dice项小类别的梯度会被大面积背景淹没。混合CrossEntropy是为了让早期训练稳定Dice则在拉高小区域的重叠率。参数说明smooth取1e-6足够CE和Dice的权重各为1.0。如果你发现外膜类别的Dice始终低可以把dice_loss乘以2.0。如果某些batch里某类完全没出现smooth能兜底但最好的做法是在dataloader里按类别比例采样保证每个batch都至少包含两个前景类别。还有一点这里的一one_hot没有包含ignore index如果你的标注里有未标记区域最好在CE里也设ignore_index并给标注-100或255。有了模型和损失下面走一遍完整训练流程。4. 训练一次颈动脉RF分割采样、增强和评估指标4.1 数据采样RF帧大用滑窗剪裁而不是整帧缩放颈动脉RF帧常见深度采样点是2048或3072扫描线192或256。直接resize到256会丢失血管壁的细节尤其是膜层之间的距离可能只有2到3个采样点。常见做法是沿着扫描线方向滑窗裁出256x256的patch。patch的尺寸要能覆盖血管壁全厚度颈总动脉直径约6到8mm如果每个像素0.1mm折算成60到80像素加上两侧外膜和周围组织256x256足够。训练时每个patch可以重叠推理时用重叠滑窗再对输出概率取均值能减少拼接缝。4.2 增强可以翻转和旋转不要加随机高斯噪声RF数据是物理采样高斯噪声增强会给网络制造假的纹理。安全的增强包括水平翻转、小角度旋转、缩放以及增益扰动。增益扰动是RF分割里特别有效的一步对整帧IQ乘一个0.9到1.1的随机因子模拟设备增益档位的小幅变化。import numpy as np from scipy import ndimage def augment_rf_train(iq, mask, angle_range15, gain_range(0.9, 1.1)): # iq: (2, H, W), mask: (H, W) if np.random.rand() 0.5: iq np.flip(iq, axis2) mask np.flip(mask, axis1) # 增益扰动整体缩放保持I/Q相位角不变 gain np.random.uniform(*gain_range) iq iq * gain # 小角度旋转不改变输出尺寸 angle np.random.uniform(-angle_range, angle_range) if np.abs(angle) 0.5: iq ndimage.rotate(iq, angle, axes(1, 2), reshapeFalse, modenearest) mask ndimage.rotate(mask, angle, axes(0, 1), reshapeFalse, order0, modenearest) return iq, mask逻辑说明翻转的axis2对应水平方向注意mask要同步翻。旋转用ndimage.rotateRF图像用线性插值标签用order0保持整数reshapeFalse避免旋转后尺寸变化导致标签错位。增益扰动要在旋转之前做这样旋转插值会在变换后的幅值上产生微小模糊等效于增加一点空间平滑。参数说明angle_range不要超过15度颈动脉在切面里基本水平或微微倾斜大角度旋转会生成不真实的体位。外膜类别的标注面积小旋转有风险如果发现旋转增强后而是伤害可以把它关掉只保留翻转和增益扰动。4.3 训练循环AdamW学习率1e-4混合损失跑50轮下面是一个最小可运行的训练循环。模型定义和损失函数沿用上一节dataloader返回的x形状是B,2,H,Wy形状是B,H,W取值范围为0到3。import torch import torch.nn as nn import torch.optim as optim from torch.utils.data import DataLoader device torch.device(cuda if torch.cuda.is_available() else cpu) model UNetRF(in_channels2, n_classes4).to(device) optimizer optim.AdamW(model.parameters(), lr1e-4, weight_decay1e-5) ce_loss nn.CrossEntropyLoss() dice_loss SoftDiceLoss(n_classes4) for epoch in range(50): model.train() total_loss 0.0 for x, y in dataloader: x, y x.to(device), y.to(device) # y.unsqueeze(1) if shape is (B,H,W) optimizer.zero_grad() logits model(x) loss ce_loss(logits, y.squeeze(1)) dice_loss(logits, y) loss.backward() optimizer.step() total_loss loss.item() avg_loss total_loss / len(dataloader) print(fepoch {epoch:02d} loss {avg_loss:.4f})逻辑说明CrossEntropy要求target是B,H,W的整数张量所以y要squeeze去掉通道维。SoftDiceLoss内部做one_hot需要target是B,H,W且是long类型。AdamW比Adam在超声高动态范围输入下更容易稳定weight_decay取1e-5不需要太大。参数说明学习率1e-4是常用起点。如果loss在10轮内震荡把学习率降到3e-5。50轮对小数据集可能不够我一般结合Early Stopping监视验证集Dice最多跑到100轮。如果训练数据来自多个设备每个batch尽量混合多个设备样本否则模型会按设备切换跳变。4.4 评估指标Dice按类别看HD95也要算训练完不能只看平均Dice。颈动脉分割里外膜和内膜中膜是两类小目标平均Dice高往往是因为背景面积大。下面这段代码分别计算每个类别的Dice和95% Hausdorff距离。import numpy as np from scipy.spatial import cKDTree def dice_one(pred, target, cls): p (pred cls) t (target cls) inter np.logical_and(p, t).sum() union p.sum() t.sum() return 2 * inter / union if union 0 else 0.0 def hd95_one(pred, target, cls): p np.argwhere(pred cls) t np.argwhere(target cls) if len(p) 0 or len(t) 0: return np.inf d1 cKDTree(t).query(p)[0] d2 cKDTree(p).query(t)[0] return np.percentile(np.concatenate([d1, d2]), 95)逻辑说明HD95用点到点的欧氏距离单位是像素要换算成毫米需要乘以深度方向的mm_per_sample。HD95比Dice更能暴露边界偏移Dice可能达到85%但HD95已经到5个像素说明边界位置肉眼可见有偏差。Dice只关心重叠面积HD95关心最大离散程度两者互补。参数说明类别标签中如果某类在图像上不存在hd95返回inf汇总时要过滤掉不要直接mean。颈动脉分割一般要求外膜类别的HD95小于2mm如果超出优先检查坐标映射和标注一致性。训练过程中遇到的很多翻车其实不来自模型而是数据侧。下面集中写几个我踩过的坑。5. 射频分割避坑清单五个让Dice从90掉到70的常见问题5.1 训练集包络清晰验证集却发虚现象模型在训练集的Dice接近0.9到了验证集只有0.7肉眼可见验证集的包络图像对比度极低血管壁边界发灰。原因训练集和验证集来自不同超声设备或不同增益档位RF包络的幅值分布不一致。我之前用全帧min-max归一化训练集里恰好有一个强反射点把上限拉高验证集没有这种点灰度整体被压暗。解决不要用固定的min-max改为每帧用第1和第99百分位数归一化lo np.percentile(log_env, 1) hi np.percentile(log_env, 99) log_env np.clip((log_env - lo) / (hi - lo 1e-8), 0, 1)这个改动牺牲了极亮点的动态范围但换来跨设备稳定性。如果还不行在数据加载时对log_env乘0.9到1.1的随机增益当作训练增强。5.2 标注明明对齐B超图分割图却整体偏移现象训练正常推理时预测的血管壁轮廓比标注整体向下偏移3到5个像素在深度方向尤其明显。原因B超显示图像做了扫描变换或像素插值标注时用的像素坐标与RF深度采样点坐标之间存在非线性映射。我早期直接按像素等比例换算忽略了显示图深度轴的插值比例。解决标注不要画在最终显示图上最好画在包络图上。如果只有B超图先做坐标校准在包络图中找两个强反射点例如血管前后壁的强回声位置在B超图上找到同名位置计算两者深度像素值的比值作为scale再按float尺度换算所有标注。不要把这个scale取整。5.3 换一个增益参数分割结果崩掉现象设备把超声增益从40dB调到55dB后同一患者的外膜区域被预测成背景。原因RF信号幅值随增益整体放大但B超图显示时又做了自动增益补偿直接喂给网络的RF包络图幅度范围漂移网络误把低幅度区域当噪声。解决在预处理中对RF包络做局部对比度归一化比如把每个像素包络值除以它附近邻域均值突出相对对比度而弱化绝对幅度。也可以把输入换成IQ双通道相位信息对增益缩放更鲁棒因为IQ的归一化共用scale已经消掉整体幅值差异。如果快速验证给log_env乘随机增益属于“后悔药”能短期有效但不能根治。5.4 小batch下Dice剧烈波动现象batch size取4时验证Dice在相邻两个epoch之间跳5个百分点学习率降到1e-5也没用。原因RF帧的标注类别不均匀有些patch几乎没有外膜区域SoftDice在这些patch上会给出很高错误梯度。小batch下这种噪声没有被平均掉。解决做一个简单的类别均衡采样。统计每个训练样本里各类别像素比例按“至少包含两类前景”的规则选择patchbatch内尽量保证每种类别都出现。如果仍波动在损失函数里对每个样本单独计算Dice后再求均值避免小batch整体Dice被一个无类别样本主导这属于处理类别不均衡的常见习惯。5.5 显存不够整帧训练不起来现象512x2048的RF帧直接塞进U-Net12GB显存也OOM。原因深度维太长即使两次下采样中间特征图仍然有128x512显存占用高。解决训练时用滑窗裁256x256 patchbatch size设8推理时重叠滑窗取平均。颈动脉壁厚度只有几十像素256x256已经包含足够上下文。如果你想保留整帧可以把模型改成深编码器加同分辨率卷积但没必要。滑窗方案的关键是patch与patch之间要有至少16像素重叠否则边界处容易出现拼接痕迹。避坑说完了最后讲一个特别实用的验证方法。6. 从分割结果测量IMT用厚度验证颈动脉壁分层是否真的可用分割模型输出了四类标签其中第2类是内膜中膜复合体。沿着每条扫描线找到这类标签在深度方向上的起始和结束位置两者之间的距离就是这段血管壁的IMT。颈动脉一条扫描线可能穿过前壁和后壁两个血管壁段所以要先把连续深度段切开再分别计算每段厚度。import numpy as np def imt_from_mask(mask, label_wall2, mm_per_sample0.026): wall_mm [] for col in range(mask.shape[1]): depths np.argwhere(mask[:, col] label_wall).ravel() if len(depths) 0: continue # 将连续深度切分成段 split np.where(np.diff(depths) 1)[0] 1 segments np.split(depths, split) for seg in segments: if len(seg) 2: wall_mm.append((seg[-1] - seg[0]) * mm_per_sample) return np.mean(wall_mm) if wall_mm else float(nan)逻辑说明mask每一列代表一条扫描线沿着深度方向自上而下搜索。正常颈动脉的IMT在0.6到1.0mm之间如果算出来平均大于2mm说明内膜中膜标签和外膜标签粘连或者类别定义串了。如果返回nan说明有的列没有分到wall标签这也是一个错误信号。验证时还有一个技巧连续取几帧超声视频IMT曲线应该是平滑的。如果相邻帧IMT突然跳0.5mm以上大概率不是生理变化而是分割边界抖动。可以用前后3帧的中值滤波把IMT曲线压平同时观察滤波前后差异差异大的样本要重点检查说明那一帧的分割不稳定。把分割图沿扫描线方向提取内膜边界和中膜边界画在包络图上人工核对比只看Dice更能发现系统性偏移。我现在每次做完颈动脉RF分割第一件事不是看平均Dice而是把IMT量一遍。这个数字太诚实了边界只要错一个采样点IMT就从正常值跳到2mm以上。希望这个习惯能帮到你。本文还有配套的精品资源点击获取