JP2011235009A - 超音波診断装置及び該装置における数値シミュレーション方法 - Google Patents

Abstract

【課題】超音波診断装置により得られる超音波ビーム方向の速度情報から、３次元的な流れの成分を考慮に入れて、流速ベクトルや圧力を求める。
【解決手段】被検体内の各計算点の超音波ビーム方向についての速度情報を初期値として、ナビエ・ストークス方程式に基づく数値計算式を時間ステップごとに反復演算することにより、各計算点について、超音波ビーム方向の速度成分と超音波ビーム方向に直交する方向の速度成分とを含む速度ベクトル及び圧力、を計算する。この反復演算において、流れの範囲の境界の近傍の計算点については、数値演算式における圧力に対応する項の入力値に対し、補正曲線Ｃに従って、基準値からその入力値までの差が閾値以上の場合にはその差を緩和する補正を行った上で数値計算式を実行する。
【選択図】図４

Description

本発明は、超音波診断装置に関し、特に超音波診断装置により得られる超音波ビーム方向の速度成分の情報から、２次元的な流速ベクトルや圧力を求めるための技術に関する。

超音波診断装置では、例えば心臓内の血流における各点からドプラ情報（速度情報）が取得され、それらを２次元マッピングすることにより２次元血流画像が構成される。これはカラードプラ法あるいはカラーフローマッピング法といわれている。一般の超音波診断装置において表示される速度は、血流の真の速度ではなく、ビーム方向に沿った速度成分である。

このような状況に対し、ビーム方向の速度成分からビームに直交する方向の速度成分を推定する方法が提案されている（特許文献１参照）。

また、特許文献２には、観測面内の各点のビーム方向速度成分から、各点のビームに直交する速度成分を推定し、この推定結果に基づき各点の圧力を推定して表示する装置が開示されている。

また、特許文献３に開示される装置は、超音波信号を体内の血管に向けて放射し、反射した超音波信号を受信する超音波計測部と、受信した信号により、血管形状と血管内の血流の流速とを得る解析処理部と、解析処理部からの血管形状から、計算格子を設定して血流速と圧力分布のシミュレーションを行うシミュレーション部と、解析処理部からの血流速と、シミュレーション部からの血流速との誤差を計算して、シミュレーション部にフィードバックするフィードバック部と、フィードバック後のシミュレーション部からの血流速と圧力分布の出力を表示する表示部とを備える。シミュレーション部は、ナビエ・ストークス方程式に基づく数値シミュレーションを行うことで、２次元流速と圧力を計算している。

特開２００５−１１０９３９号公報 特開２００２−１７７２６号公報 特開２００４−１２１７３５号公報

特許文献３の手法では、流れを超音波ビームの走査面内の２次元流れとして取り扱っており、実際の血流の３次元的な流れの成分を考慮に入れていない。

本発明は、超音波診断装置により得られる超音波ビーム方向の速度情報から、３次元的な流れの成分を考慮に入れて、流速ベクトルや圧力を求めるための手法を提供する。

本発明に係る超音波診断装置は、超音波ビームを被検体内に送受信し、被検体内からのエコー信号を取得する送受信手段と、前記エコー信号に基づき、前記被検体内の各点の超音波ビーム方向についての速度情報を求めるドプラ処理手段と、前記エコー信号に基づき、前記被検体内の流体が流れる流れ範囲の境界を検出する境界検出手段と、前記境界検出手段が検出した前記境界についての境界条件の下で、前記被検体内の各計算点の超音波ビーム方向についての速度情報を初期値として、ナビエ・ストークス方程式に基づく数値計算式を時間ステップごとに反復演算することにより、前記各計算点について、前記超音波ビーム方向の速度成分と前記超音波ビーム方向に直交する方向の速度成分とを含む速度ベクトル及び圧力、を計算する反復演算手段と、前記反復演算手段が計算した前記各計算点についての前記速度ベクトル又は前記圧力を表示するための処理を行う表示処理手段と、を備え、前記反復演算手段は、前記境界の近傍の計算点については、前記数値演算式における圧力に対応する項の入力値に対し、当該項の基準値からその入力値までの差が閾値以上の場合にはその差を緩和する補正を行った上で前記数値計算式を反復演算する、ことを特徴とする。

１つの態様では、前記境界検出手段が、前記エコー信号から求められるパワーフロー画像に基づき、前記境界を検出する、ことを特徴とする。

別の態様では、前記表示処理手段は、前記速度ベクトルを表示する速度表示画像における前記各計算点の表示色の色相を当該計算点における前記速度ベクトルの方向に基づき求め、前記速度表示画像における前記各計算点の表示色の色相以外のあらかじめ定められた色パラメータを当該計算点における前記速度ベクトルの大きさに基づき求める、ことを特徴とする。

更に別の態様では、前記反復演算手段により求められた前記計算点における前記速度ベクトルから当該計算点の渦度を計算する渦度計算手段、を更に備え、前記表示処理手段は、前記各計算点の渦度の大きさを色で表示した渦度表示画像を生成する、ことを特徴とする。

本発明に係る数値シミュレーション方法は、被検体内からのエコー信号に基づき、前記被検体内の各点の超音波ビーム方向についての速度情報を求めるステップと、前記エコー信号に基づき、前記被検体内の流体が流れる流れ範囲の境界を検出する境界検出ステップと、前記境界検出ステップで検出された前記境界についての境界条件の下で、前記被検体内の各計算点の超音波ビーム方向についての速度情報を初期値として、ナビエ・ストークス方程式に基づく数値計算式を時間ステップごとに反復演算することにより、前記各計算点について、前記超音波ビーム方向の速度成分と前記超音波ビーム方向に直交する方向の速度成分とを含む速度ベクトル及び圧力、を計算する反復演算ステップと、前記反復演算ステップで計算された前記各計算点についての前記速度ベクトル又は前記圧力を表示するための処理を行う表示処理ステップと、を備え、前記反復演算ステップでは、前記境界の近傍の計算点については、前記数値演算式における圧力に対応する項の入力値に対し、当該項の基準値からその入力値までの差が閾値以上の場合にはその差を緩和する補正を行った上で前記数値計算式を反復演算する、ことを特徴とする。

本発明によれば、超音波診断装置により得られる超音波ビーム方向の速度情報から、３次元的な流れの影響を考慮に入れて、流速ベクトルや圧力を求めることができる。

流体モデルを説明するための図である。 実施形態の数値シミュレーションの概略的な流れを示すフローチャートである。 境界の法線ベクトルの求め方の一例を説明するための図である。 ２．５Ｄ補正で用いる補正曲線の一例を示す図である。 流速ベクトルｕの方向と大きさに基づき表示色を決める色変換マップを説明するための図である。 注目点の圧力の時間変化を示す表示の例を模式的に示す図である。 ２つの注目点の圧力の時間変化を示す表示の例を模式的に示す図である。 実施形態の超音波診断装置の構成の一例を示す機能ブロック図である。

本実施形態に掛かる装置は、超音波診断装置により得られる、２次元の走査面内の各点の超音波ビーム方向についての速度成分から、数値シミュレーションにより２次元の流速ベクトルと圧力を計算するためのものである。本実施形態では、ナビエ・ストークス方程式とヘルムホルツ・ホッジ分解定理とに基づいて構築した数値計算処理を、Ｂモード画像やパワーフロー画像から求めた境界線に基づく境界条件の下で反復実行する。この数値計算において、境界近傍の圧力に対して、３次元的な影響を考慮した補正を加えることで、反復実行の間に３次元的な流れの影響を反映した流速ベクトルや圧力を求める。以下、本実施形態において超音波診断装置が行う処理の全体像、数値シミュレーションの概要及び詳細、具体的な装置構成の例を順に説明する。

＜処理の概要＞
この実施形態の処理手順は、概略以下の通りである。
１．超音波探触子から被検体内へ超音波ビームを送受波することにより、１フレーム分の超音波ビームデータ（すなわちエコー信号データ）を取得する。取得された超音波ビームデータを、この実施形態の数値シミュレーション部に読み込む。読み込んだ超音波ビームデータには、例えば、超音波ビーム走査面内の各点のエコーの強さを示すＢモードデータと、それら各点の流体速度を表すカラードプライメージング（ＣＤＩ）データとが含まれる。
２．必要に応じて、読み込んだ超音波ビームデータのスキャンコンバージョン（走査変換：例えば、セクタ走査の各ビームライン上の各点のデータを、表示に適した直交座標系に変換するなど）を行う。
３．ＣＤＩの流体速度データを、シミュレーションのための流体モデルに入力する。ここで、流体（血流など）の流れる範囲の境界の情報を、例えばＢモードデータ又はパワーフローのデータから求め、流体モデルに入力する。
４．流体モデルを用いた数値流体シミュレーションを、Ｎ回繰り返す（Ｎは、あらかじめ定められた繰り返し回数）。
５．流体シミュレーションにより、超音波ビームに直交する方向の流体の速度成分が計算される。
６．結果を表示する。
７．超音波ビームデータの次フレームに移り、ステップ１〜６を繰り返す。

＜流体モデル＞
実施形態の数値シミュレーションにおいて用いる流体モデルについて説明する。この流体モデルは、実際の流体を数学的に表したものである。この例では、流体は、図１に示すように、直交格子の格子点（ノード１０）により表現される。すなわち、この例では、一般的な超音波探触子により得られる平面内の超音波ビームデータに基づき、その平面内の流体の挙動をシミュレーションする。ノード１０の各々は、物理量として、密度、速さ及び圧力を有する。各ノード１０の位置座標は固定であるが、各ノード１０の持つ物理量は時間とともに変化し得る。直交格子のｘ方向及びｙ方向の格子間隔を、それぞれｄｘ、ｄｙとする。このように、ノード１０は、本実施形態においてシミュレーションの計算の対象となる計算点である。

流体シミュレーションの演算では、連続的な時間モデルで個々のノード１０の持つ密度、速さ、圧力を計算する。これらの計算は、時間ステップごとにナビエ・ストークス方程式を解くことにより行われる。

＜数値シミュレーションの概要＞
流体シミュレーションの対象とする流体が非圧縮性流体かつニュートン流体であると仮定する。非圧縮性流体とは、密度変化が無視できる流体のことであり、自然界の多くの流体がこれにあたる。ニュートン流体とは、流れのせん断応力（接線応力）と流れの速度勾配（せん断速度）が比例した粘性の性質を持つ流体のことである。非圧縮性は、実際は流体の「流れ」についての条件（流体自体についての条件ではない）であり、その密度が時間を通して一定のままであるという前提で流体が流れることを意味する。ニュートン流体は、一定の粘性をもつ流体であることを意味する。これらの仮定は両方とも、実施形態の主たる適用対象として想定される血流に対して合理的に妥当する。

このような非圧縮性且つニュートン流体の挙動を表す基礎方程式は、以下に示す周知のナビエ・ストークス微分方程式である。

ここで、u(x,y)は点(x,y)の持つ速度ベクトル、p(x,y)は点(x,y)での圧力、f(x,y)は点(x,y)における単位体積に作用する外力を表すベクトルである。ρは流体の密度であり、この実施形態では一定であると仮定する。また、ニュートン流体の場合粘性νは一定である。また、非圧縮性流体の場合、次式で表す質量保存則が成り立つ。

上記ナビエ・ストークス方程式（１）の左辺は速度の時間変化を示す。また、右辺の第１項は移流項、第２項は粘性拡散を表す粘性項、第３項は圧力勾配を表す圧力項、第４項は外力を表す外力項である。

なお、血流に働く外力を知ることはできないが、以下では外力は、定量化された速さ成分の中にすでに内在すると考え、外力項をなくして方程式を解く。この場合、ナビエ・ストークス方程式は以下の式（３）となる。

この方程式を解く上での主要な問題は、圧力項が知られていないということである。この問題を扱うためにヘルムホルツ・ホッジ（Helmholtz-Hodge）分解定理を用いる。ヘルムホルツ・ホッジ分解定理とは、任意のベクトル場ｗが、発散のない（すなわち質量保存の）ベクトル場ｕと、スカラー場の勾配∇ｑとに分解できるという定理である。すなわち、次式が成り立つ。

この定理は、境界条件が満たされる場合にのみ有効である。ここでの境界条件は、境界近傍ではｕが境界に平行であり、ｑが境界を横切るように変化しないということである。式で表現すると以下のようになる。

ここで、ｎは境界の法線ベクトルである。

式（４）のベクトル場ｗとスカラー場ｑを、流体の物理量と対応づけるために、式（４）の両辺を時間で微分すると、次のようになる。

従って、次の関係が成り立つ。

この式をナビエ・ストークス方程式（３）と比較することで、次の２つの関係式（６）及び（７）が成り立つことが分かる。

式（６）からはベクトル場ｗが中間的な速度場に該当し、スカラー場ｑが圧力の時間微分に該当することが分かる。この認識のもとで、式（４）に戻る。式（４）の両辺の勾配を求めると、次のようになる。

ここで、質量保存則すなわち式（２）を考慮すると、次の関係式（８）が成り立つことが分かる。

以上のことから、次のような流体シミュレーションアルゴリズムを構築することができる。

すなわち、このアルゴリズムでは、まず式（６）に基づき中間的な速度ベクトルｗを計算し、次に式（８）に基づきｑを計算し、次に式（４）に基づき質量保存の成り立つ流体の速度ベクトルｕを計算し、更に式（７）に基づき圧力ｐを計算する。

ここで、式（４）（ヘルムホルツ・ホッジ分解）は境界条件（式（５））が満たされる場合しか成り立たないので、以上のアルゴリズムの各計算式もその境界条件を満たす必要がある。境界条件を考慮に入れると、上述のアルゴリズムは次のようになる。

このアルゴリズムは、図２に示すような手順となる。すなわち、まず、超音波ビームの送受信により得られた各ノードについてのビーム方向の速度成分のみの速度情報ｕが初期値として入力される。ここで、超音波ビームの走査範囲全体の速度情報を入力してもよいが、実際のシミュレーション演算の対象となるのは、その走査範囲のうち後述する境界線により画定される流れ（血流など）の範囲内だけなので、その流れの範囲に属するノードの速度情報のみを入力するようにしてもよい。後者の場合、シミュレーションのための数値演算は、入力したノードについてのみ行う（すなわち、流れの範囲の外のノードについては計算しなくてよい）。

図２の手順では、この入力された流体の速度ベクトルｕを、式（５）の境界条件を満たすように補正する（Ｓ１）。境界条件を考慮した補正後の速度ベクトルｕ^＊を、ヘルムホルツ・ホッジ分解定理により導入した中間的な速度ベクトルｗとして用いる。なお、境界近傍以外の各ノードでは、流体の速度ベクトルｕそのものをｗとして用いる。

次に、式（６）により中間的な速度ベクトルｗを計算する（Ｓ２）。求めた中間的な速度ベクトルｗを、式（５）の境界条件を満たすように補正する（Ｓ３）。境界条件を考慮した補正後の速度ベクトルｗ^＊を式（８）に適用し、スカラー値ｑを計算する（Ｓ４）。

次に、Ｓ４で求めたスカラー値ｑを、式（５）の境界条件を満たすように補正する（Ｓ５）。また、Ｓ５では、流れの３次元的な成分を圧力に関連するスカラー値ｑに反映させるための、２．５Ｄ補正も行う。２．５Ｄとは、２次元と３次元の中間である２．５次元のことである。３次元的な流れの影響を２次元の数値シミュレーションに組み込むことを、２．５次元と表現しているのである。シミュレーションはエコー信号データで得られた２次元の面内の情報に基づき行われるが、実際の流れはその面の中から外への流出及びその面の外から中への流入を含んでいる。そこで、２．５Ｄ補正では、そのような３次元的な流出や流入の影響を、圧力に対応するスカラー値ｑに反映させるのである。この補正により、３次元的な流れの影響が数値シミュレーションに組み込まれる。なお、２．５Ｄ補正については、後で詳しく説明する。

次に、境界条件補正及び２．５Ｄ補正の後のスカラー値ｑ^＊と速度ベクトルｗ^＊とを式（４）に代入することで、流体の速度ベクトルｕを計算する（Ｓ６）。また、スカラー値場ｑ^＊から式（７）により圧力場ｐを計算する（Ｓ７）。

そして、Ｓ６で求めた流体の速度ベクトルｕを、次の反復演算におけるＳ１の入力として用いて、Ｓ１〜Ｓ７の処理を繰り返す。以上に説明したＳ１〜Ｓ７の処理ループをあらかじめ定めた回数繰り返すことで、ビーム走査面内のビーム方向の速度成分のみしかない超音波ビームデータから、３次元的な流れの影響を反映しかつビームに直交する速度成分を含んだ流体の速度ベクトルｕと圧力ｑとが求められる。

＜数値シミュレーションの具体例＞
次に、以上に説明したアルゴリズムを具体的な数値計算に実装する方法について説明する。一つの実装例は、上述の式（４）〜（８）に現れる各項を有限差分形式で時間的及び空間的な微分係数を置き換えることである。この置き換えの例を表１に示す。

ただし、単純な置き換えだけでは必ずしも十分ではない。そこで、以下のような点を考慮する。

上述の式（６）は、時間ｔに関しては線形の偏微分方程式である。したがって、この式（６）は、重ね合わせの原理により次の式（９ａ）、（９ｂ）に分解することができる。

式（９ａ）はナビエ・ストークス方程式の移流項に、式（９ｂ）は粘性項に対応する。式（９ａ）、（９ｂ）をそれぞれ個別に解くことで得られる解ｗ’とｗ”を足し合わせることで、式（６）の解ｗを得ることができる。

Ｓ１，３，５での境界条件に基づく補正では、例えばＢモード画像から境界線（例えば心臓や血管の壁）を求め、その境界線について式（５）の境界条件を満たすようにすればよい。Ｂモード画像から境界線を求めるには、例えば、Ｂモード画像を二値化すればよい。適切な二値化のしきい値は、実験などにより求めればよい。単純な例では、図３に例示するように、注目画素（図３の各例における中央の画素）の上下左右の４画素がそれぞれ白と黒のいずれであるかの組み合わせに応じて、その注目画素における境界線の法線方向（注目画素から延びる矢印）を決定する。このように求めた法線方向に基づき、速度ベクトルｕ＝（ｕ_ｘ，ｕ_ｙ）、（圧力に相当する）スカラー値ｑを以下の表に示すように補正する。

例えば、法線方向が右斜め上４５度方向の場合、境界条件補正後の速度ベクトルｕ^＊＝（０．７０７ｕ_ｘ，−０．７０７ｕ_ｙ）とする。また、スカラー値ｑ^＊＝ｑ（ｘ−Δｘ，ｙ−Δｙ）とする。表２では流体の速度ベクトルｕについての例を示したが、中間的な速度ベクトルｗについてもｕと同じ補正を行えばよい。表２に例示した速度ベクトルの補正は、大きさを変えずに方向を境界線の法線方向へと変更するものである。また、圧力に対応するスカラー値ｑの境界線補正は、境界線上のノードのスカラー値を、境界線の法線方向について１格子分手前のノードのスカラー値ｑに合わせるものであり、これにより境界線の法線方向についてｑが変化しないという境界条件を満たすことができる。

なお、以上の境界線補正の例では、注目画素の４近傍の画素の組み合わせから境界線の法線方向を求めたが、これは一例に過ぎない。もっと多数の近傍画素の情報を用いれば、境界線やその法線方向をより厳密に求めることができる。ビットマップ画像から境界線及びその法線方向を求める技術は様々に知られており、例示した方式の代わりにそのような公知の技術を用いてももちろんよい。

また、更に、Ｓ５では、境界線近傍の各ノードについて、圧力に対応するスカラー値ｑに対して２．５Ｄ補正を行う。補正の対象とする境界線の「近傍」のノードは、境界線上のノードに限定してもよいし、流れの範囲の中で境界線からあらかじめ定められた距離以内にあるノードとしてもよい。境界線からその距離以内の範囲は、境界層と捉えることができる。

この２．５Ｄ補正では、例えば図４に示すような補正曲線Ｃを用いる。図４における横軸ｑ（２Ｄ）は入力であるスカラー値ｑであり、縦軸ｑ（２．５Ｄ）は出力（補正結果）のスカラー値ｑである。このグラフは、横軸及び縦軸についてそれぞれ対称となっている。対称軸である横軸及び縦軸は、スカラー値ｑの基準とする値であり、例えば０とする。上述の数値演算では、スカラー値ｑ（及び圧力ｐ）の勾配を計算に用いており、ｑ（及びｐ）の勾配は定数を足したり引いたりしても変わらないので、簡潔化のためにｑ＝０を基準としている。ただし、あくまでこれは一例に過ぎない。

ここで、スカラー値ｑは、後述するように、ある与えられた初期値（これは後述のようにエコー信号データから得られる速度情報に基づいて仮定してもよいし、エコー信号データに関わりなくあらかじめ定めた値を用いてもよい）から数値計算を時間ステップごとに反復実行することにより、徐々に流れの性質に即した正確な値へと近づいていく。数値計算では、前回の時間ステップ（ステップ（ｋ−１）とする）についての計算で求められた速度ベクトルやスカラー値ｑ（圧力に対応する）を入力として、数値計算式を解くことで、次の時間ステップｋの速度ベクトルやスカラー値ｑを求める。

この反復演算の中の、ある時間ステップｋ（すなわちｋ番目の計算タイミング）での数値計算式の実行の際に、その数値計算式に入力するスカラー値ｑに対し、補正曲線Ｃに従った２．５Ｄ補正を施す。すなわち、時間ステップｋの計算において数値計算式に入力するスカラー値ｑは、１つ前の時間ステップ（ｋ−１）で求められたスカラー値ｑであるが、このスカラー値ｑを横軸の値ｑ（２Ｄ）とし、図４の補正曲線Ｃ上で、その横軸値に対応する縦軸の値ｑ（２．５Ｄ）を補正結果とする。そして、その補正結果を時間ステップｋの計算にてスカラー値ｑの入力値と用いて、数値計算式を解くのである。

図４の補正曲線Ｃの横方向及び縦方向の対称軸（すなわち横軸と縦軸）の交点におけるスカラー値ｑは、横軸についても、縦軸についても、あらかじめ定めた基準値ｑ_ｍである。この基準値は、上述のように例えば０としてもよい。図４の補正曲線では、入力スカラー値ｑが基準値ｑ_ｍに近い範囲では補正結果としてその入力スカラー値ｑそのものを採用し、基準値ｑ_ｍからある程度以上離れると入力スカラー値ｑよりも基準値ｑ_ｍに近い値へと補正している。特に図４の例では、入力スカラー値ｑと基準値ｑ_ｍとの差が閾値以内の範囲では、入力スカラー値ｑ（２Ｄ）の上昇（下降）に応じて補正結果のスカラー値ｑ（２．５Ｄ）も比例的に上昇（下降）するが、入力スカラー値ｑと基準値ｑ_ｍとの差がその閾値を超えると、入力スカラー値ｑが上昇（下降）しても補正結果のスカラー値ｑの上昇（下降）は鈍り、徐々に頭打ちとなる。ある１つの例では、入力スカラー値ｑと基準値ｑ_ｍとの差が閾値を超えると、補正結果のスカラー値ｑを飽和させて一定値とする。

この補正は、以下のような考え方に基づく。すなわち、一般に、血流等の流体が血管壁に当たると、その当たった部分の圧力は周囲に比べて高くなる。また、壁に当たった流体は、壁を通り抜けることはなく、壁に沿って他の方向に移動するが、このときの流体の移動は超音波ビームの走査面内に限られるものではなく、走査面の外へも３次元的に移動していく。すなわち、走査面外に流出する。逆に、壁の圧力が周囲に比べて低い場合、周囲から流体が流入するが、その流入は走査面内だけからでなく面外からも到来する。この２．５Ｄ補正では、壁上での圧力と周囲の圧力（すなわち基準値）との差が大きくなるほど、このような走査面の外との間の流出及び流入が多くなり、その圧力差が緩和されると仮定するのである。なお、上述する数値計算では、この補正は、圧力ｐそのものに対してではなく、その圧力に対応するスカラー値ｑに適用している。

図４の例では、入力スカラー値ｑがある閾値以上になると２．５Ｄ補正値を一定の値としている。これは、閾値値を超える圧力に相当する流れは超音波ビーム走査面の外に流出するとの想定に基づく。なお、図４に例示した曲線はあくまで一例に過ぎない。例えば、２．５Ｄ補正値は完全には飽和せず（すなわち一定値にせず）に少しずつ増えるようにしてもよい。

実際の超音波診断装置への実装では、２．５Ｄ補正のための変換曲線（図４）は、２Ｄでのｑ値とこれに対応する２．５Ｄでのｑ値とを記憶したテーブル、又は２Ｄでのｑ値を２．５Ｄでのｑ値に変換するための関数を表すプログラム、などの形で、超音波診断装置の制御部に記憶させておけばよい。

以上の例では、流れを規制する壁を表す境界線をＢモード画像から求めたが、これは一例に過ぎない。この代わりに、各点でのドプラ信号の強さを示すパワーフロー画像データから境界線を求めてもよい。パワーフロー画像を例えば二値化することで境界線を求めるなどである。例えば、実際の心臓内の流れでは、超音波ドプラの感度より小さい速度は測定できないので、Ｂモード画像で求められる心臓壁よりも少し内側の領域でしかドプラ速度のデータが出ない。これに対し、パワーフロー画像から境界線を求めることで、超音波ドプラの感度のある部分についてのドプラ情報に基づいた計算ができる（すなわちドプラ情報のない領域については計算をしない）ので、シミュレーションの精度が向上する。

以上の例では、境界線（壁）に該当するノードについて速度ベクトルや圧力（スカラー値ｑ）を補正したが、境界線だけでなく、境界線からある厚みを持った境界層内に属するノードについても上述の境界線補正や２．５Ｄ補正を行ってもよい。境界層の厚みはあらかじめ決定し、装置に設定しておく。この場合、補正の程度は、境界線上のノードの補正の程度を最大値とし、境界線からの距離が離れるほど補正の程度を緩和するようにしてもよい。また、境界線上のノードについての補正結果と、境界層の外側にある境界層に最も近いノードでの速度ベクトルやスカラー値ｑとから補間を行うことにより、境界層内のノードの速度ベクトルやスカラー値ｑの補正値を求めてもよい。

以上、境界近傍のノードについての速度、及び圧力に対応するスカラー値ｑについての補正について説明した。

図２の説明に戻ると、Ｓ２では、速度の補正結果を用いて移流項及び粘性項の数値計算を行う。このＳ２の計算の具体例を、以下に示す。

移流項の数値計算は、式（９ａ）の各項を単純に有限差分形式に置き換えて数値的に解くことで計算できる。また、別の例として、セミ・ラグランジュ法を用いることもできる。

セミ・ラグランジュ法では、現在の速度ベクトルｗ’が微小な時間ステップΔｔだけ後（又は前）の時刻でも保存されるものと仮定とする。そして、現在（時刻＝ｔ）のデータ点ｒ（ｒは位置ベクトル）での次の時刻（＝ｔ＋Δｔ）の速度ベクトルｗ’（ｒ，ｔ＋Δｔ）が、現在時刻（＝ｔ）においてそのデータ点ｒからその速度ベクトルｗ’×Δｔだけ戻った位置（ｒ−ｗ’Δｔ）での速度ベクトルｗ’（ｒ−ｗ’Δｔ，ｔ）と等しいと考える。すなわち、ｗ’（ｒ，ｔ＋Δｔ）＝ｗ’（ｒ−ｗ’Δｔ，ｔ）という式に従って、次の時刻ｔ＋Δｔにおける当該データ点ｒの速度ｗ’（ｒ，ｔ＋Δｔ）を計算する。ここで、データ点ｒとして流体モデルにおける各ノード１０を用いた場合、その点ｒからｗ’Δｔだけ戻った位置（ｒ−ｗ’Δｔ）は、一般に流体モデルのノード１０とは一致しないので、速度ベクトルの値が存在しない。そこで、ｗ’（ｒ−ｗ’Δｔ，ｔ）の値は、当該位置（ｒ−ｗ’Δｔ）の近傍の複数のノード１０における現在時刻ｔでの速度ベクトルｗ’を補間することにより求めればよい。補間は、公知のどのような補間法を用いてもよい。

粘性項は、式（９ｂ）を、次のように有限差分形式に書き換えて解けばよい。

この式のｘ成分を取り上げると次のようになる。

ここで、Δｘ、Δｙは図１に示した流体モデルにおけるｘ、ｙ方向についてのノード（格子点）の間隔ｄｘ、ｙのことである。ｙ成分も同様に有限差分形式で表現できる。

このような形式に直すことで、ヤコビ法による反復解法（次式）が利用可能になる。

この反復解法では、初期値ｗ_ｘ ^（０）として現在時刻ｔにおける速度場ｗ_ｘ（ｔ）を用い、この初期値から計算を開始してあらかじめ定めた回数だけ反復を繰り返す。現在時刻ｔにおける速度場ｗ_ｘ（ｔ）は、上述の境界線補正（Ｓ３）後の速度場ｗのことである。また、反復は、上記移流項において用いた時間ステップΔｔごとに行う。

なお、この反復解法の式はｘ成分についてのものである。ｙ成分についても同様の計算を行えばよい。これにより、粘性項の成分が求められる。

なお、この粘性項及び前述の移流項の反復演算は、１回ごとに、圧力に対応するスカラー値ｑや実際の速度ベクトルｕについての前回の反復演算の結果を反映した形で行う。なお、この点は、圧力に対応するスカラー値ｑや実際の速度ベクトルｕについての反復演算（後述）についても同様である。

このように計算された移流項の速度成分ｗ’と粘性項の速度成分ｗ“とを足し合わせることで、中間的な速度ベクトルｗが求められる。

なお、反復の各回で求められたベクトルｗは、Ｓ３で境界条件に基づく補正を受け、Ｓ４におけるスカラー値ｑの演算に用いられる。

次に圧力に対応するスカラー値ｑについての数値演算（Ｓ４）の具体例について説明する。既に上述のようにして中間的な速度ベクトルｗが求められているので、スカラー値ｑは上述の式（８）から計算できる。

式（８）の右辺は、有限差分方式に従って展開すると、次のようになる（表１参照）。

同様に式（８）の左辺は次のように展開される。

このように展開すると、ヤコビ法による反復解法（次式）が利用可能になる。

この反復解法では、初期値ｑ^（０）（ｘ、ｙ）して現在時刻ｔにおけるｄ（ｘ、ｙ）の値を用い、この初期値から計算を開始してあらかじめ定めた回数だけ反復を繰り返す。反復の各回の演算では、ｄ（ｘ、ｙ）の値は、その回のｗについての反復演算の結果（Ｓ３で境界条件についての補正済）に基づき求めればよい。求めた今回のｑの値は、Ｓ５で境界条件補正及び２．５Ｄ補正を施された上で、次のＳ６及びＳ７の演算に供される。

次に、速度ベクトルｕの計算（Ｓ６）の具体例について説明する。以上の演算で、ｗ及びｑが求められたので、式（４）から速度ベクトルｕを求めることができる。なお、式（４）の適用のためには、スカラー値ｑの勾配を求める必要がある。これは次式に従って行えばよい。

ここで、ｇ_ｘ、ｇ_ｙは、それぞれｑの勾配のｘ、ｙ成分である。ここで、右辺の各項のｑの値は、今回の反復におけるＳ５の補正結果である。

速度ベクトルｕのｘ、ｙ成分ｕ_ｘ、ｕ_ｙは次式で計算すればよい。

次に、実際の圧力値ｐの計算（Ｓ７）の具体例について説明する。反復演算における今回のスカラー値ｑの演算結果はＳ４及びＳ５で既に求められているので、この演算結果を式（７）に適用することで、圧力値を計算することができる。実際の数値計算では、式（７）を次式の有限差分形式に展開して計算すればよい。

以上に例示したｗ，ｑ，ｕ，ｐについての数値計算が、時間ステップΔｔごとの計算時点について反復実行される。この反復演算の初期値は、超音波診断装置により取得したあるサンプリング時刻のビーム方向の速度成分（ただし、境界付近については補正済）であるが、ナビエ・ストークス方程式及び境界条件を折り込んだ反復演算により、最終的に、流体の性質に即してビームに直交する成分を含んだ速度ベクトルｕが求められることになる。また、圧力に対応するスカラー値ｑに対して、３次元的な流れの影響を折り込むための２．５Ｄ補正を施すことで、３次元的な流れの影響を加味した速度ベクトルｕや圧力値ｐが求められる。

超音波診断装置におけるあるサンプリング時刻から次のサンプリング時刻までの間に、以上のような所定回数の反復計算を実行することで、超音波診断装置から得られる超音波ビーム方向の速度情報から、３次元的な流れの影響を加味した２次元速度ベクトルｕや圧力値ｐを、リアルタイムで計算することができる。

なお、初期値から開始して上述の所定回数の反復計算を経て得られた速度ベクトルｕの超音波ビーム方向の成分は、初期値である超音波診断装置から得た生のビーム方向速度成分とは一般に一致しない。ここで、超音波診断装置から得た生のビーム方向速度成分は実測値なので尊重するとの立場を取るならば、反復計算により得た速度ベクトルｕを生のビーム方向速度成分により補正することとなる。ここで、補正の仕方としては、次の２つが考えられる。１つは、反復計算により得た速度ベクトルｕのベクトル方向を維持したまま、その速度ベクトルｕのビーム方向の速度成分が超音波診断装置から得た生のビーム方向速度成分と等しくなるよう、速度ベクトルｕを拡大又は縮小する補正である。この補正では、補正後の速度ベクトルｕのビームに直交する方向の速度成分が補正前とは異なったものとなる。もう１つは、反復計算により得た速度ベクトルｕのうちビーム方向に直交する速度成分は保存し、ビーム方向の速度成分を超音波診断装置から得た生のビーム方向速度成分に置き換えるという補正である。この補正の場合、速度ベクトルｕの方向が、補正前とは異なったものとなる。

以上に例示した数値計算方法の例は、直交座標系の流体モデルに基づいている。したがって、セクタ走査やコンベックス走査などのような非直交座標系の超音波ビーム走査により得た超音波ビームデータについては、公知の方法で直交座標系へ座標変換してから上記の数値計算手法を適用すればよい。この場合、数値計算により求めた２次元速度ベクトルｕや圧力ｐは、その座標変換の逆変換を行うことで、元の非直交座標系での値に変換することができる。また、セクタ走査などの非直交座標系で流体モデルを構築し、その座標系での数値計算式に従って数値計算を行うようにしてももちろんよい。

＜画面表示＞
（１）流速ベクトル表示
数値シミュレーションにより求めた各ノードにおける流速ベクトルｕを、例えばそのノードから延びる線分として画面表示してもよい。この場合、線分の向き及び大きさがベクトルｕの向きと大きさを表すようにすればよい。この流速ベクトル表示は、Ｂモード画像やカラードプラ画像に重畳して表示してもよい。

（２）流体粒子移動表示
また、直感的に分かりやすくするために、数値シミュレーションにより求めた２次元の流速に従って流体粒子が移動する様子をアニメーション表示してもよい。

この表示では、例えば表示の第１フレームでは、流体モデルの各ノード（数が多すぎる場合は間引いてもよい）につき、当該ノードを表す粒子画像（所定数の画素からなる塊であってもよい）を、シミュレーションで求めた当該ノードの２次元流速ベクトルに応じた方向及び移動量だけ当該ノードの本来の位置から移動した位置に表示する。第２フレームでは、第１フレームでの移動先の位置にある各粒子画像を、第２フレームでのシミュレーションで求めた当該位置での２次元流速ベクトルに応じた方向及び移動量だけ、その位置から移動させる。なお、第１フレームでの粒子の移動先の位置における、第２フレームでの流速ベクトルは、例えば、シミュレーションにより求められた流体モデルの各ノードについての２次元流速ベクトルからの補間により求めればよい。このようにして、第３フレーム以降も同様に、前フレームでの移動先の粒子を、現フレームでのシミュレーション結果の流速ベクトルに応じて移動させる。このような処理により、流体粒子が流れていくようすを直感的に表示する動画像を生成することができる。

（３）流速マップ表示
従来の超音波診断装置の速度表示は、超音波プローブに近づく方向と遠ざかる方向の二方向についての情報しかないため、近づく方向と遠ざかる方向にそれぞれ異なる色（例えば赤と青）を割り当てて表示していた。これに対し、本実施形態では、数値シミュレーションによりビーム方向に直交する速度成分を持つ流速ベクトルｕが求められるので、より豊かな情報を表示することができる。一例として、図５に示すように、流速ベクトルｕの方向を色相で、大きさを明度で示す色変換マップを用いればよい。すなわち、各ノードの流速ベクトルｕを、この色変換マップに従って色値に変換し、この色値に応じた色を当該ノードの色として画面表示すればよい。なお、色変換マップは、例えば流速ベクトルｕのｘ、ｙ成分の組み合わせごとに、その組み合わせに対応する色値（例えば色相と明度の組、又はこの組が表す色をＲＧＢの成分で表現したデータ）を登録したデータとして実装すればよい。

以上の例では、流速ベクトルｕの大きさを表示色の明度に変換したが、明度の代わりに、濃度や彩度などのような、色相以外の色パラメータに変換してもよい。

（４）圧力表示
数値シミュレーションにより求めた各ノードの圧力ｐを、例えば色で表示することができる。例えば、圧力の平均値を黒とし、その平均より高い圧力を持つノードを青色で、平均より低い圧力を持つノードを赤色で表示するなどである。また、圧力の平均からの差の大きさに応じて、当該ノードの明度を変化させてもよい（例えば、平均からの差が大きいほど明るくするなど）。

（５）動圧表示
数値シミュレーションにより求めた各ノードの流速ベクトルｕを元に各ノードの動圧を計算し、表示してもよい。動圧は、ρｕ^２／２である。ρは既知、ｕは数値シミュレーションで計算されているので、これらから動圧を計算することができる。動圧は、上述の圧力ｐと同様の表示形態で表示することができる。

（６）渦度表示
数値シミュレーションにより求めた各ノードの流速ベクトルｕの回転（すなわち∇×ｕ。ｒｏｔ ｕとも表される）を渦度として計算し、渦度を画面表示してもよい。ベクトルｕの回転はベクトルであるが、この回転ベクトルを渦度として計算するのである。計算される渦度はベクトルであり向き（すなわち回転の方向）がある。例えば、この回転の方向が時計回りか反時計回りで色分け（例えば前者を青、後者を黄色で表すなど）し、渦度の大きさ（すなわちベクトルの絶対値）を明度などと行った色相以外のパラメータの値に反映させることで、表示色を求めてもよい。

以上の（３）〜（６）の表示方式の説明では、流体モデルにおける各ノードについての表示色（画素値）の求め方を例示したが、ノード同士の間の各画素の画素値を、当該画素の近傍のノードの画素値から補間により求め、その補間結果を表示するようにしてもよい。

（７）圧力変化表示
図６に示すように、ユーザから指定された注目点の圧力ｐ（或いは注目範囲の平均圧力）の時系列的な変化のグラフ１００を画面に表示してもよい。注目点（又は注目範囲）は、例えばＢモード表示がカラードプラ表示の画面上で、マウスやトラックボールなどのポインティングデバイスを用いてユーザが指定すればよい。図６のグラフ１００は、心房内の注目点の圧力の時間変化の例である。このグラフは、縦軸が圧力ｐ、横軸が時間ｔである。医師は、このグラフ１００の勾配（時間変化率）から、心臓の弁機能を評価することができる。

また、注目点の圧力ｐの大きさに応じた音量或いは周波数の音を出力することで、ユーザが音の時間的な変化により圧力ｐの時間的な変化を認識できるようにしてもよい。

また、図７に示すように、ユーザが指定した２つの注目点（又は２つの注目範囲）の圧力の時系列変化のグラフ１００及び１０２を同時に表示してもよい。例えば、グラフ１００は心房内の注目点の圧力の変化を示し、グラフ１０２は心室内の注目点の圧力の時間変化を示す。図７のように、心房と心室の圧力の時間変化を、同一時間・同一スケールで並べて表示することで、心房と心室との圧力差（グラフ１００と１０２の間隔）を視覚的に把握することができる。

＜超音波診断装置の構成例＞
上記の方法を実装した超音波診断装置の装置構成の例を図８に示す。図８に示した装置構成において、数値シミュレーション部２１４が前述の数値演算処理により流速ベクトルｕ及び圧力ｐを計算する機能モジュールであり、他の構成要素は、従来一般的な超音波診断装置の構成要素と同等のものでよい。

この装置では、操作パネル２００からの操作に応じた制御部２０２からの制御に従い送信部２０４にて送信信号を生成し、この送信信号によりプローブ（超音波探触子）２０６中の振動子アレイを駆動して、被検体３００内に超音波パルスのビームを送信する。そして、このビームの被検体３００内からのエコーをプローブ２０６で受信し、電気的な受信信号を生成する。この受信信号は受信部２０８で所定の信号処理を受けた上で、Ｂモード演算部２１０、カラードプラ演算部２１２に入力される。これらは、Ｂモードデータ、カラードプラデータ（すなわち各点のビーム方向速度成分）をそれぞれ計算する。

数値シミュレーション部２１４は、カラードプラ演算部２１２で求められた各点のビーム方向速度成分の値から、前述の数値シミュレーション処理を実行することにより、それら各点の２次元流速ベクトルｕと圧力ｐを計算する。すなわち、数値シミュレーション部２１４は、境界検出部２５０と反復演算部２５２を有する。境界検出部２５０が、あるフレームにおけるＢモード画像、又はカラードプラ演算部２１２が求めた速度情報から求めたパワーフロー画像から、血流範囲の境界線を特定し、境界上の各点の法線方向を求める。反復演算部２５２は、当該フレームのカラードプラ画像を出発点として、上述の移流項、粘性項、スカラー値ｑ、速度ベクトルｕ及び圧力値ｐの数値計算式を反復実行する。１回の反復ごとに、前回の反復における計算結果に対して、境界に関する補正、すなわち境界条件を満たすための補正と圧力（ｑ）についての２．５Ｄ補正とを施し、その補正結果を入力して上述した数値計算式の演算処理を実行する。この反復演算をあらかじめ定めた回数だけ繰り返すことで、当該フレームについての２次元流速ベクトルｕと圧力ｐを計算する。このような数値シミュレーション処理を１フレームごとに行うことで、２次元流速ベクトルｕと圧力ｐをリアルタイムで計算することができる。数値シミュレーション部２１４は、ソフトウエア処理で実現することもできるが、上記方法に基づく演算処理を例えばＤＳＰ（デジタルシグナルプロセッサ）等を用いて回路的に実装してもよい。

画像生成部２１６は、このように求められたＢモードデータ、カラードプラデータ、並びに２次元速度ベクトルｕ及び圧力ｐのデータを用いて診断用の表示画像を形成する。数値シミュレーションにより求めた２次元速度ベクトルｕ及び圧力ｐに基づく表示画像としては、上述の（１）〜（７）などがある。なおこの画像の生成のために、Ｂモードデータやカラードプラデータなどを必要に応じて走査変換してもよい。例えば、超音波ビーム走査がセクタ走査などのように直交座標系にしたがったものでない場合、画像表示のために直交座標系に走査変換すればよい。走査変換により表示部２２０の画像フォーマットに合わせられた各種診断画像が、診断画像モードの設定に従って適宜組み合わされることで、表示画像が形成され、表示部２２０に表示される。この表示では、上述の（１）〜（７）の画像のうちの１つ又は複数の組み合わせを、Ｂモード画像やカラードプラ画像に重ねて表示することもできる。

１０ ノード、２００ 操作パネル、２０２ 制御部、２０４ 送信部、２０６ プローブ、２０８ 受信部、２１０ Ｂモード演算部、２１２ カラードプラ演算部、２１４ 数値シミュレーション部、２１６ 画像生成部、２２０ 表示部、２５０ 境界検出部、２５２ 反復演算部。

Claims (5)

1. 超音波ビームを被検体内に送受信し、被検体内からのエコー信号を取得する送受信手段と、
   前記エコー信号に基づき、前記被検体内の各点の超音波ビーム方向についての速度情報を求めるドプラ処理手段と、
   前記エコー信号に基づき、前記被検体内の流体が流れる流れ範囲の境界を検出する境界検出手段と、
   前記境界検出手段が検出した前記境界についての境界条件の下で、前記被検体内の各計算点の超音波ビーム方向についての速度情報を初期値として、ナビエ・ストークス方程式に基づく数値計算式を時間ステップごとに反復演算することにより、前記各計算点について、前記超音波ビーム方向の速度成分と前記超音波ビーム方向に直交する方向の速度成分とを含む速度ベクトル及び圧力、を計算する反復演算手段と、
   前記反復演算手段が計算した前記各計算点についての前記速度ベクトル又は前記圧力を表示するための処理を行う表示処理手段と、
   を備え、
   前記反復演算手段は、前記境界の近傍の計算点については、前記数値演算式における圧力に対応する項の入力値に対し、基準値からその入力値までの差が閾値以上の場合にはその差を緩和する補正を行った上で前記数値計算式を反復演算する、
   ことを特徴とする超音波診断装置。
2. 請求項１に記載の超音波診断装置において、前記境界検出手段が、前記エコー信号から求められるパワーフロー画像に基づき、前記境界を検出する、ことを特徴とする超音波診断装置。
3. 請求項１又は２に記載の超音波診断装置において、前記表示処理手段は、前記速度ベクトルを表示する速度表示画像における前記各計算点の表示色の色相を当該計算点における前記速度ベクトルの方向に基づき求め、前記速度表示画像における前記各計算点の表示色の色相以外のあらかじめ定められた色パラメータを当該計算点における前記速度ベクトルの大きさに基づき求める、ことを特徴とする超音波診断装置。
4. 請求項１〜３のいずれか１項に記載の超音波診断装置において、
   前記反復演算手段により求められた前記各計算点における前記速度ベクトルから、それら各計算点の渦度を計算する渦度計算手段、を更に備え、
   前記表示処理手段は、前記各計算点の渦度が表す回転の向きを色分け表示した渦度表示画像を生成する、
   ことを特徴とする超音波診断装置。
5. 被検体内からのエコー信号に基づき、前記被検体内の各点の超音波ビーム方向についての速度情報を求めるステップと、
   前記エコー信号に基づき、前記被検体内の流体が流れる流れ範囲の境界を検出する境界検出ステップと、
   前記境界検出ステップで検出された前記境界についての境界条件の下で、前記被検体内の各計算点の超音波ビーム方向についての速度情報を初期値として、ナビエ・ストークス方程式に基づく数値計算式を時間ステップごとに反復演算することにより、前記各計算点について、前記超音波ビーム方向の速度成分と前記超音波ビーム方向に直交する方向の速度成分とを含む速度ベクトル及び圧力、を計算する反復演算ステップと、
   前記反復演算ステップで計算された前記各計算点についての前記速度ベクトル又は前記圧力を表示するための処理を行う表示処理ステップと、
   を備え、
   前記反復演算ステップでは、前記境界の近傍の計算点については、前記数値演算式における圧力に対応する項の入力値に対し、基準値からその入力値までの差が閾値以上の場合にはその差を緩和する補正を行った上で前記数値計算式を反復演算する、
   ことを特徴とする超音波診断装置における数値シミュレーション方法。
   

Images