(57)【要約】【課題】光音響イメージングにおいて、画質劣化を抑制しつつ、迅速かつ効率的に、初期音圧分布の再構成を行うことを可能にする。【解決手段】光音響波の検出信号を取得するステップと、前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を備える光音響情報処理方法において、前記音圧分布データとしてｋ空間上のデータを用い、前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う。【選択図】図４
2018033886.tif 000002
【特許請求の範囲】
【請求項１】
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を備え、
前記音圧分布データとしてｋ空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理方法。
【請求項２】
前記球殻積分にあたり、球殻を複数領域に分割するとともに、各領域で積分変数をそれぞれ選択する
請求項１に記載の光音響情報処理方法。
【請求項３】
前記初期音圧分布の再構成の結果に基づいて光音響画像を生成するとともに、生成した前記光音響画像についての画像データをデータ圧縮して記憶保持する
請求項１または２に記載の光音響情報処理方法。
【請求項４】
コンピュータに、
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を実行させるとともに、
前記音圧分布データとしてｋ空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理プログラム。
【請求項５】
光音響波の検出信号を取得する信号取得部と、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成する再構成部と、
前記再構成部での処理結果を出力する結果出力部と、を備え、
前記音圧分布データとしてｋ空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、光音響情報処理方法、光音響情報処理プログラムおよび光音響情報処理装置に関する。
【背景技術】
【０００２】
近年、例えば医療分野において、光音響効果を用いて被検体内部を画像化する光超音波イメージング法（photoacoustic imaging）を利用することが提案されている。光音響効果とは、被検体にレーザ光を照射すると、被検体内部における血液等の光吸収体に吸収され、局所的な温度上昇が熱弾性変形を引き起こし、その熱弾性変形による体積膨張で音響波（超音波）が発生する現象である。光超音波イメージング法では、発生した音響波を受信して、その音源（すなわち、被検体内部における光吸収体）の初期音圧分布を再構成することで、被検体内部を画像化する（例えば、非特許文献１参照）。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】椎名毅、「光音響を用いた新しい超音波診断装置」、メディカル＆イメージングＮｏ．３、オプトロニクス社、２０１５年、ｐ．５３−５９
【発明の概要】
【発明が解決しようとする課題】
【０００４】
光超音波イメージング法の医療適用にあたっては、画像化のために必要となる初期音圧分布の再構成を、画質劣化を抑制しつつ、迅速かつ効率的に行うことが求められる。
【０００５】
本発明は、画質劣化を抑制しつつ、迅速かつ効率的に、初期音圧分布の再構成を行うことができる光音響情報処理方法、光音響情報処理プログラムおよび光音響情報処理装置を提供することを目的とする。
【課題を解決するための手段】
【０００６】
本発明の一態様によれば、
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を備え、
前記音圧分布データとしてｋ空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理方法が提供される。
【０００７】
本発明の他の一態様によれば、
コンピュータに、
光音響波の検出信号を取得するステップと、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成するステップと、を実行させるとともに、
前記音圧分布データとしてｋ空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理プログラムが提供される。
【０００８】
本発明のさらに他の一態様によれば、
光音響波の検出信号を取得する信号取得部と、
前記検出信号から音圧分布データを求め、前記音圧分布データに基づいて初期音圧分布を再構成する再構成部と、
前記再構成部での処理結果を出力する結果出力部と、を備え、
前記音圧分布データとしてｋ空間上のデータを用い、
前記初期音圧分布の再構成を、信号のスパース性を利用した圧縮センシングを適用して行うとともに、前記圧縮センシングにおける復元計算を離散条件の球殻積分により行う
光音響情報処理装置が提供される。
【発明の効果】
【０００９】
本発明によれば、画質劣化を抑制しつつ、迅速かつ効率的に、初期音圧分布の再構成を行うことができる。
【図面の簡単な説明】
【００１０】
【図１】光音響イメージングの基本的な手順を例示する説明図である。
【図２】本発明の一実施形態に係る光音響情報処理装置を含む光音響イメージング装置の概略構成の例を模式的に示す説明図である。
【図３】順問題による繰り返し推定法を用いて再構成を行う場合の一般的な手法の概要を例示する説明図である。
【図４】本発明の一実施形態に係る光音響情報処理方法の基本的な手順を例示する説明図であり、圧縮センシングを適用して順問題を解いて再構成を行う手法の概要を示す図である。
【図５】圧縮センシングの意義を具体的に説明するための説明図である。
【図６】本発明の一実施形態に係る光音響情報処理方法における球殻積分の具体的な一態様の概要を具体的に説明するための説明図である。
【発明を実施するための形態】
【００１１】
＜本発明の一実施形態＞
以下に、本発明の実施の形態について、図面を参照しながら説明する。
【００１２】
（１）光超音波イメージング法の概要
はじめに、光音響効果を用いた光超音波イメージング法（以下、単に「光音響イメージング」ともいう。）の概要を説明する。
【００１３】
（光音響イメージングの手順）
図１は、光音響イメージングの基本的な手順を例示する説明図である。
図例のように、光音響イメージングにあたっては、以下に説明する各ステップ（以下ステップを「Ｓ」と略す。）を順に経る。
【００１４】
（Ｓ１０１：レーザ光照射）
先ず、被検体１に対して光源部１１からレーザ光を照射する。
レーザ光の波長は、検出すべき対象の光吸収特性により決まる。生体計測の場合であれば、水の吸収が少ない可視光から近赤外領域が適している。例えば、血管（血液）の描出や酸素飽和度を求める場合は、酸素化ヘモグロビンと脱酸素化ヘモグロビンのスペクトルが交差する８００ｎｍ付近の複数の波長のレーザ光を用いる。
また、光の吸収により効果的に超音波を発生させるためには、熱の拡散の時間や吸収体内を超音波が伝搬する時間に比べて、レーザ光の照射時間であるパルス幅が十分短い必要があることから、照射すべきレーザ光としてナノ秒幅のレーザ光パルスを用いる。
【００１５】
（Ｓ１０２：音響波放射）
レーザ光を照射すると、被検体１の内部では、血液等の吸収体２が光を吸収して、局所的な温度上昇が熱弾性変形を引き起こし、その熱弾性変形による体積膨張で吸収体２が音響波（初期音圧ｐ０（ｒ）)を放射する。このとき、光の吸収により発生する音響波（超音波）の大きさは、以下の（１）式で表される。
【００１６】
【数１】
2018033886.tif 000003
【００１７】
ここで、βは熱膨張係数、ｃは音速、Ｃｐは定圧比熱、μａは光の吸収係数、Ψは光量である。係数Ｇ＝βｃ２／Ｃｐは、光エネルギーから音圧への変換係数でグリューナイゼン（Gruneisen）係数と呼ばれる。Ｇは軟組織間であまり差はないため、初期音圧ｐ０は、組織の吸収係数μａと光量Ψの積に比例することになる。
【００１８】
また、吸収体２で発生した位置ｒ、時刻ｔにおける音響波ｐ（ｒ，ｔ）は、以下の（２）式で表わされる光音響波動方程式に従って、被検体１の内部を伝搬する。
【００１９】
【数２】
2018033886.tif 000004
【００２０】
（Ｓ１０３：音響波計測）
その後は、吸収体２で発生し被検体１の体表まで伝搬した音響波ｐ（ｒ，ｔ）を、体表近傍に配されて複数の圧電変換素子（トランスデューサ）を有するセンサ部１２を用いて計測する。これにより、センサ部１２では、吸収体２からの音響波ｐ（ｒ，ｔ）の検出信号（受信超音波信号）が電気信号として得られることになる。
【００２１】
（Ｓ１０４：再構成）
センサ部１２で受信超音波信号を得ると、その受信超音波信号によって特定される音圧分布ｐ（ｒ，ｔ）から、初期音圧分布ｐ０（ｒ）を再構成する。
本実施形態においては、ここで行う再構成の手法に大きな特徴があるが、この点については詳細を後述する。
【００２２】
（Ｓ１０５：光量補正）
初期音圧分布ｐ０（ｒ）を再構成したら、光源部１１の仕様から光量Ψの分布が既知であるので、その初期音圧分布ｐ０（ｒ）に対して上記の（１）式に基づき光量Ψを補正して、被検体１の内部における吸収係数の分布像を生成する。
【００２３】
このようにして得られる吸収係数分布像は、必要に応じて、いわゆる光音響画像として出力されることになる。光音響画像は、二次元画像または三次元画像の別を問わないが、その一例として最大値投影画像を含むものが知られている。最大値投影画像とは、被検体の画像をＭＩＰ（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の順問題を解いてベクトルＰｋを求め、これによりセンサ部１２で検出されたであろう音響波の推定値ｐｋ（ｒ，ｔ）を推定する（Ｓ３０３〜Ｓ３０４）。
【００５０】
このことは、伝達関数Φと基底Ψとの積を新たな伝達関数Θと考えてΦΨ＝Θとすると、以下の（３）式が成り立つことを意味する。
【００５１】
【数３】
2018033886.tif 000005
【００５２】
つまり、上述した各ステップ（Ｓ３０２〜Ｓ３０４）は、圧縮センシングを適用することで、Ｐk＝ΦＰ０kの順問題を解く代わりに、上記の（３）式に基づき、Ｐｋ＝ΘＳｋの関係を満たす信号ベクトルＳｋで最もスパースなものを解とし、その解からベクトルＰｋを求めることに相当する（Ｓ３０５）。
【００５３】
その後は、音響波の推定値ｐｋ（ｒ，ｔ）を推定したら、これを音響波ｐ（ｒ，ｔ）の計測値と比較して（Ｓ３０６）、それぞれの間の誤差が予め設定されている基準値を超えているか否かを判断する（Ｓ３０７）。その結果、誤差が基準値を超えていなければ、音響波の推定値ｐｋ（ｒ，ｔ）が音響波ｐ（ｒ，ｔ）の計測値であるものと見做す。そして、その推定値ｐｋ（ｒ，ｔ）を導き出したときの初期音圧分布の設定値ｐ０k-1（ｒ）のベクトルＰ０kを出力して処理を終了する。
【００５４】
また、誤差が基準値を超えている場合には、ｋ＝ｋ＋１にインクリメントした上で（Ｓ３０８）、誤差が小さくなるようにスパースな信号ｓk-1（ｒ）を修正する。具体的には、例えば、初期値ｓ０（ｒ）から推定した音響波の推定値ｐｋ（ｒ，ｔ）についての誤差が基準値を超えている場合であれば、その誤差が小さくなるように、初期値ｓ０（ｒ）を基にしつつ次に用いるべきスパースな信号ｓk-1（ｒ）を計算して求める。つまり、修正の都度、誤差が小さくなるように、スパースな信号ｓk-1（ｒ）を更新するのである。そして、その更新後のスパースな信号ｓk-1（ｒ）について、再び上述した一連の処理を繰り返す（Ｓ３０３〜Ｓ３０８）。これを処理が終了するまで、すなわち誤差が基準値を超えなくなるまで、繰り返し行う。
【００５５】
以上のような手順を経ることで、圧縮センシングを適用しつつ、初期音圧分布についてのベクトルＰ０kを求めることが可能となる。しかも、圧縮センシングを適用することで、信号のスパース性を利用することになるので、例えば、Ｐk＝ΦＰ０kの順問題を解く際の伝達関数Φがｍ×ｎ行列である場合にｍ＜ｎであっても、解を求めることが可能となる。つまり、再構成領域から得られる検出信号が少ない場合であっても、真の信号を推定することが可能となる。
【００５６】
図５は、圧縮センシングの意義を具体的に説明するための説明図である。
例えば、図５（ａ）に示すように、Ｐ＝ΦＰ０の関係に基づいて、ベクトルＰからベクトルＰ０を求める場合には、伝達関数Φがｍ×ｎ行列であるときにｍ＜ｎであると、他の制約条件がないと不良設定問題となってしまい、解が求まらない。
これに対して、圧縮センシング（制約条件としてＬ１最適化）を適用した場合には、例えば、図５（ｂ）に示すように、Ｐ０＝ΨＳの関係を満たすスパースな信号ベクトルＳを用いることで、上記の（３）式の関係が成り立つ。ここで、スパースな信号ベクトルＳは、Ｋ個の要素以外は「０」であると考えられる。したがって、Ｐ＝ΘＳの関係に基づいてベクトルＰから信号ベクトルＳを求める場合には、ｍ＞Ｋであれば良設定（well-posed）問題となり、解を求めることができる。つまり、Ｐ０＝ΨＳの関係を利用することで、ベクトルＰ０を求めることができる。
【００５７】
このような圧縮センシングを適用することで、再構成領域から得られる検出信号が少ない場合であっても、Ｐk＝ΦＰ０kの順問題を解くことが可能となる。したがって、逆問題を直接解きベクトルＰ０を求めて再構成を行う場合に比べると、再構成アーチファクトの発生を抑制することができ、その結果として再構成した初期音圧分布に基づいて生成する光音響画像の画質劣化を抑制して画質向上が図れるようになる。なお、生成された光音響画像は、例えば、結果出力部１３ｃが画像表示部１４へ出力して、その画像表示部１４にて表示出力されたり、あるいは画像記憶部１３ｄにより記憶保持されたりすることになる。
【００５８】
（処理の高速化の概要）
ところで、圧縮センシングでは、上記の（３）式に基づく信号ベクトルＳからベクトルＰへの生成（順問題）の過程で、伝達関数（行列）Θについての膨大な計算量を必要とする。通常、圧縮センシングを適用する場合には、信号ベクトルＳからベクトルＰへの生成を、計測信号の空間領域（実空間）で処理することが一般的である。そのため、下記の（４）式に示すように、被検体１内の全ての音源で生成され、センサ部１２における全てのトランスデューサ１２ａに到達した波形を、インパルス応答として重ね合わせる方式で計算することになり、その計算に膨大な時間がかかってしまう。
【００５９】
【数４】
2018033886.tif 000006
【００６０】
ここで、ｉは虚数単位、ｃは音速、ｋはｋ＝ω／ｃで表されるもので、ωは角周波数、ｇ（ｋ）は、トランスデューサ１２ａのインパルス応答、Ａ（ｒ´）は音源位置ｒ´の吸収係数、｜ｒ-ｒ´｜はトランスデューサ位置ｒと音源位置ｒ´との距離である。
【００６１】
上記の（４）式に基づいて計算を行う場合には、例えば、周波数帯域の分割数をＮ、トランスデューサ１２ａの数をＭ、被検体１内の音源の座標（ｘ，ｙ，ｚ）点数をＩ×Ｊ×Ｋとすると、その計算回数が積＝ＭＮＩＪＫ回、和＝ＭＮＩＪＫ回となり、膨大なものとなってしまう。
【００６２】
そこで、本実施形態においては、計算量を削減して処理の高速化を図り、その結果として圧縮センシングを適用した再構成を迅速かつ効率的に行い得るようにするために、以下に述べる（ｉ）および（ii）の特徴点を備えている。
【００６３】
すなわち、本実施形態では、（ｉ）ｋ空間（ｋ-space）でのデータを用いて計算を行うことで、計算量を削減して処理の高速化を図っている。
ｋ空間とは、空間周波数を表す空間のことであり、実空間とは互いにフーリエ変換の関係にある空間のことである。
【００６４】
さらに、本実施形態では、（ii）ｋ空間での計算を行うことを前提としたうえで、圧縮センシングにおける復元計算を離散条件の球殻積分により行う。
圧縮センシングにおける復元計算とは、トランスデューサ１２ａによる受信波のスペクトル（時間ｔについてのフーリエ変換）を求めて復元するために必要となる積分計算のことである。また、離散条件の球殻積分とは、ｋ空間に存在する球殻の表面上の離散データのみを用いて積分計算をすることである。
【００６５】
具体的には、上述した（ｉ）および（ii）の特徴点を備えることで、本実施形態では、以下の（５）式および（６）式に基づく計算を行う。
すなわち、本実施形態においては、空間を三次元フーリエ変換したｋ空間で表現することで、計測ベクトルＰのｋ空間表現である（５）式と、初期音圧ｐ０（ｒ）の空間座標ｒに関する三次元フーリエ変換（ｋ空間表現）である（６）式が得られる。
【００６６】
【数５】
2018033886.tif 000007
【００６７】
【数６】
2018033886.tif 000008
【００６８】
なお、上記の（５）式において、Ｓは、ｋ空間での半径｜ω｜／ｃの球殻上の積分を意味している。
【００６９】
このように、上記の（５）式によれば、ｋ空間での球殻積分を行うので、二重積分（球殻上積分）となっており、上記の（４）式のような実空間での三重積分を必要とする場合に比べると、計算量を削減して処理の高速化を図ることが可能となる。
【００７０】
（高速化した場合の具体的な手順）
ここで、上述した（ｉ）および（ii）の特徴点を備えつつ初期音圧分布の再構成を行う場合の具体的な手順、すなわち本実施形態における再構成の手法の具体的な手順について説明する。
【００７１】
本実施形態における具体的な手順は、図４を用いて説明した基本的な手順に当て嵌めると、以下に述べる（Ｉ）〜（VIII）の各手順を経る。
【００７２】
（Ｉ）センサ部１２で受信された音響波ｐ（ｒ，ｔ）の計測値を信号取得部１３ａが取得する（Ｓ１０１〜Ｓ１０３）。
【００７３】
（II）信号取得部１３ａが取得したｐ（ｒ，ｔ）の計測値について、再構成部１３ｂが時間ｔについてフーリエ変換し、ｐ（ｒ，ω）を求める。これにより、被検体１における音圧分布に関する計測値データとして、ｋ空間上のデータを用いることになる。
【００７４】
（III）再構成部１３ｂが初期値ｓ０（ｒ）を設定する（Ｓ３０１〜Ｓ３０２）。
【００７５】
（IV）スパースな信号ｓｋ（ｒ）の初期値ｓ０（ｒ）について、Ｐ０k＝ΨＳｋの関係に基づき、再構成部１３ｂが基底Ψで変換して、ベクトルＰ０k（ｒ）を得る（Ｓ３０２〜Ｓ３０３）。このとき、基底Ψは、予め定められたものであれば、適切なものを適宜選択して用いることが考えられる。具体的には、例えば、カーブレット基底、フーリエ基底やウェーブレット基底等のいずれかを用いることが考えられる。
【００７６】
（Ｖ）ベクトルＰ０k（ｒ）について、再構成部１３ｂが上記の（６）式に基づき、空間を三次元フーリエ変換し、ｋ空間表現のベクトルＰ０k（ｋ）を求める。これにより、被検体１における音圧分布に関するベクトルデータとして、ｋ空間上のデータを用いることになる。
【００７７】
（VI）ｋ空間表現のベクトルＰ０k（ｋ）から、ｋ空間での演算式である上記の（５）式を用いて、再構成部１３ｂが受信波のスペクトル（時間ｔについてのフーリエ変換）を求める。このとき、再構成部１３ｂは、上記の（５）式を用いることで、ｋ空間での離散条件の球殻積分を行うことになる。その結果、再構成部１３ｂは、Ｐｋ＝ΦＰ０kの順問題を解いてベクトルＰｋを求め、これによりセンサ部１２で検出されたであろう音響波の推定値ｐｋ（ｒ，ω）を推定することになる（Ｓ３０３〜Ｓ３０５）。
【００７８】
（VII）その後、再構成部１３ｂは、上記（II）で求めたｐ（ｒ，ω）と上記（VI）で求めたｐｋ（ｒ，ω）とを比較し（Ｓ３０６）、これらの間の誤差が基準値を超えている場合にはｓｋ（ｒ）を修正してｓｋ＋１（ｒ）を求める（Ｓ３０７〜Ｓ３０８）。つまり、ｋ→ｋ＋１として、上記（IV）以降の処理を繰り返し行うのである。
このとき、再構成部１３ｂは、例えば、以下の（７）式によって特定される公知のＬ１最適化アルゴリズムを用いることで、上述した繰り返し処理を効率的に収束させるようにしてもよい。
【００７９】
【数７】
2018033886.tif 000009
【００８０】
（VIII）ｐ（ｒ，ω）とｐｋ（ｒ，ω）との誤差が基準値以下の場合、再構成部１３ｂは、そのｐｋ（ｒ，ω）を導き出したときのスパースな信号ｓｋ（ｒ）を結果出力部１３ｃに出力させて処理を終了する。
【００８１】
（高速化した計算量の具体例）
本実施形態によれば、以上に説明したような処理を行うので、計算量を削減して処理の高速化を図ることができる。
【００８２】
具体的には、例えば、周波数帯域の分割数をＮ、トランスデューサ１２ａの数をＭ、被検体１内の音源の座標点数をＩ×Ｊ×Ｋとすると、上記の（６）式による計算回数がＩＪＫ・log（ＩＪＫ）回、上記の（５）式による計算回数が球殻積分で最大でもＩＪＫ回となり、本実施形態の全体ではＩＪＫ（Ｍ＋log（ＩＪＫ））回となる。
これに対して、本実施形態のようなｋ空間ではなく実空間において上記の（４）式に基づいて計算を行う場合には、既に説明したように、計算回数がＭＮＩＪＫ回となってしまう。
したがって、例えば、Ｍ，Ｎがそれぞれ１００回、Ｉ，Ｊ，Ｋがそれぞれ１０００回であると、実空間における上記（４）式の計算ではＭＮＩＪＫ＝１００２×１０００３＝１０13回の計算回数が必要であるのに対して、本実施形態ではＩＪＫ（Ｍ＋log（ＩＪＫ））＝１０００３（１００＋log（１０００３）＝１０９×１３０≒１０11回となる。つまり、本実施形態によれば、凡そ二桁程度の計算回数の削減による高速化が実現可能となる。
【００８３】
（球殻積分の具体的な一態様）
上述したように、本実施形態では、ｋ空間での計算を行うので、圧縮センシングにおける復元計算を球殻積分によって行うことが可能となる。そして、その球殻積分にあたり、ｋ空間に存在する球殻の表面上の離散データのみを用いて積分計算をすればよい。
この球殻積分については、以下に述べる態様で行えば、積分計算によって得られるデータの精度向上を図るうえで非常に好ましいものとなる。
以下に、本実施形態における球殻積分の好適な一態様について具体的に説明する。なお、以下に説明する例は、球殻積分の好適な一態様に過ぎず、本発明がこれに限定されるものではない。
【００８４】
図６は、本実施形態に係る光音響情報処理方法における球殻積分の具体的な一態様の概要を具体的に説明するための説明図である。
球殻積分にあたっては、例えば、ｋｘ，ｋｙ，ｋｚの各軸を有する座標空間（ただし、ｋｚ軸は不図示）において、ｋ＝ωｃで表される球殻を積分計算によって求める。
【００８５】
その場合において、積分計算に際しては、先ず、計算対象となる領域を｛±ｋｘ＞｜ｋｙ｜，｜ｋｚ｜｝、｛±ｋｙ＞｜ｋｘ｜，｜ｋｚ｜｝、｛±ｋｚ＞｜ｋｘ｜，｜ｋｙ｜｝の６領域に分割する。これにより、信号ベクトルＳについては、以下の（８）式および（９）式が成り立つことになる。
【００８６】
【数８】
2018033886.tif 000010
【００８７】
【数９】
2018033886.tif 000011
【００８８】
ここで、例えば、Ｓｘ＋＝｛ｋｘ＞｜ｋｙ｜，｜ｋｚ｜｝の場合に着目すると、この領域については、ｋｙ，ｋｚを積分変数として積分し、全てのグリッド上のｋy0，ｋz0について以下の（１０）式のｋｘを求める。
【００８９】
【数１０】
2018033886.tif 000012
【００９０】
そして、ｋｘの最近傍のグリッド２点ｋx0，ｋx1を線形補間して、ｐ（ｋｘ，ｋy0，ｋz0）を求める。
【００９１】
このような処理を分割した各領域の全てについて行うことで、半径｜ω｜／ｃの球殻について、その表面上の全領域についての積分結果を得ることができるようになる。
【００９２】
つまり、ここで例に挙げて説明した球殻積分の一態様では、その球殻積分にあたり、球殻を予め設定された数の複数領域（例えば６領域）に分割するとともに、予め設定された選択基準に従いつつ各領域で積分変数をそれぞれ選択する。そして、それぞれ選択した積分変数を各領域別に適用しつつ、球殻の全域についての積分計算を行うのである。
【００９３】
このような球殻積分の一態様によれば、図６に示す例において二次元に簡略化して考えると、例えば、領域Ａについてはｋｙ軸に沿った積分変数を選択し、領域Ｂについてはｋｘ軸に沿った積分変数を選択する、といったことが実現可能となる。したがって、各領域について、変動の大きな軸方向に沿って積分変数を選択し、その選択した積分変数を適用して積分計算を行うことができる。つまり、全領域に対して積分変数を一律に適用する場合に比べると、球殻の形状の特質を適切に考慮した積分計算を行うことができるので、球殻積分を行う場合であっても、その積分計算によって得られるデータの精度向上が図れる。その結果として、再構成した初期音圧分布に基づいて生成する光音響画像の画質向上が期待できる。
【００９４】
（光音響画像の記憶保持）
以上に説明した手順を経ることで、本実施形態では、被検体１における吸収体２について、初期音圧分布ｐ０（ｒ）の再構成を行う。そして、初期音圧分布ｐ０（ｒ）を再構成再構成したら、必要に応じて、光量補正を行って吸収係数分布像（すなわち光音響画像）を生成する。このようにして生成された光音響画像は、結果出力部１３ｃから画像表示部１４へ出力されたり、あるいは画像記憶部１３ｄで記憶保持されたりする。
【００９５】
画像記憶部１３ｄでの記憶保持にあたっては、光音響画像の画像データ量が膨大であることが多いため、画像データのデータ圧縮を行うことが考えられる。データ圧縮の手法については、特に限定されることはなく、公知のデータ圧縮手法を用いればよい。
【００９６】
その場合において、本実施形態で得られる光音響画像は、上述したように圧縮センシングを適用した再構成を経て生成されているため、例えばバックプロジェクション法を適用した場合に比べると、ぼけやストリークアーチファクト等の再構成アーチファクトが発生してしまうことがない。そのため、光音響画像についてデータ圧縮を行うと、再構成アーチファクトが少ない分だけ、例えばバックプロジェクション法を適用した場合に比べて、画像圧縮率が高くなる。
【００９７】
つまり、本実施形態によれば、光音響画像についてデータ圧縮を行う際の画像圧縮率を高くすることができ、画像記憶部１３ｄのために必要とするデータ記憶容量が増大してしまうのを抑制することが可能となる。
【００９８】
（３）本実施形態の効果
本実施形態によれば、以下に示す一つまたは複数の効果を奏する。
【００９９】
（ａ）本実施形態では、初期音圧分布の再構成にあたり、逆問題を直接解く代わりに、順問題を解いて音響波の推定値ｐｋ（ｒ，ｔ）を推定し、これを基準値との誤差が小さくなるまで繰り返す。したがって、逆問題を直接解く場合とは異なり、ぼけやストリークアーチファクト等の再構成アーチファクトが発生してしまうことがなく、再構成結果に基づいて光音響画像を生成する際の画質劣化を抑制することができる。
【０１００】
（ｂ）また、本実施形態では、順問題を解いて音響波の推定値ｐｋ（ｒ，ｔ）を推定するのにあたり、信号のスパース性を利用した圧縮センシングを適用する。したがって、信号のスパース性を利用することになるので、不良設定問題とはならずに解を求めることが可能となる。つまり、再構成領域から得られる検出信号が少ない場合であっても、真の信号を推定することが可能となる。このことは、光音響イメージングに適用して非常に好適であることを意味する。光音響イメージングは、再構成領域を囲む全周について検出信号を稠密に得ることが容易ではなく、また十分に広い帯域の検出信号が得されるとは限らないからである。
【０１０１】
（ｃ）また、本実施形態では、初期音圧分布の再構成にあたり、音圧分布データ（具体的には、音圧分布に関する計測値データや音圧分布に関するベクトルデータ等）としてｋ空間上のデータを用い、ｋ空間での計算を行うことを前提としたうえで、圧縮センシングにおける復元計算を離散条件の球殻積分により行う。つまり、本実施形態では、ｋ空間での球殻積分を行う。したがって、本実施形態によれば、例えば実空間での積分計算を必要とする場合に比べると、計算量を削減して処理の高速化を図ることが可能となる。具体的には、凡そ二桁程度の計算回数の削減による高速化が実現可能となる。このように、初期音圧分布の再構成に必要となる計算処理の高速化を図ることで、本実施形態では、初期音圧分布の再構成を迅速かつ効率的に行うことができる。
【０１０２】
（ｄ）また、本実施形態では、球殻積分にあたり、球殻を複数領域に分割するとともに、各領域で積分変数をそれぞれ選択する。したがって、全領域に対して積分変数を一律に適用する場合に比べると、球殻の形状の特質を適切に考慮した積分計算を行うことができるので、球殻積分を行う場合であっても、その積分計算によって得られるデータの精度向上が図れ、その結果として、再構成結果に基づいて光音響画像を生成する際の画質向上が期待できる。
【０１０３】
（ｅ）また、本実施形態では、再構成結果に基づいて生成した光音響画像の記憶保持に際し、その光音響画像の再構成アーチファクトが少ないことから、その光音響画像についてデータ圧縮を行う際の画像圧縮率を高くすることができる。
【０１０４】
＜他の実施形態＞
以上に、本発明の一実施形態を具体的に説明したが、本発明は上述の実施形態に限定されるものではなく、その要旨を逸脱しない範囲で種々変更可能である。
【０１０５】
例えば、本実施形態では、主として、被検体１の内部における血液等の吸収体２について、その吸収係数分布像を得るための光音響イメージングを例に挙げたが、本発明がこれに限定されることはない。すなわち、本発明は、血液等の吸収体以外の光音響イメージングにも適用することが可能である。
【０１０６】
また、例えば、本実施形態では、光音響イメージング装置１０を構成する各部１１〜１４、および、制御部１３を構成する各部１３ａ〜１３ｂについて具体的に説明したが、上述した実施形態で例に挙げた具体例に限定されることはなく、同等のものに代替して構成しても構わない。具体的には、その一例として、画像表示部１４に代わって、画像表示機能を備えた外部装置に有線または無線で接続するインタフェース部を設ける、といったことが考えられる。
【０１０７】
また、例えば、本実施形態では、画像の再構成にあたり、上記の（１）式〜（１０）式を基にする場合を例に挙げたが、各式はその要旨を変更しない範囲で適宜修正しても構わない。
【符号の説明】
【０１０８】
１…被検体、２…吸収体、１０…光音響イメージング装置、１１…光源部、１２…センサ部、１２ａ…圧電変換素子（トランスデューサ）、１３…制御部（光音響情報処理装置）、１３ａ…信号取得部、１３ｂ…再構成部、１３ｃ…結果出力部、１３ｄ…画像記憶部、１４…画像表示部
【図１】
2018033886.tif 000013
【図２】
2018033886.tif 000014
【図３】
2018033886.tif 000015
【図４】
2018033886.tif 000016
【図５】
2018033886.tif 000017
【図６】
2018033886.tif 000018