A phase model using the Huber norm for estimating point spread function under frozen flow hypothesis
Introduction
The observations of astronomical objects, e.g., nebula, galaxy and star cluster, from a ground-based telescope are generally interfered by the Earth’s atmosphere, i.e., the images of target objects are commonly degraded by the fickle and irregular atmospheric turbulence. To obtain the “ideal” astronomical images, astronomers are interested in quantifying the effects of atmospheric turbulence on an astronomical imaging system. Hence, both modeling and algorithmic studies on recovering the observed astronomical images are central topics in ground-based astronomy.
Typically, before the light of an astronomical object reaches the aperture of a ground-based telescope, its wavefronts are inevitably distorted by the atmospheric turbulence. The imaging system of a ground-based telescope can be modeled as where is the image of a target object in the outer space and is its observation through a ground-based telescope, is the spatially-invariant point spread function (PSF, also called the impulse response function referring to the response of an imaging system to the Dirac impulse) caused by atmospheric turbulence, “*” is the convolutional operation, and is the zero-mean white noise. The model (1) has been extensively studied in the fields, e.g., partial differential equations [1], numerical algebra [2], statistics and multi-scale functional analysis [3]. When the PSF is known, the model (1) is a trivial linear inverse problem. However, the PSF is generally unknown in practical applications, see e.g., [4], [5], [6], [7] and references therein. Thus, estimating PSF is an essential step for further processing an astronomical image. Our discussion on the model (1) is under the assumption that the PSF is unknown.
In ground-based astronomy, the PSF describes the impulse response of a long-distance imaging system to a point source. It is posed by the distortion effects of atmospheric turbulence. To tackle the ground-based astronomical images via (1), the first step is to estimate the PSF . A popular method for this purpose is attributed to the Fourier optics approach (see e.g., [8]). The basis of the Fourier optics approach is that the PSF can be expressed as the wavefront aberration of incoming light at the telescope aperture. Specifically, where denotes the Fourier transform and is its inverse, is the pupil function of telescope (1 inside the pupil and 0 otherwise), , and is the phase measuring the deviation of wavefront from planarity. Typically, the intensity of atmospheric turbulence can be quantified by the seeing condition, which is measured by the ratio with the diameter of the telescope and the Fried parameter (see e.g., [9]).
To calculate the PSF in (2), the phase should be firstly derived. One method for deriving the phase can be narrowed down to the estimate of its gradients, which are collected by the wavefront sensor (WFS) equipped on the telescope aperture. Let and denote the phase gradients in horizontal and vertical directions, respectively. Then, the phase , and its gradients and are involved in the following overdetermined linear system (see e.g., [10], [11], [12], [13]) where , are the derivative operators (see e.g., [14], [15]). Apparently, the high accurate estimates of the phase gradients can promote compelling estimates of the phase in (3) and eventually the PSF in (2). However, as the phase gradients are generally collected by the WFS under coarse grid with less information, it is nontrivial to obtain high accurate estimates of by solely one frame of phase gradient measurements. Accordingly, how to obtain the high accurate phase gradients has drawn much attention in ground-based astronomy.
Under the frozen flow hypothesis on atmospheric turbulence, Jefferies and Hart [16] suggested reconstructing high-resolution phase gradients by a sequence of low-resolution ones. Concretely, let and be the sequences of low-resolution phase gradients in horizontal and vertical directions, respectively. The high-resolution phase gradients can be derived by where is a binary matrix down-sampling the high-resolution phase gradients into low-resolution ones, is an indicator matrix representing the telescope aperture, is a structured sparse matrix describing the aliasing (i.e., the displacements of phase screen in linear constant or nonlinear velocities) of the phase gradients in the th frame, and and are additive noise induced by down-sampling and aliasing.1 Fig. 1 is the diagrams of superposing a sequence of low-resolution phase gradients as a high-resolution one (see also [18]), which shows that the phase gradients at the superposed regions admit denser measurements than those simply collected by the WFS. With the low-resolution phase gradients at different overlapping locations, the high-resolution phase gradients can be reconstructed by the fusing or image super-resolution techniques (see e.g., [19] and references therein). The method in [16] was then improved in [18], [20] by applying the Tikhonov regularizer to the ill-posed linear inverse problem (4), and then the model in [20] was generalized in [21]. Overall, a common feature of the works in [18], [20], [21] is that the phase-gradients are estimated and then the phase is calculated. Such a model that processes these two-stage procedures separately are called a “phase-gradient model”.
To estimate the phase , instead of solving first (4) and then (3) sequentially, we can derive the phase directly by consolidating (3)–(4) as For brevity, we term the estimate of phase directly from the sequences of low-resolution phase gradients as the “phase model”. In Fig. 2 (see also [22]), we draw the flowchart of phase-gradient and phase models for reconstructing an image from an astronomical image system, in which the differences between these two kinds of models are clear. As easily observed, the phase-gradient model requires solving a couple of overdetermined and ill-posed linear systems (i.e., (3), (4)), whilst the phase model only needs to solve one (i.e., (5)). Hence, the phase model is favorable for reducing computational efforts, saving storages and diminishing the error accumulations. Moreover, because the property of phase is well-defined (e.g., piece-wise smooth) whilst that of phase gradients is ambiguous, the phase model (5) is typically easier to choose a regularizer and thus develop a tailor-made solver than the phase-gradient model (4). More details on the geometrical and statistical comparisons of phase, phase gradient and natural image are discussed in Section 3.1. We proposed in [22] the first phase model by using the tight frame to regularize the ill-posed linear system (5), and developed a tight-frame-regularized model for estimating the phase directly. Numerical results reported therein demonstrate that the phase model can lead to better recovery results than those phase-gradient models.
The rest of the paper is organized as follows. In Section 2, some basic definitions and propositions are stated for the sequel discussion. In Section 3, we propose a new phase model for estimating the PSF under frozen flow hypothesis where a single-layered atmospheric turbulence is considered. By reformulating the new phase model as a saddle-point problem, we present the implementation details when a primal–dual hybrid gradient method is applied to this reformulation. Furthermore, we test the new phase model on wavefront reconstruction in the case of single-layered atmospheric turbulence. In Section 4, we extend the new phase model to the case of multi-layered atmospheric turbulence. Finally, some concluding remarks are drawn in Section 5.
Section snippets
Preliminaries
We summarize some notations and definitions which will be used in sequel analysis.
For any and , let and denote the standard definitions of - and -norm of , respectively. The Huber norm of with a parameter , denoted by , is defined as where is the Huber function defined as (see e.g., [23], [24]) It follows from, e.g., [25], that the Huber function
A new phase model for estimating the PSFs
In this section, we propose a new phase model for estimating the PSFs when single-layered atmospheric turbulence is considered in ground-based astronomy.
Extension of (12) to the case of multi-layered atmospheric turbulence
In Section 3, we estimate the PSF when only single-layered atmospheric turbulence is considered. The effects of atmospheric turbulence are generally confined to the troposphere and lower stratosphere (e.g., 0–12 km in altitude), whilst rapidly decrease and disappear above 25 km in altitude [45]. It usually comprises layers, where each layer involves a collection of lenses of varying refractive indexes, moving at different velocities with respect to ground-based telescope aperture.
Concluding remarks
We develop a new phase model for estimating the point spread functions (PSFs) in ground-based astronomy in the cases of single- and multi-layered atmospheric turbulence. The new model uses the Huber norm regularizer, which can be easily reformulated as a saddle-point problem and efficiently solved by some well-developed algorithms in optimization community. We show the effectiveness of the new phase model by some numerical experiments, including the problem of spatially-varying deblurring. It
Acknowledgments
The author is grateful to Prof. Raymond H. Chan (City University of Hong Kong), Prof. Xiaoming Yuan (The University of Hong Kong), the anonymous referees for their suggestive comments on the manuscript, and Prof. James Nagy and Dr. Qing Chu for providing their codes of Refs. [18], [20]. This work was supported by the NSFC, China (11971003) and the Fundamental Research Funds for the Central Universities, China (ZYGX2019J090).
References (53)
- et al.
Wavelet deblurring algorithms for spatially varying blur from high-resolution image reconstruction
Linear Algebra Appl.
(2003) - et al.
Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations
(2006) - et al.
Deblurring Images: Matrices, Spectra, and Filtering
(2006) - et al.
Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity
(2010) - et al.
Total variation blind deconvolution
IEEE. Trans. Image Process.
(1998) - et al.
A convergent overlapping domain decomposition method for total variation minimization
Numer. Math.
(2009) - et al.
Restoration method for spatially variant blurred images
(2011) - et al.
Iterative statistical approach to blind image deconvolution
J. Opt. Soc. Amer. A
(2000) Introduction to Fourier Optics
(1996)Adaptive Optics in Astronomy
(1999)