Zernike polynomials applied to apparent solar disk flattening for pressure profile retrievals

We present a passive method for the retrieval of atmospheric pressure profiles based on the measurement of the apparent flattening of the solar disk as observed through the atmosphere by a spaceborne imager. This method was applied to simulated sunsets. It relies on accurate representation of the solar disk, including its limb darkening, and how its image is affected by atmospheric refraction. The Zernike polynomials are used to quantify the flattening in the Sun images. The inversion algorithm relies on a transfer matrix providing the link between the atmospheric pressure profile and a sequence of Zernike moments computed on the sunset frames. The transfer matrix is determined by a training dataset of pressure profiles generated from a standard climatology. The performance and limitations of the method are assessed by two test cases. Pressure profiles similar to the training dataset show that retrieval error can be up to 10 times smaller than the natural variability in the lower mesosphere, and up to 500 times smaller in the upper troposphere. Tests with other independent profiles emphasize the need for better representativeness of the training dataset.


Introduction
The measurement of refraction angles was one of the first remote sensing techniques in the field of spaceborne atmospheric science experiments.From the first radio occultation experiment during the Mariner IV mission that successfully retrieved pressure and temperature information in Mars' atmosphere (Nicks, 1967) to the pictures of the Sun as seen through the Earth's atmosphere by SOFIE (Solar Occultation For Ice Experiment, Gordley et al., 2009a) on board AIM (Aeronomy of Ice in the Mesosphere), this technique has proved to be a valuable alternative solution with respect to other remote passive methods like the detection of CO 2 or O 2 thermal emissions.Gordley et al. (2009b) have shown that one can achieve temperature retrievals with an accuracy of ±5 K up to 60 km provided that high spatial sampling of the solar disk is available.Their method is based on accurate solar upper and lower edge detection during sunset or sunrise in order to measure refraction angles and to retrieve temperature profiles.In this paper, we follow a different strategy.By taking advantage of the asymmetrical flattening of the apparent solar disk caused by the gradient in the air density profile, we retrieve atmospheric pressure from a quantified description of the whole solar disk deformation.For this purpose, we used Zernike polynomials to guarantee a unique orthogonal decomposition of the images.By simulating a wide variety of test cases, we built a lookup table to map the link between the Zernike moments of the images and the pressure state of the atmosphere.
This theoretical study is intended to show the capabilities of a simple technique based on Sun image analysis that can be applied to any solar imager performing occultations.It does not require highly efficient pointing stabilization, and it has the advantage of using all the pixels constituting the solar disk image, allowing for measurements with very high signal-to-noise ratios.
In the following sections, we first describe the steps leading to the simulation of solar images, including the solar limb darkening and the apparent solar disk flattening.Then we present the Zernike polynomials and we apply them to simulated sunset pictures for a large and non-redundant set of atmospheric pressure states, generated by a principal Published by Copernicus Publications on behalf of the European Geosciences Union.E. Dekemper et al.: Zernike polynomials applied to solar disk flattening component analysis performed on a representative pressure dataset.Finally, we assess the limitations of the method by testing it on two different datasets.

Virtual instrument definition
The goal of this study is to exploit the relation between Sun images and atmospheric pressure profiles.To do so, we developed a virtual but realistic imaging instrument whose essential characteristics match those of ACE, the Atmospheric Chemistry Experiment on board SCISAT-1 (see Bernath et al., 2005;McElroy et al., 2007, for more details), which possesses two imagers that are used for a twofold purpose (Gilbert et al., 2007).The imagers serve primarily as a pointing support to the payload by keeping the full Sun in their field of view (FOV).Second, as they observe solar occultations through the Earth's atmosphere, they provide information about the spectral transmittance related to the atmospheric extinction profiles in two different spectral regions.The filter's central wavelengths are identical to two of the seven channels of the SAGE II instrument (Mauldin et al., 1985): a visible channel at 525 nm weakly sensitive to aerosols and ozone, and a near infrared (NIR) channel at 1020 nm sensitive to aerosols as well.
The second channel was selected as the reference for our simulations in order to minimize the effect of trace gas absorption (mainly NO 2 and O 3 ) on the transmitted irradiance in case of actual application of the method to real data.For the sake of realism, we also kept ACE's pixel size and field of view: the image is sampled by an array of 128 × 128 squared pixels dividing a total field of view of 30 mrad, which is about three times larger than the full Sun apparent angular size.Finally, we also preserved ACE's orbital parameters, and placed our virtual sensor on a 650 km circular orbit.By choosing the wavelength, the field of view and the orbit altitude, all the parameters required to simulate Sun images are defined.

Solar limb darkening
The Sun image simulation requires a description of the solar limb darkening, i.e. the apparent radial decrease of intensity across the solar disk.This darkening is also wavelengthdependent and has been empirically parameterized by Neckel (2005).Defining µ = cos β where β is the angle of emission of a Sun ray with respect to the local vertical at the Sun surface, the ratio between the light intensity emitted at this angle I (µ) and the intensity of rays emitted with β = 0 (i.e. from the center of the apparent disk) can be fitted by a 5th-order polynomial: Table 1.Values attributed to the limb darkening coefficients as given in Neckel (2005) for wavelengths belonging to the spectral range 422-1100 nm.
a 00 = 0.75267 a 30 = 2.42234 a 01 = −0.265577 a 35 = −0.017117 a 10 = 0.93874 a 40 = −1.71150 a 11 = 0.265577 a 45 = 0.011977 a 15 = −0.004095 a 50 = 0.49062 a 20 = −1.89287 a 55 = −0.003347 a 25 = 0.012582 where the A i coefficients are functions of the wavelength λ.Keeping the notation of Neckel, the six coefficients are defined by the following relations, where the values of the a ij are given in Table 1:

Atmospheric refractive index
Atmospheric refraction is the key physical process involved in the apparent flattening of the solar disk when observed through Earth's atmosphere.The solar rays that penetrate the Earth's atmosphere experience an air concentration gradient giving rise to refraction.This process depends on the refractive index n of each particular layer.In this study, we assumed hydrostatic equilibrium, making n dependent on pressure (P ) and temperature (T ).CO 2 and water vapor were also considered according to the Ciddor formula (Ciddor, 1996), which is recommended by the International Association of Geodesy (IAG) for refractive index calculations in the visible and near infrared range in air.
In our model, the P , T and H 2 O profiles are based on the US Standard 1976 atmosphere.Only CO 2 has been changed to a more up-to-date value of 400 ppmv.The division of the atmosphere into layers also follows the US Standard atmosphere: the layers are 1 km thick below 25 km, 2.5 km thick between 25 and 50 km, and 5 km thick from 50 to 100 km.A spline interpolation was performed to avoid discontinuities over multiple layers.

1-D ray tracing problem
The problem of simulating solar pictures as acquired after propagation of the light in the atmosphere mostly relies on the knowledge of three parameters (highlighted in the following sentence): there is only one path for light rays leaving the solar photosphere with a given emission angle β to graze the Earth's surface at a minimal tangent height h and hit the detector with a given viewing angle θ.The geometry is illustrated in Fig. 1 (left panel) for rays lying in the plane containing the centers of the Sun, the Earth and the spacecraft.The formalism used in the present section and the next one is an adapted version of a previous study from Vokrouhlicky et al. (1993).
The trajectory of a light ray in a refractive medium is governed by Fermat's principle, stating that light always follows the fastest path.The associated distance L traveled between time t A and time t B is where n(s) is the refractive index at the coordinate s along the path.Furthermore, assuming spherical homogeneous atmospheric layers, the following conservation formula holds for any value of s: where c st is a constant value, r denotes the geocentric distance of a point located at a coordinate s along the trajectory, and ψ is the angle between the local vertical and the ray trajectory.Applying Eq. (3) to the point of the light trajectory situated at the smallest distance from the ground, at the altitude h, one can find the relation between the tangent height h and the viewing angle θ (angle between the light ray hitting the detector and the geocentric position vector of the spacecraft): where R E is the Earth radius and d ES the geocentric spacecraft distance.
In order to find the relation between the viewing angle θ, the tangent height h and the emission angle β, one can apply some trigonometric relations to the triangles of Fig. 1: where 2δ being the total refraction angle, ω, the Sun-Earth-satellite angle, and α, the angle between the Sun geocentric position vector and the heliocentric position vector of the emission spot.The refraction angle 2δ depends on all the atmospheric layers crossed by the light ray, from the top of the atmosphere (TOA) to the lowest altitude, i.e. the tangent height.The tangent height-dependent refraction angle can be obtained by using the following expression (see for instance Dalaudier et al., 2001, for more details on the computation of refraction angles in the atmosphere): The last step to solve Eq. ( 5) with respect to θ is to remove the dependence on α by applying trigonometric relations, leading to the emission angle associated to a light ray having reached the tangent height h and being seen by the instrument at the angle θ (more details can be found in Vokrouhlicky et al., 1993): where Up to now, only emission points lying in the Sun-Earthspacecraft plane (SES plane) have been considered, as illustrated by the left panel of Fig. 1.It was possible to address them with a single viewing angle (θ).For points located outside the SES plane, one can split their associated viewing angle into the vertical viewing angle θ and a horizontal viewing angle φ, allowing to univocally associate the image pixels with the couple of coordinates (φ, θ).The origin of the frame is chosen such that (0, 0) points towards the Earth center and the θ-axis belongs to the SES plane.The solar disk center receives coordinates (0, θ c ).This is illustrated in the right panel of Fig. 1.
When no refraction occurs along the light path, the solar limb darkening is a radial effect only, i.e. the emission angle β of the detected photons is the same for points located at an equal distance of the solar disk center (φ 2 + (θ − θ c ) 2 = c st ), allowing for a 1-D treatment.But when observed through an inhomogeneous refractive medium, the radial symmetry does not hold anymore and the relation between the light ray emission angle and the viewing angle has to be computed for any combination of θ and φ.
Equation ( 8) relates β to θ in the SES plane, whose intersection with the Sun forms a disk.It allows one to compute the limb darkening along the vertical direction for this particular meridian.We will first extend the validity of Eq. ( 8) to the other adjacent slices of the Sun, then correct it for the limb darkening along the horizontal direction (which does not suffer from any atmospheric effect as no horizontal atmospheric gradient is considered).This is more easily done by dividing the Sun into vertical slices of infinitesimal width.Every point belonging to the surface of a given slice is seen with the same φ angle and each slice is actually a disk whose radius R S (φ) is approximately computed by the following relation: Obviously, R S (φ) decreases with φ and vanishes when φ is equal to half of the angular size of the Sun.The limb darkening caused by a non-zero β angle within a φ-slice is obtained by using R S (φ) instead of R S in the calculation of ρ 1 in Eq. ( 8).One obtains an updated expression for β taking into account the reduction of the vertical angular size of the φ-slices: There remains to update the parameter µ involved in Eq. ( 1) with respect to β(φ, θ), and account for the limb darkening along the horizontal direction: To illustrate the effect of the atmospheric refraction on the apparent solar disk, we simulated images of the Sun for different Sun-Earth-spacecraft arguments ω.Two of them are shown in Fig. 2. When ω increases during the occultation, the light ray bending increases also due to higher air density, resulting in an asymmetrical flattening of the Sun image as light rays from the lower part of the Sun are more bent than those from the higher part.
The Sun images were computed according to the following procedure: -For a given state of the atmosphere (e.g.pressure, temperature, humidity), the total refraction angle 2δ(h) is computed for all tangent altitudes using Eq. ( 7) (the computation grid size is 10 m, with interpolation of the geophysical quantities on this finer grid).
-The intensity received by a pixel depends on the solar limb darkening of the region of the solar disk it is looking at (µ(φ, θ )), which in turn depends on the light ray emission angle (β(φ, θ )).Both quantities are computed on a much finer grid than the final digital image resolution (900 grid nodes per pixel typically).This is particularly important for pixels capturing the solar disk edge.Equation ( 4) is solved iteratively to match θ with h using the pre-computed 2δ(h).
-The final pixel intensity is obtained by averaging the radiance computed at the 900 grid nodes and using the assumed optical parameters of the virtual imager.

Zernike moments
As the goal of this work is to retrieve atmospheric pressure information from the apparent flattening of the solar disk, a unique quantitative description of the Sun shape is preferable.It is also necessary to have a method insensitive to spacecraft pointing uncertainties responsible for rotation and translation of the Sun in the FOV.For these reasons, we selected the Zernike polynomials.
The radial function is given by the following polynomial: One immediately notices that Z m * n = Z −m n .Figure 3 shows the Zernike polynomials computed on the unit disk up to n = 4.
For our application, Cartesian coordinates are preferred over polar, i.e.Z m n (ρ, α) is equivalent to Z m n (x, y) with x 2 + y 2 ≤ 1.By orthogonality of the Zernike functions, we have x 2 +y 2 ≤1 When applying Zernike polynomials to images, one needs to discretize the functions to account for the discrete nature of digital pictures.This may lead to a loss of orthogonality.This problem was already discussed by Wang and Silva (1980), and an extensive study of the orthogonality requirement of the discrete Zernike polynomials has been performed by Allen (2010) where different grid geometries (e.g.Cartesian, polar, random) and grid resolutions (number of grid nodes) were tested.The goal was to find the best discretization strategy, close to the perfect orthogonality of the continuous case.Allen showed, as expected, that a finer grid allows for better orthogonality.In the case of a Cartesian grid like in digital images, a rule of thumb for estimating the minimal number of grid points n p ensuring "decent orthogonality" between the n z first Zernike polynomials is to satisfy log 10 n p > 1.73 log 10 n z + 0.27.According to this rule, for applications where the 10 first Zernike polynomials are needed, the picture should be made of 100 pixels at least.In our case, images are cropped such that they closely fit the left and right edges of the solar disk, which makes the final array 45 × 45 pixels large.Not all of them belong to the unit circle -the pixels close to the corners are excluded from the Zernike domain -which leaves us with a bit more than 1500 pixels concerned.As will be shown later, we will not use many Zernike polynomials, so that the uniqueness of the decomposition will be guaranteed.
The moment of order n and repetition m of the image whose pixel intensities are represented by f (x, y) is given by the following expression: There is no guarantee that the solar disk will always remain exactly at the same position in the field of view.It can be shown that Zernike moments are invariant under translation and rotation (Khotanzad and Hong, 1990) meaning that the moment modulus is preserved under rotation.A consequence of this property is that, as , it is sufficient to rely on moments with m ≥ 0 only.
The translation invariance can be achieved by preprocessing the Sun image in order to have the centroid of the solar disk co-located with the center of the unit disk.It can be done by using the regular moments (Hu, 1962).In their discrete form, they are defined by the following expression: x p y q f (x, y).
In the Cartesian frame of the initial image, the centroid is located at (x, y) where In the translated coordinates frame (x , y ), the image is centered in such a way that f (x , y ) = f (x − x, y − y).
It can be shown (Khotanzad and Hong, 1990) that Zernike moments of order n = 1, being A −1 1 and A 1 1 , are proportional to the regular moments m 10 and m 01 .Both of them vanish in the translated frame, resulting in the disappearance of the first order Zernike moments as well.
In Fig. 4, Zernike moments up to order 4 are computed for the high and low Sun cases.The first moment A 0 0 is by definition the total integrated intensity on the whole solar disk and is always dominant.

Pressure profiles retrieval
The central problem addressed by this work is retrieving pressure profiles from images of the Sun during occultation.The measurements are represented by a number of Zernike moments A m n (ω) computed for each image I (ω) acquired successively during sunset or sunrise, ω being the Sun-Earth-spacecraft argument.From a given set of Zernike moments, the measurement vector a is written as A m j n j (ω i ) being the j -th selected moment computed on the i-th image.The atmospheric pressure profile will be referred to as the state vector p, made of 46 layers of constant pressure (US Standard atmosphere).The forward model can be simply written as where F describes the ray propagation, the image acquisition and the moments computation.Equation ( 20) is not solved directly due to the nonlinearity of operator F through the refraction process itself.Ill-conditioning might also be an issue: even if the number of image acquisitions is larger than the number of retrieved parameters, the problem might be ill-conditioned due to the redundancy between successive measurements.Furthermore, iterative approaches are not suitable as all the pictures of the setting or rising Sun need to be simulated at each iteration, which is not efficient regarding the processing time.

Training dataset
The first step is to have a clear representation of the variety of realistic pressure states.We used the COSPAR International Reference Atmosphere model (Rees et al., 1990) referred in the following to CIRA-86.It consists of monthly averaged zonal pressure data on an altitude grid spacing of 5 km from 20 to 120 km between −80 • and 80 • of latitude by steps of 10 • , forming a dataset of 204 profiles.To comply with the US Standard atmosphere, the pressure profiles were extrapolated down to 0 km by setting a common value of 101 300 Pa at ground for all profiles.They were then interpolated onto the US Standard atmosphere grid using a spline method.The dataset consists of a matrix of 204 realizations (the modified CIRA-86 profiles) of 46 variables (the number of atmospheric layers).Let p = (p 1 • • • p q ) be a pressure vector where q = 46 in this case.The whole dataset is represented by the matrix with p ij being the i-th realization of the j -th variable (i.e. the i-th pressure value of the j -th layer), where i = 1, . . ., r (r = 204 for CIRA-86) and j = 1, . . ., q.The CIRA-86 climatology does reflect the annual and worldwide variability of atmospheric pressure states.However, the averaging probably made the profiles smoother and removed the extrema.This makes the dataset a good starting point for testing the method but not yet suitable for real cases.

Principal component analysis
The profiles may contain redundant structures between each other and it could be useful to describe the pressure profiles with a reduced number of parameters.We used empirical orthogonal functions (EOF) (Lorenz, 1956), often referred to as principal component analysis (PCA) (Jolliffe, 2002) to solve these issues.The technique is based on finding orthogonal vectors (the principal axes) along the directions in which most of the dispersion of repeated realizations of the same phenomenon occurs.In the following paragraphs, we focus on finding the principal axes explaining most of the variation of the atmospheric pressure profile dataset and on estimating the error of such a reduced representation.
The dataset is first centered about its mean.By averaging the realizations of each variable, we define The centered dataset is denoted by P = P − P. Then the data are normalized by their standard deviation.The covariance matrix of P is = PT • P with its elements being the σ ij , i, j = 1, . . ., q; and 0 is the same matrix but with all offdiagonal elements removed.The factor 1/(q − 1) has been neglected.The normalized dataset is then whose elements are the p * ij = p ij − p j / √ σ jj .By performing a singular value decomposition (SVD) on the normalized dataset, the eigenvectors (the principal axes) are identified and their associated singular values are sorted in a decreasing order.Defining the matrices U, L and V such that U T U = V T V = I, and L is diagonal, we have P * = U L V T .The covariance matrix of P * follows: R = P * T P * = V L L V T = V V T , where the columns of V are the eigenvectors of R, and the λ ii of the diagonal matrix are its squared singular values.From now on, the columns of V will be called the principal axes (PA) of the normalized dataset.The principal components (PC) of a given pressure profile are its coordinates in the orthogonal basis formed by the PA.Computing the PCs of the entire normalized dataset, we find the matrix of PCs C of size r × q: Table 3. Coordinates of the pivots along the five first principal axes.The median mi and the standard deviation σ i of each distribution (i = 1, . . ., 5) were used to compute the pivots.They were chosen such as to span the spectrum with a reduced number of points.See Fig. 5 for a representation of the PC spectrum.Exploring the eigenvalues of the SVD, we find that the first ones already explain a very large portion of the overall dataset variability.Table 2 shows the eigenvalues and their cumulative importance for the 10 first PA.One notices that the first four PA are responsible for more than 99 % of the variance.
As the goal is to reduce the number of dimensions of the pressure profiles, we are looking for the minimum number of PA needed to allow for a reconstruction of the dataset profiles within acceptable errors.The satisfaction criterion is based on the global mean quadratic relative error (MQRE) that we required to be lower than 1 % for the full dataset.The global MQRE is computed this way: from a number of m PA (1 ≤ m ≤ q), the profiles are reconstructed using the associated m PCs: where V m is a q × m matrix containing the m first PA, and the subscript (m) denotes quantities obtained from the truncated basis of these m first PA.The MQRE ( ) is computed by averaging the relative deviation between the initial data (p ij ) and the reconstructed ones (p The last column of Table 2 provides the MQRE computed for up to 10 PA.As can be seen, five PA are needed to go below the 1 % threshold. When the 204 profiles of the CIRA-86 database are projected onto the five first principal axes, their PCs can be seen as a spectrum.Figure 5 shows the spectrum of PCs for each principal axis.In order to build a training database that would span a wider range of atmospheric states and to avoid redundancies between profiles (as is the case in CIRA-86), we used the principal component spectra to identify a small number of pivot values from which an exhaustive but non-redundant dataset can be generated.Table 3 shows the coordinates of the pivots along the five principal axes and how they have been computed.By combining the pivots, one can create a new database of 4 2 × 3 3 = 432 different pressure profiles that constitutes the training dataset that will be used in the retrieval algorithm.

The transfer matrix
The measurement simulations consist of 23 pictures of the setting Sun.Each picture corresponds to a given value of ω, the Sun-Earth-spacecraft argument: from 113.25 • (high Sun) up to 115.45 • (low Sun) by steps of 0.1 • .The images are then projected onto the Zernike polynomial basis and the Zernike moments are computed.Only a small number k of them are kept for the training.If there are l image pixels falling within the circular domain of the Zernike functions, ordered in a column vector f (ω i , p j ) (ω i is the i-th angle during the sunset and p j is the j -th pressure profile in the dataset), then in matrix formalism, the vector containing the Zernike moments is computed by where the rows of Z are the k selected Zernike polynomials evaluated at the l pixels lying in the circular domain of the image.The Zernike moments of each image and each pressure profile are then concatenated into a matrix A of size [23 • k × 432] to form the measurement matrix of the training dataset.
The link between the measurements of matrix A and the 5 principal components describing the pressure profiles arranged in the matrix C of size [5 × 432] is expressed by the transfer matrix X of size [5 × 23 • k]: Finally, the elements of X can be easily computed through the resolution of this overconstrained problem using a leastsquare regression method:

Test cases
The training dataset was obtained from non-redundant combinations of key values (the pivots) taken in the spectrum of principal components of the CIRA-86 dataset.In that sense, not a single pressure profile of the training dataset is identical to those contained in CIRA-86.However, it is correct to state that they are similar: the principal axes of both datasets are almost parallel.
As such, the pressure profiles of the CIRA-86 climatology constitute an optimal, though non-trivial, test case.The retrieved accuracy will unveil the upper limits of the method.Results obtained with independent pressure states should exhibit larger biases, particularly when they are not accurately represented by a linear combination of the principal axes of the training dataset.In other words, the adopted retrieval strategy performs better when the sounded atmospheric state lies in the multi-dimensional space defined by the principal axes of the training dataset.Poor results mainly indicate that the training dataset is not comprehensive enough.
We illustrated this discussion with two test cases.First, we applied the method to the 204 pressure profiles of the CIRA-86 climatology.This is actually the optimal situation as they are almost perfectly represented in the basis of the training dataset principal axes.Second, we picked up 30 pressure profiles from the ACE-FTS v3.0 products of the occultations observed in 2011 (Boone et al., 2005), constituting a truly independent test sample.The results are presented and discussed below.

CIRA-86
Figure 6a shows the 204 pressure profiles of the CIRA-86 climatology (dark gray) above the 432 profiles of the training dataset (pale gray).A solar occultation has been simulated for each of the CIRA-86 atmospheric states.Zernike moments have been computed for each of the 23 pictures constituting a solar occultation measurement in order to form the measurement vector.Finally, the transfer matrix is used to retrieve the pressure profile from the Zernike moments.Figure 6b shows the relative difference between the retrieved profiles and the original ones when only the first Zernike moment is used (A 0 0 , the solar disk total intensity), while the first and the fifth moments are used together in Fig. 6c (A 0 2 better captures the solar limb darkening).
The excellent results obtained below 20 km simply follow from the extrapolation of the CIRA-86 profiles from 20 km down to a common value at ground (101 300 Pa).This extrapolation was then inherently reflected in the training dataset, so that the least-square fit cannot be mistaken in this altitude range and it should be overlooked.Between 20 and 60 km, however, the average bias is almost 0 with increasing spread though.Figure 6d shows the mean relative difference ±1σ for the retrieval performed with two Zernike moments.In this altitude range, the reasons for the discrepancy are twofold.First, not all the profiles are well described by the training dataset principal axes (though very good for most of them).Second, the solar disk flattening is decreasingly measurable with higher altitudes, given the spatial resolution of the virtual imager.This is reflected in the Zernike moments as their amplitude becomes less and less dependent on the atmospheric state as the air density decreases.Figure 7a shows the evolution of the first Zernike moment during the sunset for all 204 pressure profiles of CIRA-86.Above 60 km, all pictures give the same moment amplitude, whatever the pressure in the upper layers.Below 60 km, different pressure profiles are characterized by distinct moment values.A measure of this distinction can be done by comparing the standard deviation of the moment population with the mean amplitude (Fig. 7b).As expected, this ratio falls down to zero above 60 km; the spatial resolution of the virtual imager is too small and the refraction too weak.

ACE-FTS pressure profiles sample
We picked 30 pressure profiles among the solar occultations performed by ACE-FTS in 2011.They originate from version 3.0 of the level 2 products, sampled on an altitude grid from 5.5 to 89.5 km by steps of 1 km.In order to fit with our atmospheric model, they were extrapolated towards 0 km (101 300 Pa) and 100 km (0.03 Pa) and interpolated on the US Standard grid.The ACE-FTS profiles were selected at random among those lying between the boundaries of the training dataset (Fig. 8a).Unlike the CIRA-86 profiles, they are truly independent from the training dataset.Their coherence with the training dataset can be assessed by computing their principal components and reconstructing the profiles from the 5 first principal axes (as was done to check the number of independent variables needed to describe the CIRA-86 profiles).Figure 8b shows the relative consequent error by doing so.Clearly, the basis extracted from the CIRA-86 climatology is not suitable for any kind of profile.As expected, by simulating sunsets for these 30 profiles, and applying the retrieval algorithm, we obtained poor results (Fig. 8c).

Discussion
Two test cases have been presented to assess the limitations of this method.The first one gives a flavor of what can be achieved at best with the assumed instrument specifications.For most of these atmospheric states, the pressure profile is retrieved up to the lower mesosphere with less than 5 % relative error (less than 1 % below 30 km).Better spatial resolution should probably lift the detection limit a few kilometers upward.
The second test case revealed a critical aspect of the retrieval strategy: the suitability of the training dataset (through its basis of PA) with the considered pressure profiles must be optimized.Clearly, the smoothed profiles of the CIRA-86 climatology did not turn out to be very representative of the ACE-FTS pressure profiles.Despite this lack of representativeness, the algorithm managed to deliver acceptable results on average for the lowest layers.For instance, at 20 km, the mean relative error is close to 1 ± 1.8 % whereas the 3σvariability of CIRA-86 is close to 25 %.

Impact of error sources
In real conditions, many sources of error should be included in the development with a major contribution from shot noise and detector dark current.The error budget can be computed by starting from the covariance matrix of each image: a diagonal matrix S f i of size [l × l] containing the shot noise of each pixel belonging to the circle domain.From there, S a i , the covariance matrix of a i , of size [k × k] follows: The covariance matrix of the measurement vector (the Zernike moments of the 23 images) is formed by concatenating the S a i to form a block-diagonal matrix S a .The covariance matrix of the retrieved PCs is then given by Finally, the covariance matrix of the error on the pressure profile is found using the same procedure in Eq. ( 25): In order to estimate the expected retrieval uncertainty, we simulated a test case where the brightest pixel of the solar pictures receives 10 000 counts and suffers from a dark current of 500 counts.Taking into account the associated shot noises, we used the above formula to derive the retrieval error by taking the square root of the diagonal elements of S P .
By examining the relative importance of this uncertainty to the profile, we find values ranging typically from 0.001 % for the lower layers, to 10 % for the highest ones.These typical numbers have to be compared with the natural variability of the pressure profiles as can be computed using the CIRA-86 dataset: a few percent at lower altitudes to a hundred percent at the highest ones.This exercise, though very simple, allows to conclude that, thanks to the high signal-to-noise ratio of the measurement vector (the Zernike moments use all the pixels of the solar disk resulting in a huge amount of signal), the retrieval uncertainty is expected to be 10 to 500 times smaller than the natural variability of the retrieved quantity.
Obviously, shot noise is the easiest and more predictable source of uncertainty.Other typical imager issues are related to the optics performance (e.g.blur, ghost images, distortion) or detector features (e.g.dead or hot pixels, pixel response non-uniformity).Each of them will require careful consideration in future developments but the inherent properties of the method lead to some general comments.
-Fixed patterns in the array of pixels or systematic distortions are image features which are captured by the Zernike moments.In some situations, only specific moments are affected so that a workaround is simply not to include them in the retrieval.
-If all moments are corrupted (like in the case of large distortion), comparing the distorted moment spectrum with the clean simulated case will allow for the isolation of what depends on the refraction (flattening), and what is due to the optics.The training dataset can then be adapted to take the instrumental effects into account.
-Non-systematic distortions like the presence of a cloud (PSC or PMC) in the field of view cannot be corrected and make the image unusable for the proposed method.However, when considering the question of reliable automatic flagging of solar pictures corrupted by clouds (Dodion et al., 2007), the Zernike moments turn out to be extremely efficient as their amplitude will be considerably affected in an unusual way with respect to the cloud-free situation.

Discussion and conclusions
This theoretical study was performed in order to assess the feasibility of retrieving atmospheric pressure profiles up to the lower mesosphere by means of solar disk image analysis.It relies on an accurate simulation of the apparent solar disk (including solar limb darkening and refractive effects), the use of the Zernike polynomials as orthogonal functions to measure the apparent flattening quantitatively, the reduction of the number of variables describing the pressure profiles by using a principal component analysis, and the construction of a comprehensive dataset to link the Zernike moments of the sunset pictures to the pressure profile.Its intrinsic practical advantages are the absence of strong pointing requirements (as the Sun is a very intense light source, one only needs to take pictures with sufficiently small integration times to have sharp solar disk images and to keep it in the instrument field of view for the complete sunset/sunrise), high signal to noise ratios (all the pixels sampling the solar disk are used) and potentially applicable to any solar occultation imager dataset provided that the images have been sent to ground in their entirety.The Zernike moment computation could even be performed onboard, requiring the download of a few numbers only (the Zernike moments themselves) instead of the full pixel map intensity.This could be an interesting feature for systems having small downlink capacity.The limitations of the method have been assessed by two test cases.Pressure profiles particularly well represented by the training dataset show retrieval relative error smaller than 1 ± 5 % up to 60 km.It constitutes what can be best achieved with this method, except if an instrument with higher spatial sampling is used.In that case, results could be slightly improved.Profiles that are not well represented by the training dataset are less well retrieved, though acceptable errors are found up to 20 km.The main point here is the need for a more comprehensive dataset with its principal axes allowing for correctly reconstructing most of the real pressure profiles.
Future investigations should address the issue of a training dataset showing better representativeness of the variability of natural pressure profiles.Other developments should also focus on the vertical resolution of the retrieved profiles, consider inhomogeneous atmospheric layers and account for possible instrumental effects.
This method was studied within the framework of the development of ALTIUS (Atmospheric Limb Tracker for the Investigation of the Upcoming Stratosphere; ALTIUS, 2013), a project of a new spaceborne instrument aimed at the measurement of atmospheric trace gases using hyperspectral images of the bright limb and solar/stellar occultations.

Fig. 1 .
Fig. 1.Left panel: illustration of the 1-D ray tracing problem for photon trajectories in the Sun-Earth-spacecraft plane.Right panel: illustration of the reference frame for vertical and horizontal viewing directions as seen from the instrument.The origin of the coordinate system is at the intersection of the axes and the black box illustrates a possible instrument field of view.Relative sizes of the Earth, the Sun and the pixels are not representative of the orbit and instrument parameters used in this study.

Fig. 2 .
Fig. 2. Simulated pictures of the Sun.Left picture shows a high Sun, right picture shows a low Sun already experiencing strong flattening.Notice the asymmetrical flattening.

Fig. 3 .
Fig. 3. Usual representation of the Zernike polynomials evaluated on the unit disk up to order n = 4.If m > 0, the real part is represented.If m < 0, the imaginary part is represented.
Fig. 4. First 15 Zernike moments computed on a high Sun image and a low Sun image (same as Fig. 2, arbitrary units).In the case of the high Sun image, moments with radial features (m = 0) are the only ones to play a role in characterizing this unrefracted solar disk.For the low Sun, more moments are needed to account for the asymmetrical flattening.

Fig. 5 .
Fig. 5. Distribution of the PCs of the CIRA-86 dataset along the five first principal axes.From these distributions, pivots can be defined in order to generate a training dataset of non-redundant elements.

Fig. 6 .
Fig. 6.Retrieval algorithm applied to the simulated sunsets for each of the 204 pressure profiles of the CIRA-86 climatology.The first panel shows the CIRA-86 profiles (dark gray) above the training dataset profiles (pale gray) on a logarithmic scale.The red lines emphasize the boundaries of the training dataset.The second panel shows the relative error of the retrieved pressure in function of altitude by using information from the first Zernike polynomial only.The third panel shows the relative error obtained with the first and fifth Zernike polynomials.The last panel shows the mean relative error µ ± 1σ for the two Zernike polynomials retrieval case.

Fig. 7 .
Fig. 7. Left panel shows the moment amplitude of the first Zernike polynomial computed for every simulated solar disk image and every profile of the CIRA-86 climatology.At the beginning of the occultation, the weakness of the refraction makes the Zernike moment insensitive to the particular pressure of the upper atmospheric layers.At lower altitudes, the moment starts taking different values for different pressure profiles.This is emphasized in the second panel where the standard deviation of the moment amplitude is compared to the mean value.

Fig. 8 .
Fig. 8. First panel shows the ACE-FTS sample of 30 pressure profiles (dark gray) above the training dataset (pale gray).The red lines emphasize the boundaries of the training dataset.Second panel shows the relative difference between the original profiles and those reconstructed from a projection onto the principal axes of the training dataset.The third panel shows the relative difference between the original profiles and the retrieved ones using the first Zernike moment.

E. Dekemper et al.: Zernike polynomials applied to solar disk flattening 2.5 2-D ray tracing problem
d SE with R S being the solar radius and d SE being the Sun-Earth distance, and ρ 2 = d ES /d SE .

Table 2 .
PCA applied to the CIRA-86 pressure dataset.The eigenvalues of the 10 first principal axes and their cumulative importance for capturing the overall dataset variability are given in the first two columns.The third column gives the total relative error made when the dataset is reconstruct by the use of the first PAs only.
* MQRE = mean quadratic relative error.To circumvent this difficulty, we preferred a retrieval strategy based on a lookup table.By simulating sunsets for a large number of plausible pressure profiles, we built a large set of direct relations between pressure profiles and associated Zernike moments.The pressure profiles used to create the lookup table constitute the training dataset.If the training dataset is sufficiently exhaustive with respect to the natural variation of atmospheric pressure, then any pressure profile should be well captured by the lookup table.