実効散乱係数が未知の三次元光散乱体に対して、散乱係数μsを測定するステップと、算出した散乱係数μsに基づき、計算機シミュレーションを用いて散乱光の伝搬における非等方散乱因子gを0~1まで変化させて、散乱光からの出力ビーム径との相関関係を算出するステップと、三次元光散乱体にビーム光を照射した際の散乱光からの出力ビーム径の実測値と、前記相関関係とから、非等方散乱因子gを決定するステップと、を備えた、三次元光散乱体の実効散乱係数の算定方法。実効散乱係数が未知の三次元光散乱体に対して、散乱係数μsを測定するステップと、算出した散乱係数μsに基づき、計算機シミュレーションを用いて散乱光の伝搬における非等方散乱因子gを0~1まで変化させて、散乱光の透過/反射出力パワー比との相関関係を算出するステップと、三次元光散乱体にビーム光を照射した際の散乱光の透過/反射出力パワー比の実測値と、前記相関関係とから、非等方散乱因子gを同定するステップと、を備えた、三次元光散乱体の実効散乱係数の同定方法。散乱係数を測定するステップは、光コヒーレンストモグラフィーにより得られる振幅反射率の深さ依存性を用いて実測するものである、請求項1又は2に記載の三次元光散乱体の実効散乱係数の算定方法。前記計算機シミュレーションは、モンテカルロシミュレーションであり、散乱により変わる角度は非等方散乱因子gに依存するとされ、ある散乱から次の散乱までに進む距離は散乱係数μsに依存するとした、請求項1又は2に記載の三次元光散乱体の実効散乱係数の算定方法。非等方散乱因子を決定するステップにおいて、照射するビーム光の波長は、三次元散乱体における散乱体の粒径よりも小さくする、請求項1又は2に記載の三次元光散乱体の実効散乱係数の算定方法。本発明者は、フェムト秒レーザー加工を利用して、散乱係数制御が可能な人工的な三次元光散乱体に関する作製技術の研究を行っている。具体的には、ポリマー材料中にランダムにμmオーダーの空孔をあけ、その空孔密度により散乱係数が制御できる三次元光散乱体を研究している。また、SiO2微粒子をPMMAポリメタクリル酸メチルに拡散させる散乱体の作製技術の研究も行っている。これらの人工散乱体において、実験的に実効散乱係数を求めることが必要である。従来から、実効散乱係数を求める1つの方法として、光コヒーレンストモグラフィーを用いる方法が知られているが、これは深い領域での振幅反射率をもとに決定するため、ノイズに弱い特性を持つという欠点がある。上記状況に鑑みて、本発明は、実効散乱係数が未知の三次元光散乱体に対して、実効散乱係数を精度よく算定する方法を提供することを目的とする。上記状況に鑑みて、本発明の三次元光散乱体の実効散乱係数の算定方法は、実効散乱係数が未知の三次元光散乱体に対して下記a1~a3のステップを行う。a1散乱係数μsを測定するステップa2算出した散乱係数μsに基づき、計算機シミュレーションを用いて散乱光の伝搬における非等方散乱因子gを0~1まで変化させて、散乱光からの出力ビーム径との相関関係を算出するステップa3三次元光散乱体にビーム光を照射した際の散乱光からの出力ビーム径の実測値と、上記の相関関係とから、非等方散乱因子gを決定するステップかかるステップを備える方法によれば、散乱光のビーム径を指標として非等方散乱因子gを決定できる。光コヒーレンストモグラフィーなど既知の方法により散乱係数μsを測定し、得られた散乱係数μsを用いて、非等方散乱因子gを0から1まで変化させてモンテカルロシミュレーションを行い、散乱光の出力ビーム径を求める。例えば、非等方散乱因子gと出力ビーム径の相関グラフもしくは参照表を得る。ここで、相関グラフは、横軸が非等方散乱因子g0~1で、縦軸が出力ビーム径を示すグラフである。また、参照表は、非等方散乱因子g0~1をパラメータとする出力ビーム径の対応表である。そして、実際に散乱光にビームを照射して得られる出力ビーム径をもとに、相関グラフまたは参照表と照らし合わせて、非等方散乱因子gを決定する。本発明の三次元光散乱体の実効散乱係数の別の算定方法は、実効散乱係数が未知の三次元光散乱体に対して下記b1~b3のステップを行う。b1散乱係数μsを測定するステップb2算出した散乱係数μsに基づき、計算機シミュレーションを用いて散乱光の伝搬における非等方散乱因子gを0~1まで変化させて、散乱光の透過/反射出力パワー比との相関関係を算出するステップb3三次元光散乱体にビーム光を照射した際の散乱光の透過/反射出力パワー比の実測値と、上記の相関関係とから、非等方散乱因子gを同定するステップかかるステップを備える方法によれば、散乱光の透過/反射出力パワー比を指標として非等方散乱因子gを決定できる。光コヒーレンストモグラフィーなど既知の方法により散乱係数μsを測定し、得られた散乱係数μsを用いて、非等方散乱因子gを0から1まで変化させてモンテカルロシミュレーションを行い、散乱光の透過/反射出力パワー比を求める。例えば、非等方散乱因子gと透過/反射出力パワー比の相関グラフもしくは参照表を得る。ここで、相関グラフは、横軸が非等方散乱因子g0~1で、縦軸が透過/反射出力パワー比を示すグラフである。また、参照表は、非等方散乱因子g0~1をパラメータとする透過/反射出力パワー比の対応表である。そして、実際に散乱光にビームを照射して得られる透過/反射出力パワー比をもとに、相関グラフまたは参照表と照らし合わせて、非等方散乱因子gを決定する。ここで、上記のa1またはb1の散乱係数を測定するステップは、光コヒーレンストモグラフィーにより得られる振幅反射率の深さ依存性を用いて実測するものである。従来から用いられている光コヒーレンストモグラフィーでは、深い領域での振幅反射率をもとに決定するため、ノイズに弱い特性を持つという欠点があるが、本発明では、弱散乱領域での散乱係数を光コヒーレンストモグラフィーで測定し、透過または反射ビーム径または、透過/反射パワー比を用いて非等方散乱因子を決定するので、ノイズ特性が改善されることになる。また、上記のa2またはb2の計算機シミュレーションは、モンテカルロシミュレーションであり、散乱により変わる角度は非等方散乱因子に依存し、ある散乱から次の散乱までに進む距離は散乱係数μsに依存するものである。また、上記のa3またはb3の非等方散乱因子を決定するステップにおいて、照射するビーム光の波長は、三次元散乱体における散乱体の粒径よりも小さくする。散乱体の粒径よりもビーム光の波長が大きい場合、Mie散乱理論値との誤差が大きくなる。本発明によれば、実効散乱係数が未知の三次元光散乱体に対して、実効散乱係数を精度よく算定できるといった効果を有する。以下、本発明の実施形態について、図面を参照しながら詳細に説明していく。なお、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更及び変形が可能である。図1に、実施例1の三次元光散乱体の実効散乱係数の算定方法の処理フローを示す。図1に示すフローは、本発明の三次元光散乱体の実効散乱係数の算定方法の一実施形態であり、散乱光のビーム径を指標として非等方散乱因子gを決定するものである。実施例1の算定方法は、実効散乱係数が未知の三次元光散乱体に対して、散乱係数μsを測定するステップと、算出した散乱係数μsに基づき、計算機シミュレーションを用いて散乱光の伝搬における非等方散乱因子gを0~1まで変化させて、散乱光からの出力ビーム径との相関関係を算出するステップと、三次元光散乱体にビーム光を照射した際の散乱光からの出力ビーム径の実測値と、上記相関関係とから、非等方散乱因子gを決定するステップを備えるものである。本実施例1では、図1のフローに従う三次元光散乱体の実効散乱係数の算定方法について具体的に詳述する。上述したように、実効散乱係数は散乱係数と非等方散乱因子から決定される。先ず、実効散乱係数の決定方法を示す。次に、Mie散乱理論に基づく散乱係数および非等方散乱因子の計算方法について説明する。散乱係数を実験的に求める方法として、スペクトル領域光コヒーレンストモグラフィー以下、SD-OCTを用いる。そして、出力ビーム光径から非等方散乱因子を推定する。推定方法としては、モンテカルロシミュレーションを用いる。まず、実効散乱係数μ ́sは、散乱係数μsと非等方散乱因子gから下記の式1で与えられ、等価実効散乱係数ともいう。ここで、散乱係数μsは、媒質に入射した光が散乱するまでに進む距離である光浸達長Lの平均値の逆数で与えられ、下記の式2のようになる。また、非等方散乱因子gは、散乱の異方性を表す量であり、散乱方向の余弦平均として与えられる。gは完全後方散乱である-1と完全前方散乱である1の間の値をとり、g=0とき等方散乱を表す。光コヒーレンストモグラフィー以下、OCTは、マイケルソン型光学干渉計を応用した技術であり、低コヒーレンス光を光源として、非常に高い分解能で物体の内部構造を計測する方法である。分光技術と逆フーリエ変換を用いることにより、計算によって干渉信号を求めるFD-OCTフーリエドメインOCTは、参照光路を機械的に掃引することにより干渉信号を直接検出するTD-OCTタイムドメインOCTに対して、可動部がないため高速に計測が可能である。本実施例では、FD-OCTに分類されるSD-OCTスペクトルドメインOCTを用いるケースについて説明する。図2に、スペクトルドメインOCTの概略構成図を示す。図2に示すように、低コヒーレンス光を光源として光干渉計を構成する。光源からの光はビームスプリッタによって2分割され、一方は物体側、他方は参照鏡側に向かう。ここで対象となる物体は層構造を持つとし、その物体を干渉計の腕の一方に置く。その腕に直交する方向に反射鏡を置き、双方の腕の長さがほぼ等しくなるようにそれぞれを配置する。そして、物体にビームスプリッタ通過後の光を照射し、物体表面や内部構造の境界から反射もしくは散乱してきた光以下、物体光がビームスプリッタに戻る。そして物体光と参照鏡によって反射された光以下、参照光は、ビームスプリッタを通して重ね合わされて干渉する。その干渉光を回折格子で分光して得られたスペクトルを、レンズを通すことで波長ごとに分離し、光検出器を用いてスペクトル強度として検出する。このスペクトルでの干渉により生じた信号を干渉スペクトルと呼ぶことにする。ここで参照光と物体光の間に遅延時間が存在する場合、干渉スペクトルは遅延時間の逆数に比例した周期で変調される。検出している光スペクトルは波長成分であり、光周波数に換算して扱うことが可能であるので、干渉スペクトルに対して逆フーリエ変換を実行することにより、干渉スペクトル上の変調成分は時間軸上に変換され、参照光のパルス信号を基準とした遅延パルス信号として、反射面の振幅反射率や遅延時間が求められる。この遅延時間に対して光速を用いて距離に換算することにより、反射面の深さ位置を求めることができる。次に、モンテカルロシミュレーションについて説明する。本実施例では、出力光ビーム径から非等方散乱因子を推定するアプローチとして、モンテカルロシミュレーションを用いる。モンテカルロシミュレーションでは、光を光子に分解し、その1つ1つの散乱行程を確率に基づいて決定する。モンテカルロシミュレーションの計算に必要な情報は、散乱係数μs、非等方散乱因子g、吸収係数μaである。モンテカルロシミュレーションとその散乱行程の概略図を図3に示す。ここで、光子が一度散乱してから次の散乱までに進む距離Lは散乱係数μsと一様乱数R1を用いて下記の式3で表される。また、生体内では散乱角θは、下記の式4のような、非等方散乱因子gを用いたHenyey-Greenstein位相関数で表されるのが一般的である。モンテカルロシミュレーションにおいて、散乱角θは統計学的に位相関数pθに従わなければならない。そのために、等確率に0から1までの値をとる累積位相関数fθを下記の式5のように設定する。ここで、散乱角θは、fθと一様乱数R2を用いて、下記の式6のように表される。図4にモンテカルロシミュレーションのフローチャートを示す。まず、入射位置から光子が入射される。次に、入射された光子の移動距離を上記の式3から計算する。そして、計算された移動距離と進行角をもとに光子の新しい位置が決定される。その後、散乱角を上記の式6から計算し、新しい進行角とする。光子が媒体外に射出されると、射出位置を記録し、次の光子の散乱行程に入る。以上の行程をすべての光子において行うことで、最終的に透過光強度分布が得られる。上述のSD-OCTによる実験を行い散乱係数を求める。別の光学系による実験により得られた散乱体からの出力光強度分布からビーム径を求める。そして、求めた散乱係数を用いたモンテカルロシミュレーションによるビーム径と比較することで、非等方散乱因子を決定する。以下、詳細に説明する。散乱係数を実験的に求めるSD-OCTの光学系について、図5を参照して説明する。図5は、SD-OCTの光学系を示している。光源として、中心波長970nm、バンド幅20nm、出力3.0mWのスーパールミネッセントダイオード以下、SLDを用いた。波長範囲が780~1000nm、アイソレーションピークが39dB以上のアイソレータを、反射による戻り光をカットするために用いた。また、波長範囲が810~890nm、結合比50/50のファイバカップラを用いて、SLD光を参照側と物体側に分波し、反射してきた参照光と物体光を合波した。対物レンズL5として、光学倍率が1倍、NAが0.03、焦点距離が83mmのものを用いた。ここで、物体側は対物レンズにより波長による群速度分散が起こるので、参照側に補償板BK7を用いている。参照ミラーには、有効径が25.4mm、偏角精度が1.0秒、波面精度が0.10λの金コーティングされたミラーを用いた。コリメート/フォーカシングレンズL1~L4として、焦点距離が18.4mm、NAが0.15、ARコーティング波長が600~1500nmのものを用いた。スペクトロメーターとして、スリット幅が50μm、波長範囲が810~987nm、波長分解能が0.50nm程度のものを用いた。金コーティングミラー以下、Auミラーをサンプルとして、参照ミラーに対する光軸方向の相対位置を、それぞれ0.1mm、0.4mmに設定した状態で測定した干渉スペクトル信号を図6に示す。図6aは相対位置を0.1mmに設定した状態で測定したものであり、図6bは相対位置を0.4mmに設定した状態で測定したものである。測定条件は、スペクトロメーター露光時間が4msec、スペクトロメーターアベレージ回数が16回であり、アイソレータ直後に光学濃度以下、ODが3.0のNDフィルタを配置して減光した状態で測定した。光源のスペクトルがガウス形状であるのに対して、これらには変調成分が出ていることが確認できる。また、基準位置となる参照面から物体の反射面までの距離が離れていった場合に、変調周期が短く、変調度も小さくなることが確認できる。複数の反射面が同軸上に存在する場合、それぞれの反射面での変調周期が重なった複雑な干渉スペクトルが得られる。この干渉スペクトルに対して、波長を角周波数に換算して補間と再サンプリングを行った後、ブラックマン窓を乗じる。さらに、逆高速フーリエ変換を実行して、距離に換算し、その包絡線を検出することで、反射面の深さ位置に対する振幅反射率のプロファイルが得られる。また、Auミラーをサンプルとして、参照ミラーに対する光軸方向の相対位置を0.4mmに設定した状態で測定した反射面に対応する信号プロファイルを図7に示す。ここで、参照ミラーに対するサンプルの光軸方向での相対位置を0.4mmに設定するのは、相対位置が0mm付近ではノイズが乗りやすく、得ようとする信号に混ざってしまうためである。測定条件は、スペクトロメーター露光時間が4msec、スペクトロメーターアベレージ回数が128回であり、サンプルの同じ点に対して20回測定を行った。また、離散フーリエ変換以下、DFTデータ点数Nは65536個である。このプロファイルは20回行った測定の平均値をとり、振幅反射率を規格化したものである。これより半値全幅での深さ分解能が35.1μmであることがわかった。三次元光散乱体として、無色透明のPMMAポリメタクリル酸メチル材料の内部にSiO2粒子径が300nmと1μmの2種類異なるものをそれぞれ作製し、これをサンプルとした。参照ミラーに対するサンプル表面の光軸方向での相対位置が0.4mmになるよう設定した。また、深さ方向の信号減衰特性を考慮して、対物レンズを試料表面から内部方向に0.6mmの位置にデフォーカスした状態で測定を行った。測定条件は、スペクトロメーター露光時間が4msec、スペクトロメーターアベレージ回数が128回であり、試料横方向に50μmピッチで20点、試料縦方向に50μmピッチで20点の計400点で測定を行った。また、DFT後の再生深さ方向のデータを1μmピッチで1200点とした。DFTデータ点数Nは65536個とした。1つの試料に対して、異なる3箇所で上記の測定を行った。例として、SiO2粒子径が300nmで、体積率がそれぞれ3.38、17.87である試料に対する、3箇所の平均値をとった深さ方向反射率プロファイルを図8に示す。ここで、図8aは体積率が3.38の場合、図8bは体積率が17.87の場合を示している。共に比較的大きな山が3つ確認できるが、深さ0mm、0.7mm付近の山が、それぞれ試料表面、裏面を表している。また、深さ0.4mm付近の山は、参照ミラーに対するサンプルの光軸方向での相対位置を0.4mmに設定したことに起因する自己相関信号である。例として、SiO2粒子径が300nmで、体積率がそれぞれa3.38、b10.45、c14.48、d17.87である試料の、信号減衰特性の補正後の深さ方向反射率プロファイルを重ねたものを図9に示す。単散乱領域において、信号はexp-2μsZで減衰することから、図9に示すグラフでは、縦軸単位がdBであり、直線とみなせる部分の傾きの大きさから散乱係数が求まることがわかる。Mie散乱理論やノイズの影響を考慮したうえで、直線に最も近似する範囲である、0.1~0.3mmの範囲を選択した。そして、この範囲にあるデータに対して線形最小二乗法を適用し、得られた直線の傾きの大きさから散乱係数値を算出した。上記のような手順で算出された散乱係数を下記表1に示す。表1から、SiO2体積率が大きくなるほど、散乱係数値が大きくなることが確認できた。また、SiO2粒子径が大きい方が、散乱係数値は大きくなることも確認できた。この傾向はMie散乱理論とも一致する。この算出された散乱係数をモンテカルロシミュレーションでの入力パラメータとして用いる。作製した光散乱体の評価手法として、出力光ビーム径を用いる。図10に出力光強度を測定する光学系を示す。入力光を波長594nmのヘリウムネオンレーザーとし、2つのNDフィルタを通して光強度を調節する。レーザー光は光散乱体を通過後、散乱体後面をレンズによりCCDイメージセンサに拡大結像して記録することで出力光強度分布を得る。出力光強度分布に対して、ガウシアンフィッティングを施したグラフの1/e2全幅を出力光ビーム径とする。光散乱体の両端と中央の3箇所にて観測を行い、その平均値を出力光ビーム径として用いた。作製した8種類の散乱体に対して、出力ビーム径を調べた。その結果のグラフを図11に示す。これらから、SiO2体積率が大きくなるほど、散乱により出力ビーム径が広がることを確認できた。また、SiO2粒子径の大きい方が、出力ビーム径が大きく、より散乱していることも確認できた。この結果は、Mie散乱理論と整合性がとれている。次に、モンテカルロシミュレーションを用いて出力光ビーム径を得て、実験値との比較から、非等方散乱因子を特定する。モンテカルロシミュレーションでは、膨大なサンプル数によりデータの平滑化が行われるため、透過光強度分布は滑らかな形状となる。このため、出力光ビーム径を得る際に、出力光強度分布に直接ガウシアンフィッティングを行い、1/e2全幅を計測する。下記表2に示す条件下で、SiO2粒子径が300nmで、SiO2体積率がそれぞれa3.38、b17.87の試料に対して、モンテカルロシミュレーションを行い、透過光強度分布を得た。モンテカルロシミュレーションでは、L.Wangらが開発したMCMLMonteCarlomodelingoflighttransportinMultiLayeredtissuesのソースコードを用いて製作した。なお、他の条件として、散乱係数はSD-OCTで求めた上記の表1の値を用いた。また、Numberofpixelalongdepthdirectionは、下記表3に示す各散乱体の厚さから算出される値を用いた。非等方散乱因子を、0.1~0.9間において0.1間隔で、0.91~0.99間において、0.01間隔で変えてビーム径を算出した。透過光強度分布から得られた出力光ビーム径のグラフを図12に示す。また、下記表4に、SiO2粒子径が1μmでのサンプルでの散乱係数と非等方散乱因子の実験測定結果を示す。また、下記表5に、SiO2粒子径が1μmでのサンプルでの散乱係数と非等方散乱因子のミー散乱理論での結果を示す。上記の表4の結果は、SiO2粒子径が1μmでの散乱係数μsと非等方散乱因子gの結果である。また、上記の表5の結果は、ミー散乱理論による結果である。散乱係数は、最大2倍程度異なる結果となっているが、非等方散乱因子gは10程度の差であり、かつSiO2の体積濃度にも依らず一定値となっていることがわかる。非等方散乱因子gが一定であることは表5の結果からもわかることから、実験による散乱係数と非等方散乱因子はサンプルの作製精度から十分な結果であると考えられる。図13に示すフローは、本発明の三次元光散乱体の実効散乱係数の算定方法の他の実施形態であり、散乱光の透過/反射出力パワー比を指標として非等方散乱因子gを決定するものである。実施例2の算定方法は、実効散乱係数が未知の三次元光散乱体に対して、散乱係数μsを測定するステップと、算出した散乱係数μsに基づき、計算機シミュレーションを用いて散乱光の伝搬における非等方散乱因子gを0~1まで変化させて、散乱光の透過/反射出力パワー比との相関関係を算出するステップと、三次元光散乱体にビーム光を照射した際の散乱光の透過/反射出力パワー比の実測値と、上記相関関係とから、非等方散乱因子gを決定するステップを備えるものである。図14に、実施例2における透過/反射パワーの測定光学系の模式図を示す。透過パワーの測定では、散乱体を積分球である検出器にはりつけ、波長633nmのHe-Neレーザーからの光を入射させ、透過散乱光パワーを測定する。また、反射パワーの測定では、放物面鏡を用い、散乱光を平行光に変換し、レンズで集光することで積分球である検出器で反射パワーを測定する。この2つの比をとることで、透過/反射パワー比を得る。図15に、実施例2における透過/反射パワー比の実験値とモンテカルロシミュレーションによる結果から導出された非等方散乱因子の結果を示す。実験で得た散乱係数μsを用いて、非等方散乱因子gを0から1まで変化させて、それぞれモンテカルロシミュレーションによる光子の伝搬計算を行い、透過/反射出力パワー比を得る。そのグラフを図14に示す。実験で得られた透過/反射パワー比となるgをグラフから求めることで、非等方散乱因子gを決定した。下記表6に、SiO2粒子径が1μmでのサンプルでの透過/反射パワー比からの非等方散乱因子の実験結果とミー散乱理論での結果の比較を示す。透過/反射パワー比から求めた非等方散乱因子はMie散乱理論値と近い値となっており、また、SiO2の体積濃度に対して一定値となっていることから信頼できる値と考えられる。