以下、添付図面を参照して、本発明を実施するための形態(以下、「実施の形態」という)を説明する。
(実施の形態1)
図1は、本発明の実施の形態1に係る超音波診断装置を備えた超音波診断システムの機能構成を示すブロック図である。同図に示す超音波診断システム1は、被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波内視鏡2と、超音波内視鏡2が取得した超音波信号に基づいて超音波画像を生成する超音波診断装置3と、超音波診断装置3が生成した超音波画像を表示する表示装置4と、を備える。
超音波内視鏡2は、その先端部に、超音波診断装置3から受信した電気的なパルス信号を超音波パルス(音響パルス)に変換して被検体へ照射するとともに、被検体で反射された超音波エコーを電圧変化で表現する電気的なエコー信号に変換して出力する超音波振動子21を有する。超音波内視鏡2は、超音波振動子21をメカ的に走査させるものであってもよいし、超音波振動子21として複数の素子をアレイ状に設け、送受信にかかわる素子を電子的に切り替えたり、各素子の送受信に遅延をかけたりすることで、電子的に走査させるものであってもよい。
超音波内視鏡2は、通常は撮像光学系および撮像素子を有しており、被検体の消化管(食道、胃、十二指腸、大腸)または呼吸器(気管、気管支)へ挿入され、消化管、呼吸器やその周囲臓器(膵臓、胆嚢、胆管、胆道、リンパ節、縦隔臓器、血管等)を撮像することが可能である。また、超音波内視鏡2は、撮像時に被検体へ照射する照明光を導くライトガイドを有する。このライトガイドは、先端部が超音波内視鏡2の被検体への挿入部の先端まで達している一方、基端部が照明光を発生する光源装置に接続されている。
超音波診断装置3は、超音波診断システム1全体を制御する制御部31と、超音波診断装置3の動作に必要な各種情報を記憶する記憶部32と、超音波内視鏡2と電気的に接続され、所定の波形および送信タイミングに基づいて高電圧パルスからなる送信信号(パルス信号)を超音波振動子21へ送信するとともに、超音波振動子21から電気的な受信信号であるエコー信号を受信してデジタルの高周波(RF:Radio Frequency)信号のデータ(以下、RFデータという)を生成、出力する送受信部33と、送受信部33から受信したRFデータをもとにデジタルのBモード用受信データを生成する信号処理部34と、送受信部33から受信したRFデータに対して所定の演算を施す演算部35と、各種画像データを生成する画像処理部36と、キーボード、マウス、タッチパネル等のユーザインタフェースを用いて実現され、各種情報の入力を受け付ける入力部37と、を備える。
制御部31は、超音波診断装置3に接続される超音波内視鏡2の種別を判定する種別特定部311を有する。種別特定部311は、超音波診断装置3に接続された超音波内視鏡2内のメモリに格納されているIDを取得することによって超音波振動子21の種別を特定する。種別特定部311が特定した超音波振動子21の種別に関する情報は、後述する記憶部32の種別情報記憶部321に格納される。なお、種別特定部311は、入力部37が入力を受け付けた超音波内視鏡2の種別名に基づいて超音波振動子21の種別を特定するようにしてもよい。
制御部31は、演算および制御機能を有するCPU(Central Processing Unit)や各種演算回路等を用いて実現される。制御部31は、記憶部32が記憶、格納する情報および超音波診断装置3の作動プログラムを含む各種プログラムを記憶部32から読み出し、超音波診断装置3の作動方法に関連した各種演算処理を実行することによって超音波診断装置3を統括して制御する。
記憶部32は、種別特定部311が特定した超音波振動子21の種別情報を記憶する種別情報記憶部321と、超音波診断装置3に接続可能な超音波振動子21の種別ごとにノイズレベルデータを記憶するノイズレベルデータ記憶部322と、演算部35が行う回帰分析の対象となる周波数帯域の情報を超音波診断装置3に接続可能な超音波振動子21の種別ごとに記憶する帯域情報記憶部323とを有する。
図2は、ノイズレベルデータ記憶部322が記憶するノイズレベルデータを模式的に示す図である。図2において、ドットで示す曲面がノイズレベルデータn(f,L)を与える。ノイズレベルデータn(f,L)は、超音波エコーの周波数f、および超音波振動子21の表面と被検体(反射体)との往復距離L(受信深度の2倍に相当する距離)を離散的な変数とする関数である。ノイズレベルデータn(f,L)は、例えばノイズに相当する電圧Vを基準電圧Vcで除した量の常用対数をとってデシベル値(dB)で表現したものであり、離散的なデジタルデータである。なお、以下では、超音波振動子21の表面と被検体との往復距離Lを、単に距離Lということもある。
ノイズレベルデータ記憶部322は、超音波診断装置3に接続可能な超音波振動子21の種別ごとのノイズレベルデータを記憶している。これらのノイズレベルデータは、例えばあらかじめ工場出荷時に超音波振動子21の種別ごとに測定されて格納される。なお、送受信部33の送信を中止したときに受信したエコー信号をノイズレベルデータとして記憶するようにしてもよい。この場合には、超音波振動子21の種別だけでなく、同一種類の超音波振動子21における個体差や、同一被検体の経時変化によるノイズレベルの差異も考慮した、より正確なノイズレベルデータを取得することができる。
帯域情報記憶部323は、演算部35による回帰分析対象の帯域情報として、超音波診断装置3に接続可能な超音波振動子21の種別ごとの周波数帯域の最小値fminおよび最大値fmaxに関する情報を記憶している。帯域情報は、ノイズレベルデータと同様、あらかじめ工場出荷時に超音波振動子21の種別ごとに測定されて格納される。
記憶部32は、上記以外にも、例えば送受信部33、信号処理部34および演算部35が行う各種処理に必要な情報を記憶している。
記憶部32は、超音波診断装置3の作動方法を実行するための作動プログラムを含む各種プログラムを記憶する。各種プログラムは、ハードディスク、フラッシュメモリ、CD−ROM、DVD−ROM、フレキシブルディスク等のコンピュータ読み取り可能な記録媒体に記録して広く流通させることも可能である。なお、上述した各種プログラムは、通信ネットワークを介してダウンロードすることによって取得することも可能である。ここでいう通信ネットワークは、例えば既存の公衆回線網、LAN(Local Area Network)、WAN(Wide Area Network)などによって実現されるものであり、有線、無線を問わない。
以上の構成を有する記憶部32は、各種プログラム等が予めインストールされたROM(Read Only Memory)、および各処理の演算パラメータやデータ等を記憶するRAM(Random Access Memory)およびHDD(Hard Disk Drive)等を用いて実現される。
送受信部33は、エコー信号を増幅する信号増幅部331を有する。信号増幅部331は、受信深度が大きいエコー信号ほど高い増幅率で増幅するSTC(Sensitivity Time Control)補正を行う。図3は、信号増幅部331が行う増幅処理における受信深度と増幅率との関係を示す図である。図3に示す受信深度zは、超音波の受信開始時点からの経過時間に基づいて算出される量である。図3に示すように、増幅率β(dB)は、受信深度zが閾値zthより小さい場合、受信深度zの増加に伴ってβ0からβth(>β0)へ線型に増加する。また、増幅率β(dB)は、受信深度zが閾値zth以上である場合、一定値βthをとる。閾値zthの値は、被検体から受信する超音波信号がほとんど減衰してしまい、ノイズが支配的になるような値である。より一般に、増幅率βは、受信深度zが閾値zthより小さい場合、受信深度zの増加に伴って単調増加すればよい。なお、図3に示す関係は、予め記憶部32に記憶されている。
送受信部33は、信号増幅部331によって増幅されたエコー信号に対してフィルタリング等の処理を施した後、A/D変換することによって時間ドメインのRFデータを生成し、信号処理部34および演算部35へ出力する。なお、超音波内視鏡2が複数の素子をアレイ状に設けた超音波振動子21を電子的に走査させる構成を有する場合、送受信部33は、複数の素子に対応したビーム合成用の多チャンネル回路を有する。
送受信部33が送信するパルス信号の周波数帯域は、超音波振動子21におけるパルス信号の超音波パルスへの電気音響変換の線型応答周波数帯域をほぼカバーする広帯域にするとよい。それにより、後述する周波数スペクトルの近似処理において、精度のよい近似を行うことが可能となる。
送受信部33は、制御部31が出力する各種制御信号を超音波内視鏡2に対して送信するとともに、超音波内視鏡2から識別用のIDを含む各種情報を受信して制御部31へ送信する機能も有する。
信号処理部34は、RFデータに対してバンドパスフィルタ、包絡線検波、対数変換など公知の処理を施し、デジタルのBモード用受信データを生成する。対数変換では、RFデータを基準電圧Vcで除した量の常用対数をとってデシベル値で表現する。このBモード用受信データでは、超音波パルスの反射の強さを示す受信信号の振幅または強度が、超音波パルスの送受信方向(深度方向)に沿って並んでいる。図4は、超音波振動子21の走査領域(以下、単に走査領域ということもある)とBモード用受信データとを模式的に示す図である。図4に示す走査領域Sは扇形をなしている。これは、超音波振動子21がコンベックス振動子である場合に相当している。また、図4では、Bモード用受信データの受信深度をzとして記載している。超音波振動子21の表面から照射された超音波パルスが受信深度zにある反射体で反射し、超音波エコーとして超音波振動子21へ戻ってきた場合、その往復距離Lと受信深度zとの間には、上述したようにz=L/2の関係がある。信号処理部34は、生成したBモード用受信データを、画像処理部36のBモード画像データ生成部361へ出力する。信号処理部34は、CPU(Central Processing Unit)や各種演算用の回路等を用いて実現される。
一方、演算部35は、送受信部33が出力したRFデータに対して受信深度によらず増幅率を一定とするよう増幅補正を行う増幅補正部351と、増幅補正を行ったRFデータに高速フーリエ変換(FFT:Fast Fourier Transform)を施して周波数解析を行うことにより受信深度ごとの周波数スペクトルを算出する周波数解析部352と、ノイズレベルデータと周波数スペクトルのデータとを比較することにより演算対象とする周波数帯域を設定する帯域設定部353と、周波数解析部352が算出した周波数スペクトルにより定義される関数の距離変化率および周波数変化率を所定の順序で算出することによって2次変化率を算出する変化率算出部354と、変化率算出部354が算出した2次変化率を用いて走査領域における超音波パルスの単位距離および単位周波数あたりの減衰率を推定する減衰率推定部355と、を有する。演算部35は、CPUや各種演算用の回路等を用いて実現される。なお、演算部35を、制御部31および信号処理部34と共通のCPU等を用いて構成することも可能である。
図5は、増幅補正部351が行う増幅補正処理における受信深度と増幅率との関係を示す図である。図5に示すように、増幅補正部351が行う増幅処理における増幅率β(dB)は、受信深度zがゼロのとき最大値βth−β0をとり、受信深度zがゼロから閾値zthに達するまで線型に減少し、受信深度zが閾値zth以上のときゼロである。なお、図5に示す関係は、予め記憶部32に記憶されている。増幅補正部351が図5に示す関係に基づいてRFデータを増幅補正することにより、信号増幅部331におけるSTC補正の影響を相殺し、一定の増幅率βthの信号を出力することができる。なお、増幅補正部351が行う増幅補正処理における受信深度zと増幅率βの関係は、信号増幅部331が行う増幅補正処理における受信深度と増幅率の関係に応じて異なることは勿論である。
このような増幅補正を行う理由を説明する。STC補正は、アナログ信号波形の振幅を全周波数帯域にわたって均一に、かつ、受信深度の増加に対して単調増加する増幅率で増幅させることで、アナログ信号波形の振幅から減衰の影響を排除する補正処理である。このため、エコー信号の振幅を輝度に変換して表示するBモード画像を生成する場合、かつ、一様な組織を走査した場合には、STC補正を行うことによって深度によらず輝度値が一定になる。すなわち、Bモード画像の輝度値から減衰の影響を排除する効果を得ることができる。一方、本実施の形態1のように超音波の周波数スペクトルを算出して解析した結果を利用する場合、STC補正でも超音波の伝播に伴う減衰の影響を正確に排除できるわけではない。何故なら、減衰量は周波数によって異なるが、STC補正の増幅率は距離だけに対して変化し、周波数に対しては変化せず一定であるためである。
上述した問題、すなわち超音波の周波数スペクトルを算出して解析した結果を利用する場合、STC補正でも超音波の伝播に伴う減衰の影響を正確に排除できるわけではない、という問題を解決するには、Bモード画像を生成する際にSTC補正を施した受信信号を出力する一方、周波数スペクトルに基づいた画像を生成する際に、Bモード画像を生成するための送信とは異なる新たな送信を行い、STC補正を施していない受信信号を出力することが考えられる。ところがこの場合には、受信信号に基づいて生成される画像データのフレームレートが低下してしまうという問題がある。
そこで、本実施の形態1では、生成される画像データのフレームレートを維持しつつ、Bモード画像用にSTC補正を施した信号に対してSTC補正の影響を排除するために、増幅補正部351によって増幅率の補正を行う。
周波数解析部352は、増幅補正部351が増幅補正した各音線のRFデータ(ラインデータ)を所定の時間間隔でサンプリングし、サンプルデータを生成する。そして、周波数解析部352は、サンプルデータ群にFFT処理を施すことにより、RFデータ上の複数の箇所(データ位置)における周波数スペクトルを算出する。
図6は、超音波信号の1つの音線におけるデータ配列を模式的に示す図である。同図に示す音線SRkにおいて、白または黒の長方形は、1つのサンプル点におけるデータを意味している。また、音線SRkにおいて、右側に位置するデータほど、超音波振動子21から音線SRkに沿って計った場合の深い箇所からのサンプルデータである(図6の矢印を参照)。音線SRkは、送受信部33が行うA/D変換におけるサンプリング周波数(例えば50MHz)に対応した時間間隔で離散化されている。図6では、番号kの音線SRkの8番目のデータ位置を受信深度zの方向の初期値Z(k) 0として設定した場合を示しているが、初期値の位置は任意に設定することができる。周波数解析部352による算出結果は複素数で得られ、記憶部32に格納される。
図6に示すデータ群Fj(j=1,2,・・・,K)は、FFT処理の対象となるサンプルデータ群である。一般に、FFT処理を行うためには、サンプルデータ群が2のべき乗のデータ数を有している必要がある。この意味で、サンプルデータ群Fj(j=1,2,・・・,K−1)はデータ数が16(=24)で正常なデータ群である一方、サンプルデータ群FKは、データ数が12であるため異常なデータ群である。異常なデータ群に対してFFT処理を行う際には、不足分だけゼロデータを挿入することにより、正常なサンプルデータ群を生成する処理を行う。この点については、後述する周波数解析部352の処理を説明する際に詳述する(図15を参照)。
周波数解析部352は、RFデータから切り取ったサンプルデータ群の各々に対してFFT処理を施すことによって電圧振幅の周波数成分V(f,L)を生成する。この電圧振幅の周波数成分V(f,L)は電圧の周波数密度である。さらに、周波数解析部352は、電圧振幅の周波数成分V(f,L)を基準電圧Vcで除し、常用対数(log)をとってデシベル単位で表現する対数変換処理を施した後、適当な定数Aを乗ずることにより、次式(2)で与えられる周波数スペクトルのデータ(以下、スペクトルデータともいう)F(f,L)を生成し、帯域設定部353へ出力する。
F(f,L)=A・log{V(f,L)/Vc} ・・・(2)
ここで、logは常用対数である(以下、同じ)。
スペクトルデータF(f,L)は、サンプルデータ群の周波数fの成分である。図7は、記憶部32が記憶するスペクトルデータのデータ列を模式的に示す図である。図7において、縦方向は周波数fを示し、横方向は超音波振動子21の表面からの往復距離Lを示している。周波数fは、離散的な値0、Δf、2Δf、・・・を取る。例えば、列0のセルには、距離区間0≦L<ΔLで切り取られたサンプルデータ群に基づいて式(2)から得られるスペクトルデータF(f,0)が格納される。また、列ΔLのセルには、距離Lの区間ΔL≦L<2ΔLで切り取られたサンプルデータ群に基づいて式(2)から得られるスペクトルデータF(f,ΔL)が格納される。図7では、周波数Δf、距離ΔLのセルのみにスペクトルデータF(Δf,ΔL)のみを例示的に記載しているが、実際には全てのセルに、各セルの周波数および距離に応じたスペクトルデータが格納されていることはいうまでもない。なお、サンプルデータ群を切り取る区間の長さΔL(図6のサンプルデータ群のステップ幅Dに相当)は、例えば1.0cm程度である。また、周波数の変化量Δfは、例えば0.5MHzである。
図8は、スペクトルデータの具体例を示す図である。図8では、互いに異なる4つの距離におけるスペクトルデータF(f,L1)、F(f,L2)、F(f,L3)、F(f,L4)と周波数fとの関係を示している。ここで、4つの距離L1、L2、L3、L4は定数であり、0<L1<L2<L3<L4を満たしている。図7からも明らかなように、実際にはさらに多くのスペクトルデータF(f,L)が算出されるが、図8では、代表的な4つのスペクトルデータのみを例示している。図8に示すように、スペクトルデータF(f,L)は、距離Lが大きくなるにつれて減少していく。また、スペクトルデータF(f,L)の平均周波数は、距離Lが大きくなるにつれて低周波数側へシフトしていく。これは、被検体内での超音波が伝播する際の周波数に依存した減衰の効果によるものである。
一般に、スペクトルデータF(f,L)は、超音波が走査された組織の属性によって異なる傾向を示す。これは、スペクトルデータF(f,L)が、超音波を散乱する散乱体の大きさ、数密度、音響インピーダンス等と相関を有しているためである。ここでいう「属性」とは、例えば悪性腫瘍組織、良性腫瘍組織、内分泌腫瘍組織、粘液性腫瘍組織、正常組織、嚢胞、脈管などのことである。
帯域設定部353は、スペクトルデータF(f,L)とノイズレベルデータn(f,L)とを比較する比較部353aを有する。比較部353aは、上述した二つのデータを比較するために、超音波診断装置3に接続されている超音波振動子21に応じて、後述する回帰分析の対象とする周波数帯域U={f|fmin≦f≦fmax}の両端の周波数fmin、fmaxを帯域情報記憶部323から読み出す。周波数帯域Uは、超音波振動子21の表面(L=0)における超音波の送信波形上で比較的平坦な区間に相当し、最小値fminおよび最大値fmaxは、超音波振動子21の種別に応じて異なる。
比較部353aは、帯域情報記憶部323から読み出した周波数帯域Uにおいて、スペクトルデータF(f,L)と、ノイズレベルデータn(f,L)とを周波数f,距離Lごとにそれぞれ比較する。図8に示す4つのスペクトルデータF(f,Lp)(p=1,2,3,4)を対応するノイズレベルデータn(f,Lp)と比較した場合、例えば、p=1,2,3のとき、周波数帯域Uに含まれる任意の周波数fでF(f,Lp)>n(f,Lp)が成立する一方、p=4のとき、周波数帯域Uでは以下の2つの不等式が成立したとする。
F(f,L4)>n(f,L4) (fmin≦f<fmax’)
F(f,L4)≦n(f,L4) (fmax’≦f≦fmax)
図9は、スペクトルデータF(f,L4)とノイズレベルデータn(f,L4)との関係を示す図である。図9に示す場合、2つの周波数帯域f≦fmin’、f≧fmax’においてF(f,L4)≦n(f,L4)が成立し、ノイズレベルデータが支配的である。以下、この2つの周波数帯域のようにノイズレベルデータが支配的である周波数帯域を、ノイズ周波数帯域という。図9より、周波数fmin’、fmax’のうち、周波数帯域Uに含まれるのはfmax’である。
帯域設定部353は、比較部353aの比較結果に基づいて回帰分析対象の周波数帯域を設定する。図8および図9に示す場合、帯域設定部353は、距離Lp(p=1、2、3)に対する回帰分析対象を初期の周波数帯域Uと設定する一方、距離L=L4の回帰分析対象の周波数帯域をU’={f|fmin≦f<fmax’}と設定し、これらの周波数帯域情報を減衰率推定部355へ出力する。
図10は、帯域設定部353が設定した回帰分析対象の周波数帯域を視覚的に表現したスペクトルデータのデータ列を模式的に示す図である。図10においても、各セルのスペクトルデータF(f,L)の記載は省略している。図10では、スペクトルデータF(f,L)の種類に応じてセルごとに3種類の異なる模様のいずれかを付与している。白色のセルに格納されるスペクトルデータF(f,L)は、周波数fが周波数帯域Uに含まれるとともに、F(f,L)>n(f,L)を満たす値を有することを示している。また、斜線が記載されたセルに格納されるスペクトルデータF(f,L)は、周波数fが周波数帯域Uの外の周波数であるとともに、F(f,L)>n(f,L)を満たす値を有することを示している。さらに、ドットが記載されたセルに格納されるスペクトルデータF(f,L)は、周波数fが周波数帯域Uの外の周波数であるとともに、F(f,L)≦n(f,L)を満たす値を有することを示している。上述したように、距離L=Lp(p=1,2,3)における周波数帯域はUであり、距離L=L4における周波数帯域はU’である。なお、図10に記載されたLmax(fp)(p=1,2,3,4)については後述する。
変化率算出部354は、スペクトルデータF(f,L)により定義される数であって周波数fの1次関数である関数ψ(f,L)=F(f,L)−F(f,0)に対して周波数fおよび距離Lに対する変化率を、周波数から距離の順で順次算出する。
減衰率推定部355は、変化率算出部354が算出した2次変化率を用いて、走査領域における超音波パルスの単位距離および単位周波数あたりの減衰率を推定する。
ここで、変化率算出部354および減衰率推定部355が行う処理をより詳細に説明する。まず、スペクトルデータF(f,L)の単位距離および単位周波数あたりの減衰率の推定方法を説明する。距離Lに存在する反射体からの超音波の周波数fにおける音圧振幅P(f,L)は、正の定数μを用いて、
P(f,L)=P(f,0)・exp(−μfL) ・・・(3)
で与えられることが知られている。μfL>0なので、式(3)は、音圧振幅P(f,L)が、周波数fおよび距離Lの増加に対して指数関数的に減衰することを意味している。
一方、距離区間L〜L+ΔLの周波数fにおける超音波の減衰量をLoss(f,L)[dB]とすると、この減衰量は、
Loss(f,L)=A・log{P(f,L)/P(f,L+ΔL)}
=A・logP(f,L)−A・logP(f,L+ΔL) ・・・(4)
で定義される。ここで、右辺の定数Aは、式(2)の定数Aと同じである。したがって、単位距離および単位周波数あたりの減衰率ζは、次式で与えられる。
ζ=(∂/∂f)Lim{Loss(f,L)/ΔL}
=(∂/∂f){−A(∂/∂L)logP(f,L)}
=−A(∂2/∂f∂L)logP(f,L) ・・・(5)
ここで、Lim{Loss(f,L)/ΔL}は、関数Loss(f,L)/ΔLのΔL→0における極限を意味する。単位距離および単位周波数あたりの減衰率ζの単位は、例えば[dB/cm/MHz]である。以下、単位距離および単位周波数あたりの減衰率を単に減衰率ということもある。
なお、上記定数μと減衰率ζとの関係は次の通りである。式(5)のP(f,L)に式(3)を代入すると、減衰率ζは、
ζ=−A(∂2/∂f∂L)[log{P(f,0)・exp(−μfL)}]
=−A(∂2/∂f∂L){logP(f,0)−μfLloge}
=(loge)Aμ ・・・(6)
ここで、eは自然対数の底である。
さて、超音波振動子21の感度を周波数fの関数としてγ(f)とすると、RFデータにFFT処理を施した後の振幅成分V(f,L)は、次式(7)で与えられる。
V(f,L)=γ(f)・P(f,L) ・・・(7)
この式(7)のP(f,L)に式(3)を代入すると、
V(f,L)=γ(f)・P(f,0)・exp(−μfL)
=V(f,0)・exp(−μfL) ・・・(8)
が得られる。
式(8)を式(2)へ代入することにより、
F(f,L)=A・log{V(f,0)・exp(−μfL)/Vc}
=Alog・exp(−μfL)+Alog{V(f,0)/Vc}
=−(loge)AμfL+F(f,0) ・・・(9)
が得られる。さらに、式(9)の右辺に式(6)を代入することにより、
F(f,L)−F(f,0)=−ζfL ・・・(10)
が導出される。
式(10)の両辺に2階偏微分演算子∂2/∂L∂f、および∂2/∂f∂Lをそれぞれ作用させると、次式が得られる。
ζ=−∂2F(f,L)/∂L∂f=−∂2F(f,L)/∂f∂L ・・・(11)
ここで、∂2/∂L∂fは周波数fの偏微分を先に行う一方、∂2/∂f∂Lは距離Lの偏微分を先に行うことを意味する。
したがって、スペクトルデータF(f,L)の2次偏導関数∂2F(f,L)/∂f∂Lまたは∂2F(f,L)/∂L∂fを演算することにより、減衰率ζを推定することができる。
以上説明した減衰率の推定方法において、スペクトルデータF(f,L)の偏微分を演算することは実際には困難である場合が多い。何故なら、偏微分の定義に従うと、偏微分の演算の際には、極限Δf→0、ΔL→0(Δf、ΔLはそれぞれf、Lの微小変位)を計算する必要があるが、実際のスペクトルデータF(f,L)は離散的に定義され、これらの極限を計算することは困難であるからである。この問題を回避するために、周波数fや距離Lの隣接する離散値の差分をとることによってスペクトルデータF(f,L)の偏微分の演算を近似する方法が知られている。ところがこの方法では、偏導関数がスペクトルデータF(f,L)のゆらぎに起因するノイズを多く含んでしまう可能性がある。
本実施の形態1において、変化率算出部354は、スペクトルデータF(f,L)の関数ψ(f,L)=F(f,L)−F(f,0)に対して回帰分析を実施し、回帰直線による近似を実施する。この関数ψ(f,L)の2次偏導関数は、
∂2ψ(f,L)/∂L∂f=∂2F(f,L)/∂L∂f ・・・(12)
となるので、式(12)の右辺に式(11)を代入すると、
∂2ψ(f,L)/∂L∂f=−ζ ・・・(13)
が得られる。同様に、
∂2ψ(f,L)/∂f∂L=−ζ ・・・(14)
も得られる。式(13)、(14)は、関数ψ(f,L)を用いて減衰率ζを算出できることを示している。
ところで、関数ψ(f,L)は、式(10)により、
ψ(f,L)=−ζfL ・・・(15)
と表される。ところで、回帰直線による近似は、結局、一次関数での近似である。よって、近似する関数が一次関数に近いほど、回帰直線はその近似する関数に近づき、良い近似を与える。ここで、式(15)に示す通り、関数ψ(f,L)は周波数fの一次関数そのものである。一方、スペクトルデータF(f,L)は周波数fの一次関数に近いとは限らない。したがって、周波数fに対する関数の偏導関数を関数の回帰直線の傾き(すなわち変化率)で近似する場合、関数ψ(f,L)を用いる方が、スペクトルデータF(f,L)を用いるよりも、近似の精度が向上することとなる。
変化率算出部354は、関数ψ(f,L)の周波数fに対する偏導関数∂ψ(f,L)/∂fの近似値として、回帰分析により関数ψ(f,L)の周波数fに対する変化率(すなわち、回帰直線の傾き)を算出する。続いて、変化率算出部354は、関数ψ(f,L)の周波数fに対する変化率に対し、さらに回帰分析(第2の回帰分析)を行うことによって距離Lに対する変化率(すなわち、第2の回帰直線の傾き)を算出し、この値を2次偏導関数∂2ψ/∂L∂fの近似値とする。以下、この2次偏導関数∂2ψ/∂L∂fの近似値を2次変化率という。
変化率算出部354のさらに具体的な処理を説明する。まず、変化率算出部354は、帯域設定部353が設定した周波数帯域で、関数ψ(f,L)の周波数fを変数とする回帰直線の傾きおよび切片を得る。
図11は、例として距離L=Lpにおける関数ψ(f,Lp)と周波数fの関係、および各関数の周波数fに対する回帰直線を示す図である(p=1,2,3,4)。回帰直線J1、J2、J3は、変化率算出部354が周波数帯域Uで回帰分析を行うことによって算出した回帰直線である。また、回帰直線J4’は、変化率算出部354が周波数帯域U’で回帰分析を行うことによって算出した回帰直線である。なお、図11では、比較のため、変化率算出部354が、距離L=L4において周波数帯域Uで回帰分析を行うことによって算出した回帰直線J4も記載している。以下、回帰直線Jpの傾きをSf(Lp)とする。また、回帰直線J4’の傾きをSf’(L4)とする。
回帰直線Jpの傾きSf(Lp)は、p=1,2,3において、距離Lの増加とともに単調に減少する。これに対し、回帰直線J4の傾きSf(L4)は、回帰直線J3の傾きSf(L3)よりも大きい(Sf(L4)>Sf(L3))。一方、距離L=L4の周波数帯域U’における傾きSf’(L4)は、回帰直線J3の傾きSf(L3)よりも小さい(Sf’(L4)<Sf(L3))。これは、図8から明らかなように、関数ψ(f,L4)が周波数fmax’の近傍で極小値を取り、周波数fmaxでは極小値より大きい値を取るからである。
続いて、変化率算出部354は、傾きSf(L)を距離Lの関数とみて第2の回帰分析を行うことにより、距離Lに対する第2の回帰直線を算出する。図12は、傾きSf(L)と距離Lの関係を示す図である。変化率算出部354は、上述の通り、例に挙げた傾きSf(L1)、Sf(L2)、Sf(L3)、Sf’(L4)とそれ以外の全ての傾きを算出した。そして、変化率算出部354は、これらの傾きに基づき、さらに往復距離Lに対する第2の回帰分析を行うことで、第2の回帰直線Q1を算出する。この第2の回帰直線Q1を図12に実線で示す。ところで、変化率算出部354は、上述の通り、周波数帯域Uにおいて、傾きSf(L1)、Sf(L2)、Sf(L3)、Sf’(L4)とそれ以外の全ての傾きを算出した。比較のために、これらの傾きに基づき、さらに往復距離Lに対する第2の回帰分析を行うことによって算出した第2の回帰直線Q2を図12に破線で示す。第2の回帰直線Q1、Q2を比較すると、第2の回帰直線Q1の方が、傾きSf(L)の値によくフィットしていることがわかる。式(15)に示した通り、関数ψ(f,L)は周波数fの一次関数であり、かつ、本実施の形態1では周波数帯域U、U’を適切に設定してノイズの影響を排除したため、回帰直線J1、J2、J3、J4’の傾き(すなわち1次変化率)Sf(L1)、Sf(L2)、Sf(L3)、Sf’(L4)は、周波数fに対する関数ψ(f,L)の1次偏導関数∂ψ(f,L1)/∂f、∂ψ(f,L2)/∂f、∂ψ(f,L3)/∂f、∂ψ(f,L4)/∂fにそれぞれ良い近似を与える。さらに、式(15)に示す通り、関数ψ(f,L)は周波数fだけでなく往復距離Lの一次関数でもあるため、第2の回帰直線Q1の傾き(すなわち2次変化率)は、1次偏導関数∂ψ(f,L)/∂fの往復距離Lに対する偏導関数、すなわち式(13)の左辺の2次偏導関数∂2ψ(f,L)/∂L∂fに良い近似を与える。
減衰率推定部355は、変化率算出部354が算出した2次変化率の値を式(13)へ代入することにより、超音波の音線ごとの減衰率ζを算出する。続いて減衰率推定部355は、算出した全ての音線の減衰率ζの平均値を算出し、算出結果を走査領域の減衰率として、画像処理部36が有する合成画像データ生成部362へ出力する。なお、減衰率推定部355は、全ての音線の減衰率ζの最頻値、中央値または最大値等の統計量を走査領域の減衰率としてもよい。
画像処理部36は、エコー信号の振幅を輝度に変換して表示する超音波画像であるBモード画像データを生成するBモード画像データ生成部361と、減衰率推定部355が推定した減衰率ζの情報とBモード画像データを合成して合成画像データを生成する合成画像データ生成部362と、を有する。
Bモード画像データ生成部361は、信号処理部34からのBモード用受信データに対してゲイン処理、コントラスト処理等の公知の技術を用いた信号処理を行うとともに、表示装置4における画像の表示レンジに応じて定まるデータステップ幅に応じたデータの間引き等を行うことによってBモード画像データを生成する。Bモード画像は、色空間としてRGB表色系を採用した場合の変数であるR(赤)、G(緑)、B(青)の値を一致させたグレースケール画像である。
Bモード画像データ生成部361は、Bモード用受信データに走査領域を空間的に正しく表現できるよう並べ直す座標変換を施し、さらにBモード用受信データ間の補間処理を施して、Bモード用受信データ間の空隙を埋め、デジタルデータであるBモード画像データを生成する。Bモード画像データは図4で扇型に示された走査領域内の器官の状況を表現できるBモード画像のデジタルデータである。Bモード画像データ生成部361は、生成したBモード画像データを合成画像データ生成部362へ出力する。
合成画像データ生成部362は、走査領域の減衰率の値を示す文字データを生成し、Bモード画像に隣接して表示するようにBモード画像データと文字データとを合成することによって合成画像データを生成する。合成画像データ生成部362は、生成した合成画像データを表示装置4へ出力する。
表示装置4は、液晶または有機EL(Electro Luminescence)等からなるモニタを用いて構成される。表示装置4は、超音波診断装置3が生成した合成画像データに対応する合成画像を含む各種情報を表示する。
図13は、以上の構成を有する超音波診断装置3が実行する処理の概要を示すフローチャートである。具体的には、超音波診断装置3が超音波内視鏡2からエコー信号を受信する以降の処理の概要を示すフローチャートである。以下、図13を参照して、超音波診断装置3が行う処理を説明する。まず、超音波診断装置3は、超音波内視鏡2から超音波振動子21による測定結果としてのエコー信号を受信する(ステップS1)。
超音波振動子21からエコー信号を受信した信号増幅部331は、そのエコー信号の増幅を行う(ステップS2)。ここで、信号増幅部331は、例えば図3に示す増幅率と受信深度との関係に基づいてエコー信号の増幅(STC補正)を行う。この際、信号増幅部331におけるエコー信号の各種処理周波数帯域は、超音波振動子21による超音波エコーのエコー信号への音響電気変換の線形応答周波数帯域をほぼカバーする広帯域にするとよい。これも、後述する周波数スペクトルの近似処理において精度のよい近似を行うことを可能とするためである。
ところで、上述の通り、信号増幅部331はエコー信号を増幅し、送受信部33は増幅されたエコー信号にフィルタリング、A/D変換を施してRFデータを生成し、信号処理部34はRFデータに各種処理を施してBモード用受信データを生成した。ここで、Bモード画像データ生成部361は、信号処理部34から入力されるBモード用受信データに適切な座標変換、補間処理を施し、Bモード画像データを生成して、表示装置4へ出力する(ステップS3)。Bモード画像データを受信した表示装置4は、そのBモード画像データに対応するBモード画像を表示する。
一方、増幅補正部351は、送受信部33から出力されたRFデータに対して受信深度によらず増幅率が一定となるように増幅補正を行う(ステップS4)。ここで、増幅補正部351は、例えば図5に示す増幅率と受信深度との関係が成立するように増幅補正を行う。
この後、周波数解析部352は、増幅補正後の各音線のRFデータに対してFFTによる周波数解析を行うことによってスペクトルデータを算出する(ステップS5)。このステップS5の処理の詳細については後述する。
続いて、帯域設定部353は、回帰分析対象の周波数帯域を設定する(ステップS6)。例えば、図8および図9に示す場合、L=L1,L2,L3に対して周波数帯域Uを設定する一方、L=L4に対して周波数帯域U’を設定する。
変化率算出部354は、帯域設定部353が設定した周波数帯域に基づいて、関数ψ(f,L)の2次偏導関数∂2ψ(f,L)/∂L∂fの近似値である2次変化率を、2回の回帰分析を行うことによって算出する(ステップS7)。例えば、変化率算出部354は、2次変化率として、図12に示す回帰直線Q1の傾きを算出する。
この後、減衰率推定部355は、走査領域における超音波パルスの減衰率を推定する(ステップS8)。減衰率推定部355は、変化率算出部354が算出した2次変化率の値を式(13)の左辺へ代入することにより、音線ごとの減衰率を算出した後、算出した全ての音線の減衰率の平均値を算出し、この平均値を走査領域の減衰率として合成画像データ生成部362へ出力する。
合成画像データ生成部362は、Bモード画像データおよび走査領域の減衰率をもとに合成画像データを生成し、表示装置4へ出力する(ステップS9)。合成画像データを受信した表示装置4は、その合成画像データに対応する合成画像を表示する。図14は、表示装置4が表示する合成画像の表示例を示す図である。同図に示す合成画像101は、Bモード画像表示部102と減衰率表示部103とを有する。なお、図14では、具体的なBモード画像の表示を省略している。
ステップS9の後、超音波診断装置3は一連の処理を終了する。なお、超音波診断装置3は、ステップS1〜S9の処理を周期的に繰り返し実行する。
次に、ステップS5の周波数解析処理について、図15のフローチャートを参照して説明する。まず、周波数解析部352は、解析対象の音線を識別するカウンタkをk0とする(ステップS11)。この初期値k0は、術者等のユーザが入力部37を通じて任意に指示入力した値、もしくは、記憶部32にあらかじめ設定された値とする。
続いて、周波数解析部352は、上述の通りFFT演算用に生成した一連のデータ群(サンプルデータ群)を代表するデータ位置(受信深度に相当)Z(k)の初期値Z(k) 0を設定する(ステップS12)。例えば、図6では、上述したように、音線SRkの1番目のデータ位置を初期値Z(k) 0として設定した場合を示している。
その後、周波数解析部352は、サンプルデータ群を取得し(ステップS13)、取得したサンプルデータ群に対し、記憶部32が記憶する窓関数を作用させる(ステップS14)。このようにサンプルデータ群に対して窓関数を作用させることにより、サンプルデータ群が境界で不連続になることを回避し、アーチファクトが発生するのを防止することができる。ステップS14で適用する窓関数は、例えばHamming、Hanning、Blackman等のいずれかであり、予め記憶部32に記憶されている。
続いて、周波数解析部352は、データ位置Z(k)のサンプルデータ群が正常なデータ群であるか否かを判定する(ステップS15)。図6を参照した際に説明したように、サンプルデータ群は、2のべき乗のデータ数を有している必要がある。以下、サンプルデータ群のデータ数を2n(nは正の整数)とする。本実施の形態1では、データ位置Z(k)が、できるだけZ(k)が属するサンプルデータ群の中心になるよう設定される。具体的には、サンプルデータ群のデータ数は2nであるので、Z(k)はそのサンプルデータ群の中心に近い2n/2(=2n-1)番目の位置に設定される。この場合、サンプルデータ群が正常であるとは、データ位置Z(k)より浅部側に2n-1−1(=Nとする)個のデータがあり、データ位置Z(k)より深部側に2n-1(=Mとする)個のデータがあることを意味する。図6に示す場合、サンプルデータ群Fj(j=1,2,・・・,K−1)は正常である。なお、図6ではn=4(N=7,M=8)の場合を例示している。
ステップS15における判定の結果、データ位置Z(k)のサンプルデータ群が正常である場合(ステップS15:Yes)、周波数解析部352は、後述するステップS17へ移行する。
ステップS15における判定の結果、データ位置Z(k)のサンプルデータ群が正常でない場合(ステップS15:No)、周波数解析部352は、不足分だけゼロデータを挿入することによって正常なサンプルデータ群を生成する(ステップS16)。ステップS15において正常でないと判定されたサンプルデータ群(例えば図6のサンプルデータ群FK)は、ゼロデータを追加する前に窓関数が作用されている。このため、サンプルデータ群にゼロデータを挿入してもデータの不連続は生じない。ステップS16の後、周波数解析部352は、後述するステップS17へ移行する。
ステップS17において、周波数解析部352は、サンプルデータ群を用いてFFT演算を行うことにより、振幅の周波数分布であるスペクトルデータを得る(ステップS17)。この結果、例えば図7の各列に示すようなスペクトルデータが得られる。
続いて、周波数解析部352は、データ位置Z(k)をステップ幅Dで変化させる(ステップS18)。ステップ幅Dは、記憶部32が予め記憶しているものとする。図6では、D=15の場合を例示している。ここで、図7に示した往復距離Lの間隔ΔLは、このステップ幅Dを距離に換算したときの値(=サンプリング幅×D)の2倍で定義される。したがって、ステップ幅Dが定まると、間隔ΔLも一意に定まる。ステップ幅Dは、Bモード画像データ生成部361がBモード画像データを生成する際に利用するデータステップ幅と一致させることが望ましいが、周波数解析部352における演算量を削減したい場合には、ステップ幅Dとしてデータステップ幅より大きい値を設定してもよい。
その後、周波数解析部352は、データ位置Z(k)が音線SRkにおける最大値Z(k) maxより大きいか否かを判定する(ステップS19)。データ位置Z(k)が最大値Z(k) maxより大きい場合(ステップS19:Yes)、周波数解析部352はカウンタkを1増加させる(ステップS20)。これは、処理をとなりの音線へ移すことを意味する。一方、データ位置Z(k)が最大値Z(k) max以下である場合(ステップS19:No)、周波数解析部352はステップS13へ戻る。
ステップS20の後、周波数解析部352は、カウンタkが最大値kmaxより大きいか否かを判定する(ステップS21)。カウンタkがkmaxより大きい場合(ステップS21:Yes)、周波数解析部352は一連の周波数解析処理を終了する。一方、カウンタkがkmax以下である場合(ステップS21:No)、周波数解析部352はステップS12に戻る。この最大値kmaxは、術者等のユーザが入力部37を通じて任意に指示入力した値、もしくは、記憶部32にあらかじめ設定された値とする。
このようにして、周波数解析部352は、関心領域内の(kmax−k0+1)本の音線の各々について複数回のFFT演算を行う。
以上説明した本発明の実施の形態1によれば、周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率をこの順序で算出することによって得られる2次変化率を用いて、超音波振動子の走査領域における超音波信号の単位距離および単位周波数あたりの減衰率を推定するため、超音波の減衰率を精度よくかつ簡便に算出することができるとともに、減衰率に基づく画像の信頼性を向上させることができる。
また、本実施の形態1によれば、上述した引用文献1のように音速を得る必要はなく、送信波形がガウシアンであることを前提としていないため、減衰率を精度よく算出することができる。なお、本実施の形態1の2次変化率は、周波数と距離(または受信深度)の関数の2次変化率であり、上述した特許文献1における「位相の2次変化率」とは全く異なるものである。この点は、以下に説明する実施の形態2、3も同様である。
また、本実施の形態1によれば、スペクトルデータに応じてノイズレベルデータとの比較により演算対象の周波数帯域を設定するため、S/Nが十分高く、減衰率の推定に有効な領域のみで減衰率を算出することができる。したがって、高精度で減衰率を算出することができ、減衰率に基づく画像の信頼性を高めることができる。
また、本実施の形態1によれば、画像内の領域をマニュアルで指定するなどの煩雑な処理は不要であるため、減衰率を簡便に算出することができる。
また、本実施の形態1によれば、線型の回帰分析を用いるため、周波数変化率および距離変化率を簡便に算出することができる。
また、本実施の形態1によれば、周波数の1次関数であるスペクトルデータの関数ψ(f,L)を用いているため、回帰直線の傾きを用いた近似の精度を向上させることができる。
また、本実施の形態1によれば、ノイズレベルデータ記憶部は複数の種別または機体ごとの超音波振動子にそれぞれ対応するノイズレベルデータを記憶しているため、接続可能な全ての超音波振動子に対して超音波の減衰率を精度よく算出することができる。
なお、本実施の形態1において、比較部353aがスペクトルデータF(f,L)と比較する対象をノイズレベルデータn(f,L)そのものではなく、ノイズレベルデータn(f,L)の1次関数a・n(f,L)+b(a≧1、b≧0;a、bは定数)に代えてマージンを設けてもよい。この場合には、変化率算出部354が算出する2次変化率や、減衰率推定部355が推定する減衰率へのノイズの影響がさらに少なくなり、減衰率の算出精度が一層向上する。
また、一般にS/Nが悪いのは、周波数に依存する減衰の激しい高周波数側であることに鑑み、帯域設定部353が回帰分析対象の周波数帯域を設定する際、最小の周波数を初期値fminで固定する一方、最大の周波数の値をf>fminの範囲で初期値fmaxから変更するようにしてもよい。この場合には、変化率算出部354および減衰率推定部355が演算する際、減衰の激しい高周波数側を除外して減衰率の推定精度を高めることができる。
また、超音波診断装置3は、走査領域の減衰率を推定する代わりに、走査領域の一部の領域の減衰率を推定するようにしてもよい。この場合の領域は、入力部37を介してユーザが設定できるようしておけば好ましい。
(実施の形態2)
本発明の実施の形態2は、超音波診断装置の変化率算出部が行う2次変化率の算出方法が、上述した実施の形態1と異なる。本実施の形態2に係る超音波診断装置は、実施の形態1で説明した超音波診断装置3と同様の構成を有する。
本実施の形態2において、変化率算出部354は、回帰分析によりスペクトルデータF(f,L)の距離Lに対する変化率をまず算出して偏導関数∂ψ(f,L)/∂Lの近似値とした後、第2の回帰分析により偏導関数∂F(f,L)/∂Lの周波数fに対する変化率を算出することによって2次偏導関数∂2F(f,L)/∂f∂Lの近似値すなわち2次変化率を算出する。
この場合、距離L=0のスペクトルデータF(f,0)は距離Lの関数ではないため、
∂F(f,L)/∂L=∂ψ(f,L)/∂L
が成り立つ。したがって、本実施の形態2において、変化率算出部354は関数ψ(f,L)を算出する必要はなく、スペクトルデータF(f,L)から2次変化率を算出することができる。
変化率算出部354は、周波数fごとに回帰分析限界Lmax(f)を抽出する。回帰分析限界Lmax(f)は、周波数fにおける距離Lの最大値である。図10に示すスペクトルデータのデータ列では、周波数f=fp(p=1,2,3,4)に対する回帰分析限界Lmax(fp)を示している。図10に示す場合、回帰分析限界Lmax(fp)は、周波数fpにおける白色セルの分布範囲の右端のセルが有する距離Lの値に対応している。
変化率算出部354は、スペクトルデータF(f,L)の距離Lに対する回帰分析を実行し、距離区間0≦L≦Lmax(f)における回帰直線の傾きSL(f)を算出する。図16は、周波数fを一定にしたときのスペクトルデータF(f,L)と距離Lの関係を示す図である。具体的には、互いに異なる4つの距離におけるスペクトルデータF(f,L1)、F(f,L2)、F(f,L3)、F(f,L4)と距離Lとの関係を示している。ここで、周波数fpは定数であり、0<f1<f2<f3<f4を満たす。図16に示すように、スペクトルデータF(f,L)は、被検体内で超音波が伝搬する際の周波数依存減衰の効果により、周波数fが大きいほど距離Lの増加に伴う減衰が激しい。なお、一般にこれより多くのスペクトルデータF(f,L)が算出されるが、図16では図8と同様に、代表的な4つのスペクトルデータのみを例示している。
図16では、距離区間0≦L≦Lmax(fp)におけるスペクトルデータF(fp,L)の回帰直線Kpが示されている。スペクトルデータF(fp,L)は、L=0から回帰分析限界Lmax(fp)に達するまではほぼ線型であるが、回帰分析限界Lmax(fp)より大きくなると、ノイズの影響およびスペクトルデータF(fp,L)自体が0の近傍まで減衰することによる影響によって線型性が崩れている。
回帰直線Kpの傾きは、周波数依存減衰の効果により、周波数fが大きくなるほど急峻になる。換言すると、回帰直線Kpの傾きをSL(fp)とすると、SL(f1)>SL(f2)>SL(f3)>SL(f4)が成立する。
ここで比較のため、図16に示す4つのスペクトルデータF(fp,L)について、4つの周波数のうち最低の周波数f1の回帰分析限界Lmax(f1)を用いて回帰直線を算出する場合を考える。図16では、一例として、スペクトルデータF(f4,L)に対して回帰分析限界Lmax(f1)を適用した場合に算出される回帰直線K4’を例示している。この回帰直線K4’の傾きSL’(f4)は、スペクトルデータF(f4,L)がL=Lmax(f4)の近傍でノイズレベルまで落ち切ってしまうため、SL(f4)よりも大きい(SL’(f4)>SL(f4))。なお、図示はしないが、p=1,2,3の場合にも、同様の関係が成り立っている。
変化率算出部354は、傾きSL(f)を周波数fの関数とみて第2の回帰分析を行うことにより、周波数fに対する第2の回帰直線を算出する。図17は、傾きSL(f)と周波数fとの関係を示す図である。図17において、黒点は、変化率算出部354が距離区間0≦L≦Lmax(fp)でスペクトルデータF(fp,L)の回帰分析を行うことによって算出した傾きSL(fp)と周波数fpとの関係を示している。これに対して白点は比較用のデータであり、変化率算出部354が距離区間0≦L≦Lmax(f1)でスペクトルデータF(fp,L)の回帰分析を行うことによって算出した回帰直線Kp’の傾きSL’(fp)と周波数fpとの関係を示している。
変化率算出部354は、上述の通り、距離区間0≦L≦Lmax(f)において、例に挙げた傾きSL(f1)、SL(f2)、SL(f3)、SL(f4)と、それ以外の全ての傾きを算出した。そして、変化率算出部354は、これらの傾きに基づき、さらに周波数fに対する第2の回帰分析を行うことで、第2の回帰直線Q3を算出する。この第2の回帰直線Q3を図17に実線で示す。ところで、比較のために、変化率算出部354が、仮に距離区間0≦L≦Lmax(f1)において、傾きSL’(f1)、SL’(f2)、SL’(f3)、SL’(f4)とそれ以外をの全ての傾きを算出し、これらの傾きに基づき、さらに周波数fに対する第2の回帰分析を行った場合の第2の回帰直線Q4を図17に破線で示す。第2の回帰直線Q3、Q4を比較すると、第2の回帰直線Q3の方が、傾きSL(f)の値によくフィットしていることがわかる。第2の回帰直線Q3の傾きは、式(14)の左辺の2次偏導関数∂2ψ(f,L)/∂f∂Lに対する近似値、すなわち2次変化率を与える。
以上説明した本発明の実施の形態2によれば、周波数スペクトルにおける周波数変化率および距離変化率をこの順序で算出することによって得られる2次変化率を用いて、超音波振動子の走査領域における超音波信号の単位距離および単位周波数あたりの減衰率を推定するため、実施の形態1と同様、超音波の減衰率を精度よくかつ簡便に算出することができるとともに、減衰率に基づく画像の信頼性を向上させることができる。
また、本実施の形態2によれば、距離に対する回帰分析を周波数に対する回帰分析よりも先に実行するため、スペクトルデータをそのまま回帰分析することができ、計算量を少なくすることができる。
(実施の形態3)
図18は、本発明の実施の形態3に係る超音波診断システムの機能構成を示すブロック図である。同図に示す超音波診断システム11は、超音波内視鏡2と、超音波診断装置12と、表示装置4とを備える。
超音波診断装置12は、記憶部121および画像処理部122の構成が、上述した超音波診断装置3の記憶部32および画像処理部36とそれぞれ異なる。
記憶部121は、種別情報記憶部321、ノイズレベルデータ記憶部322、帯域情報記憶部323に加えて、減衰率の値に応じて画像に付与する視覚情報を記憶する視覚情報記憶部324を有する。ここでいう視覚情報とは、例えば輝度、色相、明度または彩度等のいずれかであり、減衰率の値に応じた値が対応付けられている。なお、視覚情報記憶部324が複数種類の視覚情報を減衰率と対応づけて記憶しておいてもよい。この場合には、ユーザが入力部37を介して所望の視覚情報を選択できるようにしておけばよい。
画像処理部122は、Bモード画像データ生成部361、合成画像データ生成部362に加えて、減衰率マップデータ生成部363を有する。減衰率マップデータ生成部363は、減衰率推定部355が推定した減衰率の値に応じた視覚情報を、視覚情報記憶部324を参照して画像に付与することによって減衰率マップデータを生成する。
図19は、以上の構成を有する超音波診断装置12が行う処理の概要を示すフローチャートである。具体的には、超音波診断装置12が超音波内視鏡2からエコー信号を受信する以降の処理の概要を示すフローチャートである。図19において、ステップS31〜ステップS36は、図13のステップS1〜S6に順次対応している。以下、ステップS36に続く処理を説明する。
ステップS36の後、変化率算出部354は、スペクトルデータの2次変化率を算出する(ステップS37)。この際、変化率算出部354は、走査領域内であらかじめ分割して設定された複数の部分領域の各々に対して2次変化率を算出する。なお、本実施の形態3において、変化率算出部354が2次変化率を算出する際の周波数変化率および距離変化率の算出順序は、特に限定されない。
図20は、部分領域の設定例を模式的に示す図である。同図に示す部分領域Rは、深度方向の長さ(深度幅)がwで、NR本の音線を含む扇形領域である。走査領域Sでは、深度幅wで送受信方向を分割するとともに、NR本の音線ごとに走査方向が分割されている。図20では、走査領域S内に複数ある送受信方向のうち、1本の送受信方向dr上で算出されたスペクトルデータのうち部分領域R内のデータを黒丸で示す一方、部分領域Rの境界に位置するスペクトルデータを白丸で示している。
変化率算出部354は、部分領域Rにおけるスペクトルデータにおける距離の基準位置(距離がゼロの位置)を、部分領域R内で超音波振動子21の表面に最も近い位置とし、この基準位置との深度の差z’の2倍(2z’)を新たな距離として演算を行う。なお、超音波振動子21の表面を距離の基準位置としたときの部分領域Rの基準位置がL=Lminである場合、上述した白丸の位置におけるスペクトルデータはF(f,Lmin)と表される。変化率算出部354は、このスペクトルデータF(f,Lmin)を式(10)のスペクトルデータF(f,0)の代わりとして、部分領域Rにおける2次変化率を算出するようにしてもよい。
減衰率推定部355は、変化率算出部354の算出結果を用いて、走査領域に含まれる複数の部分領域の各々に対し、部分領域の減衰率を推定する(ステップS38)。減衰率推定部355は、まず部分領域で算出された全ての2次変化率を用いて音線ごとの減衰率を算出する。この後、減衰率推定部355は、同一の部分領域で算出した全ての音線の減衰率の平均値を算出し、この平均値を当該部分領域の減衰率の推定値として減衰率マップデータ生成部363へ出力する。なお、減衰率推定部355は、同一の部分領域で算出した全ての音線の減衰率の最頻値、中央値または最大値等の統計量を当該部分領域の減衰率としてもよい。
この後、減衰率マップデータ生成部363は、視覚情報記憶部324を参照して、各部分領域の減衰率に応じた視覚情報を各部分領域に対して割り当てることによって減衰率マップデータを生成し、合成画像データ生成部362へ出力する(ステップS39)。なお、合成画像データ生成部362が、各部分領域の減衰率の推定値を文字情報としてさらに表示する合成画像データを生成するようにしてもよい。
合成画像データ生成部362は、Bモード画像データに減衰率マップデータを重畳することによって合成画像データを生成し、この合成画像データを表示装置4へ出力する(ステップS40)。合成画像データを受信した表示装置4は、その合成画像データに対応する合成画像を表示する。図21は、表示装置4が表示する減衰率マップデータ付きの合成画像の表示例を示す図である。図21に示す合成画像201は、各領域ごとに異なる視覚情報が割り当てられている。なお、図21では、視覚情報を模様で模式的に記載している。また、図21では、簡単のためBモード画像の具体的な表示を省略している。
本実施の形態3において、深度幅wは例えば1cm程度であるのが好ましい。深度幅wが1cm程度である場合、周波数解析部352がFFT処理を行う際のRFデータを切り取る区間の幅ΔLを2mm程度とするのが望ましい。このとき、間隔ΔLに対応する基準位置の深度z’の幅Δz’(=ΔL/2)は、1mm程度になる。
以上説明した本発明の実施の形態3によれば、周波数スペクトルまたは周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率を所定の順序で算出することによって得られる2次変化率を用いて、超音波振動子の走査領域における超音波信号の単位距離および単位周波数あたりの減衰率を推定するため、実施の形態1、2と同様、超音波の減衰率を精度よくかつ簡便に算出することができるとともに、減衰率に基づく画像の信頼性を向上させることができる。
また、本実施の形態3によれば、超音波振動子の走査領域を分割することによって得られる複数の部分領域の各々に対し、各部分領域で得られた減衰率の統計量を算出することによって各部分領域の減衰率を推定した後、各部分領域の減衰率の値に応じた視覚情報を付与することによって減衰率マップデータを生成するため、ユーザが減衰率の分布を視覚的に把握しやすい情報を提供することができる。
また、本実施の形態3において、減衰率マップの最小単位である扇形領域のうち隣接する扇形を重複させるように設定してもよい。
また、本実施の形態3において、各扇形領域の減衰率の算出を順次行う代わりに、同時並行的に行ってもよい。
また、本実施の形態3において、入力部37が入力を受け付ける設定信号に基づいて1つの扇形領域を関心領域(ROI:Region of Interest)として設定し、関心領域内の減衰率の値をさらに合成することによって合成画像データを生成してもよい。
(その他の実施の形態)
ここまで、本発明を実施するための形態を説明してきたが、本発明は、上述した実施の形態1〜3によってのみ限定されるべきものではない。例えば、超音波診断装置において、各機能を有する回路同士をバスで接続することによって構成してもよいし、一部の機能が他の機能の回路構造に内蔵されるように構成してもよい。具体的には、例えば減衰率推定部の機能を有する回路に変化率算出部の機能を内蔵してもよい。
また、超音波プローブとして、光学系のない細径の超音波ミニチュアプローブを適用してもよい。超音波ミニチュアプローブは、通常、胆道、胆管、膵管、気管、気管支、尿道、尿管へ挿入され、その周囲臓器(膵臓、肺、前立腺、膀胱、リンパ節等)を観察する際に用いられる。
また、超音波プローブとして、被検体の体表から超音波を照射する体外式超音波プローブを適用してもよい。体外式超音波プローブは、通常、腹部臓器(肝臓、胆嚢、膀胱)、乳房(特に乳腺)、甲状腺を観察する際に用いられる。
また、超音波振動子は、リニア振動子でもラジアル振動子でもコンベックス振動子でも構わない。超音波振動子がリニア振動子である場合、その走査領域は矩形(長方形、正方形)をなし、超音波振動子がラジアル振動子やコンベックス振動子である場合、その走査領域は扇形や円環状をなす。
このように、本発明は、特許請求の範囲に記載した技術的思想を逸脱しない範囲内において、様々な実施の形態等を含み得るものである。