Estimation of electron densities in the lower thermosphere from GUVI 135.6 nm tomographic inversions in support of SpreadFEx

. The SpreadFEx campaign was conducted with the goal of investigating potential neutral atmospheric dynamics inﬂuences in seeding plasma instabilities and bubbles extending to higher altitudes from September to November 2005, with primary measurements in Brazil. In this paper, we present the results of space-based UV and ground-based optical observations in support of this campaign. Speciﬁ-cally, we present multi-dimensional electron density images obtained tomographically from the 135.6 nm emissions measured by the GUVI instrument aboard the TIMED satellite that result from radiative recombination of O + and compare those with the corresponding 630.0 nm OI images recorded in the Brazilian sector. The GUVI results provide altitude vs. longitude information on depleted regions in the ionospheric plasma density that are complementary to the single-height latitude-longitude images obtained with the airglow imager.


Introduction
Irregularities resulting from plasma turbulence in the nighttime equatorial ionosphere are of considerable importance because scintillations caused by density irregularities (commonly referred to as "Equatorial Spread F ", or ESF) can result in outages of the communication and navigation systems that depend on trans-ionospheric radio links. The plasma instabilities associated with spread F have been observed since the advent of ionosondes and studied systematically since measurements were compiled at the Jicamarca Radio Observatory in Peru (Farley et al., 1970).
Correspondence to: F. Kamalabadi (farzadk@uiuc.edu) Fully developed ESF events termed Equatorial Plasma Bubbles (EPBs) refer to regions of depleted plasma density that typically originate in the bottomside post-sunset ionosphere and, while longitudinally thin, extend latitudinally along magnetic field lines. EPBs can extend vertically to altitudes above 1000 km (Kelley, 1989). This phenomenon is believed to be generated due to a variety of plasma instability processes, with the Rayleigh-Taylor instability being the primary mechanism at work. The resulting plasma irregularities manifest density depletions with latitudinal scale sizes up to several hundred kilometers. The prediction of ionospheric plasma bubbles poses challenges on the conventional observational and modeling capabilities since under seemingly identical ionospheric conditions they may occur on one day and be absent on another.
Sensing with incoherent scatter radars (ISR) is generally accepted to be the most comprehensive method of measuring plasma properties of the upper atmosphere, and hence has been the primary source of much of what has been learned about ESF. Woodman and LaHoz (1976), for example, created some of the first images of these plasma irregularities using range-time intensity plots with the Jicamarca ISR. Since that time, equatorial plasma bubbles have been observed with ground-based instruments including airglow cameras and space-borne remote sensing and in-situ observations (e.g., Kelley, 1989, and references therein).
Since ISRs generally yield only line-of-sight observations of plasma properties, they are unable to obtain a simultaneous, comprehensive view of a two-dimensional ionospheric slice (i.e., latitude vs. longitude). Airglow cameras have provided a means of effectively imaging the development and structure of EPBs from the ground. By observing an optical emission as a proxy for electron density, these cameras can provide high-resolution movies of plasma bubble formation and drift over a small area of the globe. However, the information obtained from airglow cameras is also confined to a   fixed altitude, i.e., the altitude of the source emission layer used as the density proxy. Furthermore, they are limited in their ability to quantitatively measure the depleted electron density. Nevertheless, such ground-based observations can be coordinated with space-based observations, e.g., Kelley et al. (2003), to provide complementary information on the altitude profile and density.
Tomographic imaging systems have been developed and utilized over the past decade in order to enhance observational methods for studying ionospheric processes with multi-dimensional imaging capabilities. Examples of some common techniques and some experiments are given in Bernhardt et al. (1998) and Kamalabadi et al. (2002).
With the goal of investigating potential neutral atmospheric dynamics influences in seeding equatorial spread F and plasma bubbles extending to higher altitudes, an extensive experimental campaign was performed from September to November 2005, with primary measurements in Brazil. The campaign termed the Spread F Experiment, or Spread-FEx, consisted of coordinated observations focusing on the Brazilian sector and included ground-based optical, radar, digisonde, and GPS measurements at a number of sites. Related data on plasma bubble structures were also collected by the GUVI instrument aboard the TIMED satellite.
In this paper, we present results from tomographic reconstructions of plasma bubbles from TIMED/GUVI measurements and discuss comparisons with the corresponding ground-based optical measurements. We begin by describing the space-based UV and ground-based optical observations. We then proceed to formulate the forward model relating GUVI observations to ionospheric electron densities and the inversion technique for the estimation of multi-dimensional images of plasma bubbles. Finally, we present the reconstructed results in comparison to those obtained from the ground-based red-line imager.

Observations
The observations that are the focus of this paper include UV brightness measurements from the GUVI instrument and ground-based imagers that were deployed in Brazil as part of this campaign.

Space-based UV observations: GUVI instrument
The Global Ultraviolet Imager (GUVI) (Paxton et al., 2004;Christensen et al., 2003) is one of four instruments that constitute the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) science payload, which was launched in December 2001. GUVI is a far-ultraviolet, scanning imaging spectrograph that provides horizon-tohorizon images in five selectable wavelength intervals. These wavelength intervals nominally correspond to H + 121.6 nm, O + 130.4 nm, O + 135.6 nm, and N 2 Lyman-Birge-Hopfield bands 140 to 150 nm and 165 to 180 nm. The instrument consists of a scan mirror feeding a parabolic telescope and Rowland circle spectrograph, with a wedge-and-strip detector at the focal plane (Paxton et al., 1999).
TIMED orbits at 625 km with a 74 • inclination. The GUVI instrument has a 140 • cross-track scan perpendicular to the orbit plane, as illustrated in Fig. 1. The TIMED orbit allows global coverage of the ionosphere with each orbit occurring with a local time about one minute earlier than the previous one, covering all local times every 60 days. The sensitivity of the GUVI instrument at 135.6 nm is approximately 0.5 counts/s/Rayleigh/pixel (Humm et al., 1998). The detector array has 14 spatial elements and 160 spectral elements. The instrument has 14 spatial pixels giving an 11.8 • field of view along the satellite track. There is a 127.2 • cross-track scan for disk images giving 159 cross-track samples of each of the 14 along-track pixels. Each of these pixels has an integration time of 0.064 s. There are typically 389 scans in a given orbit, with each scan taking 15 s. The tangent point is at an altitude of 519 km at a viewing angle of 80.0 • from the nadir, and the tangent point altitude is 152 km at a viewing angle of 68.8 • above the nadir. Table 1 summarized GUVI's specifications.

Ground-based airglow imager
The 630.0 nm imager was located at a site outside Brasilia, with coordinates (14.8 • S, 47.6 • W). The images span a range of 1000 km by 1000 km. Assuming an emission layer at an altitude of 250 km, the field of view ranges from 10.4 • S to 19.1 • S and from 43.1 • W to 52.1 • W. The field of view of the Brasilia imager is adjacent to the magnetic equator, and the altitude range of the emission is greater than the altitude variation along magnetic field lines over the instrument field of view. Therefore, only the longitudinal structure of the GUVI and ground-based images are comparable. The altitude profile from GUVI and the latitude profile from the ground-based imager are complementary datasets that enable imaging of the depleted regions in three dimensions.

Observation model
The UV measurements of primary concern to ionospheric tomography are passive emissions generated by recombination of oxygen ions with electrons and transition to a lower energy state. Radiative recombination of O + in the nighttime ionosphere is the principal source mechanism for a variety of radiations in the UV as well as visible and infrared wavelengths: In the UV wavelength range, the resulting excited O atom emits radiation at 135.6 and 130.4 nm, and as the recombination of O + ions to the 3 P ground state of O takes place, a narrow continuum shortward of 91.1 nm is produced. The 135.6 nm (2p 4 3 P −3s 5 S) is a prominent feature of the airglow and consists of a doublet ( 3 P 2 − 5 S 2 ) and ( 3 P 1 − 5 S 2 ) (Meier, 1991). Since the 135.6 nm emission is optically thin, its observation for a given line of sight can be modeled as an integral along that line and thus provide for a tomographic formulation. Ignoring contributions from neutralization emissions, absorption, and scattering does introduce some error into the forward model. The error from these contributions, however, is on the order of 10% or less of the observed brightness and results in a cumulative error in the reconstructed electron densities on the order of 10 4 cm −3 , and is considered negligible compared to the reconstruction error from photon noise, which is on the order of 10 5 cm −3 .
For observations made on a planar (horizontal/vertical) geometry with the spacecraft at the origin of the coordinate sys- tem, the intensity of the emission I (in Rayleighs) at a view angle θ is given by where s is the the position along the line of sight, and V denotes the volume emission rate, which can be expanded as follows: where α is the temperature-dependent recombination coefficient, n e is the electron density, and n o + is the O + density. The radiative recombination coefficient α 135.6 was taken to be 7.3×10 −13 cm 3 s −1 (Melendez-Alvira et al., 1999). This value of α 135.6 assumes a constant temperature of 1156 K, which is a reasonable approximation for the altitudes around 200 km where this emission is dominantly produced. In the F-region, [O + ] is approximately equal to [e], and thus the emission brightness is approximately equal to the line integral of the electron density squared times the recombination coefficient: Note that these integral equations make the assumption that the field of view can be approximated by a line. This description of the emission process as a line integral reduces the observed brightness to a linear function of n 2 e . Observation of the brightness of these emissions using a space-borne spectrograph can therefore provide the means to invert the integral equations above and obtain F-region volume emission rates and consequently electron densities.
Equation (4) relates the observed brightness as a function of θ at a fixed spacecraft position to an electron density field. With the satellite's motion, observations corresponding to different spacecraft positions may be collected together over a common volume to allow for a tomographic inversion.  A discrete linear observation model for the GUVI measurements of the 135.6 nm emission can be subsequently formulated by collecting multiple line-of-sight observations. The field of view of the instrument for each pixel is approximated as a single line of sight. Values from the 14 spatial pixels of the GUVI spectrograph are averaged into a single value, which results in a slight loss of spatial resolution but serves to simplify the model and enhance the counting statistics. If all pixels contribute equally to the signal, this averaging will enhance the counting statistics by a factor of √ 14. It is also assumed that the 135.6 nm emission is solely due to the radiative recombination of oxygen atoms. Alternate sources of this emission are considered to be negligible, accounting for 10 percent or less of the emission (Dymond et al., 1997).
Ionospheric structures due to plasma instabilities are generally field-aligned. The forward model uses an offset, tilted dipole magnetic coordinate system to appropriately model field-aligned instabilities. Inversion of the forward model will be greatly facilitated by assuming that the ionosphere is constant along magnetic field lines for a 10-degree segment of latitude (equivalent to 11 scans along the GUVI orbit track). Care must be taken in choosing a latitudinal segment of the ionosphere for the inversion that satisfies this requirement.

Discrete model
The three-dimensional ionosphere can be modeled as a series of two-dimensional slices in the plane perpendicular to the orbit track. Given the 74 • inclination of the orbit, these slices can be modeled as having constant latitude and varying only in altitude and longitude. A two-dimensional slice of the ionosphere can then be divided into discrete sections, with each section having constant n e . If the sections of con-stant squared electron density are arranged into a single vector x, a series of observations y can then be modeled through the following matrix equation.
where A has elements that are proportional to the length of the line of sight for the observation in each section of x and q is an additive noise term. Figure 2 shows a 3-D perspective of the discrete observation model and the resulting 2-D tomographic geometry. The measurement vector y is constructed directly from GUVI data. The 135.6 nm color is selected, and then data from each of the 14 spatial pixels are averaged together. A 10 • section of the orbit is selected, corresponding to 11 crosstrack scans of the GUVI instrument. The 159 pixels from the cross-track are then binned into 79 pixels, with the 79 pixels corresponding to evenly spaced zenith angles from 120 • to 240 • . This binning again enhances the counting statistics for the measurement and ensures that the number of measurements that compose the y vector will be roughly equal to the number of unknowns in the x vector. The resulting 2-D array can be collapsed into the 879-element y vector.
The x vector is a 754-element vector that is a collapsed two-dimensional grid of 29 altitudes and 26 longitudes. The 29 altitudes are evenly spaced from 90 km at 20 km intervals. The 26 longitudes are evenly spaced, with a range of ±5 • longitude at 0.4 • intervals, relative to the satellite longitude in the middle of the orbit section.
The additive noise q is dominated by photon counting noise. A typical integration will yield only a few counts. Statistically, this noise behaves according to a Poisson distribution, resulting in a signal-to-noise ratio (SNR, defined as y q ) equal to the square root of the number of counts. With the GUVI sensitivity of 0.5 counts/s/Rayleigh/pixel, pixel integration time of 0.064 s, and 28-pixel binning for each element in y, the estimated SNR for a 40 R source is 6 (7.78 dB). Typical 135.6 nm brightness in the equatorial anomaly is over 100 R, (Sagawa et al., 2003) yielding a SNR of approximately 10 (10 dB).
Since the ionosphere is assumed to be constant for the 11 scans used to create the y vector, calculation of the A matrix can be reduced to a two-dimensional problem by approximating the direction of the scan as being purely along the East-West direction. The longitudinal variation of the position of the TIMED satellite is retained, effectively yielding overlapping measurements from a moving sensor. The combined effect of these approximations casts the problem as a two-dimensional limited-angle tomography problem, with geometry illustrated in Fig. 2. Each row of the A matrix corresponds to one observation in the y vector. Each element in the row is the length of the portion of the line of sight that is contained by the corresponding element in the x vector. The algorithm for constructing the A matrix determines the value of these lengths for each line of sight, thereby constructing the A matrix one row at a time.

Inversion technique
Rank-deficient and ill-conditioned matrices are typical of a limited-angle tomography problem. Direct least-squares inversion is not adequate for such problems (Kamalabadi et al., 1999). Any successful inversion technique must impose additional constraints on the solution set while effectively minimizing the influence of noise.
The use of physics-based constraints can provide additional information to the brightness measurements and enhance the reconstruction. These constraints take the form of cost functionals to be minimized along with the standard least-squares cost function. The justification for the choice of the constraints are discussed in detail in Comberiate et al. (2006). This section summarizes the development of cost functionals to ensure smoothness and edge-preservation. An optimization technique is formulated for minimizing such cost functionals (Delaney and Bresler, 1998).
The smoothness constraint is similar to a quadratic regularization (Tikhonov, 1963). If the x vector is reshaped into a rectangular matrix, then a two-dimensional gradient matrix D can be formed. A standard measure of the gradient would then be the sum of the values of Dx for every pixel. Defining D in this way results in a solutionx that is globally smooth. Since the ionosphere is generally continuous (and therefore mostly smooth), this smoothness constraint greatly increases the quality of the reconstruction in preserving the background unperturbed ionosphere.
Edges in an image are the areas that would have the highest gradient values. Therefore, in order to preserve edges, high values of the gradient should have a smaller penalty, while still maintaining a penalty on low values of the gradient to maintain smoothness. The key, then, is to create a weighting function to vary the cost depending on the value of the gradient. This formulation leads to the following smoothness regularization functional: where n refers to individual elements in the vector Dx.
The following weighting function φ(t), where t is the value of the gradient, is used to enhance the quality of the reconstruction.
where T is a parameter than can be adjusted to alter the shape of the weighting function. This nonconvex weighting function will resemble a quadratic weighting function when t T and will be nearly constant when t T . Applying this φ function to our cost function gives the following expression: where λ is a regularization parameter and is adjusted to provide balance between the least-squares fit to the measurements (first term in Eq. 8) and the regularization functional (second term in Eq. 8). All that remains is then to find the optimal value of x for this cost function, given the observation matrix A and the data y.

Optimization method
In most regularization problems, the cost function is globally convex. Therefore, the goal of the minimization is to reach the "bottom" of the function. Convexity ensures that a suitably-designed method will reach the global minimum, since there are no local minima outside of the global minimum. A common approach to solving convex minimization problems is through the use of the conjugate gradient method (Moon, 2000). However, in order to reduce the penalty for large gradient values and preserve edges, we utilized a nonconvex cost function. This weighting will enhance results but will no longer allow the simple use of the conjugate gradient algorithm for minimization. An alternate method for optimization, a deterministic relaxation technique, can be used to solve this nonconvex minimization problem (Delaney and Bresler, 1998). The goal of the deterministic relaxation technique is to create a succession of convex minimization problems that will converge to the global minimum. This is done through the introduction of an additional parameter, e(x), that captures the nonconvex portion of the cost function. The cost function is then a function of both x and e(x). An initial guess, x, is used to calculate the e vector. This value e(x) is fixed and the new cost function J (x, e(x)) is minimized with respect to x. The new cost function is convex with respect to x, so the minimization is done using the conjugate gradient method. This minimum value is then used as a new This iterative progression is guaranteed to converge to the solution. The proper choice of regularization parameters is necessary to achieve the proper degree of smoothness, distinguish edges from noise, and to ensure positivity. While the application of this optimization technique to GUVI data yields good reconstructions of bubble structure, the images contain artifacts that distort the background ionosphere. The most noticeable effect was the tendency of the solutions to have unnaturally high electron density values near the satellite positions. Another problem was that the limited coverage of the satellite view produced poor reconstructions in the top corners of the image, areas which are never in the line of sight of the GUVI instrument.
An effective way to deal with this problem is to introduce additional constraints that are designed to force the solution to reside in a set of "realistic" ionospheres. The technique of projection onto convex sets (POCS) (Kamalabadi and Sharif, 2005) was incorporated into the optimization method in a way that improved the solution while still maintaining enough flexibility to conform to the real ionosphere.
Projection on convex sets was used to impose two additional constraints: a reference constraint and an amplitude constraint. The convex set C R for the reference constraint is: This constraint forces the solution x to be within a distance of the reference ionosphere, x R . The projector P R that enforces this constraint is The convex set C A for the amplitude constraint is: This constraint forces the solution to contain electron density values within a reasonable range. If a is taken to be 0, then this forces the solution to be positive. The value for b could  be obtained by a variety of means, for example obtained by the maximum electron density value in the Parameterized Ionospheric Model (PIM) (Daniell et al., 1995) simulation. The projector P A that enforces this constraint is The reconstruction is not very sensitive to the positivity parameter, provided it is large enough to enforce the constraint. The effect of the other two regularization parameters is correlated, which makes the determination of optimal values difficult. Furthermore, since the choice of regularization parameter is sensitive to the data values, a heuristic approach to regularization parameter selection is reasonable. Generally the results did not differ more than a factor of 2 for each of the regularization parameters. Typical values for the regularization parameters are λ=2.5×10 −22 and T =5.2×10 22 . This inversion technique was validated extensively in Comberiate et al. (2006) through both simulations and application to GUVI data. Analytical measures of estimation uncertainty were also provided for this algorithm in that work. Here, in the next section we present the results of the application of this technique to GUVI data during the SpreadFEx campaign and discuss comparisons with the ground-based red-line imager.

Results
GUVI and the ground-based airglow imager provide complementary observations in the sense described below. Both imaging techniques view the same three-dimensional section of the ionosphere but only produce cross-sectional images in two dimensions. The ground-based imager assumes that the emission originates only from 250 km and does not capture any altitude structure. The GUVI image assumes that the electron density in the region is constant along magnetic field 2446 F. Kamalabadi et al.: Electron densities from GUVI tomography during SpreadFEx   lines and does not capture any latitudinal structure. Used together, the two images can capture the three-dimensional structure of plasma bubbles with variations in longitude, altitude, and latitude.
GUVI reconstructions provide altitude vs. longitude electron density profiles. The ground-based airglow imager provides latitude vs. longitude images. The ground-based imager is too close to the magnetic equator to map along magnetic field lines. Therefore, validation can only be performed through a one-dimensional longitudinal comparison. This can be done by integrating the GUVI image along the altitudinal direction and integrating the ground-based image along the latitudinal direction. The shapes of the two longitudinal profiles should be comparable. The level of agreement between the two profiles can serve as a means of validation. Since the airglow images are a measure of OI intensity and the extraction of electron density from this measurement is not straight-forward and at best would require complicated modeling, here we focus on morphological rather than quantitative comparisons.
GUVI precesses through 12 h of local time approximately every 60 days. GUVI reconstructions use 3 min of data, corresponding to a 10 • segment of the orbit. The latitude value given with each image is the center point of that 10 • segment. The reconstruction is valid for ±5 • in latitude. An offset, tilted dipole magnetic field model is used and the cross-section can be projected along magnetic field lines. For the entire SpreadFEx campaign interval, GUVI viewed the nightside ionosphere on the descending (moving southward) portion of the orbit.
Each figure provides altitude vs. longitude profile of electron density values with labeled axes and a colorbar indicating the values in the image. The figure title provides the date, latitude, and universal time for the reconstruction. Examples of GUVI reconstructions during four of the campaign nights are shown in Fig. 3. As can be seen, these electron density reconstructions capture a rich spectrum of plasma structures during the formation and evolution of plasma bubble. Figure 4 shows the corresponding observations of several EPBs with GUVI and the all-sky imager on 2005 Day 274. This comparison illustrates both the capacity and the limitations of GUVI's longitudinal resolution of EPBs. GUVI observed multiple plasma bubbles at 21:18 LT (01:18 UT) that correspond to those that have drifted East into the imager field of view at 02:18 UT. GUVI sees three large depleted regions. The vertical bubble with a width of approximately 180 km (at 450 km altitude) at −65 • E corresponds to the multiple faint depletions seen at −50 • E in the all-sky image. The bubble with an approximate 30 • eastward tilt and 180 km width (at 300 km altitude) at −62 • E corresponds to the strongly depleted region seen at −48 • E in the all-sky image. The vertical bubble with a width of approximately 150 km (at 250 km altitude) at −65 • E corresponds to the thick depleted region seen at −46 • E in the all-sky image. The resolution of the GUVI reconstruction is limited by the 40 km longitudinal pixel width of the reconstruction grid, so the thinner features seen by the imager can not be resolved in the GUVI reconstruction. The depleted regions at lower altitudes in the GUVI reconstructed image are more distinct in the 630.0 nm image. Figure 5 shows a plot of the correlation between the longitudinal structures in the two images for Day 274, restricting the GUVI image to the emission layer between 230 and 270 km in altitude and assuming the structures in the GUVI image drift eastward and overlay the imagers longitudinal field of view. The correlation coefficient between the two datasets is 0.51. As described earlier, both images detect multiple bubble structures but the correlation is somewhat weakened by the difference in longitudes between the two images and the hour interval over which the bubbles evolve and drift. Figure 6 shows another set of corresponding observations from 2005 Day 275. Both images show a large area of low electron density in the eastern region and a single bubble in the west. The depletion in the GUVI image has an approximate 40 • westward tilt and a 150 km width at 450 km altitude which widens to 300 km at 560 km altitude. The bubble is centered at −47 • E at 350 km altitude and −49 • E at 550 km altitude. In the 630.0 nm image, the bubble has a width that varies from approximately 40 km to 150 km and is located between −48 • E and −50 • E. The large depleted region in the East of the GUVI reconstructed image is not an equatorial plasma bubble but appears in this case because there is substantial latitudinal variation in the background intensity over the eastern portion of the field of view of the GUVI instrument. This background gradient can be clearly seen in the 630.0 nm image. Figure 7 shows a plot of the correlation between the longitudinal structures in the two images for Day 275, again restricting the GUVI image to an emission layer between 230 and 270 km. In this case both datasets cover the same longitude span, resulting in a much higher correlation coefficient of 0.95. Both images capture a bubble structure and a substantial dropoff in electron density from West to East.

Conclusions
TIMED/GUVI data can be used to reconstruct multidimensional profiles of equatorial plasma bubbles. These reconstructions indicate the width, tilt, and depth of depletion of the plasma bubble. Coincident observations with the ground-based 6300Å airglow imager provide complementary information of altitude and high spatial resolution. The GUVI reconstructions provide a unique view for imaging and characterization of equatorial plasma bubbles.