説明

信号解析法

持続時間制約を用いることによって、ECG信号のような信号のセグメンテーションの信頼度を改善することが開示される。その信号は、隠れマルコフモデルを用いて解析される。持続時間制約は、そのモデルの特定の状態のための最小許容持続時間を指定する。持続時間制約は、モデルそのものの中に組み込むことができるか、又は従来のモデルを与えられるときの最も確からしい状態シーケンスを計算するために用いられるビタビアルゴリズム内に組み込むことができる。また、そのモデルから信頼度指標を導出することも開示され、その信頼度指標を用いて、セグメンテーションの品質及び頑強性を評価することができ、さらに、たとえば、信号内に雑音又は異常が存在することに起因して、セグメンテーションの信頼度がない信号を特定することができる。

【発明の詳細な説明】
【技術分野】
【0001】
本発明は信号を解析するための方法に関する。その解析は、持続時間制約を組み込む隠れマルコフモデリング技法を利用して、心電図(ECG)又は他の生物医学信号のような信号を、その要素波形特徴に自動的にセグメンテーションすることを含むことができ、任意選択で、それと合わせて、セグメンテーションの精度及び頑強性を指示する、関連する信頼度指標を導出することを含む。
【0002】
ECG(頭字語EKGによっても知られている)は、心臓の電気的活動を測定する重要な非侵襲性の信号である。個々の心拍は多数の異なる心臓学的段階から構成され、それらの段階は、ECG波形において1組の異なる特徴を生み出す。これらの特徴は、心臓の特定の領域における筋細胞の消極(放電)又は再分極(再充電)のいずれかを表す。図1は、人のECG波形及び関連する特徴を示す。ECG波形の標準的な特徴は、P波、QRS群及びT波である。さらに、小さなU波が(T波の後に)時折存在する。
【0003】
心周期はP波(その開始点及び終了点はPon及びPoffと呼ばれる)で開始し、それは心臓の心房消極(atrial depolarization)の周期に対応する。この後にQRS群が続いており、QRS群は一般的に、ECG波形の最もはっきりとした特徴であり、心室消極(ventricular depolarization)の周期に対応する。QRS群の開始点及び終了点はQ点及びJ点と呼ばれる。QRS群の後にT波が続いており、それは心室再分極の周期に対応する。T波の終了点はToffと呼ばれ、心周期の終了を表す(U波が存在しないものと仮定する場合)。
【0004】
ECG信号を詳細に調べることによって、特徴的なECG波形から多数の有益な測定値を導出することができる。その後、これらの測定値を用いて、患者の病状を評価することができ、心調律(cardiac rhythm)内にある任意の潜在的な異常を検出することができる。
【0005】
特に重要な測定値は「QT間隔」であり、それはQRS群の開始からT波の終了までの時間、すなわちToff−Qと定義される。この時間間隔は、心室内の電気的活動(消極及び再分極の両方)の総持続時間に対応する。
【0006】
QT間隔は、通常の間隔よりも長くなると、QT延長症候群(LQTS)の有効な指標となるので特に重要である。これは、顕在はしていないが、致命的な状態であり、その患者は、「トルサード・ド・ポワント(torsade de pointes)」として知られている非常に速く、異常な心調律(不整脈)に陥りやすくなる。このリズムが生じるとき、心臓は有効に拍動できなくなり、脳への血流が大幅に低下する。その結果、意識を突然失い、心臓死にいたる可能性もある。
【0007】
別の重要な指標はPR間隔であり、それは、P波の開始からQRS群の開始までの時間、すなわちQ−Ponと定義される。これは、心房消極の開始から心室消極の開始までの時間に対応し、房室心ブロックのような状態を検出するために用いることができる。
【0008】
QT間隔及びPR間隔の測定及び評価に加えて、その要素特徴へのECG波形の正確且つ頑強なセグメンテーションも、ST上昇のような状態を検出するために、さらには、連続した心拍内のQRS群のピーク(R点として知られている)の位置を特定することによって心拍数を計算するために重要である。
【0009】
それゆえ、ECG波形において特徴境界を正確に求めることに関する上記の重要性を踏まえて、この過程を自動化するために多大な努力がなされてきた。この分野における初期の仕事は、閾値処理及びテンプレートマッチングのような古典的な信号処理及びパターン認識手法に焦点を合わせていた。最近になって、隠れマルコフモデル(HMM)に基づく統計的な手法が、ECGセグメンテーションの問題に適用されている。
【0010】
隠れマルコフモデリングでは、根底を成す、対象となる状態シーケンスが存在するが、そのシーケンスは直に観測することはできない。一方では、さらに、この状態シーケンスに確率的に関連付けられる信号又は観測シーケンスも存在し、そこから、所与の状態シーケンスの最も確からしい値を推測することが望ましい。
【0011】
隠れマルコフモデルの中心には、2つの確率関数があり、そのうちの一方は、特定の時刻ステップの直前の時刻ステップにおける状態の値を与えられるときの、その特定の時刻ステップにおける対象となる状態の確率に関連しており、もう一方は、その特定の状態値を与えられるときの信号又は観測値の確率に関連している。
【0012】
隠れマルコフモデリング手法は、ECG波形特徴の固有の統計的な特徴を利用するので、ECG解析の場合に特に適している。さらに、そのモデルは、波形特徴の連続的な性質(すなわち、P波の後にQRS群が続き、その後にT波が続いている)も巧みに利用することができ、それは雑音のあるECG信号のセグメンテーションを改善するのを助ける。
【0013】
残念なことに、隠れマルコフモデルは、極めて信頼度の低いセグメンテーションを引き起こす傾向があるという重大な短所を抱えている。これらのセグメンテーションは、セグメント化された特徴が典型的には少数の時刻サンプルしか続かないという特徴を有する。この問題は、非常に短い持続時間を有する特徴に向かってモデルが偏ることから起こり、それゆえ、標準的なHMMは、高い度合いの頑強性を必要とするECGセグメンテーションのような応用形態には適していない。
【0014】
図2及び図4は、専門家であるECG解析者(analyst)によって求められるような波形境界(垂直な実線)及び標準的な隠れマルコフモデルを用いて推測される波形境界(垂直な破線)とともに、2つのサンプルECG波形を示す。いずれの事例でも、元の信号内にそれぞれ1つしか存在しないにもかかわらず、そのモデルは2つのQRS群、及び2つのT波を見いだしている。詳細には、モデルセグメンテーションは、非常に短いQRS群及びT波を含み、それは、P波の開始直後、及び本当のQRS群位置の直前に生じる。そのような「二重拍動」セグメンテーションは、モデルの性能及び信頼度に深刻な影響を及ぼす。
【0015】
この背景に留意して、本発明の1つの目的は、ECG信号のような信号のセグメンテーションの信頼度を高めることに関する。
【0016】
本発明のさらに別の目的は、そのモデルから信頼度指標を導出することである。この信頼度指標を用いて、セグメンテーションの品質及び頑強性を評価することができ、さらに、自動化システムによる信頼度のあるセグメンテーションが不可能であるほど雑音が多すぎるか、又は著しく異常である、いかなるECG波形をも検出することができる。
【0017】
本発明の一態様によれば、隠れマルコフモデルを用いた生物医学信号のセグメンテーションのための、コンピュータによって実施される方法であって、当該モデルは複数の状態を含み、
当該方法は、
状態のうちの少なくとも1つについて、最小持続時間制約dminを指定するステップと、
指定された最小持続時間を有するモデル内の状態毎に、当該状態を1組のサブ状態で置き換えるステップであって、当該サブ状態の総数は最小持続時間dminの値に等しい、置き換えるステップと、
1組のサブ状態を連結するステップであって、それによって、左−右マルコフ連鎖を形成し、最初のdmin−1個のサブ状態はそれぞれ、0の自己遷移確率と、その右隣の状態に遷移する1の遷移確率と、モデル内の任意の他の状態に遷移する0の遷移確率とを有する、連結するステップと、
モデルを、生物医学信号を表すデータに適用するステップであって、それによって、信号の状態へのセグメンテーションについて情報を得る、モデルを適用するステップと
を含む方法が提供される。
【0018】
本発明の別の態様によれば、ビタビアルゴリズムに対する変更形態を用いて、観測のシーケンスを含む信号を、有限状態離散時間マルコフ過程の状態のシーケンスにセグメンテーションするための、コンピュータによって実施される方法であって、変更形態は、
ビタビアルゴリズムに組み込まれる有限状態離散時間マルコフ過程の少なくとも1つの状態について持続時間変数及び持続時間制約を定義することであって、持続時間制約は当該少なくとも1つの状態について最小持続時間を指定する、定義することと、
信号に変更されたビタビアルゴリズムを適用することであって、それによって、観測のシーケンスを説明する最も確からしい持続時間制約付き状態シーケンスを計算する、適用することと
を含む方法において、
各時刻ステップにおいて、有限状態離散時間マルコフ過程内の状態毎に、その時刻ステップまでの観測のシーケンスを説明するとともに状態で終わっている、最も確からしい状態シーケンスを計算する際に、
持続時間制約を有する状態毎に、その状態のための持続時間変数を用いて、その状態だけから成るとともに直前の時刻ステップにおいてその状態で終わっている、一連の先行状態のシーケンスの長さを絶えず注視し、
その状態のための持続時間変数が、その状態のための指定された持続時間制約以上である場合には、所与の時刻ステップにおける状態シーケンス計算において、その状態からマルコフ過程内の任意の他の所与の状態への遷移が考慮され、
その状態のための持続時間変数が、その状態のための指定された持続時間制約より小さい場合には、所与の時刻ステップにおける状態シーケンス計算において、その状態からマルコフ過程内の任意の他の状態への遷移が考慮されず、
かつ、方法において、
所与の時刻ステップまで1組の最も確からしい状態シーケンスを計算した後に、その特定の状態だけから成るとともに考慮されたばかりの時刻ステップにおいてその状態で終わっている、一連の先行状態シーケンスの長さを絶えず注視するために、持続時間制約を有する状態毎に持続時間変数を更新する、方法が提供される。
【0019】
本発明のさらなる態様は、確率的なセグメンテーションアルゴリズムによってセグメンテーションされた信号を解析するための、コンピュータによって実施される方法であって、
複数のセグメンテーションされた信号特徴毎に信頼度指標を計算することと、
個々のセグメンテーションされた信号特徴長に対して信頼度指標をプロットすることと、
密度モデリング技法を適用することであって、それによって、高い信頼度特徴に関連付けられるデータ空間の適切な領域を求める、適用することと、
特定の信号特徴のための信頼度指標がこの領域の外側に該当するか否かを判定することと
を含む方法を提供する。
【0020】
ここで、一例にすぎないが、添付の図面を参照しながら、本発明の実施形態を記述する。
【0021】
本発明の以下に記述される実施形態は、ECG信号の特徴のうちの1つ又は複数を認識するようにトレーニングされる隠れマルコフモデル(HMM)の基礎を基にする。そのモデルは複数の状態から構成され、各状態はECG信号の特定の領域を表す。ECGセグメンテーションのための隠れマルコフモデルの1つの形態のアーキテクチャ又は「トポロジ」の図式的な表示が図7aに示される。このモデルは6つの固有の状態から構成され、それらの状態は、P波、P波の終了とQRS群の開始との間の基線の部分(「基線1」と呼ばれる)、QRS群、T波、U波及びT波(又はもし存在する場合にはU波)の終了と後続の心拍のP波の開始との間の基線の部分を表す。ECG波形のこれらの部分は図1に示される。ECGセグメンテーションのための隠れマルコフモデルアーキテクチャの多数の代替の形態が図7b〜図7eに示される。それぞれの場合に、そのモデルは、「隠れ」状態シーケンス(白抜きのノードによって示される)から構成され、そのシーケンスは観測される信号(陰影付きのノードによって示される)に確率論的に関連付けられる。ECGセグメンテーションの場合、隠れ状態sは、時刻tにおいて活動している特定の波形特徴に対応し、観測される信号サンプルOは、ECG波形の関連する信号サンプルに対応する。
【0022】
そのモデルは3つのパラメータ、すなわちi)初期状態分布π、ii)遷移行列A及びiii)HMMの状態k毎の観測確率モデルbによって支配される。HMMをトレーニングすることは、特定の誤差関数を最小にするために、これらのパラメータを調整することを含む。
【0023】
そのモデルパラメータのために適した1組の値を与えられるとき、HMMはECG信号のための生成確率モデルの形をとる。より厳密には、そのモデルはECG波形を生成するための確率的な手順を定義する。この手順は初期状態の選択で開始し、それは、初期状態分布πからのサンプリングによって達成される。その際、選択された特定の初期状態は、ECG信号の第1の波形特徴(すなわち、時刻サンプルt=1において活動している特徴)に対応する。その後、この波形特徴のための関連するECG信号サンプルが、その特定の状態のための対応する観測確率モデルbからのサンプリングによって生成される。
【0024】
その手順における次のステージは、次のモデル状態、すなわち時刻サンプルt=2において活動している波形特徴を確率的に選択することである。これは、遷移行列からのサンプリングによって達成される。この行列は、そのモデル内の状態毎に、次の状態が生じる確率を定義する簡単な表の形をとる。一旦、このようにして次の状態が選択されたなら、そのモデルは、適切な観測確率モデルからのサンプリングによって、関連するECG信号サンプルを生成する。この手順は、残りの時刻ステップについて連続的に繰り返される。
【0025】
特定のECG信号を与えられるとき、その信号を、モデルアーキテクチャによって定義される、対象となる種々の特徴又は領域にセグメント化すること(すなわち、各モデル状態の信号境界の位置を特定すること)は、上記の「生成手順」を逆にすることを含む。したがって、その目的は、所与のECG信号を引き起こす隠れ状態の最も確からしいシーケンスを推測することである。その際、この状態シーケンスは、ECG信号内の対象となる領域の境界を、それゆえ対応するセグメンテーションを定義する。
【0026】
このセグメンテーション手順は、ビタビアルゴリズムを用いることによって達成することができる。より厳密には、L. R. Rabiner著「A Tutorial on Hidden Markov Models and Selected Application in Speech Recognition」(Proceedings of the IEEE, 77: 257-286, 1989)に記述されるように、ビタビアルゴリズムは、トレーニングされた隠れマルコフモデルとともに所与のECG信号を入力として取り込み、最も確からしい隠れ状態シーケンスを返す。
【0027】
そのモデルがECG波形のセグメンテーションという作業を実行するのに成功するために、ECGの、頑強でかつ情報を与える符号化を用いる必要がある。特に有効な表現は、N. P. Hughes、L. Tarassenko及びS. J. Roberts著「Markov Models for Automated ECG Interval Analysis」(Advances in Neural Information Processing Systems 16, MIT Press, 2003)に記述されるような、ECG信号のアンデシメーテッド(undecimaetd)ウェーブレット変換(UWT)からの1組の係数である。
【0028】
一旦、ECGのための表現が選択されたなら、HMMをトレーニングすることは、ECGデータの確率(又は尤度)、すなわちp(O|λ)を最大にするために、そのモデルのパラメータを調整することを含む。ただし、OはECG信号のデータセットを表し、λは1組のHMMパラメータを表す。トレーニングは、教師あり(supervised)で、又は教師なし(unsupervised)で達成することができる。
【0029】
教師あり学習の場合、トレーニングデータは、多数のECG波形および対応する状態ラベルから成る。したがって、所与のECG波形からの信号サンプルO毎に、そのサンプルがどの特定の状態(すなわち、波形特徴)に属するかを指示するラベルlも入手できるようにする。その後、以下の最尤推定を用いて、初期状態分布及び遷移行列を計算することができる。
【0030】
【数1】

【0031】
ただし、mは、第1の信号サンプル(全てのECG波形にわたる)が状態iに属する全回数であり、nijは全てのラベルシーケンスにおける状態iから状態jへの全遷移数であり、Nはデータセット内の個々のECG信号の数である。
【0032】
特定の状態に属する全ての信号サンプルを抽出し、その後、観測モデルをそれらのサンプルに直に当てはめることによって、その状態k毎の観測確率モデルbを学習することができる。ガウス観測モデルの場合には、これは単に、与えられたサンプルの平均及び分散を計算することを含む。ガウス混合モデルの場合には、与えられたサンプルから観測モデルのパラメータを学習するために、EMアルゴリズム(期待値最大化アルゴリズム)を用いる必要がある。
【0033】
教師なし学習の場合、トレーニングデータは、ECG波形だけから成る。その後、バウム−ウェルチアルゴリズム(それは、より一般的なEMアルゴリズムの特殊な場合である)を用いて、隠れマルコフモデルをトレーニングすることができる。
【0034】
そのモデルが、信頼度のないセグメンテーションを引き起こすのを防ぐために、本発明の実施形態は、そのモデルに1組の持続時間制約を課す。これらの制約は、そのモデルが生理的に(関連する波形特徴の持続時間に関して)起こりそうもないセグメンテーションを推測する能力を制限する。
【0035】
詳細には、持続時間制約は、モデル内の状態毎にただ1つの数dminの形をとり、それは、その特定の状態のための(「サンプル」内の)最小許容持続時間を表す。これらの持続時間制約は、多数の異なる方法で推定することができる。ラベル付きのECGデータが入手可能である場合には、波形特徴毎の持続時間制約は、データセット内に存在するその特徴の最小持続時間の適切な割合として簡単に推定することができる。たとえば、最小T波持続時間は、データセット内に存在する最小T波持続時間の80%に設定することができる。
【0036】
ラベル付きのデータが入手できない場合に、且つ使用できるものが生のECGデータだけである場合には、持続時間制約は、そのデータから学習されなければならない。これは、トレーニングされたモデルを用いてECGデータをデータセットにセグメンテーションし、その後、結果としてもたらされたセグメンテーションから(ラベル付きのデータの場合に先に説明されたのと同じようにして)持続時間制約を推定することによって達成することができる。この場合には、最初に信頼度のない「二重拍動」セグメンテーションをすべて除去する必要があり、それを除去しなければ、推定された持続時間制約の品質に悪影響を及ぼすであろう。
【0037】
一旦、持続時間制約が計算されたなら、それらを用いて、隠れマルコフモデルのセグメンテーション性能を改善することができる。これは、持続時間制約をビタビアルゴリズムに組み込むことによって、又は持続時間制約をHMMアーキテクチャに組み込むことによって達成することができる。ここで、これらの両方の方法を詳細に検討する。
【0038】
ビタビアルゴリズムは、L. R. Rabiner著「A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition」(Proceedings of the IEEE, 77: 257-286, 1989)に詳述されるように、隠れマルコフモデルを実際にうまく利用できるようにする重要な推測手順である。ECG解析との関連では、ビタビアルゴリズムによれば、所与のECG波形におけるP波、QRS群及びT波のための最も確からしい境界位置(すなわち、オンセット及びオフセット)を特定できるようになる。
【0039】
持続時間制約をビタビアルゴリズムに組み込むためには、関連する波形特徴のそれぞれの持続時間が所与の持続時間制約に従うように、そのアルゴリズムによって返される状態シーケンスを制限する必要がある。
【0040】
より形式的には、HMMλ(K個の状態を有する)及びECG信号O=O...Oを与えられるとき、k=1...K毎に持続時間制約dmin(k)を満たす最も確からしい状態シーケンスS=s...sを推測したい。各持続時間制約は、状態kの任意の「実行(run)」の持続時間が少なくともdmin(k)サンプルでなければならないという要件を具現する。それゆえ、その目的は、(個々の状態実行の持続時間に関して)無効であると予めわかっている状態シーケンスを推測しないようにすることを確実にすることである。
【0041】
標準的なビタビアルゴリズムは、最初のt回の観測を説明するとともに状態iにおいて終了する、最も確からしい状態シーケンスの尤度を計算するための効率的な手順に基づく。
【0042】
【数2】

【0043】
そのアルゴリズムの中心には、δ(i)を計算するための以下の漸化関係がある。
【0044】
【数3】

【0045】
ただし、上記の関係は、状態シーケンスが一次マルコフであるという事実を利用している(すなわち、δ(i)を計算する際に、直前の状態しか考慮に入れる必要がない)。その際、任意の時点tにおいて、関連付けられた遷移確率が0でないなら、状態トレリスの中の最も確からしいパスは、その現在の状態から異なる状態に常に遷移することができる。したがって、持続時間制約をビタビアルゴリズムに組み込むために、このパスは、対象となる現在の状態が最小許容持続時間以上の時間にわたって占有されている場合にのみ、1つの状態から別の状態への任意の遷移が許されるように制限されなければならない。
【0046】
本発明のこの第1の実施形態による持続時間制約付きビタビアルゴリズムは以下のように進む。漸化式の各時刻ステップtにおいて、時刻t−1において終了する、状態iの現在の「実行」の持続時間d(i)を絶えず注視する。より厳密には、状態毎の持続時間変数を用いて、その特定の状態だけから成るとともに直前の時刻ステップにおいてその状態で終わっている、一連の先行状態シーケンスの長さを絶えず注視する。この持続時間が、最小許容持続時間dmin(i)以上である場合には、δ(j)を計算するときに、この状態から任意の他の状態jへの遷移を考慮する(標準的なビタビ漸化式の場合と同様)。しかしながら、現在の状態の占有が最小許容値より低い場合には、この状態から異なる状態への遷移は許されず、その状態は(その最小持続時間制約が満たされるまで)自己遷移しかできない。δ(i)がtの全ての値について計算されると、標準的な後戻り(backtracking)法を用いて、最終的な状態シーケンスを生成することができる。
【0047】
持続時間制約付きビタビアルゴリズムの手順全体が、図6の流れ図において概略的に示される。そのアルゴリズムの第1のステージは、HMMパラメータ、観測された信号、及び状態毎の最小持続時間制約を入力することである。これは、ボックス10に示される。
【0048】
ボックス20において示される、そのアルゴリズムの次のステージは、対象となる変数の初期化に関する。より厳密には、状態i毎に、δ(i)=π(O)、ψ(i)=0及びd(i)=1が設定される。標準的なビタビ手順と同じように、変数ψ(i)は、そのモデルが時刻tにおいて状態iにあった場合に、時刻t−1における最も確からしい状態を格納するために用いられる。変数dは持続時間制約付きビタビ手順に固有であり、時刻サンプルt−1において終了する、各状態の現在の「実行」の持続時間を絶えず注視するために用いられる。より厳密には、状態i毎に、持続時間変数d(i)を用いて、状態iだけから成るとともに時刻ステップt−1においてその状態iにおいて終了する、一連の先行状態シーケンスの長さが絶えず注視される。
【0049】
持続時間制約付きビタビアルゴリズムの主要部がボックス30及び40において示される。これら2つの手順は、2〜Tの全ての時刻サンプルについて連続的に繰り返される(ただし、Tは信号の長さである)。ボックス30は、所与の時刻ステップtにおける状態i毎のδ(i)及びψ(i)の値を計算する。標準的なビタビアルゴリズムでは、δ(i)及びψ(i)を計算する際にそれぞれ用いられる「max」演算子及び「argmax」演算子が、全ての可能な直前の状態(1〜K)について評価される。しかしながら、持続時間制約付きビタビアルゴリズムでは、これらの演算子は、それらの最小持続時間制約dminを満たす状態についてのみ評価される。さらに、それらの演算子は、最も確からしい直前の状態を現在考慮中の状態とすることが常に許されるように、現在の状態iについても評価される。
【0050】
これを仮定すると、それゆえ、ボックス30の最初のステップは、その最小持続時間制約を満たす1組の状態j(デフォルトで状態iを含む)を計算することである。すなわち、それらの状態kはd(k)≧dmin(k)を満たす。その後、δ(i)及びψ(i)の計算がこの組の状態jについてのみ評価される。
【0051】
持続時間制約付きビタビアルゴリズムにおける次のステージは、状態毎の持続時間変数を更新することである。より厳密には、最も確からしい直前の状態ψ(i)も状態iである場合には、この状態のための持続時間変数が1だけインクリメントされる。すなわちd(i)=d(i)+1に設定される。これは、「自己遷移」の場合に対応する。その場合、そのモデルは連続した時刻ステップにわたって同じ状態のままであり、それゆえ、時刻サンプルtにおいて終了する状態iの現在の「実行」の持続時間が1だけ増加する。しかしながら、最も確からしい直前の状態ψ(i)が状態iでない場合には、この状態の持続時間変数は1にリセットされる。すなわちd(i)=1に設定される。これは、異なる状態から状態iへの遷移の場合に対応し、それゆえ、状態iの自己遷移の現在の「実行」を終了する。
【0052】
ボックス30及び40が時刻サンプルTまで計算されたなら、標準的なビタビ後戻り手順によって、以前に計算されたδ変数及びψ変数を用いて最も確からしい状態シーケンスを抽出することができる。この後戻り手順の最初のステップは、ボックス50に示されるように、δの値を最大にする、すなわち以下の式を満たす状態を見つけることである。
【0053】
【数4】

【0054】
この値が計算されたなら、ψ(s)の値から、直前の時刻ステップ(すなわち、T−1)において最も確からしい状態を簡単に見つけることができる。その後、この帰納的な手順を用いて、全ての残りの時刻サンプルについて最も確からしい状態値を計算することができる。
【0055】
本発明の第2の実施形態によれば、HMMのセグメンテーション性能を改善する代替の方法は、持続時間制約をモデルアーキテクチャに直に組み込むことである。このようにして、一旦、制約が組み込まれたなら、その後、標準的なビタビアルゴリズムを用いて、所与のECG波形の最適な持続時間制約付きセグメンテーションを推測することができる。
【0056】
持続時間制約をHMMアーキテクチャに組み込むことは、そのモデルの各状態を1組のdmin個の新たな状態で置き換えることを含む。これらの状態は、単純に左右に連結され、共通の観測密度を共有する。
【0057】
より厳密には、そのモデルは、dmin(k)(ただし、dmin(k)>1)の最小持続時間を有する状態k毎に、元の状態kの直前にあるdmin(k)−1個の付加的な状態で拡張される。付加的な状態はそれぞれ、0の自己遷移確率と、その右隣の状態に遷移する1の遷移確率と、そのモデル内の任意の他の状態に遷移する0の遷移確率とを有する。こうして、まとめて考えると、これらの状態は単純な左−右マルコフ連鎖を形成し、その連鎖内の各状態は、(その連鎖を用いる任意の実行中に)長くとも1つの時刻サンプルの間だけ占有される。
【0058】
この連鎖の最も重要な特徴は、状態毎の観測確率モデルのパラメータが元の状態kのパラメータと同じであることである。このようにモデルパラメータを共有(share)することは一般的に、「共有化(tying)」と呼ばれる。こうして、元の状態kの代わりに用いるdmin個の新たな状態に関連付けられる観測は、元の状態と同じ組のパラメータによって支配される。
【0059】
持続時間制約をHMMアーキテクチャに組み込むための手順全体が図8に図式的に示される。その図の左側には、元の隠れマルコフモデルアーキテクチャからの状態(図7に示される)のうちの1つがある。その図の右側には、この状態が持続時間制約を組み込むことによって変形される新たな構造がある。詳細には、dmin(k)−1個の新たな状態(自己遷移ループを持たない)は、単純に左から右へと元の状態に連結され、これらの状態のための観測確率モデルのパラメータは元の状態の対応するパラメータに「共有化される」。この後者の手順は、その図の破線のボックスによって指示される。
【0060】
上記の過程は、元の隠れマルコフモデル内の元の状態(すなわち、P、B1、QRS、T等)毎に繰り返される。したがって、持続時間制約をモデルアーキテクチャに組み込む全体的な効果は、そのモデル内の元の単一の状態がそれぞれ、新たな「持続時間制約」構造によって置き換えられることである。図9は、結果として構成される持続時間制約HMMアーキテクチャの一部を示す。そのモデル内の各波形特徴(たとえば、P波、QRS群等)が、ここでは複数の状態に関連付けられる。その際、この複数の状態は、ECG波形特徴毎の最小持続時間要件をカプセル化(encapsulate)する。
【0061】
持続時間制約が、先に説明されたようにしてHMMに組み込まれたなら、その後、そのモデルは、標準的なビタビアルゴリズムを用いた新たなECG波形のセグメンテーションに用いることができる。その際、そのモデルはもはや、その持続時間が許容される持続時間より短い状態シーケンスを生成することができないので、結果として生成されるセグメンテーションは、持続時間制約に従うことを保証される。
【0062】
図3及び図5は、図2及び図4において先に検討されたのと同じECG波形を示す。図3及び図5では、専門家であるECG解析者によって求められるような波形境界は、垂直な実線によって示されており、本発明の一実施形態による、持続時間制約を利用する隠れマルコフモデルを用いて自動的に求められる状態境界は、垂直な破線として示される。そこでは、モデルと専門家である解析者とが概ね一致するだけでなく、図2及び図4に示されるような、従来モデルの信頼度のない「二重拍動」セグメンテーションが除去されており、結果的な解析は、はるかに高い頑強性を有する。それらの境界は、特定の状態の持続時間、たとえばQT間隔を求める等の、さらに別の解析に、再び自動的にかけることができ、一連のパルスを含む信号にわたる、そのような間隔の変化の統計的な解析を実行することができる。それらの結果は、生の境界情報又はさらに別の解析の結果の形で、たとえば、コンピュータシステムのディスプレイ上に表示することによって出力される。
【0063】
本発明の第3の実施形態によれば、ECG波形の正確なセグメンテーションに加えて、トレーニングされたモデルからの「信頼度指標」(又はスコア)を生成することもできる。この信頼度指標は、結果として生成されたセグメンテーションにおいて、そのモデルがどの程度信用できるかを定量化し、それを用いて、セグメンテーションされたECG波形の品質を評価することができる。それゆえ、この指標を用いるとき、雑音が多すぎるか、又は「正常でない」(波形形態に関して)ECG波形を自動的に検出し、信頼度のあるセグメンテーション、すなわちQT間隔の信頼度のある推定値を提供することができる。信頼度指標は、確率的に動作する任意のセグメンテーションアルゴリズムとともに用いることができる。
【0064】
信頼度指標は、モデルを与えられるときの信号の対数尤度(log likelihood)、すなわちlog p(O...O|λ)と定義することができるか、又は(モデルを与えられるときの)信号及び最も確からしい状態シーケンスs...sの結合(joint)対数尤度、すなわちlog p(O...O,s...s|λ)と定義することができる。対数尤度は、生の尤度とは対照的に、結果として生成される値をコンピュータのダイナミックレンジ内に保持するために用いられる。
【0065】
信頼度指標のいずれの定義も、隠れマルコフモデルの場合に、さらには隠れセミマルコフモデル(L. Thoraval、G. Carrault及びF. Mora著「Continuously variable durations HMMs for ECG segmentation」(Proceedings IEEE-EMBS, 1992, pp. 529-530)に記述される)、階乗(factorial)隠れマルコフモデル(Z. Ghahramani及びM. I. Jordan著「Factorial hidden Markov models」(Machine Learning, Vol. 29, Issue 2, Nov. 1997)に記述される)のような、これらのモデルのさらに高度な形式の場合に、効率的に計算することができる。最初の形式は、そのモデルを与えられるときの信号の対数尤度と定義され、以下の式によって与えられる、標準的なHMMフォワード変数(forwards variable)から計算することができる。
【0066】
【数5】

【0067】
結果として生成される値をコンピュータのダイナミックレンジ内に保持するために、各時刻ステップにおいてこの変数をスケーリングするのが一般的である。これは、フォワード変数を、以下の式によって与えられるスケーリング係数で割ることによって達成される。
【0068】
【数6】

【0069】
その際、そのモデルを与えられるときの信号の対数尤度は、以下のように簡単に計算することができる。
【0070】
【数7】

【0071】
第2の形式の信頼度指標は、(モデルを与えられるときの)信号と最も確からしい状態シーケンスとの結合対数尤度と定義され、ビタビアルゴリズムとの関連で先に説明されたように、デルタ変数δ(i)から計算することができる。同様に、以下の式によって与えられるスケーリング係数を用いて各時刻ステップにおいてこの変数をスケーリングするのが一般的である。
【0072】
【数8】

【0073】
その際、(モデルを与えられるときの)信号と最も確からしい状態シーケンスとの対数尤度は、以下のように簡単に計算することができる。
【0074】
【数9】

【0075】
上記のスケーリング手順を用いることへの代替形態は、L. R. Rabiner著「A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition」(Proceedings of the IEEE, 77: 257-286, 1989)に記述されるように、log α(i)又はlog δ(i)の値を直に計算することである。この場合、信頼度指標は、帰納中に直に計算され、さらに別の処理は不要である。
【0076】
多くの状況において、所与のECG信号は、対象となる多数のECG特徴(多数のECG拍動又は多数のQT間隔等)を含むことになり、それゆえ、信号内の対象となる個々の特徴毎に個別の信頼度指標を計算することが望ましい(信号全体のための集団の指標ではない)。時刻サンプルt1で開始し、時刻サンプルt2において終了する特定の特徴について、その特徴のための信頼度指標は、時刻t2までの信号部分に対する信頼度指標から、時刻t1まで(ただし、t1を含まない)の信号部分に対する信頼度指標を減算すること、すなわち以下の式によって、簡単に計算することができる。
【0077】
【数10】

【0078】
又は、信頼度指標の第2の定義を用いて、以下の式によって計算することができる。
【0079】
【数11】

【0080】
この信頼度指標を実際に用いて、特定の特徴のセグメンテーションにおいて、そのモデルが有する信頼度を評価するときには、解析中に、問題となる特徴の長さが考慮に入れられなければならない。信頼度指標の形式がいずれも特定のECG特徴の長さとともに変化するので、これが必要とされる。
【0081】
ECG特徴の長さを信頼度指標解析に組み込む最も実効的な方法は、単に、多数の正常なECG特徴(すなわち、異常と見なされないか、又は雑音によって汚染されていない特徴)の場合に、特徴長(x軸上)に対する信頼度指標(y軸上)のプロットを作成することである。図10は、そのようなプロットを示しており、対象となる特定の特徴は、個々のECG波形(すなわち、1つの心拍全体に対応するECG信号の部分)であると見なされる。
【0082】
ECG特徴長に対する信頼度指標のプロットを与えるとき、標準的な密度モデリング技法を用いることによって、高い信頼度指標に関連付けられるこの2−D空間の領域を特定することができる。詳細には、正常な(すなわち高い信頼度の)ECG特徴に対応する1組の2−Dデータ点は、ガウス分布又は確率的な線形回帰モデルのような、任意の適当な密度関数によってモデル化することができる。
【0083】
所与の密度関数は自然に、この2−Dデータ空間内の高い信頼度の領域のための境界を定義する。ガウス分布の場合、この境界は、その平均から2つの標準偏差内に存在するデータ空間の領域であると見なすことができる。確率的な線形回帰モデルの場合、この境界は、適当な統計的有意水準にある平均回帰線のための下限信頼度よりも上に存在するデータ空間の領域であると見なすことができる(R. J. Freund及びW. J. Wilson著「Regression Analysis」(Academic Press, 1998)に記述される)。この下限(99%の有意水準にある)は、図10において破線によって示される。
【0084】
セグメンテーションされた特徴のそれぞれの長さとともに、多数の新たなECG特徴のセグメンテーションのための信頼度指標を与えられるとき、任意の信頼度指標が高い信頼度領域に対する境界の外側に該当するか否かを判定することは簡単なことである。その際、それらの信頼度指標がこの領域の外側に該当する特徴のためのセグメンテーションの結果は、最終的な解析から除去することができるか、又は別法では、これらの特徴は、訓練を受けたECG解析者によってさらに手動で解析するためにハイライトすることができる。
【0085】
図11はECG信号の一部を示し、それは、人為的な筋雑音によって汚染された2つのECG波形を含む。またその図は、持続時間制約を組み入れた隠れマルコフモデルによる各波形のセグメンテーションも示す。その2つの波形は、下側のx軸上で「A」及び「B」を付される。
【0086】
残念なことに、人為的な筋雑音は、標準的なECG信号の帯域幅に類似の帯域幅を有し、それゆえ、簡単なフィルタリング技法によって除去することはできない。その雑音は、ECG信号の根底を成す波形特徴を隠し、波形の正確なセグメンテーションを非常に難しくする。このため、そのような信頼度のないセグメンテーションを検出できることが重要である。
【0087】
図12は、雑音のある2つのECG波形の場合の、波形の長さに対する信頼度指標のプロットを示す。99%下限信頼度(破線)も示される。雑音のある2つの波形のための信頼度スコアはいずれも、この線の下側にあり、それゆえ、これらの波形が、信号の自動解析中に自動的に検出されるようにする。
【0088】
本発明の上記の実施形態は、ECG信号解析との関連で記述されてきた。しかしながら、本発明はこの応用形態には限定されず、他の生物医学信号の解析のためにも用いることができる。適する生物医学信号の例には、脳波図(EEG)信号(その場合には、解析は無呼吸(apneoa)を検出するために用いることができる);筋肉の電気的活動を測定する筋電図(EMG);心臓電気生理現象の場合のような侵襲的に得られる心信号(その場合には、カテーテルが心臓に挿入され、各カテーテルは多数の電極を含む);及びペースメーカの移植前に用いられる、人体に電極を挿入することによって心臓から導出される信号を含む心内(endocardial)電位図信号が含まれる。持続時間制約付きビタビアルゴリズム及び信頼度指標技法は、隠れマルコフモデルとともに用いることには限定されない。たとえば、持続時間制約付きビタビアルゴリズムは、任意の有限状態離散時間マルコフ過程とともに用いることができ、信頼度指標は確率的に動作する任意のセグメンテーションアルゴリズムとともに用いることができる。持続時間制約付きビタビアルゴリズムは、生物医学以外の信号とともに、たとえば誤り検出/訂正符号の分野において用いることもできる。
【0089】
本発明の実施形態は、コンピュータシステム上で実行されるコンピュータプログラムによって実行される。そのコンピュータシステムには、任意のタイプのコンピュータシステムを用いることができるが、一般的には、任意の適当な言語において書かれるコンピュータプログラムを実行する従来のパーソナルコンピュータである。そのコンピュータプログラムはコンピュータ読取り可能媒体に格納することができ、そのコンピュータ読取り可能媒体は、任意のタイプであってよい。たとえば、コンピュータシステムのドライブに挿入することができ、磁気的、光学的又は光磁気的に情報を格納することができる、ディスク形媒体のような記録媒体;ハードドライブのような、コンピュータシステムの固定記録媒体;又は固体コンピュータメモリである。解析されるべき信号は、たとえばECG装置から直にコンピュータシステムに入力することができるか、又はコンピュータシステムが、予め得られている信号の記憶装置から、信号を表す情報を読み出すことができる。
【図面の簡単な説明】
【0090】
【図1】十分に標識が付けられたECG波形を示す図である。
【図2】専門家であるECG解析者によって求められるような波形境界、及び標準的な隠れマルコフモデル手法を用いて推測される波形境界とともにECG波形を示す図である。
【図3】専門家であるECG解析者によって求められるような波形境界、及び新規の持続時間制約手法を用いて推測される波形境界とともに、図2において用いられるのと同じECG波形を示す図である。
【図4】専門家であるECG解析者によって求められるような波形境界、及び標準的な隠れマルコフモデル手法を用いて推測される波形境界とともに、第2のECG波形を示す図である。
【図5】専門家であるECG解析者によって求められるような波形境界、及び新規の持続時間制約手法を用いて推測される波形境界とともに、図4において用いられるのと同じECG波形を示す図である。
【図6】新規の持続時間制約付きビタビアルゴリズムの動作を詳述する流れ図である。
【図7a】ECG信号を、P波領域、基線領域、QRS群領域、T波領域、任意選択のU波領域及び第2の基線領域にセグメンテーションするための標準的な隠れマルコフモデルのアーキテクチャを示し、P波状態のための遷移確率a及び観測モデル確率分布bも合わせて示される図である。
【図7b】図7aのモデルと同じようにECG信号をセグメンテーションするが、任意選択のU波領域を用いることなくセグメンテーションするための標準的な隠れマルコフモデルのアーキテクチャを示す図である。
【図7c】ECG信号をPR間隔領域、QRS群領域、T波領域及び基線領域にセグメンテーションするための標準的な隠れマルコフモデルのアーキテクチャを示す図である。
【図7d】ECG信号をPR間隔領域、QT間隔領域及び基線領域にセグメンテーションするための標準的な隠れマルコフモデルのアーキテクチャを示す図である。
【図7e】ECG信号をQT間隔領域、及び信号の残りの部分に対応する第2の領域(X)にセグメンテーションするための標準的な隠れマルコフモデルのアーキテクチャを示す図である。
【図8】持続時間制約の隠れマルコフモデルへの組み込みを図式的に示す図である。
【図9】ECGセグメンテーションのための組み込み式の持続時間制約を用いる隠れマルコフモデルのアーキテクチャの一部を示す図である。
【図10】99%信頼度下限とともに、多数の正常なECG波形のための、ECG波形長に対する信頼度指標のプロットを示す図である。
【図11】雑音(詳細には、「筋肉のアーテファクト(muscle artefact)」)によって汚染された2つのECG波形を含み、雑音のある2つの波形(A及びBを付される)のセグメンテーションも示される、ECG信号の一部を示す図である。
【図12】図11に示される、雑音のある2つのECG波形のための信頼度指標を示す図である。

【特許請求の範囲】
【請求項1】
隠れマルコフモデルを用いた生物医学信号のセグメンテーションのための、コンピュータによって実施される方法であって、
前記モデルは複数の状態を含み、
前記方法は、
前記状態のうちの少なくとも1つについて、最小持続時間制約dminを指定するステップと、
指定された最小持続時間を有する前記モデル内の状態毎に、前記状態を1組のサブ状態で置き換えるステップであって、前記サブ状態の総数は前記最小持続時間dminの値に等しい、置き換えるステップと、
前記1組のサブ状態を連結するステップであって、それによって、左−右マルコフ連鎖を形成し、最初のdmin−1個のサブ状態はそれぞれ、0の自己遷移確率と、その右隣の状態に遷移する1の遷移確率と、前記モデル内の任意の他の状態に遷移する0の遷移確率とを有する、連結するステップと、
前記モデルを、前記生物医学信号を表すデータに適用するステップであって、それによって、前記信号の前記状態への前記セグメンテーションについて情報を得る、前記モデルを適用するステップと
を含む方法。
【請求項2】
前記連鎖内の最後のサブ状態は、置き換えられた元の状態と同じ遷移確率を有する、請求項1に記載の方法。
【請求項3】
全dmin個のサブ状態のための観測確率モデルは共有化され、
前記マルコフ連鎖内の、前記dmin個の状態のうちの任意の状態のための前記観測モデルのパラメータは、前記連鎖内の任意の他の状態のための対応するパラメータに等しい、請求項1又は2に記載の方法。
【請求項4】
前記マルコフ連鎖内の前記サブ状態毎の前記観測モデルの前記パラメータは、置き換えられた元の状態のための前記観測モデルの前記パラメータに等しく設定される、請求項3に記載の方法。
【請求項5】
前記モデルを適用する前記ステップは、ビタビアルゴリズムに、前記モデルと、セグメンテーションされるべき前記信号とを与えることを含む、請求項1〜4のいずれか一項に記載の方法。
【請求項6】
ビタビアルゴリズムに対する変更形態を用いて、観測のシーケンスを含む信号を、有限状態離散時間マルコフ過程の状態のシーケンスにセグメンテーションするための、コンピュータによって実施される方法であって、
前記変更形態は、
前記ビタビアルゴリズムに組み込まれる前記有限状態離散時間マルコフ過程の少なくとも1つの状態について、持続時間変数及び持続時間制約を定義することであって、前記持続時間制約は前記少なくとも1つの状態について最小持続時間を指定する、定義することと、
前記信号に変更された前記ビタビアルゴリズムを適用することであって、それによって、前記観測のシーケンスを説明する最も確からしい持続時間制約付き状態シーケンスを計算する、適用することと
を含む方法において、
各時刻ステップにおいて、前記有限状態離散時間マルコフ過程内の状態毎に、その時刻ステップまでの前記観測のシーケンスを説明するとともに前記状態で終わっている、前記最も確からしい状態シーケンスを計算する際に、
持続時間制約を有する状態毎に、その状態のための前記持続時間変数を用いて、その状態だけから成るとともに直前の時刻ステップにおいてその状態で終わっている、一連の先行状態のシーケンスの長さを絶えず注視し、
その状態のための前記持続時間変数が、その状態のための前記指定された持続時間制約以上である場合には、所与の時刻ステップにおける前記状態シーケンス計算において、その状態から前記マルコフ過程内の任意の他の所与の状態への遷移が考慮され、
その状態のための前記持続時間変数が、その状態のための前記指定された持続時間制約より小さい場合には、前記所与の時刻ステップにおける前記状態シーケンス計算において、その状態から前記マルコフ過程内の任意の他の状態への遷移が考慮されず、
かつ、前記方法において、
前記所与の時刻ステップまで前記1組の最も確からしい状態シーケンスを計算した後に、その特定の状態だけから成るとともに考慮されたばかりの前記時刻ステップにおいてその状態で終わっている、前記一連の先行状態シーケンスの前記長さを絶えず注視するために、持続時間制約を有する状態毎に前記持続時間変数を更新する、方法。
【請求項7】
持続時間変数及び持続時間制約が、前記有限状態離散時間マルコフ過程の状態毎に与えられる、請求項6に記載の方法。
【請求項8】
前記有限状態離散時間マルコフ過程は隠れマルコフモデルを含む、請求項6又は7に記載の方法。
【請求項9】
確率的なセグメンテーションアルゴリズムによってセグメンテーションされた信号を解析するための、コンピュータによって実施される方法であって、
複数のセグメンテーションされた信号特徴毎に信頼度指標を計算することと、
個々の信号特徴長に対して前記信頼度指標をプロットすることと、
密度モデリング技法を適用することであって、それによって、高い信頼度特徴に関連付けられるデータ空間の適切な領域を求める、適用することと、
特定の信号特徴のための前記信頼度指標がこの領域の外側に該当するか否かを判定することと
を含む方法。
【請求項10】
前記信号特徴は、より長い信号の一部を含む、請求項9に記載の方法。
【請求項11】
前記信頼度指標は、
前記モデルを与えられるときの前記信号特徴の対数尤度と、
前記モデルを与えられるときの前記信号特徴と最も確からしい状態シーケンスとの結合対数尤度と
から成るグループから選択される指標である、請求項9又は10に記載の方法。
【請求項12】
前記信号は、心電図信号と、脳波図信号と、筋電図信号と、心臓電気生理現象信号と、心内電位図信号とから成るグループから選択される生物医学信号である、請求項1〜11のいずれか一項に記載の方法。
【請求項13】
前記信号の解析の結果についての情報を表示することをさらに含む、請求項1〜12のいずれか一項に記載の方法。
【請求項14】
前記信号は心電図信号であり、前記信号のセグメンテーションからQT間隔についての情報が得られる、請求項1〜13のいずれか一項に記載の方法。
【請求項15】
コンピュータシステム上で実行されることができるコンピュータプログラムであって、そのように実行されるときに、前記コンピュータシステムに請求項1〜14のいずれか一項に記載の方法を実行させることができるコンピュータプログラム。
【請求項16】
請求項15に記載のコンピュータプログラムを格納するコンピュータ読取り可能媒体。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate

【図6】
image rotate

【図7a】
image rotate

【図7b】
image rotate

【図7c】
image rotate

【図7d】
image rotate

【図7e】
image rotate

【図8】
image rotate

【図9】
image rotate

【図10】
image rotate

【図11】
image rotate

【図12】
image rotate


【公表番号】特表2007−536050(P2007−536050A)
【公表日】平成19年12月13日(2007.12.13)
【国際特許分類】
【出願番号】特願2007−512325(P2007−512325)
【出願日】平成17年5月6日(2005.5.6)
【国際出願番号】PCT/GB2005/001747
【国際公開番号】WO2005/107587
【国際公開日】平成17年11月17日(2005.11.17)
【出願人】(500205220)アイシス イノヴェイション リミテッド (10)
【Fターム(参考)】