CN101637388B - 呼吸运动预测方法 - Google Patents

Abstract

本发明公开了一呼吸运动预测方法，包括以下步骤：(1)输入：状态特征集；(2)获取实时的呼吸信号f(t)及其状态特征R(t)；(3)通过概率似然模型建立状态特征R(t)与状态特征集R(k)的相似性约束条件；通过概率先验模型建立相邻状态特征R(t_i)与R(t_i+1)的连续性约束条件；(4)通过最大后验概率模型从候选集中筛选满足步骤(3)所述条件的特征元素，并预测呼吸信号；(5)输出延时Δt后的呼吸信号f(t+Δt)。本发明采用最大后验概率算法构建模型，充分利用呼吸的局部和全局特征建立概率模型，与传统的线性预测模型相比：该模型能够较好地预测真实呼吸运动情况，预测的平均误差较小、提高了呼吸运动预测的精度，进而提高了放射治疗效果。

Description

呼吸运动预测方法

技术领域

本发明涉及一种放射治疗的呼吸运动预测方法，尤其是涉及一种放射治疗中呼吸信号的跟踪与预测方法。

背景技术

胸腹部器官肿瘤的放疗很大程度受呼吸运动影响。治疗中，肿瘤组织随时间动态地变化。由此产生的不确定性问题成为当前放射治疗面对的主要问题，并且解决的方法和手段还远未完善。传统的基于适形调强技术的放疗模式应对上述情况时通常采用等中心位移、呼吸训练、门控放射等技术。采用图像引导放疗应对这一情况时，需要分析肿瘤随呼吸运动的变化情况。

当前所采用的肿瘤跟踪算法，主要有：线性估计算法(LP)、线性外推算法(LE)；其中，线性估计算法(LP)算法将当前时刻肿瘤位置X_t表示为前n+1个状态的线性组合，对于一套训练样本X₁，X₂，...X_n，系数a_k可以通过解一组线性方程来最小化均方误差，不易选择最优历史状态数，线性估计只适用于延时很小的***；线性外推算法(LE)，在假设信号变化维持一个速度常量的条件下，利用最近两个时刻的样本值估计信号的变化速度，只能在极短时间内将肿瘤的运动视为常量，因此，线性外推算法(LE)算法只适用于延时很小的***。

上述传统的预测和跟踪技术的主要问题是：没有利用训练样本和先验知识的条件下，进行对象的局部线性化假设或非正则化数学建模；在对象形态发生随机变化时，无法有效逼近目标，从而导致较大的计算误差。

发明内容

针对现有技术的缺陷，本发明提供一种呼吸运动预测方法，通过最大后验概率算法构建模型，充分利用了呼吸的局部和全局特征，提高了呼吸运动预测的精度。

针对上述技术问题的解决，本发明包括以下步骤：

(1)输入：状态特征集R(k)；

(2)获取治疗中t＝t_i时刻的呼吸信号f(t)及其状态特征R(t)；

(3)通过概率似然模型建立状态特征R(t)与状态特征集R(k)的相似性约束条件；通过概率先验模型建立相邻状态特征R(t_i)与R(t_i+1)的连续性约束条件。

(4)通过最大后验概率模型从候选集中筛选满足步骤(3)所述条件的特征元素，并预测呼吸信号。

(5)输出延时Δt后的呼吸信号f(t+Δt)。

所述步骤(1)中状态特征集R(k)取得的具体步骤为：

(a)获取某个对象的一个周期T的正常呼吸信号f(k)，k＝k₁，k₂，...k_T；

(b)提取该周期T下的状态特征集R(k)，k＝k₁，k₂，...k_T；

(c)状态特征集R(k)包含了呼吸信号的方向、速率、幅度、一阶导数过零点，以及呼吸信号周期等元素。

所述步骤(3)中概率似然模型为P(Y/X)；

其中：X代表某一时刻t的状态特征R(t)中的元素，

Y代表状态特征集R(k)中对应状态的相应元素。

似然模型P(Y/X)给出了相应元素之间的相似度。

令Δt＝t_i-t_i-1，则某一t_i时刻和t_i-1时刻呼吸运动的相似性似然模型分别为：

$P\left(Y_{k_j}\right|X_{t_i})=K\·\exp (-\frac{{(r_{t_i}^x-r_{k_j}^y)}^2}{{2\mathrm{\σ}}_a^2})\·\exp (-\frac{{({\mathrm{\Δr}}_{t_i}^x-{\mathrm{\Δr}}_{k_j}^y)}^2}{{2\mathrm{\σ}}_b^2})$

$P\left(Y_{k_{j-1}}\right|X_{t_{i-1}})=K\·\exp (-\frac{{(r_{t_{i-1}}^x-r_{k_{j-1}}^y)}^2}{{2\mathrm{\σ}}_a^2})\·\exp (-\frac{{({\mathrm{\Δr}}_{t_{i-1}}^x-{\mathrm{\Δr}}_{k_{j-1}}^y)}^2}{{2\mathrm{\σ}}_b^2})$

其中：

分别为状态特征R(t_i)中的呼吸信号幅度和速率，

分别为状态特征集R(k_j)中对应的幅度和速率，

σ_a·、σ_b分别为幅度方差和速率方差，

K为与幅度方差和速率方差相关的常数。

所述步骤(3)中概率先验模型为P(X)，它给出了相邻时刻的状态特征中相应元素之间的相似度。t_i时刻呼吸运动的先验模型为：

$P\left(X_{t_i}\right)=K\·\exp (-\frac{{(r_{t_i}^x-r_{t_{i-1}}^x)}^2}{{2\mathrm{\σ}}_a^2})\·\exp (-\frac{{({\mathrm{\Δr}}_{t_i}^x-{\mathrm{\Δr}}_{t_{i-1}}^x)}^2}{{2\mathrm{\σ}}_b^2})$

其中：

为状态特征R(t_i)中两个相邻时刻的呼吸信号幅度，

为状态特征R(t_i)中两个相邻时刻的呼吸信号速率，

σ_a·、σ_b分别为幅度方差和速率方差，

K为与幅度方差和速率方差相关的常数。

所述步骤(4)中最大后验概率模型为 $P\left(X_{t_i}\right|Y_{k_j})\∝P\left(X_{t_i}\right)\·P\left(Y_{k_j}\right|X_{t_i})P\left(Y_{k_{j-1}}\right|X_{t_{i-1}}),$

它将上述两类相似度值相乘，并在“相应元素”的候选集中选择乘积最大的元素作为预测的结果输出。

所述实时呼吸信号f(t)的状态特征R(t)的获取手段包括：

(a)通过测量体表标记物相对静态参照物的位移信号变化；

(b)呼吸条件下，通过呼吸测量设备获取呼吸信号；

(c)通过计算上述(a)或(b)中呼吸信号的方向、速率、幅度、一阶导数过零点提取状态特征。

本发明的比较试验如下：1.采集11个对象在特定位置处的呼吸波信号，统一采用30Hz的采样频率分析计算相同的患者对象的呼吸信号；2.针对某一样本数据，利用相应宽度Q_w的高斯算子去除噪声；再对每一个体分析呼吸特征相似性、计算量化的呼吸状态特征集，并测量各状态阶段的动态参数；3.预测实验中，三种传统算法模型：如线性估计算法(LP)、线性外推算法(LE)和卡尔曼滤波算法(KF)，和本发明提出的最大后验概率模型(MAP)比较，并统一利用呼吸状态判别模型判断呼吸状态、引导呼吸跟踪模型。

针对每一个呼吸对象进行运动预测时，采用若干周期的数据进行试验统计，结果如表1所示：

表1左边4栏，统计了患者的个体化呼吸动态参数测量结果，其中r_L ^ex和r_R ^ex为呼气时段的首尾幅度特征，r_L ^eoe和r_R ^eoe为呼气末时段的首尾幅度特征，r_L ⁱⁿ和r_R ^ex为吸气时段的首尾幅度特征。

表1右边5栏，分别统计了无估计和利用不同算法的误差绝对值。MAP方法与真实结果之间的误差平均值为，远小于传统方法估计结果与真实结果之间的误差平均值。

表1.个体呼吸特征与呼吸预测的统计误差比较

本发明采用最大后验概率算法构建模型时，采集患者呼吸信号的曲线特征，通过最大后验概率算法在状态特征集中寻找最佳的局部信号模式与当前信号非线性匹配。提取呼吸曲线的局部特征时，用到呼吸的速率、幅度、过程相关性等标量和方向矢量的统计测量信息，充分利用呼吸的局部和全局特征建立概率模型，与传统的线性预测模型相比：：该模型能够较好地预测真实呼吸运动情况，预测的平均误差较小、提高了呼吸运动预测的精度，进而提高了放射治疗效果。

附图说明

图1是本发明的流程框图

图2是呼吸运动预测示意图

图3是呼吸运动状态判别示意图

具体实施方式

下面结合附图对本发明作进一步说明。

实施基于预测的图像引导跟踪的步骤如下：

1.如图1中的10：输入周期T下的呼吸信号样本f(k)，k＝k₁，k₂，...k_T，治疗前，一个标准呼吸周期T下：获取一个呼吸周期中的呼吸信号样本f(k)。并对该信号进行去除噪声处理。一般个体对象的正常呼吸过程伴有信号噪声，噪声信号平均周期为Q_w，用高斯算子 $G_{\mathrm{\σ}}\left(x\right)={\left(\sqrt{2\mathrm{\π}}\mathrm{\σ}\right)}^{-1}{e}^{x^2/{2\mathrm{\σ}}^2}$ 与呼吸信号卷积，并选择脉宽σ＝Q_w，可以有效去除噪声。

2.如图1中的20：提取状态特征集R(k)，k＝k₁，k₂，...k_T；治疗前，一个标准呼吸周期下，提取呼吸状态特征集R(k)，k＝k₁，k₂，...k_T，状态特征集R(k)包含了呼吸信号的方向、速率、幅度、一阶导数过零点，以及呼吸信号周期等元素。

3.如图1中的30：采集实时呼吸信号f(t)，t＝t₁，t₂，…t_i，…治疗中，获取治疗中t＝t_i时刻的呼吸信号f(t_i)，i＝1，2，...，

4.如图1中的40：提取实时状态特征R(t)，t＝t₁，t₂，…t_i，…针对治疗中实时呼吸信号f(t)，提取实时状态特征R(t)，状态特征R(t)包含了呼吸信号的方向、速率、幅度、一阶导数过零点，以及呼吸信号周期等元素。

5.如图1中的50：判断呼吸状态是否正常

本发明选择正常的呼吸状态进行判别，原因有二：第一，非正常呼吸状态如：咳嗽、哮喘、和情绪波动等，会引起患者体内的肿瘤形态的气质性变化，肿瘤运动与呼吸的相关性会变差。因此，判别呼吸状态可以提高预测的精度。第二，在相同的信号采样密度和延迟条件下，非正常呼吸会导致预测的瞬时误差增加，影响瞬时治疗效果。

呼吸状态判别建立在如下特征计算基础上：针对某一个体的呼吸过程，如果用函数 $\mathrm{\&rho;}(x_{t_i},x_{t_{i+1}}),(0<\mathrm{\&rho;}<1)$ 表示相关样点

和

之间特征矢量的相似性度量，那么特征矢量间的相似性可以用各元素间相似度的加权和表示为：

其中w_r，w_v代表幅度、速度的相关系数，上式反映了呼吸过程的局部动态相似性，只有当关联点处在相邻或近似相位条件下，相似度测量才有意义。ρ_r和ρ_v分别代表关联点之间的幅度、速度的相似性，表示为：

${\mathrm{\&rho;}}_r(x_{t_i},x_{t_{i+1}})=1-\frac{|r_x-r_y|}{r_{\max }};{\mathrm{\&rho;}}_v(x_{t_i},x_{t_{i+1}})=1-\frac{|{\mathrm{\&Delta;r}}_x-{\mathrm{\&Delta;r}}_y|}{{\mathrm{\&Delta;r}}_{\max }}$

针对特定的个体呼吸和采样频率，r_max和|Δr|_max可看作常数值。这样以来，呼吸的状态变化可以通过测量表1的呼吸特征局部参数进行判别。如表1中，如果用r^max和r^min分别代表三种呼吸状态中幅度的最大值和最小值，那么表中的参数分布范围可在个体对象的观察过程中测量得到。

表1.呼吸特征的局部参数测量依据

6.如图1中的60：概率似然模型P(Y|X)

呼吸预测过程中，似然概率模型用于测量当前局部特征与训练样本集中对应特征的相似性。令采样间隔Δt＝t_i-t_i-1，在t_i和t_i-1时刻的呼吸运动与样本信号之间的相似度可以用概率的似然模型分别表示为：

$P\left(Y_{k_j}\right|X_{t_i})=K\&CenterDot;\exp (-\frac{{(r_{t_i}^x-r_{k_j}^y)}^2}{{2\mathrm{\&sigma;}}_a^2})\&CenterDot;\exp (-\frac{{({\mathrm{\&Delta;r}}_{t_i}^x-{\mathrm{\&Delta;r}}_{k_j}^y)}^2}{{2\mathrm{\&sigma;}}_b^2})$

$P\left(Y_{k_{j-1}}\right|X_{t_{i-1}})=K\&CenterDot;\exp (-\frac{{(r_{t_{i-1}}^x-r_{k_{j-1}}^y)}^2}{{2\mathrm{\&sigma;}}_a^2})\&CenterDot;\exp (-\frac{{({\mathrm{\&Delta;r}}_{t_{i-1}}^x-{\mathrm{\&Delta;r}}_{k_{j-1}}^y)}^2}{{2\mathrm{\&sigma;}}_b^2})$

其中：

分别为状态特征R(t_i)中的呼吸信号幅度和速率，

分别为对应特征R(k_j)中的幅度和速率，

σ_a·σ_b分别为幅度方差和速率方差，

K为与幅度方差和速率方差相关的常数。

以上概率公式描述了在t_i和t_i-1时刻当前信号与观测样本信号的相似度，概率越大、相似度越强。在实际呼吸的每一个状态阶段按照采样顺序进行运动估计时，通过最大化上述概率的乘积可以使得估计的结果向平滑性、实际与测量的相似性逼近。

7.如图1中的70：概率先验模型P(X)

先验约束条件取自当前呼吸过程相邻时刻的连续性和相关性。令采样间隔Δt＝t_i-t_i-1，且t时刻被估计信号幅度为r_t ^x，变化速度为 ${\mathrm{\&Delta;r}}_t^x=r_t^x-r_{t-\mathrm{\&Delta;t}}^x$ 则t时刻呼吸运动的先验模型为：

$P\left(X_{t_i}\right)=K\&CenterDot;\exp (-\frac{{(r_{t_i}^x-r_{t_{i-1}}^x)}^2}{{2\mathrm{\&sigma;}}_a^2})\&CenterDot;\exp (-\frac{{({\mathrm{\&Delta;r}}_{t_i}^x-{\mathrm{\&Delta;r}}_{t_{i-1}}^x)}^2}{{2\mathrm{\&sigma;}}_b^2})$

其中：

为状态特征R(t_i)中两个相邻时刻的呼吸信号幅度，

为状态特征R(t_i)中两个相邻时刻的呼吸信号速率，

σ_a·、σ_b分别为幅度方差和速率方差，

K为与幅度方差和速率方差相关的常数。

8.如图1中的80：最大后验概率模型

依据贝叶斯理论：后验概率为先验和似然概率之乘积P(X|Y)＝P(X)P(Y|X)。本发明的似然模型P(Y|X)兼顾了时刻t_i和t_i-1下实时信号和观测样本的相似性，因此，最大后验概率模型表示为 $P\left(X_{t_i}\right|Y_{k_j})\&Proportional;P\left(X_{t_i}\right)\&CenterDot;P\left(Y_{k_j}\right|X_{t_i})P\left(Y_{k_j}\right|X_{t_i}),$ 它将上述两类相似度相乘，并在“相应元素”候选集中选择乘积最大的元素作为预测的结果输出。

9.呼吸预测过程中的状态判别与控制

如图2、图3所示，本发明采用了一种概率跟踪与状态控制策略，并利用有限状态模型判别呼吸状态并控制跟踪的过程。状态的特性由加权和函数

量化表示，并将呼吸状态分为四种模式：呼气EX、呼气末EOE、吸气IN、不规则状态IRR。控制过程，有限状态模型依据状态转移条件R_EOE、R_IN、R_EX、R_IRR判断呼吸各阶段的状态变化。该条件又依据了函数

和表1的呼吸局部参数统计结果。当状态参数属于事先确定的参数范围时，转移条件起作用，并完成从一个状态向对应状态的跳变。当不满足跳变条件的情况下，则继续利用概率估计模型跟踪当前状态。具体执行策略如下：

(1)从吸气末端的y₁点(如图2(a)所示)开始跟踪，起始点

作为初始种子点注册下来，有限状态模型记录当前状态类别为EX。

(2)依据当前点

的特征参数和概率候选集{Y}_EX，用概率模型预测下一点的值为

如果

仍属于样本集{Y}_EX，则继续进行当前状态下的概率跟踪。

(3)当估计到t_i+1时刻，如果出现 ${\tilde{x}}_{t_{i+1}}\&NotElement;{\left\{Y\right\}}_{\mathrm{EX}}$ 且 ${\tilde{x}}_{t_{i+1}}\&Element;{\left\{Y\right\}}_{\mathrm{EOE}},$ 则跳转到相应状态EOE，并且后续概率候选集使用{Y}_EOE；鉴于状态变化的规则，如果 ${\tilde{x}}_{t_{i+1}}\&Element;{\left\{Y\right\}}_{\mathrm{IN}},$ 则跳转到IRR。

以上执行策略中，状态判别结果用来控制概率预测。施加状态判别作用后，前述后验概率预测模型表示为 $P_c\left(X_t\right|Y_k)=T(d_{t_i},d_{t_{i-1}})\&CenterDot;P\left(X_t\right|Y_k),$ 其中

为开关函数，当相邻时刻t_i和t_i-1的呼吸状态方向

和

相同时 $T(d_{t_i},d_{t_{i-1}})=1,$ 否则为0。

10.如图1中的90：输出延时Δt后的呼吸信号f(t+Δt)

后验概率预测是一种给予当前时刻信号f(t)的特征R(t_i)。本发明完成在已知均匀采样的频率和平均延时Δt＝t_i-t_i-1条件下的状态控制和预测，并通过预测t_i+1时刻的状态特征R(t_i+1)，i＝1，2，...描述延时Δt后的呼吸信号f(t_i+1)。

Claims (9)

1.一种呼吸运动预测方法，其特征在于，包括以下内容和步骤：

(1)输入：状态特征集R(k)；

(2)获取治疗中t＝t_i时刻的呼吸信号f(t)及其状态特征R(t)；

(3)通过概率似然模型建立状态特征R(t)与状态特征集R(k)的相似性约束条件；通过概率先验模型建立相邻状态特征R(t_i)与R(t_i+1)的连续性约束条件；

(4)通过最大后验概率模型从候选集中筛选同时满足步骤(3)中两个条件的特征元素，并预测呼吸信号；

(5)输出延时Δt后的呼吸信号f(t+Δt)。

2.如权利要求1所述的呼吸运动预测方法，其特征在于，所述步骤(1)中状态特征集R(k)取得的具体步骤为：

(a)获取一个呼吸周期T的标准呼吸信号f(k)，k＝k₁，k₂，...k_T；

(b)提取该周期T下的状态特征集R(k)，k＝k₁，k₂，...k_T；

(c)状态特征集R(k)包含的元素为：呼吸信号的方向、速率、幅度、一阶导数过零点、及信号周期。

3.如权利要求1所述的呼吸运动预测方法，其特征在于，所述步骤(3)中概率似然模型为P(Y/X)，其中X代表当前时刻状态特征R(t)中的元素，Y代表状态特征集R(k)中对应状态的相应元素；概率似然模型P(Y/X)给出了“相应元素”之间的相似度。

4.如权利要求3所述的呼吸运动预测方法，其特征在于，某一t_i时刻和t_i-1时刻呼吸运动的概率似然模型分别为：

$P\left(Y_{k_j}\right|X_{t_i})=K\&CenterDot;\exp (-\frac{{(r_{t_i}^x-r_{k_j}^y)}^2}{2{\mathrm{\&sigma;}}_a^2})\&CenterDot;\exp (-\frac{{(\mathrm{\&Delta;}{r}_{t_i}^x-\mathrm{\&Delta;}{r}_{k_j}^y)}^2}{2{\mathrm{\&sigma;}}_b^2})$

$P\left(Y_{k_{j-1}}\right|X_{t_{i-1}})=K\&CenterDot;\exp (-\frac{{(r_{t_{i-1}}^x-r_{k_{j-1}}^y)}^2}{2{\mathrm{\&sigma;}}_a^2})\&CenterDot;\exp \left(\frac{{(\mathrm{\&Delta;}{r}_{t_{i-1}}^x-\mathrm{\&Delta;}{r}_{k_{j-1}}^y)}^2}{2{\mathrm{\&sigma;}}_b^2}\right)$

其中：

分别为状态特征R(t_i)中的呼吸信号幅度和速率，

分别为状态特征集R(k_j)中呼吸信号的幅度和速率，

σ_a、σ_b分别为幅度方差和速率方差，

K为与幅度方差和速率方差相关的常数。

5.如权利要求1所述的呼吸运动预测方法，其特征在于，所述步骤(3)中概率先验模型为P(X)，它给出了相邻时刻的状态特征中相应元素之间的相似度。

6.如权利要求5所述的呼吸运动预测方法，其特征在于，t_i时刻呼吸运动的先验模型为：

$P\left(X_{t_i}\right)=K\&CenterDot;\exp (-\frac{{(r_{t_i}^x-r_{t_{i-1}}^x)}^2}{2{\mathrm{\&sigma;}}_a^2})\&CenterDot;\exp (-\frac{{(\mathrm{\&Delta;}{r}_{t_i}^x-\mathrm{\&Delta;}{r}_{t_{i-1}}^x)}^2}{2{\mathrm{\&sigma;}}_b^2})$

其中：

为状态特征R(t_i)中两个相邻时刻的呼吸信号幅度，

为状态特征R(t_i)中两个相邻时刻的呼吸信号速率，

σ_a、σ_b分别为幅度方差和速率方差，

K为与幅度方差和速率方差相关的常数。

7.如权利要求1所述的呼吸运动预测方法，其特征在于，所述步骤(4)中最大后验概率模型为P(X/Y)＝P(X)P(Y/X)，它将上述概率似然模型和概率先验模型计算出的数值相乘，并在相应元素的候选集中选择乘积最大的元素作为预测的结果输出。

8.如权利要求7所述的呼吸运动预测方法，其特征在于，最大后验概率模型为：

9.如权利要求1所述的呼吸运动预测方法，其特征在于，所述状态特征R(t)的获取手段包括：

(a)通过测量体表标记物相对静态参照物的位移信号变化，提取状态特征；

(b)呼吸条件下，通过呼吸测量设备获取呼吸信号，提取状态特征。

Images