(57)【要約】【課題】層状体の各層の厚みと材質を把握することができる層状物体の測定装置および方法を提供する。【解決手段】測定装置は、層状体からの電磁波パルスの時間波形を測定した測定結果に基づいて、層状体が含んでいる複数の層のそれぞれの材質及び厚みを特定する。測定装置は、複数の材質候補のデータと、複数の厚み候補と、を記憶するデータベースと、データベースに記憶されている複数の材質候補のデータ及び複数の厚み候補の探索範囲を入力する入力部と、を有する。更に、探索範囲内の複数の材質候補のデータ及び複数の厚み候補を用いることにより層状体からの電磁波パルスの時間波形を再現して、再現した時間波形と測定結果の時間波形とを比較することで、複数の層のそれぞれの材質及び厚みを特定する処理部を備える。【選択図】図１
2014122875.tif 000002
【特許請求の範囲】
【請求項１】
層状体からの電磁波パルスの時間波形を測定した測定結果に基づいて、該層状体が含んでいる複数の層のそれぞれの材質及び厚みを特定する測定装置であって、
複数の材質候補のデータと、複数の厚み候補と、を記憶するデータベースと、
前記データベースに記憶されている前記複数の材質候補のデータ及び前記複数の厚み候補の探索範囲を入力する入力部と、
前記探索範囲内の複数の材質候補のデータ及び複数の厚み候補を用いることにより前記層状体からの電磁波パルスの時間波形を再現して、該再現した時間波形と前記測定結果の時間波形とを比較することで、前記複数の層のそれぞれの材質及び厚みを特定する処理部と、を備える、
ことを特徴とする測定装置。
【請求項２】
前記複数の材質候補のデータは、前記複数の材料候補の電磁波に対する特性を含む、
ことを特徴とする請求項１に記載の測定装置。
【請求項３】
前記複数の材料候補のデータは、前記複数の材料候補の複素屈折率スペクトルを含む、
ことを特徴とする請求項２に記載の測定装置。
【請求項４】
前記データベースに記憶されている前記複数の材料候補のデータは、ウェーブレットを用いて圧縮されている、
ことを特徴とする請求項１に記載の測定装置。
【請求項５】
前記処理部は、前記再現した時間波形と前記測定結果の時間波形とを時刻毎に比較することで、前記複数の層のそれぞれの材質及び厚みを特定する、
ことを特徴とする請求項１に記載の測定装置。
【請求項６】
前記処理部は、前記再現した時間波形をフーリエ変換することにより取得したスペクトルと前記測定結果の時間波形をフーリエ変換することにより取得したスペクトルとを周波数毎に比較することで、前記複数の層のそれぞれの材質及び厚みを特定する、
ことを特徴とする請求項１に記載の測定装置。
【請求項７】
前記層状体からの電磁波パルスは、前記層状体を透過した電磁波パルス、または、前記層状体で反射した電磁波パルスである、
ことを特徴とする請求項１に記載の測定装置。
【請求項８】
前記層状体からの電磁波パルスは、前記層状体で反射した電磁波パルスであり、
前記処理部は、前記層状体の表面からｍ番目の界面に最も早く反射されたパルスの推定検出時刻をｔｍとした場合に、前記時刻毎の比較の終了時刻をｔｍに設定して前記層状体の少なくとも（ｍ＋１）番目の層の材質及び少なくともｍ番目の層の厚みを特定するステップを、前記ｍの値を増加させながら繰り返すことにより、前記複数の層のそれぞれの材質及び厚みを特定する、
ことを特徴とする請求項５に記載の測定装置。
【請求項９】
前記処理部は、前記時刻毎の比較の開始時刻をｔｍ−Δτ０未満に設定する、
ことを特徴とする請求項８に記載の測定装置。
ただし、前記層状体に照射される電磁波パルスの時間波形のピークの半値全幅をＷとした場合に、３Ｗ≦Δτ０＜ｔｍを満たす。
【請求項１０】
電磁波パルスを前記層状体に照射する照射部と、前記層状体からの電磁波を検出する検出部と、を備える、
ことを特徴とする請求項１に記載の測定装置。
【請求項１１】
前記層状体は、電磁波パルスを透過し且つ複素屈折率及び厚みが既知の板材に、接しており、
前記照射部は、前記板材を介して前記電磁波パルスを前記層状体に照射し、
前記板材から反射された電磁波パルスと前記層状体から反射された電磁波パルスとを用いて、前記複数の層のそれぞれの材質及び厚みを特定する、
ことを特徴とする請求項１０に記載の測定装置。
【請求項１２】
前記層状体は、生体組織を含む、
ことを特徴とする請求項１に記載の測定装置。
【請求項１３】
層状体からの電磁波パルスの時間波形を測定した測定結果に基づいて、該層状体が含んでいる複数の層のそれぞれの材質及び厚みを特定する測定方法であって、
データベースに記憶されている複数の材質候補のデータ及び複数の厚み候補の探索範囲を取得するステップと、
前記探索範囲内の複数の材質候補のデータ及び複数の厚み候補を用いることにより前記層状体からの電磁波パルスの時間波形を再現するステップと、
前記再現した時間波形と前記測定結果の時間波形とを比較することで、前記複数の層のそれぞれの材質及び厚みを特定するステップと、を有する、
ことを特徴とする測定方法。
【請求項１４】
一時的でないコンピュータ読み取り可能な記憶媒体に格納され、請求項１３に記載の測定方法の各ステップをコンピュータに実行させる、
ことを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、電磁波パルスを対象に照射し、得られた電磁波パルスを解析して対象についての情報を取得する方法と、それを用いたトモグラフィ装置に関する。
【背景技術】
【０００２】
近年、周波数が30GHzから30THzの範囲におよぶ電磁波、いわゆるテラヘルツ波を用いた様々な検査技術が開発されている。テラヘルツ波の透過性を利用した非破壊検査はその一例である。特許文献１に挙げた例では、テラヘルツ帯の電磁波パルスを試料である塗装膜に照射し、反射してきた電磁波パルスを時間領域分光法（TDS法）により検出・解析する。検出信号に含まれるパルスのピークの時間差を調べることで、多層に形成された各塗装膜の膜厚を算出する。また、ピークの波形に着目し、その変化から各塗装膜の状態、例えば乾燥の様子を知ることも可能である。
【０００３】
一方、試料のテラヘルツ帯における吸収を利用し、分光の観点から試料の材質を調べる例もある。特許文献２では、テラヘルツ帯の電磁波パルスを試料に照射し、透過ないし反射してきた電磁波パルスのスペクトルを調べて試料を特定する。TDS法では電界の振幅と位相を同時に取得するため、前者のように各塗装膜の位置情報を得ることや、後者のように試料の光学特性を得ることが可能である。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開2004-28618号公報
【特許文献２】特開平10-153547号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら特許文献１の方法では、多層試料中の層の厚さが当該層中の電磁波パルスの幅と同程度に薄くなると、表裏の各界面に対応するピーク波形を分離することが困難である。また、層内での多重反射に由来するピーク波形が、たまたま他の界面由来のピーク波形と重なった場合も同様に分離が難しい。結果として、着目する層の厚みや材質の特定をうまく行えないことがあった。一方、特許文献２の方法では試料が多層化した場合に各層の材質を個別に特定することが難しい。このように、従来の多層試料の測定は、各層の厚みと材質を同時に把握することが困難であるという課題を有していた。
【課題を解決するための手段】
【０００６】
そこで本発明の一側面としての測定装置は、層状体からの電磁波パルスの時間波形を測定した測定結果に基づいて、該層状体が含んでいる複数の層のそれぞれの材質及び厚みを特定する。測定装置は、複数の材質候補のデータと、複数の厚み候補と、を記憶するデータベースと、前記データベースに記憶されている前記複数の材質候補のデータ及び前記複数の厚み候補の探索範囲を入力する入力部と、前記探索範囲内の複数の材質候補のデータ及び複数の厚み候補を用いることにより前記層状体からの電磁波パルスの時間波形を再現して、該再現した時間波形と前記測定結果の時間波形とを比較することで、前記複数の層のそれぞれの材質及び厚みを特定する処理部と、を備える。
【発明の効果】
【０００７】
本発明によれば、層状体の層の厚みと材質を把握することができる。
【図面の簡単な説明】
【０００８】
【図１】本発明に係る構成特定手法における、基本的な処理手順を示すブロック図
【図２】本発明に係る波形再構成手順を示すブロック図
【図３】本発明に係る波形の一致を評価する手順を示すブロック図
【図４】本発明に係る測定装置の一例を説明する図
【図５】本発明に係る測定装置の別の一例を説明する図
【図６】層状体の構成の一例を説明する図
【図７】データベースの基本的な構成を示す図
【図８】本発明に係る波形再構成手法を説明するための図
【図９】本発明に係る時間軸上での波形比較手法を説明する図
【図１０】実測した波形と算出した波形を比較した図
【図１１】本発明に係る周波数軸上での波形比較手法を説明する図
【図１２】主要なパルス波形の到達時刻を説明する図
【図１３】本発明に係る時間軸上での別の波形比較手法を説明する図
【発明を実施するための形態】
【０００９】
以下、図面を参照しつつ本発明の実施形態を説明する。
本実施形態では、テラヘルツ帯の電磁波パルスを層状の物体（以下、層状体）に照射し、TDS法で得た透過ないし反射後の時間波形を解析して層状体の構成、すなわち各層の材質と厚みを同時に処理部で特定する。図１に、その基本的な処理手順を示した。ここでは個々の装置やデータ処理手法の詳細には立ち入らず、層状体のある一箇所を測定する場合の手順について説明する。
【００１０】
層状体が材質の異なる複数の層を積層した構成であることと、各層を構成する材質の材料候補について予め光学特性（複素屈折率の周波数依存性）や厚みの代表値および範囲のデータ群からなるデータベースが用意されていることが前提である。データベースについては後述する。ステップＳ１０１で、これから測定する層状体に適したデータ群をデータベースからユーザーは入出力装置を使用して選択する。以後の層状体の構成の特定は、これら選択したデータの中から行う。
【００１１】
ステップＳ１０２では、参照波形eref(t)を測定する。ここで、参照波形とは層状体に照射した電磁波パルスe0(t)の時間波形であって、伝播経路（発生素子から参照試料までと参照試料から検出素子まで）の情報を含む。測定系によって取得方法は異なり、透過系では層状体無しの状態で測定した結果を用いる。一方、反射系では全反射する鏡面を試料表面と同じ位置に設置して測定した結果を用いる。次のステップＳ１０３では、対象の層状体を測定系にセットし、測定波形emeas(t)を測定する。
【００１２】
ステップＳ１０４では、先の参照波形eref(t)および測定波形emeas(t)をフーリエ変換し、それぞれのスペクトルEref(ν)およびEmeas(ν)を得ておく。ステップＳ１０５からＳ１０７は、解析のための繰り返し処理（ループ）である。まずステップＳ１０５で繰り返しの範囲を設定する。これは、先にデータ群を選択した際に絞られた材質の種類や各層の膜厚の範囲である。ついで、ループ内での解析に用いる候補のパラメータ、つまり具体的な各層の材質や厚みの値を決める。
【００１３】
つぎのステップＳ１０６では、先に定めたパラメータと参照波形eref(t)を用いて、層状体を測定した際に得られるはずの時間波形（以後、推定波形）eest(t)を再構成する。ステップＳ１０７では、再構成した推定波形eest(t)と実測して得た測定波形emeas(t)を比較する。このとき、適当な評価関数を用いて両波形の一致の程度を数値化し、パラメータの組み合わせと共に記憶する。評価関数には、例えば取りうる値の範囲（値域）が０以上で、両者が一致するほど小さな値を取る関数を選ぶ。具体的な例については後述する。
【００１４】
検討するパラメータの組み合わせが残っていればステップＳ１０５へと移り、終了したならばループを抜ける。ステップＳ１０８では、記憶しておいた評価関数の値をもとに最小値を取る、すなわち最も良く測定波形を再現するパラメータの組み合わせを抽出する。これが各層の材質と厚みを特定した結果となる。なお、この最小値が測定系のノイズレベルを考慮して設定した所定の値を超えた場合は、試料の構成が想定の範囲外であったとみなし、“NG（該当無し）”として記憶する。
【００１５】
最後のステップＳ１０９では結果を表示する。本手順は層状体の或る一点における深さ方向の構成情報を得るためのものだが、照射位置を変えつつ同様な測定を行えば２次元ないし３次元の物性の分布が得られる。これを適宜可視化することで、層状体の断面像や立体像（トモグラフィ像）を得ることができる。
【００１６】
前述したように、層状体を測定した時間波形（測定波形）から、その構成を直に求めることは難しい。一方、先のステップＳ１０６に代表される“逆の”処理、つまり既知の層状体の構成と参照波形とから測定されるはずの時間波形（推定波形）を得る処理は容易である。とはいえ、多数の層状体の構成の候補から最適な組み合わせを探す場合、特定の確度を上げるには相応の計算量（所要時間）が必要となる。そこで、測定する層状体についての事前の知識を積極的に利用し、探索の範囲を絞り込む。事前の知識としては層状体の層数、各層の材質、各層の代表的な厚み等が挙げられる。妥当な初期値を与え、探索に必要のない領域を除くことで、特定の確度を高め、所要時間の短縮を見込むことができる。
【００１７】
実際、測定する層状体が決まると、その代表的な層構成もある範囲に収まることが多い。例えば、層状体が塗装膜であれば各層の材質は各種の塗料である。各々について事前に光学特性を測定し、複数記憶しておいたものを使用する。同様に、厚みの代表値や取りうる範囲には過去の同様な層状体から得て記憶しておいた値を用いればよい。先に挙げたデータベースは、対象となる層状体以外の材質についても記憶する大規模なものでありうるが、ステップＳ１０１において候補として有望なデータ群（例えば多層の塗装膜に関連するデータ群）のみを選択し、探索範囲とする。本実施形態によれば、事前に得た様々な材質の光学特性や各層が取りうる厚みを格納したデータベースから適切なデータ群を選択し、測定結果を再現するパラメータをその探索範囲内で探索することで、層状体の各層の厚みと材質を同時に特定することができる。作業に必要な計算量（時間）も抑えられ、高速な処理を実現できる。
【００１８】
（実施例１）
ついで、本実施形態の第１の実施例について、図面を用いて詳細に説明する。代表的な装置構成や測定する層状体のモデル、データベースについて述べたあと、波形の再構成や一致の比較など処理の詳細について説明する。
【００１９】
図４に本実施例における測定系の一例を示した。層状体で反射した電磁波パルスを測定する、いわゆる反射系のテラヘルツ時間領域分光測定装置である。超短パルスレーザー１はフェムト秒オーダーの超短パルス光を出力する。ハーフミラー１２で分岐したパルス光の一方を光伝導素子２の裏面に照射し、テラヘルツ帯の電磁波パルスを発生させる。その強度は電源８より印加されるバイアス電圧の大きさに比例する。
【００２０】
光伝導素子２の表面側から取り出したテラヘルツ帯のパルスを、放物面ミラーで集光したのち層状体３に照射する。その放物面ミラーと光伝導素子２とは、照射部に相当する。層状体の構成を反映した変調を受けつつ反射されて戻ってきたテラヘルツ帯のパルスを放物面ミラーで集光し、光伝導素子５の表面に入射する。その放物面ミラーと光伝導素子５とは、検出部に相当する。一方、ハーフミラー１２で分岐した超短パルスレーザー光の片割れは、光遅延系４を経由したのち光伝導素子５の裏面に入射する。光伝導素子５は、表側に達した反射パルスの振幅強度に比例した電流を、裏側に超短パルスレーザー光が入射している間だけ出力する。
【００２１】
発振器７は電源８とロックインアンプ６に周期信号を供給する。光伝導素子２のバイアス電圧を変調しつつ、光伝導素子５で検出した信号をロックインアンプ６で位相敏感検波することで、テラヘルツ帯の反射パルスの振幅を高感度に検出する。制御系９は光遅延系４を制御して光伝導素子５の裏側に入射するパルス光の光路長を変えることで、表側に達したテラヘルツ帯の反射パルスを検出するタイミングを制御する。光路長を掃引すると同時にロックインアンプ６からの出力信号を取り込むことで、テラヘルツ帯の反射パルスの時間波形、つまり測定波形や参照波形を得ることができる。
【００２２】
参照波形の測定には、テラヘルツ帯の光を極力反射する鏡、たとえば金の薄膜を蒸着した鏡などを使用すると良い。このとき、鏡面と試料表面の位置は正確に一致、ないし所定の間隔だけ離しておく必要がある。このため、実験的には、複素屈折率が既知でかつ表面が平坦な材質（板材）を試料に押し当て、当該板材越しにテラヘルツ帯のパルスを照射・検出することも行われる。参照波形には当該材質と試料の界面からの反射光を利用するとよい。測定波形の解析は、テラヘルツ帯のパルスが当該材質を通過する際の影響を考慮したうえで行われる。
【００２３】
測定者は、入力や表示のための入出力装置（入力部）１０を介して測定の指示や画像形成などの結果の表示を表示部で行う。一方、データベース１１には予め材質の光学特性や寸法の情報が記憶されている。測定者が層状体に適したデータ群を入出力装置１０で指定すると、制御系９はデータベース１１から当該群中のデータを適宜入手し、測定波形の解析に使用する。
【００２４】
なお、本実施例ではテラヘルツ帯のパルスを発生・検出するために光伝導素子を用いるとしたが、同様な効果を持つ他の手段を用いてもかまわない。例えば、非線形光学効果を持つ材料（ＤＡＳＴ等の樹脂やＬｉＮｂＯ３等の結晶）に超短パルスレーザーを照射し、テラヘルツ帯パルスの発生源とすることができる。光伝導素子と比較して、高強度のパルス光を使える可能性がある。検出側についても、層状体からのパルス光を電気光学結晶（ＺｎＴｅ等）にプローブ光とともに照射して、偏光状態の変化から強度を読み取る構成としても良い。２次元での高速なイメージングを行える可能性がある。
【００２５】
また、処理部を含む制御系９、入出力装置１０、データベース１１は１台のコンピュータ（ＰＣ）上に構成しても、それぞれ別個の構成として必要な時に互いに接続される構成としても良い。コンピュータは、層状体の構成を特定するプログラムによって、必要なステップを実行させられる。後者の例として、大規模なデータベース１１の本体を測定系とは離れた所に置く構成が考えられる。この場合、制御系９とデータベース１１は、通信により適宜データをやり取りする。さらに、比較的計算機パワーを必要とする波形解析（再構成や一致判断）や画像処理などを、測定系外の高速な計算機で行う構成としても良い。本実施形態では、測定と解析の担当を分けるメリットが少なくないからである。例えば、測定系をコンパクトかつ高速に動作する構成とできること、測定者が層状体の候補を詳細かつ広範なデータベースから選択できること、高速・高度な処理が可能な計算機を使えること、などが挙げられる。
【００２６】
図６に本実施例で測定する層状体の一例を示す。
図６（ａ）は層状体断面の模式図である。層状体の表面から奥に向けて順に層１（６０１）、層２（６０２）、層３（６０３）…、が積層されている。各々、材質に応じた光学特性を持ち、その複素屈折率を順にｎ~1、ｎ~2、ｎ~3…、とした。各層の境界は、順に表面（６０４）、界面１（６０５）、界面２（６０６）、界面３（６０７）…、厚みはd1、d2、d3…、である。図６（ｂ）と図６（ｃ）は、層状体にテラヘルツ帯のパルスを照射した際の断面模式図である。図６（ｂ）は反射系配置、図６（ｃ）は後述する透過系配置での測定に各々対応している。層状体に照射したパルス（参照波形）eref(t)は、各界面（６０４〜６０７）に到達する度に透過光と反射光とに分岐し、各層（６０１〜６０３）を通過する度に吸収や分散の影響を受ける。測定波形emeas(t)は、層状体を反射（ｂ）ないし透過（ｃ）してきたそれらの波形の重ね合わせである。
【００２７】
本実施形態は、層状体が図６（ａ）様の構成を取ると仮定して、参照波形eref(t)と測定波形emeas(t)から各層の材質と厚みを特定する手法だということもできる。図７に本実施形態で用いるデータベースの基本的な構成を示す。本データベースは大別すると材質候補と厚み候補の２種類のデータセットから構成されている。
【００２８】
図７（ａ）は材質候補の光学特性に関するデータセットの概念図である。材質候補毎に、（１）材質に固有の材質番号、（２）材質の名称、（３）材質の複素屈折率スペクトル、の情報をそなえている。以上のデータを、少なくとも層状体を構成すると予想される材質について網羅するように揃えておく。周知のように複素屈折率ｎ~は実部の屈折率ｎと虚部の消衰係数κの和で表わすことができる。
【００２９】
【数１】
2014122875.tif 000003
【００３０】
そこで、（３）は屈折率ｎと虚部の消衰係数κに分けて記憶するとしても良い。例えば、多層の塗装膜であれば、各層に使われている塗料（名称：A、B、C等）について個別にテラヘルツ帯での複素屈折率を測定し、適宜付与した番号（１、２、３等）や属性情報（大分類、中分類等）とともにそのスペクトルをデータベースへと記録する。これらの事前測定に用いる装置は、塗装膜を評価する装置と同じでも別のものでも構わない。反射系の測定装置用に、透過系で測定・取得したスペクトルを使っても良い。ただし、（３）の帯域は測定系が発生・検出する電磁波パルスの帯域を含み、かつ波形再構成に必要な幅を最低限確保する。分解能についても同様である。繰り返し数を十分に取るなどして、正確なスペクトルを得ておくことが重要である。
【００３１】
図７（ｂ）は層状体の厚み候補に関するデータセットの概念図である。層状体の種類毎に、（４）層状体に固有の識別番号、（５）層状体の名称、（６）表面から奥へ順に付与した層番号、（７）当該層を構成する材質候補（材質番号の組）、（８）当該層が取りうる厚みの範囲、の情報をそなえている。こちらも、以上のデータを代表的な範囲を網羅するように揃えておく。
【００３２】
例えば、測定者がステップＳ１０１で層状体の候補に“多層塗装膜”を選択すると、１層から３層の各層について上記のように材質と厚みの範囲が決まる。各層の材質としては、幾つかの乾燥段階にある様々な塗料が挙げられる。以後は、これらの範囲内で測定結果の解析を進める。選択によって探索の範囲を制限することには、様々なメリットがある。まず、構成の特定に要する計算量を抑制できるため、装置の動作が高速化する。余力は各層の厚みを算出する精度の向上に振り向けても良い。また、複素屈折率は似ているが明らかに異なるし使われてもいない材質を範囲から除くことは、誤検出を避けることにもつながる。
【００３３】
他の層状体として、紙やプラスチック、布や樹脂等、広範な物体が候補である。具体的には、冊子、封筒、包装された箱やカード類、電線等の絶縁被膜などに加え、生体組織も対象の１つである。特に皮膚や（臓器の）上皮組織は層状の構成を有しており、その厚みや材質は、健康を含む様々な情報を与える。そこで、たとえば正常な皮膚組織と異常な皮膚組織（乾燥、炎症、細菌に感染、がん）について事前にスペクトルを得ておき、厚みの範囲等の情報と、“皮膚”という名称とともにデータベースに格納しておく。内臓についても同様で、それらのスペクトルは切除した切片組成の臨床解析や、内視鏡などによる組織表層の観察の際に使われる。
【００３４】
スペクトル情報の格納にあたっては、段階的・階層的な分類をしておくことが望ましい。たとえば、“皮膚”という大分類のもと、“皮膚がん”、“皮膚炎”などの中分類があり、“皮膚がん”の下には種類や進行度、悪性か良性かなどの小分類がある、といった具合である。これら分類用の情報は属性情報として付番・記憶される。
【００３５】
これらの分類は、測定の精度と速度を上げるために効果的である。たとえば、皮膚を対象とする場合、ステップＳ１０１でまず大分類の“皮膚”を選択して測定し、観察部位の概要を把握する。結果、各層の厚みの範囲や材質の候補などが分かり、探索の範囲を絞ることができる。次に中分類の“皮膚がん”を選択して測定し、がんの種類や程度、分布に関して詳しい情報を入手し、診断の助けとする。
【００３６】
ついで、既知の層状体の構成と参照波形eref(t)から推定波形eest(t)を再構成する手順、すなわち図１におけるステップＳ１０６について説明する。幾つかの手法が使えるが、ここでは伝送行列（Transfer Matrix）を用いて層状体の複素振幅反射率を求め、推定波形eest(t)を再構成する方法について説明する。
【００３７】
図２に、ステップＳ１０６の波形再構成の手順を示す。このうち、図２（ａ）が伝送行列を用いて再構成を行う手順である。なお、図２（ｂ）、図２（ｃ）は他の手法による手順なので、後の実施例中で説明する。ステップＳ１０６はループ内の１処理であり、ステップＳ１０５において層構成のパラメータが適宜決定されてから開始される。ステップＳ２０１では、再構成に必要な各層の複素屈折率と厚み{n~j,dj}j=1,2,…がこれを受けて定まる。  
【００３８】
ステップＳ２０２では、上記パラメータを用いて層状体全体の複素反射率スペクトルR(ν)が算出される（詳細は後述）。ステップＳ２０３では、これに参照波形のスペクトルEref(ν)を乗じることで、推定波形のスペクトルEest(ν)を得る。最後に、ステップＳ２０４において算出したスペクトルEest(ν)を逆フーリエ変換し、推定波形eest(t)を得る。
【００３９】
層状体全体の複素反射率スペクトルR(ν)は、各層を通過する際に受ける寄与と各界面での反射および透過からの寄与を掛け合わせることで求めることができる。図８を用いて順に説明する。簡単のため、テラヘルツ帯のパルスは層状体に垂直に入射するものとした。
【００４０】
周波数ν、時刻τ、進行方向x、波数kの進行波E+(ν)と後退波E-(ν)を次のように定義する。
【００４１】
【数２】
2014122875.tif 000004
【００４２】
ここで、波数kは光速をc、複素屈折率をｎ~としたとき以下のように求まる。
【００４３】
【数３】
2014122875.tif 000005
【００４４】
図８（ａ）は、複素屈折率ｎ~m、厚さdmの層mを、進行波E+(ν)と後退波E-(ν)が通過する様子を示している。界面mにおける電磁波にプライム’を付与して区別することにすれば、界面m-1と界面mにおける電磁波の振幅の間には以下の関係が成り立つ。以後、層mについての伝送行列をA(m)と表す。
【００４５】
【数４】
2014122875.tif 000006
【００４６】
図８（ｂ）は、複素屈折率ｎ~mの層ｍと同ｎ~m+1の層m+1との界面mを挟んで、進行波E+(ν)と後退波E-(ν)が相対する様子を示している。層m+1側の電磁波にプライム’を付与して区別すれば、層mと層m+1における電磁波の振幅の間には以下の関係が成り立つ。以後、界面mについての伝送行列をB(m)と表す。
【００４７】
【数５】
2014122875.tif 000007
【００４８】
層0を大気とし、層１から層mまである層状体全体の伝送行列Mは、各層および界面の伝送行列の積で表すことができる。
【００４９】
【数６】
2014122875.tif 000008
【００５０】
層状体を大気中に配置する場合、界面０（表面）を原点に取ると計算上便利である。その際、A(0)は単位行列となるので、式６より除いても良い。層状体前後の進行波、後退波の間には次式が成り立つ。各要素は、ステップＳ２０１で得たパラメータを式４と式５に代入し、式６を計算することで算出することができる。
【００５１】
【数７】
2014122875.tif 000009
【００５２】
反射系の場合、最終層が十分厚いなど最終界面に後退波が入射することが無い場合は、E’−=0として良い。その際の複素振幅反射率スペクトルR(ν)は次式として求めることができる。
【００５３】
【数８】
2014122875.tif 000010
【００５４】
一方、試料が鏡面（基板）の上に乗っている特別な場合ではE’+=E’−である。その際の複素振幅反射率スペクトルR(ν)は、以下のようになる。 
【００５５】
【数９】
2014122875.tif 000011
【００５６】
最後に、推定波形eest(t)と測定波形emeas(t)の一致度を評価する手順、すなわち図１におけるステップＳ１０７について、図３と図９を用いて説明する。
【００５７】
図３（ａ）は、時間軸上で推定波形eest(t)と測定波形emeas(t)を比較する場合の手順を示している。なお、図３（ｂ）については他の実施例の中で説明する。
【００５８】
図９は、実際に時間軸上で両波形を比較する際の概念図である。上から順に、推定波形、測定波形、推定波形と測定波形の差分を、時間軸上にプロットした。図９に示すように、通常、各波形は信号がほとんどない領域を含んでいる。これは主として測定上のマージンを確保するためだが、推定波形と測定波形の比較を精度よく行う目的からすると不要である。そこで、ステップ３０１により注目する領域を絞り込む。具体的な範囲は層状体の厚さから見積もれば良いが、両波形の一致度を評価している間は（範囲を）固定しておくことが望ましい。
【００５９】
ついで、ステップＳ３０２では、時間軸上での測定波形emeas(t)と推定波形eest(t)の差異δ(ｔ)を計算する。最後に、ステップＳ３０３で両パルス波形の一致度を求める。再構成に使用した材質と構成が後に再現できるよう対応付けを行ったうえで、評価関数Δに格納する。
【００６０】
【数１０】
2014122875.tif 000012
【００６１】
具体的には、ステップＳ３０１で設定した評価範囲内において、差異δ(ｔ)の累積二乗和を範囲内の測定箇所数Ｎで除し、想定する層状体の材質と構成に対応した添字で指定した先に結果を格納する。定義よりΔは常に正で、両波形が一致するほど値は小さい。評価関数Δには式１０以外にもバリエーションが考えられる。例えば、差異δ(t)の二乗和のかわりに絶対値を累積し測定箇所数Nで除してもよい。また、生体などテラヘルツ帯の光の減衰が大きな試料の場合には、試料表面に近い層からの信号ほど強度が大きい。そこで、指数関数的に変化する窓関数を適用してから式１０右辺を計算するとしてもよい。
【００６２】
ここで、上記手法の有効性を示す目的で、実際の測定結果をもとに推定波形を求めた例を紹介する。測定系は反射系のテラヘルツTDS測定装置、試料には３層からなる層状体を用いた。その構成を表１に示す。
【００６３】
層  ：  材質    厚み  （手法）
１層目：  石英  1010 um （実測）
２層目：  大気    50 um （シム厚さ）      （表１）
３層目：  樹脂  1020 um （実測）
【００６４】
１層目は厚み約1mmの平坦な石英板、３層目はこれも厚み約1mmの、テラヘルツ光を良く透過する樹脂性の板である。２層目は大気であって、上記石英板と樹脂板との間に市販の厚さ50umのドーナツ状のシムをスペーサとして挟み込むことで形成した。測定時には、この試料を大気中に保持し、石英板側からシム中央の空隙部に向けてテラヘルツ帯のパルスを照射した。
【００６５】
図１０に測定結果および本手法により求めた推定波形を示す。横軸は時間、縦軸は振幅であり、重なりを避けるように各時間波形を配置した。上から順に参照波形eref(t)、推定波形eest(t)、測定波形emeas(t)、推定波形と測定波形の差分δ(ｔ)である。参照波形eref(t)には、（表面位置が試料と共通になるよう配置した）Auミラーを測定した結果を用いた。
【００６６】
ここで、最上段の参照波形は自身の最大値、２段目以降の波形、測定波形、差分δ(ｔ)は推定波形の最大値でそれぞれ規格化している。一見して、推定波形が測定波形をほぼ再現しており、両者の差分もほとんど０に近いことが分かる。構成の特定にあたっては、事前の測定で得ていた様々な材質の複素屈折率のデータから、３種類の石英と２種類の樹脂に大気（n=1、κ=0）を加えた、合計６種類のデータ群を選択した。また、各層の厚みの範囲は実測値を参考に定めた。
【００６７】
測定波形を最も良く再現した波形から特定した、試料の構成を表２に示す。
層  ：  材質    厚み  （探索範囲）
１層目：  石英    992um（980〜1020um）
２層目：  大気    64um（40〜80um）      （表２）
３層目：  樹脂    998 um（980〜1020um）
表１と比較すると、各層の材質を正しく特定できているほか、厚みも近い値であることが分かる。
【００６８】
図１０を解釈すれば、ピーク（１）は試料表面、（３）は裏面からの最大の反射である。一方、ピーク（２）では石英と大気、大気と樹脂の両界面からの反射が重なりあっている。こうした事情は見ただけでは分かり難い。本手法を用いれば、試料各層の材質と厚みを同時に把握することができる。
【００６９】
以上、推定波形の再構成と、測定波形との一致度の評価方法について説明した。こうした比較を行うには、推定波形と測定波形の時間スケール（データ点数および時間間隔）が一致している必要がある。そのためには、データベース１１に格納した複素屈折率スペクトルｎ~(ν)の周波数スケールと、参照波形のスペクトルEref(ν)の周波数スケールとが一致していることが望ましい。毎回同じ条件で測定するならば、その時間スケールに合わせた周波数スケールで材質の複素屈折率スペクトルを用意し、データベース１１に格納しておけば良い。しかし、一般には測定点数や時間間隔といった条件は変わるため、計算の過程で適宜データの補間や近似、間引を行う必要がある。
【００７０】
そこで、例えばステップＳ２０３においてスケールの変換を行う。同じ時間スケールで測定した測定波形と参照波形のスペクトルをEmeas(ν)およびEref(ν)とする（ステップＳ１０２〜Ｓ１０４）。全系の複素振幅反射率R(ν)は、データベースに格納した複素屈折率スペクトルｎ~(ν)を用いて算出する（ステップＳ２０１〜Ｓ２０２）。ついで、複素振幅反射率R(ν)の周波数スケールを参照波形スペクトルEref(ν)のそれに一致させるよう、補間や近似を行う。調整後の反射率Rmod(ν)を参照波形スペクトルEref(ν)に乗じ、推定波形のスペクトルEest(ν)を得る（ステップＳ２０３）。以後は同様である。 
【００７１】
ここまで、データベース１１に格納した複素屈折率スペクトルのデータ（以下、スペクトルデータ）は、単純な表形式を想定していた。しかし、精度の高い生のスペクトルデータを多数の材質について記憶しておくことは、実際に使用する装置にとって、記憶容量や計算速度、費用等の面で望ましくない。そこで、データベース１１には各種材料の複素屈折率のデータを予めウェーブレット等で圧縮したのち格納しておく構成としても良い。ステップＳ１０１において、データベースからデータ群を選んだ後に、圧縮しておいたそれぞれのデータを解凍し、波形再構成用のスペクトルデータとして使用する。あるいはステップＳ２０１において、波形再構成用に必要となるその都度、必要なスペクトルデータを解凍する。解凍と同時に、測定波形と時間スケールを一致させるスケール変換を実施しても良い。その際、上記ステップＳ２０３における時間スケール合わせは不要となる。
【００７２】
また、データベース１１と制御系９は物理的に離れていても良い。たとえば、データベース１１には精度の高い、帯域も広くとったスペクトルデータを多数格納しておく。ステップＳ１０１において測定者が入出力装置１０で対象を指定すると、制御系９は必要なデータ群を適当なネットワーク経由でデータベース１１に要求する。データベース１１は、必要なスペクトルデータ（群）を、要求された時間スケールと帯域を持つよう適宜変換したうえで制御系９へと送り返す。別の例では、データベース１１は追加・変更可能な構成である。測定者は、目的に応じて外部から精度や帯域が整えられたスペクトルデータを通信や記憶媒体経由で入手し、データベース１１を更新する。
【００７３】
また、これまでは簡単のため、層状体の各層はそれぞれ均一な組成を持つかのように説明してきた。しかし、ある層におけるテラヘルツ光の散乱や、表面（界面）粗さに起因する散乱についても、伝送行列への近似により本手法へと取り込むことができる。
【００７４】
（実施例２）
本実施形態の第２の実施例について、図５を用いて説明する。先の実施例は反射系のテラヘルツ時間分光測定装置であったが、本実施例は透過系の場合である。先の実施例（図４）と比較すると、層状体３の周辺の配置が異なる。以後、共通な部分の説明は極力省き、構成や動作の差異に着目して説明する。
【００７５】
本実施例では、光伝導素子２の表面側から取り出したテラヘルツ帯のパルスを集光したのち、層状体３に照射する。層状体３を透過したパルスを再び集光して光伝導素子５の表面に入射する。時間波形の取得は先の実施例と同様に行うが、本実施例は先の反射系と参照波形eref(t)の取得方法が異なる。つまり、層状体３を装置から取り除いた“空の”状態で測定した結果を参照波形eref(t)として用いる。
【００７６】
また、推定波形eesf(t)を再構成するステップＳ１０６にも違いがある。先の実施例ではステップＳ２０２で層状体全体の複素振幅反射率スペクトルR(ν)を求めたが、本実施例では透過率スペクトルT(ν)を求める。その後、ステップＳ２０３で透過率スペクトルT(ν)を参照波形のスペクトルEref(ν)に乗じ、ステップＳ２０４でそれらに逆フーリエ変換を施すことで、最終的に推定波形eest(t)を得る。
【００７７】
層状体全体の複素振幅透過率スペクトルT(ν)は、複素振幅反射率スペクトルR(ν)と同様に、伝送行列を用いて求めることができる。具体的には、層状体最終層界面の位置における、層状体有りと無しの場合での電界強度の比を取れば良い。最終層の界面に後退波が入射することは無いとして、層状体前後の電界強度（式７）と大気中での位相変化から次式を得る。
【００７８】
【数１１】
2014122875.tif 000013
【００７９】
これ以後のステップ（Ｓ１０７〜）については、実施例１と同様なので省略する。
【００８０】
以上、本実施例によれば、透過系の測定においても本手法を用いた解析を行うことができる。対象の層状体によっては反射系より透過系の方が測定に適している。双方の配置で得た測定結果から相補的な情報が引き出せる可能性もあり、本実施例の有効性を高めている。
【００８１】
（実施例３）
先の実施例１および２では、ステップＳ１０７の波形再構成において伝送行列を用いた計算手法を使った。本実施例では、別の波形再構成手法を用いた場合について説明する。
【００８２】
層状体内のパルス光は、界面に達する毎に振幅を減じながら分岐していく（図６（ｂ）や（ｃ））。本実施例では、それらのパルスを個別に計算して足し合わせ、波形を再構成する。図２（ｂ）にその手順を示した。
【００８３】
ステップＳ２０１において、（あるループ条件内での）各層の複素屈折率と厚みの値が定まる。ついでステップＳ２０５で、各層を横断する際の位相項（吸収や分散の影響）や各界面の複素振幅透過率、ないし反射率を計算する。前者は図８（ａ）の配置において式４より、後者は図８（ｂ）の配置において次式により求まる。
【００８４】
【数１２】
2014122875.tif 000014
【００８５】
【数１３】
2014122875.tif 000015
【００８６】
あるパルスの振幅強度は、その経路が定まれば計算できる。例えば、振幅E0+の参照波形emeas(t)が層３との界面で反射され層２を往復し、再び大気に戻って来た場合の振幅E−は、以下のようになる。
【００８７】
【数１４】
2014122875.tif 000016
【００８８】
こうした層状体を反射ないし透過してきたはずの波形（以下、子波形）は、無限に考えられる。そこで、ステップＳ２０６においてこれらの子波形Eest_i(ν),i=1,2,…を求めると同時に振幅をチェックし、その最大値が所定の閾値（E0+×10-4など）以下であれば以後の計算を取りやめる。ステップＳ２０７において、これら有効な子波形を全て加算し、推定波形のスペクトルEest(ν)を算出する。最後にステップＳ２０４でフーリエ逆変換をほどこし、推定波形eest(t)を得る。
【００８９】
生体をはじめとする減衰の大きな層状体では、子波形の振幅がすぐに減るため、有効な子波形自体が少ない。このような場合には、計算量を抑える面からも本実施例を用いると都合が良い。また、界面の荒れに由来する散乱の影響等、理想からの僅かなずれを導入する場合も、本実施例の方が容易である。
【００９０】
また別の波形再構成手法として、FDTD法など既知の電磁場解析手法を用い、反射ないし透過してきた時間波形を直接算出しても良い。図２（ｃ）にその手順を示した。
【００９１】
ステップＳ２０１において、（あるループ条件内での）各層の複素屈折率と厚みの値が定まる。ステップＳ２０８において、各層の厚みの値を参考に、計算を行う領域にメッシュを切る。続くステップＳ２０９では、計算を行う領域に媒質の物性値（複素屈折率）を定め、境界条件を定める。ステップＳ２１０では電磁波解析を行い、最終的にステップＳ２１１で推定波形eest(t)を得る。電磁場解析を用いて推定波形を算出する手法のメリットは、媒質の物性値が徐々に変化している場合など、様々な分布についても柔軟に対応できる点にある。また、時間波形を直接算出するため、参照・測定波形のスペクトルを計算する必要がない（ステップＳ１０４）。
【００９２】
（実施例４）
本実施形態の第４の実施例について、図面を参照しながら説明する。先の実施例では、ステップ３０１において、測定および推定の時間波形を比較する領域（評価範囲）を固定として説明した。本実施例は反射系の構成を必須とし、評価範囲を移動させる。本件は各層界面から来る時間波形を時間分解分光をもちいて直接測定する。このため、各時間波形は早期の領域ほど層状体表面に近い層を由来とする情報を含んでいる。そこで、評価範囲を時間の経過に沿って移動ないし拡大させつつ、試料構成を表面に近い側から順に推定していく。
【００９３】
図１２に、各界面から反射されてきたパルスの代表的な到達時刻を示した。例として、図６（ａ）の試料にテラヘルツ帯のパルスを照射した場合の反射光を想定している。最上段の測定波形は各界面から反射されてきたパルスの重ね合わせである。また、次段以降は、その代表的なパルスの到着時刻を推定し界面別に示したもので、それぞれ最も早い検出時刻をt0、t1、t2、t3とした。なお、時刻tmの値は、表面からのパルスを検出した時刻t0を基準（t0=0）として、以下で与えられる。
【００９４】
【数１５】
2014122875.tif 000017
【００９５】
一般に、検出時刻tm以前に検出した時間波形は、層ｍまでの情報しか含んでいない。そこで本件は、開始から時刻t1までの時間波形で層１の構成を推定し、ついで時刻t2までで層２の構成を推定し、というように表面から奥へと順に構成の推定を進める。すなわち、処理部は、層状体の表面からｍ番目の界面に最も早く反射されたパルスの推定検出時刻をｔｍとした場合に、時刻毎の比較の終了時刻をｔｍに設定して層状体の少なくとも（ｍ＋１）番目の層の材質及び少なくともｍ番目の層の厚みを特定する。この際、ｍの値を増加させながら繰り返すことにより、複数の層のそれぞれの材質及び厚みを特定する。
【００９６】
実際の手順は以下のようになる。
（1）マージンδτと、所定の時間間隔Δτ0を定める。
（2）選択したデータ群から構成の候補を１組選び、その構成における代表パルス到達時刻（t1、t2、・・）を求める。
（3）表面を対象に評価範囲を設定し波形を比較。候補を選んでは比較することを繰り返し、表面の位置と初層（第１層）の材質を推定する。
（4）界面１を対象に評価範囲を設定し波形を比較。同様に界面１の位置と第２層の材質を推定する。
（5）順次奥へと対象とする界面を移しつつ、同様に評価範囲を設定し波形を比較。対象を界面ｋとした場合、その位置と第ｋ＋１層の材質を推定する。
（6）最終層まで推定を終えたら終了する。
【００９７】
ここで、（1）のマージンδτとは時間波形を評価する際の余裕であり、使用する測定系に応じて予め定めておく。実際のテラヘルツ帯のパルスは、図１０の測定波形で見られるように、波形のピーク時刻から前後に広がりを持つ。そこで、パルスの真の開始時刻をより分かりやすいピーク時刻を介して算出するためにマージンδτを導入する。たとえば、図１０の時間波形を、裾を含めて評価するためにはδτ≒４ps程度が適当である。
【００９８】
一方、所定の時間間隔Δτ0は評価範囲の適正な幅を定める。測定および推定の両時間波形の比較を精度よく行うには、ピークのずれを検出するに十分な間隔を確保したうえで着目する界面近傍を評価すればよい。そのような適正な幅は測定系に依存し、照射するテラヘルツ帯のパルスの幅が目安となる。ピークの時間軸上におけるずれを前後どちらであっても検出するという要請から、パルスの半値全幅（FWHM）Ｗの３倍以上、好ましくは３０倍程度の幅を所定の時間間隔Δτ0として採用する。したがって、Δτ0は、３Ｗ≦Δτ0＜ｔｍを満たす。例えば、図１０の例において、テラヘルツ帯のパルスの半値全幅は約0.3psであった。そこで、少なくともΔτ0＝1ps、望ましくは10ps程度の幅を評価用に確保する。
【００９９】
評価範囲の設定は、本実施例で重要な部分であるので、図１３を用いて詳しく説明する。図１３（ａ）は、波形比較用の基本的な窓関数の例である。１か０の値を持ち、時間軸に沿ってステップ状に変化する。このような窓関数を時間波形に乗じることで、所望の評価範囲を切り出す。切り出し開始と終了時刻、および幅（τa1、τa2、Δτa＝τa2−τa1）は、場合に応じて以下のように決定する。
【０１００】
図１３（ｂ）は、先の手順（1）と（2）を終え、手順（3）の評価を行う様子を示している。まず、表面からのパルスの到達時刻t0からマージンをδτだけ見込み、評価範囲の開始時刻τb1を決定する（τb1=t0−δτ）。ついで、界面１からの代表パルス到達時刻t1から同様にマージンをδτ見込んで評価範囲の終了時刻τb2を求める（τb2=t1−δτ）。幅Δτｂ＝τb2−τb1が所定の時間間隔Δτ0より広い（Δτｂ≧Δτ0）ことを確認できたなら、この範囲を波形比較の評価範囲とする。評価の結果、特定されるのは、到達時刻に換算された表面の位置t0と、第１層の材質n~1（複素屈折率）である。それ以外のパラメータは本段階の評価にあたり無関係なので、適当な初期値を代入しておく。ついで、界面１の位置と第２層の材質の推定へと進む。
【０１０１】
図１３（ｃ）は、手順（4）のある段階を示している。界面１および界面２からの代表パルス到達時刻t1およびt2からそれぞれマージンδτを見込み、評価範囲の開始時刻τc1および終了時刻τc2を求める。その幅はΔτc＝τc2−τc1である。このとき、対象とする層（ここでは第１層）が薄いと、幅Δτcが所定の時間間隔Δτ0より狭い（Δτc＜Δτ0）ことがありうる。そうなるとパルスの分離ができず、十分な精度で評価することができないため、評価範囲を表面から離れる方向にもう１層分拡大することを考える。
【０１０２】
図１３（ｄ）は、評価範囲を拡大した手順（4）の別の段階を示している。界面３からの代表パルス到達時刻t3からマージンδτを見込み、評価範囲の終了時刻τd2を求める。開始時刻τd1は前回と同じであり（τd1＝τc1）、幅はΔτd＝τd2−τd1となる。その幅が所定の時間間隔Δτ0より広い（Δτd≧Δτ0）ことを確認できれば、この範囲を評価範囲として採用する。逆に、所定の時間間隔Δτ0よりも狭ければ（Δτd＜Δτ0）、同様な手順でさらに範囲を拡大する。このように、処理部は比較の終了時刻（例えばτd2）を起点として所定の時間間隔Δτ0より広い評価範囲が得られるよう開始時刻（例えばτd1）を設定する。すなわち、処理部は、時刻毎の比較の開始時刻をｔｍ−Δτ０未満に設定する。こうして採用した評価範囲は複数の界面と層を含むので、評価の際にはそれらを一括して推定する。本例ならば界面１と界面２の位置（時刻に換算するとt1、t2）、および第２層と第３層の材質n~2、n~3である。界面３以降の位置および第４層以降の材質は評価にあたり無関係なので、推定の対象から外してよい。
【０１０３】
本実施例のように、評価範囲を時間の経過順に移動させながら層状体各層の構成を推定する手法は、通常の手法と比較して推定が短時間で終了する。層数が多く構成が複雑であるほどその効果は高い。通常の手法では、各層の厚みや材質といったパラメータを、選択した範囲内で網羅的に変化させる。このため、検討すべきパラメータの組み合わせは、層数の増加にともない指数的に増加する。一方、本実施例の手法は、一度に推定するパラメータ数を数種類に抑えている。したがって、範囲の割り当て工程を余分に含むとしても、総数の増加にともなう計算の増加は線形的なものにとどまる。例えば、全５層の層状体について、各パラメータ（屈折率ｎ、消衰係数κ、厚みｄ）をそれぞれ10個の候補から推定する状況を考える。通常の手法では、計（103）5＝1015個のパラメータを一挙に評価して最適な構成を推定する必要があるのに対して、本手法は、たかだか103個程度のパラメータを５回評価すればよく、計算量を抑えることができる。
【０１０４】
（実施例５）
本実施形態の第５の実施例について、図面を参照しながら説明する。ここまでの実施例では、測定波形emeas(t) と 推定波形eest(t)の一致の程度を、時間軸上でのみ評価していた。時間波形のピーク位置はパルスが通過した全光路長に大きく依存する。つまり、時間軸上で一致を評価するのは層状体の複素屈折率のうち実部の屈折率差に着目した手法だと言える。一方、テラヘルツ帯における各材質の複素屈折スペクトルを見てみると、屈折率より吸収係数（消衰係数）の方に特徴が現れていることも多い。テラヘルツ帯に特徴的なスペクトル（いわゆる指紋スペクトル）を有する材質ならば、その識別は吸収係数（消衰係数）で行う方が効果的である。
【０１０５】
そこで、本実施例では測定波形と推定波形の比較を周波数軸上で行う。すなわち、前述の実施例では、再現した時間波形と測定結果の時間波形とを時刻毎に比較したが、本実施例では、それぞれの時間波形をフーリエ変換して取得したスペクトルを周波数毎に比較することで、複数の層のそれぞれの材質及び厚みを特定する。図１１にその概念図を、図３（ｂ）に一致度の評価手順をそれぞれ示した。ステップ３０４で、比較したい周波数範囲や領域を設定する。ついで、ステップＳ３０５では測定波形のスペクトルEmeas(t)と、推定波形のスペクトルEest(t)の差異δ（ν）を計算する。最後にステップＳ３０６で両周波数スペクトルの一致度を計算し、評価関数Δに格納する。ここで、添字のiやjは想定する各層の材質と構成に対応する番号で、Ｎは差異δ（ν）の評価範囲におけるデータ点数である。
【０１０６】
【数１６】
2014122875.tif 000018
【０１０７】
以降のステップは、従来と同じである。ここで、評価範囲（候補）の材質群Ａ、Ｂ、Ｃ、… が特徴的なスペクトル（吸収線）を持つ場合、一致の評価の範囲をその近傍に限ると良い。たとえば、材質Ａが周波数ν1とν2、材質Ｂがν3、… という具合に吸収線を持ち、適当な間隔をΔνとしたとき、式１６における積算の範囲を次の様にする。
【０１０８】
【数１７】
2014122875.tif 000019
【０１０９】
以上述べたように、周波数軸上で一致の評価を行うことで、より材質の特性を反映した判別が可能である。また、ここまでの実施例では測定波形と推定波形の一致の程度を調べる際に、各層について複素屈折率（屈折率と消衰係数）と厚みといったパラメータの範囲を定め、その範囲内でパラメータを振ることを想定していた。しかし、各パラメータを総当りに振るのでは、大きく外れた領域までも調べることになり、能率が悪い。
【０１１０】
一方、式１０や式１６で表される評価関数について、値を最小化するような複数のパラメータを探すという問題は、いわゆる最適化問題に他ならない。そこで、ニュートン法、共役勾配法、シンプレックス法等、既知の最適化手法を用いて最小値を取るパラメータを探すことができる。必要な初期値はデータベースを参考に与えるものの、範囲まで指定する必要は無くなる。
【０１１１】
（他の実施形態）
上述した実施形態は、システム或いは装置のコンピュータによって実現することもできる。システム或いは装置のコンピュータは、記憶媒体（例えば、一時的でないコンピュータ読み取り可能な記憶媒体）に格納された命令を読み出し実行することによって、前述した実施形態の機能を実現する。また、こうしたコンピュータによって実行される方法により実現することもできる。コンピュータは、１以上の中央演算ユニット（ＣＰＵ）、マイクロ処理ユニット（ＭＰＵ）、中央演算装置（ＣＰＵ）、その他の回路、別個のコンピュータやコンピュータプロセッサーのネットワークなどである。コンピュータにより実行可能な命令は、ネットワークや記憶媒体などからコンピュータに供給される。記憶媒体としては、ハードディスク、ＲＡＭ、ＲＯＭ、分散コンピュータシステムの記憶部、光ディスク（ＣＤなど）、ＤＶＤ、ブルーレイディスク、フラッシュメモリデバイス、メモリカード、などがある。
【符号の説明】
【０１１２】
Ｓ１０１…対象を選択する工程、Ｓ１０２…参照試料を測定し、時間波形を得る工程、Ｓ１０３…層状体を測定し、時間波形を得る工程
【図１０】
2014122875.tif 000020
【図１】
2014122875.tif 000021
【図２】
2014122875.tif 000022
【図３】
2014122875.tif 000023
【図４】
2014122875.tif 000024
【図５】
2014122875.tif 000025
【図６】
2014122875.tif 000026
【図７】
2014122875.tif 000027
【図８】
2014122875.tif 000028
【図９】
2014122875.tif 000029
【図１１】
2014122875.tif 000030
【図１２】
2014122875.tif 000031
【図１３】
2014122875.tif 000032