CN110675960B - 基于时间序列分析的呼吸运动预测建模方法 - Google Patents

Abstract

本发明提供一种基于时间序列分析的呼吸运动预测建模方法，属于生物医学工程领域。本发明分别基于自回归积分滑动平均ARIMA模型、基于STL分解和指数平滑ETS方法和基于神经网络自回归NNAR模型计算得到呼吸运动的预测结果，定义基于RMSE和信息熵的模型评价指标，可以选择出同时具有“小误差”和“高稳定性”特征的预测模型参数，提高了预测方法的科学性。与已有技术相比大大提高了呼吸运动预测的时间长度，将以往通常只有几百毫秒的短期预测模型，升级到可以预测未来数十秒内呼吸运动幅度的长期预测模型。

Description

基于时间序列分析的呼吸运动预测建模方法

技术领域

本发明涉及肿瘤医学、核医学及生物医学工程领域，尤其涉及一种基于时间序列分析的呼吸运动预测模型的构建方法。

背景技术

放射治疗作为治疗恶性肿瘤的一项重要技术，在过去几十年得到了突飞猛进地发展。近年来，随着质子治疗技术的普及，肿瘤靶向放疗的概念被提出，放疗的精准度得到了大大提升，放疗过程中对正常器官和组织的伤害也降到较低水平。然而，当肿瘤位置位于胸腔或腹腔时，随着病人呼吸运动，肿瘤位置也随之周期性移动，导致明显的临床靶向误差并严重影响适形放射治疗的效率。此外，放射治疗系统本身的延迟也会造成照射误差，影响治疗效果。

为了解决这一问题，临床上通常采用靶区外延、呼吸门控、呼吸训练、屏气等方法，取得了一定的效果。然而，这些方法均为被动的误差补偿方式，分别有其自身的局限性。近年来，很多学者提出了呼吸波形预测的算法和数学模型，可以主动预测病人未来的呼吸运动波形，并根据预测结果及时调整照射靶区位置。精准的呼吸运动预测模型一方面可以降低病人呼吸运动带来的误差，另一方面可以补偿放疗设备自身的通信延迟。

常用的呼吸运动预测方法主要有线性外推法、卡尔曼滤波法、有限状态模型、人工神经网络、支持向量机等，这些方法对呼吸运动的短期预测卓有成效，在数百毫秒的预测尺度上的均方根误差(Root Mean Squared Error,RMSE)普遍小于1mm。然而，由于缺乏在长时间尺度上的呼吸运动预测手段，使得这些预测模型在一些新型质子治疗技术(如激光驱动质子放射治疗)难以得到应用。

发明内容

根据现有技术方案中无法实现长时间尺度上的呼吸运动预测的问题，本发明提供了一种基于时间序列分析的呼吸运动预测建模方法。

为了实现上述目的，本发明提供一种基于时间序列分析的呼吸运动预测建模方法，其具体实施步骤如下：

1)使用高精度激光测距仪获得人体表面的呼吸运动离散信号{Y_t}；

2)基于自回归积分滑动平均(ARIMA)模型计算得到呼吸运动的预测结果Y_ARIMA；

3)基于STL(Seasonal-Trend decompositionprocedure basedon Loess)分解和指数平滑(ETS)方法计算得到呼吸运动的预测结果Y_STL-ETS；

4)基于神经网络自回归(NNAR)模型计算得到呼吸运动的预测结果Y_NNAR；

5)基于上述三种方法得到的预测结果，建立等权重组合预测模型Y_combination；

6)定义基于RMSE和信息熵的模型评价指标；

7)根据步骤(6)中提出的模型评价指标，对步骤(5)中的组合预测模型进行参数优化，得到最优的组合预测模型。

其中，所述步骤2)中使用自回归积分滑动平均(ARIMA)模型计算得到呼吸运动的具体步骤为：

(a){Y_t}表示一个离散的呼吸运动时间序列,{ε_t}为白噪声，ARMA(p,q)可以写成如下形式:

Y_t＝φ₁Y_t-1+φ₂Y_t-2+…+φ_pY_t-p+ε_t+θ₁ε_t-1+…+θ_qε_t-q，

或者，φ(B)Y_t＝θ(B)ε_t。其中φ(B)＝1-φ₁B-…-φ_pB^p，θ(B)＝1+θ₁B+…+θ_qB^q，B为滞后算子。

(b)原始的呼吸运动序列并不是平稳的时间序列，对其进行一阶差分处理。差分后的时间序列为W_t＝▽^dY_t＝(1-B)^dY_t，该序列通常是一个平稳时间序列过程,d表示差分的阶数。

(c)因此，ARIMA(p,d,q)可以表示为

ARIMA(p,d,q)的阶数可以由自相关函数和偏自相关函数粗略计算，本发明中使用贝叶斯信息准则(Bayesian Information Criterion,BIC)确定模型的阶数。BIC的定义为：

BIC＝kln(n)-2ln(L)。

其中，k为模型的参数个数，L为似然函数，n为样本数目。

所述步骤3)中使用STL分解和指数平滑(ETS)方法计算呼吸运动的预测结果Y_STL-ETS的具体步骤如下：

(a)对于一个非平稳的呼吸运动时间序列{Y_t},对其进行STL分解的可加模型为：

Y_t＝S_t+T_t+R_t.

其中S_t是季节项，T_t是趋势项，R_t是剩余项。

(b)对于分解后的呼吸运动时间序列，其中的呼吸趋势项使用指数平滑模型。设{y_t}为呼吸运动时间序列分解后一个有长期趋势的项，其下一个时刻的呼吸预测值计算方法为：

其中0≤α≤1是平滑参数,

代表第t+1时刻的预测值。

所述步骤4)中使用神经网络自回归(NNAR)模型计算得到呼吸运动预测结果的模型具体如下：

(a)将呼吸运动时间序列作为神经网络的输入元素(x_i)，则神经网络的隐藏层的中间值(y_j)可以由输入元素(x_i)和偏置项通过激活函数计算得到，激活函数表达为：

其中，n_i为输出层输入元素(x_i)的个数，w_ij为输入层与隐藏层之间的权重。

(b)根据计算得到的隐藏层数值，在输出层计算呼吸运动的预测结果的公式为：

其中z_k表示使用神经网络自回归(NNAR)模型计算得到的呼吸运动预测值,权重w_jk表示隐藏层与输出层之间的权重值。

所述步骤5)中提到的权重组合预测模型Y_combination可以表示为：

Y_combination＝w_ARIMA·Y_ARIMA+w_STL-ETS·Y_STL-ETS+(1-w_ARIMA-w_STL-ETS)·Y_NNAR.

所述步骤6)中提到了基于RMSE和信息熵的模型评价指标，该评价指标的定义过程如下：

(a)定义E＝{e₁,e₂,…,e_t,…,e_m}为呼吸运动预测模型的误差序列，这里e_t表示在时刻t的预测误差。

(b)计算呼吸运动预测模型的误差序列在时刻t的相对误差比例，其计算公式为：

(c)计算预测误差序列的信息熵，其定义为：

其中k＝1/lnm,表示归一化参数。

(d)呼吸运动预测模型的评价指标定义为：

F＝RMSE·(1-Ent(E)).

其中的RMSE表示均方根偏差。

所述步骤(7)中提到的优化方法，其具体定义如下：

其中min表示优化模型的目标函数是要求F最小，s.t.表示优化模型的约束条件。通过参数扫描的方法，根据上述优化模型，优化步骤(5)中提到的组合预测模型的参数，得到最终的预测模型。

本发明的技术效果：

本发明与已有技术相比最大的突破是大大提高了呼吸运动预测的时间长度，将以往通常只有几百毫秒的短期预测模型，升级到可以预测未来数十秒内呼吸运动幅度的长期预测模型。此外，基于本发明中提出的模型评价指标，可以选择出同时具有“小误差”和“高稳定性”特征的预测模型参数，提高了预测方法的科学性。

附图说明

图1是呼吸运动预测模型的预测结果；其中，(a)为基于自回归积分滑动平均(ARIMA)模型计算得到呼吸运动的预测结果Y_ARIMA；(b)为基于STL(Seasonal-Trenddecomposition procedure based on Loess)分解和指数平滑(ETS)方法计算得到呼吸运动的预测结果Y_STL-ETS；(c)基于神经网络自回归(NNAR)模型计算得到呼吸运动的预测结果Y_NNAR；

图2是不同呼吸运动预测模型的均方根偏差比较图。

具体实施方式

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

实施基于时间序列分析的呼吸运动长期预测方法步骤如下：

1.使用高精度激光测距仪获得人体表面的呼吸运动离散信号，共选取137.25秒的离散数据进行建模和预测，其中前122.25秒的数据作为训练集，剩余数据作为测试集。

2.根据BIC最小的原则，基于自回归积分滑动平均(ARIMA)模型的参数为ARIMA(18,2,1)，具体步骤为：

(a){Y_t}表示一个离散的呼吸运动时间序列,{ε_t}为白噪声，ARMA(p,q)可以写成如下形式:

Y_t＝φ₁Y_t-1+φ₂Y_t-2+…+φ_pY_t-p+ε_t+θ₁ε_t-1+…+θ_qε_t-q，

或者，φ(B)Y_t＝θ(B)ε_t。其中φ(B)＝1-φ₁B-…-φ_pB^p，θ(B)＝1+θ₁B+…+θ_qB^q，B为滞后算子。

(b)原始的呼吸运动序列并不是平稳的时间序列，对其进行一阶差分处理。差分后的时间序列为W_t＝▽^dY_t＝(1-B)^dY_t，该序列通常是一个平稳时间序列过程,d表示差分的阶数。

(c)因此，ARIMA(p,d,q)可以表示为

ARIMA(p,d,q)的阶数可以由自相关函数和偏自相关函数粗略计算，本发明中使用贝叶斯信息准则(Bayesian Information Criterion,BIC)确定模型的阶数。BIC的定义为：

BIC＝kln(n)-2ln(L)。

其中，k为模型的参数个数，L为似然函数，n为样本数目。其预测结果如图1中(a)所示，RMSE为1.3239mm。

3.基于STL(Seasonal-Trend decomposition procedure based on Loess)分解和指数平滑(ETS)方法计算，具体步骤如下：

(a)对于一个非平稳的呼吸运动时间序列{Y_t},对其进行STL分解的可加模型为：

Y_t＝S_t+T_t+R_t.

其中S_t是季节项，T_t是趋势项，R_t是剩余项。

(b)对于分解后的呼吸运动时间序列，其中的呼吸趋势项使用指数平滑模型。设{y_t}为呼吸运动时间序列分解后一个有长期趋势的项，其下一个时刻的呼吸预测值计算方法为：

其中0≤α≤1是平滑参数,

代表第t+1时刻的预测值。

得到呼吸运动的预测结果如图1中(b)所示，RMSE为2.8339mm。

4.基于神经网络自回归(NNAR)模型计算得到呼吸运动的预测模型为NNAR(25,1,13)₂₁，

具体步骤如下：

(a)将呼吸运动时间序列作为神经网络的输入元素(x_i)，则神经网络的隐藏层的中间值(y_j)可以由输入元素(x_i)和偏置项通过激活函数计算得到，激活函数表达为：

其中，n_i为输出层输入元素(x_i)的个数，w_ij为输入层与隐藏层之间的权重。

(c)根据计算得到的隐藏层数值，在输出层计算呼吸运动的预测结果的公式为：

其中z_k表示使用神经网络自回归(NNAR)模型计算得到的呼吸运动预测值,权重w_jk表示隐藏层与输出层之间的权重值。其预测结果如图1中(c)所示，RMSE为1.4374mm。

5.基于上述三种方法得到的预测结果，建立等权重组合预测模型Y_combination，权重组合预测模型Y_combination可以表示为：

Y_combination＝w_ARIMA·Y_ARIMA+w_STL-ETS·Y_STL-ETS+(1-w_ARIMA-w_STL-ETS)·Y_NNAR.

RMSE为1.1190mm。

6.根据优化模型：

(a)定义E＝{e₁,e₂,…,e_t,…,e_m}为呼吸运动预测模型的误差序列，这里e_t表示在时刻t的预测误差。

(b)计算呼吸运动预测模型的误差序列在时刻t的相对误差比例，其计算公式为：

(c)计算预测误差序列的信息熵，其定义为：

其中k＝1/lnm,表示归一化参数。

(d)呼吸运动预测模型的评价指标定义为：

F＝RMSE·(1-Ent(E)).

其中的RMSE表示均方根偏差。

7.通过参数扫描的方法优化模型，其具体定义如下：

其中min表示优化模型的目标函数是要求F最小，s.t.表示优化模型的约束条件。计算出F最小时的模型参数，分别为w_ARIMA＝0.116，w_STL-ETS＝0.265。根据此参数建立参数优化后的组合预测模型，在15秒的预测时间长度上均方根偏差(RMSE)为0.9745mm。

优化步骤5中提到的组合预测模型的参数，得到最终的预测结果。三种单独的预测模型和两种组合预测模型的均方根偏差对比如图2所示。其中，在三种单独的预测模型中，自回归积分滑动平均(ARIMA)模型的预测效果最好，均方根偏差(RMSE)是最小的；而使用STL分解-指数平滑模型预测得到的结果最差。两种组合预测模型给出的预测结果均好于三种单独预测模型的预测结果，其中参数优化后的预测模型所得到的均方根偏差比较稳定，且均小于1毫米。

虽然本发明已以较佳实施例披露如上，然而并非用以限定本发明。任何熟悉本领域的技术人员，在不脱离本发明技术方案范围情况下，都可利用上述揭示的方法和技术内容对本发明技术方案作出许多可能的变动和修饰，或修改为等同变化的等效实施例。因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本发明技术方案保护的范围内。

Claims (5)

1.一种基于时间序列分析的呼吸运动预测建模方法，其具体实施步骤如下：

A.使用高精度激光测距仪获得人体表面的呼吸运动离散信号{Y_t}；基于自回归积分滑动平均ARIMA模型计算得到呼吸运动的预测结果Y_ARIMA；基于STL分解和指数平滑ETS方法计算得到呼吸运动的预测结果Y_STL-ETS；基于神经网络自回归NNAR模型计算得到呼吸运动的预测结果Y_NNAR；基于上述三种方法得到的预测结果，建立等权重组合预测模型Y_combination；

权重组合预测模型Y_combination表示为：

Y_combination＝w_ARIMA·Y_ARIMA+w_STL-ETS·Y_STL-ETS+(1-w_ARIMA-w_STL-ETS)·Y_NNAR，

B.定义基于RMSE和信息熵的模型评价指标；该评价指标的定义过程如下：

1)定义E＝{e₁,e₂,…,e_t,…,e_m}为呼吸运动预测模型的误差序列，e_t表示在时刻t的预测误差；

2)计算呼吸运动预测模型的误差序列在时刻t的相对误差比例，其计算公式为：

3)计算预测误差序列的信息熵，其定义为：

其中k＝1/lnm,表示归一化参数；

C.呼吸运动预测模型的评价指标定义为：

F＝RMSE·(1-Ent(E)).

其中的RMSE表示均方根偏差；

D.对上述权重组合预测模型进行参数优化，优化方法具体定义如下：

其中min表示优化模型的目标函数是要求F最小，s.t.表示优化模型的约束条件，通过参数扫描的方法，得到最终的预测模型。

2.如权利要求1所述基于时间序列分析的呼吸运动预测建模方法，其特征在于，使用自回归积分滑动平均ARIMA模型计算得到呼吸运动的具体步骤为：

2-1){Y_t}表示一个离散的呼吸运动时间序列,{ε_t}为白噪声，ARMA(p,q)写成如下形式:

Y_t＝φ₁Y_t-1+φ₂Y_t-2+…+φ_pY_t-p+ε_t+θ₁ε_t-1+…+θ_qε_t-q，

φ(B)＝1-φ₁B-…-φ_pB^p，θ(B)＝1+θ₁B+…+θ_qB^q，B为滞后算子；

2-2)原始的呼吸运动序列并不是平稳的时间序列，对其进行一阶差分处理，差分后的时间序列为

该序列是一个平稳时间序列过程,d表示差分的阶数；

2-3)ARIMA(p,d,q)表示为

3.如权利要求2所述基于时间序列分析的呼吸运动预测建模方法，其特征在于，使用贝叶斯信息准则BIC确定模型的阶数，BIC的定义为：

BIC＝kln(n)-2ln(L)，

其中，k为模型的参数个数，L为似然函数，n为样本数目。

4.如权利要求1所述基于时间序列分析的呼吸运动预测建模方法，其特征在于，使用STL分解和指数平滑ETS方法计算呼吸运动的预测结果Y_STL-ETS的具体步骤如下：

3-1)对于一个非平稳的呼吸运动时间序列{Y_t},对其进行STL分解的可加模型为：

Y_t＝S_t+T_t+R_t；

其中S_t是季节项，T_t是趋势项，R_t是剩余项；

3-2)对于分解后的呼吸运动时间序列，其中的呼吸趋势项使用指数平滑模型，设{y_t}为呼吸运动时间序列分解后一个有长期趋势的项，其下一个时刻的呼吸预测值计算方法为：

其中0≤α≤1是平滑参数,

代表第t+1时刻的预测值。

5.如权利要求1所述基于时间序列分析的呼吸运动预测建模方法，其特征在于，使用神经网络自回归NNAR模型计算得到呼吸运动预测结果的模型具体如下：

4-1)将呼吸运动时间序列作为神经网络的输入元素x_i，则神经网络的隐藏层的中间值y_j由输入元素x_i和偏置项通过激活函数计算得到，激活函数表达为：

其中，n_i为输出层输入元素x_i的个数，w_ij为输入层与隐藏层之间的权重；

4-2)根据计算得到的隐藏层数值，在输出层计算呼吸运动的预测结果的公式为：

其中z_k表示使用神经网络自回归NNAR模型计算得到的呼吸运动预测值,权重w_jk表示隐藏层与输出层之间的权重值。

Images