A phase model using the Huber norm for estimating point spread function under frozen flow hypothesis

https://doi.org/10.1016/j.cam.2021.113657Get rights and content

Abstract

In astronomical observation and space exploration, the light of an object in the outer space is inevitably distorted by the atmospheric turbulence around the Earth before it reaches the ground-based telescope aperture. The target object so observed is degraded by unknown point spread functions (PSFs). Estimating the PSFs is thus crucial for ground-based astronomy. In order to know the PSFs, we first need to estimate the phase, i.e., the wavefront aberration of incoming light at telescope aperture. In this paper, under the frozen flow hypothesis, we propose a new model for estimating the phase in the case of single- and multi-layered atmospheric turbulence. The new model, which deploys the Huber norm as the regularizer, can be cast as a saddle-point problem and solved efficiently by existing well-developed optimization algorithms. We test a spatially-varying image deblurring problem by using the estimated PSFs, which demonstrates the compelling numerical performance of the new phase model.

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 g(x,y)=k(x,y)f(x,y)+ε(x,y),where f is the image of a target object in the outer space and g is its observation through a ground-based telescope, k 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 k is known, the model (1) is a trivial linear inverse problem. However, the PSF k is generally unknown in practical applications, see e.g., [4], [5], [6], [7] and references therein. Thus, estimating PSF k is an essential step for further processing an astronomical image. Our discussion on the model (1) is under the assumption that the PSF k 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 k. 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 k can be expressed as the wavefront aberration of incoming light at the telescope aperture. Specifically, k(x,y)=|F1P(x,y)expιϕ(x,y)|2,where F denotes the Fourier transform and F1 is its inverse, P is the pupil function of telescope (1 inside the pupil and 0 otherwise), ι=1, 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 dr0 with d the diameter of the telescope and r0 the Fried parameter (see e.g., [9]).

To calculate the PSF k 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 sx and sy denote the phase gradients in horizontal and vertical directions, respectively. Then, the phase ϕ, and its gradients sx and sy are involved in the following overdetermined linear system (see e.g., [10], [11], [12], [13]) sxsy=DxDyϕ,where Dx, Dy are the derivative operators (see e.g., [14], [15]). Apparently, the high accurate estimates of the phase gradients (sx,sy) can promote compelling estimates of the phase ϕ in (3) and eventually the PSF in (2). However, as the phase gradients (sx,sy) are generally collected by the WFS under coarse grid with less information, it is nontrivial to obtain high accurate estimates of (sx,sy) by solely one frame of phase gradient measurements. Accordingly, how to obtain the high accurate phase gradients (sx,sy) 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 (sx,sy) by a sequence of low-resolution ones. Concretely, let {sxi}i=1m and {syi}i=1m be the sequences of low-resolution phase gradients in horizontal and vertical directions, respectively. The high-resolution phase gradients (sx,sy) can be derived by sxi=RWAisx+nxi,syi=RWAisy+nyi,i=1,2,,m,where RRl×n is a binary matrix down-sampling the high-resolution phase gradients into low-resolution ones, WRn×n is an indicator matrix representing the telescope aperture, AiRn×n 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 ith frame, and nxi and nyi 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 {(sxi,syi)}i=1m at different overlapping locations, the high-resolution phase gradients (sx,sy) 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 sxi=RWAiDxϕ+nxi,syi=RWAiDyϕ+nyi,i=1,2,,m.For brevity, we term the estimate of phase ϕ directly from the sequences of low-resolution phase gradients {(sxi,syi)}i=1m as the “phase model”. In Fig. 2 (see also [22]), we draw the flowchart of phase-gradient and phase models for reconstructing an image f 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 k 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 x=(x1,x2,,xn)Rn and p[1,), let xp(i=1n|xi|p)1p and xmax1in|xi| denote the standard definitions of p- and -norm of x, respectively. The Huber norm of x with a parameter μ>0, denoted by x1,μ, is defined as x1,μi=1nh(xi),xRn,where h:RR is the Huber function defined as (see e.g., [23], [24]) h(t)t22μ,if|t|<μ;|t|μ2,otherwise,tR.It follows from, e.g., [25], that the Huber function h

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 23 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)

  • ChanR.H. et al.

    Wavelet deblurring algorithms for spatially varying blur from high-resolution image reconstruction

    Linear Algebra Appl.

    (2003)
  • AubertG. et al.

    Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations

    (2006)
  • HansenP. et al.

    Deblurring Images: Matrices, Spectra, and Filtering

    (2006)
  • StarckJ.L. et al.

    Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity

    (2010)
  • ChanT.F. et al.

    Total variation blind deconvolution

    IEEE. Trans. Image Process.

    (1998)
  • FornasierM. et al.

    A convergent overlapping domain decomposition method for total variation minimization

    Numer. Math.

    (2009)
  • HadjS.B. et al.

    Restoration method for spatially variant blurred images

    (2011)
  • LamE.Y. et al.

    Iterative statistical approach to blind image deconvolution

    J. Opt. Soc. Amer. A

    (2000)
  • GoodmanJ.W.

    Introduction to Fourier Optics

    (1996)
  • RoddierF.

    Adaptive Optics in Astronomy

    (1999)
  • EllerbroekB.L.

    Efficient computation of minimum variance wavefront reconstructors with sparse matrix techniques

    J. Opt. Soc. Amer. A

    (2002)
  • GillesL.

    Order-N sparse minimun variance open-loop reconstructor for extreme adaptive optics

    Opt. Lett.

    (2003)
  • PoyneerL.A. et al.

    Fourier transform wavefront control with adaptive prediction of the atmosphere

    J. Opt. Soc. Amer. A

    (2007)
  • TalmiA. et al.

    Wavefront reconstruction from its gradients

    J. Opt. Soc. Amer. A

    (2006)
  • FriedD.L.

    Least-squares fitting a wave-front distortion estimate to an array of phase-difference measurements

    J. Opt. Soc. Amer.

    (1977)
  • HudginR.H.

    Wavefront reconstruction for compensated imaging

    J. Opt. Soc. Amer.

    (1977)
  • JefferiesS.M. et al.

    Deconvolution from wave front sensing using the frozen flow hypothesis

    Opt. Express

    (2010)
  • van DamM. et al.

    Performance of the Keck Observatory adaptive-optics system

    Appl. Opt.

    (2004)
  • ChuQ. et al.

    Iterative wavefront reconstruction for astronomical imaging

    SIAM J. Sci. Comput.

    (2012)
  • NgM.K. et al.

    Mathematical analysis of super-resolution methodology

    IEEE Signal Process. Mag.

    (2003)
  • J. Nagy, S. Jefferies, Q. Chu, Fast PSF reconstruction using the frozen flow hypothesis, in: Proceedings of the...
  • ChanR.H. et al.

    Point-spread function reconstruction in ground-based astronomy by l1lp model

    J. Opt. Soc. Amer. A

    (2012)
  • ChanR.H. et al.

    A phase model for point spread function estimation in ground-based astronomy

    Sci. China Ser. A

    (2013)
  • HuberP.J.

    Robust regression: Asymptotics, conjectures and Monte Carlo

    Ann. Statist.

    (1973)
  • LiW. et al.

    The linear 1 estimator and the Huber M-estimator

    SIAM J. Optim.

    (1998)
  • NesterovY.

    Smooth minimization of non-smooth functions

    Math. Program.

    (2005)
  • Cited by (0)

    View full text