Volume 42, number 4
OPTICS COMMUNICATIONS
15 July 1982
A PROCEDURE TO CORRECT THE IMAGES OF A S ~ A L O ~ FOR THE DISTORTIONS DUE TO ATMOSPHERIC T U ~ C E A.M.J. HUISER Ecole Polyteehnique F~d~ralede Lausanne, lnsfftut de Physique Th~orique, PHB Ecublens, CH-1015Lausanne, Switzerland -
Received 23 February 1982 Revised manuscript received27 April 1982
A method is proposed to obtain diffraction limited pictures of extended astonomieal ohmcotsin Spite of the wavefront distortions due to atmospheric turbulence, The method does not require any prior knoWle~'on ~e brightness dlstribution across the object. The presence of reference details is not necessary.
1. Introduction Images of astronomical objects obtained with earthbound telescopes are usually severely degraded due to the refractive index fluctuations in the earth's atmosphere caused e.g. by turbulence. Short time exposures, during which the atmosphere may be considered as "frozen", will yield images which are severely distorted but which still contain spacial frequencies up to the diffraction limit. Long time exposures will give images with a resolution that does not exceed that of a 10 cm telescope or worse. In this paper a method shall be presented which enables one to correct short time images for'the defects caused by the inhomogeneity of the atmosphere. In the past several techniques have been invented to reconstruct a "true", i.e. undistorted, diffraction limited picture of the object from these randomly distorted short time images. A very well known example is Labeyrie's speckle interferometry [1]. Labeyries' method exploits the randomness of the phase shifts introduced by the atmosphere to obtain a diffraction limited autoeorrelation function of the brightness distribution at the sky by summing the power spectra, i.e. the modulus of the Fourier transform of a great number of short time exposures. For simple objects this diffraction limited autocorrelation function contains enough information to retrieve a "true" image of the object. More complicated objects, however, may in practice not be uniquely defined by their autoeorrelation functions [2], although some nice results have been achieved [3]. To circumvent this problem, other techniques have been proposed by Knox and Thompson [4], Walker [5], Worden et al. [6], Bates et al. [7] and Weigelt [8]. Many of these methods require striking features such as bright or dark points in the objects to obtain an estimate of the intensity point spread function of the system atmosphere telescope during a short time exposure. An estimate of the true image is obtained by deeonvolution. The method proposed in this communication does not rely upon the presence of striking features in the object, nor does it require simple objects. It will be shown that with a suitable arrangement of the imaging system, the distorted images of completely unknown objects allow the determination of the phase shifts introduced by the atmosphere. The knowledge of these phase shifts allows one to calculate a "true" image of the object.
226
0 030-4018/82/0000-0000/$ 02.75 © 1982 North-Holland
Volume 42, number 4
OPTICS COMMUNICATIONS
15 July 1982
2. The imaging system Denoting with B(xo) the brightness distribution at the sky, one can write for the intensity distribution ~(xi) in the image plane [6] :
(xi) = f
d2xo IK/(xi - x0)] 2 B(XO),
(1)
where x i = (xi,Yi) and x 0 = (x O,YO) are the angular coordinates expressed in radians in the image plane and at the sky respectively. The function IK/(xi -x0)l 2 is the intensity point spread function of the system telescope atmosphere. Ig/(x)l 2 is the "speckle image" of a point source at the sky. While writing Ki as a function ofx i - x 0 only, it has been assumed that the image formation is isoplanatic and quasi monochromatic. One can write:
Kj( ) = f
d2u PI(u) T(u) exp(21rix • u),
(2)
where u = (u, o) are the coordinates in the plane of the entrance pupil of the telescope measured in units ?~, the mean wavelength of the recorded radiation, Pj(u) is the pupil function of the telescope and represents its aberrations; outside the aperture of the instrument P/(u) is zero. T(u) describes the phase (and amplitude) changes in the wave entering the telescope due to the random fluctuations of the refractive index of the atmosphere. As usual it is assumed that the time in which an image is recorded is sufficiently small to consider the atmosphere as frozen, i.e. T(u) as constant, in this period.
3. The method Since it will not be assumed that the object has some detail which makes it possible to obtain an estimate for IK/(x)}2, eq. (1) contains two unknowns, namely IK/(x)l 2 and B(XO). In general, one cannot expect to be able to infer both IKf(x)l 2 and B(xo) from a single image. One needs additional information. This information can be obtained by recording a second image of the object under the same atmospheric conditions 0.e. the same T(u)), but with a different pupil function. A simple arrangement with which this second image can be recorded is drawn in fig. i, and will be discussed in more detail in the next section. Fourier transforming eq. (1) and using eq. (2) one finds:
F/(u) = r(u)
f d2u' t~(u ') T(u')PT(u'- u) T*(u'- u),
/ = 1, 2,
(3)
where FI(U ) and F2(u ) are the Fourier transforms of the images obtained with the pupil functions Pl(U) and P2(u) respectively, r(u) is the Fourier transform of B(xo). One can easily eliminate F(u) by dividing eq. (3) f o r / = 1 by the same equation f o r / = 2. This leads to:
R(u) d-----efFl(u)/F2(u) = GI(U)/G2(u),
(4)
G/(u) = f d2u ' Pl(u') T(u') Pj* (u' - u)T* (u' - u)
(5)
where
is the Fourier transform of]K].(x)l 2. From (4) and (5) it is apparent that apart from T(u) also any function S(u) C exp(ia • u), with a = (a, b) a real vector and C an arbitrary complex number, is a solution of (4). So T(u) is not uniquely determined by (4). However, ifPl(U ) and P2(u) are suitably chosen, this is the only possible ambiguity. By replacing in eq. (3) T(u) by S(u), it can be seen that the only consequence of this ambiguity is that one will not Fred I'(u) itself but 227
Volume 42, number 4
OPTICS COMMUNICATIONS
15 July 1982
/ ,ELE0,
\
JECTIVE 0
I
\
I
~
IMAGES
/
/
k
\/
I
LENS~ ~o~RANSPARTEN
00:'C FL ITERS IN PIIANESOF
\
Fig, 1. A simple arrangement to obtain the two images necessaryfor the reconstsuction. ICI2 exp(ia • =)P(u). This means that instead of B(xo) one will reconstruct ICI2B(xo + a/27r), a shifted version of B. A very convenient choice of P1 (u) and P2(u) is the following @"
Pl(U,o)=P2(u,o)=l
for0~
(6a)
Pl(U, o) =A(u, v) + A(u + h, o) + A(u, - v ) + A ( u , h - o)
(6b) elsewhere,
P2(u, v) = A(u, o) - a(u + h, o) + A(u, - v ) - A(u, h - v)
(6c)
with A(u, o) = 1 i f - h - d ~ u ~ - d and - h - d ~ o ~ - d , and A(u, v) = 0 otherwise. Fig. 2 gives a visual impression of the pupil functions. The resolution that can be obtained depends on d; the strength of the turbulence fixes h (see below), d must be greater than 3h. Inserting the ddfinitions (6a, b, c) o f P 1 andP 2 into (5) one obtains for u ~ d and o ~ 0: u-d
G/(u,v)=f
du'f
v-d
dv'T(u',v')T*(u'-u,v'-v) (A(u'-u,u'-v)+_A(u'-u+h, vi-o)}, u - d - 2h v-d-h (7) where the + sign applies i f / = 1 and the - sign if] = 2. Now suppose, that for Irl ~ h/2, T(u - I') can be approximated by T(u). Such an approximation can always be justified provided that h is chosen small enough. Typically one would choose h "" 0.5 r 0 with r 0 the turbulence coherence length (see [10]). In this case, the integral on the right hand side of (7) is well approximated by a two term Riemann sum. Thus, writing: "G/(n, m ) = G/(nh , mh ) and T(n,m ) = T((n - ~)h, (m - ~)h ), eq. (7) gives for n > / l = d / h andre f> 0: "G/(n,m) = T ( n - l,m - l) "T*(-l,-l) +'T(n - l - l , n - l ) ' T * ( - I - 1 , - / ) .
(8)
Using (4), the equations (8) lead to the following recursion relation for T(n, m): T(n - l, m - l) = exp(i~) R(n, rn! 7" 1 T(n - l - 1, m - l ) ,
R(n,m)+ 1 The isoplanatic aberrations of the objective may be absorbed in T(u). 228
(9)
Volume 42, number 4
j (u,v)P
OPTICS COMMUNICATIONS
.
15 July 1982
/
+ forj= 2
/
jP(u,v)=l ,t
-d-h -d+h
////////,
~.
Fig. 2. A pair of pupil functions, or filters to be placed in the images of the aperture, which permit the determination of the wavefront distortions caused by atmosphere or otherwise.
with R(n, rn) = R(nh, mh) and exp(i~) = T*(-I - 1, -I)]T*(-I, -l). Carrying out the recursion one obtains: m)
exp(i(n
1)q0 I1 R ( t + I, m + l ) - 1 T ( 1 , m ) . t=2 R(t+ l,m +l)+ l
(10)
Since by assumption T(u) varies iitfle over distances of the order of h, T(u) can safely be approximated by smooth. ly interpolating the computed values of T((n - {)h, (m - ½)h). So, if H(u, o) smoothly interpolates the sampled values:
H((n - ½)h, (m - ½)h)= fi .~t + l,m + l ) - 1
t=2 R(t +l,m +1)+ 1
(11)
and c(o) = T(~h, o) exp(-ia/2) one has: r(u, v) = c(v) exp(iau/h ) H(u, v) .
(12)
So T(u, v) has now been determined up to an unknown function c(v) and a linear exponential. To determine c(o), consider eq. (5) for u ~ d and o < 0. In that case one finds similar to (7): u-d
o+d+2h
Gi(u,o)=/_d_h~' fa
do'T(u',o')T*(u'-u,d-@ {A(u'-u,o'-o).+-.A(u'-u,o-o'+h)}. (13)
Copying the reasoning leading from (7) to (12) one finds:
T(u, o) = d(u) exp(i[3v/h) L(u, o),
(14)
where exp(iff) = T*(-I, l + 1)/T*(-/, 1 + 2), the unknown function d(u) = T(u,d - ~h) e x p ( - i ( / - {)O) and
L(u, v) interpolates: I-1
,.,
L((n-½)h,(m-½)h)= H R(n+l,t-lt=m R ( n + l , t - l -
1)+1 1)- 1
(lS)
Equating the right hand sides of the eqs. (12) and (15) one can compute Co(O) = exp(i~/h) c(o). Inserting this into eq. (12) one gets: T(u, v) = exp(i(au + ~o)[h) Co(V) H(u, o).
(16)
Apart from the fact that the constants a and fl can be complex, the remaining unknown linear exponentials in (16) are exactly the terms expressing the ambiguity mentioned previously. Since T(u) usually varies only little 229
Volume 42, number 4
OPTICS COMMUNICATIONS
15 July 1982
across the aperture, the imaginary parts of a and 13may be accurately determined by minimizing the variation of IT(u)I across the aperture. Hence, T(u) and therefore also G/(u), have now been determined up to the anticipated inevitable ambiguity mentioned at the beginning of this section, as was to be shown.
4. Discussion The setup shown in fig. 1 is a simple arrangement with which the two required images can be recorded. It operates as follows: Using a secondary lens and a semi-transparent mirror two identical copies of the wave entering the telescope are obtained. By placing fdters, as described in the preceding section and illustrated in fig. 2, in the planes of the images of the telescope aperture, the two different pupil functions are simulated. The images of the object formed after the f'dters are those needed for the reconstruction. The arrangement described above is not the only one to which the method proposed here applies. One can use simpler pupil functions. For instance, those which obtain by omitting the terms A(u, o) and A(u,-o) in the deftnition of P1 (u, o) and removing the terms A(u + h, o) and A(u, - o ) from P2(u, o). In that case P/(u) is either 0 or 1, which means that the corresponding filters are just opaque screens with holes. The resulting calculations are almost identical to the ones given here. The cost of this simplification, however, is increased noise. Other modifications may include a change in shape of the greater transparent part in the f'dter for u i> 0 into e.g. a circle. This can be desirable when one uses separate small mirrors (or lenses) to obtain the reference waves (which are the waves emerging from the smaller transparent parts in the filters). Also many other arrangements, which shall not be mentioned here can be used. Preliminary results of computer simulations indicate that fairly high noise levels on F/(u) can be tolerated before the results of the reconstruction become meaningless. This is due to the fact that the recursion formulae for the various lines are independent. Hence, the recursion procedure is quase one-dimentional and each recursion chain in the u-direction contains only d/h steps. With d = l m (which is equivalent to a 2m telescope) and r 0 = 0.1 m, one may go up to 30% noise before the calculated phases of the endpoints of the recursion chains become meaningless. In the o-direction the chains are longer but in that case only the relative phases of lines have to be calculated. Hence, one can average over lines, which more than compensating the increased length. As usually is the case with speckle images, the signal is only a few percent of the total intensity. So, the required number of photons per image is rather high, even if noise levels up to 30% can be tolerated.
References [1] A. Labeyrie, Astron. and Astrophys. 6 (1970) 85. [2] A.M.J. Huiser and P. van Toorn, Optics Lett. 5 (1980) 499. P. van Toorn, A.M.J. Huiser and A.H. Greenaway, Solar instrumentation what's next, ed. R.B. Dunn, Sacramento Peak National Observatory (1981) p. 510. [3] J.R. Fienup, Optics Lett. 3 (1978) 27. [4] K.T. Knox and B.J. Thompson, Astrophys. J. 193 (1974) IA5. [5] J.G. Walker, Optica Acta 28 (1981) 1017. [6] S.P. Worden, M.K.Stein, G.D. Schmidt and J.R.P. Angel, Icarus 32 (1977) 450. [7] R.H.T. Bates and R.M. Lewitt, Optik 44 (1975) 1. [8] G. Weigelt, Optics Comm. 21 (1977) 55. [9] M. Born and E. Wolf, Principles of optics (Pergamon Press, Oxford, 1975) ell. 9. [10] D.L. Fried, J. Opt. Soc. Am. 55 (1965) 1427.
230