(57)【特許請求の範囲】
【請求項１】
２つの独立した測定チャネルでの、パラメータの２つの独立した測定の各々の値を予測するステップと、
２つの独立した測定チャネルでの、パラメータの２つの独立した測定によって生成された、パラメータの２つの測定値を受け取るステップと、
予測値と測定値との間の、２つの独立した測定チャネルのそれぞれについての差を計算するステップと、
それぞれの前記差に逆相関して変化する２つの重みを用いて２つの測定値を組み合わせて組み合わされた値を生成するステップと、
組み合わされた値をユーザに出力するステップと、
を含むことを特徴とする、プログラムされたコンピュータによってパラメータの測定を組み合わせる方法。
【請求項２】
予測、測定、計算、及び組み合わせのステップが連続的に繰り返され、２つの独立した測定の各々の予測値が、先行する予測値と、先行する予測値と先行する測定との間の差とに基づいていることを特徴とする請求項１に記載の方法。
【請求項３】
２つの独立した測定の各々の予測値が、線形予測モデルを使用することによって計算されることを特徴とする請求項２に記載の方法。
【請求項４】
２つの独立した測定の各々の予測値が、非線形予測モデルを使用することによって計算されることを特徴とする請求項２に記載の方法。
【請求項５】
モデルが、適応型であって測定におけるプロセスノイズの量に依存して適応することを特徴とする請求項３または４に記載の方法。
【請求項６】
２つの測定値を組み合わせるステップにおいて、各々の値の重みが、予測値と測定との間の差の二乗に逆比例して変化することを特徴とする請求項１〜５のいずれか一つに記載の方法。
【請求項７】
２つの測定値が、次の式に従って組み合わせられる（ここで、Ｍ1及びＭ2は２つの測定値であり、σ1及びσ2は２つの測定値とそれぞれの予測値との間の差である）ことを特徴とする請求項６に記載の方法。
【数１】
0004741168.tif 000004
【請求項８】
それぞれの測定の予測値が、システムのそれぞれのモデルに基づくことを特徴とする請求項１〜７のいずれか一つに記載の方法。
【請求項９】
システムのモデルが、プロセスノイズ及びセンサノイズの推定を含むことを特徴とする請求項８に記載の方法。
【請求項１０】
システムのそれぞれのモデルが相互に独立であることを特徴とする請求項８または９に記載の方法。
【請求項１１】
システムのそれぞれのモデルが、プロセスノイズ及びセンサノイズのために同じ推定を含むことを特徴とする請求項１０に記載の方法。
【請求項１２】
システムのそれぞれのモデルがカルマンフィルタであることを特徴とする請求項８、９、１０または１１のいずれか一項に記載の方法。
【請求項１３】
更に、双方の測定値とそれらの予測値との間の差が所定の時間間隔の間所定の閾値を超過する一連の測定を廃棄するステップを含むことを特徴とする請求項１〜１２のいずれか一つに記載の方法。
【請求項１４】
パラメータが心拍数であることを特徴とする請求項１〜１３のいずれか一つに記載の方法。
【請求項１５】
２つの独立した測定が心電図及びパルスオキシメトリ波形から行なわれることを特徴とする請求項１４に記載の方法。
【請求項１６】
２つの測定が多数のリード線を有する心電図の記録から行なわれることを特徴とする請求項１４に記載の方法。
【請求項１７】
３つ以上の測定が存在することを特徴とする請求項１〜１６のいずれか一つに記載の方法。
【請求項１８】
更に、双方の測定値とそれらの予測値との間の差の値に基づいて動きアーチファクトを識別するステップを含むことを特徴とする請求項１〜１７のいずれか一つに記載の方法。
【請求項１９】
請求項１〜１８のいずれか一つに記載の方法を実行するプログラム・コード手段を含むコンピュータ・プログラム。
【請求項２０】
請求項１〜１８のいずれか一つに記載の方法を実行するように構築及び配列された装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、パラメータの改善された測定を提供するため、異なったセンサからの測定を組み合わせる方法及び装置に関する。本発明は、特に、生理学的パラメータの測定に応用可能である。
【０００２】
【従来の技術】
或る種のパラメータは、複数の方法で測定されることができる。これは、同じ量の独立した測定を提供する場合に有用である。例えば、医療分野で、心拍数は、心電図（ＥＣＧ）及びパルスオキシメトリ波形（酸素の飽和を計算するために使用される）の双方から測定されることができる。
【０００３】
添付の図面において、図４はこれらの２つの波形を概略的に示す。図４ａは心拍数がＨＲ1で示された心電図であり、図４ｂは心拍数がＨＲ2で示されたパルスオキシメトリ波形である。心拍数は、突然の変化を受けることができるパラメータである。これら変化のいくつかは、有効な生理学的変化であり、例えば、早期に起こって心拍数の一時的な増加をもたらす期外収縮である。
【０００４】
図５は、心電図トレース及びパルスオキシメトリ波形の双方で発見される期外収縮５０の発生を示す。先行する鼓動と期外収縮５０との間の短縮された間隔は、心拍数の測定において心拍数の突然の増加として現れる。
【０００５】
添付図面の図１及び図２は、パルスオキシメトリ（図１）及び心電図（図２）によって測定された心拍数の時間プロットを示す。図１及び図２において、プロットの早期の部分における心拍数は、一般的に、８０ｂｐｍのオーダであるが、期外収縮によって生じ、その結果パルスオキシメトリによる測定及び心電図による測定の双方に現れている１０及び２０のように、時々、心拍数に突然の増加が存在することが分かる。
【０００６】
【発明が解決しようとする課題】
しかしながら、有効な生理学的変化から導き出された測定心拍数の変化に加えて、生理学的に有効でない他の変化、例えば、体表面に置かれたセンサの突然の動き（例えば、心電図電極がある胸部の動き）によって引き起こされる変化が生じる。図６は、パルスオキシメトリ波形において、見かけ上、鼓動の間隔を短くし、その結果、見かけ上、心拍数の増加を生じさせるアーチファクト６０の存在を示す。
【０００７】
これらの変化は、図１の１２及び図２の２２で示されるように、一方の測定には反映されるが他方の測定には反映されない。そのような変化が一方の測定に現れて他方の測定には現れないという事実は、心拍数の変化のどれが有効な生理学的変化であり、どれがアーチファクトであるかの決定を助けるために、２つの測定を組み合わせることができることを意味する。
【０００８】
しかしながら、一方の測定チャネルを他方の測定チャネルに対して有効にする通常の手法は、２つの測定の相互相関を含んで常に失敗する。なぜなら、心拍数の変化を拒絶するのではなく有効なものとして受け入れる閾値に、どのような値を与えるかを前もって（記録ごとに、患者ごとに）知ることは不可能だからである。
【０００９】
従って、図１及び図２から、２２で示されるような変化を除去する閾値を設定できるように見えるけれども、そのような閾値は、全ての患者及び全ての記録に適切であるわけではなく、パルスオキシメトリ波形を使用する場合の助けにはならない。その問題は、図１、図２、及び図３の領域ＡＦで示されるように、心拍数が急激に変化する心房細動の場合に増大する。同様の問題は、パラメータが２つ以上の測定チャネルを介して測定される他の分野で生じる。
【００１０】
【課題を解決するための手段】
本発明は、有効な変化をアーチファクトから識別できるような方法で、パラメータの２つの測定を組み合わせることによって、パラメータの測定を改善する方法及び装置を提供する。
【００１１】
従って、それは、２つの独立した測定チャネルでの、パラメータの２つの独立した測定の各々の値を予測するステップと、２つの独立した測定チャネルでの、パラメータの２つの独立した測定によって生成された、パラメータの２つの測定値を受け取るステップと、予測値と測定値との間の、２つの独立した測定チャネルのそれぞれについての差を計算するステップと、それぞれの前記差に逆相関して変化する２つの重みを用いて２つの測定値を組み合わせて組み合わされた値を生成するステップと、組み合わされた値をユーザに出力するステップと、を含む、プログラムされたコンピュータによってパラメータの測定を組み合わせる方法を提供する。
【００１２】
本発明を使用して、各々の測定のために予測が行なわれ、実際の測定がその予測と比較される。「イノベーション」と呼ばれるその差が計算され、このイノベーションは重みを計算するのに使用される。この重みは、その測定が他の測定と組み合わせられるときその測定へ与えられ、他の測定もそのイノベーションに従って重みをつけられる。
【００１３】
重みの計算は、一方の測定チャネルでのイノベーションが高く、他方の測定チャネルでのイノベーションが低い場合、低いイノベーション・チャネルからの測定が、より重く重みづけられるように計算される。これは、一方のチャネルからのイノベーションの高いレベルと他のチャネルからの低いイノベーションとが一致することは、より高いイノベーション・チャネルでのアーチファクトを示すものと考えられるからである。
【００１４】
従って、値が組み合わせられる際のそれら各値へ与えられる重みは、測定値とその予測値との間の差の２乗または絶対値に逆相関するように関連づけられる。１つの実施の形態において、その測定値は、次の式に従って組み合わせられることができる。
【００１５】
【数２】
0004741168.tif 000002
ここで、Ｍ1及びＭ2は２つの測定値であり、σ1及びσ2は２つの測定値とそれらの各予測値との間の差である。
【００１６】
予測、測定、計算、及び組み合わせのステップは、線形予測モデルに基づく測定の各々の予測値を使用して、好ましくは連続的に繰り返される。その予測値は、例えば、その先行する予測値及び先行するイノベーション（即ち、先行する予測値と先行する測定値との間の差）に基づくものである。予測値は、イノベーションの定数倍を、先行する予測値へ加えることによって得られる。その定数は、好ましくは、１より小さいか等しい正の値である。２つの独立した測定の各々の予測値については、代わりに、ニューラルネットワークのような非線形予測モデルを使用することによって計算することもできる。
【００１７】
１つの実施の形態において、予測値はシステムの数学的モデルに基づくことができる。この数学的モデルは、プロセスノイズ及びセンサ（測定）ノイズの推定を含んでよい。測定チャネルの各々のために１つのモデルが使用されるように、２つの独立したモデルが使用されてもよい。それらモデルはプロセスノイズ及びセンサノイズのための推定（２つのチャネルについて同じである場合がある）を含んでもよい。
【００１８】
１つの実施の形態において、上記モデルはカルマンフィルタである。その方法は、特に、心拍数の測定に適用できる。その場合、２つの測定チャネルは、心電図及びパルスオキシメトリ波形からのチャネルとすることができる。もっとも、その方法は、２つ以上のソースから得ることのできるパラメータの任意の他の測定に適用可能である。従って、その方法は、３つ以上の測定チャネルに適用可能であり、また、測定が独立している場合と、測定が全く独立していない場合（例えば、例えば心電図の多数のリード線に対して）との両方に適用可能である。
【００１９】
本発明は、更に、動きアーチファクトを検出することができる。この場合、動きの期間に、イノベーションの高い値が双方のチャネルで得られる。これは、その動きによって崩壊されたデータ部分を廃棄するためのトリガとして使用されることができる。
【００２０】
本発明は、コンピュータ・ソフトウェアを使用して実現できること、それにより、本発明は、方法またはその一部分を制御して実行するコンピュータ・プログラム、及びそのプログラムを搬送するコンピュータ読み取り可能な記憶媒体へ拡張されることが分かるであろう。本発明は、更に、方法を実行するための対応する装置へ拡張される。
【００２１】
【発明の実施の形態】
本発明は、添付の図面を参照しつつ、非限定的例を用いて詳細に説明される。以下に、医療分野において、心電図及びパルスオキシメトリを使用した心拍数の測定に対して適用される本発明の実施の形態を説明する。
【００２２】
図４ａで示されるように、心電図によって測定された心拍数は、２つの連続したＲ波ピークの間の間隔から導き出される。パルスオキシメトリ波形から導き出された心拍数の測定は、図４ｂで示されるように、２つの連続したピーク（または谷）の間の間隔から得られる。図１及び図２は、これら２つの測定から導き出された心拍数プロットを示す。
【００２３】
本発明のこの実施の形態において、心拍数を生成する過程のモデルが構築される。同じモデルが、各測定ソースについて（即ち、１つのモデルは心電図測定チャネルについて、及び１つのモデルはパルスオキシメトリ測定チャネルについて）独立に実行される。
【００２４】
この実施の形態において、モデルはカルマンフィルタである。一般的に、カルマンフィルタはプロセスモデル及び観測モデルを使用する。プロセスモデルは、時間ｔにおけるシステムの状態に関して、時間ｔ＋１におけるシステムの状態をモデル化する。測定または観測モデルは、時間ｔにおける測定が、時間ｔにおけるシステムの状態へどのように関連づけられるかを示す。従って、一般的な用語は次の通りである。
【００２５】
ｘt+1＝Ａｘt＋ｗ（プロセスモデル）
ｙt＝Ｃｘt＋ｖ（観測モデル）
ここで、
ｗ〜Ｎ（Ｏ，Ｑ）  −  ゼロ平均及び分散Ｑを有するガウス型プロセスノイズ
ｖ〜Ｎ（Ｏ，Ｒ）  −  ゼロ平均及び分散Ｒを有するガウス型測定ノイズ
及び、
・  ｋ  −  状態変数ｘのベクトル
・  ｎ  −  観測または測定ｙのベクトル
・  状態ｘは、単純一次マルコフ力学に従って発展する。
・  Ａは、ｋｘｋ状態遷移行列である。
・  各測定ベクトルｙは、線形観測プロセスによって現在の状態へ関連づけられる。Ｃは、ｎｘｋ観測または測定行列である。
【００２６】
この実施の形態において、一般的なカルマンフィルタはスカラー量を使用するように単純化され、同じプロセス及び測定ノイズモデル（ｗ，ｖ）が双方の測定チャネルで使用される。従って、単純化されたカルマンフィルタは次のようになる。
ｘt+1＝Ａｘt＋ｗ（プロセスモデル）
ｙt＝Ｃｘt＋ｖ（観測モデル）
【００２７】
モデルは、更に、Ｃ＝１と設定することによって単純化される。これは、心拍数が、プロセス及び測定を記述する双方の状態であることを意味する。更に、Ａ＝１と仮定される。これは、次の心拍数が、プロセスノイズモデルによって許容される可変性を有する先行する心拍数と同じであることを意味する。
【００２８】
従って、各チャネルで、このモデルを使用すると、２つの測定を組み合わせる過程は、次のステップを含む。
１．時間ｔまでの先行履歴の知識から、
（ａ）次の状態ｘｐｒｅｄを予測する。
（ｂ）ｘｐｒｅｄから、次の測定ｙｐｒｅｄを予測する。
２．測定ｙｔ＋１を行なう。
３．イノベーションを計算する。
ｙｔ＋１＝ｙｐｒｅｄ＋εｔ＋１
ここで、εｔ＋１は、実際の値と予測値との間の差、即ちイノベーションである。
４．分散σｔ＋１２を計算する。
σｔ＋１２＝εｔ＋１２
分散σ２は、予測に関連づけられた「信頼度」に逆相関している。
５．心拍数測定を、各々の測定に関連づけられた分散と逆比例させて混合する。
【数３】
0004741168.tif 000003
【００２９】
ＭＡＴＬＡＢにおけるこのモデルの実施の形態の例が、付録１に示される。その例は一般的なもので、ベクトル量に対して働くであろう。もっとも、この実施の形態において、量はスカラーである。付録１から、各々の新しい測定サイクルの予測心拍数（ｘｎｅｗ）は、前に予測された値（ｘｐｒｅｄ）にｅ（イノベーション）のＫ（カルマンゲイン）倍を加えたものに等しいことが分かる。カルマンゲインＫは、予測分散Ｖｐｒｅｄ及び測定分散Ｒから導き出される。予測分散は、前の予測分散及びプロセスノイズＱから導き出される。プロセスを開始するため、プロセスは、心拍数の初期値８０及び状態分散の初期値１００を使用して初期化される。この実施の形態におけるＱのプロ製ノイズは５に設定され、測定ノイズ分散Ｒは１０に設定される。
【００３０】
実施の形態から、カルマンフィルタを用いた通常のように、分散及びカルマンゲインは測定値に依存しないことが明らかであろう。測定値は、イノベーションｅを介する心拍数の新しい予測でのみ使用される。従って、この例で使用されるＱ及びＲの定数値について、カルマンゲインＫは０．５になる傾向を有し、状態共分散Ｖは５になる傾向を有する。しかし、遭遇されるプロセスのタイプ、例えば、健康な心拍数（その場合、低レベルのプロセスノイズが存在する）ではなく心房細動（その場合、高レベルのプロセスノイズが存在する）に従って、分散定数Ｑ及びＲ、好ましくはプロセス分散Ｑについて値を修正することによりＫを適応型にすることができる。
【００３１】
図８、図９、及び図１０は、図７に示される心拍数プロットのために（この場合、心電図測定チャネル）、推定された心拍数ｘｐｒｅｄ、イノベーションｅ、及び分散σ2の値を示す。
【００３２】
図３は、この実施の形態を使用して図１及び図２からの２つの測定プロットを組み合わせた結果を示す。図２における心電図チャネルの動きアーチファクト２２は、心房細動の期間に起こっても、組み合わせられた測定から除去されていることが分かる。
【００３３】
従って、この実施の形態において、測定とその予測値との間の差は、その測定における信頼の程度を示すために使用される。差が大きくなれば、それだけ信頼度は低くなる。前記の式（１）は、２つの測定を組み合わせるために使用される。これは、次のように要約することができる。
（１）双方のチャネルにおける有効な心拍数：  双方のチャネルにおける低いイノベーション；  双方の測定を等しく重みづける。
（２）双方のチャネルで発見される有効な突然変化（例えば、期外収縮）：  双方のチャネルにおける高いイノベーション；  従って、測定は再び等しく重みづけられる。
（３）１つのチャネルにおけるアーチファクト：  １つのチャネルだけにおける高いイノベーション；  従って、そのチャネルからの情報は、非常に低い重み（低い信頼度）を与えられることによって無視される。
【００３４】
この方法は、更に、動きアーチファクト検出器を提供するために使用されることができる。即ち、動きアーチファクトが双方のチャネルに存在し、それにより有用な情報を入手できない場合を検出するために、この方法を使用することができる。このような場合は、或る持続時間の間、双方のチャネルでイノベーションの高い値が存在することによって特徴づけられる。それを利用して、動きによって崩壊された生データ部分を廃棄し、それらの期間では、有効な心拍数の推定を導き出させないことを表示することができる。
【００３５】
（付録１）
load-ascii enghr_13
data_file = ecghr_13
time = data_file (:,1);
hr = data_file (:,2);
start = 1
stop = size (data_file), 1);
fprintf ('検出されたR-R間隔の数 = %d \n' , stop);
hr_limit = 200;
％  Ｘ（ｔ＋１）＝ＡＸ（ｔ）＋ノイズ（Ｑ）  −  ノイズｗの分散としてのＱを有するプロセスモデル
％  Ｙ（ｔ）＝ＣＸ（ｔ）＋ノイズ（Ｒ）  −  ノイズｖの分散としてのＲを有する測定モデル
ss = 1        ％  状態サイズ  −  １次元へ設定する。即ち、ルーチンはベクトルのために働くが、スカラーにする。
os = 1        ％  観測サイズ  −  １次元へ設定する。即ち、スカラーにする。
A = [1];      ％  ｘ（ｔ＋１）＝ｘ（ｔ）と仮定する。
C = [1];      ％  ｙ＝ｘと仮定する。
Q = 5.0 *eye (ss);   ％  プロセスノイズ  −  ｅｙｅはＭＡＴＬＡＢにおける恒等行列  −  ここでは単に１である。
R = 10.0 *eye (os)   ％  測定ノイズ分散
initx = [80];        ％  初期状態値（８０ｂｐｍの心拍数）
initV = 100 *eye (ss);      ％  初期状態分散
xnew = initx;        ％  初期化
Vnew = initV;
for i = start:stop - start of cycle             ％  サイクルのスタート
x = xnew;     ％  前のサイクルから更新する。
V = Vnew;     ％  前のサイクルから更新する。
xpred = A*x;         ％  状態の予測
Vpred = A*V*A' + Q;  ％  状態共分散の予測、Ａ’はＡの転置である。
ypred (i) = C*xpred; ％  測定の予測
y(i) = hr(i);        ％  「測定を行なう」
e = y(i) - ypread(i);      ％  イノベーションを計算する。
innov(i) = e;        ％  プロットのため
sigma2(i) = e*e;            ％  保存のための分散
S = C*Vpred*C' + R;  ％  イノベーションの共分散
Sinv = inv(S);              ％  Ｓを反転させる。
K = Vpred*C' *Sinv;  ％  カルマンゲインを計算する。
xnew = xpred + K*e;  ％  カルマンゲインによって制御されたイノベーションによって状態を更新する。
Vnew = (eye(ss) - K*C)*Vpred;     ％  状態共分散を更新する。
end
end
【図面の簡単な説明】
【図１】  パルスオキシメトリによって測定された心拍数のプロットを示す図である。
【図２】  心電図によって測定された心拍数のプロットを示す図である。
【図３】  本発明の実施形態に従って図１及び図２の２つのプロットを組み合わせた結果を示す図である。
【図４】  心電図上及びパルスオキシメトリ波形上の鼓動を概略的に示す図である。
【図５】  心電図及びパルスオキシメトリ波形に現れる期外収縮を概略的に示す図である。
【図６】  心電図トレース及びパルスオキシメトリ波形を、パルスオキシメトリ波形上のアーチファクトとともに示す図である。
【図７】  心電図によって測定された心拍数のプロットを示す図である。
【図８】  本発明の１つの実施の形態に従って心拍数の予測値を示す図である。
【図９】  図７及び図８から得られたイノベーションを示す図である。
【図１０】  図７及び図８から得られた分散を示す図である。
【図１】
0004741168.tif 000005
【図２】
0004741168.tif 000006
【図３】
0004741168.tif 000007
【図４】
0004741168.tif 000008
【図５】
0004741168.tif 000009
【図６】
0004741168.tif 000010
【図７】
0004741168.tif 000011
【図８】
0004741168.tif 000012
【図９】
0004741168.tif 000013
【図１０】
0004741168.tif 000014