Nanosatellites: The Next Big Chapter in Atmospheric Tomography

Abstract

Nanosatellite technology opens up new possibilities for earth observation. In the next decade, large satellite constellations will arise with hundreds, up to thousand of satellites in low earth orbit. A number of satellites will be equipped with rather low-cost sensors, such as GNSS receivers, suited for atmospheric monitoring. However, the future evolution in atmospheric science leans not only on densified observing systems but also on new, more complex analysis methods. In this regard, tomographic principles provide a unique opportunity for sensor fusion. The difficulty in performing the conversion of integral measurements into 3D images is that the signal ray path is not a straight line and the number of radio sources and detectors is limited with respect to the size of the object of interest. Therefore, the inverse problem is either solved linearly or iterative nonlinear. In this chapter, an overview about the individual solving techniques for the tomographic problem is presented, including strategies for removing deficiencies of the ill-posed problem by using truncated singular value decomposition and the L-curve technique. Applied to dense nanosatellite formations, a new quality in the reconstruction of the 3D water vapor distribution is obtained, which has the potential for leading to further advances in atmospheric science.

1. Introduction

For the reconstruction of two- or three-dimensional structures from integral measurements of atmospheric excess phase, e.g. as obtained from signals of the Global Navigation Satellite Systems (GNSS), a technique called atmospheric tomography has been invented. The basic mathematics behind was introduced by Johann Radon in 1917 and is therefore also known as the Radon transform [1]. Its first realization in form of an axial scanning computer tomograph for cross-sectional imaging of the human body was awarded in 1979 with the Nobel prize for medicine [2, 3]. Around the same time, the tomography concept was utilized for applications in geosciences. One of the very first results was communicated by [4] in 1977, who describe a three-dimensional inversion method for simultaneous reconstruction of seismic body wave velocities and epicenter coordinates.

According to [5] the basic mathematical principle of tomography is described as follows:

$$f_s={\int }_{S}g\left(s\right)\cdot \mathit{ds}$$E1

where $f_s$ is the integral function, $g\left(s\right)$ is the object property function, and $\mathit{ds}$ is a small element of the ray path $S$ along which the integral is determined.

In atmospheric tomography, $g\left(s\right)$ is typically replaced by refractivity $n$, which is connected to signal velocity $v_p$ by the constant speed of light $c$.

$$v_p=\frac{c}{n}$$E2

In vacuum $n=1$, in matter is $n\ne 1$, i.e. the signal is slower or faster than the speed of light, dependent on the electric and magnetic properties of the medium. The integral measure $f_s$ is usually replaced by excess phase¹ or signal travel time which can be converted to excess phase by multiplying with the speed of light.

One difficulty in performing the integral is that the signal path through the atmosphere depends on the object properties along the signal path and is, therefore, not a straight line. A change in atmospheric conditions leads to a change in $S$ and integral function $f_s$. Another challenge is related to the distribution of the radio sources and the number of detectors with respect to the size of the object of interest. From single satellite missions, the distribution of integral measurements is not optimal for the reconstruction of three-dimensional structures in an analytical way using the Radon transform. To overcome this limitation, in atmospheric sounding the Abel transform [6], a further simplification of the Radon transform, is generally applied. It allows for the determination of one-dimensional profiles of refractivity from measurements of excess phase, assuming spherical symmetry. In 1965, this technique was applied to measurements of the Mariner four spacecraft to study important properties of the Martian atmosphere and is nowadays commonly applied to GNSS phase measurements obtained from dedicated radio occultation missions [7, 8, 9]. However, standard processing strategies based on the Abel transform do not allow for resolving horizontal features in the atmosphere. With the advent of nanosatellite technology, the number and distribution of signals have significantly increased - leading to the situation that the assumptions made to derive the Abel transform (spherical symmetry and parallel observation paths) become a limiting factor in the analysis of space-based radio occultation observations. To overcome this limitation, the existing observations can be stacked together to solve Eq. (1) either linearly or iterative non-linearly [5]. A complete non-linear solution is difficult to achieve but also not necessary for most applications in geoscience since it can be demonstrated that the signal path is not significantly perturbed by linearization assumptions. In Section 2, common solving techniques (linear and non-linear) are presented, and in Section 3 and Section 4, it is analyzed whether they can be utilized to reconstruct refractivity fields in the neutral atmosphere from GNSS measurements of atmospheric excess phase on-board dense nanosatellite formations.

2. Solving techniques

If $f_s$ in Eq. (1) is replaced by the atmospheric excess phase ($\mathit{aep}$) and $g\left(s\right)$ by the relation $n-1$, the basic function of atmospheric tomography is obtained as follows:

$$\mathit{aep}={\int }_{S}n\cdot \mathit{ds}-{\int }_{S_0}\mathit{ds}$$E3

where $S$ is the “true” signal path and $S_0$ is the theoretical straight line signal path in a vacuum. The second term on the right-hand side of Eq. (3) stems from the definition of $\mathit{aep}$, describing only the atmospheric contribution to the excess phase. Strictly speaking also the limits of the integral have to be adapted to the relevant parts in the atmosphere, e.g. up to about $80\mathit{km}$ altitude for the neutral atmosphere.

Fermat’s principle tells us that first-order changes in the signal path cause second-order changes in signal travel time, i.e. for small variations in the object properties, the travel time is stationary. This principle is very beneficial since it allows to define two simplified versions of atmospheric tomography, the so-called linear and non-linear approach. In linear tomography, the bent signal path $S$ is replaced by a straight line $S_0$ and corrections to $n$ are made by ignoring atmospheric bending. In contrast, the iterative non-linear approach takes the signal bending into account by the definition of the ray paths but not in the inversion of $n$ along $\mathit{ds}$. This means after each processing step the signal paths are re-computed, e.g. by solving the so-called Eikonal using ray-tracing shooting techniques [10].

A numerical solution for Eq. (1) is obtained by discretizing the object of interest, e.g. the neutral atmosphere in area elements (in two-dimensions) or volume elements (in three-dimensions). Further, it is assumed that in each volume element the index of refraction is constant.² In the atmosphere, the index of refraction $n$ is close to 1, thus it can be replaced by refractivity $N=\left(n-1\right)\cdot {10}^6$. With these adaptions Eq. (1) reads:

$$\mathit{aep}=\sum _{k=1}^m{N}_k\cdot d_k$$E4

where $N_k$ is the constant refractivity and $d_k$ is the ray length in the volume element $k$. Assuming $l$ observations, indexed by $j=1,2,\dots ,l$ and $m$ volume elements (short: voxels), indexed by $k=1,2,\dots ,m$, the individual observation equations can be combined into a linear equation system. In matrix notation the resulting tomographic equation reads:

$$\mathit{AEP}=\mathit{A}\cdot \mathit{N}$$E5

where $\mathit{AEP}$ is the observation vector of size $\left(l,1\right)$ and $\mathit{N}$ is the unknown vector of size $\left(m,1\right)$ describing the properties in each volume element $k$. The $\left(l,m\right)$ matrix $\mathit{A}$ contains the spatial derivatives of the observations ${\mathit{aep}}_j$ with respect to the unknowns $N_k$.

$$\mathit{A}=\left[\begin{array}{cccc}\frac{\partial {\mathit{aep}}_1}{\partial N_1}& \frac{\partial {\mathit{aep}}_1}{\partial N_2}& \dots & \frac{\partial {\mathit{aep}}_1}{\partial N_m}\\ \frac{\partial {\mathit{aep}}_2}{\partial N_1}& \frac{\partial {\mathit{aep}}_2}{\partial N_2}& \dots & \frac{\partial {\mathit{aep}}_2}{\partial N_m}\\ ⋮& ⋮& ⋮& ⋮\\ \frac{\partial {\mathit{aep}}_l}{\partial N_1}& \frac{\partial {\mathit{aep}}_l}{\partial N_2}& \dots & \frac{\partial {\mathit{aep}}_l}{\partial N_m}\end{array}\right]$$E6

Since Eq. (4) is linear, the partial derivatives of $\mathit{aep}$ are the ray lengths ($d_k$) in each voxel. For linear tomography, $d_k$ is computed from the line of sight vector between the transmitter and receiver. In the non-linear approach, the bending of the signal is taken into account, e.g. by using ray tracing shooting techniques [13]. Therefore, a priori information about the atmospheric state (e.g. in the form of numerical weather forecasts) is needed. Dependent on the quality of the a priori model, additional iterations might be necessary. After each iteration, the estimated refractivity field is considered as input for the ray tracer. The processing is repeated until the determination of the ray path converges. This happens usually after 1–2 iterations.

2.1 The inverse problem

An analytical solution for the tomographic equation (Eq. (5)) can be found by the inversion of matrix $\mathit{A}$.

$$\mathit{N}={\mathit{A}}^{-1}\cdot \mathit{AEP}$$E7

The inverse ${\mathit{A}}^{-1}$ exists if the determinant of $\mathit{A}$ is non-zero. This requires that $\mathit{A}$ is a squared matrix $\left(l=m\right)$. Otherwise, the matrix $\mathit{A}$ is called singular, i.e. does not have a matrix inverse. Unfortunately latter is the case in most atmospheric tomography applications since the observation data, e.g. the GNSS measurements, are considered as “incomplete”. Therefore, the matrix $\mathit{A}$ has zero singular values and Eq. (7) becomes ill-posed. In literature, several strategies are described, which allow to remove the deficiencies of the ill-posed problem. They either try to solve the inverse problem or to avoid it. The most prominent ones are:

• Iterative algebraic reconstruction techniques

• Truncated singular value decomposition

• Tikhonov regularization

These three techniques have been selected since they were proven in practice as the most reliable, as described briefly in the following subsections.

2.1.1 Algebraic reconstruction techniques

The iterative algebraic reconstruction technique (ART) has been suggested in 1937 by [14] for solving linear equation systems. This technique avoids the inversion problem and initializes the matrix $\mathit{A}$ row-wise. This is very beneficial for large equation systems. Applied to Eq. (7) the ART algorithm reads:

$$N^{i+1}=N^i+\frac{\omega }{⟨A_j,A_j⟩}\cdot \left({\mathit{aep}}_j-⟨A_j^T,N^i⟩\right)\cdot A_j$$E8

where $A_j$ indicates the $j$th row of the matrix $A$, $⟨A_j,A_j⟩$ is the resulting inner product and the difference ${\mathit{aep}}_j-⟨A_j^T,N^i⟩$ is the prefit-residual from the last iteration. Based on the number of traversed volume elements and the relaxation factor $\omega$ the residual is split into multiple components and applied to $N^i$ in order to obtain the improved refractivity field $N^{i+1}$, which is again input for the next iteration. The processing is stopped once Eq. (8) converges to the solution of Eq. (7) with minimal norm if $0<\omega <2$. For ground-based GNSS networks, the best results have been obtained with a relaxation parameter of $\sim 0.175$, see [15]. Studies from [16] or [17] manifest that the algorithms of the ART family are also very successful in reconstructing the total electron content (TEC) in the ionosphere. Dependent on how the discretization is done, different realizations of ART exist. For tomographic reconstruction especially the multiplicative algebraic reconstruction techniques (MART) and the simultaneous iterative reconstruction techniques (SIRT) are worth mentioning, see [18] or [19]. In contrast to the original ART algorithm, MART leads in general to faster convergence and SIRT has the benefit of being impervious to the order of measurements (${\mathit{aep}}_j$).

2.1.2 Truncated singular value decomposition

For ill-conditioned least squares problems, [20, 21] invented a general solution, widely known as pseudo inverse or Moore-Penrose inverse ${\mathit{A}}^{+}$. A numerical solution for the pseudo inverse can be obtained by singular value decomposition [22]. This requires a split of the matrix $\mathit{A}$ into three components as follows:

$$\mathit{A}=\mathit{U}\cdot \mathit{S}\cdot {\mathit{V}}^T$$E9

with $\mathit{U}$ ($l,l$) and ${\mathit{V}}^T$($m,m$) as orthogonal matrices, containing the normalized left and right singular vectors of $\mathit{A}$, respectively. Matrix $\mathit{S}$ ($l,m$) is a diagonal matrix with singular values $s_{k,k}$ arranged in descending order. By using only the non-zero diagonal elements of $\mathit{S}$ the pseudo inverse is obtained as follows:

$${\mathit{A}}^{+}=\mathit{V}\cdot {\mathit{S}}^{-1}\cdot {\mathit{U}}^T$$E10

The 2-norm of the matrix $\mathbf{S}$ defines the condition number $\kappa \left(\mathbf{A}\right)$. It can be interpreted as the ratio between the largest and the smallest singular values.

$$\kappa \left(\mathbf{A}\right)=\frac{\left|s_{\mathit{max}}\right|}{\left|s_{\mathit{min}}\right|}$$E11

A well-conditioned matrix has a condition number $\kappa \left(\mathit{A}\right)$ near 1. The resulting tomography solution is rather insensitive to measurement errors. A large condition number indicates an ill-conditioned problem. According to Eq. (11), the condition number of $\mathit{A}$ improves by neglecting tiny singular values. This technique is known as truncated singular value decomposition (TSVD), see [23]. It allows to approximate the ill-conditioned matrix $\mathit{A}$ by a matrix $\tilde{\mathit{A}}$ of lower rank.

For atmospheric tomography, [24] suggested $s_{\mathit{lim}}=2.8\mathit{km}$ as the threshold for $s_{k,k}$, i.e. all singular values smaller than $s_{\mathit{lim}}$ are set to zero. However, in practice an optimal threshold for $s_{\mathit{lim}}$ can be determined using the L-curve technique [23]. Therefore, a set of solutions is determined with varying $s_{\mathit{lim}}$-values. For each solution, the 2-norm of the estimated parameters is plotted against the 2-norm of the residuals. By connecting all points in a log–log plot a concave L-shaped curve is obtained, whereby the corner of the curve, i.e. the point of maximum curvature, defines the optimal threshold ($s_{\mathit{opt}}$) for singular value decomposition. Figure 1 shows an L-curve for a typical GNSS tomography setting. In this example, the optimal solution is obtained by setting $s_{\mathit{opt}}$ to $0.032$. This point was found by testing various values for $s_{\mathit{lim}}$ between $3\cdot {10}^{-1}$ and $3\cdot {10}^{-5}$. After each processing step, the solution norm $\mathit{log}{‖\mathit{N}‖}_2$ is plotted against the residual norm $\mathit{log}{‖\mathit{A}\cdot \mathit{N}-\mathit{AEP}‖}_2$ and the corner point of the resulting curve is marked as the optimal solution.

Figure 1.

Representative L-curve for a typical GNSS tomography inversion problem. The red dot indicates the corner point of the L-curve and therewith the optimal tomography solution.

2.1.3 Tikhonov regularization

A more generalized solution to the regularization problem can be found in [25], who describes the minimization problem as follows:

$${\mathit{N}}_{\eta }=\mathit{arg}\mathit{min}\left\{{‖\mathit{A}\cdot \mathit{N}-\mathit{AEP}‖}_2^2+{\eta }^2{‖\mathit{L}\left({\mathit{N}}_0-\mathit{N}\right)‖}_2^2\right\}$$E12

where $\eta$ is called the regularization parameter or Tikhonov factor and ${\mathit{N}}_0$ is an approximation of $\mathit{N}$. The “size” of the solution is defined by the norm ${‖\mathit{L}\left({\mathit{N}}_0-\mathit{N}\right)‖}_2$ and the “fit” by the norm of the residual vector ${‖\mathit{A}\cdot \mathit{N}-\mathit{AEP}‖}_2$. One possibility to solve Eq. 12 is to treat it as a least squares problem. In [26] it is shown that the matrix $\mathit{L}$ can be replaced by an identity matrix $\mathit{I}$, i.e. the condition number of $\mathit{A}$ is improved by adding a small multiple of the identity to the matrix $\mathit{A}$.

$$\tilde{\mathit{A}}=\mathit{A}+\eta \cdot \mathit{I}$$E13

A possible solution for $\eta$ can be obtained by means of singular value decomposition (see Section 2.1.2). Thereby the elements of the diagonal matrix $\mathit{S}$ are replaced by the coefficients $r_{k.k}$.

$$r_{k,k}=\frac{s_{k,k}^2}{s_{k,k}^2+\eta }$$E14

If the Tikhonov factor is defined as a sharp filter

$$\eta =\left\{\begin{array}{c}1\mathit{\text{for}}{s}_{k,k}\ge s_{\mathit{lim}}\\ 0\mathit{\text{for}}{s}_{k,k}<s_{\mathit{lim}}\end{array}\right.$$E15

the resulting solution can be interpreted as smoothed TSVD solution.

2.2 The partial least squares solution

By treating Eq. (7) as a least squares problem, the basic equation of the weighted least squares estimator reads:

$$\widehat{\mathit{N}}={\left({\mathit{A}}^T\cdot \mathit{P}\cdot \mathit{A}\right)}^{-1}\cdot {\mathit{A}}^T\cdot \mathit{P}\cdot \mathit{AEP}$$E16

where $\mathit{P}$ is a weighting matrix, which allows to take the relative accuracy and possible constraints between the observations into account. The least squares solution $\widehat{\mathit{N}}$ is obtained by minimizing the 2-norm of the observation residuals. Thereby we assume, that the observations are normally distributed, i.e. free of gross errors or systematic effects.

By combining Eq. (10) with Eq. (16) the tomography solution reads:

$$\widehat{\mathit{N}}=\mathit{V}\cdot {\mathit{S}}^{-1}\cdot {\mathit{U}}^T\cdot {\mathit{A}}^T\cdot \mathit{P}\cdot \mathit{AEP}$$E17

where the columns of $\mathit{U}$ and ${\mathit{V}}^T$ are the normalized left and right singular vectors of ${\mathit{A}}^T\cdot \mathit{P}\cdot \mathit{A}$, respectively and matrix $\mathit{S}$ ($l,m$) is the diagonal singular value matrix as defined in Section 2.1.2.

2.2.1 The a priori field

In Section 2, the linear and non-linear approaches have been defined to reconstruct the GNSS signal paths through the atmosphere. While the linear approach is not dependent on any external data, the non-linear approach requires an a priori refractivity field, e.g. derived from the standard atmosphere or numerical weather forecasts, to reconstruct the bent signal path. Besides, the a priori field can be also utilized to stabilize the tomographic equation system. One possibility is to treat the additional information (${\mathit{N}}_0$) as absolute constraints:

$$\widehat{\mathit{N}}={\mathit{N}}_0+\mathit{V}\cdot {\mathit{S}}^{-1}\cdot {\mathit{U}}^T\cdot {\mathit{A}}^T\cdot \mathit{P}\cdot \left(\mathit{AEP}-\mathit{A}\cdot {\mathit{N}}_0\right)$$E18

In the following, this solution is called the constrained solution.

Another possibility to handle the extended equation system is to treat it as a system of subsets with

$${\mathit{A}}_{\mathit{ext}}=\left[\begin{array}{c}\mathit{A}\\ {\mathit{A}}_c\end{array}\right]$$E19

$${\mathit{AEP}}_{\mathit{ext}}=\left[\begin{array}{c}\mathit{AEP}\\ {\mathit{N}}_0\end{array}\right]$$E20

$${\mathit{P}}_{\mathit{ext}}=\left[\begin{array}{c}\mathit{P}\\ {\mathit{P}}_c\end{array}\right]$$E21

where ${\mathit{A}}_c$ is the design matrix and ${\mathit{P}}_c$ the weighting matrix for ${\mathit{N}}_0$. The extended equation system can be solved using Eq. (17), whereby $\mathit{A}$, $\mathit{AEP}$ and $\mathit{P}$ are replaced by its extended complements ${\mathit{A}}_{\mathit{ext}}$, ${\mathit{AEP}}_{\mathit{ext}}$ and ${\mathit{P}}_{\mathit{ext}}$. In principle, it provides the same results as the constrained solution.

A third possibility would be to solve Eq. (18) separately for each observation type using the estimates ($\widehat{\mathit{N}}$) and the variance–covariance matrix of the estimates (${\mathit{Cov}}_{\widehat{N}N}$) from the first step as a priori information for the next step. In the case of two subsets the corresponding tomography solution reads:

$${\widehat{\mathit{N}}}_1=\mathit{V}\cdot {\mathit{S}}^{-1}\cdot {\mathit{U}}^T\cdot {\mathit{A}}^T\cdot \mathit{P}\cdot \mathit{AEP}$$E22

$${\mathit{Cov}}_{\widehat{N}N}=\mathit{V}\cdot {\mathit{S}}^{-1}\cdot {\mathit{U}}^T$$E23

where $\mathit{U}$, $\mathit{V}$ and $\mathit{S}$ are obtained by singular value decomposition of the matrix ${\mathit{A}}^T\cdot \mathit{P}\cdot \mathit{A}$. For the second (final) solution both, ${\widehat{\mathit{N}}}_1$ and ${\mathit{Cov}}_{\widehat{N}N}$ are introduced into the equation system as follows:

$$\widehat{\mathit{N}}={\widehat{\mathit{N}}}_1+\mathit{V}\cdot {\mathit{S}}^{-1}\cdot {\mathit{U}}^T\cdot {\mathit{A}}_c^T\cdot {\mathit{P}}_c\cdot \left({\mathit{N}}_0-{\mathit{A}}_0\cdot {\widehat{\mathit{N}}}_1\right)$$E24

with $\mathit{S}$, $\mathit{U}$ and $\mathit{V}$ obtained by truncated singular value decomposition (see subsection 2.1.2) of the matrix ${\mathit{A}}_0^T\cdot {\mathit{P}}_0\cdot {\mathit{A}}_0+{\mathit{Cov}}_{\widehat{X}X}^{-1}.$ In the following this solution is called the partial solution. In case the matrix $\mathit{A}$ is of full rank or if only one set of observations is available, the constrained solution and the partial solution provide identical results. In the case of an ill-conditioned matrix, the partial solution has the advantage that the eigenvalue can be computed for each subset of observations. In large equation systems, this allows to reduce computational load since the matrix $\mathit{A}$ is divided into several parts.

2.2.2 Observation weights

Up to now, the individual observations were considered as uncorrelated and equally accurate. However, for varying input data it might be beneficial to set up a weighting matrix. In case the relative accuracy between observations is known, they can be directly introduced into the equation system by defining the weighting matrix

$$\mathit{P}={\sigma }_0^2\cdot {\mathit{Cov}}_{\mathit{ll}}^{-1}$$E25

where variance co-variance matrix ${\mathit{Cov}}_{\mathit{ll}}$ reflects the precision of the observations on its diagonal elements (${\sigma }_n^2$) with ${\sigma }_0^2$ as the a priori variance of the unit weight. In the case of uncorrelated observations and unit variances, the matrix ${\mathit{Cov}}_{\mathit{ll}}$ simplifies to an identity matrix of size $\left(n,n\right)$.

In case no accurate information is available, a weighting model can be utilized. For ground-based GNSS observations, an elevation-dependent weighting is common, for satellite-to-satellite observations a weighting based on carrier-to-noise density $C/N0$ seems to be more useful since the observations are usually gathered around $0deg$ the elevation angle or below. For the a priori refractivity field, a height-dependent weighting model is recommended. By focusing on the neutral atmosphere and assuming that Eq. (26) is exact

$$N=K_1\cdot \frac{p_d}{T}+K_2\cdot \frac{e}{T}+K_3\cdot \frac{e}{T^2}$$E26

the theoretical standard deviation for refractivity reads³:

$${\sigma }_N={\left[{\left(\frac{\partial N}{\partial T}\cdot {\sigma }_T\right)}^2+{\left(\frac{\partial N}{\partial p}\cdot {\sigma }_p\right)}^2+{\left(\frac{\partial N}{\partial e}\cdot {\sigma }_e\right)}^2+\dots \right]}^{\frac{1}{2}}$$E27

For in-situ measurements, [27] provides theoretical standard deviations for pressure $p$, temperature $T$, and water vapor pressure⁴ $e$ for a wide range of meteorological sensors. In addition, height-dependent error curves for the three meteorological parameters can be obtained from [28]. The resulting uncertainties, assuming standard atmospheric conditions at sea level, are listed in Table 1.

$\frac{\partial N}{\partial T}\mathrm{\Delta }T$	$\frac{\partial N}{\partial p}\mathrm{\Delta }p$	$\frac{\partial N}{\partial e}\mathrm{\Delta }e$	$\mathrm{\Delta }N$
$\pm 0.25\mathit{ppm}$	$\pm 0.08\mathit{ppm}$	$\pm 2.37\mathit{ppm}$	$\pm 2.39\mathit{ppm}$

Table 1.

Standard deviation of refractivity and its components, assuming a typical meteo sensor error of $0.3\mathit{hPa}$ for pressure, $0.2K$ for temperature and $3\%$ for relative humidity - computed for standard atmospheric conditions at sea level ($p=1013\mathit{hPa},T={15}^{\circ }C,\mathit{rh}=60\%$).

The standard deviation of refractivity is $2.39\mathit{ppm}$, which equates to a relative uncertainty of $0.75\%$. By far the largest impact ($2.37\mathit{ppm}$) is related to the uncertainty of water vapor. Consequently, the utmost care has to be taken when measuring humidity and temperature.

3. Observations of atmospheric excess phase

Satellite refractometric sounding of the atmosphere and the underlying inverse problems have been under investigation since the 1960s, see [29, 30, 31]. However, it was not until 1976 that the first radio occultation (RO) experiment was carried out to survey the earth‘s atmosphere within the Apollo-Soyuz mission [32]. Until then, the major problem noted was the lack in accuracy of refractometric measurements of phase or Doppler shift [33]. This limitation has widely been overcome with the emergence of the Global Positioning System (GPS) around the 1980s [7]. Since the proof of concept during the GPS-MET satellite mission in 1995 various satellites have been equipped with precise GNSS radio occultation receivers, leading to approximately $500-600$ globally distributed radio occultation profiles per day and satellite, assuming a 32-satellite GPS constellation.

3.1 The observation equation

The phase that a receiver obtains from a GNSS satellite can be modeled as

$$L_{r,\nu }^s={\mathrm{\varrho }}_r^s+c\cdot \delta t_r-c\cdot \delta t^s+\mathrm{\Delta }{\mathrm{\varrho }}_{r,v}^s+\nu \cdot n_{r,v}^s$$E28

with

$L_{r,\nu }^s$	... phase observation for the transmitter-receiver pair (s-r) at frequency $\nu$
${\mathrm{\varrho }}_r^s$	... geometric distance between transmitter s and receiver r
$c\cdot \delta t_r$	... correction of receiver clock at signal reception time $t$
$c\cdot \delta t^s$	... correction of satellite clock at transmission time $t-{\tau }_r^s$
$\mathrm{\Delta }{\mathrm{\varrho }}_{r,v}^s$	... all delays due to propagation effects
$n_{r,v}^s$	... unknown integer number of cycles (carrier phase ambiguity)

The objective of refractometric sounding is to extract the atmospheric propagation effects from the phase observations. Assuming that relativistic effects, satellite-specific multipath effects and antenna-specific phase center corrections are known and removed, the remaining effects in $\mathrm{\Delta }{\mathrm{\varrho }}_{r,v}^s$ can be divided into two terms:

$$\mathrm{\Delta }{\mathrm{\varrho }}_{r,v}^s=\mathrm{\Delta }{\mathrm{\varrho }}_{r,\mathit{trp}}^s+K\frac{{\mathit{TEC}}_{r,v}^s}{f_r^2}$$E29

where the first term ($\mathrm{\Delta }{\mathrm{\varrho }}_{r,\mathit{trp}}^s$) describes the delay of the carrier phase in the neutral atmosphere and the second term $\left(K\frac{{\mathit{TEC}}_{r,v}^s}{f_r^2}\right)$ the advancement of the carrier phase in the ionosphere. The integral term ${\mathit{TEC}}_{r,v}^s$ is the electron density along the ray path between transmitter s and receiver r, scaled by a constant term $K$.

A detailed description of the individual systematic effects can be found in [34]. In the following, special attention is given to the modeling and estimation of neutral atmospheric effects ($\mathrm{\Delta }{\mathrm{\varrho }}_{r,\mathit{trp}}^s$) assuming that the first-order ionospheric effect (up to $99.9\%$) can be removed by forming an ionospheric-free linear combination $L_{\mathit{\text{IF}}}$. Condition therefore is, that the receiver tracks the GNSS carrier phase simultaneous on two frequencies

$$L_{\mathit{\text{IF}}}=\frac{f_1^2\cdot L_{r,1}^s-f_2^2\cdot L_{r,2}^s}{f_1^2-f_1^2}$$E30

where the nominal frequencies $f_1$ and $f_2$ are defined by the satellite system frequency plan (e.g. $1575.42\mathit{MHz}$ for GPS $L1$ and $1227.60\mathit{MHz}$ for GPS $L2$).

3.2 Calibration of the phase signal

For the extraction of atmospheric phase excess from phase observations, first, the phase signal has to be calibrated. Therefore, the clock effects in Eq. (28) are eliminated. This can be achieved if the occulting receiver satellite can simultaneously see an occulting GNSS and a non-occulting GNSS satellite. Elimination of the clock effects is also possible if the same GNSS satellites are visible from a ground-based GNSS receiver. In addition, one needs to precisely know the position and velocity of the transmitting and receiving satellites. The orbits for the GNSS satellites can be obtained from services such as the International GNSS service (see www.igs.org). The position and velocity of the receiving satellite have to be computed using the phase measurements recorded from the GNSS radio occultation receiver or from a dedicated POD (precise orbit determination) receiver on board the receiving satellite. For further details about the procedure of POD, the reader is referred to [35]. A more detailed description of the calibration of the phase measurements is given in, e.g. [36]. Once the signal is calibrated, the atmospheric phase excess can be used to set up the observation equations for the tomographic processing. The precondition for a stable tomography solution is, that enough overlapping observations are available. Therefore, in the following the concept of a dense nanosatellite formation is introduced.

4. Concept of a dense nanosatellite formation

4.1 Introduction

In recent years, nanosatellite technology has become increasingly important for a wide variety of applications, such as communication, technology demonstration, heliophysics, astrophysics, earth science, or planetary science [37, 38, 39]. Most of the existing mission concepts are based on the CubeSat form factor established by Professor Bob Twiggs at the Department of Aeronautics and Astronautics at Stanford University in late 1999. Although small satellites have existed since the very beginning of spaceflight, the definition of the CubeSat standard “made it possible to bring production to a level of flexibility and innovation never seen before” [40]. As of September 2022, over 2000 nanosatellites were launched into orbit, with a record number of 143 satellites launched on a single rocket on board the Transporter-1 mission in January 2021. In the next decade, we expect that the number of nanosatellite launches per year will continue to rise by a factor of 4–5, leading to dense observation networks in low earth orbit. Innovations in satellite technology, such as miniaturized GNSS receivers [41] and intelligent processing strategies will further boost the realization of new observation concepts based on nanosatellite technology and the establishment of dense satellite formations in highly interesting but yet scarcely-explored regions in the earth‘s atmosphere and beyond.

4.2 The observation geometry

The multi-frequency signals from over a hundred active GNSS satellites gathered on board each nanosatellite allow for measuring the atmospheric state with unprecedented spatiotemporal resolution. For the proof of concept, we assume a formation of four nanosatellites, injected into a polar orbit. The advantage from such a configuration is that we can get simultaneous radio occultation observations that are closely located [42]. Figure 2 shows the observation geometry together with the ray paths through the lower $8\mathit{km}$ of the atmosphere.

Figure 2.

Left: The observation geometry for one GNSS satellite simultaneous observed by four nanosatellites in a string-of-pearls formation. Right: The resulting radio occultation signal paths.

The spacing between the nanosatellites is set to $\mathit{dM}=1.9deg$ (approx. $230\mathit{km}$). At an altitude of $550\mathit{km}$ this corresponds to a temporal spacing of about $30s$. Due to the limb-sounding geometry and high inclination of most nanosatellites, a global distribution can be obtained with $500-600$ radio occultation events per nanosatellite and day assuming a 32-satellite GPS constellation.

The angle under which the GNSS satellites are observed is constantly changing. In order to characterize the observation geometry, we distinguish between two scenarios. In the first scenario, the observation angle is close to $90deg$, i.e. the RO measurements are obtained in a cross-track direction. In consequence, from the four nanosatellites, we obtain ray paths that are widely parallel to each other. In the second scenario, the angle is close to $0deg$ or $180deg$. This leads to RO measurements in the flight direction or anti-flight direction. In both cases, a unique observation geometry is obtained, in which consecutive observations overlap, as shown in Figure 2.

4.3 Tomography case study

At the time of writing, real GNSS measurements from a dense nanosatellite formation were not available. Thus, for technique demonstration, a closed-loop simulation was carried out using the Weather Research and Forecasting (WRF) model to simulate the atmospheric state along the GNSS radio occultation signals shown in Figure 2.

In the first step, the signal paths through the atmosphere were reconstructed every $500\mathit{ms}$ using ray-tracing shooting techniques [13] with a step size of $1\mathit{km}$. For each ray point, wet refractivity ($N_w$) was computed from WRF temperature ($T$) and water vapor pressure ($e$) fields using Eq. (31)

$$N_w=K_2^{\prime }\cdot \frac{e}{T}+K_3\cdot \frac{e}{T^2}$$E31

where the constant $K_2^{\prime }$ is given by

$$K_2^{\prime }=K_2-K_1\cdot \frac{M_w}{M_d}$$E32

with $K_1=77.689\frac{K}{\mathit{hPa}}$, $K_2=71.2952\frac{K}{\mathit{hPa}}$ and $K_3=375463\frac{K^2}{\mathit{hPa}}$ as the refractivity constants and $M_w$ and $M_d$ as the molar mass of dry air and water vapor, respectively.

Figure 3 (left) shows the resulting wet refractivity distribution in the lowest $8\mathit{km}$ of the atmosphere, with a water vapor inversion layer at a height of approximately $2-4\mathit{km}$. By integration along the signal paths, the wet refractivity can be converted into atmospheric excess phase using Eq. (4). Figure 3 (right) shows that the inversion layer in the WRF model also propagates into the simulated observations of atmospheric phase excess.

Figure 3.

Left: Weather Research and Forecasting (WRF) model derived wet refractivity fields [ppm] with the overlaying white lines showing the tangent points of the four RO ray paths through the lower $8\mathit{km}$ of the atmosphere. Right: The resulting atmospheric excess phase observations [m] by integration through the wet refractivity field.

To reconstruct the 2D refractivity fields from the atmospheric excess phase, the area covered by the observations was discretized in area elements with a grid size of 22 km (horizontally) and 0.2 km (vertically). The tomographic processing itself was carried out with the ATom software package [13]. Table 2 summarizes the major settings.

Parameter	Settings
Case study domain	Equatorial pacific ocean ($140-150degE$)
Case study period	Late autumn 2006
Model resolution	$22\mathit{km}$ (horizontally) $\times$ $0.2\mathit{km}$ (vertically)
Tomography software	Modified version of ATom software package^*
Initial field	smooth WRF field
Inversion method	Singular value decomposition (${\mathit{\text{eigenv}}}_{\mathit{min}}=0.01{\mathit{km}}^2$)
Estimation method	Iterative weighted least squares adjustment
Convergence criteria	RMS of weighted residuals

Table 2.

Tomography settings applied for the reconstruction of refractivity fields from (simulated) RO observations.

^*

https://github.com/GregorMoeller/ATom.

The resulting refractivity fields are visualized in Figure 4. The upper left plot shows the a priori field, a smooth WRF refractivity field, used to initialize the tomography solution. For the computation of the smoothed field, a sliding window filter was applied to the WRF data to remove the inversion layer and therefore, reduce the information contained in the initial field. In the upper right plot, the actual tomography solution is shown. By comparison with the WRF reference field (lower left plot), the reconstruction capabilities of the tomography approach can be assessed. The differences between the two models are shown in the lower right plot. Overall voxels, a Root Mean Square Error (RMSE) of $0.9\mathit{ppm}$ ($21.8\%$) and a bias of $0.03\mathit{ppm}$ was received.

Figure 4.

Top left: Smooth WRF refractivity field used to initialize the tomography solution. Top right: Estimated refractivity field (tomography solution). Bottom left: WRF refractivity field (reference). Bottom right: Closed-loop validation (tomography minus WRF) to assess the performance of the tomography approach.

Overall, the best solution is obtained within the horizontal range (−250 km, 250 km) in which multiple observations overlap and therefore, help to stabilize the tomography resolution. In this core domain of the tomography model, an RMSE of $0.5\mathit{ppm}$ ($9.6\%$) was obtained, which is by a factor of two better than in the outer regions.

5. Conclusions and outlook

In this chapter, the basic aspects of the remote sensing of the lower earth’s atmosphere using tomography radio occultation methods are addressed. My motivation was to provide an overview about the current achievements in tomographic processing and its potential for the processing of radio occultation measurements collected from very light-weight and power-efficient GNSS sensors onboard dense nanosatellite formations. In a number of closed-loop validations, the expected observations have been analyzed and possible processing strategies have been evaluated. Due to the unique observation geometry, combined processing of overlapping radio occultation measurements using tomographic principles is possible and allows to generate high-resolution cross-sections of the lower atmosphere. Thus, I believe that tomography products have great potential to advance current knowledge, e.g. as a weather analysis tool or as a complementary observation technique for water vapor distribution, which can be assimilated into operational weather forecast systems. Once the required sensor technology is available, not only the communication industry but also the earth observation community will benefit from new observation concepts based on nanosatellite technology. If the proposed observation concept is also suited for the monitoring of the ionosphere has to be evaluated in future studies.

Conflict of interest

The authors declare no conflict of interest.

Abbreviations

AEP	Atmospheric Excess Phase
ART	Algebraic Reconstruction Technique
ATom	Atmospheric Tomography software package
GNSS	Global Navigation Satellite Systems
MART	Multiplicative Algebraic Reconstruction Technique
POD	Precise Orbit Determination
RMSE	Root Mean Square Error
RO	Radio Occultation
SIRT	Simultaneous Iterative Reconstruction Technique
TEC	Total Electron Content
TSVD	Truncated Singular Value Decomposition
WRF	Weather Research and Forecasting model