JP2011083393A - 睡眠ステージ自動判定の装置と方法およびそのためのコンピュータプログラム - Google Patents

Abstract

【課題】睡眠脳波データに含まれる被験者の個性によるバラツキ、電極位置の不確定性その他に起因する異常値、例外値、不規則値などを検出、低減すること、その結果として精度良く、完全自動で短時間に睡眠ステージの自動判定を行うこと。
【解決手段】収集された睡眠脳波信号を周波数解析、主成分分析を経て得られた特徴ベクトルに対して、各成分に対する百分位数統計手段の適用により１〜１０パーセンタイル値を基準とすることにより、または重み付けしたワンクラス・サポートベクトルマシンによりアウトライアを除去後のインライアデータの平均値を基準とすることにより、規格化後、多クラス識別処理を行って、判定結果を得る。
【選択図】図１

Description

本発明は睡眠時脳波を測定して睡眠ステージ（段階）を自動的に判定するための装置と方法およびそのためのコンピュータプログラムに関するものである。特に、人手を介することなく自動的に精度良く睡眠ステージを判定することのできるアルゴリズムを採用することを特徴とする自動判定装置とその方法に関する。

睡眠時脳波は医療現場では、通常、終夜睡眠ポリグラフィ検査（ＰＳＧ）によって、筋電図、眼電図、心電図、呼吸量などとともに、脳波図（ＥＥＧ）として計測されている。そして特別な訓練を受けた医療技術者による目視によって（手動で）、脳波の視察、判読が行われ、瘤状波や紡錘波といった特徴波を拠り所にしながら、２０秒ないし３０秒をエポック単位として、睡眠状態（覚醒、レム睡眠、浅い睡眠、深い睡眠など）のラベリングが行われてきた。例えば、７、８時間程度の一夜の睡眠データは８００〜１０００エポックもの膨大な脳波データを含むので全ての脳波の判読を行うには多大の労力と時間とを要することになる。このため、睡眠の測定、判定、診断を享受できる人数にはおのずから限度があり、コストもかかるので、検査方法としては一般化、普及していないのが実情である。

他方、前記のような本格的ＰＳＧ検査の代替手段として、一般人が独力で簡便に計測可能でありかつ睡眠脳波データの計測結果、評価結果も容易に入手できるような簡易脳波測定装置と睡眠状態判定方法が開発されてきており、いくつか実用化されている。
すなわち、人工知能を用いた機械学習方法による睡眠ステージの自動判定テクニックとして、代表的なのが人工ニューラルネットワーク（例えば特許文献１、非特許文献１）であり、これは波形の局所的特徴である瘤波、紡錘波などを直接的に抽出して評価するものであり、それなりに有効である。また、ステージ判定の手掛かりとして、α波、β波、θ波、δ波、σ波等の周波数帯域の強度（閾値）に基づいて判定を行う決定木手法（例えば特許文献２）が知られている。決定木手法の場合、しかしながら、被験者の個性の差異等に依っては十分な精度を自動では得難い場合がある。十分な精度を得ようとすれば一部手動での調整が必要な場合もあって、それにはスペシャリストの経験と熟練を要するなど、その分効率的でない問題があった。

そこで、本発明者らは他に代替可能で完全自動化可能な、機械学習法を用いた睡眠状態自動判定方法を模索、追究してきたところ、多クラス識別（分類）法としてサポートベクトルマシン（以下、ＳＶＭと略記）を採用した場合に、脳波データに不可避的に混入する不規則データなどを解析時に検出、排除、処理する困難と課題を克服することが可能となり、本発明に至ったものである。
尚、睡眠状態の自動検出方法として、特許文献３が謂わば、網羅的、概論的に、既存の技術である様々な機械学習テクニック、例えば前記の決定木、ニューラルネットワーク、さらにはＳＶＭなどにも言及はしているものの、具体的実施態様においても、それらテクニックを如何に適用すべきか、また前述のような脳波データに付随する様々な問題点に如何に対処するかなどの問題意識や解決について教示し示唆するものでは全くない。

特表２００１−５０３６６７号公報（対応：ＷＯ９８／２２０１９Ａ） 特開２００９−１１２４０２号公報 特表２００８−５４４７７２号公報（対応：ＷＯ２００６／１２１４５５Ａ１）

信学論（電子情報通信学会論文誌），Ｊ７６−Ａ，８，ｐｐ．１０５０−１０５８（１９９３）

脳波データの解析にとって考慮すべき問題のひとつは、先ず被験者の個性（睡眠構造）の差異の存在と日々の健康状態の変化に伴う変動のため、誤差の原因となり、延いては判定精度を劣化させることであり、それらを如何なる方法で特定、検出し排除して解析できるかが重要な鍵となる。

第２に、睡眠脳波データには、本来目的とする信号以外のノイズおよび記録における波形歪みであるアーチファクトが混入して睡眠ステージの判定に支障を来たす問題がある。アーチファクトには、生体に由来するアーチファクト、例えば心電図、筋電図、脈波、呼吸、発汗、眼球および眼瞼の動き、咳、歯軋り、体動（寝返り、トイレへの往復移動）など；電極から入るアーチファクト、例えば電極の装着不良、分極電圧の影響、電極コードの揺れ、電極コードの断線しかかりなど；外部から入るアーチファクト、例えば交流障害（ハム）、他の機器からの影響、洩れ電流、静電気など；環境的要素によるアーチファクト、例えば浅い睡眠時に周辺の音に反応する聴覚ノイズ、布ずれなどの触覚的ノイズなど、がある。この明らかに睡眠脳波とは異なる不規則データが混入する結果、脳波データには例外値（アウトライア）、異常値、不規則値がしばしば不可避的に含まれる。従って睡眠ステージの判定を精度良く行おうとすれば、そのような例外値などの不規則値を正しく特定、検出し、除去することが不可欠となる。

また、例えば携帯型１チャンネル脳波計を用いた計測の場合、従来の１０−２０法で指定される電極位置とは必ずしも一致しない位置に電極を配置する場合がある。例えば、額、耳のように髪の毛がなく皮膚が露出しているような場所であって、被験者が自ら電極をこれら位置に装着して睡眠脳波を計測しようと試みる場合にも、容易に電極の配置が可能になるような場所である。このような条件下で睡眠脳波から自動的に睡眠ステージを推定しようとすると、電極位置ずれや被験者の個性に起因する睡眠脳波周波数パワー強度の違いが、判定精度を劣化させる原因となる。従って携帯型脳波計の場合、このような電極の位置ずれや被験者の個性によるバラツキも考慮に入れる必要がある。

従って、本発明の課題は上記の問題点を解決することによって、容易に、被験者の個性等に起因する脳波データの変動を捕捉し、脳波データに混入する例外値を除去し、携帯型電極の場合の電極位置ずれによる不規則値などを検出し、補正することができるとともに精度良く、手動によらず完全自動で判定が可能な睡眠ステージ自動判定のための装置と方法およびそのためのコンピュータプログラムを提供することである。

前記課題を解決するための本発明は、第１の側面によれば、
脳波電極を通じて計測された脳波データを入力して、睡眠ステージをエポック単位の時系列で判定して出力するための睡眠ステージ判定装置であって、
前記入力された脳波データの周波数成分をエポック単位で解析して周波数パワー値を算出し、その周波数ベクトルを出力する周波数解析手段、
前記出力された周波数ベクトルの中から主成分を抽出して、固有値累積寄与率が９０％以上になるように３〜５次元の特徴ベクトルに縮約する主成分分析手段、
前記３〜５次元特徴ベクトルを入力して、その各成分に対して補正・規格化を行って不規則値を除去するための規格化手段、
および
前記正規化後の３〜５次元特徴ベクトルを入力して、覚醒、軽い睡眠、深い睡眠およびレム睡眠の少なくとも４クラスの中から、各エポックにつき１クラスを出力する多クラス識別手段を備えている睡眠ステージ自動判定装置に関する。

第２の側面によれば、
睡眠脳波信号から睡眠ステージを自動判定する方法であって、エポック単位の時系列順に次のステップ：
（１）脳波信号を周波数解析して、周波数パワー値を算出し、周波数ベクトルを出力するステップ、
（２）前記出力された周波数ベクトルを主成分分析して、３〜５次元特徴ベクトルに縮約するステップ、
（３）前記特徴ベクトルの各成分から不規則値を除去して補正・規格化するステップ、
（４）正規化された前記特徴ベクトルを多クラス識別手段に入力して、睡眠ステージとしての、覚醒、軽い睡眠、深い睡眠およびレム睡眠の少なくとも４クラスの中から１クラスを識別、出力するステップ、
を含む睡眠ステージの自動判定方法に関する。

さらに第３の側面によれば、本発明は、コンピュータに睡眠ステージの自動判定機能を付与するプログラムであって、脳波計により被験者の睡眠時脳波を計測して得られた、デジタル脳波信号のエポックごとの時系列データを時系列処理する次の工程：
（１）脳波信号を周波数解析手段に入力して、周波数パワー値を算出し、周波数ベクトルを出力する工程、
（２）前記出力された周波数ベクトルを主成分分析手段に入力して、３〜５次元特徴ベクトルに縮約する工程、
（３）得られた前記特徴ベクトルを百分位数統計手段または重み係数を付与した非線形ワンクラス・サポートベクトルマシンによって規格化する工程、
（４）正規化された前記特徴ベクトルを多クラス識別手段に入力して、睡眠ステージとしての、覚醒、軽い睡眠、深い睡眠およびレム睡眠の少なくとも４クラスの中から１クラスを識別、出力する工程、
を含む睡眠ステージの自動判定プログラムに関する。

上述したように、本発明による睡眠ステージ自動判定装置と方法およびそのためのプログラムは、例えばサポートベクトルマシンによる多クラス識別アルゴリズムを使用することおよびそのために適切な規格化手段として百分位数統計法または重み係数を付与したワンクラス・サポートベクトルマシン法を適用することを特徴とするが、このような正規化法の適用が脳波データに特有の統計的特徴、特性にとって最適かつ有効であるとの知見に基づいている。

そこで、図３に示すような典型的睡眠脳波データの統計的特徴からその理由について説明する。
一般に、睡眠ステージは３つのクラスに大別される、すなわちＮＯＮＲＥＭ（ノンレム睡眠）３０１、ＲＥＭ（レム睡眠）３０３およびＷＡＫＥ（覚醒）３０５である。レム睡眠とノンレム睡眠のデータは二峰性の山のクラスタ（３０４と３０２）に含まれ、夫々楕円状に分布しているのに対して、覚醒時データは被験者のその日の行動に依存し、偶発的にランダムに散布（３０６）している。
睡眠脳波の周波数帯域の感度に由来するバラツキは、おもに分布密度の平行移動を引き起こす。従って、被験者、測定ごとのバラツキを補正するには、適切に設けた基準点が夫々のデータで一致するように分布密度の平行移動を施す必要がある。
ところで、夫々のクラスタに属するデータの割合は、典型的には、覚醒２５％：ノンレム睡眠７０％：レム睡眠５％である（この比率は各クラスタの山の高さに比例）とされる。しかしながら、この割合はその日の睡眠状態に依存して変わり、或いは不眠症のような睡眠疾患の場合も同様である。従って全ての睡眠ステージを含む睡眠の特徴ベクトルの分布の重心を基準点にすることは出来ない（望ましくない）。
ノンレム睡眠は睡眠の深さに対してデルタ周波数帯域のパワーと正の相関を持つことが知られており、深い睡眠であればあるほどデルタパワーは原点に対して非対称に増大する。ノンレム睡眠に対応するデータ分布の広がりは日々の睡眠状態に依存する。そのため、ノンレム睡眠の分布を特徴づけるのに、例えば中央値などの指標は、被験者の個性や日々の眠りの具合によって大きく変動を受けるので、使うことができない。
他方、レム睡眠に対応する分布は睡眠中のレム頻度に応じて、分布の高さは変化するものの、分布形状が大きくは変化しないので、分布から得られる中央値などの指標は日々の睡眠深度に対して、比較的安定的に位置する特徴がある。従ってレム睡眠の分布に重きをおいた統計量を基準にすることによって安定的に睡眠状態を評価することが可能である。
以上要するに、脳波データのバラツキ補正にはレム睡眠データの扱いが鍵となるのであって、その処理方法として、レム睡眠の分布に重きをおいた統計量を基準点にすること、あるいはレム睡眠の分布に重心がくるように基準点を設計することにより、例外値を排除した後に抽出したレム睡眠クラスタの平均統計量を計算することによって、確実に睡眠状態を評価することが可能となる。

上記特徴を有する本発明の睡眠ステージ自動判定装置と方法およびそのためのプログラムは、短時間で、人手を介することなく完全自動で、睡眠状態の自動判定を精度良く行うことが可能である。また、本発明に係る睡眠ステージ自動判定は携帯型簡易脳波計から得られる脳波データのみならず、ＰＳＧから得られる脳波データ（ＥＥＧ）に対しても適用可能である。

図１は本発明に従う睡眠ステージ自動判定装置の一例を示す概要図である。 図２は睡眠ステージ自動判定方法の一例を示すフローチャートである。 図３は典型的な睡眠脳波の波形と主成分空間上のデータ散布図を示す説明図である。 図４の（A）は主成分分析後、規格化手段適用前の被験者間のバラツキの様子を表すデータ分布図、（B）は規格化手段適用後のバラツキ改善の様子を表すデータ分布図である。 図５は百分位数統計法による規格化において、パーセンタイル・パラメータを変化させた場合の正解率の変化を示すグラフである。 図６はワンクラス・サポートベクトルマシン（重み係数の導入なし）による規格化において、アウトライアを検出、除去する場合の、インライアとアウトライアの分布図である。 図７は図６において、重み係数を導入したワンクラス・サポートベクトルマシンによる規格化後のインライアとアウトライアの分布図である。 図８は本発明の睡眠ステージ自動判定装置で使用される電極と脳波計の１例であって、装着状態を示す概略図である。

以下、図面を参照して本発明の睡眠ステージ自動判定装置と方法の概要について説明する。
本発明の睡眠ステージ自動判定装置３は、図１に示すように、周波数解析器３１、主成分分析器３２、規格化器３３、多クラス識別器３４を主たる構成要素とする。その外、必要に応じて、脳波検出電極１；脳波計２；アーチファクト検出器３５；ノイズ除去・平滑化を行う後処理装置４；判定結果をビジュアル表示する表示装置５から構成されている。

図１において、脳波検出電極１と脳波計（脳波記録計）２とは脳波検出部を構成している。脳波計２は増幅器／ＡＤ変換器２１およびデータ保存器２２を備えており、判定装置３とは無線または有線で接続することができる。
脳波検出電極１としては、例えば図８に図示するように、被験者の額、耳裏部分に装着する携帯型電極８１、８２を使用することができる。両電極８１，８２はコード８３を介して脳波計２と繋がっている。また、ＰＳＧ検査で使用される導電性ペーストを皮膚に塗布して装着する皿電極タイプのものも使用できる。携帯型脳波電極としては布電極を用いた乾式電極（実用新案登録第３１５３４０９号参照）、耳電極（特開２００８−６７９１１号参照）などが挙げられる。脳波検出電極１によって収集された脳波データは増幅器／ＡＤ変換器２１とデータ保存器２２を備えた脳波計２で記録される。このデータ保存器２２は脳波計１の中に含まれていてもよく、或いは無線の送受信機を備えることによって外部記憶装置に保存されるような構成であってもよい。また、脳波データの記憶装置を設けず、睡眠ステージ判定装置３に伝送して、計測と同時的に直接、判定を行うようにすることもできる。

睡眠ステージ自動判定装置３は、伝送された前記脳波データを解析処理するための装置であって、前記脳波計２とは無線または有線で接続している。そして汎用コンピュータに判定プログラムの形で内蔵されるか、またはプログラム化された専用コンピュータの形で存在する。前記構成になる睡眠ステージ自動判定装置３では、図２に示すフローチャートに従って、判定処理工程が自動的に進行する。
すなわち、図２において、エポック単位（例えば３０秒単位のセグメント）で計測された時系列脳波信号は、時系列で睡眠ステージ自動判定装置に入力され、先ず周波数解析器３１に送られて、例えば高速フーリエ解析（ＦＦＴ）またはウエーブレット解析によって周波数パワー値が計算される（Ｓ２０）。計算された周波数パワー値は夫々の周波数成分ごとにゼロ点補正を行った後、主成分分析工程に送られる（Ｓ２１）。その過程で主成分分析器３２にもとづく主成分を表す部分空間へ射影し、非主成分であるノイズ成分が除去される。主成分空間へ射影されたデータ値は、さらに、規格化器（正規化器）３３によって成分ごとのゼロ点補正が行われた（Ｓ２２）後、多クラス識別器３４に入力されて（Ｓ２３）、最終的に睡眠ステージ（W：覚醒、R：レム睡眠、L：軽い睡眠、D：深い睡眠）の少なくとも４種類のラベルが出力される。

以下、各ステップの詳細について述べる。
（１）周波数解析工程
例えば３０秒単位の各エポックの脳波データはさらに５秒単位の６個のサブ・エポックに細分割して、サブ・エポック単位でＦＦＴ等によって周波数解析され、周波数ごとにパワー強度が計算される。６個のサブ・エポックの周波数パワー強度は平均化されて、各エポックの周波数パワー強度として算出される。パワー値は対数に換算して、次工程の特徴量として用いる。
携帯型脳波計の場合、サンプリング周波数が１２８Ｈｚであるので、最大６４Ｈｚまでの周波数に対応するパワー値を得ることはできるが、取り扱うデータ量が大きくなりすぎるので、脳波にあまり関係のないと思われる帯域を除外して、３０Ｈｚ以下の周波数成分が考慮される。５秒単位でＦＦＴ計算を行うため、周波数分解能は０．２Ｈｚ（１／５ｓｅｃ）となる。従って、３０Ｈｚまでの周波数成分は１５０（３０／０．２）の周波数点からなる。言い換えれば、１エポックの脳波データはＦＦＴ計算後、１５０次元の周波数ベクトルとして表現される。

（２）主成分分析工程
周波数ベクトルの出力後、主成分分析工程前に、主成分分析を行うためのゼロ点補正が行われる。このために、本発明では平均値でなく、百分位数（パーセンタイル）統計手段を用いて、２０〜８０パーセンタイルを基準としてゼロ点補正を行うことができる。百分位数統計法については（３）規格化工程で後述する。
ゼロ点補正後のデータに対して、主成分分析器３２により主成分ベクトルを抽出し、特徴ベクトルの固有値が大きな順番から固有値累積寄与率が９０％以上になるまで主成分ベクトルを選び出し、次元縮約が行われる。通常、おおよそ３〜５個の特徴ベクトルが選ばれる。主成分分析器３２はこれら固有ベクトルによって張られる部分空間に射影したときの座標を出力する。

（３）規格化工程 方法１
縮約された３〜５個の特徴ベクトルは、１つの方法として、その各成分について、百分位数統計量を用いてデータの１〜１０パーセンタイル値を基準として規格化することができる。
ある周波数成分に対する百分位数統計量は、その周波数特徴量の成分を小さい順（昇順）に並べ替え、先頭からｑ番目に相当する値をｑ百分位数の値とすることで計算される。この値をその周波数成分の時系列データから引き算することによりゼロ点補正が行われる。百分位数による正規化は数式１によって行われる。
本発明では、特徴ベクトルの各成分がデータの１〜１０パーセンタイル値、好ましくは１〜５パーセンタイル値を基準として調整される。
百分位数による正規化手段を用いることによって、睡眠脳波における被験者間のバラツキを吸収することができ、後続の多クラス識別手段での睡眠ステージ判定が容易になる。

（３）規格化工程 方法２
本発明では、選択的に適用できる別の規格化法として、重み付けをしたワンクラス・サポートベクトルマシン（以下、ワンクラスＳＶＭと略記）手段を用いて、正規化することもできる。

その前に、一般的ワンクラスＳＶＭの理論的側面について述べる（Scholkopf
et al, Neural Computation, 13, pp. 1443-1471 (2001); “Estimating the Support of
a High-Dimensional Distribution”参照）。
ワンクラスＳＶＭとは、数式２に示されるような二値化関数によって例外値（アウトライア）とそうでないものとを弁別する識別器である。最適な変数ρおよびαはコスト関数（目的関数）である数式２の二次計画問題を解くことによって得られる。例外値とそうでないものとを分ける識別面は変数νによって決まり、この変数は超空間内での百分位数に相当する値である。

以下、事前確率分布を導入する上で、ワンクラスＳＶＭの変更点について述べると、数式２で定められたコスト関数を最小化することによって、同様に得ることができる。
ワンクラスＳＶＭにおいては変数ｗ、ξおよび閾値ρがこのコスト関数を低減するように最適化される。変数νはインライア（想定される誤差範囲内の測定値）とアウトライア（外れ値）の割合（比）を決める意味をもつ変数である。これらの変数は線形の不等式で表される制約条件の下で解かなければならない。得られた変数の値を基に計算したｗによって計算される、αを重み係数とする非線形のカーネル関数を用いた二値化プロセスによって、インライアとアウトライアの判別が数式３によって行われる（sgnは符号関数を表す）。符号関数は引数が正の場合に＋１を、負の場合に−１を出力する。カーネル関数にはガウス関数のような非線形関数を用いることができる。

しかしながら、単純にワンクラスＳＶＭを適用しただけでは、本発明の睡眠状態解析において目的とするレム睡眠に重きを置いてクラスタ重心を計算できないので、ワンクラスＳＶＭにおける二次計画問題に対して、データの重み係数ωを導入することで、この重み係数をレム・クラスタになりそうなデータに対して多く配置することによって、そのデータの抽出を可能にしている。
下記の数式４は、数式２に対して重み係数ωを導入したときの二次計画問題を表している。重み係数は総和が１になるように規格化されており、例えば数式５のような指数関数を用いることができる。
重み係数を含む二次計画問題の解法には、データの重み係数に応じたリサンプリング処理などの方法を用いて簡便に計算することもできる。

（４）多クラス識別工程（Ｓ２３）
多クラス識別器３４には、出力するクラス数に対応する複数個のサポートベクトルマシーン（ＳＶＭ）が用いられる。夫々のＳＶＭには規格化器３３から出力された正規化後の３〜５次元の特徴ベクトルが入力される。多クラス識別器３４では、少なくとも４クラス、すなわち覚醒（Ｗ）、レム睡眠（Ｒ）、浅い睡眠（Ｌ）、深い睡眠（Ｄ）の睡眠ステージが出力されるが、ＳＶＭは本質的に２クラス識別手段であるので、４クラスの識別を行うには、２クラスずつの組み合わせごとにＳＶＭを構築する。すなわち、Ｗ対Ｒ、Ｒ対Ｌ、Ｌ対Ｄ、などのすべての組み合わせ（都合６組）のペアのＳＶＭを用意する。６個のＳＶＭは識別結果（例えば、Ｗ対Ｒを判定するＳＶＭの場合、ＷかＲ）に基づいて該当するクラスに１票の投票が行われる。多クラス識別器はこの投票結果から最も多く投票されたクラスを最終的な判定結果とする。
また、時間的な連続性を考慮するために、適当な時間窓を設定して（例えば前後７エポック）、この間に含まれるすべてのＳＶＭの結果を投票する。

以上の工程で睡眠ステージ自動判定プロセスが終了するが、必要に応じて、前記自動判定装置３に接続して、さらに時間的な睡眠ステージの遷移を考慮してノイズ除去と平滑化を行う後処理装置４が設けられる。この装置によれば、例えば隠れマルコフモデル手法による状態遷移確率から睡眠ステージの遷移を推定することによって平滑化を行うことができる。
さらに、判定後の睡眠解析結果を視覚化して表示する表示装置５を後処理装置４に接続して設けることもできる。表示装置に表される睡眠解析結果の項目には、例えば睡眠経過図（睡眠プロファイル）、総睡眠時間、入眠潜時、各睡眠ステージの時間数と睡眠全体に占める割合、覚醒指数、睡眠効率などの様々な数値、評価の外に当該被験者へのアドヴァイスなどの情報が含まれる。

［睡眠脳波の計測実験例］
以下、本発明の睡眠ステージ自動判定方法の精度、有効性を確認するために、睡眠実験を行った。この実験は、睡眠クリニック内で、睡眠に適した統制された環境下で行った。
被験者には３０歳〜５０歳の５名の健康な男性が選ばれた。被験者には、自らが受ける睡眠検査が標準的なＰＳＧ検査と同時に行われる携帯型脳波計による測定実験であることを事前に説明し、承諾を得ていた。
被験者は２０時に入床し、ベッドの上で、安静仰臥姿勢で寝るまでの間待機していた。そして所定の位置に携帯型脳波計を装着した。２２時に消灯され、携帯型脳波計による睡眠脳波測定の睡眠実験が開始された。翌朝６時頃、起床するまで脳波計測が行われた。

実験に使用した、生体電気信号計測用電極（日本光電製）は額部と耳後の２箇所に配置され、両者の電位差が１チャネル携帯型脳波計（ＯＢＩ製）によって計測された。脳波信号は低域通過特性を持つアナログ増幅器によって増幅され、１０ビットのＡＤコンバータによって毎秒１２８点でサンプリングされ、デジタル信号に変換された。デジタル信号に変換された信号は脳波計内の磁気記録媒体（ｍｉｃｒｏ−ＳＤ）によってリアルタイムに保存された。
このデータを使っての実施例を以下に述べる。

一晩の計測実験の終了後、保存された各被験者のデータは自動判定装置を内蔵するパソコンに転送され、そこで上述の判定手順による、周波数解析、主成分分析、規格化、多クラス識別の各工程を経て睡眠ステージの解析、判定が自動的に行われた。
規格化工程では、百分位数統計手段を適用した。

（i）周波数解析工程では、３０秒単位の各エポックの脳波データをさらに５秒単位のサブ・エポックに細分割して、サブ・エポック単位でＦＦＴによって周波数解析し、各周波数ごとにパワー強度が計算された。６個のサブ・エポックの周波数パワー強度は平均化し、各エポックの周波数パワー強度として算出した。パワー値は対数変換して、次工程の特徴量として用いた。
サンプリング周波数が１２８Ｈｚであるので、最大６４Ｈｚまでの周波数に対応するパワー値を得ることはできるが、脳波にあまり関係のないと思われる帯域を除外して、３０Ｈｚ以下の周波数成分を考慮した。５秒単位でＦＦＴ計算を行い、周波数分解能は０．２Ｈｚ（１／５ｓｅｃ）であるので、３０Ｈｚまでの周波数成分は１５０（３０／０．２）の周波数点からなることになる。従って、１エポックの脳波データはＦＦＴ計算後、１５０次元の周波数ベクトルとして出力された。
（ii）主成分分析工程 １５０次元の周波数ベクトルを主成分分析器２２に入力して、主成分ベクトルを抽出し、特徴ベクトルの固有値が大きな順番から固有値累積寄与率が９０％以上になるまで次元縮約し、その結果、３〜５次元の特徴ベクトルを出力した。
（iii）規格化工程では、得られた３〜５次元特徴ベクトルの各成分について、百分位数統計手段によって、１〜１０パーセンタイルの範囲から１パーセンタイルを選びその値を基準にして規格化を行った。
（iv）多クラス識別工程については、上記段落番号００２６に記載した通りである。

前記（iii）規格化工程の奏する効果について以下述べる。
主成分分析後で規格化器適用前のデータは図４（Ａ）のデータ散布図に示す通りであり、規格化器適用後のデータは図４（Ｂ）に示す通りである。各図において、各ドットは夫々のエポックにおけるデータの位置を表わし、形の異なるドットは被験者の違いを表わしている。横軸は第１主成分、縦軸は第２主成分を表わす。各ドットはエポックにおける周波数特徴ベクトルの主成分座標上の点を表わしている。一夜の計測では約６時間前後のデータが計測され、３０秒を１単位とするエポックごとに周波数特徴ベクトルが計算されたので、この図の中には夫々の被験者について６×６０×（６０／３０）＝７２０点がプロットされている。

規格化器適用前の図４（Ａ）では被験者ごとのバラツキが大きいが、これは夫々の被験者ごとに睡眠脳波を構成する周波数要素の強度が異なるために、夫々の測定データがもつ主軸とすべてのデータを考慮したデータの主成分の主軸とが必ずしも一致しないことに起因すると考えられる。
図４（Ａ）の横軸の第１主成分はおおよそすべての周波数成分の総和に対応したものであり、覚醒、体動などによって生じる大きな揺らぎに対応しており、他方縦軸の第２主成分はおおよそデルタ帯域強度とベータ帯域強度との差を表わしていると考えられる。図では、正の方向に向かってデルタ帯域強度が強くなり、負の方向に向かってベータ帯域強度が強くなる。

規格化器適用後のデータ分布を表す図４（Ｂ）の第１主成分および第２主成分の意味合いは、図４（Ａ）でのように単純ではないが、図の右上方向が、デルタパワーが増加する軸に対応し、図の左上方向が全周波数の総和に比例して増加する軸に対応すると考えられる。
図４（Ｂ）からは、データの分布が顕著に二峰性のクラスタを呈することが明らかである。即ち、上方に位置するクラスタはノンレム睡眠状態に対応しており、下方に位置するクラスタはレム睡眠状態に対応している。下方に位置するクラスタ（レム状態に対応する）は重なり整頓されているため、被験者のバラツキが小さくなっていることが分かる。

図５は主成分分析後に適用された規格化器のパーセンタイル・パラメータを変化させたときの正解率の変化を表わしたものであり、パーセンタイルが０に近づくほど正解率の向上が見られるが、１〜１０パーセンタイルの範囲で高い正解率が得られることが分かる。０パーセンタイルになる点はこのデータの例ではレム状態を示すクラスタに対応する。すなわち、上記図４（Ｂ）に示したように、レム・クラスタに重きをおくことで被験者のバラツキを低減することができることを意味している。

レム・クラスタの中心で揃えることが好適であるのは以下の理由による。
ノンレム睡眠とレム睡眠は一般に平均９０分サイクルで繰り返し、その割合はほぼ７：３であることが知られている。しかしながら、実際の個々の睡眠においてはこの知見はあてにならない。例えば、睡眠時間が短くノンレム睡眠、レム睡眠周期が終わる前に起床したような場合や、そもそも睡眠リズムを正確なリズムで繰り返している例は実際にはあまりない。このため、ノンレム睡眠とレム睡眠の割合は、個々の睡眠データによって異なる。

一方、分布の平均値はデータに含まれるノンレム睡眠とレム睡眠の割合によって変動するため、データごとに分布中心が変化する。最も出現頻度の高い睡眠ステージ２の状態はデルタ周波数パワー値でソートしたとき中ごろに位置する。従って、デルタ周波数パワー値の中央値を、分布中心とみなせば、単純に平均値を用いるよりは、被験者間のバラツキをよりよく補正することができる。

しかしながら、ノンレム睡眠は睡眠の深さによってデルタ周波数帯域のパワーと正の相関を持っているため、ノンレム睡眠の分布は深い睡眠であればあるほど、より広い分布になる。すなわち、睡眠深度はそのときの睡眠の状態に依存するため、中央値は被験者の個性や日々の眠りの具合によって大きく変動を受ける。従って、分布中ごろに位置する基準点は毎回変化するため適切でない。レム中心は睡眠深度のパワーに対して最小点に近い場所に位置するために、日々の睡眠深度のパワー変化に対してよりロバストに基準点として用いることができる。

また、最小点は計測が行われず、データがゼロ値を含む場合に対して誤作動する可能性がある。したがって、実用的には１〜１０パーセンタイル値、好ましくは１〜５パーセンタイル値程度に設定することが望ましいのはそのためである。

前記データの正解率の評価は、１つ取って置き（ｌｅａｖｅ−ｏｎｅ−ｏｕｔ）交差検証法によって、５名の被験者のデータについて行った。すなわち、５名の被験者の中から４名を学習データとして選択して識別手段を学習する。学習された識別手段に対して、残る１名の被験者データをテストデータとして用い、テストデータに付随する睡眠技師のステージ判定結果と、識別手段が出力する判定結果の睡眠クラスとが一致する数を計測する。そしてこの数をテストデータの個数で除算し、識別性能として評価した。

ワンクラスＳＶＭ手段、特に重み係数を付与したワンクラスＳＶＭ手段による正規化工程の例である。それ以外の工程は実施例１と全く同様である。
前記特徴ベクトルデータに対して重み係数を付与し、全データのうち２０％のインライアデータが検出されるように調整した（すなわち閾値ν＝２０％に設定）非線形ワンクラスＳＶＭを適用してインライアを検出し、検出したインライアデータの平均値を算出し、得られた平均値を基準点として調整した。

図６および図７は夫々、重み係数を適用しないで、単純にワンクラスＳＶＭを適用した場合および重み係数を適用したワンクラスＳＶＭ適用の場合の夫々のアウトライア検出結果を表すデータである。
図６中、黒◆はインライアデータを表し、白△はアウトライアデータを表す。点Ａはインライアデータの平均値を計算した場合の中心点を示し、点Ｂは全ての点（インライア＋アウトライア）の平均値を計算した場合の中心点を示す。単純加算平均した点である、点Ｂは中途覚醒、寝返りなどによる体動を含むデータの影響を受けるため、データが集中している範囲からずれている。
対照的に、図７ではデータの重み係数ωを導入して、この重み係数をレム・クラスタに対応するデータに対して多く配置するようにしたワンクラスＳＶＭを適用することによって、レム睡眠に対応するクラスタ（中心点Ｃ）を抽出できることが分かる。

推定された睡眠ステージの時間的な状態遷移の妥当性を考慮に入れるために、識別器の出力結果をさらに隠れマルコフモデル（ＨＭＭ）の入力としてノイズの除去と平滑化を行う。ＨＭＭでは数式６で示されるコスト関数を最大化するように睡眠状態が決められる。
を最小化するように状態が求められる。ここでＨＭＭの隠れ変数s_i ^tは４つの状態を仮定し、確率p(s_i ^t,
y_j ^t)は前段の処理結果y_j ^tと隠れ変数s_i ^tとの同時確率である。p(s_i ^t,s_j ^t+1)は現時点tでの隠れ状態iと次時点t+1での隠れ状態jへの遷移確率を表しており、ここでは、簡便的な方法として前処理結果p(y_i ^t,y_j ^t+1)を代用する。

ＨＭＭの適用によって、前段の識別器が、生理学的にあまりみられない状態遷移、例えば、深い睡眠から覚醒へ遷移する時間変化を示す場合には、ＨＭＭで定義されるコスト関数は小さくなるので、ＨＭＭによって出力される睡眠状態は、深い睡眠から浅い睡眠を経て覚醒状態へ遷移するように修正される。このように、ＨＭＭの適用によって、時間経過を考慮するノイズの除去と平滑化を行うことができる。

本発明に関わる睡眠ステージ自動判定装置、睡眠ステージ自動判定方法、そのためのプログラムによれば、睡眠脳波の知識がない一般人でも睡眠ステージ自動判定方法を組み込んだシステムを操作するだけで（脳波計を装着してスタートさせるだけで）、判定結果を得ることが可能となり、在宅での健康促進プログラムに使用したり、睡眠疾患を高速で検査するための簡易ツールとして利用することが可能である。

１ 脳波検出電極
２ 脳波計
３ 睡眠ステージ自動判定装置
４ 後処理装置
５ 表示装置
２１ 増幅器／ＡＤ変換器
２２ データ保存器
３１ 周波数解析器
３２ 主成分分析器
３３ 規格化器
３４ 多クラス識別器
３５ アーチファクト検出器

Claims (14)

1. 脳波電極を通じて計測された脳波データを入力して、睡眠ステージをエポック単位の時系列で判定して出力するための睡眠ステージ判定装置であって、
   前記入力された脳波データの周波数成分をエポック単位で解析して周波数パワー値を算出し、その周波数ベクトルを出力する周波数解析手段、
   前記出力された周波数ベクトルの中から主成分を抽出して、固有値累積寄与率が９０％以上になるように３〜５次元の特徴ベクトルに縮約する主成分分析手段、
   前記３〜５次元特徴ベクトルを入力して、その各成分に対して補正・規格化を行って不規則値を除去するための規格化手段、
   および
   前記規格化後の３〜５次元特徴ベクトルを入力して、覚醒、軽い睡眠、深い睡眠およびレム睡眠の少なくとも４クラスの中から、各エポックにつき１クラスを出力する多クラス識別手段を備えていることを特徴とする睡眠ステージ自動判定装置。
2. 前記規格化手段は百分位数統計手段を含むことを特徴とする請求項１に記載の自動判定装置。
3. 前記規格化手段は重み付けした非線形ワンクラス・サポートベクトルマシンを含むことを特徴とする請求項１に記載の自動判定装置。
4. 前記睡眠ステージ自動判定装置で得られた判定結果を視覚化して表す表示装置をさらに含む、請求項１〜３のいずれか１項に記載の自動判定装置。
5. 脳波電極および脳波計をさらに含み、前記脳波電極は被験者の額部・耳周辺部皮膚に接触して装着される電極であることを特徴とする、請求項１〜３のいずれか１項に記載の自動判定装置。
6. 前記睡眠ステージ自動判定装置に接続して、ノイズ除去と平滑化を行うための後処理装置を含む、請求項１〜３のいずれか１項に記載の自動判定装置。
7. 前記後処理装置として、隠れマルコフモデル手法による状態遷移確率から睡眠ステージの遷移（経路）を推定することによって、ノイズ除去と平滑化を行う、請求項６のいずれか１項に記載の自動判定装置。
8. アーチファクトを検出、除去するための手段をさらに含む請求項１〜３のいずれか１項に記載の自動判定装置。
9. 前記多クラス識別手段は少なくとも６個のサポートベクトルマシンを含むことを特徴とする請求項１〜３のいずれか１項に記載の自動判定装置。
10. 睡眠脳波信号から睡眠ステージを自動判定する方法であって、エポック単位の時系列順に次のステップ：
    （１）脳波信号を周波数解析して、周波数パワー値を算出し、周波数ベクトルを出力するステップ、
    （２）前記出力された周波数ベクトルを主成分分析して、３〜５次元特徴ベクトルに縮約するステップ、
    （３）前記特徴ベクトルの各成分から不規則値を除去して補正・規格化するステップ、
    （４）規格化された前記特徴ベクトルを多クラス識別手段に入力して、睡眠ステージとしての、覚醒、軽い睡眠、深い睡眠およびレム睡眠の少なくとも４クラスの中から１クラスを識別、出力するステップ、
    を含む睡眠ステージの自動判定方法。
11. 前記規格化ステップ（３）は前記特徴ベクトルの各成分を百分位数統計手段に入力して、各成分の１〜１０パーセンタイル値が基準となるように調整することを特徴とする、請求項１０に記載の睡眠ステージの自動判定方法。
12. 前記規格化ステップ（３）は、前記特徴ベクトルのデータのうち最も重要なクラスタに最も多量の重み係数を配分しかつ全データの２０〜４０％のインライアデータが抽出されるように調整した、重み付け非線形ワンクラス・サポートベクトルマシンに対し前記特徴ベクトルデータを入力し、その結果検出されたインライアデータの平均値を計算して、基準とすることを特徴とする、請求項１０に記載の睡眠ステージの自動判定方法。
13. コンピュータに睡眠ステージの自動判定機能を付与するプログラムであって、
    脳波計により被験者の睡眠時脳波を計測して得られた、デジタル脳波信号の時系列データを各エポック単位で、次の順に順次的に処理するステップ：
    （１）脳波信号を周波数解析手段に入力して、周波数パワー値を算出し、周波数ベクトルを出力するステップ、
    （２）前記出力された周波数ベクトルを主成分分析手段に入力して、３〜５次元特徴ベクトルに縮約するステップ、
    （３）得られた前記特徴ベクトルを百分位数統計手段または重み係数を付与した非線形ワンクラス・サポートベクトルマシンによって規格化するステップ、
    （４）規格化された前記特徴ベクトルを多クラス識別手段に入力して、睡眠ステージとしての、覚醒、軽い睡眠、深い睡眠およびレム睡眠の少なくとも４クラスの中から１クラスを識別、出力するステップ、
    を含むことを特徴とする睡眠ステージの自動判定プログラム。
14. 請求項１３に記載の前記睡眠ステージの自動判定プログラムを格納したことを特徴とするコンピュータ読み取り可能な記録媒体。