Mountain Top-Based Atmospheric Radio Occultation Observations with Open / Closed Loop Tracking: Experiment and Validation

: Through Global Navigation Satellite System (GNSS) occultation measurement, the global ionosphere and atmosphere can be observed. When the navigation satellites’ signal passes through the lower atmosphere, the rapid change of the atmospheric refractive index gradient will cause serious multipath phenomena in radio wave propagation. Atmospheric doppler frequency shift and amplitude signal ﬂuctuations increase drastically. Due to the attenuation of signal amplitude and the rapid change of the Doppler frequency, the general phase locked loop (PLL) cannot work properly. Hence, a more stable tracking technology is needed to track the occultation signal passing through the lower atmosphere. In this paper, a mountain-top based radio occultation experiment is performed, where we employ an open-loop receiver and remove the navigation bits by the internal demodulation.

affected by the change of the signal [27,28]. It can track full information of phases (the effects of motion, atmospheric multipath, etc.), and can detect rising occultation events more conveniently [29,30].
The mountain-based OL tracking algorithm is different from the spaceborne OL tracking algorithm in terms of the model Doppler and pseudo range. In the satellite OL tracking module, we employ the international reference atmospheric model CIRA86-Q for pseudo-range/Doppler prediction. The mountain-based OL tracking algorithm is mainly used to verify the correctness of OL tracking, not to verify the correctness of the Doppler and pseudo range model. In the verification process of the mountain-based OL tracking algorithm, we assume that the Doppler and pseudo-range are mainly changed by the movement of GNSS satellites, while the effects of the atmosphere and the ionosphere are small. In the mountain-based test, the position and speed of the navigation satellite can be obtained through the ephemeris. The geometric distance of between the station and the satellite at an occultation and the Doppler of motion obtained based on the velocity of that satellite are used as the output of the pseudo-range/Doppler model to narrow the phase search range in full OL tracking. The chip array of the mountain-based occultation receiver is expanded to ensure that the phase delay caused by the atmosphere is within the chips' range.

Open Open-Loop Data Processing
In the OL mode, the phase used to retrieve for atmosphere profiles can be expressed as sum of the modeled phase (MP) and the recovered phase (RP). The MP is obtained from the occultation geometry and climate models. The RP is restored from observations of the In-Phase(I)/Quadra-Phase(Q) components through three steps: frequency reduction conversion; NDM removal; and adjacent phase connection [28,[31][32][33][34]. The ways to eliminate the influence of NDM in the OL mode are the internal and external methods, and the biggest difference between them is that the external mode requires GPS bit data. Due to geographical restrictions, however, only about 70% of the COSMIC occultation data has corresponding GPS bit data [31][32][33][34]. Due to a lack of information in determining the initial sign of RP, there may be a different sign between the initial I/Q components and the real phase, which will negatively correct the MP. We adopt the strategy of the full OL data recovery process, assisted with closed-loop observations to avoid the above-mentioned problems. Through MPs adding or subtracting RPs, we can obtain two series phases, named corrected phases A/B(CPA/CPB). We use the correlation between the mountain-based closed-loop observation data of an elevation angle above 3 degrees and CPA/CPB to determine the initial sign of RP, and obtain a high-consistency corrected OL phase for subsequent atmospheric parameter inversion.

Mountain-Based Atmospheric Retrieval Algorithm
The GNSS signal can generally be expressed as follows: L s,r i (t) = ρ s,r (t) + c(σ s (t) + σ r (t)) + exL i (t) + m s,r (t) + n s,r (t) + N s,r (t) (1) exL s,r i (t) = Ion s,r i (t) + Atm s,r i (t) (2) where c is the speed of light, and L s,r i (t) is the observed phase of the transmitter (s) and receiver (r) of the i-th frequency at time (t). ρ s,r (t) is the geometric distance between transmitter (s) and receiver (r) at time t. σ s (t), σ r (t) are the clock difference of the transmitter (s) and receiver (r) at time (t), respectively. m s,r i (t) is the local multipath error [35,36]; n s,r i (t) is the residual phase. N s,r i (t) is the ambiguity of the whole cycle. The excess phase exL s,r i (t) of transmitter s and receiver (r) of the i-th frequency at time (t) is generally defined as the sum of neutral atmospheric delay Atm s,r i (t) and ionospheric delay Ion s,r i (t) [36][37][38][39][40]. Based on Internet GNSS Servers (IGS) products and low-orbit satellite precision orbit determination products, it is possible to obtain a sufficiently accurate distance ρ s,r (t) from the navigation satellite transmitter to the low-orbit satellite receiver and its clock difference σ s (t). The local multipath error can be ignored in mountain-based occultation. The receiver clock error σ r (t) can generally be eliminated by selecting a higher altitude navigation star as a reference and applying the satellites' phase difference. As a result, high-precision excess phase observations of the occultation can be obtained, which only include the influence of the ionosphere and atmosphere through the path. In the process of spaceborne neutral atmosphere inversion, the combination of different frequency point bending angle sequences can better eliminate the ionospheric influence. In extracting the excess phase, only the ionospheric influence of the reference star can be eliminated. In the mountain-based atmospheric occultation inversion, the atmospheric state information from the ground to the top of the mountain can be obtained through the relationship between the impact height observed by the positive and negative altitude angle. As shown in Figure 1, observation of negative elevation angle (light and dark green thick line, bending angle is α − ) can be regarded as the Abel integral (light green thick line) and positive elevation angle (light blue thick line, bending angle is α + ) from the top to the impact height, namely [39][40][41][42]: Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 17 can generally be eliminated by selecting a higher altitude navigation star as a reference and applying the satellites' phase difference. As a result, high-precision excess phase observations of the occultation can be obtained, which only include the influence of the ionosphere and atmosphere through the path. In the process of spaceborne neutral atmosphere inversion, the combination of different frequency point bending angle sequences can better eliminate the ionospheric influence. In extracting the excess phase, only the ionospheric influence of the reference star can be eliminated. In the mountain-based atmospheric occultation inversion, the atmospheric state information from the ground to the top of the mountain can be obtained through the relationship between the impact height observed by the positive and negative altitude angle. As shown in Figure 1, observation of negative elevation angle (light and dark green thick line, bending angle is - ) can be regarded as the Abel integral (light green thick line) and positive elevation angle (light blue thick line, bending angle is   ) from the top to the impact height, namely [39][40][41][42]   Let α(a) = α − (a) − α + (a) = −2a Remote Sens. 2020, 12, 4078 5 of 16 then the refractive index n of the impact height α (n 1 is the peak refractive index, ξ is the maximum impact height of the occultation) is: Thus, the refractive index profile from the ground to the top of the mountain can be obtained from the mountain-based occultation data [42,43].
In the above process, the influence of the atmosphere or ionosphere on the observation data above the top of the mountain can be eliminated in the process of searching for the symmetrical impact height and its combination. Therefore, when the clock error of the mountaintop receiver is eliminated by the reference star, only the ionospheric influence of the reference star needs to be eliminated, while the mountain-based occultation observation data does not need to consider its influence on the excess phase. Even with only single frequency observation data, the refractive index inversion results from the ground to the top of the mountain can be obtained.
As shown in Figure 2, according to the geometric relationship of the mountain-based occultation, Excess phase exL s,r i (t) and Doppler dop s,r i (t) of transmitter (s) and receiver (r) at frequency (i) at time (t), The impact height a (the superscript and subscripts and the time label are omitted here) and the bending angle α have the following relationship [44]: a(t) = n 2 r 2 sin(φ 2 ) = n 1 r 1 sin(φ 1 ) → φ 1 = arcsin( n 2 r 2 sin(φ 2 ) n 1 r 1 ) then the refractive index n of the impact height  ( 1 n is the peak refractive index, is the maximum impact height of the occultation) is: Thus, the refractive index profile from the ground to the top of the mountain can be obtained from the mountain-based occultation data [42,43].
In the above process, the influence of the atmosphere or ionosphere on the observation data above the top of the mountain can be eliminated in the process of searching for the symmetrical impact height and its combination. Therefore, when the clock error of the mountaintop receiver is eliminated by the reference star, only the ionospheric influence of the reference star needs to be eliminated, while the mountain-based occultation observation data does not need to consider its influence on the excess phase. Even with only single frequency observation data, the refractive index inversion results from the ground to the top of the mountain can be obtained.
As shown in Figure   Among these, f i is the frequency of frequency point i, and φ 1 , φ 2 , r 1 , r 2 are the signal incident angle and the center distance of the receiver and transmitter, respectively. θ is the angle between the vector from the Earth's center to the receiver and the transmitter. v 2 is the occultation velocity, and γ is the angle between the occultation velocity and the vector from the Earth's center to the transmitter.

Mountain-Based Test Results
The mountain-based occultation experiment is from 25 August to 8 September 2020. The site of the observation is the top of Wuling Mountain, the main peak of the Yanshan Mountains (located in Chengde, Hebei Province, 117.48 degrees east longitude, 40.60 degrees north latitude, and 2118 m above sea level), and the occultation antenna is deployed on the roof of Dingfeng Hotel. As shown in Figures 3 and 4, the south side of the mountain is empty. Our occultation antenna (right side) is installed on the top of the building and placed facing south, and the receiver is installed inside the room.
the mountain base test, at the top of Wuling Mountain. Coordinate point 2 represents the navigation satellite at an occultation event.
Among these, i f is the frequency of frequency point i , and 1 , 2 , 1 , 2 are the signal incident angle and the center distance of the receiver and transmitter, respectively.  is the angle between the vector from the Earth's center to the receiver and the transmitter.
2 v is the occultation velocity, and  is the angle between the occultation velocity and the vector from the Earth's center to the transmitter.

Mountain-Based Test Results
The mountain-based occultation experiment is from 25 August to 8 September 2020. The site of the observation is the top of Wuling Mountain, the main peak of the Yanshan Mountains (located in Chengde, Hebei Province, 117.48 degrees east longitude, 40.60 degrees north latitude, and 2118 m above sea level), and the occultation antenna is deployed on the roof of Dingfeng Hotel. As shown in Figures 3 and 4, the south side of the mountain is empty. Our occultation antenna (right side) is installed on the top of the building and placed facing south, and the receiver is installed inside the room. Mountain-based location, behind our experiment team is the bottom of the mountain. Our occultation antenna (the part marked in red in the picture) is mounted on a tripod and placed facing south, and our positioning antenna (the part marked in green in the picture) is also installed on the tripod. On the tripod, place it towards the zenith. The south side of the mountain is empty, which is conducive to our experiments.
The GNSS occultation receiver system used for testing is developed by Tianjin Yunyao Aerospace Technology Co., Ltd, which is compatible with GPS/BD systems. Our occultation antenna (the part marked in red in the picture) is mounted on a tripod and placed facing south, and our positioning antenna (the part marked in green in the picture) is also installed on the tripod. On the tripod, place it towards the zenith. The south side of the mountain is empty, which is conducive to our experiments.

Analysis of GPS and BD Open-Loop Tracking Results
The I and Q components of the signal have different characteristics in OL and closed-loop tracking. In the tracking process, the phases of the OL flipped by NDM and error from the OL predict model in the form of a trigonometric function on the I/Q components, which can be seen from Figure 5. The phases of the closed-loop are of huge value, with signal flipped in the I component, and the Q component is signal noise. The closed-loop tracking can detect the phase or frequency difference between the locally reproduced carrier and the received carrier through the phase detector or frequency discriminator in the carrier loop. Then, the phase alignment of the local PN code and the received code can be continuously adjusted through the code loop discriminator, so as to obtain the highest self-phase correlation peak value. It has the characteristics of small volume and low power consumption. It is also the equipment we will use in satellite networking in the future. It is placed on the anti-static pad, and ground treatment is also carried out to prevent the occurrence of static electricity and lightning strike events on the mountain.

Analysis of GPS and BD Open-Loop Tracking Results
The I and Q components of the signal have different characteristics in OL and closed-loop tracking. In the tracking process, the phases of the OL flipped by NDM and error from the OL predict model in the form of a trigonometric function on the I/Q components, which can be seen from Figure  5. The phases of the closed-loop are of huge value, with signal flipped in the I component, and the Q component is signal noise. The closed-loop tracking can detect the phase or frequency difference between the locally reproduced carrier and the received carrier through the phase detector or frequency discriminator in the carrier loop. Then, the phase alignment of the local PN code and the received code can be continuously adjusted through the code loop discriminator, so as to obtain the highest self-phase correlation peak value.  It is also the equipment we will use in satellite networking in the future. It is placed on the anti-static pad, and ground treatment is also carried out to prevent the occurrence of static electricity and lightning strike events on the mountain.
The GNSS occultation receiver system used for testing is developed by Tianjin Yunyao Aerospace Technology Co., Ltd, which is compatible with GPS/BD systems.  figure (red marking part) and verified by the company's self-researched occultation detector. It has the characteristics of small volume and low power consumption. It is also the equipment we will use in satellite networking in the future. It is placed on the anti-static pad, and ground treatment is also carried out to prevent the occurrence of static electricity and lightning strike events on the mountain.

Analysis of GPS and BD Open-Loop Tracking Results
The I and Q components of the signal have different characteristics in OL and closed-loop tracking. In the tracking process, the phases of the OL flipped by NDM and error from the OL predict model in the form of a trigonometric function on the I/Q components, which can be seen from Figure  5. The phases of the closed-loop are of huge value, with signal flipped in the I component, and the Q component is signal noise. The closed-loop tracking can detect the phase or frequency difference between the locally reproduced carrier and the received carrier through the phase detector or frequency discriminator in the carrier loop. Then, the phase alignment of the local PN code and the received code can be continuously adjusted through the code loop discriminator, so as to obtain the highest self-phase correlation peak value. This paper analyzes the mountain-based occultation observations from full OL tracking mode in the following ways. First, using the excess phase obtained from the observation data in the closedloop tracking mode as the reference, the continuity and duration of the excess phases obtained in the This paper analyzes the mountain-based occultation observations from full OL tracking mode in the following ways. First, using the excess phase obtained from the observation data in the closed-loop tracking mode as the reference, the continuity and duration of the excess phases obtained in the OL tracking verifies that the full OL tracking mode of GPS and BD are more robust than the closed-loop tracking mode. Second, the refractivity profiles retrieved from observations in the OL and closed-loop tracking modes are compared to verify their consistency.
We selected (part of) the data on September 7 for analysis. The data include rising and falling occultations. As shown in Figures 6 and 7, excess phases retrieved from observations of L1C/GPS in OL/closed-loop tracking mode are compared, as well as that of B1/BD. Figure 8 shows the difference of the excess phase retrieved from observations of L1C/GPS and B1/BD in OL/closed-loop tracking mode, which are within 0.2 m. From Figures 6-8, the following can be concluded: (1) OL tracking time is longer than closed-loop, that is, for low-angle occultation events, OL has a better tracking effect. (2) The mountain-based occultation data can only detect the peak height and below. It can be seen from Figures 6 and 7 that compared to closed-loop tracking, the OL tracking can capture a lower height, and has the same result as the closed-loop tracking. Hence, the continuity of the OL is relatively robust. From the curve, we can see that the excess phase changes obtained by GPS and BD are relatively robust, and there is no deterioration in signal capture quality due to drastic changes in the signal.       Figure 11 visualizes the difference of refractivity profiles retrieved from observations of B1/BD and L1C/GPS in OL/closed-loop mode, from which it is concluded that: 1. the bias of the refractivity retrieved from observations of L1C/GPS and B1/BD in OL/closed-loop mode is within 1N; and 2. The trend of the refractivity profiles retrieved from observations of L1C/GPS in OL/closed-loop mode is slightly better than the result of B1/BD. refractivity profiles from B1/BD signal in OL mode are about 1.2 km, and that of closed-loop mode are between 1.3-1.5 km. Figure 11 visualizes the difference of refractivity profiles retrieved from observations of B1/BD and L1C/GPS in OL/closed-loop mode, from which it is concluded that: 1. the bias of the refractivity retrieved from observations of L1C/GPS and B1/BD in OL/closed-loop mode is within 1N; and 2. The trend of the refractivity profiles retrieved from observations of L1C/GPS in OL/closed-loop mode is slightly better than the result of B1/BD.

Comparison of Closed-Loop Tracking and Open-Loop Tracking
In order to verify the effectiveness of the OL tracking algorithm, the comparison of the OL and closed-loop methods is used for verification, and the OL and closed-loop channels are used to capture and track the same satellite. Compare the start capture angle and end capture angle of the two to verify the capture effect. We selected 7 days of observation data for statistics, and the selected data are shown in Table 1:

Comparison of Closed-Loop Tracking and Open-Loop Tracking
In order to verify the effectiveness of the OL tracking algorithm, the comparison of the OL and closed-loop methods is used for verification, and the OL and closed-loop channels are used to capture and track the same satellite. Compare the start capture angle and end capture angle of the two to verify the capture effect. We selected 7 days of observation data for statistics, and the selected data are shown in Table 1: From Table 1, it can be seen that the full OL tracking method can effectively capture occultation events at lower angles. As shown in Table 1, observations from an occultation of GPS (PRN 30) at 1 September 2020 0:10:43 UT showed that its OL end pitch angle reached −1.87 deg, which is 0.88 deg smaller than its closed-loop end pitch angle. At low elevation angles, it is difficult to track in the closed-loop strategy stably, and problems such as drastic changes in the signal-to-noise ratio and loss of signal lock will occur. The fully OL tracking method effectively solves the above problems.
We have made statistics on all the data to evaluate the effect of OL tracking. One calculation method is to calculate the angles of the whole occultation tracking, and mark it as strategy A; the other is to calculate the negative angles of the occultation tracking, and mark it as strategy B. For strategy A or B, the way to evaluate the capture capability of the OL/closed loop tracking, and the effect of OL tracking are: in which α i , α e are the initial pitch angle and end pitch angle, respectively; E A o,c , E B o,c are values to evaluate the capture capability of the OL/closed-loop tracking with strategy A or B; and A, B are the effects of OL tracking evaluated with strategy A or B. According to observations from Table 1, the result of accumulating the seven-day data are: The results show that in the case of counting all angles, the OL improves the capture angle range of data volume by nearly 20%, and when only the negative angle is calculated, the OL improves the capture angle value by nearly 89%. The mountain-based situation is different from the onboard situation. We will continue to verify the actual capture effect in the future. However, there is no doubt that the fully OL capture method greatly improves the capture effect. We make the difference between the OL angle and the closed-loop angle (only the negative angle is counted), and the result is shown in Figure 12. There is no specific coupling relationship between OL tracking and closed-loop tracking, but the overall effect is better than closed-loop. The reason why there is no specific coupling relationship is due to the difference between the two capture methods. The angle difference is mainly concentrated between 1 • and 1.
in which , are the initial pitch angle and end pitch angle, respectively; E , , E , are values to evaluate the capture capability of the OL/closed-loop tracking with strategy A or B; and , are the effects of OL tracking evaluated with strategy A or B. According to observations from Table 1, the result of accumulating the seven-day data are: The results show that in the case of counting all angles, the OL improves the capture angle range of data volume by nearly 20%, and when only the negative angle is calculated, the OL improves the capture angle value by nearly 89%. The mountain-based situation is different from the onboard situation. We will continue to verify the actual capture effect in the future. However, there is no doubt that the fully OL capture method greatly improves the capture effect. We make the difference between the OL angle and the closed-loop angle (only the negative angle is counted), and the result is shown in Figure 12. There is no specific coupling relationship between OL tracking and closedloop tracking, but the overall effect is better than closed-loop. The reason why there is no specific coupling relationship is due to the difference between the two capture methods. The angle difference is mainly concentrated between 1° and 1.2°. Figure 12. Histogram of the difference of the start/end elevation angle from a rising/setting occultation event between closed-loop and OL. Because the top height of the mountain is about 2 km, the minimum elevation angle value is generally about −3° in OL mode, and the difference is mainly concentrated between 0.8° and 1.2°.

Comparison of Open-Loop Tracking Inversion Results and Meteorological Observations
The atmosphere profiles retrieved from observations in OL tracking mode are evaluated with the local observations of temperature, humidity and pressure provided by the Beijing Meteorological Bureau.
The atmospheric refractive index of the observation point can be calculated by temperature, relative humidity and pressure data, which is used to verify the accuracy of the inversion of the mountain-based experimental data. The calculation method is as follows [45][46][47][48]: Using the temperature data of the analysis field to calculate the corresponding saturated vapor pressure: Figure 12. Histogram of the difference of the start/end elevation angle from a rising/setting occultation event between closed-loop and OL. Because the top height of the mountain is about 2 km, the minimum elevation angle value is generally about −3 • in OL mode, and the difference is mainly concentrated between 0.8 • and 1.2 • .

Comparison of Open-Loop Tracking Inversion Results and Meteorological Observations
The atmosphere profiles retrieved from observations in OL tracking mode are evaluated with the local observations of temperature, humidity and pressure provided by the Beijing Meteorological Bureau. The atmospheric refractive index of the observation point can be calculated by temperature, relative humidity and pressure data, which is used to verify the accuracy of the inversion of the mountain-based experimental data. The calculation method is as follows [45][46][47][48]: Using the temperature data of the analysis field to calculate the corresponding saturated vapor pressure: E(t) = E 0 × 10 7.5t 237.3+t (14) In the above formula, E(t) is the saturated vapor pressure at the corresponding temperature, and t is the temperature in Celsius. The constant E 0 = 6.1078hpa. There is the following relationship between actual vapor pressure P w , saturated vapor pressure and relative humidity: Finally, the approximate relationship between refractive index and atmospheric parameters can be used to obtain the atmospheric refractive index: In the above formula, N is the refractive index of the atmosphere, p is the total atmospheric pressure (hpa), P w is the water vapor pressure (hpa) and T is the atmospheric temperature (K).
For each sounding profile data, find the matching occultation profile (one or more), and interpolate them to the altitude points with a distance of 0.05 km. Calculate the standard deviation of the average refractive index deviation and the absolute refractive index deviation within the height range covered by both the occultation and the sounding profile as the evaluation standard.
The results are shown in the following Table 2 (partially): It can be seen from the results that the BIAS/MSER of the OL tracking method is within 0.5% and 4%, respectively, which meet the accuracy of its applications (<0.5%, in BIAS; <5%, in MSER).

Conclusions
1. We designed and developed the occultation equipment that supports the GPS/BD systems, and tested its performance based on mountain-based experiments. It is concluded that B1/BD observation data can obtain atmospheric retrieval accuracy similar to L1C/GPS observation. Since the 2nd and 3rd generations of the BD constellations have more than 50 satellites in orbit, the equipment will be deployed on low-orbit satellites in the future to obtain more occultation observations.