WO2004010872A1

WO2004010872A1 - 超音波診断システム及び歪み分布表示方法

Info

Publication number: WO2004010872A1
Application number: PCT/JP2003/009731
Authority: WO
Inventors: Tsuyoshi Shiina; Makoto Yamakawa; Naotaka Nitta
Original assignee: Hitachi Medical Corporation
Priority date: 2002-07-31
Filing date: 2003-07-31
Publication date: 2004-02-05
Also published as: CN100475154C; EP1541090A4; US8041415B2; EP2374413A2; EP2201896A1; EP1541090A1; EP2374413B1; CN1678243A; US20110125019A1; EP2374413A3; EP1541090B1; US8382670B2; US20060052696A1

Abstract

被検体の圧縮前後に超音波探触子から出力されるRF信号を直交検波して得られる包絡線信号からなるフレームデータに計測点を設定し、計測点を前記フレームデータに対して超音波ビーム方向(1次元方向)、2次元又は3次元方向に移動して、計測点を囲む相関窓に属する圧縮前後の包絡線信号の相関係数が最大となる位置を求め、これに基づいて圧縮に伴う計測点の変位を求めるとともに、さらに圧縮前後のRF信号間の位相差を求めて圧縮に伴う計測点の変位を精度よく求める変位演算手段とを備えて構成することにより、変位量の制約を受けずに変位分布を推定することができ、演算処理時間を短縮し、横方向変位にも対応できるようにする。

Description

P T/JP2003/009731

明細書超音波診断システム及び歪み分布表示方法技術分野

本発明は、超音波診断装置を用いて、生体組織の硬さを定量的に計測することのできる超音波診断システム及び歪み分布表示方法に関する。背景技術

超音波の医療面への応用は、組織形状だけではなく組織内の音速や減衰定数などの物理量を超音波により計測 ·画像化し診断に利用しょうとする分野 (ultrasonic tissue characterization) がある。その中の 1つとして組織の硬さ、すなわち弾性特性を計測しょうとする分野があり、現在盛んに研究されている。これは、組織の弾性特性がその病理状態と深く関連しているためである。例えば、乳がんや甲状腺がんなどの硬化性がんや肝硬変、動脈硬化症などは正常組織よりも病変部分が硬くなることが知られている。これまで、生体組織の硬さ情報は触診により得ていたが、触診では客観的な情報表現が難しく、医師の経験も必要で、また計測できる領域も体表付近のある程度大きな病変に限られる。

そこで、体表から静的な圧力を加えて組織をわずかに圧縮変形させ、その際生じる組織内部の歪みを超音波により計測し、歪みから組織の弾性特性を評価する方法力 sある (J. Op ir, I. Cespedes, H. Ponnekanti , Y. Yazdi and X. Li ,

'Elastography: A quantitative method for imaging the elasticity of biological tissues", Ultrasonic Imaging, vol. 13, pp. Ill - 134, 1991)。この従来技術は、硬い組織では生じる歪みが小さく、軟らかい組織では歪みが大きくなることを基にしている。すなわち、組織を静的に圧縮してその際生じる組織内の歪み分布から組織の弾性特性を評価する。

具体的には、超音波診断装置および超音波プローブをそのまま用い、超音波プローブによって組織を圧縮する前の超音波エコー信号（圧縮前 RF信号）を計測 2 する。その後、例えば超音波プローブにより組織をわずかに（数パーセント程度）圧縮し、組織の圧縮後の超音波エコー信号（圧縮後 R F信号）を計測する。そして、計測された組織圧縮前後の R F信号から圧縮によつて組織内部の各点がどれだけ動いたかという変位量である変位分布を空間相関法により推定する。

空間相関法は、圧縮によって生じた組織内部の変位分布を組織圧縮前後の R F 信号（または R F信号の包絡線）から 2次元相関関数を用いたテンプレートマツチングにより推定する手法である。つまり、断層像データに対応する圧縮前の R F信号デ一夕に一定サイズの 2次元の相関窓（テンプレート）を設定し、その相関窓の R F信号データと最大の相関を有する圧縮後の R F信号データの 2次元位置に、着目した計測点が移動したことを自己相関処理により推定する手法である。この自己相関処理を例えば格子状に設定される各計測点について行い変位分布推定する。空間相関法の処理において、一般に、横方向（超音波ビームのスキャン方向）のサンプリング間隔は、軸方向のサンプリング間隔よりも大きいので、推定される変位成分の精度は横方向成分の方が軸方向成分よりも悪くなる。このように、空間相関法では 2次元の変位べクトル成分を推定できるという特徴がある。また、変位推定精度がサンプリング間隔によって制限されてしまうが、組織が大きく変形した場合（例えば、 5 %程度）でも変位分布を推定できる。しかし、空間相関処理にかかる計算量が膨大になるため超音波計測の利点であるリアルタイム性を損なつてしまうという問題がある。

本発明の目的は、リアルタイムに変位分布、歪み分布、弾性係数分布を求めることにある。

発明の開示

上記目的を達成するため、本発明の超音波診断システムは、超音波探触子によつて被検体の圧縮前後に計測してなる反射エコー信号（R F信号）に基づいて被検体の組織の変位を求めるにあたり、超音波探触子から出力される R F信号の包絡線信号を含む特徴量を蓄積する記憶手段と、この記憶手段に蓄積された被検体の圧縮前後の特徴量に基づき圧縮前後の特徴量の相関係数及び圧縮前後の受信信号間の位相差を求める相関演算手段と、この相関演算手段によって求められた相関係数及び圧縮前後の R F信号間の位相差に基づいて圧縮に伴う計測点の変位を求める変位演算手段とを備えたことを特徴とする。また、計測点における変位を空間微分することによつて被検体の組織の歪み分布を求める歪み演算手段と、求めた歪み分布を表示する表示手段とを備えることができる。

このように、圧縮前後の包絡線信号などの特徴量による相関を用いて計測点の変位を求めていることから、リアルタイムに変位分布を推定することができる。特に、計測点を超音波ビーム方向に超音波信号の 2分の 1波長間隔で移動して相関係数の最大位置を求めることにより、ドプラ法におけるエイリアシングの問題を解決できる。

ところで、相関演算において、計測点における圧縮前後の包絡線信号の相関係数が最大となる位置を求めるに際し、圧縮後の包絡線信号の時間軸をずらしながら圧縮後の包絡線信号の自己相関関数を計測点ごとに求めると、演算時間がかかりすぎてリアルタイム性が損なわれるおそれがある。

そこで、圧縮前後の包絡線信号の自己相関関数を先に求め、求めた自己相関関数を計測点の移動に合わせて、例えば超音波信号の 2分の 1波長間隔でシフトして相関係数を求めることが好ましい。これにより、変位演算にかかる時間を短縮して処理を高速化できる。

また、本発明の超音波診断システムは、直交検波された R F信号の包絡線信号を蓄積する記憶手段と、この記憶手段に蓄積されたスライス面に対応する前記被検体の圧縮前後の前記包絡線信号のフレームデータに計測点を設定し、該計測点を前記フレームデータに対して少なくとも 2次元方向に移動して、前記計測点を囲む 2次元相関窓に属する前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるとともに、前記圧縮前後の前記 R F信号間の位相差を求める相関演算手段と、前記相関演算手段によって求められた前記相関係数が最大となる位置及び前記位相差に基づいて前記圧縮に伴う前記計測点の少なくとも 2次元方向の変位を求める変位演算手段とを備えることを特徴とする。

すなわち、相関窓を用いた空間相関法の長所である 2次元又は 3次元の変位計測に対応でき、ドプラ法の長所であるリアルタイム性と高精度を生かして組み合わせた複合自己相関法（C A法： Comb ined Autocorrelat ion method) と称する方法を提案するものである。この複合自己相関法によれば、一定の横方向変位にも対応して変位分布を推定可能である。この場合の、 2次元方向は、超音波探触子により受信される超音波ビームの方向と、この超音波ビームの走査方向に設定する。特に、計測点を超音波ビーム方向に超音波信号の 2分の 1波長間隔で移動し、超音波ビームの走査方向に超音波ビームピッチで移動して、相関係数の最大位置を求めることが好ましい。ここで、計測点の超音波ビーム方向の移動ピッチは、超音波信号の 2分の 1波長間隔に限られるものではないが、これよりも小さいピッチが好ましい。

さらに変位演算を高速化するには、上述したと同様に、計測点を囲む 2次元相関窓に属する圧縮前後の包絡線信号の相関係数が最大となる位置を求めるに際し、圧縮前後の包絡線信号の自己相関関数を求め、この自己相関関数を計測点の移動に合わせてずらして相関係数が最大となる位置を求めることが望ましい。

また、本発明の変位演算は、 2次元に限らず 3次元に適用できる。 1次元配列の超音波探触子を用いる 3次元の場合は、記憶手段に蓄積されるフレームデータを複数のスライス面のフレームデータを有するポリュ一ムデ—夕とする。また、 2次元配列の超音波探触子を用いる 3次元の場合は、スライス方向の走査により得られる包絡線信号を加えたものとなる。そして、相関演算手段は、計測点をポリユームデータに対して 3次元方向に移動して、計測点を囲む 3次元相関窓に属する圧縮前後の包絡線信号の相関係数が最大となる位置を求めるとともに、圧縮前後の R F信号間の位相差を求める。この場合の、 3次元方向は、超音波探触子により受信される超音波ビームの方向と、超音波ビームの走査方向と、これらに直交するスライス方向である。また、相関演算手段は、圧縮前後の R F信号間の超音波ビーム方向、超音波ビームの走査方向及びこれらに直交するスライス方向の位相差を求めることが好ましい。さらに、 3次元変位を求める場合も、上述した高速処理法を適用することができる。また、計測点をスライス方向に移動する場合は、超音波ビームのスライスピッチ単位にする。

また、弾性率分布を求める方法は、被検体を有限個の要素に分割して少なくとも 2次元又は 3次元の有限要素モデルを作成し、そのモデル化の情報と求めた歪み分布を用いて弾性係数分布を演算する弾性係数演算手段を備えることができる。また、これにより求めた弾性分布を表示手段に表示することができる。この場合、弾性係数演算手段は、被検体組織を等方性弾性体及び近非圧縮性と仮定し、被検体組織を有限個の直方体要素に分割して 3次元有限要素モデル化し、各要素内では、弾性係数、応力、歪みは一様であると仮定し、弾性方程式に前記歪み分布情報を用いて弾性係数分布を演算することが好ましい。

ここで、組織を等方性弾性体と仮定するのは、外部から圧力を加えて組織を静的に圧縮した場合、応力と歪みの間の関係はほぼ線形であり、組織を弾性体として近似でき、被検体の組織はほぼ等方性が成り立つので、この発明では組織を等方性弾性体と仮定している。また、組織を近非圧縮性と仮定するのは、生体組織が非圧縮性（ポアソン比レ = 0 . 5 ) であると特殊な弾性方程式となり、有限要素法を適用することができなくなるからである。また、ポアソン比を生体内で一定とすることで弾性係数分布推定の推定パラメ一夕をヤング率のみとすることができ、逆問題を簡単化できる。また、ポアソン比はヤング率に比べ生体中であまり変化しないパラメ一夕であるため、この発明ではポアソン比を 0 . 4 9で一定とすることが好ましい。この弾性係数分布演算によれば、精度よく演算可能な軸方向の歪み分布のみから弾性係数分布を再構成することができ、安定した弾性係数分布の演算が行える。図面の簡単な説明

図 1は、超音波診断装置の原理を説明するための図である。

図 2は、静的圧縮による組織弾性計測方式の具体例及び静的圧縮による組織弹性計測方式の原理を示す図である。

図 3は、空間相関法の原理を示す図である。

図 4は、ドプラ法の原理を示す図である。

図 5は、複合自己相関法の原理を示す図である。

図 6は、複合自己相関法の基本アルゴリズムを実行する回路構成を示すブロック図である。

図 7は、本発明の一実施例である超音波診断システムの概略構成を示すブロック図である。

図 8は、 3次元複合自己相関法の基本アルゴリズムを示すフローチャート図である。

図 9は、本発明の超音波診断システムに係る 3次元複合自己相関法の基本アルゴリズムを示すフローチャート図であり、図 7の処理の一部の詳細を示すフローチヤ一卜図である。

図 1 0は、図 9のステップ S 1 5の高速化された複合自己相関法の詳細を示すフローチャート図である。

図 1 1は、本発明の超音波診断システムに係る 3次元複合自己相関法の基本ァルゴリズムを実行する回路構成を示すブロック図である。

図 1 2は、シミュレーシ aンの手順の概略を示す図である。

図 1 3は、超音波中心周波数 5 . 0 MH zの場合に用いた各点広がり関数の一例を示す図である。

図 1 4は、組織モデルの概略を示す図である。

図 1 5は、横方向変位に対する各手法の歪み推定誤差を示す図である。

図 1 6は、横方向変位が 0 . O mmの場合の各手法（複合自己相関法、拡張複合自己相関法、空間相関法）により推定した歪み分布を示す図である。

図 1 7は、横方向変位が 0 . 4 mmの場合の各手法（複合自己相関法、拡張複合自己相関法、空間相関法）により推定した歪み分布を示す図である。

図 1 8は、斜め方向圧縮の検証に使用される組織モデルの概略を示す図である。図 1 9は、組織モデルを単純に軸方向に圧縮した場合の歪み分布推定結果を示す。

図 2 0は、組織モデルを斜め方向に圧縮した場合の歪み分布推定結果を示す。図 2 1は、 3次元構造を持つ 2つの組織モデルの一例を示す図である。

図 2 2は、内包モデルにおける各推定結果を示す第 1の図である。

図 2 3は、内包モデルにおける各推定結果を示す第 2の図である。図 2 4は、層状モデルにおける各推定結果を示す第 1の図である。

図 2 5は、層状モデルにおける各推定結果を示す第 2の図である。

図 2 6は、超音波診断システムの基本構成を示す図である。

図 2 7は、超音波診断システムで用いた超音波スキャナの具体的構成を示す図である。発明を実施するための最良の形態

以下、添付図面に従って本発明に係る超音波診断システムの好ましい実施の形態について説明する。本実施例の超音波診断システムでは、包絡線相関計算の際、複合自己相関法で 1次元の相関窓で 1次元探索していた処理を 2次元の相関窓を用いて 2次元探索することにより横方向の移動にも対応した拡張複合自己相関法と呼ばれる方法を採用している。また、この拡張複合自己相関法は、軸方向には 2分の 1波長間隔、横方向にはライン間隔の格子点でのみ包絡線相関計算を行うことにより計算量を減少させて高速化を図っている。ただし、複合自己相関と同様に拡張複合自己相関法でも位相情報を利用して軸方向の変位推定精度を向上させている。しかし、横方向変位の推定はキャリアとなる信号がないため位相情報は利用できない。そのため、横方向変位推定精度は空間相関法と同様に横方向サンプリング間隔（超音波ビームライン間隔）により制限されてしまう。しかし、後で提案する弾性係数分布再構成法では軸方向の歪み（変位）分布のみから弾性係数分布を推定できるため、この実施例では横方向変位推定精度の向上は特に行わない。

本実施例の拡張複合自己相関法の具体的な構成を説明する前に、本発明の前提技術となる複合自己相関法について、図 1〜図 6を参照して説明する。図 1は、超音波診断装置の原理を説明するための図である。図から明かなように、超音波探触子である超音波プローブ 1 0は電気信号を超音波に、また超音波を電気信号に変換するものであり、この超音波プローブ 1 0を用いて生体組織 1 1内に超音波パルスを放射する。生体組織 1 1内に放射された超音波パルスは音響インピーダンスの異なる第 1の境界 1 2で一部が反射され、反射エコー 1 2 aとして超音波プローブ 1 0側に向かい、その残りは透過していく。透過した超音波パルスは次の音響インピーダンスの異なる第 2の境界 1 3で同様に一部が反射され、反射ェコ一 1 3 aとして超音波プローブ 1 0側に向かい、その残りは透過する。このようにして反射した超音波の反射エコー信号は超音波プローブ 1 0によって受波され電気信号の反射エコー信号に変換される。このとき、超音波プローブ 1 0から超音波パルスが放射されてから距離 Lの位置にある反射物体 1 4 (音響インピ一ダンスの異なる境界）からの反射エコー信号を受信するまでの時間 tは、次式 ( 1 ) となる。

1L

c … ( 1 ) ここで、 cは生体内での音速であり、軟組織では 1 5 0 0 [m/秒] にほぼ一定とみなせる。よって、反射エコー信号を受信するまでの時間 tを計測すればプローブから反射物体までの距離 Lを求めることができる。この反射エコー信号は、生体組織の音響特性情報を有することから、その反射エコー信号に基づいて Bモ ―ド断層像などの生体組織の情報をディスプレイに表示することができる。例えば、超音波診断装置を用いて、組織の硬さに関する弾性特性を計測することが行われている。弾性特性を計測する原理は、組織に機械的振動を与えてその横波の伝搬速度から硬さ情報を評価するものであり、硬い組織では横波の伝搬速度が速ぐ軟らかい組織では横波の伝搬速度が遅いことを基にしている。ただし、厳密には生体組織中を伝わる横波の伝搬速度 Vは、次式（2 ) に示したように、組織の密度 P、せん断弾性係数 ^、せん断粘性係数 ₂、および振動の角周波数 ωに関係している。

V

一方、図 2に示すように、組織を静的に圧縮してその際生じる組織内の歪み分布から弾性特性を評価する方式が提案されている。これは、硬い組織では歪みが小さく、軟らかい組織では歪みが大きくなることに基づいている。図 2 (A)は、静的圧縮による組織弾性計測方式の具体例を示す図である。図 2 (B) は、静的圧縮による組織弾性計測方式の原理を示す図である。図から明かなように、この方式は、従来からある超音波診断装置および超音波プローブ 10をそのまま用いて組織 1 1の圧縮前の超音波エコー信号（圧縮前 R F信号）を計測する。その後、超音波プローブ 10により組織 1 1をわずかに（数パーセント程度）圧縮し、組織 1 1の圧縮後の超音波エコー信号（圧縮後 RF信号）を計測する。そして、計測された組織圧縮前後の R F信号から圧縮によつて組織内部の各点がどれだけ動いたかという移動量である変位分布を推定する。この変位分布推定手法の主なものとしては、空間相関を用いる方法とドプラの原理を用いる方法がある。

図 3は、空間相関法の原理を示す図である。この方法は、圧縮によって生じた組織内部の変位分布を組織圧縮前後の RF信号（または RF信号の包絡線）から 2次元相関関数を用いたテンプレートマッチングにより推定する手法である。その具体的な処理は以下のようになる。まず、組織圧縮前後の RF信号（またはその包絡線信号）を（t， x), i ₂ (t, x) とすると、この 2つの信号の相互相関係数 C ( t， X ； n， m) は、次式（3) となる。

ここで、 tは超音波ビーム方向（軸方向）の座標、 Xはそれに直交する方向（横方向）の座標、 t。は軸方向の相関窓サイズ、 x_Qは横方向の相関窓サイズ、 Liは軸方向のサンプリング間隔、 L_x は横方向のサンプリング間隔、 n， mは整数である。そして、この相互相関関数が最大となるときの（n, m) を（k， 1 ) とすると、計測点（t， X) における軸方向の変位 u_y と横方向の変位 u_x はそれぞれ次式のようにして求められる。

u_v =kL₍ _x = \ L_X

ただし、横方向のサンプリング間隔 L_x は軸方向のサンプリング間隔よりも大きいので、推定される変位成分の精度は横方向成分の方が軸方向成分よりも悪くなる。上記の処理を各計測点について行い変位分布推定する。空間相関法は 2 次元の変位ベクトル成分を推定できるという特徴がある。また、組織が大変形した場合（5%程度）でも変位分布を推定できる。しかし、計算量が膨大になるため超音波計測の利点であるリアルタイム性を損なってしまう。また、変位推定精度もサンプリング間隔により制限されてしまうため、後に述べるドプラ法と比べると精度が悪いという問題点もある。

図 4は、ドプラ法の原理を示す図である。この方法は、圧縮によって生じた組織内部の変位分布を組織圧縮前後の R F信号から血流計測に用いられているドプラの原理を利用して推定する手法である。その具体的な処理は以下のようになる。まず、組織圧縮前後の RF信号を次式（4) のようにモデル化する。

¾(/) = ReU( e-^j(fi,oi"^e)J

i₂(t) = ReU(i -て) _e -ゾ。('- ] I

… (4) ここで、（t) は圧縮前の RF信号、 i₂ (t) は圧縮後の RF信号、 A (t) は包絡線、 ω。は超音波の中心角周波数、ては時間シフトである。そして、この 2つの RF信号をそれぞれ直交検波すると、次式のようなベースバンド信号が得られる。

sM) = A(t-T)e^j{a,oT+e)

² … (5) そして、この 2つの信号の相関関数 R₁₂ (t) は次式で表される。

… （6)

ここで、 R_A (t) は包絡線の自己相関関数、 t。は超音波ビーム軸方向の相関窓サイズである。また、 *は複素共役を表している。よって、この相関関数 R₁₂ ( t ) の位相 ψ ( t ) から圧縮による時間シフトて、軸方向変位 U_y は次式（7) のようにして求まる。て =

CT

u_v =——

y 2 … (7) ただし、 cは組織内の音速であり、生体内で一定と仮定する。

上記の処理を各計測点について行い変位分布推定する手法がドプラ法であり、ドプラの原理を基にした血流計測と同じ処理となっている。そのため、リアル夕ィム計測が可能であるという利点がある。また、位相情報を用いているので変位推定精度が空間相関法よりも良い。しかし、組織内部の移動量が大きくなると、例えば超音波中心周波数の 4分の 1波長以上となると、エイリアシングを起こしてしまい正しい変位推定ができないという問題点がある。また、上式からもわかるように軸方向の変位成分のみしか推定できない。

そこで、それら 2つの手法の長所を組み合わせた「複合自己相関法（CA法： Combined Autocorrelation Method)」を本願の発明者等は提案している。図 5 は、本願発明者等が提案している複合自己相関法の原理を示す図である。複合自己相関法は、ドプラ法におけるエイリアシングの問題を RF信号の包絡線による相関を用いることによって解決したものである。その具体的な処理は以下のようになる。

まず、組織圧縮前後の RF信号をドプラ法のときと同じように次式のようにモデル化する。

=

L ( = Re 0―て) eーゾ小

… （8) ここで、（t) は圧縮前の RF信号、 i₂ (t) は圧縮後の RF信号、 A (t) は包絡線、 ω。は超音波の中心角周波数、ては時間シフトである。そして、この 2つの RF信号をそれぞれ直交検波すると、次式のようなベ一スバンド信号が得られる。

s_l t) = A t)e^je

そして、この 2つの信号間の複素相関関数 R₁₂ ( t ； n) を次式のように定義する。

("=…，一 2, - 1，0，1,2，...） ... _{( 1}0) ここで、 Tは超音波の周期、 R_A ( t ； r) は包絡線の自己相関関数、 t。は相関窓サイズである。また、 *は複素共役を表している。ここで、 nは探索時における計測点の移動ピッチ数であり、 n=0の場合は、ドプラ法における自己相関関数の式（6) に一致する。すなわち、 n = 0の場合はドプラ法と同じであり、軸方向変位が超音波の波長の 4分の 1以上になるとエイリアシングを起こしてしまう。そこで、この問題を克服するために次式（1 1) で定義される包絡線相関係数 C ( t ； n) を用いる。

ただし、 R_u (t ； 0) は、 (t) の自己相関関数、 R₂₂ ( t ； n) は s₂ (t + ΠΤ../2) の自己相関関数である。そして、この包絡線相関係数 C (t ； n) が最大となる nの値を kとすると、その n = kのときの R₁₂ ( t ； k) の位相 φ ( t ； k) は、エイリアシングの起きていない位相となる。これは、包絡線相関を計算する間隔を 2分の 1波長（周期）に選んだためである。ちなみに、この 2 分の 1波長はエイリアシングを起こさないための最大の間隔である。よって、この φ (t ； k) を用いることにより組織圧縮による時間シフト及び軸方向変位 u_v は次式のように求まる。て =— - ^J→-k—

a>_n 2

cて

u„ =

' ² … （12) ただし、 cは組織内の音速であり、生体内で一定と仮定する。

上記の処理を各計測点について行い変位分布推定する手法が複合自己相関法であり、ドプラ法を拡張したような手法となっている。そのため、リアルタイム計測が可能な手法となっている。また、包絡線相関を用いることによってドプラ法では計測不可能であった大変形の場合（超音波の 4分の 1波長以上の変位が生じる場合）の変位分布推定にも対応している。

図 6は、複合自己相関法の基本アルゴリズムを実行する回路構成を示すプロック図である。図において、加圧前直交検波回路（QD) 131は、加圧前のェコ —信号 X ( t ) を入力し、それぞれ直交検波して、直交検波信号 I x ( t ), Qx ( t ) 信号を、第 1相関演算回路 133及び第 1相関係数演算回路 1350〜 1 35 Nに出力する。第 1加圧後直交検波回路（QD) 1320は、加圧後のェコ一信号 y (t) を入力し、それぞれ直交検波して直交検波信号 Y (t) = I y + j Qy (I y (t), Qy (t)) を、第 1相関係数演算回路 1340及び第 2相関演算回路 1350に出力する。第 1遅延回路 134は、ェコ一信号 y (t) を超音波の周期 Tだけ遅延させ、遅延したェコ一信号 =y (t -T) を第 2加圧後直交検波回路（QD) 1321に出力する。第 2遅延回路 135は、第 1遅延回路 134によって遅延されたエコー信号 _yi =y ( t -T) を同じく超音波の周期 Tだけ遅延させ、遅延したエコー信号 y₂ =y (t - 2T) を次段の第 2 加圧後直交検波回路（QD) 1322 (図示せず）に出力する。以後、 N段の遅延回路を用いて順次周期 Tの整数倍だけ信号を遅延して、遅延した信号を加圧後直交検波回路に供給する。

第 1相関演算回路 133は、信号 I x， Qxに基づいて相関値 R X Xを演算し、それを各第 2相関係数演算回路 1380〜138Nに出力する。第 2相関演算回路 1340は、加圧後直交検波回路 1320からの直交検波信号 I y (t), Qy (t) を入力し、信号 I y， Qyに基づいて相関値 Ry yを演算し、それを第 2 相関係数演算回路 1380に出力する。第 1相関係数演算回路 1350は、加圧前直交検波回路 131からの直交検波信号 I X ( t)， Qx ( t ) 及び第 1加圧後直交検波回路 1320からの直交検波信号 I y (t), Qy (t) を入力し、複素ベースバント信号 S_R , S, を求め、それを第 3相関演算回路 1360及び位相差演算回路 1370に出力する。第 3相関演算回路 1360は、第 1相関係数演算回路 1350からの複素ベースバント信号 S_R， S, を入力し、それに基づいて相関値 I Rxy Iを求め、それを第 2相関係数演算回路 1380に出力する。位相差演算回路 1370は、第 1相関係数演算回路 1350からの複素べ一スバント信号 S_R , を入力し、それに基づいて位相差 Φ。（t) を求める。第 2相関係数演算回路 1380は、第 1相関演算回路 133からの相関値 R xx、第 3 相関演算回路 1360からの相関値 I Rxy I、並びに第 2相関演算回路 1 34 0からの相関値 Ry yを入力し、これらの各相関値に基づいて相関係数 C。（t) を演算し、出力する。

第 2加圧後直交検波回路（QD) 1321は、第 1遅延回路 134によって遅延されたェコ一信号 y, =y (t -T) を入力し、それぞれ直交検波して直交検波信号 (t) = I y₁ + j Qy₁ (I y, (t), Q_yi (t)) を、第 1相関係数演算回路 1341及び第 2相関演算回路 1351に出力する。第 2相関演算回路 1341は、第 2加圧後直交検波回路（QD) 1321からの直交検波信号 I y, (t), Qy, ( t) を入力し、その信号 I y, (t), Qy, (t) に基づいて相関値 Ry iを演算し、それを第 2相関係数演算回路 1381に出力する。第 1相関係数演算回路 1351は、加圧前直交検波回路 131からの直交検波信号 I x (t)， Qx (t)、第 2加圧後直交検波回路（QD) 1 321からの直交検波信号 I y_t ( t), Qy, (t) を入力し、複素べ一スバント信号 S_K1， S_n を求め、それを第 3相関演算回路 1361及び位相差演算回路 1371に出力する。第 3相関演算回路 136 1は、第 1相関係数演算回路 1351からの複素べースバント信号 S_R1, S„を入力し、それに基づいて相関値 I Rxyi Iを求め、それを第 2相関係数演算回路 1381に出力する。位相差演算回路 1371は、 9731

15 第 1相関係数演算回路 1351からの複素べ一スバント信号 S_R1, S„を入力し、それに基づいて位相差 (t) を求める。第 2相関係数演算回路 1381は、第 1相関演算回路 133からの相関値 R xx、第 3相関演算回路 1361からの相関値 | Rxy, I、並びに第 2相関演算回路 1341からの相関値 Ry_lYlを入力し、これらの各相関値に基づいて相関係数 (t) を演算し、出力する。以下同様に、第 1遅延回路 135以降の第 2加圧後直交検波回路（QD) 13 22〜132N、第 2相関演算回路 1342〜134N、第 1相関係数演算回路 1352〜135N、第 3相関演算回路 1362〜136N、位相差演算回路 1 372〜 137 N及び第 2相関係数演算回路 1382〜 138 Nは、上述の 1段目及び 2段目の回路群と同様の処理を実行し、相関係数 C₂ (t) 〜C_N (t) 及び位相 Φ₂ (t) 〜φ_Ν (t) を出力する。このように、図 6の複合自己相関法の基本アルゴリズムを実行する回路は、加圧後のエコー信号 y ( t) を遅延回路 134〜13Nで超音波の周期 TZ2 (2分の 1波長）だけ遅延し、それを直交検波回路（QD) 1320〜132 Nを用いて個別に直交検波している。

前述のように組織庄縮に伴う変位分布が推定されたら、それを空間微分することにより歪み分布が得られる。歪み分布は定性的に組織の弾性特性を表しているものであり、歪み分布からでもかなりの弾性特性に基づいた診断は行える。しかし、肝硬変などの病変部全体が硬くなるような場合には、定量的な弾性係数によつて評価しなければ組織診断は難しい。そのため、近年、組織弾性係数分布再構成法についても研究されるようになってきた。しかし、今のところスタンダードな手法はなく、いずれの手法も研究段階であるというのが実状である。

一方、組織弾性係数分布は先にも述べたように組織内部の歪み分布と応力分布から求められる。しかし、応力分布を直接計測することは現状では困難であるため、歪み分布と組織圧縮の際の境界条件から逆問題的に弾性係数分布を再構成することになる。そのため、一般的に逆問題を解くことは難しく、現在提案されて ^る弾性係数再構成法も数少ない。従来から提案されている弾性係数再構成法を以下に説明する。

第 1に、 1次元を仮定した方法（1次元弾性体を仮定）がある。これは、 1次元弾性体を仮定して歪みの逆数を弾性係数とみなす方法である。この方法は弾性係数再構成法ではなく、歪みの逆数を求めるだけであるので、歪みにおける非定量性をそのまま残している。

第 2に、弾性方程式から応力項を消去する方法（等方性弾性体、非圧縮性、平面歪み状態を仮定）がある。これは、平面歪み状態を仮定した場合の弾性方程式を変形し、応力項を消去した方程式を用いて組織圧縮の際の境界条件（体表での外部圧力分布、または体表での変位）と歪み分布（せん断歪み成分を含む歪みテンソルの全成分) から組織弾性係数分布を再構成する手法である。ただし、絶対的な弾性係数を推定するには、弾性係数が前もってわかっている領域（参照領域）が必要となる。

第 3に、弾性微分方程式を積分する方法（等方性弾性体、非圧縮性、平面応力状態を仮定）がある。これは、平面応力状態を仮定した場合の弾性方程式を変形した応力項を含まない弾性係数に関する微分方程式を体表付近での弾性係数を基準として順次積分していくことにより、歪み分布（せん断歪み成分を含む歪みテンソルの全成分）から組織弹性係数分布を再構成する方法である。そのため、体表付近の弾性係数分布が前もって分かっている領域が必要であり、また体表付近を基準として積分を行っていくので奥に行くほど誤差が積算されるという問題点もある。

第 4に、摂動法を用いた手法（等方性弾性体、近非圧縮性、平面歪み状態を仮定）がある。これは、平面歪み状態を仮定した場合の弾性方程式を基にした摂動法により体表での外部圧力分布と超音波ビーム方向（軸方向）の歪み分布とから反復的に組織弾性係数分布を再構成する方法である。

以上、本発明の原理の基本事項及び本発明の具体的な解決課題を説明した。以下、これらの課題を解決できる本発明の実施例について説明する。図 7は、本発明の一実施例である超音波診断システムの概略構成を示すブロック図である。図において、超音波プローブ 9 1は、被検体内へ超音波を送波すると共にその反射波を受波するものであり、従来のセクタスキャンプローブ（セクタフェイズドアレイプローブ）リニアスキャンプローブ（リニアアレイプローブ）又はコンペックススキャンプローブ（コンベックスアレイプローブ）などである。

超音波プローブ 9 1からは、組織圧縮前後の R F信号が直交検波器 9 2に出力される。直交検波器 9 2は、組織圧縮前後の R F信号をそれぞれ組織圧縮前後の複素包絡線信号（I Q信号）に変換し、複素 2次元相関計算部 9 3に出力する。複素 2次元相関計算部 9 3は、組織圧縮前後の R F信号間における 2次元相関を計算し、その相関が最大となる位置を横方向変位計算部 9 4及び軸方向変位計算部 9 5に出力し、そのときの相関関数の位相を軸方向変位計算部 9 5に出力する。ただし、軸方向にはエイリアシングを起こさずに位相を検出できる最大の間隔である超音波中心周波数の 2分の 1波長間隔でのみ相関を計算するものとする。これは、超音波診断システムのリアルタイム表示を優先させるためである。従って、高精度な相関を計算するためには、この 2分の 1波長間隔に限定する必要はない。横方向変位計算部 9 4は、複素 2次元相関計算部 9 3からの横方向の相関最大位置に基づいて横方向の変位 ιι_χ を計算し、それを横方向歪み計算部 9 6に出力する。一方、軸方向変位計算部 9 5は、複素 2次元相関計算部 9 3からの軸方向の相関最大位置及びそのときの位相に基づいて軸方向の変位 u _y を計算し、それを軸方向歪み計算部 9 7に出力する。横方向歪み計算部 9 6は、横方向変位計算部 9 4からの横方向変位 u_x の分布を空間的に微分することにより横方向歪み分布 ε _χ を計算し、それを量子化部 9 8に出力する。一方、軸方向歪み計算部 9 7 は、軸方向変位計算部 9 5からの横方向変位 u_y の分布を空間的に微分することにより軸方向歪み分布 ε _γ を計算し、それを量子化部 9 8に出力する。量子化部 9 8は、横方向歪み分布 ε _χ 及び軸方向歪み分布 ε _γ をグレースケール表示（又はカラー表示）するために各歪み分布を量子化し、表示部 9 9に出力する。表示部 9 9は、量子化された各歪み分布を表示する。

次に、図 7の超音波診断システムで採用した拡張複合自己相関法の動作について説明する。まず、組織圧縮が極僅か（数パーセント以下）である場合、組織内部を局所的に見れば平行移動したと見なすことができ、組織圧縮前後の R F信号を次式のようにモデル化できる。 ¾(/,JC) = Re - 。'- j

/,ひ， x)

ここで、（t , X) は圧縮前の RF信号、 i ₂ ( t， x) は圧縮後の RF信号、 A (t， x) は包絡線、 ω。は超音波の中心角周波数、ては時間シフトであり軸方向変位の時間概念、 u_x は横方向変位、 0は初期位相である。また、ここではドプラ法や複合自己相関法のときと違い横方向の変位も考慮して圧縮前後の RF信号をモデル化している。そして、この式の中で最終的に知りたいパラメ一夕は、軸方向の変位 u_y =c て/ 2 (すなわち、時間シフトて）と横方向変位 u_x である。ただし、 cは組織内の音速であり、生体内で一定と仮定する。

そこで、これらの組織圧縮前後の RF信号を、まず直交検波器 92でそれぞれ直交検波する。すなわち、各 RF信号に超音波の中心周波数と同じ周波数の s i n波， c o s波をかけ、それぞれ低域通過フィルタをかけると、次式（14) のような複素ベースバンド信号 , s₂ が得られる。

s₁(t,x) = A(t,x)e^Jff s it, x = A(t-T,x- u_x)e^ji<OoT+0)

… (14) そして、この s_t ( t , x) と s₂ ( t + Ώ. Ύ/2 , x +mL) との間の 2次元複素相関関数 R₁₂ (t, x ； n, m) を次式（1 5) のように定義する。

*

R (t_fx n.m) +v_ix+mL + w) dvaw

T - + )

= R_A(t,x;t -n—,u_x -mL)e ² ノ

ム

("=-N_rain,"',-2,-l,0，l，2,"'，N_ma

(m = ^,-,-2,-1,0X2,-^)

··■ (1 5 ここで、 Tは超音波の周期、 Lは横方向サンプリング間隔（ライン間隔)、 R_A ( t , X ; て， U_x ) は包絡線の自己相関関数、 t。は軸方向相関窓サイズ、 X。は横方向相関窓サイズである。また、レはて時間軸方向の積分のずらし量、 Wはビームライン方向の積分のずらし量であり、 *は複素共役を表している。そして、この

2次元複素相関関数を用いて、次式（1 6)に示す 2次元包絡線相関係数 C ( t, X； n, m) を定義する。

ただし、 R„ ( t， x； 0， 0) は s , ( t， x) の自己相関関数、 R₂₂ ( t , x ； n, m) は s₂ ( t +nT/2, x+mL) の自己相関関数である。そして、この包絡線相関係数を用いて複合自己相関法の場合と同様にエイリアシングの問題を克服する。すなわち、各計測点（t , X) における C ( t , X ； n, m) と R i2 ( t , χ ； η, m) の位相 φ ( t , χ； n, m) との組 { C ( t , χ； η, m) , φ ( t , χ ; η, m)} をすベての ηと mについて求める。ここで、 nと mの範囲が十分広ければ、すなわち、包絡線相関を行う探索範囲が十分に大きければ、包絡線相関係数が最大となる（n， m) = (k, 1) に対する位相 φ (t, x； k, 1) はエイリアシングの起きていない位相となる。これは、包絡線相関 C (t, X； n， m) が最大となる（n, m) = (k, 1) のとき、 s _l ( t , x) と s ₂ ( t + kT/2, x+ 1 L) との時間シフ卜の大きさ I r-kT/2 |が TZ2 よりも小さくなる、すなわち、 ( t , ； k, 1 ) | =ω。 I て _kTZ2 Iが πよりも小さくなるためである。よって、このエイリアシングの起きていない Φ ( t , X ； k, 1 ) を用いれば、計測点（t， X) における正確な時間シフトで、軸方向変位 u_y 、そして横方向変位 ιι_χ が次式（1 7) のように求まる。て

·'·(1 7) u_x = IL

ただし、 cは組織内での音速（ここでは軟組織における一般的な音速 1 50 Om Zsで一定とする) である。したがって、組織内のすべての点で上記のように軸方向変位と横方向変位を計算すれば、軸方向変位分布 u_y (x, y) と横方向変位分布 u_x (x, y) が得られる。

また、各変位分布を次式のように空間微分することにより、軸方向歪み分布 £_y (X, y) と横方向歪み分布 ε_χ (X, y) が次式（18) のように求められる。 du (x,y)

O, ）

dy

_du_x{x,y

d x

(18) 以上のようにして、組織圧縮前後の RF信号から軸方向と横方向の変位（歪み）分布を推定することができる。ただし、上式の u_x = 1 Lからもわかるように横方向変位の精度は横方向サンプリング間隔（ライン間隔）によって制限されてしまうため、精度は若干劣るということはあるが、リアルタイムに観察できるので実用性の高いものである。

また、上述の拡張複合自己相関法は、組織の横方向移動に対応するように 2次元の相関窓と探索範囲を用いて、組織圧縮の際の組織と超音波プローブの相対的な横方向移動には対応しているが、組織圧縮の際に軸方向と横方向にそれぞれ垂直な方向（2次元超音波走査面に垂直な方向（スライス方向））の変位が生じてしまい組織移動が起こった場合には、 2次元の拡張複合自己相関法では歪みの推定を行うことができない。そのため、より安定したシステムにするには、上述の拡張複合自己相関法を 3次元の相関窓と 3次元の探索範囲を用いることにより簡単に拡張することが可能である。

図 9及び図 10は、 3次元複合自己相関法の基本アルゴリズムを示すフローチャ一ト図である。なお、図 10は、図 9の処理の一部の詳細を示すフローチヤ一卜図である。

ステップ S 11では、ステップ S 23の判定処理と組み合わせて、第 1番目の走査線から第 L番目の走査線についてそれぞれ同様の処理を行うために、走査線番号レジスタ 1に「1」を格納する。ステップ S 12では、ステップ S 18の判定処理と組み合わせて、厚み方向（フレーム）を一 Uから Uまで順次シフトする処理を実行する。

ステップ S 13では、ステップ S 17の判定処理と組み合わせて、方位方向（走査線）を一 Vから Vまで順次シフトする処理を実行する。

ステップ S 14では、ステップ S 16の判定処理と組み合わせて、距離方向（軸方向）を 0から Mまで順次シフトする処理を実行する。

ステップ S 1 5では、複合自己相関法により、距離方向（軸方向）における包絡線の相関係数 C ( 1, t ； u, V, n) を計算する。この複合自己相関法は、従来の方法でやってもいいが、それだと計算に時間を要するので、ここでは、高速化された複合自己相関法を用いて相関係数 C (1， t ； u, V, n) の計算を行う。この高速化複合自己相関法については後述する。

ステップ S 16では、前のステップ S 14と組み合わせられた処理であり、距離方向レジスタ nがその最大値 Mに達したか否かの判定を行い、達した場合にはステップ S 17に進み、そうでない場合はステップ S 14にリターンし、距離方向レジスタ nをインクリメント処理する。

ステップ S 17では、前のステップ S 13と組み合わせられた処理であり、方位方向レジス夕 Vがその最大値 Vに達したか否かの判定を行い、達した場合にはステップ S 18に進み、そうでない場合はステップ S 13にリタ一ンし、方位方向レジスタ Vをインクリメント処理する。

ステップ S 18では、前のステップ S 12と組み合わせられた処理であり、厚み方向レジスタ uがその最大値 Uに達したか否かの判定を行い、達した場合にはステップ S 1 9に進み、そうでない場合はステップ S 12にリターンし、厚み方向レジスタ uをインクリメント処理する。

ステップ S 19では、ステップ S 12〜ステップ S 18の処理によって求められた相関係数 C ( 1， t ； u, V, n) (u = -U, ...0, U) (v--V, ... 0， .„, V) (n= 0, 1， ...， N) の中から最大となる（u, v, n) を求め、それを（u。， ν。， n。）とする。なお、本実施例では、圧縮変位のみを考えて n=0, 1， ...， Nとしたが、弛緩変位をも求める場合は、 n =— Μ,···, 0, 1, „.， Nのように設定すればよいことは当業者に容易に理解されるところである。 • ステップ S 20では、ステップ S 19で求められた相関係数 C (1, t ； u₀, v。， n_Q ) について、その位相差 Φ (1， t ； u。， v_Q ， n。）を計算する。ステップ S 21では、最終的な位相差として、 φ = η。 ττ + φ ( 1 , t ； u₀, v。， n。）を計算する。

ステップ S 22では、（u。， v。， n。）の近傍の相関係数 C ( 1， t ； u, v， n) を用いて、方位方向の変位： v = v_D +Δν及び厚み方向の変位： u = u_Q +Διιを計算する。

ステップ S 23では、前のステップ S 1 1と組み合わせられた処理であり、走査線番号レジスタ 1が Lに達したか否かの判定を行い、達した場合にはステップ S 24に進み、そうでない場合はステップ S 11にリターンし、走査線番号レジスタ 1をインクリメント処理する。

ステップ S 24では、組織圧縮に伴う変位分布が推定されたら、それを空間微分することにより歪み分布を計算する。

図 10は、図 9のステップ S 15の高速化された複合己相関法の詳細を示すフローチャート図である。

ステップ S 31では、圧縮前の RF信号 Xと圧縮後の RF信号 yをそれぞれ直交検波して、以下のように I， Q信号を求める。

( t) →I x、 Qx (X ( t) = I x+ j Qxとする）

y ( t) →I y、 Qy (Y ( ) = I y+ j Qyとする）

ステップ S 32では、相関： Rxy、 Rxx、 R y yを次式に基づいて計算する。

Rxy= S X ( t +ソ） · Y* ( t + v) d v

Rx x= S X ( t + v) · X* (t +レ） dリ .

Ry y= S Y ( t + y) · Y* ( t + v) d v

ステップ S 33では、求められた相関 Rxy、 Rxx、 Ryyを用いて相関係数 C_Q を次式に基づいて計算する。

C„ = I Rxy I /7~Rxx^Ry y 2003/009731

23 ステップ S 34では、変数 nに 1をセットする。

ステップ S 3 5では、 Y_n (t) =Y (t -nT) e^iwnTを計算する。

ステップ S 36では、次式に基づいて Rxy_n ， Ry_ny_nを計算する。

Rxy_n = S X ( t + v) · Y_n* ( t + ν) d ν

= X ( t +ソ） · Υ* ( t -nT+ v) β¹'^ωηΊά v

Ry_ny_n= Y„ (t + v) Y_n* (t +ソ） dリ

= S Y ( t -nT+ v) - Y* ( t -nT+ v) d v

ステップ S 37では、求められた Rxy_n , Ry_ny_nを用いて相関係数 C_nを次式に基づいて計算する。

C„ = I Rx y_n I / Rx xT Ry_ny_s

ステップ S 38では、変数 nをインクリメント処理する。

ステップ S 39では、変数 nが最大値 Mに達したか否かを判定し、達した場合はこの処理を終了し、達していない場合は、ステップ S 3 5にリタ一ンし、同様の処理を繰り返す。

図 10のフローチャートでは、 Rxy_n ， Ry_ny_nを求めるのに、ステップ S 35で Y_nを Yから導いている。このために、回路構成の簡略化を図ることができる。以下、どのようにして Υ_η を Υから導くかについて説明する。

まず、加圧前のエコー信号 X (t) を

X ( t) =u ( t ) c o s (ω t + θ)

軸方向に加圧後のエコー信号 y (t) を

y ( t ) =x ( t + r) =u (t +て） c o s (ω (t +て） +0) とする。各信号 x， y， y_n の直交検波値は、

X (t) _→I x=0. 5 u ( t ) c o s e

Qx = - 0. 5 u ( t ) s i n Θ

(X ( t) = I x + j Q x = 0. 5 u (t) e~^ie)

y ( t ) → I y = 0. 5 u ( t + r ) c o s (ω τ + Θ)

Q y—— 0. 5 u (t +て ) s i n (ω τ + Θ)

(Y ( t) = I y + j Qy= 0. 5 u ( t +t) e-^{j(tt,r + 0)} ) y_n (t) =y (t -nT)

=u (t +t— nT) c o s (ω ( t + r -nT) + Θ) =u (t +t— nT) c o s (c t +ω (て— nT) +Θ) となる。ここで、 Tは半周期なので、

I y_n= 0. 5 u ( t + r -nT) c o s (ω (て— nT) +θ)

Qy„=- 0. 5 u (t +て _ηΤ) s i n (ω (て— nT) + θ)

(Υ = I y_n+ j Qy_n=0. 5 u ( t +て — nT) _e_" "て- )）

となる。以上の式から以下のような関係が成り立つ。

Υ_η (t) = I y_n+ j Qy_n '

=0. 5 u ( t + r - nT) e -^j(cu(r-^{nT) + e)}

=Y (t -nT) e^iMnT

これから Y_n ( t) は Y ( t) = I y+ j Qyから求まることになる。

従って、 Rxy_n, Ry_ny_nも、次式のように X、 Yから求めることができる。 Rxy_n =4 S X (t + v) · Y_n * ( t + ν) d ν

= 4 S X (t +リ） · Υ* (1; _ηΤ +ソ） e^jwnTd v

I x y_n I =R u_n

= u ( t +リ） u ( t + r -nT+ v) d v

= 4nx (t + v) - Y_n * ( t +ソ） I d v

=4 Π X ( t + y) · Y* ( t -nT+ y) e^ia,nT I d v

=4 S I X ( t + v) - Y* ( t -nT+ v) I d v Ry„y_n= u ( t + r -nT+ v) u (t +て一 nT+リ） d v

=4 Π Y„ (t + v) · Y_n * (t + ν) I d v

=4 S Y ( t -nT+ ν) · Y* (t一 nT +レ） dソ

ここで、 *は複素共役を表している。

図 1 1は、この 3次元複合自己相関法の基本アルゴリズムを実行する回路構成を示すブロック図である。複合自己相関法を実行する回路構成を従来技術の図 7 に示すようなものにすると、直交検波回路 1 320〜132 Nが多段接続されることによって、直交検波回路 1320〜1 32 Nの処理に時間を要するようになり、その処理時間は膨大なものとなってしまい、高速な演算処理の妨げとなり、リアルタイムな画像表示の妨げとなっていた。そこで、前述のような基本アルゴリズムに応じた図 11のような回路構成を採用することによって、複合自己相関法を実行する回路の処理速度を高速化している。

加圧前直交検波回路（QD) 131は、加圧前のエコー信号 X (t) を入力し、それぞれ直交検波して、直交検波信号 I x (t)， Qx (t) 信号を、第 1相関演算回路 133及び第 1相関係数演算回路 1350〜135 Nに出力する。加圧後直交検波回路（QD) 132は、加圧後のエコー信号 y (t) を入力し、それぞれ直交検波して直交検波信号 Y (t) = I y + j Qy (I y (t), Qy ( t)) を、第 1相関係数演算回路 1350、第 2相関演算回路 1340及び第 1遅延回路 134及び第 2遅延回路 135に出力する。第 1遅延回路 134及び第 2遅延 '回路 135は、直交検波信号 Y (t) をそれぞれ超音波の周期 Tだけ遅延させ、遅延した直交検波信号 Y (t -T) を第 1相関係数演算回路 1351、第 3遅延回路 136及び第 4遅延回路 137に出力する。第 3遅延回路 136及び第 4遅延回路 137は、直交検波信号 Y (t -T) をそれぞれ超音波の周期 Tだけ遅延させ、遅延した直交検波信号 Y (t -2T) を次段の第 1相関係数演算回路及び遅延回路（図示せず）に出力する。以後、 N段の遅延回路を用いて順次周期丁の整数倍だけ信号を遅延して、遅延した信号を第 1相関係数演算回路に供給する。第 1相関演算回路 133は、信号 I X, Qxに基づいて相関値 R X Xを演算し、それを各第 2相関係数演算回路 1380〜138 Nに出力する。第 2相関演算回路 1340は、加圧後直交検波回路 132からの直交検波信号 I y( t)，Qy ( t) を入力し、信号 I y， Qyに基づいて相関値 Ry yを演算し、それを第 2相関係数演算回路 1380に出力する。第 1相関係数演算回路 1350は、加圧前直交検波回路 131からの直交検波信号 I x ( t ) , Qx (t) 及び加圧後直交検波回路 132からの直交検波信号 I y (t), Qy (t) を入力し、複素ベースパント信号 S_R ， S, を求め、それを第 3相関演算回路 1360及び位相差演算回路 1 370に出力する。第 3相関演算回路 1360は、第 1相関係数演算回路 135 0からの複素べ一スバント信号 S_R， S, を入力し、それに基づいて相関値 I R xy Iを求め、それを第 2相関係数演算回路 1380に出力する。位相差演算回路 1370は、第 1相関係数演算回路 1350からの複素ベースバント信号 S_R， S, を入力し、それに基づいて位相差 Φ。 (t) を求める。第 2相関係数演算回路 1380は、第 1相関演算回路 133からの相関値 Rxx、第 3相関演算回路 1360からの相関値 I Rxy し並びに第 2相関演算回路 1340からの相関値 Ryyを入力し、これらの各相関値に基づいて相関係数 C。（t) を演算し、出力する。

第 2相関演算回路 1341は、第 1遅延回路 134及び第 2遅延回路 135からの遅延後の直交検波信号 I y (t— T), Qy (t -T)を入力し、信号 I y (t 一 T)， Qy (t -T) に基づいて相関値 Ry!y,を演算し、それを第 2相関係数演算回路 1381に出力する。第 1相関係数演算回路 1351は、加圧前直交検波回路 131からの直交検波信号 I X (t), Qx (t)、第 1遅延回路 134及び第 2遅延回路 135からの遅延後の直交検波信号 I y (t一 T)， Qy (t一 T) を入力し、複素べ一スバント信号 S_Rl, S„を求め、それを第 3相関演算回路 13 61及び位相差演算回路 1371に出力する。第 3相関演算回路 1361は、第 1相関係数演算回路 1351からの複素ベースバント信号 S_R1， S„を入力し、それに基づいて相関値 I Rxy, Iを求め、それを第 2相関係数演算回路 1381 に出力する。位相差演算回路 1371は、第 1相関係数演算回路 1351からの複素ベースバント信号 S_R1, S„を入力し、それに基づいて位相差 (t) を求める。第 2相関係数演算回路 1381は、第 1相関演算回路 133からの相関値 Rxx、第 3相関演算回路 1361からの相関値 I Rxy, I、並びに第 2相関演算回路 1341からの相関値 Ry iを入力し、これらの各相関値に基づいて相関係数 C, (t) を演算し、出力する。

第 3遅延回路 135及び第 4遅延回路 136から次段の第 2相関演算回路 13 42〜134N、第 1相関係数演算回路 1352〜135N、第 3相関演算回路 1362〜 1 36 N、位相差演算回路 1372〜； 137 N及び第 2相関係数演算回路 1382〜138Nは、上述と同様の処理を順次遅延された遅延後の直交検波信号 I y (t— 2T) · · · I y ( t -NT), Qy ( t - 2 T) · · · Qy (t 一 N T) に対して実行し、相関係数 C ₂ ( t ) 〜 C_N ( t ) 及び位相 φ ₂ ( 1；) 〜 Φ_Ν ( t ) を出力する。

次に、 3次元有限要素モデルを用いた弾性係数分布再構成法について説明する。弾性係数分布再構成逆問題を簡単化するため、この実施の形態では組織をモデル化する。これはまた、今回提案する弾性係数分布再構成法において有限要素法を用いるためでもある。この実施の形態では、組織を以下のように仮定及びモデル化する。

まず、組織を等方性弾性体と仮定する。組織歪み分布を推定する際、外部から圧力を加えて組織を静的に圧縮する。しかし、組織圧縮前後の R F信号間の相関を保っために、微小圧縮しか行わない。そのため、この場合、応力と歪みの間の関係はほぼ線形である。すなわち、組織を弾性体として近似できる。また、今回対象としている軟組織はほぼ等方性が成り立つため、この実施の形態では組織を等方性弾性体と仮定する。

さらに、組織を近非圧縮性と仮定する。生体組織は、非圧縮性（ポアソン比レ = 0 . 5 ) に近いことが知られている。そこで、ポアソン比を 0 . 4 9とし、生体内で一定とする。ここで、完全な非圧縮性を仮定しないのは、ポアソン比レ = 0 . 5とすると特殊な弹性方程式となり、今回の提案手法で用いている有限要素法が適用できなくなるためである。そして、ポアソン比を生体内で一定とすることで弾性係数分布推定の推定パラメ一夕をヤング率のみとすることができ、逆問題を簡単化できる。また、ポアソン比はヤング率に比べ生体中であまり変化しないパラメ一夕であるため、この実施の形態ではポアソン比を 0 . 4 9で一定とする。組織を 3次元有限要素モデル化する。この弾性係数分布再構成法では有限要素法を用いるため、組織を有限個の直方体要素に分割する。そして、各要素内では、弾性係数、応力、歪みは一様であると仮定する。一般的に逆問題を解くには、それに対応する順問題を理解することが重要である。今回の歪み分布と境界条件から弾性係数分布を推定する逆問題の場合、それに対応する順問題とは、弾性係数分布と境界条件から歪み分布を求めることである。そして、この順問題の数値解法の 1つが有限要素法（FEM： F ini t e El ement Me thod) である。ここで、有限要素法とは対象となる連続体を有限個の要素の集合で近似し、この集合体に対して成り立つ連立 1次方程式を数値的に解く手法のことである。なお、有限要素法の定式化については後述する。ここでは有限要素法とは「入力として物体の弹性係数分布と境界条件を与えれば、出力としてそのときの歪み（変位）分布と応力分布が得られるもの」として捉えておけば十分である。

この実施の形態では、組織を等方性弾性体で近似するため、組織内では以下のような弾性方程式（つりあい方程式 ·歪み一変位関係式 ·応力一歪み関係式）が成り立つ。

つりあい方程式は次式（19) のように表される。 =0

dx,

^(/，ゾ ⁼¹，²'³⁾ — （19) ここで、 i.、 j =1, 2, 3は、 x、 y、 zに対応する。また、歪み—変位関係式は次式（20) のように表される。 r du_t ouj

応力一歪み関係式（一般化したフックの法則）は次式（21) のように表される。び

… (21) 上記の式ではテンソル表現を用いており、実際にはつりあい方程式として 3つの式、歪み一変位方程式として 6つの式、応力—歪み関係式として 6つの式が存在する。また、座標系（x_t , x₂， x₃ ) は（x， y， z) を表しており、その他の記号は以下のことを表している。

E ：ヤング率 (弾性係数とはヤング率のことを表している)

V ：ポアソン比 T/JP2003/009731

29 ε_η：歪みテンソル

(ε_ηη= ε"+ ε₂₂ + ε ₃₃：体積歪み）

s_n：応力テンソル

djj：クロネッカーのデル夕

U; ：変位ベクトル '

f i ：体積力ベクトル（重力の影響は無視できるため、ここでは =0とする）ここで、応力一歪み関係式をについて整理すると、次のような歪み一応力関係式（22) が得られる。

ただし、 _Snn=_Sll + _S22 + s₃₃である。よって、この式の中で i = j =2とし、ヤング率 Eについて整理すると次式（23) が得られる。

E び22 - ^V (び 11 +び 33)

"22

び - び x+び J

sy - (23) 従って、軸方向（この実施の形態では、 y方向を超音波ビーム方向、すなわち軸方向とする）の歪み成分と全方向の応力成分がわかれば、ヤング率すなわち弾性係数を求めることができる。なお、上述の計算式からは、応力分布を直接計測することは現状では困難である。そこで、この実施の形態では応力分布と弾性係数分布を交互に推定 ·更新しながら、推定弹性係数分布を実際の分布に近づけていく。その弾性係数分布再構成の手順は、以下のようになる。

第 1に、未知弾性係数分布の初期値分布 {E °} として一様分布を考える。第 2に、初期弾性係数分布 {E'^Q} のときに生じる応力分布 {s' } を 3次元有限要素法により求める。具体的には、まず組織モデル内の各要素に対して歪み一変位関係式及び応力一歪み関係式をつりあい方程式に代入して得られる次式のようなつりあい方程式（24) を作る。

ただし、 u_n _ du_x du 0u₃

dx_n dx_l dx₂ dx^

- (25) である。そして、この連立方程式を以下のような境界条件のもとガウスの消去法を用いて変位について解き、弹性係数分布 {E'°} のときの変位分布 {u'°} を求める。

u I = 0 び i \y=top= Pi

σ I , =

n \x,z=si e · · · ( 2 Γ、上式において、 p；は体表における外部圧力ベクトル、 s_nは側面に垂直な方向の応力成分である。また、上段の式は底面が固定されていることを示し、中段の式は体表での応力分布は外部圧力分布に等しいことを示し、下段の式は側面が拘束されていないことをそれぞれ示している。次に、この変位分布 {iT^fl} を歪み一変位関係式に代入して、弾性係数分布 {E'^Q} のときの歪み分布 { ε '°} を求める。そして、この歪み分布 { ε '°} を応力一歪み関係式に代入することにより、弾性係数分布 {Ε } のときの応力分布 {_s °} を得る。

第 3に、 3次元有限要素法により得られた応力分布と拡張複合自己相関法により推定した軸方向（y方向）歪み分布 { ε_γ } を用いて、弾性係数分布 {E'^k} を次式（27) により更新する。， — ）

sy - (27) ただし、この式は、上述の応力一歪み関係式を ε_υについて整理し、式中の i = j =2とし、ヤング率 Eについて整理した式を書き改めたものであり、式中の k は繰り返し回数を表している。

第 4に、上述のように更新された弾性係数分布と上述の境界条件を用いて再び 3次元有限要素解析を行い、応力分布を更新する。

そして、第 3及び第 4の処理を繰り返すことにより弾性係数分布を実際の分布に近づけていく。ただし、次式（2 8 ) の条件が満たされた時点で弾性係数分布推定は収束したとみなし、推定を終了する。 it+1

丄 < Γ

N Y、 ¹

… ( 2 8 ) ここで、 1は要素番号、 Νは要素数、 Γはしきい値である。

以上が、今回提案した 3次元有限要素モデルによる弾性係数分布再構成法であり、この方法は 3次元のつりあい方程式を基に弾性係数分布を推定している。そのため、本手法は従来の手法よりもより実際の生体組織に近い仮定に基づいているので、より正確な弾性係数推定が能になる。また、本手法は精度良く推定可能な軸方向の歪み分布のみから弾性係数分布を再構成するため、安定した弾性係数分布再構成が行える。ただし、本手法は組織弾性係数の 3次元分布を推定する手法であるため、 2次元アレイ超音波プローブを用いるか、 1次元アレイ超音波プローブをスライス方向に機械的に走査することにより、対象を 3次元的に走査する必要がある。

この実施の形態の拡張複合自己相関法と 3次元有限要素モデルによる弾性係数分布再構成法の有効性をシミュレーションによって実証する。図 1 2は、このシミュレーションの手順の概略を示す図である。

第 1に、推定したい弾性係数分布を持つ組織モデルを作成する。このとき、組織モデル内には超音波エコー信号を発生させるための散乱体を分布させておく。第 2に、この組織モデルに対して外部圧力を加え、計算機上で組織圧縮を行う。そして、この圧縮による各散乱点の移動先を有限要素法などにより求める。第 3 に、組織モデル圧縮前後の散乱体分布を基に圧縮前後の R F信号を作成する。第 4に、この圧縮前後のシミュレ一ション R F信号に対して拡張複合自己相関法を適用し、組織歪み分布を推定する。第 5に、拡張複合自己相関法により推定された歪み分布と組織モデル圧縮の際に設定した境界条件（外部圧力分布など）とから 3次元有限要素モデルによる弾性係数分布再構成法により組織弾性係数分布を推定する。

今回のシミュレ一シヨンで用いた組織モデルの弾性係数分布は、各シミュレ一シヨンにおいて異なるが、いずれの場合も等方性弾性体を仮定する。なお、各シミュレーションで設定した弾性係数の値としては、今回の組織弾性計測システムで主な'対象としている乳房組織の弾性係数にほぼ即している。また、組織圧縮前後のシミュレ一ション RF信号を作成するために、各組織モデルには点散乱体を分布させた。その際、点散乱体の平均密度としては 500個/ cm3とし、組織圧縮前の散乱体の位置は一様乱数により、散乱係数は平均 1. 0、標準偏差

3の正規乱数により決めた。そして、組織圧縮後の散乱体位置は有限要素解析の結果に応じて組織圧縮前の各散乱体を移動させることにより決めている。ここで、実際の組織の散乱体に関する情報は未知であるが、シミュレーション RF信号を基に Bモード像にした際、実際の組織の Bモード像に近くなるように各パラメ一夕を設定する。

この実施の形態では、組織モデルに対する組織圧縮前後のシミュレーション R F信号を次式（29) のように組織圧縮前後の散乱体関数と超音波システムの点広がり関数との畳み込みにより作成する。

ix, y,z = fj J/z( -x y- ,ζ-ζ'^ (ΛΓ' ,γ,ζ' )dx' dy' dz

iつ (x,

ここで、 i, (x, y, z) は組織圧縮前の RF信号、 i₂ (x, y, z) は組織圧縮後の RF信号、 h ( , y, z) は超音波システムの点広がり関数（インパルス応答)、 t _L (x, y, z) は組織圧縮前の散乱体関数、 t₂ (X, y, z) は組織圧縮後の散乱体関数である。ただし、散乱体関数とは組織モデル内の散乱体が存在する位置ではその散乱係数の値をとり、その他の位置では 0であるような関数である。また、組織圧縮後の散乱体関数 ₂ (X, y, z) は組織圧縮前散乱体関数 (x, y， z) の各散乱体の位置を組織モデルの変形に応じて移動させたものである。ただし、組織圧縮に伴う各散乱体位置での変位は有限要素解析により得られる要素節点での変位べクトルを線形補間することにより求めている。

また、この実施の形態ではシミュレーション超音波システムとして無焦点、かつ減衰のないシステムを仮定する。すなわち、超音波システムの点広がり関数 h (x, y， z) は空間的に不変であると仮定する。さらに、点広がり関数は次式 (30) のように方向ごとに分離できると仮定する。 h(x,y,z) = h_x(x)h_y(y)h_z(z)

… (30) ここで、 hy (y) は超音波ビーム方向の点広がり関数、 hx (x)， h z (z) はそれぞれ超音波ビームに直交した方向の点広がり関数である。ただし、 X方向は超音波断層面内の方向（横方向）、 z方向は超音波断層面に垂直な方向（スライス方向）とする。そして、各方向の点広がり関数は実際の超音波装置により計測したワイヤー 'ターゲット（水中に張った直径 0. 13mmのワイヤ一）からの反射エコー分布を基に作成する。図 13は、超音波中心周波数 5. 0 MHzの場合に用いた各点広がり関数の一例を示す図である。図 13 (A) は軸方向の点広がり関数 hy (y) を示し、これはガウス関数に正弦波をかけたものによって実際のワイヤー .夕一ゲットからの反射エコー分布を近似し、図 13 (B) は横方向点の広がり関数 hx (x) を、図 13 (C) はスライス方向の広がり関数 h z (x) をそれぞれ示し、これらはガウス関数によって実際のワイヤー ·ターゲットからの反射ェコ一分布を近似する。また、各関数のパラメ一夕は中心周波数に応じて変えており、各シミュレ一ションの際に改めて説明する。

次に、変位（歪み）分布推定法として今回提案した拡張複合自己相関法の有効性をシミュレーションにより評価する。まず、拡張複合自己相関法の複合自己相関法に对する拡張点である組織の横方向変位に対応できる点について検証する。図 14は、組織モデルの概略を示す図である。組織モデルは、外形 60mmx 60mm (2次元）で、一様な弾性係数分布をもつモデルである。そして、この組織モデルを軸方向に一様な 3%の歪みが生じるように圧縮する。ここで、このシミュレーションに関しては拡張複合自己相関法のみの評価を行うため、組織モデルとしては単純な 1次元弾性体を仮定している。そして、組織の横方向移動（超音波プローブに対する相対的な横方向移動）に関する影響を検証するため、軸方向の圧縮と同時に横方向に 0.0mmから 1.4mmまでの横方向変位を与えた。ただし、横方向に関しては単純な平行移動であり、組織に対して超音波プローブが完全に滑つた場合を想定している。

そして、この組織モデルに対して変形前後のシミュレーション RF信号を作成する。その際用いた超音波システムに関するパラメータは、中心周波数 5. 0M Hz、パルス幅 0. 5mm、超音波ビーム幅 1. 0mm、走査ライン間隔 0. 5 mm、サンプリング周波数 30MHzである。そして、この圧縮前後のシミュレ —ション RF信号を用いて、今回提案した拡張複合自己相関法により歪み分布を推定する。また、比較のために複合自己相関法と空間相関法を用いても歪み分布推定を行った。ここで、単純に精度等を比較できるように各手法では同じサイズの相関窓と探索範囲を用いた。具体的には、拡張複合自己相関法と空間相関法では 1. 6 mm (軸方向） χ2. 5 mm (横方向）の 2次元相関窓と 5. 6 mm (軸方向） χ7. 5 mm (横方向）の 2次元探索範囲を用い、 1次元処理の複合自己相関法では軸方向だけ同じ 1. 6 mmの 1次元相関窓と 5. 6 mmの 1次元探索範囲を用いた。

このようにして各手法により歪み分布を推定した結果、図 15〜図 17のようになる。ここで、図 15は横方向変位に対する各手法の歪み推定誤差を示している。 ◊は複合自己相関法、口は拡張複合自己相関法、 △は空間相関法を示す。ただし、歪み推定誤差 eとしては次式（31) を用いた。

ここで、 ε は推定された歪み、 _{£ i}は実際の歪み（理想値）、 iは要素番号、 N は要素数である。また、図 16は横方向変位が 0. Ommの場合の各手法（複合自己相関法、拡張複合自己相関法、空間相関法）によ,り推定した歪み分布、図 1 7は横方向変位が 0 . 4 mmの場合の各手法（複合自己相関法、拡張複合自己相関法、空間相関法）により推定した歪み分布を示す図である。なお、図 1 6及び図 1 7は軸方向の深さごとに推定された歪みの平均値と標準偏差を表している。このシミュレーション結果より、従来の複合自己相関法（C A法）では組織の相対的な横方向変位が超音波ビームを越えて生じてしまう（今回の場合、横方向変位がビーム幅の片側分である 0 . 5 mmを超えてしまう）と歪み推定が急に悪くなつてしまうのに対し、拡張複合自己相関法では横方向変位の大きさにかかわらず安定して歪み推定が可能であることが理解できる。また、空間相関法も横方向変位にかかわらず安定して歪み推定が行えるものの、拡張複合自己相関法と比ベると 2倍以上誤差が大きく精度が悪いことが理解できる。また、各手法の処理時間を比較したところ、下表のように拡張複合自己相関法は複合自己相関法に比ベて 3 . 6倍時間がかかってしまうものの、空間相関法と比べると 1 Z ( 7 . 7 ) の時間しかかからなかった。そのため、拡張複合自己相関法はある程度リアル夕ィム性が保たれていることが理解できる。

次に、斜め方向圧縮に関する検証について説明する。前述のシミュレーションでは簡単な 2次元の均一組織モデルを用いたが、次は実際の生体組織と同じ 3 次元構造を持つ組織モデルを用いてシミュレーションを行う。また、超音波プローブにより組織を圧縮する際、超音波ビーム方向（軸方向）に圧縮するのが理想であるが、ここでは斜めに圧縮してしまった場合の影響を検証する。

図 1 8は、斜め方向圧縮の検証に使用される組織モデルの概略を示す図である。組織モデルは、図 1 8 (A) に示すように、外形が 6 O mmx 6 O mmx 6 O mm の 3次元構造をしており、この組織モデル中に直径 1 5 mm、長さ 6 O mmの硬い円柱形内包物が存在するようなモデルである。ここで、周辺の弹性係数（ヤング率）は 1 0 k P a、内包物の弾性係数は周辺より 3倍硬い 3 0 k P aとする。ただし、この弾性係数の値は今回主な対象としている乳房組織の弾性係数および乳がんの弾性係数を基に設定する。そして、この組織モデルを 2通りの方法で圧縮を行った。 1つ目は、図 1 8 (B ) に示すように、この組織モデルに対して上面から軸方向に一様な 2 0 0 P aの外部圧力を加えて、組織モデルを軸方向に約 2 %圧縮する。 2つ目は、図 1 8 ( C ) に示すように、この組織モデルに対してモデル上面から斜め方向の一様な外部圧力（軸方向に 2 O O P a、横方向に 3 0 P a ) を加えて、組織モデルを斜め方向に圧縮する。

そして、上記の 2通りの場合についてそれぞれ組織圧縮前後のシミュレーション R F信号を作成し、拡張複合自己相関法により歪み分布を推定する。なお、ここでも比較のために複合自己相関法と空間相関法による歪み分布推定も行う。ただし、単純に比較できるように各手法で用いた相関窓サイズと探索範囲は同じとし、そのサイズは前のシミュレ一ションと同じとする。また、シミュレーション R F信号作成の際に用いた超音波システムに関するパラメ一夕も前のシミュレ一ションと同じ、中心周波数 5 . 0 MH z、パルス幅 0 . 5 mm、横方向超音波ビ —ム幅 1 . 0 mm、スライス方向超音波ビーム幅 2 . 0 mm、走査ライン間隔 0 . 5 mm、サンプリング周波数 3 0 MH zとする。

上記のようにして行ったシミュレーションの結果は、図 1 9及び図 2 0に示すようになる。ここで、図 1 9は組織モデルを単純に軸方向に圧縮した場合の歪み分布推定結果を示し、図 2 0は組織モデルを斜め方向に圧縮した場合の歪み分布推定結果を示す。なお、各図における理想的な歪み分布とは、各条件で 3次元有限要素解析を行って得られた軸方向歪み分布であり、この歪み分布を正解としている。また、図 1 9及び図 2 0における結果は組織モデルの中央断面での結果である。ここで、図 2 0において理想歪み分布が左右対称でないのは斜め方向に圧縮した影響であり、特に今回の場合は図に向かって右斜め下に圧縮したため、図右上の部分の横方向変位が大きくなつている。

そして、このシミュレーションの結果より、軸方向に圧縮した場合は拡張複合自己相関法と複合自己相関法はほぼ同じ結果となったが、斜め方向に圧縮した場合は横方向変位が大きくなるため、複合自己相関法では推定できなくなる領域が生じてしまうのに対し、拡張複合自己相関法では先のシミュレーション同様に横方向変位の大きさにかかわらず安定して歪み推定が可能であることが改めて確認された。また、空間相関法も横方向変位の大きさには依存しないものの拡張複合自己相関法の結果と比較すると明らかに推定精度が悪いのが見てとれる。そのため、前のシミュレーションと合わせて改めて拡張複合自己相関の有効性が確認された。

前述の拡張複合自己相関法は、精度良く、かつ高速に組織歪み分布を推定できることができる。そこで、次に組織弾性映像システムの第 2段階である歪み分布から弾性係数分布を推定する手法で、今回提案した 3次元有限要素モデルによる弾性係数分布再構成法についてシミュレーシヨンにより検証を行う。

今回提案した弾性係数分布再構成法の最大の特徴は 3次元の力学的つりあい方程式に基づいて弾性係数分布を推定することである。そこで、今回提案した手法と手順的には同じであるが、 2次元の力学的つりあい方程式を基にしている点だけが異なる 2次元再構成法を用いて、提案した 3次元再構成法と比較することにより検証を行う。この 2次元再構成法では、 2次元の平面歪み状態を仮定した弾そこで、まず組織モデルとして、実際の生体組織と同じ 3次元構造を持つ図 2 1 のような 2つのモデルを用いる。図 2 1は、 3次元構造を持つ 2つの組織モデルの一例を示す図である。図 2 1 ( a ) の内包物モデルは、乳がんを模擬したモデルで、外形 1 0 0 mmx 1 0 0 mmx 1 0 0 mmのモデル中に直径 2 0 mmの硬い内包物が存在するもので、内包物の弾性係数は 3 0 k P a、周辺の弾性係数は 1 O k P aとする。これらの弾性係数の値は、先ほどのシミュレーションと同じように実際の乳房組織の弾性係数を基に決めている。また、周辺と内包物のポアソン比としては、共に非圧縮性に近いため 0 . 4 9で一様とする。そして、図 2 1 ( b ) の層状モデルは筋肉などの層状のものを模擬したモデルで、外形 1 0 0 m mx 1 0 O mmx 1 0 O mmのモデル中に厚さ 2 0 mmの硬い層がモデル中に存在するというもので、硬い層の弾性係数は 3 0 k P a、周辺の弾性係数は 1 0 k P aとする。そして、このモデルの場合もポアソン比は 049で一様とする。そして、図 21 (a) の内包物モデルの場合はモデル上,部から軸方向に 100 P aの一様な外部圧力により、図 21 (b) の層状モデルの場合はモデル上部から軸方向に 150 P aの一様な外部圧力により、各モデルをコンピュータ上で圧縮する。ここで、 2つのモデルで外部圧力の強さを変えているのは、各モデルとも同じ約 1 %の歪みが生じるようにするためである。そして、これらの組織モデル圧縮前後のシミュレーション R F信号を作成し、拡張複合自己相関法により軸方向歪み分布を推定する。そして、この推定された軸方向歪み分布と圧縮の際の境界条件とから、提案した 3次元弾性係数分布再構成法により弾性係数分布を推定する。また、同じ軸方向歪み分布と境界条件を用いて比較のために 2次元再構成法によっても弾性係数分布を推定する。ここで、シミュレーション RF信号を作成するために用いた超音波システムのパラメータとしては、中心周波数 3. 7 5MHz, パルス幅 0. 75mm、横方向超音波ビーム幅 2. 0mm, スライス方向超音波ビーム幅 2. 0mm、走査ライン間隔 2. 0mmである。また、拡張複合自己相関法における相関窓のサイズは 3. 2mm (軸方向） _x4. 0mm (横方向）、探索範囲は 1 1. 2mm (軸方向） χ14. 0mm (横方向）とする。さらに、 3次元有限要素モデルよる弾性係数分布再構成では 2 mm (軸方向） χ 2 mm (横方向） x5mm (スライス方向）の直方体要素 50000個を用いて組織モデルを構成する。

そして、このシミュレーションの結果は、図 22〜図 25に示すようになる。図 22及び図 23は内包モデルにおける各推定結果を示す、図 24及び図 25は層状モデルにおける各推定結果を示す。ただし、 3次元再構成法では弾性係数の 3次元分布を推定しているが、ここではモデル中央断面での結果のみを示す。これは、 2次元再構成法では 2次元断面での弾性係数分布しか推定できないので、 'ここでは比較できる中央断面のみ示す。また、各組織モデルにおける 3次元再構成結果と 2次元再構成結果を数値的に評価したところ次のような結果が得られた。周辺領域におけ周辺領域におモデル中心に

弹' I生 1糸数 s¾ るコノ卜

[%] [%] ラスト誤差

厂し/ <y ¾ Ί _! 内包 3次元再構成法 3. 5 15. 5 11. 0 物モ

デル 2次元再構成法 30. 9 17. 9 35. 9 層状 3次元再構成法 8. 5 26. 8 3. 1 モデ

ル 2次元再構成法 24. 9 22. 1 43. 5 ここで、用いた評価用のパラメ一夕は、周辺領域における弾性係数誤差 e_s 、周辺領域における標準偏差 SD_S 、モデル中心におけるコントラスト誤差 e_e であり、それぞれ次式のように定義する。

Ε,—ΕΛ

¾ =■

Ε.

… (32) ただし、上式における S は周辺領域における総和、 E は推定された弾性係数、

Eは実際の弾性係数、 N_s は周辺領域の要素数、 E" は周辺領域における弾性係数の平均値、 E"_eはモデル中心における推定弾性係数、 E_e はモデル中心における実際の弾性係数、 E_s は周辺領域における実際の弾性係数である。

以上のシミュレ一ション結果より、今回提案した 3次元有限要素モデルによる弾性係数分布再構成法の方が平面応力状態を仮定した 2次元再構成法よりもより正確な弾性係数を推定できることが確認された。ここで、 3次元再構成法では弹性係数の値が正確に推定されているのに対し、 2次元再構成法では実際の弾性係数の値よりも小さく推定しまっている。これは、前もって予想されていたように 2次元状態では 2次元考察面に垂直な方向の動きが制限されてしまうためである。そのため、実際の生体組織と同じ 3次元を基にした弹性係数分布再構成法が必要 2003/009731

40 であることがこのシミュレーションにより改めて示された。

次に、上述の拡張複合自己相関法と 3次元有限要素モデルによる弾性係数再構成法を実装した超音波診断システムの具体的構成について説明する。図 2 6は、この超音波診断システムの基本構成を示す図である。この超音波診断システムは、 3次元超音波スキャナ 2 8 1、超音波診断装置 2 8 2、パーソナルコンピュータ 2 8 3、パルスモ一タコントローラ 2 8 4、パルスモータ 2 8 5、圧力計 2 8 6 などから構成される。 3次元超音波スキャナ 2 8 1は超音波パルスを組織内に送信し、かつ組織からの超音波エコー信号を受信するためのものである。ただし、ここでは 3次元有限要素モデルによる弹性係数再構成法を用いるため、組織内の 3次元的なデ一夕が必要になる。そこで、この超音波診断システムでは 3次元超音波スキャナ 2 8 1は 3次元走査が可能な構成となっている。超音波診断装置 2 8 2は 3次元超音波スキャナ 2 8 1を制御したり、リアルタイムに超音波 Bモ一ド画像を表示して計測部位を決定したりするためのもので、従来の超音波診断装置をそのまま用いることができる。この超音波診断装置はフルデジタル化された装置で内部にフレームメモリを持っているため、計測した R F信号を一時的に保存しておくことが可能となっている。パーソナルコンピュータ 2 8 3は、超音波診断装置 2 8 2によって計測された R F信号を受け取り、組織弾性特性推定のための処理（前述の提案手法の処理）、および処理結果の表示を行うためのものである。パルスモータ 2 8 5は組織圧縮を制御するためのものであり、位置固定が可能なアームの先端に固定されており、かつパルスモータ 2 8 5の可動部分には 3 次元超音波スキャナ 2 8 1が取り付けてある。そして、パルスモータコントローラ 2 8 4によりこのパルスモータ 2 8 5を制御し、超音波スキャナ 2 8 1を組織表面で上下に動かすことにより数パ一セントの微小組織圧縮を精度良く行う。圧力計 2 8 6は弾性係数再構成の際に必要な境界条件である体表上での圧力を測るためのもので、超音波スキャナ 2 8 1と体表との間に置かれる。ただし、ここでは超音波スキャナ 2 8 1で組織圧縮を行った際の体表上での圧力分布は一様であると仮定し、圧力計 2 8 6で計測された値をその圧力値として用いる。

図 2 7は、この超音波診断システムで用いた超音波スキャナ 2 8 1の具体的構成を示す図である。 3次元超音波スキャナ 2 8 1は、超音波振動子が 2次元平面状に並んだ 2次元アレイ型ではなく、コンベックス型の 2次元走査プロ一ブを水中で機械的にスライス方向に振ることにより 3次元走査を行うタイプのものである。

図 2 6の超音波診断システムは、乳がん診断を主な対象としているため、超音波スキャナの特性も乳腺用に設定してある。今回用いた超音波スキャナの主な特性としては、超音波中心周波数 7 . 5 MH z , サンプリング周波数 3 0 MH z、走査ライン数 7 1本、フレーム数 4 4枚、振動子の振れ角 3 0 °、 1回の 3次元走査にかかる時間 0 . 5秒となっている。ここで、振動子の振れ角とはコンペックス型のプロ一ブをスライス方向に振る際の振れ角のことであり、フレーム数とはコンベックス型のプローブを 1回振る際に取得される走査面（フレーム）の数である。また、水中のワイヤ一·ターゲットにより超音波パルスの特性を計測したところ、パルス幅約 0 . 5 mm、横方向ビーム幅約 1 . 5 mm、スライス方向ビーム幅約 2 . 6 mmであった。

図 2 6の超音波診断システムの弾性計測の動作例について説明する。まず、超音波診断装置 2 8 2のリアルタイム Bモード像を見ながら、アームにつながった 3次元超音波スキャナ 2 8 1を動かし、計測を行いたい部位に超音波スキャナ 2 8 1を設定する。この際、超音波スキャナ 2 8 1は 3次元走査を行わず（すなわち、コンベックス型のプローブを機械的に振らず）、スキャナの中央面の Bモード像のみを超音波診断装置 2 8 2に表示する。これは、今回の超音波診断装置 2 8 2で 3次元走査をするとリアルタイムに Bモードを表示できないためである。従つて、リアルタイムに Bモ一ドを表示することのできる超音波診断装置の場合は、 3次元走査を行いながら部位の設定を行なうことができる。計測部位に超音波スキヤナ 2 8 1を移動させたら、アームの可動部を固定し、超音波スキャナ 2 8 1 を固定する。そして、組織圧縮前の 3次元 R F信号を計測する。これは、 3次元走査用のポタンを押すことにより自動的に 3次元走査される。また、 1回の 3次元走査にかかる時間はわずか 0 . 5秒である。このとき、計測された圧縮前の R Fデ一夕は超音波装置内のフレームメモリに保存しておく。次に、パルスモータ · コントローラ 2 8 4の圧縮用のポタンを 1回押すことにより、アームに取り付けられたパルスモータ 2 8 5が前もって設定しておいた量だけ超音波スキャナ 2 8 1を動かし、超音波スキャナ 2 8 1自身により組織を圧縮する。そして、パルスモー夕 2 8 5が止まって組織圧縮を行っている状態で、再び 3次元走査用のポ夕ンを押し、組織圧縮後の R Fデ一夕を計測する。ここで、組織圧縮後の R Fデ一タも圧縮前の R Fデータと同様に超音波装置 2 8 2内のフレームメモリに保存される。また、圧縮している状態で超音波スキャナ 2 8 1の端に取り付けられた圧力計の圧力を計測しておく。以上で、計測部は終わりで、組織圧縮を解除し、被験者は解放される。

被検者解放後は、パーソナルコンピュータ 2 8 3から超音波診断装置 2 8 2内のフレームメモリにアクセスし、組織圧縮前後の R Fデータをパーソナルコンビユー夕 2 8 3上のハードディスクに保存する。これは、超音波装置内のフレ一ムメモリは一時的なものであり、 1回の計測分しか容量がないためである。パーソナルコンピュータ 2 8 3は、拡張複合自己相関法と 3次元有限要素モデルによる弾性係数分布再構成法のプログラムを実行し、組織圧縮前後の R Fデータから歪み分布及び弹性係数分布を推定する。そして、パーソナルコンピュータ 2 8 3は、表示用のプログラムによってモニタ上に Bモード像、歪み像、弾性係数像をこの発明の超音波診断システム、歪み分布表示方法及び弾性係数分布表示方法によれば、横方向変位に対応して変位分布を推定することができると共に超音波ビーム方向（軸方向）の歪み分布のみから弾性係数分布を再構成することができるという効果がある。

なお、本件は、包絡線信号を例に説明したが、波形の特徴量を表すものであれば、振幅、波高値、波数を含む波形の相関を特定可能なパラメータでもよい。

Claims

請求の範囲

1 . 被検体との間で超音波信号を送受信する超音波探触子と、

この超音波探触子により検出された信号の特徴量を蓄積する記憶手段と、この記憶手段に蓄積された前記被検体の圧縮前後の前記特徴量に基づき前記圧縮前後の前記特徴量の相関係数及び前記圧縮前後の前記受信信号間の位相差を求める相関演算手段と、

この相関演算手段によって求められた前記相関係数及び前記位相差に基づいて前記圧縮に伴う前記計測点の変位並びに前記被検体の組織の歪み分布を求める演算手段と、

前記歪み分布を表示する表示手段とを備えたことを特徴とする超音波診断シスアム。

2 . 前記特徴量は、包絡線、振幅、波高値、波数を含む波形の相関を特定可能なパラメータであることを特徴とする請求項 1に記載の超音波診断システム。

3 . 前記相関演算手段は、前記記憶手段に蓄積された前記被検体の圧縮前後の前記超音波探触子により検出された信号の包絡線信号からなる超音波ビームデータに計測点を設定し、該計測点を少なくとも超音波ビーム方向に移動して、前記計測点における前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるとともに、前記圧縮前後の前記受信信号間の位相差を求め、

前記演算手段は、前記相関演算手段により求められた前記相関係数が最大となる位置及び前記位相差に基づいて前記圧縮に伴う前記計測点の変位を求める変位演算手段を含んでなることを特徴とする請求項 1に記載の超音波診断システム。

4 . 前記演算手段は、前記計測点における変位を空間微分することによつて前記被検体の組織の歪み分布を求める歪み演算手段を含んでなることを特徴とする請求項 3に記載の超音波診断システム。

5 . 前記相関演算手段は、前記計測点を前記超音波ビーム方向に前記超音波信号の 2分の 1波長間隔で移動して前記相関係数の最大位置を求めることを特徴とする請求項 3に記載の超音波診断システム。

6 . 前記相関演算手段は、前記計測点における前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるに際し、圧縮後の前記包絡線信号の自己相関関数を求め、この自己相関関数を前記計測点の移動に合わせて前記超音波信号の 2分の 1波長間隔でシフトして前記相関係数を求めることを特徴とする請求項 3に記載の超音波診断システム。

7 . 前記被検体を有限個の要素に分割して少なくとも 2次元の有限要素モデルを作成し、そのモデル化の情報と前記歪み分布を用いて弾性係数分布を演算する弾性係数演算手段を備え、前記表示手段に前記弾性分布を表示することを特徴とする請求項 6に記載の超音波診断システム。

8 . 前記相関演算手段は、前記記憶手段に蓄積されたスライス面に対応する前記被検体の圧縮前後の前記包絡線信号のフレームデータに計測点を設定し、該計測点を前記フレームデ一夕に対して少なくとも 2次元方向に移動して、前記計測点を囲む 2次元相関窓に属する前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるとともに、前記圧縮前後の前記 R F信号間の位相差を求め、

前記演算手段は、前記相関演算手段によって求められた前記相関係数が最大となる位置及び前記位相差に基づいて前記圧縮に伴う前記計測点の少なくとも 2次元方向の変位を求める変位演算手段を含んでなることを特徴とする請求項 1に記載の超音波診断システム。

9 . 前記 2次元方向は、前記超音波探触子により受信される超音波ビームの方向と、該超音波ビームの走査方向であることを特徴とする請求項 8に記載の超音波診断システム。

1 0 . 前記相関演算手段は、前記計測点を前記超音波ビーム方向に前記超音波信号の 2分の 1波長間隔で移動し、前記超音波ビームの走査方向に前記超音波ビームピッチで移動して前記相関係数の最大位置を求めることを特徴とする請求項 9に記載の超音波診断システム。

1 1 . 前記被検体を有限個の要素に分割して少なくとも 2次元の有限要素モデルを作成し、そのモデル化の情報と前記歪み分布を用いて弾性係数分布を演算する弹性係数演算手段を備え、前記表示手段に前記弾性分布を表示することを特徴とする請求項 8に記載の超音波診断システム。

1 2 . 前記相関演算手段は、前記計測点を囲む 2次元相関窓に属する前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるに際し、圧縮後の前記包絡線信号の自己相関関数を求め、この自己相関関数を前記計測点の移動に合わせてずらして前記相関係数を求めることを特徴とする請求項 8に記載の超音波診断システム。

1 3 . 前記記憶手段に蓄積されたフレームデータは、複数のスライス面のフレームデ一夕を有するボリュームデ一夕であり、

前記相関演算手段は、前記計測点を前記ポリュ一ムデータに対して 3次元方向に移動して、前記計測点を囲む 3次元相関窓に属する前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるとともに、前記圧縮前後の前記 R F信号間の位相差を求めることを特徴とする請求項 5に記載の超音波診断システム。

1 4 . 前記 3次元方向は、前記超音波探触子により受信される超音波ビームの方向と、該超音波ビームの走査方向と、これらに直交するスライス方向であることを特徴とする請求項 1 3に記載の超音波診断システム。

1 5 . 前記相関演算手段は、前記計測点を前記超音波ビーム方向に前記超音波信号の 2分の 1波長間隔で移動し、前記超音波ビームの走査方向に前記超音波ビームピッチで移動し、前記スライス方向に前記超音波ビームのスライスピッチで移動して前記相関係数の最大位置を求めることを特徴とする請求項 1 4に記載の超音波診断システム。

1 6 . 前記相関演算手段は、前記圧縮前後の前記 R F信号間の前記超音波ビーム方向、前記超音波ビームの走査方向及びこれらに直交するスライス方向の位相差を求めることを特徴とする請求項 1 3に記載の超音波診断システム。

1 7 . 前記被検体を有限個の要素に分割して少なくとも 3次元の有限要素モデルを作成し、そのモデル化の情報と前記歪み分布を用いて弾性係数分布を演算する弾性係数演算手段を備え、前記表示手段に前記弾性分布を表示することを特徴とする請求項 1 3に記載の超音波診断システム。

1 8 . 前記相関演算手段は、前記計測点を囲む 3次元相関窓に属する前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるに際し、圧縮後の前記包絡線信号の自己相関関数を求め、この自己相関関数を前記計測点の移動に合わせてずらして前記相関係数を求めることを特徴とする請求項 1 3に記載の超音波診断システム。

1 9 . 前記弾性係数演算手段は、前記被検体組織を等方性弾性体及び近非圧縮性と仮定し、前記被検体組織を有限個の直方体要素に分割して 3次元有限要素モデル化し、前記各要素内では、弾性係数、応力、歪みは一様であると仮定し、弾性方程式に前記歪み分布を用いて弾性係数分布を演算することを特徴とする請求項 1 7に記載の超音波診断システム。

2 0 . 超音波探触子によって被検体の圧縮前後に計測してなる受信信号に基づいて前記被検体の組織の変位を求め、この求めた変位に基づいて前記被検体の組織の歪み分布を求めて表示手段に表示する歪み分布表示方法において、前記圧縮前後の受信信号の特徴量を求める第 1のステップと、

前記特徴量に基づき前記圧縮前後の前記特徴量の相関係数及び前記圧縮前後の前記信号間の位相差を求める第 2のステップと、

求められた前記相関係数及び前記位相差に基づいて前記圧縮に伴う前記計測点の変位並びに前記被検体の組織の歪み分布を求める第 3のステツプと、

求めた前記歪み分布を前記表示手段に表示する第 4のステップとを有してなることを特徵とする歪み分布表示方法。

2 1 . 前記特徴量は、包絡線、振幅、波高値、波数を含む波形の相関を特定可能なパラメ一夕であることを特徴とする請求項 2 0に記載の歪み分布表示方法。

2 2 . 前記第 2ステップは、前記蓄積された前記被検体の圧縮前後の前記超音波探触子により検出された信号の包絡線信号からなる超音波ビームデータに計測点を設定し、該計測点を少なくとも超音波ビーム方向に移動して、前記計測点における前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるとともに、前記圧縮前後の前記受信信号間の位相差を求め、前記第 3ステツプは、前記求められた前記相関係数が最大となる位置及び前記位相差に基づいて前記圧縮に伴う前記計測点の変位を求めることを特徴とする請求項 2 0に記載の歪み分布表示方法。

2 3 . 前記第 3のステップは、前記計測点における変位を空間微分することによって前記被検体の組織の歪み分布を求める歪み演算手段を含んでなることを特徴とする請求項 2 2に記載の歪み分布表示方法。

2 4 . 前記第 2のステップは、前記計測点を前記超音波ビーム方向に前記超音波信号の 2分の 1波長間隔で移動して前記相関係数の最大位置を求めることを特徴とする請求項 2 2に記載の歪み分布表示方法。

2 5 . 前記第 2のステツプは、前記計測点における前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるに際し、圧縮後の前記包絡線信号の自己相関関数を求め、この自己相関関数を前記計測点の移動に合わせて前記超音波信号の 2分の 1波長間隔でシフトして前記相関係数を求めることを特徴とする請求項 2 2に記載の歪み分布表示方法。

2 6 . 前記第 2のステップは、前記被検体のスライス面に対応する圧縮前後の前記包絡線信号のフレームデータに計測点を設定し、該計測点を前記フレームデータに対して少なくとも 2次元方向に移動して、前記計測点を囲む 2次元相関窓に属する前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるとともに、前記圧縮前後の前記 R F信号間の位相差を求め、

前記第 3のステップは、前記求められた前記相関係数が最大となる位置及び前記位相差に基づいて前記圧縮に伴う前記計測点の少なくとも 2次元方向の変位を求めることを特徴とする請求項 2 0に記載の歪み分布表示方法。

2 7 . 前記 2次元方向は、前記超音波探触子により受信される超音波ビ一ムの方向と、該超音波ビームの走査方向であることを特徴とする請求項 2 6に記載の歪み分布表示方法。

2 8 . 前記第 2のステップは、前記計測点を前記超音波ビーム方向に前記超音波信号の 2分の 1波長間隔で移動し、前記超音波ビームの走査方向に前記超音波ビームピッチで移動して前記相関係数の最大位置を求めることを特徴とする請求項 2 7に記載の歪み分布表示方法。

2 9 . 前記第 2のステップは、前記計測点を囲む 2次元相関窓に属する前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるに際し、圧縮後の前記包絡線信号の自己相関関数を求め、この自己相関関数を前記計測点の移動に合わせてずらして前記相関係数を求めることを特徴とする請求項 2 6に記載の歪み分布表示方法。

3 0 . 前記第 2のステップは、前記被検体の複数のスライス面に対応する圧縮前後の前記包絡線信号のポリユームデータに計測点を設定し、該計測点を前記ボリュームデータに対して 3次元方向に移動して、前記計測点を囲む 3次元相関窓に属する前記圧縮前後の前記包絡線信号の相関係数が最大となる位置を求めるとともに、前記圧縮前後の前記 R F信号間の位相差を求め、

前記第 3のステップは、前記求められた前記相関係数が最大となる位置及び前記位相差に基づいて前記圧縮に伴う前記計測点の 3次元方向の変位を求ることを特徴とする請求項 2 0に記載の歪み分布表示方法。

3 1 . 前記 3次元方向は、前記超音波探触子により受信される超音波ビームの方向と、該超音波ビームの走査方向と、これらに直交するスライス方向であることを特徴とする請求項 3 0に記載の歪み分布表示方法。

3 2 . 前記第 2のステップは、前記計測点を前記超音波ビーム方向に前記超音波信号の 2分の 1波長間隔で移動し、前記超音波ビームの走査方向に前記超音波ビームピッチで移動し、前記スライス方向に前記超音波ビームのスライスピツチで移動して前記相関係数の最大位置を求めることを特徴とする請求項 3 1に記載の歪み分布表示方法。

3 3 . 前記被検体から受信した信号で作られた仮想モデルを有限個の要素に分割して少なくとも 2次元の有限要素モデルを作成し、そのモデル化の情報と前記求めた前記歪み分布を用いて弾性係数分布を演算する第 5のステップと、この求めた弾性係数分布を前記表示手段に表示する第 6のステップとを有してなることを特徴とする請求項 2 0に記載の歪み分布表示方法。

3 4 . 前記第 5ステツプは、前記被検体を有限個の要素に分割して 3次元の有限要素モデルを作成し、そのモデル化の情報と前記歪み分布を用いて弾性係数分布を演算する弾性係数演算手段を備え、前記表示手段に前記弾性分布を表示することを特徴とする請求項 3 3に記載の歪み分布表示方法。