(57)【要約】【課題】  特別な測定装置を用いず、光強度分布情報から光位相分布を高速且つ高精度に同定する光位相分布測定方法を提供する。【解決手段】  異なる２つの光学系とＣＣＤとを備えた光位相分布測定系を用いて、被測定光を各光学系にそれぞれ入力し、強度と位相を変調し、出力された被測定光をＣＣＤにより検出し、検出された被測定光の強度分布を画像として測定するステップと、強度分布と各光学系の光学特性とに基づいて観測方程式を設定するステップと、観測方程式から位相分布同定逆問題を設定し、設定された位相分布同定逆問題を被測定光を表す複素振幅を設計変数とする第１の非線形最適化問題として定式化するステップと、被測定光の位相分布を級数展開することにより、第１の非線形最適化問題を級数展開式の展開係数を設計変数とする第２の非線形最適化問題に変換するステップと、第２の非線形最適化問題を解くことで被測定光の位相分布を同定するステップとを有する。【選択図】  図１
2006234389.tif 000002
【特許請求の範囲】
【請求項１】
被測定光の位相分布を高速且つ高精度に測定するための光位相分布測定方法であって、
異なる２つの光学特性が既知の光学系（第１の光学系及び第２の光学系）と光波検出センサーとを備えた光位相分布測定システムを用いて、前記被測定光を前記第１の光学系と前記第２の光学系にそれぞれ入力し、強度と位相を変調し、出力された前記被測定光を前記光波検出センサーにより検出し、検出された前記被測定光の強度分布を画像として測定する強度分布測定ステップと、
強度分布測定ステップで得られた前記強度分布と、前記第１の光学系の光学特性と、前記第２の光学系の光学特性とに基づいて、前記光位相分布測定システムの観測方程式を設定する観測方程式設定ステップと、
前記観測方程式から被測定光の位相分布を同定する位相分布同定逆問題を設定し、設定された前記位相分布同定逆問題を、前記被測定光を表す複素振幅を設計変数とする第１の非線形最適化問題として定式化する第１の非線形最適化問題定式化ステップと、
前記被測定光の位相分布を級数展開することにより、前記第１の非線形最適化問題を、級数展開式の展開係数を設計変数とする第２の非線形最適化問題に変換する第２の非線形最適化問題定式化ステップと、
前記第２の非線形最適化問題を解くことで、前記被測定光の位相分布を同定する位相分布同定ステップと、
を有することを特徴とする光位相分布測定方法。
【請求項２】
2006234389.tif 000136
は前記被測定光を表す複素振幅で、[Ｈ１]は前記第１の光学系の光学特性を表すシステムマトリクスで、[Ｈ１]−１は[Ｈ１]の逆行列で、[Ｈ２]は前記第２の光学系の光学特性を表すシステムマトリクスで、
2006234389.tif 000137
は前記第１の光学系を通過した検出光を表す複素振幅で、
2006234389.tif 000138
は前記第２の光学系を通過した検出光を表す複素振幅で、
2006234389.tif 000139
と記号の入れ換えをした場合に、前記観測方程式は、次の数式で表され、
2006234389.tif 000140
前記第１の非線形最適化問題の目的関数は、
2006234389.tif 000141
で表される請求項１に記載の光位相分布測定方法。
【請求項３】
前記被測定光の位相分布φをゼルニケ多項式に展開し、ゼルニケ係数を前記第２の非線形最適化問題の設計変数とし、前記第２の非線形最適化問題の目的関数は、
2006234389.tif 000142
で表され、ただし、ｚは所定の個数のゼルニケ係数からなるベクトルである請求項２に記載の光位相分布測定方法。
【請求項４】
前記位相分布同定ステップでは、
ゼルニケ係数ｚの初期値を設定するステップ１と、
前記被測定光の位相分布φを
2006234389.tif 000143
に従って計算するステップ２と、
前記被測定光を表す複素振幅
2006234389.tif 000144
に従って計算するステップ３と、
前記被測定光の第１の仮想境界における複素振幅
2006234389.tif 000145
に従って計算するステップ４と、
ステップ４で算出された前記被測定光の第１の仮想境界における複素振幅から光強度を算出し、算出された光強度と前記光波検出センサーにより検出された光強度と比較し、比較した結果が前記第２の非線形最適化問題の目的関数が極小点に達した場合に終了し、一方、比較した結果が前記第２の非線形最適化問題の目的関数が極小点に達しない場合に、ステップ６に進むステップ５と、
前記ゼルニケ係数ｚを修正してから、ステップ２に戻るステップ６と、
を含む請求項３に記載の光位相分布測定方法。
【請求項５】
前記被測定光を表す複素振幅
2006234389.tif 000146
を２次元フーリエ変換することによって、
2006234389.tif 000147
を求めるステップＡ１と、
2006234389.tif 000148
を計算するステップＡ２と、
2006234389.tif 000149
を２次元フーリエ逆変換することによって、前記被測定光の第１の仮想境界における複素振幅
2006234389.tif 000150
を求めるステップＡ３と、
を行うことによって、前記ステップ４で計算される前記被測定光の第１の仮想境界における複素振幅
2006234389.tif 000151
を求める請求項４に記載の光位相分布測定方法。
【請求項６】
前記第１の光学系及び前記第２の光学系に対し、次のキルヒホッフの積分方程式が成立し、
2006234389.tif 000152
仮想境界Γ上の複素振幅ｕを境界条件として、前記キルヒホッフの積分方程式を直接解くことにより、前記被測定光に特別な仮定をおくことなく仮想境界Γ上のｑを導出し、仮想境界Γ上の複素振幅ｕと点Ｐの複素振幅ｕ(Ｐ)の関係を求める請求項１に記載の光位相分布測定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、光位相分布測定方法に関し、特に、簡単且つ経済的に光位相分布を高速且つ高精度に測定できるようにした光位相分布測定方法に関する。
【背景技術】
【０００２】
光位相分布（以下、光の位相分布とも称する）には、光源の波面形状、透過物体の屈折率分布、反射物体の表面形状などの情報が含まれている。光位相分布の計測は、光学システムの設計評価、半導体レーザ素子及びレンズなどの光学素子の品質管理、非接触非破壊検査や顕微鏡などの表面形状計測、光波面の整形や天体観測に利用される波面補償光学など多くの分野で必要とされている。光位相分布の定量的な測定は、工学的に非常に重要である。特に、近年ＤＶＤピックアップに代表されるレーザによる情報書き込み読み出し装置の評価においては、収差解析と呼ばれる位相分布の計測は、必要不可欠な技術である（非特許文献１及び非特許文献２を参照）。
【０００３】
従来より、ハードウェア（つまり、特別な測定装置）を用いた、光位相分布の測定方法として、例えば、干渉計により被測定光と参照光との相対的な位相分布を求める方法（以下、干渉計による光位相分布測定方法）（非特許文献３を参照）や、シャック・ハートマンセンサーにより位相勾配分布を検出し位相分布を求める方法（以下、シャック・ハートマンセンサーによる光位相分布測定方法）（非特許文献４を参照）などがある。
【非特許文献１】ケイ.エイ.ヌージェント（K.A.Nugent）・ディー.パガニン（D.paganin）・ティー.イー.グーレエフ（T.E.Gureyev）共著,「ア  フェーズ  オデッセイ（APhase Odyssey）」,フィジクス  トゥデイ（Physics Today）,第54巻,第8号,p.27-32,2001年
【非特許文献２】三好隆志著,「光応用ナノインプロセス計測技術の最新動向」,機械の研究,第54巻,第9号,p.928-933,2002年
【非特許文献３】「光測定器ガイド」,オプトロニクス社,1998年
【非特許文献４】ダブリュー.エイチ.サウスウェル（W.H.Southwell）著,「ウエーブフロント  エスティメイション  フロム  ウエーブフロント  スロープ  メジャーメントズ（Wavefront estimation from wavefrontslope measurements）」,ジェー.  オプト.  ソク.  エム.（J.Opt.Soc.Am.）,第70巻,第8号,p.998-1005,1980年
【非特許文献５】ピー.スチスケ（P.Schiske）著,「ユニークネス  アンド  アンビギュイティー  イン  フェーズ  デタミネイション  フロム  イメージ  アンド  ディフラクション  パターン  インテンシティーズ（Uniquenessand ambiguity in phase determination from image and diffraction patternintensities）」,オプティク（Optik）,第40巻,第3号,p.261-275,1974年
【非特許文献６】吉村武晃著,「光情報工学の基礎」,コロナ社,2000年
【非特許文献７】矢部・八卷共著,「非線形計画法」,朝倉書店,1999年
【非特許文献８】宮本健郎著,「光学入門」,岩波書店,1995年
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかしながら、上述した従来の光位相分布測定方法に使用されている特別な測定装置（つまり、干渉計やシャック・ハートマンセンサーなどのハードウェア）は、精密な光学機器から構成されているので非常に高価で、また、こういった特別な測定装置を設置するスペースも必要となるという問題があった。さらに、被測定光の位相分布のスケールが波長レベルになると、計測のための調整の労力や計測誤差が増大する問題もあった。
【０００５】
本発明は、上述のような事情よりなされたものであり、本発明の目的は、特別な測定装置を用いることなく、簡便に測定可能な光強度分布の情報から光位相分布を高速且つ高精度に同定することによって光位相分布を測定できるようにした光位相分布測定方法を提供することにある。
【課題を解決するための手段】
【０００６】
本発明は、被測定光の位相分布を高速且つ高精度に測定するための光位相分布測定方法に関し、本発明の上記目的は、異なる２つの光学特性が既知の光学系（第１の光学系及び第２の光学系）と光波検出センサーとを備えた光位相分布測定システムを用いて、前記被測定光を前記第１の光学系と前記第２の光学系にそれぞれ入力し、強度と位相を変調し、出力された前記被測定光を前記光波検出センサーにより検出し、検出された前記被測定光の強度分布を画像として測定する強度分布測定ステップと、強度分布測定ステップで得られた前記強度分布と、前記第１の光学系の光学特性と、前記第２の光学系の光学特性とに基づいて、前記光位相分布測定システムの観測方程式を設定する観測方程式設定ステップと、前記観測方程式から被測定光の位相分布を同定する位相分布同定逆問題を設定し、設定された前記位相分布同定逆問題を、前記被測定光を表す複素振幅を設計変数とする第１の非線形最適化問題として定式化する第１の非線形最適化問題定式化ステップと、前記被測定光の位相分布を級数展開することにより、前記第１の非線形最適化問題を、級数展開式の展開係数を設計変数とする第２の非線形最適化問題に変換する第２の非線形最適化問題定式化ステップと、前記第２の非線形最適化問題を解くことで、前記被測定光の位相分布を同定する位相分布同定ステップとを有することによって効果的に達成される。
【０００７】
また、本発明の上記目的は、
2006234389.tif 000003
は前記被測定光を表す複素振幅で、[Ｈ１]は前記第１の光学系の光学特性を表すシステムマトリクスで、[Ｈ１]−１は[Ｈ１]の逆行列で、[Ｈ２]は前記第２の光学系の光学特性を表すシステムマトリクスで、
2006234389.tif 000004
は前記第１の光学系を通過した検出光を表す複素振幅で、
2006234389.tif 000005
は前記第２の光学系を通過した検出光を表す複素振幅で、
2006234389.tif 000006
と記号の入れ換えをした場合に、前記観測方程式は、次の数式で表され、
2006234389.tif 000007
前記第１の非線形最適化問題の目的関数は、
2006234389.tif 000008
で表されることにより、或いは、前記被測定光の位相分布φをゼルニケ多項式に展開し、ゼルニケ係数を前記第２の非線形最適化問題の設計変数とし、前記第２の非線形最適化問題の目的関数は、
2006234389.tif 000009
で表され、ただし、ｚは所定の個数のゼルニケ係数からなるベクトルであることにより、或いは、前記位相分布同定ステップでは、ゼルニケ係数ｚの初期値を設定するステップ１と、前記被測定光の位相分布φを
2006234389.tif 000010
に従って計算するステップ２と、前記被測定光を表す複素振幅
2006234389.tif 000011
に従って計算するステップ３と、前記被測定光の第１の仮想境界における複素振幅
2006234389.tif 000012
に従って計算するステップ４と、ステップ４で算出された前記被測定光の第１の仮想境界における複素振幅から光強度を算出し、算出された光強度と前記光波検出センサーにより検出された光強度と比較し、比較した結果が前記第２の非線形最適化問題の目的関数が極小点に達した場合に終了し、一方、比較した結果が前記第２の非線形最適化問題の目的関数が極小点に達しない場合に、ステップ６に進むステップ５と、前記ゼルニケ係数ｚを修正してから、ステップ２に戻るステップ６とを含むことにより、或いは、前記被測定光を表す複素振幅
2006234389.tif 000013
を２次元フーリエ変換することによって、
2006234389.tif 000014
を求めるステップＡ１と、
2006234389.tif 000015
を計算するステップＡ２と、
2006234389.tif 000016
を２次元フーリエ逆変換することによって、前記被測定光の第１の仮想境界における複素振幅
2006234389.tif 000017
を求めるステップＡ３とを行うことによって、前記ステップ４で計算される前記被測定光の第１の仮想境界における複素振幅
2006234389.tif 000018
を求めることにより、或いは、前記第１の光学系及び前記第２の光学系に対し、次のキルヒホッフの積分方程式が成立し、
2006234389.tif 000019
仮想境界Γ上の複素振幅ｕを境界条件として、前記キルヒホッフの積分方程式を直接解くことにより、前記被測定光に特別な仮定をおくことなく仮想境界Γ上のｑを導出し、仮想境界Γ上の複素振幅ｕと点Ｐの複素振幅ｕ(Ｐ)の関係を求めることによって一層効果的に達成される。
【発明の効果】
【０００８】
本発明に係る光位相分布測定方法によれば、従来の光位相分布測定方法に用いられる特別な測定装置のような高価な光学機器を用いることなく、２つの異なる光学系と光波検出センサーとで構成された光位相分布測定系を用いて得られた光強度分布の情報（強度分布画像）から、被測定光の位相分布を高速且つ高精度に同定することができ、簡便かつ経済的な光位相分布の測定を可能にしたといった優れた効果を奏する。
【発明を実施するための最良の形態】
【０００９】
以下、本発明を実施するための最良の形態を図面を参照して説明する。
＜１＞光波の数理表現及び光位相分布測定系
まず、本発明で使用される光波の数理表現について説明する。
【００１０】
本発明において、位相分布測定の対象とする光（つまり、被測定光）が、例えばレーザ光などに代表されるような、単一の波長成分を持つ単色光で、且つ、異なる時刻における位相が常に一定の相関関係を持つ時間的にコヒーレントな光であることを前提とする。
【００１１】
本発明では、光学における計算において、光波ｕの時間振動項を分離し、複素振幅
2006234389.tif 000020
で光波の強度
2006234389.tif 000021
と位相θを下記数１で表現する（非特許文献６を参照）。
【００１２】
【数１】
2006234389.tif 000022
ここで、
2006234389.tif 000023
は虚数を表し、ωは角周波数を表し、ｔは時刻を表す。
【００１３】
次に、本発明において使用される、光位相分布測定系について説明する。
【００１４】
図１は、本発明で用いる光位相分布測定系の測定原理を説明するための模式図である。図１に示されるように、光位相分布測定系は、光学特性が既知である２つの異なる光学系と光波検出センサーとから構成される。光位相分布測定系の光波検出センサーとしては、ＣＣＤ撮像素子を用いるのが好ましいが、それに限定されることなく、他のセンサーを用いても良い。なお、本実施形態において、光波検出センサーとしては、ＣＣＤを用いる。
【００１５】
本発明では、このような光位相分布測定系を用いて、光波検出センサー（ＣＣＤ）により被測定光の強度分布画像を測定し、そして、得られた被測定光の強度分布画像や光学系の光学特性などの情報に基づいて、後述する光位相分布測定方法を利用して、被測定光の位相分布を同定することによって、被測定光の位相分布を測定できるようにしている。
【００１６】
前述したように、本発明では、単一の波長成分を持つ単色光を光位相分布の測定対象（つまり、被測定光）としているが、例えば、被測定光に注目する波長成分以外の波長をもつ光波が含まれてしまう場合には、光波検出センサー（ＣＣＤ）の前にフィルターを置き、注目する波長成分のみを取りだすことによって、本発明を適用することができる。
【００１７】
以下、図１を参照しながら、光位相分布測定系の測定原理及び数理モデルについて説明する。図１に示されるように、光位相分布測定系は、光学特性が既知である異なる２個の光学系（光学系１及び光学系２）と光波検出センサー（ＣＣＤ）とから構成されている。
【００１８】
光位相分布測定系の測定原理を数学的に表現すると、つまり、被測定光を表す複素数ベクトル（つまり、被測定光の複素振幅で、以下、単に、複素振幅とも称する）
2006234389.tif 000024
が光学系１による変換マトリックス[Ｈ１]により変換され、光学系１を通過した検出光を表す複素数ベクトル
2006234389.tif 000025
となり、また、被測定光を表す複素数ベクトル
2006234389.tif 000026
が光学系２による変換マトリックス[Ｈ２]により変換され、光学系２を通過した検出光を表す複素数ベクトル
2006234389.tif 000027
となることを意味する。要するに、被測定光と検出光との関係は、後述する数２及び数３で表すことができる。
【００１９】
ここで、光波検出センサー（ＣＣＤ）により検出されるのは、検出光を表す複素数ベクトル
2006234389.tif 000028
のそれぞれの成分の絶対値を並べた実数ベクトルである。
【００２０】
要するに、ＣＣＤなどの光波検出センサーにより受光検出され画像として測定されるのは、光の強度、すなわち、複素振幅の絶対値である。測定された光の強度の空間分布は、二次元の離散配列で与えられるが、本発明では、光波を表す複素振幅の分布を取扱う場合に、離散的に一次元に並べ換え複素数ベクトルで表現する。
【００２１】
なお、本発明では、光位相分布測定系に適用できる光学系は、光学特性が既知であれば、どのようなものでも構わない。簡単に利用できるものとして、例えば、レンズ、レンズの焦点はずれ、瞳関数の異なるレンズ、回折格子、空気中の回折（即ち、光学系を介さず直接入射させること）などの光学系が挙げられる。
＜２＞光位相分布測定系の観測方程式
入力光（つまり、被測定光）
2006234389.tif 000029
をｎ個に離散化し、
2006234389.tif 000030
、異なる２つの光学系（光学系１及び光学系２）を通過後の光をそれぞれｍ個に離散化し、
2006234389.tif 000031
、光学系１及び光学系２による変換マトリクスをそれぞれ
2006234389.tif 000032
とする。また、光波検出センサーであるＣＣＤにより測定された検出光の強度分布(振幅の２乗）を
2006234389.tif 000033
とすると、光位相分布測定系の観測方程式は、下記数２、数３、数４及び数５のように書くことができる。
【００２２】
【数２】
2006234389.tif 000034
【００２３】
【数３】
2006234389.tif 000035
【００２４】
【数４】
2006234389.tif 000036
【００２５】
【数５】
2006234389.tif 000037
ただし、数４と数５の形式で表現されるベクトル
2006234389.tif 000038
は、複素数ベクトル
2006234389.tif 000039
の絶対値を並べた実数ベクトルを表すものとする。
【００２６】
【数６】
2006234389.tif 000040
上記数２は光学系１による被測定光
2006234389.tif 000041
から検出光
2006234389.tif 000042
への変換式で、上記数３は光学系２による被測定光
2006234389.tif 000043
から検出光
2006234389.tif 000044
への変換式である。また、上記数４は検出光を表す複素数ベクトル
2006234389.tif 000045
の絶対値ベクトル
2006234389.tif 000046
が測定された強度
2006234389.tif 000047
として与えられる、つまり、測定量であることを表す。上記数５は、検出光を表す複素数ベクトル
2006234389.tif 000048
の絶対値ベクトル
2006234389.tif 000049
が測定された強度
2006234389.tif 000050
として与えられる。
【００２７】
ここで、既知量は、光学系１の光学特性を表す複素数マトリックス[Ｈ１]、光学系２の光学特性を表す複素数マトリックス[Ｈ２]、及び測定された検出光の強度を表す絶対値ベクトル
2006234389.tif 000051
であり、未知量は、被測定光を表す複素数ベクトル
2006234389.tif 000052
と、検出光を表す複素数ベクトル
2006234389.tif 000053
である。
＜３＞非線形最適化問題への帰着
上記数２、数３、数４及び数５で表す、光位相分布測定系の観測方程式は、下記数７、数８、数９及び数１０のように簡便化することができる。
【００２８】
【数７】
2006234389.tif 000054
【００２９】
【数８】
2006234389.tif 000055
【００３０】
【数９】
2006234389.tif 000056
【００３１】
【数１０】
2006234389.tif 000057
ここで、[Ｈ１]−１は光学系１のシステムマトリクス[Ｈ１]の逆行列である。なお、光学系のシステムマトリクス（光学系の光学特性を表す複素数マトリックス）の逆行列がない場合には、例えば、一般逆行列を導入する。また、被測定光
2006234389.tif 000058
の強度が観測されている場合には、[Ｈ１]は単位行列となる。まず、数８、数９及び数１０を解いて
2006234389.tif 000059
を求め、次に、数７を用いて
2006234389.tif 000060
を得ることができる。
【００３２】
ここで、数７、数８、数９及び数１０で表す観測方程式に対し、
2006234389.tif 000061
と記号の入れ換えをすると、下記数１１、数１２及び数１３で表す観測方程式が得られる。
【００３３】
【数１１】
2006234389.tif 000062
【００３４】
【数１２】
2006234389.tif 000063
【００３５】
【数１３】
2006234389.tif 000064
本発明の光位相分布測定方法では、被測定光の位相分布を得ることを目的としているので、被測定光の複素振幅を同定する。つまり、光波検出センサーであるＣＣＤにより測定された光の強度を観測量とし、上記数１１、数１２及び数１３で表す観測方程式を満足する被測定光を表す複素数ベクトル（複素振幅）
2006234389.tif 000065
を同定する位相分布同定逆問題を解くことになる。
【００３６】
しかしながら、上記数１１、数１２及び数１３で表す観測方程式を満たす被測定光
2006234389.tif 000066
を求めるのは、これらの観測方程式が非線形方程式であるため、直接解を求めるのは困難である。
【００３７】
そこで、本発明では、被測定光
2006234389.tif 000067
を求める位相分布同定逆問題を非線形最適化問題として定式化する。定式化された非線形最適化問題において、設計変数を被測定光を表す複素数ベクトル
2006234389.tif 000068
とし、下記数１４で表す非線形の目的関数（評価関数）を設定する。
【００３８】
【数１４】
2006234389.tif 000069
つまり、被測定光
2006234389.tif 000070
の位相分布同定問題を、設計変数を
2006234389.tif 000071
とし、数１４で表す評価関数を有する非線形最適化問題に帰着させることができる。
【００３９】
ここで、本発明で用いる光位相分布測定系の具体的な光学系の一例を図２に示して説明する。図２に示すような光学系では、
2006234389.tif 000072
を直接観測していることになり、Ｈ１は単位マトリクスになり、[Ｈ]＝[Ｈ２][Ｉ]になる。また、
2006234389.tif 000073
の観測面を図２のようにレンズの焦点面に設置すれば、下記数１５が成立する。
【００４０】
【数１５】
2006234389.tif 000074
ただし、
2006234389.tif 000075
のフーリエ変換を表す。
＜４＞被測定光の位相分布の級数展開
次に、本発明の光位相分布測定方法では、被測定光
2006234389.tif 000076
の位相分布を級数展開することにより、設計変数を
2006234389.tif 000077
とし、数１４で表す評価関数を有する非線形最適化問題（以下、単に、第１の非線形最適化問題とも称する）を後述する第２の非線形最適化問題に変換する。
【００４１】
ここで、被測定光（複素振幅）
2006234389.tif 000078
のｉ番目成分をｇｉとし、
2006234389.tif 000079
をＣＣＤにより測定された被測定光
2006234389.tif 000080
のｉ番目成分の強度とし、φｉを被測定光
2006234389.tif 000081
のｉ番目成分に対応する位相角とすると、下記数１６が成り立つ。
【００４２】
【数１６】
2006234389.tif 000082
被測定光（複素振幅）
2006234389.tif 000083
と被測定光の位相との関係を分布全体で
2006234389.tif 000084
と表現する。
【００４３】
被測定光の位相分布φを下記数１７のようにゼルニケ多項式（Zernike circle
polynomial）に展開する。
【００４４】
【数１７】
2006234389.tif 000085
ここで、ｒ,θはスクリーン上の極座標成分を表す。Ｗｊはゼルニケ（Zernike）の基底関数で、ｚｊはゼルニケ（Zernike）の係数である。数１７において、
2006234389.tif 000086
次以降の高次の項を無視し、スクリーン面を空間的に離散化したときの被測定光の位相分布のベクトル表現は、下記数１８となる。
【００４５】
【数１８】
2006234389.tif 000087
ここで、[Ｗ]は、ゼルニケ（Zernike）の基底関数を離散化したベクトルを
2006234389.tif 000088
個横に並べたマトリクスである。また、ｚは、
2006234389.tif 000089
個のゼルニケ（Zernike）の係数からなるベクトルである。
【００４６】
本発明では、被測定光の位相分布を級数展開する際に
2006234389.tif 000090
次以降の高次の項を無視できるとの理由は、実際に産業的に取り扱う、レンズやＤＶＤピックアップなどの光学素子の評価においては、位相(収差)の分布は比較的滑らかであり、ゼルニケ（Zernike）の級数の３０〜７０項ほどで、十分に位相分布を表現できることが干渉計による測定結果を調査することにより裏付けられた。要するに、被測定光の位相分布をゼルニケ級数展開する場合には、３０〜７０個の未知数を同定するだけで、被測定光の位相分布を得ることができる。
【００４７】
本発明の光位相分布測定方法では、それぞれの画素（ここで言うそれぞれの画素とは、離散化された光波の複数の成分のことを意味する）の複素振幅の成分を直接同定するのではなく、被測定光の位相分布を級数展開し、例えば、数１８のように低次の級数で展開し、少ない数の展開係数を同定するようにしている。
【００４８】
すなわち、本発明では、設計変数を級数展開の展開係数（本例では、ゼルニケ係数）とし、数１９で表す評価関数を有する非線形最適化問題（以下、第２の非線形最適化問題とも称する）に帰着させ、この第２の非線形最適化問題を解くことで、被測定光の位相分布を同定するようにしている。これにより、第２の非線形最適化問題の解を精度良く安定に求めることができ、かつ計算量を大幅に低減することができる。
【００４９】
【数１９】
2006234389.tif 000091
ここで、例えば、本発明の光位相分布測定方法を用いないで、設計変数を
2006234389.tif 000092
とし、数１４で表す評価関数を有する非線形最適化問題を解くことにすると、それぞれの画素における複素振幅の実部と虚部をそれぞれ同定する必要がある。一般的には、光波検出センサーとしてのＣＣＤの総画素数は膨大であり、例えば、画像サイズが６４×６４の粗い画像の場合でも、同定すべき未知数は、画像の総画素数ｎの２倍、２ｎ＝６４×６４×２＝８１９２となる。このような場合、数１４で表す評価関数を有する非線形最小化問題は、悪条件となり、解が安定しなかったり、また計算量が膨大になる問題が生じてしまう。
【００５０】
なお、以上では、被測定光の位相分布をゼルニケ（Zernike）の級数で展開した実施例を説明したが、本発明では、それに限定されることがなく、他にもフーリエ（Fourier）級数展開、ラグランジュ級数展開、スプライン関数による展開、やベッセル関数による展開などが利用できることは言うまでもない。また、位相分布φを展開する他にも、例えば、複素振幅
2006234389.tif 000093
を直接に展開しても、同等な効果を得ることができる。
【００５１】
さらに、被測定光の位相分布がある設計された位相分布から近い分布と予想される場合には、設計された位相分布をφ０(ｘ,ｙ)として、これを定数項とする下記数２０のような級数に展開すると、級数で展開した項（下記数２０の右辺２項目）の値を小さくすることができ、目的関数の非線形性を緩和することができる。この結果、計算が効率的になり、最適化問題の悪条件性を緩和することもできる。
【００５２】
【数２０】
2006234389.tif 000094
この適用例として、レンズの収差を測定する応用を考えたときに、レンズの収差分布の設計値をφ０(ｘ,ｙ)とし、収差分布の製品ごとのばらつきの部分だけを級数で展開することがあげられる。また、光学部品の収差測定の応用においては、ロットごとで収差分布が同じ傾向になることが知られているが、この場合、計測済みの部品の統計量である収差分布の平均値をφ０(ｘ,ｙ)とすることで、同様の効果を得ることができる。
【００５３】
以下、本発明の光位相分布測定方法の有効性を確認するために、真の位相分布が既知の被測定光の位相分布を本発明で同定することにした。即ち、既知の位相分布と強度分布を持つ光が、レンズ焦点系および回折系の２種類の光学系を伝播し、出力される光の強度分布を予め計算で求め、これに誤差を１％を含ませ、模擬データとする。この模擬データから本発明の光位相分布測定方法により、逆に位相分布を同定し、あらかじめ設定していた既知の位相分布と本発明の同定結果を比較する。なお、入力光の振幅分布は一様分布とした。
【００５４】
ここで、被測定光の既知の位相分布（真の位相分布）を図３に、本発明の光位相分布測定方法によって同定した被測定光の位相分布を図４にそれぞれ示す。また、本発明の計算結果の平均誤差は０.０１２[rad]で、計算時間は２０９[sec]である。よって、本発明を用いれば、被測定光の位相分布を高速且つ高精度に推定できることが確認された。
＜５＞被測定光の伝播の数理モデル
次に、本発明で用いる光位相分布測定系の光学系に対し、その光学特性を表すシステムマトリクス、つまり、光伝播の伝達マトリクスについて説明する。
【００５５】
ところで、均一な媒体中の単色光の定常な伝播は、下記数２１で表すヘルムホルツ方程式を満足する。
【００５６】
【数２１】
2006234389.tif 000095
ここで、ｕは光の複素振幅とし、ｋは光の波長で決まる波数とする。図５に示されるように、仮想境界Γで分かられる半無限領域Ωを考える。数２１にヘルムホルツ方程式の３次元基本解をかけて、領域Ωで積分し部分積分を２回行うと、数２２で表すキルヒホッフの積分方程式が得られる（非特許文献８を参照）。ここで、ゾンマーフェルト（Sommerfeld）の放射条件を考慮して、無限遠の境界の積分は無視している。
【００５７】
【数２２】
2006234389.tif 000096
ただし、φ＊は、ヘルムホルツ方程式の３次元基本解で、下記数２３で表される。
【００５８】
【数２３】
2006234389.tif 000097
ここで、ｒは点Ｐと点Ｑの距離である。
【００５９】
また、ｑ＊およびｑは、下記数２４で定義され、それぞれφ＊およびｕの仮想境界Γに対する外向き法線ベクトルｎ方向の微分である。
【００６０】
【数２４】
2006234389.tif 000098
領域Ω内の任意の点Ｐの複素振幅ｕ(Ｐ)を仮想境界Γ上の値ｕ(Ｑ)およびｑ(Ｑ)から求めるために、数２２のキルヒホッフの積分方程式を利用することができる。
【００６１】
一般的に、仮想境界Γ上の複素振幅ｕと点Ｐの複素振幅ｕ(Ｐ)の関係を導出するために、仮想境界Γに入射する光を球面波などに近似し、さらにその分布の勾配を計算して、仮想境界Γ上のｑを求め、数２２に代入する方法を用いる（非特許文献８を参照）。
【００６２】
ところが、実際に被測定光が球面波でないときには、光伝播の数理モデルに球面波近似を用いることによって、誤差を生じてしまう問題がある。特に、被測定光の分布が波長と同じくらいのスケールの場合には、球面波近似による誤差の問題が更に顕著になってしまう。さらに、被測定光がレーザ光である場合には、球面からの位相のずれを対象とするので、これを無視する球面波近似は、本質的な矛盾となってしまう問題もある。
【００６３】
本発明では、このような問題を解決するために、仮想境界Γ上の複素振幅ｕを境界条件として、数２２のキルヒホッフの積分方程式を直接解くことにより、入射波（被測定光）に特別な仮定をおくことなく仮想境界Γ上のｑを導出し、仮想境界Γ上の複素振幅ｕと点Ｐの複素振幅ｕ(Ｐ)の関係を得るようにしている。
【００６４】
次に、図６の模式図を参照しながら、本発明で用いられる光伝播の伝達マトリクス[Ｈ]の求め方について、以下のように具体的に説明する。
【００６５】
まず、点Ｐを仮想境界Γに近づける操作を施すと、下記数２５が得られる。
【００６６】
【数２５】
2006234389.tif 000099
様々な仮想境界Γ上の点Ｐに対し、数２５の積分方程式を空間的に離散化すると、下記数２６のマトリクス方程式が得られる。
【００６７】
【数２６】
2006234389.tif 000100
ここで、ｕ,ｑの添え字０は、第０の仮想境界Γ０上の値を表す。この第０の仮想境界Γ０が平面である場合には、上記数２６が下記数２７となる。
【００６８】
【数２７】
2006234389.tif 000101
ただし、[Ｇ０]のｉ,ｊ成分は、下記数２８で定義される。
【００６９】
【数２８】
2006234389.tif 000102
ここで、Γｊは離散化したｊ番目の仮想境界とし、
2006234389.tif 000103
はＱｊに対応した内挿関数とする。上記数２７を解くと、数２９が得られ、
2006234389.tif 000104
が求められる。
【００７０】
【数２９】
2006234389.tif 000105
同様に、第１の仮想境界Γ１上の様々な点Ｐに対し、数２２のキルヒホッフの積分方程式を空間的に離散化すると、下記数３０のマトリクス方程式が得られる。
【００７１】
【数３０】
2006234389.tif 000106
ここで、ｕ,ｑの添え字１は、第１の仮想境界Γ１上の値を表す。数３０に数２９を代入すると、下記数３１が得られる。
【００７２】
【数３１】
2006234389.tif 000107
本発明では、以上のように、入射波（被測定光）に特別な仮定をおくことなく、第０の仮想境界Γ０上の光の複素振幅
2006234389.tif 000108
と第１の仮想境界Γ１上の光の複素振幅
2006234389.tif 000109
とを関係づけることができ、よって、この関係を利用して、高精度の観測方程式を構築することが可能になる。
【００７３】
要するに、本発明では、
2006234389.tif 000110
を観測画像（つまり、ＣＣＤにより測定された強度分布画像）とすれば、本発明で用いられる光伝播の伝達マトリクス、つまり、上記数１１における光学系のシステムマトリクス[Ｈ]は、下記数３２で表すことができる。
【００７４】
【数３２】
2006234389.tif 000111
上記のように、被測定光の光の分布が波長と同じくらいのスケールの場合であっても、また、レーザ光の収差（位相）分布を測定する場合であっても、本発明の光位相分布測定方法を用いれば、被測定光を球面波で近似しないので、球面波近似で生じる誤差は発生しないため、高精度の観測方程式を構築することが可能で、よって、被測定光の位相分布を高精度に同定することができる。
＜６＞光位相分布測定方法の高速化
本発明では、フーリエ変換及び逆フーリエ変換を利用することにより、被測定光の位相分布をより高速に同定することができる。
【００７５】
まず、本発明に係る光位相分布測定方法における基本的な計算手順を以下のように説明する。なお、本実施例では、被測定光の位相分布の級数展開にゼルニケ多項式を用いる。
ステップ１：
被測定光の位相分布の級数展開の係数の初期値を設定する。本例では、ゼルニケ（Zernike）の係数ｚの初期値を設定する。
ステップ２：
被測定光の位相分布を級数展開する。本例では、被測定光の位相分布φを下記数３３に従って計算（展開）する。
【００７６】
【数３３】
2006234389.tif 000112
ステップ３：
被測定光の複素振幅
2006234389.tif 000113
を下記数３４に従って計算する。
【００７７】
【数３４】
2006234389.tif 000114
ステップ４：
被測定光の第１の仮想境界における複素振幅
2006234389.tif 000115
を下記数３５に従って計算する。
【００７８】
【数３５】
2006234389.tif 000116
ステップ５：
ステップ４で算出された被測定光の第１の仮想境界における複素振幅から光強度を算出し、算出された光強度（算出値）と、ＣＣＤにより検出された光強度（観測値）と比較する。比較した結果が終了条件を満足させた場合（つまり、第２の非線形最適化問題の目的関数が極小点に達した場合）に終了し、一方、比較した結果が終了条件を満足しない場合（つまり、第２の非線形最適化問題の目的関数が極小点に達しない場合）に、ステップ６に進む。
ステップ６：
被測定光の位相分布の級数展開の係数を修正してから、ステップ２に戻る。本例では、ゼルニケ（Zernike）の係数ｚを修正してから、ステップ２に戻る。
【００７９】
つまり、本発明では、算出値と観測値との比較結果が終了条件（つまり、第２の非線形最適化問題の目的関数）を満足するまで、以上の手順（つまり、ステップ２〜ステップ６）を繰り返し行うようになっている。
【００８０】
修正には、例えば、非特許文献７に開示されている『準ニュートン法』を利用する場合に、各ステップで未知数に対する感度の計算も必要となる。各ステップでの感度の計算は、以下のようになる。
ステップ２での位相分布の感度の算出：
被測定光の位相分布の感度を下記数３６に従って計算する。
【００８１】
【数３６】
2006234389.tif 000117
ステップ３での被測定光の複素振幅の感度の算出：
被測定光の複素振幅の感度を下記数３７に従って計算する。
【００８２】
【数３７】
2006234389.tif 000118
ステップ４での被測定光の第１の仮想境界における複素振幅の感度の算出：
被測定光の第１の仮想境界における複素振幅の感度を下記数３８に従って計算する。
【００８３】
【数３８】
2006234389.tif 000119
ステップ５での算出された強度の感度の算出：
算出された強度（算出値）の感度を下記数３９に従って計算する。
【００８４】
【数３９】
2006234389.tif 000120
ここで、
2006234389.tif 000121
は、対角成分を
2006234389.tif 000122
とする対角行列である。
【００８５】
また、マトリクス[Ｈ]が画像（強度画像）の画素数の次元ｎを持つため、ステップ４での計算量及び記憶容量が膨大となる。つまり、ステップ４において、積の計算量はＯ(ｎ２)であり、また、記憶容量もＯ(ｎ２)となり、深刻な問題となってしまう。
【００８６】
本発明の光位相分布測定方法を測定機に組み込む際には、高速計算が要求されるが、これを記憶容量が限られているＤＳＰ（Digital Signal Processor）などのマイクロプロセッサで実現する際には問題となる。
【００８７】
数３２で表すマトリクス[Ｈ]のｉ,ｊ成分は、画素ｉと画素ｊとの距離ｒijの関数である。すなわち、ステップ４での計算は、つまり、下記数４０の計算は、「畳み込み積分」で表現することができる。
【００８８】
【数４０】
2006234389.tif 000123
上記数４０に２次元フーリエ変換を施すと、下記数４１となる。
【００８９】
【数４１】
2006234389.tif 000124
ここで、＊は各分布の２次元フーリエ変換を表す。
2006234389.tif 000125
の２次元フーリエ変換で、
2006234389.tif 000126
の２次元フーリエ変換である。
2006234389.tif 000127
は畳み込み積分の積分核の関数の２次元フーリエ変換である。
2006234389.tif 000128
はｎ次のベクトルとなる。
【００９０】
次に、図７を参照しながら、本発明において、ステップ４での畳み込み積分演算をせずに、フーリエ変換及び逆フーリエ変換を利用することにより、被測定光の光位相分布をより高速に同定する手順を具体的に説明する。
【００９１】
本発明では、被測定光の光位相分布をより高速に同定するために、図７に示すように、ステップ４での計算、つまり、上記数４０の計算を畳み込み積分演算を行わずに次の手順に沿って行う。
ステップＡ１：
2006234389.tif 000129
を２次元フーリエ変換すると、下記数４２が得られ、ちなみに、数４２の計算量はＯ(ｎ)である。
【００９２】
【数４２】
2006234389.tif 000130
ステップＡ２：
前述した数４１を計算する。ちなみに、数４１の計算量もＯ(ｎ)である。
ステップＡ３：
ステップＡ２で算出された
2006234389.tif 000131
を２次元フーリエ逆変換すると、下記数４３が得られ、ちなみに、数４３の計算量はＯ(ｎ log ｎ)である。
【００９３】
【数４３】
2006234389.tif 000132
以上のように、本発明では、ステップ４を行うことで
2006234389.tif 000133
を求める代わりに、ステップＡ１〜ステップＡ３を行って、つまり、被測定光の複素振幅をフーリエ変換し（ステップＡ１）、像空間で周波数成分ごとに独立に積演算を行い（ステップＡ２）、逆フーリエ変換により実空間に戻す（ステップＡ３）ことにより、
2006234389.tif 000134
を求めることもできるので、ステップ４で必要で計算コストが高い畳み込み積分をしなくても、
2006234389.tif 000135
を求めることができる。
【００９４】
ちなみに、以上の手順（ステップＡ１〜ステップＡ３）では、合計Ｏ(ｎ)＋Ｏ(ｎ)＋Ｏ(ｎ log ｎ)の計算量となるので、ステップ４の計算量Ｏ(ｎ２)に比べて、ｎが大きな画像問題では、このようにフーリエ変換及びフーリエ逆変換を利用することにより、計算量を大幅に削減することができ、よって被測定光の光位相分布の同定をより高速化させることができる。
【００９５】
なお、本発明では、フーリエ変換及びフーリエ逆変換の代わりに、高速フーリエ変換（ＦＦＴ）と高速フーリエ逆変換（ＩＦＦＴ）を適用すれば、つまり、ステップＡ１において、２次元高速フーリエ変換（２ＤＦＦＴ）を、ステップＡ３において、２次元高速フーリエ逆変換（２ＤＩＦＦＴ）を用いれば、更なる高速化が図れることは言うまでもない。
【図面の簡単な説明】
【００９６】
【図１】本発明で用いる光位相分布測定系の測定原理及び数理モデルを説明するための模式図である。
【図２】本発明で用いる光位相分布測定系の具体的な光学系の一例を説明するための模式図である。
【図３】被測定光の真の位相分布を示す図である。
【図４】本発明の光位相分布測定方法で同定した被測定光の位相分布を示す図である。
【図５】キルヒホッフの積分方程式を説明するための模式図である。
【図６】本発明で用いられる光伝播の伝達マトリクスを説明するための模式図である。
【図７】本発明において、高速フーリエ変換及び逆高速フーリエ変換を利用することにより、被測定光の位相分布をより高速化させる原理を説明するための模式図である。
【図１】
2006234389.tif 000153
【図２】
2006234389.tif 000154
【図３】
2006234389.tif 000155
【図４】
2006234389.tif 000156
【図５】
2006234389.tif 000157
【図６】
2006234389.tif 000158
【図７】
2006234389.tif 000159