計測自動制御学会東北支部 第 313 回研究集会 (2017.12.9)
資料番号 313-3

状態空間時系列モデルに基づく肺腫瘍位置の呼吸性移動時系列予測
Prediction of respiration-induced lung tumor motion by usingstate-space representation of time series component models○佐藤 雄介*，市地 慶*，本間 経康*○Yusuke Sato*, Kei Ichiji*, Noriyasu Homma**東北大学*Tohoku University

キーワード : 状態空間時系列モデル (State-space representation of time series component models), 放射線治療 (Radiation therapy), 自己回帰モデル (Autoregressive model), トレンド成分モデル (Trend component model), カルマンフィルタ (Kalman filter)

連絡先 : 〒980-8575 仙台市青葉区星陵 2-1 東北大学大学院医学系研究科医用画像工学分野佐藤 雄介,

Tel. : (022)717-7914, E-mail : y_sato@rii.med.tohoku.ac.jp

まうため，腫瘍のみに集中した放射線照射を実現

1. はじめにするためには，腫瘍位置変動のリアルタイムな計
現代のがん治療では，がんに対するさまざまな測と位置変動に応じた照射の制御といった呼吸性治療法が提案・確立されている一方で，がん患者移動対策が必要となる 2)．の数は依然として増えており，がん治療の重要性

有力な呼吸性移動対策の 1 つとして，呼吸に合は年々増している 1)．がんの主な治療方法には外科わせて照射範囲を動かす動体追尾照射が挙げられ的療法，化学療法，放射線治療の主に 3 つがあり，

る 2)． Fig. 1 に示すような臨床で広く利用されてこのうち放射線治療は低侵襲性・副作用の局所性いる放射線治療装置（リニアック）により動体追という点で他の治療法に比べて大きな利点を持つ

尾 照 射 を 行 う 場 合 ， Fig. 2 に 示 す Multi-leaf

ことからその利用拡大が進められている．

Collimator（MLC）と呼ばれる放射線遮蔽装置を放射線治療では治療効果を最大限に高めるため腫瘍の移動に合わせて動かして照射範囲を調整すに，十分な線量を患部に照射しながらも正常組織る方式が想定されている 2)．しかし，体内腫瘍位置への放射線被ばく障害を可能な限り低減することの計測と照射制御にかかる時間により数百ミリ秒が求められる．しかし肺腫瘍のような呼吸性移動の時間遅れが発生することが報告されている 3)．

を伴う腫瘍の場合は，照射中に腫瘍が移動してし-1-

Fig. 3 肺腫瘍の呼吸性移動時系列とその成分全体として複雑な振る舞いを見せる．
呼吸性移動時系列のより正確な予測の達成には，異なる複数の時系列成分をそれぞれ正確に予測す
Fig. 1

放射線治療装置（リニアック）ることが重要と考えられる．そこで本研究では，呼吸性移動時系列に含まれる主要な成分であるベースライン成分と振動的成分をそれぞれトレンド成分モデル（Trend component model）と自己回帰（Autoregressive model, AR）成分モデルによりそれぞれモデリングし，また観測された呼吸性移動時系列よりこれらの 2 成分を分離し予測する手Fig. 2

Multi-Leaf Collimator（MLC）法について検討する．2 つの時系列成分モデルの分離・予測と統合の実現にあたり，本稿では時系列動体追尾照射システムに内在する時間遅れを補成 分 モ デ ル の 状 態 空 間 表 現 （ State-space

償し，正確な照射を実現するためには数百ミリ秒

representation of time series component models）先の腫瘍位置の予測が必要である．このため，こを導入する 6,7)．複数の時系列成分モデルをひとつれまでに腫瘍位置の時系列予測手法として，線の状態空間モデルへと統合して表現することで，形・非線形の回帰手法をはじめとしてさまざまな呼吸性移動時系列からの各時系列成分の分離・予手法が提案されている 4)．動体追尾照射の正確性は

測は，観測値から状態の推定・予測の問題として腫瘍位置の呼吸性移動の予測の正確性に依存する扱うことでカルマンフィルタなどにより効率的にため，予測性能の追求は現在も課題のひとつとな予測計算を可能とする．本稿では，呼吸性移動時っている．

系列の実データを用いた予測実験により，AR+ト

肺腫瘍位置の呼吸性移動時系列の例を Fig. 3 にレンド成分モデル，AR モデル単独，そしてトレン示す．この例に見られるように，肺腫瘍の移動にド成分モデル単独の予測性能を比較し，複数成分は，呼息と吸息の繰り返しによる周期的・振動的を想定した状態空間時系列モデルの有効性を検証な成分が観察される．くわえて，ほかの体内臓器する．

からの圧迫などにより，振動的振る舞いの中心軸であるベースラインの移動も見て取ることができ

2. 状態空間時系列モデルを用いた呼吸る．このように肺腫瘍の呼吸性移動時系列は，異性移動の予測なる性質をもつ複数の時系列成分から構成され 5)，-2-

本研究では呼吸性移動時系列をベースライン成Fig. 4

提案モデルによる予測のイメージ分とベースライン周りの振動成分の 2 つの時系列

は状態遷移のダイナミクス，システム雑音の状態に分離して考え，各成分の時系列モデルの状態空への付加，観測雑音の観測時系列への付加を規定

間表現を単一の状態空間時系列モデルに統合する．

するそれぞれ 𝑙 × 𝑙, 𝑙 × 𝑞, 𝑙 × 1 の行列である．

各成分の分離と予測にはカルマンフィルタ(2) 式の状態空間時系列モデルは，時系列解析で（Kalman filter）を用いて状態の推定を行う．Fig.

用いられる線形・ガウス型時系列モデルの多くを4 に提案モデルによる予測のイメージを示す．

表現でき，統一的に取り扱うことを可能とする 6)．

2.1. 時系列モデルの状態空間表現

2.2. 呼吸性移動時系列の状態空間モデリング

状態空間時系列モデル，または時系列成分モデ

2.2.1. トレンド成分モデル

ルの状態空間表現（State-space representation ofFig. 3 に示したように呼吸性移動時系列には振time series component models）を用いることで，動的成分の中心軸のゆっくりとした時間変化であ観測時系列に含まれる各時系列成分の分離と予測るベースライン成分が含まれる．このベースライは状態推定の問題として定式化することが可能でン成分を表すモデルとして本研究ではトレンド成ある．本稿では，1 次元の観測時系列 {𝑦𝑛 }, 𝑛 =分モデルを用いる 7)．トレンド成分モデルは長期間

1,2, … に関する次式の線形・ガウス型状態空間時系にわたって増加や減少が続くような振る舞いの表列モデルを検討する：

現に有用なモデルの 1 つであり，時系列のおおよ

𝒙𝑛 = 𝑭𝒙𝑛−1 + 𝑮𝒗𝑛(1)

𝑦𝑛 = 𝑯𝒙𝑛 + 𝑤𝑛

(2)

その傾向を表すことが可能である．
トレンド成分の表現にあたり，トレンド成分のここで， 𝒙𝑛 は直接に観測できない 𝑙 次元の状態ベ

時系列 𝑇𝑛 の 𝑘 階差分を考える．

クトルであり，観測時系列に含まれる時系列成分

𝛥𝑘 𝑇𝑛 = 𝑣𝑛

(3)により構成される． 𝒗𝑛 はシステム雑音と呼ばれ，ここで， 𝛥𝑘 𝑇𝑛 は時系列 𝑇𝑛 の 𝑘 階差分， 𝑣𝑛 は平均

平均 0，分散共分散行列 𝑸 の 𝑞 次元正規分布に従0，分散 𝜏 2 の 1 次元正規分布に従うシステム雑音

う雑音である．一方，𝑤𝑛 は平均 0，分散共分散行であり，トレンド成分に生じる微小変化を表す．

列 𝑅 の 1 次元正規分布に従う雑音である．𝑭, 𝑮, 𝑯

呼吸性移動時系列に含まれるベースライン成分-3-の表現にあたっては，式 (3) で 𝑘 = 2 とした次の 2表される正規分布に従うシステム雑音である．
AR モデル単独に従う時系列 𝑦𝑛 は式 (11)より，階差分トレンドを用いる：(4)

𝑇𝑛 = 2𝑇𝑛−1 − 𝑇𝑛−2 + 𝑣𝑛

(12)𝑦𝑛 = 𝑝𝑛 + 𝑤𝑛

上式は，局所的には一定の速度がシステム雑音にで表される．ただし 𝑤𝑛 は平均 0，分散 𝜎 2 で表さより時間変化することに対応する．すなわち，れる 1 次元正規分布に従う観測雑音である．
(5)

AR 成分モデルの状態空間表現にあたっては，式トレンド成分モデル単独に従う時系列 𝑦𝑛 は，(5)(10)より以下の 𝑚 次の状態ベクトルを定義する．

𝑇𝑛 − 𝑇𝑛−1 = 𝑇𝑛−1 − 𝑇𝑛−2 + 𝑣𝑛 ．

𝒙𝑛 = (𝑝𝑛 , … , 𝑝𝑛−𝑚+1 )T

式のトレンド成分 𝑇𝑛 に観測雑音 𝑤𝑛 を用いて次の(13)このとき 𝑚 × 𝑚 行列 𝑭 および 𝑚 × 1 行列 𝑮 ，式で表される．
𝑦𝑛 = 𝑇𝑛 + 𝑤𝑛

(6)

1 × 𝑚 行列 𝑯 を

ただし観測雑音𝑤𝑛 は平均 0, 分散𝜎 2 で表される 1𝑎11𝑭=[⋮0次元正規分布に従う雑音である．
(4)式と(6)式のトレンド成分モデルを状態空間(7)

𝑝𝑛𝑎1𝑝𝑛−11[ ⋮ ]=[⋮𝑝𝑛−𝑚+10

列𝑯を

⋯ 0]

𝑎20
⋱0⋯⋯
1𝑎𝑚 𝑝𝑛−10 𝑝𝑛−2][ ⋮ ]⋮0 𝑝𝑛−𝑚(15)(8)10+ [ ] 𝑣𝑛
⋮0𝑯 = [ 1 0]と定義すれば，次のようなトレンド成分モデルの状態空間時系列モデルが得られる．
[(14)

時系列モデルが得られる．

このとき，2 × 2行列 𝑭 および2 × 1行列 𝑮 , 1 × 2行−11] , 𝑮 = [ ],001𝑎𝑚
100] , 𝑮 = [ ],⋮⋮00と定義すると，次のような AR モデルの状態空間クトルを下記のように定義する．

2𝑭=[1⋯⋯

𝑯 = [1 0時系列モデルで表現する．式 (5)より 2 次の状態ベ

𝒙𝑛 = (𝑇𝑛 , 𝑇𝑛−1 )T𝑎20⋱𝑇𝑛
2 −1 𝑇𝑛−11]=[][] + [ ] 𝑣𝑛
𝑇𝑛−1
1 0 𝑇𝑛−20(9)𝑇𝑦𝑛 = [1 0] [ 𝑛 ] + 𝑤𝑛
𝑇𝑛−1

(10)

𝑝𝑛𝑝𝑛−1𝑦𝑛 = [1 0 ⋯ 0] [ ⋮ ] + 𝑤𝑛
𝑝𝑛−𝑚+1(16)

2.2.3. AR+トレンド成分モデル
2.2.1 と 2.2.2 節で導入したトレンド成分モデルと AR モデルを組み合わせた次のようなモデル（以2.2.2. 自己回帰（AR）モデル
ベースライン周りの振動成分の表すモデルとし下，AR+トレンド成分モデル）を考える．て自己回帰（autoregressive, AR）モデルを用いる．

𝑦𝑛 = 𝑝𝑛 + 𝑇𝑛 + 𝑤𝑛

(17)

AR モデルは過去の自分の値を回帰的に利用する

ただし 𝑤𝑛 は平均 0，分散 𝜎 2 で表される正規分布モデルで， 𝑝𝑛 を定常時系列成分としたとき，に従う雑音である．AR+トレンド成分モデルは，

𝑚𝑝𝑛 = ∑ 𝑎𝑖 𝑝𝑛−𝑖 + 𝑣𝑛

ベースライン成分とベースライン周りの振動成分(11)を組み合わせることで呼吸性移動時系列を表現す

𝑖=1と表現される．ただし 𝑚 は自己回帰の次数， 𝑎𝑖 はるモデルである．

自己回帰係数である． 𝑣𝑛 は平均 0，分散 𝜂2 としてAR+トレンド成分モデルの状態空間時系列モデ-4-ルを考える．はじめにトレンド成分モデルと AR予測を行うには，2.2 節で導入した各時系列成分かモデルの状態ベクトルを連結することで，次のよら構成される状態 𝒙𝑛 を推定・予測する必要がある．

うな 𝑙 = (𝑘 + 𝑚) 次元状態ベクトルを定義する：

今回，トレンド成分モデル，AR モデル，AR+トレ

𝒙𝑛 = (𝑇𝑛 , … , 𝑇𝑛−𝑘+1 , 𝑝𝑛 , … , 𝑝𝑛−𝑚+1 )T(18)

ンドモデルはすべて線形・ガウス型の状態空間モ

トレンド成分 𝑇𝑛 と AR 成分 𝑝𝑛 はそれぞれ独立にデルである．そこで，本研究では Fig. 5 に示すよ式 (9)，(14) に従うと考えれば，(2 + 𝑚) × (2 + 𝑚)

うにカルマンフィルタ（Kalman filter）により状行列 𝑭 および 2 × (2 + m) 行列 𝑮 ，1 × (2 + 𝑚) 行態の推定・予測を行う．以下では，カルマンフィ

列 𝑯 を以下のように定義できる．

ルタによる状態の推定・予測と観測値の予測につ2 −10 10 0𝑭=⋮⋮⋮⋮[0 000𝑎11⋮0𝑯 = [1⋯⋯𝑎20⋱⋯⋯⋯⋯

0100𝑎𝑚0,𝑮 =00⋮⋮
[00]

10 10 ⋯001,0⋮0]

いて文献 7)8)に基づき説明する．
2.3.1. カルマンフィルタによる状態推定と予測(19)

カ ル マ ン フ ィ ル タ は 時 刻 𝑛 − 1 における 状 態
𝒙𝑛−1 の推定値の平均 𝒙𝑛−1|𝑛−1 と分散共分散行列0]

𝑽𝑛−1|𝑛−1 から，時刻 𝑛 の状態 𝒙𝑛 の予測値の平均以上をまとめると，次式のように AR モデルと𝒙𝑛|𝑛−1 と分散共分散行列 𝑽𝑛|𝑛−1 を得る“1 期先予トレンド成分モデルを統合したひとつの状態空間測”と，観測値 𝑦𝑛 を含む時刻 𝑛 までに得られた情時系列モデルが得られる．

報により状態予測値 𝒙𝑛|𝑛−1 , 𝑽𝑛|𝑛−1 を 𝒙𝑛|𝑛 , 𝑽𝑛|𝑛 へ

𝑇𝑛
𝑇𝑛−1
𝑝𝑛𝑝𝑛−1⋮2 −10 10 0=⋮⋮⋮⋮[𝑝𝑛−𝑚+1 ] [0 0𝑦𝑛 = [10 100𝑎11⋮00 ⋯⋯⋯
𝑎20⋱

0 𝑇𝑛−10 𝑇𝑛−2𝑎𝑚 𝑝𝑛−10 𝑝𝑛−2⋮⋮0 ] [𝑝𝑛−𝑚 ]と修正する“フィルタリング”と呼ばれる 2 つの100+0⋮
[0る．

⋯⋯⋯⋯
10]

𝑇𝑛
𝑇𝑛−1𝑝𝑛𝑝𝑛−1⋮

処理を交互に実施することで，観測値から状態を逐次推定するアルゴリズムである．具体的に 1 期先予測とフィルタリングは以下のように計算され001𝒗 (20)0 𝑛⋮0](1) 1 期先予測：
𝒙𝑛|𝑛−1 = 𝑭𝒙𝑛−1|𝑛−1(23)

𝑽𝑛|𝑛−1 = 𝑭𝑽𝑛−1|𝑛−1 𝑭T + 𝑮𝑸𝑮T

(24)(2) フィルタリング：+ 𝑤𝑛−1(25)

𝒙𝑛|𝑛 = 𝒙𝑛|𝑛−1 + 𝑲𝑛 (𝑦𝑛 − 𝑯𝒙𝑛|𝑛−1 )(26)

𝑽𝑛|𝑛 = (𝑰 − 𝑲𝑛 𝑯)𝑽𝑛|𝑛−1(27)

𝑲𝑛 = 𝑽𝑛|𝑛−1 𝑯T (𝑯𝑽𝑛|𝑛−1 𝑯T𝑛 + 𝑹)

(21)

[𝑝𝑛−𝑚+1 ]ただし，𝒗𝑛 は 2 次元正規分布に従うシステム雑音式 (23) にあるように，1 期先予測のアルゴリズであり，その平均 𝝁 と分散共分散行列 𝑸 は，トレムにおいて 𝒙𝑛 の予測値 𝒙𝑛|𝑛−1 は， 𝒙𝑛−1 のフィルンド成分モデルと AR モデルのシステム雑音の分タ値 𝒙𝑛−1|𝑛−1 に遷移行列 𝑭 をかけることで表され

散 𝜏 2 と 𝜂2 より，以下のように与えられる：る．また，その分散共分散行列 𝑽𝑛|𝑛−1 の時間遷移𝜏20𝝁 = [ ],𝑸 = [000]𝜂2は式 (24) で表されており，𝑭𝑽𝑛−1|𝑛−1 𝑭T は 𝑭 によ(22)

る変換の影響，𝑮𝑸𝑮T はシステム雑音 𝒗𝑛 の影響を表す．

2.3. カルマンフィルタによる時系列予測
状態空間時系列モデルを利用して，観測値 𝑦𝑛 の-5-Fig. 5

カルマンフィルタによる呼吸性移動時系列からの状態推定と予測フィルタリングのアルゴリズムでは，まず式2 期先予測と同様の手順を繰り返すことで，(25)により，カルマンゲインと呼ばれる 𝑲𝑛 が求め𝑖 = 1, … , 𝑗 期先予測に関する式が次のように得らられる．ついで，式 (26) により，状態 𝒙𝑛 のフィれる．

ルタ値 𝒙𝑛|𝑛 は，予測値 𝒙𝑛|𝑛−1 と観測値の 1 期先予(4) 多段階予測：

測誤差 𝑦𝑛 − 𝑯𝒙𝑛|𝑛−1 にカルマンゲインかけたもの

𝒙𝑛+𝑖|𝑛 = 𝑭𝒙𝑛+𝑖−1|𝑛(30)により更新される．同様に，分散共分散行列予測

𝑽𝑛+𝑖|𝑛 = 𝑭𝑽𝑛+𝑖−1|𝑛 𝑭 + 𝑮𝑸𝑮T

(31)

値 𝑽𝑛|𝑛−1 もカルマンゲインを用いてフィルタ値2.3.2. 観測値の多段予測
時系列成分（状態）の多段予測値 𝒙𝑛+𝑖|𝑛 が式 (30)𝑽𝑛|𝑛 へと更新される．
カルマンフィルタでは，式 (23) より状態の推定により求まると，観測値 𝑦𝑛+𝑖 の多段予測値 𝑦𝑛+𝑖|𝑛ともに 1 期先予測値を得る．この 1 期先予測処理

は式 (2) の観測方程式に基づき，次式により求めを繰り返す事で数期先の状態の予測も可能である．

られる：
𝑦𝑛+𝑖|𝑛 = 𝑯𝒙𝑛+𝑖|𝑛

いま，時刻 𝑛 までの観測時系列 {𝑦1 , … , 𝑦𝑛 } に基づいて 𝑗 期先の状態 𝒙𝑛+𝑗 (𝑗 > 1) を推定する場合を考(32)3. 予測性能の評価実験

える．このとき，カルマンフィルタにより 𝒙𝑛+1 の1 期先予測の平均 𝒙𝑛+1|𝑛 および分散共分散行列提案する AR+トレンド成分モデルによる呼吸性𝑽𝑛+1|𝑛 が求められる．ここで観測値 𝑦𝑛+1 は未知で移動時系列の予測性能を評価するため，実データを用いた予測実験を行った．

あるため，フィルタ処理は実行不可能である．そこでこれを省略し，予測値 𝒙𝑛+1|𝑛 および 𝑽𝑛+1|𝑛 を，3.1. 予測対象データ

時刻 𝑛 + 1 のフィルタ値 𝒙𝑛+1|𝑛+1 および 𝑽𝑛+1|𝑛+1

http://signals.rob.uni-luebeck.de においてオンとみなせば，次の 2 期先予測に関する式が得られライン公開されている腫瘍の呼吸性移動時系列デる．

ータベースのうちの 1 つを用いて予測性能の評価(3) 2 期先予測：を行った 9)．Fig. 6 に本研究で用いた呼吸性移動時

𝒙𝑛+2|𝑛 = 𝑭𝒙𝑛+1|𝑛(28)

𝑽𝑛+2|𝑛 = 𝑭𝑽𝑛+1|𝑛 𝑭 + 𝑮𝑸𝑮T

(29)

系列（’p10_f33_m1_proc2_pca’）を示す．
な お ， 本 デ ー タ ベ ー ス は 米 国 Georgetown-6-

University Hospital に お い て CyberKnife

Mean Square 誤差（RMSE）を用いた. RMSE はSynchrony システムにより計測された治療中の肺次式で表される．

がんの呼吸性移動時系列（肺癌患者 31 名 304 例分）𝑁

1
𝑅𝑀𝑆𝐸 = √ ∑(𝑦𝑖 − 𝑦̂𝑖 )2𝑁から構成されている．詳細な解析については文献(33)𝑖=15)を参照されたい．

ただし，𝑦𝑖 は観測値， 𝑦̂𝑖 はその予測値， 𝑁 は予測

3.2. 比較予測モデルとパラメータ設定

性能評価区間に含まれる観測データ数である．

実験に用いた予測手法とそのパラメータ設定を3.4. 実験結果

下記に示す．パラメータ推定には予測性能評価区
間の直前 60 s 分（60 s × 26 Hz = 1560 サンプル）各モデルによる 192 ms 先（𝑗 = 5 期先）予測をの観測値を用いた．

120 秒間（120 s × 26 Hz = 3120 サンプル）にわた(1)って実施した．

トレンド成分モデル単独

システム雑音の分散共分散行列 𝑸 と観測雑音のはじめに提案法である AR+トレンド成分モデル

分散共分散行列 𝑹 には，Neldar-Mead 法（滑降シにより抽出されたトレンド成分 𝑇𝑛|𝑛 と AR 成分ンプレックス法）を用いた数値最適化により，以𝑝𝑛|𝑛 の時系列を Fig. 7 と Fig. 8 にそれぞれ示す．

下の最尤推定値を設定した．

Fig. 7 に示すように，トレンド成分では全時間にわ

𝑄 = 1.515 × 10−9 , 𝑅 = 1.159 × 10−4(2)たってゆっくりとした変動が−330 mm から−325AR モデル単独mm の範囲に確認された．これは Fig. 6 の観測時

AR モ デ ル の 次 数 は 赤 池 情 報 量 規 準 に よ り

系列の中央付近の値によく対応している．このこ𝑚 = 61 とし，Burg 法による最小二乗推定値を使とから，抽出されたトレンド成分がベースライン用して AR 係数 𝑎1 , … , 𝑎𝑚 を求めた．また，システを表現できていることがわかる．また，Fig. 8 から

ム雑音と観測雑音の分散共分散行列 𝑸 と𝑹 はトレ

は Fig. 6 の観測時系列によく似ているが，0 mm をンド成分モデル単独と同様に Neldar-Mead 法によ中心軸とした振動成分が得られている．このようる以下の最尤推定値を用いた．に AR モデルにより抽出された時系列成分は，ベ𝑄 = 1.578 × 10−3 , 𝑅 = 3.501 × 10−4(3)ースライン周りの振動成分によく対応しているこ

AR+トレンドモデル（提案法）とがわかる．以上のように，AR+トレンド成分モ

AR モデル単独と同様に，AR モデル次数 𝑚 = 61

デルは呼吸性移動に含まれる時系列成分をうまくとし，AR 係数 𝑎1 , … , 𝑎𝑚 は Burg 法による最小二乗分離できることが確認された．

推定値を使用した 7)．システム雑音と観測雑音の分

次にそれぞれの予測モデルについての結果を示

散共分散行列 𝑸 と 𝑹 はトレンド成分モデル単独・す．Fig. 9 から 11 に示す予測結果は 120 秒間の各AR モデル単独と同様に Neldar-Mead 法による最予測モデルの予測結果のうち 50 秒から 80 秒の 30尤推定値を用いた．

秒間を示している．はじめにトレンド成分モデル

𝑸 = [2.881 × 100−9単独による 192 ms 先（5 期先）予測結果を Fig. 90],1.550 × 10−3に示す．Fig. 7 の AR+トレンド成分モデルにより

𝑅 = 1.071 × 10−6

分離されたベースライン成分と異なり，トレンド成分モデル単独の予測値（赤色実線）は振動的な3.3. 性能評価指標振る舞いも表現できており，観測時系列とよく一性 能 評 価 指 標 に は 予 測 さ れ た 観 測 値 の Root-7-Fig. 6

予測対象の肺腫瘍位置の呼吸性移動時系列（観測時系列）

Fig. 7

AR+トレンドモデルにより抽出されたベースライン成分

Fig. 8AR+トレンドモデルにより抽出された振動的成分

致した変動を示した．一方，振動的な振る舞いのても観測時系列とよく一致した変動を示した．一頂点付近では予測値が実測値と大きく乖離する現方，トレンド成分モデルでも見られた振動的な振象が確認され，その周辺で予測誤差が生じているる舞いの頂点付近で予測値が実測値と乖離する現ことが確認できた．これはトレンド成分モデルが象がわずかに確認され，その周辺で予測誤差が生式(4)で表されるように局所的な速度を一定とみなじていることが確認できた．これはベースラインしたモデルであることから頂点付近での速度変化成分を表現できていないことが原因と考えられる．に対応しきれていないことが原因と考えられる．

AR モデル単独の予測誤差は RMSE=0.219 mm で

ト レ ン ド 成 分 モ デ ル 単 独 予 測 誤 差 は RMSE ＝あった．

0.272 mm となった．

最後に，提案法である AR+トレンド成分モデル

次に AR モデル単独による予測結果を Fig. 10 にによる予測結果を Fig. 11 に示す．図から見て取れ示す．AR モデル単独の予測値（緑色実線）についるように，AR＋トレンド成分モデルによる予測値-8-Fig. 9

トレンド成分モデル単独による予測結果

Fig. 10AR モデル単独による予測結果

Fig. 11

AR+トレンド成分モデルによる予測結果（赤色実線）は，AR 単独モデルとよく似た変動をAR 成分を分離・予測することで時系列モデル単独示し，観測値ともよく一致した変動を示した．まを用いる場合よりも予測性能が向上することが確た，Fig. 10 の AR モデル単独の予測値に見られた認された．

ような振幅の頂点付近での誤差が減少しているこ

3.5. 考察とが確認できた． AR モデルとトレンド成分モデ

本研究で提案した提案法の AR+トレンド成分モルの併用により，AR モデル単独では表現しきれなデルは AR モデル単独，トレンドモデル単独よりかったベースライン成分の変動が表現され，誤差も優れた予測性能を示した．動体追尾照射におい減少の効果があったと考えられる．AR+トレンドて呼吸性移動の予測性能は，放射線照射の正確成分モデルの予測誤差は RMSE=0.216 mm と今回度・精度に直結するため，今後も高い予測性能が評価したモデルの中で最小であった．以上のよう求められると考えられる．提案法のように呼吸性に，状態空間時系列モデルによりトレンド成分と移動時系列を主要な時系列成分に分離して予測す-9-

るアプローチは，今後の更なる予測性能追求に有

参考文献

効であると考えられる．

1)

本研究では，呼吸性移動時系列をガウス分布に

国立がん研究センターがん情報サービス :がん登録・統計, 年次推移，従う線形・時不変の状態空間モデルで表現し，カ

http://ganjoho.jp/reg_stat/statistics/stat/an

ルマンフィルタによる推定・予測を行った．一方，

nual.html

実際の肺腫瘍の呼吸性移動のダイナミクスには非

2)

P. J. Keall et al: The management of

線形性・時変性が含まれ，またその分布はガウス

respiratory motion in radiation oncology型とは限らないと考えられる．したがって，非線

report of AAPM Task Group 76, Med. Phys.,

形・時変・非ガウス型のモデルにより時系列成分

33-10, pp. 3874-900 (2006)

を表現し，これを拡張カルマンフィルタや粒子フ

3)

P. R. Poulsen et al.: Detailed analysis of

ィルタといった手法により推定・予測することでlatencies in image-based dynamic MLC

更なる予測性能向上が期待される．

tracking, Med. Phys., 37-9, pp. 4998-5005(2010)

4. おわりに4)

本研究では，放射線治療の肺腫瘍の動体追尾照P. S. Verma et al: Survey: Real-Time Tumor
MotionPrediction

forImage-Guided

射において必須技術である呼吸性移動の時系列予

Radiation Treatment, Comput. Sci. Eng.,

測の性能向上のために, 状態空間時系列モデルを

13-5, pp.24-35 (2011)

用いて呼吸性移動時系列に含まれる異なる性質の5)

Y. Suh et al: An analysis of thoracic and複数の時系列成分への分離・予測を試みた．単独

abdominal tumour motion for stereotactic

成分モデルとの予測性能比較実験の結果より，

body radiotherapy patients, Phys. Med.

AR+トレンド成分モデルの誤差が一番小さく，単独成分モデルによる予測と比べて予測性能が向上

Biol., 53-13, pp. 3623-40 (2008)6)北川源四郎, 佐藤整尚 : 一般状態空間モデルしたことが確認された．このことは，状態空間時に よ る 分 散 変 動 時 系 列 の 解 析 , IMES

系列モデルにより，呼吸性移動時系列を複数の異Discussion Paper Series, Discussion Paperなる性質を持つ時系列成分への分離・予測する手

No. 98-J-22 (1998)

法の有効性を示唆する．

7)

公開データベースの全データを用いて予測を行い，本手法の信頼性を確認すること，ならびに同(2005)
8)

一データベースによる検証を行っている先行研究の予測手法との性能比較が今後の課題である．

北川源四郎 : 時系列解析入門, 岩波書店

足立修一, 丸田一郎 : カルマンフィルタの基礎, 東京電機大学出版局 (2012)

9)

F. Ernst: Compensating for Quasi-periodic
MotioninRobotic

Radiosurgery,

Springer-Verlag New York (2012)- 10 -