JP2014100249A - 血管解析装置、医用画像診断装置、血管解析方法、及び血管解析プログラム - Google Patents

Abstract

【課題】血管の構造流体解析の精度の向上ができる血管解析装置を提供する。
【解決手段】記憶装置は、被検体の血管に関する時系列の医用画像のデータを記憶する。領域設定部５１は、時系列の医用画像に含まれる血管領域に解析対象領域を設定し、解析対象領域に潜在変数同定領域を設定する。画像解析・追尾処理部５３は、時系列の医用画像を画像処理して解析対象領域の時系列の形態指標と時系列の形状変形指標とを算出する。力学モデル構築部５５は、時系列の形態指標と時系列の形状変形指標と時系列の医用画像とに基づいて解析対象領域の構造流体解析に関する力学モデルを暫定的に構築する。統計的同定部は、暫定的に構築された力学モデルに基づく血管形態指標の予測値及び血液流量指標の予測値の少なくとも一方が予め計測された血管形態指標の観測値及び血液流量指標の観測値の少なくとも一方に整合するように潜在変数同定領域に関する潜在変数を同定する。
【選択図】図４

Description

本発明の実施形態は、血管解析装置、医用画像診断装置、血管解析方法、及び血管解析プログラムに関する。

三大疾患の一つである心疾患の原因となる冠動脈の狭窄や、脳動脈瘤、あるいはそれらの予兆となる頸動脈のプラークによる狭窄を、非侵襲あるいは低侵襲に予防・診断するための技術が望まれている。

冠動脈の狭窄は、虚血性心疾患に至る重大な病変である。冠動脈の狭窄診断としては、カテーテルによる冠動脈造影検査（ＣＡＧ：Coronary Angiography）が主流である。冠動脈の器質的病変の診断指標として、心筋血流予備量比（ＦＦＲ：Fractional Flow Reserve）がある。ＦＦＲは、狭窄非存在下の最大冠血流に対する狭窄存在下の最大冠血流の比率として定義される。ＦＦＲは、狭窄近位部冠内圧に対する狭窄遠位部冠内圧の比率に略一致する。カテーテル先端に設けられた圧力センサが測定される。すなわち、ＦＦＲの測定は、カテーテル手術を必要とする。

冠動脈の狭窄解析が心臓ＣＴで可能であれば、カテーテル手術によるＦＦＲの測定に比して、低侵襲、患者の負担低減、及び医療コストが節約できる。しかしながら、心臓ＣＴにおいては、ＣＴ画像に含まれるプラーク領域あるいは血栓領域の大きさに基づく指標しか低侵襲に計測することしかできない。ＣＴ画像に基づいて血栓前後の圧力差などを構造流体解析により計測できれば、血栓（あるいはプラーク）が及ぼす影響の定量化が期待できる。

冠循環の動態評価として臨床的には、超高速ＣＴ、シネアンギオグラム、超音波法、ＳＰＥＣＴ（シングルフォトンエミッショントモグラフィ）やＰＥＴ（ポジトロンエミッショントモグラフィ）を含む核医学イメージング、ＭＲＩ（核磁気共鳴画像法）などが開発導入され、診断や治療法の評価に役立っている。

しかしながら、冠微小血管を医用画像診断装置で正確に捉えることは困難である。また、血管形状が鮮明であっても医用画像にノイズが含まれる場合や、生体組織の境界の閾値設定に曖昧性が存在する場合も多い。このように、医用画像診断装置から得られる血管形状は不確定性を有している。

臨床応用で活用する場合、冠微小血管より上流の大動脈起始部から冠動脈の太い領域のみを対象として解析が行われる場合も多い。冠動脈の血流が冠微小血管の緊張性（トーヌス）にも大きく影響をうけるため、太い領域の冠動脈の出口における流量あるいは圧力もしくはそれらの変化率といった流体解析の境界条件を適切に設定することが課題となる。また、冠動脈の血流は、心臓の拍動による機械的因子（拍動による全体的な動き、局所的な伸縮、ねじり、せん断変形による強制変位あるいは外力）を受ける。流体解析のみでは、心臓の拍動等の機械的因子の影響を考慮できないため、血流の流量分布や内圧分布を精確に計測することができない。一方、画像で捉えられる心臓および血管系を対象とし、機械的因子の影響を考慮した構造-流体連成解析も実施されている。しかし、構造−流体連成解析を行う場合でも、血液（造影剤を含む）の流体解析における血管の入口や出口の境界条件や血管やプラークの材料モデルを正しく設定することが困難な場合も多い。また、画像に描出されない微小血管が存在する場合には、微小血管が血流に与える影響を考慮できない場合もある。そのため、構造-流体連成解析の解析結果は、実際の血流や血管変形を再現できていない恐れがある。また、境界条件や負荷条件、材料モデルが適切でない場合や、血管が大きな動きを伴う場合、収束性や解析安定性に問題がある場合もある。このように、従来の血管の構造流体解析は、多大な解析リソースと解析時間とが必要となる場合や、解析結果の誤差が大きくなる場合があり、現実的に臨床の現場で活用するのに問題が生じる場合がある。

特開２００８−２４１４３２号公報

門岡ら（ITUジャーナル，心臓シミュレータが拓くテーラーメード医療〜世界最先端の心臓シミュレータとその適用例のご紹介〜，Vol.41,No.6）

実施形態の目的は、血管（血液を含む）の構造流体解析および画像解析・追尾処理に関する精度の向上を可能とする血管解析装置、医用画像診断装置、血管解析方法、及び血管解析プログラムを提供することにある。

本実施形態に係る血管解析装置は、被検体の血管（血液を含む）に関する時系列の医用画像のデータを記憶する記憶部と、前記時系列の医用画像に含まれる血管領域に解析対象領域を設定し、前記解析対象領域に潜在変数同定領域を設定する設定部と、前記時系列の医用画像を画像処理して前記解析対象領域の時系列の形態指標と時系列の形状変形指標とを算出する算出部と、前記時系列の形態指標と前記時系列の形状変形指標と前記時系列の医用画像とに基づいて前記解析対象領域の構造流体解析に関する力学モデルを暫定的に構築する構築部と、前記暫定的に構築された力学モデルに基づく血管形態指標の予測値及び血液流量指標の予測値の少なくとも一方が予め計測された血管形態指標の観測値及び血液流量指標の観測値の少なくとも一方に整合するように前記潜在変数同定領域に関する潜在変数を同定する同定部と、を具備する。ここで、潜在変数とは、材料特性をあらわす材料モデルパラメータ，血管内腔に負荷される内圧分布などの負荷条件パラメータ、構造解析や流体解析の境界条件パラメータ、時系列の形態指標や形状変形指標の不確定性に関連する血管組織や血液および造影剤の空間座標あるいは幾何学的構造の不確定性に関連するばらつき分布パラメータあるいは医用画像データに内在するばらつき分布パラメータ（補足：各ＣＴ値のばらつき分布や生体組織の境界閾値のばらつき分布など）、の少なくともひとつを意味する。

本実施形態に係る医用画像診断装置（Ｘ線コンピュータ断層撮影装置）の概略的なブロック構成を示す図である。 本実施形態に係る構造流体解析の対象領域に関する力学モデルの一例を示す図。 図１のシステム制御部の制御のもとに行われる構造流体解析処理の典型的な流れを示す図。 図１の画像処理装置のブロック構成を示す図。 血管の芯線に直交する断面を模式的に示す図。 図４の画像解析・追尾処理による画像追尾処理に利用される血管芯線の形態の時間的変化を示す図。 図４の画像解析・追尾処理による画像追尾処理を説明するための図であり、時刻ｔと時刻ｔ＋Δｔとの間における追尾処理の一例を示す図。 図４の画像解析・追尾処理による血管芯線の曲げ変形や回転変位の算出例を示す図。 図４の力学モデル構築部により構築される形状モデルの芯線に直交する断面を示す図。 図４の力学モデル構築部による形状モデルへの強制変位履歴の割り当てを説明するための図。 図４の力学モデル構築部による形状モデルへの強制変位履歴の他の割り当て方法を示す図であり、血管形状変形指標がねじれの場合における割り当て例を示す図。 図４の力学モデル構築部による形状モデルへの強制変位履歴の他の割り当て方法を示す図であり、血管形状変形指標が曲げの場合における割り当て例を示す。 図４の統計的同定部による、階層ベイズモデル及びマルコフ連鎖モンテカルロ法による負荷条件（血管内の平均圧力）に関する事後分布算出と平均内圧の同定とを説明するための図。 図４の統計的同定部による、階層ベイズモデル及びマルコフ連鎖モンテカルロ法による材料モデルパラメータに関する事後分布算出と材料モデルパラメータ（血管壁の等価弾性率）の同定とを説明するための図。 図１の表示機器による力学的指標の一つである内圧の空間分布の表示例を示す図。 図１の表示機器による血液流量指標の一つである流速値の空間分布の表示例を示す図。 図１の表示機器により表示される左冠動脈起始部の血圧に関するグラフ。 図１の表示機器により表示されるＬＣＸとＬＤＡとの分岐点付近の血圧に関するグラフ。 図１の表示機器により表示される芯線方向に関する血圧変化に関するグラフ。 図４の力学モデル構築部による強制変位履歴の他の割り当て例を示す図。 図４の力学モデル構築部による強制変位履歴の他の割り当て例を示す図。 図４の力学モデル構築部による強制変位履歴の他の割り当て例を示す図。 本実施形態に係る繊維グループの変形を示す図。 本実施形態に係る肉厚円筒の力学モデルの直交断面を示す図。 図２０の微小扇形要素の拡大図。

以下、図面を参照しながら本実施形態に係わる血管解析装置、医用画像診断装置、血管解析方法、及び血管解析プログラムを説明する。

本実施形態に係る血管解析装置は、医用画像診断装置により発生された医用画像に含まれる血管領域を構造流体解析するためのコンピュータ装置である。本実施形態に係る血管解析装置は、医用画像診断装置に組み込まれていても良いし、医用画像診断装置とは別体のワークステーション等のコンピュータ装置であっても良い。以下、説明を具体的に行うため本実施形態に係る血管解析装置は、医用画像診断装置に組み込まれているものとする。

本実施形態に係る医用画像診断装置は、被検体をスキャンするための撮像機構を装備する如何なる種類の画像診断装置にも適用可能である。本実施形態に係る医用画像診断装置としては、例えば、Ｘ線コンピュータ断層撮影装置（Ｘ線ＣＴ装置）、磁気共鳴診断装置、超音波診断装置、ＳＰＥＣＴ装置、ＰＥＴ装置、放射線治療装置等が適宜利用可能である。以下、説明を具体的に行うため本実施形態に係る医用画像診断装置は、Ｘ線コンピュータ断層撮影装置であるものとする。

図１は、本実施形態に係る医用画像診断装置（Ｘ線コンピュータ断層撮影装置）の概略的なブロック構成図である。図１に示すように、Ｘ線コンピュータ断層撮影装置は、ＣＴ架台１０とコンソール２０とを有する。ＣＴ架台１０は、コンソール２０の架台制御部２３からの制御に従ってＸ線で被検体の撮像部位をスキャンする。撮像部位は、例えば、心臓である。ＣＴ架台１０は、Ｘ線管１１、Ｘ線検出器１３、及びデータ収集装置１５を有している。Ｘ線管１１とＸ線検出器１３とは、回転軸Ｚ回りに回転可能にＣＴ架台１０に装備されている。Ｘ線管１１は、造影剤が注入された被検体にＸ線を照射する。Ｘ線検出器１３は、Ｘ線管１１から発生され被検体を透過したＸ線を検出し、検出されたＸ線の強度に応じた電気信号を発生する。データ収集装置１５は、Ｘ線検出器１３から電気信号を読み出してデジタルデータに変換する。１ビュー毎のデジタルデータのセットは、生データセットと呼ばれている。複数のスキャン時刻に関する時系列の生データセットは、非接触データ伝送装置（図示しない）によりコンソール２０に伝送される。

コンソール２０は、システム制御部２１を中枢として、架台制御部２３、再構成装置２５、画像処理装置２７、入力機器２９、表示機器３１、及び記憶装置３３を有している。

架台制御部２３は、ユーザにより入力機器２９を介して設定されたスキャン条件に応じてコンソール２０内の各装置を制御する。

再構成装置２５は、生データセットに基づいて被検体に関するＣＴ画像のデータを発生する。具体的には、まず、再構成装置２５は、生データセットに前処理を施して投影データセットを発生する。前処理としては、対数変換や不均一補正、キャリブレーション補正等が含まれる。次に、再構成装置２５は、投影データセットに画像再構成処理を施してＣＴ画像のデータを発生する。画像再構成アルゴリズムとしては、ＦＢＰ（filtered back projection）法等の解析学的画像再構成法や、ＭＬ−ＥＭ（maximum likelihood expectation maximization）法やＯＳ−ＥＭ（ordered subset expectation maximization）法等の逐次近似画像再構成等の既存のアルゴリズムが採用可能である。本実施形態において再構成装置２５は、時系列の投影データセットに基づいて時系列のＣＴ画像のデータを発生する。ＣＴ画像は、造影剤により造影された血管に関する画素領域（以下、血管領域と呼ぶことにする。）を含んでいる。なお、ＣＴ画像は、ＣＴ値の２次元空間分布を表現するスライスデータであっても良いし、ＣＴ値の３次元空間分布を表現するボリュームデータであっても良い。以下、ＣＴ画像はボリュームデータであるとする。時系列のＣＴ画像のデータは、記憶装置３３に記憶される。

画像処理装置２７は、時系列のＣＴ画像に基づいて力学モデルを構築して構造流体解析を実行する。画像処理装置２７の処理の詳細については後述する。

入力機器２９は、ユーザからの各種指令や情報入力を受け付ける。入力機器２９としては、キーボードやマウス、スイッチ等が利用可能である。

表示機器３１は、ＣＴ画像や構造流体解析結果等の種々の情報を表示する。表示機器３１としては、例えばＣＲＴディスプレイや、液晶ディスプレイ、有機ＥＬディスプレイ、プラズマディスプレイ等が適宜利用可能である。

記憶装置３３は、ハードディスク装置等の種々の記憶媒体により構成される。記憶装置３３は、時系列の投影データや時系列のＣＴ画像データ等の種々のデータを記憶する。例えば、記憶装置３３は、時系列のＣＴ画像データをＤＩＣＯＭ（digital imaging and communications in medicine）規格に準拠した医用画像ファイル形式で記憶する。また、記憶装置３３は、外部機器により収集された医用データを時系列のＣＴ画像データに医用画像ファイル内において関連付けて記憶しても良い。

システム制御部２１は、中央演算処理装置（ＣＰＵ：central processing unit）や読み出し専用メモリ（ＲＯＭ：read only memory）、ランダムアクセスメモリ（ＲＡＭ：random access memory）を有する。システム制御部２１は、Ｘ線コンピュータ断層撮影装置の中枢として機能する。システム制御部２１は、ＲＯＭやＲＡＭに記憶されている血管解析プログラムを実行して本実施形態に係る血管構造解析処理を実行する。

なお、システム制御部２１、画像処理装置２７、入力機器２９、表示機器３１、及び記憶装置３３は、血管解析装置５０を構成する。本実施形態のように血管解析装置５０は、医用画像診断装置（Ｘ線コンピュータ断層撮影装置）に組み込まれていても良いし、医用画像診断装置とは別体のコンピュータ装置であっても良い。血管解析装置５０が医用画像診断装置とは別体の場合、血管解析装置５０は、医用画像診断装置やＰＡＣＳ（picture archiving and communication systems）からネットワークを介して時系列のＣＴ画像等の医用データを収集すれば良い。

次に、本実施形態に係る動作例を詳細に説明する。なお、本実施形態に係る血管解析装置、医用画像診断装置、血管解析方法、及び血管解析プログラムは、心臓血管や頸動脈、脳動脈等の人体のあらゆる部位の血管を解析対象とすることができる。しかしながら、以下、説明を具体的に行うため本実施形態に係る解析対象は、心臓の血管であるとする。

心臓の血管としては、例えば、冠動脈と大動脈とが挙げられる。冠動脈は、大動脈の冠動脈起始部から始まり心筋表面を走行し、心外膜側から内膜側に入り込む。冠動脈は、心筋の内膜において無数の毛細管に分岐する。分岐後、無数の毛細管は、再び統合して大心静脈を形成し、冠静脈洞に接続する。冠血管系は、他の臓器と異なり、心筋の収縮及び弛緩という力学的変化のなかで、灌流が保障されなければならないという点で特徴的である。

冠血流の特徴としては、心筋収縮による機械的血流阻害作用で冠動脈起始部の内圧が高くなる収縮期よりも、左心室拡張期に灌流圧が低下したときに多く流れることである。そのため、正常の冠動脈血流速波形は収縮期と拡張期の二峰性であり、拡張期血流が優位である。肥大型心筋症や大動脈弁狭窄症では収縮期に逆行性波を認め、大動脈逆流症では収縮期順行波が大きくなるなど疾患によって特異的な血流波形を呈することが知られている。また、拡張期の順行性波形は左室拡張機能、特に左室弛緩と密接な関係がある。左室弛緩遅延例では拡張期波形のピークが後ろにずれ、また減速脚がゆるやかになる傾向がある。また、このような症例では、頻拍時には拡張期の冠血流は十分に増大できず、心筋虚血を助長すると考えられる。

解剖学的に大動脈起始部から分岐する左右冠動脈に大動脈圧に等しい冠灌流圧（すなわち、冠動脈が分枝する大動脈起始部の圧力）がかかることにより、冠血流が生じる。冠血流を決定するのは大動脈圧である駆動圧とともに冠血管抵抗が重要である。１４０〜１８０μｍ以上の太い冠血管には冠血管低抗の２０％程度が存在するのに対し、１００〜１５０μｍ以下の微小血管には抵抗成分の残りの多くが存在するといわれる。従って、いわゆる冠狭窄などのない場合には抵抗値は冠微小血管の緊張性（トーヌス）に左右される。

血管抵抗因子としては、血管特性、動脈硬化、管狭窄、血液粘性、機械的因子があげられる。冠微小血管のトーヌスは、血管特性、心筋代謝（心筋酸素消費）、神経体液性因子、機械的因子、体液因子としての各種の血管作動性物質、血液粘性に規定され、さらに、心肥大、冠動脈硬化などを含めた様々な病変によっても影響され冠循環障害を起こす。

冠動脈血流拍動は、冠動脈血流の拍動パターン、心筋収縮による心筋内血流の制御、機械的刺激に対する心筋内血管の反応に影響される。心筋収縮が血流を阻害する機序としては、心筋内圧の上昇、心筋内血管容量の変化、心筋内血管の圧迫が挙げられる。心筋拡張期の血流規定因子には、拡張期の冠動脈圧、拡張期の血管外力、心拍数、心周期に占める拡張期の割合、心筋弛緩が存在する。

本実施形態に係る血管解析装置５０は、時系列のＣＴ画像に基づいて力学モデルを構築し、力学モデルを利用して心臓の血管についての構造流体解析を実行し、血管内の力学的指標や血管流量指標を高精度に算出する。精度の高い力学的指標や血管流量指標を算出するためには、力学モデルに精度の高い潜在変数を割り当てる必要がある。血管解析装置５０は、力学モデルを構築する際、初期的な力学モデルに逆解析を施して潜在変数を統計的に同定する。これにより血管解析装置５０は、高精度の潜在変数を決定することができる。力学的指標は、血管壁や血液に関する力学的な指標を意味する。血管壁に関する力学的指標としては、例えば、血管壁の変位に関する指標、血管壁に生じる応力やひずみに関する指標、血管内腔に負荷される内圧分布に関する指標、血管の硬さなどをあらわす材料特性に関する指標等に分類される。血管の硬さなどをあらわす材料特性に関する指標は、血管組織の応力とひずみの関係を示す曲線の平均的な傾き等が挙げられる。血液に関する力学的な指標における血管流量指標は、血管を流れる血液に関する血行動態の指標を意味する。血管流量指標としては、例えば、血液の流量、血液の流速、血液の粘性等が挙げられる。なお、潜在変数は、例えば、血管の材料構成式、あるいは血液の材料構成式といった材料モデルのパラメータ（例えばヤング率やポアソン比など）、血管内腔に負荷される内圧分布などの負荷条件パラメータ、構造解析や流体解析の境界条件パラメータ、時系列の形態指標や形状変形指標の不確定性に関連するばらつき分布パラメータの少なくとも一つを含む。ここで、時系列の形態指標や形状変形指標の不確定性に関連するばらつき分布パラメータとは、医用画像データには、各ＣＴ値のノイズに起因したばらつき分布や、生体組織の境界閾値の曖昧性に起因する確率分布などが存在することを考慮して、血管組織や血液の境界座標および特徴点（血管分岐部や造影剤分散配置など）の空間座標における不確定性、あるいは幾何学的構造パラメータ（芯線に垂直な断面の内腔半径など）の不確定性、あるいは医用画像データ自体（ＣＴ値や境界閾値など）の不確定性を、確率分布として表現したものである。

力学モデルは、血管や血液の挙動を表現するための数値モデルである。力学モデルは、構造流体解析の手法に応じて異なるタイプを有している。例えば、力学モデルは、連続体力学モデルと簡易的力学モデルとに分類される。連続体力学モデルは、例えば、有限要素法（ＦＥＭ：finite element method）や境界要素法に用いられる。簡易的力学モデルは、例えば、材料力学に基づく材料力学モデルと流れ学に基づく流体力学モデルとに分類される。なお、以下の説明において特に言及しない場合、力学モデルのタイプについては特に限定しないものとする。初期的な力学モデルは、潜在変数の確率分布や変数範囲から得られる潜在変数のパラメータに関するサンプリング集合（各パラメータの組み合わせの集合）が割り当てられた力学モデルを意味するものとする。

図２は、構造流体解析の対象領域（以下、解析対象領域と呼ぶ）に関する力学モデルＭ１の一例を示す図である。図２に示すように、力学モデルＭ１は、大動脈領域Ｒ１と右冠動脈領域Ｒ２と左冠動脈領域Ｒ３とを有している。血液は、大動脈から右冠動脈または左冠動脈へ流れる。

図２に示すように、力学モデルＭ１の大動脈起始部側の末端が血流の入口に設定され、右冠動脈領域の末端と左冠動脈領域の末端とが血流の出口に設定される。入口と出口との各々に境界条件が設定される。入口に関する境界条件は、例えば、入口における血流の流速あるいは血流による圧力、もしくはそれらの変化率を含む。出口に関する境界条件は、例えば、出口における血流の流速あるいは血流による圧力、もしくはそれらの変化率を含む。大動脈、右冠動脈、及び左冠動脈の変形は、血流に起因する血管壁への力学的作用、心臓の拍動による血管壁への力学的作用（外力）、血管断面境界の負荷条件、血管壁の材料モデル、血管の無応力状態、及び血管壁の幾何学的形状等の様々な要因に依存する。ここで、血流に起因する血管壁への力学的作用は、例えば、血流に起因する内圧と血流に起因するせん断力とを含む。血流に起因する内圧により、血管半径方向あるは血管内腔面の垂直方向に変形する。心臓の拍動による血管壁への力学的作用と、血流に起因するせん断力は、血管芯線方向に関する伸縮、ねじり、曲げといった血管壁への力学的作用による変形を生じさせる。血管芯線方向に関する伸縮、ねじり、曲げといった血管の全体的および局所的な変形挙動は、大動脈領域Ｒ１、右冠動脈領域Ｒ２、及び左冠動脈領域Ｒ３、あるいはその一部、あるいはその周辺組織に、強制変位（移動ベクトルや回転変位）あるいは荷重ベクトルの時間的変化を与えることで、負荷条件を設定する。また、血流に起因する内圧に基づく血管半径方向あるいは内腔面に垂直な方向への変形は、血管内腔への圧力分布の時間的変化として与えることで、負荷条件を設定する。

大動脈領域Ｒ１、右冠動脈領域Ｒ２、及び左冠動脈領域Ｒ３には、構造流体解析において、強制変位による変位拘束条件が割り当てられる。これにより、構造流体解析における血管壁の変形自由度を縮小することができ、計算収束性の安定化や、解析時間の短縮を実現できる。

また、例えば、血管の形状の変形度合は、血管壁の材料に依存する。そのため、大動脈領域Ｒ１、右冠動脈領域Ｒ２、及び左冠動脈領域Ｒ３に材料モデルが割り当てられる。

また、例えば、血管の形状の変形度合は、血管の無応力状態に依存する。負荷条件の初期値として、血管の残留応力分布を負荷条件として与えても良い。ここで、血管解析対象の数値計算用力学モデルの空間を離散化した節点集合および節点から構成される要素集合において、材料モデルや境界条件や負荷条件などの解析条件を同定しない領域と、同定する領域にわけ、解析条件を同定しない領域内の節点には強制変位履歴の変位拘束条件を与え、材料モデルを同定する領域では、血管壁表面（外表面）の節点のみ強制変位履歴の変位拘束条件を与え、血管壁内部は変位自由度を確保し変位拘束を与えない。これにより、構造流体解析の変形自由度をおさえることができ、安定的かつ効率的に解析を行うことができる。ただし、血管壁表面上に緩衝用のダミーの要素集合を設け、その表面の節点に強制変位を与えてもよい。これにより、血管内腔にプラークなどによる突起がある場合や、血管分岐部など、内圧による負荷が、芯線方向の断面外の変形に影響を与える場合には、血管内腔だけではなく血管壁の形態指標も参照して血管に負荷される荷重ベクトルと内圧を分離して同定することができる。また、生理学的には壁表面の脂肪層を模擬しており、一方、数値計算上は、血管壁表面に強制変位を与えることで血管壁内部に局所的に実際とは異なる高い応力が発生することを避ける効果もある。

これら材料モデル、境界条件、及び負荷条件等の潜在変数に関するパラメータは、後述する力学モデルに基づく逆解析（統計的同定処理）により同定される。逆解析により同定された精確な潜在変数は、力学モデルに割り当てられる。精確な潜在変数が割り当てられた力学モデルにより、解析対象血管領域外の血管や心臓等の外部要因による当該解析対象血管領域への影響を加味した構造流体解析あるいは流体解析あるいは構造解析あるいは画像解析に基づく血行動態解析を実行することができる。血管解析装置５０は、力学モデルの構築に関し、逆解析による潜在変数の同定により、従来例に係る次の４点の困難を解決することができる。困難１．冠動脈の材料モデルの同定方法。困難２．心臓の形状の変形の冠動脈への影響の組み込み。困難３．冠動脈の境界条件の同定方法。困難４．医用画像データの不確定性に起因したばらつきを有する血管形状による画像解析や構造流体解析。この４点の困難の克服により、血管解析装置５０は、逆解析による潜在変数の同定を行わない従来の血管構造流体解析に比して、解析精度の向上を実現する。

次に、本実施形態に係る構造流体解析処理の詳細について説明する。図３は、本実施形態に係るシステム制御部２１の制御のもとに行われる構造流体解析処理の典型的な流れを示す図である。図４は、画像処理装置２７のブロック構成を示す図である。

図３に示すように、構造流体解析処理において、まず、システム制御部２１により記憶装置３３から処理対象の医用画像ファイルが読み出され、画像処理装置２７に供給される。医用画像ファイルは、時系列のＣＴ画像のデータの他に、当該被検体に関する血管内腔に関する圧力値のデータ、血液流量指標の観測値のデータ、及びプラーク指標を含んでいる。ＣＴ画像以外ではＭＲＩ画像や超音波エコー画像であってもよい。時系列のＣＴ画像のデータは、時系列のＣＴ値の３次元空間分布を表現するデータである。時系列のＣＴ画像は、例えば、１心拍で約２０枚、すなわち、約２０心位相分のＣＴ画像を含んでいる。

図３に示すように、システム制御部２１は、画像処理装置２７に領域設定処理を行わせる（ステップＳ１）。ステップＳ１において画像処理装置２７の領域設定部５１は、時系列のＣＴ画像に含まれる血管領域に構造流体解析の解析対象領域を設定する。解析対象領域は、冠動脈に関する血管領域の任意の一部分に設定される。例えば、領域設定部５１は、ユーザによる入力機器２９を介した指示、または、画像処理により血管領域に解析対象領域と同定対象領域を設定する。

ここで、図５を参照しながら、血管の構造について説明する。図５は、血管の芯線に直交する断面（以下、血管断面と呼ぶ）を模式的に示す図である。図５に示すように、血管は、管状の血管壁を有している。血管壁の中心軸は芯線と呼ばれている。血管壁の内側は内腔と呼ばれている。内腔に血液が流れる。内腔と血管壁との境は血管内壁と呼ばれている。血管壁の外側には心筋等の血管周辺組織が分布している。血管壁と血管周辺組織との境は血管外壁と呼ばれている。血管壁内部にプラークが発生することがある。図５に示すように、プラークは、例えば、石灰化した石灰化プラーク、粥状プラーク等に分類される。粥状プラークは、やわらかく、血管内壁が破裂して血栓として血管内部に染み出す危険性があり、不安定プラークと呼ばれることもある。従って、プラークの性状を把握することは臨床的に有用である。プラークの性状や存在領域は、医用画像ファイルに含まれるプラーク指標により特定可能である。プラーク指標は、例えば骨のＣＴ値を基準に正規化したＣＴ値の大きさなどにより相対的に判別することができる。しかしながら、血管内部のプラークの変形特性や硬さを解析するのは容易ではない。

ステップＳ１が行われるとシステム制御部２１は、画像処理装置２７に画像解析・追尾処理を行わせる（ステップＳ２）。ステップＳ２において画像処理装置２７の画像解析・追尾処理部５３は、時系列のＣＴ画像に画像処理を施して時系列の血管形態指標と時系列の血管形状変形指標とを算出する。具体的には、画像解析・追尾処理部５３は、時系列のＣＴ画像に画像解析処理を施すことにより時系列の血管形態指標を算出し、時系列のＣＴ画像に追尾処理を施すことにより時系列の血管形状変形指標を算出する。

より詳細には、画像解析・追尾処理部５３は、画像解析処理において、各ＣＴ画像から血管領域を抽出し、血管の内腔に関する画素領域（以下、血管内腔領域と呼ぶ）と血管壁に関する画素領域（以下、血管壁領域と呼ぶ）とを特定する。画像解析・追尾処理部５３は、血管形態指標として、血管の芯線に垂直な断面、あるいは血管内腔面に垂直な面が、血管内腔、血管壁、プラーク領域に交わる領域上の複数の画素の３次元座標を特定する。なお、血管形態指標は、３次元座標だけでなく、芯線に垂直な断面における一定角度ごとの血管内腔の半径や直径及び０°の方向ベクトル、あるいは断面における全角度に対する平均面積や平均半径、あるいは、芯線方向に垂直な複数の断面で囲まれた血管内腔容積、あるいは内腔面に垂直な複数断面で囲まれた血管壁容積やプラーク容積等の幾何学的指標でも良い。

追尾処理において、画像解析・追尾処理部５３は、ユーザからの入力機器２９を介した指示または画像処理により、血管領域や血液や造影剤やプロトンにおける特徴点や特徴形状、代表点、画素等の複数の追跡点を設定する。例えば、画像解析・追尾処理部５３は、血管分岐部や表面の特徴形状などの追跡点集合を設定する。各時刻（各心位相）における画像解析・追尾処理部５３の追尾処理により得られた追跡点集合の変位データから、力学モデルの血管壁表面あるいは血管壁内部あるいは血管内腔における節点の変位の時間的変化を補間処理などにより算出し、強制変位として与える。また、例えば、画像解析・追尾処理部５３は、力学モデルに血管芯線上の節点を定義する。画像解析・追尾処理部５３は、力学モデルの血管壁表面あるいは血管壁内部あるいは血管内腔における節点の変位の時間的変化から、血管の芯線方向に関する伸縮やねじりや曲げに関する変形を抽出し、血管芯線と芯線に垂直な断面における節点の強制変位として与えることで表現しても良い。このように、血管形状変形指標としては、力学モデルにおける各時刻の節点の強制変位データ（強制変位履歴）を特定する。

以下、図６及び図７を参照しながら、画像解析・追尾処理を説明する。図６は、血管芯線の形態の時間的変化を示す図である。図６に示すように、例えば、時系列の医用画像は、１心拍につき２０枚のＣＴ画像を含んでいるものとする。すなわち、心位相０％から９５％まで５％おきにＣＴ画像が得られているものとする。画像解析・追尾処理部５３により血管領域の芯線が抽出される。図６に示すように、芯線の形態は、心位相の経過に従って変化している。

図７は、時刻ｔと時刻ｔ＋Δｔとの間における追尾処理の一例を示す図である。図７に示すように、血管芯線上にＰ１からＰ１０の力学モデルにおける節点が設定されており、各断面上の血管の力学モデルの節点と力学的につながっている。ただし、血液の力学モデルの節点とは独立である。血管の追跡点の変位データをもとに、血管芯線上のＰ１からＰ２０の節点の変位データを補間などの処理により算出し、各節点に強制変位が設定されているものとする。血管形状変形指標と血管形態指標について説明するため、節点Ｐ１３と節点Ｐ１４とにより規定される局所血管領域ＲＡを考える。時刻ｔにおいて、芯線方向に関する節点Ｐ１３と節点Ｐ１４との間の距離がＬであり、血管領域の半径がｒであるとする。画像解析・追尾処理部５３から、節点Ｐ１３と節点Ｐ１４の血管芯線方向の伸縮やねじりや曲げといった強制変位を抽出することにより、節点Ｐ１３の強制変位（３次元空間における移動変位と、芯線方向の回転変位）と節点Ｐ１４の強制変位（３次元空間における移動変位と、芯線方向の回転変位）を算出する。

図８は、血管芯線の曲げ変形や回転変位の算出例を示す図である。図８に示すように、例えば、ねじり角は、線ａ及びｂから構成される面の法線方向ベクトルcの変化から算出しても良い。

図７や図８に示すように、画像解析・追尾処理部５３は、追跡点の座標と移動ベクトルとに基づいて、芯線上の各節点の強制変位（３次元空間における移動変位と、芯線方向の回転変位）を算出し、血管形状変形指標を算出する。例えば、画像解析・追尾処理部５３は、隣り合う２つの節点の座標差の時間変化を芯線方向に関する伸縮距離ΔＬとして算出する。また、画像解析・追尾処理部５３は、各芯線上の節点について、当該節点と当該節点を含む血管領域断面上の他の節点（血管内腔あるいは血管壁あるいはプラーク領域における節点）との間の距離の時間変化を半径方向に関する伸縮距離Δｒとして算出する。また、画像解析・追尾処理部５３は、各追跡点について、当該追跡点の近傍の複数の追跡点の座標と移動ベクトルとに基づいて、芯線上の当該節点の芯線方向のねじれ角度Δθを算出する。

また、画像解析・追尾処理部５３は、血液領域の造影剤やプロトンの画像追尾により、血液流量指標として、流速あるいは芯線方向断面の平均流速あるいは平均流量を算出してもよい。

血管形状変形指標は、力学モデルにおける強制変位として利用される。以下、時系列の血管形態指標を形状履歴と呼び、時系列の血管形状変形指標を強制変位履歴と呼ぶことにする。

ステップＳ２が行われるとシステム制御部２１は、画像処理装置２７に構築処理を行わせる（ステップＳ３）。ステップＳ３において画像処理装置２７の力学モデル構築部５５は、形状履歴（時系列の血管形態指標）と強制変位履歴（時系列の血管形状変形指標）と時系列の医用画像（ＣＴ画像，ＭＲＩ画像，超音波エコー画像などのＤＩＣＯＭデータ）とに基づいて解析対象領域に関する力学モデルを暫定的に構築する。力学モデルは、構造流体解析を行うための解析対象領域に関する数値モデルである。

以下、ステップＳ３について詳細に説明する。力学モデル構築部５５は、まず、医用画像と形状履歴とに基づいて、力学モデル（数理モデル）を解くための形状モデルを構築する。形状モデルは、各時刻における血管領域の幾何学的構造を模式的に表現したものである。形状モデルは、例えば、複数の離散化領域に区分されている。各離散化領域の頂点は、節点と呼ばれる。力学モデル構築部５５は、時刻毎の医用画像に含まれる血管領域と血管形態指標とに基づいて時刻毎の形状モデルを構築しても良いし、特定の時相の医用画像に含まれる血管領域と血管形態指標とに基づいて時刻毎の形状モデルを構築しても良い。また、例えば、初期の負荷状態として、解析対象領域に対応する血管に残留応力が存在しないと仮定する場合、無応力状態の時相として、解析対象領域に対応する血管が最も収縮した時相を無応力状態であると仮定する。

図９は、形状モデルの芯線に直交する断面を示す図である。図９に示すように、形状モデルは、芯線から外側に向けて血管内腔領域、血管壁領域を有している。プラークが存在する場合、血管壁領域にプラーク領域を設けても良い。また、血管周辺組織による血管への影響を考慮する場合、血管周辺組織のダミー要素を血管壁領域の外側に設けても良い。

形状モデルが構築されると力学モデル構築部５５は、各潜在変数の確率分布や変数範囲から得られる潜在変数のパラメータに関するサンプリング値（例えばマルコフ連鎖モンテカルロ法などによる、各パラメータの組み合わせの集合からのサンプリング）を力学モデルに設定する。例えば、力学モデル構築部５５は、図２に示すように、大動脈領域Ｒ１の大動脈起始部側の末端に入口に関する境界条件の同定対象の領域（以下、境界条件同定領域）を設定し、右冠動脈領域Ｒ２の末端と左冠動脈領域Ｒ３の末端とに出口に関する境界条件同定領域を設定する。力学モデル構築部５５は、各境界条件同定領域に境界条件の確率分布や変数範囲から得られる境界条件のパラメータに関するサンプリング値を割り当てる。また、力学モデル構築部５５は、大動脈領域Ｒ１、右冠動脈領域Ｒ２、及び左冠動脈領域Ｒ３に材料モデルの同定対象の領域（以下、材料モデル同定領域と呼ぶ）及び負荷条件の同定対象の領域（以下、負荷条件同定領域と呼ぶ）を設定する。力学モデル構築部５５は、各材料モデル同定領域に材料モデルの確率分布や変数範囲から得られる材料モデルのパラメータに関するサンプリング値を割り当て、各負荷条件同定領域に負荷条件の確率分布や変数範囲から得られる負荷条件のパラメータに関するサンプリング値を割り当てる。血管は、流量が０でも残留応力が存在すると言われている。例えば、力学モデル構築部５５は、流量が０の場合の残留応力を負荷条件の初期値として解析対象領域に割り当てても良い。また、力学モデル構築部５５は、幾何学的構造に不確定性がある部位に、幾何学的構造の同定対象の領域（以下、幾何学的構造同定領域と呼ぶ）を設定しても良い。なお、幾何学的構造のパラメータは、幾何学的構造の不確定性に関連するばらつき分布パラメータ、あるいは医用画像データに内在するばらつき分布パラメータであり、各ＣＴ値のばらつき分布や生体組織の境界閾値のばらつき分布などであってよい。詳細は、後述するが、力学モデル構築部５５は、プラーク領域に材料モデルを設定しても良い。材料モデルの詳細については後述する。

形状モデルが構築されると力学モデル構築部５５は、ステップＳ２において算出された時系列の血管形状変形指標、すなわち、強制変位履歴を形状モデルに割り当てる。潜在変数及び強制変位履歴が割り当てられた形状モデルを力学モデルと呼ぶことにする。

図１０は、形状モデルＭ２は血管や血液の力学モデルの一部を示しており、力学モデルの節点への強制変位履歴の割り当てを説明するための図である。図１０は、形状モデルＭ２の一部分を示している。ただし、図１０では芯線がＭ２内に位置する場合を示しているが、芯線がＭ２外に位置する場合でも良い。図１０に示すように、形状モデルＭ２には複数の節点ＰＮ（ＰＮ１，ＰＮ２）が設定されている。芯線上の節点をＰＮ１と称し、血管や血液の力学モデルにおける節点をＰＮ２と称することにする。形状モデルＭ２は、ダミー要素表面、血管外壁、血管内壁、プラーク領域表面、プラーク領域内部、あるいは血液部、に設定される。力学モデル構築部５５は、形状モデルＭ２の各節点ＰＮ１に強制変位、すなわち、血管形状変化指標を時刻毎に割り当てる。

具体的には、力学モデル構築部５５は、芯線上の隣り合う節点ＰＮ１と節点ＰＮ１とをビーム要素（あるいはリジッド要素）ＥＢで結ぶ。また、力学モデル構築部５５は、節点ＰＮ１と当該節点ＰＮ１を通る直交断面に含まれる他の節点ＰＮ２とをビーム要素ＥＢで結ぶ。力学モデル構築部５５は、節点ＰＮ１及びビーム要素ＥＢに各血管形状変形指標の形状変位方向に関する拘束条件を割り当てる。材料モデルや血管内腔の内圧を同定する領域では、強制変位としては、芯線方向に関する血管壁（あるいはダミー要素）表面の伸縮、血管壁（あるいはダミー要素）表面のねじれ、血管壁（あるいはダミー要素）表面の曲げ変形が挙げられる。例えば、材料モデルや血管内腔の内圧を同定しない領域では、芯線方向に関する強制変位だけではなく、半径方向に関する血管壁の時系列の伸縮（変位）も強制変位履歴として割り当てられる。また、血管内腔に突起がある場合や血管分岐部などのように、内圧が、芯線方向断面外への変形に寄与する場合は、その領域には強制変位を与えずに、その領域の周辺部（たとえばダミー要素の表面節点）のみに強制変位履歴を割り当てる。また、力学モデル構築部５５、節点ＰＮ１及びビーム要素ＥＢに時系列の血管形状変形指標を強制変位履歴として割り当てられる。これにより、血管の全体及び局所に関する伸縮変形やねじれ変形や曲げ変形が表現される。

なお、強制変位履歴の割当対象は、芯線上の節点及びビーム要素に限定されない。図１１及び図１２は、形状モデルへの強制変位履歴の他の割り当て方法を示す図である。図１１は、血管形状変形指標がねじれの場合における割り当て例を示し、図１２は、血管形状変形指標が曲げの場合における割り当て例を示している。図１１及び図１２に示すように、力学モデル構築部５５は、形状モデルの表面あるいは内部の節点に直接的に強制変位履歴を割り当てても良い。例えば、ステップＳ２において画像解析・追尾処理部５３は、特徴点についての伸縮量やねじり量や曲げ量等の血管形状変形指標を算出する。力学モデル構築部５５は、算出された血管形状変形指標を特徴点の周囲の節点へ補間（内挿や外挿）することにより直接的に割り当てる。

本実施形態に係る画像処理装置２７は、ステップＳ３において暫定的に構築された力学モデルを用いて逆解析を施し、力学モデルに設定される潜在変数を統計的に同定する。統計的同定処理は、後述のステップＳ６において行われる。ステップＳ４及びＳ５は、それぞれ統計的同定処理に用いる血管形態指標及び血液流量指標を算出するために設けられている。

ステップ３が行われるとシステム制御部２１は、画像処理装置に２７に血管応力解析処理を行わせる（ステップＳ４）。ステップＳ４において画像処理装置２７の血管応力解析部５７は、現段階の力学モデルに血管応力解析を施して時系列の血管形態指標の予測値を算出する。血管形態指標は、既述の血管形態指の何れであっても良いが、例えば、血管芯線方向に関する内腔領域の断面形状指標や血管壁の断面形状指標が用いられると良い。具体的には、内腔領域の断面形状指標は、内腔領域の注目画素の座標値、内腔領域の幾何学的構造パラメータ（内腔領域の半径、内腔領域の直径など）の少なくとも一つである。また、血管壁領域の断面形状指標は、具体的には、血管壁領域の注目画素の座標値、血管壁領域の幾何学的構造パラメータ（血管壁領域の半径、壁領域の直径など）の少なくとも一つである。なお、予測値は、力学モデルに血管応力解析を施して算出された血管形態指標の算出値を意味する。

また、ステップ３が行われるとシステム制御部２１は、画像処理装置２７に血液流体解析処理を行わせる（ステップＳ５）。ステップＳ５において画像処理装置２７の血液流体解析部５９は、暫定的に構築された力学モデルに血液流体解析を施して時系列の血液流量指標の予測値を算出する。血液流量指標は、血流量または流速、あるいはその空間的・時間的な平均値の少なくとも一つである。なお、予測値は、力学モデルに血液流体解析を施して算出された血液流体指標の算出値を意味する。

ステップＳ４及びＳ５が行われるとシステム制御部２１は、画像処理装置２７に同定処理を行わせる（ステップＳ６）。ステップＳ６において画像処理装置２７の統計的同定部６１は、ステップＳ４において算出された血管形態指標の予測値とステップＳ５において算出された血液流量指標の予測値とが、事前に収集された血管形態指標の観測値と血液流量指標の観測値とに整合するように力学的モデルの潜在変数のパラメータを統計的に同定する。

図４に示すように、統計的同定部６１は、第１統計的同定部６１−１と第２統計的同定部６１−２とを有している。第１統計的同定部６１−１は、血管形態指標の予測値が血管形態指標の観測値に整合するように潜在変数のパラメータを統計的に同定する。第２統計的同定部６１−２は、血液流量指標の予測値が血液流量指標の観測値に整合するように潜在変数のパラメータを統計的に同定する。以下、第１統計的同定部６１−１と第２統計的同定部６１−２とを順番に説明する。

具体的には、ステップＳ６において第１統計的同定部６１−１は、ステップＳ４において算出された血管形態指標の予測値と観測値とに基づくデータ分布を設定する。データ分布は、例えば、血管形態指標の予測値と観測値との誤差に関する多変量正規分布関数を示す。具体的には、力学モデルにおける各節点あるいは各要素での予測値と観測値の誤差に関する正規分布関数値を算出し、それらの積をデータ分布として統計的同定処理部において設定する。データ分布は、時刻毎に個別に設定されても良いし、複数時刻まとめて設定されても良い。次に、第１統計的同定部６１−１は、力学モデルの潜在変数に事前分布（事前確率分布）を割り当てる。具体的には、材料モデル、境界条件、負荷条件、及び時系列の形態指標や形状変形指標の不確定性に関連するパラメータの各々に事前分布が割り当てられる。例えば、負荷条件のパラメータの一つである血管内腔に関する圧力値に関する事前分布が割り当てられる。圧力値の取り得る値の範囲（想定範囲）は、経験的に予め限定することができる。第１統計的同定部６１−１は、これら想定範囲内に限定して内圧値に関するモンテカルロシミュレーションを実行することにより、各離散化領域について内圧値の確率分布、すなわち、事前分布を算出する。また、力学モデル構築部５５は、事前分布として、芯線方向の圧力分布は滑らかであること、また、時間経過に伴う圧力変化が滑らかであること、血流の逆流がないことが観測されている場合には芯線方向の平均的な圧力変化の傾きが負であることを、例えば、多変量正規分布関数により数学的に表現した確率分布を事前分布として設定しても良い。これら想定範囲内に限定して、想定した確率分布に従って、負荷条件のパラメータに関するモンテカルロシミュレーションを実行することができ、力学モデルへ設定するための負荷条件（潜在変数）のサンプリング値を得ることができる。次に、第１統計的同定部６１−１は、各潜在変数について、事前分布とデータ分布とに統計的同定処理を施すことにより事後分布（事後確率分布）を算出する。統計的同定処理は、例えば、階層ベイズモデルやマルコフ連鎖モデルが挙げられる。そして、第１統計的同定部６１−１は、各潜在変数について、事後分布の最頻値や平均値等の統計値から各潜在変数のパラメータを同定する。例えば、上述の例の場合、血管内腔圧力値に関する事後分布が算出され、この事後分布から血管内腔圧力値の同定値が算出される。

図１３は、階層ベイズモデル及びマルコフ連鎖モンテカルロ法による負荷条件（血管内の平均圧力）に関する事後分布算出と平均内圧の同定とを説明するための図である。図１３に示すように、血管起始部から延びる血管に石灰化プラーク領域と粥状プラーク領域とが設定されているものとする。石灰化プラーク領域は材料モデル同定領域に設定され、粥状プラーク領域に材料モデル同定領域に設定される。血管起始部から血管芯線方向に沿って進むにつれて血管内圧が降下する。血管芯線に沿って複数の節点が設定される。各節点を含む直交断面（節点断面）において内腔内圧の事後分布が算出され、事後分布の最頻値が同定される。

なお、血管形態指標の観測値としては、例えば、ステップＳ２において算出された血管形態指標が用いられる。

第２統計的同定部６１−２による処理は、データ分布の算出に用いる指標が異なるだけで第１統計的同定部６１−１による処理と同様である。すなわち、第２統計的同定部６１−２は、まず、ステップＳ５において算出された血液流量指標の予測値と観測値とに基づくデータ分布を設定する。次に、第２統計的同定部６１−２は、力学モデルの潜在変数に事前分布を割り当てる。例えば、血管に関する材料モデルのパラメータや血液に関する材料モデルのパラメータ、プラークに関する材料モデルのパラメータに関する事前分布が割り当てられる。これら材料モデルのパラメータとしては、例えば、弾性率などの材料モデルパラメータや、血液の構成式における粘性に関するパラメータが挙げられる。材料モデルのパラメータの想定範囲や確率分布は、経験的に予め設定することができる。第２統計的同定部６１−２は、各離散化領域について材料モデルのパラメータの確率分布、すなわち、事前分布を設定し、これら想定範囲内に限定して、想定した確率分布に従って、材料モデルのパラメータに関するモンテカルロシミュレーションを実行することができ、力学モデルへ設定するための材料モデルパラメータ（潜在変数）のサンプリング値を得ることができる。次に、第２統計的同定部６１−２は、各潜在変数について、事前分布とデータ分布とに統計的同定処理を施すことにより事後分布を算出し、算出された事後分布の統計値から各潜在変数のパラメータを同定する。例えば、上述の例の場合、材料モデルのパラメータに関する事後分布が算出され、この事後分布から材料モデルのパラメータの同定値が算出される。

図１４は、階層ベイズモデル及びマルコフ連鎖モンテカルロ法による材料モデルパラメータに関する事後分布算出と材料モデルパラメータ（血管壁の等価弾性率）の同定とを説明するための図である。図１４に示すように、血管モデルは、図１３と同様であるとする。材料モデル同定領域に限定して、血管壁の材料モデルのパラメータ（例えば、等価弾性率）の事後分布が算出され、事後分布の最頻値が同定される。

なお、血液流量指標の観測値は、例えば、大動脈に送り出される血流量変化であると仮定し、血管形態指標の観測値を、時系列のＣＴ画像から画像処理により計測される左心室の容積変化値（ＣＦＡ）として用いることができる。造影剤の冠動脈内注入後の造影剤の画像追尾により特徴点の移動量の時間的変化を算出することにより、流速や流量を算出してもよい。また、造影剤の血管芯線方向あるいは時間的な特定領域の濃度変化量を取得し、その濃度変化を各領域の芯線方向距離間隔で除した値や、濃度変化の時間的変化率から、流速や流量を算出してもよい。ＭＲＩの場合はプロトンの画像追尾を用い、超音波エコーの場合には、コントラストエコー図法などにより流量を算出する。

また、解析対象領域の各画素の座標値が確定値であることを前提としない場合、すなわち、解析対象領域の幾何学的構造に不確定性があると仮定する場合、潜在変数に幾何学的構造を含めても良い。この場合、統計的同定部６１は、各節点の座標値の芯線方向に関する所定範囲内の変動、あるいは、解析対象領域の径の所定範囲内の変動を表現した正規分布などの確率分布を事前分布として設定すると良い。この場合、解析対象領域の形状が滑らかであること、また、芯線における節点の順番が不変であるという制約も事前分布として設定しても良い。これら想定範囲内に限定して、想定した確率分布に従って、幾何学的構造のパラメータに関するモンテカルロシミュレーションを実行することができ、力学モデルへ設定するための幾何学的構造の不確定性パラメータ（潜在変数）のサンプリング値を得ることができる。

なお、各ステップＳ６において第１統計的同定部６１−１による統計的同定処理と第２統計的同定部６１−２による統計的同定処理との両方が行われなくても良い。すなわち、各ステップＳ６においては、第１統計的同定部６１−１による統計的同定処理と第２統計的同定部６１−２による統計的同定処理との何れか一方が行われても良い。

また、上記の例においては、第１統計的同定部６１−１は、血管形態指標の予測値が血管形態指標の観測値に整合するように潜在変数のパラメータを統計的に同定し、第２統計的同定部６１−２は、血液流量指標の予測値が血液流量指標の観測値に整合するように潜在変数のパラメータを統計的に同定するとした。しかしながら、統計的同定部６１は、構造−流体連成解析に基づいて、血管形態指標の予測値と血液流量指標の予測値とが血管形態指標の観測値と血液流量指標の観測値とに整合するように潜在変数のパラメータを統計的に同定しても良い。統計的同定部６１による統計的同定処理のさらなる詳細については後述する。

ステップＳ６が行われるとシステム制御部２１は、画像処理装置２７に設定処理を行わせる（ステップＳ７）。ステップＳ７において画像処理装置２７の力学モデル構築部５５は、ステップＳ６において算出された潜在変数のパラメータを力学モデルに設定する。

ステップＳ７が行われるとシステム制御部２１は、同定終了条件が満たされたか否かを判定する（ステップＳ８）。ステップＳ８において同定終了条件が満たされていないと判定した場合（ステップＳ８：ＮＯ）、システム制御部２１は、ステップＳ４、Ｓ５、Ｓ６、Ｓ７、及びＳ８を繰り返す。ここで、同定終了条件は、同定終了を判定するための指標（以下、同定終了指標と呼ぶ）が規定値に達するか否かにより表現される。同定終了指標としては、例えば、血管形態指標の予測値と観測値との差分値が挙げられる。この場合、システム制御部２１は、この差分値が既定値よりも大きい場合、同定終了条件が満たされていないと判定し、差分値が既定値よりも小さい場合、同定終了条件が満たされたと判定する。また、同定終了指標は、例えば、モンテカルロ法のサンプリング点の数でも良い。この場合、システム制御部２１は、このサンプリング点の数が既定値よりも小さい場合、同定終了条件が満たされていないと判定し、サンプリング点の数が既定値よりも大きい場合、同定終了条件が満たされたと判定する。同定終了条件が満たされた場合、力学モデル構築部５５は、その時点の最新の力学モデルを最終的な力学モデルに設定する。

最終的な力学モデルが構築されると、力学モデル構築部５５は、血管形状変形指標の観測値、最終的な力学モデルにおける負荷条件のパラメータ、及び材料モデルのパラメータを関連付けたモデル（以下、関連モデルと呼ぶ）を算出する。関連モデルは、記憶装置３３に記憶される。関連モデルは、検索の容易性等のため患者情報や検査情報等に関連付けて記憶されると良い。なお、血管形態指標や血液流量指標の観測値、最終的な力学モデルにおける負荷条件のパラメータ、及び材料モデルのパラメータは、必ずしもモデルの形態で関連付けられる必要はなく、例えば、テーブルあるいはデータベースであっても良い。

上記のステップＳ４、Ｓ５、Ｓ６、Ｓ７、及びＳ８は、同一の同定法で反復しても良いし、異なる同定法で反復しても良い。異なる同定法で反復する場合、例えば、まず、簡易的力学モデルを利用して潜在変数を暫定的に同定し、次に、連続体力学モデルを利用して潜在変数を正確に同定しても良い。このように統計的同定処理を異なる手法で２段階に分けて行うことにより、潜在変数のパラメータを短時間で収束させることができる。簡易的力学モデルを利用する方法としては、内圧及び外圧を厚肉円筒の材料力学の式と、ハーゲン・ポアズイユ流れ及び修正ベルヌーイの式とを用いる方法が挙げられる。連続体力学モデルを利用する方法としては、ＦＥＭ構造流体解析が挙げられる。簡易的力学モデルを利用する同定法と連続体力学モデルを利用する同定法との詳細については後述する。

ステップＳ８において同定終了受件が満たされたと判定した場合（ステップＳ８：ＹＥＳ）、システム制御部２１は、画像解析・追尾処理部５３に修正処理を行わせても良い（ステップＳ９）。ステップＳ９において画像解析・追尾処理部５３は、統計的同定法による逆解析で得られた潜在変数のもとで実施した構造流体解析結果（力学的指標の予測値及び血液流体指標の予測値）が観測値（力学的指標の観測値及び血液流体指標の観測値）に整合するように、時系列の医用画像に含まれる血管領域の形状を修正しても良い。表示機器３１は、修正後の時系列の医用画像に基づく診断結果を表示する。これにより、血管解析装置５０は、最終的な力学モデルを考慮した診断結果を表示することができる。あるいは、表示機器３１は、逆解析による同定とその構造流体解析により観測結果とが整合しない血管箇所・領域を画面に表示しても良い。例えば、血管の挙動の動きの速い心位相の画像はぼやけることが多く、医用画像をもとにした画像解析により観測した血管形状には誤差が大きい箇所や領域が存在する。比較的、血管の挙動が安定した心位相のデータでは、ノイズが少ない画像から得られる。そのような誤差分布の小さい血管形状データに基づき、誤差が大きな心位相での血管形状を、力学モデルを用いることで、正しく内挿することができ、そのような誤差が大きい血管箇所や領域について、正しく内挿した形状とともに、元データからのばらつき分布とともに、表示することができる。これにより、血管形状表示の安定性を確保できるとともに、形状の不確定性をユーザが認識できる。

ステップＳ９が行われるとシステム制御部２１は、画像処理装置２７に血管応力解析処理を行わせる（ステップＳ１０）。ステップＳ１０において画像処理装置２７の血管応力解析部５７は、最終的な力学モデルに血管応力解析を施して時系列の力学的指標の予測値の空間分布を算出する。具体的には、離散化領域毎に力学的指標の予測値が算出される。

また、ステップＳ９が行われるとシステム制御部２１は、画像処理装置２７に血液流体解析処理を行わせる（ステップＳ１１）。ステップＳ１１において画像処理装置２７の血液流体解析部５９は、暫定的に構築された力学モデルに血液流体解析を施して時系列の血液流量指標の予測値の空間分布を算出する。具体的には、離散化領域毎に血液流量指標の予測値が算出される。

なお、力学的指標または血液流量指標として、ＦＦＲが算出されても良い。

ステップＳ９及びＳ１０が行われるとシステム制御部２１は、表示機器３１に表示処理を行わせる（ステップＳ１１）。ステップＳ１１において表示機器３１は、ステップＳ９において算出された時系列の力学的指標の予測値とステップＳ１０において算出された時系列の血液流量指標の予測値とを表示する。例えば、表示機器３１は、時系列の力学的指標または時系列の血管流量指標を、時系列の力学モデルを当該予測値に応じた色で動画的に表示する。このため、表示機器３１は、各種の予測値とカラー値（例えば、ＲＧＢ）との関係を示すカラーテーブルを保持している。表示機器３１は、予測値に応じたカラー値をカラーテーブルを利用して特定し、特定されたカラー値に応じた色で当該予測値に対応する離散化領域を表示する。

図１５は、力学的指標の一つである内圧の空間分布の表示例を示す図である。図１５に示すように、表示機器３１は、力学モデルを構成する各離散化領域を当該離散化領域に関する内圧値に応じた色で動画的に表示する。ユーザは、力学モデルを観察することにより、経時的且つ空間的に変化する力学的指標を色で把握することができる。

図１６は、血液流量指標の一つである流速値の空間分布の表示例を示す図である。図１６に示すように、表示機器３１は、力学モデルを構成する各離散化領域を当該離散化領域に関する流速値に応じた色で動画的に表示する。ユーザは、力学モデルを観察することにより、経時的且つ空間的に変化する血液流量指標を色で把握することができる。

例えば、血管内が完全に狭窄している場合、狭窄部位の内圧は、非狭窄部位の内圧よりも小さい。力学的指標として内圧が指定された場合、ユーザは、力学モデル上での局所的な色の違いにより狭窄の有無を判断することができる。また、狭窄部位の流量は、非狭窄部位の流用よりも小さい。血液流量指標として流量が指定された場合、ユーザは、力学モデル上での局所的な色の違いにより狭窄の有無を判断することができる。

また、血管応力解析部５７は、プラーク領域の材料モデルパラメータの同定結果に基づいて力学的指標として硬さ値の空間分布を算出しても良い。この場合、表示機器３１は、プラーク領域に関する硬さ値の空間分布を力学モデル上において表示しても良い。また、表示機器３１は、プラーク領域周辺の内圧分布や応力分布、ひずみ分布を表示して良い。ユーザは、これらの表示をプラークの性状と破綻しやすさとを推定することに活用することができる。

力学的指標及び血液流量指標の予測値は、力学モデルの離散化領域の色で表現する方法のみに限定されない。例えば、図１７、図１８、及び図１９に示すように、グラフで表示しても良い。なお、図１７は、左冠動脈起始部の血圧に関するグラフである。図１７のグラフの縦軸は正規化した血圧に規定され、横軸は心位相［％］に規定される。図１８は、ＬＣＸとＬＤＡとの分岐点付近の血圧に関するグラフである。図１８のグラフの縦軸は正規化した血圧に規定され、横軸は心位相［％］に規定される。図１９は、芯線方向に関する血圧変化に関するグラフである。図１９のグラフの縦軸は血圧比に規定され、横軸は大動脈からの距離［ｍｍ］に規定される。表示機器３１は、力学的指標及び血液流量指標の予測値をグラフで表示することにより、これら値をユーザに簡便に把握させることができる。

ステップＳ１１が行われると構造流体解析処理が終了する。

なお、図１０において、強制変位履歴は、形状モデルの芯線部と外壁部とに設定されるとしたが、強制変位履歴の設定箇所は、これに限定されない。例えば、強制変位履歴は、芯線部と外壁部との間の血管壁領域に設定されても良い。

なお、強制変位履歴の拘束条件の割り当て対象は、境界条件及び材料モデルを同定するか否かに応じて切り分けられても良い。図２０、強制変位履歴の他の割り当て例を示す図であり、形状モデルの断面を示している。例えば、図２０の（ａ）に示すように、境界条件及び材料モデルを同定する場合、形状モデルの外壁部ＯＷ上の節点ＰＮ２のみに強制変位履歴が割り当てられ、血管壁領域ＲＶの節点ＰＮ３に強制変位履歴を割り当てないと良い。図２０の（ｂ）に示すように、境界条件及び材料モデルを同定しない場合、形状モデルの外壁部の節点ＰＮ２と血管壁領域ＲＶの節点ＰＮ３との両方に強制変位履歴が割り当てられると良い。この場合、芯線上の節点ＰＮ１に強制変位履歴が割り当てられる。また、外壁部ＯＷの節点ＰＮ２と節点ＰＮ１とをビーム要素ＥＢで結び、ビーム要素ＥＢ上の節点ＰＮ２及びＰＮ３にも強制変位履歴が割り当てられても良い。このとき、半径方向に関する収縮及び膨張は、ビーム要素ＥＢの伸縮変位で表現する。なお、血管内腔領域ＲＩには、強制変位履歴が割り当てられなくて良い。

図２１は、強制変位履歴の割り当ての他の例を示す図であり、血管周辺組織のダミー要素ＲＤを含む形状モデルの断面を示している。図２１に示すように、ダミー要素ＲＤは、血管壁領域ＲＮの外側に設定される。形状モデルがダミー要素ＲＤを含む場合、血管壁領域ＲＮに加え、ダミー要素ＲＤにも節点ＰＮ４が設定される。節点ＰＮ４にも強制変位履歴が割り当てられる。力学モデル構築部５５は、境界条件及び材料モデルを同定する場合、血管壁領域ＲＶに含まれる節点ＰＮ３に強制変位履歴を割り当て、境界条件及び材料モデルを同定しない場合、強制変位履歴を割り当てなくても良い。節点ＰＮ３に強制変位履歴を割り当てる場合、内腔領域ＲＩに関する形状指標以外にも血管壁領域ＲＶに関する形状指標も参照して材料モデルの同定が行われる。

図２２は、強制変位履歴の割り当ての他の例を示す図であり、プラーク領域ＲＰを含む形状モデルの断面を示している。図２２に示すように、プラーク領域ＲＰは、血管壁領域ＲＶに含まれる。プラーク領域ＲＰは、材料モデル同定領域に設定される。プラーク領域ＲＰについては、内腔形状指標、血管壁形状指標、及びプラーク指標を考慮して材料モデルが同定される。既述のように、プラーク指標は、例えば、超音波診断装置による組織性状診断により得られたプラークの性状に関するデータである。力学モデル構築部５５は、性状に応じてプラーク領域を複数の部分領域に区分し、複数の部分領域に個別に材料モデル同定領域を設定する。各部分領域には、予め、当該部分領域の性状に応じたパラメータ範囲が設定されると良い。既述の統計的同定処理により、各部分領域についての材料モデルパラメータが統計的同定部６１により同定される。そして、ステップＳ１０において表示機器３１が力学的指標として血管の材料特性に関する指標を表示することにより、ユーザは、プラークの性状を正確且つ容易に把握することができる。

次に、潜在変数の一つである材料モデルについて詳細に説明する。血管の材料モデルとしては、弾性モデル、超弾性モデル、異方性超弾性モデル、粘性特性を考慮した超弾性モデルなどが適用可能である。異方性超弾性モデルとしては、例えば、Y.C.Funにより提案された数理モデルや、Holzapfel-Gasser構成式と呼ばれる数理モデルが適用できる。単位参照体積あたりのひずみエネルギーは以下の（１）式で表わされる。（１）式の第１項は、コラーゲンを含まない等方性基礎材料のせん断変形に関するエネルギーを表し、第２項は、コラーゲンを含まない等方性基礎材料の体積変形に関するエネルギーを表し、第３項は、コラーゲン繊維各グループの寄与（繊維方向の分散を考慮）を表している。

図２３は、繊維グループの変形を示す図である。図２３に示すように、円筒形上の外膜を仮定する。芯線方向ｚと周方向θとにより規定される面における平均方向Ａでの繊維グループの変形は、下記の（２）式で表される。

（１）及び（２）式中の材料モデルに関するパラメータとしては、以下の表１に示すように、材料パラメータと繊維分散パラメータとがある。材料パラメータとしては、Ｃ１０、Ｄ，Ｋ１、Ｋ２などが用いられ、繊維分散パラメータとしては、Ｋａｐｐａやγなどが用いられる。各パラメータのデフォルト値と制約条件とは表１に示す通りである。

血液の材料モデルは、以下の（３）式のようなＣａｓｓｏｎ構成式や（４）式のようなＨＢ構成式が好適である。

これら材料モデルのパラメータは、上述のステップＳ６において統計的同定部６１により、血管形態指標及び血液流体指標を用いた統計的同定処理により同定される。

次に、統計的同定部６１により行われる統計的同定処理の詳細について説明する。

時系列の医用画像から計測される血管形態指標及び血液流量指標のような観測変数は、不確定性を有している。統計的同定部６１は、このような不確定性が存在する状況下における潜在変数の統計的同定法として、階層ベイズモデルとマルコフ連鎖モンテカルロ法とに基づく統計的手法を活用する。

既述のように、ステップＳ６において統計的同定部６１は、ステップＳ４において算出された血管形態指標または血液流量指標の予測値と観測値とに基づくデータ分布を設定する。データ分布は、例えば、血管形態指標または血液流量指標の予測値と観測値との誤差に関する多変量正規分布関数を示す。データ分布は、時刻毎に個別に設定されても良いし、複数時刻まとめて設定されても良い。次に、統計的同定部６１は、形状モデルの強制変位と潜在変数とに事前分布を割り当てる。事前分布は、取り得る値の確率分布を示す。次に、統計的同定部６１は、潜在変数に関する数値シミュレーションのパラメータサーベイを実行し、潜在変数と血管形態指標または血液流量指標との関係を表現するモデルを構築する。例えば、材料モデルパラメータと内圧分布パラメータと血管形態指標または血液流量指標との関係がモデルに規定される。なお、血管形態指標または血液流量指標と潜在変数との関係は、モデルという形態ではなく、データベースあるいはテーブルにより規定されても良い。これらモデル、データベース、または、テーブルは、記憶装置３３に記憶される。このモデル、データベース、または、テーブルを利用して統計的同定部６１は、事前分布から血管形態指標または血液流量指標の確率分布を算出する。統計的同定部６１は、階層ベイズモデルとマルコフ連鎖モンテカルロ法とにより得られる事後分布から潜在変数を統計的に同定する。

この同定問題は、具体的には、次の３条件を満足しない不良設定問題となる。３条件は、（１）解の存在が保証される、（２）解が唯一に定まる、（３）解がデータに対して連続的に変化し測定誤差に対して解が安定している、である。不良設定問題は、正規化理論とその拡張とからなる枠組み内で捉えると扱いやすい。不良設定問題を解くことは、標準正規化理論では不十分である。不良設定問題を解くためには、内部状態の不連続点を検出し、検出された不連続点を内部状態の推定に役立てる方法論が必要となる。この点で、本実施形態に係る統計的同定処理においてマルコフ確率場理論が有効となる。

血管形態指標と血液流量指標とに不確定性が存在する環境下において統計的同定部６１は、適切な制約条件の下において潜在変数の確率分布パラメータを同定する。適切な制約条件を決定するためには、解の性質を事前に知っている必要がある。統計的同定部６１は、解空間の制約条件に関するデータベースをシミュレーションと観測値とに基づいて発生する。統計的同定部６１は、発生されたデータベースを利用して、超多自由度大規模問題に対してマルコフ確率場理論と階層ベイズモデルとに基づく統計的同定処理を実行する。制約条件となる事前分布の設定では、多くの数値実験結果に基づいてこれらの要因に関するパラメータの確率分布が並列的に個別に構成される。統計的同定部６１は、複数の確率分布を統合しデータの欠損を補間することで潜在変数のパラメータを同定する。この処理のために、統計的同定部６１は、マルコフ確率場理論を用いたモデルに基づく階層ベイズ法による推定を行う。解析対象とする構造の変形状態の実測結果から、同定した中間変数をもとに、ある負荷条件と境界条件における圧力や流量分布が推定できるという仕組みである。

冠動脈の構造流体解析における材料モデル、境界条件、及び負荷条件の同定問題は、非線形逆解析と位置付けられ、解の一意性及び安定性が保証されない場合が多い。生体組織の材料特性や血圧の現実的に取り得る範囲は先験情報として想定できるため、これらを事前分布の確率分布として設定できる。また、圧力や変位は空間的時間的に滑らかであることも想定できるためこの情報も先験情報として事前分布の確率分布として設定できる。あるいは、血液の流れに逆流が生じていない事実を考慮できる場合は血管芯線方向の全体的な圧力分布の傾きは負（圧力降下が存在）であることも制約条件として用いて良い。負荷条件（内圧分布など）、境界条件、及び材料モデルに対して、時系列のＣＴ画像に基づく血管形状変形指標の観測値と力学モデルに基づく血管形状変形指標の予測値との２乗誤差分布をデータ分布として設定できる。観測可能な平均流量に関する２乗誤差分布もデータ分布として追加してもよい。これらの事前分布とデータ分布とに基づいて階層ベイズモデルとモンテカルロ法とを利用して事後分布を算出することができる。事後分布の発生確率や分散により、潜在変数のパラメータの同定値を得ることができる。発生確率が高く、分散が小さいほうが確信度合の高い同定値といえる。事後分布が多峰性分布となる場合でも、複数の同定値のうちの分散が小さい同定値を選択すれば良い。あるいは、複数の同定値が存在し得る場合、それぞれの同定条件で、構造流体解析を実施し、それぞれの可能性を認識して、同定値や解析結果を、診断や予防の指針情報として活用することができる。時系列のＣＴ画像にも誤差が含まれていることから、力学モデルの各節点に関する血管形態指標にも誤差が含まれる。このため、各血管形態指標を、例えば、時系列のＣＴ画像から計測された血管形態指標の予測値を平均値とした正規分布の確率変数として扱い、位置の空間的順序を保つという制約を含めた上で、事前分布を設定してもよい。また、潜在変数のパラメータの同定において、一意性がなく、複数の候補が考えられる場合がある。この場合、時系列のＣＴ画像から計測された血管形態指標の不確定性に従う乱数のサンプリング点に対する、潜在変数の同定値のサンプル集合の変動幅をチェックすることで、各同定値の候補のロバスト性（安定性）を判定する。各同定値の候補のロバスト性に基づいて最終的な同定値を決定しても良い。

次に、力学モデルの詳細について説明する。既述のように、力学モデル構築部５５は、力学モデルの種類に応じて異なるタイプの力学モデルを構築することができる。連続体力学に基づくＦＥＭを用いる場合、力学モデル構築部５５は、血管壁の応力解析用のための形状モデル（ＦＥＭモデル）と血液の流体解析用のための形状モデル（ＦＥＭモデル）との両方を構築する。

材料力学に基づく簡易的な同定法を用いる場合、材料力学における内圧を受ける厚肉円筒の式から近似的に圧力と弾性率と変位との関係を求める。この場合、芯線方向に配列された複数の離散化領域の各々について厚肉円筒近似が形状モデルとして用いられる。具体的には、力学モデル構築部５５は、芯線上に離散的に配列された節点を通る断面上の血管内腔形状と血管壁表面形状と断面中心とを特定する。次に力学モデル構築部５５は、血管内腔形状と血管壁表面形状とに基づいて平均面積、内腔の平均半径、及び平均壁厚を算出する。そして力学モデル構築部５５は、平均面積、内腔の平均半径、及び平均壁厚に基づいて、各離散化領域の血管領域に厚肉円筒近似を施して形状モデルを構築する。

流れ力学に基づく簡易的な同定法を用いる場合、血流の平均圧力と平均流量とを近似的に求めるため、流体力学における修正ベルヌーイの式、あるいは、Hagen-Poiseuille Flow（ハーゲン・ポアズイユ流れ）の式を用いる。この場合、複数の離散化領域の各々の血圧差と流量との関係を近似的に求めるための形状モデルが構築される。ここで、血圧差は、入口の血圧と出口の圧力との圧力差であり、流量は、単位時間あたり入口流量（あるいは流速）と出口流量（あるいは流速）とを意味する。ただし、力学モデル構築部５５は、ステップＳ２において算出された芯線方向に関する伸縮距離とねじれ角とに基づいて、各節点に関する移動ベクトルと当該節点に隣接する節点に関する回転変位とを各時刻で対応させると良い。

次に、連続体力学モデルを用いる構造流体解析について説明する。

血管及び血液（血管と血液を構成する物質）についての運動学は、その運動を引き起こしている力とは無関係である。血管及び血液についての運動学の基本概念は、位置、時間、物体、運動、および変形し得る物質の集合についての直感的概念を数学的用語へ抽象化したものである。血管及び血液に関する変形及び運動の局所的解析を支配する基礎的運動学テンソルは、変形勾配テンソルFと速度勾配テンソルLとである。変形勾配テンソルFは、運動する血管及び血液の物質要素に生じる大きさ及び形の変化を規定する。変形勾配テンソルFは、回転テンソル（正格直交テンソル）Rとストレッチングテンソル（正値対称テンソル）U,Vとの積で表わされる。ストレッチングテンソルU,Vは、まず、基本形態での正規直交ベクトルRによって定められる方向に関するストレッチを課し、次に、正規直交ベクトルRによって与えられる剛体回転を課すことによりもたらされる。なお、ストレッチを課す順番と剛体回転を課す順番とは逆でも良い。速度勾配テンソルLは、基準形態に依存せず現在形態にのみ依存する。速度勾配テンソルLは、運動する血管及び血液の物質要素に生じる大きさ及び形の変化が生じる速度を規定する。速度勾配テンソルLは、ひずみ速度テンソルD（対称テンソル）とスピンテンソル（反対称テンソル）とに分離できる。ひずみ速度テンソルDは、物体がその現在形態をちょうど通過するときのストレッチの変化率を表す。スピンテンソルは、物体がその現在形態をちょうど通過するときの回転の変化率を表す。血管の力学モデルの外表面節点（あるいは積分点）の一部に関する変形勾配テンソルと速度勾配テンソルとに変位拘束条件（時間的変化を含む）を割り当てる。力学モデルの内腔節点（あるいは積分点）に関する予測値（変形勾配テンソル、速度勾配テンソル、あるいはそれらの関数値（例えば、変位や面積でも良い））と、観測値（観測データから得られた変形勾配テンソル、速度勾配テンソル、あるいはそれらの関数値）とが整合するように、材料モデルのパラメータ、負荷条件（力学モデルの内腔における表面力ベクトル）、及び境界条件（血管境界における力ベクトル）を同定する。ここで、血管内部の応力の初期状態は、予め仮定しておいてもよいし、同定してもよい。

連続体力学に基づく力学モデルは、運動中の血管および血液における質量、運動学、運動量、角運動量、及びエネルギーの平衡を表す方程式を基礎としている。質量、力、熱、及び内部エネルギーという概念が基本である。平衡則とは、運動量の時間全微分が物体力と接触力との和に等しくなり、また、角運動量の時間全微分が物体トルクと接触トルクとの和に等しくなり、また、運動エネルギーと内部エネルギーとの時間的変化が、仕事率（力学的エネルギー）と単位時間あたりの熱供給と熱流束との和に等しくなる、ということを意味している。平衡則、構成式、及び跳躍条件から、変形勾配テンソル、速度勾配テンソル、及び応力テンソルといった場の方程式を導き、血管及び血液の力学モデルを記述できる。また、血管及び血液におけるひずみ場は、適合条件を満たす。ここで、構成式とは、密度、内部エネルギー、速度ベクトル、応力テンソル、熱流束ベクトル、及び温度からなる１０個のスカラ方程式の組の関係を与えるものである。場の釣り合い式の中の１７個のスカラ場、つまり、密度、内部エネルギー、速度ベクトル、応力テンソル、及び熱流束ベクトルのうち、場の釣り合いにより８個のスカラ関係が与えられ、温度を加えた残りの未知量の関係を与える。ただし、物体力ｂと熱源ｒとは既知としている。これらのスカラ場を与える方程式のパラメータが材料モデルパラメータである。

連続体力学に基づく力学モデルは、有限要素法又は境界要素法による数値解析法により変位ベクトル、応力テンソル、ひずみテンソル、及び速度ベクトル等の場の近似解を、与えられた境界条件、負荷条件、及び材料モデルのもとに算出することができる。

構造−流体連成解析において、構造及び流体の方程式を解く方法は、一体型解法（monolithic method）と分離型解法（partitioned method）とのいずれでも良い。また、構造と流体との間の境界面での連成は、弱連成でも強連成でもよい。また、流体解析においては、ＡＬＥ法に代表されるような境界面追跡型の手法で移動境界を扱ってもよいし、Immersed Boundary method、Immersed Finite Element Method、Fictitious Domain Methodなどの境界面補足型の手法でもよい。

次に、簡易的な力学モデルの例として、内圧と外圧とを受ける厚肉円筒の材料力学の式と、ハーゲン・ポアズイユ流れ及び修正ベルヌーイの式とについて詳細に説明する。

まず、図２４と図２５とを参照しながら、厚肉円筒の材料力学の式について説明する。図２４は、肉厚円筒の力学モデルの直交断面を示す図であり、図２５は、図２４の微小扇形要素の拡大図である。内半径r_a、外半径r_bの厚肉円筒に内圧p_aと外圧p_bとが作用する場合の応力やひずみ、変位などの式を説明する。E及びνは材料モデルパラメータである。Eは弾性率、νはポアソン比を表している。厚肉円筒では，半径応力σ_rも考慮し、円周応力σ_θの半径方向分布も考慮する必要がある。以降では、軸方向のひずみε_zは、断面の位置及び向きに関して一様であるとする。また円筒断面は、軸対称である場合について説明するが任意形状であっても良い。円筒断面が軸対称である場合、平衡条件は、任意の断面の半径方向についてのみ考えれば良い。任意の断面上で半径rとr+drの同心円筒と中心角dθで切り取られた単位厚さ1の微小扇形要素について、半径方向に関する力の平衡を考える。変形も軸対称であるため、ab面及びbc面には、せん断応力が生じないので、垂直応力のみ作用する。従って、半径方向に関する力の平衡は、以下の（５）式のように表現することができる。

σ_r rdθ+2σθ_r sin(dθ/2)-（σ_r+(dσ_r/dr)dr）(r+dr)dθ=0 …（５）
ここで、drはrより小さく、dσ_rはσ_rより小さいので、（５）式に含まれる高次の微小項を省略し、sin(dθ/2)≒dθ/2とすれば、（５）式は、以下の（６）式のように表現することができる。

rdσ_r/dr+σ_r-σ_θ=0 …（６）
半径rにおける半径方向の変位をuとすれば、u+drでの同方向の変位はu+(du/dr)drとなるので、半径rにおける半径方向のひずみε_rはε_r=du/drとなる。また、半径方向の変位uによって、半径rの円は半径r+uの円になる。従って、円周ひずみε_θは、以下の（７）式のように表現することができる。

ε_θ=(2π(r+u)-2πr)/2πr=u/r …（７）
また、応力とひずみの関係式から、以下の（８）式または（９）式が得られる。

d²u/dr²+(1/r)(du/dr)-u/r²=0 …（８）
d²u/dr²+d(u/r)/dr=0 …（９）
（８）式または（９）式を積分すると以下の（１０）式を得ることができる。

u=c₁r+c₂/r …（１０）
これにより以下の（１１）式、（１２）式、（１３）式が得られる。

σ_r=(E/((1+ν)(1-2ν)))(c₁-(1-2ν)(c₂/r²)+νε_z) …（１１）
σ_θ=(E/((1+ν)(1-2ν)))(c₁+(1-2ν)(c₂/r²)+νε_z) …（１２）
σ_z=(Eν/((1+ν)(1-2ν)))(2c₁+((1-ν)/ν)ε_z) …（１３）
（１１）式、（１２）式、（１３）式における定数c₁，c₂は、周辺条件、すなわち、円筒の内周r=raでσ_r=-p_a、外周r=r_bでσ_r=-p_bから定めることができる。この周辺条件により、（１１）式、（１２）式、及び（１３）式から以下の（１４）式、（１５）式、（１６）式をそれぞれ得ることができる。また、変位uは、以下の（１７）式のように表現することができる。

σ_r=(1/(r_b ²-r_a ²))(r_a ²(1-r_b ²/r²)p_a-r_b ²(1-r_a ²/r²)p_b) …（１４）
σ_θ=(1/(r_b ²-r_a ²))(r_a ²(1+r_b ²/r²)p_a-r_b ²(1+r_a ²/r²)p_b) …（１５）
σ_z=2ν(r_a ²p_a-r_b ²p_b)/(r_b ²-r_a ²)+Eε_z=ν(σ_r+σ_θ)+Eε_z …（１６）
u=((1+ν)(1-2ν)/E)((r_a ²p_a-r_b ²p_b)/(r_b ²-r_a ²))r
+((1+ν)/E)((r_a ²r_b ²)/((r_b ²-r_a ²)r))(p_a-p_b)-νε_zr …（１７）
上記の（１６）式及び（１７）式は、ε_zの項が含んでいる。従って、（１６）式は、対象とする円筒の境界における拘束条件に応じて異なる。例えば、対象とする円筒の両端が拘束される場合、ε_z=0となり、両端が開放され、σ_z=0となる。

rθ面内のせん断応力τ_rは、以下の（１８）式のように表現することができる。

τ_r=(1/2)(σ_r-σ_θ)=((r_a ²r_b ²)/((r_b ²-r_a ²)r²))p_b …（１８）
θz面内のせん断応力τ_r’は、σ_z=0のとき最大になるのでσ_z=0とすると、以下の（１９）式のように表現することができる。

τ_r’=(1/2)|σ_θ|=(1/2)((r_b ²)/(r_b ²-r_a ²))(1+r_a ²/r²)p_b …（１９）
次に、ハーゲン・ポアズイユ流れ及び修正ベルヌーイの式について説明する。

流体が円管内へ流入すると、下流に進むにつれて圧力は降下するとともに、流れの速度分布も徐々に変化する。流れが管内へ流入すると管壁から境界層が発達し、下流へ進むにつれて境界層は厚さを増し、ついには管内の流れは境界層に覆われる。このため、速度分布は管入口のほぼ平らな分布から下流の放物形分布へと変化し、それ以降の速度分布は変化しなくなる。この状態を完全に発達した流れと呼び、管摩擦損失による圧力降下も一定の割合になる。流れが管入口から発達した流れに達する区間を助走区間または入口区間といい、その区間の長さを助走距離または入口長さと呼ぶ。完全に発達した流れの速度分布は下流方向へ変化しないから、以下の（２０）式に示すように、管摩擦損失によって生じる圧力損失ΔPの作用力と流体の粘性によって生じるせん断応力τの摩擦力は釣り合うことになる。円管内流れの場合、レイノルズ数Reがおよそ2300以下のとき層流となるといわれている。

速度分布uは層流の場合、以下の（２１）式のように軸対称な回転放物面で表わされる。

流量Qは、速度分布uを管断面全体にわたって積分することにより、以下の（２２）式のように表現することができる。

断面平均流速vは、以下の（２３）式のように表現することができる。

圧力勾配は、流れ方向へ一定で、圧力は減少する。圧力勾配は、管長lの圧力降下Δpを用いると、以下の（２４）式のように表現できる。（２４）式に（２２）式を代入すると（２５）式が得られる。（２５）式は、流量Qが圧力損失Δpに比例することを意味する。この関係を満たす流れをハーゲン・ポアズイユ流れという。

損失がある場合の修正ベルヌーイの式は、以下の（２５）式のように表現することができる。

ここで、lは管の長さ、dは管内径、vは管内平均速度、λは管摩擦係数である。管摩擦係数λは、流れが層流の場合にはレイノルズ数Reによって、乱流の場合にはレイノルズ数Reと表面粗さとによって定まる値になる。流体の粘性によって摩擦抵抗が必ず作用する。この摩擦抵抗は、流れを駆動する動力あるいはエネルギーを消費することになるのでエネルギー損失になる。

ハーゲンポワズイユ流れの式をダルシー・ワイズバッハの式のように変形すると、圧力差と流速および管内径とに関する以下の（２７）式が得られる。

以上でハーゲン・ポアズイユ流れ及び修正ベルヌーイの式の説明を終了する。

上述のように、画像処理装置２７は、この材料力学の式とハーゲン・ポアズイユ流れ及び修正ベルヌーイの式とを利用して潜在変数を同定することができる。例えば、力学モデルとして、血管の変形を厚肉管の材料力学の式を用い、管径変化を内圧変化と弾性率とにより表現する場合について考える。無応力状態を初期形状（例えば、血管が最も収縮する状態）と仮定した場合、血管壁及びプラークの弾性率をある値に設定すると、血管内腔の平均半径などの血管形状変形指標の観測値の時間的変化量と内圧の変化量との関係式が得られる。血管形状変形指標の観測値は、時系列のＣＴ画像から計測される。この血管形状変形指標の観測値の時間的変化量に合致するように血管の内圧分布の時間的変化が決定される。この内圧分布の下に血液の流体解析を行うことで血管流量指標の予測値が計測される。この血管流量指標の予測値が観測値に一致しない場合、画像処理装置２７は、最初に決めた血管壁あるいはプラークの弾性率を変更して、さらに同様の解析を行う。これを繰り返すことにより、画像処理装置２７は、血管形状変形指標の観測値と血液流量指標の観測値とに整合する血管壁およびプラークの弾性率、内圧分布、流体解析の圧力境界条件などの潜在変数を決定することができる。この決定方法をより効率的かつ安定的に行うために、階層ベイズモデルとマルコフ連鎖モンテカルロ法とによる統計的同定手法を用いてもよい。

上記のように、本実施形態に係る血管解析装置５０は、記憶装置３３、領域設定部５１、画像解析・追尾処理部５３、力学モデル構築部５５、及び統計的同定部６１を有している。記憶装置３３は、被検体の血管に関する時系列の医用画像のデータを記憶する。領域設定部５１は、時系列の医用画像に含まれる血管領域に解析対象領域を設定する。画像解析・追尾処理部５３は、時系列の医用画像を画像処理して解析対象領域の時系列の形態指標と時系列の形状変形指標とを算出する。力学モデル構築部５５は、時系列の形態指標と時系列の形状変形指標と時系列の医用画像とに基づいて解析対象領域の構造流体解析に関する力学モデルを暫定的に構築する。統計的同定部６１は、暫定的に構築された力学モデルに基づく血管形態指標及び血液流量指標の予測値が予め計測された血管形態指標及び血液流量指標の観測値に整合するように解析対象領域に関する力学モデルの潜在変数を同定する。

上記構成により、本実施形態に係る血管解析装置５０は、材料モデル、境界条件、負荷条件、及び幾何学的構造等の潜在変数を血管形状変形指標と血液流量指標とを用いた逆解析により同定することができる。血管解析装置５０は、潜在変数を変更しながら逆解析を反復して行うことにより、上述の４点の困難、すなわち、１．冠動脈の材料モデルの同定方法、２．心臓の形状の変形の冠動脈への影響の組み込み、３．冠動脈の境界条件の同定方法、４．不確定性を有する血管形状を利用した材料モデルや負荷条件、境界条件の同定方法、を全て加味した潜在変数を同定することができる。従って、血管解析装置５０は、ＣＴ画像に描出されない血管や心臓等の外部要因による影響を加味した構造流体解析を実行することができる。

かくして、本実施形態によれば、血管の構造流体解析の精度の向上が可能となる。

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

１０…ＣＴ架台、１１…Ｘ線管、１３…Ｘ線検出器、１５…データ収集装置、２０…コンソール、２１…システム制御部、２３…架台制御部、２５…再構成装置、２７…画像処理装置、２９…入力機器、３１…表示機器、３３…記憶装置、５０…血管解析装置、５１…領域設定部、５３…画像解析・追尾処理部、５５…力学モデル構築部、５７…血管応力解析部、５９…血液流体解析部、６１…統計的同定部

Claims (16)

1. 被検体の血管に関する時系列の医用画像のデータを記憶する記憶部と、
   前記時系列の医用画像に含まれる血管領域に解析対象領域を設定し、前記解析対象領域に潜在変数同定領域を設定する設定部と、
   前記時系列の医用画像を画像処理して前記解析対象領域の時系列の形態指標と時系列の形状変形指標とを算出する算出部と、
   前記時系列の形態指標と前記時系列の形状変形指標と前記時系列の医用画像とに基づいて前記解析対象領域の構造流体解析に関する力学モデルを暫定的に構築する構築部と、
   前記暫定的に構築された力学モデルに基づく血管形態指標の予測値及び血液流量指標の予測値の少なくとも一方が予め計測された血管形態指標の観測値及び血液流量指標の観測値の少なくとも一方に整合するように前記潜在変数同定領域に関する潜在変数を同定する同定部と、
   を具備する血管解析装置。
2. 前記同定部は、前記潜在変数として、前記潜在変数同定領域に関する材料モデル、前記潜在変数同定領域の血液流入口及び血液流出口に関する境界条件、前記潜在変数同定領域の血流による内圧分布を含む負荷条件、及び前記潜在変数同定領域の幾何学的構造の少なくとも一つを同定する、請求項１記載の血管解析装置。
3. 前記同定部は、
   前記暫定的に構築された力学モデルに応力解析を施して前記血管形態指標の予測値を推定する血管応力解析部と、
   前記暫定的に構築された力学モデルに流体解析を施して前記血液流量指標の予測値を推定する血液流体解析部と、
   前記血管形態指標の予測値が前記血管形態指標の観測値に整合するように前記潜在変数を統計的に同定する第１統計的同定部と、
   前記血液流量指標の予測値が前記液流量指標の観測値に整合するように前記潜在変数を統計的に同定する第２統計的同定部と、
   を有する、請求項２記載の血管解析装置。
4. 前記構築部と前記同定部とを制御する制御部をさらに備え、
   前記制御部は、前記同定部により潜在変数が同定される毎に所定の同定終了指標が既定値に到達しているか否かを判定し、前記同定終了指標が前記既定値に到達していないと判定した場合、前記構築部を制御して前記同定された潜在変数を利用して力学モデルを再構築させ、前記同定部を制御して前記再構築された力学モデルの潜在変数を同定させ、前記同定終了指標が前記既定値に到達した判定した場合、前記構築部を制御して最新の力学モデルを最終的な力学モデルに設定させる、
   請求項１記載の血管解析装置。
5. 前記最終的な力学モデルに応力解析を施して時系列の力学的指標を算出する血管応力解析部と、
   前記時系列の力学的指標を表示する表示部と、
   をさらに備える請求項４記載の血管解析装置。
6. 前記最終的な力学モデルに流体解析を施して時系列の血液流体指標を算出する血液流体解析部と、
   前記時系列の血液流体指標を表示する表示部と、
   をさらに備える請求項４記載の血管解析装置。
7. 前記最終的な力学モデルに基づく力学的指標と血液流量指標との少なくとも一方の予測値が力学的指標と血液流量指標との少なくとも一方の観測値に整合するように前記時系列の医用画像に含まれる血管領域を修正する修正部、をさらに備える請求項４記載の血管解析装置。
8. 前記構築部は、前記時系列の形態指標と前記時系列の医用画像とに基づいて前記解析対象領域の形状モデルを構築し、前記形状モデルに前記時系列の形状変形指標を強制変位として割り当て、前記潜在変数を割り当てることにより前記力学モデルを構築する、請求項２記載の血管解析装置。
9. 前記構築部は、前記力学モデルを材料モデル同定領域と非同定領域とに区分し、前記材料モデル同定領域の表面部に強制変位の拘束条件を割り当て、前記材料モデル同定領域の内部に強制変位の拘束条件を割り当てず、前記非同定領域の表面部及び内部に強制変位の拘束条件を割り当てる、請求項８記載の血管解析装置。
10. 前記力学モデルに基づく、血圧による内圧の時系列の空間分布、材料モデルパラメータ、及び血管形態指標を関連付けて記憶する関連付け記憶部をさらに備える、請求項１記載の血管解析装置。
11. 前記時系列の形態指標と前記時系列の形状変形指標とは確率分布で与えられる、請求項１記載の血管解析装置。
12. 前記構築部は、前記暫定的に構築された力学モデルを超音波診断装置からのプラーク指標に応じて複数の部分領域に区分し、前記複数の部分領域の各々に材料モデルのためのパラメータを割り当てる、請求項１記載の血管解析装置。
13. 前記血液流量指標の観測値は、前記時系列の医用画像に基づく血管内腔の容積の時間変化、あるいは、超音波ドプラ計測に基づく血流量の時間変化である、請求項１記載の血管解析装置。
14. 被検体の血管に関する時系列の医用画像のデータを発生する撮像機構と、
    前記時系列の医用画像に含まれる血管領域に解析対象領域を設定し、前記解析対象領域に潜在変数同定領域を設定する設定部と、
    前記時系列の医用画像を画像処理して前記解析対象領域の時系列の形態指標と時系列の形状変形指標とを算出する算出部と、
    前記時系列の形態指標と前記時系列の形状変形指標と前記時系列の医用画像とに基づいて前記解析対象領域の構造流体解析に関する力学モデルを暫定的に構築する構築部と、
    前記暫定的に構築された力学モデルに基づく血管形態指標の予測値及び血液流量指標の予測値の少なくとも一方が予め計測された血管形態指標の観測値及び血液流量指標の観測値の少なくとも一方に整合するように前記潜在変数同定領域に関する潜在変数を同定する同定部と、
    を具備する医用画像診断装置。
15. 時系列の医用画像に含まれる血管領域に解析対象領域を設定し、前記解析対象領域に潜在変数同定領域を設定し、
    前記時系列の医用画像を画像処理して前記解析対象領域の時系列の形態指標と時系列の形状変形指標とを算出し、
    前記時系列の形態指標と前記時系列の形状変形指標と前記時系列の医用画像とに基づいて前記解析対象領域の構造流体解析に関する力学モデルを暫定的に構築し、
    前記暫定的に構築された力学モデルに基づく血管形態指標の予測値及び血液流量指標の予測値の少なくとも一方が予め計測された血管形態指標の観測値及び血液流量指標の観測値の少なくとも一方に整合するように前記潜在変数同定領域に関する潜在変数を同定する、
    ことを具備する血管解析方法。
16. コンピュータに、
    時系列の医用画像に含まれる血管領域に解析対象領域を設定し、前記解析対象領域に潜在変数同定領域を設定する機能と、
    前記時系列の医用画像を画像処理して前記解析対象領域の時系列の形態指標と時系列の形状変形指標とを算出する機能と、
    前記時系列の形態指標と前記時系列の形状変形指標と前記時系列の医用画像とに基づいて前記解析対象領域の構造流体解析に関する力学モデルを暫定的に構築する機能と、
    前記暫定的に構築された力学モデルに基づく血管形態指標の予測値及び血液流量指標の予測値の少なくとも一方が予め計測された血管形態指標の観測値及び血液流量指標の観測値の少なくとも一方に整合するように前記潜在変数同定領域に関する潜在変数を同定する機能と、
    を実現させる血管解析プログラム。

Images