面向放疗的人体胸腹表面运动预测方法与流程

本发明面向放疗的人体胸腹表面运动预测方法涉及精密测量、肿瘤医学、统计数学等技术领域。

背景技术：

放射治疗(后文简称为放疗)是治疗癌症最主要的手段之一。大约70％的癌症患者在治疗癌症的过程中会用到放疗，约有40％的癌症可以通过放疗治愈。放射治疗效果与放射线照射到肿瘤区域的精度和剂量有关。消化系统恶性肿瘤位于胸腔腹腔内，受呼吸运动影响，肿瘤的位置和体积都会随时间而改变，造成放射线照射精度降低，从而影响放疗效果。

为了解决呼吸运动降低放疗效果的问题，很多学者开展了靶区扩边、屏气、呼吸门控、被动加压、四维放疗和实时跟踪等技术或方法的研究工作。这些工作有效降低了呼吸运动对放疗效果的影响。然而，由于这些工作都是在呼吸运动发生后开展的滞后补偿工作，因此无法从根本上克服呼吸运动降低放疗效果的问题。

如果更够开展超前补偿方法，就有希望能够从根本上克服呼吸运动降低放疗效果的问题。要想实现对呼吸运动超前补偿，首先需要对呼吸运动进行预测。目前，呼吸运动方法可以分为以下两大类：第一、模型预测方法，这种方法通过监测前期呼吸运动数据，建立呼吸运动模型，再基于模型来推测未来的呼吸运动；第二、无模型预测方法，这种方法通过对呼吸运动数据进行训练，形成启发式的学习算法，通过学习来预测未来的呼吸运动。

虽然呼吸运动不断重复吸气和呼气动作，但是呼吸运动却并不是以吸气和呼气为一个周期进行简单重复，在每次呼吸过程中，受自身以及外界的影响，运动规律都会发生变化。模型预测方法无法照顾到这种变化，因此预测精度会受到制约，所以，更多学者开始在无模型预测方法中进行尝试。

技术实现要素：

为了实现对呼吸运动进行准确预测，同时为了照顾到呼吸运动规律的变化，本发明提出了一种面向放疗的人体胸腹表面运动预测方法，该方法的最大特点是并不是利用先验知识对未知量进行直接预测；而是采用对特征量和发生时间分别进行预测的方式，并结合特征量和发生时间与未知量发生时间的关系，来实现对未知量的预测，该方法能够充分利用特征量的有序性和惯性，实现对未知量的更准确预测。

本发明的目的是这样实现的：

面向放疗的人体胸腹表面运动预测方法，包括以下步骤：

步骤a、利用fasttrack，采集胸腹表面呼吸运动数据flow-high(x)；

步骤b、利用高通滤波器，滤除步骤a得到的胸腹表面呼吸运动数据flow-high(x)中的低频成分，去除放疗过程中人体自身除呼吸运动外的其他运动，得到胸腹表面呼吸运动数据fhigh(x)；

步骤c、利用低通滤波器，滤除步骤b得到的胸腹表面呼吸运动数据fhigh(x)中的高频成分，去除胸腹表面呼吸运动数据采集过程中的噪声，得到胸腹表面呼吸运动数据f(x)；

步骤d、从步骤c得到的胸腹表面呼吸运动数据f(x)中获取特征数据，作为训练集，包括以下步骤：

步骤d1、采集胸腹表面呼吸运动数据的极大值，这些极大值为：fmax(i),i＝1,2,3,…,n；

步骤d2、确定步骤d1所述极大值发生的时间，这些极大值发生时间为：tmax(i),i＝1,2,3,…,n；

步骤d3、采集胸腹表面呼吸运动数据的极小值，这些极小值为：fmin(i),i＝1,2,3,…,n；

步骤d4、确定步骤d3所述极小值发生的时间，这些极小值发生时间为：tmin(i),i＝1,2,3,…,n；

步骤e、根据步骤d得到的特征数据，预测下一组特征数据，包括以下步骤：

步骤e1、预测下一个极大值fmax(n+1)，包括以下步骤：

步骤e11、构造分辨率为n×n的矩阵k，在矩阵k中，第i行第j列元素k(i,j)为：

步骤e12、按照如下公式，计算fmax(n+1)

fmax(n+1)＝k(n,:)k^-1(fmax')

其中，k(n,:)表示矩阵k的第n行，k^-1表示矩阵k的逆矩阵，fmax'表示胸腹表面呼吸运动数据的极大值fmax的转置；

步骤e2、预测下一个极大值发生时间，包括以下步骤：

步骤e21、按照如下公式，计算in1

步骤e22、按照如下公式，计算in2

步骤e23、按照如下公式，计算tmax

步骤e24、按照如下公式，计算intmax

步骤e25、按照如下公式，计算下一个极大值发生时间tmax(n+1)

步骤e3、预测下一个极小值fmin(n+1)，包括以下步骤：

步骤e31、构造分辨率为n×n的矩阵k，在矩阵k中，第i行第j列元素k(i,j)为：

步骤e32、按照如下公式，计算fmin(n+1)

fmin(n+1)＝k(n,:)k^-1(fmin')

其中，k(n,:)表示矩阵k的第n行，k^-1表示矩阵k的逆矩阵，fmin'表示胸腹表面呼吸运动数据的极小值fmin的转置；

步骤e4、预测下一个极大值发生时间，包括以下步骤：

步骤e41、按照如下公式，计算in1

步骤e42、按照如下公式，计算in2

步骤e24、按照如下公式，计算t

步骤e24、按照如下公式，计算intmin

步骤e25、按照如下公式，计算下一个极大值发生时间tmin(n+1)

步骤f、根据给定的时间t，预测胸腹表面呼吸运动数据f(t)，包括以下步骤：

步骤f1、判断时间t位于最后一个从极小值向极大值变化的区间，还是位于最后一个从极大值向极小值变化的区间；

如果，时间t位于最后一个从极小值向极大值变化的区间，计算每一个从极小值向极大值变化的区间内，对应时间下的对应幅值比例；

如果，时间t位于最后一个从极大值向极小值变化的区间，计算每一个从极大值向极小值变化的区间内，对应时间下的对应幅值比例；

步骤f2、将计算得到的对应幅值比例做平均；

步骤f3、根据步骤f2得到的比例平均值，计算得到胸腹表面呼吸运动数据f(t)。

有益效果：

第一、本发明面向放疗的人体胸腹表面运动预测方法的最大特点就是并不是利用先验知识对未知量进行直接预测，而是采用对特征量和发生时间分别进行预测的方式，并结合特征量和发生时间与未知量发生时间的关系，来实现对未知量的预测，这种通过预测特征量的间接预测方法能够充分利用特征量的有序性和惯性，实现对未知量的更准确预测。

第二、在本发明中，极值的预测和发生时间的预测采用了两种不同的方法，这是分别针对极值的发生规律与发生时间的规律不同而设计，这样的设计，可以保证提高任何一个参数的预测精度都不会以牺牲另一个预测精度为代价，有利于提高最终未知量的预测精度。

第三、在对极值进行预测的时候，相关性矩阵的构造是影响预测精度的关键，在本发明中，矩阵k中的元素计算公式的构造也是本发明的创新点，该计算公式直接决定了本发明具有较高的预测精度。

除此之外，还有以下两点需要说明：

第一、本发明虽然依附肿瘤医学的背景，但是本发明是从技术角度出发，是对胸腹表面呼吸运动数据进行提取和预测，并且从本发明的技术目的和效果来看，无法实现对疾病进行诊断和治疗，因此不属于专利法第二十五条所述的疾病的诊断和治疗方法，即本发明不存在客体问题。

第二、本发明公开足够充分，只要将说明书具体实施例二所记载的源代码在matlab软件上运行，即可得到运行结果，无论是本领域技术人员还是非本领域非技术人员，都能够实现本申请。

附图说明

图1是本发明方法仿真结果。

图2是执行本发明方法的matlab软件运行界面截图。

具体实施例

下面结合附图对本发明具体实施例作进一步详细描述。

实施例一

本实施例是面向放疗的人体胸腹表面运动预测方法的实施例。

本实施例的面向放疗的人体胸腹表面运动预测方法，包括以下步骤：

步骤a、利用fasttrack，采集胸腹表面呼吸运动数据flow-high(x)；

步骤b、利用高通滤波器，滤除步骤a得到的胸腹表面呼吸运动数据flow-high(x)中的低频成分，去除放疗过程中人体自身除呼吸运动外的其他运动，得到胸腹表面呼吸运动数据fhigh(x)；

步骤c、利用低通滤波器，滤除步骤b得到的胸腹表面呼吸运动数据fhigh(x)中的高频成分，去除胸腹表面呼吸运动数据采集过程中的噪声，得到胸腹表面呼吸运动数据f(x)；

步骤d、从步骤c得到的胸腹表面呼吸运动数据f(x)中获取特征数据，作为训练集，包括以下步骤：

步骤d1、采集胸腹表面呼吸运动数据的极大值，这些极大值为：fmax(i),i＝1,2,3,…,n；

步骤d2、确定步骤d1所述极大值发生的时间，这些极大值发生时间为：tmax(i),i＝1,2,3,…,n；

步骤d3、采集胸腹表面呼吸运动数据的极小值，这些极小值为：fmin(i),i＝1,2,3,…,n；

步骤d4、确定步骤d3所述极小值发生的时间，这些极小值发生时间为：tmin(i),i＝1,2,3,…,n；

步骤e、根据步骤d得到的特征数据，预测下一组特征数据，包括以下步骤：

步骤e1、预测下一个极大值fmax(n+1)，包括以下步骤：

步骤e11、构造分辨率为n×n的矩阵k，在矩阵k中，第i行第j列元素k(i,j)为：

步骤e12、按照如下公式，计算fmax(n+1)

fmax(n+1)＝k(n,:)k^-1(fmax')

其中，k(n,:)表示矩阵k的第n行，k^-1表示矩阵k的逆矩阵，fmax'表示胸腹表面呼吸运动数据的极大值fmax的转置；

步骤e2、预测下一个极大值发生时间，包括以下步骤：

步骤e21、按照如下公式，计算in1

步骤e22、按照如下公式，计算in2

步骤e23、按照如下公式，计算tmax

步骤e24、按照如下公式，计算intmax

步骤e25、按照如下公式，计算下一个极大值发生时间tmax(n+1)

步骤e3、预测下一个极小值fmin(n+1)，包括以下步骤：

步骤e31、构造分辨率为n×n的矩阵k，在矩阵k中，第i行第j列元素k(i,j)为：

步骤e32、按照如下公式，计算fmin(n+1)

fmin(n+1)＝k(n,:)k^-1(fmin')

其中，k(n,:)表示矩阵k的第n行，k^-1表示矩阵k的逆矩阵，fmin'表示胸腹表面呼吸运动数据的极小值fmin的转置；

步骤e4、预测下一个极大值发生时间，包括以下步骤：

步骤e41、按照如下公式，计算in1

步骤e42、按照如下公式，计算in2

步骤e24、按照如下公式，计算t

步骤e24、按照如下公式，计算intmin

步骤e25、按照如下公式，计算下一个极大值发生时间tmin(n+1)

步骤f、根据给定的时间t，预测胸腹表面呼吸运动数据f(t)，包括以下步骤：

步骤f1、判断时间t位于最后一个从极小值向极大值变化的区间，还是位于最后一个从极大值向极小值变化的区间；

如果，时间t位于最后一个从极小值向极大值变化的区间，计算每一个从极小值向极大值变化的区间内，对应时间下的对应幅值比例；

如果，时间t位于最后一个从极大值向极小值变化的区间，计算每一个从极大值向极小值变化的区间内，对应时间下的对应幅值比例；

步骤f2、将计算得到的对应幅值比例做平均；

步骤f3、根据步骤f2得到的比例平均值，计算得到胸腹表面呼吸运动数据f(t)。

实施例二

本实施例是以具体仿真程序来实现本发明目的的仿真实施例。

本实施例所记载的程序在matlab软件上可以直接运行。在本实施例中，呼吸运动被简化成了三角波，特征点为三角波的极大值、极小值以及这两个极值发生的时间。相邻两个极值之间直线连接。共有十组极大值和极小值数据，其中，前九组作为训练集，用于预测第九组和第十组之间某一时刻的未知量，而第十组数据作为标准，对未知量的预测进行对比。

程序如下：

程序运行结果如图1所示，程序运行界面如图2所示。从程序运行结果可以看出，预测数据与真实数据之间的绝对误差仅为0.0172，相对误差仅为2.3％，说明本发明具有较高的预测精度。