Free Access
Issue
A&A
Volume 616, August 2018
Article Number A95
Number of page(s) 21
Section Celestial mechanics and astrometry
DOI https://doi.org/10.1051/0004-6361/201732537
Published online 31 August 2018

© ESO 2018

1. Introduction

Astrometry, which deals with the accurate and precise measurement of the positions and motions of celestial objects, is the oldest branch of observational astronomy, dating back at least to Hipparchus of Nicaea in 190 BC. Since, from its very beginnings, this branch of astronomy has required measurements over time to fulfil its goals, it could be considered the precursor of the currently fashionable “time-domain astronomy”, preceding it by at least 20 centuries. In recent years, astrometry has experienced a “coming of age” motivated by the rapid increase in positional precision allowed by the use of all-digital techniques and space observatories1.

A number of techniques have been proposed to estimate the location and flux of celestial sources as recorded on digital detectors, such as charge-coupled devices (CCD). In this context, estimators based on the use of a least squares (LS) error principle have been widely adopted (King 1983; Stetson 1987; Alard & Lupton 1998). The use of this type of decision rule has been traditionally justified through heuristic reasons. First, LS methods are conceptually straightforward to formulate based on the observation model of these problems. Second, they offer computationally efficient implementations and have shown reasonable performance (Lee & Van Altena 1983; Stone1989; Vakili & Hogg 2016). Finally, the LS approach was the classical method used when the observations were obtained with analog devices (Van Altena & Auer 1975; Auer & Van Altena 1978), which are well characterized by a Gaussian noise model for the observations. In the Gaussian case the LS is equivalent to the maximum likelihood (ML) solution (Chun-Lin 1993; Gray & Davisson 2004; Cover & Thomas 2012), and, consequently, the LS method was taken from the analog to the digital observational (Poisson noise model) setting.

Considering astrometry as an inference problem (of, usually, point sources), the astrometric community has been interested for a long time in understanding the fundamental performance limits (or information bounds) of this task (Lindegren 2010 and references therein). It is well understood by the community that the characterization of this precision limit offers the possibility of understanding the complexity of the task and how it depends on key attributes of the problem, like the quality of the observational site, the performance of the instrument (CCD), and the details of the experimental conditions (Mendez et al. 2013, 2014). On the other hand, it provides meaningful benchmarks to define the optimality of practical estimators in the process of comparing their performance with the bounds (Lobos et al. 2015).

Concerning the characterization and analysis of fundamental performance bounds, we can mention some works on the use of the parametric Cramér–Rao (CR) bound by Lindegren (1978), Jakobsen et al. (1992), Zaccheo et al. (1995), Adorf (1996), and Bastian (2004). The CR bound is a minimum variance bound (MVB) for the family of unbiased estimators (Rao, 1945; Cramér, 1946). In astrometry and joint photometry and astrometry, Mendez et al. (2013, 2014) have recently studied the structure of this bound, and have analyzed its dependency with respect to important observational parameters under realistic (ground-based) astronomical observing conditions. In this context, closed-form expressions for the Cramér–Rao bound were derived in a number of important settings (high pixel resolution and low and high signal-to-noise (S/N) regimes), and their trends were explored across different CCD pixel resolutions and the position of the object in the CCD array. As an interesting outcome of those studies, the analysis of the CR bound has allowed us to predict the optimal pixel resolution of the array, as well as providing a formal justification for some heuristic techniques commonly used to improve performance in astrometry, like dithering for undersampled images (Mendez et al. 2013, Sect. 3.3). Recently, an application of the CR bound to moving sources has been done by Bouquillon et al. (2017), indicating excellent agreement between our theoretical predictions, simulations, and actual ground-based observations of the Gaia satellite, in the context of the Gaia Ground Based Optical Tracking (GBOT) program (Altmann et al. 2014). The use of the CR bound on other applications is also of interest, for example, in assessing the performance of star trackers to guide satellites with demanding pointing constraints (Zhang et al. 2016), or to meaningfully compare positional differences from different catalogs (for example involving the Sloan Digital Sky Survey (SDSS) and Gaia, see Lemon et al. 2017). Finally, a formulation of the (non-parametric) Bayesian CR bound in astrometry, using the so-called “Van-Trees inequality” (Van Trees, 2004), has been presented by our group in Echeverria et al. (2016). This approach is particularly well suited for objects at the edge of detectability, and where some prior information is available, and has been proposed for the analysis of Gaia data for faint sources, or for those with a poor observational history (Michalik et al. 2015; Michalik & Lindegren, 2016).

From the perspective of astrometric estimators, Lobos et al. (2015) have studied in detail the performance of the widely adopted LS estimator. In particular extending the result in So et al. (2013), Lobos and collaborators derived lower and upper bounds for the mean square error (MSE) of the LS estimator. Using these bounds, the optimality of the LS estimator was analyzed, demonstrating that for high S/N there is a considerable gap between the CR bound and the performance of the LS estimator (indicating a lack of optimality of this estimator). This work showed that for the very low S/N observational regime (weak astronomical sources), the LS estimator is near optimal, as its performance closely follows the CR bound. The limitations of the LS method in the medium to high S/N regime proved in that work opens up the question of studying alternative estimators that could achieve the CR bound on these regimes, which is the main focus of this paper, as outlined below.

2. Contribution and organization

In this work we study the ML estimator in astrometry, motivated by its well-known optimality properties in a classical parametric estimation setting with independent and identically distributed measurements (i.i.d.; Kendall et al. 1999). We know that in the i.i.d. case this estimator is efficient with respect to the CR limit (Kay 1993), but it is important to emphasize that the observational setting of astrometry deviates from the classical i.i.d. case and, consequently, the analysis of its optimality is still an open problem. In particular, we face the technical challenge of evaluating its performance, a problem that, to the best of our knowledge, has not been addressed by the astrometric community. Concerning the independent but not identically distributed case, Bradley & Gart (1962) and Hoadley (1971) give conditions under which ML estimators are consistent2 and asymptotically normal3. Those conditions, however, are technically difficult to proof in this context.

The main challenge here is the fact that, as in the case of the LS estimator (Lobos et al. 2015), the ML estimator is the solution of an optimization problem with a nonconvex cost function of the data. This implies that it is not possible to directly compute the performance of the method. To address this technical issue, we extend the approach proposed by Fessler (1996) to approximate the variance and the mean of an implicit estimator solution of a generic optimization problem of the data through the use of a Taylor approximation around the mean measurement (see Theorem 1 below). Our extension considers high order approximations of the function that allow us not only to estimate the performance of the ML estimator through an explicit nominal value, but also provides a confidence interval around it. With this result we revisit the more general weighted least square (WLS) and ML methods providing specific upper and lower bounds for both methods (see Theorems 2 and 3). The main findings from our analysis of the bounds are two-fold: first we show that the WLS exhibits a sub-optimality similar to that of the LS method for medium to high S/N regimes discovered by Lobos et al. (2015) and, second, that the ML estimator achieves the CR limit for medium to high S/N and, consequently, it is optimal on those regimes. This last result is remarkable because, in conjunction with the result presented in Lobos et al. (2015), we are able to identify estimators that achieve the fundamental performance limits of astrometry in all the S/N regimes for the problem.

The paper is organized as follows: Sect. 3 introduces the background, preliminaries, and notation of the problem. Section 4 presents the main methodological contribution of this work. Section 5 presents the application to astrometry considering the WLS and the ML schemes. Finally, numerical analysis of these performance bounds are presented in Sect. 6 and the final remarks and conclusion are given in Sect. 7.

3. Preliminaries and background

We begin by introducing the problem of astrometry. For simplicity, we focus on the one-dimensional (1D) scenario of a linear array detector, as it captures the key conceptual elements of the problem4.

3.1. Relative astrometry as a parameter estimation problem

The main problem at hand is the inference of the relative position (in the array) of a point source. This source is modeled by two scalar quantities, the position of object x c ∈ ℝ in the array5, and its intensity (or brightness, or flux) that we denote by ∈ ℝ+ . These two parameters induce a probability distribution over an observation space that we denote by 𝕏. Formally, given a point source represented by the pair (x c, ), it creates a nominal intensity profile in a photon integrating device (PID), typically a CCD, which can be expressed by(1)

where ϕ(xx c,σ) denotes the one-dimensional normalized point spread function (PSF) and where σ is a generic parameter that determines the width (or spread) of the light distribution on the detector (typically a function of wavelength and the quality of the observing site, see Sect. 6; Mendez et al. 2013, 2014).

The profile in Eq. (1) is not measured directly, but it is observed through three sources of perturbations. First, an additive background noise that captures the photon emissions of the open (diffuse) sky, and the noise of the instrument itself (the read-out noise and dark-current Janesick 2001, 2007; Howell 2006; McLean 2008), modeled by i in Eq. (2). Second, an intrinsic uncertainty between the aggregated intensity (the nominal object brightness plus the background) and actual measurements, which is modeled by independent random variables that follow a Poisson probability law. Finally, we need to account for the spatial quantization process associated with the pixel resolution of the PID as specified in Eqs. (2) and (3). Modeling these effects, we have a countable collection of independent and non-identically distributed random variables (observations or counts) {I i : i ∈ ℤ, where I i ∼ Poisson(λ i(x c, )), driven by the expected intensity at each pixel element i, given by(2)

and(3)

where 𝔼{} is the expectation value of the argument and {x i : i ∈ ℤ} denotes the standard uniform quantization of the real line array with resolution ∆x > 0, that is, x i+1x i = Δx for all i ∈ ℤ. In practice, the PID has a finite collection of measurement elements (or pixels) I 1,…,In, then a basic assumption here is that we have a good coverage of the object of interest, in the sense that for a given position x c (4)

At the end, the likelihood (probability) of the joint observations In = (I1,…,In) (with values in ℕ n ) given the source parameters (x c, ) is given by(5)

where denotes the probability mass function (PMF) of the Poisson law (Gray & Davisson, 2004).

Finally, if is assumed to be known6, the astrometric estimation is the task of defining a decision rule τn() : ℕn → Θ with Θ = ℝ being the parameter space, where given an observation I n the estimated position is given by c(In ) = τn (In ).

3.2. The Cramér–Rao bound

In astrometry the Cramér–Rao bound has been used to bound the variance (estimation error) of any unbiased estimator (Mendez et al. 2013, 2014). In general, let In be a collection of independent observations that follow a parametric PMF defined on ℕ. The parameters to be estimated from In will be denoted in general by the vector . Let τn(In ) : ℕ n → Θ be an unbiased estimator7 of , and be the likelihood of the observation given . Then, the Cramér–Rao bound (Rao 1945; Cramér 1946) establishes that if(6)

then, the variance (denoted by Var), satisfies that(7)

where is the Fisher information matrix given by, ∈i, j ∈ {1,…,m},(8)

In particular, for the scalar case (m = 1), we have that for all θ ∈ Θ, (9)

where is the collection of all unbiased estimators and . For astrometry, Mendez et al. (2013, 2014) have characterized and analyzed the Cramér–Rao bound, leading to

3.3. Achievability and performance of the LS estimator

Concerning the achievability of the CR bound with a practical estimator, Lobos et al. (2015, Proposition 2) have demonstrated that this bound cannot be attained, meaning that for any unbiased estimator τn (·) we have that(11)

where In follows the Poisson PMF in Eq. (5).

This finding should be interpreted with caution, considering its pure theoretical meaning. This is because Eq. (11) does not exclude the possibility that the CR bound could be approximated arbitrarily closely by a practical estimation scheme. Motivated by this refined conjecture, Lobos et al. (2015) proposed to study the performance of the widely adopted LS estimator8 with the goal of deriving operational upper and lower performance bounds of its performance that could be used to determine how far this scheme could depart from the CR limit. Then, from this result, it was possible to evaluate the goodness of the LS estimator for concrete observational regimes. To bound the performance of the LS estimator, the challenge was that τ LS(In ) is an implicit function of the data (where no close-form expression is available) and, consequently, Lobos et al. (2015) derived a result to bound the estimation error and the variance of τ LS(In ). We can briefly summarize the main result presented in Lobos et al. (2015), Theorem 1) by saying that under certain mild sufficient conditions (that were shown to be realistic for astrometry), there is a constant δ > 0 (that depends on the observational regime, in particular the S/N) and a nominal variance , which is determined in closed-form in the result, from which it is possible to bound (τ LS(In ) by the simple expression(12)

where(13)

We note that when δ is small, tightly determines the performance of the LS estimator, and its comparison with can be used to evaluate the goodness of the LS estimator for astrometry. Based on a careful comparison, it was shown in Lobos et al. (2015, Sect. 4) that, in general, is close to for the small S/N regime of the problem. However, for moderate to hig h S/N regimes, the gap between and becomes quite significant9.

These unfavorable findings for the LS method have motivated us to study alternative schemes that could potentially approach better the Cramér–Rao bound for the rich observational context of medium to high S/N regimes. This will be the focus of the following sections, where in particular we explore the performance of the ML and WLS estimators, thus extending and generalizing the analysis done for the LS estimator by our group presented in Lobos et al. (2015).

4. Bounding the performance of an implicit estimator

Before we go to the case of the WLS and the ML estimators, we present a general result that bounds the performance of any estimator that is the solution of a generic optimization problem. Let us consider a vector of observations In = (I1,…In) ∈ ℝn and a general so-called cost function J(α, In ). Then the estimation of xc from the data is the solution of the optimization problem(14)

where α represents the position of the object in the context of astrometry. As in our previous work (Lobos et al. 2015), the challenge here is that this estimator is implicit because no closed-form expression of the data which solves Eq. (14) is assumed. In particular, this implies that both the variance and the estimation error of τJ(In) cannot be determined directly. To address this technical issue, we extend the approach proposed by Fessler (1996) to approximate the variance and the mean of an implicit estimator solution of a problem described by Eq. (14) through the use of a Taylor approximation around the mean measurement, that is, .

More precisely, we assume that J(α,In ) has a unique global minimum at τJ (In ), and that it has a regular behavior, so its partial derivatives are zero, that is,(15)

Then we can obtain τJ (In ), by a first order Taylor expansion around the mean by(16) (17)

where t ∈ [0,1] is fixed but unknown10. For simplicity, Eq. (17) can be written in matrix form as(18)

where denotes the row gradient operator and is the residual error of the Taylor expansion. From Eq. (18) we can readily obtain the expression for its variance(19) (20) (21)

In Eq. (21) we recognize two terms: that captures the linear behavior of τJ (·) around and γJ (n) which reflects the deviation from this linear trend. It should be noted that the above expression does not depend on τ(In) itself, but on its partial derivatives evaluated at the mean vector of observations. Then in the adoption of this approach to estimate Var{τJ (In )}, a key task is to determine .

In an effort to analyze both the linear and non-linear aspects of a general intrinsic estimator solution to Eq. (14), the following result offers sufficient conditions to determine in closed-form, and to bound the magnitude of the residual term γJ (n) in Eq. (21).

The proof of this result and the expression for in Eqs. (24) and (25) are presented in detail in Appendix A.

Revisiting the equality in Eq. (21), Theorem 1 provides general sufficient conditions to bound the residual term γJ (n) and by doing that, a way of bounding the variance of τJ (I) which is the solution of Eq. (14). In particular, it is worth noting that if the ratio , then Eq. (25) offers a tight bound for . In this last context, (called the nominal value of the result) provides a very good approximation for .

On the application of this result to the WLS and ML estimators, we will see that the main assumption in Eq. (22) is satisfied in both cases (see Eqs. (B.8) and (C.6)), and from that plays an important role to approximate the performance of ML and WLS in a wide range of observational regimes. In addition, the analysis of the bias in Eq. (24) shows that these estimators are unbiased for any practical purpose and, consequently, contrasting their performance (estimation error ) with the CR bound is a meaningful way to evaluate optimality.

5. Application to astrometry

In this section we apply Theorem 1 to bound the variances of the ML and WLS estimators in the context of astrometry. Following the model presented in Sect. 3.1, In = (I1,...,In)T denotes the measurements acquired by each pixel of the array and where each of them follows a Poisson distribution given by(26)

as expressed by Eqs. (2) and (3).

5.1. Bounding the variance of the WLS estimator

The WLS estimator, denoted by τWLS(In) in Eq. (28), is implicitly defined through a cost function given by(27)

where is a weight vector, and α is a general source position parameter. Specifically we have that(28)

Applying Theorem 1 we obtain the following result:

The proof of this result and the expressions for ϵWLS(n) and βWLS(n) in Eqs. (29) and (30), respectively, are elaborated in Appendix B. This result offers concrete expressions to bound the bias as well as the variance of the WLS estimator. For the bias bound in Eq. (29), it will be shown that ϵWLS(n) is very small (orders of magnitude smaller than xc) for all the observational regimes explored in this work and, consequently, the WLS can be considered an unbiased estimator in astrometry, as would be expected. Concerning the bounds for the variance in Eq. (30), we will show that for high and moderate S/N regimes the ratio , and consequently in this context is a precise estimator of . For the very small S/N the results offers an admissible interval around the nominal value . Therefore in any context shows to be a meaningful approximation for the performance of the WLS.

5.2. Bounding the variance of the ML estimator

The ML estimator, denoted by τWLS(In) in Eq. (34), is implicitly defined through a cost function(33)

where α is a general source position parameter. Specifically, given an observation In we have that(34)

Applying Theorem 1 we obtain the following result:

The proof of this result and the expressions for ϵML(n) and βML(n) in Eqs. (35) and (36), respectively, are elaborated in Appendix C.

6. Numerical analysis

In this section we evaluate numerically the performance bounds obtained in Sect. 5 for the WLS and ML estimators, and compare them with the astrometric CR bound in Proposition 1. The idea is to consider some realistic astrometric conditions to evaluate the expressions developed in Theorems 2 and 3 and their dependency on important observational conditions and regimes. As we shall see, key variables in this analysis are the trade-off between the intensity of the object and the noise represented by the S/N value, and the pixel resolution of the CCD.

6.1. Experimental setting

We adopt some realistic design and observing variables to model the problem (Mendez et al. 2013, 2014). For the PSF, analytical and semi-empirical forms have been introduced (see for instance the ground-based model in King 1971 and the space-based models by King 1983 or Bendinelli et al. 1987). In this work we will adopt a Gaussian PSF, that is, in Eq. (3), and where σ is the width of the PSF and is assumed to be known. This PSF has been found to be a good representation for typical astrometric-quality ground-based data (Méndez et al. 2010). In terms of nomenclature, measured in arcsec, denotes the full-width at half-maximum parameter, which is an overall indicator of the image quality at the observing site (Chromey, 2016).

The background profile, represented by {i, i = 1,…,n} in Eq. (2), is a function of several variables, like the wavelength of the observations, the moon phase (which contributes significantly to the diffuse sky background), the quality of the observing site, and the specifications of the instrument itself. We will consider a uniform background across pixels underneath the PSF, that is, i = for all i. To characterize the magnitude of , it is important to first mention that the detector does not measure photon counts [e] directly, but a discrete variable in “Analog to Digital Units (ADUs)” of the instrument, which is a linear proportion of the photon counts (McLean 2008). This linear proportion is characterized by the gain of the instrument G (scaling value) in units of [e/ADU]. We can define F/G and B/G as the intensity of the object and noise, respectively, in the specific ADUs of the instrument. Then, the background (in ADUs) depends on the pixel size Δx arcsec as follows:(38)

where fs is the (diffuse) sky background in ADU arcsec−1, while D and RON2, both measured in e, model the dark-current and read-out-noise of the detector on each pixel, respectively. The first component in Eq. (38) is attributed to the site, and its effect is proportional to the pixel size. On the other hand, the second component is attributed to errors of the PID (detector), and it is pixel-size independent. This distinction is central when analyzing the performance as a function of the pixel resolution of the array (see details in Mendez et al. 2013, Sect. 4). More important is the fact that in typical ground-based astronomical observation, long exposure times are considered, which implies that the background is dominated by diffuse light coming from the sky, and not from the detector Mendez et al. (2013), Sect. 4).

For the experimental conditions, we consider the scenario of a ground-based station located at a good site with clear atmospheric conditions and the specification of current science-grade CCDs, where fs = 1502.5 ADU arcsec−1, D = 0, RON = 5 e arcsec (equivalent to arcsec), and G = 2, e ADUe−1 (with these values we will have B = 313 ADU for Δx = 0.2 arcsec using Eq. (38)). In terms of scenarios of analysis, we explore different pixel resolutions for the CCD array Δx ∈ [0.1,0.65] measured in arcsec, and different signal strengths ∈ {1080; 3224; 20 004; 60 160}, measured in e, which corresponds to S/N ∈ ∼ {12,32,120,230}. Increasing implies increasing the S/N of the problem, which can be approximately measured by the ratio F̃/B̃. On a given detector plus telescope setting, these different S/N scenarios can be obtained by appropriately changing the exposure time (open shutter) that generates the image.

thumbnail Fig. 1.

Relative performance of the bias (as measured by ) stipulated by Theorem 1 for the WLS estimator (left side, Eq. (29)) and the ML estimator (right side, Eq. (35)). Results are reported for different values of the source flux ∈ {1080, 3224, 20 004, 60, 160}, all in e (top to bottom symbols respectively), as a function of the detector pixel size. The 0% level corresponds to having achieved no bias.

Open with DEXTER

6.2. Bias analysis

Considering the upper bound terms ϵWLS(n) and ϵML(n) for the bias error obtained from Theorems 2 and 3 for the WLS and ML, respectively, Fig. 1 presents the relative bias error for different S/N regimes and pixel resolutions. In the case of the ML estimator, the bounds for relative bias error are very small in all the explored S/N regimes and pixel resolutions meaning that for any practical purposes this estimator is unbiased as expected from theory (Gray & Davisson, 2004). For the case of the WLS, we observe that from medium to high S/N the relative error bound obtained is very small, meanwhile at low S/N unbiasedness cannot be fully guaranteed from the bound in Eq. (29). In general, our results show that both WLS and ML are unbiased estimators for astrometry in a wide range of relevant observational regimes (in particular from medium to high S/N) and, consequently, it is meaningful to analyze the optimality of these estimators in comparison with the CR bound in those regimes.

In the following sections, we move to the analysis of the variance of the WLS and ML with particular focus on the medium to high S/N regimes across all pixel resolutions using the performance bounds derived in Eqs. (30) and (36), respectively.

6.3. Performance analysis of the WLS estimator

In this section, we evaluate numerically the expression derived in Theorem 2 to bound the variance of the WLS estimator in Eq. (30). For that we characterize the admissible regime predicted for the variance of the WLS estimator, that is, the interval(39)

for S/N ∈ ∼ {12,32,120,230} and Δx ∈ [0.01,0.65] arcsec. In these bounds, we recognize its central value (or nominal value) in Eq. (31) and the length of the interval 2βWLS(n) that is determined in closed form for its numerical evaluation in Eqs. (A.2) and (B.11). We note that 2βWLS(n) can be considered an indicator of the precision of our result to approximate the variance of the WLS in astrometry.

thumbnail Fig. 2.

Range of the square root of the variance performance (in miliarcsecond = mas) for the WLS method in astrometry using uniform weights (equivalent to the LS method) predicted by Theorems 1 and 2, Eq. (30). Results are reported for different representative values of and across different pixel sizes (top-left to bottom-right): ∈ {1080; 3224; 20 004; 60 160} e.

Open with DEXTER

6.3.1. Revisiting the uniform weight case

To begin the analysis, we consider the setting of uniform weights across pixels, that is, the case of the LS estimator and, without loss of generality, we locate the object in the center of the field of view11, which can be considered a reasonable scenario to represent the complexity of the astrometry task. At this point, it is important to remind the reader that from the analysis of in Remark 2, the nominal value is equivalent to the CR bound when the wi are selected as a function of the true position in Eq. (32). In view of this observation, the selection of non-uniform weights can be interpreted as biasing the estimation to a particular area of the angular space, which goes in contradiction to the essence of the inference problem that estimates the position with no prior information, and only relies on the measured counts. From this interpretation, revisiting the LS estimator is an important first step in the analysis of the WLS framework.

On the specifics, the boundaries of the interval in Eq. (39) and its nominal values are illustrated in Fig. 2 for the different observational regimes. In addition, Fig. 2 shows the CR bound as a reference to evaluate the optimality of the LS scheme across settings. We observe that for the low S/N ∼ 12 regime, the nominal values precisely match the CR bound, however, our result is not conclusive as the interval around the nominal performance is significantly large. This is the regime where our result is not conclusive regarding the performance of the LS estimator. In the regime S/N ∈ (30,50) (top right panel on Fig. 2), we notice an important reduction in the range of admissible performance, and our result becomes more informative and meaningful. In this context, the nominal values are very close to the CR bound, and we could assert that the LS estimator offers sufficiently good performance in the sense that it is very competitive with the MVB. Importantly, when we move to the regime of relatively high S/N and very high S/N (from 100 to 300), our results are very accurate in predicting the performance of the LS method, and we find that the gap between the CR and the nominal value is very significant (the deviation from the MVB is 16% and 30% for S/N 120 and 230 at Δx = 0.2 arcsec, respectively). This last result confirms one of the main findings presented in Lobos et al. (2015), who showed that for medium to high S/N the LS estimator is suboptimal with respect to the MVB. In Fig. 2 we also show the square root of the empirical variance () with respect to the empirically-determined mean position c (using the WLS estimator), all as derived from the simulations, showing good consistency with our predictions.

thumbnail Fig. 3.

Worst-case discrepancies in Eq. (40) for the WLS estimator using the weights set indexed by the positions . Results are reported for two S/N scenarios, namely = 20 004 e (left) and = 60 160 e (right) (right), and across different pixel sizes.

Open with DEXTER

6.3.2. Non-uniform weight case

The sub-optimality of the WLS scheme from moderate to high S/N seen in Fig. 2 can be extended for any arbitrary selection of a fixed set of weights, as would be expected. Given that the space of weights’ selection is literally unlimited, we use the insight obtained from Remark 2 that states that a selection of weights can be interpreted as an specific prior on the position of the object where the optimum, but unfortunately unknown, selection (achieving the CR bound) of weights in Eq. (32) is an explicit function of the unknown position of the object. Then, we consider a finite set of positions {xc,1,…,xc,M} that uniformly partition the field of view, and their corresponding weight sets {wi(xc,1) : 1 = 1,…,n}, {wi(xc,1) : 1 = 1,…,n}, {wi(xc,M) : 1 = 1,…,n} using Eq. (32) to cover a representative collections of weights for the problem of astrometry.

Then, for a specific selection of weights in our admissible collection (attributed to a prior belief regarding the position of the object in the field of view), we evaluate the worse discrepancy between (which is the most favorable expression for the variance predicted from Theorem 2), and the CR bound in Proposition 1, across a collection of presumed object positions in the range of positions

where denotes the center of the array (which, as indicated at the beginning of Sect. 6.3.1, is equal to the true object position xc) and is the dispersion parameter of the PSF. The idea of using this worst-case difference is justified from the fact that in this parameter estimation problem we do not know the position of the object and, consequently, the optimality of any WLS estimator should be evaluated in the worse-case situation, as the scheme should be able to estimate the position of the object in any scenario (position). More precisely, for a given weight selection {wi, i=1,…,n}, we use the worst-case discrepancy(40)

For this analysis, we note that both and are functions of the position xc, Δx, and S/N.

Figure 3 illustrates the worst-case discrepancy in Eq. (40) for the medium and high S/N regimes where Theorem 2 provides an accurate and meaningful prediction of the performance of the WLS method, that is, for S/N ∈ {120,230}, and across Δx ∈ [0.05,0.7] arcsec. The discrepancy is quite significant, on the order of 37% and 60% in the range for Δx ∈ [0.1,0.3] arcsec for S/N 120 and 230, respectively.

To refine the worst-case analysis presented in Fig. 3, and to evaluate in more detail the sensitivity of the discrepancy indicator given by Eq. (40), we evaluate the discrepancy of WLS using the weights associated to (the center position of the array) with respect to the CR bound associated to the positions , to study how the discrepancy (measuring the non-optimality of the method) evolves when the adopted position moves far from the prior imposed by the WLS in the center of the array. Figure 4 illustrates this behavior, where it is possible to see that the discrepancy is very sensitive and proportional to the misassumption of the object position, where the worst-case discrepancy in the maximum achievable location precision is on the order of 40% for pixel sizes in the range [0.1, 0.6] arcsec for S / N ~ 120, and about 60% for pixel sizes in the range [0.1, 0.6] arcsec for S / N ~ 230. These worst-case scenarios happens in both cases when the object is located the farthest from the prior assumption, as would be expected.

The main conclusion derived form this CR bound analysis is that, independent of the weight selection adopted, as long as the weights are fixed, the WLS estimator is not able to achieve the CR bound in all observational regimes. More precisely, the discrepancy (measuring the non-optimality) in the less favorable case of an hypothetical and feasible position of the object is very significant, in the range of 40%–60% for the important regime of high and very high S/N.

thumbnail Fig. 4.

Performance discrepancies (measuring the non-optimality) of the WLS estimator using the center position as a prior for the weight selection with respect to the CR bound obtained for the true object positions . Results are reported for two S/N scenarios, namely = 20 004 e (left) and = 60 160 e (right), and across different pixel sizes.

Open with DEXTER

6.4. Performance analysis of the ML estimator

In this section we perform the same analysis done for the WLS in Sect. 6.3, but using the result in Theorem 3. In particular, we consider the admissible regime for the variance of the ML estimator given by

in Eq. (36), where the nominal value in this case, in Eq. (37), is precisely the CR bound, while the length of the interval 2βML(n), is given by Eqs. (A.2) and (C.9) in closed form.

Considering an object located in the center of the array, that is, xc = x0 the performance curves are presented in Fig. 5 for S/N ∈ {12, 32, 120, 230} and Δx ∈ [0.1, 0.65] arcsec. First, we note that there is a significant difference in the predictions of our methodology for the ML estimator in comparison with what we predict in the WLS case. In fact, the results of our approach are very precise for the determination of the variance of the ML estimator in all the regimes, from small to high S/N, which is remarkable. More important it is the fact that, from these findings, we observe that the performance deviation from the MVB in the worst case (small S/N) is very small (see Table 1, first row), while for any practical purposes the variance of the ML estimator achieves the CR limit for all the other regimes, from medium to high S/N, which is a numerical corroboration of the optimality of the ML estimator in astrometry, as predicted theoretically by Theorems 1 and 3. In Fig. 5 we also show the square root of the empirical variance (Var(τM̂L)) with respect to the empirically-determined mean position x̂c (using the ML estimator), all as derived from the simulations, showing good consistency with our predictions.

Complementing this analysis, we conducted the same comparison considering different positions for the object within the array obtaining the same trends and conclusions. To summarize these results, Fig. 6 shows the value , which is an indicator of the quality of the estimator (the smaller the better) for different scenarios of the position of the object . In particular, for all the evaluated positions, the relative discrepancy is bounded (relative to the CR bound) by values in the range of 0.025% and 0.012% for pixel resolution in the range Δx ∈ [0.1, 0.2] arcsec for S/N = 120 and S/N = 230, respectively. Finally, Table 1 presents the relative discrepancy for all the range of S/N values considered in this study for the case Δx = 0.2 arcsec.

We conclude from this analysis that the ML estimator in nearly optimal for the medium, high, and very high S/N regimes across pixel resolutions, achieving the MVB for astrometry. This is an interesting result, since it lends further support to the adoption of these types of estimators for very demanding astrometric applications, as has been done in the case of Gaia (Lindegren, 2008). We note that Vakili & Hogg (2016) reach the same conclusion regarding the optimality of the ML method in comparison with the MVB, through simulations of two-dimensional (2D) CCD images using a broad set of Moffat PSF stellar profiles. While their results are purely empirical, it is interesting that they test the ML using a different PSF from ours, and in a 2D scenario, and yet they reach the same conclusions. More recently, Gai et al. (2017) have tested (also empirically) the ML method in a 1D scenario (similar to ours), but in the context of a Gaia-like PSF. They find that the ML is unbiased (in agreement with our results, see Fig. 1) and by comparing two implementations of the ML they conclude that they predict self-consistent and reliable results over a broad range of flux, background, and instrument response variations. It would still be interesting to compare the performance of those implementations against the CR MVB in order to further test our theoretical predictions.

thumbnail Fig. 5.

Range of the square root of the variance performance (in miliarcsecond = mas) for the ML method in astrometry as predicted by Theorems 1 and 3, Eq. (36). Results are reported for different representative values of F̃ and across different pixel sizes (top-left to bottom-right): F̃ ∈ {1080, 3224, 20 004, 60 160} e.

Open with DEXTER
Table 1.

Performance quality of the ML estimator relative to the Cramér–Rao bound expressed in terms of the indicator from the result in Theorem 3.

6.5. Comments on an adaptive WLS estimator for astrometry

In Sect. 5.1 we presented results that offer a nominal prediction for the variance of the WLS method through Eq. (30), which turns out to be very accurate in the regime from medium to high S/N as shown in Sect. 6.3. Interestingly, Remark 2 tells us that these nominal values precisely match the CR limit for an optimal selection of weights given in Eq. (32), namely wi ∼ 1/λi(xc) for all i = 1, . . . , n (compare Eqs. (31) and (37)). As this selection is unfeasible, because it requires the knowledge of xc (see the expression in Eq. (2)), we can approximate this value by a noisy version of it, considering the fact that the expected value of the observations Ii that we measure at pixel i is λi(xc) using our model in Eq. (4). Therefore Ii can be interpreted as a noisy version of λi(xc) and(41)

can be seen as a noisy version of the ideal weight . Adopting this data-driven weighting approach, we would have an adaptive WLS method as the weights are not fixed but instead they are a function of the data {I1 : i = 1, . . . , n}. This selection of weights can be interpreted as an empirical version of the optimal weights that achieves the CR bound. Then, the problem reduces to solve(42)

where(43)

Figure 7 presents the performance of this scheme for the same regimes we have been exploring in this work, supporting the conjecture that this selection of weights resembles the optimal weight selection and in fact achieves MSE performances that are surprisingly close to the CR bound in all the observational regimes.

Our results in this sub-section show that some commonly adopted weighting schemes, using the analogous idea to Eq. (41), are a very good data-driven choice for methods that employ a WLS scheme, such is the case, for example, of the well-know PSF stellar fitting program (including astrometry) Dominion Astrophysical Observatory Photometry (DAOPHOT), described in Stetson (1987, Eq. (10)).

thumbnail Fig. 6.

Performance optimality of the ML estimator (computed as ) for different positions of the target object in the array, as a function of pixel resolution. The left panel shows the case = 20 004 e, the right panel the case = 60 160 e.

Open with DEXTER
thumbnail Fig. 7.

Performance comparison between the of the adaptive WLS estimator and σCR(n), both in mas. Results are reported for different F̂ and across different pixel sizes: (top-left to bottom-right) F̃ ∈ {1080, 3224, 60 160}e.

Open with DEXTER

6.6. Undersampled case

So far we have analyzed the performance of the ML and WLS estimators in regimes where the PSF (assumed to have a FWHM = 1 arcsec) is well sampled by the detector. In this section we analyze the undersampled case, where we have only a few pixels that capture most the flux of the source. This scenario would be typical of nearly diffraction-limited images taken, for example, from space, or with adaptive optics (AO) system from the ground.

Table 2 shows the performance of the bounds presented in Theorems 2 and 3 in the undersampled regime, where the PSF is concentrated in only a few pixels (order of five to nine pixels), and where the flux is relatively large. For the space-based case, the background is reduced almost to the pure instrument noise (left numbers in the table), while in the ground-based case the background would be larger, but we would still have a very large Strehl ratio. From the table it can be seen that the performance is nearly optimal with the ML estimator, even in the severely undersampled regime (small FWHM values). On the other hand, the WLS deteriorates its performance very quickly as the image size decreases, reaching up to nearly 45% when the FWHM is on the order of the detector pixel size.

We emphasize that the above results are computed in the ideal situation of a uniform pixel response function (and a perfect flat-field correction), whereas in reality a non-uniform pixel response has a critical impact on the Cramér–Rao bound as demonstrated by Adorf (1996), Fig. 1). In addition, in the case of stare-mode space-based images acquired under very small background, charge transfer inefficiency effects due to time-varying charge traps in the detector would also blur (and systematically shift) the images, deteriorating the ultimate astrometric precision above the Cramér–Rao limit (for details see, e.g., Bristow et al. 2006, especially their Fig. 10).

7. Summary, conclusions, and outlook

In this paper we study the performance of the WLS and ML estimators for relative astrometry on digital detectors subject to Poisson noise, in comparison with the best possible attainable precision given by the CR bound. Our study includes analytical results and numerical simulations under realistic observational conditions to help us to corroborate our theoretical findings.

By generalizing the proposal by Fessler (1996) we are able to obtain, for the first time, close-form expressions for the variance and the mean of implicit estimators (as is the particular case of the WLS and ML schemes), which can be computed directly from the data (see Theorem 1, in particular Eqs. (24) and (25), and Appendix A). When applying this result to astrometry with digital detectors, we are able to bound both the bias and the variance of the relative position of a celestial source on a CCD array as a function of all the relevant parameters of the problem (see Eqs. (29) and (30) or Eqs. (35) and (36) for the WLS and ML estimators, respectively). We verified that the bias of the WLS and ML methods are negligible in all the observational regimes explored in this paper (see Fig. 1).

A careful analysis of our predictions confirms earlier results by (Lobos et al. 2015; for the LS method) in that the WLS method is, in general, sub-optimal (in comparison with the MVB given by the CR result), especially at high and very high S/N (see the two bottom panels on Fig. 2). However, a judicious data-driven selection of weights (called “adaptive” WLS method by us, Sect. 6.5), improves the performance of the WLS substantially (see Fig. 7). This is an interesting result, given the widespread use and simple numerical implementation of the WLS method.

The ML method is found to have both a smaller bias than the WLS method (compare left and right panels of Fig. 1), although the bias on both methods is already quite small, and a tight correspondence to the MVB throughout the entire range of S/N regimes explored in this paper (Fig. 5). Therefore, the ML estimator for astrometry is consistently optimal, and should be the estimator of choice for high-precision applications. This paper, along with Lobos et al. (2015), completes an in-depth study of the performance of commonly used estimators in astrometry using PIDs, and sets the stage for the development of codes that could efficiently implement astrometric ML estimators on 2D detectors, incorporating also the simultaneous measurement of fluxes, as explored by Gai et al. (2017).

Table 2.

Performance quality of the ML and WLS estimators relative to the nominal bound expressed in terms of the indicator and , respectively, in the undersampled regime (space-based and ground-based are on the left and right, respectively).


1

See, e.g., Reffert (2009), Fig. 1a in Høg (2017) for an overview spanning more than 2000 years of astrometry, and Benedict et al. (2016) for a summary of the contributions from the Hubble Space Telescope (HST) (fine guide sensors), and, of course, the exquisite prospects from Prusti et al. (2016), with applications ranging from fundamental astrophysics (Van Altena, 2013; Cacciari et al. 2016), to cosmology (Lattanzi, 2012).

2

As the sample size increases, the sampling distribution of the estimator becomes increasingly concentrated at the true parameter value.

3

More precisely, the distribution around the true parameter approaches a normal distribution as the sample size grows.

4

The analysis can be extended to the two-dimensional case as presented in Mendez et al. (2013).

5

This captures the angular position in the sky and it is measured in seconds of arc (arcsec thereafter), through the “plate-scale”, which is an optical design feature of the instrument plus telescope configuration.

6

The joint estimation of photometry and astrometry is the task of estimating both from the observations (see Mendez et al. 2014).

7

In the sense that, for all , .

8

This is the solution of , with λi(α) = giα + i, (α) being a generic variable representing the astrometric position, gi(·) is given by Eq. (3), and where arg min represents the argument that minimizes the expression. More details are presented in Lobos et al. (2015).

9

In particular, for the very high S/N regime and assuming , Lobos et al. (2015), Proposition 3) show that this gap reaches the condition

10

It follows that .

11

It is important to remember that the CR bound is a function of the value of the parameter to be estimated, in this case the position xc (see the Fisher information in Eq. (10)).

12

Considering that .

13

The derivation of this identity is presented in Appendix B.2

14

Considering that .

15

The derivation of this result is presented in Appendix C.2

Acknowledgments

SE acknowledges support from CONICYT-PCHA/MagisterNacional/2016 – 22160840. The authors want to thank the support of the Advanced Center for Electrical and Electronic Engineering, AC3E, Basal Project FB0008, PIA ACT1405, and the Chilean Centro de Excelencia en Astrofisica y Tecnologias Afines (CATA) BASAL PFB/06 from CONICYT. The authors also acknowledge support by CONICYT/FONDECYT grants # 1151213, 1170044, and 1170854. RAM also acknowledges support from the Project IC120009 Millennium Institute of Astrophysics (MAS) of the Iniciativa Cientifica Milenio del Ministerio de Economia, Fomento y Turismo de Chile. Finally, we thank the anonymous referee for his or her careful reading of our manuscript, and for pointing out the exploration of the relevant space-based undersampled regime that led to the discussion in Sect. 6.6.

References

  1. Adorf, H.-M. 1996, in Astronomical Data Analysis Software and Systems V, 101, 13 [NASA ADS] [Google Scholar]
  2. Alard, C., & Lupton, R. H. 1998, ApJ, 503, 325 [NASA ADS] [CrossRef] [Google Scholar]
  3. Altmann, M., Bouquillon, S., & Taris, F. 2014, Proc. SPIE, 9149, 91490P [CrossRef] [Google Scholar]
  4. Auer, L., & Van Altena, W. 1978, AJ, 83, 531 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bastian, U. 2004, GAIA Technical Note, 2004 BASNOCODE [Google Scholar]
  6. Bendinelli, O., Parmeggiani, G., Piccioni, A., & Zavatti, F. 1987, AJ, 94, 1095 [NASA ADS] [CrossRef] [Google Scholar]
  7. Benedict, G. F., McArthur, B. E., Nelan, E. P., & Harrison, T. E. 2016, PASP, 129, 012001 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bouquillon, S., Mendez, R., Altmann, M., et al. 2017, A&A, 606, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Bradley, R. A., & Gart, J. J. 1962, Biometrika, 49, 205 [CrossRef] [Google Scholar]
  10. Bristow, P., Kerber, F., & Rosa, M. 2006, in The 2005 HST Calibration Workshop: Hubble After the Transition to Two-Gyro Mode, 299 [Google Scholar]
  11. Cacciari, C., Pancino, E., & Bellazzini, M. 2016, Astron. Nachr., 337, 899 [NASA ADS] [CrossRef] [Google Scholar]
  12. Chromey, F. R. 2016, To Measure the Sky: An Introduction to Observational Astronomy (Cambridge: Cambridge University Press) [Google Scholar]
  13. Chun-Lin, L. 1993, IAU Symp., 156, 113 [NASA ADS] [Google Scholar]
  14. Cover, T. M., & Thomas, J. A. 2012, Elements of Information Theory (New York: Wiley) [Google Scholar]
  15. Cramér, H. 1946, Scand. Actuar. J., 1946, 85 [CrossRef] [Google Scholar]
  16. Echeverria, A., Silva, J. F., Mendez, R. A., & Orchard, M. 2016, A&A, 594, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Fessler, J. A. 1996, IEEE Trans. Image Process., 5, 493 [NASA ADS] [CrossRef] [Google Scholar]
  18. Gai, M., Busonero, D., & Cancelliere, R. 2017, PASP, 129, 054502 [NASA ADS] [CrossRef] [Google Scholar]
  19. Gray, R. M., & Davisson, L. D. 2004, An Introduction to Statistical Signal Processing (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
  20. Hoadley, B. 1971, Ann. Math. Stat., 1977 [CrossRef] [Google Scholar]
  21. Høg, E. 2017, ArXiv e-prints [arXiv:1707.01020] [Google Scholar]
  22. Howell, S. B. 2006, Handbook of CCD Astronomy (Cambridge: Cambridge University Press), 5 [Google Scholar]
  23. Jakobsen, P., Greenfield, P., & Jedrzejewski, R. 1992, A&A, 253, 329 [NASA ADS] [Google Scholar]
  24. Janesick, J. R. 2001, Scientific Charge-Coupled Devices (Bellingham: SPIE Press), 83 [Google Scholar]
  25. Janesick, J. R. 2007, Photon Transfer (San Jose: SPIE Press) [CrossRef] [Google Scholar]
  26. Kay, S. M. 1993, Fundamentals of Statistical Signal Processing. Vol 1, Estimation Theory (Englewood Cliffs: Prentice-Hall) [Google Scholar]
  27. Kendall, M., Stuart, A., Ord, J., & Arnold, S. 1999, Kendall’s Advanced Theory of Statistics. Vol. 2A (London: Hodder Arnold Publication) [Google Scholar]
  28. King, I. R. 1971, PASP, 83, 199 [NASA ADS] [CrossRef] [Google Scholar]
  29. King, I. R. 1983, PASP, 95, 163 [NASA ADS] [CrossRef] [Google Scholar]
  30. Lattanzi, M. 2012, Mem. Soc. Astron. It. 83, 1033 [NASA ADS] [Google Scholar]
  31. Lee, J.-F., & Van Altena, W. 1983, AJ, 88, 1683 [NASA ADS] [CrossRef] [Google Scholar]
  32. Lemon, C. A., Auger, M. W., McMahon, R. G., & Koposov, S. E. 2017, MNRAS, 472, 5023 [NASA ADS] [CrossRef] [Google Scholar]
  33. Lindegren, L. 1978, IAU Colloq., 48, 197 [NASA ADS] [Google Scholar]
  34. Lindegren, L. 2008, Gaia DPAC Public Document GAIA-C3-TN-LU-LL-078 [Google Scholar]
  35. Lindegren, L. 2010, ISSI Scientific Reports Series, 9, 279 [NASA ADS] [Google Scholar]
  36. Lobos, R. A., Silva, J. F., Mendez, R. A., & Orchard, M. 2015, PASP, 127, 1166 [NASA ADS] [CrossRef] [Google Scholar]
  37. McLean, I. S. 2008, Electronic Imaging in Astronomy: Detectors and Instrumentation (New York: Springer Science & Business Media) [Google Scholar]
  38. Méndez, R. A., Costa, E., Pedreros, M. H., et al. 2010, PASP, 122, 853 [NASA ADS] [CrossRef] [Google Scholar]
  39. Mendez, R. A., Silva, J. F., & Lobos, R. 2013, PASP, 125, 580 [NASA ADS] [CrossRef] [Google Scholar]
  40. Mendez, R. A., Silva, J. F., Orostica, R., & Lobos, R. 2014, PASP, 126, 798 [NASA ADS] [Google Scholar]
  41. Michalik, D., & Lindegren, L. 2016, A&A, 586, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Michalik, D., Lindegren, L., Hobbs, D., & Butkevich, A. G. 2015, A&A, 583, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Prusti, T., De Bruijne, J., Brown, A. G., et al. 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Rao, C. R. 1945, Bull. Calcutta Math. Soc., 37, 81 [Google Scholar]
  45. Reffert, S. 2009, New Astron. Rev., 53, 329 [NASA ADS] [CrossRef] [Google Scholar]
  46. So, H. C., Chan, Y. T., Ho, K., & Chen, Y. 2013, IEEE Signal Process. Mag., 30, 162 [NASA ADS] [CrossRef] [Google Scholar]
  47. Stetson, P. B. 1987, PASP, 99, 191 [NASA ADS] [CrossRef] [Google Scholar]
  48. Stone, R. C. 1989, AJ, 97, 1227 [NASA ADS] [CrossRef] [Google Scholar]
  49. Vakili, M., & Hogg, D. W. 2016 ArXiv e-prints [arXiv:1610.05873] [Google Scholar]
  50. Van Altena, W. F. 2013, Astrometry for Astrophysics: Methods, Models, and Applications (Cambridge: Cambridge University Press) [Google Scholar]
  51. Van Altena, W. F., & Auer, L. 1975, in Image Processing Techniques in Astronomy (Berlin: Springer), 411 [NASA ADS] [CrossRef] [Google Scholar]
  52. Van Trees, H. L. 2004, Detection, Estimation, and Modulation Theory, Part I: Detection, Estimation, and Linear Modulation Theory (New York: Wiley) [Google Scholar]
  53. Zaccheo, T., Gonsalves, R., Ebstein, S., & Nisenson, P. 1995, ApJ, 439, L43 [NASA ADS] [CrossRef] [Google Scholar]
  54. Zhang, J., Hao, Y. C., Wang, L., & Long, Y. 2016, PASP, 128, 035003 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A

Proof of Theorem 1

We begin presenting the expressions for to complete the statement of the result.(A.1) (A.2)

where(A.3) (A.4)

and, finally,(A.5)

Proof of Theorem 1: using the chain rule in the cost function J(α, In) and taking the partial derivative of both sides in Eq. (15), we have that(A.6)

Thus, we have n equations with one unknown, and it holds for any In. Defining the operators(A.7)

of dimensions 1 × 1 and 1 × n, respectively, we can express Eq. (A.6) in matrix form as(A.8)

Assuming that the matrix ∇20 J(τ(In), In) is non singular, we can calculate ∇τ(In) from Eq. (A.8),(A.9)

Finally, using Eq. (A.9), evaluating at Īn, and then replacing in Eq. (21), we have that(A.10)

Moving into the residual term γJ(n) in Eq. (21) captured by βJ(n), we must consider the variance of the error function Var{e(I̅n, InI̅n)} and the covariance . For the first, we have that(A.11)

On the other hand, for the covariance, using the main assumption in Eq. (22), it is clear that(A.12)

From this,(A.13)

Finally, replacing Eqs. (A.11) and (A.13) in the definition of γJ(n), we have that(A.14)

For the bias expression of the result in Eq. (24), using the hypothesis in Eq. (23), we can take expectation at both sides of Eq. (18) to obtain that(A.15)

Appendix B

Proof of Theorem 2

Proof: the proof and, in particular, the derivation of , βWLS(n) and ϵWLS(n) simply reduces to a straightforward application of Theorem 1. For that we need to first validate the assumptions of Theorem 1. If we begin with Eq. (A.9),(B.1)

and then we calculate the gradient terms on the right hand side (RHS) of Eq. (B.1) for our WLS context, it follows that(B.2) (B.3) (B.4)

Following Eq. (B.1), we need to evaluate Eqs. (B.3) and (B.4) at α = τWLS(I̅n). For that, we have the following12:(B.5)

Then we will use the following result.

Proposition 2

Under the assumption of a Gaussian PSF, τWLS(I̅n) = xc.

We note that this proposition is the second assumption used in Theorem 1. Using this proposition, we obtain that (B.6) (B.7)

Finally, applying Eqs. (B.6) and (B.7) in Eq. (A.9), we have that(B.8)

which offers the decomposition needed for the application of Theorem 1 (Eq. (22)). For the value of in Eq. (A.5), since the observations are independent and follow a Poisson distribution, we have that(B.9)

Then if we replace Eqs. (B.6), (B.7), and (B.9) in Eq. (A.5), we have that(B.10)

On the other hand, the expression for βWLS(n) and ϵWLS(n) can be determined from the evaluation of Eqs. (A.2) and (A.1), respectively. Looking at them, the problem reduces to determine the key term . For that, if we use Fessler (1996), Eq. (17) we can obtain the following identity13,(B.11)

which concludes the result.

B.1. Proof of Proposition 2

Proof: using the function , we need to show that the minimum is reached only at α = xc. From this, we have that h(α) ≥ 0 and it achieves its minimum at xc. To prove uniqueness, let us assume that there is another position at which h is zero. Then(B.12)

The last identity is not possible, because if we use a Gaussian PSF there is at least one i ∈ {1, . . . , n} such that .

Proof of Eq. (B.11)

Proof: if we recall Fessler (1996), Eq. (17) and consider JWLS(α, In) as the cost function, we have that(B.13)

where from the definition of JML (α, In) we have that(B.14) (B.15) (B.16) (B.17)

Concerning is just the i-th component of the gradient in Eq. (B.1), then we use Eqs. (B.3) and (B.4)(B.18)

Finally, replacing Eqs. (B.14)–(B.18) in Eq. (B.13), and evaluating in α = τWLS(I̅n + t(InI̅n)), we obtain the desired result(B.19)

Appendix C

Proof of Theorem 3

Proof: again the proof and the derivation of , βML(n) and ϵML(n) reduce to apply Theorem 1. First, we need to validate the assumption of Theorem 1. Beginning with the equality in Eq. (B.1), it follows that(C.1) (C.2)

To evaluate these two expression at α = τML(I̅n) as required in Eq. (B.1), we use that14 (C.3)

Then we will use the following result.

Proposition 3

Under the assumption of a Gaussian PSF, α = τML(I̅n) = xc.

We notice again that this proposition is the second assumption used in Theorem 1. From this proposition, it follows that(C.4) (C.5)

Finally, we apply Eqs. (C.4) and (C.5) in Eq. (A.9) to obtain that(C.6)

which shows that the sufficient condition in Eq. (22) of Theorem 1 is satisfied. To compute the value in Eq. (A.5), we have that(C.7)

since the observations are independent and follow a Poisson distribution. Then, replacing Eqs. (C.4), (C.5), and (C.7) in Eq. (A.5), we have that(C.8)

which resolves the identity in Eq. (37). Finally βML(n) and ϵML(n) come from evaluating Eqs. (A.2) and (A.1) in this ML context. For that we only need to determine . Using Fessler (1996), Eq. (17), we can obtain the identity15 (C.9)

which concludes the result.

Proof of Proposition 3

Proof: let us consider the function given by(C.10)

We note that(C.11)

where applying the first order condition y = λi, ∀i ∈ {1, . . . , n}. Returning to our problem in Eq. (C.3) where λi = I̅i(xc) and yi = λi(α), it is clear, considering the Gaussian profile in PSF, that(C.12)

which concludes the result.

Proof of Eq. (C.9)

Proof: if we recall Fessler (1996, Eq. (17)) and consider JML(α, In) as the cost function, we have that(C.13)

where from the definition of JML(α, In) we have that(C.14) (C.15) (C.16) (C.17)

Concerning , this expression is the i-th component of the gradient in Eq. (B.1), then we use Eqs. (C.1) and (C.2)(C.18)

Finally, replacing Eqs. (C.14)–(C.18) in Eq. (C.13), and evaluating in α = τML(I̅n + t(InI̅n)), we obtain the desired result(C.19)

All Tables

Table 1.

Performance quality of the ML estimator relative to the Cramér–Rao bound expressed in terms of the indicator from the result in Theorem 3.

Table 2.

Performance quality of the ML and WLS estimators relative to the nominal bound expressed in terms of the indicator and , respectively, in the undersampled regime (space-based and ground-based are on the left and right, respectively).

All Figures

thumbnail Fig. 1.

Relative performance of the bias (as measured by ) stipulated by Theorem 1 for the WLS estimator (left side, Eq. (29)) and the ML estimator (right side, Eq. (35)). Results are reported for different values of the source flux ∈ {1080, 3224, 20 004, 60, 160}, all in e (top to bottom symbols respectively), as a function of the detector pixel size. The 0% level corresponds to having achieved no bias.

Open with DEXTER
In the text
thumbnail Fig. 2.

Range of the square root of the variance performance (in miliarcsecond = mas) for the WLS method in astrometry using uniform weights (equivalent to the LS method) predicted by Theorems 1 and 2, Eq. (30). Results are reported for different representative values of and across different pixel sizes (top-left to bottom-right): ∈ {1080; 3224; 20 004; 60 160} e.

Open with DEXTER
In the text
thumbnail Fig. 3.

Worst-case discrepancies in Eq. (40) for the WLS estimator using the weights set indexed by the positions . Results are reported for two S/N scenarios, namely = 20 004 e (left) and = 60 160 e (right) (right), and across different pixel sizes.

Open with DEXTER
In the text
thumbnail Fig. 4.

Performance discrepancies (measuring the non-optimality) of the WLS estimator using the center position as a prior for the weight selection with respect to the CR bound obtained for the true object positions . Results are reported for two S/N scenarios, namely = 20 004 e (left) and = 60 160 e (right), and across different pixel sizes.

Open with DEXTER
In the text
thumbnail Fig. 5.

Range of the square root of the variance performance (in miliarcsecond = mas) for the ML method in astrometry as predicted by Theorems 1 and 3, Eq. (36). Results are reported for different representative values of F̃ and across different pixel sizes (top-left to bottom-right): F̃ ∈ {1080, 3224, 20 004, 60 160} e.

Open with DEXTER
In the text
thumbnail Fig. 6.

Performance optimality of the ML estimator (computed as ) for different positions of the target object in the array, as a function of pixel resolution. The left panel shows the case = 20 004 e, the right panel the case = 60 160 e.

Open with DEXTER
In the text
thumbnail Fig. 7.

Performance comparison between the of the adaptive WLS estimator and σCR(n), both in mas. Results are reported for different F̂ and across different pixel sizes: (top-left to bottom-right) F̃ ∈ {1080, 3224, 60 160}e.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.