WO2016103839A1 - 超音波診断装置、超音波診断装置の作動方法および超音波診断装置の作動プログラム - Google Patents

Abstract

　本発明にかかる超音波診断装置は、超音波信号の周波数を解析することによって超音波信号の受信深度ごとの周波数スペクトルを算出する周波数解析部と、周波数解析部が算出した周波数スペクトルまたは該周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率を所定の順序で算出することによって周波数スペクトルまたは関数の２次変化率を算出する変化率算出部と、変化率算出部が算出した２次変化率を用いて、超音波振動子の走査領域内の所定領域における超音波信号の単位距離および単位周波数あたりの減衰率を推定する減衰率推定部と、を備える。

Description

超音波診断装置、超音波診断装置の作動方法および超音波診断装置の作動プログラム

　本発明は、被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波診断装置、超音波診断装置の作動方法および超音波診断装置の作動プログラムに関する。

　従来、被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波診断装置において、被検体内での超音波の減衰率を算出する技術が知られている（例えば、特許文献１を参照）。この技術では、超音波探触子から被検体へ向けてガウシアンパルスを送信するとともに、超音波探触子から出力される受信信号を直交検波して、超音波エコーの位相を示す位相φ(t）を取得し、スペックル成分を取り除いた時間ｔの２次微分値ｄ²φ／ｄｔ²を求める。特許文献１では、単位距離当たりの減衰（特許文献１における減衰率）α、２次微分値ｄ²φ／ｄｔ²、音速ｖ、および既知の帯域幅Δｗが満たす関係式
　　ｄ²φ／ｄｔ²＝－２πα×（Δｗ）²×ｖ　　　・・・（１）
を用いることによって減衰率αを算出する（πは円周率）。ここで、特許文献１では、音速ｖを別の方法で得る、とだけ記載されており、その詳細については開示されていない。

　また、ノイズ領域を低Ｓ／Ｎ領域とみなして区別し、この低Ｓ／Ｎ領域の情報を、減衰率に基づく画像である減衰画像とともに表示する超音波診断装置も知られている（例えば、特許文献２を参照）。この技術では、中心周波数が４ＭＨｚの超音波信号を送受信させた場合と同一の条件下で超音波の送信を停止させ、被検体の各位置から受信したノイズ信号に基づいてノイズ画像を生成する。検査中は、ノイズ画像と４ＭＨｚのＢモード画像との輝度を比較し、輝度が同一であるピクセルを低Ｓ／Ｎ領域として抽出し、この低Ｓ／Ｎ領域におけるピクセルの情報と、別に生成した減衰画像とをモニタに表示させる。

特開２０１０－８２２３０号公報 特開２０１３－５８７６号公報

　特許文献１に記載の技術では、音速を別の方法で得るという記載のみが開示されているものの、一般に超音波の経路に応じて異なる音速を正確に求めることは難しい。このため、式（１）を用いて減衰率を精度よく算出することは難しいといわざるを得なかった。

　また、式（１）では、送信波形がガウシアンパルスであることが前提となっているが、超音波振動子のダンパーなどに相当の工夫を施しても、送信波形には尾引きが生じるのが通例である。このため、送信波形がガウシアンパルスであることを前提とするのは不適切であり、式（１）から減衰率を算出しても、高精度の減衰率を算出できるとはいい難かった。なお、特許文献１では、送信波形がガウシアンパルス以外である超音波を被検体へ向けて送信してもよいことが記載されているが、その場合の減衰率の具体的な算出方法は一切開示されていない。

　また、特許文献１に記載の技術では、式（１）を導出する過程で、断層面の注目領域における中心値を採用したり、断層面をスペックル領域と非スペックル領域とに分けたりすることにより、ノイズ成分としてスペックル成分を除くことが記載されている。しかしながら、実際にこの技術を実施する場合、注目領域やスペックル領域等、画像内の領域をその都度マニュアルで指定しなければならず、処理が煩雑であった。また、特許文献１には、２次位相差にてスペックル成分を除く構成および作用効果について詳細な記載がなされているわけではなく、当業者から見ても不明な点が多かった。

　以上のように、特許文献１に記載の技術では、超音波の減衰率を精度よく、かつ簡便に、算出することは困難であった。

　次に、特許文献２に記載の技術が有する問題点を説明する。特許文献２に記載の技術では、送受信条件として中心周波数４ＭＨｚが与えられている。超音波診断装置の送信波形は、中心周波数の周囲成分を含むのが通例である。特に、距離分解能を向上させる場合、パルス幅を短くすべく、送受信帯域を広げる必要がある。この場合、送受信帯域内の全ての周波数成分でＳ／Ｎが十分高いとは限らない。そのため、特許文献２に記載の技術では、Ｓ／Ｎが十分高くない周波数成分を用いることで減衰率の算出精度が低下し、減衰率に基づく画像の信頼性が損なわれてしまうおそれがあった。

　本発明は、上記に鑑みてなされたものであって、超音波の減衰率を精度よくかつ簡便に算出することができるとともに、減衰率に基づく画像の信頼性を向上させることができる超音波診断装置、超音波診断装置の作動方法および超音波診断装置の作動プログラムを提供することを目的とする。

　上述した課題を解決し、目的を達成するために、本発明に係る超音波診断装置は、被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波診断装置であって、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度ごとの周波数スペクトルを算出する周波数解析部と、前記周波数解析部が算出した周波数スペクトルまたは該周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率を所定の順序で算出することによって前記周波数スペクトルまたは前記関数の２次変化率を算出する変化率算出部と、前記変化率算出部が算出した２次変化率を用いて、前記超音波振動子の走査領域内の所定領域における前記超音波信号の単位距離および単位周波数あたりの減衰率を推定する減衰率推定部と、を備えたことを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、周波数および受信深度に応じたノイズレベルデータを記憶するノイズレベルデータ記憶部と、前記ノイズレベルデータ記憶部が記憶するノイズレベルデータと、前記周波数スペクトルのデータとを比較することにより、前記変化率算出部が算出対象とする前記周波数スペクトルの周波数帯域を設定する帯域設定部と、をさらに備え、前記変化率算出部は、前記帯域設定部が設定した周波数帯域における前記周波数変化率を算出することを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記帯域設定部は、同一周波数における前記周波数スペクトルのデータと前記ノイズレベルデータとを比較することによってノイズに相当する周波数の帯域であるノイズ周波数帯域を抽出し、該ノイズ周波数帯域を除外することによって前記周波数帯域を設定することを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記帯域設定部は、所定の周波数より高い周波数帯域を前記ノイズ周波数帯域として抽出することを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記帯域設定部は、前記周波数スペクトルのデータの値が閾値以下の領域を前記ノイズ周波数帯域として抽出することを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記閾値は、前記ノイズレベルデータ以上の値であることを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記変化率算出部は、前記超音波の音線ごとに前記２次変化率を算出し、前記減衰率推定部は、前記音線ごとの前記２次変化率を用いて前記音線ごとの減衰率を算出し、前記走査領域における全ての前記音線ごとの減衰率の統計量を算出することによって前記減衰率を推定することを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記変化率算出部は、前記周波数スペクトルと基準の受信深度における前記周波数スペクトルとの差によって定義される関数に対し、前記周波数帯域で周波数を変数とする回帰分析を行うことによって回帰直線を算出し、前記回帰直線の傾きに対して前記超音波振動子と前記被検体との往復距離を変数とする第２の回帰分析を行うことによって第２の回帰直線を算出し、該第２の回帰直線の傾きを前記２次変化率とすることを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記変化率算出部は、前記周波数帯域に基づいて周波数ごとに定まる距離区間で前記周波数スペクトルに対して前記超音波振動子と前記被検体との往復距離を変数とする回帰分析を行うことによって回帰直線を算出し、該回帰直線の傾きに対して周波数を変数とする第２の回帰分析を行うことによって第２の回帰直線を算出し、該第２の回帰直線の傾きを前記２次変化率とすることを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記超音波画像として信号の振幅を輝度に変換して表示するＢモード画像データを生成するＢモード画像データ生成部と、前記減衰率推定部が推定した減衰率に関する情報と前記Ｂモード画像データ生成部が生成したＢモード画像データとを用いて合成画像データを生成する合成画像データ生成部と、をさらに備えたことを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記減衰率推定部は、前記走査領域を分割することによって得られる複数の部分領域の各々に対し、各部分領域で得られた減衰率の統計量を算出することによって各部分領域の減衰率を推定し、各部分領域の減衰率の値に応じた視覚情報を付与することによって減衰率マップデータを生成する減衰率マップデータ生成部をさらに備えたことを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記超音波画像として信号の振幅を輝度に変換して表示するＢモード画像データを生成するＢモード画像データ生成部と、前記減衰率マップデータ生成部が生成した減衰率マップデータと前記Ｂモード画像データ生成部が生成したＢモード画像データとを用いて合成画像データを生成する合成画像データ生成部と、をさらに備えたことを特徴とする。

　本発明に係る超音波診断装置は、上記発明において、前記ノイズレベルデータ記憶部は、当該超音波診断装置に接続可能な前記超音波プローブが備える前記超音波振動子に応じたノイズレベルデータを記憶することを特徴とする。

　本発明に係る超音波診断装置の作動方法は、被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波診断装置の作動方法であって、周波数解析部が前記超音波信号の周波数を解析することによって前記超音波信号の受信深度ごとの周波数スペクトルを算出する周波数解析ステップと、変化率算出部が前記周波数スペクトルまたは該周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率を所定の順序で算出することによって前記周波数スペクトルまたは前記関数の２次変化率を算出する変化率算出ステップと、減衰率推定部が前記２次変化率を用いて、前記超音波振動子の走査領域内の所定領域における前記超音波信号の単位距離および単位周波数あたりの減衰率を推定する減衰率推定ステップと、を有することを特徴とする。

　本発明に係る超音波診断装置の作動プログラムは、被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波診断装置に、周波数解析部が前記超音波信号の周波数を解析することによって前記超音波信号の受信深度ごとの周波数スペクトルを算出する周波数解析ステップと、変化率算出部が前記周波数スペクトルまたは該周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率を所定の順序で算出することによって前記周波数スペクトルまたは前記関数の２次変化率を算出する変化率算出ステップと、減衰率推定部が前記２次変化率を用いて、前記超音波振動子の走査領域内の所定領域における前記超音波信号の単位距離および単位周波数あたりの減衰率を推定する減衰率推定ステップと、を実行させることを特徴とする。

　本発明によれば、超音波の減衰率を精度よくかつ簡便に算出することができるとともに、減衰率に基づく画像の信頼性を向上させることができる。

図１は、本発明の実施の形態１に係る超音波診断装置を備えた超音波診断システムの機能構成を示すブロック図である。 図２は、本発明の実施の形態１に係る超音波診断装置が備えるノイズレベルデータ記憶部が記憶するノイズレベルデータを模式的に示す図である。 図３は、本発明の実施の形態１に係る超音波診断装置が備える信号増幅部が行う増幅処理における受信深度と増幅率との関係を示す図である。 図４は、超音波振動子の走査領域とＢモード用受信データとを模式的に示す図である。 図５は、本発明の実施の形態１に係る超音波診断装置が備える増幅補正部が行う増幅補正処理における受信深度と増幅率との関係を示す図である。 図６は、超音波信号の１つの音線におけるデータ配列を模式的に示す図である。 図７は、周波数スペクトルのデータ列を模式的に示す図である。 図８は、周波数スペクトルのデータの具体例を示す図である。 図９は、周波数スペクトルのデータとノイズレベルデータとの関係を示す図である。 図１０は、本発明の実施の形態１に係る超音波診断装置が備える帯域設定部が設定した回帰分析対象の周波数帯域を視覚的に表現した周波数スペクトルのデータを模式的に示す図である。 図１１は、代表的な距離の関数と、各関数の回帰直線を示す図である。 図１２は、回帰直線の傾きと距離の関係を示す図である。 図１３は、本発明の実施の形態１に係る超音波診断装置が実行する処理の概要を示すフローチャートである。 図１４は、表示装置が表示する合成画像の表示例を示す図である。 図１５は、本発明の実施の形態１に係る超音波診断装置が行う周波数解析処理の概要を示すフローチャートである。 図１６は、周波数が一定の周波数スペクトルのデータと距離の関係を示す図である。 図１７は、周波数スペクトルのデータの回帰直線の傾きと周波数との関係を示す図である。 図１８は、本発明の実施の形態３に係る超音波診断システムの機能構成を示すブロック図である。 図１９は、本発明の実施の形態３に係る超音波診断装置が実行する処理の概要を示すフローチャートである。 図２０は、超音波振動子の走査領域における部分領域の設定例を示す図である。 図２１は、表示装置が表示する減衰率マップデータ付きの合成画像の表示例を示す図である。

　以下、添付図面を参照して、本発明を実施するための形態（以下、「実施の形態」という）を説明する。

（実施の形態１）
　図１は、本発明の実施の形態１に係る超音波診断装置を備えた超音波診断システムの機能構成を示すブロック図である。同図に示す超音波診断システム１は、被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波内視鏡２と、超音波内視鏡２が取得した超音波信号に基づいて超音波画像を生成する超音波診断装置３と、超音波診断装置３が生成した超音波画像を表示する表示装置４と、を備える。

　超音波内視鏡２は、その先端部に、超音波診断装置３から受信した電気的なパルス信号を超音波パルス（音響パルス）に変換して被検体へ照射するとともに、被検体で反射された超音波エコーを電圧変化で表現する電気的なエコー信号に変換して出力する超音波振動子２１を有する。超音波内視鏡２は、超音波振動子２１をメカ的に走査させるものであってもよいし、超音波振動子２１として複数の素子をアレイ状に設け、送受信にかかわる素子を電子的に切り替えたり、各素子の送受信に遅延をかけたりすることで、電子的に走査させるものであってもよい。

　超音波内視鏡２は、通常は撮像光学系および撮像素子を有しており、被検体の消化管（食道、胃、十二指腸、大腸）または呼吸器（気管、気管支）へ挿入され、消化管、呼吸器やその周囲臓器（膵臓、胆嚢、胆管、胆道、リンパ節、縦隔臓器、血管等）を撮像することが可能である。また、超音波内視鏡２は、撮像時に被検体へ照射する照明光を導くライトガイドを有する。このライトガイドは、先端部が超音波内視鏡２の被検体への挿入部の先端まで達している一方、基端部が照明光を発生する光源装置に接続されている。

　超音波診断装置３は、超音波診断システム１全体を制御する制御部３１と、超音波診断装置３の動作に必要な各種情報を記憶する記憶部３２と、超音波内視鏡２と電気的に接続され、所定の波形および送信タイミングに基づいて高電圧パルスからなる送信信号（パルス信号）を超音波振動子２１へ送信するとともに、超音波振動子２１から電気的な受信信号であるエコー信号を受信してデジタルの高周波（ＲＦ：Radio　Frequency）信号のデータ（以下、ＲＦデータという）を生成、出力する送受信部３３と、送受信部３３から受信したＲＦデータをもとにデジタルのＢモード用受信データを生成する信号処理部３４と、送受信部３３から受信したＲＦデータに対して所定の演算を施す演算部３５と、各種画像データを生成する画像処理部３６と、キーボード、マウス、タッチパネル等のユーザインタフェースを用いて実現され、各種情報の入力を受け付ける入力部３７と、を備える。

　制御部３１は、超音波診断装置３に接続される超音波内視鏡２の種別を判定する種別特定部３１１を有する。種別特定部３１１は、超音波診断装置３に接続された超音波内視鏡２内のメモリに格納されているＩＤを取得することによって超音波振動子２１の種別を特定する。種別特定部３１１が特定した超音波振動子２１の種別に関する情報は、後述する記憶部３２の種別情報記憶部３２１に格納される。なお、種別特定部３１１は、入力部３７が入力を受け付けた超音波内視鏡２の種別名に基づいて超音波振動子２１の種別を特定するようにしてもよい。

　制御部３１は、演算および制御機能を有するＣＰＵ（Central　Processing　Unit）や各種演算回路等を用いて実現される。制御部３１は、記憶部３２が記憶、格納する情報および超音波診断装置３の作動プログラムを含む各種プログラムを記憶部３２から読み出し、超音波診断装置３の作動方法に関連した各種演算処理を実行することによって超音波診断装置３を統括して制御する。

　記憶部３２は、種別特定部３１１が特定した超音波振動子２１の種別情報を記憶する種別情報記憶部３２１と、超音波診断装置３に接続可能な超音波振動子２１の種別ごとにノイズレベルデータを記憶するノイズレベルデータ記憶部３２２と、演算部３５が行う回帰分析の対象となる周波数帯域の情報を超音波診断装置３に接続可能な超音波振動子２１の種別ごとに記憶する帯域情報記憶部３２３とを有する。

　図２は、ノイズレベルデータ記憶部３２２が記憶するノイズレベルデータを模式的に示す図である。図２において、ドットで示す曲面がノイズレベルデータｎ（ｆ，Ｌ）を与える。ノイズレベルデータｎ（ｆ，Ｌ）は、超音波エコーの周波数ｆ、および超音波振動子２１の表面と被検体（反射体）との往復距離Ｌ（受信深度の２倍に相当する距離）を離散的な変数とする関数である。ノイズレベルデータｎ（ｆ，Ｌ）は、例えばノイズに相当する電圧Ｖを基準電圧Ｖ_cで除した量の常用対数をとってデシベル値（ｄＢ）で表現したものであり、離散的なデジタルデータである。なお、以下では、超音波振動子２１の表面と被検体との往復距離Ｌを、単に距離Ｌということもある。

　ノイズレベルデータ記憶部３２２は、超音波診断装置３に接続可能な超音波振動子２１の種別ごとのノイズレベルデータを記憶している。これらのノイズレベルデータは、例えばあらかじめ工場出荷時に超音波振動子２１の種別ごとに測定されて格納される。なお、送受信部３３の送信を中止したときに受信したエコー信号をノイズレベルデータとして記憶するようにしてもよい。この場合には、超音波振動子２１の種別だけでなく、同一種類の超音波振動子２１における個体差や、同一被検体の経時変化によるノイズレベルの差異も考慮した、より正確なノイズレベルデータを取得することができる。

　帯域情報記憶部３２３は、演算部３５による回帰分析対象の帯域情報として、超音波診断装置３に接続可能な超音波振動子２１の種別ごとの周波数帯域の最小値ｆminおよび最大値ｆmaxに関する情報を記憶している。帯域情報は、ノイズレベルデータと同様、あらかじめ工場出荷時に超音波振動子２１の種別ごとに測定されて格納される。

　記憶部３２は、上記以外にも、例えば送受信部３３、信号処理部３４および演算部３５が行う各種処理に必要な情報を記憶している。

　記憶部３２は、超音波診断装置３の作動方法を実行するための作動プログラムを含む各種プログラムを記憶する。各種プログラムは、ハードディスク、フラッシュメモリ、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ、フレキシブルディスク等のコンピュータ読み取り可能な記録媒体に記録して広く流通させることも可能である。なお、上述した各種プログラムは、通信ネットワークを介してダウンロードすることによって取得することも可能である。ここでいう通信ネットワークは、例えば既存の公衆回線網、ＬＡＮ（Local　Area　Network）、ＷＡＮ（Wide　Area　Network）などによって実現されるものであり、有線、無線を問わない。

　以上の構成を有する記憶部３２は、各種プログラム等が予めインストールされたＲＯＭ（Read　Only　Memory）、および各処理の演算パラメータやデータ等を記憶するＲＡＭ（Random　Access　Memory）およびＨＤＤ（Hard　Disk　Drive）等を用いて実現される。

　送受信部３３は、エコー信号を増幅する信号増幅部３３１を有する。信号増幅部３３１は、受信深度が大きいエコー信号ほど高い増幅率で増幅するＳＴＣ（Sensitivity　Time　Control）補正を行う。図３は、信号増幅部３３１が行う増幅処理における受信深度と増幅率との関係を示す図である。図３に示す受信深度ｚは、超音波の受信開始時点からの経過時間に基づいて算出される量である。図３に示すように、増幅率β（ｄＢ）は、受信深度ｚが閾値ｚ_thより小さい場合、受信深度ｚの増加に伴ってβ_０からβ_th（＞β₀）へ線型に増加する。また、増幅率β（ｄＢ）は、受信深度ｚが閾値ｚ_th以上である場合、一定値β_thをとる。閾値ｚ_thの値は、被検体から受信する超音波信号がほとんど減衰してしまい、ノイズが支配的になるような値である。より一般に、増幅率βは、受信深度ｚが閾値ｚ_thより小さい場合、受信深度ｚの増加に伴って単調増加すればよい。なお、図３に示す関係は、予め記憶部３２に記憶されている。

　送受信部３３は、信号増幅部３３１によって増幅されたエコー信号に対してフィルタリング等の処理を施した後、Ａ／Ｄ変換することによって時間ドメインのＲＦデータを生成し、信号処理部３４および演算部３５へ出力する。なお、超音波内視鏡２が複数の素子をアレイ状に設けた超音波振動子２１を電子的に走査させる構成を有する場合、送受信部３３は、複数の素子に対応したビーム合成用の多チャンネル回路を有する。

　送受信部３３が送信するパルス信号の周波数帯域は、超音波振動子２１におけるパルス信号の超音波パルスへの電気音響変換の線型応答周波数帯域をほぼカバーする広帯域にするとよい。それにより、後述する周波数スペクトルの近似処理において、精度のよい近似を行うことが可能となる。

　送受信部３３は、制御部３１が出力する各種制御信号を超音波内視鏡２に対して送信するとともに、超音波内視鏡２から識別用のＩＤを含む各種情報を受信して制御部３１へ送信する機能も有する。

　信号処理部３４は、ＲＦデータに対してバンドパスフィルタ、包絡線検波、対数変換など公知の処理を施し、デジタルのＢモード用受信データを生成する。対数変換では、ＲＦデータを基準電圧Ｖ_cで除した量の常用対数をとってデシベル値で表現する。このＢモード用受信データでは、超音波パルスの反射の強さを示す受信信号の振幅または強度が、超音波パルスの送受信方向（深度方向）に沿って並んでいる。図４は、超音波振動子２１の走査領域（以下、単に走査領域ということもある）とＢモード用受信データとを模式的に示す図である。図４に示す走査領域Ｓは扇形をなしている。これは、超音波振動子２１がコンベックス振動子である場合に相当している。また、図４では、Ｂモード用受信データの受信深度をｚとして記載している。超音波振動子２１の表面から照射された超音波パルスが受信深度ｚにある反射体で反射し、超音波エコーとして超音波振動子２１へ戻ってきた場合、その往復距離Ｌと受信深度ｚとの間には、上述したようにｚ＝Ｌ／２の関係がある。信号処理部３４は、生成したＢモード用受信データを、画像処理部３６のＢモード画像データ生成部３６１へ出力する。信号処理部３４は、ＣＰＵ（Central　Processing　Unit）や各種演算用の回路等を用いて実現される。

　一方、演算部３５は、送受信部３３が出力したＲＦデータに対して受信深度によらず増幅率を一定とするよう増幅補正を行う増幅補正部３５１と、増幅補正を行ったＲＦデータに高速フーリエ変換（ＦＦＴ：Fast　Fourier　Transform）を施して周波数解析を行うことにより受信深度ごとの周波数スペクトルを算出する周波数解析部３５２と、ノイズレベルデータと周波数スペクトルのデータとを比較することにより演算対象とする周波数帯域を設定する帯域設定部３５３と、周波数解析部３５２が算出した周波数スペクトルにより定義される関数の距離変化率および周波数変化率を所定の順序で算出することによって２次変化率を算出する変化率算出部３５４と、変化率算出部３５４が算出した２次変化率を用いて走査領域における超音波パルスの単位距離および単位周波数あたりの減衰率を推定する減衰率推定部３５５と、を有する。演算部３５は、ＣＰＵや各種演算用の回路等を用いて実現される。なお、演算部３５を、制御部３１および信号処理部３４と共通のＣＰＵ等を用いて構成することも可能である。

　図５は、増幅補正部３５１が行う増幅補正処理における受信深度と増幅率との関係を示す図である。図５に示すように、増幅補正部３５１が行う増幅処理における増幅率β（ｄＢ）は、受信深度ｚがゼロのとき最大値β_th－β₀をとり、受信深度ｚがゼロから閾値ｚ_thに達するまで線型に減少し、受信深度ｚが閾値ｚ_th以上のときゼロである。なお、図５に示す関係は、予め記憶部３２に記憶されている。増幅補正部３５１が図５に示す関係に基づいてＲＦデータを増幅補正することにより、信号増幅部３３１におけるＳＴＣ補正の影響を相殺し、一定の増幅率β_thの信号を出力することができる。なお、増幅補正部３５１が行う増幅補正処理における受信深度ｚと増幅率βの関係は、信号増幅部３３１が行う増幅補正処理における受信深度と増幅率の関係に応じて異なることは勿論である。

　このような増幅補正を行う理由を説明する。ＳＴＣ補正は、アナログ信号波形の振幅を全周波数帯域にわたって均一に、かつ、受信深度の増加に対して単調増加する増幅率で増幅させることで、アナログ信号波形の振幅から減衰の影響を排除する補正処理である。このため、エコー信号の振幅を輝度に変換して表示するＢモード画像を生成する場合、かつ、一様な組織を走査した場合には、ＳＴＣ補正を行うことによって深度によらず輝度値が一定になる。すなわち、Ｂモード画像の輝度値から減衰の影響を排除する効果を得ることができる。一方、本実施の形態１のように超音波の周波数スペクトルを算出して解析した結果を利用する場合、ＳＴＣ補正でも超音波の伝播に伴う減衰の影響を正確に排除できるわけではない。何故なら、減衰量は周波数によって異なるが、ＳＴＣ補正の増幅率は距離だけに対して変化し、周波数に対しては変化せず一定であるためである。

　上述した問題、すなわち超音波の周波数スペクトルを算出して解析した結果を利用する場合、ＳＴＣ補正でも超音波の伝播に伴う減衰の影響を正確に排除できるわけではない、という問題を解決するには、Ｂモード画像を生成する際にＳＴＣ補正を施した受信信号を出力する一方、周波数スペクトルに基づいた画像を生成する際に、Ｂモード画像を生成するための送信とは異なる新たな送信を行い、ＳＴＣ補正を施していない受信信号を出力することが考えられる。ところがこの場合には、受信信号に基づいて生成される画像データのフレームレートが低下してしまうという問題がある。

　そこで、本実施の形態１では、生成される画像データのフレームレートを維持しつつ、Ｂモード画像用にＳＴＣ補正を施した信号に対してＳＴＣ補正の影響を排除するために、増幅補正部３５１によって増幅率の補正を行う。

　周波数解析部３５２は、増幅補正部３５１が増幅補正した各音線のＲＦデータ（ラインデータ）を所定の時間間隔でサンプリングし、サンプルデータを生成する。そして、周波数解析部３５２は、サンプルデータ群にＦＦＴ処理を施すことにより、ＲＦデータ上の複数の箇所（データ位置）における周波数スペクトルを算出する。

　図６は、超音波信号の１つの音線におけるデータ配列を模式的に示す図である。同図に示す音線ＳＲ_kにおいて、白または黒の長方形は、１つのサンプル点におけるデータを意味している。また、音線ＳＲ_kにおいて、右側に位置するデータほど、超音波振動子２１から音線ＳＲ_kに沿って計った場合の深い箇所からのサンプルデータである（図６の矢印を参照）。音線ＳＲ_kは、送受信部３３が行うＡ／Ｄ変換におけるサンプリング周波数（例えば５０ＭＨｚ）に対応した時間間隔で離散化されている。図６では、番号ｋの音線ＳＲ_kの８番目のデータ位置を受信深度ｚの方向の初期値Ｚ^(k) ₀として設定した場合を示しているが、初期値の位置は任意に設定することができる。周波数解析部３５２による算出結果は複素数で得られ、記憶部３２に格納される。

　図６に示すデータ群Ｆ_j（ｊ＝１，２，・・・，Ｋ）は、ＦＦＴ処理の対象となるサンプルデータ群である。一般に、ＦＦＴ処理を行うためには、サンプルデータ群が２のべき乗のデータ数を有している必要がある。この意味で、サンプルデータ群Ｆ_j（ｊ＝１,２，・・・，Ｋ－１）はデータ数が１６（＝２⁴）で正常なデータ群である一方、サンプルデータ群Ｆ_Kは、データ数が１２であるため異常なデータ群である。異常なデータ群に対してＦＦＴ処理を行う際には、不足分だけゼロデータを挿入することにより、正常なサンプルデータ群を生成する処理を行う。この点については、後述する周波数解析部３５２の処理を説明する際に詳述する（図１５を参照）。

　周波数解析部３５２は、ＲＦデータから切り取ったサンプルデータ群の各々に対してＦＦＴ処理を施すことによって電圧振幅の周波数成分Ｖ（ｆ，Ｌ）を生成する。この電圧振幅の周波数成分Ｖ（ｆ，Ｌ）は電圧の周波数密度である。さらに、周波数解析部３５２は、電圧振幅の周波数成分Ｖ（ｆ，Ｌ）を基準電圧Ｖ_cで除し、常用対数（ｌｏｇ）をとってデシベル単位で表現する対数変換処理を施した後、適当な定数Ａを乗ずることにより、次式（２）で与えられる周波数スペクトルのデータ（以下、スペクトルデータともいう）Ｆ（ｆ，Ｌ）を生成し、帯域設定部３５３へ出力する。
　　Ｆ（ｆ，Ｌ）＝Ａ・ｌｏｇ｛Ｖ（ｆ，Ｌ）／Ｖ_c｝　　・・・（２）
ここで、ｌｏｇは常用対数である（以下、同じ）。

　スペクトルデータＦ（ｆ，Ｌ）は、サンプルデータ群の周波数ｆの成分である。図７は、記憶部３２が記憶するスペクトルデータのデータ列を模式的に示す図である。図７において、縦方向は周波数ｆを示し、横方向は超音波振動子２１の表面からの往復距離Ｌを示している。周波数ｆは、離散的な値０、Δｆ、２Δｆ、・・・を取る。例えば、列０のセルには、距離区間０≦Ｌ＜ΔＬで切り取られたサンプルデータ群に基づいて式（２）から得られるスペクトルデータＦ（ｆ，０）が格納される。また、列ΔＬのセルには、距離Ｌの区間ΔＬ≦Ｌ＜２ΔＬで切り取られたサンプルデータ群に基づいて式（２）から得られるスペクトルデータＦ（ｆ，ΔＬ）が格納される。図７では、周波数Δｆ、距離ΔＬのセルのみにスペクトルデータＦ（Δｆ，ΔＬ）のみを例示的に記載しているが、実際には全てのセルに、各セルの周波数および距離に応じたスペクトルデータが格納されていることはいうまでもない。なお、サンプルデータ群を切り取る区間の長さΔＬ（図６のサンプルデータ群のステップ幅Ｄに相当）は、例えば１．０ｃｍ程度である。また、周波数の変化量Δｆは、例えば０．５ＭＨｚである。

　図８は、スペクトルデータの具体例を示す図である。図８では、互いに異なる４つの距離におけるスペクトルデータＦ（ｆ，Ｌ１）、Ｆ（ｆ，Ｌ２）、Ｆ（ｆ，Ｌ３）、Ｆ（ｆ，Ｌ４）と周波数ｆとの関係を示している。ここで、４つの距離Ｌ１、Ｌ２、Ｌ３、Ｌ４は定数であり、０＜Ｌ１＜Ｌ２＜Ｌ３＜Ｌ４を満たしている。図７からも明らかなように、実際にはさらに多くのスペクトルデータＦ（ｆ，Ｌ）が算出されるが、図８では、代表的な４つのスペクトルデータのみを例示している。図８に示すように、スペクトルデータＦ（ｆ，Ｌ）は、距離Ｌが大きくなるにつれて減少していく。また、スペクトルデータＦ（ｆ，Ｌ）の平均周波数は、距離Ｌが大きくなるにつれて低周波数側へシフトしていく。これは、被検体内での超音波が伝播する際の周波数に依存した減衰の効果によるものである。

　一般に、スペクトルデータＦ（ｆ，Ｌ）は、超音波が走査された組織の属性によって異なる傾向を示す。これは、スペクトルデータＦ（ｆ，Ｌ）が、超音波を散乱する散乱体の大きさ、数密度、音響インピーダンス等と相関を有しているためである。ここでいう「属性」とは、例えば悪性腫瘍組織、良性腫瘍組織、内分泌腫瘍組織、粘液性腫瘍組織、正常組織、嚢胞、脈管などのことである。

　帯域設定部３５３は、スペクトルデータＦ（ｆ，Ｌ）とノイズレベルデータｎ（ｆ，Ｌ）とを比較する比較部３５３ａを有する。比較部３５３ａは、上述した二つのデータを比較するために、超音波診断装置３に接続されている超音波振動子２１に応じて、後述する回帰分析の対象とする周波数帯域Ｕ＝｛ｆ｜ｆmin≦ｆ≦ｆmax｝の両端の周波数ｆmin、ｆmaxを帯域情報記憶部３２３から読み出す。周波数帯域Ｕは、超音波振動子２１の表面（Ｌ＝０）における超音波の送信波形上で比較的平坦な区間に相当し、最小値ｆminおよび最大値ｆmaxは、超音波振動子２１の種別に応じて異なる。

　比較部３５３ａは、帯域情報記憶部３２３から読み出した周波数帯域Ｕにおいて、スペクトルデータＦ（ｆ，Ｌ）と、ノイズレベルデータｎ（ｆ，Ｌ）とを周波数ｆ，距離Ｌごとにそれぞれ比較する。図８に示す４つのスペクトルデータＦ（ｆ，Ｌｐ）（ｐ＝１，２，３，４）を対応するノイズレベルデータｎ（ｆ，Ｌｐ）と比較した場合、例えば、ｐ＝１，２，３のとき、周波数帯域Ｕに含まれる任意の周波数ｆでＦ（ｆ，Ｌｐ）＞ｎ（ｆ，Ｌｐ）が成立する一方、ｐ＝４のとき、周波数帯域Ｕでは以下の２つの不等式が成立したとする。
　　Ｆ（ｆ，Ｌ４）＞ｎ（ｆ，Ｌ４）　（ｆmin≦ｆ＜ｆmax’）
　　Ｆ（ｆ，Ｌ４）≦ｎ（ｆ，Ｌ４）　（ｆmax’≦ｆ≦ｆmax）

　図９は、スペクトルデータＦ（ｆ，Ｌ４）とノイズレベルデータｎ（ｆ，Ｌ４）との関係を示す図である。図９に示す場合、２つの周波数帯域ｆ≦ｆmin’、ｆ≧ｆmax’においてＦ（ｆ，Ｌ４）≦ｎ（ｆ，Ｌ４）が成立し、ノイズレベルデータが支配的である。以下、この２つの周波数帯域のようにノイズレベルデータが支配的である周波数帯域を、ノイズ周波数帯域という。図９より、周波数ｆmin’、ｆmax’のうち、周波数帯域Ｕに含まれるのはｆmax’である。

　帯域設定部３５３は、比較部３５３ａの比較結果に基づいて回帰分析対象の周波数帯域を設定する。図８および図９に示す場合、帯域設定部３５３は、距離Ｌｐ（ｐ＝１、２、３）に対する回帰分析対象を初期の周波数帯域Ｕと設定する一方、距離Ｌ＝Ｌ４の回帰分析対象の周波数帯域をＵ’＝｛ｆ｜ｆmin≦ｆ＜ｆmax’｝と設定し、これらの周波数帯域情報を減衰率推定部３５５へ出力する。

　図１０は、帯域設定部３５３が設定した回帰分析対象の周波数帯域を視覚的に表現したスペクトルデータのデータ列を模式的に示す図である。図１０においても、各セルのスペクトルデータＦ（ｆ，Ｌ）の記載は省略している。図１０では、スペクトルデータＦ（ｆ，Ｌ）の種類に応じてセルごとに３種類の異なる模様のいずれかを付与している。白色のセルに格納されるスペクトルデータＦ（ｆ，Ｌ）は、周波数ｆが周波数帯域Ｕに含まれるとともに、Ｆ（ｆ，Ｌ）＞ｎ（ｆ，Ｌ）を満たす値を有することを示している。また、斜線が記載されたセルに格納されるスペクトルデータＦ（ｆ，Ｌ）は、周波数ｆが周波数帯域Ｕの外の周波数であるとともに、Ｆ（ｆ，Ｌ）＞ｎ（ｆ，Ｌ）を満たす値を有することを示している。さらに、ドットが記載されたセルに格納されるスペクトルデータＦ（ｆ，Ｌ）は、周波数ｆが周波数帯域Ｕの外の周波数であるとともに、Ｆ（ｆ，Ｌ）≦ｎ（ｆ，Ｌ）を満たす値を有することを示している。上述したように、距離Ｌ＝Ｌｐ（ｐ＝１，２，３）における周波数帯域はＵであり、距離Ｌ＝Ｌ４における周波数帯域はＵ’である。なお、図１０に記載されたＬｍａｘ（ｆｐ）（ｐ＝１，２，３，４）については後述する。

　変化率算出部３５４は、スペクトルデータＦ（ｆ，Ｌ）により定義される数であって周波数ｆの１次関数である関数ψ（ｆ，Ｌ）＝Ｆ（ｆ，Ｌ）－Ｆ（ｆ，０）に対して周波数ｆおよび距離Ｌに対する変化率を、周波数から距離の順で順次算出する。

　減衰率推定部３５５は、変化率算出部３５４が算出した２次変化率を用いて、走査領域における超音波パルスの単位距離および単位周波数あたりの減衰率を推定する。

　ここで、変化率算出部３５４および減衰率推定部３５５が行う処理をより詳細に説明する。まず、スペクトルデータＦ（ｆ，Ｌ）の単位距離および単位周波数あたりの減衰率の推定方法を説明する。距離Ｌに存在する反射体からの超音波の周波数ｆにおける音圧振幅Ｐ（ｆ，Ｌ）は、正の定数μを用いて、
　Ｐ（ｆ，Ｌ）＝Ｐ（ｆ，０）・ｅｘｐ（－μｆＬ）　　　・・・（３）
で与えられることが知られている。μｆＬ＞０なので、式（３）は、音圧振幅Ｐ（ｆ，Ｌ）が、周波数ｆおよび距離Ｌの増加に対して指数関数的に減衰することを意味している。

　一方、距離区間Ｌ～Ｌ＋ΔＬの周波数ｆにおける超音波の減衰量をＬｏｓｓ（ｆ，Ｌ）［ｄＢ］とすると、この減衰量は、
　　Ｌｏｓｓ（ｆ，Ｌ）＝Ａ・ｌｏｇ｛Ｐ（ｆ，Ｌ）／Ｐ（ｆ，Ｌ＋ΔＬ）｝
　　　　＝Ａ・ｌｏｇＰ（ｆ，Ｌ）－Ａ・ｌｏｇＰ（ｆ，Ｌ＋ΔＬ）　　・・・（４）
で定義される。ここで、右辺の定数Ａは、式（２）の定数Ａと同じである。したがって、単位距離および単位周波数あたりの減衰率ζは、次式で与えられる。
　　ζ＝（∂／∂ｆ）Ｌｉｍ｛Ｌｏｓｓ（ｆ，Ｌ）／ΔＬ｝
　　　＝(∂／∂ｆ）｛－Ａ(∂／∂Ｌ）ｌｏｇＰ（ｆ，Ｌ）｝
　　　＝－Ａ（∂²／∂ｆ∂Ｌ）ｌｏｇＰ（ｆ，Ｌ）　　　・・・（５）
ここで、Ｌｉｍ｛Ｌｏｓｓ（ｆ，Ｌ）／ΔＬ｝は、関数Ｌｏｓｓ（ｆ，Ｌ）／ΔＬのΔＬ→０における極限を意味する。単位距離および単位周波数あたりの減衰率ζの単位は、例えば［ｄＢ／ｃｍ／ＭＨｚ］である。以下、単位距離および単位周波数あたりの減衰率を単に減衰率ということもある。

　なお、上記定数μと減衰率ζとの関係は次の通りである。式（５）のＰ（ｆ，Ｌ）に式（３）を代入すると、減衰率ζは、
　　ζ＝－Ａ（∂²／∂ｆ∂Ｌ）［ｌｏｇ｛Ｐ（ｆ，０）・ｅｘｐ（－μｆＬ）｝］
　　　＝－Ａ（∂²／∂ｆ∂Ｌ）｛ｌｏｇＰ（ｆ，０）－μｆＬｌｏｇｅ｝
　　　＝（ｌｏｇｅ）Ａμ　　　・・・（６）
ここで、ｅは自然対数の底である。

　さて、超音波振動子２１の感度を周波数ｆの関数としてγ（ｆ）とすると、ＲＦデータにＦＦＴ処理を施した後の振幅成分Ｖ（ｆ，Ｌ）は、次式（７）で与えられる。
　　Ｖ（ｆ，Ｌ）＝γ（ｆ）・Ｐ（ｆ，Ｌ）　　　・・・（７）
この式（７）のＰ（ｆ，Ｌ）に式（３）を代入すると、
　　Ｖ（ｆ，Ｌ）＝γ（ｆ）・Ｐ（ｆ，０）・ｅｘｐ（－μｆＬ）
　　　　　　　　＝Ｖ（ｆ，０）・ｅｘｐ（－μｆＬ）　　　・・・（８）
が得られる。

　式（８）を式（２）へ代入することにより、
　　Ｆ（ｆ，Ｌ）＝Ａ・ｌｏｇ｛Ｖ（ｆ，０）・ｅｘｐ(－μｆＬ）／Ｖ_c｝
　　＝Ａｌｏｇ・ｅｘｐ（－μｆＬ）＋Ａｌｏｇ｛Ｖ（ｆ，０）／Ｖ_c｝
　　＝－（ｌｏｇｅ）ＡμｆＬ＋Ｆ（ｆ，０）　　　・・・（９）
が得られる。さらに、式（９）の右辺に式（６）を代入することにより、
　　Ｆ（ｆ，Ｌ）－Ｆ（ｆ，０）＝－ζｆＬ　　　・・・（１０）
が導出される。

　式（１０）の両辺に２階偏微分演算子∂²／∂Ｌ∂ｆ、および∂²／∂ｆ∂Ｌをそれぞれ作用させると、次式が得られる。
　ζ＝－∂²Ｆ（ｆ，Ｌ）／∂Ｌ∂ｆ＝－∂²Ｆ（ｆ，Ｌ）／∂ｆ∂Ｌ　　・・・（１１）
ここで、∂²／∂Ｌ∂ｆは周波数ｆの偏微分を先に行う一方、∂²／∂ｆ∂Ｌは距離Ｌの偏微分を先に行うことを意味する。

　したがって、スペクトルデータＦ（ｆ，Ｌ）の２次偏導関数∂²Ｆ（ｆ，Ｌ）／∂ｆ∂Ｌまたは∂²Ｆ（ｆ，Ｌ）／∂Ｌ∂ｆを演算することにより、減衰率ζを推定することができる。

　以上説明した減衰率の推定方法において、スペクトルデータＦ（ｆ，Ｌ）の偏微分を演算することは実際には困難である場合が多い。何故なら、偏微分の定義に従うと、偏微分の演算の際には、極限Δｆ→０、ΔＬ→０（Δｆ、ΔＬはそれぞれｆ、Ｌの微小変位）を計算する必要があるが、実際のスペクトルデータＦ（ｆ，Ｌ）は離散的に定義され、これらの極限を計算することは困難であるからである。この問題を回避するために、周波数ｆや距離Ｌの隣接する離散値の差分をとることによってスペクトルデータＦ（ｆ，Ｌ）の偏微分の演算を近似する方法が知られている。ところがこの方法では、偏導関数がスペクトルデータＦ（ｆ，Ｌ）のゆらぎに起因するノイズを多く含んでしまう可能性がある。

　本実施の形態１において、変化率算出部３５４は、スペクトルデータＦ（ｆ，Ｌ）の関数ψ（ｆ，Ｌ）＝Ｆ（ｆ，Ｌ）－Ｆ（ｆ，０）に対して回帰分析を実施し、回帰直線による近似を実施する。この関数ψ（ｆ，Ｌ）の２次偏導関数は、
　　∂²ψ（ｆ，Ｌ）／∂Ｌ∂ｆ＝∂²Ｆ（ｆ，Ｌ）／∂Ｌ∂ｆ　　・・・（１２）
となるので、式（１２）の右辺に式（１１）を代入すると、
　　∂²ψ（ｆ，Ｌ）／∂Ｌ∂ｆ＝－ζ　　・・・（１３）
が得られる。同様に、
　　∂²ψ（ｆ，Ｌ）／∂ｆ∂Ｌ＝－ζ　　・・・（１４）
も得られる。式（１３）、（１４）は、関数ψ（ｆ，Ｌ）を用いて減衰率ζを算出できることを示している。

　ところで、関数ψ（ｆ，Ｌ）は、式（１０）により、
　　ψ（ｆ，Ｌ）＝－ζｆＬ　　・・・（１５）
と表される。ところで、回帰直線による近似は、結局、一次関数での近似である。よって、近似する関数が一次関数に近いほど、回帰直線はその近似する関数に近づき、良い近似を与える。ここで、式（１５）に示す通り、関数ψ（ｆ，Ｌ）は周波数ｆの一次関数そのものである。一方、スペクトルデータＦ（ｆ，Ｌ）は周波数ｆの一次関数に近いとは限らない。したがって、周波数ｆに対する関数の偏導関数を関数の回帰直線の傾き（すなわち変化率）で近似する場合、関数ψ（ｆ，Ｌ）を用いる方が、スペクトルデータＦ（ｆ，Ｌ）を用いるよりも、近似の精度が向上することとなる。

　変化率算出部３５４は、関数ψ（ｆ，Ｌ）の周波数ｆに対する偏導関数∂ψ（ｆ，Ｌ）／∂ｆの近似値として、回帰分析により関数ψ（ｆ，Ｌ）の周波数ｆに対する変化率（すなわち、回帰直線の傾き）を算出する。続いて、変化率算出部３５４は、関数ψ（ｆ，Ｌ）の周波数ｆに対する変化率に対し、さらに回帰分析（第２の回帰分析）を行うことによって距離Ｌに対する変化率（すなわち、第２の回帰直線の傾き）を算出し、この値を２次偏導関数∂²ψ／∂Ｌ∂ｆの近似値とする。以下、この２次偏導関数∂²ψ／∂Ｌ∂ｆの近似値を２次変化率という。

　変化率算出部３５４のさらに具体的な処理を説明する。まず、変化率算出部３５４は、帯域設定部３５３が設定した周波数帯域で、関数ψ（ｆ，Ｌ）の周波数ｆを変数とする回帰直線の傾きおよび切片を得る。

　図１１は、例として距離Ｌ＝Ｌｐにおける関数ψ（ｆ，Ｌｐ）と周波数ｆの関係、および各関数の周波数ｆに対する回帰直線を示す図である（ｐ＝１，２，３，４）。回帰直線Ｊ１、Ｊ２、Ｊ３は、変化率算出部３５４が周波数帯域Ｕで回帰分析を行うことによって算出した回帰直線である。また、回帰直線Ｊ４’は、変化率算出部３５４が周波数帯域Ｕ’で回帰分析を行うことによって算出した回帰直線である。なお、図１１では、比較のため、変化率算出部３５４が、距離Ｌ＝Ｌ４において周波数帯域Ｕで回帰分析を行うことによって算出した回帰直線Ｊ４も記載している。以下、回帰直線Ｊｐの傾きをＳｆ（Ｌｐ）とする。また、回帰直線Ｊ４’の傾きをＳｆ’（Ｌ４）とする。

　回帰直線Ｊｐの傾きＳｆ（Ｌｐ）は、ｐ＝１，２，３において、距離Ｌの増加とともに単調に減少する。これに対し、回帰直線Ｊ４の傾きＳｆ（Ｌ４）は、回帰直線Ｊ３の傾きＳｆ（Ｌ３）よりも大きい（Ｓｆ（Ｌ４）＞Ｓｆ（Ｌ３））。一方、距離Ｌ＝Ｌ４の周波数帯域Ｕ’における傾きＳｆ’（Ｌ４）は、回帰直線Ｊ３の傾きＳｆ（Ｌ３）よりも小さい（Ｓｆ’（Ｌ４）＜Ｓｆ（Ｌ３））。これは、図８から明らかなように、関数ψ（ｆ，Ｌ４）が周波数ｆmax’の近傍で極小値を取り、周波数ｆmaxでは極小値より大きい値を取るからである。

　続いて、変化率算出部３５４は、傾きＳｆ（Ｌ）を距離Ｌの関数とみて第２の回帰分析を行うことにより、距離Ｌに対する第２の回帰直線を算出する。図１２は、傾きＳｆ（Ｌ）と距離Ｌの関係を示す図である。変化率算出部３５４は、上述の通り、例に挙げた傾きＳｆ（Ｌ１）、Ｓｆ（Ｌ２）、Ｓｆ（Ｌ３）、Ｓｆ’（Ｌ４）とそれ以外の全ての傾きを算出した。そして、変化率算出部３５４は、これらの傾きに基づき、さらに往復距離Ｌに対する第２の回帰分析を行うことで、第２の回帰直線Ｑ₁を算出する。この第２の回帰直線Ｑ₁を図１２に実線で示す。ところで、変化率算出部３５４は、上述の通り、周波数帯域Ｕにおいて、傾きＳｆ（Ｌ１）、Ｓｆ（Ｌ２）、Ｓｆ（Ｌ３）、Ｓｆ’（Ｌ４）とそれ以外の全ての傾きを算出した。比較のために、これらの傾きに基づき、さらに往復距離Ｌに対する第２の回帰分析を行うことによって算出した第２の回帰直線Ｑ₂を図１２に破線で示す。第２の回帰直線Ｑ₁、Ｑ₂を比較すると、第２の回帰直線Ｑ₁の方が、傾きＳｆ（Ｌ）の値によくフィットしていることがわかる。式（１５）に示した通り、関数ψ（ｆ，Ｌ）は周波数ｆの一次関数であり、かつ、本実施の形態１では周波数帯域Ｕ、Ｕ’を適切に設定してノイズの影響を排除したため、回帰直線Ｊ１、Ｊ２、Ｊ３、Ｊ４’の傾き（すなわち１次変化率）Ｓｆ（Ｌ１）、Ｓｆ（Ｌ２）、Ｓｆ（Ｌ３）、Ｓｆ’（Ｌ４）は、周波数ｆに対する関数ψ（ｆ，Ｌ）の１次偏導関数∂ψ（ｆ，Ｌ１）／∂ｆ、∂ψ（ｆ，Ｌ２）／∂ｆ、∂ψ（ｆ，Ｌ３）／∂ｆ、∂ψ（ｆ，Ｌ４）／∂ｆにそれぞれ良い近似を与える。さらに、式（１５）に示す通り、関数ψ（ｆ，Ｌ）は周波数ｆだけでなく往復距離Ｌの一次関数でもあるため、第２の回帰直線Ｑ₁の傾き（すなわち２次変化率）は、１次偏導関数∂ψ（ｆ，Ｌ）／∂ｆの往復距離Ｌに対する偏導関数、すなわち式（１３）の左辺の２次偏導関数∂²ψ（ｆ，Ｌ）／∂Ｌ∂ｆに良い近似を与える。

　減衰率推定部３５５は、変化率算出部３５４が算出した２次変化率の値を式（１３）へ代入することにより、超音波の音線ごとの減衰率ζを算出する。続いて減衰率推定部３５５は、算出した全ての音線の減衰率ζの平均値を算出し、算出結果を走査領域の減衰率として、画像処理部３６が有する合成画像データ生成部３６２へ出力する。なお、減衰率推定部３５５は、全ての音線の減衰率ζの最頻値、中央値または最大値等の統計量を走査領域の減衰率としてもよい。

　画像処理部３６は、エコー信号の振幅を輝度に変換して表示する超音波画像であるＢモード画像データを生成するＢモード画像データ生成部３６１と、減衰率推定部３５５が推定した減衰率ζの情報とＢモード画像データを合成して合成画像データを生成する合成画像データ生成部３６２と、を有する。

　Ｂモード画像データ生成部３６１は、信号処理部３４からのＢモード用受信データに対してゲイン処理、コントラスト処理等の公知の技術を用いた信号処理を行うとともに、表示装置４における画像の表示レンジに応じて定まるデータステップ幅に応じたデータの間引き等を行うことによってＢモード画像データを生成する。Ｂモード画像は、色空間としてＲＧＢ表色系を採用した場合の変数であるＲ（赤）、Ｇ（緑）、Ｂ（青）の値を一致させたグレースケール画像である。

　Ｂモード画像データ生成部３６１は、Ｂモード用受信データに走査領域を空間的に正しく表現できるよう並べ直す座標変換を施し、さらにＢモード用受信データ間の補間処理を施して、Ｂモード用受信データ間の空隙を埋め、デジタルデータであるＢモード画像データを生成する。Ｂモード画像データは図４で扇型に示された走査領域内の器官の状況を表現できるＢモード画像のデジタルデータである。Ｂモード画像データ生成部３６１は、生成したＢモード画像データを合成画像データ生成部３６２へ出力する。

　合成画像データ生成部３６２は、走査領域の減衰率の値を示す文字データを生成し、Ｂモード画像に隣接して表示するようにＢモード画像データと文字データとを合成することによって合成画像データを生成する。合成画像データ生成部３６２は、生成した合成画像データを表示装置４へ出力する。

　表示装置４は、液晶または有機ＥＬ（Electro　Luminescence）等からなるモニタを用いて構成される。表示装置４は、超音波診断装置３が生成した合成画像データに対応する合成画像を含む各種情報を表示する。

　図１３は、以上の構成を有する超音波診断装置３が実行する処理の概要を示すフローチャートである。具体的には、超音波診断装置３が超音波内視鏡２からエコー信号を受信する以降の処理の概要を示すフローチャートである。以下、図１３を参照して、超音波診断装置３が行う処理を説明する。まず、超音波診断装置３は、超音波内視鏡２から超音波振動子２１による測定結果としてのエコー信号を受信する（ステップＳ１）。

　超音波振動子２１からエコー信号を受信した信号増幅部３３１は、そのエコー信号の増幅を行う（ステップＳ２）。ここで、信号増幅部３３１は、例えば図３に示す増幅率と受信深度との関係に基づいてエコー信号の増幅（ＳＴＣ補正）を行う。この際、信号増幅部３３１におけるエコー信号の各種処理周波数帯域は、超音波振動子２１による超音波エコーのエコー信号への音響電気変換の線形応答周波数帯域をほぼカバーする広帯域にするとよい。これも、後述する周波数スペクトルの近似処理において精度のよい近似を行うことを可能とするためである。

　ところで、上述の通り、信号増幅部３３１はエコー信号を増幅し、送受信部３３は増幅されたエコー信号にフィルタリング、Ａ／Ｄ変換を施してＲＦデータを生成し、信号処理部３４はＲＦデータに各種処理を施してＢモード用受信データを生成した。ここで、Ｂモード画像データ生成部３６１は、信号処理部３４から入力されるＢモード用受信データに適切な座標変換、補間処理を施し、Ｂモード画像データを生成して、表示装置４へ出力する（ステップＳ３）。Ｂモード画像データを受信した表示装置４は、そのＢモード画像データに対応するＢモード画像を表示する。

　一方、増幅補正部３５１は、送受信部３３から出力されたＲＦデータに対して受信深度によらず増幅率が一定となるように増幅補正を行う（ステップＳ４）。ここで、増幅補正部３５１は、例えば図５に示す増幅率と受信深度との関係が成立するように増幅補正を行う。

　この後、周波数解析部３５２は、増幅補正後の各音線のＲＦデータに対してＦＦＴによる周波数解析を行うことによってスペクトルデータを算出する（ステップＳ５）。このステップＳ５の処理の詳細については後述する。

　続いて、帯域設定部３５３は、回帰分析対象の周波数帯域を設定する（ステップＳ６）。例えば、図８および図９に示す場合、Ｌ＝Ｌ１，Ｌ２，Ｌ３に対して周波数帯域Ｕを設定する一方、Ｌ＝Ｌ４に対して周波数帯域Ｕ’を設定する。

　変化率算出部３５４は、帯域設定部３５３が設定した周波数帯域に基づいて、関数ψ（ｆ，Ｌ）の２次偏導関数∂²ψ（ｆ，Ｌ）／∂Ｌ∂ｆの近似値である２次変化率を、２回の回帰分析を行うことによって算出する（ステップＳ７）。例えば、変化率算出部３５４は、２次変化率として、図１２に示す回帰直線Ｑ₁の傾きを算出する。

　この後、減衰率推定部３５５は、走査領域における超音波パルスの減衰率を推定する（ステップＳ８）。減衰率推定部３５５は、変化率算出部３５４が算出した２次変化率の値を式（１３）の左辺へ代入することにより、音線ごとの減衰率を算出した後、算出した全ての音線の減衰率の平均値を算出し、この平均値を走査領域の減衰率として合成画像データ生成部３６２へ出力する。

　合成画像データ生成部３６２は、Ｂモード画像データおよび走査領域の減衰率をもとに合成画像データを生成し、表示装置４へ出力する（ステップＳ９）。合成画像データを受信した表示装置４は、その合成画像データに対応する合成画像を表示する。図１４は、表示装置４が表示する合成画像の表示例を示す図である。同図に示す合成画像１０１は、Ｂモード画像表示部１０２と減衰率表示部１０３とを有する。なお、図１４では、具体的なＢモード画像の表示を省略している。

　ステップＳ９の後、超音波診断装置３は一連の処理を終了する。なお、超音波診断装置３は、ステップＳ１～Ｓ９の処理を周期的に繰り返し実行する。

　次に、ステップＳ５の周波数解析処理について、図１５のフローチャートを参照して説明する。まず、周波数解析部３５２は、解析対象の音線を識別するカウンタｋをｋ₀とする（ステップＳ１１）。この初期値ｋ₀は、術者等のユーザが入力部３７を通じて任意に指示入力した値、もしくは、記憶部３２にあらかじめ設定された値とする。

　続いて、周波数解析部３５２は、上述の通りＦＦＴ演算用に生成した一連のデータ群（サンプルデータ群）を代表するデータ位置（受信深度に相当）Ｚ^(k)の初期値Ｚ^(k) ₀を設定する（ステップＳ１２）。例えば、図６では、上述したように、音線ＳＲ_kの１番目のデータ位置を初期値Ｚ^(k) ₀として設定した場合を示している。

　その後、周波数解析部３５２は、サンプルデータ群を取得し（ステップＳ１３）、取得したサンプルデータ群に対し、記憶部３２が記憶する窓関数を作用させる（ステップＳ１４）。このようにサンプルデータ群に対して窓関数を作用させることにより、サンプルデータ群が境界で不連続になることを回避し、アーチファクトが発生するのを防止することができる。ステップＳ１４で適用する窓関数は、例えばHamming、Hanning、Blackman等のいずれかであり、予め記憶部３２に記憶されている。

　続いて、周波数解析部３５２は、データ位置Ｚ^(k)のサンプルデータ群が正常なデータ群であるか否かを判定する（ステップＳ１５）。図６を参照した際に説明したように、サンプルデータ群は、２のべき乗のデータ数を有している必要がある。以下、サンプルデータ群のデータ数を２ⁿ（ｎは正の整数）とする。本実施の形態１では、データ位置Ｚ^(k)が、できるだけＺ^(k)が属するサンプルデータ群の中心になるよう設定される。具体的には、サンプルデータ群のデータ数は２ⁿであるので、Ｚ^(k)はそのサンプルデータ群の中心に近い２ⁿ／２（＝２^n-1）番目の位置に設定される。この場合、サンプルデータ群が正常であるとは、データ位置Ｚ^(k)より浅部側に２^n-1－１（＝Ｎとする）個のデータがあり、データ位置Ｚ^(k)より深部側に２^n-1（＝Ｍとする）個のデータがあることを意味する。図６に示す場合、サンプルデータ群Ｆ_j（ｊ＝１,２，・・・，Ｋ－１）は正常である。なお、図６ではｎ＝４（Ｎ＝７，Ｍ＝８）の場合を例示している。

　ステップＳ１５における判定の結果、データ位置Ｚ^(k)のサンプルデータ群が正常である場合（ステップＳ１５：Ｙｅｓ）、周波数解析部３５２は、後述するステップＳ１７へ移行する。

　ステップＳ１５における判定の結果、データ位置Ｚ^(k)のサンプルデータ群が正常でない場合（ステップＳ１５：Ｎｏ）、周波数解析部３５２は、不足分だけゼロデータを挿入することによって正常なサンプルデータ群を生成する（ステップＳ１６）。ステップＳ１５において正常でないと判定されたサンプルデータ群（例えば図６のサンプルデータ群Ｆ_K）は、ゼロデータを追加する前に窓関数が作用されている。このため、サンプルデータ群にゼロデータを挿入してもデータの不連続は生じない。ステップＳ１６の後、周波数解析部３５２は、後述するステップＳ１７へ移行する。

　ステップＳ１７において、周波数解析部３５２は、サンプルデータ群を用いてＦＦＴ演算を行うことにより、振幅の周波数分布であるスペクトルデータを得る（ステップＳ１７）。この結果、例えば図７の各列に示すようなスペクトルデータが得られる。

　続いて、周波数解析部３５２は、データ位置Ｚ^(k)をステップ幅Ｄで変化させる（ステップＳ１８）。ステップ幅Ｄは、記憶部３２が予め記憶しているものとする。図６では、Ｄ＝１５の場合を例示している。ここで、図７に示した往復距離Ｌの間隔ΔＬは、このステップ幅Ｄを距離に換算したときの値（＝サンプリング幅×Ｄ）の２倍で定義される。したがって、ステップ幅Ｄが定まると、間隔ΔＬも一意に定まる。ステップ幅Ｄは、Ｂモード画像データ生成部３６１がＢモード画像データを生成する際に利用するデータステップ幅と一致させることが望ましいが、周波数解析部３５２における演算量を削減したい場合には、ステップ幅Ｄとしてデータステップ幅より大きい値を設定してもよい。

　その後、周波数解析部３５２は、データ位置Ｚ^(k)が音線ＳＲ_kにおける最大値Ｚ^(k) _maxより大きいか否かを判定する（ステップＳ１９）。データ位置Ｚ^(k)が最大値Ｚ^(k) _maxより大きい場合（ステップＳ１９：Ｙｅｓ）、周波数解析部３５２はカウンタｋを１増加させる（ステップＳ２０）。これは、処理をとなりの音線へ移すことを意味する。一方、データ位置Ｚ^(k)が最大値Ｚ^(k) _max以下である場合（ステップＳ１９：Ｎｏ）、周波数解析部３５２はステップＳ１３へ戻る。

　ステップＳ２０の後、周波数解析部３５２は、カウンタｋが最大値ｋ_maxより大きいか否かを判定する（ステップＳ２１）。カウンタｋがｋ_maxより大きい場合（ステップＳ２１：Ｙｅｓ）、周波数解析部３５２は一連の周波数解析処理を終了する。一方、カウンタｋがｋ_max以下である場合（ステップＳ２１：Ｎｏ）、周波数解析部３５２はステップＳ１２に戻る。この最大値ｋ_maxは、術者等のユーザが入力部３７を通じて任意に指示入力した値、もしくは、記憶部３２にあらかじめ設定された値とする。

　このようにして、周波数解析部３５２は、関心領域内の（ｋ_max－ｋ₀＋１）本の音線の各々について複数回のＦＦＴ演算を行う。

　以上説明した本発明の実施の形態１によれば、周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率をこの順序で算出することによって得られる２次変化率を用いて、超音波振動子の走査領域における超音波信号の単位距離および単位周波数あたりの減衰率を推定するため、超音波の減衰率を精度よくかつ簡便に算出することができるとともに、減衰率に基づく画像の信頼性を向上させることができる。

　また、本実施の形態１によれば、上述した引用文献１のように音速を得る必要はなく、送信波形がガウシアンであることを前提としていないため、減衰率を精度よく算出することができる。なお、本実施の形態１の２次変化率は、周波数と距離（または受信深度）の関数の２次変化率であり、上述した特許文献１における「位相の２次変化率」とは全く異なるものである。この点は、以下に説明する実施の形態２、３も同様である。

　また、本実施の形態１によれば、スペクトルデータに応じてノイズレベルデータとの比較により演算対象の周波数帯域を設定するため、Ｓ／Ｎが十分高く、減衰率の推定に有効な領域のみで減衰率を算出することができる。したがって、高精度で減衰率を算出することができ、減衰率に基づく画像の信頼性を高めることができる。

　また、本実施の形態１によれば、画像内の領域をマニュアルで指定するなどの煩雑な処理は不要であるため、減衰率を簡便に算出することができる。

　また、本実施の形態１によれば、線型の回帰分析を用いるため、周波数変化率および距離変化率を簡便に算出することができる。

　また、本実施の形態１によれば、周波数の１次関数であるスペクトルデータの関数ψ（ｆ，Ｌ）を用いているため、回帰直線の傾きを用いた近似の精度を向上させることができる。

　また、本実施の形態１によれば、ノイズレベルデータ記憶部は複数の種別または機体ごとの超音波振動子にそれぞれ対応するノイズレベルデータを記憶しているため、接続可能な全ての超音波振動子に対して超音波の減衰率を精度よく算出することができる。

　なお、本実施の形態１において、比較部３５３ａがスペクトルデータＦ（ｆ，Ｌ）と比較する対象をノイズレベルデータｎ（ｆ，Ｌ）そのものではなく、ノイズレベルデータｎ（ｆ，Ｌ）の１次関数ａ・ｎ（ｆ，Ｌ）＋ｂ（ａ≧１、ｂ≧０；ａ、ｂは定数）に代えてマージンを設けてもよい。この場合には、変化率算出部３５４が算出する２次変化率や、減衰率推定部３５５が推定する減衰率へのノイズの影響がさらに少なくなり、減衰率の算出精度が一層向上する。

　また、一般にＳ／Ｎが悪いのは、周波数に依存する減衰の激しい高周波数側であることに鑑み、帯域設定部３５３が回帰分析対象の周波数帯域を設定する際、最小の周波数を初期値ｆminで固定する一方、最大の周波数の値をｆ＞ｆminの範囲で初期値ｆmaxから変更するようにしてもよい。この場合には、変化率算出部３５４および減衰率推定部３５５が演算する際、減衰の激しい高周波数側を除外して減衰率の推定精度を高めることができる。

　また、超音波診断装置３は、走査領域の減衰率を推定する代わりに、走査領域の一部の領域の減衰率を推定するようにしてもよい。この場合の領域は、入力部３７を介してユーザが設定できるようしておけば好ましい。

（実施の形態２）
　本発明の実施の形態２は、超音波診断装置の変化率算出部が行う２次変化率の算出方法が、上述した実施の形態１と異なる。本実施の形態２に係る超音波診断装置は、実施の形態１で説明した超音波診断装置３と同様の構成を有する。

　本実施の形態２において、変化率算出部３５４は、回帰分析によりスペクトルデータＦ（ｆ，Ｌ）の距離Ｌに対する変化率をまず算出して偏導関数∂ψ（ｆ，Ｌ）／∂Ｌの近似値とした後、第２の回帰分析により偏導関数∂Ｆ（ｆ，Ｌ）／∂Ｌの周波数ｆに対する変化率を算出することによって２次偏導関数∂²Ｆ（ｆ，Ｌ）／∂ｆ∂Ｌの近似値すなわち２次変化率を算出する。

　この場合、距離Ｌ＝０のスペクトルデータＦ（ｆ，０）は距離Ｌの関数ではないため、
　　∂Ｆ（ｆ，Ｌ）／∂Ｌ＝∂ψ（ｆ，Ｌ）／∂Ｌ
が成り立つ。したがって、本実施の形態２において、変化率算出部３５４は関数ψ（ｆ，Ｌ）を算出する必要はなく、スペクトルデータＦ（ｆ，Ｌ）から２次変化率を算出することができる。

　変化率算出部３５４は、周波数ｆごとに回帰分析限界Ｌmax（ｆ）を抽出する。回帰分析限界Ｌmax（ｆ）は、周波数ｆにおける距離Ｌの最大値である。図１０に示すスペクトルデータのデータ列では、周波数ｆ＝ｆｐ（ｐ＝１，２，３，４）に対する回帰分析限界Ｌmax（ｆｐ）を示している。図１０に示す場合、回帰分析限界Ｌmax（ｆｐ）は、周波数ｆｐにおける白色セルの分布範囲の右端のセルが有する距離Ｌの値に対応している。

　変化率算出部３５４は、スペクトルデータＦ（ｆ，Ｌ）の距離Ｌに対する回帰分析を実行し、距離区間０≦Ｌ≦Ｌmax（ｆ）における回帰直線の傾きＳＬ（ｆ）を算出する。図１６は、周波数ｆを一定にしたときのスペクトルデータＦ（ｆ，Ｌ）と距離Ｌの関係を示す図である。具体的には、互いに異なる４つの距離におけるスペクトルデータＦ（ｆ，Ｌ１）、Ｆ（ｆ，Ｌ２）、Ｆ（ｆ，Ｌ３）、Ｆ（ｆ，Ｌ４）と距離Ｌとの関係を示している。ここで、周波数ｆｐは定数であり、０＜ｆ１＜ｆ２＜ｆ３＜ｆ４を満たす。図１６に示すように、スペクトルデータＦ（ｆ，Ｌ）は、被検体内で超音波が伝搬する際の周波数依存減衰の効果により、周波数ｆが大きいほど距離Ｌの増加に伴う減衰が激しい。なお、一般にこれより多くのスペクトルデータＦ（ｆ，Ｌ）が算出されるが、図１６では図８と同様に、代表的な４つのスペクトルデータのみを例示している。

　図１６では、距離区間０≦Ｌ≦Ｌmax（ｆｐ）におけるスペクトルデータＦ（ｆｐ，Ｌ）の回帰直線Ｋｐが示されている。スペクトルデータＦ（ｆｐ，Ｌ）は、Ｌ＝０から回帰分析限界Ｌmax（ｆｐ）に達するまではほぼ線型であるが、回帰分析限界Ｌmax（ｆｐ）より大きくなると、ノイズの影響およびスペクトルデータＦ（ｆｐ，Ｌ）自体が０の近傍まで減衰することによる影響によって線型性が崩れている。

　回帰直線Ｋｐの傾きは、周波数依存減衰の効果により、周波数ｆが大きくなるほど急峻になる。換言すると、回帰直線Ｋｐの傾きをＳＬ（ｆｐ）とすると、ＳＬ（ｆ１）＞ＳＬ（ｆ２）＞ＳＬ（ｆ３）＞ＳＬ（ｆ４）が成立する。

　ここで比較のため、図１６に示す４つのスペクトルデータＦ（ｆｐ，Ｌ）について、４つの周波数のうち最低の周波数ｆ１の回帰分析限界Ｌmax（ｆ１）を用いて回帰直線を算出する場合を考える。図１６では、一例として、スペクトルデータＦ（ｆ４，Ｌ）に対して回帰分析限界Ｌmax（ｆ１）を適用した場合に算出される回帰直線Ｋ４’を例示している。この回帰直線Ｋ４’の傾きＳＬ’（ｆ４）は、スペクトルデータＦ（ｆ４，Ｌ）がＬ＝Ｌmax（ｆ４）の近傍でノイズレベルまで落ち切ってしまうため、ＳＬ（ｆ４）よりも大きい（ＳＬ’（ｆ４）＞ＳＬ（ｆ４））。なお、図示はしないが、ｐ＝１，２，３の場合にも、同様の関係が成り立っている。

　変化率算出部３５４は、傾きＳＬ（ｆ）を周波数ｆの関数とみて第２の回帰分析を行うことにより、周波数ｆに対する第２の回帰直線を算出する。図１７は、傾きＳＬ（ｆ）と周波数ｆとの関係を示す図である。図１７において、黒点は、変化率算出部３５４が距離区間０≦Ｌ≦Ｌmax（ｆｐ）でスペクトルデータＦ（ｆｐ，Ｌ）の回帰分析を行うことによって算出した傾きＳＬ（ｆｐ）と周波数ｆｐとの関係を示している。これに対して白点は比較用のデータであり、変化率算出部３５４が距離区間０≦Ｌ≦Ｌmax（ｆ１）でスペクトルデータＦ（ｆｐ，Ｌ）の回帰分析を行うことによって算出した回帰直線Ｋｐ’の傾きＳＬ’（ｆｐ）と周波数ｆｐとの関係を示している。

　変化率算出部３５４は、上述の通り、距離区間０≦Ｌ≦Ｌmax（ｆ）において、例に挙げた傾きＳＬ（ｆ１）、ＳＬ（ｆ２）、ＳＬ（ｆ３）、ＳＬ（ｆ４）と、それ以外の全ての傾きを算出した。そして、変化率算出部３５４は、これらの傾きに基づき、さらに周波数ｆに対する第２の回帰分析を行うことで、第２の回帰直線Ｑ₃を算出する。この第２の回帰直線Ｑ₃を図１７に実線で示す。ところで、比較のために、変化率算出部３５４が、仮に距離区間０≦Ｌ≦Ｌmax（ｆ１）において、傾きＳＬ’（ｆ１）、ＳＬ’（ｆ２）、ＳＬ’（ｆ３）、ＳＬ’（ｆ４）とそれ以外をの全ての傾きを算出し、これらの傾きに基づき、さらに周波数ｆに対する第２の回帰分析を行った場合の第２の回帰直線Ｑ₄を図１７に破線で示す。第２の回帰直線Ｑ₃、Ｑ₄を比較すると、第２の回帰直線Ｑ₃の方が、傾きＳＬ（ｆ）の値によくフィットしていることがわかる。第２の回帰直線Ｑ₃の傾きは、式（１４）の左辺の２次偏導関数∂²ψ（ｆ，Ｌ）／∂ｆ∂Ｌに対する近似値、すなわち２次変化率を与える。

　以上説明した本発明の実施の形態２によれば、周波数スペクトルにおける周波数変化率および距離変化率をこの順序で算出することによって得られる２次変化率を用いて、超音波振動子の走査領域における超音波信号の単位距離および単位周波数あたりの減衰率を推定するため、実施の形態１と同様、超音波の減衰率を精度よくかつ簡便に算出することができるとともに、減衰率に基づく画像の信頼性を向上させることができる。

　また、本実施の形態２によれば、距離に対する回帰分析を周波数に対する回帰分析よりも先に実行するため、スペクトルデータをそのまま回帰分析することができ、計算量を少なくすることができる。

（実施の形態３）
　図１８は、本発明の実施の形態３に係る超音波診断システムの機能構成を示すブロック図である。同図に示す超音波診断システム１１は、超音波内視鏡２と、超音波診断装置１２と、表示装置４とを備える。

　超音波診断装置１２は、記憶部１２１および画像処理部１２２の構成が、上述した超音波診断装置３の記憶部３２および画像処理部３６とそれぞれ異なる。

　記憶部１２１は、種別情報記憶部３２１、ノイズレベルデータ記憶部３２２、帯域情報記憶部３２３に加えて、減衰率の値に応じて画像に付与する視覚情報を記憶する視覚情報記憶部３２４を有する。ここでいう視覚情報とは、例えば輝度、色相、明度または彩度等のいずれかであり、減衰率の値に応じた値が対応付けられている。なお、視覚情報記憶部３２４が複数種類の視覚情報を減衰率と対応づけて記憶しておいてもよい。この場合には、ユーザが入力部３７を介して所望の視覚情報を選択できるようにしておけばよい。

　画像処理部１２２は、Ｂモード画像データ生成部３６１、合成画像データ生成部３６２に加えて、減衰率マップデータ生成部３６３を有する。減衰率マップデータ生成部３６３は、減衰率推定部３５５が推定した減衰率の値に応じた視覚情報を、視覚情報記憶部３２４を参照して画像に付与することによって減衰率マップデータを生成する。

　図１９は、以上の構成を有する超音波診断装置１２が行う処理の概要を示すフローチャートである。具体的には、超音波診断装置１２が超音波内視鏡２からエコー信号を受信する以降の処理の概要を示すフローチャートである。図１９において、ステップＳ３１～ステップＳ３６は、図１３のステップＳ１～Ｓ６に順次対応している。以下、ステップＳ３６に続く処理を説明する。

　ステップＳ３６の後、変化率算出部３５４は、スペクトルデータの２次変化率を算出する（ステップＳ３７）。この際、変化率算出部３５４は、走査領域内であらかじめ分割して設定された複数の部分領域の各々に対して２次変化率を算出する。なお、本実施の形態３において、変化率算出部３５４が２次変化率を算出する際の周波数変化率および距離変化率の算出順序は、特に限定されない。

　図２０は、部分領域の設定例を模式的に示す図である。同図に示す部分領域Ｒは、深度方向の長さ（深度幅）がｗで、ＮＲ本の音線を含む扇形領域である。走査領域Ｓでは、深度幅ｗで送受信方向を分割するとともに、ＮＲ本の音線ごとに走査方向が分割されている。図２０では、走査領域Ｓ内に複数ある送受信方向のうち、１本の送受信方向ｄｒ上で算出されたスペクトルデータのうち部分領域Ｒ内のデータを黒丸で示す一方、部分領域Ｒの境界に位置するスペクトルデータを白丸で示している。

　変化率算出部３５４は、部分領域Ｒにおけるスペクトルデータにおける距離の基準位置（距離がゼロの位置）を、部分領域Ｒ内で超音波振動子２１の表面に最も近い位置とし、この基準位置との深度の差ｚ’の２倍（２ｚ’）を新たな距離として演算を行う。なお、超音波振動子２１の表面を距離の基準位置としたときの部分領域Ｒの基準位置がＬ＝Ｌminである場合、上述した白丸の位置におけるスペクトルデータはＦ（ｆ，Ｌmin）と表される。変化率算出部３５４は、このスペクトルデータＦ（ｆ，Ｌmin）を式（１０）のスペクトルデータＦ（ｆ，０）の代わりとして、部分領域Ｒにおける２次変化率を算出するようにしてもよい。

　減衰率推定部３５５は、変化率算出部３５４の算出結果を用いて、走査領域に含まれる複数の部分領域の各々に対し、部分領域の減衰率を推定する（ステップＳ３８）。減衰率推定部３５５は、まず部分領域で算出された全ての２次変化率を用いて音線ごとの減衰率を算出する。この後、減衰率推定部３５５は、同一の部分領域で算出した全ての音線の減衰率の平均値を算出し、この平均値を当該部分領域の減衰率の推定値として減衰率マップデータ生成部３６３へ出力する。なお、減衰率推定部３５５は、同一の部分領域で算出した全ての音線の減衰率の最頻値、中央値または最大値等の統計量を当該部分領域の減衰率としてもよい。

　この後、減衰率マップデータ生成部３６３は、視覚情報記憶部３２４を参照して、各部分領域の減衰率に応じた視覚情報を各部分領域に対して割り当てることによって減衰率マップデータを生成し、合成画像データ生成部３６２へ出力する（ステップＳ３９）。なお、合成画像データ生成部３６２が、各部分領域の減衰率の推定値を文字情報としてさらに表示する合成画像データを生成するようにしてもよい。

　合成画像データ生成部３６２は、Ｂモード画像データに減衰率マップデータを重畳することによって合成画像データを生成し、この合成画像データを表示装置４へ出力する（ステップＳ４０）。合成画像データを受信した表示装置４は、その合成画像データに対応する合成画像を表示する。図２１は、表示装置４が表示する減衰率マップデータ付きの合成画像の表示例を示す図である。図２１に示す合成画像２０１は、各領域ごとに異なる視覚情報が割り当てられている。なお、図２１では、視覚情報を模様で模式的に記載している。また、図２１では、簡単のためＢモード画像の具体的な表示を省略している。

　本実施の形態３において、深度幅ｗは例えば１ｃｍ程度であるのが好ましい。深度幅ｗが１ｃｍ程度である場合、周波数解析部３５２がＦＦＴ処理を行う際のＲＦデータを切り取る区間の幅ΔＬを２ｍｍ程度とするのが望ましい。このとき、間隔ΔＬに対応する基準位置の深度ｚ’の幅Δｚ’（＝ΔＬ／２）は、１ｍｍ程度になる。

　以上説明した本発明の実施の形態３によれば、周波数スペクトルまたは周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率を所定の順序で算出することによって得られる２次変化率を用いて、超音波振動子の走査領域における超音波信号の単位距離および単位周波数あたりの減衰率を推定するため、実施の形態１、２と同様、超音波の減衰率を精度よくかつ簡便に算出することができるとともに、減衰率に基づく画像の信頼性を向上させることができる。

　また、本実施の形態３によれば、超音波振動子の走査領域を分割することによって得られる複数の部分領域の各々に対し、各部分領域で得られた減衰率の統計量を算出することによって各部分領域の減衰率を推定した後、各部分領域の減衰率の値に応じた視覚情報を付与することによって減衰率マップデータを生成するため、ユーザが減衰率の分布を視覚的に把握しやすい情報を提供することができる。

　また、本実施の形態３において、減衰率マップの最小単位である扇形領域のうち隣接する扇形を重複させるように設定してもよい。

　また、本実施の形態３において、各扇形領域の減衰率の算出を順次行う代わりに、同時並行的に行ってもよい。

　また、本実施の形態３において、入力部３７が入力を受け付ける設定信号に基づいて１つの扇形領域を関心領域（ＲＯＩ：Region　of　Interest）として設定し、関心領域内の減衰率の値をさらに合成することによって合成画像データを生成してもよい。

（その他の実施の形態）
　ここまで、本発明を実施するための形態を説明してきたが、本発明は、上述した実施の形態１～３によってのみ限定されるべきものではない。例えば、超音波診断装置において、各機能を有する回路同士をバスで接続することによって構成してもよいし、一部の機能が他の機能の回路構造に内蔵されるように構成してもよい。具体的には、例えば減衰率推定部の機能を有する回路に変化率算出部の機能を内蔵してもよい。

　また、超音波プローブとして、光学系のない細径の超音波ミニチュアプローブを適用してもよい。超音波ミニチュアプローブは、通常、胆道、胆管、膵管、気管、気管支、尿道、尿管へ挿入され、その周囲臓器（膵臓、肺、前立腺、膀胱、リンパ節等）を観察する際に用いられる。

　また、超音波プローブとして、被検体の体表から超音波を照射する体外式超音波プローブを適用してもよい。体外式超音波プローブは、通常、腹部臓器（肝臓、胆嚢、膀胱）、乳房（特に乳腺）、甲状腺を観察する際に用いられる。

　また、超音波振動子は、リニア振動子でもラジアル振動子でもコンベックス振動子でも構わない。超音波振動子がリニア振動子である場合、その走査領域は矩形（長方形、正方形）をなし、超音波振動子がラジアル振動子やコンベックス振動子である場合、その走査領域は扇形や円環状をなす。

　このように、本発明は、特許請求の範囲に記載した技術的思想を逸脱しない範囲内において、様々な実施の形態等を含み得るものである。

　以上のように、本発明にかかる超音波診断装置、超音波診断装置の作動方法および超音波診断装置の作動プログラムは、超音波の減衰率を精度よくかつ簡便に算出するとともに、減衰率に基づく画像の信頼性を向上させるのに有用である。

　１、１１　超音波診断システム
　２　超音波内視鏡
　３、１２　超音波診断装置
　４　表示装置
　２１　超音波振動子
　３１　制御部
　３２、１２１　記憶部
　３３　送受信部
　３４　信号処理部
　３５　演算部
　３６、１２２　画像処理部
　３７　入力部
　１０１、２０１　合成画像
　１０２　Ｂモード画像表示部
　１０３　減衰率表示部
　３１１　種別特定部
　３２１　種別情報記憶部
　３２２　ノイズレベルデータ記憶部
　３２３　帯域情報記憶部
　３２４　視覚情報記憶部
　３３１　信号増幅部
　３５１　増幅補正部
　３５２　周波数解析部
　３５３　帯域設定部
　３５３ａ　比較部
　３５４　変化率算出部
　３５５　減衰率推定部
　３６１　Ｂモード画像データ生成部
　３６２　合成画像データ生成部
　３６３　減衰率マップデータ生成部

Claims (15)

1. 　被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波診断装置であって、
   　前記超音波信号の周波数を解析することによって前記超音波信号の受信深度ごとの周波数スペクトルを算出する周波数解析部と、
   　前記周波数解析部が算出した周波数スペクトルまたは該周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率を所定の順序で算出することによって前記周波数スペクトルまたは前記関数の２次変化率を算出する変化率算出部と、
   　前記変化率算出部が算出した２次変化率を用いて、前記超音波振動子の走査領域内の所定領域における前記超音波信号の単位距離および単位周波数あたりの減衰率を推定する減衰率推定部と、
   　を備えたことを特徴とする超音波診断装置。
2. 　周波数および受信深度に応じたノイズレベルデータを記憶するノイズレベルデータ記憶部と、
   　前記ノイズレベルデータ記憶部が記憶するノイズレベルデータと、前記周波数スペクトルのデータとを比較することにより、前記変化率算出部が算出対象とする前記周波数スペクトルの周波数帯域を設定する帯域設定部と、
   　をさらに備え、
   　前記変化率算出部は、前記帯域設定部が設定した周波数帯域における前記周波数変化率を算出することを特徴とする請求項１に記載の超音波診断装置。
3. 　前記帯域設定部は、同一周波数における前記周波数スペクトルのデータと前記ノイズレベルデータとを比較することによってノイズに相当する周波数の帯域であるノイズ周波数帯域を抽出し、該ノイズ周波数帯域を除外することによって前記周波数帯域を設定することを特徴とする請求項２に記載の超音波診断装置。
4. 　前記帯域設定部は、所定の周波数より高い周波数帯域を前記ノイズ周波数帯域として抽出することを特徴とする請求項３に記載の超音波診断装置。
5. 　前記帯域設定部は、前記周波数スペクトルのデータの値が閾値以下の領域を前記ノイズ周波数帯域として抽出することを特徴とする請求項３に記載の超音波診断装置。
6. 　前記閾値は、前記ノイズレベルデータ以上の値であることを特徴とする請求項５に記載の超音波診断装置。
7. 　前記変化率算出部は、前記超音波の音線ごとに前記２次変化率を算出し、
   　前記減衰率推定部は、前記音線ごとの前記２次変化率を用いて前記音線ごとの減衰率を算出し、前記走査領域における全ての前記音線ごとの減衰率の統計量を算出することによって前記減衰率を推定することを特徴とする請求項１に記載の超音波診断装置。
8. 　前記変化率算出部は、
   　前記周波数スペクトルと基準の受信深度における前記周波数スペクトルとの差によって定義される関数に対し、前記周波数帯域で周波数を変数とする回帰分析を行うことによって回帰直線を算出し、前記回帰直線の傾きに対して前記超音波振動子と前記被検体との往復距離を変数とする第２の回帰分析を行うことによって第２の回帰直線を算出し、該第２の回帰直線の傾きを前記２次変化率とすることを特徴とする請求項２に記載の超音波診断装置。
9. 　前記変化率算出部は、
   　前記周波数帯域に基づいて周波数ごとに定まる距離区間で前記周波数スペクトルに対して前記超音波振動子と前記被検体との往復距離を変数とする回帰分析を行うことによって回帰直線を算出し、該回帰直線の傾きに対して周波数を変数とする第２の回帰分析を行うことによって第２の回帰直線を算出し、該第２の回帰直線の傾きを前記２次変化率とすることを特徴とする請求項２に記載の超音波診断装置。
10. 　前記超音波画像として信号の振幅を輝度に変換して表示するＢモード画像データを生成するＢモード画像データ生成部と、
    　前記減衰率推定部が推定した減衰率に関する情報と前記Ｂモード画像データ生成部が生成したＢモード画像データとを用いて合成画像データを生成する合成画像データ生成部と、
    　をさらに備えたことを特徴とする請求項１に記載の超音波診断装置。
11. 　前記減衰率推定部は、前記走査領域を分割することによって得られる複数の部分領域の各々に対し、各部分領域で得られた減衰率の統計量を算出することによって各部分領域の減衰率を推定し、
    　各部分領域の減衰率の値に応じた視覚情報を付与することによって減衰率マップデータを生成する減衰率マップデータ生成部をさらに備えたことを特徴とする請求項１に記載の超音波診断装置。
12. 　前記超音波画像として信号の振幅を輝度に変換して表示するＢモード画像データを生成するＢモード画像データ生成部と、
    　前記減衰率マップデータ生成部が生成した減衰率マップデータと前記Ｂモード画像データ生成部が生成したＢモード画像データとを用いて合成画像データを生成する合成画像データ生成部と、
    　をさらに備えたことを特徴とする請求項１１に記載の超音波診断装置。
13. 　前記ノイズレベルデータ記憶部は、当該超音波診断装置に接続可能な前記超音波プローブが備える前記超音波振動子に応じたノイズレベルデータを記憶することを特徴とする請求項２に記載の超音波診断装置。
14. 　被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波診断装置の作動方法であって、
    　周波数解析部が前記超音波信号の周波数を解析することによって前記超音波信号の受信深度ごとの周波数スペクトルを算出する周波数解析ステップと、
    　変化率算出部が前記周波数スペクトルまたは該周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率を所定の順序で算出することによって前記周波数スペクトルまたは前記関数の２次変化率を算出する変化率算出ステップと、
    　減衰率推定部が前記２次変化率を用いて、前記超音波振動子の走査領域内の所定領域における前記超音波信号の単位距離および単位周波数あたりの減衰率を推定する減衰率推定ステップと、
    　を有することを特徴とする超音波診断装置の作動方法。
15. 　被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波診断装置に、
    　周波数解析部が前記超音波信号の周波数を解析することによって前記超音波信号の受信深度ごとの周波数スペクトルを算出する周波数解析ステップと、
    　変化率算出部が前記周波数スペクトルまたは該周波数スペクトルを用いて定義される関数における距離変化率および周波数変化率を所定の順序で算出することによって前記周波数スペクトルまたは前記関数の２次変化率を算出する変化率算出ステップと、
    　減衰率推定部が前記２次変化率を用いて、前記超音波振動子の走査領域内の所定領域における前記超音波信号の単位距離および単位周波数あたりの減衰率を推定する減衰率推定ステップと、
    　を実行させることを特徴とする超音波診断装置の作動プログラム。

Images