SLR troposphere delay modeling method and system based on ray tracing

Abstract

The invention discloses an SLR troposphere delay modeling method based on ray tracing, and the method comprises the steps: obtaining basic meteorological parameters which are defined on a multi-layer air pressure surface, and have a corresponding value at each horizontal grid node; respectively carrying out vertical interpolation, high altitude continuation and refractive index calculation based on the obtained basic meteorological parameters to obtain a refractive index field defined on the three-dimensional grid; and calculating integral delay along the observation path based on the refractive index field in combination with an SLR ray tracing algorithm. Specifically, a numerical weather model (NWM) is combined with a ray tracing technology to carry out three-dimensional modeling on an atmospheric refractive index field so as to realize delay correction with higher precision.

Description

Technical Field

[0001] The present invention specifically relates to a SLR tropospheric delay modeling method and system based on ray tracing. Background Art

[0002] Satellite laser ranging (SLR) is one of the four primary measurement techniques used to maintain the International Terrestrial Reference Frame (ITRF). Thanks to its high-precision measurements from ground stations to satellites, SLR has been widely used for satellite orbit verification and estimation of Earth's center of mass and scale parameters. To establish standards and ensure data consistency, the International Laser Ranging Service (ILRS) was established to promote research and applications related to SLR technology.

[0003] For decades, SLRs have been used to observe a growing number of satellites, including GNSS satellites, spherical satellites, and low-Earth orbit (LEO) satellites. Historically, SLR observations primarily targeted spherical satellites like LAGEOS, which are equipped only with retroreflectors and no other sensors, relying on Earth's gravity for orbital stability. Despite decades in orbit, LAGEOS satellites remain operational, providing critical center-of-mass support. Today, SLR observations are being used on a growing range of satellite types, such as GNSS satellites equipped with planar laser retroreflector arrays (LRAs). Within the GPS system, only GPS35 and GPS36 have SLR data, both of which have been retired. Most satellites in the GLONASS, BeiDou, and Galileo systems have SLR observations and are in good operating condition, with SLRs frequently used for orbit verification. ILRS has also deployed several observation missions for LEO satellites, which are typically equipped with hemispherical LRA arrays. For example, there are the GRACE and Swarm satellites equipped with tetrahedral LRAs, the Sentinel-3 satellite equipped with 7 cube reflectors arranged in a spherical pattern, and the Sentinel-6A and Jason satellites equipped with 9 cube reflectors. Due to the different orbital altitudes of different types of satellites, their tropospheric delays exhibit different characteristics.

[0004] Tropospheric delay is a significant source of error in SLR observations because it interferes with ranging bias estimation and affects the resolution of SLR data. The commonly used tropospheric correction model for SLR is based on the Mendes model. This model uses weather information from the observing station to calculate the zenithal tropospheric delay and then uses a mapping function to project this delay onto the corresponding elevation angle. This algorithm performs well at higher elevation angles, but at lower elevation angles, the mapping function struggles to accurately reflect the true tropospheric delay, potentially leading to errors in the low-elevation-angle estimates and compromising the quality of SLR observations. Summary of the Invention

[0005] To overcome the above-mentioned shortcomings of the prior art, the present invention provides a SLR tropospheric delay modeling method based on ray tracing. By adopting a numerical weather model (NWM) combined with ray tracing technology, the atmospheric refractive index field is three-dimensionally modeled to achieve higher-precision delay correction.

[0006] According to one aspect of the present invention, a ray tracing-based SLR tropospheric delay modeling method is provided, comprising:

[0007] Obtaining basic meteorological parameters, wherein the basic meteorological parameters are defined on a multi-layer pressure surface and have corresponding values ​​at each horizontal grid node;

[0008] Based on the basic meteorological parameters obtained, vertical interpolation, high-altitude extension and refractive index calculation are performed to obtain the refractive index field defined on the three-dimensional grid;

[0009] Based on the refractive index field, combined with the SLR ray tracing algorithm, the integrated delay along the observation path is calculated.

[0010] As a further technical solution, after calculating the integrated delay along the observation path, the tropospheric delay is calculated according to the following formula: , Where L represents the path length of the ray tracing algorithm integration, G represents the geometric path length of the electromagnetic wave propagating in a vacuum, represents the radial distance, represents the co-latitude, represents longitude, and t represents travel time.

[0011] As a further technical solution, the basic meteorological parameters include atmospheric pressure P, temperature T and water vapor pressure e.

[0012] As a further technical solution, vertical interpolation calculation is performed based on the acquired basic meteorological parameters, including: , in , and Represents the height values ​​of the lower and upper layers respectively. and Represent the meteorological parameter values ​​of the lower and upper layers respectively.

[0013] As a further technical solution, for areas above the highest layer of ECMWF, an exponential model is used for extrapolation, and the high-altitude extension expression of the pressure is: in Reference height The air pressure at is the atmospheric scale height.

[0014] As a further technical solution, refractive index calculation is performed based on the acquired basic meteorological parameters, including: , in, Represents the dry refractive index, calculated based on dry air; Represents the wet refractive index, calculated based on water vapor.

[0015] As a further technical solution, the method further includes:

[0016] The numerical meteorological model ERA5 is used to provide basic meteorological data.

[0017] According to one aspect of the present invention, a ray tracing-based SLR tropospheric delay modeling system is provided, comprising:

[0018] The first main module is used to obtain basic meteorological parameters, which are defined on a multi-layer pressure surface and have corresponding values ​​at each horizontal grid node;

[0019] The second main module is used to perform vertical interpolation, high-altitude extension and refractive index calculation based on the basic meteorological parameters obtained, and obtain the refractive index field defined on the three-dimensional grid;

[0020] The third main module is used to calculate the integrated delay along the observation path based on the refractive index field in combination with the SLR ray tracing algorithm.

[0021] As a further technical solution, it also includes:

[0022] The fourth main module is used to calculate the tropospheric delay according to the following formula: , Where L represents the path length of the ray tracing algorithm integration, G represents the geometric path length of the electromagnetic wave propagating in a vacuum, represents the radial distance, represents the co-latitude, represents longitude, and t represents travel time.

[0023] According to one aspect of the present invention, a non-transitory computer-readable storage medium is provided, wherein the non-transitory computer-readable storage medium stores computer instructions, wherein the computer instructions enable the computer to execute the SLR tropospheric delay modeling method based on ray tracing.

[0024] Compared with the prior art, the present invention has the following beneficial effects:

[0025] The method of constructing a three-dimensional refractive index field for SLR based on ECMWF meteorological data covers key links such as meteorological data interpolation, standard atmosphere extension, and refractive index sub-item modeling. Compared with traditional mapping functions, it has higher spatiotemporal adaptability and accuracy, especially at low altitude angles and when the weather changes drastically. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, a brief introduction will be given below to the drawings used in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0027] Figure 1 A schematic flow chart of a ray tracing-based SLR tropospheric delay modeling method provided in an embodiment of the present invention.

[0028] Figure 2 Schematic diagrams of three-dimensional (b) and two-dimensional (a) ray tracing algorithms provided by embodiments of the present invention.

[0029] Figure 3 This is a distribution map of measurement stations provided by an embodiment of the present invention.

[0030] Figure 4 Schematic diagram comparing the absolute values ​​of SLR residuals between the MF group and the Ray-3d-1h group provided by an embodiment of the present invention.

[0031] Figure 5 This is a boxplot of the SLR residuals provided by an embodiment of the present invention.

[0032] Figure 6 The embodiments of the present invention provide a tropospheric delay zenith map and a difference zenith map.

[0033] Figure 7 (a) MF group SLR residual zenith map (b) Ray-3d-1h group SLR residual zenith map (c) Ray-3d-6h group SLR residual zenith map (d) Ray-2d-1h group SLR residual zenith map provided by the embodiment of the present invention.

[0034] Figure 8(a) MF group SLR residual zenith map (b) Ray-3d-1h group SLR residual zenith map (c) Ray-3d-6h group SLR residual zenith map (d) Ray-2d-1h group SLR residual zenith map provided by the embodiment of the present invention.

[0035] Figure 9 (a) MF group SLR residual zenith map (b) Ray-3d-1h group SLR residual zenith map (c) Ray-3d-6h group SLR residual zenith map (d) Ray-2d-1h group SLR residual zenith map provided by the embodiment of the present invention.

[0036] Figure 10 A schematic diagram of SLR residual statistics for elevation angles below 20 degrees provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0037] Ray tracing, based on fundamental principles of geometric optics, calculates tropospheric delay by determining the signal path from a space transmitter to a ground receiver. Ray tracing is considered the most accurate method for measuring tropospheric delay and has been widely used in VLBI data processing. For VLBI tropospheric delay estimation, 3D ray tracing algorithms offer the highest accuracy but also require significant computational resources. 2D ray tracing methods significantly improve computational efficiency by restricting the algorithm to the vertical plane through a fixed azimuth angle, but their effectiveness in SLR observations requires further verification.

[0038] ERA5 is the latest numerical weather model (NWM) product from the European Centre for Medium-Range Weather Forecasts (ECMWF), replacing ERA-Interim since 2019. It offers higher spatial resolution (approximately 31 kilometers) and finer temporal resolution (1 hour). However, the 1-hour resolution NWM requires significant computer storage resources, leading many studies to interpolate 6-hour weather products for tropospheric delay processing. Using a 6-hour resolution NWM may not accurately capture rapid changes in water vapor activity, potentially causing problems in tropospheric delay calculations. Whether a 6-hour resolution NWM can be effectively used for SLR observations also requires further study. To further clarify the objectives, technical solutions, and advantages of the embodiments of the present invention, the technical solutions of the embodiments will be described below in detail, in conjunction with the accompanying drawings. It should be understood that the described embodiments represent only a portion of the embodiments, not all of them. All other embodiments derived by persons of ordinary skill in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention. In addition, the technical features in the various embodiments or single embodiments provided by the present invention may be arbitrarily combined with each other to form a new technical solution. Such combination is not restricted by the sequence of steps and / or structural composition pattern, but must be based on the ability of ordinary technicians in this field to implement it. When the combination of technical solutions is mutually contradictory or cannot be implemented, it should be deemed that such a combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0039] The embodiment of the present invention provides a SLR tropospheric delay modeling method based on ray tracing, such as Figure 1 As shown, first, basic meteorological parameters are obtained. The basic meteorological parameters are defined on a multi-layer pressure surface and have corresponding values ​​at each horizontal grid node. Then, based on the obtained basic meteorological parameters, vertical interpolation, high-altitude extension and refractive index calculation are performed to obtain a refractive index field defined on a three-dimensional grid. Subsequently, based on the refractive index field, the integral delay along the observation path is calculated in combination with the SLR ray tracing algorithm.

[0040] Satellite laser ranging (SLR) determines distance by measuring the time of flight of a laser. Laser is a specific type of electromagnetic wave, and its observation equation is as follows:

[0041]

[0042] in, is the SLR observation value; c is the speed of light; It represents the round-trip propagation time of the laser between the station and the satellite; is the geometric distance between the station and the satellite; the systematic error in SLR observations—the ranging bias—is given by The time deviation caused by the time tag error is represented by The atmospheric delay caused by tropospheric refraction is represented by , while the delay caused by relativistic effects is given by The movement of the station due to factors such as solid tide and tidal load is represented by express; Indicates noise.

[0043] like Figure 2 As shown in (a) and (b), ray tracing algorithms can be divided into three-dimensional (3D) and two-dimensional (2D) algorithms, depending on the method used to approximate the ray propagation path. In 3D algorithms, rays propagate in three-dimensional space, affected by both horizontal and vertical refractive index gradients. In contrast, 2D algorithms restrict rays to a specific azimuthal plane, where they are affected only by the vertical refractive index gradient.

[0044] The three-dimensional ray tracing algorithm is considered the most accurate method for integrating the actual propagation path of electromagnetic waves. The tropospheric delay calculated using this method can be described as the difference between the path length L integrated by the ray tracing algorithm and the geometric path length G of the electromagnetic wave propagating in a vacuum, expressed by the following formula:

[0045]

[0046]

[0047] Where L represents the path length of the ray tracing algorithm integration, G represents the geometric path length of the electromagnetic wave propagating in a vacuum, represents the radial distance, represents the colatitude (ranging from 0 to π), represents longitude (ranging from 0 to 2π), and t represents the travel time.

[0048] The above 3D ray tracing algorithm determines the path as a three-dimensional space curve, which is very close to the actual path of electromagnetic waves propagating in a neutral atmosphere. By restricting the ray path to a vertical plane with a fixed azimuth angle, the three-dimensional ray tracing system can be simplified to a two-dimensional system. In this case, the refractive index at the colatitude is and longitude The gradient component in the direction is set to zero.

[0049]

[0050] This method assumes that when an electromagnetic wave passes through a neutral atmosphere, the refractive index path has only radial (vertical) gradients and lacks horizontal gradients. Existing research shows that three-dimensional ray tracing algorithms have the highest accuracy, while two-dimensional algorithms are more computationally efficient.

[0051] Preferably, the embodiment of the present invention uses observation data from 26 ILRS stations, such as Figure 3 These stations are primarily located in Oceania, Europe, East Asia, and North America. The observation period covers 060 to 151 DOY (Day of the Year) in 2023. The data used is in the conventional point format. It is worth noting that the ILRS announced on January 25, 2023, that it would no longer accept V1 CRD and CPF format files. Therefore, all observation data used in the examples of this invention are derived from V2 np2 files, rather than the previously common npt files.

[0052] Figure 3 In the station distribution diagram, the size of the circles represents the number of observations, with the largest circles representing the YARL stations with the most observations. This example uses SLR observation data from multiple types of satellites, including spherical satellites, GNSS satellites, and low-Earth orbit (LEO) satellites, as shown in Table 1. Each type of satellite has different characteristics.

[0053] Table 1 Satellite information table

[0054] Satellite type name Altitude (km) LRA Type Track source Spherical satellite Lageos1 5850 Covered with 426 corner reflectors ILRSB Spherical satellite Lageos2 5625 Covered with 426 corner reflectors ILRSB LEO GRACE-C 500 Pyramid LRA, including 4 corner reflectors NASA / JPL LEO GRACE-D 500 Pyramid LRA, including 4 corner reflectors NASA / JPL LEO JASON3 1336 Hemispherical LRA, including 9 corner reflectors CNES LEO Sentinel-6A 1336 Hemispherical LRA, including 9 corner reflectors ESA / CPOD GNSS BDS 20000-36000 Planar array composed of multiple corner reflectors WHU GNSS GALILEO 23200-23650 Planar array composed of multiple corner reflectors WHU

[0055] GNSS satellites typically orbit at high altitudes, ranging from 20,000 to 40,000 kilometers. Their laser retroreflector arrays (LRAs) typically consist of a planar array of multiple corner reflectors. This LRA ensures the feasibility of GNSS satellite SLR observations. However, the difficulty in determining the laser reflection position within the LRA can introduce systematic errors, making ranging corrections challenging and affecting observation results.

[0056] Spherical satellites typically operate at altitudes of several thousand kilometers, but there are exceptions, such as the Etalon satellite, which is close to 20,000 kilometers (Appleby 1998). These spherical satellites are typically covered with a large number of corner reflectors, and the SLR observations can be corrected to the satellite's optical center using the products published by (Rodríguez, Appleby, and Otsubo 2019), effectively reducing most of the systematic errors associated with ranging corrections.

[0057] LEO satellites orbit at altitudes ranging from a few hundred kilometers to over a thousand kilometers. They are equipped with LRAs of various shapes, typically consisting of several corner reflectors mounted on a hemispherical mount. Generally speaking, the lower the satellite's altitude, the smaller its LRA. For example, the GRACE satellites are equipped with pyramidal LRAs consisting of four corner reflectors. In contrast, the JASON-3 and Sentinel-6A satellites are equipped with hemispherical LRAs consisting of nine corner reflectors. Due to the higher orbital altitude of the JASON satellites, their LRAs are significantly larger than those of the GRACE satellites. Ranging corrections for these LEO satellites can be adjusted using the azimuth and elevation models proposed by (Montenbruck and Neubert 2011).

[0058] The data processing strategy used in the embodiment of the present invention is shown in Table 2. SLRF2020 is an extension of ITRF2020, which includes additional SLR stations not covered by ITRF and provides high-quality coordinates for SLR stations. SLR observations with residual values ​​greater than 0.15m are considered abnormal and are discarded. The LRA offset is a vector pointing from the satellite's centroid to the LRA reference point and is used for correction. Because the Galileo and BeiDou systems are still under development and their satellite information is updated frequently, it may be difficult for ILRS to collect statistics in a timely manner. Therefore, the LRA offsets of Galileo and BeiDou satellites are derived from their official metadata websites. For LEO satellites, ILRS has well collected LRA offsets, which can be obtained directly from their websites.

[0059] Table 2 Data processing strategy table

[0060] project Description or source Reference Frame SLRF2020 (Pavlis et al. 2021) Observational data Normal points Residual elimination threshold 0.15m Galileo LRA offset Galileo Metadata Website BDS LRA offset BDServices website LEO LRA offset ILRS Ray tracing software RADIATE (Hofmeister 2016) Ray tracing method 3D Numerical Weather Model (NWM) ERA5 NWM horizontal resolution 0.25° x 0.25° (approximately 31 km) NWM temporal resolution One hour Refractive index calculation Based on dry air and water vapor contributions (Mendes and Pavlis 2004) Vertical interpolation method Linear interpolation for pressure, temperature, and water vapor; logarithmic interpolation for water vapor where applicable Height limit 84000m Supplementary atmospheric models 1976 U.S. Standard Atmosphere

[0061] The ray tracing software used in the present invention is based on RADIATE, an open-source tool developed by the Vienna University of Technology (TU Wien) (Hofmeister and Böhm 2017). The software supports ray tracing for various spatial geodetic techniques, including GNSS, VLBI, and SLR, and can use both 2D and 3D ray tracing algorithms.

[0062] NWM uses ERA5, which has a spatial resolution of 0.25° x 0.25° (approximately 31 kilometers) and a temporal resolution of one hour (Hersbach et al. 2020). In ray tracing, ERA5's high spatial and temporal resolution provides high-quality meteorological data, significantly improving the accuracy of the results (Zhang et al. 2023). The refractive index algorithm used in this example follows the method proposed by (Mendes and Pavlis 2004), splitting the calculation into dry air and water vapor. The NWM has an upper altitude limit of 84,000 meters. For altitudes above this limit, the standard atmosphere model provided by the 1976 U.S. Standard Atmosphere is applied.

[0063] The results are divided into four groups, as shown in Table 3. The first group uses the most classic mapping function method provided by (Mendes and Pavlis 2004), designated the Mapping Function Group (MF). The second group employs a 3D ray tracing algorithm with a 1-hour ERA5 product, designated the Ray-3D-1h Group. The third group employs a 3D ray tracing algorithm with a numerical weather model (NWM) interpolated from a 6-hour ERA5 product, designated the Ray-3D-6h Group. The fourth group employs a 2D ray tracing algorithm with a 1-hour ERA5 product. The results from these groups are compared to verify the effectiveness of the ray tracing algorithm in correcting tropospheric delay estimates from SLR observations of different satellite types. This comparison also explores the impact of the 2D and 3D algorithms on the results and assesses the differences caused by different meteorological product time intervals.

[0064] Table 3 Results grouping

[0065] Comparison of the absolute values ​​of SLR residuals between the MF group and the Ray-3d-1h group Figure 4 As shown in Figure 2, red scatter points indicate that the absolute value of the SLR residuals is larger in the MF group, while blue scatter points indicate that the SLR residuals are larger in the Ray-3d-1h group. The larger the difference, the further the point deviates from the middle auxiliary line.

[0066] The present invention compares the SLR residuals of the MF group and the Ray-3d-1h group, using 0 as the reference value, as shown in Figure 4This figure shows the difference between the two sets of SLR residuals, using their absolute values. The red scattered points in the figure represent instances where the absolute SLR residuals of the MF group exceed those of the Ray-3d-1h group. Clearly, for most satellites, the ray tracing method reduces the absolute values ​​of the SLR residuals, demonstrating its effectiveness in minimizing the impact of tropospheric delay on SLR observations and reducing their systematic errors.

[0067] Comparing the performance of different satellites reveals that the ray tracing method is particularly effective for Low Earth Orbit (LEO) satellites. For example, for the GRACE-C satellite, the Ray-3d-1h group had 2,741 observations with large residual errors, while the MF group had 3,814 observations, a reduction of approximately 28%. Similarly, the Ray-3d-1h group reduced residual errors for the GRACE-D satellite by 22%, for the Jason-3 satellite by 11%, and for the Sentinel-6A satellite by 14%. These results demonstrate that the ray tracing method effectively improves observational results when processing SLR data from LEO satellites.

[0068] like Figure 5 The SLR residual boxplot shown uses the interquartile range (IQR) as an evaluation metric to avoid outliers that could lead to gross errors. Research has found that Global Navigation Satellite System (GNSS) satellites have large SLR residuals, which is attributed to the low orbital accuracy and the low accuracy of the Laser Retroreflector Array (LRA) range corrections used for GNSS satellites. Because GNSS satellites consist of a planar array of multiple corner reflectors, it is difficult to determine the precise location of the reflected light based on the incident azimuth and elevation angles of the laser, making correction difficult and resulting in large SLR residuals.

[0069] Comparing the MF group with the other three groups, little difference in the SLR residuals was observed for high-orbit satellites such as GNSS and LAGEOS. However, for LEO satellites, the use of ray tracing significantly reduced the IQR of the SLR residuals. For example, for the GRACE-C satellite, the IQR for the MF group was 18.6 mm, for the Ray-3d-1h group was 16.2 mm, for the Ray-3d-6h group was 16.7 mm, and for the Ray-2d-1h group was 17.1 mm, representing reductions of 13%, 10%, and 8%, respectively, after ray tracing. For the GRACE-D satellite, the IQR for the MF group was 15.7 mm, for the Ray-3d-1h group was 13.4 mm, for the Ray-3d-6h group was 14.5 mm, and for the Ray-2d-1h group was 14.6 mm, representing reductions of 15%, 8%, and 7%, respectively. In addition, the IQR of the Jason3 satellite is reduced by 2-4% after applying ray tracing, and the IQR of the Sentinel-6a satellite is reduced by 4-8%.

[0070] Overall, the improvement brought by ray tracing is mainly reflected in LEO satellites, and the IQR difference of medium earth orbit (MEO) satellites is very small, which is consistent with the Figure 4 The results are consistent with those shown. However, the ray tracing method itself does not differentiate between MEO and LEO satellites for adjustment, indicating that SLR observations of LEO satellites exhibit distinct characteristics compared to those of MEO satellites. Furthermore, among the three ray tracing results, the Ray-3d-1h group generally performs best for most LEO satellites, which is consistent with our expectations.

[0071] To determine why the ray tracing method performs differently for LEO and MEO satellites, this example compares the tropospheric delay values ​​of the MF and Ray-3d-1h groups. For example, the YARL station, located in western Oceania, has the most SLR observation data. It is the most commonly used SLR station and its observation data is relatively high quality.

[0072] Table 4 Statistics of the number of SLR observations and the proportion of low-altitude angle observations

[0073] Since satellites with similar orbital altitudes have similar performance in SLR observations, four satellites were selected for analysis: BeiDou (BDS) satellites with orbital altitudes exceeding 20,000 km; LAGEOS-1 satellites with orbital altitudes of several thousand kilometers; Sentinel-6A satellites with altitudes above 1,000 km but below 2,000 km; and GRACE-C satellites with altitudes below 1,000 km. Figure 6 In each sub-figure, the left figure shows the total delay (STD) of the slant path obtained by the mapping function group (MF), in meters; the right figure shows the difference in STD between the ray tracing group and the mapping function group, in millimeters.

[0074] As can be seen from the comparison in the left figure, the STD increases with decreasing altitude. This is because at lower altitudes, the laser travels a longer path through the troposphere and encounters greater variations in the refractive index, resulting in increased tropospheric delay. The distribution of the SLR residuals also differs due to differences in satellite orbital inclination and altitude. For BDS satellites, their altitude and inclination characteristics result in fewer low-altitude observations, while LEO satellites like GRACE-C and Sentinel-6A have relatively more low-altitude observations.

[0075] Comparing the right figure reveals that the differences between the ray tracing algorithm and the mapping function algorithm primarily occur at low altitude angles. At higher altitude angles, the STDs obtained by the two methods are nearly identical. However, when the altitude angle drops below 30°, a clear difference between the MF and ray tracing groups begins to emerge, becoming significant below 20°. As shown in Table 4, since the proportion of low-altitude SLR observations for LEO satellites is much higher than for MEO satellites, and the differences between the ray tracing and mapping function methods are concentrated at low altitude angles, the ray tracing algorithm demonstrated superior performance for LEO satellites in the previous statistical analysis.

[0076] Table 5. SLR residual statistics of the LAGEOS-1 satellite at the YARL station (unit: mm) (low-elev refers to statistics with elevation angles below 20 degrees)

[0077] like Figure 7Figures (a) and (b) show the results for the MF group, while (c) and (d) show the results for the ray tracing group. For the LAGEOS-1 satellite, due to the limited number of low-altitude observations, the overall differences between the groups were not significant. As shown in Table [tab-lageos1-7090], the IQR differences among the four groups were small. However, after applying the ray tracing method, the overall mean was closer to zero. For example, the mean for the Ray-3d-1h group was 2.54 mm, a decrease of 1.72 mm from the MF group's 4.26 mm.

[0078] In previous studies, we found that ray tracing methods primarily improve low-altitude SLR observations. A comparison showed that the Ray-3d-1h group achieved the best results, with an IQR of 16.98 mm at low altitudes, a 13.5% reduction compared to the MF group's 19.62 mm. The Ray-2d-1h and Ray-3d-6h groups failed to achieve the same level of correction as the Ray-3d-1h method, possibly due to lower algorithm accuracy or the coarser resolution of the NWM. Nevertheless, they still provide some improvement in tropospheric delay estimates for low-altitude observations from the LAGEOS-1 satellite.

[0079] Table 6. SLR residual statistics of Sentinel-6A satellite at YARL station (unit: mm) (low-elev refers to statistics with elevation angles below 20 degrees)

[0080] like Figure 8 As shown, the significantly deviated SLR residuals in the MF group are concentrated in the low-elevation angle regions in the northwest and northeast regions. Differences between the other three groups and the MF group are also concentrated in these observations. The key difference is that in the Ray-3d-1h group, the SLR residuals in these regions approach zero after ray tracing. However, in the Ray-3d-6h group, the SLR residuals decrease excessively after ray tracing, resulting in significantly negative values, indicating that untimely NWM updates can lead to poor ray tracing results.

[0081] The results for the Ray-2d-1h group appear more unstable than those for the Ray-3d-1h group. The low-elevation angle SLR residuals become more negative in the northeast, while they become more positive in the southeast, both moving further away from zero. This suggests that insufficient algorithmic precision may also introduce errors.

[0082] As shown in Table 6, the IQR for the low-elevation angle was 8.22 mm for the MF group, 6.81 mm for the Ray-3d-1h group, 8.58 mm for the Ray-3d-6h group, and 8.37 mm for the Ray-2d-1h group. The Ray-3d-1h group performed best, achieving a 17% reduction compared to the MF group. In contrast, the Ray-3d-6h and Ray-2d-1h groups actually showed an increase relative to the MF group, highlighting the importance of accurate NWM and ray tracing algorithms.

[0083] Table 7. SLR residual statistics of GRACEC satellite at YARL station (unit: mm) (low-elev refers to the statistics with elevation angles below 20 degrees)

[0084] like Figure 9 As shown, the GRACEC satellite is a polar-orbiting satellite with an orbital inclination of 89 degrees (close to 90 degrees). Therefore, its residuals appear as stripes in the zenith diagram. Due to its relatively low orbital altitude, a large number of low-altitude SLR observations exist, accounting for approximately 18% of the total observations, with elevations below 20 degrees.

[0085] Comparing the MF group with the other three groups reveals similar results at higher elevation angles. However, after applying the ray tracing algorithm, low-elevation angle observations showed significant improvement. Because low-elevation angle observations account for a large proportion of the total observations, their impact on the overall statistical results is more significant. As shown in Table 7, the overall IQR for the MF group was 12.18 mm, compared with 11.16 mm for the Ray-3d-1h group, 11.60 mm for the Ray-3d-6h group, and 11.23 mm for the Ray-2d-1h group, representing reductions of 8.4%, 4.8%, and 7.8%, respectively.

[0086] Due to the presence of numerous low-angle observations distributed across multiple regions, the Ray-3d-1h group's results show significant improvement compared to the MF group in the west, north, east, southwest, and southeast directions. The results of Ray-3d-6h and Ray-2d-1h are generally similar to those of Ray-3d-1h, with some exceptions. For example, in the southwest, the low-angle observations of Ray-3d-6h and Ray-2d-1h show not only no improvement but actually deteriorate, with larger SLR residuals after ray tracing. These issues are likely due to insufficient precision in the NWM and algorithm.

[0087] Overall, the IQR for the low-elevation angle was 14.30 mm for the MF group, 10.18 mm for the Ray-3d-1h group, 11.59 mm for the Ray-3d-6h group, and 11.34 mm for the Ray-2d-1h group. The Ray-3d-1h group performed best, achieving a 29% reduction compared to the MF group. While the Ray-3d-6h and Ray-2d-1h groups also showed improvement relative to the MF group, anomalies caused their overall performance to fall short of that of the Ray-3d-1h group.

[0088] Figure 10 The SLR residual statistics for elevation angles below 20 degrees are given. Figure 10 As shown in Figure 2, this example lists the SLR residuals for all elevation angles below 20 degrees and uses the IQR to verify the improvement achieved by the ray tracing method. The SLR residuals for the Lageos1 satellite are relatively large because its orbital accuracy is weaker than that of other satellites, and its higher orbital altitude results in fewer low-elevation angle observations, increasing randomness. For the Lageos1 satellite, the IQR for the MF group is 19.5 mm; for the Ray-3d-1h group it is 17.9 mm; for the Ray-3d-6h group it is 17.7 mm; and for the Ray-2d-1h group it is 17.4 mm, indicating that the IQR is reduced by approximately 2 mm after using ray tracing.

[0089] For both GRACEC and GRACED satellites, the Ray-3d-1h group showed significant improvements. The IQR for GRACEC decreased by approximately 3.8 mm (approximately 27%), while the IQR for GRACED decreased by 2.4 mm (approximately 21%). However, the improvements for the Ray-3d-6h and Ray-2d-1h groups were less significant, likely due to lower algorithmic and NWM accuracy.

[0090] For Jason3 and Sentinel-6a, the ray tracing algorithm did not show the same significant improvement as for GRACE. The Ray-3d-6h and Ray-2d-1h groups showed significantly worse results, indicating that 2D ray tracing and 6-hour interval NWM products still have issues in SLR data processing. It is recommended that 1-hour interval NWM and 3D ray tracing be used preferentially in SLR observation processing.

[0091] It is worth noting that the number of low-elevation angle observations from Lageos2, BDS, and Galileo are 48, 33, and 27, respectively. These satellites have very few low-elevation angle observations, which may lead to greater randomness and are therefore not shown in the figure.

[0092] This example evaluates the performance of tropospheric delay correction using classical mapping functions and various ray tracing methods in SLR observations from three perspectives: statistical performance comparison, physical explanation of the advantage of ray tracing at low altitude angles, and the impact of dimensionality and meteorological model resolution on ray tracing accuracy, with a focus on the differences in satellite altitudes and observation geometries.

[0093] First, statistical analysis shows that ray tracing—particularly the 3D method using 1-hour ERA5 data (Ray-3d-1h)—significantly improves tropospheric delay corrections at low altitudes. For Low Earth Orbit (LEO) satellites like GRACE-C and GRACE-D (which have a high proportion of low-altitude observations, approximately 18% of the total observations), the overall SLR residual IQR is reduced by 13–15%, while the IQR for low-altitude observations alone is reduced by 21–27%. For Medium Earth Orbit (MEO) satellites like LAGEOS-1, applying the ray tracing algorithm reduces the IQR for low-altitude observations by approximately 8%. However, since low-altitude observations only account for 2% of the total LAGEOS-1 observations, their impact on the overall SLR residual is limited. In summary, ray tracing significantly improves tropospheric delay estimates for low-altitude SLR observations.

[0094] Secondly, the superiority of the ray tracing algorithm over the classical mapping function stems primarily from its physical modeling of the laser pulse propagation path in a stratified, refractive atmosphere. The mapping function estimates the STD by empirically projecting the zenith total delay (ZTD), while ray tracing obtains the STD by directly integrating the three-dimensional refractive index field. Benefiting from the high-precision NWM, ray tracing can more accurately capture the true structure of tropospheric delays. This distinction becomes particularly critical at low altitudes, where the atmospheric path is longer and the refractive index gradient is more complex. SLR residual analysis shows that the difference between the mapping function and ray tracing begins to increase at altitudes below 30° and is most significant below 20°. Because LEO satellites are frequently involved in such low-altitude observations, the benefits of ray tracing are significantly amplified.

[0095] Third, comparisons between the three ray-tracing configurations demonstrate that algorithm accuracy and meteorological input resolution are equally critical. 3D ray tracing using 1-hour interval ERA5 data (Ray-3d-1h) outperforms both the 2D method (Ray-2d-1h) and the 3D configuration using 6-hour interval interpolated data (Ray-3d-6h) in most cases. The 2D algorithm simplifies the signal path to a flat surface, failing to capture atmospheric asymmetries—particularly under azimuth variations or complex weather conditions. This simplification can introduce local errors and lead to unstable residual corrections. Similarly, using 6-hour interval interpolated data reduces the temporal resolution of the refractivity field, potentially missing rapid atmospheric variations. In some cases, Ray-3d-6h and Ray-2d-1h not only fail to improve the results but can even worsen them, particularly in low-altitude observations from Sentinel-6A and GRACE-C. These findings clearly demonstrate that the spatial dimensionality and temporal resolution of the atmospheric model are equally important for achieving accurate ray-tracing corrections.

[0096] Table 8 Statistics of time and storage resources required for tropospheric delay correction using different ray tracing methods (daily average)

[0097] The time and storage required for tropospheric delay correction using ray tracing are summarized in Table [tab-resource]. Ray tracing is performed hourly, so the total processing time is calculated as 24 times the time required for each hour of data. Among these methods, Ray-3d-6h is the most time-consuming because it requires two NWM products per hour—one corresponding to the first and last hour of the six-hour window. For example, to process SLR observations between 03:00 and 04:00, interpolation is required using the NWMs from 00:00–01:00 and 06:00–07:00, both of which require the construction of a refractivity field. Due to the small number of SLR observations, the actual time spent on ray tracing is negligible; the majority of the time is spent constructing the refractivity field, which is not particularly time-consuming. Therefore, the processing time for Ray-3d-1h remains acceptable. Regarding storage, ERA5 NWMs are stored hourly, with each file approximately 1.3 GB. Therefore, Ray-3d-1h and Ray-2d-1h require approximately 1.3 GB × 24 per day, while Ray-3d-6h only requires NWM at 00:00, 06:00, 12:00, and 18:00. Although storage usage is a consideration, we still recommend using Ray-3d-1h because SLR processing results may occasionally deteriorate when using Ray-3d-6h.

[0098] In summary, a 3D ray tracing method using high-resolution (1-hour interval) meteorological data offers significant advantages over traditional mapping functions for SLR tropospheric delay correction, particularly for LEO satellites and low-altitude observations. The results also suggest that simplifying ray tracing—either by reducing dimensionality or using data with coarser temporal resolution—can negate these benefits and introduce new biases. Therefore, future implementations of tropospheric correction in SLR data processing should prioritize full 3D ray tracing using high-resolution, timely meteorological inputs.

[0099] The various embodiments of the present invention are implemented through programmed processing by a device with processor functionality. Therefore, in practical engineering, the technical solutions and functions of the various embodiments of the present invention are encapsulated into various modules. Based on this reality, and in addition to the aforementioned embodiments, embodiments of the present invention provide a ray tracing-based SLR tropospheric delay modeling system. This system is used to implement the ray tracing-based SLR tropospheric delay modeling method described in the aforementioned method embodiments.

[0100] The system includes: a first main module, used to obtain basic meteorological parameters, which are defined on a multi-layer pressure surface and have corresponding values ​​on each horizontal grid node; a second main module, used to perform vertical interpolation, high-altitude extension and refractive index calculation based on the obtained basic meteorological parameters to obtain a refractive index field defined on a three-dimensional grid; a third main module, used to calculate the integral delay along the observation path based on the refractive index field in combination with the SLR ray tracing algorithm.

[0101] An embodiment of the present invention provides a ray tracing-based SLR tropospheric delay modeling system, which adopts the aforementioned modules and uses a numerical weather model (NWM) combined with ray tracing technology to perform three-dimensional modeling of the atmospheric refractive index field to achieve higher-precision delay correction.

[0102] It should be noted that the system embodiments provided by the present invention are not only used to implement the methods in the above-mentioned method embodiments, but also used to implement the methods in other method embodiments provided by the present invention. The only difference lies in the setting of corresponding functional modules, and the principles thereof are basically the same as the principles of the above-mentioned system embodiments provided by the present invention. As long as those skilled in the art refer to the specific technical solutions in other method embodiments on the basis of the above-mentioned system embodiments, obtain corresponding technical means and technical solutions composed of these technical means by combining technical features, and on the premise of ensuring the practicality of the technical solutions, improve the modules in the above-mentioned system embodiments to obtain corresponding system class embodiments for implementing the methods in other method class embodiments.

[0103] Based on the same inventive concept as the aforementioned embodiment, an embodiment of the present invention further provides a non-transitory computer-readable storage medium, wherein the non-transitory computer-readable storage medium stores computer instructions, and the computer instructions enable the computer to execute the SLR tropospheric delay modeling method based on ray tracing.

[0104] In summary, the present invention uses high temporal and spatial resolution meteorological data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) as modeling input to construct a three-dimensional refractive index field that varies with space and altitude. The modeling process consists of three main stages: vertical interpolation of meteorological parameters, high-altitude extension, and refractivity calculation. An in-depth analysis of the constant factors and formula structure used in the Fortran code provides a refined modeling foundation for high-precision SLR delay correction.

[0105] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the technical solutions of the embodiments of the present invention.

Claims

1. A ray tracing-based SLR tropospheric delay modeling method, characterized in that: include: Obtaining basic meteorological parameters, wherein the basic meteorological parameters are defined on a multi-layer pressure surface and have corresponding values ​​at each horizontal grid node; Based on the basic meteorological parameters obtained, vertical interpolation, high-altitude extension and refractive index calculation are performed to obtain the refractive index field defined on the three-dimensional grid; Based on the refractive index field, combined with the SLR ray tracing algorithm, the integrated delay along the observation path is calculated.

2. The SLR tropospheric delay modeling method based on ray tracing according to claim 1, characterized in that: After calculating the integrated delay along the observation path, the tropospheric delay is calculated as follows: ， Where L represents the path length of the ray tracing algorithm integration, G represents the geometric path length of the electromagnetic wave propagating in a vacuum, represents the radial distance, represents the co-latitude, represents longitude, and t represents travel time.

3. The SLR tropospheric delay modeling method based on ray tracing according to claim 1, characterized in that: The basic meteorological parameters include atmospheric pressure P, temperature T and water vapor pressure e.

4. The SLR tropospheric delay modeling method based on ray tracing according to claim 3, characterized in that: Based on the basic meteorological parameters obtained, vertical interpolation calculations are performed, including: ， in , and Represents the height values ​​of the lower and upper layers respectively. and Represent the meteorological parameter values ​​of the lower and upper layers respectively.

5. The SLR tropospheric delay modeling method based on ray tracing according to claim 3, characterized in that: For the area above the highest layer of ECMWF, the exponential model is used for extrapolation, and the high-altitude extension expression of the pressure is: in Reference height The air pressure at is the atmospheric scale height.

6. The SLR tropospheric delay modeling method based on ray tracing according to claim 3, characterized in that: Based on the basic meteorological parameters obtained, refractive index calculation is performed, including: ， in, Represents the dry refractive index, calculated based on dry air; Represents the wet refractive index, calculated based on water vapor.

7. The SLR tropospheric delay modeling method based on ray tracing according to claim 1, characterized in that: The method further comprises: using the numerical meteorological model ERA5 to provide basic meteorological data.

8. A ray tracing based SLR tropospheric delay modeling system, characterized in that: include: The first main module is used to obtain basic meteorological parameters, which are defined on a multi-layer pressure surface and have corresponding values ​​at each horizontal grid node; The second main module is used to perform vertical interpolation, high-altitude extension and refractive index calculation based on the basic meteorological parameters obtained, and obtain the refractive index field defined on the three-dimensional grid; The third main module is used to calculate the integrated delay along the observation path based on the refractive index field in combination with the SLR ray tracing algorithm.

9. The SLR tropospheric delay modeling system based on ray tracing according to claim 8, characterized in that: Also includes: The fourth main module is used to calculate the tropospheric delay according to the following formula: ， Where L represents the path length of the ray tracing algorithm integration, G represents the geometric path length of the electromagnetic wave propagating in a vacuum, represents the radial distance, represents the co-latitude, represents longitude, and t represents travel time.

10. A non-transitory computer-readable storage medium, characterized in that The non-transitory computer-readable storage medium stores computer instructions, and the computer instructions enable the computer to execute the SLR tropospheric delay modeling method based on ray tracing according to any one of claims 1 to 7.

Images