1. Introduction

Phase-shifting interferometry (PSI) has been widely used for precision measurements of surface profile [1, 2]. Today the most common phase-shifting interferometer is in Fizeau style due to its simplicity, i.e. no other components than the reference surface and the test surface are involved in establishing the optical path difference [3, 4]. To implement high precision spherical surface testing with the PSI technique, much efforts have been done by former researchers.

As early as in 1980, Moore and Slaymaker found that the phase shift caused by the axial translation of the spherical reference surface is not uniform but field dependent [5]. Fortunately, as the well-designed modern phase shifting algorithms have high tolerance to the scaling error in the nominal phase shift increment, in most cases the accuracies achieved in spherical surface testing are acceptable under laboratory controlled conditions [4, 6, 7]. Modifications on the interferometer were also suggested to implement the uniform phase shift, by replacing the mechanical phase shifting with the wavelength tuning phase shifting or inserting a flat reference surface isolated from the transmission sphere [6, 7].

Misalignment of the spherical test surfaces is another error source. Traditionally, misalignment aberrations are removed by subtracting Zernike tilt and defocus terms from the wavefront data [8], which is only valid for the test surface with small numerical aperture (NA) or small misalignment amounts [9]. The computer aided alignment methods [10] may be adapted to tackle with the residual high-order aberrations due to misalignment. Recently, Wang et al. suggested a wavefront difference algorithm to remove the misalignment aberrations, based on the high-order approximation of the optical path difference [9]. Similarly, Yuan et al. proposed a (n + 4)-frame strategy to remove the modified defocus aberrations in the wavefront data, which was retrieved by the standard phase-shifting algorithms [11]. With their idea it is unnecessary to adjust the interferogram to the null-fringe pattern, which makes the spherical surface testing much easier.

From the conventional view, vibration is boring for precision temporal phase-shifting interferometry, as it would introduce uncertainty in the phase shifts. To solve this problem, Novak et al. proposed simultaneous phase shifting interferometry schemes [12], which made the testing systems more expensive and complex. While some other researchers focused on design of phase-shifting method, so as to retrieve wavefront data from interferogram sequence with large unknown tilt phase shifts [13–17]. With such methods it is feasible to realize precision surface testing without the use of any phase-shifting devices, but by effectively utilizing the phase shifts induced from vibration. In other words, like every coin has two sides, sometimes we can even benefit from vibration for surface testing. However, the methods proposed in [13–17] are only applicable to testing of flat surfaces or spherical surfaces with small NA. Recently, Deck put forward a vibration-insensitive iterative phase-shifting method based on a physical model of the cavity, and it has been validated in spherical surface testing with NA values as high as 0.77 [1]. In that method, the initial estimation for the wavefront data is acquired with the standard phase shifting algorithms. However, we find that the convergence performance of that method relies on the initial estimation values. Under severe vibration conditions the phase shifts would deviate too far from the nominal status, and the quality of the initial estimations would be low. As a result, it may fail to retrieval the wavefront, as shown in our simulation results in section 4.

Imaging distortion of the interferometer is also a factor to influence the testing accuracy, which is often underestimated in spherical surface testing. Selberg showed that the radial imaging distortion would create errors in the spherical surface testing results, appearing as third order coma, third order spherical and fifth order spherical terms [18]. If the fringe patterns are not well-nulled, such errors should not be neglected in the case of precision fast spherical surface testing [3, 18]. ZYGO Corporation has patented the technique for in situ calibration of the imaging distortion [19]. As the commercial interferometers are well designed and the polished spherical surfaces always have small fabrication errors, it may be adequate to model the imaging distortion with the primary aberration (third-order aberration) theory [20], in the same way as in [21, 22].

In this paper a general method is proposed for interferometric spherical surface testing with unknown phase shifts. With modifications on the physical model of the interferometer cavity proposed by Deck [1], the influence of imaging distortion on the wavefront retrieval is also taken into account. Meanwhile, the proposed method has good convergence stability, as a result of the much emphasis paid on finding the reasonable initial values.

The retrieved wavefront results by the proposed method would have indetermination in the global sign if no prior knowledge about the phase shifts is available, which is common to any asynchronous approaches [15]. Such indetermination also exists in any single-frame phase extraction methods [23–27] given no prior information. However, in practice we need not concern too much on it, as elimination of the global sign ambiguity is always much easier than retrieval of the wavefront shape.

The rest of the paper is organized as follows: the physical model of the interferometer cavity is given in section 2, the principle of the proposed method is provided in detail in section 3, while the simulation and experimental results are reported in section 4 and section 5, respectively. Finally, discussions are drawn in section 6.

2. Physical model of the interferometer cavity

Recently, Deck has proposed a physical model of the Fizeau interferometer cavity for spherical surface testing [1] as follows

Equation (1) can also be rewritten as

In Deck’s model, the imaging distortion of the interferometer has been ignored. However, as shown in [3, 18] the errors induced by the image distortion cannot be neglected in precision fast spherical surface testing, excepting that the fringe patterns are well-nulled. Therefore, in this work the model is modified to include the influence of imaging distortion. Considering that the commercial interferometers are well designed, we simply model the imaging distortion as $\rho \text{'}=\rho +\Delta p\approx \rho \left(1+\epsilon {\rho }^2\right)$. Substituting it into Eq. (5), we will get

3. Principle of the proposed method

The method we proposed includes three parts. Firstly, the interferogram sequence is normalized, and processed frame by frame with the single-frame phase extraction technique to retrieve the individual phase maps as well as the phase shifts between the frames. These phase shift data are subsequently fitted. The initial estimations for the wavefront, direct current and interference contrast terms are calculated with the least-squares method. Secondly, the wavefront data is further refined in an iterative way, by fitting the sequence of interferograms to the physical model of the interferometer cavity with the linear regression technique (which is also named as MPSI technique in [1]). Thirdly, the aberrations introduced by the surface misalignment and the imaging distortion are removed, so as to achieve the wavefront result related to the actual surface profile.

Figure 1 shows a block diagram of the proposed method.

 

Fig. 1 Block diagram of the proposed method

Download Full Size | PDF

In this section, all the involved parameters will be investigated mainly in the polar coordinate system for illustration. Actually, in our programs all the data are processed in the Cartesian coordinate system.

3.1 Initial estimation of the phase shifts and wavefront

In this section we just aim to solve the initial estimations for the phase shifts and wavefront data $\Theta \left(t_1\right)$. It is adequate to adopt Deck’s model of cavity. Besides, precision prior knowledge of the parameters$g,\epsilon$and NA is not required for the moment, i.e. we can simply assume$g=0,\epsilon =0$ (considering that typical laser Fizeau glass reference surfaces have an intensity reflectivity of $r_{ref}\approx 4\%$, $g\ll 1$always holds); even if $NA$is unknown, we can still find a reasonable initial guess for it.

The captured interferograms are processed as follows. Firstly, all the interferograms are normalized [28] individually and phase maps are extracted from them with the single-frame phase extraction technique. Specifically, the regularized frequency-stabilizing method proposed in [27] is chosen, allowing for its simplicity, and ability to deal with single interferogram with both open and closed fringe patterns. Then the ambiguities in the global sign of these phase maps are eliminated with the aid of the least-squares method [29]. Meanwhile, the initial guess for NA can be achieved by inspecting the fitting results of these phase shift data, if it is unknown. Hence we can effectively extract the unknown phase shift information between the interferograms, and further make a more reliable initial estimation for the wavefront data.

(a) Processing the interferogram sequence with the regularized frequency-stabilizing method

After processing the whole interferogram sequence frame by frame with the regularized frequency-stabilizing method [27], we will obtain a series of estimated unwrapped phase maps $\widehat{\Theta }\left(x,y,t_i\right),i=1,2,\dots N$. As pointed out in section 1, all these phase maps would have ambiguities in the global sign, i.e. $\widehat{\Theta }\left(x,y,t_i\right)+2n_i\pi \approx \pm \Theta \left(x,y,t_i\right),i=1,2,\dots N$, where $n_i$ is an integer.

(b) The ambiguities of the phase shift data due to the sign ambiguities in the estimated phase maps

We can rewrite Eq. (2) as

It would be a direct way to get the phase shift data, by taking the difference between the estimated phase maps $\widehat{\Theta }\left(x,y,t_i\right),i=1,2,\dots N$. However, due to the existence of global sign ambiguities in these phase maps, there will also exist ambiguities in the phase shift data.

Considering the case of three phase-shifted interferograms $I_1,I_i,I_j(1<i<j\le N)$, if we assume $\widehat{\Theta }\left(\rho ,\theta ,t_1\right)\approx \Theta \left(\rho ,\theta ,t_1\right)$, there are totally four candidate solutions for $\delta \left(\rho ,\theta ,t_i\right)$ and$\delta \left(\rho ,\theta ,t_j\right)$, i.e.

(c) Elimination of the global sign ambiguities in the estimated phase maps

Considering that $g\ll 1$always holds, the interferogram intensity can be approximately written as

When there are more than three phase-shifted interferograms, the global signs of all these estimated phase maps$\widehat{\Theta }\left(\rho ,\theta ,t_i\right)$, $i=1,2,\dots ,N$, can be determined in the same way. As a result, the ambiguities in the phase shift data are also eliminated. It is worth noting that if we alternatively assume$\widehat{\Theta }\left(\rho ,\theta ,t_1\right)\approx -\Theta \left(\rho ,\theta ,t_1\right)$, the other set of phase shift estimations with inverse sign will be obtained.

(d) Fitting of the phase shift data and initial estimation for NA

In the proposed method the phase shift data should be fitted according to the physical model of the interferometer cavity, in order to improve their accuracy and facilitate the subsequent refinement process. In some cases NA of the test surface may be unknown. Then we have to make reasonable initial guess for it, so as to ensure the validity of the above fitting procedure. Here we put forward a simple way to implement it.

Referring to Eq. (6) and Eq. (8), $\delta \left(\rho ,\theta ,t\right)$ can be further expanded in the following form

When the values of NA and $\epsilon$are both unknown, the phase shift data $\widehat{\delta }\left(\rho ,\theta ,t\right)$ can be approximately fitted as

We can extract out partial phase shift data$\widehat{\delta }\left(\rho ,\theta ,t\right)$ within a square sub-region, e.g.$-r\le \rho \cos \theta \le r,-r\le \rho \sin \theta \le r$ ($r$ is recommended to be larger than 0.45). The extracted data can be denoted as ${\widehat{\delta }}_{sub}\left(\rho ,\theta ,t\right)$. Fit it with Eq. (13) over a set of sampled NA values in the potential range$\left[N{A}_{\min },N{A}_{\max }\right]$, and define a merit function

The proposed way to find $\widehat{NA}$ is quite straightforward to be understood, i.e. when $N{A}_{samp}\approx NA$we will have${\widehat{p}}_{\mathrm{mod}}\left(N{A}_{samp},t\right)\approx 0$, thus$C\left(N{A}_{samp}\right)\approx 0$; otherwise, we will have${\widehat{p}}_{\mathrm{mod}}\left(N{A}_{samp},t\right)\gg 0$, thus $C\left(N{A}_{samp}\right)\gg 0$, as a result of the nonorthogonality between the piston term and the term$\sqrt{1-{\rho }^{2}N{A}_{samp}{}^2}$. The proposed way is always effective unless$\Delta k_z\left(t\right)\approx 0$ holds for all the interference patterns, which is a very rare case with multiple frames and can be easily avoided.

Obviously, once $\widehat{NA}$ is determined, the estimation results for the fitted phase shifted data ${\widehat{\delta }}_{fit}\left(\rho ,\theta ,t\right)$, as well as the related fitting parameters$\Delta k_x(t)$, $\Delta k_y(t)$ and$\Delta k_z(t)$, are also obtained simultaneously during the above fitting process.

Up to now, we can get updated estimation results for$A\left(\rho ,\theta \right),V\left(\rho ,\theta \right),$and $\Theta \left(\rho ,\theta ,t_1\right)$with the least-squares method [29], by taking the whole interferogram sequence into account. All these estimation results (including the fitting parameters$\Delta k_x(t)$, $\Delta k_y(t)$ and$\Delta k_z(t)$for the phase shift data) will be adopted as initial values for the subsequent iterative refinement procedure, as illustrated in section 3.2.

3.2 Refinement of wavefront data

To refine the wavefront data ${\widehat{\Theta }}_{\text{w}}\left(\rho ,\theta ,t_1\right)$ the MPSI technique proposed in [1] is adapted, incorporating with the modified physical model of interferometer cavity (see Eq. (6)). The refinement operation relates to an iterative optimization process, and each iteration consists of two steps.

At the first step, the estimates for the time-independent variables$A\left(\rho ,\theta \right)$,$V\left(\rho ,\theta \right)$, and$\Theta \left(\rho ,\theta ,t_1\right)$ (see Eq. (4)) are treated as constants. Fitting each interferogram to a linearized form of the mathematical model can provide new estimates for the time-dependent variables $\Delta k_x(t)$, $\Delta k_y(t)$, and $\Delta k_z(t)$, as well as the constants $\epsilon$and $NA$(if they are unknown). We can define

While at the second step the time-dependent variables $\Delta k_x\left(t\right)$, $\Delta k_y\left(t\right)$, and$\Delta k_z\left(t\right)$ are treated as constants. New estimates for the time-independent variables$A\left(\rho ,\theta \right)$, $V\left(\rho ,\theta \right)$, and $\Theta \left(\rho ,\theta ,t_1\right)$, as well the constant $g$(if it is unknown) can be solved by the least squares method as indicated in [1]. Here, the estimate for$\Theta \left(\rho ,\theta ,t_1\right)$ is wrapped within the$[0,2\pi )$ range.

The iterations continue until a termination condition is satisfied.

It is worth noting that in [1] $\Phi$ is selected as the interested unknown wavefront quantity instead of $\Theta \left(t_1\right)$. This minor difference may have influence on the final wavefront results, as shown in section 4.1.

3.3 Removement of aberration errors due to misalignment and imaging distortion

We can unwrap the refined wavefront data${\widehat{\Theta }}_{\text{w}}\left(\rho ,\theta ,t_1\right)$ with some specific unwrapping method, such as the Goldstein method [30], and the outcome is denoted as$\tilde{\Theta }\left(\rho ,\theta ,t_1\right)$. When we test spherical surfaces (especially the ones with high NA values), adjusting the interferogram to the null-fringe pattern usually require much skill and patience [9]. Besides, it would be difficult and subjective to judge the best alignment status, if the test surface has large fabrication errors. In other words, it is almost unavoidable to include some misalignment aberrations in$\widehat{\Theta }\left(\rho ,\theta ,t_1\right)$, i.e. $\Delta \left(\rho ,\theta ,t_1\right)\ne 0$ (see Eqs. (4) and (5)). Meanwhile, as stated in section 1, the imaging distortion will introduce extra aberrations in presence of misalignment errors.

In [9] Wang et al. have analyzed the misalignment aberrations, based on the high-order approximation of the optical path difference. However, in that paper the imaging system of the interferometer is assumed distortionless. In this section, we will extend their analysis to incorporate the influence of imaging distortion. Additionally, in [9] two separate measurements with the test surface located at different defocus positions are required to eliminate the misalignment aberrations, if no prior knowledge of NA is available. While in our proposed method it would not be a problem, as NA can be solved with the technique illustrated in section 3.2.

As proposed in [9], the aberrations purely induced by misalignment (see Eq. (5)) can be well approximated with a series of Zernike polynomials [9, 31] as follows

We find that when the factors of misalignment and imaging distortion are both considered, the induced aberrations (see Eq. (6)) can be well approximated in a similar way as follows

Besides, based on Eq. (25) it is able to retrieve the values of$k_x\left(t\right)$, $k_y\left(t\right)$and$k_z\left(t\right)$$k_z\left(t\right)$ from the coefficients$a_1\text{'}$,$a_2\text{'}$ and$a_3\text{'}$, respectively. However, the linear ratios between$a_1\text{'},a_2\text{'},a_3\text{'}$ and $k_x\left(t\right),k_y\left(t\right),k_z\left(t\right)$ actually vary when changing the data sampling densities, i.e. number of CCD pixels, as a result of the discrete sampling and the non-uniform mapping relationship between the Cartesian and polar coordinates. Therefore, in this paper these linear ratios are determined numerically with the practical data sampling configuration.

Now,$\tilde{\Theta }\left(\rho ,\theta ,t_1\right)$ is fitted with Zernike polynomials [9, 31]. Based on the obtained fitting coefficients$a_1\text{'}\left(t_1\right)$,$a_2\text{'}\left(t_1\right)$ and$a_3\text{'}\left(t_1\right)$, the values of$k_x\left(t_1\right),k_y\left(t_1\right)$, and$k_z\left(t_1\right)$ can be retrieved. Then the aberrations due to misalignment and imaging distortion can be calculated via Eq. (6), and separated from$\tilde{\Theta }\left(\rho ,\theta ,t_1\right)$, i.e. $\widehat{\Phi }\left(\rho ,\theta ,t_1\right)=\tilde{\Theta }\left(\rho ,\theta ,t_1\right)-\widehat{\Delta }\left(\rho ,\theta ,t_1\right)$. Meanwhile, the values of$k_x\left(t_i\right),k_y\left(t_i\right),k_z\left(t_i\right),1<i\le N$can be determined as by-products, since$\Delta k_x\left(t\right)$, $\Delta k_y\left(t\right)$, and $\Delta k_z\left(t\right)$ have been solved in section 3.2.

If the test surface has large fabrication errors, its ancillary aberrations due to imaging distortion should also be considered to ensure the measurement accuracy. For the double pass of test beam in the Fizeau cavity,$0.5\widehat{\Phi }\left(\rho ,\theta ,t_1\right)$, i.e.$0.5\widehat{\Phi }\left(x,y,t_1\right)$ in the Cartesian coordinate system would be determined as the wavefront result related to the actual surface profile.

4. Simulation and discussion

To validate the effectiveness of the proposed method, a series of computation simulations have been carried out.

4.1 Performance in comparison with Deck’s method

(a) Comparisons in convergence performance

Thirteen frames of interferograms with pixels of $500\times 500$ are generated according to Eqs. (1), (2), (4) and (6) by setting the parameters as follows:$A\left(x,y\right)=145\exp \left(x^2+y^2\right)$, $V\left(x,y\right)=100\exp \left(x^2+y^2\right)$and$\Phi \left(x,y\right)=1.5\left(x^2+y^2\right)y-2y^3$ (i.e. a scaled trefoil term in$y$ direction), where$x^2+y^2\le 1$; $\epsilon =0.02$, $g=0.1$ ($r_{ref}=0.4\%$, $r_{test}=25\%$), NA = 0.67 (F/0.75), and values for the time-dependent variables$k_x\left(t_i\right)$,$k_y\left(t_i\right)$,$k_z\left(t_i\right)$,$1\le i\le 13$ are listed in Table 1. The first three frames of simulated interferograms, as well as the simulated surface profile (simulated fabrication errors), i.e. $0.5\Phi \left(x,y\right)$ are shown in Fig. 2.

Table 1. The simulated time-dependent variables [k_x(t_i)_,k_y (t_i)_,k_z(t_i)], i = 1,2,…13, and the retrieved ones by the proposed method (units in radians)

View Table | View all tables in this article

 

Fig. 2 The simulated interferograms and surface profile. (a-c) The first three frames of simulated interferograms. (d) The simulated surface profile. The data shown in (d) are in radians.

Download Full Size | PDF

We have retrieved the surface profile with the proposed method assuming that the constants$\epsilon$, g, NA are all unknowns. The retrieval results are given in Fig. 3. As shown in Fig. 3(a), the initial guess value of NA is 0.674, where the cost function (see Eq. (14)) achieves its minimum normalized value. The initial estimation result for $\Theta \left(x,y,t_1\right)$is shown in Fig. 3(b), and its related error is shown in Fig. 3(c), from which it can be seen that the initial estimation result has been very close to the truth. When ignoring the imaging distortion, the errors of the retrieved surface profile would be $0.24rad$ and $0.03rad$ in the PV (peak-to-valley) and RMS (root-mean-square) measures, respectively, as shown in Fig. 3(d). From it we can find substantial residual coma aberrations, which is in agreement with the predictions given in section 3.3. Figure 3(e) shows a plot of the PV errors in retrieved surface profile versus iteration number by the proposed method. Specifically, after 200 iterations the errors of the retrieved surface profile are as small as$8.2\times {10}^{-4}rad$and$2.3\times {10}^{-4}rad$in PV and RMS measures, respectively, as shown in Fig. 3(f). Meanwhile, the final retrieved values of $\epsilon$, g, NA by our proposed method are 0.02, 0.10 and 0.67, respectively. The retrieved values of the time-dependent variables$k_x\left(t_i\right)$, $k_y\left(t_i\right)$, $k_z\left(t_i\right)$are given in Table 1, and it demonstrates that the retrieval errors are all below $0.001rad$.

 

Fig. 3 The retrieved results by the proposed method. (a) The initial guess for NA. (b) The initial estimation of the wavefront data related to time t₁. (c) The related residual errors in the wavefront data shown in (b). (d) The retrieved surface profile when imaging distortion is neglected. (e) A plot of the PV and RMS errors in retrieved surface profile versus iteration number by the proposed method. (f) The residual errors in the retrieved surface profile by the proposed method. The data shown in (b-c) and (d, f) are in radians.

Download Full Size | PDF

Besides, we find that the initial guess for NA is nontrivial. For instance, if we arbitrarily choose the initial guess value as smaller than 0.60 (or larger than 0.73) in this simulation, it even fails to retrieval the surface profile. We think that like many other iterative optimization techniques, the convergence performance of MPSI technique also relies on the initial values. As the quality of initial guess for NA becomes worse, the quality of initial estimations for wavefront data as well as phase shift data would degenerate until they are outside of the convergence domain. However, it would not be a problem to the proposed method, as we can always obtain high-quality initial guess for NA by the way proposed in section 3.1.

For comparisons, we have also provided the retrieval results by Deck's method [1]. Here, the number of iterations used in MPSI technique is set as 2000, and the problem has been relaxed, i.e. the constants g and NA are assumed to be known quantities. The initial estimations for the direct current and interference contrast terms are acquired in the same way as in [1], while the initial estimation for the wavefront data is solved by a standard 13-frame algorithm with$\pi /4$ on-axis phase increment [32]. The initial estimation results are shown in Figs. 4(a)-4(c), and the erroneously retrieved surface profile is given in Fig. 4(d). In this case, it never converges to the right solution even with more iterations. Anyway, this instance reveals again the importance of reasonable initial estimations to the reliable convergence.

 

Fig. 4 The retrieved results by the method proposed in [1]. (a-c) The initial estimations of the direct current and interference contrast terms, as well as the wavefront data. (d) The erroneously retrieved surface profile in wrapped form. The data shown in (c-d) are in radians.

Download Full Size | PDF

The above simulation results demonstrate that our proposed method can well retrieve the spherical surface profile with unknown phase shifts, and exhibits much more enhanced convergence performance compared to the original MPSI technique. Additionally, as explained in section 1, the retrieved wavefront results by the proposed method would have indetermination in the global sign if no prior knowledge about the phase shifts is available. To reasonably evaluate the estimation performance of the proposed method, the simulation results (as well as the experimental results shown in section 5) are all assigned with the correct signs.

(b) The non-uniqueness of wavefront retrieval results by Deck’s method

At the first glance, Deck's method [1] seems to be more compact in data processing flow than ours, i.e. in [1] the remove of misalignment aberrations seemingly has been integrated in MPSI optimization procedure, while in our proposed method compulsory post processing is required on the refined wavefront data. Actually, it is not definitely true.

For simplicity, we assume the imaging system of the interferometer has none distortion, and$\left\{\widehat{\Phi }\left(x,y\right),k_x\left(t_i\right),k_y\left(t_i\right),k_z\left(t_i\right)\right\},1\le i\le N$is one solution provided by the method in [1]. Then it is easy to infer that$\left\{\widehat{\Phi }\text{'}\left(x,y\right),k_x\text{'}\left(t_i\right),k_y\text{'}\left(t_i\right),k_z\text{'}\left(t_i\right)\right\}$ would also be the possible solution for the same problem (on the premise of successful convergence), with

To validate it, a simulation example is given below by setting the parameters as follows. $A\left(x,y\right),B\left(x,y\right)$, and$\Phi \left(x,y\right)$are defined as same as in the above simulations; $g=0.1,\epsilon =0,NA=0.67$(here NA is assumed as a known quantity); $k_x\left(t_i\right)=k_y\left(t_i\right)=0$,$k_z\left(t_i\right)=\pi /4\times \left(i-7\right)+\pi \times \left[rand(1)-0.5\right]$,$1\le i\le 13$, where rand() is a built-in function of MATLAB, to generate uniformly distributed pseudorandom numbers. Then one set of retrieved surface profile by the method in [1] is shown in Fig. 5. As shown in Fig. 5(b), the retrieved surface profile contains misalignment-style aberration errors, which equals to $0.082rad$in PV measure. To pursue super-high precision measurement results (i.e. at nanometer scale), the wavefront data should be further calibrated, as shown in Fig. 5(c).

 

Fig. 5 The retrieved results. (a-b) The retrieved surface profile and the corresponding error map, by the method proposed in [1]. (c) The resultant error map after removement of the numerical misalignment-style aberrations. All the data shown are in radians.

Download Full Size | PDF

In contrast, in our proposed method$\Theta \left(x,y,t_1\right)$is selected as the interested unknown wavefront quantity instead of$\Phi \left(x,y\right)$, during the MPSI optimization procedure. Therefore, the outcome of the iterative optimization procedure would be unique, so is the final retrieved surface profile.

4.2 Some considerations for stability of the proposed method

The stability is an important aspect of the phase extraction methods. This section will make some discussions about it.

(a) Sensitivity to the amplitude of the unknown phase shifts

Vibration perturbation is one of the main factors to introduce uncertainties in phase shift amounts. Some references [32, 33] have investigated the stability performance of the specific phase-shifting methods under vibration conditions, by calculating a spectrum of vibration sensitivity. However, in those papers the vibrations are regarded as ideal single-frequency events. Such assumptions may be violated to some degree in real-world applications. On the other hand, the iterative nature of the proposed method makes it difficult to derive analytic formulas for vibration sensitivity. Additionally, in some modern interferometers camera shuttering has been employed to eliminate motion-induced contrast blurring [1], then the contrast variations among interferograms are always negligible. Therefore, in more general sense we may model the effect of low-frequency environment perturbations as random unknown phase shifts with certain amplitudes. Based on that, some simulation results are presented here to compare the stability performance of the proposed method with Deck’s method [1].

The setting of simulation parameters excepting the variables$k_x\left(t_i\right),k_y\left(t_i\right),k_z\left(t_i\right)$is the same as in section 4.1(b). While

 

Fig. 6 Comparisons of the stability performance between Deck’s method and our proposed method. (a) The successful convergence rate of Deck’s method with different amplitudes of random unknown phase shifts along z-axis. (b) The successful convergence rate of our proposed method.

Download Full Size | PDF

Such phenomenon is straightforward to be understood. As$r_z$get larger, the phase retrieval performance of the standard phase-shifting algorithms would decline, i.e. the quality of the initial estimation results for Deck’s method would decrease, so the stability of Deck’s method is prone to be influenced. In contrast, in the proposed method the initial estimation results are obtained based on the single-frame phase extraction technique. So its stability performance has no direct collection with the values of $r_z$, but depends on the effectiveness of the adopted single-frame phase extraction technique. For simplicity, here we have only shown the cases where random phase shifts exist merely along the z-axis. Similar phenomena can be observed if random phase shifts exist also along the other axes.

(b) The requirements on the captured interferograms

As stated in the previous paragraph, the effectiveness of the adopted single-frame phase extraction technique is crucial for the stability of the proposed method. As a variant of the frequency-guided sequential demodulation method [25], the regularized frequency-stabilizing method has poor performance in the area where the fringe pattern is sparse. Since critical points have low total local frequencies, the frequency-guided strategy has been proposed to reduce the possibility of error propagating from critical point [26]. In this work, we have fused the frequency-guided strategy into the regularized frequency-stabilizing method. However, in the extreme cases, i.e. the captured interferograms are at near-null status, such improved regularized frequency-stabilizing method would be insufficient to provide reliable phase extraction results.

In view of the above consideration, in this work the interferograms with unknown phase shifts are all captured at status with noticeable misalignment amounts, so as to ensure the effectiveness of the single-frame phase extraction procedure.

(c) The applicability range when NA is unknown

As stated in section 4.1, when NA is unknown we need to provide a reasonable initial guess for it, so as to ensure the successful convergence of the proposed method. We find that by the way proposed in section 3.1, reliable initial guess for NA can always be obtained (the typical guess errors are less than 0.01), if its truth value is no less than 0.30. Such initial guess precision for NA is far more than adequate, in consideration of the stability of the proposed method.

When the truth value of NA is quite small, e.g. NA = 0.15, the initial guess for it would be unreliable. It is because that in this case the essential difference between the piston term and $\sqrt{1-{\rho }^{2}N{A}^2}$(see section 3.1) would be minor and the fitting coefficients for them would be sensitive to the residual errors in the initially estimated phase shift data. On the other hand, the test spherical surface with quite small NA has been very close to the flat surface. Therefore, in such cases we can simply set the initial guess for NA as 0, and the stability of the proposed method can be maintained.

Another direct solution to the spherical surface testing problems with unknown NA would be trials of several different initial guess values (i.e. the phase extraction results are determined by observing the convergence process), but at the cost of increasing computation load. From the above discussion, we can conclude that the lack of NA information would not be an obstacle for high precision surface testing.

5. Experimental results

We have also carried out experiments with the 4-inch ZYGO Fizeau interferometer, to validate the experimental performance of the proposed method. The cavity consisted of an F/0.68 spherical surface (also provided by ZYGO, with reflectivity of 4%) illuminated through an F/0.75 transmission sphere (with reflectivity of 4%). To introduce unknown phase shifts in the recorded interferograms, the active vibration isolation workstation was turned off, and the 5-axis mount (loading the test surface) was manually adjusted in random directions along with randomly unknown adjustment amounts. Twelve frames of interferograms are captured. The central 92% subregion of the measured data is selected to analyze in the experimental validation. The average filter (4 × 4) is applied on the retrieved surface profile, to restrain the random system noise [9, 11].

The first frame of interferograms (with unknown phase shifts) is shown in Fig. 7(a). By our proposed method, the retrieved distortion coefficient $\epsilon$ is equal to$-0.013$, and the retrieved values for the time-dependent variables are listed in Table 2. The retrieved surface profile is shown in Fig. 7(b), which equals to $0.268rad$ and $0.028rad$ in the PV and RMS measures, respectively. At the same time, the retrieval result when neglecting the imaging distortion is also shown in Fig. 7(c), which equals to $0.558rad$ and $0.087rad$ in the PV and RMS measures, respectively. From Fig. 7(c), the dominant coma components can be visually perceived. Although not shown here, in this case it fails to converge by the original MPSI technique [1], since the initial estimates deviates too much from the truth.

 

Fig. 7 (a) The first frame of the captured interferograms with unknown phase shifts. (b-c) The retrieved surface profile by the proposed method, with (b) or without (c) considering the imaging distortion. (d) The first frame of interferograms with well controlled phase shifts. (e) The retrieved surface profile based on the results given by ZYGO’s MetroPro software. (f) The difference between (b) and (e). The data shown in (b-c) and (e-f) are in radians.

Download Full Size | PDF

Table 2. The retrieved time-dependent variables$\left[k_x\left(t_i\right),k_y\left(t_i\right),k_z\left(t_i\right)\right]$, $i=1,2,\mathrm{...},12$ by the proposed method (units in radians)

View Table | View all tables in this article

We also tested the same surface under vibration-isolated conditions with a calibrated PZT, and the outcome is taken as the reference result. The first frame of the thirteen captured interferograms (with well controlled phase shifts) is shown in Fig. 7(d). With the results provided by the ZYGO Metropro software, we can obtain the retrieved surface profile as shown in Fig. 7(e) after removing the aberrations due to misalignment and imaging distortion.

The difference between the retrieval results given by the proposed method and the MetroPro software is shown in Fig. 7(f), which equals to $0.262rad$ and $0.020rad$ in the PV and RMS measures, respectively, i.e. an experimental accuracy of $0.0032\lambda$ RMS has been achieved. As the test surface has high fabrication precision (which is comparable to the experimental accuracy of our method), the visual discrepancy between Figs. 7(b) and 7(e) is evident (as shown in Fig. 7(f)), however, we can still notice some existent similarity between them.

To check repeatability of the proposed method, another two independent measurements have been carried out with the same configuration. The retrieved results are shown in Fig. 8, where the first column represents the first frame of the captured interferograms including unknown phase shifts, respectively, the second column provides the retrieved surface profiles by the proposed method, and the third column shows the residual differences between these retrieved results and the retrieval result shown in Fig. 7(e). The residual difference shown in Fig. 8(c) equals to $0.281rad$ and $0.051rad$ in the PV and RMS measures, respectively, while the residual difference shown in Fig. 8(f) equals to $0.245rad$ and $0.021rad$ in the PV and RMS measures, respectively. In summary, the residual differences for all the three independent measurements equal to$0.0032\lambda$, $0.0033\lambda$and$0.0040\lambda$in the RMS measures, respectively.

 

Fig. 8 The retrieved results by the proposed method in another two independent measurements. (a, d) show the first frame of the captured interferograms including unknown phase shifts in the two independent measurements, respectively. (b, e) show the retrieved surface profiles by the proposed method. (c, f) show the residual differences between these retrieved results and the retrieval result shown in Fig. 7(e). The data shown in (b-c) and (e-f) are in radians.

Download Full Size | PDF

Some factors may account for the residual difference shown in Figs. 7(f), 8(c), and 8(f). As in our experiment the interferograms containing unknown phase shifts were captured in a relative long time duration, the laser beam would meet some temporal intensity fluctuations, i.e. parameters$A(x,y)$, $B(x,y)$would be no longer time-independent in the strict sense, which will introduce some errors into the retrieved wavefront results. As suggested in [1], the addition of an intensity monitor may provide enough information to correct such errors. Moreover, the dynamic variations of the testing environment, such as the air turbulence will also lead to some differences between measurements.

6. Conclusion

We have proposed a general method for spherical surface testing with unknown tilt phase shifts, based on an improved physical model of the interferometer cavity taking into account both the rigid cavity motions and the imaging distortion. With the regularized frequency-stabilizing method and the least-squares method, we can obtain reasonable initial estimations for the wavefront and phase shift data. Such wavefront result is further refined in an iterative way, by fitting the sequence of interferograms to the physical model of the interferometer cavity with the linear regression technique. As reasonable initial estimations are available, the convergence stability of the mentioned iterative refinement procedure is greatly enhanced compared with the original MPSI technique. After removing the aberrations due to the surface misalignment and the imaging distortion, the surface profile can be retrieved. Both simulations and experiments have demonstrated the validity of the proposed method. The proposed method is applicable in a testing environment with low frequency and high amplitude vibration. With this method costly and accurate phase-shifting devices are no longer required for steady-state measurement, by effectively utilizing the phase shifts induced from environment perturbations.

In fact, the proposed method can be further extended for the real applications with the current framework. The first concern is about the imaging distortion. In some self-developed interferometer, distortion is no longer a one-dimensional function of the radius, but depends in addition on the azimuth angle [3]. In this case, the calibration for the image distortion is needed, and the physical model of the interferometer cavity should be modified accordingly. Secondly, the selection of single-frame phase extraction method is flexible, i.e. some other reliable and faster ones can also be considered, as a result of the rapid development on it. Additionally, in the proposed method reliable phase unwrapping operations are essential for fitting the phase shift data and removing the misalignment aberrations. However, if the captured interferograms contain much noise, phase unwrapping may become difficult. Then more advanced wavefront data filtering and phase unwrapping techniques [34] should be adopted.