(57)【要約】【課題】  特別な測定装置を用いることなく、簡便に測定可能な光強度分布の情報から光位相分布を同定することによって光位相分布を測定できるようにした光位相分布測定方法及び測定システムを提供する。【解決手段】  異なる複数の光学特性が既知の光学系と光波検出センサーとを備える光位相分布測定システムを用いて、被測定光の強度分布を測定する強度分布測定ステップと、強度分布測定ステップで得られた前記強度分布と前記各光学系の光学特性とに基づいて、前記光位相分布測定システムの観測方程式を設定する観測方程式設定ステップと、前記観測方程式から被測定光の位相分布を同定する位相分布同定逆問題を設定する位相分布同定逆問題設定ステップと、設定された前記位相分布同定逆問題を非線形の最適化問題として定式化する最適化問題定式化ステップと、前記最適化問題を解くことで、前記被測定光の位相分布を同定する位相分布同定ステップとを有する。【選択図】  図２
2005331440.tif 000002
【特許請求の範囲】
【請求項１】
異なる複数の光学特性が既知の光学系と光波検出センサーとを備え、
被測定光を前記各光学系にそれぞれ入力し、強度と位相を変調し、出力された前記被測定光を前記光波検出センサーにより検出し、検出された前記被測定光の強度分布を画像として測定することを特徴とする光位相分布測定システム。
【請求項２】
前記光波検出センサーはＣＣＤ撮像素子である請求項１に記載の光位相分布測定システム。
【請求項３】
前記複数の光学系は回折系とレンズ焦点系とから構成される請求項１に記載の光位相分布測定システム。
【請求項４】
被測定光の位相分布を測定するための光位相分布測定方法であって、
請求項１乃至請求項３の何れかの光位相分布測定システムを用いて、前記被測定光の強度分布を測定する強度分布測定ステップと、
強度分布測定ステップで得られた前記強度分布と前記各光学系の光学特性とに基づいて、前記光位相分布測定システムの観測方程式を設定する観測方程式設定ステップと、
前記観測方程式から被測定光の位相分布を同定する位相分布同定逆問題を設定する位相分布同定逆問題設定ステップと、
設定された前記位相分布同定逆問題を非線形の最適化問題として定式化する最適化問題定式化ステップと、
前記最適化問題を解くことで、前記被測定光の位相分布を同定する位相分布同定ステップと、
を有することを特徴とする光位相分布測定方法。
【請求項５】
前記最適化問題において、設計変数を被測定光を表す複素数ベクトルとし、目的関数に前記被測定光の位相分布に関する先験情報に基づいて設定される適切化関数が含まれる請求項４に記載の光位相分布測定方法。
【請求項６】
前記先験情報は、前記被測定光の位相は滑らかに分布しているであろうという情報である請求項５に記載の光位相分布測定方法。
【請求項７】
前記最適化問題の目的関数の最小化は準ニュートン法を用いて行う請求項４に記載の光位相分布測定方法。
【請求項８】
前記最適化問題の設計変数の隣り合う所定の点を代表するものを第２の設計変数とし、前記第２の設計変数で表現されたサイズの小さな第２の最適化問題の推定解を、前記最適化問題の前記設計変数の初期推定解として用い、前記最適化問題を解くようにする請求項４に記載の光位相分布測定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、光位相分布測定方法及び測定システムに関し、特に、簡単且つ経済的に光位相分布を測定できるようにした光位相分布測定方法及び測定システムに関する。
【背景技術】
【０００２】
光位相分布（以下、光の位相分布とも称する）には、光源の波面形状、透過物体の屈折率分布、反射物体の表面形状などの情報が含まれている。光位相分布の計測は、光学システムの設計評価、半導体レーザ素子及びレンズなどの光学素子の品質管理、非接触非破壊検査や顕微鏡などの表面形状計測、光波面の整形や天体観測に利用される波面補償光学など多くの分野で必要とされている。光位相分布の定量的な測定は、工学的に非常に重要である（非特許文献１、非特許文献２、非特許文献３を参照）。
【０００３】
従来より、ハードウェア（つまり、特別な測定装置）を用いた、光位相分布の測定方法として、例えば、干渉計により被測定光と参照光との相対的な位相分布を求める方法（以下、干渉計による光位相分布測定方法）（非特許文献４及び非特許文献５を参照）や、シャック・ハートマンセンサーにより位相勾配分布を検出し位相分布を求める方法（以下、シャック・ハートマンセンサーによる光位相分布測定方法）（非特許文献６及び非特許文献７を参照）などがある。
【０００４】
上記の従来の光位相分布測定方法を次のようにより詳細に説明する。
【０００５】
まず、従来の干渉計による光位相分布測定方法とは、被測定光を２つ以上の光に分割し、別々の光路を通ったあと再び重ね合わせ、光路差により発生する干渉縞を捉え、これを解析して光の位相分布を求めるものである。
【０００６】
つまり、図５３に示すように、従来の干渉計による光位相分布測定方法では、被測定光を２つに分割し、一方の波面（位相分布）に対して他方の波面がどれだけ異なっているかを調べる。基準となる光は参照光と呼ばれる。参照光はリファレンスレンズとスペイシャルフィルタによって均一な波面の光として取りだされ、被測定光と干渉する。参照光と被測定光の光路差が波長の整数倍であれば干渉した光は明るくなり、光路差が波長の整数倍から波長の半分だけずれていれば干渉光は暗くなる。場所により光路差が一定でない場合には、明暗の干渉縞が観察される。
【０００７】
干渉縞画像から位相分布を求めるには、例えば、図５４に示されるような光路長の異なる複数の干渉縞画像を解析する必要がある。干渉により出来た干渉縞の正弦波状の強度分布を走査し、強度の変化から位相を求める、非特許文献４に開示されているフリンジスキャン法が代表的な解析手法である。
【０００８】
また、従来のシャック・ハートマンセンサーによる光位相分布測定方法では、シャック・ハートマンセンサー（ＳＨＷＳ:Shack-Hartmann  wavefront  sensor）が多数の小型レンズとＣＣＤとから構成される。図５５に示されるように、小型レンズに入射した光は、その部分の局所的な光の位相勾配によって、焦点を結ぶ位置が焦点面上でΔｓずれる。このずれを全ての小型レンズにおいて検出する。図５６に示されたような焦点面上の結像パターンは、ハートマンパターン（Hartmann  pattern）と呼ばれ、ＣＣＤで検出される。
【０００９】
非特許文献８に開示されるように、Δｓと位相の勾配は比例関係にあり、このことから局所的な位相の勾配情報を知ることができ、総合することで全体の位相分布を求めることができる。
【非特許文献１】ケイ．エイ．ヌージェント（K.A.Nugent），デイ．パガニーニ（D.paganin），ティー．イー．グリエーフ（T.E.Gureyev）共著,「フェイズ  オデッセイ（A  Phase  Odyssey）」、フィジクス  トゥデイ（Physics  Today），第54巻，第8号,2001年,ｐ27−32
【非特許文献２】エル．エイ．トンプスン（L.A.Thompson）著,「アダプティブ  オプティクス  イン  アストロノミー（Adaptive  Optics  in  Astronomy）」，フィジクス  トゥデイ（Physics  Today），第47巻，第12号,1994年,ｐ24−31
【非特許文献３】三好隆志  著,「光応用ナノインプロセス計測技術の最新動向」,機械の研究,第54巻，第9号、2002年,ｐ928−933
【非特許文献４】「光測定器ガイド」,オプトロニクス社,1998年
【非特許文献５】高田元弘  著,「干渉計によるレーザ光の波面収差測定」,光アライアンス，第2版、2001年,ｐ19−22
【非特許文献６】ダブリュー．エイチ．サウスウェル（W.H.Southwell）著，「ウェイブ  フロント  エスティメイション  フロム  ウェイブフロント  スロープ  メジャーメント（Wave  front  estimation  from  wavefront  slope  measurements）」,ジェイ.  オーピーティー.  エスオーシー.  エイエム.（J.Opt.Soc.Am.）,第70巻,第8号,1980年,ｐ998−1005
【非特許文献７】アール．アーウォン（R.Irwan），アール．ジー．レイン（R.G.Lane）共著,「アナリシス  オブ  オプティマル  セントロイド  エスティメイション  アプライド  トゥー  シャック  ハートマン  センシング（Analysis  of  optimal  centroid  estimation  applied  to  Shack-Hartmann  sensing）」,アプライド  オプティクス（Applied  Optics）,第38巻,第32号,1999年,ｐ6737−6743
【非特許文献８】イリット  ラス（ILIT  RAS）,「シャック−ハートマン  ウェイブフロント  センサー  フォア  レーザ  ビーム  アナリシス（Shack-Hartmann  Wavefront  Sensor  for  Laser  Beam  Analysis）」、インターネット＜http://www.laser.ru/adopt/Science/sens/sensor.htm＞
【非特許文献９】谷田貝豊彦  著,「光とフーリエ変換」,朝倉書店,1992年
【非特許文献１０】吉村武晃  著,「光情報工学の基礎」,コロナ社,2000年
【非特許文献１１】矢部・八卷  共著,「非線形計画法」,朝倉書店,1999年
【発明の開示】
【発明が解決しようとする課題】
【００１０】
しかしながら、上述した従来の光位相分布測定方法に使用されている特別な測定装置（つまり、干渉計やシャック・ハートマンセンサーなどのハードウェア）は、精密な光学機器から構成されているので非常に高価で、また、こういった特別な測定装置を設置するスペースも必要となるという問題があった。
【００１１】
本発明は、上述のような事情よりなされたものであり、本発明の目的は、特別な測定装置を用いることなく、簡便に測定可能な光強度分布の情報から光位相分布を同定することによって光位相分布を測定できるようにした光位相分布測定方法及び測定システムを提供することにある。
【課題を解決するための手段】
【００１２】
本発明は、光位相分布測定システムに関し、本発明の上記目的は、異なる複数の光学特性が既知の光学系と光波検出センサーとを備え、被測定光を前記各光学系にそれぞれ入力し、強度と位相を変調し、出力された前記被測定光を前記光波検出センサーにより検出し、検出された前記被測定光の強度分布を画像として測定することにより、或いは、前記光波検出センサーはＣＣＤ撮像素子であるようにすることにより、或いは、前記複数の光学系は回折系とレンズ焦点系とから構成されることによって効果的に達成される。
【００１３】
また、本発明は光位相分布測定方法に関し、本発明の上記目的は、被測定光の位相分布を測定するための光位相分布測定方法であって、本発明の光位相分布測定システムを用いて、前記被測定光の強度分布を測定する強度分布測定ステップと、強度分布測定ステップで得られた前記強度分布と前記各光学系の光学特性とに基づいて、前記光位相分布測定システムの観測方程式を設定する観測方程式設定ステップと、前記観測方程式から被測定光の位相分布を同定する位相分布同定逆問題を設定する位相分布同定逆問題設定ステップと、設定された前記位相分布同定逆問題を非線形の最適化問題として定式化する最適化問題定式化ステップと、前記最適化問題を解くことで、前記被測定光の位相分布を同定する位相分布同定ステップとを有することにより、或いは、前記最適化問題において、設計変数を被測定光を表す複素数ベクトルとし、目的関数に前記被測定光の位相分布に関する先験情報に基づいて設定される適切化関数が含まれるようにすることにより、或いは、前記先験情報は、前記被測定光の位相は滑らかに分布しているであろうという情報であることにより、或いは、前記最適化問題の目的関数の最小化は準ニュートン法を用いて行うことにより、或いは、前記最適化問題の設計変数の隣り合う所定の点を代表するものを第２の設計変数とし、前記第２の設計変数で表現されたサイズの小さな第２の最適化問題の推定解を、前記最適化問題の前記設計変数の初期推定解として用い、前記最適化問題を解くようにすることによって効果的に達成される。
【発明の効果】
【００１４】
本発明に係る光位相分布測定方法及び光位相分布測定システムによれば、従来の光位相分布測定方法に用いられる特別な測定装置のような高価な光学機器を用いることなく、簡便かつ経済的に測定可能な光強度分布の情報から光位相分布を安定的に推定することができ、簡便かつ経済的な光位相分布の測定を可能にしたといった優れた効果を奏する。
【発明を実施するための最良の形態】
【００１５】
以下、本発明を実施するための最良の形態を図面を参照して説明する。
＜１＞数理モデル
まず、本発明において使用される、光波や光位相分布測定系の数理モデルについて説明する。
＜１−１＞光波の数理モデル
本発明において、位相分布測定の対象とする光（つまり、被測定光）が、例えばレーザ光などに代表されるような、単一の波長成分を持つ単色光で、且つ、異なる時刻における位相が常に一定の相関関係を持つ時間的にコヒーレントな光であることを前提とする。
＜１−１−１＞光波の複素数表示
光は電磁波の一種であり、電場、磁場が振動しながら伝搬する波である。電場Ｅ、磁場Ｈ、誘電率ε、透磁率μを関係づけるマックスウェルの方程式は，下記数１、数２、数３及び数４といった４つの式から構成される。
【００１６】
【数１】
2005331440.tif 000003
【００１７】
【数２】
2005331440.tif 000004
【００１８】
【数３】
2005331440.tif 000005
【００１９】
【数４】
2005331440.tif 000006
これらの式より、以下のような波動方程式を得る。
【００２０】
【数５】
2005331440.tif 000007
【００２１】
【数６】
2005331440.tif 000008
電界Ｅに関する波動方程式である上記数５の解の一つは、下記数７のように表される。
【００２２】
【数７】
2005331440.tif 000009
ここで、ａは振幅で、ｅｐは電界の振動方向を表す単位偏光ベクトルで、ωは角周波数で、ｔは時刻で、ｋは波数ベクトルで、ｒは位置ベクトルで、φは初期位相である。また、ｉは虚数を表す。なお、磁場Ｈについても同様であるので以下省略する。
【００２３】
ここで、位相の時間に依存しない成分をまとめて、下記数８と置きなおし、時間振動項を分離して表現すると、下記数９になる。
【００２４】
【数８】
2005331440.tif 000010
【００２５】
【数９】
2005331440.tif 000011
本発明では、光学における計算において、このように光波の時間振動項を分離し、複素数ｖで光波の強度ａと位相θを表現する。複素数ｖは複素振幅と呼ばれ、下記数１０のようになる。
【００２６】
【数１０】
2005331440.tif 000012
また、複素振幅ｖを実数部ｖｒと虚数部ｖｉで分離して表現した場合に、下記数１１、数１２、数１３のようになる。
【００２７】
【数１１】
2005331440.tif 000013
【００２８】
【数１２】
2005331440.tif 000014
【００２９】
【数１３】
2005331440.tif 000015
本発明では、上述したように、光を複素数で表現するものとすることを前提とする。
＜１−１−２＞光波の離散表示
後述するように、本発明では、光位相分布測定系において、被測定光を検出するための光波検出センサーとして、ＣＣＤ撮像素子（以下、単にＣＣＤとも呼ばれる）を用いることが好ましい。光波検出センサー（本実施形態では、ＣＣＤ）により、検出されて画像として測定されるのは、被測定光の強度、すなわち、上述した複素振幅の絶対値である。測定された強度の空間分布は、二次元の離散配列で与えられるが、本発明では、光波を表す複素振幅の分布を取扱う場合に、図１に示されるような複素数の二次元配列を離散的に一次元に並べ換え、下記数１４で表される複素数ベクトルｖで表現する。
【００３０】
【数１４】
2005331440.tif 000016
ここで、光位相分布測定系で測定されるのは、複素数ベクトルのそれぞれの成分の絶対値である。以降、複素数ベクトルｖのそれぞれ成分の絶対値を並べた実数ベクトルを下記数１５のように表記するものとする。
【００３１】
【数１５】
2005331440.tif 000017
＜１−２＞光位相分布測定系（以下、単に測定系とも称する）の数理モデル
＜１−２−１＞光位相分布測定系
光位相分布測定系は、光学特性が既知である複数（少なくとも２つ以上）の異なる光学系と光波検出センサーとから構成される。光位相分布測定系の光波検出センサーとしては、ＣＣＤ撮像素子を用いるのが好ましいが、それに限定されることなく、他のセンサーを用いても良い。なお、本実施形態において、光波検出センサーとしては、ＣＣＤを用いる。
【００３２】
本発明では、このような光位相分布測定系を用いて、光波検出センサー（ＣＣＤ）により被測定光の強度分布画像を測定し、そして、得られた被測定光の強度分布画像や光学系の光学特性などの情報に基づいて、光位相分布測定方法を利用して、被測定光の位相分布を同定することによって、被測定光の位相分布を測定できるようにしている。
【００３３】
前述したように、本発明では、単一の波長成分を持つ単色光を光位相分布の測定対象（つまり、被測定光）としているが、例えば、被測定光に注目する波長成分以外の波長をもつ光波が含まれてしまう場合には、光波検出センサー（ＣＣＤ）の前にフィルターを置き、注目する波長成分のみを取りだすことによって、本発明を適用することができる。
【００３４】
以下、図２を参照しながら、光位相分布測定系の測定原理及び数理モデルについて説明する。図２に示されるように、光位相分布測定系は、光学特性が既知である異なるｋ個の光学系と光波検出センサー（ＣＣＤ）とから構成されている。
【００３５】
光位相分布測定系では、まず、被測定光、すなわち、位相分布を測定したい光を光学特性が既知の光学系に入力し、強度と位相を変調する。そして、出力された光波を光波検出センサー（ＣＣＤ）により検出し、検出光の強度分布を画像として測定する。これらの手順を同じ被測定光に対し、複数の光学系（ここではｋ個の光学系、つまり、光学系１、光学系２、…、光学系ｋ）を用いて行う。
【００３６】
光位相分布測定系の測定原理を数学的に表現すると、被測定光を表す複素数ベクトルｇが光学系による変換マトリックスＨにより変換され、検出光を表す複素数ベクトルｆとなることを表している。すなわち、下記数１６と表すことができる。
【００３７】
【数１６】
2005331440.tif 000018
ここで、光波検出センサー（ＣＣＤ）により検出されるのは、検出光を表す複素数ベクトルｆのそれぞれの成分の絶対値を並べた実数ベクトルである。
【００３８】
本発明の光位相分布測定方法では、検出光の絶対値ベクトル｜ｆ｜を観測量とし、被測定光を表す複素数ベクトルｇを同定することによって、被測定光の位相分布の測定を可能にしている。
【００３９】
本発明の光位相分布測定系では、使用する光学系の種類および数は任意である。適切な光学系を適用することで、光位相分布測定系の観測方程式の悪条件性を克復し、唯一解を求めることが可能となる。
＜１−２−２＞光学系による光波の変換の数理モデル
本発明の光位相分布測定系に適用できる光学系は、光学特性が既知であれば、どのようなものでも構わない。簡単に利用できるものとして、例えば、レンズ、レンズの焦点はずれ、瞳関数の異なるレンズ、回折格子、空気中の回折（即ち、光学系を介さず直接入射させること）などの光学系が挙げられる。
【００４０】
以下では、三つの光学系の例をあげ、光学系による光波の変換の数理モデルについて説明する。詳細については、非特許文献９、非特許文献１０などを参照されたい。
＜Ａ＞回折系による光波の変換の数理モデル
図３に示すように、ｘｙ面内に光波の複素振幅分布ｇ(ｘ,ｙ)が局在しており、ｚ軸方向に距離Ｒだけ伝搬した後のＸＹ面内の複素振幅分布ｆ(Ｘ,Ｙ)とする。このとき、距離Ｒが波長に比べ十分大きく、回折波がｚ軸近傍に存在する場合を考えると、下記数１７のような光波の伝搬前後の関係を表すフレネル回折の式を得る。
【００４１】
【数１７】
2005331440.tif 000019
ここで、ｋは波数であり、波長をλとすると、ｋ＝２π／λとなる。下記数１８とすると、下記数１９で表すたたみ込み積分の定義より、上記数１７は、下記数２０のようにたたみ込み積分で表すことができる。
【００４２】
【数１８】
2005331440.tif 000020
【００４３】
【数１９】
2005331440.tif 000021
【００４４】
【数２０】
2005331440.tif 000022
上記数２０を離散化し、行列、ベクトル表示で書くと、下記数２１のようになる。
【００４５】
【数２１】
2005331440.tif 000023
ここで、簡単のために、１次元で考え、具体的に書き出してみる。上記数２１において、下記数２２、数２３とすると、フレネル回折の変換マトリックスは、下記数２４のようになる。
【００４６】
【数２２】
2005331440.tif 000024
【００４７】
【数２３】
2005331440.tif 000025
【００４８】
【数２４】
2005331440.tif 000026
ここで、下記数２５が成立する。
【００４９】
【数２５】
2005331440.tif 000027
ｈ(Ｘ,Ｙ)は、点光源入力に対する光学システムのインパルス応答を表し、点像応答関数（ＰＳＦ:point  spread  function）と呼ばれる。ＰＳＦは入力関数に依存せずシステム固有の性質を示す。光波の回折現象も上記数１８のＰＳＦを持つ光学システムと見なせる。この他の光学システムにおいても、ＰＳＦが前もって分かれば、入力ｇに対するシステムの出力ｆを求めることができる。
＜Ｂ＞結像レンズ系による光波の変換の数理モデル
図４に示すように、ｘｙ面内に光波の複素振幅分布ｇ(ｘ,ｙ)が局在しており、距離ｌだけ伝搬後、レンズを透過し、再び距離ｒだけ伝搬した後のＸＹ面内の複素振幅分布をｆ(Ｘ,Ｙ)とする。
【００５０】
光波ｇはレンズの直前まで、距離ｌだけ伝搬しフレネル回折する。レンズ直前の光波ｕ−(ｘ,ｙ)は、下記数２６のようになる。
【００５１】
【数２６】
2005331440.tif 000028
この光波が焦点距離ｆのレンズによる位相変調および瞳関数ｐ(ｘ,ｙ)による振幅変調を受けると、レンズ直後の光波ｕ＋(ｘ,ｙ)は、下記数２７のようになる。
【００５２】
【数２７】
2005331440.tif 000029
ここで、下記数２８で表す瞳関数は、レンズの開口を表す関数である。
【００５３】
【数２８】
2005331440.tif 000030
この光波が観測面まで、距離ｒだけ伝搬し、再びフレネル回折する。観測される出力光波ｆ(ｘ,ｙ)は、下記数２９のようになる。
【００５４】
【数２９】
2005331440.tif 000031
以上の三つの関係から、観測面での光波を求めることができる。観測位置は結像条件を満足し、距離ｒ,ｌおよび焦点距離ｆについてレンズの公式、結像倍率Ｍ、瞳関数のフーリエ変換形をそれぞれ下記数３０、数３１、数３２とする。
【００５５】
【数３０】
2005331440.tif 000032
【００５６】
【数３１】
2005331440.tif 000033
【００５７】
【数３２】
2005331440.tif 000034
これらの関係式を用いて整理し、定数項を無視すれば、入力光と出力光の関係は、下記数３３のようなたたみ込み積分で書ける。詳細は非特許文献１０を参照する。
【００５８】
【数３３】
2005331440.tif 000035
ここで、レンズの直径をＤとし、１次元で扱うと、関数Ｐは矩形関数のフーリエ変換で表される。すなわち、図５に示されるような結像光学系の点像応答関数ＰＳＦは、下記数３４で表すことが分かる。
【００５９】
【数３４】
2005331440.tif 000036
具体的な変換マトリックスの形は、下記数３５とすれば、数２４と同様に結像光学系の変換マトリックスを書くことができる。
【００６０】
【数３５】
2005331440.tif 000037
上記数３５式は、Ｘ＝ｘ近傍で値を持つので、結像光学系の変換行列は、帯行列になることがわかる。
【００６１】
レンズの直径Ｄ→∞のときは、ＰＳＦはｈ(ｘ)＝δ(ｘ)となり、点像に広がりのない、ボケのない完全な画像が観察面に出力される。つまり、変換マトリックスは単位行列となる。
＜Ｃ＞回折格子による光波の変換の数理モデル
図６に示すように、ｘｙ面内に光波の複素振幅分布ｇ(ｘ,ｙ)が局在しており、距離ｌだけ伝搬後、回折格子を透過し、再び距離ｒだけ伝搬した後のＸＹ面内の複素振幅分布をｆ(Ｘ,Ｙ)とする。
【００６２】
入力となる光波ｇは回折格子の直前まで、距離ｌだけ伝搬しフレネル回折する。回折格子の直前の光波ｕ−(ｘ,ｙ)は、下記数３６のようになる。
【００６３】
【数３６】
2005331440.tif 000038
次に、回折格子の複素透過率をｔ(ｘ,ｙ)とすると、回折格子を透過した直後の光波ｕ＋(ｘ,ｙ)は、下記数３７のようになる。
【００６４】
【数３７】
2005331440.tif 000039
この光波が観測面まで、距離ｒだけ伝搬し、再びフレネル回折する。観測される出力光波ｆ(ｘ,ｙ)は、下記数３８のようになる。
【００６５】
【数３８】
2005331440.tif 000040
以上の処理をまとめると、下記数３９のようになる。
【００６６】
【数３９】
2005331440.tif 000041
上記数３９を離散化し、行列ベクトル表示で書くと、下記数４０、数４１のようになる。
【００６７】
【数４０】
2005331440.tif 000042
【００６８】
【数４１】
2005331440.tif 000043
ここで、１次元で考え、具体的に書き出してみる。行列[ｈｒ],[ｈｌ]は、数２４と同様である。また、行列[ｔ]は、対角行列となり、下記数４２のようになる。
【００６９】
【数４２】
2005331440.tif 000044
＜２＞光位相分布測定方法（以下、位相分布同定手法とも呼ばれる）
以下では、位相分布同定逆問題を設定し、具体的な位相分布同定手法について説明する。
＜２−１＞位相分布同定逆問題の設定
＜２−１−１＞光位相分布測定系の観測方程式
入力光（つまり、被測定光）ｇをｎ個に離散化し、ｇ∈Ｃｎ、異なるｋ個の光学系を通過後の光をそれぞれｍ個に離散化しｆ１,ｆ２,…,ｆｋ∈Ｃｍ、これらの光学系による変換マトリックスをそれぞれ[Ｈ１],[Ｈ２],…,[Ｈｋ]∈Ｃｍ×ｎとする。また、光波検出センサーであるＣＣＤにより測定された強度分布を

2005331440.tif 000045
とすると、光位相分布測定系の観測方程式は、下記数４３、数４４のように書くことができる。
【００７０】
【数４３】
2005331440.tif 000046
【００７１】
【数４４】
2005331440.tif 000047
上記数４３は、光学系による被測定光ｇから検出光ｆへの変換式であり、上記数４４は、複素数ベクトルｆの絶対値ベクトル|ｆ|が測定された強度
2005331440.tif 000048
として与えられる、つまり、測定量であることを表す。ここで、既知量は、光学系の光学特性を表す複素数マトリックス[Ｈ]、および測定強度を表す絶対値ベクトル
2005331440.tif 000049
であり、未知量は、被測定光を表す複素数ベクトルｇと検出光を表す複素数ベクトルｆである。
【００７２】
本発明の光位相分布測定方法では、被測定光の位相分布を得ることを目的としているので、被測定光の複素振幅を同定する。つまり、複数の測定強度
2005331440.tif 000050
を観測量とし、光位相分布測定系の観測方程式を満足する被測定光の複素数ベクトルｇを同定する位相分布同定逆問題を解くこととなる。
＜２−１−２＞最適化問題
上記数４３と数４４で表す、光位相分布測定系の観測方程式を満たす入力光（つまり、被測定光）ｇを求めるのが、本発明の光位相分布測定方法の目的となるが、これらの方程式系は、非線形方程式のため、直接解を求めるのは困難である。
【００７３】
そこで、本発明の光位相分布測定方法では、入力光（つまり、被測定光）ｇを求める位相分布同定逆問題を非線形の最適化問題（以下、単に、最適化問題とも称する）として定式化する。定式化された最適化問題において、設計変数を被測定光を表す複素数ベクトルｇとし、下記数４５のような非線形の目的関数（つまり、評価関数）を設定する。
【００７４】
【数４５】
2005331440.tif 000051
この位相分布同定逆問題では、上記数４５を０とする、つまり、最小とする複素数ベクトルｇが推定解となる。
【００７５】
しかし、光位相分布測定系の測定値の誤差により、光位相分布測定系の観測方程式は、厳密には成立しない。このことが位相分布同定逆問題を悪条件化させることが考えられる。そこで、位相分布同定逆問題の悪条件性を克復するため、被測定光に関する先験情報を利用し、位相分布同定逆問題を適切化する。すなわち、上記数４５に適切化項をつけて得られた下記数４６を目的関数として用いる。
【００７６】
【数４６】
2005331440.tif 000052
ここで、ｗは適切化の重み係数で、Λ(ｇ)は適切化のために新たに加えられた関数である。
＜２−２＞最適化の手法
この最適化問題の数値解法および、具体的な先験情報の利用法については、以下のように説明する。
＜２−２−１＞準ニュートン法
本実施形態では、光位相分布測定方法において、上記数４６で表す目的関数の最小化をするにあたり、最適化問題に対する数値解法の中で、非特許文献１１に開示されている『準ニュートン法』を用いる。
【００７７】
まず、準ニュートン法について、以下のように簡単な説明をする。
【００７８】
ここで、ｎ変数の目的関数ｆ(ｘ)を最小にするｘ∈Ｒｎを見つけるという、最適化問題があるとする。
【００７９】
ニュートン法では、目的関数を２次モデルで近似し、そのモデル関数を局所的に最小化することを繰り返し行う。ヘッセ行列∇２ｆ(ｘ)、勾配ベクトル∇ｆ(ｘ)、探索方向ｄとすると、モデル関数は、下記数４７で表される。
【００８０】
【数４７】
2005331440.tif 000053
ヘッセ行列が正定値対称行列ならば、ｑ(ｄ)を最小にする点は、∇ｑ(ｄ)＝０を満たすベクトルｄを求めればよい。つまり、下記数４８で表す連立１次方程式を解けばよいことになる。
【００８１】
【数４８】
2005331440.tif 000054
上記数４８で表すこの方程式をニュートン方程式といい、解ｄをニュートン方向という。
【００８２】
しかし、いつでもヘッセ行列が正定値である保証はなく、ニュートン方向がｆ(ｘ)の降下方向になるとは限らない。そこで、ヘッセ行列を適当な正定値対称行列Ｂｋで近似し、探索方向を決定しようというのが、準ニュートン法である。
【００８３】
準ニュートン法には超一次の速い収束性を持ち、２階微分の計算を必要としないという利点がある。
【００８４】
次に、準ニュートン法のアルゴリズムを以下のように示す。
ステップ１  初期設定
近似解の初期点ｘ０、正定値対称な初期行列Ｂ０を与え、ｋ＝０とおく。
ステップ２  探索方向ｄｋの決定
下記数４９で表す連立一次方程式を解き、探索方向ｄｋを求める。
【００８５】
【数４９】
2005331440.tif 000055
ステップ３  直線探索
ｄｋ方向でのステップ幅αｋを決定する
ステップ４  近似解の更新
下記数５０のように、近似解を更新する。
【００８６】
【数５０】
2005331440.tif 000056
ステップ５  収束判定
収束条件が満たされていれば、ｘｋ＋１を解とみなし、終了する。さもなくば、ステップ６へ行く。
ステップ６  Ｂｋの更新
更新公式によりＢｋ＋１を計算する。
ステップ７  ｋ＋１をｋに代入して（つまり、ｋ：＝ｋ＋１とし）、ステップ２へ
また、直線探索法には様々な手法が存在するが、実用上には次のようなArmijoの基準を用いた簡単な探索法で十分である。
【００８７】
０＜ξ＜１であるような定数ξに対して、下記数５１を満たすαを選ぶ。
【００８８】
【数５１】
2005331440.tif 000057
Armijoの基準の幾何学的意味は図７で示される。つまり、図７は目的関数ｆ(ｘ)を探索方向ｄｋで切った断面図であり、横軸がステップ幅αに対応する。原点におけるｆ(ｘｋ＋αｄｋ)の接線ｙ＝ｆ(ｘｋ)＋α∇ｆ(ｘｋ)Ｔｄｋの傾きを緩和して得られた直線ｙ＝ｆ(ｘｋ)＋ξα∇ｆ(ｘｋ)Ｔｄｋよりも関数値が低くなる区間からステップ幅αを選ぶことになる。
【００８９】
Armijoの基準の直線探索アルゴリズムは、次のようになる。
ステップＡ  初期設定
現在の近似解ｘｋ、パラメータ０＜ξ＜１、０＜τ＜１を与える。
ステップＢ  探索
探索方向ｄｋで、Armijoの基準を満たすステップ幅αｋを求める。
【００９０】
ステップＢ１  βｋ,０＝１,ｉ＝０とおく。
【００９１】
ステップＢ２  下記数５２で表すArmijoの基準を満たすなら、ステップＣへ、さもなくばステップＢ３へいく。
【００９２】
【数５２】
2005331440.tif 000058
ステップＢ３  βｋ,ｉ＋１＝τβｋ,ｉ、ｉ：＝ｉ＋１とおいて（つまり、ｉ＋１をｉに代入して）、ステップＢ２へいく。
ステップＣ  ステップ幅決定
αｋ＝βｋ,ｉとおく。
また、Ｂｋの更新には種々の公式が考えられているが、代表的なものは、下記数５３で表すＢ公式のＢＦＧＳ（Broyden  Fletcher  Goldfarb  Shanno）更新公式である。
【００９３】
【数５３】
2005331440.tif 000059
なお、Ｂ０は単位行列とするのが一般的である。すなわち、初回の探索方向は最急降下方向である。また、ｓｋ,ｙｋは、下記数５４で表されるベクトルである。
【００９４】
【数５４】
2005331440.tif 000060
ヘッセ行列の逆行列を直接行列Ｈｋで近似することも考えられている。このとき、探索方向は連立方程式を解くこと無く、下記数５５として求めることができる。
【００９５】
【数５５】
2005331440.tif 000061
Ｈｋの更新公式は、上記数５３の逆行列を計算することによって得られる。具体的には、Ｈ公式のＢＦＧＳ更新公式は、下記数５６のような形になる。
【００９６】
【数５６】
2005331440.tif 000062
Ｂ公式およびＨ公式を用い、実際にそれぞれ最適化の計算を行ったところ、収束性に違いはほとんど見られなかった。Ｈ公式は連立方程式を解く必要がないため、計算コスト上有利である。本実施形態では、最適化問題ではＨ公式を採用した。
＜２−２−２＞被測定光に関する先験情報の利用
被測定光の光波の強度および位相の空間分布は、滑らかに分布していることが予測された場合に、つまり、「複素振幅（つまり、位相）は滑らかに分布しているであろう」ということが先験情報として利用できる。この先験情報を具体的に表現すると、『ある点における複素振幅とその隣り合う点における複素振幅との差は、小さいはずである』という情報である。これは最適化問題の目的関数に隣り合う点の差を最小化するという目的を加えることで実現できる。また、被測定光の複素ベクトルの分布が滑らかであるという先験情報は、最適化問題の目的関数に被測定光の複素ベクトルｇのノルムを最小化する目的を加えても実現できる。
【００９７】
一次元に分布する被測定光ｇをｎ個に離散化した場合、つまり、設計変数がｇ∈Ｃｎの場合に、隣り合う点の差の最小化を表す目的関数は、下記数５７のように表現できる。
【００９８】
【数５７】
2005331440.tif 000063
上記数５７のような目的関数を先験情報として、数４６の適切化関数とする。つまり、光波が一次元で分布している、設計変数がｇ∈Ｃｎの場合では、下記数５８のような適切化関数を用いる。
【００９９】
【数５８】
2005331440.tif 000064
ここで、マトリックスＣ∈Ｒｎ×ｎは、下記数５９で表される。
【０１００】
【数５９】
2005331440.tif 000065
同様に、被測定光ｇが二次元（ｎ×ｎ）で分布している場合には、隣り合う点の差の最小化を表す目的関数は、次のように表現できる。まず、簡単のため、設計変数を二次元配列[ｇ]∈Ｃｎ×ｎのまま表記すると、下記数６０のようになる。
【０１０１】
【数６０】
2005331440.tif 000066
上記数６０を設計変数がベクトル
2005331440.tif 000067
の表記で書き直すと、下記数６１のようになる。
【０１０２】
【数６１】
2005331440.tif 000068
ここで、上記数６１において、第一項は図８(Ａ)で表される行方向の、第二項は図８(Ｂ)で表される列方向の隣り合う点の差を表している。図８はｎ＝４の場合を示している。よって、上記数５８と同様に表せば、被測定光が二次元（ｎ×ｎ）で分布している場合の適切化関数は、下記数６２のようになる。
【０１０３】
【数６２】
2005331440.tif 000069
Ｃ１は、図８で表されるように、二次元の複素振幅分布の行方向の差をとるマトリックスであり、Ｃ２は列方向の差をとるマトリックスである。つまり、マトリックス
2005331440.tif 000070
は、数５９で表すマトリックスＣを用い、下記数６３で表される。
【０１０４】
【数６３】
2005331440.tif 000071
また、マトリックス
2005331440.tif 000072
は、１≦ｉ≦ｎ(ｎ−１)としたときに、（ｉ,ｉ）成分および（ｉ,ｉ＋ｎ）成分は、下記数６４のようになり、その他の成分は０となる。
【０１０５】
【数６４】
2005331440.tif 000073
具体的に、マトリックスで書き表してみると、下記数６５のような形となる。
【０１０６】
【数６５】
2005331440.tif 000074
また、被測定光の複素ベクトルｇのノルムを最小化するという先験情報を利用する場合には、マトリックスＣは単位マトリックスで表される。
【０１０７】
また、光位相分布の計測を表面形状計測などの分野に利用する場合に、表面形状の概形は、前もって分かっていることがある。この表面形状の概形といった情報は、光位相分布の概形が既知であることと等価であり、被測定光に関する先験情報として利用できる。すなわち、設計変数の参照解ｇｒｅｆが与えられているときに、適切化関数として、下記数６６を用いることができる。
【０１０８】
【数６６】
2005331440.tif 000075
本発明では、上述したような被測定光に関する先験情報を利用することが可能である。位相分布同定逆問題の設定により、適切化関数を使い分けたり、組み合わせることで、位相分布同定逆問題の悪条件性を克復することができる。
【０１０９】
なお、本発明で利用可能な先験情報として、上述した２種類の情報に限定されることなく、被測定光に関する他の先験情報を利用することが可能であることは言うまでもない。
＜２−２−３＞最適化問題のサイズ変換方法
本発明では、観測量として光波検出センサー（ＣＣＤ）により得られる光波の強度分布情報は、膨大なデータとなる。このため、非常に大規模な最適化問題を解くことになり、膨大な計算時間が必要となってしまう問題が生じる。この問題を克復するため、本発明の光位相分布測定方法では、設定した大規模な最適化問題に対応するために、次のような最適化問題のサイズ変換方法を用いる。
【０１１０】
つまり、設計変数ｇ∈ＣＮ、観測量|ｆ|∈ＲＮ、光学系の変換マトリックス[Ｈ]∈ＣＮ×Ｎの場合に、下記数６７で表す最適化問題があるとき、この最適化問題（以下、元の最適化問題と称する）を直接解くことはせず、図９に示されるように、元の最適化問題の設計変数ｇを隣り合ういくつかの点を代表する設計変数ｇ´∈Ｃｎ（ただしｎ＜Ｎとする）で表現する。
【０１１１】
【数６７】
2005331440.tif 000076
図９では、設計変数ｇの隣り合う４点を代表するものを設計変数ｇ´としている。まず、この設計変数ｇ´で表現されたサイズの小さな最適化問題を解く。観測量は元の最適化問題の|ｆ|の隣り合ういくつかの点を平均化して、サイズの小さな最適化問題に対応する観測量|ｆ´|∈Ｒｎを作るものとする。図９では、|ｆ|の隣り合う４点づつの平均をとって並べたものを|ｆ´|∈Ｒｎとしている。また、光学系の変換マトリックスは、サイズの小さな最適化問題に対応するように、前述した光学系による変換マトリックスの計算式により、新たに[Ｈ´]∈Ｃｎ×ｎを用意する。つまり、元の最適化問題より、サイズの小さな下記数６８で表す最適化問題を解く。
【０１１２】
【数６８】
2005331440.tif 000077
上記数６８で表すこのサイズの小さな最適化問題の推定解を、元の最適化問題の設計変数の初期値（つまり、初期推定解）として用い、元の最適化問題を解く。
【０１１３】
上述した最適化問題のサイズ変換方法の手順を図を用いて説明する。例えば、一次元の１００点の最適化問題（つまり、元の最適化問題）があるときに、１０点ずつをまとめて、まず、１０点の最適化問題（つまり、サイズの小さな最適化問題）を解くことによって、図１０に示されるような設計変数ｇ´の推定解が得られる。次に、この推定解を元のサイズの１００点に拡張し、図１１に示されるように、１００点の最適化問題の設計変数ｇの初期推定解（初期値）とする。最後に、この初期推定解（初期値）を用い、この１００点の最適化問題を解くことによって、図１２に示されるように、１００点の最適化問題の設計変数ｇの推定解が得られる。
【０１１４】
要するに、本発明では、このような最適化問題のサイズ変換方法を用いることで、元の最適化問題において、解の探索を真の解のより近くから開始することができ、直接解くのに比べ少ない反復回数で収束することが期待できる。
【０１１５】
そして、本発明では、このような最適化問題のサイズ変換方法を用いることで、計算時間がかからないサイズの小さな最適化問題で、おおまかな推定を行い、大きな計算コストがかかる元の最適化問題での計算時間を短縮することができる。
【０１１６】
さらに、本発明の最適化問題のサイズ変換方法では、上述したような二段階の推定だけでなく、サイズの異なる幾つかの最適化問題を用意しておき、サイズの小さな最適化問題から徐々にサイズの大きな最適化問題へと数段階の推定を行うことも可能である。
＜３＞本発明を適用した光位相分布測定結果
本発明に係る光位相分布測定系及び光位相分布測定方法の有効性を確認するために、数値実験を行った。本発明の光位相分布測定系及び光位相分布測定方法を用いて、光位相分布の幾つかの測定結果（同定例）を以下のように示す。
＜３−１＞被測定光が一次元分布の場合
被測定光が一次元分布する場合での測定結果（つまり、例１〜例８)を示し、本発明における先験情報の利用および最適化問題のサイズ変換方法の有効性を確認する。
＜３−１−１＞位相分布同定逆問題設定
まず、本発明の光位相分布測定系に用いる光学系は、二種類とし、つまり、回折系とレンズ焦点系を用いるものとする。ここで、レンズは直径無限大の理想的なレンズであるとする。
【０１１７】
入力光および観測光をｎ個に離散化した場合に、回折系の変換マトリックス[Ｈ１]の(ｉ,ｊ)成分は、＜１＞数理モデルで説明したように、下記数６９で表される。
【０１１８】
【数６９】
2005331440.tif 000078
ここで、ｋは波数で、Ｒは光の入力面と観測面との距離である。被測定光の波長λ＝０.５×１０−６(ｍ)、つまり、波数ｋ＝２π／λ＝１.２５６×１０−７(１／ｍ)とする。また、Ｒ＝０.１(ｍ)とし、入力面ｘおよび観測面Ｘの離散点の間隔Δｘ,ΔＸ＝１.０×１０−５(ｍ)とする。
【０１１９】
レンズ焦点系の変換マトリックス[Ｈ２]は、直径が無限大の理想的なレンズを用いた場合に、下記数７０に表されるように、焦点でボケのない完全な像を結ぶため単位行列となる。
【０１２０】
【数７０】
2005331440.tif 000079
次に、ある正解の入力光ｇを設定し、これを元に、数４３で表す光位相分布測定系の観測方程式を用いて、模擬的に測定された強度分布|ｆ１|（回折系）、|ｆ２|（レンズ焦点系）を作成する。これを有効数字３桁に丸め、誤差を含ませたものを観測量とする。
＜３−１−２＞本発明の光位相分布測定方法を用いた測定結果（同定結果）
まず、本発明において、先験情報の利用の有効性を確認する。
【０１２１】
入力面および観測面を１００点に離散化した例１を図１３〜図１５、例２を図１６〜図１８、例３を図１９〜図２１、例４を図２２〜図２４にそれぞれ示す。図１３、図１６、図１９及び図２２は、被測定光の強度である観測量|ｆ１|,|ｆ２|を表し、図１４、図１７、図２０及び図２３は、先験情報を利用せずに、本発明を用いて同定した結果を示し、図１５、図１８、図２１及び図２４は、『複素振幅分布（つまり、位相）は滑らかに分布する』という先験情報を利用して、本発明を用いて同定した結果を示している。
【０１２２】
例１、例２及び例３は、被測定光の強度分布がガウス分布の例である。先験情報を用いないで本発明を用いて同定した結果である図１４、図１７、図２０を見ると、中央部では同定精度が比較的に良いが、被測定光の強度が微弱な両端部では精度があまり良くないことが分かる。続いて、『複素振幅分布は滑らかに分布する』という先験情報を用いて本発明を用いて同定した結果である図１５、図１８を見ると、全体に渡りうまく同定できていることがよく分かり、本発明において先験情報の利用の有効性の確認ができた。
【０１２３】
しかし、例３において、『複素振幅分布は滑らかに分布する』という先験情報を用いて本発明を用いて同定した結果である図２１に示されるように、被測定光の強度が微弱な領域で大きな位相差が存在するときは両端部までうまく同定できていない。強度微弱領域における設計変数値の挙動は、目的関数値にあまり影響を与えないため、解が収束しづらいものと推測できる。
【０１２４】
そこで、例４のように、図２２に示されるような被測定光の強度分布を一様なものとして、本発明の光位相分布測定方法を用いて、光位相分布同定を行った。すると、図２３、図２４のように全体に渡り精度良く同定することができた。先験情報を用いて本発明を用いて同定した結果である図２４の方が僅かに精度が良いようである。
【０１２５】
このように、被測定光の強度を一様な分布としても、例５の被測定光および観測光を１０００点に離散化し、位相分布が複雑な形状をしている場合の同定例では、図２６を見ると、先験情報を用いないで本発明を用いて同定した場合に、全く見当はずれな解が出てきていることが分かる。このような悪条件な位相分布同定逆問題においても、先験情報を用いて本発明を用いて同定した結果である図２７を見ると、解の変動を抑え、うまく同定できていることが良く分かる。光の強度が一様な分布の場合でも、本発明において、先験情報の利用の有効性の確認ができた。
【０１２６】
次に、最適化の計算が収束するまでの反復回数の面から、本発明において、先験情報の利用の有効性を確認する。下記表１は、例１〜例５において、先験情報を用いない場合と先験情報を用いた場合の最適化の計算が収束するまでにかかった反復回数を示している。
【０１２７】
【表１】
2005331440.tif 000080
表１から、例１、例２、例３及び例５では、先験情報の利用は、同定精度の向上と計算時間の短縮に貢献していることが分かる。また、例４では、先験情報を用いた場合も先験情報を用いない場合も最終的な同定精度はほぼ同程度であるが、やはり先験情報の利用は計算時間を短縮させている。このことから、計算時間の面でも、本発明において、先験情報の利用の有効性の確認ができた。
【０１２８】
以上より、『複素振幅が滑らかに分布する』という先験情報は、解を安定化させ、同定精度の向上及び計算時間の短縮に有効であることが確認できた。
【０１２９】
次に、本発明において、最適化問題のサイズ変換方法（以下、二段階の同定手法とも称する）の利用の有効性を確認する。
【０１３０】
二段階の同定手法を適用した例６を図２８〜図３０、例７を図３１〜図３３、例８を図３４〜図３６にそれぞれ示す。例６、例７及び例８において、入力面および観測面は、それぞれ１０００点ずつに離散化する。例６、例７及び例８において、二段階の同定手法は、次のように適用される。
【０１３１】
まず、１０００点の最適化問題を１００点の最適化問題に変換し、そして、サイズの小さな１００点の最適化問題を解く。次に、この１００点の最適化問題での推定解を１０００点の最適化問題の初期推定解として用い、１０００点の最適化問題を解くようにする。
【０１３２】
図２８、図３１、及び図３４は、観測量である被測定光の強度分布を示し、図２９、図３２、及び図３５は、１００点の最適化問題の推定解（つまり、二段階の同定手法の一段階目で得られた位相分布）を示し、図３０、図３３、及び図３６は、元の１０００点の最適化問題の推定解（つまり、二段階の同定手法の二段階目で得られた位相分布）を示す。
例６、例７及び例８からは、１００点の最適化問題で解をおおまかに推定した後に、元の１０００点の最適化問題で詳細に解を同定している様子がよく分かる。
【０１３３】
下記表２は、二段階の同定手法を用いた場合と二段階の同定手法を用いない場合の最適化の計算が収束するまでにかかる反復回数を示している。
【０１３４】
【表２】
2005331440.tif 000081
表２では、二段階の同定手法を用いた場合の最適化の計算が収束するまでにかかる反復回数は、一段階目(１００点の最適化問題)と二段階目(１０００点の最適化問題)と分けて表示してある。なお、設計変数１００点の１回の反復計算にかかる時間は、１０００点のそれと比べ非常に短いため無視できる程度である。
【０１３５】
表２より、二段階の同定手法を用いることで、最適化の計算が収束するまでの反復回数を大幅に減少させたことがよく分かる。本発明の光位相分布測定方法において、二段階の同定手法の利用は、計算時間の短縮に有効であることが確認できた。勿論、この二段階の同定手法を多段階に拡張して用いることで、大規模な最適化問題にも対応することができる。
＜３−２＞被測定光が二次元分布の場合
被測定光が二次元分布の場合でも、本発明に係る光位相分布測定方法が有効であることを測定結果（つまり、例９〜例１２）により確認する。
＜３−２−１＞位相分布同定逆問題設定
まず、本発明の光位相分布測定系に用いる光学系は、二種類とし、つまり、回折系とレンズ焦点系を用いるものとする。ここで、レンズは直径無限大の理想的なレンズであるとする。
【０１３６】
入力光および観測光をｎ×ｎ個に離散化した場合に、回折系の変換マトリックス[Ｈ１]の（ｎ(ｉ−１)＋ｊ,ｎ(ｋ−１)＋ｌ）成分は、下記数７１で表される。
【０１３７】
【数７１】
2005331440.tif 000082
ここで、ｋは波数で、Ｒは光の入力面と観測面との距離である。被測定光の波長λ＝０.５×１０−６(ｍ)、つまり、波数ｋ＝２π／λ＝１.２５６×１０−７(１／ｍ)とする。また、Ｒ＝０.１(ｍ)とし、入力面(ｘ,ｙ)および観測面(Ｘ,Ｙ)の離散点の間隔Δｘ,Δｙ,ΔＸ,ΔＹ＝１.０×１０−５(ｍ)とする。
【０１３８】
また、レンズ焦点系の変換マトリックス[Ｈ２]も、一次元分布の場合と同様に、直径が無限大の理想的なレンズを用いた場合に、数７０に表されるように、焦点でボケのない完全な像を結ぶため単位行列となる。
【０１３９】
次に、ある正解の入力光ｇを設定し、これを元に、数４３で表す光位相分布測定系の観測方程式を用いて、模擬的に測定された強度分布|ｆ１|（回折系）、|ｆ２|（レンズ焦点系）を作成する。これを有効数字３桁に丸め、誤差を含ませたものを観測量とする。
＜３−２−２＞本発明の光位相分布測定方法を用いた測定結果（同定結果）
入力面および観測面を４０×４０点に離散化した例９を図３７〜図４０、例１０を図４１〜図４４、例１１を図４５〜図４８、例１２を図４９〜図５２にそれぞれ示す。例９、例１０、例１１及び例１２は、二段階の同定手法を適用せず、『複素振幅が滑らかに分布する』という先験情報を利用して、本発明の光位相分布測定方法を用いた測定結果（同定結果）である。
【０１４０】
図３７、図４１、図４５及び図４９は、被測定光の強度である観測量|ｆ１|の分布を示す。なお、観測量|ｆ２|は、観測量|ｆ１|と一様分布となるため、その分布図を省略する。図３８、図４２、図４６及び図５０は、被測定光の正解の位相分布である。図３９、図４３、図４７及び図５１は、本発明の光位相分布測定方法を用いて同定された被測定光の位相分布である。図４０、図４４、図４８及び図５２は、正解および本発明による同定結果の位相分布のｙ＝２０(点)における断面を示している。
【０１４１】
例９や例１０のような被測定光の正解の位相が滑らかに分布している場合に、本発明を用いた同定結果を見ると、全領域で良く正解と一致していることがよく分かる。また、例１１や例１２のような被測定光の正解の位相が急激に変化するような分布をしている場合においても、本発明を用いた同定結果を見ると、多少のばらつきが見られるが全領域でうまく同定できていることがよく分かる。また、離散点数を増やし空間周波数を上げることで、このようなばらつきを抑えることができる。
【０１４２】
以上のように、例１〜例１２を通して、本発明に係る光位相分布測定系及び光位相分布測定方法の有効性が確認された。
＜４＞本発明のまとめ
上述したように、本発明に係る光位相分布測定方法は、異なる複数の光学系を通過した被測定光の強度画像から、被測定光の位相分布を同定することを最大の特徴としている。
【０１４３】
要するに、本発明では、まず、異なる複数の光学系と光波検出センサーとから構成される光位相分布測定系を用いて、同一の被測定光を光学特性が既知であるこれら複数の光学系にそれぞれ通過させ、強度と位相を変調し、被測定光の強度分布を複数測定する。そして、得られた被測定光の強度分布の情報及び光学系の光学特性に基づいて、光位相分布測定系の観測方程式をたて、この観測方程式から被測定光の位相分布を同定する位相分布同定逆問題を設定し、設定された位相分布同定逆問題を非線形の最適化問題として定式化し、この最適化問題を解くことで、被測定光の位相分布を同定するようにしている。
【０１４４】
なお、上述したように、準ニュートン法を用いて、本発明の最適化問題の目的関数の最小化を行うことが好ましいが、それに限定されることなく、最適化問題に対する数値解法であれば、他の方法を用いて本発明の最適化問題を解くことができる。
【０１４５】
また、本発明では、光位相分布測定系において、適切な光学系を任意の個数設定することや、光位相分布測定方法において、被測定光の位相分布に関する先験情報を与えることによって、光位相分布測定系の観測方程式の悪条件性を克復し、解を安定化させ、唯一解を求めることが可能となる。
【０１４６】
更に、本発明の光位相分布測定方法で設定する位相分布同定逆問題において、観測量として、光波検出センサー（ＣＣＤ）により得られる被測定光の強度分布情報は、膨大なデータになるため、非常に大規模な最適化問題を解くことになる。本発明の光位相分布測定方法では、最適化問題のサイズ変換方法を用いて、このような大規模な最適化問題に対応するようにしている。
【図面の簡単な説明】
【０１４７】
【図１】光波の複素数の二次元配列を説明するための模式図である。
【図２】本発明に係る光位相分布測定系の測定原理及び数理モデルを説明するための模式図である。
【図３】回折系による光波の変換の数理モデルを説明するための模式図である。
【図４】結像レンズ系による光波の変換の数理モデルを説明するための模式図である。
【図５】結像光学系の点像応答関数ＰＳＦを示す図である。
【図６】回折格子による光波の変換の数理モデルを説明するための模式図である。
【図７】Armijoの基準の幾何学的意味を説明するための図である。
【図８】光波が二次元（４×４）で分布している場合の行方向又は列方向の隣り合う点を説明するための模式図である。
【図９】最適化問題のサイズ変換を説明するための模式図である。
【図１０】サイズの小さな最適化問題の設計変数ｇ´の推定解を示す図である。
【図１１】元の最適化問題の設計変数ｇの初期値を示す図である。
【図１２】元の最適化問題の設計変数ｇの推定解を示す図である。
【図１３】例１において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図１４】例１において、先験情報を利用せずに本発明を用いて同定された被測定光の位相分布図である。
【図１５】例１において、先験情報を利用して本発明を用いて同定された被測定光の位相分布図である。
【図１６】例２において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図１７】例２において、先験情報を利用せずに本発明を用いて同定された被測定光の位相分布図である。
【図１８】例２において、先験情報を利用して本発明を用いて同定された被測定光の位相分布図である。
【図１９】例３において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図２０】例３において、先験情報を利用せずに本発明を用いて同定された被測定光の位相分布図である。
【図２１】例３において、先験情報を利用して本発明を用いて同定された被測定光の位相分布図である。
【図２２】例４において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図２３】例４において、先験情報を利用せずに本発明を用いて同定された被測定光の位相分布図である。
【図２４】例４において、先験情報を利用して本発明を用いて同定された被測定光の位相分布図である。
【図２５】例５において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図２６】例５において、先験情報を利用せずに本発明を用いて同定された被測定光の位相分布図である。
【図２７】例５において、先験情報を利用して本発明を用いて同定された被測定光の位相分布図である。
【図２８】例６において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図２９】例６において、二段階の同定手法の一段階目で得られた位相分布図である。
【図３０】例６において、二段階の同定手法の二段階目で得られた位相分布図である。
【図３１】例７において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図３２】例７において、二段階の同定手法の一段階目で得られた位相分布図である。
【図３３】例７において、二段階の同定手法の二段階目で得られた位相分布図である。
【図３４】例８において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図３５】例８において、二段階の同定手法の一段階目で得られた位相分布図である。
【図３６】例８において、二段階の同定手法の二段階目で得られた位相分布図である。
【図３７】例９において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図３８】例９において、被測定光の正解の位相分布図である。
【図３９】例９において、本発明を用いて同定された被測定光の位相分布図である。
【図４０】例９において、本発明を用いて同定された被測定光の位相分布図のｙ＝２０(点)における断面図である。
【図４１】例１０において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図４２】例１０において、被測定光の正解の位相分布図である。
【図４３】例１０において、本発明を用いて同定された被測定光の位相分布図である。
【図４４】例１０において、本発明を用いて同定された被測定光の位相分布図のｙ＝２０(点)における断面図である。
【図４５】例１１において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図４６】例１１において、被測定光の正解の位相分布図である。
【図４７】例１１において、本発明を用いて同定された被測定光の位相分布図である。
【図４８】例１１において、本発明を用いて同定された被測定光の位相分布図のｙ＝２０(点)における断面図である。
【図４９】例１２において、光位相分布測定系を用いて測定された被測定光の強度分布図である。
【図５０】例１２において、被測定光の正解の位相分布図である。
【図５１】例１２において、本発明を用いて同定された被測定光の位相分布図である。
【図５２】例１２において、本発明を用いて同定された被測定光の位相分布図のｙ＝２０(点)における断面図である。
【図５３】従来の干渉計による光位相分布測定方法を説明するための模式図である。
【図５４】干渉縞画像の例を示す図である。
【図５５】従来のシャック・ハートマンセンサーによる光位相分布測定方法を説明するための模式図である。
【図５６】ハートマンパターンの例を示す図である。
【図１】
2005331440.tif 000083
【図２】
2005331440.tif 000084
【図３】
2005331440.tif 000085
【図４】
2005331440.tif 000086
【図５】
2005331440.tif 000087
【図６】
2005331440.tif 000088
【図７】
2005331440.tif 000089
【図８】
2005331440.tif 000090
【図９】
2005331440.tif 000091
【図１０】
2005331440.tif 000092
【図１１】
2005331440.tif 000093
【図１２】
2005331440.tif 000094
【図１３】
2005331440.tif 000095
【図１４】
2005331440.tif 000096
【図１５】
2005331440.tif 000097
【図１６】
2005331440.tif 000098
【図１７】
2005331440.tif 000099
【図１８】
2005331440.tif 000100
【図１９】
2005331440.tif 000101
【図２０】
2005331440.tif 000102
【図２１】
2005331440.tif 000103
【図２２】
2005331440.tif 000104
【図２３】
2005331440.tif 000105
【図２４】
2005331440.tif 000106
【図２５】
2005331440.tif 000107
【図２６】
2005331440.tif 000108
【図２７】
2005331440.tif 000109
【図２８】
2005331440.tif 000110
【図２９】
2005331440.tif 000111
【図３０】
2005331440.tif 000112
【図３１】
2005331440.tif 000113
【図３２】
2005331440.tif 000114
【図３３】
2005331440.tif 000115
【図３４】
2005331440.tif 000116
【図３５】
2005331440.tif 000117
【図３６】
2005331440.tif 000118
【図３７】
2005331440.tif 000119
【図３８】
2005331440.tif 000120
【図３９】
2005331440.tif 000121
【図４０】
2005331440.tif 000122
【図４１】
2005331440.tif 000123
【図４２】
2005331440.tif 000124
【図４３】
2005331440.tif 000125
【図４４】
2005331440.tif 000126
【図４５】
2005331440.tif 000127
【図４６】
2005331440.tif 000128
【図４７】
2005331440.tif 000129
【図４８】
2005331440.tif 000130
【図４９】
2005331440.tif 000131
【図５０】
2005331440.tif 000132
【図５１】
2005331440.tif 000133
【図５２】
2005331440.tif 000134
【図５３】
2005331440.tif 000135
【図５４】
2005331440.tif 000136
【図５５】
2005331440.tif 000137
【図５６】
2005331440.tif 000138