# Atmospheric Sounding from Fengyun-3C GPS Radio Occultation Observations: First Results and Validation

Authors: Shuanggen Jin; Chao Gao; Junhai Li
Journal: Advances in Meteorology
 (2019)

---

## Abstract

Carrying global positioning system (GPS) radio occultation (RO) receiver, Chinese meteorological satellite Fengyun-3C (FY-3C) was launched on September 23, 2013, which provides new observation data for observations and studies of weather and climate change. In this paper, the results of FY-3C GPS RO atmospheric sounding are presented for the first time, including high-order ionospheric correction, atmospheric parameters estimation, and evaluation by COSMIC and radiosonde observations as well as applications in estimating gravity wave activities. It is found that the effect of the ionospheric correction residual on the phase delay is below 20 mm, which has minimal impact on bending angle estimation and generates differences of about 1 K in the average temperature profile. The difference between FY-3C and COSMIC temperatures at all heights is within 1°C, and the tropopause temperature and height have a good consistency. Deviations from Radiosonde measurements are within 2°C, and the tropopause temperature and height results also have a strong consistency. Furthermore, global gravity wave potential energy is estimated from FY-3C GPS RO, exhibiting similar behavior to results derived from COSMIC radio occultation measurements. The mean value of the gravity wave potential energy near the equator is 10 J/kg and decreases toward the two poles while in the northern hemisphere, it is stronger than that in the southern hemisphere.


---

## Body

## 1. Introduction

GPS radio occultation measurements can provide abundant atmospheric parameters like the atmospheric refractive index, pressure, and temperature, which can be used to study the atmospheric structure, variations, and dynamics of the Earth. In 1995, the low-orbit satellite Microlab1 with GPS receiver was launched, and its measurement results showed that the Earth’s atmospheric parameters could be inverted by receiving GPS radio occultation signals, and the accuracy of the retrieving temperature is about 1 K [1, 2]. Subsequently, a variety of low-orbit satellites carrying Global Navigation Satellite System (GNSS) radio occultation receivers have been widely used for meteorology and atmospheric studies. For example, the satellite CHAMP (CHAllenging Minisatellite Payload) launched in July 2000 provided GNSS occultation data as well as Argentinian satellite SAC-C (Satellite de Aplicaciones Cientificas-C) and US/German satellite GRACE (Gravity Recovery and Climate Experiment), and the deviation of the inverted temperature profiles is less than 1 K [3, 4]. In addition, Hajj and Romans pointed out that GPS/MET can be used for analysis of ionospheric E and F layers’ electron density variations [5]. Jakowski et al. obtained the electron density profile results of the ionosphere from CHAMP’s radio occultation data [6]. In 2006, the United States and Taiwan, China, jointly developed the COSMIC (Constellation Observing System for Meteorology Ionosphere and Climate) mission with six small satellites, which can provide more than 2,500 occultation events per day. This greatly enhances the geographic coverage of the Earth’s atmosphere, which also provides rich data for the ionospheric and tropospheric applications [7–14].On September 23, 2013, Chinese meteorological satellite FY-3C with a track altitude of 836 km and an orbital inclination of 98.75° was launched successfully. FY-3C is the second generation of China’s polar-orbit meteorological satellites and belongs to the near-polar solar synchronous orbit. It is equipped with a variety of payloads, such as the space environment detector, the solar radiation detector, the Earth radiation detector, the microwave radiometer, the microwave thermometer, the infrared spectrometer, and GNSS radio occultation sounder. Therefore, this new FY-3C satellite will further improve the spatial and temporal resolution of the GNSS radio occultation dataset and provide abundant data for studies of the ionosphere, troposphere, and climate change. In this paper, the first atmospheric results are presented from FY-3C GPS RO data, including high-order ionospheric effects, atmospheric temperature and pressure, and evaluation with COSMIC and radiosonde observations as well as applications in estimating gravity wave potential energy.


## 2. Data Processing and Methods

### 2.1. FY-3C GPS RO Data and Radiosonde Data

The radio occultation receiver GNOS (GNSS occultation sounder) carried on FY-3C is the first satellite occultation receiver of China. It can simultaneously receive the navigation signals of both GPS and BeiDou satellites for occultation observations. The GNOS occultation receiver is equipped with the latest open-loop technology, which can reach 100 Hz in neutral atmospheric occultation detection, and the tracked signal can extend down to 1 to 2 km above the surface [15, 16]. FY-3C can provide GPS occultation events about 500 times per day. The GPS occultation data were released in August 2014, which can be obtained from the Fengyun meteorological satellites’ official website (http://satellite.nsmc.org.cn), including the occurring time and duration of occultation events as well as dual-frequency atmospheric phase data. Detailed precise orbits for the FY-3C spacecraft and corresponding attitude data are described in Li et al. [17].The radiosonde data are from the IGRA (Integrated Global Radiosonde Archive) project of NOAA (National Oceanic and Atmospheric Administration). The project provided more than 2,700 radiosonde stations around the world from 1905 to present. The data contain the location of stations, the observation time, the sounding level above stations, the temperature, the atmospheric pressure, etc. Figure1 shows FY-3C occultation events and co-located radiosonde stations from September 21th to 28th, 2014.Figure 1
FY-3C occultation events and co-located radiosonde stations from September 21th to 28th, 2014.


### 2.2. GPS Occultation Data Processing

Firstly, the occultation bending angle profile parameters are obtained by Doppler shift or additional phase method. Secondly, the bending angle profile is transformed into the refractive index profile by the Abel transform. Using the relationship between the atmospheric refractive index and the temperature as well as the pressure, the neutral atmospheric temperature and pressure can be calculated. The most commonly used method to obtain the occultation bending angle profile through LEO occultation Doppler shift or additional phase data is geometrical optics or wave optical method.The relationship between atmospheric Doppler shiftΔf and the motions of GNSS satellite and LEO satellite can be described by the following formula [1]:(1)λΔf=vt⇀·kt⇀−vr⇀·kr⇀−vt⇀−vr⇀⋅k⇀,where λ is the signal wavelength; vt⇀ and vr⇀ are the velocity vectors of GNSS satellite and LEO satellite, respectively; and k⇀, kt⇀ and kr⇀ are the unit vector of GNSS satellite to LEO satellite, the unit vector of GNSS satellite signal broadcasting, and the unit vector of LEO satellite signal received, respectively.As can be seen from Figure2, if the occultation event is limited in the occultation plane, the bending angle α is the sum of the angles δt and δr. Then, we assume that the refractive indexes near the tangent point are the same value and follow the spherical symmetry, δt and δr also obey the shell rule [18], namely,(2)a=nrtrtsinθt+δt=nrrrrsinθr+δr,where a is the impact parameter; nrt and nrr are the refractive indexes of LEO satellite and GNSS satellite, respectively; rt and rr are the orbital radius of LEO satellite and GNSS satellite, respectively; θt is the angle between the unit vector of GNSS-LEO line of sign and GNSS satellite orbit vector; and θr is the angle between the unit vector of GNSS-LEO line of sign and LEO satellite orbit vector.Figure 2
Schematic diagram for inversion of atmospheric bending angle using geometrical optics.Under the assumption that the refractive index exhibits stable stratification near the tangent point, we can get the relation between the bending angleα and the refractive index [1]:(3)αa=−2a∫a∞1x2−a2dlnndxdx,where n is the refractive index, and x is the integral variable of the impact parameter. Through the Abel transformation, we can obtain the relation between the refractive index n and the bending angle. In the neutral atmosphere, the refractive index is a function of the atmospheric temperature, pressure, and water vapor. So by putting certain restrictions and conditions on the refractive index, we can obtain the neutral atmospheric parameters. Nowadays, there are several different RO processing systems [19], and here, we used the Radio Occultation Processing Package (ROPP).


### 2.3. High-Order Correction Algorithm

When the GNSS signal propagates in the ionosphere and the neutral atmosphere, the phase delay caused by the atmospheric refraction can be expressed as(4)ΔL=L−ρ=∫nneu+niondl−ρ,where ΔL is the phase delay, L is the phase measurement, ρ is the geometric distance between GNSS satellite and LEO satellite, nneu and nion are the refractive indexes in the neutral atmosphere and the ionosphere, respectively. And the ionospheric refractive index nion is associated with the electron density and the magnetic field intensity of the GNSS signal propagation path, which can be described by Appleton–Hartree formula [20]:(5)nion2=1−2X1−X21−X−Yt2±Yt4+4Yl21−X2,(6)X=e24π2ε0mNef,Y=e2πmB0f,Yt=Ysinθ,Yl=Ycosθ,where B0 is the magnetic field intensity, Ne is the electron density, f is the frequency of the electromagnetic wave, m is the electron mass, ε0 is the dielectric constant, e is the electronic charge, and θ is the angle between the spread direction of the signal and the direction of the geomagnetism. Simplifying equation (5), we get the following:(7)nion=1−12X−12XYcosθ−14X12X+Y21+cosθ.The ionospheric delay can be further expressed as [21](8)Iϕ=∫niondl=−qf2−s2f3−r3f4.The quantities in equation (8) are calculated as follows:(9)q=Cx2∫Nedl,s=CxCy∫NeB0cosθdl,r=r1+r2+r3,r1=3Cx28∫Ne2dl,r2=3CxCy22∫NeB02cos2θdl,r3=3CxCy24∫NeB02dl,Cx=e24π2ε0m,Cy=e2πm.The second term on the right side of equation (8) is the second-order ionospheric delay, and the third term is the third-order ionospheric delay. Assessing and correcting the effect of the ionospheric correction residual on an ionospheric or atmospheric occultation involves solving the second-order and third-order ionospheric delays (the fourth-order term is very small and is therefore negligible).To estimate the high-order ionospheric delay of the GNSS signal, the path information between the signal’s emission position and the received position is needed. Born and Wolf proposed a signal propagation model in the nonvacuum space—the optical path formula [18]:(10)ddsndr→ds=nd2r→ds2+dndsdr→ds=∇n,where r→ is the Cartesian coordinate of the signal path, s is the length of the signal path, n is the refractive index, and ∇n is the gradient component of the refractive index in the Cartesian coordinate system. Then, we can get(11)d2r→ds2=1n∇n−dndsdr→ds,where r→ contains three components x,y,z, and if the refractive index of each point and the gradient of the refractive index ∇n are known, we can use the fourth-order Runge–Kutta numerical integration method to numerically integrate each component and find the specific coordinates of each point in the signal path.By the refractive index formula, the refractive index can be obtained according to the meteorological elements of each point in the signal path. In order to carry out the ray tracing, we choose the neutral atmospheric model MSIS90 to obtain neutral atmospheric parameters such as the atmospheric pressure and temperature. The MSIS90 model integrates the rocket mass spectrometer, the ground noncoherent scattering radar, and satellite data to fit the global temperature and material density model at different locations, altitudes, and times [22]. The NeQuick2 model is used to obtain the ionospheric electron density information under different solar activity cycles [23]. The NeQuick2 is an ionospheric model developed by ICTP (International Center for Theoretical Physics). In addition, the magnetic field strength and three-dimensional components of arbitrary point on the global are used from IGRF12 (International Geomagnetic Reference Field, 12th generation), which uses the mass data of various observation methods to model the global geomagnetic data [24]. After the pressure, temperature, electron density, geomagnetic intensity, and other parameters are obtained through the models, we can use numerical integration to obtain the GNSS signal path according to equation (10) and then obtain the values of the second-order and the third-order ionospheric delays from the GNSS signal path according to equation (8).


## 2.1. FY-3C GPS RO Data and Radiosonde Data

The radio occultation receiver GNOS (GNSS occultation sounder) carried on FY-3C is the first satellite occultation receiver of China. It can simultaneously receive the navigation signals of both GPS and BeiDou satellites for occultation observations. The GNOS occultation receiver is equipped with the latest open-loop technology, which can reach 100 Hz in neutral atmospheric occultation detection, and the tracked signal can extend down to 1 to 2 km above the surface [15, 16]. FY-3C can provide GPS occultation events about 500 times per day. The GPS occultation data were released in August 2014, which can be obtained from the Fengyun meteorological satellites’ official website (http://satellite.nsmc.org.cn), including the occurring time and duration of occultation events as well as dual-frequency atmospheric phase data. Detailed precise orbits for the FY-3C spacecraft and corresponding attitude data are described in Li et al. [17].The radiosonde data are from the IGRA (Integrated Global Radiosonde Archive) project of NOAA (National Oceanic and Atmospheric Administration). The project provided more than 2,700 radiosonde stations around the world from 1905 to present. The data contain the location of stations, the observation time, the sounding level above stations, the temperature, the atmospheric pressure, etc. Figure1 shows FY-3C occultation events and co-located radiosonde stations from September 21th to 28th, 2014.Figure 1
FY-3C occultation events and co-located radiosonde stations from September 21th to 28th, 2014.


## 2.2. GPS Occultation Data Processing

Firstly, the occultation bending angle profile parameters are obtained by Doppler shift or additional phase method. Secondly, the bending angle profile is transformed into the refractive index profile by the Abel transform. Using the relationship between the atmospheric refractive index and the temperature as well as the pressure, the neutral atmospheric temperature and pressure can be calculated. The most commonly used method to obtain the occultation bending angle profile through LEO occultation Doppler shift or additional phase data is geometrical optics or wave optical method.The relationship between atmospheric Doppler shiftΔf and the motions of GNSS satellite and LEO satellite can be described by the following formula [1]:(1)λΔf=vt⇀·kt⇀−vr⇀·kr⇀−vt⇀−vr⇀⋅k⇀,where λ is the signal wavelength; vt⇀ and vr⇀ are the velocity vectors of GNSS satellite and LEO satellite, respectively; and k⇀, kt⇀ and kr⇀ are the unit vector of GNSS satellite to LEO satellite, the unit vector of GNSS satellite signal broadcasting, and the unit vector of LEO satellite signal received, respectively.As can be seen from Figure2, if the occultation event is limited in the occultation plane, the bending angle α is the sum of the angles δt and δr. Then, we assume that the refractive indexes near the tangent point are the same value and follow the spherical symmetry, δt and δr also obey the shell rule [18], namely,(2)a=nrtrtsinθt+δt=nrrrrsinθr+δr,where a is the impact parameter; nrt and nrr are the refractive indexes of LEO satellite and GNSS satellite, respectively; rt and rr are the orbital radius of LEO satellite and GNSS satellite, respectively; θt is the angle between the unit vector of GNSS-LEO line of sign and GNSS satellite orbit vector; and θr is the angle between the unit vector of GNSS-LEO line of sign and LEO satellite orbit vector.Figure 2
Schematic diagram for inversion of atmospheric bending angle using geometrical optics.Under the assumption that the refractive index exhibits stable stratification near the tangent point, we can get the relation between the bending angleα and the refractive index [1]:(3)αa=−2a∫a∞1x2−a2dlnndxdx,where n is the refractive index, and x is the integral variable of the impact parameter. Through the Abel transformation, we can obtain the relation between the refractive index n and the bending angle. In the neutral atmosphere, the refractive index is a function of the atmospheric temperature, pressure, and water vapor. So by putting certain restrictions and conditions on the refractive index, we can obtain the neutral atmospheric parameters. Nowadays, there are several different RO processing systems [19], and here, we used the Radio Occultation Processing Package (ROPP).


## 2.3. High-Order Correction Algorithm

When the GNSS signal propagates in the ionosphere and the neutral atmosphere, the phase delay caused by the atmospheric refraction can be expressed as(4)ΔL=L−ρ=∫nneu+niondl−ρ,where ΔL is the phase delay, L is the phase measurement, ρ is the geometric distance between GNSS satellite and LEO satellite, nneu and nion are the refractive indexes in the neutral atmosphere and the ionosphere, respectively. And the ionospheric refractive index nion is associated with the electron density and the magnetic field intensity of the GNSS signal propagation path, which can be described by Appleton–Hartree formula [20]:(5)nion2=1−2X1−X21−X−Yt2±Yt4+4Yl21−X2,(6)X=e24π2ε0mNef,Y=e2πmB0f,Yt=Ysinθ,Yl=Ycosθ,where B0 is the magnetic field intensity, Ne is the electron density, f is the frequency of the electromagnetic wave, m is the electron mass, ε0 is the dielectric constant, e is the electronic charge, and θ is the angle between the spread direction of the signal and the direction of the geomagnetism. Simplifying equation (5), we get the following:(7)nion=1−12X−12XYcosθ−14X12X+Y21+cosθ.The ionospheric delay can be further expressed as [21](8)Iϕ=∫niondl=−qf2−s2f3−r3f4.The quantities in equation (8) are calculated as follows:(9)q=Cx2∫Nedl,s=CxCy∫NeB0cosθdl,r=r1+r2+r3,r1=3Cx28∫Ne2dl,r2=3CxCy22∫NeB02cos2θdl,r3=3CxCy24∫NeB02dl,Cx=e24π2ε0m,Cy=e2πm.The second term on the right side of equation (8) is the second-order ionospheric delay, and the third term is the third-order ionospheric delay. Assessing and correcting the effect of the ionospheric correction residual on an ionospheric or atmospheric occultation involves solving the second-order and third-order ionospheric delays (the fourth-order term is very small and is therefore negligible).To estimate the high-order ionospheric delay of the GNSS signal, the path information between the signal’s emission position and the received position is needed. Born and Wolf proposed a signal propagation model in the nonvacuum space—the optical path formula [18]:(10)ddsndr→ds=nd2r→ds2+dndsdr→ds=∇n,where r→ is the Cartesian coordinate of the signal path, s is the length of the signal path, n is the refractive index, and ∇n is the gradient component of the refractive index in the Cartesian coordinate system. Then, we can get(11)d2r→ds2=1n∇n−dndsdr→ds,where r→ contains three components x,y,z, and if the refractive index of each point and the gradient of the refractive index ∇n are known, we can use the fourth-order Runge–Kutta numerical integration method to numerically integrate each component and find the specific coordinates of each point in the signal path.By the refractive index formula, the refractive index can be obtained according to the meteorological elements of each point in the signal path. In order to carry out the ray tracing, we choose the neutral atmospheric model MSIS90 to obtain neutral atmospheric parameters such as the atmospheric pressure and temperature. The MSIS90 model integrates the rocket mass spectrometer, the ground noncoherent scattering radar, and satellite data to fit the global temperature and material density model at different locations, altitudes, and times [22]. The NeQuick2 model is used to obtain the ionospheric electron density information under different solar activity cycles [23]. The NeQuick2 is an ionospheric model developed by ICTP (International Center for Theoretical Physics). In addition, the magnetic field strength and three-dimensional components of arbitrary point on the global are used from IGRF12 (International Geomagnetic Reference Field, 12th generation), which uses the mass data of various observation methods to model the global geomagnetic data [24]. After the pressure, temperature, electron density, geomagnetic intensity, and other parameters are obtained through the models, we can use numerical integration to obtain the GNSS signal path according to equation (10) and then obtain the values of the second-order and the third-order ionospheric delays from the GNSS signal path according to equation (8).


## 3. Results and Discussion

### 3.1. Effect of High-Order Ionospheric Delays

In the estimation of the ionospheric correction residual on occultation events, we mainly account for the second-order and the third-order ionospheric delays, and the ionospheric correction residuals. We evaluated all the second-order ionospheric delays of FY-3C satellite neutral atmospheric occultation events in September 2014. In the vicinity of the surface, the second-order delays are within ±15 mm. As the heights increase, the second-order delays also increase. At altitudes above 90 km, the second-order ionospheric delays are about ±20 mm, and the average of the second-order ionospheric delays at all heights is about ±7 mm. The second-order ionospheric delays at 15 km, 45 km, and 90 km are shown in Figure3. At the height of 15 km, almost all of the second-order delays are less than 15 mm; at the height of 45 km, there are less than 5% of the second-order delays greater than 15 mm; and when the height rises to 90 km, the ratio of the second-order delays over 15 mm is about 10%.Figure 3
The statistics histograms of the second-order ionospheric delays at 15 km (a), 45 km (b), and 90 km (c).
(a)
(b)
(c)The residual effect of ionospheric delay after using the dual-frequency correction is called the ionospheric correction residual, and the statistic histograms are shown in Figure4. At the height of 15 km, 98% of the residuals are less than ±10 mm, at the height of 45 km, there are 95% of the residuals less than ±10 mm, but the ratio of the second-order ionospheric correction residuals within ±10 mm is less than 70% at the height of 90 km, which indicates that the effect of the ionospheric correction residual on the upper atmosphere is more pronounced. Since the top-to-bottom method is always used when we invert the refraction rate and the temperature, the ionospheric correction residual of the upper atmosphere will affect the accuracy of the lower atmosphere. Therefore, we usually limit the maximum height of inversion to weaken the effect of the upper atmosphere in the data processing.Figure 4
The statistic histograms of the second-order ionospheric correction residuals at 15 km (a), 45 km (b), and 90 km (c).
(a)
(b)
(c)In addition, the third-order ionospheric delays in the measurement of neutral atmospheric occultation are also estimated, and the values are very small. The maximum of which is 0.5 mm and the average is only about 0.15 mm, which is only 1/40 to 1/30 of the second-order ionospheric delay. Therefore, the second-order ionospheric delay is the main source of the ionospheric correction residual.We further estimated the effect of the ionospheric correction residual on the bending angle and temperature of the neutral atmosphere. Here, the second-order ionospheric delay is directly eliminated from the original observation, and then the inversion is performed according to the ionospheric-free linear combination of both GPS frequencies to obtain the bending angle, the refractive index, and the temperature of the atmosphere, which are then compared with the parameters obtained from the original data without eliminating the high-order ionospheric delay.In order to verify the feasibility of correcting the high-order ionospheric term by directly subtracting the second-order ionospheric delay from the original observation data, we compared the results from colocated occultation and radiosonde measurements for 9 different radiosonde stations (Figure5). It was found that most temperature profiles obtained through the correction of high-order ionospheric delays are closer to the results of radiosonde, and the temperature measured by radiosonde is not affected by the ionosphere. Therefore, it is feasible to use this method to effectively eliminate the effect of the high-order ionospheric delays.Figure 5
Effect of the ionospheric correction residual on the neutral atmospheric temperature estimation at 9 radiosonde stations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Furthermore, 7,144 occultation events were estimated with and without eliminating the high-order ionospheric delays, which were divided into three groups according to the maximum residuals of temperature profiles between 2 K, 4 K, and 6 K. There are 6,508 events located within 2 K, accounting for 91.1% of the total sample, 375 residuals of the temperature are between 2 K and 4 K, accounting for 5.2%, and the residuals of the temperature between 4 K and 6 K are 261 times, accounting for 3.7% of the estimated total sample. And the average temperature residual of ionospheric correction at all altitudes is about 1 K.


### 3.2. Comparison with COSMIC and Radiosonde

#### 3.2.1. Comparison between COSMIC and FY-3C Data

In order to verify the validity of FY-3C occultation data, we obtained COSMIC occultation data and looked for the overlapping occultation events from COSMIC and FY-3C within 3° in latitude, 3° in longitude, and 3 hours. There are totally 255 overlapping events from September 21th to 28th, 2014. The profiles of four overlapping occultation events are shown in Figure6. The inversion results of the temperature profiles from FY-3C occultation data are similar to those from COSMIC occultation data, but the temperature profiles of FY-3C are slightly smoother than those of COSMIC.Figure 6
Comparison of the temperature profiles obtained from FY-3C and COSMIC occultation data.
(a)
(b)
(c)
(d)From the statistics of the difference and the standard deviation of the temperature profiles at the height of 10 km, 15 km, 20 km, 25 km, 30 km, and 35 km (Figure7), we find that the average differences between the FY-3C and COSMIC temperature profiles are within 1°C, and the standard deviation is about 2 to 3°C. In the refractive index curve, we can also see the consistency of the refractive index obtained by FY-3C and COSMIC occultation data.Figure 7
Statistics of temperature profiles and the relative refraction index obtained from FY-3C and COSMIC occultation data.
(a)
(b)In addition, we obtained and compared the tropopause height and tropopause temperature from COSMIC and FY-3C occultation data, respectively (Figure8). It was found that there is a strong correlation of the tropopause height between FY-3C data and COSMIC data, and the correlation coefficient is 0.94, indicating that they have similar capability to retrieve the tropopause height. The tropopause temperatures from FY-3C data and COSMIC data also have a very strong correlation, with a correlation coefficient of 0.96.Figure 8
Comparison of the tropopause height (a) and temperature (b) obtained from FY-3C and COSMIC occultation data.
(a)
(b)


#### 3.2.2. Comparison between Radiosonde and FY-3C Data

In order to verify the accuracy of FY-3C occultation data, the independent radiosonde data are further used. We chose radiosonde data from the IGRA program for the experiment, and the criteria for co-located station observations is also set as 3° in latitude, 3° in longitude, and 3 hours. Figure9 shows the temperature profiles of four repeated observations. We find that the temperature profiles are very similar in the region below 30 km, which indicates that FY-3C occultation data have a strong reliability.Figure 9
Comparison of the temperature profiles obtained from FY-3C GPS radio occultation data and radiosonde data.
(a)
(b)
(c)
(d)From the statistics of the difference and the standard deviation of the temperature profiles at the heights of 10 km, 15 km, 20 km, 25 km, and 30 km (Figure10), we find that the difference of the temperature between FY-3C and radiosonde is within 2°C, and the standard deviation is about 3°C. The difference of the refractive index calculated by radiosonde and FY-3C is very small with a relative error of less than 0.5%, and the standard deviation is less than 1.4%, which indicates that FY-3C occultation data have the same accuracy as radiosonde data in terms of statistics.Figure 10
Statistics of the temperature profiles and the relative refraction index obtained from FY-3C occultation data and radiosonde data.
(a)
(b)In addition, we further compared the tropopause height and tropopause temperature obtained from radiosonde and FY-3C (Figure11). The correlation coefficient of the tropopause height is 0.93, and the correlation coefficient of the tropopause temperature is 0.95. It indicates that the tropopause height and temperature estimated from FY-3C are equivalent to that from radiosonde, which is consistent with the conclusion of Liao et al. [25].Figure 11
Comparison of the tropopause height (a) and temperature (b) obtained from FY-3C occultation data and radiosonde data.
(a)
(b)


### 3.3. Gravity Wave Potential Energy from FY-3C GPS Radio Occultation

Considering the volume and spatial and temporal distribution of FY-3C occultation data, the Earth is meshed into a 90 × 90 × 1 grid with a spatial resolution of 2° in latitude, 4° in longitude, and a temporal resolution of 1 month. The temperature profile is interpolated with a resolution of 200 m, and the average is taken as the value of the grid node.First, we performed S-transform on each height of every latitude circle and extracted the components whose wavenumber are of 0–6 to eliminate the interferences of some zonal disturbances, like the planetary wave [26, 27]. The method to separate disturbance temperature field from the background temperature field is to conduct continuous wavelet transform and then use the Butterworth filter with the interception frequency of 0.1–0.5 × 103 MHz for bandpass filtering, and finally, the inverse transformation is performed to get the background temperature field.The atmospheric gravity wave is produced by tiny air mass through reciprocating vibration under the influence of gravity and buoyancy. We call this oscillation frequency buoyant frequency (Brumt–Vaisala frequency):(12)Nz2=gT¯∂T¯∂z+gcp,where Nz is the buoyant frequency, g is the gravitational acceleration, T¯ is the background temperature field, and cp is the heat capacity at constant pressure. The gravity wave potential energy can be described as(13)Epz=12gNz2T′T¯2,where Ep is the gravity wave potential energy and T′ is the gravity wave disturbance temperature. As for the processing mode of T′2, we adopted the method of Tsuda et al. and take 2 km as the window and 200 m as the step size to calculate the smoothing values in the window [28]:(14)T′2=1z2−z1∫z1z2Tz′2dz.Figure12 shows the global distributions of gravity wave potential energy obtained from FY-3C occultation data at the height of 20–30 km in northern hemisphere winter from 2015 to 2016. The white patches in the map are missing data due to the insufficient number of FY-3C occultation events. The results show that the average of the gravity wave potential energy near the equator is the maximum with a value of 10 J/kg due to strong convection near the equator and decreases toward the two poles, while the potential energy in the northern hemisphere is stronger than that in the southern hemisphere. Also, the gravity wave potential energy over the Eurasian continent is relatively higher. These features are consistent with the excitation mechanism of gravity wave, such as the terrain, wind shear, and so on. Moreover, it is also consistent with other results [26, 27]. In addition, the filtering methods or the density of the grid also affects the accuracy of the gravity wave parameter estimation [29].Figure 12
Global distribution of gravity wave potential energy from FY-3C occultation data at the height of 20–30 km in northern hemisphere winter from 2015 to 2016.Figure13 shows the seasonal distributions of global gravity wave potential energy at the height of 15–20 km. High values of the gravity wave potential energy are basically distributed near the equator and it decreases toward the two poles, which may be related to the convection near the equator.Figure 13
Global distributions of gravity wave potential energy at the height of 15–20 km in northern hemisphere. Spring (a), summer (b), autumn (c), and winter (d).
(a)
(b)
(c)
(d)Furthermore, the gravity wave potential energy in the eastern hemisphere (0°E–180°E) is more obvious than that in the western hemisphere, which may be due to the greater land area of the Eurasia in the eastern hemisphere. At the same time, we can also find high values of gravity wave potential energy at the western side of the Rocky Mountains, the northern side of the Himalayas in Eurasia in northern hemisphere summer, and on the west of the Andes in northern hemisphere winter. To sum up, the gravity wave at the height of 15–20 km is basically excited by the terrain, and there are extreme values of the gravity wave potential energy at the topographic shear.Terrain is one of the main factors that excite gravity wave potential energy [30]. Figure 14 shows the zonal distribution of gravity wave potential energy at the latitude of 0°S–30°S in northern hemisphere summer. There are high values of the gravity wave potential energy below 25 km at the longitude of 120°W, 90°W, 50°W, and 10°E. Examining terrain heights suggests that the high gravity wave potential energy values at longitude 120°W are likely caused by the topography of islands in the ocean, which confirms the conclusion of Alexander et al. [30]. The high values at about 90°W are due to the Andes of South America on the east, the high values at 50°W are due to the continental terrain of the South American continent, and the high values near 10°E are due to the presence of the African continent on the east. Therefore, most gravity waves are excited by the landform and then spread up to the height of about 25 km.Figure 14
Zonal distribution of gravity wave potential energy at the latitude of 0°S–30°S in summer of the northern hemisphere (a) and maximum terrain height (b).
(a)
(b)


## 3.1. Effect of High-Order Ionospheric Delays

In the estimation of the ionospheric correction residual on occultation events, we mainly account for the second-order and the third-order ionospheric delays, and the ionospheric correction residuals. We evaluated all the second-order ionospheric delays of FY-3C satellite neutral atmospheric occultation events in September 2014. In the vicinity of the surface, the second-order delays are within ±15 mm. As the heights increase, the second-order delays also increase. At altitudes above 90 km, the second-order ionospheric delays are about ±20 mm, and the average of the second-order ionospheric delays at all heights is about ±7 mm. The second-order ionospheric delays at 15 km, 45 km, and 90 km are shown in Figure3. At the height of 15 km, almost all of the second-order delays are less than 15 mm; at the height of 45 km, there are less than 5% of the second-order delays greater than 15 mm; and when the height rises to 90 km, the ratio of the second-order delays over 15 mm is about 10%.Figure 3
The statistics histograms of the second-order ionospheric delays at 15 km (a), 45 km (b), and 90 km (c).
(a)
(b)
(c)The residual effect of ionospheric delay after using the dual-frequency correction is called the ionospheric correction residual, and the statistic histograms are shown in Figure4. At the height of 15 km, 98% of the residuals are less than ±10 mm, at the height of 45 km, there are 95% of the residuals less than ±10 mm, but the ratio of the second-order ionospheric correction residuals within ±10 mm is less than 70% at the height of 90 km, which indicates that the effect of the ionospheric correction residual on the upper atmosphere is more pronounced. Since the top-to-bottom method is always used when we invert the refraction rate and the temperature, the ionospheric correction residual of the upper atmosphere will affect the accuracy of the lower atmosphere. Therefore, we usually limit the maximum height of inversion to weaken the effect of the upper atmosphere in the data processing.Figure 4
The statistic histograms of the second-order ionospheric correction residuals at 15 km (a), 45 km (b), and 90 km (c).
(a)
(b)
(c)In addition, the third-order ionospheric delays in the measurement of neutral atmospheric occultation are also estimated, and the values are very small. The maximum of which is 0.5 mm and the average is only about 0.15 mm, which is only 1/40 to 1/30 of the second-order ionospheric delay. Therefore, the second-order ionospheric delay is the main source of the ionospheric correction residual.We further estimated the effect of the ionospheric correction residual on the bending angle and temperature of the neutral atmosphere. Here, the second-order ionospheric delay is directly eliminated from the original observation, and then the inversion is performed according to the ionospheric-free linear combination of both GPS frequencies to obtain the bending angle, the refractive index, and the temperature of the atmosphere, which are then compared with the parameters obtained from the original data without eliminating the high-order ionospheric delay.In order to verify the feasibility of correcting the high-order ionospheric term by directly subtracting the second-order ionospheric delay from the original observation data, we compared the results from colocated occultation and radiosonde measurements for 9 different radiosonde stations (Figure5). It was found that most temperature profiles obtained through the correction of high-order ionospheric delays are closer to the results of radiosonde, and the temperature measured by radiosonde is not affected by the ionosphere. Therefore, it is feasible to use this method to effectively eliminate the effect of the high-order ionospheric delays.Figure 5
Effect of the ionospheric correction residual on the neutral atmospheric temperature estimation at 9 radiosonde stations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Furthermore, 7,144 occultation events were estimated with and without eliminating the high-order ionospheric delays, which were divided into three groups according to the maximum residuals of temperature profiles between 2 K, 4 K, and 6 K. There are 6,508 events located within 2 K, accounting for 91.1% of the total sample, 375 residuals of the temperature are between 2 K and 4 K, accounting for 5.2%, and the residuals of the temperature between 4 K and 6 K are 261 times, accounting for 3.7% of the estimated total sample. And the average temperature residual of ionospheric correction at all altitudes is about 1 K.


## 3.2. Comparison with COSMIC and Radiosonde

### 3.2.1. Comparison between COSMIC and FY-3C Data

In order to verify the validity of FY-3C occultation data, we obtained COSMIC occultation data and looked for the overlapping occultation events from COSMIC and FY-3C within 3° in latitude, 3° in longitude, and 3 hours. There are totally 255 overlapping events from September 21th to 28th, 2014. The profiles of four overlapping occultation events are shown in Figure6. The inversion results of the temperature profiles from FY-3C occultation data are similar to those from COSMIC occultation data, but the temperature profiles of FY-3C are slightly smoother than those of COSMIC.Figure 6
Comparison of the temperature profiles obtained from FY-3C and COSMIC occultation data.
(a)
(b)
(c)
(d)From the statistics of the difference and the standard deviation of the temperature profiles at the height of 10 km, 15 km, 20 km, 25 km, 30 km, and 35 km (Figure7), we find that the average differences between the FY-3C and COSMIC temperature profiles are within 1°C, and the standard deviation is about 2 to 3°C. In the refractive index curve, we can also see the consistency of the refractive index obtained by FY-3C and COSMIC occultation data.Figure 7
Statistics of temperature profiles and the relative refraction index obtained from FY-3C and COSMIC occultation data.
(a)
(b)In addition, we obtained and compared the tropopause height and tropopause temperature from COSMIC and FY-3C occultation data, respectively (Figure8). It was found that there is a strong correlation of the tropopause height between FY-3C data and COSMIC data, and the correlation coefficient is 0.94, indicating that they have similar capability to retrieve the tropopause height. The tropopause temperatures from FY-3C data and COSMIC data also have a very strong correlation, with a correlation coefficient of 0.96.Figure 8
Comparison of the tropopause height (a) and temperature (b) obtained from FY-3C and COSMIC occultation data.
(a)
(b)


### 3.2.2. Comparison between Radiosonde and FY-3C Data

In order to verify the accuracy of FY-3C occultation data, the independent radiosonde data are further used. We chose radiosonde data from the IGRA program for the experiment, and the criteria for co-located station observations is also set as 3° in latitude, 3° in longitude, and 3 hours. Figure9 shows the temperature profiles of four repeated observations. We find that the temperature profiles are very similar in the region below 30 km, which indicates that FY-3C occultation data have a strong reliability.Figure 9
Comparison of the temperature profiles obtained from FY-3C GPS radio occultation data and radiosonde data.
(a)
(b)
(c)
(d)From the statistics of the difference and the standard deviation of the temperature profiles at the heights of 10 km, 15 km, 20 km, 25 km, and 30 km (Figure10), we find that the difference of the temperature between FY-3C and radiosonde is within 2°C, and the standard deviation is about 3°C. The difference of the refractive index calculated by radiosonde and FY-3C is very small with a relative error of less than 0.5%, and the standard deviation is less than 1.4%, which indicates that FY-3C occultation data have the same accuracy as radiosonde data in terms of statistics.Figure 10
Statistics of the temperature profiles and the relative refraction index obtained from FY-3C occultation data and radiosonde data.
(a)
(b)In addition, we further compared the tropopause height and tropopause temperature obtained from radiosonde and FY-3C (Figure11). The correlation coefficient of the tropopause height is 0.93, and the correlation coefficient of the tropopause temperature is 0.95. It indicates that the tropopause height and temperature estimated from FY-3C are equivalent to that from radiosonde, which is consistent with the conclusion of Liao et al. [25].Figure 11
Comparison of the tropopause height (a) and temperature (b) obtained from FY-3C occultation data and radiosonde data.
(a)
(b)


## 3.2.1. Comparison between COSMIC and FY-3C Data

In order to verify the validity of FY-3C occultation data, we obtained COSMIC occultation data and looked for the overlapping occultation events from COSMIC and FY-3C within 3° in latitude, 3° in longitude, and 3 hours. There are totally 255 overlapping events from September 21th to 28th, 2014. The profiles of four overlapping occultation events are shown in Figure6. The inversion results of the temperature profiles from FY-3C occultation data are similar to those from COSMIC occultation data, but the temperature profiles of FY-3C are slightly smoother than those of COSMIC.Figure 6
Comparison of the temperature profiles obtained from FY-3C and COSMIC occultation data.
(a)
(b)
(c)
(d)From the statistics of the difference and the standard deviation of the temperature profiles at the height of 10 km, 15 km, 20 km, 25 km, 30 km, and 35 km (Figure7), we find that the average differences between the FY-3C and COSMIC temperature profiles are within 1°C, and the standard deviation is about 2 to 3°C. In the refractive index curve, we can also see the consistency of the refractive index obtained by FY-3C and COSMIC occultation data.Figure 7
Statistics of temperature profiles and the relative refraction index obtained from FY-3C and COSMIC occultation data.
(a)
(b)In addition, we obtained and compared the tropopause height and tropopause temperature from COSMIC and FY-3C occultation data, respectively (Figure8). It was found that there is a strong correlation of the tropopause height between FY-3C data and COSMIC data, and the correlation coefficient is 0.94, indicating that they have similar capability to retrieve the tropopause height. The tropopause temperatures from FY-3C data and COSMIC data also have a very strong correlation, with a correlation coefficient of 0.96.Figure 8
Comparison of the tropopause height (a) and temperature (b) obtained from FY-3C and COSMIC occultation data.
(a)
(b)


## 3.2.2. Comparison between Radiosonde and FY-3C Data

In order to verify the accuracy of FY-3C occultation data, the independent radiosonde data are further used. We chose radiosonde data from the IGRA program for the experiment, and the criteria for co-located station observations is also set as 3° in latitude, 3° in longitude, and 3 hours. Figure9 shows the temperature profiles of four repeated observations. We find that the temperature profiles are very similar in the region below 30 km, which indicates that FY-3C occultation data have a strong reliability.Figure 9
Comparison of the temperature profiles obtained from FY-3C GPS radio occultation data and radiosonde data.
(a)
(b)
(c)
(d)From the statistics of the difference and the standard deviation of the temperature profiles at the heights of 10 km, 15 km, 20 km, 25 km, and 30 km (Figure10), we find that the difference of the temperature between FY-3C and radiosonde is within 2°C, and the standard deviation is about 3°C. The difference of the refractive index calculated by radiosonde and FY-3C is very small with a relative error of less than 0.5%, and the standard deviation is less than 1.4%, which indicates that FY-3C occultation data have the same accuracy as radiosonde data in terms of statistics.Figure 10
Statistics of the temperature profiles and the relative refraction index obtained from FY-3C occultation data and radiosonde data.
(a)
(b)In addition, we further compared the tropopause height and tropopause temperature obtained from radiosonde and FY-3C (Figure11). The correlation coefficient of the tropopause height is 0.93, and the correlation coefficient of the tropopause temperature is 0.95. It indicates that the tropopause height and temperature estimated from FY-3C are equivalent to that from radiosonde, which is consistent with the conclusion of Liao et al. [25].Figure 11
Comparison of the tropopause height (a) and temperature (b) obtained from FY-3C occultation data and radiosonde data.
(a)
(b)


## 3.3. Gravity Wave Potential Energy from FY-3C GPS Radio Occultation

Considering the volume and spatial and temporal distribution of FY-3C occultation data, the Earth is meshed into a 90 × 90 × 1 grid with a spatial resolution of 2° in latitude, 4° in longitude, and a temporal resolution of 1 month. The temperature profile is interpolated with a resolution of 200 m, and the average is taken as the value of the grid node.First, we performed S-transform on each height of every latitude circle and extracted the components whose wavenumber are of 0–6 to eliminate the interferences of some zonal disturbances, like the planetary wave [26, 27]. The method to separate disturbance temperature field from the background temperature field is to conduct continuous wavelet transform and then use the Butterworth filter with the interception frequency of 0.1–0.5 × 103 MHz for bandpass filtering, and finally, the inverse transformation is performed to get the background temperature field.The atmospheric gravity wave is produced by tiny air mass through reciprocating vibration under the influence of gravity and buoyancy. We call this oscillation frequency buoyant frequency (Brumt–Vaisala frequency):(12)Nz2=gT¯∂T¯∂z+gcp,where Nz is the buoyant frequency, g is the gravitational acceleration, T¯ is the background temperature field, and cp is the heat capacity at constant pressure. The gravity wave potential energy can be described as(13)Epz=12gNz2T′T¯2,where Ep is the gravity wave potential energy and T′ is the gravity wave disturbance temperature. As for the processing mode of T′2, we adopted the method of Tsuda et al. and take 2 km as the window and 200 m as the step size to calculate the smoothing values in the window [28]:(14)T′2=1z2−z1∫z1z2Tz′2dz.Figure12 shows the global distributions of gravity wave potential energy obtained from FY-3C occultation data at the height of 20–30 km in northern hemisphere winter from 2015 to 2016. The white patches in the map are missing data due to the insufficient number of FY-3C occultation events. The results show that the average of the gravity wave potential energy near the equator is the maximum with a value of 10 J/kg due to strong convection near the equator and decreases toward the two poles, while the potential energy in the northern hemisphere is stronger than that in the southern hemisphere. Also, the gravity wave potential energy over the Eurasian continent is relatively higher. These features are consistent with the excitation mechanism of gravity wave, such as the terrain, wind shear, and so on. Moreover, it is also consistent with other results [26, 27]. In addition, the filtering methods or the density of the grid also affects the accuracy of the gravity wave parameter estimation [29].Figure 12
Global distribution of gravity wave potential energy from FY-3C occultation data at the height of 20–30 km in northern hemisphere winter from 2015 to 2016.Figure13 shows the seasonal distributions of global gravity wave potential energy at the height of 15–20 km. High values of the gravity wave potential energy are basically distributed near the equator and it decreases toward the two poles, which may be related to the convection near the equator.Figure 13
Global distributions of gravity wave potential energy at the height of 15–20 km in northern hemisphere. Spring (a), summer (b), autumn (c), and winter (d).
(a)
(b)
(c)
(d)Furthermore, the gravity wave potential energy in the eastern hemisphere (0°E–180°E) is more obvious than that in the western hemisphere, which may be due to the greater land area of the Eurasia in the eastern hemisphere. At the same time, we can also find high values of gravity wave potential energy at the western side of the Rocky Mountains, the northern side of the Himalayas in Eurasia in northern hemisphere summer, and on the west of the Andes in northern hemisphere winter. To sum up, the gravity wave at the height of 15–20 km is basically excited by the terrain, and there are extreme values of the gravity wave potential energy at the topographic shear.Terrain is one of the main factors that excite gravity wave potential energy [30]. Figure 14 shows the zonal distribution of gravity wave potential energy at the latitude of 0°S–30°S in northern hemisphere summer. There are high values of the gravity wave potential energy below 25 km at the longitude of 120°W, 90°W, 50°W, and 10°E. Examining terrain heights suggests that the high gravity wave potential energy values at longitude 120°W are likely caused by the topography of islands in the ocean, which confirms the conclusion of Alexander et al. [30]. The high values at about 90°W are due to the Andes of South America on the east, the high values at 50°W are due to the continental terrain of the South American continent, and the high values near 10°E are due to the presence of the African continent on the east. Therefore, most gravity waves are excited by the landform and then spread up to the height of about 25 km.Figure 14
Zonal distribution of gravity wave potential energy at the latitude of 0°S–30°S in summer of the northern hemisphere (a) and maximum terrain height (b).
(a)
(b)


## 4. Summary

In this paper, we analyzed the influence of high-order ionospheric delay on atmospheric occultation and proposed a corresponding correction algorithm to eliminate the effect. The effect of the ionospheric correction residual on the phase delay is less than 20 mm, which has little effect on bending angle and about 1 K on average temperature profile. By comparing with COSMIC occultation data and radiosonde data, the difference of the temperature between FY-3C occultation data and COSMIC occultation data is within 1°C, and the difference between FY-3C occultation data and the radiosonde data is within 2°C. This proves that the data from Chinese meteorological satellite FY-3C are highly reliable and can be used for in-depth study of atmospheric science. Finally, the global gravity wave potential energy is obtained from FY-3C occultation data, which are consistent with those obtained by COSMIC occultation data. The high values of gravity wave potential energy are basically distributed near the equator and decrease toward the two poles. Most of the gravity waves below 25 km are mainly excited by the terrain. In the future, with further development and improvement of more BeiDou constellations and GNOS RO receiver, novel BeiDou RO observations and applications are expected in atmospheric sounding [31].



---

*Source: 4780143-2019-08-04.xml*