Revision of the atmospheric modeling for SNR observations in ground-based GNSS Reflectometry

: The application of signal-to-noise (SNR) observations from ground-based GNSS Reflectometry is becoming an operational tool for coastal sea-level altimetry. As in all data analyses, systematic influences must be reduced here too, to achieve reliable results. A prominent influence results from atmospheric refraction. Different approaches exist to describe or to correct for this influence. In our contribution we will revise the latest developments and suggest a simple atmospheric interferometric delay model that takes into account ray bending as well as along-path propagation delay. The findings are double-checked by numerical experiments based on a step-by-step raytracing procedure. of the reflected signal passes through the specular point and the end of the mirrored antenna below the water surface. The bending angle  e sp at the specular point is spanned by the bent ray path of the reflected signal and the connection from the specular point to the satellite. The bent ray interferometric delay  bent is derived similar to the vacuum case as the difference between the lengths of bent ray path of the reflected signal and that of the rotated bent ray path of the


Introduction
In 1993, Martin-Neira [1] proposed first to use GNSS reflectometry (GNSS-R) for the observation of sea surface heights. Since then, the basic concepts have been adopted for ground-based applications by many groups, which developed various strategies to analyze in particular the oscillating structure of the signal-to-noise ratio (SNR). The wide manifold of approaches to estimate the height of a GNSS antenna above a water surface reach from the frequency analysis by means of Lomb-Scargle Periodogram [2] over Wavelet Analysis Periodogram [3] to inverse modeling [4] of full models even for moving platforms [5] or in real-time [6].
The latest results from an inter-comparison campaign [7] demonstrated an excellent agreement and the capability to derive sea surface heights with a quality of better than 5 cm. Hence, groundbased GNSS-R by means of SNR-analysis seems to be developing into an operational tool for coastal sea-level altimetry, that could possibly reach the quality level of well-established observation methods used by conventional tide gauges. This goal could only be achieved if the full potential of the method can be utilized. For this, it is necessary to take into account as many systematic effects as possible.
The analysis of SNR-data bases on the fact that the direct and the reflected signal from a GNSS satellite interfere at the antenna. The relative phase between the direct and the reflected signal yield an oscillation of the SNR for a moving satellite, that is a function of the interferometric delay between both signals and the wave length of the signal [8]. In most approaches a multipath delay model according to [9] is used (Figure 1).

Figure 1.
Simple sketch of the interferometric delay geometry. The signal from the satellite is reflected at specular point (sp) and received at the antenna together with the direct signal. Dashed lines below the water surface represent the mirrored antenna height H and the reflected signal for an elevation angle e. The interferometric delay  is derived from the length difference between the reflected signal and the orthogonal projection of the direct signal to the reflected signal.
From the geometry of Figure 1 the familiar equation for the relation between the height of an GNSS antenna over the water surface, the reflector height H, and the interferometric delay  with respect to the elevation angle e can easily be derived =2Hsin(e) (1) The interferometric delay from eq. (1) governs the frequency of the SNR oscillation. Neglecting trends, signal amplitudes and attenuations, the oscillation is commonly expressed as a function of a phase offset  0 and the wave length of the GNSS signal λ by It is clear that Figure 1 presents a simplified model of the real geometry. Hence, eq. (1) will only hold under the assumption of different approximations that might yield systematic errors in the estimation of the reflector height. In the past, some of these approximations were examined and more realistic models were derived, yielding an increased quality of the results.
As one of the major approximations, Figure 1 assumes a horizontal plane reflector. This seems to be an adequate approximation for lower reflector heights and larger elevation angles. If low elevation angles (for example less than 5°) should be used, at least a spherical reflecting surface should be applied [10]. In this case, a more realistic incidence angle at the specular point can be calculated from an iteration [5].
The interferometric delay is calculated from eq. (1) as the difference of the length of the reflected signal path minus the one of the direct signal path, projected orthogonally to the direction of the reflected signal. This is only valid under the assumption of parallel signal paths, what would require a satellite at an infinite distance. In reality, the distance is finite and therefore the paths are not parallel. Although the elevation angle at the antenna and the incidence angle at the specular point are similar, they are not equal even for a plane reflector. The difference might be negligible for many applications but it can become important if another major approximation, the assumption of a rectilinear geometry, should be overcome. Figure 1 shows the geometry under pure vacuum conditions only. In reality, both paths pass through the neutral atmosphere and are therefore refracted while the propagation velocity along the paths is retarded. Hence, atmospheric corrections are mandatory to derive high-quality results. In the past, several studies considered atmospheric corrections in GNSS-R for example by using raytracing [11], [12], by application of an adaptive mapping function [13] or simply by taking into account the atmospheric effects below the antenna [14]. These corrections are not easy to integrate into the afore mentioned SNR analysis based on eq. (2). That might be a reason why some groups using these methods do not consider or even mention the atmospheric effect.
The authors of [15] suggested to correct for the atmospheric influence in a model based on eq. (1) using a correction of the elevation angle due to the bending of the refracted ray. In [16] it was remarked that this approach does not take into account the propagation delay and the authors suggested to use an interferometric atmospheric delay derived from mapping functions for the wet and the dry atmosphere together with a sophisticated atmospheric model. Lately, the authors of [17] used a rigorous atmospheric raytracing procedure and found that the atmospheric interferometric delay should be decomposed into an along-path delay, which results from the propagation delay along the bent ray paths and a vacuum atmospheric delay, that accounts for the difference of the interferometric delays from a vacuum ray path and a bent ray path. From numerical simulations, the authors found similar values for both components. In order to derive a formalism that could be used in practical SNR analysis, they developed a method based on rectilinear approximation of some of the bent ray paths [18].
The aim of this investigation is to evaluate the latest developments with regard to atmospheric corrections for SNR analyses that make use of eq. (2). Based on that, a simple model for the atmospheric interferometric delay will be suggested that takes into account ray bending as well as along-path propagation delay. In Section 2, we will revise the vacuum interferometric delay model presented in Figure 1 under the assumption of intersecting direct and reflected ray paths. Based on this, we will evaluate the geometric atmospheric delay from [17] in Section 3 by deriving a closed formulation for the bent ray interferometric delay model and by comparing it to the vacuum interferometric delay model. A simple interferometric delay model that incorporates ray bending and propagation delay will be deduced in Section 4. Our model will be tested using numerical experiments by means of a simple raytracing procedure in Section 5. Section 6 will conclude our findings.

Vacuum Ray Paths
To derive the vacuum interferometric delay assuming a satellite at a finite distance, we will keep to the idea of signals travelling in a pure vacuum. For the ease of understanding we will examine the interferometric delay for ground-based GNSS-R with low reflector heights and elevation angles larger than 5°, which allows us to deal with a planar reflector. Figure 1 must be modified so that the paths of the direct and reflected signal intersects at the satellite. It is clear that the elevation angle at the antenna e and the incidence angle at the specular point esp are not equal anymore. Therefore, we cannot use the orthogonal projection of the direct single path to that one of the reflected signal to calculate the interferometric delay. Instead, we have to rotate the direct signal by the intersecting angle  of the paths (Figure 2), that is the difference of esp and e.
Using this angle, the vacuum interferometric delay  is simply calculated as sin( ) ee 2H sin(e) 1 cot(e) tan 2 (4) From eq. (4), it is obvious that eq. (1) only holds, if the elevation angle at the antenna and the incidence angle at the specular point are equal because only in that case the factor in the bracket becomes 1.
To quantify the influence of , we can use typical values for the geometric elements. Let us assume a satellite at a distance of 25,000 km at an elevation angle of 5° from the antenna and a reflector height of 10 m. For these values we calculate an incidence angle esp of about 5.000045662° and from eq. (4) a vacuum interferometric delay of 1.743123 m. If we would apply eq. (1), the resulting interferometric delay would differ from the correct one by less than 1 m. Even for a reflector height of 100 m and an elevation angle of 1° the difference would be smaller than 1 mm. Therefore, it is reasonable to replace eq. (4) by the approximation from eq. (1) in the case of a pure vacuum.
For further evaluations it might be helpful to split the length of the reflected signal path into components above (Dr,a) and below (Dr,b) the horizon of the antenna. From Figure 2 we find = r ,a d sp sin(e) DD sin(e ) (5) Here, we cannot neglect the difference between the elevation angle and the incidence angle because the relation of the sin would become 1 if both angles would be equal. This could only be the case for an satellite at an infinite distance, since only in that case Dr,a would be equal to Dd. For the component below the antenna horizon we find Hence, we can express the vacuum interferometric delay in very good approximation from eq. (1) and eq. (6) as

Bent Ray Paths
The signals do not travel in a vacuum but in an atmosphere with variable density and will therefore experience a refraction according to Snell's law. If we assume a typical atmospheric structure where the density decreases with an increasing altitude, the signals will travel along curved paths. The length of the bent ray paths could be used to define and calculate a bent ray interferometric delay bent, too. The authors of [17] defined a geometric atmospheric delay as the difference of this delay and the vacuum interferometric delay from eq. (1) or (4). Since we assume an atmosphere for the ray bending but no propagation delay along the ray path, a bent ray interferometric delay is physical impossible [19], but might be useful in understanding the geometry.
The geometry of the bent ray paths is presented in Figure 3. The elevation angle and the incidence angle refer to the tangent of the bent ray paths As described by [15], the elevation angle is changed with respect to the vacuum conditions by a bending angle e of the direct signal. According to [17], the difference between the bending angle at the antenna and the specular point esp is a thousand times smaller than die bending angle itself and might be neglected so that esp could possibly be replaced by e. Figure 3. The bent ray paths from the satellite to the antenna and the specular point above the antenna horizon are not shown. For the component below the antenna horizon we assume straight lines instead of bent rays since the curvature is very small. The elevation angle e of the direct path between the antenna and the satellite defers from the bent ray path by the bending angle e. The bent ray path of the reflected signal passes through the specular point and the end of the mirrored antenna below the water surface. The bending angle esp at the specular point is spanned by the bent ray path of the reflected signal and the connection from the specular point to the satellite. The bent ray interferometric delay bent is derived similar to the vacuum case as the difference between the lengths of bent ray path of the reflected signal and that of the rotated bent ray path of the direct signal.
From Figure 3 it can be seen that the elevation angle e and the incidence angle esp refer to the straight lines connecting the satellite with the antenna and the specular point respectively. The bent ray paths defer from these angles by the bending angle. Since the curvature of the bent rays is rather similar because they pass through almost the same part of the atmosphere, we can derive the bent ray interferometric delay from a rotation of the bent ray of the direct signal as in the vacuum case. Again, we have to use the intersecting angle  of the strait line connections as the angle of rotation, but the intersection of the arc with the bent reflected ray path and the straight line defer now. The point of intersection of the arc with the bent reflected ray path can be found in very good approximation by the elongation of the chord to that path line.
From Figure 3 and this assumption, we can find the relation Considering eq. (4), this relation yields the bent ray interferometric delay as cos( e ) ee 1 tan( e ) tan 2 (9) Again, we can state as in the vacuum case that the influence of the difference between e and esp is of minor order and the approximated bent ray interferometric delay reads   bent sp sin(e) 2H cos( e ) (10) Here, we can calculate the component of the ray path of the reflected signal below the antenna horizon from From eq. (10) and eq. (11) we find the relation between this component and the bent ray interferometric delay as For a numeric evaluation we again use the setting from Section 2 and assume that the bending angle esp could possibly be replaced by e. We used Bennett's formula [20] and calculated the bending angle for a temperature of 23° C and a pressure of 1013 hPa as 0.1596° for an elevation angle of 5°. With these values, the bent ray interferometric delay from eq. (9) will become 1.743129 m. The difference between the vacuum interferometric delay from Section 2 and the bent ray interferometric delay is less than 7 m. A comparison of the interferometric delays from the approximation formulae eq. (1) and (9) yields a very similar value. Even for a reflector height of 100 m and an elevation angle of 1° the interferometric delay difference from eq. (1) and (9) is less than a tenth of a mm.
We can validate this result by a rough approximation. Let us replace the bent ray paths by arcs. The arc's chords should be the direct ray path Dd and the reflected ray path Dr respectively from the vacuum case. Let us further assume an angle between the arc's chord and the tangent at the end of the arc of twice the bending angle. The length of these arcs will be much larger than that of the bent ray paths due to the larger curvature of the arcs. If we calculate the bent ray interferometric delay from the approximating arcs, we get about 10 m for an elevation angle of 5° and a reflector height of 10 m and less than 1 mm for an elevation angle of 1° and a reflector height of 100 m Hence, we cannot confirm the results from [17] that shows values of about 6.5 cm for their geometric atmospheric delay for a reflector height of 10 m and an elevation angle of 5°. It is remarkable that a value of about 5.5 cm results with our setting, if we would apply the aforementioned simplification of parallel ray paths. In that inaccurate case, we would have to project the bent ray path of the direct signal orthogonal to that one of the reflected signal. To do so, we have to use e+e instead of e in eq. (1). Taking the bending angle of about 0.185° from [17], we end up with a difference of 6.4 cm for parallel ray paths. However, it is important to take into account the intersection of the bent ray paths at the satellite, although it has no major impact on the vacuum interferometric delay.

Atmospheric Ray Paths
In Section 3, we derived the component of the bent ray path of the reflected signal below the antenna horizon while neglecting the propagation delay. Because the component is commonly small for low antenna heights, we already approximated them by a straight line. Hence, we can simply account for the retardation of the propagation velocity of radio waves be multiplying this component by an index of refraction nb for the atmosphere below the antenna horizon.
For the component above the antenna horizon as well as for the bent ray path of the direct signal, we cannot use this simplification. As can be seen from Figure 3, the piercing point of the vacuum path of the reflected signal at the antenna horizon defers from that one of the bent ray path of the reflected signal. To avoid a break in the ray path, we have to follow the tangent of the bent ray path of the reflected signal until it intersects with the tangent of the bent path of the direct signal. Again, we have to project the direct path to the reflected path by a rotation. If we assume that the bending angles at the antenna and the specular point are almost equal, we can once more use the intersecting angle  of the strait line connections as the angle of rotation, but the center of rotation defers. Here, we have to use the intersection of the tangents of the bent rays rather than the satellite. After that, both path lengths might be multiplied by an average index of refraction and used for calculating the atmospheric interferometric delay atmo.
As we have seen from Section 2, for typical geometrical settings in GNSS-R we can map the direct path by an orthogonal projection in very good approximation of rotational projection to the reflected path and end up with an almost same value for the atmospheric interferometric delay atmo.
For the geometry below the antenna horizon, the change of the center of rotation or the orthogonal projection is likewise important. Figure 4 shows that we can calculate the interferometric delay in the same way as in the vacuum case, but we do have to take into account the change of the elevation angle and replace e by e+e and in esp by esp+esp in eq. (4), (6) and (7). To avoid a break in ray paths above the antenna horizon, the bent paths above the antenna horizon have to be rotated, whereby the center of rotation results from the intersection of the tangents of the bent paths. Again, the interferometric delay can be calculated from the triangle spanned by the chord (green line) and the reflector height H. It should be mentioned, that the propagation delay is still neglected in this figure.
Hence, the atmospheric delay can be expressed by the component below the antenna horizon from an adaptation of eq. (7) n D sin(e e )sin(e e) (13) The component below the antenna horizon can likewise be derived from a modification of eq. (6), but that was already done in eq. (11) in Section 3. Together with eq. (13) we can finally calculate the atmospheric interferometric delay from This is the well-known formulation from [15] but accounting for the propagation delay in addition to the bending model. The modification might seem small but is important. Let us imagine a theoretical observation with an elevation angle of 90°, at which the bending angle vanishes. Although we would not observe any reflection in reality, we could calculate the theoretical interferometric delays. The formulation from [15] would end up in the vacuum interferometric delay from eq. (1) as simply twice the reflector height. Because the reflected signal would pass the atmosphere twice, the atmospheric interferometric delay must be twice the reflector height, multiplied by the index of refraction in the lower part of the atmosphere. That is exactly what eq. (14) yields.

Numerical Experiment
We examined our findings a by numerical experiment based on a simple step-by-step raytracing procedure as described in [19]. We assumed a spherical earth with a radius of 6378137 m and a satellite at an altitude of 20,000 km above the earth. We defined a 2D coordinate system starting at the center of the sphere. The vertical axis was set to pass through the antenna. The reflector height was set to 10 m and the elevation angles range from 1° to 90°. From the altitude of the satellite and the radius of the spherical earth we calculated the coordinates of the satellite for all elevation angles.
The atmosphere was approximated as a layered spherical structure with a layer increment of 10 m. The index of refraction was taken from the same CIRA-86 model [21] and calculated in the same manner as in [17], whereby the pressure below 20 km was log-linear interpolated, so that it fits to the CIRA-86 value at an altitude of 20 km and to 1013.15 hPa at the ground.
We applied an iterative computation of the raytracing. The rays were calculated in the inverse direction, what means that we started at the antenna or specular point and computed step-by-step the piercing point of the ray with the upper-nearest atmospheric layer, taking into account Snell's law to derive the deflection of the ray at this layer limit. Above the top layer of the atmospheric model in an altitude of 120 km we assumed a constant index of refraction of 1, and therefore, a straight line as the last ray.
The initial elevation angle of the bent path was set to the vacuum elevation angle. The perpendicular distance of the satellite from the last ray was used to derive a correction for the initial elevation angle and applied in the next iteration step. The iteration stopped when the last ray passed the satellite within a range of a tenth of a millimeter.
For the ray path of the reflected signal we calculated the incidence angle according to [5] and combined it with the bending angle from the raytracing of the direct signal to compute the coordinates of the specular point. The 2D coordinate system was rotated so that the vertical axis passes through the specular point and the same iteration as for the direct signal was conducted. After the raytracing iteration was finished, the coordinate system was rotated back. The resulting elevation angle of the bent path of the reflected signal was used to recalculate the position of the specular point and the iterative raytracing was repeated. The iteration of the coordinates of the specular point was stopped when the change of coordinates was smaller than 1 millimeter.
The raytracing yields the geometric length of the paths. For the computation of the radio length, we used the ray parts between two layers as finite differences and the mean index of refraction to derive the radio length of the paths from numerical integration.
Hence, we obtained from the raytracing the bending angle, the vacuum interferometric delay, the bent ray interferometric delay and the atmospheric interferometric delay. In a first step, the bending angle from our raytracing was compared to that one from Bennett's formula. Figure 5 depicts that the bending angles agree very well for elevation angles above 5°. Larger discrepancy for lower elevation angles may result from different atmospheric models applied here and in Bennett's development. Hence, it can be stated that the raytracing procedure yields reliable results. Next, we compared the bent ray interferometric delay and the vacuum interferometric delay, both derived from raytracing. The differences (blue line presented in Figure 6) are less than 1 cm for all elevation angles. This comparison confirms our findings from Section 3. For the lowest elevation angle, this difference is more than ten times smaller than the difference between atmospheric interferometric delay and the vacuum interferometric delay (orange line in Figure 6). The latter can be compared to the along-path-delay from [17]. For an elevation angle of 5° the authors found an along-path-delay of about 6.9 cm. The difference of the atmospheric interferometric delay and the vacuum interferometric delay from our raytracing is about 5.2 for that elevation angle. The discrepancy results from the difference in the bending angles. From our raytracing we derived a bending angle of 0.149. The bending angle from [17] is about 0.185°. Applying eq. (14) for both bending angles results in a difference of the atmospheric interferometric delays of about 1.3 cm. Figure 6. Differences of the bend ray interferometric delay and atmospheric interferometric delay respectively to the vacuum interferometric delay from raytracing.
Finally, we compared the atmospheric interferometric delay from raytracing and that derived from eq. (14) for the bending angle from raytracing as well as from Bennett's formula. Figure 7 demonstrates the high quality of eq. (14) since the absolute differences are less than 1 mm for all elevation angles if we use the bending angle from raytracing. The differences for the case when the bending angle from Bennett's formula is uses in eq. (14) are quite large for lower elevation angles. They become less than 1 mm only for elevation angles larger than about 12°.

Figure 7.
Absolute differences of the atmospheric interferometric delay from raytracing and from the approximation from eq (14) using the bending angle from raytracing and Bennett's formula respectively.
The variation of the index of refraction of the atmosphere below the antenna is likewise important. Figure 8 shows the difference between the atmosphere for an index of refraction from the CIRA-96 model and pure vacuum for that part of the atmosphere only. The differences increase with an increasing elevation angle with a maximum of 2H(nb-1) for an elevation angle of 90°. In the typical range of the elevation angles used in ground-based GNSS-R of about 5° to 30°, the differences are almost larger than that from Figure 7 for the bending angle from Bennett's formula. This leads to the conclusion that both the upper and the lower part of the atmosphere should be modeled well.

Conclusions
We examined the atmospheric modeling in relation to the analysis of SNR data from groundbased GNSS-R observations from a geometric point of view. This was completed by a numerical test applying a simple raytracing.
We revised the vacuum interferometric delay model and extended it for intersecting vacuum ray paths. The quantification of the influence of the intersecting angle shows that the assumption of parallel rays in a vacuum is appropriate for the typical settings of ground-based GNSS-R observations.
We used the deductions from the vacuum case to evaluate the case of physically impossible refracted rays in a vacuum. A closed formula for the bent ray interferometric delay, what is the difference of the geometric lengths of the bent ray of direct and the reflected signal, was derived. The comparison with the vacuum interferometric delay showed even for larger reflector heights and low elevation angles non-significant differences. Hence, we cannot confirm the results from other groups.
Taking into account the retardation of the propagation velocity of radio waves in non-vacuum conditions yields atmospheric ray paths. Based on the preceding findings we derived a relation between the atmospheric interferometric delay and the component of the refracted ray path of the reflected signal below the antenna horizon. The final formulation of the atmospheric interferometric delay is an extension of a well-known formula.
We compared the theoretical results by calculating the various path lengths and delays from a simple raytracing, using a typical atmospheric model. The comparison of the bending angle of the direct signal path from this raytracing with a standard formula showed good agreement. Hence, it seems reasonable to assume that the results from our raytracing are reliable. The atmospheric interferometric delay from our formulation agrees very well with that resulting from the raytracing for all elevation angles.
The evaluation of our formula for the atmospheric interferometric delay shows that both the atmospheric layer above and below the antenna horizon should be modeled well. Since the layer above the antenna horizon influence only the bending angle, the modelling of that part might be less important for larger elevation angles. The layer below the antenna horizon influences the atmospheric interferometric delay also for larger elevation angle. Hence, we recommend to include the humidity besides temperature and pressure in the computation of the index of refraction similar to [16] since it might show a strong variability, especially over water. In the future, studies on the behavior of atmosphere over water surface as suggested by [22] might benefit from eq. (14), too.