Top > Search of Japanese Patents > PHOTOACOUSTIC INFORMATION PROCESSING METHOD, PHOTOACOUSTIC INFORMATION PROCESSING PROGRAM AND PHOTOACOUSTIC INFORMATION PROCESSING DEVICE > Specification

Specification :(In Japanese)光音響情報処理方法、光音響情報処理プログラムおよび光音響情報処理装置

Country (In Japanese)日本国特許庁(JP)
Gazette (In Japanese)特許公報(B2)
Patent Number P6739099
Publication number P2018-033886A
Date of registration Jul 27, 2020
Date of issue Aug 12, 2020
Date of publication of application Mar 8, 2018
Title of the invention, or title of the device (In Japanese)光音響情報処理方法、光音響情報処理プログラムおよび光音響情報処理装置
IPC (International Patent Classification) A61B   8/13        (2006.01)
FI (File Index) A61B 8/13 ZDM
Number of claims or invention 5
Total pages 18
Application Number P2016-172204
Date of filing Sep 2, 2016
Date of request for substantive examination Jun 3, 2019
Patentee, or owner of utility model right (In Japanese)【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
Inventor, or creator of device (In Japanese)【氏名】椎名 毅
【氏名】近藤 健悟
Representative (In Japanese)【識別番号】100145872、【弁理士】、【氏名又は名称】福岡 昌浩
【識別番号】100091362、【弁理士】、【氏名又は名称】阿仁屋 節雄
【識別番号】100187632、【弁理士】、【氏名又は名称】橘高 英郎
Examiner (In Japanese)【審査官】冨永 昌彦
Document or reference (In Japanese)特開2013-005871(JP,A)
特表2015-507947(JP,A)
米国特許出願公開第2013/0190595(US,A1)
特表2002-536041(JP,A)
米国特許第06216025(US,B1)
国際公開第2014/207440(WO,A1)
Field of search A61B 8/00 - 8/15
JSTPlus/JMEDPlus/JST7580(JDreamIII)
Scope of claims (In Japanese)【請求項1】
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を備え、
前記音圧分布データとしてk空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理方法。
【請求項2】
前記球殻積分にあたり、球殻を複数領域に分割するとともに、各領域で積分変数をそれぞれ選択する
請求項1に記載の光音響情報処理方法。
【請求項3】
前記初期音圧分布の再構成の結果に基づいて光音響画像を生成するとともに、生成した前記光音響画像についての画像データをデータ圧縮して記憶保持する
請求項1または2に記載の光音響情報処理方法。
【請求項4】
コンピュータに、
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を実行させるとともに、
前記音圧分布データとしてk空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理プログラム。
【請求項5】
光音響波の検出信号を取得する信号取得部と、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成する再構成部と、
前記再構成部での処理結果を出力する結果出力部と、を備え、
前記音圧分布データとしてk空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理装置。
Detailed description of the invention (In Japanese)【技術分野】
【0001】
本発明は、光音響情報処理方法、光音響情報処理プログラムおよび光音響情報処理装置に関する。
【背景技術】
【0002】
近年、例えば医療分野において、光音響効果を用いて被検体内部を画像化する光超音波イメージング法(photoacoustic imaging)を利用することが提案されている。光音響効果とは、被検体にレーザ光を照射すると、被検体内部における血液等の光吸収体に吸収され、局所的な温度上昇が熱弾性変形を引き起こし、その熱弾性変形による体積膨張で音響波(超音波)が発生する現象である。光超音波イメージング法では、発生した音響波を受信して、その音源(すなわち、被検体内部における光吸収体)の初期音圧分布を再構成することで、被検体内部を画像化する(例えば、非特許文献1参照)。
【先行技術文献】
【0003】

【非特許文献1】椎名毅、「光音響を用いた新しい超音波診断装置」、メディカル&イメージングNo.3、オプトロニクス社、2015年、p.53-59
【発明の概要】
【発明が解決しようとする課題】
【0004】
光超音波イメージング法の医療適用にあたっては、画像化のために必要となる初期音圧分布の再構成を、画質劣化を抑制しつつ、迅速かつ効率的に行うことが求められる。
【0005】
本発明は、画質劣化を抑制しつつ、迅速かつ効率的に、初期音圧分布の再構成を行うことができる光音響情報処理方法、光音響情報処理プログラムおよび光音響情報処理装置を提供することを目的とする。
【課題を解決するための手段】
【0006】
本発明の一態様によれば、
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を備え、
前記音圧分布データとしてk空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理方法が提供される。
【0007】
本発明の他の一態様によれば、
コンピュータに、
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を実行させるとともに、
前記音圧分布データとしてk空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理プログラムが提供される。
【0008】
本発明のさらに他の一態様によれば、
光音響波の検出信号を取得する信号取得部と、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成する再構成部と、
前記再構成部での処理結果を出力する結果出力部と、を備え、
前記音圧分布データとしてk空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理装置が提供される。
【発明の効果】
【0009】
本発明によれば、画質劣化を抑制しつつ、迅速かつ効率的に、初期音圧分布の再構成を行うことができる。
【図面の簡単な説明】
【0010】
【図1】光音響イメージングの基本的な手順を例示する説明図である。
【図2】本発明の一実施形態に係る光音響情報処理装置を含む光音響イメージング装置の概略構成の例を模式的に示す説明図である。
【図3】順問題による繰り返し推定法を用いて再構成を行う場合の一般的な手法の概要を例示する説明図である。
【図4】本発明の一実施形態に係る光音響情報処理方法の基本的な手順を例示する説明図であり、圧縮センシングを適用して順問題を解いて再構成を行う手法の概要を示す図である。
【図5】圧縮センシングの意義を具体的に説明するための説明図である。
【図6】本発明の一実施形態に係る光音響情報処理方法における球殻積分の具体的な一態様の概要を具体的に説明するための説明図である。
【発明を実施するための形態】
【0011】
<本発明の一実施形態>
以下に、本発明の実施の形態について、図面を参照しながら説明する。

【0012】
(1)光超音波イメージング法の概要
はじめに、光音響効果を用いた光超音波イメージング法(以下、単に「光音響イメージング」ともいう。)の概要を説明する。

【0013】
(光音響イメージングの手順)
図1は、光音響イメージングの基本的な手順を例示する説明図である。
図例のように、光音響イメージングにあたっては、以下に説明する各ステップ(以下ステップを「S」と略す。)を順に経る。

【0014】
(S101:レーザ光照射)
先ず、被検体1に対して光源部11からレーザ光を照射する。
レーザ光の波長は、検出すべき対象の光吸収特性により決まる。生体計測の場合であれば、水の吸収が少ない可視光から近赤外領域が適している。例えば、血管(血液)の描出や酸素飽和度を求める場合は、酸素化ヘモグロビンと脱酸素化ヘモグロビンのスペクトルが交差する800nm付近の複数の波長のレーザ光を用いる。
また、光の吸収により効果的に超音波を発生させるためには、熱の拡散の時間や吸収体内を超音波が伝搬する時間に比べて、レーザ光の照射時間であるパルス幅が十分短い必要があることから、照射すべきレーザ光としてナノ秒幅のレーザ光パルスを用いる。

【0015】
(S102:音響波放射)
レーザ光を照射すると、被検体1の内部では、血液等の吸収体2が光を吸収して、局所的な温度上昇が熱弾性変形を引き起こし、その熱弾性変形による体積膨張で吸収体2が音響波(初期音圧p0(r))を放射する。このとき、光の吸収により発生する音響波(超音波)の大きさは、以下の(1)式で表される。

【0016】
【数1】
JP0006739099B2_000002t.gif

【0017】
ここで、βは熱膨張係数、cは音速、Cpは定圧比熱、μaは光の吸収係数、Ψは光量である。係数G=βc2/Cpは、光エネルギーから音圧への変換係数でグリューナイゼン(Gruneisen)係数と呼ばれる。Gは軟組織間であまり差はないため、初期音圧p0は、組織の吸収係数μaと光量Ψの積に比例することになる。

【0018】
また、吸収体2で発生した位置r、時刻tにおける音響波p(r,t)は、以下の(2)式で表わされる光音響波動方程式に従って、被検体1の内部を伝搬する。

【0019】
【数2】
JP0006739099B2_000003t.gif

【0020】
(S103:音響波計測)
その後は、吸収体2で発生し被検体1の体表まで伝搬した音響波p(r,t)を、体表近傍に配されて複数の圧電変換素子(トランスデューサ)を有するセンサ部12を用いて計測する。これにより、センサ部12では、吸収体2からの音響波p(r,t)の検出信号(受信超音波信号)が電気信号として得られることになる。

【0021】
(S104:再構成)
センサ部12で受信超音波信号を得ると、その受信超音波信号によって特定される音圧分布p(r,t)から、初期音圧分布p0(r)を再構成する。
本実施形態においては、ここで行う再構成の手法に大きな特徴があるが、この点については詳細を後述する。

【0022】
(S105:光量補正)
初期音圧分布p0(r)を再構成したら、光源部11の仕様から光量Ψの分布が既知であるので、その初期音圧分布p0(r)に対して上記の(1)式に基づき光量Ψを補正して、被検体1の内部における吸収係数の分布像を生成する。

【0023】
このようにして得られる吸収係数分布像は、必要に応じて、いわゆる光音響画像として出力されることになる。光音響画像は、二次元画像または三次元画像の別を問わないが、その一例として最大値投影画像を含むものが知られている。最大値投影画像とは、被検体の画像をMIP(maximum intensity projection:最大値投影法)で表示した、いわゆるMIP画像と呼ばれるものである。

【0024】
なお、ここでは、光量補正を行って光音響画像を出力する場合を例に挙げているが、光量補正を行わずに、初期音圧分布p0(r)を最終的な再構成結果として出力したり、あるいは求めた吸収係数分布を用いて例えば酸素飽和度のような他の指標を求め、吸収係数分布像以外のものを生成して出力したりするようにしてもよい。

【0025】
(光音響イメージングのための装置構成)
続いて、上述した手順の光音響イメージングを行うために必要となる装置(以下「光音響イメージング装置」という。)の構成について説明する。

【0026】
図2は、本実施形態に係る光音響情報処理装置を含む光音響イメージング装置の概略構成の例を模式的に示す説明図である。
図例のように、光音響イメージング装置10は、大別すると、光源部11と、センサ部12と、制御部13と、画像表示部14と、を備えて構成されている。

【0027】
(光源部)
光源部11は、被検体1に照射するレーザ光を発するもので、公知の光源(例えばレーザ光源装置)を用いて構成されたものである。具体的には、例えば、Nd-YAGレーザ、Tiサファイヤレーザ、アレキサンドライトレーザ、GaAs半導体レーザが用いられる。レーザ光のパルス幅としては、例えば100ns以下が好適な条件である。レーザ光のパルス出力としては、照射面積や繰り返し周波数に依存するが、例えば数十マイクロジュール/パルスから100ミリジュール/パルスが好適な範囲である。レーザ光の被検体1への照射強度は、例えば人体の最大許容露光量を超えない範囲で行う。レーザ光のパルス波長は、例えば400nm以上、1600nm以下が好適な範囲である。さらには、生体内において吸収が少ない700nmから1100nmの範囲がより好適な範囲である。

【0028】
(センサ部)
センサ部12は、光源部11のレーザ光照射に応じて吸収体2で発生した音響波(以下「光音響波」ともいう。)が被検体1の体表まで伝搬すると、その伝搬した光音響波(超音波)を受信して計測するもので、その光音響波の検出信号を電気信号として得る圧電変換素子(トランスデューサ)12aを複数有して構成されたものである。なお、センサ部12の具体的な構成(トランスデューサ12aの種類、数、配置等)については、公知技術を利用したものであればよく、ここではその詳細な説明を省略する。

【0029】
(制御部)
制御部13は、光音響イメージング装置10の全体の動作制御を行うものであり、CPU(Central Processing Unit)、ROM(Read Only Memory)、RAM(Random Access Memory)等の組み合わせからなる演算部、フラッシュメモリやHDD(Hard Disk Drive)等の記憶部、外部インタフェース等のデータ入出力部といったハードウエア資源を備えて構成されたものである。つまり、制御部13は、コンピュータ装置としてのハードウエア資源を備えて構成されており、演算部が記憶部に記憶されたプログラムを実行することにより、そのプログラム(ソフトウエア)とハードウエア資源とが協働して、光音響イメージング装置10の動作を制御するようになっている。

【0030】
また、制御部13は、演算部が記憶部に記憶された所定プログラムを実行することにより、信号取得部13a、再構成部13b、結果出力部13cおよび画像記憶部13dとして機能するようになっている。
信号取得部13aは、センサ部12からの光音響波の検出信号を取得するものである。
再構成部13bは、信号取得部13aが取得した検出信号から吸収体2についての音圧分布データを求め、その音圧分布データに基づいて吸収体2おける初期音圧分布p0を再構成して、必要に応じて光音響画像を生成するものである。なお、再構成部13bで行う再構成については詳細を後述する。
結果出力部13cは、光音響画像に代表される再構成部13bでの処理結果を外部装置、具体的には画像表示部14へ出力して、その画像表示部14に表示出力させるものである。
画像記憶部13dは、再構成部13bが生成した光音響画像を必要に応じて記憶保持するものである。なお、画像記憶部13dが光音響画像を記憶保持する際には、そのデータ量を削減する画像圧縮を経ているものとする。また、画像記憶部13dが記憶保持する光音響画像は、結果出力部13cで未出力のものであってもよいし、既に出力されたものであってもよい。

【0031】
これらの各部13a~13dを備える制御部13は、光音響波の検出信号から初期音圧分布p0を再構成する光音響情報処理装置として機能することになる。つまり、制御部13は、本発明に係る「光音響情報処理装置」の一実施形態に相当する。

【0032】
このような制御部13において、上述した各部13a~13dとしての機能は、記憶部に記憶された所定プログラムを演算部が実行することによって実現される。つまり、上述した各部13a~13dとしての機能を実現する所定プログラムは、本発明に係る「光音響情報処理プログラム」の一実施形態に相当する。
その場合に、光音響情報処理プログラムは、コンピュータ装置としての制御部13にインストール可能なものであれば、当該コンピュータ装置で読み取り可能な記録媒体(例えば、磁気ディスク、光ディスク、光磁気ディスク、半導体メモリ等)に格納されて提供されるものであってもよいし、インターネットや専用回線等のネットワークを通じて外部から提供されるものであってもよい。

【0033】
(画像表示部)
画像表示部14は、画像表示を行うディスプレイパネルを備えて構成されたもので、制御部13が生成した光音響画像を含む各種情報の表示出力を行うものである。

【0034】
(2)再構成の手法の詳細
次に、光音響画像を得るための初期音圧分布の再構成の手法について詳しく説明する。なお、以下に説明する処理は、主として制御部13において実行される。

【0035】
(逆問題→順問題)
再構成は、検出点で得られた音響波p(r,t)から初期音圧分布p0(r)を導き出すことであり、数学的には逆問題と呼ばれる。

【0036】
例えば、初期音圧分布p0(r)のベクトルP0から音響波p(r,t)のベクトルPへの伝達関数をΦとすると、ベクトルP0とベクトルPは、P=ΦP0の関係となる。その場合に、再構成にあたっては、逆行列Φ-1が存在すれば、その逆行列Φ-1を用いて逆問題を解くことで、ベクトルPからベクトルP0を求めることができる。

【0037】
このような逆問題を直接解く手法としては、例えば、バックプロジェクション(Back Projection)法が知られている。しかしながら、バックプロジェクション法のような逆問題を直接解く手法では、ベクトルP0を正しく再構成するために、その再構成領域を囲む全周について検出信号を稠密に得る必要があり、また十分に広い帯域の検出信号を得る必要があり、このような検出信号が得られないと正しく再構成できず、再構成アーチファクトが発生してしまうおそれがある。

【0038】
このことから、再構成にあたっては、逆問題を直接解く代わりに、例えば逐次近似法のように、解を設定して順問題で音響波の推定値pk(r,t)を推定し、これを誤差が所定の基準値以下に小さくなるまで繰り返す、といった手法が用いられることがある。

【0039】
(順問題の繰り返し推定法を用いる場合の一般的な手法)
ここで、再構成に順問題の繰り返し推定法を用いる場合の一般的な手法を説明する。
図3は、順問題による繰り返し推定法を用いて再構成を行う場合の一般的な手法の概要を例示する説明図である。

【0040】
再構成を行う場合には、これに先立ち、被検体1へのレーザ光の照射に応じて(S101)、被検体1内の吸収体2から光音響波(初期音圧p0(r))が放射され(S102)、被検体1の体表まで伝搬した音響波p(r,t)の計測値が計測される(S103)。

【0041】
その一方で、図例の手法では、初期音圧分布の設定値p0k-1(r)について、k=1のときの初期値p00(r)を用意する(S201)。k=1のときの初期値p00(r)は、その値が所定値として予め設定されたもの(すなわち固定値であるもの)を用いることが考えられる。ただし、これに限定されることはなく、必要に応じて逐次的に発生または修正させたもの(すなわち可変値であるもの)を初期値として用いてもよいし、予め準備された複数の中から選択されたもの(すなわち選択値であるもの)を初期値として用いてもよいし、別の手法(例えばバックプロジェクション法)によって求めた値を用いることも可能である。そして、初期値p00(r)を用意したら、その初期値p00(r)について、初期音圧分布の設定値p0k-1(r)のベクトルP0kからPk=ΦP0kの順問題を解いてベクトルPkを求め、これによりセンサ部12で検出されたであろう音響波の推定値pk(r,t)を推定する(S202~S204)。

【0042】
音響波の推定値pk(r,t)を推定したら、これを音響波p(r,t)の計測値と比較して(S205)、それぞれの間の誤差が予め設定されている基準値を超えているか否かを判断する(S206)。その結果、誤差が基準値を超えていなければ、音響波の推定値pk(r,t)が音響波p(r,t)の計測値であるものと見做す。そして、その推定値pk(r,t)を導き出したときの初期音圧分布の設定値p0k-1(r)のベクトルP0kを出力して処理を終了する。

【0043】
また、誤差が基準値を超えている場合には、k=k+1にインクリメントした上で(S207)、誤差が小さくなるように初期音圧分布の設定値p0k-1(r)を修正する(S202)。具体的には、例えば、初期値p00(r)から推定した音響波の推定値pk(r,t)についての誤差が基準値を超えている場合であれば、その誤差が小さくなるように、初期値p00(r)を基にしつつ次に用いるべき設定値p0k-1(r)を計算して求める。つまり、修正の都度、誤差が小さくなるように、設定値p0k-1(r)を更新するのである。そして、その更新後の設定値p0k-1(r)について、再び上述した一連の処理を繰り返す(S203~S207)。これを処理が終了するまで、すなわち誤差が基準値を超えなくなるまで、繰り返し行う。

【0044】
以上のような手順を経ることで、逆問題を直接解くことを必要とせずに、初期音圧分布の設定値p0k-1(r)から音響波の推定値pk(r,t)を得る順問題を解くことにより、初期音圧分布についてのベクトルP0kを求めることが可能となる。

【0045】
ただし、上述した手順の一般的な手法では、例えば、Pk=ΦP0kの順問題を解く際の伝達関数Φがm×n行列であると、m<nのときには、他の制約条件を入れない限り、不良設定(ill-posed)問題となってしまい、解が不定となって一意に求まらない。つまり、再構成領域から得られる検出信号が少ない場合には、真の信号を推定することができないおそれがある。

【0046】
(本実施形態における再構成の手法の基本的な手順)
以上のことから、本実施形態においては、順問題を解いて音響波の推定値pk(r,t)を得るのにあたり、上述した手順の一般的な手法ではなく、信号のスパース性を利用した圧縮センシング(Compressed Sensing)を適用する。圧縮センシングとは、スパース性(零成分が多いという性質)を持つ高次元ベクトルで表される信号を少ない観測から復元する技術である。

【0047】
図4は、本実施形態に係る光音響情報処理方法の基本的な手順を例示する説明図であり、圧縮センシングを適用して順問題を解いて再構成を行う手法の概要を示す図である。

【0048】
図例のように、圧縮センシングを適用する場合においても、被検体1へのレーザ光の照射に応じて(S101)、被検体1内の吸収体2から光音響波(初期音圧p0(r))が放射され(S102)、被検体1の体表まで伝搬した音響波p(r,t)の計測値が計測される(S103)。

【0049】
その一方で、図例の手法では、予めスパースな信号sk-1(r)について、k=1のときの初期値s0(r)を用意する(S301)。この初期値s0(r)については、上述した初期値p00(r)の場合と同様に、予め設定された固定値を用いることが考えられるが、これに限定されることはなく、可変値、選択値、別の手法で求めた値等を用いることも可能である。そして、初期値s0(r)を用意したら、その初期値s0(r)について、スパースな信号sk(r)についての信号ベクトルSkを、予め定められた所定の基底Ψで変換し、初期音圧分布の設定値p0k(r)のベクトルP0kを得る(S302~S303)。つまり、P0k=ΨSkの関係に基づき、信号ベクトルSkをある基底Ψで変換してベクトルP0kを求める。ベクトルP0kを求めたら、次いで、そのベクトルP0kからPk=ΦP0kの順問題を解いてベクトルPkを求め、これによりセンサ部12で検出されたであろう音響波の推定値pk(r,t)を推定する(S303~S304)。

【0050】
このことは、伝達関数Φと基底Ψとの積を新たな伝達関数Θと考えてΦΨ=Θとすると、以下の(3)式が成り立つことを意味する。

【0051】
【数3】
JP0006739099B2_000004t.gif

【0052】
つまり、上述した各ステップ(S302~S304)は、圧縮センシングを適用することで、Pk=ΦP0kの順問題を解く代わりに、上記の(3)式に基づき、Pk=ΘSkの関係を満たす信号ベクトルSkで最もスパースなものを解とし、その解からベクトルPkを求めることに相当する(S305)。

【0053】
その後は、音響波の推定値pk(r,t)を推定したら、これを音響波p(r,t)の計測値と比較して(S306)、それぞれの間の誤差が予め設定されている基準値を超えているか否かを判断する(S307)。その結果、誤差が基準値を超えていなければ、音響波の推定値pk(r,t)が音響波p(r,t)の計測値であるものと見做す。そして、その推定値pk(r,t)を導き出したときの初期音圧分布の設定値p0k-1(r)のベクトルP0kを出力して処理を終了する。

【0054】
また、誤差が基準値を超えている場合には、k=k+1にインクリメントした上で(S308)、誤差が小さくなるようにスパースな信号sk-1(r)を修正する。具体的には、例えば、初期値s0(r)から推定した音響波の推定値pk(r,t)についての誤差が基準値を超えている場合であれば、その誤差が小さくなるように、初期値s0(r)を基にしつつ次に用いるべきスパースな信号sk-1(r)を計算して求める。つまり、修正の都度、誤差が小さくなるように、スパースな信号sk-1(r)を更新するのである。そして、その更新後のスパースな信号sk-1(r)について、再び上述した一連の処理を繰り返す(S303~S308)。これを処理が終了するまで、すなわち誤差が基準値を超えなくなるまで、繰り返し行う。

【0055】
以上のような手順を経ることで、圧縮センシングを適用しつつ、初期音圧分布についてのベクトルP0kを求めることが可能となる。しかも、圧縮センシングを適用することで、信号のスパース性を利用することになるので、例えば、Pk=ΦP0kの順問題を解く際の伝達関数Φがm×n行列である場合にm<nであっても、解を求めることが可能となる。つまり、再構成領域から得られる検出信号が少ない場合であっても、真の信号を推定することが可能となる。

【0056】
図5は、圧縮センシングの意義を具体的に説明するための説明図である。
例えば、図5(a)に示すように、P=ΦP0の関係に基づいて、ベクトルPからベクトルP0を求める場合には、伝達関数Φがm×n行列であるときにm<nであると、他の制約条件がないと不良設定問題となってしまい、解が求まらない。
これに対して、圧縮センシング(制約条件としてL1最適化)を適用した場合には、例えば、図5(b)に示すように、P0=ΨSの関係を満たすスパースな信号ベクトルSを用いることで、上記の(3)式の関係が成り立つ。ここで、スパースな信号ベクトルSは、K個の要素以外は「0」であると考えられる。したがって、P=ΘSの関係に基づいてベクトルPから信号ベクトルSを求める場合には、m>Kであれば良設定(well-posed)問題となり、解を求めることができる。つまり、P0=ΨSの関係を利用することで、ベクトルP0を求めることができる。

【0057】
このような圧縮センシングを適用することで、再構成領域から得られる検出信号が少ない場合であっても、Pk=ΦP0kの順問題を解くことが可能となる。したがって、逆問題を直接解きベクトルP0を求めて再構成を行う場合に比べると、再構成アーチファクトの発生を抑制することができ、その結果として再構成した初期音圧分布に基づいて生成する光音響画像の画質劣化を抑制して画質向上が図れるようになる。なお、生成された光音響画像は、例えば、結果出力部13cが画像表示部14へ出力して、その画像表示部14にて表示出力されたり、あるいは画像記憶部13dにより記憶保持されたりすることになる。

【0058】
(処理の高速化の概要)
ところで、圧縮センシングでは、上記の(3)式に基づく信号ベクトルSからベクトルPへの生成(順問題)の過程で、伝達関数(行列)Θについての膨大な計算量を必要とする。通常、圧縮センシングを適用する場合には、信号ベクトルSからベクトルPへの生成を、計測信号の空間領域(実空間)で処理することが一般的である。そのため、下記の(4)式に示すように、被検体1内の全ての音源で生成され、センサ部12における全てのトランスデューサ12aに到達した波形を、インパルス応答として重ね合わせる方式で計算することになり、その計算に膨大な時間がかかってしまう。

【0059】
【数4】
JP0006739099B2_000005t.gif

【0060】
ここで、iは虚数単位、cは音速、kはk=ω/cで表されるもので、ωは角周波数、g(k)は、トランスデューサ12aのインパルス応答、A(r´)は音源位置r´の吸収係数、|r-r´|はトランスデューサ位置rと音源位置r´との距離である。

【0061】
上記の(4)式に基づいて計算を行う場合には、例えば、周波数帯域の分割数をN、トランスデューサ12aの数をM、被検体1内の音源の座標(x,y,z)点数をI×J×Kとすると、その計算回数が積=MNIJK回、和=MNIJK回となり、膨大なものとなってしまう。

【0062】
そこで、本実施形態においては、計算量を削減して処理の高速化を図り、その結果として圧縮センシングを適用した再構成を迅速かつ効率的に行い得るようにするために、以下に述べる(i)および(ii)の特徴点を備えている。

【0063】
すなわち、本実施形態では、(i)k空間(k-space)でのデータを用いて計算を行うことで、計算量を削減して処理の高速化を図っている。
k空間とは、空間周波数を表す空間のことであり、実空間とは互いにフーリエ変換の関係にある空間のことである。

【0064】
さらに、本実施形態では、(ii)k空間での計算を行うことを前提としたうえで、圧縮センシングにおける復元計算を離散条件の球殻積分により行う。
圧縮センシングにおける復元計算とは、トランスデューサ12aによる受信波のスペクトル(時間tについてのフーリエ変換)を求めて復元するために必要となる積分計算のことである。また、離散条件の球殻積分とは、k空間に存在する球殻の表面上の離散データのみを用いて積分計算をすることである。

【0065】
具体的には、上述した(i)および(ii)の特徴点を備えることで、本実施形態では、以下の(5)式および(6)式に基づく計算を行う。
すなわち、本実施形態においては、空間を三次元フーリエ変換したk空間で表現することで、計測ベクトルPのk空間表現である(5)式と、初期音圧p0(r)の空間座標rに関する三次元フーリエ変換(k空間表現)である(6)式が得られる。

【0066】
【数5】
JP0006739099B2_000006t.gif

【0067】
【数6】
JP0006739099B2_000007t.gif

【0068】
なお、上記の(5)式において、Sは、k空間での半径|ω|/cの球殻上の積分を意味している。

【0069】
このように、上記の(5)式によれば、k空間での球殻積分を行うので、二重積分(球殻上積分)となっており、上記の(4)式のような実空間での三重積分を必要とする場合に比べると、計算量を削減して処理の高速化を図ることが可能となる。

【0070】
(高速化した場合の具体的な手順)
ここで、上述した(i)および(ii)の特徴点を備えつつ初期音圧分布の再構成を行う場合の具体的な手順、すなわち本実施形態における再構成の手法の具体的な手順について説明する。

【0071】
本実施形態における具体的な手順は、図4を用いて説明した基本的な手順に当て嵌めると、以下に述べる(I)~(VIII)の各手順を経る。

【0072】
(I)センサ部12で受信された音響波p(r,t)の計測値を信号取得部13aが取得する(S101~S103)。

【0073】
(II)信号取得部13aが取得したp(r,t)の計測値について、再構成部13bが時間tについてフーリエ変換し、p(r,ω)を求める。これにより、被検体1における音圧分布に関する計測値データとして、k空間上のデータを用いることになる。

【0074】
(III)再構成部13bが初期値s0(r)を設定する(S301~S302)。

【0075】
(IV)スパースな信号sk(r)の初期値s0(r)について、P0k=ΨSkの関係に基づき、再構成部13bが基底Ψで変換して、ベクトルP0k(r)を得る(S302~S303)。このとき、基底Ψは、予め定められたものであれば、適切なものを適宜選択して用いることが考えられる。具体的には、例えば、カーブレット基底、フーリエ基底やウェーブレット基底等のいずれかを用いることが考えられる。

【0076】
(V)ベクトルP0k(r)について、再構成部13bが上記の(6)式に基づき、空間を三次元フーリエ変換し、k空間表現のベクトルP0k(k)を求める。これにより、被検体1における音圧分布に関するベクトルデータとして、k空間上のデータを用いることになる。

【0077】
(VI)k空間表現のベクトルP0k(k)から、k空間での演算式である上記の(5)式を用いて、再構成部13bが受信波のスペクトル(時間tについてのフーリエ変換)を求める。このとき、再構成部13bは、上記の(5)式を用いることで、k空間での離散条件の球殻積分を行うことになる。その結果、再構成部13bは、Pk=ΦP0kの順問題を解いてベクトルPkを求め、これによりセンサ部12で検出されたであろう音響波の推定値pk(r,ω)を推定することになる(S303~S305)。

【0078】
(VII)その後、再構成部13bは、上記(II)で求めたp(r,ω)と上記(VI)で求めたpk(r,ω)とを比較し(S306)、これらの間の誤差が基準値を超えている場合にはsk(r)を修正してsk+1(r)を求める(S307~S308)。つまり、k→k+1として、上記(IV)以降の処理を繰り返し行うのである。
このとき、再構成部13bは、例えば、以下の(7)式によって特定される公知のL1最適化アルゴリズムを用いることで、上述した繰り返し処理を効率的に収束させるようにしてもよい。

【0079】
【数7】
JP0006739099B2_000008t.gif

【0080】
(VIII)p(r,ω)とpk(r,ω)との誤差が基準値以下の場合、再構成部13bは、そのpk(r,ω)を導き出したときのスパースな信号sk(r)を結果出力部13cに出力させて処理を終了する。

【0081】
(高速化した計算量の具体例)
本実施形態によれば、以上に説明したような処理を行うので、計算量を削減して処理の高速化を図ることができる。

【0082】
具体的には、例えば、周波数帯域の分割数をN、トランスデューサ12aの数をM、被検体1内の音源の座標点数をI×J×Kとすると、上記の(6)式による計算回数がIJK・log(IJK)回、上記の(5)式による計算回数が球殻積分で最大でもIJK回となり、本実施形態の全体ではIJK(M+log(IJK))回となる。
これに対して、本実施形態のようなk空間ではなく実空間において上記の(4)式に基づいて計算を行う場合には、既に説明したように、計算回数がMNIJK回となってしまう。
したがって、例えば、M,Nがそれぞれ100回、I,J,Kがそれぞれ1000回であると、実空間における上記(4)式の計算ではMNIJK=1002×10003=1013回の計算回数が必要であるのに対して、本実施形態ではIJK(M+log(IJK))=10003(100+log(10003)=109×130≒1011回となる。つまり、本実施形態によれば、凡そ二桁程度の計算回数の削減による高速化が実現可能となる。

【0083】
(球殻積分の具体的な一態様)
上述したように、本実施形態では、k空間での計算を行うので、圧縮センシングにおける復元計算を球殻積分によって行うことが可能となる。そして、その球殻積分にあたり、k空間に存在する球殻の表面上の離散データのみを用いて積分計算をすればよい。
この球殻積分については、以下に述べる態様で行えば、積分計算によって得られるデータの精度向上を図るうえで非常に好ましいものとなる。
以下に、本実施形態における球殻積分の好適な一態様について具体的に説明する。なお、以下に説明する例は、球殻積分の好適な一態様に過ぎず、本発明がこれに限定されるものではない。

【0084】
図6は、本実施形態に係る光音響情報処理方法における球殻積分の具体的な一態様の概要を具体的に説明するための説明図である。
球殻積分にあたっては、例えば、kx,ky,kzの各軸を有する座標空間(ただし、kz軸は不図示)において、k=ωcで表される球殻を積分計算によって求める。

【0085】
その場合において、積分計算に際しては、先ず、計算対象となる領域を{±kx>|ky|,|kz|}、{±ky>|kx|,|kz|}、{±kz>|kx|,|ky|}の6領域に分割する。これにより、信号ベクトルSについては、以下の(8)式および(9)式が成り立つことになる。

【0086】
【数8】
JP0006739099B2_000009t.gif

【0087】
【数9】
JP0006739099B2_000010t.gif

【0088】
ここで、例えば、Sx={kx>|ky|,|kz|}の場合に着目すると、この領域については、ky,kzを積分変数として積分し、全てのグリッド上のky0,kz0について以下の(10)式のkxを求める。

【0089】
【数10】
JP0006739099B2_000011t.gif

【0090】
そして、kxの最近傍のグリッド2点kx0,kx1を線形補間して、p(kx,ky0,kz0)を求める。

【0091】
このような処理を分割した各領域の全てについて行うことで、半径|ω|/cの球殻について、その表面上の全領域についての積分結果を得ることができるようになる。

【0092】
つまり、ここで例に挙げて説明した球殻積分の一態様では、その球殻積分にあたり、球殻を予め設定された数の複数領域(例えば6領域)に分割するとともに、予め設定された選択基準に従いつつ各領域で積分変数をそれぞれ選択する。そして、それぞれ選択した積分変数を各領域別に適用しつつ、球殻の全域についての積分計算を行うのである。

【0093】
このような球殻積分の一態様によれば、図6に示す例において二次元に簡略化して考えると、例えば、領域Aについてはky軸に沿った積分変数を選択し、領域Bについてはkx軸に沿った積分変数を選択する、といったことが実現可能となる。したがって、各領域について、変動の大きな軸方向に沿って積分変数を選択し、その選択した積分変数を適用して積分計算を行うことができる。つまり、全領域に対して積分変数を一律に適用する場合に比べると、球殻の形状の特質を適切に考慮した積分計算を行うことができるので、球殻積分を行う場合であっても、その積分計算によって得られるデータの精度向上が図れる。その結果として、再構成した初期音圧分布に基づいて生成する光音響画像の画質向上が期待できる。

【0094】
(光音響画像の記憶保持)
以上に説明した手順を経ることで、本実施形態では、被検体1における吸収体2について、初期音圧分布p0(r)の再構成を行う。そして、初期音圧分布p0(r)を再構成再構成したら、必要に応じて、光量補正を行って吸収係数分布像(すなわち光音響画像)を生成する。このようにして生成された光音響画像は、結果出力部13cから画像表示部14へ出力されたり、あるいは画像記憶部13dで記憶保持されたりする。

【0095】
画像記憶部13dでの記憶保持にあたっては、光音響画像の画像データ量が膨大であることが多いため、画像データのデータ圧縮を行うことが考えられる。データ圧縮の手法については、特に限定されることはなく、公知のデータ圧縮手法を用いればよい。

【0096】
その場合において、本実施形態で得られる光音響画像は、上述したように圧縮センシングを適用した再構成を経て生成されているため、例えばバックプロジェクション法を適用した場合に比べると、ぼけやストリークアーチファクト等の再構成アーチファクトが発生してしまうことがない。そのため、光音響画像についてデータ圧縮を行うと、再構成アーチファクトが少ない分だけ、例えばバックプロジェクション法を適用した場合に比べて、画像圧縮率が高くなる。

【0097】
つまり、本実施形態によれば、光音響画像についてデータ圧縮を行う際の画像圧縮率を高くすることができ、画像記憶部13dのために必要とするデータ記憶容量が増大してしまうのを抑制することが可能となる。

【0098】
(3)本実施形態の効果
本実施形態によれば、以下に示す一つまたは複数の効果を奏する。

【0099】
(a)本実施形態では、初期音圧分布の再構成にあたり、逆問題を直接解く代わりに、順問題を解いて音響波の推定値pk(r,t)を推定し、これを基準値との誤差が小さくなるまで繰り返す。したがって、逆問題を直接解く場合とは異なり、ぼけやストリークアーチファクト等の再構成アーチファクトが発生してしまうことがなく、再構成結果に基づいて光音響画像を生成する際の画質劣化を抑制することができる。

【0100】
(b)また、本実施形態では、順問題を解いて音響波の推定値pk(r,t)を推定するのにあたり、信号のスパース性を利用した圧縮センシングを適用する。したがって、信号のスパース性を利用することになるので、不良設定問題とはならずに解を求めることが可能となる。つまり、再構成領域から得られる検出信号が少ない場合であっても、真の信号を推定することが可能となる。このことは、光音響イメージングに適用して非常に好適であることを意味する。光音響イメージングは、再構成領域を囲む全周について検出信号を稠密に得ることが容易ではなく、また十分に広い帯域の検出信号が得されるとは限らないからである。

【0101】
(c)また、本実施形態では、初期音圧分布の再構成にあたり、音圧分布データ(具体的には、音圧分布に関する計測値データや音圧分布に関するベクトルデータ等)としてk空間上のデータを用い、k空間での計算を行うことを前提としたうえで、圧縮センシングにおける復元計算を離散条件の球殻積分により行う。つまり、本実施形態では、k空間での球殻積分を行う。したがって、本実施形態によれば、例えば実空間での積分計算を必要とする場合に比べると、計算量を削減して処理の高速化を図ることが可能となる。具体的には、凡そ二桁程度の計算回数の削減による高速化が実現可能となる。このように、初期音圧分布の再構成に必要となる計算処理の高速化を図ることで、本実施形態では、初期音圧分布の再構成を迅速かつ効率的に行うことができる。

【0102】
(d)また、本実施形態では、球殻積分にあたり、球殻を複数領域に分割するとともに、各領域で積分変数をそれぞれ選択する。したがって、全領域に対して積分変数を一律に適用する場合に比べると、球殻の形状の特質を適切に考慮した積分計算を行うことができるので、球殻積分を行う場合であっても、その積分計算によって得られるデータの精度向上が図れ、その結果として、再構成結果に基づいて光音響画像を生成する際の画質向上が期待できる。

【0103】
(e)また、本実施形態では、再構成結果に基づいて生成した光音響画像の記憶保持に際し、その光音響画像の再構成アーチファクトが少ないことから、その光音響画像についてデータ圧縮を行う際の画像圧縮率を高くすることができる。

【0104】
<他の実施形態>
以上に、本発明の一実施形態を具体的に説明したが、本発明は上述の実施形態に限定されるものではなく、その要旨を逸脱しない範囲で種々変更可能である。

【0105】
例えば、本実施形態では、主として、被検体1の内部における血液等の吸収体2について、その吸収係数分布像を得るための光音響イメージングを例に挙げたが、本発明がこれに限定されることはない。すなわち、本発明は、血液等の吸収体以外の光音響イメージングにも適用することが可能である。

【0106】
また、例えば、本実施形態では、光音響イメージング装置10を構成する各部11~14、および、制御部13を構成する各部13a~13bについて具体的に説明したが、上述した実施形態で例に挙げた具体例に限定されることはなく、同等のものに代替して構成しても構わない。具体的には、その一例として、画像表示部14に代わって、画像表示機能を備えた外部装置に有線または無線で接続するインタフェース部を設ける、といったことが考えられる。

【0107】
また、例えば、本実施形態では、画像の再構成にあたり、上記の(1)式~(10)式を基にする場合を例に挙げたが、各式はその要旨を変更しない範囲で適宜修正しても構わない。
【符号の説明】
【0108】
1…被検体、2…吸収体、10…光音響イメージング装置、11…光源部、12…センサ部、12a…圧電変換素子(トランスデューサ)、13…制御部(光音響情報処理装置)、13a…信号取得部、13b…再構成部、13c…結果出力部、13d…画像記憶部、14…画像表示部
Drawing
(In Japanese)【図1】
0
(In Japanese)【図2】
1
(In Japanese)【図3】
2
(In Japanese)【図4】
3
(In Japanese)【図5】
4
(In Japanese)【図6】
5