一种时间强度曲线的呼吸运动补偿及双峰拟合的方法
【技术领域】
[0001] 本发明涉及超声造影成像技术领域，特别涉及一种时间强度曲线的呼吸运动补偿 及双峰拟合的方法。
【背景技术】
[0002] 超声造影剂在临床B超成像中的引入，有效增强血流回波信号的强度，极大的提 高超声图像中灌注组织的对比度，弥补了原有普通超声检测信号的不足，并使利用超声造 影技术来定量估计血流灌注成为可能。在人体微小血管和组织灌注检测与成像方面，与CT、 MRI等其他医学检测方法相比，具有成像效果好、实时、操作简便、无离子辐射、无损性、使用 面广等优点，因此在医学超声成像中发挥着越来越重要的作用。在具体操作中，医生通常利 用时间强度曲线（time-intensity curve, TIC)来反映局部组织血液中造影剂含量或浓度 随时间的变化，来定量分析局部组织的血流灌注信息。TIC的获取通过采集造影剂在组织血 管灌注过程中的图像帧序列，选定感兴趣区（region of interesting, ROI)，测量造影后各 帧图像感兴趣区的声学强度，绘制声学强度随时间变化曲线而获得的。TIC可反映微泡背 向散射强度随时间的变化，即局部造影剂浓度随时间的变化。目前，医生通常通过提取TIC 的曲线参数，如峰值、峰值时间、峰面积、冲入率等来进行病理评价，这种操作技术已在飞利 浦HD15造影定量分析软件中出现；并已广泛用于估计心肌血流灌注、区分肝脏的良性恶性 病变、研宄肿瘤的新生血管的血液供给情况等。
[0003] 然而在进行超声造影诊断时，病人不自主的呼吸将使待观察区域发生运动偏移， 这在肝脏造影成像中，尤为明显。而对于TIC的提取则基于某一 ROI的声学强度，一般此ROI 在序列B超造影图像中的位置是固定不变的，这样，人体的呼吸运动会使前后帧ROI区域不 是同一位置，而发生了偏移。这对于医生进行定量分析及处理时造成了干扰，并导致结果出 现错误，影响超声造影定量评价的准确性。在实际操作中，医生虽要求病人屏住呼吸，减少 呼吸运动的干扰，但在长达数分钟的造影检查过程中病人并不能持久屏住呼吸。所以，一般 情况下，超声造影成像是在病人自由呼吸过程中进行的，但呼吸运动的干扰，降低了 TIC提 取的准确性。目前在文献中用到的呼吸运动补偿方法包括图像配准方法、呼吸门控算法等。 图像配准方法，如块匹配方法，算法简单、易于实现，但其遍历搜索方法耗时太久，所选用的 平均误差函数等作为评判标准并不全面。而呼吸门控算法中目前最新文献中以主成分分析 方法为代表，但其运算复杂度较大，不易实现；同时需人工剔除错误的肝脏造影图像，不仅 耗时，而且存在主观因素的影响。
[0004] 另外，目前文献、专利中所用到的基于短时团式灌注模型分析TIC时，关注点主要 是造影剂的首次通过峰拟合曲线及参数的提取，而忽略了因心脏周期搏动带来的造影剂在 观察区域再循环峰的客观存在。而再循环峰包含了丰富的生理信息。在微泡的浓度达到峰 值浓度之后，随后发生的是相对缓慢的微泡浓度下降，直至下降趋势受到微泡再循环的干 扰，微泡再循环尖峰表明微泡已经完成整个循环并重新回到成像目标部位。首次峰和再循 环峰之间的拐点具有重要的生理学意义。该拐点是微泡冲出相与再次冲入相的冲入冲出量 相等的时间点，该点的时刻可以作为判断微循环及局部血供性能的重要依据。因此，在对 TIC进行曲线拟合时，要同时关注首次通过峰和再循环峰的拟合，即双峰拟合，只有这样才 能充分提取TIC曲线所包含的生理学信息，充分利用TIC曲线来进行病理评价。
[0005] 为此，对TIC进行运动补偿，并进行双峰拟合，对于准确提取TIC趋势、获取准确的 灌注信息以及更好地对血流灌注进行评价具有重要意义。
【发明内容】
[0006] 针对肝脏超声造影成像使用TIC进行灌注评价时TIC受到强烈的呼吸运动干扰以 及再次循环峰被忽略的实际问题，本发明的目的在于提供一种时间强度曲线的呼吸运动补 偿及双峰拟合的方法。
[0007] 为实现上述目的，本发明采用了以下技术方案：
[0008] (1)基于以累积绝对误差为标准的块匹配方法，在轴向以及横向对提取自序列B 超造影图像的时间强度曲线进行呼吸运动补偿；
[0009] (2)对呼吸运动补偿后的时间强度曲线进行预处理得到预处理后的时间强度曲 线，预处理依次包括频域滤波、时域滤波以及基线归零；
[0010] (3)通过强度阈值法和局部最小值法确认拟合区间的起始点和终止点；
[0011] (4)对预处理后的时间强度曲线进行再循环峰拟合判定，如存在再循环峰，进行再 循环峰拟合，然后输出再循环峰拟合与首次通过峰拟合结果相加得到的双峰拟合结果，并 输出首次通过峰与再循环峰的拐点强度值和时刻，拐点是指首次通过峰下降支与再循环峰 上升支的交叉点。
[0012] 所述步骤（1)具体包括以下步骤：
[0013] (1. 1)设置当前帧、当前块的大小以及搜索路径及搜索区域；
[0014] (1.2)以下一帧为匹配帧，当前块所在的中心位置为初始位移，在搜索区域内依次 按轴向、横向选定匹配块；
[0015] (1. 3)计算匹配块与当前块的累积绝对误差，找到最小累积绝对误差所对应的匹 配块的中心位置，计算该匹配块相对当前块的轴向偏移位移以及横向偏移位移；
[0016] (1. 4)重复步骤（1. 2)~（1. 3)，得到序列B超造影图像中最小累积绝对误差所对 应的匹配块的轴向偏移位移以及横向偏移位移；
[0017] 按帧依次将需提取时间强度曲线的感兴趣区的行、列位置与所述块匹配方法得到 的横向偏移位移、轴向偏移位移对应叠加，实现运动补偿；按帧计算运动补偿后的对应感兴 趣区内的灰度均值以及所述当前帧感兴趣区内的灰度均值，得到呼吸运动补偿后的时间强 度曲线。
[0018] 所述步骤（2)具体包括以下步骤：
[0019] (2. 1)对呼吸运动补偿后的时间强度曲线进行巴特沃斯低通滤波，截止频率为 0. 2~0. 7Hz，滤波器的阶数为2~3 ;然后进行中值时域滤波，邻域窗长取值为：时间强度 曲线总数据点数/20 ;
[0020] (2. 2)对于经过步骤（2. 1)处理后的呼吸运动补偿后的时间强度曲线，选取造影 微泡进入成像平面的时刻之前的时间强度曲线幅值均值作为基线值，将基线值从整个时间 强度曲线序列信号中减去，然后将时间强度曲线中所有的负强度值的点赋值为极小的数。
[0021] 所述步骤（3)具体包括以下步骤：
[0022] (3. 1)找出时间强度曲线的峰值Ip和峰值时间t p;沿着时间强度曲线的上升阶段 按时间增加的方向，找到时间强度曲线强度首次高于IpIO~15%的数据点对应的时刻，作 为拟合区间的起始点；
[0023] (3. 2)设定强度阈值为Ip的40~50%;从tp开始，沿着时间强度曲线下降阶段找 到所有局部极小值及对应的时刻，保存为局部极小值集合；若局部极小值集合为空，则在时 间强度曲线下降阶段按时间增加的方向找到强度首次低于I P20~25%的数据点对应的时 亥IJ，作为拟合区间的终止点；
[0024] (3. 3)若局部极小值集合不为空，从距离tp最近的局部极小值对应的时刻开始，沿 着时间增加的方向进行判定，如该局部极小值小于所述强度阈值，则将该局部极小值对应 时刻作为拟合区间的终止点，否则循环至下一个局部极小值点进行判定；如遍历t p后所有 局部极小值对应的时刻，找不到满足判定条件的局部极小值点，则在时间强度曲线下降阶 段按时间增加的方向找到强度首次低于I P20~25 %的数据点对应的时刻，作为拟合区间的 终止点。
[0025] 所述步骤（4)具体包括以下步骤：
[0026] (4. 1)对预处理后的时间强度曲线进行首次通过峰拟合：对预处理后的时间强度 曲线，确定拟合区间起始点与终止点，并进行线性GAMM团注模型拟合；GAMM团注模型的 尺度因子以及曲线形状参数通过多元线性回归分析方法确定；
[0027] (4. 2)将预处理后的时间强度曲线与首次通过峰拟合结果相减得剩余时间强度曲 线，若剩余时间强度曲线的峰值大于预处理后的时间强度曲线峰值的20~25%，则进行再 循环峰拟合，然后输出双峰拟合结果，否则直接输出首次通过峰拟合结果。
[0028] 若剩余时间强度曲线的峰值大于预处理后的时间强度曲线峰值的20~25%，则：
[0029] (1)对剩余时间强度曲线进行再循环峰拟合：对剩余时间强度曲线，确认拟合区 间的起始点与终止点，然后进行线性GAMM团注模型拟合；GA