胸部数字合成体层成像时的呼吸信号提取

李 仟，刘 佳，莫 英，于佳弘，张昆鹏，张 华

南方医科大学生物医学工程学院，广东 广州 510000

肺癌是最常见的致命癌症之一［1］。胸部数字X射线平片摄影（DR）是肺部病变检测［2］最常用的影像学方式。然而，DR产生的二维图像包含重叠的解剖结构，大大降低了DR对于肺部结节尤其是肺部微小结节或隔膜下结节检测［3-4］的特异性。胸部数字合成体层成像（CTS）是一种新兴的能够通过去除部分重叠组织来提高解剖的可视性的成像技术，该技术利用有限角度下的扫描投影重建伪3D图像［5-7］，因此，CTS检测肺结节的灵敏度至少是DR的3倍［8-9］。肺癌胸部合成的检出率与低剂量CT相当［10］。与标准胸部CT相比，胸部数字合成体层成像分辨率更高，对患者的辐射剂量更低。可作为一种增强肺实质检测和描述肺结节等异常的影像学工具［11］。

呼吸运动引起的伪影在CTS检查中较为常见。在某些情况下，由于呼吸运动的存在，导致CTS对肺部病灶的特异性降至DR水平［9,12-13］。在临床实践中，放射技师很难直接从X线投影图像上判断病人是否屏住呼吸。如果能够在扫描过程中自动识别是否屏气良好（存在呼吸运动），将会有非常重要的临床意义。在CTS成像过程中，X射线管扫描在平行于探测器的平行直线路径上［14］，投影域中的运动物体的运动模式与CT的不同，因此通常用于从胸锥束CT 投影中提取呼吸信号的方法［15］，例如Amsterdam Shroud（AS）［16］方法，可能不适用于CTS。

目前，对CTS中呼吸信号提取算法的研究较少，sin-quadratic模型［17］对于某些特殊状态，如扫描过程中病人存在不规律的呼吸，可能不适合提取呼吸信号。在CTS成像中，管沿患者的头尾方向直线移动，且每个视图的X射线投影都可以清楚地显示横膈膜的边缘线。当前X射线球管的速度有限，故对呼吸信号提取影响不大。在此基础上，本文提出了一个模型，以横膈肌的物理位置作为替代信号来分析胸部数字合成体层成像时的呼吸运动。该模型不仅包含了患者的呼吸运动，而且考虑了CTS系统几何结构引起的基础运动。本文依据投影数据本身，通过对目标对象呼吸运动的实时分析，不需要任何其他外部呼吸检测装置作为参考，直接从投影图中提取呼吸信号。

1 方法

1.1 CTS的运动分析

1.1.1 静止物体CTS运动分析 在CTS扫描过程中，X射线球管与探测器是连续平行相对运动的［18-20］，其运动连续且均匀。当X射线管匀速移动时，单个点的运动轨迹将被很好地限制在一定范围内。如图1所示，D是X射线管平面到探测器所在平面的距离，d是探测器与X射线管同步移动支点平面的高度。假设在某个时刻，X射线管位于坐标c处，探测器的中心位于坐标b处，那么位于（x，y）的物体的投影脉冲相对于探测器中点的相对位移m则有：m=b-a

图1 CTS的几何形状Fig.1 Geometry of chest tomosynthesis.

由CTS系统几何可以得到：

由上文可知，一旦CTS的几何形状固定，m就与c线性相关。对于使用固定探测器的系统（例如GE Volume RAD），m与c也存在线性关系。如果X射线管匀速移动，那么相邻投影图中单个点的移动距离几乎是恒定的［21］，这是CTS模式所固有的。

然而，横膈膜的运动比单个点的运动要复杂的多。从解剖学的角度来看，人体横膈膜肌表现为丘陵状。因此，即使是屏住呼吸的病人，每次曝光时X射线和横膈膜的交点是不同的。式（1）所示的线性关系不足以描述横膈膜的整体的运动轨迹。图2显示了有无呼吸运动的横膈膜的CTS扫描示意图。

1.1.2 CTS运动物体的运动分析 除了CTS系统成像几何引起的基础运动之外，在CTS成像过程中无法屏住呼吸的患者，其横膈膜投影位置也会受到呼吸运动的影响。所以投影X线片中横膈肌的运动可以看作是成像系统引起的基础运动和患者呼吸运动的结合。同时，本文中不再关注横膈肌上单点的运动状态，而是将横膈肌看作整体来构建运动模型。具体来说，我们可以将横膈膜的运动建模为：

其中M表示在投影中横膈膜的总运动，M呼吸表示患者呼吸引起的运动，M基础表示基础运动。

1.2 运动模型

如图2B所示，X线投影图中横膈膜的运动可以概括为由系统几何引起基础运动加上病人呼吸引起的呼吸运动。考虑到如图2所示的，探测器、放射源异界扫描位置的关系，可以推定其线性相关。根据本文2.1分析，采用一次多项式估计横膈膜的由CTS系统几何引起的基础运动分量，并从整体运动中减去基本运动，呼吸信号就可以被检测到。

图2 CTS扫描横膈膜分别在（A）静息状态（B）呼吸状态Fig.2 CTS scan of the diaphragms in stationary state(A)and when breathing(B).

1.3 呼吸信号提取

基于以上分析，我们在图3中给出了提取呼吸信号工作流程图，该流程图包括以下步骤：

步骤一：从投影中提取横膈膜信息（图3A～C）

通过提取膈膜的运动来产生呼吸信号。按照投影顺序，依次对每个投影图像（图像尺寸为1440*1440）进行差分。本文中的差分图像由每个投影图像自身两两相对应像素值相减（同一列中的第n行的像素值减去第n-2行的像素值）得到，以突出图像中的变化部分。这种方法差分后的图像往往能够检测出运动目标的轮廓，得到竖直方向上也就是横膈膜运动方向上的像素值的变化程度。一方面由于横膈膜位于胸腔中，在投影图像的中下部分；另一方面由于左侧横膈膜的运动受心脏跳动的影响；所以在投影图中选择右边的横膈膜运动来提取信号。在实现过程中，我们根据多幅图像分析获得构造横膈膜位置图所需处理的图像区域。再通过将该区域水平方向求和，获取该幅图像横膈膜位置信息。

步骤二：将逐帧横膈膜位置信息组合获得横膈膜位置图（图3D）

横膈膜位置图是将信息列按照其投影顺序组合在一起，得到一个的大小为nb*np的矩阵。其中nb是每列的像素值的数目，np是投影图的数目。

步骤三：从横膈膜位置图中提取横膈膜运动轨迹（图3E）

本文中，我们利用基于动态规划算法从横膈膜位置图中提取横膈膜的运动轨迹。限制条件是轮廓线精确地通过每列仅一次。基于DP理论［22-24］，可以建立成本函数：

其中（x1，x2，...xn）是每列中路径的大小，n是列数（也是本研究中投影的数量）。i数据项c（xi）是通过路径xi的成本，先前项d（xi-1,xi）是xi-1到xi之间部分路径的成本。这里，用横膈膜位置图的梯度大小去构造数据项：。先前项是xi-1与xi之间的绝对距离同时满足d（xi-1，xi）≤k。所以轮廓线也就是横膈膜的运动轨迹是通过求解等式（3）得到的。

步骤四：拟合横膈膜运动（图3F）

对步骤三中获取的膈肌运动轨迹，我们采用一次函数y=ax+b来拟合横膈膜的基本运动分量，拟合方法为最小二乘法。最小二乘法通过等式（4）最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得拟合函数的参数，并使得这些求得的数据与实际数据之间误差的平方和为最小。

步骤五：获取呼吸信号（图3G）

图3 应用一次模型提取呼吸信号流程图Fig.3 Flowchart of applying model for respiratory signal extraction.(A)are the CTS projections.The projections of the diaphragm in the right pulmonary field of the red rectangle in(B)selected to construct the diaphragma location map.The orange circle marked curve in(F)is the diaphragm moving trajectory extracted from the diaphragm location map(E).The solid blue curve in(F)and the black circular marked curve in(G)are the model fitted diaphragm trajectory caused by the system respectively.The circle marked red curve in(G)is the estimated respiratory signal.

呼吸信号可以通过等式（5），利用横膈膜的整体运动减去拟合的一次函数就可以得到。

其中M表示在投影中横膈膜的总运动，M呼吸表示患者呼吸引起的运动，M拟合表示拟合的基础运动。

1.4 数据采集

1.4.1 物理体模数据 物理体模实验在LUNGMAN上进行。该体模是一个精确的真人大小的人体躯干解剖模型，胸壁的厚度是基于临床数据的测量。软组织替代材料和合成骨的X射线吸收率与人体组织非常接近。体模的骨骼和血管在图像上显示出与真人一样的对比度和电子管电压由纵隔、肺血管和腹部块组成的内部成分很容易分离，允许插入模拟肿瘤或其他病变。物体体模数据采集过程与下文中患者数据采集过程相同。图4给出了物理体模的图片。

图4 物理体模Fig.4 Multipurpose chest phantom.

1.4.2 XCAT 体模数据 采用动态数字胸部XCAT 体模［25］进行实验仿真。XCAT体模可以提供真实的人体解剖模型。在人体解剖模四维胸部体模中设计10个呼吸阶段。按照设定的次序进行排列10个呼吸阶段。设定数字体模的体素大小为1 mm×1 mm×1 mm。仿真中使用了在下节中描述的Shimadzu系统的几何结构。我们将大约一个半呼吸周期的4D体模投影到模拟CTS系统上，得到74张尺寸为400×400的投影。采用快速射线追踪算法进行仿真。不考虑X射线散射、光束质量和探测器特性等因素的影响，因为这些因素不会影响所提出模型的适用性。图5为模拟投影X线照片示例。

图5 XCAT数据仿真CTS投影图像Fig.5 Digital XCAT phantom-based simulation.

1.4.3 病人数据 有关病人的研究得到了南方医科大学南方医院伦理委员会的批准，在病人知情的情况下收集并使用数据。病人数据是通过Sonialvision Safire II射线照相系统（Shimadzu,Japan）获取的，该系统配备了一个直接转换的数字平板探测器。震源到等中心点的距离为998 mm，震源到检测器的距离为1098毫米。数据采集时设置管电压为110 kVp，管电流［26］为160 mA，采集管扫角为±20°、尺寸为1440×1440的X线投影图片74张。采集3组患者数据，分别标记为患者a,b,c。其中患者b为明显的脊柱侧弯病人。图6显示了3组患者数据的X线投影图。

图6 3个病人的X线投影图Fig.6 Projection radiographs of 3 patients.A:Patient a.B:Patient b.C:Patient c.

2 结果

2.1 物理体模数据的结果

从图7中可以看出，横隔膜的运动轨迹是线性的。

图7 物理体模数据结果Fig.7 Results of multipurpose chest phantom data.A:Projections of multipurpose chest phantom.B:Extracted diaphragm trajectory in stationary state.C:The trajectory of the diaphragm in stationary state is linear.

2.2 XCAT数字体模数据的结果

图8显示了XCAT数字体模横隔膜位置图和提取的横隔膜运动轨迹。从结果可以看出，横膈膜的边缘可以准确的代表运动。基于动态规划的技术可以从横膈膜位置图很好的提取横膈膜运动轨迹信息。图9C红色带圈标记曲线为本文模型提取出的呼吸运动轨迹。图10为设计的呼吸信号（圆形标记红色实线）与本文方法提取的呼吸信号（圆形标记黑色实线）以及陶等人方法提取的呼吸信号（‘+’标记绿色实线）三者之间的视觉对比。通过计算得到陶等人提出的方法提取的呼吸信号曲线与设计的呼吸信号曲线之间的相关系数为0.9683，而本文方法提取的呼吸信号曲线与设计的呼吸信号曲线之间的相关系数为0.9797。

图8 数字体模数据的膜片位置图和提取的膜片运动轨迹（红色曲线）Fig.8 Diaphragm location map of the phantom data and the extracted diaphragm motion trajectory(red curve).

图9 数字体模数据结果Fig.9 Results of digital phantom data.A:Diaphragm motion trajectory.B:Model fitted trajectory (blue curve).C:Diaphragm trajectory(circle-marked black curve)and respiratory motion trajectory(circle-marked red curve).

图10 本文方法提取的呼吸信号曲线Fig.10 Comparison of the respiratory signal curve extracted by the method in this study (black curve marked by diamond),the respiratory signal curve extracted by Tao et al.'s method (green curve marked by '+') and the respiratory signal curve designed by volume model(red curve marked by circle).

2.3 患者数据结果

3例病患利用本文所提出的一次模型提取呼吸信号的结果如图11～12所示。可以看出，三位患者都存在不同程度的呼吸。患者a和患者b都是有规律的呼吸，而患者c扫描前期屏息良好，后期突然剧烈呼吸。

图11 3位患者数据利用一次模型拟合基本运动分量结果Fig.11 Model fitting of basic motion components based on data from 3 patients.A:Patient a.B:Patient b.C:Patient c.

3 讨论

本文在充分考虑到CTS系统的扫描几何形状引起的基础运动以及患者自身的呼吸运动后，提出了一种新的呼吸运动分析方法。在实验中，用一次多项式模拟了由CTS系统扫描几何形状引起的线性基础运动。由于在实际应用过程中病人呼吸状态的不确定性，我们并未用某一种特定的模型去拟合呼吸运动。呼吸运动信息是由实际的横膈膜的整体运动信号减去拟合的横膈膜基础运动得到。

具体地说，呼吸信号是通过从每个二维投影中提取横膈膜运动信息估计得到的。这是通过将投影图进行差分后手动选择大致位于右隔膜中央区域的像素水平求和得到信息列，作为替代呼吸信号来完成的。因为在这个位置的运动比其他位置更明显。提取的信息列按照投影顺序一起排列，以产生一个隔膜位置图。采用基于动态规划边缘提取方法提取光阑运动轨迹。然后利用一次模型拟合运动轨迹，估计横膈膜由系统几何引起的基础运动分量，并从整体运动中减去基础运动，可以很容易地估计呼吸信号。整个过程的计算时间远小于1 s，因此本文方法可以看作是一种实时的方法。

本研究中，物理体模数据提取结果可以看出，利用物理体模无呼吸运动的情况下提取横膈膜运动轨迹来测定用一次线性模型标定基础运动的可行性，本文模型能够很好地估计横膈膜总体的由系统几何引起的线性基础运动分量。定量分析表明，XCAT仿真实验中本文提出的方法和2018年Tao等［17］提出的sin-quadratic模型都能有效地提取设计的呼吸信号，本文模型所提取的呼吸信号与原设计的呼吸信号之间的相关系数为0.9797，陶等人提出的模型所提取的呼吸信号与原设计的呼吸信号之间的相关系数为0.9683，与Tao等［17］的方法相比本文所提出的模型精确度更高。定性和定量结果均可表明，该模型能有效地从CTS投影中提取呼吸信号。3组患者数据结果也表明，在不同的呼吸状态下，本文模型都可以有效地对呼吸信号进行估计，表现得较为稳健。在病人有规律的呼吸模式下，本文提出的方法和2018年Tao等［17］提出的sin-quadratic模型都能有效地估计患者的呼吸信号。但有一点值得注意：对于不规则呼吸模式的患者，相比Tao等［17］利用正弦函数去拟合呼吸运动，本文并未用某一种特定的模型去拟合呼吸运动。

图12 3位患者数据提取呼吸信号结果Fig.12 Results of respiratory signal extraction from the datasets of 3 patients.A:Patient a.B:Patient b.C:Patient c.

因为CTS成像的扫描时间为6～12 s，也就是大约有1～1.5个呼吸周期。在真实的临床病例中，部分患者在整个数据采集过程中可能无法屏住呼吸如同患者数据中的患者c，在采集过程中突然剧烈呼吸，很明显这种不规则的呼吸状态并不能用正弦函数来表达［27-28］。在这种情况下，Tao等［17］提出的sin-quadratic模型表现的不够稳健不适合提取呼吸信号。然而本文并未用某一种特定的模型去拟合呼吸运动，不规则呼吸信号仍然可以检测到，本文模型表现得更为稳健，适合大多数真实的临床病例扫描状态下。

总的来说，我们开发了一种实时的呼吸信号运动分析方法。放射技师很难直接从X线投影图像上判断病人是否屏住呼吸，本文的方法可以在扫描过程中自动识别是否屏气良好（存在呼吸运动），直接从CTS投影中提取横膈肌运动轨迹，而不需要任何额外的呼吸检测设备或其他参考图像，具有重要的临床意义。据所获得的呼吸信号，医生可立即对患者进行重新扫描或指导无运动CTS图像重建。一旦检测到强烈紊乱的呼吸信号，医生可立即对存在呼吸运动的病人进行重新扫描。在呼吸运动存在的情况下模糊的肋骨和脊柱在无呼吸运动4D-CTS重建［29-30］中可以清晰地呈现出来，提高医生的诊断精确度。