JP2013500757A - 心弾道図信号解析の方法および装置 - Google Patents

Abstract

心弾道図信号の解析のための方法および装置が提供される。本方法は、BCG信号中のユーザーの心拍の典型的な特徴を位置特定することによってBCG信号中の心拍を検出することを含み、前記心拍の典型的な特徴はトレーニング段階の間に得られたものである。

Description

本発明は、心弾道図（ballistocardiogram）信号の解析のための方法および装置に、詳細には心弾道図信号における単独心拍イベントの検出のための方法および装置に関する。

心弾道計（BCG: ballistocardiograph）は、心臓によって拍出される際の血液の勢いに起因する人体の動きを測定する。

BCGは、体に電極を貼り付けたり、ベルト、布などのような特別なセンサーを身につけたりする必要なしに身体生命徴候〔バイタルサイン〕の測定が可能であるという点で、心電計（ECG: electrocardiograph）に対して利点がある。この邪魔にならない性質のため、BCGは、夜間、長時間にわたって人の心臓活動をモニタリングするのに最も適している。同様に、BCGシステムは、病院の一般病棟で患者の睡眠の質を下げることなく患者をモニタリングするために追加的な安全施策として使うことができる。現代のBCGシステムはベッドに完全統合でき、単一スイッチによって作動させられるので、そのようなシステムによって提供される追加的な安全性は、健康ケア専門家に要求する超過努力が最小限ですむ。

Kortelainen,J.M. and Virkkala,J.、"FFT averaging of multichannel BCG signals from bed mattress sensor to improve estimation of heart beat interval"、Engineering in Medicine and Biology Society, 2007, EMBS2007、29th Annual International Conference of the IEEE, 22-26 August 2007、pp.6685-6688 J.H.Shin, B.H.Choi, Y.G.Lim, D.U.Joeng and K.S.Part, "Automatic Ballistocardiogram (BCG) Beat Detection Using a Template Matching Approach"、Engineering in Medicine and Biology Society, 2008, EMBS2008、30th Annual International Conference of the IEEE, 21-24 August 2008 D.H.Phan, S.Bonnet, R.Guillemaud, E.Castelli, N.Y.Pham Thi、"Estimation of Respiratory Waveform and Heart Rate Using an Accelerometer"、Engineering in Medicine and Biology Society, 2008, EMBS2008、30th Annual International Conference of the IEEE, 21-24 August 2008 J.N.Edson, R.Dickes, G.H.Flamm and M.Tobin、"Higher Frequency Phenomena in the Normal Ballistocardiogram"、Circulation Research、1958、Vol.1、pp.405-409

現在のところ、心拍数を決定するために心弾道図信号を解析するためのアルゴリズムは、スペクトル法または時間領域法を使い、たとえば信号の自己相関関数を評価することによってある種のパターンの再発生を検出する。こうしたアプローチのすべてにおいて、複数の心拍をカバーするよう数秒間持続する信号セグメントを考慮する必要がある。結果として、ある時間期間にわたる平均心拍が得られるが、心拍から心拍の間の情報は得られない。

心弾道図信号からの心拍間推定のためのいくつかのアルゴリズムが提起されているが、これらは適正に動作するために大型で高価なセンサー・アレイ（非特許文献１）または人間の対話（非特許文献２）を必要とするか、あるいは異なる複数のセンサー・モダリティーを使い、精確さを欠く（非特許文献３）。

これらのアルゴリズムが市場に投入できるかどうか、あるいはこれらのアルゴリズムが、特に不整脈の患者について、心弾道図の患者内および患者間多様性に対処できるかどうかには疑問がある。

不整脈は広範な問題であり、危機的な健康状態の前兆となりうる。米国心臓協会によれば、推定220万人の米国人が心房細動を患って暮らしている。これは、いくつかのアルゴリズムにとっての規則的な心拍という前提が、たとえば病院の一般病棟においてモニタリングされるべき患者のかなりの部分について成り立たないということを意味している。

不整脈の深刻さに依存して、BCG信号における心拍を検出しようとする既存のアルゴリズムは遅かれ早かれ失敗する。通例、心臓がほぼ規則的に拍動していることを前提とするからである。よって、この規則性の前提に基づかず、よって深刻な不整脈の事例でも機能できる心弾道図信号処理アルゴリズムが必要とされている。BCG信号の信頼できる心拍間の解析を提供できるアルゴリズムも望ましい。

さらに、既存のBCG解析アルゴリズムの多くはオフラインで機能する。よって、望ましいアルゴリズムは、オンラインでの（すなわち連続的な）使用に好適である、よって病院、介護施設または自宅でのモニタリングおよび警告タスクに好適であるものである。

本発明の第一の側面によれば、心弾道図（BCG）信号中のユーザーの心拍を検出する方法であって、前記BCG信号中のユーザーの心拍の典型的な特徴を位置特定することによって前記BCG信号中の心拍を検出することを含み、前記心拍の典型的な特徴はトレーニング段階の間に得られたものである、方法が提供される。

本発明の第二の側面によれば、ユーザーの心弾道図信号を測定するデバイスとともに使うための装置であって、前記デバイスから心弾道図信号を受け取る手段と；受け取った心弾道図信号に対して上記方法を実行する処理手段とを有する装置が提供される。

本発明の第三の側面によれば、コンピュータまたはプロセッサ上で実行されたときに、前記コンピュータまたはプロセッサに上記方法を実行させるよう構成されたコンピュータ・プログラム・コードを含むコンピュータ・プログラム・プロダクトが提供される。

本発明について、単に例として、図面を参照しつつ、以下の詳細な説明において記述する。
本発明に基づく、心拍を検出する原理を例解するフローチャートである。 フィルタリングされていない心弾道図信号を示すグラフである。 本発明に基づく、アルゴリズムをトレーニングする方法を示すフローチャートである。 呼吸成分を除去するために心弾道図信号をフィルタリングした結果を示すグラフである。 フィルタリングされた心弾道図信号を平滑化した結果を示すグラフである。 本発明のある実施形態に基づく、心弾道図信号から導出されるパラメータの図解である。 図３のステップ１０７の実装をより詳細に示すフローチャートである。 心弾道図信号の高周波数フィルタ応答包絡線を示すグラフである。 心弾道図信号中の心拍を検出するアルゴリズムを使う方法を示すフローチャートである。 信頼性スコア三値組（triplet）を示すグラフである。 ECGを使って検出された心拍に対する、本アルゴリズムの有効性を示すグラフである。 本アルゴリズムによって検出された心拍から心拍までの間隔とECGを使って検出された心拍から心拍までの間隔を比較するグラフである。 ECGに対して、不整脈があるときの心拍検出における本アルゴリズムの有効性を示すグラフである。 本発明のある実施形態に基づく装置のブロック図である。

図１のフローチャートは、本発明に基づくアルゴリズムによって使用される一般的な原理を示している。第一段階（ステップ１０）であるトレーニング段階では、心弾道図（BCG）信号の短いセグメントが解析され、心拍の典型的な特徴が判別される。第二段階（ステップ１２）である割り当て段階では、ライブBCG信号の新たに記録されたサンプルが心拍を求めてスキャンされる。換言すれば、本アルゴリズムは、トレーニング段階で学習された特徴を、新たに記録された信号中に再発見しようとする。トレーニング段階と割り当て段階の動作について簡単に概観しておく。

トレーニング期間（ステップ１０）中に心拍の典型的な特徴を学習するために、トレーニングに使われる心弾道図のセグメントはまず、呼吸成分があればそれを除去し、該信号を平滑化するために、フィルタリングされる。結果として得られる信号中で諸特定点が同定され、それが信号をパラメータ化するために使用できる。各特定点において、後続のn個のパラメータが、特徴ベクトルを形成するために使われる。結果として得られる特徴ベクトルは、類似性によってクラスター化される。ここでの原理は、個々の心拍をカバーする特徴ベクトルはまとめてクラスター化されるのに十分類似しており、他の特徴ベクトルは異なるクラスターに見出されるであろうということである。

この処理とは別に、もとの（フィルタリングされていない）トレーニング・セグメントは、呼吸信号を除去するためにフィルタリングされるが、このときは平滑化は行われない。この信号において高周波数成分が決定され、これは通例、各心拍に一致する。このプロセスにおいて、高周波数成分を検出するために使用されるフィルタの通過帯域周波数は自動的に調整される。

心拍に関係する特徴ベクトルを含むクラスターは、決定された高周波数成分と一緒に個々のクラスターの特性を使って識別される。

割り当て段階（ステップ１２）の間に、新たに測定または受領された信号がやはり上記のようにしてフィルタリングされ、平滑化され、パラメータ化される。この信号のパラメータが、トレーニング段階（ステップ１０）において心拍の特徴を含むと識別されたクラスターと比較される。さらに、この新たに到着した信号において高周波数成分が検出される。心拍を含むクラスターとの類似性および高周波数成分の位置に基づいて、心拍の位置の第一の推定がなされる。第一の推定によって与えられる心拍の位置は、二つの連続する心拍信号の間の自己相関に基づく洗練段階において較正される。

本発明に基づくアルゴリズムの動作について、ここで付属の図面の図２ないし図１０を参照しつつより詳細に述べる。

図２は、低周波数の擬似正弦波信号として現れる呼吸成分と鋭いバーストにつながる心拍成分からなるBCG信号の典型的なセグメントを示している。提起されるアルゴリズムの目標は、心拍に関係したバーストを同定し、それらの間のパルスからパルスの（pulse-to-pulse）（心拍から心拍の（beat-to-beat）ともいう）時間距離を決定することである。この目的のため、上で要約したように、本アルゴリズムはまず、心拍を特徴付ける特徴を検出し、学習する。

これらの特徴は個人によって変わり、モニタリングされる人物のBCGセンサーに対する位置にも強く依存することは理解されるであろう。

〈学習段階〉
本アルゴリズムの学習段階（ステップ１０）が図３に描かれている。ステップ１００では、心弾道図信号のセグメントを含む、s_trainと表されるトレーニング信号が得られる。トレーニング・セグメントは、いくつかの心拍を含むほど十分長い必要があり、アーチファクトがない必要がある。このセグメントは好ましくは10秒から60秒の間の長さである。

ステップ１０１では、1Hzのカットオフ周波数をもつ高域通過フィルタによって、心弾道図信号から呼吸成分が除去される。結果として得られるフィルタリングされた信号s_filt（図４に描かれる）は、のちにさらに述べるように、二つの異なる処理経路において使用される。

本アルゴリズムの後続の諸ステップでは、上記の前処理された信号のすべての極値（すなわち極大または極小）が解析される。これらのステップにおいて、無視できるほど小さい極値が本アルゴリズムの成功を妨害することを考えつつ、有意な、よく現れた極値のみが考慮に入れられるべきである。この理由で、セグメントはまず、本アルゴリズムの第一の経路（ステップ１０２ないし１０６）では、10Hzのカットオフ周波数をもつ低域通過フィルタによって平滑化される（図３のフローチャートのステップ１０２）。図５は平滑化された信号s_smoothを示している。見て取れるように、図４においてまだ見えた高周波数ノイズが、よって小さな極値も除去されて、信号の基本的な形だけが残っている。

BCG信号s_smoothにおいて、呼吸成分が除去されたのちは、心拍は、類似した振幅およびピーク間距離をもつピークのグループとして見分けられる。

ステップ１０３では、フィルタリングされ、平滑化されたBCG信号の特徴点が検出される。ある好ましい実施形態では、これらの特徴点は信号の極大および極小である。

ステップ１０４では、信号s_smooth中の検出された各極大がパラメータ化される。ある好ましい実施形態では、各極大kは次の特徴を使ってパラメータ化される。
ｉ）極大の振幅（a_max）；
ｉｉ）その極大と右側の次の極小との間の距離（d_max）；
ｉｉｉ）その極大の右側の次の極小の振幅（a_min）；および
ｉｖ）その極大の右側の次の極小と、右側の次の極大との間の距離（d_min）。

これらのパラメータは図６に示されており、これらは信号の形のおおまかな記述を与える。

ステップ１０５では、n個（nは通例n＝3ないし10の範囲である）の連続するピークのパラメータが集められ、特徴ベクトルf_kを形成するために使われる。特徴ベクトルf_kは次式によって与えられる。

f_k＝[a_max,k, d_max,k, a_min,k, d_min,k, a_max,k+1, a_min,k+1, d_min,k+1,…,
a_max,k+n-1, d_max,k+n-1, a_min,k+1-1, d_min,k+n-1] (1)
ここで、k＝1,2,…,K−n＋1であり、Kはトレーニング・セグメントs_trainにおけるピークの数を表す。各特徴ベクトルf_k中に集められるピークの数nは、ほぼBCG信号中の各心拍の間に現れるピークの数に対応する。よって、ステップ１０５はK−n＋1個の相異なる特徴ベクトルf_kを与える。

当業者は、BCG信号s_smoothの他のパラメータ化が可能であることを理解するであろう。たとえば、信号は、スプラインによって近似されることができ、該スプライン近似の節点が極大および極小ならびに可能性としてはそれらの中間の少数の（通例一つだけの）他の点である。あるいはまた、信号の最も有意な諸フーリエ係数に関する該信号のパラメータ化も可能である。

ステップ１０６では、特徴ベクトルf_kは類似性によってクラスター化される。クラスタリング方法は通例高次元ベクトル空間より低次元ベクトル空間で最もよく機能するので、本発明のいくつかの実施形態では、これに先立って、主成分解析によって特徴ベクトルの次元を低下させることができる。

BCG信号s_smoothにおいて、心拍はピークの反復パターンとして見える。したがって、心拍を記述する特定の特徴ベクトルは類似するはずであり、心拍の中間の信号のランダムなゆらぎを記述する特徴ベクトルとは顕著に区別されるはずである。

したがって、ステップ１０６において、類似性による特徴ベクトルのグループ化は、心拍のみに関係した特徴ベクトルを含むクラスター（単数または複数）につながる。

当業者は、特徴ベクトルf_kをクラスターにグループ化するためには、さまざまな方法、たとえばｋ平均クラスタリング、階層的クラスタリング、自己組織化マップなどが好適であることを認識するであろう。また、クラスタリング段階において、さまざまな距離指標d、たとえばユークリッド距離、個々の特徴ベクトル間の角度などが使用できることも理解されるであろう。

本発明のある好ましい実施形態では、標準的なクラスタリング方法とは別に、各特徴ベクトルについて、毎回、ある距離指標（好ましくはユークリッド距離）に関して次のP個の近隣（Pは10程度）を含むクラスターを形成することが有利である。

いずれのクラスタリング技法も、ある数M個のクラスターにつながり、心拍パターンを記述する特徴ベクトルを含むクラスターを同定することが必要である（ステップ１０７）。

図７のフローチャートがステップ１０７で実行される段階を示している。

ステップ１０７１では、M個のクラスターのそれぞれについて、クラスター中心f_cmが決定される（ここで、mはm番目のクラスターを示し、m＝1,…,M）。クラスター中心を決定する技法は当業者にはよく知られており、本稿ではさらに述べることはしない。

ステップ１０７２では、ステップ１０７１において同定されたM個のクラスター中心のそれぞれについて、クラスター中心f_cmに最も近い、集合{f_k|k＝1,2,…,K−n＋1}内の特徴ベクトルf_kが同定される。この特徴ベクトルは、m番目のクラスターによって表される信号部分の典型（archetype）またはモデルである特徴ベクトルなので、f_archmと表される。よって、心拍信号を表すクラスターについては、典型特徴ベクトルf_archまたはそのクラスターについてのクラスター中心f_cは、パラメータ化された心拍についてのモデルと見ることができる。

次に、ステップ１０７３において、f_archmに最もよく対応するサブセグメントs_archmを求めてトレーニング信号s_smoothが探索される。有効的には、特徴ベクトルf_archmについて、学習手順のステップ１０３および１０４が逆転されることになる。よって、位置特定されたサブセグメントs_archmは、f_archmに最も近いパラメータ化をもつs_smoothのn個のピークを含む。

以下に述べるステップ１０７４〜１０７８は、ステップ１０７２で同定されたM個のサブセグメントs_archおよびM個のクラスター中心から、心拍信号を表すクラスターおよびクラスター中心がどのようにして同定されるかを示している。以下の段階ではf_cが使われているものの、代わりにf_archが使われることができることは理解されるであろう。

まず、ステップ１０７４において、各クラスター中心f_cmと個々の特徴ベクトルf_kとの間の距離関数が計算される。好ましくは、クラスタリング段階（ステップ１０６）において使われたのと同じ距離指標が使われる。

この距離関数dは、下記の条件
d(f_c,f_p-1)＞d(f_c,f_p)＜d(f_c,f_p+1)
が満たされるとき、f_cmに関して、位置pに極小をもつことになる。ここで、p∈[1,2,…,K−n＋1]である。

f_cmと、そのクラスターに関係したまたはそのクラスターの一部である任意の特徴ベクトルとの間の距離関数dが極小をもつことは理解されるであろう。

次いで、ステップ１０７５において、各サブセグメントs_archmとトレーニング信号s_smoothとの間の相互相関が決定される。相互相関の結果は、サブセグメントs_archmがトレーニング信号s_smoothと最も類似するところ（すなわち、m番目のクラスター中の特徴ベクトルが表している信号部分）に最大をもつことは理解されるであろう。

心拍を表すクラスターおよびクラスター中心（これらはこの段階ではM個のクラスターのうちからまだ同定されていない）を考えると、BCG信号s_smooth内の各心拍のところにおいて、距離関数における極小および相互相関における極大があるであろう。

このように、ステップ１０７６では、距離関数の極小および相互相関関数の極大を使って、心拍位置の第一の推定が、M個のクラスターおよびクラスター中心のそれぞれについて別個に決定される。位置を同定するためにこれら二つの関数がどのように使われるかについてのさらなる詳細は、下記の割り当てセクションにおいて与えられる。

心拍信号に関係するクラスターおよびクラスター中心が心拍の位置を与えることに加えて、他のクラスターも距離関数における極小および相互相関関数における極大をもつ、s_smooth内の位置を同定するであろう。しかしながら、これらの極値は一般に、BCG信号において散発的に現れるのみである。

よって、BCG信号中の二つの個々の極値の間の時間を検査することによって、心拍に関係しないクラスターが識別できる（ステップ１０７７）。連続する「心拍」どうしの間の間隔がたとえば3秒より大きいおよび／または0.25秒より小さい場合、対応するクラスターは心拍に関係していないと想定でき、さらなる解析において無視できる。

上記のステップ１０７７における検査は、クラスターの多くを解析から排除することにおいて有用であるが、二つ以上のクラスターがこの検査に合格することもしばしばである。よって、さらなるノックアウト基準が必要とされる。

心拍がBCGにおける比較的高い周波数をもつセグメントにつながるというのは一般的な観察である。非特許文献４では、これらの高周波数成分が、心臓弁の開閉によって引き起こされる機械的な力に関係していると言われている。これらの成分の周波数は個々に、またBCGセンサーの位置にも依存して変動しうるので、この方法においてそれらを検出するためには、動的に調整される帯域通過（bandpass）フィルタが使用される（図３のステップ１１０参照）。

最適パラメータは、信号セグメントs_filt（図３のステップ１０１より）を、さまざまな帯域通過フィルタを用いてフィルタリングすることによって決定される。帯域通過フィルタのセットは、個々のフィルタがそれぞれ2Hzの帯域幅をもち、セット全体が4Hzから20Hzまでの周波数範囲を、0.1〜0.5Hzの増分でカバーするよう構築される。よって、帯域通過フィルタのセットは、カットオフ周波数＝{[4Hz,6Hz],[4.1Hz,6.1Hz],[4.2Hz,6.2Hz],…,[20Hz,22Hz]}をもつ諸フィルタを含むことができる。

上記セット中の各フィルタについて、フィルタ応答が二乗され、3.5Hzのカットオフ周波数をもつ低域通過フィルタが、（図８に示すような）フィルタ応答包絡線を得るために適用される。

フィルタ応答包絡線のそれぞれについて、その包絡線における各ピークに、
ａ）そのピークの振幅と、平均ピーク振幅にピーク振幅標準偏差を加えたものとの比、および
ｂ）ピークの対称性
によって、有意性（relevance）スコアが割り当てられる。

次いで、最高の平均有意性スコアにつながるフィルタが選択され、ステップ１１０において、フィルタリングされたBCG信号s_filtの高周波数成分を検出するために使われる。有意性スコア計算のより詳細な記述は下記の割り当てセクションで与えられる。

最適なフィルタを使ってステップ１１０で検出された高周波数成分は、今や、ステップ１０７において心拍に関係したクラスターを同定するのを助けることができる。ここで図７、特にステップ１０７８に戻ると、距離関数における極小および相互相関関数における極大がステップ１１０において決定された高周波数成分の位置と最もよく一致するクラスターが、心拍に最もよく関係していると考えられる。さらに、距離関数の極小の深さ、相互相関関数の極大の振幅ならびに高周波数成分の最大の振幅が、心拍（heart beat）パターンを最もよく記述するクラスターを決定するために考慮に入れられることができる。このクラスターは本稿ではc_HBと記される。この段階での解析に残っていた他のすべてのクラスターは今や破棄される。

このクラスターのクラスター中心f_cは今やc_HBcentreと表され（図３のステップ１０８）、c_HBcentreに最も近く一致するs_smooth中のセグメントが同定される。このセグメントはs_HBarchと記される（ステップ１０９）。これでアルゴリズムの学習部分が終わる。

〈割り当て段階〉
ひとたび中心c_HBcentreおよび心拍信号の対応する典型s_HBarchが決定されると、新しくはいってくるBCG信号（オンラインの場合）またはBCG信号の残り（オフラインの場合）が心拍を求めてスキャンされることができる（図１のステップ１２）。以下では、オンラインの場合について述べる。当業者は、同様の考察がオフラインの場合にも成り立つことを理解するであろう。

図９は、割り当て段階の間に使われる方法を示している。

BCGセンサーから新しいサンプル（s_real）がはいってくるときは常に、図３に示したトレーニング・アルゴリズムのステップ１０１、１０２と同様に、フィルタリングされ（ステップ２０１）、平滑化される（ステップ２０２）。同様に、フィルタリングされ、平滑化された信号における特徴点（すなわち、極大および極小）が決定される（ステップ１０３に対応するステップ２０３）。この信号は次いでこれらの特徴点においてパラメータ化され、特徴ベクトルが構成される（ステップ２０４、２０５）。

ステップ２０６（これは図７のステップ１０７４に対応する）では、特徴ベクトルと心拍クラスター中心c_HBcentreとの間の距離が計算され、極小を求めてスキャンされる。同様に、s_HBarchと新しくはいってくる信号（フィルタリングおよび平滑化後）との間の相互相関およびその極大が決定される（ステップ２０７）。

さらに、学習アルゴリズムのステップ１１０において決定された最適帯域通過フィルタを適用して、高周波数包絡線における極大を求めて、新しいBCG信号（s_filt）をスキャンする（ステップ２０８）。

ステップ２０６、２０７および２０８における解析の結果は、BCG信号s_real中のどこに心拍が現れるかを識別するために使うことができる。しかしながら、距離関数の極小、相互相関の極大および高周波数成分の極大はやや異なる位置（時刻）に現れうることを注意しておく。さらに、偽の高周波数成分のみならず距離および相互相関関数の誤った極値が現れることがある。そのような場合、個々の基準の間の合理的な妥協を見出すためには、個々の基準の信頼性を評価することが重要である。よって、ステップ２０９では、ステップ２０６、２０７および２０８において見出された極小および極大の信頼性が評価される。

距離関数の極小の信頼性は、二つの隣り合う極大の高さに対するその深さによって評価される。この目的のために、左の極大と当該極小との間の振幅の差（h_l）および右の極大と当該極小の間の振幅の差（h_r）が計算される。すると、振幅a_minをもつ極小の信頼性r_dは
r_d＝(h_l＋h_r)／(2(a_min＋max(h_r,h_l)))
として評価される。

構成により、r_dは規格化されている（すなわち0≦r_d≦1）ことが示せる。最大スコア1は、隣接する両ピークの高さに対して対称的であり、振幅0をもつ極小に割り当てられる。

同様に、振幅a_maxをもつ相互相関ピークは
r_x＝(h_l＋h_r)／a_max
を使って評価される。ここで、h_lおよびh_rはそれぞれ当該ピークと左および右の極小との間の振幅差である。r_dと同様、r_xも規格化されており、最大スコアは、左および右の極小より上の高さに関して対称的であるピークに割り当てられる。

最後に、信号の高周波数成分におけるピークの信頼性は、
r_hf＝(h_l＋h_r)／(2a_max)・(a_max/a_ref)
を使って評価される。

この場合、公式の第一の因子はr_xの計算と同一であり、よってピークの相対的な対称性を定量化する。一方、第二の因子はピーク振幅a_maxと参照振幅a_refとの比によってスコアを調整する。参照振幅a_refはトレーニング・シーケンスの間に、高周波数ピークの平均振幅に振幅の標準偏差を加えたものとして計算される。

最終的に、どの時点で心拍が実際に起こったかを判別するために、上記で計算された信頼性r_x、r_dおよびr_hfについての情報がどのようにして有利に組み合わされることができるかを判別しなければならない（ステップ２１０）。相互相関基準、距離基準および高周波数成分基準がそれぞれ独自に心拍を検出する時点は、それぞれt_x、t_dおよびt_hfと表される。理想的な場合では、これら三つの時点は同一である（すなわちt_x＝t_d＝t_hf）。しかしながら、現実世界の場面では、これらは、BCG信号の異なる特徴に焦点を当てることによって心拍を検出しようとしているので、互いにわずかに異なっている。心拍が三つの基準（距離、相互相関および高周波数成分）すべてによって同時によく検出される場合、時刻t_d、t_xおよびt_hfは明確に見分けられる三値組をなすであろう（図１０参照）。以下の手順のねらいは、これらの三値組を同定し、該三値組を、個々の心拍を同定するために使うことのできる代表的な値によって置き換えることである。

一般性を失うことなく、時刻t_hfを本アルゴリズムのさらなる記述のための出発点として考える。高周波数基準はこの時点t_hfにおいて心拍を示唆する。その場合、本アルゴリズムはt_xおよびt_dが区間[t_hf−c,t_hf＋c]にあるかどうかを検出する。ここで、cは定数であり、通例約0.3秒に設定される。t_xおよびt_dがこの区間内にあれば、三つの時刻t_x、t_dおよびt_hfとおよび三つの信頼性r_x、r_dおよびr_hfの代表値R_hf(r^* _hf,t^* _hf)が構築される。これは時刻t^* _hfおよび信頼性r^* _hfによって定義される。代表時刻t^* _hfは対応する信頼性r_x、r_dおよびr_hfで重みをかけた時刻t_x、t_dおよびt_hfの和として、次のように計算される。

t^* _hf＝(r_xt_x＋r_dt_d＋r_hft_hf)／(r_x＋r_d＋r_hf)
代表信頼性は単に、相互相関、距離および高周波数基準の信頼性の和として計算される。すなわち、
r^* _hf＝(r_x＋r_d＋r_hf)／N N＝3
となる。

これまでは、距離および相互相関基準はそれぞれ区間[t_hf−c,t_hf＋c]内の唯一の時刻t_dおよびt_xにつながると想定してきた。もしこれらの時刻の一方しかこの区間内に見出されない、あるいはどちらもこの区間内に見出されない場合でも、本手順は上記のとおりに機能する。ただし、上記の公式において欠けている時点の信頼性は0に設定される。したがって、Nの値は0でない信頼性の数に設定される。

同様に、距離および／または相関基準が区間[t_hf−c,t_hf＋c]内の二つ以上の振幅を示唆することもありうる。距離基準が区間[t_hf−c,t_hf＋c]内の時刻t_d1およびt_d2に二つの振幅を示すとすると、対応する信頼性r_d1およびr_d2は、t_hfを中心とし、t_hfからの距離とともに線形もしくは非線形に減少する対称的な窓関数によって重みをかけられる。ある好ましい実施形態では、区間[t_hf−c,t_hf＋c]の両端で（ほとんど）0にまで減少するガウス窓関数が使われる。重みをかけられた信頼性が大きいほうの時点t_d1またはt_d2のみが、本アルゴリズムの後続の諸ステップにおいて考慮され、他方は破棄される。

これまでは、上記の手順は出発点としての時点t_hfに基づいていた。この出発点は任意に選んだのであり、好ましくは、実際上、上で概説した手順は出発点としてのt_dおよびt_xに関してそれぞれ繰り返される。これは新しい代表R_d(t^* _d,r^* _d)およびR_x(t^* _x,r^* _x)につながる。図１０に示されるような明確に見分けられる三値組の場合には、代表t^* _d、t^* _xおよびt^* _hfの時間成分は同一である。

心拍が起こったと考えられる各時点tについて、本アルゴリズムはスコアを計算する。時刻tに代表が見出されない（すなわち、どの代表も、tに等しい時間成分t^* _d、t^* _xまたはt^* _hfをもたない）場合、スコアは0である。それ以外の場合には、スコアは、tに等しい時間成分をもつ代表の信頼性成分の和に等しい。

たとえば、t^* _d＝t^* _x＝t^* _hf＝tで、三つの代表R_d(t^* _d,r^* _d)、R_x(t^* _x,r^* _x)およびR_hf(t^* _hf,r^* _hf)がある場合、時刻tでの対応するスコアS(t)は
S(t)＝r^* _d＋r^* _x＋r^* _hf
に等しい。

S(t)が所定の閾値S_threshより大きい場合には，本アルゴリズムは時刻tにおいて心拍を検出する。そうでない場合には、心拍は検出されず、本アルゴリズムは次のサンプルを待つ（オンラインの場合）。当業者は、閾値S_threshがトレーニングおよび／または割り当て段階の間に、最後に検出されたq個の心拍の信頼性値に依存して調整されることができることを理解するであろう。

さらに、モニタリングされる人物がかなり規則的な心拍を示す場合、最近決定された心拍から心拍までの間隔が、次の心拍がいつ起こる可能性が最も高いかを予測するために使われることができる。これらの時点付近では、心拍を見逃す数を減らすために、閾値S_threshは低くされる。

上記で同定された点はすでに、BCG信号において個々の心拍がどこで起こったかをよく示している。しかしながら、個々の心拍から心拍までの時間間隔をより精密に決定するためには、さらなる洗練ステップが必要とされる。このステップでは、ピークのパターンが信号中で初めて反復する時間期間が微調整される。

具体的には、Nをこれまで見出された心拍の数として、t₁,t₂,…,t_Nを心拍が検出された時点であるとする。さらに、1≦p≦N−1として、s(t_p,t_p+1)が時点t_pとt_p+1の間に記録された信号s_filtのセグメントを表すとする。

信号中でp番目のピークがいつ反復するかを微調整するために、セグメントs(t_p,t_p+1)と次のセグメントs(t_p+1＋dt,2t_p+1−t_p−1＋dt)との間の相互相関xcorrを最大化する引数t_shift
t_shift＝argmax_dt xcorr(s(t_p,t_p+1), s(t_p+1＋dt,2t_p+1−t_p＋dt)) dt∈[−0.15s,0.15s]
が決定される。p番目のピークとp＋1番目のピークとの間の心拍から心拍までの間隔（interval）i_p,p+1は最終的には
i_p,p+1＝t_p+1−t_p＋dt
として計算できる。

このステップ後、心拍の検出および心拍から心拍までの間隔の計算が終了し、本アルゴリズムは次のサンプルが到着するのを待つ。

〈結果〉
図１１は、本発明に基づくアルゴリズムを、心拍がかなり規則的に現れる、すなわち不整脈が存在しないBCGデータの短いセグメントに適用したときに、該アルゴリズムによって得られた結果を示すグラフである。参照として、BCG信号と同時に得られたECG信号がプロットされている。二つのプロットを比べることによって、本アルゴリズムがBCG信号中の個々の心拍を明確に同定していることが見て取れる。

さらに、図１２は約８分にわたって同時に記録されたBCG信号とECG信号の解析を示している。それぞれの新しい心拍について、直前の心拍までの時間距離がECGにおいて決定される。同様に、本発明に基づくアルゴリズムが、BCG信号に基づいて対応する心拍から心拍の距離を計算するために使われる。ECGおよびBCG解析両方の対応する心拍から心拍までの距離は対にされ、図１２に示される散布図にプロットされている。BCGおよびECG解析の両方が同一の推定につながるのであれば、図１２のプロットはxy対角線（破線で示される）に沿ったエントリーのみを含むはずである。この対角線は明らかに図１２において優勢であり、非常に少数の逸脱する対のみが見出せる。セグメント全体にわたって平均されると、ECGベースのピークからピークまでの間隔とBCG間隔との間の絶対的な偏差は7msであった。全体として、すべてのピークのうち見逃されたのは0.6%だけであり、0.2%の擬陽性が検出された。

最後に、図１３は非周期的に拍動する心臓のECG信号および対応するBCG信号の両方を示している。規則的に拍動する心臓の場合（図１１）と同様、本発明に基づくアルゴリズムはこのセグメントにおけるすべての心拍を検出する。これは、心拍が規則的に拍動するという前提に基づくアルゴリズムでは可能ではない。

最新技術のBCGセンサーが、モニタリングされるべき人物のベッド中に見えないように統合でき、よって心拍数、また呼吸数の全く邪魔にならないモニタリングを提供できることが当業者には理解されるであろう。さらに、BCGベースのモニタリング・システムは、確立された集中治療室（ICU: intensive care unit）モニタリング・システムよりもずっと安価である。これら二つの特徴のため、BCGベースのソリューションは、患者がもはやフルスケールのICUモニタリング・システムは必要としないが、医師がまだ患者の生命機能に関心がある病院の一般病棟にとって、理想的な解決策になる。特に、この場面において、不整脈はしばしば他の深刻な健康問題の指標または前兆となりうるので、不整脈の正確な検出は重要である。

BCGベースの技術は、心拍数および呼吸数の長期のモニタリングのために家庭でも使用されることができる。本願は、たとえば、心不全を患う患者にとって多大な利益になるものである。その場合、心拍数変動の低下が代償不全の前兆と見られるのからである。

さらに、特にその邪魔にならない性質のため、BCGは、夜間に、モニタリングされる人物を乱すことなく、睡眠の質を評価するために使用できる。

本発明は方法またはアルゴリズムの面で記載してきたが、本発明はBCGシステム（すなわち、BCG信号を測定するための装置と組み合わせたコンピュータ装置）においてまたはスタンドアローンのコンピュータ・システムまたはプログラムとして実装されることができることは理解されるであろう。BCGシステムは、アナログまたはデジタル形式の心弾道図信号を本発明の装置に提供することができる、本発明の装置はこの信号をしかるべく受信するよう適応されることができることは理解されるであろう。たとえば、BCGシステムは心弾道図信号をアナログ形式で当該装置に提供でき、当該装置は、心弾道図信号のデジタル表現を当該装置内の好適にプログラムされたデジタル信号プロセッサに提供するために、アンチエイリアス・フィルタおよびアナログ‐デジタル変換器を有することができる。あるいはまた、心弾道図信号が当該装置に（特に当該装置内のデジタル信号プロセッサに）デジタル形式で提供されるよう、BCGシステムがアナログ‐デジタル変換器を実装することもできる。当該装置は、BCGシステムへの有線または無線接続など、いかなる適切な手段を使って心弾道図信号を受信することもできる。

本発明を実装するための装置の一つの実施形態が図１４に示されている。心弾道図信号はBCGセンサー３０２から当該装置３０４に提供される。装置３０４はBCG信号を入力ポート３０６において受信し、BCG信号をプロセッサ３０８を使って、上記の記述で述べたように処理する。プロセッサ３０８に当該方法を実行させるための命令はメモリ３１０に記憶されることができる。

したがって、心弾道図信号における単独心拍イベントを検出するための改善された方法および装置が提供される。

本発明は、図面および上記の記述において詳細に図示され、説明されてきたが、そのような図示および説明は、制約するものではなく、例解または例示するものと考えられるものである。本発明は開示される実施形態に限定されるものではない。

図面、本開示および付属の請求項を吟味することから、特許請求される発明を実施する当業者は、開示される実施形態への変形を理解し、実施することができる。請求項において、「有する」「含む」の語は他の要素やステップを排除するものではなく、単数形の表現は複数を排除するものではない。

単一のプロセッサまたは他のユニットが請求項に記載されるいくつかの項目の機能を充足してもよい。ある種の施策が互いに異なる従属請求項に記載されているというだけの事実が、それらの施策の組み合わせが有利に使用できないことを示すものではない。請求項において参照符号があったとしても、範囲を限定するものと解釈すべきではない。コンピュータ・プログラムは、他のハードウェアと一緒にまたは他のハードウェアの一部として供給される、光記憶媒体または半導体媒体といった好適な媒体上で記憶／配布されてもよいが、インターネットまたは他の有線または無線の遠隔通信システムを介してなど、他の形で配布されてもよい。

本発明の第一の側面によれば、心弾道図（BCG）信号中のユーザーの心拍を検出する方法であって、請求項１記載の諸段階を含む方法が提供される。

Claims (15)

1. 心弾道図（BCG）信号中のユーザーの心拍を検出する方法であって：
   前記BCG信号においてユーザーの心拍の典型的な特徴を位置特定することによって前記BCG信号中の心拍を検出する段階を含み、前記心拍の典型的な特徴はトレーニング段階の間に得られたものである、
   方法。
2. 請求項１記載の方法であって、前記心拍の典型的な特徴は、典型的な心拍のモデル特徴ベクトルc_HBcentreを含み、心拍を検出する前記段階は：
   前記BCG信号中に特徴点を識別する段階と；
   識別された各特徴点において前記BCG信号のパラメータを決定する段階と；
   決定されたパラメータから複数の特徴ベクトルを構成する段階と；
   前記複数の特徴ベクトルおよびc_HBcentreを使って前記BCG信号中の心拍を検出する段階とを含む、
   方法。
3. 請求項２記載の方法であって、前記心拍の典型的な特徴はさらに、前記トレーニング段階において使われたBCG信号の、前記モデル特徴ベクトルに対応する部分s_HBarchを含み、前記複数の特徴ベクトルおよびc_HBcentreを使って前記BCG信号中の心拍を検出する前記段階は：
   前記複数の特徴ベクトル、c_HBcentreおよびs_HBarchを使って前記BCG信号中の心拍を検出する段階を含む、
   方法。
4. 請求項３記載の方法であって、前記複数の特徴ベクトル、c_HBcentreおよびs_HBarchを使って前記BCG信号中の心拍を検出する前記段階は：
   前記複数の特徴ベクトルのそれぞれとc_HBcentreとの間の距離を計算する段階（２０６）と；
   結果として得られる複数の距離において極小を識別する段階（２０６）と；
   s_HBarchと前記BCG信号との間の相互相関を決定する段階（２０７）と；
   前記相互相関において極大を識別する段階（２０７）と；
   前記複数の距離における極小および前記相互相関における極大の位置から前記BCG信号中の心拍を識別する段階（２０８、２０９、２１０）とを含む、
   方法。
5. 請求項４記載の方法であって、心拍を検出する前記段階がさらに：
   前記BCG信号をフィルタ処理して高周波数成分包絡線を得る段階（２０８）と；
   前記高周波数成分包絡線において極大を識別する段階（２０８）とを含み、
   前記BCG信号中の心拍を識別する前記段階が：
   前記複数の距離における極小、前記相互相関における極大および前記高周波数成分包絡線における極大の位置から、前記BCG信号中の心拍を識別する段階（２０９、２１０）を含む、
   方法。
6. 請求項５記載の方法であって、前記BCG信号中の心拍を識別する前記段階がさらに：
   前記複数の距離における極小、前記相互相関における極大および前記高周波数成分包絡線における極大のそれぞれの信頼性を評価する段階を含む、
   方法。
7. 請求項６記載の方法であって、前記BCG信号中の心拍を識別する前記段階がさらに：
   前記複数の距離における極小、前記相互相関における極大および前記高周波数成分包絡線における極大から三値組を形成する段階（２１０）と；
   評価された信頼性から各三値組について、前記BCG信号において心拍が起こった時刻を示す代表値を決定する段階（２１０）とを含む、
   方法。
8. 請求項１ないし７のうちいずれか一項記載の方法であって、心拍の典型的な特徴がBCG信号s_trainの一部から得られるトレーニング段階をさらに含み、前記トレーニング段階は：
   s_train中に特徴点を識別する段階と；
   識別された各特徴点においてs_trainのパラメータを決定する段階と；
   決定されたパラメータから複数の特徴ベクトルを構成する段階と；
   前記複数の特徴ベクトルをそれらの特徴ベクトルの類似性に従って複数のクラスターにグループ分けする段階と；
   心拍に関係した前記複数のクラスター中のクラスターを識別する段階と；
   前記クラスター中の特徴ベクトルから、典型的な心拍についてのモデル特徴ベクトルc_HBcentreを決定し、前記モデル特徴ベクトルに対応するBCG信号の部分s_HBarchを決定することによって心拍の典型的な特徴を得る段階とを含む、
   方法。
9. 請求項８記載の方法であって、心拍に関係した前記複数のクラスター中のクラスターを識別する前記段階が：
   前記複数のクラスター中の各クラスターについて：
   クラスター中心f_cmを決定する段階（１０７１）と；
   前記クラスター中心f_cmに最も近く一致する特徴ベクトルf_archmを同定する段階（１０７２）と；
   前記特徴ベクトルf_archmに対応する前記BCG信号s_trainの部分s_archmを位置特定する段階（１０７３）と；
   前記クラスター中心f_cmと前記クラスター中の各特徴ベクトルとの間の距離関数を計算する段階（１０７４）と；
   s_archmの前記BCG信号との相互相関を決定する段階（１０７５）と；
   前記距離関数についての極小および前記相互相関についての極大をもつ前記クラスター中の特徴ベクトルから前記BCG信号中の心拍位置の推定を決定する段階（１０７６）と；
   前記BCG信号strainをフィルタ処理して高周波数成分の位置を同定する段階と；
   心拍に関係したクラスターを：
   （ｉ）指定された時間窓内にはいる量だけ離間している心拍位置の諸推定（１０７７）；および
   （ｉｉ）前記高周波数成分の位置と最もよく一致する、前記距離関数中の極小および前記相互相関中の極大（１０７８）
   をもつクラスターとして同定する段階とを有する、
   方法。
10. 請求項２ないし９のうちいずれか一項記載の方法であって、前記BCG信号中の前記特徴点が前記BCG信号中の極大を含む、方法。
11. 請求項１０記載の方法であって、前記BCG信号中の各極大点の前記パラメータが：
    ｉ）前記極大の振幅（a_max）；
    ｉｉ）その極大と右側の次の極小との間の距離（d_max）；
    ｉｉｉ）その極大の右側の次の極小の振幅（a_min）；および
    ｉｖ）その極大の右側の次の極小と、右側の次の極大との間の距離（d_min）
    を含む、方法。
12. 請求項２ないし１１のうちいずれか一項記載の方法であって：
    主成分解析を使って前記特徴ベクトルの次元を減らす段階をさらに含む、
    方法。
13. 請求項１ないし１２のうちいずれか一項記載の方法であって：
    N個の心拍がそれぞれ時刻t₁,t₂,…,t_Nに検出されたとして、1≦p≦N−1として、時刻t_pとt_p+1の二つの検出された心拍の間の前記BCG信号のセグメントと、時刻(t_p+1＋dt)と時刻(2t_p+1−t_p＋dt)との間の前記BCG信号ののちのセグメントとの間の相互相関を最大にするようなパラメータdtの値を同定することによって、前記BCG信号中に検出された心拍の位置を洗練する段階をさらに含む、
    方法。
14. ユーザーの心弾道図信号を測定するデバイスとともに使うための装置であって：
    前記デバイスから心弾道図信号を受け取る手段と；
    受け取った心弾道図信号に対して請求項１ないし１３のうちいずれか一項記載の方法を実行する処理手段とを有する、
    装置。
15. コンピュータまたはプロセッサ上で実行されたときに、前記コンピュータまたはプロセッサに請求項１ないし１３のうちいずれか一項記載の方法を実行させるよう構成されたコンピュータ・プログラム・コードを含むコンピュータ・プログラム。

Images