Discussion and Comparison of the Mapping Functions in Radio Frequencies

Some mapping functions in radio frequency are reviewed at the begin­ ning of this paper. From the viewpoint of statistics and mathematics, the error sources of refractive delay computation are analyzed together with the ray tracing technique and radiosonde data application. The NMF model (Niell 1996) used globally distributed radiosonde network data and tabular coefficients in its continued fraction expression and, therefore, reached a high accuracy with wide latitude and meteorological parameter coverage (MacMillan and Ma 1997). The UNSW931 model (Yan and Ping 1995) in­ troduced an improved continued fraction form, which was derived from the generator function method and had a good convergency at low eleva­ tion observations. The influences of ground meteorological records on the mapping functions are generally discussed. Comparisons among different mapping functions are performed together with radiosonde data. In order to fit the requirements of more accurate and lower elevation observations in future space technique applications, we suggest constructing a new map­ ping function, which will involve the main advantages of present models and improve or avoid inadequate expressions as far as possible. (


INTRODUCTION
In modern space techniques, atmospheric refraction was shown to be one of the most important error sources.Improvement of computational accuracy of refractive delay is one of the critical tasks for the geo-and space science researches in both the present and the future.
The refractive delay of signal transferring between source and receiver is defined as the differ ence of optical and geometric �istances: L\.cr = J nqf -J ds , , (1) e where 1 is the track of the signal, sis the geometric straight line connecting source and receiver, n is the refractive index of atmosphere.The concept of the mapping function was firstly pro posed by Marini (1972), which nowadays has become a basic formulation used in atmospheric refraction delay correction,.The refractive delay L\.cr is then formally written as: where e is the elevation angle of the relevant object, L\.cr z is defined as the zenith delay which can be accurately solved by ground meteorological records (Saastamoinen 1972), and the con tinued fraction mapping function m(e) with a set of constant coefficients a,b,c ... has the form In comparison with the series expansion method, used for instance by Saastamoinen (1972), the above continued fraction expression may reach higher accuracy and lower elevation cov erage by a simpler computation procedure.
Noticing the convergency requirement at zenith and incorporated Chao's (1972) formulation, Davis et al. (1985) replaced the second sin( e)of Eq.(3) by tan( e ).In order to describe the influences of atmospheric profiles on the mapping function, Davis et al. (1985) developed a parameterized continued fraction model, which was named the Cf A2. et al. (1972)).Another contribution made by the CfA2.2 model is that the angular argument e in Eq. ( 4) was changed from the arrival direction of the object in Marini (1972) to the true direction in the CfA2.2.The true direction can be calculated more accurately from the ephem eris of the object rather than the arrival direction, which is usually measured by a local instrument.
Instead of the simplest exponential atmospheric profile used by Marini model, the CfA2.2 model furthermore adopted.theStandard Atmosphere (Allen 1973) in its coefficient resolution.
The CfA2.2 model is, at present, still widely used in the Global Positioning System (GPS) and in other space techniques because of its noticeable accuracy.
Based on data sets from ten radiosonde sites, Herring (1992) developed another expres sion form of continued fraction mapping function for hydrostatic and wet delays, respectively: (5) In the above expression the coefficients a,b,c for both hydrostatic and wet components weie shown as a function of latitude and the height of site, as well as of the surface temperature of the station (see Eqs. ( 8) and ( 9) in Herring (1992)).With its coefficient sets, Eq. ( 5) was named the MTT mapping function.After Herring's model, Niell (1996) further used more widely distributed radiosonde net work data sets, and developed a global mapping function for atmospheric delay in radio frequencies, which was named as the NMF model.The mathematical structure of the NMF model formally inherites the normalized continued fraction of Eq. ( 5) used by the MTT model (Herring 1992).But the coefficients a,b,c of the NMF are purely related to the station geo graphic location and the Day-of-Year (DOY), and are totally independent of meteorological records: (6) where D0 is the adopted phase DOY 28; d the DOY time; </Ji the latitude of station.Parameters b and c have similar structure to a.All parameters in Eq. ( 6) are tabulated with respect to latitudes (see Tables 3 and 4 of Niell (1996) for hydrostatic and wet components, respectively).Linear interpolation is suggested for the values of a,b,c at latitude of station.(In this paper, some printing mistakes in Eqs. ( 4) and ( 5) of Niell's contribution have been carefully taken into account.) Based on the generator function theory, Yan and Ping (1995) mathematically proved that the complementary error function is an approximate generator function of atmospheric refrac tive integrals of both delay and bending, and that any expansion form of the complementary error function can be used to construct the mapping function of atmospheric refraction.Ray tracking procedures can then be used to adjust the coefficients of the mapping functions.A new improved continued fraction mapping function form is then deduced from an expression of the complementary error function (Press et al. 1991 where D; are the coefficients to be adjusted, �o is the true zenith angle of object, I is defined as the normalized effective zenith argument (see Eq. ( 12) in Yan and Ping 1995).A further de scription of coefficients D; can be found in Eq. (21) of Yan and Ping (1995) for the Standard Atmosphere, which is named as the UNSW931 model.For the Hopfield atmospheric profile (Hopfield 1969(Hopfield , 1971)), the coefficients Di of the improved continued fraction mapping func tion is listed in Eq. (23) of Yan and Ping (1995), which is named the UNSW932 model.Yan (1996) further applied the generator function method to the astronomical refraction, and de veloped the improved continued fraction to construct the mapping function of the astronomi cal refraction.The introduction of the mapping function to the astronomical refraction may improve the computational accuracy to a great extent in comparison with the conventional formulas listed in the Astronomical Almanac (2001), especially in low elevation observations (Yan 1998).In optical wavelength facilities, Yan and Wang (1999) gave a frequency-related mapping function for laser ranging application.They also discussed in detail some subsidiary corrections to the mapping function, for instance the correction of a finite distance object, which is suitable to observations of low orbit spacecraft.

ERROR ANALYSIS OF THE MAPPING FUNCTION
From Eqs. (1-3), we may realize that the computational errors of atmospheric refraction are mainly related to the following sources.The first is the selection of the mathematical expression of the map ping function.The second is the deviation of the employed atmospheric profile in the model from the actual atmospheric distribution above the stations.The first error source involves two factors: the choice of the continued fraction expression of the mapping function, and the choice of the expression for the coefficients in the continued fraction.The development of the primary modes of Chao (1972), Marini (1972) for the comprehensive models of CfA2.2, MTT, UNSW931, NMF reflects this discussion which is gradually going into details.Except for simpler models, such as Marini, the choice• of the continued fraction expression plays a role only in the observations of rather low elevations.
From simulation computation, we can find that the coefficient expression becomes one of the main problems in atmospheric refractive delay (see Section 5).There are basically two sorts of the coefficient expressions of the continued fraction mapping function.The first is the linear coefficient presentation of meteorological and geophysical parameters with respect to their nominal values, as is typically being used in the CfA2.2,UNSW931 and UNSW932 models.In comparison with the primary constant coefficients in Marini (1972), the influences from the meteorological parameters are compensated to a certain extent in these models.But, when the meteorological records are far apart from the nominal values of the model, or the observation data cover wide meteorological conditions, this sort of coefficient expression may produce rather large errors.The second is the tabular coefficients as typically used in the NMF model.This simple mathematical presentation effectively keeps substantial tolerance with respect to the change of the observational seasons and of the geographic locations of stations.It may reach wider seasonal and latitude coverage and maintain high accuracy, even in severe climate conditions near polar regions.
The Standard Atmosphere is composed of a troposphere with constant temperature lapse and a stratosphere with constant temperature (Allen 1973).It is a kind of simplified global average of earth atmosphere.This model can be used to construct the coefficients of the map ping function, for instance the CfA2.2 and UNSW931.The measurement data from radio sonde stations strongly represent local atmospheric profile characteristics, and can be used to estimate the deviations of the chosen atmosphere model from the actual local atmospheric profile.Therefore, the mapping functions constructed from radiosonde data sets involve fewer errors in the atmospheric refraction computations for sites nearby these radiosonde stations.Table 1 summarizes the main features of different mapping functions.

METEOROLOGICAL AND GEOPHYSICAL PARAMETERS
The ray tracing technique can be used to quantitatively estimate the influences of various meteorological and geophysical parameters on the mapping functions.The ray tracing tech nique is a numerical integral procedure based on Snell's law.It provides accurate results of refraction for a given theoretical atmospheric model or for radiosonde data sets.The refractive delay and bending angle can be solved together.Presently, this technique has become a basic method for planetary atmospheric research, not only in atmospheric refraction, but also in atmospheric retrieval.We developed a counter-trace integral method in our ray tracing com putation software SHAOMF (Yan 1999).The departure point of the integral is set on the receiver, and the integral will end at the top of atmosphere, or at the source of signal.The arrival elevation is taken as an input parameter of the integral.The true elevation angle of the object can be exactly computed from an integral procedure.A least-squares solution can be used to get a coefficient solution of the mapping function.
Because the dispersive effect of signals in radio frequencies is not as significant as in optical frequencies, the frequency-related term is not considered in radio frequency application.Theoretically, we know that the mapping function depends on the detailed information of the atmospheric profile above the observational site.But, the earth's atmosphere is much more complicated with respect to changing space and time.In atmospheric refraction research, we have to implement a ldnd of average atmospheric profile to describe the actual earth atmosphere, which comes from the local and the global atmospheric soundings.Linear expressions of me teorological and geophysical parameters for the coefficients of the mapping function, for in stance the CfA2.2 or UNSW93 l, can, to a certain extent, eliminate the influences from atmo spheric variation.But, they will suffer accuracy corruption in severe climate conditions when the meteorological records are largely deviate from the nominal values used in the model.
In fact, the NMF model used average meteorological parameters from radiosonde data to replace the individual meteorological values at sites, and adopted a coefficient expression related to the site geographic position and the observational date.This approach still involves some remnant errors, because the meteorological and atmospheric parameters are quite differ ent in the same DOY or for different areas in the same latitude zone.
From recent knowledge of atmospheric sciences, it is clear that the earth's atmosphere above the tropopause is rather stable, but the part below the tropopause is in largely changing with space and time.The main part of water vapor in the atmosphere is below 8km.The low layer atmosphere influences mainly the climate variation globally and locally.In the same latitude zone the physical parameters of the earth atmosphere, such as temperature, pressure and humidity, are not only presented as the time functions of the height, but they also signifi cantly depend on the geographic location (sea or land, inland or coastal), geomorphologic landform (plane or mountain, forest or desert), weather condition (rain, snow or sunny), global circulation (warm or cold), daytime (day or night), human activities (which recently became more and more significant), and so on (Lamb 1972(Lamb , 1977)).From global meteorological records, 15 degree deviation of the ground temperature is normal for different sites in the same latitude and on the same DOY.Larger deviations are also possible under extreme weather conditions.20 mbar surface pressure deviation can be found in summer and autumn hurricane season records.Especially geographic location of stations and weather variation can cause significant differences in humidity, which produce noticeable influence on both the hydrostatic and the wet components of the mapping function.
Differentiating Eq. ( 2), we have: We may write the relation between the mapping function and the meteorological and geophysical parameters deviations: where p is the meteorological and geophysical parameter vector taken into consideration.
With the ray tracing technique, we can estimate the second term of Eq. ( 8) which comes from the variations of surface meteorological records.Numerical computations show that the ground temperature deviations may produce more significant influences on the mapping func tion than ground pressure deviations.This might be the reason why Herring (1992) adds a temperature-related correction to the coefficients of his mapping function MTT.

HEIGHT CORRECTION OF STATION IN THE MAPPING FUNCTION
In atmospheric refraction computation, the gravitational acceleration correction of the location of station was primarily discussed by Saastamoinen (1972).The station height and latitude correction terms were involved through adding a multiplicative factor w-1( </J ,h0) W(</J, h0) = 1-0.00266cos2¢-0.00028ho, (10) to the expression of the zenith delay L).a z of Eq. ( 2).Herring (1992) considered the influence of geographic location of site to the mapping function.He introduced a station height term to the coefficient expression of the MTT model Eq. ( 5).But, a rather general theoretical discus sion on the influences of station height on the mapping function has not yet been found or referenced.It might become a new topic in mapping function research under higher accuracy requirement of atmospheric delay in future.In the NMF model empirical formulas of station height correction were adopted, see Eqs. ( 6) and ( 7) ofNiell (1996).Simulation computations proved that these empirical formulas might involve rather larger errors, because of their inad equate mathematical expression and related dimension errors inside.Mathematically, there are obviously rather big differences between the height correction of the NMF model and the direct derivative of Eqs. ( 2) and ( 5).
By the ray tracing technique, we can check the height correction of the NMF model.The right of Fig. 1 gives the results of the NMF model, and the left of Fig. 1 originates from the ray tracing technique which might be considered as the true values of the height correction.A numerical comparison between both sides of Fig. 1 demonstrates a rather large incompatibility in the NMF height correction.We believe that the height correction of the NMF mapping function might be reconsidered either by theoretical deduction or by an appropriate empirical expression.Perhaps, the expression of MTT model (Herring 1992) could also be taken into account.

COMPARISON WITH RADIOSONDE DATA
Today there are about 1000 radiosonde stations distributed over the Earth surface, and two observations per day are carried out at every station at QhGT and 12hGT, respectively.By the ray tracing technique the atmospheric profiles of temperature, pressure and humidity ob tained from radiosonde measurements can be used not only to construct the mapping function (Herring 1992;Niell 1996), but also to estimate the accuracy of the different mapping functions.We select two radiosonde data sets coming from significant different geographic zones: Fairbanks, Alaska (Station 70261) during February 1-28, 1990 (in total 39 groups), and West Palm Beach, Florida (Station 72203) February 10-28, 1990 (in total 29 groups).Figs. 2 and 3 illustrate the comparisons of the different mapping functions with respect to the ray tracing results of the first several days' measurements from these two stations.The results obtained from other day data have a similar tendency.In comparison with radiosonde data, Fig. 3 shows a general closer tendency for the NMF than the CfA2.2,MTT, UNSW931 at Station 70261.But Fig. 2 displays that the NMF does not dominate the UNSW931 or MIT models at Station 72203.
The CfA2.2 and UNSW93 l model used the same atmospheric profile and the same coef ficient expression, but different continued fraction forms in the mapping function.From the comparisons of the CfA2.2 and UNSW931 in Figs. 2 and 3, we may formulate that: (a) The UNSW93 l is better than the CfA2.2 in 72203 station; that might be explained by the adoption of the new continued fraction by the UNSW93 1. But, in 70261 the CfA2.2 is slightly better.There might be two reasons.The first is that the UNSW93 l gives rather large weight to the nominal point.The constant terms of the coefficient in the UNSW931 are directly taken from the nominal point.The errors of the UNSW931 will become more significant than the CfA2.2 while the ground meteorological records run far from the nomi nal point.The second is that, in coefficient resolution procedure of the UNSW93 l, the temperature coverage was taken from +35 °C to -20 °C, which does not cover the tempera ture range of station 70261.
(b) In station 7026 1 the results of the UNSW93 l still present a convergency tendency when elevation is going down to horizon.That means UNSW931 has a better convergency than other in very low elevation coverage.
Some radiosonde stations were built near the main space observational instruments.They offer a possibility to construct a local or quasi-local mapping function with high accuracy on top of the instruments.A local mapping function from radiosonde data can still keep high accuracy under severe climate conditions.For instance, the radiosonde data of Station 70261 are taken from typical harsh winter weather samples (the ground temperature was between -30 ° and -40 °C.The geophysical parameters of this station greatly deviated from the nominal values adopted in the UNSW931 and CfA2.2.In this case, the NMF model is better.The reasons are: (1) the NMF model uses tabular coefficient forms in its mapping function; (2) the NMF model involves the radiosonde data of Fairbanks in coefficient resolution procedure.
Because there are many geophysical factors which may influence the atmospheric pro files globally or locally (see Section 3), it will be very difficult to find out a simple mathemati cal expression for a global mapping function by radiosonde data, no matter how accurate are the radiosonde data, and how much radiosonde station data are used in resolution.As an example, the data of Station 72203 are in a normal weather condition, and the ground meteorological records are close to the nominal values of the UNSW93 l model.the UNSW93 l has a better coincident result than the NMF model, especially at low elevation.� 0 F--.=!'"'F=l""=!="il'l=i�t+-+.-:-�.... .. .. This expression of the continued fraction may help us to compare the MTT and NMF models with the other mapping functions.We believe, there is no reason to consider the ex pression of Eq. ( 12) superior than Eqs.(4) or (7).
From Fig. 2, we can find that if the observation elevation angles are constrained above 5°C, the choice of the con tinued fraction of the mapping function is not a critical problem.
When the elevation is going below 5°C, the improved continued fraction obtained from the generator function method might have better con vergency, because it is based on a spherical earth atmospheric approach (Yan and Ping 1995).From Fig. 3, we may state that the structure of the coefficients of the mapping fun ction is a dominant problem within a wide coverage of climate conditions.The tabular parameter form of the NMF (Niell 1996) is a simple and an effective way to eliminate the nonlin ear characteristics of the coefficients of the mapp ing function in a very wide variable range of the meteorological and the geophysical factors of earth atmosphere.
Tabular parameter form and application of widely distributed radio network data make the NMF model to be, in our opinion, a kind of quasi-local mapping function, not a global one.
The high accuracy and wide coverage of meteorological conditions of the NMF model mainly come from its tabular parameter form.Because of the variability of earth atmosphere with time and sp ace, it mig ht be efficient to introduce a ground temperature related correction term to the NMF model.The second possible improvement to the NMF model might be the recon sideration of its height correction formula.For an expression of empirical formulas similar to Eqs. (8) and (9) of Niell ( 1996), a dimensionless normalized height argument is strongly required.
Rigorous deduction of the station height correction of the mapping function is a very difficult mathematical problem.A doubtful choice of station height correction will fin ally con taminate the accuracy of the mapping function (see Fig. 1).Perhaps the correction method of station height in the MTT model (Herring 1992) can be further studied.
In the list of names of the radiosonde stations used in the NMF model, we find Fairbanks.
But the station of West Palm Beach is not yet included.In the latitude belt of Fairbanks there are only few radiosonde stations involved in the NMF.Therefore, the radiosonde data from the Fairbanks station have noticeable weights.In the middle latitude belt, in which West Palm Beach is located, there are much more radioson de station s included in the NMF model.From Fig. 2, the results of the NMF in West Palm Beach station, which can be recognized as the average of other radiosonde stations in this latitude belt, are not significantly over other map ping functions.The first possible reason of this fact might be the inappropriate height correc tion formulation of the NMF model.The second reason might be that the ground meteorologi cal records of West Palm Beach are near the nominal values of the meteorological parameters in the CfA2.2 and UNSW93 l models.And the third might be that the NMF does not consider variations of meteorological parameters, especially ground temperature.

DISCUSSION
In modem space techniques, low elevation observations are required in some observa tional versions.VLBI observations have reached 3° or lower elevation.Observations below 10° in SLR are considered in a "synchronous" mode.As a special application of the atmo spheric retrieval, a mountain-based or airborne GPS receiver involves the low observations near and below horizon (Zuffada et al. 1999).As analyzed in previous sections, we think the tabular coefficient presentation used by NMF is the best form for further new mapping func tion coefficient construction.Radiosonde data can offer most accurate atmospheric profile above stations.The new continued fraction form of Eq. ( 7) has a relatively good convergent tendency of convergency at low elevation.By elevation cut-off tests of VLBI observations, Ping et al. (1997) found that the UNSW931 model performed good repetitions of base-line measurements.For PRARE measurement data, Zhang et al. (2001) give results that the UNSW931 and NMF models can reach almost the same accuracy, and they are significantly better than other mapping functions.In order to improve the accuracy of the mapping function in low elevation, we also suggest adding more ray tracing integrals below 10° in coefficient resolution.Nine elevation integrals from 3° to 90° used by Marini, MTT and NMF might be not enough.
) === ______ a __ _ sine+ -----b- tan e + --sine+ c(4) in which the coefficients a, b, c are presented as the linear function of the total surface pressure P0, the partial pressure of water vapor at the surface e0, the surface temperature T0, the tropo spheric temperature lapse rate /3 and the height of tropopause h (see Eqs.(10-12) of Davist.

Fig. 1 .
Fig. 1.Site height correction with respect to sea level.

Table 1 .
Comparisons of different mapping functions.