特許 5852331 データ処理装置、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】5852331

(24)【登録日】2015年12月11日

(45)【発行日】2016年2月3日

(54)【発明の名称】データ処理装置、及びプログラム

(51)【国際特許分類】

   A61B 5/05 20060101AFI20160114BHJP

【ＦＩ】

   A61B5/05 A

【請求項の数】5

【全頁数】19

(21)【出願番号】特願2011-124688(P2011-124688)

(22)【出願日】2011年6月2日

(65)【公開番号】特開2012-249817(P2012-249817A)

(43)【公開日】2012年12月20日

【審査請求日】2014年5月30日

(73)【特許権者】

【識別番号】513155677

【氏名又は名称】株式会社新領域技術研究所

(74)【代理人】

【識別番号】100122275

【弁理士】

【氏名又は名称】竹居 信利

(72)【発明者】

【氏名】武田 常広

(72)【発明者】

【氏名】宇野 裕

(72)【発明者】

【氏名】天野 薫

【審査官】 樋熊 政一

(56)【参考文献】

【文献】 国際公開第２０１０／１０７０９３（ＷＯ，Ａ１）

【文献】 特開平０８−１１２２６４（ＪＰ，Ａ）

【文献】 特開平１１−００９５６７（ＪＰ，Ａ）

【文献】 特開平１０−２９０７９２（ＪＰ，Ａ）

【文献】 特開２００９−２９７３８２（ＪＰ，Ａ）

【文献】 特開２００８−１７８５４６（ＪＰ，Ａ）

【文献】 特開２００５−２８７６７５（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ａ６１Ｂ ５／００−５／０５３

Ａ６１Ｂ ５／０６−５／２２

(57)【特許請求の範囲】

【請求項1】

測定された、ノイズと自発脳磁場とを含んだ時系列の安静時データと、ノイズと自発脳磁場と誘発脳磁場とを含んだ時系列の刺激時データとを受け入れる受入手段と、
少なくとも一つのノイズパラメータを用いてノイズを模式的に表すノイズ演算式に含まれる前記ノイズパラメータを、過去に得られたノイズのみのデータを用いて定められる基準に基づいて決定し、ノイズの時系列推定データを得る手段と、
前記安静時データの対応する時刻のデータから前記ノイズの時系列推定データを差引きして、ノイズがないと仮定される状態のデータを得る手段と、
自発脳磁場に含まれる信号を模式的に、自発脳磁場の成分の和として近似的に表すときの、当該自発脳磁場の成分を表すモデルパラメータを、前記ノイズがないと仮定される状態のデータに基づいて求め、当該自発脳磁場の成分を表すモデルパラメータと、自発脳磁場に含まれる信号を模式的に、自発脳磁場の成分の和として近似的に表すモデルとに基づいて、自発脳磁場の時系列推定データを得る推定手段と、
前記刺激時データからの対応する時刻のデータから、前記ノイズの時系列推定データと、前記自発脳磁場の時系列推定データとを差引きしてノイズ及び自発脳磁場がないと仮定される状態のデータを得て出力する手段と、
を含むデータ処理装置。

【請求項2】

請求項１記載のデータ処理装置であって、
前記推定手段は、時刻ｔにおける安静時データがｍtであったとき、時刻ｔ以前の自発脳磁場を表すデータ及び時刻ｔ以前のノイズを表すデータとを含めた状態ベクトルがｘtである確率分布ｐ（ｘt｜ｍt）を最大とする前記モデルパラメータを、時刻ｔで状態ベクトルがｘtである確率ｐ（ｘt）を事前確率とし、時刻ｔで状態ベクトルがｘtであったときに、前記自発脳磁場に含まれる信号を模式的に、自発の磁場の成分の和として近似的に表すモデルが表す時刻ｔでの安静時データの確率分布ｐ（ｍt｜ｘt）を事後確率とし、当該事前確率と事後確率の積を被積分関数とする積分結果を用いたベイズ推定により推定する推定手段であるデータ処理装置。

【請求項3】

請求項２記載のデータ処理装置であって、
複数の演算装置要素を含み、
前記推定手段は、前記ベイズ推定における、前記事前確率と前記事後確率との積を被積分関数とする積分を行うにあたり、積分範囲内からランダムに選択される積分変数の値を代入したときの前記被積分関数の値の演算を、複数の演算装置要素に分散的に行わせ、モンテカルロ法により前記積分を遂行するデータ処理装置。

【請求項4】

請求項１から３のいずれか一項に記載のデータ処理装置であって、
前記ノイズの時系列推定データを得る手段は、前記安静時データまたは刺激時データの実際の測定に先立って、前記受入手段が受け入れるデータを、前記ノイズのみのデータとして用いるデータ処理装置。

【請求項5】

コンピュータを、
測定された、ノイズと自発脳磁場とを含んだ時系列の安静時データと、ノイズと自発脳磁場と誘発脳磁場とを含んだ時系列の刺激時データとを受け入れる受入手段と、
少なくとも一つのノイズパラメータを用いてノイズを模式的に表すノイズ演算式に含まれる前記ノイズパラメータを、過去に得られたノイズのみのデータを用いて定められる基準に基づいて決定し、ノイズの時系列推定データを得る手段と、
前記安静時データの対応する時刻のデータから前記ノイズの時系列推定データを差引きして、ノイズがないと仮定される状態のデータを得る手段と、
自発脳磁場に含まれる信号を模式的に、自発脳磁場の成分の和として近似的に表すときの、当該自発脳磁場の成分を表すモデルパラメータを、前記ノイズがないと仮定される状態のデータに基づいて求め、当該自発脳磁場の成分を表すモデルパラメータと、自発脳磁場に含まれる信号を模式的に、自発脳磁場の成分の和として近似的に表すモデルとに基づいて、自発脳磁場の時系列推定データを得る推定手段と、
前記刺激時データからの対応する時刻のデータから、前記ノイズの時系列推定データと、前記自発脳磁場の時系列推定データとを差引きしてノイズ及び自発脳磁場がないと仮定される状態のデータを得て出力する手段と、
として機能させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、データ処理装置、及びプログラムに係り、特にノイズの影響除去に関する。

【背景技術】

【0002】

脳磁計（ＭＥＧ）等、微小信号を計測する装置では、外来の種々のノイズに信号が埋もれてしまうため、ノイズの除去が課題となっている。かかるノイズ除去の一手法として、従来より、複数回の測定結果を加算平均する方法が知られている。また特許文献１には、脳磁場信号に対して、心拍の発生時刻に同期させて加算平均を取ることにより、心磁によるノイズ信号を抽出する技術が開示されている。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特開２００９−１９５５７１号公報

【発明の概要】

【発明が解決しようとする課題】

【0004】

しかしながら、上記加算平均を用いるノイズの除去方法では、測定回数をなるべく多くしなければならない。ところが測定の対象によっては多数回の測定が困難である場合もある。また脳磁計等の場合、測定回数を増大させるほど被験者の負担が大きくなってしまう。

【0005】

本発明は上記実情に鑑みて為されたもので、単一の測定結果からでも信号の推定を可能とするデータ処理装置を提供することを、その目的の一つとする。

【課題を解決するための手段】

【0006】

上記従来例の問題点を解決するための本発明は、データ処理装置であって、測定されたデータを受け入れる受入手段と、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、を含むこととしたものである。

【発明の効果】

【0007】

本発明によれば、ノイズ及び信号を、パラメータを含んだ模式的な演算式として表現し、そのパラメータを推定するようにしたことで、単一の測定結果からでも信号の推定を可能とした。

【図面の簡単な説明】

【0008】

【図1】本発明の実施の形態に係るデータ処理装置の例を表す構成ブロック図である。

【図2】本発明の実施の形態に係るデータ処理装置の例を表す機能ブロック図である。

【図3】本発明の実施の形態に係るデータ処理装置による処理例を表すフローチャート図である。

【図4】従来の方法で得られた誘発脳磁場の波形と、これらの波形から得られた等磁界線図の例を表す説明図である。

【図5】ノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図の例を表す説明図である。

【図6】本発明の実施例に係るデータ処理装置が出力する誘発脳磁場の波形と、これらの波形から得られた等磁界線図の例を表す説明図である。

【図7】ノイズを含めて実際に得られる時系列データの別の例と、当該時系列データを基に得られた等磁界線図の別の例を表す説明図である。

【図8】本発明の実施例に係るデータ処理装置で推定して得た等磁界線図の例を表す説明図である。

【図9】本発明の実施例に係るデータ処理装置で推定して得た複数の等磁界線図と、従来の方法で推定された等磁界線図とのコヒーレンスの例を表す説明図である。

【図10】従来の方法によって推定された等磁界線図の例を表す説明図である。

【発明を実施するための形態】

【0009】

本発明の実施の形態に係るデータ処理装置を図面を参照しつつ説明する。本実施の形態のデータ処理装置１は、図１に例示するように、制御部１１、記憶部１２、操作部１３、出力部１４、及び入力インタフェース部１５を含んで構成される。以下では、具体的に、脳磁計（ＭＥＧ）の出力するデータ（時間変化する一連のデータ）を入力インタフェース部１５を介して受け入れる例について説明するが、本実施の形態のデータ処理装置１は、脳磁計のデータを処理するものに限られない。

【0010】

脳磁計のデータには、被験者のいない状態で脳磁計が出力するデータ（ノイズのみのデータ）と、被験者が安静な状態にあるときに脳磁計が出力するデータ（外的刺激を受けていない状態：安静時データと呼ぶ）と、被験者に対して外的刺激を与えたときに脳磁計が出力するデータ（ノイズ、自発脳磁場及び、刺激による誘発脳磁場が累算されたデータ：刺激時データと呼ぶ）とがある。

【0011】

制御部１１は、ＣＰＵ（Central Processing Unit）等のプログラム制御デバイスであり、記憶部１２に格納されたプログラムに従って動作する。この制御部１１は、入力インタフェース部１５から入力される、外部の測定装置（ここでは例えば脳磁計とする）により測定されたデータを受け入れる。また、この制御部１１は、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、ノイズパラメータを決定する。

【0012】

さらに制御部１１は、測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、測定されたデータと、ノイズ演算式とに基づいて、信号パラメータを推定し、当該推定した信号パラメータと、信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する。この制御部１１の詳しい処理の内容は後に説明する。

【0013】

記憶部１２は、メモリデバイス等であり、制御部１１によって実行されるプログラムを保持している。このプログラムは、例えばＤＶＤ−ＲＯＭ（Digital Versatile Disc Read Only Memory）等のコンピュータ可読な記録媒体に格納されて提供され、この記憶部１２に複写されたものであってもよい。また、この記憶部１２は制御部１１のワークメモリとしても動作する。

【0014】

操作部１３は、キーボードやマウスなどを含み、利用者の指示操作を受け入れて、制御部１１に当該指示操作の内容を出力する。出力部１４は、ディスプレイや、プリンタ、ネットワークインタフェース等であり、制御部１１から入力される指示に従って、情報を出力する。

【0015】

入力インタフェース部１５は、ＵＳＢ（Universal Serial Bus）やパラレルポート等のデータの入力を受け入れるデバイスであり、ここでの例では脳磁計に接続されて、脳磁計が出力するデータを受け入れ、当該受け入れたデータを制御部１１に出力する。

【0016】

本実施の形態では、脳磁計から入力されるデータは時々刻々と変化する（時間変化する）。そして時刻ｔにおけるデータｍ_tが、測定の対象たる信号ｓ_tとノイズとの和により表現されているものとする。すなわち、

【数1】

ここにノイズは、要因ごとに複数の成分の和として表現されており、（１）式ではτ_tがベースラインドリフト成分、ｑ_tが交流電源成分、ｅ_tがベースラインドリフトや交流電源成分以外の定常ノイズ成分、ｎ_tは、このようにモデル化したことによるデータからの差を表すモデル化誤差成分であるものとしている。また本実施形態では、不定期に生じる測定値のとび（例えばグリッチノイズ）、眼球運動や心拍等の生体から生じるノイズはないものとする。

【0017】

ノイズの各成分は、具体的に次のような形の演算式で表されるものとする。ベースラインドリフト成分τ_tの時間発展を、平滑化事前分布を用いた確率的トレンドモデルを採用し、

【数2】

と表す。

【0018】

ここでv_{τ ,t}はトレンドの変化を起こす外乱であり、平均ゼロ分散σ_τのガウス分布と仮定しておく。なお、この（２）式、及び以下の説明において時刻ｔ＋ｊは、ｔから予め定めた単位時間（例えばＴとする）を用いてｊ×Ｔだけ後の時刻を意味するものである。

【0019】

さらに（１）式におけるｑ_tの時間発展を、サンプリング周波数と交流周波数とから決まる自己回帰モデル（ARモデル）を用いて、

【数3】

と表す。

【0020】

ここでも同様に、v_{q ,t}は平均ゼロ分散σ_qのガウス分布と仮定する。またΔTはサンプリング周波数と交流電源の周波数から決定される、離散化された周期を表す。

【0021】

定常ノイズ成分ｅ_tの時間発展は次数ｌのＡＲモデルを仮定する。すなわち、

【数4】

とする。

【0022】

これにおいても、ｖ_{e ,t}は平均ゼロ分散σ_eのガウス分布と仮定する。モデル化誤差成分ｎ_tは、平均ゼロ分散σ_nのガウス分布と仮定しておく。

【0023】

本実施の形態の制御部１１は、ノイズ演算式として、（１）から（４）式の和、

【数5】

を用いる。このノイズ演算式に含まれるノイズパラメータを成分とするベクトルθ（モデルパラメータθベクトル）は、従って、

【数6】

となる。

【0024】

本実施の形態の制御部１１は、機能的には、図２に例示するように、データ記録部２１、モデル学習部２２、ノイズパラメータ管理部２３、信号推定部２４、及び信号出力部２５を含んで構成される。ここにデータ記録部２１は、入力インタフェース部１５を介して受け入れられるデータを記憶部１２に蓄積して保持する。本実施の形態では、このデータ記録部２１は、操作部１３から入力される指示に従い、受け入れたデータを、ノイズのみのデータ、ノイズと自発脳磁場とを含んだデータ、ノイズ、自発脳磁場、及び誘発脳磁場を含んだデータのいずれかに分類して記憶部１２に格納する。

【0025】

モデル学習部２２は、データ記録部２１によって記憶部１２に蓄積されたデータを参照し、過去に得られたデータ（被験者のいない状態で得られたノイズのみのデータ）を用いて定められる基準に基づいて、（６）式のモデルパラメータθベクトルの成分（各ノイズパラメータ）を決定する。

【0026】

具体的に、このモデル学習部２２は、過去に得られたノイズのみの時系列のデータｍ₁，ｍ₂，…，ｍ_tを記憶部１２から読み出す。そしてベクトルθが与えられたときにデータがこの順に観測される確率ｐ（ｍ₁，ｍ₂，…，ｍ_t｜θ）を最大とするベクトルθを求める。すなわち、現実の観測変数ｍ_tを生成する真の確率モデルは未知であるものの、当該真の確率モデルを、近似的に（５）式で与えられるものとし、真の確率モデルと近似的なモデルとの擬距離であるＫＬ（カルバック・ライブラー）ダイバージェンス（ＫＬ情報量）を最小とする（５）式のパラメータ（ベクトルθの成分）を求めることとするのである。

【0027】

例えば上記ＫＬ情報量を最小とする赤池情報量基準を用いてベクトルθを決定すればよい。このためモデル学習部２２は、上記確率ｐ（ｍ₁，ｍ₂，…，ｍ_t｜θ）の対数log（ｐ（ｍ₁，ｍ₂，…，ｍ_t｜θ））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｍ₁，ｍ₂，…，ｍ_t｜θ）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここではθの次元である、（ｌ＋４）となる。

【0028】

モデル学習部２２は、このＡＩＣの値を、複数の予め定めたθである値の組ごとに演算する。具体的にここで予め定めたθの値の組は、
ａ₁，ａ₂，…，ａ_l，σ_τ，σ_q，σ_e，σ_n
で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したＡＩＣの値が最小となるθを、上記複数の予め定めたθである値の組のうちから選択する。つまりモデル学習部２２は、ＡＩＣの値を最小とするベクトルθを、まずはグリッドサーチにより大域的に粗く探索する。

【0029】

モデル学習部２２は、さらにグリッドサーチにより求めたベクトルθ（θ₀とする）を初期値として、非線形最適化により最適化されたベクトル

【数7】

を求める。

【0030】

ここで非線形最適化の一例としては、上記対数尤度log（ｐ（ｍ₁，ｍ₂，…，ｍ_t｜θ））をカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθ₀を初期値として、対数尤度を最大とするベクトルθを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθを、

【数7】

とする。ここで終了条件は、更新前のθと更新後のθとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。

【0031】

ノイズパラメータ管理部２３は、モデル学習部２２が演算により求めた

【数7】

を、ノイズパラメータとして記憶部１２に格納する。また、このノイズパラメータ管理部２３は、信号推定部２４からノイズパラメータの要求を受けて、記憶部１２に格納したノイズパラメータを読み出し、信号推定部２４に対して出力する。

【0032】

本実施の形態の信号推定部２４は、自発脳磁場推定部２４ａと、誘発脳磁場推定部２４ｂとを含んで構成される。自発脳磁場推定部２４ａは、記憶部１２に格納されている、ノイズと自発脳磁場とを含んだ安静時データ（安静状態にある被験者を測定した脳磁場のデータ）を読み出す。ここで自発脳磁場推定部２４ａは、自発脳磁場が

【数8】

によって表される確率モデルで近似的に表されるものとする。ここにｒ_tは時刻ｔの自発脳磁場成分を意味し、δ波(約２Ｈｚ)，θ波(約４Ｈｚ)，α波(約１０Ｈｚ)，β波(約２０Ｈｚ)，low-γ波(約４０Ｈｚ)，high-γ波(約８０Ｈｚ)の成分から構成されるものとし、各成分を、

【数9】

と表したものである。なお、ここでのθは、θ波に対応する自発脳磁場成分であることを意味するものとして便宜的に利用しているもので、パラメータのベクトルθとは異なる。

【0033】

また、これらの成分の各々はいずれも、次の準周期振動モデルで表現されるものとする。

【0034】

【数10】

ここでωは、δ，θ，α，β，γ_l，γ_hのそれぞれをとる。また、ｆ_ωは、成分ω（δ，θ，α，β，γ_l，γ_hのそれぞれをとる）の中心周波数であり、ｖ_t^(ω)は、平均が０、分散がσ_ω²であるガウス分布に属するとの条件、

【数11】

を満足する確率的入力である。

【0035】

以上より、自発脳磁場の成分を表すモデルパラメータであるθベクトルは、

【数12】

となる。ここにおいて、左辺のθ_pのθはθベクトルを表すθであり、右辺中、σ_θ²とあるθは、θ波に対応する自発脳磁場成分を表す便宜的なものである。

【0036】

自発脳磁場推定部２４ａは、記憶部１２に格納されているノイズと自発脳磁場とを含んだ時系列の安静時データｍ₁，ｍ₂，…，ｍ_tを記憶部１２から読み出す。

【0037】

また自発脳磁場推定部２４ａは、ノイズパラメータ管理部２３に対してノイズパラメータを要求し、この要求に応答してノイズパラメータ管理部２３が出力するノイズパラメータと（５）式とを用いて、ノイズの時系列推定データν₁，ν₂，…を得る。自発脳磁場推定部２４ａは、ノイズと自発脳磁場とを含んだ時系列の安静時データｍ₁，ｍ₂，…，ｍ_tのそれぞれ対応する時刻のデータからノイズの時系列推定データの対応する時刻のデータを差引きして、ｍ′_i＝ｍ_i−ν_iを求め、ノイズがないと仮定される状態のデータ、ｍ′₁，ｍ′₂，…，ｍ′_tを得る。

【0038】

そして自発脳磁場推定部２４ａは、（１０）式で表されるベクトルθ_pが与えられたときに安静時データが実際に測定された時系列順に観測される確率ｐ（ｍ′₁，ｍ′₂，…，ｍ′_t｜θ_p）を最大とするベクトルθ_pを求める。この場合も、現実の観測変数ｍ′_tを生成する真の確率モデルは未知であるものの、当該真の確率モデルを、近似的に（７）式で与えられるものとし、真の確率モデルと近似的なモデルとの擬距離であるＫＬ（カルバック・ライブラー）ダイバージェンス（ＫＬ情報量）を最小とする（７）式のパラメータ（（１０）式のベクトルθ_pの成分）を求めることとするのである。

【0039】

例えば上記ＫＬ情報量を最小とする赤池情報量基準を用いて（１０）式のベクトルθを決定すればよい。このため自発脳磁場推定部２４ａは、上記確率ｐ（ｍ′₁，ｍ′₂，…，ｍ′_t｜θ_p）の対数log（ｐ（ｍ′₁，ｍ′₂，…，ｍ′_t｜θ_p））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｍ′₁，ｍ′₂，…，ｍ′_t｜θ_p）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここでは（１０）式のベクトルθ_pの次元である、６となる。

【0040】

自発脳磁場推定部２４ａは、このＡＩＣの値を、複数の予め定めたθ_pである値の組ごとに演算する。具体的にここで予め定めたθ_pの値の組は、σ_ω²（ただしωはδ，θ，α，β，γ_l，γ_hのそれぞれをとる）で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したＡＩＣの値が最小となるθの値の組を、上記複数の予め定めたθ_pとしての値の組のうちから選択する。つまり自発脳磁場推定部２４ａは、ＡＩＣの値を最小とするベクトルθ_pを、まずはグリッドサーチにより大域的に粗く探索する。

【0041】

自発脳磁場推定部２４ａは、さらにグリッドサーチにより求めたベクトルθ_p（θ_p0とする）を初期値として、非線形最適化により最適化された自発脳磁場のモデルのためのパラメータのベクトル

【数13】

を求める。

【0042】

ここで非線形最適化の一例としては、上記対数尤度log（ｐ（ｍ′₁，ｍ′₂，…，ｍ′_t｜θ_p））をカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθ_p0を初期値として、対数尤度を最大とするベクトルθ_pを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθ_pを、

【数13】

とする。ここで終了条件は、更新前のθ_pと更新後のθ_pとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。

【0043】

自発脳磁場推定部２４ａは、こうして求めた

【数13】

を、誘発脳磁場推定部２４ｂに出力する。

【0044】

誘発脳磁場推定部２４ｂは、記憶部１２に格納されている、ノイズと自発脳磁場とを含んだ刺激時データ（刺激を受けた状態にある被験者を測定した脳磁場のデータ）を読み出す。ここで誘発脳磁場推定部２４ｂは、誘発脳磁場を次の２次の線型モデル

【数14】

のインパルス応答で近似できるものとする。

【0045】

本実施の形態の誘発脳磁場推定部２４ｂは、誘発脳磁場を、確率的入力を持つ２次の定常線形モデルとして（１１）式で表現しておく。この（１１）式においてモデルパラメータはシステム入力強度と、システム入力に対するモデルの応答を決めるシステムパラメータｂ₁，ｂ₂である。

【0046】

ここでは単一の線型モデルのインパルス応答で誘発脳磁場を近似できるものとして、この刺激（システム入力）の確率分布を次のようにして与える。すなわち、誘発脳磁場推定部２４ｂは、システムパラメータとしてｂ₁，ｂ₂をもつシステムのインパルス応答ｈの時系列ｈ₁，ｈ₂，…と、得られている刺激時データｍ₁，ｍ₂，…との相互相関関数を求める。そして誘発脳磁場推定部２４ｂは、刺激の確率分布が、上記相互相関の絶対値が最大となる時間的ラグｔ（誘発脳磁場が現れる時刻ｔ）を中心とする二項分布
Ｂ（２ｔ，０．５）
であるとする。ここでは時間が離散化されているので二項分布とした。また、入力強度はラグｔでの相互相関関数の値を平均値Iとするガウス分布
Ｎ（Ｉ，σ_I²）
としておく。

【0047】

誘発脳磁場推定部２４ｂは、ノイズパラメータ管理部２３に対してノイズパラメータを要求し、この要求に応答してノイズパラメータ管理部２３が出力するノイズパラメータと（５）式とを用いて、ノイズの時系列推定データν₁，ν₂，…を得る。誘発脳磁場推定部２４ｂは、また自発脳磁場推定部２４ａから入力される自発脳磁場に関するパラメータと（７）式とを用いて、自発脳磁場の時系列推定データｇ₁，ｇ₂，…を得る。そして誘発脳磁場推定部２４ｂは、刺激時データｍ₁，ｍ₂，…，ｍ_tのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、ｍ″_i＝ｍ_i−ν_i−ｇ_iを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、ｍ″₁，ｍ″₂，…，ｍ″_tを得る。

【0048】

（１１）式では、モデルから直接にデータが得られるものではないので、上述のように刺激の確率分布を、ガウス分布を用いた演算式（これには上記パラメータｂ₁，ｂ₂が含まれる）で表現しておき、この演算式における時刻ｔでの期待値μ（ｔ）と分散σ²とを用いた確率差分方程式により、刺激に対する推定される誘発脳磁場の時刻ｔでの信号ｓ_t及び時刻ｔ−１での信号ｓ_t-1を、その成分に含む状態ベクトルｘの時間発展の演算式を得ておく。

【0049】

この演算式は、時刻ｔ，ｔ−１において推定される誘発脳磁場ｓ_tから、時刻ｔ＋１において推定される誘発脳磁場ｓ_t+1が得られる時間発展の確率（運動のモデル）ｐ（ｓ_t+1｜ｓ_t，ｓ_t-1）を表す。

【0050】

誘発脳磁場推定部２４ｂは、システムパラメータｂ₁，ｂ₂を用いて表される信号ｓ_t，ｓ_t-1を成分として含む状態ベクトルｘ_tについて、初期状態ｐ（ｘ₀）（これは誘発脳磁場が発生していないものとして既知である）と、上記時間発展の確率とを用いて、時刻ｔで観測されるデータがｍ″_tであったとき、状態ベクトルがｘ_tである確率ｐ（ｘ_t｜ｍ″_t）を演算する。本実施の形態の誘発脳磁場推定部２４ｂは、この確率ｐ（ｘ_t｜ｍ″_t）を、次のように求める。

【0051】

すなわち、ベイズの定理より、

【数15】

が成立するので、これを変形して、

【数16】

を得る。ここでｐ（ｍ″_t）は、

【数17】

と表すことができ、またｐ（ｘ_t）は、時間発展の確率を用いて、

【数18】

と表すことができるから、結局（１３）式は、

【数19】

とすることができる。

【0052】

この（１６）式は、時刻ｔで観測されるデータがｍ″_tであったとき、状態ベクトルがｘ_tである確率密度ｐ（ｘ_t｜ｍ″_t）が、観測のモデルｐ（ｍ″_t｜ｘ_t）と、運動のモデルｐ（ｘ_t+1｜ｘ_t）と、初期状態ｐ（ｘ₀）とから推定できることを表す。

【0053】

誘発脳磁場推定部２４ｂは、システムパラメータｂ₁，ｂ₂を含むベクトルθsが与えられたときに確率ｐ（ｘ_t｜ｍ″_t）を最大にするベクトルθsを求める。具体的に誘発脳磁場推定部２４ｂは、ＫＬ情報量を最小とする赤池情報量基準を用いてベクトルθsを決定する。このため誘発脳磁場推定部２４ｂは、上記確率ｐ（ｘ_t｜ｍ″_t）の対数log（ｐ（ｘ_t｜ｍ″_t））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｘ_t｜ｍ″_t）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここではθsの次元である、２となる。

【0054】

誘発脳磁場推定部２４ｂは、このＡＩＣの値を、複数の予め定めたθsである値の組ごとに演算する。具体的にここで予め定めたθsの値の組は、
ｂ₁，ｂ₂
で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したＡＩＣの値が最小となるθsを、上記複数の予め定めたθsとしての値の組のうちから選択する。つまり誘発脳磁場推定部２４ｂは、ＡＩＣの値を最小とするベクトルθsを、まずはグリッドサーチにより大域的に粗く探索する。

【0055】

誘発脳磁場推定部２４ｂは、さらにグリッドサーチにより求めたベクトルθs（θs₀とする）を初期値として、非線形最適化により最適化されたベクトル

【数20】

を求める。

【0056】

ここで非線形最適化の一例としては、上記対数尤度log（ｐ（ｘ_t｜ｍ″_t））をカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθs₀を初期値として、対数尤度を最大とするベクトルθsを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθsを、

【数20】

とする。ここで終了条件は、更新前のθsと更新後のθsとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。

【0057】

信号出力部２５は、誘発脳磁場推定部２４ｂが演算により得た

【数20】

と、（１１）式で表されるモデルの式とを用いて、誘発脳磁場の時間変化波形を演算して出力する。

【0058】

さらに本実施の形態の制御部１１は、以上の演算を、複数のセンサ（脳磁計の場合はSQUIDs）から得られた各データについて行い、センサごとに

【数20】

のパラメータを得てもよい。このときには、複数のセンサの一つを注目センサとするとき、注目センサに係るパラメータを得るためには、ノイズのみのデータや安静時データ、刺激時データは、当該注目センサによって得られたものを用いてもよい。

【0059】

そして制御部１１は、複数のセンサの各々に対応して得られたパラメータと（１１）式とを用い、各センサの位置における誘発脳磁場の時間変化波形を演算し、等磁界線図を作成して出力してもよい。この等磁界線図の作成処理は広く知られているので、ここでの詳細な説明を省略する。

【0060】

［モンテカルロ法による演算］
ところで制御部１１における（１６）式の演算では、右辺分母の積分（ベイズ推定における、事前確率と事後確率との積を被積分関数とする積分）が解析的な解が得られないので困難になる。そこで本実施の形態の制御部１１は、複数の演算素子を含み、次のように演算を行うこととしてもよい。

【0061】

ここで複数の演算素子は、ＣＰＵでなく、ＧＰＵ（Graphics Processing Unit）等の演算素子を用いてもよい。このように複数の演算素子を利用できる場合、制御部１１は、各演算素子に対し、積分範囲内からランダムにｘ_tとｘ_t-1との値を選択させ、当該ランダムに選択した積分変数の値（複数ある場合はそれぞれ積分範囲からランダムに選択する）を被積分関数に代入して得た被積分関数値を、並列的に演算させる。そして制御部１１は、各演算素子で得られた複数の被積分関数値を累算してモンテカルロ法により積分を遂行する。

【0062】

［ノイズのみのデータを得るタイミング］
また、以上の説明において、ノイズパラメータの決定に用いた、ノイズのみのデータは、演算を行う時点より過去であれば、いつ得られたものであってもよかった。例えば刺激データを得る直前（刺激データを実際に測定する被験者が入る直前に、測定に先立って得られたデータ）であってもよいし、刺激データを測定した後、被験者が退席してから得られたデータであってもよい。

【0063】

［自発脳磁場をベイズ推定で求める例］
ここまでの説明において、自発脳磁場推定部２４ａは、（１０）式で表されるベクトルθが与えられたときに、安静時データが、測定された時系列順に観測される確率ｐ（ｍ′₁，ｍ′₂，…，ｍ′_t｜θ）を最大とするベクトルθを求めていた。しかしながら、本実施の形態はこれに限られない。

【0064】

すなわち、自発脳磁場推定部２４ａは、時刻ｔでの状態ベクトルｘ_tとして、時刻ｔ以前の自発脳磁場を表すデータｇ_t，ｇ_t-1，…と、時刻ｔ以前の各ノイズを表すデータτ_t，ｑ_t,ｑ_t-1, …, ｑ_t-T+2, ｅ_t, …, ｅ_t-l+1，ｎ_t…とを含めた
x_t=[τ_t，ｑ_t,ｑ_t-1, …, ｑ_t-T+2, ｅ_t, …, ｅ_t-l+1，ｎ_t，ｇ_t，ｇ_t-1，…]^T
を生成する。なお、肩のＴは、転置を表す。

【0065】

ここで、各ノイズのデータや自発脳磁場のデータを、時刻ｔよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻ｔ＋１における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なり、例えば自発脳磁場でいえば、ｇ_t一つだけの場合もあり得る。

【0066】

自発脳磁場推定部２４ａは、誘発脳磁場推定部２４ｂが行う（１６）式の演算と同様の演算により、（１０）式のベクトルθが与えられたときに確率ｐ（ｘ_t｜ｍ′_t）を最大にするベクトルθを求める。具体的に自発脳磁場推定部２４ａは、ＫＬ情報量を最小とする赤池情報量基準を用いてベクトルθを決定する。このため自発脳磁場推定部２４ａは、上記確率ｐ（ｘ_t｜ｍ′_t）の対数log（ｐ（ｘ_t｜ｍ′_t））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｘ_t｜ｍ′_t）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここではθの次元である、６となる。

【0067】

そして自発脳磁場推定部２４ａは、このＡＩＣを最小とするθを、既に述べたグリッドサーチや非線形最適化を組み合わせた方法で見出すこととしてもよい。

【0068】

この際の（１６）式の演算においても、モンテカルロ法による演算を行うことができる。

【0069】

［自発脳磁場もベイズ推定に含める例］
さらに、ここまでの説明では、誘発脳磁場推定部２４ｂは、自発脳磁場推定部２４ａから入力される自発脳磁場に関するパラメータと（７）式とを用いて、自発脳磁場の時系列推定データｇ₁，ｇ₂，…を得て、刺激時データｍ₁，ｍ₂，…，ｍ_tのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、ｍ″_i＝ｍ_i−ν_i−ｇ_iを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、ｍ″₁，ｍ″₂，…，ｍ″_tを得ていた。

【0070】

しかしながら本実施の形態はこれに限られない。すなわち誘発脳磁場推定部２４ｂは、刺激時データｍ₁，ｍ₂，…，ｍ_tのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータを差引きして、ｍ″_i＝ｍ_i−ν_iを求め、ノイズがないと仮定される状態のデータ、ｍ″₁，ｍ″₂，…，ｍ″_tを得ておき、次のように処理を遂行してもよい。

【0071】

この場合、誘発脳磁場推定部２４ｂは、時刻ｔでの状態ベクトルｘ_tに、時刻ｔ以前の自発脳磁場を表すデータｇ_t，ｇ_t-1，…を含めて、
x_t=[s_t,s_t-1，ｇ_t，ｇ_t-1，…]^T
とする。なお、肩のＴは、転置を表す。

【0072】

ここで、自発脳磁場を時刻ｔよりどれだけ前のデータから含めるかは、ｇ_t+1を予測するに十分な数とする。これは採用するモデルによっても異なり、例えばｇ_t一つだけの場合もあり得る。

【0073】

誘発脳磁場推定部２４ｂが行う（１６）式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部２４ｂは、既に述べたと同様にして

【数20】

を求めることができる。

【0074】

［ノイズもベイズ推定に含める例］
さらに、ここまでの説明では、誘発脳磁場推定部２４ｂは、自発脳磁場の時系列推定データｇ₁，ｇ₂，…を得て、刺激時データｍ₁，ｍ₂，…，ｍ_tのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、ｍ″_i＝ｍ_i−ν_i−ｇ_iを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、ｍ″₁，ｍ″₂，…，ｍ″_tを得ていた。

【0075】

しかしながら本実施の形態はこれに限られない。すなわち誘発脳磁場推定部２４ｂは、刺激時データｍ₁，ｍ₂，…，ｍ_tのそれぞれ対応する時刻のデータから、自発脳磁場の時系列推定データｇ₁，ｇ₂，…の対応する時刻のデータを差引きして、ｍ″_i＝ｍ_i−ｇ_iを求め、自発脳磁場がないと仮定される状態のデータ、ｍ″₁，ｍ″₂，…，ｍ″_tを得ておき、次のように処理を遂行してもよい。

【0076】

この場合、誘発脳磁場推定部２４ｂは、時刻ｔでの状態ベクトルｘ_tに、時刻ｔ以前の各ノイズを表すデータτ_t，ｑ_t,ｑ_t-1, …, ｑ_t-T+2, ｅ_t, …, ｅ_t-l+1，ｎ_t…を含め、
x_t=[s_t, s_t-1，τ_t，ｑ_t, ｑ_t-1, …, ｑ_t-T+2,ｅ_t, …, ｅ_t-l+1，ｎ_t]^T
とする。なお、肩のＴは、転置を表す。

【0077】

ここで、各ノイズのデータを、時刻ｔよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻ｔ＋１における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なる。

【0078】

誘発脳磁場推定部２４ｂが行う（１６）式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部２４ｂは、既に述べたと同様にして

【数20】

を求めることができる。

【0079】

［ノイズも自発脳磁場もベイズ推定に含める例］
さらに誘発脳磁場推定部２４ｂは、刺激時データｍ₁，ｍ₂，…，ｍ_tのそれぞれを、そのまま上記ｍ″₁，ｍ″₂，…，ｍ″_tとして、次のように処理を遂行してもよい。

【0080】

すなわち誘発脳磁場推定部２４ｂは、時刻ｔでの状態ベクトルｘ_tに、時刻ｔ以前の各ノイズを表すデータτ_t，q_t,ｑ_t-1, …, ｑ_t-T+2, ｅ_t, …, ｅ_t-l+1，ｎ_t…を含め、
x_t=[s_t,s_t-1，τ_t，ｑ_t, ｑ_t-1, …, ｑ_t-T+2,ｅ_t, …, ｅ_t-l+1，ｎ_t，ｇ_t，ｇ_t-1，…]^T
とする。なお、肩のＴは、転置を表す。

【0081】

ここで、各ノイズのデータを、時刻ｔよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻ｔ＋１における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なる。

【0082】

また誘発脳磁場推定部２４ｂは、時刻ｔ以前の自発脳磁場を表すデータｇ_t，ｇ_t-1，…も時刻ｔでの状態ベクトルｘ_tに含める。ここで、自発脳磁場を時刻ｔよりどれだけ前のデータから含めるかは、ｇ_t+1を予測するに十分な数とする。これは採用するモデルによっても異なり、例えばｇ_t一つだけの場合もあり得る。

【0083】

誘発脳磁場推定部２４ｂが行う（１６）式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部２４ｂは、既に述べたと同様にして

【数20】

を求めることができる。

【0084】

［グリッドサーチにおけるＡＲモデルの変数変換］
さらにここまでの説明では定常ノイズ成分ｅ_tの時間発展を表すＡＲモデルにおいては、安定性の条件が見にくいので、次の方法により繰り返し変数変換を行って、１次のＡＲモデルに変換してもよい。

【0085】

すなわち、Ｌ次のＡＲモデルをＬ−１次のＡＲモデルに変換する方法として、

【数21】

を用いる。ただし、ｉ＝１，２，…，ｌ−１である。ここにａの肩にある添字のｌは、ｌ次の変数であることを意味する。この（１７）式による変換を、１次となるまで繰り返して行うと、各変換の段階で、ａ_j^jが得られる（ここにｊ＝１，２，…ｌ）。すると、安定なＡＲモデルとなるための必要十分条件が、

【数22】

と得られる。そこで、グリッドサーチにおいては、この（１８）式で表される範囲を定義域として、この定義域を均等に分割した格子点を用いればよい。

【0086】

［他のモデル］
なお、一般に活動部位や刺激の種類に応じて、誘発脳磁場モデルの入力の確率分布は異なると考えられるうえ、線形モデルに限らず非線形なモデルが適当な場合が想定される。また本実施の形態は脳磁計のデータに限らず、一般的な測定データの処理に用いることができるものである。そこで利用者がそれぞれの仮説に従い、測定対象となるデータのモデル（ここでの例では誘発脳磁場のモデル）の確率分布を作成してもよい。本実施の形態の処理は、種々の確率分布を仮定した場合でも実行可能なものである。

【0087】

また、自発脳磁場のモデルについても、ここでは（８）式によって表されるものとしたが、振幅の時間変化が大きい場合には、確率的入力項のパワーが時間的に変化するものと仮定し、当該パワーをパラメータを含む何らかの関数形で記述し、当該パラメータをモデルパラメータに含めて処理を行ってもよい。さらに、被験者の各自発脳磁場帯域のピーク周波数が典型的なものと大きく異なる場合は、中心周波数もモデルパラメータに含めてもよい。これらの場合も、上述の処理と同様にモデルパラメータの一つとして推定を行うことができる。

【0088】

［動作］
本実施の形態は、以上の構成を有し、次のように動作する。すなわち、本実施の形態のデータ処理装置１は、ノイズのみのデータ、安静時データ、刺激時データを受け入れて記憶部１２に保持しており、図３に例示するように、まずノイズのみのデータを用いて予め定めたノイズのモデルに含まれるノイズパラメータを推定する（Ｓ１）。

【0089】

またデータ処理装置１は、安静時データに含まれるノイズを、上記推定したノイズパラメータと予め定めたノイズのモデルとを用いて除去し、ノイズを除去した安静時データを参照して、予め定めた安静時データのモデルに含まれるパラメータを推定する（Ｓ２）。

【0090】

そしてデータ処理装置１は、刺激時データから、処理Ｓ１で推定したノイズパラメータと予め定めたノイズのモデルとを用いて除去し、刺激時データに含まれる自発脳磁場を、処理Ｓ２で推定したパラメータと、予め定めた安静時データのモデルとを用いて除去する。そしてデータ処理装置１は、こうしてノイズ及び自発脳磁場を除去したデータを参照して、予め定めた刺激時データのモデル（誘発脳磁場のモデル）に含まれるパラメータを推定する（Ｓ３）。

【0091】

データ処理装置１は、処理Ｓ３で推定したパラメータと誘発脳磁場のモデルとを用いて、誘発脳磁場のデータを推定して出力する（Ｓ４）。

【実施例】

【0092】

本発明の一実施例について説明する。ここでは被験者に、「ボタンを押してもらう」という刺激を与えたときの誘発脳磁場を得たときの例について述べる。

【0093】

図４は、ＬＰＦ（Low Pass Filter）を経た１６９回の測定データを加算平均した、従来の方法で得られた誘発脳磁場の波形（複数のセンサの信号を重ね合わせて表している、以下も同様）と、これらの波形から得られた等磁界線図である。図５は、ノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図である。この図５のようなデータを１６９回得て加算平均したものが図４に示した従来の結果である。

【0094】

図６は、この図５に示したデータ（単一測定のデータ）に基づき、上記本実施の形態で説明した処理を経て本発明の実施例に係るデータ処理装置１が出力する誘発脳磁場の波形（複数のセンサの信号を重ね合わせて表している）と、これらの波形から得られた等磁界線図である。

【0095】

図４に示した従来の解析図と、図６に示した本実施例の解析図とを比較すると、略一致していることが理解される。

【0096】

図７は、自発脳磁場の強度が大きい場合のノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図である。この図７に例示したようなデータを複数回得て加算平均した、従来の方法によって推定された等磁界線図と、この図７に例示したような複数のデータのそれぞれから、本実施例のデータ処理装置１で推定して得た複数の等磁界線図（図８に例示する）とのコヒーレンスを演算して表したものを図９に示す。また、従来の方法で推定された等磁界線図を図１０に示す。図９，図１０から見られるように、図１０で得られた等磁界線図のうち、信号の強度が大きい領域でのコヒーレンスは、図９に示されるように０．８を超えている。つまり、本実施例のデータ処理装置１での推定は、ロバストネスの面からも十分であることが理解される。

【符号の説明】

【0097】

１ データ処理装置、１１ 制御部、１２ 記憶部、１３ 操作部、１４ 出力部、１５ 入力インタフェース、２１ データ記録部、２２ モデル学習部、２３ ノイズパラメータ管理部、２４ 信号推定部、２５ 信号出力部。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】