(57)【要約】【課題】  光音響波の受信信号に基づいて取得される初期音圧分布を用いて、被検体情報を精度良く取得できる光音響装置を提供することを目的とする。【解決手段】  本明細書が開示する光音響装置は、光が被検体に照射されることにより発生する光音響波を受信して受信信号を出力する音響波受信部と、受信信号を用いて被検体内の初期音圧分布を取得し、初期音圧分布と前記被検体の背景吸収係数分布を含む補正項とを含む光の伝播を表す方程式を用いて被検体情報を取得する処理部とを有する。【選択図】  図３
2015134155.tif 000002
【特許請求の範囲】
【請求項１】
光が被検体に照射されることにより発生する光音響波を受信して受信信号を出力する音響波受信部と、
前記受信信号を用いて前記被検体内の初期音圧分布を取得し、
前記初期音圧分布と前記被検体の背景吸収係数分布を含む補正項とを含む光の伝播を表す方程式を用いて、被検体情報を取得する処理部とを有することを特徴とする光音響装置。
【請求項２】
前記処理部は、
前記光の伝播を表す方程式を用いて前記被検体内の光フルエンス分布を取得し、
前記初期音圧分布と前記光フルエンス分布とを用いて前記被検体情報を取得することを特徴とする請求項１に記載の光音響装置。
【請求項３】
前記音響波受信部は、
複数の照射状態のそれぞれで光が前記被検体に照射されることにより発生した光音響波を受信して前記複数の照射状態に対応する複数の前記受信信号を出力し、
前記処理部は、
複数の前記受信信号のそれぞれを用いて前記複数の照射状態に対応する複数の前記初期音圧分布を取得し、
前記複数の照射状態に対応する複数の前記光の伝播を表す方程式を用いて、前記被検体情報を取得することを特徴とする請求項１に記載の光音響装置。
【請求項４】
前記処理部は、
前記光の伝播を表す方程式を満たすことを拘束条件として、複数の照射状態に対応する前記被検体内の複数の光フルエンス分布と、前記被検体内の光吸収体の吸収係数分布と、前記被検体内のグルナイゼン係数分布と、前記複数の照射状態に対応する複数の前記初期音圧分布とに基づく関係式を含む目的関数が最小となるように前記背景吸収係数分布と前記光吸収体の吸収係数分布とを取得し、
前記背景吸収係数分布と前記光吸収体の吸収係数分布との和を取得することを特徴とする請求項３に記載の光音響装置。
【請求項５】
前記処理部は、
照射状態がｉであるときの位置ｒにおける前記初期音圧分布の値をδｐ０（ｉ）（ｒ）、照射状態がｉであるときの位置ｒにおける前記光フルエンス分布の値をφ（ｉ）（ｒ）、位置ｒにおける前記背景吸収係数分布の値をμａＢ（ｒ）、位置ｒにおける拡散係数分布の値をＤ（ｒ）、位置ｒにおける前記被検体のグルナイゼン係数分布の値をΓ（ｒ）、位置ｒにおける光のソース項をｓ（ｒ）、位置ｒにおける前記光吸収体の吸収係数分布をδμａ（ｒ）とすると、
式（１）で示される前記光の伝播を表す方程式を満たすことを拘束条件として、式（２）および式（３）で示される前記関係式ｆを含む目的関数が最小とする前記背景吸収係数分布と前記拡散係数分布と前記光吸収体の吸収係数分布とを取得することを特徴とする請求項４に記載の光音響装置。
【数１】
2015134155.tif 000024
【数２】
2015134155.tif 000025
【数３】
2015134155.tif 000026
【請求項６】
前記処理部は、
式（１）に基づき式（３）の光フルエンス分布を消去し、前記関係式を含む目的関数を最小とする前記背景吸収係数分布および前記光吸収体の吸収係数分布を取得することを特徴とする請求項５に記載の光音響装置。
【請求項７】
前記処理部は、
位置ｒにおける前記初期音圧分布の値をδｐ０（ｒ）、位置ｒにおける前記光フルエンス分布の値をφ（ｒ）、位置ｒにおける前記背景吸収係数分布の値をμａＢ（ｒ）、位置ｒにおける拡散係数分布の値をＤ（ｒ）、位置ｒにおける前記被検体のグルナイゼン係数分布の値をΓ（ｒ）、位置ｒにおける光のソース項をｓ（ｒ）とするときに、
次の式を満たす前記光フルエンス分布を取得することを特徴とする請求項２に記載の光音響装置。
【数４】
2015134155.tif 000027
【請求項８】
前記処理部は、
前記背景吸収係数分布として位置によらず一様な値の吸収係数を用いることを特徴とする請求項１から７のいずれか１項に記載の光音響装置。
【請求項９】
前記処理部は、
前記背景吸収係数分布として前記被検体の平均的な吸収係数を用いることを特徴とする請求項８に記載の光音響装置。
【請求項１０】
前記処理部は、
前記背景吸収係数分布として空間分布を有する吸収係数を用いることを特徴とする請求項１から７のいずれか１項に記載の光音響装置。
【請求項１１】
前記背景吸収係数分布を取得する吸収係数取得部を更に有することを特徴とする請求項１から１０のいずれか１項に記載の光音響装置。
【請求項１２】
光が被検体に照射されることにより発生する光音響波を受信して受信信号を出力する音響波受信部と、
前記受信信号を用いて前記被検体内の初期音圧分布を取得し、
前記初期音圧分布と初期音圧分布の補正とを含む光の伝播を表す方程式を用いて、被検体情報を取得する処理部とを有することを特徴とする光音響装置。
【請求項１３】
前記処理部は、
前記音響波受信部が測定することのできなかった光音響波に由来する初期音圧分布を前記初期音圧分布の補正として用いることを特徴とする請求項１２に記載の光音響装置。
【請求項１４】
光が被検体に照射されることにより発生する光音響波を受信して得られた受信信号を用いて前記被検体内の初期音圧分布を取得する工程と、
前記初期音圧分布と前記被検体の背景吸収係数分布を含む補正項とを含む光の伝播を表す方程式を用いて、被検体情報を取得する工程とを有することを特徴とする信号処理方法。
【請求項１５】
請求項１４に記載の信号処理方法をコンピュータに実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は光が被検体に照射されることにより発生する光音響波を受信して被検体情報を取得する光音響装置、信号処理方法、およびプログラムに関するものである。
【背景技術】
【０００２】
生体などの被検体に光が照射されると、被検体内で伝播・拡散した光が吸収され圧力が生じ、光音響波が発生する。近年、この光音響効果を用いた装置が盛んに研究、開発されている。特に、発生した光音響波をさまざまな方向から受信し、受信した光音響波を解析処理することで、被検体内部の被検体情報の空間分布を取得することができる光音響イメージング（Ｐｈｏｔｏ−Ａｃｏｕｓｔｉｃ  Ｉｍａｇｉｎｇ：ＰＡＩ）が注目されている。
【０００３】
ＰＡＩで得られる被検体情報の中でも期待されているのが吸収係数分布である。光の吸収は主に血管内のヘモグロビン（オキシヘモグロビン、デオキシヘモグロビン）由来であり、光の波長を変えて吸収係数分布を測定することで、各ヘモグロビンの濃度分布を求めることができる。また、各ヘモグロビンの濃度分布に基づいて、総ヘモグロビン中の酸化ヘモグロビンの割合、すなわち酸素飽和度分布を求めることができる。このような被検体情報を精度良く取得するために各波長における吸収係数分布を精度良く取得することが望まれている。
【０００４】
以下、吸収係数分布の取得方法を説明する。理解を助けるためにまずは実際に起きる現象を順番に解説し、その次に現象を逆にたどるように吸収係数分布の取得方法を説明する。上述のように、光音響効果は、被検体に光が照射され、光が拡散・伝播し、光を吸収して圧力が生じ、光音響波が発生するという現象であった。生体内の可視から近赤外の波長の光の伝播・拡散は非常に浅いところを除けば光拡散方程式でよく近似できることが知られており、定常状態の光拡散方程式は式（１）のように書くことができる。
【０００５】
【数１】
2015134155.tif 000003
ここで、Ｄは拡散係数分布、φは被検体内部の光フルエンス分布、μａは被検体内部の吸収係数分布（真の吸収係数分布）、ｓは光のソース項、ｒは位置である。また圧力と吸収係数と光フルエンスには式（２）のような関係がある。
【０００６】
【数２】
2015134155.tif 000004
ここで、ｐ０は真の初期音圧分布、Γはグルナイゼン係数分布である。この発生した真の初期音圧ｐ０の光音響波が被検体を伝播し、音響波受信器で測定される。
【０００７】
これを逆にたどるように吸収係数分布を取得する。まず、被検体中を伝播してきた光音響波の音圧の時間変化を測定し、非特許文献１に記載されたような方法でその測定結果から被検体に生じた音圧分布を取得する。これを初期音圧分布δｐ０と呼ぶこととする。
【０００８】
この初期音圧分布δｐ０を真の初期音圧分布ｐ０だと考え、式（２）に基づいて、初期音圧分布δｐ０を光フルエンス分布φおよびグルナイゼン係数分布Γで除することにより吸収係数分布μａを得ることができる。ところがこの手続きから理解されるように、初期音圧分布δｐ０から吸収係数分布μａを得るためには、光フルエンス分布φを求めることが必要である。そこで式（１）を用いて光フルエンス分布φを求める。例えば式（１）中のμａφを、式（２）に基づきｐ０／Γへと変形して式（３）を作る。
【０００９】
【数３】
2015134155.tif 000005
よって、初期音圧分布δｐ０を真の初期音圧分布ｐ０だと考える、つまり式（４）のようだと考えると、式（４）を光フルエンス分布φについて解き、求めることができる。
【００１０】
【数４】
2015134155.tif 000006
このように、初期音圧分布δｐ０を含む光伝搬方程式を解くことにより光フルエンス分布φを取得することができるのである。
【先行技術文献】
【非特許文献】
【００１１】
【非特許文献１】Ｃ．Ｌｉ，Ｌ．Ｖ．Ｗａｎｇ，“Ｐｈｏｔｏａｃｏｕｓｔｉｃ  ｔｏｍｏｇｒａｐｈｙ  ａｎｄ  ｓｅｎｓｉｎｇ  ｉｎ  ｂｉｏｍｅｄｉｃｉｎｅ”，Ｐｈｙｓ．Ｍｅｄ．Ｂｉｏｌ．５４（２００９）Ｒ５９．
【非特許文献２】Ｄ．Ｃｏｎｔｉｎｉ，ｅｔ．ａｌ，“Ｐｈｏｔｏｎ  ｍｉｇｒａｔｉｏｎ  ｔｈｒｏｕｇｈ  ａ  ｔｕｒｂｉｄ  ｓｌａｂ  ｄｅｓｃｒｉｂｅｄ  ｂｙ  ａ  ｍｏｄｅｌ  ｂａｓｅｄ  ｏｎ  ｄｉｆｆｕｓｉｏｎ  ａｐｐｒｏｘｉｍａｔｉｏｎ”．Ｉ．Ｔｈｅｏｒｙ”，Ａｐｐｌ．Ｐｈｙｓ．３６（１９９７）４５８７．
【発明の概要】
【発明が解決しようとする課題】
【００１２】
ところで、音響波受信器には受信可能な周波数帯域の制限があり、すべての周波数成分の光音響波を受信することが困難である。そのため音響波受信器では測定する周波数帯域を決める必要がある。例えばあるサイズの血管を観察対象とする場合、そのサイズの血管で発生する光音響波の周波数分布に合わせて音響波受信器の周波数帯域を決める。しかしながら、被検体から発生する光音響波は観察対象となる血管で主に発生する光音響波の周波数成分とは異なる周波数成分も含むため、上述のように周波数帯域を決めると、測定困難な周波数成分が存在する。そのため測定することができなかった周波数成分の分だけ、初期音圧分布は精度が低くなってしまう。
【００１３】
すなわち、初期音圧分布δｐ０は、真の初期音圧分布ｐ０とは異なる場合がある。そのため、初期音圧分布δｐ０を用いて、例えば式（４）に示す光拡散方程式を解いて取得される光フルエンス分布は精度が低い場合がある。
【００１４】
そこで本発明は、初期音圧分布δｐ０を用いて、被検体情報を精度良く取得できる光音響装置を提供することを目的とする。
【課題を解決するための手段】
【００１５】
本明細書が開示する光音響装置は、光が被検体に照射されることにより発生する光音響波を受信して受信信号を出力する音響波受信部と、受信信号を用いて被検体内の初期音圧分布を取得し、初期音圧分布と前記被検体の背景吸収係数分布を含む補正項とを含む光の伝播を表す方程式を用いて被検体情報を取得する処理部とを有する。
【発明の効果】
【００１６】
本発明に係る光音響装置によれば、光音響波の受信信号に基づいて得られる初期音圧分布を用いて、被検体情報を精度良く取得することができる。
【図面の簡単な説明】
【００１７】
【図１】実施形態１及び実施形態２に係る光音響装置の模式図
【図２】コンピュータとその他の構成との接続を示す図
【図３】実施形態１に係る光音響装置の作動フローを示す図
【図４】実施形態１において光を照射する様子を示す図
【図５】実施形態１における各種分布を示す図
【図６】実施形態２に係る光音響装置の作動フローを示す図
【図７】実施例１における数値実験に使用されたモデルを示す図
【図８】実施例１において観測結果を模擬計算する際に得られた光フルエンス分布を示す図
【図９】実施例１において得られた光フルエンス分布を示す図
【図１０】実施例１において得られた吸収係数分布を示す図
【図１１】比較例において得られた吸収係数分布を示す図
【図１２】実施例２における第二の照射状態を示す図
【図１３】実施例２におけるＧａｕｓｓ−Ｎｅｗｔｏｎ法の反復回数と目的関数の値を示す図
【図１４】実施例２において得られた吸収係数分布を示す図
【発明を実施するための形態】
【００１８】
前述したように、初期音圧分布δｐ０は、真の初期音圧分布ｐ０とは異なる場合がある。そのため式（４）を解くことによって得られる光フルエンス分布φは、真の初期音圧分布ｐ０とは異なる初期音圧分布δｐ０を用いて取得されているため、精度は低くなってしまう。
【００１９】
そのため式（２）に基づいて、初期音圧分布δｐ０を光フルエンス分布φおよびグルナイゼン係数分布Γで除することによって吸収係数分布を得ると、吸収係数分布の精度も低くなってしまう。
【００２０】
そこで本発明者は、上記課題に鑑みて、初期音圧分布δｐ０に初期音圧分布の補正ｐ０Ｂを加えて真の初期音圧分布ｐ０に近づけた上で、光フルエンス分布φを取得することを着想した。さらに本発明者は、初期音圧分布の補正ｐ０Ｂは音響波受信器が測定することのできなかった光音響波が原因であるとして光フルエンス分布φを取得することを見出した。すなわち、本発明者は初期音圧分布δｐ０と音響波受信器が測定することのできなかった光音響波に由来する初期音圧分布の補正ｐ０Ｂとを含む光伝搬を表す式を満たす光フルエンス分布φを取得することを見出した。
【００２１】
例えば、式（５）に示すように、初期音圧分布δｐ０と初期音圧分布の補正ｐ０Ｂとを含む光拡散方程式を解くことにより、精度の高い光フルエンス分布φを取得することができる。
【００２２】
【数５】
2015134155.tif 000007
【００２３】
既に述べたように、式（５）に含まれる初期音圧分布δｐ０と初期音圧分布の補正ｐ０Ｂの和が真の初期音圧分布ｐ０に近づくようにｐ０Ｂを設定することを考える。そうすれば、式（３）に示すような真の初期音圧分布を用いて光フルエンス分布φを取得することと近い効果を得ることができる。すなわち、式（５）により得られる光フルエンス分布φは、初期音圧分布の補正ｐ０Ｂを考慮せずに取得された光フルエンス分布、例えば式（４）により得られる光フルエンス分布と比べて真の値に近づく。
【００２４】
そして式（６）に基づいて、初期音圧分布δｐ０と初期音圧分布の補正ｐ０Ｂとの和を光フルエンス分布φおよびグルナイゼン係数分布Γで除することによって得られる吸収係数分布μａは、差分ｐ０Ｂを考慮しない場合と比べて真の値に近づく。なお式（６）は式（２）の初期音圧分布ｐ０を初期音圧分布δｐ０と初期音圧分布の補正ｐ０Ｂとの和に置き換えたものである。
【００２５】
【数６】
2015134155.tif 000008
また、このように取得された吸収係数分布μａに基づいて得られる他の被検体情報についても真の値に近づく。
【００２６】
本発明に係る光音響装置により得られる被検体情報としては、被検体の吸収係数分布や被検体を構成する物質の濃度分布などがある。物質の濃度とは、酸素飽和度、オキシヘモグロビン濃度、デオキシヘモグロビン濃度、および総ヘモグロビン濃度などである。総ヘモグロビン濃度とは、オキシヘモグロビン濃度およびデオキシヘモグロビン濃度の和である。
【００２７】
なお、以下、初期音圧分布の補正ｐ０Ｂとして、音響波受信器が測定することのできなかった光音響波に由来すると考える例を説明する。ただし本発明によれば、初期音圧分布の補正はこれ以外のあらゆる要因に起因して生じるものであってかまわない。例えば、光音響波の受信信号から初期音圧分布を取得する際に生じた精度の低下を補正するとしてもよい。
【００２８】
また、以下、本発明を光拡散方程式に適用する例を説明するが、本発明を適用可能な光伝搬を表す式はこれに限らない。真の初期音圧分布ｐ０を用いて光フルエンス分布を取得することのできる光伝搬を表す式である限り、本発明を適用することができる。例えば、この場合、光伝搬を表す式中の真の初期音圧分布ｐ０を、初期音圧分布δｐ０と初期音圧分布の補正ｐ０Ｂとの和に置き換えることにより、本発明を適用することができる。
【００２９】
［実施形態１］
実施形態１では、本実施形態に係る光音響装置が取得した初期音圧分布と、音響波受信器が測定できなかった光音響波に由来する初期音圧分布の補正とを用いて光フルエンス分布を取得する例を説明する。
【００３０】
本実施形態では、音響波受信器が測定できなかった光音響波にはそれに対応する吸収係数が存在すると考える。この吸収係数を背景吸収係数分布μａＢとする。すなわち、本実施形態においては、初期音圧分布の補正ｐ０Ｂは式（７）と書けることとする。
【００３１】
【数７】
2015134155.tif 000009
【００３２】
なお、背景吸収係数分布μａＢは、位置によらず一様な値であってもよいし、被検体内の複数の位置で異なる値であってもよい。例えば、被検体内の平均的な吸収係数の値を、被検体内の各位置の背景吸収係数分布の値とすることができる。または、被検体内を複数の領域に分割し、それぞれの領域内の平均的な吸収係数の値を各領域内の各位置の背景吸収係数分布の値とすることもできる。
【００３３】
図１に示す模式図を用いて本実施形態に係る光音響装置を説明する。本実施形態に係る光音響装置は、光源１１０、光学系１２０、音響波受信部としてのトランスデューサアレイ１３０、処理部としてのコンピュータ１５０、表示装置１６０、入力部１７０、吸収係数取得部１８０を有する。
【００３４】
以下、本実施形態に係る光音響装置の各構成について説明する。
【００３５】
（被検体１００及び光吸収体１０１）
これらは本発明の光音響装置の一部を構成するものではないが、以下に説明する。本発明の光音響装置は、人や動物の悪性腫瘍や血管疾患などの診断や化学治療の経過観察などを主な目的とする。よって、被検体１００としては生体、具体的には人体や動物の乳房や頸部、腹部などの診断の対象部位が想定される。
【００３６】
また、被検体１００の内部にある光吸収体１０１が観察対象となる。光吸収体１０１は、被検体１００の内部で相対的に吸収係数が高いものが好ましい。また、光吸収体１０１は、被検体１００よりも十分に小さいものが想定される。例えば、被検体１００が人体である場合、オキシヘモグロビン、デオキシヘモグロビン、メトヘモグロビン、あるいはそれらを含む血管や新生血管が光吸収体１０１の対象となる。その他、頸動脈壁のプラークなどもその対象となる。また、造影剤を用いる場合、造影剤が光吸収体１０１の対象となる。
【００３７】
（光源１１０）
光源は、効率的に光音響波を発生させるため、数ナノから数マイクロ秒オーダーのパルス光を発生可能なパルス光源である好ましい。具体的に光源は５ｎｓｅｃ以上５００ｎｓｅｃ以下のパルス幅の光を発生させることが好ましい。
【００３８】
また、光の波長は、被検体内部まで光が伝搬する波長であることが好ましい。具体的に、被検体が生体の場合は水による光の吸収を避けるために、光源は４００ｎｍ以上１６００ｎｍ以下の波長の光を発生させることが好ましい。また、波長依存性を有する光吸収体を測定するために、光源は単一波長の光のみを発生するものではなく、異なる波長を発生できるものであってもよい。
【００３９】
例えば、光源としてはレーザや発光ダイオードなどを用いることも可能である。レーザとしては、固体レーザ、ガスレーザ、色素レーザ、半導体レーザ、ＯＰＯ（Ｏｐｔｉｃａｌ  Ｐａｒａｍｅｔｒｉｃ  Ｏｓｃｉｌｌａｔｏｒｓ）を用いたレーザなど様々なレーザを使用することができる。
【００４０】
（光学系１２０）
光学系１２０は、光源から出射された光を所望の光分布形状に加工しながら被検体１００に導く部材である。光学系１２０は、レンズやミラーなどの光学部品や光ファイバなどの光導波路などにより構成される。例えば光学部品としては、光を反射するミラーや、光を集光したり拡大したり形状を変化させるレンズ、光を拡散させる拡散板などである。このような光学系１２０は、光源から発せられた光が被検体に所望の形状で照射されれば、どのようなものを用いてもかまわない。また強度については生体の安全性を確保した上で、ノイズに対するシグナルの割合を確保できるように選べばよい。
【００４１】
また被検体の被検体情報を得たい領域が広いため、一度にすべての領域を所望の強度で照らすことが出来ない場合には、複数回照射してスキャンすることなどのよって、所望の領域を十分な強度で照らすことができるようにしてもよい。
【００４２】
（トランスデューサアレイ１３０）
トランスデューサアレイに１３０に用いられる音響波受信器は、音響波を受信して電気信号に変換できるものであればどのようなトランスデューサを用いてもよい。トランスデューサアレイ１３０が観察対象から発生する光音響波の周波数成分をできる限り受信できるように、トランスデューサの受信帯域は観察対象の大きさに合わせて設計することが好ましい。
【００４３】
例えば、トランスデューサとしては、圧電現象、光の共振、静電容量の変化等を用いたものなどを用いることができる。
【００４４】
また、被検体の形状などの事前情報を用いずに、被検体情報の空間分布を取得するためには、複数の位置で音響波を受信する必要がある。そこで本実施形態のようにトランスデューサをアレイ状に配置し、複数の位置で同時受信してもよい。また、移動部がトランスデューサを移動させることにより、トランスデューサが複数の位置で音響波を受信してもよい。
【００４５】
（音響整合材１４０）
音響整合材１４０は本発明の光音響装置の一部を構成するものではないが説明する。音響整合材１４０は、被検体１００とトランスデューサアレイ１３０との間に設けられることにより、被検体１００とトランスデューサアレイ１３０との音響整合を図る部材である。音響整合材１４０としては、水、ならびに、水、オイル、およびアルコールなどから構成されるジェルなどが用いられる。
【００４６】
（コンピュータ１５０）
コンピュータ１５０には典型的にはワークステーションや大規模並列クラスタなどが用いられ、補正処理や画像再構成処理などがあらかじめプログラミングされたソフトウェアにより行われる。なお、ワークステーションで行うようなソフトウェア処理ではなく、ハードウェア処理により被検体情報を取得することもできる。また、本実施形態においてコンピュータ１５０が実行するそれぞれの処理を、それぞれ別の装置で行ってもよい。
【００４７】
コンピュータ１５０内の演算部１５２は、トランスデューサアレイ１３０から出力された電気信号に対して所定の処理を施すことができる。
【００４８】
また、制御部としての演算部１５２は、図２に示すようにバス２００を介して光音響装置を構成する各構成の作動を制御することができる。例えば、演算部１５２は、光源から発生する光の発光タイミングの制御し、光をトリガ信号として光音響波の受信のタイミングを制御することができる。
【００４９】
演算部１５２は、典型的にはＣＰＵ、ＧＰＵ、Ａ／Ｄ変換器などの素子や、ＦＰＧＡ、ＡＳＩＣなどの回路から構成される。なお、演算部１５２は、１つの素子や回路から構成されるだけではなく、複数の素子や回路から構成されていてもよい。また、演算部１５２が行う各処理をいずれの素子や回路が実行してもよい。
【００５０】
また、コンピュータ１５０内の記憶部１５１は、典型的にはＲＯＭ、ＲＡＭ、およびハードディスクなどの記憶媒体から構成される。なお、記憶部１５１は、１つの記憶媒体から構成されるだけでなく、複数の記憶媒体から構成されていてもよい。
【００５１】
また、コンピュータ１５０は、同時に複数の信号をパイプライン処理できるように構成されていることが好ましい。これにより、被検体情報を取得するまでの時間を短縮することができる。
【００５２】
なお、コンピュータ１５０が行う信号処理や光音響装置の作動制御を実行させるプログラムを記憶部１５１に保存しておくことができる。ただし、プログラムが保存される記憶部１５１は、非一時的な記録媒体である。
【００５３】
（表示装置１６０）
表示装置１６０は、コンピュータ１５０から出力された被検体情報の画像や数値情報などを表示する装置である。表示装置１６０には、液晶ディスプレイ、プラズマディスプレイ、有機ＥＬディスプレイ、およびＦＥＤなどあらゆる方式のディスプレイを採用することができる。なお、表示装置１６０は、本発明の光音響装置とは別に提供されていてもよい。
【００５４】
（入力部１７０）
入力部１７０は、コンピュータ１５０に情報を入力するために、ユーザーが情報を指定できるように構成された部材である。入力部１７０としては、キーボード、マウス、タッチパネル、ダイヤル、およびボタンなどを用いることができる。入力部１７０としてタッチパネルを採用する場合、表示装置１６０が入力部１７０を兼ねるタッチパネルであってもよい。
【００５５】
（吸収係数取得部１８０）
吸収係数取得部１８０は、被検体１００の背景吸収係数分布を取得する装置である。例えば、吸収係数取得部１８０としては、近赤外分光法に時間分解分光法を導入して被検体１００の背景吸収係数分布を取得する時間分解分光システムを採用することができる。その他、吸収係数取得部１８０としては、ＭＲＩ、ＣＴ、超音波診断装置などの種々の画像形成装置を採用することができる。このとき取得される背景吸収係数分布は、空間的に一様な値（例えば、被検体１００内のの平均な値）であってもよいし、空間分布を有するものであってもよい。
【００５６】
なお、吸収係数取得部１８０として時間分解分光システムなどを採用する場合、吸収係数取得部１８０が還元散乱係数分布（ｒｅｄｕｃｅｄ  ｓｃａｔｔｅｒｉｎｇ  ｃｏｅｆｆｉｃｉｅｎｔ）分布μＳ’や拡散係数分布Ｄを併せて取得してもよい。このとき取得される還元散乱係数分布や拡散係数分布は、空間的に一様な値（例えば、被検体１００内のの平均な値）であってもよいし、空間分布を有するものであってもよい。なお、還元散乱係数分布μｓ’および拡散係数分布Ｄは式（８）および式（９）で定義される係数である。ｇは散乱異方性、μｓは散乱係数を示す。
【００５７】
【数８】
2015134155.tif 000010
【００５８】
【数９】
2015134155.tif 000011
【００５９】
また、吸収係数取得部１８０としての画像形成装置で得られた被検体１００の内部の画像に基づいて、演算部１５２が各領域の背景吸収係数分布の値を推定してもよい。この場合、演算部１５２を吸収係数取得部１８０の一部の構成とすることができる。
【００６０】
次に、図３に示すフローを用いて本実施形態に係る光音響装置の作動を説明する。
【００６１】
（Ｓ１００：光音響波の受信信号を取得する工程）
まず光源１１０は光１２１を発生し、光１２１が光学系１２０を介して被検体１００に照射される。そして、光１２１が被検体１００の内部に位置する光吸収体１０１に吸収される。光１２１を吸収した光吸収体１０１が瞬間的に膨張することにより、光音響波１０２が発生する。
【００６２】
トランスデューサアレイ１３０は、光音響波１０２を受信し、時系列の受信信号に変換する。コンピュータ１５０は、トランスデューサアレイ１３０から出力された時系列の受信信号に増幅、Ａ／Ｄ変換の処理を施し、コンピュータ１５０内の記憶部１５１に格納する。本実施形態においてコンピュータ１５０は、光源１１０が光１２１を発生したタイミングに、トランスデューサアレイ１３０から出力された受信信号に対して上記処理を開始する。
【００６３】
なお、本工程で被検体１００に照射された光１２１の強度情報は記憶部１５１に格納され、後述するＳ４００で説明する照射光の境界条件として用いられる。
【００６４】
（Ｓ２００：光音響波の受信信号を用いて初期音圧分布を取得する工程）
コンピュータ１５０は、記憶部１５１に格納された時系列の受信信号を用いて非特許文献１に記載されたような再構成処理を行い、被検体１００内の初期音圧分布δｐ０を形成することができる。再構成アルゴリズムとしては、トモグラフィー技術で通常に用いられるタイムドメインあるいはフーリエドメインでの逆投影、Ｄ＆Ｓ（Ｄｅｌａｙ＆Ｓｕｍ）などが用いられる。また、再構成に多くの時間を要することが可能な場合は、逆問題などの分野で用いられる繰り返し手法などの画像再構成手法を用いることもできる。
【００６５】
なお、Ｓ１００においてフォーカスしたトランスデューサを順次移動して時系列の受信信号を取得した場合、再構成なしに高次元の初期音圧分布を形成することができる。この場合、コンピュータ１５０は再構成処理を行う必要はない。
【００６６】
（Ｓ３００：被検体の背景吸収係数分布を取得する工程）
本実施形態では、吸収係数取得部１８０としての時間分解分光システムが、背景吸収係数分布を取得する。また本実施形態において吸収係数取得部１８０は、還元散乱係数分布も取得する。このとき取得される背景吸収係数分布や還元散乱係数分布は、空間的に一様な値（例えば、被検体１００内のの平均な値）であってもよいし、空間分布を有するものであってもよい。
【００６７】
なお、本工程においてユーザーが背景吸収係数分布を選択してもよい。例えば、記憶部１５１に格納された複数の背景吸収係数分布の中からユーザーが入力部１７０を用いて所望の背景吸収係数分布を選択してもよい。
【００６８】
また、被検体１００が生体である場合、まずユーザーが入力部１７０を用いて被検体１００の年齢や性別などの情報を入力し、入力された情報がコンピュータ１５０に出力されてもよい。続いて、演算部１５２が記憶部１５１に格納された複数の背景吸収係数分布の統計値の中からユーザーにより入力された情報に対応する背景吸収係数分布の統計値を選択することができる。
【００６９】
また、ユーザーが画像形成装置で得られた被検体１００の内部の画像を確認し、入力部１７０を用いて各領域の背景吸収係数分布の値を入力してもよい。
【００７０】
また、拡散係数分布もしくは還元散乱係数分布についても、同様にユーザーが入力部１７０を用いて入力することができる。
【００７１】
また、後述するＳ４００において、被検体１００の背景吸収係数分布として予め設定された値を用いる場合、本工程を行う必要はない。
【００７２】
（Ｓ４００：初期音圧分布と背景吸収係数分布を含む補正項とを含む光伝搬を表す方程式を用いて光フルエンス分布を取得する工程）
演算部１５２は、Ｓ２００で取得した初期音圧分布と、Ｓ３００で取得した被検体１００の背景吸収係数分布を含む補正項とを含む光伝搬を表す式を用いて、被検体１００内の光フルエンス分布を取得する。この光伝搬を表す式について、以下詳細に説明する。
【００７３】
図４Ａは本実施形態に係る光１２１の照射状態を示す。ここでは、光１２１を被検体１００に照射し、光吸収体１０１ａおよび１０１ｂで強い光吸収が生じる場合を考える。図４Ｂは光１２１を当てた時の被検体１００内のＡ−Ａ’上の光フルエンス分布を模式的に示したものである。
【００７４】
図５Ａは、図４Ａに示す場合のときにＡ−Ａ’上で発生する光音響波の真の初期音圧分布ｐ０を模式的に示したものである。すなわち、光吸収体１０１ａおよび１０１ｂで発生する光音響波に対応する真の初期音圧分布は４１０および４２０である。同様に、被検体１００の全体で発生する低周波成分の光音響波に対応する４３０も、真の初期音圧分布である。
【００７５】
ところが、通常トランスデューサアレイ１３０には周波数帯域の制限があるため、発生した光音響波のすべての周波数成分を測定することは困難である。また、本実施形態に係るトランスデューサ１３０は、光吸収体１０１ａおよび１０１ｂで発生する光音響波の周波数成分を主に受信するように構成されている。そのため、このトランスデューサ１３０からの受信信号を用いて得られる初期音圧分布δｐ０は模式的には図５Ｂに示すような形となり、初期音圧４３０が再現されているとは言い難い。
【００７６】
図５Ａと図５Ｂとを比較して理解されるように、真の初期音圧分布ｐ０と初期音圧分布δｐ０とは異なる。ここで、真の初期音圧分布ｐ０と初期音圧分布δｐ０との差を図５Ｃに示す。初期音圧分布の補正ｐ０Ｂが図５Ｃに示す差と同じになるように設定できれば、初期音圧分布δｐ０と初期音圧分布の補正ｐ０Ｂとの和は真の初期音圧分布ｐ０になるだろう。
【００７７】
前述したように本実施形態では初期音圧分布の補正ｐ０Ｂは式（７）によって表わされるものとした。すなわち本発明者は、図５Ｃに示す真の初期音圧分布ｐ０と初期音圧分布δｐ０との差は背景吸収係数分布μａＢに由来すると考えた。
【００７８】
式（５）に示す光拡散方程式に式（７）を代入すると、式（１０）を得ることができる。
【００７９】
【数１０】
2015134155.tif 000012
ここで、式（１０）中のΓμａＢφが本実施形態における「背景吸収係数分布を含む補正項」に相当する。本実施形態において、初期音圧分布δｐ０はＳ２００で取得され、背景吸収係数分布μａＢおよび拡散係数分布ＤはＳ３００で取得されている。また、グルナイゼン係数分布Γは例えば血管を想定して０．２５など、被検体を想定してある一定の値を設定するか、あらかじめ空間分布を求めることができればその値を入れればよい。ユーザーが入力部１７０を用いてグルナイゼン係数分布Γを入力することとしてもよい。式（１０）においては光フルエンス分布φのみが変数となるため、演算部１５２は式（１０）を解いて、光フルエンス分布φを取得することができる。
【００８０】
このようにして得られた光フルエンス分布は真の値に近いと想定される。なぜならそれは初期音圧分布には含まれなかった背景吸収係数分布の寄与（ΓμａＢφ）によって補正されているからである。すなわち、初期音圧分布に加えて、背景吸収係数分布も用いることにより光フルエンス分布を精度良く取得することができる。なお、式（１０）を変形した式を用いたとしても、等価な変形であれば、変形後の式を満たすように光フルエンス分布を取得する形態も本発明に包含される。
【００８１】
（Ｓ５００：初期音圧分布と光フルエンス分布とを用いて被検体情報を取得する工程）
演算部１５２は、式（６）及び式（７）により、被検体情報としての被検体１００内の吸収係数分布μａを取得する。本工程では、Ｓ４００において取得された精度の高い光フルエンス分布φを用いて吸収係数分布μａを取得しているため、精度の高い吸収係数分布μａを取得することができる。
【００８２】
そして、演算部１５２は、以上のように得られた吸収係数分布の画像や数値情報などを表示装置１６０に表示させることができる。
【００８３】
以上、本実施形態によれば、音響波受信器が測定することのできない光音響波があった場合であっても、初期音圧分布に加えて背景吸収係数分布も用いて光フルエンス分布を取得することにより、被検体内の光フルエンス分布を精度良く取得することができる。
【００８４】
なお、Ｓ１００からＳ５００までの工程を異なる波長の光を用いて実行し、異なる波長の光に対応する吸収係数分布を取得してもよい。また、演算部１５２は、複数の波長のそれぞれに対応する複数の吸収係数分布を用いて、被検体情報として被検体１００を構成する物質の濃度分布を取得してもよい。本実施形態によれば、各波長の光に対応する吸収係数分布のそれぞれが精度良く取得されるため、これらを用いて取得される物質の濃度分布も精度良く取得することができる。
【００８５】
（実施形態２）
実施形態１の方法では、用いられる背景吸収係数分布の精度が十分でない場合は、精度の高い光フルエンス分布を取得することができず、吸収係数分布についても精度良く取得することができない。そこで、実施形態２では互いに異なる複数の照射状態のそれぞれで光を被検体に照射する。そして、複数の照射状態に対応する複数の初期音圧分布と、複数の照射状態に対応する複数の光フルエンス分布と、被検体内の光吸収体の吸収係数分布と、被検体内のグルナイゼン係数分布とに基づく関係式を含む目的関数を最小にする最適化問題を解くことにより背景吸収係数分布を取得する。
【００８６】
ここでいう光吸収体の吸収係数分布とは、式（１１）で示されるδμａのことを指し、音響波受信器が測定できた光音響波に対応する吸収係数を意味している。またその定義かもしくは式（６）から式（７）の辺々を引けば理解できるように、光吸収体の吸収係数分布δμａは吸収係数分布μａから背景吸収係数分布μａＢを引いたものとなる。
【００８７】
【数１１】
2015134155.tif 000013
【００８８】
本実施形態によれば、背景吸収係数分布の精度が十分でない場合や、背景吸収係数分布を取得することができない場合であっても、背景吸収係数分布を精度良く取得することができる。
【００８９】
また、本実施形態によれば、音響波受信器が測定することのできない光音響波があった場合であっても、吸収係数分布を精度良く取得することができる。
【００９０】
なお、本明細書において異なる照射状態とは、被検体に対する光の照射位置が異なることを指す。照射状態間の初期音圧分布が有意に異なる程度に光の照射位置を変えて測定する必要がある。
【００９１】
また、光フルエンス分布を取得するために用いる背景吸収係数分布以外のパラメータについても背景吸収係数分布と同様に、そのパラメータを繰り返し演算により精度良く求めることができる。なお、本実施形態では、背景吸収係数分布に加えて、拡散係数分布および光吸収体の吸収係数分布の値についても最適化問題を解くことにより求める例を説明する。
【００９２】
なお本実施形態に係る光音響装置の構成は実施形態１で示した構成と同一である。また光学系１２０とトランスデューサアレイ１３０は上述したように異なる照射状態の初期音圧分布を取得できるように動作させる必要がある。
【００９３】
以下、図６に示すフローを用いて本実施形態に係る光音響装置の作動を説明する。なお、本実施形態においても、コンピュータ１５０が図１に示す光音響装置を構成する各構成の作動を制御している。
【００９４】
（Ｓ２１００：複数の照射状態で初期音圧分布を取得する工程）
まず、第１の照射状態で実施形態１のＳ１００およびＳ２００の工程を行い、第１の初期音圧分布δｐ０（１）を取得する。
【００９５】
続いて、照射状態を変えて、第２の照射状態でＳ１００およびＳ２００の工程を行い、第２の初期音圧分布δｐ０（２）を取得する。本実施形態では説明の簡単のために２パターンの光の照射状態の例を示しているが、３パターン以上の光の照射状態を用いても構わない。
【００９６】
（Ｓ２２００：複数の照射状態に対応する複数の光伝搬を表す方程式を用いて、背景吸収係数分布および吸収係数分布を取得する工程）
続いて演算部１５２は、Ｓ２１００で取得した初期音圧分布と、被検体１００の背景吸収係数分布を含む補正項とを含む複数の照射状態に対応する複数の光伝搬を表す式を用いて、背景吸収係数分布および吸収係数分布を取得する。
【００９７】
本実施形態では上記光伝播を表す式として式（１２）を用いることとする。式（１２）は式（１０）を光の照射状態ごとに書きなおしただけである。また式（１２）に示す光拡散方程式を満たすことを拘束条件として、式（１３）及び式（１４）で示す関係式ｆを含む目的関数を最小とするような、背景吸収係数分布μａＢと光吸収体の吸収係数分布δμａを求めるという最適化問題を解くこととする。上記最適化問題を解いた結果求まる背景吸収係数分布μａＢと光吸収体の吸収係数分布δμａとを足すことにより、吸収係数分布μａが求まる。
【００９８】
すなわち、演算部１５２は、式（１２）に示す初期音圧分布δｐ０（ｉ）と背景吸収係数分布μａＢを含む補正項とを含む光伝搬を表す式を満たすことを拘束条件として、各照射状態のそれぞれにおける光フルエンス分布φ（ｉ）と、光吸収体の吸収係数分布δμａと、被検体内のグルナイゼン係数分布Γと、各照射状態のそれぞれにおける初期音圧分布δｐ０（ｉ）とに基づく式（１３）及び式（１４）に示す関係式ｆを含む目的関数が最小となるように背景吸収係数分布μａＢと光吸収体の吸収係数分布δμａとを取得する。
【００９９】
【数１２】
2015134155.tif 000014
【０１００】
【数１３】
2015134155.tif 000015
【０１０１】
【数１４】
2015134155.tif 000016
【０１０２】
なお、数値計算によって上記最適化問題を解く場合には、厳密に目的関数の最小値に達することは非常に困難であるので、必要な精度に応じて計算を終了させてもよい。この場合も通例、目的関数が最小となるまで計算したと言える。
【０１０３】
式（１２）は初期音圧分布δｐ０（ｉ）と背景吸収係数分布μａＢを含む補正項とを含む光伝搬を表す式の一形態に過ぎない。また、式（１３）及び式（１４）は各照射状態のそれぞれにおける光フルエンス分布φ（ｉ）と、光吸収体の吸収係数分布δμａと、被検体内のグルナイゼン係数分布Γと、各照射状態のそれぞれにおける初期音圧分布δｐ０（ｉ）との関係式の一形態に過ぎない。よって異なる光伝播を表す式、異なる関係式を採用しても構わない。
【０１０４】
また、演算部１５２は、吸収係数分布μａに基づいてその他の被検体情報を取得することができる。
【０１０５】
以上の本実施形態に係る方法によれば、十分な精度の背景吸収係数分布を得られない場合であっても、背景吸収係数分布を精度良く取得することができる。また、初期音圧分布と背景吸収係数分布とを用いて光フルエンス分布を取得するため、背景吸収係数分布に起因して初期音圧の値が低下してしまう場合であっても、被検体の吸収係数分布を精度良く取得することができる。
【０１０６】
（実施例１）
本実施例では、実施形態１に基づいた数値実験を行った例を説明する。
【０１０７】
被検体１００は５ｃｍ角の立方体の形状をしているものとした。また、波長７００ｎｍの光１２１を被検体１００の上部全面から入射するものとした。なお、光１２１の強度を１０ｍＪ／ｍ２、グルナイゼン係数分布の値を１とした。
【０１０８】
光吸収体１０１は血管を模して円環体をなしており、大半径を１．５ｃｍ、小半径を１ｍｍとした。円環体の光吸収体１０１の中心が立方体の被検体１００の中心と一致するように光吸収体１０１を配置した。また、円環体の輪の中心軸が立方体の被検体１００のある側面と直交するように、光吸収体１０１を配置した。光吸収体１０１はヘモグロビンを想定し、人体の血管に標準的に存在する程度の密度を想定した。また酸素飽和度は９０％を想定した。そのため、波長７００ｎｍの光に対する光吸収体１０１の吸収係数が約０．２４ｍｍ−１となった。また、背景吸収係数分布としては位置によらず一様な値の吸収係数を設定し、光吸収体１０１の吸収係数の１％の大きさである０．００２４ｍｍ−１とした。また、還元散乱係数についても位置によらず一様に０．９６ｍｍ−１とした。拡散係数分布は式（９）から求めた。
【０１０９】
図７に被検体の中心を通り、上下に切った、ある断面での吸収係数の空間分布の様子を示す。なお、図７、図１０及び図１１における吸収係数の単位はｍ−１である。図７における明暗は吸収係数の大きさを表す。図７の白い円環が比較的に吸収係数の高い光吸収体１０１である。また、図７の黒い部分が比較的吸収係数の低い光吸収体１０１の背景構造であり、音響波受信器では受信が困難である低周波成分の光音響波を発生する。なお光１２１は紙面上部から入射する。
【０１１０】
なお数値計算には離散化が必要だが、被検体１００は直交等間隔格子で離散化され、各軸の格子数は１００個とした。
【０１１１】
今回は数値実験のため観測結果も模擬計算を行う。この計算方法にも任意性があるが、本実施例では以下の方法を用いた。まずは光の伝播だが式（１）を用いて光フルエンスの空間分布を求めた。光フルエンスの境界条件は非特許文献２で提案された方法を用いた。光の入射は被検体上部全面の境界条件のソース項に値を設定する方法で入射を表現することとし、被検体の屈折率は１．４とした。式（１）の離散化方法には中心差分を用い、連立一次方程式は共役勾配法を用いて、残差の２乗ノルムが１０−７以下となるまで解いた。求めた光フルエンスの空間分布の様子を図８に示す。上部より入射した光が拡散、吸収されながら下部へと伝わる様子が見受けられるだろう。
【０１１２】
次は光音響波の受信信号を取得する工程Ｓ１００が実施される。つまり真の初期音圧分布ｐ０は式（２）で求めることができ、これが観測装置まで伝播し、観測される。しかし既に述べたように、観測装置の測定可能な周波数帯域には制限があり、測定されるのは一部にとどまる。本実施例ではそれを再現するために、音響波受信器が測定することのできなかった光音響波に由来する初期音圧分布を式（７）の右辺、即ちグルナイゼン係数分布と背景吸収係数分布と光フルエンスとの積とした。そして、工程Ｓ１００で観測される信号はｐ０と音響波受信器が測定することのできなかった光音響波に由来する初期音圧分布との差に相当する光音響波が伝播して、信号として観測されるものとした。
【０１１３】
次は光音響波の受信信号を用いて初期音圧分布を取得する工程Ｓ２００が実施される。非特許文献１に記載されたような再構成処理を行い、被検体１００内の初期音圧分布δｐ０を取得することができる。本実施例では理想的に再構成が行われたことを想定し、Ｓ１００で用いたｐ０と音響波受信器が測定することのできなかった光音響波に由来する初期音圧分布との差が初期音圧分布δｐ０として求まったとした。
【０１１４】
次は被検体の背景吸収係数分布を取得する工程Ｓ３００が実施される。ここでは近赤外分光法に時間分解分光法を導入した時間分解分光システムが理想的に動作したと想定して、位置によらず一様な値である背景吸収係数分布μａＢと     還元散乱係数分布μＳ’が正しく求まったとした。またそれらに基づいて、位置によらず一様な値である拡散係数分布Ｄが得られた。
【０１１５】
次は初期音圧分布と背景吸収係数分布とを用いて光フルエンス分布を取得する工程Ｓ４００が実施される。式（１０）で示される方程式をφについて解けばよい。光フルエンスの境界条件は非特許文献２で提案された方法を用いた。光の入射は被検体上部全面の境界条件のソース項に値を設定する方法で入射を表現することとし、被検体の屈折率は１．４とした。式（１０）の離散化方法には中心差分を用い、連立一次方程式は共役勾配法を用いて、残差の２乗ノルムが１０−７以下となるまで解いた。なお、式（１０）を解く際の離散化の方法や連立一次方程式の解法には任意性があり、公知のあらゆる方法により解くことができる。すなわち、本実施例の方法は一例に過ぎない。また、連立一次方程式を解く際の残差の計算を終了する閾値については、ユーザーが入力部１７０を用いて適宜設定することができる。図９に本工程で求まった光フルエンス分布φを示した。
【０１１６】
次は初期音圧分布と光フルエンス分布とを用いて被検体情報を取得する工程Ｓ５００が実施される。つまり式（６）及び式（７）により、被検体情報としての被検体１００内の吸収係数分布μａを取得する。図１０に本実施例で算出された吸収係数分布μａを示す。
【０１１７】
一方、比較例として、初期音圧分布δｐ０を用いて式（４）の光拡散方程式により光フルエンス分布φを算出し、さらに初期音圧分布δｐ０を光フルエンス分布φおよびグルナイゼン係数分布Γで除することにより得た吸収係数分布を図１１に示す。
【０１１８】
図７のように光吸収体１０１内の真の吸収係数の値は均一である。それに対して図１０の本実施例で求めた光吸収体１０１内の吸収係数は均一に分布している。一方、図１１の比較例で求めた光吸収体１０１の吸収係数は下部の吸収係数が低くなっており、不均一な吸収係数分布となっていることが理解される。また図１０に示す光吸収体１０１の吸収係数のほうが、図１１に示す光吸収体１０１の吸収係数と比べて真の値（０．２４ｍｍ−１）に近いことが理解される。すなわち、本実施例の方法により得られる吸収係数分布の方が、比較例の方法により得られる吸収係数分布よりも精度が高いことが理解される。
【０１１９】
本実施例によれば、音響波受信器が測定することのできない光音響波があった場合であっても、精度の高い吸収係数分布を取得することができる。
【０１２０】
（実施例２）
本実施例では、実施形態２に基づいた数値実験を行った例を説明する。対象となる被検体１００及び光吸収体１０１は実施例１と同様のものを用いる。
【０１２１】
本実施例では実施例１とは異なり、二つの照射状態を用いる。その様子を図７と図１２に示す。本実施例では、図７に示したように被検体１００の上部（紙面上部）全面に光１２１を照射する第１の照射状態と、図１２に示したように被検体１００の下部（紙面下部）全面に光１２１を照射する第２の照射状態との２つの照射状態を用いた。
【０１２２】
今回は数値実験のため観測結果も模擬計算を行う。これは実施例１と同様の方法を実行すればよいが、照射状態が２種類あるので、それぞれで実行し、信号として観測されるものとした。
【０１２３】
次は複数の照射状態で初期音圧分布を取得する工程Ｓ２１００が実施される。これも実施例１と同様の方法を実行し、各照射状態における初期音圧分布δｐ０（１）およびδｐ０（２）を取得する。
【０１２４】
次は複数の照射状態に対応する複数の光伝搬を表す方程式を用いて、背景吸収係数分布および吸収係数分布を取得する工程Ｓ２２００が実施される。ここでは式（１２）を満たすことを拘束条件として、式（１３）および式（１４）で示される関係式ｆを含む目的関数を最小とするような背景吸収係数分布μａＢ、拡散係数分布Ｄ、および光吸収体の吸収係数分布δμａを求めた。そして背景吸収係数分布μａＢと光吸収体１０１の吸収係数分布δμａとの和をとることで、真の吸収係数分布μａを求めた。
【０１２５】
なお本実施例では背景吸収係数分布μａＢおよび拡散係数分布Ｄは空間的に一様な値とした。背景吸収係数分布μａＢも拡散係数分布Ｄも空間的に分布を持つとしても構わないが、空間分布を表現する変数の数が多くなれば多くなるほど、精度は向上するものの解く手法は複雑になり、計算時間も膨大になるので、なるべく少なく抑えることが望ましい。
【０１２６】
本実施例では目的関数を最小化する方法にはＧａｕｓｓ−Ｎｅｗｔｏｎ法を用いた。また、拘束条件付きの最適化問題のため、式（１２）に示す拘束条件を用いて各照射状態における光フルエンス分布φ（１）およびφ（２）を消去した上で目的関数を最小とする変数を求めるという方法を用いた。以下で詳述する。
【０１２７】
Ｇａｕｓｓ−Ｎｅｗｔｏｎ法自体は公知であるため詳述しないが、その実施のためには目的関数と目的関数の勾配と目的関数の近似ヘッセ行列が必要である。そしてそれらは式（１４）に示す残差ｈ（ｉ）と残差のヤコビアン∂ｈ（ｉ）／∂ξがあれば求めることができる。ここでξは式（１５）で定義されるように表記の便利のために変数をまとめたものであり、ρによってどの変数を指すかを示している。
【０１２８】
【数１５】
2015134155.tif 000017
【０１２９】
残差ｈ（ｉ）は式（１４）で示されるが、これを求めるためにはφ（ｉ）が必要である。そのためには式（１２）を解いてφ（ｉ）を求めればよい。また残差のヤコビアン∂ｈ（ｉ）／∂ξは式（１６）で示されるが、これを求めるためには∂φ（ｉ）／∂ξが必要である。
【０１３０】
【数１６】
2015134155.tif 000018
そこで∂φ（ｉ）／∂ξを求めるために式（１２）を変数で微分する。すると式（１７）、式（１８）、式（１９）となる。
【０１３１】
【数１７】
2015134155.tif 000019
【０１３２】
【数１８】
2015134155.tif 000020
【０１３３】
【数１９】
2015134155.tif 000021
よって、式（１７）、式（１８）、式（１９）の方程式を解くことにより∂φ（ｉ）／∂ξを求め、式（１６）に代入し、残差のヤコビアン∂ｈ（ｉ）／∂ξを求めればよい。なお、式（１７）は簡単に解けて式（２０）となる。
【０１３４】
【数２０】
2015134155.tif 000022
このことにより、式（１７）を解くために必要な膨大な計算時間は必要なくなり、また目的関数の近似ヘッセ行列は疎行列となり、メモリの制限をさほど気にすることなく大規模な問題に対してもＧａｕｓｓ−Ｎｅｗｔｏｎ法によって解くことができる。これは式（１２）から変数をなるべく取り除くように変形することを本発明者が着想した結果である。
【０１３５】
なお背景吸収係数分布μａＢおよび拡散係数分布Ｄが、空間的に一様ではなく分布を持つ場合も同様の手続きで式変形を行えば∂φ（ｉ）／∂ξを求めることができる。但し背景吸収係数分布μａＢの空間分布をｍ個の変数で取り扱えば式（１８）に相当する連立一次方程式をｍ回解かねば∂φ（ｉ）／∂ξを求めることができない。また拡散係数分布Ｄの空間分布をｍ個の変数で取り扱えば式（１９）に相当する連立一次方程式をｍ回解かねば∂φ（ｉ）／∂ξを求めることができない。そのため空間分布を表現する変数の数が多くなれば多くなるほど、計算時間は膨大となるので、なるべく少なく抑えることが望ましい。
【０１３６】
次に最適化問題における変数の初期値について述べる。背景吸収係数分布μａＢおよび拡散係数分布Ｄの初期値としては被検体で想定されうる値を使えば良い。背景吸収係数分布μａＢおよび拡散係数分布Ｄはこの問題を解くことで求めることができるので、正確である必要はないが、おおよそ想定されうる値を使う方が好ましい。例えば実施形態１における被検体の背景吸収係数分布を取得する工程Ｓ３００を実施し、求めた値を用いてもよい。なお本実施例では想定されうる値が、それぞれ正しい値の２倍と０．５倍の値であったという例を示す。また、光吸収体１０１の吸収係数分布δμａの初期値は式（１２）および式（２１）に基づいて求めた。
【０１３７】
【数２１】
2015134155.tif 000023
すなわち、初期値として与えた背景吸収係数分布μａＢおよび拡散係数分布Ｄを用いて、式（１２）をφ（ｉ）について解き、式（２１）に基づいて光吸収体１０１の吸収係数分布δμａの初期値を求めた。なお本実施例の場合、照射状態は２種類であったので、ｉに関する和はｉ＝１とｉ＝２との和を取った。
【０１３８】
なお、式（１２）の離散化方法は中心差分を用いた。また、式（１９）、式（２０）の離散化方法もそれに準じた。式（１２）、式（１８）、式（１９）の連立一次方程式の解法には共役勾配法を用いて、残差の２乗ノルムが１０−７以下となるまで解いた。目的関数の降下方向を定める連立一次方程式の解法も共役勾配法を用いて、残差の２乗ノルムが１０−７以下となるまで解いた。直線探索は２分割法を用い目的関数が減少することを保証するようにした。なお、離散化の方法、連立一次方程式の解法、直線探索の方法には任意性があり、公知のあらゆる方法により解くことができる。すなわち、本実施例の方法は一例に過ぎない。また式（１２）で示される目的関数には最適化問題に対して用いられる様々な手法を適用することができる。例えば、Ｔｉｋｈｏｎｏｖ正則化やＴｏｔａｌ  Ｖａｒｉａｔｉｏｎ項を加えること、解の範囲制限を加えることなどである。
【０１３９】
図１３に横軸をＧａｕｓｓ−Ｎｅｗｔｏｎ法の反復回数、縦軸を目的関数の値としたときの本実施例におけるグラフを示す。このようにおおよそ６回程度で目的関数の値は十分に下がっている。
【０１４０】
図１４に本実施例で求めた真の吸収係数分布μａを示す。なお図１４における吸収係数の単位はｍ−１である。図１４によれば、求められた真の吸収係数分布μａは、背景吸収係数分布（０．００２４ｍｍ−１）と光吸収体１０１の吸収係数（０．２４ｍｍ−１）と近い値になっていることが理解される。
【０１４１】
本実施例によれば、本実施例のように背景吸収係数分布の初期値の精度が十分でない場合であっても、被検体の吸収係数分布を精度良く取得することができる。
【０１４２】
以上、特定の実施形態を説明したが、本発明はこれらの実施形態に限らず、特許請求の範囲を逸脱しない限りにおいて、種々の変形例、応用例も包含するものである。
【符号の説明】
【０１４３】
１００  被検体
１１０  光源
１３０  トランスデューサアレイ
１５０  コンピュータ
【図１】
2015134155.tif 000028
【図２】
2015134155.tif 000029
【図３】
2015134155.tif 000030
【図４】
2015134155.tif 000031
【図５】
2015134155.tif 000032
【図６】
2015134155.tif 000033
【図７】
2015134155.tif 000034
【図８】
2015134155.tif 000035
【図９】
2015134155.tif 000036
【図１０】
2015134155.tif 000037
【図１１】
2015134155.tif 000038
【図１２】
2015134155.tif 000039
【図１３】
2015134155.tif 000040
【図１４】
2015134155.tif 000041