JP6739099B2 - 光音響情報処理方法、光音響情報処理プログラムおよび光音響情報処理装置 - Google Patents

Description

本発明は、光音響情報処理方法、光音響情報処理プログラムおよび光音響情報処理装置に関する。

近年、例えば医療分野において、光音響効果を用いて被検体内部を画像化する光超音波イメージング法（photoacoustic imaging）を利用することが提案されている。光音響効果とは、被検体にレーザ光を照射すると、被検体内部における血液等の光吸収体に吸収され、局所的な温度上昇が熱弾性変形を引き起こし、その熱弾性変形による体積膨張で音響波（超音波）が発生する現象である。光超音波イメージング法では、発生した音響波を受信して、その音源（すなわち、被検体内部における光吸収体）の初期音圧分布を再構成することで、被検体内部を画像化する（例えば、非特許文献１参照）。

椎名毅、「光音響を用いた新しい超音波診断装置」、メディカル＆イメージングＮｏ．３、オプトロニクス社、２０１５年、ｐ．５３−５９

光超音波イメージング法の医療適用にあたっては、画像化のために必要となる初期音圧分布の再構成を、画質劣化を抑制しつつ、迅速かつ効率的に行うことが求められる。

本発明は、画質劣化を抑制しつつ、迅速かつ効率的に、初期音圧分布の再構成を行うことができる光音響情報処理方法、光音響情報処理プログラムおよび光音響情報処理装置を提供することを目的とする。

本発明の一態様によれば、
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を備え、
前記音圧分布データとしてｋ空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理方法が提供される。

本発明の他の一態様によれば、
コンピュータに、
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を実行させるとともに、
前記音圧分布データとしてｋ空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理プログラムが提供される。

本発明のさらに他の一態様によれば、
光音響波の検出信号を取得する信号取得部と、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成する再構成部と、
前記再構成部での処理結果を出力する結果出力部と、を備え、
前記音圧分布データとしてｋ空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理装置が提供される。

本発明によれば、画質劣化を抑制しつつ、迅速かつ効率的に、初期音圧分布の再構成を行うことができる。

光音響イメージングの基本的な手順を例示する説明図である。 本発明の一実施形態に係る光音響情報処理装置を含む光音響イメージング装置の概略構成の例を模式的に示す説明図である。 順問題による繰り返し推定法を用いて再構成を行う場合の一般的な手法の概要を例示する説明図である。 本発明の一実施形態に係る光音響情報処理方法の基本的な手順を例示する説明図であり、圧縮センシングを適用して順問題を解いて再構成を行う手法の概要を示す図である。 圧縮センシングの意義を具体的に説明するための説明図である。 本発明の一実施形態に係る光音響情報処理方法における球殻積分の具体的な一態様の概要を具体的に説明するための説明図である。

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

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

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

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

（Ｓ１０２：音響波放射）
レーザ光を照射すると、被検体１の内部では、血液等の吸収体２が光を吸収して、局所的な温度上昇が熱弾性変形を引き起こし、その熱弾性変形による体積膨張で吸収体２が音響波（初期音圧ｐ_０（ｒ）)を放射する。このとき、光の吸収により発生する音響波（超音波）の大きさは、以下の（１）式で表される。

ここで、βは熱膨張係数、ｃは音速、Ｃ_ｐは定圧比熱、μ_ａは光の吸収係数、Ψは光量である。係数Ｇ＝βｃ^２／Ｃ_ｐは、光エネルギーから音圧への変換係数でグリューナイゼン（Gruneisen）係数と呼ばれる。Ｇは軟組織間であまり差はないため、初期音圧ｐ_０は、組織の吸収係数μ_ａと光量Ψの積に比例することになる。

また、吸収体２で発生した位置ｒ、時刻ｔにおける音響波ｐ（ｒ，ｔ）は、以下の（２）式で表わされる光音響波動方程式に従って、被検体１の内部を伝搬する。

（Ｓ１０３：音響波計測）
その後は、吸収体２で発生し被検体１の体表まで伝搬した音響波ｐ（ｒ，ｔ）を、体表近傍に配されて複数の圧電変換素子（トランスデューサ）を有するセンサ部１２を用いて計測する。これにより、センサ部１２では、吸収体２からの音響波ｐ（ｒ，ｔ）の検出信号（受信超音波信号）が電気信号として得られることになる。

（Ｓ１０４：再構成）
センサ部１２で受信超音波信号を得ると、その受信超音波信号によって特定される音圧分布ｐ（ｒ，ｔ）から、初期音圧分布ｐ_０（ｒ）を再構成する。
本実施形態においては、ここで行う再構成の手法に大きな特徴があるが、この点については詳細を後述する。

（Ｓ１０５：光量補正）
初期音圧分布ｐ_０（ｒ）を再構成したら、光源部１１の仕様から光量Ψの分布が既知であるので、その初期音圧分布ｐ_０（ｒ）に対して上記の（１）式に基づき光量Ψを補正して、被検体１の内部における吸収係数の分布像を生成する。

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

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

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

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

（光源部）
光源部１１は、被検体１に照射するレーザ光を発するもので、公知の光源（例えばレーザ光源装置）を用いて構成されたものである。具体的には、例えば、Ｎｄ−ＹＡＧレーザ、Ｔｉサファイヤレーザ、アレキサンドライトレーザ、ＧａＡｓ半導体レーザが用いられる。レーザ光のパルス幅としては、例えば１００ｎｓ以下が好適な条件である。レーザ光のパルス出力としては、照射面積や繰り返し周波数に依存するが、例えば数十マイクロジュール／パルスから１００ミリジュール／パルスが好適な範囲である。レーザ光の被検体１への照射強度は、例えば人体の最大許容露光量を超えない範囲で行う。レーザ光のパルス波長は、例えば４００ｎｍ以上、１６００ｎｍ以下が好適な範囲である。さらには、生体内において吸収が少ない７００ｎｍから１１００ｎｍの範囲がより好適な範囲である。

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

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

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

これらの各部１３ａ〜１３ｄを備える制御部１３は、光音響波の検出信号から初期音圧分布ｐ_０を再構成する光音響情報処理装置として機能することになる。つまり、制御部１３は、本発明に係る「光音響情報処理装置」の一実施形態に相当する。

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

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

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

（逆問題→順問題）
再構成は、検出点で得られた音響波ｐ（ｒ，ｔ）から初期音圧分布ｐ_０（ｒ）を導き出すことであり、数学的には逆問題と呼ばれる。

例えば、初期音圧分布ｐ_０（ｒ）のベクトルＰ_０から音響波ｐ（ｒ，ｔ）のベクトルＰへの伝達関数をΦとすると、ベクトルＰ_０とベクトルＰは、Ｐ＝ΦＰ_０の関係となる。その場合に、再構成にあたっては、逆行列Φ^-1が存在すれば、その逆行列Φ^-1を用いて逆問題を解くことで、ベクトルＰからベクトルＰ_０を求めることができる。

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

このことから、再構成にあたっては、逆問題を直接解く代わりに、例えば逐次近似法のように、解を設定して順問題で音響波の推定値ｐ^ｋ（ｒ，ｔ）を推定し、これを誤差が所定の基準値以下に小さくなるまで繰り返す、といった手法が用いられることがある。

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

再構成を行う場合には、これに先立ち、被検体１へのレーザ光の照射に応じて（Ｓ１０１）、被検体１内の吸収体２から光音響波（初期音圧ｐ_０（ｒ）)が放射され（Ｓ１０２）、被検体１の体表まで伝搬した音響波ｐ（ｒ，ｔ）の計測値が計測される（Ｓ１０３）。

その一方で、図例の手法では、初期音圧分布の設定値ｐ_０ ^k-1（ｒ）について、ｋ＝１のときの初期値ｐ_０ ^０（ｒ）を用意する（Ｓ２０１）。ｋ＝１のときの初期値ｐ_０ ^０（ｒ）は、その値が所定値として予め設定されたもの（すなわち固定値であるもの）を用いることが考えられる。ただし、これに限定されることはなく、必要に応じて逐次的に発生または修正させたもの（すなわち可変値であるもの）を初期値として用いてもよいし、予め準備された複数の中から選択されたもの（すなわち選択値であるもの）を初期値として用いてもよいし、別の手法（例えばバックプロジェクション法）によって求めた値を用いることも可能である。そして、初期値ｐ_０ ^０（ｒ）を用意したら、その初期値ｐ_０ ^０（ｒ）について、初期音圧分布の設定値ｐ_０ ^k-1（ｒ）のベクトルＰ_０ ^kからＰ^ｋ＝ΦＰ_０ ^kの順問題を解いてベクトルＰ^ｋを求め、これによりセンサ部１２で検出されたであろう音響波の推定値ｐ^ｋ（ｒ，ｔ）を推定する（Ｓ２０２〜Ｓ２０４）。

音響波の推定値ｐ^ｋ（ｒ，ｔ）を推定したら、これを音響波ｐ（ｒ，ｔ）の計測値と比較して（Ｓ２０５）、それぞれの間の誤差が予め設定されている基準値を超えているか否かを判断する（Ｓ２０６）。その結果、誤差が基準値を超えていなければ、音響波の推定値ｐ^ｋ（ｒ，ｔ）が音響波ｐ（ｒ，ｔ）の計測値であるものと見做す。そして、その推定値ｐ^ｋ（ｒ，ｔ）を導き出したときの初期音圧分布の設定値ｐ_０ ^k-1（ｒ）のベクトルＰ_０ ^kを出力して処理を終了する。

また、誤差が基準値を超えている場合には、ｋ＝ｋ＋１にインクリメントした上で（Ｓ２０７）、誤差が小さくなるように初期音圧分布の設定値ｐ_０ ^k-1（ｒ）を修正する（Ｓ２０２）。具体的には、例えば、初期値ｐ_０ ^０（ｒ）から推定した音響波の推定値ｐ^ｋ（ｒ，ｔ）についての誤差が基準値を超えている場合であれば、その誤差が小さくなるように、初期値ｐ_０ ^０（ｒ）を基にしつつ次に用いるべき設定値ｐ_０ ^k-1（ｒ）を計算して求める。つまり、修正の都度、誤差が小さくなるように、設定値ｐ_０ ^k-1（ｒ）を更新するのである。そして、その更新後の設定値ｐ_０ ^k-1（ｒ）について、再び上述した一連の処理を繰り返す（Ｓ２０３〜Ｓ２０７）。これを処理が終了するまで、すなわち誤差が基準値を超えなくなるまで、繰り返し行う。

以上のような手順を経ることで、逆問題を直接解くことを必要とせずに、初期音圧分布の設定値ｐ_０ ^k-1（ｒ）から音響波の推定値ｐ^ｋ（ｒ，ｔ）を得る順問題を解くことにより、初期音圧分布についてのベクトルＰ_０ ^kを求めることが可能となる。

ただし、上述した手順の一般的な手法では、例えば、Ｐ^k＝ΦＰ_０ ^kの順問題を解く際の伝達関数Φがｍ×ｎ行列であると、ｍ＜ｎのときには、他の制約条件を入れない限り、不良設定（ill-posed）問題となってしまい、解が不定となって一意に求まらない。つまり、再構成領域から得られる検出信号が少ない場合には、真の信号を推定することができないおそれがある。

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

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

図例のように、圧縮センシングを適用する場合においても、被検体１へのレーザ光の照射に応じて（Ｓ１０１）、被検体１内の吸収体２から光音響波（初期音圧ｐ_０（ｒ）)が放射され（Ｓ１０２）、被検体１の体表まで伝搬した音響波ｐ（ｒ，ｔ）の計測値が計測される（Ｓ１０３）。

その一方で、図例の手法では、予めスパースな信号ｓ^k-1（ｒ）について、ｋ＝１のときの初期値ｓ^０（ｒ）を用意する（Ｓ３０１）。この初期値ｓ^０（ｒ）については、上述した初期値ｐ_０ ^０（ｒ）の場合と同様に、予め設定された固定値を用いることが考えられるが、これに限定されることはなく、可変値、選択値、別の手法で求めた値等を用いることも可能である。そして、初期値ｓ^０（ｒ）を用意したら、その初期値ｓ^０（ｒ）について、スパースな信号ｓ^ｋ（ｒ）についての信号ベクトルＳ^ｋを、予め定められた所定の基底Ψで変換し、初期音圧分布の設定値ｐ_０ ^k（ｒ）のベクトルＰ_０ ^kを得る（Ｓ３０２〜Ｓ３０３）。つまり、Ｐ_０ ^k＝ΨＳ^ｋの関係に基づき、信号ベクトルＳ^ｋをある基底Ψで変換してベクトルＰ_０ ^kを求める。ベクトルＰ_０ ^kを求めたら、次いで、そのベクトルＰ_０ ^kからＰ^ｋ＝ΦＰ_０ ^kの順問題を解いてベクトルＰ^ｋを求め、これによりセンサ部１２で検出されたであろう音響波の推定値ｐ^ｋ（ｒ，ｔ）を推定する（Ｓ３０３〜Ｓ３０４）。

このことは、伝達関数Φと基底Ψとの積を新たな伝達関数Θと考えてΦΨ＝Θとすると、以下の（３）式が成り立つことを意味する。

つまり、上述した各ステップ（Ｓ３０２〜Ｓ３０４）は、圧縮センシングを適用することで、Ｐ^k＝ΦＰ_０ ^kの順問題を解く代わりに、上記の（３）式に基づき、Ｐ^ｋ＝ΘＳ^ｋの関係を満たす信号ベクトルＳ^ｋで最もスパースなものを解とし、その解からベクトルＰ^ｋを求めることに相当する（Ｓ３０５）。

その後は、音響波の推定値ｐ^ｋ（ｒ，ｔ）を推定したら、これを音響波ｐ（ｒ，ｔ）の計測値と比較して（Ｓ３０６）、それぞれの間の誤差が予め設定されている基準値を超えているか否かを判断する（Ｓ３０７）。その結果、誤差が基準値を超えていなければ、音響波の推定値ｐ^ｋ（ｒ，ｔ）が音響波ｐ（ｒ，ｔ）の計測値であるものと見做す。そして、その推定値ｐ^ｋ（ｒ，ｔ）を導き出したときの初期音圧分布の設定値ｐ_０ ^k-1（ｒ）のベクトルＰ_０ ^kを出力して処理を終了する。

また、誤差が基準値を超えている場合には、ｋ＝ｋ＋１にインクリメントした上で（Ｓ３０８）、誤差が小さくなるようにスパースな信号ｓ^k-1（ｒ）を修正する。具体的には、例えば、初期値ｓ^０（ｒ）から推定した音響波の推定値ｐ^ｋ（ｒ，ｔ）についての誤差が基準値を超えている場合であれば、その誤差が小さくなるように、初期値ｓ^０（ｒ）を基にしつつ次に用いるべきスパースな信号ｓ^k-1（ｒ）を計算して求める。つまり、修正の都度、誤差が小さくなるように、スパースな信号ｓ^k-1（ｒ）を更新するのである。そして、その更新後のスパースな信号ｓ^k-1（ｒ）について、再び上述した一連の処理を繰り返す（Ｓ３０３〜Ｓ３０８）。これを処理が終了するまで、すなわち誤差が基準値を超えなくなるまで、繰り返し行う。

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

図５は、圧縮センシングの意義を具体的に説明するための説明図である。
例えば、図５（ａ）に示すように、Ｐ＝ΦＰ_０の関係に基づいて、ベクトルＰからベクトルＰ_０を求める場合には、伝達関数Φがｍ×ｎ行列であるときにｍ＜ｎであると、他の制約条件がないと不良設定問題となってしまい、解が求まらない。
これに対して、圧縮センシング（制約条件としてＬ１最適化）を適用した場合には、例えば、図５（ｂ）に示すように、Ｐ_０＝ΨＳの関係を満たすスパースな信号ベクトルＳを用いることで、上記の（３）式の関係が成り立つ。ここで、スパースな信号ベクトルＳは、Ｋ個の要素以外は「０」であると考えられる。したがって、Ｐ＝ΘＳの関係に基づいてベクトルＰから信号ベクトルＳを求める場合には、ｍ＞Ｋであれば良設定（well-posed）問題となり、解を求めることができる。つまり、Ｐ_０＝ΨＳの関係を利用することで、ベクトルＰ_０を求めることができる。

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

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

ここで、ｉは虚数単位、ｃは音速、ｋはｋ＝ω／ｃで表されるもので、ωは角周波数、ｇ（ｋ）は、トランスデューサ１２ａのインパルス応答、Ａ（ｒ´）は音源位置ｒ´の吸収係数、｜ｒ-ｒ´｜はトランスデューサ位置ｒと音源位置ｒ´との距離である。

上記の（４）式に基づいて計算を行う場合には、例えば、周波数帯域の分割数をＮ、トランスデューサ１２ａの数をＭ、被検体１内の音源の座標（ｘ，ｙ，ｚ）点数をＩ×Ｊ×Ｋとすると、その計算回数が積＝ＭＮＩＪＫ回、和＝ＭＮＩＪＫ回となり、膨大なものとなってしまう。

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

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

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

具体的には、上述した（ｉ）および（ii）の特徴点を備えることで、本実施形態では、以下の（５）式および（６）式に基づく計算を行う。
すなわち、本実施形態においては、空間を三次元フーリエ変換したｋ空間で表現することで、計測ベクトルＰのｋ空間表現である（５）式と、初期音圧ｐ_０（ｒ）の空間座標ｒに関する三次元フーリエ変換（ｋ空間表現）である（６）式が得られる。

なお、上記の（５）式において、Ｓは、ｋ空間での半径｜ω｜／ｃの球殻上の積分を意味している。

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

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

本実施形態における具体的な手順は、図４を用いて説明した基本的な手順に当て嵌めると、以下に述べる（Ｉ）〜（VIII）の各手順を経る。

（Ｉ）センサ部１２で受信された音響波ｐ（ｒ，ｔ）の計測値を信号取得部１３ａが取得する（Ｓ１０１〜Ｓ１０３）。

（II）信号取得部１３ａが取得したｐ（ｒ，ｔ）の計測値について、再構成部１３ｂが時間ｔについてフーリエ変換し、ｐ（ｒ，ω）を求める。これにより、被検体１における音圧分布に関する計測値データとして、ｋ空間上のデータを用いることになる。

（III）再構成部１３ｂが初期値ｓ^０（ｒ）を設定する（Ｓ３０１〜Ｓ３０２）。

（IV）スパースな信号ｓ^ｋ（ｒ）の初期値ｓ^０（ｒ）について、Ｐ_０ ^k＝ΨＳ^ｋの関係に基づき、再構成部１３ｂが基底Ψで変換して、ベクトルＰ_０ ^k（ｒ）を得る（Ｓ３０２〜Ｓ３０３）。このとき、基底Ψは、予め定められたものであれば、適切なものを適宜選択して用いることが考えられる。具体的には、例えば、カーブレット基底、フーリエ基底やウェーブレット基底等のいずれかを用いることが考えられる。

（Ｖ）ベクトルＰ_０ ^k（ｒ）について、再構成部１３ｂが上記の（６）式に基づき、空間を三次元フーリエ変換し、ｋ空間表現のベクトルＰ_０ ^k（ｋ）を求める。これにより、被検体１における音圧分布に関するベクトルデータとして、ｋ空間上のデータを用いることになる。

（VI）ｋ空間表現のベクトルＰ_０ ^k（ｋ）から、ｋ空間での演算式である上記の（５）式を用いて、再構成部１３ｂが受信波のスペクトル（時間ｔについてのフーリエ変換）を求める。このとき、再構成部１３ｂは、上記の（５）式を用いることで、ｋ空間での離散条件の球殻積分を行うことになる。その結果、再構成部１３ｂは、Ｐ^ｋ＝ΦＰ_０ ^kの順問題を解いてベクトルＰ^ｋを求め、これによりセンサ部１２で検出されたであろう音響波の推定値ｐ^ｋ（ｒ，ω）を推定することになる（Ｓ３０３〜Ｓ３０５）。

（VII）その後、再構成部１３ｂは、上記（II）で求めたｐ（ｒ，ω）と上記（VI）で求めたｐ^ｋ（ｒ，ω）とを比較し（Ｓ３０６）、これらの間の誤差が基準値を超えている場合にはｓ^ｋ（ｒ）を修正してｓ^ｋ＋１（ｒ）を求める（Ｓ３０７〜Ｓ３０８）。つまり、ｋ→ｋ＋１として、上記（IV）以降の処理を繰り返し行うのである。
このとき、再構成部１３ｂは、例えば、以下の（７）式によって特定される公知のＬ１最適化アルゴリズムを用いることで、上述した繰り返し処理を効率的に収束させるようにしてもよい。

（VIII）ｐ（ｒ，ω）とｐ^ｋ（ｒ，ω）との誤差が基準値以下の場合、再構成部１３ｂは、そのｐ^ｋ（ｒ，ω）を導き出したときのスパースな信号ｓ^ｋ（ｒ）を結果出力部１３ｃに出力させて処理を終了する。

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

具体的には、例えば、周波数帯域の分割数をＮ、トランスデューサ１２ａの数をＭ、被検体１内の音源の座標点数をＩ×Ｊ×Ｋとすると、上記の（６）式による計算回数がＩＪＫ・log（ＩＪＫ）回、上記の（５）式による計算回数が球殻積分で最大でもＩＪＫ回となり、本実施形態の全体ではＩＪＫ（Ｍ＋log（ＩＪＫ））回となる。
これに対して、本実施形態のようなｋ空間ではなく実空間において上記の（４）式に基づいて計算を行う場合には、既に説明したように、計算回数がＭＮＩＪＫ回となってしまう。
したがって、例えば、Ｍ，Ｎがそれぞれ１００回、Ｉ，Ｊ，Ｋがそれぞれ１０００回であると、実空間における上記（４）式の計算ではＭＮＩＪＫ＝１００^２×１０００^３＝１０¹³回の計算回数が必要であるのに対して、本実施形態ではＩＪＫ（Ｍ＋log（ＩＪＫ））＝１０００^３（１００＋log（１０００^３）＝１０^９×１３０≒１０¹¹回となる。つまり、本実施形態によれば、凡そ二桁程度の計算回数の削減による高速化が実現可能となる。

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

図６は、本実施形態に係る光音響情報処理方法における球殻積分の具体的な一態様の概要を具体的に説明するための説明図である。
球殻積分にあたっては、例えば、ｋｘ，ｋｙ，ｋｚの各軸を有する座標空間（ただし、ｋｚ軸は不図示）において、ｋ＝ωｃで表される球殻を積分計算によって求める。

その場合において、積分計算に際しては、先ず、計算対象となる領域を｛±ｋｘ＞｜ｋｙ｜，｜ｋｚ｜｝、｛±ｋｙ＞｜ｋｘ｜，｜ｋｚ｜｝、｛±ｋｚ＞｜ｋｘ｜，｜ｋｙ｜｝の６領域に分割する。これにより、信号ベクトルＳについては、以下の（８）式および（９）式が成り立つことになる。

ここで、例えば、Ｓ_ｘ ^＋＝｛ｋｘ＞｜ｋｙ｜，｜ｋｚ｜｝の場合に着目すると、この領域については、ｋｙ，ｋｚを積分変数として積分し、全てのグリッド上のｋy0，ｋz0について以下の（１０）式のｋｘを求める。

そして、ｋｘの最近傍のグリッド２点ｋx0，ｋx1を線形補間して、ｐ（ｋｘ，ｋy0，ｋz0）を求める。

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

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

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

（光音響画像の記憶保持）
以上に説明した手順を経ることで、本実施形態では、被検体１における吸収体２について、初期音圧分布ｐ_０（ｒ）の再構成を行う。そして、初期音圧分布ｐ_０（ｒ）を再構成再構成したら、必要に応じて、光量補正を行って吸収係数分布像（すなわち光音響画像）を生成する。このようにして生成された光音響画像は、結果出力部１３ｃから画像表示部１４へ出力されたり、あるいは画像記憶部１３ｄで記憶保持されたりする。

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

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

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

（３）本実施形態の効果
本実施形態によれば、以下に示す一つまたは複数の効果を奏する。

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

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

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

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

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

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

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

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

また、例えば、本実施形態では、画像の再構成にあたり、上記の（１）式〜（１０）式を基にする場合を例に挙げたが、各式はその要旨を変更しない範囲で適宜修正しても構わない。

１…被検体、２…吸収体、１０…光音響イメージング装置、１１…光源部、１２…センサ部、１２ａ…圧電変換素子（トランスデューサ）、１３…制御部（光音響情報処理装置）、１３ａ…信号取得部、１３ｂ…再構成部、１３ｃ…結果出力部、１３ｄ…画像記憶部、１４…画像表示部

Claims (5)

1. 光音響波の検出信号を取得するステップと、
   前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を備え、
   前記音圧分布データとしてｋ空間上のデータを用い、
   前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
   光音響情報処理方法。
2. 前記球殻積分にあたり、球殻を複数領域に分割するとともに、各領域で積分変数をそれぞれ選択する
   請求項１に記載の光音響情報処理方法。
3. 前記初期音圧分布の再構成の結果に基づいて光音響画像を生成するとともに、生成した前記光音響画像についての画像データをデータ圧縮して記憶保持する
   請求項１または２に記載の光音響情報処理方法。
4. コンピュータに、
   光音響波の検出信号を取得するステップと、
   前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を実行させるとともに、
   前記音圧分布データとしてｋ空間上のデータを用い、
   前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
   光音響情報処理プログラム。
5. 光音響波の検出信号を取得する信号取得部と、
   前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成する再構成部と、
   前記再構成部での処理結果を出力する結果出力部と、を備え、
   前記音圧分布データとしてｋ空間上のデータを用い、
   前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
   光音響情報処理装置。