(57)【特許請求の範囲】
【請求項１】
波面パラメータを推定する装置であって、
光源と、
レンズレットアレイと、
前記光源によって生成され、前記レンズレットアレイを通過した光を検出するためのセンサアレイを含む検出器と、
前記レンズレットアレイと前記光源の間に配置される波面補正エレメントと、
前記波面補正エレメントの前記光源側に位置する所定の面における波面の少なくとも１つの波面パラメータを推定するように構成されるデータ解析部とを、備え、
前記レンズレットアレイおよび前記センサアレイは波面センサを形成するように配置され、
前記波面補正エレメントは負の平凹レンズであり、前記負の平凹レンズの平坦面は前記レンズレットアレイに接触し、
前記波面補正エレメントは、前記波面補正エレメントに入射する光の波面の収差を少なくとも部分的に減少させて、前記収差が少なくとも部分的に減少した波面が前記波面センサのダイナミックレンジ内に入るように構成されることを特徴とする装置。
【請求項２】
前記レンズレットアレイは複数のレンズレットを含み、前記光源が発する前記光は前記レンズレットを通って進み、
前記センサアレイは、前記複数のレンズレットを透過した光によって形成される複数の焦点スポットを検出するように構成され、
前記データ解析部は、前記検出器から前記複数の焦点スポットの光強度分布を取得するように構成されるデータ取得系を含み、
前記データ解析部は、前記取得系によって取得された前記複数の焦点スポットの前記光強度分布を前記波面の前方伝播モデル計算値と比較することによって前記波面の収差を推定することを特徴とする請求項１に記載の装置。
【請求項３】
前記波面センサは非干渉測定センサであり、前記センサアレイは前記レンズレットアレイの焦点面に配置されることを特徴とする請求項１に記載の装置。
【請求項４】
前記波面センサはシャック−ハルトマン波面センサであることを特徴とする請求項１に記載の装置。
【請求項５】
前記波面補正エレメントは、前記波面補正エレメントに入射する光の波面の収差項を部分的にまたは全体的に補正することにより前記波面補正エレメントに入射する光の波面の収差を少なくとも部分的に減少させることを特徴とする請求項１乃至４のいずれか１項に記載の装置。
【請求項６】
波面パラメータを推定する方法であって、
波面センサと光源の間に検査エレメントを配置することと、
前記光源によって生成される光で、前記光の波面が前記検査エレメントを通過し前記波面センサによって検出されるように、前記検査エレメントを照明することと、
前記波面が前記波面センサのダイナミックレンジ内に留まるように、前記波面センサと前記検査エレメントの間に配置される波面補正エレメントを用いて前記波面の収差を部分的に減少させることと、
データ解析部を用いて、前記波面補正エレメントの前記光源側に位置する所定の面における前記波面の少なくとも１つのパラメータを推定すること、とを含み、
前記波面センサは、レンズレットアレイと、前記光源によって生成され、前記レンズレットアレイを通過した光を検出するためのセンサアレイを含む検出器と、を含み、
前記波面補正エレメントは、負の平凹レンズであり、前記レンズレットアレイは、前記負の平凹レンズの平坦面上に形成されていることを特徴とする方法。
【請求項７】
前記波面補正エレメントは、前記波面補正エレメントに入射する光の波面の収差項を部分的にまたは全体的に補正することにより前記波面補正エレメントに入射する光の波面の収差を少なくとも部分的に減少させることを特徴とする請求項６に記載の方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本出願は、参照により開示全体が本明細書に組み込まれている、２０１１年１０月１４日に出願した仮出願第６１／５４７５８３号の優先権を主張するものである。
【０００２】
本発明は、波面形状や収差などの光学波面パラメータを推定する装置および方法に関する。
【背景技術】
【０００３】
波面測定は通常、光学的な表面の質を検査し、光学エレメントの光学性能を評価するために用いられる。波面測定は、高精度光学系のアライメントにも用いられる。波面は、すべての光線の位相が同じである点の集まり（線、または３次元の波の伝播では面）である。最も簡単な形態の波面は平面波の波面であり、この場合、光線は互いに平行であり、平らな波面がセンサに当たる。レンズやミラーなどの光学エレメントを通る波面の伝播では、一般に、レンズの厚さ、表面形状の不完全さ、屈折率のばらつき、および他の要因により波面の形状が変化する。波面形状の望まれない変化は収差として知られている。そのため、波面プロファイルを知り波面収差を補正することは、光学エレメントを設計する際や、新たに設計した光学系の性能を評価する際に極めて重要である。例えば、結像系一式を組み立てる前に、その系に含まれる各光学単位（単位レンズ）の性能を確認することが必要である。単位レンズまたは単レンズ自体はそれぞれある種の収差を有し得るので、結像レンズの質を高精度に制御することが必要である。大きな波面収差を測定する一応用例は、レンズやミラーなどの非球面光学エレメントの曲率中心の検査である。
【０００４】
光ビームの波面の質を測定する従来の方法では干渉測定波面センサを用いる。干渉測定波面センサでは、光源ビームの小部分を空間フィルタリングすることによって球面参照波を生成し、続いてそれを元の波面と合成してインターフェログラムを生成する。当技術分野では周知のように、インターフェログラムの干渉縞を解析して光ビームの質を評価することができる。しかし、元のビームを分割しそれを参照ビームと再合成すると、例えば、光路に誤差があったり、光学構成エレメントのアライメントが不適切な場合には、システム収差が発生しやすい。従来の波面測定法に代わるものとして、元のビームを分割し再合成する必要のないシャック−ハルトマン（Ｓｈａｃｋ−Ｈａｒｔｍａｎ）波面センサなどの非干渉測定波面センサ（ＮＩＷＦＳ）を用いる方法がある。
【０００５】
シャック−ハルトマン波面センサ（ＳＨＷＦＳ）などのＮＩＷＦＳは、干渉測定センサよりもダイナミックレンジが広い。ＳＨＷＦＳセンサの基本的かつ通常用いられている一構成は、レンズレットアレイと、レンズレットアレイの後側焦点面に配置される光検出器（典型的にはＣＣＤカメラ）とからなる。シャック−ハルトマン波面センサは、２次元（２Ｄ）レンズレットアレイを用いることによって測定対象の入射ビームの波面を複数の小ビームに分割する。このレンズレットアレイのレンズレットはそれぞれ、光検出器の表面に別々の独立した焦点（スポット）を生成する。各スポットの重心位置は、参照ビームと歪んだビームの間の波面収差を示す距離および方向によって変位する。したがって、ＳＨＷＦＳによる波面測定は、参照（平面）波面に対する収差を含む波面の局所的な勾配を推定することに基づいている。一般に、波面推定の手順は、位相が波面の局所的な勾配の集合として示されるか、開口全体にわたって求められるある種のモーダル関数の係数の集合として示されるかに応じて、ゾーナル式およびモーダル式のいずれかに分類し得る。後者では、焦点スポットの変位をゼルニケ多項式の形で表し得る。
【０００６】
ＳＨＷＦＳは干渉計測波面センサに対して以下のいくつかの利点を有する。ＳＨＷＦＳは干渉計測センサよりもダイナミックレンジが広い。入射する放射はコヒーレントでなくてよい。ＳＨＷＦＳは単一像から波面情報のすべてを取得し得るので、露光時間を短くすることができ、そのため、振動の影響を受けにくい。より重要には、ＳＨＷＦＳでは放射照度分布および位相分布がともに得られる。
【０００７】
図１に、ＳＨＷＦＳを利用する波面測定装置１００の構成例を示す。図１に示すように、波面測定装置１００は、典型的には、光源１０２、ニュートラルデンシティ（ＮＤ）フィルタ１０４、ビームエキスパンダ１０６、投影レンズ１０８、検査光学素子（検査エレメントまたは試料）１１０、レンズ群１１２、レンズレットアレイ（マイクロレンズアレイ）１１４、光検出器１１６（ＣＣＤセンサ）、およびデータ解析部１１８を含み得る。このように配置された波面測定装置１００は、試料（検査光学素子１１０）がそれを通過する光の波面に及ぼす影響を特徴付けるように構成されている。
【０００８】
図２に、ＳＨＷＦＳの光学構成をより詳細に示す。図２では、レンズレットアレイ１１４を通過した光が入射するセンサアレイ２１０が光検出器１１６に含まれている。センサアレイ２１０上の焦点スポット２５０、２５１、２５２、２５３、および２５４の場所は、入射する波面１５０の局所的な傾きによって決まる。波面１５０の局所的な傾きは、図１の検査光学素子１１０による収差によって生じる。この局所的な傾きは、焦点スポットの場所の変化によって計算し得る。レンズレットアレイ１１４の各レンズレットから得られる局所的な傾きの情報を用いることによって波面を再構成することができる。
【０００９】
波面の偏差量がＳＨＷＦＳのダイナミックレンジ内であれば、センサアレイ２１０上の各焦点スポットの位置を別々に検出し、対応するレンズレットに割り当てることができ、波面プロファイルを容易に求めることができる。しかし、図３に示すように波面の偏差がＳＨＷＦＳのダイナミックレンジを越える場合、センサアレイＳＨＷＦＳはもはや波面を解析することができない。具体的には、図３では、レンズレットアレイ１１４は、焦点スポット３５０、３５１、３５２、３５３、および３５４を形成する。しかし、波面１６０が大きな収差を含むために、焦点スポット３５０はセンサアレイ２１０の表面の範囲からはみ出てしまい、焦点スポット３５３と３５４は互いに交差してしまう（スポット交差）。つまり、焦点スポット３５３および３５４を形成するビームは互いに重なり合い、交差した後でセンサアレイ２１０上のそれぞれの場所に結像する。図３に示す状況では、焦点スポット３５０はセンサアレイ２１０の範囲からはみ出しているので、焦点スポット３５０を検出することができず、焦点スポット３５３および３５４はスポット交差のために明確に割り当てることができない。
【００１０】
図４に、収差を含む波面を検出したセンサアレイ２１０からの出力データの例を示す。図４では、波面の外側領域の焦点スポットを示す関心領域Ｓが拡大されて、焦点スポットが互いに交差したり、センサアレイ２１０の有効区域からはみ出したりすることの影響がよりよく示されている。
【００１１】
ＳＨＷＦＳのダイナミックレンジを拡大する技術や波面収差を解析的に補償する技術が以前に提案されており、その開発が継続中である。こうした比較的注目されている技術をいくつか以下にまとめる。
（１）ヌルレンズ
ヌルレンズは、予想される波面収差を補償または除去するように特に設計された１組の光学素子を含む。ヌルレンズ技術は、（製造誤差がないという前提で）検査光学素子の波面収差を完全に補償し得るので、レンズレットアレイ上で波面偏差を効果的に打ち消すことができる。しかし、この技術では、正確な測定を行うためには極めて正確なヌルレンズを製作することが必要である。したがって、ヌルレンズの製作コストが途方もなく高額になり得る。さらに、このようなヌルレンズは波面収差を予想し得る特定の検査光学素子に対して設計されるので、形状または特性が異なる検査光学素子によって形成される他の波面にはこの技術を適用し得ないことがある。ヌルレンズ技術の例が、Ｚｍｅｋに付与された米国特許第５，２３３，１７４号に記載されている。この特許を参照により本明細書に組み込む。
（２）推定技術
ヌルレンズ技術の代わりに、収差を含む波面を測定する波面推定技術が提案されており、典型的には補正アルゴリズムを用いて収差を補正または補償する。推定技術の一例が、参照文献１：Ｍｉｃｈａｅｌ  Ｃ．Ｒｏｇｇｅｍａｎｎ、Ｔｉｍｏｔｈｙ  Ｊ．Ｓｃｈｕｌｚ、Ｃｈｅｅ  Ｗ．Ｎｇａｉ、およびＪａｓｏｎ  Ｔ．Ｋｒａｆｔ、「Ｊｏｉｎｔ  ｐｒｏｃｅｓｓｉｎｇ  ｏｆ  Ｈａｒｔｍａｎｎ  ｓｅｎｓｏｒ  ａｎｄ  ｃｏｎｖｅｎｔｉｏｎａｌ  ｉｍａｇｅ  ｍｅａｓｕｒｅｍｅｎｔｓ  ｔｏ  ｅｓｔｉｍａｔｅ  ｌａｒｇｅ  ａｂｅｒｒａｔｉｏｎｓ：ｔｈｅｏｒｙ  ａｎｄ  ｅｘｐｅｒｉｍｅｎｔａｌ  ｒｅｓｕｌｔｓ」、Ａｐｐｌ．Ｏｐｔ．３８、２２４９〜２２５５（１９９９）に開示されている。
【００１２】
別の波面推定技術では、最尤推定法（Ｍａｘｉｍｕｍ  Ｌｉｋｅｌｉｈｏｏｄ  Ｅｓｔｉｍａｔｉｏｎ、ＭＬＥ）を用いて波面を再構築する。ＭＬＥ技術の例が、参照文献２：Ｈａｒｒｉｓｏｎ  Ｈ．Ｂａｒｒｅｔｔ、Ｃｈｒｉｓｔｏｐｈｅｒ  Ｄａｉｎｔｙ、およびＤａｖｉｄ  Ｌａｒａ、「Ｍａｘｉｍｕｍ−ｌｉｋｅｌｉｈｏｏｄ  ｍｅｔｈｏｄｓ  ｉｎ  ｗａｖｅｆｒｏｎｔ  ｓｅｎｓｉｎｇ：ｓｔｏｃｈａｓｔｉｃ  ｍｏｄｅｌｓ  ａｎｄ  ｌｉｋｅｌｉｈｏｏｄ  ｆｕｎｃｔｉｏｎｓ」、Ｊ．Ｏｐｔ．Ｓｏｃ．Ａｍ．Ａ．、２４、３９１〜４１４（２００７）に開示されている。
（３）ステッチング
米国特許出願公開第２００９／０２８４７５３号には、数学的方法を用いて一連の波面測定値を合わせて「ステッチ」する技術が記載されている。各測定毎に、異なる焦点、波面の傾き、または参照収差をダイナミックレンジ限界開口と合わせて用いる。この技術は、センサのダイナミックレンジを効果的に拡大し得ると言われているが、一連の波面測定を行い、数学的方法を用いて得られた測定値を合わせてステッチすることは計算集約的プロセスであることが、当業者には容易に理解されよう。
【先行技術文献】
【特許文献】
【００１３】
【特許文献１】米国特許第５，２３３，１７４号
【特許文献２】米国特許出願公開第２００９／０２８４７５３号
【非特許文献】
【００１４】
【非特許文献１】Ｍｉｃｈａｅｌ  Ｃ．Ｒｏｇｇｅｍａｎｎ、Ｔｉｍｏｔｈｙ  Ｊ．Ｓｃｈｕｌｚ、Ｃｈｅｅ  Ｗ．Ｎｇａｉ、およびＪａｓｏｎ  Ｔ．Ｋｒａｆｔ、「Ｊｏｉｎｔ  ｐｒｏｃｅｓｓｉｎｇ  ｏｆ  Ｈａｒｔｍａｎｎ  ｓｅｎｓｏｒ  ａｎｄ  ｃｏｎｖｅｎｔｉｏｎａｌ  ｉｍａｇｅ  ｍｅａｓｕｒｅｍｅｎｔｓ  ｔｏ  ｅｓｔｉｍａｔｅ  ｌａｒｇｅ  ａｂｅｒｒａｔｉｏｎｓ：ｔｈｅｏｒｙ  ａｎｄ  ｅｘｐｅｒｉｍｅｎｔａｌ  ｒｅｓｕｌｔｓ」、Ａｐｐｌ．Ｏｐｔ．３８、２２４９〜２２５５（１９９９）
【非特許文献２】Ｈａｒｒｉｓｏｎ  Ｈ．Ｂａｒｒｅｔｔ、Ｃｈｒｉｓｔｏｐｈｅｒ  Ｄａｉｎｔｙ、およびＤａｖｉｄ  Ｌａｒａ、「Ｍａｘｉｍｕｍ−ｌｉｋｅｌｉｈｏｏｄ  ｍｅｔｈｏｄｓ  ｉｎ  ｗａｖｅｆｒｏｎｔ  ｓｅｎｓｉｎｇ：ｓｔｏｃｈａｓｔｉｃ  ｍｏｄｅｌｓ  ａｎｄ  ｌｉｋｅｌｉｈｏｏｄ  ｆｕｎｃｔｉｏｎｓ」、Ｊ．Ｏｐｔ．Ｓｏｃ．Ａｍ．Ａ．、２４、３９１〜４１４（２００７）
【非特許文献３】Ｋｒｉｅｇら、「Ｉｎｖｅｒｓｅ  ｐｒｏｐａｇａｔｉｏｎ  ａｌｇｏｒｉｔｈｍ  ｆｏｒ  ａｎｇｓｔｒｏｍ  ａｃｃｕｒａｃｙ  ｉｎｔｅｒｆｅｒｏｍｅｔｅｒ」、Ｉｎｔｅｒｆｅｒｏｍｅｔｒｙ  ＸＩＩ：Ｔｅｃｈｎｉｑｕｅｓ  ａｎｄ  Ａｎａｌｙｓｉｓ、５５３１巻、Ｐｒｏｃｅｅｄｉｎｇｓ  ｏｆ  ｔｈｅ  Ｓｏｃｉｅｔｙ  ｏｆ  Ｐｈｏｔｏ−Ｏｐｔｉｃａｌ  Ｉｎｓｔｒｕｍｅｎｔａｔｉｏｎ  Ｅｎｇｉｎｅｅｒｓ（ＳＰＩＥ）、Ｃｒｅａｔｈら監修、２００４
【非特許文献４】Ｊａｅ  Ｍｙｕｎｇ、「Ｔｕｔｏｒｉａｌ  ｏｎ  Ｍａｘｉｍｕｍ  Ｌｉｋｅｌｉｈｏｏｄ  Ｅｓｔｉｍａｔｉｏｎ」、Ｊｏｕｒｎａｌ  ｏｆ  Ｍａｔｈｅｍａｔｉｃａｌ  Ｐｓｙｃｈｏｌｏｇｙ、４７巻、９０〜１００（２００３）
【非特許文献５】Ｌ．Ｒ．Ｆｕｒｅｎｌｉｄ、Ｊ．Ｙ．Ｈｅｓｔｅｍａｎ、およびＨ．Ｈ．Ｂａｒｒｅｔｔ、「Ｆａｓｔ  ｍａｘｉｍｕｍ−ｌｉｋｅｌｉｈｏｏｄ  ｅｓｔｉｍａｔｉｏｎ  ｍｅｔｈｏｄｓ  ｆｏｒ  ｓｃｉｎｔｉｌｌａｔｉｏｎ  ｃａｍｅｒａｓ  ａｎｄ  ｏｔｈｅｒ  ｏｐｔｉｃａｌ  ｓｅｎｓｏｒｓ」、Ｐｒｏｃ．ＳＰＩＥ、６７０７、２００７
【発明の概要】
【発明が解決しようとする課題】
【００１５】
当技術分野の上記状況に鑑みて、本発明者らは本発明で、波面パラメータを推定する新規技術を開発した。この技術の利点は、中でもとりわけ、波面誤差を補正して打ち消すと同時に、波面パラメータの測定および解析に必要とされる時間を短縮し得ることである。
【課題を解決するための手段】
【００１６】
本発明の一態様によれば、波面パラメータを推定する装置は、光源と、レンズレットアレイと、光源によって生成され、レンズレットアレイを通過した光を検出するための検出器と、レンズレットアレイと光源の間に配置される波面補正エレメントと、補正エレメントの光源側に位置する所定の面における波面の少なくとも１つの波面パラメータを推定するように構成されるデータ解析部とを備え、レンズレットアレイおよびセンサアレイは非干渉測定波面センサを形成するように配置され、前記波面補正エレメントは負の平凹レンズであり、前記負の平凹レンズの平坦面は前記レンズレットアレイに接触し、波面補正エレメントは波面の収差を少なくとも部分的に減少させて、前記収差が少なくとも部分的に減少した波面が前記波面センサのダイナミックレンジ内に入るように構成される。
【００１８】
本発明の別の態様によれば、波面パラメータを推定する方法は、波面センサと光源の間に検査エレメントを配置することと、前記光源によって生成される光で、前記光の波面が前記検査エレメントを通過し前記波面センサによって検出されるように、前記検査エレメントを照明することと、前記波面が前記波面センサのダイナミックレンジ内に留まるように、前記波面センサと前記検査エレメントの間に配置される波面補正エレメントを用いて前記波面の収差を部分的に減少させることと、データ解析部を用いて、前記波面補正エレメントの前記光源側に位置する所定の面における前記波面の少なくとも１つのパラメータを推定すること、とを含み、前記波面センサは、レンズレットアレイと、前記光源によって生成され、前記レンズレットアレイを通過した光を検出するためのセンサアレイを含む検出器と、を含み、前記波面補正エレメントは、負の平凹レンズであり、前記レンズレットアレイは、前記負の平凹レンズの平坦面上に形成されている。
【００１９】
本発明のさらなる特徴は、例示実施形態の以下の説明を添付の図面と併せ読めば、当業者には明らかである。
【図面の簡単な説明】
【００２０】
【図１】シャック−ハルトマン波面センサ（ＳＨＷＦＳ）を利用する波面測定装置の光学構成を示す図。
【図２】収差を含まない波面が入射した状態のＳＨＷＦＳの光学構成を示す図。
【図３】大きな収差を含む波面が入射したときのＳＨＷＦＳの動作を示す図。
【図４】収差を含む波面が入射したときにＳＨＷＦＳのセンサアレイから得られる出力データを示す図。
【図５】本発明の一態様による、波面パラメータを推定する波面測定装置の例を示す図。
【図６】波面補正エレメントが間に介在する状態での検査エレメントから光検出器までの光線追跡図を示す図。
【図７Ａ】検査エレメントがある場合とない場合の光検出器によって検出される像の例を示す図。
【図７Ｂ】検査エレメントがある場合とない場合の光検出器によって検出される像の例を示す図。
【図８】平行平面板が間に介在する状態での検査エレメントから光検出器までの光線追跡図を示す図。
【図９】収縮グリッド探索アルゴリズムを用いる最適化プロセスをグラフィカルに示す図。
【図１０】パラメータ｛ａ１，ａ２｝の探索空間におけるコスト関数面プロファイルを示す図。
【図１１】パラメータ｛ａ３，ａ４｝の探索空間におけるコスト関数面プロファイルを示す図。
【図１２】パラメータ｛ａ１，ａ３｝の組合せのコスト関数面を示す図。
【図１３】本発明の第２態様による波面パラメータを推定する波面測定装置の例を示す図。
【図１４】データ解析部の概略図を示す図。
【図１５】光学構成を示す図。
【図１６】光学構成を示す図。
【図１７】光学構成を示す図。
【図１８】光学構成を示す図。
【図１９】シミュレーションされた検出器データの例を示す図。
【図２０】ゼルニケ係数に対するコスト関数の変化を示す図。
【図２１Ａ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２１Ｂ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２１Ｃ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２１Ｄ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２１Ｅ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２１Ｆ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２１Ｇ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２１Ｈ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２１Ｉ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２１Ｊ】様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２２】最適化手順のフローチャートを示す。
【図２３Ａ】対象なコスト関数面プロファイルの９点を示す図。
【図２３Ｂ】対象なコスト関数面プロファイルの９点を示す図。
【図２４Ａ】非対象なコスト関数面プロファイルの９点を示す図。
【図２４Ｂ】非対象なコスト関数面プロファイルの９点を示す図。
【図２５Ａ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２５Ｂ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２５Ｃ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２５Ｄ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２５Ｅ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２５Ｆ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２５Ｇ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２５Ｈ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２５Ｉ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２５Ｊ】例示の波面の様々なパラメータの組合せの２次元コスト関数面を示す図。
【図２６】収縮グリッド探索アルゴリズムに基づいてゼルニケ係数を推定するフロープロセスを示す図。
【図２７】最尤（ＭＬ）推定プロセスに用いられるゼルニケ係数の推定値とターゲット値の差を示す図。
【図２８】係数｛Ｃ７，Ｃ８｝の組合せのコスト関数面を示す図。
【図２９Ａ】係数｛Ｃ３，Ｃ８｝の組合せのコスト関数面を示す図。
【図２９Ｂ】パラメータ｛Ｃ３，Ｃ８｝のフィッティング面を示す図。
【図３０】フィッティング方法を用いる収縮グリッド探索アルゴリズムに基づいてゼルニケ係数を推定する改善された最適化プロセスを示す図。
【図３１】フィッティング方法を用いるＭＬ推定中の、推定されたゼルニケ係数とターゲット係数の差を示す図。
【発明を実施するための形態】
【００２１】
本発明による例示実施形態を以下に説明する。本明細書では、波面は、この波面を構成するすべての光線がほぼ同じ位相を有する点の集まり（線、または３次元の波伝播では面）とみなし得る。「波面」は一般に光学的波面の形で記述される。当業者には理解されるように、「光学の」、「光の」などの用語は、明確にことわらない限り可視光に限定されるものではない。具体的には、ほぼ同じ波面特性を有する波を生成し得る電磁スペクトルは、その波が赤外、可視、紫外、またはより大きいエネルギー範囲（例えば、Ｘ線およびガンマ線）の電磁スペクトルのいずれであろうが、かなりの範囲のエネルギーにわたっていることを当業者は理解しているはずである。さらに、波面収差について言及する場合、「収差」という用語は、予想、推定、または確定された参照値からの任意の偏差であり、これにはいくつかの（最小許容差を超える）偏差の和も含まれることを理解されたい。したがって、波面収差の「補正」とは、予想、推定、または確定された参照値からの任意の偏差を減少させる、最小限に抑える、または補償することを指し得る。さらに、波面「パラメータ」について言及する場合、波面の傾き、形状、位相、表面、場所、速さなどを特徴付ける任意の値または値の群をパラメータとみなし得る。
【００２２】
本発明の一態様は、波面形状や波面収差などの光学系の波面パラメータを推定する装置および方法に関するものである。図５に、本発明の一実施形態による、波面パラメータを推定する装置５００の例示セットアップを示す。
【００２３】
図５では、装置５００は、光源５０１、波面センサ５１０、およびデータ解析部５１８を備える。光源５０１は、コリメート光源または点光源とし得る単色光源として実装し得る。光学検査エレメント（すなわち、被検査光学エレメント）５０３は、光源５０１と波面センサ５１０の間に光路に沿って配置される。必要に応じて、光源５０１と検査エレメント５０３の間にレンズ５０９およびビームエキスパンダ５０２を配置し得る。ＮＤフィルタ（図示せず）もこのセットアップに適用し得る。
【００２４】
波面センサ５１０は、好ましくは非干渉測定波面センサ（ＮＩＷＦＳ）として実装される。レンズレットアレイ５０５および光検出器５０６がこの非干渉測定波面センサを形成するように配置される。本実施形態では、ＮＩＷＦＳの例としてシャック−ハルトマン波面センサ（ＳＨＷＦＳ）を用いる。ただし、本実施形態はＳＨＷＦＳに限定されるものではない。多くのタイプの非干渉測定波面センサがあり、それぞれがそれ自体の特定の配置を有する傾向がある。添付の特許請求の範囲に含まれる限り、任意の機器を使用して波面パラメータを測定し得る。本発明の一態様では、負の平凹レンズなどの補正エレメント５０４（波面補正エレメント）がＳＨＷＦＳの光源側に導入される。具体的には、補正エレメント５０４の例として負の平凹レンズが検査エレメント５０３とレンズレットアレイ５０５の間に配置される。光源５０１によって生成される光は、レンズ５０９およびビームエキスパンダ５０２を通って進み、検査エレメント５０３に入射するビームを形成する。このビームは、検査エレメント５０３を通過した後で補正エレメント５０４に進み、レンズレットアレイ５０５を通過して、ＣＣセンサアレイなどの光検出器５０６に至る。光検出器５０６に入射する光の強度が検出され、データ解析部５１８が光検出器５０６からの信号を電子的に処理（解析）し得る。データ解析部５１８は、コンピュータ制御下でデータを取得し記憶する。データ解析部５１８は、特殊なソフトウェアを用いてデータフィッティングを行い、様々な計算を行い、データをグラフィカルな形式で表示する。
【００２５】
補正エレメント５０４（負の平凹レンズなど）は、未補正の波面をＳＨＷＦＳのダイナミックレンジよりも大きくしてしまういくつかの収差項を部分的または全体的に補正することによってデータ取得ステップおよび実際の収差の推定ステップをより正確かつより実用的にする単一または一群の高精度光学構成エレメントを含む。補正エレメントは光源５０１とレンズレットアレイ５０５の間であればどこにでも配置し得るが、レンズレットアレイと検査エレメントの間に補正エレメントを配置することが好ましい。つまり、補正エレメント５０４は、レンズレットアレイ５０５の光源側に配置される。単一の平凹レンズの代わりに、検査エレメントの特徴（特性）に基づいて、単一の凹レンズ、単一の凸レンズ、単一の平凸レンズ、またはこれらの組合せを補正エレメント５０４として使用してもよい。実際、補正エレメント５０４は、補正エレメント５０４の表面がレンズレットアレイ５０５と接触し得るようにレンズレットアレイ５０５と一体化することさえ可能である。
【００２６】
図６に、検査エレメント５０３の出射面から像面ＩＰ（すなわち、光検出器５０６のセンサアレイが配置される面）までの光線追跡図の例を示す。検査エレメント５０３の最終面から発する光を、像面に向かって凹である軌跡を有する波面として、図の最も左側に示す。ある実験配置では、この波面の近軸焦点は、検査エレメント５０３の最終面から右に１１０．６４ｍｍ離れていた。波面収差を補正するか、または少なくとも部分的に補償するために、補正エレメント５０４（例えば、ｆが−７５ｍｍで、φが５０．８ｍｍの平凹レンズ）を検査エレメント５０３の後ろに配置した。２５×２５個の焦点距離４０ｍｍのレンズレットをピッチ１ｍｍで配列したアレイ（図示せず）を補正エレメント５０４の後ろ（図６では右）に配置した。実装を容易にするためにレンズレットアレイを負レンズ（補正エレメント）の平坦面上に製作することもできるし、補正レンズに隣接してその検出器側に配置することもできる。３６×３６ｍｍ２の正方形センサアレイを含む光検出器（図５の５０６）をレンズレットアレイの焦点面に配置する。レンズレットアレイの焦点面で、光検出器のセンサアレイの各エレメントは、レンズレットアレイ内の対応するレンズレットによって形成される光ビームを受け取る。原理的には、センサアレイ全体からのデータを用いてレンズレットアレイに入射する光ビームの波面を電子的に再構築し得る。
【００２７】
検出器アレイによって検出される波面は多くの異なる方法で表現し得る。波面は、例えば、任意の平面または曲面における光の場の位相に関する１組のゼルニケ係数の形で表現し得る。グリッド上のサンプルとして波面の位相の形で波面を表現することもできる。あるいは、波面は、電場の位相に関する解析的表現によって特徴付けることができる。実際、各単一レンズレットに関連するセンサアレイから得られるデータの副集合をいくつか用いて、波面の局所的な傾きを計算し得る。次いで、これらの局所的な傾きを用いて光学系の他の面での波面パラメータを推定する。天文計測に応用する適応光学素子の場合、波面は大気中を進むが、局所的な傾きの観察値を求め、大気状態を知ることにより、波面パラメータを簡単かつ正確に推定し得る。しかし、天文測定に応用する適応光学素子の場合と異なり、製造または実験的なセットアップにおける波面は大きな収差を含んでおり、大気中で歪んだ典型的な波から大きく逸脱している。そのため、高精度光学素子では、局所的な傾きの値を観察し、大気状態を知るだけでは、波面パラメータを簡単かつ正確に推定し得ない。
【００２８】
具体的には、図５および図６の配置では、光は、検査エレメント５０３から空気（または別の既知の均質な媒質）中を補正エレメント５０４まで伝播し、次いでレンズレットアレイ５０５に至ると仮定する。図６に示す構成では、検査エレメント５０３からの波面が完全な球面波であれば、補正エレメント５０４からコリメートされた光として出射されるように補正エレメント５０４により光が曲げられる。補正エレメント５０４の後ろにはレンズレットアレイ５０５が配置される（ただし、直接接触する必要はない）。レンズレットアレイ５０５の各レンズレットは局所的な波面をサンプリングし、光検出器５０６の対応する領域に像を生成する。
【００２９】
検査エレメント５０３から出射する波がほぼ球面波の場合、光検出器５０６上の像は、図７Ａに示すもののように見えるはずである。一方、検査エレメント５０３から出射する波に少なくともいくらか収差が含まれている場合、光検出器５０６上の像は図７Ｂに示すもののように見えることがある。図６ならびに図７Ａおよび図７Ｂを見れば、レンズレットアレイ５０５のレンズレットはそれぞれ光検出器５０６上にスポットを生成し得ることが理解されよう。このスポットパターンおよびスポット自体の特徴、例えば勾配や重心位置は、検査エレメント５０３から出射する波面の収差によって、かつ、補正エレメント５０４およびレンズレットアレイ５０５の光学的かつ物理的な特性によって決まる。しかし、本発明のこの態様では、検査エレメント５０３から出射する波面の収差を大きく減少させるように補正エレメント５０４が特に設計（構成）されており、そのため、波面推定をより容易に実施し得る。
【００３０】
より具体的には、補正エレメント５０４の有用性は、図８に示す光線追跡を調べ、それを図６と比較することによってよりよく理解し得る。図８には、（図６に示す）補正エレメント５０４が取り除かれ、その代わりに光学的に透明な平行平面板５２４が置かれた場合を示す。レンズレットアレイ５０５および光検出器５０６は図６に示すものから変えていない。図８から、補正エレメント５０４がないと、個々のレンズレットからの光線が光検出器５０６の中心領域に集束されることが理解されよう。さらに、図８の場合、図６の場合に反して、個々のレンズレットからの焦点スポットが交差して隣接するレンズレットの検出器領域に入りやすい（スポット交差）。
【００３１】
言い換えると、波面に大きな収差が含まれる状態では、シャック−ハルトマン波面センサにおける焦点スポットは、ほぼ規則的なグリッドに整列される単純なスポットにならない。そうではなく、スポットが複雑に分散し、おそらくは隣接するレンズレットによって生成されるスポットと重なり合い干渉し合う。原理上は、これらの複雑なパターンは測定中の波面についての情報を多く含んでいるが、実際には、これらのパターンを正確に妥当な時間内で解釈するのは極めて難しいことがある。
【００３２】
レンズレットアレイの光源側に負レンズなどの適切に設計された補正エレメントを追加すると、収差を含む波における非合焦や低次の球面収差からの大きな寄与を相殺または減少させるように作用し、残余の収差は、検出器の表面全体にわたって、ただしその範囲内できれいな輪郭の焦点スポットを得るのに十分に小さいことが図６と図８の比較から理解されよう。言い換えると、負の凹レンズなどの形態の補正エレメント５０４は、検査エレメント５０３によって生じる波面収差を少なくとも部分的に補正し、波面センサのダイナミックレンジを広げることができ、それによって、波面推定装置内で検査物体のセンサ側のほぼすべての面における波面パラメータの迅速な推定が容易になる。
【００３３】
図６に示すセットアップにおける光線追跡は、検査エレメント５０３から出射する波面のパラメータを推定する反復手順における前方モデルとして効果的に使用し得る。この推定プロセスは、波面パラメータの最尤（ＭＬ）推定によって実施し得る。最尤推定問題は一般に、パラメータ空間全体にわたる探索として定式化され、式（１）によって与えられる。
【００３４】
【数１】
0005564545.tif 000002
・・・（１）
ここで、θは複数のパラメータからなるベクトルであり、ｇは複数の観察値（データ）からなるベクトル、
【００３５】
【数２】
0005564545.tif 000003
【００３６】
は複数の推定モデルパラメータからなるベクトルである。解くべき一般問題は、１組の観察値ｇが与えられると仮定して、これらの観察データを生成する最も高い確率（最大尤度）を有するパラメータθの組はどれかという問いの形で提示される。
【００３７】
本例では、ＭＬ推定アルゴリズムに入力されるのは、光検出器で観察される１組のスポット特徴である。各スポット毎に考え得る特徴には、重心、尖度、歪度、分散、ピッチなどが含まれる。この１組の特徴は、ノイズの多い測定（観察）データを含み、ｇと表す。ｇは光検出器によって収集される像から計算し得ることに留意されたい。特徴は、図６に示すように光線追跡図を計算し、これらを用いて検出器上の放射照度パターンを計算することによって、波面パラメータベクトルθからも計算し得る。この場合、検出器で得られる波面パラメータの１組の特徴を
【００３８】
【数３】
0005564545.tif 000004
と表す。この表記を用いてＭＬ推定を式（２）によって数学的に表現し得る。
【００３９】
【数４】
0005564545.tif 000005
・・・（２）
ここで、
【数５】
0005564545.tif 000006
はθのＭＬ推定値を示し、
【数６】
0005564545.tif 000007
【００４０】
で表される量は特徴ベクトルｇを測定する確率である。特徴ベクトルｇは、「実際の」特徴ベクトルが
【数７】
0005564545.tif 000008
であるという仮定で条件付けられる。
【００４１】
上記の最大化問題を解くのに、共役勾配探索、最急降下探索、疑似焼きなまし法、収縮グリッド探索など多くの効率的なアルゴリズムが利用可能である。例として、収縮グリッド探索の詳細な説明を、この探索法が前方伝播モデルに適用される場合について以下に示す。収縮グリッド探索アルゴリズムは、当業者なら容易に実施し得るものである。可能な応用例として、補正エレメント５０４の前の様々な面（検査中の光学エレメントの射出面など）における波面パラメータを、対象とする面までレンズレットアレイにおける場を後方伝播させることによって推定することができる。つまり、本発明のこの態様では、検査エレメント５０３によって生じる波面収差は、像面から補正エレメント５０４と検査エレメント５０３の間の面まで波面パラメータを後方伝播させることによって推定し得る。後方伝播は数値計算によって実施される。このような数値計算は光線追跡またはフーリエ光学技術を用いて実施し得る。ただし、後方伝播を実施する際、後方伝播アルゴリズムの正確さおよび処理時間を慎重に解析すべきである。処理時間の検討については、式（５）〜（１５）を参照して以下で論じる。後方伝播アルゴリズムの例が参照文献４：Ｋｒｉｅｇら、「Ｉｎｖｅｒｓｅ  ｐｒｏｐａｇａｔｉｏｎ  ａｌｇｏｒｉｔｈｍ  ｆｏｒ  ａｎｇｓｔｒｏｍ  ａｃｃｕｒａｃｙ  ｉｎｔｅｒｆｅｒｏｍｅｔｅｒ」、Ｉｎｔｅｒｆｅｒｏｍｅｔｒｙ  ＸＩＩ：Ｔｅｃｈｎｉｑｕｅｓ  ａｎｄ  Ａｎａｌｙｓｉｓ、５５３１巻、Ｐｒｏｃｅｅｄｉｎｇｓ  ｏｆ  ｔｈｅ  Ｓｏｃｉｅｔｙ  ｏｆ  Ｐｈｏｔｏ−Ｏｐｔｉｃａｌ  Ｉｎｓｔｒｕｍｅｎｔａｔｉｏｎ  Ｅｎｇｉｎｅｅｒｓ（ＳＰＩＥ）、Ｃｒｅａｔｈら監修、２００４に記載されている。
【００４２】
上記で論じたように、非干渉測定センサ、特にシャック−ハルトマン波面センサは一般に、大きな合焦／非合焦項（大きな収差を含む波面）を伴う光学エレメントの検査には効果的ではないことに留意されたい。というのは、波面偏差がこうしたセンサのダイナミックレンジを越えやすいからである。上記の説明では、適切に設計された補正エレメントを検査エレメントと従来のＳＨＷＦＳの間に配置して使用することによって、解析中の波面の大きな合焦／非合焦または他の大振幅収差を補正エレメントが効果的に、少なくとも部分的に補償することができ、残った小振幅の収差を解明するだけであることが示された。さらに、補正エレメントをレンズレットアレイに接触させることは厳密には必要とされないので、補正エレメントとレンズレットアレイの間の隙間は、特定の推定問題に応じて、例えば、解析中の波面の収差の予想されるレベルに基づいて変更することができる。
【００４３】
一実施形態では、波面をゼルニケ多項式で表すことができ、検査エレメントの出射面の検出器側の所定の面からの前方伝播モデルを用いてゼルニケ係数を最適化することによって波面パラメータを推定し得る。所定の面の光源側の他の面において波面を決定しなければならない場合、前記所定の面からの後方伝播によって波面パラメータを推定し得る。
【００４４】
波面パラメータの推定には最尤推定（ＭＬＥ）技術を用いる。ＭＬＥは周知の推定方法であるが、高次元問題に実際に用いるには演算集約的過ぎると考えられている。上記で論じたように、参照文献２にＭＬＥ技術の一例が詳細に論じられている。別の例が参照文献３：Ｊａｅ  Ｍｙｕｎｇ、「Ｔｕｔｏｒｉａｌ  ｏｎ  Ｍａｘｉｍｕｍ  Ｌｉｋｅｌｉｈｏｏｄ  Ｅｓｔｉｍａｔｉｏｎ」、Ｊｏｕｒｎａｌ  ｏｆ  Ｍａｔｈｅｍａｔｉｃａｌ  Ｐｓｙｃｈｏｌｏｇｙ、４７巻、９０〜１００（２００３）に提供されている。以下で説明するように、本願では、収縮グリッド法を用いることによってＭＬＥを波面推定で実用的なものにすることを紹介する。
【００４５】
一般に、ＭＬＥ技術を用いて推定を行うには、推定するパラメータが与えられると仮定して、データの確率を知ることが必要とされる。確率分布関数（ＰＤＦ）はｐｒ（ｇ｜θ）で示され、ここで、ｇはＭ×１ベクトルによって特徴付けられる観察データベクトルであり、θはＰ×１ベクトルによって特徴付けられるパラメータベクトルである。
【００４６】
θの尤度は式（３）によって表される。
【００４７】
Ｌ（θ｜ｇ）＝ｐｒ（ｇ｜θ）
・・・（３）
θの最尤（ＭＬ）推定値は、観察データｇが与えられると仮定して、式（４）によって定義される。
【００４８】
【数８】
0005564545.tif 000009
・・・（４）
【００４９】
実際には、解析的にＭＬ推定値を得ることは通常可能でなく、数値的にＭＬ推定を実行し得る。尤度の数値最適化は、非線形最適化法において反復的に実施される。適切なパラメータ（本例ではゼルニケ係数）を見つけるために、パラメータの値をそれらが収束する方向に変化させながら計算を繰り返し実行し得る。推測対象の
【００５０】
【数９】
0005564545.tif 000010
【００５１】
の計算を前方モデルと呼ぶことがある。この数値最適化では、いくつかのモデルは、それらのアルゴリズムにおいて導関数を用いる。例えば、共役勾配探索法または最急降下探索法が最適化アルゴリズムによく用いられる。
【００５２】
導関数を用いる計算にかかる演算時間は極めて長くなることがあり、リソース集約的であるため、導関数が不要なアルゴリズムがより効率的であり、そのため好ましい。
【００５３】
導関数が不要なアルゴリズムがいくつか提案されている。収縮グリッド探索アルゴリズムは、並列処理用に単純かつ効率的なアルゴリズムとして提案されたものである。例えば、参照文献５：Ｌ．Ｒ．Ｆｕｒｅｎｌｉｄ、Ｊ．Ｙ．Ｈｅｓｔｅｍａｎ、およびＨ．Ｈ．Ｂａｒｒｅｔｔ、「Ｆａｓｔ  ｍａｘｉｍｕｍ−ｌｉｋｅｌｉｈｏｏｄ  ｅｓｔｉｍａｔｉｏｎ  ｍｅｔｈｏｄｓ  ｆｏｒ  ｓｃｉｎｔｉｌｌａｔｉｏｎ  ｃａｍｅｒａｓ  ａｎｄ  ｏｔｈｅｒ  ｏｐｔｉｃａｌ  ｓｅｎｓｏｒｓ」、Ｐｒｏｃ．ＳＰＩＥ、６７０７、２００７を参照されたい。収縮グリッド探索アルゴリズムによれば、所望の精度レベルによって決まる所定数の反復において、関数の最大（または最小）値が特定される。この探索法における最適化プロセスが、参照文献５に詳細に説明されているが、ここでは便宜上、それを以下のようにまとめる。このアルゴリズムの実施は、複数のパラメータθからなるＭ×１ベクトルを含むＮ次元空間の探索に適用される。各パラメータ毎に、
１．物理的に妥当なパラメータ空間の周囲の領域を、データベクトル値にかかわらず、無視し得るほど小さい（ほぼゼロに等しい）尤度をもたらす面積（パラメータの数、したがって探索の次元数によっては体積または超体積）で埋める。こうすると、ステップ３ａ（下記）での検査場所に、許容し得る推定値ドメインの外側を、条件検査（ソフトウェアにおける「ｉｆ／ｔｈｅｎ文」）を用いずに、最大尤度として戻される可能性なしに探査させることができる。（あるいは、実行速度をわずかに犠牲にして、物理的に妥当なパラメータ空間内の場所のみが検査されるように条件文を用いる。）
２．較正データ（または予備知識）のみによって決まる尤度表現におけるすべての項を事前に演算する。
【００５４】
３．ｉについて、１から所望の精度で推定値を解くのに必要とされる反復回数まで、
ａ．被検査パラメータ領域を覆う正方形（立方体または超立方体）の中心として定義される検査場所における尤度または対数尤度を演算し、
ｂ．次の反復の検査領域の中心として、尤度または対数尤度が最大になる（または最小限に抑える場合には、最小になる）場所を選択し、
ｃ．各次元における収縮率（例えば、１／２）によって検査領域のサイズを小さくする。
【００５５】
図９に、収縮グリッド探索法のグラフィカルな例を、ＩからＶＩまでの漸進ステップで示す。ここでの目標は、様々なゼルニケ多項式の係数のどの集合が最も大きな（最大）尤度を有しているかを判断することである。図９では、ステップ（Ｉ）で、区域９２０（探索空間）内で、１６個の点９１０からなるグリッドが所定の探索範囲ε内に設定されている。点９１０の数および間隔は、所望の精度レベルに基づいて定義し得る。初期探索範囲ε内に点９１０を定義した後で、これら１６個の点において尤度を演算する。尤度が最も大きい（最大である）適切な点を選択する。ステップ（ＩＩ）で、ステップ（Ｉ）で見つけた最大値を中心として新たなグリッドを定義する。新たなグリッド内の点の間隔（収縮率）は、前のグリッドの間隔の半分とする。図９に示すように、尤度が最も大きい（最大である）適切な点を含む探索範囲を、段階（Ｉ）から段階（ＶＩ）まで、収縮率によって徐々に狭める。
【００５６】
ここで、この収縮グリッド探索法にかかる演算時間を考える。そのために、以下の表記を定義する。
Ｐ：最適化されるパラメータの数
εｐ：最適化を行う初期探索範囲、ｐ＝｛１，２，．．．，Ｐ｝
ｍｐ：各パラメータの更新（尤度分布の計算）後の縮小比（縮小率）、ｐ＝｛１，２，．．．，Ｐ｝
Ｎｐ：尤度分布の計算の反復回数（図９ではＮｐ＝６）
ｌｐ：探索区域の１辺上の尤度分布計算点数、ｐ＝｛１，２，．．．，Ｐ｝（図９ではｌ＝４）
導関数不要の最適化方法では、演算時間の大部分が前方モデルの計算に費やされる。したがって、ＭＬＥのための全演算時間ｔＭＬＥは、式（５）によって以下のように推定される。
【００５７】
ｔＭＬＥ＝ｔｆＮｔｏｔａｌ
・・・（５）
ここで、Ｎｔｏｔａｌは前方モデル計算数の合計であり、ｔｆは１回の前方モデル計算にかかる演算時間である。
【００５８】
すべてのパラメータを同時に最適化する場合、前方モデルの計算点の合計数Ｎｔｏｔａｌは式（６）によって与えられる。
【００５９】
Ｎｔｏｔａｌ＝Ｎｍａｘｌｐ
ただし（Ｎｔｏｔａｌ＝Ｎｍａｘｌ１ｌ２．．．ｌｉ．．．ｌｐ）
・・・  （６）
ここで、すべてのｌｐがｌに等しく、すべてのＮｐがＮｍａｘに等しいと考える。ｌ＝４、Ｐ＝９、Ｎｍａｘ＝６、およびｔｆ＝１０秒の場合、ｔＭＬＥは式（７）によって計算される。
【００６０】
ｔＭＬＥ＝ｔｆＮｍａｘｌｐ＝１０×６×４９／３６００（時間）≒４３６９（時間）
・・・（７）
この時間（約４３６９時間）は、ほとんどの応用例で非現実的であり得る。したがって、収縮グリッド探索技術の場合、探査空間の次元が小さいことが好ましい。
【００６１】
例えば、図９では、Ｎｔｏｔａｌ＝６×４２＝９６点である。同時に最適化しない場合のＮｔｏｔａｌを考えると、パラメータの数Ｐを群の数Ｇによって分割し得る。各群のパラメータ数はｋｇ＝｛ｋ１，ｋ２，．．．，ｋＧ｝であり、各群毎の反復回数はＮｇ＝｛Ｎ１，Ｎ２，．．．，ＮＧ｝である。こうすると、Ｎｔｏｔａｌは式（８）になる。
【００６２】
【数１０】
0005564545.tif 000011
・・・（８）
ここでｋ０＝０
ｌｇ＝｛ｌ１，ｌ２，．．．，ｌＧ｝の場合、式（８）は式（９）になる。
【００６３】
【数１１】
0005564545.tif 000012
・・・（９）
各群毎の反復回数がＮｍａｘである場合、式（８）は式（１０）になる。
【００６４】
【数１２】
0005564545.tif 000013
・・・（１０）
各パラメータ毎の順次最適化では、それぞれの最適化における反復回数はＮｐになる。式（１０）を用いて、Ｎｔｏｔａｌは式（１１）になる。
【００６５】
【数１３】
0005564545.tif 000014
・・・（１１）
ここで、
【００６６】
【数１４】
0005564545.tif 000015
・・・（１２）
各パラメータにおいてＮｐおよびｌｐが一定の場合、式（８）を簡略化することができ、それによって式（１３）が得られる。
Ｎｔｏｔａｌ＝ＮｍａｘＰｌ
・・・（１３）
「ｎ」個の並列プロセッサを使用し得る場合、数ｎＣＰＵを考えることができる。そうすると、式（９）を用いて、各プロセッサ毎の前方モデル計算の合計数は式（１４）になる。
【００６７】
【数１５】
0005564545.tif 000016
・・・（１４）
式（１４）は、各プロセッサ毎の前方モデル計算数が減少し得ることを示している。例えば、９個のパラメータ（すなわちＰ＝９）は３つの群に分けられ、各群は３つのパラメータを有する。
【００６８】
この場合、Ｇ＝３およびｋｇ＝｛３，３，３｝である。式（５）および式（１４）を用いて、ｌ＝４、Ｐ＝９、Ｎｍａｘ＝６、およびｔｆ＝１０秒の場合、ｔＭＬＥは式（１５）になる。
【００６９】
【数１６】
0005564545.tif 000017
・・・（１５）
したがって、並列処理を用いる収縮グリッド探索アルゴリズムは、探索空間の次元が小さい場合には効率的な探索方法となり得る。例えば、この最適化方法において、４つのパラメータを以下のように最適化し得る。
【００７０】
まず、パラメータａの集合により一群のパラメータＡが形成されるとし、Ａ＝｛ａ１，ａ２，ａ３，ａ４｝とする。
【００７１】
次いで、パラメータの集合Ａを２つの群Ａ１＝｛ａ１，ａ２｝とＡ２＝｛ａ３，ａ４｝に分割する。Ａ１群のパラメータおよびＡ２群のパラメータは同時に最適化し得る。パラメータの初期値は、例えば、｛ａ１，ａ２，ａ３，ａ４｝＝｛０，０，０．５，０｝と設定される。真値はそれぞれ｛ａ１，ａ２，ａ３，ａ４｝＝｛０，０，０，０｝である。
【００７２】
正規化コスト関数は、対数尤度関数の否定に対応する。コスト関数の最小値を与えるゼルニケ係数は最大尤度推定値である。コスト関数は、光検出器の各画素における実験データと計算データの差の二乗に関係している。本実施形態では、図示の便宜上、光検出器が含むセンサアレイは、１２８０×１０２４画素の行列を形成すると仮定する。｛ａ１，ａ２｝の探索空間におけるコスト関数面プロファイルが図１０に示すものになると、｛ａ１，ａ２｝の最適値は、収縮グリッドアルゴリズムなどの最適化方法を用いて｛−０．２５，０．０｝に収束する。
【００７３】
その後、｛ａ３，ａ４｝も、図１１に示すように｛０．２５，０．０｝と最適化される。その結果、｛ａ１，ａ２，ａ３，ａ４｝は｛−０．２５，０，０．２５，０｝と最適化される。
【００７４】
最適化された値と真値の間には差が存在する。この前方モデルが｛ａ１，ａ２，ａ３，ａ４｝の線形関数である場合、この最適化結果はほぼ最適である。
【００７５】
しかし、典型的には、前方モデルは非線形関数である。この場合、前方モデルは｛ａ１，ａ３｝の非線形関数である。そのため、｛ａ１，ａ３｝の組合せのコスト関数面は図１２に示すものになる。この例では、ａ３の初期値は０．５であり、これは、真値０．０からの偏りが０．５であることを示している。したがって、｛ａ１，ａ２｝がａ３＝０．５の状態で最適化されると、ａ１は−０．２５に最適化され、これが偏りになる。その結果、最適化点は、真のターゲット値と異なる極小値に向かって収束するように見える。
【００７６】
本発明の第２実施形態では、真のターゲット値と異なる極小値に向かって収束する最適化パラメータの問題の解法を提案する。それによって、推定の正確さが改善される。
【００７７】
図１３に、本発明の第２実施形態による波面測定装置１１１１の構成を示す。波面測定装置１１１１は、光源１１００、波面センサ１１１５、およびデータ解析部１１１８を備える。光源１１００は光ビームＢを発し、この光ビームはフィルタリング光学素子１１０１（例えばＮＤフィルタ）を通過し、照明光学素子１１０２によって整形される。次いで、整形されたビームは検査エレメント１１０３を照明する。検査エレメント１１０３は、レンズ、プリズム、導波路などの単一の光学エレメントとすることもできるし、光学エレメントの組合せとすることもできる。検査エレメント１１０３によって導入される波面収差を補償（低減）するために、レンズレットアレイ１１０４と検査エレメント１１０３の間に（第１実施形態で説明した）補正エレメント５０４を任意選択で配置し得る。
【００７８】
この光ビームは、検査エレメント１１０３を通過し、レンズレットアレイ１１０４に進み、次いで、センサアレイ１１０５を含む光検出器１１０９に至る。センサアレイ１１０５は、この光学系によって生成される光強度分布を検出する。
【００７９】
センサアレイ１１０５を備えた光検出器１１０９は、この光学系全体の情報を有する信号を出力し得る。光検出器１１０９からの出力は、既知の方法で、例えば、データ取得系１１２０のアナログデジタル（ＡＤ）コンバータを用いてデジタル化し得る。次いで、デジタル化されたデータはデータ解析部１１１８によって解析される。データ解析部１１１８は、検査エレメント１１０３を含むこの光学系に関する様々な光学パラメータを解析する。データ解析部１１１８は、コンピュータ制御下でデータを取得し記憶する。データ解析部１１１８では、特殊なソフトウェアを用いてデータフィッティングを行い、様々な計算を行い、データをグラフィカルな形式で表示する。
【００８０】
光源１１００として、パルスレーザまたは連続波（ＣＷ）レーザなどのコヒーレント光源、あるいは水銀ランプまたはタングステンランプなどの非コヒーレント光源を使用し得る。ＣＷレーザまたはランプの場合、光ビームを変調するために光チョッパなどの変調装置が必要とされることがある。光ビームを整形する照明光学素子は当業者には周知であると考えられ、本明細書では説明しない。レンズレットアレイ１１０４については、（１次元および２次元の）様々なレンズレットアレイが当技術分野で周知の商業供給業者から容易に入手可能である。アレイ内のレンズレットの形状、ピッチ、および焦点距離などの特定のパラメータは一般に、測定装置の特定の用途によって決まる。同様に、センサアレイ１１０５は、当技術分野では周知のように、単一のＣＣＤカメラとして実装し得る。あるいは、センサアレイ１１０５は、例えば大規模集積（ＶＬＳＩ）電子技術を用いて専用に設計し得る１次元または２次元のＣＭＯＳセンサアレイとして実装し得る。
【００８１】
図１４に、データ解析部１１１８の概略ブロック図を示す。データ解析部１１１８は、以下に示すステップＳ１１〜Ｓ１８を含む波面推定プロセスを実行し得る。そのために、データ解析部１１１８は、較正モジュール１１０７、パラメータ推定モジュール１１０８、およびハードディスクドライブまたはポータブルメモリデバイスなどの記憶ユニット１１１２（メモリ）に記憶されるデータベース１１１０を備える。データ解析部１１１８は、ステップＳ１１〜Ｓ１８のプロセスシーケンスを実行するように特に構成（プログラム）し得る一般的なプログラマブルコンピュータまたはデータ処理ユニット（マイクロプロセッサ）として実装し得る。データ解析部１１１８およびその構成エレメント（モジュール）は、個々のハードウェアエレメントまたは分散処理ユニットとして実装することもできる。本実施形態では、中央処理装置（ＣＰＵ）などのデータ処理ユニットを含むプログラマブルコンピュータとしてデータ解析部１１１８を実装し得る。一般的なプログラマブルコンピュータの代わりに、あるいはそれに加えて、専用のマイクロプロセッサまたは分散型コンピュータネットワークなどの他のハードウェアデバイスによって、ステップＳ１１〜Ｓ１８のプロセスシーケンスを実現し得る。
【００８２】
データ解析部１１１８は、図１３の波面測定装置の動作中に図１４に示す波面推定プロセスを実施する。ステップＳ１１で、レンズレットパラメータおよび測定された照明分布が手作業または自動的に記憶ユニット１１１２に入力される。レンズレットパラメータは、レンズレットアレイ１１０４の焦点距離およびピッチ、レンズレットアレイ内のレンズレットの数などの光学情報を含み得る。測定された照明分布は、センサアレイ１１０５からデータ取得系１１２０を介して得ることもできるし、局所記憶ユニット１１１２または遠隔記憶装置（図示せず）内のあらかじめ記憶した測定値（例えば、最後に測定した照明分布）から提供されることもある。次いで、入力されたレンズレットパラメータおよび照明分布は、較正モジュール１１０７に提供される。較正モジュール１１０７の機能には、レンズレットアレイと検出器アレイの位置関係についての情報を得る機能、あるいは、これらのアレイの位置を適切にアライメント（較正）する機能がある。つまり、以下でさらに説明するように、較正モジュール１１０７を用いて、例えば当技術分野では周知のフィードバック回路を介して実際の光学ユニットを対話形式でアライメントし得る。これに加えて、またはこの代わりに、較正モジュールを用いて、推定プロセス用の入力情報として使用し得る情報（データ）として光学ユニット間の位置関係を得ることもできる。そのために、ステップＳ１２で、レンズレットアレイおよび他の光学エレメントの位置が、検査エレメント１１０３なしで生成される参照データを用いて定義される。それによって、センサアレイ１１０５は、スクリーンショット１１２１に示すようなほぼ平面波のプロファイルを検出する。データ取得系１１２０は、この参照データをデータ解析部１１１８の較正モジュール１１０７に提供する。これにより、レンズレットアレイの位置を含めて、測定装置の初期光学構成がステップＳ１２で設定される。較正モジュール１１０７によって初期光学構成が設定されると、フロープロセスは、パラメータ推定モジュール１１０８で実施されるシーケンスステップに進む。
【００８３】
データ解析部１１１８によって参照データが取得され、測定装置の初期光学構成が確立されると、検査エレメント１１０３が、光源と波面センサの間の光路に位置決めされ測定が行われる。つまり、図１３で、光源１１００が発する光が検査エレメント１１０３に向けられて、この光の波面が検査光学素子を通過し、波面センサ１１１５のセンサアレイ１１０５によって検出される。ステップＳ１３で、検査エレメント１１０３によって形成されると予想される波面を記述する初期パラメータがパラメータ推定モジュール１１０８に設定される。具体的には、初期パラメータは、デフォルト値としてゼロとし得るゼルニケ係数として設定される。言い換えると、ステップＳ１３では、パラメータ推定モジュール１１０８は、ほぼ平面波（ゼルニケ係数の値がゼロ）を表し得る初期パラメータを設定する。ステップＳ１１〜Ｓ１３が完了すると、ステップＳ１４で、ステップＳ１１〜Ｓ１２の少なくとも１つから受け取った値およびステップＳ１３で設定された初期パラメータに基づいて前方モデルが計算される。つまり、ステップＳ１４では、参照データ（検査光学素子なしで得られたデータ）、測定装置の光学構成、および検査光学素子のデータに基づいて、前方伝播モデルが、ステップＳ１３で設定された初期パラメータに基づいて計算される。以下でより詳細に論ずるように、前方伝播モデルとして、小ビームに基づく伝播モデルまたは光線に基づく伝播モデルを使用し得る。最尤推定（式（１）参照）に関しては、ステップＳ１４で計算された前方伝播モデルを複数の推定モデルパラメータからなるベクトルとみなし得る。
【００８４】
前方モデル（シミュレーションデータ）が確立されると、プロセスはステップＳ１５に進む。ステップＳ１５で、データ取得系１１２０は、センサアレイ１１０５から実際の実験データを取得する。具体的には、ステップＳ１５で、図１３の光学装置において、検査レンズ（検査エレメント１１０３）を用いて実際の実験データが測定される。このようにして、図１４のスクリーンショット１１２２に示すように、センサアレイ１１０５が、検査エレメント１１０３によって形成される波面プロファイルを検出する。データ取得系１１２０は、この実際の実験データ（測定データ）をデータ解析部１１１８のパラメータ推定モジュール１１０８に提供する。したがって、最尤推定（式（１）参照）に関しては、光路内に検査エレメント１１０３が配置された状態で得られた測定データを、複数の観察結果からなるベクトル（複数の観察値からなるベクトル）とみなし得る。
【００８５】
ステップＳ１５で、データ解析部１１１８は、実際の実験データ（例えば、データ取得系１１２０のスクリーンショット１１２２に示すもの）と、ステップＳ１４で前方モデルを用いることによって得られたシミュレーションデータとを用いてパラメータの尤度を計算する。つまり、ステップＳ１５での目的は、複数の観察結果（測定データ）からなるベクトルが与えられると仮定して、観察データを生成（と合致）する最も高い尤度（最大尤度）を有する前方モデルにおけるパラメータからなるベクトルを見つけることである。パラメータの最大尤度値は、実際の実験データとシミュレーションデータの差がどのくらい小さくなり得るかを示し得る。つまり、パラメータと前方モデルが、実験データを生成した実際の光学系をどのくらい正確にモデル化し得るかを示している。ここで、尤度は、実際の実験データとシミュレーションデータが収斂（合致）し得るかどうかを決定し得る（図２２、図２６、および図３０ならびに式１６〜１９を参照して以下で論じる）推定アルゴリズムに基づくものとし得る。
【００８６】
そのため、ステップＳ１６で、実際の実験データ（測定データ）とシミュレーションデータ（前方モデル）が収斂したかどうかが検査される。つまり、ステップＳ１６で、前方モデルのゼルニケ係数（パラメータ）が、測定された波面プロファイルを表す同じパラメータに収斂するかどうかが検査される。ステップＳ１６で、計算結果が収斂している場合、ステップＳ１７で推定されたパラメータ（前方モデルによって確立されたパラメータ）が出力される。しかし、計算結果が収斂しなかった場合、ステップＳ１８で推定されたパラメータを変更することができ、ステップＳ１４とＳ１６の間のプロセスがこの結果が収斂するまで反復して繰り返される。具体的には、本発明では、前方モデルが測定データにほぼ収斂するまでゼルニケ係数の値を反復して変更することによって推定されたパラメータを変更し得る。例えば、式（１６〜１９）を参照して以下で論ずるように、前方モデルによって得られる複数のパラメータＭからなるベクトルにおけるデータｇの値を反復して変更し、最大尤度が得られるか検査して複数の観察値θからなるベクトルを求める。実際の実験データと、推定されたパラメータを用いることによって得られるシミュレーションデータとの差が閾値以下のとき、上記検査結果を収斂したとみなし得る。このようにして、初期ゼルニケ係数（理想波面）と、測定データ（観察された波面パラメータ）にほぼ収斂（合致）する推定されたゼルニケ係数との差を数学的に得ることによって、被検査光学エレメントによって生じた波面収差を求めることができる。
【００８７】
データ解析部１１１８の較正モジュール１１０７により、解析部は、光学素子の位置、検出器パラメータ、および照明分布などいくつかの光学パラメータを較正することができる。具体的には、ステップＳ１１でレンズパラメータおよび測定された照明分布が入力されると、較正モジュール１１０７は、既知の電子技術により、光学系を、その中で所望のパラメータが確立されるように較正する。例えば、較正モジュール１１０７は、電子信号を送信して光源１１００を、センサアレイ１１０５において所望の照明分布を得るのに十分な出力で点灯し得る。さらに、較正モジュール１１０７は、検査エレメント１１０３、レンズレットアレイ１１０４、および検出器１１０９の少なくとも１つの所望の位置を示す電子信号を出力し得る。較正モジュール１１０７から出力される較正信号に応答して、検査エレメント１１０３、レンズレットアレイ１１０４、または検出器１１０９は、所望の位置に移動し（またはアライメントされ）て、センサアレイ１１０５において最適な照明分布が得られる。
【００８８】
パラメータ推定モジュール１１０８は、検査エレメント１１０３が光学系内に存在するときの関心の対象である光学（波面）パラメータを推定する推定機能を実施する。これらのパラメータは、検査エレメント１１０３を含む各光学エレメントの位置ならびに検査エレメント１１０３の特定の平面また球面における波面偏差または強度分布を含み得る。
【００８９】
図１４のフロープロセスは、データベース１１１０に記憶されたデータをうまく利用できるように以下のように改変し得る。具体的には、ステップＳ１３で初期パラメータを設定した後で、データ解析部１１１８は（初期パラメータに基づいて）、データベース１１１０内を探索して、以前に確立した前方モデルに関するデータがデータベース内に存在するかを判断し得る。以前に確立した前方モデルを示す有用なデータがデータベース１１１０内に存在する場合、このデータを用いてステップＳ１５で尤度を計算する。一方、データベース１１１０内にこのようなデータが存在しない場合、先に述べたようにステップＳ１４で前方モデルを計算し、計算結果をデータベース１１１０に記憶し得る。
【００９０】
このデータベース１１１０は、以前に計算された複数の前方モデル出力データを記憶し得る。これらの前方モデル出力データは、前方モデルにおいて定義される計算パラメータとともに記憶され、次いで、このデータベースが索引付けされる。同じまたは類似の検査光学素子（検査エレメント）が波面測定装置によって測定されると、データ解析部におけるパラメータ推定プロセスは、索引付けされたデータベース内の類似パラメータを用いて計算されたデータを探索する。データベース１１１０内に類似のデータが見つかると、（ステップＳ１５で）このデータを用いて尤度を計算する。見つからない場合には、新たな前方モデルが被検査光学エレメントに基づいて計算され、次いで、計算されたデータがデータベースに記憶される。類似の検査光学素子が繰り返し測定される間に、データベース１１１０は類似のパラメータを用いてすでに計算されたデータを記憶し得るので、必要な前方計算数を減らすことができる。
【００９１】
言い換えると、繰返し測定中に処理時間が短くなり得る。したがって、波面測定装置における処理を加速することができる。こうして尤度計算が加速されることにより、より短時間でより多くの反復回数が実現し得る。
【００９２】
上記で述べたように、この測定系についてのＭＬＥにおいて前方モデル計算を使用し得る。前方モデルを計算するためにデータ解析部は光検出器からの出力信号を用い、この目的のために光伝播計算が実行される。
【００９３】
例えば、Ｏｐｔｉｃａｌ  Ｒｅｓｅａｒｃｈ  Ａｓｓｏｃｉａｔｅｓ社によって開発された市販ソフトウェアパッケージであるＣＯＤＥ  Ｖ（登録商標）におけるＢＳＰ（ビーム合成伝播、Ｂｅａｍ  Ｓｙｎｔｈｅｓｉｓ  Ｐｒｏｐａｇａｔｉｏｎ）を用いて前方伝播モデルを実施し得る。ＢＳＰは、小ビームに基づく伝播方法に基づくものである。この小ビームに基づく方法は、高速フーリエ変換（ＦＦＴ）に基づく伝播モデルのいくつかの制限を克服することができる。
【００９４】
前方モデル構成の例として、図１５および図１６に光学構成を示す。図１５は、射出瞳平面２２００を用いた前方伝播についての光学構成を示し、図１６は、射出瞳球面２３００を用いた前方モデルについての光学構成を示す。２２０５および２３０５はレンズレットアレイを示し、２２１０および２３１０はＣＣＤイメージセンサなどの光検出器を示す。
【００９５】
射出瞳の平面２２００または球面２３００において波面または複雑な振幅分布を定義する。射出瞳の平面２２００または球面２３００における光線を、図示しない光源から（像面ＩＰとも称する）ＣＣＤイメージセンサ面まで追跡する。瞳における波面特性および被検査光学エレメントの他の光学設計パラメータに基づいて光線の場所および方向が計算される。
【００９６】
小ビームに基づく伝播では、各光線はガウシアンビームなどの細い小ビームである。イメージセンサ面上で多数の小ビームが重ね合わされ、次いで、像面上の強度プロファイルが計算される。データ解析部は、光学系またはその構成エレメントの、対象とする光学パラメータを簡単な光学系を用いて推定し得る。したがって、本実施形態では、従来の測定技術では測定し得なかった大きな収差を含む波面を測定することができる。その結果、簡単な光学構成で、レンズ、レンズユニット、または光学構成エレメントの光学性能を検証し得る。
【００９７】
瞳上で波面が極めて大きな偏差を有している場合、この偏差を表すにはより多くの小ビームが必要である。小ビームの数は演算リソースによって制限されるので、波面偏差量があるレベルを越えると、小ビームに基づく伝播の正確さは急激に悪くなる。
【００９８】
図１７に、前方モデルの光線追跡図を示す。図１７では、光源側から波面センサ側への順で、光線が検査エレメント２５０４およびレンズレットアレイ２５０５を通って進み、ＣＣＤセンサアレイ２５０６によって検出される。所定の面２５００が、検査エレメント２５０４のすぐ隣の光源側位置に確立され、レンズ設計データ（検査エレメントパラメータ）があらかじめ既知であると仮定される。このモデルにおいて検査エレメントのレンズ設計データが定義されると、所定の面２５００における波面偏差を、レンズ設計データを知らずに所定の面が確立される場合よりもはるかに簡単に推定することができる。レンズ設計データおよび所定の面が確立されているという前提に基づいて、前方伝播モデルが計算される。この前方モデルを用いて、データ解析部において波面偏差または最適化されたレンズ設計パラメータが推定される。次いで、これらの推定結果を用いて波面収差を補正し得る。
【００９９】
さらに、対象とする面に対して測定パラメータを外挿（投影）することによって、対象とする任意の光学パラメータを推定することもできる。図１８に、レンズユニット（レンズ群）２６０８を測定するための光学構成を示す。この場合、表面プロファイル、レンズの厚さ、レンズ間距離、偏心、傾き、または屈折率分布などの光学パラメータを、データ解析部で推定するターゲットパラメータとし得る。
【０１００】
例として、図１８に示すレンズユニット２６０８の波面測定を考え得る。表１に示すように、便宜上、フリンジゼルニケ多項式で波面偏差を表すが、他の関数も使用し得る。データ解析部で推定するパラメータは、これらの多項式の係数である。まず、波面の初期パラメータを、図１８のレンズユニット２６０８（検査光学素子または検査エレメント）の光源側の所定の面２６０９において定義する。ＳＨＷＦＳを形成するように配置されたレンズレットアレイ２６１０およびセンサアレイ２６１５を用いてレンズユニット２６０８を通過する光を検出する。初期レンズ設計データおよび所定の面があらかじめ既知なので、所定の面２６０９から対象とする面（表面）まで初期波面パラメータを前方伝播させることによって、レンズユニット２６０８における対象とするパラメータを推定し得る。
表１：図１８の所定の面２６０９において定義され、対象とする面または表面まで前方伝播されるゼルニケ多項式
【０１０１】
【表１】
0005564545.tif 000018
【０１０２】
最初に、所望の面における波面パラメータを推定するプロセスを開始するために、図１４のステップＳ１１に従ってレンズ設計パラメータおよび検出器データを示す情報が入力される。検出器データに関しては、検出器自体の特定のデータ、例えば、ＣＣＤイメージセンサの画素数（例えば１２８０×１０２４）、ＰＳＮＲ（ピーク信号とノイズの比、Ｐｅａｋ  Ｓｉｇｎａｌ  ｔｏ  Ｎｏｉｓｅ  Ｒａｔｉｏ）＝４００、およびＡＤコンバータの８ビットのデジタル化分解能を入力し得る。より重要には、検出器アレイのところで受け取られる特定の測定照明分布を入力し得る。図１９に、図１８に示すように配置された光学系から得られる検出器データの例を示す。検出器データはデータ解析部１１１８に渡される。
【０１０３】
図１４のデータ解析部１１１８は、較正機能およびパラメータ推定機能を有する。較正機能により、光学素子の位置、検出器パラメータ、照明分布などいくつかの光学パラメータが較正される。
【０１０４】
パラメータ推定機能により、波面収差などの光学パラメータが推定され、それによって、対象とする典型的なパラメータである表面形状または波面偏差が求められる。図１８の例では、まず、既知のレンズ設計データに基づいて所定の面２６０９における波面パラメータを推定し得るが、波面伝播を用いることによって、検査エレメント２６０８に含まれる各光学エレメントの位置または形状に対応する任意のパラメータを、対象とする特定の平面または球面における波面偏差または強度分布を推定することによって推定し得る。
【０１０５】
上記で説明したように、確率分布関数（ＰＤＦ）をｐｒ（ｇ｜θ）で表す。ここで、ｇは観察データベクトルであり、Ｍ×１ベクトルとして特徴付けられ、θはパラメータベクトルであり、Ｐ×１ベクトルとして特徴付けられる。上記で説明したように、ｇが与えられると仮定してθの尤度が式（１）によって定義され、ｇが与えられると仮定してθの最大尤度（ＭＬ）の推定値が式（２）によって定義される。
【０１０６】
実際のケースでは、通常、解析的にＭＬ推定値を得ることができないので、ＭＬ推定値は数値的に求めなければならない。
【０１０７】
例えば、ｇが与えられると仮定してθの確率分布関数（ＰＤＦ）がガウシアン分布によって定義され、尤度が式（１６）によって与えられる。
【０１０８】
【数１７】
0005564545.tif 000019
・・・（１６）
ここで、ｇはパラメータデータであり、ｍはベクトルＭの指数であり、
【０１０９】
【数１８】
0005564545.tif 000020
【０１１０】
はパラメータデータの平均値であり、前方モデルを用いて計算され、σ２はノイズによる分散である。したがって、対数尤度を以下のように記述し得る。
【０１１１】
【数１９】
0005564545.tif 000021
・・・（１７）
その結果、ＭＬ推定値は以下のように記述し得る。
【０１１２】
【数２０】
0005564545.tif 000022
・・・（１８）
図１４に示すパラメータ推定プロセスでは、ＭＬ推定値は反復的に計算される。そのために、平均値データ
【０１１３】
【数２１】
0005564545.tif 000023
【０１１４】
が、ある検査場所における尤度を計算するのに必要である。データ
【０１１５】
【数２２】
0005564545.tif 000024
【０１１６】
は、このプロセスの「前方モデルを計算する」ステップで計算され、次いで、対数尤度ｌｎＬ（θ｜ｇ）がこのプロセスの「尤度を計算する」ステップで計算される。ＭＬ推定値
【０１１７】
【数２３】
0005564545.tif 000025
【０１１８】
は最適化探索法によって求められる。そのため、この探索を最適化するために、式（１９）によって定義するように、前方モデルにおいて計算されるパラメータの各値にコスト関数が割り当てられる。
【０１１９】
【数２４】
0005564545.tif 000026
・・・（１９）
【０１２０】
式（１８）は、
【数２５】
0005564545.tif 000027
【０１２１】
がコスト関数の最小値を有するθであることを意味する。
【０１２２】
したがって、コスト関数プロファイル（あるいは、パラメータの数によっては表面または超表面）内で
【数２６】
0005564545.tif 000028
を探索区域Ａ全体にわたって探索する。
【０１２３】
図２０に、ゼルニケ係数のコスト関数の例を示す。図２０では、Ｃ２は、表１における２次のゼルニケ多項式の係数を示す。この場合、所定の面における理想表面には波面偏差は存在せず、そのため、最適値はゼロになる。最適化プロセスでは、解析部がコスト関数面のプロファイル上でこれらの最適係数を探す。例えば、パラメータＣ２が計算されると、他のパラメータを指定値に固定し得る。
【０１２４】
図２１Ａ〜図２１Ｊに、様々なパラメータ（Ｃｎはｎ次ゼルニケ多項式の係数を意味する）の組合せの２次元コスト関数面を示す。本発明者らは本発明で、図２１Ｂに示す｛Ｃ４，Ｃ９｝、図２１Ｆに示す｛Ｃ２，Ｃ７｝、および図２１Ｇに示す｛Ｃ３，Ｃ８｝のコスト関数が密接に相関している（すなわち、Ｃｎが軸対称のケース）ことを見いだした。したがって、これらの組合せは同時に推定しなければならない。
【０１２５】
言い換えると、これらの組合せのコスト関数面におけるこれらの最適パラメータの探索を実施する。これらの組合せのコスト関数を用いて、前方モデル計算数を減らすことができ、そのため、効率的なＭＬＥを実現し得る。
【０１２６】
コスト関数面の対称性を調べて相関するパラメータを見つけるには、以下の手順を用いることができる。
【０１２７】
（１）対称性を調べる区域Ａを決める。図２１Ａ〜図２１Ｊでは、区域Ａは、ある閾値よりも小さいコスト関数を有する係数の幅として定義される。
【０１２８】
（２）区域Ａにおいてコスト関数面の対称性を調べる。対称性は、図２１Ａ〜図２１Ｊに示すコスト関数面を用いて調べることができる。
【０１２９】
（３）同時に最適化すべきパラメータ（係数）の組合せを決める。図２１Ａ〜図２１Ｊでは、これらの組合せは、｛Ｃ４，Ｃ９｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ９，Ｃ１６｝と決定し得る。同時に最適化すべきパラメータの組合せは、上記手順（２）の下で確立される対称性に基づいて決定（判断）し得る。
【０１３０】
（４）その後、データ解析部で最適化を実施する。
【０１３１】
（５）最適化がうまく収斂しない場合、ステップ（４）の最新の出力パラメータを用いてステップ（１）を繰り返す。
【０１３２】
図２２に、コスト関数面の対称性を調べるプロセスのフローチャートの例を示す。ステップＳ２００で、所定の面における理想パラメータが設定される。要するに、理想（平面）波面に対して、例えば、図１７に示す前方モデルの場合、所定の面においてゼルニケ係数をゼロに設定し得る。所定の面における理想パラメータは、例えば図１６に示すように、理想球面波面の関数の形で送り出すこともできる。つまり、理想値は、前方モデルに用いられる光学構成の設計値である。図２０に示し図２１Ａ〜図２１Ｊで解析される対称性を調べるための区域ＡをステップＳ２０１で求める（決定する）。次いで、ステップＳ２０２で、区域Ａ内のコスト関数面の対称性を調べる。調べた対称性に基づいて、ステップＳ２０３で、同時に最適化すべきパラメータの組合せを求める（決定する）。次いで、ステップＳ２０４で初期値として入力される実際の設計パラメータを用いて、ステップＳ２０５で最適化を実行する。理論的には、任意の初期パラメータを受け取ることができるが、実際には、初期パラメータには（例えば、所定の面におけるゼルニケ係数によって定義される）実際の設計値が用いられる。さらに、以前に測定したデータなどの他のパラメータが利用可能な場合、図２２のＳ２０４でこのような以前のデータも初期パラメータとして使用し得る。最後に、ステップＳ２０６で、コスト関数に基づいて最適化されたパラメータが出力される。パラメータが良好に最適化されない場合（Ｓ２０５のＮＯ）、ステップＳ２０１〜Ｓ２０３で、同時の最適化するパラメータの数を増やし得る。
【０１３３】
図２２のステップＳ２０２に対応する手順（２）では、コスト関数の対称性、すなわち、最適化されるパラメータ間の相関の程度を調べるのにいくつかの方法を用い得る。例えば、対称性を調べるために用いるコスト関数の点数は９点と少なくし得る。
【０１３４】
図２３Ａ〜２３Ｂおよび図２４Ａ〜２４Ｂに、最適化される例示パラメータのコスト関数対称性を調べる簡略化された方法を示す。図２３Ａおよび図２３Ｂは、コスト関数最小化プロセスにおいてパラメータを独立に探索し得る対称な場合を示す。Ｐ１１、Ｐ１２、Ｐ１３、Ｐ２１、Ｐ２２、Ｐ２３、Ｐ３１、Ｐ３２、およびＰ３３におけるコスト関数の９点は図２３Ａに示すように定義される。区域Ａ内の各点の位置に基づいて、点Ｐｘｘの位置に対応するコスト値Ｌｘｘをコスト関数に割り当てることができる。こうすると、図２３Ｂで、Ｌ２２＜Ｌ３２＜Ｌ１２、Ｌ２１＜Ｌ３１＜Ｌ１１、およびＬ２３＜Ｌ３３＜Ｌ１３になる。次いで、コスト関数の３つの値からなる組をそれぞれ、図２３Ｂに示すように、放物線関数に対して曲線フィッティングし得る。３つの曲線フィッティングされた放物線関数の最小値は、点線によって示されるパラメータａ１のほぼ同じ値のところで得られる。したがって、これら３つの放物線関数の軸の場所および向きがほぼ同じであり、コスト関数の値を最小値から最大値まで昇順に並べ替えることができる場合、コスト関数のパラメータは対称であると言え、独立に最適化し得る。
【０１３５】
一方、図２４Ａおよび図２４Ｂには非対称な場合を示す。図２４Ｂでは、放物線関数の３つの最小値がａ１に沿った異なる場所で得られる。さらに、区域Ａ内の各点の位置に基づくコスト関数の値から、Ｌ２２＜Ｌ３２＜Ｌ１２、Ｌ１１＜Ｌ２１＜Ｌ３１、およびＬ３３＜Ｌ２３＜Ｌ１３が得られる。非対称の場合には、コスト関数を最小値から最大値まで昇順に並べ替えることができない。そのため、コスト関数の全体の最小値を見つけるには、パラメータａ１およびａ３を一緒に最適化する必要がある。したがって、このコスト関数は所定の区域Ａ内では非対称であると言って間違いない。これらの基準（コスト関数値のランク順序付け）は、対称（非相関）または非対称（相関）の結論を引き出すには十分ではない。しかし、最小値が単独で共通の位置で生じることは、これらのパラメータを独立に探索し得ることを保証するのに十分であり、ランク順序付けと最小値の場所を組み合わせると、対称性のありなしを適切に判断することが保証され得る。
【０１３６】
所定の面にある量の波面偏差が存在する場合、この所定の面における波面のゼルニケ係数を表にすることができ、これを表２に示す。この場合、コスト関数面は図２５Ａ〜図２５Ｊに示すようになる。
表２：所定の面における収差を含む波面のゼルニケ係数
【０１３７】
【表２】
0005564545.tif 000029
【０１３８】
図２５Ａ〜図２５Ｊのコスト関数の面プロファイルは、パラメータ（ゼルニケ係数）の真値が異なるために、図２１Ａ〜図２１Ｊに示すものと異なる。具体的には、図２５Ｉの｛Ｃ４，Ｃ５｝と図２５Ｊの｛Ｃ４，Ｃ７｝のコスト関数は、図２１Ｉおよび図２１Ｊにそれぞれ示す同類のパラメータのコスト関数よりも大きな相関を示している。ゼルニケ多項式の対称性によれば、｛Ｃ４，Ｃ６｝と｛Ｃ４，Ｃ８｝も相関するはずである。したがって、ゼルニケ係数の数が少なく収斂が実現されない場合、｛Ｃ４，Ｃ５｝、｛Ｃ４，Ｃ６｝、｛Ｃ４，Ｃ７｝、および｛Ｃ４，Ｃ８｝の組合せを追加する必要があり得る。
【０１３９】
上記で説明したように、同時に最適化しなければならない組合せは、｛Ｃ４，Ｃ９｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ９，Ｃ１６｝である。さらに、以下の組合せ（｛Ｃ４，Ｃ５｝、｛Ｃ４，Ｃ６｝、｛Ｃ４，Ｃ７｝、および｛Ｃ４，Ｃ８｝）を追加し得る。その結果、｛Ｃ２，Ｃ４，Ｃ７｝、｛Ｃ３，Ｃ４，Ｃ８｝、または｛Ｃ４，Ｃ９，Ｃ１６｝の組合せも、より効率的な収斂が得られるように働き得る。最適化の開始時において、図１９に示すような検出器データは、倍率および位置のせいで、前方モデルの結果と異なることがある。この場合、先端、傾き、および焦点の成分をあらかじめ最適化する必要がある。
【０１４０】
先端、傾き、および焦点は、｛Ｃ２，Ｃ３，Ｃ４｝に対応する。したがって、この場合のゼルニケ係数は、図２６のフローチャートに示す手順を用いて推定し得る。図２６では、ステップＳ２３００で、各パラメータ毎の初期値が設定される。ステップＳ２３０１で、ゼルニケ多項式の係数Ｃ２、Ｃ３、およびＣ４が１回り目の収縮グリッドグローバル探索として推定される。ステップＳ２３０２、Ｓ２３０３、Ｓ２３０４、およびＳ２３０５で、引き続き１回り目の収縮グリッドグローバル探索中に、｛Ｃ５，Ｃ６｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ４，Ｃ９，Ｃ１６｝がそれぞれ推定される。
【０１４１】
１回り目の収縮グリッドグローバル探索後、２回り目〜６回り目の収縮グリッドローカル探索を順次実行し得る。ステップＳ２３０６、Ｓ２３０７、Ｓ２３０８、Ｓ２３０９、Ｓ２３１０、およびＳ２３１１で、｛Ｃ２，Ｃ３｝、｛Ｃ４，Ｃ９｝、｛Ｃ５，Ｃ６｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ９，Ｃ１６｝がそれぞれ推定される。ステップＳ２３１２で、所望の反復回数Ｎが完了したかどうかの判断がなされる。この反復回数が完了していない場合、収縮グリッドローカル探索が繰り返される。この反復回数が完了した場合、ステップＳ２３１３で、尤度が最大である値が推定結果として出力される。
【０１４２】
収斂の検査として、図２７に、ＭＬ推定中のゼルニケ係数の推定値と既知のターゲット値との差を示す。ターゲット値は、表２にまとめられた設計係数を意味する。図２７のデータは以下のようにして得られる。まず、前方モデルにターゲット値（表２のゼルニケ係数）を入力することによって、光分布のシミュレーションデータ（図１９など）を得ることができる。次いで、このシミュレーションデータに基づいて、ＭＬ推定下でゼルニケ係数を推定し得る。最後に、ゼルニケ係数の推定値とターゲット値の差を得ることができる。
【０１４３】
最適化後、最適化された（推定された）係数とターゲット係数の間に差異がある場合には、以下のプロセスを用いることができる。
【０１４４】
図２８に、｛Ｃ７，Ｃ８｝の組合せのコスト関数面を示す。ターゲット係数は｛５．０λ，５．０λ｝であるが、最適化された値は｛５．０８３４λ，３．２３９５λ｝である。このように、コスト関数面上での探索は極小値に収束したことがわかる。具体的には、本発明者らは本発明で、収差が大きい場合には、探索が誤ってグローバル最小値近傍の極小値に収束し得ることを見いだした。
【０１４５】
図２９Ａに、上記の収縮グリッド探索を用いたグローバル最適化後の｛Ｃ３，Ｃ８｝の組合せのコスト関数面を示す。この場合、探索点は、真のターゲット値と異なる極小値に向かって収束し得る。正確さを改善するために、本発明者らは本発明で、収縮グリッド探索法に対するさらなる改善を考案した。具体的には、
（１）下記の式（２０）によって定義される単純な２次関数を、例えば最小二乗法を用いて、コスト関数の離散サンプルに対してフィッティングする。
ｆ（ｘ，ｙ）＝ａ１ｘ２＋ａ２ｙ２＋ａ３ｘｙ＋ａ４ｘ＋ａ５ｙ＋ａ６
・・・（２０）
（２）この２次関数の最小値が存在する場所を見つけ、この最小値が生じる場所を、収縮グリッド探索の次の反復のための新たな中心として用いる。この方法によって、グローバル最小値に向かって収束するように改善されることが示された。
【０１４６】
図２９Ａに｛Ｃ３，Ｃ８｝尤度の離散データを示し、図２９Ｂに｛Ｃ３，Ｃ８｝空間のフィッティングされた面を示す。この場合、｛Ｃ３，Ｃ８｝の値＝｛９．８９９７８λ，４．６４９２６λ｝が、フィッティングされた面から最適係数として得られた。この方法を用いて、図２７に示す最適化手順が大きく改善された。
【０１４７】
図３０に、収縮グリッド探索アルゴリズムについての改変されたフロープロセスを示す。具体的には、図３０は、フィッティング方法を用いて改善されたゼルニケ係数を推定する最適化プロセスを示している。この改善された最適化では、１回り目の収縮グリッドグローバル探索（ステップＳ２７０１〜Ｓ２７０５）と２回り目〜６回り目の収縮グリッドローカル探索（ステップＳ２７１０〜Ｓ２７１５）との間にフィッティングプロセスが実行される。ステップＳ２７００で初期値が設定された後で、ステップＳ２７０１からＳ２７０５で、｛Ｃ２，Ｃ３，Ｃ４｝、｛Ｃ５，Ｃ６｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ４，Ｃ９，Ｃ１６｝のパラメータの組合せに対して収縮グリッドグローバル探索が実行される。
【０１４８】
ステップＳ２７０６からＳ２７０９で、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、｛Ｃ４，Ｃ９｝、および｛Ｃ５，Ｃ６｝に対するフィッティングプロセスが実行される。このフィッティングプロセスの後で、ステップＳ２７１０からＳ２７１５で、｛Ｃ２，Ｃ３｝、｛Ｃ４，Ｃ９｝、｛Ｃ５，Ｃ６｝、｛Ｃ２，Ｃ７｝、｛Ｃ３，Ｃ８｝、および｛Ｃ９，Ｃ１６｝のパラメータの組合せに対して収縮グリッド探索が実行される。ステップＳ２７１６で、所望の反復回数Ｎが完了したかどうかの判断がなされる。この反復回数が完了していない場合、収縮グリッドローカル探索が繰り返される。この反復回数が完了した場合、ステップＳ２７１７で推定値が出力される。図３０で説明した改善されたプロセスは、上記の第１実施形態の後方伝播計算に適用し得る。
【０１４９】
図３１に、このフィッティング方法を用いたＭＬ推定中の、推定されたゼルニケ係数とターゲット係数の差を示す。この結果から、ターゲット係数と推定された係数の差が０．１λよりも小さいくなり得る。したがって、本発明者らは本発明で、フィッティング方法を用いる収縮グリッド探索が収束性能を劇的に改善し得ることを確認した。
【０１５０】
本明細書で説明したように、最適化にかかる演算時間が短くなり、推定結果の正確さが改善されることにより、有利には、ここで説明した原理から逸脱することなく、波面測定技術が改善され得る。
【０１５１】
第１および第２の実施形態の概念は、複数の検出器または検出器アレイから得られるデータの処理に適用可能である。複数の検出器から得られるデータは、検出器位置、照明モード、偏光状態、または検査エレメント位置を変えて取得し得る。
【０１５２】
本発明の少なくとも１つの実施形態により、従来の測定技術では測定し得なかった大きな収差を含む波面を測定することができる。したがって、単レンズやレンズユニットなどの検査エレメントの光学性能の検査または検証を簡単な光学構成で実現し得る。
【０１５３】
第１実施形態では、補正エレメントは、検査エレメントによって生じるある種の大きな低次波面収差を補償することができる。第２実施形態では、収縮グリッド探索のフィッティング方法により、推定値の収束が効率的かつ正確に計算されることが保証される。
【０１５４】
例示の実施形態を参照して本発明を説明してきたが、本発明は開示した例示の実施形態に限定されるものではないことを理解されたい。添付の特許請求の範囲は、最も広い解釈に従うものとし、そのため、あらゆる改変形態ならびに等価な構造および機能が含まれるものとする。
【符号の説明】
【０１５５】
１００  波面測定装置
１０２  光源
１０４  ＮＤフィルタ
１０６  ビームエキスパンダ
１０８  投影レンズ
１１０  検査光学素子
１１２  レンズ群
１１４  レンズレットアレイ
１１６  光検出器
１１８  データ解析部
１５０  波面
１６０  波面
２１０  センサアレイ
２５０〜２５４  焦点スポット
３５０〜３５４  焦点スポット
５００  波面測定装置
５０１  光源
５０２  ビームエキスパンダ
５０３  光学検査エレメント
５０４  補正エレメント
５０５  レンズレットアレイ
５０６  光検出器
５０９  レンズ
５１０  波面センサ
５１８  データ解析部
５２４  平行平面板
１１００  光源
１１０１  フィルタリング光学素子
１１０２  照明光学素子
１１０３  検査エレメント
１１０４  レンズレットアレイ
１１０５  センサアレイ
１１０９  光検出器
１１１１  波面測定装置
１１１５  波面センサ
１１０７  較正モジュール
１１０８  パラメータ推定モジュール
１１１０  データベース
１１１２  記憶ユニット
１１１８  データ解析部
１１２０  データ取得系
１１２１  スクリーンショット
１１２２  スクリーンショット
２２００  射出瞳平面
２２０５  レンズレットアレイ
２２１０  光検出器
２３００  射出瞳球面
２３０５  レンズレットアレイ
２３１０  光検出器
２５００  所定の面
２５０４  検査エレメント
２５０５  レンズレットアレイ
２５０６  ＣＣＤセンサアレイ
２６０８  レンズユニット
２６０９  所定の面
２６１０  レンズレットアレイ
２６１５  センサアレイ
ａ１〜ａ４  パラメータ
Ａ  区域
Ｂ  光ビーム
Ｃ２〜Ｃ１６  係数
ＩＰ  像面
Ｌ１１〜Ｌ３３  コスト値
Ｐ１１〜Ｐ３３  コスト関数の点
Ｓ１１〜Ｓ１８  波面推定プロセスステップ
Ｓ２００〜Ｓ２０６  対称性調査プロセスステップ
Ｓ２３０１〜Ｓ２３１１  係数推定プロセスステップ
Ｓ２７０１〜Ｓ２７０５  係数推定プロセスステップ
Ｓ２７０６〜Ｓ２７０９  フィッティングプロセスステップ
Ｓ２７１０〜Ｓ２７１５  係数推定プロセスステップ
【図１】
0005564545.tif 000030
【図２】
0005564545.tif 000031
【図３】
0005564545.tif 000032
【図４】
0005564545.tif 000033
【図５】
0005564545.tif 000034
【図６】
0005564545.tif 000035
【図７Ａ】
0005564545.tif 000036
【図７Ｂ】
0005564545.tif 000037
【図８】
0005564545.tif 000038
【図９】
0005564545.tif 000039
【図１０】
0005564545.tif 000040
【図１１】
0005564545.tif 000041
【図１２】
0005564545.tif 000042
【図１３】
0005564545.tif 000043
【図１４】
0005564545.tif 000044
【図１５】
0005564545.tif 000045
【図１６】
0005564545.tif 000046
【図１７】
0005564545.tif 000047
【図１８】
0005564545.tif 000048
【図１９】
0005564545.tif 000049
【図２０】
0005564545.tif 000050
【図２１Ａ】
0005564545.tif 000051
【図２１Ｂ】
0005564545.tif 000052
【図２１Ｃ】
0005564545.tif 000053
【図２１Ｄ】
0005564545.tif 000054
【図２１Ｅ】
0005564545.tif 000055
【図２１Ｆ】
0005564545.tif 000056
【図２１Ｇ】
0005564545.tif 000057
【図２１Ｈ】
0005564545.tif 000058
【図２１Ｉ】
0005564545.tif 000059
【図２１Ｊ】
0005564545.tif 000060
【図２２】
0005564545.tif 000061
【図２３Ａ】
0005564545.tif 000062
【図２３Ｂ】
0005564545.tif 000063
【図２４Ａ】
0005564545.tif 000064
【図２４Ｂ】
0005564545.tif 000065
【図２５Ａ】
0005564545.tif 000066
【図２５Ｂ】
0005564545.tif 000067
【図２５Ｃ】
0005564545.tif 000068
【図２５Ｄ】
0005564545.tif 000069
【図２５Ｅ】
0005564545.tif 000070
【図２５Ｆ】
0005564545.tif 000071
【図２５Ｇ】
0005564545.tif 000072
【図２５Ｈ】
0005564545.tif 000073
【図２５Ｉ】
0005564545.tif 000074
【図２５Ｊ】
0005564545.tif 000075
【図２６】
0005564545.tif 000076
【図２７】
0005564545.tif 000077
【図２８】
0005564545.tif 000078
【図２９Ａ】
0005564545.tif 000079
【図２９Ｂ】
0005564545.tif 000080
【図３０】
0005564545.tif 000081
【図３１】
0005564545.tif 000082