方園

安徽省第二人民醫(yī)院 （安徽合肥 230000）

時(shí)間序列為按照時(shí)間順序獲得的一系列觀測(cè)值的組合。隨著科學(xué)技術(shù)的不斷發(fā)展，對(duì)于時(shí)間序列的研究被應(yīng)用于金融、氣象、農(nóng)業(yè)生產(chǎn)等領(lǐng)域[1-3]。時(shí)間序列數(shù)據(jù)挖掘工作的主要任務(wù)是獲取數(shù)據(jù)中蘊(yùn)含的與時(shí)間相關(guān)的信息，掌握事物發(fā)展的規(guī)律。通過(guò)對(duì)時(shí)間歷史數(shù)據(jù)的挖掘，我們可以預(yù)測(cè)未來(lái)一段時(shí)間可能發(fā)生的數(shù)據(jù)，并以此分析事物未來(lái)的發(fā)展趨勢(shì)[4]。

在放射治療中，醫(yī)用直線加速器是不可或缺的設(shè)備，其輸出劑量直接決定了患者治療的預(yù)期效果，因此，保證輸出劑量的穩(wěn)定性非常重要[5]。由于醫(yī)用直線加速器是一類結(jié)構(gòu)復(fù)雜的高精度大型醫(yī)療設(shè)備，某些元器件的故障或老化均會(huì)導(dǎo)致加速器輸出劑量發(fā)生突發(fā)的或緩慢的變化，因此，對(duì)加速器輸出劑量的日常安全監(jiān)測(cè)成了放射治療質(zhì)量保證（quality assurance，QA）工作的重點(diǎn)[6-7]。很多學(xué)者利用統(tǒng)計(jì)學(xué)處理方法開展了大量關(guān)于醫(yī)用直線加速器劑量穩(wěn)定性的研究，且取得了較好的研究成果，為加速器劑量穩(wěn)定性的監(jiān)測(cè)工作提供了參考依據(jù)[8-10]。也有學(xué)者利用時(shí)間序列分析方法對(duì)醫(yī)用直線加速器質(zhì)量穩(wěn)定性進(jìn)行了分析和預(yù)測(cè)，并證實(shí)了此方法在加速器劑量預(yù)測(cè)工作中的應(yīng)用可行性[11]。鑒于此，本研究對(duì)醫(yī)用直線加速器質(zhì)控工作中的劑量監(jiān)測(cè)數(shù)據(jù)進(jìn)行了研究，并利用時(shí)間序列挖掘知識(shí)進(jìn)行了數(shù)據(jù)分析，以達(dá)到對(duì)加速器劑量偏移進(jìn)行預(yù)測(cè)的目的。

1 數(shù)據(jù)來(lái)源

本研究的數(shù)據(jù)來(lái)源于我院放射治療中心醫(yī)科達(dá)precise 直線加速器，選取3個(gè)調(diào)整周期內(nèi)共計(jì)66條加速器劑量實(shí)測(cè)數(shù)據(jù)作為實(shí)驗(yàn)樣本。

2 研究方法

首先，選取1個(gè)調(diào)整周期內(nèi)的前19條劑量監(jiān)測(cè)數(shù)據(jù)作為原始研究數(shù)據(jù)，對(duì)其進(jìn)行穩(wěn)定性分析；然后，利用時(shí)間序列研究中常用的差分自回歸移動(dòng)平均（auto regressive integrated moving average，ARIMA）模型對(duì)其進(jìn)行建模并預(yù)測(cè)后5次的測(cè)量結(jié)果；最后，將預(yù)測(cè)結(jié)果與另外兩個(gè)測(cè)量周期用同樣預(yù)測(cè)方法得出的預(yù)測(cè)結(jié)果進(jìn)行對(duì)比分析。

2.1 線性模型的建立

ARIMA 模型于1976年由Box 和Jenkins 等學(xué)者提出，之后被廣泛應(yīng)用于經(jīng)濟(jì)、金融和醫(yī)療等領(lǐng)域[12-14]。ARIMA（p，d，q）模型的通式為：

該模型分為前后兩部分，前部分?1yt-1+?2yt-2+...+?pyt-p是自回歸方程，后部分et-(θ1et-1+θ2et-2+...+θqet-q）是誤差移動(dòng)方程，其中，yt是當(dāng)前值，?1、?2...是自相關(guān)系數(shù)，θ1、θ2...是誤差系數(shù)，yt-1、yt-2...及et-1、et-2...是滯后項(xiàng)，et是誤差，p 是自回歸項(xiàng)數(shù)，q 為移動(dòng)平均項(xiàng)數(shù)，d 為時(shí)間序列成為平穩(wěn)時(shí)所作的差分次數(shù)。

2.2 應(yīng)用ARIMA 模型對(duì)數(shù)據(jù)進(jìn)行分析

存在一定發(fā)展趨勢(shì)的時(shí)間序列均是非平穩(wěn)的，而ARIMA模型必須建立在平穩(wěn)的基礎(chǔ)上才有意義[15]。因此，時(shí)間序列建模首先應(yīng)去除確定性因素，包括趨勢(shì)和季節(jié)性因素，然后檢測(cè)剩下數(shù)據(jù)的平穩(wěn)性，最終通過(guò)檢驗(yàn)的數(shù)據(jù)才可進(jìn)行ARIMA建模，建模的具體參數(shù)可以通過(guò)殘差自相關(guān)函數(shù)（autocorrelation function，ACF）和偏自相關(guān)函數(shù)（partial autocorrelation function，PACF）值來(lái)確定，R語(yǔ)言中有1個(gè)被命名為“auto.arima”的函數(shù)，可以用來(lái)直接定參，然后用此參數(shù)進(jìn)行模型的擬合和預(yù)測(cè)[16]。

auto.arima 函數(shù)是Hyndman-Khandakar 算法的1個(gè)變種，其結(jié)合了單位根檢驗(yàn)、最小化赤池信息準(zhǔn)則（Akaike information criterion，AIC）和極大似然估計(jì)（maximum likelihood estimate，MLE）等評(píng)價(jià)標(biāo)準(zhǔn)來(lái)獲得1個(gè)ARIMA 模型。

Hyndman-Khandakar 自動(dòng)ARIMA 建模算法步驟如下：

步驟1：通過(guò)重復(fù)地進(jìn)行單位根檢驗(yàn)（Kwiatkowski Phillips Schmidt Shin test，KPSS）測(cè)試來(lái)確定差分階數(shù)d（0≤d≤2）。

步驟2：對(duì)數(shù)據(jù)差分d 次之后，通過(guò)最小化AIC來(lái)選擇最優(yōu)的p，q；算法通過(guò)stepwise search 而不是遍歷所有可能的p，q 組合來(lái)尋找最優(yōu)的p，q 組合。

2.3 數(shù)據(jù)驗(yàn)證

我們通過(guò)對(duì)周期內(nèi)的數(shù)據(jù)進(jìn)行觀察，發(fā)現(xiàn)其呈現(xiàn)逐步上升的趨勢(shì)且不存在季節(jié)性因素，劑量實(shí)測(cè)數(shù)據(jù)序列見圖1。

圖1 劑量實(shí)測(cè)

經(jīng)過(guò)一階差分和R 語(yǔ)言中的auto.arima 函數(shù)，我們得到1個(gè)p，q 的最優(yōu)組合，即ARIMA（1，1，0）。劑量預(yù)測(cè)模型的擬合預(yù)測(cè)見圖2。

圖2 劑量預(yù)測(cè)模型的擬合預(yù)測(cè)

通過(guò)對(duì)預(yù)測(cè)模型的殘差A(yù)CF 和PACF 值的觀察，發(fā)現(xiàn)取值均在置信區(qū)間內(nèi)，因此殘差檢驗(yàn)可通過(guò)。預(yù)測(cè)模型的殘差A(yù)CF 和PACF 值見圖3；預(yù)測(cè)模型的參數(shù)估計(jì)見表1，由表可知，ARIMA（1，1，0）模型的參數(shù)在統(tǒng)計(jì)學(xué)上顯著。

表1 ARIMA 模型的參數(shù)

圖3 預(yù)測(cè)模型的殘差A(yù)CF 和PACF 值

3 結(jié)果與分析

經(jīng)過(guò)ARIMA 模型參數(shù)的選取、檢驗(yàn)和確認(rèn)，發(fā)現(xiàn)運(yùn)用此模型獲得的在1個(gè)預(yù)測(cè)周期內(nèi)的預(yù)測(cè)值能夠較好地反映醫(yī)用直線加速器劑量變化的范圍和趨勢(shì)，預(yù)測(cè)結(jié)果見圖4和表2。表2是周期內(nèi)第20～24次的劑量預(yù)測(cè)值與實(shí)測(cè)值的對(duì)比，通過(guò)對(duì)比可知，預(yù)測(cè)值對(duì)實(shí)測(cè)值有較好的體現(xiàn)，預(yù)測(cè)誤差為0.18%～0.80%。

表2 劑量預(yù)測(cè)數(shù)據(jù)與實(shí)際數(shù)據(jù)的對(duì)比

圖4 劑量預(yù)測(cè)值

在對(duì)數(shù)據(jù)結(jié)果進(jìn)行總結(jié)時(shí)，我們發(fā)現(xiàn)預(yù)測(cè)是基于1個(gè)固定周期內(nèi)的數(shù)據(jù)進(jìn)行的，為了明確ARIMA模型在其他周期內(nèi)的數(shù)據(jù)表現(xiàn)，本研究就此增加了兩個(gè)測(cè)量周期的數(shù)據(jù)并用上述方法進(jìn)行了預(yù)測(cè)和對(duì)比，發(fā)現(xiàn)ARIMA（1，1，0）模型對(duì)其他測(cè)量周期的數(shù)據(jù)預(yù)測(cè)也具有較好的表現(xiàn)。周期二、三的預(yù)測(cè)結(jié)果見圖5和表3。表3中周期二、三分別代表的是兩個(gè)測(cè)量周期第20～24次的實(shí)測(cè)值和預(yù)測(cè)值，由表可知，運(yùn)用此模型兩組預(yù)測(cè)值與實(shí)測(cè)值的對(duì)比表現(xiàn)良好，誤差為-1.28%～0.50%。

表3 周期二、三劑量預(yù)測(cè)數(shù)據(jù)與實(shí)際數(shù)據(jù)的對(duì)比

4 小結(jié)

本研究提出的利用時(shí)間序列進(jìn)行醫(yī)用直線加速器劑量預(yù)測(cè)方法，具有較好的預(yù)測(cè)精度，對(duì)加速器劑量參數(shù)的及時(shí)調(diào)整起到了參考和提示作用。但是，ARIMA 模型也存在一定的局限性，僅適用于正常檢測(cè)和調(diào)整內(nèi)的數(shù)據(jù)預(yù)測(cè)，當(dāng)醫(yī)用直線加速器出現(xiàn)維修和保養(yǎng)，且維修和保養(yǎng)的項(xiàng)目影響到劑量輸出時(shí)，此模型將不再適用，檢測(cè)周期應(yīng)以維修和保養(yǎng)后的時(shí)間開始計(jì)算。