特開 2022-26116 超音波ＣＴ装置、超音波画像生成方法、および、超音波画像生成装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2022026116

(43)【公開日】2022-02-10

(54)【発明の名称】超音波ＣＴ装置、超音波画像生成方法、および、超音波画像生成装置

(51)【国際特許分類】

   A61B 8/15 20060101AFI20220203BHJP

【ＦＩ】

A61B8/15

【審査請求】未請求

【請求項の数】13

【出願形態】ＯＬ

(21)【出願番号】P 2020129431

(22)【出願日】2020-07-30

(71)【出願人】

【識別番号】320011683

【氏名又は名称】富士フイルムヘルスケア株式会社

(74)【代理人】

【識別番号】110000888

【氏名又は名称】特許業務法人 山王坂特許事務所

(72)【発明者】

【氏名】鈴木 敦郎

(72)【発明者】

【氏名】坪田 悠史

(72)【発明者】

【氏名】寺田 崇秀

(72)【発明者】

【氏名】山中 一宏

(72)【発明者】

【氏名】川畑 健一

【テーマコード（参考）】

4C601

【Ｆターム（参考）】

4C601BB03

4C601BB13

4C601DD08

4C601DD20

4C601DE18

4C601EE01

4C601EE04

4C601GB05

4C601GC02

4C601GC10

(57)【要約】

【課題】ＦＷＩ法においてシミュレーションにより模擬音源から発生する音波と、振動子が照射する音波との波形のずれを低減し、高解像度の再構成画像を得る。
【解決手段】複数の振動子が配列された振動子アレイの１以上の振動子から音波を被写体に対して送信する。被写体の撮像領域を透過した音波を複数の振動子が測定した測定音圧を複数の振動子から受け取る。測定音圧を処理して、撮像領域の透過波画像を生成する。模擬音源から時間的に音圧が変化する模擬音波を発生させ、模擬音波が撮像領域を透過して複数の模擬検出器に到達する際の音圧を計算音圧として、計算により求め、透過波画像を初期画像として、計算音圧を用いて透過波画像を逐次的に修正する。このとき、逐次更新部は、計算音圧を用いて、模擬音源の発生する模擬音波の波形を、振動子が送信する音波の波形に近づける演算処理を行う。
【選択図】図５

【特許請求の範囲】

【請求項1】

複数の振動子が配列された振動子アレイの１以上の前記振動子から音波を被写体に対して送信させる送信部と、
前記被写体の撮像領域を透過した前記音波の音圧を複数の振動子により測定した測定音圧を前記複数の振動子から受け取る受信部と、
前記測定音圧を処理して、前記撮像領域の透過波画像を生成する画像再構成部と、
模擬音源から時間的に音圧が変化する模擬音波を発生させ、当該模擬音波が前記撮像領域を透過して複数の模擬検出器に到達する際の音圧を計算音圧として、計算により求め、前記透過波画像を初期画像として、前記計算音圧を用いて前記透過波画像を逐次的に修正する逐次更新部とを有し、
前記逐次更新部は、前記計算音圧を用いて、前記模擬音源の発生する前記模擬音波の波形を、前記振動子が送信する前記音波の波形に近づける演算処理を行うことを特徴とする超音波ＣＴ装置。

【請求項2】

請求項１に記載の超音波ＣＴ装置において、前記逐次更新部は、前記演算処理として、前記模擬音源の発生する前記模擬音波の前記時間的な音圧の変化の波形を複数種類に異ならせ、当該複数種類の前記波形についてそれぞれ前記計算音圧を求め、前記測定音圧と前記計算音圧との差分から、前記透過波画像の逐次的な修正に用いる前記模擬音波の波形を、前記複数種類の波形の中から選択することを特徴とする超音波ＣＴ装置。

【請求項3】

請求項２に記載の超音波ＣＴ装置であって、前記模擬音波の前記複数種類の波形は、位相が異なることを特徴とする超音波ＣＴ装置。

【請求項4】

請求項２に記載の超音波ＣＴ装置であって、前記逐次更新部は、前記模擬音波の前記複数種類の波形について、前記測定音圧と前記計算音圧との差分を表す所定の関数の値を算出し、その関数の値が最小となる波形を選択することを特徴とする超音波ＣＴ装置。

【請求項5】

請求項４に記載の超音波ＣＴ装置であって、前記逐次更新部は、前記所定の関数としては、前記測定音圧と前記計算音圧との差の二乗和を算出することを特徴とする超音波ＣＴ装置。

【請求項6】

請求項２に記載の超音波ＣＴ装置であって、前記模擬音波の前記複数種類の波形の音圧は、a(k)f（ただし、fは、任意の複素数、a(k)は、a(k)=exp{i・θ_k}、θ_k=（2π/N）・k、k=0, 1, …N-1、Nは任意の整数）で表される音源から計算することを特徴とする超音波ＣＴ装置。

【請求項7】

請求項２に記載の超音波ＣＴ装置であって、前記逐次更新部は、前記透過波画像を逐次的に修正する際に、前記測定音圧と前記計算音圧との差分を最小にする音速の補正係数を求め、前記音速の補正係数を用いて前記透過波画像を逐次的に修正することを特徴とする超音波ＣＴ装置。

【請求項8】

請求項１に記載の超音波ＣＴ装置であって、前記逐次更新部は、前記模擬音源を、前記送信部が音波を送信させる前記振動子の位置に配置し、前記模擬検出器を前記受信部が前記測定音圧を受け取る前記振動子の位置に配置することを特徴とする超音波ＣＴ装置。

【請求項9】

請求項１に記載の超音波ＣＴ装置であって、前記逐次更新部は、前記被写体を配置せずに水について測定した測定音圧と、前記模擬音波の波形を前記振動子が送信する前記音波の波形に近づける前記演算処理後の前記模擬音波について求めた前記計算音圧との位相差および音圧差の少なくとも一方を算出することを特徴とする超音波ＣＴ装置。

【請求項10】

請求項９に記載の超音波ＣＴ装置であって、前記逐次更新部は、前記位相差および前記音圧差の少なくとも一方が予め定めた閾値以上の場合、キャリブレーションの実行を促す表示を表示部に表示させることを特徴とする超音波ＣＴ装置。

【請求項11】

請求項９に記載の超音波ＣＴ装置であって、前記逐次更新部は、前記位相差および前記音圧差の少なくとも一方が予め定めた閾値以上の場合、前記閾値以上の位相差および音圧差となった測定音圧をマスク処理して、前記透過波画像の逐次的な修正に用いないことを特徴とする超音波ＣＴ装置。

【請求項12】

複数の振動子が配列された振動子アレイの１以上の前記振動子から音波を被写体に対して送信し、
前記被写体の撮像領域を透過した前記音波の音圧を複数の振動子により測定した測定音圧を受け取り、
前記測定音圧を処理して、前記撮像領域の透過波画像を生成し、
模擬音源から時間的に音圧が変化する模擬音波を発生させ、当該模擬音波が前記撮像領域を透過して複数の模擬検出器に到達する際の音圧を計算音圧として、計算により求め、
前記計算音圧を用いて、前記模擬音源の発生する前記模擬音波の波形を、前記振動子が送信する前記音波の波形に近づけ、
前記音波の波形に近づけた前記模擬音波の前記計算音圧を用いて、前記透過波画像を初期画像として、前記透過波画像を逐次的に修正する
ことを特徴とする超音波画像生成方法。

【請求項13】

複数の振動子が配列された振動子アレイの１以上の前記振動子から音波を被写体に対して送信し、前記被写体の撮像領域を透過した前記音波の音圧を複数の振動子により測定した測定音圧を受け取って、前記測定音圧を処理して、前記撮像領域の透過波画像を生成する画像再構成部と、
模擬音源から時間的に音圧が変化する模擬音波を発生させ、当該模擬音波が前記撮像領域を透過して複数の模擬検出器に到達する際の音圧を計算音圧として、計算により求め、前記透過波画像を初期画像として、前記計算音圧を用いて前記透過波画像を逐次的に修正する逐次更新部とを有し、
前記逐次更新部は、前記計算音圧を用いて、前記模擬音源の発生する前記模擬音波の波形を、前記振動子が送信する前記音波の波形に近づける演算処理を行うことを特徴とする超音波画像生成装置。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、超音波ＣＴの画像再構成法に関するものである。

【背景技術】

【0002】

超音波測定を乳がんの検出に応用した医療用診断装置として、乳房専用超音波ＣＴ(Computed tomography)装置が知られている（特許文献１）。超音波ＣＴ装置では、水中に挿入された乳房の周囲に、超音波の送信および受信器である振動子アレイを配置し、乳房を透過した超音波を全周方向について受信し、受信信号から音速や減衰量の分布を示す画像を再構成する。一般に、組織の音速および減衰量は定量値であり、腫瘍における音速および減衰量は、周囲の乳腺、脂肪等の正常組織に比べて高いため、透過波画像から腫瘍を定量的に検出することが可能である。

【0003】

透過波画像の撮像方法として、一つの音源（振動子）から所定の角度で広がる拡散波を被写体に照射し、透過信号を取得する方法が知られている（非特許文献１）。

【0004】

透過信号から透過波画像を生成する画像再構成法としては、ストレートレイ(straight-ray)法やベントレイ（bent-ray法)が知られている。straight-ray法は超音波の軌跡を直線に近似して画像再構成を行う方法であり、計算は高速であるが空間分解能は低い。bent-ray法は超音波の屈折を考慮して画像再構成を行う方法であり、空間分解能はstraight-ray法よりも高い。

【0005】

特許文献１には、ストレートレイ法やベントレイ法よりも、高い空間分解能の音速画像を再構成する方法として、ＦＷＩ（full waveform inversion）法が開示されており、臨床データへも適用されている。

【0006】

ＦＷＩ法では、測定された透過信号から従来の再構成法で音速画像を得て、この音速画像を初期画像とし、ある音源と初期画像からシミュレーションにより音圧の分布を計算し、この作成したシミュレーションデータ（音圧の分布）と測定データとの誤差が小さくなるように音速画像を更新してゆく。特許文献１には、ＦＷＩ法の計算アルゴリズムが開示されている。

【先行技術文献】

【特許文献】

【0007】

【特許文献1】米国特許出願公開第２０１６／００３００００号明細書

【非特許文献1】Pratt, R., ”Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model.” GEOPHYSICS, 1999, 64(3), 888-901.

【発明の概要】

【発明が解決しようとする課題】

【0008】

ＦＷＩ法では、撮像領域の音圧分布のシミュレーションに用いる音源は、実際の音源（振動子）と同様に、時間的に音圧が変化する波を発生する。シミュレーションの音源の発する波の信号を、実際の音源の発する波の信号に一致させるため、特許文献１では、その段落００４５に記載されているように、信号源スケーリング係数γを算出している。具体的には、特許文献１では、信号源スケーリング係数γの推定を線形推定問題として扱い、複素数のスケーリング係数γを、実測値の受信信号とシミュレーション結果の波の信号とを掛け合わせたものを、シミュレーション結果の波の信号の二乗で除することにより算出することを開示している。

【0009】

しかしながら、上記スケーリング係数γの算出には、実際の音源が発生する波の位相と、シミュレーション時の音源が発する波の位相が一致しているかどうかについては考慮されていない。

【0010】

シミュレーション時の音源から発生する波の位相が、実際の音源が発生する波の位相と異なると、受信信号の波の到達時刻の推定に最大で１周期分のずれが生じる可能性がある。そのため、両者の位相を精度よく一致させることがFWIの画質向上のために重要である。

【0011】

しかしながら、実際の振動子から発せられる時点の波の位相を高精度に制御することは、現実的には難しいため、シミュレーションにおける音源が発する波の位相を、実際の音源の波の位相に一致するように精度よく推定する必要がある。

【0012】

本発明の目的は、ＦＷＩ法においてシミュレーションにより模擬音源から発生する音波と、振動子が照射する音波との波形のずれを低減し、高解像度の再構成画像を得ることにある。

【課題を解決するための手段】

【0013】

上記目的を達成するために、本発明によれば、複数の振動子が配列された振動子アレイの１以上の振動子から音波を被写体に対して送信させる送信部と、被写体の撮像領域を透過した音波を複数の振動子が測定した測定音圧を複数の振動子から受け取る受信部と、測定音圧を処理して、撮像領域の透過波画像を生成する画像再構成部と、模擬音源から時間的に音圧が変化する模擬音波を発生させ、模擬音波が撮像領域を透過して複数の模擬検出器に到達する際の音圧を計算音圧として、計算により求め、透過波画像を初期画像として、計算音圧を用いて透過波画像を逐次的に修正する逐次更新部とを有する。逐次更新部は、計算音圧を用いて、模擬音源の発生する模擬音波の波形を、振動子が送信する音波の波形に近づける演算処理を行う。

【発明の効果】

【0014】

本発明によれば、ＦＷＩ法においてシミュレーションにより模擬的な音源から発生する波の音圧の変化を、実際の振動子が照射する音圧の変化に近づけることができるため、高解像度の再構成画像を得ることができる。

【図面の簡単な説明】

【0015】

【図1】実施形態１の超音波ＣＴ装置の構成を示すブロック図である。

【図2】実施形態１の超音波ＣＴ装置の振動子アレイ３０の斜視図である。

【図3】振動子３ａから照射される照射波を示す説明図である。

【図4】（ａ）超音波ＣＴ装置で再構成された透過波画像（初期画像）を示す説明図であり、（ｂ）透過波画像に対して模擬音源５３を設定し模擬音波５０を照射して模擬検出器５４の位置で計算音圧を求めることを示す説明図であり、（ｃ）計算音圧と測定音圧の差を示すグラフであり、（ｄ）計算音圧を測定音圧の差から求めた補正係数Δｃの画像であり、（ｅ）補正係数Δｃを差し引いた修正後の画像を示す図であり、（ｆ）最終画像を示す説明図である。

【図5】実施形態１の超音波ＣＴ装置の動作の概要を示すフローチャートである。

【図6】実施形態１の超音波ＣＴ装置の詳しい動作を示すフローチャートである。

【図7】実施形態２の超音波ＣＴ装置のキャリブレーション判定機能を示すフローチャートである。

【図8】実施形態２の超音波ＣＴ装置のキャリブレーション判定機能の判定結果を示す表示画面例である。

【図9】実施形態２の超音波ＣＴ装置のキャリブレーション判定機能により素子の劣化判定の結果を示す表示画面例である。

【発明を実施するための形態】

【0016】

以下、本発明を一実施形態について、図面を参照して説明する。

【0017】

＜＜実施形態１＞＞
＜概要＞
まず、実施形態１の超音波ＣＴ装置の概要について図１等を用いて説明する。

【0018】

図１は、超音波ＣＴ装置の全体構造を示す図であり、図２は、水槽４に挿入された被写体１の乳房と振動子アレイ３の斜視図であり、図３は、振動子アレイ３０から超音波４０を照射することを示す図である。ここでは、乳癌検診に適した超音波ＣＴ装置について説明するため、乳房の画像を生成する例について説明するが、撮影対象は、乳房に限定されない。また、本発明は生体の音響情報の断層像を取得する超音波CT装置に限定するものではなく、例えば、非特許文献１に示されるような、地層の音響情報の断層像を取得する場合にも適用可能である。したがって、本発明では、超音波CT装置を物質の音響情報の断層像を取得する装置と定義する。

【0019】

図１のように、実施形態１の超音波ＣＴ装置は、送信部６と、受信部７と、画像再構成部８と、逐次更新部９と、入力受付部１０と、記憶部１１と、全体の動作を制御する制御部１２とを備えて構成される。送信部６および受信部７には、図２に示したように、複数の振動子３が、所定の形状（例えばリング状）に配列された振動子アレイ３０が接続されている。逐次更新部９には、表示部１１が接続されている。

【0020】

振動子アレイ３０は、図２に示すように円柱状の水槽４の内部または外部に配置されている。水槽４は、被写体１を乗せるベッド２に設けられた開口の下部に配置されている（図１参照）。被写体１は、ベッド２上でうつぶせになり、ベッド２の開口に乳房を挿入する。図示していないが、水槽４の軸方向に振動子アレイを平行移動することができる駆動機構が備えられていてもよい。振動子３は、超音波の送受信器であり、例えば圧電素子を用いる。水槽４には、温水が満たされ、予備タンク５が接続されている。予備タンク５は、水槽４の温水を浄化、過熱、脱気する機能を備えている。

【0021】

図３に示すように、送信部６は、振動子アレイ３０の１以上の振動子３ａに送信信号（電気信号）を出力し、１以上の振動子３ａから音波（ここでは超音波）４０を送信させる。被写体１の撮像領域を透過した音波４０は、振動子アレイ３０の複複数の振動子に到達し、その音圧が測定される。受信部７は、複数の振動子３が測定した測定音圧（電気信号）を受け取る。リング状の振動子アレイ３０の場合、撮像領域は被写体１のリング内に位置する領域である。受信部７は、測定音圧を受け取ってＡ／Ｄ変換する。

【0022】

画像再構成部８は、測定音圧を処理して、撮像領域の透過波画像を再構成する。透過波画像として、撮像領域における物性値分布（例えば、音速分布および／または減衰量分布）を示す画像を公知の方法により生成する。

【0023】

逐次更新部９は、画像再構成部８が生成した透過波画像を初期画像（図４（ａ）参照）として、計算により求めた音圧を用いて透過波画像を逐次的に修正する。すなわち、逐次更新部９は、模擬的な音源（以下模擬音源と呼ぶ）５３から時間的に音圧が変化する模擬的な音波（以下模擬音波）５０（音速ｃ）を発生させ（図４（ｂ）参照）、この模擬音波５０が撮像領域を透過して複数の模擬的な検出器（模擬検出器と呼ぶ）５４に到達する際の音圧（以下、計算音圧と呼ぶ）を、演算により求める。逐次更新部９は、例えば、図４（ｃ）のように、計算音圧と測定音圧との差分を求め、この差分を最小にする音速の補正係数Δｃ（図４（ｄ）参照）を求め、音速の補正係数Δｃを用いて透過波画像を逐次的に修正する（図４（ｅ），（ｆ）参照）。

【0024】

このとき、本実施形態１では、逐次更新部９は、計算音圧を用いて、逐次修正に用いる模擬音源５３の発生する模擬音波５０の波形を、振動子３ａが送信する音波４０の波形に近づける処理を行う。

【0025】

これにより、ＦＷＩ法においてシミュレーションにより模擬的な音源５３から発生する音波５０の音圧の波形を、実際の振動子３ａが照射する音波４０の音圧の波形に近づけることができるため、高精細な再構成画像を得ることができる。

【0026】

例えば、逐次更新部９は、模擬音波５０の時間的な音圧の変化の波形を複数種類に異ならせ、複数種類の波形についてそれぞれ計算音圧を求め、測定音圧と計算音圧との差分に基づいて、透過波画像の逐次的な修正に用いる模擬音源５３の発生する模擬音波５０の波形を選択する。

【0027】

この、模擬音波５０の波形の選択は、上述の透過波画像の逐次的な修正を行う前に１回のみ行ってもよいし、逐次的な修正処理における画像の更新のたびに行ってもよい。
なお、ここでいう模擬音波５０の波形とは、音圧の時間変化を表す形状であればいかなるものも含み、例えば、位相、最大振幅、周波数の他、波の時間変化が直線（矩形波や三角波等）か曲線（正弦波等の関数）か、曲線である場合には曲線形状等が異なるものは異なる波形とする。

【0028】

例えば、逐次更新部９は、模擬音波５０の複数種類の波形として、位相が異なる波形を用意し、その中から特定の位相の波形を選択する構成とする。

【0029】

選択方法としては例えば、逐次更新部９は複数種類の波形について、測定音圧と計算音圧との差分を表す所定の関数の値を算出し、その関数の値が最小となる波形を選択する。関数としては、例えば測定音圧と計算音圧との差の二乗または二乗和を算出する関数を用いる。

【0030】

なお、逐次更新部９は、模擬音源５３を、送信部６が音波を送信させる振動子３ａの位置に配置し、模擬検出器５４を受信部７が測定音圧を受け取る振動子３の位置に配置することが望ましい。

【0031】

以下、超音波ＣＴ装置の各部の動作について、図５等のフローを用いて説明する。

【0032】

なお、本実施形態では、画像再構成部８および逐次更新部９は、ＣＰＵ（Central Processing Unit）やＧＰＵ（Graphics Processing Unit）等のプロセッサーと、メモリとを備えたコンピュータ等によって構成され、ＣＰＵが、メモリに格納されたプログラムを読み込んで実行することにより画像再構成部８および逐次更新部９の機能をソフトウエアにより実現する構成とする。ただし、画像再構成部８および逐次更新部８は、その一部または全部をハードウエアによって実現することも可能である。例えば、ＡＳＩＣ（Application Specific Integrated Circuit）のようなカスタムＩＣや、ＦＰＧＡ（Field-Programmable Gate Array）のようなプログラマブルＩＣを用いて画像再構成部８および逐次更新部８を構成し、その動作を実現するように回路設計を行えばよい。

【0033】

図５のフローを用いて、超音波ＣＴ装置の動作の概要を説明する。

【0034】

まず、ステップ５０１において、送信部６は、振動子３ａから音波４０を被写体１に対して送信する。被写体１の撮像領域を透過した音波４０は、振動子３により受信される。受信部７は、振動子３から測定音圧ｕ_ｍを受け取る。

【0035】

ステップ５０２において、画像再構成部８は、測定音圧ｕ_ｍを用いて、透過波画像を生成する。ここでは、透過波画像として、音速分布画像（以下、音速画像と呼ぶ）を生成する例について以下説明する。

【0036】

ステップ５０３において、逐次更新部９は、模擬音源５３から模擬音波５０を被写体１に対して照射した場合に模擬検出器５４によって検出される計算音圧ｕ_ｓを算出し、計算音圧ｕ_ｓを用いて模擬音波５０の波形を音波４０に近づける処理を行う。

【0037】

ステップ５０４において、逐次更新部９は、所定の目的関数の値が閾値以上である場合、ステップ５０５に進み、ステップ５０３で音波４０に近づけた模擬音波５０を用いて、音速画像を修正する。

【0038】

逐次更新部９は、ステップ５０３～５０５をステップ５０４で目的関数の値が閾値より小さくなるまで繰り返す。

【0039】

これにより、波形を音波４０に近づけた模擬音波５０を用いて、音速画像を修正できるため、精度よく高解像度の再構成画像を得ることができる。

【0040】

図５のステップ５０１～５０５の処理の詳細を図６のフローを用いて説明する。

【0041】

＜ステップ５０１＞
ステップ５０１の処理を、図６のステップ５１１～５１３により詳しく説明する。

【0042】

ステップ５１１において、本実施形態の超音波ＣＴ装置において、入力受付部１０がユーザから撮像開始の指示を受けた場合、送信部６は、記憶部１１に予め格納しておいた振動子３ａをビューごとに特定し、振動子３ａに対して、予め定めた周波数（数ＭＨｚ程度）の送信信号を出力する。送信信号を受け取った振動子３ａは、音波４０を被写体１に対して送信する。

【0043】

ステップ５１２において、被写体１に照射された音波４０の一部は、被写体１を透過し、振動子３により受信され、電気信号（測定音圧ｕ_ｍ）に変換される。受信部７は、振動子３から測定音圧ｕ_ｍを受け取って、Ａ／Ｄ変換等の処理を行う。

【0044】

ステップ５１１，５１２を超音波の送信および受信を、全ビューから行う（ステップ５１３）。

【0045】

＜ステップ５０２＞
ステップ５０２の処理を、図６のステップ５１４、５１５により詳しく説明する。

【0046】

ステップ５１４において、画像再構成部８は、測定音圧ｕ_ｍを用いて画像再構成処理を行うことにより、被写体１内における音速画像を生成する。ステップ５１５において、音速画像は、逐次更新部９が行う逐次更新処理の初期画像ｃ^ｎ（ｎ＝０）となる。

【0047】

具体的な画像再構成処理例について簡単に説明する。例えば、画像再構成部８は、ストレートレイ法により音速画像を求めることができる。すなわち、画像再構成部８は、各ビューにおける各振動子３の出力した測定音圧に対して、時間方向にヒルベルト変換を実施し、受信波の最大振幅の受信タイミングを求める。画像再構成部８は、被写体１の挿入前に予め受信しておいた各振動子３の測定音圧についても同様に最大振幅の受信タイミングを求める。画像再構成部８は、被写体１の挿入前後の受信タイミングの差を、各ビュー、各受信チャネルについてそれぞれ計算し、それらデータの集合であるサイノグラムを得る。画像再構成部８は、受信タイミングの差のサイノグラムをフィルタ補正逆投影法（Filtered Back Projection, FBP）等で処理することにより、断層画像を再構成する。この断層画像は、被写体１の挿入前後の、超音波の「遅さ（Slowness）」の差の分布画像である。「遅さ」は、音速の逆数である。画像再構成部８は、水の音速値（推定値）を用いて、「遅さ（Slowness）」の差の分布画像から、被写体１の音速の分布画像（音速画像）を生成する。

【0048】

＜ステップ５０３～５０５＞
ステップ５０３において、ステップ５０４、５０５における音速画像の逐次更新処理に用いる模擬音波５０の波形を、実際に送信した音波４０の波形に近づける処理を図６のステップ５１６～５２０により行う。その説明のために、まず、ＦＷＩ法による逐次更新処理の原理について説明する。

【0049】

＜ＦＷＩ法の逐次更新の原理＞
ＦＷＩ法は、計算音圧ｕ_ｓを算出する音圧計算アルゴリズムと、目的関数を最小にするように透過波画像を修正する逆問題のアルゴリズムとからなる。

【0050】

音圧計算アルゴリズムは、音速画像が表す空間の振動子３ａの位置に模擬音源５３を配置し、模擬音源５３から模擬音波５０を発生し、音速画像が表す空間を透過し、振動子３の位置に配置した模擬検出器５４に到達した模擬音波５０の音圧（計算音圧ｕ_ｓ）をシミュレーションにより求める。

【0051】

ＦＷＩを周波数領域で実行する場合には、模擬音源５３を以下の式（１）のように複素数で表す。

【0052】

【数1】

【0053】

ここで、Aは模擬音源の強度および位相を調整する複素数の係数、iは虚数、θは位相である。模擬音源の調整係数Aは、例えば公知のPrattの信号推定法により求めることができる（非特許文献１）。

【0054】

音圧は、例えば、周波数領域のヘルムホルツ方程式を有限差分法で解く方法を用いることができる。周波数領域におけるヘルムホルツ方程式は次の式で表わされる。

【0055】

【数2】

【0056】

ここで、ωは角周波数、rは位置である。c(r)は、音速画像（撮像領域）の位置ｒの画素の画素値（音速）である。u(r,ω)は、音速画像の位置ｒの音圧を表すベクトルであり行列である。f(r)は、式（１）の関数で表される模擬音源であり、ｒは、音源の位置である。

【0057】

式（２）において

【数3】

はラプラシアンである。

【0058】

空間を離散化し、式（２）を有限差分法による行列表現で表わすと以下のようになる。

【0059】

【数4】

【0060】

式（４）において、S(r)は、インピーダンス行列と呼ばれ、有限差分法における音圧u(r)の係数行列を表す。すなわち、S(r)は、式（５）に示すように、式（２）における音圧u(r)の左側の括弧内の行列を表している。なお、式（４）において、周波数ωは一定であると仮定し、ωは省略している。

【数5】

【0061】

式（４）の右側の第２式から明らかなように、式（１）で表される音源f（r）と、インピーダンス行列Ｓ（r）の逆行列Ｓ（r）^－１から、音速画像の位置ｒの音圧ｕ(ｒ)を計算により求めることができる。計算により求められた位置ｒの音圧を行列で計算音圧ｕ_ｓ(ｒ)と表す。

【0062】

つぎに、逆問題のアルゴリズムの解法として、式（４）の計算で得られる計算音圧ｕ_ｓ(ｒ)と、ステップ５１２で測定した測定音圧ｕ_mの残差二乗和である目的関数E(c)を式（６）から求め、その最小化を行う。

【0063】

【数6】

式（６）において、Hはエルミート転置を表す。

【0064】

式（６）の目的関数E(c)を最小化するために、音速画像の画素値ｃは、例えば以下の式（７）で表される最急降下法により逐次的に修正することができる。

【0065】

【数7】

【0066】

式（７）において、ｎは反復回数であり、ｃ^ｎは、修正前（ｎ回目）の音速画像上の画素値（音速または減衰量）、ｃ^ｎ＋１は、修正後（ｎ＋１回目）の透過波画像上の画素値であり、αはステップ長と呼ばれる修正量を調整するパラメータである。

【0067】

ここでは、式（７）の右辺第2項を式（８）のように以下Δｃで表す。

【0068】

【数8】

【0069】

なお、上記説明においては、最急降下法により音速画像を修正しているが、共役勾配法などの他のアルゴリズムを用いることも可能である。また、本実施形態のＦＷＩは、周波数領域、時間領域のいずれの領域でも実行することが可能である。

【0070】

＜ステップ５０３＞
上記原理を用いて、ステップ５０３において模擬音波５０の波形を音波４０の波形に近づける処理を図６のステップ５１６～５２０により具体的に説明する。

【0071】

まず、ステップ５１６において、逐次更新部９は、ステップ５１５の音速画像の画素値の示す音速c(r)（ｒは、画素の位置）から、式（５）によってインピーダンス行列S(r)を算出し、その逆行列S(r)^－１を算出する。
ステップ５１７において、模擬音源５３を式（１）によりｆで表す。逐次更新部９は、算出したS(r)^－１と、模擬音源５３であるｆと、式（４）の右側の第２式から、模擬検出器５４の位置の計算音圧ｕ_ｓを算出する。

【0072】

つぎに、ステップ５１８において、逐次更新部９は、上記模擬音源５３の計算音圧ｕ_ｓに、式（９）により表されるＮ種類の係数ａ(k)（ｋ＝０，１，・・・，Ｎ－１）をそれぞれ掛けることにより、Ｎ種類の計算音圧ａ(k)ｕ_ｓを生成する。式（９）から明らかなようにＮ種類の係数ａ(k)は、それぞれ異なる角度θｋの複素数であり、計算音圧ｕ_ｓの位相をずらす作用をする係数である。

【0073】

【数9】

【0074】

つぎに、ステップ５１９において、逐次更新部９は、計算音圧ａ(k)ｕ_ｓを式（６）の計算音圧ｕ_ｓ(c)と置き換えた式（６’）により、ステップ５１８で算出したＮ種類の計算音圧ａ(k)ｕ_ｓと、ステップ５１２で測定した測定音圧ｕ_ｍから、目的関数Ｅ（k）を算出する。これにより、Ｎ種類の目的関数Ｅ（k）が得られる。

【0075】

【数10】

【0076】

ステップ５２０において、逐次更新部９は、Ｎ種類の目的関数Ｅ（k）のうち最小値となる目的関数Ｅ（k'）（ｋ’は、最小値となる目的関数Ｅ（k）の時のｋの値）を選択し、その最小値の目的関数Ｅ（k'）が得られた係数a(k')を選択する。

【0077】

上述の式（４）の右側の第２式より、インピーダンス行列S(r)が同一である場合、すなわち、同一被写体１である場合、式（４）の第２式の両辺にa(k')を掛けても式（４）は成り立つ。よって、計算音圧ｕ_ｓにa(k')を掛けたa(k')ｕ_ｓは、fで表される模擬音源５３にa(k')を掛けた音源a(k')fによって得られることがわかる。すなわち、最小の目的関数E（k'）となる音圧a(k')ｕ_ｓは、模擬音源５３の位相を調整したa(k')fによって得られる。

【0078】

このように、上記ステップ５１６～５２０により、模擬音源５３から発生する音圧波形５０を振動子３ａの発生する音波４０の位相に最も近づけることのできる模擬音源５３の位相を調整したa(k')fを求めることができる。

【0079】

＜ステップ５０４＞
図５のステップ５０４の目的関数を求める処理を、図６のステップ５２１において説明する。

【0080】

ステップ５２１において、逐次更新部９は、計算音圧a(k')u_ｓと、ステップ５１２で得た測定音圧ｕ_ｍと、式（６）から、目的関数E(k')を算出する。逐次更新部９は、算出した目的関数E(k')が予め定めた閾値よりも小さい場合、処理を終了し、ステップ５１５の音速画像ｃ^ｎが最終画像となる。目的関数E(k')が予め定めた閾値よりも大きい場合、ステップ５０５（図６のステップ５２２～５２３）に進み、音速画像を修正する。

【0081】

＜ステップ５０５＞
図５のステップ５０５の音速画像の修正処理を、図６のステップ５２２～５２３において説明する。

【0082】

ステップ５２２において、逐次更新部９は、例えば最急降下法によりE(ｃ)を最小化する音速ｃを算出する。具体的には、音速画像の画素ごとに、式（８）により、補正係数Δｃを算出する。
ステップ５２３において、逐次更新部９は、式（７）のように、画素値ｃ^ｎからΔｃを差し引くことにより修正後の画素値cⁿ⁺¹を算出する。これにより、修正後の音速画像cⁿ⁺¹を得る。

【0083】

逐次更新部９は、ステップ５１５に戻り、ステップ５２３で算出した修正後の音速画像cⁿ⁺¹に音速画像ｃ^ｎを置き換え、ステップ５２１で目的関数E(k')が閾値thより小さくなるまで、ステップ５１５～ステップ５２３を繰り返すことにより逐次的に音速画像を修正する。

【0084】

逐次更新部９は、ステップ５２１で目的関数E(k')が閾値thより小さくなった場合、逐次修正が収束したと判断して、透過波画像を表示部１３に表示させる。

【0085】

上述してきたように、本実施形態では、FWI法において、模擬音源５３の発生する音波５０の波形（例えば位相）を、振動子３ａが発生する音波４０に逐次的に近づけることができるため、高い空間分解能を有する音速画像を生成することができる。

【0086】

なお、図５および図６のフローでは、逐次修正のたびに、ステップ５０３（ステップ５１６～５２０）を行って係数a(k')を求める構成について説明したが、１回目だけ係数a(k')を算出し、２回目以降は、従来と同様に逐次修正を行う構成としてもよい。

【0087】

なお、ステップ５２１において、上述の実施形態では、目的関数E(ｃ)が予め定めた閾値th以下になった場合に逐次修正が収束したと判断しているが、これ以外の判断基準を用いることも可能である。例えば、上記ステップ５１５～ステップ５２３を予め定めた回数繰り返した場合には、逐次修正が収束したと判断することができる。

【0088】

また、上述の実施形態では、画像再構成部８は、透過波画像として、音速画像を生成し、音速画像を逐次修正する構成について説明したが、減衰量分布画像（減衰画像）を生成し、同様に減衰画像を逐次修正することも可能である。

【0089】

なお、振動子アレイ３０を上下（水槽４の軸方向）に移動させながら、受信部７が測定音圧を収集することにより、画像再構成部８は、被写体１の透過波画像を三次元画像として取得することもできる。

【0090】

＜＜実施形態２＞＞
実施形態２の超音波CT装置について説明する。実施形態２の超音波CT装置は、実施形態１と同様の構成であり、実施形態１と同様の機能に加えて、キャリブレーション判定機能を備えている。

【0091】

FWIのシミュレーションでは、振動子アレイ３０の複数の振動子３ごとの位置（座標）と、振動子３を構成する素子の機能のばらつきにより生じる音波４０の送受信タイミングを較正するための振動子３ごとの遅延時間を予め求めておく必要がある。これらの振動子位置と遅延時間は、キャリブレーション情報と呼ばれ、キャリブレーション測定によって予め求められている。キャリブレーション情報（振動子位置、遅延時間）が正確である場合、測定音圧と計算音圧は、精度よく一致する。

【0092】

実施形態２の超音波CT装置は、キャリブレーション判定機能により、被写体１を配置せずに、水に対して、音波４０を送受信して測定音圧を測定し、計算音圧と比較する。このとき、実施形態１と同様に、模擬音源５３から発生する模擬音波５０の波形を、音波４０に近づける処理を行うことにより、キャリブレーション情報が正確かどうかを精度よく判定することができる。

【0093】

キャリブレーション情報の判定を行う際の超音波CT装置の各部の動作について図７のフローを用いて説明する。図７のフローにおいて、図６のフローと同様の動作のステップについては同じ番号を付し、簡単に説明する。

【0094】

ステップ５１１おいて、水槽４内に被写体１を配置せず、水のみに対して振動子３から音波４０を送信し、ステップ７１２において、測定音圧ｕ_ｍを得るとともに、水温を記録する。全ビューについて、ステップ５１１，７１２を行う。

【0095】

ステップ７１４において、ステップ７１２で測定した水温から、Del Grossの式等により水の音速を算出し、算出した音速値を一様にもつ音速画像を生成する。

【0096】

生成した水の音速画像を用いて、実施形態１と同様にステップ５１７～５２０を行うことにより、模擬音源５３から発生する音波５０の波形を音波４０に近づける係数a(k')を選択する。

【0097】

ステップ７２１において、計算音圧a(k'）u_ｓと、ステップ７１２で水を測定したときの測定音圧u_ｍの位相差を算出する。この位相差が０に近いほどキャリブレーション情報の精度が高く、位相差が大きくなるほどキャリブレーション情報の精度が低いと判定することができる。よって、位相差が予め定めて閾値ｔｈ’以上である場合、ステップ７２２に進み、算出した計算音圧a(k'）u_ｓと測定音圧u_ｍの位相差と、振幅差（音圧差）を例えば、図８のように表示部１３に表示するとともに、キャリブレーションの実行を促す表示をする。

【0098】

このように、実施形態２の装置では、水を測定して上述のようにキャリブレーション情報の精度判定を行うことができる。これを定期的（例えば月１回）に行うことにより、経年劣化により、再キャリブレーションが必要かどうかをユーザーに知らせることが可能である。

【0099】

実施形態２では、計算音圧a(k'）u_ｓと測定音圧u_ｍの位相差だけでなく、振幅差（音圧差）が閾値以上であるかどうかも判定することができるため、位相差および／または振幅差が閾値以上である振動子３を、異常値を出力する劣化した不良素子であるとして、図９のように、素子交換を促す表示を行うこともできる。

【0100】

さらに、素子交換までの間、不良素子である振動子３の出力（測定音圧ｕ_ｍ）をマスクし、画像の生成に用いないようにすることができる。具体的には、図６のフローにおいて、測定音圧ｕ_ｍを用いるステップ５１４、５１９および５２１では不良素子の振動子３の測定音圧u_ｍをマスク処理し、演算に用いない。

【0101】

また、超音波CT装置を据え付けた後のキャリブレーション測定の後で行うことにより、キャリブレーション情報の精度判定を行うことができる。

【符号の説明】

【0102】

１…被写体、２…ベッド、３…振動子、４…水槽、５…予備タンク、６…送信部、７…受信部、８…画像再構成部、９…逐次更新部、１０…入力受付部、１１…記憶部、１２…制御部、１３…表示部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】