JP2014130124A - 面内空間分解能の向上方法、向上装置及び向上プログラム - Google Patents

Abstract

【課題】テラヘルツ時間領域分光法による空間計測の分解能を向上させることができるようにする。
【解決手段】計測対象物にテラヘルツパルス波を照射して計測対象物からの反射波の電界の計測を行い（Ｓ１）、計測によって得られた電界データに周波数フィルタを適用して電界データの周波数を限定し（Ｓ２）、周波数が限定された電界データを用いてスペクトル強度を計算する（Ｓ３）ようにした。
【選択図】図１

Description

本発明は、面内空間分解能の向上方法、向上装置及び向上プログラムに関する。さらに詳述すると、本発明は、テラヘルツ時間領域分光イメージングによる断面構造の計測に用いて効果的な深度方向の空間分解能を向上させるための技術に関する。

テラヘルツ波技術の進展により、テラヘルツ波イメージングによる非破壊検査は様々な分野で用いられ始めている。これまでに、遮熱材の欠陥検出，隠匿物検出，薬剤のコーティング厚さ分布計測，木材やコンクリートなどの建材中の欠陥検出，絵画の塗料分析などが報告されている。テラヘルツ波はプラスチック，樹脂，セラミック，木材などを比較的良く透過するため、対象物の断面構造（言い換えると、深度方向の分布）を計測することができ、これを面内方向の二次元分布計測と組み合わせることで三次元計測が可能となる。なお、０．１〜１０ 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年）。したがって、従来のテラヘルツ時間領域分光法は空間分解能が十分に高いとは言えない。

そこで、本発明は、テラヘルツ時間領域分光法による空間計測の分解能を向上させることができる面内空間分解能の向上方法、向上装置及び向上プログラムを提供することを目的とする。

本発明者は、テラヘルツ時間領域分光法による空間計測の分解能の向上についての検討を重ね、反射波の電界データ（具体的には、反射波の電界の時間変化を表す波形）の周波数の範囲を限定することによって空間分解能が変化することを突き止めた。

本発明の面内空間分解能の向上方法は、上記知見に基づくものであり、計測対象物にテラヘルツパルス波を照射して計測対象物からの反射波の電界の計測を行い、計測によって得られた電界データに周波数フィルタを適用して電界データの周波数を限定し、周波数が限定された電界データを用いてスペクトル強度を計算するようにしている。

また、本発明の面内空間分解能の向上装置は、計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置から読み込む手段と、電界データに周波数フィルタを適用して電界データの周波数を限定する手段と、周波数が限定された電界データを用いてスペクトル強度を計算する手段とを有するようにしている。

また、本発明の面内空間分解能の向上プログラムは、計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置から読み込む手段、電界データに周波数フィルタを適用して電界データの周波数を限定する手段、周波数が限定された電界データを用いてスペクトル強度を計算する手段としてコンピュータを機能させるようにしている。

これら本発明の面内空間分解能の向上方法、向上装置及び向上プログラムによると、計測によって得られた反射波の電界データの周波数を適切に限定した上でスペクトル強度を計算するようにしているので、空間分解能の向上に対して不必要なデータが排除され、空間分解能の向上に必要なデータのみによってスペクトル強度が計算される。

また、本発明の面内空間分解能の向上方法、向上装置及び向上プログラムは、周波数フィルタとして数式１に示す矩形フィルタ若しくは数式２に示すガウスフィルタを用いるようにし、数式１について０．５≦ａ≦１．２５且つ０．７５≦ｂ≦１．５（ただし、ａ＜ｂ）、或いは、数式２について０．５≦μ≦１．５且つ０．２５≦σ≦１．０であるようにしても良い。

ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
をそれぞれ表す。

ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
σ：フィルタ関数としての分散値
をそれぞれ表す。

これらの場合には、空間分解能の向上に対して不必要なデータが適切に排除されて空間分解能の向上に必要なデータのみが適切に抽出される。

また、本発明の面内空間分解能の向上方法、向上装置及び向上プログラムは、計算されたスペクトル強度に対して線形フィルタ，逆ガウス型フィルタ，上下左右方向エッジ強調マスク，上下左右斜め方向エッジ強調マスクのうちのいずれか一つを更に適用するようにしても良い。この場合には、コントラスト比を向上させ、空間分解能の更なる向上が図られる。

本発明の面内空間分解能の向上方法、向上装置及び向上プログラムによれば、空間分解能の向上に対して不必要なデータを排除し、空間分解能の向上に必要なデータのみによってスペクトル強度を計算することができるので、テラヘルツ波イメージングの空間分解性能を高め、テラヘルツパルス波を用いたイメージングによる非破壊検査技術としての有用性、延いては信頼性の向上を図ることが可能になる。具体的には、検査対象物の断面構造（即ち、深度方向の分布）の空間分解能を向上させて検査処理に対して有用な情報を提供することを可能にし、検査精度の向上が可能になる。

本発明の面内空間分解能の向上方法の実施形態の一例を説明するフローチャートである。 実施形態の面内空間分解能の向上方法を面内空間分解能の向上プログラムを用いて実施する場合の当該プログラムによって実現される面内空間分解能の向上装置の機能ブロック図である。 実施形態のテラヘルツ波イメージング装置の概略構成図である。 本発明で用いられ得る矩形フィルタ（周波数帯域：０．５〜１．０ THz）を示す図である。 本発明で用いられ得る矩形フィルタ（周波数帯域：１．０〜１．５ THz）を示す図である。 本発明で用いられ得るガウスフィルタ（中心周波数：０．７５ THz）を示す図である。 本発明で用いられ得るガウスフィルタ（中心周波数：１．２５ THz）を示す図である。 実施例１で用いた解像度ターゲット（テストパターン）を示す図である。 実施例１におけるテストパターンからの反射波の全周波数帯域におけるスペクトル強度を積算した結果を示す図である。 実施例１における線群に対して垂直な直線上の強度分布の一例を示す図である。 実施例１における周波数帯域０．５〜１．０ THz の矩形フィルタを適用した場合のスペクトル強度の計算結果を示す図である。 実施例１における周波数帯域１．０〜１．５ THz の矩形フィルタを適用した場合のスペクトル強度の計算結果を示す図である。 実施例１における中心周波数０．７５ THz のガウスフィルタを適用した場合のスペクトル強度の計算結果を示す図である。 実施例１における中心周波数１．２５ THz のガウスフィルタを適用した場合のスペクトル強度の計算結果を示す図である。 実施例１における線群に対して垂直な直線上の強度分布から求めたコントラスト比を計算した結果を示す図である。 実施例１における周波数フィルタの適用結果に空間周波数フィルタ（逆ガウス型フィルタ）を更に適用した場合のスペクトル強度の計算結果を示す図である。 実施例１における周波数フィルタの適用結果に上下左右方向エッジ強調マスクを更に適用した場合のスペクトル強度の計算結果を示す図である。

以下、本発明の構成を図面に示す実施の形態の一例に基づいて詳細に説明する。

図１から図７に、本発明の面内空間分解能の向上方法、向上装置及び向上プログラムの実施形態の一例を示す。

本実施形態の面内空間分解能の向上方法は、図１に示すように、計測対象物にテラヘルツパルス波を照射して計測対象物からの反射波の電界の計測を行い（Ｓ１）、計測によって得られた電界データに周波数フィルタを適用して電界データの周波数を限定し（Ｓ２）、周波数が限定された電界データを用いてスペクトル強度を計算する（Ｓ３）ようにしている。

また、本実施形態の面内空間分解能の向上装置は、計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置としてのデータサーバ(１６)から読み込む手段(１１ａ)と、電界データに周波数フィルタを適用して電界データの周波数を限定する手段(１１ｂ)と、周波数が限定された電界データを用いてスペクトル強度を計算する手段(１１ｃ)とを有する。

さらに、本実施形態の面内空間分解能の向上プログラムは、計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置としてのデータサーバ(１６)から読み込む手段(１１ａ)、電界データに周波数フィルタを適用して電界データの周波数を限定する手段(１１ｂ)、周波数が限定された電界データを用いてスペクトル強度を計算する手段(１１ｃ)としてコンピュータを機能させる。

そして、面内空間分解能の向上方法の実行にあたっては、まず、計測対象物にテラヘルツパルス波を照射することによる計測対象物からの反射波の電界の計測が行われる（Ｓ１）。

本発明で用いられ得るテラヘルツ波イメージング装置の一例を図３に示す。このテラヘルツ波イメージング装置は、レーザ２から出射された超短パルスレーザ光をハーフミラー８によって二つに分岐する。そして、分岐したうちの一方を送信用光伝導アンテナ３に入射して幅約１ ps のテラヘルツ波パルスを発生させ、このテラヘルツ波はレンズ４によって計測対象物１の表面上に集束され、計測対象物１から反射されたテラヘルツ波はレンズ４及び偏光子９を介して受信用光伝導アンテナ５に導かれて超短パルスレーザ光を入射した時点における電界を計測するようにしている。

テラヘルツ波の電界の時間変化を得るには、ハーフミラー８によって分岐したうちの一方のレーザ光であって送信用光伝導アンテナ３へ導入するレーザ光に対し、前記分岐したうちの他方のレーザ光であって受信用光伝導アンテナ５へ導入するレーザ光の時間遅延を走査する。この時間遅延の走査は可動ミラー６を用いて行う。本実施形態の装置では、可動ミラー６は１秒間に１００回走査されるため、サンプリング周波数は１００ Hz である。

受信用光伝導アンテナ５によって受信された反射波のテラヘルツ波信号は、必要に応じてアンプ（図示していない）によって増幅された後にＡ／Ｄコンバータでデジタル化され、保存される。可動ミラー６の位置を変えながら上記動作を繰り返すことによって、テラヘルツ波の時間波形として取得される。

また、送受信器（具体的にはアンテナ３，５など）は一体型のヘッドに搭載されて可動ステージ７に取り付けられている。そして、ヘッドを平面方向に走査することによってテラヘルツ波を照射する計測対象物１上の平面位置（ｘ，ｙ）を変えることができる。個々の平面位置（ｘ，ｙ）において可動ミラー６の移動による時間遅延の走査を行い、テラヘルツ波の時間波形を取得し、最終的には時間波形の２次元分布が得られる。

なお、本実施形態では図３に示すテラヘルツ波イメージング装置を用いるようにしているが、本発明で用いられ得る計測装置は図３に示すものには限られない。すなわち、テラヘルツ時間領域分光法において用いられ得るものであって、具体的には、テラヘルツ波を出射して反射波の電界（テラヘルツパルス波の振動電場の強さ）を計測し得るものであれば、どのような計測装置であっても構わない。

次に、Ｓ１の処理によって計測され取得された反射波の電界データに対して周波数フィルタを適用することによる電界データの周波数の限定が行われる（Ｓ２）。

ここで、本発明の面内空間分解能の向上方法におけるＳ２以降の処理は本発明の面内空間分解能の向上装置によって実行され得る。

そして、本発明の面内空間分解能の向上方法におけるＳ２以降の処理及びこれら処理を実行する面内空間分解能の向上装置は、本発明の面内空間分解能の向上プログラムをコンピュータ上で実行することによっても実現され得る。本明細書では、面内空間分解能の向上プログラムをコンピュータ上で実行することによってＳ２以降の処理を実行する面内空間分解能の向上装置が実現されると共に面内空間分解能の向上方法におけるＳ２以降の処理が実行される場合を説明する。

面内空間分解能の向上プログラム１７を実行するためのコンピュータ１０（本実施形態では、面内空間分解能の向上装置１０でもある）の全体構成を図２に示す。このコンピュータ１０（面内空間分解能の向上装置１０）は、制御部１１、記憶部１２、入力部１３、表示部１４及びメモリ１５を備え相互にバス等の信号回線によって接続されている。また、コンピュータ１０には記憶装置としてのデータサーバ１６がバス等の信号回線によって接続されており、その信号回線を介してデータや制御指令等の信号の送受信（即ち出入力）が相互に行われる。

制御部１１は記憶部１２に記憶されている面内空間分解能の向上プログラム１７によってコンピュータ１０全体の制御並びに面内空間分解能の向上に係る演算を行うものであり、例えばＣＰＵ（中央演算処理装置）である。

記憶部１２は少なくともデータやプログラムを記憶可能な装置であり、例えばハードディスクである。

メモリ１５は制御部１１が種々の制御や演算を実行する際の作業領域であるメモリ空間となるものであり、例えばＲＡＭ(Ｒandom Ａccess Ｍemory の略)である。

入力部１３は少なくとも作業者の命令を制御部１１に与えるためのインターフェイスであり、例えばキーボードである。

表示部１４は制御部１１の制御によって文字や図形等の描画・表示を行うものであり、例えばディスプレイである。

そして、本実施形態では、上述のＳ１の処理において計測され取得された反射波の電界データが反射波電界データベース１８としてデータサーバ１６に格納(保存)される。具体的には、原点(任意)からの横方向平面位置ｘ mm、縦方向平面位置ｙ mm における周波数ｆ THz 別の反射波の電界(電場)の強さＦ(ｘ，ｙ，ｆ) a.u. の値が反射波電界データベース１８に記録される。なお、データサーバ１６に格納される反射波の電界の強さＦ(ｘ，ｙ，ｆ) のデータは全ての周波数(ｆ)成分を含むものである。また、Ｓ１の処理において計測されたデータは、例えば、計測装置から適当な記憶媒体に保存されて作業者によってデータサーバ１６内の反射波電界データベース１８に記録・蓄積されるようにしても良いし、計測装置から通信手段（無線・有線）を介してデータサーバ１６内の反射波電界データベース１８に自動的に記録・蓄積されるようにしても良い。

そして、コンピュータ１０（本実施形態では、面内空間分解能の向上装置１０でもある）の制御部１１には、面内空間分解能の向上プログラム１７を実行することにより、Ｓ１の処理において計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置としてのデータサーバ１６から読み込む処理を行うデータ読込部１１ａと、電界データに周波数フィルタを適用して電界データの周波数を限定する処理を行う周波数フィルタ部１１ｂと、周波数が限定された電界データを用いてスペクトル強度を計算する処理を行うスペクトル強度計算部１１ｃとが構成される。

面内空間分解能の向上プログラム１７が実行されることによる具体的な処理としては、まず、コンピュータ１０（面内空間分解能の向上装置１０）の制御部１１に構成されたデータ読込部１１ａが反射波の電界データの読み込みを行う（Ｓ２−１）。

具体的には、データ読込部１１ａは、Ｓ１の処理において計測され取得されてデータサーバ１６に格納されている反射波電界データベース１８に記録されている反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータをデータサーバ１６から読み込む。

そして、データ読込部１１ａは、読み込んだ反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータをメモリ１５に記憶させる。

続いて、制御部１１の周波数フィルタ部１１ｂが、反射波の電界データに対して周波数フィルタを適用することによる電界データの周波数の限定を行う（Ｓ２−２）。

具体的には、周波数フィルタ部１１ｂは、Ｓ２−１の処理においてメモリ１５に記憶された反射波の電界の強さＦ(ｘ，ｙ，ｆ) のデータをメモリ１５から読み込む。

そして、周波数フィルタ部１１ｂは、周波数フィルタを適用して反射波の電界データの周波数帯域を限定する。本発明における周波数フィルタとしては、例えば、周波数軸(周波数の並び)に対して関数値が矩形の分布型を有する矩形フィルタ（数式３）やガウス分布型を有するガウスフィルタ（数式４）が用いられ得る。

ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
をそれぞれ表す。

ここに、ｗ(ｆ)：フィルタ関数，
ｆ：周波数 [THz]，
μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
σ：フィルタ関数としての分散値 [THz]
をそれぞれ表す。

なお、矩形フィルタにおける、フィルタ関数によって限定する周波数帯域の下限周波数ａ THz の値は０．５〜１．２５ THz の範囲で適宜選択されると共に上限周波数ｂ THz の値は０．７５〜１．５ THz の範囲で適宜選択され、０．７５≦ａ≦１．２５且つ１．０≦ｂ≦１．５が好ましく、１．００≦ａ≦１．２５且つ１．２５≦ｂ≦１．５がより一層好ましい（ただし、いずれについてもａ＜ｂ）。また、ｂ−ａ≧０．２５が好ましく、ｂ−ａ＝０．５がより一層好ましい。以上を踏まえて具体的には例えば、ａ＝０．５且つｂ＝０．７５，ａ＝０．７５且つｂ＝１．０が選択され得て、ａ＝０．５且つｂ＝１．０が好ましく（図４）、ａ＝１．０且つｂ＝１．２５，ａ＝１．２５且つｂ＝１．５がより一層好ましく、ａ＝１．０且つｂ＝１．５が最も好ましい（図５）。

また、ガウスフィルタにおける、フィルタ関数によって着目する周波数帯域の中心周波数σ THz の値は、０．５〜１．５ THz の範囲で適宜選択され、０．７５≦σ≦１．５が好ましく（σ＝０．７５：図６）、１．０≦σ≦１．５がより一層好ましく、σ＝１．２５が更に好ましい（図７）。

また、分散値σの値は、０．２５〜１．０ THzの範囲で適宜選択され、０．２５≦σ≦０．７５が好ましく、０．４≦σ≦０．６がより一層好ましく、σ＝０．５が更に好ましい。

そして、周波数フィルタ部１１ｂは、周波数フィルタを適用した後の反射波の電界の強さＦ_flt(ｘ，ｙ，ｆ)を数式５によって計算する。

ここに、Ｆ_flt(ｘ，ｙ，ｆ)：周波数フィルタ適用後の反射波の電界の強さ [a.u.]，
ｗ(ｆ)：フィルタ関数，
Ｆ(ｘ，ｙ，ｆ)：反射波の電界の強さ(計測値) [a.u]，
ｘ：原点(任意)からの横方向平面位置 [mm]，
ｙ：原点(任意)からの縦方向平面位置 [mm]，
ｆ：周波数 [THz]
をそれぞれ表す。

そして、周波数フィルタ部１１ｂは、計算した周波数フィルタ適用後の反射波の電界の強さＦ_flt(ｘ，ｙ，ｆ)の値をメモリ１５に記憶させる。

次に、Ｓ２−２の処理によって計算された周波数フィルタ適用後の周波数が限定された反射波の電界の強さＦ_flt(ｘ，ｙ，ｆ)のデータを用いてのスペクトル強度の計算を行う（Ｓ３）。

具体的には、スペクトル強度計算部１１ｃは、Ｓ２−２の処理においてメモリ１５に記憶された周波数フィルタ適用後の反射波の電界の強さＦ_flt(ｘ，ｙ，ｆ)の値をメモリ１５から読み込み、平面位置(ｘ，ｙ)別の全周波数成分を積算したスペクトル強度Ｓ(ｘ，ｙ)を、数式６によって計算する。

ここに、Ｓ(ｘ，ｙ)：スペクトル強度 [a.u]，
Ｆ_flt(ｘ，ｙ，ｆ)：周波数フィルタ適用後の反射波の電界の強さ [a.u.]，
ｘ：原点(任意)からの横方向平面位置 [mm]，
ｙ：原点(任意)からの縦方向平面位置 [mm]，
ｆ：周波数 [THz]
をそれぞれ表す。

そして、スペクトル計算部１１ｂは、Ｓ３までの処理による計算結果として、計算した平面位置(ｘ，ｙ)別のスペクトル強度Ｓ(ｘ，ｙ)の値を、表示部１４に表示したり、例えば記憶部１２やデータサーバ１６に計算結果データファイルとして保存したりする。

そして、制御部１１は、面内空間分解能の向上の処理を終了する（ＥＮＤ）。

以上の構成を有する本発明の面内空間分解能の向上方法、向上装置及び向上プログラムによれば、計測によって得られた反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータの周波数ｆを適切に限定した上でスペクトル強度Ｓ(ｘ，ｙ)を計算するようにしているので、空間分解能の向上に対して不必要なデータを排除し、空間分解能の向上に必要なデータのみによってスペクトル強度を計算することができ、テラヘルツ波イメージングの空間分解性能を高め、テラヘルツパルス波を用いたイメージングによる非破壊検査技術としての有用性、延いては信頼性の向上を図ることが可能になる。

なお、上述の形態は本発明の好適な形態の一例ではあるが本発明の実施の形態がこれに限定されるものではなく、本発明の要旨を逸脱しない範囲において種々変形実施可能である。例えば、上述の実施形態では反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータが蓄積される記憶手段をデータサーバ１６としているが、記憶部１２でも良いし、他の記憶装置を用いるようにしても良い。また、反射波の電界の強さＦ(ｘ，ｙ，ｆ)のデータを記憶装置に一旦蓄積することなく、計測装置からコンピュータ１０（面内空間分解能の向上装置１０）に直接入力し、入力されたデータから処理を順次行うようにしても良い。

また、上述の実施形態における周波数フィルタの適用結果に、他のフィルタ・他の画像処理を更に適用するようにしても良い。ここで、更に適用するフィルタ・画像処理は、特定のものに限定されるものではなく、具体的には例えば空間周波数フィルタや画像フィルタなどが考えられる。

（空間周波数フィルタの適用）
空間周波数フィルタは、「投影」を行った後の二次元分布に対して二次元フーリエ変換を行い、空間周波数（ｋ_x，ｋ_y）特性に対してフィルタ処理を行うものである。すなわち、もとの強度分布をＩ(ｘ，ｙ)，その二次元フーリエ変換をＦ(ｋ_x，ｋ_y)，空間周波数フィルタをＨ(ｋ_x，ｋ_y)とすると、処理後の強度分布Ｉ'(ｘ，ｙ)は数式７で表される。

空間分解能の向上に有効な空間周波数フィルタとしてはハイパスフィルタが挙げられる。具体的には、急峻なカットオフを用いると処理画像が劣化するため、数式８に示す線形フィルタＨ_L、或いは、数式９に示す逆ガウス型フィルタＨ_Gなどが具体的には挙げられる。ここで、ｋ＝(ｋ_x＋ｋ_y)^1/2である。

数式８における定数ｋ₁及びｋ₂は、画像化に用いる空間周波数成分の範囲を決めるものであり、特定の値に限定されるものではなく、フィルタを適用するデータの特性などを考慮して適宜設定される。例えば、テラヘルツ波イメージング結果のサンプリング間隔やフィルタ適用前の解像度などを考慮して適宜設定される。

具体的には例えば、テラヘルツ波イメージング結果がサンプリング間隔０．２ mm にて２５０×２５０点（つまり、計測範囲は５０ mm ×５０ mm）測定したものである場合には、二次元フーリエ変換における空間周波数の最低値は１／５０ mm ＝０．０２ mm^-1，最高値は１２５／５０ mm ＝２．５ mm^-1 になる。そして、フィルタ適用前の解像度が例えば１．６ mm^-1 程度である場合には、この値がｋ₁とｋ₂との中間付近に来るようにｋ₁＝０．４ mm^-1，ｋ₂＝２．５ mm^-1 に設定されることが例として挙げられる。

なお、画像の適切なコントラストを保つため、ｋ＜ｋ₁の場合について、Ｈ_L＝０ではなく、Ｈ_Lがゼロよりも大きい値（具体的には例えば、Ｈ_L＝０．０５）に設定されても良い。

一方、数式９における定数ａは、フィルタ適用後に残留する平均輝度の割合を示すものであり、特定の値に限定されるものではなく、フィルタを適用するデータの特性など考慮して適宜設定される。例えば、フィルタ適用前の解像度などを考慮して適宜設定される。

具体的には例えば、画像の適切なコントラストを保つことを考慮してａ＝０．９５に設定されることが例として挙げられる（なお、コントラストを保つ必要がない場合は、ａ＝１．０に設定される場合もある）。

また、数式９における定数Δｋは、空間周波数フィルタの幅を示すものであり、特定の値に限定されるものではなく、フィルタを適用するデータの特性などを考慮して適宜設定される。例えば、フィルタ適用前の解像度に相当するｋ＝１．６ mm^-1 においてＨ_G＝０．５になるようにΔｋが設定される、すなわち、０．５＝１−０．９５exp[−(１．６／Δｋ)²]からΔｋ＝２．０ mm^-1 に設定されることが例として挙げられる。

そして、上述の実施形態の周波数フィルタの適用結果への空間周波数フィルタの更なる適用として、具体的には、Ｓ３までの処理による計算結果としての平面位置(ｘ，ｙ)別のスペクトル強度Ｓ(ｘ，ｙ)に対し、数式８で表される線形フィルタＨ_Lを適用したり、数式９で表される逆ガウス型フィルタＨ_Gを適用したりする。つまり、面内空間分解能の向上装置であれば、数式８で表される線形フィルタＨ_Lを適用する手段を更に有するようにしたり、数式９で表される逆ガウス型フィルタＨ_Gを適用する手段を更に有するようにしたりし、或いは、面内空間分解能の向上プログラムであれば、数式８で表される線形フィルタＨ_Lを適用する手段としてもコンピュータを機能させるようにしたり、数式９で表される逆ガウス型フィルタＨ_Gを適用する手段としてもコンピュータを機能させるようにしたりする。

（画像フィルタの適用）
空間分解能の向上に有効な画像フィルタとしては原画像の２次元行列とマスクと呼ばれる２次元行列Ｋとの畳み込みを行うものが挙げられる。具体的には、原画像の画素(ｉ，ｊ)における輝度をＩ_ijとした場合、処理後の輝度Ｉ'_ijは数式１０で表される。

数式１０において、Ｉ_ij;klは画素(ｉ，ｊ)を中心とする、大きさｋ×ｌの副行列であり、マスクと同じ大きさである。一例として、エッジ強調に用いられるラプラシアン(Laplacian)フィルタの場合のマスクＫ_Lは数式１１で与えられる。以下では、「ラプラシアン」はこのフィルタを適用した結果を示す。

「高域強調フィルタ」（アンシャープマスク処理）は、原画像からラプラシアンを差し引いた形、即ち数式１２で表される。

数式１２におけるαは、原画像に対するラプラシアンの割合を示すパラメータであり、α＝０は処理なし（即ち、原画像のまま）の場合に、α＝１は原画像とラプラシアンとを重畳した場合に、α>>１は原画像の平均輝度を維持したラプラシアンに、それぞれ相当する。

数式１１，１２は左右又は上下方向におけるエッジ強調を行うためのマスクであるが、これに斜め方向におけるエッジ強調を加えたマスクは数式１３で表される。

そして、上述の実施形態の周波数フィルタの適用結果への画像フィルタの更なる適用として、具体的には、Ｓ３までの処理による計算結果としての平面位置(ｘ，ｙ)別のスペクトル強度Ｓ(ｘ，ｙ)に対し、数式１２で表される上下左右方向エッジ強調マスクＫ_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を表す数式１２に関し、α＝３．０にされた。

図１７に示す結果から、線群２．５の５本の線を明確に判別することができると共に、コントラスト比が大幅に向上することが確認された。

以上の結果から、周波数フィルタの適用結果に空間周波数フィルタ若しくは画像フィルタ（エッジ強調マスク）を更に適用することにより、コントラスト比が大幅に向上し、空間分解能の更なる向上が可能であることが確認された。

１０ 面内空間分解能の向上装置
１７ 面内空間分解能の向上プログラム

Claims (9)

1. 計測対象物にテラヘルツパルス波を照射して計測対象物からの反射波の電界の計測を行い、前記計測によって得られた電界データに周波数フィルタを適用して前記電界データの周波数を限定し、前記周波数が限定された電界データを用いてスペクトル強度を計算することを特徴とする面内空間分解能の向上方法。
2. 前記周波数フィルタとして以下の数式１に示す矩形フィルタ
   

   

   ここに、ｗ(ｆ)：フィルタ関数，
   ｆ：周波数 [THz]，
   ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
   ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
   をそれぞれ表す。
   若しくは以下の数式２に示すガウスフィルタ
   

   

   ここに、ｗ(ｆ)：フィルタ関数，
   ｆ：周波数 [THz]，
   μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
   σ：フィルタ関数としての分散値
   をそれぞれ表す。
   を用いるようにし、数式１について０．５≦ａ≦１．２５且つ０．７５≦ｂ≦１．５（ただし、ａ＜ｂ）、或いは、数式２について０．５≦μ≦１．５且つ０．２５≦σ≦１．０であることを特徴とする請求項１記載の面内空間分解能の向上方法。
3. 計算された前記スペクトル強度に対して線形フィルタ，逆ガウス型フィルタ，上下左右方向エッジ強調マスク，上下左右斜め方向エッジ強調マスクのうちのいずれか一つを更に適用することを特徴とする請求項１記載の面内空間分解能の向上方法。
4. 計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置から読み込む手段と、前記電界データに周波数フィルタを適用して前記電界データの周波数を限定する手段と、前記周波数が限定された電界データを用いてスペクトル強度を計算する手段とを有することを特徴とする面内空間分解能の向上装置。
5. 前記周波数フィルタとして以下の数式３に示す矩形フィルタ
   

   

   ここに、ｗ(ｆ)：フィルタ関数，
   ｆ：周波数 [THz]，
   ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
   ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
   をそれぞれ表す。
   若しくは以下の数式４に示すガウスフィルタ
   

   

   ここに、ｗ(ｆ)：フィルタ関数，
   ｆ：周波数 [THz]，
   μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
   σ：フィルタ関数としての分散値
   をそれぞれ表す。
   を用いるようにし、数式３について０．５≦ａ≦１．２５且つ０．７５≦ｂ≦１．５（ただし、ａ＜ｂ）、或いは、数式４について０．５≦μ≦１．５且つ０．２５≦σ≦１．０であることを特徴とする請求項４記載の面内空間分解能の向上装置。
6. 計算された前記スペクトル強度に対して線形フィルタ，逆ガウス型フィルタ，上下左右方向エッジ強調マスク，上下左右斜め方向エッジ強調マスクのうちのいずれか一つを適用する手段を更に有することを特徴とする請求項４記載の面内空間分解能の向上装置。
7. 計測対象物にテラヘルツパルス波を照射して計測された計測対象物からの反射波の電界データを記憶装置から読み込む手段、前記電界データに周波数フィルタを適用して前記電界データの周波数を限定する手段、前記周波数が限定された電界データを用いてスペクトル強度を計算する手段としてコンピュータを機能させることを特徴とする面内空間分解能の向上プログラム。
8. 前記周波数フィルタとして以下の数式５に示す矩形フィルタ
   

   

   ここに、ｗ(ｆ)：フィルタ関数，
   ｆ：周波数 [THz]，
   ａ：フィルタ関数によって限定する周波数帯域の下限周波数 [THz]，
   ｂ：フィルタ関数によって限定する周波数帯域の上限周波数 [THz]
   をそれぞれ表す。
   若しくは以下の数式６に示すガウスフィルタ
   

   

   ここに、ｗ(ｆ)：フィルタ関数，
   ｆ：周波数 [THz]，
   μ：フィルタ関数によって着目する周波数帯域の中心周波数 [THz]，
   σ：フィルタ関数としての分散値
   をそれぞれ表す。
   を用いるようにし、数式５について０．５≦ａ≦１．２５且つ０．７５≦ｂ≦１．５（ただし、ａ＜ｂ）、或いは、数式６について０．５≦μ≦１．５且つ０．２５≦σ≦１．０であることを特徴とする請求項７記載の面内空間分解能の向上プログラム。
9. 計算された前記スペクトル強度に対して線形フィルタ，逆ガウス型フィルタ，上下左右方向エッジ強調マスク，上下左右斜め方向エッジ強調マスクのうちのいずれか一つを適用する手段としてもコンピュータを機能させることを特徴とする請求項７記載の面内空間分解能の向上プログラム。

Images