CN116361714B - A Classification Method of Horizontal Tropospheric Delay Considering Non-isotropy - Google Patents

Abstract

The invention belongs to the technical field of global satellite navigation and positioning, and specifically discloses a horizontal tropospheric delay classification method in consideration of non-isotropy. The method comprises the steps of: estimating the tropospheric delay data at different azimuths at each altitude within the altitude range to be classified; estimating the non-isotropic values at different azimuths at each altitude according to the tropospheric delay data; Based on the medium error of the tropospheric delay, the classification thresholds at different altitude angles are constructed based on IGG-3; the size relationship between the anisotropy value and the absolute value of the classification threshold is judged, so as to realize the horizontal classification of the tropospheric delay. The present invention takes into account the non-isotropy of the tropospheric delay in the horizontal direction, realizes quantitative analysis by defining the non-isotropic value, and introduces IGG-3 to establish a classification threshold, thus well realizing the horizontal classification of the tropospheric delay .

Description

A Classification Method of Horizontal Tropospheric Delay Considering Non-isotropy

technical field

The invention belongs to the technical field of global satellite navigation and positioning, and in particular relates to a method for classifying tropospheric delay in the horizontal direction in consideration of non-isotropy.

Background technique

Tropospheric delay refers to the effect of changing the speed and path of electromagnetic wave signals when they penetrate the neutral atmosphere. It is one of the important sources of error in the high-precision navigation and positioning of the Global Navigation Satellite System (GNSS). The estimation of high-precision tropospheric delay is of great significance to precise positioning and GNSS meteorology. Generally, the tropospheric delay in the zenith direction of the reference station is called Zenith Tropospheric Delay (ZTD), and the other directions are called Slant Path Delay (SPD). In GNSS data processing, SPD is generally obtained by multiplying ZTD by a mapping function. Among them, the tropospheric delay mapping function is a key factor to improve the accuracy of GNSS tropospheric delay estimation. The function of the mapping function is to link the ZTD with the SPD in any direction to simplify the calculation.

NMF, GMF models and VMF series models are currently the most widely used and most accurate tropospheric mapping function models. The modeling of the above models is based on the principle of azimuth symmetry of meteorological elements, that is, such models are established on the basis of troposphere isotropy. However, since the atmosphere is flowing and the water vapor movement also has high temporal and spatial variation characteristics, the SPD will vary with the azimuth angle, that is, the tropospheric delay in different directions above the station will be different. It can be seen that the scheme of using the mapping function to estimate the tropospheric delay ignores the difference of the tropospheric delay in different horizontal directions.

In the high-precision positioning algorithm processing, the horizontal gradient term is generally introduced to correct the error caused by the mapping function ignoring the difference in tropospheric delay at different azimuth angles. However, the essence of horizontal gradient correction is to simply consider the atmosphere to be anisotropic, expand the first-order Fourier series of SPD in the azimuth domain, ignore higher-order terms, and then use the parameter estimation method to solve the east-west direction and north-south direction and multiply by the corresponding horizontal gradient mapping function to obtain the anisotropy value. The method of combining mapping function and horizontal gradient function to estimate tropospheric delay takes into account the difference of tropospheric delay in different horizontal directions, but it does not match the actual characteristics of tropospheric delay in the horizontal direction.

It can be seen that neither isotropy nor anisotropy can accurately represent the actual characteristics of tropospheric delay in the horizontal direction.

Before the present invention appeared, a method of estimating the tropospheric delay with high precision by using the integral method (ray tracing method) was also proposed. Based on the NWP data, this method uses the ray tracing method to estimate the SPD of all azimuth angles at a certain altitude angle at DOY196 UTC18 at the BJFS station in 2018, and makes the difference between the SPD and the average SPD corresponding to each azimuth angle at this altitude angle. Among them, the SPD values at the azimuth angles of 180-360° are similar, but the SPD values at the azimuth angles of 30-150° are significantly different. The lower the angle, the greater the difference. It can be seen from this that under the conditions of the fluidity of the atmosphere and the high temporal and spatial variation of water vapor movement, the tropospheric delay is neither isotropic nor anisotropic, but non-isotropic. Although the integration method can estimate high-precision tropospheric delay, it is difficult to apply this method to real-time positioning because of the difficulty in obtaining real-time meteorological data.

Glossary:

The isotropy of tropospheric delay means that the physical and chemical properties of objects will not change due to different measurement directions, also known as homogeneity.

The anisotropy of tropospheric delay means that all or part of the physical and chemical properties of an object change with the change of the measurement direction, showing differences in different directions.

The non-isotropy of tropospheric delay means that at a certain altitude angle, the tropospheric delay at some azimuths can be considered to be approximately equal, but at other azimuths, the tropospheric delay is quite different.

Contents of the invention

The object of the present invention is to propose a kind of horizontal direction tropospheric delay classification method that takes account of non-isotropy, quantify and express the non-isotropy of tropospheric delay by the method of proposed non-isotropic value, combine the tropospheric delay at each height angle simultaneously A classification threshold function is established based on IGG-3, so as to realize the classification of different properties (isotropic and non-isotropic) of tropospheric delay in the horizontal direction.

In order to achieve the above object, the present invention adopts the following technical solutions:

A method for classifying tropospheric delays in the horizontal direction in consideration of non-isotropy, comprising the following steps:

Step 1. Estimate the tropospheric delay data at different azimuth angles at each altitude angle within the altitude angle range to be classified;

Step 2. According to the tropospheric delay data obtained in step 1, estimate the non-isotropy value Δ at different azimuth angles at each altitude angle, where there is a corresponding anisotropy value Δ at each azimuth angle at each altitude angle;

Step 3. Combining the σ at each altitude angle, determine the threshold function k σ at different altitude angles; where σ is the median error of the tropospheric delay at the altitude angle where Δ is located, and k is the three-stage sliding window based on IGG-3 function;

Step 4. Determine the size relationship between the non-isotropic value Δ and the absolute value of the threshold function k·σ at different azimuth angles at each altitude, and classify the different properties of the tropospheric delay in the horizontal direction according to the size relationship.

The present invention has the following advantages:

As mentioned above, the present invention relates to a method for classifying tropospheric delay in the horizontal direction in consideration of non-isotropy, which breaks the original theory of isotropy and anisotropy, and takes into account the The non-isotropy in the horizontal direction, in addition to quantitative analysis by defining the non-isotropic value, also introduces IGG-3 to establish the two different properties of isotropy and non-isotropy used to distinguish tropospheric delay The classification threshold of , thereby realizing the horizontal direction classification of the tropospheric delay of the present invention. Among them, the non-isotropic value is the basis for establishing the tropospheric delay classification method in the present invention, and the classification threshold established based on IGG-3 can well ensure the continuity between the azimuth angles of the same nature after classification. The tropospheric delay classification method proposed in the present invention is conducive to establishing a tropospheric delay correction model that is closer to the real situation of the atmosphere, improves the estimation accuracy of the tropospheric delay, and reduces the influence of the tropospheric delay on high-precision GNSS positioning. The classification method proposed in the present invention will provide technical support for GNSS high-precision tropospheric delay estimation and extreme weather forecasting.

Description of drawings

FIG. 1 is a flow chart of a method for classifying tropospheric delays in the horizontal direction considering non-isotropy in an embodiment of the present invention.

Fig. 2 is a curve diagram of the value of the coefficient k in the embodiment of the present invention.

Fig. 3 is a schematic diagram of a non-isotropic classification method for tropospheric delay in an embodiment of the present invention.

Fig. 4 is the classification plan view of 5°-40° elevation angle and 10°-360° azimuth angle classification of BAKE station by using the method of the present invention.

Fig. 5 is a schematic diagram of a horizontal gradient model based on a non-isotropic classification method in an embodiment of the present invention.

Figure 6 is a schematic diagram of RT-SPD, MF-SPD with horizontal gradient correction, and SPD estimated based on the horizontal gradient model estimation of the classification method of the present invention at 17.5° elevation angle SPD at BAKE station 2019 DOY001 UTC18 time.

Detailed ways

Below in conjunction with accompanying drawing and specific embodiment the present invention is described in further detail:

As shown in Figure 1, the present embodiment describes a method for classifying tropospheric delay in the horizontal direction considering non-isotropy, which includes the following steps:

Step 1. Estimate the tropospheric delay data of different azimuth angles at each altitude angle (within the range of 0°-360°) within the altitude angle range to be classified. The needs are self-determined.

In this embodiment, for example, the ray tracing method is used to estimate the tropospheric delay data using meteorological parameters.

Step 2. Based on the tropospheric delay data obtained in step 1, estimate the anisotropy value Δ at different azimuths at each altitude, where there is a corresponding Δ at each azimuth at each altitude.

The method of determining the non-isotropic value Δ is: make the difference between the SPD of each azimuth angle at a certain elevation angle and the meanspd of the SPD of each azimuth angle at the elevation angle, and define the difference as non-isotropic value Δ.

The non-isotropic value Δ is due to the asymmetry, and its expression is shown in Equation (1):

Δ=SPD - mean spd (1)

where SPD represents the tropospheric delay, and meanspd is the average value of the SPD at each azimuth at the altitude angle to be classified, that is, the isotropic value of the tropospheric delay due to the symmetry of the atmosphere.

The anisotropy value Δ is a quantitative representation of the anisotropy of the tropospheric delay. Moreover, after the Lilliefors hypothesis test, the non-isotropic value at a certain height angle is in line with the normal distribution.

The non-isotropic value Δ inherits some of the characteristics of SPD, such as Δ increases exponentially with the decrease of the altitude angle. According to the experiment, the value of Δ is about 0.1mm when the altitude angle is 40°, and at 5° The maximum height can reach dm level.

Since the non-isotropic value Δ defined in this embodiment removes the isotropic part, the non-isotropic value Δ can better represent the change characteristics of the horizontal direction of the SPD at a certain elevation angle.

In order to realize the present invention, in addition to realizing quantitative analysis by defining the non-isotropy value Δ, another key issue is how to determine a threshold to distinguish the two different properties of isotropy and non-isotropy of tropospheric delay.

Step 3. Combining the σ at each altitude angle, determine the threshold function k σ at different altitude angles; where σ is the median error of the tropospheric delay at the altitude angle where Δ is located, and k is the three-stage sliding window based on IGG-3 function.

The present invention will use a function with continuity as the threshold value, and the threshold value in the form of a function with continuity can effectively detect abnormal values in the data, and avoid a certain azimuth in the continuous range of azimuth angles showing the same nature of the tropospheric delay. When the tropospheric delay at the angle exhibits other properties, the azimuths where the classified tropospheric delays with the same properties are kept continuous, which is conducive to establishing a tropospheric delay estimation model that is closer to the real situation of the atmosphere.

For this reason, the present invention uses the classification method of the IGG-3 weight scheme to construct a three-stage sliding window function k as a part of the threshold. The IGG-3 rights scheme divides data into three categories: effective information, usable information, and harmful information through data quality, and makes full use of effective information, limits the impact of available information, and excludes harmful information. It is a scheme suitable for measurement data processing.

A three-stage sliding window function k is constructed based on IGG-3, and its expression is shown in formula (2).

(2)

Among them, |u| is the standardized residual, , e ₀ is the first limit threshold, and e ₁ is the second limit threshold.

According to a large number of repeated experiments, it is shown that within the range of 40°-5° altitude angle, the value range of |u| is about 0.3~3.8, so in this embodiment, the first limit threshold e ₀ and the second limit threshold e ₁ They are respectively set as: e ₀ =0.5, e ₁ =2.0.

The graph of the value range of k is shown in Figure 2. In order to make the classification results close to the actual characteristics of the tropospheric delay and facilitate the establishment of a more accurate tropospheric delay estimation model by using the classification results, the maximum value of k is set to 0.8, the minimum value is set to 0.7, and the rest are set to It has a continuous function form, so it can make full use of effective information, limit available information, and exclude harmful information, so as to realize the continuity between azimuth angles of the same nature after classification.

Under the conditions of a certain spatial position, time, and altitude angle, the k times of the medium error σ of the SPD at the altitude angle is used as the threshold value, and the tropospheric delay non-isotropy value is used for quantification, and combined with the non-isotropic value of the tropospheric delay same-sex classification.

Step 4. Determine the size relationship between the non-isotropic value Δ and the absolute value of the threshold function k·σ at different azimuth angles at each altitude, and classify the different properties of the tropospheric delay in the horizontal direction according to the size relationship.

When Δ>|k·σ|, the tropospheric delay at the azimuth angle corresponding to Δ that satisfies this condition is positively anisotropic; The tropospheric delay at is reverse anisotropy.

When -|k·σ|≤Δ≤|k·σ|, the tropospheric delay at the azimuth corresponding to Δ that satisfies this condition is isotropic.

As shown in Figure 3, the vertical axis represents different properties, the horizontal axis represents the non-isotropic value, the dotted gray stripe area represents the forward part, the solid gray stripe area represents the reverse part, and the solid black stripe area in the center represents Same-sex section.

At a certain moment, for any station, the present invention can classify the tropospheric delay at each azimuth at different elevation angles of the station, so as to provide technical support for establishing a more accurate tropospheric delay estimation model that is more in line with the actual characteristics of the atmosphere.

The following uses the SPD of the BAKE station in the IGS station at DOY00158 UTC18 in 2020, with an elevation angle of 5°–40° and an azimuth angle of 10°–360° as an example data, and is classified according to the above classification method.

Fig. 4 shows the result of classifying example data according to the classifying method and steps of the present invention. Among them, the polar coordinate system takes 0 m as the center, its polar axis represents SPD (in m), and the polar angle represents the azimuth angle.

The different circular curves in the figure represent the SPD at different altitude angles, the innermost curve is the SPD at each azimuth angle at an altitude angle of 40°, and the outermost circle is the SPD at each azimuth angle at an altitude angle of 5°. The dotted portion of each curve represents the inverse SPD, the dotted portion with a circle mark represents the orthotropic SPD, and the solid line represents the isotropic SPD.

It can be seen from Fig. 4 that the SPD between the 150° azimuth angle and the 270° azimuth angle at each altitude angle of the BAKE station is positive, while the SPD is reversed between the 330° and 90° azimuth angle, and the other azimuth angles are negative. The SPD presents the same sex.

After applying the classification method proposed in the present invention to complete the classification of SPD, targeted modeling can be aimed at SPDs of different properties (forwardness, homosexuality, reverseness), so that the established model takes into account the non-variety of SPD in the horizontal direction. Isotropy, rather than building a simple model based on isotropy or anisotropy, is conducive to estimating SPDs that are closer to the real characteristics of tropospheric delay and have higher accuracy, and then in the fields of GNSS high-precision positioning, GNSS meteorology, etc. provide technical support.

In addition, in order to verify the effectiveness of the classification method proposed in the present invention, based on the classification method of the present invention, combined with the horizontal gradient correction model, a non-isotropic horizontal gradient correction method based on tropospheric delay is established.

As shown in Figure 5, the non-isotropic horizontal gradient correction method based on tropospheric delay includes the following steps:

Step 1. Make a difference between the MF–SPD with horizontal gradient correction and RT–SPD between the same azimuth angles at each altitude angle, and determine the maximum difference between the two at each altitude angle (hereinafter referred to as the maximum difference) Whether it is greater than 0.

For the convenience of description, the RT-SPD scheme is taken as scheme 1, and the MF-SPD scheme is taken as scheme 2.

The mapping function used in Scheme 2 is VMF1.

Second step. The data of scheme 1 is classified according to the proposed classification method of the present invention, wherein the azimuth of the positive SPD is called the forward part, and the azimuth of the reverse SPD is called the reverse sex part.

Step 3. If the maximum difference at a certain altitude is greater than 0, it means that the estimated value of scheme 2 is greater than the estimated value of scheme 1, that is, in the positive part, the horizontal gradient leads to the difference between MF–SPD and RT–SPD The difference increases. In this case, only VMF1 is used to estimate the SPD in the directivity part, and the SPD in the remaining azimuths will be estimated according to scheme 2.

If the maximum difference at a certain altitude is less than 0, it means that the estimated value of scheme 2 is smaller than that of scheme 1, that is, in the inverse part, the horizontal gradient causes the difference between MF–SPD and RT–SPD to increase, and in In this case, only VMF1 is used to estimate the SPD of the inverse part, while the other azimuths will be estimated according to scheme 2.

In addition, the boundary between SPDs of different properties will be smoothed.

Taking the SPD of 5-40° elevation angle and 10-360° azimuth angle at BAKE station 2019 DOY001 UTC18 as a calculation example, construct scheme 3 (SPD estimated based on the horizontal gradient model of the classification method of the present invention), and combine the implementation steps of this model Option 1 and Option 2 mentioned in , verify the validity of the model and evaluate the accuracy of the model, and Option 1 will be used as a reference.

Fig. 6 shows the estimated SPD (17.5° elevation angle) of the BAKE station according to the three schemes. In the figure, the polar axis of the polar coordinate system represents the SPD value (in m, 7.590 m at the pole and 7.620 m at the outermost boundary), and the polar angle represents the azimuth.

Combining three kinds of points with different densities into a graph is the result of scheme 1, and the density of the points from dense to sparse indicates reverse, homosexuality, and forwardness respectively; the oblique line graph represents the result of scheme 2; the graph filled with black represents the scheme 3 results.

For the altitude angle of 17.5°, the maximum difference between scheme 2 and scheme 1 is less than 0, so the reverse part (10°–120° azimuth) is only used to estimate SPD using VMF1, and the remaining azimuth angles are estimated according to scheme 2.

It can be seen from Figure 6 that, compared with Scheme 2, Scheme 3 is obviously closer to Scheme 1, and its estimated value is more accurate.

Known by quantitative analysis, when using scheme 1 as a reference, the RMS of scheme 2 is 6.1mm, and the RMS of scheme 3 is 3.9mm, it can be seen that scheme 3 based on the proposed classification method of the present invention has improved about 36% compared with scheme 2 precision. %.

This experiment proves that the non-isotropic classification method of tropospheric delay proposed by the present invention is effective.

The invention provides new modeling ideas and technical support for establishing a high-precision tropospheric delay correction model that is closer to the real situation of the atmosphere, and has important research significance and application value in GNSS meteorology, GNSS high-precision positioning, and the like.

Of course, the above descriptions are only preferred embodiments of the present invention, and the present invention is not limited to the above-mentioned embodiments. It should be noted that all equivalent substitutions made by any person skilled in the art under the teaching of this specification , obvious deformation forms, all fall within the essential scope of this specification, and should be protected by the present invention.

Claims (2)

1. A method for classifying retardation of a troposphere in a horizontal direction considering non-isotropy is characterized in that,

the method comprises the following steps:

step 1, estimating troposphere delay data of different azimuth angles at each altitude angle in an altitude angle range to be classified;

step 2, estimating non-isotropy values delta of different azimuth angles at each altitude angle according to the troposphere delay data obtained in the step 1; wherein a corresponding non-isotropic value delta exists at each azimuth at each altitude;

in the step 2, the expression of the non-isotropy value delta is as shown in formula (1):

Δ=SPD﹣meanspd （1）

wherein SPDs represent tropospheric delay, and meanspd is the average value of SPDs of azimuth angles at altitude to be classified;

step 3, combining sigma at each height angle to determine a threshold function k.sigma at different height angles; wherein sigma is a medium error of tropospheric delay at a height angle where delta is located, and k is a three-section sliding window function constructed based on IGG-3;

in the step 3, the expression of the three-section sliding window function k is shown in formula (2);

（2）

where |u| is the normalized residual,，e ₀ for a first limiting threshold, e ₁ Is a second limiting threshold;

step 4, judging the magnitude relation between the non-isotropy value delta of each azimuth angle at each altitude angle and the absolute value of the threshold function k.sigma, and realizing classification of different properties of tropospheric delay in the horizontal direction according to the magnitude relation;

in the step 4, the classification process of the different properties of the tropospheric delay in the horizontal direction is as follows:

when delta > |k.sigma|, the troposphere delay at the azimuth angle corresponding to delta meeting the condition is positive; when delta is minus k-sigma, the tropospheric delay at the azimuth angle corresponding to delta meeting the condition is reverse;

when the absolute value of the k and the absolute value of the sigma is less than or equal to the absolute value of the delta, the troposphere delay at the azimuth angle corresponding to the delta meeting the condition shows the same polarity.

2. The method of classifying a horizontal tropospheric delay in consideration of non-isotropy as claimed in claim 1,

in the step 1, tropospheric delay data is estimated by using meteorological parameters by a ray tracing method.