超音波の医療面への応用もエレクトロニクス技術の進歩と相まってさまざまな臨床領域へと広がっている。その例としては、超音波を生体情報の取得手段として利用する超音波断層像(Bモード像)やドプラ血流計測、そして超音波のエネルギーを直接利用する超音波ハイパーサーミア(温熱治療)や体外衝撃波結石破砕装置などがある。これらの中でも特に超音波Bモード像は計測のリアルタイム性、手軽さ、安全性のため臨床の場で幅広く利用されている。ここで、超音波Bモード像とは体内に超音波を放射して音響インピーダンスが異なる組織境界での反射エコーを輝度変調しながら、これを2次元断面的に走査することによって組織の形状を画像化したものである。
これに対し、組織形状だけではなく組織内の音速や減衰定数などの物理量を超音波により計測・画像化し診断に利用しようとするウルトラソニック ティシュー キャラクタライゼイション(ultrasonic tissue characterization)と呼ばれる分野がある。これは、組織の物理量を計測して組織診断に利用とするものである。そして、その中の1つとして組織の硬さ、すなわち弾性特性を計測しようとする分野があり、現在盛んに研究されている。これは、組織の弾性特性がその病理状態と深く関連しているためである。例えば、乳がんや甲状腺がんなどの硬化性がんや肝硬変、動脈硬化症などは正常組織よりも病変部分が硬くなることが知られている。そして、これまでこれらの硬さ情報は触診により得ていた。しかし、触診では客観的な情報表現が難しく、医師の経験も必要で、また計測できる領域も体表付近のある程度大きな病変に限られる。
そこで、超音波やMRIを利用して組織の弾性特性を定量的に計測・画像化しようとする研究が行われるようになった。まず、体表から機械的振動を与えその横波の伝播速度を超音波により計測し、横波の伝播速度から組織の硬さを評価する試みが行われた(非特許文献1)。これを第1の従来技術とする。この第1の従来技術は、硬い組織では横波の伝播速度が速く、軟らかい組織では横波の伝播速度が遅いということを基にしている。しかし、この方法は分解能が低いという問題点があった。
これに対し、体表から静的な圧力を加えて組織をわずかに圧縮変形させ、その際生じる組織内部の歪みを超音波により計測し、歪みから組織の弾性特性を評価する方法が1990年頃から始まった(非特許文献2)。これを第2の従来技術とする。この第2の従来技術は、硬い組織では生じる歪みが小さく、軟らかい組織では歪みが大きくなることを基にしている。
そして、その後MRIを用いて同様の原理により組織歪みから弾性特性を評価する方法が研究されるようになった。しかし、MRIを利用した方法はその特性上リアルタイム計測が困難であり、また体表から組織変形を加えることが難しいという問題点がある。従って、現在ではリアルタイム計測が可能で軽便な超音波を利用して、静的組織圧縮下における歪み推定原理に基づいた組織弾性特性評価を採用するようになってきた。
一般に、超音波とは「人間の可聴音域(約20Hz〜20kHz)より周波数の高い音」と定義されている。しかし、使う場面によっては人間の耳に聞こえる音も超音波と呼ばれることがある。そこで、最近では「超音波とは聞くことを目的としない音」と定義されるようになってきた。ただし、超音波診断装置で用いられる超音波の周波数は1MHz〜10MHzが主流である。現在、超音波は医学をはじめとして様々な分野で利用されているが、特に生体計測の分野では以下のような性質のため超音波が広く用いられ、超音波診断装置として利用されている。
・超音波は生体を媒質として伝播できる。
・超音波が生体中を進む速度(約1500m/s)は光(電磁波)に比べ桁違いに遅い。
・超音波には指向性があるため、音のビームとして利用できる。
・弱いパワーであれば生体に対して無侵襲である。
・生体の組織によって音響特性が異なるため、組織の境界で反射エコーが得られる。
図1は、超音波診断装置の原理を説明するための図である。図から明らかなように、超音波プローブ10は電気信号を超音波に、また超音波を電気信号に変換するものであり、この超音波プローブ10を用いて生体組織11内に超音波パルスを放射する。生体組織11内に放射された超音波パルスは音響インピーダンスの異なる第1の境界12で一部が反射され、反射エコー12aとして超音波プローブ10側に向かい、その残りは透過していく。透過した超音波パルスは次の音響インピーダンスの異なる第2の境界13で同様に一部が反射され、反射エコー13aとして超音波プローブ10側に向かい、その残りは透過する。このようにして反射した反射エコー(超音波エコー信号)は超音波プローブ10によって受波され電気信号に変換される。受波された反射エコー信号は、受信エコー信号のようになっている。このとき、超音波プローブ10から超音波パルスが放射されてから距離Lの位置にある反射物体14(音響インピーダンスの異なる境界)からのエコー信号を受信するまでの時間tは、
となる。ここで、cは生体内での音速であり、軟組織では1500[m/秒]にほぼ一定とみなせる。よって、超音波エコー信号を受信するまでの時間tを計測すればプローブから反射物体までの距離Lを求めることができる。
そして、電気信号に変換された超音波信号(受信エコー信号)をディスプレイに表示する方法としては、図2に示すような3種類の方法がある。図2(A)は、Aモード方式であり、表示用ディスプレイの横軸にプローブからの距離、縦軸に受信した反射エコーの強度(振幅)をとり反射エコー信号をグラフ状に表示するものである。図2(B)は、Bモード方式であり、超音波プローブを2次元断層的に走査したときに得られる反射エコー信号の強度を輝度変調し、走査位置に応じてディスプレイに表示するものである。この方式を用いると生体内の断層像が得られるため、今日最も広く用いられている。また、このとき得られる断層像をBモード像という。図2(C)は、Mモード方式であり、対象となる物体が運動している場合、超音波プローブ10の位置を固定しても時々刻々異なったAモード波形が観測される。このAモード波形を輝度変調してディスプレイの縦方向に表示し、さらに、時間に伴って横方向に走査する方式をMモードという。この方式を用いると組織の動く様子が画像化されるため、心臓の弁や壁の運動を調べるのに利用されている。
図3は、超音波プローブの種類を示す図である。Bモード像の走査方式・走査形状の違いにより現在、様々な超音波プローブが利用されている。まず、超音波ビームの走査方式の種類としては以下に示す3通りの方式がある。
第1は、手動走査方式である。これは、振動子(圧電素子)を先端に1つだけ装着したプローブを手で体表に沿わせて走査し、そのプローブの位置や角度をアームの検出機構により検出して、プローブの動きに対応した画像を表示する方式である。第2は、機械走査方式である。これは、振動子を先端に1つだけ装着したプローブをモーター等により動かし、そのプローブの位置や角度を検出機構により検出して、その動きに応じた画像を表示する方式である。第3は、電子走査方式である。これは、短冊状の振動子を先端に多数装着したプローブを用い、駆動する振動子を電子スイッチ等により制御し、走査を行う方式である。これらの走査方式の中で、現在広く用いられている方式は、電子走査方式であり、機械走査方式は一部の特殊な用途に用いられているのみである。
次に、電子走査方式のプローブでも走査形状の違いにより以下のように分けられている。第1は、セクタ走査方式である。この方式は、図3(A)に示すように、超音波ビームを扇状に走査するもので、このような走査を行うプローブをセクタスキャンプローブ(セクタフェイズドアレイプローブ)という。浅部の視野は狭いが、深部では広い範囲を観測することが可能であるため、肋骨やガス像の合間からの観察に優れている。第2は、リニア走査方式である。この方式は、図3(B)に示すように、超音波ビームを直線状に走査するもので、このような走査を行うプローブをリニアスキャンプローブ(リニアアレイプローブ)という。浅部で広い視野が得られるため、腹部検査で用いられている。第3は、オフセットセクタ走査方式である。この方式は、図3(C)に示すように、超音波ビームを扇状に走査するが、要の部分を表示しないもので、このような走査を行うプローブをコンベックススキャンプローブ(コンベックスアレイプローブ)という。浅部から深部まで広い範囲を観測できるため、腹部検査で広く用いられている。このような走査形状を持った電子走査方式の超音波プローブが現在、主に用いられている。その他、特殊なものとしては血管内部から血管周辺を観察するためのカテーテルプローブや超音波顕微鏡用の超高周波超音波プローブなどもある。また、最近では3次元の超音波像を得るための2次元アレイプローブの開発も行われている。
図4は、超音波診断装置を用いて、組織の硬さに関する情報(組織の弾性特性)を計測する手法(機械的振動下における横波伝播速度からの弾性特性評価)の一例を示す図である。これは、前述の第1の従来技術に相当するものであり、超音波を用いて組織の硬さに関する情報を計測する方式であり、組織に機械的振動を与えてその横波の伝播速度から硬さ情報を評価する方式である。この方式は、硬い組織では横波の伝播速度が速く、軟らかい組織では横波の伝播速度が遅いことを基にしている。ただし、厳密には生体組織中を伝わる横波の伝播速度は次式のように組織の密度、せん断弾性係数、せん断粘性係数、および振動の周波数に関係している。
ここで、vは横波の伝播速度、μ
1 はせん断弾性係数、μ
2 はせん断粘性係数、ρは組織の密度、ωは機械振動の角周波数である。
この方式では、まず低周波(数百ヘルツ)で振動する低周波振動子41を生体組織11の体表に接触させ、組織内部に振動を伝播させる。この振動により誘起された横波の振幅と位相の分布を血流計測に用いられるドプラ法を用いて計測する。そして、横波の振幅と位相の分布から組織の弾性特性(横波の伝播速度)を推定することになる。ただし、その際、組織の粘弾性特性は無視し、また組織の密度は一様であると仮定する。このように仮定すると組織のせん断弾性係数μ1 は、μ1=ρv2のように横波の伝播速度の2乗に比例する。
しかし、組織の粘弾性特性を無視することは難しく、組織の密度も生体内で変化するため、この方法により組織の弾性特性を定量的に評価することは難しい。また、横波の伝播速度分布も機械振動の波長程度の分解能でしか得られない。
そこで、機械的振動を与えて組織の弾性特性を評価するものに対して、前述の第2の従来技術のように、組織を静的に圧縮してその際生じる組織内の歪み分布から弾性特性を評価する方式が提案されている。これは、硬い組織では歪みが小さく、軟らかい組織では歪みが大きくなることに基づいている。
図5(A)は、静的圧縮による組織弾性計測方式の具体例を示す図である。図5(B)は、静的圧縮による組織弾性計測方式の原理を示す図である。図から明らかなように、この方式は、従来の超音波診断装置および超音波プローブ10をそのまま用いる。まず、超音波プローブ10によって組織11の圧縮前の超音波エコー信号(圧縮前RF信号)を計測する。その後、超音波プローブ10自身により組織11をわずかに(数パーセント程度)圧縮し、組織11の圧縮後の超音波エコー信号(圧縮後RF信号)を計測する。そして、計測された組織圧縮前後のRF信号から圧縮によって組織内部の各点がどれだけ動いたかという移動量である変位分布を推定する。この変位分布推定手法の主なものとしては、空間相関を用いるものとドプラの原理を用いるものがある。
図6は、空間相関法の原理を示す図である。この方法は、圧縮によって生じた組織内部の変位分布を組織圧縮前後のRF信号(またはRF信号の包絡線)から2次元相関関数を用いたテンプレートマッチングにより推定する手法である。その具体的な処理は以下のようになる。まず、組織圧縮前後のRF信号(またはその包絡線信号)をi
1(t,x),i
2(t,x)とすると、この2つの信号の相互相関係数C(t,x;n,m)は、
となる。ここで、tは超音波ビーム方向(軸方向)の座標、xはそれに直交する方向(横方向)の座標、t
0 は軸方向の相関窓サイズ、x
0 は横方向の相関窓サイズ、L
t は軸方向のサンプリング間隔、L
x は横方向のサンプリング間隔、n,mは整数である。そして、この相互相関関数が最大となるときの(n,m)を(k,l)とすると、計測点(t,x)における軸方向の変位u
y と横方向の変位u
x はそれぞれ次式のようにして求められる。
u
y =kL
t
u
x =lL
x
ただし、横方向のサンプリング間隔L
x は軸方向のサンプリング間隔L
t よりも大きいので、推定される変位成分の精度は横方向成分の方が軸方向成分よりも悪くなる。上記の処理を各計測点について行い変位分布推定する手法が空間相関法である。そのため、空間相関法では2次元の変位ベクトル成分を推定できるという特徴がある。また、組織が大変形(5%程度)した場合でも変位分布を推定できる。しかし、計算量が膨大になるため超音波計測の利点であるリアルタイム性を損なってしまう。また、変位推定精度もサンプリング間隔により制限されてしまうため、後に述べるドプラ法と比べると精度が悪いという問題点もある。
図7は、ドプラ法の原理を示す図である。この方法は、圧縮によって生じた組織内部の変位分布を組織圧縮前後のRF信号から血流計測に用いられているドプラの原理を利用して推定する手法である。その具体的な処理は以下のようになる。まず、組織圧縮前後のRF信号を次式のようにモデル化する。
ここで、i
1 (t)は圧縮前のRF信号、i
2 (t)は圧縮後のRF信号、A(t)は包絡線、ω
0 は超音波の中心角周波数、τは時間シフトである。そして、この2つのRF信号をそれぞれ直交検波すると、次式のようなベースバンド信号が得られる。
そして、この2つの信号の自己相関関数R
12(t)(本来は相互相関関数であるが共に同じ部位からの信号であるためドプラ計測では自己相関関数と呼ぶ)は次式で表される。
ここで、R
A (t)は包絡線の自己相関関数、t
0 は相関窓サイズである。また、*は複素共役を表している。よって、この自己相関関数R
12(t)の位相φ(t)から圧縮による時間シフトτ、軸方向変位u
y が次式のようにして求まる。
ただし、cは組織内の音速であり、生体内で一定と仮定する。
(本発明の原理の説明)
本発明は、生体組織の硬さを定量的に計測するにあたり、2次元又は3次元変位分布の計算時間を短縮でき、かつ、変位分布の計算精度を向上させるため、生体組織の変位を推定する公知のドプラ法を改良したことを特徴とする。ここで、まず、本発明の原理について説明する。
前述したように、空間相関法とドプラ法はそれぞれに一長一短があり、共に臨床応用に耐えられるものではない。そこで、この2つの手法の長所を組み合わせた「複合自己相関法(CA法:Combined Autocorrelation Method)」を本願の発明者等は提案している。
図8は、本願発明者等が既に提案している複合自己相関法の原理を示す図である。複合自己相関法は、ドプラ法におけるエイリアシングの問題をRF信号の包絡線による相関を用いることによって解決したものである。その具体的な処理は以下のようになる。
まず、組織圧縮前後のRF信号をドプラ法のときと同じように次式のようにモデル化する。
ここで、i
1 (t)は圧縮前のRF信号、i
2 (t)は圧縮後のRF信号、A(t)は包絡線、ω
0 は超音波の中心角周波数、τは時間シフトである。そして、この2つのRF信号をそれぞれ直交検波すると、次式のようなベースバンド信号が得られる。
そして、この2つの信号間の複素相関関数R
12(t;n)を次式のように定義する。
ここで、Tは超音波の周期、R
A (t;τ)は包絡線の自己相関関数、t
0 は相関窓サイズである。また、*は複素共役を表している。ここで、n=0の場合は、ドプラ法における自己相関関数の式(数6)に一致する。すなわち、n=0の場合はドプラ法と同じであり、軸方向変位が超音波の波長の4分の1以上になるとエイリアシングを起こしてしまう。そこで、この問題を克服するために次式で定義される包絡線相関係数C(t;n)を用いる。
ただし、R
11(t;0)は、s
1 (t)の自己相関関数、R
22(t;n)はs
2 (t+nT/2)の自己相関関数である。そして、この包絡線相関係数が最大となるnの値をkとすると、そのとき(n=k)のR
12(t;k)の位相φ(t;k)はエイリアシングの起きていない位相となる。これは、包絡線相関を計算する間隔を2分の1波長(周期)に選んだためである。ちなみに、この2分の1波長はエイリアシングを起こさないための最大の間隔である。よって、このφ(t;k)を用いることにより組織圧縮による時間シフトτ及び軸方向変位u
y は次式のように求まる。
ただし、cは組織内の音速であり、生体内で一定と仮定する。
上記の処理を各計測点について行い変位分布を推定する手法が複合自己相関法であり、ドプラ法を拡張したような手法となっている。そのため、リアルタイム計測が可能な手法となっている。また、包絡線相関を用いることによってドプラ法では計測不可能であった大変形の場合(超音波の4分の1波長以上の変位が生じる場合)の変位分布推定にも対応している。
図9は、前述の複合自己相関法の基本アルゴリズムを実行する回路構成を示すブロック図である。加圧前直交検波回路(QD)131は、加圧前のエコー信号x(t)を入力し、それぞれ直交検波して、直交検波信号Ix(t),Qx(t)信号を、第1相関演算回路133及び第1相関係数演算回路1350〜135Nに出力する。第1加圧後直交検波回路(QD)1320は、加圧後のエコー信号y(t)を入力し、それぞれ直交検波して直交検波信号Y(t)=Iy+jQy(Iy(t),Qy(t))を、第1相関係数演算回路1340及び第2相関演算回路1350に出力する。第1遅延回路134は、エコー信号y(t)を超音波の半期Tの1/2だけ遅延させ、遅延したエコー信号y1 =y(t−T/2)を第2加圧後直交検波回路(QD)1321に出力する。第2遅延回路135は、第1遅延回路134によって遅延されたエコー信号y1 =y(t−T/2)を同じく超音波の周期Tの1/2だけ遅延させ、遅延したエコー信号y2 =y(t−T)を次段の第2加圧後直交検波回路(QD)1322(図示せず)に出力する。以後、N段の遅延回路を用いて順次周期Tの整数倍だけ信号を遅延して、遅延した信号を加圧後直交検波回路に供給する。
第1相関演算回路133は、信号Ix,Qxに基づいて相関値Rxxを演算し、それを各第2相関係数演算回路1380〜138Nに出力する。第2相関演算回路1340は、加圧後直交検波回路1320からの直交検波信号Iy(t),Qy(t)を入力し、信号Iy,Qyに基づいて相関値Ryyを演算し、それを第2相関係数演算回路1380に出力する。第1相関係数演算回路1350は、加圧前直交検波回路131からの直交検波信号Ix(t),Qx(t)及び第1加圧後直交検波回路1320からの直交検波信号Iy(t),Qy(t)を入力し、複素ベースバント信号SR ,SI を求め、それを第3相関演算回路1360及び位相差演算回路1370に出力する。第3相関演算回路1360は、第1相関係数演算回路1350からの複素ベースバント信号SR ,SI を入力し、それに基づいて相関値|Rxy|を求め、それを第2相関係数演算回路1380に出力する。位相差演算回路1370は、第1相関係数演算回路1350からの複素ベースバント信号SR ,SI を入力し、それに基づいて位相差φ0 (t)を求める。第2相関係数演算回路1380は、第1相関演算回路133からの相関値Rxx、第3相関演算回路1360からの相関値|Rxy|、並びに第2相関演算回路1340からの相関値Ryyを入力し、これらの各相関値に基づいて相関係数C0 (t)を演算し、出力する。
第2加圧後直交検波回路(QD)1321は、第1遅延回路134によって遅延されたエコー信号y1 =y(t−T/2)を入力し、それぞれ直交検波して直交検波信号Y1 (t)=Iy1 +jQy1 (Iy1 (t),Qy1 (t))を、第1相関係数演算回路1341及び第2相関演算回路1351に出力する。第2相関演算回路1341は、第2加圧後直交検波回路(QD)1321からの直交検波信号Iy1 (t),Qy1 (t)を入力し、その信号Iy1 (t),Qy1 (t)に基づいて相関値Ry1y1を演算し、それを第2相関係数演算回路1381に出力する。第1相関係数演算回路1351は、加圧前直交検波回路131からの直交検波信号Ix(t),Qx(t)、第2加圧後直交検波回路(QD)1321からの直交検波信号Iy1 (t),Qy1 (t)を入力し、複素ベースバント信号SR1,SI1を求め、それを第3相関演算回路1361及び位相差演算回路1371に出力する。第3相関演算回路1361は、第1相関係数演算回路1351からの複素ベースバント信号SR1,SI1を入力し、それに基づいて相関値|Rxy1 |を求め、それを第2相関係数演算回路1381に出力する。位相差演算回路1371は、第1相関係数演算回路1351からの複素ベースバント信号SR1,SI1を入力し、それに基づいて位相差φ1 (t)を求める。第2相関係数演算回路1381は、第1相関演算回路133からの相関値Rxx、第3相関演算回路1361からの相関値|Rxy1 |、並びに第2相関演算回路1341からの相関値Ry1y1を入力し、これらの各相関値に基づいて相関係数C1 (t)を演算し、出力する。
以下同様に、第1遅延回路135以降の第2加圧後直交検波回路(QD)1322〜132N、第2相関演算回路1342〜134N、第1相関係数演算回路1352〜135N、第3相関演算回路1362〜136N、位相差演算回路1372〜137N及び第2相関係数演算回路1382〜138Nは、上述の1段目及び2段目の回路群と同様の処理を実行し、相関係数C2 (t)〜CN (t)及び位相φ2 (t)〜φN (t)を出力する。
上述の複合自己相関法の基本アルゴリズムを実行する回路は、加圧後のエコー信号y(t)を遅延回路134〜13Nで超音波の周期T/2(2分の1波長)だけ遅延し、それを直交検波回路(QD)1320〜132Nを用いて個別に直交検波している。
前述のように組織圧縮に伴う変位分布が推定されたら、それを空間微分することにより歪み分布が得られる。歪み分布は定性的に組織の弾性特性を表しているものであり、歪み分布からでもかなりの弾性特性に基づいた診断は行える。しかし、肝硬変などの病変部全体が硬くなるような場合には、定量的な弾性係数によって評価しなければ組織診断は難しい。そのため、近年、組織弾性係数分布再構成法についても研究されるようになってきた。しかし、今のところスタンダードな手法はなく、いずれの手法も研究段階であるというのが実状である。
組織弾性係数分布は先にも述べたように組織内部の歪み分布と応力分布から求められる。しかし、応力分布を直接計測することは現状では困難であるため、歪み分布と組織圧縮の際の境界条件から逆問題的に弾性係数分布を再構成することになる。そのため、一般的に逆問題を解くことは難しく、現在提案されている弾性係数再構成法も数少ない。従来から提案されている弾性係数再構成法を以下に説明する。
第1に、1次元を仮定した方法(1次元弾性体を仮定)がある。これは、1次元弾性体を仮定して歪みの逆数を弾性係数とみなす方法である。この方法は弾性係数再構成法ではなく、歪みの逆数を求めるだけであるので、歪みにおける非定量性をそのまま残している。
第2に、弾性方程式から応力項を消去する方法(等方性弾性体、非圧縮性、平面歪み状態を仮定)がある。これは、平面歪み状態を仮定した場合の弾性方程式を変形し、応力項を消去した方程式を用いて組織圧縮の際の境界条件(体表での外部圧力分布、または体表での変位)と歪み分布(せん断歪み成分を含む歪みテンソルの全成分)から組織弾性係数分布を再構成する手法である。ただし、絶対的な弾性係数を推定するには、弾性係数が前もってわかっている領域(参照領域)が必要となる。
第3に、弾性微分方程式を積分する方法(等方性弾性体、非圧縮性、平面応力状態を仮定)がある。これは、平面応力状態を仮定した場合の弾性方程式を変形した応力項を含まない弾性係数に関する微分方程式を体表付近での弾性係数を基準として順次積分していくことにより、歪み分布(せん断歪み成分を含む歪みテンソルの全成分)から組織弾性係数分布を再構成する方法である。そのため、体表付近の弾性係数分布が前もって分かっている領域が必要であり、また体表付近を基準として積分を行っていくので奥に行くほど誤差が積算されるという問題点もある。
第4に、摂動法を用いた手法(等方性弾性体、近非圧縮性、平面歪み状態を仮定)がある。これは、平面歪み状態を仮定した場合の弾性方程式を基にした摂動法により体表での外部圧力分布と超音波ビーム方向(軸方向)の歪み分布とから反復的に組織弾性係数分布を再構成する方法である。
前述の複合自己相関法は、ドプラ法と同じ1次元の処理を基にしているため、超音波プローブが相対的に横方向に移動してしまい、超音波ビーム方向(軸方向)に直交する方向(横方向)の変位が生じてしまった場合(横方向変位が超音波ビーム幅を超えてしまう場合)には、組織圧縮後のRF信号が無相関となってしまい変位推定に失敗してしまうという問題点がある。すなわち、横方向の変位に対応することができずに、超音波ビーム方向(軸方向)の変位成分のみしか推定できないという問題がある。
前述の組織弾性係数分布再構成法は、いずれの方法も問題を簡単化するために2次元状態(平面歪み状態や平面応力状態)を仮定している。しかし、本来3次元構造をなす生体組織を2次元で近似すると弾性係数を過小評価する恐れがあり好ましくない。これは、2次元状態では考察面に垂直な方向の歪みや応力を考慮に入れていないことが原因である。
本発明は、生体組織の硬さを定量的に計測するにあたり、2次元又は3次元変位分布の計算時間を短縮でき、かつ、変位分布の計算精度を向上させること、及び生体組織の弾性係数を精度よく計測することを特徴とする。
(実施の形態1)
以下に、本発明の実施形態1について、図面を参照して説明する。
図10は、本発明に係る超音波診断システムの好ましい一例の実施形態1の超音波診断システムの概略構成を示すブロック図である。この超音波診断システムでは、包絡線相関計算の際、複合自己相関法で1次元の相関窓で1次元探索していた処理を2次元の相関窓を用いて2次元探索することにより横方向の移動にも対応した拡張複合自己相関法と呼ばれる方法を採用している。この拡張複合自己相関法は、軸方向には2分の1波長間隔、横方向にはライン間隔の格子点でのみ包絡線相関計算を行うことにより計算量を減少させて高速化を図っている。ただし、複合自己相関と同様に拡張複合自己相関法でも位相情報を利用して軸方向の変位推定精度を向上させている。しかし、横方向変位の推定はキャリアとなる信号がないため位相情報は利用できない。そのため、横方向変位推定精度は空間相関法と同様に横方向サンプリング間隔(ライン間隔)により制限されてしまう。しかし、後で提案する弾性係数分布再構成法では軸方向の歪み(変位)分布のみから弾性係数分布を推定できるため、ここでは横方向変位推定精度の向上は特に行わない。この拡張複合自己相関法の具体的な構成について図10を用いて説明する。
図10において、超音波プローブ91は、被検体内へ超音波を送波すると共にその反射波を受波するものであり、従来のセクタスキャンプローブ(セクタフェイズドアレイプローブ)、リニアスキャンプローブ(リニアアレイプローブ)又はコンベックススキャンプローブ(コンベックスアレイプローブ)などである。超音波プローブ91からは、組織圧縮前後のRF信号が直交検波器92に出力される。直交検波器92は、組織圧縮前後のRF信号をそれぞれ組織圧縮前後の複素包絡線信号(IQ信号)に変換し、複素2次元相関計算部93に出力する。複素2次元相関計算部93は、組織圧縮前後のRF信号間における2次元相関を計算し、その相関が最大となる位置を横方向変位計算部94及び軸方向変位計算部95に出力し、そのときの相関関数の位相を軸方向変位計算部95に出力する。ただし、軸方向にはエイリアシングを起こさずに位相を検出できる最大の間隔である超音波中心周波数の2分の1波長間隔でのみ相関を計算するものとする。これは、超音波診断システムのリアルタイム表示を優先させるためである。従って、高精度な相関を計算するためには、この2分の1波長間隔に限定する必要はない。
横方向変位計算部94は、複素2次元相関計算部93からの横方向の相関最大位置に基づいて横方向の変位ux を計算し、それを横方向歪み計算部96に出力する。一方、軸方向変位計算部95は、複素2次元相関計算部93からの軸方向の相関最大位置及びそのときの位相に基づいて軸方向の変位uy を計算し、それを軸方向歪み計算部97に出力する。横方向歪み計算部96は、横方向変位計算部94からの横方向変位ux の分布を空間的に微分することにより横方向歪み分布εx を計算し、それを量子化部98に出力する。一方、軸方向歪み計算部97は、軸方向変位計算部95からの横方向変位uy の分布を空間的に微分することにより軸方向歪み分布εy を計算し、それを量子化部98に出力する。量子化部98は、横方向歪み分布εx 及び軸方向歪み分布εy をグレースケール表示(又はカラー表示)するために各歪み分布を量子化し、表示部99に出力する。表示部99は、量子化された各歪み分布を表示する。
次に、図10の超音波診断システムで採用した拡張複合自己相関法の動作について説明する。まず、組織圧縮が極僅か(数パーセント以下)である場合、組織内部を局所的に見れば平行移動したと見なすことができ、組織圧縮前後のRF信号を次式のようにモデル化できる。
ここで、i
1 (t,x)は圧縮前のRF信号、i
2 (t,x)は圧縮後のRF信号、A(t,x)は包絡線、ω
0 は超音波の中心角周波数、τは時間シフト、u
x は横方向変位である。また、ここではドプラ法や複合自己相関法のときと違い横方向の変位も考慮して圧縮前後のRF信号をモデル化している。そして、この式の中で最終的に知りたいパラメータは、軸方向の変位u
y =cτ/2(すなわち、時間シフトτ)と横方向変位u
x である。ただし、cは組織内の音速であり、生体内で一定と仮定する。
そこで、まずこれらの組織圧縮前後のRF信号を直交検波器92でそれぞれ直交検波する。すなわち、各RF信号に超音波の中心周波数と同じ周波数のsin波,cos波をかけ、それぞれ低域通過フィルタをかける。すると、以下のような複素ベースバンド信号s
1 ,s
2 が得られる。
そして、このs
1 (t,x)とs
2 (t+nT/2,x+mL)との間の2次元複素相関関数R
12(t,x;n,m)を次式のように定義する。
ここで、Tは超音波の周期、Lは横方向サンプリング間隔(ライン間隔)、R
A (t,x;τ,u
x )は包絡線の自己相関関数、t
0 は軸方向相関窓サイズ、x
0 は横方向相関窓サイズである。また、*は複素共役を表している。そして、この2次元複素相関関数を用いて2次元包絡線相関係数C(t,x;n,m)を以下のように定義する。
ただし、R
11(t,x;0,0)はs
1 (t,x)の自己相関関数、R
22(t,x;n,m)はs
2 (t+nT/2,x+mL)の自己相関関数である。そして、この包絡線相関係数を用いて複合自己相関法の場合と同様にエイリアシングの問題を克服する。すなわち、各計測点(t,x)におけるC(t,x;n,m)とR
12(t,x;n,m)の位相φ(t,x;n,m)との組{C(t,x;n,m),φ(t,x;n,m)}をすべてのnとmについて求める。ここで、nとmの範囲が十分広ければ、すなわち、包絡線相関を行う探索範囲が十分に大きければ、包絡線相関係数が最大となる(n,m)=(k,l)に対する位相φ(t,x;k,l)はエイリアシングの起きていない位相となる。これは、包絡線相関C(t,x;n,m)が最大となる(n,m)=(k,l)のとき、s
1 (t,x)とs
2 (t+kT/2,x+lL)との時間シフトの大きさ|τ−kT/2|がT/2よりも小さくなる、すなわち、|φ(t,x;k,l)|=ω
0 |τ−kT/2|がπよりも小さくなるためである。よって、このエイリアシングの起きていないφ(t,x;k,l)を用いれば、計測点(t,x)における正確な時間シフトτ、軸方向変位u
y 、そして横方向変位u
x が次式のように求まる。
ただし、cは組織内での音速(ここでは軟組織における一般的な音速1500m/sで一定とする)である。したがって、組織内のすべての点で上記のように軸方向変位と横方向変位を計算すれば、軸方向変位分布u
y (x,y)と横方向変位分布u
x (x,y)が得られる。
また、各変位分布を次式のように空間微分することにより、軸方向歪み分布ε
y (x,y)と横方向歪み分布ε
x (x,y)が次式のように求められる。
以上のようにして、実施形態1によれば、組織圧縮前後のRF信号から軸方向と横方向の変位(歪み)分布を推定することができる。ただし、上式のu
x =lLからもわかるように横方向変位の精度は横方向サンプリング間隔(ライン間隔)によって制限されてしまうため、精度は若干劣るということはあるが、リアルタイムに観察できるので実用性の高いものである。
(実施の形態2)
実施の形態1の拡張複合自己相関法は、組織の横方向移動に対応するように2次元の相関窓と探索範囲を用いて、組織圧縮の際の組織と超音波プローブの相対的な横方向移動には対応している。しかし、組織圧縮の際に軸方向と横方向にそれぞれ垂直な方向(2次元超音波走査面に垂直な方向(スライス方向))の変位が生じてしまい組織移動が起こった場合には、2次元の拡張複合自己相関法では歪みの推定を行うことができない。つまり、実施の形態1では、問題を簡単化するために2次元状態(平面歪み状態や平面応力状態)を仮定しているが、本来3次元構造をなす生体組織を2次元で近似すると弾性係数を過小評価する恐れがあり好ましくない。これは、2次元状態では考察面に垂直な方向の歪みや応力を考慮に入れていないことが原因である。
そのため、より安定したシステムにするには、上述の拡張複合自己相関法を、以下に述べる実施の形態2に示すように、3次元の相関窓と3次元の探索範囲を用いることにより簡単に拡張することが可能である。
図11及び図12は、本発明の実施の形態2の3次元複合自己相関法の基本アルゴリズムを示すフローチャート図である。なお、図12は、図11の処理の一部の詳細を示すフローチャート図である。
ステップS11では、ステップS23の判定処理と組み合わせて、第1番目の走査線から第L番目の走査線についてそれぞれ同様の処理を行うために、走査線番号レジスタlに「1」を格納する。
ステップS12では、ステップS18の判定処理と組み合わせて、厚み方向(フレーム)を−UからUまで順次シフトする処理を実行する。
ステップS13では、ステップS17の判定処理と組み合わせて、方位方向(走査線)を−VからVまで順次シフトする処理を実行する。
ステップS14では、ステップS16の判定処理と組み合わせて、距離方向(軸方向)を0からMまで順次シフトする処理を実行する。
ステップS15では、複合自己相関法により、距離方向(軸方向)における包絡線の相関係数C(l,t;u,v,n)を計算する。この複合自己相関法は、従来の方法でやってもいいが、それだと計算に時間を要するので、ここでは、高速化された複合自己相関法を用いて相関係数C(l,t;u,v,n)の計算を行う。この高速化複合自己相関法については後述する。
ステップS16では、前のステップS14と組み合わせられた処理であり、距
離方向レジスタnがその最大値Mに達したか否かの判定を行い、達した場合には
ステップS17に進み、そうでない場合はステップS14にリターンし、距離方向レジスタnをインクリメント処理する。
ステップS17では、前のステップS13と組み合わせられた処理であり、方位方向レジスタvがその最大値Vに達したか否かの判定を行い、達した場合にはステップS18に進み、そうでない場合はステップS13にリターンし、方位方向レジスタvをインクリメント処理する。
ステップS18では、前のステップS12と組み合わせられた処理であり、厚み方向レジスタuがその最大値Uに達したか否かの判定を行い、達した場合にはステップS19に進み、そうでない場合はステップS12にリターンし、厚み方向レジスタuをインクリメント処理する。
ステップS19では、ステップS12〜ステップS18の処理によって求めれた相関係数C(l,t;u,v,n)(u=−U,…0,…,U)(v=−V,…0,…,V)(n=0,1,…,N)の中から最大となる(u,v,n)を求め、それを(u0 ,v0 ,n0 )とする。
ステップS20では、ステップS19で求められた相関係数C(l,t;u0 ,v0 ,n0 )について、その位相差φ(l,t;u0 ,v0 ,n0 )を計算する。
ステップS21では、最終的を位相差として、φ=n0 π+φ(l,t;u0 ,v0 ,n0 )を計算する。
ステップS22では、(u0 ,v0 ,n0 )の近傍の相関係数C(l,t;u,v,n)を用いて、方位方向の変位:v=v0 +Δv及び厚み方向の変位:u=u0 +Δuを計算する。
ステップS23では、前のステップS11と組み合わせられた処理であり、走査線番号レジスタlがLに達したか否かの判定を行い、達した場合にはステップS24に進み、そうでない場合はステップS11にリターンし、走査線番号レジスタlをインクリメント処理する。
ステップS24では、組織圧縮に伴う変位分布が推定されたら、それを空間微分することにより歪み分布を計算する。
図13は、図12のステップS15の高速化された複合自己相関法の詳細を示すフローチャート図である。
ステップS31では、圧縮前のRF信号の包絡線xと、圧縮後のRF信号の包絡線をそれぞれ直交検波して、以下のようにI,Q信号を求める。
x(t)→Ix、Qx(X(t)=Ix+jQxとする)
y(t)→Iy、Qy(Y(t)=Iy+jQyとする)
ステップS32では、相関:Rxy、Rxx、Ryyを次式に基づいて計算する。
Rxy=∫X(t+ν)・Y* (t+ν)dν
Rxx=∫X(t+ν)・X* (t+ν)dν
Ryy=∫Y(t+ν)・Y* (t+ν)dν
ステップS33では、求められた相関Rxy、Rxx、Ryyを用いて相関係数C0 を次式に基づいて計算する。
C0 =|Rxy|/√Rxx√Ryy
ステップS34では、変数nに1をセットする。
ステップS35では、Yn (t)=Y(t−nT)ejωnTを計算する。
ステップS36では、次式に基づいてRxyn ,Rynynを計算する。
Rxyn =∫X(t+ν)・Yn * (t+ν)dν
=∫X(t+ν)・Y* (t−nT+ν)ejωnTdν
Rynyn=∫Yn(t+ν)Yn * (t+ν)dν
=∫Y(t−nT+ν)・Y* (t−nT+ν)dν
ステップS37では、求められたRxyn ,Rynynを用いて相関係数Cn を次式に基づいて計算する。
Cn =|Rxyn |/√Rxx√Rynyn
ステップS38では、変数nをインクリメント処理する。
ステップS39では、変数nが最大値Mに達したか否かを判定し、達した場合はこの処理を終了し、達していない場合は、ステップS35にリターンし、同様の処理を繰り返す。
図13のフローチャートでは、Rxyn ,Rynynを求めるのに、ステップS35でYnをYから導いている。このために、回路構成の簡略化を図ることができる。以下、どのようにしてYn をYから導くかについて説明する。
まず、加圧前のエコー信号x(t)を
x(t)=u(t)cos(ωt+θ)
軸方向に加圧後のエコー信号y(t)を
y(t)=x(t+τ)=u(t+τ)cos(ω(t+τ)+θ)
とする。
各信号x,y,yn の直交検波値は、
x(t)→Ix=0.5u(t)cosθ
Qx=−0.5 u(t)sinθ
(X(t)=Ix+jQx=0.5u(t)e-jθ)
y(t)→Iy=0.5u(t+τ)cos(ωτ+θ)
Qy=−0.5u(t+τ)sin(ωτ+θ)
(Y(t)=Iy+jQy=0.5u(t+τ)e-j(ωτ+θ) )
yn (t)=y(t−nT)
=u(t+τ−nT)cos(ω(t+τ−nT)+θ)
=u(t+τ−nT)cos(ωt+ω(τ−nT)+θ)
となる。ここで、Tは半周期なので、
Iyn=0.5u(t+τ−nT)cos(ω(τ−nT)+θ)
Qyn=−0.5u(t+τ−nT)sin(ω(τ―nT)+θ)
(Yn=Iyn+jQyn=0.5u(t+τ−nT)e-j(ω(τ-nT)+θ))
となる。以上の式から以下のような関係が成り立つ。
Yn (t)=Iyn+jQyn
=0.5u(t+τ−nT)e-j(ω(τ-nT)+θ)
=Y(t−nT)ejωnT
これからYn (t)はY(t)=Iy+jQyから求まることになる。
従って、Rxyn,Rynynも、次式のようにX、Yから求めることができる。
Rxyn =4∫X(t+ν)・Yn * (t+ν)dν
=4∫X(t+ν)・Y*(t−nT+ν)ejωnTdν
|Rxyn |=Run
=∫u(t+ν)u(t+τ−nT+ν)dν
=4∫|X(t+ν)・Yn * (t+ν)|dν
=4∫|X(t+ν)・Y*(t−nT+ν)ejωnT|dν
=4∫|X(t+ν)・Y*(t−nT+ν)|dν
Rynyn=∫u(t+τ−nT+ν)u(t+τ−nT+ν)dν
=4∫|Yn(t+ν)・Yn * (t+ν)|dν
=4∫Y(t−nT+ν)・Y*(t−nT+ν)dν
ここで、*は複素共役を表している。
図14は、この3次元複合自己相関法の基本アルゴリズムを実行する回路構成を示すブロック図である。複合自己相関法を実行する回路構成を図9に示すようなものにすると、直交検波回路1320〜132Nが多段接続されることによって、直交検波回路1320〜132Nの処理に時間を要する。そのため、その処理時間が膨大なものとなってしまい、高速な演算処理の妨げとなり、リアルタイムな画像表示の妨げとなる。そこで、前述のような基本アルゴリズムに応じた図14のような回路構成を採用することによって、複合自己相関法を実行する回路の処理速度を高速化している。
加圧前直交検波回路(QD)131は、加圧前のエコー信号x(t)を入力し、それぞれ直交検波して、直交検波信号Ix(t),Qx(t)信号を、第1相関演算回路133及び第1相関係数演算回路1350〜135Nに出力する。加圧後直交検波回路(QD)132は、加圧後のエコー信号y(t)を入力し、それぞれ直交検波して直交検波信号Y(t)=Iy+jQy(Iy(t),Qy(t))を、第1相関係数演算回路1350、第2相関演算回路1340及び第1遅延回路134及び第2遅延回路135に出力する。第1遅延回路134及び第2遅延回路135は、直交検波信号Y(t)をそれぞれ超音波の周期Tの1/2だけ遅延させ、遅延した直交検波信号Y(t−T/2)を第1相関係数演算回路1351、第3遅延回路136及び第4遅延回路137に出力する。第3遅延回路136及び第4遅延回路137は、直交検波信号Y(t−T/2)をそれぞれ超音波の周期Tの1/2だけ遅延させ、遅延した直交検波信号Y(t−T)を次段の第1相関係数演算回路及び遅延回路(図示せず)に出力する。以後、N段の遅延回路を用いて順次周期Tの1/2の整数倍だけ信号を遅延して、遅延した信号を第1相関係数演算回路に供給する。
第1相関演算回路133は、信号Ix,Qxに基づいて相関値Rxxを演算し、それを各第2相関係数演算回路1380〜138Nに出力する。第2相関演算回路1340は、加圧後直交検波回路132からの直交検波信号Iy(t),Qy(t)を入力し、信号Iy,Qyに基づいて相関値Ryyを演算し、それを第2相関係数演算回路1380に出力する。第1相関係数演算回路1350は、加圧前直交検波回路131からの直交検波信号Ix(t),Qx(t)及び加圧後直交検波回路132からの直交検波信号Iy(t),Qy(t)を入力し、複素ベースバント信号SR ,SI を求め、それを第3相関演算回路1360及び位相差演算回路1370に出力する。第3相関演算回路1360は、第1相関係数演算回路1350からの複素ベースバント信号SR ,SI を入力し、それに基づいて相関値|Rxy|を求め、それを第2相関係数演算回路1380に出力する。
位相差演算回路1370は、第1相関係数演算回路1350からの複素ベースバント信号SR ,SI を入力し、それに基づいて位相差φ0 (t)を求める。第2相関係数演算回路1380は、第1相関演算回路133からの相関値Rxx、第3相関演算回路1360からの相関値|Rxy|、並びに第2相関演算回路1340からの相関値Ryyを入力し、これらの各相関値に基づいて相関係数C0 (t)を演算し、出力する。
第2相関演算回路1341は、第1遅延回路134及び第2遅延回路135からの遅延後の直交検波信号Iy(t−T/2),Qy(t−T/2)を入力し、信号Iy(t−T/2),Qy(t−T/2)に基づいて相関値Ry1y1を演算し、それを第2相関係数演算回路1381に出力する。第1相関係数演算回路1351は、加圧前直交検波回路131からの直交検波信号Ix(t),Qx(t)、第1遅延回路134及び第2遅延回路135からの遅延後の直交検波信号Iy(t−T/2),Qy(t−T/2)を入力し、複素ベースバント信号SR1,SI1を求め、それを第3相関演算回路1361及び位相差演算回路1371に出力する。第3相関演算回路1361は、第1相関係数演算回路1351からの複素ベースバント信号SR1,SI1を入力し、それに基づいて相関値|Rxy1 |を求め、それを第2相関係数演算回路1381に出力する。位相差演算回路1371は、第1相関係数演算回路1351からの複素ベースバント信号SR1,SI1を入力し、それに基づいて位相差φ1 (t)を求める。第2相関係数演算回路1381は、第1相関演算回路133からの相関値Rxx、第3相関演算回路1361からの相関値|Rxy1 |、並びに第2相関演算回路1341からの相関値Ry1y1を入力し、これらの各相関値に基づいて相関係数C1 (t)を演算し、出力する。
第3遅延回路135及び第4遅延回路136から次段の第2相関演算回路1342〜134N、第1相関係数演算回路1352〜135N、第3相関演算回路1362〜136N、位相差演算回路1372〜137N及び第2相関係数演算回路1382〜138Nは、上述と同様の処理を順次遅延された遅延後の直交検波信号Iy(t−2T/2)・・・Iy(t−NT/2),Qy(t−2T/2)・・・Qy(t−NT/2)に対して実行し、相関係数C2 (t)〜CN (t)及び位相φ2 (t)〜φN (t)を出力する。すなわち、言い換えれば、実施の形態2は、被検体との間で超音波信号を送受信する超音波探触子により検出された超音波信号のフレームデータからなるボリュームデータを記憶しておき、記憶された被検体の圧縮前後のボリュームデータに計測点を設定し、この計測点をボリュームデータに対して3次元方向に移動して、計測点を囲む3次元相関窓に属する圧縮前後の相関係数が最大となる位置及び位相差を求め、求められた相関係数が最大となる位置及び位相差に基づいて圧縮に伴う計測点の変位を求めることを特徴とする。さらに、求めた変位に基づいて、被検体の歪み分布を求めることができる。
次に、3次元有限要素モデルを用いた弾性係数分布再構成法について説明する。弾性係数分布再構成逆問題を簡単化するため、この実施の形態では組織をモデル化する。これはまた、本実施の形態の弾性係数分布再構成法において有限要素法を用いるためでもある。この実施の形態では、組織を以下のように仮定及びモデル化する。
まず、組織を等方性弾性体と仮定する。組織歪み分布を推定する際、外部から圧力を加えて組織を静的に圧縮する。しかし、組織圧縮前後のRF信号間の相関を保つために、微小圧縮しか行わない。そのため、この場合、応力と歪みの間の関係はほぼ線形である。すなわち、組織を弾性体として近似できる。また、今回対象としている軟組織はほぼ等方性が成り立つため、この実施の形態では組織を等方性弾性体と仮定する。
さらに、組織を近非圧縮性と仮定する。生体組織は、非圧縮性(ポアソン比ν=0.5)に近いことが知られている。そこで、ポアソン比を0.49とし、生体内で一定とする。ここで、完全な非圧縮性を仮定しないのは、ポアソン比ν=0.5とすると特殊な弾性方程式となり、本実施の形態で用いている有限要素法が適用できなくなるためである。そして、ポアソン比を生体内で一定とすることで弾性係数分布推定の推定パラメータをヤング率のみとすることができ、逆問題を簡単化できる。また、ポアソン比はヤング率に比べ生体中であまり変化しないパラメータであるため、この実施の形態ではポアソン比を0.49で一定とする。
組織を3次元有限要素モデル化する。この弾性係数分布再構成法では有限要素法を用いるため、組織を有限個の直方体要素に分割する。そして、各要素内では、弾性係数、応力、歪みは一様であると仮定する。一般的に逆問題を解くには、それに対応する順問題を理解することが重要である。今回の歪み分布と境界条件から弾性係数分布を推定する逆問題の場合、それに対応する順問題とは、弾性係数分布と境界条件から歪み分布を求めることである。そして、この順問題の数値解法の1つが有限要素法(FEM:Finite Element Method)である。
ここで、有限要素法とは対象となる連続体を有限個の要素の集合で近似し、この集合体に対して成り立つ連立1次方程式を数値的に解く手法のことである。なお、有限要素法の定式化については後述する。ここでは有限要素法とは「入力として物体の弾性係数分布と境界条件を与えれば、出力としてそのときの歪み(変位)分布と応力分布が得られるもの」として捉えておけば十分である。
この実施の形態では、組織を等方性弾性体で近似するため、組織内では以下のような弾性方程式(つりあい方程式・歪み−変位関係式・応力−歪み関係式)が成り立つ。
つりあい方程式は次式のように表される。
歪み−変位関係式は次式のように表される。
応力−歪み関係式(一般化したフックの法則)は次式のように表される。
上記の式ではテンソル表現を用いており、実際にはつりあい方程式として3式、歪み−変位方程式として6式、応力−歪み関係式として6式が存在する。また、座標系(x
1 ,x
2 ,x
3 )は(x,y,z)を表しており、その他の記号は以下のことを表している。
E:ヤング率(弾性係数とはヤング率のことを表している)
ν:ポアソン比
ε
ij:歪みテンソル
(ε
nn=ε
11+ε
22+ε
33:体積歪み)
σ
ij:応力テンソル
δ
ij:クロネッカーのデルタ
u
i :変位ベクトル
f
i :体積力ベクトル(重力の影響は無視できるため、ここではf
i =0とする)
ここで、応力−歪み関係式をε
ijについて整理すると、次のような歪み−応力関係式が得られる。
ただし、σ
nn=σ
11+σ
22+σ
33である。よって、この式の中でi=j=2とし、ヤング率Eについて整理すると次式が得られる。
従って、軸方向(この実施の形態では、y方向を超音波ビーム方向、すなわち軸方向とする)の歪み成分と全方向の応力成分がわかれば、ヤング率すなわち弾性係数を求めることができる。なお、上述の計算式からは、応力分布を直接計測することは現状では困難である。そこで、この実施の形態では応力分布と弾性係数分布を交互に推定・更新しながら、推定弾性係数分布を実際の分布に近づけていく。その弾性係数分布再構成の手順は、以下のようになる。
第1に、未知弾性係数分布の初期値分布{E^
0}として一様分布を考える。第2に、初期弾性係数分布{E^
0}のときに生じる応力分布{σ^
0}を3次元有限要素法により求める。具体的には、まず組織モデル内の各要素に対して歪み−変位関係式及び応力−歪み関係式をつりあい方程式に代入して得られる次式のようなつりあい方程式を作る。
ただし、
である。そして、この連立方程式を以下のような境界条件のもとガウスの消去法を用いて変位について解き、弾性係数分布{E^
0}のときの変位分布{u^
0}を求める。
上式において、piは体表における外部圧力ベクトル、σnは側面に垂直な方向の応力成分である。また、上段の式は底面が固定されていることを示し、中段の式は体表での応力分布は外部圧力分布に等しいことを示し、下段の式は側面が拘束されていないことをそれぞれ示している。次に、この変位分布{u^
0}を歪み−変位関係式に代入して、弾性係数分布{E^
0}のときの歪み分布{ε^
0}を求める。そして、この歪み分布{ε^
0}を応力−歪み関係式に代入することにより、弾性係数分布{E^
0}のときの応力分布{σ^
0}を得る。
第3に、3次元有限要素法により得られた応力分布と拡張複合自己相関法により推定した軸方向(y方向)歪み分布{ε
y }を用いて、弾性係数分布{E^
k}を次式により更新する。
ただし、この式は、上述の応力−歪み関係式をε
ijについて整理し、式中のi=j=2とし、ヤング率Eについて整理した式を書き改めたものであり、式中のkは繰り返し回数を表している。
第4に、上述のように更新された弾性係数分布と上述の境界条件を用いて再び3次元有限要素解析を行い、応力分布を更新する。
そして、第3及び第4の処理を繰り返すことにより弾性係数分布を実際の分布に近づけていく。ただし、次式の条件が満たされた時点で弾性係数分布推定は収束したとみなし、推定を終了する。
ここで、lは要素番号、Nは要素数、Γはしきい値である。
以上が、実施の形態2の3次元有限要素モデルによる弾性係数分布再構成法であり、この方法は3次元のつりあい方程式を基に弾性係数分布を推定している。そのため、本手法は従来の手法よりもより実際の生体組織に近い仮定に基づいているので、より正確な弾性係数推定が可能になる。また、本手法は精度良く推定可能な軸方向の歪み分布のみから弾性係数分布を再構成するため、安定した弾性係数分布再構成が行える。ただし、本手法は組織弾性係数の3次元分布を推定する手法であるため、2次元アレイ超音波プローブを用いるか、1次元アレイ超音波プローブをスライス方向に機械的に走査することにより、対象を3次元的に走査する必要がある。
(実施形態1,2のシミュレーション)
本発明の実施の形態1,2の拡張複合自己相関法と3次元有限要素モデルによる弾性係数分布再構成法の有効性をシミュレーションによって実証する。図15は、このシミュレーションの手順の概略を示す図である。
第1に、推定したい弾性係数分布を持つ組織モデルを作成する。このとき、組織モデル内には超音波エコー信号を発生させるための散乱体を分布させておく。第2に、この組織モデルに対して外部圧力を加え、計算機上で組織圧縮を行う。そして、この圧縮による各散乱点の移動先を有限要素法などにより求める。第3に、組織モデル圧縮前後の散乱体分布を基に圧縮前後のRF信号を作成する。第4に、この圧縮前後のシミュレーションRF信号に対して拡張複合自己相関法を適用し、組織歪み分布を推定する。第5に、拡張複合自己相関法により推定された歪み分布と組織モデル圧縮の際に設定した境界条件(外部圧力分布など)とから3次元有限要素モデルによる弾性係数分布再構成法により組織弾性係数分布を推定する。
今回のシミュレーションで用いた組織モデルの弾性係数分布は、各シミュレーションにおいて異なるが、いずれの場合も等方性弾性体を仮定する。なお、各シミュレーションで設定した弾性係数の値としては、本実施の形態の組織弾性計測システムで主な対象としている***組織の弾性係数にほぼ即している。また、組織圧縮前後のシミュレーションRF信号を作成するために、各組織モデルには点散乱体を分布させた。その際、点散乱体の平均密度としては500個/cm3とし、組織圧縮前の散乱体の位置は一様乱数により、散乱係数は平均1.0、標準偏差0.3の正規乱数により決めた。そして、組織圧縮後の散乱***置は有限要素解析の結果に応じて組織圧縮前の各散乱体を移動させることにより決めている。ここで、実際の組織の散乱体に関する情報は未知であるが、シミュレーションRF信号を基にBモード像にした際、実際の組織のBモード像に近くなるように各パラメータを設定する。
この実施の形態では、組織モデルに対する組織圧縮前後のシミュレーションRF信号を次式のように組織圧縮前後の散乱体関数と超音波システムの点広がり関数との畳み込みにより作成する。
ここで、i
1 (x,y,z)は組織圧縮前のRF信号、i
2 (x,y,z)は組織圧縮後のRF信号、h(x,y,z)は超音波システムの点広がり関数(インパルス応答)、t
1 (x,y,z)は組織圧縮前の散乱体関数、t
2 (x,y,z)は組織圧縮後の散乱体関数である。ただし、散乱体関数とは組織モデル内の散乱体が存在する位置ではその散乱係数の値をとり、その他の位置では0であるような関数である。また、組織圧縮後の散乱体関数t
2 (x,y,z)は組織圧縮前散乱体関数t
1 (x,y,z)の各散乱体の位置を組織モデルの変形に応じて移動させたものである。ただし、組織圧縮に伴う各散乱***置での変位は有限要素解析により得られる要素節点での変位ベクトルを線形補間することにより求めている。
また、この実施の形態ではシミュレーション超音波システムとして無焦点、かつ減衰のないシステムを仮定する。すなわち、超音波システムの点広がり関数h(x,y,z)は空間的に不変であると仮定する。さらに、点広がり関数は次式のように方向ごとに分離できると仮定する。
ここで、hy(y)は超音波ビーム方向の点広がり関数、hx(x),hz(z)はそれぞれ超音波ビームに直交した方向の点広がり関数である。ただし、x方向は超音波断層面内の方向(横方向)、z方向は超音波断層面に垂直な方向(スライス方向)とする。そして、各方向の点広がり関数は実際の超音波装置により計測したワイヤー・ターゲット(水中に張った直径0.13mmのワイヤー)からの反射エコー分布を基に作成する。図16は、超音波中心周波数5.0MHzの場合に用いた各点広がり関数の一例を示す図である。図16(A)は軸方向の点広がり関数hy(y)を示し、これはガウス関数に正弦波をかけたものによって実際のワイヤー・ターゲットからの反射エコー分布を近似し、図16(B)は横方向点の広がり関数hx(x)を、図16(C)はスライス方向の広がり関数hz(x)をそれぞれ示し、これらはガウス関数によって実際のワイヤー・ターゲットからの反射エコー分布を近似する。また、各関数のパラメータは中心周波数に応じて変えており、各シミュレーションの際に改めて説明する。
次に、変位(歪み)分布推定法として本実施の形態の拡張複合自己相関法の有効性をシミュレーションにより評価する。まず、拡張複合自己相関法の複合自己相関法に対する拡張点である組織の横方向変位に対応できる点について検証する。
図17は、組織モデルの概略を示す図である。組織モデルは、外形60mm×60mm(2次元)で、一様な弾性係数分布をもつモデルである。そして、この組織モデルを軸方向に一様な3%の歪みが生じるように圧縮する。ここで、このシミュレーションに関しては拡張複合自己相関法のみの評価を行うため、組織モデルとしては単純な1次元弾性体を仮定している。そして、組織の横方向移動(超音波プローブに対する相対的な横方向移動)に関する影響を検証するため、軸方向の圧縮と同時に横方向に0.0mmから1.4mmまでの横方向変位を与えた。ただし、横方向に関しては単純な平行移動であり、組織に対して超音波プローブが完全に滑った場合を想定している。
そして、この組織モデルに対して変形前後のシミュレーションRF信号を作成する。その際用いた超音波システムに関するパラメータは、中心周波数5.0MHz、パルス幅0.5mm、超音波ビーム幅1.0mm、走査ライン間隔0.5mm、サンプリング周波数30MHzである。そして、この圧縮前後のシミュレーションRF信号を用いて、本発明に係る拡張複合自己相関法により歪み分布を推定する。また、比較のために複合自己相関法と空間相関法を用いても歪み分布推定を行った。ここで、単純に精度等を比較できるように各手法では同じサイズの相関窓と探索範囲を用いた。具体的には、拡張複合自己相関法と空間相関法では1.6mm(軸方向)×2.5mm(横方向)の2次元相関窓と5.6mm(軸方向)×7.5mm(横方向)の2次元探索範囲を用い、1次元処理の複合自己相関法では軸方向だけ同じ1.6mmの1次元相関窓と5.6mmの1次元探索範囲を用いた。
このようにして各手法により歪み分布を推定した結果、図18〜図20のようになる。ここで、図18は横方向変位に対する各手法の歪み推定誤差を示している。◇は複合自己相関法、□は拡張複合自己相関法、△は空間相関法を示す。ただし、歪み推定誤差eとしては次式を用いた。
ここで、ε^
iは推定された歪み、ε
i は実際の歪み(理想値)、iは要素番号、Nは要素数である。また、図19は横方向変位が0.0mmの場合の各手法(複合自己相関法、拡張複合自己相関法、空間相関法)により推定した歪み分布、図20は横方向変位が0.4mmの場合の各手法(複合自己相関法、拡張複合自己相関法、空間相関法)により推定した歪み分布を示す図である。なお、図19及び図20は軸方向の深さごとに推定された歪みの平均値と標準偏差を表している。
このシミュレーション結果より、1次元の複合自己相関法(CA法)では組織の相対的な横方向変位が超音波ビームを越えて生じてしまう(今回の場合、横方向変位がビーム幅の片側分である0.5mmを超えてしまう)と歪み推定が急に悪くなってしまうのに対し、2次元の拡張複合自己相関法では横方向変位の大きさにかかわらず安定して歪み推定が可能であることが理解できる。また、空間相関法も横方向変位にかかわらず安定して歪み推定が行えるものの、拡張複合自己相関法と比べると2倍以上誤差が大きく精度が悪いことが理解できる。また、各手法の処理時間を比較したところ、下表のように拡張複合自己相関法は複合自己相関法に比べて3.6倍時間がかかってしまうものの、空間相関法と比べると1/(7.7)の時間しかかからなかった。そのため、拡張複合自己相関法はある程度リアルタイム性が保たれていることが理解できる。
┌−−−−−−−−−−−┬−−−−−−−−−−−┬−−−−−−−−─┐
│ 手 法 │ 処理時間 │ 処理時間比較 │
├−−−−−−−−−−−┼−−−−−−−−−−−┼−−−−−−−−─┤
│ 複合自己相関法 │ 26秒 │ 1/(3.6) │
├−−−−−−−−−−−┼−−−−−−−−−−−┼−−−−−−−−─┤
│ 拡張複合自己相関法 │ 1分34秒 │ 1.0 │
├−−−−−−−−−−−┼−−−−−−−−−−−┼−−−−−−−−─┤
│ 空間相関法 │ 12分5秒 │ 7.7 │
└−−−−−−−−−−−┴−−−−−−−−−−−┴−−−−−−−−─┘
次に、斜め方向圧縮に関する検証について説明する。前述のシミュレーションでは簡単な2次元の均一組織モデルを用いたが、次は実際の生体組織と同じ3次元構造を持つ組織モデルを用いてシミュレーションを行う。また、超音波プローブにより組織を圧縮する際、超音波ビーム方向(軸方向)に圧縮するのが理想であるが、ここでは斜めに圧縮してしまった場合の影響を検証する。
図21は、斜め方向圧縮の検証に使用される組織モデルの概略を示す図である。組織モデルは、図21(A)に示すように、外形が60mm×60mm×60mmの3次元構造をしており、この組織モデル中に直径15mm、長さ60mmの硬い円柱形内包物が存在するようなモデルである。ここで、周辺の弾性係数(ヤング率)は10kPa、内包物の弾性係数は周辺より3倍硬い30kPaとする。ただし、この弾性係数の値は今回主な対象としている***組織の弾性係数および乳がんの弾性係数を基に設定する。そして、この組織モデルを2通りの方法で圧縮を行った。1つ目は、図21(B)に示すように、この組織モデルに対して上面から軸方向に一様な200Paの外部圧力を加えて、組織モデルを軸方向に約2%圧縮する。2つ目は、図21(C)に示すように、この組織モデルに対してモデル上面から斜め方向の一様な外部圧力(軸方向に200Pa、横方向に30Pa)を加えて、組織モデルを斜め方向に圧縮する。
そして、上記の2通りの場合についてそれぞれ組織圧縮前後のシミュレーションRF信号を作成し、拡張複合自己相関法により歪み分布を推定する。なお、ここでも比較のために複合自己相関法と空間相関法による歪み分布推定も行う。ただし、単純に比較できるように各手法で用いた相関窓サイズと探索範囲は同じとし、そのサイズは前のシミュレーションと同じとする。また、シミュレーションRF信号作成の際に用いた超音波システムに関するパラメータも前のシミュレーションと同じ、中心周波数5.0MHz、パルス幅0.5mm、横方向超音波ビーム幅1.0mm、スライス方向超音波ビーム幅2.0mm、走査ライン間隔0.5mm、サンプリング周波数30MHzとする。
上記のようにして行ったシミュレーションの結果は、図22及び図23に示すようになる。ここで、図22は組織モデルを単純に軸方向に圧縮した場合の歪み分布推定結果を示し、図23は組織モデルを斜め方向に圧縮した場合の歪み分布推定結果を示す。なお、各図における理想的な歪み分布とは、各条件で3次元有限要素解析を行って得られた軸方向歪み分布であり、この歪み分布を正解としている。また、図22及び図23における結果は組織モデルの中央断面での結果である。ここで、図23において理想歪み分布が左右対称でないのは斜め方向に圧縮した影響であり、特に今回の場合は図に向かって右斜め下に圧縮したため、図右上の部分の横方向変位が大きくなっている。
そして、このシミュレーションの結果より、軸方向に圧縮した場合は拡張複合自己相関法と複合自己相関法はほぼ同じ結果となったが、斜め方向に圧縮した場合は横方向変位が大きくなるため、複合自己相関法では推定できなくなる領域が生じてしまうのに対し、拡張複合自己相関法では先のシミュレーション同様に横方向変位の大きさにかかわらず安定して歪み推定が可能であることが改めて確認された。また、空間相関法も横方向変位の大きさには依存しないものの拡張複合自己相関法の結果と比較すると明らかに推定精度が悪いのが見てとれる。そのため、前のシミュレーションと合わせて改めて拡張複合自己相関の有効性が確認された。
前述の拡張複合自己相関法は、精度良く、かつ高速に組織歪み分布を推定できることができる。そこで、次に組織弾性映像システムの第2段階である歪み分布から弾性係数分布を推定する手法で、実施の形態2の3次元有限要素モデルによる弾性係数分布再構成法についてシミュレーションにより検証を行う。
実施の形態2の弾性係数分布再構成法の最大の特徴は3次元の力学的つりあい方程式に基づいて弾性係数分布を推定することである。そこで、実施の形態2の手法と手順的には同じであるが、2次元の力学的つりあい方程式を基にしている点だけが異なる実施形態1の2次元再構成法を用いて、実施の形態2の3次元再構成法と比較することにより検証を行う。実施形態1の2次元再構成法では、2次元の平面歪み状態を仮定した弾性方程式および有限要素法を用いて弾性係数分布を推定する。
そこで、まず組織モデルとして、実際の生体組織と同じ3次元構造を持つ図24のような2つのモデルを用いる。図24は、3次元構造を持つ2つの組織モデルの一例を示す図である。図24(a)の内包物モデルは、乳がんを模擬したモデルで、外形100mm×100mm×100mmのモデル中に直径20mmの硬い内包物が存在するもので、内包物の弾性係数は30kPa、周辺の弾性係数は10kPaとする。これらの弾性係数の値は、先ほどのシミュレーションと同じように実際の***組織の弾性係数を基に決めている。また、周辺と内包物のポアソン比としては、共に非圧縮性に近いため0.49で一様とする。そして、図24(b)の層状モデルは筋肉などの層状のものを模擬したモデルで、外形100mm×100mm×100mmのモデル中に厚さ20mmの硬い層がモデル中に存在するというもので、硬い層の弾性係数は30kPa、周辺の弾性係数は10kPaとする。そして、このモデルの場合もポアソン比は049で一様とする。
そして、図24(a)の内包物モデルの場合はモデル上部から軸方向に100Paの一様な外部圧力により、図24(b)の層状モデルの場合はモデル上部から軸方向に150Paの一様な外部圧力により、各モデルをコンピュータ上で圧縮する。ここで、2つのモデルで外部圧力の強さを変えているのは、各モデルとも同じ約1%の歪みが生じるようにするためである。そして、これらの組織モデル圧縮前後のシミュレーションRF信号を作成し、拡張複合自己相関法により軸方向歪み分布を推定する。そして、この推定された軸方向歪み分布と圧縮の際の境界条件とから、提案した3次元弾性係数分布再構成法により弾性係数分布を推定する。また、同じ軸方向歪み分布と境界条件を用いて比較のために2次元再構成法によっても弾性係数分布を推定する。ここで、シミュレーションRF信号を作成するために用いた超音波システムのパラメータとしては、中心周波数3.75MHz、パルス幅0.75mm、横方向超音波ビーム幅2.0mm、スライス方向超音波ビーム幅2.0mm、走査ライン間隔2.0mmである。また、拡張複合自己相関法における相関窓のサイズは3.2mm(軸方向)×4.0mm(横方向)、探索範囲は11.2mm(軸方向)×14.0mm(横方向)とする。さらに、3次元有限要素モデルよる弾性係数分布再構成では2mm(軸方向)×2mm(横方向)×5mm(スライス方向)の直方体要素50000個を用いて組織モデルを構成する。
そして、このシミュレーションの結果は、図25〜図28に示すようになる。図25及び図26は内包モデルにおける各推定結果を示す、図27及び図28は層状モデルにおける各推定結果を示す。ただし、3次元再構成法では弾性係数の3次元分布を推定しているが、ここではモデル中央断面での結果のみを示す。これは、2次元再構成法では2次元断面での弾性係数分布しか推定できないので、ここでは比較できる中央断面のみ示す。また、各組織モデルにおける3次元再構成結果と2次元再構成結果を数値的に評価したところ次のような結果が得られた。
┌−−−−−−−−−−┬−−−−−−−┬−−−−−−┬−−−−−−−┐
│ │ 周辺領域に │周辺領域に │モデル中心に │
│ │おける弾性係数│おける標準 │おけるコントラ│
│ │ 誤差[%] │偏差[%] │スト誤差[%]│
├−−┬−−−−−−−┼−−−−−−−┼−−−−−−┼−−−−−−−┤
│内包│3次元再構成法│ 3.5 │ 15.5│ 11.0 │
│物モ├−−−−−−−┼−−−−−−−┼−−−−−−┼−−−−−−−┤
│デル│3次元再構成法│ 30.9 │ 17.9│ 35.9 │
├−−┼−−−−−−−┼−−−−−−−┼−−−−−−┼−−−−−−−┤
│層状│3次元再構成法│ 8.5 │ 26.8│ 3.1 │
│モデ├−−−−−−−┼−−−−−−−┼−−−−−−┼−−−−−−−┤
│ル │3次元再構成法│ 24.9 │ 22.1│ 43.5 │
└−−┴−−−−−−−┴−−−−−−−┴−−−−−−┴−−−−−−−┘
ここで、用いた評価用のパラメータは、周辺領域における弾性係数誤差e
s 、周辺領域における標準偏差SD
s 、モデル中心におけるコントラスト誤差e
c であり、それぞれ次式のように定義する。
ただし、上式におけるΣは周辺領域における総和、E^ は推定された弾性係数、Eは実際の弾性係数、N
s は周辺領域の要素数、E~
sは周辺領域における弾性係数の平均値、E^
cはモデル中心における推定弾性係数、E
c はモデル中心における実際の弾性係数、E
s は周辺領域における実際の弾性係数である。
以上のシミュレーション結果より、実施形態2の3次元有限要素モデルによる弾性係数分布再構成法の方が平面応力状態を仮定した2次元再構成法よりもより正確な弾性係数を推定できることが確認された。ここで、3次元再構成法では弾性係数の値が正確に推定されているのに対し、2次元再構成法では実際の弾性係数の値よりも小さく推定しまっている。これは、前もって予想されていたように2次元状態では2次元考察面に垂直な方向の動きが制限されてしまうためである。そのため、実際の生体組織と同じ3次元を基にした弾性係数分布再構成法が必要であることがこのシミュレーションにより改めて示された。
次に、上述の拡張複合自己相関法と3次元有限要素モデルによる弾性係数再構成法を実装した超音波診断システムの具体的構成について説明する。図29は、この超音波診断システムの基本構成を示す図である。この超音波診断システムは、3次元超音波スキャナ281、超音波診断装置282、パーソナルコンピュータ283、パルスモータコントローラ284、パルスモータ285、圧力計286などから構成される。3次元超音波スキャナ281は超音波パルスを組織内に送信し、かつ組織からの超音波エコー信号を受信するためのものである。ただし、ここでは3次元有限要素モデルによる弾性係数再構成法を用いるため、組織内の3次元的なデータが必要になる。そこで、この超音波診断システムでは3次元超音波スキャナ281は3次元走査が可能な構成となっている。超音波診断装置282は3次元超音波スキャナ281を制御したり、リアルタイムに超音波Bモード画像を表示して計測部位を決定したりするためのもので、従来の超音波診断装置をそのまま用いることができる。この超音波診断装置はフルデジタル化された装置で内部にフレームメモリを持っているため、計測したRF信号を一時的に保存しておくことが可能となっている。パーソナルコンピュータ283は、超音波診断装置282によって計測されたRF信号を受け取り、組織弾性特性推定のための処理(前述の提案手法の処理)、および処理結果の表示を行うためのものである。パルスモータ285は組織圧縮を制御するためのものであり、位置固定が可能なアームの先端に固定されており、かつパルスモータ285の可動部分には3次元超音波スキャナ281が取り付けてある。そして、パルスモータコントローラ284によりこのパルスモータ285を制御し、超音波スキャナ281を組織表面で上下に動かすことにより数パーセントの微小組織圧縮を精度良く行う。圧力計286は弾性係数再構成の際に必要な境界条件である体表上での圧力を測るためのもので、超音波スキャナ281と体表との間に置かれる。ただし、ここでは超音波スキャナ281で組織圧縮を行った際の体表上での圧力分布は一様であると仮定し、圧力計286で計測された値をその圧力値として用いる。
図30は、この超音波診断システムで用いた超音波スキャナ281の具体的構成を示す図である。3次元超音波スキャナ281は、超音波振動子が2次元平面状に並んだ2次元アレイ型ではなく、コンベックス型の2次元走査プローブを水中で機械的にスライス方向に振ることにより3次元走査を行うタイプのものである。
図29の超音波診断システムは、乳がん診断を主な対象としているため、超音波スキャナの特性も乳腺用に設定してある。今回用いた超音波スキャナの主な特性としては、超音波中心周波数7.5MHz、サンプリング周波数30MHz、走査ライン数71本、フレーム数44枚、振動子の振れ角30°、1回の3次元走査にかかる時間0.5秒となっている。ここで、振動子の振れ角とはコンベックス型のプローブをスライス方向に振る際の振れ角のことであり、フレーム数とはコンベックス型のプローブを1回振る際に取得される走査面(フレーム)の数である。また、水中のワイヤー・ターゲットにより超音波パルスの特性を計測したところ、パルス幅約0.5mm、横方向ビーム幅約1.5mm、スライス方向ビーム幅約2.6mmであった。
図29の超音波診断システムの弾性計測の動作例について説明する。まず、超音波診断装置282のリアルタイムBモード像を見ながら、アームにつながった3次元超音波スキャナ281を動かし、計測を行いたい部位に超音波スキャナ281を設定する。この際、超音波スキャナ281は3次元走査を行わず(すなわち、コンベックス型のプローブを機械的に振らず)、スキャナの中央面のBモード像のみを超音波診断装置282に表示する。これは、実施の形態の超音波診断装置282で3次元走査をするとリアルタイムにBモードを表示できないためである。従って、リアルタイムにBモードを表示することのできる超音波診断装置の場合は、3次元走査を行いながら部位の設定を行なうことができる。計測部位に超音波スキャナ281を移動させたら、アームの可動部を固定し、超音波スキャナ281を固定する。そして、組織圧縮前の3次元RF信号を計測する。これは、3次元走査用のボタンを押すことにより自動的に3次元走査される。また、1回の3次元走査にかかる時間はわずか0.5秒である。このとき、計測された圧縮前のRFデータは超音波装置内のフレームメモリに保存しておく。次に、パルスモータ・コントローラ284の圧縮用のボタンを1回押すことにより、アームに取り付けられたパルスモータ285が前もって設定しておいた量だけ超音波スキャナ281を動かし、超音波スキャナ281自身により組織を圧縮する。そして、パルスモータ285が止まって組織圧縮を行っている状態で、再び3次元走査用
のボタンを押し、組織圧縮後のRFデータを計測する。ここで、組織圧縮後のRFデータも圧縮前のRFデータと同様に超音波装置282内のフレームメモリに保存される。また、圧縮している状態で超音波スキャナ281の端に取り付けられた圧力計の圧力を計測しておく。以上で、計測部は終わりで、組織圧縮を解除し、被験者は解放される。
被検者解放後は、パーソナルコンピュータ283から超音波診断装置282内のフレームメモリにアクセスし、組織圧縮前後のRFデータをパーソナルコンピュータ283上のハードディスクに保存する。これは、超音波装置内のフレームメモリは一時的なものであり、1回の計測分しか容量がないためである。パーソナルコンピュータ283は、拡張複合自己相関法と3次元有限要素モデルによる弾性係数分布再構成法のプログラムを実行し、組織圧縮前後のRFデータから歪み分布及び弾性係数分布を推定する。そして、パーソナルコンピュータ283は、表示用のプログラムによってモニタ上にBモード像、歪み像、弾性係数像を並べて表示する。
以上説明したように、本発明の各実施形態の超音波診断システム、歪み分布表示方法及び弾性係数分布表示方法によれば、横方向変位に対応して変位分布を推定することができると共に超音波ビーム方向(軸方向)の歪み分布のみから弾性係数分布を再構成することができる。