- Email: [email protected]
Numerical simulation of Shack–Hartmann wavefront sensor to characterize the turbulence phase statistics
ARTICLE IN PRESS
Optics and Lasers in Engineering 44 (2006) 466–478
Numerical simulation of Shack–Hartmann wavefront sensor to characterize the turbulence phase statistics Zhiling Jiang, Yang Dai, Faquan Li, Xuewu Cheng, Shunsheng Gong Wuhan Institute of Physics and Mathematics, The Chinese Academy of Sciences, P.O. Box 71010, Wuhan, Hubei, 430071, China Received 17 November 2004; received in revised form 8 March 2005; accepted 5 April 2005
Abstract Shack–Hartmann wavefront sensor (SHWS) is one of the important parts in an adaptive optics system and its detection accuracy has considerable influences on the performances of the whole system. In this paper, a simulator of such type of wavefront sensor is built and tested by measuring the slope structure function of Kolmogorov turbulence simulated by random phase screens. The numerical results show that it works well to characterize the turbulence phase statistics. In combination with the noise simulation, this simulator will be a very useful tool to further investigate and optimize the system parameters to improve the detection accuracy and the performance of the whole system. r 2005 Elsevier Ltd. All rights reserved. Keywords: Numerical simulation; Hartmann wavefront sensor; Turbulence statistics characterization; Slope structure function
Corresponding author.
E-mail address: [email protected] (Z. Jiang). 0143-8166/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlaseng.2005.04.012
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
467
1. Introduction It is well known that the wavefront distortion resulting from atmospheric turbulence is the crucial factor limiting the resolution of the optical system. The adaptive optics technology is usually utilized to actively correct the received wavefront and improve the performance of the system [1]. The wavefront shearing interferometer and Shack–Hartmann wavefront sensor (SHWS) are the most popular devices used to measure the distorted wavefront. The latter has the advantages such as short response time, easy operation, and convenient to comprehend the wavefront modes and has been widely used in real-time adaptive optical systems, beam quality diagnosis and the test of the optical components. Characterization of the turbulence statistics is very important to understand its influences on optical system. Based on the wavefront slope measurement, Fried [2] presented the theoretical expression of slope structure function for circular aperture, which could be applied to the long-time exposure imaging by using a single aperture or differential image motion measurement by using a pair of subapertures [3,4]. Later, Dayton [5] and Nicholls [6], etc. introduced and discussed the measurement of turbulence statistics by using SHWS; Silbaugh [7] derived the expression of slope structure function of arbitrary power-law turbulence power spectrum and applied them to analyze the experimental results; Rao [8] further extended that theory and showed the comparison of the theoretical results and experimental measurements. All those research work showed that the SHWS could be conveniently applied to characterize the atmospheric turbulence statistics; however, few of studies have been carried out to investigate the influences of the sensor parameters such as size of the subaperture and sampling points as well as lenslet array on such characterization. Thus, such wavefront sensor simulator is essentially required. Due to the advantages of low cost and easy realization, numerical simulation technology has been widely applied to study the problems of the atmosphere turbulence [9–15], such as analysis of the performance of the adaptive optical systems, design and optimization of the system parameters, and simulation of the experiments that could not be easily or impossibly realized in practice. Recently, we have applied this technology to numerically analyze the influences of noise on centroid detection accuracy [16,17]. In this paper, we initially set up a SHWS simulator, and test its performance by using Kolmogorov phase screens generated by random numbers. The simulated results show that this simulator could be used to characterize the turbulence statistics and will be a useful tool to further investigate the influences of system parameters on measurement results. The main contents of this paper include: in Section 2, the theoretical considerations of the numerical simulation such as generation of the phase screens, model of wavefront sensor, automatic detection of spot centroid and characterization of turbulence parameters are presented; in Section 3, the simulation results and the discussion are given and the conclusion is outlined in Section 4.
ARTICLE IN PRESS 468
Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
2. Theoretical analysis and numerical treatment Fig. 1 shows the principle of atmospheric turbulence measurement by using SHWS. The waves affected by refractive index fluctuation of atmosphere incident on the entrance pupil and enter into the lenslet array after an optical system. The wavefront is split into many subapertures and imaged on the focal plane of the lenslet. A CCD is used to measure the centroid of the spot in each subaperture. Under the parabolic axis approximation, the equation describing the wave propagation along z direction in random medium could be simplified as [9] 2 qc q c q2 c 2ik þ 2k2 Dnc ¼ 0, qz q2 x q2 y
(1)
where c is the complex amplitude of the wave, k ¼ 2pn=l0 is the wave number, l0 is the wavelength in vacuum and n ¼ n0 þ Dn is the refractive index of air, n0 and Dn are the constant and the fluctuation of the refractive index, respectively. The popular method to solve the above equation is splitting-step Fourier transform algorithm, i.e. it splits the propagation distance into many segments with interval of Dzn and in the nth segment the wave solution could be represented by "
( Fn ¼ F F
1
Fn1 exp i
k2x þ k2y 4k
!# Dzn
) expðikDnDzn Þ exp i
k2x þ k2y 4k
! Dzn ,
(2) where F is the spectrum amplitude of the wave field, kx and ky are the transverse components of the wave vector, F and F 1 are Fourier and inverse Fourier transform operations. From this expression, we could see that the first and the third items represent the wave propagation in vacuum, and a phase item is added (the second item) to show the influence of the turbulence on wave field in this distance. Hence, the propagation of wave in random medium could be described by the propagation in vacuum and a series of discrete phase screens. If the nonisoplanatic effects are not taken into account, a single phase screen located at the entrance pupil of the optical system will be adequate to approximate the static turbulence [1].
Fig. 1. Principle scheme of atmosphere turbulence measurement.
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
469
Simulation of phase screens
Statistics of the phase screens
Wave propagation and imaging
Spot centroid detection
Turbulence characterization analysis Fig. 2. Flowing chart of turbulence simulation, Hartmann spots array generation and atmospheric turbulence phase statistics measurement.
According to the scheme shown in Fig. 1, the aspects that should be considered in numerical simulation are summarized as: generation of random phase screen, describing of the lenslet array, centroid detection and analysis of the simulation results. We will discuss each of them in the following sections. Fig. 2 shows the flowing chart of simulation. 2.1. Generation of phase screen The phase screen generated from the random numbers under certain conditions could simulate the statistical character of the atmospheric turbulence. Fast Fourier transform algorithm [9] is the most basic approach to generate the phase screens. However, in order to sample adequately low frequencies, the width of the phase screens needs to be in the order of the outer scale, which may be hundreds of meters, and due to the fact that the limited sampling interval and the coherent periodic property of discrete Fourier transform, the conventional algorithm needs a large number of sampling points and much computational time. This disadvantage could be overcome by adding high-order harmonics as shown in Refs. [12,14]. The random phase screen could also be well generated from Zernike functions [18], which has the advantage that the remaining statistical variance could be estimated from this simulation and provides useful information for correction system; but the number of Zernike functions increases rapidly with the size of phase screen, which raises the difficulty for simulation. Another limitation of this algorithm is that it only simulates the turbulence in circular region due to the definition of the Zernike function. In our simulation, we use the eigenvalue method and random mid-point displacement technique to generate the phase screen with a certain power spectrum [13].
ARTICLE IN PRESS 470
Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
The general phase structure function of turbulence statistics could be defined as [6] b2 x1 x2 ; y1 y2 , (3) Dðx1 ; y1 ; x2 ; y2 Þ ¼ gb r0 where gb is the generalized coefficient, r0 is an analogous quantity to describe the turbulence and b is the power-law exponent of the phase spectrum. If b is equal to 11 3, r0 is the coherent length of Kolmogorov turbulence r0 . If fðx; yÞ represents the phase screen (assuming that the mean component in the whole aperture has been removed), the covariance of the phase screen is given by Cðx1 ; y1 ; x2 ; y2 Þ 1 ¼ Dðx1 ; y1 ; x2 ; y2 Þ 2ZZ 1 0 0 0 0 0 0 Dðx1 ; y1 ; x2 ; y2 ÞW ðx1 ; y1 Þ dx1 dy1 þ 2 ZZ 1 0 0 0 0 0 0 Dðx1 ; y1 ; x2 ; y2 ÞW ðx2 ; y2 Þ dx2 dy2 þ 2 Z Z Z Z 1 0 0 0 0 0 0 0 0 0 0 0 0 Dðx1 ; y1 ; x2 ; y2 ÞW ðx1 ; y1 ÞW ðx2 ; y2 Þ dx1 dy1 dx2 dy2 . 2
ð4Þ
In practice, it requires to rearrange the screen into a one-dimensional vector to calculate the covariance: FðmM þ nÞ ¼ T½fðmDx; nDyÞ;
m; n ¼ 1; 2; . . . M,
(5)
where M is sampling points in rows and columns, T is the transform function. The covariance matrix could be written as C ¼ hFFT i.
(6)
The eigen-system of C could be used to describe the phase screen. Assuming that the eigenvalues of C are l ¼ ðl1 ; l2 . . . lMM Þ, and the eigenvectors are E ¼ ðE 1 ; E 2 . . . E MM Þ, phase screen could be simulated by using the random vector X where each element is independent and has its variance given by the corresponding element of l: F ¼ EX .
(7)
Due to the limitation of computation load, the number of sampling points is small and the random mid-point displacement technique is used to improve the resolution of the phase screen. According to the self-similarity of turbulence, the phase screen with different sizes could be achieved by adding a modified factor [12]. 2.2. The description of lenslet array The lenslet array is an important part of SHWS, which transfers the wavefront measurement to centroid detection. According to the wave optics, the lens could be
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
treated as a phase device, and the mathematic expression is k FL ðx; yÞ ¼ wðx; yÞ exp i ðx2 þ y2 Þ , 2f
471
(8)
where k is the wave number, f is the focal length and wðx; yÞ is the aperture function defined as ( 1; jxjp d2 ; yp d2 ; wðx; yÞ ¼ (9) 0; others; where d is the size of the lens. It is easy to obtain the intensity distribution in detection plane as Z 1 Z 1 2 k Iðxf ; yf Þ / c1 ðx; yÞwðx; yÞ exp i ðxxf þ yyf Þ dxdy f 1 1 2 ¼ FTfc ðx; yÞwðx; yÞg , ð10Þ 1
where c1 is the wave field distribution at the front of the lenslet array. From Eq. (10), we can see that the intensity distribution could be obtained by using fast Fourier transform. It should be noted that the optical system in front of the lenslet array is not taken into account in our present model because it is easily solved by changing the absolute size of the wave field and no other changes are required for the sampling points and the calculation matrix. Of course, a phase modification should be added into the wave field if the influence of the aberration resulting from the optical system is investigated. 2.3. Automatic detection of the spot centroid How to accurately measure the centroid is important for SHWS, which will determine the performance of the system. Different to the simulation, the detection region for each subaperture is not known in advance in practical measurement. To simulate the practical case, we adapt the FFT technique [19] to realize the division of subapertures and the detection of centroids. For a periodic interferometric image, its intensity could be written as, for example, gðx; yÞ ¼ aðx; yÞ þ bðx; yÞ cos½2pf 0 x þ fðx; yÞ.
(11)
For simpleness, we only discuss the case along x direction. Usually, the phase information could be obtained by using FFT technique (the typical instrument is shearing interferometer). In our simulation, we use this method to determine the frequency f 0 , which represents the average interval of the fringes. Eq. (11) could be re-written as gðx; yÞ ¼ aðx; yÞ þ cðx; yÞ expði2pf 0 xÞ þ c ðx; yÞ expði2pf 0 xÞ, where cðx; yÞ ¼ 12bðx; yÞ exp½ifðx; yÞ.
(12)
ARTICLE IN PRESS 472
Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
Fig. 3. Typical simulation result using FFT technique.
The spectrum distribution could be obtained if we apply FFT to Eq. (12): Gðf x ; yÞ ¼ Aðf x ; yÞ þ Cðf x f 0 ; yÞ þ C ðf x þ f 0 ; yÞ,
(13)
where f x is the frequency along x direction. We could find that A is located at the origin and C and C are located at f ¼ f 0 and f ¼ f 0 , respectively. f 0 could be determined and the spot interval in spatial could be known according to Fourier transform relationship. In our simulation, we firstly determine the centroid position of the whole wave field in detection region, and then divide the subapertures according to the transverse and longitudinal interval. A typical simulation result is shown in Fig. 3, where the left shows the interferometric image and the right shows the spectrum distribution, which has been amplified. 2.4. Character of turbulence statistics The statistical character could be achieved by using detection information of SHWS [5–8]. The slope structure function for SHWS is defined as [7] DaS^ ðX ; X 0 Þ ¼ h½S a^ ðX Þ Sa^ ðX 0 Þ2 i,
(14)
where S is the wavefront slope, and a^ is the normal direction (x or y). After a certain mathematical calculation, the final result is DaS^ ðDx; DyÞ ¼ gb d 2
b2 Z d du triðuÞf2½Dx2 þ ðu þ DyÞ2 ðb2Þ=2 r0
½ðDx 1Þ2 þ ðu þ DyÞ2 ðb2Þ=2 ½ðDx þ 1Þ2 þ ðu þ DyÞ2 ðb2Þ=2 2½jujb2 ð1 þ u2 Þðb2Þ=2 g,
ð15Þ
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
473
where Dx and Dy (unit is the size of lenslet d) are the distance of two lenslets along x direction and y direction, respectively; the definitions of other parameters could be found in Eq. (3). triðxÞ is the triangle function defined as ( jxjp1; 1 jxj; triðxÞ ¼ (16) 0; elsewhere: Consider the slope structure function for slopes in the x direction. The ratio of transverse to longitudinal covariance is used to characterize the value of b [6]: RðDÞ ¼
DX S ð0; DÞ , DX S ðD; 0Þ
(17)
where D is the separation between any two lenslets. Based on this formula, we could compare the simulation with theoretical results.
3. Analysis and discussion of the simulation results To test the feasibility of the simulator, we do not consider the optical system in front of the SHWS and the optical wave will directly incident into the lenslet array. For simplification, the plane wave passing through only one phase screen located at the front of the lenslet array is simulated by using our program. The phase screens follow the Kolmogorov power spectrum [1]: 5=3
FðkÞ ¼ 0:023r0
jkj11=3 .
The corresponding phase structure function is 5=3 j rj 2 Df ðrÞ ¼ hðfðr1Þ fðr1 þ rÞ i ¼ 6:88 . r0
(18)
(19)
The adequate frame number of phase screens is necessarily required in simulation. One hand, it should ensure to meet the statistical behavior of the turbulence and, on the other hand, it should be limited to a certain amount to reduce the computation load and improve the work efficiency. We firstly compare the statistics of the simulated phase screen with different parameters (such as sampling points and turbulence character parameter). Fig. 4 gives the statistical results of the simulation phase screen with sampling points 144 144 and D=r0 ¼ 1. We can see that the behavior of the simulation will agree with the theoretical prediction with increase of amount of phase screen frames. In this case, the error will be small enough when we use 500 frames to simulate the turbulence. The statistical variance of phase screen with a certain character parameter (D=r0 ¼ 1) and different sampling points (80 80, 144 144 and 240 240) is shown in Fig. 5, and with a certain number of sampling points (80 80) and different parameters (D=r0 ¼ 1,2 and 3) is shown in Fig. 6. From these figures, we could find that the statistical error is in the descent tendency with increase of frames of phase
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
474 7
4 3
Phase structure function (Dφ )
Solid line: theoretical results, and plus symbol: simulation 1:50, 2:100, 3:300,and 4:500,respectively. 6
2
5
1
4
3
2
1
0 0.0
0.2
0.4
0.6
0.8
1.0
r/D Fig. 4. Statistics of the simulated phase screens as function of normalized distance: the simulated data: plus symbol and the theory data: solid line for D=r0 ¼ 1, and the sampling points: 144 144.
Simulated accuracy (variance σ2 )
1:Rectangle:80∗80; 2:Circle:144∗144; 3:dot:240∗240
0.01
3 1E-3
1
1E-4
2 100
Number of simulated phase screens Fig. 5. Simulated accuracy changed with the simulated frames of the phase screens for D=r0 ¼ 1.
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
475
Simulated accuracy (variance σ2 )
10
D/r0=3
1
0.1
D/r0=2
0.01
1E-3
D/r0=1
1E-4 100
1000
Number of simulated phase screens Fig. 6. Simulated accuracy changed with the simulated number of the phase screens for the sampling number of 80 80.
screen. The statistical fluctuation for small sampling points will be larger than that of the big ones. And we also find that the bigger the value of D=r0 , the worse the simulated accuracy. All those results show that the proper sampling points and enough frames are required to achieve the expected statistical behavior. In the following, we show the simulated results for a plane wave passing through the phase screen with D=r0 ¼ 3. The interval of sampling points is 7 mm, each lenslet with focal length of 25 mm is 16 16 pixels and 15 15 lenslet array is used in simulation, the frame number of phase screen is 1000. A typical spot image obtained from SHWS simulator is shown in Fig. 7. The size of subaperture is obtained according to the method outlined in Section 2 and is 16 16 pixels (those values are same as our simulation parameters but are obtained from FFT) along transverse and longitudinal directions, respectively. Fig. 8 gives the centroid of the spot for about 500 frames in one of the subapertures along x and y directions (unit is the size of pixel), respectively. The mean slope structure function for the SHWS between any two subapertures could be achieved according to the following formula: N X M 1 X ¯X D ðDÞ ¼ ½SX ðX Þ S X ðX DÞ2 , S MN n¼1 m¼1
(20)
where M is the frame number of phase screens, and N is the number of subaperture pairs with distance of D in a single frame. Fig. 9 shows the comparison of the calculation results with theoretical prediction. We can see that the simulator works
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
476
Fig. 7. 15 15 Hartmann spot array with a D=r0 ¼ 3 turbulence screen.
1.0
1.0
Centroid difference along y axis (pixel)
Centroid difference along x axis (pixel)
Position difference of centroid (0,0) along x axis
0.5
0.0
-0.5
0.5
0.0
-0.5
-1.0
-1.0 0
(a)
Position difference of centroid (0,0) along y axis
100
200
300
400
Part number of simulated phase screens
0
500
(b)
100
200
300
400
500
Part number of simulated phase screens
Fig. 8. Centroid position variance of spot (0,0): (a) along x-axis; (b) along y-axis (240 240 and D=r0 ¼ 3).
well to characterize the turbulence statistics. It also shows that our simulation is feasible. Finally, we should note that our simulation is initial and simplified. In practice, we should simulate the wave propagation according to Eq. (2); Also, the influence of the noise resulting from CCD and aberration arising from optical system and lenslet
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
477
Transverse over longitudinal differential variance
1.0
0.9
0.8
0.7
0.6
0.5
0.4 0
1
2
3
4
5
6
7
8
9
10
11
Subaperture separation/subaperture size
Fig. 9. Ratio of transverse and longitudinal differential variances in simulated data (symbols) compared with theoretical curves in the case of b ¼ 44 12, in which the turbulence spectrum obeys Kolmogorov statistics.
have not been taken into account in our simulation; however, those effects could be easily added in our simulator.
4. Conclusion In summary, a SHWS simulator has been realized and tested by using random phase screen following Kolmogorov power spectrum. The statistics of generated phase screen with different parameters have been compared and the simulation results show that the proper sampling points and enough frames are required to achieve the expected statistical behavior. The FFT technique is used in our simulation to automatically divide the detection region and the slope structure function is calculated according to the output of the simulator. The simulation results show that the simulator works well to characterize the turbulence phase statistics. It will be a useful tool to further investigate the influences of different parameters on performance of the whole system. References [1] Roggeman MC, Welsh B. Imaging Through Turbulence. Boca Raton, FL: CRC Press; 1996. [2] Fried DL. Differential angle of arrival: theory evaluation and measurement feasibility. Radio Sci 1975;10(1):71–6. [3] Fried DL. Optical resolution through a randomly inhomogeneous medium for very long and very short exposures. J Opt Soc Am 1966;56(10):1372–9.
ARTICLE IN PRESS 478
Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
[4] Sarazin M, Roddier F. The ESO differential image motion monitor. Astron Astrophys 1990;227:294–300. [5] Dayton D, Pierson B, Spielbursch B, Gonglewski J. Atmospheric structure function measurements with a Shack–Hartmann wavefront sensor. Opt Lett 1992;17(24):1737–9. [6] Nicholls TW, Boreman GD, Dainty JC. Use of a Shack–Hartmann wavefront sensor to measure deviations from a Kolmogrov phase spectrum. Opt Lett 1995;20(24):2460–2. [7] Silbaugh E, Welsh BM, Roggermann MC. Characterization of atmospheric turbulence phase statistics using wavefront slope measurements. J Opt Soc Am A 1996;13(12):2453–60. [8] Rao CH, Jiang WH, Ling N. Measuring the power-law exponent of an atmospheric turbulence phase spectrum with a Shack–Hartmann wavefront sensor. Opt Lett 1999;24(5):1008–10. [9] Fleck Jr JA, Mcrris JR, Feit MD. Time-dependent propagation of high energy laser beams through the atmosphere. Appl Phys 1976;10:129–60. [10] Martin JM, Flatte SM. Intensity images and statistics from numerical simulation of wave propagation in 3-D random media. Appl Opt 1988;27:2111–26. [11] Coles WA, Filice JP, Frehlich RG, Yadlowsky M. Simulation of wave propagation in threedimensional random media. Appl Opt 1995;12:2089–101. [12] Lane RG, Glindemann A, Dainty JC. Simulation of a Kolmogorov phase screen. Wave Random Media 1992;2:209–24. [13] Harding CM, Johnston RA, Lane RG. Fast simulation of a Kolmogorov phase screen. Appl Opt 1999;38:2161–70. [14] Rod Frehlich. Simulation of laser propagation in a turbulence atmosphere. Appl Opt 2000;39:393–7. [15] Yan Hx, Li SS, Zhang DL, Chen S. Numerical simulation of an adaptive optics system with laser propagation in the atmosphere. Appl Opt 2000;39(18):3023–31. [16] Zhiling Jiang, Shunsheng Gong, Yang Dai. Monte-Carlo analysis of centroid detected accuracy for wavefront sensor. Opt Laser Technol, in press. [17] Zhiling Jiang, Shunsheng Gong, Yang Dai. Numerical Study of centroid detection accuracy for Shack–Hartmann Wavefront Sensor. Opt Laser Technol, in press. [18] Nicolas Roddier. Atmospheric wavefront simulation using Zernike polynomials. Opt Eng 1990;29(10):1174–80. [19] Mitsuo Takeda, Hideki Ina, Seiji Kobayashi. Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry. J Opt Soc Am 1982;72(1):156–60.
Optics and Lasers in Engineering 44 (2006) 466–478
Numerical simulation of Shack–Hartmann wavefront sensor to characterize the turbulence phase statistics Zhiling Jiang, Yang Dai, Faquan Li, Xuewu Cheng, Shunsheng Gong Wuhan Institute of Physics and Mathematics, The Chinese Academy of Sciences, P.O. Box 71010, Wuhan, Hubei, 430071, China Received 17 November 2004; received in revised form 8 March 2005; accepted 5 April 2005
Abstract Shack–Hartmann wavefront sensor (SHWS) is one of the important parts in an adaptive optics system and its detection accuracy has considerable influences on the performances of the whole system. In this paper, a simulator of such type of wavefront sensor is built and tested by measuring the slope structure function of Kolmogorov turbulence simulated by random phase screens. The numerical results show that it works well to characterize the turbulence phase statistics. In combination with the noise simulation, this simulator will be a very useful tool to further investigate and optimize the system parameters to improve the detection accuracy and the performance of the whole system. r 2005 Elsevier Ltd. All rights reserved. Keywords: Numerical simulation; Hartmann wavefront sensor; Turbulence statistics characterization; Slope structure function
Corresponding author.
E-mail address: [email protected] (Z. Jiang). 0143-8166/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlaseng.2005.04.012
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
467
1. Introduction It is well known that the wavefront distortion resulting from atmospheric turbulence is the crucial factor limiting the resolution of the optical system. The adaptive optics technology is usually utilized to actively correct the received wavefront and improve the performance of the system [1]. The wavefront shearing interferometer and Shack–Hartmann wavefront sensor (SHWS) are the most popular devices used to measure the distorted wavefront. The latter has the advantages such as short response time, easy operation, and convenient to comprehend the wavefront modes and has been widely used in real-time adaptive optical systems, beam quality diagnosis and the test of the optical components. Characterization of the turbulence statistics is very important to understand its influences on optical system. Based on the wavefront slope measurement, Fried [2] presented the theoretical expression of slope structure function for circular aperture, which could be applied to the long-time exposure imaging by using a single aperture or differential image motion measurement by using a pair of subapertures [3,4]. Later, Dayton [5] and Nicholls [6], etc. introduced and discussed the measurement of turbulence statistics by using SHWS; Silbaugh [7] derived the expression of slope structure function of arbitrary power-law turbulence power spectrum and applied them to analyze the experimental results; Rao [8] further extended that theory and showed the comparison of the theoretical results and experimental measurements. All those research work showed that the SHWS could be conveniently applied to characterize the atmospheric turbulence statistics; however, few of studies have been carried out to investigate the influences of the sensor parameters such as size of the subaperture and sampling points as well as lenslet array on such characterization. Thus, such wavefront sensor simulator is essentially required. Due to the advantages of low cost and easy realization, numerical simulation technology has been widely applied to study the problems of the atmosphere turbulence [9–15], such as analysis of the performance of the adaptive optical systems, design and optimization of the system parameters, and simulation of the experiments that could not be easily or impossibly realized in practice. Recently, we have applied this technology to numerically analyze the influences of noise on centroid detection accuracy [16,17]. In this paper, we initially set up a SHWS simulator, and test its performance by using Kolmogorov phase screens generated by random numbers. The simulated results show that this simulator could be used to characterize the turbulence statistics and will be a useful tool to further investigate the influences of system parameters on measurement results. The main contents of this paper include: in Section 2, the theoretical considerations of the numerical simulation such as generation of the phase screens, model of wavefront sensor, automatic detection of spot centroid and characterization of turbulence parameters are presented; in Section 3, the simulation results and the discussion are given and the conclusion is outlined in Section 4.
ARTICLE IN PRESS 468
Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
2. Theoretical analysis and numerical treatment Fig. 1 shows the principle of atmospheric turbulence measurement by using SHWS. The waves affected by refractive index fluctuation of atmosphere incident on the entrance pupil and enter into the lenslet array after an optical system. The wavefront is split into many subapertures and imaged on the focal plane of the lenslet. A CCD is used to measure the centroid of the spot in each subaperture. Under the parabolic axis approximation, the equation describing the wave propagation along z direction in random medium could be simplified as [9] 2 qc q c q2 c 2ik þ 2k2 Dnc ¼ 0, qz q2 x q2 y
(1)
where c is the complex amplitude of the wave, k ¼ 2pn=l0 is the wave number, l0 is the wavelength in vacuum and n ¼ n0 þ Dn is the refractive index of air, n0 and Dn are the constant and the fluctuation of the refractive index, respectively. The popular method to solve the above equation is splitting-step Fourier transform algorithm, i.e. it splits the propagation distance into many segments with interval of Dzn and in the nth segment the wave solution could be represented by "
( Fn ¼ F F
1
Fn1 exp i
k2x þ k2y 4k
!# Dzn
) expðikDnDzn Þ exp i
k2x þ k2y 4k
! Dzn ,
(2) where F is the spectrum amplitude of the wave field, kx and ky are the transverse components of the wave vector, F and F 1 are Fourier and inverse Fourier transform operations. From this expression, we could see that the first and the third items represent the wave propagation in vacuum, and a phase item is added (the second item) to show the influence of the turbulence on wave field in this distance. Hence, the propagation of wave in random medium could be described by the propagation in vacuum and a series of discrete phase screens. If the nonisoplanatic effects are not taken into account, a single phase screen located at the entrance pupil of the optical system will be adequate to approximate the static turbulence [1].
Fig. 1. Principle scheme of atmosphere turbulence measurement.
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
469
Simulation of phase screens
Statistics of the phase screens
Wave propagation and imaging
Spot centroid detection
Turbulence characterization analysis Fig. 2. Flowing chart of turbulence simulation, Hartmann spots array generation and atmospheric turbulence phase statistics measurement.
According to the scheme shown in Fig. 1, the aspects that should be considered in numerical simulation are summarized as: generation of random phase screen, describing of the lenslet array, centroid detection and analysis of the simulation results. We will discuss each of them in the following sections. Fig. 2 shows the flowing chart of simulation. 2.1. Generation of phase screen The phase screen generated from the random numbers under certain conditions could simulate the statistical character of the atmospheric turbulence. Fast Fourier transform algorithm [9] is the most basic approach to generate the phase screens. However, in order to sample adequately low frequencies, the width of the phase screens needs to be in the order of the outer scale, which may be hundreds of meters, and due to the fact that the limited sampling interval and the coherent periodic property of discrete Fourier transform, the conventional algorithm needs a large number of sampling points and much computational time. This disadvantage could be overcome by adding high-order harmonics as shown in Refs. [12,14]. The random phase screen could also be well generated from Zernike functions [18], which has the advantage that the remaining statistical variance could be estimated from this simulation and provides useful information for correction system; but the number of Zernike functions increases rapidly with the size of phase screen, which raises the difficulty for simulation. Another limitation of this algorithm is that it only simulates the turbulence in circular region due to the definition of the Zernike function. In our simulation, we use the eigenvalue method and random mid-point displacement technique to generate the phase screen with a certain power spectrum [13].
ARTICLE IN PRESS 470
Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
The general phase structure function of turbulence statistics could be defined as [6] b2 x1 x2 ; y1 y2 , (3) Dðx1 ; y1 ; x2 ; y2 Þ ¼ gb r0 where gb is the generalized coefficient, r0 is an analogous quantity to describe the turbulence and b is the power-law exponent of the phase spectrum. If b is equal to 11 3, r0 is the coherent length of Kolmogorov turbulence r0 . If fðx; yÞ represents the phase screen (assuming that the mean component in the whole aperture has been removed), the covariance of the phase screen is given by Cðx1 ; y1 ; x2 ; y2 Þ 1 ¼ Dðx1 ; y1 ; x2 ; y2 Þ 2ZZ 1 0 0 0 0 0 0 Dðx1 ; y1 ; x2 ; y2 ÞW ðx1 ; y1 Þ dx1 dy1 þ 2 ZZ 1 0 0 0 0 0 0 Dðx1 ; y1 ; x2 ; y2 ÞW ðx2 ; y2 Þ dx2 dy2 þ 2 Z Z Z Z 1 0 0 0 0 0 0 0 0 0 0 0 0 Dðx1 ; y1 ; x2 ; y2 ÞW ðx1 ; y1 ÞW ðx2 ; y2 Þ dx1 dy1 dx2 dy2 . 2
ð4Þ
In practice, it requires to rearrange the screen into a one-dimensional vector to calculate the covariance: FðmM þ nÞ ¼ T½fðmDx; nDyÞ;
m; n ¼ 1; 2; . . . M,
(5)
where M is sampling points in rows and columns, T is the transform function. The covariance matrix could be written as C ¼ hFFT i.
(6)
The eigen-system of C could be used to describe the phase screen. Assuming that the eigenvalues of C are l ¼ ðl1 ; l2 . . . lMM Þ, and the eigenvectors are E ¼ ðE 1 ; E 2 . . . E MM Þ, phase screen could be simulated by using the random vector X where each element is independent and has its variance given by the corresponding element of l: F ¼ EX .
(7)
Due to the limitation of computation load, the number of sampling points is small and the random mid-point displacement technique is used to improve the resolution of the phase screen. According to the self-similarity of turbulence, the phase screen with different sizes could be achieved by adding a modified factor [12]. 2.2. The description of lenslet array The lenslet array is an important part of SHWS, which transfers the wavefront measurement to centroid detection. According to the wave optics, the lens could be
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
treated as a phase device, and the mathematic expression is k FL ðx; yÞ ¼ wðx; yÞ exp i ðx2 þ y2 Þ , 2f
471
(8)
where k is the wave number, f is the focal length and wðx; yÞ is the aperture function defined as ( 1; jxjp d2 ; yp d2 ; wðx; yÞ ¼ (9) 0; others; where d is the size of the lens. It is easy to obtain the intensity distribution in detection plane as Z 1 Z 1 2 k Iðxf ; yf Þ / c1 ðx; yÞwðx; yÞ exp i ðxxf þ yyf Þ dxdy f 1 1 2 ¼ FTfc ðx; yÞwðx; yÞg , ð10Þ 1
where c1 is the wave field distribution at the front of the lenslet array. From Eq. (10), we can see that the intensity distribution could be obtained by using fast Fourier transform. It should be noted that the optical system in front of the lenslet array is not taken into account in our present model because it is easily solved by changing the absolute size of the wave field and no other changes are required for the sampling points and the calculation matrix. Of course, a phase modification should be added into the wave field if the influence of the aberration resulting from the optical system is investigated. 2.3. Automatic detection of the spot centroid How to accurately measure the centroid is important for SHWS, which will determine the performance of the system. Different to the simulation, the detection region for each subaperture is not known in advance in practical measurement. To simulate the practical case, we adapt the FFT technique [19] to realize the division of subapertures and the detection of centroids. For a periodic interferometric image, its intensity could be written as, for example, gðx; yÞ ¼ aðx; yÞ þ bðx; yÞ cos½2pf 0 x þ fðx; yÞ.
(11)
For simpleness, we only discuss the case along x direction. Usually, the phase information could be obtained by using FFT technique (the typical instrument is shearing interferometer). In our simulation, we use this method to determine the frequency f 0 , which represents the average interval of the fringes. Eq. (11) could be re-written as gðx; yÞ ¼ aðx; yÞ þ cðx; yÞ expði2pf 0 xÞ þ c ðx; yÞ expði2pf 0 xÞ, where cðx; yÞ ¼ 12bðx; yÞ exp½ifðx; yÞ.
(12)
ARTICLE IN PRESS 472
Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
Fig. 3. Typical simulation result using FFT technique.
The spectrum distribution could be obtained if we apply FFT to Eq. (12): Gðf x ; yÞ ¼ Aðf x ; yÞ þ Cðf x f 0 ; yÞ þ C ðf x þ f 0 ; yÞ,
(13)
where f x is the frequency along x direction. We could find that A is located at the origin and C and C are located at f ¼ f 0 and f ¼ f 0 , respectively. f 0 could be determined and the spot interval in spatial could be known according to Fourier transform relationship. In our simulation, we firstly determine the centroid position of the whole wave field in detection region, and then divide the subapertures according to the transverse and longitudinal interval. A typical simulation result is shown in Fig. 3, where the left shows the interferometric image and the right shows the spectrum distribution, which has been amplified. 2.4. Character of turbulence statistics The statistical character could be achieved by using detection information of SHWS [5–8]. The slope structure function for SHWS is defined as [7] DaS^ ðX ; X 0 Þ ¼ h½S a^ ðX Þ Sa^ ðX 0 Þ2 i,
(14)
where S is the wavefront slope, and a^ is the normal direction (x or y). After a certain mathematical calculation, the final result is DaS^ ðDx; DyÞ ¼ gb d 2
b2 Z d du triðuÞf2½Dx2 þ ðu þ DyÞ2 ðb2Þ=2 r0
½ðDx 1Þ2 þ ðu þ DyÞ2 ðb2Þ=2 ½ðDx þ 1Þ2 þ ðu þ DyÞ2 ðb2Þ=2 2½jujb2 ð1 þ u2 Þðb2Þ=2 g,
ð15Þ
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
473
where Dx and Dy (unit is the size of lenslet d) are the distance of two lenslets along x direction and y direction, respectively; the definitions of other parameters could be found in Eq. (3). triðxÞ is the triangle function defined as ( jxjp1; 1 jxj; triðxÞ ¼ (16) 0; elsewhere: Consider the slope structure function for slopes in the x direction. The ratio of transverse to longitudinal covariance is used to characterize the value of b [6]: RðDÞ ¼
DX S ð0; DÞ , DX S ðD; 0Þ
(17)
where D is the separation between any two lenslets. Based on this formula, we could compare the simulation with theoretical results.
3. Analysis and discussion of the simulation results To test the feasibility of the simulator, we do not consider the optical system in front of the SHWS and the optical wave will directly incident into the lenslet array. For simplification, the plane wave passing through only one phase screen located at the front of the lenslet array is simulated by using our program. The phase screens follow the Kolmogorov power spectrum [1]: 5=3
FðkÞ ¼ 0:023r0
jkj11=3 .
The corresponding phase structure function is 5=3 j rj 2 Df ðrÞ ¼ hðfðr1Þ fðr1 þ rÞ i ¼ 6:88 . r0
(18)
(19)
The adequate frame number of phase screens is necessarily required in simulation. One hand, it should ensure to meet the statistical behavior of the turbulence and, on the other hand, it should be limited to a certain amount to reduce the computation load and improve the work efficiency. We firstly compare the statistics of the simulated phase screen with different parameters (such as sampling points and turbulence character parameter). Fig. 4 gives the statistical results of the simulation phase screen with sampling points 144 144 and D=r0 ¼ 1. We can see that the behavior of the simulation will agree with the theoretical prediction with increase of amount of phase screen frames. In this case, the error will be small enough when we use 500 frames to simulate the turbulence. The statistical variance of phase screen with a certain character parameter (D=r0 ¼ 1) and different sampling points (80 80, 144 144 and 240 240) is shown in Fig. 5, and with a certain number of sampling points (80 80) and different parameters (D=r0 ¼ 1,2 and 3) is shown in Fig. 6. From these figures, we could find that the statistical error is in the descent tendency with increase of frames of phase
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
474 7
4 3
Phase structure function (Dφ )
Solid line: theoretical results, and plus symbol: simulation 1:50, 2:100, 3:300,and 4:500,respectively. 6
2
5
1
4
3
2
1
0 0.0
0.2
0.4
0.6
0.8
1.0
r/D Fig. 4. Statistics of the simulated phase screens as function of normalized distance: the simulated data: plus symbol and the theory data: solid line for D=r0 ¼ 1, and the sampling points: 144 144.
Simulated accuracy (variance σ2 )
1:Rectangle:80∗80; 2:Circle:144∗144; 3:dot:240∗240
0.01
3 1E-3
1
1E-4
2 100
Number of simulated phase screens Fig. 5. Simulated accuracy changed with the simulated frames of the phase screens for D=r0 ¼ 1.
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
475
Simulated accuracy (variance σ2 )
10
D/r0=3
1
0.1
D/r0=2
0.01
1E-3
D/r0=1
1E-4 100
1000
Number of simulated phase screens Fig. 6. Simulated accuracy changed with the simulated number of the phase screens for the sampling number of 80 80.
screen. The statistical fluctuation for small sampling points will be larger than that of the big ones. And we also find that the bigger the value of D=r0 , the worse the simulated accuracy. All those results show that the proper sampling points and enough frames are required to achieve the expected statistical behavior. In the following, we show the simulated results for a plane wave passing through the phase screen with D=r0 ¼ 3. The interval of sampling points is 7 mm, each lenslet with focal length of 25 mm is 16 16 pixels and 15 15 lenslet array is used in simulation, the frame number of phase screen is 1000. A typical spot image obtained from SHWS simulator is shown in Fig. 7. The size of subaperture is obtained according to the method outlined in Section 2 and is 16 16 pixels (those values are same as our simulation parameters but are obtained from FFT) along transverse and longitudinal directions, respectively. Fig. 8 gives the centroid of the spot for about 500 frames in one of the subapertures along x and y directions (unit is the size of pixel), respectively. The mean slope structure function for the SHWS between any two subapertures could be achieved according to the following formula: N X M 1 X ¯X D ðDÞ ¼ ½SX ðX Þ S X ðX DÞ2 , S MN n¼1 m¼1
(20)
where M is the frame number of phase screens, and N is the number of subaperture pairs with distance of D in a single frame. Fig. 9 shows the comparison of the calculation results with theoretical prediction. We can see that the simulator works
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
476
Fig. 7. 15 15 Hartmann spot array with a D=r0 ¼ 3 turbulence screen.
1.0
1.0
Centroid difference along y axis (pixel)
Centroid difference along x axis (pixel)
Position difference of centroid (0,0) along x axis
0.5
0.0
-0.5
0.5
0.0
-0.5
-1.0
-1.0 0
(a)
Position difference of centroid (0,0) along y axis
100
200
300
400
Part number of simulated phase screens
0
500
(b)
100
200
300
400
500
Part number of simulated phase screens
Fig. 8. Centroid position variance of spot (0,0): (a) along x-axis; (b) along y-axis (240 240 and D=r0 ¼ 3).
well to characterize the turbulence statistics. It also shows that our simulation is feasible. Finally, we should note that our simulation is initial and simplified. In practice, we should simulate the wave propagation according to Eq. (2); Also, the influence of the noise resulting from CCD and aberration arising from optical system and lenslet
ARTICLE IN PRESS Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
477
Transverse over longitudinal differential variance
1.0
0.9
0.8
0.7
0.6
0.5
0.4 0
1
2
3
4
5
6
7
8
9
10
11
Subaperture separation/subaperture size
Fig. 9. Ratio of transverse and longitudinal differential variances in simulated data (symbols) compared with theoretical curves in the case of b ¼ 44 12, in which the turbulence spectrum obeys Kolmogorov statistics.
have not been taken into account in our simulation; however, those effects could be easily added in our simulator.
4. Conclusion In summary, a SHWS simulator has been realized and tested by using random phase screen following Kolmogorov power spectrum. The statistics of generated phase screen with different parameters have been compared and the simulation results show that the proper sampling points and enough frames are required to achieve the expected statistical behavior. The FFT technique is used in our simulation to automatically divide the detection region and the slope structure function is calculated according to the output of the simulator. The simulation results show that the simulator works well to characterize the turbulence phase statistics. It will be a useful tool to further investigate the influences of different parameters on performance of the whole system. References [1] Roggeman MC, Welsh B. Imaging Through Turbulence. Boca Raton, FL: CRC Press; 1996. [2] Fried DL. Differential angle of arrival: theory evaluation and measurement feasibility. Radio Sci 1975;10(1):71–6. [3] Fried DL. Optical resolution through a randomly inhomogeneous medium for very long and very short exposures. J Opt Soc Am 1966;56(10):1372–9.
ARTICLE IN PRESS 478
Z. Jiang et al. / Optics and Lasers in Engineering 44 (2006) 466–478
[4] Sarazin M, Roddier F. The ESO differential image motion monitor. Astron Astrophys 1990;227:294–300. [5] Dayton D, Pierson B, Spielbursch B, Gonglewski J. Atmospheric structure function measurements with a Shack–Hartmann wavefront sensor. Opt Lett 1992;17(24):1737–9. [6] Nicholls TW, Boreman GD, Dainty JC. Use of a Shack–Hartmann wavefront sensor to measure deviations from a Kolmogrov phase spectrum. Opt Lett 1995;20(24):2460–2. [7] Silbaugh E, Welsh BM, Roggermann MC. Characterization of atmospheric turbulence phase statistics using wavefront slope measurements. J Opt Soc Am A 1996;13(12):2453–60. [8] Rao CH, Jiang WH, Ling N. Measuring the power-law exponent of an atmospheric turbulence phase spectrum with a Shack–Hartmann wavefront sensor. Opt Lett 1999;24(5):1008–10. [9] Fleck Jr JA, Mcrris JR, Feit MD. Time-dependent propagation of high energy laser beams through the atmosphere. Appl Phys 1976;10:129–60. [10] Martin JM, Flatte SM. Intensity images and statistics from numerical simulation of wave propagation in 3-D random media. Appl Opt 1988;27:2111–26. [11] Coles WA, Filice JP, Frehlich RG, Yadlowsky M. Simulation of wave propagation in threedimensional random media. Appl Opt 1995;12:2089–101. [12] Lane RG, Glindemann A, Dainty JC. Simulation of a Kolmogorov phase screen. Wave Random Media 1992;2:209–24. [13] Harding CM, Johnston RA, Lane RG. Fast simulation of a Kolmogorov phase screen. Appl Opt 1999;38:2161–70. [14] Rod Frehlich. Simulation of laser propagation in a turbulence atmosphere. Appl Opt 2000;39:393–7. [15] Yan Hx, Li SS, Zhang DL, Chen S. Numerical simulation of an adaptive optics system with laser propagation in the atmosphere. Appl Opt 2000;39(18):3023–31. [16] Zhiling Jiang, Shunsheng Gong, Yang Dai. Monte-Carlo analysis of centroid detected accuracy for wavefront sensor. Opt Laser Technol, in press. [17] Zhiling Jiang, Shunsheng Gong, Yang Dai. Numerical Study of centroid detection accuracy for Shack–Hartmann Wavefront Sensor. Opt Laser Technol, in press. [18] Nicolas Roddier. Atmospheric wavefront simulation using Zernike polynomials. Opt Eng 1990;29(10):1174–80. [19] Mitsuo Takeda, Hideki Ina, Seiji Kobayashi. Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry. J Opt Soc Am 1982;72(1):156–60.