基于PCA和TSSM模型的膈肌预测系统及方法与流程

本发明属于图像处理技术领域，更进一步涉及医学影像辅助干预技术领域中的一种基于主成分分析pca(principalcomponentsanalysis)和分步子空间映射tssm(two-stepsubspacemapping)模型的膈肌预测系统及方法。本发明在不使用标记物的情况下，通过体外胸腹部表面的位移预测人体膈肌的运动。

背景技术：

动态目标跟踪法是目前最为有效的呼气管理技术，直接方式是植入体内的标志物进行成像，进而进行目标跟踪。但实验对象需要进行昂贵的、创伤性手术。植入标志物的替代方法，是使用呼吸信息对内部解剖结构的运动进行建模。通过肺活量测定和实时位置跟踪系统，研究肺肿瘤运动与呼吸运动的相关性。相比于其它解剖结构，膈肌与呼吸运动具有更强的相关性。另一方面，人体内部结构的运动情况，可以通过机体外部来进行观察和预测。

在malinowski等人在其发表的论文"mitigatingerrorsinexternalrespiratorysurrogate-basedmodelsoftumorposition"(jinternationaljournalofradiationoncology,biology,physics,"vol.82,no.5,2012)中公开公布了一种减轻外呼吸替代肿瘤位置模型的误差的方法。该方法的主要实现步骤是：在患者身上贴附了多个标志物，并建立了标志物-肿瘤的位移关系模型。该方法通过评估肿瘤位置预测模型的潜在误差来源，最终选择了偏最小二乘回归(pls)进行模型拟合。该方法存在的不足之处是，模型的准确性与患者数量显著相关，采用外标志物建立的肿瘤位置预测模型的准确性受患者间差异、肿瘤-替代物相关性、训练数据选择的影响。

在南方医科大学申请的专利“一种胸部数字合成x射线体层成像中呼吸运动分析方法”(公开号：cn109146842a，申请号：cn201810745977.4，申请日：2018年07月09日)中公开了一种胸部数字合成x射线体层成像中呼吸运动分析方法，其制备方法包括以下步骤：s1，通过胸部数字合成体层成像，获取cts投影数据；s2,从s1中获取的cts投影数据提取膈肌的运动轨迹；s3,从s2获取的膈肌运动进行进行正弦二次模型拟合。该方法存在的不足之处是：部分患者在整个数据采集过程中可能无法屏住呼吸，在采集结束时开始呼吸。且质心仅以最大振幅的点来表示，忽略了膈肌在左右方向上的位移，导致振幅最大的点正下方振幅最小时的点并不是同一个点从而使得简化对象不准确的问题。

技术实现要素：

本发明的目的是针对上述已有技术的不足，提出了一种基于pca和tssm模型的膈肌预测系统及方法，用于解决pls模型外标志物建立的肿瘤位置预测模型的准确性受患者个体差异、测量精度、肿瘤-替代物相关性、训练数据选择和正弦二次模型对患者数据采集时可能无法屏住呼吸而产生数据采集不全的问题。

实现本发明目的的思路是：对于pls模型和正弦二次模型都存在的初始采集的数据就会受到患者差异影响的问题，本发明采用了主成分分析法对患者的数据进行降维处理，保留了不同患者之间的共性，极大的消除个体差异性；pls模型中标志物的移动会导致模型误差随时间越来越大。本发明采用了icp算法，通过对不同相位的胸腹表面图像的像素点进行配准计算出胸腹表面位移，而人体本身的胸腹表面的相对位置是相对固定的，使得模型不会随时间变差，也减少了患者植入标志物的痛苦。

本发明的系统包括图像分割模块，膈肌位移计算模块、胸腹表面位移模块、图像输出模块、pca模块、模型拟合模块；其中，

所述膈肌图像分割模块，用于选取所选患者的人体三维图像，将所选相位的人体三维图像以像素值差值法处理为仅有右肺下膈肌的三维图像；

所述胸腹表面图像分割模块，用于选取所选患者的所选相位的人体三维图像，将所选相位的人体三维图像以掩膜法处理为仅有胸腹表面的三维图像；

所述膈肌位移计算模块，将膈肌三维图像的全部膈肌像素点的空间三维坐标值取均值，将该均值作为右肺下膈肌的质心坐标值；将该坐标值与所选患者的右肺下膈肌的初始质心坐标值相减，得到所选相位的膈肌位移；

所述胸腹表面位移计算模块，用于使用icp算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标与呼吸周期的0％相位的胸腹表面图的空间三维坐标进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移；

所述pca模块，用于将每个患者全部相位的膈肌位移连接为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；用于将每个患者的归一化后的30维的向量，进行主成分分析法处理，降维并映射到d维子空间r^v中；用于将每个患者每个相位的胸腹表面位移归一化为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；用于将每个患者归一化后的30维向量，进行主成分分析法处理，降维并映射到d维子空间r^u中；用于将经过tssm模型计算后得到的膈肌位移按照pca特征子空间的系数和基，即计算得到归一化的膈肌位移数据按照膈肌位移的归一化公式，将归一化的膈肌位移数据逆归一化，得到最终的膈肌位移v^out；

所述模型拟合模块，用于利用线性岭回归优化法将映射到r^v中的膈肌位移数据与映射到r^u中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式β^opt，即为tssm模型。

所述图像输出模块在计算机上显示出横坐标为相位，纵坐标为膈肌位移的曲线图。

本发明方法的具体步骤包括如下：

(1)呼吸时相数据采集：

(1a)选取患有肺部癌症或胰腺癌的至少20个患者，利用4dct扫描仪获取每个患者连续的至少5个呼吸周期的时相数据，计算每个患者的平均呼吸周期，从5个呼吸周期中选取最接近平均呼吸周期的一个时相数据，作为该患者的呼吸周期；

(1b)用与4dct扫描仪连接的计算机打开每个患者的平均呼吸周期时相数据，将该患者呼吸周期从0％开始，每隔10％处选取一次时相数据重建成10个相位的人体三维图像；

(1c)从20个患者中随机选取16个患者，将所选16个患者的人体三维图像包含的全部时相数据集组成初始训练数据，剩余4个患者的人体三维图像作为测试数据；

(2)计算每个患者的4dct图像的膈肌位移和胸腹表面位移：

(2a)从16个患者中选取一个未选过的患者：

(2b)膈肌图像分割模块选取所选患者呼吸周期0％相位的患者体三维图像，利用像素值差值法，将所选相位的人体三维图像处理为仅有右肺下膈肌的三维图像；

(2c)胸腹表面图像分割模块选取所选患者的呼吸周期的0％相位的人体三维图像，将该图像以掩膜法处理为仅有胸腹表面的三维图像；

(2d)选取所选患者的一个未选过相位的人体三维图像，采用与步骤(2b)相同的处理方法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有右肺下膈肌的膈肌三维图像，膈肌位移计算模块将该膈肌三维图像的全部膈肌像素点的空间三维坐标取均值，将该均值作为右肺下膈肌的质心坐标；将质心坐标值与所选患者的右肺下膈肌的初始质心坐标相减，得到所选相位的膈肌位移；

(2e)采用与步骤(2c)中相同的处理方法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有的胸腹表面的胸腹部三维图像，胸腹表面位移计算模块用icp算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标值与呼吸周期的0％相位的胸腹表面图的空间三维坐标值进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移；

(2f)判断是否选完所选患者的所有相位的人体三维图像，若是，则执行步骤(2g)，否则，执行步骤(2d)；

(2g)判断是否选完16个患者，若是，则执行步骤(3)，否则，执行步骤(2a)；

(3)将膈肌位移数据和胸腹表面位移数据映射到各自的子空间中：

(3a)模块将每个患者全部相位的膈肌位移连接为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；

(3b)pca模块将每个患者的归一化后的30维的向量，进行主成分分析法处理，降维并映射到d维子空间r^v中；

(3c)pca模块将每个患者每个相位的胸腹表面位移归一化为一个10个x轴上的位移值+10个y轴上的位移值+10个z轴上的位移的30维的向量，并进行归一化；

(3d)pca模块将每个患者归一化后的30维向量，进行主成分分析法处理，降维并映射到d维子空间r^u中；

(4)模型拟合模块利用线性岭回归优化法将映射到r^v中的膈肌位移数据与映射到r^u中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式β^opt，即为tssm模型。

(5)实时输出测试患者的相位-膈肌位移的曲线图：

(5a)用4dct扫描仪获取测试患者的呼吸时相数据；

(5b)测试患者的呼吸时相数据依次经过：胸腹表面图像分割模块、胸腹表面位移计算模块、pca模块、tssm模型、pca模块的逆归一化处理、最后输出测试患者的膈肌位移；

(5c)图像输出模块在计算机上显示出横坐标为膈肌位移相位，纵坐标为膈肌位移的曲线图。

(6)进一步将tssm推广到基于核映射的非线性子空间：

将胸腹表面位移数据用经过核函数变换后的线性组合来表示，基于核函数的非线性岭回归优化变换ktssm模型。

本发明与现有技术相比具有以下优点：

第一，由于本发明系统中的胸腹表面图像分割模块用于从每个相位的人体三维中分割出胸腹表面图像，通过计算胸腹表面位移预测出膈肌位移，克服了现有技术在偏最小二乘回归模型中，通过计算外标志物位移计算出肿瘤位移的技术中需要在患者胸腹表面植入标志物从而对患者造成创伤的问题，使得本发明具有对患者的伤害更少的优点。

第二，由于本发明中采用了膈肌图像分割，并以所有膈肌像素点的平均值来表示膈肌质心的方法，克服了现有技术在正弦二次模型中，仅以膈肌上呼吸振幅最大的点表示膈肌质心，但膈肌的运动在左右方向上有一定的偏移，导致在振幅最高处和振幅最低处的点实际上并不是膈肌上同一个点的缺点，使得本发明具有选取的膈肌质心更加准确的优点。

第三，由于本发明方法中采用了主成分分析法，将患者的膈肌位移和胸腹表面位移降维并映射到仅包涵数据本质特征的子空间中，克服了现有技术中采集呼吸时相数据时就易受患者的个体差异、噪声的影响的问题，使得本发明具有更抗干扰的优点。

附图说明

图1是本发明系统的模块示意图；

图2是本发明方法的流程图；

图3是本发明仿真图；

图4是本发明的性能评价箱线图。

具体实施方式

下面结合附图对本发明做进一步的详细描述。

参照附图1，对本发明的系统做进一步的详细描述。

所述膈肌图像分割模块，用于选取所选患者的所选相位的患者体三维图像，将所选相位的人体三维图像以像素值差值法处理为仅有右肺下膈肌的三维图像。

所述胸腹表面图像分割模块，用于选取所选患者的所选相位的人体三维图像，将将所选相位的人体三维图像以掩膜法处理为仅有胸腹表面的三维图像。

所述膈肌位移计算模块，将膈肌三维图像的全部膈肌像素点的空间三维坐标取均值，将该均值作为右肺下膈肌的质心坐标；将该坐标与所选患者的右肺下膈肌的初始质心坐标相减，得到所选相位的膈肌位移。

所述胸腹表面位移计算模块，用于用icp算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标与呼吸周期的0％相位的胸腹表面图的空间三维坐标进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移。

所述pca模块，用于将每个患者的10个相位归一化后的向量，进行主成分分析法处理，得到降维并映射到d维子空间r^v中的；用于将每个患者的10个相位，归一化后的向量，进行主成分分析法处理，得到降维并映射到d维子空间r^v中的；用于将经过tssm模型计算得到的膈肌位移根据pca特征子空间的系数和基来计算膈肌位移的预测数据，即得到归一化的膈肌位移数据根据膈肌位移的归一化公式，将归一化的膈肌位移数据逆归一化，得到最终的膈肌位移v^out。

所述模型拟合模块，用于利用线性岭回归优化法将映射到r^v中的膈肌位移数据与映射到r^u中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式β^opt，即为tssm模型。

所述图像输出模块在计算机上显示出横坐标为相位，纵坐标为膈肌位移的曲线图。

参照附图2，对本发明的方法做进一步的详细描述。

步骤1，呼吸时相数据采集。

选取患有肺部癌或胰腺癌或肝癌的至少20个患者，利用4dct扫描仪获取每个患者连续的至少5个呼吸周期的时相数据，计算每个患者的平均呼吸周期，从5个呼吸周期中选取最接近平均呼吸周期的一个时相数据，作为该患者的呼吸周期。

用与4dct扫描仪连接的计算机打开每个患者的平均呼吸周期时相数据，将该患者呼吸周期从0％开始，每隔10％处选取一次时相数据重建成10个相位的人体三维图像。

从20个患者中随机选取16个患者，将所选16个患者的人体三维图像包含的全部时相数据集组成初始训练集；剩余4个患者的全部时相数据集组成测试集。

步骤2，计算出每个患者的4dct图像的膈肌位移和胸腹表面位移：

(2.1)从16个患者中选取一个未选过的患者：

(2.2)膈肌图像分割模块选取所选患者呼吸周期0％相位的患者体三维图像，将所选相位的人体三维图像以像素值差值法处理为仅有右肺下膈肌的三维图像，对右肺下膈肌的三维图像中全部膈肌像素点的空间三维坐标取均值，将该均值作为右肺下膈肌的初始质心坐标。

所述的像素值差值法的步骤如下：

(2.2.1)将所选相位的三维人体图像中的最大目标区域作为身体区域；通过形态学方法中的开运算，去除每个相位的三维人体图像中背景中不属于人体，但像素不为0的孤立区域；利用形态学方法中的闭运算填充身体区域中的孔，得到分割后的身体区域。

(2.2.2)从分割后的身体区域中，分割出肺部区域；对三维人体图像，将背景像素值均置为0，分割后的身体区域像素值均置为1，肺部区域像素值均置为3；得到三维人体分割图像。

(2.2.3)用三维roi掩膜对三维人体分割图像进行分割，得到从前向后看的肺部三维图像。

(2.2.4)对肺部三维图像，从左到右在矢状平面上计算肺部区域的各个切片的面积，生成横坐标为左右方向的三维roi掩膜切面，纵坐标为肺部面积的曲线图；选择曲线两个峰之间的最低点，作为中位点，定位出该点横坐标的坐标值。

(2.2.5)在三维人体分割图像中，将中位点以左的人体和肺部的像素点全部置为0；再从下到上找到像素值差值为2的分界面，该边界面上的全部像素组成右肺下表面三维图像。

(2.2.6)对该右肺下表面，从上方，先通过形态学运算中的腐蚀去除掉一些属于肺部，但不属于膈肌部分的像素点，再在经过腐蚀后的右肺下表面上填充被过度腐蚀掉的膈肌像素点。

(2.3)胸腹表面图像分割模块选取所选患者的呼吸周期的0％相位的人体三维图像，将该图像以掩膜法处理为仅有胸腹表面的三维图像。

所述的掩膜法的步骤如下：

(2.3.1)对人体三维图像中背景区域的所有像素值均设置为0、身体区域中的所有像素值均设置为1；沿x轴负方向找出像素值差值为1的边界面，将该边界面称为初始胸腹部表面。

(2.3.2)在初始胸腹部表面上沿y轴正方向找到最宽的位置，最宽位置对应人体的肩膀，去除肩膀之上的区域，进行腐蚀操作，得到了每个相位的胸腹表面图。

(2.3.3)将所选患者的所有相位的胸腹表面图取交集，得到胸腹部掩膜，将该掩膜与初始胸腹部表面取交集，得到最终胸腹部表面。

(2.4)选取所选患者的一个未选过相位的人体三维图像，采用像素值差值法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有右肺下膈肌的膈肌三维图像，膈肌位移计算模块将该膈肌三维图像的全部膈肌像素点的空间三维坐标取均值，将该均值作为右肺下膈肌的质心坐标。

(2.5)将所选患者所选相位的右肺下膈肌的质心坐标与所选患者的右肺下膈肌的初始质心坐标相减，得到所选相位的膈肌位移，其表达式如下：

其中，i＝1,2,……16表示训练集内的患者的序号，k＝1,2,……,10表示4dct图像的相位的序号，表示第i个患者的第k个相位的右肺下膈肌全部的像素点的质心，表示第i个患者呼吸周期的0％相位的右肺下膈肌全部的像素点的质心；表示第i个患者第k个相位右肺下膈肌全部的像素点在x轴上的坐标值的平均值，表示第i个患者呼吸周期的0％相位的右肺下膈肌全部的像素点在x轴上坐标值的平均值；表示第i个患者第k个相位右肺下膈肌全部的像素点在y轴上的坐标值的平均值，表示第i个患者的呼吸周期0％相位右肺下膈肌全部的像素点在y轴上的坐标值的平均值，表示第i个患者第k个相位右肺下膈肌全部的像素点在z轴上的坐标值的平均值，表示第i个患者的呼吸周期0％相位右肺下膈肌全部的像素点在z轴上的坐标值的平均值；表示第i个患者第k个相位右肺下膈肌质心的位移；表示第i个患者第k个相位右肺下膈肌的质心在x轴上的位移，表示第i个患者第k个相位右肺下膈肌的质心在y轴上的位移，表示第i个患者第k个相位右肺下膈肌的质心在z轴上的位移。

为了便于表达，本发明后续将忽略平均符号和差值符号，使用vi(k)＝{oi(k),pi(k),qi(k)},k＝1,2,……,10表示第i个病人第k个相位的膈肌位移。

(2.6)采用掩膜法，对所选相位的人体三维图像进行处理，得到所选相位的人体三维图像中仅有的胸腹表面的胸腹部三维图像，胸腹表面位移计算模块用icp算法将所选患者的所选相位的胸腹部三维图像的空间三维坐标与呼吸周期的0％相位的胸腹表面图的空间三维坐标进行配准，将配准后得到的平移向量作为所选相位的胸腹表面位移，其表达式如下：

其中，表示第i个患者第k个相位的胸腹表面的位移，表示第i个患者第k个相位的沿x轴方向上的位移，表示第i个患者第k个相位的沿y轴方向上的位移，表示第i个患者第k个相位的沿z轴方向上的位移。

为了便于表达，本发明后续将忽略平均符号和插值符号，使用ui(k)＝{ri(k),si(k),ti(k)},k＝1,2,……,10表示第i个病人第k个相位的胸腹部表面位移。

所述icp算法是指：

(2.6.1)将呼吸周期0％相位的胸腹表面图的全部胸腹表面像素点的三维坐标作为源点云m，取点集mj∈m。

(2.6.2)将所选相位的胸腹表面图的全部胸腹表面像素点的三维坐标作为目标点云n，取点集nj∈n。

(2.6.3)计算映射矩阵h，使误差函数e(g,l)最小，映射矩阵h和误差函数e(g,l)如下：

其中，j表示点集mj和nj中每个点的序列号，j表示j的最大值，gj×j表示j阶旋转矩阵，lj×1表示j维平移向量，o1×j表示j维零向量，表示使{·}取得最小值时的g和l的取值。

(2.6.4)对点集nj进行旋转和平移，计算点集nj'＝{nj'＝gnj+l,nj∈n}。

(2.6.5)计算映射矩阵h'，使误差函数e(g',l')最小，映射矩阵h'和e(g',l')如下：

第6步，判断误差函数e(g',l')是否小于阈值或者大于预设的最大迭代次数，若是，则停止迭代更新，对所有平移向量求和，即为将所选相位的胸腹表面图的全部胸腹表面像素点与呼吸周期0％相位的胸腹表面图的全部胸腹表面像素点，进行配准后的平移向量，并执行(2.7)；否则，执行本步骤的(2.6.4)。

(2.7)判断是否选完所选患者的所有相位的人体三维图像，若是，执行本步骤的(2.8)，否则，执行本步骤的(2.4)。

(2.8)判断是否选完16个患者，若是，执行步骤3，否则，执行本步骤的(2.1)。

步骤3，将膈肌位移数据和胸腹表面位移数据映射到各自的子空间中。

pca模块利用主成分分析法，将膈肌位移vi和胸腹表面位移ui，进行归一化，降维并映射到各自的d维子空间r^v，r^u中。

所述的主成分分析法是指：

(3.1)分别将每个患者的10个相位的膈肌位移值中每个相位的x轴、y轴、z轴的位移值，分别组成一个10维向量，再将这3个向量，按照x轴向量+y轴向量+z轴向量，组成一个30维的向量vi；

vi＝{oi(1),oi(2),...,oi(10),pi(1),pi(2),...,pi(10),qi(1),qi(2),...qi(10)}

采用与本步骤的(3.1)相同的处理方法，将每个患者的10个相位的胸腹表面位移值组成一个30维的向量ui；

ui＝{ri(1),ri(2),...,ri(10),si(1),si(2),...,si(10),ti(1),ti(2),...ti(10)}

(3.2)对向量vi，ui分别进行归一化处理，得到和归一化公式如下：

其中，表示第i个患者膈肌位移归一化后的30维向量，i表示所选患者的序号，vi(n)表示向量vi内每一个元素，表示第i个患者膈肌位移向量vi(n)中30个元素中最小的元素，表示第i个患者膈肌位移向量vi(n)中30个元素中最大的元素；表示第i个患者胸腹表面位移归一化后的30维向量，ui(n)表示向量ui内每一个元素，第i个患者胸腹表面位移3向量ui(n)30个元素中最小的元素，表示第i个患者胸腹表面位移向量ui(n)30个元素中最大的元素。

(3.3)按照下式，计算协方差矩阵c和b：

其中，表示第i个患者膈肌位移归一化后的30维向量，表示的转置，表示第i个患者胸腹表面位移归一化后的30维向量，表示的转置。

(3.4)按照下式，计算协方差矩阵b的特征值及对应的特征向量：

det(λi-b)＝0

(λi-b)e＝0

其中，det(·)表示行列式符号，λ表示协方差矩阵b的特征值，i表示和协方差矩阵b阶数相同的单位矩阵，e表示特征值λ对应的特征向量。

按照与本步骤的(3.4)相同的处理方式，计算协方差矩阵c的特征值及对应的特征向量。

(3.5)对协方差矩阵b的全部特征值，从大到小选择大于阈值ε的d个特征值，阈值ε在构造域的过程中，需要在数据完整性和数据冗余之间进行平衡。

按照与本步骤的(3.5)相同的处理方式，选择协方差矩阵c的d个特征值。

(3.6)根据d个特征值，得到d个对应的特征向量，该d个特征向量分别作为膈肌位移数据和胸腹部表面位移数据的基，表示为ωv和ωu，ωu,ωv∈r^30×d，分别构建出d维子空间r^v和r^u。

(3.7)和是正交化的(且其中id是d的单位矩阵)；分别通过运算和将膈肌位移数据和胸腹表面位移数据映射到各自的子空间中，其中，是在子空间r^u中的对应点，是在子空间r^v中的对应点。

步骤4，模型拟合模块利用线性岭回归优化法将映射到r^v中的膈肌位移数据与映射到r^u中的胸腹部表面位移数据连接起来，得到变换矩阵β，求出β的最优化表达式β^opt，即为tssm模型。

所述的变换矩阵β的公式如下：

其中，表示当{·}取最小值时的β的值，||·||2表示2范数，λ≥0是控制收缩量的规则化参数。

所述的最优表达式β^opt操作是指：对β求导，并使等式等于0，则可以获得最优表达式，如下所示：

其中，opt表示最优，i表示单位矩阵。

tssm模型表达式如下：

其中，表示测试数据经过tssm模型的输出结果，表示测试患者的胸腹表面位移在d维子空间r^u的表示。

步骤5，实时输出测试患者的相位-膈肌位移的坐标图。

用4dct扫描仪获取测试患者的呼吸时相数据。

将测试患者的呼吸时相数据依次经过：胸腹表面图像分割模块、胸腹表面位移计算模块、pca模块、tssm模型、pca模块的逆归一化处理，最后输出测试患者的膈肌位移。

在计算机上显示出横坐标为相位，纵坐标为膈肌位移的曲线图。

步骤6，推导出基于核函数的非线性子空间的ktssm模型。

用替换掉tssm中的用k(·,·)替换掉两个高维子空间w^u中的两个向量的内积，得到非线性岭回归优化的变换矩阵β^'opt和ktssm模型。

其中，k(·,·)表示映射的核函数，表示测试数据中d维子空间σ中的胸腹表面位移，表示训练数据中d维子空间r^u中的胸腹表面位移，φ表示非线性子空间r^u到高维子空间w^u的映射。

利用替换掉非线性子空间r^u的胸腹表面位移，得到基于核函数的非线性岭回归优化变换表达式如下：

其中，β'表示非线性岭回归优化变换矩阵。

岭回归优化方程的最优解表示为：

ktssm模型表达式如下：

其中，表示使用ktssm模型计算后得到的d维子空间r^v中的膈肌位移，

表示训练数据中的d维子空间r^v中的膈肌位移。

本发明的效果可通过以下仿真进一步说明。

1.仿真实验条件：

本发明的仿真实验的硬件平台为：处理器为interxeon3.6ghz，32g内存。

本发明的仿真实验的软件配置为windows10操作系统和matlabr2019b。

2.仿真内容：

本发明的仿真是用本发明的系统和本发明的方法，采用了来自20位实验对象的4dct数据集，每个实验对象包括10个相位(具体为平均呼吸周期的0％、10％、…、90％)。随机选出16个患者的呼吸时相数据作为训练数据，剩余4名患者的呼吸时相数据作为测试数据。

本发明仿真实验是采用本发明的tssm模型和2个ktssm模型和一个现有技术(偏最小二乘回归模型)分别对输入的训练患者的呼吸时相数据和测试患者的呼吸时相数据进行计算，获得相位-膈肌位移的坐标图。

在仿真实验中采用的两个ktssm模型是指：

两个ktssm模型是指本发明步骤6中推导出基于核函数的非线性子空间的ktssm模型：多项式核的ktssm模型和高斯核的ktssm模型。

在仿真实验中采用的一个现有技术偏最小二乘法模型是指：

malinowski等人在“mitigatingerrorsinexternalrespiratorysurrogate-basedmodelsoftumorposition，internationaljournalofradiationoncology,biology,physics,vol.82,no.5,2012”中提到的偏最小二乘回归模型，简称pls模型。

下面结合图3仿真图对本发明的效果做进一步的描述。

图3对应四种算法的膈肌位移的预测结果和实际结果的曲线图，图3中有4个子图，分别对应4种模型的预测结果曲线和实际结果的曲线图，在每个子图中，又有4张图，每张图对应一名测试患者的预测结果和实际结果的曲线图。图3中的x轴表示十个相位，y轴表示膈肌位移值，单位是毫米(mm)。图3中以标示的曲线表示实际结果的曲线，即通过膈肌图像分割，膈肌位移计算得到的值，图3中以“米”字标示的曲线表示预测结果的曲线，即将胸腹表面位移输入预测模型后得到的值。图3(a)为对4名患者采用pls模型计算出的预测结果和实际结果的曲线图，图3(b)为对4名患者采用线性tssm模型计算出的预测结果和实际结果的曲线图，图3(c)为4名患者采用多项式核的ktssm模型计算出的预测结果和实际结果的曲线图，图3(d)为4名患者采用高斯核的ktssm模型计算出的预测结果和实际结果的曲线图。两条曲线重合的越好，表明该模型的预测效果越好。图3(a)中的4名患者的预测值与实际值之间的误差都明显大于图3(b)，说明tssm模型相比pls模型具有更好的预测效果。

3.仿真结果分析。

图4表示在膈肌位移在z轴方向上每个相位的预测性能评价的箱线图，每个相位都包含100次独立运行的结果。图4中，在每个方框上，中心的标记表示中值，而方框的底部和顶部边缘分别表示100次绝对误差结果中的从小到大排列的第25％和第75％的绝对误差结果，“+”符号表示异常值。图4中,x轴表示十个相位，y轴表示预测值和实际值之间的估计误差绝对值，单位是毫米(mm)。图4(a)为对pls模型进行预测性能评价的箱线图，图4(b)为对线性tssm模型进行预测性能评价的箱线图，图4(c)为对多项式核的tssm模型进行预测性能评价的箱线图，图4(d)为对高斯核的tssm模型进行预测性能评价的箱线图。在图4(a)中的预测误差值比图4(b)、(c)、(d)中的误差值大得多，说明pls模型相比tssm模型的预测误差更大。

利用三个指标进行定量评估四种模型的性能：均方误差mse、r²误差和平均绝对百分比误差mape，利用下面公式计算mse、r²、mape，将所有计算结果绘制成表1：

其中，f＝1,2,3,4表示测试集内患者序号，是膈肌实际的位移值，即将膈肌图像从人体三维图像中分割出来计算得到的膈肌位移值，是模型的预测值，即将胸腹表面位移值输入进模型后经过计算得到的膈肌位移值，σ表示求和。mse衡量模型的预测均方误差，mse的值越低，预测结果越好。r²比较模型的优劣，消除数据值域对实验结果的影响。r²的范围为[-∞,1]，值越大预测效果越好；mape衡量模型的预测值与实际值的平均偏差，平均偏差值越低预测结果越好。

表1主方向预测结果的统计结果一览表

如表1所示，每个表格的前一个值代表100次独立运行的平均值，(·)中的值代表100次独立运行的标准偏差，本发明仿真实验记录mse和mape的最小值，记录r²的最大值。从表1中可见由顶部到底部，噪声系数依次增加，tssm模型相比pls模型，mse值更小，r²值更大，mape值更小，显示了tssm算法优于pls模型；通过tssm算法获得的mse值和r²值非常令人满意；由于向量值的不平衡性，mape值大于100％。