Back to INAG Homepage Back to UAG-104 Contents page

Simplified Numerical Procedures Applied to the True-height Analysis of Ionograms

L.-C. Tsai
Space Dynamics Laboratory and Department of Electrical Engineering, Utah State University, Logan, UT 84322-4120, USA

F. T. Berkey
Center for Atmospheric and Space Sciences, Utah State University, Logan, UT 84322-4405, USA

G. S. Stiles
Department of Electrical Engineering, Utah State University, Logan, UT 84322-4120, USA

Abstract

An efficient numerical process can enable true height calculations to be readily carried out on small computer systems. We have developed a numerical integration method termed "m't-fitting", where m' is the group refractive index and t is defined as Ö(1-fN 2/f 2) for the O-wave and Ö[1-fN 2/(f 2-f fH)] for the X-wave. This method is best suited to the true-height analysis of digitally recorded ionograms. It can be used to analyse either the O- or X-waves in an ionogram and includes the effect of the Earth's magnetic field. For second order polynomial profile analysis, the m't-fitting technique utilises 38-62% fewer numerical operations than the Gaussian quadrature integration method. Applying 5 and 12 term m't-fitting to 100 data points, the method completes a true-height analysis in 0.77 and 1.43 seconds, respectively, using a 80486/33 computer system. Furthermore, applied to the overlapping Chapman model profile adopted by U.R.S.I. Working Group G.6.2 [McNamara and Titheridge, 1977], the test results show rms errors of 8.3 and 7.0 meters for the O and X wave components (12-term m't-fitting), over a set of 22 selected frequencies.

Introduction

Ionograms have provided the basis for much of our understanding of the physics of the ionosphere. Even more information can be obtained from these records by inverting the virtual height-frequency trace to obtain the ionospheric electron density profile, a process called true-height analysis. The proliferation of digitally recorded ionograms and the development of powerful microcomputers now enable true-height analysis methods to be readily implemented.

Titheridge [1985] developed one of the most widely used true-height analysis methods using polynomial equations to represent h( fN ), where h is the true height and fN is commonly expressed as plasma density. The POLynomial ANalysis program (POLAN) is a computer method that implements a generalised, polynomial analysis, the execution time of which is primarily dependent on the group refractive index calculation. In this paper, we describe an efficient method that can be used to derive the integral of the group refractive index m', and enable true-height analysis of digital ionograms to execute in less than one second for 100 data points. A monotonic ionosphere is used for the purpose of a simplified comparison of the two true-height analysis methods described previously.

"m't-fitting" technique

The electromagnetic field in the ionosphere is governed by Maxwell's equations. It is assumed that there is no permanent space charge and that all field variables vary sinusoidally with time at the same angular frequency w. Ionospheric sounding utilizes radio waves transmitted vertically from the ground, which travel upward into the ionosphere, where they are reflected, and return to a ground-based receiver. The measured virtual height h' and the true height h are related by h' = ò0hr m'dh, where hr is the height of the reflection, and m' is the group refractive index, a function of the radio frequency f, the plasma frequency fN, the magnetic gyrofrequency fH, and the magnetic dip angle I.

The group refractive index m' is a complicated function, as expressed in Shinn and Whale [1952], and approaches infinity at the reflection point. In order to avoid singular points in the integrations of m', Shinn [1955] and Becker [1960] introduced the reduced frequency fR. This quantity is defined to be f for the O-wave and (f 2-f´fH)0.5 for the X-wave. By changing variables, we make the following substitution:

At the reflection point, the product m't has a finite value, except for the ordinary component at 90° magnetic dip angle. Numerical methods permit approximation of the product m't by an interpolating polynomial in t which fits a selected set of points (ti, m'i ti), where i = 0, 1, 2, .., n. Such a polynomial is unique, because there is only one polynomial of degree n passing through n+1 points. The Newton-Gregory interpolating polynomial [§5 of Cheney and Kincaid, 1980] can be used to find the polynomial that passes through a group of equispaced points for t Î [0,1]. Hence, the product m't can be derived as a function of t as follows

where s = t/d, and d is the uniform difference between the t-values. With this m't-fitting technique, the series terms of the group refractive index m' can be calculated by the Shinn-Whale formula, and each coefficient of the Newton-Gregory forward interpolating polynomial can be derived. We assume that the variation of plasma frequency fN with true height h is given by

where the point (fa, ha) is the origin for the real-height region under consideration, and the value of NT defines the order of the analysis. This model yields

where h" is the reduced virtual height, equal to h' minus the group retardation due to those parts of the profile below the point (fa, ha). Therefore, by the m't-fitting technique and replacing fN2 with f, we can express the integral in equation (4) in terms of this known polynomial as

 

and

Complexity analysis

Titheridge [1985] developed the FORTRAN program POLAN that implements a generalised polynomial method that can be applied to the true-height analysis of ionograms. Ten standard modes of analysis are provided in POLAN, dependent on: 1) NT, the order of the polynomial invoked for real-height analysis; 2) NV, the number of virtual heights; 3) NR, the number of known true heights and; 4) NH, the number of new true heights calculated in each step. The coefficients of true-height polynomials yield an exact solution if NT = NV + _NR_, or a least squares solution if NT < NV + _NR_. Calculations normally proceed in a stepwise fashion to obtain successive, overlapping sections of the real-height profile. POLAN invokes a five- or twelve-term Gaussian quadrature integration, depending on the desired accuracy. In the case of the NG-term used in the Gaussian quadrature method, since all NT values of integration for equation (4) are obtained from the same NG values of m', there are a total of {(NT + 11) NG - 1} additions and {(NT + 19) NG + 2 NT} multiplications required to calculate these integration terms.

/

For the m't-fitting technique, since the gyrofrequency fH, the magnetic dip angle I, and all radio frequencies are known, the coefficients of the Newton-Gregory interpolating polynomial for the product m't as a function of t can be precalculated. If we assume NG terms (n = NG - 1) for the Newton-Gregory interpolating polynomial, there are a total of {3 NTNG + 0.5 NT2 + 1} additions and {(2 NT + 2) NG + 2.5 NT2 + 8 NT - 1} multiplications required to calculate NT values of integration for m'.

Figure 1 shows the number of mathematical operations (addition and multiplication) required for the two integration techniques, as a function of NG between 3 and 20. With increasing the polynomial order, this diagram shows that the number of operations increases proportionally more for the m't-fitting technique (100% for NT = 4) than for the Gaussian quadrature method (13% for NT = 4) at a NG value of 12. The smaller increase in the number of numerical operations for the Gaussian quadrature method results from using the same NG values of the group refractive index m' for every integration in equation (17). None-the-less, except for small values of NG in the fourth order case, more mathematical operations are used by the Gaussian quadrature method since more complex calculations are required to derive the group refractive index m'. For the parabolic lamination method (NT = 2), the m't-fitting technique utilizes 38-62% fewer numerical operations than the Gaussian quadrature method. For second order polynomial profile analysis, Figure 5 shows the execution times of the 5-term and 12-term m't-fitting techniques as a function of echo number N, as implemented on an 80486/33 computer. As expected, the execution times of these methods are approximately proportional to the square of the number of echoes N. They range from 0.77 (1.43) seconds for 100 data points to 3.14 (5.63) seconds for 200 data points when the 5(12)-term m't-fitting techniques are applied.

Accuracy analysis

In general, the Newton-Gregory interpolating polynomial does not produce results which are identical with the real product m't derived from the Shinn-Whale formula, even though common values can be found. The non-matching values obtained from the interpolating polynomial are subject to error, which decreases (as expected) as the order of the interpolating polynomial is increased. A convergence test for various orders of the interpolating polynomial reveals the accuracy which can be expected and provides a basis for selection of the order of polynomial. The errors are computed relative to results derived using an 18-term Gaussian quadrature integration with 103 intervals over t Î [0,1], in which all abscissas and weights were obtained from Love [1966] and have an error less than 1.5´10-23. From the results, the worst case error of the X-wave is two orders of magnitude smaller than that of the O-wave; i.e. the interpolating polynomial for the mx't curve has smaller error than that derived for the mo't curve. The magnitude of relative error is less than 1.2´10-3 for a 5-term m't-fitting integration and less than 2.2´10-6 for a 12-term m't-fitting integration in the case of the X-wave, I=65°, and Y³0.1. In order to decrease the error associated with the O-wave, we simulate the product mo't by using three non-overlapping polynomials instead of one over the interval t Î [0,1]. The points of intersection can be obtained by trial and error or found approximately from the curve for I = 65°, which varies by a factor of 2 between t=0.1 and 0.3. In this analysis, we chose the points of intersection at 0.1 and 0.3. The relative errors for this approach are less than 8.5´10-4 for 5 terms in each polynomial section and less than 2.5´10-6 for 12 terms in the case of I=65° and Y³0.1. The associated execution time, which is the same as that of adding two more data points to the single-m't-fitting method, is increased slightly but is still a factor of 2 faster relative to the Gaussian quadrature integration method, when the number of data points exceeds 10.

Inadequate starting height data is a main source of error in most methods of true-height analysis. The simplest analysis makes the assumption that the virtual height of the lowest observed frequency is the first point of the electron density profile. Thus, the presence of underlying ionization is ignored. This is generally a good approximation, especially when the first observed frequency occurs below the critical frequency of the E layer and when the D layer is not present. A different method can be applied when there is no critical frequency below the lowest observed frequency. Here, we assume that the true heights of the starting region in an ionogram vary linearly with f and that the higher order derivative is equal to zero. Therefore, the starting equations can be expressed as

where z' is constant over the starting range. Since there are two unknowns, ho and z', in equation (7), the value of i must be at least two. If i is greater than two, we have more than two equations and must solve for the two unknowns by a least squares method.

A model profile adopted by U.R.S.I. Working Group G.6.2 [McNamara and Titheridge, 1977] consists of two overlapping Chapman layers, in which the Chapman tail has been cut off at 0.8 MHz. Virtual heights were calculated using a 40-point Gaussian integration applied over five successive height ranges, yielding errors of less than 1 m. We have used the 5- and 12-term m't-fitting techniques to analyze 150 data points in which the frequency interval is ³ 0.02 MHz as recommended by McNamara and Titheridge [1977]. Table 1 shows true-height errors (in meters) obtained for the O- and X-wave components at 22 selected frequencies. Application of the 12-term m't-fitting technique to these profiles yields a difference error for the O-wave component that is a factor of ~5 smaller than that obtained with the 5-term m't-fitting technique. However, when either the 5- or 12-term m't-fitting technique is applied to the X-mode data, the difference error is essentially the same for each frequency-height data pair. This result is consistent with a smaller error in the interpolating fitting polynomial for the mx't curve than that derived for the mo't curve. The rms deviations of 23.5 and 7.0 meters for the X-wave component are the result of the approximation techniques that are employed.

Table 1. The errors in true-height analysis of the overlapping Chapman profiles as presented by McNamara and Titheridge [1977].

/

Discussion and summary

In POLAN, a polynomial representation of the true-height profile is used to the inversion of ionograms. Although Gaussian quadrature integration is a very accurate numerical integration method, the computational overhead of this method is large due to the complicated mathematical calculations required to derive the various group refractive indices. However, in most instances of ionospheric sounding, the transmitted frequencies are well known. Using the characteristics of polynomial analysis, we can also express the product of the group refractive index and the variable t (defined in equation (1)) as a polynomial of t at a given radio frequency. Since all of the integrations can be expressed as polynomial integrations and because the precalculation step for the coefficients of the m't fitting polynomial saves significant computational time, the m't-fitting technique can utilise 38-62% fewer numerical operations than the Gaussian quadrature integration method for second order polynomial profile analysis. We note that echoes are not obtained at every transmitted frequency for a variety of reasons (interference, absorption, etc.). The Gaussian quadrature method cannot readily utilize a precalculated step, since it requires knowledge of the start and end frequencies of each interval. However, in the m't fitting technique, we can precalculate the polynomial at every sounding frequency and calculate the integral of the product of m' and t later, at any interval or for any combination of echo frequency. We have demonstrated that the Newton-Gregory interpolating polynomial can be used to derive a polynomial that passes through a group of equispaced points. Furthermore, utilising a technique whereby three interpolating polynomials intersect, the inherent error of the derived product m't can be decreased. Applying this method to model ionospheric data, the test results indicate that an error of a few meters can be obtained in the true-height analysis of an overlapping Chapman profile.

Acknowledgements

This work has been supported by the Space Dynamics Laboratory at Utah State University. NSF grant DPP91-19382 also provided funding for one of the authors (FTB).

References

Becker, W. (1960), Tables of ordinary and extraordinary refractive indices, group refractive indices and h'o, x(f) curves for standard ionospheric layer model, in Mitteilungen aus dem Max-Planck-Institut für Aeronomie, 4, 25 pp., Springer-Verlag, New York.

Cheney, Ward and David Kincaid (1980), Numerical mathematics and computing, Brooks/Cole Publishing Company, California.

Love, Carl H. (1966), Abscissas and weights for Gaussian quadrature, National Bureau of Standards, Monograph 98.

McNamara, L. F. and J. E. Titheridge (1977), Numerical ionograms for comparing methods of N(h) analysis, IPS Series R Reports X-5.

Shinn, D. H. and H. A. Whale (1952), Group velocities and group heights from the magneto-ionic theory, J. Atmos. Terr. Phys., 2, 85-105.

Shinn, D. H. (1955), Tables of group refractive index for the ordinary ray in the ionosphere, in Physics of the Ionosphere, pp. 402-406, Physical Society, London.

Titheridge, J. E. (1985), Ionogram analysis with the generalised program POLAN, World Data Center A for Solar-Terrestrial Physics, Report UAG-93.

Titheridge, J. E. (1988), The real height analysis of ionograms: A generalized formulation, Radio Sci., 23, 831-849.

Back to INAG Homepage Back to UAG-104 Contents page