血管状態評価装置、血管状態評価方法および血管状態評価プログラム
【課題】動脈硬化の度合いをより高い精度で評価することが可能な血管状態評価装置、血管状態評価方法および血管状態評価プログラムを提供する。
【解決手段】位相線傾き算出部（実測）３２は、周波数変換部３０ａおよび３０ｂから出力される位相特性Ｐａ（ｆ）およびＰｂ（ｆ）を受信し、両者間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出する。一方、位相線傾き算出部（モデル）３８は、伝達関数算出部３６で算出された伝達関数Ｇａ（ｆ）と伝達関数Ｇｂ（ｆ）との間の位相差特性を算出し、算出した位相差特性を探索部４０へ出力する。探索部４０は、変数ｋをフィッティングして、傾きｇ（ｋ）と傾きｇｅｘｐとが互いに略一致する変数ｋｏｐｔ（最適解）を決定する。この変数ｋｏｐｔが被験者の動脈硬化の度合いを示す指標となる。
【解決手段】位相線傾き算出部（実測）３２は、周波数変換部３０ａおよび３０ｂから出力される位相特性Ｐａ（ｆ）およびＰｂ（ｆ）を受信し、両者間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出する。一方、位相線傾き算出部（モデル）３８は、伝達関数算出部３６で算出された伝達関数Ｇａ（ｆ）と伝達関数Ｇｂ（ｆ）との間の位相差特性を算出し、算出した位相差特性を探索部４０へ出力する。探索部４０は、変数ｋをフィッティングして、傾きｇ（ｋ）と傾きｇｅｘｐとが互いに略一致する変数ｋｏｐｔ（最適解）を決定する。この変数ｋｏｐｔが被験者の動脈硬化の度合いを示す指標となる。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は生体を構成する血管の状態を評価する血管状態評価装置、血管状態評価方法および血管状態評価プログラムに関し、特に血管経路を伝達関数としてモデル化した上で、血管の動脈硬化の度合いを評価する技術に関するものである。
【背景技術】
【０００２】
近年、動脈硬化に起因する循環器系疾患が増加しており、これを受けて血管の動脈硬化の度合いを評価する評価装置が実用化されている。動脈硬化の度合いを評価する代表的な方法として、脈波速度法が知られている。この脈波速度法は、心臓の拍動に伴う血圧変化が血管上を伝播する速度（脈波速度）と血管の弾力度合い（剛性）との間に相関があることを利用したものである。すなわち、弾性管である血管内を脈波が進む際、管壁が硬く、内径が細く、管厚が厚いほど脈波速度は増大するため、脈波速度を測定することで血管の動脈硬化の度合いを知ることができる。特に、両上腕および両足首における血圧の時間波形を用いる、ｂａＰＷＶ法（brachial-ankle Pulse Wave Velocity method）による評価装置が実用化されている。
【０００３】
このような脈波速度の測定方法として、たとえば、特開２００６−３２６３３４号公報（特許文献１）には、複数の電圧波形のうち、隣り合う電圧電極対から得られた電圧波形間の時間的なひずみを検出する検出手段と、複数の電圧波形の全てについて、隣り合う電圧電極間の距離及び／又は時間的なずれとを用い、隣り合う電圧電極間における脈波伝播速度又は脈波伝播時間の変化率を求める算出手段とを有する脈波伝播速度測定装置が開示されている。
【特許文献１】特開２００６−３２６３３４号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
特開２００６−３２６３３４号公報（特許文献１）に記載されるように、脈波速度を測定する方法として、血管経路上の複数の点で測定した時間波形の間の時間差（遅れ時間）を求め、各点の心臓からの経路差を当該時間差で除することで脈波速度が算出される。
【０００５】
しかしながら、実際の脈波速度は伝播経路や周波数に依存するため、単純に経路差を時間差で除するだけでは、厳密な脈波速度を算出することはできない。すなわち、被評価者の血管径や血管長、および脈波に含まれる周波数成分などに依存して、本来の脈波速度からずれた値が算出される場合があった。そのため、動脈硬化の度合いの評価精度をあまり高められないという問題があった。
【０００６】
そこで、この発明は、かかる問題を解決するためになされたものであり、その目的は、動脈硬化の度合いをより高い精度で評価することが可能な血管状態評価装置、血管状態評価方法および血管状態評価プログラムを提供することである。
【課題を解決するための手段】
【０００７】
この発明のある局面に従う血管状態評価装置は、生体を構成する血管を複数の区間に分割してモデル化した循環系モデルを格納する記憶手段を備え、循環系モデルは、複数の区間の各々を代表する形状値を含む。さらに、血管状態評価装置は、生体の第１測定部位に装着されて第１生体信号の時間波形を測定する第１測定手段と、生体の第２測定部位に装着され、第１測定手段と同期して第２生体信号の時間波形を測定する第２測定手段と、第１生体信号と第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出する第１算出手段と、第１測定部位までの血管経路に対応して循環系モデルに基づいて規定された第１伝達関数と、第２測定部位までの血管経路に対応して循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出する第２算出手段とを備え、第１伝達関数および第２伝達関数は、血管の弾力度合いを示す弾力度変数を含む。さらに、血管状態評価装置は、第１算出手段によって算出される実測の位相差特性に基づいて、第２算出手段によって算出される位相差特性をフィッティングすることにより、弾力度変数を決定する探索手段を備える。
【０００８】
この発明によれば、循環系モデルに基づいて、第１および第２測定部位までの血管経路に対応してそれぞれ第１および第２伝達関数が規定された上で、第１伝達関数と第２伝達関数との間の位相差特性を算出する。ここで、第１および第２伝達関数は、血管の弾力度合いを示す弾力度変数を含む。一方で、第１および第２測定部位からそれぞれ第１および第２生体信号の時間波形が測定された上で、第１生体信号と第２生体信号との間の各周波数成分についての位相差に基づいて実測の位相差特性を算出する。さらに、実測の位相差特性に基づいて、第１伝達関数と第２伝達関数との間の位相差特性をフィッティングすることにより、弾力度変数を決定する。
【０００９】
ここで、第１および第２伝達関数は、脈波の周波数特性を含んで規定されるので、フィッティングによって、実測の測定信号の周波数特性を反映した弾力度変数の最適解が得られることになる。よって、周波数特性を考慮した上で、血管の状態（動脈硬化の度合い）をより高い精度で評価することができる。
【００１０】
好ましくは、血管状態評価装置は、第１および第２測定部位までの血管経路にそれぞれに対応する各区間の形状値に基づいて、第１および第２伝達関数を算出する伝達関数算出手段をさらに備える。
【００１１】
好ましくは、伝達関数算出手段は、各区間に対応付けられた、血管の圧力および血液流量を入力変数とする分布定数モデルを用いて、第１および第２伝達関数を算出し、分布定数モデルの各々は、対応する区間における血液の流れやすさに応じた縦インピーダンスと、弾力度変数を有する横インピーダンスとを含む。
【００１２】
好ましくは、血管状態評価装置は、探索手段によってフィッティングされた弾力度変数に基づいて、血管における脈波速度を算出する脈波速度算出手段をさらに備える。
【００１３】
さらに好ましくは、脈波速度算出手段は、第１測定部位までの血管経路に対応する各区間の形状値および第２測定部位までの血管経路に対応する各区間の形状値に基づいて脈波速度を算出する。
【００１４】
好ましくは、循環系モデルは、形状値として、血管径および血管長を含む。
好ましくは、循環系モデルは、生体を構成する血管を複数の区分に分類した上で、複数の区分の少なくとも１つの区分に属する血管をモデル化したものである。
【００１５】
さらに好ましくは、生体を構成する血管は、血管径の大きさに基づいて複数の区分に分類される。
【００１６】
また好ましくは、伝達モデル算出手段は、各区間に含まれる血管のうち、循環系モデルにおいてモデル化されていない血管をモデル化した末梢部モデルを、各区間に対応する循環系モデルに付加した上で、伝達関数を算出する。
【００１７】
さらに好ましくは、伝達モデル算出手段は、各区間の循環系モデルを血管の形状差に基づいて変換することで、当該区間の末梢部モデルを算出する。
【００１８】
さらに好ましくは、伝達モデル算出手段は、末梢部モデルの終端を無反射条件の下で伝達関数を算出する。
【００１９】
好ましくは、血管状態評価装置は、第１生体信号から各周波数成分についての位相を示す第１位相特性を算出する第１周波数変換手段と、第２生体信号から各周波数成分についての位相を示す第２位相特性を算出する第２周波数変換手段とをさらに備え、第１算出手段は、第１位相データと第２位相データとを差分して差分位相データを算出する。さらに、血管状態評価装置は、差分位相データにおける周期遅れに起因する位相ずれを、１または２以上の周期に相当する位相の単位で補正することで、実測の位相差特性を算出する。
【００２０】
好ましくは、第１算出手段は、第１生体信号と第２生体信号との間のコヒーレンス値が予め定められたしきい値より高い周波数成分を用いて実測の位相差特性を算出する。
【００２１】
この発明の別の局面に従えば、生体を構成する血管を複数の区間に分割してモデル化した循環系モデルを用いて、生体を構成する血管の状態を評価する血管状態評価方法であって、循環系モデルは、複数の区間の各々を代表する形状値を含む。血管状態評価方法は、生体の第１測定部位から第１生体信号の時間波形を測定するとともに、生体の第２測定部位から第２生体信号の時間波形を測定するステップと、第１生体信号と第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するステップと、第１測定部位までの血管経路に対応して循環系モデルに基づいて規定された第１伝達関数と、第２測定部位までの血管経路に対応して循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出するステップとを備え、第１伝達関数および第２伝達関数は、血管の弾力度合いを示す弾力度変数を含む。さらに、血管状態評価方法は、実測の位相差特性に基づいて、第１伝達関数と第２伝達関数との間の位相差特性をフィッティングすることにより、弾力度変数を決定するステップを備える。
【００２２】
この発明のさらに別の局面に従えば、演算処理部を有するコンピュータに、生体を構成する血管を複数の区間に分割して予めモデル化した循環系モデルを用いて、生体を構成する血管の状態の評価を実行させるための血管状態評価プログラムであって、循環系モデルは、複数の区間の各々を代表する形状値を含む。血管状態評価プログラムは、演算処理部が、生体の第１測定部位で測定される第１生体信号の時間波形を取得するとともに、生体の第２測定部位で測定される第２生体信号の時間波形を取得するステップと、演算処理部が、第１生体信号と第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するステップと、演算処理部が、第１測定部位までの血管経路に対応して循環系モデルに基づいて規定された第１伝達関数と、第２測定部位までの血管経路に対応して循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出するステップとを備え、第１伝達関数および第２伝達関数は、血管の弾力度合いを示す弾力度変数を含む。さらに、血管状態評価プログラムは、演算処理部が、演算処理部が、実測の位相差特性に基づいて、第１伝達関数と第２伝達関数との間の位相差特性をフィッティングすることにより、弾力度変数を決定するステップを備える。
【発明の効果】
【００２３】
この発明によれば、動脈硬化の度合いをより高い精度で評価することが可能な血管状態評価装置、血管状態評価方法および血管状態評価プログラムを実現できる。
【発明を実施するための最良の形態】
【００２４】
この発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。
【００２５】
［実施の形態１］
（装置構成）
図１は、この発明の実施の形態１に従う血管状態評価装置１００の概略構成図である。
【００２６】
図１を参照して、この発明の実施の形態１に従う血管状態評価装置（以下、単に「評価装置」とも記す）１００は、制御部２と、表示部４と、操作部６と、測定部２０ａ，２０ｂとを含む。
【００２７】
制御部２は、評価装置１００全体の制御を行なう装置であり、代表的に、ＣＰＵ（Central Processing Unit）１０と、ＲＯＭ（Read Only Memory）１２と、ＲＡＭ（Random Access Memory）１４とを含むコンピュータで構成される。
【００２８】
ＣＰＵ１０は、演算処理部に相当し、ＲＯＭ１２に予め格納されているプログラムを読出して、ＲＡＭ１４をワークメモリとして使用しながら、当該プログラムを実行する。また、ＲＯＭ１２には、少なくとも後述する循環系モデルが予め格納されており、ＣＰＵ１０は、本実施の形態に従う血管状態評価方法を格納したプログラムの実行に際して、当該循環系モデルを参照する。
【００２９】
また、制御部２には、表示部４および操作部６が接続されている。表示部４は、ユーザによる各種設定の入力を促したり、制御部２からの演算結果を表示したりする。これに対して、ユーザは、表示部４に表示される内容を確認しながら操作部６を操作して、所望の設定入力や操作を行なう。なお、表示部４は、一例として、ＬＥＤ（Light Emitting Diode）やＬＣＤ（Liquid Crystal Display）などからなる。
【００３０】
より具体的には、制御部２は、測定部２０ａ，２０ｂに対して測定指令を与えるとともに、当該測定指令に応答して測定された測定信号Ｐａ（ｔ），Ｐｂ（ｔ）を受信し、当該測定信号Ｐａ（ｔ），Ｐｂ（ｔ）に基づいて、本実施の形態に従う血管状態評価方法を実行する。
【００３１】
測定部２０ａ，２０ｂは、被験者２００の所定の測定部位に装着された押圧カフ（空気袋）２４ａ，２４ｂの内圧（以下、「カフ圧」という）を加圧して、それぞれの測定部位における生体信号（一例として、脈波）の時間波形を測定する。なお、後述するように、制御部２は、測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するので、制御部２からは、測定部２０ａおよび２０ｂが互いに同期して生体信号を測定できるように、測定指令が同時に与えられる。
【００３２】
より詳細には、たとえば、押圧カフ２４ａおよび２４ｂは、それぞれ被験者２００の足首部および上腕部に装着され、それぞれ配管２２ａおよび２２ｂを介して測定部２０ａおよび２０ｂから供給される空気によって加圧される。そして、この加圧によって、押圧カフ２４ａおよび２４ｂは対応の測定部位に押圧され、当該測定部位の脈波に応じた圧力変化がそれぞれ配管２２ａおよび２２ｂを介して測定部２０ａおよび２０ｂへ伝達される。測定部２０ａ，２０ｂは、この伝達される圧力変化を検出すことで、測定部位の脈波の時間波形を測定する。なお、測定信号Ｐａ（ｔ）およびＰｂ（ｔ）の所定の周波数成分（一例として、０〜２０［Ｈｚ］）に対して演算処理を行なうことが好ましいので、測定信号Ｐａ（ｔ）およびＰｂ（ｔ）の測定周期（サンプリング周期）は、この周波数成分に応じた時間間隔（一例として、２５ｍｓｅｃ）より短くすることが好ましい。
【００３３】
このような測定動作を実行するために、測定部２０ａは、圧力センサ２８ａと、調圧弁２６ａと、圧力ポンプ２５ａと、配管２７ａとを含む。圧力センサ２８ａは、配管２２ａを介して伝達される圧力変動を検出するための検出部位であり、一例として、単結晶シリコンなどからなる半導体チップに所定間隔に配列された複数のセンサエレメントを含む。調圧弁２６ａは、圧力ポンプ２５ａと押圧カフ２４ａとの間に介挿され、測定時に押圧カフ２４ａを加圧に用いられる圧力を所定の範囲に維持する。圧力ポンプ２５ａは、制御部２からの測定指令に応じて作動し、押圧カフ２４ａを加圧するための加圧空気を供給する。
【００３４】
同様に、測定部２０ｂは、圧力センサ２８ｂと、調圧弁２６ｂと、圧力ポンプ２５ｂと、配管２７ｂとを含む。各部の構成については、測定部２０ａと同様であるので、詳細な説明は繰返さない。
【００３５】
なお、本実施の形態では、生体信号の一例として、脈波によって生じる圧力変化を、圧力カフを用いて測定する構成について説明するが、たとえば、被験者２００の測定部位に微少の一定電流を流すとともに、脈波の伝播に応じて生じるインピーダンス（生体インピーダンス）の変化によって生じる電圧変化を測定するようにしてもよい。
【００３６】
なお、図１に示す評価装置１００と本願発明との対応関係については、測定部２０ａ，２０ｂと、配管２２ａ，２２ｂと、押圧カフ２４ａ，２４ｂとが「第１および第２の測定手段」に相当する。
【００３７】
（機能ブロック図）
制御部２は、予め格納する循環系モデルに基づいて、それぞれ押圧カフ２４ａおよび２４ｂが装着された測定部位までの血管経路に対応して規定される２つの伝達関数を算出する。ここで、それぞれの伝達関数は、血管の弾力度合いを示す弾力度変数を含む。すなわち、弾力度変数は、血管の動脈硬化の度合いを示す指標である。本実施の形態では、弾力度変数の代表例として「ヤング率」を用いるが、血管の剛性や柔軟度を示す他の変数を用いてもよい。そして、制御部２は、測定信号Ｐａ（ｔ）およびＰｂ（ｔ）を周波数領域の信号に変換した上で、両者間の実測の位相差特性を算出し、この実測の位相差特性が上述した２つの伝達関数の間の位相差特性と一致するように、弾力度変数をフィッティング（同定）する。このフィッティングされた弾力度変数が被験者２００の動脈硬化の度合いを示す値となる。以下、このような制御部２における処理動作を実現するための機能ブロックについて説明する。
【００３８】
図２は、この発明の実施の形態１に従う血管状態評価装置１００の制御部２で実行される機能を模式的に示す機能ブロック図である。
【００３９】
図２を参照して、制御部２は、周波数変換部（ＦＦＴ）３０ａ，３０ｂと、位相線傾き算出部（実測）３２と、記憶部３４と、伝達関数算出部３６と、位相線傾き算出部（モデル）３８と、探索部４０と、評価部４２とを実現する。
【００４０】
周波数変換部３０ａおよび３０ｂは、それぞれ時間波形である測定信号Ｐａ（ｔ）およびＰｂ（ｔ）を所定期間にわたって受信し、当該受信した測定信号Ｐａ（ｔ）およびＰｂ（ｔ）を周波数領域の関数に変換する周波数変換部である。代表的に、周波数変換部３０ａおよび３０ｂは、高速フーリエ変換（ＦＦＴ：Fast Fourier Transformer）を用いて、周波数変換を実行する。なお、高速フーリエ変換に限らず、時間領域の関数をフーリエ級数などの周波数領域の関数に変換するものであれば、いずれのロジックを用いてもよい。
【００４１】
そして、周波数変換部３０ａは、測定信号Ｐａ（ｔ）の各周波数成分についての位相を示す位相特性Ｐａ（ｆ）を算出し、算出した位相特性Ｐａ（ｆ）を位相線傾き算出部（実測）３２へ出力する。同様に、周波数変換部３０ｂは、測定信号Ｐｂ（ｔ）の各周波数成分についての位相を示す位相特性Ｐｂ（ｆ）を算出し、算出した位相特性Ｐｂ（ｆ）を位相線傾き算出部（実測）３２へ出力する。
【００４２】
位相線傾き算出部（実測）３２は、ユーザによる操作部６（図１）などの操作に応答して、測定指令を測定部２０ａ，２０ｂへ与える。この測定指令を与えた後、周波数変換部３０ａおよび３０ｂから出力される位相特性Ｐａ（ｆ）および位相特性Ｐｂ（ｆ）を受信し、両者間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出する。具体的には、位相線傾き算出部（実測）３２は、周波数成分毎に位相特性Ｐａ（ｆ）およびＰｂ（ｆ）の値を比較し、両者の位相差を算出する。後述するように、このように算出された位相差は、周波数成分についての一次関数として近似できるので、この近似された一次関数（位相線）の傾きｇｅｘｐ［ｄｅｇ／Ｈｚ］を、実測の位相差特性として探索部４０へ出力する。すなわち、偏角φｅｘｐ＝∠位相特性Ｐａ（ｆ）／位相特性Ｐｂ（ｆ）として算出される偏角φを用いて、傾きｇｅｘｐ＝ｔａｎ（φｅｘｐ）として定義できる。
【００４３】
一方、伝達関数算出部３６は、心臓からそれぞれ押圧カフ２４ａおよび２４ｂが装着された２つの測定部位までの血管経路の伝達特性を示す２つの伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）を算出し、位相線傾き算出部（モデル）３８へ出力する。より具体的には、伝達関数算出部３６は、記憶部３４に予め格納されている循環系モデルに基づいて、心臓を入力端とした全身についての脈波伝播モデル（伝達関数）を算出した上で、全身の脈波伝播モデルの中で２つの測定部位までの血管経路に対応する伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）を算出する。この際、伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）には、変数ｋを含む形でヤング率が組込まれるが、変数ｋには探索部４０によって具体な値が設定される。
【００４４】
記憶部３４は、被験者２００を構成する血管を複数の区間に分割してモデル化した循環系モデルを格納している。この循環系モデルは、各区間に対応付けて、各区間を代表する形状値が規定されている。このような形状値の一例として、本実施の形態では、各区間の血管径、血管長、血管壁の厚さなどが含まれる。なお、循環系モデルについては、後に詳述する。
【００４５】
位相線傾き算出部（モデル）３８は、伝達関数Ｇａ（ｆ）と伝達関数Ｇｂ（ｆ）との間の位相差特性を算出し、算出した位相差特性を探索部４０へ出力する。具体的には、位相線傾き算出部（モデル）３８は、周波数領域における位相特性Ｐａ（ｆ）と位相特性Ｐｂ（ｆ）との間の位相差である位相線の傾きｇ（ｋ）［ｄｅｇ／Ｈｚ］を、位相差特性として探索部４０へ出力する。ここで、傾きｇ（ｋ）は、偏角φｍｏｄｅｌ＝∠伝達関数Ｇａ（ｆ）／伝達関数Ｇｂ（ｆ）として算出される偏角φｍｏｄｅｌを用いて、傾きｇ（ｋ）＝ｔａｎ（φｍｏｄｅｌ）として定義できる。
【００４６】
探索部４０は、位相線傾き算出部（実測）３２によって算出される傾きｇｅｘｐに基づいて、位相線傾き算出部（モデル）３８によって算出される傾きｇ（ｋ）をフィッティングすることにより、変数ｋを決定する。すなわち、傾きｇ（ｋ）と傾きｇｅｘｐとが互いに略一致するまで、変数ｋが初期値ｋ０から順次変更され、伝達関数算出部３６および位相線傾き算出部（モデル）３８における演算処理が繰返し実行される。このようにして、傾きｇ（ｋ）と傾きｇｅｘｐとが互いに略一致する変数ｋｏｐｔ（最適解）が決定されると、探索部４０は、決定した変数ｋの値を評価部４２へ出力する。この決定された変数ｋの最適解が、被験者２００の動脈硬化の度合いを示す指標となる。
【００４７】
評価部４２は、探索部４０で決定された最適解ｋｏｐｔ（または、最適解ｋｏｐｔを用いて換算されたヤング率）を予め定められた基準値と比較し、動脈硬化の度合いについての評価を、表示部４（図１）などに出力する。
【００４８】
以下、このような主要機能の動作や構成について詳述する。
（物理モデル）
上述したように、伝達関数算出部３６は、心臓を入力端（始点）とする血管経路の伝達特性を示す伝達関数を算出するが、この算出される伝達関数は、脈波が血管を伝播する力学的なモデルから解析的に算出される。本実施の形態では、血管の各区間を一次元線形分布定数モデル化して、伝達関数を算出する構成について例示する。
【００４９】
まず、血管を軸対象の微小変形する薄肉円管とし、内部の血液の流れを非粘性流体の層流とし、かつ反射波はないと仮定してモデル化すると、脈波速度Ｃｐと血管壁のヤング率Ｅとの関係は、Moens-Kortewegの式と呼ばれる（１）式で表される。なお、脈波速度Ｃｐは、心臓の拍動に伴う血圧変化が血管上を伝播する速度を意味する。
【００５０】
【数１】
【００５１】
この（１）式から、血管が硬く、内腔が細く、または血管壁が厚いほど脈波速度Ｃｐが増加することがわかる。
【００５２】
図３は、血管内の血液の一次元流れモデルを示す図である。
一般に、血管に比べて血液の体積弾性率は十分高いので、血管を弾性円管とし、血液を非圧縮性流体として考えることができる。このような弾性管内における一次元流れの支配方程式は、以下のように導かれる。
【００５３】
図３を参照して、一次元流れモデルの断面ＣＳ１−ＣＳ２間の検査体積５０に関する質量の保存について考える。断面ＣＳ１の内腔の面積をＡ（＝πｒｉ２），流体（血液）の密度をρ，圧力をｐ，断面平均流速をＵとし、断面ＣＳ１−ＣＳ２間にある分枝血管へ単位時間に漏れ出す流体の体積を単位長さおよび単位圧力あたりＧとすると，質量保存則より（２）式が成り立つ。ここで、非圧縮性流体では密度ρは一定であるので（２）式は（３）式のように簡略化できる。
【００５４】
【数２】
【００５５】
図４は、図３に示す検査体積５０に働く力５２および出入りする運動量５４を示す図である。
【００５６】
図４を参照して、検査体積５０内における運動量５４の単位時間当たりの変化は、流入する正味の運動量５４と検査体積５０に及ぼされる力５２に等しいので，高次の微少項を省略して（４）式が導かれる。
【００５７】
【数３】
【００５８】
（４）式を連続の式を用いて整理すると（５）式に示す運動方程式が得られる．
【００５９】
【数４】
【００６０】
次に、血管を一次元線形分布定数モデル化するために、（３）式および（５）式において非線形項を省略し、変数を圧力ｐと体積流量ｑ（＝ＡＵ）に置き換えることにより、（６）式および（７）式が得られる。
【００６１】
【数５】
【００６２】
ここで、（６）式および（７）式の４個の係数の物理的な意味については、Ｒは血液が流れるときの粘性抵抗を示し、Ｌは流れが変化するときに急な変化を妨げようとする血液の慣性を示し、Ｇは血管外もしくは分岐管に流れ出す血液の流れやすさを示し、Ｃは圧力変化に応じて血管が伸縮する際に血管内に血液を蓄える能力を示す。
【００６３】
図５は、血管を一次元線形分布定数モデル化した模式図である。図５（ａ）は、（６）式および（７）式と血管の物理モデルとを対応付ける図である。図５（ｂ）は、図５（ａ）に示す物理モデルを電気的な等価回路に置き換えた図である。
【００６４】
すなわち、（６）式および（７）式は、図５（ａ）に示すような物理モデルと対応付けることができる。さらに、（６）式および（７）式において、圧力ｐを電圧ｖに置き換え、流量ｖを電流ｉに置き換えることで、図５（ｂ）に示すような電気的な等価回路（分布定数回路）に置き換えることができる。ここで、Ｒは抵抗を示し、Ｌはインダクタンスを示し、Ｇはアドミタンスを示し、Ｃはキャパシタンスを示す。
【００６５】
ここで、（６）式についてみれば、血管系では運動方程式に対応する一方、電気系ではオームの法則に対応する。そして、血管系において、断面ＣＳ１と断面ＣＳ２との間の圧力勾配によって流体が加速される現象が、電気系において、インダクタンスの両端に印加される電位差が電流を引き起こす現象に対応することを意味している。
【００６６】
また、（７）式についてみれば、血管系では連続の式（質量保存則）に対応する一方、電気系では電荷の保存則に対応する。そして、血管系において、断面ＣＳ１から断面ＣＳ２へ進めない質量の滞留分が血管を押し広げて圧力の上昇を引き起こす現象が、キャパシタに溜まった電荷が電圧の上昇を引き起こす現象に対応することを意味している。
【００６７】
さらに、（６）式および（７）式において、ｐ＝Ｐｅｊｗｔとし、ｑ＝Ｑｅｊｗｔとすると、それぞれ（８）式および（９）式に示す関係式が導出される。
【００６８】
【数６】
【００６９】
本明細書では、図５（ｂ）および（８）式に示すＺｌ（＝ｒ＋ｊωＬ）を「縦インピーダンス」と称し、図５（ｂ）および（９）式に示すＺｔ（＝（Ｇ＋ｊωＣ）−１）を「横インピーダンス」と称す。（８）式および（９）式の一般解は、ｘ＝０における圧力の進行波の振幅値をＰｆとし、後退波の振幅値をＰｒとして、それぞれ（１０）式および（１１）式となる。なお、角周波数ωと周波数ｆとの間には、ω＝２πｆの関係が成立する。
【００７０】
【数７】
【００７１】
また、伝播定数γは、減衰定数βと位相速度（脈波速度）Ｃｐとを用いて、（１２）式のように表される。
【００７２】
【数８】
【００７３】
ここで、位相速度Ｃｐは脈波が単位時間に進む距離を示す量であり、減衰定数βは脈波の振幅が単位距離を進むごとにｅ−β倍になることを示す。また、特性インピーダンスＺ０は（１３）式のように表せ、単位体積の脈波を進行方向へ進めるために必要な圧力を示す。
【００７４】
【数９】
【００７５】
さらに、距離ｌｓｅを隔てた二点での圧力Ｐｓ，Ｐｅと体積流量Ｑｓ，Ｑｅは、（１４）式の伝達行列で接続される。
【００７６】
【数１０】
【００７７】
本実施の形態では、（１４）式に示す伝達行列を血管の各区間に対応付けて算出するとともに、対象とする血管経路に応じて、各区間に対応する伝達行列を縦続接続することで、伝達関数を算出する。このとき、任意の境界より下流の条件はその境界での圧力Ｐｘ，体積速度Ｑｘとの比である（１５）式のインピーダンスＺｘで表される。
【００７８】
【数１１】
【００７９】
また、進行波と後退波の振幅の比である反射率Ｓｐは、（１６）式で表される。
【００８０】
【数１２】
【００８１】
（縦インピーダンスの算出）
縦インピーダンスＺｌは、流体の粘性抵抗と慣性の項とからなり、血管断面内の流速分布をモデル化することで求められる。
【００８２】
本実施の形態では、ウマースリーモデルに基づいて、縦インピーダンスを算出する。このウマースリーモデルは、ニュートン流体の円管内脈波流が十分に発達した状態での流速分布を表したものである。このウマースリーモデルに基づく縦インピーダンスは、第１種ベッセル関数Ｊｎを用いて（１７）式で表せる。
【００８３】
【数１３】
【００８４】
ここで、（１７）式中のαは「ウマースリーのアルファ」と称され、脈波流の粘性項と慣性項の比を示す量であり、定常流におけるレイノルズ数に相当する。また、血液の密度ρは、代表的に１．０３×１０３［ｋｇ／ｍ３］とし、血液の粘性係数μは、代表的に４×１０−３［Ｐａ・ｓ］とした。
【００８５】
なお、（１７）式に示すウマースリーモデルに代えて、非粘性モデルを用いてもよい。このモデルは、血液を非粘性流体とし、断面内流速を一定とするものである。この非粘性モデルに基づく縦インピーダンスは、（１８）式で表せる。
【００８６】
【数１４】
【００８７】
さらに、上記のモデルに代えてポアゼイユモデルを用いてもよい。このモデルは、ニュートン流体の円管内定常流が十分に発達した状態での流速分布を表したものである。このポアゼイユモデルに基づく縦インピーダンスは、（１９）式で表せる。
【００８８】
【数１５】
【００８９】
（横インピーダンスの算出）
横インピーダンスは、遺漏または分岐の項Ｇと、管のコンプライアンスの項Ｃとからなる。
【００９０】
遺漏または分岐の項については、血管壁から周囲組織への遺漏も分岐もない場合には、Ｇ＝０とする。これに対して、分岐がある場合には、分岐管のアドミタンスをＧとする。
【００９１】
次に、管のコンプライアンスの項については、厚肉円管をモデル化したコンプライアンスを用いることができる。外圧および軸方向ひずみ一定の条件における厚肉円管の軸対称微小変形のコンプライアンスは、（２０）式で表される。
【００９２】
【数１６】
【００９３】
ここで、血管壁のポアソン比νは、代表的に０．５とした。
なお、（２０）式に示す厚肉円管をモデル化したコンプライアンスに代えて、薄肉円管をモデル化したコンプライアンスを用いてもよい。外圧および軸方向ひずみ一定の条件における薄肉円管の軸対称微小変形のコンプライアンスは、（２１）式で表される。
【００９４】
【数１７】
【００９５】
（循環系モデル）
この発明の実施の形態１に従う血管状態評価装置１００で用いる循環系モデルは、生体を構成する血管を複数の区間に分割してモデル化したものである。このような循環系モデルの代表的なものとして、参考文献１「Avolio,A.P, Multi-branched Model of Human Arterial System, 1980, Med. & Biol. Engng. & Comp., 18,796」に記載されている、いわゆる「Ａｖｏｌｉｏモデル」が知られており、本実施の形態においても、循環系モデルとして、このＡｖｏｌｉｏモデルを採用する。
【００９６】
図６は、Ａｖｏｌｉｏモデルの模式図である。
図６を参照して、Ａｖｏｌｉｏモデルでは、全身の動脈を１２８の血管要素（区間）に分割し、各区間を代表する形状値を規定している。Ａｖｏｌｉｏモデルの各区間に対応付けられた形状値の一部を付表に示す。このＡｖｏｌｉｏモデルは、形状値として、各区間に対応付けられた、長さ・半径・管壁の厚さ・ヤング率を含む。なお、Ａｖｏｌｉｏモデル中のヤング率は一応の基準値であり、後述するフィッティング処理においては、この基準値に変数ｋを乗じた値が用いられる。
【００９７】
この循環系モデルは、生体を構成するさまざまな血管を複数の区分に分類した上で、当該複数の区分の少なくとも１つの区分に属する血管をモデル化したものである。代表的に、血管は、その血管径の大きさに基づいて、血管径の大きいものから順に大動脈、中動脈（φ３．２ｍｍ以上）、小動脈（φ０．５ｍｍ以上）、細動脈（φ０．０３ｍｍ以上）、毛細血管などに区分される。そして、Ａｖｏｌｉｏモデルは、これらの区分のうち、大動脈および中動脈に区分される血管についてモデル化したものである。
【００９８】
なお、血管の区分の方法については、血管径の大きさに限らず、別の指標に基づいて区分してもよい。
【００９９】
伝達関数算出部３６（図２）は、記憶部３４（図２）に予め格納されているこのような循環系モデルを参照して、（１７）式および（２０）式に従って、各区間の縦インピーダンスおよび横インピーダンスを算出する。さらに、伝達関数算出部３６は、算出した縦インピーダンスおよび横インピーダンスを用いて、（１２）式，（１３）式，（１４）式に従って各区間の伝達行列を算出し、各区間の実際の接続関係に対応付けて、各伝達関数を縦続接続および並列接続することで、心臓を基準点とした全身についての脈波伝播モデル（伝達関数）を算出する。より具体的には、（１４）式に示す２行×２列の伝達行列が、各区間の接続関係（連続、分岐、および終端など）に応じて、順次接続されていく。
【０１００】
そして、伝達関数算出部３６は、算出された全身についての脈波伝播モデル（伝達関数）の中で、２つの測定部位までの血管経路に対応する伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）を算出する。また、終端部では、その反射率に応じて（１５）式の制約が加えられる。
【０１０１】
なお、心臓から吐出される圧力（圧力Ｐｓ）および体積流量（体積流量Ｑｓ）は未知数となるが、本実施の形態では、伝達関数Ｇａ（ｆ）と伝達関数Ｇｂ（ｆ）との間の位相差特性を算出すれば十分であるので、このような未知数が存在していても、互いにキャンセルできる。
【０１０２】
（末梢部モデル）
このような、全身についての脈波伝播モデル（伝達関数）ならびに伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）の算出に際して、上述のＡｖｏｌｉｏモデルに末梢部モデルを加えることが好ましい。これは、Ａｖｏｌｉｏモデルは、比較的太い血管（大動脈および中動脈）について厳密に形状値を与えている一方で、終端条件については末梢血管を模擬した一定の反射率を規定するにとどまっているからである。そのため、より高い評価精度を得るためには、Ａｖｏｌｉｏモデルにおいてモデル化されていない血管（小動脈、細動脈、毛細血管）を考慮することが好ましい。このようなＡｖｏｌｉｏモデルにおいてモデル化されていない血管をモデル化したモデル（以下、「末梢部モデル」と称す）を、Ａｖｏｌｉｏモデルから算出される伝達行列に付加した上で、全身についての脈波伝播モデル（伝達関数）などを算出する構成について説明する。
【０１０３】
このような末梢部モデルは、末梢血管の形状値と、当該末梢血管の上流側に接続されている区間の形状値との間の形状差を利用して算出する。本実施例では、形状差として、代表的に各血管の総断面積の差を用いる。
【０１０４】
図７は、末梢部モデルの模式図である。
図７を参照して、各血管（動脈）が循環系全体で有する総断面積は、分岐を経て下流の細い血管へ向かうにつれて増加する。特に、中動脈から細動脈までの分岐に伴う血管総断面積増加率は、参考文献２「William F.Ganong, Review of Medical Physiology 15ed.」によれば、２０倍程度と報告されている。
【０１０５】
そこで、本実施の形態では、小動脈および細動脈を末梢部モデルの対象とする。そして、一例として、分岐に伴う血管総断面積増加率は、中動脈から小動脈の過程で４倍とし、小動脈から細動脈の過程で５倍であるとし、小動脈および細動脈の長さはそれぞれ１０ｃｍおよび５ｃｍとする。また、小動脈および細動脈の血管径については、参考文献２に記載の一般的な値を採用し、それぞれの血管壁の厚さは、上流に接続されている中動脈の血管径と血管壁の厚さとの比に従って決定する。さらに、小動脈および細動脈のヤング率については、上流に接続されている中動脈のヤング率と同じ値を採用する。
【０１０６】
たとえば、橈骨動脈（図６に示すＡｖｏｌｉｏモデルの区間番号８８または９３）と、それに付加すべき末梢部モデルの形状値を以下の表に示す。
【０１０７】
【表１】
【０１０８】
（末梢部モデルの終端条件）
上述した末梢部モデルにおける細動脈の終端条件は、以下に説明するように任意に設定できる。これは、末梢血管によって構成される中動脈終端での反射率（以下、末梢反射率と称する）が細動脈の終端条件に依存しないためである。そこで、本実施の形態では、末梢部モデルにおける終端を無反射条件の下で、全身についての脈波伝播モデル（伝達関数）などを算出する。
【０１０９】
図８は、表１に示す形状値を用いて、図７に示す細動脈終端での反射率ＳＴを変化させた場合の末梢反射率の変化を示す図である。
【０１１０】
図８（ａ）は、反射率ＳＴ＝０（無反射）の場合を示し、図８（ｂ）は、反射率ＳＴ＝１（閉鎖端）の場合を示し、図８（ｃ）は、反射率ＳＴ＝−１（開放端）の場合を示す。
【０１１１】
図８（ａ）〜図８（ｃ）を参照して、細動脈終端での反射率ＳＴにかかわらず、振幅特性および周波数特性のいずれについても、略同一の挙動を示していることが分かる。すなわち、細動脈の終端条件は末梢反射率に寄与していないといえる。
【０１１２】
この現象を物理的に解釈できるように、各動脈内における脈波の伝播性状についても説明する。
【０１１３】
図９は、中動脈・小動脈・細動脈を伝播する脈波の減衰定数を数値計算によって求めた結果を示す図である。
【０１１４】
図９を参照して、脈波が伝播する血管の管径が小さくなるにつれて、減衰定数が指数的に大きくなっていることが分かる。これは、細い管内では、脈波（流体）の管壁摩擦抵抗の影響が大きくなるためであると考えられる。すなわち、細動脈終端に到達する脈波の振幅は、中動脈を伝播する脈波の振幅に比較して極めて小さくなっており、細動脈終端からの反射波は十分に減衰していると考えることができる。
【０１１５】
図１０は、図７に示す中動脈の始端（座標ｘ）から所定の周波数で加振（加圧）した場合の血管内圧力分布を示す図である。なお、図１０の座標系は、図７に示す中動脈の始端を原点（ｘ＝０）とし、下流方向をｘの正方向としたものである。
【０１１６】
図１０には、加振周波数ｆを１Ｈｚ，５Ｈｚ，１０Ｈｚ，２０Ｈｚの４種類とした場合の結果を示す。図１０から分かるように、いずれの加振周波数においても細動脈内部を伝播する脈波はその終端付近において十分に減衰している。これは、上述したように、細動脈内における脈波の距離減衰が大きいためである。
【０１１７】
したがって、細動脈の終端にどのような境界条件を設けても、細動脈の終端からの反射波が上流の循環系に影響を及ぼすことはないと考えられる。すなわち、末梢反射率は、末梢部モデルにおける終端の境界条件に依存することなく、末梢血管の形状値のみによって決定される。
【０１１８】
（位相差特性）
図１１は、一様な管路における脈波伝播の様子を示す模式図である。
【０１１９】
図１１を参照して、反射波が存在せず、脈波速度Ｃｐが（１）式で与えられる周波数に依存しない定数であるとする。すると、測定部位Ｍｐａの測定部位Ｍｐｂに対する脈波の位相遅れφは（２２）式で表される。
【０１２０】
【数１８】
【０１２１】
（２２）式を脈波速度Ｃｐおよび周波数ｆを用いて書き直せば、（２３）式が得られる。
【０１２２】
【数１９】
【０１２３】
（２２）式より、測定部位Ｍｐａ−測定部位Ｍｐｂ間の位相線図（位相差特性）は、周波数ｆの一次関数となり、その勾配は脈波速度Ｃｐに応じた値になることが分かる。さらに、（１）式および（２３）式を用いて（２４）式が得られる。
【０１２４】
【数２０】
【０１２５】
（２４）式より、血管壁のヤング率Ｅが大きくなるほど位相線図の勾配が緩やかになることが分かる。
【０１２６】
（フィッティング）
再度、図２を参照して、上述したように２つの測定部位間の位相差特性は周波数ｆの一次関数となるので、探索部４０は、この一次関数の傾きが実測値と一致するように、モデルをフィッティングする。
【０１２７】
より具体的には、伝達関数算出部３６は、Ａｖｏｌｉｏモデル中の各区間のヤング率Ｅｎ（ｎ＝１〜１２８）に変数ｋを乗じて得られる、仮のヤング率ｋ・Ｅｎを用いて伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）を算出する。そして、探索部４０は、位相線傾き算出部（実測）３２が算出する位相線の傾きｇｅｘｐと、位相線傾き算出部（モデル）３８が算出する位相線の傾きｇ（ｋ）との偏差Δ（＝｜ｇｅｘｐ−ｇ（ｋ）｜）が最小化するように、変数ｋを最適化する。この最適化処理は、代表的に数理計画法（たとえば、最小二乗法）を用いて行なわれるが、数理計画法については公知であるので、詳細な説明は行なわない。
【０１２８】
（実測の位相差特性の算出）
以下では、本実施の形態に従う血管状態評価装置１００を用いて、２人の被験者２００ａおよび２００ｂについて、実際に測定を行なった結果を示す。
【０１２９】
図１２は、被験者２００ａ，２００ｂの上腕および足関節に押圧カフを装着して測定した圧力の時間波形を示す図である。なお、測定信号Ｐａ（ｔ）が足関節における圧力を示し、測定信号Ｐｂ（ｔ）が上腕における圧力を示す。
【０１３０】
図１３は、図１２に示す上腕の圧力波形と足関節の圧力波形との間のコヒーレンスを示す図である。
【０１３１】
図１２（ａ）および図１３（ａ）は、被験者２００ａの測定結果を示し、図１２（ｂ）および図１３（ｂ）は、被験者２００ｂの測定結果を示す。
【０１３２】
ここで、コヒーレンスとは、波形間の周波数領域での相関性を示す指標であり、コヒーレンスが大きいほど（１に近いほど）両波形の相関性が高いことになる。なお、図１３に示すコヒーレンスは、（２５）式に従って算出した。
【０１３３】
【数２１】
【０１３４】
図１３（ａ）および図１３（ｂ）に示すように、特定の周波数においてコヒーレンスが低くなっていることが分かる。このようなコヒーレンス低下の力学的な要因としては、様々なものが考えられるが、代表的に、人体が有する非線形性や、測定部位での圧力が特定の周波数において節になってしまうことなどが影響していると考えられる。また、被験者の姿勢や測定中のわずかな動きといった人為的な要因によってもコヒーレンスが低下することも考えられる。
【０１３５】
このようなコヒーレンスの低いデータは、解析上の誤差を拡大してしまう恐れがあるため、測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間のコヒーレンスが予め定められたしきい値（たとえば、０．７）より小さいデータを排除することが望ましい。そのため、位相線傾き算出部（実測）３２（図２）は、測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間のコヒーレンスが予め定められたしきい値より高い周波数成分のみを用いて、位相線の傾きｇｅｘｐを算出する。
【０１３６】
図１４は、図１２に示す測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間の各周波数成分についての位相差をプロットした位相線図である。図１４（ａ）は、図１２（ａ）に対応する位相線図を示し、図１４（ｂ）は、図１２（ｂ）に対応する位相線図を示す。なお、図１４（ａ）および図１４（ｂ）では、コヒーレンスが「０．７」より低いデータを排除している。
【０１３７】
図１４（ａ）および図１４（ｂ）を参照して、それぞれの位相線図は、±１８０°を境界として不連続点を有している。これは、所定の周波数以上の周波数成分において、１周期（３６０°）以上の位相差が生じていることを意味する。そこで、位相線傾き算出部（実測）３２（図２）は、このような位相線図の不連続点に対して、１または２以上の周期に相当する単位（ｎ×３６０°）で補正を行なった上で、実測の位相差特性を算出する。
【０１３８】
図１５および図１６は、位相線傾き算出部（実測）３２が実施する位相線図の補正処理を説明するための模式図である。
【０１３９】
図１５（ａ）を参照して、測定信号Ｐａ（ｔ）を周波数変換して得られる位相特性Ｐａ（ｆ）と、測定信号Ｐｂ（ｔ）を周波数変換して得られる位相特性Ｐｂ（ｆ）とを比較し、周波数ｆｉに対応する位相差Ａｉを位相線図上にプロットする。なお、周波数ｆｉは、低周波側から数えてｉ番目の周波数成分である。そして、位相線図上にプロットされる位相差Ａｉのうち、不連続点が存在しない範囲のｎ個の位相差｛Ａ１，Ａ２，・・・，Ａｎ｝を用いて、初期回帰直線ｌ０を算出する（図１５（ｂ））。
【０１４０】
次に、ｎ＋１番目の位相差Ａｎ＋１と、初期回帰直線ｌ０の周波数ｆｎ＋１に対応する位相とが比較される。図１５（ｃ）に示すように、両者の偏差Δφ（０，ｉ＋１）が１８０°未満であれば、初期基準回帰ｌ０の算出に用いたｎ個の位相差｛Ａ１，Ａ２，・・・，Ａｎ｝にｎ＋１番目の位相差Ａｎ＋１を加えた位相差群｛Ａ１，Ａ２，・・・，Ａｎ，Ａｎ＋１｝を用いて、回帰直線ｌ１を算出する（図１５（ｄ））。
【０１４１】
一方、図１６（ｃ）に示すように、両者の偏差Δφ（０，ｉ＋１）が１８０°以上であれば、不連続点があると判断される。そして、初期回帰直線ｌ０に対する偏差が１８０°未満となるように、位相差Ａｎ＋１に３６０°×ｍ（ｍは、１以上の整数）を減算して、位相差Ａｎ＋１を位相差＃Ａｎ＋１に遷移させる。すなわち、測定されたデータが有する見かけ上の誤差を補正する。
【０１４２】
そして、初期基準回帰ｌ０の算出に用いたｎ個の位相差｛Ａ１，Ａ２，・・・，Ａｎ｝に、この補正された位相差＃Ａｎ＋１を加えた位相差群｛Ａ１，Ａ２，・・・，Ａｎ，＃Ａｎ＋１｝を用いて、回帰直線ｌ１を算出する（図１６（ｄ））。
【０１４３】
以下、同様にして、すべての位相差Ａｉについて、プロットおよび回帰直線の更新が繰返される。
【０１４４】
図１７は、図１４に示す位相線図を補正して連続化した結果を示す図である。図１７（ａ）は、図１４（ａ）に対応する位相線図を示し、図１７（ｂ）は、図１４（ｂ）に対応する位相線図を示す。
【０１４５】
図１７（ａ）および図１７（ｂ）を参照して、上述の方法によって、各位相差が補正されて位相線図が連続化されていることが分かる。また、図１７（ａ）および図１７（ｂ）には、プロットした位相データの回帰直線についてもそれぞれ示されており、これらの課回帰直線の傾きが図２に示す傾きｇｅｘｐに相当する。
【０１４６】
（フローチャート）
図１８は、この発明の実施の形態１に従う血管状態評価装置１００において実行される処理の手順を示すフローチャートである。図１８のフローチャートに示される各処理は、制御部２のＣＰＵ１０がＲＯＭ１２に予め格納されているプログラムを読出し、ＲＡＭ１３上に展開して実行することで、図２に示す各機能を実現する。
【０１４７】
図１８を参照して、ユーザによる操作部６などの操作に応答して、ＣＰＵ１０は、測定部２０ａ，２０ｂに対して測定指令を与え、測定部２０ａ，２０ｂが被験者２００の所定の測定部位における生体信号の測定を開始する（ステップＳ１００）。
【０１４８】
次に、ＣＰＵ１０は、測定部２０ａ，２０ｂで測定される時間波形である測定信号Ｐａ（ｔ），Ｐｂ（ｔ）を周波数領域の位相特性Ｐａ（ｆ），Ｐｂ（ｆ）に変換する（ステップＳ１０２）。そして、ＣＰＵ１０は、位相特性Ｐａ（ｆ）と位相特性Ｐｂ（ｆ）との間の各周波数成分についての位相差に基づいて、実測の位相差特性（傾きｇｅｘｐ）を算出する（ステップＳ１０４）。
【０１４９】
また、ＣＰＵ１０は、変数ｋを初期値ｋ０に設定する（ステップＳ１０６）。そして、ＲＯＭ１２などに格納されている循環系モデルを参照して、心臓からそれぞれ押圧カフ２４ａ，２４ｂが装着された２つの測定部位までの血管経路の伝達特性を示す２つの伝達関数を算出する（ステップＳ１０８）。この伝達関数の算出に際して、各区間のヤング率は、循環系モデル中の規定された基準のヤング率に変数ｋを乗じた値が用いられる。さらに、ＣＰＵ１０は、ステップＳ１０６で算出した２つの伝達関数の間の位相差特性（傾きｇ（ｋ））を算出する（ステップＳ１１０）。
【０１５０】
その後、ＣＰＵ１０は、ステップＳ１０４において算出された実測の位相差特性（傾きｇｅｘｐ）と、ステップＳ１１０において算出された伝達関数の間の位相差特性（傾きｇ（ｋ））との偏差Δ（＝｜ｇｅｘｐ−ｇ（ｋ）｜）を算出する（ステップＳ１１２）。そして、ＣＰＵ１０は、偏差Δが予め定められた収束条件を満足しているか否かを判断する（ステップＳ１１４）。この収束条件としては、代表的に、偏差Δが予め定められたしきい値より小さいか否かが判断される。
【０１５１】
偏差Δが予め定められた収束条件を満足していなければ（ステップＳ１１４においてＮＯ）、ＣＰＵ１０は、偏差Δを小さくする方向に、変数ｋを所定値だけ増減させる（ステップＳ１１６）。そして、ステップＳ１０８以降の処理を再度実行する。
【０１５２】
一方、偏差Δが予め定められた収束条件を満足していれば（ステップＳ１１４においてＹＥＳ）、ＣＰＵ１０は、当該時点の変数ｋの値を最適解ｋｏｐｔとして決定する（ステップＳ１１８）。そして、ＣＰＵ１０は、決定した最適解ｋｏｐｔ、最適解ｋｏｐｔを用いて換算されたヤング率、最適解ｋｏｐｔの評価結果などを、表示部４へ出力する（ステップＳ１２０）。そして、評価処理は終了する。
【０１５３】
なお、上述の説明においては、心臓を入力端とした全身についての脈波伝播モデル（伝達関数）、末梢部モデル、および伝達関数Ｇａ（ｆ），Ｇｂ（ｆ）を算出する方法について詳述したが、このようなモデルや伝達関数を必ずしも評価処理毎に算出する必要はない。すなわち、評価処理前に算出したモデルや伝達関数を記憶部３４に予め格納しておいてもよい。
【０１５４】
この発明の実施の形態１によれば、実測の測定信号の周波数特性を反映した弾力度変数ｋの最適解を得ることができる。これにより、周波数特性を考慮した上で、血管の状態（動脈硬化の度合い）をより高い精度で評価することができる。
【０１５５】
［実施の形態２］
上述のこの発明の実施の形態１では、血管の弾力度合いを示す弾力度変数を算出する構成について例示したが、従来から実用化されている脈波速度法と同様の評価基準を用いることができるように、脈波速度を算出することも有効である。そのため、この発明の実施の形態２では、予め定められた２点間の脈波速度を算出する構成について説明する。
【０１５６】
この発明の実施の形態２に従う血管状態評価装置１００＃の構成は、制御部２＃で実行される処理を除いて、図１に示すこの発明の実施の形態１に従う血管状態評価装置１００と同様であるので、詳細な説明は繰返さない
図１９は、この発明の実施の形態２に従う血管状態評価装置１００＃の制御部２＃で実行される機能を模式的に示す機能ブロック図である。
【０１５７】
図１９を参照して、制御部２＃は、図２に示すこの発明の実施の形態１に従う制御部２で実行される機能において、評価部４２に代えて、脈波速度モデル算出部４４と、平均脈波速度算出部４６と、評価部４８とを設けたものである。その他の機能については、図２と同様であるので、詳細な説明は繰返さない。
【０１５８】
脈波速度モデル算出部４４は、探索部４０によってフィッティングされた最適解ｋｏｐｔに基づいて、予め定められた２点間を伝播する脈波速度を算出するための数学的なモデルを算出する。ここで、２点間を伝播する脈波速度とは、２つの測定部位間において空間的に平均化した脈波速度を意味する。すなわち、上述した循環系モデルの各区間の脈波速度は、管径および管長に応じて変化するため、２つの測定部位間の伝播経路の形状に応じて、脈波速度は増速したり減速したりする。そこで、従来の脈波速度法との整合性をとるために、以下に示すような空間的な平均化処理を行なって、測定部位間の脈波速度（以下、「平均脈波速度」とも称す）を算出する。脈波速度モデル算出部４４は、このような平均化処理を含んだモデルを算出する。
【０１５９】
平均脈波速度算出部４６は、脈波速度モデル算出部４４で算出された数学モデルに基づいて演算を行ない、平均脈波速度Ｃａｖｅを算出する。
【０１６０】
評価部４８は、平均脈波速度算出部４６で算出された平均脈波速度Ｃａｖｅを予め定められた基準値と比較し、動脈硬化の度合いについての評価を、表示部４（図１）などに出力する。
【０１６１】
（平均脈波速度の算出）
記憶部３４に予め格納されている循環系モデル、および探索部４０によってフィッティングされた最適解ｋｏｐｔに基づいて、２つの測定部位間の伝播経路における血管の形状値は既知となっている。
【０１６２】
図２０は、２つの測定部位ＭｐＡ，ＭｐＢの間の経路を模式的に示す管路モデルの図である。
【０１６３】
図２０を参照して、２つの測定部位ＭｐＡ，ＭｐＢの間の経路に、ｎ個の要素管路（区間）が直列的に接続されているとする。区間ｉについての区間長をＬｉとし、脈波速度をＣｉとし、脈波の伝播に要する時間をｔｉとすれば、測定部位ＭｐＡ−ＭｐＢ間の平均脈波速度は、（２６）式で与えられる。
【０１６４】
【数２２】
【０１６５】
ここで、ｔｉ＝Ｌｉ／Ｃｉであることを用いれば、（２６）式は（２７）式にように表される。
【０１６６】
【数２３】
【０１６７】
なお、各区間の脈波速度Ｃｉは、（１２）式を用いて、（２８）式のように表される。
【０１６８】
【数２４】
【０１６９】
このように、各区間の脈波速度を厳密に評価することで、測定部位間の平均脈波速度について周波数特性を反映した上で正確に算出することができる。さらに、各区間の脈波速度が経路によって異なる場合においても、測定部位間の平均的な脈波速度として有意な値を保証することができる。
【０１７０】
（評価例）
図２１は、図１２に示す被験者２００ａ，２００ｂから実測した測定信号Ｐａ（ｔ），Ｐｂ（ｔ）に基づいて、平均脈波速度を算出した結果を示す図である。図２１（ａ）は、被験者２００ａの算出結果を示し、図２１（ｂ）は、被験者２００ｂの算出結果を示す。
【０１７１】
図２１（ａ）および図２１（ｂ）を参照して、被験者２００ａおよび２００ｂのそれぞれについて、平均脈波速度が周波数特性を反映した上で同定されていることが分かる。なお、被験者２００ａおよび２００ｂの算出結果を比較すると、被験者２００ｂの脈波速度がより大きく、相対的に動脈硬化が進行していることを示している。
【０１７２】
図２２は、この発明の実施の形態２に従う血管状態評価方法による算出結果を従来の脈波速度法（ｂａＰＷＶ法）によって得られた測定結果と比較した図である。なお、図２２は、２３名の被験者を対象とした結果を示す。また、上述したように、各区間の脈波速度Ｃｉは周波数依存性を有しているため、各測定体の拍動周期に対応した周波数に対応する脈波速度Ｃｉおよび平均脈波速度Ｃａｖｅを採用した。
【０１７３】
図２２を参照して、この発明の実施の形態２に従う血管状態評価方法による算出結果と、従来の脈波速度法によって得られた測定結果との間の相関係数は０．９３となった。この結果から、本実施の形態に従う血管状態評価方法は、従来の脈波速度法による測定結果と比較的高い相関を有しており、脈波速度法によって蓄積された評価基準と同様の評価基準を用いて、動脈硬化の度合いを評価することができる。
【０１７４】
（フローチャート）
図２３は、本実施の形態２に従う血管状態評価装置１００＃において実行される処理の手順を示すフローチャートである。図２３のフローチャートに示される各処理は、制御部２＃のＣＰＵ１０がＲＯＭ１２に予め格納されているプログラムを読出し、ＲＡＭ１３上に展開して実行することで、図１９に示す各機能を実現する。
【０１７５】
図２３を参照して、まず、ＣＰＵ１０は、図１８に示すステップＳ１００〜ステップＳ１１８と同様の処理を実行する。なお、これらの処理は、図１８と同様であるので、詳細な説明は繰返さない。
【０１７６】
そして、ＣＰＵ１０は、ステップＳ１１８において決定した最適解ｋｏｐｔと、ＲＯＭ１２などに格納されている循環系モデルの形状値とに基づいて、予め定められた２つの測定部位間の脈波速度モデルを算出する（ステップＳ２０２）。さらに、ＣＰＵ１０は、算出した脈波速度モデルに従って、平均脈波速度を算出する（ステップＳ２０４）。
【０１７７】
さらに、ＣＰＵ１０は、ステップＳ２０４で算出した平均脈波速度の評価結果などを、表示部４へ出力する（ステップＳ２０６）。そして、評価処理は終了する。
【０１７８】
この発明の実施の形態２によれば、実測の測定信号の周波数特性を反映した平均脈波速度を高い精度で算出できる。この平均脈波速度は、従来の脈波速度法を用いて測定した結果とも比較的高い相関値をもつので、従来から蓄積されてきた脈波速度法における評価基準を用いて、血管の状態を判定することもできる。
【０１７９】
［その他の形態］
さらに、本実施の形態に従う血管状態評価装置における評価方法を実現させるためのプログラムを提供することもできる。このようなプログラムは、コンピュータに付属するフレキシブルディスク、ＣＤ−ＲＯＭ（Compact Disk-Read Only Memory）、ＲＯＭ、ＲＡＭおよびメモリカードなどのコンピュータ読取り可能な記録媒体に記録させた上、プログラム製品として提供することもできる。あるいは、コンピュータに内蔵するハードディスクなどの記録媒体にて記録させて、プログラムとして提供することもできる。また、ネットワークを介したダウンロードによって、プログラムとして提供することもできる。
【０１８０】
なお、本発明にかかるプログラムは、コンピュータのオペレーションシステム（ＯＳ）の一部として提供されるプログラムモジュールのうち、必要なモジュールを所定の配列で所定のタイミングで呼出して処理を実行させるものであってもよい。その場合、プログラム自体には上記モジュールが含まれずＯＳと協働して処理が実行される。このようなモジュールを含まないプログラムも、本発明にかかるプログラムに含まれ得る。
【０１８１】
また、本発明にかかるプログラムは、たとえば通常の血圧測定を実行するためのプログラムといった他のプログラムの一部に組込まれて提供されるものであってもよい。その場合にも、プログラム自体には当該他のプログラムに含まれるモジュールが含まれず、他のプログラムと協働して処理が実行される。このような他のプログラムに組込まれたプログラムも、本発明にかかるプログラムに含まれ得る。
【０１８２】
提供されるプログラム製品は、ハードディスクなどのプログラム格納部にインストールされて実行される。なお、プログラム製品は、プログラム自体と、プログラムが記録された記録媒体とを含む。
【０１８３】
今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は、上記した説明ではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。
【０１８４】
（付表）
Ａｖｏｌｉｏモデルの主要な形状値を以下の付表に示す。
【０１８５】
【表２】
【０１８６】
【表３】
【０１８７】
【表４】
【０１８８】
【表５】
【０１８９】
【表６】
【０１９０】
【表７】
【図面の簡単な説明】
【０１９１】
【図１】この発明の実施の形態１に従う血管状態評価装置の概略構成図である。
【図２】この発明の実施の形態１に従う血管状態評価装置の制御部で実行される機能を模式的に示す機能ブロック図である。
【図３】血管内の血液の一次元流れモデルを示す図である。
【図４】図３に示す検査体積に働く力および出入りする運動量を示す図である。
【図５】血管を一次元線形分布定数モデル化した模式図である。
【図６】Ａｖｏｌｉｏモデルの模式図である。
【図７】末梢部モデルの模式図である。
【図８】表１に示す形状値を用いて、図７に示す細動脈終端での反射率ＳＴを変化させた場合の末梢反射率の変化を示す図である。
【図９】中動脈・小動脈・細動脈を伝播する脈波の減衰定数を数値計算によって求めた結果を示す図である。
【図１０】図７に示す中動脈の始端（座標ｘ）から所定の周波数で加振（加圧）した場合の血管内圧力分布を示す図である。
【図１１】一様な管路における脈波伝播の様子を示す模式図である。
【図１２】被験者の上腕および足関節に押圧カフを装着して測定した圧力の時間波形を示す図である。
【図１３】図１２に示す上腕の圧力波形と足関節の圧力波形との間のコヒーレンスを示す図である。
【図１４】図１２に示す測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間の各周波数成分についての位相差をプロットした位相線図である。
【図１５】位相線傾き算出部（実測）が実施する位相線図の補正処理を説明するための模式図である。
【図１６】位相線傾き算出部（実測）が実施する位相線図の補正処理を説明するための模式図である。
【図１７】図１４に示す位相線図を補正して連続化した結果を示す図である。
【図１８】この発明の実施の形態１に従う血管状態評価装置において実行される処理の手順を示すフローチャートである。
【図１９】この発明の実施の形態２に従う血管状態評価装置の制御部で実行される機能を模式的に示す機能ブロック図である。
【図２０】２つの測定部位ＭｐＡ，ＭｐＢの間の経路を模式的に示す管路モデルの図である。
【図２１】図１２に示す被験者から実測した測定信号Ｐａ（ｔ），Ｐｂ（ｔ）に基づいて、平均脈波速度を算出した結果を示す図である。
【図２２】この発明の実施の形態２に従う血管状態評価方法による算出結果を従来の脈波速度法（ｂａＰＷＶ法）によって得られた測定結果と比較した図である。
【図２３】本実施の形態２に従う血管状態評価装置において実行される処理の手順を示すフローチャートである。
【符号の説明】
【０１９２】
２，２＃ 制御部、４ 表示部、６ 操作部、１０ ＣＰＵ、１２ ＲＯＭ、１４ ＲＡＭ、２０ａ，２０ｂ 測定部、２２ａ，２２ｂ，２７ａ，２７ｂ 配管、、２４ａ，２４ｂ 押圧カフ、２５ａ，２５ｂ 圧力ポンプ、２６ａ，２６ｂ 調圧弁、２８ａ，２８ｂ 圧力センサ、３０ａ，３０ｂ 周波数変換部（ＦＦＴ）、３２ 位相線傾き算出部（実測）、３４ 記憶部、３６ 伝達関数算出部、３８ 位相線傾き算出部（モデル）、４０ 探索部、４２，４８ 評価部、４４ 脈波速度モデル算出部、４６ 平均脈波速度算出部、５０ 検査体積、１００，１００＃ 血管状態評価装置（評価装置）。
【技術分野】
【０００１】
この発明は生体を構成する血管の状態を評価する血管状態評価装置、血管状態評価方法および血管状態評価プログラムに関し、特に血管経路を伝達関数としてモデル化した上で、血管の動脈硬化の度合いを評価する技術に関するものである。
【背景技術】
【０００２】
近年、動脈硬化に起因する循環器系疾患が増加しており、これを受けて血管の動脈硬化の度合いを評価する評価装置が実用化されている。動脈硬化の度合いを評価する代表的な方法として、脈波速度法が知られている。この脈波速度法は、心臓の拍動に伴う血圧変化が血管上を伝播する速度（脈波速度）と血管の弾力度合い（剛性）との間に相関があることを利用したものである。すなわち、弾性管である血管内を脈波が進む際、管壁が硬く、内径が細く、管厚が厚いほど脈波速度は増大するため、脈波速度を測定することで血管の動脈硬化の度合いを知ることができる。特に、両上腕および両足首における血圧の時間波形を用いる、ｂａＰＷＶ法（brachial-ankle Pulse Wave Velocity method）による評価装置が実用化されている。
【０００３】
このような脈波速度の測定方法として、たとえば、特開２００６−３２６３３４号公報（特許文献１）には、複数の電圧波形のうち、隣り合う電圧電極対から得られた電圧波形間の時間的なひずみを検出する検出手段と、複数の電圧波形の全てについて、隣り合う電圧電極間の距離及び／又は時間的なずれとを用い、隣り合う電圧電極間における脈波伝播速度又は脈波伝播時間の変化率を求める算出手段とを有する脈波伝播速度測定装置が開示されている。
【特許文献１】特開２００６−３２６３３４号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
特開２００６−３２６３３４号公報（特許文献１）に記載されるように、脈波速度を測定する方法として、血管経路上の複数の点で測定した時間波形の間の時間差（遅れ時間）を求め、各点の心臓からの経路差を当該時間差で除することで脈波速度が算出される。
【０００５】
しかしながら、実際の脈波速度は伝播経路や周波数に依存するため、単純に経路差を時間差で除するだけでは、厳密な脈波速度を算出することはできない。すなわち、被評価者の血管径や血管長、および脈波に含まれる周波数成分などに依存して、本来の脈波速度からずれた値が算出される場合があった。そのため、動脈硬化の度合いの評価精度をあまり高められないという問題があった。
【０００６】
そこで、この発明は、かかる問題を解決するためになされたものであり、その目的は、動脈硬化の度合いをより高い精度で評価することが可能な血管状態評価装置、血管状態評価方法および血管状態評価プログラムを提供することである。
【課題を解決するための手段】
【０００７】
この発明のある局面に従う血管状態評価装置は、生体を構成する血管を複数の区間に分割してモデル化した循環系モデルを格納する記憶手段を備え、循環系モデルは、複数の区間の各々を代表する形状値を含む。さらに、血管状態評価装置は、生体の第１測定部位に装着されて第１生体信号の時間波形を測定する第１測定手段と、生体の第２測定部位に装着され、第１測定手段と同期して第２生体信号の時間波形を測定する第２測定手段と、第１生体信号と第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出する第１算出手段と、第１測定部位までの血管経路に対応して循環系モデルに基づいて規定された第１伝達関数と、第２測定部位までの血管経路に対応して循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出する第２算出手段とを備え、第１伝達関数および第２伝達関数は、血管の弾力度合いを示す弾力度変数を含む。さらに、血管状態評価装置は、第１算出手段によって算出される実測の位相差特性に基づいて、第２算出手段によって算出される位相差特性をフィッティングすることにより、弾力度変数を決定する探索手段を備える。
【０００８】
この発明によれば、循環系モデルに基づいて、第１および第２測定部位までの血管経路に対応してそれぞれ第１および第２伝達関数が規定された上で、第１伝達関数と第２伝達関数との間の位相差特性を算出する。ここで、第１および第２伝達関数は、血管の弾力度合いを示す弾力度変数を含む。一方で、第１および第２測定部位からそれぞれ第１および第２生体信号の時間波形が測定された上で、第１生体信号と第２生体信号との間の各周波数成分についての位相差に基づいて実測の位相差特性を算出する。さらに、実測の位相差特性に基づいて、第１伝達関数と第２伝達関数との間の位相差特性をフィッティングすることにより、弾力度変数を決定する。
【０００９】
ここで、第１および第２伝達関数は、脈波の周波数特性を含んで規定されるので、フィッティングによって、実測の測定信号の周波数特性を反映した弾力度変数の最適解が得られることになる。よって、周波数特性を考慮した上で、血管の状態（動脈硬化の度合い）をより高い精度で評価することができる。
【００１０】
好ましくは、血管状態評価装置は、第１および第２測定部位までの血管経路にそれぞれに対応する各区間の形状値に基づいて、第１および第２伝達関数を算出する伝達関数算出手段をさらに備える。
【００１１】
好ましくは、伝達関数算出手段は、各区間に対応付けられた、血管の圧力および血液流量を入力変数とする分布定数モデルを用いて、第１および第２伝達関数を算出し、分布定数モデルの各々は、対応する区間における血液の流れやすさに応じた縦インピーダンスと、弾力度変数を有する横インピーダンスとを含む。
【００１２】
好ましくは、血管状態評価装置は、探索手段によってフィッティングされた弾力度変数に基づいて、血管における脈波速度を算出する脈波速度算出手段をさらに備える。
【００１３】
さらに好ましくは、脈波速度算出手段は、第１測定部位までの血管経路に対応する各区間の形状値および第２測定部位までの血管経路に対応する各区間の形状値に基づいて脈波速度を算出する。
【００１４】
好ましくは、循環系モデルは、形状値として、血管径および血管長を含む。
好ましくは、循環系モデルは、生体を構成する血管を複数の区分に分類した上で、複数の区分の少なくとも１つの区分に属する血管をモデル化したものである。
【００１５】
さらに好ましくは、生体を構成する血管は、血管径の大きさに基づいて複数の区分に分類される。
【００１６】
また好ましくは、伝達モデル算出手段は、各区間に含まれる血管のうち、循環系モデルにおいてモデル化されていない血管をモデル化した末梢部モデルを、各区間に対応する循環系モデルに付加した上で、伝達関数を算出する。
【００１７】
さらに好ましくは、伝達モデル算出手段は、各区間の循環系モデルを血管の形状差に基づいて変換することで、当該区間の末梢部モデルを算出する。
【００１８】
さらに好ましくは、伝達モデル算出手段は、末梢部モデルの終端を無反射条件の下で伝達関数を算出する。
【００１９】
好ましくは、血管状態評価装置は、第１生体信号から各周波数成分についての位相を示す第１位相特性を算出する第１周波数変換手段と、第２生体信号から各周波数成分についての位相を示す第２位相特性を算出する第２周波数変換手段とをさらに備え、第１算出手段は、第１位相データと第２位相データとを差分して差分位相データを算出する。さらに、血管状態評価装置は、差分位相データにおける周期遅れに起因する位相ずれを、１または２以上の周期に相当する位相の単位で補正することで、実測の位相差特性を算出する。
【００２０】
好ましくは、第１算出手段は、第１生体信号と第２生体信号との間のコヒーレンス値が予め定められたしきい値より高い周波数成分を用いて実測の位相差特性を算出する。
【００２１】
この発明の別の局面に従えば、生体を構成する血管を複数の区間に分割してモデル化した循環系モデルを用いて、生体を構成する血管の状態を評価する血管状態評価方法であって、循環系モデルは、複数の区間の各々を代表する形状値を含む。血管状態評価方法は、生体の第１測定部位から第１生体信号の時間波形を測定するとともに、生体の第２測定部位から第２生体信号の時間波形を測定するステップと、第１生体信号と第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するステップと、第１測定部位までの血管経路に対応して循環系モデルに基づいて規定された第１伝達関数と、第２測定部位までの血管経路に対応して循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出するステップとを備え、第１伝達関数および第２伝達関数は、血管の弾力度合いを示す弾力度変数を含む。さらに、血管状態評価方法は、実測の位相差特性に基づいて、第１伝達関数と第２伝達関数との間の位相差特性をフィッティングすることにより、弾力度変数を決定するステップを備える。
【００２２】
この発明のさらに別の局面に従えば、演算処理部を有するコンピュータに、生体を構成する血管を複数の区間に分割して予めモデル化した循環系モデルを用いて、生体を構成する血管の状態の評価を実行させるための血管状態評価プログラムであって、循環系モデルは、複数の区間の各々を代表する形状値を含む。血管状態評価プログラムは、演算処理部が、生体の第１測定部位で測定される第１生体信号の時間波形を取得するとともに、生体の第２測定部位で測定される第２生体信号の時間波形を取得するステップと、演算処理部が、第１生体信号と第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するステップと、演算処理部が、第１測定部位までの血管経路に対応して循環系モデルに基づいて規定された第１伝達関数と、第２測定部位までの血管経路に対応して循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出するステップとを備え、第１伝達関数および第２伝達関数は、血管の弾力度合いを示す弾力度変数を含む。さらに、血管状態評価プログラムは、演算処理部が、演算処理部が、実測の位相差特性に基づいて、第１伝達関数と第２伝達関数との間の位相差特性をフィッティングすることにより、弾力度変数を決定するステップを備える。
【発明の効果】
【００２３】
この発明によれば、動脈硬化の度合いをより高い精度で評価することが可能な血管状態評価装置、血管状態評価方法および血管状態評価プログラムを実現できる。
【発明を実施するための最良の形態】
【００２４】
この発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。
【００２５】
［実施の形態１］
（装置構成）
図１は、この発明の実施の形態１に従う血管状態評価装置１００の概略構成図である。
【００２６】
図１を参照して、この発明の実施の形態１に従う血管状態評価装置（以下、単に「評価装置」とも記す）１００は、制御部２と、表示部４と、操作部６と、測定部２０ａ，２０ｂとを含む。
【００２７】
制御部２は、評価装置１００全体の制御を行なう装置であり、代表的に、ＣＰＵ（Central Processing Unit）１０と、ＲＯＭ（Read Only Memory）１２と、ＲＡＭ（Random Access Memory）１４とを含むコンピュータで構成される。
【００２８】
ＣＰＵ１０は、演算処理部に相当し、ＲＯＭ１２に予め格納されているプログラムを読出して、ＲＡＭ１４をワークメモリとして使用しながら、当該プログラムを実行する。また、ＲＯＭ１２には、少なくとも後述する循環系モデルが予め格納されており、ＣＰＵ１０は、本実施の形態に従う血管状態評価方法を格納したプログラムの実行に際して、当該循環系モデルを参照する。
【００２９】
また、制御部２には、表示部４および操作部６が接続されている。表示部４は、ユーザによる各種設定の入力を促したり、制御部２からの演算結果を表示したりする。これに対して、ユーザは、表示部４に表示される内容を確認しながら操作部６を操作して、所望の設定入力や操作を行なう。なお、表示部４は、一例として、ＬＥＤ（Light Emitting Diode）やＬＣＤ（Liquid Crystal Display）などからなる。
【００３０】
より具体的には、制御部２は、測定部２０ａ，２０ｂに対して測定指令を与えるとともに、当該測定指令に応答して測定された測定信号Ｐａ（ｔ），Ｐｂ（ｔ）を受信し、当該測定信号Ｐａ（ｔ），Ｐｂ（ｔ）に基づいて、本実施の形態に従う血管状態評価方法を実行する。
【００３１】
測定部２０ａ，２０ｂは、被験者２００の所定の測定部位に装着された押圧カフ（空気袋）２４ａ，２４ｂの内圧（以下、「カフ圧」という）を加圧して、それぞれの測定部位における生体信号（一例として、脈波）の時間波形を測定する。なお、後述するように、制御部２は、測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するので、制御部２からは、測定部２０ａおよび２０ｂが互いに同期して生体信号を測定できるように、測定指令が同時に与えられる。
【００３２】
より詳細には、たとえば、押圧カフ２４ａおよび２４ｂは、それぞれ被験者２００の足首部および上腕部に装着され、それぞれ配管２２ａおよび２２ｂを介して測定部２０ａおよび２０ｂから供給される空気によって加圧される。そして、この加圧によって、押圧カフ２４ａおよび２４ｂは対応の測定部位に押圧され、当該測定部位の脈波に応じた圧力変化がそれぞれ配管２２ａおよび２２ｂを介して測定部２０ａおよび２０ｂへ伝達される。測定部２０ａ，２０ｂは、この伝達される圧力変化を検出すことで、測定部位の脈波の時間波形を測定する。なお、測定信号Ｐａ（ｔ）およびＰｂ（ｔ）の所定の周波数成分（一例として、０〜２０［Ｈｚ］）に対して演算処理を行なうことが好ましいので、測定信号Ｐａ（ｔ）およびＰｂ（ｔ）の測定周期（サンプリング周期）は、この周波数成分に応じた時間間隔（一例として、２５ｍｓｅｃ）より短くすることが好ましい。
【００３３】
このような測定動作を実行するために、測定部２０ａは、圧力センサ２８ａと、調圧弁２６ａと、圧力ポンプ２５ａと、配管２７ａとを含む。圧力センサ２８ａは、配管２２ａを介して伝達される圧力変動を検出するための検出部位であり、一例として、単結晶シリコンなどからなる半導体チップに所定間隔に配列された複数のセンサエレメントを含む。調圧弁２６ａは、圧力ポンプ２５ａと押圧カフ２４ａとの間に介挿され、測定時に押圧カフ２４ａを加圧に用いられる圧力を所定の範囲に維持する。圧力ポンプ２５ａは、制御部２からの測定指令に応じて作動し、押圧カフ２４ａを加圧するための加圧空気を供給する。
【００３４】
同様に、測定部２０ｂは、圧力センサ２８ｂと、調圧弁２６ｂと、圧力ポンプ２５ｂと、配管２７ｂとを含む。各部の構成については、測定部２０ａと同様であるので、詳細な説明は繰返さない。
【００３５】
なお、本実施の形態では、生体信号の一例として、脈波によって生じる圧力変化を、圧力カフを用いて測定する構成について説明するが、たとえば、被験者２００の測定部位に微少の一定電流を流すとともに、脈波の伝播に応じて生じるインピーダンス（生体インピーダンス）の変化によって生じる電圧変化を測定するようにしてもよい。
【００３６】
なお、図１に示す評価装置１００と本願発明との対応関係については、測定部２０ａ，２０ｂと、配管２２ａ，２２ｂと、押圧カフ２４ａ，２４ｂとが「第１および第２の測定手段」に相当する。
【００３７】
（機能ブロック図）
制御部２は、予め格納する循環系モデルに基づいて、それぞれ押圧カフ２４ａおよび２４ｂが装着された測定部位までの血管経路に対応して規定される２つの伝達関数を算出する。ここで、それぞれの伝達関数は、血管の弾力度合いを示す弾力度変数を含む。すなわち、弾力度変数は、血管の動脈硬化の度合いを示す指標である。本実施の形態では、弾力度変数の代表例として「ヤング率」を用いるが、血管の剛性や柔軟度を示す他の変数を用いてもよい。そして、制御部２は、測定信号Ｐａ（ｔ）およびＰｂ（ｔ）を周波数領域の信号に変換した上で、両者間の実測の位相差特性を算出し、この実測の位相差特性が上述した２つの伝達関数の間の位相差特性と一致するように、弾力度変数をフィッティング（同定）する。このフィッティングされた弾力度変数が被験者２００の動脈硬化の度合いを示す値となる。以下、このような制御部２における処理動作を実現するための機能ブロックについて説明する。
【００３８】
図２は、この発明の実施の形態１に従う血管状態評価装置１００の制御部２で実行される機能を模式的に示す機能ブロック図である。
【００３９】
図２を参照して、制御部２は、周波数変換部（ＦＦＴ）３０ａ，３０ｂと、位相線傾き算出部（実測）３２と、記憶部３４と、伝達関数算出部３６と、位相線傾き算出部（モデル）３８と、探索部４０と、評価部４２とを実現する。
【００４０】
周波数変換部３０ａおよび３０ｂは、それぞれ時間波形である測定信号Ｐａ（ｔ）およびＰｂ（ｔ）を所定期間にわたって受信し、当該受信した測定信号Ｐａ（ｔ）およびＰｂ（ｔ）を周波数領域の関数に変換する周波数変換部である。代表的に、周波数変換部３０ａおよび３０ｂは、高速フーリエ変換（ＦＦＴ：Fast Fourier Transformer）を用いて、周波数変換を実行する。なお、高速フーリエ変換に限らず、時間領域の関数をフーリエ級数などの周波数領域の関数に変換するものであれば、いずれのロジックを用いてもよい。
【００４１】
そして、周波数変換部３０ａは、測定信号Ｐａ（ｔ）の各周波数成分についての位相を示す位相特性Ｐａ（ｆ）を算出し、算出した位相特性Ｐａ（ｆ）を位相線傾き算出部（実測）３２へ出力する。同様に、周波数変換部３０ｂは、測定信号Ｐｂ（ｔ）の各周波数成分についての位相を示す位相特性Ｐｂ（ｆ）を算出し、算出した位相特性Ｐｂ（ｆ）を位相線傾き算出部（実測）３２へ出力する。
【００４２】
位相線傾き算出部（実測）３２は、ユーザによる操作部６（図１）などの操作に応答して、測定指令を測定部２０ａ，２０ｂへ与える。この測定指令を与えた後、周波数変換部３０ａおよび３０ｂから出力される位相特性Ｐａ（ｆ）および位相特性Ｐｂ（ｆ）を受信し、両者間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出する。具体的には、位相線傾き算出部（実測）３２は、周波数成分毎に位相特性Ｐａ（ｆ）およびＰｂ（ｆ）の値を比較し、両者の位相差を算出する。後述するように、このように算出された位相差は、周波数成分についての一次関数として近似できるので、この近似された一次関数（位相線）の傾きｇｅｘｐ［ｄｅｇ／Ｈｚ］を、実測の位相差特性として探索部４０へ出力する。すなわち、偏角φｅｘｐ＝∠位相特性Ｐａ（ｆ）／位相特性Ｐｂ（ｆ）として算出される偏角φを用いて、傾きｇｅｘｐ＝ｔａｎ（φｅｘｐ）として定義できる。
【００４３】
一方、伝達関数算出部３６は、心臓からそれぞれ押圧カフ２４ａおよび２４ｂが装着された２つの測定部位までの血管経路の伝達特性を示す２つの伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）を算出し、位相線傾き算出部（モデル）３８へ出力する。より具体的には、伝達関数算出部３６は、記憶部３４に予め格納されている循環系モデルに基づいて、心臓を入力端とした全身についての脈波伝播モデル（伝達関数）を算出した上で、全身の脈波伝播モデルの中で２つの測定部位までの血管経路に対応する伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）を算出する。この際、伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）には、変数ｋを含む形でヤング率が組込まれるが、変数ｋには探索部４０によって具体な値が設定される。
【００４４】
記憶部３４は、被験者２００を構成する血管を複数の区間に分割してモデル化した循環系モデルを格納している。この循環系モデルは、各区間に対応付けて、各区間を代表する形状値が規定されている。このような形状値の一例として、本実施の形態では、各区間の血管径、血管長、血管壁の厚さなどが含まれる。なお、循環系モデルについては、後に詳述する。
【００４５】
位相線傾き算出部（モデル）３８は、伝達関数Ｇａ（ｆ）と伝達関数Ｇｂ（ｆ）との間の位相差特性を算出し、算出した位相差特性を探索部４０へ出力する。具体的には、位相線傾き算出部（モデル）３８は、周波数領域における位相特性Ｐａ（ｆ）と位相特性Ｐｂ（ｆ）との間の位相差である位相線の傾きｇ（ｋ）［ｄｅｇ／Ｈｚ］を、位相差特性として探索部４０へ出力する。ここで、傾きｇ（ｋ）は、偏角φｍｏｄｅｌ＝∠伝達関数Ｇａ（ｆ）／伝達関数Ｇｂ（ｆ）として算出される偏角φｍｏｄｅｌを用いて、傾きｇ（ｋ）＝ｔａｎ（φｍｏｄｅｌ）として定義できる。
【００４６】
探索部４０は、位相線傾き算出部（実測）３２によって算出される傾きｇｅｘｐに基づいて、位相線傾き算出部（モデル）３８によって算出される傾きｇ（ｋ）をフィッティングすることにより、変数ｋを決定する。すなわち、傾きｇ（ｋ）と傾きｇｅｘｐとが互いに略一致するまで、変数ｋが初期値ｋ０から順次変更され、伝達関数算出部３６および位相線傾き算出部（モデル）３８における演算処理が繰返し実行される。このようにして、傾きｇ（ｋ）と傾きｇｅｘｐとが互いに略一致する変数ｋｏｐｔ（最適解）が決定されると、探索部４０は、決定した変数ｋの値を評価部４２へ出力する。この決定された変数ｋの最適解が、被験者２００の動脈硬化の度合いを示す指標となる。
【００４７】
評価部４２は、探索部４０で決定された最適解ｋｏｐｔ（または、最適解ｋｏｐｔを用いて換算されたヤング率）を予め定められた基準値と比較し、動脈硬化の度合いについての評価を、表示部４（図１）などに出力する。
【００４８】
以下、このような主要機能の動作や構成について詳述する。
（物理モデル）
上述したように、伝達関数算出部３６は、心臓を入力端（始点）とする血管経路の伝達特性を示す伝達関数を算出するが、この算出される伝達関数は、脈波が血管を伝播する力学的なモデルから解析的に算出される。本実施の形態では、血管の各区間を一次元線形分布定数モデル化して、伝達関数を算出する構成について例示する。
【００４９】
まず、血管を軸対象の微小変形する薄肉円管とし、内部の血液の流れを非粘性流体の層流とし、かつ反射波はないと仮定してモデル化すると、脈波速度Ｃｐと血管壁のヤング率Ｅとの関係は、Moens-Kortewegの式と呼ばれる（１）式で表される。なお、脈波速度Ｃｐは、心臓の拍動に伴う血圧変化が血管上を伝播する速度を意味する。
【００５０】
【数１】
【００５１】
この（１）式から、血管が硬く、内腔が細く、または血管壁が厚いほど脈波速度Ｃｐが増加することがわかる。
【００５２】
図３は、血管内の血液の一次元流れモデルを示す図である。
一般に、血管に比べて血液の体積弾性率は十分高いので、血管を弾性円管とし、血液を非圧縮性流体として考えることができる。このような弾性管内における一次元流れの支配方程式は、以下のように導かれる。
【００５３】
図３を参照して、一次元流れモデルの断面ＣＳ１−ＣＳ２間の検査体積５０に関する質量の保存について考える。断面ＣＳ１の内腔の面積をＡ（＝πｒｉ２），流体（血液）の密度をρ，圧力をｐ，断面平均流速をＵとし、断面ＣＳ１−ＣＳ２間にある分枝血管へ単位時間に漏れ出す流体の体積を単位長さおよび単位圧力あたりＧとすると，質量保存則より（２）式が成り立つ。ここで、非圧縮性流体では密度ρは一定であるので（２）式は（３）式のように簡略化できる。
【００５４】
【数２】
【００５５】
図４は、図３に示す検査体積５０に働く力５２および出入りする運動量５４を示す図である。
【００５６】
図４を参照して、検査体積５０内における運動量５４の単位時間当たりの変化は、流入する正味の運動量５４と検査体積５０に及ぼされる力５２に等しいので，高次の微少項を省略して（４）式が導かれる。
【００５７】
【数３】
【００５８】
（４）式を連続の式を用いて整理すると（５）式に示す運動方程式が得られる．
【００５９】
【数４】
【００６０】
次に、血管を一次元線形分布定数モデル化するために、（３）式および（５）式において非線形項を省略し、変数を圧力ｐと体積流量ｑ（＝ＡＵ）に置き換えることにより、（６）式および（７）式が得られる。
【００６１】
【数５】
【００６２】
ここで、（６）式および（７）式の４個の係数の物理的な意味については、Ｒは血液が流れるときの粘性抵抗を示し、Ｌは流れが変化するときに急な変化を妨げようとする血液の慣性を示し、Ｇは血管外もしくは分岐管に流れ出す血液の流れやすさを示し、Ｃは圧力変化に応じて血管が伸縮する際に血管内に血液を蓄える能力を示す。
【００６３】
図５は、血管を一次元線形分布定数モデル化した模式図である。図５（ａ）は、（６）式および（７）式と血管の物理モデルとを対応付ける図である。図５（ｂ）は、図５（ａ）に示す物理モデルを電気的な等価回路に置き換えた図である。
【００６４】
すなわち、（６）式および（７）式は、図５（ａ）に示すような物理モデルと対応付けることができる。さらに、（６）式および（７）式において、圧力ｐを電圧ｖに置き換え、流量ｖを電流ｉに置き換えることで、図５（ｂ）に示すような電気的な等価回路（分布定数回路）に置き換えることができる。ここで、Ｒは抵抗を示し、Ｌはインダクタンスを示し、Ｇはアドミタンスを示し、Ｃはキャパシタンスを示す。
【００６５】
ここで、（６）式についてみれば、血管系では運動方程式に対応する一方、電気系ではオームの法則に対応する。そして、血管系において、断面ＣＳ１と断面ＣＳ２との間の圧力勾配によって流体が加速される現象が、電気系において、インダクタンスの両端に印加される電位差が電流を引き起こす現象に対応することを意味している。
【００６６】
また、（７）式についてみれば、血管系では連続の式（質量保存則）に対応する一方、電気系では電荷の保存則に対応する。そして、血管系において、断面ＣＳ１から断面ＣＳ２へ進めない質量の滞留分が血管を押し広げて圧力の上昇を引き起こす現象が、キャパシタに溜まった電荷が電圧の上昇を引き起こす現象に対応することを意味している。
【００６７】
さらに、（６）式および（７）式において、ｐ＝Ｐｅｊｗｔとし、ｑ＝Ｑｅｊｗｔとすると、それぞれ（８）式および（９）式に示す関係式が導出される。
【００６８】
【数６】
【００６９】
本明細書では、図５（ｂ）および（８）式に示すＺｌ（＝ｒ＋ｊωＬ）を「縦インピーダンス」と称し、図５（ｂ）および（９）式に示すＺｔ（＝（Ｇ＋ｊωＣ）−１）を「横インピーダンス」と称す。（８）式および（９）式の一般解は、ｘ＝０における圧力の進行波の振幅値をＰｆとし、後退波の振幅値をＰｒとして、それぞれ（１０）式および（１１）式となる。なお、角周波数ωと周波数ｆとの間には、ω＝２πｆの関係が成立する。
【００７０】
【数７】
【００７１】
また、伝播定数γは、減衰定数βと位相速度（脈波速度）Ｃｐとを用いて、（１２）式のように表される。
【００７２】
【数８】
【００７３】
ここで、位相速度Ｃｐは脈波が単位時間に進む距離を示す量であり、減衰定数βは脈波の振幅が単位距離を進むごとにｅ−β倍になることを示す。また、特性インピーダンスＺ０は（１３）式のように表せ、単位体積の脈波を進行方向へ進めるために必要な圧力を示す。
【００７４】
【数９】
【００７５】
さらに、距離ｌｓｅを隔てた二点での圧力Ｐｓ，Ｐｅと体積流量Ｑｓ，Ｑｅは、（１４）式の伝達行列で接続される。
【００７６】
【数１０】
【００７７】
本実施の形態では、（１４）式に示す伝達行列を血管の各区間に対応付けて算出するとともに、対象とする血管経路に応じて、各区間に対応する伝達行列を縦続接続することで、伝達関数を算出する。このとき、任意の境界より下流の条件はその境界での圧力Ｐｘ，体積速度Ｑｘとの比である（１５）式のインピーダンスＺｘで表される。
【００７８】
【数１１】
【００７９】
また、進行波と後退波の振幅の比である反射率Ｓｐは、（１６）式で表される。
【００８０】
【数１２】
【００８１】
（縦インピーダンスの算出）
縦インピーダンスＺｌは、流体の粘性抵抗と慣性の項とからなり、血管断面内の流速分布をモデル化することで求められる。
【００８２】
本実施の形態では、ウマースリーモデルに基づいて、縦インピーダンスを算出する。このウマースリーモデルは、ニュートン流体の円管内脈波流が十分に発達した状態での流速分布を表したものである。このウマースリーモデルに基づく縦インピーダンスは、第１種ベッセル関数Ｊｎを用いて（１７）式で表せる。
【００８３】
【数１３】
【００８４】
ここで、（１７）式中のαは「ウマースリーのアルファ」と称され、脈波流の粘性項と慣性項の比を示す量であり、定常流におけるレイノルズ数に相当する。また、血液の密度ρは、代表的に１．０３×１０３［ｋｇ／ｍ３］とし、血液の粘性係数μは、代表的に４×１０−３［Ｐａ・ｓ］とした。
【００８５】
なお、（１７）式に示すウマースリーモデルに代えて、非粘性モデルを用いてもよい。このモデルは、血液を非粘性流体とし、断面内流速を一定とするものである。この非粘性モデルに基づく縦インピーダンスは、（１８）式で表せる。
【００８６】
【数１４】
【００８７】
さらに、上記のモデルに代えてポアゼイユモデルを用いてもよい。このモデルは、ニュートン流体の円管内定常流が十分に発達した状態での流速分布を表したものである。このポアゼイユモデルに基づく縦インピーダンスは、（１９）式で表せる。
【００８８】
【数１５】
【００８９】
（横インピーダンスの算出）
横インピーダンスは、遺漏または分岐の項Ｇと、管のコンプライアンスの項Ｃとからなる。
【００９０】
遺漏または分岐の項については、血管壁から周囲組織への遺漏も分岐もない場合には、Ｇ＝０とする。これに対して、分岐がある場合には、分岐管のアドミタンスをＧとする。
【００９１】
次に、管のコンプライアンスの項については、厚肉円管をモデル化したコンプライアンスを用いることができる。外圧および軸方向ひずみ一定の条件における厚肉円管の軸対称微小変形のコンプライアンスは、（２０）式で表される。
【００９２】
【数１６】
【００９３】
ここで、血管壁のポアソン比νは、代表的に０．５とした。
なお、（２０）式に示す厚肉円管をモデル化したコンプライアンスに代えて、薄肉円管をモデル化したコンプライアンスを用いてもよい。外圧および軸方向ひずみ一定の条件における薄肉円管の軸対称微小変形のコンプライアンスは、（２１）式で表される。
【００９４】
【数１７】
【００９５】
（循環系モデル）
この発明の実施の形態１に従う血管状態評価装置１００で用いる循環系モデルは、生体を構成する血管を複数の区間に分割してモデル化したものである。このような循環系モデルの代表的なものとして、参考文献１「Avolio,A.P, Multi-branched Model of Human Arterial System, 1980, Med. & Biol. Engng. & Comp., 18,796」に記載されている、いわゆる「Ａｖｏｌｉｏモデル」が知られており、本実施の形態においても、循環系モデルとして、このＡｖｏｌｉｏモデルを採用する。
【００９６】
図６は、Ａｖｏｌｉｏモデルの模式図である。
図６を参照して、Ａｖｏｌｉｏモデルでは、全身の動脈を１２８の血管要素（区間）に分割し、各区間を代表する形状値を規定している。Ａｖｏｌｉｏモデルの各区間に対応付けられた形状値の一部を付表に示す。このＡｖｏｌｉｏモデルは、形状値として、各区間に対応付けられた、長さ・半径・管壁の厚さ・ヤング率を含む。なお、Ａｖｏｌｉｏモデル中のヤング率は一応の基準値であり、後述するフィッティング処理においては、この基準値に変数ｋを乗じた値が用いられる。
【００９７】
この循環系モデルは、生体を構成するさまざまな血管を複数の区分に分類した上で、当該複数の区分の少なくとも１つの区分に属する血管をモデル化したものである。代表的に、血管は、その血管径の大きさに基づいて、血管径の大きいものから順に大動脈、中動脈（φ３．２ｍｍ以上）、小動脈（φ０．５ｍｍ以上）、細動脈（φ０．０３ｍｍ以上）、毛細血管などに区分される。そして、Ａｖｏｌｉｏモデルは、これらの区分のうち、大動脈および中動脈に区分される血管についてモデル化したものである。
【００９８】
なお、血管の区分の方法については、血管径の大きさに限らず、別の指標に基づいて区分してもよい。
【００９９】
伝達関数算出部３６（図２）は、記憶部３４（図２）に予め格納されているこのような循環系モデルを参照して、（１７）式および（２０）式に従って、各区間の縦インピーダンスおよび横インピーダンスを算出する。さらに、伝達関数算出部３６は、算出した縦インピーダンスおよび横インピーダンスを用いて、（１２）式，（１３）式，（１４）式に従って各区間の伝達行列を算出し、各区間の実際の接続関係に対応付けて、各伝達関数を縦続接続および並列接続することで、心臓を基準点とした全身についての脈波伝播モデル（伝達関数）を算出する。より具体的には、（１４）式に示す２行×２列の伝達行列が、各区間の接続関係（連続、分岐、および終端など）に応じて、順次接続されていく。
【０１００】
そして、伝達関数算出部３６は、算出された全身についての脈波伝播モデル（伝達関数）の中で、２つの測定部位までの血管経路に対応する伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）を算出する。また、終端部では、その反射率に応じて（１５）式の制約が加えられる。
【０１０１】
なお、心臓から吐出される圧力（圧力Ｐｓ）および体積流量（体積流量Ｑｓ）は未知数となるが、本実施の形態では、伝達関数Ｇａ（ｆ）と伝達関数Ｇｂ（ｆ）との間の位相差特性を算出すれば十分であるので、このような未知数が存在していても、互いにキャンセルできる。
【０１０２】
（末梢部モデル）
このような、全身についての脈波伝播モデル（伝達関数）ならびに伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）の算出に際して、上述のＡｖｏｌｉｏモデルに末梢部モデルを加えることが好ましい。これは、Ａｖｏｌｉｏモデルは、比較的太い血管（大動脈および中動脈）について厳密に形状値を与えている一方で、終端条件については末梢血管を模擬した一定の反射率を規定するにとどまっているからである。そのため、より高い評価精度を得るためには、Ａｖｏｌｉｏモデルにおいてモデル化されていない血管（小動脈、細動脈、毛細血管）を考慮することが好ましい。このようなＡｖｏｌｉｏモデルにおいてモデル化されていない血管をモデル化したモデル（以下、「末梢部モデル」と称す）を、Ａｖｏｌｉｏモデルから算出される伝達行列に付加した上で、全身についての脈波伝播モデル（伝達関数）などを算出する構成について説明する。
【０１０３】
このような末梢部モデルは、末梢血管の形状値と、当該末梢血管の上流側に接続されている区間の形状値との間の形状差を利用して算出する。本実施例では、形状差として、代表的に各血管の総断面積の差を用いる。
【０１０４】
図７は、末梢部モデルの模式図である。
図７を参照して、各血管（動脈）が循環系全体で有する総断面積は、分岐を経て下流の細い血管へ向かうにつれて増加する。特に、中動脈から細動脈までの分岐に伴う血管総断面積増加率は、参考文献２「William F.Ganong, Review of Medical Physiology 15ed.」によれば、２０倍程度と報告されている。
【０１０５】
そこで、本実施の形態では、小動脈および細動脈を末梢部モデルの対象とする。そして、一例として、分岐に伴う血管総断面積増加率は、中動脈から小動脈の過程で４倍とし、小動脈から細動脈の過程で５倍であるとし、小動脈および細動脈の長さはそれぞれ１０ｃｍおよび５ｃｍとする。また、小動脈および細動脈の血管径については、参考文献２に記載の一般的な値を採用し、それぞれの血管壁の厚さは、上流に接続されている中動脈の血管径と血管壁の厚さとの比に従って決定する。さらに、小動脈および細動脈のヤング率については、上流に接続されている中動脈のヤング率と同じ値を採用する。
【０１０６】
たとえば、橈骨動脈（図６に示すＡｖｏｌｉｏモデルの区間番号８８または９３）と、それに付加すべき末梢部モデルの形状値を以下の表に示す。
【０１０７】
【表１】
【０１０８】
（末梢部モデルの終端条件）
上述した末梢部モデルにおける細動脈の終端条件は、以下に説明するように任意に設定できる。これは、末梢血管によって構成される中動脈終端での反射率（以下、末梢反射率と称する）が細動脈の終端条件に依存しないためである。そこで、本実施の形態では、末梢部モデルにおける終端を無反射条件の下で、全身についての脈波伝播モデル（伝達関数）などを算出する。
【０１０９】
図８は、表１に示す形状値を用いて、図７に示す細動脈終端での反射率ＳＴを変化させた場合の末梢反射率の変化を示す図である。
【０１１０】
図８（ａ）は、反射率ＳＴ＝０（無反射）の場合を示し、図８（ｂ）は、反射率ＳＴ＝１（閉鎖端）の場合を示し、図８（ｃ）は、反射率ＳＴ＝−１（開放端）の場合を示す。
【０１１１】
図８（ａ）〜図８（ｃ）を参照して、細動脈終端での反射率ＳＴにかかわらず、振幅特性および周波数特性のいずれについても、略同一の挙動を示していることが分かる。すなわち、細動脈の終端条件は末梢反射率に寄与していないといえる。
【０１１２】
この現象を物理的に解釈できるように、各動脈内における脈波の伝播性状についても説明する。
【０１１３】
図９は、中動脈・小動脈・細動脈を伝播する脈波の減衰定数を数値計算によって求めた結果を示す図である。
【０１１４】
図９を参照して、脈波が伝播する血管の管径が小さくなるにつれて、減衰定数が指数的に大きくなっていることが分かる。これは、細い管内では、脈波（流体）の管壁摩擦抵抗の影響が大きくなるためであると考えられる。すなわち、細動脈終端に到達する脈波の振幅は、中動脈を伝播する脈波の振幅に比較して極めて小さくなっており、細動脈終端からの反射波は十分に減衰していると考えることができる。
【０１１５】
図１０は、図７に示す中動脈の始端（座標ｘ）から所定の周波数で加振（加圧）した場合の血管内圧力分布を示す図である。なお、図１０の座標系は、図７に示す中動脈の始端を原点（ｘ＝０）とし、下流方向をｘの正方向としたものである。
【０１１６】
図１０には、加振周波数ｆを１Ｈｚ，５Ｈｚ，１０Ｈｚ，２０Ｈｚの４種類とした場合の結果を示す。図１０から分かるように、いずれの加振周波数においても細動脈内部を伝播する脈波はその終端付近において十分に減衰している。これは、上述したように、細動脈内における脈波の距離減衰が大きいためである。
【０１１７】
したがって、細動脈の終端にどのような境界条件を設けても、細動脈の終端からの反射波が上流の循環系に影響を及ぼすことはないと考えられる。すなわち、末梢反射率は、末梢部モデルにおける終端の境界条件に依存することなく、末梢血管の形状値のみによって決定される。
【０１１８】
（位相差特性）
図１１は、一様な管路における脈波伝播の様子を示す模式図である。
【０１１９】
図１１を参照して、反射波が存在せず、脈波速度Ｃｐが（１）式で与えられる周波数に依存しない定数であるとする。すると、測定部位Ｍｐａの測定部位Ｍｐｂに対する脈波の位相遅れφは（２２）式で表される。
【０１２０】
【数１８】
【０１２１】
（２２）式を脈波速度Ｃｐおよび周波数ｆを用いて書き直せば、（２３）式が得られる。
【０１２２】
【数１９】
【０１２３】
（２２）式より、測定部位Ｍｐａ−測定部位Ｍｐｂ間の位相線図（位相差特性）は、周波数ｆの一次関数となり、その勾配は脈波速度Ｃｐに応じた値になることが分かる。さらに、（１）式および（２３）式を用いて（２４）式が得られる。
【０１２４】
【数２０】
【０１２５】
（２４）式より、血管壁のヤング率Ｅが大きくなるほど位相線図の勾配が緩やかになることが分かる。
【０１２６】
（フィッティング）
再度、図２を参照して、上述したように２つの測定部位間の位相差特性は周波数ｆの一次関数となるので、探索部４０は、この一次関数の傾きが実測値と一致するように、モデルをフィッティングする。
【０１２７】
より具体的には、伝達関数算出部３６は、Ａｖｏｌｉｏモデル中の各区間のヤング率Ｅｎ（ｎ＝１〜１２８）に変数ｋを乗じて得られる、仮のヤング率ｋ・Ｅｎを用いて伝達関数Ｇａ（ｆ）およびＧｂ（ｆ）を算出する。そして、探索部４０は、位相線傾き算出部（実測）３２が算出する位相線の傾きｇｅｘｐと、位相線傾き算出部（モデル）３８が算出する位相線の傾きｇ（ｋ）との偏差Δ（＝｜ｇｅｘｐ−ｇ（ｋ）｜）が最小化するように、変数ｋを最適化する。この最適化処理は、代表的に数理計画法（たとえば、最小二乗法）を用いて行なわれるが、数理計画法については公知であるので、詳細な説明は行なわない。
【０１２８】
（実測の位相差特性の算出）
以下では、本実施の形態に従う血管状態評価装置１００を用いて、２人の被験者２００ａおよび２００ｂについて、実際に測定を行なった結果を示す。
【０１２９】
図１２は、被験者２００ａ，２００ｂの上腕および足関節に押圧カフを装着して測定した圧力の時間波形を示す図である。なお、測定信号Ｐａ（ｔ）が足関節における圧力を示し、測定信号Ｐｂ（ｔ）が上腕における圧力を示す。
【０１３０】
図１３は、図１２に示す上腕の圧力波形と足関節の圧力波形との間のコヒーレンスを示す図である。
【０１３１】
図１２（ａ）および図１３（ａ）は、被験者２００ａの測定結果を示し、図１２（ｂ）および図１３（ｂ）は、被験者２００ｂの測定結果を示す。
【０１３２】
ここで、コヒーレンスとは、波形間の周波数領域での相関性を示す指標であり、コヒーレンスが大きいほど（１に近いほど）両波形の相関性が高いことになる。なお、図１３に示すコヒーレンスは、（２５）式に従って算出した。
【０１３３】
【数２１】
【０１３４】
図１３（ａ）および図１３（ｂ）に示すように、特定の周波数においてコヒーレンスが低くなっていることが分かる。このようなコヒーレンス低下の力学的な要因としては、様々なものが考えられるが、代表的に、人体が有する非線形性や、測定部位での圧力が特定の周波数において節になってしまうことなどが影響していると考えられる。また、被験者の姿勢や測定中のわずかな動きといった人為的な要因によってもコヒーレンスが低下することも考えられる。
【０１３５】
このようなコヒーレンスの低いデータは、解析上の誤差を拡大してしまう恐れがあるため、測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間のコヒーレンスが予め定められたしきい値（たとえば、０．７）より小さいデータを排除することが望ましい。そのため、位相線傾き算出部（実測）３２（図２）は、測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間のコヒーレンスが予め定められたしきい値より高い周波数成分のみを用いて、位相線の傾きｇｅｘｐを算出する。
【０１３６】
図１４は、図１２に示す測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間の各周波数成分についての位相差をプロットした位相線図である。図１４（ａ）は、図１２（ａ）に対応する位相線図を示し、図１４（ｂ）は、図１２（ｂ）に対応する位相線図を示す。なお、図１４（ａ）および図１４（ｂ）では、コヒーレンスが「０．７」より低いデータを排除している。
【０１３７】
図１４（ａ）および図１４（ｂ）を参照して、それぞれの位相線図は、±１８０°を境界として不連続点を有している。これは、所定の周波数以上の周波数成分において、１周期（３６０°）以上の位相差が生じていることを意味する。そこで、位相線傾き算出部（実測）３２（図２）は、このような位相線図の不連続点に対して、１または２以上の周期に相当する単位（ｎ×３６０°）で補正を行なった上で、実測の位相差特性を算出する。
【０１３８】
図１５および図１６は、位相線傾き算出部（実測）３２が実施する位相線図の補正処理を説明するための模式図である。
【０１３９】
図１５（ａ）を参照して、測定信号Ｐａ（ｔ）を周波数変換して得られる位相特性Ｐａ（ｆ）と、測定信号Ｐｂ（ｔ）を周波数変換して得られる位相特性Ｐｂ（ｆ）とを比較し、周波数ｆｉに対応する位相差Ａｉを位相線図上にプロットする。なお、周波数ｆｉは、低周波側から数えてｉ番目の周波数成分である。そして、位相線図上にプロットされる位相差Ａｉのうち、不連続点が存在しない範囲のｎ個の位相差｛Ａ１，Ａ２，・・・，Ａｎ｝を用いて、初期回帰直線ｌ０を算出する（図１５（ｂ））。
【０１４０】
次に、ｎ＋１番目の位相差Ａｎ＋１と、初期回帰直線ｌ０の周波数ｆｎ＋１に対応する位相とが比較される。図１５（ｃ）に示すように、両者の偏差Δφ（０，ｉ＋１）が１８０°未満であれば、初期基準回帰ｌ０の算出に用いたｎ個の位相差｛Ａ１，Ａ２，・・・，Ａｎ｝にｎ＋１番目の位相差Ａｎ＋１を加えた位相差群｛Ａ１，Ａ２，・・・，Ａｎ，Ａｎ＋１｝を用いて、回帰直線ｌ１を算出する（図１５（ｄ））。
【０１４１】
一方、図１６（ｃ）に示すように、両者の偏差Δφ（０，ｉ＋１）が１８０°以上であれば、不連続点があると判断される。そして、初期回帰直線ｌ０に対する偏差が１８０°未満となるように、位相差Ａｎ＋１に３６０°×ｍ（ｍは、１以上の整数）を減算して、位相差Ａｎ＋１を位相差＃Ａｎ＋１に遷移させる。すなわち、測定されたデータが有する見かけ上の誤差を補正する。
【０１４２】
そして、初期基準回帰ｌ０の算出に用いたｎ個の位相差｛Ａ１，Ａ２，・・・，Ａｎ｝に、この補正された位相差＃Ａｎ＋１を加えた位相差群｛Ａ１，Ａ２，・・・，Ａｎ，＃Ａｎ＋１｝を用いて、回帰直線ｌ１を算出する（図１６（ｄ））。
【０１４３】
以下、同様にして、すべての位相差Ａｉについて、プロットおよび回帰直線の更新が繰返される。
【０１４４】
図１７は、図１４に示す位相線図を補正して連続化した結果を示す図である。図１７（ａ）は、図１４（ａ）に対応する位相線図を示し、図１７（ｂ）は、図１４（ｂ）に対応する位相線図を示す。
【０１４５】
図１７（ａ）および図１７（ｂ）を参照して、上述の方法によって、各位相差が補正されて位相線図が連続化されていることが分かる。また、図１７（ａ）および図１７（ｂ）には、プロットした位相データの回帰直線についてもそれぞれ示されており、これらの課回帰直線の傾きが図２に示す傾きｇｅｘｐに相当する。
【０１４６】
（フローチャート）
図１８は、この発明の実施の形態１に従う血管状態評価装置１００において実行される処理の手順を示すフローチャートである。図１８のフローチャートに示される各処理は、制御部２のＣＰＵ１０がＲＯＭ１２に予め格納されているプログラムを読出し、ＲＡＭ１３上に展開して実行することで、図２に示す各機能を実現する。
【０１４７】
図１８を参照して、ユーザによる操作部６などの操作に応答して、ＣＰＵ１０は、測定部２０ａ，２０ｂに対して測定指令を与え、測定部２０ａ，２０ｂが被験者２００の所定の測定部位における生体信号の測定を開始する（ステップＳ１００）。
【０１４８】
次に、ＣＰＵ１０は、測定部２０ａ，２０ｂで測定される時間波形である測定信号Ｐａ（ｔ），Ｐｂ（ｔ）を周波数領域の位相特性Ｐａ（ｆ），Ｐｂ（ｆ）に変換する（ステップＳ１０２）。そして、ＣＰＵ１０は、位相特性Ｐａ（ｆ）と位相特性Ｐｂ（ｆ）との間の各周波数成分についての位相差に基づいて、実測の位相差特性（傾きｇｅｘｐ）を算出する（ステップＳ１０４）。
【０１４９】
また、ＣＰＵ１０は、変数ｋを初期値ｋ０に設定する（ステップＳ１０６）。そして、ＲＯＭ１２などに格納されている循環系モデルを参照して、心臓からそれぞれ押圧カフ２４ａ，２４ｂが装着された２つの測定部位までの血管経路の伝達特性を示す２つの伝達関数を算出する（ステップＳ１０８）。この伝達関数の算出に際して、各区間のヤング率は、循環系モデル中の規定された基準のヤング率に変数ｋを乗じた値が用いられる。さらに、ＣＰＵ１０は、ステップＳ１０６で算出した２つの伝達関数の間の位相差特性（傾きｇ（ｋ））を算出する（ステップＳ１１０）。
【０１５０】
その後、ＣＰＵ１０は、ステップＳ１０４において算出された実測の位相差特性（傾きｇｅｘｐ）と、ステップＳ１１０において算出された伝達関数の間の位相差特性（傾きｇ（ｋ））との偏差Δ（＝｜ｇｅｘｐ−ｇ（ｋ）｜）を算出する（ステップＳ１１２）。そして、ＣＰＵ１０は、偏差Δが予め定められた収束条件を満足しているか否かを判断する（ステップＳ１１４）。この収束条件としては、代表的に、偏差Δが予め定められたしきい値より小さいか否かが判断される。
【０１５１】
偏差Δが予め定められた収束条件を満足していなければ（ステップＳ１１４においてＮＯ）、ＣＰＵ１０は、偏差Δを小さくする方向に、変数ｋを所定値だけ増減させる（ステップＳ１１６）。そして、ステップＳ１０８以降の処理を再度実行する。
【０１５２】
一方、偏差Δが予め定められた収束条件を満足していれば（ステップＳ１１４においてＹＥＳ）、ＣＰＵ１０は、当該時点の変数ｋの値を最適解ｋｏｐｔとして決定する（ステップＳ１１８）。そして、ＣＰＵ１０は、決定した最適解ｋｏｐｔ、最適解ｋｏｐｔを用いて換算されたヤング率、最適解ｋｏｐｔの評価結果などを、表示部４へ出力する（ステップＳ１２０）。そして、評価処理は終了する。
【０１５３】
なお、上述の説明においては、心臓を入力端とした全身についての脈波伝播モデル（伝達関数）、末梢部モデル、および伝達関数Ｇａ（ｆ），Ｇｂ（ｆ）を算出する方法について詳述したが、このようなモデルや伝達関数を必ずしも評価処理毎に算出する必要はない。すなわち、評価処理前に算出したモデルや伝達関数を記憶部３４に予め格納しておいてもよい。
【０１５４】
この発明の実施の形態１によれば、実測の測定信号の周波数特性を反映した弾力度変数ｋの最適解を得ることができる。これにより、周波数特性を考慮した上で、血管の状態（動脈硬化の度合い）をより高い精度で評価することができる。
【０１５５】
［実施の形態２］
上述のこの発明の実施の形態１では、血管の弾力度合いを示す弾力度変数を算出する構成について例示したが、従来から実用化されている脈波速度法と同様の評価基準を用いることができるように、脈波速度を算出することも有効である。そのため、この発明の実施の形態２では、予め定められた２点間の脈波速度を算出する構成について説明する。
【０１５６】
この発明の実施の形態２に従う血管状態評価装置１００＃の構成は、制御部２＃で実行される処理を除いて、図１に示すこの発明の実施の形態１に従う血管状態評価装置１００と同様であるので、詳細な説明は繰返さない
図１９は、この発明の実施の形態２に従う血管状態評価装置１００＃の制御部２＃で実行される機能を模式的に示す機能ブロック図である。
【０１５７】
図１９を参照して、制御部２＃は、図２に示すこの発明の実施の形態１に従う制御部２で実行される機能において、評価部４２に代えて、脈波速度モデル算出部４４と、平均脈波速度算出部４６と、評価部４８とを設けたものである。その他の機能については、図２と同様であるので、詳細な説明は繰返さない。
【０１５８】
脈波速度モデル算出部４４は、探索部４０によってフィッティングされた最適解ｋｏｐｔに基づいて、予め定められた２点間を伝播する脈波速度を算出するための数学的なモデルを算出する。ここで、２点間を伝播する脈波速度とは、２つの測定部位間において空間的に平均化した脈波速度を意味する。すなわち、上述した循環系モデルの各区間の脈波速度は、管径および管長に応じて変化するため、２つの測定部位間の伝播経路の形状に応じて、脈波速度は増速したり減速したりする。そこで、従来の脈波速度法との整合性をとるために、以下に示すような空間的な平均化処理を行なって、測定部位間の脈波速度（以下、「平均脈波速度」とも称す）を算出する。脈波速度モデル算出部４４は、このような平均化処理を含んだモデルを算出する。
【０１５９】
平均脈波速度算出部４６は、脈波速度モデル算出部４４で算出された数学モデルに基づいて演算を行ない、平均脈波速度Ｃａｖｅを算出する。
【０１６０】
評価部４８は、平均脈波速度算出部４６で算出された平均脈波速度Ｃａｖｅを予め定められた基準値と比較し、動脈硬化の度合いについての評価を、表示部４（図１）などに出力する。
【０１６１】
（平均脈波速度の算出）
記憶部３４に予め格納されている循環系モデル、および探索部４０によってフィッティングされた最適解ｋｏｐｔに基づいて、２つの測定部位間の伝播経路における血管の形状値は既知となっている。
【０１６２】
図２０は、２つの測定部位ＭｐＡ，ＭｐＢの間の経路を模式的に示す管路モデルの図である。
【０１６３】
図２０を参照して、２つの測定部位ＭｐＡ，ＭｐＢの間の経路に、ｎ個の要素管路（区間）が直列的に接続されているとする。区間ｉについての区間長をＬｉとし、脈波速度をＣｉとし、脈波の伝播に要する時間をｔｉとすれば、測定部位ＭｐＡ−ＭｐＢ間の平均脈波速度は、（２６）式で与えられる。
【０１６４】
【数２２】
【０１６５】
ここで、ｔｉ＝Ｌｉ／Ｃｉであることを用いれば、（２６）式は（２７）式にように表される。
【０１６６】
【数２３】
【０１６７】
なお、各区間の脈波速度Ｃｉは、（１２）式を用いて、（２８）式のように表される。
【０１６８】
【数２４】
【０１６９】
このように、各区間の脈波速度を厳密に評価することで、測定部位間の平均脈波速度について周波数特性を反映した上で正確に算出することができる。さらに、各区間の脈波速度が経路によって異なる場合においても、測定部位間の平均的な脈波速度として有意な値を保証することができる。
【０１７０】
（評価例）
図２１は、図１２に示す被験者２００ａ，２００ｂから実測した測定信号Ｐａ（ｔ），Ｐｂ（ｔ）に基づいて、平均脈波速度を算出した結果を示す図である。図２１（ａ）は、被験者２００ａの算出結果を示し、図２１（ｂ）は、被験者２００ｂの算出結果を示す。
【０１７１】
図２１（ａ）および図２１（ｂ）を参照して、被験者２００ａおよび２００ｂのそれぞれについて、平均脈波速度が周波数特性を反映した上で同定されていることが分かる。なお、被験者２００ａおよび２００ｂの算出結果を比較すると、被験者２００ｂの脈波速度がより大きく、相対的に動脈硬化が進行していることを示している。
【０１７２】
図２２は、この発明の実施の形態２に従う血管状態評価方法による算出結果を従来の脈波速度法（ｂａＰＷＶ法）によって得られた測定結果と比較した図である。なお、図２２は、２３名の被験者を対象とした結果を示す。また、上述したように、各区間の脈波速度Ｃｉは周波数依存性を有しているため、各測定体の拍動周期に対応した周波数に対応する脈波速度Ｃｉおよび平均脈波速度Ｃａｖｅを採用した。
【０１７３】
図２２を参照して、この発明の実施の形態２に従う血管状態評価方法による算出結果と、従来の脈波速度法によって得られた測定結果との間の相関係数は０．９３となった。この結果から、本実施の形態に従う血管状態評価方法は、従来の脈波速度法による測定結果と比較的高い相関を有しており、脈波速度法によって蓄積された評価基準と同様の評価基準を用いて、動脈硬化の度合いを評価することができる。
【０１７４】
（フローチャート）
図２３は、本実施の形態２に従う血管状態評価装置１００＃において実行される処理の手順を示すフローチャートである。図２３のフローチャートに示される各処理は、制御部２＃のＣＰＵ１０がＲＯＭ１２に予め格納されているプログラムを読出し、ＲＡＭ１３上に展開して実行することで、図１９に示す各機能を実現する。
【０１７５】
図２３を参照して、まず、ＣＰＵ１０は、図１８に示すステップＳ１００〜ステップＳ１１８と同様の処理を実行する。なお、これらの処理は、図１８と同様であるので、詳細な説明は繰返さない。
【０１７６】
そして、ＣＰＵ１０は、ステップＳ１１８において決定した最適解ｋｏｐｔと、ＲＯＭ１２などに格納されている循環系モデルの形状値とに基づいて、予め定められた２つの測定部位間の脈波速度モデルを算出する（ステップＳ２０２）。さらに、ＣＰＵ１０は、算出した脈波速度モデルに従って、平均脈波速度を算出する（ステップＳ２０４）。
【０１７７】
さらに、ＣＰＵ１０は、ステップＳ２０４で算出した平均脈波速度の評価結果などを、表示部４へ出力する（ステップＳ２０６）。そして、評価処理は終了する。
【０１７８】
この発明の実施の形態２によれば、実測の測定信号の周波数特性を反映した平均脈波速度を高い精度で算出できる。この平均脈波速度は、従来の脈波速度法を用いて測定した結果とも比較的高い相関値をもつので、従来から蓄積されてきた脈波速度法における評価基準を用いて、血管の状態を判定することもできる。
【０１７９】
［その他の形態］
さらに、本実施の形態に従う血管状態評価装置における評価方法を実現させるためのプログラムを提供することもできる。このようなプログラムは、コンピュータに付属するフレキシブルディスク、ＣＤ−ＲＯＭ（Compact Disk-Read Only Memory）、ＲＯＭ、ＲＡＭおよびメモリカードなどのコンピュータ読取り可能な記録媒体に記録させた上、プログラム製品として提供することもできる。あるいは、コンピュータに内蔵するハードディスクなどの記録媒体にて記録させて、プログラムとして提供することもできる。また、ネットワークを介したダウンロードによって、プログラムとして提供することもできる。
【０１８０】
なお、本発明にかかるプログラムは、コンピュータのオペレーションシステム（ＯＳ）の一部として提供されるプログラムモジュールのうち、必要なモジュールを所定の配列で所定のタイミングで呼出して処理を実行させるものであってもよい。その場合、プログラム自体には上記モジュールが含まれずＯＳと協働して処理が実行される。このようなモジュールを含まないプログラムも、本発明にかかるプログラムに含まれ得る。
【０１８１】
また、本発明にかかるプログラムは、たとえば通常の血圧測定を実行するためのプログラムといった他のプログラムの一部に組込まれて提供されるものであってもよい。その場合にも、プログラム自体には当該他のプログラムに含まれるモジュールが含まれず、他のプログラムと協働して処理が実行される。このような他のプログラムに組込まれたプログラムも、本発明にかかるプログラムに含まれ得る。
【０１８２】
提供されるプログラム製品は、ハードディスクなどのプログラム格納部にインストールされて実行される。なお、プログラム製品は、プログラム自体と、プログラムが記録された記録媒体とを含む。
【０１８３】
今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は、上記した説明ではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。
【０１８４】
（付表）
Ａｖｏｌｉｏモデルの主要な形状値を以下の付表に示す。
【０１８５】
【表２】
【０１８６】
【表３】
【０１８７】
【表４】
【０１８８】
【表５】
【０１８９】
【表６】
【０１９０】
【表７】
【図面の簡単な説明】
【０１９１】
【図１】この発明の実施の形態１に従う血管状態評価装置の概略構成図である。
【図２】この発明の実施の形態１に従う血管状態評価装置の制御部で実行される機能を模式的に示す機能ブロック図である。
【図３】血管内の血液の一次元流れモデルを示す図である。
【図４】図３に示す検査体積に働く力および出入りする運動量を示す図である。
【図５】血管を一次元線形分布定数モデル化した模式図である。
【図６】Ａｖｏｌｉｏモデルの模式図である。
【図７】末梢部モデルの模式図である。
【図８】表１に示す形状値を用いて、図７に示す細動脈終端での反射率ＳＴを変化させた場合の末梢反射率の変化を示す図である。
【図９】中動脈・小動脈・細動脈を伝播する脈波の減衰定数を数値計算によって求めた結果を示す図である。
【図１０】図７に示す中動脈の始端（座標ｘ）から所定の周波数で加振（加圧）した場合の血管内圧力分布を示す図である。
【図１１】一様な管路における脈波伝播の様子を示す模式図である。
【図１２】被験者の上腕および足関節に押圧カフを装着して測定した圧力の時間波形を示す図である。
【図１３】図１２に示す上腕の圧力波形と足関節の圧力波形との間のコヒーレンスを示す図である。
【図１４】図１２に示す測定信号Ｐａ（ｔ）と測定信号Ｐｂ（ｔ）との間の各周波数成分についての位相差をプロットした位相線図である。
【図１５】位相線傾き算出部（実測）が実施する位相線図の補正処理を説明するための模式図である。
【図１６】位相線傾き算出部（実測）が実施する位相線図の補正処理を説明するための模式図である。
【図１７】図１４に示す位相線図を補正して連続化した結果を示す図である。
【図１８】この発明の実施の形態１に従う血管状態評価装置において実行される処理の手順を示すフローチャートである。
【図１９】この発明の実施の形態２に従う血管状態評価装置の制御部で実行される機能を模式的に示す機能ブロック図である。
【図２０】２つの測定部位ＭｐＡ，ＭｐＢの間の経路を模式的に示す管路モデルの図である。
【図２１】図１２に示す被験者から実測した測定信号Ｐａ（ｔ），Ｐｂ（ｔ）に基づいて、平均脈波速度を算出した結果を示す図である。
【図２２】この発明の実施の形態２に従う血管状態評価方法による算出結果を従来の脈波速度法（ｂａＰＷＶ法）によって得られた測定結果と比較した図である。
【図２３】本実施の形態２に従う血管状態評価装置において実行される処理の手順を示すフローチャートである。
【符号の説明】
【０１９２】
２，２＃ 制御部、４ 表示部、６ 操作部、１０ ＣＰＵ、１２ ＲＯＭ、１４ ＲＡＭ、２０ａ，２０ｂ 測定部、２２ａ，２２ｂ，２７ａ，２７ｂ 配管、、２４ａ，２４ｂ 押圧カフ、２５ａ，２５ｂ 圧力ポンプ、２６ａ，２６ｂ 調圧弁、２８ａ，２８ｂ 圧力センサ、３０ａ，３０ｂ 周波数変換部（ＦＦＴ）、３２ 位相線傾き算出部（実測）、３４ 記憶部、３６ 伝達関数算出部、３８ 位相線傾き算出部（モデル）、４０ 探索部、４２，４８ 評価部、４４ 脈波速度モデル算出部、４６ 平均脈波速度算出部、５０ 検査体積、１００，１００＃ 血管状態評価装置（評価装置）。
【特許請求の範囲】
【請求項１】
生体を構成する血管を複数の区間に分割してモデル化した循環系モデルを格納する記憶手段を備え、
前記循環系モデルは、前記複数の区間の各々を代表する形状値を含み、さらに
生体の第１測定部位に装着されて第１生体信号の時間波形を測定する第１測定手段と、
生体の第２測定部位に装着され、前記第１測定手段と同期して第２生体信号の時間波形を測定する第２測定手段と、
前記第１生体信号と前記第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出する第１算出手段と、
前記第１測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第１伝達関数と、前記第２測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出する第２算出手段とを備え、
前記第１伝達関数および前記第２伝達関数は、血管の弾力度合いを示す弾力度変数を含み、さらに、
前記第１算出手段によって算出される実測の位相差特性に基づいて、前記第２算出手段によって算出される位相差特性をフィッティングすることにより、前記弾力度変数を決定する探索手段を備える、血管状態評価装置。
【請求項２】
前記第１および第２測定部位までの血管経路にそれぞれに対応する各区間の前記形状値に基づいて、前記第１および第２伝達関数を算出する伝達関数算出手段をさらに備える、請求項１に記載の血管状態評価装置。
【請求項３】
前記伝達関数算出手段は、各区間に対応付けられた、血管の圧力および血液流量を入力変数とする分布定数モデルを用いて、前記第１および第２伝達関数を算出し、
前記分布定数モデルの各々は、対応する区間における血液の流れやすさに応じた縦インピーダンスと、前記弾力度変数を有する横インピーダンスとを含む、請求項２に記載の血管状態評価装置。
【請求項４】
前記探索手段によってフィッティングされた前記弾力度変数に基づいて、血管における脈波速度を算出する脈波速度算出手段をさらに備える、請求項１〜３のいずれか１項に記載の血管状態評価装置。
【請求項５】
前記脈波速度算出手段は、前記第１測定部位までの血管経路に対応する各区間の前記形状値および前記第２測定部位までの血管経路に対応する各区間の前記形状値に基づいて前記脈波速度を算出する、請求項４に記載の血管状態評価装置。
【請求項６】
前記循環系モデルは、前記形状値として、血管径および血管長を含む、請求項１〜５のいずれか１項に記載の血管状態評価装置。
【請求項７】
前記循環系モデルは、生体を構成する血管を複数の区分に分類した上で、前記複数の区分の少なくとも１つの区分に属する血管をモデル化したものである、請求項１〜６のいずれか１項に記載の血管状態評価装置。
【請求項８】
前記生体を構成する血管は、血管径の大きさに基づいて複数の区分に分類される、請求項７に記載の血管状態評価装置。
【請求項９】
前記伝達モデル算出手段は、各区間に含まれる血管のうち、前記循環系モデルにおいてモデル化されていない血管をモデル化した末梢部モデルを、各区間に対応する前記循環系モデルに付加した上で、前記伝達関数を算出する、請求項２または３に記載の血管状態評価装置。
【請求項１０】
前記伝達モデル算出手段は、各区間の前記循環系モデルを血管の形状差に基づいて変換することで、当該区間の前記末梢部モデルを算出する、請求項９に記載の血管状態評価装置。
【請求項１１】
前記伝達モデル算出手段は、前記末梢部モデルの終端を無反射条件の下で前記伝達関数を算出する、請求項９または１０に記載の血管状態評価装置。
【請求項１２】
前記第１生体信号から各周波数成分についての位相を示す第１位相特性を算出する第１周波数変換手段と、
前記第２生体信号から各周波数成分についての位相を示す第２位相特性を算出する第２周波数変換手段とをさらに備え、
前記第１算出手段は、
前記第１位相データと前記第２位相データとを差分して差分位相データを算出し、さらに
前記差分位相データにおける周期遅れに起因する位相ずれを、１または２以上の周期に相当する位相の単位で補正することで、前記実測の位相差特性を算出する、請求項１〜１１のいずれか１項に記載の血管状態評価装置。
【請求項１３】
前記第１算出手段は、前記第１生体信号と前記第２生体信号との間のコヒーレンス値が予め定められたしきい値より高い周波数成分を用いて前記実測の位相差特性を算出する、請求項１〜１２のいずれか１項に記載の血管状態評価装置。
【請求項１４】
生体を構成する血管を複数の区間に分割してモデル化した循環系モデルを用いて、生体を構成する血管の状態を評価する血管状態評価方法であって、
前記循環系モデルは、前記複数の区間の各々を代表する形状値を含み、
前記血管状態評価方法は、
生体の第１測定部位から第１生体信号の時間波形を測定するとともに、生体の第２測定部位から第２生体信号の時間波形を測定するステップと、
前記第１生体信号と前記第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するステップと、
前記第１測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第１伝達関数と、前記第２測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出するステップとを備え、
前記第１伝達関数および前記第２伝達関数は、血管の弾力度合いを示す弾力度変数を含み、さらに、
前記実測の位相差特性に基づいて、前記第１伝達関数と前記第２伝達関数との間の位相差特性をフィッティングすることにより、前記弾力度変数を決定するステップを備える、血管状態評価方法。
【請求項１５】
演算処理部を有するコンピュータに、生体を構成する血管を複数の区間に分割して予めモデル化した循環系モデルを用いて、生体を構成する血管の状態の評価を実行させるための血管状態評価プログラムであって、
前記循環系モデルは、前記複数の区間の各々を代表する形状値を含み、
前記血管状態評価プログラムは、
前記演算処理部が、生体の第１測定部位で測定される第１生体信号の時間波形を取得するとともに、生体の第２測定部位で測定される第２生体信号の時間波形を取得するステップと、
前記演算処理部が、前記第１生体信号と前記第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するステップと、
前記演算処理部が、前記第１測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第１伝達関数と、前記第２測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出するステップとを備え、
前記第１伝達関数および前記第２伝達関数は、血管の弾力度合いを示す弾力度変数を含み、さらに、
前記演算処理部が、前記実測の位相差特性に基づいて、前記第１伝達関数と前記第２伝達関数との間の位相差特性をフィッティングすることにより、前記弾力度変数を決定するステップを備える、血管状態評価プログラム。
【請求項１】
生体を構成する血管を複数の区間に分割してモデル化した循環系モデルを格納する記憶手段を備え、
前記循環系モデルは、前記複数の区間の各々を代表する形状値を含み、さらに
生体の第１測定部位に装着されて第１生体信号の時間波形を測定する第１測定手段と、
生体の第２測定部位に装着され、前記第１測定手段と同期して第２生体信号の時間波形を測定する第２測定手段と、
前記第１生体信号と前記第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出する第１算出手段と、
前記第１測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第１伝達関数と、前記第２測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出する第２算出手段とを備え、
前記第１伝達関数および前記第２伝達関数は、血管の弾力度合いを示す弾力度変数を含み、さらに、
前記第１算出手段によって算出される実測の位相差特性に基づいて、前記第２算出手段によって算出される位相差特性をフィッティングすることにより、前記弾力度変数を決定する探索手段を備える、血管状態評価装置。
【請求項２】
前記第１および第２測定部位までの血管経路にそれぞれに対応する各区間の前記形状値に基づいて、前記第１および第２伝達関数を算出する伝達関数算出手段をさらに備える、請求項１に記載の血管状態評価装置。
【請求項３】
前記伝達関数算出手段は、各区間に対応付けられた、血管の圧力および血液流量を入力変数とする分布定数モデルを用いて、前記第１および第２伝達関数を算出し、
前記分布定数モデルの各々は、対応する区間における血液の流れやすさに応じた縦インピーダンスと、前記弾力度変数を有する横インピーダンスとを含む、請求項２に記載の血管状態評価装置。
【請求項４】
前記探索手段によってフィッティングされた前記弾力度変数に基づいて、血管における脈波速度を算出する脈波速度算出手段をさらに備える、請求項１〜３のいずれか１項に記載の血管状態評価装置。
【請求項５】
前記脈波速度算出手段は、前記第１測定部位までの血管経路に対応する各区間の前記形状値および前記第２測定部位までの血管経路に対応する各区間の前記形状値に基づいて前記脈波速度を算出する、請求項４に記載の血管状態評価装置。
【請求項６】
前記循環系モデルは、前記形状値として、血管径および血管長を含む、請求項１〜５のいずれか１項に記載の血管状態評価装置。
【請求項７】
前記循環系モデルは、生体を構成する血管を複数の区分に分類した上で、前記複数の区分の少なくとも１つの区分に属する血管をモデル化したものである、請求項１〜６のいずれか１項に記載の血管状態評価装置。
【請求項８】
前記生体を構成する血管は、血管径の大きさに基づいて複数の区分に分類される、請求項７に記載の血管状態評価装置。
【請求項９】
前記伝達モデル算出手段は、各区間に含まれる血管のうち、前記循環系モデルにおいてモデル化されていない血管をモデル化した末梢部モデルを、各区間に対応する前記循環系モデルに付加した上で、前記伝達関数を算出する、請求項２または３に記載の血管状態評価装置。
【請求項１０】
前記伝達モデル算出手段は、各区間の前記循環系モデルを血管の形状差に基づいて変換することで、当該区間の前記末梢部モデルを算出する、請求項９に記載の血管状態評価装置。
【請求項１１】
前記伝達モデル算出手段は、前記末梢部モデルの終端を無反射条件の下で前記伝達関数を算出する、請求項９または１０に記載の血管状態評価装置。
【請求項１２】
前記第１生体信号から各周波数成分についての位相を示す第１位相特性を算出する第１周波数変換手段と、
前記第２生体信号から各周波数成分についての位相を示す第２位相特性を算出する第２周波数変換手段とをさらに備え、
前記第１算出手段は、
前記第１位相データと前記第２位相データとを差分して差分位相データを算出し、さらに
前記差分位相データにおける周期遅れに起因する位相ずれを、１または２以上の周期に相当する位相の単位で補正することで、前記実測の位相差特性を算出する、請求項１〜１１のいずれか１項に記載の血管状態評価装置。
【請求項１３】
前記第１算出手段は、前記第１生体信号と前記第２生体信号との間のコヒーレンス値が予め定められたしきい値より高い周波数成分を用いて前記実測の位相差特性を算出する、請求項１〜１２のいずれか１項に記載の血管状態評価装置。
【請求項１４】
生体を構成する血管を複数の区間に分割してモデル化した循環系モデルを用いて、生体を構成する血管の状態を評価する血管状態評価方法であって、
前記循環系モデルは、前記複数の区間の各々を代表する形状値を含み、
前記血管状態評価方法は、
生体の第１測定部位から第１生体信号の時間波形を測定するとともに、生体の第２測定部位から第２生体信号の時間波形を測定するステップと、
前記第１生体信号と前記第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するステップと、
前記第１測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第１伝達関数と、前記第２測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出するステップとを備え、
前記第１伝達関数および前記第２伝達関数は、血管の弾力度合いを示す弾力度変数を含み、さらに、
前記実測の位相差特性に基づいて、前記第１伝達関数と前記第２伝達関数との間の位相差特性をフィッティングすることにより、前記弾力度変数を決定するステップを備える、血管状態評価方法。
【請求項１５】
演算処理部を有するコンピュータに、生体を構成する血管を複数の区間に分割して予めモデル化した循環系モデルを用いて、生体を構成する血管の状態の評価を実行させるための血管状態評価プログラムであって、
前記循環系モデルは、前記複数の区間の各々を代表する形状値を含み、
前記血管状態評価プログラムは、
前記演算処理部が、生体の第１測定部位で測定される第１生体信号の時間波形を取得するとともに、生体の第２測定部位で測定される第２生体信号の時間波形を取得するステップと、
前記演算処理部が、前記第１生体信号と前記第２生体信号との間の各周波数成分についての位相差に基づいて、実測の位相差特性を算出するステップと、
前記演算処理部が、前記第１測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第１伝達関数と、前記第２測定部位までの血管経路に対応して前記循環系モデルに基づいて規定された第２伝達関数との間の位相差特性を算出するステップとを備え、
前記第１伝達関数および前記第２伝達関数は、血管の弾力度合いを示す弾力度変数を含み、さらに、
前記演算処理部が、前記実測の位相差特性に基づいて、前記第１伝達関数と前記第２伝達関数との間の位相差特性をフィッティングすることにより、前記弾力度変数を決定するステップを備える、血管状態評価プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【公開番号】特開２００８−２４６０１０（Ｐ２００８−２４６０１０Ａ）
【公開日】平成２０年１０月１６日（２００８．１０．１６）
【国際特許分類】
【出願番号】特願２００７−９２５２２（Ｐ２００７−９２５２２）
【出願日】平成１９年３月３０日（２００７．３．３０）
【出願人】（５０４１３２２７２）国立大学法人京都大学 (1,269)
【出願人】（５０３２４６０１５）オムロンヘルスケア株式会社 (584)
【Ｆターム（参考）】
【公開日】平成２０年１０月１６日（２００８．１０．１６）
【国際特許分類】
【出願日】平成１９年３月３０日（２００７．３．３０）
【出願人】（５０４１３２２７２）国立大学法人京都大学 (1,269)
【出願人】（５０３２４６０１５）オムロンヘルスケア株式会社 (584)
【Ｆターム（参考）】
[ Back to top ]