(57)【要約】【課題】心拍のゆらぎを解析することによって、体が健康かどうか、逆に問題あるかどうかを、出来るだけ簡便に計測判断する。【解決手段】心拍数を測定し、心拍データを取得して、その時系列から平均値を計算し、時系列のそれぞれの要素とこの平均値との差をさらに差分し、新しい時系列を算出し、この新しい時系列を所定の長さｌのボックスに分割し、これをランダムに並べなおし、そのデータを足し合わせ、これからＤＦＡを用いて、スケーリング指数とする心拍ゆらぎの解析方法とし、スケーリング指数を基に、健康上問題があると判断する。【選択図】図５
2008173160.tif 000002
【特許請求の範囲】
【請求項１】
心拍数を測定し、心拍データを取得して、その時系列から平均値を計算し、
時系列のそれぞれの要素とこの平均値との差をさらに差分し、新しい時系列を算出し、
この新しい時系列を所定の長さｌのボックスに分割し、
これをランダムに並べなおし、そのデータを足し合わせ、
これからＤＦＡを用いて、スケーリング指数とする心拍ゆらぎの解析方法。
【請求項２】
前記所定の長さのボックスの長さｌは、３０〜１５０拍である請求項１の心拍ゆらぎの解析方法。
【請求項３】
前記並べなおしは、ランダムルーチン  メルセンヌ・ツイスター法によって行われる請求項１または２の心拍ゆらぎの解析方法。
【請求項４】
前記並べなおしは、時系列の長さｌのＮ／ｌ個のボックス（ここに、Ｎはｌの倍数の整数である。）に分割される請求項１〜３のいずれかの心拍ゆらぎの解析方法。
【請求項５】
前記ボックスが４〜１０００個の整数に分割される請求項４の心拍ゆらぎの解析方法。
【請求項６】
前記ボックスの長さｌが３０〜６０である請求項２〜６のいずれかの心拍ゆらぎの解析方法。
【請求項７】
前記スケーリング指数は、４次関数へのフィッティングで求める請求項１〜６のいずれかの心拍ゆらぎの解析方法。
【請求項８】
前記スケーリング指数は、ｌｏｇ−ｌｏｇプロットの指数である請求項１〜７のいずれかの心拍ゆらぎの解析方法。
【請求項９】
請求項１〜８のいずれかの解析方法で得られたスケーリング指数を基に、１．１超であったら、健康上問題があると判断する健康状態の判断方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、心拍のゆらぎの解析方法およびそれを用いた健康状態の判断方法に関するものである。
【背景技術】
【０００２】
従来、ゆらぎそのものは知られており、フーリエ解析（線形）での結果が病院や医学現場で論議され始めている。昨今の北米神経科学会議（参加者４万人、世界最大の神経の会議）でこの周波数解析が沢山発表されている。例えば、いろいろなゆらぎ成分のうちで、高周波成分（ＨＦ）は抑制性神経活動（副交感神経の心臓作用効果）を反映するものであり、低い周波数の成分（ＬＦ）は交感・副交感の両方を反映しており、自律神経の「バランス」を見るにはＨＦとＬＦの比を数値的に計算する、というような技術が、近年、具体的になってきている。
【０００３】
このゆらぎを診断しようという提案に関しては、いわば、我々に最も近い既存の技術として、２４時間、心電図を使い、ゆらぎ解析をして心臓病の患者と健常者とが違う指数になることを発見した報告がある。これはフラクタルとかマルチフラクタルなどとも呼ばれ、非線形テクニックに該当する。しかし、非線形解析法を実際に組み込んだ診断装置が完成しているわけではなく、また、診断方法概念として利用可能なものとして完成したか、あるいは器械や基盤やプログラムのような完成した形で発表したかどうか、と言うと、完成されたものとは言い難い。周期性のある生体信号を統計物理学的観点で数えなおし、統計生理学的演算により、体の状態と指数との関係に生理学的に意味を持たせたことが新規で優位であるといえる。
【０００４】
心拍のゆらぎが隠れた情報を含んでいる可能性があるという考え方は、１９８０年代に基礎理論が発表されているが、一昼夜等の長時間の心電データを適切に処理するには当時のコンピュータでは計算時間も膨大になり、適切な解析結果を得るまでには至っていなかったと推察される。そこで、我々は、昨今の現実に立ち返って、日常生活の中で、健康指針として簡単且つ短時間に計測評価できる方法を研究し、数理学的原理に立ち返り、理論を再考し、新たなアルゴリズムに反映した。その結果、心拍脈波の取得時間を短縮し、長時間の脈波データを必要としない、短時間の解析技術に仕上げることができた。この解析モデルで実際のヒトの健康状態との相関を積み上げてきたところ、短時間のデータでも、脈波データがあれば、そのゆらぎ方が、健康か否かを生理学的に判断できる確信を得た。
【０００５】
心臓のゆらぎで計算をしてみせた具体的論文（非特許文献１）や、心臓など生理学にゆらぎ解析が役立つことを主張し総説にまとめている論文  （非特許文献２）や、心拍変動  heart rate variability俗にＨＲＶというが、これを検査する最近主流のフーリエ解析法について沢山の論文から例示されている（非特許文献３）。
しかし、非特許文献１等の方法のＤＦＡ（Ｄｅｔｒｅｎｄｅｄ  Ｆｌｕｃｔｕａｔｉｏｎ  Ａｎａｌｉｓｉｓ：トレンド除去心拍分析）では、心拍ゆらぎの測定にあたり、数十分から２４時間の時間を要し、迅速且つ簡便・的確な判断ができないという問題点があった。
【非特許文献１】シー・ケイ・ペン，エス・ハブリン，エイチ・イー・スタンレイ，エー・エル・ゴールドバーガー(C.-K.Peng, S. Havlin, H. E. Stanley, and A. L. Goldberger), “クオンティフィケイション  オブ  スケーリング  エクスポーネンツ  アンド  クロスオーバー  フェノミナ  イン  ノンステショナリイ  ハートビート  タイム  シリーズ(Quantification of scaling exponents and crossover phenomenain nonstationary heartbeat time series)”,カオス( Chaos), 第５巻(Vol.5), 1995, pp.82-87
【非特許文献２】スタンレイ  エイチイー，等(Stanley HE, et al.) (1999)  スタティカル  フィジックス  アンド  フィジオロジイ：モノフラクタクル  アンド  マルチフラクタクル  アプローチイズ．フィジカ  エー(Statistical physics and physiology: Monofractal andmultifractal approaches. Physica A), 270:309-324.
【非特許文献３】ヴィラ，エムピー，等(Villa, MP, et al.) (2000)  エフェクツ  オブ  スリープ  ステージ  アンド  エージ  オン  ショート−ターム  ハート  レイト  ヴァリアビィリティ  ディアリング  スリープ  イン  ヘルシイ  インファンツ  アンド  チルドレン( Effects of sleep stage and age on short-term heart ratevariability during sleep in healthy infants and children).
【発明の開示】
【発明が解決しようとする課題】
【０００６】
本発明は、心拍のゆらぎを解析することによって、体が健康かどうか、逆に問題あるかどうかを、出来るだけ簡便に計測判断することを目的とする。さらに詳しくは、温度計で体温をはかったり、血圧計で血圧をはかるように簡便な方法で、例えば指先等の脈波から得られるゆらぎを解析指数化して、心臓の状態を把握することを目的とする。そこで、本発明者らは、健康指針として簡単且つ短時間に計測評価できる方法を研究し、数理学的原理に立ち返り、理論を再考し、新たなアルゴリズムに反映させた。その結果、心拍脈波の取得時間を短縮し、長時間の脈波データを必要としない、短時間の解析技術に仕上げることができた。この解析モデルで実際のヒトの健康状態との相関を積み上げてきたところ、短時間のデータでも、脈波データがあれば、そのゆらぎ方が、健康か否かを生理学的に判断できる確信を得た。
本発明者らは、ノコギリガザミに付いて、長期間の心拍を測定したところ、非特許文献１と同様に、指数が１．０を超えると、突然の死が訪れることを発見した。モデル動物でも非特許文献１と同じような結果が得られたのである。そこで、本発明者らは、ヒトに適用して、短時間のデータで健康状態の判断が出来る心拍ゆらぎの解析方法をについて研究したのである。
【課題を解決するための手段】
【０００７】
本発明は、心拍数を測定し、心拍データを取得して、その時系列から平均値を計算し、時系列のそれぞれの要素とこの平均値との差をさらに差分し、新しい時系列を算出し、この新しい時系列を所定の長さｌのボックスに分割し、これをランダムに並べなおし、そのデータを足し合わせ、これからＤＦＡを用いて、スケーリング指数とする心拍ゆらぎの解析方法である。
この場合、前記所定の長さのボックスの長さｌは、３０〜１５０拍であることが好ましい。
また、前記並べなおしは、ランダムルーチン  メルセンヌ・ツイスター法によって行われることが好ましい。
また、前記並べなおしは、時系列の長さｌのＮ／ｌ個のボックス（ここに、Ｎはｌの倍数の整数である。）に分割されることが好ましい。
そして、前記ボックスが４〜１０００個の整数に分割されることが好ましい。
さらに、前記ボックスの長さｌが３０〜６０であることが好ましい。
また、前記スケーリング指数は、４次関数へのフィッティングで求めることが好ましい。
前記スケーリング指数は、ｌｏｇ−ｌｏｇプロットの指数であることが好ましい。
さらに、本発明は、前記のいずれかの解析方法で得られたスケーリング指数を基に、１．１超であったら、健康上問題があると判断する健康状態の判断方法である。
【発明の効果】
【０００８】
本発明によれば、例えば、指先から例えば、２５０拍程度の心拍を取得するだけで、３００００拍以上にも及ぶデータに拡張できるデータが得られるので、簡便な方法で、ゆらぎ指数を解析指数化できるので、健康状態を迅速に判断できる。我々が得た成果は、２４時間の長い記録でなく、時間を短縮したことと、生理学意義付けを得た点が優位性である。
【発明を実施するための最良の形態】
【０００９】
本発明では、従来の心拍測定と同様に、先ず、心拍データの取得を行う。取得に当たっては、指先等に、電極を取り付け、心拍信号を取得すればよい。
図１には、取得された信号ｘｉが示される。
【００１０】
今、長さNの時系列ｘｋを考える。ｋは時系列の添え字で自然数とする（ｋ＝１，２，３，・・・・・，N）。この時系列は、得られた揺らぎのデータである。
【００１１】
最初にｘｋの平均値〈ｘ〉を計算する。そして時系列のそれぞれの要素と平均値との差を足し合わせた新しい時系列ｙｉを定義する（ステップ１）。
【００１２】
【数１】
2008173160.tif 000003
【００１３】
時系列ｙｉを長さｌのＮ／ｌ個のボックスに分割する。各々のボックスは重ならないように分割する（ステップ２）。
【００１４】
それぞれのボックスのトレンド（局所トレンド）ｙｖ（ｉ）を決定する。ここで、局所トレンドを最小自乗法によって決定された近似曲線によって定義する。最小自乗法で用いる関数は１次関数から４次関数を用いる。決定されたトレンドを差し引き、時系列ｚｉを定義する（ステップ３）。
ｚｉ＝ｙｉ−ｙｖ（ｉ）
【００１５】
図１（a）は、長さ５００の時系列ｘｉである。縦線はボックスの区切りを表す。ここではボックス幅は１００である。図１（ｂ）は、黒点は足し合わされた時系列ｙｉである。斜線はそれぞれの局所トレンドで４次関数へのフィッティングである。（図１（ｃ）は、局所トレンドを除いた時系列ｚｉである。
【００１６】
図５は、ｌogｌのｌogＳ（ｌ）に対するプロットである。
【００１７】
次に、ステップ４
【００１８】
【数２】
2008173160.tif 000004
【００１９】
を計算する。ここでｚｉｊ＋ｌ−ｚｉｊの値は長さｌのボックスにおける、１番目の要素ｚｉｊ＋１と最後の要素ｚｉｊ＋ｌとの差でｊ番目のボックスにおける変位に相当する（ステップ４）。
【００２０】
ボックス幅ｌの値を変え、ステップ２からステップ４の計算を繰り返し、それぞれのｌにおけるＳ（ｌ）の値を計算する。通常、ｌ＝１０〜Ｎ／４の範囲で行う（但し、Ｎ＞４０の場合）（ステップ５）。
【００２１】
Ｓ（ｌ）とｌの関係によりスケーリング指数を決定する。中心極限定理によって解釈すると、時系列の確率密度関数がガウス分布であった場合、それを足し合わせた時系列はランダムウォークを示しＳ（ｌ）∝ｌ０．５の関係を持つ。しかし、確率密度関数がガウス分布ではなく、分布の裾がベキ型の減衰を示す安定分布の場合Ｓ（ｌ）∝ｌαの関係を持つ。このαの値がスケーリング指数を示す。
ここで、α＝０．５の場合はガウス分布の場合と同様である。また、α＝１．０の場合、確率密度関数はコーシー分布となる（ステップ６）。
【００２２】
このような場合、図２〜４には、より好ましい態様が示される。図２は、ボックス数５０で行っている。そして、１〜５のボックスをメルセンヌ・ツイスター法のランダムルーチンによって、並べなおしをしている（図３）。さらに、ボックスの継ぎ目でズレが生じているので、これを解決するために、データを差分する。即ち、ボックスの平均値を０にして、ランダムに並べなおす。そして、データを足していく（図４）。
以上の方法により２５０拍データを３００００拍に拡張しＤＦＡを行えばよい。
【００２３】
ＤＦＡの方法を図を使って説明すると、図５のようになる。図５は、ｌogｌの図であり、この図の傾きが即ちＳ（ｌ）∝ｌαで決定されるスケーリング指数である。局所トレンドを決定する際に用いる関数の次数によってスケーリング指数は変化するが、次数を上げていくと殆ど変化しなくなる。これのようになって初めてトレンドが除かれたと考えられ、その次数でのスケーリング指数を以って時系列ｘｉのスケーリング指数として採用する。図５の場合、４次関数のフィッティングに最も近似しているので、スケーリング指数αは０．７５と決定した。
【００２４】
ペンら（Peng et aｌ.）（１９９４）のＤＦＡの方法と我々のＤＦＡの方法で異なる点はステップ４における計算のみである。PengらはＳ（ｌ）の計算のかわりに分散それぞれのボックスにおいてＦｊ（ｌ）を計算している。
【００２５】
【数３】
2008173160.tif 000005
【００２６】
但し、ｙｖ（ｘ）はディトレンドで使用したフィッティング関数である。
そして各ボックスでの分散を平均した
【００２７】
【数４】
2008173160.tif 000006
【００２８】
の値をもってスケーリング指数を決定している。スケーリング指数がＦ（ｌ）∝ｌαのようにｌに依存すると言う。
【００２９】
次に、このような計算の流れをフロ−チャートによって説明する（図６、図７）。
図６では、先ず、心拍データを取得する。次いで、心拍数に変換した後、差分処理を行う。
そして、小区間のデータ数ｎを指定する。ｎは、２０より大であって、通常２０〜１５０、好ましくは、３０〜８０、さらには、３０〜６０、特に５０が好ましい。
その際、各区間番号を付ける。ボックスの数は、４〜１０００が好ましい。
そして、各小区間の平均値を０に変換する。
この後、データ数が３００００を超えているか否かを判断する。
３００００未満のときには、ランダムルーチンによる番号を決定する。ランダムルチーンとしては、公知のメルセンヌ・ツイスター法を用いる。
当該番号の小区間データを追加する。その後、再びデータ数が３００００未満か否かを判断する。そして、３００００超であったら、加算処理に移る。
ＤＦＡでは、図７に示されるように、先ず、データの読み込みを行う。その後、ＤＦＡ、Ｂｏｘサイズｌの初期値をセットする。ｌが例えば、１０００未満であるか否かを判断する。１０００未満であれば、データを長さｌの小区間に分割する。各区間にて、最小二乗法によるディトレンドを行う。そして、各区間の最初と最後のデータの差の二乗平均の平方根ｆを記録する。そして、再度、Ｂｏｘサイズの指定を行う。これにより、１０００を超えたときには、フィッティングＢｏｘサイズ区間を指定する。そして、Ｂｏｘサイズｌと各二乗平均平方根とをｌｏｇ−ｌｏｇスケールでフィッティングする。
フィッティングの一次の係数によりスケーリング指数が決定される。
【実施例】
【００３０】
本発明の発明者の一人、矢澤の心拍の測定結果を図８に示す。
心拍ゆらぎの算出には、図２〜４に示されるボックス数５０とこのボックスを１〜５個並べなおして行った。この結果、基準の０．７８に極めて近い０．６６が得られた。この場合は、注意（突然死の危険はないが、低すぎる。制御系に不健康な問題が内包されている。）との判断がなされる。
その他、同様にして、何人もの被験者のゆらぎを算出した結果を図９〜１２に示す。
これらの例から、本発明による算出方法が優れた一致を見せることがわかる。
なお、図１１に示される被験者は、スケーリング指数αが１．４を超えており、その後死亡した例である。
このようにして、スケーリング指数を算出することができるが、この指数を用いて、健康状態を判断することができる。１に近い値なら、心配無用のシグナルを出すようにすればよい。これに対し、低い数値なら、ストレスや運動やサプリメントを要するので、注意（突然死の危険はないが、低すぎる。制御系に不健康な問題が内包されている。）のシグナルを発すればよい。異常に高い数値のときは、心臓疾患の可能性があるので、警告を発する。このように、健康を指向する利用者には好都合である。また、医療現場では、新しい診断項目の追加ができるだけでなく、末期医療における生体システムの動作状態を数値化できることも可能となる。その他、医学領域のみならず、ペット、競馬馬、動物園等の獣医学等関係領域での診断、運輸交通機関運転士（航空機、電車、トラック、バス、タクシーなど）の始業開始前、運転中診断、警備保障会社の遠隔地診断や生命保険会社での診断技術としての利用、携帯電話端末等に代表される移動端末への心拍監視診断警告発信機構の組み込み、ゲーム、クイズ等の遊戯器具や、心電計、血圧計、便座等の既存の脈波検出可能な装置への接続、この他、パソコンのマウス、自動車、椅子、ベット、眼鏡、腕輪、腕時計、入れ歯、健康器具のハンドル、イヤホン等のおよそ身体に接触して数分間の脈波の検出に供せられる状況にあるものへの接続などによる診断に適用することができ、高齢化社会対策、高密度交通往来危険対策、深刻な病気になる前の予防予知の支援が可能となる。
【産業上の利用可能性】
【００３１】
本発明によれば、簡易に且つ迅速に健康上の問題あるかどうかの判断をすることができる。
【図面の簡単な説明】
【００３２】
【図１】本発明の算出方法を説明するためのグラフである。
【図２】本発明の算出方法によって得られるＤＦＡのグラフである。
【図３】本発明の算出方法によって得られるＤＦＡのグラフである。
【図４】本発明の算出方法によって得られるＤＦＡのグラフである。
【図５】実施例の算出方法によって得られたＤＦＡのグラフである。
【図６】本発明の算出方法による計算の流れを示すフロ−チャートである。
【図７】本発明の算出方法による計算の流れを示すフロ−チャートである。
【図８】実施例の算出方法によって得られたＤＦＡのグラフである。
【図９】実施例の算出方法によって得られたＤＦＡのグラフである。
【図１０】実施例の算出方法によって得られたＤＦＡのグラフである。
【図１１】実施例の算出方法によって得られたＤＦＡのグラフである。
【図１２】実施例の算出方法によって得られたＤＦＡのグラフである。
【図１】
2008173160.tif 000007
【図２】
2008173160.tif 000008
【図３】
2008173160.tif 000009
【図４】
2008173160.tif 000010
【図５】
2008173160.tif 000011
【図６】
2008173160.tif 000012
【図７】
2008173160.tif 000013
【図８】
2008173160.tif 000014
【図９】
2008173160.tif 000015
【図１０】
2008173160.tif 000016
【図１１】
2008173160.tif 000017
【図１２】
2008173160.tif 000018