Application of generalized singular value decomposition to ionospheric tomography

. The electron density distribution of the low- and mid-latitude ionosphere has been investigated by the computerized tomography technique using a Generalized Singular Value Decomposition (GSVD) based algorithm. Model ionospheric total electron content (TEC) data obtained from the International Reference Ionosphere 2001 and slant relative TEC data measured at a chain of three stations receiving transit satellite transmissions in Alaska, USA are used in this analysis. The issue of optimum efﬁciency of the GSVD algorithm in the reconstruction of ionospheric structures is being addressed through simulation of the equatorial ionization anomaly (EIA), in addition to its application to investigate complicated ionospheric density irregularities. Results show that the Generalized Cross Validation approach to ﬁnd the regularization parameter and the corresponding solution gives a very good reconstructed image of the low-latitude ionosphere and the EIA within it. Provided that some minimum norm is fulﬁlled, the GSVD solution is found to be least affected by considerations, such as pixel size and number of ray paths. The method has also been used to investigate the behaviour of the mid-latitude ionosphere under magnetically quiet and disturbed conditions.


Introduction
Computerized tomography (CT) technique has had a revolutionary impact in the field of medical diagnostic imaging and is now routinely used to produce high quality images of a section of the human body. The mathematical foundation of CT was laid as early as 1917 by Radon (Deans, 1983). The technique has also been applied to the study of ocean structures (Munk and Wunch, 1979) and geological formations (Dines and Lyttle, 1979). CT technique can be used for investigation of the vertical and horizontal structure of the ionosphere when the line integral of the parameter (e.g. electron density) is the measured data through the region of interest.
Correspondence to: P. K. Bhuyan (bhuyan@dibru.ernet.in) A radio signal passing through the ionosphere suffers time delay, refraction and polarization rotations during its passage through the ionized plasma. Since the beginning of the satellite era, radio beacons on board artificial satellites have been routinely used to measure the total ionospheric electron content (TEC) using techniques such as Faraday rotation, group delay, Doppler shift (de Mendonca and Garriott, 1962;Davies, 1980Davies, , 1989. TEC is a measure of the total number of electrons contained in a tube of unit cross sectional area along the line of sight from the transmitter to the receiver. Individual and coordinated TEC measurements provide information about the temporal and horizontal distribution of electron content on various time scales. However, since TEC is the line integral of the electron density along the ray path, it cannot be used for investigating the vertical structure of the ionosphere using conventional techniques. Austen et al. (1986Austen et al. ( , 1988 first proposed the application of computerized tomography technique for reconstruction of the electron density distribution below and above the F2 peak by utilizing TEC observations. However, for utilizing TEC obtained from absolute phase measurements, an additional procedure of initial phase estimation is necessary, which leads to significant difficulties and reduces the reconstruction accuracy. Pryse and Kersley (1992) reported preliminary experimental results based on the method proposed by Austen et al. (1986Austen et al. ( , 1988. Andreeva et al. (1990) and Kunitsyn and Tereshchenko (1992) proposed a phase-difference approach making use of relative TEC derivatives and obtained experimental ray tomographic (RT) images of the ionosphere for the first time. The genre of computerized ionospheric tomography (CIT) that followed Austen et al. (1986Austen et al. ( , 1988 basically adheres to the principle of measurement of total electron content at a meridional chain of ground receivers using a transit satellite and then inverting the same to reconstruct a two-dimensional electron density versus height profile. Raymund et al. (1990) applied the computerized ionospheric tomography (CIT) to realistically simulate the ionospheric electron density variations over 16 • in latitude within height range of 50 to 1000 km.
A number of theoretical as well as experimental CIT schemes have been reported so far in literature (Raymund et al., 1990;Fremouw et al., 1992;Na and Lee, 1990, 1994Afraimovitch et al., 1992;Andreeva et al., 1990Andreeva et al., , 1992Andreeva et al., , 2001Tereshchenko, 1991, 1992;Klobuchar et al., 1992;Pryse and Kersley, 1992;Kunitake et al., 1995;Kersley et al., 1997;Zhou et al., 1999;Tsai et al., 2002;Kamalabadi et al., 2002;Bhuyan et al., 2002). The Algebraic Reconstruction Technique, ART (Austen et al., 1988), the Multiplicative Algebraic Reconstruction Technique, MART (Raymund et al., 1990;Kersley et al., 1993) and Maximum Entropy Method, MEM (Fougere, 1995) are some of the widely used CIT algorithms. Kunitsyn et al. (1994) proposed the DART algorithm (Decomposed Algebraic Reconstruction Technique) and obtained good results in ionospheric tomography problems. It contains the first two modified members of MART Taylor Decomposition. Yeh and Raymund (1991) and Raymund et al. (1994) discussed some of the limitations of ionospheric imaging by tomography. Raymund (1995) has made a comparative assessment of various CIT algorithms and observed that no single algorithm can be considered as the best. Some algorithms that construct well in one case can do very poorly in other cases. Some algorithms do well in parts (e.g. bottom-side gradients) but not as well at others. It was also observed that the fundamental assumptions make significant differences and proper choice of the algorithm will be the prime requisite for meaningful data reconstruction. The ART and MART algorithms are conceptually simple and computationally efficient, which make them choices as general purpose CIT algorithms. The ART algorithm used to reconstruct the images seeks to minimize the root mean square difference between the observed TEC data and TEC data computed from the reconstruction. The ART algorithm is iterative and requires some starting image as an initial guess. The algorithm computes the root mean squared difference at each step of iteration and then makes additive changes that seek to minimize the difference. In many ways MART is similar to ART. It is iterative and starts from an initial guess. It compares TEC computed from the initial guess with the measured TEC. At each iteration step, changes to the image are based on the difference between the initial guess and measured data. Of course, the difference of MART with ART is that the changes are multiplicative rather than additive. More detailed discussion of algorithms used in tomography can be found in Raymund (1994), Pryse et al. (1998) and Kunitsyn and Tereshchenko (2003).
In this paper, we investigate the electron density distribution of the low-and mid-latitude ionosphere through tomographic reconstruction using the Generalized Singular Value Decomposition algorithm reported by Bhuyan et al. (2002). The TEC data used in this analysis for reconstruction of the low-latitude ionosphere are obtained from the International Reference Ionosphere 2001. Further, slant TEC data obtained from a chain of three receiving stations of LEO satellites in the Alaska region, USA, and made available through the Internet by International Ionospheric Tomography Community, are used to study the dynamics of the mid-latitude ionosphere. We seek to address the issue of optimum efficiency of the GSVD technique in the reconstruction of largescale ionospheric structures, such as the equatorial ionization anomaly in addition to its application to investigate smallscale ionospheric density irregularities. The algorithm is further used to investigate the growth and development of the mid-latitude trough in electron density.

Theory
In a CIT experiment, the phase difference suffered by two coherent radio signals transmitted from radio beacons on board transit satellites is measured using a meridional chain of ground receivers. For example, a number of CIT experiments (Kersely et al., 1993(Kersely et al., , 1997Foster et al., 1994, etc.) have been performed using 150 MHz and 400 MHz transmissions from NNSS, Cicada satellites. The dispersive phase φ (phase difference between the two signals reduced to a reference frequency) is proportional to the Total Electron Content (TEC) along the path S connecting the satellite and receiver and is given by (Andreeva et al., 1992) where λ is the probing wavelength, r e is the classical electron radius, k is proportionality constant and S Ndσ is the TEC along the path S. k is of the order of unity and is connected to the reference frequency at which the dispersive phase is measured. This absolute phase approach to the ray tomography of the ionosphere suffers from the problem of initial phase determination known as phase ambiguity or 2nπ ambiguity. This ambiguity in phase measurement causes substantial error in the reconstruction process (Kunitsyn and Tereshchenko, 2003). However, the problem of initial phase determination can be avoided, if the phase-difference tomography as suggested by Andreeva et al. (1992) is adopted.
Subject to appropriate discretization, the ray tomography problem transforms into a problem of solving a system of linear equations. The system of linear equations can be written as where A M×N is the operator matrix for transition of the sought function (i.e. the function to be reconstructed, x N×1 in our notation) to the linear integrals (y M ×1). The matrix A takes different forms depending on whether the approach is based on phase or phase difference. Its form also depends on the approximation used for the function sought. A detailed discussion can be found in literature . The different orders of approximations for the function sought, that was studied by Kunitsyn et al. (1994), are the piecewise constant, piecewise planar, the product of linear approximations, the product of cubic approximations, etc. However, it should be noted that as the order of approximation increases, the matrix properties for solving the inverse problem worsen.
In this paper, the matrix A has been calculated using piecewise planar approximation of the function, sought following Andreeva et al. (1992). The piecewise constant approximation has been avoided, as it leads to the largest error in reconstruction . The discretization process and the intersection of a single ray path with only a very small fraction of the total number of the pixels makes the geometry matrix largely sparse and hence highly ill conditioned. Kunitake et al. (1995) used a method known as Modified Truncated Singular Value Decomposition (MTSVD) for solving such ill-conditioned problems. This method is based on the truncation of singular values causing larger perturbations in the solution. The choice of the truncation parameter is arbitrary. Zhou et al. (1999) derived an optimal truncation criterion for using SVD for reconstruction purposes. The approach is based on the inversion of a weighted geometry matrix constructed by incorporating a priori information in the form of a priori and a posteriori covariance for both the data and the model parameters (function sought). In the inversion of the weighted geometry matrix, the series of the weighted singular values has been truncated at unity. The truncation criterion has been chosen due to the fact that singular values greater than unity decrease the trace of the a posteriori covariance relative to the a priori covariance, while smaller ones do not. The smaller the trace of the a posteriori covariance, the lower the uncertainty is in the solution. Bhuyan et al. (2002) reported a new CIT algorithm based on the Generalized Singular Value Decomposition (GSVD). For the sake of completeness of the paper the algorithm is discussed below in brief.

The algorithm
Ionospheric tomography is a part of the family of ill-posed problems, which can be written in the form of a system of linear equations given by Eq. (2). A common way to solve Eq. (2) is to generate the estimated value of x as the leastsquare solution to the set of equations where ||·|| denotes l 2 norm of a vector. When A is ill conditioned, the least-square reconstructed object x ls obtained through Eq.
(3) will be corrupted by amplified noise. There are many methods for regularizing such problems in order to generate reasonable estimate. Tikhonov regularization is one such method (Tikhonov, 1963). The GSVD is a generalization of the singular value decomposition (SVD). The SVD of a matrix A is a decomposition of the form where the sets {u i } and {v i } are orthonormal and λ i satisfy λ 1 ≥λ 2 ≥λ 3 . . . ≥λ m ≥0.
The GSVD can be used to solve the damped least-squares problem, as proposed by Tikhonov (1963). This approach amounts to finding out the x that solves where α is the regularization parameter and α>0. L is a positive definite matrix. L takes different forms in accordance with the order of regularization. For zero order Tikhonov regularization, L=I, the identity matrix and for first order Tikhonov regularization where L is a P ×N matrix. Under GSVD technique, the matrices A M×N and L P ×N with M≥N ≥P , and rank (L)=P , can be written as and where U is M×N and orthogonal, V is P ×P and orthogonal, and B is N×N and non-singular. is P ×P diagonal matrix with elements λ i such that 0≤λ 1 ≤λ 2 ≤λ 3 ≤λ 4 . . . ≤λ P ≤1 C is a P ×P diagonal matrix with elements µ i such that The values of λ i and µ i are normalized so that The generalized singular values are such that 0≤γ 1 ≤γ 2 ≤γ 3 ≤ . . .≤γ p . It can be easily seen that and when P <N , the matrix L has a nontrivial null space. The vectors B P +1 , B P +2 , . . . , B N form a basis for the null space of L.
Equation (3) can now be written as The normal equations are or, The Tikhonov regularized solution can be written from the above equation as The solution x is a function of the regularization parameter α, as well as L. The solution varies with the regularization parameter quite strongly as compared to L and unfortunately, in the damped least-square approach, the value of α is not specified. A key element of any Tikhonov regularized formulation is the proper choice of the regularization parameter α. If α is chosen to be too small, the reconstruction will be dominated by large, high frequency components. If α is chosen to be too large, important information in the data will be suppressed. The indeterminacy of α can be eliminated by the method of Generalized Cross Validation (GCV). Let Here V(α) measures the misfit. As α increases, V(α) also increases. T(α) is a slowly increasing function of α. A regularization parameter α is selected, for which G(α) has a minimum value. Bhuyan et al. (2002) proposed the determination of the regularization parameter α using the method of Generalized Cross Validation (GCV). GCV has been used for the purpose of blur estimation and resolution enhancement in image processing (Nguyen et al., 2001). The solution also varies with the form of L, i.e. the order of regularization. But, as the regularization approach tries to find a solution which is vertically as well as horizontally smooth, one has to choose an optimum regularization. In case of second order regularization, the large gradients are heavily penalized (Kamalabadi et al., 2002), thus tending to hide important aspects of the image of a highly disturbed ionosphere. A first order regularization should be a better option in such cases. In this paper, first order regularization has been used in all the cases.
In contrast to the MTSVD technique, the GSVD approach is based on the method of regularization of the problem, which eliminates the need for truncation of the singular values. It is noteworthy that the matrix L incorporates information regarding the smoothness of the solution. Lx for first order regularization is a finite difference approximation to the first derivative of x and for the second order; it is a finite difference approximation to the second derivative of x.

CIT reconstruction of the low-latitude ionosphere
The well-known equatorial ionization or "Appleton" anomaly (EIA) over Indian low latitudes was investigated by reconstructing the same using the GSVD algorithm developed for this purpose and described in Sect. 2. The ionosphere within the altitude range of 100 km to 1000 km and latitude range of 18 • N to 28 • N was divided into pixels of 55 km×100 km dimensions (details in Table 1). Bhuyan et al. (2002) used lower resolution 110 km× 100 km pixel size for reconstructing the low-latitude ionosphere within a latitude span of 6.5 • . The TEC along 220 different ray paths within the region of interest was obtained from the International Reference Ionosphere (IRI-2001). Figure 1a illustrates the electron density distribution for 20 January 2004 at 12:00 LT along 75 • E meridian. The reconstructed electron density distribution is plotted in Fig. 1b. The regularization parameter α has a value of 5.85. It is seen that the equatorial ionization anomaly (EIA) is reconstructed particularly well. The tilt of the ionosphere is also well preserved. Andreeva et al. (2000), Yeh et al. (2001) have shown by radio tomographic investigations that EIA crest can be quite well pronounced. Also, the anomaly core was found to be oriented in the direction of the geomagnetic field. However, the IRI is a climatological model that gives monthly averages of ionospheric parameters under quiet solar activity conditions and since the TEC data used for the present reconstruction of the EIA, is obtained from the IRI, investigation of structural peculiarities of EIA as obtained by tomographic experiments such as the ones mentioned above, could not be undertaken.

Effect of geometry on reconstruction
The computations carried out with the GSVD-based algorithm require large computer disc space, which, in turn, adversely affects the resolution one tries to achieve. In order to investigate the effect of the geometry on the resolution of the reconstructed image, another geometry with increased number of sampling ray paths was considered. The number of sampling ray paths in this case was increased to 250 from the initial 220 within the same altitude latitude range described in Sect. 3.1. The reconstruction is shown in Fig. 1c. The regularization parameter α has a value of 14.0737 in this case. As the number of ray paths increases and the number of pixels increases the sparsity of the geometry matrix, the solution should have become more non-unique and unstable. However, the Generalized Cross Validation approach to find the regularization parameter and the corresponding solution gives a very good reconstructed image, with the prominent features reconstructed particularly well. The reconstruction is poor only near the edges of the image. It may also be noted that the increase in resolution keeping the same latitudinal span increases the sparsity of the geometry matrix, which can, in turn, adversely affect the reconstruction.

Reconstruction of complicated irregularities
The ionosphere has a complicated structure where irregularities of different dimensions cohabit the regular structures, such as EIA, mid latitude trough, etc. For the reconstruction of such complicated ionospheric structures, a model ionosphere was first constructed by superimposing three irregular structures on a background ionosphere obtained from the Chapman function, as shown in Fig. 2a. Figure 2b shows the reconstructed ionosphere using the geometry described in Sect. 3.1. All the irregularities preserve their positions in latitude and altitude in the reconstructed image. It should be noted that the regularization approach introduces a smoothening effect in the gradients of electron density (Kamalabadi et al., 2002). It has been noticed during the simulation of different test cases that in case of a highly disturbed ionosphere, first order regularization is optimum. Second order regularization has a much stronger smoothening effect. The regularization parameter α has a value of 0.39131 in this case.

Reconstruction of mid-latitude ionosphere
Earlier simulation studies (Andreeva et al., 1992) have shown that a minimum of three receiving stations could form a good platform for performing reconstruction of the ionosphere from total electron content data. The GSVD algorithm is   Figure 4a represents the reconstruction for 26 February. The midlatitude trough is reconstructed reasonably well. Figure 4b represents the reconstruction for 14 July. Departures from the average quiet day behavior of the midlatitude ionosphere are quite evident from the reconstruction and it matches with what is expected from Figure 3b. The variation of relative TEC with arrival angle of the signal at all the receivers suggests the presence of such a structure in the ionosphere. It should be noted that the latitudinal span of the two reconstructions is different, owing to the availability of slant TEC data.

Discussion
The formulation of phase-difference tomography makes it free from the phase ambiguity problem. The reconstruction result is affected only by the presence of noise and random recording errors (Andreeva et al., 1992). These errors are averaged out during processing of data and as such, the reconstruction becomes more effective. A reconstruction error can be characterized by relative errors ρ l 2 and ρ c as being discrete analogue for the deviation The values of the parameters ρ l 2 and ρ c for reconstruction of EIA for the two cases are 0.4251, 0.035 and 0.6148, 0.0093, respectively. In case of reconstruction of complicated irregularities, the respective values are 0.0789 and 0.019. Thus, comparison of ρ l 2 and ρ c for the first two cases shows that while the overall reconstruction error is lower in the first test case compared to the second test case, the reconstruction error of the maximum value of the sought function (x) is higher. The low value of both the parameters in the last case points to better reconstruction. The indeterminacy of the regularization parameter α has been removed by using GCV, which makes the algorithm self-consistent. It must be noted that the GSVD approach requires the number of ray paths to equal or exceed the number of pixels. The algorithm is capable of reconstructing regular as well as complicated structures.

Conclusion
A radio tomographic algorithm based on the Generalized Singular Value Decomposition technique is used to investigate the regular and irregular structure of the F region ionosphere at equatorial and mid latitudes. Model and measured TEC data are used to reconstruct the equatorial ionization anomaly and the mid-latitude ionosphere, respectively. The method has also been used to reconstruct complicated density irregularities embedded in a simulated regular ionosphere. Results indicate that this GSVD based CIT algorithm is capable of recreating the electron density distribution of the ionosphere satisfactorily. Another advantage of the method has been its ability to reconstruct without the help of any initial guess regarding the state of the ionosphere in the region of interest.