WO2013076984A1

WO2013076984A1 - 光音響信号処理装置及び方法

Info

Publication number: WO2013076984A1
Application number: PCT/JP2012/007492
Authority: WO
Inventors: 覚入澤; 和弘広田
Original assignee: 富士フイルム株式会社
Priority date: 2011-11-22
Filing date: 2012-11-21
Publication date: 2013-05-30
Also published as: JP2013128758A; US11213209B2; JP5779567B2; US20140257079A1

Abstract

【課題】実用的な装置において、検出信号から吸収分布を求める。【解決手段】被検体に対して光が照射されることで、被検体内で光音響信号が生じる。プローブ（１１）で検出された光音響信号は、受信回路を介して、ＡＤ変換手段（２２）でサンプリングされる。光微分波形逆畳込み手段（２５）は、被検体に照射された光の時間波形の微分波形を光微分波形として、光音響信号から光微分波形をデコンボリューションする。このデコンボリューションにより、吸収分布が得られる。

Description

光音響信号処理装置及び方法

　本発明は、光音響信号処理装置及び方法に関し、更に詳しくは、被検体に照射された光により被検体内で生じた光音響信号を信号処理する光音響信号処理装置及び方法に関する。

　生体内部の状態を非侵襲で検査できる画像検査法の一種として、超音波検査法が知られている。超音波検査では、超音波の送信及び受信が可能な超音波探触子を用いる。超音波探触子から被検体（生体）に超音波を送信させると、その超音波は生体内部を進んでいき、組織界面で反射する。超音波探触子でその反射音波を受信し、反射超音波が超音波探触子に戻ってくるまでの時間に基づいて距離を計算することで、内部の様子を画像化することができる。

　また、光音響効果を利用して生体の内部を画像化する光音響イメージングが知られている。一般に光音響イメージングでは、パルスレーザ光を生体内に照射する。生体内部では、例えば生体組織がパルスレーザ光のエネルギーを吸収し、そのエネルギーによる断熱膨張により超音波（光音響信号）が発生する。この光音響信号を超音波プローブなどで検出し、検出信号に基づいて光音響画像を構成することで、光音響信号に基づく生体内の可視化が可能である。

　ここで、通常の再構成方法（例えばフーリエドメイン法（ＦＴＡ法）や遅延加算法（Delay&Sum法））は、実質的に圧力分布を画像化しており、生成された光音響画像は光吸収体の分布画像にはなっていない。圧力分布を画像化した光音響画像では、１本の血管が二重に表示されることがある。この現象は、以下のように説明できる。すなわち、微小要素（微小吸収体）が発生する圧力波形をミクロ波形とすると、マクロ構造物としての血管から発生する圧力波形は、ミクロ波形の重ね合わせとして捉えることができる。微小吸収体で発生したミクロ波形を重ね合わせていくと、血管内部で発生したミクロの圧力波形は近接した前後のミクロ波形と重ね合わされてキャンセルされる。血管から出たミクロの圧力波形を全て重ね合わせると、最終的に血管表面で発生した圧力波形の一部（前端）と裏面で発生した圧力波形の一部（後端）とがキャンセルされずに残り、これらが観測されることになる。そのため、血管が二重に表示される。１本の血管が二重に表示されると、画像判定上、血管の位置の確認が困難になるという問題が生じる。また、複数の異なる波長で取得した画像同士を計算する際に、位置がずれやすく適切な結果が得られないなどの問題も生じる。

　圧力分布画像に代えて、吸収分布画像を生成する技術がこれまでにいくつか知られている。例えば特許文献１には、試料の熱的インパルス応答から光音響画像の分解能劣化を修復する逆フィルタを求め、得られた光音響画像に逆フィルタを作用させることで、理想的な光音響画像、すなわち試料表面の点光源によって励起され、検出されるその点（無限小）における熱的インピーダンス情報（＝無限小なる点熱源の集合）を得ることが記載されている。

　特許文献１では、より詳細には、まず、試料の熱インパルス応答ｈ（ｘ，ｙ）を計算し、次いで光音響画像ｐ（ｘ，ｙ）を構成する。熱的インパルス応答は、無限小なる一点の温度変化が試料表面の微小変位に変換されるまでの伝達関数と定義されている。その後、熱インパルス応答ｈ（ｘ，ｙ）と光音響画像ｐ（ｘ，ｙ）とをそれぞれフーリエ変換し、フーリエ変換像Ｈ（μ，υ）、Ｐ（μ，υ）を得る。１／Ｈ（μ，υ）を逆フィルタとして用い、Ｑ（μ，υ）＝Ｐ（μ，υ）・（１／Ｈ（μ，υ））により、Ｑ（μ，υ）を計算する。そのようにして計算されたＱ（μ，υ）をフーリエ逆変換することで、理想的な光音響画像ｑ（ｘ，ｙ）が得られる。

　上記特許文献の他にも、非特許文献１には、論理的には、有限の時間幅を持つ光パルスη（ｔ）をフーリエ変換したものをη（ｋ）としたとき、その微分をｉη（ｋ）として考慮することが記載されている。実験的には、超音波検出素子であるＰＺＴ（チタン酸ジルコン酸鉛）での検出帯域以内に励起光パルス波形が入るようにパルス幅を長くしたマイクロ波を被検体に照射し、通常のＰＺＴプローブで光音響信号を検出し、吸収分布を再構成している。

　また、非特許文献２には、被写体中の微小要素からの圧力波形として光パルス微分関数と装置インパルス応答関数とを合わせたミクロ波形と、吸収分布とを、観測圧力波形に関連付けることが記載されている。吸収像再構成は、光微分とシステム応答とを不可分な状態で含むｐ_ｄ０を測定して、各素子の圧力波形からｐ_ｄ０をデコンボリューションしたのちに、フィルタ補正逆投影法（Filtered Backprojection法）を用いて行う。実験的には、パルス幅の短いパルスレーザ光で励起し、超音波の検出帯域を通常の超音波診断装置よりも広げて、ハイドロホンとオシロスコープとで光音響信号を検出し、吸収分布を再構成する。

特開平３－１５６３６２号公報

Yuan Xu, et al., IEEE Transactions on Medical Imaging, Volume 21 (2002), pp.823-828. Yi Wang, et al., Physics in Medicine and Biology, Volume 49 (2004), pp.3117-3124

　非特許文献１では、光パルスをフーリエ変換したη（ｋ）を位置依存の関数η（ｒ，ｋ）として扱っていない。このため、ｔ＝０に再構成できる（検出した波形を基に時間を遡って光が入射した瞬間（ｔ＝０）に発生した圧力分布を推定する計算を行うことができる）場合には正確な吸収分布を得ることができるものの、ｔ≠０に再構成した場合（光が入射した瞬間（ｔ＝０）の圧力分布を再構成できずに、時刻ｔ＝０からしばらく時間が経った後の圧力分布を推定することになる場合）には、光パルス幅成分が除去できず、圧力分布となる。

　また、非特許文献１及び２では、励起レーザと超音波検出装置とのどちらかを実用的な範囲から外してレーザ発光時間と超音波検出時とを合わせて再構成している。このため、非特許文献１及び２では問題は明確には現れないものの、実用的な装置構成を考えた場合、非特許文献１及び２の手法では、ｔ＝０とする再構成が困難である。すなわち、例えば実用的な装置構成として、
　・サンプリング周波数１００ＭＨｚ以下で、ＰＺＴなどを利用した狭帯域プローブを用いた超音波検出装置
　・強い光音響信号が出る１－１００ｎｓオーダーの光パルス幅を持つ励起レーザ
を用いた場合、レーザ発光が超音波検出時間と比較して短時間の現象のため、ｔ＝０に相当する状態（吸収分布と圧力分布とが比例する時間帯）に正確に再構成できない。

　ここで、「時刻ｔ＝０の圧力分布」は「吸収分布」を表すので、時刻ｔ＝０の圧力分布が求められれば吸収分布を得ることができる。しかし、一般的な超音波検出装置のサンプリング間隔は２５ｎｓ程度であり、光が当たった瞬間の時刻ｔ＝０のつもりで計算しても、実際はｔ＝±１２．５ｎｓ程度の時間幅でずれが生じる。例えば光パルス幅が１００ｎｓと長い場合には、上記のずれ（±１２．５ｎｓ）は誤差と考えればよいものの、光パルス幅が１０ｎｓであれば、上記のずれは誤差とは呼べなくなり、光が当たった瞬間の圧力分布というよりも圧力波の伝播過程の圧力分布に移行することとなる。その「圧力波の伝播過程の圧力分布」は、「吸収分布」に一致しない。

　更に、実験的な生体等のサンプルにおいては、ｔ＝０の圧力分布を定義することが困難である。生体内の音速を例えば１５３０ｍ／ｓと仮定し、検出時刻とレーザ照射時刻との差を伝播時間とすると、伝播時間から伝播距離が求められる。生体内の音速が１５３０ｍ／ｓで一定であれば、伝播時間から求めた伝播距離は実際の伝播距離と一致する。しかし、実際には、音速は生体内で一様ではなく、計算上の伝播距離と実際の伝播距離とにずれが生じる。従って、検出信号から伝播距離を推定する場合には、音速差に起因する伝搬距離のあいまいさが残る。生体内の伝播距離のあいまいさを伝播時間のあいまいさと捉えると、時刻ｔ＝０も曖昧性を持つことになり、ｔ＝０の圧力分布というのもあいまいになり、定義が困難となる。ｔ＝０の分布が吸収分布であるのに対し、ｔ＞０の分布は伝播時の圧力分布であり、これらが混ざると吸収分布とは言えなくなる。

　本発明は、上記に鑑み、実用的な装置においても検出信号から吸収分布を求めることができる光音響信号処理装置及び方法を提供することを目的とする。

　上記課題を解決するために、本発明は、被検体に対する出射光により被検体内で生じた光音響信号をサンプリングするサンプリング手段と、サンプリングされた光音響信号から被検体への出射光の光強度の時間波形の微分波形である光微分波形をデコンボリューションした信号を生成する光微分波形逆畳込み手段とを備えたことを特徴とする光音響信号処理装置を提供する。

　本発明の光音響信号処理装置は、被検体への出射光の光強度の時間波形の微分波形を取得する光微分波形取得手段を更に備えていてもよい。

　本発明では、光微分波形逆畳込み手段が、サンプリングされた光音響信号をフーリエ変換する第１のフーリエ変換手段と、光微分波形を所定のサンプリングレートでサンプリングした信号をフーリエ変換する第２のフーリエ変換手段と、フーリエ変換された光微分波形の逆数を逆フィルタとして求める逆フィルタ演算手段と、フーリエ変換された光音響信号に逆フィルタを適用するフィルタ適用手段と、逆フィルタが適用された光音響信号をフーリエ逆変換するフーリエ逆変換手段とを有する構成を採用することができる。

　光音響信号が第１のサンプリングレートでサンプリングされ、光微分波形が第１のサンプリングレートよりも高い第２のサンプリングレートでサンプリングされており、第１のサンプリングレートでサンプリングされた光音響信号を、第２のサンプリングレートでリサンプルするリサンプル手段を更に備え、第１のフーリエ変換手段が、リサンプル手段でリサンプルされた光音響信号をフーリエ変換するものとしてもよい。この場合、第１のフーリエ変換手段と第２のフーリエ変換手段とが、おなじデータ点数でフーリエ変換を行うこととしてもよい。

　上記に代えて、光音響信号が第１のサンプリングレートでサンプリングされ、光微分波形が第１のサンプリングレートよりも高い第２のサンプリングレートでサンプリングされ、かつ、第１のフーリエ変換手段が第１のデータ点数でフーリエ変換を行い、第２のフーリエ変換手段が、第１のデータ点数よりも多い第２のデータ点数でフーリエ変換を行うものであり、フーリエ変換された光音響信号に対して、第１のデータ点数と第２のデータ点数との差の分だけ中央に０を付加するゼロパディングを行うゼロパディング手段を更に備え、フィルタ適用手段が、ゼロパディング手段でゼロパディングが施された信号に対して逆フィルタを適用するものとしてもよい。

　あるいは、光音響信号が第１のサンプリングレートでサンプリングされ、光微分波形が第１のサンプリングレートよりも高い第２のサンプリングレートでサンプリングされ、かつ、第１のフーリエ変換手段が第１のデータ点数でフーリエ変換を行い、第２のフーリエ変換手段が、第１のデータ点数よりも多い第２のデータ点数でフーリエ変換を行うものであり、フーリエ変換された光微分波形から、第１のデータ点数と第２のデータ点数の差の分だけ高周波成分サンプル点を除去する高周波成分サンプル点除去手段を更に備え、逆フィルタ演算手段が、光微分波形をフーリエ変換した信号から高周波成分サンプル点を除去した信号の逆数を逆フィルタとして求めるものであることとしてもよい。

　第２のデータ点数は、第１のデータ点数に、第２のサンプリングレートと第１のサンプリングレートとの比を乗じたデータ点数と等しいか、又はそれよりも多いデータ点数であるものとすることができる。

　本発明の光音響信号処理装置では、複数の検出素子で検出され、サンプリング手段でサンプリングされた光音響信号に基づいて光音響信号を再構成する光音響信号再構成手段を更に備え、第１のフーリエ変換手段が、光音響信号再構成手段で再構成された光音響信号をフーリエ変換する構成を採用することが好ましい。

　光微分波形をデコンボリューションした信号から、光音響信号を検出する検出器の受信角度依存特性の影響を除去するように、光微分波形をデコンボリューションした信号を補正する補正手段を更に備える構成をとしてもよい。

　また、補正手段が、検出器の受信角度依存特性に代えて、又は加えて、光微分波形をデコンボリューションした信号から被検体における光の入射光分布の影響を除去するように、光微分波形をデコンボリューションした信号を補正してもよい。

　光微分波形逆畳込み手段が、デコンボリューション時に、更に、ノイズ増幅周波数帯をフィルタリングすることとしてもよい。

　本発明では、被検体に対する出射光が複数の波長の光を含み、サンプリング手段が各波長の光に対応した光音響信号をサンプリングし、光微分波形逆畳込み手段が各波長の光に対応した光音響信号から光微分波形をデコンボリューションした信号を生成するものであってもよく、その場合、各波長の光に対応した、光音響信号から光微分波形をデコンボリューションした信号を演算処理する２波長データ演算手段を更に備える構成とすることができる。

　本発明の光音響信号処理装置は、光微分波形をデコンボリューションした信号に基づいて光音響画像を生成する光音響画像生成手段を更に備える構成を採用することができる。

　サンプリング手段が、被検体に対して送信された音響波に対する反射音響波をサンプリングするものであり、サンプリングされた反射音響波に基づいて反射音響波画像（超音波画像）を生成する反射音響波画像生成手段と、光音響画像と反射音響波画像とを合成する画像合成手段とを更に備える構成を採用してもよい。

　画像合成手段が、光音響画像と反射音響波画像とを重畳することで画像合成を行うこととしてもよい。

　本発明は、また、被検体に対する出射光に起因して被検体内で発生した光音響信号を検出するステップと、検出された光音響信号から、被検体への出射光の光強度の時間波形の微分波形である光微分波形をデコンボリューションするステップとを有する光音響信号処理方法を提供する。

　本発明の光音響信号処理装置では、デコンボリューションするステップが、検出された光音響信号をフーリエ変換するステップと、光微分波形をフーリエ変換するステップと、フーリエ変換された光微分波形の逆数を逆フィルタとして求めるステップと、フーリエ変換された光音響信号に逆フィルタを適用するステップと、逆フィルタが適用された光音響信号をフーリエ逆変換するステップとを含むものであってもよい。

　本発明の光音響信号処理装置及び方法では、被検体内で生じた光音響信号から、被検体に照射された光の光強度の時間波形を微分した光微分波形をデコンボリューションする。このようにすることで、光微分項がコンボリューションされた検出信号から光微分項を除去することができ、検出信号から吸収分布を求めることができる。そのような吸収分布を画像化した場合には、吸収分布画像を示す光音響画像が得られる。

光パルス微分波形デコンボリューションの基本アルゴリズムを示すブロック図。再構成後の光音響信号を示す波形図。ＦＦＴ後の光音響信号ＦＦＴを示す波形図。光パルス微分波形（ｈ）を示す波形図。ＦＦＴ後の光パルス微分波形ＦＦＴ（ｆｆｔ＿ｈ）を示す波形図。光パルス微分波形ＦＦＴフィルタを示す波形図。デコンボリューション後のＦＦＴ波形を示す波形図。逆変換された光音響信号を示す波形図。再構成後の光音響信号に基づいて生成した光音響画像を示す図。デコンボリューション後の光音響信号に基づいて生成した光音響画像を示す図。本発明の第１実施形態の光音響信号処理装置を含む光音響画像生成装置を示すブロック図。光微分波形逆畳込み手段を示すブロック図。光音響画像生成の動作手順を示すフローチャート。本発明の第２実施形態における光微分波形逆畳込み手段を示すブロック図。４００ＭＨｚのサンプリングレートでサンプリングした光パルス微分波形を示す波形図。４０ＭＨｚのサンプリングレートでサンプリングした光パルス微分波形を示す波形図。本発明の第３実施形態における光微分波形逆畳込み手段を示すブロック図。光音響信号（周波数領域）を示すグラフ。ゼロパディング後の光音響信号を示すグラフ。本発明の第４実施形態における光微分波形逆畳込み手段を示すブロック図。光パルス微分波形（周波数領域）を示すグラフ。高周波成分サンプル点が除去された光パルス微分波形を示すグラフ。本発明の第５実施形態の光音響画像生成装置を示すブロック図。本発明の第６実施形態の光音響画像生成装置を示すブロック図。

　本発明の実施の形態の説明に先立って、本発明の概要について説明する。光吸収体であるミクロ吸収粒子を考え、このミクロ吸収粒子がパルスレーザ光を吸収して圧力波（光音響圧力波）が生じることを考える。時刻をｔとして、位置ｒにあるあるミクロ吸収粒子から発生する光音響圧力波を、位置Ｒで観測した場合の圧力波形ｐ_{ｍｉｃｒｏ}（Ｒ，ｔ）は、[Phys. Rev. Lett. 86(2001)3550.]より、以下の球面波となる。

　ここで、Ｉ（ｔ）は励起光の光強度の時間波形であり、係数ｋは、粒子が光を吸収して音響波を出力する際の変換係数であり、ｖ_ｓは被検体の音速である。また、位置ｒ、Ｒは、空間上の位置を示すベクトルである。ミクロ吸収粒子から発生する圧力は、上記式に示すように、光パルス微分波形に比例した球面波となる。

　実際にイメージングする対象から得られる圧力波形は、よりマクロな吸収体のサイズを有しているため、上記のミクロ吸収波形を重ね合わせた波形になると考える（重ね合わせの原理）。ここで、マクロな光音響波を発する粒子の吸収分布をＡ（ｒ－Ｒ）とし、そのマクロな吸収体からの圧力の観測波形をｐ_{ｍａｃｒｏ}（Ｒ，ｔ）とする。観測位置Ｒでは、各時刻において、観測位置Ｒから半径ｖ_ｓｔに位置する吸収粒子からの光音響波が観測されることになるため、観測波形ｐ_{ｍａｃｒｏ}（Ｒ，ｔ）は、以下の圧力波形の式で示される。

　上記式（１）からわかるように、観測波形は、光パルス微分のコンボリューション型を示す。

　非特許文献２では、上記式に更に装置インパルス応答をコンボリューションした式を基本に考え、観測波形ｐ_{ｍａｃｒｏ}から、光微分とシステム応答とを不可分な状態で含むｐ_ｄ０をデコンボリューションした後に、吸収分布Ａ（ｒ－Ｒ）の再構成をフィルタ補正逆投影法で行うことを提案している。非特許文献２では、光パルス微分の影響よりも装置インパルス応答を考慮することを重視しており、それゆえに、装置として十分Ｓ／Ｎ比（Signal to Noise Ratio）が取れていない周波数帯域まで強調され、処理後の画像ノイズが増加する。そのため、非特許文献２では、高周波フィルタを含めて処理することが必要になる。

　非特許文献２のように帯域が広い超音波プローブを用いるならば上記の方法でもよい。しかし、実用的な狭帯域プローブを用いた場合には、装置のインパルス応答に対して検出する超音波信号の周波数が低いため、通常の超音波プローブで検出する信号（低周波）に対してデコンボリューションする波形の帯域が広くなり、適切にデコンボリューションすることができず、画像として破たんが生じる。そこで、吸収分布を得る上で重要なのは光パルス微分項を考慮することであるため、本発明では、デコンボリューション処理において、光パルス微分項のみを考慮してデコンボリューションを行うこととする。

　更に、本発明では、従来の超音波システムにおいても用いられている圧力分布を求める再構成（ＦＴＡ法、ＤｎＳ法、ＢＰ法など）を適用後に、再構成後の画像がｔ≠０の圧力分布、すなわち圧力波の伝播過程の圧力分布であることを認識の上で、これを吸収分布に変換することを考えた。圧力分布再構成の基本的な考え方としては、検出位置Ｒ＝（ｘ，ｙ，０）の再構成後の圧力分布ｐ_ｒｅｃ（Ｒ，ｔ）は、各時刻ｔにおける、Ｒの検出軸（ｒ－Ｒ）上の｜ｒ－Ｒ｜位置に存在する吸収体から発生する球面波を、周囲の圧電素子の信号も含めて足し合わせてその位置における圧力強度を計算して得られる。従って、検出軸（ｒ－Ｒ）に存在するミクロの吸収体から発生し、伝播する光音響波を重ねあわせたｐ_ｒｅｃ（Ｒ，ｔ）は以下のように表記できる。

このように吸収分布を１次元で考えて良くなることで、上記式のような圧力表記が可能となる。上記式（２）は、検出軸（ｒ－Ｒ）をｚ軸、検出素子からの距離｜ｒ－Ｒ｜をｚとすると、下記のように表すことができる。

　更に、式（３）から積分とは関係ないｘ，ｙは表記を省略し、ｚ軸を時間で表記すると、上記式は、下記式のように表すことができる。

　このように、（ｘ，ｙ，０）に位置する検出素子における１軸（時間軸又はｚ’軸）のコンボリューション表記が可能となる。

　上記式（４）の両辺をフーリエ変換し、周波数軸において、圧力分布のフーリエ係数を光パルスの時間微分のフーリエ係数で割ることで、光パルス微分をデコンボリューションすることができる。

　デコンボリューション後、得られた式を、フーリエ逆変換することで、Ａ（ｘ，ｙ，ｖ_ｓｔ）を求め、吸収分布を画像化することができる。ここで求めたＡ（ｘ，ｙ，ｖ_ｓｔ）には、検出素子受信角度依存性Ｄ（ｘ，ｙ，ｚ）や、プローブ帯域の固有振動が重畳されている可能性がある。例えば、装置関数Ｄ（ｘ，ｙ，ｚ）を事前に求めておいた上でその逆数をＡ（ｘ，ｙ，ｖ_ｓｔ）に掛けることで、検出素子受信角度依存性の影響を除去できる。また、帯域の固有振動に関してはヒルベルト変換や直交検波処理により強度画像化すれば、その影響を除去することができる。更に、検体に入射する光空間分布Ｌ（ｘ，ｙ，ｚ）を別に観測やシミュレーションにより求め、吸収係数に比例する画素値μ（ｘ，ｙ，ｚ）を、μ（ｘ，ｙ，ｚ）＝Ａ（ｘ，ｙ，ｖ_ｓｔ）／Ｌ（ｘ，ｙ，ｚ）により求めてもよい。この場合、より生体組織と密接な関係のある物理量である吸収係数の分布画像を得ることができる。

　図１は、光パルス微分波形デコンボリューションの基本アルゴリズムを示す。再構成後の光音響信号を入力し、再構成後の光音響信号をＦＦＴ（Fast Fourier Transform）によりフーリエ変換する（ステップＳ１０１）。図２Ａに再構成後の光音響信号を示し、図２ＢにＦＦＴ後の光音響信号ＦＦＴを示す。フーリエ変換することで、図２Ａに示す時間領域の信号が、図２Ｂに示すような周波数領域の信号に変換される。なお、図２Ｂでは、光音響信号ＦＦＴの絶対値を示しているが、実際の処理では複素数のまま処理される。

　光パルス微分波形ｈをＦＦＴによりフーリエ変換する（ステップＳ１０２）。図２Ｃに光パルス微分波形（ｈ）を示し、図２ＤにＦＦＴ後の光パルス微分波形ＦＦＴ（ｆｆｔ＿ｈ）を示す。フーリエ変換することで、図２Ｃに示す時間領域の信号（波形）が、図２Ｄに示す周波数領域の信号に変換される。なお、図２Ｃにおける黒丸は、光パルス微分波形におけるサンプリング点を表している。また、図２Ｄでは、光パルス微分波形ＦＦＴの絶対値を示しているが、実際の処理では複素数のまま処理される。

　ステップＳ１０２で得られたＦＦＴ後の光パルス微分波形ＦＦＴ（ｆｆｔ＿ｈ）の逆数を、光パルス微分波形ＦＦＴフィルタ（逆フィルタ）として求める（ステップＳ１０３）。光パルス微分波形ＦＦＴフィルタは、具体的にはconj(fft_h)/abs(fft_h)²で求めることができる。ここで、conj(fft_h)はfft_hの共役複素数、abs(fft_h)はfft_hの絶対値を表す。図２Ｅに、光パルス微分波形ＦＦＴフィルタを示す。図２Ｄに示す光パルス微分波形ＦＦＴの逆数を求めることで、図２Ｅに示すような光パルス微分波形ＦＦＴフィルタを得ることができる。

　上記のようにして求めた光パルス微分ＦＦＴフィルタと、再構成後の光音響信号ＦＦＴとを要素ごとに乗算し、光音響信号ＦＦＴから光パルス微分波形をデコンボリューションする（ステップＳ１０４）。図２Ｆに、デコンボリューション後のＦＦＴ波形を示す。図２Ｂに示す光音響信号ＦＦＴと図２Ｅに示す光パルス微分波形ＦＦＴフィルタとの乗算を行うことで、図２Ｆに示すＦＦＴ波形が得られる。

　ステップＳ１０４で光パルス微分波形をデコンボリューションしたＦＦＴ波形を、逆ＦＦＴによりフーリエ逆変換し、周波数領域の信号を時間領域の信号に戻す（ステップＳ１０５）。図２Ｇは、逆変換された光音響信号を示す。図２Ｆに示すＦＦＴ波形（周波数領域の信号）を逆ＦＦＴすることで、図２Ｇに示すデコンボリューション後の光音響信号（時間領域の信号）が得られる。このデコンボリューション後の光音響信号は、光吸収分布に光パルス微分波形（図２Ｃ）がコンボリューションされた再構成後の光音響信号（図２Ａ）から、光パルス微分波形をデコンボリューションした吸収分布に相当する。

　図３Ａに、再構成後の光音響信号（図２Ａ）に基づいて生成した光音響画像を示し、図３Ｂに、デコンボリューション後の光音響信号（図２Ｇ）に基づいて生成した光音響画像を示す。図３Ａに示す、再構成後の光音響信号に基づいて生成した光音響画像は、実質的に圧力分布を画像化したものであり、１本の血管が二重に表示されるなど、画像判定上、血管の位置が確認しづらい。これに対し、図３Ｂに示すデコンボリューション後の光音響信号に基づいて生成した光音響画像は、光パルス微分波形をデコンボリューションしていることで吸収体の分布を画像化できており、血管の位置を確認しやすくなっている。

　以下、図面を参照し、本発明の実施の形態を詳細に説明する。図４は、本発明の第１実施形態の光音響信号処理装置を含む光音響画像生成装置を示す。光音響画像生成装置（光音響画像診断装置）１０は、超音波探触子（プローブ）１１、超音波ユニット１２、及び光源ユニット（レーザユニット）１３を備える。

　レーザユニット１３は、被検体に照射するレーザ光を生成する。レーザ光の波長は、観察対象物に応じて適宜設定すればよい。レーザユニット１３が出射するレーザ光は、例えば光ファイバなどの導光手段を用いてプローブ１１まで導光され、プローブ１１から被検体に照射される。プローブ１１は、レーザユニット１３から出射した光が被検体に照射された後に、被検体内の光吸収体がレーザ光を吸収することで生じた超音波（光音響信号）を検出する。プローブ１１は、例えば一次元配列された複数の超音波振動子を有する。

　超音波ユニット１２は、光音響信号処理装置に相当する。超音波ユニット１２は、受信回路２１、ＡＤ変換手段２２、受信メモリ２３、光音響画像再構成手段２４、光微分波形逆畳込み手段２５、補正手段２６、検波・対数変換手段２７、光音響画像構築手段２８、トリガ制御回路２９、及び制御手段３０を有する。受信回路２１は、プローブ１１で検出された光音響信号を受信する。ＡＤ変換手段２２は、サンプリング手段であり、受信回路２１が受信した光音響信号をサンプリングしてデジタル信号に変換する。ＡＤ変換手段２２は、例えば、外部から入力する所定周波数のＡＤクロック信号に基づいて、所定のサンプリング周期で光音響信号をサンプリングする。受信メモリ２３は、ＡＤ変換手段２２でサンプリングされた光音響信号を記憶する。

　光音響画像再構成手段２４は、受信メモリ２３から光音響信号を読み出し、プローブ１１の複数の超音波振動子で検出された光音響信号に基づいて、光音響画像の各ラインのデータを生成する。光音響画像再構成手段２４は、例えばプローブ１１の６４個の超音波振動子からのデータを、超音波振動子の位置に応じた遅延時間で加算し、１ライン分のデータを生成する（遅延加算法）。光音響画像再構成手段２４は、遅延加算法に代えて、ＢＰ法（Back Projection）により再構成を行ってもよい。あるいは光音響画像再構成手段２４は、ハフ変換法又はフーリエ変換法を用いて再構成を行ってもよい。

　光微分波形逆畳込み手段２５は、再構成された光音響信号から被検体に照射された光の光強度の時間波形の微分波形である光パルス微分波形（光微分波形）をデコンボリューションした信号を生成する。光パルス微分波形をデコンボリューションすることで、ｔ≠０に再構成した圧力分布から、ｔ＝０に再構成した圧力分布、すなわち吸収分布を求めることができる。光微分波形逆畳込み手段２５は、再構成前の光音響信号に対してデコンボリューションを行ってもよい。

　補正手段２６は、光パルス微分波形がデコンボリューションされた信号を補正し、光パルス微分波形がデコンボリューションされた信号から、プローブ１１における超音波振動子の受信角度依存特性の影響を除去する。また、補正手段２６は、受信角度依存特性に加えて、又はこれらに代えて、光パルス微分波形がデコンボリューションされた信号から被検体における光の入射光分布の影響を除去する。補正手段２６を省き、これらの補正を行わずに、光音響画像の生成を行ってもよい。

　検波・対数変換手段２７は、補正後の各ラインのデータの包絡線を求め、求めた包絡線を対数変換する。包絡線を求める検波手段としては、ヒルベルト変換や直交検波など従来から用いられている手法を用いることができる。これにより、超音波振動子の固有振動による帯域の影響が除去できる。光音響画像構築手段２８は、対数変換が施された各ラインのデータに基づいて、光音響画像を生成する。光音響画像構築手段２８は、例えば光音響信号（ピーク部分）の時間軸方向の位置を光音響層画像における深さ方向の位置に変換して光音響画像を生成する。

　制御手段３０は、超音波ユニット１２内の各部を制御する。トリガ制御回路２９は、光音響画像生成に際して、レーザユニット１３にフラッシュランプトリガ信号を送る。また、フラッシュランプトリガ信号の出力後に、Ｑスイッチトリガ信号を送る。レーザユニット１３は、フラッシュランプ３１とＱスイッチ３２とを含む。レーザユニット１３は、フラッシュランプトリガ信号を受けてフラッシュランプ３１を点灯し、レーザ励起を開始する。レーザユニット１３は、Ｑスイッチトリガ信号が入力されるとＱスイッチをＯＮにし、パルスレーザ光を出射する。トリガ制御回路２９は、被検体に対するレーザ光照射と同期してＡＤ変換手段２２にサンプリングトリガ信号を送り、ＡＤ変換手段２２における光音響信号のサンプリング開始タイミングを制御する。

　図５に、光微分波形逆畳込み手段２５を示す。光微分波形逆畳込み手段２５は、フーリエ変換手段４１、４２と、逆フィルタ演算手段４３と、フィルタ適用手段４４と、フーリエ逆変換手段４５と、光微分波形取得手段５１とを有する。光微分波形取得手段５１は、光パルス微分波形を取得する。光微分波形取得手段５１は、例えばメモリから光パルス微分波形を読み出す。これに代えて、メモリから被検体に照射された光の光強度の時間波形を読み出し、それを時間微分してもよい。また、被検体に照射された光の光強度の時間波形を測定し、その測定結果を時間微分することで光パルス微分波形を取得してもよい。更には、被検体に照射されたパルス光の光強度の時間波形又はその微分波形を表す関数であって、パルス光のパルス時間幅を独立変数とする関数と、被検体に照射されたパルス光のパルス時間幅の測定結果とを用いて、光パルス微分波形を取得するようにしてもよい。

　フーリエ変換手段（第１のフーリエ変換手段）４１は、離散フーリエ変換により、再構成された光音響信号を時間領域の信号から周波数領域の信号へと変換する。フーリエ変換手段（第２のフーリエ変換手段）４２は、離散フーリエ変換により、光パルス微分波形を所定のサンプリングレートでサンプリングした信号を時間領域の信号から周波数領域の信号へと変換する。フーリエ変換のアルゴリズムにはＦＦＴを用いることができる。

　本実施形態においては、光音響信号のサンプリングレートと光パルス微分波形のサンプリングレートとは等しいものとする。例えば光音響信号はＦｓ＝４０ＭＨｚのサンプリングクロックに同期してサンプリングされており、光微分パルスも、Ｆｓ＿ｈ＝４０ＭＨｚのサンプリングレートでサンプリングされている。フーリエ変換手段４１は、４０ＭＨｚでサンプリングされた光音響信号を、例えば１０２４点のフーリエ変換でフーリエ変換する。また、フーリエ変換手段４２は、４０ＭＨｚでサンプリングされた光パルス微分波形を１０２４点のフーリエ変換でフーリエ変換する。

　逆フィルタ演算手段４３は、フーリエ変換された光パルス微分波形の逆数を逆フィルタとして求める。例えば逆フィルタ演算手段４３は、光パルス微分波形ｈをフーリエ変換した信号をｆｆｔ＿ｈとしたとき、conj(fft_h)/abs(fft_h)²を逆フィルタとして求める。フィルタ適用手段４４は、フーリエ変換手段４１でフーリエ変換された光音響信号に、逆フィルタ演算手段４３で求められた逆フィルタを適用する。フィルタ適用手段４４は、例えば、要素ごとに、光音響信号のフーリエ係数と逆フィルタのフーリエ係数とを乗算する。逆フィルタが適用されることで、周波数領域の信号において、光パルス微分波形がデコンボリューションされる。フーリエ逆変換手段４５は、フーリエ逆変換により、逆フィルタが適用された光音響信号を、周波数領域の信号から時間領域の信号へと変換する。フーリエ逆変換により、時間領域の吸収分布信号が得られる。

　図６は、動作手順を示す。トリガ制御回路２９は、レーザユニット１３に対してフラッシュランプトリガ信号を出力する。レーザユニット１３は、フラッシュランプトリガ信号を受けてフラッシュランプ３１を点灯する。トリガ制御回路２９は、所定のタイミングでＱスイッチトリガ信号を出力する。レーザユニット１３は、Ｑスイッチトリガ信号が入力されると、Ｑスイッチ３２をＯＮにし、パルスレーザ光を出射する。出射したパルスレーザ光は、例えばプローブ１１まで導光され、プローブ１１から被検体に照射される（ステップＳ１）。

　プローブ１１は、レーザ光の照射後、レーザ光の照射により被検体内で発生した光音響信号を検出する（ステップＳ２）。超音波ユニット１２の受信回路２１は、プローブ１１で検出された光音響信号を受信する。トリガ制御回路２９は、被検体に対する光照射のタイミングに合わせてＡＤ変換手段２２にサンプリングトリガ信号を送る。ＡＤ変換手段２２は、サンプリングトリガ信号を受けて光音響信号のサンプリングを開始し、光音響信号のサンプリングデータを受信メモリ２３に格納する。

　光音響画像再構成手段２４は、受信メモリ２３から光音響信号のサンプリングデータを読み出し、読み出した光音響信号のサンプリングデータに基づいて、光音響信号を再構成する（ステップＳ３）。光微分波形逆畳込み手段２５は、再構成された光音響信号から、被検体に照射されたパルスレーザ光の光強度の時間波形を微分した光パルス微分波形をデコンボリューションする（ステップＳ４）。このデコンボリューションにより、吸収分布を示す光音響信号が得られる。

　補正手段２６は、光パルス微分波形がデコンボリューションされた信号を、検出素子受信角度依存性や被検体における光の入射分布で補正する。検波・対数変換手段２７は、補正手段２６で補正された光音響信号の包絡線を求め、求めた包絡線を対数変換する。光音響画像構築手段２８は、対数変換が施された各ラインのデータに基づいて、光音響画像を生成する（ステップＳ５）。この光音響信号は、吸収分布を画像化した吸収分布画像である。画像表示手段１４は、表示画面上に、吸収分布画像である光音響画像を表示する（ステップＳ６）。

　本実施形態では、ひとまず、光音響画像再構成手段２４にて、通常の再構成法により発光時刻（ｔ＝０）の圧力分布として光音響信号（光音響画像）を再構成する。次に、光の発光時間は実際には有限の長さを持つことから、再構成時にｔ＝０としていた時刻を有限の長さの時間と考え、光微分波形逆畳込み手段２５にて、再構成後の光音響画像から光パルス微分波形をデコンボリューションする。光パルス微分波形をデコンボリューションすることで、吸収分布を得ることができ、吸収分布画像を生成することができる。このような手法を採用することで、実用的な光パルス幅と実用的な超音波システム、或いは実際の生体を観測した場合でも、吸収分布を画像化することができる。これは、現状システムの検出器の帯域やＡＤサンプリングを使用できる利点がある。また、本実施形態においては光音響画像の再構成で圧力分布を一度出しているため、既存の超音波アルゴリズム、装置との親和性が高い。

　非特許文献２との比較では、非特許文献２では、光微分関数と装置インパルス応答関数とが不可分な状態でデコンボリューションしている。狭帯域の通常の超音波プローブ、例えば８ＭＨｚの超音波プローブを用いた場合、その超音波プローブにより４～１２ＭＨｚの信号を検出できるが、帯域の端となる４ＭＨｚや１２ＭＨｚでは感度が低いため、４ＭＨｚや１２ＭＨｚのＳ／Ｎは、８ＭＨｚのＳ／Ｎとよりも低い。非特許文献２では、装置インパルス応答を考慮（補正）することを重視しており、感度の低い４ＭＨｚや１２ＭＨｚの信号が強調され、結果としてＳ／Ｎの悪い周波数成分が強調された画像が生成される。これに対し、本願発明では、光パルス微分波形のみをデコンボリューションしており、光パルス微分波形でＳ／Ｎを下げずに処理しつつ、装置インパルス応答に相当する成分をデコンボリューション後に除去することができ、Ｓ／Ｎを下げずに、光音響画像を生成できる。

　次いで、本発明の第２実施形態を説明する。第１実施形態では、光音響信号のサンプリングレートと光パルス微分波形のサンプリングレートとが一致しており、双方の信号を同じデータ点数でフーリエ変換した。本実施形態では、光音響信号を低速サンプリングする一方で、光パルス微分波形を高速サンプリングする。つまり、光パルス微分波形のサンプリングレートを光音響信号のサンプリングレートよりも高く設定する。例えば光音響信号のサンプリング間隔（サンプリングレートの逆数）は、被検体に照射される光のパルス時間幅よりも長く設定される。フーリエ変換に際しては、低サンプリングレートでサンプリングされた光音響信号を、光パルス微分波形のサンプリングレートと同じサンプリングレートでリサンプル（アップサンプル）した上で、フーリエ変換を行う。その他の点は、第１実施形態と同様でよい。

　図７は、本実施形態における光微分波形逆畳込み手段２５ａを示す。本実施形態における光微分波形逆畳込み手段２５ａは、図５に示す第１実施形態における光微分波形逆畳込み手段２５の構成に加えて、リサンプル手段４６及び４７を有する。リサンプル手段４６は、アップサンプル手段であり、低いサンプリングレートでサンプリングされた光音響信号のサンプリングデータを、光パルス微分波形のサンプリングレートと同じサンプリングレートでアップサンプルする。リサンプル手段４６は、例えば、低サンプリングレートでサンプリングされた光音響信号のサンプル点間にゼロを付加し、アップサンプル前のナイキスト周波数でカットするローパスフィルタをかけることでアップサンプルを行う。

　例えば、ＡＤ変換手段２２（図４）における光音響信号のサンプリングレート（第１のサンプリングレート）が４０ＭＨｚであり、光パルス微分波形のサンプリングレート（第２のサンプリングレート）が４００ＭＨｚであったとする。この場合、リサンプル手段４６は、４０ＭＨｚの光音響信号を４００ＭＨｚの信号にアップサンプルする。フーリエ変換手段４１は、リサンプル手段４６でアップサンプルされた光音響信号をフーリエ変換する。光音響信号をフーリエ変換するフーリエ変換手段４１と、光パルス微分波形をフーリエ変換するフーリエ変換手段４２とは、同じデータ点数でフーリエ変換を行う。例えばフーリエ変換手段４１は光音響信号を８１９２点の周波数領域の信号に変換し、フーリエ変換手段４２は光パルス微分波形を８１９２点の周波数領域の信号に変換する。

　フィルタ適用手段４４は、アップサンプルされた光音響信号をフーリエ変換した信号に対して逆フィルタを適用する。フーリエ逆変換手段４５は、逆フィルタが適用された信号を、周波数領域の信号から時間領域の信号（吸収分布）へと変換する。時間領域の信号に戻された吸収分布信号は、例えば４００ＭＨｚにアップサンプルされた状態の信号となっている。リサンプル手段４７は、吸収分布信号が、光音響信号の元のサンプルリングレートでサンプリングされた信号となるように、吸収信号をダウンサンプルする。リサンプル手段４７は、例えば４００ＭＨｚの吸収信号を４０ＭＨｚの吸収信号にダウンサンプルする。ダウンサンプリングは、例えばダウンサンプル後のナイキスト周波数でカットするローパスフィルタをかけた後に、サンプル点を間引くことで行う。

　図８Ａに、４００ＭＨｚのサンプリングレートでサンプリングした光パルス微分波形を示し、図８Ｂに、４０ＭＨｚのサンプリングレートでサンプリングした光パルス微分波形を示す。サンプリングレート４００ＭＨｚでは、図８Ａに示すように、光パルス微分波形を正確に再現できる。一方、光パルス微分波形のサンプリングレートを光音響信号のサンプリングレートに合わせ、４０ＭＨｚでサンプリングすると、図８Ｂに示すように、光パルス微分波形を正確に再現できなくなる。

　フィルタ適用手段４４にて光音響信号をフーリエ変換した信号に逆フィルタを適用する際には、双方のデータ点数が揃っている必要がある。光音響信号のサンプリングレートに合わせて光パルス微分波形のサンプリングレートを設定すると、図８Ｂに示したように、波形変化に対してサンプリング周波数が低すぎ、光パルス微分波形が正確に再現できない。このような光パルス微分波形から求めた逆フィルタを適用した場合、光パルス微分項を正確にデコンボリューションできずに、吸収分布を正しく求められないこともある。

　一方、光パルス微分波形を正確に再現するために光パルス微分波形のサンプリングレートを例えば４００ＭＨｚに設定し、光音響信号のサンプリングレートを４００ＭＨｚに合わせるとした場合は、光パルス微分項を正確にデコンボリューションでき、吸収分布を正しく求めることができる。しかしながら、その場合、ＡＤ変換手段２２には高速なＡＤ変換器が要求され、また、サンプリングデータの総数が増えることから、受信メモリ２３（図４）に要求されるメモリ容量が増大する。更に、光音響画像再構成手段２４で取り扱うデータが増えるため、再構成に要する時間も長くなる。

　本実施形態では、リサンプル手段４６で、事後的に光音響信号のサンプリングデータをリサンプルする。本実施形態では、検出後の光音響信号を信号処理でアップサンプルしているため、光音響の検出から再構成までは低速サンプリングしつつも、光パルス微分項を正確にデコンボリューションすることができる。本実施形態では、ＡＤ変換手段２２に高速なＡＤ変換器は不要であり、受信メモリ２３に必要なメモリ容量も増大しない。また、光音響信号の再構成に要する時間も増大せず、光音響信号の検出時に高いサンプリングレートでサンプリングする場合に比して、処理時間を短縮することができる。

　続いて、本発明の第３実施形態を説明する。本実施形態では、第２実施形態と同様に、光パルス微分波形のサンプリングレートを光音響信号のサンプリングレートよりも高く設定する。第２実施形態では、低サンプリングレートでサンプリングされた光音響信号をアップサンプルし、双方の信号を同じデータ点数でフーリエ変換した。本実施形態では、光パルス微分波形のフーリエ変換を、光音響信号のフーリエ変換のデータ点数よりも多いデータ点数で行い、フーリエ変換された光音響信号に対して、データ点数の差の分だけ中央（高周波成分領域）にゼロ点を付加する。その他の点は、第１実施形態と同様でよい。

　図９は、本実施形態における光微分波形逆畳込み手段２５ｂを示す。本実施形態における光微分波形逆畳込み手段２５ｂは、図５に示す第１実施形態における光微分波形逆畳込み手段２５の構成に加えて、ゼロパディング手段４８とゼロ点除去手段４９とを有する。例えば、光音響信号のサンプリングレート（第１のサンプリングレート）は４０ＭＨｚであり、光パルス微分波形のサンプリングレート（第２のサンプリングレート）は３２０ＭＨｚであるとする。フーリエ変換手段４１は、例えば４０ＭＨｚの光音響信号を１０２４点（第１のデータ点数）の周波数領域の信号に変換し、フーリエ変換手段４２は、３２０ＭＨｚの光パルス微分波形を８１９２点（第２のデータ点数）の周波数領域の信号に変換する。第２のデータ点数は、第１のデータ点数に、第２のサンプリングレートと第１のサンプリングレートとの比を乗じたデータ点数と等しいか、又はそれよりも多い。

　ゼロパディング手段４８は、フーリエ変換手段４１から周波数領域の信号に変換された光音響信号を入力する。ゼロパディング手段４８は、フーリエ変換された光音響信号に対して、フーリエ変換後の光音響信号と光パルス微分波形のデータ点数の差の分だけ中央にゼロ点（信号値ゼロの点）を付加する。ゼロパディング手段４８は、例えばデータ点数１０２４点の光音響信号（周波数領域）を、周波数帯域の中心周波数で２つに分割し、分割した２つの周波数領域の間にデータ点数の差の分だけゼロ点を付加し、光パルス微分波形（周波数領域）のデータ点数と同じデータ点数８１９２点の光音響信号を生成する。ゼロ点の付加は、周波数領域におけるアップサンプリングに相当する。

　フィルタ適用手段４４は、ゼロパディング手段４８でゼロパディングが施された信号に対して逆フィルタを適用する。ゼロ点除去手段４９は、逆フィルタが適用された信号からゼロパディング手段４８で“０”が付加された周波数帯域を除去する。例えばゼロパディング手段４８にてデータ点数１０２４点の光音響信号（周波数領域）がデータ点数８１９２点の信号に変換されていたとき、ゼロ点除去手段４９は、フィルタ適用後の信号（データ点数８１９２点）をデータ点数１０２４点の信号に戻す。ゼロ点の除去は、周波数領域におけるダウンサンプリングに相当する。フーリエ逆変換手段４５は、データ点数１０２４点に戻された信号を、周波数領域の信号から時間領域の信号へと変換する。

　図１０Ａに、フーリエ変換された光音響信号を示し、図１０Ｂに、ゼロパディング後の光音響信号を示す。例えば、ＡＤ変換手段２２（図４）における光音響信号のサンプリングレートが４０ＭＨｚであるとき、その光音響信号をフーリエ変換した信号は、図１０Ａに示すように、０ＭＨｚから４０ＭＨｚまでの周波数帯域の信号となる。この信号を、ナイキスト周波数（サンプリング周波数の１／２）である２０ＭＨｚを境に２つの領域Ａ、Ｂに２分割する。ゼロパディング手段４８は、図１０Ｂに示すように、２つの領域の間にゼロ点を８１９２－１０２４＝７１６８個挿入する。ゼロ点が付加された結果、領域Ｂの信号は、３００ＭＨｚから３２０ＭＨｚの周波数領域に対応した信号となる。

　本実施形態では、低サンプリングレートでサンプリングされた光音響信号を周波数領域の信号に変換し、変換された周波数領域の信号の高周波成分の領域のゼロ点を付加する。本実施形態と第２実施形態との相違点は、第２実施形態では、光音響信号をアップサンプルするのに対し、本実施形態では、光音響信号を周波数領域でアップサンプルする点である。時間領域に代え、周波数領域において、双方の信号の帯域差を埋めるようにリサンプル（アップサンプル）を行う場合も、第２実施形態と同様に、光音響の検出から再構成までは低速サンプリングしつつも、光パルス微分項を正確にデコンボリューションすることができる。

　引き続いて、本発明の第４実施形態を説明する。本実施形態においても、第２及び第３実施形態と同様に、光パルス微分波形のサンプリングレートを光音響信号のサンプリングレートよりも高く設定する。本実施形態では、光パルス微分波形を、光音響信号のフーリエ変換のデータ点数よりも多いデータ点数で行い、フーリエ変換された光微分波形から高周波成分サンプル点を除去し、その逆数を逆フィルタとして求める。その他の点は、第１実施形態と同様でよい。

　図１１は、本実施形態における光微分波形逆畳込み手段２５ｃを示す。本実施形態における光微分波形逆畳込み手段２５ｃは、図５に示す第１実施形態における光微分波形逆畳込み手段２５の構成に加えて、高周波成分サンプル点除去手段５０を有する。例えば、光音響信号のサンプリングレート（第１のサンプリングレート）は４０ＭＨｚであり、光パルス微分波形のサンプリングレート（第２のサンプリングレート）は３２０ＭＨｚであるとする。フーリエ変換手段４１は、例えば４０ＭＨｚの光音響信号を１０２４点（第１のデータ点数）の周波数領域の信号に変換し、フーリエ変換手段４２は、３２０ＭＨｚの光パルス微分波形を８１９２点（第２のデータ点数）の周波数領域の信号に変換する。第２のデータ点数は、第１のデータ点数に、第２のサンプリングレートと第１のサンプリングレートとの比を乗じたデータ点数と等しいか、又はそれよりも多い。

　高周波成分サンプル点除去手段５０は、フーリエ変換手段４２から周波数領域の信号に変換された光パルス微分波形を入力する。高周波成分サンプル点除去手段５０は、フーリエ変換された光パルス微分波形から、フーリエ変換後の光音響信号と光パルス微分波形のデータ点数の差の分だけ高周波成分サンプル点を除去する。高周波成分サンプル点除去手段５０は、例えばデータ点数８１９２点の光パルス微分波形（周波数領域）から高周波成分に相当する中央のデータ点を削除し、光音響信号（周波数領域）のデータ点数と同じデータ点数１０２４点の光パルス微分波形を生成する。高周波成分サンプル点の除去は、周波数領域における光パルス微分波形のダウンサンプリングに相当する。

　図１２Ａに、フーリエ変換された光パルス微分波形を示し、図１２Ｂに、高周波成分サンプル点が除去された光パルス微分波形を示す。例えば、光パルス微分波形のサンプリングレートが３２０ＭＨｚであるとき、その光パルス微分波形をフーリエ変換した信号（データ点数８１９２点）は、図１２Ａに示すように、０ＭＨｚから３２０ＭＨｚまでの周波数帯域の信号となる。この信号を、１番目のデータ点から５１２番目までの領域（領域Ａ）、５１３番目のデータ点から７６８０番目のデータ点までの領域（領域Ｂ）、及び、７６８１番目のデータ点から８１９２番目のデータ点までの領域（領域Ｃ）の３つの領域に分け、領域Ｂのデータ点を除去する。図１２Ｂに示すように、領域Ａと領域Ｃとをつなげることで、０ＭＨｚから４０ＭＨｚまでの周波数帯域に対応したデータ点数１０２４点の光パルス微分波形（周波数領域）が得られる。

　逆フィルタ演算手段４３は、高周波成分サンプル点が除去された光パルス微分波形（周波数領域）の逆数を逆フィルタとして求める。逆フィルタ演算手段４３は、例えばデータ点が８１９２点から１０２４点に削減された光パルス微分波形の逆数を逆フィルタとして求める。フィルタ適用手段４４は、例えばデータ点数１０２４点の光音響信号（周波数領域）と逆フィルタとを要素ごとに乗算する。フーリエ逆変換手段４５は、逆フィルタが適用された信号を、周波数領域の信号から時間領域の信号へと変換する。

　ここで、第３実施形態では、フィルタ適用手段４４は、図１０Ｂに示す高周波成分の領域にゼロ点が付加された光音響信号（周波数領域）と、図１２Ａに示す光パルス微分波形（周波数領域）の逆数とを乗算する。光音響信号の高周波成分領域の値は“０”であるため、光パルス微分波形の高周波成分（図１２Ａの領域Ｂ）は、逆フィルタ適用後の光音響信号に影響を与えない。従って、本実施形態のように、光パルス微分波形の周波数領域の信号から高周波成分サンプル点を除去し、高周波成分を除去した光パルス微分波形から逆フィルタを求め、求めた逆フィルタを光音響信号（周波数領域）に適用しても、得られる結果は第３実施形態と同じ結果となる。つまり、本実施形態においても、第３実施形態と同様な効果が得られる。

　続いて、本発明の第５実施形態を説明する。図１３は、本発明の第５実施形態の光音響画像生成装置を示す。本実施形態の光音響画像生成装置１０ａにおける超音波ユニット１２ａは、図４に示す第１実施形態の光音響画像生成装置１０における超音波ユニット１２の構成に加えて、送信制御回路３３、データ分離手段３４、超音波画像再構成手段３５、検波・対数変換手段３６、超音波画像構築手段３７、及び画像合成手段３８を備える。本実施形態の光音響画像生成装置１０ａは、光音響画像に加えて、超音波画像の生成を行う点で、第１実施形態の光音響画像生成装置１０と相違する。なお、本実施形態においては音響波として超音波を用いるが、被検対象や測定条件等に応じて適切な周波数を選択することにより、可聴周波数の音響波であっても良い。また、本実施形態を第２から第４実施形態の何れかと組み合わせ、それら実施形態において超音波画像の生成を行うこととしてもよい。

　本実施形態では、プローブ１１は、光音響信号の検出に加えて、被検体に対する音響波（超音波）の出力（送信）、及び送信した超音波に対する被検体からの反射超音波の検出（受信）を行う。トリガ制御回路２９は、超音波画像（反射音響波画像）の生成時は、送信制御回路３３に超音波送信を指示する旨の超音波送信トリガ信号を送る。送信制御回路３３は、トリガ信号を受けると、プローブ１１から超音波を送信させる。プローブ１１は、超音波の送信後、被検体からの反射超音波を検出する。超音波の送受信は分離してもよい。例えばプローブ１１とは異なる位置から超音波の送信を行い、その送信された超音波に対する反射超音波をプローブ１１で受信してもよい。

　プローブ１１が検出した反射超音波は、受信回路２１を介してＡＤ変換手段２２に入力される。トリガ制御回路２９は、超音波送信のタイミングに合わせてＡＤ変換手段２２にサンプリグトリガ信号を送り、反射超音波のサンプリングを開始させる。ここで、反射超音波はプローブ１１と超音波反射位置との間を往復するのに対し、光音響信号はその発生位置からプローブ１１までの片道である。反射超音波の検出には、同じ深さ位置で生じた光音響信号の検出に比して２倍の時間がかかるため、ＡＤ変換手段２２のサンプリングクロックは、光音響信号サンプリング時の半分、例えば２０ＭＨｚとしてもよい。ＡＤ変換手段２２は、反射超音波のサンプリングデータを受信メモリ２３に格納する。光音響信号の検出（サンプリング）と、反射超音波の検出（サンプリング）とは、どちらを先に行ってもよい。

　データ分離手段３４は、受信メモリ２３に格納された光音響信号のサンプリングデータと反射超音波のサンプリングデータとを分離する。データ分離手段３４は、分離した光音響信号のサンプリングデータを光音響画像再構成手段２４に入力する。光パルス微分波形のデコンボリューションを含む光音響画像（吸収分布画像）の生成は、第１実施形態と同様である。データ分離手段３４は、分離した反射超音波のサンプリングデータを、超音波画像再構成手段３５に入力する。

　超音波画像再構成手段３５は、プローブ１１の複数の超音波振動子で検出された反射超音波（そのサンプリングデータ）に基づいて、超音波画像の各ラインのデータを生成する。各ラインのデータの生成には、光音響画像再構成手段２４における各ラインのデータの生成と同様に、遅延加算法などを用いることができる。検波・対数変換手段３６は、超音波画像再構成手段３５が出力する各ラインのデータの包絡線を求め、求めた包絡線を対数変換する。

　超音波画像構築手段３７は、対数変換が施された各ラインのデータに基づいて、超音波画像を生成する。超音波画像再構成手段３５、検波・対数変換手段３６、及び超音波画像構築手段３７は、反射超音波に基づいて超音波画像を生成する超音波画像生成手段（反射音響波画像生成手段）を構成する。画像合成手段３８は、光音響画像と超音波画像とを合成する。画像合成手段３８は、例えば光音響画像と超音波画像とを重畳することで画像合成を行う。合成された画像は、画像表示手段１４に表示される。画像合成を行わずに、画像表示手段１４に、光音響画像と超音波画像とを並べて表示し、或いは光音響画像と超音波画像とを切り替えてすることも可能である。

　本実施形態では、光音響画像生成装置は、光音響画像に加えて超音波画像を生成する。超音波画像を参照することで、光音響画像では画像化することができない部分を観察することができる。光パルス微分波形をデコンボリューションすることで吸収分布を画像化できる点は、第１実施形態と同様である。また、超音波画像の生成と光音響画像の生成とで、画像再構成や検波・対数変換などのアルゴリズムの大部分を共通化でき、ＦＰＧＡ回路構成やソフトの簡略化が可能であるなどの実用上のメリットを有する。

　引き続いて、本発明の第６実施形態を説明する。図１４は、本発明の第６実施形態の光音響画像生成装置を示す。本実施形態は、被検体に対して複数の波長の光を照射する点で、第１実施形態と相違する。本実施形態の光音響画像生成装置１０ｂにおける超音波ユニット１２ｂは、図４に示す第１実施形態の光音響画像生成装置１０における超音波ユニット１２の構成に加えて、複数波長の光に対する光音響信号（光音響画像）の演算を行う２波長データ演算手段３９を備える。なお、本実施形態を第２から第５実施形態の何れかと組み合わせ、それら実施形態において複数波長の光を照射し、複数波長に対する光音響信号（光音響画像）の演算することとしてもよい。

　本実施形態では、レーザユニット１３は、複数の波長の光を切り替えて出射可能に構成されている。レーザユニット１３は、例えば波長７５０ｎｍのパルスレーザ光と波長８００ｎｍのパルスレーザ光を切り替えて出射する。プローブ１１は、各波長のパルスレーザ光の出射後に被検体からの光音響信号を検出し、受信メモリ２３には、各波長に対応した光音響信号のサンプリングデータが格納される。格納された各波長に対応した光音響信号は、それぞれ光音響画像再構成手段で再構成される。

　光微分波形逆畳込み手段２５は、光音響画像再構成手段２４による再構成後、各波長に対応した光音響信号（光音響画像）から、被検体に照射された各波長の光の光強度の時間波形の微分波形（光微分波形）をそれぞれデコンボリューションする。各波長に対応した光微分波形がデコンボリューションされた光音響信号は、補正手段２６による補正の後、２波長データ演算手段３９にて処理される。

　ここで、生体組織の多くは光吸収特性が光の波長に応じて変わり、また一般に、その光吸収特性も組織ごとに特有のものとなっている。例えば、ヒトの動脈に多く含まれる酸素化ヘモグロビン（酸素と結合したヘモグロビン:oxy-Hb）の波長７５０ｎｍにおける分子吸収係数は波長８００ｎｍにおけるそれも低い。また、静脈に多く含まれる脱酸素化ヘモグロビン（酸素と結合していないヘモグロビンdeoxy-Hb)の波長７５０ｎｍにおける分子吸収係数は波長８００ｎｍにおけるそれよりも高い。この性質を利用し、波長８００ｎｍで得られた光音響信号に対して、波長７５０ｎｍで得られた光音響信号が相対的に大きいのか小さいのかを調べることで、動脈からの光音響信号と静脈からの光音響信号とを判別することができる。

　２波長データ演算手段３９は、例えば複数波長に対応した光音響信号間の相対的な大小関係を比較する。具体的には、２波長データ演算手段３９は、波長７５０ｎｍの光が照射されたときに検出された光音響信号と、波長８００ｎｍの光が照射されたときに検出された光音響信号とを比較し、どちらがどれだけ大きいかを調べる。画像表示に際しては、波長７５０ｎｍの光が照射されたときに検出された光音響信号が大きければ静脈からの光音響信号と判断できるため、その部分を青色で表示するとよい。また、波長８００ｎｍの光が照射されたときに検出された光音響信号が大きければ動脈からの光音響信号と判断できるため、その部分を赤色で表示するとよい。なお、第１実施形態においても説明したように、補正手段２６は省いてもよい。

　本実施形態では、２波長データ演算手段３９は、光微分波形のデコンボリューション後に、複数波長に対応した光音響信号の演算を行う。複数の波長の光を被検体に照射する場合、例えば１つ目の波長の光を被検体に照射した光音響信号を検出した後、２つ目の光を被検体に照射して光音響信号を検出するとき、体動などの影響で、波長間で位置ずれが生じることがある。複数波長間での光音響信号を比較する際には、同じ場所から発生した光音響信号同士を比較することが好ましい。光微分波形をデコンボリューションしない場合は、図３Ａに示したように１本の血管が二重に表示されるなど、画像判定上、血管の位置が確認しづらく、位置ずれ補正がしにくい。光微分波形をデコンボリューションすることで、図３Ｂに示したように光吸収分布を画像化でき、血管の位置が確認しやすくなり、位置ずれ補正も容易となる。

　また、図３Ａに示したように１本の血管が二重に表示される場合は血管の内部に相当する部分には信号が存在せず（信号レベルが所定レベルよりも低い）、複数波長の画像間で位置ずれが生じていると、信号が存在する部分すなわち血管の重複部分が少なくなる。この場合、複数波長間で光音響信号の相対的な大小関係を適切に比較することが難しくなる。これに対し、図３Ｂに示したように光吸収分布を画像化した場合は、血管内の多くの部分が信号が存在する部分となり、多少位置ずれが生じたとしても、複数波長の画像間で多くの部分が重複する。従って、本実施形態において、双方の画像を位置合わせずに比較する場合でも、位置ずれの影響を軽減することができる。

　なお、上記各実施形態では、光音響信号及び光パルス微分波形を周波数領域の信号に変換し、周波数領域でデコンボリューション後に時間領域の信号に戻しているが、これには限定されない。光パルス微分波形のデコンボリューションを時間領域で行うことも可能である。また、光微分波形逆畳込み手段２５は、デコンボリューション時に、光音響信号に対して何らかのフィルタをかける処理をおこなってもよい。例えば光微分波形逆畳込み手段２５が、デコンボリューション時に、ノイズ増幅周波数帯をフィルタリングするようにしてもよい。

　上記各実施形態では、光音響信号から光微分波形をデコンボリューションした後に光音響画像（吸収分布画像）を生成することとしているが、これに加えて、又はこれに代えて、光微分波形をデコンボリューションせずに光音響画像（圧力分布画像）を生成してもよい。例えば、ユーザが、スイッチや表示モニタ上で操作を行うことで、デコンボリューション処理の有無を選択できるようにしておき、ユーザがデコンボリューション処理の実施を選択したときには光微分波形のデコンボリューションを行った上で光音響画像を生成し、ユーザがデコンボリューション処理の不実施を選択したときは光微分波形のデコンボリューションを行わずに光音響画像を生成してもよい。例えば、光微分波形のデコンボリューションを行ったときは、光音響信号を赤・黒の色に対応付けて表示し、デコンボリューションなしのときは、光音響信号を青・黒の色に対応付けて表示してもよい。

　また、デコンボリューションなしの場合の光音響画像を生成し、コンピュータがその光音響画像を解析することで、血管部分が２本に分かれているか否かを判定し、血管が２本に分かれていると判定されたときに、その血管部分のみを対象に光微分波形のデコンボリューション処理を行うようにしてもよい。その際、デコンボリューション処理を実施した血管部分の表示色を、他の未処理の血管部分の表示色とは異なる色とし、デコンボリューション処理が行われた血管と、他の未処理の血管とが容易に判別可能になるようにしてもよい。

　以上、本発明をその好適な実施形態に基づいて説明したが、本発明の光音響信号処理装置及び方法は、上記実施形態にのみ限定されるものではなく、上記実施形態の構成から種々の修正及び変更を施したものも、本発明の範囲に含まれる。

Claims

　被検体に対する光出射光により被検体内で生じた光音響信号をサンプリングするサンプリング手段と、
　前記サンプリングされた光音響信号から前記被検体への出射光の光強度の時間波形の微分波形である光微分波形をデコンボリューションした信号を生成する光微分波形逆畳込み手段とを備えたことを特徴とする光音響信号処理装置。
　前記被検体への出射光の光強度の時間波形の微分波形を取得する光微分波形取得手段を更に備えたことを特徴とする請求項１に記載の光音響信号処理装置。
　前記光微分波形逆畳込み手段が、
　前記サンプリングされた光音響信号をフーリエ変換する第１のフーリエ変換手段と、
　前記光微分波形を所定のサンプリングレートでサンプリングした信号をフーリエ変換する第２のフーリエ変換手段と、
　前記フーリエ変換された光微分波形の逆数を逆フィルタとして求める逆フィルタ演算手段と、
　前記フーリエ変換された光音響信号に前記逆フィルタを適用するフィルタ適用手段と、
　前記逆フィルタが適用された光音響信号をフーリエ逆変換するフーリエ逆変換手段と、
を有するものであることを特徴とする請求項１又は２に記載の光音響信号処理装置。
　前記光音響信号が第１のサンプリングレートでサンプリングされ、前記光微分波形が前記第１のサンプリングレートよりも高い第２のサンプリングレートでサンプリングされており、
　前記第１のサンプリングレートでサンプリングされた光音響信号を、前記第２のサンプリングレートでリサンプルするリサンプル手段を更に備え、
　前記第１のフーリエ変換手段が、前記リサンプル手段でリサンプルされた光音響信号をフーリエ変換するものであることを特徴とする請求項３に記載の光音響信号処理装置。
　前記第１のフーリエ変換手段と前記第２のフーリエ変換手段とが、おなじデータ点数でフーリエ変換を行うものであることを特徴とする請求項４に記載の光音響信号処理装置。
　前記光音響信号が第１のサンプリングレートでサンプリングされ、前記光微分波形が前記第１のサンプリングレートよりも高い第２のサンプリングレートでサンプリングされ、かつ、前記第１のフーリエ変換手段が第１のデータ点数でフーリエ変換を行い、前記第２のフーリエ変換手段が、前記第１のデータ点数よりも多い第２のデータ点数でフーリエ変換を行うものであり、
　前記フーリエ変換された光音響信号に対して、前記第１のデータ点数と前記第２のデータ点数との差の分だけ中央に０を付加するゼロパディングを行うゼロパディング手段を更に備え、
　前記フィルタ適用手段が、前記ゼロパディング手段でゼロパディングが施された信号に対して前記逆フィルタを適用するものであることを特徴とする請求項３に記載の光音響信号処理装置。
　前記光音響信号が第１のサンプリングレートでサンプリングされ、前記光微分波形が前記第１のサンプリングレートよりも高い第２のサンプリングレートでサンプリングされ、かつ、前記第１のフーリエ変換手段が第１のデータ点数でフーリエ変換を行い、前記第２のフーリエ変換手段が、前記第１のデータ点数よりも多い第２のデータ点数でフーリエ変換を行うものであり、
　前記フーリエ変換された光微分波形から、前記第１のデータ点数と前記第２のデータ点数の差の分だけ高周波成分サンプル点を除去する高周波成分サンプル点除去手段を更に備え、
　前記逆フィルタ演算手段が、前記光微分波形をフーリエ変換した信号から高周波成分サンプル点を除去した信号の逆数を逆フィルタとして求めるものであることを特徴とする請求項３に記載の光音響信号処理装置。
　前記第２のデータ点数は、前記第１のデータ点数に、前記第２のサンプリングレートと前記第１のサンプリングレートとの比を乗じたデータ点数と等しいか、又はそれよりも多いデータ点数であることを特徴とする請求項６又は７に記載の光音響信号処理装置。
　前記光音響信号のサンプリング間隔が、前記被検体への出射光のパルス時間幅よりも長いことを特徴とする請求項４から８何れかに記載の光音響信号処理装置。
　複数の検出素子で検出され、前記サンプリング手段でサンプリングされた光音響信号に基づいて前記光音響信号を再構成する光音響信号再構成手段を更に備え、
　前記光微分波形逆畳込み手段が、前記光音響信号再構成手段で再構成された光音響信号から前記光微分波形をデコンボリューションするものであることを特徴とする請求項１から９の何れかに記載の光音響信号処理装置。
　前記光微分波形をデコンボリューションした信号から、前記光音響信号を検出する検出器の受信角度依存特性の影響を除去するように、前記光微分波形をデコンボリューションした信号を補正する補正手段を更に備えることを特徴とする請求項１から１０何れかに記載の光音響信号処理装置。
　前記補正手段が、前記検出器の受信角度依存特性に代えて、又は加えて、前記光微分波形をデコンボリューションした信号から前記被検体における光の入射光分布の影響を除去するように、前記光微分波形をデコンボリューションした信号を補正するものであることを特徴とする請求項１１に記載の光音響信号処理装置。
　前記光微分波形逆畳込み手段が、前記デコンボリューション時に、更に、ノイズ増幅周波数帯をフィルタリングすることを特徴とする請求項１から１２何れかに記載の光音響信号処理装置。
　被検体に対する出射光が複数の波長の光を含み、前記サンプリング手段が各波長の光に対応した光音響信号をサンプリングし、前記光微分波形逆畳込み手段が各波長の光に対応した光音響信号から光微分波形をデコンボリューションした信号を生成するものであり、
　各波長の光に対応した、光音響信号から光微分波形をデコンボリューションした信号を演算処理する２波長データ演算手段を更に備えたことを特徴とする請求項１から１３何れかに記載の光音響信号処理装置。
　前記光微分波形をデコンボリューションした信号に基づいて光音響画像を生成する光音響画像生成手段を更に備えること特徴とする請求項１から１３何れかに記載の光音響信号処理装置。
　前記サンプリング手段が、前記被検体に対して送信された音響波に対する反射音響波を更にサンプリングするものであり、
　前記サンプリングされた反射音響波に基づいて反射音響波画像を生成する反射音響波画像生成手段と、
　前記光音響画像と前記反射音響波画像とを合成する画像合成手段とを更に備えることを特徴とする請求項１５に記載の光音響信号処理装置。
　前記画像合成手段が、前記光音響画像と反射音響波画像とを重畳することで画像合成を行うものであることを特徴とする請求項１６に記載の光音響信号処理装置。
　被検体に対する出射光に起因して被検体内で発生した光音響信号を検出するステップと、
　前記検出された光音響信号から、前記被検体への出射光の光強度の時間波形の微分波形である光微分波形をデコンボリューションするステップとを有する光音響信号処理方法。
　前記デコンボリューションするステップが、
　前記検出された光音響信号をフーリエ変換するステップと、
　前記光微分波形をフーリエ変換するステップと、
　前記フーリエ変換された光微分波形の逆数を逆フィルタとして求めるステップと、
　前記フーリエ変換された光音響信号に前記逆フィルタを適用するステップと、
　前記逆フィルタが適用された光音響信号をフーリエ逆変換するステップと、
を含むことを特徴とする請求項１８に記載の光音響信号処理方法。