JP2019051011A - 瞬時心拍の時系列データの補完装置、補完方法及びそのプログラム - Google Patents

Abstract

【課題】 計測異常などによって生じた欠損区間のある瞬時心拍データであっても、適切なスペクトル解析を実現することができる、瞬時心拍の時系列データの補完装置を提供する。【解決手段】 実施形態の瞬時心拍の時系列データの補完装置は、被験者の心電を計測する心電計測部から出力される信号に基づいて算出される、時系列で隣接する２つのＲ波の間隔である瞬時心拍の内、正常な値と判断した複数の瞬時心拍からなる時系列データを入力とし、前記入力された前記時系列データの周波数成分に基づく補完関数によって補完値を算出し（Ｓ６１〜Ｓ６２）、前記算出した前記補完値により、前記入力された前記時系列データにおいて欠損した瞬時心拍を補完する（Ｓ６３〜Ｓ６６）、瞬時心拍補完部を具備する。【選択図】図６

Description

本発明は、瞬時心拍の時系列データの補完装置、補完方法及びそのプログラムに関する。

１ 心拍変動解析に関する概要
自律神経には、交感神経と迷走神経の二つがある。両自律神経は各臓器などに広く分布し、循環や代謝をはじめとする不随意な身体機能を制御する。多くの場合、両自律神経がひとつの臓器を拮抗的に支配すると言われている。

自律神経活動のひとつである交感神経活動は、暗算負荷などのストレス刺激によって亢進することが知られている。

もう一つの自律神経である迷走神経は、当該神経が支配する各臓器において主に副交換性の神経活動を担うことから、副交感神経活動と同等に理解されることも多い。なお、「迷走神経」とは、厳密には脳神経のひとつである第Ｘ神経の名称であり、脳から各臓器などに至る当該神経すべてを指す。このため、支配対象となっている臓器の名称を付記することで、対象臓器における副交感神経活動を示す場合がある（例：心臓迷走神経）。

自律神経が支配する臓器のひとつに心臓がある。心臓は交感神経、迷走神経によって拮抗的に支配されており、両自律神経活動を反映すると言われている（参考文献［ｉ］参照）。

特に、隣接する二つのＲ波の間隔である瞬時心拍（ＲＲＩ：Ｒ−Ｒ ｉｎｔｅｒｖａｌ）のゆらぎは両自律神経活動によって変化することが知られている。なお、Ｒ波とは、心電図計測によって得られる心電波形のひとつであり、心臓の脱分極活動を反映している（参考文献［ｉｉ］参照）。図９は、Ｒ波と瞬時心拍（ＲＲＩ）との関係を示す図である。

瞬時心拍は、直流成分、低周波数成分、高周波数成分の三つの成分を含むと言われている（参考文献［ｉ］参照）。なお、低周波数成分、高周波数成分については、それぞれの帯域である０．０４〜０．１５Ｈｚ、０．１５〜０．４０Ｈｚの各周波数成分の合計からなる。

実環境で自律神経活動を推定する手法として、瞬時心拍変動の周波数スペクトル解析がある。この手法によれば、不等間隔である瞬時心拍を周波数スペクトルで解析した際の低周波数成分（以降ＨＲＶ_ＬＦ）は交感神経活動と心臓迷走神経活動、高周波数成分（以降ＨＲＶ_ＨＦ）は心臓迷走神経活動を反映する指標として解釈される（参考文献［ｉ］参照）。

心電図を計測する手段のひとつとして、ホルター心電計などのウェアラブルなデバイスがある。これらのデバイスを用いて取得する心電図は、電極の変形やズレをはじめとする電極異常、あるいは、体動、発汗、静電気など様々な要因によって計測異常が生じる（参考文献［ｉｉｉ］参照）。この計測異常は、心電図では図１０に示すようなアーチファクトやノイズという形で確認できる。

計測異常のうち、アーチファクトとして観測される波形はＲ波と類似しており、ひとつ以上連続して観測される。このため、心電図を解析してＲ波を抽出するアルゴリズムが、アーチファクトをＲ波と誤判断して抽出してしまう場合もある（以降では、このようなものを「計測異常Ｒ波」、計測異常Ｒ波によって誤算出した瞬時心拍を「計測異常瞬時心拍」と表記する）。

ＨＲＶ_ＬＦおよびＨＲＶ_ＨＦは、解析対象となるすべてのデータが正常な瞬時心拍である場合でのみ、自律神経活動を反映する。ここでいう正常な状態とは、計測対象と計測器両方において異常がない状態を意味する。計測対象の異常とは不整脈などであり、計測器の異常とは心電図で計測異常が生じている状態を指す。

ノイズやアーチファクトをＲ波と誤判断したものは、その発生機序から心臓の脱分極活動を一切反映しない。このため、解析対象となる瞬時心拍を構成するＲ波のうち、少なくともひとつがアーチファクトをＲ波と誤判断したものである場合、ＨＲＶ_ＬＦおよびＨＲＶ_ＨＦのいずれも自律神経活動を反映するとは言えない。

ノイズやアーチファクトを誤判定したものを含む瞬時心拍と、正常に計測した瞬時心拍とを、周波数解析で区別することは原理的に不可能である。そのため、周波数解析を行う前に、誤判定した瞬時心拍を除去する必要がある。

２ 心拍の周波数特徴量の算出
心拍特徴量のうち、周波数特徴量を求める際には、不等間隔である瞬時心拍を補間関数によって等間隔なデータにリサンプリングしてからパワースペクトル密度を求める必要がある（参考文献［ｉ］参照）。リサンプリング手法には線形補間やスプライン補間などがあるが、リサンプリング後のデータに心拍と似たゆらぎを残せるスプライン補間が用いられる場合が多い。

なお、スプライン補間においては、心拍のゆらぎである心拍変動から心拍特徴量を計算するため、一般に３次スプライン関数が用いられる。すなわち、スプライン関数を用いた補間そのものとしては、例えば平滑化スプライン関数など他の手法を使うことも手段としては考えられる。しかし、平滑化スプライン関数は、全計測データと補間関数との間の誤差が最も少なくなるように補間関数を決定するため、補間関数が実際の計測データを通らない場合がある。このような場合、平滑化スプライン関数によって求まった補間関数は心拍変動を正しく反映しているとはいえないため、心拍特徴量の算出には３次スプライン関数が好まれる傾向にある。

出来る限り正常な状態の瞬時心拍のみからリサンプリングデータを生成するため、一般的には、補間処理の前に、計測異常瞬時心拍つまり瞬時心拍の異常値の除外を行う。瞬時心拍の異常値を除外する手法として、瞬時心拍の時間特徴量を使用するものがある。具体的な例としては、瞬時心拍の下限値・上限値、および前後の瞬時心拍との差分値に閾値を設定し、当該閾値を逸脱したものを除外する手法（参考文献［ｉｖ］参照）や、瞬時心拍の正規分布から外れたものを除外する手法（参考文献［ｖ］参照）がある。後者の手法（参考文献［ｖ］参照）の中では、瞬時心拍の平均±標準偏差による異常値検出が最も簡便なものであり、一般的には２σあるいは３σルールが用いられることが多い。

周波数特徴量は欠損値の影響を受けて大きく変動することが知られている（参考文献［ｉ］参照）。例えば、高周波数成分と低周波数成分の比であるＨＲＶ_{ＬＦ／ＨＦ}は、欠損率６％で１．５倍、１２％で２倍近くまで増加する場合があるとの報告がある。

上記のスペクトル解析を行う場合は、各瞬時心拍データを後方のＲ波の時間的位置にプロットしたタコグラムが必要となる（参考文献［ｉ］参照）。このタコグラムの縦軸は計測した瞬時心拍の値であるが、横軸の取り方は過去の瞬時心拍の累積時間（手法ａ）を使う場合と、瞬時心拍の生起時間を使う場合（手法ｂ）の二つがある。手法ａ、ｂの時間は、計測した瞬時心拍がすべて正常であり、異常値除外もない場合は等しくなる。

すなわち、異常値除外が生じ得る状況において、ＲＲＩタコグラムを上記手法ａに基づいて作成すると、本来解析対象とすべきタコグラムを得られなくなる恐れがある。このため、基本的には、理論通りのタコグラムを得ることが可能な上記手法ｂが用いられる場合が多い。

３ ウェアラブル心電計を用いた心拍変動解析
ウェアラブル心電計を用いる場合、一般には長時間に及ぶ心拍変動解析を行う場合が多い。このような場合には、一定長の解析窓幅と解析シフト幅を設定し、対象データを時間経過に応じて変更しながら解析を行う。

従来、参考文献［ｉｖ］に示した手法や参考文献［ｖ］に示した手法で瞬時心拍の異常値除外を行なったデータをスプライン補間してパワースペクトル密度を算出する方法がある（非特許文献１）。

さらに別の従来手法として、参考文献［ｉｖ］に示した手法や参考文献［ｖ］に示した手法で瞬時心拍の異常値除外を行った後、データの補間処理をする前に、正常に計測できた瞬時心拍の平均値を用いて欠損区間を補完する方法がある（非特許文献２、非特許文献３）。

神田尚子、佐久間大輝、吉永努、入江英嗣，色彩環境下での心拍変動と作業能率との相関に関する検討，インタラクティブシステムとソフトウェアに関するワークショップ，ｐｐ．２３１−２３２，２０１２ 青木良輔、江口佳那、島内末廣、千葉昭宏、武田十季，日常環境における人の内部状態推定のための着衣型ウェアラブル心電計のデザイン検討，信学技報，Ｖｏｌ．１１７，Ｎｏ．１０９，ｐｐ．１２７−１３２，２０１７（公知日２０１７年６月１６日、発表日２０１７年６月２４日） 江口佳那、青木良輔、吉田和広、山田智広，ウェアラブル心電計から計測された心拍の周波数特徴量解析を対象とした欠損ＲＲＩの補完手法，マルチメディア、分散、協調とモバイル（ＤＩＣＯＭＯ２０１７）シンポジウム，ｐｐ．８８８−８９７，２０１７（公知日２０１７年６月２１日、発表日２０１７年６月２９日）

上記手法ｂでＲＲＩタコグラムを作成し、従来手法のあと欠損値がある状態でスプライン補間を行う場合、欠損区間の長さによってはスプライン曲線が生理的に取りえない値となってしまう場合がある（例えば、２５０［ｍｓｅｃ］以下の値や負値など）。このようなリサンプリングデータから得たパワースペクトル密度は、本来取得し得ないタコグラムから算出されたものであるため、心臓の活動を正しく反映しているとは言えない。

さらに、ウェアラブル心電計で取得した瞬時心拍には、計測異常瞬時心拍が含まれることもある。このような計測異常瞬時心拍を異常値として除外すると、瞬時心拍の時間特徴量のみに基づく異常値除外よりも多くの瞬時心拍を異常値として除外することになる。したがって、通常の計測機器以上に解析対象データに欠損値が含まれる可能性が高いと考えられる。

また、上記非特許文献２、３に記載の手法によれば、上述のようなスプライン曲線が生理的に取りえない値域に発振することを防ぐことができる。しかし、上記非特許文献２、３に記載の手法では、補完に用いる値が一定であり、正常に計測した瞬時心拍が持つゆらぎ成分がないため、補完後の瞬時心拍を解析して得られる周波数特徴量が、本来取るべき値から乖離する場合がある。

本発明は、上記実情に鑑みてなされたものであり、計測異常などによって生じた欠損区間のある瞬時心拍データであっても、適切なスペクトル解析を実現することができる、瞬時心拍の時系列データの補完装置、補完方法及びそのプログラムを提供することを目的とする。

本発明の第１の態様は、被験者の心電を計測する心電計測手段から出力される信号に基づいて算出される、時系列で隣接する２つのＲ波の間隔である瞬時心拍の内、正常な値と判断した複数の瞬時心拍からなる時系列データを入力とし、前記入力された前記時系列データの周波数成分に基づく補完関数によって補完値を算出し、前記算出した前記補完値により、前記入力された前記時系列データにおいて欠損した瞬時心拍を補完する瞬時心拍補完部を具備するようにしたものである。

本発明の第２の態様によれば、第１の態様において、前記正常な値と判断した複数の瞬時心拍からなる前記時系列データのうち連続した瞬時心拍のデータ長に基づいて、前記時系列データにおいて欠損した前記瞬時心拍を補完するための前記補完関数の項を決定する瞬時心拍補完関数判別部を更に具備し、前記瞬時心拍補完部は、前記瞬時心拍補完関数判別部で決定した前記補完関数の前記項に応じて前記補完関数の各項について算出されたパラメータを適用した前記補完関数によって前記補完値を算出し、前記算出した前記補完値により、前記入力された前記時系列データにおいて欠損した瞬時心拍を補完する。

本発明の第３の態様によれば、第１または第２の態様において、前記瞬時心拍補完部における前記補完関数を構成する項は、直流成分、低周波数成分、高周波数成分のいずれか一つ以上である。

本発明の第４の態様によれば、瞬時心拍の時系列データの補完方法として、被験者の心電を計測する心電計測手段から出力される信号に基づいて算出される、時系列で隣接する２つのＲ波の間隔である瞬時心拍の内、正常な値と判断した複数の瞬時心拍からなる時系列データを入力し、前記入力した前記時系列データの周波数成分に基づく補完関数によって補完値を算出し、前記算出した前記補完値により、前記入力した前記時系列データにおいて欠損した瞬時心拍を補完するようにしたものである。

本発明の第５の態様によれば、第４の態様において、前記正常な値と判断した複数の瞬時心拍からなる前記時系列データのうち連続した瞬時心拍のデータ長に基づいて、前記時系列データにおいて欠損した前記瞬時心拍を補完するための前記補完関数の項を決定し、前記瞬時心拍を補完することは、前記決定した前記補完関数の前記項に応じて前記補完関数の各項について算出されたパラメータを適用した前記補完関数によって前記補完値を算出し、前記算出した前記補完値により、前記入力した前記時系列データにおいて欠損した瞬時心拍を補完する。

本発明の第６の態様によれば、第４または第５の態様において、前記正常な値と判断した複数の瞬時心拍からなる前記時系列データのうち連続した瞬時心拍のデータ長に基づいて、前記時系列データにおいて欠損した前記瞬時心拍を補完するための前記補完関数の項を決定し、
前記瞬時心拍を補完することは、前記決定した前記補完関数の前記項に応じて前記補完関数の各項について算出されたパラメータを適用した前記補完関数によって、前記入力した前記時系列データにおいて欠損した瞬時心拍を補完する。

本発明の第７の態様によれば、コンピュータによって実行されたときに、前記コンピュータを、第１の態様乃至第３の態様のいずれか一つの態様の前記瞬時心拍の時系列データの補完装置として機能させるためのプログラムである。

本発明によれば、計測異常などによって生じた欠損区間のある瞬時心拍データであっても、適切なスペクトル解析を実現することができる、瞬時心拍の時系列データの補完装置、補完方法及びそのプログラムを提供することができる。

本発明の実施形態に係る瞬時心拍の時系列データの補完装置１２を適用した補完システム１０の構成を示す図である。 補完装置１２の瞬時心拍評価部２３による瞬時心拍評価及び異常値除外と瞬時心拍再計算部２６による瞬時心拍再計算の一例を示す図である。 補完システム１０の動作を説明するためのフローチャートを示す図である。 Ｓ５における瞬時心拍補完関数判別部２７の補完関数の項数または項判別処理を説明するためのフローチャートである。 Ｓ５１における対象データ取得処理における区間データの取得方法を説明するための図である。 Ｓ６における瞬時心拍補完部２８の瞬時心拍の時系列データに対する補完処理を説明するためのフローチャートである。 瞬時心拍電位の時系列データと内挿された補完値との関係を示す図である。 次の瞬時心拍との時間間隔が短くなる場合に補完値を内挿しないことを説明するための図である。 図９は、Ｒ波と瞬時心拍（ＲＲＩ）との関係を示す図である。 図１０は、心電図における計測異常（アーチファクト、ノイズ）の例を示す図である。

以下、図面を参照して本発明の実施形態に係る瞬時心拍の時系列データの補完装置を適用した補完システムについて説明する。

図１は、瞬時心拍の時系列データの補完システム１０の構成を示す図である。瞬時心拍の時系列データの補完システム１０は、心電計測部１１と、本発明の実施形態に係る瞬時心拍の時系列データの補完装置１２と、を具備する。

一例として、瞬時心拍の補完システム１０は、心電計測部１１を被験者（ユーザ）に装着可能なウェアラブルデバイスとし、瞬時心拍の時系列データの補完装置１２をスマートフォン、タブレット型端末、パーソナルコンピュータ（ＰＣ）などのコンピュータデバイスとしたシステムにより実現される。例えば、コンピュータデバイスは、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）などのプロセッサと、プロセッサに接続されるメモリと、心電計測部１１と（例えば無線で）通信するための通信インタフェースと、を備える。

なお、瞬時心拍の補完システム１０の実現形態は、この例に限定されるものではない。
例えば、瞬時心拍の補完システム１０は１つのデバイスとして実現されても良い。また、心電計測部１１は瞬時心拍の補完システム１０の外部に設けられても良い。言い換えると、瞬時心拍の補完システム１０は、心電計測部１１に相当する外部の心電計測装置から被験者の心電を計測した結果を取得しても良い。

心電計測部１１は、被験者の心電を計測し、計測結果をＲ波抽出部２１に送る。心電は、循環器系の生体信号であり、例えば、心室の収縮と同期した周期的な信号を含む。心電計測部１１は、少なくとも２極の電極によって心電の計測を行う。計測結果は、心電図におけるＲ波相当の心電を抽出可能な時系列データを含む。

例えば、計測結果は心電図の時系列データを含む。心電計測部１１は、Ｒ波相当の心電を計測することができれば良く、その実現形態は問わない。例えば、心電計測部１１はホルター心電計からなる。

瞬時心拍の時系列データの補完装置１２は、Ｒ波抽出部２１、Ｒ波関連情報記録部２２、瞬時心拍算出部２３、瞬時心拍記録部２４、瞬時心拍評価部２５、瞬時心拍再計算部２６、瞬時心拍補完関数判別部２７、瞬時心拍補完部２８及びスペクトル解析処理部２９を具備する。これらＲ波抽出部２１、Ｒ波関連情報記録部２２、瞬時心拍算出部２３、瞬時心拍記録部２４、瞬時心拍評価部２５、瞬時心拍再計算部２６、瞬時心拍補完関数判別部２７、瞬時心拍補完部２８及びスペクトル解析処理部２９の機能は、例えば、プロセッサがメモリに格納されているプログラムを読み出して実行することにより実現される。なお、これらの機能の一部または全部は、特定用途向け集積回路（ＡＳＩＣ）などの回路によって実現されても良い。

Ｒ波抽出部２１は、心電計測部１１によって計測された心電の時系列データを解析し、Ｒ波を抽出する。実施形態では、具体的なＲ波の抽出方法は問わない。後続処理で必要がある場合については、抽出したＲ波に関連する情報を、Ｒ波関連情報記録部２２に記録する。

Ｒ波関連情報記録部２２は、瞬時心拍評価部２５において、公知の技術によって瞬時心拍以外の異常値除外処理を行う場合などにおいて、Ｒ波抽出部２１で抽出したＲ波に関する情報を記録する。実施形態では、具体的な記録方法については特に指定しない。

瞬時心拍算出部２３は、Ｒ波抽出部２１で抽出したＲ波に基づき、瞬時心拍を算出し、算出した瞬時心拍情報を瞬時心拍記録部２４に保存する。

瞬時心拍記録部２４は、瞬時心拍算出部２３で算出した瞬時心拍を記録する。具体的な記録形式については特に指定しないが、例えば、瞬時心拍の行列や、瞬時心拍を構成する１つ目のＲ波の時刻情報と瞬時心拍の二つから構成されるデータ行列が考えられる。なお、本機能は実施形態において必須の機能ではない。Ｒ波の電位情報に加え、瞬時心拍の時間情報も考慮して瞬時心拍を評価する場合のみ必要となる。

瞬時心拍評価部２５は、Ｒ波関連情報記録部２２の情報を基に、瞬時心拍算出部２３で算出した瞬時心拍を評価し、異常値除外処理を行う。実施形態では具体的な処理内容については規定しない。例えば、参考文献［ｉｖ］、［ｖ］の手法ように、瞬時心拍の時間特徴のみに基づく異常値除外処理のみでも良い。

また、瞬時心拍再計算部２６は、実施形態においては任意の処理であるが、例えば、以下の処理が含まれる。

具体的には、瞬時心拍再計算部２６は、瞬時心拍評価部２５による評価及び異常値除外処理の結果に基づいて瞬時心拍の再計算を行う。図２は本発明の実施形態における補完装置１２の瞬時心拍評価部２３による瞬時心拍評価及び異常値除外と瞬時心拍再計算部２６による瞬時心拍再計算の一例を示す図である。

例えば、瞬時心拍評価部２５は、表１で示す＃２、＃３および＃４のいずれかに該当すると評価された瞬時心拍を除去し、このような異常値除外処理を行った瞬時心拍の時系列データを瞬時心拍補完関数判別部２７及び瞬時心拍補完部２８に出力する。また、瞬時心拍再計算部２６へは、評価結果を出力する。瞬時心拍再計算部２６は、この瞬時心拍評価部２５による評価結果に基づいて、瞬時心拍評価部２５によって除去された瞬時心拍を構成する、隣接する２つのＲ波における、アーチファクトとして判別されたＲ波を棄却し、残りの正常計測状態と判別されたＲ波のみを用いて、隣接する２つのＲ波である瞬時心拍を構成する。

図２に示した例では、再計算前のＲ波は、時系列に沿って、Ｒ波番号「１」、「２」の「Ｒ（正常計測状態）」、Ｒ波番号「３」の「Ａ（アーチファクト）」、Ｒ波番号「４」の「Ｒ」、Ｒ波番号「５」、「６」、「７」の「Ａ」、Ｒ波番号「８」、「９」の「Ｒ」と判別されたＲ波である。
これに対し、瞬時心拍再計算部２６は、瞬時心拍評価部２５によってＲ波番号「３」、「５」、「６」、「７」の「Ａ」と判別されたＲ波を棄却し、残りのＲ波番号「１」、「２」、「４」、「８」、「９」の「Ｒ」と判別されたＲ波を用いて、１つ目の瞬時心拍（隣接する、Ｒ波番号「１」、「２」に対応する２つのＲ波で構成）、２つ目の瞬時心拍（隣接する、Ｒ波番号「２」、「４」に対応する２つのＲ波で構成）、３つ目の瞬時心拍（隣接する、Ｒ波番号「４」、「８」に対応する２つのＲ波で構成）、４つ目の瞬時心拍（隣接する、Ｒ波番号「８」、「９」に対応する２つのＲ波で構成）を構成する。これら再計算後の、瞬時心拍評価部２５による、瞬時心拍の計測状態の評価は、表１に示す「＃１（２つとも正常計測状態）」のみとなる。

なおこのとき、再計算後の瞬時心拍が洞調律の範囲以下（例えば１５００ｍｓ以下）であることを、再計算結果を以降の処理において使用する最低条件としても良い。すなわち、瞬時心拍再計算部２６は、Ｒ波とＲ波の間に瞬発的に生じた異常のみを再計算対象とし、それ以上の値については計測異常として除外する。

瞬時心拍補完関数判別部２７は、瞬時心拍評価部２５で異常値除外処理を行った或いは、瞬時心拍再計算部２６で再計算が行なわれた瞬時心拍の時系列データに基づいて、補完処理に用いる補完関数の項数または項を判別する。具体的な手法については後述する。

瞬時心拍補完部２８は、瞬時心拍補完関数判別部２７で判定した補完関数の項数または項に応じて補完関数を決定し、瞬時心拍評価部２５で瞬時心拍の異常値除外処理を行った或いは、瞬時心拍再計算部２６で再計算が行なわれた瞬時心拍の時系列データを対象として、この補完関数を用いた補完処理を行う。具体的な手法については後述する。

スペクトル解析処理部２９は、瞬時心拍補完部２８で補完処理を行った瞬時心拍の時系列データを対象としてスペクトル解析を行う。実施形態では具体的な処理については特に指定しないが、例えば参考文献［ｉ］に記載のように、スプライン補間によってリサンプリングデータを取得した後、自己回帰モデルによってパワースペクトル密度を求める処理などが考えられる。

次に、実施形態に係る瞬時心拍の補完システムの動作について、図３のフローチャートを参照して説明する。なお、実施形態では、参考文献［ｉｖ］、［ｖ］の手法に基づく瞬時心拍評価のみを行なう場合を仮定する。

まず、心電計測部１１は、被験者の心電を計測して、被験者の心電図のための時系列の心電データをＲ波抽出部２１に送る（Ｓ１）。Ｒ波抽出部２１は、心電計測部１１にて計測された時系列の心電データからＲ波を抽出する（Ｓ２）。

瞬時心拍算出部２３は、Ｒ波抽出部２１で取得したＲ波に基づき、隣接する２つのＲ波から瞬時心拍を算出する（Ｓ３）。実施形態では、参考文献［ｉｖ］、［ｖ］の手法により、瞬時心拍評価部２５で瞬時心拍の時間情報を用いた評価を行うため、算出した瞬時心拍に関する情報を瞬時心拍記録部２４に記録するが、瞬時心拍の時間情報を異常値除外処理に用いない場合は記録しなくても良い。

瞬時心拍評価部２５は、Ｒ波関連情報記録部２２に保存したＲ波の情報や、瞬時心拍記録部２４の時間情報に基づいて、瞬時心拍を評価する（Ｓ４）。実施形態では、参考文献［ｉｖ］、［ｖ］の手法に基づき２５０［ｍｓｅｃ］未満あるいは１５００［ｍｓｅｃ］を越える瞬時心拍を除外した後、平均±３×標準偏差の範囲から外れた瞬時心拍を異常値として除外する。

瞬時心拍補完関数判別部２７は、瞬時心拍評価部２５または瞬時心拍再計算部２６の処理で正常な値と判断された時系列データに基づいて、補完処理に用いる補完関数の項数または項を判別する（Ｓ５）。具体的な判別処理の例は、後述する。

瞬時心拍補完部２８は、瞬時心拍評価部２５で瞬時心拍の異常値除外処理を行った或いは、瞬時心拍再計算部２６で再計算が行なわれた瞬時心拍の時系列データを対象として、補完関数を用いた補完処理を行う（Ｓ６）。このとき、補完に用いる補完関数の項数または項は、Ｓ５で判別した項数または項とする。具体的な補完処理の例は、後述する。

その後、スペクトル解析処理部２９により、所定の条件のもとに補完処理が行なわれた瞬時心拍の時系列データに対してスペクトル解析が行なわれる（Ｓ７）。

図４は、Ｓ５における瞬時心拍補完関数判別部２７の補完関数の項数または項判別処理を説明するためのフローチャートである。実施形態では、正常な値と判断された複数の瞬時心拍からなる時系列データの内、最も長いデータを用いる場合について説明する。

同図に示すように、瞬時心拍補完関数判別部２７は、まず、瞬時心拍評価部２５または瞬時心拍再計算部２６の処理で正常な値と判断された複数の瞬時心拍からなる時系列データの内、最も長いデータを対象データとして取得する（Ｓ５１）。ここで、対象データとは、複数個の連続した正常な瞬時心拍によって構成される瞬時心拍の総和（累積時間であって、瞬時心拍の個数ではない）が最も大きい区間のデータのことである。例えば、心電計測部１１によって計測される心電の時系列データが図５８のようであった場合には、補完対象区間は「異常データ区間＃１」及び「異常データ区間＃２」となり、これら異常データ区間の前後の「正常データ区間＃１」、「正常データ区間＃２」、「正常データ区間＃３」が、瞬時心拍評価部２５または瞬時心拍再計算部２６の処理で正常な値と判断された複数の瞬時心拍からなるデータ区間、つまり正常データ区間となる。この例では、瞬時心拍補完関数判別部２７は、最も長いデータである「正常データ区間＃３」を、対象データとして取得する。

なお、ここでは、最も長いデータを対象データとして取得するものとしたが、本実施形態はこの例に限定されるものではない。データの選択取得の方法としては、補完対象区間の直近の正常計測区間のデータを取得する方法もある。例えば、補完対象区間である異常データ区間の直前の正常計測区間のデータを採用しても良い。すなわち、補完対象区間「異常データ区間＃１」に対しては「正常データ区間＃１」のデータ、補完対象区間「異常データ区間＃２」に対しては「正常データ区間＃２」のデータ、を採用することができる。あるいは、補完対象区間の直後の正常データ区間のデータを取得するようにしても良いし、補完対象区間の前後の正常データ区間のうち最も長いデータを採用するようにしても良い。また、補完対象の区間毎に上記データ選択方法を組合せたり、上記のデータ選択方法を用いて決定された複数の補完関数から適宜選択したりするようにしても良い。

そして、この取得した対象データのデータ長が、高周波数成分が必要とするデータ長以上である、という第１の条件を充足するか判別する（Ｓ５２）。本実施形態では、この第１の条件の充足判別の具体的な手法については特に規定しないが、例えば、高周波数成分の必要データ長を０．１５Ｈｚの１／２周期である３秒、また後述する低周波数成分の必要データ長を０．０４Ｈｚの１／２周期である１２．５秒とすることなどが考えられる。

上記第１の条件を充足していると判別した場合には、瞬時心拍補完関数判別部２７は、補完に用いる関数である補完関数の項数または項を決定する。まず、対象データのデータ長が、低周波数成分が必要とするデータ長以上である、という第２の条件を充足するか判別し（Ｓ５３）、その第２の条件を充足している場合には、瞬時心拍補完部２８が補完処理に用いる補完関数は、直流成分、低周波数成分、高周波数成分の三つの成分を項に持つものと決定する（Ｓ５４）。また、Ｓ５３で上記第２の条件を充足していないと判別した場合には、その第２の条件を満たさない対象データについて、そのデータ長が、低周波数成分が必要とするデータ長未満、且つ、高周波数成分が必要とするデータ長以上である、という第３の条件を充足するか判別する（Ｓ５５）。この第３の条件を充足している場合には、瞬時心拍補完部２８が補完処理に用いる補完関数は、直流成分、高周波数成分の二つの成分を項に持つものと決定する（Ｓ５６）。そして、Ｓ５５で上記第３の条件を充足していないと判別した場合或いは、上記Ｓ５２で上記第１の条件を充足していないと判別した場合には、その第３または第１の条件を満たさない対象データについては、瞬時心拍補完部２８が補完処理に用いる補完関数は、直流成分のみを項に持つものと決定する（Ｓ５７）。

なお、上記Ｓ５２における上記第１の条件として、対象データのデータ長が高周波数成分が必要とするデータ長以上であることを判別している。そのため、上記Ｓ５３において上記第２の条件（対象データのデータ長が低周波数成分が必要とするデータ長以上である）を満たしていなければ、必然的に、対象データのデータ長は、低周波数成分が必要とするデータ長未満且つ高周波数成分が必要とするデータ長以上である。よって、Ｓ５５での上記第３の条件の充足判別は省略し、Ｓ５３で上記第２の条件を充足していないと判別した場合にはＳ５６へと進んで、瞬時心拍補完部２８が補完処理に用いる補完関数は、直流成分、高周波数成分の二つの成分を項に持つものと決定するようにしても構わない。

図６は、Ｓ６における瞬時心拍補完部２８の欠損区間を有する瞬時心拍の時系列データに対する補完処理を説明するためのフローチャートである。

同図に示すように、瞬時心拍補完部２８は、まず、補完に用いる補完関数を算出する（Ｓ６１）。この補完関数の算出処理の詳細について、以下、瞬時心拍補完関数判別部２７がＳ５４の補完処理に用いる補完関数は、直流成分、低周波数成分、高周波数成分の三つの成分を項に持つものと決定した場合を例にして説明する。

この場合、瞬時心拍補完部２８は、瞬時心拍が直流成分、低周波数成分、高周波数成分の三つの成分からなることに着目し、正常な値と判断された複数の瞬時心拍からなる時系列データ、すなわち正常に計測した瞬時心拍から、これらの成分を算出する。本実施形態では、各成分算出に用いるデータ区間については具体的に規定しないが、例えば、正常に計測した瞬時心拍の中で最も長いデータなどが考えられる。

なお、直流成分、低周波数成分、高周波数成分の各成分の構成方法についても、本実施形態で詳細には規定しないが、一例として、各成分の代表値を用いる方法がある。この方法によれば、各成分に含まれ得る誤判定した瞬時心拍の影響を低減しつつ、各成分の大局的な変動を捉えることが可能となる。

また、本実施形態では、各成分の算出方法についても具体的に規定しない。一例としては、対象データの直流成分（０Ｈｚ）または対象データの平均値を直流成分とし、対象データから直流成分を除いたデータをケプストラム解析して得られるピーク周波数を用いる方法がある。この場合、少なくとも二つの大きさの異なるケプストラム次数を用いてケプストラム解析を行い、大きいケプストラム次数を用いた場合のピーク周波数を低周波数成分、小さいケプストラム次数を用いた場合のピーク周波数を高周波数成分とする。これ以外の例としては、対象データの直流成分を直流成分とし、対象データからこの直流成分を除いたデータをフーリエ変換して、低周波数成分領域（０．０４〜０．１５Ｈｚ）、高周波数成分領域（０．１５〜０．４０Ｈｚ）それぞれにおけるピーク周波数を各周波数成分とし、フーリエ変換における直流成分（０Ｈｚ）または対象データの平均値を直流成分とする方法がある。なお、本実施形態では、いずれの場合についても、低周波数成分、高周波数成分については、算出した周波数における位相成分を周波数領域で同時に算出し、補完関数に用いることを想定する。

また、本実施形態では、対象データの瞬時心拍の平均値、または、低周波数成分・高周波数成分の算出に用いた周波数スペクトルの直流成分を、直流成分とする。

このような、補完関数の各項について算出されたパラメータについて、算出した直流成分をＲＲＩ_ＤＣと表し、低周波数成分の振幅をＣ_ＬＦ、周波数をｆ_ＬＦ、位相をφ_ＬＦと表し、高周波数成分の振幅をＣ_ＨＦ、周波数をｆ_ＨＦ、位相をφ_ＨＦと表すとき、瞬時心拍補完部２８は、ｓｉｎ関数またはｃｏｓ関数を用いて、以下の式を補完関数として算出する。
estimatedＲＲＩ＝ＲＲＩ_ＤＣ＋Ｃ_ＬＦｓｉｎ（２πｆ_ＬＦｔ＋φ_ＬＦ）＋Ｃ_ＨＦｓｉｎ（２πｆ_ＨＦｔ＋φ_ＨＦ）
（または、estimatedＲＲＩ＝ＲＲＩ_ＤＣ＋Ｃ_ＬＦｃｏｓ（２πｆ_ＬＦｔ＋φ_ＬＦ）＋Ｃ_ＨＦｃｏｓ（２πｆ_ＨＦｔ＋φ_ＨＦ））

そして、瞬時心拍補完部２８は、このように瞬時心拍補完関数判別部２７で決定した補完関数の項に応じて補完関数の各項について算出されたパラメータを適用した補完関数によって、瞬時心拍評価部２５で瞬時心拍の異常値除外処理を行った或いは、瞬時心拍再計算部２６で再計算が行なわれた、少なくとも一つの欠損区間を有する瞬時心拍の時系列データの内、一つの未処理の欠損区間に対して、補完値を算出する（Ｓ６２）。この場合、補完値の算出は、未処理の欠損区間の内、時系列的に古いものから処理を行っても良いし、新しいものから処理を行っても良いし、任意のものに対して行っても構わない。

そして、瞬時心拍補完部２８は、この算出した補完値で欠損区間を補完（内挿）することとなる。なお、本実施形態では具体的な内挿間隔については特に規定しない。ＲＲＩタコグラムとして妥当と考えられる値のみを不等間隔で補完しても良いし、正常計測データのリサンプリングを行うサンプリングレートにしたがって等間隔に補完しても良い。なお、不等間隔で補完するときについては、任意時刻における補完関数の推定値が、直前のＲＲＩ観測時刻との経過時間と必ずしも合致しない時がある。このような場合について、例えば、これら二つの値が合致する実数を補完値としても良いし、これら二つの値の差が最も小さくなる値を補完値としても良い。

また特に、この欠損区間補完（内挿）処理において、瞬時心拍補完部２８が、ＲＲＩタコグラムとして妥当と考えられる値のみを不等間隔で補完する場合については、上記非特許文献２や上記非特許文献３に開示されているように、補完した瞬時心拍とその次の瞬時心拍との時間間隔が短くならないようにすることで、結果としてスプライン曲線が生理的に取り得ない値域に発振することを防止することが望ましい。そのため、算出した補完値で、直ちに欠損区間を補完（内挿）するのではなく、以下のような処理を実行する。

すなわち、瞬時心拍補完部２８は、欠損区間の時間を算出する（Ｓ６３）。

次に、欠損区間の時間が補完時間以上であるかが判断される（Ｓ６４）。Ｓ６４において、欠損区間の時間が補完時間以上であると判断された場合、瞬時心拍補完部２８は、欠損区間の時間から補完時間を差し引いた時間を算出する。そして、欠損区間の時間から補完時間を差し引いた時間（欠損区間の時間−補完時間）が補完対象時間以上であるかが判断される（Ｓ６５）。ここで、補完時間＜補完対象時間である。

Ｓ６５において、欠損区間の時間から補完時間を差し引いた時間（欠損区間の時間−補完時間）が補完対象時間以上であると判断された場合、瞬時心拍補完部２８は欠損区間を補完値で補完（内挿）し（Ｓ６６）、Ｓ６２の処理に戻る。なお、実施形態において、補完値の内挿方法は問わない。図７は、瞬時心拍電位の時系列データと内挿された補完値との関係を示す図である。同図においては、瞬時心拍の時系列データｄ１〜ｄ６に、算出された欠損区間ｄｐにおいて補完値ｄｃ_１が補完された例を示している。

なお、実施形態では、補完対象時間については詳細に規定しないが、例えば、２５０［ｍｓｅｃ］未満の値へのスプライン曲線の発振を防止するため、補完時間の１．２５倍とすることも考えられる。

こうして、Ｓ６２〜Ｓ６６の処理を繰り返すことで、欠損区間ｄｐにおいて補完値が内挿されていく。そして、例えば図８に示すように、３つめの補完値ｄｃ_３について、Ｓ６５において、欠損区間の時間から補完時間を差し引いた時間が、補完対象時間未満であると判定される。つまり、Ｓ６２で算出された補完値ｄｃ_３は次の瞬時心拍の時系列データｄ４との時間間隔が短いとして、この場合には、未処理の欠損区間があるか否かを確認する（Ｓ６７）。そして、未だ処理していない欠損区間がある場合には、瞬時心拍補完部２８は上記Ｓ６２の処理に戻る。補完関数の算出を都度行う場合は、Ｓ６１に戻っても良い。

また、Ｓ６４において、欠損区間の時間が補完時間未満の場合においても、瞬時心拍補完部２８は、上記Ｓ６７の処理に進んで、未処理の欠損区間があるか否かを確認し、未だ処理していない欠損区間がある場合には、Ｓ６２の処理に戻る。

こうして、Ｓ６２〜Ｓ６６の処理を繰り返し、全ての欠損区間に対する処理が終了したならば、瞬時心拍補完部２８は、このＳ６における補完処理を終了する。

なお、実施形態においては、図６に示したように、欠損区間を有する瞬時心拍の時系列データを補完する場合と、補完しない場合とがあるが、欠損区間を補完した瞬時心拍の時系列データに対しては、所定の識別子を付しても良い。従って、実施形態の瞬時心拍の時系列データの補完装置によれば、タコグラムの横軸を手法ｂによって生成し、データに識別子を付加することにより、欠損が生じていないデータと、欠損のあるデータとの区別が可能になる。

また、補完関数により補完値を求め、欠損区間をこの補完値で内挿することにより、スプライン曲線が生理的に取りえない値域に発振することを防止し、生理的に妥当なパワースペクトル密度関数を算出可能となる。

さらに、欠損区間の補完を行う前に、補完時間と補完対象時間の二つの時間を用いて欠損区間を評価することで、補完した瞬時心拍とその次の瞬時心拍との時間間隔が解析に影響を及ぼすほど短縮してしまうことを防止する。なお、欠損区間の評価を補完時間のみで行った場合、補完した瞬時心拍とその次の瞬時心拍との時間間隔が短くなり、結果としてスプライン曲線が生理的に取り得ない値域に発振する場合がある。

なお、Ｓ６１の補完関数の算出処理においては、直流成分、低周波数成分、高周波数成分のそれぞれについて、各成分の代表周波数のみを用いる場合について説明したが、各成分については複数の周波数成分を用いても良い。また、各成分において用いる周波数成分の数は一致しなくても良い。例えば、直流成分は一つ、低周波数成分は二つの周波数、高周波数成分は三つの周波数を組み合わせることとして、式の対応項数を増やしても良い。

また、Ｓ６１の補完関数の算出処理において、欠損区間の前の、正常に計測したデータが任意長であるとき、瞬時心拍補完部２８は、その正常計測データの変化傾向を踏まえて補完関数の位相をさらに調整しても良い。実施形態では具体的な手法については特に規定しないが、例えば、欠損区間に対応する補完関数が正常計測データの最終データに最も近い値からはじまるようにしたり、正常計測データの変曲点に合わせて補完関数の位相を調整したりすることなどが考えられる。

また、Ｓ６１の補完関数の算出処理において、直流成分、高周波数成分の二つの項を用いる場合は、低周波数成分を除く二つの成分について、Ｓ６１の補完関数の算出処理と同様にして補完関数を算出し、補完処理を行う。なお、この場合についても、上述のように、各成分において用いる周波数の数は一つ以上であれば良く、実施形態ではその数については規定しない。

また、Ｓ６１の補完関数の算出処理において、直流成分のみを用いる場合は、直流成分のみをＳ６１の補完関数の算出処理と同様にして補完関数を算出し、補完処理を行う。なお、この場合についても、各成分において用いる周波数の数は一つ以上であれば良く、実施形態ではその数については規定しない。

また、Ｓ６１の補完関数の算出処理において、瞬時心拍補完部２８は、複数の補完関数を用いても良い。例えば、ひとかたまりの正常計測データに対して一つの補完関数を算出し、ある欠損区間の補完には、そのデータ観測時刻に最も近い正常計測データから得た補完関数を用いることとしても良い。

以上のように、本発明の実施形態に係る瞬時心拍の時系列データの補完装置１２は、正常に計測した瞬時心拍の変動特性に基づいた補完関数によって瞬時心拍の欠損区間を補完することで、スプライン補間を用いた心拍変動のスペクトル解析を高精度に実現可能にすることができる。

すなわち、本発明の実施形態に係る瞬時心拍の時系列データの補完装置１２によれば、入力された瞬時心拍の周波数成分に基づいた補完関数によって欠損区間を補完することにより、スプライン曲線が生理的に取りえない値域に発振することを防止し、生理的に妥当なパワースペクトル密度関数を算出可能となる。

また、補完に用いる値を、瞬時心拍の特性に基づいて直流成分、低周波数成分、高周波数成分から算出した補完関数にしたがって算出することで、欠損が生じた場合でも、実態に近い周波数特徴量の算出が可能となる。このとき、低周波数帯域や高周波数帯域のスペクトルをすべて補完関数に反映するのではなく、代表周波数のみを用いることによって、ノイズやアーチファクトを誤判定した瞬時心拍のゆらぎの影響を低減することができる。

また、本発明の実施形態に係る瞬時心拍の時系列データの補完装置１２では、正常に計測したデータのデータ長に基づいて補完に用いる関数を変更することで、可能な限りもっともらしい値による欠損区間の補完を可能にすることができる。

以上、本発明の実施形態を説明したが、この実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。この新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。この実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。

・ 参考文献
（ｉ） 井上博，循環器疾患と自律神経機能，第２版，医学書院，２００１
（ｉｉ） 奥出潤，これならわかる！かんたんポイント心電図，第２版，医学書院，２０１１
（ｉｉｉ） 日本光電，雑音混入のメカニズムと対策：きれいな心電図を記録するポイント 〜ホルター心電図編〜，（２０１７年３月２３日確認），ｈｔｔｐ：／／ｗｗｗ．ｎｉｈｏｎｋｏｈｄｅｎ．ｃｏ．ｊｐ／ｉｒｙｏ／ｐｏｉｎｔ／ｈｏｌｔｅｒ／ｍｅｃｈａｎｉｓｍ．ｈｔｍｌ
（ｉｖ） 佐久間大輝ら，座位状態での心拍測定を用いたリアルタイムなストレス緩和システム，マルチメディア、分散協調とモバイルシンポジウム２０１３論文集，ｐｐ．１１８８−１１９５、２０１３
（ｖ） 横田康成ら，心拍変動時系列変化を用いた敗血症の前駆症状モニタリング，第５４回自動制御連合講演会，ｐｐ．１２５８−１２６１，２０１１

１０…瞬時心拍の補完システム、１１…心電計測部、１２…補完装置、２１…Ｒ波抽出部、２２…Ｒ波関連情報記録部、２３…瞬時心拍算出部、２４…瞬時心拍記録部、２５…瞬時心拍評価部、２６…瞬時心拍再計算部、２７…瞬時心拍補完関数判別部、２８…瞬時心拍補完部、２９…スペクトル解析処理部。

Claims (7)

1. 被験者の心電を計測する心電計測手段から出力される信号に基づいて算出される、時系列で隣接する２つのＲ波の間隔である瞬時心拍の内、正常な値と判断した複数の瞬時心拍からなる時系列データを入力とし、前記入力された前記時系列データの周波数成分に基づく補完関数によって補完値を算出し、前記算出した前記補完値により前記入力された前記時系列データにおいて欠損した瞬時心拍を補完する瞬時心拍補完部、
   を具備する、瞬時心拍の時系列データの補完装置。
2. 前記正常な値と判断した複数の瞬時心拍からなる前記時系列データのうち連続した瞬時心拍のデータ長に基づいて、前記時系列データにおいて欠損した前記瞬時心拍を補完するための前記補完関数の項を決定する瞬時心拍補完関数判別部を更に具備し、
   前記瞬時心拍補完部は、前記瞬時心拍補完関数判別部で決定した前記補完関数の前記項に応じて前記補完関数の各項について算出されたパラメータを適用した前記補完関数によって前記補完値を算出し、前記算出した前記補完値により前記入力された前記時系列データにおいて欠損した瞬時心拍を補完する、
   請求項１記載の瞬時心拍の時系列データの補完装置。
3. 前記瞬時心拍補完部における前記補完関数を構成する項は、直流成分、低周波数成分、高周波数成分のいずれか一つ以上である、請求項１または２記載の瞬時心拍の時系列データの補完装置。
4. 被験者の心電を計測する心電計測手段から出力される信号に基づいて算出される、時系列で隣接する２つのＲ波の間隔である瞬時心拍の内、正常な値と判断した複数の瞬時心拍からなる時系列データを入力し、
   前記入力した前記時系列データの周波数成分に基づく補完関数によって補完値を算出し、
   前記算出した前記補完値により、前記入力した前記時系列データにおいて欠損した瞬時心拍を補完する、
   瞬時心拍の時系列データの補完方法。
5. 前記正常な値と判断した複数の瞬時心拍からなる前記時系列データのうち連続した瞬時心拍のデータ長に基づいて、前記時系列データにおいて欠損した前記瞬時心拍を補完するための前記補完関数の項を決定し、
   前記瞬時心拍を補完することは、前記決定した前記補完関数の前記項に応じて前記補完関数の各項について算出されたパラメータを適用した前記補完関数によって、前記補完値を算出し、前記算出した前記補完値により、前記入力した前記時系列データにおいて欠損した瞬時心拍を補完する、
   請求項４記載の瞬時心拍の時系列データの補完方法。
6. 前記補完関数を構成する項は、直流成分、低周波数成分、高周波数成分のいずれか一つ以上である、請求項４または５記載の瞬時心拍の時系列データの補完方法。
7. コンピュータによって実行されたときに、前記コンピュータを、請求項１乃至請求項３いずれか１項に記載の前記瞬時心拍の時系列データの補完装置として機能させるためのプログラム。