特許 7161756 生体信号解析装置及び生体信号解析方法

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

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

(11)【特許番号】

(24)【登録日】2022-10-19

(45)【発行日】2022-10-27

(54)【発明の名称】生体信号解析装置及び生体信号解析方法

(51)【国際特許分類】

   A61B 5/374 20210101AFI20221020BHJP

   A61B 10/00 20060101ALI20221020BHJP

【ＦＩ】

A61B5/374

A61B10/00 H

A61B5/374 ZDM

【請求項の数】 6

(21)【出願番号】P 2018231225

(22)【出願日】2018-12-10

(65)【公開番号】P2020092753

(43)【公開日】2020-06-18

【審査請求日】2021-10-12

(73)【特許権者】

【識別番号】504136568

【氏名又は名称】国立大学法人広島大学

(74)【代理人】

【識別番号】100196380

【弁理士】

【氏名又は名称】森 匡輝

(72)【発明者】

【氏名】辻 敏夫

(72)【発明者】

【氏名】曽 智

(72)【発明者】

【氏名】古居 彬

(72)【発明者】

【氏名】大西 亮太

(72)【発明者】

【氏名】秋山 倫之

(72)【発明者】

【氏名】竹内 章人

【審査官】外山 未琴

(56)【参考文献】

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

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

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

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

【文献】Furui, Akira, Hideaki Hayashi, and Toshio Tsuji，An EMG pattern classification method based on a mixture of variance distribution models，2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)，2018年10月29日

【文献】古居 彬、辻 敏夫，表面筋電位の分散分布モデルと生体ゆらぎ評価への応用，生理心理学と精神生理学，2018年，36巻 2号，66頁

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

Ａ６１Ｂ ５／００－５／３９８

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

【請求項1】

生体信号を取得する生体信号取得部と、
前記生体信号取得部で取得された複数の生体信号のゆらぎを表す１次元の特徴パラメータを算出する演算部と、を備え、
前記演算部は、
前記複数の生体信号を混合した多次元分散分布モデルである尺度混合モデルを用いて、前記特徴パラメータを算出し、
前記特徴パラメータは、
前記多次元分散分布モデルを表す分散共分散行列の分布として用いられる逆ウィシャート分布の自由度を含む、
ことを特徴とする生体信号解析装置。

【請求項2】

前記生体信号は脳波である、
ことを特徴とする請求項１に記載の生体信号解析装置。

【請求項3】

前記生体信号に基づいて、異常脳波を伴う疾患を検出する、
ことを特徴とする請求項１又は２に記載の生体信号解析装置。

【請求項4】

前記生体信号を、複数の周波数フィルタを用いて、所定の周波数帯域の信号に分解するフィルタ部を備え、
前記演算部は、
各周波数帯域に分解された前記生体信号について特徴パラメータを算出する、
ことを特徴とする請求項１から３のいずれか一項に記載の生体信号解析装置。

【請求項5】

前記演算部は、
所定の時間幅を有する窓を設定し、窓の開始時刻を所定のずらし幅でずらしながら、設定された窓の時間幅の前記生体信号について特徴パラメータを算出する、
ことを特徴とする請求項１から４のいずれか一項に記載の生体信号解析装置。

【請求項6】

複数の生体信号を取得するステップと、
取得された前記複数の生体信号のゆらぎを表す１次元の特徴パラメータを算出するステップと、を含み、
前記特徴パラメータは、前記複数の生体信号を混合した多次元分散分布モデルである尺度混合モデルを用いて算出され、前記多次元分散分布モデルを表す分散共分散行列の分布として用いられる逆ウィシャート分布の自由度を含む、
ことを特徴とする生体信号解析方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、生体信号解析装置及び生体信号解析方法に関する。

【背景技術】

【0002】

従来、脳波、脈拍等の生体信号を用いて、被験者の身体の状態を分析する装置が開発されている。例えば、特許文献１の診断支援装置では、被験者の頭部上の所定部位で測定された脳波の時系列データを解析することにより、被験者が起こしたけいれんが、けいれん重積型脳症によるけいれんか熱性けいれんかを判別する。

【先行技術文献】

【特許文献】

【0003】

【文献】特開２０１７－２０９２６７号公報

【発明の概要】

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

【0004】

特許文献１の診断支援装置では、複数の部位で測定された脳波を個々に分析して、判別情報を生成する。したがって、個々の脳波に含まれるノイズの影響を受けやすく、判定精度は低いものとなる。

【0005】

また、特許文献１の診断支援装置では脳波を周波数解析して被験者の状態を判断する。特許文献１のように、生体信号の周波数又は振幅の大きさに着目した解析方法では、被験者の身体的特徴による個体差や生体信号に含まれるノイズの影響を受けやすいため解析精度が低い。

【0006】

本発明は、上述の事情に鑑みてなされたものであり、複数の生体信号を混合して被験者の身体の状態を分析することにより、判定精度の高い生体信号解析装置及び生体信号解析方法を提供することを目的とする。

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

【0007】

上記目的を達成するために、この発明の第１の観点に係る生体信号解析装置は、
生体信号を取得する生体信号取得部と、
前記生体信号取得部で取得された複数の生体信号のゆらぎを表す１次元の特徴パラメータを算出する演算部と、を備え、
前記演算部は、
前記複数の生体信号を混合した多次元分散分布モデルである尺度混合モデルを用いて、前記特徴パラメータを算出し、
前記特徴パラメータは、
前記多次元分散分布モデルを表す分散共分散行列の分布として用いられる逆ウィシャート分布の自由度を含む。

【0010】

また、前記生体信号は脳波である、
こととしてもよい。

【0011】

また、前記生体信号に基づいて、異常脳波を伴う疾患を検出する、
こととしてもよい。

【0012】

また、前記生体信号を、複数の周波数フィルタを用いて、所定の周波数帯域の信号に分解するフィルタ部を備え、
前記演算部は、
各周波数帯域に分解された前記生体信号について特徴パラメータを算出する、
こととしてもよい。

【0013】

また、前記演算部は、
所定の時間幅を有する窓を設定し、窓の開始時刻を所定のずらし幅でずらしながら、設定された窓の時間幅の前記生体信号について特徴パラメータを算出する、
こととしてもよい。

【0014】

また、本発明の第２の観点に係る生体信号解析方法は、
複数の生体信号を取得するステップと、
取得された前記複数の生体信号のゆらぎを表す１次元の特徴パラメータを算出するステップと、を含み、
前記特徴パラメータは、前記複数の生体信号を混合した多次元分散分布モデルである尺度混合モデルを用いて算出され、前記多次元分散分布モデルを表す分散共分散行列の分布として用いられる逆ウィシャート分布の自由度を含む。

【発明の効果】

【0015】

本発明の生体信号解析装置及び生体信号解析方法によれば、複数の生体信号を混合して、生体信号のゆらぎを表す１次元パラメータを算出するので、容易にかつ精度よく被験者の状態変化を分析することが可能である。

【図面の簡単な説明】

【0016】

【図1】本発明の実施の形態に係る生体信号解析装置のブロック図である。

【図2】実施の形態に係る生体信号解析の流れを示すフローチャートである。

【図3】実施の形態に係る脳波の測定箇所を示す概念図である。

【図4】焦点発作を生じた場合の脳波の例を示す図であり、（Ａ）が左半球の脳波、（Ｂ）が右半球の脳波である。

【図5】図４の焦点発作を生じた場合の脳波の解析結果を示すコンター図であり、（Ａ）が左半球の脳波の解析結果、（Ｂ）が右半球の脳波の解析結果である。

【図6】欠神発作を生じた場合の脳波の例を示す図であり、（Ａ）が左半球の脳波、（Ｂ）が右半球の脳波である。

【図7】図６の欠神発作を生じた場合の脳波の解析結果を示すコンター図であり、（Ａ）が左半球の脳波の解析結果、（Ｂ）が右半球の脳波の解析結果である。

【図8】ミオクロニー発作を生じた場合の脳波の例を示す図であり、（Ａ）が左半球の脳波、（Ｂ）が右半球の脳波である。

【図9】図８のミオクロニー発作を生じた場合の脳波の解析結果を示すコンター図であり、（Ａ）が左半球の脳波の解析結果、（Ｂ）が右半球の脳波の解析結果である。

【図10】本発明と二乗平均平方根を用いた従来法との癲癇発作の識別率を示すグラフである。

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

【0017】

以下、図を参照しつつ、本発明の実施の形態に係る生体信号解析装置１について説明する。図１のブロック図に示すように、生体信号解析装置１は、センサ１１、制御ユニット２０を備える。

【0018】

センサ１１は、取得する生体信号の種類に合わせて選択される検出器である。取得する生体信号は、例えば脳波、脈拍等であり、特に限定されない。

【0019】

制御ユニット２０は、例えばコンピュータ装置であり、図１に示すように、制御部２１、記憶部２２、表示部２３、入力部２４を備える。

【0020】

制御部２１は、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、水晶発振器等から構成されており、生体信号解析装置１の動作を制御するとともに、センサ１１で取得された生体信号に基づいて特徴パラメータを算出し、被験者の身体状態を判定する。制御部２１は、制御部２１のＲＯＭ、記憶部２２等に記憶されている各種動作プログラム及びデータをＲＡＭに読み込んでＣＰＵを動作させることにより、図１に示される制御部２１の各機能を実現させる。これにより、制御部２１は、生体信号取得部２１１、フィルタ部２１２、演算部２１３として動作する。

【0021】

生体信号取得部２１１は、センサ１１を制御して、被験者の生体信号を取得する。また、取得した生体信号をフィルタ部２１２へ送信する。本実施の形態では、複数のセンサ１１を用いて、同時に複数の生体信号が取得される。

【0022】

フィルタ部２１２は、特定の周波数帯域の信号を通過させる周波数フィルタを備え、生体信号取得部２１１で取得した生体信号を特定の周波数帯域の信号へ変換する。

【0023】

演算部２１３は、フィルタ部２１２から受信した複数の生体信号に基づいて、被験者の身体の状態を示す１次元の特徴パラメータを算出する。特徴パラメータは、複数の生体信号のゆらぎを示すものであり、ゆらぎの大きさ、すなわち時間的な変化の大きさによって被験者の身体状態の変化を表す。特徴パラメータの詳細な算出方法については、後述する。

【0024】

記憶部２２は、ハードディスク、フラッシュメモリ等の不揮発性メモリであり、複数の生体信号から特徴パラメータを算出する演算アルゴリズム、算出された特徴パラメータ等を記憶する。

【0025】

表示部２３は、コンピュータ装置である制御ユニット２０に備えられた表示用デバイスであり、例えば液晶パネルである。表示部２３は、センサ１１で取得された生体信号、演算部２１３で算出された特徴パラメータ等を表示する。

【0026】

入力部２４は、生体信号取得の開始、終了指示、特徴パラメータ算出のための条件変更等を入力するための入力デバイスである。入力部２４は、制御ユニット２０に備えられたキーボード、マウス等である。

【0027】

続いて、生体信号解析装置１を用いた生体信号解析方法について、図２のフローチャートを参照しつつ、具体的に説明する。本実施の形態では、複数チャンネルの脳波を生体信号として取得し、被験者が癲癇の発作を生じているか否かを解析する場合を例として説明する。

【0028】

医師等のユーザは、図３に示すように、被験者の頭部に複数のセンサ１１を取り付ける。そして、入力部２４を操作して、脳波の取得を開始する。生体信号取得部２１１は、脳波の取得開始指示により、センサ１１を制御して、被験者の脳波を取得する（ステップＳ１１）。生体信号取得部２１１は、取得した脳波データをフィルタ部２１２へ送信する（ステップＳ１２）。

【0029】

フィルタ部２１２は、受信した各チャンネルの脳波データを、周波数フィルタによって、所定の周波数帯域の信号に変換する（ステップＳ１３）。本実施の形態に係るフィルタ部２１２は、取得された脳波データを、α帯域（８～１２Ｈｚ）、β帯域（１３～２４Ｈｚ）、γ帯域（２５～１００Ｈｚ）、δ帯域（１～３Ｈｚ）、θ帯域（４～７Ｈｚ）の５つの帯域の信号に分解する。そして、演算部２１３は、フィルタ部２１２から受信した各周波数帯域の信号について、脳波のゆらぎを表す１次元パラメータである特徴パラメータを算出する。これにより、癲癇発作の種類等について、より詳細な分析が可能となる。

【0030】

以下、特徴パラメータの算出方法について説明する。演算部２１３は、分解された各帯域の脳波について、同一の時刻から所定の時間幅を持つ窓を設定する（ステップＳ１４）。

【0031】

演算部２１３は、分解された帯域ごとに、信号のゆらぎを表す特徴パラメータを算出する（ステップＳ１５）。より具体的には、演算部２１３は、まず、ステップＳ１４で設定した窓内の、各チャンネルの脳波を混合して解析するための尺度混合モデルを生成する。本実施の形態に係る尺度混合モデルでは、センサ１１であるＬ極（Ｌは自然数）の電極で計測された脳波ｘをｘ∈Ｒ^Ｌ、平均ゼロの確率変数と考え、分散共分散行列Σ∈Ｒ^Ｌ×Ｌの多変量正規分布を用いて表現する。ここで、Σは確率変数であり、Σの分布は自由度ν∈Ｒ^１と尺度行列Ψ∈Ｒ^Ｌ×Ｌによって特徴付けられる。

【0032】

上記の分散共分散行列Σが与えられた下での脳波ｘの条件付き分布Ｐ（ｘ｜Σ）は、以下の式（１）に示す多変量正規分布で表される。

【数1】

【0033】

分散共分散行列Σの分布として、以下の式（２）に示す、多変量正規分布の共役事前分布である逆ウィシャート分布Ｐ（Σ；ν，Ψ）を仮定する。

【数2】

【0034】

ここで、ν、Ψは、逆ウィシャート分布を決定するパラメータであり、νは自由度、Ψは尺度行列と呼ばれる。

【0035】

ｘの周辺分布Ｐ（ｘ）は、以下の式（３）、（４）で表される。

【数3】

【0036】

ただし、Δはマハラノビス距離の平方で、以下の式（５）で表される。

【数4】

【0037】

上記の式（１）、（３）より、Ｐ（ｘ）は、分散共分散行列の異なる多変量正規分布を無限個足し合わせた分布である。本実施の形態では、この多次元分散分布モデルＰ（ｘ）を尺度混合モデルとして用いる。

【0038】

続いて、計測された脳波ｘから、分散共分散行列Σのパラメータとして定義される特徴パラメータを推定する方法について説明する。

【0039】

本実施の形態では、観測されたＮサンプルの脳波データＸ＝｛ｘ_ｎ∈Ｒ^Ｌ；ｎ＝１，・・・、Ｎ｝に基づいて分散共分散行列Σのパラメータν、Ψを推定する。観測されたサンプル系列が尤も生起され易いパラメータは、次式で表す周辺分布の尤度を最大化することにより推定できる。

【数5】

【0040】

本実施の形態では、ＥＭ（Expectation-Maximization）アルゴリズムを用いて上記の推定を行う。

【0041】

簡単のため、式（４）において、パラメータをν’＝ν－Ｌ＋１，Ψ’＝Ψ／（ν－Ｌ＋１）と置き換えると、周辺尤度Ｐ（ｘ_ｎ）は以下の式（７）のように書き換えることができる。

【数6】

【0042】

ただし、Δ’は以下の式（８）で表される。

【数7】

ここで、式（７）はStudentのｔ分布と等価である。

【0043】

式（７）、（８）を用いて、以下の手順（ａ）～（ｄ）で周辺尤度Ｐ（ｘ_ｎ）が最大となるパラメータν’、Ψ’を推定する。
（ａ）各パラメータに任意の初期値を与える。
（ｂ）現在のパラメータ値を使って、以下の式（９）に示す重みパラメータｗ_ｎを計算する（Ｅステップ）。

【数8】

【0044】

（ｃ）手順（ｂ）で求めたｗ_ｎを用いて、以下の式（１０）により、Ψ’を更新する（Ｍステップ）。

【数9】

【0045】

ν’については、対数周辺尤度を最大化するν’を線形探索によって求める。

【数10】

【0046】

（ｄ）更新後のパラメータを用いて対数周辺尤度を算出し、収束基準を満たすまで、上記の手順（ｂ）及び（ｃ）を反復する。最後に、推定されたｔ分布のパラメータであるν’、Ψ’を、元の分散共分散行列Σのパラメータν、Ψに変換する。

【0047】

上記の手順（ａ）～（ｄ）により、計測した脳波ｘから分散共分散行列Σのパラメータを推定する。

【0048】

推定されたパラメータであるνは、ｔ分布の自由度に対応し、分布のガウス性を特徴付ける。本実施の形態の尺度混合モデルにおいては、分散に相当するパラメータの確率的ゆらぎによって、分布のガウス性が変化すると考えられる。すなわち、観測された脳波ｘから自由度νを推定することにより、脳波ｘの分散共分散行列Σのゆらぎを評価することができる。以上のようにして、脳波のゆらぎを表す１次元パラメータである特徴パラメータとして自由度νを算出する。

【0049】

特徴パラメータとしての自由度νを算出した後、解析を継続する場合（ステップＳ１６のＹＥＳ）、窓の開始時刻を所定のずらし幅で遅らせて（ステップＳ１７）、ステップＳ１５に戻り、特徴パラメータである自由度νの算出を繰り返す。

【0050】

解析を継続しない場合（ステップＳ１６のＮＯ）、制御部２１は、算出された自由度νを時系列データとして記憶部２２へ記憶させる。また、制御部２１は、算出された自由度νの時系列データを、被験者の身体状態の変化を表す特徴パラメータとして、表示部２３に表示させる（ステップＳ１８）。そして、生体信号解析装置１は、解析処理を終了する。

【0051】

図４（Ａ）、（Ｂ）は、癲癇患者の脳波データの例であり、図４（Ａ）は左半球の脳波データ、図４（Ｂ）は右半球の脳波データを表している。図５（Ａ）、（Ｂ）は、図４（Ａ）、（Ｂ）の各データを上述の解析方法によって解析した解析結果を示す図である。図５（Ａ）、（Ｂ）では、自由度νの逆数である１／νを特徴パラメータとして、周波数帯域ごとに時系列のコンター図として示している。計測チャンネル数は、図３に示す１６チャンネルであり、計測時間は５０±２秒、サンプリング周波数は５００Ｈｚ、窓幅は１０秒、ずらし幅は１秒である。

【0052】

図４（Ａ）、（Ｂ）に示す脳波は、被験者が焦点発作を発生している場合のデータであり、図５（Ａ）、（Ｂ）に示すように、発作の発生を明確に判別できることがわかる。

【0053】

図６（Ａ）、（Ｂ）及び図７（Ａ）、（Ｂ）は別の癲癇患者の脳波の例である。図６（Ａ）、（Ｂ）は、欠神発作を発生している場合のデータであり、図７（Ａ）、（Ｂ）に示すように、発作の種類が異なる場合でも、本発明に係る生体信号解析により、発作の発生を明確に判別できることがわかる。

【0054】

図８（Ａ）、（Ｂ）及び図９（Ａ）、（Ｂ）は、さらに別の癲癇患者の脳波の例である。図８（Ａ）、（Ｂ）は、ミオクロニー発作を発生している場合のデータであり、図９（Ａ）、（Ｂ）に示すように、発作が上記２つの発作の場合と異なる周波数帯域、異なる大きさの信号として表れる場合であっても、複数の周波数フィルタを適用することで、明確に判別できることがわかる。

【0055】

以上、説明したように、本発明に係る生体信号解析装置及び生体信号解析方法では、複数の生体信号を混合して、生体信号のゆらぎを表す１次元パラメータを算出するので、容易に、かつ精度よく被験者の身体状態の変化を分析することが可能である。

【0056】

すなわち、本発明の解析方法によって１次元パラメータとして解析できるので、被験者の頭部に取り付けられた複数のセンサ１１で取得される各チャンネルの脳波を医師が観察し、発作が生じているか否かを経験的に判定する方法に比べて、癲癇発作の発生等、被験者の身体状態の変化をより容易に分析することができる。

【0057】

また、複数の生体信号のゆらぎを解析できるので、個々の生体信号の振幅変化、周波数変化を解析する従来の方法に比べて、より精度よく癲癇発作等被験者の身体状態の変化を分析することができる。

【0058】

例として、γ帯域の信号に着目して、本発明の方法によって発作を識別した場合の識別率と、二乗平均平方根（ＲＭＳ：Root Mean Square）により振幅変化を特徴量として発作を識別した場合の識別率とを比較する。図１０に示すように、本発明の方法による識別率は、二重平均平方根を用いた場合の識別率よりも高いことがわかる。

【0059】

また、本実施の形態では、センサ１１で取得された生体信号をフィルタ部２１２の周波数フィルタによって、異なる周波数帯域を持つ複数の信号に分割して解析している。これにより、身体状態の変化の種類によって、その変化の表れやすい周波数帯域が異なる場合でも、適切に状態変化の発生を分析することができる。

【0060】

本実施の形態では、式（７）～（１１）に示すＥＭアルゴリズムによって特徴パラメータを推定することとしたが、これに限られない。特徴パラメータは、式（６）に示すモデルの周辺尤度を最大化可能な任意のアルゴリズムで推定できる。また、周辺尤度最大化以外のハイパーパラメータ推定法を適用することも可能である。

【0061】

本実施の形態では、被験者の脳波を生体信号として解析したが、これに限られない。例えば、脈拍、呼吸、発汗、心電図等の生体信号であってもよい。また、脳波と脈拍、呼吸と心電図等、複数の種類の生体信号を混合して解析することとしてもよい。これにより、身体の様々な部位の状態変化や複合的な状態変化を分析することができる。

【0062】

また、本実施の形態では、生体信号である被験者の脳波に基づいて、癲癇の発作を検出することとしたが、これに限られない。例えば、脳波を含む生体信号に基づいて、意識障害を伴う脳炎、糖尿病性昏睡、脳血管障害、脳腫瘍、アルツハイマー型痴呆等の異常脳波を伴う疾患を、検出することとしてもよい。また、被験者の覚醒状態、睡眠状態の別や睡眠深度などの脳状態を分析することとしてもよい。

【0063】

また、本実施の形態では、センサ１１で取得した生体信号である脳波を、制御ユニット２０で逐次解析することとしたが、これに限られない。例えば、過去に取得された生体信号を制御ユニット２０の生体信号取得部２１１に読み込んで、解析を実行し、現在の生体信号の解析結果と比較することとしてもよい。

【産業上の利用可能性】

【0064】

本発明は、複数の生体信号に基づいて身体状態の変化を分析する生体信号解析装置に好適である。特に、複数チャンネルの脳波を長時間観察することが求められる癲癇患者の発作の発生検出を行う生体信号解析装置に好適である。

【符号の説明】

【0065】

１ 生体信号解析装置、１１ センサ、２０ 制御ユニット、２１ 制御部、２１１ 生体信号取得部、２１２ フィルタ部、２１３ 演算部、２２ 記憶部、２３ 表示部、２４ 入力部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】