WO2017090732A1 - 呼吸推定方法および装置 - Google Patents

Abstract

呼吸推定方法は、被験者の心機能に関する第１のデータの時系列信号と被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から被験者の呼吸数をそれぞれ推定する第１のステップ（Ｓ１０７）と、第１のデータから推定した呼吸数と第２のデータから推定した呼吸数の各々についてカルマンフィルタによりノイズを濾した呼吸数をそれぞれ推定する第２のステップ（Ｓ１０８）と、この第２のステップで得られた複数の呼吸数の推定値の重み付け平均化処理を実行する第３のステップ（Ｓ１０９）とを含む。

Description

呼吸推定方法および装置

　本発明は、生体情報センサを用いて被験者の呼吸状態を推定する呼吸推定方法および装置に関するものである。

　人の呼吸を継続的に観察する呼吸連続モニタリングは、呼吸器系疾患の診断や呼吸器系疾患患者のリハビリテーション、胎児や患者の遠隔監視、ストレス診断等への応用が期待されている。一般的な呼吸計測の方法としては、マスクによる呼吸フロー計、鼻腔に設置したサーミスターによる熱計測、胸部での歪み計測バンド、胸部の電流インピーダンス計などが用いられるが、これらの方法には、被験者に装置装着時の違和感を与えるという問題がある。

　一方、近年、ウェアラブル機器の進展により、人体への装着感がより軽減されたデバイスが開発されている。例えば、図２２に示すようなシャツ等の衣類１００に電極１０１，１０２を埋め込んだウェアラブルセンサが提案されている（特許文献１）。このウェアラブルセンサは、導電性繊維からなる電極１０１，１０２を心臓１０３の近くに配置し、配線を介して電極１０１，１０２と端末（不図示）とを接続することにより、ウェアラブル心電計として機能する。

　以前より、心電位から人の呼吸の状態を推定することは試みられているが、心電位のＲ波と１つ前のＲ波の間隔であるＲＲ間隔を用いて被験者の呼吸の状態を推定しようとした場合、被験者の自律神経系の影響が顕著に現れ、測定時や測定前後の被験者の精神状態や被験者の年齢により推定結果が変化する。また、心電位のＲ波振幅を用いて被験者の呼吸の状態を推定しようとした場合、被験者の体動や皮膚状態による電極と皮膚との接触インピーダンスの変化、および被験者の体格や皮膚状態の個人差によって測定誤差が生じる。

　このようなことから、例えば非特許文献１では、従来の心電計のデータに被験者の加速度のデータを加えて、被験者の呼吸の状態（呼吸数）を推定する方法が開示されている。この非特許文献１に開示された方法では、心電位のＲＲ間隔とＲ波振幅と加速度の３つの指標の重み付け平均をとることによって、推定精度を向上させる。非特許文献１に開示された方法では、ＲＲ間隔から推定した呼吸レートをＢＲ_RSA、Ｒ波振幅から推定した呼吸レートをＢＲ_QRSa、加速度データから推定した呼吸レートをＢＲ_accelとすると、その重み付け平均をとった呼吸レートＢＲ_combinedは、次式のようになる。
　ＢＲ_combined＝（Ｑ_RSAＢＲ_RSA＋Ｑ_QRSaＢＲ_QRSa＋Ｑ_accelＢＲ_QRSa）
　　　　　　　　／（Ｑ_RSA＋Ｑ_QRSa＋Ｑ_accel）　　　　　　　・・・（１）

　式（１）の重み付け定数Ｑ_RSA，Ｑ_QRSa，Ｑ_accelは、ＲＲ間隔、Ｒ波振幅、加速度のそれぞれの指標によって異なり、信号品質によって０～１の範囲で調整される。時系列データのピーク検出により、ピーク振幅の標準偏差、ピーク振幅の平均値、極小値間の時間変化定数、極大値と極小値の数と全ての極大値と極小値の数の比の４つのパラメータの線形和として誤差定数Ｅが算出され、信号品質Ｑが式（２）で計算される。
　Ｑ＝ｅｘｐ｛－（Ｅ／τ）｝　　　　　　　　　　　　　　・・・（２）
ここでτは経験的に定められた定数である。なお、信号振幅がノイズレベルの場合、Ｑはゼロに設定される。

　また、このような複数のセンサを利用したロバストな統合演算処理については、ロボット制御やバイオ医療、装置モニタリング等の分野において、カルマンフィルタを用いた方法が知られている（非特許文献２）。カルマンフィルタに基づく統合処理は、バイオ医療分野では、心電波形と動脈血圧による心拍数のロバスト計測に適用されている（非特許文献３）。

　しかしながら、非特許文献１に開示された方法では、心拍周期によって呼吸サンプリング間隔が制限されるため、周期の早い呼吸曲線を計測することが難しいという問題点があった。また、非特許文献１に開示された方法では、心電位のＲ波振幅の大小や加速度計の信号波形振幅の大小が被験者の体動によるアーチファクトからも発生し、また個人差により、呼吸レートの算出に用いる定数τを変えなければならないため、必ずしもＳＮ（signal-to-noise）比が向上せず、呼吸推定精度が向上しないという問題点があった。

　また、被験者の呼吸状態を推定する分野においては、従来、カルマンフィルタに基づく信号の統合処理は実現されていなかった。

特開２０１５－８３０４５号公報

A.M.Chan，N.Ferdosi，and R.Narasimhan，"Ambulatory Respiratory Rate Detection using ECG and a Triaxial Accelerometer"，35th Annual International Conference of the IEEE EMBS，pp.4058-4061，July，2013 R.C.Luo，C.C.Yih，and K.L.Su，"Multisensor fusion and integration：approaches，applications，and future research directions"，IEEE Journal Sensor，pp.107-119，Vol.2，no.2，Apr.2002 Q Li，R G Mark，and G D Clifford，"Robust heart rate estimation from multiple asynchronous noisy sources using signal quality indices and a Kalman filter"，Physiological Measurement，29，pp.15-32，2008

　本発明は、上記課題を解決するためになされたもので、様々な状況下での呼吸推定精度を高めることができる呼吸推定方法および装置を提供することを目的とする。

　本発明の呼吸推定方法は、被験者の心機能に関する第１のデータの時系列信号と前記被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から前記被験者の呼吸数をそれぞれ推定する第１のステップと、前記第１のデータから推定した呼吸数と前記第２のデータから推定した呼吸数の各々についてカルマンフィルタによりノイズを濾した呼吸数をそれぞれ推定する第２のステップと、この第２のステップで得られた複数の呼吸数の推定値の重み付け平均化処理を実行する第３のステップとを含むことを特徴とするものである。

　また、本発明の呼吸推定装置は、被験者の心機能に関する第１のデータの時系列信号と前記被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から前記被験者の呼吸数をそれぞれ推定する呼吸数推定部と、前記第１のデータから推定した呼吸数と前記第２のデータから推定した呼吸数の各々についてノイズを濾した呼吸数をそれぞれ推定するカルマンフィルタと、このカルマンフィルタで得られた複数の呼吸数の推定値の重み付け平均化処理を実行する統合処理部とを備えることを特徴とするものである。

　また、本発明の呼吸推定方法は、被験者の心機能に関する第１のデータの時系列信号と前記被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から位相情報をそれぞれ抽出する第１のステップと、前記第１のデータの位相情報と前記第２のデータの位相情報の各々についてカルマンフィルタによりノイズを濾した位相情報をそれぞれ推定する第２のステップと、この第２のステップで得られた複数の推定位相値の重み付け平均化処理を実行する第３のステップと、この第３のステップで統合された位相値を周波数に変換することにより、前記被験者の呼吸周波数を求める第４のステップとを含むことを特徴とするものである。

　また、本発明の呼吸推定装置は、被験者の心機能に関する第１のデータの時系列信号と前記被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から位相情報をそれぞれ抽出する特徴量抽出部と、この特徴量抽出部で得られた前記第１のデータの位相情報と前記第２のデータの位相情報の各々についてノイズを濾した位相情報をそれぞれ推定するカルマンフィルタと、このカルマンフィルタで得られた前記第１のデータの推定位相値と前記第２のデータの推定位相値の重み付け平均化処理を実行する統合処理部と、この統合処理部で統合された位相値を周波数に変換することにより、前記被験者の呼吸周波数を求める呼吸周波数変換部とを備えることを特徴とするものである。

　本発明によれば、被験者の心機能に関する第１のデータの時系列信号と被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から被験者の呼吸数を推定し、第１のデータから推定した呼吸数と第２のデータから推定した呼吸数の各々についてノイズを濾した呼吸数をカルマンフィルタで推定し、このカルマンフィルタで得られた複数の呼吸数の推定値を重み付け平均化処理することにより、個人差や幅広い測定状態に対して被験者の呼吸数を推定可能であり、かつ被験者の体動のアーチファクトの影響も低減することができ、呼吸推定精度を高めることができる。

　また、本発明によれば、被験者の心機能に関する第１のデータの時系列信号と被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から位相情報を抽出し、これらの位相情報の各々をカルマンフィルタで処理し、第１のデータの推定位相値と第２のデータの推定位相値とを重み付け平均化処理することにより、個人差や幅広い測定状態に対して被験者の呼吸周波数を推定可能であり、かつ被験者の体動のアーチファクトの影響も低減することができ、呼吸推定精度を高めることができる。

図１は、本発明の第１実施例に係る呼吸推定装置の１構成例を示すブロック図である。 図２は、本発明の第１実施例に係る呼吸推定装置の動作を説明するフローチャートである。 図３は、本発明の第１実施例に係る呼吸推定装置の呼吸数推定部の詳細な動作を説明するフローチャートである。 図４は、本発明の第１実施例に係る呼吸推定装置のカルマンフィルタのブロック線図である。 図５は、本発明の第１実施例に係る呼吸推定装置の重み付け定数生成部と重み付け演算部と適応推定部の動作を説明するフローチャートである。 図６は、本発明の第１実施例に係る呼吸推定装置の統合処理部の構成例を示すブロック図である。 図７は、本発明の第１実施例に係る呼吸推定装置の統合処理部の詳細な動作を説明するフローチャートである。 図８は、本発明の第１実施例に係る呼吸推定装置による呼吸数の統合結果を示す図である。 図９は、本発明の第２実施例に係る呼吸推定装置の構成を示すブロック図である。 図１０は、本発明の第２実施例に係る呼吸推定装置の動作を説明するフローチャートである。 図１１は、本発明の第２実施例において特徴量抽出処理以降の動作を詳細に説明するフローチャートである。 図１２Ａは、本発明の第２実施例に係る呼吸推定装置の特徴量抽出部の構成を示すブロック図、図１２Ｂ、図１２Ｃ、図１２Ｄは、特徴量抽出部の各部の信号波形を示す図である。 図１３は、本発明の第２実施例に係る呼吸推定装置の特徴量抽出部の詳細な動作を説明するフローチャートである。 図１４は、本発明の第２実施例に係る呼吸推定装置のカルマンフィルタの処理を説明する図である。 図１５は、本発明の第２実施例に係る呼吸推定装置のカルマンフィルタのブロック線図である。 図１６は、本発明の第２実施例におけるカルマンフィルタ処理後の推定位相を示す図である。 図１７は、角度変位の位相情報を基に推定した呼吸周波数信号列、および角度変位の位相情報にカルマンフィルタ処理を適用した後の位相情報を基に推定した呼吸周波数信号列を示す図である。 図１８は、本発明の第２実施例に係る呼吸推定装置の統合処理部の構成を示すブロック図である。 図１９は、本発明の第２実施例に係る呼吸推定装置の統合処理部の詳細な動作を説明するフローチャートである。 図２０は、本発明の第２実施例に係る呼吸推定装置で得られた呼吸周波数信号列、およびカルマンフィルタ処理後の角度変位、ＲＲ間隔、Ｒ波振幅の位相情報を基に推定した呼吸周波数信号列を示す図である。 図２１は、本発明の第１、第２実施例に係る呼吸推定装置を実現するコンピュータの構成例を示すブロック図である。 図２２は、ウェアラブルセンサを人体に装着した様子を示す模式図である。

［第１実施例］
　以下、本発明の実施例について図面を参照して説明する。図１は本発明の第１実施例に係る呼吸推定装置の構成を示すブロック図である。呼吸推定装置は、被験者の心電位波形を測定する心電計１と、被験者の呼吸運動による３軸加速度を検出する３軸加速度計２と、心電計１と３軸加速度計２で得られたデータを記憶する記憶部３，４と、心電位波形からＲ波振幅を検出するＲ波振幅検出部５と、心電位波形からＲＲ間隔を検出するＲＲ間隔検出部６と、３軸加速度信号から加速度ベクトルの角度変位を検出する加速度変位検出部７とを備えている。また、呼吸推定装置は、Ｒ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々をサンプリングする標本化部８と、Ｒ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々を帯域制限するバンドパスフィルタ９とを備えている。また、呼吸推定装置は、バンドパスフィルタ９で得られたＲ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々をフーリエ変換して、Ｒ波振幅、ＲＲ間隔、角度変位の各々の信号の周波数スペクトルを求め、これら周波数スペクトルの各々から被験者の呼吸数を推定する呼吸数推定部１０を備えている。また、呼吸推定装置は、Ｒ波振幅の時系列信号から推定した呼吸数とＲＲ間隔の時系列信号から推定した呼吸数と角度変位の時系列信号から推定した呼吸数の各々についてノイズを濾した呼吸数を推定するカルマンフィルタ１１と、カルマンフィルタ１１で得られた呼吸数の推定値を重み付け平均化することにより、データを統合する統合処理部１２とを備えている。さらに、呼吸推定装置は、Ｒ波振幅、ＲＲ間隔、角度変位の各々の信頼度を示すエラー係数を測定するエラー測定部１３と、エラー係数に基づいて重み付け行列の成分である重み付け定数を求める重み付け定数生成部１４と、カルマンフィルタ１１の推定誤差が最小となるようにカルマンフィルタ１１の共分散推定誤差行列を更新する適応推定部１５と、重み付け行列と共分散推定誤差行列に基づいてカルマン利得を更新する重み付け演算部１６とを備えている。

　図２は呼吸推定装置の動作を説明するフローチャートである。心電計１は、被験者の心電位波形を測定し、心電位波形の時系列信号列を出力する（図２ステップＳ１００）。記憶部３は、心電計１から出力された心電位波形の時系列信号列を記憶する。

　周知のとおり、心電位波形は、連続した心拍波形からなり、１つの心拍波形は、心房や心室の活動を反映したＰ波、Ｑ波、Ｒ波、Ｓ波、Ｔ波等の成分からなっている。Ｒ波振幅検出部５は、記憶部３に格納された心電位波形の信号からＲ波の振幅を検出する（図２ステップＳ１０１）。

　具体的には、Ｒ波振幅検出部５は、ウェイブレット変換を用いて、Ｒ波のピーク、Ｑ波のピーク、Ｓ波のピークを検出する。そして、Ｒ波振幅検出部５は、Ｒ波のピーク値からＳ波のピーク値までの振幅をＲ波振幅とする。このように、Ｒ波とＳ波との差を用いることにより、心電位波形のベースラインの変動の影響を低減することができる。Ｒ波振幅検出部５は、心電位波形のＲ波ごとに振幅検出を行う。

　ＲＲ間隔検出部６は、記憶部３に格納された心電位波形の信号からＲＲ間隔を検出する（図２ステップＳ１０２）。ＲＲ間隔検出部６は、上記と同様にウェイブレット変換を用いてＲ波のピークを検出し、Ｒ波と１つ前のＲ波の時間間隔であるＲＲ間隔を求める。ＲＲ間隔検出部６は、心電位波形のＲ波ごとにＲＲ間隔検出を行う。

　一方、３軸加速度計２は、例えば被験者の胸部に装着され、被験者の呼吸運動による３軸加速度を検出して、３軸加速度の時系列信号列を出力する（図２ステップＳ１０３）。記憶部４は、３軸加速度計２から出力された３軸加速度の時系列信号列を記憶する。なお、本実施例の３軸加速度計２は、回転速度検出手段、例えばジャイロセンサー（不図示）を具備しており、Ｘ方向、Ｙ方向、Ｚ方向の３軸の各々について回転速度を検出できるようになっている。

　加速度変位検出部７は、記憶部４に格納された３軸加速度の信号から加速度ベクトルの角度変位を検出する（図２ステップＳ１０４）。この角度変位を検出するには、Ｘ、Ｙ、Ｚの３軸方向の加速度変位の平均から加速度の変化面を規定した後に、Ｘ、Ｙ、Ｚの３軸の加速度データからなる加速度ベクトルの前記変化面への射影ベクトルの角度を角度変位として算出すればよい。

　このような方法は、例えば文献「A.Bates，M.J.Ling，J.Mann and D.K.Arvind，”Respiratory rate and flow waveform estimation from tri-axial accelerometer data”，International Conference on Body Sensor Network，pp.144-150，June 2010」に開示されている。加速度変位検出部７は、このような角度変位の検出を加速度のサンプリング周期ごとに行う。

　続いて、標本化部８は、Ｒ波振幅検出部５から出力されるＲ波振幅の時系列信号、ＲＲ間隔検出部６から出力されるＲＲ間隔の時系列信号、加速度変位検出部７から出力される角度変位の時系列信号の各々を、心電計１のサンプリング周波数および３軸加速度計２のサンプリング周波数よりも遅いサンプリング周波数（例えば１５Ｈｚ間隔）でダウンサンプリングする（図２ステップＳ１０５）。さらに、標本化部８は、Ｒ波振幅の時系列信号、ＲＲ間隔の時系列信号については、１５Ｈｚ間隔へのダウンサンプリングを行った後に、５秒ステップサイズで３０秒間の時間領域にサンプリングする。ここで、ＤＣ（直流）成分は削除する。

　バンドパスフィルタ９は、標本化部８が取得したＲ波振幅の時系列信号、ＲＲ間隔の時系列信号、および角度変位の時系列信号の各々を帯域制限する（図２ステップＳ１０６）。バンドパスフィルタ９を用いる理由は、人の呼吸周波数が低周波のみに限られるためである。このバンドパスフィルタ９の典型的な通過帯域は、例えば０．１５～０．４Ｈｚである。

　なお、Ｒ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々から得られる呼吸数の推定値に相違があるため、これら時系列信号の各々に対するバンドパスフィルタ９の通過帯域を異なる値に設定してもよい。例えば角度変位から得られる呼吸数は６～６０ＢｒＰＭ（０．１～１．０Ｈｚ）であり、Ｒ波振幅から得られる呼吸数は６ＢｒＰＭ～心拍数の半分程度であり、ＲＲ間隔から得られる呼吸数は９ＢｒＰＭ～心拍数の半分程度である。

　次に、呼吸数推定部１０は、バンドパスフィルタ９によって帯域制限されたＲ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々から被験者の呼吸数を推定する（図２ステップＳ１０７）。本実施例では、人の呼吸の動きを波形で表したときに、呼吸曲線の理想モデルにおいては、人の呼吸の動きが正弦波になると仮定し、周波数解析により呼吸数情報を抽出する。図３は呼吸数推定部１０の詳細な動作を説明するフローチャートである。

　まず、呼吸数推定部１０は、Ｒ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々にハミング窓関数を掛け合わせた後に、Ｒ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々を高速フーリエ変換して、Ｒ波振幅、ＲＲ間隔、角度変位の各々の信号の周波数スペクトルを求める（図３ステップＳ１０７０）。周知のとおり、ハミング窓関数は、所望の入力データを切り出すために使用される。

　続いて、呼吸数推定部１０は、高速フーリエ変換で得られたＲ波振幅の周波数スペクトル、ＲＲ間隔の周波数スペクトル、角度変位の周波数スペクトルの各々について、ピークの周波数を被験者の呼吸周波数として検出する。

　このとき、呼吸数推定部１０は、周波数スペクトルに複数のピークが現れた場合には、強度が最大のピークの周波数を当該周波数スペクトルから求まる呼吸周波数とする。呼吸数は呼吸周波数の逆数である。こうして、呼吸数推定部１０は、Ｒ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々から被験者の呼吸数を推定することができる（図３ステップＳ１０７１）。

　次に、カルマンフィルタ１１は、呼吸数推定部１０がＲ波振幅の時系列信号から推定した呼吸数と呼吸数推定部１０がＲＲ間隔の時系列信号から推定した呼吸数と呼吸数推定部１０が角度変位の時系列信号から推定した呼吸数の各々についてノイズを濾した呼吸数を推定する（図２ステップＳ１０８）。

　図４はカルマンフィルタ１１の構成を示すブロック線図である。ベイジアン推定と自動再帰推定に基づくカルマンフィルタは、ノイズが無相関で、ゼロ平均、白色雑音の際に、最適線形推定器であることが知られている。このような条件が真でなくても、最適ではないが最良の線形推定器と考えられる。呼吸数推定部１０がＲ波振幅の時系列信号から推定した呼吸数の値とＲＲ間隔の時系列信号から推定した呼吸数の値と角度変位の時系列信号から推定した呼吸数の値とは、それぞれカルマンフィルタ１１に入力される。生体システムを記述する物理量ｘ_kは再帰的に決定されており，次式で表される。

　ここで、ｕ_kはシステム入力、ｘ_kは雑音無しの物理量である。また、Ａはシステムモデルを示すｎ×ｍ行列（ｍは測定数、ｎは生体システムからの信号数で、本実施例ではｎ＝３）、Ｈは測定系モデルを示すｍ×ｎ行列である。式（３）に示すように、物理量ｘ_k+1には、前時刻における物理量Ａｘ_kと、ｗ_kの生体システム雑音とが含まれる。式（４）に示すように、システム入力ｕ_kには、測定値Ｈｘ_kと、ｖ_kの測定系雑音とが含まれる。本実施例の場合、測定値Ｈｘ_kは、呼吸数推定部１０から入力される呼吸数の値である。システムの更新は下記式により行う。

　ここで、ｄ_kはカルマンイノベーション、ハットｘ_kは物理量ｘ_kの推定値（本実施例では呼吸数の推定値）である。以下、同様に文字上に付した「∧」をハットと呼ぶ。ここで、Ｓ_kは従来のｎ×ｍ行列のカルマン利得とは異なる。変更カルマン利得Ｓ_kは次式より求める。

　ここで、Ｐ_k ^-はハットｘ_kの共分散推定誤差行列であり、Ｈ^Tは行列Ｈの転置行列であり、Ｇ_kは測定系雑音ｖ_kの共分散行列であるＲ_kに基づくパラメータ（行列）である。式（７）によれば、既知の情報を用いて、共分散行列Ｒ_kを調整することによりパラメータＧ_kを求め、変更カルマン利得Ｓ_kを決定することができる。既知の情報は、例えば、外部からのインパルス雑音やセンサ固有の雑音であり、重み付け演算部１６からの定数および定常パラメータからなる。

　パラメータＧ_kには従来のカルマンフィルタで用いることができない情報が含まれているので、共分散推定誤差行列Ｐ_k ^-を更新するのに用いるのは適切ではない。パラメータＧ_kの推定方法は後述するが、本実施例では、共分散推定誤差行列Ｐ_k ^-の更新には、共分散行列Ｒ_kを用いた標準的な算定法を採用する。従来のカルマン利得行列Ｋ_kは次のように定義される。

　ここで、Ｐ_k ⁺はＡハットｘ_kの共分散推定誤差行列であり、Ｑ_kは生体システム雑音ｗ_kの共分散行列、Ｉは単位行列、Ａ^Tは行列Ａの転置行列である。式（８）、式（１０）に用いる共分散行列Ｑ_k、Ｒ_kの値を適応的に求める手法は複数存在する。本実施例では、適応推定部１５において、次式を用い、共分散行列Ｒ_kの値を推定する。

　ここで、ｍは推定に用いる窓長さである。式（１１）の方法では、カルマンイノベーションｄ_kが窓のなかで定常的であることが必要である。この前提は、非定常的な生体システムでは近似となる。また、カルマンフィルタが最適の際に、カルマンイノベーションｄ_kはゼロ平均、白色雑音となる。もしカルマンイノベーションｄ_kが非ゼロ平均で有色であるならば、共分散行列Ｒ_kを調整することで、白色化を試みる例がある。この際に、定常システムでは、最適な共分散行列Ｒ_kは最終的に一定値へ到達する。本実施例では、非定常的な生体システムを前提とするため、ある限られた窓長さｍで、定常状態となることを仮定する。本実施例では、同様の近似をカルマン利得等の他の定常パラメータにも適用する。

　生体システム雑音ｗ_kの共分散行列Ｑ_kに関しては、共分散行列Ｒ_kと同様に推定手法が知られているが、定常状態が仮定できないため、誤差を含むことが考えられる。そこで、本実施例では、共分散行列Ｑ_kを、予め人工的に高い値に設定された成分を有する所定の行列とし、モデルの誤差が大きいものという前提を課すことで、新たな推定値の信頼性が古い推定値よりも高くなるようにする。

　次に、重み付け定数生成部１４について説明する。前述したとおり、既知の情報をカルマンフィルタ１１に注入する必要がある。本実施例では、心電計１および３軸加速度計２から得られた信号のエラーを測定することにより、重み付け行列Ｗ_kの成分である重み付け定数を求める。情報の信頼度に応じて、システムにフィードバックする。

　このようなフィードバックする定数（エラー係数）を－１から＋１の間で定める。情報の信頼度が低い場合、エラー係数を－１とし、情報の信頼度が高い場合、エラー係数を＋１とする。エラーが未知の場合は、エラー係数を０とする。エラー測定部１３によるエラー係数の具体的な決定方法については後述する。

　エラー係数はカルマンフィルタ状態と関連する重み付け行列Ｗ_kの対角成分に配置される。例えば加速度のエラー係数を重み付け行列Ｗ_kの対角成分Ｗ₁₁として配置し、ＲＲ間隔のエラー係数を重み付け行列Ｗ_kの対角成分Ｗ₂₂として配置し、Ｒ波振幅のエラー係数を重み付け行列Ｗ_kの対角成分Ｗ₃₃として配置する。

　本実施例では、信号間の重み付け共分散も必要である。２つの信号は同様に誤差が小さい、誤差が大きい場合には似たようなエラー係数となるが、２つの信号が依存していることを意味していない。しかしながら、それらの信号の誤差の大小は、独立の雑音源が関与しているはずである。

　本実施例では、エラー係数間の差を用いて、重み付け行列Ｗ_kの非対角成分Ｗ_ij，Ｗ_ji（ｉ，ｊは重み付け定数Ｗの次元を示すインデックスであり、ｉ≠ｊ）として配置される重み付け定数を次式により計算する。

　なお、時間インデックスを示すｋは式（１２）に記述していない。こうして、重み付け定数生成部１４は、エラー測定部１３から出力されるエラー係数に基づいて重み付け行列Ｗ_kを決定する。
　次に、式（７）に用いるパラメータＧ_kを求めるためには、単純に重み付け行列Ｗ_kを用いて、式（１３）の計算を行う。

　式（１３）における記号「○」はアダマール積を意味する。αはエラー係数の強さを調整するスケールファクタである（０≦α≦１）。αは収束性を高めるために、例えば、０．５と設定しておく。

　図５は重み付け定数生成部１４と適応推定部１５と重み付け演算部１６の動作を説明するフローチャートである。上記のように、ある窓長さｍ（測定数）におけるエラー測定で得られたエラー係数に基づいて、重み付け定数生成部１４は、短期的なシステム推定を行い、重み付け行列Ｗ_kを決定する（図５ステップＳ２００）。

　適応推定部１５は、カルマンフィルタ状態と過去の状態から、理想的な雑音状態（カルマンイノベーションが定常状態）を仮定して、周知の手法により長期的なシステム推定を行い、共分散行列Ｒ_kを推定し（図５ステップＳ２０１）、さらに式（８）～式（１０）により定常パラメータである共分散推定誤差行列Ｐ_k ^-を決定する（図５ステップＳ２０２）。式（８）～式（１０）はハットｘ_kの推定誤差が最小となるようにカルマン利得Ｋ_kを再帰的に決定する周知の手法を表しており、この手法により共分散推定誤差行列Ｐ_k ^-，Ｐ_k ⁺を更新することができる。

　そして、重み付け演算部１６は、重み付け定数生成部１４が決定した重み付け行列Ｗ_kを用いて式（１３）によりパラメータＧ_kを算出し、このパラメータＧ_kと適応推定部１５が決定した共分散推定誤差行列Ｐ_k ^-とを用いて式（７）により変更カルマン利得Ｓ_kを算出する（図５ステップＳ２０３）。以上により、カルマン利得Ｓ_kを適応的に変更することができる。

　本実施例ではシステムを下記のとおり定める。

　式（１４）において、ｘ_accは呼吸数推定部１０が角度変位（３軸加速度）の時系列信号から推定した呼吸数、ｘ_RRIは呼吸数推定部１０がＲＲ間隔の時系列信号から推定した呼吸数、ｘ_RAは呼吸数推定部１０がＲ波振幅の時系列信号から推定した呼吸数である。ｕ_acc，ｕ_RRI，ｕ_RAはそれぞれｘ_acc，ｘ_RRI,ｘ_RAに対応するシステム入力である。

　本実施例は、Ｒ波振幅の時系列信号から推定した呼吸数とＲＲ間隔の時系列信号から推定した呼吸数と角度変位（３軸加速度）の時系列信号から推定した呼吸数とを統合することを特徴としている。すなわち、パラメータＧ_k、変更カルマン利得Ｓ_kおよび式（１６）の行列Ａに基づいて、３つの呼吸数の推定値は式（６）により１つの推定値に結合される。

　なお、測定系モデルは下記で表すことができる。

　本実施例では、上記のとおり、共分散行列Ｒ_kについては適応推定部１５で推定し、一方、共分散行列Ｑ_kについては推定せずに固定値とするため、共分散行列Ｑ_kをあらかじめ定めておく必要がある。システムを単純化するために、プロセス雑音は独立と仮定する。すなわち、共分散行列Ｑ_kの対角成分のみが値を有し、非対角成分が０となるようにする。

　例えば加速度に対応する成分を共分散行列Ｑ_kの対角成分Ｑ₁₁として配置し、ＲＲ間隔に対応する成分を共分散行列Ｑ_kの対角成分Ｑ₂₂として配置し、Ｒ波振幅に対応する成分を共分散行列Ｑ_kの対角成分Ｑ₃₃として配置する。これら対角成分の値は、それぞれの信号（加速度、Ｒ波振幅、ＲＲ間隔）で同じ値とする。

　式（１６）の行列Ａに基づいて、３つの呼吸数の推定値は式（６）により１つの推定値に結合される。ここで、行列Ａを式（１８）のようにした場合、統合処理部１２により呼吸数の推定値を統合する。

　統合処理部１２は、カルマンフィルタ１１から出力された、Ｒ波振幅の時系列信号に基づく呼吸数の推定値とＲＲ間隔の時系列信号に基づく呼吸数の推定値と角度変位（３軸加速度）の時系列信号に基づく呼吸数の推定値とを、カルマンフィルタ１１の自乗推定誤差に基づく重みを用いて重み付け平均化処理することにより、これら呼吸数の推定値を統合する（図２ステップＳ１０９）。

　図６は統合処理部１２の構成を示すブロック図、図７は統合処理部１２の詳細な動作を説明するフローチャートである。統合処理部１２は、Ｒ波振幅の時系列信号に基づく呼吸数、ＲＲ間隔の時系列信号に基づく呼吸数、角度変位（３軸加速度）の時系列信号に基づく呼吸数の各々について、カルマンフィルタ１１の自乗推定誤差から統合処理のための重み付け定数を算出する重み付け定数生成処理部１２０と、この重み付け定数を用いて、複数の呼吸数の推定値を重み付け平均化する重み付け平均化処理部１２１とから構成される。上記のカルマンフィルタ処理では、時刻ごとに呼吸数の推定値が得られると共に変更カルマン利得Ｓ_kおよび共分散推定誤差行列Ｐ_k ^-，Ｐ_k ⁺が更新される。

　カルマンフィルタ１１の自乗推定誤差は、共分散推定誤差行列Ｐ_k ⁺から得ることができる。本実施例では、角度変位（３軸加速度）の時系列信号に基づく呼吸数の自乗推定誤差をσ１、ＲＲ間隔の時系列信号に基づく呼吸数の自乗推定誤差をσ２、Ｒ波振幅の時系列信号に基づく呼吸数の自乗推定誤差をσ３とする。

　重み付け定数生成処理部１２０は、Ｒ波振幅、ＲＲ間隔、角度変位のそれぞれの自乗推定誤差σｉから重み付け定数αｉを式（１９）により算出する（図７ステップＳ１０９０）。

　重み付け定数生成処理部１２０は、式（１９）の算出をＲ波振幅、ＲＲ間隔、角度変位の各々について行う。そして、重み付け平均化処理部１２１は、重み付け定数生成処理部１２０が算出した重み付け定数αｉ（ｉ＝１，２，３）を用いて、角度変位（３軸加速度）の時系列信号に基づく呼吸数の推定値ハットｘ１とＲＲ間隔の時系列信号に基づく呼吸数の推定値ハットｘ２とＲ波振幅の時系列信号に基づく呼吸数の推定値ハットｘ３とを重み付け平均化処理した呼吸数の値である統合出力値ハットｘｆ_kを式（２０）のように算出する（図７ステップＳ１０９１）。統合処理部１２は、以上のような統合処理を時刻ごとに行う。

　次に、エラー測定部１３によるエラー係数の決定方法について説明する。加速度のエラー係数に関しては、３軸加速度計２に具備されたジャイロセンサーにより、Ｘ、Ｙ、Ｚの３軸の各々について回転速度を検出する。エラー測定部１３は、Ｘ、Ｙ、Ｚの３軸の回転速度のうち、全ての回転速度が所定の第１の回転速度閾値（例えば０．２ｒａｄ／ｓｅｃ）以上であれば、加速度のエラー係数を＋１とし、全ての回転速度が所定の第２の回転速度閾値（例えば０．００３ｒａｄ／ｓｅｃ）以下であれば、加速度のエラー係数を－１とする（図５ステップＳ１９０）。

　また、エラー測定部１３は、Ｘ、Ｙ、Ｚの３軸の回転速度の代表値（例えば平均値）が第１の回転速度閾値と第２の回転速度閾値の間にある場合には、この代表値に応じて－１～＋１の間の値を割り当てる線形スケーリングの手法により、加速度のエラー係数を決定する（ステップＳ１９０）。

　また、エラー測定部１３は、Ｒ波振幅のエラー係数に関しては、各ビートに合わせて、エラーを抽出する（ステップＳ１９０）。具体的には、エラー測定部１３は、心電計１が測定した心電位波形からＱ波、Ｒ波、Ｓ波を除き、これらを除いた波形に関して、ビート間の相関係数を算出する。エラー測定部１３は、相関係数が所定の相関係数閾値以上の高い値であれば、Ｒ波振幅のエラー係数を－１とし、相間係数が０以下であれば、Ｒ波振幅のエラー係数を＋１とする。相関係数が高いことは雑音源からの雑音が大きく、誤差が大きいことを意味している。

　また、エラー測定部１３は、算出した相関係数が所定の相関係数閾値と０との間にある場合には、相関係数に応じて－１～＋１の間の値を割り当てる線形スケーリングの手法により、Ｒ波振幅のエラー係数を決定する。

　一方、ＲＲ間隔についてはエラー係数を算定する適当な方法がなく、適切なエラー係数を求めることが困難である。このため、エラー測定部１３は、ＲＲ間隔のエラー係数を固定値、具体的には０とする（ステップＳ１９０）。

　本実施例の呼吸推定装置によるデータ再構成結果（呼吸数）のＲＭＳ誤差（ＢｒＰＭ）を表１に示す。

　表１は、加速度（角度変位）から推定した呼吸数の誤差、ＲＲ間隔から推定した呼吸数の誤差、Ｒ波振幅から推定した呼吸数の誤差、これらを統合した呼吸数（統合処理部１２の統合出力）の誤差を、被験者の異なる姿勢や動作において比較した表である。ここでは、臥位、座位、歩行、走行のそれぞれの場合について推定した呼吸数の誤差を示している。

　表１によれば、加速度、ＲＲ間隔、Ｒ波振幅から個別に推定した呼吸数の推定値よりも、統合出力の方がより精度が良いことがわかる。また、歩行や走行のように大きいアーチファクトが存在する動作においても統合出力の誤差が小さく、推定精度が向上していることが分かる。

　統合処理部１２の統合出力において、異なる個人差における推定誤差率を、臥位（Laying down）、座位（Sitting）、歩行（Walking）、走行（Running）の各場合についてパーセンタイルグラフで示したのが、図８である。図８によれば、歩行している複数の被験者について呼吸数を推定した場合でも、３０％以下の誤差率が実現できていることが分かる。

　本実施例は、被験者の年齢、自律神経、皮膚状態、体格の個人差に対応するため、適応的にカルマン利得を変更しながら、複数のセンサデータを統合し，最良のデータを抽出することで被験者の呼吸数を推定するものである。

　３軸加速度のデータを呼吸数の推定に使用する場合、被験者の体動、体格の個人差等により算定式の最適化が困難であった。ＲＲ間隔のデータを呼吸数の推定に使用する場合、心拍周期による呼吸周期算出の限界があり、また被験者の精神状態や加齢により変化する自律神経系の影響により算定式の最適化が困難であった。Ｒ波振幅のデータを呼吸数の推定に使用する場合、被験者の体動や皮膚状態による接触インピーダンスの変化、皮膚状態の個人差等により算定式の最適化が困難であった。これに対して、本実施例では、複数のセンサデータを適応的に統合することで、信号のＳＮ劣化が生じた場合でも、呼吸数の雑音を低減することができ、安定した呼吸数の推定が可能である。また、本実施例では、エラー測定部１３と重み付け定数生成部１４と適応推定部１５と重み付け演算部１６とを設けることにより、カルマンフィルタのカルマン利得に、心機能に関するデータおよび加速度データの信頼度を反映させることができ、呼吸数の推定精度を高めることができる。

　本実施例で説明した記憶部３，４とＲ波振幅検出部５とＲＲ間隔検出部６と加速度変位検出部７と標本化部８とバンドパスフィルタ９とカルマンフィルタ１１と統合処理部１２とエラー測定部１３と重み付け定数生成部１４と適応推定部１５と重み付け演算部１６とは、ＣＰＵ（Central Processing Unit）、記憶装置及びインタフェースを備えたコンピュータと、これらのハードウェア資源を制御するプログラムによって実現することができる。ＣＰＵは、記憶装置に格納されたプログラムに従って本実施例で説明した処理を実行する。

　なお、心電計１は、シャツ等の衣類に取り付けられた電極と、同じく衣類に取り付けられたモニタ装置内の心電位波形信号処理部とから構成され、電極と心電位波形信号処理部との間が配線によって接続されている。同様に、３軸加速度計２は、衣類に取り付けられたセンサ部と、上記のモニタ装置内に設けられた加速度信号処理部とから構成され、センサ部と加速度信号処理部との間が配線によって接続されている。

　本実施例は、このような心電計１と３軸加速度計２とが衣類に取り付けられたウェアラブル機器と一緒に設けてもよいし、別に設けてもよい。すなわち、前記のモニタ装置内に、記憶部３，４とＲ波振幅検出部５とＲＲ間隔検出部６と加速度変位検出部７と標本化部８とバンドパスフィルタ９とカルマンフィルタ１１と統合処理部１２とエラー測定部１３と重み付け定数生成部１４と適応推定部１５と重み付け演算部１６とを設けてもよいし、これらを別の装置内に設けてもよい。

　記憶部３，４とＲ波振幅検出部５とＲＲ間隔検出部６と加速度変位検出部７と標本化部８とバンドパスフィルタ９とカルマンフィルタ１１と統合処理部１２とエラー測定部１３と重み付け定数生成部１４と適応推定部１５と重み付け演算部１６とを、別の装置内に設ける場合には、この装置に対して、心電計１で得られた心電位波形の信号と３軸加速度計２で得られた３軸加速度信号および回転速度信号の情報を例えば無線送信することになる。

［第２実施例］
　次に、本発明の第２実施例について説明する。図９は本発明の第２実施例に係る呼吸推定装置の構成を示すブロック図である。呼吸推定装置は、心電計２１と、３軸加速度計２２と、記憶部２３，２４と、Ｒ波振幅検出部２５と、ＲＲ間隔検出部２６と、加速度変位検出部２７とを備えている。また、呼吸推定装置は、標本化部２８と、バンドパスフィルタ２９と、バンドパスフィルタ２９で得られたＲ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々から位相情報を抽出する特徴量抽出部３０とを備えている。さらに、呼吸推定装置は、Ｒ波振幅の位相情報、ＲＲ間隔の位相情報、角度変位の位相情報の各々についてノイズを濾した位相情報を推定するカルマンフィルタ３１と、カルマンフィルタ３１で得られたＲ波振幅の推定位相値、ＲＲ間隔の推定位相値、角度変位の推定位相値を重み付け平均化することによりデータを統合する統合処理部３２と、統合処理部３２で統合された位相値を周波数に変換した結果を、被験者の呼吸周波数とする呼吸周波数変換部３３とを備えている。

　図１０は呼吸推定装置の動作を説明するフローチャートである。心電計２１は、被験者の心電位波形を測定し、心電位波形の時系列信号列を出力する（図１０ステップＳ３００）。記憶部２３は、心電計２１から出力された心電位波形の時系列信号列を記憶する。

　Ｒ波振幅検出部２５は、記憶部２３に格納された心電位波形の信号からＲ波の振幅を検出する（図１０ステップＳ３０１）。心電位波形のＲ波を認識する方法としては、例えば特開２００３－５６１号公報に開示された技術がある。この技術は、心電位波形の山と谷との振幅に基づいた閾値でＲ波を認識するものである。Ｒ波振幅検出部２５は、心電位波形のＲ波ごとに振幅検出を行う。

　ＲＲ間隔検出部２６は、記憶部２３に格納された心電位波形の信号からＲＲ間隔を検出する（図１０ステップＳ３０２）。心電位波形のＲＲ間隔を検出する技術としては、例えば文献「“ECG Implementation on the TMS320C5515 DSP Medical Development Kit (MDK) with the ADS1298 ECG-FE”，Texas Instruments Incorporated，＜http://www.ti.com/lit/an/sprabj1/sprabj1.pdf＞，2011」に開示された技術がある。この技術は、心電位波形を時間差分した値の変化を基にＲＲ間隔を求めるものである。ＲＲ間隔検出部２６は、心電位波形のＲ波ごとにＲＲ間隔検出を行う。

　一方、３軸加速度計２２は、被験者の胸部に装着され、被験者の呼吸運動による３軸加速度を検出して、３軸加速度の時系列信号列を出力する（図１０ステップＳ３０３）。記憶部２４は、３軸加速度計２２から出力された３軸加速度の時系列信号列を記憶する。

　加速度変位検出部７と同様に、加速度変位検出部２７は、記憶部２４に格納された３軸加速度の信号から加速度ベクトルの角度変位を検出する（図１０ステップＳ３０４）。加速度変位検出部２７は、角度変位の検出を加速度のサンプリング周期ごとに行う。

　続いて、標本化部２８は、Ｒ波振幅検出部２５から出力されるＲ波振幅の時系列信号、ＲＲ間隔検出部２６から出力されるＲＲ間隔の時系列信号、加速度変位検出部２７から出力される角度変位の時系列信号の各々を、心電計２１のサンプリング周波数および３軸加速度計２２のサンプリング周波数よりも遅いサンプリング周波数（例えば１Ｈｚ間隔）でサンプリングする（図１０ステップＳ３０５）。

　バンドパスフィルタ２９は、標本化部２８が取得したＲ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々を帯域制限する（図１０ステップＳ３０６）。バンドパスフィルタ２９を用いる理由は、人の呼吸周波数が低周波のみに限られるためである。このバンドパスフィルタ２９の通過帯域は、例えば０．１５～０．４Ｈｚである。

　特徴量抽出部３０に係る呼吸推定装置の処理の流れを、図１１のフローチャートを用いてより詳細に説明する。特徴量抽出部３０は、バンドパスフィルタ２９で得られたＲ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々から瞬時的な線形化された位相情報及び瞬時的な振幅情報及び瞬時的な周波数情報を抽出する（図１０、図１１ステップＳ３０７）。本実施例では、人の呼吸の動きを波形で表したときに、呼吸曲線の理想モデルにおいては、短時間での呼吸数の変化が無視できるほど小さいことと呼吸の動きが正弦波になると仮定し、ヒルベルト変換により瞬時的な位相情報を抽出する。

　図１２Ａは特徴量抽出部３０の構成を示すブロック図、図１３は特徴量抽出部３０の詳細な動作を説明するフローチャートである。特徴量抽出部３０は、ヒルベルト変換部３１０と、角度算出部３１１と、アンラップ処理部３１２とから構成される。

　まず、ヒルベルト変換部３１０は、バンドパスフィルタ２９によって帯域制限されたＲ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号をヒルベルト変換して、位相がπ／２異なる２つの信号成分（実部成分および虚部成分）を生成する（図１３ステップＳ３０７０）。ヒルベルト変換部３１０に入力される信号をＡｅｘｐ（－ｉθ）の正弦波で表現すると（図１２Ｂ）、生成される実部成分はＡｃｏｓθで表現される正弦波であり、虚部成分はｉＡｓｉｎθで表現される正弦波の信号である。Ａは振幅、θは角度、ｉは虚数単位である。

　続いて、角度算出部３１１は、ヒルベルト変換部３１０によって生成された実部成分Ａｃｏｓθと虚部成分ｉＡｓｉｎθとから角度θ（－πから＋π）を算出する（図１３ステップＳ３０７１）。

　最後に、アンラップ処理部３１２は、角度算出部３１１が算出した角度θを連続した位相値に線形化する位相アンラップを行う（図１３ステップＳ３０７２）。角度算出部３１１が算出する角度θは－πから＋πの間の値となるので、図１２Ｃに示すように隣り合う点に２πの位相飛びが生じる場合がある。そこで、アンラップ処理部３１２は、例えば２πを足したり引いたりすることで、位相を繋ぎ合わせるようにする。これにより、図１２Ｄに示すような連続した位相が得られる。

　特徴量抽出部３０は、以上のような処理を、バンドパスフィルタ２９によって帯域制限されたＲ波振幅の時系列信号、ＲＲ間隔の時系列信号、角度変位の時系列信号の各々について行う。すなわち、ヒルベルト変換部３１０と角度算出部３１１とアンラップ処理部３１２の組を３組設け、Ｒ波振幅、ＲＲ間隔、角度変位の３つの時系列信号の処理を同時並行に行うようにすればよい。

　次に、カルマンフィルタ３１は、Ｒ波振幅の位相情報、ＲＲ間隔の位相情報、角度変位の位相情報の各々についてノイズを濾した位相情報を推定する（図１０、図１１ステップＳ３０８）。図１４に示すように、カルマンフィルタ３１は、Ｒ波振幅、ＲＲ間隔、角度変位の各々について設けられ、特徴量抽出部３０で得られたＲ波振幅の位相情報、ＲＲ間隔の位相情報、角度変位の位相情報がそれぞれ対応するカルマンフィルタ３１－１，３１－２，３１－３に入力される。

　図１５はカルマンフィルタ３１のブロック線図である。なお、ここでは３つのカルマンフィルタ３１－１，３１－２，３１－３のうちの１つのカルマンフィルタについてブロック線図を示す。カルマンフィルタは、ベイジアン推定と自動再帰推定に基づくものである。システムへの入力はガウシアン雑音を有するという前提の下、カルマン利得Ｋは最小自乗近似により求める。生体システムを記述する物理量ｘ（ｋ）は再帰的に決定されており、次式で表される。

　ここで、ｕ（ｋ）はシステム入力、ｘ（ｋ）は雑音無しの物理量である。また、ｍは測定数、ｎは生体システムからの信号数、Ａはシステムモデルを示すｎ×ｍ行列、Ｈは測定系モデルを示すｍ×ｎ行列である。式（２１）に示すように、物理量ｘ（ｋ＋１）には、前時刻における物理量Ａｘ（ｋ）と、ｗ_s（ｋ）の生体システム雑音とが含まれる。式（２２）に示すように、システム入力ｕ（ｋ）には、測定値Ｈｘ（ｋ）と、ｗ_m（ｋ）の測定系雑音とが含まれる。本実施例の場合、測定値Ｈｘ（ｋ）は特徴量抽出部３０から入力される位相値である。

　ここで、Ｋ（ｋ）はカルマン利得を示すｎ×ｍ行列、ハットｘ（ｋ）は物理量ｘ（ｋ）の推定値（本実施例では位相の推定値）である。カルマン利得Ｋ（ｋ）は次の式（２４）～式（２６）により求めることができる。

　ここで、Ｒはセンサ雑音に関する共分散マトリクス、Ｑは生体システム雑音に関する共分散マトリクス、Ｐは推定値に関する共分散マトリクスである。Ｈ^T，Ａ^Tはそれぞれ行列Ｈ，Ａの転置行列である。カルマン利得Ｋ（ｋ）は、測定系雑音ｗ_m（ｋ）を最小化するように再帰的に決定される。共分散マトリクスＰ（ｋ｜ｋ）の対角成分はフィルタ処理により最小化される自乗推定誤差である。

　センサ雑音に関する共分散マトリクスＲを予め正確に設定するために、被験者の呼吸レートが安定している時刻域を利用し、特徴量抽出部３０で得られたＲ波振幅の位相情報、ＲＲ間隔の位相情報、角度変位の位相情報の各々について、当該時刻域における標準偏差σを計算する。

　Ｒ波振幅、ＲＲ間隔、角度変位のそれぞれの入力が互いに独立であるという前提では、共分散マトリクスＲは対角行列である。そこで、Ｒ波振幅の位相情報について得られた標準偏差σの平均値を、Ｒ波振幅用のカルマンフィルタ３１－１の共分散マトリクスＲの対角成分として予め設定し、ＲＲ間隔の位相情報について得られた標準偏差σの平均値を、ＲＲ間隔用のカルマンフィルタ３１－２の共分散マトリクスＲの対角成分として予め設定し、角度変位の位相情報について得られた標準偏差σの平均値を、角度変位用のカルマンフィルタ３１－３の共分散マトリクスＲの対角成分として予め設定する。センサ雑音に関する共分散マトリクスＲは、測定環境に依存するものの、個人差による違いはさほど大きくないと考えられる。

　一方、生体システム雑音に関する共分散マトリクスＱは、生体システム雑音であり、個人差が反映される。この共分散マトリクスＱは対角行列である。共分散マトリクスＱの対角成分を変更可能なパラメータとして、被験者毎に共分散マトリクスＱの対角成分の最適値を数値的にテストし、被験者に応じた対角成分の値を予め設定しておくようにすればよい。以上のように、予め被験者に呼吸をさせつつデータを計測し、計測したデータに基づいて共分散マトリクスＲ，Ｑを決定しておくことで、位相の推定精度を高めることができる。ここで、各パラメータの例示を行う。ｘ（ｋ）は位相θｋとドットθｋのベクトルとして表現できる。ドットθｋは位相θｋの微分である。ここでは、文字上に付した「・」をドットと呼ぶ。共分散マトリクスＱ、Ｒの定数ｑ１およびｒ１は設定可能なパラメータであり、上述の通りに設定を行う。例えば、ｑ１＝１×１０^-3、ｒ１＝２．４である。

　カルマンフィルタ３１（３１－１～３１－３）は、フィルタ処理を標本化部２８のサンプリング周期ごと（本実施例の例では１秒ごと）に行う。
　図１６に本実施例におけるカルマンフィルタ処理後の推定位相を示す。７０はカルマンフィルタ３１に入力される測定値（位相値）を示し、７１はカルマンフィルタ３１から出力される推定値（位相値）を示し、７２は位相の理想直線を示している。測定値に対して推定値は理想直線との誤差が小さくなっており，測定値では標準誤差０．５６６３（ｒａｄ）に対して、フィルタ処理後の推定値では標準誤差が０．２７１１（ｒａｄ）となっている。

　図１７の（ａ）は角度変位の位相情報を基に推定した呼吸周波数信号列、図１７の（ｂ）はカルマンフィルタ処理後の角度変位の位相情報を基に推定した呼吸周波数信号列を示す図である。図１７の（ａ）と図１７の（ｂ）を比較すると、カルマンフィルタ処理を適用することにより、被験者の体動等による呼吸レートの揺らぎや雑音が低減され、呼吸周波数の推定値が安定していることが分かる。

　次に、統合処理部３２は、カルマンフィルタ３１（３１－１～３１－３）で得られたＲ波振幅の推定位相値、ＲＲ間隔の推定位相値、角度変位の推定位相値を、カルマンフィルタ３１（３１－１～３１－３）の自乗推定誤差に基づく重みを用いて重み付け平均化処理することにより、Ｒ波振幅、ＲＲ間隔、角度変位の推定位相値のデータを統合する（図１０、図１１ステップＳ３０９）。

　図１８は統合処理部３２の構成を示すブロック図、図１９は統合処理部３２の詳細な動作を説明するフローチャートである。統合処理部３２は、重み付け定数生成処理部３２０と、重み付け平均化処理部３２１とから構成される。上記のカルマンフィルタ処理では、時刻ごと（本実施例の例では１秒ごと）に推定位相値が得られると共に自乗推定誤差が更新される。

　本実施例では、Ｒ波振幅用のカルマンフィルタ３１－１におけるカルマンフィルタ処理の自乗推定誤差をσ１、ＲＲ間隔用のカルマンフィルタ３１－２におけるカルマンフィルタ処理の自乗推定誤差をσ２、角度変位用のカルマンフィルタ３１－３におけるカルマンフィルタ処理の自乗推定誤差をσ３とする。Ｒ波振幅、ＲＲ間隔、角度変位のそれぞれのカルマンフィルタ処理の自乗推定誤差σｉ（ｉ＝１，２，３）は、次式のように表される。

　重み付け定数生成処理部３２０は、Ｒ波振幅、ＲＲ間隔、角度変位のそれぞれの自乗推定誤差σｉから重み付け定数αｉを式（２８）により算出する（図１９ステップＳ３０９０）。

　重み付け定数生成処理部３２０は、式（２８）の算出をＲ波振幅、ＲＲ間隔、角度変位の各々について行う。そして、重み付け平均化処理部３２１は、重み付け定数生成処理部３２０が算出した重み付け定数αｉ（ｉ＝１，２，３）を用いて、Ｒ波振幅の推定位相値ハットｘ１、ＲＲ間隔の推定位相値ハットｘ２、角度変位の推定位相値ハットｘ３を重み付け平均化処理した位相値である統合出力値ハットｘｆ（ｋ）を式（２９）のように算出する（図１９ステップＳ３０９１）。

　統合処理部３２は、以上のような統合処理を標本化部２８のサンプリング周期ごと（本実施例の例では１秒ごと）に行う。

　次に、呼吸周波数変換部３３は、統合処理部３２によって統合された位相値を周波数に変換して呼吸周波数信号を出力する（図１０、図１１ステップＳ３１０）。統合処理部３２から出力される位相値を時間微分すれば、瞬時角周波数が得られるので、この瞬時角周波数を２πで割ることにより、呼吸周波数を求めることができる。こうして、呼吸周波数変換部３３は、呼吸周波数のデータを時系列順に並べた呼吸周波数信号列を生成して出力する。

　図２０の（ａ）は本実施例の呼吸推定装置で得られた呼吸周波数信号列、図２０の（ｂ）はカルマンフィルタ処理後の角度変位の位相情報を基に推定した呼吸周波数信号列、図２０の（ｃ）はカルマンフィルタ処理後のＲＲ間隔の位相情報を基に推定した呼吸周波数信号列、図２０の（ｄ）はカルマンフィルタ処理後のＲ波振幅の位相情報を基に推定した呼吸周波数信号列を示す図である。

　本実施例は、被験者の年齢、自律神経、皮膚状態、体格の個人差に対応するため、複数のセンサデータを統合し、最良のデータを抽出することで被験者の呼吸周波数を推定するものである。

　３軸加速度のデータを呼吸周波数の推定に使用する場合、体動、体格の個人差等による測定誤差を除くことが困難であった。ＲＲ間隔のデータを呼吸周波数の推定に使用する場合、心拍周期による呼吸周期算出の限界があり、また精神状態や加齢により変化する自律神経系の影響による測定誤差を除くことが困難であった。Ｒ波振幅のデータを呼吸周波数の推定に使用する場合、体動や皮膚状態による接触インピーダンスの変化、皮膚状態の個人差等による測定誤差を除くことが困難であった。これに対して、本実施例では、複数のセンサデータを統合することで、ＳＮ比劣化が生じた場合でも、呼吸周波数のノイズを低減することができ、安定した呼吸周波数の推定が可能である。

　本実施例で説明した記憶部２３，２４とＲ波振幅検出部２５とＲＲ間隔検出部２６と加速度変位検出部２７と標本化部２８とバンドパスフィルタ２９と特徴量抽出部３０とカルマンフィルタ３１と統合処理部３２と呼吸周波数変換部３３とは、ＣＰＵ、記憶装置及びインタフェースを備えたコンピュータと、これらのハードウェア資源を制御するプログラムによって実現することができる。ＣＰＵは、記憶装置に格納されたプログラムに従って本実施例で説明した処理を実行する。

　心電計２１は、シャツ等の衣類に取り付けられた電極と、同じく衣類に取り付けられたモニタ装置内の心電位波形信号処理部とから構成され、電極と心電位波形信号処理部との間が配線によって接続されている。同様に、３軸加速度計２２は、衣類に取り付けられたセンサ部と、上記のモニタ装置内に設けられた加速度信号処理部とから構成され、センサ部と加速度信号処理部との間が配線によって接続されている。

　本実施例は、このような心電計２１と３軸加速度計２２とが衣類に取り付けられたウェアラブル機器と一緒に設けてもよいし、別に設けてもよい。すなわち、前記のモニタ装置内に、記憶部２３，２４とＲ波振幅検出部２５とＲＲ間隔検出部２６と加速度変位検出部２７と標本化部２８とバンドパスフィルタ２９と特徴量抽出部３０とカルマンフィルタ３１と統合処理部３２と呼吸周波数変換部３３とを設けてもよいし、これらを別の装置内に設けてもよい。

　記憶部２３，２４とＲ波振幅検出部２５とＲＲ間隔検出部２６と加速度変位検出部２７と標本化部２８とバンドパスフィルタ２９と特徴量抽出部３０とカルマンフィルタ３１と統合処理部３２と呼吸周波数変換部３３とを、別の装置内に設ける場合には、この装置に対して、心電計２１で得られた心電位波形の信号と３軸加速度計２２で得られた３軸加速度信号とを例えば無線送信することになる。

　ここで、第１実施例、第２実施例で説明したコンピュータの構成例を図２１に示す。コンピュータは、ＣＰＵ４０と、記憶装置４１と、インターフェイス装置（以下、Ｉ／Ｆと略する）４２とを備えている。Ｉ／Ｆ４２には、心電計１，２１と３軸加速度計２，２２などが接続される。このようなコンピュータにおいて、本発明の呼吸推定方法を実現させるためのプログラムは、フレキシブルディスク、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ、メモリカードなどの記録媒体に記録された状態で提供され、記憶装置４１に格納される。ＣＰＵ４０は、記憶装置４１に記憶されたプログラムに従って第１実施例または第２実施例で説明した処理を実行する。

　本発明は、３軸加速度計、心電計の組み合わせのみに関わらず、磁場センサやマイクロ波センサなど、その他のセンサから得られるデータの組み合わせにも適用が可能である。

　また、本発明では、被験者の心機能に関するデータの例として、Ｒ波振幅とＲＲ間隔を例に挙げて説明したが、Ｒ波振幅とＲＲ間隔を用いることは本発明において必須の構成要件ではない。心機能に関するデータの別の例としては、Ｓ波振幅の時系列信号、あるいはＲ波のピーク値からＳ波のピーク値までのＲＳ振幅の時系列信号がある。なお、第２実施例ではＲ波そのものの振幅を検出しているが、第１実施例ではＲＳ振幅を検出して、この検出結果をＲ波振幅としている。したがって、第１実施例のＲ波振幅検出部５は、ＲＳ振幅検出部と見なすことも可能である。第１実施例においてＲ波振幅検出部５の代わりに、第２実施例のＲ波振幅検出部２５を用いてもよい。また、心電位波形の信号からＳ波振幅を検出する場合には、Ｓ波振幅検出部を設けるようにすればよい。Ｓ波振幅検出部は、Ｒ波の場合と同様に、閾値を用いることでＳ波を検出することが可能である。Ｓ波振幅検出部は、心電位波形のＳ波ごとに振幅検出を行う。

　また、本発明では、Ｒ波振幅、ＲＲ間隔、角度変位の各々から求めた３つの変数（第１実施例では呼吸数、第２実施例では位相値）をカルマンフィルタ１１，３１に入力しているが、これに限るものではない。被験者の心機能に関するデータから求めた変数と加速度データ（角度変位）から求めた変数のうち、２つの変数または４つ以上の変数をカルマンフィルタ１１，３１に入力するように拡張することも可能である。

　本発明は、人の呼吸を継続的に観察する呼吸連続モニタリングに適用することができる。

　１，２１…心電計、２，２２…３軸加速度計、３，４，２３，２４…記憶部、５，２５…Ｒ波振幅検出部、６，２６…ＲＲ間隔検出部、７，２７…加速度変位検出部、８，２８…標本化部、９，２９…バンドパスフィルタ、１０…呼吸数推定部、１１，３１…カルマンフィルタ、１２，３２…統合処理部、１３…エラー測定部、１４…重み付け定数生成部、１５…適応推定部、１６…重み付け演算部、３０…特徴量抽出部、３３…呼吸周波数変換部、１２０…重み付け定数生成処理部、１２１…重み付け平均化処理部、３１０…ヒルベルト変換部、３１１…角度算出部、３１２…アンラップ処理部、３２０…重み付け定数生成処理部、３２１…重み付け平均化処理部。

Claims (13)

1. 　被験者の心機能に関する第１のデータの時系列信号と前記被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から前記被験者の呼吸数をそれぞれ推定する第１のステップと、
   　前記第１のデータから推定した呼吸数と前記第２のデータから推定した呼吸数の各々についてカルマンフィルタによりノイズを濾した呼吸数をそれぞれ推定する第２のステップと、
   　この第２のステップで得られた複数の呼吸数の推定値の重み付け平均化処理を実行する第３のステップとを含むことを特徴とする呼吸推定方法。
2. 　請求項１記載の呼吸推定方法において、
   　前記第１のステップは、
   　前記第１のデータの時系列信号と前記第２のデータの時系列信号の各々をフーリエ変換することにより、前記第１のデータと前記第２のデータの各々の信号の周波数スペクトルをそれぞれ求めるステップと、
   　前記第１のデータから求めた周波数スペクトルと前記第２のデータから求めた周波数スペクトルの各々から前記被験者の呼吸数をそれぞれ推定するステップとを含むことを特徴とする呼吸推定方法。
3. 　請求項１または２記載の呼吸推定方法において、
   　前記第１のデータと前記第２のデータの各々の信頼度を示すエラー係数をそれぞれ測定する第４のステップと、
   　前記エラー係数に基づいて重み付け行列の成分である重み付け定数を求める第５のステップと、
   　前記カルマンフィルタの推定誤差が最小となるように前記カルマンフィルタの共分散推定誤差行列を更新する第６のステップと、
   　前記重み付け行列と前記共分散推定誤差行列とに基づいて前記カルマンフィルタのカルマン利得を更新する第７のステップとをさらに含むことを特徴とする呼吸推定方法。
4. 　請求項３記載の呼吸推定方法において、
   　前記第７のステップは、測定系雑音の共分散行列を前記重み付け行列に応じて調整し、調整後の共分散行列を用いて前記カルマン利得を算出するステップを含むことを特徴とする呼吸推定方法。
5. 　請求項１乃至４のいずれか１項に記載の呼吸推定方法において、
   　前記第３のステップは、
   　前記第１のデータと前記第２のデータの各々に関する前記カルマンフィルタの自乗推定誤差から前記重み付け平均化処理のための重み付け定数を算出するステップと、
   　このステップで算出した重み付け定数を用いて、前記第１のデータに基づく呼吸数の推定値と前記第２のデータに基づく呼吸数の推定値とを重み付け平均化するステップとを含むことを特徴とする呼吸推定方法。
6. 　請求項１乃至５のいずれか１項に記載の呼吸推定方法において、
   　前記第１のデータとして、前記被験者の心電位波形からＲ波の振幅を検出する第１の検出ステップと、
   　前記第１のデータとして、前記被験者の心電位波形からＲ波と１つ前のＲ波との間隔であるＲＲ間隔を検出する第２の検出ステップと、
   　前記第２のデータとして、前記被験者の呼吸運動による３軸加速度の信号から加速度ベクトルの角度変位を検出する第３の検出ステップとをさらに含むことを特徴とする呼吸推定方法。
7. 　被験者の心機能に関する第１のデータの時系列信号と前記被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から前記被験者の呼吸数をそれぞれ推定する呼吸数推定部と、
   　前記第１のデータから推定した呼吸数と前記第２のデータから推定した呼吸数の各々についてノイズを濾した呼吸数をそれぞれ推定するカルマンフィルタと、
   　このカルマンフィルタで得られた複数の呼吸数の推定値の重み付け平均化処理を実行する統合処理部とを備えることを特徴とする呼吸推定装置。
8. 　被験者の心機能に関する第１のデータの時系列信号と前記被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から位相情報をそれぞれ抽出する第１のステップと、
   　前記第１のデータの位相情報と前記第２のデータの位相情報の各々についてカルマンフィルタによりノイズを濾した位相情報をそれぞれ推定する第２のステップと、
   　この第２のステップで得られた複数の推定位相値の重み付け平均化処理を実行する第３のステップと、
   　この第３のステップで統合された位相値を周波数に変換することにより、前記被験者の呼吸周波数を求める第４のステップとを含むことを特徴とする呼吸推定方法。
9. 　請求項８記載の呼吸推定方法において、
   　前記第１のステップは、
   　前記第１のデータと前記第２のデータの各々について、時系列信号をヒルベルト変換して実部成分と虚部成分の２つの信号成分をそれぞれ生成するステップと、
   　前記第１のデータと前記第２のデータの各々について、前記実部成分と前記虚部成分とから角度をそれぞれ算出するステップと、
   　前記第１のデータと前記第２のデータの各々について、前記角度を連続した位相値にする位相アンラップをそれぞれ行うステップとを含むことを特徴とする呼吸推定方法。
10. 　請求項８または９記載の呼吸推定方法において、
    　前記第３のステップは、
    　前記第１のデータと前記第２のデータの各々に関する前記カルマンフィルタの自乗推定誤差から重み付け定数を算出するステップと、
    　このステップで算出した重み付け定数を用いて、前記第１のデータの推定位相値と前記第２のデータの推定位相値の重み付け平均化処理を実行するステップとを含むことを特徴とする呼吸推定方法。
11. 　請求項８乃至１０のいずれか１項に記載の呼吸推定方法において、
    　前記第１のデータの時系列信号と前記第２のデータの時系列信号の各々を、心電位波形のサンプリング周波数および３軸加速度信号のサンプリング周波数よりも遅いサンプリング周波数でそれぞれサンプリングする第５のステップと、
    　前記第５のステップで取得した前記第１のデータの時系列信号と前記第２のデータの時系列信号の各々をそれぞれ帯域制限する第６のステップとをさらに含み、
    　前記第１のステップは、前記第６のステップで帯域制限した前記第１のデータの時系列信号と前記第２のデータの時系列信号の各々から位相情報をそれぞれ抽出するステップを含むことを特徴とする呼吸推定方法。
12. 　請求項８乃至１１のいずれか１項に記載の呼吸推定方法において、
    　前記第１のデータとして、前記被験者の心電位波形からＲ波の振幅を検出する第１の検出ステップと、
    　前記第１のデータとして、前記被験者の心電位波形からＲ波と１つ前のＲ波との間隔であるＲＲ間隔を検出する第２の検出ステップと、
    　前記加速度データとして、前記被験者の呼吸運動による３軸加速度の信号から加速度ベクトルの角度変位を検出する第３の検出ステップとをさらに含むことを特徴とする呼吸推定方法。
13. 　被験者の心機能に関する第１のデータの時系列信号と前記被験者の呼吸運動による加速度に関する第２のデータの時系列信号の各々から位相情報をそれぞれ抽出する特徴量抽出部と、
    　この特徴量抽出部で得られた前記第１のデータの位相情報と前記第２のデータの位相情報の各々についてノイズを濾した位相情報をそれぞれ推定するカルマンフィルタと、
    　このカルマンフィルタで得られた前記第１のデータの推定位相値と前記第２のデータの推定位相値の重み付け平均化処理を実行する統合処理部と、
    　この統合処理部で統合された位相値を周波数に変換することにより、前記被験者の呼吸周波数を求める呼吸周波数変換部とを備えることを特徴とする呼吸推定装置。

Images