A method for calculating atmospheric time delay based on ray tracing

Abstract

The application discloses a kind of atmospheric time delay calculation method based on ray tracing, the present application relates to atmospheric time delay calculation technical field, the specific content of the present application is: obtaining multi-source meteorological data, extracting temperature, specific humidity, three kinds of key meteorological parameters such as air pressure, by optimizing interpolation to complete vertical layering blank and correct grid deviation, to construct three-dimensional atmospheric parameter field of time and space continuity;Based on atmospheric physical characteristics and electromagnetic wave propagation principle, calculate atmospheric refractive index and its gradient;With station location as starting point, set the angle of incidence of electromagnetic wave, adopt nonlinear dynamic step length strategy to track propagation path, and compare actual and vacuum propagation time to obtain time delay result.The present application considers data integrity and timeliness, avoids experience model deviation, and improves the calculation precision of atmospheric time delay.

Description

Technical Field

[0001] This invention relates to the field of atmospheric time delay calculation technology, specifically to an atmospheric time delay calculation method based on ray tracing. Background Technology

[0002] Atmospheric time delay is the time difference between the propagation of electromagnetic waves in the atmosphere and in a vacuum. Its calculation accuracy directly affects the application effects in fields such as GNSS high-precision positioning and aerospace telemetry and control. However, existing technologies still have the following shortcomings:

[0003] First, relying heavily on low spatiotemporal resolution reanalysis data or local real-time data makes it difficult to balance the integrity of global coverage with the timeliness of short-term atmospheric changes, resulting in biases in data input.

[0004] Secondly, refractive index calculations often rely on empirical fitting models, which do not fully incorporate atmospheric physical properties and electromagnetic wave propagation principles. This makes it difficult to accurately capture the nonlinear changes in atmospheric parameters, thus affecting the accuracy of refractive index and its gradient calculations.

[0005] Third, the gaps in the vertical layering are mostly filled by simple linear interpolation, and the grid data deviation is not effectively corrected, resulting in poor spatiotemporal continuity and low reliability of the constructed three-dimensional atmospheric parameter field;

[0006] Fourth, the fixed step size strategy is often used in ray tracing, which makes it difficult to adapt to the spatial changes in refractive index gradient, resulting in a difficulty in balancing the accuracy and real-time performance of atmospheric time delay calculation.

[0007] Therefore, there is an urgent need for a method for calculating atmospheric time delay based on ray tracing. Summary of the Invention

[0008] To address the shortcomings of existing technologies, this invention provides an atmospheric delay calculation method based on ray tracing, which solves the problems of low accuracy, insufficient real-time performance, and inadequate scene adaptability in existing atmospheric delay calculation methods.

[0009] To achieve the above objectives, the present invention provides the following technical solution: a method for calculating atmospheric time delay based on ray tracing, comprising:

[0010] Step 1: Acquire global high spatiotemporal resolution reanalysis data and short-term high-frequency tropospheric real-time data, extract three key meteorological parameters: temperature, specific humidity, and air pressure, and use interpolation to fill in the meteorological parameters of each layer from the ground to the top of the effective influence layer of the atmosphere, and correct grid data deviations, and finally form a spatiotemporally continuous and fully covered three-dimensional atmospheric parameter field.

[0011] Step 2: Based on the three-dimensional atmospheric parameter field, and according to the atmospheric physical properties and the principle of electromagnetic wave propagation, calculate the atmospheric refractive index and its gradient.

[0012] Step 3: Starting from the geographical location of the target station, set the zenith angle and azimuth angle of the electromagnetic wave incident. Based on the atmospheric refractive index and its gradient, use a dynamic step size strategy to continuously update the propagation position and direction of the electromagnetic wave. Analyze the actual propagation path obtained by tracking and the propagation of the electromagnetic wave in the vacuum to calculate the atmospheric time delay result.

[0013] As a further aspect of the present invention, the interpolation method includes logarithmic interpolation and cubic spline interpolation: logarithmic interpolation is used for two key meteorological parameters, specific humidity and air pressure; cubic spline interpolation is used for temperature.

[0014] As a further aspect of the present invention, the specific operation for correcting grid data deviation is as follows:

[0015] Extract the mean value M of the same height parameter of 8 grid points in the 3×3 adjacent grid point set of the target grid point, obtain the parameter U of the height layer directly above the target grid point and the parameter L of the height layer directly below the target grid point, and calculate the correction value = (M+U+L) / 3 for the three types of parameters: temperature, air pressure and specific humidity.

[0016] Extract the height parameter R of the target grid point from the global high spatiotemporal resolution reanalysis data, and check whether R is within the fluctuation range of the height parameter in the same latitude zone of the data;

[0017] Extract the target grid point height parameter Q from the short-term high-frequency tropospheric real-time data, and verify whether Q is within the nominal accuracy range of the data observation equipment;

[0018] The final parameter value is determined according to the following rules: if both R and Q meet or do not meet the range requirements, the final parameter is the calculated correction value; if only R meets the range requirements, the final parameter = (calculated correction value + R) / 2; if only Q meets the range requirements, the final parameter value = (calculated correction value + Q) / 2.

[0019] As a further aspect of the present invention, according to the formula Calculate the atmospheric refractive index n, where q is the specific humidity, p is the total atmospheric pressure, e is the water vapor pressure, and p d Let T be the dry air pressure, T be the atmospheric temperature, and k1, k2, and k3 be constant coefficients.

[0020] As a further aspect of the present invention, according to the formula Update the propagation position and direction of the electromagnetic wave, where, h is the refractive index gradient. m Let m be the normal vector of the m-th iteration. Let r be the tangential vector of the m-th iteration. m Let r be the position vector of the m-th iteration, Δs be the iteration step size, and r be the position vector of the m-th iteration. m+1Let m+1 be the position vector of the (m+1)th iteration. h is the tangential vector of the (m+1)th iteration, and h0 is the normal unit vector. It is a tangential unit vector.

[0021] As a further aspect of the present invention, according to the formula Calculate the dynamic step size Δs, where Δs min The minimum step size is Δs. max This represents the maximum step size. α is the reference gradient, and α is the power coefficient.

[0022] As a further aspect of the present invention, the specific operation for calculating the atmospheric time delay result is as follows:

[0023] The complete actual path obtained by ray tracing is divided into several continuous line segments according to the iteration step size. For each line segment, its corresponding atmospheric refractive index is extracted. The propagation time of a single line segment is calculated by dividing the line segment length by the electromagnetic wave propagation speed. Then, the propagation times of all line segments are summed to obtain the actual total propagation time.

[0024] Calculate the straight-line distance between the target station and the satellite using their geographical location and spatial position; divide this straight-line distance by the speed of light in a vacuum to obtain the vacuum propagation time.

[0025] The atmospheric time delay is obtained by subtracting the vacuum propagation time from the calculated total propagation time.

[0026] This invention provides a method for calculating atmospheric time delay based on ray tracing, which has the following advantages compared with the prior art:

[0027] (1) This invention combines multi-source data fusion, optimized interpolation and grid deviation correction to take into account both the global coverage integrity and short-term timeliness of the data, effectively fill the gaps in vertical height layering and correct data deviations, making the constructed three-dimensional atmospheric parameter field more spatiotemporally continuous and reliable, and reducing the time delay calculation error caused by the data level from the source.

[0028] (2) This invention calculates atmospheric refractive index and its gradient based on atmospheric physical properties and electromagnetic wave propagation principle, abandons traditional empirical fitting model, accurately captures nonlinear changes and spatial distribution laws of atmospheric parameters, and improves the reliability of path tracking.

[0029] (3) The present invention uses a nonlinear dynamic step size strategy to perform ray tracing and adjusts the step size adaptively through the refractive index gradient, which not only ensures the path tracing accuracy in the high gradient region, but also improves the computational efficiency in the low gradient region, effectively balancing the accuracy and real-time performance of time delay calculation. Attached Figure Description

[0030] Figure 1 This is a flowchart of the steps of the present invention. Detailed Implementation

[0031] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0032] like Figure 1 This invention provides a method for calculating atmospheric time delay based on ray tracing;

[0033] As an embodiment of this application, the specific steps include the following:

[0034] Step 1: Acquire global high spatiotemporal resolution reanalysis data and short-term high-frequency tropospheric real-time data, extract three key meteorological parameters: temperature, specific humidity, and air pressure, and use interpolation to fill in the meteorological parameters of each layer from the ground to the top of the effective influence layer of the atmosphere, and correct grid data deviations, and finally form a spatiotemporally continuous and fully covered three-dimensional atmospheric parameter field.

[0035] Step 2: Based on the three-dimensional atmospheric parameter field, and according to the atmospheric physical properties and the principle of electromagnetic wave propagation, calculate the atmospheric refractive index and its gradient.

[0036] Step 3: Starting from the geographical location of the target station, set the zenith angle and azimuth angle of the electromagnetic wave incident. Based on the atmospheric refractive index and its gradient, use a dynamic step size strategy to continuously update the propagation position and direction of the electromagnetic wave. Analyze the actual propagation path obtained by tracking and the propagation of the electromagnetic wave in the vacuum to calculate the atmospheric time delay result.

[0037] As a second embodiment of this application, it is implemented based on the first embodiment, except that this embodiment includes:

[0038] Step 1: Acquire global high spatiotemporal resolution reanalysis data and short-term high-frequency tropospheric real-time data, and extract three key meteorological parameters: temperature, specific humidity, and air pressure.

[0039] The advantage of global high spatiotemporal resolution reanalysis data (such as reanalysis data obtained by the European Centre for Medium-Range Weather Forecasts ECMWF) is that it covers the entire globe and has complete vertical stratification, but its temporal resolution is relatively low and it is usually updated hourly. It is suitable for building an overall framework of three-dimensional atmosphere.

[0040] The core value of short-term, high-frequency tropospheric real-time data (such as tropospheric data obtained from the Crustal Dynamics Data Information System CDDIS) is its high temporal resolution. It is usually updated every minute and can accurately capture short-term atmospheric changes, making up for the time lag of reanalysis data.

[0041] If only global high spatiotemporal resolution reanalysis data is used, it cannot respond to short-term atmospheric changes and parameter distortion will occur in scenarios such as rainfall; if only short-term high-frequency tropospheric real-time data is used, upper atmospheric parameters are missing and a complete vertical parameter field cannot be constructed; therefore, combining the two types of data can take into account both global integrity and short-term timeliness.

[0042] Temperature, specific humidity, and air pressure are the three key meteorological parameters that form the basis for calculating atmospheric refractive index. Refractive index reflects the atmosphere's ability to refract electromagnetic waves, and its essence is determined by the atmosphere's thermal state, water vapor content, and molecular density, which correspond to temperature, specific humidity, and air pressure, respectively.

[0043] To address the gaps in the data at vertical height, an interpolation method is used to fill in the meteorological parameters of each layer within the effective atmospheric influence layer from the ground to the top of the effective atmospheric influence layer. The top of the effective atmospheric influence layer refers to the critical height at which the influence of the atmosphere on the refraction of electromagnetic waves is negligible. When the height exceeds this critical value, the atmosphere is extremely thin, the refractive index approaches 1, and the propagation path of electromagnetic waves is no longer deflected by atmospheric refraction. Subsequent ray tracing does not need to consider the atmospheric influence above this height.

[0044] Ground weather stations can only observe near-surface parameters, while upper-air data rely on satellite remote sensing or radiosondes. Among these, radiosondes have a low observation frequency, and satellite remote sensing data has a fixed grid interval, resulting in the lack of intermediate layer parameters between adjacent grid points.

[0045] Since the original observation data is discrete, in order to construct a continuous three-dimensional parameter field, it is necessary to divide it into vertical layers at fixed intervals. The layers between discrete observations naturally form gaps, which are filled by interpolation.

[0046] The interpolation methods include logarithmic interpolation and cubic spline interpolation.

[0047] Logarithmic interpolation targets two key meteorological parameters: specific humidity and air pressure. The vertical distribution of these two parameters is not linear, but rather exhibits an exponential pattern. Logarithmic interpolation can first transform the exponential relationship into a linear relationship before performing interpolation calculations, making it more effective than linear interpolation in correcting errors.

[0048] Cubic spline interpolation is used for temperature. The temperature at the bottom of the stratosphere may fluctuate nonlinearly, first decreasing and then increasing. Using linear interpolation will cause abrupt changes in parameters, while cubic spline interpolation constructs a smooth polynomial curve to ensure continuous parameter changes between adjacent layers. It can accurately capture the nonlinear trend of temperature while avoiding amplifying local observation noise.

[0049] By correcting grid data biases, a spatiotemporally continuous and fully covered three-dimensional atmospheric parameter field is finally formed.

[0050] The specific operation for correcting grid data deviation is as follows:

[0051] Extract the mean value M of the same height parameter of 8 grid points in the 3×3 neighboring grid point set of the target grid point. This value reflects the atmospheric consistency of the local area around the target grid point and avoids interference from a single adjacent abnormal grid point.

[0052] Extract the height layer parameter U immediately above and the height layer parameter L immediately below the target grid point;

[0053] Atmospheric parameters follow fixed physical laws in the vertical direction. For example, the temperature in the lower troposphere decreases with increasing altitude, and the air pressure decreases with increasing altitude. This trend is objective and stable. The vertical distance between the target grid point and the adjacent altitude layer directly above and below is the shortest, and the parameters have the strongest correlation, providing the most direct vertical constraints.

[0054] For the three parameters of temperature, air pressure, and specific humidity, calculate the correction value as (M+U+L) / 3 respectively.

[0055] The M value focuses on horizontal local consistency, while the U and L values ​​focus on vertical physical rationality. These three values ​​belong to two key dimensions of constraints and are equally important. Therefore, the arithmetic mean can balance the needs of both dimensions.

[0056] Extract the target grid point height parameter R from the global high spatiotemporal resolution reanalysis data, and verify whether R is within the fluctuation range of the same height parameter in the same latitude zone of the data. This range is clearly given in the description document of the global high spatiotemporal resolution reanalysis data.

[0057] Extract the target grid point height parameter Q from the short-term high-frequency tropospheric real-time data, and verify whether Q is within the nominal accuracy range of the data observation equipment, which is clearly marked in the equipment technical manual;

[0058] The final parameter is determined according to the following rules: if both R and Q meet or do not meet the range requirements, the final parameter is the calculated correction value; if only R meets the range requirements, the final parameter = (calculated correction value + R) / 2; if only Q meets the range requirements, the final parameter = (calculated correction value + Q) / 2.

[0059] Step 2: Based on the three-dimensional atmospheric parameter field, and according to the atmospheric physical properties and the principle of electromagnetic wave propagation, calculate the atmospheric refractive index and its gradient.

[0060] The essence of ray tracing is to simulate the actual propagation path of electromagnetic waves in the atmosphere. The direction and speed of electromagnetic wave propagation are entirely determined by the atmospheric refractive index: the refractive index reflects the atmosphere's ability to refract and slow down electromagnetic waves, while the gradient reflects the spatial variation trend of this ability. Together, they constitute the core input for path simulation.

[0061] According to the formula Calculate the atmospheric refractive index n, where q is the specific humidity, corresponding to the corrected specific humidity, representing the ratio of the mass of water vapor in the air to the mass of dry air; p is the total atmospheric pressure, corresponding to the corrected atmospheric pressure, representing the pressure of the atmosphere on a unit area; e is the vapor pressure, representing the pressure generated by pure water vapor in the air; p d is the dry air pressure, representing the pressure of pure air after removing water vapor from the total air pressure; T is the atmospheric temperature, corresponding to the corrected temperature, representing the atmospheric thermal state; k1, k2, and k3 are constant coefficients.

[0062] The atmosphere is composed of dry air and moist air (i.e., water vapor). The two contribute to the refractive index through different mechanisms. The formula calculates the contributions of dry air and moist air separately, then sums them up, which reflects the actual composition of the atmosphere. The contribution of dry air is defined in the formula above. The contribution of the moist air is given by the formula above. .

[0063] Step 3: Starting from the geographical location of the target station, set the zenith angle and azimuth angle of the electromagnetic wave incident, and continuously update the propagation position and direction of the electromagnetic wave using a dynamic step size strategy based on the atmospheric refractive index and its gradient. The geographical location of the target station includes latitude, longitude and altitude.

[0064] According to the formula Update the propagation position and direction of the electromagnetic wave, where, The refractive index gradient characterizes the trend and rate of change of the refractive index in space; h m This is the normal vector for the m-th iteration, perpendicular to the current propagation direction of the electromagnetic wave, used to correct the turning direction of the path; Let r be the tangential vector of the m-th iteration, parallel to the current propagation direction of the electromagnetic wave, representing the instantaneous propagation direction; m Let be the position vector of the m-th iteration, representing the current spatial coordinates of the electromagnetic wave (composed of three-dimensional data: longitude, latitude, and altitude); Δs is the iteration step size, a dynamic value; r m+1 Let be the position vector for the (m+1)th iteration, which is the updated electromagnetic wave spatial coordinate; h is the tangential vector of the (m+1)th iteration, which is the updated electromagnetic wave propagation direction; h0 is the normal unit vector. These are tangential unit vectors, all of which are known quantities;

[0065] The turning point in the propagation path of electromagnetic waves is driven by the refractive index gradient: the larger the gradient, the more pronounced the path deflection. The formula is derived from h... m Quantify the deflection direction to ensure that the correction logic conforms to the physical essence;

[0066] Position updates need to consider both the current directional displacement and the cornering correction displacement, in r m+1 In the formula, It is the basic displacement along the current direction. It is the additional displacement caused by turning, and only by combining the two can the propagation of the curve be simulated more accurately.

[0067] The dynamic step size Δs needs to be dynamically adjusted according to the refractive index gradient. The specific formula is as follows:

[0068] , where △s min The minimum step size is Δs. max This represents the maximum step size. The reference gradient is taken as a critical value of a common gradient in the atmosphere, and α is a power coefficient used to adjust the nonlinear decay rate, with a value greater than 0.

[0069] The range of the squared cosine function in the formula is [0,1], which strictly limits the step size to [Δs]. min ,△s max This avoids the problem of the step size approaching 0 when the gradient is too large or the step size increasing infinitely when the gradient is too small.

[0070] The power function in the formula is used to enhance the overall nonlinearity. Since the dynamic range of the atmospheric refractive index gradient is extremely wide, a power function with respect to the refractive index gradient is introduced. This allows for a gradual increase in step size when the gradient is small, while the gradient approaches its maximum value when it is close to its maximum value. The nonlinear response with a rapidly decreasing time step size better reflects the actual atmospheric gradient distribution characteristics than a linear relationship.

[0071] By analyzing the actual propagation path obtained from the tracking and the propagation of electromagnetic waves in a vacuum, the atmospheric time delay is calculated. The specific operation is as follows:

[0072] The complete actual path obtained by ray tracing is divided into several continuous line segments according to the iteration step size. For each line segment, its corresponding atmospheric refractive index is extracted. The propagation time of a single line segment is calculated by dividing the line segment length by the electromagnetic wave propagation speed. Then, the propagation times of all line segments are summed to obtain the actual total propagation time.

[0073] Different line segments correspond to different refractive indices, and the propagation speed of electromagnetic waves = speed of light / refractive index. By splitting and calculating, the actual speed of each path segment can be accurately matched, avoiding errors caused by the overall average speed.

[0074] Calculate the straight-line distance between the target station and the satellite using their geographical location and spatial position; divide this straight-line distance by the speed of light in a vacuum to obtain the vacuum propagation time.

[0075] The speed of light is constant in a vacuum, and electromagnetic waves propagate in straight lines. This calculation does not depend on any atmospheric parameters, and the results obtained are relatively objective.

[0076] The atmospheric time delay is obtained by subtracting the vacuum propagation time from the calculated total actual propagation time.

[0077] The physical essence of atmospheric time delay is the extra time electromagnetic waves take to propagate in the atmosphere compared to in a vacuum; this can be directly obtained by calculating the difference between the two.

[0078] Some of the data in the above formulas are numerical calculations with dimensions removed, and the contents not described in detail in this specification are all prior art known to those skilled in the art.

[0079] The above embodiments are only used to illustrate the technical methods of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical methods of the present invention without departing from the spirit and scope of the technical methods of the present invention.

Claims

1. A method for calculating atmospheric time delay based on ray tracing, characterized in that, include: Step 1: Acquire global high spatiotemporal resolution reanalysis data and short-term high-frequency tropospheric real-time data, extract three key meteorological parameters: temperature, specific humidity, and air pressure. Addressing the vertical gaps in the data, interpolation methods are used to fill in the meteorological parameters for each layer from the ground to the top of the effective atmospheric influence layer, and grid data bias is corrected. This ultimately forms a spatiotemporally continuous and fully covered three-dimensional atmospheric parameter field. The interpolation methods include logarithmic interpolation and cubic spline interpolation: logarithmic interpolation is used for specific humidity and air pressure, while cubic spline interpolation is used for temperature. The method for correcting grid data bias is... The operation is as follows: Extract the mean value M of the same height parameter of 8 grid points in the 3×3 adjacent grid point set of the target grid point, obtain the parameter U of the adjacent height layer directly above the target grid point and the parameter L of the adjacent height layer directly below the target grid point, and calculate the correction value = (M+U+L) / 3 for the three types of parameters: temperature, air pressure and specific humidity; extract the same height parameter R of the target grid point in the global high spatiotemporal resolution reanalysis data, and check whether R is within the fluctuation range of the same height parameter in the same latitude zone of the data; extract the same height parameter Q of the target grid point in the short-term high frequency tropospheric real-time data, and check whether Q is within the nominal accuracy range of the data observation equipment; The final parameter value is determined according to the following rules: if both R and Q meet or do not meet the range requirements, the final parameter is the calculated correction value; if only R meets the range requirements, the final parameter = (calculated correction value + R) / 2; if only Q meets the range requirements, the final parameter value = (calculated correction value + Q) / 2. Step 2: Based on the three-dimensional atmospheric parameter field, and according to the atmospheric physical properties and electromagnetic wave propagation principle, calculate the atmospheric refractive index and its gradient; Step 3: Starting from the geographical location of the target station, set the zenith angle and azimuth angle of the electromagnetic wave incident. Based on the atmospheric refractive index and its gradient, continuously update the propagation position and direction of the electromagnetic wave using a dynamic step-size strategy. Analyze the actual propagation path obtained through ray tracing and the propagation of the electromagnetic wave in a vacuum to calculate the atmospheric time delay. Specifically, the complete actual path obtained through ray tracing is divided into several continuous line segments according to the iteration step size. For each line segment, its corresponding atmospheric refractive index is extracted. The propagation time of a single line segment is calculated by dividing the line segment length by the electromagnetic wave propagation speed. The propagation times of all line segments are then summed to obtain the total actual propagation time. The straight-line distance between the target station's geographical location and the satellite's spatial position is calculated. This straight-line distance is divided by the speed of light in a vacuum to obtain the vacuum propagation time. The difference between the calculated total actual propagation time and the vacuum propagation time is the atmospheric time delay result.

2. The atmospheric time delay calculation method based on ray tracing according to claim 1, characterized in that, According to the formula Calculate the atmospheric refractive index n, where q is the specific humidity, p is the total atmospheric pressure, e is the water vapor pressure, and p d Let T be the dry air pressure, T be the atmospheric temperature, and k1, k2, and k3 be constant coefficients.

3. The atmospheric time delay calculation method based on ray tracing according to claim 1, characterized in that, According to the formula Update the propagation position and direction of the electromagnetic wave, where, h is the refractive index gradient. m Let m be the normal vector of the m-th iteration. Let r be the tangential vector of the m-th iteration. m Let r be the position vector of the m-th iteration, Δs be the iteration step size, and r be the position vector of the m-th iteration. m+1 Let m+1 be the position vector of the iteration. Let h be the tangential vector of the (m+1)th iteration, and h0 be the normal unit vector. It is a tangential unit vector.

4. The atmospheric time delay calculation method based on ray tracing according to claim 3, characterized in that, According to the formula Calculate the dynamic step size Δs, where Δs min The minimum step size is Δs. max This represents the maximum step size. α is the reference gradient, and α is the power coefficient.

Images