Method for generating and setting weights of regional enhancement observation based on gim model

Abstract

The application provides a regional enhancement observation generation and weight setting method based on a GIM model, comprising the following steps: generating regional constraint points and extracting VTEC values from IONEX data; combining latitude, distance and density factor to calculate the geographical distribution weight of GNSS stations; assigning the weight of the nearest station to the constraint points and calculating the maximum weight; constructing an enhanced design matrix containing the constraint points, and solving spherical harmonic coefficients by using weighted least squares. The application has the beneficial effects of improving the VTEC estimation accuracy in the sparse area of GNSS stations, reducing the condition number of the design matrix to enhance the stability of the model, adaptively adjusting the weight, supporting the configuration of any regional enhancement, and being efficient in calculation.

Description

Technical Field

[0001] This invention belongs to the field of GNSS ionospheric modeling and regional enhancement technology, and in particular relates to a method for generating and weighting regional enhancement observations based on the GIM model. Background Technology

[0002] In GNSS ionospheric modeling, spherical harmonic function expansion is commonly used to parametrically model the global ionospheric VTEC. The spherical harmonic function model represents VTEC as: Where t represents epoch time; Indicates latitude and longitude; VTEC This indicates that the latitude and longitude coordinates at time t are ( The total number of vertical electrons at position ( ); N represents the order of the spherical harmonic function; n and m represent the order and degree of the spherical harmonic function, respectively; , , These are the coefficients, cosine coefficients, and sine coefficients of the nth-order m-degree associated Legendre polynomial. The solution is obtained using weighted least squares estimation. , A spherical harmonic function model of the total vertical electron volume of the globe is obtained. However, due to the uneven geographical distribution of GNSS stations, the observation data in some areas such as the ocean, polar regions, and remote areas is sparse, which leads to an increase in the condition number of the design matrix A in sparse regions. The global ionospheric model based on the spherical harmonic function is not accurate enough in these regions. The existing technology has the following problems: (1) The selection of regional observation points lacks specificity and does not consider the distribution characteristics of GNSS stations; (2) The weight allocation mechanism is simple, only considering the distance factor and ignoring the latitude and density factors; (3) It cannot effectively balance the observation contributions of different regions, resulting in a decrease in accuracy in sparse regions. Summary of the Invention

[0003] In view of this, the present invention aims to propose a method for generating and weighting regional enhanced observations based on the GIM model, so as to solve at least one of the problems existing in the above-mentioned prior art.

[0004] To achieve the above objectives, the technical solution of the present invention is implemented as follows: The method for generating and weighting regional augmented observations based on the GIM model includes the following steps: S1. Based on the IONEX format GIM data file, parse the file header information to obtain the mesh parameters; S2. Determine the geographic boundary of the target augmentation area based on grid parameters; S3. Within the target area, randomly generate a set number of constraint points and obtain the positions of the constraint points; S4. For each constraint point, use interpolation to extract the VTEC value from the IONEX data and store the constraint point information as a constraint point set. S5. Based on the global information set of GNSS stations, calculate the geographical distribution weight of each GNSS station and generate a set of station weights. S6. Based on the set of constraint points and the set of station weights, for each constraint point, find the station closest to the constraint point, assign the geographical distribution weight of the station to the constraint point, and then combine it with the global constraint weight factor to calculate the constraint point weight. S7. Construct the design matrix: For each GNSS observation point, occultation observation point, and constraint point, calculate the row vector of the spherical harmonic function basis function and fill the row vector into the design matrix; S8. Based on the VTEC values ​​of each GNSS observation point, occultation observation point, and constraint point, construct an observation vector; based on the geographical distribution weights of each GNSS observation point, occultation observation point, and constraint point, construct a weight matrix. S9. Based on the design matrix, weight matrix, and observation vector, construct a weighted least squares objective function; S10. Based on the objective function, obtain and solve the normal equation, and then solve for the spherical harmonic coefficients.

[0005] Furthermore, in step S3, the expression for the constraint point location is as follows: ; ; In the formula, , Let the latitude and longitude of the i-th generated constraint point be... and In the interval Random numbers that are uniformly distributed within the range.

[0006] Furthermore, in step S5, the expression for calculating the geographic distribution weights is as follows: ; In the formula, For geographical distribution weights, , , These are the importance coefficients for latitude, distance, and density weighting factors, respectively. Latitude weighting factor Distance weighting factor Density weighting factor; ; In the formula, The latitude of the measuring station; ; In the formula, For the station The minimum great circle distance to other stations. This is the distance attenuation scale parameter; ; In the formula, For the purpose of measuring station Center, search radius The number of stations within the area This is the density threshold parameter.

[0007] Furthermore, in step S6, the expression for the constraint point weights is as follows: ; ; ; In the formula, The index of the station closest to constraint point i. Let i be the distance from point i to station j. This represents the total number of stations. For the nearest neighbor station j ∗ Geographic weight, For the constraint point weights, The geographical weights inherited by constraint point i, As a global constraint weighting factor, Distance constraint point Recent station index Surrounding search radius The number of stations within the area This is the density threshold parameter.

[0008] Furthermore, in step S7, for each GNSS observation point, the expression for the row vector of the spherical harmonic function basis function is calculated as follows: ; In the formula, To design the k-th row vector of matrix A, Let n be a Legendre polynomial, where n is the order and m is the degree. , Let be the latitude and longitude of the k-th GNSS observation point.

[0009] Furthermore, in step S9, the expression for the weighted least squares objective function is as follows: ; ; In the formula, This represents the total number of observations. This refers to the total number of GNSS vertical electron observations. This refers to the total number of observed vertical electrons during occultation. The number of GIM constraint points. The objective function value, Let be the vector of spherical harmonic coefficients to be solved. To design the matrix, For the observation vector, This is the weight matrix. denoted as the total number of spherical harmonic coefficients, and N as the maximum order of the spherical harmonic function.

[0010] Furthermore, in step S10, the expression for the normal equation is as follows: ; In the formula, Let A be the transpose of matrix A.

[0011] Compared with existing technologies, the region augmentation observation generation and weight setting method based on the GIM model described in this invention has the following advantages: This invention improves regional accuracy, increasing VTEC estimation accuracy by 20-40% in sparse GNSS station areas, particularly in marine, polar, and remote regions. It also enhances model stability by reducing the condition number of the design matrix, improving numerical stability, and decreasing the variance of spherical harmonic coefficient estimation. Furthermore, it improves adaptability by automatically adjusting constraint point weights based on station distribution without manual intervention, adapting to different regional characteristics. It offers regional flexibility, supporting enhancement configurations for any geographical area, allowing selection of enhancement regions based on actual needs. Finally, it improves computational efficiency, enabling multiple reuse of pre-calculated weights, reducing redundant computation overhead, and supporting parallel processing. Attached Figure Description

[0012] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings: Figure 1 This is a schematic diagram of the overall process described in an embodiment of the present invention; Figure 2 This is a schematic diagram of the geographic distribution weight calculation process described in an embodiment of the present invention; Figure 3 This is a schematic diagram of the constraint point weight allocation process according to an embodiment of the present invention. Detailed Implementation

[0013] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0014] In the description of this invention, it should be understood that the terms "center," "longitudinal," "lateral," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicating orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. Furthermore, the terms "first," "second," etc., are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, a feature defined with "first," "second," etc., may explicitly or implicitly include one or more of that feature. In the description of this invention, unless otherwise stated, "a plurality of" means two or more.

[0015] In the description of this invention, it should be noted that, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal connection of two components. Those skilled in the art will understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0016] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0017] like Figures 1 to 3 As shown, the method for generating and setting weights for region-enhanced observations based on the GIM model includes the following specific steps: S1. Based on the IONEX format GIM data file, parse the file header information to obtain the mesh parameters; S2. Determine the geographic boundary of the target augmentation area based on grid parameters; S3. Within the target area, randomly generate a set number of constraint points and obtain the positions of the constraint points; S4. For each constraint point, use interpolation to extract the VTEC value from the IONEX data and store the constraint point information as a constraint point set. S5. Based on the global information set of GNSS stations, calculate the geographical distribution weight of each GNSS station and generate a set of station weights. S6. Based on the set of constraint points and the set of station weights, for each constraint point, find the station closest to the constraint point, assign the geographical distribution weight of the station to the constraint point, and then combine it with the global constraint weight factor to calculate the constraint point weight. S7. Construct the design matrix: For each GNSS observation point, occultation observation point, and constraint point, calculate the row vector of the spherical harmonic function basis function and fill the row vector into the design matrix; S8. Based on the VTEC values ​​of each GNSS observation point, occultation observation point, and constraint point, construct an observation vector; based on the geographical distribution weights of each GNSS observation point, occultation observation point, and constraint point, construct a weight matrix. S9. Based on the design matrix, weight matrix, and observation vector, construct a weighted least squares objective function; S10. Based on the objective function, obtain and solve the normal equation, and then solve for the spherical harmonic coefficients.

[0018] In a preferred embodiment of the present invention, in step S1, the header information of the GIM data file in IONEX format is parsed to obtain the mesh parameters.

[0019] Specifically, the IONEX (IONosphere map exchange) format is an ionospheric data exchange format released by the International GNSS Service (IGS). GIM (Global Ionosphere Map) data files contain gridded data of the global ionosphere's VTEC (Vertical Total Electron Content). The file header information includes key parameters such as the grid's latitude and longitude range, time resolution, altitude information, and spherical harmonic order.

[0020] In this embodiment, the following mesh parameters are obtained by parsing the IONEX file header: Latitude range: to Grid spacing Longitude range: to Grid spacing Time resolution: Typically 2 hours; Altitude information: Height of the ionospheric thin layer, typically 450 km; Maximum order of the spherical harmonic function. It is usually of order 15.

[0021] In a preferred embodiment of the present invention, in step S2, the geographical boundary of the target enhancement area is determined based on the grid parameters.

[0022] The target enhancement region refers to the specific geographic area where VTEC estimation accuracy needs to be improved. Based on the grid parameters obtained in step S1, the geographic boundaries of the target region are determined, including: Minimum latitude and maximum latitude Minimum longitude and maximum longitude .

[0023] The target enhancement area can be of any shape and size, such as ocean areas, polar areas, remote areas, and other areas with sparse GNSS stations.

[0024] In a preferred embodiment of the present invention, in step S3, a set number of constraint points are randomly generated within the target area, and the positions of the constraint points are obtained.

[0025] Constraint points are virtual observation points used to enhance the amount of observation information in the target area. In this embodiment, the expression for the location of the constraint point is as follows: ; ; In the formula, , For the first The latitude and longitude of each generated constraint point. and In the interval Random numbers that are uniformly distributed within the range.

[0026] Number of constraint points This is a configuration parameter; the default value is 5183. The number of constraint points can be adjusted based on the size of the target region and the desired level of enhancement.

[0027] In a preferred embodiment of the present invention, in step S4, for each constraint point, the VTEC value is extracted from the IONEX data using interpolation, and the constraint point information is stored as a constraint point set.

[0028] For each constraint point Extract a specified time from IONEX data VTEC value The specific interpolation method is as follows: 1. Time interpolation: If IONEX data is interpolated in time... If no exact match exists, linear interpolation along the time dimension is performed. Let... and To and Two adjacent IONEX data moments ( ),but: ; 2. Spatial interpolation: If at position If no exact match exists, bilinear interpolation is used to obtain the match from surrounding grid points. Let... Located at four grid points , , , Between, then: ; in, , , , These are the VTEC values ​​for the four grid points.

[0029] If the extracted VTEC value is less than 0 or is NaN (Not a Number), then skip the constraint point and do not include it in the constraint point set.

[0030] Store valid constraint point information as a set of constraint points: ; in, This represents the number of valid constraint points.

[0031] In a preferred embodiment of the present invention, in step S5, based on the global information set of GNSS stations, the geographical distribution weight of each GNSS station is calculated, and a station weight set is generated.

[0032] This step corresponds to the technical solution described in claim 3. For example... Figure 2 As shown, the calculation process of geographic distribution weights includes: obtaining information from all GNSS stations, traversing each station, calculating latitude weight, distance weight, and density weight respectively, combining the weights and normalizing them, and finally storing the station weights.

[0033] Specifically, the first step is to obtain a global information set from all GNSS stations: ;in, For the number of stations, For the first Latitude and longitude of each station, in radians.

[0034] For each GNSS station Calculate its geographical distribution weight This weight is composed of three factors: (1) Latitude weighting factor Considering the weaker ionospheric activity in high-latitude regions, the square of the cosine function is used as the latitude weight, as shown in the following expression: This weight is at the equator ( The value is 1 at point ) and in the polar regions ( The value at point ) is 0. This design reflects the physical property that the electron density of the ionosphere varies with latitude.

[0035] (2) Distance weighting factor Calculation of stations Minimum great circle distance to all other stations The great circle distance is calculated using the Haversine formula, as shown below: In the formula, The average radius of the Earth , .

[0036] The distance weights use a Gaussian function, expressed as follows: In the formula, For the station The minimum great circle distance to other stations. This is the distance attenuation scale parameter, with a default value of 500 km.

[0037] (3) Density weighting factor Statistics based on monitoring stations Center, search radius Number of stations within (default value is 1000 km) The density weights use the Sigmoid function, expressed as follows: In the formula, For the purpose of measuring station Center, search radius The number of stations within the area This is the density threshold parameter, with a default value of 5 stations.

[0038] (4) Combination of geographical distribution weights: The three weighting factors are combined by weighted average to obtain the station The geographical distribution weights are expressed as follows: In the formula, , , These are the importance coefficients for latitude, distance, and density weighting factors, with default values ​​of 1.0, 1.0, and 0.5 respectively. The normalized weights... Generate a set of station weights: That is, the mapping relationship between the location of the monitoring station and the weight.

[0039] In a preferred embodiment of the present invention, in step S6, based on the set of constraint points and the set of station weights, for each constraint point, the station closest to the constraint point is found, and the geographical distribution weight of the station is assigned to the constraint point. Then, the constraint point weight is calculated by combining the global constraint weight factor.

[0040] like Figure 3 As shown, the constraint point weight allocation process includes: obtaining the constraint point location, initializing the weight to 1.0, traversing the station weight list, calculating the great circle distance, determining if the distance is smaller, updating the minimum distance and weight, applying the global weight factor, and obtaining the final weight. Specifically: (1) Find the nearest station: For each constraint point In the set of station weights Find the station closest to the constraint point. Calculate the great circle distances from the constraint point to all stations using the Haversine formula. Find the station corresponding to the minimum distance, as shown in the following expression: In the formula, Distance constraint point The latest station index, Constraint points to the station distance, This represents the total number of monitoring stations.

[0041] (2) Inherit geographical weights: Assign the geographical distribution weights of the nearest station to the constraint points, as shown in the following expression: In the formula, Distance constraint point Recent station index Surrounding search radius The number of stations within the area This is the density threshold parameter. Constraint points Inherited geographical weights, the physical meaning of which is in areas with sparse monitoring stations, The exponent term is negative. Approaching 1, constraint points are assigned high weights to enhance observational contributions; in densely populated areas of stations, The exponent term is positive. Approaching 0, the constraint points are assigned low weights to avoid redundant observations overwhelming the true observations.

[0042] (3) Calculate the final weights: combine the global constraint weight factors Configuration parameter value range The default value is 2.0, which is used to calculate the final weights of the constraint points. The expression is as follows: Weight The role of the weight matrix used in weighted least squares estimation and the global constraint weight factors: This factor is used to improve the overall contribution of constraint points in weighted least squares estimation. Combined with the design of high weights for sparse regions in geographical weights, this factor can ensure that constraint points play a significant enhancing role in sparse regions, effectively reducing the condition number of the design matrix and improving model stability.

[0043] In a preferred embodiment of the present invention, in step S7, a design matrix is ​​constructed: for each GNSS observation point, occultation observation point, and constraint point, the row vector of the spherical harmonic function basis function is calculated, and the row vector is filled into the design matrix.

[0044] The spherical harmonic function model represents VTEC as: ; in, For Legendre polynomials, For order, For the number of times, and is the spherical harmonic coefficient.

[0045] For each GNSS observation point Calculate its spherical harmonic basis function row vector, as shown in the following expression: In the formula, For designing a matrix The Row vectors.

[0046] Similarly, for occultation observation points and GIM constraint points, the row vectors of the spherical harmonic basis functions are calculated using the same method. All row vectors are then filled into the design matrix. ,in The total number of spherical harmonic coefficients. is the maximum order of the spherical harmonic function.

[0047] In a preferred embodiment of the present invention, in step S8, an observation vector is constructed based on the VTEC value of each GNSS observation point, occultation observation point, and constraint point, and a weight matrix is ​​constructed based on the geographical distribution weight of each GNSS observation point, occultation observation point, and constraint point.

[0048] Construct observation vectors For GNSS observation points, the observed value is the VTEC value; for occultation observation points, the observed value is the VTEC value corresponding to the peak electron density location; for GIM-constrained points, the observed value is... .

[0049] Constructing the weight matrix For GNSS and occultation observation points, the weights are determined based on the observation accuracy; for GIM-constrained points, the weights are... .

[0050] In a preferred embodiment of the present invention, in step S9, a weighted least squares objective function is constructed based on the design matrix, the weight matrix, and the observation vector.

[0051] First, count the number of observation data: ; In the formula, This refers to the total number of GNSS vertical electron observations. This refers to the total number of observed vertical electrons during occultation. This represents the number of GIM constraint points.

[0052] Construct the weighted least squares objective function: ; In the formula, Let be the vector of spherical harmonic coefficients to be solved. To design the matrix, For the observation vector, This is the weight matrix. It is the square root of the weight matrix.

[0053] In a preferred embodiment of the present invention, in step S10, the normal equation is obtained and solved based on the objective function, and then the spherical harmonic coefficients are obtained.

[0054] For the objective function Differentiating and setting it equal to zero, we obtain the normal equation: In the formula, For matrix The transpose of .

[0055] Use numerically stable solution methods, such as SVD decomposition or Cholesky decomposition, to solve for the spherical harmonic coefficients. .

[0056] Output estimated spherical harmonic coefficients This is used for subsequent VTEC prediction and ionospheric modeling.

[0057] Example 2: Detailed Implementation of Constraint Point Generation This embodiment details the method for generating the constraint point positions in step S3.

[0058] Constraint points are randomly generated virtual observation points within the target enhancement region, used to increase the amount of observational information in sparse areas. The location of the constraint points is determined by the following formula: ; ; in: and The minimum and maximum latitudes of the target area; and The minimum and maximum longitudes of the target area; and In the interval A random number that is uniformly distributed within the range; and For the first The latitude and longitude of each constraint point. The number of constraint points. This is a configuration parameter with a default value of 5183. This value can be adjusted based on the size of the target area and the desired level of enhancement. For example, for a smaller area such as 10°×10°, fewer constraint points, such as 500, can be used; for a larger area such as 60°×60°, more constraint points, such as 5000, can be used.

[0059] The process of generating constraint points is as follows: 1. Determine the geographical boundaries of the target area. ; 2. Generate For random numbers ,in ; 3. Calculate the position of each constraint point according to the above formula. ; 4. Extract the VTEC value for each constraint point from the IONEX data; 5. Filter out constraint points with invalid VTEC values, values ​​less than 0, or values ​​of NaN to obtain the set of valid constraint points. .

[0060] Example 3: Detailed Implementation of Geographic Distribution Weight Calculation This embodiment details the method for calculating the geographical distribution weight in step S5.

[0061] The geographic distribution weight reflects the distribution characteristics of GNSS stations in geographic space and is composed of three factors: latitude weight factor, distance weight factor, and density weight factor.

[0062] 3.1 Latitude Weighting Factor: Latitude Weighting Factor Considering the physical characteristics of weaker ionospheric activity in high-latitude regions, the square of the cosine function is used as the latitude weight: ; in, For the station Latitude, unit: radians.

[0063] The characteristic of this weight: at the equator ( ) place, The weight is the largest; in the polar regions ( ) place, In the equatorial region, the weight is the smallest; in the mid-latitude region, the weight is between 0 and 1. This design reflects the physical characteristics of ionospheric electron density varying with latitude: ionospheric activity is strongest in the equatorial region and weakest in the polar region.

[0064] 3.2 Distance Weighting Factor: Distance Weighting Factor This reflects the spatial distribution density between stations, calculated using a Gaussian function: ;in: For the station The minimum great circle distance to other stations; This is the distance attenuation scale parameter, with a default value of 500 km.

[0065] The great circle distance is calculated using the Haversine formula: ;in: The average radius of the Earth; It is the difference in latitude; This is due to the difference in longitude.

[0066] Characteristics of distance weighting factors: When hour, The one with the highest weight; when hour, ;when hour, The weight of a single station is the lowest. This design reflects the local density characteristics of the station distribution: the denser the stations, the lower the weight of a single station; the sparser the stations, the higher the weight of a single station.

[0067] 3.3 Density Weighting Factor: Density Weighting Factor This reflects the number of stations within a certain range around the current station, calculated using the Sigmoid function: ;in: For the purpose of measuring station Center, search radius The default value is the number of stations within 1000 km; This is the density threshold parameter, with a default value of 5 stations.

[0068] Characteristics of density weighting factor: When hour, The weight is the smallest; when hour, ;when hour, The higher the station density, the higher the weight; the lower the station density, the lower the weight. This design reflects the regional density characteristics of the station distribution: the higher the station density, the higher the weight; the lower the station density, the lower the weight.

[0069] 3.4 Combination of Geographical Distribution Weights: The three weighting factors are combined by weighted average to obtain the station's weighted average. Geographical distribution weights: ;in, , , These are the importance coefficients for latitude, distance, and density weighting factors, respectively, with default values ​​of 1.0, 1.0, and 0.5.

[0070] The principle for setting importance coefficients: Latitude and distance weights have higher importance coefficients, which can be 1.0, because they directly reflect the physical characteristics of the ionosphere and the spatial distribution of stations; density weights have lower importance coefficients, which can be 0.5, because they are a supplementary description of the station distribution density. Normalized weights. .

[0071] Example 4: Detailed Implementation of Constraint Point Weight Allocation This embodiment details the method for calculating the constraint point weights in step S6.

[0072] The constraint point weight reflects the importance of the constraint point in the target area. It is calculated by finding the station closest to the constraint point and inheriting the geographical distribution weight of that station.

[0073] 4.1 Find the nearest station: For each constraint point In the set of station weights Find the station closest to the constraint point. Calculate the great circle distances from the constraint point to all stations using the Haversine formula. Find the station corresponding to the minimum distance: ;in: Distance constraint point Recent station index; Constraint points to the station The distance; This represents the total number of monitoring stations.

[0074] 4.2 Inherit Geographic Weights: Assign the geographic distribution weights of the nearest monitoring station to the constraint points: ; in: Nearest neighbor station Geographic weight; Constraint points Inherited geographical weights.

[0075] This design allows the weight of constraint points to reflect the distribution characteristics of stations in their respective areas: in densely populated areas, the weight of constraint points is higher; in sparsely populated areas, the weight of constraint points is lower.

[0076] 4.3 Calculate the final weights: Combine global constraint weight factors Calculate the final weights of the constraint points: ;in, This is the global constraint weight factor, and the default value of the configuration parameter is 1.0.

[0077] The role of global constraint weighting factors: when When the constraint point's weight is equal to its inherited geographical weight, when... When, enhance the effect of constraint points; when At this time, the effect of the constraint point is weakened. Weight The weight matrix used in weighted least squares estimation.

[0078] Example 5: Detailed Implementation of Row Vectors for Spherical Harmonic Function Basis Functions This embodiment details the method for calculating the row vectors of the spherical harmonic function basis functions in step S7.

[0079] The spherical harmonic function model represents VTEC as: ; in: For Legendre polynomials; For order, For the number of times; and These are the spherical harmonic coefficients; This is the maximum order of the spherical harmonic function, usually 15.

[0080] For each GNSS observation point Calculate its spherical harmonic basis function row vectors: ; in: For designing a matrix The Row vectors; , For the first Latitude and longitude of each GNSS observation point.

[0081] The dimension of a row vector is This corresponds to the number of all spherical harmonic coefficients. For example, when hour, .

[0082] Legendre polynomials The calculation uses a recursive formula: ; ; Similarly, for occultation observation points and GIM constraint points, the spherical harmonic basis function row vectors are calculated using the same method and then filled into the design matrix. middle.

[0083] Example 6: Detailed Implementation of the Weighted Least Squares Objective Function This embodiment details the method for constructing the weighted least squares objective function in step S9.

[0084] 6.1 Statistics of the total number of observations: First, the number of observation data points is counted: ;in: This refers to the total number of GNSS vertical electron observations. The number of observations of the total vertical electron count of the occultation; This represents the number of GIM constraint points.

[0085] 6.2 Weighted Least Squares Objective Function: Constructing the weighted least squares objective function: ;in: Let be the vector of spherical harmonic coefficients to be solved; For designing the matrix; For observation vectors; This is the weight matrix; The square root of the weight matrix; The total number of spherical harmonic coefficients; is the maximum order of the spherical harmonic function.

[0086] The physical meaning of the objective function: The VTEC value is predicted based on the spherical harmonic coefficient; The actual observed VTEC value; The residual between the predicted and observed values; For weighted residuals; This is the sum of squares of the weighted residuals. The objective function is minimized. The optimal spherical harmonic coefficients can be obtained. This minimizes the weighted sum of squared residuals between the predicted and observed values.

[0087] Example 7: Detailed Implementation of Solving Normal Equations This embodiment details the method for solving the normal equation in step S10.

[0088] 7.1 Derivation of the normal equation: For the objective function Differentiate and set it equal to zero: The normal equation is obtained by rearranging the equations: ;in, For matrix The transpose of .

[0089] 7.2 Solving the Normal Equations: The normal equations are a system of linear equations that can be solved using various numerical methods. (1) SVD decomposition method: for the matrix Perform Singular Value Decomposition (SVD): ;in, and It is an orthogonal matrix. Let be a diagonal matrix. Then the solution for the spherical harmonic coefficients is: The advantage of the SVD decomposition method is its good numerical stability and applicability to ill-conditioned matrices.

[0090] (2) Cholesky decomposition method: due to It is a symmetric positive definite matrix, and can be decomposed using Cholesky's method. ;in, Let be a lower triangular matrix. Then the solution for the spherical harmonic coefficients is: The advantages of the Cholesky decomposition method are its high computational efficiency and suitability for matrices with good condition numbers.

[0091] 7.3 Output of spherical harmonic coefficients: Output the estimated spherical harmonic coefficients. This is used for subsequent VTEC prediction and ionospheric modeling. Spherical harmonic coefficients It contains spatial distribution information of ionospheric VTEC, which can be used to predict VTEC values ​​at any location; analyze the spatiotemporal variation characteristics of the ionosphere; and provide ionospheric correction information for GNSS positioning.

[0092] Example 8: Detailed Implementation of Enhanced Design Matrix Construction This embodiment details the method for constructing the enhanced design matrix in step S4.

[0093] 8.1 Pre-assignment of the design matrix: Pre-assignment of the design matrix ,in: This represents the total number of observations. The number of spherical harmonic coefficients; This is the maximum order of the spherical harmonic function.

[0094] 8.2 GNSS Observation Filling Design Matrix: For the first GNSS observation points Calculate its spherical harmonic function basis function row vectors And fill it into the design matrix. The OK.

[0095] 8.3 Occultation Observation Filling Design Matrix: For the first... An occultation observation, using the location of its electron density peak. ,in, , These represent the longitude and latitude at the peak of the ionospheric electron density, respectively. The design matrix rows are calculated and filled according to the method in step 8.2.

[0096] 8.4 GIM Constraint Point Filling Design Matrix: For the first GIM constraint points Calculate the row vectors of the spherical harmonic function basis functions according to the method in step 8.2. And fill it into the design matrix. The Okay. It's important to note that adding GIM constraint points as virtual observations to the design matrix increases the amount of observation information in sparse regions and improves the condition number of the design matrix.

[0097] 8.5 Construction of Observation Vectors and Weight Matrix: Constructing Observation Vectors : For GNSS observation points, the observed value is the VTEC value; for occultation observation points, the observed value is the VTEC value corresponding to the peak electron density location; for GIM-constrained points, the observed value is... .

[0098] Constructing the weight matrix (Diagonal matrix): For GNSS and occultation observation points, the weights are determined based on the observation accuracy; for GIM constraint points, the weights are... .

[0099] Specific implementation of the present invention in practical applications: Application scenarios: This invention is applicable to the field of GNSS ionospheric modeling and regional enhancement technology, especially for improving the VTEC estimation accuracy in sparse GNSS stations such as oceans, polar regions, and remote areas.

[0100] Implementation steps: Step 1: Obtain the GIM data file in IONEX format, parse the file header information, and obtain the mesh parameters.

[0101] Step 2: Determine the geographical boundaries of the target enhancement area. For example, select the Pacific region, with latitudes from -60° to 60° and longitudes from 120° to 240°.

[0102] Step 3: Randomly generate 5183 constraint points within the target area and calculate the position of each constraint point.

[0103] Step 4: Extract the VTEC value of each constraint point from the IONEX data, filter out invalid values, and obtain the set of valid constraint points.

[0104] Step 5: Obtain the global information set of all GNSS stations and calculate the geographical distribution weight of each station.

[0105] Step 6: For each constraint point, find the nearest station, inherit its geographical distribution weight, and calculate the final weight.

[0106] Step 7: Construct the design matrix, including the row vectors of spherical harmonic basis functions for GNSS observation points, occultation observation points, and GIM constraint points.

[0107] Step 8: Construct the observation vector and weight matrix.

[0108] Step 9: Construct the weighted least squares objective function.

[0109] Step 10: Solve the normal equations to obtain the spherical harmonic coefficients.

[0110] The method of this invention improves the accuracy of VTEC estimation in the Pacific region by about 30%, reduces the condition number of the design matrix by about 50%, and significantly enhances model stability.

[0111] The beneficial effects of this invention are: This invention improves regional accuracy, increasing VTEC estimation accuracy by 20-40% in sparse GNSS station areas, particularly in marine, polar, and remote regions. It also enhances model stability by reducing the condition number of the design matrix, improving numerical stability, and decreasing the variance of spherical harmonic coefficient estimation. Furthermore, it improves adaptability by automatically adjusting constraint point weights based on station distribution without manual intervention, adapting to different regional characteristics. It offers regional flexibility, supporting enhancement configurations for any geographical area, allowing selection of enhancement regions based on actual needs. Finally, it improves computational efficiency, enabling multiple reuse of pre-calculated weights, reducing redundant computation overhead, and supporting parallel processing.

[0112] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for generating and weighting regional augmented observations based on the GIM model, characterized in that: Includes the following steps: S1. Based on the IONEX format GIM data file, parse the file header information to obtain the mesh parameters; S2. Determine the geographic boundary of the target augmentation area based on grid parameters; S3. Within the target enhancement area, randomly generate a set number of constraint points and obtain the positions of the constraint points; S4. For each constraint point, use interpolation to extract the global ionospheric VTEC value from the IONEX data and store the constraint point information as a constraint point set. S5. Based on the global information set of GNSS stations, calculate the geographical distribution weight of each GNSS station and generate a set of station weights. S6. Based on the set of constraint points and the set of station weights, for each constraint point, find the station closest to the constraint point, assign the geographical distribution weight of the station to the constraint point, and then combine it with the global constraint weight factor to calculate the constraint point weight. S7. Construct the design matrix: For each GNSS observation point, occultation observation point, and constraint point, calculate the row vector of the spherical harmonic function basis function and fill the row vector into the design matrix; S8. Based on the VTEC values ​​of each GNSS observation point, occultation observation point, and constraint point, construct an observation vector; based on the geographical distribution weights of each GNSS observation point, occultation observation point, and constraint point, construct a weight matrix. S9. Based on the design matrix, weight matrix, and observation vector, construct a weighted least squares objective function; S10. Based on the objective function, obtain and solve the normal equation, and then solve for the spherical harmonic coefficients; In step S5, the expression for calculating the geographic distribution weights is as follows: ； In the formula, For geographical distribution weights, , , These are the importance coefficients for latitude, distance, and density weighting factors, respectively. Latitude weighting factor Distance weighting factor Density weighting factor; ； In the formula, The latitude of the measuring station; ； In the formula, For the station The minimum great circle distance to other stations. This is the distance attenuation scale parameter; ； In the formula, For the purpose of measuring station Center, search radius The number of stations within the area This is the density threshold parameter; In step S6, the expression for the constraint point weights is as follows: ； ； ； In the formula, The index of the station closest to constraint point i. Let i be the distance from point i to station j. This represents the total number of stations. For the constraint point weights, The geographical weights inherited by constraint point i As a global constraint weighting factor, Distance constraint point Recent station index Surrounding search radius The number of stations within the area This is the density threshold parameter; In step S7, for each GNSS observation point, the expression for the row vector of the spherical harmonic function basis function is calculated as follows: ； In the formula, To design the k-th row vector of matrix A, Let n be a Legendre polynomial, where n is the order and m is the degree. , Let be the latitude and longitude of the k-th GNSS observation point.

2. The method for generating and weighting regional augmented observations based on the GIM model according to claim 1, characterized in that: In step S3, the expression for the position of the constraint point is as follows: ； ； In the formula, , Let the latitude and longitude of the i-th generated constraint point be... and The minimum and maximum latitudes of the target area; and Given the minimum and maximum longitudes of the target area, and In the interval Random numbers that are uniformly distributed within the range.

3. The method for generating and weighting regional enhanced observations based on the GIM model according to claim 1, characterized in that: In step S9, the expression for the weighted least squares objective function is as follows: ； ； In the formula, This represents the total number of observations. This refers to the total number of GNSS vertical electron observations. This refers to the total number of observed vertical electrons during occultation. The number of GIM constraint points. The objective function value, Let be the vector of spherical harmonic coefficients to be solved. To design the matrix, For the observation vector, This is the weight matrix. denoted as the total number of spherical harmonic coefficients, and N as the maximum order of the spherical harmonic function.

4. The method for generating and weighting regional enhanced observations based on the GIM model according to claim 1, characterized in that: In step S10, the expression for the normal equation is as follows: ； In the formula, To design the transpose of matrix A.

Images