JP4258015B2 - 超音波診断システム、歪み分布表示方法及び弾性係数分布表示方法 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は、超音波診断装置を用いて、生体組織の硬さを定量的に計測することのできる超音波診断システム、歪み分布表示方法及び弾性係数分布表示方法に関する。
【０００２】
【従来の技術】
超音波の医療面への応用もエレクトロニクス技術の進歩と相まってさまざまな臨床領域へと広がっている。その例としては、超音波を生体情報の取得手段として利用する超音波断層像（Ｂモード像）やドプラ血流計測、そして超音波のエネルギーを直接利用する超音波ハイパーサーミア（温熱治療）や体外衝撃波結石破砕装置などがある。これらの中でも特に超音波Ｂモード像は計測のリアルタイム性、手軽さ、安全性のため臨床の場で幅広く利用されている。ここで、超音波Ｂモード像とは体内に超音波を放射して音響インピーダンスが異なる組織境界での反射エコーを輝度変調しながら、これを２次元断面的に走査することによって組織の形状を画像化したものである。
【０００３】
これに対し、組織形状だけではなく組織内の音速や減衰定数などの物理量を超音波により計測・画像化し診断に利用しようとするウルトラソニック ティシュー キャラクタライゼイション（ｕｌｔｒａｓｏｎｉｃ ｔｉｓｓｕｅ ｃｈａｒａｃｔｅｒｉｚａｔｉｏｎ）と呼ばれる分野がある。これは、組織の物理量を計測して組織診断に利用とするものである。そして、その中の１つとして組織の硬さ、すなわち弾性特性を計測しようとする分野があり、現在盛んに研究されている。これは、組織の弾性特性がその病理状態と深く関連しているためである。例えば、乳がんや甲状腺がんなどの硬化性がんや肝硬変、動脈硬化症などは正常組織よりも病変部分が硬くなることが知られている。そして、これまでこれらの硬さ情報は触診により得ていた。しかし、触診では客観的な情報表現が難しく、医師の経験も必要で、また計測できる領域も体表付近のある程度大きな病変に限られる。
【０００４】
そこで、超音波やＭＲＩを利用して組織の弾性特性を定量的に計測・画像化しようとする研究が行われるようになった。まず、体表から機械的振動を与えその横波の伝播速度を超音波により計測し、横波の伝播速度から組織の硬さを評価する試みが行われた（Ｒ．Ｍ．Ｌｅｒｎｅｒ，Ｓ．Ｒ．Ｈｕａｎｇ ａｎｄ Ｋ．Ｊ．Ｐａｒｋｅｒ，”Ｓｏｎｏｅｌａｓｔｉｃｉｔｙ ｉｍａｇｅｓ ｄｅｒｉｖｅｄ ｆｒｏｍ ｕｌｔｒａｓｏｕｎｄ ｓｉｇｎａｌｓ ｉｎ ｍｅｃｈａｎｉｃａｌｌｙ ｖｉｂｒａｔｅｄ ｔｉｓｓｕｅｓ”， Ｕｌｔｒａｓｏｕｎｄ Ｍｅｄ．Ｂｉｏｌ．，ｖｏｌ．１６，ｎｏ．３，ｐｐ．２３１−２３９，１９９０）。これを第１の従来技術とする。この第１の従来技術は、硬い組織では横波の伝播速度が速く、軟らかい組織では横波の伝播速度が遅いということを基にしている。しかし、この方法は分解能が低いという問題点があった。
【０００５】
これに対し、体表から静的な圧力を加えて組織をわずかに圧縮変形させ、その際生じる組織内部の歪みを超音波により計測し、歪みから組織の弾性特性を評価する方法が１９９０年頃から始まった（Ｊ．Ｏｐｈｉｒ， Ｉ．Ｃｅｓｐｅｄｅｓ，Ｈ．Ｐｏｎｎｅｋａｎｔｉ，Ｙ．Ｙａｚｄｉ ａｎｄ Ｘ．Ｌｉ，”Ｅｌａｓｔｏｇｒａｐｈｙ： Ａ ｑｕａｎｔｉｔａｔｉｖｅ ｍｅｔｈｏｄ ｆｏｒ ｉｍａｇｉｎｇ ｔｈｅ ｅｌａｓｔｉｃｉｔｙ ｏｆ ｂｉｏｌｏｇｉｃａｌ ｔｉｓｓｕｅｓ”，Ｕｌｔｒａｓｏｎｉｃ Ｉｍａｇｉｎｇ， ｖｏｌ．１３，ｐｐ．１１１−１３４，１９９１）これを第２の従来技術とする。この第２の従来技術は、硬い組織では生じる歪みが小さく、軟らかい組織では歪みが大きくなることを基にしている。
【０００６】
そして、その後ＭＲＩを用いて同様の原理により組織歪みから弾性特性を評価する方法が研究されるようになった。しかし、ＭＲＩを利用した方法はその特性上リアルタイム計測が困難であり、また体表から組織変形を加えることが難しいという問題点がある。従って、現在ではリアルタイム計測が可能で軽便な超音波を利用して、静的組織圧縮下における歪み推定原理に基づいた組織弾性特性評価を採用するようになってきた。
【０００７】
一般に、超音波とは「人間の可聴音域（約２０Ｈｚ〜２０ｋＨｚ）より周波数の高い音」と定義されている。しかし、使う場面によっては人間の耳に聞こえる音も超音波と呼ばれることがある。そこで、最近では「超音波とは聞くことを目
的としない音」と定義されるようになってきた。ただし、超音波診断装置で用いられる超音波の周波数は１ＭＨｚ〜１０ＭＨｚが主流である。現在、超音波は医学をはじめとして様々な分野で利用されているが、特に生体計測の分野では以下のような性質のため超音波が広く用いられ、超音波診断装置として利用されている。
・超音波は生体を媒質として伝播できる。
・超音波が生体中を進む速度（約１５００ｍ／ｓ）は光（電磁波）に比べ桁違いに遅い。
・超音波には指向性があるため、音のビームとして利用できる。
・弱いパワーであれば生体に対して無侵襲である。
・生体の組織によって音響特性が異なるため、組織の境界で反射エコーが得られる。
【０００８】
図１は、超音波診断装置の原理を説明するための図である。図から明らかなように、超音波プローブ１０は電気信号を超音波に、また超音波を電気信号に変換するものであり、この超音波プローブ１０を用いて生体組織１１内に超音波パルスを放射する。生体組織１１内に放射された超音波パルスは音響インピーダンスの異なる第１の境界１２で一部が反射され、反射エコー１２ａとして超音波プローブ１０側に向かい、その残りは透過していく。透過した超音波パルスは次の音響インピーダンスの異なる第２の境界１３で同様に一部が反射され、反射エコー１３ａとして超音波プローブ１０側に向かい、その残りは透過する。このようにして反射した反射エコー（超音波エコー信号）は超音波プローブ１０によって受波され電気信号に変換される。受波された反射エコー信号は、受信エコー信号のようになっている。このとき、超音波プローブ１０から超音波パルスが放射されてから距離Ｌの位置にある反射物体１４（音響インピーダンスの異なる境界）からのエコー信号を受信するまでの時間ｔは、
【数０１】

となる。ここで、ｃは生体内での音速であり、軟組織では１５００［ｍ／秒］にほぼ一定とみなせる。よって、超音波エコー信号を受信するまでの時間ｔを計測すればプローブから反射物体までの距離Ｌを求めることができる。
【０００９】
そして、電気信号に変換された超音波信号（受信エコー信号）をディスプレイに表示する方法としては、図２に示すような３種類の方法がある。図２（Ａ）は、Ａモード方式であり、表示用ディスプレイの横軸にプローブからの距離、縦軸に受信した反射エコーの強度（振幅）をとり反射エコー信号をグラフ状に表示するものである。図２（Ｂ）は、Ｂモード方式であり、超音波プローブを２次元断層的に走査したときに得られる反射エコー信号の強度を輝度変調し、走査位置に応じてディスプレイに表示するものである。この方式を用いると生体内の断層像が得られるため、今日最も広く用いられている。また、このとき得られる断層像をＢモード像という。図２（Ｃ）は、Ｍモード方式であり、対象となる物体が運動している場合、超音波プローブ１０の位置を固定しても時々刻々異なったＡモード波形が観測される。このＡモード波形を輝度変調してディスプレイの縦方向に表示し、さらに、時間に伴って横方向に走査する方式をＭモードという。この方式を用いると組織の動く様子が画像化されるため、心臓の弁や壁の運動を調べるのに利用されている。
【００１０】
図３は、超音波プローブの種類を示す図である。Ｂモード像の走査方式・走査形状の違いにより現在、様々な超音波プローブが利用されている。まず、超音波ビームの走査方式の種類としては以下に示す３通りの方式がある。
【００１１】
第１は、手動走査方式である。これは、振動子（圧電素子）を先端に１つだけ装着したプローブを手で体表に沿わせて走査し、そのプローブの位置や角度をアームの検出機構により検出して、プローブの動きに対応した画像を表示する方式である。第２は、機械走査方式である。これは、振動子を先端に１つだけ装着したプローブをモーター等により動かし、そのプローブの位置や角度を検出機構により検出して、その動きに応じた画像を表示する方式である。第３は、電子走査方式である。これは、短冊状の振動子を先端に多数装着したプローブを用い、駆動する振動子を電子スイッチ等により制御し、走査を行う方式である。これらの走査方式の中で、現在広く用いられている方式は、電子走査方式であり、機械走査方式は一部の特殊な用途に用いられているのみである。
【００１２】
次に、電子走査方式のプローブでも走査形状の違いにより以下のように分けられている。第１は、セクタ走査方式である。この方式は、図３（Ａ）に示すように、超音波ビームを扇状に走査するもので、このような走査を行うプローブをセクタスキャンプローブ（セクタフェイズドアレイプローブ）という。浅部の視野は狭いが、深部では広い範囲を観測することが可能であるため、肋骨やガス像の合間からの観察に優れている。第２は、リニア走査方式である。この方式は、図３（Ｂ）に示すように、超音波ビームを直線状に走査するもので、このような走査を行うプローブをリニアスキャンプローブ（リニアアレイプローブ）という。浅部で広い視野が得られるため、腹部検査で用いられている。第３は、オフセットセクタ走査方式である。この方式は、図３（Ｃ）に示すように、超音波ビームを扇状に走査するが、要の部分を表示しないもので、このような走査を行うプローブをコンベックススキャンプローブ（コンベックスアレイプローブ）という。浅部から深部まで広い範囲を観測できるため、腹部検査で広く用いられている。このような走査形状を持った電子走査方式の超音波プローブが現在、主に用いられている。その他、特殊なものとしては血管内部から血管周辺を観察するためのカテーテルプローブや超音波顕微鏡用の超高周波超音波プローブなどもある。また、最近では３次元の超音波像を得るための２次元アレイプローブの開発も行われている。
【００１３】
図４は、超音波診断装置を用いて、組織の硬さに関する情報（組織の弾性特性）を計測する手法（機械的振動下における横波伝播速度からの弾性特性評価）の一例を示す図である。これは、前述の第１の従来技術に相当するものであり、超音波を用いて組織の硬さに関する情報を計測する方式であり、組織に機械的振動を与えてその横波の伝播速度から硬さ情報を評価する方式である。この方式は、硬い組織では横波の伝播速度が速く、軟らかい組織では横波の伝播速度が遅いことを基にしている。ただし、厳密には生体組織中を伝わる横波の伝播速度は次式のように組織の密度、せん断弾性係数、せん断粘性係数、および振動の周波数に関係している。
【数０２】

ここで、ｖは横波の伝播速度、μ₁ はせん断弾性係数、μ₂ はせん断粘性係数、ρは組織の密度、ωは機械振動の角周波数である。
【００１４】
この方式では、まず低周波（数百ヘルツ）で振動する低周波振動子４１を生体組織１１の体表に接触させ、組織内部に振動を伝播させる。この振動により誘起された横波の振幅と位相の分布を血流計測に用いられるドプラ法を用いて計測する。そして、横波の振幅と位相の分布から組織の弾性特性（横波の伝播速度）を推定することになる。ただし、その際、組織の粘弾性特性は無視し、また組織の密度は一様であると仮定する。このように仮定すると組織のせん断弾性係数μ₁ は、μ₁＝ρｖ²のように横波の伝播速度の２乗に比例する。
【００１５】
しかし、組織の粘弾性特性を無視することは難しく、組織の密度も生体内で変化するため、この方法により組織の弾性特性を定量的に評価することは難しい。また、横波の伝播速度分布も機械振動の波長程度の分解能でしか得られない。
【００１６】
そこで、機械的振動を与えて組織の弾性特性を評価するものに対して、前述の第２の従来技術のように、組織を静的に圧縮してその際生じる組織内の歪み分布から弾性特性を評価する方式が提案されている。これは、硬い組織では歪みが小さく、軟らかい組織では歪みが大きくなることに基づいている。
【００１７】
図５（Ａ）は、静的圧縮による組織弾性計測方式の具体例を示す図である。図５（Ｂ）は、静的圧縮による組織弾性計測方式の原理を示す図である。図から明らかなように、この方式は、従来の超音波診断装置および超音波プローブ１０をそのまま用いる。まず、超音波プローブ１０によって組織１１の圧縮前の超音波エコー信号（圧縮前ＲＦ信号）を計測する。その後、超音波プローブ１０自身により組織１１をわずかに（数パーセント程度）圧縮し、組織１１の圧縮後の超音波エコー信号（圧縮後ＲＦ信号）を計測する。そして、計測された組織圧縮前後のＲＦ信号から圧縮によって組織内部の各点がどれだけ動いたかという移動量である変位分布を推定する。この変位分布推定手法の主なものとしては、空間相関を用いるものとドプラの原理を用いるものがある。
【００１８】
図６は、空間相関法の原理を示す図である。この方法は、圧縮によって生じた組織内部の変位分布を組織圧縮前後のＲＦ信号（またはＲＦ信号の包絡線）から２次元相関関数を用いたテンプレートマッチングにより推定する手法である。その具体的な処理は以下のようになる。まず、組織圧縮前後のＲＦ信号（またはその包絡線信号）をｉ₁（ｔ，ｘ），ｉ₂（ｔ，ｘ）とすると、この２つの信号の相互相関係数Ｃ（ｔ，ｘ；ｎ，ｍ）は、
【数０３】

となる。ここで、ｔは超音波ビーム方向（軸方向）の座標、ｘはそれに直交する方向（横方向）の座標、ｔ₀ は軸方向の相関窓サイズ、ｘ₀ は横方向の相関窓サイズ、Ｌ_t は軸方向のサンプリング間隔、Ｌ_x は横方向のサンプリング間隔、ｎ，ｍは整数である。そして、この相互相関関数が最大となるときの（ｎ，ｍ）を（ｋ，ｌ）とすると、計測点（ｔ，ｘ）における軸方向の変位ｕ_y と横方向の変位ｕ_x はそれぞれ次式のようにして求められる。
ｕ_y ＝ｋＬ_t
ｕ_x ＝ｌＬ_x
ただし、横方向のサンプリング間隔Ｌ_x
は軸方向のサンプリング間隔Ｌ_t よりも大きいので、推定される変位成分の精度は横方向成分の方が軸方向成分よりも悪くなる。上記の処理を各計測点について行い変位分布推定する手法が空間相関法である。そのため、空間相関法では２次元の変位ベクトル成分を推定できるという特徴がある。また、組織が大変形（５％程度）した場合でも変位分布を推定できる。しかし、計算量が膨大になるため超音波計測の利点であるリアルタイム性を損なってしまう。また、変位推定精度もサンプリング間隔により制限されてしまうため、後に述べるドプラ法と比べると精度が悪いという問題点もある。
【００１９】
図７は、ドプラ法の原理を示す図である。この方法は、圧縮によって生じた組織内部の変位分布を組織圧縮前後のＲＦ信号から血流計測に用いられているドプラの原理を利用して推定する手法である。その具体的な処理は以下のようになる。まず、組織圧縮前後のＲＦ信号を次式のようにモデル化する。
【数０４】

ここで、ｉ₁ （ｔ）は圧縮前のＲＦ信号、ｉ₂ （ｔ）は圧縮後のＲＦ信号、Ａ（ｔ）は包絡線、ω₀ は超音波の中心角周波数、τは時間シフトである。そして、この２つのＲＦ信号をそれぞれ直交検波すると、次式のようなベースバンド信号が得られる。
【数０５】

そして、この２つの信号の自己相関関数Ｒ₁₂（ｔ）（本来は相互相関関数であるが共に同じ部位からの信号であるためドプラ計測では自己相関関数と呼ぶ）は次式で表される。
【数０６】

ここで、Ｒ_A （ｔ）は包絡線の自己相関関数、ｔ₀ は相関窓サイズである。また、＊は複素共役を表している。よって、この自己相関関数Ｒ₁₂（ｔ）の位相φ（ｔ）から圧縮による時間シフトτ、軸方向変位ｕ_y
が次式のようにして求まる。
【数０７】

ただし、ｃは組織内の音速であり、生体内で一定と仮定する。
【００２０】
上記の処理を各計測点について行い変位分布を推定する手法がドプラ法であり、ドプラの原理を基にした血流計測と同じ処理となっている。そのため、リアルタイム計測が可能であるという利点がある。また、位相情報を用いているので変位推定精度が空間相関法よりも良い。しかし、組織内部の移動量が大きい（超音波中心周波数の４分の１波長以上となる）とエイリアシングを起こしてしまい正しい変位推定ができないという問題点がある。また、上式からもわかるように、ドプラ法は、軸方向の変位成分のみしか推定できないから、２次元の変位を推定できないという問題がある。
【００２１】
【発明が解決しようとする課題】
上述したように、空間相関法によると、２次元の変位を推定できるが、計算量が膨大になるため超音波計測の利点であるリアルタイム性を損なってしまうという問題がある。また、変位推定精度もサンプリング間隔により制限されてしまうため、ドプラ法と比べると精度が悪いという問題もある。
【００２２】
一方、ドプラ法によると、計算量が少ないのでリアルタイム計測が可能であり、また、位相情報を用いているので変位推定精度が空間相関法よりも良いが、２次元の変位を推定できないという問題がある。
【００２３】
本発明が解決しようとする課題は、生体組織の硬さを定量的に計測するにあたり、２次元変位分布の計算時間を短縮でき、かつ、変位分布の計算精度を向上させることにある。
【００２４】
【課題を解決するための手段】
上記課題を解決するため、本発明の超音波診断システムの第1の態様は、被検体組織との間で超音波信号を送受信する超音波探触子と、前記被検体組織に対する前記超音波探触子による圧縮前後に前記超音波探触子により受信された前記超音波エコー信号をそれぞれ直交検波して包絡線信号を作成する直交検波手段と、前記圧縮前後における前記被検体組織内の複数の計測点における前記包絡線信号間の相関係数を、各計測点を囲む２次元相関窓を用いて予め設定された２次元探索範囲内で、軸方向及びこの軸方向に直交する横方向にそれぞれ所定値ずつ離散した複数の格子点毎に計算し、前記相関係数が最大となる相関最大位置を求め、前記相関最大位置における前記圧縮前後の前記包絡線信号間の位相差を求める相関演算手段と、前記相関演算手段によって求められた前記相関最大位置及び前記位相差に基づいて前記圧縮に伴う前記被検体組織内の前記各計測点の前記軸方向及び前記横方向における変位を求める変位演算手段と、前記被検体組織内の前記軸方向及び前記横方向の前記各計測点における変位を空間微分することによって歪み分布情報を求める歪み演算手段と、前記歪み分布情報を表示する表示手段とを備えて構成される。特に、前記変位演算手段は、前記相関演算手段により求められた前記相関最大位置及び該相関最大位置における前記位相差のうち、前記横方向における前記相関最大位置に基づいて前記各計測点の前記横方向の変位を求め、前記相関演算手段により求められた前記軸方向の前記相関最大位置及び該相関最大位置における前記位相差の双方に基づいて前記各計測点の前記軸方向の変位を求めて前記歪み演算手段に出力することを特徴とする。
【００２５】
このように、本発明の相関演算手段は、直交検波手段から出力される包絡線信号を用いて被検体組織の圧縮前後の超音波エコー信号間で相関を計算する。このときに、サンプル間隔毎に相関計算を行うと計算量が膨大となるので、この発明では、相関演算手段は、各計測点を囲む２次元相関窓を用いて予め設定された２次元探索範囲内で、軸方向及びこの軸方向に直交する横方向にそれぞれ所定値ずつ離散した複数の格子点毎に計算する。これによって、計算量が大幅に減少し、２次元変位分布の計算時間を短縮でき、かつ、変位分布の計算精度を向上させることができるから、生体組織の硬さを定量的に計測でき、弾性情報のリアルタイム表示が可能となる。
【００２６】
この場合において、前記相関演算手段は、前記軸方向には２分の１波長間隔、前記横方向には前記超音波探触子のライン間隔の格子点でのみ相関係数の計算を行うようにすることができる。これは、相関演算手段が２次元相関窓内の軸方向及び横方向にそれぞれ計算する離散値を具体的に示したものである。この場合、さらに、前記相関演算手段は、前記直交検波手段によって直交検波された包絡線信号を用いて、前記軸方向に２分の１波長の整数倍だけシフトした直交検波包絡線信号を作成し、前記圧縮前後における前記包絡線信号間の相関係数を各格子点について計算するように構成できる。これは、直交検波された包絡線信号を用いて、各格子点における直交検波包絡線信号を作成するようにしたものであり、各格子点の相関係数の演算時に直交検波を繰り返す必要がなくなり、高速化及び回路の簡略化を図ることができる。
【００２７】
また、本発明の超音波診断システムの第２の態様は、第１の態様に加えて、前記被検体組織を有限個の要素に分割して少なくとも２次元の有限要素モデル化し、そのモデル化の情報と前記歪み分布情報を用いて弾性係数分布情報を演算する弾性係数情報演算手段と、前記弾性係数分布情報を表示する表示手段とを備えて構成することを特徴とする。これは、第１の態様の超音波診断システムによって得られた歪み分布情報を用いて弾性係数分布情報を演算するようにしたものである。組織を等方性弾性体と仮定するのは、外部から圧力を加えて組織を静的に圧縮した場合、応力と歪みの間の関係はほぼ線形であり、組織を弾性体として近似でき、被検体の組織はほぼ等方性が成り立つので、この発明では組織を等方性弾性体と仮定している。また、組織を近非圧縮性と仮定するのは、生体組織が非圧縮性（ポアソン比ν＝０．５）であると特殊な弾性方程式となり、有限要素法を適用することができなくなるからである。また、ポアソン比を生体内で一定とすることで弾性係数分布推定の推定パラメータをヤング率のみとすることができ、逆問題を簡単化できる。また、ポアソン比はヤング率に比べ生体中であまり変化しないパラメータであるため、この発明ではポアソン比を０．４９で一定とすることが好ましい。
そして、組織を少なくとも２次元の有限要素モデル化、すなわち、組織を有限個の要素に分割し、各要素内で弾性方程式に歪み分布情報を適用して弾性係数分布情報を演算する。この弾性係数分布演算によれば、精度よく演算可能な軸方向の歪み分布のみから弾性係数分布を再構成することができ、安定した弾性係数分布の演算が行える。
【００２８】
第２の態様において、前記相関演算手段が、前記軸方向には２分の１波長間隔、前記横方向には前記超音波探触子のライン間隔の格子点でのみ相関係数の計算を行うように構成できるのは、第１の態様と同様である。また、この場合において、第１の態様と同様に、前記相関演算手段が、前記直交検波手段によって直交検波された包絡線信号を用いて、前記軸方向に２分の１波長の整数倍だけシフトした直交検波包絡線信号を作成し、前記圧縮前後における前記包絡線信号間の相関係数を各格子点について計算するように構成することができる。
【００２９】
また、第２の態様において、前記弾性係数情報演算手段は、前記被検体組織を等方性弾性体及び近非圧縮性と仮定し、前記被検体組織を有限個の直方体要素に分割して３次元有限要素モデル化し、前記各要素内では、弾性係数、応力、歪みは一様であると仮定し、弾性方程式に前記歪み分布情報を用いて弾性係数分布情報を演算するように構成することができる。
【００３０】
本発明に係る歪み分布表示方法は、第１の態様の超音波診断システムの各構成要素である超音波探触子、直交検波手段、相関演算手段、変位演算手段、歪み演算手段及び表示手段の処理内容を各ステップとすることにより構成できる。
【００３１】
また、本発明に係る弾性係数分布表示方法は、第２の態様の超音波診断システムの各構成要素である超音波探触子、直交検波手段、相関演算手段、変位演算手段、歪み演算手段、弾性係数情報演算手段及び表示手段の処理内容を各ステップとすることにより構成できる。
【００３２】
【発明の実施の形態】
（本発明の原理の説明）
本発明は、生体組織の硬さを定量的に計測するにあたり、２次元変位分布の計算時間を短縮でき、かつ、変位分布の計算精度を向上させるため、生体組織の変位を推定する公知のドプラ法を改良したことを特徴とする。ここで、まず、本発明の原理について説明する。
【００３３】
前述したように、空間相関法とドプラ法はそれぞれに一長一短があり、共に臨床応用に耐えられるものではない。そこで、この２つの手法の長所を組み合わせた「複合自己相関法（ＣＡ法：Ｃｏｍｂｉｎｅｄ Ａｕｔｏｃｏｒｒｅｌａｔｉｏｎ Ｍｅｔｈｏｄ）」を本願の発明者等は提案している。
【００３４】
図８は、本願発明者等が既に提案している複合自己相関法の原理を示す図である。複合自己相関法は、ドプラ法におけるエイリアシングの問題をＲＦ信号の包絡線による相関を用いることによって解決したものである。その具体的な処理は以下のようになる。
【００３５】
まず、組織圧縮前後のＲＦ信号をドプラ法のときと同じように次式のようにモデル化する。
【数０８】

ここで、ｉ ₁ （ｔ）は圧縮前のＲＦ信号、ｉ ₂ （ｔ）は圧縮後のＲＦ信号、Ａ（ｔ）は包絡線、ω ₀ は超音波の中心角周波数、τは時間シフトである。そして、この２つのＲＦ信号をそれぞれ直交検波すると、次式のようなベースバンド信号が得られる。
【数０９】

そして、この２つの信号間の複素相関関数Ｒ ₁₂ （ｔ；ｎ）を次式のように定義する。
【数１０】

ここで、Ｔは超音波の周期、Ｒ _A （ｔ；τ）は包絡線の自己相関関数、ｔ ₀ は相関窓サイズである。また、＊は複素共役を表している。ここで、ｎ＝０の場合は、ドプラ法における自己相関関数の式（数６）に一致する。すなわち、ｎ＝０の場合はドプラ法と同じであり、軸方向変位が超音波の波長の４分の１以上になるとエイリアシングを起こしてしまう。そこで、この問題を克服するために次式で定義される包絡線相関係数Ｃ（ｔ；ｎ）を用いる。
【数１１】

ただし、Ｒ ₁₁ （ｔ；０）は、ｓ ₁ （ｔ）の自己相関関数、Ｒ ₂₂ （ｔ；ｎ）はｓ ₂ （ｔ＋ｎＴ／２）の自己相関関数である。そして、この包絡線相関係数が最大となるｎの値をｋとすると、そのとき（ｎ＝ｋ）のＲ ₁₂ （ｔ；ｋ）の位相φ（ｔ；ｋ）はエイリアシングの起きていない位相となる。これは、包絡線相関を計算する間隔を２分の１波長（周期）に選んだためである。ちなみに、この２分の１波長はエイリアシングを起こさないための最大の間隔である。よって、このφ（ｔ；ｋ）を用いることにより組織圧縮による時間シフトτ及び軸方向変位ｕ _y は次式のように求まる。
【数１２】

ただし、ｃは組織内の音速であり、生体内で一定と仮定する。
【００３６】
上記の処理を各計測点について行い変位分布を推定する手法が複合自己相関法であり、ドプラ法を拡張したような手法となっている。そのため、リアルタイム計測が可能な手法となっている。また、包絡線相関を用いることによってドプラ法では計測不可能であった大変形の場合（超音波の４分の１波長以上の変位が生じる場合）の変位分布推定にも対応している。
【００３７】
図９は、前述の複合自己相関法の基本アルゴリズムを実行する回路構成を示すブロック図である。加圧前直交検波回路（ＱＤ）１３１は、加圧前のエコー信号ｘ（ｔ）を入力し、それぞれ直交検波して、直交検波信号Ｉｘ（ｔ），Ｑｘ（ｔ）信号を、第１相関演算回路１３３及び第１相関係数演算回路１３５０〜１３５Ｎに出力する。第１加圧後直交検波回路（ＱＤ）１３２０は、加圧後のエコー信号ｙ（ｔ）を入力し、それぞれ直交検波して直交検波信号Ｙ（ｔ）＝Ｉｙ＋ｊＱｙ（Ｉｙ（ｔ），Ｑｙ（ｔ））を、第１相関係数演算回路１３４０及び第２相関演算回路１３５０に出力する。第１遅延回路１３４は、エコー信号ｙ（ｔ）を超音波の半期Ｔの１／２だけ遅延させ、遅延したエコー信号ｙ ₁ ＝ｙ（ｔ−Ｔ／２）を第２加圧後直交検波回路（ＱＤ）１３２１に出力する。第２遅延回路１３５は、第１遅延回路１３４によって遅延されたエコー信号ｙ ₁ ＝ｙ（ｔ−Ｔ／２）を同じく超音波の周期Ｔの１／２だけ遅延させ、遅延したエコー信号ｙ ₂ ＝ｙ（ｔ−Ｔ）を次段の第２加圧後直交検波回路（ＱＤ）１３２２（図示せず）に出力する。以後、Ｎ段の遅延回路を用いて順次周期Ｔの整数倍だけ信号を遅延して、遅延した信号を加圧後直交検波回路に供給する。
【００３８】
第１相関演算回路１３３は、信号Ｉｘ，Ｑｘに基づいて相関値Ｒｘｘを演算し、それを各第２相関係数演算回路１３８０〜１３８Ｎに出力する。第２相関演算回路１３４０は、加圧後直交検波回路１３２０からの直交検波信号Ｉｙ（ｔ），Ｑｙ（ｔ）を入力し、信号Ｉｙ，Ｑｙに基づいて相関値Ｒｙｙを演算し、それを第２相関係数演算回路１３８０に出力する。第１相関係数演算回路１３５０は、加圧前直交検波回路１３１からの直交検波信号Ｉｘ（ｔ），Ｑｘ（ｔ）及び第１加圧後直交検波回路１３２０からの直交検波信号Ｉｙ（ｔ），Ｑｙ（ｔ）を入力し、複素ベースバント信号Ｓ _R ，Ｓ _I を求め、それを第３相関演算回路１３６０及び位相差演算回路１３７０に出力する。第３相関演算回路１３６０は、第１相関係数演算回路１３５０からの複素ベースバント信号Ｓ _R ，Ｓ _I を入力し、それに基づいて相関値｜Ｒｘｙ｜を求め、それを第２相関係数演算回路１３８０に出力する。位相差演算回路１３７０は、第１相関係数演算回路１３５０からの複素ベースバント信号Ｓ _R ，Ｓ _I を入力し、それに基づいて位相差φ ₀ （ｔ）を求める。第２相関係数演算回路１３８０は、第１相関演算回路１３３からの相関値Ｒｘｘ、第３相関演算回路１３６０からの相関値｜Ｒｘｙ｜、並びに第２相関演算回路１３４０からの相関値Ｒｙｙを入力し、これらの各相関値に基づいて相関係数Ｃ ₀ （ｔ）を演算し、出力する。
【００３９】
第２加圧後直交検波回路（ＱＤ）１３２１は、第１遅延回路１３４によって遅延されたエコー信号ｙ ₁ ＝ｙ（ｔ−Ｔ／２）を入力し、それぞれ直交検波して直交検波信号Ｙ ₁ （ｔ）＝Ｉｙ ₁ ＋ｊＱｙ ₁ （Ｉｙ ₁ （ｔ），Ｑｙ ₁ （ｔ））を、第１相関係数演算回路１３４１及び第２相関演算回路１３５１に出力する。第２相関演算回路１３４１は、第２加圧後直交検波回路（ＱＤ）１３２１からの直交検波信号Ｉｙ ₁ （ｔ），Ｑｙ ₁ （ｔ）を入力し、その信号Ｉｙ ₁ （ｔ），Ｑｙ ₁ （ｔ）に基づいて相関値Ｒｙ ₁ ｙ ₁ を演算し、それを第２相関係数演算回路１３８１に出力する。第１相関係数演算回路１３５１は、加圧前直交検波回路１３１からの直交検波信号Ｉｘ（ｔ），Ｑｘ（ｔ）、第２加圧後直交検波回路（ＱＤ）１３２１からの直交検波信号Ｉｙ ₁ （ｔ），Ｑｙ ₁ （ｔ）を入力し、複素ベースバント信号Ｓ _R1 ，Ｓ _I1 を求め、それを第３相関演算回路１３６１及び位相差演算回路１３７１に出力する。第３相関演算回路１３６１は、第１相関係数演算回路１３５１からの複素ベースバント信号Ｓ _R1 ，Ｓ _I1 を入力し、それに基づいて相関値｜Ｒｘｙ ₁ ｜を求め、それを第２相関係数演算回路１３８１に出力する。位相差演算回路１３７１は、第１相関係数演算回路１３５１からの複素ベースバント信号Ｓ _R1 ，Ｓ _I1 を入力し、それに基づいて位相差φ ₁ （ｔ）を求める。第２相関係数演算回路１３８１は、第１相関演算回路１３３からの相関値Ｒｘｘ、第３相関演算回路１３６１からの相関値｜Ｒｘｙ ₁ ｜、並びに第２相関演算回路１３４１からの相関値Ｒｙ ₁ ｙ ₁ を入力し、これらの各相関値に基づいて相関係数Ｃ ₁ （ｔ）を演算し、出力する。
【００４０】
以下同様に、第１遅延回路１３５以降の第２加圧後直交検波回路（ＱＤ）１３２２〜１３２Ｎ、第２相関演算回路１３４２〜１３４Ｎ、第１相関係数演算回路１３５２〜１３５Ｎ、第３相関演算回路１３６２〜１３６Ｎ、位相差演算回路１３７２〜１３７Ｎ及び第２相関係数演算回路１３８２〜１３８Ｎは、上述の１段目及び２段目の回路群と同様の処理を実行し、相関係数Ｃ ₂ （ｔ）〜Ｃ _N （ｔ）及び位相φ ₂ （ｔ）〜φ _N （ｔ）を出力する。
【００４１】
上述の複合自己相関法の基本アルゴリズムを実行する回路は、加圧後のエコー信号ｙ（ｔ）を遅延回路１３４〜１３Ｎで超音波の周期Ｔ／２（２分の１波長）だけ遅延し、それを直交検波回路（ＱＤ）１３２０〜１３２Ｎを用いて個別に直交検波している。
【００４２】
前述のように組織圧縮に伴う変位分布が推定されたら、それを空間微分することにより歪み分布が得られる。歪み分布は定性的に組織の弾性特性を表しているものであり、歪み分布からでもかなりの弾性特性に基づいた診断は行える。しかし、肝硬変などの病変部全体が硬くなるような場合には、定量的な弾性係数によって評価しなければ組織診断は難しい。そのため、近年、組織弾性係数分布再構成法についても研究されるようになってきた。しかし、今のところスタンダードな手法はなく、いずれの手法も研究段階であるというのが実状である。
【００４３】
組織弾性係数分布は先にも述べたように組織内部の歪み分布と応力分布から求められる。しかし、応力分布を直接計測することは現状では困難であるため、歪み分布と組織圧縮の際の境界条件から逆問題的に弾性係数分布を再構成することになる。そのため、一般的に逆問題を解くことは難しく、現在提案されている弾性係数再構成法も数少ない。従来から提案されている弾性係数再構成法を以下に説明する。
【００４４】
第１に、１次元を仮定した方法（１次元弾性体を仮定）がある。これは、１次元弾性体を仮定して歪みの逆数を弾性係数とみなす方法である。この方法は弾性係数再構成法ではなく、歪みの逆数を求めるだけであるので、歪みにおける非定量性をそのまま残している。
【００４５】
第２に、弾性方程式から応力項を消去する方法（等方性弾性体、非圧縮性、平面歪み状態を仮定）がある。これは、平面歪み状態を仮定した場合の弾性方程式を変形し、応力項を消去した方程式を用いて組織圧縮の際の境界条件（体表での外部圧力分布、または体表での変位）と歪み分布（せん断歪み成分を含む歪みテンソルの全成分）から組織弾性係数分布を再構成する手法である。ただし、絶対的な弾性係数を推定するには、弾性係数が前もってわかっている領域（参照領域）が必要となる。
【００４６】
第３に、弾性微分方程式を積分する方法（等方性弾性体、非圧縮性、平面応力状態を仮定）がある。これは、平面応力状態を仮定した場合の弾性方程式を変形した応力項を含まない弾性係数に関する微分方程式を体表付近での弾性係数を基準として順次積分していくことにより、歪み分布（せん断歪み成分を含む歪みテンソルの全成分）から組織弾性係数分布を再構成する方法である。そのため、体表付近の弾性係数分布が前もって分かっている領域が必要であり、また体表付近を基準として積分を行っていくので奥に行くほど誤差が積算されるという問題点もある。
【００４７】
第４に、摂動法を用いた手法（等方性弾性体、近非圧縮性、平面歪み状態を仮定）がある。これは、平面歪み状態を仮定した場合の弾性方程式を基にした摂動法により体表での外部圧力分布と超音波ビーム方向（軸方向）の歪み分布とから反復的に組織弾性係数分布を再構成する方法である。
【００４８】
前述の複合自己相関法は、ドプラ法と同じ１次元の処理を基にしているため、超音波プローブが相対的に横方向に移動してしまい、超音波ビーム方向（軸方向）に直交する方向（横方向）の変位が生じてしまった場合（横方向変位が超音波ビーム幅を超えてしまう場合）には、組織圧縮後のＲＦ信号が無相関となってしまい変位推定に失敗してしまうという問題点がある。すなわち、横方向の変位に対応することができずに、超音波ビーム方向（軸方向）の変位成分のみしか推定できないという問題がある。
【００４９】
（実施の形態１）
以下に、生体組織の硬さを定量的に計測するにあたり、２次元変位分布の計算時間を短縮でき、かつ、変位分布の計算精度を向上させることができる本発明の実施形態１について、図面を参照して説明する。
【００５０】
図１０は、本発明に係る超音波診断システムの好ましい一例の実施形態１の超音波診断システムの概略構成を示すブロック図である。この超音波診断システムでは、包絡線相関計算の際、複合自己相関法で１次元の相関窓で１次元探索していた処理を２次元の相関窓を用いて２次元探索することにより横方向の移動にも対応した拡張複合自己相関法と呼ばれる方法を採用している。この拡張複合自己相関法は、軸方向には２分の１波長間隔、横方向にはライン間隔の格子点でのみ包絡線相関計算を行うことにより計算量を減少させて高速化を図っている。ただし、複合自己相関と同様に拡張複合自己相関法でも位相情報を利用して軸方向の変位推定精度を向上させている。しかし、横方向変位の推定はキャリアとなる信号がないため位相情報は利用できない。そのため、横方向変位推定精度は空間相関法と同様に横方向サンプリング間隔（ライン間隔）により制限されてしまう。しかし、後で提案する弾性係数分布再構成法では軸方向の歪み（変位）分布のみから弾性係数分布を推定できるため、ここでは横方向変位推定精度の向上は特に行わない。この拡張複合自己相関法の具体的な構成について図１０を用いて説明する。
【００５１】
図１０において、超音波プローブ９１は、被検体内へ超音波を送波すると共にその反射波を受波するものであり、従来のセクタスキャンプローブ（セクタフェイズドアレイプローブ）、リニアスキャンプローブ（リニアアレイプローブ）又はコンベックススキャンプローブ（コンベックスアレイプローブ）などである。超音波プローブ９１からは、組織圧縮前後のＲＦ信号が直交検波器９２に出力される。直交検波器９２は、組織圧縮前後のＲＦ信号をそれぞれ組織圧縮前後の複素包絡線信号（ＩＱ信号）に変換し、複素２次元相関計算部９３に出力する。複素２次元相関計算部９３は、組織圧縮前後のＲＦ信号間における２次元相関を計算し、その相関が最大となる位置を横方向変位計算部９４及び軸方向変位計算部９５に出力し、そのときの相関関数の位相を軸方向変位計算部９５に出力する。ただし、軸方向にはエイリアシングを起こさずに位相を検出できる最大の間隔である超音波中心周波数の２分の１波長間隔でのみ相関を計算するものとする。これは、超音波診断システムのリアルタイム表示を優先させるためである。従って、高精度な相関を計算するためには、この２分の１波長間隔に限定する必要はない。
【００５２】
横方向変位計算部９４は、複素２次元相関計算部９３からの横方向の相関最大位置に基づいて横方向の変位ｕ_x を計算し、それを横方向歪み計算部９６に出力する。一方、軸方向変位計算部９５は、複素２次元相関計算部９３からの軸方向の相関最大位置及びそのときの位相に基づいて軸方向の変位ｕ_y を計算し、それを軸方向歪み計算部９７に出力する。横方向歪み計算部９６は、横方向変位計算部９４からの横方向変位ｕ_x の分布を空間的に微分することにより横方向歪み分布ε_x を計算し、それを量子化部９８に出力する。一方、軸方向歪み計算部９７は、軸方向変位計算部９５からの横方向変位ｕ_y の分布を空間的に微分することにより軸方向歪み分布ε_y を計算し、それを量子化部９８に出力する。量子化部９８は、横方向歪み分布ε_x 及び軸方向歪み分布ε_y をグレースケール表示（又はカラー表示）するために各歪み分布を量子化し、表示部９９に出力する。表示部９９は、量子化された各歪み分布を表示する。
【００５３】
次に、図１０の超音波診断システムで採用した拡張複合自己相関法の動作について説明する。まず、組織圧縮が極僅か（数パーセント以下）である場合、組織内部を局所的に見れば平行移動したと見なすことができ、組織圧縮前後のＲＦ信号を次式のようにモデル化できる。
【数１３】

ここで、ｉ₁ （ｔ，ｘ）は圧縮前のＲＦ信号、ｉ₂ （ｔ，ｘ）は圧縮後のＲＦ信号、Ａ（ｔ，ｘ）は包絡線、ω₀ は超音波の中心角周波数、τは時間シフト、ｕ_x は横方向変位である。また、ここではドプラ法や複合自己相関法のときと違い横方向の変位も考慮して圧縮前後のＲＦ信号をモデル化している。そして、この式の中で最終的に知りたいパラメータは、軸方向の変位ｕ_y ＝ｃτ／２（すなわち、時間シフトτ）と横方向変位ｕ_x である。ただし、ｃは組織内の音速であり、生体内で一定と仮定する。
【００５４】
そこで、まずこれらの組織圧縮前後のＲＦ信号を直交検波器９２でそれぞれ直交検波する。すなわち、各ＲＦ信号に超音波の中心周波数と同じ周波数のｓｉｎ波，ｃｏｓ波をかけ、それぞれ低域通過フィルタをかける。すると、以下のような複素ベースバンド信号ｓ₁ ，ｓ₂ が得られる。
【数１４】

そして、このｓ₁ （ｔ，ｘ）とｓ₂ （ｔ＋ｎＴ／２，ｘ＋ｍＬ）との間の２次元複素相関関数Ｒ₁₂（ｔ，ｘ；ｎ，ｍ）を次式のように定義する。
【数１５】

ここで、Ｔは超音波の周期、Ｌは横方向サンプリング間隔（ライン間隔）、Ｒ_A （ｔ，ｘ；τ，ｕ_x ）は包絡線の自己相関関数、ｔ₀ は軸方向相関窓サイズ、ｘ₀ は横方向相関窓サイズである。また、＊は複素共役を表している。そして、この２次元複素相関関数を用いて２次元包絡線相関係数Ｃ（ｔ，ｘ；ｎ，ｍ）を以下のように定義する。
【数１６】

ただし、Ｒ₁₁（ｔ，ｘ；０，０）はｓ₁ （ｔ，ｘ）の自己相関関数、Ｒ₂₂（ｔ，ｘ；ｎ，ｍ）はｓ₂ （ｔ＋ｎＴ／２，ｘ＋ｍＬ）の自己相関関数である。そして、この包絡線相関係数を用いて複合自己相関法の場合と同様にエイリアシングの問題を克服する。すなわち、各計測点（ｔ，ｘ）におけるＣ（ｔ，ｘ；ｎ，ｍ）とＲ₁₂（ｔ，ｘ；ｎ，ｍ）の位相φ（ｔ，ｘ；ｎ，ｍ）との組｛Ｃ（ｔ，ｘ；ｎ，ｍ），φ（ｔ，ｘ；ｎ，ｍ）｝をすべてのｎとｍについて求める。ここで、ｎとｍの範囲が十分広ければ、すなわち、包絡線相関を行う探索範囲が十分に大きければ、包絡線相関係数が最大となる（ｎ，ｍ）＝（ｋ，ｌ）に対する位相φ（ｔ，ｘ；ｋ，ｌ）はエイリアシングの起きていない位相となる。これは、包絡線相関Ｃ（ｔ，ｘ；ｎ，ｍ）が最大となる（ｎ，ｍ）＝（ｋ，ｌ）のとき、ｓ₁ （ｔ，ｘ）とｓ₂ （ｔ＋ｋＴ／２，ｘ＋ｌＬ）との時間シフトの大きさ｜τ−ｋＴ／２｜がＴ／２よりも小さくなる、すなわち、｜φ（ｔ，ｘ；ｋ，ｌ）｜＝ω₀ ｜τ−ｋＴ／２｜がπよりも小さくなるためである。よって、このエイリアシングの起きていないφ（ｔ，ｘ；ｋ，ｌ）を用いれば、計測点（ｔ，ｘ）における正確な時間シフトτ、軸方向変位ｕ_y 、そして横方向変位ｕ_x が次式のように求まる。
【数１７】

ただし、ｃは組織内での音速（ここでは軟組織における一般的な音速１５００ｍ／ｓで一定とする）である。したがって、組織内のすべての点で上記のように軸方向変位と横方向変位を計算すれば、軸方向変位分布ｕ_y （ｘ，ｙ）と横方向変位分布ｕ_x （ｘ，ｙ）が得られる。すなわち、言い換えれば、実施の形態１は、被検体との間で超音波信号を送受信する超音波探触子により受信された受信信号を記憶しておき、記憶された被検体の圧縮前後のフレームデータに計測点を設定し、その計測点を囲む２次元相関窓を用いて予め設定された２次元探索範囲内で移動して、計測点における前記圧縮前後の受信信号の包絡線信号の相関係数が最大となる相関最大位置及び該相関最大位置における前記圧縮前後の前記受信信号間の位相差を求め、これによって求められた相関係数が最大となる位置及び位相差に基づいて圧縮に伴う計測点の変位を求めることにより、２次元の変位を求めることができる。
【００５５】
また、各変位分布を次式のように空間微分することにより、軸方向歪み分布ε_y （ｘ，ｙ）と横方向歪み分布ε_x （ｘ，ｙ）が次式のように求められる。
【数１８】

以上のようにして、実施形態１によれば、組織圧縮前後のＲＦ信号から軸方向と横方向の変位（歪み）分布を推定することができる。ただし、上式のｕ_x ＝ｌＬからもわかるように横方向変位の精度は横方向サンプリング間隔（ライン間隔）によって制限されてしまうため、精度は若干劣るということはあるが、リアルタイムに観察できるので実用性の高いものである。
【００５６】
（実施の形態２）
実施の形態１の拡張複合自己相関法は、組織の横方向移動に対応するように２次元の相関窓と探索範囲を用いて、組織圧縮の際の組織と超音波プローブの相対的な横方向移動には対応している。しかし、組織圧縮の際に軸方向と横方向にそれぞれ垂直な方向（２次元超音波走査面に垂直な方向（スライス方向））の変位が生じてしまい組織移動が起こった場合には、２次元の拡張複合自己相関法では歪みの推定を行うことができない。つまり、実施の形態１では、問題を簡単化するために２次元状態（平面歪み状態や平面応力状態）を仮定しているが、本来３次元構造をなす生体組織を２次元で近似すると弾性係数を過小評価する恐れがあり好ましくない。これは、２次元状態では考察面に垂直な方向の歪みや応力を考慮に入れていないことが原因である。
【００５７】
そのため、より安定したシステムにするには、上述の拡張複合自己相関法を、以下に述べる実施の形態２に示すように、３次元の相関窓と３次元の探索範囲を用いることにより簡単に拡張することが可能である。
【００５８】
図１１及び図１２は、本発明の実施の形態２の３次元複合自己相関法の基本アルゴリズムを示すフローチャート図である。なお、図１２は、図１１の処理の一部の詳細を示すフローチャート図である。
【００５９】
ステップＳ１１では、ステップＳ２３の判定処理と組み合わせて、第１番目の走査線から第Ｌ番目の走査線についてそれぞれ同様の処理を行うために、走査線番号レジスタｌに「１」を格納する。
【００６０】
ステップＳ１２では、ステップＳ１８の判定処理と組み合わせて、厚み方向（フレーム）を−ＵからＵまで順次シフトする処理を実行する。
【００６１】
ステップＳ１３では、ステップＳ１７の判定処理と組み合わせて、方位方向（走査線）を−ＶからＶまで順次シフトする処理を実行する。
【００６２】
ステップＳ１４では、ステップＳ１６の判定処理と組み合わせて、距離方向（軸方向）を０からＭまで順次シフトする処理を実行する。
【００６３】
ステップＳ１５では、複合自己相関法により、距離方向（軸方向）における包絡線の相関係数Ｃ（ｌ，ｔ；ｕ，ｖ，ｎ）を計算する。この複合自己相関法は、従来の方法でやってもいいが、それだと計算に時間を要するので、ここでは、高速化された複合自己相関法を用いて相関係数Ｃ（ｌ，ｔ；ｕ，ｖ，ｎ）の計算を行う。この高速化複合自己相関法については後述する。
【００６４】
ステップＳ１６では、前のステップＳ１４と組み合わせられた処理であり、距
離方向レジスタｎがその最大値Ｍに達したか否かの判定を行い、達した場合には
ステップＳ１７に進み、そうでない場合はステップＳ１４にリターンし、距離方向レジスタｎをインクリメント処理する。
【００６５】
ステップＳ１７では、前のステップＳ１３と組み合わせられた処理であり、方位方向レジスタｖがその最大値Ｖに達したか否かの判定を行い、達した場合にはステップＳ１８に進み、そうでない場合はステップＳ１３にリターンし、方位方向レジスタｖをインクリメント処理する。
【００６６】
ステップＳ１８では、前のステップＳ１２と組み合わせられた処理であり、厚み方向レジスタｕがその最大値Ｕに達したか否かの判定を行い、達した場合にはステップＳ１９に進み、そうでない場合はステップＳ１２にリターンし、厚み方向レジスタｕをインクリメント処理する。
【００６７】
ステップＳ１９では、ステップＳ１２〜ステップＳ１８の処理によって求めれた相関係数Ｃ（ｌ，ｔ；ｕ，ｖ，ｎ）（ｕ＝−Ｕ，…０，…，Ｕ）（ｖ＝−Ｖ，…０，…，Ｖ）（ｎ＝０，１，…，Ｎ）の中から最大となる（ｕ，ｖ，ｎ）を求め、それを（ｕ₀ ，ｖ₀ ，ｎ₀ ）とする。
【００６８】
ステップＳ２０では、ステップＳ１９で求められた相関係数Ｃ（ｌ，ｔ；ｕ₀ ，ｖ₀ ，ｎ₀ ）について、その位相差φ（ｌ，ｔ；ｕ₀ ，ｖ₀ ，ｎ₀ ）を計算する。
【００６９】
ステップＳ２１では、最終的を位相差として、φ＝ｎ₀ π＋φ（ｌ，ｔ；ｕ₀ ，ｖ₀ ，ｎ₀ ）を計算する。
【００７０】
ステップＳ２２では、（ｕ₀ ，ｖ₀ ，ｎ₀ ）の近傍の相関係数Ｃ（ｌ，ｔ；ｕ，ｖ，ｎ）を用いて、方位方向の変位：ｖ＝ｖ₀ ＋Δｖ及び厚み方向の変位：ｕ＝ｕ₀ ＋Δｕを計算する。
【００７１】
ステップＳ２３では、前のステップＳ１１と組み合わせられた処理であり、走査線番号レジスタｌがＬに達したか否かの判定を行い、達した場合にはステップＳ２４に進み、そうでない場合はステップＳ１１にリターンし、走査線番号レジスタｌをインクリメント処理する。
【００７２】
ステップＳ２４では、組織圧縮に伴う変位分布が推定されたら、それを空間微分することにより歪み分布を計算する。
【００７３】
図１３は、図１２のステップＳ１５の高速化された複合自己相関法の詳細を示すフローチャート図である。
【００７４】
ステップＳ３１では、圧縮前のＲＦ信号の包絡線ｘと、圧縮後のＲＦ信号の包絡線をそれぞれ直交検波して、以下のようにＩ，Ｑ信号を求める。
ｘ（ｔ）→Ｉｘ、Ｑｘ（Ｘ（ｔ）＝Ｉｘ＋ｊＱｘとする）
ｙ（ｔ）→Ｉｙ、Ｑｙ（Ｙ（ｔ）＝Ｉｙ＋ｊＱｙとする）
ステップＳ３２では、相関：Ｒｘｙ、Ｒｘｘ、Ｒｙｙを次式に基づいて計算する。
Ｒｘｙ＝∫Ｘ（ｔ＋ν）・Ｙ^* （ｔ＋ν）ｄν
Ｒｘｘ＝∫Ｘ（ｔ＋ν）・Ｘ^* （ｔ＋ν）ｄν
Ｒｙｙ＝∫Ｙ（ｔ＋ν）・Ｙ^* （ｔ＋ν）ｄν
ステップＳ３３では、求められた相関Ｒｘｙ、Ｒｘｘ、Ｒｙｙを用いて相関係数Ｃ₀ を次式に基づいて計算する。
Ｃ₀ ＝｜Ｒｘｙ｜／√Ｒｘｘ√Ｒｙｙ
ステップＳ３４では、変数ｎに１をセットする。
【００７５】
ステップＳ３５では、Ｙ_n （ｔ）＝Ｙ（ｔ−ｎＴ）ｅ^{j ω nT}を計算する。
【００７６】
ステップＳ３６では、次式に基づいてＲｘｙ_n ，Ｒｙ_nｙ_nを計算する。
Ｒｘｙ_n ＝∫Ｘ（ｔ＋ν）・Ｙ_n ^* （ｔ＋ν）ｄν
＝∫Ｘ（ｔ＋ν）・Ｙ^* （ｔ−ｎＴ＋ν）ｅ^{j ω nT}ｄν
Ｒｙ_nｙ_n＝∫Ｙ_n（ｔ＋ν）Ｙ_n ^* （ｔ＋ν）ｄν
＝∫Ｙ（ｔ−ｎＴ＋ν）・Ｙ^* （ｔ−ｎＴ＋ν）ｄν
ステップＳ３７では、求められたＲｘｙ_n ，Ｒｙ_nｙ_nを用いて相関係数Ｃ_n を次式に基づいて計算する。
Ｃ_n ＝｜Ｒｘｙ_n ｜／√Ｒｘｘ√Ｒｙ_nｙ_n
ステップＳ３８では、変数ｎをインクリメント処理する。
ステップＳ３９では、変数ｎが最大値Ｍに達したか否かを判定し、達した場合はこの処理を終了し、達していない場合は、ステップＳ３５にリターンし、同様の処理を繰り返す。
【００７７】
図１３のフローチャートでは、Ｒｘｙ_n ，Ｒｙ_nｙ_nを求めるのに、ステップＳ３５でＹ_nをＹから導いている。このために、回路構成の簡略化を図ることができる。以下、どのようにしてＹ_n をＹから導くかについて説明する。
まず、加圧前のエコー信号ｘ（ｔ）を
ｘ（ｔ）＝ｕ（ｔ）ｃｏｓ（ωｔ＋θ）
軸方向に加圧後のエコー信号ｙ（ｔ）を
ｙ（ｔ）＝ｘ（ｔ＋τ）＝ｕ（ｔ＋τ）ｃｏｓ（ω（ｔ＋τ）＋θ）
とする。
各信号ｘ，ｙ，ｙ_n の直交検波値は、
ｘ（ｔ）→Ｉｘ＝０．５ｕ（ｔ）ｃｏｓθ
Ｑｘ＝−０．５ ｕ（ｔ）ｓｉｎθ
（Ｘ（ｔ）＝Ｉｘ＋ｊＱｘ＝０．５ｕ（ｔ）ｅ^{-j θ}）
ｙ（ｔ）→Ｉｙ＝０．５ｕ（ｔ＋τ）ｃｏｓ（ωτ＋θ）
Ｑｙ＝−０．５ｕ（ｔ＋τ）ｓｉｎ（ωτ＋θ）
（Ｙ（ｔ）＝Ｉｙ＋ｊＱｙ＝０．５ｕ（ｔ＋τ）ｅ^{-j( ωτ + θ )} ）
ｙ_n （ｔ）＝ｙ（ｔ−ｎＴ）
＝ｕ（ｔ＋τ−ｎＴ）ｃｏｓ（ω（ｔ＋τ−ｎＴ）＋θ）
＝ｕ（ｔ＋τ−ｎＴ）ｃｏｓ（ωｔ＋ω（τ−ｎＴ）＋θ）
となる。ここで、Ｔは半周期なので、
Ｉｙ_n＝０．５ｕ（ｔ＋τ−ｎＴ）ｃｏｓ（ω（τ−ｎＴ）＋θ）
Ｑｙ_n＝−０．５ｕ（ｔ＋τ−ｎＴ）ｓｉｎ（ω（τ―ｎＴ）＋θ）
（Ｙ_n＝Ｉｙ_n＋ｊＱｙ_n＝０．５ｕ（ｔ＋τ−ｎＴ）ｅ^{-j( ω ( τ -nT)+ θ )}）
となる。以上の式から以下のような関係が成り立つ。
Ｙ_n （ｔ）＝Ｉｙ_n＋ｊＱｙ_n
＝０．５ｕ（ｔ＋τ−ｎＴ）ｅ^{-j( ω ( τ -nT)+ θ )}
＝Ｙ（ｔ−ｎＴ）ｅ^{j ω nT}
これからＹ_n （ｔ）はＹ（ｔ）＝Ｉｙ＋ｊＱｙから求まることになる。
従って、Ｒｘｙ_n，Ｒｙ_nｙ_nも、次式のようにＸ、Ｙから求めることができる。
Ｒｘｙ_n ＝４∫Ｘ（ｔ＋ν）・Ｙ_n ^* （ｔ＋ν）ｄν
＝４∫Ｘ（ｔ＋ν）・Ｙ^*（ｔ−ｎＴ＋ν）ｅ^{j ω nT}ｄν
｜Ｒｘｙ_n ｜＝Ｒｕ_n
＝∫ｕ（ｔ＋ν）ｕ（ｔ＋τ−ｎＴ＋ν）ｄν
＝４∫｜Ｘ（ｔ＋ν）・Ｙ_n ^* （ｔ＋ν）｜ｄν
＝４∫｜Ｘ（ｔ＋ν）・Ｙ^*（ｔ−ｎＴ＋ν）ｅ^{j ω nT}｜ｄν
＝４∫｜Ｘ（ｔ＋ν）・Ｙ^*（ｔ−ｎＴ＋ν）｜ｄν
Ｒｙ_nｙ_n＝∫ｕ（ｔ＋τ−ｎＴ＋ν）ｕ（ｔ＋τ−ｎＴ＋ν）ｄν
＝４∫｜Ｙ_n（ｔ＋ν）・Ｙ_n ^*
（ｔ＋ν）｜ｄν
＝４∫Ｙ（ｔ−ｎＴ＋ν）・Ｙ^*（ｔ−ｎＴ＋ν）ｄν
ここで、＊は複素共役を表している。
【００７８】
図１４は、この３次元複合自己相関法の基本アルゴリズムを実行する回路構成を示すブロック図である。複合自己相関法を実行する回路構成を図９に示すようなものにすると、直交検波回路１３２０〜１３２Ｎが多段接続されることによって、直交検波回路１３２０〜１３２Ｎの処理に時間を要する。そのため、その処理時間が膨大なものとなってしまい、高速な演算処理の妨げとなり、リアルタイムな画像表示の妨げとなる。そこで、前述のような基本アルゴリズムに応じた図１４のような回路構成を採用することによって、複合自己相関法を実行する回路の処理速度を高速化している。
【００７９】
加圧前直交検波回路（ＱＤ）１３１は、加圧前のエコー信号ｘ（ｔ）を入力し、それぞれ直交検波して、直交検波信号Ｉｘ（ｔ），Ｑｘ（ｔ）信号を、第１相関演算回路１３３及び第１相関係数演算回路１３５０〜１３５Ｎに出力する。加圧後直交検波回路（ＱＤ）１３２は、加圧後のエコー信号ｙ（ｔ）を入力し、それぞれ直交検波して直交検波信号Ｙ（ｔ）＝Ｉｙ＋ｊＱｙ（Ｉｙ（ｔ），Ｑｙ（ｔ））を、第１相関係数演算回路１３５０、第２相関演算回路１３４０及び第１遅延回路１３４及び第２遅延回路１３５に出力する。第１遅延回路１３４及び第２遅延回路１３５は、直交検波信号Ｙ（ｔ）をそれぞれ超音波の周期Ｔの１／２だけ遅延させ、遅延した直交検波信号Ｙ（ｔ−Ｔ／２）を第１相関係数演算回路１３５１、第３遅延回路１３６及び第４遅延回路１３７に出力する。第３遅延回路１３６及び第４遅延回路１３７は、直交検波信号Ｙ（ｔ−Ｔ／２）をそれぞれ超音波の周期Ｔの１／２だけ遅延させ、遅延した直交検波信号Ｙ（ｔ−Ｔ）を次段の第１相関係数演算回路及び遅延回路（図示せず）に出力する。以後、Ｎ段の遅延回路を用いて順次周期Ｔの１／２の整数倍だけ信号を遅延して、遅延した信号を第１相関係数演算回路に供給する。
【００８０】
第１相関演算回路１３３は、信号Ｉｘ，Ｑｘに基づいて相関値Ｒｘｘを演算し、それを各第２相関係数演算回路１３８０〜１３８Ｎに出力する。第２相関演算回路１３４０は、加圧後直交検波回路１３２からの直交検波信号Ｉｙ（ｔ），Ｑｙ（ｔ）を入力し、信号Ｉｙ，Ｑｙに基づいて相関値Ｒｙｙを演算し、それを第２相関係数演算回路１３８０に出力する。第１相関係数演算回路１３５０は、加圧前直交検波回路１３１からの直交検波信号Ｉｘ（ｔ），Ｑｘ（ｔ）及び加圧後直交検波回路１３２からの直交検波信号Ｉｙ（ｔ），Ｑｙ（ｔ）を入力し、複素ベースバント信号Ｓ_R ，Ｓ_I を求め、それを第３相関演算回路１３６０及び位相差演算回路１３７０に出力する。第３相関演算回路１３６０は、第１相関係数演算回路１３５０からの複素ベースバント信号Ｓ_R ，Ｓ_I を入力し、それに基づいて相関値｜Ｒｘｙ｜を求め、それを第２相関係数演算回路１３８０に出力する。
位相差演算回路１３７０は、第１相関係数演算回路１３５０からの複素ベースバント信号Ｓ_R ，Ｓ_I を入力し、それに基づいて位相差φ₀ （ｔ）を求める。第２相関係数演算回路１３８０は、第１相関演算回路１３３からの相関値Ｒｘｘ、第３相関演算回路１３６０からの相関値｜Ｒｘｙ｜、並びに第２相関演算回路１３４０からの相関値Ｒｙｙを入力し、これらの各相関値に基づいて相関係数Ｃ₀ （ｔ）を演算し、出力する。
【００８１】
第２相関演算回路１３４１は、第１遅延回路１３４及び第２遅延回路１３５からの遅延後の直交検波信号Ｉｙ（ｔ−Ｔ／２），Ｑｙ（ｔ−Ｔ／２）を入力し、信号Ｉｙ（ｔ−Ｔ／２），Ｑｙ（ｔ−Ｔ／２）に基づいて相関値Ｒｙ₁ｙ₁を演算し、それを第２相関係数演算回路１３８１に出力する。第１相関係数演算回路１３５１は、加圧前直交検波回路１３１からの直交検波信号Ｉｘ（ｔ），Ｑｘ（ｔ）、第１遅延回路１３４及び第２遅延回路１３５からの遅延後の直交検波信号Ｉｙ（ｔ−Ｔ／２），Ｑｙ（ｔ−Ｔ／２）を入力し、複素ベースバント信号Ｓ_R1，Ｓ_I1を求め、それを第３相関演算回路１３６１及び位相差演算回路１３７１に出力する。第３相関演算回路１３６１は、第１相関係数演算回路１３５１からの複素ベースバント信号Ｓ_R1，Ｓ_I1を入力し、それに基づいて相関値｜Ｒｘｙ₁ ｜を求め、それを第２相関係数演算回路１３８１に出力する。位相差演算回路１３７１は、第１相関係数演算回路１３５１からの複素ベースバント信号Ｓ_R1，Ｓ_I1を入力し、それに基づいて位相差φ₁ （ｔ）を求める。第２相関係数演算回路１３８１は、第１相関演算回路１３３からの相関値Ｒｘｘ、第３相関演算回路１３６１からの相関値｜Ｒｘｙ₁ ｜、並びに第２相関演算回路１３４１からの相関値Ｒｙ₁ｙ₁を入力し、これらの各相関値に基づいて相関係数Ｃ₁ （ｔ）を演算し、出力する。
【００８２】
第３遅延回路１３５及び第４遅延回路１３６から次段の第２相関演算回路１３４２〜１３４Ｎ、第１相関係数演算回路１３５２〜１３５Ｎ、第３相関演算回路１３６２〜１３６Ｎ、位相差演算回路１３７２〜１３７Ｎ及び第２相関係数演算回路１３８２〜１３８Ｎは、上述と同様の処理を順次遅延された遅延後の直交検波信号Ｉｙ（ｔ−２Ｔ／２）・・・Ｉｙ（ｔ−ＮＴ／２），Ｑｙ（ｔ−２Ｔ／２）・・・Ｑｙ（ｔ−ＮＴ／２）に対して実行し、相関係数Ｃ₂ （ｔ）〜Ｃ_N （ｔ）及び位相φ₂ （ｔ）〜φ_N （ｔ）を出力する。
【００８３】
次に、３次元有限要素モデルを用いた弾性係数分布再構成法について説明する。弾性係数分布再構成逆問題を簡単化するため、この実施の形態では組織をモデル化する。これはまた、本実施の形態の弾性係数分布再構成法において有限要素法を用いるためでもある。この実施の形態では、組織を以下のように仮定及びモデル化する。
【００８４】
まず、組織を等方性弾性体と仮定する。組織歪み分布を推定する際、外部から圧力を加えて組織を静的に圧縮する。しかし、組織圧縮前後のＲＦ信号間の相関を保つために、微小圧縮しか行わない。そのため、この場合、応力と歪みの間の関係はほぼ線形である。すなわち、組織を弾性体として近似できる。また、今回対象としている軟組織はほぼ等方性が成り立つため、この実施の形態では組織を等方性弾性体と仮定する。
【００８５】
さらに、組織を近非圧縮性と仮定する。生体組織は、非圧縮性（ポアソン比ν＝０．５）に近いことが知られている。そこで、ポアソン比を０．４９とし、生体内で一定とする。ここで、完全な非圧縮性を仮定しないのは、ポアソン比ν＝０．５とすると特殊な弾性方程式となり、本実施の形態で用いている有限要素法が適用できなくなるためである。そして、ポアソン比を生体内で一定とすることで弾性係数分布推定の推定パラメータをヤング率のみとすることができ、逆問題を簡単化できる。また、ポアソン比はヤング率に比べ生体中であまり変化しないパラメータであるため、この実施の形態ではポアソン比を０．４９で一定とする。
【００８６】
組織を３次元有限要素モデル化する。この弾性係数分布再構成法では有限要素法を用いるため、組織を有限個の直方体要素に分割する。そして、各要素内では、弾性係数、応力、歪みは一様であると仮定する。一般的に逆問題を解くには、それに対応する順問題を理解することが重要である。今回の歪み分布と境界条件から弾性係数分布を推定する逆問題の場合、それに対応する順問題とは、弾性係数分布と境界条件から歪み分布を求めることである。そして、この順問題の数値解法の１つが有限要素法（ＦＥＭ：Ｆｉｎｉｔｅ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ）である。
【００８７】
ここで、有限要素法とは対象となる連続体を有限個の要素の集合で近似し、この集合体に対して成り立つ連立１次方程式を数値的に解く手法のことである。なお、有限要素法の定式化については後述する。ここでは有限要素法とは「入力として物体の弾性係数分布と境界条件を与えれば、出力としてそのときの歪み（変位）分布と応力分布が得られるもの」として捉えておけば十分である。
【００８８】
この実施の形態では、組織を等方性弾性体で近似するため、組織内では以下のような弾性方程式（つりあい方程式・歪み−変位関係式・応力−歪み関係式）が成り立つ。
【００８９】
つりあい方程式は次式のように表される。
【数１９】

歪み−変位関係式は次式のように表される。
【数２０】

応力−歪み関係式（一般化したフックの法則）は次式のように表される。
【数２１】

上記の式ではテンソル表現を用いており、実際にはつりあい方程式として３式、歪み−変位方程式として６式、応力−歪み関係式として６式が存在する。また、座標系（ｘ₁ ，ｘ₂ ，ｘ₃ ）は（ｘ，ｙ，ｚ）を表しており、その他の記号は以下のことを表している。
Ｅ：ヤング率（弾性係数とはヤング率のことを表している）
ν：ポアソン比
ε_ij：歪みテンソル
（ε_nn＝ε₁₁＋ε₂₂＋ε₃₃：体積歪み）
σ_ij：応力テンソル
δ_ij：クロネッカーのデルタ
ｕ_i ：変位ベクトル
ｆ_i ：体積力ベクトル（重力の影響は無視できるため、ここではｆ_i ＝０とする）
ここで、応力−歪み関係式をε_ijについて整理すると、次のような歪み−応力関係式が得られる。
【数２２】

ただし、σ_nn＝σ₁₁＋σ₂₂＋σ₃₃である。よって、この式の中でｉ＝ｊ＝２とし、ヤング率Ｅについて整理すると次式が得られる。
【数２３】

従って、軸方向（この実施の形態では、ｙ方向を超音波ビーム方向、すなわち軸方向とする）の歪み成分と全方向の応力成分がわかれば、ヤング率すなわち弾性係数を求めることができる。なお、上述の計算式からは、応力分布を直接計測することは現状では困難である。そこで、この実施の形態では応力分布と弾性係数分布を交互に推定・更新しながら、推定弾性係数分布を実際の分布に近づけていく。その弾性係数分布再構成の手順は、以下のようになる。
【００９０】
第１に、未知弾性係数分布の初期値分布｛Ｅ^⁰｝として一様分布を考える。第２に、初期弾性係数分布｛Ｅ^⁰｝のときに生じる応力分布｛σ^⁰｝を３次元有限要素法により求める。具体的には、まず組織モデル内の各要素に対して歪み−変位関係式及び応力−歪み関係式をつりあい方程式に代入して得られる次式のようなつりあい方程式を作る。
【数２４】

ただし、
【数２５】

である。そして、この連立方程式を以下のような境界条件のもとガウスの消去法を用いて変位について解き、弾性係数分布｛Ｅ^⁰｝のときの変位分布｛ｕ^⁰｝を求める。
【数２６】

上式において、ｐｉは体表における外部圧力ベクトル、σｎは側面に垂直な方向の応力成分である。また、上段の式は底面が固定されていることを示し、中段の式は体表での応力分布は外部圧力分布に等しいことを示し、下段の式は側面が拘束されていないことをそれぞれ示している。次に、この変位分布｛ｕ^⁰｝を歪み−変位関係式に代入して、弾性係数分布｛Ｅ^⁰｝のときの歪み分布｛ε^⁰｝を求める。そして、この歪み分布｛ε^⁰｝を応力−歪み関係式に代入することにより、弾性係数分布｛Ｅ^⁰｝のときの応力分布｛σ^⁰｝を得る。
【００９１】
第３に、３次元有限要素法により得られた応力分布と拡張複合自己相関法により推定した軸方向（ｙ方向）歪み分布｛ε_y ｝を用いて、弾性係数分布｛Ｅ^^k｝を次式により更新する。
【数２７】

ただし、この式は、上述の応力−歪み関係式をε_ijについて整理し、式中のｉ＝ｊ＝２とし、ヤング率Ｅについて整理した式を書き改めたものであり、式中のｋは繰り返し回数を表している。
【００９２】
第４に、上述のように更新された弾性係数分布と上述の境界条件を用いて再び３次元有限要素解析を行い、応力分布を更新する。
【００９３】
そして、第３及び第４の処理を繰り返すことにより弾性係数分布を実際の分布に近づけていく。ただし、次式の条件が満たされた時点で弾性係数分布推定は収束したとみなし、推定を終了する。
【数２８】

ここで、ｌは要素番号、Ｎは要素数、Γはしきい値である。
【００９４】
以上が、実施の形態２の３次元有限要素モデルによる弾性係数分布再構成法であり、この方法は３次元のつりあい方程式を基に弾性係数分布を推定している。そのため、本手法は従来の手法よりもより実際の生体組織に近い仮定に基づいているので、より正確な弾性係数推定が可能になる。また、本手法は精度良く推定可能な軸方向の歪み分布のみから弾性係数分布を再構成するため、安定した弾性係数分布再構成が行える。ただし、本手法は組織弾性係数の３次元分布を推定する手法であるため、２次元アレイ超音波プローブを用いるか、１次元アレイ超音波プローブをスライス方向に機械的に走査することにより、対象を３次元的に走査する必要がある。
【００９５】
（実施形態１，２のシミュレーション）
本発明の実施の形態１，２の拡張複合自己相関法と３次元有限要素モデルによる弾性係数分布再構成法の有効性をシミュレーションによって実証する。図１５は、このシミュレーションの手順の概略を示す図である。
【００９６】
第１に、推定したい弾性係数分布を持つ組織モデルを作成する。このとき、組織モデル内には超音波エコー信号を発生させるための散乱体を分布させておく。第２に、この組織モデルに対して外部圧力を加え、計算機上で組織圧縮を行う。そして、この圧縮による各散乱点の移動先を有限要素法などにより求める。第３に、組織モデル圧縮前後の散乱体分布を基に圧縮前後のＲＦ信号を作成する。第４に、この圧縮前後のシミュレーションＲＦ信号に対して拡張複合自己相関法を適用し、組織歪み分布を推定する。第５に、拡張複合自己相関法により推定された歪み分布と組織モデル圧縮の際に設定した境界条件（外部圧力分布など）とから３次元有限要素モデルによる弾性係数分布再構成法により組織弾性係数分布を推定する。
【００９７】
今回のシミュレーションで用いた組織モデルの弾性係数分布は、各シミュレーションにおいて異なるが、いずれの場合も等方性弾性体を仮定する。なお、各シミュレーションで設定した弾性係数の値としては、本実施の形態の組織弾性計測システムで主な対象としている乳房組織の弾性係数にほぼ即している。また、組織圧縮前後のシミュレーションＲＦ信号を作成するために、各組織モデルには点散乱体を分布させた。その際、点散乱体の平均密度としては５００個／ｃｍ３とし、組織圧縮前の散乱体の位置は一様乱数により、散乱係数は平均１．０、標準偏差０．３の正規乱数により決めた。そして、組織圧縮後の散乱体位置は有限要素解析の結果に応じて組織圧縮前の各散乱体を移動させることにより決めている。ここで、実際の組織の散乱体に関する情報は未知であるが、シミュレーションＲＦ信号を基にＢモード像にした際、実際の組織のＢモード像に近くなるように各パラメータを設定する。
【００９８】
この実施の形態では、組織モデルに対する組織圧縮前後のシミュレーションＲＦ信号を次式のように組織圧縮前後の散乱体関数と超音波システムの点広がり関数との畳み込みにより作成する。
【数２９】

ここで、ｉ₁ （ｘ，ｙ，ｚ）は組織圧縮前のＲＦ信号、ｉ₂ （ｘ，ｙ，ｚ）は組織圧縮後のＲＦ信号、ｈ（ｘ，ｙ，ｚ）は超音波システムの点広がり関数（インパルス応答）、ｔ₁ （ｘ，ｙ，ｚ）は組織圧縮前の散乱体関数、ｔ₂ （ｘ，ｙ，ｚ）は組織圧縮後の散乱体関数である。ただし、散乱体関数とは組織モデル内の散乱体が存在する位置ではその散乱係数の値をとり、その他の位置では０であるような関数である。また、組織圧縮後の散乱体関数ｔ₂ （ｘ，ｙ，ｚ）は組織圧縮前散乱体関数ｔ₁ （ｘ，ｙ，ｚ）の各散乱体の位置を組織モデルの変形に応じて移動させたものである。ただし、組織圧縮に伴う各散乱体位置での変位は有限要素解析により得られる要素節点での変位ベクトルを線形補間することにより求めている。
【００９９】
また、この実施の形態ではシミュレーション超音波システムとして無焦点、かつ減衰のないシステムを仮定する。すなわち、超音波システムの点広がり関数ｈ（ｘ，ｙ，ｚ）は空間的に不変であると仮定する。さらに、点広がり関数は次式のように方向ごとに分離できると仮定する。
【数３０】

ここで、ｈｙ（ｙ）は超音波ビーム方向の点広がり関数、ｈｘ（ｘ），ｈｚ（ｚ）はそれぞれ超音波ビームに直交した方向の点広がり関数である。ただし、ｘ方向は超音波断層面内の方向（横方向）、ｚ方向は超音波断層面に垂直な方向（スライス方向）とする。そして、各方向の点広がり関数は実際の超音波装置により計測したワイヤー・ターゲット（水中に張った直径０．１３ｍｍのワイヤー）からの反射エコー分布を基に作成する。図１６は、超音波中心周波数５．０ＭＨｚの場合に用いた各点広がり関数の一例を示す図である。図１６（Ａ）は軸方向の点広がり関数ｈｙ（ｙ）を示し、これはガウス関数に正弦波をかけたものによって実際のワイヤー・ターゲットからの反射エコー分布を近似し、図１６（Ｂ）は横方向点の広がり関数ｈｘ（ｘ）を、図１６（Ｃ）はスライス方向の広がり関数ｈｚ（ｘ）をそれぞれ示し、これらはガウス関数によって実際のワイヤー・ターゲットからの反射エコー分布を近似する。また、各関数のパラメータは中心周波数に応じて変えており、各シミュレーションの際に改めて説明する。
【０１００】
次に、変位（歪み）分布推定法として本実施の形態の拡張複合自己相関法の有効性をシミュレーションにより評価する。まず、拡張複合自己相関法の複合自己相関法に対する拡張点である組織の横方向変位に対応できる点について検証する。
【０１０１】
図１７は、組織モデルの概略を示す図である。組織モデルは、外形６０ｍｍ×６０ｍｍ（２次元）で、一様な弾性係数分布をもつモデルである。そして、この組織モデルを軸方向に一様な３％の歪みが生じるように圧縮する。ここで、このシミュレーションに関しては拡張複合自己相関法のみの評価を行うため、組織モデルとしては単純な１次元弾性体を仮定している。そして、組織の横方向移動（超音波プローブに対する相対的な横方向移動）に関する影響を検証するため、軸方向の圧縮と同時に横方向に０．０ｍｍから１．４ｍｍまでの横方向変位を与えた。ただし、横方向に関しては単純な平行移動であり、組織に対して超音波プローブが完全に滑った場合を想定している。
【０１０２】
そして、この組織モデルに対して変形前後のシミュレーションＲＦ信号を作成する。その際用いた超音波システムに関するパラメータは、中心周波数５．０ＭＨｚ、パルス幅０．５ｍｍ、超音波ビーム幅１．０ｍｍ、走査ライン間隔０．５ｍｍ、サンプリング周波数３０ＭＨｚである。そして、この圧縮前後のシミュレーションＲＦ信号を用いて、本発明に係る拡張複合自己相関法により歪み分布を推定する。また、比較のために複合自己相関法と空間相関法を用いても歪み分布推定を行った。ここで、単純に精度等を比較できるように各手法では同じサイズの相関窓と探索範囲を用いた。具体的には、拡張複合自己相関法と空間相関法では１．６ｍｍ（軸方向）×２．５ｍｍ（横方向）の２次元相関窓と５．６ｍｍ（軸方向）×７．５ｍｍ（横方向）の２次元探索範囲を用い、１次元処理の複合自己相関法では軸方向だけ同じ１．６ｍｍの１次元相関窓と５．６ｍｍの１次元探索範囲を用いた。
【０１０３】
このようにして各手法により歪み分布を推定した結果、図１８〜図２０のようになる。ここで、図１８は横方向変位に対する各手法の歪み推定誤差を示している。◇は複合自己相関法、□は拡張複合自己相関法、△は空間相関法を示す。ただし、歪み推定誤差ｅとしては次式を用いた。
【数３１】

ここで、ε^_iは推定された歪み、ε_i
は実際の歪み（理想値）、ｉは要素番号、Ｎは要素数である。また、図１９は横方向変位が０．０ｍｍの場合の各手法（複合自己相関法、拡張複合自己相関法、空間相関法）により推定した歪み分布、図２０は横方向変位が０．４ｍｍの場合の各手法（複合自己相関法、拡張複合自己相関法、空間相関法）により推定した歪み分布を示す図である。なお、図１９及び図２０は軸方向の深さごとに推定された歪みの平均値と標準偏差を表している。
【０１０４】
このシミュレーション結果より、１次元の複合自己相関法（ＣＡ法）では組織の相対的な横方向変位が超音波ビームを越えて生じてしまう（今回の場合、横方向変位がビーム幅の片側分である０．５ｍｍを超えてしまう）と歪み推定が急に悪くなってしまうのに対し、２次元の拡張複合自己相関法では横方向変位の大きさにかかわらず安定して歪み推定が可能であることが理解できる。また、空間相関法も横方向変位にかかわらず安定して歪み推定が行えるものの、拡張複合自己相関法と比べると２倍以上誤差が大きく精度が悪いことが理解できる。また、各手法の処理時間を比較したところ、下表のように拡張複合自己相関法は複合自己相関法に比べて３．６倍時間がかかってしまうものの、空間相関法と比べると１／（７．７）の時間しかかからなかった。そのため、拡張複合自己相関法はある程度リアルタイム性が保たれていることが理解できる。
┌−−−−−−−−−−−┬−−−−−−−−−−−┬−−−−−−−−─┐
│ 手 法 │ 処理時間 │ 処理時間比較 │
├−−−−−−−−−−−┼−−−−−−−−−−−┼−−−−−−−−─┤
│ 複合自己相関法 │ ２６秒 │ １／（３．６） │
├−−−−−−−−−−−┼−−−−−−−−−−−┼−−−−−−−−─┤
│ 拡張複合自己相関法 │ １分３４秒 │ １．０ │
├−−−−−−−−−−−┼−−−−−−−−−−−┼−−−−−−−−─┤
│ 空間相関法 │ １２分５秒 │ ７．７ │
└−−−−−−−−−−−┴−−−−−−−−−−−┴−−−−−−−−─┘
次に、斜め方向圧縮に関する検証について説明する。前述のシミュレーションでは簡単な２次元の均一組織モデルを用いたが、次は実際の生体組織と同じ３次元構造を持つ組織モデルを用いてシミュレーションを行う。また、超音波プローブにより組織を圧縮する際、超音波ビーム方向（軸方向）に圧縮するのが理想であるが、ここでは斜めに圧縮してしまった場合の影響を検証する。
【０１０５】
図２１は、斜め方向圧縮の検証に使用される組織モデルの概略を示す図である。組織モデルは、図２１（Ａ）に示すように、外形が６０ｍｍ×６０ｍｍ×６０ｍｍの３次元構造をしており、この組織モデル中に直径１５ｍｍ、長さ６０ｍｍの硬い円柱形内包物が存在するようなモデルである。ここで、周辺の弾性係数（ヤング率）は１０ｋＰａ、内包物の弾性係数は周辺より３倍硬い３０ｋＰａとする。ただし、この弾性係数の値は今回主な対象としている乳房組織の弾性係数および乳がんの弾性係数を基に設定する。そして、この組織モデルを２通りの方法で圧縮を行った。１つ目は、図２１（Ｂ）に示すように、この組織モデルに対して上面から軸方向に一様な２００Ｐａの外部圧力を加えて、組織モデルを軸方向に約２％圧縮する。２つ目は、図２１（Ｃ）に示すように、この組織モデルに対してモデル上面から斜め方向の一様な外部圧力（軸方向に２００Ｐａ、横方向に３０Ｐａ）を加えて、組織モデルを斜め方向に圧縮する。
【０１０６】
そして、上記の２通りの場合についてそれぞれ組織圧縮前後のシミュレーションＲＦ信号を作成し、拡張複合自己相関法により歪み分布を推定する。なお、ここでも比較のために複合自己相関法と空間相関法による歪み分布推定も行う。ただし、単純に比較できるように各手法で用いた相関窓サイズと探索範囲は同じとし、そのサイズは前のシミュレーションと同じとする。また、シミュレーションＲＦ信号作成の際に用いた超音波システムに関するパラメータも前のシミュレーションと同じ、中心周波数５．０ＭＨｚ、パルス幅０．５ｍｍ、横方向超音波ビーム幅１．０ｍｍ、スライス方向超音波ビーム幅２．０ｍｍ、走査ライン間隔０．５ｍｍ、サンプリング周波数３０ＭＨｚとする。
【０１０７】
上記のようにして行ったシミュレーションの結果は、図２２及び図２３に示すようになる。ここで、図２２は組織モデルを単純に軸方向に圧縮した場合の歪み分布推定結果を示し、図２３は組織モデルを斜め方向に圧縮した場合の歪み分布推定結果を示す。なお、各図における理想的な歪み分布とは、各条件で３次元有限要素解析を行って得られた軸方向歪み分布であり、この歪み分布を正解としている。また、図２２及び図２３における結果は組織モデルの中央断面での結果である。ここで、図２３において理想歪み分布が左右対称でないのは斜め方向に圧縮した影響であり、特に今回の場合は図に向かって右斜め下に圧縮したため、図右上の部分の横方向変位が大きくなっている。
【０１０８】
そして、このシミュレーションの結果より、軸方向に圧縮した場合は拡張複合自己相関法と複合自己相関法はほぼ同じ結果となったが、斜め方向に圧縮した場合は横方向変位が大きくなるため、複合自己相関法では推定できなくなる領域が生じてしまうのに対し、拡張複合自己相関法では先のシミュレーション同様に横方向変位の大きさにかかわらず安定して歪み推定が可能であることが改めて確認された。また、空間相関法も横方向変位の大きさには依存しないものの拡張複合自己相関法の結果と比較すると明らかに推定精度が悪いのが見てとれる。そのため、前のシミュレーションと合わせて改めて拡張複合自己相関の有効性が確認された。
【０１０９】
前述の拡張複合自己相関法は、精度良く、かつ高速に組織歪み分布を推定できることができる。そこで、次に組織弾性映像システムの第２段階である歪み分布から弾性係数分布を推定する手法で、実施の形態２の３次元有限要素モデルによる弾性係数分布再構成法についてシミュレーションにより検証を行う。
【０１１０】
実施の形態２の弾性係数分布再構成法の最大の特徴は３次元の力学的つりあい方程式に基づいて弾性係数分布を推定することである。そこで、実施の形態２の手法と手順的には同じであるが、２次元の力学的つりあい方程式を基にしている点だけが異なる実施形態１の２次元再構成法を用いて、実施の形態２の３次元再構成法と比較することにより検証を行う。実施形態１の２次元再構成法では、２次元の平面歪み状態を仮定した弾性方程式および有限要素法を用いて弾性係数分布を推定する。
【０１１１】
そこで、まず組織モデルとして、実際の生体組織と同じ３次元構造を持つ図２４のような２つのモデルを用いる。図２４は、３次元構造を持つ２つの組織モデルの一例を示す図である。図２４（ａ）の内包物モデルは、乳がんを模擬したモデルで、外形１００ｍｍ×１００ｍｍ×１００ｍｍのモデル中に直径２０ｍｍの硬い内包物が存在するもので、内包物の弾性係数は３０ｋＰａ、周辺の弾性係数は１０ｋＰａとする。これらの弾性係数の値は、先ほどのシミュレーションと同じように実際の乳房組織の弾性係数を基に決めている。また、周辺と内包物のポアソン比としては、共に非圧縮性に近いため０．４９で一様とする。そして、図２４（ｂ）の層状モデルは筋肉などの層状のものを模擬したモデルで、外形１００ｍｍ×１００ｍｍ×１００ｍｍのモデル中に厚さ２０ｍｍの硬い層がモデル中に存在するというもので、硬い層の弾性係数は３０ｋＰａ、周辺の弾性係数は１０ｋＰａとする。そして、このモデルの場合もポアソン比は０４９で一様とする。
【０１１２】
そして、図２４（ａ）の内包物モデルの場合はモデル上部から軸方向に１００Ｐａの一様な外部圧力により、図２４（ｂ）の層状モデルの場合はモデル上部から軸方向に１５０Ｐａの一様な外部圧力により、各モデルをコンピュータ上で圧縮する。ここで、２つのモデルで外部圧力の強さを変えているのは、各モデルとも同じ約１％の歪みが生じるようにするためである。そして、これらの組織モデル圧縮前後のシミュレーションＲＦ信号を作成し、拡張複合自己相関法により軸方向歪み分布を推定する。そして、この推定された軸方向歪み分布と圧縮の際の境界条件とから、提案した３次元弾性係数分布再構成法により弾性係数分布を推定する。また、同じ軸方向歪み分布と境界条件を用いて比較のために２次元再構成法によっても弾性係数分布を推定する。ここで、シミュレーションＲＦ信号を作成するために用いた超音波システムのパラメータとしては、中心周波数３．７５ＭＨｚ、パルス幅０．７５ｍｍ、横方向超音波ビーム幅２．０ｍｍ、スライス方向超音波ビーム幅２．０ｍｍ、走査ライン間隔２．０ｍｍである。また、拡張複合自己相関法における相関窓のサイズは３．２ｍｍ（軸方向）×４．０ｍｍ（横方向）、探索範囲は１１．２ｍｍ（軸方向）×１４．０ｍｍ（横方向）とする。さらに、３次元有限要素モデルよる弾性係数分布再構成では２ｍｍ（軸方向）×２ｍｍ（横方向）×５ｍｍ（スライス方向）の直方体要素５００００個を用いて組織モデルを構成する。
【０１１３】
そして、このシミュレーションの結果は、図２５〜図２８に示すようになる。図２５及び図２６は内包モデルにおける各推定結果を示す、図２７及び図２８は層状モデルにおける各推定結果を示す。ただし、３次元再構成法では弾性係数の３次元分布を推定しているが、ここではモデル中央断面での結果のみを示す。これは、２次元再構成法では２次元断面での弾性係数分布しか推定できないので、ここでは比較できる中央断面のみ示す。また、各組織モデルにおける３次元再構成結果と２次元再構成結果を数値的に評価したところ次のような結果が得られた。
┌−−−−−−−−−−┬−−−−−−−┬−−−−−−┬−−−−−−−┐
│ │ 周辺領域に │周辺領域に │モデル中心に │
│ │おける弾性係数│おける標準 │おけるコントラ│
│ │ 誤差［％］ │偏差［％］ │スト誤差［％］│
├−−┬−−−−−−−┼−−−−−−−┼−−−−−−┼−−−−−−−┤
│内包│３次元再構成法│ ３．５ │ １５．５│ １１．０ │
│物モ├−−−−−−−┼−−−−−−−┼−−−−−−┼−−−−−−−┤
│デル│３次元再構成法│ ３０．９ │ １７．９│ ３５．９ │
├−−┼−−−−−−−┼−−−−−−−┼−−−−−−┼−−−−−−−┤
│層状│３次元再構成法│ ８．５ │ ２６．８│ ３．１ │
│モデ├−−−−−−−┼−−−−−−−┼−−−−−−┼−−−−−−−┤
│ル │３次元再構成法│ ２４．９ │ ２２．１│ ４３．５ │
└−−┴−−−−−−−┴−−−−−−−┴−−−−−−┴−−−−−−−┘
ここで、用いた評価用のパラメータは、周辺領域における弾性係数誤差ｅ_s 、周辺領域における標準偏差ＳＤ_s 、モデル中心におけるコントラスト誤差ｅ_c であり、それぞれ次式のように定義する。
【数３２】

ただし、上式におけるΣは周辺領域における総和、Ｅ^ は推定された弾性係数、Ｅは実際の弾性係数、Ｎ_s は周辺領域の要素数、Ｅ~_sは周辺領域における弾性係数の平均値、Ｅ^_cはモデル中心における推定弾性係数、Ｅ_c はモデル中心における実際の弾性係数、Ｅ_s は周辺領域における実際の弾性係数である。
【０１１４】
以上のシミュレーション結果より、実施形態２の３次元有限要素モデルによる弾性係数分布再構成法の方が平面応力状態を仮定した２次元再構成法よりもより正確な弾性係数を推定できることが確認された。ここで、３次元再構成法では弾性係数の値が正確に推定されているのに対し、２次元再構成法では実際の弾性係数の値よりも小さく推定しまっている。これは、前もって予想されていたように２次元状態では２次元考察面に垂直な方向の動きが制限されてしまうためである。そのため、実際の生体組織と同じ３次元を基にした弾性係数分布再構成法が必要であることがこのシミュレーションにより改めて示された。
【０１１５】
次に、上述の拡張複合自己相関法と３次元有限要素モデルによる弾性係数再構成法を実装した超音波診断システムの具体的構成について説明する。図２９は、この超音波診断システムの基本構成を示す図である。この超音波診断システムは、３次元超音波スキャナ２８１、超音波診断装置２８２、パーソナルコンピュータ２８３、パルスモータコントローラ２８４、パルスモータ２８５、圧力計２８６などから構成される。３次元超音波スキャナ２８１は超音波パルスを組織内に送信し、かつ組織からの超音波エコー信号を受信するためのものである。ただし、ここでは３次元有限要素モデルによる弾性係数再構成法を用いるため、組織内の３次元的なデータが必要になる。そこで、この超音波診断システムでは３次元超音波スキャナ２８１は３次元走査が可能な構成となっている。超音波診断装置２８２は３次元超音波スキャナ２８１を制御したり、リアルタイムに超音波Ｂモード画像を表示して計測部位を決定したりするためのもので、従来の超音波診断装置をそのまま用いることができる。この超音波診断装置はフルデジタル化された装置で内部にフレームメモリを持っているため、計測したＲＦ信号を一時的に保存しておくことが可能となっている。パーソナルコンピュータ２８３は、超音波診断装置２８２によって計測されたＲＦ信号を受け取り、組織弾性特性推定のための処理（前述の提案手法の処理）、および処理結果の表示を行うためのものである。パルスモータ２８５は組織圧縮を制御するためのものであり、位置固定が可能なアームの先端に固定されており、かつパルスモータ２８５の可動部分には３次元超音波スキャナ２８１が取り付けてある。そして、パルスモータコントローラ２８４によりこのパルスモータ２８５を制御し、超音波スキャナ２８１を組織表面で上下に動かすことにより数パーセントの微小組織圧縮を精度良く行う。圧力計２８６は弾性係数再構成の際に必要な境界条件である体表上での圧力を測るためのもので、超音波スキャナ２８１と体表との間に置かれる。ただし、ここでは超音波スキャナ２８１で組織圧縮を行った際の体表上での圧力分布は一様であると仮定し、圧力計２８６で計測された値をその圧力値として用いる。
【０１１６】
図３０は、この超音波診断システムで用いた超音波スキャナ２８１の具体的構成を示す図である。３次元超音波スキャナ２８１は、超音波振動子が２次元平面状に並んだ２次元アレイ型ではなく、コンベックス型の２次元走査プローブを水中で機械的にスライス方向に振ることにより３次元走査を行うタイプのものである。
【０１１７】
図２９の超音波診断システムは、乳がん診断を主な対象としているため、超音波スキャナの特性も乳腺用に設定してある。今回用いた超音波スキャナの主な特性としては、超音波中心周波数７．５ＭＨｚ、サンプリング周波数３０ＭＨｚ、走査ライン数７１本、フレーム数４４枚、振動子の振れ角３０°、１回の３次元走査にかかる時間０．５秒となっている。ここで、振動子の振れ角とはコンベックス型のプローブをスライス方向に振る際の振れ角のことであり、フレーム数とはコンベックス型のプローブを１回振る際に取得される走査面（フレーム）の数である。また、水中のワイヤー・ターゲットにより超音波パルスの特性を計測したところ、パルス幅約０．５ｍｍ、横方向ビーム幅約１．５ｍｍ、スライス方向ビーム幅約２．６ｍｍであった。
【０１１８】
図２９の超音波診断システムの弾性計測の動作例について説明する。まず、超音波診断装置２８２のリアルタイムＢモード像を見ながら、アームにつながった３次元超音波スキャナ２８１を動かし、計測を行いたい部位に超音波スキャナ２８１を設定する。この際、超音波スキャナ２８１は３次元走査を行わず（すなわち、コンベックス型のプローブを機械的に振らず）、スキャナの中央面のＢモード像のみを超音波診断装置２８２に表示する。これは、実施の形態の超音波診断装置２８２で３次元走査をするとリアルタイムにＢモードを表示できないためである。従って、リアルタイムにＢモードを表示することのできる超音波診断装置の場合は、３次元走査を行いながら部位の設定を行なうことができる。計測部位に超音波スキャナ２８１を移動させたら、アームの可動部を固定し、超音波スキャナ２８１を固定する。そして、組織圧縮前の３次元ＲＦ信号を計測する。これは、３次元走査用のボタンを押すことにより自動的に３次元走査される。また、１回の３次元走査にかかる時間はわずか０．５秒である。このとき、計測された圧縮前のＲＦデータは超音波装置内のフレームメモリに保存しておく。次に、パルスモータ・コントローラ２８４の圧縮用のボタンを１回押すことにより、アームに取り付けられたパルスモータ２８５が前もって設定しておいた量だけ超音波スキャナ２８１を動かし、超音波スキャナ２８１自身により組織を圧縮する。そして、パルスモータ２８５が止まって組織圧縮を行っている状態で、再び３次元走査用
のボタンを押し、組織圧縮後のＲＦデータを計測する。ここで、組織圧縮後のＲＦデータも圧縮前のＲＦデータと同様に超音波装置２８２内のフレームメモリに保存される。また、圧縮している状態で超音波スキャナ２８１の端に取り付けられた圧力計の圧力を計測しておく。以上で、計測部は終わりで、組織圧縮を解除し、被験者は解放される。
【０１１９】
被検者解放後は、パーソナルコンピュータ２８３から超音波診断装置２８２内のフレームメモリにアクセスし、組織圧縮前後のＲＦデータをパーソナルコンピュータ２８３上のハードディスクに保存する。これは、超音波装置内のフレームメモリは一時的なものであり、１回の計測分しか容量がないためである。パーソナルコンピュータ２８３は、拡張複合自己相関法と３次元有限要素モデルによる弾性係数分布再構成法のプログラムを実行し、組織圧縮前後のＲＦデータから歪み分布及び弾性係数分布を推定する。そして、パーソナルコンピュータ２８３は、表示用のプログラムによってモニタ上にＢモード像、歪み像、弾性係数像を並べて表示する。
【０１２０】
以上説明したように、本発明の各実施形態の超音波診断システム、歪み分布表示方法及び弾性係数分布表示方法によれば、横方向変位に対応して変位分布を推定することができると共に超音波ビーム方向（軸方向）の歪み分布のみから弾性係数分布を再構成することができる。
【０１２１】
【発明の効果】
本発明によれば、生体組織の硬さを定量的に計測するにあたり、２次元変位分布の計算時間を短縮でき、かつ、変位分布の計算精度を向上させることができる。
【図面の簡単な説明】
【図１】 超音波診断装置の原理を説明するための図である。
【図２】 超音波信号（受信エコー信号）をディスプレイに表示する方法を示す図であり、図２（Ａ）はＡモード方式、図２（Ｂ）はＢモード方式、図２（Ｃ）はＭモード方式を示す。
【図３】 超音波プローブの種類を示す図である。
【図４】 超音波診断装置を用いて、組織の硬さに関する情報（組織の弾性特性）を計測する手法（機械的振動下における横波伝播速度からの弾性特性評価）の一例を示す図である。
【図５】 静的圧縮による組織弾性計測方式の具体例及び静的圧縮による組織弾性計測方式の原理を示す図である。
【図６】 空間相関法の原理を示す図である。
【図７】 ドプラ法の原理を示す図である。
【図８】 複合自己相関法の原理を示す図である。
【図９】 複合自己相関法の基本アルゴリズムを実行する回路構成を示すブロック図である。
【図１０】 本発明の一実施の形態である超音波診断システムの概略構成を示すブロック図である。
【図１１】 ３次元複合自己相関法の基本アルゴリズムを示すフローチャート図である。
【図１２】 本発明の他の実施の形態の超音波診断システムに係る３次元複合自己相関法の基本アルゴリズムを示すフローチャート図であり、図１１の処理の一部の詳細を示すフローチャート図である。
【図１３】 図１２のステップＳ１５の高速化された複合自己相関法の詳細を示すフローチャート図である。
【図１４】 本発明の他の実施の形態の超音波診断システムに係る３次元複合自己相関法の基本アルゴリズムを実行する回路構成を示すブロック図である。
【図１５】 シミュレーションの手順の概略を示す図である。
【図１６】 超音波中心周波数５．０ＭＨｚの場合に用いた各点広がり関数の一例を示す図である。
【図１７】 組織モデルの概略を示す図である。
【図１８】 横方向変位に対する各手法の歪み推定誤差を示す図である。
【図１９】 横方向変位が０．０ｍｍの場合の各手法（複合自己相関法、拡張複合自己相関法、空間相関法）により推定した歪み分布を示す図である。
【図２０】 横方向変位が０．４ｍｍの場合の各手法（複合自己相関法、拡張複合自己相関法、空間相関法）により推定した歪み分布を示す図である。
【図２１】 斜め方向圧縮の検証に使用される組織モデルの概略を示す図である。
【図２２】 組織モデルを単純に軸方向に圧縮した場合の歪み分布推定結果を示す。
【図２３】 組織モデルを斜め方向に圧縮した場合の歪み分布推定結果を示す。
【図２４】 ３次元構造を持つ２つの組織モデルの一例を示す図である。
【図２５】 内包モデルにおける各推定結果を示す第１の図である。
【図２６】 内包モデルにおける各推定結果を示す第２の図である。
【図２７】 層状モデルにおける各推定結果を示す第１の図である。
【図２８】 層状モデルにおける各推定結果を示す第２の図である。
【図２９】 超音波診断システムの基本構成を示す図である。
【図３０】 超音波診断システムで用いた超音波スキャナの具体的構成を示す図である。
【符号の説明】
９１ 超音波プローブ
９２ 直交検波器
９３ 複素２次元相関計算部
９４ 横方向変位計算部
９５ 軸方向変位計算部
９６ 横方向歪み計算部
９７ 軸方向歪み計算部
９８ 量子化部
９９ 表示部
１３１ 加圧前直交検波回路
１３２ 加圧後直交検波回路
１３３ 第１相関演算回路
１３４ 第１遅延回路
１３５ 第２遅延回路
１３６ 第３遅延回路
１３７ 第４遅延回路
１３４０〜１３４Ｎ 第２相関演算回路
１３５０〜１３５Ｎ 第１相関係数演算回路
１３６０〜１３６Ｎ 第３相関演算回路
１３７０〜１３７Ｎ 位相差演算回路
１３８０〜１３８Ｎ 第２相関係数演算回路
２８１ ３次元超音波スキャナ
２８２ 超音波診断装置
２８３ パーソナルコンピュータ
２８４ パルスモータコントローラ
２８５ パルスモータ
２８６ 圧力計

Claims (19)

1. 被検体組織との間で超音波信号を送受信する超音波探触子と、
   前記被検体組織に対する前記超音波探触子による圧縮前後に前記超音波探触子により受信された前記超音波エコー信号をそれぞれ直交検波して包絡線信号を作成する直交検波手段と、
   前記圧縮前後における前記被検体組織内の複数の計測点における前記包絡線信号間の相関係数を、各計測点を囲む２次元相関窓を用いて予め設定された２次元探索範囲内で、軸方向及びこの軸方向に直交する横方向にそれぞれ所定値ずつ離散した複数の格子点毎に計算し、前記相関係数が最大となる相関最大位置を求め、前記相関最大位置における前記圧縮前後の前記包絡線信号間の位相差を求める相関演算手段と、
   前記相関演算手段によって求められた前記相関最大位置及び前記位相差に基づいて前記圧縮に伴う前記被検体組織内の前記各計測点の前記軸方向及び前記横方向における変位を求める変位演算手段と、
   前記被検体組織内の前記軸方向及び前記横方向の前記各計測点における変位を空間微分することによって歪み分布情報を求める歪み演算手段と、
   前記歪み分布情報を表示する表示手段とを備え、
   前記変位演算手段は、前記相関演算手段により求められた前記相関最大位置及び該相関最大位置における前記位相差のうち、前記横方向における前記相関最大位置に基づいて前記各計測点の前記横方向の変位を求め、前記相関演算手段により求められた前記軸方向の前記相関最大位置及び該相関最大位置における前記位相差の双方に基づいて前記各計測点の前記軸方向の変位を求めて前記歪み演算手段に出力することを特徴とする超音波診断システム。
2. 請求項１において、前記２次元探索範囲は、前記２次元相関窓のサイズより大きいことを特徴とする超音波診断システム。
3. 請求項１において、前記相関演算手段は、前記軸方向には２分の１波長間隔、前記横方向には前記超音波探触子のライン間隔の格子点でのみ相関係数の計算を行うことを特徴とする超音波診断システム。
4. 請求項３において、前記相関演算手段は、前記直交検波手段によって直交検波された包絡線信号を用いて、前記軸方向に２分の１波長の整数倍だけシフトした直交検波包絡線信号を作成し、前記圧縮前後における前記包絡線信号間の相関係数を各格子点について計算することを特徴とする超音波診断システム。
5. 被検体組織との間で超音波信号を送受信する超音波探触子と、
   前記被検体組織に対する前記超音波探触子による圧縮前後に前記超音波探触子により受信された前記超音波エコー信号をそれぞれ直交検波して包絡線信号を作成する直交検波手段と、
   前記圧縮前後における前記被検体組織内の複数の計測点における前記包絡線信号間の相関係数を、各計測点を囲む２次元相関窓を用いて予め設定された２次元探索範囲内で、軸方向及びこの軸方向に直交する横方向にそれぞれ所定値ずつ離散した複数の格子点毎に計算し、前記相関係数が最大となる相関最大位置を求め、前記相関最大位置における前記圧縮前後の前記包絡線信号間の位相差を求める相関演算手段と、
   前記相関演算手段によって求められた前記相関最大位置及び前記位相差に基づいて前記圧縮に伴う前記被検体組織内の前記各計測点の前記軸方向及び前記横方向における変位を求める変位演算手段と、
   前記被検体組織内の前記軸方向及び前記横方向の前記各計測点における変位を空間微分することによって歪み分布情報を求める歪み演算手段と、
   前記被検体組織を有限個の要素に分割して少なくとも２次元の有限要素モデル化し、そのモデル化の情報と前記歪み分布情報を用いて弾性係数分布情報を演算する弾性係数情報演算手段と、
   前記弾性係数分布情報を表示する表示手段とを備え、
   前記変位演算手段は、前記相関演算手段により求められた前記相関最大位置及び該相関最大位置における前記位相差のうち、前記横方向における前記相関最大位置に基づいて前記各計測点の前記横方向の変位を求め、前記相関演算手段により求められた前記軸方向の前記相関最大位置及び該相関最大位置における前記位相差の双方に基づいて前記各計測点の前記軸方向の変位を求めて前記歪み演算手段に出力することを特徴とする超音波診断システム。
6. 請求項５において、前記相関演算手段は、前記軸方向には２分の１波長間隔、前記横方向には前記超音波探触子のライン間隔の格子点でのみ相関係数の計算を行うことを特徴とする超音波診断システム。
7. 請求項６において、前記相関演算手段は、前記直交検波手段によって直交検波された包絡線信号を用いて、前記軸方向に２分の１波長の整数倍だけシフトした直交検波包絡線信号を作成し、前記圧縮前後における前記包絡線信号間の相関係数を各格子点について計算することを特徴とする超音波診断システム。
8. 請求項５において、前記弾性係数情報演算手段は、前記被検体組織を等方性弾性体及び近非圧縮性と仮定し、前記被検体組織を有限個の直方体要素に分割して３次元有限要素モデル化し、前記各要素内では、弾性係数、応力、歪みは一様であると仮定し、弾性方程式に前記歪み分布情報を用いて弾性係数分布情報を演算することを特徴とする超音波診断システム。
9. 被検体組織に接触する超音波探触子によって取得された超音波エコー信号を入力する第１のステップと、
   前記被検体組織に対する前記超音波探触子による圧縮前後における前記超音波エコー信号をそれぞれ直交検波して包絡線信号を作成する第２のステップと、
   前記圧縮前後における前記被検体組織内の複数の計測点における前記包絡線信号間の相関係数を、各計測点を囲む２次元相関窓を用いて予め設定された２次元探索範囲内で、軸方向及びこの軸方向に直交する横方向にそれぞれ所定値ずつ離散した複数の格子点毎に計算し、前記相関係数が最大となる相関最大位置を求め、前記相関最大位置における前記圧縮前後の前記包絡線信号間の位相差を求める第３のステップと、
   前記相関演算手段によって求められた前記相関最大位置及び前記位相差に基づいて前記圧縮に伴う前記被検体組織内の前記各計測点の前記軸方向及び前記横方向における変位を求める第４のステップと、
   前記被検体組織内の前記軸方向及び前記横方向の前記各計測点における変位を空間微分することによって歪み分布情報を求める第５のステップと、
   前記歪み分布情報を表示する第６のステップとを含んで構成され、
   前記第４のステップは、前記第３のステップで求められた前記相関最大位置及び該相関最大位置における前記位相差のうち、前記横方向における前記相関最大位置に基づいて前記各計測点の前記横方向の変位を求め、前記第３のステップで求められた前記軸方向の前記相関最大位置及び該相関最大位置における前記位相差の双方に基づいて前記各計測点の前記軸方向の変位を求めて前記第５のステップに出力することを特徴とする歪み分布表示方法。
10. 請求項９において、前記第３のステップで、前記軸方向には２分の１波長間隔、前記横方向には前記超音波探触子のライン間隔の格子点でのみ相関係数の計算を行うことを特徴とする歪み分布表示方法。
11. 請求項１０において、前記第３のステップで、前記第２のステップによって直交検波された包絡線信号を用いて、前記軸方向に２分の１波長の整数倍だけシフトした直交検波包絡線信号を作成し、前記圧縮前後における前記包絡線信号間の相関係数を各格子点について計算することを特徴とする歪み分布表示方法。
12. 被検体組織に接触する超音波探触子によって取得された超音波エコー信号を入力する第１のステップと、
    前記被検体組織に対する前記超音波探触子の圧縮前後における前記超音波エコー信号をそれぞれ直交検波して包絡線信号を作成する第２のステップと、
    前記圧縮前後における前記被検体組織内の複数の計測点における前記包絡線信号間の相関係数を、各計測点を囲む２次元相関窓を用いて予め設定された２次元探索範囲内で、軸方向及びこの軸方向に直交する横方向にそれぞれ所定値ずつ離散した複数の格子点毎に計算し、前記相関係数が最大となる相関最大位置を求め、前記相関最大位置における前記圧縮前後の前記包絡線信号間の位相差を求める第３のステップと、
    前記相関演算手段によって求められた前記相関最大位置及び前記位相差に基づいて前記圧縮に伴う前記被検体組織内の前記各計測点の前記軸方向及び前記横方向における変位を求める第４のステップと、
    前記被検体組織内の前記軸方向及び前記横方向の前記各計測点における変位を空間微分することによって歪み分布情報を求める第５のステップと、
    前記被検体組織を有限個の要素に分割して少なくとも２次元の有限要素モデル化し、そのモデル化の情報と前記歪み分布情報を用いて弾性係数分布情報を演算する第６のステップと、
    前記弾性係数分布情報を表示する第７のステップとを含んで構成され、
    前記第４のステップは、前記第３のステップで求められた前記相関最大位置及び該相関最大位置における前記位相差のうち、前記横方向における前記相関最大位置に基づいて前記各計測点の前記横方向の変位を求め、前記第３のステップで求められた前記軸方向の前記相関最大位置及び該相関最大位置における前記位相差の双方に基づいて前記各計測点の前記軸方向の変位を求めて前記第５のステップに出力することを特徴とする弾性係数分布表示方法。
13. 請求項１２において、前記第３のステップで、前記軸方向には２分の１波長間隔、前記横方向には前記超音波探触子のライン間隔の格子点でのみ相関係数の計算を行うことを特徴とする弾性係数分布表示方法。
14. 請求項１３において、前記第３のステップで、前記第２のステップによって直交検波された包絡線信号を用いて、前記軸方向に２分の１波長の整数倍だけシフトした直交検波包絡線信号を作成し、前記圧縮前後における前記包絡線信号間の相関係数を各格子点について計算することを特徴とする弾性係数分布表示方法。
15. 請求項１２において、前記第６のステップで、前記被検体組織を等方性弾性体及び近非圧縮性と仮定し、前記被検体組織を有限個の直方体要素に分割して３次元有限要素モデル化し、前記各要素内では、弾性係数、応力、歪みは一様であると仮定し、弾性方程式に前記歪み分布情報を用いて弾性係数分布情報を演算することを特徴とする弾性係数分布表示方法。
16. 被検体との間で超音波信号を送受信する超音波探触子と、
    前記超音波探触子により受信された受信信号を記憶する記憶手段と、
    前記記憶手段に記憶された前記被検体の圧縮前後のフレームデータに計測点を設定し、該計測点を前記フレームデータに対して前記計測点を囲む２次元相関窓を用いて予め設定された２次元探索範囲内で移動して、前記計測点における前記圧縮前後の受信信号の包絡線信号の相関係数が最大となる相関最大位置及び該相関最大位置における前記圧縮前後の前記包絡線信号間の位相差を求める相関演算手段と、
    前記相関演算手段によって求められた前記相関最大位置及び前記位相差に基づいて前記圧縮に伴う前記計測点の変位並びに前記被検体の組織の歪み分布を求める演算手段と、
    前記歪み分布を表示する表示手段とを備え、
    前記歪み分布を求める演算手段は、前記相関演算手段により求められた前記相関最大位置及び該相関最大位置における前記位相差のうち、前記横方向における前記相関最大位置に基づいて前記各計測点の前記横方向の変位を求め、前記相関演算手段により求められた前記軸方向の前記相関最大位置及び該相関最大位置における前記位相差の双方に基づいて前記各計測点の前記軸方向の変位を求めることを特徴とする超音波診断システム。
17. 請求項１６において、前記相関演算手段は、前記計測点を前記超音波ビーム方向に前記超音波信号の２分の１波長間隔で移動し、前記超音波ビームの走査方向に前記超音波ビームピッチで移動して前記相関係数の最大位置を求めることを特徴とする超音波診断システム。
18. 請求項１６において、前記被検体を有限個の要素に分割して少なくとも２次元の有限要素モデルを作成し、そのモデル化の情報と前記歪み分布を用いて弾性係数分布を演算する弾性係数演算手段を備え、前記表示手段に前記弾性分布を表示することを特徴とする超音波診断システム。
19. 請求項１６において、前記相関演算手段は、前記計測点を囲む２次元相関窓を用いて前記包絡線信号の相関係数が最大となる位置を求めるに際し、圧縮後の前記包絡線信号の自己相関関数を求め、この自己相関関数を前記計測点の移動に合わせてずらして前記相関係数を求めることを特徴とする超音波診断システム。

Images