計測自動制御学会産業論文集
Vol. 8, No. 13, 108/112 (2009)

位相勾配検出による干渉縞画像の２次元周波数推定
北 川 克 一*Two-Dimensional Frequency Estimation for Fringe Analysis by Phase Gradient Detection
Katsuichi KITAGAWA*
A new frequency estimation technique is proposed, which enables us to estimate two-dimensional frequencies of interferometricfringes with high accuracy and low computational cost. It is accomplished by phase gradient detection, where phases are calculated
by a local model fitting algorithm for carrier pattern analysis. The algorithms used and experimental results are presented.

Key words: frequency estimation, phase gradient, interferometry, fringe analysis, surface profiler, single-shot
1.

はじめに近年，半導体や液晶など様々な産業分野において，ナノメートルオーダの表面凹凸形状を精度良く測定したいという要求が高まってきている．光干渉を用いた表面形状測定法は，速度や測定精度，保守性の観点から最も有望な計測手法である．
代表的な光干渉計測法である位相シフト法 1) では，干渉計の測定面と参照面の相対距離を変えながら複数枚の干渉画像を撮像し，その情報から表面形状を推定する．この方法では，複数Fig. 1 Optics of spatial carrier interferometry.の画像を撮像する必要があるため，振動などの外乱のある環境下では精度が大きく低下するという問題がある．この解決策として，一枚の画像から表面形状を求めるワンショット計測法が提案されている．その代表的なものは，参照面を傾斜させてキャリア縞を生成させる方法であり，空間キャリア縞法 2)〜8)と呼ばれている（Fig. 1）．この方法によって得られる１枚の干渉縞画像（Fig. 2）からフーリエ変換法 3) ，空間位相同期法 4)〜6) ，Fig. 2 Interferogram with carrier fringes.

局所モデル適合法(Local Model Fitting 法；LMF 法) 7)8) などによ

約 470nm が測定可能な最大段差となる．

り表面形状が求められる．
しかし，いずれの方法においても，位相シフト法と同様，隣この方法を延長して，さらなる測定レンジ拡大を目的に，緑接画素間に光源波長の 1/4 以上の段差が存在する場合には，正色（中心波長 530nm）の LED 照明を追加し，カメラの G 信号しい位相アンラッピング（位相接続 ）ができないという問題がも利用する 3 波長ワンショット計測法を検討した 11)．この方式ある．この問題解決のために，筆者らは，青色（中心波長 470nm）の実現には，RGB 信号間のクロストーク補正，高精度な周波数と赤色（中心波長 627nm）の 2 個の LED 照明とカラーカメラを推定，3 波長アンラッピングなど，新たな課題の解決が必要で用いる 2 波長同時撮像系を実現し，350nm の段差の 2 波長ワンあった．得られた最終的な概略フローを Fig. 3 に示す．本報は，．この方法では，撮像したカラーこのフローのなかの周波数推定について述べるものである．ここ

画像を B,R 成分に分離し，各画素における位相を局所モデル適で言う周波数とは，試料面と参照面の相対傾斜角と波長により決合法により求め，等価波長法，あるいは，それを拡張した次数まる x 方向，y 方向のキャリア縞周波数である．

ショット測定に成功した

9) 10)

決定法を用いて，2 波長アンラッピング（注 1）を実施し，高さにここで，周波数推定の目的を述べる．干渉縞画像の位相計算変換している．使用した 2 波長の等価波長は，1877nm であり，には，先に筆者らが開発した局所モデル適合法（付録Ⅰに概要を述べる）を利用する．この方法は，正弦波状モデル関数に含まれ

東レエンジニアリング（株）エレクトロニクス事業本部
開発センター 滋賀県大津市大江 1‑1‑45* R&D Center, Electronics Division, Toray Engineering Co.,Ltd.
*るパラメータ（振幅，周波数，位相，直流成分）のうち，周波数を既知とし，残る 3 個のパラメータを最小自乗法で求める．よって，位相計算の前に周波数の推定が必要である．
さらに，周波数推定の必要精度について考察する．周波数に誤(注 1)ここで言うアンラッピングとは，複数波長を利用するものであり，通常の隣接画素情報を利用するものとは異なる．本報では，後者を「位相接続」と表記して区別する．

差があると，2.1 節で述べるように位相が直線的に変化する．この誤差が測定結果に及ぼす影響を考えると，1 波長法，および，アンラッピングに等価波長法を用いた 2 波長法の場合には，表(Received April 16, 2009)

面形状が傾斜して測定されるだけであって，問題にならない．し108

f = f 0 + (dφ ′ ( x) / dx) / 2π(5)′ここで，位相勾配 dφ / dx は，最低 2 点のデータから求められるので，最も簡易的には，基準平面内の異なる 2 点(x1,x2)における位相 φ ′1 ， φ ′2 から，f = f 0 + (1 / 2π )(φ '2 −φ '1 ) /( x2 − x1 )(6)

により，周波数 f が得られることになる．得られた周波数を初期値にして，推定を繰り返すことにより推定精度を向上することもできる．なお，正しい位相勾配を求めるためには，通常，位相接続 が必要であるが，この問題については 4.2 節に述べる．

2.22 次元推定

つぎに，この方法を 2 次元に拡張すると，(1)，(2)，(4) 式はそれぞれ次式のようになる．
Fig. 3 Flowchart of three-wavelength
single-shot interferometry.
かし，3 波長法の場合には，アンラッピングに合致法を用いて縞次数を決定するため，僅かの位相誤差が次数の変化として拡g ( x, y ) = A cos(2πf x x + 2πf y y + φ ( x, y ))(7)g ' ( x, y ) = A′ cos(2πf x 0 x + 2πf y 0 y + φ ' ( x, y ))(8)

φ ' ( x, y ) = 2π ( f x − f x 0 ) x + 2π ( f y − f y 0 ) y + φ ( x, y )(9)また，(5) 式に相当する周波数推定式は，次式のようになり，大する可能性がある．付録Ⅱに述べる考察によれば，周波数推x 方向と y 方向の周波数が同時に推定できる．

f x = f x 0 + (dφ ′ ( x, y ) / dx) / 2πf y = f y 0 + (dφ ′ ( x, y ) / dy ) / 2π定の必要精度は 0.4%程度であり，かなりの高精度が要求される．
正弦波信号の周波数推定法としては，自己相関法，フーリエ変換法，Prony 法 12) など多くの推定手法が提案されている．しかし，精度が高く，計算負荷の低い実用的な手法は見当たらな

3.

い．さらに，Fig. 2 に示すように縞が座標軸に平行な場合には，

(11)

計算機実験

本提案手法の検証のため，計算機実験を行なった．

その軸方向の縞周波数がゼロに近く，1 次元的な推定法では高精度な推定が困難である．

3.1

筆者らは，位相勾配を利用した周波数推定法を考案し，上記

1 次元推定実験

3.1.1 実験条件の問題を解決した 11) 13)．本報では，その推定原理，計算機実験観測データを関数 g ( x) = cos(2πfx) から生成した．ここで，真結果，実装上の問題点と解決策，実試料実験結果について述べの周波数 f =0.020 とし，g(x)には (-0.1,+0.1)の一様雑音を付加しる．

た．初期周波数 f0 =0.019（相対誤差 5%）とし，位相計算用デ2.

2.1(10)周波数推定原理ータサイズは，25 画素とした．

1 次元推定3.1.2 実験結果

本提案手法（位相勾配法と呼ぶ）は，周波数誤差と位相勾配

Fig. 4 に観測データとモデル関数（ただし，振幅 A′ を 1，位との関係を利用する．本手法は 2 次元の周波数推定に適用可能相φ′ (0)を 0 としている）を示す．フィッティングにより得られであるが，先ず 1 次元の周波数推定で原理を述べる．
Fig. 2 の水平方向を x 軸として，観測信号が次式で表されるとする．

g ( x) = A cos(2πfx + φ ( x))(1)

ここで，A は振幅，f が求める周波数，φ(x)は点 x における位相であって高さにより変化する．
この観測信号に，初期周波数推定値 f0 のモデル関数

g ' ( x) = A′ cos( 2πf 0 x + φ ' ( x))(2)をフィッティングして，位相φ'(x) を求める．この位相計算には，Fig. 4 Observed data and model function.

先に筆者らが開発した局所モデル適合法（付録Ⅰ参照）を利用する．
(1)(2) 式から次式が成立する．

2πfx + φ ( x ) = 2πf 0 x + φ ' ( x )よって，得られる位相 φ ′ (x) は，次式で表される．
φ ' ( x ) = 2π ( f − f 0 ) x + φ ( x )

(3)
(4)

ここで， φ (x) =一定と見なせる領域（基準平面領域と呼ぶ）を考えると，位相 φ ′ (x ) は x の 1 次式となり，その勾配が周波数誤差 f - f0 に比例する．よって，(4) 式の両辺を微分して得ら
Fig. 5 Phase gradient.

れる次式により，周波数推定ができる．

109

た位相φ′ (x)と回帰式を Fig. 5 に示す．位相勾配 0.00628 から，周波数推定値は，0.019 + 0.00628/2π = 0.02000 となり，誤差が0.00001 以下で真値と一致した．また，初期周波数の広い範囲で同様の結果が得られた．
3.22 次元推定実験

y 方向の周波数がゼロの場合の 2 次元周波数推定を行なう．
従来の 1 次元信号処理では周波数推定が極めて困難なケースである．

(a) Observed image

3.2.1 実験条件(b) Model image

Fig. 6 Two-dimensional frequency estimation.

観測データを関数 g ( x, y ) = cos( 2πf x x + 2πf y y ) + 1 から生成した．ここで，真の周波数を fx =0.2，fy =0.0 とし，g(x)には(-0.1,+0.1)の一様雑音を付加した．Fig. 6(a)に観測画像を示す．
また，初期周波数を fx0 =0.19（相対誤差 5%）, fy0 =0.01 とした．
Fig. 6(b)にモデル画像（ただし，振幅 A′ =1，位相φ′ (x,y)=0 としている）を示す．位相計算用データは，3×3 画素とした．
3.2.2 実験結果
フィッティングにより得られた位相φ′ (x,y)を Fig. 7 に，そのプロファイルを Fig. 8(a)(b)に示す．
x 方向, y 方向の位相勾配は，位相データに，平面 Z = aX + bY + c をフィッティングし，係数

Fig. 7 Two-dimensional phase map.

a,b から求められ,それぞれ 0.06275，-0.06298 (rad/pixel) であった．これから，x 方向周波数推定値 fx =0.19999，y 方向周波数推定値 fy =-0.00002 が得られた．推定を 10 回繰り返した時の平均値±標準偏差は, fx = 0.20000±0.00003, fy = -0.00002±0.00002 であった．また，初期周波数の広い範囲で同様の結果が得られることを確認した．干渉縞の周波数や方向に依らず，高精度で 2次元周波数推定のできることが示された．
4.
4.1

実装上の問題点と解決策(a) x-direction(b) y-direction

Fig. 8 Phase profiles along three selected lines.

推定フローと問題点

推定フローを Fig. 9 に示す．点線は，得られた推定値を初期値として推定を繰り返すループを示している．周波数補正値 ∆fx，∆fyの絶対値が収束判定しきい値ε 以下になることを収束条件とする．
しきい値ε は，周波数推定の必要精度が 0.4%程度であること（付録Ⅱ参照）を考慮して設定する．
4.2

実装上の問題点

2.1 節で述べたように，正しい位相勾配を求めるためは，通常，位相接続が必要である．必要な例を Fig. 10(a)に示す．位相勾配は小さいが，位相値が±π近傍にあるため，不連続点が発生している．位相接続には多くの提案があるが，ノイズに強いロバストな位相接続は計算負荷が高い．そこで，位相接続を不要化する方法を検討した．
4.3

解 決 策 (1)− 原 点 シ フ ト 法

位相接続 を不要化するには，(9) 式により得られる位相値をゼロに近づけることが有効であり，これは以下に述べる「原点シフト」で実現できる．すなわち，Fig. 10(b)に示すように，モデル信号の位相を観測信号の位相と一致させると，求められる

Fig. 9 Flowchart of frequency estimation.

位相はゼロ近傍となる．実際には，観測データ信号の最大個所を求め，その点を座標原点として位相計算する．この方法によ数の小領域で粗推定を行なう．このとき，前節で述べた「原点り，Fig. 10(b)に示すように，位相接続 が不要となる．

シフト法」を併用することにより，位相接続 を不要化する．Fig.
11(a)における黄色実線が，原点と座標軸を示している．つぎに，

4.4

解 決 策 (2)− 2 ス テ ッ プ 法

第 2 ステップとして，粗推定で得られた周波数の平均値を初期初期値の誤差が大きいと，位相勾配が大きくなり，さらに推値とし，Fig. 11(b)の点線で示すように，広い領域の周波数推定定領域が大きいと，位相が±πを超える可能性がある．そこで，を実施する．ここでも，「原点シフト法」を採用する．この方先ず第 1 ステップとして，Fig. 11(a)に青色点線で示すように複法により，位相接続 が不要となった．

110

(a) With non-optimized origin(b) With optimized origin

Fig. 10 Effects of origin optimization.

Fig. 12 Experimental setup.

(a) 1st step: Rough estimation(b) 2nd step: Fine estimation

Fig. 11 Two-step estimation.
5.
5.1

実試料実験

実験方法

Fig. 13 Observed red image of 1µm step.

実験装置を Fig. 12 に示す．３波長ワンショット測定用に製作されたもので，光源は 3 色(RGB)LED 照明装置であり，干渉画像はカラーカメラで撮像される．1μm 標準段差試料を撮像したカラー干渉縞画像の R 成分を Fig. 13 に示す．中央下部の矩形部が段差である．
5.2

実験結果

Fig. 13 の上部に白線で囲んだ矩形領域（400×200 画素）を基準平面領域とし（注 2），初期値を変化させながら，周波数推定した．位相計算用データは，25×5 画素とした．その結果，広い範囲の初期値に対して，安定に収束し，x,y 方向の周波数推定値
Fig. 14 Estimated frequencies versus initial values.

として，fx =0.030064， fy =0.000361 が得られた．これは，水平方向が 512 画素の画面内の縞本数に換算すると，x 方向 15.393本, y 方向 0.185 本となる．Fig. 14 は 1 回の推定における（すな次元周波数の同時推定が可能という特徴がある．また，周波数わち，推定を繰り返さない場合の）初期値と推定結果との関係がゼロ近傍の場合でも適用可能である．計算機実験と実試料実を示す．回帰係数が 0.0029 であることから，1 回の推定で周波験により，提案手法の妥当性・有効性を確認した．さらに，実数誤差が約 3/1,000 に縮小されることが分かる．収束判定しき

装上の問題点として，位相接続 を取り上げ，これを不要化するい値εを 0.4%と設定した場合には，1 回の推定で収束し，推定方法を考案した．の繰り返しは不要になっている．推定に要する時間は，PC（Pentium 1.6GHz）において，23ms であった．また，基準平面領域の位相計算を全画素ではなく，100×100 画素に間引いても参 考 文 献

十分な精度が得られ，この場合の計算時間は 8ms であった．
6.

1)

まとめ

2)

本報告では，縞画像の新しい周波数推定法として，位相勾配3)法を提案した．画像の局所位相値が周波数誤差に依存することを利用して，周波数を推定する．精度が高く，計算負荷が軽く，24)(注 2)基準平面領域は，2.1 節で述べたように，位相が一定（すなわち，高さが一定）と見なせる領域に設定する．単一領域の必要はなく，画面内の複数個所に散在していても良い．

5)
6)
7)

111

J. H. Brunning et al.: Digital wavefront measuring interferometer for
testing optical surfaces and lenses, Appl. Opt., 13, 2693/2703 (1974).
加藤純一：実時間干渉じま解析とその応用,精密工学会誌，
64(9),1289/1293 (1998).
M. Takeda, H. Ina, and S. Kobayashi: Fourier-transform method offringe-pattern analysis for computer-based topography and interferometry,
J. Opt. Soc. Am., 72, 156/160 (1982).
S. Toyooka and M. Tominaga: Spatial fringe scanning for optical phasemeasurement, Opt. Commun. , 51, 68/70 (1984).
K. H. Womack: Interferometric phase measurement using spatialsynchronous detection, Opt. Eng., 23, 391/395 (1984).
J. Kato et al.: Video-rate fringe analyzer based on phase-shifting electronicmoire patterns, Appl. Opt., 36, 8403/8412 (1997).
M. Sugiyama, H. Ogawa, K. Kitagawa and K. Suzuki: Single-shot surface

profiling by local model fitting, Appl. Opt., 45, 7999/8005 (2006).
杉山将，松坂拓哉，小川英光，北川克一，鈴木一嘉：急峻な段差を持つ表面のワンショット形状計測法，精密工学会 2007 年度春季大会学術講演会講演論文集, 585/586 (2007).
9) K. Kitagawa, M. Sugiyama, T. Matsuzaka, H. Ogawa, and K. Suzuki:Two-wavelength single-shot interferometry, Proc. of SICE Annual
Conference 2007 in Takamatsu (計測自動制御学会学術講演会予稿集),
724/728 (2007).
10) 北川克一, 杉山将, 松坂拓哉, 小川英光, 鈴木一嘉：2 波長ワンショ
ット干渉計測，精密工学会誌, 75(2), 273/277 (2009).
11) 北川克一：3 波長ワンショット干渉計測，ViEW2008 ビジョン技術の実利用ワークショップ講演論文集, 5/10 (2008).
12) 灘谷演，久保和良：正弦波周波数推定法の評価，第 39 回計測自動制御学会学術講演会予稿集，109D-2 (2000).
13) 北川克一：位相勾配検出による干渉縞画像の２次元周波数推定，2009 年度精密工学会春季大会学術講演会論文集，169/170(2009).のゼロ点を画面中央に置くと x の最大値が 256 画素となるので，

8)∆h の概略最大値は次式で表される．
∆h ≅ (300∗256)(∆f)(A8)

合致法においては，各波長の位相から縞次数を変えて高さ候補値を求め，3 つの波長の高さ候補値が最も良く合致する組み合わせを求める．誤った次数の選択を避けるために必要な高さ誤差は，経験上，10nm 程度である．よって，許容される高さ誤差を 10nmとすると，許容される周波数誤差は(A8)式より，約 0.00013(1/pix)となる．これは，512 画素の画面内の縞本数に換算すると，0.067本に相当する．本報告の実験では，縞本数が 15 本程度なので，許容相対誤差は約 0.4%となる．

付録Ⅰ：位相計算法
縞画像の各点の位相計算に使用した局所モデル適合法（Local--------------------------------------------------------------［著 者 紹 介］

Model Fitting 法）7))8) の概要を以下に述べる．
各点の輝度 g ( x, y ) が次式で表されるとする．

北川 克一（正会員）g ( x, y ) = a( x, y ) + b( x, y ) cos(φ (x, y) + 2πf x x + 2πf y y )(A1)1964 年東京大学計数工学科卒．同年，東レ（株）
入社．1989 年より画像処理を応用した半導体検査機器の研究開発に従事．2000 年より東レエンジニアリング（株）技監．2001 年度計測
自動制御学会技術賞，ViEW2003 小田原賞，手島記念財団発明賞を受賞．計測制御エンジニア．

ここで， a( x, y ) は直流成分， b( x, y ) は振幅， φ ( x, y ) が位相である． f x , f y はそれぞれ x 方向，y 方向のキャリア縞周波数であり，既知とする．
局所モデル適合法では，各点の近傍で a ( x, y ) , b( x, y ) ,φ ( x, y ) が局所的に一定と仮定して，正弦波状モデル関数を次式で定義する．
g ( x, y ) = a + b cos(φ + 2πf x x + 2πf y y )(A2)

このモデル関数を，各点の近傍 n 点(n≧3)の輝度データに最小自乗適合する．しかし，上記のモデルは非線形であるので，その適合は計算負荷が高い．そこで， ξ c = b cos φ ， ξ s = b sin φの変数変換により，次式のように線形化する．
g ( x, y ) = a + ξ c φ c ( x, y ) + ξ s φ s ( x, y )ここで，()(A3)(φ c (x, y) = cos 2πf x x + 2πf y y ， φ s (x, y) = − sin 2πf x x + 2πf y y)である．
この結果，a,b,φ を求める問題は a, ξc, ξs を求める線形最小自乗問題に変換され，3 元連立 1 次方程式を解くことにより，a，ξc，ξs を算出することができる．位相φ は次式により求められる．

φ = arctan(ξs /ξc)

(A4)

付録Ⅱ：周波数推定の必要精度の考察
3 波長アンラッピングに合致法を用いて縞次数を決定する場合に必要とされる周波数推定精度について考察する．単純化のために，x 軸方向の 1 次元の計測を考える．また，誤差の絶対値に着目し，正負符号を無視する．すると，周波数誤差 ∆f による位相誤差∆φ は，(4) 式の変形により，次式で表される．
∆φ = 2π (∆f)x(A5)位相誤差による高さ誤差 ∆h は，波長をλとして，∆h = λ (∆φ)/4π(A6)で表される．(A5) (A6) 式から∆h = (λ x/2) (∆f)(A7)が得られる．
本報告の実験条件では，最大波長λが約 600nm であり，座標系112