(57)【要約】【課題】テラヘルツ時間領域分光法による空間計測の分解能を向上させることができるようにする。【解決手段】計測対象物にテラヘルツパルス波を照射して計測対象物からの反射波の電界の計測を行い（Ｓ１）、計測によって得られた電界データに周波数フィルタを適用して電界データの周波数を限定し（Ｓ２）、周波数が限定された電界データを用いてスペクトル強度を計算する（Ｓ３）ようにした。【選択図】図１
2014130124.tif 000002
【特許請求の範囲】
【請求項１】
計測対象物にテラヘルツパルス波を照射して計測対象物からの反射波の電界の計測を行い、前記計測によって得られた電界データに周波数フィルタを適用して前記電界データの周波数を限定し、前記周波数が限定された電界データを用いてスペクトル強度を計算することを特徴とする面内空間分解能の向上方法。
【請求項２】
前記周波数フィルタとして以下の数式１に示す矩形フィルタ
【数１】
2014130124.tif 000016
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
をそれぞれ表す。
若しくは以下の数式２に示すガウスフィルタ
【数２】
2014130124.tif 000017
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
σ：フィルタ関数としての分散値
をそれぞれ表す。
を用いるようにし、数式１について０．５≦ａ≦１．２５且つ０．７５≦ｂ≦１．５（ただし、ａ＜ｂ）、或いは、数式２について０．５≦μ≦１．５且つ０．２５≦σ≦１．０であることを特徴とする請求項１記載の面内空間分解能の向上方法。
【請求項３】
計算された前記スペクトル強度に対して線形フィルタ，逆ガウス型フィルタ，上下左右方向エッジ強調マスク，上下左右斜め方向エッジ強調マスクのうちのいずれか一つを更に適用することを特徴とする請求項１記載の面内空間分解能の向上方法。
【請求項４】
計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置から読み込む手段と、前記電界データに周波数フィルタを適用して前記電界データの周波数を限定する手段と、前記周波数が限定された電界データを用いてスペクトル強度を計算する手段とを有することを特徴とする面内空間分解能の向上装置。
【請求項５】
前記周波数フィルタとして以下の数式３に示す矩形フィルタ
【数３】
2014130124.tif 000018
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
をそれぞれ表す。
若しくは以下の数式４に示すガウスフィルタ
【数４】
2014130124.tif 000019
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
σ：フィルタ関数としての分散値
をそれぞれ表す。
を用いるようにし、数式３について０．５≦ａ≦１．２５且つ０．７５≦ｂ≦１．５（ただし、ａ＜ｂ）、或いは、数式４について０．５≦μ≦１．５且つ０．２５≦σ≦１．０であることを特徴とする請求項４記載の面内空間分解能の向上装置。
【請求項６】
計算された前記スペクトル強度に対して線形フィルタ，逆ガウス型フィルタ，上下左右方向エッジ強調マスク，上下左右斜め方向エッジ強調マスクのうちのいずれか一つを適用する手段を更に有することを特徴とする請求項４記載の面内空間分解能の向上装置。
【請求項７】
計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置から読み込む手段、前記電界データに周波数フィルタを適用して前記電界データの周波数を限定する手段、前記周波数が限定された電界データを用いてスペクトル強度を計算する手段としてコンピュータを機能させることを特徴とする面内空間分解能の向上プログラム。
【請求項８】
前記周波数フィルタとして以下の数式５に示す矩形フィルタ
【数５】
2014130124.tif 000020
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
をそれぞれ表す。
若しくは以下の数式６に示すガウスフィルタ
【数６】
2014130124.tif 000021
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
σ：フィルタ関数としての分散値
をそれぞれ表す。
を用いるようにし、数式５について０．５≦ａ≦１．２５且つ０．７５≦ｂ≦１．５（ただし、ａ＜ｂ）、或いは、数式６について０．５≦μ≦１．５且つ０．２５≦σ≦１．０であることを特徴とする請求項７記載の面内空間分解能の向上プログラム。
【請求項９】
計算された前記スペクトル強度に対して線形フィルタ，逆ガウス型フィルタ，上下左右方向エッジ強調マスク，上下左右斜め方向エッジ強調マスクのうちのいずれか一つを適用する手段としてもコンピュータを機能させることを特徴とする請求項７記載の面内空間分解能の向上プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、面内空間分解能の向上方法、向上装置及び向上プログラムに関する。さらに詳述すると、本発明は、テラヘルツ時間領域分光イメージングによる断面構造の計測に用いて効果的な深度方向の空間分解能を向上させるための技術に関する。
【背景技術】
【０００２】
テラヘルツ波技術の進展により、テラヘルツ波イメージングによる非破壊検査は様々な分野で用いられ始めている。これまでに、遮熱材の欠陥検出，隠匿物検出，薬剤のコーティング厚さ分布計測，木材やコンクリートなどの建材中の欠陥検出，絵画の塗料分析などが報告されている。テラヘルツ波はプラスチック，樹脂，セラミック，木材などを比較的良く透過するため、対象物の断面構造（言い換えると、深度方向の分布）を計測することができ、これを面内方向の二次元分布計測と組み合わせることで三次元計測が可能となる。なお、０．１〜１０ THz の範囲がテラヘルツ波帯やテラヘルツ領域などと呼ばれる。
【０００３】
深度方向の空間分解能を得るには極めて短いパルスを用いて断面構造における各界面からの反射波を計測する。このようなパルスは１０〜１００ fs オーダのレーザ光を光伝導アンテナに照射することで発生させることができる。これらのパルスは、時間幅が短い分、広帯域の周波数特性を持つ。テラヘルツ時間領域分光法（ＴＨｚ−ＴＤＳ：Ｔerahertz Ｔime Ｄomain Ｓpectroscopy の略）ではこの特性を活かし、時間波形から対象物の吸収特性や屈折率などを得るようにしている（非特許文献１）。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】Ｙ．Ｌee：”Ｐrinciples of terahertz science and technology”，Ｓpringer，Ｎew Ｙork，２００９年
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、従来のテラヘルツ時間領域分光法では、これまでの実験結果から、反射波の振幅分布に基づく、つまり全周波数成分を含む場合の空間分解能は０．６ mm 程度であり、厚さ約１ mm のセラミック層を透視した場合には１ mm 程度が空間分解能の上限(限界)であることが分かっている（福地哲生 他：「テラヘルツ波を用いた遮熱コーティングの非破壊検査における空間分解能評価」，光応用・視覚，計測合同研究会資料，LAV-12-3/IM-12-11，pp.11-16，2012年；福地哲生 他：「テラヘルツ波を用いたガスタービン用遮熱コーティングの非破壊検査技術の開発（第2報）」，電力中央研究所報告H11004，2012年）。したがって、従来のテラヘルツ時間領域分光法は空間分解能が十分に高いとは言えない。
【０００６】
そこで、本発明は、テラヘルツ時間領域分光法による空間計測の分解能を向上させることができる面内空間分解能の向上方法、向上装置及び向上プログラムを提供することを目的とする。
【課題を解決するための手段】
【０００７】
本発明者は、テラヘルツ時間領域分光法による空間計測の分解能の向上についての検討を重ね、反射波の電界データ（具体的には、反射波の電界の時間変化を表す波形）の周波数の範囲を限定することによって空間分解能が変化することを突き止めた。
【０００８】
本発明の面内空間分解能の向上方法は、上記知見に基づくものであり、計測対象物にテラヘルツパルス波を照射して計測対象物からの反射波の電界の計測を行い、計測によって得られた電界データに周波数フィルタを適用して電界データの周波数を限定し、周波数が限定された電界データを用いてスペクトル強度を計算するようにしている。
【０００９】
また、本発明の面内空間分解能の向上装置は、計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置から読み込む手段と、電界データに周波数フィルタを適用して電界データの周波数を限定する手段と、周波数が限定された電界データを用いてスペクトル強度を計算する手段とを有するようにしている。
【００１０】
また、本発明の面内空間分解能の向上プログラムは、計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置から読み込む手段、電界データに周波数フィルタを適用して電界データの周波数を限定する手段、周波数が限定された電界データを用いてスペクトル強度を計算する手段としてコンピュータを機能させるようにしている。
【００１１】
これら本発明の面内空間分解能の向上方法、向上装置及び向上プログラムによると、計測によって得られた反射波の電界データの周波数を適切に限定した上でスペクトル強度を計算するようにしているので、空間分解能の向上に対して不必要なデータが排除され、空間分解能の向上に必要なデータのみによってスペクトル強度が計算される。
【００１２】
また、本発明の面内空間分解能の向上方法、向上装置及び向上プログラムは、周波数フィルタとして数式１に示す矩形フィルタ若しくは数式２に示すガウスフィルタを用いるようにし、数式１について０．５≦ａ≦１．２５且つ０．７５≦ｂ≦１．５（ただし、ａ＜ｂ）、或いは、数式２について０．５≦μ≦１．５且つ０．２５≦σ≦１．０であるようにしても良い。
【００１３】
【数１】
2014130124.tif 000003
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
をそれぞれ表す。
【００１４】
【数２】
2014130124.tif 000004
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
σ：フィルタ関数としての分散値
をそれぞれ表す。
【００１５】
これらの場合には、空間分解能の向上に対して不必要なデータが適切に排除されて空間分解能の向上に必要なデータのみが適切に抽出される。
【００１６】
また、本発明の面内空間分解能の向上方法、向上装置及び向上プログラムは、計算されたスペクトル強度に対して線形フィルタ，逆ガウス型フィルタ，上下左右方向エッジ強調マスク，上下左右斜め方向エッジ強調マスクのうちのいずれか一つを更に適用するようにしても良い。この場合には、コントラスト比を向上させ、空間分解能の更なる向上が図られる。
【発明の効果】
【００１７】
本発明の面内空間分解能の向上方法、向上装置及び向上プログラムによれば、空間分解能の向上に対して不必要なデータを排除し、空間分解能の向上に必要なデータのみによってスペクトル強度を計算することができるので、テラヘルツ波イメージングの空間分解性能を高め、テラヘルツパルス波を用いたイメージングによる非破壊検査技術としての有用性、延いては信頼性の向上を図ることが可能になる。具体的には、検査対象物の断面構造（即ち、深度方向の分布）の空間分解能を向上させて検査処理に対して有用な情報を提供することを可能にし、検査精度の向上が可能になる。
【図面の簡単な説明】
【００１８】
【図１】本発明の面内空間分解能の向上方法の実施形態の一例を説明するフローチャートである。
【図２】実施形態の面内空間分解能の向上方法を面内空間分解能の向上プログラムを用いて実施する場合の当該プログラムによって実現される面内空間分解能の向上装置の機能ブロック図である。
【図３】実施形態のテラヘルツ波イメージング装置の概略構成図である。
【図４】本発明で用いられ得る矩形フィルタ（周波数帯域：０．５〜１．０ THz）を示す図である。
【図５】本発明で用いられ得る矩形フィルタ（周波数帯域：１．０〜１．５ THz）を示す図である。
【図６】本発明で用いられ得るガウスフィルタ（中心周波数：０．７５ THz）を示す図である。
【図７】本発明で用いられ得るガウスフィルタ（中心周波数：１．２５ THz）を示す図である。
【図８】実施例１で用いた解像度ターゲット（テストパターン）を示す図である。
【図９】実施例１におけるテストパターンからの反射波の全周波数帯域におけるスペクトル強度を積算した結果を示す図である。
【図１０】実施例１における線群に対して垂直な直線上の強度分布の一例を示す図である。
【図１１】実施例１における周波数帯域０．５〜１．０ THz の矩形フィルタを適用した場合のスペクトル強度の計算結果を示す図である。
【図１２】実施例１における周波数帯域１．０〜１．５ THz の矩形フィルタを適用した場合のスペクトル強度の計算結果を示す図である。
【図１３】実施例１における中心周波数０．７５ THz のガウスフィルタを適用した場合のスペクトル強度の計算結果を示す図である。
【図１４】実施例１における中心周波数１．２５ THz のガウスフィルタを適用した場合のスペクトル強度の計算結果を示す図である。
【図１５】実施例１における線群に対して垂直な直線上の強度分布から求めたコントラスト比を計算した結果を示す図である。
【図１６】実施例１における周波数フィルタの適用結果に空間周波数フィルタ（逆ガウス型フィルタ）を更に適用した場合のスペクトル強度の計算結果を示す図である。
【図１７】実施例１における周波数フィルタの適用結果に上下左右方向エッジ強調マスクを更に適用した場合のスペクトル強度の計算結果を示す図である。
【発明を実施するための形態】
【００１９】
以下、本発明の構成を図面に示す実施の形態の一例に基づいて詳細に説明する。
【００２０】
図１から図７に、本発明の面内空間分解能の向上方法、向上装置及び向上プログラムの実施形態の一例を示す。
【００２１】
本実施形態の面内空間分解能の向上方法は、図１に示すように、計測対象物にテラヘルツパルス波を照射して計測対象物からの反射波の電界の計測を行い（Ｓ１）、計測によって得られた電界データに周波数フィルタを適用して電界データの周波数を限定し（Ｓ２）、周波数が限定された電界データを用いてスペクトル強度を計算する（Ｓ３）ようにしている。
【００２２】
また、本実施形態の面内空間分解能の向上装置は、計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置としてのデータサーバ(１６)から読み込む手段(１１ａ)と、電界データに周波数フィルタを適用して電界データの周波数を限定する手段(１１ｂ)と、周波数が限定された電界データを用いてスペクトル強度を計算する手段(１１ｃ)とを有する。
【００２３】
さらに、本実施形態の面内空間分解能の向上プログラムは、計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置としてのデータサーバ(１６)から読み込む手段(１１ａ)、電界データに周波数フィルタを適用して電界データの周波数を限定する手段(１１ｂ)、周波数が限定された電界データを用いてスペクトル強度を計算する手段(１１ｃ)としてコンピュータを機能させる。
【００２４】
そして、面内空間分解能の向上方法の実行にあたっては、まず、計測対象物にテラヘルツパルス波を照射することによる計測対象物からの反射波の電界の計測が行われる（Ｓ１）。
【００２５】
本発明で用いられ得るテラヘルツ波イメージング装置の一例を図３に示す。このテラヘルツ波イメージング装置は、レーザ２から出射された超短パルスレーザ光をハーフミラー８によって二つに分岐する。そして、分岐したうちの一方を送信用光伝導アンテナ３に入射して幅約１ ps のテラヘルツ波パルスを発生させ、このテラヘルツ波はレンズ４によって計測対象物１の表面上に集束され、計測対象物１から反射されたテラヘルツ波はレンズ４及び偏光子９を介して受信用光伝導アンテナ５に導かれて超短パルスレーザ光を入射した時点における電界を計測するようにしている。
【００２６】
テラヘルツ波の電界の時間変化を得るには、ハーフミラー８によって分岐したうちの一方のレーザ光であって送信用光伝導アンテナ３へ導入するレーザ光に対し、前記分岐したうちの他方のレーザ光であって受信用光伝導アンテナ５へ導入するレーザ光の時間遅延を走査する。この時間遅延の走査は可動ミラー６を用いて行う。本実施形態の装置では、可動ミラー６は１秒間に１００回走査されるため、サンプリング周波数は１００ Hz である。
【００２７】
受信用光伝導アンテナ５によって受信された反射波のテラヘルツ波信号は、必要に応じてアンプ（図示していない）によって増幅された後にＡ／Ｄコンバータでデジタル化され、保存される。可動ミラー６の位置を変えながら上記動作を繰り返すことによって、テラヘルツ波の時間波形として取得される。
【００２８】
また、送受信器（具体的にはアンテナ３，５など）は一体型のヘッドに搭載されて可動ステージ７に取り付けられている。そして、ヘッドを平面方向に走査することによってテラヘルツ波を照射する計測対象物１上の平面位置（ｘ，ｙ）を変えることができる。個々の平面位置（ｘ，ｙ）において可動ミラー６の移動による時間遅延の走査を行い、テラヘルツ波の時間波形を取得し、最終的には時間波形の２次元分布が得られる。
【００２９】
なお、本実施形態では図３に示すテラヘルツ波イメージング装置を用いるようにしているが、本発明で用いられ得る計測装置は図３に示すものには限られない。すなわち、テラヘルツ時間領域分光法において用いられ得るものであって、具体的には、テラヘルツ波を出射して反射波の電界（テラヘルツパルス波の振動電場の強さ）を計測し得るものであれば、どのような計測装置であっても構わない。
【００３０】
次に、Ｓ１の処理によって計測され取得された反射波の電界データに対して周波数フィルタを適用することによる電界データの周波数の限定が行われる（Ｓ２）。
【００３１】
ここで、本発明の面内空間分解能の向上方法におけるＳ２以降の処理は本発明の面内空間分解能の向上装置によって実行され得る。
【００３２】
そして、本発明の面内空間分解能の向上方法におけるＳ２以降の処理及びこれら処理を実行する面内空間分解能の向上装置は、本発明の面内空間分解能の向上プログラムをコンピュータ上で実行することによっても実現され得る。本明細書では、面内空間分解能の向上プログラムをコンピュータ上で実行することによってＳ２以降の処理を実行する面内空間分解能の向上装置が実現されると共に面内空間分解能の向上方法におけるＳ２以降の処理が実行される場合を説明する。
【００３３】
面内空間分解能の向上プログラム１７を実行するためのコンピュータ１０（本実施形態では、面内空間分解能の向上装置１０でもある）の全体構成を図２に示す。このコンピュータ１０（面内空間分解能の向上装置１０）は、制御部１１、記憶部１２、入力部１３、表示部１４及びメモリ１５を備え相互にバス等の信号回線によって接続されている。また、コンピュータ１０には記憶装置としてのデータサーバ１６がバス等の信号回線によって接続されており、その信号回線を介してデータや制御指令等の信号の送受信（即ち出入力）が相互に行われる。
【００３４】
制御部１１は記憶部１２に記憶されている面内空間分解能の向上プログラム１７によってコンピュータ１０全体の制御並びに面内空間分解能の向上に係る演算を行うものであり、例えばＣＰＵ（中央演算処理装置）である。
【００３５】
記憶部１２は少なくともデータやプログラムを記憶可能な装置であり、例えばハードディスクである。
【００３６】
メモリ１５は制御部１１が種々の制御や演算を実行する際の作業領域であるメモリ空間となるものであり、例えばＲＡＭ(Ｒandom Ａccess Ｍemory の略)である。
【００３７】
入力部１３は少なくとも作業者の命令を制御部１１に与えるためのインターフェイスであり、例えばキーボードである。
【００３８】
表示部１４は制御部１１の制御によって文字や図形等の描画・表示を行うものであり、例えばディスプレイである。
【００３９】
そして、本実施形態では、上述のＳ１の処理において計測され取得された反射波の電界データが反射波電界データベース１８としてデータサーバ１６に格納(保存)される。具体的には、原点(任意)からの横方向平面位置ｘ mm、縦方向平面位置ｙ mm における周波数ｆ THz 別の反射波の電界(電場)の強さＦ(ｘ，ｙ，ｆ) a.u. の値が反射波電界データベース１８に記録される。なお、データサーバ１６に格納される反射波の電界の強さＦ(ｘ，ｙ，ｆ) のデータは全ての周波数(ｆ)成分を含むものである。また、Ｓ１の処理において計測されたデータは、例えば、計測装置から適当な記憶媒体に保存されて作業者によってデータサーバ１６内の反射波電界データベース１８に記録・蓄積されるようにしても良いし、計測装置から通信手段（無線・有線）を介してデータサーバ１６内の反射波電界データベース１８に自動的に記録・蓄積されるようにしても良い。
【００４０】
そして、コンピュータ１０（本実施形態では、面内空間分解能の向上装置１０でもある）の制御部１１には、面内空間分解能の向上プログラム１７を実行することにより、Ｓ１の処理において計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置としてのデータサーバ１６から読み込む処理を行うデータ読込部１１ａと、電界データに周波数フィルタを適用して電界データの周波数を限定する処理を行う周波数フィルタ部１１ｂと、周波数が限定された電界データを用いてスペクトル強度を計算する処理を行うスペクトル強度計算部１１ｃとが構成される。
【００４１】
面内空間分解能の向上プログラム１７が実行されることによる具体的な処理としては、まず、コンピュータ１０（面内空間分解能の向上装置１０）の制御部１１に構成されたデータ読込部１１ａが反射波の電界データの読み込みを行う（Ｓ２−１）。
【００４２】
具体的には、データ読込部１１ａは、Ｓ１の処理において計測され取得されてデータサーバ１６に格納されている反射波電界データベース１８に記録されている反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータをデータサーバ１６から読み込む。
【００４３】
そして、データ読込部１１ａは、読み込んだ反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータをメモリ１５に記憶させる。
【００４４】
続いて、制御部１１の周波数フィルタ部１１ｂが、反射波の電界データに対して周波数フィルタを適用することによる電界データの周波数の限定を行う（Ｓ２−２）。
【００４５】
具体的には、周波数フィルタ部１１ｂは、Ｓ２−１の処理においてメモリ１５に記憶された反射波の電界の強さＦ(ｘ，ｙ，ｆ) のデータをメモリ１５から読み込む。
【００４６】
そして、周波数フィルタ部１１ｂは、周波数フィルタを適用して反射波の電界データの周波数帯域を限定する。本発明における周波数フィルタとしては、例えば、周波数軸(周波数の並び)に対して関数値が矩形の分布型を有する矩形フィルタ（数式３）やガウス分布型を有するガウスフィルタ（数式４）が用いられ得る。
【００４７】
【数３】
2014130124.tif 000005
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
をそれぞれ表す。
【００４８】
【数４】
2014130124.tif 000006
ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
σ：フィルタ関数としての分散値 [THz]
をそれぞれ表す。
【００４９】
なお、矩形フィルタにおける、フィルタ関数によって限定する周波数帯域の下限周波数ａ THz の値は０．５〜１．２５ THz の範囲で適宜選択されると共に上限周波数ｂ THz の値は０．７５〜１．５ THz の範囲で適宜選択され、０．７５≦ａ≦１．２５且つ１．０≦ｂ≦１．５が好ましく、１．００≦ａ≦１．２５且つ１．２５≦ｂ≦１．５がより一層好ましい（ただし、いずれについてもａ＜ｂ）。また、ｂ−ａ≧０．２５が好ましく、ｂ−ａ＝０．５がより一層好ましい。以上を踏まえて具体的には例えば、ａ＝０．５且つｂ＝０．７５，ａ＝０．７５且つｂ＝１．０が選択され得て、ａ＝０．５且つｂ＝１．０が好ましく（図４）、ａ＝１．０且つｂ＝１．２５，ａ＝１．２５且つｂ＝１．５がより一層好ましく、ａ＝１．０且つｂ＝１．５が最も好ましい（図５）。
【００５０】
また、ガウスフィルタにおける、フィルタ関数によって着目する周波数帯域の中心周波数σ THz の値は、０．５〜１．５ THz の範囲で適宜選択され、０．７５≦σ≦１．５が好ましく（σ＝０．７５：図６）、１．０≦σ≦１．５がより一層好ましく、σ＝１．２５が更に好ましい（図７）。
【００５１】
また、分散値σの値は、０．２５〜１．０ THzの範囲で適宜選択され、０．２５≦σ≦０．７５が好ましく、０．４≦σ≦０．６がより一層好ましく、σ＝０．５が更に好ましい。
【００５２】
そして、周波数フィルタ部１１ｂは、周波数フィルタを適用した後の反射波の電界の強さＦflt(ｘ，ｙ，ｆ)を数式５によって計算する。
【数５】
2014130124.tif 000007
ここに、Ｆflt(ｘ，ｙ，ｆ)：周波数フィルタ適用後の反射波の電界の強さ [a.u.]，
ｗ(ｆ)：フィルタ関数，
Ｆ(ｘ，ｙ，ｆ)：反射波の電界の強さ(計測値) [a.u]，
ｘ：原点(任意)からの横方向平面位置 [mm]，
ｙ：原点(任意)からの縦方向平面位置 [mm]，
ｆ：周波数 [THz]
をそれぞれ表す。
【００５３】
そして、周波数フィルタ部１１ｂは、計算した周波数フィルタ適用後の反射波の電界の強さＦflt(ｘ，ｙ，ｆ)の値をメモリ１５に記憶させる。
【００５４】
次に、Ｓ２−２の処理によって計算された周波数フィルタ適用後の周波数が限定された反射波の電界の強さＦflt(ｘ，ｙ，ｆ)のデータを用いてのスペクトル強度の計算を行う（Ｓ３）。
【００５５】
具体的には、スペクトル強度計算部１１ｃは、Ｓ２−２の処理においてメモリ１５に記憶された周波数フィルタ適用後の反射波の電界の強さＦflt(ｘ，ｙ，ｆ)の値をメモリ１５から読み込み、平面位置(ｘ，ｙ)別の全周波数成分を積算したスペクトル強度Ｓ(ｘ，ｙ)を、数式６によって計算する。
【００５６】
【数６】
2014130124.tif 000008
ここに、Ｓ(ｘ，ｙ)：スペクトル強度 [a.u]，
Ｆflt(ｘ，ｙ，ｆ)：周波数フィルタ適用後の反射波の電界の強さ [a.u.]，
ｘ：原点(任意)からの横方向平面位置 [mm]，
ｙ：原点(任意)からの縦方向平面位置 [mm]，
ｆ：周波数 [THz]
をそれぞれ表す。
【００５７】
そして、スペクトル計算部１１ｂは、Ｓ３までの処理による計算結果として、計算した平面位置(ｘ，ｙ)別のスペクトル強度Ｓ(ｘ，ｙ)の値を、表示部１４に表示したり、例えば記憶部１２やデータサーバ１６に計算結果データファイルとして保存したりする。
【００５８】
そして、制御部１１は、面内空間分解能の向上の処理を終了する（ＥＮＤ）。
【００５９】
以上の構成を有する本発明の面内空間分解能の向上方法、向上装置及び向上プログラムによれば、計測によって得られた反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータの周波数ｆを適切に限定した上でスペクトル強度Ｓ(ｘ，ｙ)を計算するようにしているので、空間分解能の向上に対して不必要なデータを排除し、空間分解能の向上に必要なデータのみによってスペクトル強度を計算することができ、テラヘルツ波イメージングの空間分解性能を高め、テラヘルツパルス波を用いたイメージングによる非破壊検査技術としての有用性、延いては信頼性の向上を図ることが可能になる。
【００６０】
なお、上述の形態は本発明の好適な形態の一例ではあるが本発明の実施の形態がこれに限定されるものではなく、本発明の要旨を逸脱しない範囲において種々変形実施可能である。例えば、上述の実施形態では反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータが蓄積される記憶手段をデータサーバ１６としているが、記憶部１２でも良いし、他の記憶装置を用いるようにしても良い。また、反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータを記憶装置に一旦蓄積することなく、計測装置からコンピュータ１０（面内空間分解能の向上装置１０）に直接入力し、入力されたデータから処理を順次行うようにしても良い。
【００６１】
また、上述の実施形態における周波数フィルタの適用結果に、他のフィルタ・他の画像処理を更に適用するようにしても良い。ここで、更に適用するフィルタ・画像処理は、特定のものに限定されるものではなく、具体的には例えば空間周波数フィルタや画像フィルタなどが考えられる。
【００６２】
（空間周波数フィルタの適用）
空間周波数フィルタは、「投影」を行った後の二次元分布に対して二次元フーリエ変換を行い、空間周波数（ｋx，ｋy）特性に対してフィルタ処理を行うものである。すなわち、もとの強度分布をＩ(ｘ，ｙ)，その二次元フーリエ変換をＦ(ｋx，ｋy)，空間周波数フィルタをＨ(ｋx，ｋy)とすると、処理後の強度分布Ｉ'(ｘ，ｙ)は数式７で表される。
【数７】
2014130124.tif 000009
【００６３】
空間分解能の向上に有効な空間周波数フィルタとしてはハイパスフィルタが挙げられる。具体的には、急峻なカットオフを用いると処理画像が劣化するため、数式８に示す線形フィルタＨL、或いは、数式９に示す逆ガウス型フィルタＨGなどが具体的には挙げられる。ここで、ｋ＝(ｋx＋ｋy)1/2である。
【数８】
2014130124.tif 000010
【数９】
2014130124.tif 000011
【００６４】
数式８における定数ｋ1及びｋ2は、画像化に用いる空間周波数成分の範囲を決めるものであり、特定の値に限定されるものではなく、フィルタを適用するデータの特性などを考慮して適宜設定される。例えば、テラヘルツ波イメージング結果のサンプリング間隔やフィルタ適用前の解像度などを考慮して適宜設定される。
【００６５】
具体的には例えば、テラヘルツ波イメージング結果がサンプリング間隔０．２ mm にて２５０×２５０点（つまり、計測範囲は５０ mm ×５０ mm）測定したものである場合には、二次元フーリエ変換における空間周波数の最低値は１／５０ mm ＝０．０２ mm-1，最高値は１２５／５０ mm ＝２．５ mm-1 になる。そして、フィルタ適用前の解像度が例えば１．６ mm-1 程度である場合には、この値がｋ1とｋ2との中間付近に来るようにｋ1＝０．４ mm-1，ｋ2＝２．５ mm-1 に設定されることが例として挙げられる。
【００６６】
なお、画像の適切なコントラストを保つため、ｋ＜ｋ1の場合について、ＨL＝０ではなく、ＨLがゼロよりも大きい値（具体的には例えば、ＨL＝０．０５）に設定されても良い。
【００６７】
一方、数式９における定数ａは、フィルタ適用後に残留する平均輝度の割合を示すものであり、特定の値に限定されるものではなく、フィルタを適用するデータの特性など考慮して適宜設定される。例えば、フィルタ適用前の解像度などを考慮して適宜設定される。
【００６８】
具体的には例えば、画像の適切なコントラストを保つことを考慮してａ＝０．９５に設定されることが例として挙げられる（なお、コントラストを保つ必要がない場合は、ａ＝１．０に設定される場合もある）。
【００６９】
また、数式９における定数Δｋは、空間周波数フィルタの幅を示すものであり、特定の値に限定されるものではなく、フィルタを適用するデータの特性などを考慮して適宜設定される。例えば、フィルタ適用前の解像度に相当するｋ＝１．６ mm-1 においてＨG＝０．５になるようにΔｋが設定される、すなわち、０．５＝１−０．９５exp[−(１．６／Δｋ)2]からΔｋ＝２．０ mm-1 に設定されることが例として挙げられる。
【００７０】
そして、上述の実施形態の周波数フィルタの適用結果への空間周波数フィルタの更なる適用として、具体的には、Ｓ３までの処理による計算結果としての平面位置(ｘ，ｙ)別のスペクトル強度Ｓ(ｘ，ｙ)に対し、数式８で表される線形フィルタＨLを適用したり、数式９で表される逆ガウス型フィルタＨGを適用したりする。つまり、面内空間分解能の向上装置であれば、数式８で表される線形フィルタＨLを適用する手段を更に有するようにしたり、数式９で表される逆ガウス型フィルタＨGを適用する手段を更に有するようにしたりし、或いは、面内空間分解能の向上プログラムであれば、数式８で表される線形フィルタＨLを適用する手段としてもコンピュータを機能させるようにしたり、数式９で表される逆ガウス型フィルタＨGを適用する手段としてもコンピュータを機能させるようにしたりする。
【００７１】
（画像フィルタの適用）
空間分解能の向上に有効な画像フィルタとしては原画像の２次元行列とマスクと呼ばれる２次元行列Ｋとの畳み込みを行うものが挙げられる。具体的には、原画像の画素(ｉ，ｊ)における輝度をＩijとした場合、処理後の輝度Ｉ'ijは数式１０で表される。
【数１０】
2014130124.tif 000012
【００７２】
数式１０において、Ｉij;klは画素(ｉ，ｊ)を中心とする、大きさｋ×ｌの副行列であり、マスクと同じ大きさである。一例として、エッジ強調に用いられるラプラシアン(Laplacian)フィルタの場合のマスクＫLは数式１１で与えられる。以下では、「ラプラシアン」はこのフィルタを適用した結果を示す。
【数１１】
2014130124.tif 000013
【００７３】
「高域強調フィルタ」（アンシャープマスク処理）は、原画像からラプラシアンを差し引いた形、即ち数式１２で表される。
【数１２】
2014130124.tif 000014
【００７４】
数式１２におけるαは、原画像に対するラプラシアンの割合を示すパラメータであり、α＝０は処理なし（即ち、原画像のまま）の場合に、α＝１は原画像とラプラシアンとを重畳した場合に、α>>１は原画像の平均輝度を維持したラプラシアンに、それぞれ相当する。
【００７５】
数式１１，１２は左右又は上下方向におけるエッジ強調を行うためのマスクであるが、これに斜め方向におけるエッジ強調を加えたマスクは数式１３で表される。
【数１３】
2014130124.tif 000015
【００７６】
そして、上述の実施形態の周波数フィルタの適用結果への画像フィルタの更なる適用として、具体的には、Ｓ３までの処理による計算結果としての平面位置(ｘ，ｙ)別のスペクトル強度Ｓ(ｘ，ｙ)に対し、数式１２で表される上下左右方向エッジ強調マスクＫUを適用したり、数式１３で表される上下左右斜め方向エッジ強調マスクＫVを適用したりする。つまり、面内空間分解能の向上装置であれば、数式１２で表される上下左右方向エッジ強調マスクＫUを適用する手段を更に有するようにしたり、数式１３で表される上下左右斜め方向エッジ強調マスクＫVを適用する手段を更に有するようにしたりし、或いは、面内空間分解能の向上プログラムであれば、数式１２で表される上下左右方向エッジ強調マスクＫUを適用する手段としてもコンピュータを機能させるようにしたり、数式１３で表される上下左右斜め方向エッジ強調マスクＫVを適用する手段としてもコンピュータを機能させるようにしたりする。
【実施例１】
【００７７】
本発明の面内空間分解能の向上方法によるテラヘルツ波イメージングにおける空間分解能の向上効果を検証した実施例を図８から図１７を用いて説明する。
【００７８】
（１）検証方法
本実施例では、図３に示すテラヘルツ波イメージング装置（具体的には、Picometrix社製・T-Ray4000）を用いると共に、図８に示す光学機器用の解像度ターゲット（テストパターン）を用いて検証を行った。なお、０〜２ THz の周波数帯域で計測を行った。
【００７９】
ここで、既存の検討により、反射波の振幅分布に基づく、全周波数成分を含む場合の空間分解能評価を行い、０．６ mm 程度の空間分解能であることが確認されている（福地哲生 他：「テラヘルツ波を用いた遮熱コーティングの非破壊検査における空間分解能評価」，光応用・視覚，計測合同研究会資料，ＬＡＶ−１２−３／ＩＭ−１２−１１，pp.１１−１６，２０１２年）。
【００８０】
本実施例で用いたテストパターン（図８）は、ソーダ石灰ガラスに縦横方向に５本の線群から成るパターンをクロム蒸着したものであり、米国国立標準技術研究所(ＮＩＳＴ)の規格ＮＢＳ１９６３Ａに準拠している。線群の間隔は最大１．０ mm から最小０．０５５ mm であり、線幅は間隔の１／２である。テストパターン中の各線群に付されている数値は１ mm あたりの線数であり、線間隔 [mm]はその逆数になる。なお、以下においては、各線群に付されている数値を用いて、例えば「線群１．０」や「線群１．１」のように各線群のことを表記する。
【００８１】
ここでは、図４〜７に示す周波数フィルタを用いて反射波の電界データの周波数帯域を限定し、イメージング結果を比較した。
【００８２】
そして、数式６によって、矩形フィルタを用いた場合は、矩形範囲として指定した周波数帯域における反射波のスペクトル強度の積分値を計算し、一方、ガウスフィルタを用いた場合は、反射波のスペクトル強度分布とガウス分布との畳み込み積分値を計算した。
【００８３】
（２）検証結果
まず、本発明の面内空間分解能の向上方法を適用しないで（即ち、周波数フィルタを適用しないで）、テストパターンからの反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータを用いて全周波数帯域におけるスペクトル強度を積算して図９に示す反射波の電界強度分布が得られた。
【００８４】
図９に示す結果から、個別の線として識別できる線群は１．６以下（即ち、線間隔０．６２５ mm 以上）であることが確認された。なお、識別の可否は、線群の各線方向に対して直交する直線上の強度分布に５つの有意な極小値が確認できるか否かで判断した（図１０に一例を示す）。
【００８５】
次に、テストパターンからの反射波の電界データに、帯域０．５〜１．０ THz の矩形フィルタ（図４参照）を適用した上でスペクトル強度を積算して図１１に示す結果が得られ、また、帯域１．０〜１．５ THz の矩形フィルタ（図５参照）を適用した上でスペクトル強度を積算して図１２に示す結果が得られた。
【００８６】
図１１における解像度は図９と大きく変わらないことが確認された。一方、図１２における解像度は図９と比べて明らかに改善されており、縦線は線群２．５（線間隔０．４０ mm）、横線は線群２．０（線間隔０．５０ mm）の場合でも識別可能であることが確認された。なお、横方向と縦方向との分解能が異なるのは、送受信器を搭載した可動ステージ７のたわみが影響しているためと考えられた。
【００８７】
また、テストパターンからの反射波の電界データに、中心周波数０．７５ THz のガウスフィルタ（図６参照）を適用した上でスペクトル強度を積算して図１３に示す結果が得られ、中心周波数１．２５ THz のガウスフィルタ（図７参照）を適用した上でスペクトル強度を積算して図１４に示す結果が得られた。なお、ガウスフィルタの半値全幅は０．３５ THz とした。
【００８８】
図１３における解像度は図９や図１１と大きく変わらないことが確認された。一方、図１４においては、縦線は線群２．５（線間隔０．４０ mm）、横線は線群２．０（線間隔０．５０ mm）の場合でも識別可能であり、図１２における解像度と殆ど同じであることが確認された。
【００８９】
また，図９と図１１〜１４とを比較し、周波数フィルタの適用によって線群の濃淡(コントラスト)が増加していることが確認された。これは、各線群に対して直交する方向の強度分布の変調度が向上していることによると考えられた。図９，１２，１４中の線群１．６（線間隔０．６２５ mm）に対して直交する直線上の強度分布を求めて図１０に示す結果が得られた。
【００９０】
図１０に示す結果から、周波数フィルタ無しの場合でも５つの極小値が確認でき、個別の線が読み取れると判断できるものの、変調度（即ち、極大値と極小値との差を平均値で除算した指標）は４％程度であることが確認された。一方で、１．０〜１．５ THz の矩形フィルタ又は中心周波数１．２５ THz のガウスフィルタを用いた場合の変調度は４０％程度であり、フィルタ無しの場合と比べて１桁向上することが確認された。
【００９１】
また、周波数フィルタの適用による濃淡(コントラスト)の向上を定量的に評価するため、図９，１１〜１４中の線群１．０〜２．０（線間隔０．５〜１．０ mm）について、線群に対して直交する直線上の強度分布から５つの極大値の平均値Ｉmaxと５つの極小値の平均値Ｉminとを求め、Ｃ＝(Ｉmax−Ｉmin)／(Ｉmax＋Ｉmin)で与えられるコントラスト比を算出して図１５に示す結果が得られた。なお、コントラスト比Ｃの値が大きいほど濃淡が明確であり、完全な灰色の場合はコントラスト比Ｃ＝０、完全な白黒の場合はコントラストＣ＝１になる。
【００９２】
図１５に示す結果から、周波数帯域１．０〜１．５ THz の矩形フィルタ又は中心周波数１．２５ THz のガウスフィルタを用いることにより、コントラスト比を４倍以上向上させることが可能であることが確認された。
【００９３】
以上の結果から、テラヘルツ波イメージングにおいて、周波数フィルタを用いて反射波の電界データの周波数帯域を適切に限定することによって空間分解能とコントラスト比との向上が可能であることが確認された。具体的には、本発明を適用することにより、本発明を適用しない場合（即ち、全周波数帯域を用いた場合）には０．６ mm 程度である空間分解能を０．４０ mm まで向上させることが可能であり、また、本発明を適用しない場合と比べてコントラスト比を４倍以上向上させることが可能であることが確認された。なお、空間分解能の向上効果は周波数フィルタの形状が異なっても大きな差違はなく、矩形フィルタでもガウスフィルタでも同様の結果が得られることが確認された。
【００９４】
（３）更なるフィルタ・画像処理の適用
次に、周波数フィルタの適用結果に対して他のフィルタや画像処理を更に適用した場合の空間分解能の向上効果を検証した。
【００９５】
まず、周波数フィルタ（ガウスフィルタ：中心周波数１．２５ THz，半値全幅０．４０ THz）の適用結果に空間周波数フィルタである逆ガウス型フィルタＨGを更に適用して図１６に示す結果が得られた。なお、図１６の右図は、左図(画像)中の破線内の輝度を横方向に積算した結果である。また、逆ガウス型フィルタＨGを表す数式９に関し、ａ＝０．９５，Δｋ＝２．０ mm-1 にされた。
【００９６】
図１６に示す結果から、線群２．５の５本の線を明確に判別することができると共に、コントラスト比が大幅に向上することが確認された。
【００９７】
また、周波数フィルタ（ガウスフィルタ：中心周波数１．２５ THz，半値全幅０．４０ THz）の適用結果に上下左右方向エッジ強調マスクＫUを更に適用して図１７に示す結果が得られた。なお、図１７の右図は、左図(画像)中の破線内の輝度を横方向に積算した結果である。また、上下左右方向エッジ強調マスクＫUを表す数式１２に関し、α＝３．０にされた。
【００９８】
図１７に示す結果から、線群２．５の５本の線を明確に判別することができると共に、コントラスト比が大幅に向上することが確認された。
【００９９】
以上の結果から、周波数フィルタの適用結果に空間周波数フィルタ若しくは画像フィルタ（エッジ強調マスク）を更に適用することにより、コントラスト比が大幅に向上し、空間分解能の更なる向上が可能であることが確認された。
【符号の説明】
【０１００】
１０  面内空間分解能の向上装置
１７  面内空間分解能の向上プログラム
【図１】
2014130124.tif 000022
【図２】
2014130124.tif 000023
【図３】
2014130124.tif 000024
【図４】
2014130124.tif 000025
【図５】
2014130124.tif 000026
【図６】
2014130124.tif 000027
【図７】
2014130124.tif 000028
【図８】
2014130124.tif 000029
【図９】
2014130124.tif 000030
【図１０】
2014130124.tif 000031
【図１１】
2014130124.tif 000032
【図１２】
2014130124.tif 000033
【図１３】
2014130124.tif 000034
【図１４】
2014130124.tif 000035
【図１５】
2014130124.tif 000036
【図１６】
2014130124.tif 000037
【図１７】
2014130124.tif 000038