JP5808836B2

JP5808836B2 - 波面を推定する方法

Info

Publication number: JP5808836B2
Application number: JP2014086774A
Authority: JP
Inventors: 大久保　彰律; 彰律大久保; バレットハリソン; フレンリドラース; カウッチルカ
Original assignee: University of Arizona
Current assignee: University of Arizona
Priority date: 2011-10-14
Filing date: 2014-04-18
Publication date: 2015-11-10
Anticipated expiration: 2032-10-12
Also published as: US9182289B2; JP2014132280A; JP5564545B2; US20130092816A1; JP2013088432A

Description

本出願は、参照により開示全体が本明細書に組み込まれている、２０１１年１０月１４日に出願した仮出願第６１／５４７５８３号の優先権を主張するものである。

本発明は、波面を推定する方法に関する。

波面測定は通常、光学的な表面の質を検査し、光学エレメントの光学性能を評価するために用いられる。波面測定は、高精度光学系のアライメントにも用いられる。波面は、すべての光線の位相が同じである点の集まり（線、または３次元の波の伝播では面）である。最も簡単な形態の波面は平面波の波面であり、この場合、光線は互いに平行であり、平らな波面がセンサに当たる。レンズやミラーなどの光学エレメントを通る波面の伝播では、一般に、レンズの厚さ、表面形状の不完全さ、屈折率のばらつき、および他の要因により波面の形状が変化する。波面形状の望まれない変化は収差として知られている。そのため、波面プロファイルを知り波面収差を補正することは、光学エレメントを設計する際や、新たに設計した光学系の性能を評価する際に極めて重要である。例えば、結像系一式を組み立てる前に、その系に含まれる各光学単位（単位レンズ）の性能を確認することが必要である。単位レンズまたは単レンズ自体はそれぞれある種の収差を有し得るので、結像レンズの質を高精度に制御することが必要である。大きな波面収差を測定する一応用例は、レンズやミラーなどの非球面光学エレメントの曲率中心の検査である。

光ビームの波面の質を測定する従来の方法では干渉測定波面センサを用いる。干渉測定波面センサでは、光源ビームの小部分を空間フィルタリングすることによって球面参照波を生成し、続いてそれを元の波面と合成してインターフェログラムを生成する。当技術分野では周知のように、インターフェログラムの干渉縞を解析して光ビームの質を評価することができる。しかし、元のビームを分割しそれを参照ビームと再合成すると、例えば、光路に誤差があったり、光学構成エレメントのアライメントが不適切な場合には、システム収差が発生しやすい。従来の波面測定法に代わるものとして、元のビームを分割し再合成する必要のないシャック−ハルトマン（Ｓｈａｃｋ−Ｈａｒｔｍａｎ）波面センサなどの非干渉測定波面センサ（ＮＩＷＦＳ）を用いる方法がある。

シャック−ハルトマン波面センサ（ＳＨＷＦＳ）などのＮＩＷＦＳは、干渉測定センサよりもダイナミックレンジが広い。ＳＨＷＦＳセンサの基本的かつ通常用いられている一構成は、レンズレットアレイと、レンズレットアレイの後側焦点面に配置される光検出器（典型的にはＣＣＤカメラ）とからなる。シャック−ハルトマン波面センサは、２次元（２Ｄ）レンズレットアレイを用いることによって測定対象の入射ビームの波面を複数の小ビームに分割する。このレンズレットアレイのレンズレットはそれぞれ、光検出器の表面に別々の独立した焦点（スポット）を生成する。各スポットの重心位置は、参照ビームと歪んだビームの間の波面収差を示す距離および方向によって変位する。したがって、ＳＨＷＦＳによる波面測定は、参照（平面）波面に対する収差を含む波面の局所的な勾配を推定することに基づいている。一般に、波面推定の手順は、位相が波面の局所的な勾配の集合として示されるか、開口全体にわたって求められるある種のモーダル関数の係数の集合として示されるかに応じて、ゾーナル式およびモーダル式のいずれかに分類し得る。後者では、焦点スポットの変位をゼルニケ多項式の形で表し得る。

ＳＨＷＦＳは干渉計測波面センサに対して以下のいくつかの利点を有する。ＳＨＷＦＳは干渉計測センサよりもダイナミックレンジが広い。入射する放射はコヒーレントでなくてよい。ＳＨＷＦＳは単一像から波面情報のすべてを取得し得るので、露光時間を短くすることができ、そのため、振動の影響を受けにくい。より重要には、ＳＨＷＦＳでは放射照度分布および位相分布がともに得られる。

図１に、ＳＨＷＦＳを利用する波面測定装置１００の構成例を示す。図１に示すように、波面測定装置１００は、典型的には、光源１０２、ニュートラルデンシティ（ＮＤ）フィルタ１０４、ビームエキスパンダ１０６、投影レンズ１０８、検査光学素子（検査エレメントまたは試料）１１０、レンズ群１１２、レンズレットアレイ（マイクロレンズアレイ）１１４、光検出器１１６（ＣＣＤセンサ）、およびデータ解析部１１８を含み得る。このように配置された波面測定装置１００は、試料（検査光学素子１１０）がそれを通過する光の波面に及ぼす影響を特徴付けるように構成されている。

図２に、ＳＨＷＦＳの光学構成をより詳細に示す。図２では、レンズレットアレイ１１４を通過した光が入射するセンサアレイ２１０が光検出器１１６に含まれている。センサアレイ２１０上の焦点スポット２５０、２５１、２５２、２５３、および２５４の場所は、入射する波面１５０の局所的な傾きによって決まる。波面１５０の局所的な傾きは、図１の検査光学素子１１０による収差によって生じる。この局所的な傾きは、焦点スポットの場所の変化によって計算し得る。レンズレットアレイ１１４の各レンズレットから得られる局所的な傾きの情報を用いることによって波面を再構成することができる。

波面の偏差量がＳＨＷＦＳのダイナミックレンジ内であれば、センサアレイ２１０上の各焦点スポットの位置を別々に検出し、対応するレンズレットに割り当てることができ、波面プロファイルを容易に求めることができる。しかし、図３に示すように波面の偏差がＳＨＷＦＳのダイナミックレンジを越える場合、センサアレイＳＨＷＦＳはもはや波面を解析することができない。具体的には、図３では、レンズレットアレイ１１４は、焦点スポット３５０、３５１、３５２、３５３、および３５４を形成する。しかし、波面１６０が大きな収差を含むために、焦点スポット３５０はセンサアレイ２１０の表面の範囲からはみ出てしまい、焦点スポット３５３と３５４は互いに交差してしまう（スポット交差）。つまり、焦点スポット３５３および３５４を形成するビームは互いに重なり合い、交差した後でセンサアレイ２１０上のそれぞれの場所に結像する。図３に示す状況では、焦点スポット３５０はセンサアレイ２１０の範囲からはみ出しているので、焦点スポット３５０を検出することができず、焦点スポット３５３および３５４はスポット交差のために明確に割り当てることができない。

図４に、収差を含む波面を検出したセンサアレイ２１０からの出力データの例を示す。図４では、波面の外側領域の焦点スポットを示す関心領域Ｓが拡大されて、焦点スポットが互いに交差したり、センサアレイ２１０の有効区域からはみ出したりすることの影響がよりよく示されている。

ＳＨＷＦＳのダイナミックレンジを拡大する技術や波面収差を解析的に補償する技術が以前に提案されており、その開発が継続中である。こうした比較的注目されている技術をいくつか以下にまとめる。
（１）ヌルレンズ
ヌルレンズは、予想される波面収差を補償または除去するように特に設計された１組の光学素子を含む。ヌルレンズ技術は、（製造誤差がないという前提で）検査光学素子の波面収差を完全に補償し得るので、レンズレットアレイ上で波面偏差を効果的に打ち消すことができる。しかし、この技術では、正確な測定を行うためには極めて正確なヌルレンズを製作することが必要である。したがって、ヌルレンズの製作コストが途方もなく高額になり得る。さらに、このようなヌルレンズは波面収差を予想し得る特定の検査光学素子に対して設計されるので、形状または特性が異なる検査光学素子によって形成される他の波面にはこの技術を適用し得ないことがある。ヌルレンズ技術の例が、Ｚｍｅｋに付与された米国特許第５，２３３，１７４号に記載されている。この特許を参照により本明細書に組み込む。
（２）推定技術
ヌルレンズ技術の代わりに、収差を含む波面を測定する波面推定技術が提案されており、典型的には補正アルゴリズムを用いて収差を補正または補償する。推定技術の一例が、参照文献１：ＭｉｃｈａｅｌＣ．Ｒｏｇｇｅｍａｎｎ、ＴｉｍｏｔｈｙＪ．Ｓｃｈｕｌｚ、ＣｈｅｅＷ．Ｎｇａｉ、およびＪａｓｏｎＴ．Ｋｒａｆｔ、「ＪｏｉｎｔｐｒｏｃｅｓｓｉｎｇｏｆＨａｒｔｍａｎｎｓｅｎｓｏｒａｎｄｃｏｎｖｅｎｔｉｏｎａｌｉｍａｇｅｍｅａｓｕｒｅｍｅｎｔｓｔｏｅｓｔｉｍａｔｅｌａｒｇｅａｂｅｒｒａｔｉｏｎｓ：ｔｈｅｏｒｙａｎｄｅｘｐｅｒｉｍｅｎｔａｌｒｅｓｕｌｔｓ」、Ａｐｐｌ．Ｏｐｔ．３８、２２４９〜２２５５（１９９９）に開示されている。

別の波面推定技術では、最尤推定法（ＭａｘｉｍｕｍＬｉｋｅｌｉｈｏｏｄＥｓｔｉｍａｔｉｏｎ、ＭＬＥ）を用いて波面を再構築する。ＭＬＥ技術の例が、参照文献２：ＨａｒｒｉｓｏｎＨ．Ｂａｒｒｅｔｔ、ＣｈｒｉｓｔｏｐｈｅｒＤａｉｎｔｙ、およびＤａｖｉｄＬａｒａ、「Ｍａｘｉｍｕｍ−ｌｉｋｅｌｉｈｏｏｄｍｅｔｈｏｄｓｉｎｗａｖｅｆｒｏｎｔｓｅｎｓｉｎｇ：ｓｔｏｃｈａｓｔｉｃｍｏｄｅｌｓａｎｄｌｉｋｅｌｉｈｏｏｄｆｕｎｃｔｉｏｎｓ」、Ｊ．Ｏｐｔ．Ｓｏｃ．Ａｍ．Ａ．、２４、３９１〜４１４（２００７）に開示されている。
（３）ステッチング
米国特許出願公開第２００９／０２８４７５３号には、数学的方法を用いて一連の波面測定値を合わせて「ステッチ」する技術が記載されている。各測定毎に、異なる焦点、波面の傾き、または参照収差をダイナミックレンジ限界開口と合わせて用いる。この技術は、センサのダイナミックレンジを効果的に拡大し得ると言われているが、一連の波面測定を行い、数学的方法を用いて得られた測定値を合わせてステッチすることは計算集約的プロセスであることが、当業者には容易に理解されよう。

米国特許第５，２３３，１７４号米国特許出願公開第２００９／０２８４７５３号

ＭｉｃｈａｅｌＣ．Ｒｏｇｇｅｍａｎｎ、ＴｉｍｏｔｈｙＪ．Ｓｃｈｕｌｚ、ＣｈｅｅＷ．Ｎｇａｉ、およびＪａｓｏｎＴ．Ｋｒａｆｔ、「ＪｏｉｎｔｐｒｏｃｅｓｓｉｎｇｏｆＨａｒｔｍａｎｎｓｅｎｓｏｒａｎｄｃｏｎｖｅｎｔｉｏｎａｌｉｍａｇｅｍｅａｓｕｒｅｍｅｎｔｓｔｏｅｓｔｉｍａｔｅｌａｒｇｅａｂｅｒｒａｔｉｏｎｓ：ｔｈｅｏｒｙａｎｄｅｘｐｅｒｉｍｅｎｔａｌｒｅｓｕｌｔｓ」、Ａｐｐｌ．Ｏｐｔ．３８、２２４９〜２２５５（１９９９）ＨａｒｒｉｓｏｎＨ．Ｂａｒｒｅｔｔ、ＣｈｒｉｓｔｏｐｈｅｒＤａｉｎｔｙ、およびＤａｖｉｄＬａｒａ、「Ｍａｘｉｍｕｍ−ｌｉｋｅｌｉｈｏｏｄｍｅｔｈｏｄｓｉｎｗａｖｅｆｒｏｎｔｓｅｎｓｉｎｇ：ｓｔｏｃｈａｓｔｉｃｍｏｄｅｌｓａｎｄｌｉｋｅｌｉｈｏｏｄｆｕｎｃｔｉｏｎｓ」、Ｊ．Ｏｐｔ．Ｓｏｃ．Ａｍ．Ａ．、２４、３９１〜４１４（２００７）Ｋｒｉｅｇら、「Ｉｎｖｅｒｓｅｐｒｏｐａｇａｔｉｏｎａｌｇｏｒｉｔｈｍｆｏｒａｎｇｓｔｒｏｍａｃｃｕｒａｃｙｉｎｔｅｒｆｅｒｏｍｅｔｅｒ」、ＩｎｔｅｒｆｅｒｏｍｅｔｒｙＸＩＩ：ＴｅｃｈｎｉｑｕｅｓａｎｄＡｎａｌｙｓｉｓ、５５３１巻、ＰｒｏｃｅｅｄｉｎｇｓｏｆｔｈｅＳｏｃｉｅｔｙｏｆＰｈｏｔｏ−ＯｐｔｉｃａｌＩｎｓｔｒｕｍｅｎｔａｔｉｏｎＥｎｇｉｎｅｅｒｓ（ＳＰＩＥ）、Ｃｒｅａｔｈら監修、２００４ＪａｅＭｙｕｎｇ、「ＴｕｔｏｒｉａｌｏｎＭａｘｉｍｕｍＬｉｋｅｌｉｈｏｏｄＥｓｔｉｍａｔｉｏｎ」、ＪｏｕｒｎａｌｏｆＭａｔｈｅｍａｔｉｃａｌＰｓｙｃｈｏｌｏｇｙ、４７巻、９０〜１００（２００３）Ｌ．Ｒ．Ｆｕｒｅｎｌｉｄ、Ｊ．Ｙ．Ｈｅｓｔｅｍａｎ、およびＨ．Ｈ．Ｂａｒｒｅｔｔ、「Ｆａｓｔｍａｘｉｍｕｍ−ｌｉｋｅｌｉｈｏｏｄｅｓｔｉｍａｔｉｏｎｍｅｔｈｏｄｓｆｏｒｓｃｉｎｔｉｌｌａｔｉｏｎｃａｍｅｒａｓａｎｄｏｔｈｅｒｏｐｔｉｃａｌｓｅｎｓｏｒｓ」、Ｐｒｏｃ．ＳＰＩＥ、６７０７、２００７

当技術分野の上記状況に鑑みて、本発明者らは本発明で、波面パラメータを推定する新規技術を開発した。この技術の利点は、中でもとりわけ、波面誤差を補正して打ち消すと同時に、波面パラメータの測定および解析に必要とされる時間を短縮し得ることである。

本発明の一態様は、ゼルニケ多項式で偏差が表される波面を最尤推定法を用いることによって推定する方法であって、レンズレットアレイおよび検出器を含む波面センサを用いることによって得られるデータを準備する工程と、ゼルニケ多項式の互いに異なる２つのゼルニケ係数により前記検出器上の強度分布をコスト関数として求め、前記２つのゼルニケ係数をそれぞれ座標軸とする座標面において当該コスト関数を表したコスト関数分布を生成する工程と、前記座標軸に対する前記コスト関数分布の対称性を調べるための区域を選択する工程と、前記区域において前記コスト関数分布の前記座標軸に対する対称性を調べる工程と、同時に最適化すべきゼルニケ多項式の係数の組合せとして、前記対称性を調べた結果に基づいて、前記座標軸に対して非対称なゼルニケ係数の組み合わせを選択する工程と、前記コスト関数の最小値を探索することによって前記組合せに従うゼルニケ多項式の係数について最大尤度最適化を実行する工程と、を含む。

本発明のさらなる特徴は、例示実施形態の以下の説明を添付の図面と併せ読めば、当業者には明らかである。

シャック−ハルトマン波面センサ（ＳＨＷＦＳ）を利用する波面測定装置の光学構成を示す図。収差を含まない波面が入射した状態のＳＨＷＦＳの光学構成を示す図。大きな収差を含む波面が入射したときのＳＨＷＦＳの動作を示す図。収差を含む波面が入射したときにＳＨＷＦＳのセンサアレイから得られる出力データを示す図。本発明の一態様による、波面パラメータを推定する波面測定装置の例を示す図。波面補正エレメントが間に介在する状態での検査エレメントから光検出器までの光線追跡図を示す図。検査エレメントがある場合とない場合の光検出器によって検出される像の例を示す図。検査エレメントがある場合とない場合の光検出器によって検出される像の例を示す図。平行平面板が間に介在する状態での検査エレメントから光検出器までの光線追跡図を示す図。収縮グリッド探索アルゴリズムを用いる最適化プロセスをグラフィカルに示す図。パラメータ｛ａ１，ａ２｝の探索空間におけるコスト関数面プロファイルを示す図。パラメータ｛ａ３，ａ４｝の探索空間におけるコスト関数面プロファイルを示す図。パラメータ｛ａ１，ａ３｝の組合せのコスト関数面を示す図。本発明の第２態様による波面パラメータを推定する波面測定装置の例を示す図。データ解析部の概略図を示す図。光学構成を示す図。光学構成を示す図。光学構成を示す図。光学構成を示す図。シミュレーションされた検出器データの例を示す図。ゼルニケ係数に対するコスト関数の変化を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。様々なパラメータの組合せの２次元コスト関数面を示す図。最適化手順のフローチャートを示す。対象なコスト関数面プロファイルの９点を示す図。対象なコスト関数面プロファイルの９点を示す図。非対象なコスト関数面プロファイルの９点を示す図。非対象なコスト関数面プロファイルの９点を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。収縮グリッド探索アルゴリズムに基づいてゼルニケ係数を推定するフロープロセスを示す図。最尤（ＭＬ）推定プロセスに用いられるゼルニケ係数の推定値とターゲット値の差を示す図。係数｛Ｃ７，Ｃ８｝の組合せのコスト関数面を示す図。係数｛Ｃ３，Ｃ８｝の組合せのコスト関数面を示す図。パラメータ｛Ｃ３，Ｃ８｝のフィッティング面を示す図。フィッティング方法を用いる収縮グリッド探索アルゴリズムに基づいてゼルニケ係数を推定する改善された最適化プロセスを示す図。フィッティング方法を用いるＭＬ推定中の、推定されたゼルニケ係数とターゲット係数の差を示す図。

本発明による例示実施形態を以下に説明する。本明細書では、波面は、この波面を構成するすべての光線がほぼ同じ位相を有する点の集まり（線、または３次元の波伝播では面）とみなし得る。「波面」は一般に光学的波面の形で記述される。当業者には理解されるように、「光学の」、「光の」などの用語は、明確にことわらない限り可視光に限定されるものではない。具体的には、ほぼ同じ波面特性を有する波を生成し得る電磁スペクトルは、その波が赤外、可視、紫外、またはより大きいエネルギー範囲（例えば、Ｘ線およびガンマ線）の電磁スペクトルのいずれであろうが、かなりの範囲のエネルギーにわたっていることを当業者は理解しているはずである。さらに、波面収差について言及する場合、「収差」という用語は、予想、推定、または確定された参照値からの任意の偏差であり、これにはいくつかの（最小許容差を超える）偏差の和も含まれることを理解されたい。したがって、波面収差の「補正」とは、予想、推定、または確定された参照値からの任意の偏差を減少させる、最小限に抑える、または補償することを指し得る。さらに、波面「パラメータ」について言及する場合、波面の傾き、形状、位相、表面、場所、速さなどを特徴付ける任意の値または値の群をパラメータとみなし得る。

本発明の一態様は、波面形状や波面収差などの光学系の波面パラメータを推定する装置および方法に関するものである。図５に、本発明の一実施形態による、波面パラメータを推定する装置５００の例示セットアップを示す。

図５では、装置５００は、光源５０１、波面センサ５１０、およびデータ解析部５１８を備える。光源５０１は、コリメート光源または点光源とし得る単色光源として実装し得る。光学検査エレメント（すなわち、被検査光学エレメント）５０３は、光源５０１と波面センサ５１０の間に光路に沿って配置される。必要に応じて、光源５０１と検査エレメント５０３の間にレンズ５０９およびビームエキスパンダ５０２を配置し得る。ＮＤフィルタ（図示せず）もこのセットアップに適用し得る。

波面センサ５１０は、好ましくは非干渉測定波面センサ（ＮＩＷＦＳ）として実装される。レンズレットアレイ５０５および光検出器５０６がこの非干渉測定波面センサを形成するように配置される。本実施形態では、ＮＩＷＦＳの例としてシャック−ハルトマン波面センサ（ＳＨＷＦＳ）を用いる。ただし、本実施形態はＳＨＷＦＳに限定されるものではない。多くのタイプの非干渉測定波面センサがあり、それぞれがそれ自体の特定の配置を有する傾向がある。添付の特許請求の範囲に含まれる限り、任意の機器を使用して波面パラメータを測定し得る。本発明の一態様では、負の平凹レンズなどの補正エレメント５０４（波面補正エレメント）がＳＨＷＦＳの光源側に導入される。具体的には、補正エレメント５０４の例として負の平凹レンズが検査エレメント５０３とレンズレットアレイ５０５の間に配置される。光源５０１によって生成される光は、レンズ５０９およびビームエキスパンダ５０２を通って進み、検査エレメント５０３に入射するビームを形成する。このビームは、検査エレメント５０３を通過した後で補正エレメント５０４に進み、レンズレットアレイ５０５を通過して、ＣＣセンサアレイなどの光検出器５０６に至る。光検出器５０６に入射する光の強度が検出され、データ解析部５１８が光検出器５０６からの信号を電子的に処理（解析）し得る。データ解析部５１８は、コンピュータ制御下でデータを取得し記憶する。データ解析部５１８は、特殊なソフトウェアを用いてデータフィッティングを行い、様々な計算を行い、データをグラフィカルな形式で表示する。

補正エレメント５０４（負の平凹レンズなど）は、未補正の波面をＳＨＷＦＳのダイナミックレンジよりも大きくしてしまういくつかの収差項を部分的または全体的に補正することによってデータ取得ステップおよび実際の収差の推定ステップをより正確かつより実用的にする単一または一群の高精度光学構成エレメントを含む。補正エレメントは光源５０１とレンズレットアレイ５０５の間であればどこにでも配置し得るが、レンズレットアレイと検査エレメントの間に補正エレメントを配置することが好ましい。つまり、補正エレメント５０４は、レンズレットアレイ５０５の光源側に配置される。単一の平凹レンズの代わりに、検査エレメントの特徴（特性）に基づいて、単一の凹レンズ、単一の凸レンズ、単一の平凸レンズ、またはこれらの組合せを補正エレメント５０４として使用してもよい。実際、補正エレメント５０４は、補正エレメント５０４の表面がレンズレットアレイ５０５と接触し得るようにレンズレットアレイ５０５と一体化することさえ可能である。

図６に、検査エレメント５０３の出射面から像面ＩＰ（すなわち、光検出器５０６のセンサアレイが配置される面）までの光線追跡図の例を示す。検査エレメント５０３の最終面から発する光を、像面に向かって凹である軌跡を有する波面として、図の最も左側に示す。ある実験配置では、この波面の近軸焦点は、検査エレメント５０３の最終面から右に１１０．６４ｍｍ離れていた。波面収差を補正するか、または少なくとも部分的に補償するために、補正エレメント５０４（例えば、ｆが−７５ｍｍで、φが５０．８ｍｍの平凹レンズ）を検査エレメント５０３の後ろに配置した。２５×２５個の焦点距離４０ｍｍのレンズレットをピッチ１ｍｍで配列したアレイ（図示せず）を補正エレメント５０４の後ろ（図６では右）に配置した。実装を容易にするためにレンズレットアレイを負レンズ（補正エレメント）の平坦面上に製作することもできるし、補正レンズに隣接してその検出器側に配置することもできる。３６×３６ｍｍ^２の正方形センサアレイを含む光検出器（図５の５０６）をレンズレットアレイの焦点面に配置する。レンズレットアレイの焦点面で、光検出器のセンサアレイの各エレメントは、レンズレットアレイ内の対応するレンズレットによって形成される光ビームを受け取る。原理的には、センサアレイ全体からのデータを用いてレンズレットアレイに入射する光ビームの波面を電子的に再構築し得る。

検出器アレイによって検出される波面は多くの異なる方法で表現し得る。波面は、例えば、任意の平面または曲面における光の場の位相に関する１組のゼルニケ係数の形で表現し得る。グリッド上のサンプルとして波面の位相の形で波面を表現することもできる。あるいは、波面は、電場の位相に関する解析的表現によって特徴付けることができる。実際、各単一レンズレットに関連するセンサアレイから得られるデータの副集合をいくつか用いて、波面の局所的な傾きを計算し得る。次いで、これらの局所的な傾きを用いて光学系の他の面での波面パラメータを推定する。天文計測に応用する適応光学素子の場合、波面は大気中を進むが、局所的な傾きの観察値を求め、大気状態を知ることにより、波面パラメータを簡単かつ正確に推定し得る。しかし、天文測定に応用する適応光学素子の場合と異なり、製造または実験的なセットアップにおける波面は大きな収差を含んでおり、大気中で歪んだ典型的な波から大きく逸脱している。そのため、高精度光学素子では、局所的な傾きの値を観察し、大気状態を知るだけでは、波面パラメータを簡単かつ正確に推定し得ない。

具体的には、図５および図６の配置では、光は、検査エレメント５０３から空気（または別の既知の均質な媒質）中を補正エレメント５０４まで伝播し、次いでレンズレットアレイ５０５に至ると仮定する。図６に示す構成では、検査エレメント５０３からの波面が完全な球面波であれば、補正エレメント５０４からコリメートされた光として出射されるように補正エレメント５０４により光が曲げられる。補正エレメント５０４の後ろにはレンズレットアレイ５０５が配置される（ただし、直接接触する必要はない）。レンズレットアレイ５０５の各レンズレットは局所的な波面をサンプリングし、光検出器５０６の対応する領域に像を生成する。

検査エレメント５０３から出射する波がほぼ球面波の場合、光検出器５０６上の像は、図７Ａに示すもののように見えるはずである。一方、検査エレメント５０３から出射する波に少なくともいくらか収差が含まれている場合、光検出器５０６上の像は図７Ｂに示すもののように見えることがある。図６ならびに図７Ａおよび図７Ｂを見れば、レンズレットアレイ５０５のレンズレットはそれぞれ光検出器５０６上にスポットを生成し得ることが理解されよう。このスポットパターンおよびスポット自体の特徴、例えば勾配や重心位置は、検査エレメント５０３から出射する波面の収差によって、かつ、補正エレメント５０４およびレンズレットアレイ５０５の光学的かつ物理的な特性によって決まる。しかし、本発明のこの態様では、検査エレメント５０３から出射する波面の収差を大きく減少させるように補正エレメント５０４が特に設計（構成）されており、そのため、波面推定をより容易に実施し得る。

より具体的には、補正エレメント５０４の有用性は、図８に示す光線追跡を調べ、それを図６と比較することによってよりよく理解し得る。図８には、（図６に示す）補正エレメント５０４が取り除かれ、その代わりに光学的に透明な平行平面板５２４が置かれた場合を示す。レンズレットアレイ５０５および光検出器５０６は図６に示すものから変えていない。図８から、補正エレメント５０４がないと、個々のレンズレットからの光線が光検出器５０６の中心領域に集束されることが理解されよう。さらに、図８の場合、図６の場合に反して、個々のレンズレットからの焦点スポットが交差して隣接するレンズレットの検出器領域に入りやすい（スポット交差）。

言い換えると、波面に大きな収差が含まれる状態では、シャック−ハルトマン波面センサにおける焦点スポットは、ほぼ規則的なグリッドに整列される単純なスポットにならない。そうではなく、スポットが複雑に分散し、おそらくは隣接するレンズレットによって生成されるスポットと重なり合い干渉し合う。原理上は、これらの複雑なパターンは測定中の波面についての情報を多く含んでいるが、実際には、これらのパターンを正確に妥当な時間内で解釈するのは極めて難しいことがある。

レンズレットアレイの光源側に負レンズなどの適切に設計された補正エレメントを追加すると、収差を含む波における非合焦や低次の球面収差からの大きな寄与を相殺または減少させるように作用し、残余の収差は、検出器の表面全体にわたって、ただしその範囲内できれいな輪郭の焦点スポットを得るのに十分に小さいことが図６と図８の比較から理解されよう。言い換えると、負の凹レンズなどの形態の補正エレメント５０４は、検査エレメント５０３によって生じる波面収差を少なくとも部分的に補正し、波面センサのダイナミックレンジを広げることができ、それによって、波面推定装置内で検査物体のセンサ側のほぼすべての面における波面パラメータの迅速な推定が容易になる。

図６に示すセットアップにおける光線追跡は、検査エレメント５０３から出射する波面のパラメータを推定する反復手順における前方モデルとして効果的に使用し得る。この推定プロセスは、波面パラメータの最尤（ＭＬ）推定によって実施し得る。最尤推定問題は一般に、パラメータ空間全体にわたる探索として定式化され、式（１）によって与えられる。

・・・（１）
ここで、θは複数のパラメータからなるベクトルであり、ｇは複数の観察値（データ）からなるベクトル、

は複数の推定モデルパラメータからなるベクトルである。解くべき一般問題は、１組の観察値ｇが与えられると仮定して、これらの観察データを生成する最も高い確率（最大尤度）を有するパラメータθの組はどれかという問いの形で提示される。

本例では、ＭＬ推定アルゴリズムに入力されるのは、光検出器で観察される１組のスポット特徴である。各スポット毎に考え得る特徴には、重心、尖度、歪度、分散、ピッチなどが含まれる。この１組の特徴は、ノイズの多い測定（観察）データを含み、ｇと表す。ｇは光検出器によって収集される像から計算し得ることに留意されたい。特徴は、図６に示すように光線追跡図を計算し、これらを用いて検出器上の放射照度パターンを計算することによって、波面パラメータベクトルθからも計算し得る。この場合、検出器で得られる波面パラメータの１組の特徴を

と表す。この表記を用いてＭＬ推定を式（２）によって数学的に表現し得る。

・・・（２）
ここで、
はθのＭＬ推定値を示し、

で表される量は特徴ベクトルｇを測定する確率である。特徴ベクトルｇは、「実際の」特徴ベクトルが
であるという仮定で条件付けられる。

上記の最大化問題を解くのに、共役勾配探索、最急降下探索、疑似焼きなまし法、収縮グリッド探索など多くの効率的なアルゴリズムが利用可能である。例として、収縮グリッド探索の詳細な説明を、この探索法が前方伝播モデルに適用される場合について以下に示す。収縮グリッド探索アルゴリズムは、当業者なら容易に実施し得るものである。可能な応用例として、補正エレメント５０４の前の様々な面（検査中の光学エレメントの射出面など）における波面パラメータを、対象とする面までレンズレットアレイにおける場を後方伝播させることによって推定することができる。つまり、本発明のこの態様では、検査エレメント５０３によって生じる波面収差は、像面から補正エレメント５０４と検査エレメント５０３の間の面まで波面パラメータを後方伝播させることによって推定し得る。後方伝播は数値計算によって実施される。このような数値計算は光線追跡またはフーリエ光学技術を用いて実施し得る。ただし、後方伝播を実施する際、後方伝播アルゴリズムの正確さおよび処理時間を慎重に解析すべきである。処理時間の検討については、式（５）〜（１５）を参照して以下で論じる。後方伝播アルゴリズムの例が参照文献４：Ｋｒｉｅｇら、「Ｉｎｖｅｒｓｅｐｒｏｐａｇａｔｉｏｎａｌｇｏｒｉｔｈｍｆｏｒａｎｇｓｔｒｏｍａｃｃｕｒａｃｙｉｎｔｅｒｆｅｒｏｍｅｔｅｒ」、ＩｎｔｅｒｆｅｒｏｍｅｔｒｙＸＩＩ：ＴｅｃｈｎｉｑｕｅｓａｎｄＡｎａｌｙｓｉｓ、５５３１巻、ＰｒｏｃｅｅｄｉｎｇｓｏｆｔｈｅＳｏｃｉｅｔｙｏｆＰｈｏｔｏ−ＯｐｔｉｃａｌＩｎｓｔｒｕｍｅｎｔａｔｉｏｎＥｎｇｉｎｅｅｒｓ（ＳＰＩＥ）、Ｃｒｅａｔｈら監修、２００４に記載されている。

上記で論じたように、非干渉測定センサ、特にシャック−ハルトマン波面センサは一般に、大きな合焦／非合焦項（大きな収差を含む波面）を伴う光学エレメントの検査には効果的ではないことに留意されたい。というのは、波面偏差がこうしたセンサのダイナミックレンジを越えやすいからである。上記の説明では、適切に設計された補正エレメントを検査エレメントと従来のＳＨＷＦＳの間に配置して使用することによって、解析中の波面の大きな合焦／非合焦または他の大振幅収差を補正エレメントが効果的に、少なくとも部分的に補償することができ、残った小振幅の収差を解明するだけであることが示された。さらに、補正エレメントをレンズレットアレイに接触させることは厳密には必要とされないので、補正エレメントとレンズレットアレイの間の隙間は、特定の推定問題に応じて、例えば、解析中の波面の収差の予想されるレベルに基づいて変更することができる。

一実施形態では、波面をゼルニケ多項式で表すことができ、検査エレメントの出射面の検出器側の所定の面からの前方伝播モデルを用いてゼルニケ係数を最適化することによって波面パラメータを推定し得る。所定の面の光源側の他の面において波面を決定しなければならない場合、前記所定の面からの後方伝播によって波面パラメータを推定し得る。

波面パラメータの推定には最尤推定（ＭＬＥ）技術を用いる。ＭＬＥは周知の推定方法であるが、高次元問題に実際に用いるには演算集約的過ぎると考えられている。上記で論じたように、参照文献２にＭＬＥ技術の一例が詳細に論じられている。別の例が参照文献３：ＪａｅＭｙｕｎｇ、「ＴｕｔｏｒｉａｌｏｎＭａｘｉｍｕｍＬｉｋｅｌｉｈｏｏｄＥｓｔｉｍａｔｉｏｎ」、ＪｏｕｒｎａｌｏｆＭａｔｈｅｍａｔｉｃａｌＰｓｙｃｈｏｌｏｇｙ、４７巻、９０〜１００（２００３）に提供されている。以下で説明するように、本願では、収縮グリッド法を用いることによってＭＬＥを波面推定で実用的なものにすることを紹介する。

一般に、ＭＬＥ技術を用いて推定を行うには、推定するパラメータが与えられると仮定して、データの確率を知ることが必要とされる。確率分布関数（ＰＤＦ）はｐｒ（ｇ｜θ）で示され、ここで、ｇはＭ×１ベクトルによって特徴付けられる観察データベクトルであり、θはＰ×１ベクトルによって特徴付けられるパラメータベクトルである。

θの尤度は式（３）によって表される。

Ｌ（θ｜ｇ）＝ｐｒ（ｇ｜θ）
・・・（３）
θの最尤（ＭＬ）推定値は、観察データｇが与えられると仮定して、式（４）によって定義される。

・・・（４）

実際には、解析的にＭＬ推定値を得ることは通常可能でなく、数値的にＭＬ推定を実行し得る。尤度の数値最適化は、非線形最適化法において反復的に実施される。適切なパラメータ（本例ではゼルニケ係数）を見つけるために、パラメータの値をそれらが収束する方向に変化させながら計算を繰り返し実行し得る。推測対象の

の計算を前方モデルと呼ぶことがある。この数値最適化では、いくつかのモデルは、それらのアルゴリズムにおいて導関数を用いる。例えば、共役勾配探索法または最急降下探索法が最適化アルゴリズムによく用いられる。

導関数を用いる計算にかかる演算時間は極めて長くなることがあり、リソース集約的であるため、導関数が不要なアルゴリズムがより効率的であり、そのため好ましい。

導関数が不要なアルゴリズムがいくつか提案されている。収縮グリッド探索アルゴリズムは、並列処理用に単純かつ効率的なアルゴリズムとして提案されたものである。例えば、参照文献５：Ｌ．Ｒ．Ｆｕｒｅｎｌｉｄ、Ｊ．Ｙ．Ｈｅｓｔｅｍａｎ、およびＨ．Ｈ．Ｂａｒｒｅｔｔ、「Ｆａｓｔｍａｘｉｍｕｍ−ｌｉｋｅｌｉｈｏｏｄｅｓｔｉｍａｔｉｏｎｍｅｔｈｏｄｓｆｏｒｓｃｉｎｔｉｌｌａｔｉｏｎｃａｍｅｒａｓａｎｄｏｔｈｅｒｏｐｔｉｃａｌｓｅｎｓｏｒｓ」、Ｐｒｏｃ．ＳＰＩＥ、６７０７、２００７を参照されたい。収縮グリッド探索アルゴリズムによれば、所望の精度レベルによって決まる所定数の反復において、関数の最大（または最小）値が特定される。この探索法における最適化プロセスが、参照文献５に詳細に説明されているが、ここでは便宜上、それを以下のようにまとめる。このアルゴリズムの実施は、複数のパラメータθからなるＭ×１ベクトルを含むＮ次元空間の探索に適用される。各パラメータ毎に、
１．物理的に妥当なパラメータ空間の周囲の領域を、データベクトル値にかかわらず、無視し得るほど小さい（ほぼゼロに等しい）尤度をもたらす面積（パラメータの数、したがって探索の次元数によっては体積または超体積）で埋める。こうすると、ステップ３ａ（下記）での検査場所に、許容し得る推定値ドメインの外側を、条件検査（ソフトウェアにおける「ｉｆ／ｔｈｅｎ文」）を用いずに、最大尤度として戻される可能性なしに探査させることができる。（あるいは、実行速度をわずかに犠牲にして、物理的に妥当なパラメータ空間内の場所のみが検査されるように条件文を用いる。）
２．較正データ（または予備知識）のみによって決まる尤度表現におけるすべての項を事前に演算する。

３．ｉについて、１から所望の精度で推定値を解くのに必要とされる反復回数まで、
ａ．被検査パラメータ領域を覆う正方形（立方体または超立方体）の中心として定義される検査場所における尤度または対数尤度を演算し、
ｂ．次の反復の検査領域の中心として、尤度または対数尤度が最大になる（または最小限に抑える場合には、最小になる）場所を選択し、
ｃ．各次元における収縮率（例えば、１／２）によって検査領域のサイズを小さくする。

図９に、収縮グリッド探索法のグラフィカルな例を、ＩからＶＩまでの漸進ステップで示す。ここでの目標は、様々なゼルニケ多項式の係数のどの集合が最も大きな（最大）尤度を有しているかを判断することである。図９では、ステップ（Ｉ）で、区域９２０（探索空間）内で、１６個の点９１０からなるグリッドが所定の探索範囲ε内に設定されている。点９１０の数および間隔は、所望の精度レベルに基づいて定義し得る。初期探索範囲ε内に点９１０を定義した後で、これら１６個の点において尤度を演算する。尤度が最も大きい（最大である）適切な点を選択する。ステップ（ＩＩ）で、ステップ（Ｉ）で見つけた最大値を中心として新たなグリッドを定義する。新たなグリッド内の点の間隔（収縮率）は、前のグリッドの間隔の半分とする。図９に示すように、尤度が最も大きい（最大である）適切な点を含む探索範囲を、段階（Ｉ）から段階（ＶＩ）まで、収縮率によって徐々に狭める。

ここで、この収縮グリッド探索法にかかる演算時間を考える。そのために、以下の表記を定義する。
Ｐ：最適化されるパラメータの数
ε_ｐ：最適化を行う初期探索範囲、ｐ＝｛１，２，．．．，Ｐ｝
ｍ_ｐ：各パラメータの更新（尤度分布の計算）後の縮小比（縮小率）、ｐ＝｛１，２，．．．，Ｐ｝
Ｎ_ｐ：尤度分布の計算の反復回数（図９ではＮ_ｐ＝６）
ｌ_ｐ：探索区域の１辺上の尤度分布計算点数、ｐ＝｛１，２，．．．，Ｐ｝（図９ではｌ＝４）
導関数不要の最適化方法では、演算時間の大部分が前方モデルの計算に費やされる。したがって、ＭＬＥのための全演算時間ｔ_ＭＬＥは、式（５）によって以下のように推定される。

ｔ_ＭＬＥ＝ｔ_ｆＮ_{ｔｏｔａｌ}
・・・（５）
ここで、Ｎ_{ｔｏｔａｌ}は前方モデル計算数の合計であり、ｔ_ｆは１回の前方モデル計算にかかる演算時間である。

すべてのパラメータを同時に最適化する場合、前方モデルの計算点の合計数Ｎ_{ｔｏｔａｌ}は式（６）によって与えられる。

Ｎ_{ｔｏｔａｌ}＝Ｎ_ｍａｘｌ^ｐ
ただし（Ｎ_{ｔｏｔａｌ}＝Ｎ_ｍａｘｌ_１ｌ_２．．．ｌ_ｉ．．．ｌ_ｐ）
・・・（６）
ここで、すべてのｌ_ｐがｌに等しく、すべてのＮ_ｐがＮ_ｍａｘに等しいと考える。ｌ＝４、Ｐ＝９、Ｎ_ｍａｘ＝６、およびｔ_ｆ＝１０秒の場合、ｔ_ＭＬＥは式（７）によって計算される。

ｔ_ＭＬＥ＝ｔ_ｆＮ_ｍａｘｌ^ｐ＝１０×６×４^９／３６００（時間）≒４３６９（時間）
・・・（７）
この時間（約４３６９時間）は、ほとんどの応用例で非現実的であり得る。したがって、収縮グリッド探索技術の場合、探査空間の次元が小さいことが好ましい。

例えば、図９では、Ｎ_{ｔｏｔａｌ}＝６×４^２＝９６点である。同時に最適化しない場合のＮ_{ｔｏｔａｌ}を考えると、パラメータの数Ｐを群の数Ｇによって分割し得る。各群のパラメータ数はｋ_ｇ＝｛ｋ_１，ｋ_２，．．．，ｋ_Ｇ｝であり、各群毎の反復回数はＮ_ｇ＝｛Ｎ_１，Ｎ_２，．．．，Ｎ_Ｇ｝である。こうすると、Ｎ_{ｔｏｔａｌ}は式（８）になる。

・・・（８）
ここでｋ_０＝０
ｌ_ｇ＝｛ｌ_１，ｌ_２，．．．，ｌ_Ｇ｝の場合、式（８）は式（９）になる。

・・・（９）
各群毎の反復回数がＮ_ｍａｘである場合、式（８）は式（１０）になる。

・・・（１０）
各パラメータ毎の順次最適化では、それぞれの最適化における反復回数はＮ_ｐになる。式（１０）を用いて、Ｎ_{ｔｏｔａｌ}は式（１１）になる。

・・・（１１）
ここで、

・・・（１２）
各パラメータにおいてＮ_ｐおよびｌ_ｐが一定の場合、式（８）を簡略化することができ、それによって式（１３）が得られる。
Ｎ_{ｔｏｔａｌ}＝Ｎ_ｍａｘＰｌ
・・・（１３）
「ｎ」個の並列プロセッサを使用し得る場合、数ｎ_ＣＰＵを考えることができる。そうすると、式（９）を用いて、各プロセッサ毎の前方モデル計算の合計数は式（１４）になる。

・・・（１４）
式（１４）は、各プロセッサ毎の前方モデル計算数が減少し得ることを示している。例えば、９個のパラメータ（すなわちＰ＝９）は３つの群に分けられ、各群は３つのパラメータを有する。

この場合、Ｇ＝３およびｋ_ｇ＝｛３，３，３｝である。式（５）および式（１４）を用いて、ｌ＝４、Ｐ＝９、Ｎ_ｍａｘ＝６、およびｔ_ｆ＝１０秒の場合、ｔ_ＭＬＥは式（１５）になる。

・・・（１５）
したがって、並列処理を用いる収縮グリッド探索アルゴリズムは、探索空間の次元が小さい場合には効率的な探索方法となり得る。例えば、この最適化方法において、４つのパラメータを以下のように最適化し得る。

まず、パラメータａの集合により一群のパラメータＡが形成されるとし、Ａ＝｛ａ１，ａ２，ａ３，ａ４｝とする。

次いで、パラメータの集合Ａを２つの群Ａ１＝｛ａ１，ａ２｝とＡ２＝｛ａ３，ａ４｝に分割する。Ａ１群のパラメータおよびＡ２群のパラメータは同時に最適化し得る。パラメータの初期値は、例えば、｛ａ１，ａ２，ａ３，ａ４｝＝｛０，０，０．５，０｝と設定される。真値はそれぞれ｛ａ１，ａ２，ａ３，ａ４｝＝｛０，０，０，０｝である。

正規化コスト関数は、対数尤度関数の否定に対応する。コスト関数の最小値を与えるゼルニケ係数は最大尤度推定値である。コスト関数は、光検出器の各画素における実験データと計算データの差の二乗に関係している。本実施形態では、図示の便宜上、光検出器が含むセンサアレイは、１２８０×１０２４画素の行列を形成すると仮定する。｛ａ１，ａ２｝の探索空間におけるコスト関数面プロファイルが図１０に示すものになると、｛ａ１，ａ２｝の最適値は、収縮グリッドアルゴリズムなどの最適化方法を用いて｛−０．２５，０．０｝に収束する。

その後、｛ａ３，ａ４｝も、図１１に示すように｛０．２５，０．０｝と最適化される。その結果、｛ａ１，ａ２，ａ３，ａ４｝は｛−０．２５，０，０．２５，０｝と最適化される。

最適化された値と真値の間には差が存在する。この前方モデルが｛ａ１，ａ２，ａ３，ａ４｝の線形関数である場合、この最適化結果はほぼ最適である。

しかし、典型的には、前方モデルは非線形関数である。この場合、前方モデルは｛ａ１，ａ３｝の非線形関数である。そのため、｛ａ１，ａ３｝の組合せのコスト関数面は図１２に示すものになる。この例では、ａ３の初期値は０．５であり、これは、真値０．０からの偏りが０．５であることを示している。したがって、｛ａ１，ａ２｝がａ３＝０．５の状態で最適化されると、ａ１は−０．２５に最適化され、これが偏りになる。その結果、最適化点は、真のターゲット値と異なる極小値に向かって収束するように見える。

本発明の第２実施形態では、真のターゲット値と異なる極小値に向かって収束する最適化パラメータの問題の解法を提案する。それによって、推定の正確さが改善される。

図１３に、本発明の第２実施形態による波面測定装置１１１１の構成を示す。波面測定装置１１１１は、光源１１００、波面センサ１１１５、およびデータ解析部１１１８を備える。光源１１００は光ビームＢを発し、この光ビームはフィルタリング光学素子１１０１（例えばＮＤフィルタ）を通過し、照明光学素子１１０２によって整形される。次いで、整形されたビームは検査エレメント１１０３を照明する。検査エレメント１１０３は、レンズ、プリズム、導波路などの単一の光学エレメントとすることもできるし、光学エレメントの組合せとすることもできる。検査エレメント１１０３によって導入される波面収差を補償（低減）するために、レンズレットアレイ１１０４と検査エレメント１１０３の間に（第１実施形態で説明した）補正エレメント５０４を任意選択で配置し得る。

この光ビームは、検査エレメント１１０３を通過し、レンズレットアレイ１１０４に進み、次いで、センサアレイ１１０５を含む光検出器１１０９に至る。センサアレイ１１０５は、この光学系によって生成される光強度分布を検出する。

センサアレイ１１０５を備えた光検出器１１０９は、この光学系全体の情報を有する信号を出力し得る。光検出器１１０９からの出力は、既知の方法で、例えば、データ取得系１１２０のアナログデジタル（ＡＤ）コンバータを用いてデジタル化し得る。次いで、デジタル化されたデータはデータ解析部１１１８によって解析される。データ解析部１１１８は、検査エレメント１１０３を含むこの光学系に関する様々な光学パラメータを解析する。データ解析部１１１８は、コンピュータ制御下でデータを取得し記憶する。データ解析部１１１８では、特殊なソフトウェアを用いてデータフィッティングを行い、様々な計算を行い、データをグラフィカルな形式で表示する。

光源１１００として、パルスレーザまたは連続波（ＣＷ）レーザなどのコヒーレント光源、あるいは水銀ランプまたはタングステンランプなどの非コヒーレント光源を使用し得る。ＣＷレーザまたはランプの場合、光ビームを変調するために光チョッパなどの変調装置が必要とされることがある。光ビームを整形する照明光学素子は当業者には周知であると考えられ、本明細書では説明しない。レンズレットアレイ１１０４については、（１次元および２次元の）様々なレンズレットアレイが当技術分野で周知の商業供給業者から容易に入手可能である。アレイ内のレンズレットの形状、ピッチ、および焦点距離などの特定のパラメータは一般に、測定装置の特定の用途によって決まる。同様に、センサアレイ１１０５は、当技術分野では周知のように、単一のＣＣＤカメラとして実装し得る。あるいは、センサアレイ１１０５は、例えば大規模集積（ＶＬＳＩ）電子技術を用いて専用に設計し得る１次元または２次元のＣＭＯＳセンサアレイとして実装し得る。

図１４に、データ解析部１１１８の概略ブロック図を示す。データ解析部１１１８は、以下に示すステップＳ１１〜Ｓ１８を含む波面推定プロセスを実行し得る。そのために、データ解析部１１１８は、較正モジュール１１０７、パラメータ推定モジュール１１０８、およびハードディスクドライブまたはポータブルメモリデバイスなどの記憶ユニット１１１２（メモリ）に記憶されるデータベース１１１０を備える。データ解析部１１１８は、ステップＳ１１〜Ｓ１８のプロセスシーケンスを実行するように特に構成（プログラム）し得る一般的なプログラマブルコンピュータまたはデータ処理ユニット（マイクロプロセッサ）として実装し得る。データ解析部１１１８およびその構成エレメント（モジュール）は、個々のハードウェアエレメントまたは分散処理ユニットとして実装することもできる。本実施形態では、中央処理装置（ＣＰＵ）などのデータ処理ユニットを含むプログラマブルコンピュータとしてデータ解析部１１１８を実装し得る。一般的なプログラマブルコンピュータの代わりに、あるいはそれに加えて、専用のマイクロプロセッサまたは分散型コンピュータネットワークなどの他のハードウェアデバイスによって、ステップＳ１１〜Ｓ１８のプロセスシーケンスを実現し得る。

データ解析部１１１８は、図１３の波面測定装置の動作中に図１４に示す波面推定プロセスを実施する。ステップＳ１１で、レンズレットパラメータおよび測定された照明分布が手作業または自動的に記憶ユニット１１１２に入力される。レンズレットパラメータは、レンズレットアレイ１１０４の焦点距離およびピッチ、レンズレットアレイ内のレンズレットの数などの光学情報を含み得る。測定された照明分布は、センサアレイ１１０５からデータ取得系１１２０を介して得ることもできるし、局所記憶ユニット１１１２または遠隔記憶装置（図示せず）内のあらかじめ記憶した測定値（例えば、最後に測定した照明分布）から提供されることもある。次いで、入力されたレンズレットパラメータおよび照明分布は、較正モジュール１１０７に提供される。較正モジュール１１０７の機能には、レンズレットアレイと検出器アレイの位置関係についての情報を得る機能、あるいは、これらのアレイの位置を適切にアライメント（較正）する機能がある。つまり、以下でさらに説明するように、較正モジュール１１０７を用いて、例えば当技術分野では周知のフィードバック回路を介して実際の光学ユニットを対話形式でアライメントし得る。これに加えて、またはこの代わりに、較正モジュールを用いて、推定プロセス用の入力情報として使用し得る情報（データ）として光学ユニット間の位置関係を得ることもできる。そのために、ステップＳ１２で、レンズレットアレイおよび他の光学エレメントの位置が、検査エレメント１１０３なしで生成される参照データを用いて定義される。それによって、センサアレイ１１０５は、スクリーンショット１１２１に示すようなほぼ平面波のプロファイルを検出する。データ取得系１１２０は、この参照データをデータ解析部１１１８の較正モジュール１１０７に提供する。これにより、レンズレットアレイの位置を含めて、測定装置の初期光学構成がステップＳ１２で設定される。較正モジュール１１０７によって初期光学構成が設定されると、フロープロセスは、パラメータ推定モジュール１１０８で実施されるシーケンスステップに進む。

データ解析部１１１８によって参照データが取得され、測定装置の初期光学構成が確立されると、検査エレメント１１０３が、光源と波面センサの間の光路に位置決めされ測定が行われる。つまり、図１３で、光源１１００が発する光が検査エレメント１１０３に向けられて、この光の波面が検査光学素子を通過し、波面センサ１１１５のセンサアレイ１１０５によって検出される。ステップＳ１３で、検査エレメント１１０３によって形成されると予想される波面を記述する初期パラメータがパラメータ推定モジュール１１０８に設定される。具体的には、初期パラメータは、デフォルト値としてゼロとし得るゼルニケ係数として設定される。言い換えると、ステップＳ１３では、パラメータ推定モジュール１１０８は、ほぼ平面波（ゼルニケ係数の値がゼロ）を表し得る初期パラメータを設定する。ステップＳ１１〜Ｓ１３が完了すると、ステップＳ１４で、ステップＳ１１〜Ｓ１２の少なくとも１つから受け取った値およびステップＳ１３で設定された初期パラメータに基づいて前方モデルが計算される。つまり、ステップＳ１４では、参照データ（検査光学素子なしで得られたデータ）、測定装置の光学構成、および検査光学素子のデータに基づいて、前方伝播モデルが、ステップＳ１３で設定された初期パラメータに基づいて計算される。以下でより詳細に論ずるように、前方伝播モデルとして、小ビームに基づく伝播モデルまたは光線に基づく伝播モデルを使用し得る。最尤推定（式（１）参照）に関しては、ステップＳ１４で計算された前方伝播モデルを複数の推定モデルパラメータからなるベクトルとみなし得る。

前方モデル（シミュレーションデータ）が確立されると、プロセスはステップＳ１５に進む。ステップＳ１５で、データ取得系１１２０は、センサアレイ１１０５から実際の実験データを取得する。具体的には、ステップＳ１５で、図１３の光学装置において、検査レンズ（検査エレメント１１０３）を用いて実際の実験データが測定される。このようにして、図１４のスクリーンショット１１２２に示すように、センサアレイ１１０５が、検査エレメント１１０３によって形成される波面プロファイルを検出する。データ取得系１１２０は、この実際の実験データ（測定データ）をデータ解析部１１１８のパラメータ推定モジュール１１０８に提供する。したがって、最尤推定（式（１）参照）に関しては、光路内に検査エレメント１１０３が配置された状態で得られた測定データを、複数の観察結果からなるベクトル（複数の観察値からなるベクトル）とみなし得る。

ステップＳ１５で、データ解析部１１１８は、実際の実験データ（例えば、データ取得系１１２０のスクリーンショット１１２２に示すもの）と、ステップＳ１４で前方モデルを用いることによって得られたシミュレーションデータとを用いてパラメータの尤度を計算する。つまり、ステップＳ１５での目的は、複数の観察結果（測定データ）からなるベクトルが与えられると仮定して、観察データを生成（と合致）する最も高い尤度（最大尤度）を有する前方モデルにおけるパラメータからなるベクトルを見つけることである。パラメータの最大尤度値は、実際の実験データとシミュレーションデータの差がどのくらい小さくなり得るかを示し得る。つまり、パラメータと前方モデルが、実験データを生成した実際の光学系をどのくらい正確にモデル化し得るかを示している。ここで、尤度は、実際の実験データとシミュレーションデータが収斂（合致）し得るかどうかを決定し得る（図２２、図２６、および図３０ならびに式１６〜１９を参照して以下で論じる）推定アルゴリズムに基づくものとし得る。

そのため、ステップＳ１６で、実際の実験データ（測定データ）とシミュレーションデータ（前方モデル）が収斂したかどうかが検査される。つまり、ステップＳ１６で、前方モデルのゼルニケ係数（パラメータ）が、測定された波面プロファイルを表す同じパラメータに収斂するかどうかが検査される。ステップＳ１６で、計算結果が収斂している場合、ステップＳ１７で推定されたパラメータ（前方モデルによって確立されたパラメータ）が出力される。しかし、計算結果が収斂しなかった場合、ステップＳ１８で推定されたパラメータを変更することができ、ステップＳ１４とＳ１６の間のプロセスがこの結果が収斂するまで反復して繰り返される。具体的には、本発明では、前方モデルが測定データにほぼ収斂するまでゼルニケ係数の値を反復して変更することによって推定されたパラメータを変更し得る。例えば、式（１６〜１９）を参照して以下で論ずるように、前方モデルによって得られる複数のパラメータＭからなるベクトルにおけるデータｇの値を反復して変更し、最大尤度が得られるか検査して複数の観察値θからなるベクトルを求める。実際の実験データと、推定されたパラメータを用いることによって得られるシミュレーションデータとの差が閾値以下のとき、上記検査結果を収斂したとみなし得る。このようにして、初期ゼルニケ係数（理想波面）と、測定データ（観察された波面パラメータ）にほぼ収斂（合致）する推定されたゼルニケ係数との差を数学的に得ることによって、被検査光学エレメントによって生じた波面収差を求めることができる。

データ解析部１１１８の較正モジュール１１０７により、解析部は、光学素子の位置、検出器パラメータ、および照明分布などいくつかの光学パラメータを較正することができる。具体的には、ステップＳ１１でレンズパラメータおよび測定された照明分布が入力されると、較正モジュール１１０７は、既知の電子技術により、光学系を、その中で所望のパラメータが確立されるように較正する。例えば、較正モジュール１１０７は、電子信号を送信して光源１１００を、センサアレイ１１０５において所望の照明分布を得るのに十分な出力で点灯し得る。さらに、較正モジュール１１０７は、検査エレメント１１０３、レンズレットアレイ１１０４、および検出器１１０９の少なくとも１つの所望の位置を示す電子信号を出力し得る。較正モジュール１１０７から出力される較正信号に応答して、検査エレメント１１０３、レンズレットアレイ１１０４、または検出器１１０９は、所望の位置に移動し（またはアライメントされ）て、センサアレイ１１０５において最適な照明分布が得られる。

パラメータ推定モジュール１１０８は、検査エレメント１１０３が光学系内に存在するときの関心の対象である光学（波面）パラメータを推定する推定機能を実施する。これらのパラメータは、検査エレメント１１０３を含む各光学エレメントの位置ならびに検査エレメント１１０３の特定の平面また球面における波面偏差または強度分布を含み得る。

図１４のフロープロセスは、データベース１１１０に記憶されたデータをうまく利用できるように以下のように改変し得る。具体的には、ステップＳ１３で初期パラメータを設定した後で、データ解析部１１１８は（初期パラメータに基づいて）、データベース１１１０内を探索して、以前に確立した前方モデルに関するデータがデータベース内に存在するかを判断し得る。以前に確立した前方モデルを示す有用なデータがデータベース１１１０内に存在する場合、このデータを用いてステップＳ１５で尤度を計算する。一方、データベース１１１０内にこのようなデータが存在しない場合、先に述べたようにステップＳ１４で前方モデルを計算し、計算結果をデータベース１１１０に記憶し得る。

このデータベース１１１０は、以前に計算された複数の前方モデル出力データを記憶し得る。これらの前方モデル出力データは、前方モデルにおいて定義される計算パラメータとともに記憶され、次いで、このデータベースが索引付けされる。同じまたは類似の検査光学素子（検査エレメント）が波面測定装置によって測定されると、データ解析部におけるパラメータ推定プロセスは、索引付けされたデータベース内の類似パラメータを用いて計算されたデータを探索する。データベース１１１０内に類似のデータが見つかると、（ステップＳ１５で）このデータを用いて尤度を計算する。見つからない場合には、新たな前方モデルが被検査光学エレメントに基づいて計算され、次いで、計算されたデータがデータベースに記憶される。類似の検査光学素子が繰り返し測定される間に、データベース１１１０は類似のパラメータを用いてすでに計算されたデータを記憶し得るので、必要な前方計算数を減らすことができる。

言い換えると、繰返し測定中に処理時間が短くなり得る。したがって、波面測定装置における処理を加速することができる。こうして尤度計算が加速されることにより、より短時間でより多くの反復回数が実現し得る。

上記で述べたように、この測定系についてのＭＬＥにおいて前方モデル計算を使用し得る。前方モデルを計算するためにデータ解析部は光検出器からの出力信号を用い、この目的のために光伝播計算が実行される。

例えば、ＯｐｔｉｃａｌＲｅｓｅａｒｃｈＡｓｓｏｃｉａｔｅｓ社によって開発された市販ソフトウェアパッケージであるＣＯＤＥＶ（登録商標）におけるＢＳＰ（ビーム合成伝播、ＢｅａｍＳｙｎｔｈｅｓｉｓＰｒｏｐａｇａｔｉｏｎ）を用いて前方伝播モデルを実施し得る。ＢＳＰは、小ビームに基づく伝播方法に基づくものである。この小ビームに基づく方法は、高速フーリエ変換（ＦＦＴ）に基づく伝播モデルのいくつかの制限を克服することができる。

前方モデル構成の例として、図１５および図１６に光学構成を示す。図１５は、射出瞳平面２２００を用いた前方伝播についての光学構成を示し、図１６は、射出瞳球面２３００を用いた前方モデルについての光学構成を示す。２２０５および２３０５はレンズレットアレイを示し、２２１０および２３１０はＣＣＤイメージセンサなどの光検出器を示す。

射出瞳の平面２２００または球面２３００において波面または複雑な振幅分布を定義する。射出瞳の平面２２００または球面２３００における光線を、図示しない光源から（像面ＩＰとも称する）ＣＣＤイメージセンサ面まで追跡する。瞳における波面特性および被検査光学エレメントの他の光学設計パラメータに基づいて光線の場所および方向が計算される。

小ビームに基づく伝播では、各光線はガウシアンビームなどの細い小ビームである。イメージセンサ面上で多数の小ビームが重ね合わされ、次いで、像面上の強度プロファイルが計算される。データ解析部は、光学系またはその構成エレメントの、対象とする光学パラメータを簡単な光学系を用いて推定し得る。したがって、本実施形態では、従来の測定技術では測定し得なかった大きな収差を含む波面を測定することができる。その結果、簡単な光学構成で、レンズ、レンズユニット、または光学構成エレメントの光学性能を検証し得る。

瞳上で波面が極めて大きな偏差を有している場合、この偏差を表すにはより多くの小ビームが必要である。小ビームの数は演算リソースによって制限されるので、波面偏差量があるレベルを越えると、小ビームに基づく伝播の正確さは急激に悪くなる。

図１７に、前方モデルの光線追跡図を示す。図１７では、光源側から波面センサ側への順で、光線が検査エレメント２５０４およびレンズレットアレイ２５０５を通って進み、ＣＣＤセンサアレイ２５０６によって検出される。所定の面２５００が、検査エレメント２５０４のすぐ隣の光源側位置に確立され、レンズ設計データ（検査エレメントパラメータ）があらかじめ既知であると仮定される。このモデルにおいて検査エレメントのレンズ設計データが定義されると、所定の面２５００における波面偏差を、レンズ設計データを知らずに所定の面が確立される場合よりもはるかに簡単に推定することができる。レンズ設計データおよび所定の面が確立されているという前提に基づいて、前方伝播モデルが計算される。この前方モデルを用いて、データ解析部において波面偏差または最適化されたレンズ設計パラメータが推定される。次いで、これらの推定結果を用いて波面収差を補正し得る。

さらに、対象とする面に対して測定パラメータを外挿（投影）することによって、対象とする任意の光学パラメータを推定することもできる。図１８に、レンズユニット（レンズ群）２６０８を測定するための光学構成を示す。この場合、表面プロファイル、レンズの厚さ、レンズ間距離、偏心、傾き、または屈折率分布などの光学パラメータを、データ解析部で推定するターゲットパラメータとし得る。

例として、図１８に示すレンズユニット２６０８の波面測定を考え得る。表１に示すように、便宜上、フリンジゼルニケ多項式で波面偏差を表すが、他の関数も使用し得る。データ解析部で推定するパラメータは、これらの多項式の係数である。まず、波面の初期パラメータを、図１８のレンズユニット２６０８（検査光学素子または検査エレメント）の光源側の所定の面２６０９において定義する。ＳＨＷＦＳを形成するように配置されたレンズレットアレイ２６１０およびセンサアレイ２６１５を用いてレンズユニット２６０８を通過する光を検出する。初期レンズ設計データおよび所定の面があらかじめ既知なので、所定の面２６０９から対象とする面（表面）まで初期波面パラメータを前方伝播させることによって、レンズユニット２６０８における対象とするパラメータを推定し得る。
表１：図１８の所定の面２６０９において定義され、対象とする面または表面まで前方伝播されるゼルニケ多項式

最初に、所望の面における波面パラメータを推定するプロセスを開始するために、図１４のステップＳ１１に従ってレンズ設計パラメータおよび検出器データを示す情報が入力される。検出器データに関しては、検出器自体の特定のデータ、例えば、ＣＣＤイメージセンサの画素数（例えば１２８０×１０２４）、ＰＳＮＲ（ピーク信号とノイズの比、ＰｅａｋＳｉｇｎａｌｔｏＮｏｉｓｅＲａｔｉｏ）＝４００、およびＡＤコンバータの８ビットのデジタル化分解能を入力し得る。より重要には、検出器アレイのところで受け取られる特定の測定照明分布を入力し得る。図１９に、図１８に示すように配置された光学系から得られる検出器データの例を示す。検出器データはデータ解析部１１１８に渡される。

図１４のデータ解析部１１１８は、較正機能およびパラメータ推定機能を有する。較正機能により、光学素子の位置、検出器パラメータ、照明分布などいくつかの光学パラメータが較正される。

パラメータ推定機能により、波面収差などの光学パラメータが推定され、それによって、対象とする典型的なパラメータである表面形状または波面偏差が求められる。図１８の例では、まず、既知のレンズ設計データに基づいて所定の面２６０９における波面パラメータを推定し得るが、波面伝播を用いることによって、検査エレメント２６０８に含まれる各光学エレメントの位置または形状に対応する任意のパラメータを、対象とする特定の平面または球面における波面偏差または強度分布を推定することによって推定し得る。

上記で説明したように、確率分布関数（ＰＤＦ）をｐｒ（ｇ｜θ）で表す。ここで、ｇは観察データベクトルであり、Ｍ×１ベクトルとして特徴付けられ、θはパラメータベクトルであり、Ｐ×１ベクトルとして特徴付けられる。上記で説明したように、ｇが与えられると仮定してθの尤度が式（１）によって定義され、ｇが与えられると仮定してθの最大尤度（ＭＬ）の推定値が式（２）によって定義される。

実際のケースでは、通常、解析的にＭＬ推定値を得ることができないので、ＭＬ推定値は数値的に求めなければならない。

例えば、ｇが与えられると仮定してθの確率分布関数（ＰＤＦ）がガウシアン分布によって定義され、尤度が式（１６）によって与えられる。

・・・（１６）
ここで、ｇはパラメータデータであり、ｍはベクトルＭの指数であり、

はパラメータデータの平均値であり、前方モデルを用いて計算され、σ^２はノイズによる分散である。したがって、対数尤度を以下のように記述し得る。

・・・（１７）
その結果、ＭＬ推定値は以下のように記述し得る。

・・・（１８）
図１４に示すパラメータ推定プロセスでは、ＭＬ推定値は反復的に計算される。そのために、平均値データ

が、ある検査場所における尤度を計算するのに必要である。データ

は、このプロセスの「前方モデルを計算する」ステップで計算され、次いで、対数尤度ｌｎＬ（θ｜ｇ）がこのプロセスの「尤度を計算する」ステップで計算される。ＭＬ推定値

は最適化探索法によって求められる。そのため、この探索を最適化するために、式（１９）によって定義するように、前方モデルにおいて計算されるパラメータの各値にコスト関数が割り当てられる。

・・・（１９）

式（１８）は、

がコスト関数の最小値を有するθであることを意味する。

したがって、コスト関数プロファイル（あるいは、パラメータの数によっては表面または超表面）内で
を探索区域Ａ全体にわたって探索する。

図２０に、ゼルニケ係数のコスト関数の例を示す。図２０では、Ｃ２は、表１における２次のゼルニケ多項式の係数を示す。この場合、所定の面における理想表面には波面偏差は存在せず、そのため、最適値はゼロになる。最適化プロセスでは、解析部がコスト関数面のプロファイル上でこれらの最適係数を探す。例えば、パラメータＣ２が計算されると、他のパラメータを指定値に固定し得る。

図２１Ａ〜図２１Ｊに、様々なパラメータ（Ｃｎはｎ次ゼルニケ多項式の係数を意味する）の組合せの２次元コスト関数面を示す。本発明者らは本発明で、図２１Ｂに示す｛Ｃ４，Ｃ９｝、図２１Ｆに示す｛Ｃ２，Ｃ７｝、および図２１Ｇに示す｛Ｃ３，Ｃ８｝のコスト関数が密接に相関している（すなわち、Ｃｎが軸対称のケース）ことを見いだした。したがって、これらの組合せは同時に推定しなければならない。

言い換えると、これらの組合せのコスト関数面におけるこれらの最適パラメータの探索を実施する。これらの組合せのコスト関数を用いて、前方モデル計算数を減らすことができ、そのため、効率的なＭＬＥを実現し得る。

コスト関数面の対称性を調べて相関するパラメータを見つけるには、以下の手順を用いることができる。

（１）対称性を調べる区域Ａを決める。図２１Ａ〜図２１Ｊでは、区域Ａは、ある閾値よりも小さいコスト関数を有する係数の幅として定義される。

（２）区域Ａにおいてコスト関数面の対称性を調べる。対称性は、図２１Ａ〜図２１Ｊに示すコスト関数面を用いて調べることができる。

（３）同時に最適化すべきパラメータ（係数）の組合せを決める。図２１Ａ〜図２１Ｊでは、これらの組合せは、｛Ｃ４，Ｃ９｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ９，Ｃ１６｝と決定し得る。同時に最適化すべきパラメータの組合せは、上記手順（２）の下で確立される対称性に基づいて決定（判断）し得る。

（４）その後、データ解析部で最適化を実施する。

（５）最適化がうまく収斂しない場合、ステップ（４）の最新の出力パラメータを用いてステップ（１）を繰り返す。

図２２に、コスト関数面の対称性を調べるプロセスのフローチャートの例を示す。ステップＳ２００で、所定の面における理想パラメータが設定される。要するに、理想（平面）波面に対して、例えば、図１７に示す前方モデルの場合、所定の面においてゼルニケ係数をゼロに設定し得る。所定の面における理想パラメータは、例えば図１６に示すように、理想球面波面の関数の形で送り出すこともできる。つまり、理想値は、前方モデルに用いられる光学構成の設計値である。図２０に示し図２１Ａ〜図２１Ｊで解析される対称性を調べるための区域ＡをステップＳ２０１で求める（決定する）。次いで、ステップＳ２０２で、区域Ａ内のコスト関数面の対称性を調べる。調べた対称性に基づいて、ステップＳ２０３で、同時に最適化すべきパラメータの組合せを求める（決定する）。次いで、ステップＳ２０４で初期値として入力される実際の設計パラメータを用いて、ステップＳ２０５で最適化を実行する。理論的には、任意の初期パラメータを受け取ることができるが、実際には、初期パラメータには（例えば、所定の面におけるゼルニケ係数によって定義される）実際の設計値が用いられる。さらに、以前に測定したデータなどの他のパラメータが利用可能な場合、図２２のＳ２０４でこのような以前のデータも初期パラメータとして使用し得る。最後に、ステップＳ２０６で、コスト関数に基づいて最適化されたパラメータが出力される。パラメータが良好に最適化されない場合（Ｓ２０５のＮＯ）、ステップＳ２０１〜Ｓ２０３で、同時の最適化するパラメータの数を増やし得る。

図２２のステップＳ２０２に対応する手順（２）では、コスト関数の対称性、すなわち、最適化されるパラメータ間の相関の程度を調べるのにいくつかの方法を用い得る。例えば、対称性を調べるために用いるコスト関数の点数は９点と少なくし得る。

図２３Ａ〜２３Ｂおよび図２４Ａ〜２４Ｂに、最適化される例示パラメータのコスト関数対称性を調べる簡略化された方法を示す。図２３Ａおよび図２３Ｂは、コスト関数最小化プロセスにおいてパラメータを独立に探索し得る対称な場合を示す。Ｐ１１、Ｐ１２、Ｐ１３、Ｐ２１、Ｐ２２、Ｐ２３、Ｐ３１、Ｐ３２、およびＰ３３におけるコスト関数の９点は図２３Ａに示すように定義される。区域Ａ内の各点の位置に基づいて、点Ｐｘｘの位置に対応するコスト値Ｌｘｘをコスト関数に割り当てることができる。こうすると、図２３Ｂで、Ｌ２２＜Ｌ３２＜Ｌ１２、Ｌ２１＜Ｌ３１＜Ｌ１１、およびＬ２３＜Ｌ３３＜Ｌ１３になる。次いで、コスト関数の３つの値からなる組をそれぞれ、図２３Ｂに示すように、放物線関数に対して曲線フィッティングし得る。３つの曲線フィッティングされた放物線関数の最小値は、点線によって示されるパラメータａ１のほぼ同じ値のところで得られる。したがって、これら３つの放物線関数の軸の場所および向きがほぼ同じであり、コスト関数の値を最小値から最大値まで昇順に並べ替えることができる場合、コスト関数のパラメータは対称であると言え、独立に最適化し得る。

一方、図２４Ａおよび図２４Ｂには非対称な場合を示す。図２４Ｂでは、放物線関数の３つの最小値がａ１に沿った異なる場所で得られる。さらに、区域Ａ内の各点の位置に基づくコスト関数の値から、Ｌ２２＜Ｌ３２＜Ｌ１２、Ｌ１１＜Ｌ２１＜Ｌ３１、およびＬ３３＜Ｌ２３＜Ｌ１３が得られる。非対称の場合には、コスト関数を最小値から最大値まで昇順に並べ替えることができない。そのため、コスト関数の全体の最小値を見つけるには、パラメータａ１およびａ３を一緒に最適化する必要がある。したがって、このコスト関数は所定の区域Ａ内では非対称であると言って間違いない。これらの基準（コスト関数値のランク順序付け）は、対称（非相関）または非対称（相関）の結論を引き出すには十分ではない。しかし、最小値が単独で共通の位置で生じることは、これらのパラメータを独立に探索し得ることを保証するのに十分であり、ランク順序付けと最小値の場所を組み合わせると、対称性のありなしを適切に判断することが保証され得る。

所定の面にある量の波面偏差が存在する場合、この所定の面における波面のゼルニケ係数を表にすることができ、これを表２に示す。この場合、コスト関数面は図２５Ａ〜図２５Ｊに示すようになる。
表２：所定の面における収差を含む波面のゼルニケ係数

図２５Ａ〜図２５Ｊのコスト関数の面プロファイルは、パラメータ（ゼルニケ係数）の真値が異なるために、図２１Ａ〜図２１Ｊに示すものと異なる。具体的には、図２５Ｉの｛Ｃ４，Ｃ５｝と図２５Ｊの｛Ｃ４，Ｃ７｝のコスト関数は、図２１Ｉおよび図２１Ｊにそれぞれ示す同類のパラメータのコスト関数よりも大きな相関を示している。ゼルニケ多項式の対称性によれば、｛Ｃ４，Ｃ６｝と｛Ｃ４，Ｃ８｝も相関するはずである。したがって、ゼルニケ係数の数が少なく収斂が実現されない場合、｛Ｃ４，Ｃ５｝、｛Ｃ４，Ｃ６｝、｛Ｃ４，Ｃ７｝、および｛Ｃ４，Ｃ８｝の組合せを追加する必要があり得る。

上記で説明したように、同時に最適化しなければならない組合せは、｛Ｃ４，Ｃ９｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ９，Ｃ１６｝である。さらに、以下の組合せ（｛Ｃ４，Ｃ５｝、｛Ｃ４，Ｃ６｝、｛Ｃ４，Ｃ７｝、および｛Ｃ４，Ｃ８｝）を追加し得る。その結果、｛Ｃ２，Ｃ４，Ｃ７｝、｛Ｃ３，Ｃ４，Ｃ８｝、または｛Ｃ４，Ｃ９，Ｃ１６｝の組合せも、より効率的な収斂が得られるように働き得る。最適化の開始時において、図１９に示すような検出器データは、倍率および位置のせいで、前方モデルの結果と異なることがある。この場合、先端、傾き、および焦点の成分をあらかじめ最適化する必要がある。

先端、傾き、および焦点は、｛Ｃ２，Ｃ３，Ｃ４｝に対応する。したがって、この場合のゼルニケ係数は、図２６のフローチャートに示す手順を用いて推定し得る。図２６では、ステップＳ２３００で、各パラメータ毎の初期値が設定される。ステップＳ２３０１で、ゼルニケ多項式の係数Ｃ２、Ｃ３、およびＣ４が１回り目の収縮グリッドグローバル探索として推定される。ステップＳ２３０２、Ｓ２３０３、Ｓ２３０４、およびＳ２３０５で、引き続き１回り目の収縮グリッドグローバル探索中に、｛Ｃ５，Ｃ６｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ４，Ｃ９，Ｃ１６｝がそれぞれ推定される。

１回り目の収縮グリッドグローバル探索後、２回り目〜６回り目の収縮グリッドローカル探索を順次実行し得る。ステップＳ２３０６、Ｓ２３０７、Ｓ２３０８、Ｓ２３０９、Ｓ２３１０、およびＳ２３１１で、｛Ｃ２，Ｃ３｝、｛Ｃ４，Ｃ９｝、｛Ｃ５，Ｃ６｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ９，Ｃ１６｝がそれぞれ推定される。ステップＳ２３１２で、所望の反復回数Ｎが完了したかどうかの判断がなされる。この反復回数が完了していない場合、収縮グリッドローカル探索が繰り返される。この反復回数が完了した場合、ステップＳ２３１３で、尤度が最大である値が推定結果として出力される。

収斂の検査として、図２７に、ＭＬ推定中のゼルニケ係数の推定値と既知のターゲット値との差を示す。ターゲット値は、表２にまとめられた設計係数を意味する。図２７のデータは以下のようにして得られる。まず、前方モデルにターゲット値（表２のゼルニケ係数）を入力することによって、光分布のシミュレーションデータ（図１９など）を得ることができる。次いで、このシミュレーションデータに基づいて、ＭＬ推定下でゼルニケ係数を推定し得る。最後に、ゼルニケ係数の推定値とターゲット値の差を得ることができる。

最適化後、最適化された（推定された）係数とターゲット係数の間に差異がある場合には、以下のプロセスを用いることができる。

図２８に、｛Ｃ７，Ｃ８｝の組合せのコスト関数面を示す。ターゲット係数は｛５．０λ，５．０λ｝であるが、最適化された値は｛５．０８３４λ，３．２３９５λ｝である。このように、コスト関数面上での探索は極小値に収束したことがわかる。具体的には、本発明者らは本発明で、収差が大きい場合には、探索が誤ってグローバル最小値近傍の極小値に収束し得ることを見いだした。

図２９Ａに、上記の収縮グリッド探索を用いたグローバル最適化後の｛Ｃ３，Ｃ８｝の組合せのコスト関数面を示す。この場合、探索点は、真のターゲット値と異なる極小値に向かって収束し得る。正確さを改善するために、本発明者らは本発明で、収縮グリッド探索法に対するさらなる改善を考案した。具体的には、
（１）下記の式（２０）によって定義される単純な２次関数を、例えば最小二乗法を用いて、コスト関数の離散サンプルに対してフィッティングする。
ｆ（ｘ，ｙ）＝ａ_１ｘ^２＋ａ_２ｙ^２＋ａ_３ｘｙ＋ａ_４ｘ＋ａ_５ｙ＋ａ_６
・・・（２０）
（２）この２次関数の最小値が存在する場所を見つけ、この最小値が生じる場所を、収縮グリッド探索の次の反復のための新たな中心として用いる。この方法によって、グローバル最小値に向かって収束するように改善されることが示された。

図２９Ａに｛Ｃ３，Ｃ８｝尤度の離散データを示し、図２９Ｂに｛Ｃ３，Ｃ８｝空間のフィッティングされた面を示す。この場合、｛Ｃ３，Ｃ８｝の値＝｛９．８９９７８λ，４．６４９２６λ｝が、フィッティングされた面から最適係数として得られた。この方法を用いて、図２７に示す最適化手順が大きく改善された。

図３０に、収縮グリッド探索アルゴリズムについての改変されたフロープロセスを示す。具体的には、図３０は、フィッティング方法を用いて改善されたゼルニケ係数を推定する最適化プロセスを示している。この改善された最適化では、１回り目の収縮グリッドグローバル探索（ステップＳ２７０１〜Ｓ２７０５）と２回り目〜６回り目の収縮グリッドローカル探索（ステップＳ２７１０〜Ｓ２７１５）との間にフィッティングプロセスが実行される。ステップＳ２７００で初期値が設定された後で、ステップＳ２７０１からＳ２７０５で、｛Ｃ２，Ｃ３，Ｃ４｝、｛Ｃ５，Ｃ６｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ４，Ｃ９，Ｃ１６｝のパラメータの組合せに対して収縮グリッドグローバル探索が実行される。

ステップＳ２７０６からＳ２７０９で、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、｛Ｃ４，Ｃ９｝、および｛Ｃ５，Ｃ６｝に対するフィッティングプロセスが実行される。このフィッティングプロセスの後で、ステップＳ２７１０からＳ２７１５で、｛Ｃ２，Ｃ３｝、｛Ｃ４，Ｃ９｝、｛Ｃ５，Ｃ６｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ９，Ｃ１６｝のパラメータの組合せに対して収縮グリッド探索が実行される。ステップＳ２７１６で、所望の反復回数Ｎが完了したかどうかの判断がなされる。この反復回数が完了していない場合、収縮グリッドローカル探索が繰り返される。この反復回数が完了した場合、ステップＳ２７１７で推定値が出力される。図３０で説明した改善されたプロセスは、上記の第１実施形態の後方伝播計算に適用し得る。

図３１に、このフィッティング方法を用いたＭＬ推定中の、推定されたゼルニケ係数とターゲット係数の差を示す。この結果から、ターゲット係数と推定された係数の差が０．１λよりも小さいくなり得る。したがって、本発明者らは本発明で、フィッティング方法を用いる収縮グリッド探索が収束性能を劇的に改善し得ることを確認した。

本明細書で説明したように、最適化にかかる演算時間が短くなり、推定結果の正確さが改善されることにより、有利には、ここで説明した原理から逸脱することなく、波面測定技術が改善され得る。

第１および第２の実施形態の概念は、複数の検出器または検出器アレイから得られるデータの処理に適用可能である。複数の検出器から得られるデータは、検出器位置、照明モード、偏光状態、または検査エレメント位置を変えて取得し得る。

本発明の少なくとも１つの実施形態により、従来の測定技術では測定し得なかった大きな収差を含む波面を測定することができる。したがって、単レンズやレンズユニットなどの検査エレメントの光学性能の検査または検証を簡単な光学構成で実現し得る。

第１実施形態では、補正エレメントは、検査エレメントによって生じるある種の大きな低次波面収差を補償することができる。第２実施形態では、収縮グリッド探索のフィッティング方法により、推定値の収束が効率的かつ正確に計算されることが保証される。

例示の実施形態を参照して本発明を説明してきたが、本発明は開示した例示の実施形態に限定されるものではないことを理解されたい。添付の特許請求の範囲は、最も広い解釈に従うものとし、そのため、あらゆる改変形態ならびに等価な構造および機能が含まれるものとする。

１００波面測定装置、１０２光源、１０４ＮＤフィルタ、１０６ビームエキスパンダ、１０８投影レンズ、１１０検査光学素子、１１２レンズ群、１１４レンズレットアレイ、１１６光検出器、１１８データ解析部、１５０波面、１６０波面

Claims

ゼルニケ多項式で偏差が表される波面を最尤推定法を用いることによって推定する方法であって、
レンズレットアレイおよび検出器を含む波面センサを用いることによって得られるデータを準備する工程と、
ゼルニケ多項式の互いに異なる２つのゼルニケ係数により前記検出器上の強度分布をコスト関数として求め、前記２つのゼルニケ係数をそれぞれ座標軸とする座標面において当該コスト関数を表したコスト関数分布を生成する工程と、
前記座標軸に対する前記コスト関数分布の対称性を調べるための区域を選択する工程と、
前記区域において前記コスト関数分布の前記座標軸に対する対称性を調べる工程と、
同時に最適化すべきゼルニケ多項式の係数の組合せとして、前記対称性を調べた結果に基づいて、前記座標軸に対して非対称なゼルニケ係数の組み合わせを選択する工程と、
前記コスト関数の最小値を探索することによって前記組合せに従うゼルニケ多項式の係数について最大尤度最適化を実行する工程と、を含むことを特徴とする方法。
２次関数に対してコスト関数面をフィッティングして、前記組合せに従うゼルニケ係数の前記最適化を改善するフィッティングプロセスをさらに含むことを特徴とする請求項１に記載の方法。
前記データは、前記レンズレットアレイの前に配置される波面補正レンズを用いて、前記波面の収差を前記波面センサのダイナミックレンジに入るまで低減することによって得られることを特徴とする請求項１に記載の方法。
前記座標軸に対する対称性を調べる工程は、
前記２つのゼルニケ係数の一方を一定値として、他方のゼルニケ係数を変化させた際の前記コスト関する分布を２次関数でフィッティングし、前記２次関数からコスト関数が最小値となる前記他方のゼルニケ係数の値を求めることを、前記一方の互いに異なる複数の値の各々について行い、
前記一方の互いに異なる複数の値の各々について求められたコスト関数が最小値となる前記他方のゼルニケ係数の値の一致度から対称性を求める工程を含み、
前記座標軸に対して非対称なゼルニケ係数の組み合わせを選択する工程は、前記複数のゼルニケ係数が一致しない前記２つのゼルニケ係数の組み合わせを選択する工程を含む、ことを特徴とする請求項１に記載の方法。