特許 7607881 超音波検査データから波形インバージョンによって高速かつロバストに体組織をイメージングするための方法

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

(12)【公報種別】特許公報(B1)

(11)【特許番号】

(24)【登録日】2024-12-20

(45)【発行日】2025-01-06

(54)【発明の名称】超音波検査データから波形インバージョンによって高速かつロバストに体組織をイメージングするための方法

(51)【国際特許分類】

   A61B 8/13 20060101AFI20241223BHJP

【ＦＩ】

A61B8/13

【請求項の数】 4

(21)【出願番号】P 2023205720

(22)【出願日】2023-12-05

【審査請求日】2023-12-22

【早期審査対象出願】

(73)【特許権者】

【識別番号】592036391

【氏名又は名称】株式会社地球科学総合研究所

(73)【特許権者】

【識別番号】516247100

【氏名又は名称】ユニヴェルシテ・グルノーブル・アルプ

【氏名又は名称原語表記】Ｕｎｉｖｅｒｓｉｔｅ Ｇｒｅｎｏｂｌｅ Ａｌｐｅｓ

(74)【代理人】

【識別番号】110000051

【氏名又は名称】弁理士法人共生国際特許事務所

(72)【発明者】

【氏名】加藤 政史

(72)【発明者】

【氏名】寺西 慶裕

(72)【発明者】

【氏名】疋田 朗

(72)【発明者】

【氏名】淺川 栄一

(72)【発明者】

【氏名】ロマン ブロシエ

(72)【発明者】

【氏名】ルードヴィック メチビエ

(72)【発明者】

【氏名】ジャン ヴィリュー

(72)【発明者】

【氏名】ペンリン ヤン

【審査官】井海田 隆

(56)【参考文献】

【文献】特開２０２０－０１８７８９（ＪＰ，Ａ）

【文献】中国特許出願公開第１１５７３６９８６（ＣＮ，Ａ）

【文献】特表２０２０－５０１７３５（ＪＰ，Ａ）

【文献】米国特許出願公開第２０１５／０２７２５０６（ＵＳ，Ａ１）

【文献】米国特許出願公開第２０１８／０２１７８７３（ＵＳ，Ａ１）

【文献】特開２０２２－０２６１１６（ＪＰ，Ａ）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ａ６１Ｂ ８／００－８／１５

Ｇ０１Ｎ ２９／００－２９／５２

Ｇ０１Ｖ １／００－１／５２

(57)【特許請求の範囲】

【請求項1】

超音波検査データから波形インバージョンによる体組織のイメージング方法であって、
（ａ）超音波検査データを取得する段階、
（ｂ）粘性付き音響波動方程式について、複数の計算機（ＣＰＵ）のｘ方向とｙ方向とにそれぞれ独立に等間隔に分割されたメモリ空間において、ｘ方向に隣り合うＣＰＵ間及びｙ方向に隣り合うＣＰＵ間で、変位速度場及び音圧場のデータ交換を行い、１つの波動場を複数の計算機（ＣＰＵ）で生成する段階、
（ｃ）波動場から受信点で記録される波形について、前記（ａ）段階で取得した観測波形との間の誤差２乗、及び下記式（５）のｆｃを全てのＣＰＵから統合する段階、
（ｄ）前記（ｂ）段階で分割されたメモリ空間において計算した修正量勾配と疑似ヘッセ行列とに、ＣＰＵ数に等しくｘ方向に等間隔にメモリ空間を分割した後、ｙ方向にガウシアンフィルターを適用する段階、
（ｅ）前記（ｄ）段階でｘ方向のみに等間隔に分割したメモリ空間からｙ方向にＣＰＵ数に等しく等間隔なメモリ空間に分割した後に、ｘ方向にガウシアンフィルターを適用する段階、
（ｆ）準ニュートン法によって、全波形インバージョンのモデル更新に必要なステップ長を、前記（ｅ）段階でｙ方向に等間隔に分割されたメモリ空間において計算する段階、及び
（ｇ）前記（ｂ）～（ｆ）段階を繰り返し、音波伝播速度に加えて、粘性Ｑ値のイメージングを得る段階からなり、
前記粘性付き音響波動方程式が下記式（１）～（４）の連立微分方程式によって表現されることを特徴とする高速かつロバストな体組織のイメージング方法。
〔式１〕

〔式２〕

〔式３〕

〔式４〕

ここで、上付きの・は時間方向の微分を表し、ρおよびＶはそれぞれ媒質の密度と音波伝播速度であり、σ、ｖ_ｘ、ｖ_ｙはそれぞれ音圧、ｘ方向の変位速度、ｙ方向の変位速度であり、γｌは記憶変数であり、τｌは粘性Ｑ値を定義するＬ個の周波数（ωｌ）から換算される緩和時間であり、Ｌは緩和メカニズムの個数（ここでは、３個とする）であり、
〔式５〕

ここで、ｆｃ（スカラー量）はコスト関数で、全受信点の全時刻において、観測データ（ｄ＿ｏｂｓ）と初期モデルで計算された合成データ（ｄ＿ｓｙｎ）の残差２乗和、Ｎｓ，Ｎｒ，Ｎｔはそれぞれ発信点数、受信点数、時間サンプル数である。

【請求項2】

前記（ｂ）段階で、各時間ステップにおいてデータの交換を行い、分割したメモリ空間の分割数はＣＰＵ数及びメモリ空間からなる計算機の環境に応じて設定し、前記分割したメモリ空間は全波形インバージョンの順伝播と逆伝播および順伝播の巻き戻し計算の３つの波動場計算全てに適用し、前記分割したメモリ空間のまま、修正量勾配と疑似ヘッセ行列を求めることを特徴とする請求項１に記載の高速かつロバストな体組織のイメージング方法。

【請求項3】

前記（ｆ）段階に計算される修正量勾配の下記（９）に示す２乗和と修正量勾配と疑似ヘッセ行列の下記（８）に示す相乗和は、全てのＣＰＵから統合することを特徴とする請求項２に記載の高速かつロバストな体組織のイメージング方法。
〔８〕

〔９〕

ここでｉはグリッド番号でＮは総グリッド数、Ｈ^－１は、疑似ヘッセ行列の逆行列でＧは修正量勾配である。

【請求項4】

前記超音波検査データは、超音波による乳がん検査データであり、
前記体組織は、皮膚、乳腺、癌化した組織を含むことを特徴とする請求項３に記載の高速かつロバストな体組織のイメージング方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、超音波医療機器で取得されたデータについて、従来の反射面（あるいは散乱点）のイメージ（エコーイメージ）とは異なり、体内を音波伝播速度と減衰Ｑ値でイメージする全波形インバージョンに関するものである。本発明が請求する内容は、超音波検査データに対して、全波形インバージョンを実用化するために必要不可欠である数値演算の高速計算を実現するための手法である。

【背景技術】

【0002】

全波形インバージョンは、地震学の分野で開発された波動論を応用した地下の可視化技術である（非特許文献１、２）。この全波形インバージョンは、原理的に超音波医療検査機器で取得されたデータに適用することが可能で、検査装置と併せた形でいくつかの既往特許が発表されている（特許文献１～６）。

【0003】

超音波で体内をイメージングすることは、Ｘ線を用いるレントゲン検査やＣＴ検査、マンモグラフィー検査に比べて、被曝リスクが無い点で望ましい。超音波医療機器の代表として、エコー検査機がある。これは、リアルタイムで画像を提供することができる。しかしながら、エコー検査は熟練した医師や検査技師が繰り返し検査を実施することで、問題箇所を推定していく作業が必要である。

【0004】

乳がん検査の代表であるマンモグラフィー検査は、デンスブレスト乳房の場合、密集した乳腺と癌化した組織の識別が困難となる。全波形インバージョンでは音波伝播速度と減衰Ｑ値の２つの物理量で体組織をイメージングできることから、デンスブレスト乳房でも、乳腺と癌化した組織の識別を数値的に行うことが可能となる。
ＭＲＩは、Ｘ線ＣＴよりも解像度でやや劣るものの、被曝リスクは無く検査当日中にはイメージング結果が得られる程度の計算コストである。

【0005】

全波形インバージョンは、特許文献５の記載のように、ＭＲＩに匹敵するためにはリアルタイム画像の提供が必要で、リアルタイム画像を提供するためには桁単位でターンアラウンドタイムを低減する必要がある。
また、全波形インバージョンは限られた周波数帯（特にターゲットの大きさに対して数波長以上の伝播距離を要する高周波数帯）のデータを使用する際は、結果の初期モデル依存性が高く、ロバスト性を確保することが困難である。
時間領域の全波形インバージョンは周波数領域のそれに比べて、ロバスト性を大幅に改善するが、計算時間に問題があり十分な検討は為されてこなかった。
以上のような問題点によって、全波形インバージョンは、超音波での体組織のイメージング技術としては現在のところ実用化に至っていない。

【先行技術文献】

【特許文献】

【0006】

【文献】特開２０２０－１８７８９号公報

【文献】特開２０２１－１９８３９号公報

【文献】特開２０１８－５２８８２９号公報

【文献】特開２０２０－５０１６４８号公報

【文献】特表２０２０－５０１７３５号公報

【文献】国際公開第２０１３－１１６８５４号

【非特許文献】

【0007】

【文献】Ｌａｉｌｌｙ， Ｐ．，“ Ｔｈｅ ｓｅｉｓｍｉｃ ｉｎｖｅｒｓｅ ｐｒｏｂｌｅｍ ａｓ ａ ｓｅｑｕｅｎｃｅ ｏｆ ｂｅｆｏｒｅ ｓｔａｃｋ ｍｉｇｒａｔｉｏｎｓ”， Ｃｏｎｖｅｒｅｎｃｅ ｏｎ Ｉｎｖｅｒｓｅ Ｓｃａｔｔｅｒｉｎｇ， Ｔｈｅｏｒｙ ａｎｄ ａｐｐｌｉｃａｔｉｏｎ， Ｓｏｃｉｅｔｙ ｆｏｒ Ｉｎｄｕｓｔｒｉａｌ ａｎｄ Ａｐｐｌｉｅｄ Ｍａｔｈｅｍａｔｉｃｓ， １９８３， ｐ２０６－２２０

【文献】Ｔａｒａｎｔｏｌａ， Ａ．，“ Ｉｎｖｅｒｓｉｏｎ ｏｆ ｓｅｉｓｍｉｃ ｒｅｆｌｅｃｔｉｏｎ ｄａｔａ ｉｎ ｔｈｅ ａｃｏｕｓｔｉｃ ａｐｐｒｏｘｉｍａｔｉｏｎ”， Ｇｅｏｐｈｙｓｉｃｓ， １９８４， ｖｏｌ． ４９， ｎｏ． ８， ｐｐ． １２５９－１２６６

【文献】加藤政史、“列車走行に起因する地盤振動の粘弾性波動論に基づく数値シミュレーションに関する研究“、京都大学大学院工学研究科博士学位論文、２００８

【文献】Ｖｉｒｉｅｕｘ Ｊ．，”Ｐ－ＳＶ ｗａｖｅ ｐｒｏｐａｇａｔｉｏｎ ｉｎ ｈｅｔｅｒｏｇｅｎｅｏｕｓ ｍｅｄｉａ： Ｖｅｌｏｃｉｔｙ－ｓｔｒｅｓｓ ｄｉｎｉｔｅ－ｄｉｆｆｅｒｅｎｃｅ ｍｅｔｈｏｄ” Ｇｅｏｐｈｙｓｉｃｓ， １９８６， ｖｏｌ．５１， Ｐ８８９－９０１

【文献】Ｙａｎｇ， Ｐ．，他３名，“ Ｗａｖｅｆｉｅｌｄ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｉｎ ａｔｔｅｎｕａｔｉｎｇ ｍｅｄｉａ： Ａ ｃｈｅｃｋｐｏｉｎｔｉｎｇ－ａｓｓｉｓｔｅｄ ｒｅｖｅｒｓｅ－ｆｏｒｗａｒｄ ｓｉｍｕｌａｔｉｏｎ ｍｅｔｈｏｄ“，Ｇｅｐｈｙｓｉｃｓ，２０１６， ｖｏｌ． ８１， ｎｏ． ６， ｐｐ． Ｒ３４９－３６２

【文献】Ｙａｎｇ， Ｐ．，他３名，“ Ａ ｒｅｖｉｅｗ ｏｎ ｔｈｅ ｓｙｓｔｅｍａｔｉｃ ｆｏｒｍｕｌａｔｉｏｎ ｏｆ ３－Ｄ ｍｕｌｔｉｐａｒａｍｅｔｅｒ ｆｕｌｌ ｗａｖｅｆｏｒｍ ｉｎｖｅｒｓｉｏｎ ｉｎ ｖｉｓｃｏｅｌａｓｔｉｃ ｍｅｄｉｕｍ”，Ｇｅｏｐｈｙｓｉｃａｌ Ｊｏｕｒｎａｌ Ｉｎｔｅｒｎａｔｉｏｎａｌ， ２０１７， ｖｏｌ． ２０７， ｐ１２９－１４９

【文献】Ｓｈｉｎ Ｃ．， 他２名， “Ｉｍｐｒｏｖｅｄ ａｍｐｌｉｔｕｄｅ ｐｒｅｓｅｒｖａｔｉｏｎ ｆｏｒ ｐｒｅｓｔａｃｋ ｄｅｐｔｈ ｍｉｇｒａｔｉｏｎ ｂｙ ｉｎｖｅｒｓｅ ｓｃａｔｔｅｒｉｎｇ ｔｈｅｏｒｙ．”Ｇｅｏｐｈｙｓｉｃａｌ Ｐｒｏｓｐｅｃｔｉｎｇ， ２００１， ｖｏｌ． ４９．５， ｐ５９２－６０６

【発明の概要】

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

【0008】

本発明の目的とするところは、上記従来技術の課題等に鑑み、超音波検査データから波形インバージョンによって高速かつロバストに体組織をイメージングするための方法を提供することである。

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

【0009】

まず、全波形インバージョンのロバスト性を高めるために、時間領域の計算手法を採用する。時間領域の全波形インバージョンでは、画像ピクセルサイズは１波長に対して３０ピクセル以上を必要とする。これは計算時間の増大につながるが、画像解像度が高まるという利点を生む。
既往特許が引用する音響波動方程式ではなく、粘性付き音響波動方程式を採用する。音波伝播速度に加えて粘性Ｑ値をイメージング結果として提供することができるようになる。

【0010】

次に、計算を高速化するために、多数の計算機において並列計算を行う。個々の計算機は汎用ＣＰＵで構成されたもので良いため、多数の計算機を同時使用したとしても費用は安価である。ただし、全波形インバージョンは膨大なメモリ空間が必要とされるため、メモリの分割を行う必要がある。すなわち、全波形インバージョンの計算を高速化するための手段は、以下のように、メモリ空間の分割方法を全波形インバージョンに組み込むことに要約される。

【0011】

（１）波動方程式（順伝播・逆伝播ともに）を解く際はｘ方向とｙ方向を等間隔にメモリを分割する。分割されたメモリ空間において、初期モデルについての波動方程式を有限差分法で解く。有限差分法のタイムステップ毎に、分割された領域の端部のデータ交換を行う。同時に、観測データとの誤差であるコスト関数を各領域で求めた後、分割されたメモリ領域から１つに集めて足し合わせておく。
（２）（１）で分割された領域のまま修正量勾配（Ｇｒａｄｉｅｎｔ）と疑似ヘッセ行列（Ｐｓｅｕｄｏ Ｈｅｓｓｉａｎ）を計算する。
（３）修正量勾配と疑似ヘッセ行列にｙ方向にガウシアンスムージングを掛ける。この際はｘ方向のみ等間隔に分割する。
（４）（３）につづきｘ方向にガウシアンスムージングを掛ける。この際はｙ方向にのみ等間隔に分割する。
（５）（４）のｙ方向で等間隔に分割されたメモリ領域のまま、疑似ヘッセ行列による修正量勾配を補正する。
（６）（４）のｙ方向で等間隔に分割したまま、修正量勾配のノルムを計算し、分割されたメモリ領域から集めて１つに足し合わせる。
（７）（５）で求めたノルムから修正量の方向ベクトルを計算し、修正量勾配と掛け合わせることで、初期モデルの修正量を得る。得られた修正量を初期モデルとして更新されたモデルを得る。
（８）メモリ分割領域を（１）の形態に戻し、（１）～（７）を繰り返す。
（９）コスト関数が一定の閾値以下に収束するか、繰り返し回数が設定した上限数に達したら、更新作業を終了し、最終的なイメージング結果とする。

【0012】

すなわち、超音波検査データから波形インバージョンによる体組織のイメージング方法は、
（ａ）超音波検査データを取得する段階、
（ｂ）粘性付き音響波動方程式について、複数の計算機（ＣＰＵ）のｘ方向とｙ方向とにそれぞれ独立に等間隔に分割されたメモリ空間において、ｘ方向に隣り合うＣＰＵ間及びｙ方向に隣り合うＣＰＵ間で、変位速度場及び音圧場のデータ交換を行い、１つの波動場を複数の計算機（ＣＰＵ）で生成する段階、
（ｃ）波動場から受信点で記録される波形について、前記（ａ）段階で取得した観測波形との間の誤差２乗、及び下記式（５）のｆｃを全てのＣＰＵから統合する段階、
（ｄ）前記（ｂ）段階で分割されたメモリ空間において計算した修正量勾配と疑似ヘッセ行列とに、ＣＰＵ数に等しくｘ方向に等間隔にメモリ空間を分割した後、ｙ方向にガウシアンフィルターを適用する段階、
（ｅ）前記（ｄ）段階でｘ方向のみに等間隔に分割したメモリ空間からｙ方向にＣＰＵ数に等しく等間隔なメモリ空間に分割した後に、ｘ方向にガウシアンフィルターを適用する段階、
（ｆ）準ニュートン法によって、全波形インバージョンのモデル更新に必要なステップ長を、前記（ｅ）段階でｙ方向に等間隔に分割されたメモリ空間において計算する段階、及び
（ｇ）前記（ｂ）～（ｆ）段階を繰り返し、音波伝播速度に加えて、粘性Ｑ値のイメージングを得る段階からなり、
前記粘性付き音響波動方程式が下記式（１）～（４）の連立微分方程式によって表現されることを特徴とする。
〔式１〕

〔式２〕

〔式３〕

〔式４〕

ここで、上付きの・は時間方向の微分を表し、ρおよびＶはそれぞれ媒質の密度と音波伝播速度であり、σ、ｖ_ｘ、ｖ_ｙはそれぞれ音圧、ｘ方向の変位速度、ｙ方向の変位速度であり、γｌは記憶変数であり、τｌは粘性Ｑ値を定義するＬ個の周波数（ωｌ）から換算される緩和時間であり、Ｌは緩和メカニズムの個数（ここでは、３個とする）であり、
〔式５〕

ここで、ｆｃ（スカラー量）はコスト関数で、全受信点の全時刻において、観測データ（ｄ＿ｏｂｓ）と初期モデルで計算された合成データ（ｄ＿ｓｙｎ）の残差２乗和、Ｎｓ，Ｎｒ，Ｎｔはそれぞれ発信点数、受信点数、時間サンプル数である。

【発明の効果】

【0013】

本発明によれば、超音波検査データから波形インバージョンによって高速かつロバストに体組織をイメージングするための方法を提供できる。

【図面の簡単な説明】

【0014】

【図1】送受信点配置の図である。

【図2】数値モデルの説明図（ａ）と実験に基づいた体組織の物理量（ｂ）（粘性Ｑ値と音波伝播速度Ｖとの関係）である。

【図3】ウェーブレットの図（ａ）波形、（ｂ）スペクトルである。

【図4】スタガードグリッド配置の説明図である。

【図5】波動場のスナップショット（発信時刻から６４μｓ後）の図である。

【図6】 メモリー空間分割の説明図（ａ）分割無し、（ｂ）２分割、（ｃ）４分割、（ｄ）８分割、（ｅ）１６分割である。

【図7】ｘ方向に隣り合うＣＰＵ間のデータ交換ダイアグラムの説明図である。

【図8】ｙ方向に隣り合うＣＰＵ間のデータ交換ダイアグラムの説明図である。

【図9】観測データの図である。

【図10】全波形インバージョンの初期モデルの説明図（ａ）Ｖ音波伝播速度、（ｂ）粘性Ｑ値、（ｃ）ρ媒質の密度である。

【図11】初期モデルで計算された合成データの図である。

【図12】観測データと合成データの残差の図である。

【図13】残差を逆伝播させた波動場のスナップショット（発信時刻から６４μｓ後）の図である。

【図14】修正量勾配の図（ａ）ρＶ^２、（ｂ）Ｑ^－１である。

【図15】計算するメモリー空間の分割方法の説明図（ａ）送信点１～Ｎ／４の４分割、（ｂ）計算するメモリー空間のｘ方向分割、（ｃ）計算するメモリー空間のｙ方向分割である。

【図16】結果イメージの図（ａ）Ｖ音波伝播速度、（ｂ）粘性Ｑ値である。

【図17】コスト関数のグラフである。

【図18】モデル更新の説明図である。

【図19】本発明のフローチャートである。

【図20】ベンチマークのグラフである。

【図21】周波数領域波動場のスナップショットの図（ａ）順伝播、（ｂ）逆伝播である。

【図22】周波数領域の修正量勾配の説明図（ａ）Ｖ音波伝播速度、（ｂ）粘性Ｑ値である。

【図23】本発明のフローチャートである。

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

【0015】

形態を説明するため、図１に示す数値モデルおよび圧電素子の配置について、数値シミュレーションを行う。図２に実験に基づいた体組織の物理量を示す。水と腫瘍の音波伝播速度が近いが、粘性Ｑ値は大きな差がある。これは、体組織の水分が多い場所（血管や乳腺）に腫瘍が存在する場合、音波伝播速度のイメージングのみでは腫瘍の検出が困難で、粘性Ｑ値のイメージング結果が腫瘍検出に大きな役割を果たすことが示唆される。

【0016】

図３に示す超音波信号は汎用の圧電素子で生成できる２オクターブ程度かそれ以下の狭帯域（図２の場合、おおよそ２００ｋＨｚ～８００ｋＨｚの範囲であるが、帯域は任意の周波数を選択することができる）の超音波信号である。

【0017】

図１においてジオメトリは２５６個の圧電素子を１つのセグメントグループとして８つのセグメント（時計回りに都合２０４８個の圧電素子）によって、体組織を取り囲んでいることとする。実用上の観点から、セグメント間にはある程度のギャップが存在し、このギャップ区間においてはデータ取得できないこととする（半径１１．７ｃｍ、受信器間隔は０．３３６ｍｍ、セグメント間のギャップは５．７１ｍｍ～６．３８ｍｍ）。

【0018】

全波形インバージョンを理解するために、図１の矢印（ＴＲ１）で示す位置の圧電素子から図５で示す波形の超音波信号が発信された場合を考える。

【0019】

圧電素子から発信された圧力波は、式（１）フックの法則、式（２）記憶変数（γ）および式（３）、（４）運動方程式の４つの連立微分方程式によって表現される粘性付き音響波動方程式を有限差分法で解くことで音圧（σ）と変位速度（ｖ_ｘ、ｖ_ｙ）で構成される全波動場が計算される。
〔式１〕

〔式２〕

〔式３〕

〔式４〕

【0020】

上付きの・は時間方向の微分を表す。ρおよびＶはそれぞれ媒質の密度と音波伝播速度である。σ、ｖ_ｘ、ｖ_ｙはそれぞれ音圧、ｘ方向の変位速度、ｙ方向の変位速度である。Ｌは緩和メカニズムの個数（３個とする）である。厳密な導出方法および粘性Ｑ値を定義するＬ個の周波数（ω_ｌ）から緩和時間（τ_ｌ）に換算する方法は発明者の（非特許文献３）に示される。

【0021】

有限差分法は、上記微分方程式（１）～（４）を直接離散化（∂ｘ→Δｘ；∂ｔ→Δｔ）することで解く手法である。有限差分法の計算精度を高めるために図４に示すスタガードグリッド配置を採用する。スタガードグリッド配置は発明者の指導教官であるＶｉｒｉｅｕｘ， Ｊ．（非特許文献４）によって考案された。離散化精度は空間方向に４次、時間方向に２次を採用する。

【0022】

図５に順伝播計算の発信時刻から６４μｓ後の音圧場のスナップショットを示す。

【0023】

計算を高速化するために、計算するメモリ空間は図６の区切り位置において、複数の計算機に分割されている（図６の場合は分割なしの場合から１６台の計算機で分割する場合を示したが、ＣＰＵとネットワークの能力に応じて分割数は自由に設定できる）。図７および図８にＣＰＵ間でのデータ交換のダイアグラムを示す。

【0024】

図９に２０４８個の圧電素子で受信された音圧の時系列を示す。

【0025】

次に、未知のターゲット（説明のため、図１および図２のモデルをターゲットとする。）で観測された模擬の観測波形（図８）から全波形インバージョンによってイメージングを行う方法ついて説明する。全波形インバージョンには初期モデルが必要であるため、一般的な初期モデルの取得方法である走時断層撮影法（Ｔｒａｖｅｌ ｔｉｍｅ ｔｏｍｏｇｒａｐｈｙ）から音波速度の初期モデルを得る。粘性Ｑ値は走時断層撮影法からは得られないため、図２から腫瘍と皮膚を除いた、水と脂肪から成るものを初期モデルとする。密度は水の値で一定とし、全波形インバージョンでは解かないこととする。

【0026】

図１０に走時断層撮影法で得られる初期モデルを示す。この「初期モデルでの順伝播波動場」を計算すると、観測点において、図９に合成観測データが作成される（図１１）。
図１２に真のモデルで取得された観測データ（図９）と初期モデルから得た合成観測データ（図１１）の残差を示す。初期モデルでは送受信点の間を体組織がさえぎらない場合の音波は、観測データと同じように合成できる。一方で、初期モデルは明瞭な体組織の境界がないため、反射波がほとんど合成されていないことが分かる。

【0027】

残差は全受信点の全時刻において、観測データ（ｄ＿ｏｂｓ）と初期モデルで計算された合成データ（ｄ＿ｓｙｎ）の残差２乗和を下記式５で計算しておく。
〔式５〕

Ｎｓ，Ｎｒ，Ｎｔはそれぞれ発信点数、受信点数、時間サンプル数である。ｆｃはスカラー量で、あとの準ニュートン法で極小化されるべき目的関数であるため、コスト関数ともよぶ。コスト関数は、分割された計算するメモリ空間の値をまとめておく必要があるため、全ての計算機の値を統合し、全ての計算機で同じ値を持つようにしておく。

【0028】

この観測波形残差を逆伝播させる（図１３に残差を逆伝播させた波動場のスナップショット（発信時刻から６４μｓ後）を示す）。「初期モデルでの順伝播波動場」と、「観測波形残差の逆伝播により生成される波動場」の二つについて、イメージングさせる地点において時間方向に相互相関を計算する。この相互相関のゼロラグ値が全波形インバージョンにおけるモデル修正量の勾配である。順伝播の波動場は計算機メモリに記憶するのではなく、発明者の非特許文献５に示される、順伝播波動場の逆戻し計算を行うことで、計算効率良く修正量勾配総量を求める。同時に疑似ヘッセ行列も求める。

【0029】

本発明では音波伝播速度に加えて粘性Ｑ値のイメージングが行われる。そのためには、修正量勾配の総量を音波伝播速度と粘性Ｑ値の修正量勾配に分解する必要がある。この分解方法は、非特許文献１、２では示されていない。そこで発明者は、その方法を非特許文献６に示した。

【0030】

上述の、順伝播、逆伝播、修正勾配量の計算について、全波形インバージョンを（多くの先行技術が使用している）周波数領域ではなく、時間領域で行う技術的特徴を説明するために以下の比較例と共に示す。

【0031】

音波伝播速度と粘性Ｑ値の修正量勾配が求まったら、それぞれの修正量勾配と疑似ヘッセ行列に対してガウシアンフィルター（スムージング）を掛ける。ガウシアンスムージングは、分割されたメモリ空間のままでは、個々の計算機内で計算することはできないため、計算するメモリ空間をｘ方向に分割し、ｙ方向のガウシアンスムージングを行う。次に計算するメモリ空間をｙ方向に分割しｘ方向にガウシアンスムージングを行う（図１５）。ガウシアンスムージングの特性上、この（ｘ、ｙの）順番が逆でも結果は同じである。

【0032】

ガウシアンスムージングを掛けた後、音波伝播速度と粘性Ｑ値に対する２つの疑似ヘッセ行列の逆行列（Ｈ^－１）で修正量勾配（Ｇ）を割ることで、修正量勾配を補正する。この疑似ヘッセ行列による補正操作は、修正量勾配への影響が支配的なヘッセ行列の対角成分のみを使用することで、膨大な計算量が必要なヘッセ行列の逆行列を陽に解くことなく、修正量勾配の解への収束性を高める効果を生むと言えるが（非特許文献７）、数値的に計算を可能にしているだけで、理論的には不完全であり経験的にはうまく行くというものである。また、ジオメトリの不均質性に起因する局所的に誤ったイメージの生成を抑制する効果が期待できる。

【0033】

〔式６〕
Δｍ＝－Ｈ^－１Ｇ （６）
こうして得られた上記式６の修正量（Δｍ）を初期モデルに加えることで修正されたモデルを得ることが出来るが、この修正方法は線形近似であるため、修正されたモデルで数１０回から数１００回の計算を繰り返して、最終的なイメージングを得る。モデル修正の際は、一般的に修正量勾配にステップ長（α）を掛けることで収束を早める操作が為される。繰り返しｋ回目のモデル（ｍ_ｋ）からｋ＋１回目のモデル（ｍ_ｋ＋１）に更新する計算は、下記式７で表される。
〔式７〕
ｍ_ｋ＋１ ＝ ｍ_ｋ ＋ αΔｍ_ｋ （７）
ステップ長はコスト関数（受信点位置のみで評価する）から求めるという考えもある（特許文献５が採用している）が、本来的には系全体で評価した方が良いため修正量勾配と補正後の下記（８）の修正量の積（スカラー）

を下記（９）の修正量勾配の２乗（スカラー）

で割った値を用いる。ｉはグリッド番号でＮは総グリッド数である。コスト関数は収束したかどうかの判定にのみ利用する。ここまでは最急降下法であり、αは目的関数を接線方向に極小解へと導く。ところが、時間領域で算出された修正量勾配には局所的な誤差が少なく、ＬＢＦＧＳ（Ｌｉｍｉｔｅｄ－ｍｅｍｏｒｙ Ｂｒｏｙｄｅｎ Ｆｌｅｔｃｈｅｒ Ｇｏｌｄｆａｒｂ Ｓｈａｎｎｏ）法という準ニュートン法を採用して、αを目的関数の２次関数近似の接線方向に修正することで、収束をさらに早めることが可能である。

【0034】

ただし、繰り返し回数ｋ＝１のときは、コスト関数の履歴が存在しないので最急降下法でモデル更新を行う。さらに繰り返し回数ｋ＝２においては共役勾配法でモデル更新を行う。これにより、ｋ＝３以降はＬＢＦＧＳ法に必要な最低２回のコスト関数および補正量が得られる。

【0035】

繰り返しインバージョンを行う際の計算機ごとの分割メモリ空間は、ガウシアンフィルター後のｙ方向に分割されたままで良い。このインバージョンテクニックを応用したモデル更新の計算においては、修正量勾配の２乗の総和、修正量勾配と補正後の修正量の積の総和を、それらの履歴と比較することで実行される。そのため、ｙ方向に分割されたメモリ空間からＬ２－ノルムを統合して

【0036】

以上のようにして得られたイメージング結果を図１６に示す。

【0037】

図１７にコスト関数を繰り返し回数についてプロットする。繰り返し回数に伴ってコスト関数が減少してゆくことが分かる。

【0038】

図１８に、繰り返し１０回目、２０回目、３０回目の更新モデルを並べる。腫瘍のフォーカスは繰り返し回数１０回以内で捉えられるが、皮膚がフォーカスするまでには、繰り返し回数は２０～３０回が必要である。

【0039】

音波伝播速度は精度よくイメージされる。粘性Ｑ値には一定の誤差がある。この誤差の原因の一部は密度モデルを一定にしていることである。しかし問題の大半は、皮膚の周りを見ても分かるように、イメージが振動を起こしているため、ターゲットの大きさに対して音波の波長が短すぎると言える。しかしながら、粘性Ｑ値イメージの相対的な大小関係は真のモデルを保存しており、画像として診断する上では十分な精度である。

【0040】

図１９にフローチャートをまとめる。本願発明の計算するメモリ空間の分解方法を、全波形インバージョンのフローチャートに沿って図示した。

【0041】

図２０に計算時間およびメモリについてベンチマーク試験を行った結果を示す。ベンチマーク試験ではｘ方向のピクセル数３００５、ｙ方向のピクセル数３００５、時間方向のサンプル数５６２５、発信点数２４点である。計算機はＩｎｔｅｌ社製ＸＥＯＮ（Ｅ５－２６９０ ｖ４クロック周波数２．６ＧＨｚ、１４コア）ＣＰＵが２個、メモリ２５６ＧＢｙｔｅで１機のサーバーを構成し、２４個のサーバーをＬＡＮケーブルで接続したＰＣクラスタを使用した。並列数は２４、４８、９６、１９２、３８４について並列数２４について正規化した値を表示する。それぞれ有限差分法計算時のメモリ空間の分割数が１、２、４、８、１６を意味する。メモリ空間の分割数が２の時は、計算時間は半分以下となり、メモリの削減による計算機内の効率の改善が、通信によるロスを上回っている。メモリ空間の分割数２、４の時は、通信が必要な計算するメモリ空間の辺の長さがすべての計算機で均質であるが、８、１６は中央部を計算する計算機ほど通信する辺の長さが長くなる。この通信量の不均質性によるロスが大きくなることがわかる。しなしながら、計算時間とメモリはともに並列数にほぼ比例して減少させることが出来る。

【0042】

〔比較例〕
〔周波数領域の全波形インバージョン〕
全波形インバージョンを時間領域と周波数領域で行う際の違いを説明する。
図２（ｂ）に示した５つの周波数サンプル（２５１７７０、３７５７４８、５０１６３３、６２５６１０、７５１４９５Ｈｚ）を使用して周波数領域の場合の全波形インバージョンについて、修正量勾配の計算を行った。

【0043】

図２１は、周波数５０１６３３Ｈｚの周波数（コサイン変換）領域の波動場である。周波数領域の波動場は単振動応答であるため、縞模様をしている。注目すべき点として、体組織の外側や、送受信点配置の円の外側に強い応答が現れる。ただし、発信点から開口幅が９０°以内の受信点の記録について逆伝播を行うことにした。これは開口幅１８０°として全受信点の記録を逆伝播させると、修正量勾配の結果が悪化したためである。

【0044】

図２２に修正量勾配を示す。周波数領域の場合、修正量勾配は、順伝播波動場と逆伝播波動場の複素積の実部であるが、体組織の外側は送受信点配置の円の外側も大きな値が残る結果になってしまう。これは、図１４の時間領域で求めた修正量勾配と全く異なる様相を呈する。これは全波形インバージョンで旧来から問題視されているサイクルスキップの問題である。この問題を低減するために、各種のスムージングや正則化計算が考案されているが、発明者が試行した限りでは解決できる方法は見つからなかった。

【0045】

体組織内にも一定の周期で修正量勾配が大きくなり、同心円のような模様を示している。これは使用した５つの周波数サンプルの位相が揃う周期を意味する。残念ながら、上記の範囲で周波数サンプルを変更しても、ターゲットの大きさ（１６ｃｍ）に対して、波長が非常に短い（３ｍｍ程度）ため、位相が揃う（強めあったり弱めあったりする）周期は必ず生じる。これを避けるためには、全ての周波数サンプルを計算すれば良いが、これを試行することは時間領域の全波形インバージョンと同じことで、時間領域の計算と比べて、非常に重くなる。
サイクルスキップの問題は、信号の周波数帯を下げることで解決できるが、それにはターゲットの大きさに対して１～１０波長で透過する音波を生成する必要がある。機械的に低周波数を発生させるためには圧電素子を大型化する必要があるが、これは超音波医療機器のメリットである低費用な検査機器という意義に反する。

【0046】

また、セグメント間のギャップが生成する２次的な波動場が存在する。この人工的に発生するノイズに対して、修正量勾配の反応が大きいことも分かる。
従って、特許文献では周波数領域と時間領域のいずれでも実施可能という記述があるが、周波数領域に関しては、全ての周波数の計算を行うことができれば可能であるという点において不十分である。例えば、特許文献４では、全波形インバージョンを境界領域の特徴付けと称し、Ｘ線ＣＴや走時断層撮影イメージ（初期モデル）の補強を行うことに留めており、初期モデルである程度イメージできていないものが全波形インバージョンでイメージされることはない。これは周波数領域の全波形インバージョンのロバスト性が低いということを暗に示ものである。従って、出願人は、時間領域でのみロバストにイメージング結果を得ることが可能であると主張する（図２３）。さらに、音波伝播速度に加えて、粘性Ｑ値のイメージを得ることで、腫瘍の検出確率を向上されることができると考えられる。

【産業上の利用可能性】

【0047】

本発明によれば、超音波検査データから波形インバージョンによって高速かつロバストに体組織をイメージングするための方法を提供できるため、既往の検査システムのＸ線ＣＴやＭＲＩよりも安全、低コストかつ高解像度に体組織のイメージングができる。特に、従来、時間領域の波形インバージョンの計算時間の問題を大量のＣＰＵに効果的な並列計算をさせることで、既往の検査システムと同程度の計算時間で解析結果が得られるという産業上の利用の可能性は大きいといえる。

【符号の説明】

【0048】

Ｔ／Ｒ１ スナップショット（図５）の発信点、受信点１番目
Ｔ／Ｒ 送受信点
（Ａ） 腫瘍
（Ｂ） 脂肪
（Ｃ） 皮膚
（Ｄ） 水
Ｖ 音波伝播速度
Ｑ 粘性値
ρ 媒質の密度

【要約】      （修正有）

【課題】高速かつロバストな体組織をイメージングするための方法を提供する。
【解決手段】（ａ）超音波検査データ取得、（ｂ）複数のＣＰＵで１つの波動場生成、（ｃ）波動場から受信点で記録される波形について（ａ）段階で取得した観測波形との間の誤差２乗、式（５）のｆｃを全ＣＰＵから統合、（ｄ）修正勾配量と疑似ヘッセ行列をＣＰＵ数の均等な領域にｘ方向に分割しなおす、（ｅ）ｙ方向に分割しなおす、（ｆ）準ニュートン法によって全波形インバージョンのモデル更新に必要なステップ長を計算する、（ｇ）（ｂ）～（ｆ）を繰り返し、音波伝播速度に加えて粘性Ｑ値のイメージングを得る段階からなる。

ここでｆｃはコスト関数で、全受信点の全時刻の観測データ（ｄ＿ｏｂｓ）と初期モデルで計算された合成データ（ｄ＿ｓｙｎ）の残差２乗和、Ｎｓ，Ｎｒ，Ｎｔは発信点数、受信点数、時間サンプル数である。
【選択図】図１９

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】

【図20】

【図21】

【図22】

【図23】