計測⾃動制御学会東北⽀部 第 341 回研究集会 （2023.3.23）
資料番号 341-6

状態空間時系列モデルに基づく粒⼦フィルタを⽤いた
腫瘍位置移動の適応的予測の試み
Adaptive prediction of intrafractional tumor motion based on state spacemodel of time-series components using particle filter◯篠原 唯*, 市地 慶*, ⽯井 万結*, 淡路 樹*, 本間 経康*◯Yui Shinohara*, Kei Ichiji*, Mayu Ishii, Itsuki Awaji, Noriyasu Homma**東北⼤学*Tohoku University
キーワード：放射線治療（radiation therapy）, 肺がん（lung cancer）, ⾃⼰組織化状態空間モデル（self-organized state-space model）, 粒⼦フィルタ（particle filter）

連絡先：〒980-8575 仙台市⻘葉区星陵 2-1 東北⼤学⼤学院医学系研究科 医⽤画像⼯学分野篠原 唯, Tel.： 022-717-8190, E-mail：yui.shinohara.p8@dc.tohoku.ac.jpの腫瘍位置変動は照射精度に影響するため，対策が必要である 4)．治療中の腫瘍位置変動対策として
は，Fig.2 のように治療装置に取り付けられている
透視装置を⽤いてリアルタイムに腫瘍の位置を計
測 し ， マ ル チ リ ー フ コ リ メ ー タ （ multi-leaf
collimator; MLC）を腫瘍に合わせて変形させながら照射範囲を制御する動態追尾照射などが開発・利⽤されている 5)．しかし，腫瘍位置の計測・照射範囲の制御システムには数⼗から数百ミリ秒程度の時間遅れが発⽣することが報告されている 5)．
この時間遅れの補償には，腫瘍位置の計測に基づく将来の腫瘍位置変動の時系列予測（以下，腫瘍
位置変動予測）を⾏い，予測位置に基づき照射の先回り制御する⽅法が有効である 6)．腫瘍位置変動予測は，直近の腫瘍位置変動とその先の腫瘍位置の変化の関係を数理モデルで表現することで実現され，これまでにさまざまな腫瘍位置変動の予測モ
デルが開発・提案されている 7-12)．例えば，Support
vector machine や Kernel density estimation といった ⼀ 般 的 な 機 械 学 習 ⼿ 法 7) の ほ か 近 年 はTransformer モデルなど深層学習モデルを応⽤した⼿法 8)9)が登場している．特に深層学習に基づく予測モデルは⾼い予測性能を達成可能と期待されている 13)．これら学習に基づくモデルは，通常，治療開始以前に取得されたデータを⽤いて腫瘍位置変動の履歴を⼊⼒とし，その先の腫瘍位置変動の予測誤差を最⼩化するように訓練が⾏われ，予測のための数理表現は必ずしも明⽩ではない．他⽅，腫瘍位置変動は Fig. 1 で⽰したように異なる

1. はじめに⾼齢化社会の進展によりがんの罹患数は増加し続けており 1)，がんの治療法の治療効果向上と副作⽤低減のための研究・開発が継続的に⾏われている．がんの治療法のうち，放射線治療は，低侵襲かつ副作⽤の局所性から，⾼齢者などの体⼒の低下した患者でも利⽤しやすい治療として⽤いられている 2)．放射線治療の治療効果を⾼めるためにはがん細胞を死滅させるのに⼗分な線量を治療対象へ照射すると同時に，副作⽤を極⼒低減するために治療対象周辺の正常組織への被曝を抑制することが求められる．
このような治療の実現には，3 次元的な照射を⾏う体幹部定位放射線治療や照射野内の放射線強度を変化させる強度変調放射線治療などに代表されるように腫瘍周辺のみに限局した⾼精度な照射が必要となる 2)．しかし，肺がんや肝がんなどでは腫瘍の位置が変動する腫瘍位置変動が⽣じる 3)ため，腫瘍への線量不⾜や正常組織に余分に照射してしまう可能性がある．Fig.1 に肺がんの腫瘍位置変動の例を⽰す．
この例では主要な成分として 5〜10 秒
程度の呼吸による繰り返しの動きがみられる．⼀⽅，1 秒ほどの周期で⽣じている拍動や，呼吸性の周期的変動の中⼼基軸のずれ（ベースライン変動）
として他臓器からの圧迫による呼吸や拍動とは異なるもっとゆっくりとした動きも含まれており，
腫瘍位置変動には性質の異なる複数の動きが含まれていることがわかる．このような放射線治療中

1

繰り返しの動き
拍動

ベースラインの変化

Fig. 1 肺がんの腫瘍位置変動の例
ビームの制御腫瘍位置の計測

Fig. 2 腫瘍位置変動に応じた照射システムの例（動体追尾照射）
性質の複数の動きの成分により構成される．学習に基づくモデルが良好な予測性能を達成するためには，これら複数の動きの成分の膨⼤な組み合わせからなる腫瘍位置変動の複雑なパターンを網羅的に表現することが求められる．しかし，限られた臨床データから適切な量・質のデータを準備し，それらをもとに許容可能な計算コスト・時間で実⽤上⼗分な予測性能を発揮可能な機械学習系モデルを構築・訓練する⽅法は確⽴されておらず，いまだ発展途上にある．
複雑な腫瘍位置変動パターンを先験知識に基づいてその成分・要素ごとに明⽰的に分けて取り扱うことができ，かつ訓練データへ強く依存せずに構築・利⽤可能な⼿法として，従来から統計的時系
列モデルに基づく腫瘍位置変動予測法が提案されている 10-12)．たとえば，先⾏研究 11)では，状態空間時系列モデルにより腫瘍位置変動時系列に含まれるベースライン成分と振動的成分をそれぞれトレ ン ド 成 分 モ デ ル と ⾃ ⼰ 回 帰 （ Autoregressivemodel, AR）成分モデルで表現し，成分ごとにカルマンフィルタにより分離・予測し，それらを統合することで腫瘍位置変動を予測する⼿法が提案されている．この⼿法は性能良好な機械学習系モデル
14)15)に匹敵する予測性能を達成可能と報告されており，状態空間時系列モデルに基づき複数成分を推定・予測するアプローチの有効性が⽰唆されている．⼀⽅，この⼿法では各時系列成分モデルのパラメータは固定されていることから，呼吸周期や振幅のゆらぎをはじめとする各成分の性質の時間変化によって予測誤差が増⼤するおそれがある．
腫瘍位置変動時系列に含まれる時系列成分の性質変化へ対応するためには，状態空間時系列モデルのパラメータを適応的に時間変化させることが有効と考えられる．
そこで本研究では，より正確な腫瘍位置変動予測の実現を⽬的とし，統計的時系列モデルに基づ
き，時変の時系列成分を適応的に分離・推定し，腫瘍位置変動を予測する⼿法について提案する．具体的に各時系列成分モデルのパラメータの時間変化を表現・推定するため，提案法では時変パラメー
タも状態として取り込む⾃⼰組織化状態空間モデ
ルを導⼊する．これと粒⼦フィルタと組み合わせることで，各時系列成分とそのパラメータを同時推定し，適応的な時系列成分分離とそれに基づく
腫瘍位置変動予測を試みる．性質変化を伴う呼吸波形を模した⼈⼯時系列データと腫瘍位置変動時
系列の実データを⽤いた予測実験により，提案法の有効性を検証する．

2. ⾃⼰組織化状態空間モデルを⽤いた腫瘍位置変動の予測
本研究では，先⾏研究に倣って腫瘍位置変動時
系列をベースライン成分とベースライン付近の振動的な成分の 2 つの時系列成分から構成されるものと考え，各成分の時系列モデルの状態空間表現
を 1 つの状態空間時系列モデルに統合する．
また，時間変化への対応のために，各時系列成分モデルのパラメータについても状態を定義し，時系列成分モデルの状態と組み合わせた拡⼤状態ベクトルを定義した⾃⼰組織化状態空間モデルによって予測する．

2

次元の観測時系列 {𝑦! }, 𝑛 = 1,2, … に関する次式の線形・ガウス型状態空間時系列モデルを検討する．

2.1. 腫瘍位置変動の時系列表現：AR＋トレンド成分モデル
本研究では，先⾏研究 11)と同様に腫瘍位置変動
時系列をトレンド成分モデルと⾃⼰回帰（Autoregressive; AR）モデルを組み合わせた次のようなモデル（以下，AR+トレンド成分モデル）で表現する．
𝑦! = 𝑇! + 𝑝! + 𝑤!
(1)
ここで 𝑦! は観測値の時系列であり，𝑝! は AR 成
分，𝑇! はトレンド成分である．また，𝑤! は平均 0，
分散 𝜎 " で表される正規分布に従う観測ノイズである．
AR+トレンド成分モデルを構成する各時系列成分モデルについて以下で詳しく説明する．

𝑦! = 𝑯𝒙! + 𝑤!

(5)

うノイズである．⼀⽅，𝑤! は平均 0，分散共分散
⾏列 𝑅 = 𝜎 " の 1 次元正規分布に従うノイズである．
𝑭, 𝑮, 𝑯 はそれぞれ状態遷移のダイナミクス，シス
テムノイズの状態への付加，観測ノイズの観測時
系列への付加を規定する 𝑙 × 𝑙, 𝑙 × 𝑘, 𝑙 × 1 の係数⾏列である．
2.2.1. AR 成分モデルの状態空間表現

呼吸性変動など腫瘍位置変動時系列に含まれる
周期的・振動的な振る舞いを表すモデルとして⾃⼰回帰（steady autoregressive, AR）モデルを⽤いる．AR モデルは過去の時系列値を回帰的に利⽤するモデルであり，

式(2)の 𝑞 次 AR 成分モデルの状態空間表現を考える．AR 成分時系列値より状態ベクトル𝒙(,! を𝒙(,! = [𝑝! , … , 𝑝!$%+' ],(6)とし，係数⾏列 𝑭( , 𝑮( , 𝑯( をそれぞれ次のように定義する：
𝑎' 𝑎" ⋯ 𝑎%1
1 0 ⋯ 00𝑭( = D
H , 𝑮( = D H,
⋮
⋱
⋮
⋮
(7)
0
1 0
0
𝑯( = [1 0 ⋯ 0]

%

(2)

#&'

と表現される．ただし，𝑝! は AR 成分の時系列値，
𝑞 は⾃⼰回帰モデルの次数，𝑎# , 𝑖 = 1,2, … , 𝑞 は⾃⼰回帰係数である． 𝑣(,! は平均 0，分散 𝜂" の正規分布に従うシステムノイズである．

上記の状態ベクトル・係数⾏列を⽤いると次のよ
うな AR 成分モデルの状態空間時系列モデルが得られる．

2.1.2. トレンド成分モデル
腫瘍位置変動時系列には呼吸や拍動のほかにもそれらの周期よりもゆっくりとした時間変化であるベースライン成分が含まれる．このベースライ
ン成分は周辺臓器の蠕動，筋緊張・弛緩や重⼒からの影響を受けていると考えられる．ベースライン
成分を表すモデルとしてトレンド成分モデルを⽤いる．トレンド成分モデルは⽐較的⻑期にわたって緩やかに増加・減少が続くような振る舞いの表
現に有⽤なモデルの 1 つであり，時系列のおおよ
その傾向を表すことが可能である 16)．
トレンド成分は 1 時刻につきわずかに変化すると仮定すると，トレンド成分の時系列値 𝑇! は次のような１次トレンドモデルで表現できる．
𝑇! = 𝑇!$' + 𝑣*,!

(4)

ここで， 𝒙! は直接に観測できない 𝑙 次元の状態ベ
クトルであり，観測時系列に含まれる時系列成分
により構成される． 𝒗! はシステムノイズと呼ばれ，
平均 0，分散共分散⾏列 𝑸 の 𝑘 次元正規分布に従2.1.1. AR 成分モデル

𝑝! = ) 𝑎# 𝑝!$# + 𝑣(,!

𝒙! = 𝑭𝒙!$' + 𝑮𝒗!

𝒙(,! = 𝑭( 𝒙(,!$' + 𝑮( 𝒗(,!
𝑎' 𝑎" ⋯ 𝑎% 𝑝!$'1
1 0 ⋯ 0 𝑝!$"0
=D
H D ⋮ H + D H 𝒗(,! ,⋮⋱⋮⋮0 0 1 0 𝑝!$%0𝒗(,! ~𝒩(0, 𝜂" )(8)𝑦! = 𝑯( 𝒙(,! + 𝑤(,!
𝑝!
𝑝!$'
= [1 0 ⋯ 0] D ⋮ H + 𝑤!
𝑝!$%+'(9)2.2.2. トレンド成分モデルの状態空間表現
前節で⽰した１次トレンド成分モデルの状態空間表現について考える．式(3)より，トレンド成分
モデルの状態 𝒙*,! および係数⾏列 𝑭* , 𝑮* , 𝑯* をそ
れぞれ

(3)

ここで 𝑣*,! は平均 0，分散 𝜏 " の 1 次元正規分布に従うシステムノイズである．
2.2. AR＋トレンド成分モデルの状態空間表現
前節の２成分から構成される時系列モデルを状態空間モデルとして表現することで，観測時系列に含まれる各時系列成分の分離・予測を状態推定の問題として考えることができる 16)．ここでは 1

𝒙*,! = [𝑇! ](10)

𝑭* = 𝑮* = 𝑯* = [1](11)と定義する．すると，次のようなトレンド成分モデ
ルの状態空間時系列表現が得られる．

3

𝒙*,! = 𝑭* 𝒙*,!$' + 𝑮* 𝒗*,!
= 𝑇!$' + 𝒗*,! ,𝒗*,! ~𝒩(0, 𝜏 " )(12)𝑦! = 𝑯* 𝒙*,! + 𝑤!
= 𝑇! + 𝑤!

(13)を時間変化させることを考える．前節の式(16)(17)のモデルには，AR 係数 𝑎# ，AR 成分のシステムノ
イズ分散 𝜂" ，トレンド成分のシステムノイズ分散
𝜏 " が時系列成分モデルのパラメータとして存在す
る1．これらはそれぞれ周期的成分の周期性，振幅，
そしてベースライン成分の緩急の変化に関連する．
以下では，2 次 AR＋トレンド成分モデルを例に
これらパラメータ時間変化の導⼊を考える．2 次
AR＋トレンド成分モデルの係数⾏列 𝑭, 𝑮, 𝑯 は2.2.3. AR＋トレンド成分モデル
AR+トレンド成分モデルの状態空間時系列モデ
ルを考える．はじめにトレンド成分モデルと AR モ
デルの状態ベクトルを連結することで，次のような 𝑙 = (1 + 𝑞) 次元状態ベクトル𝒙! を定義する．
,

𝒙! = M𝑇! , … , 𝑇!$-+' , 𝑝! , … , 𝑝!$%+' N

1 0 0
1 0
𝑭 = ]0 𝑎' 𝑎" ^ , 𝑮 = ]0 1^,
(19)
0 1 0
0 0
𝑯 = [1 1 0]のようにそれぞれ定義でき，2 次 AR＋トレンド成分モデルは次のように状態空間モデルとして表現
できる：

(14)

さらに係数⾏列 𝑭, 𝑮, 𝑯 は式(7)，式(11)で定義した
各時系列成分の係数⾏列を組み合わることでそれぞれ以下のように定義できる．
𝑭
𝑭=O *

𝑮=O

𝑮*

𝐻 = [𝑯 *1 0 ⋯ ⋯ 0
⎡0 𝑎 𝑎 ⋯ 𝑎 ⎤
'
"
%
⎢
⎥
P = ⎢⋮ 1 0 ⋯ 0 ⎥ ,
𝑭(
⎢⋮ ⋮
⋱
⋮⎥
⎣0 0
1 0⎦
1 0
⎡0 1⎤
⎢
⎥
P = ⎢0 0⎥ ,
𝑮(
⎢⋮ ⋮⎥
⎣0 0⎦
𝑯( ] = [1 1 0 ⋯ 0]𝑇!
1 0
] 𝑝! ^ = ]0 𝑎'
𝑝!$'
0 1

(15)

0
1⎤
⎥
0⎥ 𝒗!
⋮⎥
0⎦

𝑇!
⎡ 𝑝
⎤
!
⎢
⎥
𝑝
[
]
𝑦! = 1 1 0 ⋯ 0 ⎢ !$' ⎥ + 𝑤!
⎢ ⋮ ⎥⎣𝑝!$.+' ⎦

(21)

𝑇!
1
0
] 𝑝! ^ = ]0 𝑎',!
𝑝!$'
0
1

(16)

𝑇!$'
0
1 0
𝑎",! ^ ]𝑝!$' ^ + ]0 1^ 𝒗!
0 0
0 𝑝!$"

𝜏"𝑸=[ !
0𝑦! = 𝑯𝒙! + 𝑤!

𝑇!
⎡ 𝑝 ⎤
⎢ ! ⎥
= [1 1 0 ⋯ 0] ⎢ 𝑝!$' ⎥ + 𝑤!
⎢ ⋮ ⎥⎣𝑝!$%+'⎦(20)

さらに，2 次 AR＋トレンド成分モデルのパラメータ (𝑎' , 𝑎" , 𝜏 " , 𝜂" ) を時刻 𝑛 により変化するもの
として，_𝑎',! , 𝑎",! , 𝜏!" , 𝜂!" ` と書き換えると以下の時変状態空間モデルを得る．

以上をまとめると，次式のように AR モデルとトレ
ンド成分モデルを統合したひとつの状態空間時系
列モデルが得られる．
𝒙! = 𝑭𝒙!$' + 𝑮𝒗!
1 0 ⋯ ⋯ 0 𝑇!$'1
⎡0 𝑎 𝑎 ⋯ 𝑎 ⎤ ⎡𝑝 ⎤ ⎡0
!$'
'
"
%
⎥ ⎢
⎢
⎥⎢
= ⎢ ⋮ 1 0 ⋯ 0 ⎥ ⎢𝑝!$" ⎥ + ⎢0
⋮
⎢⋮ ⋮
⋱
⋮ ⎥⎢
⎥ ⎢⋮
⎣0 0
1 0 ⎦ ⎣𝑝!$% ⎦ ⎣0

0 𝑇!$'
1 0
𝑎" ^ ]𝑝!$' ^ + ]0 1^ 𝒗!
0 𝑝!$"
0 0

0
\
𝜂!"

(22)

(23)

次に，パラメータ _a𝑎#,! b, 𝜏!" , 𝜂!" ` の具体的な時
間変化を時系列モデルとして表現することを考える．いま，それぞれのパラメータの単位時間当たりの変化量は微⼩と考えると，パラメータ _a𝑎#,! b, 𝜏!" ,𝜂!" ` は１次トレンドモデルによって以下の式のように表現することができる．

(17)

ただし，𝒗! は 2 次元正規分布に従うシステムノイ
ズであり，その平均 𝝁 と分散共分散⾏列 𝑸 は，トレンド成分モデルと AR モデルのシステムノイズの分散 𝜏 " と 𝜂" より，以下のように与えられる．

腫瘍位置変動時系列に含まれる成分の性質変化への適応のために，状態空間モデルのパラメータ𝑎',! = 𝑎',!$' + 𝑣/,!
(24)
(25)
𝑎",! = 𝑎",!$' + 𝑣/,!
"
(26)
𝜏!" = 𝜏!$'+ 𝑣0!,!
"
𝜂!" = 𝜂!$'+ 𝑣1!,!
(27)
ここで， 𝑣/,! , 𝑣0!,! , 𝑣1!,! はそれぞれ平均 0，分散
σ"/ , σ"0! , σ"1! の正規分布に従うシステムノイズである．

1

分散𝜎 ! も存在するが本研究ではこれを時不変とする．

"

𝜏0𝝁 = Y Z,𝑸 = [
0
0

0
\
𝜂"

(18)

2.2.4. 適応的 AR＋トレンド成分モデル

状態空間モデル全体のパラメータとして観測ノイズの4パラメータ _𝑎',! , 𝑎",! , 𝜏!" , 𝜂!" `を状態ベクトル

さらに拡⼤状態ベクトルに対応するように係数⾏列および分散共分散⾏列をそれぞれ以下のように定義する．

*
𝜽! = M𝑎',! , 𝑎",! , 𝜏!" , 𝜂!" N とし，式(24)から式(27)をまとめて表現すると次の状態遷移⽅程式を得る：
𝜽! = 𝑭2,! 𝜽!$' + 𝑮2 𝒗2,!
1 0 0 01 0 0 00 1 0 00 1 0 0
(28)
=D
+D
H𝜽
H𝒗
0 0 1 0 !$'0 0 1 0 2,!
0 0 0 10 0 0 1
ただし， 𝒗2,! は 4 次元正規分布に従うシステムノイズであり，その平均 𝝁2 と分散共分散⾏列 𝚺 は，
パラメータの時系列モデルのシステムノイズ分散σ"/ , σ"0! , σ"1! より，以下のように与えられる．
σ"
0
⎡ /
⎤
σ"/
⎢
⎥
0
𝝁2 = D H , 𝚺 = ⎢
(29)
⎥
σ"0!
0
⎢
⎥
"
0
σ1! ⎦
⎣

1
⎡
𝑎',! 𝑎",!
⎢1
0⎢
𝑭! = ⎢
1 0 0
⎢
0 1 0
⎢
0 0 1
⎣
0 0 0
1
⎡
⎤
1
⎢
⎥
0
⎢
⎥
𝑮=⎢
1
⎥,
1
⎢
⎥
⎢
⎥
1
⎣
1⎦
𝑯 = [1 1 0 0 0 0 0]式(16)の AR＋トレンド成分モデルの状態に，式(28)のパラメータ状態を組み合わせることで，次のように拡⼤状態ベクトルを定義することができ
る.
𝑇!
⎡ 𝑝 ⎤
⎢ ! ⎥
⎢𝑝!$' ⎥
𝒙! = ⎢ 𝑎',! ⎥⎢ 𝑎",! ⎥⎢ " ⎥⎢ 𝜏! ⎥⎣ 𝜂!" ⎦𝑇!
⎡ 𝑝 ⎤ ⎡1
⎢ ! ⎥ ⎢
⎢𝑝!$' ⎥ ⎢⎢ 𝑎',! ⎥ = ⎢⎢ 𝑎",! ⎥⎢ " ⎥ ⎢⎢ 𝜏! ⎥ ⎢⎣ 𝜂!" ⎦ ⎣

𝑎',!
1

⎤
⎥
⎥
0⎥,
0⎥
0⎥
1⎦

(31)

式(30)(31)をまとめると，次式のように AR モデ
ルとトレンド成分モデル，それぞれのパラメータを統合したひとつの状態空間時系列モデルが得ら
れる．

(30)

𝑎",!
0

1
0
0
0

𝑇!$'
1
⎤ ⎡ 𝑝!$' ⎤ ⎡⎥ ⎢ 𝑝!$" ⎥ ⎢⎥ ⎢⎥ ⎢𝑎
0 ⎥ ⎢ ',!$' ⎥ + ⎢⎢𝑎⎥
0 ⎥ ⎢ ",!$' ⎥ ⎢"𝜏
0 ⎥ ⎢ !$' ⎥ ⎢"1 ⎦ ⎣ 𝜂!$'
⎦ ⎣

0 0
1 0
0 1
0 0
𝑇!
⎡ 𝑝 ⎤
⎢ ! ⎥
⎢𝑝!$' ⎥
𝑦! = [1 1 0 0 0 0 0] ⎢ 𝑎',! ⎥ + 𝑤!
⎢ 𝑎",! ⎥⎢ " ⎥⎢ 𝜏! ⎥⎣ 𝜂!" ⎦

1
0

1

1

1

⎤
⎥
⎥
⎥ 𝒗! ,
⎥
⎥
1⎦

(32)

𝜏"
⎡
⎤
𝜂"
⎢
⎥
"
σ/3
⎢
⎥
(33)
𝑸! = ⎢
σ"/3
⎥
⎢
⎥
σ"0!
⎢
⎥
"
σ1! ⎦
⎣
式(33)のように元の状態空間モデルのパラメー
タを拡⼤状態ベクトルとして推定対象に含める時
変状態空間モデルを⾃⼰組織化状態空間モデルと呼ぶ 17)．

ただし，𝒗! は平均 𝝁 = [0,0,0,0,0,0]* と時変の分散
共分散⾏列 𝑸! の 6 次元正規分布に従うシステム
ノイズである．ここで 𝑸! は，トレンド成分モデル，
AR 成分モデル，パラメータの時系列モデルのシス
テムノイズ分散 𝜏 " , 𝜂" , σ"/ , σ"0! , σ"1! より，以下のよ
うに与えられる：5

初期分布の⽣成

(a)

𝑛=𝑛+1
システムノイズの⽣成

粒⼦ 𝑗 = 1, … , 𝑚 について乱数を⽣成し，時刻 𝑛
で状態 𝒙! に⽣じうる摂動の 𝑚 個の独⽴な実現値とみなせるシステムノイズを得る．

1
𝑚
𝒙0|0
, … , 𝒙0|0

1

𝑚

𝑗

𝒗𝑛 , … , 𝒗𝑛

𝑗

, 𝒗𝑛 ~𝑝 𝒗 𝒙 𝑛−1|𝑛−1

システムノイズの⽣成(')

現時刻の状態の予測
(𝑗 )

𝑗

𝑗

𝑗

𝒙𝑛|𝑛−1 = 𝑭𝑛−1 𝒙 𝑛−1|𝑛−1 + 𝑮𝒗𝑛

(b)

(7)

(.)

現時刻の状態の予測式(33)の状態遷移⽅程式に従い，現時刻 𝑛 の状
態 𝒙! の予測分布 𝑝(𝒙! |𝑦! , 𝑦!$' , … ) の実現値とみなせる状態値を粒⼦ 𝑗 = 1, … , 𝑚 ごとに得る．

観測値に基づく尤度計算
𝑗

(7)

(")

j𝒗! , 𝒗! , … , 𝒗! , … , 𝒗! k , 𝒗! ~𝒩(0, 𝑸! )(𝑗 )𝑤𝑛 = 𝑝 𝒙𝑛|𝑛−1 𝑦𝑛

リサンプリングによる現時刻の状態の推定

(7)

(7)

(7)

(7)

𝒙!|!$' = 𝑭!$' 𝒙 !$'|!$' + 𝑮𝒗!
予測分布に従う粒⼦は 1 期前のフィルタの実現値から⽣成できるとみなすことができる．

1
𝑚
𝒙𝑛|𝑛
, … , 𝒙𝑛|𝑛Fig. 3 粒⼦フィルタでの状態推定アルゴリズム(c)

観測値に基づく尤度計算(7)2.3. ⾃⼰組織化状態空間モデルによる時系列予測

現時刻の観測値 𝑦! と各粒⼦の状態 𝒙!|!$' から

状態空間時系列モデルを利⽤して，観測値 𝑦! の将来値を予測するには，2.2.4.で導⼊した各時系列
成分とそのパラメータから構成される拡⼤状態
𝒙! を推定・予測する必要がある．先⾏研究 11)は時不変かつ線形・ガウス型の状態空間モデルを対象とすることから，状態推定・予測には⼀般的なカル
マンフィルタが⽤いられた．⼀⽅，本研究では時変パラメータを状態に取り込んだ⾃⼰組織化型の状態空間モデルを対象とした推定・予測が求められる．そこで本研究では⾃⼰組織化状態空間モデル
を⽐較的容易に扱うことのできる粒⼦フィルタ（Particle filter）を⽤いて時系列成分とそのパラメータの状態推定・予測する．

式(33)の観測⽅程式より予測された観測値 𝑦!|!$'に基づき，各粒⼦の尤度を計算する．

(7)

(7)

(d)

(#)

基 づ き ， ス テ ッ プ (b) の 予 測 分 布 実 現 値(#)

j𝒙!|!$' k , 𝑖 = 1,2, … , 𝑚 から 𝑚 回のリサンプリング（復元抽出）を⾏い，フィルタ分布の実現値

(7)

(&)

($)
:"

(')

;:" +⋯+:" =(#)

𝒙!|!$' がの確率で割り当てられる．

Fig.3 にアルゴリズムを図⽰したものを⽰す．
2.3.2. 時系列成分（状態）の多段予測
現時刻 𝑛 よりも 1 期先の予測分布に従う粒⼦(7)(7)
𝒙!+'|! は，時刻 𝑛 のフィルタの実現値 𝒙 !|! に遷移(7)
⾏列 𝑭! をかけることで表される．
(7)(7) (7)(34)
𝒙!+'|! = 𝑭! 𝒙 !|!
1 期先予測を繰り返す事で数期先の状態の予測も可能である．ただし，1 期先の観測値 𝑦!+' は未知であるため，前節の(c)観測値に基づく尤度計算と，それに続く(d)リサンプリングを⾏うことはで(7)

きない．そこで 1 期先予測分布の実現値 𝒙!+'|! を
時刻 𝑛 + 1 のフィルタの実現値とみなすことで 2
期先予測値を次式により得る．

粒⼦ 𝑗 = 1, … , 𝑚 について適当な乱数を⽣成し，各状態の初期分布の実現値とみなせる 𝑚 個の値を
得る．
(7)

(7)

(7) (7)(35)
𝒙!+"|! = 𝑭! 𝒙!+'|!
同様に，ℎ 期先の状態の予測値は以下のように表
せる．

(.)

j𝒙4 , 𝒙4 , … , 𝒙4 , … , 𝒙4 k
(2)

(.)

ィ ル タ 実 現 値 𝒙!|! に は 予 測 値

初期分布の⽣成

(")

(7)

(')

j𝒙!|! , … , 𝒙!|! , … , 𝒙!|! kを得る．リサンプリングにより，𝑗 番⽬の粒⼦のフ粒⼦フィルタは，時刻 𝑛 − 1 までの観測済みの値を も と に 状 態 𝒙! を 予 測 し た 分 布𝑝(𝒙! |𝑦!$' , 𝑦!$" , … ) と，時刻 𝑛 で新たに観測された値をもとに状態 𝒙! の予測分布を更新したフィ
ルタ分布 𝑝(𝒙! |𝑦! , 𝑦!$' , … ) とを独⽴に⽣成した 𝑚
個の粒⼦を⽤いて近似する．このことは，それぞれの分布関数を 𝑚 個の粒⼦から定義される経験分布関数によって近似していることに相当する．以下
では，2.2.4 で説明した適応的 2 次 AR＋トレンド
成分モデルを例として粒⼦フィルタでの状態推定のアルゴリズムを説明する 17)．

(')

リサンプリングステップ(c)で求めた尤度 j𝑤! k , 𝑖 = 1,2, … , 𝑚 に

2.3.1. 粒⼦フィルタによる状態推定

(1)

(7)

𝑤! = 𝑝 m𝒙!|!$' n𝑦! o

𝑛 > 0 の各時刻で以下の(a)から(d)のステップを繰り返す．

(7)(7) (7)𝒙!+>|! = 𝑭! 𝒙!+>$'|!

6

(36)

Fig. 4 ⼈⼯時系列データ

Fig 5. 呼吸運動データ
Fig.5 に本研究で⽤いた⼈⼯時系列と腫瘍位置変動
時系列をそれぞれ⽰す．

2.3.3. 観測値の予測
将来の時刻 𝑛 + ℎ の観測値 𝑦!+> の多段予測値(')

3.2. ⽐較予測モデルと条件(.)の分布 j𝑦!+>|! , … , 𝑦!+>|! k は式(37)の多段予測に(')

実験に⽤いた予測⼿法とその設定を下記に⽰す．

(.)

よる状態予測値 j𝒙!+>|! , … , 𝒙!+>|! k へ観測⾏列をかけることで求めることができる．
(7)

(7)

𝑦!+>|! = 𝑯𝒙!+>|!

(1)

適応的 AR+トレンド成分モデル（提案法）

パラメータを時間変化させる⾃⼰組織化状態空間モデルに基づき粒⼦フィルタで予測を⾏う提案法である．使⽤する粒⼦数は 𝑚 = 1000 個とした．
モデルパラメータ初期値には予測性能評価区間(37)

3. 予測実験
提案する適応的 AR+トレンド成分モデルによる
腫瘍位置変動時系列の予測性能を評価するため，
⼈⼯時系列と実データの 2 つを⽤いた予測実験を⾏った．

の直前 30 秒（ 30 s × 10 Hz = 300 サンプル）の観測値を⽤いた．AR 成分モデルの次数 𝑞 は最⼤次数𝑞./? = 10として AIC により決定し，Burg 法による最⼩⼆乗推定値を使⽤して AR 係数の初期値
𝑎',4 , … , 𝑎%,4 を求めた．また，システムノイズの分
散共分散⾏列の初期値 𝑸4 と観測ノイズの分散共分散⾏列 𝑹 は，Burg 法による最⼩⼆乗推定値の残差項の分散𝜎3" をもとに以下のように設定した．

3.1. 予測対象データ
はじめに，呼吸波形を模した振幅と周期がランダムに変化する⼈⼯時系列データを⽤いて予測性能の評価を⾏った．⼈⼯時系列の⽣成には呼吸波
形の模擬に使われる cos4(⋅) 関数を⽤いた．デー
タの⻑さは 200 秒，サンプリング周波数は 10 Hz
である．次に，http://signals.rob.uni-luebeck.de においてオンライン公開されている腫瘍位置変動時
系列データベースのうちの 1 つを⽤いて予測性能の評価を⾏った 19)．サンプリング周波数を 26 Hz
から 10 Hz にダウンサンプリングし，そのうち692.3 秒分を使⽤した．なお，本データベースは⽶
国 Georgetown University Hospital に お い てCyberKnife Synchrony システムにより計測された治療中の肺がんの腫瘍位置変動時系列（肺癌患者31 名 304 例分）から構成されている．Fig.4 および

𝜎"
⎡ 3
⎢2
⎢
⎢
𝑸4 = ⎢
⎢
⎢
⎢
⎣

𝜎3"
2

10$@

10$@

10$@

⎤
⎥
⎥
⎥
⎥,
⎥
⎥
⎥
10$@ ⎦

(38)

𝜎3"
P
4
AR+トレンド成分モデル（従来法）
𝑹=O

(2)

時不変の状態空間モデルに基づき粒⼦フィルタ
で予測を⾏う従来法 11)に相当する⽅法である．使7

⽤する粒⼦数は 𝑚 = 1000 個とした．
モデルパラメータは提案法と同様に予測性能評
価区間の直前 30 秒の観測値より設定した．すなわ
ち，AR 成分モデルの次数は AIC に基づき𝑚 ≤ 10 ．
で設定し，AR 係数は Burg 法による最⼩⼆乗推定値を⽤いた．システムノイズの分散共分散⾏列 𝑸と観測ノイズの分散共分散⾏列 𝑹 は以下のとおり
である．
𝜎"
⎡ 3
𝑸=⎢2
⎢
⎣

(3)

⎤
𝜎3"
⎥
,
𝑹
=
O
P
4
𝜎3" ⎥
2⎦
予測による遅れ補償なし

最⼩であった．Fig. 7 に拡⼤して⽰した 35 秒間は周期・振幅の変動が強く⾒られる区間である．この区間に着⽬すると，従来法は頂点付近では予測値が観測値と乖離があり，変動に対応できていない
傾向がみてとれる．⼀⽅，提案法は，頂点付近の予測値は従来法と⽐べて観測値に近い値となっていることが確認できる．Fig. 8 は上から AR 係数
𝑎' , 𝑎" , 𝑎B の時間変化，AR 成分のシステムノイズ分
散 𝜂" の時間変化，トレンド成分のシステムノイズ
分散 𝜏 " の時間変化を⽰している．周期・振幅が変化した点に応じてパラメータの値が変化していることが確認される．このことから，⾃⼰組織化状態
空間モデルのパラメータ同時推定により，時系列の性質変化への適応がなされていると考えること
ができる．

(39)

予測による遅れ補償なしに照射を実施したときに相当するものとして零次ホールド法（zero-order
hold: ZOH）による予測も⽐較対象とした．最新の値をそのまま⽤いることで，ℎ 期分ずれた同じ波形が出⼒される．
𝑦x!+>|! = 𝑦!

3.4.2. 腫瘍位置変動時系列データでの予測結果
腫瘍位置変動時系列における各予測モデルによる予測結果を Fig. 9-12 に⽰す．Fig. 9 は予測結果
全体を，Fig. 10 は予測結果のうち 600 秒から 635
秒までの 35 秒間をそれぞれ⽰している．図中，⿊の細線は観測値，⾚の実線と⻘の波線は提案法と(40)3.3. 予測性能評価指標
予測性能の評価には予測された観測値の Root
Mean Square 誤差（RMSE）を⽤いた．RMSE は次式で表される．

(7)

従来法の予測値（予測分布j𝑦!+>|! kの平均），緑の点
線は遅れ補償なしの予測値である．また，Fig. 11は 600 秒から 635 秒までの 35 秒間における提案
法のパラメータの時間変化を⽰している．図中，⽔⾊の実線は提案法のパラメータの推定値の平均，
灰⾊の破線はパラメータの推定値の標準偏差である．
Fig. 9 より，予測波形全体ではいずれの⽅法の予測波形に⼤きな差はみられないが，予測誤差は従来法が RMSE=0.636, 提案法が RMSE＝0.549,
遅れ補償なしが RMSE=0.867 となり，提案法が最⼩であった．
Fig. 10 の 600 秒から 635 秒までの 35 秒間は
特に周期の変動が⾒られる区間である．この区間に着⽬しても，従来法と提案法の予測波形に⼤きな差は⾒られず，どちらも頂点と⾕を捉えることができていないことがわかる．
Fig. 11 は上からこの区間における AR 係数
𝑎' , 𝑎" , 𝑎B の時間変化，AR 成分のシステムノイズ
分散 𝜂" の時間変化，トレンド成分のシステムノ
イズ分散 𝜏 " の時間変化を⽰している．624 秒付
近の周期が急激に短くなるような点でいずれのパラメータの値も⼤きく変動していることが⾒て取れ，時系列の性質が変化するような点でパラメータも時間変化していると⾔える．
さらに，Fig. 12 に Fig. 9-10 の 611 秒から 616
秒までの 5 秒間を拡⼤して⽰す．この区間に着⽬すると，腫瘍変動時系列には 1 秒程度の周期をも

A

1
RMSE = } )(𝑦# − 𝑦x# )"
𝑁

(41)

#&'

ただし，𝑦# は観測値， 𝑦x# はその予測値， 𝑁 は予測
性能評価区間に含まれる観測データ数である．
3.4. 予測実験結果
各モデルによる 0.3 秒先（ℎ = 3 期先）予測を時系列の後半 60%にわたって実施した．
3.4.1. ⼈⼯時系列データでの予測結果
それぞれの予測モデルについての予測結果をFig. 6-8 に⽰す．Fig. 6 は予測結果全体を，Fig. 7
は予測結果のうち 110 秒から 145 秒の 35 秒間を⽰している．図中，⿊の細線は観測値，⾚の実線と⻘の破線は提案法と従来法の予測値（予測分布(7)j𝑦!+>|! k の平均 𝑦•!+>|! ，緑の点線は遅れ補償なしの予測値である．また，Fig. 8 は 110 秒から 145
秒の 35 秒間における提案法のパラメータの時間
変 化を⽰している．図中，⽔⾊の実線は提案法のパラメータの推定値の平均，灰⾊の破線はパラメ
ータの推定値の標準偏差である．
Fig. 6 より，予測波形全体ではいずれの⽅法の予測波形に⼤きな差はみられないが，予測誤差は提案法が RMSE＝0.094, 従来法が RMSE=0.106，
遅れ補償なしが RMSE=0.149 となり，提案法が8

Fig 6. ⼈⼯時系列データの予測波形（全体図）

𝜏2

𝜂2

Fig 7. ⼈⼯時系列データの予測波形（110 秒〜145 秒までの拡⼤図）

Fig 8. ⼈⼯時系列データにおける提案法のパラメータの時間変化の推定値
上から，AR 係数𝑎' , 𝑎" , 𝑎B ，AR 成分のシステムノイズ分散 𝜂" ，トレンド成分のシステムノイズ分散 𝜏 "9Fig 9. 呼吸運動データの予測波形（全体図）

𝜏2

𝜂2

Fig 10. 呼吸運動データの予測波形（600 秒〜635 秒までの拡⼤図）

Fig 11. 呼吸運動データにおける提案法のパラメータの時間変化の推定値
上から，AR 係数𝑎' , 𝑎" , 𝑎B ，AR 成分のシステムノイズ分散 𝜂" ，トレンド成分のシステムノイズ分散 𝜏 "Fig 12. 呼吸運動データの予測波形（611 秒〜616 秒までの拡⼤図）10

つ拍動による振動的な動きも含まれていることがわかる．従来法の予測値にも，1 秒程度の周期の振動的な動きがみられるが，観測値に⽐べてはっきりせず，またピークの位置もずれていることが観察できる．⼀⽅，提案法による予測値の拍動ピークは，観測値の拍動ピークとタイミングがほぼ⼀致しており，より適切に拍動成分を表現できている．

ベースに含まれる他のデータを⽤いて，近年の深層学習に基づく⼿法など他の予測⼿法との性能⽐
較が今後の課題である．

参考⽂献
1)

年次推移：[国⽴がん研究センター がん統計].”
https://ganjoho.jp/reg_stat/statistics/stat/an
nual.html (accessed Nov. 24,
2) K. Takahiro and I. Kiyonari, 診療放射線技術
下巻, 改訂第 14 版. 南光堂, 2019.
3) Y. Suh, S. Dieterich, B. Cho, and P. J. Keall, “Ananalysis of thoracic and abdominal tumour
motion for stereotactic body radiotherapy
patients,” Phys Med Biol, vol. 53, no. 13, p.
3623, Jun. 2008, doi: 10.1088/00319155/53/13/016.
4) P. J. Keall et al., “The management of
respiratory motion in radiation oncology reportof AAPM Task Group 76,” Med Phys, vol. 33,
no. 10, pp. 3874‒3900, 2006, doi:
10.1118/1.2349696.
5) P. R. Poulsen, B. Cho, A. Sawant, D. Ruan, and
P. J. Keall, “Detailed analysis of latencies inimage-based dynamic MLC tracking,” MedPhys, vol. 37, no. 9, p. 4998, 2010, doi:
10.1118/1.3480504.
6) P. J. Keall et al., “AAPM Task Group 264: The
safe clinical implementation of MLC tracking
in radiotherapy,” Med Phys, vol. 48, no. 5, May
2021, doi: 10.1002/mp.14625.
7) A. Jöhl et al., “Performance comparison ofprediction filters for respiratory motion
tracking in radiotherapy,” Med Phys, vol. 47,
no. 2, pp. 643‒650, Feb. 2020, doi:
10.1002/MP.13929.
8) N. Fukumitsu et al., “Reproducibility of image
quality for moving objects using respiratorygated computed tomography: a study using aphantom model,” J Radiat Res, vol. 53, no. 6,
pp.
945‒953,
Nov.
2012,
doi:
10.1093/JRR/RRS039.
9) ⻄⼭ほか, 注意機構に基づく肺腫瘍位置変動
予測, 第 55 回⽇本⽣体医⼯学会東北⽀部⼤会,
2021
10) Toftegaard et al., “Potential improvements oflung and prostate MLC tracking investigatedby treatment simulations,” Med. Phys, vol. 34,
no. 1, p. 7, 2018, doi: 10.1002/mp.12868.
11) 佐藤ほか, 腫瘍位置変動対策のための肺腫瘍

3.5. 考察⼈⼯時系列データにおいては提案した適応的
AR+トレンド成分モデルは従来法である AR+トレ
ンドモデル，遅れ補償なしと⽐較して優れた予測性能を⽰し，特に振幅や周期が変化するような区間での予測誤差低減効果がみられた．実際の腫瘍
位置変動データにおいても，提案した適応的 AR+
トレンド成分モデルの予測誤差指標は最⼩となったが，呼吸振幅や周期が変化するような区間において従来法との⼤きな差は⾒られなかった．⼀⽅，
提案法は拍動による振動的な動きを従来法よりも適切に予測・表現することができていた．これらのことからは，呼吸性に加えて拍動性を含む場合には，パラメータ時間変化を導⼊しても AR＋トレン
ド成分モデルのみでは，これら２種類の周期的変動成分を⼗分に表現しきれなかった可能性が考えられる．今回，提案法の AR 成分モデルが拍動性の変動を従来法よりもよく捉えられていた点を考慮すると，呼吸性変動の表現・予測を AR 成分モデルとは別の時系列成分モデルに分担させる⽅策が考えられる．呼吸性変動の時系列表現としては，時変
季節性⾃⼰回帰モデル 10)などが提案されており，
このような時系列成分を導⼊して時変状態空間時
系列モデルを構築することで，より適切に腫瘍位置変動の性質変化を表現でき，予測性能が向上する可能性があると考えられる．

4. おわりに本研究では，⾼精度放射線治療のための腫瘍位
置変動の予測性能向上のために，パラメータを時間変化させる時変状態空間モデルによる予測を試みた．パラメータを時間変化させない従来法との予測性能⽐較実験の結果より，⼈⼯時系列データ
ではパラメータの時間変化によって予測性能向上が確認されたが，実際の腫瘍位置運動データでは予測性能向上はわずかであった．呼吸と拍動の２種類の周期的変動を少ない成分モデルで表現しきれなかったことが原因と考えらえる．
AR＋トレンド成分に加えて，時変季節性成分など今回考慮していない時系列成分を導⼊すること
で更なる予測性能向上を図ることと，公開データ

11

位置の時系列成分分離に基づく予測,⽇本放射線技術学会（JSRT）
第 46 回秋季学術⼤会, 201812) K. C. McCall and R. Jeraj, “Dual-component
model of respiratory motion based on the
periodic autoregressive moving average(periodic ARMA) method,” Phys Med Biol, vol.
52, no. 12, p. 3455, May 2007, doi:
10.1088/0031-9155/52/12/009.
13) A. Mylonas, J. Booth, and D. T. Nguyen, “A
review of artificial intelligence applications formotion tracking in radiotherapy,” J Med
Imaging Radiat Oncol, vol. 65, no. 5, pp. 596‒
611, Aug. 2021, doi: 10.1111/17549485.13285.
14) F. Ernst, A. Schlaefer, and A. Schweikard,“Prediction of Respiratory Motion with
Wavelet-Based Multiscale Autoregression,”
Lecture Notes in Computer Science, vol. 4792
LNCS, no. PART 2, pp. 668‒675, Jan. 2007,
doi: 10.1007/978-3-540-75759-7_81.
15) F. Ernst, A. Schweikard, F. Ernst, · A
Schweikard, and A. Schweikard, “Forecasting
respiratory motion with accurate online
support vector regression (SVRpred),” Int J
CARS, vol. 4, pp. 439‒447, 2009, doi:
10.1007/s11548-009-0355-5.
16) 北川源四郎, 時系列解析⼊⾨. 岩波書店, 2005.
17) G. Kitagawa, “A Self-Organizing State-Space
Model A Self-Organizing State-Space Model
Genshiro KITAGAWA,” Source: Journal of theAmerican Statistical Association, vol. 93, no.
443, pp. 1203‒1215, 1998.
18) 樋⼝知之, 予測に活かす統計モデリングの基
本, 改訂第 2 版. 講談社, 2022.
19) F. Ernst, Compensating for Quasi-periodic
Motion in Robotic Radiosurgery. New York,
NY: Springer New York, 2012. doi:
10.1007/978-1-4614-1912-9.

12