基于4D CT的肺部运动呼吸模型

梁 睿，张建勋，冯昌利，赵汝琛

(1. 南开大学机器人与信息自动化研究所，天津 300071；2. 广州雅敦微创科技有限公司，广州 510110)

基于4D CT的肺部运动呼吸模型

梁 睿1，张建勋1，冯昌利1，赵汝琛2

(1. 南开大学机器人与信息自动化研究所，天津 300071；2. 广州雅敦微创科技有限公司，广州 510110)

肺部的呼吸运动会在图像导航穿刺手术中引入误差，并且该误差是穿刺导航手术的主要误差．为解决该问题，提出了一个接近亚体素级精度、方便易用、运算简便的呼吸模型．首先根据肺内部各标记点的真实运动信息推算运动状态参数；然后根据运动状态参数的运动趋势，提取各点运动的总体规律；最后结合已获得的样本数据，得到一个一般化的呼吸模型．交叉实验结果表明，运用该模型可以将由呼吸引起的原始误差降低 75%左右，在呼吸幅度较小的情况下，模型可以达到亚体素级的精度，并且具有运算简便、运用简单的特点．

呼吸模型；曲线拟合；运动状态参数；误差

随着手术导航技术的发展成熟，手术导航定位越来越精准．但是，与此同时，人们对手术导航的精度要求也越来越高．

手术导航系统的误差主要来自于两方面：一方面来源于图像与病体组织的配准误差，主要原因是病体标记点图像扫描产生的伪影、位置传感器的定位精度不高，目前已经有很多非常成熟的伪影消除方法[1]、补偿方法[2]等解决上述问题；另一方面来自于病体自然的生理运动，尤其是人体的呼吸．呼吸是人体内部器官最大的生理运动，根据文献记载，其运动幅度一般有0.5～4.0,cm[3]，对于4.0,cm以下的病灶，呼吸的研究显得尤为重要．呼吸一刻不能停止，它引起的肺部移动量也最大，连带着其他脏器做较大运动，所以研究呼吸时肺部产生的三维变化和位移，建立相应的运动模型预测其运动趋势，对手术导航系统精确度的改进具有重大意义．

此外，CT扫描技术也日益发展，逐渐由 2D CT发展为3D CT，进而出现了4D CT．4D CT能同时监控呼吸信号，为人体呼吸模型的研究开辟了一个新途径[4]．

以往呼吸模型的研究根据指导思想的不同分为4个类别：基于生物学的模型、基于信号的模型、基于图像配准的模型[5]和基于统计学的模型[6]．

基于生物学的模型是根据生理学和解剖学的理论建立力学模型，大都利用有限元分析的方法．但是该方法的假设过于理想化，无法真正诠释肺的运动，实验结果的准确度不高．

基于信号的模型，其外贴标记点无法测知内部运动，而内嵌标记点由于个数有限，只能感知局部个别点．该方法虽然比生物学模型准确，但是研究对象有限的特点限制了该模型的研究进展．

基于图像配准的模型是近年来的研究热点．但是，图像配准尤其是 3D图像配准耗时长，通常达到20～30,h，并且其配准误差不可知．利用带有误差的数据进行研究，容易产生累计误差．

基于统计学的模型是提取若干个样本的呼吸运动信息，再归纳出一个普遍的规律，因此，该方法普遍应用性强，但没有真正表达其个性特征．根据个性特征调整呼吸模型以适应每一个人的呼吸运动，是目前研究发展的方向．

本文提出了一种基于4D CT的肺部运动呼吸模型，结合基于图像配准模型和基于统计学模型的思想，同时避免了基于图像配准模型由图像配准所带来的缺点．首先对4个样本的4D CT的图像数据进行处理，为了避免配准带来的误差，利用医学专家针对整个肺部特征做的标记点进行运动信息分析；然后获得每一个呼吸相位各标记点的位置信息，归纳每一个呼吸相位的运动特征；最后，综合 4个样本的肺部运动特征，拟合出一个普遍的统计学模型．

1 呼吸模型

1.1 研究材料

本文所用第 3方数据以及信息来自文献[7]．该论文及其网站上提供了研究图像配准算法时所用的实验材料、实验过程以及实验结果．笔者利用其中 4个样本的4D CT图像以及医学专家针对这4个样本所做的特征标记点，其中3个样本有100个标记点，1个样本有 41个标记点(标记点位置分布如图 1所示)，每个样本在整个呼吸周期的每个相位均做了标记点[8]．4个样本均患有非小细胞肺癌，4D CT扫描采用螺旋模式，分辨率为 512×512×141，所有切片间距均为 2,mm，每张切片两个方向的像素间距为0.976,6,mm．呼吸信号侦测利用 Lafayette仪器公司生产的气动胸部压力带进行．标记过程采用了一个半自动化的软件，该软件具有良好的交互界面，能够自动计算点与点的相似度，协助观察者找到对应点的位置．

图1 41个标记点空间分布投影Fig.1 Projection of 41 fiducials’ space distribution

1.2 研究方法

由于呼吸运动的速度和幅度具有个体内和个体间的差异性，病患的每个呼吸周期都没有重复性，不同病患的呼吸运动都不同，因此无法研究时间与呼吸运动的关系，但是多项研究表明，呼吸周期的相位与呼吸运动之间有一定的规律[9]．目前确定呼吸相位的方法有很多[10]，其中应用比较普遍的有两种：一种是根据呼吸通气设备或者腹部压力传感器得到信号的振幅来对应呼吸相位；另一种是根据这些设备得到信号的相位来对应呼吸相位．针对不同呼吸相位对 CT图像的肺部呼吸进行建模．临床应用时，可以一方面利用相关设备监测呼吸相位，另一方面运用模型对图像进行调整，使其实时对应该病患体内的肺部运动状态．

该模型的建立基于第1.1节的CT图像．每一个样本的4D CT图像包含10组，分别是呼吸相位0、10%、20%、30%、40%、50%、60%、70%、80%和90%．其中，按照4D CT的惯例标准，0是吸气末端，50%或者 60%是呼气末端．为了更好地找出肺部每个点的普遍运动规律，根据肺部运动规律，假设从 0到50%(或 60%)的变化是肺部运动的最大幅度，并且肺部的每一个点都是从0到50%的位置上做往复运动，那么该点某个时刻的位置可以描述为

式中：t为相位；k为标记点；i为运动状态参数，表示呼吸运动幅度的百分比，取值为 0～100；为 t相位下、i运动幅度时，第 k个标记点的理论估计坐标值；为0相位下第k个标记点的坐标；为50%相位下第k个标记点的坐标；ti为t相位下各点的运动幅度．

根据式(1)需要确定 t与 i的对应关系．本文采用误差最小化法找到这个对应关系，即 t相位下，找到一个i值令误差最小，且有

式中,tkP 为 t相位下第 k个标记点的坐标值．该误差的标准差为

呼吸模型研究过程流程如图2所示．

图2 呼吸模型研究过程流程Fig.2 Flow chart of breathing model research process

1.3 呼吸模型的研究

根据上述研究过程，针对4个样本的材料分别进行研究．每一个相位都会得到一个误差最小时的运动状态参数，称为最优运动参数．得出 4个样本的整个呼吸周期中的最优运动参数的散点图，如图 3所示.

图3 4个样本的最优运动参数散点图Fig.3 Plot of optimal motion parameters for the four samples

2 多样本一般化的呼吸模型

根据以上4个样本的各相位最优运动参数，分别采用不同方法进行曲线拟合，对拟合结果进行比较，得出最优的一般化呼吸运动趋势．选取几种效果较好的结果进行对比，如图4所示．其拟合方法有傅里叶曲线拟合、高斯曲线拟合、多项式曲线拟合和正弦曲线拟合．

傅里叶曲线拟合式为

拟合结果为

高斯曲线拟合式为

拟合结果为

多项式曲线拟合式为

拟合结果为

正弦曲线拟合式为

拟合结果为

何为优质？韩光曙分析称，医院是救死扶伤的场所，必须尊重生命，对生命有敬畏感。“乞丐与王子平等”的济世情怀在医院传承百余年。“但若仅强调平等，机械地完成工作也可以。鼓楼医院希望的是，在提供治病救人服务的同时，让患者感到温馨。”

图4 4种方法的拟合结果Fig.4 Fitting results with four kinds of methods

运用不同函数进行曲线拟合，得出结果如表1所示，误差单位为该点运动幅度的百分比．其中SSE为和方差，RMSE为均方根，则有

SSE、RMSE越小，说明模型选择和拟合越好．

R-square为确定系数，其计算式为

R-square越接近1，方程变量对y的解释能力越强．

R-square-a为调整后的确定系数，它是嵌套模型的最优参考指标，其计算式为

表1 4种方法结果比较Tab.1 Comparison of fitting results with four kinds of methods %

经过比较发现高斯曲线拟合和正弦曲线拟合产生了局部畸变，这种运动状态在肺的自然生理运动中是不存在的，所以只在傅里叶曲线拟合和多项式曲线拟合的方法中进行取舍．从表1中，可以看到傅里叶曲线拟合算法的精度较高，故选择傅里叶拟合结果来进行验证．

3 实验及模型误差

由于4D CT辐射很大，国内该类型仪器很少，而且实验成本高昂，样本取材很受限制，所以采取样本间的交叉验证方法对仅有的 4组样本数据来进行如下4组数值实验，分别计算了其在模型作用下的各相位各点的平均误差、标准差以及分解到x、y、z方向上的平均误差．

该交叉验证实验结果如表2和表3所示．

表2 交叉验证实验结果1Tab.2 Result 1 of the cross validation test

表3 交叉验证实验1的xyz方向上的误差Tab.3 Comparison of the errors in xyz directions in result 1

计算的整个呼吸运动中各点平均运动位移为7.670,8,mm，标准差为 5.054,7，所以，根据模型计算的结果误差减小为无模型计算结果的 16.87%．从表3可以看出，该实验结果误差达到了亚体素级．

3.2 用第1、2、4个样本拟合，用第3个样本验证

该交叉验证的实验结果如表4和表5所示．

表4 交叉验证实验结果2Tab.4 Result 2 of the cross validation test

表5 交叉验证实验2的xyz方向上的误差Tab.5 Comparison of the errors in xyz directions in result 2

计算的整个呼吸运动各点平均运动位移为14.044,5,mm，标准差为 7.204,3，所以，根据模型计算的结果误差减小为无模型计算结果的 26.47%．该组数据的呼吸运动幅度较大，所以误差较大．

3.3 用第1、3、4个样本拟合，用第2个样本验证

该交叉验证的实验结果如表6和表7所示．

表6 交叉验证结果3Tab.6 Result 3 of the cross validation test

表7 交叉验证实验3的xyz方向上的误差Tab.7 Comparison of the errors in xyz directions in result 3

计算的整个呼吸运动各点平均运动位移为5.904,8,mm，标准差为 2.732,7，所以，根据模型计算的结果误差减小为无模型计算结果的 22.87%．从表7可以看出，该组实验的误差也几乎达到了亚体素级别．

3.4 用第2、3、4个样本拟合，用第1个样本验证

该交叉验证的实验结果如表8和表9所示．

表8 交叉验证结果4Tab.8 Result 4 of the cross validation test

表9 交叉验证实验4的xyz方向上的误差Tab.9 Comparison of the errors in xyz directions in result 4

计算的整个呼吸运动各点平均运动位移为5.747,4,mm，标准差为 2.773,9，所以，根据模型计算的结果误差减小为无模型计算结果的 25.81%．从表9中可以看出，该组实验的误差几乎达到了亚体素级别．

综合 4个实验可以看出，运用该呼吸模型进行估计可以把呼吸运动引起的误差降低为 25%左右，而且在呼吸运动较小的情况下可以达到亚体素别．

3.5 讨 论

根据表10的无模型下各相位的平均误差可以看出，本文模型虽然在整体的误差消除作用比较显著，达到了 25%左右，但在首相位(10%)和末相位(90%)时，误差消除作用比较微弱，甚至偶有增长．这是因为，呼吸运动在运动开始和末端时候的位移相对于整个呼吸过程来讲比较小，而且运动增长速率也比较小，这时的呼吸运动存在迟滞效应，运动模式也比较混沌[11]．该阶段下肺脏各点运动具有非线性不规则的特点，每个点运动差异较大，这种现象的分析和建模非常复杂．而本文为了增强临床的实用性和运算的简便性，对该现象的建模进行了统一简化，使这两个相位下模型的运动参数几乎为 0，也就是弱化了此时模型的作用，从而使模型对于这两个相位下的误差消除作用比较小．但是由于该相位下各点的运动幅度相对较小，在正常的呼吸幅度较小的情况下，呼吸误差已经达到了亚体素级别，在图像的呈现上看不出具体差别．即使该相位下对呼吸误差消除作用较弱，也保证了较小的误差，使其在呼吸幅度较小的情况下也达到了亚体素级别．所以此相位下的误差消除作用较弱并不影响整体的误差消除效果．这样既保证了文中模型的有效性，也保证了其实用性和简便性．

表10 无模型下每个呼吸相位各点的平均误差Tab.10 Average movement of each point in each phase without the model

4 结 语

为了处理穿刺手术中呼吸运动引入的误差，本文根据已有的肺部呼吸数据，找到普遍呼吸运动规律，简化为一个运算简捷、方便易用的呼吸模型．首先根据肺部标记点运动数据进行参数寻优，然后针对参数进行多项式拟合，最后根据多个样本得出一般化的呼吸模型．

该模型很容易运用到穿刺手术场景，而且运算量很小，不会影响实时性．在呼吸运动较小的情况下，该模型的误差可以达到亚体素级别．

该模型也有其局限性．该模型假设同一个呼吸相位时每个点的运动参数相同，这只是比较一般的规律，如果能够根据位置的不同具有不同的运动参数变化规律，就会使精度更好，这需要进一步进行研究，同时需要增加样本参数进行验证．

［1］ Bal M，Spies L. Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering[J]. Medical Physics，2006，33(8)：2852-2859.

［2］ Yu Hengyong，Zeng Kai，Bharkhada D K，et al. A segmentation-based method for metal artifact reduction[J]. Academic Radiology，2007，14(4)：495-504.

［3］ Li Xing，Thorndyke B，Schreibmann E，et al. Overview of image-guided radiation therapy[J]. Medical Dosimetry，2006，31(2)：91-112.

［4］ Keall P，Yamamoto T，Suh Y. Introduction to 4D Motion Modeling and 4D Radiotherapy[M]. Berlin Heidelberg：Springer Berlin Heidelberg，2013.

［5］ McClelland J R，Blackall J M，Tarte S，et al. A continuous 4D motion model from multiple respiratory cycles for use in lung radiotherapy[J]. Medical Physics，2006，33(9)：3348-3358.

［6］ George R，Vedam S S，Chung T D，et al. The application of the sinusoidal model to lung cancer patient respiratory motion[J]. Medical Physics，2005，32(9)：2850-2861.

［7］ Vandemeulebroucke J. Lung motion modelling and estimation for image guided radiation therapy[EB/OL]. http：//www. creatis. insa-lyon. fr/rio/popi-model，2012-12-20.

［8］ Murphy K，van Ginneken B，Klein S，et al. Semiautomatic construction of reference standards for evaluation of image registration[J]. Medical Image Analysis，2011，15(1)：71-84.

［9］ Yamamoto T，Kabus S，Lorenz C，et al. 4D CT lung ventilation images are affected by the 4D CT sorting method[J]. Medical Physics，2013，40(10)：101907.

［10］ Heinrich H P，Jenkinson M，Brady S，et al. MRF-based deformable registration and ventilation estimation of lung CT[J]. IEEE Transactions on Medical Imaging， 2013，32(7)：1239-1248.

［11］ Rahni A A A，Lewis E，Wells K. Characterisation of respiratory motion extracted from 4D MRI[C]//SPIE Medical Imaging International Society for Optics and Photonics. Florida，USA，2013：86692Z.

(责任编辑：金顺爱)

Model of Lung Respiratory Movement Based on 4D CT

Liang Rui1，Zhang Jianxun1，Feng Changli1，Zhao Roger2
(1. Institute of Robotics & Automatic Information System，Nankai University，Tianjin 300071，China；2. Acton Medical Device Corporation of Guangzhou，Guangzhou 510110，China)

Lung respiratory movement can cause errors in the operation of image navigation surgery and they are the main errors in the navigation system. To solve this problem，a convenient，simple and easy-to-use breathing model whose precision was close to the sub-voxel was proposed. Motion parameters were first calculated according to the actual lung movement information of each point. The general movement law of each point was then obtained from the movement trend of the motion parameters. Finally a generic lung breathing model was constructed based on all the sample data. Results of cross experiments show that the proposed model can reduce the original errors caused by lung breath by about 75%. Besides，in circumstances of lung breath with small amplitude，the proposed model is simple and easy-to-use and can reach the sub-voxel precision.

breathing model；curve fitting；motion parameter；error

TP391.4

A

0493-2137(2015)04-0298-07

10.11784/tdxbz201311065

2013-11-22；

2014-03-10.

国家科技支撑计划资助项目(2012BAI14B00).

梁 睿（1987— ），女，博士研究生，lrcalr@163.com.

张建勋，zhjx@robot.nankai.edu.cn.

时间：2014-03-13.

http：//www.cnki.net/kcms/doi/10.11784/z201311065.html.