Atmospheric delay correction method based on three-dimensional ray tracing

Abstract

The invention discloses an atmospheric delay correction method based on three-dimensional ray tracing. The method comprises the following steps: obtaining an incident angle of laser on an atmosphere top by using a position vector of the laser atmosphere top and a unit incident vector of a first layer of laser; layering the atmosphere; the atmospheric refractive index of the actual track is calculated hierarchically, and the actual track of the laser passing through the atmosphere is calculated hierarchically until the laser reaches the ground, so that the position coordinates of the laser spot on the ground are obtained; calculating a track of laser passing through the atmosphere under theoretical conditions; and by using the obtained atmospheric refractive index and the actual laser track, subtracting the theoretical distance from the atmospheric refractive index integral along the actual laser track to obtain the atmospheric delay. For a satellite-borne laser radar system under the condition of a large incident angle, a large error can be generated when the atmospheric delay is calculated by adopting the relation between a mapping function and zenith delay, and the method can accurately calculate the light track and the atmospheric refractive index in the process that laser passes through the atmospheric layer, so that higher calculation precision can be obtained, and high-precision atmospheric delay correction can be provided.

Description

Technical Field

[0001] The present invention relates to an atmospheric delay correction method based on three-dimensional ray tracing, and in particular to a ranging error correction method for atmospheric delay generated when a satellite-borne laser radar observes the earth under large incident angle conditions, and belongs to the field of satellite remote sensing. Background Art

[0002] Satellite-borne ocean lidar can obtain subsurface information of the ocean and realize ocean remote sensing from two-dimensional plane to three-dimensional stereoscopic observation. It is not limited by sunlight and can achieve long-term continuous observation. Compared with traditional remote sensing methods, it has irreplaceable advantages.

[0003] At present, satellite-borne laser radars for earth observation mainly use near-zenith (small incident angle) directions to emit laser pulses in order to obtain higher ground spot geometric positioning accuracy. The ranging error caused by atmospheric refraction is one of the main errors in laser radar geometric positioning, and the resulting error seriously affects the geometric positioning accuracy of ground spots. The atmospheric delay generated by the laser passing through the atmosphere can be calculated using the relationship between the mapping function and the zenith delay. However, for ocean laser radar satellites, laser radars receive signals through backscattering. When the laser is emitted in the near-zenith (small incident angle) direction, the reflected signal from the ocean surface will be 1-2 orders of magnitude higher than the water body scattered signal. The surface signal may saturate the detector and affect the inversion of the subsurface water body signal. Therefore, the incident angle of satellite-borne ocean laser radar satellites is generally large. For this large incident angle tilt, using the relationship between the mapping function and the zenith delay to calculate the atmospheric delay will bring large errors, and a more accurate model is needed to correct the error. When the satellite's orbital altitude is 400km and the laser incident angle is 15°, the ground spot can be offset by hundreds of kilometers relative to the subsatellite point. If the temperature, humidity and atmospheric pressure of the subsatellite point coordinates are still used, the calculation of atmospheric delay may introduce large errors. There are two types of ray tracing models: two-dimensional and three-dimensional. The calculation principles of the two-dimensional ray tracing model and the three-dimensional ray tracing model are the same. In the two-dimensional ray tracing model, it is believed that the light reaches the ground from the satellite along the near-zenith direction. In the process of the laser reaching the ground from the top of the atmosphere, only the changes in the atmospheric parameters in the vertical direction are considered. The three-dimensional ray tracing introduces vector calculation, which can continuously update the three-dimensional coordinates of the end of the light vector. The layer-by-layer update of the laser three-dimensional coordinates can simultaneously consider the changes in the atmospheric parameters in the vertical and horizontal directions, and thus more accurately calculate the atmospheric refractive index. Therefore, the three-dimensional ray tracing model can accurately calculate the light trajectory and the atmospheric refractive index, and the three-dimensional ray tracing model can obtain more accurate calculation accuracy. Summary of the invention

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

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

[0006] Step 1, calculate the incident angle θ1 of the laser at the top of the atmosphere;

[0007] Step 2, stratify the atmosphere:

[0008] Using the altitude of the laser at the top of the atmosphere and the uncorrected altitude on the ground, the atmosphere is layered every 30m, and the number of layers obtained is recorded as K;

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

[0010] a is the position vector (x i ,y i ,z i ), convert it into latitude and longitude height (B i ,L i ,H i ), i = 1, ..., K-1, using B i and L i Find 8 adjacent grid points in the ERA5 data.

[0011] b position vector (x i ,y i ,z i ) and its eight adjacent grid points, using spline interpolation in the vertical direction and Kriging interpolation in the horizontal direction.

[0012] c The height H of the i-th layer i Convert to altitude h i According to the interpolation of the first three items in step 3b and the atmospheric refractive index calculation formula, the atmospheric refractive index n of the laser at the i+1 layer is obtained. i+1 ;

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

[0014] a Using the vector form of the refraction law, according to the laser unit incident vector P of the i-th layer 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 of the i-th layer of atmosphere n i , and obtain the laser unit incident vector P of the i+1th layer i+1 ;

[0015] b Using vector calculation, according to the laser unit incident vector P of the i+1th layer i+1 , i-th layer position vector (x i ,y i ,z i ) and the distance from the i+1th layer of atmosphere to the Earth's center of mass R0-(i+1)h x , we get the distance ρ of the laser in the i-th layer of atmosphere i , and then the position vector (x i+1 ,y i+1 ,z i+1 ), where R0 is the distance from the center of the earth to the top of the atmosphere, h x =30m;

[0016] c uses the laser unit incident vector P of the i+1th layer i+1 and the position vector (x i+1 ,y i+1 ,z i+1 ), and the incident angle θ of the laser at the i+1 layer is obtained 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 surface;

[0018] Step 6, calculate the trajectory of the laser passing through the atmosphere under theoretical conditions, and obtain the theoretical distance from the top of the atmosphere to the ground according to 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;

[0019] Step 7, using the atmospheric refractive index obtained in steps 3 and 4 and the actual laser trajectory, integrate the atmospheric refractive index along the actual laser trajectory minus the theoretical distance to obtain the atmospheric delay.

[0020] The step 1 is specifically as follows:

[0021] Obtain the initial laser unit incident vector P0 of the satellite-borne laser radar and the satellite orbit position vector (x0, y0, z0), and use vector calculation to obtain the position vector (x1, y1, z1) of the laser top of the atmosphere and the first layer laser unit incident vector P1, and use the position vector (x1, y1, z1) and the laser unit incident vector P1 to obtain the incident angle θ1 of the laser at the top of the atmosphere;

[0022] The step 2 stratifies the atmosphere, as follows:

[0023] Convert the position vector (x1, y1, z1) of the laser at the top of the atmosphere in step 1 to latitude and longitude (B1, L1, H1), and further convert the elevation H1 to the altitude h1. end The atmosphere is stratified every 30 m, and the number of layers obtained is recorded as K.

[0024] Beneficial Effects

[0025] For satellite-borne lidar systems under large incident angle conditions, calculating the atmospheric delay using the relationship between the mapping function and the zenith delay may result in large errors. The present invention can solve this problem. The present invention corrects the atmospheric delay based on three-dimensional ray tracing, and can accurately calculate the light trajectory and atmospheric refractive index of the laser during its passage through the atmosphere, thereby achieving higher calculation accuracy and providing a high-precision atmospheric delay correction amount. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Figure 1 It is the overall flow chart of the present invention.

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

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

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

[0030] Figure 5 Figure 3. The atmospheric pressure, temperature, relative humidity and saturated water vapor pressure change with the potential height.

[0031] in, Figure 5 (a) Figure 5 (b) Figure 5 (c) are graphs showing the changes in atmospheric pressure, temperature, relative humidity and refractivity with geopotential height. Figure 5 (d) is a graph showing the variation of saturated water vapor pressure with potential height.

[0032] Figure 6 Graph showing atmospheric refractivity changing with potential height.

[0033] Figure 7 Comparison of atmospheric delay calculated using 2D and 3D ray tracing. DETAILED DESCRIPTION

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

[0035] Step 1, obtain the initial laser unit incident vector P0 of the satellite-borne laser radar and the satellite orbit position vector (x0, y0, z0), use vector calculation to obtain the position vector (x1, y1, z1) of the laser top of the atmosphere and the first layer laser unit incident vector P1, and use the position vector (x1, y1, z1) and the laser unit incident vector P1 to obtain the incident angle θ1 of the laser at the top of the atmosphere.

[0036] Step 2, stratify the atmosphere:

[0037] Convert the position vector (x1, y1, z1) of the laser at the top of the atmosphere in step 1 to latitude and longitude (B1, L1, H1), and further convert the elevation H1 to the altitude h1. end The atmosphere is stratified every 30 m, and the number of layers obtained is recorded as K.

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

[0039] a is the position vector (x i ,y i ,z i ), convert it into latitude and longitude height (B i ,L i ,H i ), i = 1, ..., K-1, using B i and L i Find 8 adjacent grid points in the ERA5 data.

[0040] b position vector (x i ,y i ,z i ) and its eight adjacent grid points, using spline interpolation in the vertical direction and Kriging interpolation in the horizontal direction.

[0041] c The height H of the i-th layer i Convert to altitude h i According to the interpolation of the first three items in step 3b and the atmospheric refractive index calculation formula, the atmospheric refractive index n of the laser at the i+1 layer is obtained. i+1 .

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

[0043] a Using the vector form of the refraction law, according to the laser unit incident vector P of the i-th layer i , position vector (xi ,y i ,z i ), laser incident angle θ i , the atmospheric refractive index n in step 3c i+1 and the refractive index of the i-th layer of atmosphere n i , and obtain the laser unit incident vector P of the i+1th layer i+1 .

[0044] b Using vector calculation, according to the laser unit incident vector P of the i+1th layer i+1 , i-th layer position vector (x i ,y i ,z i ) and the distance from the i+1th layer of atmosphere to the Earth's center of mass R0-(i+1)h x , we get the distance ρ of the laser in the i-th layer of atmosphere i , and then the position vector (x i+1 ,y i+1 ,z i+1 ), where R0 is the distance from the center of the earth to the top of the atmosphere, h x =30m.

[0045] c uses the laser unit incident vector P of the i+1th layer i+1 and the position vector (x i+1 ,y i+1 ,z i+1 ), and the incident angle θ of the laser at the i+1 layer is obtained 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 surface.

[0047] Step 6, calculate the trajectory of the laser through the atmosphere under theoretical conditions, and obtain the theoretical distance from the top of the atmosphere to the ground 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 surface in step 5.

[0048] Step 7, using the atmospheric refractive index obtained in steps 3 and 4 and the actual laser trajectory, integrate the atmospheric refractive index along the actual laser trajectory minus the theoretical distance to obtain the atmospheric delay.

[0049] The following is based on theoretical derivation to illustrate the rationality of this method.

[0050] 1. Theoretical formula of satellite-borne laser radar ranging error

[0051] When the laser passes through the atmosphere, the change in atmospheric refractive index will cause additional optical path difference and produce atmospheric delay error. The one-way atmospheric delay of the lidar system is expressed as

[0052]

[0053] The first integral term represents the integration of the actual atmospheric refractive index along the actual transmission path of the laser pulse in the atmosphere, and the second integral term represents the integration of the theoretical atmospheric refractive index along the theoretical transmission path of the laser pulse in the atmosphere.

[0054] In the real atmosphere, the atmospheric refractive index is very close to 1. In order to facilitate analysis and calculation, the atmospheric refractive index is introduced: atmospheric refractive index = atmospheric refractive index * 10 -6 +1. Under standard atmospheric conditions, atmospheric pressure P = 1013.25 hPa, temperature T = 373.15 K, CO2 content is 0.00375%, water vapor pressure P w =0.0hPa, the expression of atmospheric refractivity is:

[0055]

[0056] For a laser radar system with a wavelength of λ=532nm, under standard atmospheric conditions, the atmospheric refractive index is 293.5405ppm, that is, the atmospheric refractive index is 1.000293.

[0057] The expression of actual atmospheric refractivity is:

[0058]

[0059] The water vapor partial pressure P is not given in the meteorological data set. w The value of , which is related to the gas temperature T and relative humidity rh, is obtained by iteratively using the Chebyshev polynomial to obtain the saturated water vapor pressure P 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, and its coefficients are 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 is the saturated water vapor pressure, and the expression of water vapor partial pressure is as follows:

[0062] P w =RhP s (6)

[0063] When doing ray tracing, you need to use the altitude h and the surface elevation h e The relationship with the altitude h is as follows:

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

[0065] Where N is the geoid difference.

[0066] The ERA5 global reanalysis data used to calculate the atmospheric refractivity gives the data on the corresponding pressure layer through the potential height. In order to facilitate the calculation, the potential height H needs to be converted into the altitude h.

[0067]

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

[0069] 2. Assumptions about the atmosphere

[0070] Because the atmosphere is affected by the earth's gravity, most of the gases are concentrated on the earth's surface, and 90% of the atmosphere is concentrated below 30km above sea level. For places above 30km above sea level, the atmosphere is very thin and its effect on the laser can be ignored. When performing three-dimensional ray tracing, it is assumed that the height of the top of the atmosphere is the position where the atmospheric pressure is 1mb (the altitude is much greater than 30km), and the atmosphere is layered every 30m. The laser propagates in a straight line within the same layer, and the laser is only refracted at the interface of adjacent atmospheric layers, obeying the law of refraction. Ignore the changes in the altitude of the atmosphere, and assume that the distance from each layer of the atmosphere to the center of mass of the earth is the same.

[0071] 3. The actual trajectory of the laser through the atmosphere

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

[0073] The laser travels from the laser emission point M on the satellite to the top point N of the atmosphere in a straight line. The optical path diagram is as follows: 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 i), the distance from the top of the atmosphere to the center of mass of the earth O is R0, let |MN|=ρ0, then MN=ρ0P0, according to the calculation of vector

[0074]

[0075] From formula (9), we can get a quadratic equation about ρ0, whose value is as follows:

[0076]

[0077] According to P0 and ρ0, the laser 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 by the vector ON and MN. The formula is as follows:

[0078]

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

[0080] The position vector (x i ,y i ,z i ) in the ECEF coordinate system, transform it to the latitude and longitude coordinate system, that is, (B i , L i , H i ).

[0081] According to the laser B at the i layer i and L i , the atmospheric pressure P, temperature T, relative humidity rh and potential height H data of the 8 grid points closest to the laser point are obtained in the ERA5 data, and the potential height H is converted to the altitude h. The atmospheric pressure P, temperature T, and relative humidity rh are interpolated vertically using the spline interpolation method. The atmospheric refractive index calculation formula is used to obtain the relationship between the atmospheric refractive index of the 8 neighboring points and the altitude. Then, the Kriging interpolation method is used in space to obtain the relationship between the atmospheric refractive index of the laser point and the altitude. According to the altitude h of the laser point at the i-th layer, the atmospheric refractive index of the laser point is calculated. i , and get the atmospheric refractive index n of the i+1th layer i+1 (The refractive index of the i+1 layer of atmosphere is expressed as the refractive index of the bottom of the i layer of atmosphere).

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

[0083] like Figure 4 As shown, the laser light is refracted at the interface between the two atmospheric layers, θ i and θ i+1 denote the incident angles of the i-th and i+1-th atmospheric layers, respectively, and hx is the width of each atmospheric layer, h x =30m, the unit incident vector of the laser at the i+1 layer is P i+1 (i.e., the laser unit refraction vector of 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+1th layer of atmosphere to the Earth's center of mass O is R0-(i+1)h x , the unit incident vector P of the laser at the i+1 layer i+1 It can be calculated using the vector form of the law of refraction, and the calculation formula is as follows:

[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] By getting a i A quadratic equation in one variable has the following value:

[0088]

[0089] From the above, we can know that the distance ρ that the laser travels in the i+1th layer is i+1 , and then we can get the position vector (x i+1 ,y i+1 ,z i+1 ).

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

[0091]

[0092] The expression for calculating atmospheric delay by three-dimensional ray tracing of satellite-borne lidar system is:

[0093]

[0094] Specific Examples

[0095] Calculate the atmospheric path delay at a position 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 orbit position vector (-5967486.65, 2134782.03, -2454362.63). The distance from the laser at the top of the atmosphere to the center of mass of the earth is 6423754.71 m. Using formulas (9) and (10), we can obtain the laser unit incident vector (0.9506738, -0.070408, 0.3020961) and position vector (-5599883.23, 2107557.01, -2337549.08) in the atmosphere. The distance from the top of the atmosphere to the center of mass of the earth is 6423754.71 m.

[0096] Convert the position vector of the laser at the top of the atmosphere into longitude and latitude (-21.2910, 160.3161, 48416.81). The elevation of the laser at the top of the atmosphere is 48416.81, so the altitude of the laser at the top of the atmosphere is 48406.05m. The uncorrected longitude and latitude coordinates of the ground are (-21.4935, 159.2450, 53.3370), and the altitude of the laser on the ground is -2.33m. Using the altitudes of the laser at the top of the atmosphere and the ground, the atmosphere can be divided into 1613 layers at 30-meter intervals.

[0097] The atmospheric pressure, temperature, relative humidity and geopotential height of the 8 grid points near the ERA5 meteorological data are found by using the longitude and latitude of the laser at the top of the atmosphere. Since the ERA5 data gives the corresponding temperature, relative humidity and geopotential height data in 37 pressure layers, the geopotential height on these 8 grid points is converted to altitude for the convenience of calculation. Therefore, the altitude can be used to give the corresponding atmospheric pressure, temperature and relative humidity data. For these 8 grid points, in the vertical direction, the altitude is used to take Kriging interpolation for atmospheric pressure, temperature and relative humidity in the horizontal direction, and the atmospheric refractive index calculation formula is used to obtain the atmospheric refractive index at different altitudes on the 8 grid points. In the horizontal direction, the atmospheric refractive index of the 8 grid points is interpolated by Kriging using the longitude and latitude of the laser in the first layer of the atmosphere, and then the altitude of the first layer of the atmosphere is used to obtain the atmospheric refractive index of the laser in the first layer. The atmospheric refractive index outside the atmosphere is 1 by default.

[0098] The incident angle at the first layer of the atmosphere is 15.8927°, the atmospheric refractive index is 0.3018ppm, and the atmospheric refractive index is 1.0000003. The laser unit incident vector and position vector and the atmospheric refractive index of the 0th layer are 1. Using the vector form of the refraction law, the laser unit refraction vector (0.950674, -0.070407, 0.302096) is obtained, which is the laser unit incident vector at the second layer of the atmosphere.

[0099] By using vector calculation, according to the laser unit incident vector of the second layer, the position vector of the first layer and the distance from the second layer of atmosphere to the center of mass of the earth, which is 6423724.70m, the laser advance distance of 31.1923m in the first layer of atmosphere is obtained, and then the position vector of the laser in the second layer is obtained (-5599853.58, 2107554.81, -2337539.67). Using the laser unit incident vector and position vector of the second layer, the incident angle of the laser in the second layer is 15.8927.

[0100] Repeat the calculation of the distance traveled by the light from the first layer of the atmosphere to the 1613th layer and the atmospheric refractive index. The atmospheric refractive index of the 1613th layer is 265.9154ppm, the atmospheric refractive index is 1.0002659, the laser incident angle is 16.0121°, the distance traveled by the light is 31.2109m, and the position vector of the light at the top of the atmosphere and the position vector of the ground are used to obtain the theoretical distance of the laser, which is 50328.6429765252m. The final atmospheric delay is 2.2439m using formula (16).

[0101] Using the meteorological data provided by ERA5, the meteorological data of the upper layer at 21.4539°S, 159.2450°E were obtained, such as Figure 5 The atmospheric pressure, temperature, relative humidity and saturated water vapor pressure change with the potential height as shown in the graph, where: Figure 5 (a) Figure 5 (b) Figure 5 (c) shows the curves of atmospheric pressure, temperature, relative humidity and refractivity as a function of geopotential height. The saturated water vapor pressure is obtained by iteratively using Chebyshev polynomials. Figure 5 (d) is a graph showing the variation of saturated water vapor pressure with potential height.

[0102] According to the above meteorological data and the calculation formula of atmospheric refractivity, the curve of atmospheric refractivity changing with potential height is obtained. Figure 6 As shown.

[0103] The calculation of atmospheric delay requires the use of the latitude and longitude, elevation, laser pointing and satellite orbit position data of the laser footprint point without atmospheric delay correction and the meteorological data of ERA5.

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

Claims

1. Atmospheric delay correction method based on three-dimensional ray tracing, characterized by The steps include: Step 1: Calculate the incident angle of the laser at the top of the atmosphere θ 1; Step 2, stratify the atmosphere: Using the altitude of the laser at the top of the atmosphere and the uncorrected altitude on the ground, the atmosphere is layered every 30m, and the number of layers obtained is recorded as K; Step 3, calculate the atmospheric refractive index of the actual trajectory layer by layer: aFor the laser i The position vector of the atmosphere ( x i , y i , z i ), converting it to latitude and longitude height ( B i , L i , H i ), i =1,…,K-1, using B i and L i Find 8 adjacent grid points in the ERA5 data; b position vector ( x i , y i , z i ) and its neighboring 8 grid points, using spline interpolation in the vertical direction and Kriging interpolation in the horizontal direction; c will i Layer elevation H i Convert to altitude h i According to the interpolation of the first three items in step 3b and the atmospheric refractive index calculation formula, the laser is obtained at the i +1 layer of atmospheric refractive index n i+1 ; Step 4: Calculate the actual trajectory of the laser through the atmosphere in layers: a Using the vector form of the law of refraction, according to 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 Refractive index of the atmosphere n i , get the i +1 stack of laser unit incident vector P i+1 ; b Using vector calculation, according to i +1 stack of laser unit incident vector P i+1 , i-layer position vector ( x i , y i , z i ) and i +1 layer of atmosphere to the Earth's center of mass R 0-( i +1) h x , get the laser in the i Distance of advance in the atmosphere ρ i , and then the laser is obtained in the first i +1 layer position vector ( x i+1 , y i+1, z i+1 ),in R 0 is the distance from the center of the earth to the top of the atmosphere, h x =30m; c Using the i +1 layer of laser unit incident vector P i+1 and the position vector ( x i+1 , y i+1, z i+1 ), and the laser is obtained at 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 surface; Step 6, calculate the trajectory of the laser passing through the atmosphere under theoretical conditions: according to 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 obtained in steps 3 and 4 and the actual laser trajectory, integrate the atmospheric refractive index along the actual laser trajectory minus the theoretical distance to obtain the atmospheric delay.

2. The atmospheric delay correction method based on three-dimensional ray tracing as claimed in claim 1, characterized in that The step 1 is specifically as follows: Get the initial laser unit incident vector of the satellite-borne laser radar P 0. Satellite orbit position vector ( x 0, y 0, z 0), using vector calculation, we get the position vector of the laser top of the atmosphere ( x 1, y 1, z 1) and the first layer laser unit incident vector P 1. Using the position vector ( x 1, y 1, z 1) and the laser unit incident vector P 1 Get 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 claimed in claim 1, characterized in that The step 2 stratifies the atmosphere, as follows: 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 ( B 1, L 1, H 1), further H 1 Convert to altitude h 1. Use lasers at high altitudes at the top of the atmosphere h 1 and the uncorrected altitude at ground level h end The atmosphere is stratified every 30 m, and the number of layers obtained is recorded as K.

Images