1. Introduction

The distorting effects of atmospheric turbulence on ground-based imagery of astronomical and space objects are well-documented, and numerous methods have been proposed, implemented, and tested for their correction¹. These methods can be broadly classified into two categories: post-processing through signal processing whereby the image distortions are identified and compensated by computer processing of the recorded imagery; and pre-processing through adaptive optics whereby the distortions are identified and compensated with wavefront sensors and deformable mirrors so that undistorted imagery is, in principle, recorded directly. With current technology the cost of adaptive compensation can be prohibitive, particularly if full compensation is desired. Because of this, many systems in use today benefit only from partial compensation or none at all, and post-processing methods such as the one discussed in this paper provide an economical means for restoring a telescope’s lost resolving power.

It is widely accepted that only modest improvements in image quality can be attained through the measurement and processing of long exposure (on the order of minutes to hours) telescope imagery. This difficulty is often quantified through Fried’s seeing parameter r ₀, which provides a measure of the largest spatial frequency that is present in the long-exposure imagery. Whereas a diffraction-limited telescope records spatial frequencies as high as the ratio D/λd_o , where D is the telescope diameter, λ is the nominal wavelength being sensed, and d_o is the distance from the telescope to the object, long-exposure imagery acquired through atmospheric turbulence will only contain frequencies to the turbulence-limited cut-off frequency r ₀/λd ₀, where r ₀ can take values under 10cm. For a 2m telescope with r ₀ = 10cm, then, recovery of diffraction-limited resolution requires one to super-resolve the long-exposure imagery by a factor of 20. Such improvements are rarely, if ever, accomplished with real data.

Exposure times on the order of a few milliseconds, however, provide imagery with spatial frequency information out to the diffraction limit. This fact was first exploited by Labeyrie², and has since spawned a number of signal-processing methods broadly classified as speckle imaging techniques. A review of these is found in Ref. 1 and the many references within. Many of these, such as Labeyrie’s speckle interferometry, the Knox-Thompson method³, and the triple-correlation method⁴, rely on the existence of a data transformation that is invariant to the effects of atmospheric turbulence. After performing this transformation along with any required calibration, a non-trivial image-recovery problem must be solved to obtain an estimate of the underlying object. As an alternative to these methods, the data can be processed to directly estimate and correct for the effects of atmospheric turbulence. Methods based on this approach are broadly classified as multiframe blind deconvolution ⁵ techniques. The application of a multiframe blind deconvolution method to real telescope data is the focus of this paper.

2. Mathematical model

For vertical path, ground-to-space imaging, the turbulence-induced distortions are commonly modeled as a time-varying phase screen that induces a point-spread function given by:

where: u and y are two-dimensional spatial variables in the telescope-pupil and image planes, respectively; A(·) is the telescope pupil function; θ_t (·) is the time-varying phase distortion caused by turbulence-induced changes in the refractive index of the atmosphere; λ is the nominal wavelength; and d_o is the distance from the object to the telescope. This is the isoplanatic approximation and results in an imaging model that is linear and spatially invariant. In the absence of other effects, then, the time-varying intensity recorded by the telescope’s camera is:

where o(·) is the object’s incoherent scattering function, after magnification and inversion⁶. The notation i(·; θ_t , o) shows the image intensity’s explicit dependence on the unknown object o(·) and turbulence-induced phase aberrations θ_t (·).

Many telescope imaging systems in use today use charge-coupled-device (CCD) cameras to record the image intensity. As discussed by Snyder et al. in Ref. 7, however, the data acquired by these cameras include other effects such as internal and external background counts, nonuniform camera response, electronic read-out noise, and other sources of camera bias. Accordingly, a reasonable model for the data acquired by the nth detector element in a CCD camera array during the kth exposure is:

where N_k [n] is a Poisson distributed random variable that models the number of object-induced photocounts recorded by the detector, M_k [n] is a Poisson distributed random variable that models the number of background-induced (both internal and external) photocounts recorded by the detector, g_k [n] is a Gaussian distributed random variable that models the CCD read-out noise, and b[n] is a deterministic bias induced by the electronic circuitry used to acquire and record the data. All of these variables are usually independent from each other and across detector elements.

The expected number of object-induced photocounts is modeled as:

where t_k is the time and T the duration of the kth exposure, y_n is the spatial region over which the nth detector element integrates the image intensity, and γ[n] is a spatially-varying scale factor that accounts for detector efficiency. Pixels that record meaningless data are accommodated with this model by setting γ[n] = 0. The phase aberrations are usually constant over a short exposure, and the detector regions are often small in size relative to the fluctuations in the image intensity so that the integrating operation can be modeled by the sampling operation:

where y_n is the location of the nth detector element and |y| is its integration area. Furthermore, the object and point-spread functions are often approximated on a discrete grid so that:

and

where Δ_y, Δ_x, and Δ_u are the image-, object-, and pupil-plane grid sizes, respectively. When the pupil-plane array size is N × N and Δ_uΔ_y/λd_o = 1/N, the computations required to create the point-spread function can be accomplished efficiently by using a fast Fourier transform (FFT) algorithm.

The expected number of background-induced photocounts M_k [n] is denoted by I_b [n] and assumed to be the same for each exposure. The CCD read-out noise g_k [n] are typically modeled as zero-mean, independent random variables with variance σ ²[n]. Appropriate values for the gain function γ[·], background mean I_b [·], read-noise variance σ ²[·], and bias b[·] are usually determined through a carefully controlled study of the data acquisition system. This typically involves the analysis of data acquired both without illumination (dark frames) and with uniform illumination (flat fields).

3. Maximum-likelihood estimation

Consistent with the noise models developed in the previous section, the data recorded by each detector element in a CCD camera are a mixture of independent Poisson and Gaussian random variables. Accordingly, the probability of receiving N object- and background-induced photocounts in the nth detector element during the kth exposure is:

where

Here, o[m] = T|y|${\mathrm{\Delta }}_x^2$ o(x_m ) contains the dependence on the unknown object, and h(·; θ_tk) contains the dependence on the unknown phase aberrations. Furthermore, the probability density for the read-out noise is

and the density for the measured data at each detector element is:

For a given set of data {d_k [·]}, the maximum-likelihood estimate of the unknown object and phase aberrations is selected to maximize the likelihood:

or, as is commonly done, its logarithm (the log-likelihood):

The complicated form for the measurement density (Eq. 11) makes the maximization of 𝐿(o, θ) an overly complicated problem. When the read-out noise variance is large (greater than 50 or so), however, the data can be modified according to:

and these modified data are then similar (in distribution) to the sum of three Poisson-distributed random variables⁸. Accordingly, the modified data are approximately Poisson distributed with the mean value I_k [n] + I_b [n] + σ ²[n]. The probability mass function for these data is then modeled as:

so that the log-likelihood is:

where terms not dependent on the object or phase aberrations have been omitted. The maximum-likelihood phase aberration and object estimates are then those that maximize 𝐿(o, θ) subject to any physical constraints that are imposed on these parameters. For all problems we have considered, the object intensity o[·] is constrained to be a nonnegative function. Known object support is another constraint that is commonly used; however, implementation of this constraint in an optimization procedure usually results in a trivial reduction of the number of unknown parameters.

3.1 Comments

In recent years, many methods for the blind deconvolution of astronomical and space-object imagery have been proposed. These can be broadly classified into two categories: i) the Ayers-Dainty⁹ class of iterative algorithms; and ii) constrained numerical optimization algorithms. The methodology described here is one of constrained numerical optimization, and is a straightforward extension of the method proposed by Schulz⁵ with important modifications to accommodate real sensor artifacts such as read-noise, internal and external background, and non-uniform gain. Other methods of constrained optimization have been proposed^10,11,12, and some have relied on essentially the same methodology as the one proposed in Ref. 5 with the differences being the method used for numerical optimization¹², and the utilization (or lack thereof) of realistic sensor models. Key points regarding the methodology presented in this paper are:

3.2 Numerical optimization

Many methods are available for performing the optimization required to solve for the maximum-likelihood aberration and object parameters. Gradient-based methods and the expectation-maximization (EM) method proposed in Ref. 5 are two popular choices. However, one should not confuse the method used for solving the optimization problem with the methodology used to pose the problem. Provided the problem is successfully solved, the optimization method will only influence algorithm performance (computation time). The methodology used to define the problem, which is greatly dependent upon the use of sound mathematical and physical models, will determine estimator performance (bias and variance).

With straightforward modifications to accommodate sensor effects such as non-uniform gain and backgrounds, the EM method of Ref. 5 can be applied to this problem. Additionally, the derivatives of the log-likelihood with respect to the object and phase-aberration parameters are easily derived and can be used in a gradient-based optimization method. Both methods have performed well for this problem, and the determination of an optimal method of optimization is still an open problem.

4. Results with real data

Telescope data were acquired at the Air Force Maui Optical Station (AMOS) on March 9, 1995 and subsequently processed to determine maximum-likelihood object estimates. The telescope and imaging-system parameters are listed in Table 1. The object under observation was the Hubble Space Telescope in its 600 km orbit, and 656 short exposure images were acquired, each with an 8 ms exposure time. Four of these images are shown in Fig. 1. The 656 images were partitioned into 41 sets of 16, and each set in the partition was used to obtain one object restoration. The optimization method used was a straightforward extension of the expectation-maximization method described in Ref. 5. Four of these restorations are shown in Fig. 2. Each restoration required approximately 15 minutes of processing time on 9 nodes of an IBM SP2 at the Maui High Performance Computing Center. Whereas very little detail is available in the raw data, the restored images allow one to clearly identify the satellite’s solar panels along with the cover (lower left) for the space telescope’s 2.4 m primary mirror. An mpeg video showing all data and restorations is available at the URL listed as Ref. 14.

5. Summary

For ground-based surveillance of astronomical and space objects, post-detection processing through multiframe blind deconvolution is a viable method for restoring image resolution. To support this claim, a maximum-likelihood estimation method has been implemented and tested on real telescope data, and preliminary results of this study are presented in this paper. Much work still needs to be done, however, and important topics for consideration include: the use of prior models for the object and phase aberration parameters; the optimal use of measurement diversity as supplied by a wavefront sensor or phase-diversity system¹⁵; clever methods for generating initial parameter estimates; and the development and implementation of improved methods of optimization.

Table 1:. Telescope and imaging system parameters for the AMOS imagery of the Hubble Space Telescope.

View Table

 

Figure 1: Four short-exposure images of the Hubble Space Telescope as acquired by the 1.6 m AMOS telescope.

Download Full Size | PDF

 

Figure 2: Four restored images of the Hubble Space Telescope.

Download Full Size | PDF