LS-SVM / MPA-BP combined model fitting GNSS elevation anomaly method

Abstract

The invention discloses a method for fitting GNSS elevation anomaly by using an LS-SVM / MPA-BP combined model. The method comprises the following steps: S1, partitioning a measurement region; s2, selecting a sub-region as a current subregion; s3, an LS-SVM model is trained; s4, obtaining a fitting elevation abnormal value of a known point of the sub-region; s5, the MPA-BP model is trained; s6, obtaining fitting elevation abnormal values of to-be-measured points of the sub-regions; s7, obtaining a fitting elevation anomaly residual error of the to-be-measured point of the sub-region; s8, acquiring an elevation abnormal value of the current partition; and S9, selecting the sub-regions which are not selected as the current subregion, repeating the steps S3-S8 until all the sub-regions are selected, and obtaining the elevation anomaly of the measurement region. According to the method, the combined model is adopted to carry out partition fitting on the GNSS elevation anomaly, the method has good extrapolation capability and high generalization capability, the applicable geographical scene is improved, the nonlinear characteristic in GNSS elevation anomaly fitting can be better captured, and the problems of difficult initialization, low precision and the like are solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of engineering surveying, and particularly relates to a method for fitting GNSS height anomaly by using an LS-SVM / MPA-BP combined model. Background Art

[0002] In recent years, the global navigation satellite system (GNSS) measurement technology has become the main means of providing high-precision three-dimensional positioning data in the field of engineering surveying. Its height measurement mainly provides geodetic height, which is based on the reference ellipsoid. Among them, the English name of the global navigation satellite system is Global Navigation Satellite System, abbreviated as GNSS. However, the geodetic height cannot directly reflect the actual terrain undulation and gravity change, and normal height is usually required in engineering applications in China, which is based on the quasi-geoid. Therefore, in order to convert the geodetic height measured by GNSS into the normal height required in practice, it is necessary to calculate the difference between the geodetic height and the normal height, that is, the height anomaly.

[0003] At present, the GNSS height anomaly fitting methods include function model fitting method, earth gravity field model method, machine learning method, etc. The function model fitting methods all fit a trend surface similar to the height anomaly to replace the quasi-geoid of the fitting area. However, the quasi-geoid is a physical surface and it is difficult to be completely approximated by a mathematical function, so the fitting error generated is relatively large. The earth gravity field model method has a best accuracy of about 20 cm because it is necessary to obtain gravity data and the spatial resolution of the model is limited. As a new technology for processing non-linear data fitting, the machine learning model has developed rapidly in recent years, but its initial parameters are random, which greatly affects the accuracy of fitting the height anomaly.

[0004] Generally speaking, various fitting methods have their own advantages and disadvantages. However, the research of different scholars has shown that the elevation anomaly fitting accuracy of machine learning models is significantly better than the first two. The most commonly used machine learning models are the Support Vector Machine (SVM) model and the Back Propagation Neural Network (BPNN) model. Among them, the English name of the Support Vector Machine model is SupportVectorModel, abbreviated as SVM; the English name of the Back Propagation Neural Network model is BackPropagationNeuralNetwork, abbreviated as BPNN. Due to the flexibility of its kernel function, the SVM model can effectively capture the complex non-linear relationship between the elevation anomaly and its features such as Gaussian plane coordinates and Digital Elevation Model (DEM), improving the fitting accuracy. Among them, the English name of the Digital Elevation Model is Digital Elevation Model, abbreviated as DEM. And a regularization term is added to the objective function, which can balance the model complexity and error, reducing the risk of overfitting. However, it also has many disadvantages. For example: ① It is difficult to process large-scale data. For large-scale GNSS data, it may cause the kernel matrix to exceed the memory limit of the computer, resulting in inefficient processing; ② Its adaptability to non-uniform data distribution is limited. The sample-dense area will cause overfitting, and the sample-sparse area will cause prediction deviation; ③ The limitation of kernel function selection. When the regional elevation anomaly has multi-scale features, the local terrain undulation is uneven or the overall terrain is relatively flat, and a single kernel function may not be able to take into account both local and global features. For the BPNN, due to its simple network topology structure and strong self-adaptability, it is widely used in GNSS elevation anomaly fitting. However, since its initial parameters are randomly generated, it has a significant impact on the fitting accuracy. At the same time, the BPNN is optimized through the gradient descent algorithm and is prone to falling into local optimal solutions. In elevation anomaly fitting, if the data distribution is complex, the model may not be able to find the global optimal solution. In view of the defects of the BPNN, some scholars have proposed to use some optimization algorithms to determine the optimal initial parameters of the BPNN, thereby improving the fitting accuracy. Currently, the optimization algorithms that have been applied to elevation conversion in strip areas mainly include genetic algorithms and sparrow algorithms. These algorithms can effectively improve the elevation conversion accuracy compared with the traditional BPNN and have achieved good results. However, for the above two models, using a single machine learning model to fit GNSS elevation anomalies does not take into account the residual information of the fitted elevation anomalies. As a result, the accuracy of elevation anomaly fitting conversion is not high enough, and the capture of non-linear changes is not precise enough. Summary of the Invention

[0005] In order to solve the above technical problems, the present invention provides a method for fitting GNSS elevation anomalies using a LS-SVM / MPA-BP combined model, wherein a least squares algorithm is used to optimize the SVM model in part of the processing; and a marine predator algorithm is used to optimize the BP neural network model in part of the processing to achieve elevation conversion. The English name of the least squares algorithm is Least Squares, abbreviated as LS; the English name of the marine predator algorithm is Marine Predators Algorithm, abbreviated as MPA. At the same time, combined with the idea of ​​zoning fitting, the combined model is used to perform fine fitting on each area, better handle geographical complexity and uneven data distribution, and achieve high-precision elevation conversion of line engineering.

[0006] The present invention is achieved through the following technical solutions.

[0007] The present invention provides a LS-SVM / MPA-BP combined model fitting GNSS elevation anomaly method, comprising the following steps:

[0008] S1. Obtain a measurement area, partition the measurement area using an RF model, and obtain several sub-areas;

[0009] S2, selecting a sub-region as a current partition, wherein the current partition includes sub-region known points and sub-region test points;

[0010] S3, training an LS-SVM model according to the characteristic variables of the sub-region known points, wherein the LS-SVM model is used to calculate the fitted elevation anomaly value, wherein the characteristic variables of the sub-region known points include Gaussian projection coordinates of the sub-region known points, DEM data of the sub-region known points and measured elevation anomaly values;

[0011] S4, inputting the characteristic variables of the known points in the sub-region into the LS-SVM model to obtain the fitted elevation anomaly values ​​of the known points in the sub-region;

[0012] S5, training an MPA-BP model according to characteristic variables of sub-region known points and fitted elevation anomaly values ​​of sub-region known points, wherein the MPA-BP model is used to calculate fitted elevation anomaly residuals;

[0013] S6, inputting characteristic variables of the sub-region test points into the LS-SVM model to obtain the fitted elevation anomaly values ​​of the sub-region test points, wherein the characteristic variables of the sub-region test points include Gaussian projection coordinates of the sub-region test points and DEM data;

[0014] S7, inputting the characteristic variables of the sub-region test points and the fitted elevation anomaly values ​​of the sub-region test points into the MPA-BP model to obtain the fitted elevation anomaly residuals of the sub-region test points;

[0015] S8. Calculate the height anomaly value of the sub-region measurement points according to the height anomaly value fitted by the sub-region measurement points to be measured and the height anomaly residual of the sub-region measurement points to be measured, and obtain the height anomaly value of the current partition;

[0016] S9. Select the unselected sub-region as the current partition, and repeat steps S3 - S8 until all sub-regions have been selected, and obtain the height anomaly of the measurement area.

[0017] Preferably, the steps for obtaining the measurement area and partitioning the measurement area through the RF model to obtain several sub-regions are as follows:

[0018] S11. Obtain the measurement area and set the parameters of the RF model; among them, the measurement area includes several known points, and the parameters of the RF model include the number of decision trees, the number of training samples K, and the number of sampling sets M;

[0019] S12. Randomly and repeatedly draw K training samples from the known points in the measurement area with replacement through the RF model;

[0020] S13. Repeat step S12 M times to obtain M sampling sets;

[0021] S14. Train M classification results corresponding to the M sampling sets by constructing decision trees;

[0022] S15. Obtain the model output by combining the M classification results through a combination strategy to obtain several sub-regions.

[0023] Preferably, the steps for training the LS - SVM model according to the characteristic variables of the known points in the sub-region are as follows:

[0024] S31. Select points from the known points in the sub-region to obtain a training set and a validation set, and set the training error stop condition. The training set includes several training samples;

[0025] S32. Use the Gaussian projection coordinates of the known points in the sub-region and the DEM data of the known points in the sub-region as the characteristic variables of the training samples, and use the measured height anomaly value as the target variable;

[0026] S33. Determine the SVM objective function according to the characteristic variables and target variables of the training samples, determine the corresponding constraint conditions, and obtain the SVM model;

[0027] S34. Optimize the SVM model through the least squares algorithm and update the SVM model;

[0028] S35. Calculate the training error of the SVM model according to the validation set, and judge whether the training error of the SVM model meets the training error stop condition. If so, obtain the LS - SVM model; otherwise, return to step S34.

[0029] Preferably, the expression of the SVM objective function is as follows:

[0030]

[0031] where D is the number of training samples, m is the dimension of the training samples, W ∈ R m is the weight vector, C is the regularization parameter, and ε i is the tolerance deviation;

[0032] The expression of the corresponding constraint condition is:

[0033] y i = W T g(x i ) + b + ε i i = 1,..., D

[0034] where D is the number of training samples, g(x i ) is the mapping function, b is the bias, and ε i is the tolerance deviation.

[0035] Preferably, the SVM model is optimized by the least squares algorithm, and updating the SVM model includes the following steps:

[0036] S341. Obtain the optimization conditions and define the Lagrangian function according to the SVM model;

[0037] S342. Calculate the Lagrangian function according to the optimization conditions to obtain the SVM optimization equations;

[0038] S343. Rearrange the SVM optimization equations to obtain the SVM model;

[0039] The expression of the Lagrangian function is:

[0040]

[0041] where D is the number of training samples, W ∈ R m is the weight vector, C is the regularization parameter, ε i is the tolerance deviation, g(x i ) is the mapping function, b is the bias, and a is the Lagrange multiplier;

[0042] The expression of the SVM optimization equations is:

[0043]

[0044] where m is the dimension of the training samples, W ∈ R m is the weight vector, C is the regularization parameter, ε i is the tolerance deviation, g(x i) is the mapping function, b is the bias, and a i is the Lagrange multiplier;

[0045] The expression of the SVM model is:

[0046]

[0047] where D is the number of training samples, g(x i ) is the mapping function, b is the bias, and a i is the Lagrange multiplier.

[0048] Preferably, the training of the MPA-BP model according to the characteristic variables of the known points in the sub-region and the fitting height anomaly values of the known points in the sub-region includes the following steps:

[0049] S51. Initialize the BP neural network structure, set the error precision requirement, and the BP neural network structure includes an input layer, a hidden layer, an output layer, weights, and bias values;

[0050] S52. Calculate the fitting height anomaly residual according to the measured height anomaly value and the fitting height anomaly value of the known points in the sub-region;

[0051] S53. Use the Gaussian coordinates projection of the known points in the sub-region, the DEM data of the known points in the sub-region, and the fitting height anomaly values of the known points in the sub-region as the input layer, and the fitting height anomaly residual as the expected output layer to train the BP neural network structure;

[0052] S54. Use the MPA algorithm to correct the weights and bias values and update the BP neural network structure;

[0053] S55. Calculate the training error of the BP neural network structure, determine whether it meets the error precision requirement. If so, obtain the MPA-BP neural network model; otherwise, return to step S54.

[0054] Preferably, the use of the MPA algorithm to correct the weights and bias values and update the BP neural network structure includes the following steps:

[0055] S541. Set the fitness precision requirement, and initialize the prey matrix and predator matrix according to the weights and bias values of the BP neural network model;

[0056] S542. Optimize the prey matrix by simulating the best encounter rate strategy between predators and prey, and update and record the prey matrix;

[0057] S543. Calculate the fitness of the prey matrix, detect the top predator, and update the predator matrix;

[0058] S544. Determine whether the fitness of the predator matrix meets the fitness accuracy requirement. If so, proceed to step S545; otherwise, return to step S542.

[0059] S545. Obtain the optimal weight position and the optimal bias value position according to the predator matrix, and use the weight at the optimal weight position and the bias value at the optimal bias value position as the initialization parameters of the BP neural network structure.

[0060] Preferably, the method for optimizing the prey matrix by simulating the best encounter rate strategy between predators and prey, updating and recording the prey matrix includes the following steps:

[0061] S5421. Simulate the high-speed escape behavior of the prey, obtain the high-speed escape optimization strategy, and optimize the prey matrix according to the high-speed escape.

[0062] S5422. Simulate the random swimming behavior when the speeds of the prey and the predator are comparable, obtain the escape strategy when the speeds are comparable, and optimize the prey matrix according to the escape strategy when the speeds are comparable.

[0063] S5423. Simulate the low-speed escape behavior of the prey, obtain the low-speed escape strategy, and optimize the prey matrix according to the low-speed escape strategy.

[0064] S5424. Simulate the formation of eddies and the effect of the fish school aggregation device, obtain the eddy and dispersion effect strategy, optimize the prey matrix according to the eddy and dispersion effect strategy, and record the prey matrix.

[0065] Preferably, the expression of the high-speed escape optimization strategy is:

[0066]

[0067] where represents the moving step size of the i-th individual, is a random number containing a normal distribution based on Brownian motion, n is the population size, P is a constant, is a vector of uniform random numbers in [0,1], prey is the prey matrix, and Elite is the predator matrix;

[0068] The escape strategy when the speeds are comparable specifically means: dividing the population into two parts, one part is the exploration population and the other part is the exploitation population;

[0069] The expression of the exploration population is:

[0070]

[0071] where represents the moving step size of the i-th individual, is a random number vector based on the Lévy distribution, n is the number of populations, P is a constant, is a uniform random number vector in [0, 1], prey is the prey matrix, and Elite is the predator matrix;

[0072] The expression for the developed population is:

[0073]

[0074]

[0075] where, represents the movement step size of the i-th individual, is a random number containing a normal distribution representing Brownian motion, n is the number of populations, P is a constant, prey is the prey matrix, Elite is the predator matrix, Iter is the current iteration number, Max_Iter is the maximum iteration number, and CF is the scaling factor;

[0076] The expression for the low-speed escape strategy is:

[0077]

[0078] where, represents the movement step size of the i-th individual, is a random number containing a normal distribution representing Brownian motion, n is the number of populations, P is a constant, prey is the prey matrix, Elite is the predator matrix, and CF is the scaling factor.

[0079] Preferably, the expression for the vortex and dispersion effect strategy is:

[0080]

[0081] where, FADs is the probability of the influence of the fish school aggregation device effect on the optimization process, is a binary vector including 0 and 1, r is a random number in [0, 1], X max and X min represent vectors of the lower dimension limit, the subscripts r1 and r2 represent random indices of the prey matrix, prey is the prey matrix, CF is the scaling factor, is a uniform random number vector in [0, 1].

[0082] The beneficial effects of the present invention are as follows:

[0083] 1. A method for partitioning and fitting GNSS height anomalies in a strip area using an LS-SVM / MPA-BP neural network combined model, which makes full use of the advantages of different height anomaly fitting models. Without the need to determine the geometric surface, it can solve the problems of local optimal values generated by the training set and uncertain initial parameters of the model in the space of a small number of GNSS / leveling sample points, thereby achieving high-precision conversion from geodetic height to normal height;

[0084] 2. Compared with the function model fitting method and the earth gravity field method, it can better capture the non-linear characteristics of height anomalies. Compared with the overall area fitting, it takes into account the influence of the undulation of the regional terrain and the unevenness of some data points on the fitting of GNSS height anomalies in the strip area using the LS-SVM model. Selecting a suitable kernel function for each area for refined fitting can more effectively reduce the deviation of height anomalies and meet the engineering requirements;

[0085] 3. Compared with the standard BP neural network, the MPA optimization algorithm can optimize the weights and thresholds of the BP neural network, making the weights and thresholds of the BP neural network better values during fitting, so as to improve the fitting accuracy of the fitting residuals of height anomalies. And compared with the existing optimized BP neural network algorithms, the MPA algorithm shows a better fitting effect in height conversion;

[0086] 4. Then, combining the GNSS height anomalies fitted by the LS-SVM model with the height anomaly residuals fitted by the MPA-BP neural network achieves the correction of the fitted values of GNSS height anomalies, obtaining a better height conversion result than fitting height anomalies with a single LS-SVM model;

[0087] 5. This method of partitioning and combined model is also applicable to the case where the height anomaly changes of long-distance and large-range strip terrains are complex, increasing the applicable range of the model for fitting GNSS height anomalies. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 is the flowchart of the method provided by the embodiment of the present invention;

[0089] Figure 2 is the flowchart of the LS-SVM / MPA-BP combined model for fitting sub-region height anomalies provided by the embodiment of the present invention;

[0090] Figure 3 is the flowchart of the LS-SVM model for fitting sub-region height anomalies provided by the embodiment of the present invention;

[0091] Figure 4 is the flowchart of the MPA optimized BP neural network structure provided by the embodiment of the present invention;

[0092] Figure 5It is a schematic diagram of the BP neural network fitting geoid height residual structure provided by an embodiment of the present invention. Detailed implementation manners

[0093] The technical solutions of the present invention will be further described below, but the scope of protection is not limited thereto.

[0094] As Figure 1 shown, a method for fitting GNSS geoid height using an LS-SVM / MPA-BP combined model includes the following steps:

[0095] S1. Obtain a measurement area, partition the measurement area through an RF model, and obtain a number of sub-areas;

[0096] The RF model, full name Random Forest model, that is, a random forest model. A random forest refers to a classifier that uses multiple decision trees to train and predict samples. It contains a classifier with multiple decision trees, and the output category is determined by the mode of the categories output by individual trees. Random forest is a flexible and easy-to-use machine learning algorithm that can generally obtain good results even without hyperparameter tuning.

[0097] This embodiment is applied to long-distance line projects, and the measurement area is a long-distance strip-shaped area. Among them, specific long-distance line projects include railways, roads, pipelines, etc.

[0098] This step includes the following content:

[0099] S11. Obtain a measurement area and set the parameters of the RF model; wherein, the measurement area includes a number of known points, and the parameters of the RF model include the number of decision trees, the number of training samples K, and the number of sampling sets M;

[0100] In this embodiment, setting the parameters of the RF model includes:

[0101] The number of clustering categories: a machine learning parameter related to the number of input data, which is set by oneself. It should be noted that setting it too small will result in insufficient training samples, and setting it too large cannot capture the influence of terrain features well;

[0102] The number of trees in the random forest: the default value is 100. Generally speaking, if the number of trees is too small, it is easy to underfit. If the number of trees is too large, it will increase its computational amount, and when the number of trees reaches a certain amount, increasing the number of trees further will result in little improvement in the obtained model. Therefore, a moderate value is usually selected, and it is recommended to be between 100 and 150;

[0103] The maximum number of features: When selecting the most suitable attribute, the features for partitioning cannot exceed this value, which is related to the number of feature quantities M of the input data. Generally, the default "auto" is used, which means that each node randomly considers when partitioning For a feature, if there are a relatively large number of features, the following ratio can be used to divide the number of feature values to control the generation time of the decision tree: log 2 M, M,

[0104] Maximum depth: By default, it can be left blank. If not entered, the decision tree will not limit the depth of the subtree when building the subtree, which will result in each leaf node having only one class. Generally, this value can be ignored when there is little data or few features. If the model has a large number of samples and many features, it is recommended to limit this maximum depth. The specific value depends on the data distribution and can usually be between 10 and 100;

[0105] Minimum number of samples required for further partitioning of internal nodes: When partitioning nodes based on attributes, the minimum number of samples for each partition. If the number of samples at a certain node is less than this value, no further partitioning will be performed, and the default value is 2. When the sample order of magnitude is very large, it is recommended to increase this value.

[0106] Minimum number of samples in leaf nodes: This value limits the minimum number of samples in leaf nodes. If the number of samples in a certain leaf node is less than the number of samples, it will be pruned together with other sibling nodes, and only the parent node will be retained. The default value is 1. If the sample order of magnitude is very large, it is recommended to increase this value.

[0107] S12. Randomly and repeatedly draw K training samples with replacement from the known points in the measurement area through the RF model;

[0108] S13. Repeat step S12 M times to obtain M sampling sets;

[0109] S14. Train M classification results corresponding to the M sampling sets by constructing decision trees;

[0110] S15. Obtain the model output through the combination strategy for the M classification results to obtain several sub-regions.

[0111] S2. Select a sub-region as the current partition, and the current partition includes sub-region known points and sub-region points to be measured;

[0112] S3. Train the LS-SVM model based on the feature variables of the sub-region known points. The LS-SVM model is used to calculate the fitted elevation anomaly value. The feature variables of the sub-region known points include the Gaussian projection coordinates of the sub-region known points, the DEM data of the sub-region known points, and the measured elevation anomaly value;

[0113] The LS-SVM model is the SVM model optimized by the LS algorithm. In this embodiment, the measured elevation anomaly value is denoted as ξ.

[0114] On the basis of completing the zoning, to achieve high-precision elevation anomaly fitting with the combined model, it is necessary to construct the fitting model for each area separately. The specific training steps are illustrated with a sub-area. For example, Figure 2 As shown, the steps of training the LS-SVM model according to the characteristic variables of the known points in the sub-area include the following:

[0115] S31. Select points from the known points in the sub-area, obtain the training set and the validation set, and set the training error stop condition. The training set includes several training samples;

[0116] For example, Figure 3 As shown, take the point selection result as the training set D=(x i ,y i ), where x i is the characteristic variable of the training sample; i is the number of training samples; y i is the target variable of the training sample. The remaining known points are used as the validation set to evaluate the model performance.

[0117] S32. Use the Gaussian projection coordinates of the known points in the sub-area and the DEM data of the known points in the sub-area as the characteristic variables of the training samples, and use the measured elevation anomaly value as the target variable;

[0118] S33. Determine the SVM objective function according to the characteristic variables and target variables of the training samples, determine the corresponding constraint conditions, and obtain the SVM model;

[0119] According to the principle of minimizing the structural risk, convert the function fitting problem into a function constraint optimization problem, that is, the expression of the SVM objective function is obtained as:

[0120]

[0121] where D is the number of training samples, m is the dimension of the training samples, W∈R m is the weight vector, C is the regularization parameter, ε i is the tolerance deviation;

[0122] The expression of the corresponding constraint condition is:

[0123] y i =W T g(x i )+b+ε i i = 1,..., D

[0124] where D is the number of training samples, g(x i ) is the mapping function, b is the bias, and ε i is the tolerance deviation.

[0125] S34. Optimize the SVM model through the least squares algorithm and update the SVM model;

[0126] The step of optimizing the SVM model through the least squares algorithm and updating the SVM model includes the following steps:

[0127] S341. Obtain the optimization conditions and define the Lagrangian function according to the SVM model;

[0128] The expression of the Lagrangian function is:

[0129]

[0130] where D is the number of training samples, W ∈ R m is the weight vector, C is the regularization parameter, ε i is the tolerance deviation, g(x i ) is the mapping function, b is the bias, a i is the Lagrange multiplier;

[0131] S342. Calculate the Lagrangian function according to the optimization conditions to obtain the SVM optimization equations;

[0132] Specifically, according to the optimization conditions, respectively find and set the result to zero, then the expression of the SVM optimization equations obtained is:

[0133]

[0134] where m is the dimension of the training samples, W ∈ R m is the weight vector, C is the regularization parameter, ε i is the tolerance deviation, g(x i ) is the mapping function, b is the bias, a i is the Lagrange multiplier;

[0135] S343. Rearrange the SVM optimization equations to obtain the SVM model;

[0136] The expression of the SVM model is:

[0137]

[0138] where D is the number of training samples, g(x i ) is the mapping function, b is the bias, a i is the Lagrange multiplier.

[0139] S35. Calculate the training error of the SVM model according to the validation set, and determine whether the training error of the SVM model satisfies the training error stop condition. If so, obtain the LS-SVM model; otherwise, return to step S34.

[0140] If the training error of the SVM model satisfies the training error stop condition, then use this SVM model as the LS-SVM model.

[0141] After obtaining the LS-SVM model, introduce the kernel function K(x, x i ) to replace g(x)g(x i ), to solve the high-dimensional problem, that is:

[0142]

[0143] where D is the number of training samples, K(x, x i ) is the kernel function, b is the bias, and a i is the Lagrange multiplier.

[0144] S4. Input the characteristic variables of the known points in the sub-region into the LS-SVM model to obtain the fitting elevation anomaly values of the known points in the sub-region;

[0145] In this embodiment, the following parameters also need to be set when inputting the characteristic variables of the known points in the sub-region into the LS-SVM model:

[0146] Regularization hyperparameter: Controls the tolerance of the model to errors. If set to a larger value, the loss function is larger and overfitting is likely to occur, which is suitable for fitting elevation anomalies in flat areas. Because the terrain distribution is relatively uniform and the data is less affected by the terrain, high-precision elevation anomalies can be fitted; if set to a smaller value, the algorithm complexity is smaller, emphasizing the smoothness of the model, which is suitable for areas with uneven terrain distribution, such as mountains and river valleys; but it may lead to underfitting and unable to fully learn the non-linear characteristics of the data. The default value is 1. In the case where the data characteristics are not clear, it is recommended to search for the optimal value in a relatively wide range: 1 - e3 to 1 + e3.

[0147] Kernel function hyperparameter: For the fitting of non-linear elevation anomaly values, it is most applicable to use the radial basis kernel function, abbreviated as RBF, which is also called the Gaussian kernel function. Therefore, the setting of the kernel function hyperparameter here is mainly for the kernel width parameter of the RBF kernel function. A larger kernel function width ignores the detailed changes and emphasizes the global characteristics, which is suitable for areas with gentle elevation anomaly changes and relatively uniform terrain distribution. A smaller kernel function width emphasizes the local changes more and is suitable for processing areas with rapid elevation anomaly changes. It is recommended to determine the range of the hyperparameter between 0.1 and 100, or it can also be determined by the following formula:

[0148]

[0149] In the formula: σ is the kernel function width, n_distance represents the distance from all sample points in the two-dimensional space to the point (0, 0), and k is a constant between [0.5, 2].

[0150] In this embodiment, the fitted height anomaly value of the known points in the sub-region is denoted as ξ'.

[0151] S5. Train the MPA-BP model according to the characteristic variables of the known points in the sub-region and the fitted height anomaly value of the known points in the sub-region. The MPA-BP model is used to calculate the fitted height anomaly residual.

[0152] The MPA-BP model is a BP model optimized by the MPA algorithm.

[0153] As Figure 4 shown, training the MPA-BP model according to the characteristic variables of the known points in the sub-region and the fitted height anomaly value of the known points in the sub-region includes the following steps:

[0154] S51. Initialize the BP neural network structure and set the error precision requirement. The BP neural network structure includes an input layer, a hidden layer, an output layer, weights, and bias values.

[0155] The weights and bias values are necessary parameters connecting each layer, and the number of hidden layers can be set according to requirements.

[0156] The GNSS height anomaly fitting belongs to a simple regression problem. Therefore, in this embodiment, only a simple BP neural network with a single hidden layer structure is required, as Figure 5 shown. Assume that the input layer sample vector is the output layer sample vector is the hidden layer output vector is the expected output sample vector is the parameters from the input layer to the hidden layer are the weights and the bias value the parameters from the hidden layer to the output layer are the weights W (2) and the bias value b (2) , and f(x) is selected as the activation function for each layer.

[0157] Then, the output vector from the input layer to the hidden layer can be obtained by the following formula:

[0158]

[0159] o j = f(z j )

[0160] where z j represents the input value of the j-th neuron in the hidden layer, represents the weight from the i-th x input vector in the input layer to the j-th neuron in the hidden layer, represents the bias value of the j-th neuron in the hidden layer, and o j represents the output value of the j-th neuron in the hidden layer.

[0161] Similarly, the output vector from the hidden layer to the output layer is as follows:

[0162]

[0163] y n = f(z n )

[0164] where z n represents the input value of the nth neuron in the output layer, represents the weight value from the jth neuron in the hidden layer to the nth neuron in the output layer, represents the bias value of the nth neuron in the output layer, and y n represents the nth output layer sample vector.

[0165] S52. Calculate the fitting elevation anomaly residual according to the measured elevation anomaly value and the fitting elevation anomaly value of the known points in the sub-region;

[0166] The expression for the fitting elevation anomaly residual is:

[0167] Δξ = ξ - ξ′

[0168] where ξ' is the fitting elevation value of the known points in the sub-region, and ξ is the measured elevation anomaly value.

[0169] S53. Use the Gaussian coordinate projection of the known points in the sub-region, the DEM data of the known points in the sub-region, and the fitting elevation anomaly value of the known points in the sub-region as the input layer, and the fitting elevation anomaly residual as the expected output layer to train the BP neural network structure;

[0170] Using the Gaussian coordinate projection of the known points in the sub-region, the DEM data of the known points in the sub-region, and the fitting elevation anomaly value of the known points in the sub-region as the input layer, the error function can be obtained according to the output vector from the input layer to the hidden layer and the output vector from the hidden layer to the output layer:

[0171]

[0172] where y n represents the nth output layer sample vector, E represents the training error value, and d n represents the nth expected output layer sample vector;

[0173] Then, according to the error function, the gradient descent method is used to correct the weight value and the bias value. The adjustment magnitude of the weight value is proportional to the gradient descent of the error:

[0174]

[0175]

[0176] where, It represents the weight correction value from the $i$-th input vector $\mathbf{x}$ in the input layer to the $j$-th neuron in the hidden layer; It represents the weight correction value from the $j$-th neuron in the hidden layer to the $n$-th neuron in the output layer; the negative sign represents gradient descent, and the proportionality coefficient $\eta\in(0,1)$, which is often called the learning efficiency of the BP neural network.

[0177] In this embodiment, the following BP neural network related parameters also need to be set:

[0178] Number of nodes in the hidden layer: A machine learning parameter, which can be determined by the trial-and-error method or the formula method. It is recommended to be between 5 and 20. The following formula can be referred to:

[0179] $S = \log$ 2 $n$

[0180] $S = 2n + 1$

[0181] $S < n - 1$

[0182]

[0183] where $S$ is the number of neurons in the hidden layer, $n$ is the dimension of the input data, $m$ is the dimension of the output data, and $\alpha$ is a constant between $[0,10]$.

[0184] Number of nodes in the input layer: A machine learning parameter, which is related to the dimension of the input data. If only the planar coordinates are input, it is 2. If the input data includes planar coordinates and DEM data, it is 3;

[0185] Number of nodes in the output layer: A machine learning parameter, which is related to the dimension of the output data. Generally, the output value is the height anomaly value, so the number of nodes in the output layer is 1;

[0186] Momentum factor: A machine learning parameter, which is between 0 and 1, generally 0.2;

[0187] Target accuracy: A machine learning parameter, which is set by oneself. When the accuracy reaches the target accuracy, the BP neural network stops iterating and outputs the final result;

[0188] Number of iterations: The number of iterations of the BP neural network. When the predetermined number of iterations is reached, the calculation stops and the result is output.

[0189] After iteration, the BP neural network structure is obtained.

[0190] S54. Use the MPA algorithm to correct the weights and bias values, and update the BP neural network structure;

[0191] The MPA optimization algorithm can optimize the weights and thresholds of the BP neural network, so that the weights and thresholds are relatively optimal values during the fitting of the BP neural network, thereby improving the fitting accuracy of the fitting height anomaly residuals.

[0192] In this embodiment, the following relevant parameters of the MPA algorithm also need to be set:

[0193] Maximum number of generations: The number of times the MPA algorithm evolves, greater than 0. The larger the value, the more iteration times;

[0194] Population size: A machine learning parameter, greater than 0, default is 0. The larger the value, the larger the population size and the more individuals participating in the optimization. However, if the number is too large, the calculation efficiency will be reduced. Generally, it is set to about 20.

[0195] The method of using the MPA algorithm to correct the weight and bias values and update the BP neural network structure includes the following steps:

[0196] S541. Set the fitness accuracy requirement, and initialize the prey matrix and predator matrix according to the weight and bias values of the BP neural network model;

[0197] Similar to most metaheuristic methods, MPA is a population-based method. The initialization method for each element is a uniform distribution over the search space:

[0198] X 0 = X(Xmin max ) min

[0199] where X max and X min are the lower and upper bounds of the variable, and rand is a uniform random vector in the range of 0 to 1.

[0200] After each element is initialized, the first prey matrix is formed:

[0201]

[0202] where X n,d represents the d-th dimension of the n-th prey.

[0203] Calculate the fitness of each prey individual to obtain the optimal individual vector Copy this vector n times to construct the Elite matrix.

[0204]

[0205] where n is the size of the population and d is the number of dimensions of the problem.

[0206] At the end of each iteration, if there is a top predator with better fitness, then the current top predator will be replaced by it, and the Elite matrix will be updated.

[0207] S542. Optimize the prey matrix by simulating the optimal encounter rate strategy between predators and prey, and update and record the prey matrix;

[0208] The step of optimizing the prey matrix by simulating the optimal encounter rate strategy between predators and prey, and updating and recording the prey matrix includes the following steps:

[0209] S5421. Simulate the high-speed escape behavior of the prey, obtain the high-speed escape optimization strategy, and optimize the prey matrix according to the high-speed escape;

[0210] Simulating the high-speed escape behavior of the prey is the first optimization stage. When the high-speed ratio appears or the prey moves faster than the predator, this generally occurs in the first third of the iteration. At this time, obtain the high-speed escape optimization strategy and optimize the prey matrix according to the high-speed escape. The expression of the high-speed escape optimization strategy is:

[0211]

[0212]

[0213] Where, represents the moving step of the i-th individual, is a random number containing a normal distribution based on Brownian motion, n is the population size, P is a constant, is a vector of uniform random numbers in [0,1], prey is the prey matrix, and Elite is the predator matrix;

[0214] S5422. Simulate the random swimming behavior when the speeds of the prey and the predator are equal, obtain the speed-equivalent escape strategy, and optimize the prey matrix according to the speed-equivalent escape strategy;

[0215] Simulating the random swimming behavior when the speeds of the prey and the predator are equal is the second optimization stage. When the unit speed ratio or both the predator and the prey move at the same speed, both are looking for prey. This usually occurs when the iteration number is greater than one-third of the maximum iteration number and less than two-thirds of it. At this time, obtain the speed-equivalent escape strategy and optimize the prey matrix according to the speed-equivalent escape strategy. The specific meaning of the speed-equivalent escape strategy is: divide the population into two parts, one part is the exploration population and the other part is the exploitation population;

[0216] The expression of the exploration population is:

[0217]

[0218] Where, represents the moving step of the i-th individual, is a vector of random numbers based on the Lévy distribution, n is the population size, P is a constant, is a vector of uniformly random numbers in [0, 1], prey is the prey matrix, and Elite is the predator matrix;

[0219] The expression of the developed population is:

[0220]

[0221]

[0222] where, represents the moving step of the i-th individual, is a random number containing a normal distribution based on Brownian motion, n is the population number, P is a constant, prey is the prey matrix, Elite is the predator matrix, Iter is the current iteration number, Max_Iter is the maximum iteration number, and CF is the scaling factor;

[0223] S5423. Simulate the low-speed escape behavior of the prey, obtain the low-speed escape strategy, and optimize the prey matrix according to the low-speed escape strategy;

[0224] Simulating the low-speed escape behavior of the prey is the third optimization stage. When the low-speed ratio or the predator's moving speed is faster than the prey, this usually occurs in the last stage of optimization, that is, when the iteration number is greater than two-thirds of the maximum iteration number. At this time, obtain the low-speed escape strategy and optimize the prey matrix according to the low-speed escape strategy. The expression of the low-speed escape strategy is:

[0225]

[0226] where, represents the moving step of the i-th individual, is a random number containing a normal distribution based on Brownian motion, n is the population number, P is a constant, prey is the prey matrix, Elite is the predator matrix, and CF is the scaling factor.

[0227] S5424. Simulate the formation of eddies and the effect of the fish aggregation device, obtain the eddy and dispersion effect strategy, optimize the prey matrix according to the eddy and dispersion effect strategy, and record the prey matrix.

[0228] In the last stage of iteration, utilize the eddy and dispersion effect to further improve the search ability of the algorithm through random perturbation, and avoid the BP neural network hovering near the local optimal solution. The expression of the eddy and dispersion effect strategy is:

[0229]

[0230] where, FADs is the probability of the influence of the fish aggregation device effect on the optimization process, is a binary vector including 0 and 1, r is a random number in [0, 1], X max and X min is a vector representing the lower bound of the dimension, the subscripts r1 and r2 represent random indices of the prey matrix, prey is the prey matrix, CF is the scaling factor, is a vector of uniformly random numbers in [0, 1].

[0231] Specifically, FADs is the probability of the influence of the fish aggregating device effect on the optimization process, generally 0.2; is a binary vector including 0 and 1, which is constructed by generating a random vector in [0, 1], if the array is less than 0.2, change its array to 0, if the array is greater than 0.2, change its array to 1. r is a random number in [0, 1]; X max and X min is a vector representing the lower bound of the dimension; the subscripts r1 and r2 represent random indices of the prey matrix. When r ≤ FADs, the predator will make longer jumps in different dimensions to find other optimal predation ranges, so as to achieve the effect of jumping out of the local optimum. When r > FADs, the predator will move randomly within the current predation range.

[0232] Among them, the rule for recording the prey matrix is the ocean memory, specifically: calculate the fitness value of each individual in the prey matrix, and use the ocean memory to update the historical positions of each prey in the prey matrix, and save the positions with higher fitness to ensure that the individuals will not degenerate. If the fitness of an individual prey is better than the fitness of the corresponding position in the Elite matrix, then replace the corresponding individual in the original Elite matrix with this individual, and there will be random mutations and re-explorations during the iteration process to prevent early convergence of the population and help find better positions.

[0233] S543. Calculate the fitness of the prey matrix, detect the top predator, and update the predator matrix;

[0234] S544. Judge whether the fitness of the predator matrix meets the fitness accuracy requirement. If so, go to step S545; otherwise, return to step S542;

[0235] S545. Obtain the optimal weight position and the optimal bias value position according to the predator matrix, and use the weight at the optimal weight position and the bias value at the optimal bias value position as the initialization parameters of the BP neural network structure.

[0236] S55. Calculate the training error of the BP neural network structure, and judge whether it meets the error accuracy requirement. If so, obtain the MPA - BP neural network model; otherwise, return to step S54.

[0237] S6. Input the characteristic variables of the points to be measured in the sub-region into the LS-SVM model to obtain the fitted elevation anomaly values of the points to be measured in the sub-region. The characteristic variables of the points to be measured in the sub-region include the Gauss projection coordinates of the points to be measured in the sub-region and the DEM data;

[0238] In this embodiment, the fitted elevation anomaly value of the points to be measured in the sub-region is denoted as ξ”.

[0239] S7. Input the characteristic variables of the points to be measured in the sub-region and the fitted elevation anomaly values of the points to be measured in the sub-region into the MPA-BP model to obtain the fitted elevation anomaly residuals of the points to be measured in the sub-region;

[0240] In this embodiment, the fitted elevation anomaly residual of the points to be measured in the sub-region is denoted as Δξ'.

[0241] S8. Calculate the elevation anomaly values of the points to be measured in the sub-region based on the fitted elevation anomaly values of the points to be measured in the sub-region and the fitted elevation anomaly residuals of the points to be measured in the sub-region, and obtain the elevation anomaly values of the current partition;

[0242]

[0243] Among them, H GNSS is the geodetic height of the point to be measured, H 0 is the normal height of this point, and the difference between the two is the elevation anomaly.

[0244] S9. Select the unselected sub-regions as the current partition, and repeat steps S3 to S8 until all sub-regions have been selected, to obtain the elevation anomaly of the measurement area.

[0245] So far, each sub-region uses the LS-SVM model to fit the GNSS elevation anomaly, and uses the MPA-BP neural network model to fit its residual as the correction term of the elevation anomaly to improve the fitting accuracy of the elevation anomaly.

[0246] The present invention adopts a method of partitioning and fitting the GNSS height anomaly in a strip area using an LS-SVM / MPA-BP neural network combined model. By making full use of the advantages of different height anomaly fitting models, without the need to determine the geometric surface, it can solve the problems of local optimal values generated by the training set and uncertain initial parameters of the model in the space of a small number of GNSS / level sample points, thus achieving a high-precision conversion from geodetic height to normal height. Compared with the function model fitting method and the earth gravity field method, it can better capture the non-linear characteristics of the height anomaly. Compared with the overall area fitting, it takes into account the influence of the terrain undulation and uneven distribution of data points in the strip area on the fitting of the GNSS height anomaly by using the LS-SVM model. Selecting a suitable kernel function for each area for refined fitting can more effectively reduce the deviation of the height anomaly and meet the engineering requirements. Compared with the standard BP neural network, the MPA optimization algorithm can optimize the weights and thresholds of the BP neural network, making the weights and thresholds of the BP neural network optimal values during fitting, thereby improving the fitting accuracy of the fitting residuals of the height anomaly. And compared with the existing optimized BP neural network algorithms, the MPA algorithm shows a better fitting effect in height conversion. Then, combining the GNSS height anomaly fitted by the LS-SVM model with the height anomaly residuals fitted by the MPA-BP neural network achieves the correction of the fitted value of the GNSS height anomaly, obtaining a better height conversion result than fitting the height anomaly with a single LS-SVM model. This method of partitioning and combined model is also applicable to the case where the height anomaly of long-distance and large-scale strip terrain changes complexly, increasing the applicable range of the model for fitting the GNSS height anomaly.

Claims

1. A method for fitting GNSS elevation anomaly using a LS-SVM / MPA-BP combined model, characterized in that: The following steps are involved: S1. Obtain a measurement area, partition the measurement area using an RF model, and obtain several sub-areas; S2, selecting a sub-region as a current partition, wherein the current partition includes sub-region known points and sub-region test points; S3, training an LS-SVM model according to the characteristic variables of the sub-region known points, wherein the LS-SVM model is used to calculate the fitted elevation anomaly value, wherein the characteristic variables of the sub-region known points include Gaussian projection coordinates of the sub-region known points, DEM data of the sub-region known points and measured elevation anomaly values; S4, inputting the characteristic variables of the known points in the sub-region into the LS-SVM model to obtain the fitted elevation anomaly values ​​of the known points in the sub-region; S5, training an MPA-BP model according to characteristic variables of sub-region known points and fitted elevation anomaly values ​​of sub-region known points, wherein the MPA-BP model is used to calculate fitted elevation anomaly residuals; S6, inputting characteristic variables of the sub-region test points into the LS-SVM model to obtain the fitted elevation anomaly values ​​of the sub-region test points, wherein the characteristic variables of the sub-region test points include Gaussian projection coordinates of the sub-region test points and DEM data; S7, inputting the characteristic variables of the sub-region test points and the fitted elevation anomaly values ​​of the sub-region test points into the MPA-BP model to obtain the fitted elevation anomaly residuals of the sub-region test points; S8, calculating the elevation anomaly value of the sub-region to be tested point according to the fitted elevation anomaly value of the sub-region to be tested point and the fitted elevation anomaly residual of the sub-region to be tested point, and obtaining the elevation anomaly value of the current partition; S9, select the sub-area that has not been selected as the current partition, repeat steps S3 to S8 until all sub-areas have been selected, and obtain the elevation anomaly of the measurement area.

2. The method for fitting GNSS elevation anomaly according to claim 1, characterized in that: The obtaining of the measurement area, partitioning the measurement area by the RF model, and obtaining a plurality of sub-areas comprises the following steps: S11, obtaining a measurement area and setting parameters of an RF model; wherein the measurement area includes a number of known points, and the parameters of the RF model include the number of decision trees, the number of training samples K, and the number of sampling sets M; S12, repeatedly randomly extracting K training samples with replacement from known points in the measurement area through the RF model; S13, repeat step S12 M times to obtain M sampling sets; S14, forming corresponding M classification results by constructing a decision tree training for the M sampling sets; S15. The M classification results are combined with a strategy to obtain a model output and a number of sub-regions.

3. The method for fitting GNSS elevation anomaly according to claim 1, characterized in that: The training of the LS-SVM model according to the characteristic variables of the known points in the sub-region comprises the following steps: S31, selecting known points in the sub-region, obtaining a training set and a validation set, and setting a training error stop condition, wherein the training set includes a plurality of training samples; S32, using the Gaussian projection coordinates of the known points in the sub-region and the DEM data of the known points in the sub-region as the characteristic variables of the training samples, and using the measured elevation anomaly values ​​as the target variables; S33, determining the SVM objective function according to the characteristic variables and target variables of the training samples, determining the corresponding constraint conditions, and obtaining the SVM model; S34, optimizing the SVM model by using a least squares algorithm, and updating the SVM model; S35. Calculate the training error of the SVM model based on the validation set, and determine whether the training error of the SVM model meets the training error stopping condition. If yes, obtain the LS-SVM model, otherwise return to step S34.

4. The method for fitting GNSS elevation anomaly according to claim 3, characterized in that: The expression of the SVM objective function is: Where D is the number of training samples, m is the dimension of training samples, W∈R m is the weight vector, C is the regularization parameter, ε i Tolerance for deviations; The expression of the corresponding constraint condition is: y i =W T g(x i )+b+ε i i=1,...,D Where D is the number of training samples, g(x i ) is the mapping function, b is the bias, ε i Tolerance deviation.

5. The method for fitting GNSS elevation anomaly according to claim 3, characterized in that: The method of optimizing the SVM model by the least squares algorithm and updating the SVM model comprises the following steps: S341, obtaining optimization conditions, and defining a Lagrangian function according to the SVM model; S342, calculating the Lagrangian function according to the optimization conditions to obtain the SVM optimization equation group; S343, sort out the SVM optimization equation group and obtain the SVM model; The expression of the Lagrangian function is: Where D is the number of training samples, W∈R m is the weight vector, C is the regularization parameter, ε i is the tolerance deviation, g(x i ) is the mapping function, b is the offset, a i is the Lagrange multiplier; The expression of the SVM optimization equation group is: Where m is the dimension of the training sample, W∈R m is the weight vector, C is the regularization parameter, ε i is the tolerance deviation, g(x i ) is the mapping function, b is the offset, a i is the Lagrange multiplier; The expression of the SVM model is: Where D is the number of training samples, g(x i ) is the mapping function, b is the offset, a i is the Lagrange multiplier.

6. The method for fitting GNSS elevation anomaly according to claim 1, characterized in that: The training of the MPA-BP model according to the characteristic variables of the sub-region known points and the sub-region known points fitting elevation anomaly values ​​comprises the following steps: S51, initializing the BP neural network structure and setting the error accuracy requirement, wherein the BP neural network structure includes an input layer, a hidden layer, an output layer, a weight and a bias value; S52, calculating the fitted elevation anomaly residual according to the measured elevation anomaly value and the fitted elevation anomaly value of the sub-region known point; S53, using the Gaussian coordinate projection of the sub-region known points, the DEM data of the sub-region known points and the fitted elevation anomaly values ​​of the sub-region known points as the input layer, and the fitted elevation anomaly residual as the expected output layer, to train the BP neural network structure; S54, using the MPA algorithm to correct weights and bias values, and update the BP neural network structure; S55, calculate the training error of the BP neural network structure and determine whether it meets the error accuracy requirement. If yes, obtain the MPA-BP neural network model, otherwise return to step S54.

7. The method for fitting GNSS elevation anomaly according to claim 6, characterized in that: The method of using the MPA algorithm to correct weights and bias values ​​and update the BP neural network structure includes the following steps: S541, setting the fitness accuracy requirement, and initializing the prey matrix and the predator matrix according to the weights and bias values ​​of the BP neural network model; S542, optimizing the prey matrix by simulating the optimal encounter rate strategy between the predator and the prey, and updating and recording the prey matrix; S543, calculating the fitness of the prey matrix, detecting the top predator, and updating the predator matrix; S544, determine whether the fitness of the predator matrix meets the fitness accuracy requirement, if yes, proceed to step S545, otherwise return to step S542; S545, obtaining the optimal weight position and the optimal bias position according to the predator matrix, and using the weight of the optimal weight position and the bias value of the optimal bias position as initialization parameters of the BP neural network structure.

8. The method for fitting GNSS elevation anomaly according to claim 7, characterized in that: The method of optimizing the prey matrix by simulating the optimal encounter rate strategy between the predator and the prey, and updating and recording the prey matrix comprises the following steps: S5421, simulating the high-speed escape behavior of prey, obtaining the high-speed escape optimization strategy, and optimizing the prey matrix according to the high-speed escape; S5422, simulate the random swimming behavior when the speeds of prey and predator are similar, obtain the speed-equivalent escape strategy, and optimize the prey matrix according to the speed-equivalent escape strategy; S5423, simulating the low-speed escape behavior of the prey, obtaining the low-speed escape strategy, and optimizing the prey matrix according to the low-speed escape strategy; S5424. Simulate the vortex formation and fish aggregation device effects, obtain the vortex and dispersion effect strategies, optimize the prey matrix according to the vortex and dispersion effect strategies, and record the prey matrix.

9. The method for fitting GNSS elevation anomaly according to claim 8, characterized in that: The expression of the high-speed escape optimization strategy is: in, represents the moving step length of the i-th individual, is a random number based on the normal distribution representing Brownian motion, n is the population size, P is a constant, is a uniform random number vector in [0,1], prey is the prey matrix, and Elite is the predator matrix; The speed-equal escape strategy specifically refers to: dividing the population into two parts, one part is the exploration population, and the other part is the development population; The expression of the exploration population is: in, represents the moving step length of the i-th individual, is a random number vector based on Lévy distribution, n is the population size, P is a constant, is a uniform random number vector in [0,1], prey is the prey matrix, and Elite is the predator matrix; The expression of the development population is: in, represents the moving step length of the i-th individual, is a random number based on the normal distribution representing Brownian motion, n is the population size, P is a constant, prey is the prey matrix, Elite is the predator matrix, Iter is the current iteration number, Max_Iter is the maximum iteration number, and CF is the scaling factor; The expression of the low-speed escape strategy is: in, represents the moving step length of the i-th individual, is a random number based on the normal distribution representing Brownian motion, n is the population size, P is a constant, prey is the prey matrix, Eltie is the predator matrix, and CF is the scaling factor.

10. The method for fitting GNSS elevation anomaly according to claim 8, characterized in that: The expression of the eddy current and dispersion effect strategy is: Among them, FADs is the probability of the fish aggregation device effect affecting the optimization process, is a binary vector containing 0 and 1, r is a random number in [0,1], X max and X min The vector representing the lower limit of the dimension, the subscripts r1 and r2 represent random indices of the prey matrix, prey is the prey matrix, CF is the scaling factor, is a vector of uniform random numbers in [0,1].

Images