Atmospheric delay correction method based on three-dimensional ray tracing

Abstract

The method for atmospheric delay correction based on three-dimensional light ray tracing includes obtaining the incidence angle of laser at the top of the atmosphere by using the position vector of the top of the atmosphere and the first layer laser unit incident vector; layering the atmosphere; layering the atmospheric refractive index of the actual trajectory, layering the actual trajectory of laser passing through the atmosphere until the laser reaches the ground, thereby obtaining the position coordinates of the laser spot on the ground; calculating the trajectory of laser passing through the atmosphere under the theoretical condition; using the obtained atmospheric refractive index and the actual trajectory of laser, integrating the atmospheric refractive index along the actual trajectory of laser and subtracting the theoretical distance, thereby obtaining the atmospheric delay. For the spaceborne laser radar system under the condition of large incidence angle, the relationship between the mapping function and the zenith delay is used to calculate the atmospheric delay, which may cause a larger error, the present application can accurately calculate the light ray trajectory and the atmospheric refractive index in the process of laser passing through the atmosphere, thereby obtaining higher calculation precision and providing high-precision atmospheric delay correction.

Description

Technical Field

[0001] This invention relates to an atmospheric delay correction method based on three-dimensional ray tracing, and more particularly to a ranging error correction method for atmospheric delay generated during Earth observation by spaceborne lidar under large incident angle conditions, belonging to the field of satellite remote sensing. Background Technology

[0002] Spaceborne marine lidar can acquire subsurface information of the ocean, enabling ocean remote sensing to move from two-dimensional to three-dimensional observation. It is not limited by sunlight and can achieve long-term continuous observation, giving it irreplaceable advantages over traditional remote sensing methods.

[0003] Currently, spaceborne lidar for Earth observation primarily employs near-zenith (small incident angle) laser pulse emission to achieve high geometric positioning accuracy of the ground spot. Atmospheric refraction is one of the main errors in lidar geometric positioning, significantly impacting its accuracy. Atmospheric delay during laser transmission can be calculated using a mapping function and the relationship between zenith delay and atmospheric delay. However, for ocean lidar satellites, which receive signals through backscattering, the reflected signal from the ocean surface is 1-2 orders of magnitude higher than the scattered signal from the water body when the laser is emitted near-zenith (small incident angle). This surface signal may saturate the detector, affecting the inversion of subsurface water signals. Therefore, spaceborne ocean lidar satellites typically have large incident angles. For such large incident angles, calculating atmospheric delay using a mapping function and the relationship between zenith delay introduces significant errors, necessitating a more accurate model to correct these errors. When a satellite orbits at an altitude of 400 km and the laser incident angle is 15°, the offset of the ground spot relative to the nadir point can reach hundreds of kilometers. If the temperature, humidity, and atmospheric pressure at the nadir point coordinates are still used, calculating atmospheric delay may introduce significant errors. Ray tracing models exist in two dimensions (two-dimensional) and three dimensions. The calculation principles of the two-dimensional and three-dimensional ray tracing models are consistent. In the two-dimensional ray tracing model, it is assumed that light travels from the satellite to the ground along a near-zenith direction, and only the vertical changes in atmospheric parameters are considered during the laser's journey from the top of the atmosphere to the ground. Three-dimensional ray tracing introduces vector calculations, continuously updating the three-dimensional coordinates of the ray vector's endpoint. By updating the laser's three-dimensional coordinates layer by layer, changes in atmospheric parameters in both the vertical and horizontal directions can be considered simultaneously, leading to a more accurate calculation of atmospheric refractive index. Therefore, the three-dimensional ray tracing model can accurately calculate the ray trajectory and atmospheric refractive index, achieving higher computational precision. Summary of the Invention

[0004] The purpose of this invention is to provide an atmospheric delay correction method based on three-dimensional ray tracing, which is applied to atmospheric delay correction of spaceborne lidar under large incident angle operating conditions, so as to improve the calculation accuracy of atmospheric delay correction.

[0005] The atmospheric delay correction method based on three-dimensional ray tracing is characterized by the following steps:

[0006] Step 1: Calculate the incident angle θ1 of the laser at the top of the atmosphere;

[0007] Step 2, divide the atmosphere into layers:

[0008] Using the altitude of the laser at the top of the atmosphere and the uncorrected altitude at the ground, the atmosphere is divided into layers every 30m, and the resulting number of layers is denoted as K.

[0009] Step 3: Calculate the atmospheric refractive index of the actual trajectory layer by layer:

[0010] a is the position vector of the laser in the i-th atmospheric layer (x) i ,y i ,z i ), convert it to latitude and longitude (B) i ,L i H i ), i = 1, ..., K-1, using B i and L i Find the 8 nearest grid points in the ERA5 data.

[0011] b Position vector (x) i ,y i ,z i The atmospheric pressure P, temperature tmp, and relative humidity rh of the grid points and its eight adjacent grid points are calculated using spline interpolation in the vertical direction and kriging interpolation in the horizontal direction.

[0012] c will determine the elevation H of the i-th layer. i Convert to altitude h i Based on the first three interpolation terms in step 3b and the atmospheric refractive index calculation formula, the atmospheric refractive index n of the laser in layer i+1 is obtained. i+1 ;

[0013] Step 4: Calculate the actual trajectory of the laser as it travels through the atmosphere in layers:

[0014] a uses the vector form of the law of refraction, based on the unit incident vector P of the i-th layer laser. i Position vector (x) i ,y i ,z i Laser incident angle θ i The atmospheric refractive index n in step 3c i+1and the refractive index n of the i-th layer of atmosphere i The unit incident vector P of the laser at layer i+1 is obtained. i+1 ;

[0015] b. Using vector calculations, based on the (i+1)th layer laser unit incident vector P i+1 i-layer position vector (x) i ,y i ,z i ) and the distance R0-(i+1)h from the (i+1)th atmospheric layer to the Earth's center of mass. x The travel distance ρ of the laser in the i-th layer of the atmosphere is obtained. i This allows us to obtain the position vector (x) of the laser in the (i+1)th layer. i+1 ,y i+1 ,z i+1 ), where R0 is the distance from the Earth's center to the top of the atmosphere, h x =30m;

[0016] c utilizes the laser unit incident vector P of the (i+1)th layer i+1 and position vector (x) i+1 ,y i+1 ,z i+1 ), to obtain the incident angle θ of the laser at layer i+1. i+1 ;

[0017] Step 5: Repeat steps 3 and 4 until the laser reaches the ground, thereby obtaining the position coordinates of the laser spot on the ground.

[0018] Step 6: Calculate the trajectory of the laser through the atmosphere under theoretical conditions. Based on the position vector of the laser point at the top of the atmosphere in Step 1 and the position coordinates of the laser spot on the ground in Step 5, obtain the theoretical distance of the laser from the top of the atmosphere to the ground.

[0019] Step 7: Using the atmospheric refractive index and the actual laser trajectory obtained in Steps 3 and 4, integrate the atmospheric refractive index along the actual laser trajectory and subtract the theoretical distance to obtain the atmospheric delay.

[0020] Step 1 is described in detail as follows:

[0021] The initial laser unit incident vector P0 and satellite orbit position vector (x0, y0, z0) of the spaceborne lidar are obtained. Using vector calculation, the position vector (x1, y1, z1) of the laser at the top of the atmosphere and the laser unit incident vector P1 of the first layer are obtained. The incident angle θ1 of the laser at the top of the atmosphere is obtained using the position vector (x1, y1, z1) and the laser unit incident vector P1.

[0022] Step 2 involves dividing the atmosphere into layers, as detailed below:

[0023] The position vector (x1, y1, z1) of the laser at the top of the atmosphere in step 1 is converted into latitude, longitude, and altitude (B1, L1, H1). The altitude H1 is then converted into elevation h1. The elevation h1 of the laser at the top of the atmosphere and the uncorrected elevation h1 at the ground level are then used to determine the final elevation. end The atmosphere is divided into layers every 30m, and the resulting number of layers is denoted as K.

[0024] Beneficial effects

[0025] For spaceborne lidar systems with large incident angles, calculating atmospheric delay using the relationship between the mapping function and zenith delay may produce significant errors. This invention can solve this problem. Based on three-dimensional ray tracing, this invention corrects atmospheric delay and can accurately calculate the ray trajectory and atmospheric refractive index during the laser's passage through the atmosphere, thus achieving higher calculation accuracy and providing a high-precision atmospheric delay correction. Attached Figure Description

[0026] Figure 1 This is the overall flowchart of the present invention.

[0027] Figure 2 This is a flowchart of the steps for "calculating the refractive index of the i+1 layer atmosphere using ERA5 data".

[0028] Figure 3 This is a diagram showing the optical path of a laser beam from a satellite to the top of the atmosphere.

[0029] Figure 4 This is the optical path diagram of a laser beam traveling from the i-th layer to the i+1-th layer within the atmosphere.

[0030] Figure 5 Figure: Curves showing the variation of atmospheric pressure, temperature, relative humidity, and saturated water vapor pressure with geopotential height.

[0031] in, Figure 5 (a) Figure 5 (b) Figure 5 (c) The graphs show the changes in atmospheric pressure, temperature, relative humidity, and refractive index with geopotential height. Figure 5 (d) shows the curve of saturated water vapor pressure as a function of potential height.

[0032] Figure 6 Atmospheric refractive index as a function of geopotential height.

[0033] Figure 7 Comparison of atmospheric delay calculations using 2D and 3D ray tracing. Detailed Implementation

[0034] like Figure 1 , 2As shown, the atmospheric delay correction method based on three-dimensional ray tracing is characterized by the following steps:

[0035] Step 1: Obtain the initial laser unit incident vector P0 and satellite orbit position vector (x0, y0, z0) of the spaceborne lidar. Using vector calculation, obtain the position vector (x1, y1, z1) of the laser at the top of the atmosphere and the first layer laser unit incident vector P1. Using the position vector (x1, y1, z1) and the laser unit incident vector P1, obtain the incident angle θ1 of the laser at the top of the atmosphere.

[0036] Step 2, divide the atmosphere into layers:

[0037] The position vector (x1, y1, z1) of the laser at the top of the atmosphere in step 1 is converted into latitude, longitude, and altitude (B1, L1, H1). The altitude H1 is then converted into elevation h1. The elevation h1 of the laser at the top of the atmosphere and the uncorrected elevation h1 at the ground level are then used to determine the final elevation. end The atmosphere is divided into layers every 30m, and the resulting number of layers is denoted as K.

[0038] Step 3: Calculate the atmospheric refractive index of the actual trajectory layer by layer:

[0039] a is the position vector of the laser in the i-th atmospheric layer (x) i ,y i ,z i ), convert it to latitude and longitude (B) i ,L i H i ), i = 1, ..., K-1, using B i and L i Find the 8 nearest grid points in the ERA5 data.

[0040] b Position vector (x) i ,y i ,z i The atmospheric pressure P, temperature tmp, and relative humidity rh of the grid points and its eight adjacent grid points are calculated using spline interpolation in the vertical direction and kriging interpolation in the horizontal direction.

[0041] c will determine the elevation H of the i-th layer. i Convert to altitude h i Based on the first three interpolation terms in step 3b and the atmospheric refractive index calculation formula, the atmospheric refractive index n of the laser in layer i+1 is obtained. i+1 .

[0042] Step 4: Calculate the actual trajectory of the laser as it travels through the atmosphere in layers:

[0043] a uses the vector form of the law of refraction, based on the unit incident vector P of the i-th layer laser. i Position vector (x)i ,y i ,z i Laser incident angle θ i The atmospheric refractive index n in step 3c i+1 and the refractive index n of the i-th layer of atmosphere i The unit incident vector P of the laser at layer i+1 is obtained. i+1 .

[0044] b. Using vector calculations, based on the (i+1)th layer laser unit incident vector P i+1 i-layer position vector (x) i ,y i ,z i ) and the distance R0-(i+1)h from the (i+1)th atmospheric layer to the Earth's center of mass. x The travel distance ρ of the laser in the i-th layer of the atmosphere is obtained. i This allows us to obtain the position vector (x) of the laser in the (i+1)th layer. i+1 ,y i+1 ,z i+1 ), where R0 is the distance from the Earth's center to the top of the atmosphere, h x =30m.

[0045] c utilizes the laser unit incident vector P of the (i+1)th layer i+1 and position vector (x) i+1 ,y i+1 ,z i+1 ), to obtain the incident angle θ of the laser at layer i+1. i+1 .

[0046] Step 5: Repeat steps 3 and 4 until the laser reaches the ground, thereby obtaining the position coordinates of the laser spot on the ground.

[0047] Step 6: Calculate the trajectory of the laser through the atmosphere under theoretical conditions. Based on the position vector of the laser point at the top of the atmosphere in Step 1 and the position coordinates of the laser spot on the ground in Step 5, obtain the theoretical distance of the laser from the top of the atmosphere to the ground.

[0048] Step 7: Using the atmospheric refractive index and the actual laser trajectory obtained in Steps 3 and 4, integrate the atmospheric refractive index along the actual laser trajectory and subtract the theoretical distance to obtain the atmospheric delay.

[0049] The rationale for this method will first be explained based on theoretical derivation.

[0050] 1. Theoretical Formula for Ranging Error of Spaceborne LiDAR

[0051] During the transmission of laser light through the atmosphere, changes in the atmospheric refractive index cause an additional optical path difference, resulting in atmospheric delay error. The single-pass atmospheric delay of a lidar system is expressed as...

[0052]

[0053] The first integral term represents the integral along the actual propagation path of the laser pulse in the atmosphere with respect to the actual atmospheric refractive index, while the second integral term represents the integral along the theoretical propagation path of the laser pulse in the atmosphere with respect to the theoretical atmospheric refractive index.

[0054] In a real atmospheric environment, the atmospheric refractive index is very close to 1. For ease of analysis and calculation, we introduce atmospheric refractive index: Atmospheric refractive index = Atmospheric refractive index * 10-1 -6 +1. Under standard atmospheric conditions, the atmospheric pressure P = 1013.25 hPa, the temperature T = 373.15 K, the CO2 content is 0.00375%, and the water vapor pressure P... w At 0.0 hPa, the atmospheric refractive index is expressed as follows:

[0055]

[0056] For a lidar system with a wavelength λ = 532 nm, the atmospheric refractive index is 293.5405 ppm under standard atmospheric conditions, which is equivalent to an atmospheric refractive index of 1.000293.

[0057] The expression for actual atmospheric refractive index is:

[0058]

[0059] The meteorological dataset does not provide the partial pressure of water vapor, P. w The value of P is related to the gas temperature T and relative humidity rh. The saturated vapor pressure P is obtained using Chebyshev polynomial iteration. s The calculation formula is as follows:

[0060]

[0061] Among them, P b =1000Pa, T max =648K,T min =273K, E i (x) is a Chebyshev polynomial with coefficients a and b respectively. i ={2794.027, 1430.604, -18.234, 7.674, -0.022, 0.263, 0.146, 0.055, 0.033, 0.015, 0.013}, P s The expression for the partial pressure of water vapor, where saturation water vapor pressure is given, is as follows:

[0062] P w =RhP s (6)

[0063] When performing ray tracing, altitude h and ground elevation h are required. e The relationship between altitude h and the altitude is as follows:

[0064] h = h e -N(7)

[0065] Where N represents the geoid difference.

[0066] The ERA5 global reanalysis data used to calculate atmospheric refractive index provides data for the corresponding pressure layers via geopotential height. For ease of calculation, the geopotential height H needs to be converted to altitude h.

[0067]

[0068] in, Where is the geographical latitude, R is the average radius of the Earth, R = 6371009 m, and the constant term g eq = 9.7803267715 m / s 2 k = 0.001931851353, e 2 =0.00669438002290.

[0069] 2. Assumptions about the atmosphere

[0070] Due to the Earth's gravity, the vast majority of gases are concentrated on the Earth's surface, with 90% of the atmosphere located below 30 km altitude. Above 30 km, the atmosphere is very thin and its effect on lasers is negligible. In three-dimensional ray tracing, we assume the top of the atmosphere is at an atmospheric pressure of 1 mb (altitude much greater than 30 km), and that the atmosphere is divided into layers every 30 m. The laser propagates in a straight line within the same layer, and refraction occurs only at the interface between adjacent atmospheric layers, obeying the law of refraction. We ignore variations in atmospheric altitude and assume the distance from each atmospheric layer to the Earth's center of mass is the same.

[0071] 3. The actual trajectory of a laser as it travels through the atmosphere

[0072] 3.1 Calculation of the optical path of laser light from satellite to the top of the atmosphere

[0073] The laser travels in a straight line from laser emission point M on the satellite to point N at the top of the atmosphere, as shown in the optical path diagram. Figure 3 As shown, the unit pointing vector of the laser (in the ECEF coordinate system) is P0, and the orbital position vector of the satellite (in the ECEF coordinate system) is OM = (x i ,y i ,z iLet R0 be the distance from the top of the atmosphere to the Earth's center of mass O. Let |MN| = ρ0, then MN = ρ0P0. According to vector calculations...

[0074]

[0075] Equation (9) yields a quadratic equation in ρ0, the value of which is as follows:

[0076]

[0077] Based on P0 and ρ0, the unit incident vector MN = P1 and the position vector NO = (x1, y1, z1) of the laser at the top of the atmosphere are obtained. The incident angle θ1 of the laser at the top of the atmosphere is calculated from the vectors ON and MN, as follows:

[0078]

[0079] 3.2 Calculation of the refractive index of the laser in the (i+1)th atmospheric layer

[0080] The position vector of the laser in the i-th layer of atmosphere (x) i ,y i ,z i In the ECEF coordinate system, transform it to the latitude, longitude, and altitude coordinate system, i.e., (B i L i H i ).

[0081] According to the laser in layer i, B i and L i In ERA5 data, atmospheric pressure P, temperature T, relative humidity rh, and geopotential height H of the eight nearest grid points to the laser point were obtained. The geopotential height H was then converted to altitude h. Using spline interpolation of atmospheric pressure P, temperature T, and relative humidity rh in the vertical direction, and employing the formula for calculating atmospheric refractive index, the relationship between atmospheric refractive index and altitude at the eight nearest points was obtained. Then, Kriging interpolation was used spatially to obtain the relationship between atmospheric refractive index and altitude at the laser point. Based on the altitude h of the laser point at layer i... i The atmospheric refractive index n of the (i+1)th layer is obtained. i+1 (The refractive index of the i+1 layer of atmosphere is expressed as the refractive index at the bottom of the i-th layer of atmosphere).

[0082] 3.3 Calculation of the laser beam trajectory in the (i+1)th layer of the atmosphere

[0083] like Figure 4 As shown, the laser beam is refracted at the interface between the two atmospheres, θ i and θ i+1 Let h represent the incident angles of the i-th and (i+1)-th atmospheric layers, respectively.x h is the width of each atmospheric layer. x =30m, the unit incident vector of the laser at layer i+1 is P i+1 (i.e., the unit refraction vector of the laser in the i-th layer), the position vector (normal vector) of the laser in the i-th layer is (x i y i , z i The distance from the (i+1)th atmospheric layer to the Earth's center of mass O is R0-(i+1)h. x The unit incident vector P of the laser at layer i+1 i+1 The calculation can be performed using the vector form of the law of refraction, as shown in the following formula:

[0084]

[0085] The distance the laser travels in the i-th layer of the atmosphere is ρ. i , ρ i+1 The calculation formula is as follows:

[0086]

[0087] From this, we obtain a value about ρ i A quadratic equation in one variable has the following value:

[0088]

[0089] From the above, we can see that the distance ρ the laser travels in the (i+1)th layer is... i+1 Therefore, the position vector (x) of the laser in the (i+1)th layer can be obtained. i+1 ,y i+1 ,z i+1 ).

[0090] The formula for calculating the incident angle of the laser at the (i+1)th moment is as follows:

[0091]

[0092] The atmospheric delay expression for three-dimensional ray tracing calculation of a 4-satellite-borne lidar system is as follows:

[0093]

[0094] Specific examples

[0095] Calculate the atmospheric path delay at a location of 21.4539°S, 159.2450°E, the initial laser unit incident vector of the spaceborne lidar (0.9506738, -0.070408, 0.3020961), and the satellite orbital position vector (-5967486.65, 2134782.03, -2454362.63). The distance from the top of the atmosphere to the center of mass of the Earth is 6,423,754.71 m. Using formulas (9) and (10), the unit incident vector (0.9506738, -0.070408, 0.3020961) and position vector (-5599883.23, 2107557.01, -2337549.08) of the laser in the atmosphere are obtained. The distance from the top of the atmosphere to the center of mass of the Earth is 6,423,754.71 m.

[0096] Converting the laser's position vector at the top of the atmosphere to latitude, longitude, and altitude (-21.2910, 160.3161, 48416.81), the laser's altitude at the top of the atmosphere is 48416.81, resulting in an altitude of 48406.05m. The uncorrected ground latitude, longitude, and altitude coordinates are (-21.4935, 159.2450, 53.3370), giving the laser's altitude at ground level as -2.33m. Using the laser's altitude at the top of the atmosphere and the ground level, and dividing the atmosphere into 1613 layers at 30-meter intervals, we can achieve this.

[0097] Using a laser, the atmospheric pressure, temperature, relative humidity, and geopotential height of eight neighboring grid points in ERA5 meteorological data are located at the longitude and latitude of the top of the atmosphere. Since ERA5 data uses 37 pressure layers to provide corresponding temperature, relative humidity, and geopotential height data, for ease of calculation, the geopotential height at these eight grid points is converted to altitude. Therefore, the corresponding atmospheric pressure, temperature, and relative humidity data can be obtained using the altitude. For each of these eight grid points, in the vertical direction, Kriging interpolation is performed on the atmospheric pressure, temperature, and relative humidity using the altitude. The atmospheric refractive index at different altitudes at the eight grid points is then obtained using the atmospheric refractive index calculation formula. In the horizontal direction, Kriging interpolation is performed on the atmospheric refractive index of the eight grid points using the longitude and latitude of the laser at the first layer of the atmosphere. Then, the atmospheric refractive index of the laser at the first layer of the atmosphere is obtained using the altitude of the first layer. The atmospheric refractive index outside the atmosphere is assumed to be 1.

[0098] With an incident angle of 15.8927° in the first layer of the atmosphere, an atmospheric refractive index of 0.3018 ppm, and an atmospheric refractive index of 1.0000003, the laser unit incident vector, the position vector, and the refractive index of the 0th layer are used to obtain the laser unit refractive vector (0.950674, -0.070407, 0.302096) using the vector form of the law of refraction. This is the laser unit incident vector in the second layer of the atmosphere.

[0099] Using vector calculations, based on the laser unit incident vector in the second layer, the position vector in the first layer, and the distance from the second layer of atmosphere to the Earth's center of mass (6423724.70 m), the laser's travel distance in the first layer of atmosphere is calculated to be 31.1923 m. This leads to the laser's position vector in the second layer: (-5599853.58, 2107554.81, -2337539.67). Using the laser unit incident vector and position vector in the second layer, the laser's incident angle in the second layer is calculated to be 15.8927°.

[0100] The light travel distance and atmospheric refractive index from the first layer of the atmosphere to the 1613th layer were repeatedly calculated. The atmospheric refractive index of the 1613th layer was 265.9154 ppm, the atmospheric refractive index was 1.0002659, the laser incident angle was 16.0121°, and the light travel distance was 31.2109 m. Using the position vector of the light at the top of the atmosphere and the position vector of the ground, the theoretical distance of the laser was obtained as 50328.6429765252 m. Using formula (16), the final atmospheric delay was obtained as 2.2439 m.

[0101] Using meteorological data provided by ERA5, meteorological data for the upper atmosphere at 21.4539°S, 159.2450°E were obtained, such as... Figure 5 The graph shows the changes in atmospheric pressure, temperature, relative humidity, and saturated vapor pressure with geopotential height. Figure 5 (a) Figure 5 (b) Figure 5 (c) The graphs show the changes in atmospheric pressure, temperature, relative humidity, and refractive index with geopotential height. The saturated vapor pressure is obtained using Chebyshev polynomial iteration. Figure 5 (d) shows the curve of saturated water vapor pressure as a function of potential height.

[0102] Based on the above meteorological data and the formula for calculating atmospheric refractive index, the curve of atmospheric refractive index as a function of geopotential height is obtained. Figure 6 As shown.

[0103] Calculating atmospheric delay requires using uncorrected latitude, longitude, elevation, laser pointing, and satellite orbital position data of the laser footprint point, as well as ERA5 meteorological data.

[0104] The correction value for atmospheric delay was calculated using satellite-borne data and ERA5 data, such as... Figure 7 As shown. The initial incident angle of the laser is 8°. There is a certain difference between the atmospheric delay calculated by two-dimensional ray tracing and three-dimensional ray tracing. The calculation result of three-dimensional ray tracing is on average 0.0242m larger than that of two-dimensional ray tracing.

Claims

1. An atmospheric delay correction method based on three-dimensional ray tracing, characterized in that... The steps include the following: Step 1: Calculate the incident angle of the laser at the top of the atmosphere. θ 1; Step 2, divide the atmosphere into layers: Using the altitude of the laser at the top of the atmosphere and the uncorrected altitude at the ground, the atmosphere is divided into layers every 30m, and the resulting number of layers is denoted as K. Step 3: Calculate the atmospheric refractive index of the actual trajectory layer by layer: a For laser in the first i The position vector of the layer of atmosphere ( x i , y i , z i ), convert it to latitude and longitude ( B i , L i , H i ), i =1,…,K-1, using B i and L i Find the 8 nearest grid points in the ERA5 data; b position vector ( x i , y i , z i The atmospheric pressure P, temperature tmp, and relative humidity rh of the grid points and its eight adjacent grid points are obtained by spline interpolation in the vertical direction and kriging interpolation in the horizontal direction. c will be the i Elevation of the layer H i Convert to altitude h i Based on the first three interpolation terms in step 3b and the atmospheric refractive index calculation formula, the laser beam at the [missing term] is obtained. i Atmospheric refractive index of +1 layer n i+1 ; Step 4: Calculate the actual trajectory of the laser as it travels through the atmosphere in layers: a uses the vector form of the law of refraction, according to the first... i Layer laser unit incident vector P i Position vector ( x i , y i , z i ), laser incident angle θ i Atmospheric refractive index in step 3c n i+1 and i Atmospheric refractive index n i , obtained the i +1 layer laser unit incident vector P i+1 ; b uses vector calculations, according to the first... i +1 layer laser unit incident vector P i+1 i-layer position vector ( x i , y i , z i ) and the i The distance from +1 layer of atmosphere to the Earth's center of mass R 0-( i +1) h x The laser was obtained in the first... i The distance traveled by the layer of atmosphere ρ i Thus, the laser is obtained in the first... i +1 layer position vector ( x i+1 , y i+1, z i+1 ),in R 0 represents the distance from the Earth's center to the top of the atmosphere. h x =30m; c. Using the first i Laser unit incident vector at +1 layer P i+1 and position vector ( x i+1 , y i+1, z i+1 ), to obtain the laser in the first i +1 layer incident angle θ i+1 ; Step 5: Repeat steps 3 and 4 until the laser reaches the ground, thereby obtaining the position coordinates of the laser spot on the ground. Step 6, calculate the trajectory of the laser through the atmosphere under theoretical conditions: based on the position vector of the laser point at the top of the atmosphere in Step 1 and the position coordinates of the laser spot on the ground in Step 5, obtain the theoretical distance of the laser from the top of the atmosphere to the ground. Step 7: Using the atmospheric refractive index and the actual laser trajectory obtained in Steps 3 and 4, integrate the atmospheric refractive index along the actual laser trajectory and subtract the theoretical distance to obtain the atmospheric delay.

2. The atmospheric delay correction method based on three-dimensional ray tracing as described in claim 1, characterized in that... Step 1 is described in detail as follows: Obtain the initial laser unit incident vector of the spaceborne lidar P 0. Satellite orbital position vector ( x 0, y 0, z 0), using vector calculations, the position vector of the laser's atmospheric top is obtained ( x 1, y 1, z 1) and the unit incident vector of the first layer laser P 1. Using position vectors ( x 1, y 1, z 1) and laser unit incident vector P 1. Obtain the incident angle of the laser at the top of the atmosphere. θ 1。 3. The atmospheric delay correction method based on three-dimensional ray tracing as described in claim 1, characterized in that... Step 2 involves dividing the atmosphere into layers, as detailed below: The position vector of the laser at the top of the atmosphere in step 1 ( x 1, y 1, z 1) Convert to latitude and longitude (latitude and longitude) B 1, L 1, H 1) Further elevation H 1 converted to high altitude h 1. Utilizing lasers at high altitudes above the atmosphere h 1 and the uncorrected altitude on the ground h end The atmosphere is divided into layers every 30m, and the resulting number of layers is denoted as K.

Images