JP2022016414A - 超音波減衰パラメータを推定するための方法及びシステム - Google Patents

Abstract

【課題】媒体中の領域における超音波減衰パラメータを推定するための方法を提供する。
【解決手段】本方法は、少なくとも１つの超音波トランスデューサ（２）に関連付けられた処理装置（８）によって実施される。本方法は、（ａ）媒体において少なくとも１つのパルスがトランスデューサによって送信される送信ステップ（１０１）と、（ｂ）パルスに応答してデータがトランスデューサによって取得される受信ステップ（１０２）と、（ｃ）処理装置によってデータが処理されて領域の後方散乱された取得データを提供する処理ステップ（１０３）と、（ｄ）時空間領域における深度の関数である、後方散乱された取得データの自己相関関数が決定される関数決定ステップ（１０４）と、（ｅ）自己相関関数に基づいて超音波減衰パラメータが推定される減衰推定ステップ（１０５）とを含む。本発明は、対応するシステムにも関連する。
【選択図】図３

Description

本発明は、特に医療用イメージング（画像生成）のための、イメージング方法及び上記方法を実現する装置に関する。

本開示は、特に、媒体中の領域における超音波減衰パラメータを推定するための方法に関し、より正確には、媒体の内部の領域における局所的な超音波減衰パラメータを推定して画像生成するための方法に関する。本方法は、上記媒体に対して所定関係を有する複数のトランスデューサ（例えば、ライン又はアレイ）に関連付けられた処理装置によって実現されてもよい。

古典的な超音波イメージングは、所与の点に合焦する円筒波を用いて媒体に対して超音波照射（insonification）を行うことからなる。この単一の超音波照射の後方散乱エコーを用いて、画像のライン全体が動的受信ビームフォーミング処理を用いて計算される。画像全体を構築するために、この手順は、所与の深度において横方向のラインに沿って走査する、合焦された一組の波（焦点面と呼ぶ）を送信することで繰り返される。合焦された各波について、動的ビームフォーミングが実行され、画像全体が１ラインずつ構築されて取得される。動的ビームフォーミングは、受信モードでは一様な合焦を保証するが、送信モードでは、焦点は所与の深度に固定される。最終的な画像は、焦点面において、また、焦点軸長に対応する媒体の限られた領域において最適である。しかしながら、回析の法則によって課されるこのエリアの外部では、画像品質は、（合焦されたビームの近傍界及び遠方界における）他の深度では急速に劣化する。

この制限を克服するために、古典的な解決方法は、多焦点イメージングを実行するというものであり、画像全体にわたって均質な品質を得るために、異なる送信焦点深度が使用される。所与の焦点深度における各送信は、軸焦点距離によって定められた範囲を有する領域における部分画像の実施を可能にする最終的な画像は、様々な深度に対応するこれらの部分画像を再結合したものを用いて取得される。最適な多焦点画像は、典型的には、数十個の焦点面を必要とする。このことは、典型的には１０フレーム／秒未満のフレームレート制限をもたらし、この制限は超音波イメージングでは許容できない。画像品質とフレームレートとの間での良好な妥協点は、約４つの焦点深度画像である。

合成動的送信合焦を実行することで、画像品質の改善を期待することができる。そのようなアプローチは、互いに異なる一組の超音波照射のビームフォーミングを行って、次いで結合することにより、動的送信合焦（すなわち、画像中の画素数と同じ個数の焦点深度）を再合成することにある。

さらに、本願出願人によって出願された特許文献１によれば、改善された合成超音波イメージング方法が知られ、これは、例えば特許文献２に開示される従来技術の平面波合成超音波イメージング方法を改善することを可能にする。

特許文献１は、少なくとも以下のステップを含む超音波イメージングのための方法を提案する。
ａ）画像生成される領域に対して複数の超音波が送信され、各超音波に応答してトランスデューサのアレイによって各組の生データが取得される送信ステップであって、上記複数の超音波が互いに異なる空間周波数内容を有し、領域における画像生成される複数の場所のそれぞれについて、各組の生データが、対応する超音波に応答してトランスデューサによって受信された時間信号を表す送信ステップ。
ｂ）画像生成される領域における複数の仮想送信焦点ゾーンのそれぞれについて、複数組の生データから少なくとも一組のコヒーレントデータが合成されるコヒーレンス向上ステップ。
ｃ）仮想送信焦点ゾーンのそれぞれに含まれる複数の場所のそれぞれについて、一組のコヒーレントデータを用いてビームフォーミングによって画像の画素が計算されるビームフォーミングステップ。

これらの性質のおかげで、生データの空間的コヒーレンスは、ビームフォーミング前にステップｂ）において回復され、それにより、様々な超音波の送信から受信されたデータを正確に結合することを可能にする。空間的コヒーレンスを回復する必要性は、画像生成される領域を空間的に広く拡散された波動場で照射する場合、媒体から戻るエコーを、画像生成される領域においてランダムに分散されたインコヒーレントなソース（散乱器）から発生された波動場とみなすことができるという事実に起因し、したがって、波動場の空間的コヒーレンスは生データでは失われる（又は非常に低くなる）。次いで、ビームフォーミングは、コヒーレンス回復ステップから得られるコヒーレントデータに対して実行されてもよく、その結果、より正確な画像をもたらす。

合成ビームフォーミングも呼ばれるこの技術に基づいて、特許文献３は、低減されたスペックルノイズを有して、媒体の内部の領域の画像を生成するもう１つのイメージング方法を開示する。このため、本方法は以下のステップを含む。
（ａ）媒体の内部においてトランスデューサによって第１の複数の波が送信される送信ステップ。
（ｂ）波に応答して上記トランスデューサによって一組のデータが取得される受信ステップ。
（ｃ）第２の複数のビームフォーミング処理によって一組のデータが処理されて、画像の少なくとも一部のビームフォーミングされた画素値を提供するビームフォーミングステップであって、各ビームフォーミング処理は、送信重み付けベクトルで生成された波に対応する一組のデータを用いるか、又は、ビームフォーミングされた画素値の計算において送信重み付けベクトルを用いるビームフォーミングステップ。
（ｄ）上記第２の複数のものに係るビームフォーミングされた画素値が互いに結合されて、画像の内部の画素の画素値を提供する結合ステップであって、送信重み付けベクトルが互いに異なりかつ直交する結合ステップ。

これらの特徴のおかげで、各送信重み付けベクトルは、無相関のスペックルノイズを生成し、重み付けられたデータの結合は、低減されたスペックルノイズを有する領域の画像を計算することを可能にする。

欧州特許出願公開第２１０１１９１号公報 米国特許第６５５１２４６号明細書 国際公開第２０１７／０９３７７８号

LYMBERIS ET AL: "Estimation of frequency-dependent attenuation based on parametric spectral analysis and correlation lags of the demodulated echo signal", ULTRASONIC IMAGING, DYNAMEDIA INC., SILVER SPRING, MD, US, vol. 13, no. 1, 1 January 1991 (1991-01-01), pages 1-26, XP026409816, ISSN: 0161-7346. HYUNGSUK KIM ET AL: "Attenuation estimation using spectral cross-correlation", IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS AND FREQUENCY CONTROL, IEEE, US, vol. 54, no. 3, 1 March 2007 (2007-03-01), pages 510-519, XP011175817, ISSN: 0885-3010, DOI: 10.1109/TUFFC.2007.274.

しかしながら、超音波イメージングにおいて望ましくは考慮されるべき別の現象は、検査される媒体の内部における超音波減衰である。超音波減衰は、送信された超音波への応答に対して直接的に影響する。それによって、減衰は、敏感な、周波数及び深度に依存する現象を構成する。したがって、組織減衰を考慮して例えば時間利得補償によって従来行われるように、結果的に得られる計算された画像に対する減衰のいかなる影響も補償することが望ましい。

一方、超音波減衰の信頼できる推定は、他の超音波診断目的で使用されてもよい。例えば、検査された媒体について、例えばヒトの肝臓又は筋肉について推定された大域的減衰パラメータは、その脂肪含有量を決定するために使用されてもよい。さらに、例えば画像を計算するために処理された、局所的な減衰の分布は、例えばヒトの乳房における、癌を検出するために使用されてもよい。

しかしながら、超音波減衰を推定するための従来の方法は、通常、周波数（又はスペクトル）分析を使用する。しかしながら、そのような分析はいくつかの欠点を意味する。

減衰が、敏感な、周波数及び深度に依存する現象であるので、正確な評価には、十分に大きな空間及び周波数の観測範囲を必要とする。このことは、大きな周波数及び空間的範囲の使用が望まれるかもしれないことを意味する。しかしながら、減衰を理由として、空間的範囲が大きくなり、また、周波数内容の変動が大きくなると、バイアスが大きくなる。一方、空間的範囲が小さくなると、スペクトル解析が不十分になる。さらに、空間的範囲が大きくなると、回析の影響が大きくなる。さらに、スペックルノイズは、減衰推定の変動を増大させる。

従って、この必要なトレードオフを考慮すると、周波数（又はスペクトル）分析によって超音波減衰を推定することは、不正確な結果をもたらし、それと同時に、特に、必要とされるフーリエ変換に起因して、非常に大きな計算上のコストを必要とする。結果として、可能な結果の変動が増大し、従って、推定の信頼性が低下する、すなわち再現可能性が低下するという問題がある。従って、取得時間と精度のレベルとの間のトレードオフが考慮されなければならない。Ａ．Ｌｙｍｂｅｒｉｓ他は、周波数に依存する減衰を推定する方法を提案する。特に、提案された方法２において、平均周波数推定量は、復調信号の利用可能なすべての自己相関遅延（ＡＣｎ）から導出された。非特許文献１を参照。さらに、Ｈｙｕｎｇｓｕｋ他は、スペクトル相互相関を用いた減衰推定方法を開示する。非特許文献２を参照。

現在、前述の問題を克服し、特に、媒体中の領域における超音波減衰パラメータを高い信頼性で推定するための方法及びシステムであって、より高速であり、計算上のコストがより低く、したがって、必要な処理パワーがより小さいという利点を有する方法及びシステムを提供することが引き続き望まれている。

従って、本開示の実施形態によれば、媒体中の領域における超音波減衰パラメータを推定するための方法が提供される。上記方法は、少なくとも１つの超音波トランスデューサ（これは上記媒体と所定関係にあってもよい）に関連付けられた処理装置によって実施される。本方法は、
（ａ）上記媒体において少なくとも１つのパルスがトランスデューサによって送信される送信ステップと、
（ｂ）上記パルスに応答してデータがトランスデューサによって取得される受信ステップと、
（ｃ）上記処理装置によって上記データが処理されて上記領域の後方散乱された取得データを提供する処理ステップと、
（ｄ）時空間領域における深度の関数である、上記後方散乱された取得データの自己相関関数がゼロ（０）遅延で決定される関数決定ステップと、
（ｅ）上記自己相関関数に基づいて超音波減衰パラメータが推定される減衰推定ステップとを含む。

上記方法は、例えば、検査された領域全体に係る、例えば超音波画像全体に係る１つの減衰パラメータの大域的減衰推定方法と呼ばれてもよい。

そのような方法を提供することによって、周波数（又はスペクトル）領域の代わりに時空間領域において分析を実行することで超音波減衰パラメータを推定することが可能になる。従って、推定は、さらにずっと正確になるという利点があり、より少ない計算コストをともない、これにより、例えばリアルタイムの、又は少なくとも準リアルタイムの計算モードを改善する。さらに、増大した精度に起因して、変動の減少を達成することができ、したがって、増大した再現性を達成することができる。

処理ステップ（ｃ）は、データがビームフォーミング処理によって処理されて領域のビームフォーミングされた取得データを提供するビームフォーミングステップを含んでもよい。上記ビームフォーミングされた取得データ又は後方散乱された取得データは、例えば画素値又はボクセル値の形式で提示されてもよい。

ビームフォーミング処理に起因して、取得データにおける回析パターンを低減することが可能になる。

ビームフォーミング処理は、例えば、合成ビームフォーミング処理であってもよい。これは、回析パターンをさらに低減することを可能にする。

さらに、スペックルノイズをさらに低減するようにステップ（ｃ）のビームフォーミング処理を構成することが可能である。例えば、ステップ（ｃ）は、
（ｃ１）複数のビームフォーミング処理によってデータが処理されて、少なくとも領域のビームフォーミングされた取得データを提供するビームフォーミングステップであって、各ビームフォーミング処理は、送信重み付けベクトルで生成されたパルスに対応するデータを用いるか、又は、ビームフォーミングされた取得データの計算において送信重み付けベクトルを用いるビームフォーミングステップと、
（ｃ２）上記複数のものに係るビームフォーミングされた取得データが互いに結合されて、生成された超音波画像の内部の画素の画素値を提供する結合ステップであって、送信重み付けベクトルは互いに異なりかつ直交する結合ステップとを含む。

しかしながら、後方散乱された取得データを提供するために、上述したビームフォーミング処理の代替処理が使用されてもよい。例えば、音響レンズを介して送信パルス（及びオプションで受信信号）を音響的にビーム整形することが使用されてもよく、又は、機械的な掃引／走査が使用されてもよい。そのような場合、処理ステップ（ｃ）においてビームフォーミング処理は省略されてもよい。

減衰推定ステップ（ｅ）において、超音波減衰パラメータは、深度方向の所定区間Δｚにわたって推定されてもよい。例えば、上記区間Δｚは、媒体における指定された第１の深度及び第２の深度の間で予め定義されてもよい。媒体における深度は、プローブの表面と関心対象点（空間的場所）との間の距離として定義されてもよいことに注意すべきである。

減衰決定ステップにおいて、超音波減衰パラメータは、例えば区間Δｚに係る、自己相関関数を解くことで決定されてもよい。

後方散乱された取得データの自己相関関数は、深度及び時間の関数であってもよい。

関数決定ステップ（ｄ）において、後方散乱された取得データにおける残留回折パターンは、予め決められた回析基準データに基づいて補償されてもよい。

後方散乱された取得データにおける残留回析パターンは、予め決められた基準超音波減衰パラメータを有する基準媒体に対してステップ（ａ’）～（ｅ’）を含む方法を適用し、自己相関関数を比較することで補償される。

言いかえれば、推定方法は、予め決められた、すなわち既知の超音波減衰パラメータを用いて、そのような基準媒体、例えばファントムに適用されてもよい。従って、決定された自己相関関数は基準自己相関関数と呼ばれてもよく、この自己相関関数を解くことで、後方散乱された取得データにおける残留回析パターンを決定してもよく、及び／又は、検査される関心対象領域を含む媒体の自己相関関数においてそれを補償してもよい。

自己相関関数を比較することは、両方の関数の間のデルタ関数を計算してデルタ関数を解き、媒体の超音波減衰パラメータを決定することを含んでもよい。

例えば、基準超音波減衰パラメータは、検査される媒体における区間Δｚに類似した基準媒体における区間Δｚ’に係るパラメータであってもよい。

予め決められた回析基準データは、基準媒体における異なる深度に係る予め決められたパラメータを含んでもよい。

パルスは、広帯域パルス及び／又は複数のパルスであってもよい。例えば、１つ又は複数の広帯域パルスを使用することは、空間的な限定、すなわちより正確な区間決定を可能にするので、推定の精度を増大させうる。

本開示はさらに、媒体中の領域における局所的超音波減衰パラメータを推定するための方法に関する。上記方法は上述した大域的減衰推定方法（すなわちステップ（ａ）～（ｅ））を含んでもよい。媒体内の深度方向ｚにおける複数の軸線（例えば、ステップ（ａ）～（ｃ）において走査されたライン）のそれぞれについて、
・上記自己相関関数を、駆動入力を有するノイズありの状態空間モデルとしてモデリングすることと、
・上記自己相関関数を、チホノフ正則化法、リッジ回帰法、及び／又はベイジアン回帰法を含む正則化法によって正則化することと、
・サビツキー・ゴーレイ（Savitzky-Golay）フィルタを用いて上記自己相関関数をフィルタリングすることと、モンテカルロ法と、
特異スペクトル分析と、
圧縮された検出及び／又は疎サンプリング法と、
のうちの少なくとも１つによって、Ｅｙを自己相関関数としかつｚを深度とするとき、数式Ｆ（Ｅｙ（ｚ），Ｂ）＝０における局所的超音波減衰Ｂを推定することを含む。

上記方法は、例えば、検査された領域全体に係る、例えば１つの超音波画像に係る複数の減衰パラメータの局所的減衰推定方法と呼ばれてもよい。上記パラメータは、超音波画像における画素又はボクセルに割り当てられてもよい。さらに、それらは、検査された領域の追加の減衰マップの超音波画像におけるセグメント又はスーパー画素又はクラスタを形成してもよい。

この方法のおかげで、関心対象領域に係る局所的超音波減衰パラメータを推定することが可能になり、領域にわたる減衰特性のマップ及び／又は１つ又は複数の画像を生成できるようになる。

数式Ｆ（Ｅｙ（ｚ），Ｂ）を解くために使用されうる従来のインバージョン法と比較して、提案した方法は、より速く収束するので、より少ない計算負荷を意味する。同時に、それらは、例えば、取得した局所推定値のより高い空間分解能について、より正確な結果をもたらす。

自己相関関数は、局所的減衰関数を取得するために、深度方向（ｚ）において離散化されてもよい。

駆動入力は、基準媒体の予め決められた減衰パラメータを含んでもよい。

モデリングされた状態空間モデルは、カルマンフィルタ又は予め定義されたベイジアン法を用いて、局所的超音波減衰パラメータを決定してもよい。

カルマンフィルタを使用することは、最小平均二乗誤差基準を最適化するという利点を有しうる。

媒体の反射率パターンは、予め定義された収束しきい値が達成されるまで、モデリングされた状態空間モデルに対してＥＭアルゴリズムを繰り返すことで補償されてもよい。例えば、ＥＭアルゴリズムは、反射率パラメータを出力する第１のカルマンフィルタを使用してもよく、駆動入力として超音波減衰パラメータを使用してもよい。第２のカルマンフィルタは、超音波減衰パラメータを出力してもよく、駆動入力として反射率パラメータを使用してもよい。両方のフィルタは、例えば、予め定義された収束しきい値が達成されるまで、駆動入力として他方の出力をそれぞれ用いるように接続されてもよい。

本開示はさらに、媒体の領域内における局所的超音波減衰の画像生成を行うための方法に関する。本方法は、前述した局所的減衰推定方法を含んでもよい。減衰画像は、例えば各軸線について、局所的減衰パラメータに基づいて構築されてもよい。

媒体の領域内における局所的超音波減衰の画像生成を行うための方法は、少なくとも１つの予め定義されたカテゴリにそれぞれ関連付けられた１つ又は複数の関心対象領域に画像をセグメント化するようにトレーニング（学習）された機械学習モデルに減衰画像を供給することをさらに含んでもよい。

上記機械学習モデルは、例えば人工ニューラルネットワーク、例えば畳み込みニューラルネットワークを含んでもよい（モデルは例えばコンピュータで実装されてもよい）。モデルは、トレーニングデータとして例えばそれぞれ注釈された画像を用いる教師ありトレーニング方法によって、又は、教師なしトレーニング方法によってトレーニングされてもよい。

最後に、本開示は、（例えば上記媒体と所定関係にある）少なくとも１つの超音波トランスデューサ（２）に関連付けられているか接続されている処理装置を備える、媒体中の領域における超音波減衰パラメータを推定するためのシステムに関する。本システム又は処理装置は、
（ａ）上記媒体においてトランスデューサにより少なくとも１つのパルスを送信し、
（ｂ）上記パルスに応答してトランスデューサによりデータを取得し、
（ｃ）上記処理装置によりデータを処理して上記領域の後方散乱取得データを提供し、
（ｄ）時空間領域における深度の関数である、上記後方散乱された取得データの自己相関関数をゼロ（０）遅延で決定し、
（ｅ）上記自己相関関数に基づいて超音波減衰パラメータを推定するように構成される。

本システムは、オプションで、少なくとも１つのトランスデューサを備えてもよい。

少なくとも１つのトランスデューサは、パルスを送信して組織の応答を受信するように構成された単一のトランスデューサであってもよい。例えば、合焦されたトランスデューサは、例えば、凹面形式を有してもよく、又は、各レンズを有してもよい。また、単一のトランスデューサを掃引することも可能である。

複数のトランスデューサ及び／又はトランスデューサアレイ２を用いることも可能である。例えば、典型的には、軸Ｘ（水平又はアレイ方向Ｘ）に沿って並置された数十個のトランスデューサ（例えば１００～３００）を含む、リニアアレイが提供されてもよい。本開示の実施例では、３Ｄプローブが使用されてもよい。

パルスを送信して応答を受信するために１つ又は複数の同じトランスデューサが使用されてもよく、又は、送信及び受信のために異なるトランスデューサが使用されてもよい。

本開示はさらに、コンピュータによって実行されるとき、上述した方法のうちの少なくとも１つに係るステップを実行するための命令を含むコンピュータプログラムに関してもよい。

最後に、本開示はさらに、コンピュータによって読み取り可能な記録媒体であって、コンピュータによって実行されるとき、上述した方法のうちの少なくとも１つに係るステップを実行するための命令を含むコンピュータプログラムをその上に記録した記録媒体に関してもよい。

矛盾しない限り、上述した構成要素及び本明細書内で説明したものが組み合わされもよいことが意図される。

先述の概要説明及び後述の詳細説明の両方は例示及び説明にすぎず、実例の提示を目的とし、請求項に記載された開示の限定ではないことは理解されるべきである。

添付図面は、本明細書に組み込まれ、本明細書の一部を構成し、その説明とともに本開示の実施形態を示し、その原理をサポート及び例証として提示される。

本開示の実施形態に係る超音波装置を示す概略図である。 図１の装置の一部を示すブロック図である。 本開示に係る超音波減衰パラメータを推定するための方法であって、図１の装置で実装された方法を示す図である。 本開示に係る方法の第１の例示的な実施形態（大域的減衰推定）を示す図である。 図３又は図４の方法によって取得された超音波画像の例を示す。 本開示に係る方法の第２の例示的な実施形態（局所的減衰推定）を示す図である。

ここで、本開示の例示的な実施形態が詳細に参照され、その実施例は添付図面に図示されている。図面の全体にわたって可能であればどこでも、同じ部分又は同様の部分を示すために同じ参照番号が使用される。

図１に示す装置は、媒体の領域１、例えば、生体組織及び特に患者のヒト組織の画像生成を行うように適応化される。本装置は、例えば、以下の構成要素を含んでもよい。

・少なくとも１つのトランスデューサは、例えば、パルスを送信して組織の応答を受信するように構成された単一のトランスデューサ。また、複数のトランスデューサ及び／又はトランスデューサアレイ２を用いることも可能である。例えば、典型的には、通常のプローブで既知のように、軸Ｘ（水平又はアレイ方向Ｘ）に沿って並置された数十個のトランスデューサ（例えば１００～３００）を含む、リニアアレイが提供されてもよい。この場合、アレイ２は、領域１の二次元（２Ｄ）イメージングを実行するように適応化されるが、アレイ２は、領域１の３Ｄイメージングを実行するように適応化された二次元アレイであってもよい。トランスデューサアレイ２は、曲線に沿って整列された複数のトランスデューサを含む凸面のアレイであってもよい。パルスを送信して応答を受信するために１つ又は複数の同じトランスデューサが使用されてもよく、又は、送信及び受信のために異なるトランスデューサが使用されてもよい。
・トランスデューサアレイを制御してそこから信号を取得する電子的ベイ３。
・電子的ベイ３を制御し、例えば、電子的ベイから取得された画像を見るためのマイクロコンピュータ４（変形例では、単一の電子装置が、電子的ベイ３及びマイクロコンピュータ４のすべての機能を満たしてもよい）。

図１における軸Ｚは、軸Ｘに垂直な軸であり、それは、通常、アレイのトランスデューサによって生成された超音波ビームの方向であり、例えば、検査された媒体の深度方向にある。この方向は、本開示では、垂直方向又は軸方向と呼ばれる。

図２に示すように、電子的ベイ３は、例えば、以下の構成要素を含んでもよい。

・トランスデューサアレイ２のＬ個のトランスデューサ（Ｔｌ－ＴＬ）に個々に接続されたＬ個のアナログ／ディジタル変換器５（Ａ／Ｄｉ－Ａ／ＤＬ）；
・ｎ個のアナログ／ディジタル変換器５にそれぞれ接続されたＬ個のバッファメモリー６（Ｂｉ－Ｂｎ）；
・バッファメモリー６及びマイクロコンピュータ４と通信する中央処理装置８（ＣＰＵ）；
・中央処理装置８に接続されたメモリ９（ＭＥＭ）；
・中央処理装置８に接続されたディジタル信号プロセッサ１０（ＤＳＰ）。

本願で開示した装置は超音波イメージングのための装置であり、トランスデューサは超音波トランスデューサであり、実現された方法は、領域１に係る超音波減衰パラメータを推定し、また、オプションで、領域１の超音波画像を生成してもよい。

しかしながら、本装置は、超音波とは異なる波（超音波波長とは異なる波長を有する波）を用いる任意のイメージング装置であってもよく、この場合、トランスデューサ及び電子的ベイ構成要素は上記波に適応化される。

図３は、本開示に係る超音波減衰パラメータを推定するための方法であって、図１の装置で実装された方法を示す図である。

方法ステップは、主として中央処理装置８によって、最終的には、ディジタル信号プロセッサ１０の支援を受けて、又は他の任意の手段によって制御される。本方法は、
（ａ）上記媒体において少なくとも１つのパルスがトランスデューサによって送信される送信ステップ（１０１）と、
（ｂ）上記パルスに応答してデータがトランスデューサによって取得される受信ステップ（１０２）と、
（ｃ）上記処理装置によって上記データが処理されて上記領域の後方散乱取得データを提供する処理ステップ（１０３）と、
（ｄ）空間－時間領域における深度の関数である、上記後方散乱取得データの自己相関関数が決定される関数決定ステップ（１０４）と、
（ｅ）上記自己相関関数に基づいて超音波減衰パラメータが推定される減衰推定ステップ（１０５）とを含む。

オプションで、超音波画像は、ステップ１０３の後方散乱された取得データに基づいて生成されてもよい。

例えば、推定された超音波減衰パラメータは、以下の目的で使用されてもよい。

・それが大域的推定量を提供する場合、肝臓脂肪症を評価すること、
・それが大域的推定量を提供する場合、デュシェンヌ型筋ジストロフィーキャリアを検出すること、及び／又は
・それが減衰のマップを提供する場合、１つ又は複数の乳房腫瘤を特徴づけること。

それらの実施例は、本開示で説明する方法の実証的なアプリケーションとして提供されるが、他のアプリケーション、例えば動物の身体及び／又は材料分析に関連したものもまた同様に関心対象であってもよい。

ステップ１０１～１０５はオプションでループを形成してもよく、例えば、ステップ１０５からステップ１０１に戻ることでループ１０７を形成してもよい。このことで、複数の減衰パラメータを推定できるようになる可能性があり、オプションで、各サイクルのステップ１０３において、さらに、超音波画像が生成されてもよい。追加的又は代替的に、ステップ１０５からステップ１０４に戻ることで、ステップ１０４及び１０５をループするループ１０８を追加することが可能である。このことは、図６のコンテキストでさらに詳細に説明するように、複数の局所的推定量（オプションで、関心対象領域の異なるエリアにおける推定量）を推定できるようになる可能性がある。

いくつかの実施形態では、上記方法は、図４及び図５のコンテキストでさらに詳細に説明するように、大域的超音波減衰パラメータを推定するために使用されてもよく、又は、図６のコンテキストで説明するように、局所的超音波減衰パラメータを推定するために使用されてもよい。しかしながら、本開示の詳細な実施形態について詳しく説明する前に、以下、いくつかの基本について提示する。

以下の説明では、次の用語を用いる。

ＡＣＦ：自己相関関数（Auto Correlation Function）
β：超音波減衰（ｄＢ／ｃｍ／ＭＨｚ又はＮｐ／ｃｍ／ＭＨｚ）
ｃ：媒体における音速
ｆ_０：パルスの中心周波数
λ：パルスの波長（λ＝ｃ／ｆ_０）
ＭＭＳＥ：最小平均二乗誤差（Minimum Mean Square Error）
ＰＳＤ：パワースペクトル密度（Power Spectral Density）
ＲＯＩ：関心対象領域（Region Of Interest）
ＲＦデータ：プローブ上で取得されたディジタル化信号を格納する二次元アレイ（トランスデューサ対時間）
σ_ｆ：ガウス形状スペクトラムの標準偏差。

以下では、媒体が均質であり、複数の弱い散乱器を含むと仮定する。送信パルスの解析的信号は、ガウシアン包絡線及び中心周波数ｆ_０を有すると仮定する。

以下では、本開示で使用されうる一般的な信号モデリングについて説明する。それによって、信号モデリングは、モデルに対して超音波エコー（すなわちＰＳＤ）、超音波減衰、回析、反射率を次第に追加することで、次第に向上される。

まず、画素において測定されたパワー信号密度（Power Signal Density：ＰＳＤ）がどのようにモデリングされうるかについて説明する。この説明セクションでは、簡単化のため、超音波減衰は考慮しない。合焦された開口について、時刻ｔにおいて等時性の体積によって後方散乱された信号は、次式のようにモデリングされてもよい。

ここで、ｎ（ｔ）は、微視的な自己相関関数、すなわち次式を満たす、ゼロ平均のガウス確率過程である。

従って、後方散乱された信号の自己相関関数の表現は、次式のように定式化されてもよい。

ウィーナー・ヒンチン定理を用いると、後方散乱された信号のＰＳＤであるＳ（ｆ）は、次式のように表すことができる。

ビームフォーミング処理は、関心対象の画素からトランスデューサへの伝搬時間に対応する遅延を有してトランスデューサにおいて取得された複数の信号を蓄積することにある。従って、画素において測定されたＰＳＤは、数式（５）によって与えられる。

次に、超音波減衰がどのようにモデリングされうるかについて説明する。特に、超音波減衰は、次式の伝達関数を有するフィルタとしてモデリングされてもよい。

係数２は往復の伝搬を表わす。軸方向の送信及び受信ビームを仮定して、（５）及び（６）の組み合わせは、深度ｚにおいてビームフォーミングされた画素のＰＳＤの表現を導く。

ここで、次式を用いる。

数式（７）で与えられるＰＳＤにおいてスペクトル領域から時空間領域へのウィーナー・ヒンチン定理を呼び出すことで、次式が得られる。

遅延τ＝０において得られたＡＣＦは、画素値の平均パワーを与える。数式（８）を用いて次式が得られる。

（９）から（１０）への変形は、（９）の右辺がガウス関数の－∞から＋∞までの積分であり、従って、その積分値が１になるという事実によって正当化される。

次に、画素信全体号がどのようにモデリングされうるかについて説明する。トランスデューサの感度、ビームフォーミングアンテナの利得、回析を考慮し、数式（１０）と組み合わせることで、次式の表現を導く。

数式（１１）は、その文献で説明された方法の核心を表す。そのフレームワークにおいて、ｃｏｎｓｔａｎｔ（ｚ）が推定されなければならない。

数式（１１）に基づく超音波減衰推定方法が以下の利点を提示するということが本開示の発明者によって発見された。

・以下の理由により、回析効果にはそれほど敏感でない：
－合成ビームフォーミングデータに対して動作すること；従って、任意の空間的場所において、送信及び受信の両方の焦点が保証されること。
－時空間領域において動作し、回析効果はパルスの周波数範囲にわたって平均化されること。
－送信パルスが広帯域であること。
・周波数に基づく方法とは異なり、時間（又は空間的）周波数の妥協を必要としない。
・ビームフォーミングされたデータにおける局所的エネルギー推定は簡単かつ頑健であり、その結果、減衰推定量はより頑健になる。

図４を参照して、大域的超音波減衰パラメータを推定するための方法について説明する。従って、図４は、本開示に係る方法の第１の例示的な実施形態（大域的減衰推定）を示す図である。特に、図４の方法は、図３の方法の特定の実施形態であってもよい。

図４の方法は、受信ステップｂ、ここでは参照数字’１０２及び１０２’で開始し、このとき、データは、パルスに応答してトランスデューサによって取得される。このデータは、ここでは、検査された媒体に係る「ＲＦデータ集合」と呼ばれ、また、基準媒体、例えばファントムに係る「ＲＦトレーニングデータ集合」と呼ばれる。ステップ１０２～１０４は、検査された媒体のデータの処理を示し、一方、同様のステップ１０２’～１０４’は、基準媒体のデータの処理を示す。

次のステップ１０３及び１０３’「合成ビームフォーミング」はそれぞれ、本開示及び図３のステップ（ｃ）、すなわち、処理装置によってデータが処理されて領域の後方散乱された取得データを提供する処理ステップに対応する。

次のステップ１０４、１０４’「平均パワー推定」はそれぞれ、本開示及び図３のステップ（ｄ）、すなわち、時空間領域における深度の関数である、後方散乱取得データ（特に、ステップ１０３のビームフォーミングされた取得データ）の自己相関関数が決定される関数決定ステップに対応する。

次のステップ１０４ａにおいて、２つの自己相関関数が互いに比較され、後方散乱された取得データにおける回析パターン及び／又は未知のトランスデューサ感度を補償する。特に、両方の関数の間のデルタ関数が決定されてもよい。上記比較ステップの例示的な実施形態を以下で説明する。

ステップ１０５において、２つの自己相関関数の比較、すなわち、特にそれらのデルタ関数を解くことで、超音波減衰パラメータが推定される。より詳しくは、以下の通りである。

以下の説明では、次の基本的仮定を用いる。

・媒体における減衰が一定であること。
・数式（１１）のｃｏｎｓｔａｎｔ（ｚ）はｚに依存しない、すなわち、それは関心対象領域ＲＯＩにわたって一定であること。

図３及び図４の方法の目標は、特定のＲＯＩにおいてβを推定することである。ＲＯＩは、図５に示すような肝臓の均質部分であってもよい。

ｚ_０及びｚ_１は（それぞれ）、ＲＯＩの頂部及び底部の深度を（それぞれ）示すものとする。深度は、プローブの表面と関心対象点（空間的場所）との間の距離を意味する。

深度ｚ_０及びｚ_１における平均パワーの対数の差を評価することで、次式が得られる。

数式（１２）（これはステップ１０４の結果として理解されてもよい）は、２つの未知数、すなわち、ｃｏｎｓｔａｎｔ（ｚ_１，ｚ_０）及びβを含む。

βを推定するために、ｃｏｎｓｔａｎｔ（ｚ_１，ｚ_０）を推定する必要がある。ｃｏｎｓｔａｎｔ（ｚ_１，ｚ_０）はいくつかのファクターに依存するが、媒体には依存せず、従って、それは、既知のβに関連付けられたトレーニングデータ集合（すなわちステップ１０２’の基準データ）において学習されてもよい。このことは、減衰について較正されたファントム（すなわち基準媒体）において超音波データを取得することで達成可能である。δＥ_ｒｅｆ（ｚ_０，ｚ_１）は、較正されたファントムにおける深度ｚ_０及びｚ_１で測定された平均パワーの対数の差を示し、β_ｒｅｆは、その較正された（既知の）減衰を示し、ΔＥ（ｚ_０，ｚ_１）は、δＥ（ｚ_０，ｚ_１）及びδＥ_ｒｅｆ（ｚ_０，ｚ_１）の差を示すものとする。次いで、それは次式で決定可能である。

数式（１３）は、比較ステップ１０４ａの例示的な結果であってもよい。数式（１３）（すなわちステップ１０５）の検査は、ΔＥ（ｚ_０，ｚ_１）＝０を解くことがβの推定値を与えることを明らかにする。

数式（１３）は次式で近似されうることに注意する。

このことは、数式（１３）の二次の項が比較的小さく、したがって無視できるという状況に起因する。

この方法の重要なステップは、深度ｚ_０及びｚ_１における平均パワーを推定することにあるということに注意する。その推定処理の第１ステージは、ステップ１０２の取得された超音波ＲＦデータに対して、ステップ１０３で合成ビームフォーミングを行うことであってもよく、その目標は、回析効果を最小化することにある。第２ステージは、特定の深度におけるこれらのビームフォーミングされたデータの平均パワーを推定することであってもよい。しかしながら、ビームフォーミングされたデータは、超音波イメージングにおける公知の問題、すなわち、それらがスペックルノイズによって損なわれるという問題の影響を受ける。スペックルノイズは、空間分解能を犠牲にして空間的平均化を行うことにより平滑化されうる。代替として、ビームフォーミングされた取得データを得るためにビームフォーミングを使用することが可能である。この場合、平均パワー推定値の変動を減少させるために、合成ビームフォーミングステージにおいて開口における直交アポディゼーションを使用することが可能になる。

本開示の方法、特に、上述した例示的な方法は、次の利点を有する。

・方程式ソルバ（方程式解計算器）は、（例えば、比較的小さな二次の項に起因して）数式（１３）及び（１４）の両方に基づくものであってもよい。
・スペクトル解析は実行されず、従って、時間／空間対周波数分解能のトレードオフを考慮する必要はない。
・回析を最小化するために、広帯域のパルスが使用されてもよい。
・送信及び受信の両方の合焦を保証するために、従って、回析を低減するために、合成ビームフォーミングが実行されてもよい。
・スペックルを平滑化して平均パワー変動推定値を減少させるために、従って、超音波減衰推定値の変動を低減するために、直交アポディゼーションが使用されてもよい。
・回析又はプローブ感度のような、媒体に依存しないパラメータが、較正されたファントムにおいて学習される。

図６は、本開示に係る方法の第２の例示的な実施形態（局所的減衰推定）を示す図である。

図６の方法は図４のもの、すなわち、ステップ１０２～１０４及びステップ１０２’～１０４’に対応する。しかしながら、（図４のステップ１０４ａ及び１０５の代わりに）後のステップ１０６において、局所的減衰パラメータを推定するために、両方の自己相関関数がカルマンフィルタを用いて比較される。

図６の方法の目標は、１つ又は複数の局所的減衰パラメータを、例えば局所的推定値のマップの形式で取得することにある。このことは、基本的な仮定が図４の方法のものと比較してわずかに異なることを意味する。特に、それは、ここでは、超音波減衰パラメータが空間的場所に応じて変化するということを考慮しなければならない。

ここで、深度方向ｚにおいて軸線に沿った局所的超音波減衰を推定し、ラインに沿って走査すること、すなわち、Ｂモードイメージングで通常行われるものと同様の方法で走査することが提案される。

従って、この処理の主な目標は、特定の走査線における超音波減衰パラメータを、深度の関数、すなわちβ（ｚ）として推定することにある。このため、数式（１１）の逆関数を求めることを提案する。この数式から直接的に逆関数を求めることは容易ではない。したがって、ノイズを含む状態空間系を観測（又は測定）することによって、平均パワーの対数を深度の関数としてモデリングすることが提案される。

深度の関数としての平均パワーの対数は離散化されてもよく、次式の線形確率差分方程式（１５）によって支配される、離散時間制御された処理の測定としてモデリングされてもよい。

ここで、
・Ｘｋは（隠れた）状態ベクトルである。
・Ａは遷移行列である。
・ｗ_ｋ及びｖ_ｋは処理及び測定のノイズである。
・Ｈは測定行列である。
・ＬＥ_ｋは、テスト対象媒体及び深度ｋΔｚにおける既知の減衰を有する基準媒体のインデックスｋにおける平均パワーの対数と、深度ｋΔｚにおける基準媒体の平均パワーの対数との差を含むベクトルである。
・Δｚは深度サンプリングステップである。
・Ｂは制御行列である。
・ｕ_ｋは駆動入力である。

Ｅ［ｋ］は、深度ｋΔｚにおいて平均パワー推定値の対数を示すものとする。数式（１１）のテーラー展開を行うことにより、次式が得られる。

ここで、次式のように近似されてもよい。

ΔＥ［ｋ］は、テスト対象媒体及び既知の減衰β_ｒｅｆを有する基準媒体、例えば較正されたファントムのインデックスｋにおける平均パワーの対数の差を示すものとする。この場合、次式が得られる。

ｖ_ｋにおいて（ｃｏｎｓｔａｎｔ（ｋ）－ｃｏｎｓｔａｎｔ（ｋ－１））が含まれる場合、次式が得られる。

ここで、
・Ｅ_ｒｅｆは、基準媒体についてインデックスｋにおいて推定された平均パワーの対数である。
・ｎ_１［ｋ］、ｎ_２［ｋ］、ｎ_３［ｋ］は、状態ノイズ、すなわち、ゼロ平均の白色ガウス雑音である。

次式を満たす場合、数式（１７）は、数式（１５）の形式にすることができる。

カルマンフィルタが、ＭＭＳＥの意味で、数式（１５）における状態ベクトルｘ_ｋの最適な推定量を提供することがわかった。その結果、β［ｋ］の最適な推定量（ＭＭＳＥの意味で）を提供する。

請求項を含む本開示の全体にわたって、用語「１つの～を備える（含む）」は、他の記載がなければ、「少なくとも１つの～を備える（含む）」と同義のものとして理解されるべきである。さらに、請求項を含む本開示で示した任意の範囲は、他の記載がなければ、その端の値を含むものとして理解されるべきである。説明した構成要素に係る特定の値は、当業者に既知である、許容された製造時又は産業上の公差内にあるものと理解されるべきであり、用語「実質的に」及び／又は「約」及び／又は「概して」のいかなる使用も、そのような許容された公差内にあることを意味するものと理解されるべきである。

本開示は特定の実施形態を参照してここに説明したが、これらの実施形態は、本開示の原理及び応用の単なる例示であることが理解されるべきである。

明細書及び実施例が例示的なものとしてのみ考慮され、本開示の真の範囲は添付の特許請求の範囲によって示されることが意図される。

要約すると、上述した本開示に係る方法は、より正確な減衰推定を可能にし、より少ない計算上のコストを意味し、このことは、特に、リアルタイム計算モードを改善する。さらに、増大した精度に起因して、変動の減少を達成することができ、したがって、増大した再現性を達成することができる。

Claims (17)

1. 媒体中の領域における超音波減衰パラメータを推定するための方法であって、上記方法は、少なくとも１つの超音波トランスデューサ（２）に関連付けられた処理装置（８）によって実施され、上記方法は、
   （ａ）上記媒体において少なくとも１つのパルスがトランスデューサによって送信される送信ステップ（１０１）と、
   （ｂ）上記パルスに応答してデータがトランスデューサによって取得される受信ステップ（１０２）と、
   （ｃ）上記処理装置によって上記データが処理されて上記領域の後方散乱された取得データを提供する処理ステップ（１０３）と、
   （ｄ）時空間領域における深度の関数である、上記後方散乱された取得データの自己相関関数がゼロ遅延で決定される関数決定ステップ（１０４）と、
   （ｅ）上記自己相関関数に基づいて超音波減衰パラメータが推定される減衰推定ステップ（１０５）とを含む、
   方法。
2. 上記処理ステップ（ｃ）は、上記データがビームフォーミング処理によって処理されて上記領域のビームフォーミングされた取得データを提供するビームフォーミングステップ（１０３）を含む、
   請求項１記載の方法。
3. 上記減衰推定ステップ（ｅ）において、超音波減衰パラメータは、深度方向の所定区間（Δｚ）にわたって推定される、
   請求項１又は２記載の方法。
4. 上記減衰推定ステップ（ｅ）において、上記超音波減衰パラメータは、上記自己相関関数を（区間（Δｚ）について）解くことで決定される、
   請求項１～３のうちの１つに記載の方法。
5. 上記後方散乱された取得データの自己相関関数は、深度及び時間の関数である、
   請求項１～４のうちの１つに記載の方法。
6. 上記関数決定ステップ（ｄ）において、上記後方散乱された取得データにおける残留回折パターンは、予め決められた回析基準データに基づいて補償される、
   請求項１～５のうちの１つに記載の方法。
7. 上記後方散乱された取得データにおける残留回析パターンは、予め決められた基準超音波減衰パラメータを有する基準媒体に対してステップ（ａ’）～（ｅ’）を含む方法を適用し、上記自己相関関数を比較することで補償される、
   請求項１～６のうちの１つに記載の方法。
8. 上記自己相関関数を比較することは、上記媒体及び上記基準媒体における２つの異なる深度である自己相関関数及びｚ_０及びｚ_１の間のデルタ関数ΔＥ（ｚ_０，ｚ_１）を計算することと、上記デルタ関数ΔＥ（ｚ_０，ｚ_１）＝０（０１０２）を解いて上記媒体の超音波減衰パラメータを決定することとを含む、
   請求項７記載の方法。
9. 上記予め決められた回析基準データは、基準媒体における少なくとも２つの異なる深度（ｚ_０，ｚ_１）に係る予め決められたパラメータを含む、
   請求項６～８のうちの１つに記載の方法。
10. 上記パルスは、広帯域パルス及び／又は複数のパルスである、
    請求項１～９のうちの１つに記載の方法。
11. 媒体中の領域における局所的超音波減衰パラメータを推定するための方法であって、上記方法は、請求項１～１０のうちの１つに記載の方法を含み、
    上記方法は、上記媒体内の深度方向（ｚ）における複数の軸線のそれぞれについて、
    上記自己相関関数を、駆動入力を有するノイズありの状態空間モデルとしてモデリングすることと、
    上記自己相関関数を、チホノフ正則化法、リッジ回帰法、及び／又はベイジアン回帰法を含む正則化法によって正則化することと、
    サビツキー・ゴーレイ（Ｓａｖｉｔｚｋｙ－Ｇｏｌａｙ）フィルタを用いて上記自己相関関数をフィルタリングすることと、モンテカルロ法と、特異スペクトル分析と、
    圧縮された検出及び／又は疎サンプリング法と、
    のうちの少なくとも１つを実行することにより、Ｅｙを上記自己相関関数としかつｚを上記深度とするとき、数式Ｆ（Ｅｙ（ｚ），β）＝０における局所的超音波減衰βを推定することを含む、
    方法。
12. 上記自己相関関数が、ノイズありの状態空間モデルとしてモデリングされる場合、上記駆動入力は、基準媒体の予め決められた減衰パラメータを含む、
    請求項１１記載の方法。
13. 上記モデリングされた状態空間モデルは、カルマンフィルタ又は予め定義されたベイジアン法を用いて、局所的超音波減衰パラメータを決定する、
    請求項１１又は１２記載の方法。
14. 上記媒体の反射率パターンは、予め定義された収束しきい値が達成されるまで、上記モデリングされた状態空間モデルに対してＥＭアルゴリズムを繰り返すことで補償され、
    上記ＥＭアルゴリズムは、反射率パラメータを出力しかつ駆動入力として超音波減衰パラメータを用いる第１のカルマンフィルタと、超音波減衰パラメータを出力しかつ駆動入力として反射率パラメータを用いる第２のカルマンフィルタとを使用し、
    両方のフィルタは、駆動入力として他方の出力をそれぞれ用いるように接続される、
    請求項１１～１３のうちの１つに記載の方法。
15. 媒体の領域内における局所的超音波減衰の画像生成を行うための方法であって、請求項１１～１４のうちの１つに記載の方法を含み、
    局所的減衰パラメータに基づいて減衰画像が構築される、
    方法。
16. 媒体の領域内における局所的超音波減衰の画像生成を行うための請求項１５記載の方法であって、
    少なくとも１つの予め定義されたカテゴリにそれぞれ関連付けられた１つ又は複数の関心対象領域に画像をセグメント化するようにトレーニングされた機械学習モデルに減衰画像を供給することを含む、
    方法。
17. 媒体中の領域における超音波減衰パラメータを推定するためのシステムであって、上記システムは、少なくとも１つの超音波トランスデューサ（２）に関連付けられた処理装置（８）を備え、
    上記システムは、
    （ａ）上記媒体においてトランスデューサにより少なくとも１つのパルスを送信し、
    （ｂ）上記パルスに応答してトランスデューサによりデータを取得し、
    （ｃ）上記処理装置によりデータを処理して上記領域の後方散乱された取得データを提供し、
    （ｄ）時空間領域における深度の関数である、上記後方散乱された取得データの自己相関関数をゼロ遅延で決定し、
    （ｅ）上記自己相関関数に基づいて超音波減衰パラメータを推定するように構成された、
    システム。

Images