JP5584300B2 - 画像再構成方法及び装置 - Google Patents

Description

本発明は、計測対象とする計測領域との間で信号を送受信し、その実受信信号に基づいて計測領域の物性を定量する画像再構成方法及び装置に関する。本発明は、計測対象内部の密度などの物性値が不均一に分布している場合にも、不均一性による屈折や多重反射などの現象に起因する画質劣化なく計測領域内の物性値の分布を実現可能な時間で算出可能とする技術に関する。

計測装置とは計測対象の内部を非侵襲に画像化する装置で、計測対象の境界から計測対象の内部に測定信号（例えば超音波）を送信し、測定信号の伝播現象を記述する物理方程式に基づいたアルゴリズムによって、境界で測定された測定値を内部の物性値分布に変換する。測定信号は物理方程式の変数に、物性値は物理方程式の係数や外力項に相当し、物理方程式は測定信号と物性値の関係をあらわす式として使われる。

物理方程式に基づいて測定信号を物性値に変換するアルゴリズムは、科学計算の分野で係数同程問題や境界値問題として研究されているが、装置に実装されるアルゴリズムは、計算量を実装可能な範囲に減らすために方程式を近似して変換の精度を落としている。現在製品化されている装置のアルゴリズムは、アンテナ技術を基本とした指向性を合成する方法で、測定対象内の物性値を均一だと近似している。CTと同様な投影写像のアルゴリズムを搭載した先端的な試みである乳房専用装置でも、測定信号の直線性を仮定する意味で物性値が均一だと近似している。

また、特許文献１にあるように、生体の物性値の不均一を考慮し、ヘルムホルツ方程式などの波の伝播を考慮して、送信波形から受信波形を計算する。実測の受信波形と計算された受信波形の差が小さくなるように、生体の物性値を修正するステップを繰返すことで、生体の物性値を推定する方法もあった。

特表2006-511298

Breast Imaging using transmission Ultrasound: Reconstructing Tissue Parameters of Sound Speed and Attenuation ;Cuiping Li; Duric, N.; Lianjie Huang; BioMedical Engineering and Informatics, 2008. BMEI 2008. International Conference on ;Vol.2 pp.708-712（2008）

しかし、上述した非特許文献１に代表される方法は近似した方程式に基づいているため、処理時間が短いが、計測対象の内部構造が不均一の場合には、不均一性による屈折や多重反射などの現象に起因する画質劣化を取り除くことが原理的にできないという問題がある。

また、特許文献１に代表される方法では、計測対象の内部での波の伝播を計算するのに時間が掛るのみでなく、この計算を実測の受信データと計算された受信データのずれが収束するまで繰返し計算が必要なため、リアルタイム性が求められる医療装置としては、実用的でないという課題があった。

本発明者は、計測対象の境界から対象内に測定信号を送受信し、境界での測定値から計測対象内部の物性値を算出する方法及び装置であって、計測対象の内部が不均一の場合にも、不均一性による屈折や多重反射などの現象に起因する画質劣化なく、実装可能な計算規模で、計測対象内の物性値の分布を算出することができる技術を提案する。

特に、実測の受信データから、物性値の分布を推定する演算の規模を大幅に短縮する技術を提案する。

本発明は、近似していない物理方程式から導出される評価量を使用し、測定対象の物性値を求める。本発明では、仮定した物性値を真の物性値に近づけるであろうずれ量を推定し、当該ずれ量の更新を通じて真の物性値を発見する手法を採用する。なお、画像再構成に寄与の大きい信号のみを取得する送信シーケンスを用いることが好ましい。具体的には、(a) 本発明は、計測対象とする計測領域との間で信号を送受信し、計測領域の物性を定量する画像再構成方法として、以下のステップを有するものを提案する。

(1) 計測領域との間で信号を送受信する送受信データ取得ステップ（手段）
(2) 計測領域から受信された実受信信号から物性値を算出する物性値算出ステップ（手段）
(3) 物性値を表示値に換算する物性値−表示値換算ステップ（手段）
例えば物性値算出ステップは、(i) 運動方程式の演算子項と外力項の差分の残差である方程式残差量と、(ii)２種類の境界条件による運動方程式の解の一致度により物性値の不均一性を検出する式の残差である不均一検出式残差量と、(iii) 制約条件の残差である条件式残差量との線形和を評価量とし、評価量を極値にする物性値を物性値−表示値換算ステップ（手段）に出力する。

(b) なお、前述の送受信データ取得ステップは、(i) 送受信素子の送受信回数と、(ii) 送受信素子の空間配置と、(iii) 送受信波形とを決定して送受信結果を取得する送受信シーケンス決定ステップを有することが望ましい。

ここで、受信シーケンス決定ステップが決定する送受信素子の送受信回数は２回以上であり、送信波形は0.5MHz以上10MHz以下のキャリア周波数が異なるバースト長の等しい複数の平面バースト波を少なくとも含み、送信素子の空間配置はキャリア周波数に逆比例する間隔であり、受信素子の空間配置は計測領域に対する見込み角を最大にする配置であることが望ましい。

(c) または、前述の送受信シーケンス決定ステップが決定する送受信回数は２回以上であり、送受信時刻が早いほど送受信の周波が低い送受信シーケンスを出力することが望ましい。

(d) または、前述の不均一検出式残差量は、(i) 境界値が実送受信データの場合のノイマン境界条件による解とディリクレ境界条件による解の差分である差分ベクトルと、(ii)境界値がゼロの場合のノイマン境界条件に対する運動方程式のグリーン関数とディリクレ境界条件に対する運動方程式のグリーン関数の差分である演算子行列と物性値の不均一性を表すベクトルの積との差分の残差であることが望ましい。

(e) または、前述の評価量のうち条件式残差量は、評価量物性値分布の性質を表す事前確率である事前確率量、及び、指向性を合成する方法若しくは投影写像により再構成された画像との一致度を表す低次近似一致度量の少なくとも一方を含むことが望ましい。

(f) または、前述の物性値算出ステップは、評価量を極値にする物性値を、確率的手法の１つであるメトロポリスサンプリングを用いて算出することが望ましい。

(g) 評価量または、前述の方程式残差は、動弾性方程式の周波数空間表示、動弾性方程式の時間応答、ヘルムホルツポテンシャル表示値のいずれかの方程式の残差であることが望ましい。また、物性値は、ラメ定数λ、ラメ定数μ、慣性密度ρ、減衰係数ｃ、これらの比、力源ｆのうちの１つ以上であることが望ましい。また、表示量は、λ、μ、ρ、ｃ又はこれらの関数である体積弾性率、縦波速度、横波速度、ポアソン比、ヤング率、インピーダンス、力源ｆのうち１つ以上であることが望ましい。

本発明は、近似していない物理方程式を用いるため、不均一性による屈折や多重反射などの現象に起因する画質劣化のない画像再構成が可能になる。また、画像再構成の処理速度は、信号取得時間と計算時間によって決まるが、画像再構成に寄与の大きい信号のみを取得する送信シーケンスを用いることで信号取得時間を短縮し、かつ、仮定した物性値の真の物性値からのずれ量を推定する式を用いて計算時間を短縮することで、画像再構成の処理時間を実現可能な長さにすることができる。

実施例に係る画像再構成装置のブロック構成を示す図。 実施例に係る画像再構成方法の処理の流れを示すフローチャート。 実施例に係る画像再構成方法の処理の流れを示すフローチャート。 実施例に係る画像再構成装置の探触子と計測対象の関係を概念的に説明する図。 実施例に係る画像再構成装置の探触子と計測対象の関係を概念的に説明する図。 実施例に係る画像再構成装置の物性値算出部の処理の概要を説明する図。 実施例に係る画像再構成装置の物性値算出部の処理の概要を説明する図。 実施例に係る画像再構成装置の表示部に表示される画面を概念的に示す図。 実施例に係る画像再構成装置のデータフォーマットの一例を概念的に示す図。 実施例に係る画像再構成装置の評価量の違いと物性値の求め方の関係を概念的に示す図。 実施例に係る画像再構成装置の送受信データ取得部が出力する送信シーケンスの基本シーケンス例を概念的に示す図。 実施例に係る画像再構成装置の送受信データ取得部が出力する送信シーケンスの一例を概念的に示す図。 実施例に係る画像再構成装置の送受信データ取得部が出力する送信シーケンスの一例を概念的に示す図。 実施例に係る画像再構成装置の送受信データ取得部が出力する送信シーケンスの他の一例を概念的に示す図。 データの処理フローの詳細図。 引例と本発明の処理ブロックの比較図。

以下、図面に基づいて、本発明の実施の形態を説明する。なお、後述する装置構成や処理動作の内容は発明を説明するための一例である。本発明は、後述する装置構成同士や処理動作同士の任意の組み合わせ、後述する装置構成や処理動作に既知の技術を追加する組み合わせ、後述する装置構成や処理動作の一部を既知の技術で置換する組み合わせも包含する。

（実施例）
まず本発明の特長を図１３を用いて、概略を説明する。特許文献１の方法は、生体の物性値の不均一を考慮し、ヘルムホルツ方程式などの波の伝播を考慮して、送信波形から受信波形を計算する。実測の受信波形と計算された受信波形の差が小さくなるように、生体の物性値を修正するステップを繰返すことで、生体の物性値を推定する方法であり、装置構成と処理フローは図１３（ａ）に示されるとおりである。

一方、本発明では、図１３（ｂ）に示す構成となり、実際の送受信データと、仮定された物性値の分布から確からしい物性値と仮定された物性値のずれ量を一回の計算ステップで計算する。このずれ量と、仮定された物性値の和を確からしい物性として出力する。この演算では、特許文献１では大規模であった演算規模を大きく改善する事ができる。もちろん、ノイズ等の影響によって、一度の計算ステップでは信頼性が課題となる場合には、リアルタイム性を損ねない範囲で、演算を繰返し、信頼性を向上することも有効な方法である。

（画像再構成装置の構成）
以下、図１〜図１１を参照し、実施列を説明する。図１は、本発明の画像再構成装置の構成を示すブロック図である。実施例に係る画像再構成装置は、本体１と探触子２から構成される。実施例に係る本体１はコンピュータ構成を有し、数値計算その他の信号処理はコンピュータプログラムを通じて実現する。

探触子２は複数の送受信素子２１を含む。本体１は送信部３、受信部４、信号処理部５、メモリ６、表示部７、入力手段８、制御部９を含む。信号処理部５は更に、計測領域との間で計測信号を送受信し、受信された信号である実受信信号を出力する送受信データ取得部５１、前記送受信データ取得部から前記実受信信号を入力し、計測領域の物性値に変換して出力する物性値算出部５２、前記物性値算出部から前記物性値を入力して表示量に換算する物性値−表示量換算部５３を含む。図中太枠で示す物性値算出部５２が、今回提案する処理機能の中でも特に重要な処理機能を提供する。

（画像再構成装置による計測動作）
実施例に係る画像再構成装置は、次の順序で信号処理を実行する。
まず、複数の送受信素子２１のそれぞれに対応する送信波形である送信シーケンスが、信号処理部５の送受信データ取得部５１により決定される。決定された送信シーケンスは、送受信データ取得部５１から送信部３に直接出力される、又はメモリ６に書き込まれる。メモリ６に書き込まれた送信シーケンスは、所定のタイミングでメモリ６から読み出され、送信部３に出力される。

次に、表示部７の表示画面を目視しながら計測作業を行う作業者が、計測領域の撮像の開始を、マウス、キーボード、操作ボタンその他の入力手段８を通じて操作入力する。当該操作入力を検知した入力手段８は、撮像開始信号を制御部９に出力する。制御部９は、撮像開始信号を入力すると、撮像処理の開始を信号処理部５に通知する。

送信部３は、信号処理部５又はメモリ６から送信シーケンスを入力すると、それぞれ対応する送受信素子２１から計測対象に信号を送信する。例えば生体等の撮像対象に超音波信号を送信する。

送信シーケンスの送信開始と同時に、受信部４は、撮像対象から反射又は撮像対象を透過した信号に対応する受信信号を送受信素子２１により受信する。受信部４は、個々の送受信素子２１に対応する受信信号を信号処理部５に出力する。

まず、信号処理部５は、送受信データ取得部５１において、送受信素子２１のそれぞれに対応する送受信データを入力する。次に、信号処理部５は、物性値算出部５２において、送受信素子２１のそれぞれに対応する送受信データに基づいて、ラメ定数λ、μ、慣性密度ρなどの物性値を算出する。次に、信号処理部５は、物性値−表示量換算部５３において、ラメ定数λ、μ、慣性密度ρ等の物性値を、物性値そのもの、縦波音速（λ+2μ）/ρ、横波音速μ/ρ、体積弾性率（λ+2μ）/ρその他の表示量に換算し、メモリ６に出力する。

メモリ６は書き込まれた表示量をモニタその他の表示部７に出力する。表示部７は、各物性値に対応する表示量を画面上に表示する。

（処理動作の概要）
図２−１に、本実施列に係る画像再構成方法の処理の流れを説明するフローチャートを示す。画像再構成方法を構成する各ステップは、計測装置（例えば計算機）の数値演算部における数値計算を通じて実行される。各ステップで使用する数値及び計算結果は記憶領域（メモリ）に保存される。これらの数値は、数値演算の進行に応じ、プログラムの指示に従い、数値演算部と記憶領域（メモリ）との間で受け渡しされる。

図２−１の（ａ）は画像再構成方法の処理の概要を示し、図２−１の（ｂ）はステップＳ１００の具体的な処理内容を示している。

図２−１の（ａ）に示すように、本実施例に係る画像再構成方法においては、送受信データ取得ステップ（Ｓ１００）と、物性値算出ステップ（Ｓ２００）と、物性値−表示量換算ステップ（Ｓ３００）とで構成される。送受信データ取得ステップ（Ｓ１００）は、計測領域との間で計測信号を送受信し、受信された信号である実受信信号を物性値算出ステップに出力する。物性値算出ステップ（Ｓ２００）は、前述した実受信信号から計測領域の物性値を算出し、物性値−表示量換算ステップに出力する。物性値−表示量換算ステップ（Ｓ３００）は、前述した物性値を表示量に換算する処理を実行する。図中太枠で示す物性値算出ステップ（Ｓ２００）が、今回提案する処理機能の中でも特に重要な処理機能を提供する。

（送受信データ取得ステップの内容）
図２−１の（ｂ）に、送受信データ取得ステップ（Ｓ１００）の詳細な処理内容を示す。送受信データ取得ステップ（Ｓ１００）は、送受信シーケンス決定ステップ（Ｓ１０１）と、送信ステップ（Ｓ１０２）と、受信ステップ（Ｓ１０３）で構成される。送受信シーケンス決定ステップ（Ｓ１０１）では、送受信シーケンスの条件の決定処理が実行される。送信ステップ（Ｓ１０２）では、計測領域に計測信号（例えば超音波）を送信する処理が実行され、受信ステップ（Ｓ１０３）では、計測領域から受信信号を受信する処理が実行される。図中太線で示す送受信シーケンス決定ステップ（Ｓ１０１）では、物性値算出ステップ（Ｓ２００）における数値計算を効率化するための数値条件が決定される。

（物性値算出ステップの内容）
図２−２の（ａ）に、物性値算出ステップ（Ｓ２００）の詳細な処理内容を示す。物性値算出ステップ（Ｓ２００）は、物性値の数値計算に後述する評価量Ｊを使用する。本実施例の場合、計測対象内部の物性値は、評価量Ｊの極値として求められる。

物性値算出ステップ（Ｓ２００）が開始すると、まず、初期物性値設定ステップ（Ｓ２０１）が物性値の初期値を設定する。次に、送受信データ算出ステップ（Ｓ２０２）が物性値の初期値に対する受信信号である実受信信号を計算する。次の収束判定ステップ（Ｓ２０３）は、評価量Ｊが極値か否かにより計算結果が収束するか否かを判定する。ステップ（Ｓ２０３）の判定結果が否定（Ｎｏ）である場合、物性値更新ステップ（Ｓ２０４）は物性値を更新し、更新後の物性値について計算受信データ算出ステップ（Ｓ２０２）に戻る。このループ処理は、ステップ（Ｓ２０３）の判定結果が肯定（Ｙｅｓ）になるまで繰り返される。ステップ（Ｓ２０３）の判定結果が肯定（Ｙｅｓ）である場合の物性値は、物性値−表示量換算ステップ（Ｓ３００）に出力される。図２−２の（ａ）に示す処理ステップのうち太枠で示す収束判定ステップ（Ｓ２０３）と物性値更新ステップ（Ｓ２０４）が新規な処理を含むステップである。このステップは、評価量Ｊとその極値の求め方に関する処理を含む。

この実施例の場合、「評価量」は、(1) 運動方程式の演算子項と外力項との差分の残差である方程式残差量と、(2) ２種類の境界条件による運動方程式の解の一致度により物性値の不均一性を検出する式の残差である不均一検出式残差量と、(3) 制約条件の残差である条件式残差量との線形和として与えられる。なお、評価量は、これらの残差量のうちの１つ又は任意の複数を含む線形和若しくはそれらの積若しくはそれらを指数に含む指数関数として与えることもできる。

（評価量の内容）
以下では、評価量の内容について説明する。最初に、評価量Ｊが収束するか否かを、評価量Ｊが極値（この実施例ではゼロ）になるか否かにより判定する場合について説明する。しかし、この判定は、「不均一検出式残差量」Ｊの指数関数P＝exp（-J^2）が最大になることと同値である。すなわち、評価量Ｊによる変分法は、確率密度を表現する指数関数Ｐの最大化処理として表現することもできる。

運動方程式は、計測信号が超音波である場合、動弾性方程式の周波数表示、又は、動弾性方程式の時間応答表示、又は、動弾性方程式のヘルムホルツポテンシャル表示で与えられる。この場合、物性値は、ラメ定数λ、ラメ定数μ、慣性密度ρである。因みに、計測領域が人体の場合、初期物性値は、計測領域の内部全域でλは2e9［kg/m/sec^2］、ρは1e3［kg/m^3］でそれぞれ一様として与えられる。

まず、基本的な評価量について説明する。基本的な評価量は、運動方程式の演算子項と外力項の差分の残差で、方程式残差量Ｊ_ｅｑと呼ぶ。運動方程式の演算子を行列Ｇ（ε）［ｉ］［ｊ］、変位をベクトルｐ、外力をベクトルｆで表すとする。Ｇは求めるべき物性値εの関数である。このとき、これらの量の間には関係Ｇ＊ｐ＝ｆが成り立つ。このため、方程式残差量Ｊ_ｅｑをＪ_ｅｑ＝Ｇ＊ｐ−ｆと定義すると、Ｊ_ｅｑは物性値εの関数となり、ε_Ａが真値の場合に方程式残差量Ｊ_ｅｑが０になる。方程式残差量を０にする物性値ε_Ａを求める方法には、一般的に知られている方法を使えばよい。例えばＪ_ｅｑ＝−（Ｇ＊ｐ−ｆ）^２として最急降下法を使うなどが例として挙げられる。

この評価量Ｊ_ｅｑを用いれば、理想的な（近似していない）物理方程式を満たす解を、評価量Ｊの極値（この実施例ではゼロ）を与える物性値として求めることができる。しかも、求められる物性値からは、屈折や散乱による影響を除去することができる。

次に、物性値を高速に求めるための評価量について説明する。
運動方程式の境界条件は、ノイマン条件とディリクレ条件の２種類である。ここで、設定した物性値を係数とする運動方程式のノイマン条件（境界値＝実送受信データ）下の解をｈ１、ディリクレ条件（境界値＝実送受信データ）下の解をｄ１、設定した物性値を係数とする運動方程式のノイマン条件（境界値＝０；非一様性考慮）下の解をｈ２、ディリクレ条件（境界値＝０；非一様性考慮）下の解をｄ２とする。この場合、実送受信データを境界値とする非一様性を考慮した解は、ノイマン条件でｈ１＋ｄ１、ディリクレ条件でｈ２＋ｄ２として与えられる。非一様性を正確に推定できた場合、２つの条件は一致する。すなわち、ｈ１＋ｄ１＝ｈ２＋ｄ２が成立する。

そこで、本明細書においては、運動方程式の解の一致度により物性値の不均一性を検出する式として、非一様性を正確に推定できた場合に解がゼロになる次式の評価量Ｊを使用する。

Ｊ_ｄｈ＝（−ｈ１＋ｈ２）−（ｄ１−ｄ２） （１）式

この評価量Ｊ（Ｊ_ｄｈ）は、特許請求の範囲における「不均一検出式残差量」に対応する。

ここで、ラメ定数λの設定値と真値との差を表す量をσ、慣性密度ρの設定値と真値との差を表す量をμで表し、ノイマン条件に対する運動方程式のグリーン関数の行列表示（これは有限要素法などの手法と同様に構成される）ＧＮ［ｉ］［ｊ］と、ディリクレ条件に対する運動方程式のグリーン関数の行列表示ＧＤ［ｉ］［ｊ］とを用いると、（１）式の右辺第２の量（すなわち、ｄ１−ｄ２）は、行列（ＧＤ^−１−ＧＮ^−１）とベクトル（σ、μ）との積として表すことができる。

本実施例では、前述した関係により、σについて導き出される以下の関係式
ＧＮ（ｄε_ｉ）＊（ｈ１＋ｄ１）＝σ （２）式
を用いて物性値のずれ量ｄε_Ｉを計算し（ずれ量推定ステップ）、物性値の設定値ε_ｉをε_ｉ＋１＝ε_ｉ＋ｄε_ｉに逐次更新する（物性値更新ステップ）。ここでのずれ量ｄε_Ｉは、物性値の初期値ε_Ｉ０に対する値である。なお、（２）式の導出についての詳細は後述する。

この物性値の設定と更新（Ｓ２０４）を、設定された微少量ｅについて、ｉ番目に更新した物性値ε_i を変数とする評価量Ｊ_ｄｈ（ε_i ）が０とみなせる条件（すなわち、｜Ｊ_ｄｈ（ε_i ）｜＜ｅ）、あるいは評価量の指数関数ｅｘｐ（−｜Ｊ_ｄｈ（ε_i ）｜^２）が最大とみなせる条件（すなわち、｜ｅｘｐ（−｜Ｊ_ｄｈ（ε_i＋１ ）｜^２）−ｅｘｐ（−｜Ｊ_ｄｈ（ε_i ）｜^２）｜＜ｅ）が満たされるまで繰り返し（Ｓ２０３）、条件（収束条件）が満たされたときの物性値の設定値を推定値として物性値算出ステップの出力とする（Ｓ２０５）。なお、評価量Ｊは、送受信データにより算出される（Ｓ２０２）。

この評価量Ｊを用いれば、理想的な（近似していない）物理方程式を満たす解を、評価量Ｊの極値（この実施例ではゼロ）を与える物性値として求めることができる。しかも、求められる物性値からは、屈折や散乱による影響を除去することができる。さらに、短い計算時間で物性値を求めることができる。

次に、指数関数Ｐの最大化処理として、収束判定ステップ（Ｓ２０３）を実現する場合について説明する。拘束条件に事前確率を含ませる場合、収束判定ステップ（Ｓ２０３）は、図２−２の（ｂ）のように表すことができる。この判定処理は、事前確率算出ステップ（Ｓ２０３１）、条件つき確率算出ステップ（Ｓ２０３２）、事後確率算出ステップ（Ｓ２０３３）、事後確率最大化判定（Ｓ２３４）から構成される。

ここで、例えば計測領域内の物性値を“ｃ”、不連続部で値“１”を採り、連続部で値 “０”、を採る線過程を“lc”で表すと、空間的に一定値が連続する領域とその境界からなるという仮定を表す事前確率量Jcは、次式のように表すことができる。

Jc＝(∇c)^2(1−lc)＋lc （３）式

この事前確率量Jcが、特許請求の範囲における「条件式残差量」Jcに対応する。なお、当該事前確率量Jcについても、「不均一検出式残差量」Ｊの場合と同様に、指数関数Pc＝exp(−Jc^2)として表すことができる。このように、物性値の算出に事前確率量を導入す
ると、物性値の算出に必要な情報量を補うことができる。このため、計測点数が不足している場合にも計測領域の物性値を計算することができる。なお、事前確率量Jcで評価量Ｊを補助する場合の評価量は、J^2＋Jc^2又はP＊Pc（＝exp(−(J^2＋Jc^2))で与えられる。

ここで、評価量Ｊは受信データの１次式で与えられ、指数関数Ｐの指数（−J^2）は受信データの２次式で与えられる。一般に逆行列の計算に長時間を要する。しかし、各関数が１次又は２次式であることで、ステップＳ２０４における物性値の微小な変化に伴う行列（ＧＮ−ＧＤ）の微小な変化に対してSherman-Morrison公式やWoodbury公式その他の行列の摂動展開公式を適用でき、評価量Ｊの計算結果の更新に要する演算量を低減することができる。具体的には、逆行列の計算時間の大幅な低減により、評価量Ｊ全体の計算時間を短縮できる。

以上の説明では、物性値の計算を補助する条件式を反映した量（条件式残差量）を、事前確率量Jcを用いて表現する場合について説明した。しかし、従来の画像再構成方法によって計算された画像を補助的に用いた低次近似一致度量を定義し、前述した「条件式残差量」Jcを、事前確率量と低次近似一致度量の和と定義することもできる。

具体的には、例えば事前確率算出ステップＳ２０３１を、（１）事前確率量目算出ステップと、（２）低次近似一致度量計算ステップとで構成しても良い。ここで、低次近似一致度量計算ステップは、（１）指向性を合成する方法、投影写像その他の従来方法により物性値分布ｃ’を画像化する画像作成ステップと、（２）物性値更新ステップＳ２０４で仮定する物性値分布ｃ（画像）と従来方法による物性値分布ｃ’（画像）との差分（＝c-c’）の２次式(c-c’)^2である低次近似一致度量を計算する低次近似一致度量計算ステップとで構成し、「条件式残差量」Jcを、事前確率量と低次近似一致度量の和として算出しても良い。

Jc＝(∇c)^2(1−lc)＋lc＋(c−c’)^2 （４）式

以上においては、「不均一検査式残差量」Ｊと「条件式残差量」Jcで規定される２次式の線形和J^2＋Jc^2を最小とする物性値又はそれら２つの残差量で規定される指数関数P＊Pc＝exp(−(J^2＋Jc^2))を最大にする物性値を求める方法について説明したが、さらに運動方程式の右辺と左辺の差分の２次式で与えられる方程式残差量Jeを評価量に含めても良い。この場合には、２次式の線形和J^2＋Jc^2＋Je^2を最小とする又はその指数関数P＊Pc＊Pe＝exp(−（J^2＋Jc^2＋Je^2））を最大とする物性値を求めても良い。

本発明において、評価量の条件式残差量が、物性値分布の性質を表す事前確率である事前確率量及び指向性を合成する方法若しくは投影写像により受信信号をモデル化する方法による再構成画像との一致度を表す低次近似一致度量のうちの少なくとも１つ以上を含む場合には、事前確率の導入により情報量を補うことができる。このため、計測点数が不足する場合にも、計測領域の物性値を計算することができる。

また、前述した行列（ＧＮ−ＧＤ）はフルランク行列でない。このため、行列（ＧＮ−ＧＤ）×ベクトル（σ、μ）＝０の解である（σ＊、μ＊）が不定量となる。従って、前述した物性値更新ステップ（Ｓ２０４）においては、この不定量を確率的にサンプリングすることにより物性値を更新する例が考えられる。

確率的に物性値の更新値の候補（σ’、μ’）をサンプリングする場合には、例えば次式により、現在の物性値（σ、μ）を更新値に更新することができる。

P(σ,μ→σ’,μ’)＝min(1,P(σ,μ)/P(σ’,μ’）） （５）式

なお、P(σ,μ)はexp(−Jσ,μ^2−J）であり、P(σ’,μ’)はexp(−Jσ’,μ’^2−J）である。この更新規則はメトロポリスサンプリング法と呼ばれる。以上の処理内容の実行により、計算規模を小さくすることができる。

（具体例１）
（探触子の配置例と計測結果の表示例）
以下、前述した画像再構成装置を人体内部の計測に使用する場合について説明する。図３−１及び図３−２は、探触子２を用いた計測対象の測定態様を概念的に説明する図である。図３−１は、計測領域がヒトの頭蓋である場合を示し、図３−２の（ａ）は、計測領域がヒトの腹部である場合を示し、図３−２の（ｂ）は計測領域がヒトの乳房である場合を示している。

各図における「００１」〜「０１３」が不均一性の原因となる内部構造である。因みに、「００１」は大脳灰白質、「００２」は大脳白質、「００３」は頭蓋骨、「００４」は上皮、「００５」は脂肪、「００６」は上皮、「００７」は筋肉、「００８」は体腔、「００９」は臓器、「０１０」は上皮、「０１１」は脂肪、「０１２」は繊維腫、「０１３」は癌組織である。探触子２の空間上の配置位置は、図３−１及び図３−２の（ａ）のように計測領域の境界の一部でも良いし、図３−２の（ｂ）のように計測領域の境界全体でも良い。なお、図中の符号「０」は計測対象を示す。

図４−１は、物性値算出部５２の処理の概要を説明する図である。図４−１の（ａ）は、探触子２が計測対象０の境界の一部に接している状態を示す。前述したように、計測対象０の内部構造の評価に使用する評価量は、運動方程式を用いて計算する。運動方程式は、有限要素法や境界要素法と同様、計測対象の空間が微小なセルに分割されていると考え、方程式の変数をセルの節点に対して定義し、方程式の係数を各セルに対して定義する。

物性値算出部５２は、図４−１の（ｂ）に示すように、計測対象の空間内に数値計算用の計測領域４１１を設定する。図中、計測領域４１１は、探触子２の送受信面に対して一定の大きさを有する断面矩形形状の領域として与えられる。次に、物性値算出部５２は、計測領域４１１のうち実受信データが得られる境界４１２を設定する。図４−１の（ｂ）では、探触子２のうち送受信素子２１の配列範囲に境界４１２を設定する。さらに、物性値算出部５２は、計測領域４１１の内部を微小なセル４１３に仮想的に分割する。図では、セル４１３の水平方向の一辺長を送受信素子２１の長さに揃えている。また、図では、セル４１３の深さ方向の一辺長も水平方向の一辺長に揃えている。もっとも、各辺の長さは、計測能力に応じて定めることができる。図４−１の（ｃ）は、各セル４１３の評価量を極値にする物性値４２の違いを濃度の違いによってセル毎に表現した図である。

図４−２の（ａ）〜（ｃ）は、それぞれ図４−１の（ａ）〜（ｃ）に対応する。ただし、図４−２は、計測領域の境界全体を取り囲むように、送受信素子２１を円環状に配置する例を表している。図４−２の（ｂ）の場合、計測領域４１１は円形状に設定され、境界４１２は計測領域４１１を取り囲む正多角形形状として設定される。このように計測領域４１１が非矩形形状である場合、セル４１３は三角形状に設定することが好ましい。セル４１３の基本形状を三角形状とすることにより、計測領域４１１の内部に過不足無くセルを配置することができる。因みに、図４−２の場合には、物性値４２を４段階で表している。

図５に、物性値−表示量換算部５３による出力画面例を示す。図５の（ａ）及び（ｂ）は、図４−１の（ｃ）に対応する。ここでも、計測対象は生体である。評価量で評価される運動方程式が動弾性方程式である場合、物性値は、ラメ定数λ、ラメ定数μ、慣性密度ρ、減衰係数ｃとなる。表示量は、λ、μ、ρ、ｃと、これらの関数である体積弾性率、縦波速度、横波速度、ポアソン比、ヤング率、インピーダンスなどが考えられる。表示部７は、これらのうち１つ以上（好ましくは複数個）の表示量を同一画面上に同時に表示しても良い。なお、図５の（ａ）は縦波速度の表示例であり、図５の（ｂ）は減衰係数の表示例である。例えばこれら２つの画面を、表示部７の画面上に同時に表示することが望ましい。複数の表示量を同一画面上に表示することにより、インピーダンス（形状）が同じであるが体積弾性率（硬さ）が異なる部位である腫瘍など診断に有用な情報を容易に認識することができる。

（データフォーマットの例）
図６にデータフォーマットの一例を示す。なお、図６に示すデータフォーマットは、有限要素法の一般の教科書に記載されているフォーマットと同様のものである。図６のデータを使ったときの処理フローの詳細を図１２に示す。図６の（ａ）は図４−１（ｂ）で示した、撮像対象領域を仮想的に分割した微小なセルである。図６の（ａ１）に示すようにセルに番号を振り、図６の（ａ２）に示すようにセルの辺にも番号を振る。本例ではセルの数が９、辺の数が２４の例を示した。図６の（ａ２）中、２重線で示した辺は撮像領域の境界で、実送受信データ（音圧値、音圧の空間勾配値）に相当し、１重線の辺では値が不明である。物性値はセルごとに定義され、測定値は辺ごとに定義されるとする。

これらの設定に対して、物性値算出ステップＳ２００内では、図６の（ｂ）から（ｅ）のデータフォーマットを用いて物性値ε_Aiを算出する。（ｂ１）は物性値の真値ε_Ａ、（ｂ２）は物性値の仮定値ε_ｉ、（ｂ３）はi番目の仮定値とｉ＋１番目の仮定値の差分値ｄεであり、それぞれサイズ９（セルの数）のベクトルで表される。（ｃ１）は各辺での音圧の値ｐ（ディリクレ条件）、（ｃ２）は各辺での音圧の勾配の値ｋｄ_ｎｐ（ノイマン条件）である。（ｃ１）及び（ｃ２）とも、サイズ２４（辺の数）のベクトルで、図６（ａ２）の２重線で示した、境界上の辺、別の言い方をすると実送受信データが得られている辺（本例では１、３、６、８、９、１１、１６、１８、２０、２２、２３、２４番目の辺）のみで値を持ち、他の辺での値は０とする。

（ｄ１）はディリクレ条件に対する、運動方程式のグリーン関数を表す行列ＧＤ［ｉ］［ｊ］、（ｄ２）はノイマン条件に対する、運動方程式のグリーン関数を表す行列ＧＮ［ｉ］［ｊ］であり、サイズは２４（辺の数）行２４列、物性値の設定値の関数になっている。すなわち、物性値の仮定値を変えると行列の値が変わる。物性値と行列の関係は、一般の教科書を参照されたい。これらのデータを用いて、図２−１及び図２−２で述べた評価量を構成する。以下では、不均一検出式残差量を用いる例を説明する。具体的には、まず、観測対象の境界上で音圧および音圧の空間勾配の送受信データをそれぞれ（ｃ１）及び（ｃ２）に格納する。次に、物性値の仮定値（ｂ２）を設定して、（ｄ１）行列ＧＤ［ｉ］［ｊ］、（ｄ２）ＧＮ［ｉ］［ｊ］を計算する。次に、ｇを境界上の辺を１、その他の辺を０とする対角行列、Ｎｇを単位行列からｇを差し引いた行列、σを物性値の不均一性を外力に換算した仮想のベクトルとして、ディリクレ条件下の解（境界値＝実送受信データ）ｈ１はＧＤ^−１＊（ｐ＋Ｎｇ＊ｐ＊ｇ）、ディリクレ条件下の解（境界値＝０；非一様性考慮）ｄ１はＧＤ^−１＊（σ）、ノイマン条件下の解（境界値＝０；非一様性考慮）ｈ２はＧＮ^−１＊（ｋｄ_ｎｐ）、ノイマン条件下の解（境界値＝０；非一様性考慮）ｄ２はＧＮ^−１＊（σ）で表されるので、不均一検出式残差量Ｊ_ｄｈ＝（−ｈ１＋ｈ２）−（ｄ１−ｄ２）を−ＧＤ^−１＊（ｐ＋Ｎｇ＊ｐ＊ｇ）＋ＧＮ^−１＊（ｋｄ_ｎｐ）−（ＧＤ^−１−ＧＮ^−１）＊（σ）、により計算する。

このとき、σはＪ_ｄｈ＝０から値を計算し、（２）式に示す関係式（すなわち、ＧＮ（ｄε_ｉ）＊（ｈ１＋ｄ１）＝σ）から物性値のずれ量ｄε_ｉを計算して、物性値の設定値ε_ｉをε_ｉ＋１＝ε_ｉ＋ｄε_ｉに更新する。

なお、図６の送受波データｃ１、ｃ２は、時系列データを周波数解析した結果の、ある周波数成分を表しており、周波数毎に解析した結果から、周波数成分ごとの物性値の分布を表示してもよいし、周波数成分ごとの物性値を送受信の有効な帯域内で積算した物性値の分布を表示に使ってもよい。また、ここまで記述を行った式を、時系列データに対応するように適宜変形を行い、過渡的な物性値分布を表示することも、勿論本発明の思想の範囲内である。

（評価量と物性値の求め方の関係）
図７に、評価量と物性値の求め方の関係を説明する概念図を示す。図７の（ａ１）は方程式残差量Ｊ_ｅｑを評価量とし、評価量を０にする物性値を真値としてＳ２００の出力とする場合、（ａ２）は方程式残差量Ｊ_ｅｑを評価量とし、評価量を最大にする物性値を真値としてＳ２００の出力とする場合、（ｂ）は不均一検出式残差量Ｊ_ｈｄを評価量とし、評価量を最大にする物性値を真値とする場合、（ｃ）は方程式残差量Ｊ_ｅｑと条件式残差量Ｊ_ｃの積を評価量とし、評価量を最大にする物性値を真値とする場合を示す。方程式残差量Ｊ_ｅｑを評価量とする（ａ１）及び（ａ２）の場合、物性値εの設定値ε_ｉ（ｉ＝１、２、…）を従来行われている方法でサンプルして更新し、物性値の真値をε_Ａに収束させる。従来のサンプル方法は、ランダムにεを更新するか、現設定値ε_ｉにおける評価量の局所的な傾きｄＪ／ｄεに比例した差分ｄε_ｉ＝ｃ＊ｄＪ／ｄε（ε＝ε_Ｉ）を用いて次の設定値ｄε_ｉ＋１＝ε_ｉ＋ｄε_ｉを決める。不均一検出式残差量Ｊ_ｈｄを評価量とする（ｂ）の場合、真値の推定値ε_Ａ と現在の設定値との差分を計算して次の設定値ｄε_ｉ＋１＝ε_ｉ＋ｄε_ｉを決める。真値と設定値との差（ずれ量）を推定して物性値を更新する（ｂ）の方法は、真値の情報は使わず設定値の情報のみを用いる（ａ）の方法より、早く真値に収束させることができる。方程式残差量Ｊ_ｅｑと条件式残差量Ｊ_ｃの積を評価量とする（ｃ）の場合は、条件式を課すことにより、図７の（ａ２）及び（ｂ）より幅の狭い曲線で示したように、条件式を満たさない物性値が評価範囲外になり、評価範囲を狭めることで、測定値が十分でない場合に物理的に不適切な解に収束することを防ぎ、また、図７の（ａ２）よりも早く真値に収束させることができる。

（送信シーケンスの送信例）
図８に、送受信データ取得部５１が出力する送信シーケンスの基本シーケンスの一例を示す。この実施例の場合、送受信素子２１は、送信シーケンスを２回以上送受信する。図８の（ａ）は１回目の送受信例であり、図８の（ｂ）は２回目の送受信例を示している。図中、「２１」は送信シーケンスを出力する送信素子、「２２」は送波波形、「２３」は波数ベクトルを示している。本例では、１回目と２回目のそれぞれにおいて、複数の送受信素子２１からバースト長の等しい信号を同時に送出する。これにより、位相面が揃った平面バースト波２２が生成される。

図８の場合、１回目の送信シーケンスと２回目の送信シーケンスではキャリア周波数を変えている。ここでのキャリア周波数は、好ましくは0.5MHz以上10MHzの範囲で可変する。また、送受信素子２１の空間配置をキャリア周波数に逆比例する間隔に設定する。さらに、１回目と２回目において、対応する波数ベクトル２３が同様の形状になるように指向性を合成する。すなわち、キャリア周波数の低い１回目の送波のメインローブ２３ＬＭとグレーティングローブ２３ＬＧ１、２３ＬＧ２のそれぞれと、キャリア周波数が高い２回目の送波のメインローブ２３ＨＭとグレーティングローブ２３ＨＧ１、２３ＨＧ２とが同じ角度方向になるように指向性が合成される。

なお、図８の（ａ）の２１Ｌ１（網掛けで示す）は素子配置の重心を示したもので、１回目の送波に使用する送受信素子２１の個数と２回目の送波に使用する送受信素子２１の個数が等しい必要はない。因みに、１回目の送信のキャリア周波数がと２回目の送信のキャリア周波数が２倍異なる場合、１回目の送信においては２個の送受信素子２１からなるサブアレイを前述した重心位置に配置しても良い。

以上説明した送信波形及び送受信素子２１の空間配置を有する基本シーケンスの採用により、キャリア周波数のみが異なる類似形状、かつ、平易な形状の音場を生成することができる。この結果、評価量を用いたラメ定数と慣性密度の比である音速の同定精度を向上できる。また、キャリア周波数の違いを利用することにより、慣性質量の同定精度を向上できる。また、本例のように平面バースト波を利用することにより、１つの送受信素子２１だけから信号を送信する場合に比べて測定ノイズに起因するＳ／Ｎを向上させることができる。

図９及び図１０に、送受信データ取得部５１による送信シーケンスの出力例を示す。ｚｙ７は探触子２が計測対象０の境界の一部に接している場合を示し、図９は探触子２が計測対象０の境界全体に接している場合を示す。

図９の（ａ１）と（ａ２）は、それぞれ送信素子Ｔと受信素子Ｒの配置例を示している。送信素子Ｔはグレーティング（grating）の生成が問題にならない程度に疎であり、平面波が形成できる程度に密な配置である。一方、受信素子Ｒの空間配置は計測領域に対する見込み角を最大にする配置である。受信素子Ｒの空間配置は計測領域に対する見込み角を最大にする配置である結果、物性値の計算精度が向上する。

図９の（ｂ１）〜（ｂ３）は、図８で説明した基本シーケンスにおいて、計測領域内の波数ベクトルが可能な限り異なる方向に横切るように複数回繰り返す送波シーケンスを示している。受信素子Ｒの空間配置は、計測領域に対する見込み角を最大とする配置とする結果、評価量に含まれる行列の固有空間のうち、装置構成の制約範囲内で最大限独立な成分に相当する情報を入手可能になり、計測点数が不足している場合の係数同定精度が向上する。

図９のように探触子２が計測対象０の境界の一部に接している場合は、送信素子Ｔの配置位置を固定し、各送信素子Ｔから信号が送信される時刻に位相差を設けることにより、波数ベクトルの方向を任意に可変できる。

一方、図１０に示すように、探触子２が計測対象０の境界全体に接している場合は、送信素子Ｔの配置位置を変化させることにより、波数ベクトルの方向を変化させることができる。

例えば図９の例は、探触子２が計測対象０の境界の一部に接している場合において、６回の送受信を行う例を示している。偶数回目の送受信時で使用する送信素子Ｔの素子間隔は図９の（ａ１）に示すようにある程度密に配置する。一方、奇数回目の送受信時で使用する送信素子Ｔの素子間隔は、偶数回目の送受信時における素子間隔の２倍に設定する。受信素子Ｒは１回目の送受信時から６回目の送受信時まで図９の（ａ２）に示すように固定である。

１回目の送受信時には図９の（ｂ１）の方向に低周波バースト平面波を送信し、２回目の送受信時には図９の（ｂ１）の方向に高周波バースト平面波を送信し、３回目の送受信時には図９の（ｂ２）の方向に低周波バースト平面波を送信し、４回目の送受信時には図９の（ｂ２）の方向に高周波バースト平面波を送信し、５回目の送受信時には図９の（ｂ３）の方向に低周波バースト平面波を送信し、６回目の送受信時には図９の（ｂ３）の方向に高周波バースト平面波を送信する。勿論、各バースト平面波の受信波を受信素子Ｒがそれぞれ受信する。そして、これら６本の運動方程式を連立して評価量を作る。

図１０の例は、探触子２が計測対象０の境界全体に接している場合において、８回の送受信を行う例を示している。偶数回目の送受信で使用する送信素子Ｔの配置範囲が偏在するように、かつ、素子間隔がある程度密になるように配置する。一方、奇数回目の送受信時で使用する送信素子Ｔの素子間隔は、偶数回目の送受信時における素子間隔の２倍に設定する。受信素子Ｒは１回目の送受信時から８回目の送受信時まで図１０の（ａ２）に示すように均等かつ固定的に配置する。

１回目の送受信時には図中上側に素子間隔が２ｄとなるように配置した送信素子Ｔから図１０の（ｂ１）の方向に低周波バースト平面波を送信し、２回目の送受信時には図中上側に素子間隔がｄとなるように配置した送信素子Ｔから図１０の（ｂ１）の方向に高周波バースト平面波を送信し、３回目の送受信時には図中右側に素子間隔が２ｄとなるように配置した送信素子Ｔから図１０の（ｂ２）の方向に低周波バースト平面波を送信し、４回目の送受信時には図中右側に素子間隔がｄとなるように配置した送信素子Ｔから図１０の（ｂ２）の方向に高周波バースト平面波を送信し、５回目の送受信時には図中下側に素子間隔が２ｄとなるように配置した送信素子Ｔから図１０の（ｂ３）の方向に低周波バースト平面波を送信し、６回目の送受信時には図中下側に素子間隔がｄとなるように配置した送信素子Ｔから図１０の（ｂ３）の方向に高周波バースト平面波を送信し、７回目の送受信時には図中左側に素子間隔が２ｄとなるように配置した送信素子Ｔから図１０の（ｂ４）の方向に低周波バースト平面波を送信し、８回目の送受信時には図中左側に素子間隔がｄとなるように配置した送信素子から図１０の（ｂ４）の方向に高周波バースト平面波を送信する。勿論、各バースト平面波の受信波を受信素子Ｒがそれぞれ受信する。そして、これら８本の運動方程式を連立して評価量を作る。

図１１に、前記送受信データ取得部５１が出力する送信シーケンスの他の例を示す。図９及び図１０の場合には、個々の送受信回に関する限り、送信素子Ｔの配置を切り替え、バースト平面波の送信方向が特定の一方向に向かうように設定した。しかし、計測領域に形成される音場が対称性を持つという基準で送受信シーケンスを設定する場合、図１１の（ａ１）及び（ａ２）に示すように送受信素子２１が円周に沿って配置できるのであれば、図１１の（ｂ）に示すように、円筒波を生成することも考えられる。

この場合、送信素子Ｔと受信素子Ｒを、図１１の（ａ１）及び（ａ２）に示すように円周上に等間隔に配置された同一の送受信素子２１を切り替えて用いても良い。この場合、送信素子Ｔと受信素子Ｒとで使用する送受信素子２１の切り替えをせずに、同一素子を送受両方に用いることができる。このため、素子切り替えスイッチからノイズが混入されるのを防ぐことができ、物性値の計算精度を向上させることができる。

なお、図９では送受信データを６回、図１０では送受信データを８回送受信する例を説明したが、評価量Ｊは、Ｎ＝（素子数）＊（送受信時間）／（時間サンプル間隔）とするＮ行Ｎ列の行列で表されるため基底数はＮとなる。そのため、Ｎ回分の送受信データと、Ｎ回以上分の送受信データの情報量は同じである。一方、撮像装置において撮像時間は短いことが要請される。再構成画像の精度と撮像時間を両立させるためには、画像の精度に直結する情報量を増やさず、撮像時間を長くするＮ＋１回以上の送受信は避けるべきであり、送受信回数はＮ＝（素子数）＊（送受信時間）／（時間サンプル間隔）以下とすることが好ましい。

１：探触子、２：本体、３：送信部、４：受信部、５：信号処理部、５１：送受信データ取得部、５２：物性値算出部、５３：物性値−表示量換算部、６：メモリ、７：表示部、８：入力手段、９：制御部。

Claims (12)

1. 計測領域との間で信号を送受信し、計測領域の物性を、数値計算を通じて定量する画像再構成方法であって、
   計測領域との間で信号を送受信する送受信データ取得ステップと、
   計測領域から受信された実受信信号から物性値を算出する物性値算出ステップと、
   物性値を表示値に換算する物性値−表示値換算ステップとを有し、
   前記物性値算出ステップは、
   前記送受信データにおいて取得された送受信データと、仮定された物性値の分布と波の伝播を考慮した方程式から導出した、前記仮定された物性値のずれ量の推定を計算するステップを有し、前記仮定された物性値と前記ずれ量の和を算出された物性値として出力し、
   さらに、前記物性値算出ステップは、
   前記方程式から導出した、物性値が真値ε _Ａ の場合に０又は極値になる評価量Ｊが、物性値の設定値について収束条件を満たすか否か判定し、前記評価量Ｊが収束条件を満たさない場合には、物性値の真値ε _Ａ に対する設定値のずれ量ｄεの推定値の更新を通じて物性値の設定値を更新して前記判定を再実行し、前記評価量Ｊが収束条件を満たす場合には更新後の物性値を算出結果として出力する
   ことを特徴とする画像再構成方法。
2. 請求項１に記載の画像再構成方法において、波の伝播を考慮した方程式から前記ずれ量を推定するステップは、前記計測領域を微小な空間セルに分割して計算することを特徴する画像再構成方法。
3. 請求項１に記載の画像再構成方法において、
   前記物性値算出ステップは、
   ｉを１以上ｎ−１（ｎは２以上の整数であり、物性値ε_ｉの設定回数の最大値）の整数とし、ｉ回目の物性値の設定値ε_ｉと真値ε_Ａの差分であるずれ量ｄε_ｉを推定するずれ量推定ステップと、
   ｉ回目の物性値の設定値がε_ｉ、推定されたずれ量がｄε_ｉの場合に、ｉ＋１回目の物性値の設定値ε_ｉ＋１をε_ｉ＋ｄε_ｉで与える物性値更新ステップと、
   ｅを予め設定した微小量とする場合に、物性値ε_ｉを変数とする評価量Ｊが、｜Ｊ（ε_ｉ）｜＜ｅ又は｜Ｊ（ε_ｉ＋１）−Ｊ（ε_ｉ）｜＜ｅを満たすとき、物性値ε_ｉが真値ε_ｉに収束したと判定する収束判定ステップとを有する
   ことを特徴とする画像再構成方法。
4. 請求項３に記載の画像再構成方法は、
   送受信シーケンスは、ずれ量ｄε_ｊの推定への寄与が大きい送受信シーケンスである
   ことを特徴とする画像再構成方法。
5. 請求項１に記載の画像再構成方法は、
   運動方程式の演算子項と外力項の差分の残差である方程式残差量と、
   ２種類の境界条件による運動方程式の解の一致度により物性値の不均一性を検出する式の残差である不均一検出式残差量と、
   制約条件の残差である条件式残差量のうちの１つ又は任意の複数を含む線形和若しくはそれらの積若しくはそれらを指数に含む指数関数を前記評価量とする
   ことを特徴とする画像再構成方法。
6. 請求項１に記載の画像再構成方法であって、
   前記送受信データ取得ステップは、
   送受信素子の送受信回数と、送受信素子の空間配置と、送受信波形とを決定して送受信結果を取得する送受信シーケンス決定ステップを有し、
   前記送受信シーケンス決定ステップは、送受信素子の送受信回数が２回以上であり、送信波形は0.5MHz以上10MHz以下のキャリア周波数が異なるバースト長の等しい複数の平面バースト波を少なくとも含み、送信素子の空間配置は前記キャリア周波数に逆比例する間隔であり、受信素子の空間配置は計測領域に対する見込み角を最大とする配置である送受信シーケンスを出力する
   ことを特徴とする画像再構成方法。
7. 請求項６に記載の画像再構成方法であって、
   前記送受信シーケンス決定ステップは、送受信回数は２回以上であり、送受信時刻が早いほど送受信の周波が低い送受信シーケンスを出力する
   ことを特徴とする画像再構成方法。
8. 請求項５に記載の画像再構成方法であって、
   前記物性値算出ステップは、
   境界値が実送受信データの場合のノイマン境界条件による解とディリクレ境界条件による解の差分である差分ベクトルと、
   境界値がゼロの場合のノイマン境界条件に対する運動方程式のグリーン関数とディリクレ境界条件に対する運動方程式のグリーン関数の差分である演算子行列と物性値の不均一性を表すベクトルの積と
   の差分の残差で与えられる前記不均一検出式残差量を線形和の１項あるいは指数関数積の１因数とする量を極値にする物性値を出力する
   ことを特徴とする画像再構成方法。
9. 請求項１に記載の画像再構成方法であって、
   前記物性値算出ステップは、
   評価量物性値分布の性質を表す事前確率である事前確率項と、指向性を合成する方法若しくは投影写像により再構成された画像との一致度を表す低次近似一致度量の少なくとも一方を含む前記評価量のうち条件式残差量を、線形和の１項あるいは指数関数積の１因数とする量を極値とする物性値を出力する
   ことを特徴とする画像再構成方法。
10. 請求項１に記載の画像再構成方法であって、
    前記物性値算出ステップは、
    前記評価量を極値にする物性値を、確率的手法の１つであるメトロポリスサンプリングを用いて算出する
    ことを特徴とする画像再構成方法。
11. 請求項５に記載の画像再構成方法であって、
    前記方程式残差量は、動弾性方程式の周波数空間表示、動弾性方程式の時間応答、動弾性方程式のヘルムホルツポテンシャル表示のうちのいずれかの方程式の残差であり、
    前記物性値は、ラメ定数λ、ラメ定数μ、慣性密度ρ、減衰係数ｃ、ラメ定数と慣性密度の比、ラメ定数と減衰係数の比、外力ｆのうちの１つ以上であり、
    前記表示量は、前記λ、前記μ、前記ρ、前記ｃ又はこれらの関数である体積弾性率、縦波速度、横波速度、ポアソン比、ヤング率、インピーダンス、外力のうちの１つ以上である
    ことを特徴とする画像再構成方法。
12. 計測領域との間で信号を送受信し、計測領域の物性を定量する計測装置であって、
    計測領域との間で信号を送受信する送受信データ取得部と、
    計測領域から受信した実受信信号から物性値を算出する物性値算出部と、
    物性値を表示値に換算する物性値−表示値換算部とを有し、
    前記物性値算出部は、
    前記送受信データにおいて取得された送受信データと、仮定された物性値の分布と波の伝播を考慮した方程式から導出した、前記仮定された物性値のずれ量の推定を計算する処理部を有し、前記仮定された物性値と前記ずれ量の和を算出された物性値として出力し、
    さらに、前記物性値算出部は、
    前記方程式から導出した、物性値が真値ε _Ａ の場合に０又は極値になる評価量Ｊが、物性値の設定値について収束条件を満たすか否か判定し、前記評価量Ｊが収束条件を満たさない場合には、物性値の真値ε _Ａ に対する設定値のずれ量ｄεの推定値の更新を通じて物性値の設定値を更新して前記判定を再実行し、前記評価量Ｊが収束条件を満たす場合には更新後の物性値を算出結果として出力する
    ことを特徴とする計測装置。

Images