以下、添付図面を参照して、本発明を実施するための形態(以下、「実施の形態」という)を説明する。
図1は、本発明の一実施の形態に係る超音波観測装置を備えた超音波診断システムの機能構成を示すブロック図である。同図に示す超音波診断システム1は、観測対象である被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波内視鏡2と、超音波内視鏡2が取得した超音波信号に基づいて超音波画像を生成する超音波観測装置3と、超音波観測装置3が生成した超音波画像を表示する表示装置4と、を備える。
超音波内視鏡2は、その先端部に、超音波観測装置3から受信した電気的なパルス信号を超音波パルス(音響パルス)に変換して被検体へ照射するとともに、被検体で反射された超音波エコーを電圧変化で表現する電気的なエコー信号に変換して出力する超音波振動子21を有する。超音波振動子21は、コンベックス振動子、リニア振動子およびラジアル振動子のいずれでも構わない。超音波内視鏡2は、超音波振動子21をメカ的に走査させるものであってもよいし、超音波振動子21として複数の素子をアレイ状に設け、送受信にかかわる素子を電子的に切り替えたり、各素子の送受信に遅延をかけたりすることで、電子的に走査させるものであってもよい。以下、本実施の形態では、説明の便宜上、超音波振動子21がリニア振動子であるものとする。
超音波内視鏡2は、通常は撮像光学系および撮像素子を有しており、被検体の消化管(食道、胃、十二指腸、大腸)または呼吸器(気管、気管支)へ挿入され、消化管、呼吸器やその周囲臓器(膵臓、胆嚢、胆管、胆道、リンパ節、縦隔臓器、血管等)を撮像することが可能である。また、超音波内視鏡2は、撮像時に被検体へ照射する照明光を導くライトガイドを有する。このライトガイドは、先端部が超音波内視鏡2の被検体への挿入部の先端まで達している一方、基端部が照明光を発生する光源装置に接続されている。
超音波観測装置3は、超音波内視鏡2と電気的に接続され、所定の波形および送信タイミングに基づいて高電圧パルスからなる送信信号(パルス信号)を超音波振動子21へ送信するとともに、超音波振動子21から電気的な受信信号であるエコー信号を受信してデジタルの高周波(RF:Radio Frequency)信号のデータ(以下、RFデータという)を生成、出力する送受信部31と、送受信部31から受信したRFデータをもとにデジタルのBモード用受信データを生成する信号処理部32と、送受信部31から受信したRFデータに対して所定の演算を施す演算部33と、各種画像データを生成する画像処理部34と、キーボード、マウス、タッチパネル等のユーザインタフェースを用いて実現され、各種情報の入力を受け付ける入力部35と、超音波診断システム1全体を制御する制御部36と、超音波観測装置3の動作に必要な各種情報を記憶する記憶部37と、を備える。
送受信部31は、エコー信号を増幅する信号増幅部311を有する。信号増幅部311は、受信深度が大きいエコー信号ほど高い増幅率で増幅するSTC(Sensitivity Time Control)補正を行う。図2は、信号増幅部311が行う増幅処理における受信深度と増幅率との関係を示す図である。図2に示す受信深度zは、超音波の受信開始時点からの経過時間に基づいて算出される量である。図2に示すように、増幅率β(dB)は、受信深度zが閾値zthより小さい場合、受信深度zの増加に伴ってβ0からβth(>β0)へ線型に増加する。また、増幅率β(dB)は、受信深度zが閾値zth以上である場合、一定値βthをとる。閾値zthの値は、観測対象から受信する超音波信号がほとんど減衰してしまい、ノイズが支配的になるような値である。より一般に、増幅率βは、受信深度zが閾値zthより小さい場合、受信深度zの増加に伴って単調増加すればよい。図2に示す関係は、予め記憶部37に記憶されている。
送受信部31は、信号増幅部311によって増幅されたエコー信号に対してフィルタリング等の処理を施した後、A/D変換することによって時間ドメインのRFデータを生成し、信号処理部32および演算部33へ出力する。超音波内視鏡2が複数の素子をアレイ状に設けた超音波振動子21を電子的に走査させる構成を有する場合、送受信部31は、複数の素子に対応したビーム合成用の多チャンネル回路を有する。
送受信部31が送信するパルス信号の周波数帯域は、超音波振動子21におけるパルス信号の超音波パルスへの電気音響変換の線型応答周波数帯域をほぼカバーする広帯域にするとよい。また、信号増幅部311におけるエコー信号の各種処理周波数帯域は、超音波振動子21による超音波エコーのエコー信号への音響電気変換の線型応答周波数帯域をほぼカバーする広帯域にするとよい。これらにより、後述する周波数スペクトルの近似処理を実行する際、精度のよい近似を行うことが可能となる。
送受信部31は、制御部36が出力する各種制御信号を超音波内視鏡2に対して送信するとともに、超音波内視鏡2から識別用のIDを含む各種情報を受信して制御部36へ送信する機能も有する。
信号処理部32は、RFデータに対してバンドパスフィルタ、包絡線検波、対数変換など公知の処理を施し、デジタルのBモード用受信データを生成する。対数変換では、RFデータを基準電圧で除した量の常用対数をとってデシベル値で表現する。信号処理部32は、生成したBモード用受信データを、画像処理部34へ出力する。信号処理部32は、CPU(Central Processing Unit)や各種演算回路等を用いて実現される。
演算部33は、送受信部31が出力したRFデータに対して受信深度によらず増幅率を一定とするよう増幅補正を行う増幅補正部331と、増幅補正を行ったRFデータに高速フーリエ変換(FFT:Fast Fourier Transform)を施して周波数解析を行うことにより超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析部332と、各周波数スペクトルの特徴量を算出する特徴量算出部333と、超音波画像を複数の領域に分割したときの各分割領域において、超音波が観測対象を伝播する際の互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値の各々を用いることにより、各周波数スペクトルの特徴量に対して超音波の影響を排除する減衰補正を施すことによって減衰率候補値ごとの各周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに複数の減衰率候補値の中から観測対象に最適な減衰率を設定する減衰率設定部334と、減衰率設定部334が分割領域ごとに設定した最適な減衰率のうち超音波振動子の表面とサンプリング点との間に存在する分割領域の最適な減衰率を用いてサンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて特徴量の減衰補正を行うことによって補正特徴量を算出する特徴量補正部335と、を有する。演算部33は、CPUや各種演算回路等を用いて実現される。
図3は、増幅補正部331が行う増幅補正処理における受信深度と増幅率との関係を示す図である。図3に示すように、増幅補正部331が行う増幅処理における増幅率β(dB)は、受信深度zがゼロのとき最大値βth−β0をとり、受信深度zがゼロから閾値zthに達するまで線型に減少し、受信深度zが閾値zth以上のときゼロである。なお、図3に示す関係は、予め記憶部37に記憶されている。増幅補正部331が図3に示す関係に基づいてデジタルRF信号を増幅補正することにより、信号増幅部311におけるSTC補正の影響を相殺し、一定の増幅率βthの信号を出力することができる。なお、増幅補正部331が行う増幅補正処理における受信深度zと増幅率βの関係は、信号増幅部311が行う増幅処理における受信深度と増幅率の関係に応じて異なることは勿論である。
このような増幅補正を行う理由を説明する。STC補正は、アナログ信号波形の振幅を全周波数帯域にわたって均一に、かつ、深度に対しては単調増加する増幅率で増幅させることで、アナログ信号波形の振幅から減衰の影響を排除する補正処理である。このため、エコー信号の振幅を輝度に変換して表示するBモード画像を生成する場合、かつ、一様な組織を走査した場合には、STC補正を行うことによって深度によらず輝度値が一定になる。すなわち、Bモード画像の輝度値から減衰の影響を排除する効果を得ることができる。
一方、本実施の形態のように超音波の周波数スペクトルを算出して解析した結果を利用する場合、STC補正でも超音波の伝播に伴う減衰の影響を正確に排除できるとは限らない、という問題がある。なぜなら、一般に減衰量は周波数によって異なるが(後述する式(1)を参照)、STC補正の増幅率は距離だけに応じて変化し、周波数依存性がないためである。
上述した問題を解決するには、Bモード画像を生成する際にSTC補正を施した受信信号を出力する一方、周波数スペクトルに基づいた画像を生成する際に、Bモード画像を生成するための送信とは異なる新たな送信を行い、STC補正を施していない受信信号を出力することが考えられる。ところがこの場合には、受信信号に基づいて生成される画像データのフレームレートが低下してしまうという問題がある。
そこで、本実施の形態では、生成される画像データのフレームレートを維持しつつ、Bモード画像用にSTC補正を施した信号に対してSTC補正の影響を排除するために、増幅補正部331によって増幅率の補正を行う。
周波数解析部332は、増幅補正部331が増幅補正した各音線のRFデータ(ラインデータ)を所定の時間間隔でサンプリングし、サンプルデータを生成する。周波数解析部332は、サンプルデータ群にFFT処理を施すことにより、RFデータ上の複数の箇所(データ位置)における周波数スペクトルを算出する。
図4は、超音波信号の1つの音線におけるデータ配列を模式的に示す図である。同図に示す音線SRkにおいて、白または黒の長方形は、1つのサンプル点におけるデータを意味している。また、音線SRkにおいて、右側に位置するデータほど、超音波振動子21から音線SRkに沿って計った場合の深い箇所からのサンプルデータである(図4の矢印を参照)。音線SRkは、送受信部31が行うA/D変換におけるサンプリング周波数(例えば50MHz)に対応した時間間隔で離散化されている。図4では、番号kの音線SRkの8番目のデータ位置を受信深度zの方向の初期値Z(k) 0として設定した場合を示しているが、初期値の位置は任意に設定することができる。周波数解析部332による算出結果は複素数で得られ、記憶部37に格納される。
図4に示すデータ群Fj(j=1、2、・・・、K)は、FFT処理の対象となるサンプルデータ群である。一般に、FFT処理を行うためには、サンプルデータ群が2のべき乗のデータ数を有している必要がある。この意味で、サンプルデータ群Fj(j=1、2、・・・、K−1)はデータ数が16(=24)で正常なデータ群である一方、サンプルデータ群FKは、データ数が12であるため異常なデータ群である。異常なデータ群に対してFFT処理を行う際には、不足分だけゼロデータを挿入することにより、正常なサンプルデータ群を生成する処理を行う。この点については、周波数解析部332の処理を説明する際に詳述する(図10を参照)。
図5は、周波数解析部332が算出した周波数スペクトルの例を示す図である。ここでいう「周波数スペクトル」とは、サンプルデータ群にFFT処理を施すことによって得られた「ある受信深度zにおける強度の周波数分布」を意味する。また、ここでいう「強度」とは、例えばエコー信号の電圧、エコー信号の電力、超音波エコーの音圧、超音波エコーの音響エネルギー等のパラメータ、これらパラメータの振幅や時間積分値やその組み合わせのいずれかを指す。
図5では、横軸に周波数fを取っている。また、図5では、縦軸に、強度I0を基準強度Ic(定数)で除した量の常用対数(デシベル表現)I=10log10(I0/Ic)を取っている。図5において、受信深度zは一定である。図5に示す直線L10については後述する。なお、本実施の形態において、曲線および直線は、離散的な点の集合からなる。
図5に示す周波数スペクトルC1において、以後の演算に使用する周波数帯域の下限周波数fLおよび上限周波数fHは、超音波振動子21の周波数帯域、送受信部31が送信するパルス信号の周波数帯域などをもとに決定されるパラメータであり、例えばfL=3MHz、fH=10MHzである。以下、図5において、下限周波数fLおよび上限周波数fHによって定まる周波数帯域を「周波数帯域U」という。
一般に、周波数スペクトルは、観測対象が生体組織である場合、超音波が走査された生体組織の性状によって異なる傾向を示す。これは、周波数スペクトルが、超音波を散乱する散乱体の大きさ、数密度、音響インピーダンス等と相関を有しているためである。ここでいう「生体組織の性状」とは、例えば悪性腫瘍(癌)、良性腫瘍、内分泌腫瘍、粘液性腫瘍、正常組織、嚢胞、脈管などのことである。
特徴量算出部333は、所定周波数帯域における周波数スペクトルに対して回帰分析を行って一次式で近似することにより、この近似した一次式を特徴付ける特徴量を算出する。例えば、図5に示す周波数スペクトルC1の場合、特徴量算出部333は、周波数帯域Uで回帰分析を行うことによって近似直線L10を得る。以下、近似直線L10を周波数fの一次式I=a0f+b0で表すと、特徴量算出部333は、直線L10に対応する特徴量として、傾きa0、切片b0、および周波数帯域Uの中心周波数fM=(fs+fe)/2における強度Iの値であるミッドバンドフィット(Mid-band fit)c0=a0fM+b0を算出する。なお、特徴量算出部333は、回帰分析によって二次以上の多項式で周波数スペクトルを近似するようにしてもよい。
3つの補正前特徴量のうち、傾きa0は、超音波の散乱体の大きさと相関を有し、一般に散乱体が大きいほど傾きが小さな値を有すると考えられる。また、切片b0は、散乱体の大きさ、音響インピーダンスの差、散乱体の数密度(濃度)等と相関を有している。具体的には、切片b0は、散乱体が大きいほど大きな値を有し、音響インピーダンスの差が大きいほど大きな値を有し、散乱体の数密度が大きいほど大きな値を有すると考えられる。ミッドバンドフィットc0は、傾きa0と切片b0から導出される間接的なパラメータであり、有効な周波数帯域内の中心におけるスペクトルの強度を与える。このため、ミッドバンドフィットc0は、散乱体の大きさ、音響インピーダンスの差、散乱体の数密度に加えて、Bモード画像の輝度とある程度の相関を有していると考えられる。
減衰率設定部334は、超音波画像を複数の領域に分割したときの各分割領域において、単位長さおよび単位周波数あたりの超音波の減衰量を与える減衰率を用いて特徴量の減衰補正を行う。図6は、超音波画像における分割領域の設定例を模式的に示す図である。図6では、矩形状をなす超音波画像101を30個の分割領域Pij(i=1〜5、j=1〜6)に分割した場合を示している。なお、図6の破線は音線を表している。図6に示す場合、分割領域Pijは、すべて同じ面積を有している。このため、分割領域Pijにおける超音波振動子21の表面位置102からの深度方向の長さH(以下、分割領域Pijの高さHという)は一定である。分割領域Pijに関する情報は、記憶部37が有する分割領域情報記憶部371に格納されている。なお、音線の数はあくまでも一例に過ぎず、図6に示す場合に限定されるものではない。また、図6は超音波振動子21がリニア振動子である場合の超音波画像であるため矩形状をなしているが、他のタイプの超音波振動子21の場合に超音波画像の外形状が異なることはいうまでもない。
一般に、超音波の減衰量A(f,z)は、超音波が受信深度0と受信深度zとの間を往復する間に生じる減衰であり、往復する前後の強度変化(デシベル表現での差)として定義される。減衰量A(f,z)は、一様な組織内では周波数に比例することが経験的に知られており、以下の式(1)で表現される。
A(f,z)=2αzf ・・・(1)
ここで、比例定数αは減衰率と呼ばれる量であり、単位長さおよび単位周波数あたりの超音波の減衰量を与える。また、zは超音波の受信深度であり、fは周波数である。減衰率αの具体的な値は、観測対象が生体である場合、生体の部位に応じて定まる。減衰率αの単位は、例えばdB/cm/MHzである。
減衰率設定部334は、複数の減衰率候補値の中から最適な減衰率を設定する。この際、減衰率設定部334は、減衰率候補値αを用いて、特徴量算出部333が算出した特徴量(傾きa0、切片b0、ミッドバンドフィットc0)に対し、以下に示す式(2)〜(4)にしたがって減衰補正を行うことにより、予備補正特徴量a、b、cを算出する。
a=a0+2αz ・・・(2)
b=b0 ・・・(3)
c=c0+A(fM,z)=c0+2αzfM(=afM+b) ・・・(4)
式(2)、(4)からも明らかなように、減衰率設定部334は、超音波の受信深度zが大きいほど、補正量が大きい補正を行う。また、式(3)によれば、切片に関する補正は恒等変換である。これは、切片が周波数0(Hz)に対応する周波数成分であって減衰の影響を受けないためである。
図7は、減衰率設定部334が補正した予備補正特徴量a、b、cをパラメータとして有する直線を示す図である。直線L1の式は、
I=af+b=(a0+2αz)f+b0 ・・・(5)
で表される。この式(5)からも明らかなように、直線L1は、減衰補正前の直線L10と比較して、傾きが大きく(a>a0)、かつ切片が同じ(b=b0)である。
減衰率設定部334は、各分割領域において、減衰率候補値ごとに算出した予備補正特徴量の統計的なばらつきが最小である減衰率候補値を最適な減衰率として設定する。本実施の形態では、統計的なばらつきを示す量として分散を適用する。この場合、減衰率設定部334は、分散が最小となる減衰率候補値を最適な減衰率として設定する。上述した3つの予備補正特徴量a、b、cのうち独立なのは2つである。加えて、予備補正特徴量bは減衰率に依存しない。したがって、予備補正特徴量a、cに対して最適な減衰率を設定する場合、減衰率設定部334は、予備補正特徴量aおよびcのいずれか一方の分散を算出すればよい。
ただし、減衰率設定部334が、予備補正特徴量aを用いて最適な減衰率を設定する場合は予備補正特徴量aの分散を適用し、予備補正特徴量cを用いて最適な減衰率を設定する場合は予備補正特徴量cの分散を適用するのがより好ましい。これは、減衰量A(f,z)を与える式(1)があくまで理想的なものに過ぎず、現実には以下の式(6)の方が適切であることによる。
A(f,z)=2αzf+2α1z ・・・(6)
式(6)の右辺第2項のα1は、超音波の受信深度zに比例して信号強度が変化する大きさを表す係数であり、観測対象の組織が不均一であることや、ビーム合成時のチャンネル数の変更などに起因して発生する信号強度の変化を表す係数である。式(6)の右辺第2項が存在するため、予備補正特徴量cを用いて最適な減衰率を設定する場合は、予備補正特徴量cの分散を適用した方が正確に減衰を補正することができる(式(4)を参照)。一方、周波数fに比例する係数である予備補正特徴量aを用いて最適な減衰率を設定する場合は、予備補正特徴量aの分散を適用した方が、式(6)の右辺第2項の影響を排除して正確に減衰を補正することができる。
ここで、統計的なばらつきに基づいて最適な減衰率を設定することができる理由を説明する。観測対象に最適な減衰率を適用した場合、観測対象と超音波振動子21との距離に関わらず、特徴量は観測対象に固有の値へ収束し、統計的なばらつきが小さくなると考えられる。その一方で、観測対象に適合しない減衰率候補値を最適な減衰率とした場合、減衰補正が過剰であるかまたは不足するため、超音波振動子21との距離に応じて特徴量にずれが生じ、特徴量の統計的なばらつきが大きくなると考えられる。したがって、統計的なばらつきが最も小さい減衰率候補値が、観測対象にとって最適な減衰率であるということができる。
図8は、同じ観測対象に対して2つの異なる減衰率候補値に基づいてそれぞれ減衰補正された予備補正特徴量の分布例を模式的に示す図である。図8では、横軸を予備補正特徴量とし、縦軸を頻度としている。図8に示す2つの分布曲線N1、N2は、頻度の総和が同じである。図8に示す場合、分布曲線N1は、分布曲線N2と比較して特徴量の統計的なばらつきが小さく(分散が小さく)、山が急峻な形状をなす。したがって、減衰率設定部334は、この2つの分布曲線N1、N2に対応する2つの減衰率候補値から最適な減衰率を設定する場合、分布曲線N1に対応する減衰率候補値を最適な減衰率として設定する。
特徴量補正部335は、分割領域ごとの最適な減衰率を用いてサンプリング点における単位周波数あたりの累積減衰率(以下、単に累積減衰率ともいう)を算出し、該累積減衰率を用いて特徴量の減衰補正を行う。任意のサンプリング点における累積減衰率は、超音波振動子21の表面からの距離とその間に存在する分割領域における最適な減衰率を用いて算出される。分割領域PIJ(1≦I≦imax、1≦J≦jmax;I,Jは定数)で超音波振動子21に近い側の境界からの距離がh(0<h≦H)であるようなサンプリング点SIJ(h)における累積減衰率γIJ(h)は、
γIJ(h)=[Σj=1,2,…,J-1(2H・α(PIj))]+2h・α(PIJ)
・・・(7)
と表される。ここで、Σj=1,2,…,J-1はj=1〜J−1の和を取ることを意味し、α(PIj)は分割領域PIjにおける最適な減衰率を表している。式(7)の右辺の第1項は、分割領域ごとの超音波の往復距離2Hに対して最適な減衰率α(PIj)を乗じたものの和である。また、式(7)の右辺の第2項は、分割領域PIJにおけるサンプリング点SIJ(h)までの超音波の往復距離2hに分割領域PIJにおける最適な減衰率α(PIJ)を乗じたものである。このようにして、特徴量補正部335は、超音波振動子21の表面からの減衰率を累積していくことによって累積減衰率γIJ(h)を算出する。最適な減衰率の単位をdB/cm/MHzとする場合、累積減衰率の単位はdB/MHzである。
例えば図6では、式(7)でI=2、J=3(ただしimax=5、jmax=6)の場合を例示している。この場合、サンプリング点S23(h)における累積減衰率γ23(h)は、式(7)より、
γ23(h)=2H・α(P21)+2H・α(P22)+2h・α(P23)
・・・(8)
となる。
特徴量補正部335は、累積減衰率γIJ(h)を用いてサンプリング点SIJ(h)における特徴量を以下のように減衰補正する。
aIJ(h)=a0+2γIJ(h) ・・・(9)
bIJ(h)=b0 ・・・(10)
cIJ(h)=c0+2fMγIJ(h) ・・・(11)
画像処理部34は、エコー信号の振幅を輝度に変換して表示する超音波画像であるBモード画像データを生成するBモード画像データ生成部341と、特徴量算出部333が算出した特徴量に関する情報を表示する特徴量画像データを生成する特徴量画像データ生成部342と、を有する。
Bモード画像データ生成部341は、信号処理部32から受信したBモード用受信データに対してゲイン処理、コントラスト処理等の公知の技術を用いた信号処理を行うとともに、表示装置4における画像の表示レンジに応じて定まるデータステップ幅に応じたデータの間引き等を行うことによってBモード画像データを生成する。Bモード画像は、色空間としてRGB表色系を採用した場合の変数であるR(赤)、G(緑)、B(青)の値を一致させたグレースケール画像である。
Bモード画像データ生成部341は、Bモード用受信データに走査範囲を空間的に正しく表現できるよう並べ直す座標変換を施した後、Bモード用受信データ間の補間処理を施すことによってBモード用受信データ間の空隙を埋め、Bモード画像データを生成する。Bモード画像データ生成部341は、生成したBモード画像データを特徴量画像データ生成部342へ出力する。
特徴量画像データ生成部342は、特徴量算出部333が算出した特徴量に関連する視覚情報をBモード画像データにおける画像の各画素に対して重畳することによって特徴量画像データを生成する。特徴量画像データ生成部342は、例えば図4に示す1つの振幅データ群Fj(j=1、2、・・・、K)のデータ量に対応する画素領域に対し、その振幅データ群Fjから算出される周波数スペクトルの特徴量に対応する視覚情報を割り当てる。特徴量画像データ生成部342は、例えば上述した傾き、切片、ミッドバンドフィットのいずれか一つに視覚情報としての色相を対応付けることによって特徴量画像データを生成する。なお、特徴量画像データ生成部342が、傾き、切片、ミッドバンドフィットから選択される2つの特徴量の一方に色相を対応付けるとともに、他方に明暗を対応付けることによって特徴量画像データを生成するようにしてもよい。特徴量に関連する視覚情報としては、色相や明暗(明度)のほか、例えば彩度、輝度値、R(赤)、G(緑)、B(青)などの所定の表色系を構成する色空間の変数を挙げることができる。
制御部36は、演算および制御機能を有するCPU(Central Processing Unit)や各種演算回路等を用いて実現される。制御部36は、記憶部37が記憶、格納する情報を記憶部37から読み出し、超音波観測装置3の作動方法に関連した各種演算処理を実行することによって超音波観測装置3を統括して制御する。なお、制御部36を信号処理部32および演算部33と共通のCPU等を用いて構成することも可能である。
記憶部37は、分割領域に関する情報を記憶する分割領域情報記憶部371と、周波数解析部332が算出した周波数スペクトルの情報を受信深度および受信方向とともに記憶するスペクトル情報記憶部372と、特徴量算出部333が算出する特徴量および特徴量補正部335が補正した補正特徴量に関する情報を記憶する特徴量情報記憶部373と、減衰率設定部334が分割領域ごとに設定する最適な減衰率および特徴量補正部335が算出するサンプリング点毎の累積減衰率に関する情報を記憶する減衰率情報記憶部374と、を有する。
記憶部37は、上記以外にも、例えば増幅処理に必要な情報(図2に示す増幅率と受信深度との関係)、増幅補正処理に必要な情報(図3に示す増幅率と受信深度との関係)、減衰補正処理に必要な情報(式(1)参照)、周波数解析処理に必要な窓関数(Hamming、Hanning、Blackman等)の情報等を記憶する。
また、記憶部37は、超音波観測装置3の作動方法を実行するための作動プログラムを含む各種プログラムを記憶する。作動プログラムは、ハードディスク、フラッシュメモリ、CD−ROM、DVD−ROM、フレキシブルディスク等のコンピュータ読み取り可能な記録媒体に記録して広く流通させることも可能である。なお、上述した各種プログラムは、通信ネットワークを介してダウンロードすることによって取得することも可能である。ここでいう通信ネットワークは、例えば既存の公衆回線網、LAN(Local Area Network)、WAN(Wide Area Network)などによって実現されるものであり、有線、無線を問わない。
以上の構成を有する記憶部37は、各種プログラム等が予めインストールされたROM(Read Only Memory)、および各処理の演算パラメータやデータ等を記憶するRAM(Random Access Memory)等を用いて実現される。
図9は、以上の構成を有する超音波観測装置3が行う処理の概要を示すフローチャートである。具体的には、超音波観測装置3が超音波内視鏡2からエコー信号を受信する以降の処理の概要を示すフローチャートである。以下、図9を参照して、超音波観測装置3が行う処理を説明する。まず、超音波観測装置3は、超音波内視鏡2から超音波振動子21による観測対象の測定結果としてのエコー信号を受信する(ステップS1)。
超音波振動子21からエコー信号を受信した信号増幅部311は、そのエコー信号の増幅を行う(ステップS2)。ここで、信号増幅部311は、例えば図2に示す増幅率と受信深度との関係に基づいてエコー信号の増幅(STC補正)を行う。
続いて、Bモード画像データ生成部341は、信号増幅部311が増幅したエコー信号を用いてBモード画像データを生成して、表示装置4へ出力する(ステップS3)。Bモード画像データを受信した表示装置4は、そのBモード画像データに対応するBモード画像を表示する。
増幅補正部331は、送受信部31から出力されたRFデータに対して受信深度によらず増幅率が一定となるように増幅補正を行う(ステップS4)。ここで、増幅補正部331は、例えば図3に示す増幅率と受信深度との関係に基づいて増幅補正を行う。
この後、周波数解析部332は、増幅補正後の各音線のRFデータに対してFFTによる周波数解析を行うことによって全てのサンプルデータ群に対する周波数スペクトルを算出し、スペクトル情報記憶部372へ格納する(ステップS5)。図10は、ステップS5において周波数解析部332が行う処理の概要を示すフローチャートである。以下、図10に示すフローチャートを参照して、周波数解析処理を詳細に説明する。
まず、周波数解析部332は、解析対象の音線を識別するカウンタkをk0とする(ステップS21)。
続いて、周波数解析部332は、FFT演算用に生成した一連のデータ群(サンプルデータ群)を代表するデータ位置(受信深度に相当)Z(k)の初期値Z(k) 0を設定する(ステップS22)。例えば、図4では、上述したように、音線SRkの8番目のデータ位置を初期値Z(k) 0として設定した場合を示している。
その後、周波数解析部332は、サンプルデータ群を取得し(ステップS23)、取得したサンプルデータ群に対し、記憶部37が記憶する窓関数を作用させる(ステップS24)。このようにサンプルデータ群に対して窓関数を作用させることにより、サンプルデータ群が境界で不連続になることを回避し、アーチファクトが発生するのを防止することができる。
続いて、周波数解析部332は、データ位置Z(k)のサンプルデータ群が正常なデータ群であるか否かを判定する(ステップS25)。図4を参照した際に説明したように、サンプルデータ群は、2のべき乗のデータ数を有している必要がある。以下、正常なサンプルデータ群のデータ数を2n(nは正の整数)とする。本実施の形態では、データ位置Z(k)が、できるだけZ(k)が属するサンプルデータ群の中心になるよう設定される。具体的には、サンプルデータ群のデータ数は2nであるので、Z(k)はそのサンプルデータ群の中心に近い2n/2(=2n-1)番目の位置に設定される。この場合、サンプルデータ群が正常であるとは、データ位置Z(k)より浅部側に2n-1−1(=Nとする)個のデータがあり、データ位置Z(k)より深部側に2n-1(=Mとする)個のデータがあることを意味する。図4に示す場合、サンプルデータ群Fj(j=1、2、・・・、K−1)は正常である。なお、図4ではn=4(N=7,M=8)の場合を例示している。
ステップS25における判定の結果、データ位置Z(k)のサンプルデータ群が正常である場合(ステップS25:Yes)、周波数解析部332は、後述するステップS27へ移行する。
ステップS25における判定の結果、データ位置Z(k)のサンプルデータ群が正常でない場合(ステップS25:No)、周波数解析部332は、不足分だけゼロデータを挿入することによって正常なサンプルデータ群を生成する(ステップS26)。ステップS25において正常でないと判定されたサンプルデータ群(例えば図4のサンプルデータ群FK)は、ゼロデータを追加する前に窓関数が作用されている。このため、サンプルデータ群にゼロデータを挿入してもデータの不連続は生じない。ステップS26の後、周波数解析部332は、後述するステップS27へ移行する。
ステップS27において、周波数解析部332は、サンプルデータ群を用いてFFT演算を行うことにより、振幅の周波数分布である周波数スペクトルを得る(ステップS27)。
続いて、周波数解析部332は、データ位置Z(k)をステップ幅Dで変化させる(ステップS28)。ステップ幅Dは、記憶部37が予め記憶しているものとする。図4では、D=15の場合を例示している。ステップ幅Dは、Bモード画像データ生成部341がBモード画像データを生成する際に利用するデータステップ幅と一致させることが望ましいが、周波数解析部332における演算量を削減したい場合には、ステップ幅Dとしてデータステップ幅より大きい値を設定してもよい。
その後、周波数解析部332は、データ位置Z(k)が音線SRkにおける最大値Z(k) maxより大きいか否かを判定する(ステップS29)。データ位置Z(k)が最大値Z(k) maxより大きい場合(ステップS29:Yes)、周波数解析部332はカウンタkを1増加させる(ステップS30)。これは、処理をとなりの音線へ移すことを意味する。一方、データ位置Z(k)が最大値Z(k) max以下である場合(ステップS29:No)、周波数解析部332はステップS23へ戻る。
ステップS30の後、周波数解析部332は、カウンタkが最大値kmaxより大きいか否かを判定する(ステップS31)。カウンタkがkmaxより大きい場合(ステップS31:Yes)、周波数解析部332は一連の周波数解析処理を終了する。一方、カウンタkがkmax以下である場合(ステップS31:No)、周波数解析部332はステップS22に戻る。この最大値kmaxは、術者等のユーザが入力部35を通じて任意に指示入力した値、もしくは、記憶部37にあらかじめ設定された値とする。
このようにして、周波数解析部332は、解析対象領域内の(kmax−k0+1)本の音線の各々について複数回のFFT演算を行う。FFT演算の結果として得られる周波数スペクトルは、受信深度、受信方向とともにスペクトル情報記憶部372に格納される。
なお、以上の説明では、周波数解析部332が超音波信号を受信したすべての領域に対して周波数解析処理を行うものとしたが、入力部35が特定の深度幅および音線幅で区切られる部分領域の設定入力を受け付け可能な構成とし、設定された部分領域内においてのみ周波数解析処理を行うようにすることも可能である。
以上説明したステップS5の周波数解析処理に続いて、特徴量算出部333は、各分割領域に含まれるサンプリング点で周波数スペクトルの特徴量を算出する(ステップS6)。具体的には、特徴量算出部333は、所定の周波数帯域の周波数スペクトルに対して回帰分析を行うことによって一次式I=a0f+b0で近似し、特徴量として傾きa0、切片b0、ミッドバンドフィットc0を算出する。例えば、図5に示す直線L10は、特徴量算出部333が周波数帯域Uの周波数スペクトルC1に対し回帰分析によって近似した回帰直線である。
この後、減衰率設定部334は、後述する減衰補正を行う際に適用する減衰率候補値αの値を所定の初期値α0に設定する(ステップS7)。この初期値α0の値は、予め減衰率情報記憶部374が記憶しておくようにすればよい。
続いて、減衰率設定部334は、特徴量算出部333が各周波数スペクトルに対して近似した特徴量に対し、減衰率候補値をαとして減衰補正を行うことにより、予備補正特徴量を算出し、減衰率候補値αとともに特徴量情報記憶部373に格納する(ステップS8)。図7に示す直線L1は、減衰率設定部334が減衰補正処理を行うことによって得られる直線の例である。
ステップS8において、減衰率設定部334は、上述した式(2)、(4)における受信深度zに、超音波信号の音線のデータ配列を用いて得られるデータ位置Z=(fsp/2vs)Dnを代入することによって算出する。ここで、fspはデータのサンプリング周波数、vsは音速、Dはデータステップ幅、nは処理対象の振幅データ群のデータ位置までの音線の1番目のデータからのデータステップ数である。例えば、データのサンプリング周波数fspを50MHzとし、音速vsを1530m/secとし、図4に示すデータ配列を採用してステップ幅Dを15とすると、z=0.2295n(mm)となる。
減衰率設定部334は、減衰率設定部334が各周波数スペクトルに対して減衰補正することによって得られた複数の予備補正特徴量のうち代表となる予備補正特徴量の分散を算出し、減衰率候補値αと対応づけて特徴量情報記憶部373へ格納する(ステップS9)。例えば、予備補正特徴量が傾きa、ミッドバンドフィットcである場合、上述したように、減衰率設定部334は、予備補正特徴量aおよびcのいずれか一方の分散を算出する。このステップS9において、予備補正特徴量aを用いて特徴量画像を生成する場合は各分割領域における予備補正特徴量aの分散を適用し、予備補正特徴量cを用いて特徴量画像を生成する場合は各分割領域における予備補正特徴量cの分散を適用するのが好ましい。
この後、減衰率設定部334は、減衰率候補値αの値をΔαだけ増加させ(ステップS10)、増加後の減衰率候補値αと所定の最大値αmaxとの大小を比較する(ステップS11)。ステップS11における比較の結果、減衰率候補値αが最大値αmaxより大きい場合(ステップS11:Yes)、超音波観測装置3はステップS12へ移行する。一方、ステップS11における比較の結果、減衰率候補値αが最大値αmax以下である場合(ステップS11:No)、超音波観測装置3はステップS8へ戻る。
ステップS12において、減衰率設定部334は、分割領域ごとに、特徴量情報記憶部373が記憶する減衰率候補値ごとの予備補正特徴量の分散を参照し、該分散が最小である減衰率候補値を分割領域における最適な減衰率として設定する(ステップS12)。
図11は、減衰率設定部334が行う処理の概要を示す図である。α0=0(dB/cm/MHz)、αmax=1.0(dB/cm/MHz)、Δα=0.2(dB/cm/MHz)とした場合の減衰率候補値αと分散S(α)との関係の例を示す図である。図11に示す場合、減衰率候補値αが0.2(dB/cm/MHz)のときに分散が最小値S(α)minをとる。したがって、図11に示す場合、減衰率設定部334は、α=0.2(dB/cm/MHz)を最適な減衰率として設定する。
この後、特徴量補正部335は、減衰率設定部334が各分割領域で設定した最適な減衰率を用いて、サンプリング点における累積減衰率を算出する(ステップS13)。例えば、図6に示す分割領域PIJのサンプリング点SIJ(h)における累積減衰率γIJ(h)は、式(7)で表される。
続いて、特徴量補正部335は、累積減衰率を用いて特徴量を減衰補正することによって補正特徴量を算出する(ステップS14)。例えば、図6に示す分割領域PIJのサンプリング点SIJ(h)における傾きa0、切片b0、ミッドバンドフィットc0の補正特徴量aIJ(h)、bIJ(h)、cIJ(h)は、式(9)〜(11)を用いてそれぞれ算出される。
特徴量画像データ生成部342は、Bモード画像データ生成部341が生成したBモード画像データにおける各画素に対し、ステップS14で算出された補正特徴量と関連づけた視覚情報(例えば色相)を重畳することによって特徴量画像データを生成する(ステップS15)。特徴量画像データ生成部342は、生成した特徴量画像データを表示装置4へ送信する。特徴量画像データを受信した表示装置4は、受信した特徴量画像データに対応する特徴量画像を表示する。
ステップS15の後、超音波観測装置3は一連の処理を終了する。なお、超音波観測装置3は、ステップS1〜S15の処理を周期的に繰り返し実行する。
以上説明した本発明の一実施の形態によれば、超音波画像を複数の領域に分割した各分割領域において、互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値に応じた周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに複数の減衰率候補値の中から観測対象に最適な減衰率を分割領域ごとに設定し、超音波振動子の表面とサンプリング点との間に存在する分割領域の最適な減衰率を用いてサンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて特徴量の減衰補正を行うことによって補正特徴量を算出するため、減衰率が不均一な観測対象であっても不均一さを考慮した補正特徴量を算出することができる。したがって、本発明によれば、減衰率が不均一である観測対象の組織性状を精度よく鑑別することが可能となる。
また、本実施の形態によれば、超音波振動子の表面位置とサンプリング点との間に存在する分割領域ごとの最適な減衰率を各分割領域における深度方向の往復距離で重み付けしたものを累積加算することによってサンプリング点における累積減衰率を算出するため、超音波振動子の表面からサンプリング点に至るまでの途中の減衰率を適切に設定することができる。
ここまで、本発明を実施するための形態を説明してきたが、本発明は上述した実施の形態によってのみ限定されるべきものではない。例えば、超音波画像の分割領域の設定方法は、図6に示すものに限られるわけではない。図12は、超音波画像における分割領域の別な設定例を示す図である。図12では、超音波画像201を30個の分割領域Qij(i=1〜5、j=1〜6)に分割した場合を示している。分割領域Qijは、深度が同じ領域すなわちjの値が共通な領域はすべて同じ面積を有している。一方、分割領域Qijは、深度が異なる領域では、超音波振動子21の表面位置202からの深度が大きいほど高さHj(j=1〜6)が大きくなっている(H1<H2<・・・<H6)。一般に、超音波の減衰は深度が大きいほど大きくなるので、超音波振動子21から遠方の分割領域の面積を大きくすることにより、遠方におけるS/N比を向上させることができる。
この場合、分割領域QIJ(1≦I≦imax,1≦J≦jmax;I,Jは定数)で超音波振動子21に近い側の境界からの距離がh(0<h≦HJ)であるようなサンプリング点SIJ(h)における累積減衰率γIJ(h)は、
γIJ(h)=[Σj=1,2,…,J-1(2Hj・α(QIj))]+2h・α(QIJ)
・・・(12)
と表される。この式(12)は、式(7)と比較して右辺第1項の分割領域QIjの深度方向の長さHjが領域に依存している点が異なる。
また、減衰率設定部334が行う最適な減衰率の設定方法は上述した方法に限定されるわけではない。図13は、減衰率設定部334が行う最適な減衰率の別な設定方法の概要を示す図である。図13では、α0=0(dB/cm/MHz)、αmax=1.0(dB/cm/MHz)、Δα=0.2(dB/cm/MHz)とした場合の減衰率候補値αと分散S(α)との関係の例を示しており、減衰率候補値α=0、0.2、0.4、0.6、0.8、1.0(いずれもdB/cm/MHz)における分散S(α)の値は、図11とそれぞれ同じである。この場合には、減衰率設定部334が最適な減衰率を設定する前に、特徴量算出部333が回帰分析を行うことによって減衰率候補値αにおける分散S(α)の値を補間する曲線Rを算出する。その後、減衰率設定部334は、この曲線Rに対し、0(dB/cm/MHz)≦α≦1.0(dB/cm/MHz)における最小値S’(α)minを算出し、そのときの減衰率候補値の値α’を最適な減衰率として設定する。したがって、図13に示す場合、最適な減衰率α’は、0(dB/cm/MHz)と0.2(dB/cm/MHz)の間の値となる。
また、減衰率設定部334は、超音波画像の全てのフレームで最適な減衰率に相当する最適減衰率相当値をそれぞれ算出し、最新のフレームにおける最適減衰率相当値を含む所定数の最適減衰率相当値の平均値、中央値または最頻値を最適な減衰率として設定してもよい。この場合には、各フレームで最適な減衰率を設定する場合と比較して、最適な減衰率の変化が少なくなってその値を安定させることができる。
また、減衰率設定部334は、超音波画像の所定のフレーム間隔で最適な減衰率を設定するようにしてもよい。これにより、計算量を大幅に削減することができる。この場合には、次に最適な減衰率を設定するまでの間、最後に設定した最適な減衰率の値を使用すればよい。
また、入力部35が減衰率候補値の初期値α0の設定変更の入力を受け付け可能な構成としてもよい。
また、統計的なばらつきを与える量として、例えば標準偏差、母集団における特徴量の最大値と最小値の差、特徴量の分布の半値幅のいずれかを適用することも可能である。なお、統計的なばらつきを与える量として分散の逆数を適用する場合も考えられるが、この場合には、その値が最大となる減衰率候補値が最適な減衰率となることはいうまでもない。
また、減衰率設定部334が複数種類の予備補正特徴量の統計的なばらつきをそれぞれ算出し、統計的なばらつきが最小である場合の減衰率候補値を最適な減衰率として設定することも可能である。
また、減衰率設定部334が複数の減衰率候補値を用いて周波数スペクトルを減衰補正し、減衰補正後の周波数スペクトルに対して回帰分析を行うことによって予備補正特徴量を算出するようにしてもよい。
また、超音波内視鏡以外の超音波プローブに対しても適用することが可能である。超音波プローブとして、例えば光学系のない細径の超音波ミニチュアプローブを適用してもよい。超音波ミニチュアプローブは、通常、胆道、胆管、膵管、気管、気管支、尿道、尿管へ挿入され、その周囲臓器(膵臓、肺、前立腺、膀胱、リンパ節等)を観察する際に用いられる。また、超音波プローブとして、被検体の体表から超音波を照射する体外式超音波プローブを適用してもよい。体外式超音波プローブは、通常、腹部臓器(肝臓、胆嚢、膀胱)、乳房(特に乳腺)、甲状腺を観察する際に用いられる。
このように、本発明は、請求の範囲に記載した技術的思想を逸脱しない範囲内において、様々な実施の形態を含みうるものである。