Space-Borne GNSS-R Ionospheric Delay Error Elimination by Optimal Spatial Filtering

Global Navigation Satellite System Reflectometry (GNSS-R) technology is a new and promising remote sensing technology, especially satellite-based GNSS-R remote sensing, which has broad application prospects. In this work, the ionospheric impacts on space-borne GNSS-R sea surface altimetry were investigated. An analysis of optimal values for spatial filtering to remove ionospheric delays in space-borne GNSS-R altimetry was conducted. Considering that there are few satellite-borne GNSS-R orbit observations to date, simulated high-resolution space-borne GNSS-R orbital data were used for a comprehensive global and applicable study. The curves of absolute bias in relation to the bilateral filtering points were verified to achieve the minimum absolute bias. The optimal filtering points were evaluated in both statistical probability density and quantile analysis to show the reliability of the selected values. The proposed studies are helpful and valuable for the future implementation of high-accuracy space-borne GNSS-R sea surface altimetry.


Introduction
Global Navigation Satellite System Reflectometry (GNSS-R) technology is a promising remote sensing technology, especially satellite-based GNSS-R remote sensing, which has broad application prospects [1,2]. Compared with traditional microwave remote sensing technology, GNSS-R has the advantages of low power consumption, light weight, small size, low cost, and passivity. Since it can receive reflection signals from multiple satellites at the same time, its spatial resolution is higher than that of traditional microwave remote sensing technology, and its data can fill some areas that cannot be covered by traditional remote sensing technology. Currently, GNSS-R has been widely used for the remote sensing of sea surface winds and heights, soil assessments, and ice status monitoring.
GNSS-R sea surface altimetry is currently a research focus. At present, a large number of ground and airborne GNSS-R altimetry experiments have been carried out, which have fully verified the effectiveness of this technology [3][4][5][6]. The development of ground and airborne GNSS-R altimetry will promote the application of this technology in satellite platforms. Compared to traditional sea surface altimetry, this technology can cover a vast region of the sea and help measure the sea surface height from a global view. Space-borne GNSS-R sea surface altimetry uses the single difference between the reflected navigation signal from the sea surface and the direct signal from the navigation satellite simultaneously received by the GNSS-R receiver loaded on a low earth orbit (LEO) satellite to retrieve the sea surface height. For the first time, researchers can use TDS-1 data to retrieve sea surface height experiments [7]. For the two different regions considered, the difference between themean sea surface height (MSSH) of DTU10 [8] is 8.1 and 7.4 m, which is within a tolerable level. The corresponding influence that affects sea surface altimetry accuracy was provided after the performance analysis of Existing GNSS-R satellites, such as TDS-1 and CYGNSS, use only single frequency (L1) measurements. Under such conditions, the ionospheric model is necessary to retrieve ionospheric delays for single frequency measurements. A series of models were proposed and accepted for ionospheric studies. Here, the GIM model is applied in this work to estimate the ionospheric delay for space-borne GNSS-R observations. The GIM is a common model used to study global ionospheric changes and provides global ionosphere total electron content (TEC) [18][19][20]. The global ionosphere has been divided with a longitude step of 5° and a latitude step of 2.5°, making for a total of 5183 grids globally. The GIM is provided by IGS data analysis centres with temporal intervals of 2 hours or even less. The VTEC at the ionospheric pierce points shown in Figure 1 was calculated with the GIM to obtain the VTEC from the GNSS transmitter to the sea surface reflecting point (VTEC1), the sea surface reflecting point to the space-borne GNSS-R receiver (VTEC2) and the STEC from the GNSS transmitter to the space-borne GNSS-R receiver (STEC3). In further calculations, the effect of STEC3 is very minor and can be ignored; thus, only STEC1 and STEC2 were used. For the path of STEC3, the altitude of ionosphere is larger than 600 km, according to the vertical structure of ionosphere, the electron density is mainly concentrated below 450 km in the F2 layer, as revealed by many electron density observations [21]. In high altitude of ionosphere (600-1000 km and above), the impact of ionosphere on GNSS signal is relatively weak.
Once VTEC1 and VTEC2 were calculated, a mapping function was used to project the VTEC to STEC along the signal propagation path, the calculation was represented by Equations (1)-(4) [8]. The mapping function is given as: where is the elevation of the satellite at the sea surface reflection point, is the radius of the Earth, and h is the height of the ionosphere used in the GIM model. The total ionospheric delay was computed with the combination of STEC1 and STEC2 and is given as: Figure 1. In the GNSS-R ionospheric error model, the ranging error of the ionosphere is estimated for the signals reflected from the Earth using the global ionosphere map (GIM) model with data from International GNSS Service (IGS), and the vertical total electron content (VTEC) value is obtained and then mapped to slant total electron content (STEC) based on different elevations for a single point.
Existing GNSS-R satellites, such as TDS-1 and CYGNSS, use only single frequency (L1) measurements. Under such conditions, the ionospheric model is necessary to retrieve ionospheric delays for single frequency measurements. A series of models were proposed and accepted for ionospheric studies.
Here, the GIM model is applied in this work to estimate the ionospheric delay for space-borne GNSS-R observations. The GIM is a common model used to study global ionospheric changes and provides global ionosphere total electron content (TEC) [18][19][20]. The global ionosphere has been divided with a longitude step of 5 • and a latitude step of 2.5 • , making for a total of 5183 grids globally. The GIM is provided by IGS data analysis centres with temporal intervals of 2 h or even less. The VTEC at the ionospheric pierce points shown in Figure 1 was calculated with the GIM to obtain the VTEC from the GNSS transmitter to the sea surface reflecting point (VTEC1), the sea surface reflecting point to the space-borne GNSS-R receiver (VTEC2) and the STEC from the GNSS transmitter to the space-borne GNSS-R receiver (STEC3). In further calculations, the effect of STEC3 is very minor and can be ignored; thus, only STEC1 and STEC2 were used. For the path of STEC3, the altitude of ionosphere is larger than 600 km, according to the vertical structure of ionosphere, the electron density is mainly concentrated below 450 km in the F2 layer, as revealed by many electron density observations [21]. In high altitude of ionosphere (600-1000 km and above), the impact of ionosphere on GNSS signal is relatively weak.
Once VTEC1 and VTEC2 were calculated, a mapping function was used to project the VTEC to STEC along the signal propagation path, the calculation was represented by Equations (1)-(4) [8]. The mapping function is given as: where El is the elevation of the satellite at the sea surface reflection point, R E is the radius of the Earth, and h is the height of the ionosphere used in the GIM model. The total ionospheric delay was computed with the combination of STEC1 and STEC2 and is given as: ( The total ionospheric delay is finally represented as: In the space-borne GNSS-R altimetry model, the height error corresponding to the ionospheric delay is expressed as: where θ is the incident angle of the sea surface reflecting point. To use the simulated orbital observation, the ionospheric delay was simulated as: where ∆ρ IonoNoise is Gaussian noise with standard deviation varying in accordance with the incident angle, which is given as: Here, supposing σ h = 20 cm, that is, the influence of ionospheric delay observation noise on height accuracy is 20 cm, then the standard deviation of Gaussian noise will vary with the incident angle of the specular reflection point.

Spatial Filtering of Ionospheric Delay
In space-borne GNSS-R altimetry, the calculation was represented from Equations (7)-(11) [13], the single difference of measurement is given as: where h means the height measured from GNSS-R altimetry. I means the ionospheric delay related to θ, the incidence angle of the single point, f means the frequency of the signal. This formula can be rewritten as: where y = ρ/(2cosθ), x = I /(2cosθ). If the ionospheric delay is corrected by a dual frequency method (assuming L1 and L5 bands are used), the above equation can be expressed as follows at different frequencies: The error propagation function was derived by solving the above equation, which is given as: Sensors 2020, 20, 5535

of 17
Since the ionospheric delay is spatially related to itself and not related to the topography of the sea, it can be considered that the ionospheric delay changes little (or is smooth) on a limited spatial scale. Therefore, the accuracy of the ionospheric delay can be improved by a regression of multiple continuous sampling points. If the regression is accurate, the accuracy of ionospheric delay estimation will be improved by √ n [13], signal n means the number of points to be used for filtering. If the dual frequency method is used to correct the ionospheric delay and spatial filtering is carried out, the final altimetry measurement accuracy can be expressed as follows: This spatial filtering method is based on two assumptions: (1) the variation of ionospheric delay along an orbit is gentle; (2) since the variation of ionospheric delay is predictable, the accuracy can be improved by effective regression. The proposed spatial filtering method can help reduce the ionospheric delay for space-borne GNSS-R propagating signals, resulting in improvements in altimetry measurement accuracy.
The GNSS-R simulated trajectory data were analyzed by the spatial filtering method in the statistical aspect, the simulated data were provided by the National Space Science Center of China, with over 2 million GNSS-R reflection trajectories in global coverage. The simulated GNSS-R trajectories were quite similar to the real condition, to supplement the vacancy of real space-borne GNSS-R observations for analysis. In the proposed method, the GNSS-R ionospheric delays were derived by Equation (5) with simulated trajectories, then the mean value filter method was used to compare the ratio of MAE (mean absolute error). The optimal filtering performance was decided with the variation curve of filtering points f n. The experiments were conducted with over 2 million global simulated GNSS-R trajectories to decide the distribution of optimal filtering points. The results were evaluated with two criteria: to find the peak value of f nmin and to use 99% quantile for f n. Compared with polynomial regression, the mean value filtering method has good computation efficiency over the polynomial method, and has better accuracy under small and medium spatial scales.

Ionospheric Delay from Simulated Orbital Data
Simulated orbital data are used in this experiment, including the Global Positioning System (GPS) satellite trajectory, reflection point trajectory, space-borne receiver trajectory, elevation, azimuth, and incident angle of the reflection point. The ionospheric delay of each sea surface reflection point signal propagating path can be calculated according to the method in Section 2.1. Then, a single frequency ionospheric delay correction of the reflection point trajectory is obtained from the GIM model, and noise is added to the ionospheric delay correction to simulate ionospheric delay measurements. Figure 2 shows the trajectories of reflection points distributed in eight different regions of the world and the corresponding single frequency ionospheric delay measurements. The red curves in Figure 2 show the single frequency ionospheric delay corrections from the GIM, and the blue curves show the ionospheric delay measurements with observation noise.
Then, a single frequency ionospheric delay correction of the reflection point trajectory is obtained from the GIM model, and noise is added to the ionospheric delay correction to simulate ionospheric delay measurements. Figure 2 shows the trajectories of reflection points distributed in eight different regions of the world and the corresponding single frequency ionospheric delay measurements. The red curves in Figure 2 show the single frequency ionospheric delay corrections from the GIM, and the blue curves show the ionospheric delay measurements with observation noise.

Spatial Filtering Method Results for Ionospheric Delay
The simulated ionospheric delay series for different reflection point trajectories were first analyzed by filtering the mean values. The steps of this analysis were introduced as follows: (1) select a trajectory point named SP0, the point to be estimated or filtered, on a certain reflection trace; (2) take 100 continuous ionospheric delay measurements before and after the point to form an analysis sample with a total of 201 reflection points; (3) take the point SP0 as the center and select f n samples before and after the center separately, where the mean value of these samples is considered the estimation of the center value; f n samples means the number of points before and after the center separately to be used for filtering, called filter n( f n); (4) study the variation of absolute bias in relation to the changes of f n, the absolute bias is the difference between the estimation and the real value derived from the single frequency ionospheric delay correction model. The results are shown in Figure 3. In the results, Figure 3a shows the sample trajectory of SP0; Figure 3b shows the ionospheric delay measurements and correction from the GIM model; Figure 3c shows the variation curve with bilateral filtering distance f n of the deviation σ, the difference between the estimated value of the mean filter and the model value (ρ IonoMod ); it demonstrates the sequence of 30 dual-frequency ionospheric delay observations for the same orbit trajectory. The observations were derived from model plus noise, containing certain randomity. Figure 3d shows the sequence of 30 ionospheric delay dual frequency observation values for the same orbit trajectory because the observation values are obtained by the modeled value plus noise, so there is certain randomness. Through multiple filtering analysis, the randomness is eliminated; that is, the filtering curve in Figure 3c is the average value of 30 ionospheric delay observation values passing through the same point.
Sensors 2020, 20, 5535 7 of 17 ionospheric delay observations for the same orbit trajectory. The observations were derived from model plus noise, containing certain randomity. Figure 3d shows the sequence of 30 ionospheric delay dual frequency observation values for the same orbit trajectory because the observation values are obtained by the modeled value plus noise, so there is certain randomness. Through multiple filtering analysis, the randomness is eliminated; that is, the filtering curve in Figure 3c is the average value of 30 ionospheric delay observation values passing through the same point. 30 times with the same method to decrease randomity brought by the noise, 30 times were efficient for computation consumption; (d) deviation plot (filtering plot) with between the estimated and modeled ionospheric error (true value for comparison), with a parameter value σ / 0 . and 0 indicate the mean absolute error after and before filtering.
During the mean filtering process, the standard deviation of absolute bias was denoted as σ, which was computed by the following equation:  During the mean filtering process, the standard deviation of absolute bias was denoted as σ which was computed by the following equation: where ρ ionoEstimated ( f n) means the value of ionospheres delay after noise being filtered using f n ( 2 f n + 1 points around the single point to be analyzed) with the filtering method mentioned above, ρ ionoMod ( f n) is the value of ionospheres delay as a contrast. The step of summing means the average of 30 times estimating. Figure 4 shows two typical mean filtering curves; with an increment of the sampling points f n, the σ has two variation trends: (a) σ decreases with fast increment of f n; after reaching to a bottom value, σ gradually increases. Thus, an optimal f n can be found to approach the best estimation between the real value and the estimated value, called f nmin. (b) σ decreases drastically with increasing f n and then decreases gradually. Under this condition, an optimal f n can be found to approach the best estimation between the real value and the estimated value. To summarize, an optimal sampling point f nmin can be obtained for ideal filtering.
the has two variation trends: (a) decreases with fast increment of ; after reaching to a bottom value, gradually increases. Thus, an optimal can be found to approach the best estimation between the real value and the estimated value, called min. (b) decreases drastically with increasing and then decreases gradually. Under this condition, an optimal can be found to approach the best estimation between the real value and the estimated value. To summarize, an optimal sampling point can be obtained for ideal filtering. It can be seen from the above results that the accuracy of ionospheric delay can be greatly improved if is found after mean filtering. Without loss of generality, a global view of spatial filtering estimation was conducted as follows to obtain the relationship between spatial parameters of longitude, latitude, incident angle, and the optimal for spatial filtering.

Statistical Analysis of the Spatial Filtering Method
For this purpose, with the applicable GIM global grids, the reflection points were analyzed from two aspects. First, the global distribution analysis of was performed by replacing the mean value of in each grid with the values for all regions in each grid; second, a statistical distribution of was performed for all the analysed reflection points. Each spatial filtering curve was studied, as shown in Figure 5.
All the reflection points were analysed with the GIM global grids. The average value of for each grid point replaced all the values in the grid, and the global distribution for the mean value of was then analyzed. The statistical distribution for the optimal filtering points and of all the reflection points were studied. The statistical distribution for was analyzed when different were considered for analysis.  It can be seen from the above results that the accuracy of ionospheric delay can be greatly improved if f nmin is found after mean filtering. Without loss of generality, a global view of spatial filtering estimation was conducted as follows to obtain the relationship between spatial parameters of longitude, latitude, incident angle, and the optimal f nmin for spatial filtering.

Statistical Analysis of the Spatial Filtering Method
For this purpose, with the applicable GIM global grids, the reflection points were analyzed from two aspects. First, the global distribution analysis of f nmin was performed by replacing the mean value of f nmin in each grid with the values for all regions in each grid; second, a statistical distribution of f nmin was performed for all the analysed reflection points. Each spatial filtering curve was studied, as shown in Figure 5.   All the reflection points were analysed with the GIM global grids. The average value of f nmin for each grid point replaced all the values in the grid, and the global distribution for the mean value of f nmin was then analyzed. The statistical distribution for the optimal filtering points f nmin and σ( f nmin) of all the reflection points were studied. The statistical distribution for σ( f n) was analyzed when different f n were considered for analysis. Figure 5 shows the global distribution of the mean f nmin. Each grid value represents the mean value of all reflection points in the grid. The up plot demonstrates the global VTEC provided by IGS, the VTEC values were differentiated along the latitudinal direction, and the absolute values were considered. The bottom plot shows the global distribution of the mean f nmin. The absolute differentiated VTEC and mean f nmin were negatively correlated, especially in the peak regions of ∂VTEC ∂lat . Each spatial filtering curve was studied, as shown in Figure 6. Most values of f nmin (80%) were concentrated in a range of f n < 55, BDS constellation data for GNSS-R f n < 55 ; with peak probability when f n = 30 for GPS constellation data for GNSS-R, f n = 28 for BDS-R data; and fewer data were distributed in a range between 55 and 100 for both GPS-R and BDS-R data. However, the probability increased drastically when f n = 100, this is because the overall filtering length is 201, and f n = 100 is the maximum median filtering length. At the same time, the standard deviation of absolute bias decreased greatly with the increment of f n and then decreased gradually, in agreement with Figure 3b. For the statistical analysis of σ min , it was found that 99% of the σ min was less than 0.17σ 0 , 95% of the σ min was less than 0.14σ 0 , and the peak frequency values of the σ min were 0.09σ 0 .  Each spatial filtering curve was studied, as shown in Figure 6. Most values of (80%) were concentrated in a range of < 55 , BDS constellation data for GNSS-R < 55 ; with peak probability when = 30 for GPS constellation data for GNSS-R, = 28 for BDS-R data; and fewer data were distributed in a range between 55 and 100 for both GPS-R and BDS-R data. However, the probability increased drastically when = 100, this is because the overall filtering length is 201, and = 100 is the maximum median filtering length. At the same time, the standard deviation of absolute bias decreased greatly with the increment of and then decreased gradually, in agreement with Figure 3b. For the statistical analysis of , it was found that 99% of the was less than 0.17 , 95% of the was less than 0.14 , and the peak frequency values of the were 0.09 . From the above analysis, it can be concluded that when filtering ionospheric delay observations, there is an optimal for each reflection point that provides the best filtering. However, is a random variable and cannot be predicted through the relevant parameters of reflection points. To achieve a better filtering effect, the peak value of should be selected. When ∆h = 20 cm, = 29; under such conditions, a greater probability can be obtained to obtain a better filtering effect. Here, if all ionospheric delay samples are filtered with = 29 on both sides, the results From the above analysis, it can be concluded that when filtering ionospheric delay observations, there is an optimal f nmin for each reflection point that provides the best filtering. However, f nmin is a random variable and cannot be predicted through the relevant parameters of reflection points. To achieve a better filtering effect, the peak value of f nmin should be selected. When ∆h = 20 cm, f nmin = 29; under such conditions, a greater probability can be obtained to obtain a better filtering effect. Here, if all ionospheric delay samples are filtered with f nmin = 29 on both sides, the results are shown in Figure 7. Moreover, for ∆h = 10, 20, 30 cm, the f nmin distributions were demonstrated, and the corresponding peak values of f nmin were 21, 29, and 31. Figure 7 shows the statistical probability distribution of σ( f nmin) with different ∆h values. Each σ( f nmin) value is divided by 2 to compare the change in noise accuracy before and after filtering, and 2 represents the noise accuracy before filtering. Figure 8 shows the distribution σ( f nmin)/σ 0 when ∆h = 10, 20, 30 cm. From the above analysis, it can be concluded that when filtering ionospheric delay observations, there is an optimal for each reflection point that provides the best filtering. However, is a random variable and cannot be predicted through the relevant parameters of reflection points. To achieve a better filtering effect, the peak value of should be selected. When ∆h = 20 cm, = 29; under such conditions, a greater probability can be obtained to obtain a better filtering effect. Here, if all ionospheric delay samples are filtered with = 29 on both sides, the results are shown in Figure 7. Moreover, for ∆h = 10, 20, 30 cm , the distributions were demonstrated, and the corresponding peak values of were 21, 29, and 31. Figure 7 shows the statistical probability distribution of σ fnmin with different ∆h values. Each value is divided by 2∆h cos θ to compare the change in noise accuracy before and after filtering, and 2∆h cos θ represents the noise accuracy before filtering. Figure 8 shows the distribution / when ∆h = 10, 20, 30 cm. If the peak probability for the distribution is selected during filtering, that is, all points adopt bilateral filtering distance = 29, then the distribution for all filtered = 29 / is shown in Figure 9.  former is 0.1 σ ; the latter is 0.4 σ . That is, although the filtering result falling within 0.1 σ is the largest, the filtering result is greater than 0.2 σ in nearly 5% of cases, and the filtering result is greater than 0.4 σ in 1% of cases. The fluctuation range of the filtering effect is large. Therefore, according to the method of directly taking , = 15~40 can be selected to conduct filtering analysis on all reflection points and observe the distribution of / when takes different values to search for a more stable value of .  Figure 9. The probability distribution for = 29 / in all samples when the noise accuracy is ∆ℎ = 20 and = 29 for all samples to be filtered. Figure 10 shows the probability distribution for = 29 / when the accuracy is 20 cm, and If the peak probability for the f nmin distribution is selected during filtering, that is, all points adopt bilateral filtering distance f nmin = 29, then the distribution for all filtered σ( f n = 29)/σ 0 is shown in Figure 9. Here, the peak values of the probability distribution are 0.1 σ 0 , 0.95 quantile 0.21 σ 0 , and 0.99 quantile 0.43 σ 0 . The 0.95 and 0.99 quantiles can be compared with the values of 2σ and 3σ in the Gaussian distribution. It was found that when the peak value of f nmin = 29 is selected to filter the ionospheric delay observations for all reflection points, the filtering effect is not ideal because the peak for σ( f n = 29)/σ 0 is nearly 3 times different from that for the 0.99 quantile, and the former is 0.1 σ 0 ; the latter is 0.4 σ 0 . That is, although the filtering result falling within 0.1 σ 0 is the largest, the filtering result is greater than 0.2 σ 0 in nearly 5% of cases, and the filtering result is greater than 0.4 σ 0 in 1% of cases. The fluctuation range of the filtering effect is large. Therefore, according to the method of directly taking f nmin, f n = 15 ∼ 40 can be selected to conduct filtering analysis on all reflection points and observe the distribution of σ( f n)/σ 0 when f n takes different values to search for a more stable value of f n.    Figure 11 shows the 99% quantile variation curve in relation to different values. This result shows that this sampling point selection method is more stable.  Figure 10 shows the probability distribution for σ( f n = 29)/σ 0 when the accuracy is 20 cm, and f n is selected for different values in a range from 13 to 21. It can be seen that the 0.99 quantile of σ( f n = 29)/σ 0 varies when f n takes different values. That is, under the confidence level of 99%, the accuracy improvement after filtering is different. When f n = 17, the 99% quantile of σ( f n = 29)/σ 0 reaches a minimum value of 0.195, corresponding to an accuracy improvement of 0.195. Figure 11 shows the 99% quantile variation curve in relation to different f n values. This result shows that this sampling point selection method is more stable.  In this subsection, two methods of filtering point selection were discussed to obtain the optimal filtering. One is to select the point with the highest probability of ( peak), which creates a greater probability of obtaining the optimal filtering; the other is to select the appropriate to get a better 99% quantile of σ fn /σ , which makes the filtering effect more robust. A large number of epochs of the GIM data were used for modeling the ionosphere delay, the data were selected with 300 days span in 2017. The spatial mean value filtering method was applied with two optimal filtering criteria suggested above, results were demonstrated in Figure 12 and Figure 13. Figure 12a shows peak value is almost stable for the different time span with minor fluctuation. The histogram of peak value is described in Figure 12b. The optional value of the filtering point is = 20,28, and 32 corresponding to ∆ℎ = 10, 20, and 30 . In this subsection, two methods of filtering point selection were discussed to obtain the optimal filtering. One is to select the point with the highest probability of f nmin ( f nmin peak), which creates a greater probability of obtaining the optimal filtering; the other is to select the appropriate f n to get a better 99% quantile of σ( f n)/σ 0 , which makes the filtering effect more robust.
A large number of epochs of the GIM data were used for modeling the ionosphere delay, the data were selected with 300 days span in 2017. The spatial mean value filtering method was applied with two optimal filtering criteria suggested above, results were demonstrated in Figures 12 and 13. Figure 12a shows f nmin peak value is almost stable for the different time span with minor fluctuation. The histogram of f nmin peak value is described in Figure 12b. The optional value of the filtering point is f nmin peak = 20, 28, and 32 corresponding to ∆h = 10, 20, and 30 cm.
to get a better 99% quantile of σ fn /σ , which makes the filtering effect more robust. A large number of epochs of the GIM data were used for modeling the ionosphere delay, the data were selected with 300 days span in 2017. The spatial mean value filtering method was applied with two optimal filtering criteria suggested above, results were demonstrated in Figure 12 and Figure 13. Figure 12a shows peak value is almost stable for the different time span with minor fluctuation. The histogram of peak value is described in Figure 12b. The optional value of the filtering point is = 20,28, and 32 corresponding to ∆ℎ = 10, 20, and 30 . In Figure 13, it is found that the value of . varies with different epoch, indicating that this method is unstable in temporal domain. However, when the filtering point was limited less than 10, and ∆ℎ = 10, 20, and 30 , the values of σ fn /σ were almost the same.  In Figure 13, it is found that the value of . varies with different epoch, indicating that this method is unstable in temporal domain. However, when the filtering point was limited less than 10, and ∆ℎ = 10, 20, and 30 , the values of σ fn /σ were almost the same.   In Figure 13, it is found that the value of f n minQ0.99 varies with different epoch, indicating that this method is unstable in temporal domain. However, when the filtering point was limited less than 10, and ∆h = 10, 20, and 30 cm, the values of σ( f n)/σ 0 were almost the same. f n = 10 was used ass the optional value to get a better 99% quantile of σ( f n)/σ 0 . Under such condition, the valve of σ Q0.99 was same with different ionospheric epoch, indicating that the distribution of σ( f n = 10) was almost same; the value of σ Q0.99 decreased greatly when f n < 10. In the result, when f n = 10, σ Q0.99 = 0.3, it means that the error after filtering has 99% probability less than 0.3 times of the error before filtering.
Two methods of filtering point selection are discussed and optimized to obtain the optimal filtering effect with enough validity. The optional filtering point with method are 20,28, and 32 (about 459 km, 638.2 km, and 706 km along the track of single point, the sample frequency is 2 second along the track. In addition, f n is bilateral) corresponding to ∆h = 10, 20, and 30 cm. The elective filtering point with method is 10 (about 236 km along the SP track) corresponding to ∆h = 10, 20, and 30 cm.

Discussion
In the global distribution of f nmin, an obvious correlation was noticed between the mean grid value of f nmin and the absolute value of the zonal differential of the vertical total electron content in this grid. This is because most GNSS-R satellites are inclined orbit satellites, which leads to a single orbit on the earth reference ellipsoid surface. The variation in VTEC along latitude directly determines the change in ionospheric delay in the latitude direction. The correlation coefficient between the mean f nmin and the absolute zonal differential of VTEC was computed to be 0.38, which is not as large as expected. The probable explanation is that f nmin was determined by many factors that cannot be determined only by the zonal differential of VTEC. A better method should be applied to explore the inner mechanism for the correlation between f nmin and the zonal differential of VTEC.
For the statistical distribution of optimal f nmin, the ∆h = 20cm noise accuracy was set as the background. The statistical distribution of f nmin resembles the Gaussian distribution. Though the GNSS ranging error is timely correlated mainly due to atmospheric effects and multipath, as denoted by [22], it is still reasonable to make Gaussian noise assumption for the ionospheric delay in space-borne GNSS-R condition, since in this work, the observations were considered in 2 seconds, and the estimation of f nmin depends only on the number of observations rather than the time interval between observations. For the real space-borne GNSS-R observations, the effects of temporal correlation for the observations should be carefully studied. Under Gaussian noise assumption of ionospheric delay, 80% of the values were less than 55, corresponding to a distance along the orbit track of 1100 km. In other words, the ionospheric delay was filtered with good effect on a 1000-km spatial scale. The peak value of f nmin was obtained when f nmin was 29, corresponding to a distance along the orbit trace of 600 km. This means that under experimental assumptions, the optimal spatial filter statistical results for ionospheric delay lie in a distance range from 400 to 1000 km, which is determined by the spatial distribution of the ionosphere. Varied noise backgrounds were tested for statistical analysis to exclude occasionality.
In traditional altimetry technology, such as Jason series satellites and other altimetry satellites, the altimetry accuracy is relatively high at the centimeter level. For such high-precision altimetry, ionospheric delay correction was implemented, but there were still observation errors in the dual frequency calibration of ionospheric delay. To further improve the accuracy, the observation error needs to be eliminated, and smooth filtering is generally applied. The proposed spatial filtering distances of the Jason satellite were 100~150 km and 150~200 km, corresponding to 6:00-24:00 and 0-6:00 local time, respectively; that is, the spatial filtering reference distance of ionosphere delay in different local times is given. At present, there are no ionospheric delay observation data in space-borne GNSS-R satellites. Therefore, the simulation method was used to analyse the spatial filtering of ionospheric delay for space-borne GNSS-R altimetry. These experimental results can provide a useful reference for the spatial filtering method of ionospheric delay observation errors in GNSS-R altimetry and set the basis for follow-up research in related fields. More accuracy analysis will be conducted for the distribution of optimal spatial filtering points in different local times.

Conclusions
In this paper, a spatial filtering method of ionospheric delay observations in space-borne GNSS-R altimetry was studied. The ionospheric delay observations were computed with the GIM model provided by IGS. The spatial filtering analysis of these ionospheric delay observations was conducted with the absolute bias between the estimated value and the true value as the criterion to evaluate the filtering quality. The curves of absolute bias in relation to the bilateral filtering points f n were verified to achieve the minimum absolute bias. The global distribution of f nmin was investigated and compared with the ionospheric correction derived from the GIM, and it showed that f nmin and ionospheric correction were negatively correlated. The optimal filtering points f nmin were evaluated with two methods. One method used a statistical analysis to find the maximum probability density at f nmin = 20, 28, 32, corresponding to an altimetry accuracy of 10, 20, and 30 cm. The other method selected filtering points with a minimum value of the 99% quantile of σ( f n)/σ 0 to obtain a more stable filtering effect at f n = 10 corresponding to an altimetry accuracy of 10, 20, and 30 cm, with less computation assumption when f n = 10.