JP7236312B2 - 超音波診断装置、信号処理装置、および、プログラム - Google Patents

Description

本発明は、超音波診断装置に係り、被検体内においてせん断波を発生させ、その伝搬速度を計測することにより、生体組織の性状を評価する技術に関する。

超音波やＭＲＩ（Magnetic Resonance Imaging）、Ｘ線ＣＴ（Computed Tomography）に代表される医療用の画像表示装置は、目視できない生体内の情報を数値または画像の形態で提示する装置として広く利用されている。中でも超音波を利用して画像を表示する超音波撮像装置は、他の装置と比較して高い時間分解能を備えており、例えば拍動下の心臓を滲みなく画像化できる性能を持つ。

被検体である生体内を伝搬する波は、主に縦波と横波に区別される。超音波撮像装置の製品に搭載されている組織形状を映像化する技術や、血流速度を計測する技術は、主に縦波（音速約１５４０ｍ／ｓ）の情報を利用している。

近年、生体内を伝搬する横波（以降、せん断波）を利用して組織の弾性率を評価する技術が注目されており、慢性肝疾患や癌に対する臨床利用が進められている。この技術では、計測対象となる組織内部にせん断波を発生させ、その伝搬速度から弾性率等の弾性を表す評価指標を算出する。せん断波を発生させる手法は、機械方式と放射圧方式に大別される。機械方式は、バイブレータ等を利用して体表面に１ｋＨｚ程度の振動を与えてせん断波を発生させる方式で、振動源となる駆動装置が必要である。一方、放射圧方式は、超音波を組織内の局所に集中させる集束超音波を利用して生体内に音響放射圧を加え、瞬時的に発生する組織変位を利用してせん断波を発生させる。いずれの方式も、発生したせん断波による組織変位を超音波により計測することにより、伝搬速度を算出する。そして、算出されたせん断波の伝搬速度から、組織性状を表す弾性率等の特性値を計算によって求める。

このように、せん断波を利用して組織の弾性を評価する方法は、弾性を定量的に計測できるため、腫瘍診断において極めて重要であり、臨床価値が高い。しかし、せん断波を用いて組織弾性を計測する際、せん断波が組織構造による反射、屈折、回折および減衰を起こし、計測精度や再現率を低下させ、診断性能を劣化させることが知られている。

例えば、特許文献１には、検査対象を伝搬するせん断波の波面振幅の分布を計測し、このデータをフーリエ変換し、フーリエ空間上でフィルタ処理をすることで、せん断波の計測対象となる主成分と、反射、屈折および回折成分とを分け、主成分を抽出する方法が開示されている。

具体的には、特許文献１の技術では、検査対象のある深さｚにおける、伝搬するせん断波の波面振幅データを、方位方向－時間（ｘ－ｔ面）の２次元平面にプロットし、これを２次元フーリエ変換することにより、空間周波数ｋと時間周波数ｆの２次元のフーリエ空間（ｋ－ｆ面）の２次元平面における強度分布に変換する。せん断波の速度は、ｋ－ｆ面において、波面振幅データの点と原点を結ぶ直線とｋ軸とのなす角度θに比例するので、主成分近傍の速度成分のみを取り出すように、所定の角度範囲の波面振幅データを取り出すフィルタ処理を行う。フィルタ処理後のｋ－ｆ面上の波面振幅データを、フーリエ逆変換することにより、実空間（ｘ－ｔ面）に変換し、速度算出を行なう。この処理をすべての深さｚに対して行なうことにより、ｘ－ｚ面全体のせん断波速度マップを作成することができる。

特開２０１８－９９１８０号公報

上述した放射圧方式によって、せん断波の伝搬速度から、組織性状を表す弾性率を計算によって精度よく求めるためには、伝搬速度を正確に算出する必要がある。例えば、予め伝搬方向を仮定して２以上の計測点を設け、両計測点の通過に要する時間を正確に計測することが考えられるが、せん断波は、計測ラインに対し水平に伝搬するとは限らず、体内の組織構造や音波伝搬の物理的な原理に起因した反射、屈折、回折および減衰などの影響を受けるため、様々な周波数や進行方向成分を含む。このため、反射波等の影響を受けずに、せん断波の速度を精度よく、かつ、再現性よく計測することは難しい。

特許文献１に記載の技術では、フーリエ変換を応用し、フーリエ空間（ｋ－ｆ面）で所定の角度範囲の波面振幅データを取り出すフィルタ処理によりせん断波の主成分を抽出することができる。

しかしながら、ｋ－ｆ面における波面振幅データの点と原点を結ぶ直線とｋ軸とのなす角度θと、伝搬速度との対応関係は、角度θが大きい（せん断波の速度が大きい）領域では、わずかな角度変化によって、対応する伝搬速度が大きく異なるという特性があるため、フィルタを掛けて主成分のみを精度よく抽出するのは容易ではない。

本発明の目的は、せん断波の速度を精度よく計測することのできる超音波診断装置を提供することである。

上記目的を達成するために、本発明によれば、検査対象に生じたせん断波による組織の変位の時間変化データを、検査対象に対して超音波を送信し、反射波を受信することにより得た受信信号から算出する計測部と、変位の時間変化データを、空間周波数と時間周波数を２軸とする周波数空間の変位分布を示すスペクトルデータに変換し、周波数空間のスペクトルデータにフィルタを掛けることにより所定の領域のスペクトルデータを抽出する抽出部と、抽出部が抽出した所定の領域のスペクトルデータに基づいてせん断波の速度を算出する速度算出部とを有する超音波診断装置が提供される。抽出部は、周波数空間においてスペクトルデータを所定の角度だけ回転させるスペクトル回転部を備え、回転後のスペクトルデータに対してフィルタを掛けることにより所定の領域のスペクトルデータを抽出する。

本発明によれば、せん断波の速度を精度よく計測することができる。

本発明の実施形態の超音波診断装置の一構成例のブロック図である。 実施形態１における実施形態の超音波診断装置の動作を示すフローチャートである。 実施形態の超音波の送受信の方法をそれぞれ示す説明図である。 （ａ）はせん断波の伝搬する様子、（ｂ）はせん断波の伝搬を時系列情報、（ｃ）はある深さＺiにおけるｘ－ｔ平面の変位分布、をそれぞれ示す説明図である。 （ａ）はある深さにおけるｘ－ｔ平面の変位分布を、（ｂ）は（ａ）を２次元フーリエ変換することによって得られたｋ－ｆ面のパワースペクトルを、（ｃ）は（ｂ）に対して行うフィルタ処理をしたものを、それぞれ示した図である。 （ａ）はフィルタの第１象限における速度分布の等高線表示したものを、（ｂ）は（ａ）の速度分布から算出した速度勾配分布をベクトル表示したものを、（ｃ）は（ｂ）の速度勾配の絶対値を等高線表示したものを、（ｄ）はフィルタ上の第１象限の位置を明確にするため、フィルタ全体の速度分布を等高線表示したものを、（ｅ）は一つのメッシュ当たりの速度範囲を、それぞれ示す図である。 （ａ）はある深さにおけるｘ－ｔ平面の変位分布を、（ｂ）は（ａ）を２次元フーリエ変換してｋ－ｆ面のパワースペクトルにしたものを、（ｃ）は（ｂ）をαだけ回転させたパワースペクトルを、それぞれ示す図である。 （ａ）は、所定の角度範囲を抽出するフィルタを示す図である。（ｂ）は、フィルタによる抽出半径範囲の決定方法を示す図である。（ｃ）は、（ａ）、（ｂ）で決定された抽出角度、抽出半径に基づき作成されたフィルタを示す図である。 円弧上のピーク探索によって主成分を算出する手法を示した処理フローを示す図である。 （ａ）および（ｂ）は、円弧上のピーク探索によって主成分を算出する手法の概要を示す図である。 Ｒａｄｏｎ変換によって主成分を算出する手法を示した処理フローを示す図である。 （ａ）および（ｂ）は、Ｒａｄｏｎ変換によって主成分を算出する手法の概要を示す図である。 速度勾配の絶対値の分布によってｋ－ｆ面のパワースペクトルの回転角を決める手法の処理フローを示す図である。 （ａ）および（ｂ）は、速度勾配の絶対値の分布によってｋ－ｆ面のパワースペクトルの回転角を決める手法の概要を示す図である。 閾値処理によって、ｋ－ｆ面のパワースペクトルの回転角を決める手法の処理フローを示す図である。 （ａ）および（ｂ）は、閾値処理によって、ｋ－ｆ面のパワースペクトルの回転角を決める手法における、閾値による許容範囲の決定方法を示す図である。 （ａ）および（ｂ）は、閾値処理によって、ｋ－ｆ面のパワースペクトルの回転角を決める手法の概要を示す図である。 （ａ）は回転ｋ－ｆ面のパワースペクトルに対するフィルタ処理を、（ｂ）はフィルタ処理後の回転ｋ－ｆ面のパワースペクトルを、（ｃ）は（ｂ）をフーリエ逆変換し、ｘ－ｔ面の変位分布を、それぞれ示す図である。 フィッティングによるせん断波速度の算出方法の処理フローを示す図である。 （ａ）および（ｂ）は、フィッティングによるせん断波速度の算出方法の概要を示す図である。 実施形態２における実施形態の超音波診断装置の動作を示すフローチャートである。 ｋ－ｆ空間分布の回転角度及びフィルタを生成し、メモリに記憶する処理フローを示す図である。 中心角度を求める手法を示した図である。 フィルタの中心角度を求める手法を示した図である。 超音波診断装置の表示部における表示画面例であって、Ｂ像と、弾性またはせん断波速度マップが表示され、さらにフィルタ処理の前後の波面の勾配分布が表示された画面例を示す図である。 超音波診断装置の表示部における表示画面例であって、Ｂ像と、弾性またはせん断波速度マップが表示され、フィルタ処理前後のＲＯＩ中の速度分布と、主成分、フィルタの速度範囲が表示された画面例を示す図である。 超音波診断装置の表示部における表示画面例であって、ユーザが速度範囲を入力する入力画面領域２７０１－２７０３の表示画面の例を示す図である。 超音波診断装置のメモリに格納されたテーブルの例であり、選択可能な臓器や疾患と、それらに適した速度範囲の関係を示すテーブルの例を示す図である。

＜＜実施形態１＞＞
以下、本発明の実施形態を図面に従い説明する。
＜超音波診断装置の全体構成＞
図１に実施形態の超音波診断装置の一構成例のブロック図を示す。本実施形態の超音波診断装置は、送受信制御部２０と、制御部（信号処理装置）３０とを備えている。また、超音波診断装置には、探触子１０と、外部入力デバイス１３と、表示部１６が接続されている。

送信制御部２０は、探触子１０を構成する各振動子に受け渡す送信信号を生成する送信ビームフォーマ２１と、探触子１０の各振動子の出力から、検査対象１００内の所定の点についての受信信号を生成する受信ビームフォーマ２２を備えている。また、制御部３０は、計測部３１と、フィルタ生成部３２と、主成分抽出部３３と、速度算出部３４と、画像生成部３５を備えている。

探触子１０は、音源となる振動子（素子）が規則的に並べられている。各々の素子に対し、送信ビームフォーマ２１は、所定の遅延時間ずつ遅延させた送信信号を出力する。振動子は、この送信信号によって振動することにより、所望の超音波のビームを形成する。送信された超音波ビームは、検査対象内部で反射等し、振動子１０に戻る。振動子１０は戻ってきた超音波を信号化し、受信ビームフォーマ２２に送る。受信ビームフォーマ２２は、受信走査線上の複数の点について振動子の出力信号を整相加算することにより、受信信号（ＲＦ（Radio Frequency）信号）を生成する。

計測部３１は、ＲＦ信号を用いて、検査対象１００の内部の組織の変位を時系列に計測する。

フィルタ生成部３２は、送受信制御部２０から制御部３０に送られた信号や、送受信制御部２０で用いたパラメータ、外部入力デバイス１３から入力された情報に基づいて、主成分抽出部３３で用いるフィルタを生成する。フィルタ生成部３２には、せん断波の主成分抽出の精度を向上させるためのスペクトル回転部３７が備えられている。

主成分抽出部は、せん断波データにフィルタを適用することにより、せん断波の主成分４０２を抽出する。

速度算出部３４は、主成分抽出部３３で得られたせん断波の主成分から、せん断波の速度を算出する。

画像生成部３５は、速度算出部３４で得たせん断波の速度や、せん断波の速度を弾性率に変換した画像データを作成し、表示部１６に送る。

制御部３０の計測部３１、フィルタ生成部３２、主成分抽出部３３、速度算出部３４、画像生成部３５は、ソフトウェアによって実現することも可能であるし、その一部または全部をハードウェアによって実現することも可能である。ソフトウェアによって実現する場合、制御部３０をＣＰＵ（Central Processing Unit）やＧＰＵ（Graphics Processing Unit）等のプロセッサにより構成し、制御部３０に予め格納されたプログラムを読み込んで実行することにより、計測部３１、フィルタ生成部３２、主成分抽出部３３、速度算出部３４、画像生成部３５の機能を実現する。また、ハードウェアによって実現する場合には、ＡＳＩＣ（Application Specific Integrated Circuit）のようなカスタムＩＣやＦＰＧＡ（Field-Programmable Gate Array）のようなプログラマブルＩＣを用い、計測部３１、フィルタ生成部３２、主成分抽出部３３、速度算出部３４、画像生成部３５の動作を少なくとも実現するように回路設計を行なえばよい。

＜超音波診断装置の各部の動作＞
以下、上述した各部の動作を図２から図５を用いて具体的に説明する。ここでは、制御部３０をソフトウェアにより実現する場合を例に説明する。図２は、装置全体の動作を示すフローチャートである。図３は、探触子１０から送信し、受信される超音波の位置と送信方向を模式的に示す説明図である。図４（ａ）はせん断波による組織の変位分布と伝搬方向をｘ－ｚ面において示し、図４（ｂ）はせん断波による組織の変位分布を時系列に取得することを示し、図４（ｃ）は、図４（ｂ）のせん断波による組織の変位分布を検査対象のある深さＺiにおけるｘ－ｔ面において示す図である。図５（ａ）は、図４（ｃ）のｘ－ｔ面の変位分布の強度の大きな方向とｔ方向のなす角を示し、図５（ｂ）はｘ－ｔ面の変位分布を２次元フーリエ変換した空間周波数－時間周波数（ｋ－ｆ）面での変位分布（パワースペクトル）を示し、図５（ｃ）は、（ｋ－ｆ）面のパワースペクトルにおいて特定の範囲の変位データを抽出するフィルタの例を示している。

（（ステップ２０１））
まず、ステップ２０１において、制御部３０は、送受信制御部２０に指示し、検査対象１００内にせん断波を発生させることのできる強度の第１超音波３０１を探触子１０から送信させる。第１超音波３０１が照射された検査対象１００には、第１超音波３０１の焦点付近に音響放射圧が生じて局所的に圧力が加わるため、焦点を中心にせん断波が発生し、放射状に伝搬する。これにより、検査対象１００内に設定されたＲＯＩ３００内へせん断波を伝搬させることができる。

具体的には、制御部３０は、図３に示した検査対象１００内のＲＯＩ（Region of Interest）３００の位置を送受信制御部２０に指示する。送受信制御部２０は、送信信号を送信ビームフォーマ２１に生成させる。生成された送信信号は、探触子１０を構成する振動子にそれぞれ出力される。これにより、探触子１０から、所定の音響強度を有し、所定の深さの焦点に集束する第１超音波がＲＯＩ３００もしくはその近傍の所定の位置に送信される。対象物１００内の第１超音波３０１の焦点付近には音響放射圧が加わり、第１超音波３０１の送信を停止すると圧力負荷が排除されて復元力が働くため、対象部１００にはせん断波が発生する。せん断波は、第１超音波３０１が照射された位置を基点に放射状に伝搬する。図２では、第１超音波３０１の焦点から右方向に伝搬するせん断波３０４がＲＯＩ３００を通過するように第１超音波３０１を送信している。

検査対象１００が均質で無限に広がっていると仮定すると、放射圧によって発生したせん断波は、図３に示すように、放射圧が加わった方向（深さ方向）４０１に対して垂直な方向に検査対象１００内を伝搬する。しかし、実際のせん断波は、生体内の不均質な組織構造や音波伝搬の物理的な性質により、屈折、反射、回折等が発生するため、図４（ａ）に示すように様々な方向へ伝搬する。結果として、せん断波は、主成分の他に、屈折・反射・回折波４０２等の成分を持ち、これらが重なり合って、様々な周波数成分と速度成分を持つ。そのため、本実施形態では後述するステップ２１０，２１１において、フィルタ生成部３２が、せん断波の様々な成分から、真の計測対象である主成分４０３を抽出するためのフィルタを生成する。主成分抽出部３３は、フィルタ処理後のせん断波データから主成分４０３の速度を算出する。ステップ２１０、２１１の動作については、後で詳しく説明する。

（（ステップ２０２））
次に、ステップ２０２において、制御部３０は、送受信制御部２０に指示し、図３のように、せん断波速度計測用の第２超音波３０２を探触子１０からＲＯＩ３００に送信し、その反射波等を受信させる。第２超音波３０２は、ＲＯＩ３００内のせん断波が伝搬する方向（例えばｘ方向）に例えば等間隔に設定された計測点３０５に向かって順次送信される。具体的には、送受信制御部２０は、送信ビームフォーマ２１に送信信号を生成させ、送信信号は、探触子１０の各振動子に出力される。これにより、探触子１０から複数の計測点３０４に所定のタイミングで順次第２超音波３０２が送信される。

第２超音波３０２は、計測点３０５において反射等され、探触子１０に戻り、探触子１０の振動子によってそれぞれ受信される。送受信制御部２０は、第２超音波３０２の複数の計測点３０５をそれぞれ通り、深さ方向（ｚ方向）に延びる複数の受信走査線を設定する。受信ビームフォーマ２２は、各振動子の受信信号に整相加算等の受信ビームフォーミング処理を行うことにより、受信走査線上に設定された複数の深さｚの点（受信焦点）にそれぞれ焦点を合わせた整相後受信信号を得る。これにより、整相後受信信号が、受信走査線方向に連結されたＲＦ信号が生成される。

送受信制御部２０は、第２超音波３０２の送信と反射波等の受信を、所定の時間間隔で繰り返し、時間経過ごとのＲＦ信号を複数の受信走査線についてそれぞれ生成する。

（（ステップ２０３））
ステップ２０３において、制御部３０の計測部３１は、ＲＦ信号に基づいて、受信走査線上のｚ方向（深さ方向）の複数の受信焦点ごとに変位（せん断波の振幅）を計測する。具体的には、同一の受信走査線について時系列に得たＲＦ信号同士の相互相関演算により、複数の計測点３０５を含む各受信焦点について組織変位を求める。これにより、計測部３１は、ｚ－ｘ面におけるせん断波の変位（振幅）の分布を、時系列に得ることができる（図４（ｂ）参照）。

（（ステップ２０４））
ステップ２０４では、計測部３１は、図４（ｂ）のような、ｚ－ｘ面におけるせん断波の変位（振幅）の分布の時系列な情報から、ある深さＺiにおける変位の時間変化（ｘ－ｔ面の変位分布）を図４（ｃ）のように取得する。

（（ステップ２０５））
ステップ２０５において、フィルタ生成部３２は、図４（ｃ）のｘ－ｔ面の変位分布中に含まれる、分離が容易なノイズを除去する。例えば患者の体動や術者の手ぶれに起因するノイズは超音波の信号に比べ非常に低周波なため、フィルタ生成部３２は、公知のノイズ除去手法により容易に除去することができる。また、機器のシステムノイズのような、恒常的かつ周期的にｘ－ｔ面の変位分布にあらわれるノイズに関しても、フィルタ生成部３２は、公知のノイズ除去手法により、容易に除去することができる。

（（ステップ２０６））
ステップ２０６では、フィルタ生成部３２は、図５（ａ）、（ｂ）に示すように、ｘ－ｔ面の変位分布を、２次元フーリエ変換することにより、空間周波数ｋと時間周波数ｆを２軸とする平面（ｋ－ｆ）面における変位分布（以下、パワースペクトルと呼ぶ）に変換する。

図５（ａ）に示すようにｘ－ｔ面の変位分布においては、区間Δｘにおけるせん断波速度Ｖ_ｓは、パワースペクトルの変位の大きい軸方向５１と時間ｔ軸とのなす角度θにより表される。これを数式で表すと、式（１）のようになる。

これに対し、周波数空間ｋ－ｆ面のパワースペクトルにおいては、速度Ｖ_ｓは、図５（ｂ）に示すように、変位（振幅）の大きい軸方向（主成分４０３の方向）５１と、空間周波数ｋ軸とのなす角度θにより表される。これを数式で表すと、式（２）のようになる。

式（１）、（２）から明らかなようにｘ－ｔ面の変位分布及びｋ－ｆ面のパワースペクトルにおけるせん断波速度は共通の角度θを用いて表すことができる。

（（ステップ２０７～２１３））
フィルタ生成部３２は、ステップ２０７～２１０において、図５（ｃ）に示すように、ｋ－ｆ面のパワースペクトルの変位の大きい軸方向（主成分４０３の方向）５１の周囲にある所定の角度範囲の成分を抽出するフィルタ５０２を生成する。主成分抽出部３３は、ステップ２１１において、フィルタ５０２をｋ－ｆ面のパワースペクトルに適用することにより、主成分４０３とその周囲の成分を抽出し、それ以外の成分を除去する。

（フィルタ５０２による主成分の抽出精度の説明）
ここで、フィルタ生成部３２が生成するフィルタ５０２の主成分４０３の抽出の精度について説明する。

上述の式（２）からわかるように、ｋ－ｆ面のパワースペクトルにおいて、空間周波数ｋ軸と変位の大きい軸方向（主成分４０３の方向）５１とのなす角度θと速度Ｖ_ｓは、線形関係にない。具体的には、速度Ｖ_ｓが大きく（角度θが大きく）なるほど、角度θの変化は小さくなる。よって、ｋ－ｆ面において等速度線６０１を示すと図６（ａ）、（ｄ）のようになり、主成分方向５１と空間周波数ｋ軸とのなす角θが小さい低速域では、等速度線６０１は疎であるが、主成分方向５１と空間周波数ｋ軸とのなす角θが大きい（９０度に近い）高速域では、等速度線６０１は密になる。なお、図６（ａ）は、図５（ｂ）の周波数空間ｋ－ｆにおけるパワースペクトルの第１象限のみの等速度線６０１を示し、図６（ｄ）は、４つの象限の等速度線６０１を示している。

この速度分布の勾配∇・Ｖ_ｓは、式（３）で表され、例えば図６（ｂ）のように図示することができる。また、速度分布の勾配の絶対値は、式（４）のように表され、図６（ｃ）のように図示することができる。

上述したようにフィルタ生成部３２は、図５（ｃ）に示すように、主成分方向５１の周囲にあたる特定の角度成分をフィルタにより抽出することにより、主成分４０３とその周囲の成分を抽出することができる。このとき、図６（ａ）、（ｄ）に示したように、周波数空間ｋ－ｆ面における等速度線６０１は、空間周波数ｋ軸からの角度（すなわち速度）が大きいほど密になる。一方、周波数空間ｋ－ｆ面におけるパワースペクトル（変位分布）は離散信号であり、周波数空間ｋ－ｆ面は一様なサイズのメッシュ（画素）に分割され、メッシュごとに変位の値が割り当てられる。このため、メッシュの位置によって、一つのメッシュの速度分解能が異なる。これを図６（ｅ）を用いて具体的に説明する。

図６（ｅ）は、ｋ－ｆ面を構成するメッシュのうち、空間周波数ｋ軸に近い位置にあるメッシュ６０４と、時間周波数ｆ軸に近い位置にあるメッシュ６０５とを示している。メッシュ６０４とメッシュ６０５のサイズは同一である。なお、図６（ｅ）では、わかりやすくするためにメッシュ６０４，６０５のサイズは、実際より大きく表している。メッシュ６０４は、速度勾配の小さな領域に位置するのに対し、メッシュ６０５は速度勾配の大きな領域に位置するため、メッシュ６０５の速度分解能は、メッシュ６０４の速度分解能よりも低い。すなわち、メッシュ６０４を抽出するフィルタを生成した場合、速度分解能が高く所望の速度を抽出することができるが、メッシュ６０５を選択するフィルタを生成した場合、速度勾配の大きい部分にあるため、広い範囲の速度が抽出され、速度分解能が低下する。

そこで本実施形態では、フィルタ生成部３２は、図１に示すようにスペクトル回転部３７を備え、スペクトル回転部３７は、ｋ－ｆ面におけるパワースペクトルをｋ－ｆ面の原点に対して回転させる。これにより、スペクトル回転部３７は、パワースペクトル（スペクトルデータ）を回転させることにより、抽出したい所定の領域のスペクトルデータ（主成分方向５１のデータ）を、周波数空間ｋ－ｆ面上の速度分解能が高い領域に移動させる。具体的には、パワースペクトルの主成分方向５１を、空間周波数ｋ軸とのなす角θが小さい、速度分解能が高い低速度領域まで移動させる。

（ステップ２０７～２０９）
具体的には、ステップ２０７～２０９において、スペクトル回転部３７は、図７（ａ）のｘ－ｔ面の変位分布を２次元フーリエ変換したｋ－ｆ面におけるパワースペクトル（図７（ｂ））５０１をｋ－ｆ面の原点を中心にして反時計回りに角度α回転させる。これにより、ｋ－ｆ面の速度分解能の低い領域（速度の大きい領域）に存在している主成分方向５１の変位データを、図７（ｃ）のように速度分解能の高い領域（速度の小さい領域）に移動させることができる。ステップ２０７～２０９の各ステップの詳細は後述する。

（ステップ２１０）
フィルタ生成部３２は、ステップ２１０において、図８（ａ）に示すように回転移動後のパワースペクトル７０１の主成分方向５１を中心とする所定角度範囲を抽出するフィルタを生成する。

（ステップ２１１）
主成分抽出部３３は、フィルタ生成部３２が生成したフィルタ５０２を用いて、ｋ－ｆ面のパワースペクトルのフィルタ処理を行なって、主成分４０３およびその周辺の成分の変位データを抽出する。

（ステップ２１２～２１３）
速度算出部３４は、ステップ２１１において抽出された主成分４０３およびその周辺の成分の変位データから速度を算出する。このとき、主成分４０３の軸方向が、ステップ２０８において角度αで回転していることを考慮し、角度αの回転角に相当する速度を除去して、主成分４０３の変位データの速度を算出する。これにより、高い分解能で速度抽出を行なうことができるため、せん断波の計測精度を向上させることができる。

（（ステップ２１４））
制御部３０は、ステップ２０４～２１３の処理を、すべての深さＺについて速度が算出されるまで繰り返し行う（ステップ２１４）。

＜ステップ２０７～２０９の詳細＞
上述のステップ２０７～２０９の詳細な処理について説明する。

（ステップ２０７の詳細）
ステップ２０７において、スペクトル回転部３７は、図７（ｂ）のようなパワースペクトル５０１における主成分の方向５１を求め、求めた主成分の方向とｋ軸との角度θmainを求める。

ステップ２０７の処理の具体例を、図９の処理フロー（ステップ８０１～８０４）と図１０の説明図を用いて説明する。

まず、図９のステップ８０１では、スペクトル回転部３７は、図１０（ａ）に示すように、ｋ－ｆ面のｋ軸となす角度θが０からπ／２の範囲（第１象限）において、半径ｒiの円弧ｃiを設定し、円弧ｃi上の変位データのうちピーク値（最大値）９０１の変位データを探索する。

次に、図９のステップ８０２において、スペクトル回転部３７は、ステップ８０１で求めたピーク値９０１の変位データの位置と原点を結んだ線分と、ｋ軸となす角度θiを算出する。

スペクトル回転部３７は、上記ステップ８０１及び８０２を、探索半径ｒiが予め定めた最大値ｒ_maxとなるまで繰り返す（ステップ８０３）。最大値ｒ_maxは、ｋ－ｆ空間の対角線の半分の長さなど、測定の度に自動的に決定してよいし、予め値を定め、定めた値をメモリに記憶しておき、測定時に呼び出して用いてもよい。

図９のステップ８０３において、探索半径が最大値ｒ_maxに達したならば、スペクトル回転部３７は、ステップ８０４に進み、ステップ８０２で算出した角度θi群の中央値または平均値等の代表となる値を算出し、その値を主成分４０３の角度θmainとする。図１０（ｂ）に、探索された複数のピーク値９０１の位置と、主成分４０３の角度θmainの位置関係が示す。

また、別の手法として、例えば図１１のような処理フロー（ステップ１００１～１００３）により、主成分４０３の角度θmainを求めることもできる。図１２は、図１１のフローに対応する説明図である。

この手法では、まずステップ１００１において、図１２（ａ）のｘ－ｔ面の変位分布をラドン変換する。これにより、ｘ－ｔ面の角度θごとに、その角度θ上にある変位データを投影面ｓに投影（積算）した積算値がθ－ｓ面上に算出される。θ－ｓ面において、最も積算値が大きくなる角度を主成分４０３の角度θmainとする。

（ステップ２０８、２０９の詳細）
次に、ステップ２０８について説明する。ステップ２０８では、スペクトル回転部３７は、図７（ｃ）に示すパワースペクトルの回転角αを求める。回転角αは、例えば図１３に示すような処理フローによって算出される。図１４は、図１３の処理フローに対応する説明図である。

まず、図１３のステップ１２０１～１２０３では、スペクトル回転部３７は、図６（ａ）～（ｃ）に示すような、ｋ－ｆ面のパワースペクトルの速度分布を算出し（ステップ１２０１）、ｋ－ｆ面のパワースペクトル内の速度勾配∇・Ｖ_ｓを算出し（ステップ１２０２）、その絶対値｜∇・Ｖ_ｓ｜を算出する（ステップ１２０３）。

次に、スペクトル回転部３７は、ステップ１２０４において、回転前のｋ－ｆ面のパワースペクトルＧ（ｋ，ｆ）を，回転行列Ｒ（α_tmp）を用いて、式（５）の演算により、角度α_tmp回転させた後、式（６）にように速度勾配の絶対値｜∇・Ｖ_ｓ｜を掛けることにより、コスト関数Ψ_gradを算出する（ステップ１２０５）。

速度勾配の絶対値の小さい部分に、より大きなパワー（変位データ）が存在していれば、より高分解能に速度を抽出できるため、図１３のステップ１２０６では、スペクトル回転部３７は、コスト関数Ψ_gradが最小となる角度α_tmpを０＜角度α_tmp＜π／２の範囲で探索する。求められた角度α_tmpをステップ２０８で使用する回転角αとして決定する。

また、スペクトル回転部３７は、図１５のフローに示す別の方法により、回転角αを決定してもよい。まず。図１５のステップ１４０１～１４０３では、スペクトル回転部３７は、ステップ１２０１～１２０３と同様に、ｋ－ｆ面のパワースペクトル内の速度勾配の絶対値｜∇・Ｖ_ｓ｜の分布を算出する（図１６（ａ））。

次に、ステップ１４０４では、スペクトル回転部３７は、算出した速度勾配の絶対値｜∇・Ｖ_ｓ｜に対し、予め定めた閾値Ｔにより閾値処理を施し、速度勾配の絶対値｜∇・Ｖ_ｓ｜が閾値以下の範囲を許容範囲として設定する（図１６（ｂ））。これにより、スペクトル回転部３７は、許容範囲の領域内のメッシュの値を１とするテーブルＰ（ｋ，ｆ）を生成する（図１７（ａ））。

なお、上述の閾値Ｔは、測定の都度自動的に決めてもよいし、予めメモリに記憶させたものを呼び出してもよい。閾値の自動的な決定方法としては、例えば速度勾配の絶対値の最大値に対して一定の割合を閾値とする方法がある。

次に、ステップ１４０５、１４０６は、式（５）と同様に、回転前のｋ－ｆ面のパワースペクトルＧ（ｋ，ｆ）を，回転行列Ｒ（α_tmp）を用いて、式（５）の演算により、角度α_tmp回転させた後、式（７）のようにテーブルＰ（ｋ，ｆ）を掛けることにより、コスト関数Ψ_thを算出する。

許容範囲内に、より大きなパワー（変位データ）が存在していれば、より高分解能に速度を抽出できるため、図１５のステップ１４０７、１４０８では、スペクトル回転部３７は、コスト関数Ψ_thが最大となる角度α_tmpを０＜角度α_tmp＜π／２の範囲で探索する。求められた角度α_tmpをステップ２０９で使用する回転角αとして決定する。

スペクトル回転部３７は、ステップ２０９において、ステップ２０８で求めた回転角αだけパワースペクトルを回転させる。

＜ステップ２１０の詳細＞
ステップ２１０の詳細について説明する。ステップ２１０において、フィルタ生成部３２は、速度分離フィルタ５０２の生成を行なう。図８は、速度分離フィルタの作成方法を示した図である。図８（ａ）は、速度分離フィルタの抽出角度範囲の決定方法を、図８（ｂ）は、速度分離フィルタの抽出半径の決定方法を、図８（ｃ）は作成されたフィルタを示している。図８（ａ）において、θ_０は、主成分の速度θmainをα回転させた角度である。この角度θ_０を中心に、抽出する速度範囲に対応する角度範囲であるθsup及びθinfを決定する。θsup及びθinfのθ_０に対する角度の差（角度幅）は、測定の都度自動的に決定してもよいし、予めメモリに記憶されたものを呼び出して用いてもよい。次に、フィルタ生成部３２は、図８（ｂ）のように、抽出する半径範囲を決定する。この処理は、ｘ－ｔ面の変位分布全体にバンドパスフィルタをかけることに相当するため、ｋ－ｆ面のパワースペクトル上に、明確な周波数に依存する特性を持ち、かつ抽出角度範囲の指定のみでは除去しきれない不要成分が存在している場合に有効である。これにより、図８（ｃ）のように、図８（ａ）、（ｂ）の両方の条件を考慮したフィルタが生成される。

ステップ２１０においてフィルタ生成部３２がフィルタを作成した後、ステップ２１１において、主成分抽出部３３は、ｋ－ｆ面のパワースペクトルの変位データのフィルタ処理を行って（図１８（ａ）参照）、主成分４０３およびその周辺の成分を抽出する（図１８（ｂ））。

ステップ２１２にて、速度算出部３４は、逆フーリエ変換処理によりｋ－ｆ面のパワースペクトルをｘ－ｔの変位分布へ変換する（図１８（ｃ）参照）。

ステップ２１１、２１２の処理を表式化すると、フィルタ処理前のｘ－ｔ面の変位分布をｇ_before（ｘ、ｔ）、それを２次元フーリエ変換によってｋ－ｆ面のパワースペクトルに変換したものをG（ｋ、ｆ）、回転行列をR（α）、フィルタをH（ｋ、ｆ）、フィルタ処理後のｘ－ｔ面の変位分布をｇ_after（ｘ、ｔ）とすると、ｇ_before（ｘ、ｔ）とｇ_after（ｘ、ｔ）は、式（８）、（９）の関係にある。

＜ステップ２１３の詳細＞
ステップ２１３の詳細について説明する。ステップ２１３において、速度算出部３４は、速度算出を行なう。図１９は、速度算出を行なう処理フローを示したもので、図２０はその説明図である。

まず、図１９のステップ１９０１では、速度算出部３４は、ｘ－ｔの変位分布において、ｘ軸について計測範囲２００１を設定する。

ステップ１９０２では、速度算出部３４は、設定した計測範囲内の各ｘの点において、ｘ－ｔの変位分布の変位がピーク値を示す時間ｔ（ピーク時間２００２）を求める（図２０（ａ）参照）。

次に、ステップ１９０３では、速度算出部３４は、計測範囲２００１において算出した複数のピーク時間２００２に対し、図２０（ｂ）のようにフィッティングを行なう。フィッティングの手法は、例えば、最小二乗法による線形回帰や、Ｍ推定による重みづけ線形回帰がある。線形回帰を行なった際の直線の傾きを求めることにより、速度算出部３４は、速度を求めることができる。

このとき、傾きから、ステップ２０９で回転させた角度αに相当する角度を差し引くことにより、回転前の速度を求めることができる。

ステップ１９０１から１９０４を、計測範囲を全ｘ対して移動させて行なうことにより、各計測範囲について速度を求めることができる。

（ステップ２１４）
ステップ２１４では、ステップ２０４からステップ２１３までの処理をすべてのｚに対して行なう。これによって、ｘ－ｚ平面の、せん断波速度のマップを作成することができる。このせん断波速度のマップは、そのまま表示部１６に表示してもよいし、式（１０）の関係（Ｅは弾性率、ρは媒質の密度）から、弾性率に変換してから表示部１６に表示してもよい。なお、表示方式に関しては実施形態３で詳細に述べる。

上述してきたように、本実施形態１によれば、ｋ－ｆ面のパワースペクトルの主成分の方向を角度α回転させてからフィルタを生成しているため、速度分解能の高いフィルタを生成することができる。よって、フィルタ処理により、主成分を精度よく抽出することができるため、主成分の速度の算出（計測）精度が向上するという効果が得られる。

＜＜実施形態２＞＞
上述した実施形態１では、図２のフローに示したように、深さＺi毎にフィルタを作成し、回転角度αを求めた。この方法は、適応的で正確な処理が可能であるが、計算量が大きく、処理に時間がかかる。

そこで、実施形態２では、せん断波速度の範囲が対象臓器によってある程度わかっていることを利用して、この速度レンジを元に予めフィルタや回転角度αを算出してメモリに格納しておき、図２１に示したように、ステップ２１０３においてフィルタや回転角度を主成分抽出部３３がこれを読み出して使用する。このとき、全ての深さＺiにおいて、同じ形状のフィルタと同じ回転角度αを適用する構成とする。これにより、実施形態２では、高速な処理を実現することができる。図２１のステップ２１０３以外のステップ２１０１～２１１２は、図２の各ステップと同様であるので説明を省略する。

図２２は、フィルタ生成部３２がフィルタと回転角度αを予め算出し、メモリに記憶するまでの処理を示したフローである。

まず、図２２のステップ２２０１では、フィルタ生成部３２は、対象臓器や疾患等によって、せん断波の速度範囲（最高速度V_H及び最低速度V_L）を指定する。これは、ユーザが指定してもよいし、接続されている探触子１０の種類や、ＲＦ信号から自動的に判別したものでもよい。

次に、図２２のステップ２２０２では、フィルタ生成部３２は、ｋ－ｆ面のパワースペクトル上の、中心角度θ_centerを例えば、式（１１）により算出する。

次に、ステップ２２０３では、フィルタ生成部３２は、フィルタの回転角度θ_Ｒを設定する。なお、θ_Ｒは、－π／２から+π／２の範囲で変化させる。次に、ステップ２２０４にて、予め定めておいたフィルタの速度範囲に対応する角度α_Ｌ、α_Ｈ用いて、図２３のように、フィルタ生成部３２は、フィルタを作成する。

ステップ２２０５では、フィルタ生成部３２は、コスト関数Ψ_filterを式（１２）により算出する（図２４参照）。

ここで、H（ｋ，ｆ）は作成したフィルタ、Q（ｋ，ｆ）は、例えば図１４（ａ）に示した速度勾配の絶対値分布を示す。

Q（ｋ，ｆ）として速度勾配の絶対値分布を使用した場合、フィルタ生成部３２は、ステップ２２０７において、コスト関数Ψ_filterが最小となる角度θ_Ｒを、回転後の中心角度θ_０と設定する。次に、ステップ２２０８において、回転角度αを式（１３）により決定する。

によって決定する。次に、ステップ２２０９において、式（１２）のコスト関数を用いて決定した中心角度θ_０のフィルタを作成する。

ステップ２２１０では、これらのフィルタと回転角度をメモリに記憶する。

このように、実施形態２では、予めフィルタや回転角度αを算出してメモリに格納しておき、ステップ２１０３においてフィルタや回転角度αを主成分抽出部３３がこれを読み出して使用するため、速度計測中の計算量を低減でき、少ない計算量で素早く、しかも精度よく主成分の速度を算出することができる。

＜＜実施形態３＞＞
実施形態１や、実施形態２において、図１の表示部１６における表示例を図２５及び図２６に示す。

図２５では、装置の表示部に、Ｂモード画像と弾性マップ２５０１を表示している。弾性マップは、ユーザが指定したＲＯＩ２５０２中に表示される。また、計測した弾性率は、表示部の一部に数値として表示してもよい。

これに加え、図２５の表示例では、処理前のせん断波の波面の勾配分布２５０３と、実施形態１または実施形態２の処理方法により処理後のせん断波の主成分の波面の勾配分布２５０４が示される。ユーザは、処理前と処理後の波面の勾配分布を比較することで、実施形態１または２による処理の効果を確認することができる。

図２６は、図２５と同じように、Ｂモード画像と弾性マップ２５０１を表示している。これに加え、図２６では、ＲＯＩ２５０２中の速度のヒストグラムが表示されている。実施形態１または実施形態２の処理方法により処理前と後のヒストグラムを比較することで、処理による効果を確認することができる。

図２７では、装置の表示部に、ユーザが速度範囲を入力する入力画面領域２７０１－２７０３を表示している。

画面領域２７０１には、フィルタを適用するか否か、選択するボタンが表示されている。このボタンをＯＮにするとフィルタを適用したせん断波速度の計測が行われる。

また、画面領域２７０２には、フィルタにおける最高、最低速度の手動指定を受け付ける画面が表示されている。ユーザは、この部分に最高、最低速度を入力すると、それに応じたフィルタ処理が行われる。

画面領域２７０３は、フィルタで抽出する最高、最低速度を自動指定するために、臓器や疾患の選択を受け付ける画面が表示されている。この領域には臓器や、疾患名が表示されており、ユーザは、診断したい臓器や、疾患を選択することができる。

選択可能な臓器や疾患に適した速度範囲が、予め定められ、図２８に示すようにテーブル形式等で装置内のメモリに保持されている。フィルタ生成部３２は、ユーザが選択した臓器や疾患に応じて、テーブル等から速度範囲を呼び出し、それに応じたフィルタを生成し、このフィルタによりフィルタ処理が行われる。

１０：探触子
１３：外部入力デバイス
１６：表示部
２０：送受信制御部
２１：送信ビームフォーマ
２２：受信ビームフォーマ
３０：制御部
３１：計測部
３２：フィルタ生成部
３３：主成分抽出部
３４：速度算出部
３５：画像生成部
１００：検査対象
３００：ＲＯＩ
３０１：第１超音波
３０２：第２超音波
３０３：せん断波
３０４：計測点
４０１：放射圧の方向
４０２：反射・屈折・回折波
４０３：主成分
４０４：ＲＯＩ
５０１：ｋ－ｆ面のパワースペクトル
５０２：フィルタ
６０１：等速度線
６０２：勾配ベクトル
６０４，６０５：メッシュ
７０１：ｋ－ｆ面の回転後のパワースペクトル
９０１：変位のピーク値
１１０１：投影エネルギー分布
１１０２：探索範囲
２００１：計測範囲
２００２：ピーク時間
２００３：線形回帰直線
２５０１：Ｂ像＋弾性マップまたはせん断波速度マップ
２５０２：ＲＯＩ
２５０３：処理前のせん断波の波面の勾配分布
２５０４：処理後のせん断波の波面の勾配分布
２５０５：計測数値
２６０１：せん断波速度のヒストグラム

Claims (13)

1. 検査対象に生じたせん断波による組織の変位の時間変化データを、前記検査対象に対して超音波を送信し、反射波を受信することにより得た受信信号から算出する計測部と、
   前記変位の時間変化データを、空間周波数と時間周波数を２軸とする周波数空間の変位分布を示すスペクトルデータに変換し、前記周波数空間のスペクトルデータにフィルタを掛けることにより所定の領域のスペクトルデータを抽出する抽出部と、
   前記抽出部が抽出した所定の領域のスペクトルデータに基づいて前記せん断波の速度を算出する速度算出部とを有し、
   前記抽出部は、前記周波数空間において前記スペクトルデータを所定の角度だけ回転させるスペクトル回転部を備え、回転後のスペクトルデータに対してフィルタを掛けることにより前記所定の領域のスペクトルデータを抽出することを特徴とする超音波診断装置。
2. 請求項１に記載の超音波診断装置であって、前記抽出部が抽出する前記所定の領域のスペクトルデータは、前記せん断波の主成分のデータであることを特徴とする超音波診断装置。
3. 請求項１に記載の超音波診断装置であって，前記スペクトル回転部は、前記スペクトルデータを回転させることにより、前記抽出部が抽出する前記所定の領域のスペクトルデータを、前記周波数空間上の速度分解能が高い領域に移動させることを特徴とする超音波診断装置。
4. 請求項１に記載の超音波診断装置であって、前記スペクトル回転部は、前記スペクトルデータの変位のピーク値のデータを、空間周波数軸に近づける方向に、周波数空間の中心を回転中心として前記スペクトルデータを回転させることを特徴とする超音波診断装置。
5. 請求項４に記載の超音波診断装置であって、前記スペクトル回転部は、前記スペクトルデータの変位のピーク値のデータの位置を算出し、前記スペクトルデータを回転させる角度を決定することを特徴とする超音波診断装置。
6. 請求項４に記載の超音波診断装置であって、予め求めておいた回転角が格納されたメモリをさらに有し、
   前記スペクトル回転部は、前記回転角を前記メモリから読み出して、前記回転角だけ前記スペクトルデータを回転させることを特徴とする超音波診断装置。
7. 請求項６に記載の超音波診断装置であって、ユーザから検査対象とする臓器および疾患名の少なくとも一方の選択を受け付ける入力部をさらに有し、
   前記メモリには、前記検査対象の複数の臓器および複数の疾患の少なくとも一方について、それぞれ回転角が格納され、
   前記スペクトル回転部は、前記入力部において選択された臓器および／または疾患に対応する前記回転角を前記メモリから読み出して、前記回転角だけ前記スペクトルデータを回転させることを特徴とする超音波診断装置。
8. 請求項１に記載の超音波診断装置であって、前記検査対象に探触子から超音波を送信させ、探触子が受信した反射波の受信信号を受け取る送受信部をさらに有し、
   前記送受信部は、前記検査対象に超音波を送信することにより前記せん断波を生じさせることを特徴とする超音波診断装置。
9. 請求項１に記載の超音波診断装置であって、前記抽出部は、前記周波数空間において、前記変位がピーク値を示す方向を中心に所定の角度範囲で前記スペクトルデータを抽出するフィルタを用いることを特徴とする超音波診断装置。
10. 請求項１に記載の超音波診断装置であって、前記スペクトル回転部は、前記周波数空間上で前記せん断波の速度勾配を算出することを特徴とする超音波診断装置。
11. 請求項１０に記載の超音波診断装置であって、前記スペクトル回転部は、前記せん断波の速度勾配から、前記回転の角度を決定することを特徴とした超音波診断装置。
12. せん断波が伝搬している検査対象に対して超音波を送信し、反射波を受信した受信信号から得られた変位の時間変化データを、空間周波数と時間周波数を２軸とする周波数空間の変位分布を示すスペクトルデータに変換し、前記周波数空間のスペクトルデータにフィルタを掛けることにより所定の領域のスペクトルデータを抽出する抽出部と、
    前記抽出部が抽出した所定の領域のスペクトルデータに基づいて前記せん断波の速度を算出する速度算出部とを有し、
    前記抽出部は、前記周波数空間において前記スペクトルデータを所定の角度だけ回転させるスペクトル回転部を備え、回転後のスペクトルデータに対してフィルタを掛けることにより前記所定の領域のスペクトルデータを抽出することを特徴とする信号処理装置。
13. コンピュータを、
    せん断波が伝搬している検査対象の組織の変位の時間変化データを、空間周波数と時間周波数を２軸とする周波数空間の変位分布を示すスペクトルデータに変換し、前記周波数空間のスペクトルデータにフィルタを掛けることにより所定の領域のスペクトルデータを抽出する抽出手段と、
    前記抽出手段が抽出した所定の領域のスペクトルデータに基づいて前記せん断波の速度を算出する速度算出手段と
    して機能させるプログラムであって、
    前記抽出手段は、前記周波数空間において前記スペクトルデータを所定の角度だけ回転させ、回転後のスペクトルデータに対してフィルタを掛けることにより前記所定の領域のスペクトルデータを抽出することを特徴とするプログラム。
    
    

Images