JP7585655B2 - 波形発生特定方法およびプログラム - Google Patents

Description

本発明は、波形発生特定方法およびプログラムに関する。

従来、脳磁計や脳波計の解析において、特徴的な波形を解析することは非常に重要である。ただし、技術の進歩によって計測時のサンプリング周波数やセンサの数は増加傾向にあり、結果として目視での波形の探索にかかる時間も増加する傾向がある。

例えば、臨床における脳磁計を用いたてんかんの診断では、等価電流双極子法と呼ばれる手法を用いててんかん病変部の局在性を評価する。等価電流双極子法では、頭皮上で計測された磁場を生み出す電流源（ダイポール）を推定する。ダイポールを推定するためには、複数センサの時系列の中から、特徴的な波形情報（ＩＥＤ：Interictal Epileptiform Discharge）が発生した時刻（起始部）とその波形情報が出現しているセンサの絞り込みが必要となる。

現状では、医師は手動でＩＥＤの探索・起始部の決定を行っているが、脳磁計のデータは膨大であるため、個別のＩＥＤについて手動で正確にセンサおよびＩＥＤの時刻の抽出を行うのは困難である。

そこで、特許文献１には、精度良く特徴波を検出する目的で、単一センサもしくは複数センサにおける脳電磁波の時系列データをウェーブレット変換して、ＩＥＤをはじめとした特徴的な波形を検出する技術が開示されている。

しかしながら、特許文献１に開示の技術によれば、等価電流双極子法を実行するためのＩＥＤの時刻の決定やセンサの絞り込みが適切に行えないという問題があった。

本発明は、上記に鑑みてなされたものであって、特徴的な波形情報（ＩＥＤ）が現れている時刻の決定とセンサの抽出をより正確に行うことを目的とする。

上述した課題を解決し、目的を達成するために、本発明の波形発生特定方法は、複数のセンサにより取得した個々の波形データについて、少なくとも１つ以上の特徴的な波形情報と比較する比較ステップと、前記波形データの区間と前記特徴的な波形情報との相関の割合に基づいて、前記波形データの少なくとも一定区間における特徴的な波形情報の出現確度を決定する決定ステップと、前記決定された特徴的な波形情報の出現確度に対し、閾値を用いてセンサの領域を抽出し、該領域に対してフィルタを適用して特徴的な波形情報の時刻と前記センサの領域とを絞り込む絞込ステップと、前記絞込ステップで絞り込まれた時刻及びセンサの領域を用いて、特徴的な波形情報に一致する区間の出現した時刻及びセンサのインデックスを特定する特定ステップと、を含むことを特徴とする。

本発明によれば、特徴的な波形情報（ＩＥＤ）が現れている時刻の決定とセンサの抽出をより正確に行うことができる。

図１は、実施形態に係る生体信号計測システムの構成を示す概略図である。 図２は、サーバの機能構成を概略的に示すブロック図である。 図３は、情報処理装置のハードウェア構成の一例を示す図である。 図４は、情報処理装置の機能ブロック構成の一例を示す図である。 図５は、ダイポール推定処理の流れを示すフローチャートである。 図６は、機械学習を用いたＩＥＤの確率マップの生成の際の学習処理の流れを示すフローチャートである。 図７は、ＩＥＤの確率マップの一例を示す図である。 図８は、閾値処理後のＩＥＤ確率マップの一例を示す図である。 図９は、後処理後のＩＥＤ確率マップの一例を示す図である。 図１０は、ＩＥＤ確率マップを用いたＵＩ例を示す図である。 図１１は、機械学習以外の手法を用いたＩＥＤの確率マップの算出例を示す図である。 図１２は、センサ数の拡張方法の一例を示す概念図である。

以下に添付図面を参照して、波形発生特定方法およびプログラムの実施の形態を詳細に説明する。また、以下の実施形態によって本発明が限定されるものではなく、以下の実施形態における構成要素には、当業者が容易に想到できるもの、実質的に同一のもの、およびいわゆる均等の範囲のものが含まれる。さらに、以下の実施の形態の要旨を逸脱しない範囲で構成要素の種々の省略、置換、変更および組み合わせを行うことができる。

本実施形態は、ＩＥＤ（Interictal Epileptiform Discharge）の確率マップを算出するに際して、以下の特徴を有する。要するに、ＩＥＤの確率マップを用いることで、等価電流双極子法に必要なＩＥＤの起始時刻の検出とセンサの選択が可能になることが特徴になっている。

医師が解析を行う際には、ＩＥＤの起始部を等価電流双極子法で解析することが理想ではあるが、起始部はＳ／Ｎが低いことが知られているため、実際には起始部からＩＥＤのピークまでの中でＩＥＤの時刻を決定している。本実施形態では起始部に限定しないで解析を行う。

（生体信号計測システムの概略）
図１は、実施形態に係る生体信号計測システムの構成を示す概略図である。図１を参照しながら、本実施形態に係る生体信号計測システム１の概略について説明する。

生体信号計測システム１は、特定の発信源（生体部位）からの被検者の複数種類の生体信号（例えば、脳磁（ＭＥＧ：Magneto－encephalography）信号、および脳波（ＥＥＧ：Electro－encephalography）信号を計測し、表示するシステムである。なお、本発明において、測定対象となる生体信号は、脳磁信号および脳波信号に限られるものではなく、例えば、心臓の活動に応じて発生する電気信号（心電図として表現可能な電気信号）などであってもよい。

図１に示すように、生体信号計測システム１は、被検者の１以上の生体信号を測定する測定装置３と、測定装置３で測定された１種類以上の生体信号を記録するサーバ４０と、サーバ４０に記録された１種類以上の生体信号を解析する生体信号表示装置である情報処理装置５０と、を含む。測定装置３は、例えば脳磁場や刺激を与えたタイミング等の計測値を収集する脳磁計である。なお、図１では、サーバ４０と情報処理装置５０とが別々に記載されているが、例えば、サーバ４０が有する機能の少なくとも一部が情報処理装置５０に組み込まれる形態であってもよい。

図１の例では、被検者（被測定者）は、頭に脳波測定用の電極（またはセンサ）を付けた状態で測定テーブル４に仰向けで横たわり、測定装置３のデュワ３１の窪み３２に頭部を入れる。デュワ３１は、液体ヘリウムを用いた極低温環境の保持容器であり、デュワ３１の窪み３２の内側には脳磁測定用の多数の磁気センサが配置されている。測定装置３は、電極からの脳波信号と、磁気センサからの脳磁信号とを収集し、収集した脳波信号および脳磁信号を含むデータ（以下、「測定データ」と称する場合がある）をサーバ４０へ出力する。サーバ４０へ出力された測定データは、情報処理装置５０に読み出されて表示され、解析される。一般的に、磁気センサを内蔵するデュワ３１および測定テーブル４は、磁気シールドルーム内に配置されているが、図１では便宜上、磁気シールドルームの図示を省略している。

情報処理装置５０は、複数の磁気センサからの脳磁信号の波形データと、複数の電極からの脳波信号の波形データを、同じ時間軸上に同期させて表示する装置である。脳波信号とは、神経細胞の電気的な活動（シナプス伝達の際のニューロンの樹状突起で起きるイオン電荷の流れ）を電極間の電圧値として表される信号である。脳磁信号とは、脳の電気活動により生じた微小な電場変動を表す信号である。脳磁場は、高感度の超伝導量子干渉計（ＳＱＵＩＤ：Superconducting Quantum Interference Device）センサで検知される。これらの脳波信号および脳磁信号は、「生体信号」の一例である。

ここで、図２はサーバ４０の機能構成を概略的に示すブロック図である。図２に示すように、サーバ４０は、データ取得部４１と、データ記憶部４２とを有する。

データ取得部４１は、測定装置３から測定データを定期的に取得する。ここで、測定データは、測定装置３のデュワ３１の複数の磁気センサで計測した個々の波形データである。

データ記憶部４２は、測定装置３から取得した測定データを記憶する。

（情報処理装置のハードウェア構成）
図３は、情報処理装置５０のハードウェア構成の一例を示す図である。図３を参照しながら、情報処理装置５０のハードウェア構成について説明する。

図３に示すように、情報処理装置５０は、ＣＰＵ（Central Processing Unit）１０１と、ＲＡＭ（Random Access Memory）１０２と、ＲＯＭ（Read Only Memory）１０３と、補助記憶装置１０４と、ネットワークＩ／Ｆ１０５と、入力装置１０６と、表示装置１０７と、を有し、これらがバス１０８で相互に接続されている。

ＣＰＵ１０１は、情報処理装置５０の全体の動作を制御し、各種の情報処理を行う演算装置である。ＣＰＵ１０１は、ＲＯＭ１０３または補助記憶装置１０４に記憶された情報表示プログラムを実行して、測定収集画面および解析画面の表示動作を制御する。

ＲＡＭ１０２は、ＣＰＵ１０１のワークエリアとして用いられ、主要な制御パラメータおよび情報を記憶する揮発性の記憶装置である。ＲＯＭ１０３は、基本入出力プログラム等を記憶する不揮発性の記憶装置である。例えば、上述の情報表示プログラムがＲＯＭ１０３に記憶されているものとしてもよい。

補助記憶装置１０４は、ＨＤＤ（Hard Disk Drive）またはＳＳＤ（Soliｄ State Drive）等の記憶装置である。補助記憶装置１０４は、例えば、情報処理装置５０の動作を制御する制御プログラム、ならびに、情報処理装置５０の動作に必要な各種のデータおよびファイル等を記憶する。

ネットワークＩ／Ｆ１０５は、サーバ４０等のネットワーク上の機器と通信を行うための通信インターフェースである。ネットワークＩ／Ｆ１０５は、例えば、ＴＣＰ（Transmission Control Protocol）／ＩＰ（Internet Protocol）に準拠したＮＩＣ（Network Interface Card）等によって実現される。

入力装置１０６は、タッチパネルの入力機能、キーボード、マウスおよび操作ボタン等のユーザインターフェース等である。表示装置１０７は、各種の情報を表示するディスプレイ装置である。表示装置１０７は、例えば、タッチパネルの表示機能、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）または有機ＥＬ（Electro－Luminescence）等によって実現される。表示装置１０７は、測定収集画面および解析画面を表示し、入力装置１０６を介した入出力操作に応じて画面が更新される。

なお、図３に示す情報処理装置５０のハードウェア構成は一例であり、これ以外の装置が備えられるものとしてもよい。また、図３に示す情報処理装置５０は、例えば、ＰＣ（Personal Computer）を想定したハードウェア構成であるが、これに限定されるものではなく、タブレット等のモバイル端末であってもよい。この場合、ネットワークＩ／Ｆ１０５は、無線通信機能を有する通信インターフェースであればよい。

（情報処理装置の機能ブロック構成）
図４は、情報処理装置５０の機能ブロック構成の一例を示す図である。図４を参照しながら、情報処理装置５０の機能ブロック構成について説明する。

図４に示すように、情報処理装置５０は、前処理部５０１と、ＩＥＤ確率マップ算出部５０２と、閾値処理部５０３と、後処理部５０４と、ダイポール推定部５０５と、を備える。

前処理部５０１は、センサの抽出・拡張、ダウンサンプリングや周波数フィルタの適用、アーチファクト除外、不良チャネル処理、時間窓切り出し、磁場データの標準化などの前処理を行う。

センサの抽出・拡張は、後述のＩＥＤ確率マップ算出の際に、全センサを使って算出する方法以外にも、予めグルーピングした一定数センサ群のみを使用して算出することも可能である。センサのグルーピングには、側頭葉や前頭葉など解剖学的な基準に即して設定したものや、単純に近くにある任意の数のセンサをグルーピングしたものなどが考えられる。また、後述の学習過程において、使用したセンサ数よりも、算出の際のセンサ数が少ない場合には、架空のセンサを拡充することも可能である。

ダウンサンプリングは、学習時に使用したサンプリング周波数に合わせる目的で適用される。周波数フィルタも同様で学習時に適用したフィルタリングと同じものを適用する。よく使われるフィルタは、３５Ｈｚのローパスフィルタ、３Ｈｚ－３５Ｈｚのバンドパスフィルタなどである。

アーチファクト除外は、心拍性のアーチファクトや瞬き・体動によるアーチファクトを除外する目的で、ＩＣＡ（非特許文献１を参照）やＤＳＳＰ（非特許文献２を参照）などが適用される。

不良チャネル処理は、予め設定した閾値を超えるような磁場変化が観測されたセンサを除外したり、周囲のセンサ値を使用して補間したりする処理を指す。

時間窓の切り出し方は、時間窓の長さ分だけずらしていくオーバーラップなしのものや、時間窓の半分の長さをオーバーラップさせる方法、１／４の長さをオーバーラップさせる方法などがある。オーバーラップさせる場合には、後述のＩＥＤ確率マップ算出の際にオーバーラップ分の加算平均を行う。

磁場データの標準化は、切り出した時間窓の中で平均が０、分散が１になるような標準化が適用される。標準化の他に、予め設定した磁場の範囲内を－１～１となるように正規化する手法などを用いてもよい。

ＩＥＤ確率マップ算出部５０２は、特徴的な波形（ＩＥＤ：Interictal Epileptiform Discharge）の確率マップを算出する。

ＩＥＤとは、てんかん病変部に特徴的なSpike Wave、Spike and Wave、Poly Spike and Wave、Sharp Waveなどの波形を指す。本実施形態では、ＩＥＤ確率マップ算出部５０２は、事前に学習した機械学習のモデルを用いて確率マップを算出する。機械学習のモデルの作成方法については後述する。

閾値処理部５０３は、閾値用いて、ＩＥＤ確率マップ算出部５０２において得られたＩＥＤの確率マップからＩＥＤの確率が高い時間とセンサの領域の絞り込みを行う。

後処理部５０４は、閾値処理部５０３で閾値処理したＩＥＤの確率マップからダイポール推定に使用するためのサンプリング時点およびセンサの抽出を行う後処理を実行する。また、後処理部５０４は、ピーク検出手法やマップが複数のピークを持つ場合（病変部が２つ存在する場合などにみられるケース）には、マップを分離するための後処理を行う。

ダイポール推定部５０５は、ダイポール推定部５０５以前の処理で得られたセンサとＩＥＤの起始部のサンプリング時点を用いてダイポール推定（非特許文献３を参照）を行う。

次に、ダイポール推定処理の流れについて説明する。

ここで、図５はダイポール推定処理の流れを示すフローチャートである。処理に先立ち、測定装置３は、脳磁計計測を実施し、デュワ３１の複数の磁気センサで計測した個々の波形データをサーバ４０に出力する。サーバ４０は、測定装置３のデュワ３１の複数の磁気センサで計測した個々の波形データをデータ記憶部４２に記憶する。

そして、図５に示すように、まず、前処理部５０１は、測定装置３のデュワ３１の複数の磁気センサで計測した個々の波形データをサーバ４０のデータ記憶部４２から取得する（ステップＳ１）。

次に、前処理部５０１は、取得したデータに対し、センサの抽出・拡張、ダウンサンプリングや周波数フィルタの適用、アーチファクト除外、不良チャネル処理、時間窓切り出し、磁場データの標準化などの前処理を行う（ステップＳ２）。

次に、ＩＥＤ確率マップ算出部５０２は、特徴的な波形（ＩＥＤ）の確率マップを算出する（ステップＳ３）。

ＩＥＤの確率マップの算出については、機械学習を用いて算出したモデルを適用する方法のほかにも、従来のスパイク位置の検出アルゴリズム（非特許文献４を参照）を用いたものでもよい。機械学習を用いる場合には、医師が手動で検出したＩＥＤの時点を正解として学習させるため、従来のスパイク位置を検出する手法と比較して、医師が解析した場合と類似した結果が得られるという効果がある。従来の手法では、ＩＥＤの起始部もしくはＩＥＤのピークを検出することを目的として開発されているが、実際の解析では、医師は波形や所見などからダイポール推定を行う時刻を起始部からピークまでの間で調整しながら解析している。

ここで、ＩＥＤ確率マップ算出部５０２が、ＩＥＤの確率マップの生成において機械学習を用いる場合について説明する。

図６は、機械学習を用いたＩＥＤの確率マップの生成の際の学習処理の流れを示すフローチャートである。図６に示すように、ＩＥＤ確率マップ算出部５０２は、前処理部５０１で前処理を実行した前処理済み脳磁図データを取得すると（ステップＳ１１）、前処理済み脳磁図データに対して、ＩＥＤの確率マップを出力するための学習を行う（ステップＳ１２）。なお、学習の際に用いる正解データには、ダイポール推定を行った際に記録されているサンプリング時点とその時使用したセンサの情報を用いる。

具体的には、非特許文献５に記載のＵ－Ｎｅｔに代表されるSemantic Segmentationというタスクでよく使用されているネットワークを応用することで実現される。Semantic Segmentationは、入力となる全配列の各要素についてラベリングを行うタスクのことで、医療画像においては、脳腫瘍や癌の領域を自動で推定する際に用いられる。ここでは、医師がダイポール推定を行った際に選択されたセンサと推定を行った時刻から作成したマスクデータを正解として、入力の脳磁図データがマスクデータとなるように学習を行う。このとき、汎化性能を高める目的で、センサの順番をランダムに並び替えたり、特定のセンサのある範囲の値を０にしたりといったデータ拡張手法を適用してもよい。

また、Semantic Segmentation以外のタスクであるObject DetectionやInstance Segmentationといったタスクで使われるネットワークでも同様に実現可能である。Object Detectionは、入力となる配列中のどこに目的となるobjectが存在しているかを矩形などで推定する手法で、この手法を使用する場合には、ＩＥＤ確率マップは算出できないが、検出された矩形からダイポール推定を行う時刻（例えば矩形の中点）とセンサを直接取得することができる。Instance Segmentationは、Semantic SegmentationとObject Detectionを混ぜたようなタスクで、検出したObjectの分類と個数に加えてそのObjectの領域も算出することが可能である。Instance Segmentationでは、ＩＥＤの確率マップを算出することができるので、そのまま応用可能である。

次いで、ＩＥＤ確率マップ算出部５０２は、ステップＳ２の学習が終わると、学習済みモデルを作成する（ステップＳ１３）。作成された学習済みモデルは、機械学習の推論、つまり未知データに対してＩＥＤを検出する時に用いられる。

ここで、ＩＥＤ確率マップ算出部５０２が算出するＩＥＤの確率マップについて説明する。

図７は、ＩＥＤの確率マップの一例を示す図である。図７においては、横軸が時間、縦軸はセンサのインデックスを表している。図７に示す例の場合、時間のインデックスが１００、センサのインデックスが４０～８０あたりにＩＥＤが出現したと推定される。なお、図７においては、ＩＥＤの確率に応じて白黒の濃淡で表現したが、ＩＥＤの確率に応じて色が変化するようなカラー表示でもよい。

図５に戻り、次に、閾値処理部５０３は、ＩＥＤ確率マップ算出部５０２において得られたＩＥＤの確率マップに対して閾値処理を行う（ステップＳ４）。具体的には、閾値処理部５０３は、閾値用いて、ＩＥＤ確率マップ算出部５０２において得られたＩＥＤの確率マップからＩＥＤの確率が高い時間とセンサの領域の絞り込みを行う。

なお、閾値処理部５０３においては、閾値を高くすればするほど、よりＩＥＤらしい波形のみを選択できるが、ダイポール推定を行うＩＥＤそのものの数が少なくなってしまう。一方、閾値処理部５０３においては、閾値を低くした場合、誤検出は増えるが、より多くのＩＥＤを検出することができる。さらに、脳磁図のデータは個人差の影響が強いため、閾値を下げなければ、ＩＥＤがほとんど検出できない可能性もある。

そこで、閾値処理部５０３においては、ＵＩ（User Interface）と連携させてインタラクティブに検出するＩＥＤの数を制御する構成としてもよい。

閾値処理部５０３は、閾値として、基本的には「０．８」などのプリセット値を用いる。

ここで、図８は閾値処理後のＩＥＤ確率マップの一例を示す図である。図８に示す例においては、閾値処理部５０３は、閾値として「０．５」を用いた。図８に示す２次元マップであるＩＥＤ確率マップ６０５では、５０％以上の確率でＩＥＤが存在している箇所を抽出したことに相当する。

図５に戻り、次に、後処理部５０４は、閾値処理部５０３で閾値処理したＩＥＤの確率マップ６０５からダイポール推定に使用するためのサンプリング時点およびセンサの抽出を行う後処理を実行する（ステップＳ５）。

具体的には、閾値処理部５０３を適用すると、所定の確率以上の値を持つ時刻とセンサの領域とが抽出される。ピーク検出、もしくはＩＥＤの検出時点を明確にする目的で、抽出されたセンサの領域に対して、ガウシアンフィルタ等のフィルタをピーク検出の前に適用するようにしてもよい。抽出されたセンサの領域に対して、ガウシアンフィルタ等のフィルタをピーク検出の前に適用した後にピーク検出を行うことで、ＩＥＤの時点が決定される。ピーク検出を行う際には、抽出されたセンサの領域を含む確率マップをセンサ方向に加算平均した１次元の確率マップに対して行ってもよい。加算平均を行わない場合は、センサごとにピーク検出を行い、その時点を平均することでＩＥＤの時刻を決定するようにしてもよい。

閾値処理部５０３は、検出されたピーク時点における所定の確率以上の値をもつセンサをダイポール推定に使用するセンサとして抽出する。このとき、ダイポール推定を行う際のセンサの数が少なくなると、ダイポール推定解の安定性が低くなることが知られている。そのためここでは、選択されたセンサ数が規定の数に満たない場合には、そのＩＥＤ時点ではダイポール推定を行わないようにすることも可能である。

また、別のセンサ抽出方法として、予めセンサのグループを定義し、抽出されたセンサをそのセンサが属するグループのすべてのセンサへと拡張することでダイポール推定解の安定性を高めることもできる。ここでグループとは、側頭葉付近や前頭葉付近といった解剖学的な基準と照らして設定することを基本とするが、単純に任意の数を含む近傍領域を一つのグループとしてもよい。

ここで、図１２はセンサ拡張方法の一例を示す概念図である。図１２（ａ）は、被検者の頭部２００と、頭部２００を覆うように配置された脳磁図のセンサ２１０と、予め設定されたセンサ２１０のグループ２２０と、を表している。図１２（ｂ）において黒丸で表されたセンサ２１０ａは、閾値処理された段階で抽出されたセンサを表している。これらの閾値処理された段階で抽出されたセンサ２１０ａのみを用いてダイポール推定を行うことも可能だが、前述の通り、ダイポール推定解の安定性を高めるために、センサの拡張を実施することも可能である。図１２（ｃ）は、センサ数の拡張を行った結果である。すなわち、図１２（ｃ）に示すように、閾値処理された段階で抽出されたセンサ２１０ａをそのセンサ２１０ａが属するグループ２２０のすべてのセンサ２１０へと拡張する。これにより、ダイポール推定解の安定性を高めることができる。

また、後処理部５０４は、ＩＥＤ確率マップが複数のピークを持つ場合（病変部が２つ存在する場合などにみられるケース）には、マップを分離するための後述の後処理を行うこともできる。

近傍領域の設定の別の例として、選択されたセンサ座標から凸包を構成し、次に凸包を所定の閾値でモルフォロジー変換を行って空間を拡張し、拡張された空間に含まれるセンサを一つのグループとするなどが挙げられる。

また、選択センサデータの磁場の最大値を持つセンサと（湧き出し）、最小値を持つセンサ（吸い込み）の空間座標の中点を基準に最近傍の点から順に取得してもよい。

また、センサデータから取得した磁場情報から、所定の磁場閾値で空間的な等磁場線図を再構成し、選択されたセンサの磁場情報から、既に選択されたセンサの磁場情報の近傍でかつ同じ等磁場線に属するセンサを選択しても良い。

また、近傍領域の磁場の空間分布を用いた決定手段と解剖的な基準との組み合わせによるグループを構成することによって、空間的かつ解剖的な視点でのセンサグループを構成しても良い。

最後に、ダイポール推定部５０５は、ダイポール推定部５０５以前の処理で得られたセンサとＩＥＤの起始部のサンプリング時点を用いてダイポール推定（非特許文献３を参照）を行う（ステップＳ６）。ここでは、ダイポール推定以外にも、最小ノルム法（非特許文献６を参照）やLCMV Beamformer法（非特許文献７を参照）といった空間フィルタ法による電流源解析法で代替することが可能である。空間フィルタ法は、予め脳内に数万程度のダイポールを配置し、各ダイポールのおける電流の経時変化を求めることができる手法である。

ここで、図９は後処理後のＩＥＤ確率マップの一例を示す図である。後処理部５０４は、図８で得られた閾値処理後のＩＥＤ確率マップに対して、ピーク検出と閾値処理とを実行することで、ダイポール推定のサンプリング時点およびセンサの抽出とを行う。図９に示す例では、後処理部５０４は、ピーク検出と閾値処理とを実行することで、指定されたダイポール推定のサンプリング時点５１１と、センサを抽出したチャネル５１２とを得る。

ダイポール推定部５０５は、ダイポール推定のサンプリング時点５１１と、センサを抽出したチャネル５１２とを用いてダイポール推定を行う。

なお、図９に示す例は、スパイクが１つしかない場合を模擬しているが、二つ以上ある場合には、ｗａｔｅｒｓｈｅｄアルゴリズム（非特許文献８を参照）やランダムウォーク分割アルゴリズム（非特許文献９を参照）などを用いてマップを分離する処理を後処理部５０４で適用することで、単一のＩＥＤのみを対象としてダイポール推定を行うことが可能となる。

ここで、図１０はＩＥＤ確率マップを用いたＵＩ例を示す図である。図１０に示すように、上述のようにして算出された２次元マップであるＩＥＤ確率マップを元の波形に重畳することで、学習済みネットワークが波形のどこに注目したか俯瞰することが可能である。これによって、推定てんかん病変部の局在性の診断に悪影響を及ぼす誤検出などの推定結果に悪影響を与えるデータを除外できる。表示の際には、全センサをそのまま表示してもよいが、手動で選択したセンサのみの表示や、予め設定したセンサグループ単位で表示することも可能である。

このように表示することにより、推定結果と元々の波形と見比べることで、より正確なＩＥＤの時刻およびセンサの選択が可能となるので、より正確なてんかん病変の局在性診断が可能となるという効果がある。

なお、図１０に示す例においては、推定された確率に応じて濃淡が変化するように描画しているが、確率に応じてカラーマップを設定することも可能である。例えば、確率０であれば黒、高くなるにつれて青から赤、確率１であれば白となるように線形に色を変化させるようなカラーマップを適用してもよい。

特に、ＩＥＤ確率マップ算出部５０２が、機械学習の際に等価電流双極子法において理想的なＩＥＤの時点と絞り込んだセンサを正解データとして学習することにより、ダイポール推定の高精度化を図ることができる。図１０に示す例は、機械学習を用いて、てんかんのＩＥＤの起始部とそのＩＥＤが出現しているセンサを抽出した例である。図１０中で濃く表現されている箇所が、機械学習が推定した箇所を表している。

なお、上記は脳磁計の例で説明したが、これに限られることなく、本発明は脳波計についても適用することが可能である。

このように本実施形態によれば、ＩＥＤの確率マップを用いることで、等価電流双極子法に必要なＩＥＤの起始時刻の検出とセンサの選択が可能になるので、より正確に特徴的な波形情報（ＩＥＤ）が現れている時刻の決定とセンサの抽出が可能になり、等価電流双極子法を用いたてんかん病変部の局在性診断が可能となる。

また、機械学習を用いた自動ＩＥＤ探索手法を用いることで、ＩＥＤの探索・起始部の決定を自動化させることができる。

ここで、ＩＥＤ確率マップ算出部５０２が、ＩＥＤの確率マップの生成において機械学習以外の手法を用いる場合について説明する。

図１１は、機械学習以外の手法を用いたＩＥＤの確率マップの算出例を示す図である。図１１に示す例においては、ＩＥＤ確率マップ算出部５０２は、テンプレートマッチングを用いている。ＩＥＤ確率マップ算出部５０２は、脳磁図におけるＩＥＤの形態的な特徴について調査された文献（非特許文献１０を参照）より作成したＩＥＤの形状を一般化したテンプレートを用いて、ＩＥＤの検出を行う。

図１１に示すように、ＩＥＤ確率マップ算出部５０２は、脳磁図データ６０４の一つのセンサの脳磁データの時系列に対してテンプレート６０１～６０３のようなテンプレートを走査していき、類似度よって２次元マップであるＩＥＤ確率マップ６０５を作成する。ここで作られるマップは厳密には確率マップではないが、目的上マップとしては利用できるので便宜上確率マップとする。

類似度の算出には、正規化相互相関（ＮＣＣ：Normalized Cross-Correlation）や、零平均正規化相互相関（ＺＮＣＣ：Zero-mean Normalized Cross-Correlation）など画像におけるテンプレートマッチングの類似度算出方法を１次元化したものや、単純にフィルタとの畳み込みを正規化したものなどを用いることで計算できる。ＩＥＤの検出にはテンプレートマッチング以外にも、ＩＥＤのベースラインからの急峻な遷移を検出することで導出する手法などいくつか適用することが可能である。

なお、ＩＥＤ確率マップ算出部５０２において機械学習を用いない場合では、ピーク位置でのダイポール推定は好ましくない結果をもたらす可能性があるため、少しサンプル時点を前にずらした時点で行うなどの後処理が必要である。

このようにアルゴリズムを用いた自動ＩＥＤ探索手法を用いることで、ＩＥＤの探索・起始部の決定を自動化させることができる。

また、上述の実施形態において、生体信号計測システム１の各機能部の少なくともいずれかがプログラムの実行によって実現される場合、そのプログラムは、ＲＯＭ等に予め組み込まれて提供される。また、上述の実施形態に係る生体信号計測システム１で実行されるプログラムは、インストール可能な形式または実行可能な形式のファイルでＣＤ－ＲＯＭ、フレキシブルディスク（ＦＤ）、ＣＤ－Ｒ（Compact Disc Recordable）、ＤＶＤ（Digital Versatile Disc）等のコンピュータで読み取り可能な記録媒体に記録して提供するように構成してもよい。

また、上述の実施形態の生体信号計測システム１で実行されるプログラムを、インターネット等のネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることにより提供するように構成してもよい。

また、上述の実施形態の生体信号計測システム１で実行されるプログラムを、インターネット等のネットワーク経由で提供または配布するように構成してもよい。また、上述の実施形態の生体信号計測システム１で実行されるプログラムは、上述した各機能部のうち少なくともいずれかを含むモジュール構成となっており、実際のハードウェアとしてはＣＰＵがＲＯＭ等からプログラムを読み出して実行することにより、上述の各機能部が主記憶装置上にロードされて生成されるようになっている。

６０５ ２次元マップ

特許第３１０７０２０号公報

E. Javier, H. Roberto, A. Daniel, F. Alberto, and L. C. Miguel, "Artifact removal in magnetoencephalogram background activity with independent component analysis," IEEE Trans Biomed Eng, vol. 54, no. 11, pp. 1965-1973, 2007. K. Sekihara, Y. Kawabata, S. Ushio, S. Sumiya, S. Kawabata, Y. Adachi, and S. S. Nagarajan, "Dual signal subspace projection (DSSP): a novel algorithm for removing large interference in biomagnetic measurements," Journal of Neural Engineering, vol. 13, no. 3, p. 036007, 2016. M. Scherg, "Fundamentals of dipole source potential analysis," in Auditory Evoked Magnetic Fields and Potentials, M. Hoke, F. Grandori, and G. L. Romani, Eds. Basel, Switzerland: Karger, 1989, vol. 6 A. Ossadtchi, S. Baillet, J. Mosher, D. Thyerlei, W. Sutherling, and R. Leahy, "Automated interictal spike detection and source localization in magnetoencephalography using independent components analysis and spatio-temporal clustering," Clinical Neurophysiology, vol. 115, no. 3, pp. 508-522, 2004. O. Ronneberger P. Fischer and T. Brox "U-net: Convolutional networks for biomedical image segmentation," Proc. Int. Conf. Medical Image Comput. Comput.-Assisted Intervention, pp. 234-241 2015. K. Sekihara, M. Sahani, and S. S. Nagarajan, "Localization bias and spatial resolution of adaptive and non-adaptive spatial filters for MEG source reconstruction," NeuroImage, vol. 25, no. 4, pp. 1056-1067, 2005. B. V. Veen, W. V. Drongelen, M. Yuchtman, and A. Suzuki, "Localization of brain electrical activity via linearly constrained minimum variance spatial filtering," IEEE Transactions on Biomedical Engineering, vol. 44, no. 9, pp. 867-880, 1997. J. Cousty, G. Bertrand, L. Najman, and M. Couprie, "Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 31, no. 8, pp. 1362-1374, 2009. L. Grady, "Random Walks for Image Segmentation," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 28, no. 11, pp. 1768-1783, 2006. R. Nowak, M. Santiuste, and A. Russi, "Toward a definition of MEG spike: Parametric description of spikes recorded simultaneously by MEG and depth electrodes", Seizure, vol. 18, no. 9, pp. 652-655, 2009.

Claims (14)

1. 複数のセンサにより取得した個々の波形データについて、少なくとも１つ以上の特徴的な波形情報と比較する比較ステップと、
   前記波形データの区間と前記特徴的な波形情報との相関の割合に基づいて、前記波形データの少なくとも一定区間における特徴的な波形情報の出現確度を決定する決定ステップと、
   前記決定された特徴的な波形情報の出現確度に対し、閾値を用いてセンサの領域を抽出し、該領域に対してフィルタを適用して特徴的な波形情報の時刻と前記センサの領域とを絞り込む絞込ステップと、
   前記絞込ステップで絞り込まれた時刻及びセンサの領域を用いて、特徴的な波形情報に一致する区間の出現した時刻及びセンサのインデックスを特定する特定ステップと、
   を含むことを特徴とする波形発生特定方法。
2. 前記決定ステップは、前記波形データの少なくとも一定区間における特徴的な波形情報の出現確度を表す２次元マップを表示する、
   ことを特徴とする請求項１に記載の波形発生特定方法。
3. 前記決定ステップは、前記２次元マップを前記複数のセンサにより取得した個々の波形データに重畳表示する、
   ことを特徴とする請求項２に記載の波形発生特定方法。
4. 前記決定ステップは、事前に学習した機械学習のモデルを用いて特徴的な波形情報の確率マップを算出する、
   ことを特徴とする請求項１に記載の波形発生特定方法。
5. 前記決定ステップは、機械学習の際に、等価電流双極子法において理想的な特徴的な波形情報の時点と、絞り込んだセンサのインデックスとを正解データとして学習させる、
   ことを特徴とする請求項４に記載の波形発生特定方法。
6. 前記特定ステップは、予めセンサのグループを定義し、前記確率マップにおいて所定の確率以上の値をもつセンサをそのセンサが属するグループのすべてのセンサへと拡張して、センサを絞り込む、
   ことを特徴とする請求項４に記載の波形発生特定方法。
7. 前記決定ステップは、前記２次元マップを、特徴的な波形情報の確率に応じて色が変化するようなカラー表示とする、
   ことを特徴とする請求項２に記載の波形発生特定方法。
8. 複数のセンサにより取得した個々の波形データについて、少なくとも１つ以上の特徴的な波形情報と比較する比較ステップと、
   前記波形データの区間と前記特徴的な波形情報との相関の割合に基づいて、前記波形データの少なくとも一定区間における特徴的な波形情報の出現確度を決定する決定ステップと、
   前記決定された特徴的な波形情報の出現確度に対し、閾値を用いてセンサの領域を抽出し、該領域に対してフィルタを適用して特徴的な波形情報の時刻と前記センサの領域とを絞り込む絞込ステップと、
   前記絞込ステップで絞り込まれた時刻及びセンサの領域を用いて、特徴的な波形情報に一致する区間の出現した時刻及びセンサのインデックスを特定する特定ステップと、
   をコンピュータに実行させるためのプログラム。
9. 前記決定ステップは、前記波形データの少なくとも一定区間における特徴的な波形情報の出現確度を表す２次元マップを表示する、
   ことを特徴とする請求項８に記載のプログラム。
10. 前記決定ステップは、前記２次元マップを前記複数のセンサにより取得した個々の波形データに重畳表示する、
    ことを特徴とする請求項９に記載のプログラム。
11. 前記決定ステップは、事前に学習した機械学習のモデルを用いて特徴的な波形情報の確率マップを算出する、
    ことを特徴とする請求項８に記載のプログラム。
12. 前記決定ステップは、機械学習の際に、等価電流双極子法において理想的な特徴的な波形情報の時点と、絞り込んだセンサのインデックスとを正解データとして学習させる、
    ことを特徴とする請求項１１に記載のプログラム。
13. 前記特定ステップは、予めセンサのグループを定義し、前記確率マップにおいて所定の確率以上の値をもつセンサをそのセンサが属するグループのすべてのセンサへと拡張して、センサを絞り込む、
    ことを特徴とする請求項１１に記載のプログラム。
14. 前記決定ステップは、前記２次元マップを、特徴的な波形情報の確率に応じて色が変化するようなカラー表示とする、
    ことを特徴とする請求項９に記載のプログラム。

Images