MEGデータの解析において,従来は100回程度の試行を加算平均することで意味のある信号を得ることが一般的であったが,近年のBrian computer interface(BCI)技術の進展に伴い,単一試行データからの信号抽出の重要性が高まっている.BCIへの応用のみならず,MEGデータのノイズ除去はMEG解析の前処理として非常に重要なものでもある.

本論文は,MEG計測データの時間構造に着目し,MEGノイズを推定・除去する方法の提案と検討を行った.第一章では,MEG計測の背景を説明し,ノイズに対処するための主要な既存手法の概略と問題点を概説した.第二章はMEG計測におけるノイズの性質について説明し,これまでの研究を踏まえ提案した新しいノイズ除去手法の概略を説明した.さらにその手法の有効性を示唆するための数値実験について記述した.第三章では,現実に即したより自然な信号モデルの導入とその予備的数値実験を示した.第四章では,実ノイズと2次の線形モデルで生成した誘発脳磁場信号を合成したデータを用いて提案手法の有効性を検討した.さらに実ノイズとMEGデータの同期加算によって得られた信号を合成したデータを用いて精度のよい信号成分の推定ができることを示した.第五章は,人為的な手を全く加えない単一試行MEGデータに提案手法を実行した結果を示した.第六章では,単一試行MEG解析に必要な計算資源の量について検討した.また,我々の手法は効率の良い並列処理が可能であることを示し,それを実現するためのGraphics Processing Unit(GPU)を用いた実装について概説した.第七章では,これまでのまとめと考察,将来の展望について述べた.以下各章について詳しく述べていく.

MEGデータのノイズ除去法は大きく4つに分類できる.参照チャネル等の異なる観測データを利用する手法,MEG計測系の幾何的配置の情報を利用する空間フィルタ,MEGデータの持つ時間構造を利用した手法,多変量解析の応用である.第一章ではこれら4クラスでの代表的な手法を挙げて,その有効性,問題点を概説した.

第二章では,これまでの既存のMEGノイズの除去法を踏まえて,それらの手法の問題点を解決する新しい手法を提案した.その手法は,MEG計測データの生成モデルを観測データから推定し,その最適推定モデルを用いて,MEGデータの推定値を得る方法である.すなわち,まずMEGデータの生成モデルを次のように与える.

ここではあるチャンネルにおける時刻でのMEG計測値であり,は誘発脳磁場成分,は1Hz未満のベースラインドリフト成分,は交流電源由来の定常周期成分,はls次元の線形システムから生成される背景雑音成分は確率的に揺らぐセンサ雑音及びモデル化誤差成分であり,平均 ,分散のガウス過程である.ここで,は次に示すような自己回帰過程であると仮定する.

ただし,vtは平均 ,分散σ2e ,のガウス過程である. tt,Pt 成分は特殊な自己回帰過程で記述できる.提案手法は,ターゲットとする脳活動が含まれていないデータのプレ区間(stが0の区間)を利用して,背景雑音モデルを学習する.最適モデルを決定するときの評価関数として,赤池情報量基準(AIC)を用いる.その定義式は次式である.

ここでp(m1,...,mt|Θ)は,モデルパラメータΘが与えられたときの対数尤度である.また,MはモデルパラメータΘである.ここでのモデルパラメータベクトルは具体的には, etの回帰係数{ai,...al}とその入力となるガウス過程の分散σ2e ,センサノイズの分散σ2nから構成される.一般に条件付き同時確率密度p(m1,...,m|Θ)の評価は難しいが,本手法のような時間構造を仮定することにより成立する観測値mtのマルコフ性を利用できて,次のように書き下すことができる.

ここで, Miは観測時系列{mi,...,mi}を表す.この式の右辺は,カルマンフィルタを用いることで非常に効率的に評価できる[1].最適な背景雑音の近似モデルを探索するには,AICを評価関数としてグリッドサーチ,非線形最適化を順次行い,モデルパラメータとその次数を決定すればよい.これらの背景雑音モデルを学習した後に,信号を含む区間データにおいて,同様に信号モデルを推定する.ここでは単純にAR過程を仮定し,で与える.ここで,wtは平均0 ,分散σ2sのガウス過程である.

この手法を確かめるためにシミュレーションを行った.利用したデータセットは実ノイズデータに信号を加算したものである.信号は20Hzと70Hzにピークを持つようなARモデルに白色ガウス雑音を入力として与えた時の出力を,電流双極子のモーメントとし,その磁束を求めたものを使用した.また,本手法の性能を評価するために,Fast ICAを用いたノイズ除去法と比較を行った [2].この結果としては,提案手法はICAと比較してSN比の改善がよいことが示された(図1).また,本手法は各チャンネルを独立に処理する線形フィルタであるのでICAモデルを変化させないことから,本手法の後処理としてICAを適用できる.これにより,さらなるSN比の改善ができることを示した.さらに推定されたノイズモデルのパワースペクトルのconsistencyから提案した背景ノイズモデルが適切であることを示した.

第三章では,実際の誘発脳磁場の近似モデルとして,定常AR過程は適切でないということを示し誘発脳磁場モデルの拡張を行った.定常AR過程のパワースペクトルp(f)とモデルパラメータとの関係はで書けるので,定常AR過程では誘発脳磁場のパワーの時間変化を適切に扱えないことがわかる.そこで,誘発脳磁場の近似モデルを, ARモデルへの入力の確率構造が時変である非定常モデルに拡張した.ここで導入した動的に変化するより複雑なモデルを推定するために,演算効率の観点からパーティクルフィルタ(PF)を利用した[3][4].PFは一般の確率分布を、その分布の標本の集まりで近似を行い,その時間発展と新たな観測値を用いたフィルタ分布を逐次的に計算するアルゴリズムであり,並列計算に向いたものである.そこで,このPFのGPU上での実装を行った.

第四章では,本手法の有効性を検証するため,2次の線形モデルで生成した誘発脳磁場信号と実ノイズを合成したデータから,信号を抽出出来ることを示した(図2,3).続いて被験者がボタン押しをしている際に計測されたMEGの同期加算信号と実ノイズを合成したデータから信号を抽出出来ることを示した.先行する手法との比較も行った.

第五章では,実際のMEG計測によって得られた単一試行データに対して本手法を適用した.本手法においては自発活動のモデルを導入していないため,誘発活動の推定には条件が悪いが,それでも10試行程度の同期加算を行うことで,明瞭な誘発脳磁場を得ることができた.

第六章では,本手法をMEGデータに適用するときに問題となる,演算の効率性について考察を行った.GPUを用いて非常に高度な並列化を行うことで,Matlab上での実装と比較して約1000倍高速に実行する.ことに成功した.この実装において重要な,累積分布関数の算出,尤度による重み付けに基づくリサンプリングを行う部分での並列実装の要点を記述した.

第七章では以上二つの数値実験のまとめと総合考察,今後の展望を述べた. MEGデータの時間構造を活用することで,単一試行データからでも誘発脳磁場の検出が可能であることを数値実験により示した.また,提案手法と既存の手法を組み合わせることで,ノイズ除去性能の向上が図れることを数値実験により示した.同期加算平均によって得られた誘発脳磁場を用いた数値実験でも単一試行から誘発脳磁場の検出が十分に可能である.ことを示した.また,自発活動のモデルを未導入にも関わらず,本手法を適用後,10回程度の同期加算を行えば,100試行を用いて求められた同期加算波形成分を検出することができた.

本手法は,ノイズや誘発脳磁場に対する解析者の信念をモデルの構造に込めることにより,観測データを最もよく説明するモデルを得ることで,各成分の分解を得る方法である.我々の研究の意義は,完全とは言えないモデルでも,単一試行MEGデータから微小な誘発脳磁場成分を抽出できる可能性が十分にあることを示した点にあると考える.単一試行MEGデータ中の脳磁場については未知な部分が多いが,MEG研究者が経験的に得ている仮説を今後種々取り入れることで,それらの仮説の検証が可能になるはずである.断続的なモデルの改良を行っていくことが,MEG研究の目標である脳という未知システムを理解する有力な手段であることを示せた点も,我々の研究の成果である.

図1 解析結果(28ch分を表示,全ての図で同一スケール),左上:原データ,右上:用いた信号,左下:提案手法による処理結果,右下:ICAによる処理結果

図2 定常ARモデルと時変分散ARモデルの推定

緑線が定常ARモデルの推定,青線が時変分散ARモデルの推定

図3 磁場マップにおける比較

左上:作成した信号 左下:真の信号 中:推定に用いたデータ 右上:推定されたノイズ成分の和 左下:推定された信号

図4 コヒーレンスマップ

左図は今回用いた100試行に渡る,真の信号とのコヒーレンス関数を信号の帯域で平均したもの.右図は真の信号の試行間平均のパワーを正規化したもの.信号のパワーが強い領域では非常に高いコヒーレンスが得られている.

本論文は全7章よりなり,MEG計測データの時間構造に着目し,MEGノイズを推定・除去する方法の検討と提案を行った。

第1章では,MEG計測の背景を説明し,MEG計測におけるノイズの性質について説明し,これらのノイズに対処するべく提案された主な既存手法の概略とその有効性と問題点を概説した。MEGデータのノイズ除去法は大きく4つに分類できる。参照チャネルを利用する手法,空間フィルタ,時間構造を利用したフィルタ,多変量解析の応用である。第一章ではこれら4クラスでの代表的な手法を説明し,その有効性,問題点を概説した。

第2章はこれまでの研究を踏まえ提案した新しい手法を述べた。その手法は,MEG計測データの生成モデルを観測データから推定し,その最適推定モデルを用いて,MEGデータの推定値を得る方法である。すなわち,まずMEGデータの生成モデルを次のように与える。

ここでmtはあるチャンネルにおける時刻tでのMEG計測値であり,は誘発脳磁場成分,etはls次元の線形システムから生成される背景雑音成分,ntは確率的に揺らぐセンサ雑音及びモデル化誤差成分であり,平均0 ,分散σ2nのガウス過程である。ここで, etは次に示すような自己回帰過程であると仮定する。

ただし, vtは平均0 ,分散σ2g ,のガウス過程である。提案手法は,ターゲットとする脳活動が含まれていないデータのプレ区間(stが0の区間)を利用して,背景雑音モデルを学習する。最適モデルを決定するときの評価関数として,赤池情報量基準(AIC)を用いる。その定義式は次式である。

ここでP(m1,...,mt | Θ)を,モデルパラメータΘ が与えられたときに観測時系列{m1,...,mt}が観測される条件付き同時確率密度関数とし,これは時系列モデルの対数尤度と等しい。また,Mはモデルパラメータ数である。ここでのモデルパラメータベクトルΘは具体的には,etの回帰係数{ai,...,al}とその入力となるガウス過程の分散σ2e ,センサノイズの分散σ2nから構成される。最適な背景雑音の近似モデルを探索するには,AICを評価関数として非線形最適化を行い,モデルパラメータとその次数を決定すればよい。これらの背景雑音モデルを学習した後に,脳磁場を含むデータにおいて,同様に誘発脳磁場のモデルを推定する。誘発脳磁場の構造はAR過程を仮定し,で与える。ここで, Wtは平均0 ,分散σ2sのガウス過程である。

第3章は数値実験により提案手法の有効性を検討した。計算機で作成した誘発脳磁場成分を,MEGノイズデータと合成したデータセットを作成し,提案手法でのノイズ除去の結果と,独立成分分析(ICA)との結果を比較した。作成したデータセットは、実ノイズデータに擬似生体信号(以後,信号と記述する)を加算して作成した。信号は20Hzと70Hzにピークを持つようなARモデルに白色ガウス雑音を入力として与えた時の出力を,電流双極子のモーメントとし,その磁束を求めたものを使用した。この結果,提案手法はICAと比較して全般にSN比の改善がよいことが示された。

第4章では,実際のMEG計測で得られた誘発脳磁場を利用して,提案手法の有効性を検討した。その結果から我々の手法の問題点を指摘し,それを解決するために手法の拡張を行った。一般に誘発脳磁場のパワーは大きく時間変化している。定常AR過程のパワースペクトルとモデルパラメータとの関係は 〓で書けるので,定常AR過程では誘発脳磁場のパワーの時間変化を適切に扱えないことがわかる。そこで,誘発脳磁場の近似モデルとして, ARモデルの入力のガウス過程の分散が時変であるモデルσ2s,tに拡張した。さらに,この分散σ2s,tが急変することは現実的でないので、モデルに平滑化制約を課した。

第5章では,提案手法を実行するうえで問題となる演算時間について検討した。上述のような動的に変化する、より複雑なモデルパラメータを推定するために,演算効率の観点からパーティクルフィルタ(PF)を利用した。PFは一般の確率分布を、その分布の標本の集まりで近似を行い,その時間発展と新たな観測値を用いたフィルタ分布を逐次的に計算するアルゴリズムであり,並列計算に向いたものである。このPFの演算の処理時間を短縮させるために,Graphics Processing Unit(GPU)を用いたプログラムを作成し,超並列演算の概要について説明した。

第6章は,新たな提案手法を用いて,数値実験を行いその性能を評価し,その結果について検討を行った。推定に用いた信号にはパワーの時間変動が大きい二つの電流双極子によって生成される誘発脳磁場様の信号を用いた。その結果から磁場マップを再構成しダイポール推定を行った結果,真値(-4,5,4)cm,(-4,-5,4)cmに対し,(-4.59,5.43,4.34)cm,(-4.90,-5.69,4.72)cmにGOF=79.89で推定できた。推定位置の誤差はそれぞれ0.84cm,1.15cmとなり,単一試行データでもかなり良い推定ができることを示した。

第7章では,研究のまとめと考察,将来の展望について述べた。MEGデータの時間構造を活用することで,単一試行データからでも誘発脳磁場の検出が可能であることを数値実験により示した。また,提案手法と既存の手法を組み合わせることで,ノイズ除去性能の向上が図れることを数値実験により示した。今後の展望としては,加算平均によって得られた誘発脳磁場を用いた数値実験を行い,その有効性を確認すること,その次に実際のMEGデータを単一試行で解析することが挙げられる。

なお、本論文第2章と第6章は、武田常広との共同研究であるが、論文提出者が主体となって分析及び検証を行ったもので、論文提出者の寄与が十分であると判断する。したがって、博士(科学)の学位を授与できると認める。