特開 2022-65 血管剛性推定方法、血管剛性推定装置及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】特開2022-65(P2022-65A)

(43)【公開日】2022年1月4日

(54)【発明の名称】血管剛性推定方法、血管剛性推定装置及びプログラム

(51)【国際特許分類】

   A61B 5/02 20060101AFI20211203BHJP

【ＦＩ】

   A61B5/02 A

【審査請求】未請求

【請求項の数】9

【出願形態】ＯＬ

【全頁数】13

(21)【出願番号】特願2020-105880(P2020-105880)

(22)【出願日】2020年6月19日

【国等の委託研究の成果に係る記載事項】（出願人による申告）令和元年度、国立研究開発法人科学技術振興機構研究成果展開事業センター・オブ・イノベーションプログラム「精神的価値が成長する感性イノベーション拠点」委託研究開発、産業技術力強化法第１７条の適用を受ける特許出願

(71)【出願人】

【識別番号】504136568

【氏名又は名称】国立大学法人広島大学

(74)【代理人】

【識別番号】100196380

【弁理士】

【氏名又は名称】森 匡輝

(72)【発明者】

【氏名】辻 敏夫

(72)【発明者】

【氏名】古居 彬

(72)【発明者】

【氏名】曽 智

(72)【発明者】

【氏名】棟安 俊文

(72)【発明者】

【氏名】笹岡 貴史

(72)【発明者】

【氏名】山脇 成人

【テーマコード（参考）】

4C017

【Ｆターム（参考）】

4C017AA07

4C017AA08

4C017AA09

4C017AA19

4C017AB03

4C017BC11

4C017FF05

(57)【要約】

【課題】心拍より高いサンプリング周波数で血管剛性を推定可能な血管剛性推定方法、血管剛性推定装置及びプログラムを提供する。
【解決手段】本発明に係る血管剛性推定方法は、容積脈波と血圧との関係で表され、容積脈波の二階微分を含む慣性項、容積脈波の一階微分を含む粘性項及び容積脈波の指数関数を含む剛性項を含む血管剛性モデルに基づいて、容積脈波の二階微分の係数である慣性パラメータ、容積脈波の一階微分の係数である粘性パラメータ、剛性項を線形近似して得られる容積脈波の係数である剛性パラメータ及び定数項を最小二乗法によって推定する、第１推定工程と、慣性パラメータ及び粘性パラメータを血管剛性モデルに代入して、最小二乗法を用いることにより、血管剛性を推定する、第２推定工程と、を含む。
【選択図】図２

【特許請求の範囲】

【請求項1】

容積脈波と血圧との関係で表され、前記容積脈波の二階微分を含む慣性項、前記容積脈波の一階微分を含む粘性項及び前記容積脈波の指数関数を含む剛性項を含む血管剛性モデルに基づいて、前記容積脈波の二階微分の係数である慣性パラメータ、前記容積脈波の一階微分の係数である粘性パラメータ、前記剛性項を線形近似して得られる前記容積脈波の係数である剛性パラメータ及び定数項を最小二乗法によって推定する、第１推定工程と、
前記慣性パラメータ及び前記粘性パラメータを前記血管剛性モデルに代入して、最小二乗法を用いることにより、血管剛性を推定する、第２推定工程と、を含む、
ことを特徴とする血管剛性推定方法。

【請求項2】

前記血管剛性モデルは、以下に示す対数線形化末梢血管粘弾性モデル

【数1】

ただし、μは血管壁の慣性、ηは血管壁の粘性、βは血管壁の剛性、Ｐ_ｂ（ｔ）は時刻ｔにおける動脈血圧、Ｐ_ｌ（ｔ）は時刻ｔにおける容積脈波、Ｐ_ｂβ０は基準血圧、Ｐ_ｂβｎｌ（Ｐ_ｌ（ｔ））は静脈影響
である、
ことを特徴とする請求項１に記載の血管剛性推定方法。

【請求項3】

前記血管剛性は、心拍より高いサンプリング周波数で推定される、
ことを特徴とする請求項１又は２に記載の血管剛性推定方法。

【請求項4】

前記第２推定工程では、
平均血圧以上の範囲で最小二乗法を適用することにより、
前記血管剛性を推定する、
ことを特徴とする請求項１から３のいずれか一項に記載の血管剛性推定方法。

【請求項5】

前記第１推定工程及び前記第２推定工程の最小二乗法の窓幅は、
心拍間隔よりも長い、
ことを特徴とする請求項１から４のいずれか一項に記載の血管剛性推定方法。

【請求項6】

前記血圧と前記容積脈波との時間差を補正する時間差補正工程を含む、
ことを特徴とする請求項１から５のいずれか一項に記載の血管剛性推定方法。

【請求項7】

前記時間差補正工程では、
心拍ごとに、対応する個別時間差を算出し、前記第１推定工程の最小二乗法の窓幅に対する各心拍の影響度合いに基づいた重みを前記個別時間差に乗じて得られる調整後個別時間差に基づいて、前記時間差を算出する、
ことを特徴とする請求項６に記載の血管剛性推定方法。

【請求項8】

容積脈波と血圧との関係で表され、前記容積脈波の二階微分を含む慣性項、前記容積脈波の一階微分を含む粘性項及び前記容積脈波の指数関数を含む剛性項を含む血管剛性モデルに基づいて、血管剛性を推定する演算部を備え、
前記演算部は、
前記血管剛性モデルのうち、前記容積脈波の二階微分の係数である慣性パラメータ、前記容積脈波の一階微分の係数である粘性パラメータ、前記剛性項を線形近似して得られる前記容積脈波の係数である剛性パラメータ及び定数項を最小二乗法によって推定し、
推定された前記慣性パラメータ及び前記粘性パラメータを前記血管剛性モデルに代入して、最小二乗法を用いることにより、前記血管剛性を推定する、
ことを特徴とする血管剛性推定装置。

【請求項9】

コンピュータを、
容積脈波と血圧との関係で表され、前記容積脈波の二階微分を含む慣性項、前記容積脈波の一階微分を含む粘性項及び前記容積脈波の指数関数を含む剛性項を含む血管剛性モデルに基づいて、
前記容積脈波の二階微分の係数である慣性パラメータ、前記容積脈波の一階微分の係数である粘性パラメータ、前記剛性項を線形近似して得られる前記容積脈波の係数である剛性パラメータ及び定数項を最小二乗法によって推定し、
前記慣性パラメータ及び前記粘性パラメータを前記血管剛性モデルに代入して、最小二乗法を用いることにより、血管剛性を推定する演算部、
として機能させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、血管剛性推定方法、血管剛性推定装置及びプログラムに関する。

【背景技術】

【0002】

従来、血管剛性の変化を指標として、自律神経活動などの身体の状態を把握する方法が開発されている（例えば、非特許文献１）。非特許文献１では、血管剛性と筋交感神経活動との相関関係に着目し、血管剛性から自律神経活動の働きを評価する方法が開示されている。

【先行技術文献】

【非特許文献】

【0003】

【非特許文献1】平野博大，他９名、「対数線形化末梢血管粘弾性インデックスの提案と交感神経遮断術への応用」、計測自動制御学会論文集、Vol.48，No.11，p.731-739、2012年

【発明の概要】

【発明が解決しようとする課題】

【0004】

非特許文献１では、血管剛性として、血圧と指尖容積脈波との関係に基づく対数線形化末梢血管粘弾性モデルを用いて推定される末梢の血管剛性を用いることとしている。また、非特許文献１の血管剛性推定方法では、指数関数となる指尖容積脈波を含む剛性項を消去して血管剛性を算出するため、心拍を基準とした差分計算を行っているので、心拍１拍ごとに血管剛性が推定される。

【0005】

しかしながら、身体の状態を精度よく把握するためには、血管剛性をより細かく、すなわち高いサンプリング周波数で推定することが求められる。

【0006】

本発明は、上述の事情に鑑みてなされたものであり、心拍より高いサンプリング周波数で血管剛性を推定可能な血管剛性推定方法、血管剛性推定装置及びプログラムを提供することを目的とする。

【課題を解決するための手段】

【0007】

上記目的を達成するために、この発明の第１の観点に係る血管剛性推定方法は、
容積脈波と血圧との関係で表され、前記容積脈波の二階微分を含む慣性項、前記容積脈波の一階微分を含む粘性項及び前記容積脈波の指数関数を含む剛性項を含む血管剛性モデルに基づいて、前記容積脈波の二階微分の係数である慣性パラメータ、前記容積脈波の一階微分の係数である粘性パラメータ、前記剛性項を線形近似して得られる前記容積脈波の係数である剛性パラメータ及び定数項を最小二乗法によって推定する、第１推定工程と、
前記慣性パラメータ及び前記粘性パラメータを前記血管剛性モデルに代入して、最小二乗法を用いることにより、血管剛性を推定する、第２推定工程と、を含む。

【0008】

また、前記血管剛性モデルは、以下に示す対数線形化末梢血管粘弾性モデル

【数1】

ただし、μは血管壁の慣性、ηは血管壁の粘性、βは血管壁の剛性、Ｐ_ｂ（ｔ）は時刻ｔにおける動脈血圧、Ｐ_ｌ（ｔ）は時刻ｔにおける容積脈波、Ｐ_ｂβ０は基準血圧、Ｐ_ｂβｎｌ（Ｐ_ｌ（ｔ））は静脈影響
である、
こととしてもよい。

【0009】

また、前記血管剛性は、心拍より高いサンプリング周波数で推定される、
こととしてもよい。

【0010】

また、前記第２推定工程では、
平均血圧以上の範囲で最小二乗法を適用することにより、
前記血管剛性を推定する、
こととしてもよい。

【0011】

また、前記第１推定工程及び前記第２推定工程の最小二乗法の窓幅は、
心拍間隔よりも長い、
こととしてもよい。

【0012】

また、前記血圧と前記容積脈波との時間差を補正する時間差補正工程を含む、
こととしてもよい。

【0013】

また、前記時間差補正工程では、
心拍ごとに、対応する個別時間差を算出し、前記第１推定工程の最小二乗法の窓幅に対する各心拍の影響度合いに基づいた重みを前記個別時間差に乗じて得られる調整後個別時間差に基づいて、前記時間差を算出する、
こととしてもよい。

【0014】

また、本発明の第２の観点に係る血管剛性推定装置は、
容積脈波と血圧との関係で表され、前記容積脈波の二階微分を含む慣性項、前記容積脈波の一階微分を含む粘性項及び前記容積脈波の指数関数を含む剛性項を含む血管剛性モデルに基づいて、血管剛性を推定する演算部を備え、
前記演算部は、
前記血管剛性モデルのうち、前記容積脈波の二階微分の係数である慣性パラメータ、前記容積脈波の一階微分の係数である粘性パラメータ、前記剛性項を線形近似して得られる前記容積脈波の係数である剛性パラメータ及び定数項を最小二乗法によって推定し、
推定された前記慣性パラメータ及び前記粘性パラメータを前記血管剛性モデルに代入して、最小二乗法を用いることにより、前記血管剛性を推定する。

【0015】

また、本発明の第３の観点に係るプログラムは、
コンピュータを、
容積脈波と血圧との関係で表され、前記容積脈波の二階微分を含む慣性項、前記容積脈波の一階微分を含む粘性項及び前記容積脈波の指数関数を含む剛性項を含む血管剛性モデルに基づいて、
前記容積脈波の二階微分の係数である慣性パラメータ、前記容積脈波の一階微分の係数である粘性パラメータ、前記剛性項を線形近似して得られる前記容積脈波の係数である剛性パラメータ及び定数項を最小二乗法によって推定し、
前記慣性パラメータ及び前記粘性パラメータを前記血管剛性モデルに代入して、最小二乗法を用いることにより、血管剛性を推定する演算部、
として機能させる。

【発明の効果】

【0016】

本発明の血管剛性推定方法、血管剛性推定装置およびプログラムによれば、サンプリング時刻ごとに、最小二乗法を用いて慣性パラメータ及び粘性パラメータを推定し、推定された慣性パラメータ及び粘性パラメータを用いて血管剛性を推定する。したがって、心拍より高いサンプリング周波数で血管剛性を推定することが可能である。

【図面の簡単な説明】

【0017】

【図1】本発明の実施の形態に係る血管剛性推定装置の機能ブロック図である。

【図2】実施の形態１に係る血管剛性推定の流れを示すフローチャートである。

【図3】実施の形態２に係る血管剛性推定の流れを示すフローチャートである。

【図4】血圧と容積脈波との時間差を示す概念図である。

【図5】数値例に係る計測波形を示すグラフである。

【図6】数値例に係るＰＳＤを示すグラフであり、（Ａ）は第１解析区間のグラフ、（Ｂ）は第２解析区間のグラフである。

【発明を実施するための形態】

【0018】

以下、図を参照しつつ、本発明の実施の形態に係る血管剛性推定装置１について説明する。図１のブロック図に示すように、血管剛性推定装置１は、センサ１０、推定装置本体２０を備える。

【0019】

センサ１０は、取得する生体信号の種類に合わせて選択される検出器である。取得する生体信号は、血管剛性βを算出するために用いられる時系列信号であり、例えば、心臓の活動を示す心電、血圧、血管の容積変化を示す容積脈波等である。

【0020】

推定装置本体２０は、例えばコンピュータ装置であり、制御部２１、記憶部２２、表示部２３、入力部２４を備える。

【0021】

制御部２１は、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、水晶発振器等から構成されており、血管剛性推定装置１の動作を制御する。また、制御部２１は、センサ１０で取得された生体信号と、予め記憶部２２に記憶されている血管剛性βの算出モデル（以下、血管剛性モデルという。）とに基づいて、血管剛性βを算出する。制御部２１は、制御部２１のＲＯＭ、記憶部２２等に記憶されている各種動作プログラム及びデータをＲＡＭに読み込んでＣＰＵを動作させることにより、図１に示す制御部２１の各機能を実現させる。これにより、制御部２１は、生体信号取得部２１１、演算部２１２として動作する。

【0022】

生体信号取得部２１１は、センサ１０を制御して、被験者の生体信号を取得する。また、取得した生体信号を記憶部２２へ送信し、記憶させる。本実施の形態では、複数のセンサ１０を用いて、同時に複数の生体信号、すなわち心電、血圧、容積脈波等が取得される。

【0023】

演算部２１２は、生体信号取得部２１１でセンサ１０から取得した生体信号と血管剛性モデルとに基づいて、血管剛性βを算出する。また、演算部２１２は、算出された血管剛性βを推定結果として記憶部２２に記憶させる。血管剛性βの推定方法の詳細については、後述する。

【0024】

記憶部２２は、ハードディスク、フラッシュメモリ等の不揮発性メモリであり、血管剛性モデル、演算部２１２で算出された血管剛性βの推定結果等を記憶する。

【0025】

表示部２３は、コンピュータ装置である推定装置本体２０に備えられた表示用デバイスであり、例えば液晶パネルである。表示部２３は、センサ１０で取得された生体信号、演算部２１２で推定された血管剛性β等を表示する。

【0026】

入力部２４は、生体信号取得の開始、終了指示、血管剛性βを推定するためのパラメータ等を入力するための入力デバイスである。入力部２４は、推定装置本体２０に備えられたキーボード、マウス等である。

【0027】

続いて、血管剛性推定装置１を用いた血管剛性推定方法について、図２のフローチャートを参照しつつ、具体的に説明する。本実施の形態では、心電、動脈血圧及び指尖容積脈波を生体信号として取得して、生体信号と血管剛性算出モデルとに基づいて血管剛性βを推定する場合を例として説明する。

【0028】

まず、測定工程として、医師等の測定者は、心電、動脈血圧及び指尖容積脈波を測定するためのセンサ１０を、被験者に装着する。そして、測定者は、入力部２４を操作して、心電、動脈血圧及び指尖容積脈波の取得を開始する。生体信号取得部２１１は、生体信号の取得開始指示により、センサ１０を制御して、被験者の心電、動脈血圧及び指尖容積脈波を取得する（ステップＳ１１）。生体信号取得部２１１は、取得した生体信号のデータを記憶部２２へ送信し、記憶部２２に記憶させる（ステップＳ１２）。

【0029】

続いて、第１推定工程として、演算部２１２は、生体信号取得部２１１で取得した生体信号のデータと、記憶部２２に予め記憶されている血管剛性βの算出モデルに基づいて、パラメータ推定を行う（ステップ１３）。

【0030】

ここで、第１推定工程に係るパラメータ推定について詳細に説明する。本実施の形態に係る血管剛性βの算出モデルは、以下の式（１）に示す対数線形化末梢血管粘弾性モデルである。

【数2】

ここで、μは血管壁の慣性、ηは血管壁の粘性、βは血管壁の剛性（血管剛性）を表す。また、Ｐ_ｂ（ｔ）は時刻ｔにおける動脈血圧、Ｐ_ｌ（ｔ）は時刻ｔにおける血管の容積変化（指尖容積脈波）、Ｐ_ｂβ０は基準血圧、Ｐ_ｂβｎｌ（Ｐ_ｌ（ｔ））は静脈影響を表す。

【0031】

静脈影響項Ｐ_ｂβｎｌ（Ｐ_ｌ（ｔ））は、被験者の動脈で計測される動脈血圧、指尖容積脈波に対する影響を考慮するためのパラメータであり、例えば、高血圧、うっ血性心不全などの心臓の影響、駆血、重力などの静脈の影響等が考慮される。静脈影響Ｐ_ｂβｎｌ（Ｐ_ｌ（ｔ））は、被験者及び計測状況に応じて適宜調整される。

【0032】

式（１）の血管剛性モデルは、容積脈波と血圧との関係で表され、容積脈波の二階微分を含む慣性項、容積脈波の一階微分を含む粘性項及び容積脈波の指数関数を含む剛性項を含む。

【0033】

式（１）の指数項について、Ｐ_ｌ（ｔ）＝０まわりでマクローリン展開を行い、線形近似すると、以下の式（２）となる。

【数3】

ただし、β_Ａは以下の式のとおりである。

【数4】

【0034】

慣性μ、粘性η、exp｛Ｐ_ｂβｎｌ＋Ｐ_ｂβｎｌ（０）｝、剛性パラメータβ_Ａは、任意の窓幅ｗの血圧Ｐ_ｂ（ｔ）、容積脈波Ｐ_ｌ（ｔ）を用いて、式（２）に基づき線形最小二乗法を実行することによって推定することができる。

【0035】

本実施の形態に係る演算部２１２は、窓幅ｗとして、センサ１０で測定された動脈血圧Ｐ_ｂ（ｔ）、容積脈波Ｐ_ｌ（ｔ）を用い、上記で推定した容積脈波の二階微分の係数である慣性パラメータのμ、容積脈波の一階微分の係数である粘性パラメータのη、剛性項を線形近似して得られる容積脈波の係数である剛性パラメータのβ_Ａ、定数項であるexp｛Ｐ_ｂβｎｌ＋Ｐ_ｂβｎｌ（０）｝を最小二乗法によって推定する。窓幅ｗ及びサンプリング周波数は、特に限定されないが、上記複数のパラメータを最小二乗法で推定するのに適した値に設定することが好ましい。

【0036】

これにより、血管剛性推定装置１は、任意の窓幅ｗでパラメータ推定を行うことができ、心拍より高いサンプリング周波数でパラメータ推定を行うことが可能である。すなわち、以下の式に示す従来のパラメータ推定方法のように、心拍（Ｒ波）でタイミングを調整し、異なる時刻の差分から、指数項を消去してパラメータ推定を行う必要がないので、本発明の方法では、より高いサンプリング周波数でパラメータ推定を行うことができる。

【数5】

【0037】

図２のフローチャートに戻り、第２推定工程として、演算部２１２は、第１推定工程で推定された慣性μ、粘性ηを用いて、血管剛性βを推定する（ステップＳ１４）。具体的な推定方法は以下の通りである。

【0038】

式（１）を指数項とそれ以外の項で整理し、両辺に自然対数をとると、以下の式（３）となる。

【数6】

【0039】

ここで、右辺の各項は、第１推定工程で推定された慣性μ、粘性ηと測定データから数値として求めることができる。演算部２１２は、式（３）に慣性μ、粘性ηを代入して、最小二乗法を適用することにより血管剛性βを推定する。

【0040】

また、最小二乗法は、静脈影響を受けない範囲で行う。具体的には、血圧Ｐ_ｂ（ｔ）が、平均血圧以上となる範囲で最小二乗法による推定を行う。これにより、式（３）の左辺第３項を０として、血管剛性β、基準血圧Ｐ_ｂβ０を推定することができる。そして、演算部２１２は、推定された血管剛性βと対応するサンプリング時刻とを記憶部２２に記憶させる。

【0041】

演算部２１２は、血管剛性推定処理が、予め定められた終了条件に達するまで（ステップＳ１５のＮＯ）、サンプリングステップを進めて窓幅ｗをずらしながら（ステップＳ１６）、ステップＳ１３からＳ１４の血管剛性βの推定処理を繰り返す。

【0042】

また、制御部２１は、血管剛性推定処理が、予め定められた所定の終了条件に達した場合（ステップＳ１５のＹＥＳ）、血管剛性βの推定を終了する。そして、制御部２１は、推定結果としての血管剛性βを表示部２３に表示させるとともに、記憶部２２に記憶させる（ステップＳ１７）。

【0043】

以上説明したように、本実施の形態に係る血管剛性推定方法では、第１推定工程において、心拍数に拘束されることなく、最小二乗法によって各項のパラメータを推定する。また、第２推定工程として、第１推定工程で推定された慣性μ、粘性ηを用いて血管剛性βを推定するので、心拍より高いサンプリング周波数で血管剛性を推定することが可能である。

【0044】

（実施の形態２）
実施の形態１では、式（１）の血管剛性モデルで用いられる血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との時刻の差は考慮せず、計測された時刻ｔにおける血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）とを用いて血管剛性βを推定することとした。一般的に、血圧Ｐ_ｂ（ｔ）の測定箇所と容積脈波Ｐ_ｌ（ｔ）の測定箇所とは同一の箇所ではなく、血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）とは、被験者の体に別個の測定器具を取り付けて測定される。

【0045】

したがって、血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との間には血流による変化に時間的な誤差が生じる。血管剛性βの推定間隔を小さくし、サンプリング周波数を高めると、血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との間の時間的な誤差の影響が大きくなると考えられる。よって、精度よく血管剛性βを推定するためには、血圧Ｐ_ｂ（ｔ）の変化と容積脈波Ｐ_ｌ（ｔ）の変化との間の時間差を補正することが求められる。本実施の形態では、血管剛性βの推定を行う前に、血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との間の時間差を補正する場合について説明する。

【0046】

本実施の形態に係る血管剛性推定装置１では、演算部２１２が、測定された血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）とに基づいて、血管剛性βを推定する前に、時間差補正工程として、血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との時間差を補正する点で、実施の形態１と異なる。その他、血管剛性推定装置１に係るセンサ１０、推定装置本体２０の記憶部２２、表示部２３等の構成については、実施の形態１と同様であるので同じ符号を付して、詳細な説明は省略する。

【0047】

図３のフローチャートを参照しつつ、本実施の形態に係る血管剛性推定方法について具体的に説明する。ステップＳ３１〜Ｓ３２に係る測定工程は、実施の形態１に係るステップＳ１１〜Ｓ１２の計測工程と同様である。続いて、演算部２１２は、時間差補正工程として、記憶部２２に記憶されている測定データに基づいて、血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との時間差を補正する。時間差補正に用いる測定データは、例えば、図４に示すように、心電、血圧Ｐ_ｂ（ｔ）、容積脈波Ｐ_ｌ（ｔ）とする。

【0048】

本実施の形態に係る窓幅ｗ、すなわち第１推定工程、第２推定工程の最小二乗法で用いられる窓時間は、１．２秒としている。したがって、通常、窓幅ｗは、心拍間隔よりも長いので、窓幅ｗ内には、心電データから得られる心拍、すなわちＲ波は、１つ又は２つ存在する。演算部２１２は、時間差補正工程として、血管剛性推定を行う対象の窓幅ｗ内にＲ波がいくつあるかを示す心拍回数Ｐ_ｋを算出する（ステップＳ３３）。

【0049】

続いて、演算部２１２は、窓幅ｗにおいて、窓の開始時刻から、窓幅ｗ内の最初のＲ波ピークまでの時間である重みｗ_１を算出する（ステップＳ３４）。

【0050】

心拍回数Ｐ_ｋが１である場合（ステップＳ３５のＹＥＳ）、演算部２１２は、最初のＲ波ピークから窓の終了時刻までの時間である重みｗ_２を算出するとともに、重みｗ_３＝０と設定する（ステップＳ３６）。

【0051】

心拍回数Ｐ_ｋが２である場合（ステップＳ３５のＮＯ）、演算部２１２は、最初のＲ波ピークから、２つ目のＲ波ピークまでの時間である重みｗ_２を算出する。さらに、演算部２１２は、２つ目のＲ波ピークから窓の終了時刻までの時間である重みｗ_３を算出する（ステップＳ３７）。

【0052】

続いて、演算部２１２は、血圧Ｐ_ｂ（ｔ）の変曲点の時刻を算出する。変曲点は、血圧Ｐ_ｂ（ｔ）の１階微分の値が極大となる時刻である。また、演算部２１２は、変曲点の時刻として、窓の開始時刻以前で、窓の開始時刻から最も近い変曲点の時刻（φｂ_１）、窓の開始時刻以降で窓の開始時刻から最も近い変曲点の時刻（φｂ_２）、及びφｂ_２の次の変曲点の時刻（φｂ_３）の３つを算出する（ステップＳ３８）。

【0053】

また、演算部２１２は、容積脈波Ｐ_ｌ（ｔ）の変曲点の時刻を算出する。変曲点は、容積脈波Ｐ_ｌ（ｔ）の１階微分の値が極大となる時刻である。また、演算部２１２は、変曲点の時刻として、窓の開始時刻以前で、窓の開始時刻から最も近い変曲点の時刻（φｌ_１）、窓の開始時刻以降で窓の開始時刻から最も近い変曲点の時刻（φｌ_２）、及びφｌ_２の次の変曲点の時刻（φｌ_３）の３つを算出する（ステップＳ３９）。

【0054】

演算部２１２は、上記で算出した変曲点ごとの血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との時間差を算出する（ステップＳ４０）。具体的には、φｂ_１とφｌ_１との差を第１の個別時間差φ１、φｂ_２とφｌ_２との差を第２の個別時間差φ２、φｂ_３とφｌ_３との差を第３の個別時間差φ３とする。

【0055】

演算部２１２は、算出した個別時間差φ１〜φ３、重みｗ_１〜ｗ_３に基づいて、以下の式（４）から、補正値τを算出する。

【数7】

【0056】

すなわち、演算部２１２は、心拍ごとに、対応する個別時間差φ１〜φ３を算出する。そして、演算部２１２は、窓幅ｗに対する各心拍の影響度合いに基づいた重みｗ_１〜ｗ_３を個別時間差に乗じて、式（４）の各項に対応する調整後個別時間差を算出する。演算部２１２は、調整後個別時間差に基づいて、血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との時間差である補正値τを算出する。

【0057】

以下、血管剛性推定装置１は、実施の形態１と同様の第１推定工程（ステップＳ４１）、第２推定工程（ステップＳ４２）で、血管剛性βの推定を行う。この際、演算部２１２は、ステップＳ４０で算出した補正値τを用いて、血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との時間差を補正して用いる。具体的には、以下の式（５）に示すように、Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ−τ）とを用いて血管剛性βを推定する。

【数8】

【0058】

血管剛性推定装置１は、予め定められた所定の終了条件に達するまで（ステップＳ４３のＮＯ）、サンプリングステップを進めて窓幅ｗをずらす（ステップＳ４４）。そして、各サンプリングステップで、時間差補正工程、第１推定工程及び第２推定工程を実施して、血管剛性βの推定処理を行う。予め定められた所定の終了条件に達すると（ステップＳ４３のＹＥＳ）、血管剛性推定装置１は、推定結果を表示部２３に表示させるとともに、記憶部２２に記憶させる（ステップＳ４５）。

【0059】

以上説明したように、本実施の形態に係る血管剛性推定装置１では、血圧Ｐ_ｂ（ｔ）と容積脈波Ｐ_ｌ（ｔ）との時間差を補正して血管剛性βを推定するので、血圧Ｐ_ｂ（ｔ）の測定箇所と容積脈波Ｐ_ｌ（ｔ）の測定箇所との距離にかかわらず、精度よく血管剛性βを推定することが可能である。

【0060】

本実施の形態では、第１推定工程に係る窓幅ｗと第２推定工程の窓幅ｗは同じであることとしたが、これに限られない。例えば、第１推定工程の窓幅ｗと第２推定工程の窓幅ｗとが異なる場合、いずれかの工程の窓幅ｗを適宜選択して、時間差を表す補正値τを算出すればよい。

【0061】

また、本実施の形態では、３つの個別時間差に基づいて補正値τを算出することとしたが、これに限られない。例えば、窓幅ｗをより大きく設定した場合、被験者の心拍数が大きい場合等、窓幅ｗに係る心拍数に応じて個別時間差の数を増加させることとしてもよい。

【0062】

（数値例）
以下、実施の形態２の血管剛性推定方法を用いて、血管剛性βを推定し、血管剛性βのパワースペクトル密度（ＰＳＤ：Power Spectral Density）を算出した例について説明する。本例では、胸腔鏡下交感神経遮断術を受ける患者１５名を被験者とし、サンプリング周波数を５Ｈｚとして生体情報を取得して、右手の血管剛性βを推定した。解析区間は、図５に示す期間のうち右交感神経幹遮断直前の３０秒間（第１解析区間）と、右交感神経幹遮断６０秒後からの３０秒間（第２解析区間）である。また、ＰＳＤ算出法として、Ｂｕｒｇ法を用いる。

【0063】

図６（Ａ）は、右交感神経刺激中である第１解析区間のＰＳＤを示している。より詳細には、図６（Ａ）に示すグラフは、全体のパワー値の合計で除することにより、正規化されたＰＳＤに基づいて、被験者１５名分の正規化ＰＳＤを示している。また、被験者１５名分の正規化ＰＳＤの平均を太線で示している。従来の血管剛性推定方法では、心拍１拍ごとに、血管剛性βを推定していたので、血管剛性βの時間変化に含まれる周波数成分を詳細に分析することは難しかった。しかしながら、図６（Ａ）、（Ｂ）のＰＳＤのグラフに示すように、本発明に係る血管剛性推定方法によれば、詳細な周波数成分分析を行うことができる。例えば、図６（Ａ）、（Ｂ）のグラフから、１〜１．５Ｈｚ付近にピークを確認することができ、血管剛性βが、心拍リズムを反映していると推測できる。

【0064】

また、０．２Ｈｚ付近のピークから、血管剛性βが、呼吸リズムを反映していると推測できる。さらに、０．１Ｈｚ付近のピークから、血管剛性βと収縮期血圧の低周波成分との相関関係を分析することができる。

【0065】

以上、説明したように、本発明に係る血管剛性推定方法及び血管剛性推定装置によれば、心拍のタイミングに拘束されず、心拍より高いサンプリング周波数で血管剛性βを推定することができるので、より詳細に身体状態の分析を行うことができる。

【0066】

また、本実施の形態に係る血管剛性推定方法は、通常のコンピュータシステムを用いて実現可能である。例えば、上記の動作を実行するためのコンピュータプログラムを、ＵＳＢメモリ、ＤＶＤ−ＲＯＭ等のコンピュータが読み取り可能な記録媒体に格納して配布し、当該コンピュータプログラムをコンピュータにインストールすることにより、コンピュータ装置を上記の血管剛性推定方法を実行する推定装置として機能させることができる。

【0067】

また、本実施の形態では、センサ１０で取得した生体信号である動脈血圧、指尖容積脈波から、推定装置本体２０で血管剛性βを推定することとしたが、これに限られない。例えば、過去に取得され、記憶部２２に記憶されている生体信号を、演算部２１２に読み込んで、血管剛性βの推定を行うこととしてもよい。また、遠隔地の被験者に装着されたセンサ１０から、推定装置本体２０へネットワークを介して生体信号を読み込んで、血管剛性βを推定することとしてもよい。

【産業上の利用可能性】

【0068】

本発明は、身体の状態を把握するための血管剛性推定に好適である。特に、血管剛性の時間変化から、様々な周波数の信号に由来する身体状態の変化を推定する場合の血管剛性推定に好適である。

【符号の説明】

【0069】

１ 血管剛性推定装置、１０ センサ、２０ 推定装置本体、２１ 制御部、２１１ 生体信号取得部、２１２ 演算部、２２ 記憶部、２３ 表示部、２４ 入力部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】