CN101794357A - 呼吸同步跟踪系统的位置信号的模糊预测方法及预测器 - Google Patents

Abstract

呼吸同步跟踪系统在治疗射线投放时，可以监测病人的呼吸，并对肿瘤目标运动作出补偿。在治疗射线投放系统的指令和响应之间存在时间延迟，这种延迟可能会导致意外的、甚至危险的振动。为了提高同步跟踪性能，可以用信号预测技术来弥补时间延迟。本发明提供了一种新的同步跟踪系统的位置信号的模糊预测方法以及预测器，该预测器能够处理随时间变化的幅度和频率问题，比较好的解决了放射治疗中同步跟踪的问题。

Description

呼吸同步跟踪系统的位置信号的模糊预测方法及预测器

技术领域

本发明涉及一种位置信号的预测方法及预测器，特别是呼吸同步跟踪系统的位置信号的模糊预测方法以及按照这种方法设计的预测器。

背景技术

放射治疗(简称放疗)已经在医学临床得到广泛应用，成为治疗恶性肿瘤的主要手段之一。精确放射治疗的目标是最大限度地杀灭肿瘤细胞并有效保护周围正常组织和重要器官。研究表明，剂量精度每提高1％，治愈率便可提高2％，因此，提高放射治疗的精确度是非常重要的。

呼吸在引起体外(胸部/腹部)表面的运动的同时，还会导致体内肿瘤(肺部，肝部和胰腺部位的肿瘤)的运动。体内肿瘤由于受呼吸和器官本身的运动影响，其位置有很大的不确定性。通过影像检查显示胸部和腹部肿瘤和器官的移动度在10～30mm，有时甚至超过30mm，对由于呼吸引起的运动肿瘤，比如肺部和肝脏肿瘤，图像引导系统需要对肿瘤运动实时追踪，从而保证在整个治疗射线投放过程中，能够动态地移动或控制治疗射线束，使其准确地投放到病灶区，这需要通过呼吸运动跟踪系统来实现。

图像引导系统拍摄患者体部高清晰的射线影像，通过事先植入肿瘤内部或附近的金属标记物，构造空间立体结构，通过迭代算法精确算出肿瘤的空间位置；治疗中，通过实时的金标影像与定位CT数字生成的DRR中的金标影像位置的对比，可自动进行校正，这称为金标追踪。对于随呼吸运动的肿瘤，记录患者皮肤运动，形成相应的呼吸曲线。在这个动态连续的呼吸曲线过程中，通过摄影系统在不同的呼吸时相拍摄人体影像，并通过其中的金标确立一个呼吸模型关系。静态金标影像与动态的外在呼吸曲线相结合确立的关系能准确的反映体内随呼吸运动的靶区的位移。这样就可以通过病人体表的位置计算出体内肿瘤的位置。

带有图像引导的放疗系统的治疗射线投放装置的控制指令和实际的响应之间存在时间延迟。处理不好，这种时间延迟将可能在治疗期间导致意想不到的、甚至危险的振动。为了获得更好的跟踪性能，位置信号的预测可以用来补偿时间延迟。在当前的国际和国内流行的商用放疗系统里，呼吸运动跟踪系统中的同步跟踪子系统采用的预测方法通常是模式匹配法(Pattern-Matching)。和最小均方(Least-Mean-Square或LMS)。模式匹配算法的总体思路是：根据到当前位置为止的一段历史位置记录，定义一个模式(即以前的一些位置点)，在该段历史位置记录历史上搜索与预先定义的模式最佳匹配的模式；再从该匹配模式获取最好的“预测”。也就是说，Y’_p＝Y_0+p，其中Y′p是我们想要的预测值，Y₀是最好的匹配模式，而p是预测步骤。LMS算法是一种最速下降算法地近似算法，它使用了一个瞬时估计来近似代价函数的梯度矢量。梯度的估计数是根据对输入向量和误差信号的采样值。该算法遍历每个滤波器的系数，并沿近似梯度方向移动。

在实际的呼吸同步跟踪系统中，由于人体呼吸作用的影响，胸/腹部位置以及体内肿瘤的位置信号是随时间周期性变化的；更关键的是，其振幅和周期(或频率)也随时间变化。如果目前的呼吸模式在以前的时间框架里发生过并在历史记录里保存下来，那么模式匹配能够提供完善的多步以后的位置信号预测。一旦上述条件不成立，它的预测精度将恶化。在某些情况下可能还会导致预测间断问题。

因此，我们需要一种更好的预测方法，来解决以上存在的问题。

发明内容

本发明针对现有技术的不足，提供了一种新的基于模糊算法的同步跟踪系统的位置信号的预测方法，以解决现有技术中存在的一些困难问题，为放射治疗的同步跟踪提供了一种新的途径。

考虑到如图2所示的肿瘤位置随呼吸而变化，而且其随时间曲线的幅度和周期也是变化的，本发明的具体步骤如下：采集病人体表某点的位置Po并计算与其对应的肿瘤位置P(x，y，z)在一段时间内的位置参数(x，y，z)并且存入计算机，根据位置参数计算出相应的速度参数(dx，dy，dz)并存入计算机：

${\mathrm{dx}}_{i+1}=\frac{x_{i+1}-x_i}{t_{i+1}-t_i}$

${\mathrm{dy}}_{i+1}=\frac{y_{i+1}-y_i}{t_{i+1}-t_i}$

${\mathrm{dz}}_{i+1}=\frac{z_{i+1}-z_i}{t_{i+1}-t_i};$

根据计算机内存储的位置参数(x，y，z)，找出最大和最小的位置(x_max，y_max，z_max)和(x_min，y_min，z_min)，根据计算机内存储的速度参数(dx，dy，dz)，找出最大和最小的速度(dx_max，dy_max，dz_max)和(dx_min，dy_min，dz_min)；

制定模糊规则：根据最大和最小的位置参数，把x模糊划分出n₁个子集，这n₁个子集的每一个子集再被根据最大和最小的速度参数模糊划分为n₂个子集，得到M＝n₁*n₂个子模型，每个子模型都对应一个规则，在这里，我们以n₁取3为例加以说明，n₁＝3，即得到大(Big)，中(Medium)，小(Small)三个不同的子集，如图3所示。接下来，每个子集可以根据速度再分为n₂个子集，这里我们n₂也取3，即负(Negative)，零(Zero)和正(Positive)。这样我们一共就有M＝n₁*n₂个子模型。每个子模型都给出了一个线性方程组，其广义方程是这样的：

${\hat{y}}_m(k+p)=a_{m,0}y\left(k\right)+...+a_{m,r}y(k-r)+c_m$

是第m规则的模式输出，m的范围是在1~M，而M是方程的总数，p是预测步骤，r是预测阶数，可以选择2-15中任意一个数字，a和c是参数。在n₁取3，n₂取3的情况下，这些规则是：

R1：如果y(k)是Big而且dy(k)是Positive那么

${\hat{y}}_1(k+p)=a_{\mathrm{1,0}}y\left(k\right)+...+a_{1,r}y(k-r)+c_1$

R2：如果y(k)是Big而且dy(k)是Zero那么

${\hat{y}}_2(k+p)=a_{\mathrm{2,0}}y\left(k\right)+...+a_{2,r}y(k-r)+c_2$

……

Rm：如果y(k)是Medium而且dy(k)是Positive那么

${\hat{y}}_m(k+p)=a_{m,0}y\left(k\right)+...+a_{m,r}y(k-r)+c_m$

……

RM：如果y(k)是Small而且dy(k)是Negative那么

${\hat{y}}_M(k+p)=a_{M,0}y\left(k\right)+...+a_{M,r}y(k-r)+c_M$

那么，最后的预测结果即是以上子模型输出的综合结果的加权平均值

$\hat{y}(k+p)=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_i\·{\hat{y}}_i(k+p)}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_i}$

这里μ_i是模糊模型的隶属度，M＝9是模糊规则的数目。

采用同样的步骤，我们可以得到我们所想要的x和z的预测。

这样，我们就得到了位置P(x，y，z)的预测值。

参数的训练

所有的参数可以用梯度下降算法来训练，根据各个子模型所对应的计算机内已经存储的已知的位置参数(x，y，z)、速度参数(dx，dy，dz)以及该体表位置在相应时间点上的所对应的体内病灶的位置参数采用梯度下降算法，训练出对应子模型的参数a_m和c_m；

根据存储的点的历史记录，对参数a_m和c_m进行训练。

先定义成本函数J如下：

$J=\frac{1}{2}{e}^2=\frac{1}{2}{(y(k+p)-y_m(k+p))}^2$

每个参数的更新如下：

a_i，new＝a_i，old+Δa_i

c_new＝c_old+Δc

a和c的初始值可以在-0.001到0.001之间随机选择，但不为0这里，

${\mathrm{\Δa}}_i={\mathrm{\η}}_a\frac{\∂J}{\∂a_i}=-{\mathrm{\η}}_{a}e\frac{{\∂y}_m(k+p)}{{\∂a}_i}=-{\mathrm{\η}}_a\mathrm{ey}(k-i)$

$\mathrm{\Δc}={\mathrm{\η}}_c\frac{\∂J}{\∂c}=-{\mathrm{\η}}_{c}e\frac{{\∂y}_m(k+p)}{\∂c}=-{\mathrm{\η}}_{c}e$

η_a和η_c是学习率.惯量β可以被添加到Δa_i和Δc来加速训练的收敛。

我们η_a和η_c取0.3，β取0.01

Δa_i(k)＝-η_aey(k-i)+βΔa_i(k-1)

Δc(k)＝-η_ce+βΔc(k-1)

根据之前存储的数据，对a和c进行训练，得到a和c的值。

对随机选来的10组病人数据进行验证后，我们得到了类似的结果。表1显示了不同预测器对其中一个病人的比较结果。这里p＝3，η_a和η_c为0.3，β为0.01。

表1不同预测器对同一病人的比较结果

(单位为mm。ex，ey和ez分别是x，Y和Z方向的预测误差，|ex|表示ex的绝对值，Max是最大值，Mean是平均值，#|ex|＞1则表示|ex|中大于1mm的样本点数目，样本总数是12.5万)

	模式匹配	LMS	自适应模糊预测
				Max(|ex|)	3.42	3.63	2.59
Max(|ey|)	4.96	4.80	3.15
				Max(|ez|)	5.78	4.03	3.44
Mean(ex)	-0.0039	0.0006	0.0002
				Mean(ey)	0.0076	-0.0014	0.0019
Mean(ez)	0.0080	-0.0015	-0.0017
				#|ex|＞4	0	0	0
#|ey|＞4	5	1	0
				#|ez|＞4	22	0	0
#|ex|＞2	69	19	41
				#|ey|＞2	192	22	38
#|ez|＞2	336	21	105
				#|ex|＞1	461	38	123
#|ey|＞1	1321	61	301
				#|ez|＞1	1869	113	591

预测的位置将被用来指挥治疗射线投放装置作相应的移动。如果预测误差较大，人体器官中好的一部分就有可能将受射线的影响。如果预测误差有太多的高频成分，可能会导致治疗射线投放装置振动。此外，一般期望射线能够均匀地治疗肿瘤。因此，最好是使得平均预测误差接近于零。

从表1所列的测试结果和详细的数据观测，我们发现该自适应模糊预测器能够处理随时间变化的幅度和频率问题。具体而言：

Ο与模式匹配相比，它的预测中的大值误差的数目更小。

Ο与模式匹配和LMS相比，它的预测的最大误差最小。

Ο与模式匹配和LMS相比，它的预测的平均误差最小。

Ο不象LMS，它从没有导致预测误差的大的突变。

Ο它的预测误差具有更少量的高频噪音。这对模式匹配是个问题，对LMS更是个大问题。

当采用n₁＝5，n₂＝5时，并没有观察到明显的差异。不同的应用可以相应调整或比较几种情况后选取最佳值。

附图说明

图1是同步跟踪子系统的结构图

图2是人体胸部位置信号由于呼吸作用而随时间变化示意图

图3是模糊位置模型的模糊划区图

具体实施方式

下面给出本发明的一个较好的实施例，并予以详细描述，以使能更好的本发明的结构特征、功能，而不是用来限制本发明的权利要求保护范围。

如图1所示，本发明装置包括数据采集模块，计算机，功率放大器，驱动电路以及动作机构。数据采集器将采集到的前30000个点的位置，以及在该时刻对应的体内的肿瘤的位置数据存入计算机，根据位置参数计算出相应的速度参数(dx，dy，dz)并存入计算机：

${\mathrm{dx}}_{i+1}=\frac{x_{i+1}-x_i}{t_{i+1}-t_i}$

${\mathrm{dy}}_{i+1}=\frac{y_{i+1}-y_i}{t_{i+1}-t_i}$

${\mathrm{dz}}_{i+1}=\frac{z_{i+1}-z_i}{t_{i+1}-t_i};$

根据计算机内存储的位置参数(x，y，z)，找出最大和最小的位置(x_max，y_max，z_max)和(x_min，y_min，z_min)；根据计算机内存储的速度参数(dx，dy，dz)，找出最大和最小的速度(dx_max，dy_max，dz_max)和(dx_min，dy_min，dz_min)，根据模糊划分的特征以及已经记录的最大和最小值，把x分为5个子集，它们可以是非常小(Very Small)，小(Small)，零(Zero)，大(Big)，非常大(Very Big)；y和z也各被分成5个子集；同样的，根据模糊划分的特征以及已经记录的最大和最小值，把速度也分为5个子集，可以是负大(Negative Big)，负小(Negative Small)，零(Zero)，正小(Positive Small)和正大(Positive Big)。这样对于x一共就有5*5＝25个子模型，y，z同理。所有的子模型都可以用以下方程表示：

${\hat{y}}_m(k+p)=a_{m,0}y\left(k\right)+...+a_{m,r}y(k-r)+c_m$

是第m规则的模式输出，m的范围是在1~M，而M是方程的总数，p是预测步骤，r是预测阶数，可以选择2-15中任意一个数字，这里我们选5，a_m和c_m是参数。

Claims (4)

1.一种呼吸同步跟踪系统的位置信号的模糊预测方法，其特征在于包括以下步骤：

(1)在病人体表确定一体表位置P₀并计算出相对应的体内病灶的目标位置参数P(x，y，z)；

(2)采集该位置P(x，y，z)在一段时间内的位置参数(x，y，z)存入计算机；

(3)根据P(x，y，z)的位置参数计算出相应的速度参数(dx，dy，dz)并存入计算机：

${\mathrm{dx}}_{i+1}=\frac{x_{i+1}-x_i}{t_{i+1}-t_i}$

${\mathrm{dy}}_{i+1}=\frac{y_{i+1}-y_i}{t_{i+1}-t_i}$

${\mathrm{dz}}_{i+1}=\frac{z_{i+1}-z_i}{t_{i+1}-t_i};$

(4)根据计算机内存储的位置参数(x，y，z)，找出最大和最小的位置(x_max，y_max，z_max)和(x_min，y_min，z_min)；

(5)根据计算机内存储的速度参数(dx，dy，dz)，找出最大和最小的速度(dx_max，dy_max，dz_max)和(dx_min，dy_min，dz_min)；

(6)制定模糊规则：根据最大和最小的位置参数，把x模糊划分出n1个子集，这n1个子集的每一个子集再被根据最大和最小的速度参数模糊划分为n2个子集，得到M＝n1*n2个子模型，每个子模型都对应一个规则，该规则用线性方程

${\hat{x}}_m(k+p)=a_{m,0}x\left(k\right)+...+a_{m,r}x(k-r)+c_m$

表示，

是第m规则的模式输出，m的范围是在1～M，而M是规则的总数，p是预测步骤，r是预测阶数，a_m和c_m是在第m规则的参数；对于y，z作同样处理；

(7)根据各个子模型所对应的计算机内已经存储的已知的位置参数(x，y，z)、速度参数(dx，dy，dz)，采用梯度下降算法，训练出对应子模型的参数a_m和c_m；

(8)采集下一时间点体表目标的位置参数，并计算出该位置参数对应的体内病灶的位置参数，根据已经制定的模糊规则，把计算出的参数映射到相应的子模型里，得到输出

对相应映射的子模型的输出结果加权平均，即为所需要的预测结果，

$\hat{x}(k+p)=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_i\·{\hat{x}}_i(k+p)}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_i}$

$\hat{y}(k+p)=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_i\·{\hat{y}}_i(k+p)}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_i}$

$\hat{z}(k+p)=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_i\·{\hat{z}}_i(k+p)}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_i}$

μ_i是子模型的隶属度。

2.根据权利要求1所述的模糊规则，其特征在于：把x模糊划分出5个子集，这5个子集再被根据速度模糊划分为5个子集，得到M＝25个子模型。

3.根据权利要求1所述的线性方程，其特征在于所述预测阶数r为5。

4.一种用于放射治疗的同步跟踪定位系统，包括：数据采集装置，同步跟踪子系统，执行机构，数据采集装置将采集到的数据送到同步子系统，处理放大后传送到执行机构，其特征在于：所述同步跟踪子系统包括有内建有如权利要求1所述的模糊规则以及所述子模型的计算机。

Images