Combined observation fusion method and system for slope radar and GNSS three-dimensional deformation data

Abstract

The invention discloses a slope radar and GNSS three-dimensional deformation data combined observation fusion method and system. The method comprises the following steps: S1, acquiring slope radar data, performing discrete value elimination processing, acquiring GNSS three-dimensional deformation data, and performing discrete value elimination processing and spatial interpolation processing; s2, calculating displacement data of the slope radar data and converting the displacement data to a BLH coordinate system; s3, constructing a combined observation model of the slope radar data and the GNSS three-dimensional deformation data under the same BLH coordinate system; and S4, carrying out data fusion on the slope radar data and the GNSS three-dimensional deformation data according to the weight by taking the iteratively optimized weight of the combined observation model. According to the method, the combined observation model of the multi-source data is jointly constructed based on the slope radar observation equation and the GNSS observation equation, the weighted residual algorithm based on the least square method is introduced, the weighted residual sum is minimized, the high-precision three-dimensional deformation fusion result is generated, and reliable technical support is provided for subsequent slope deformation analysis.

Description

Technical Field

[0001] The present invention relates to the field of slope deformation monitoring in mining areas, and particularly to a combined observation and fusion method and system for three-dimensional deformation data of slope radar and GNSS. Background Art

[0002] Both slope radar and GNSS technologies can be applied to slope deformation detection. The slope radar emits microwave signals, and the microwave signals will be reflected after encountering the slope surface, and the radar receives the reflected signals; according to the time difference, frequency change and other information between the transmitted signal and the received signal, the distance change between the slope surface and the radar can be calculated. The GNSS technology uses satellite signals emitted by GPS, Beidou satellite navigation systems, etc., and determines the three-dimensional position of the receiver (installed at the slope monitoring point) by receiving the signals emitted by the satellites, including longitude, latitude and elevation.

[0003] In the safety monitoring of open-pit mine slopes, the slope radar realizes the efficient monitoring of large-scale displacement and deformation of slopes through coherent synthetic aperture radar imaging technology; this technology can obtain high-resolution displacement information in a short time and is suitable for high-frequency, large-scale monitoring tasks, but the noise interference in the mine environment may reduce its accuracy. In addition, radar data usually can only provide one-dimensional displacement information along the line of sight direction and it is difficult to directly reflect the deformation characteristics of the three-dimensional space. The GNSS technology can provide high-precision and continuous displacement data and can also be applied to slope monitoring (in slope deformation monitoring, multiple GNSS receivers are installed at key positions on the slope; when the slope deforms, the positions of these receivers will change, and by comparing the position information at different time points, the deformation data of the slope can be obtained). The GNSS technology is suitable for long-term monitoring tasks and can accurately reflect the change of the three-dimensional coordinates of the monitoring point; however, GNSS monitoring points are usually relatively scattered, with a limited coverage range and it is difficult to provide comprehensive regional displacement information.

[0004] Both slope radar and GNSS technologies can be applied to slope deformation detection. How to combine slope radar and GNSS technologies for slope deformation detection is the current research direction; however, how to utilize the data obtained by slope radar and GNSS technologies respectively to realize the fusion of the two data (the technical advantages of the two can be complementary); moreover, constructing the complementary fusion of multi-source data of slope radar and GNSS technologies is of great significance for slope deformation detection and can provide more comprehensive technical support for slope safety monitoring in complex environments. Summary of the Invention

[0005] The purpose of the present invention is to overcome the data fusion problem of existing slope radar and GNSS deformation monitoring, and provide a joint observation fusion method of slope radar and GNSS three-dimensional deformation data, which utilizes the data complementarity of slope radar data and GNSS three-dimensional deformation data. The slope radar data and GNSS three-dimensional deformation data respectively constitute the slope radar observation equation and GNSS observation equation. Based on the two types of observation equations, a joint observation model of multi-source data is jointly constructed, and a weighted residual algorithm based on the least squares method is introduced. During the fusion process, the weights of the slope radar data and the GNSS three-dimensional deformation data are adjusted and the weighted residual sum is minimized, thereby finally generating a high-precision three-dimensional deformation fusion result. The present invention combines multi-source observation data, effectively utilizes the complementary advantages of different source data and fused observations, and can obtain high-precision three-dimensional displacement data of the target, providing reliable technical support for subsequent slope deformation analysis.

[0006] The purpose of the present invention is achieved through the following technical solutions:

[0007] A joint observation and fusion method for slope radar and GNSS three-dimensional deformation data, the method comprising:

[0008] S1. Obtain slope radar data and perform discrete value elimination processing, and obtain GNSS three-dimensional deformation data and perform discrete value elimination processing and spatial interpolation processing;

[0009] S2. Calculate the displacement data of the slope radar data and convert it into the BLH coordinate system;

[0010] S3. Construct a joint observation model of slope radar data and GNSS 3D deformation data in the BLH coordinate system. The weighted iterative optimization expression of the two types of observation data in the joint observation model is as follows: i represents the category of observation data, W i ′ represents the weight corresponding to the observed data category i after iteration, W i Represents the weight corresponding to the observed data category i before iteration, until the weight variances of the two types of observed data are approximately equal, that is,

[0011] S4. The joint observation model takes the weights after iterative optimization to fuse the slope radar data and the GNSS three-dimensional deformation data according to the weights.

[0012] In order to better implement the present invention, in method S3, the joint observation model solves the residual minimization process based on the least squares principle, and the result expression is as follows:

[0013] X=(A T PA) -1 A T PL, where A is the coefficient matrix, AT is the transpose of the coefficient matrix, P is the weight matrix, and L is the observation matrix of the joint observation model.

[0014] Preferably, the joint observation model consists of a slope radar observation equation and a GNSS observation equation.

[0015] The expression of the slope radar observation equation is as follows:

[0016] Among them, V los is the deformation value in the radar line-of-sight direction, S e , S n , S u are the projection vectors in the longitude, latitude, and elevation directions of the slope radar data respectively, D e , D n , D u are the deformation values in the longitude, latitude, and elevation directions to be calculated respectively, D′ los is the initial displacement in the radar line-of-sight direction, D′ e , D′ n , D′ u are the initial deformation values in the longitude, latitude, and elevation directions respectively;

[0017] The expression of the GNSS observation equation is as follows:

[0018] Among them, L e，gnss , L n，gnss , L u，gnss are the deformation values in the longitude, latitude, and elevation directions decomposed from the GNSS three-dimensional deformation data respectively, D e , D n , D u are the deformation values in the longitude, latitude, and elevation directions to be calculated respectively, V e，gnss , V n，gnss , V u，gnss are the residuals in the east-west, north-south, and vertical directions respectively;

[0019] The expression of the observation matrix L of the joint observation model combining the slope radar observation equation and the GNSS observation equation is as follows:

[0020]

[0021] is the fusion result X of the joint observation model.

[0022] Preferably, in method S1, the method for removing discrete values from GNSS three-dimensional deformation data is as follows:

[0023] Standardize the GNSS three-dimensional deformation data, calculate the mean of the GNSS three-dimensional deformation data within the sliding window, and preliminarily determine the data points deviating from the mean of the sliding window as discrete values;

[0024] Classify the preliminarily separated discrete values through a support vector machine with the goal of maximizing the classification margin. Use the label +1 to represent normal data points and -1 to represent discrete data points; the constraint condition expression of the support vector machine is: d n，svm (w svm T g n，svm +b svm )>1, g n，svm represents the combination of displacement, velocity, and acceleration vectors input for each data point, d n，svm represents the category of the samples after rough separation, w svm is the normal vector of the hyperplane, perpendicular to the decision hyperplane; b svm is the bias term and is used to adjust the position of the decision hyperplane; the constraint condition of the support vector machine classifies each data point on one side of the hyperplane and at least 1 away from the decision hyperplane;

[0025] Through the combination of the hyperplane normal vector w svm and the bias term b svm , calculate the distance of each data point to the hyperplane. If the distance is less than or equal to -1, determine the data point as a discrete value and perform elimination processing.

[0026] Preferably, for the discrete values that cannot be completely separated in the preliminary separation, introduce a soft margin by introducing a relaxation variable ξ n，svm , with as the goal, and the constraint condition expression is as follows:

[0027] d n，svm (w svm T g i，svm +b svm )>1 - ξ n，svm , where C svm is a hyperparameter used to control the maximization of the soft margin and the error penalty, and further separate and eliminate the discrete values that cannot be completely separated.

[0028] Preferably, the method for eliminating discrete values from slope radar data is as follows:

[0029] Preliminarily detect discrete values using proximity value analysis, then calculate the neighborhood set of each data point and the mean and standard deviation of the line-of-sight direction displacement within the neighborhood using the Euclidean distance, and judge suspected discrete values according to the following formula:

[0030] |d LOS，j -u N(j)| > k ld ·σ N(j) , where u N(j) , σ N(j) are the mean and variance in the field respectively, and d LOS，j is the displacement of data point j; k ld is a parameter preset based on the characteristics of slope radar historical data;

[0031] When constructing a random forest decision tree model, the splitting node selects the optimal splitting feature and threshold through the Gini index to maximize the reduction of the Gini index after splitting. t L , t R are the normal value node and discrete value node after splitting the slope radar data respectively, N L , N R are the sample numbers of the normal value and discrete value nodes respectively, N is the total number of samples of the slope radar data, G(t L ) is the Gini index of the normal value, and G(t R ) is the Gini index of the discrete value;

[0032] The discrete values in the slope radar data are determined by the majority voting of all decision trees in the random forest decision tree model and removed.

[0033] Preferably, in method S1, the spatial interpolation processing of GNSS three-dimensional deformation data is performed by the Kriging method using spatial correlation for interpolation, and the Kriging interpolation formula is as follows:

[0034] where [[ID=,41]]is the predicted value of the unknown position x0, [[ID=4,3]]is the value on the known data point x i with number i, and λ i is the Kriging weight with number i.

[0035] Preferably, in method S2, the expression of the three-dimensional displacement conversion relationship is as follows:

[0036] where S e , S n , S u are the projection vectors in the longitude, latitude, and elevation directions of the slope radar data respectively, (B i , L i , H i ) are the coordinates in the geocentric rectangular coordinate system respectively, and (B0, L0, H0) is the BLH coordinate of the radar equipment center point of the slope radar data.

[0037] A joint observation and fusion system for slope radar and GNSS three-dimensional deformation data, including a slope radar data processing module, a GNSS three-dimensional deformation data processing module, and a joint observation model. The slope radar data processing module is used to obtain slope radar data and perform discrete value elimination processing. The GNSS three-dimensional deformation data processing module is used to obtain GNSS three-dimensional deformation data and perform discrete value elimination processing and spatial interpolation processing. The joint observation model includes a coordinate transformation calculation module. The coordinate transformation calculation module is used to calculate the displacement data of the slope radar data and convert it to the BLH coordinate system. The coordinate transformation calculation module calibrates the slope radar data and the GNSS three-dimensional deformation data to the same BLH coordinate system. The weight iteration optimization expression of the two types of observation data, namely the slope radar data and the GNSS three-dimensional deformation data, in the joint observation model is as follows: i represents the category of observation data, W i ′ represents the weight corresponding to the i-th iteration of the observation data category, W i represents the weight corresponding to the i-th iteration of the observation data category before iteration, until the weight variances of the two types of observation data are approximately equal, that is The joint observation model takes the weights after iterative optimization and performs data fusion on the slope radar data and the GNSS three-dimensional deformation data according to the weights.

[0038] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0039] (1) The present invention combines threshold classification, the characteristic of maximizing the classification interval of the support vector machine, and the decision tree integration ability of the random forest to achieve precise identification and elimination of abnormal points in the data. It shows stronger accuracy and classification efficiency when dealing with complex data, significantly improving the accuracy and reliability of outlier processing; performs three-dimensional spatial modeling and decomposition on the slope radar data, obtains the displacement in the line-of-sight direction of the slope radar and converts it to the BLH coordinate system, and combines the maximum likelihood optimization method to optimize the variogram parameters. By optimizing the key parameters such as nugget value, sill value, and range in the variogram, the accuracy of the Kriging interpolation result is significantly improved.

[0040] (2) The present invention utilizes the data complementarity of the slope radar data and the GNSS three-dimensional deformation data. The slope radar data and the GNSS three-dimensional deformation data respectively constitute a slope radar observation equation and a GNSS observation equation. Based on the two types of observation equations, a joint observation model of multi-source data is jointly constructed, and a weighted residual algorithm based on the least squares method is introduced. During the fusion process, the weights of the slope radar data and the GNSS three-dimensional deformation data are adjusted and the minimization of the weighted residual sum is achieved, finally generating a high-precision three-dimensional deformation fusion result; the present invention combines multi-source observation data, effectively utilizes the complementary advantages and fusion observations of different source data, can obtain high-precision three-dimensional displacement data of the target, and provides reliable technical support for subsequent slope deformation analysis. Brief Description of the Drawings

[0041] Figure 1 It is the method flow chart of the combined observation and fusion method of the present invention;

[0042] Figure 2 It is the schematic diagram of the principle of the combined observation and fusion method in the embodiment;

[0043] Figure 3 It is the graph of the displacement data changing with time in the GNSS three-dimensional deformation data from January to April in a certain mining area in Xinjiang in the embodiment;

[0044] Figure 4 It is the effect diagram of screening and classification of the GNSS three-dimensional deformation data from January to April in a certain mining area in Xinjiang by the support vector machine in the embodiment;

[0045] Figure 5 It is the schematic diagram of the grid point data of the slope radar data of a 100x100 grid of the slope radar in a certain mining area in Xinjiang in the embodiment;

[0046] Figure 6 For Figure 5 It is the schematic diagram of the grid point data after the standardized processing of the slope radar data in

[0047] Figure 7 For Figure 6 It is the schematic diagram of the adjacent value analysis and separation result of the slope radar data in

[0048] Figure 8 For Figure 7 It is the schematic diagram of the random forest separation result of the slope radar data in

[0049] Figure 9 It is the schematic diagram of the principle of the decomposition and conversion of the slope radar data;

[0050] Figure 10 It is the graph of the change of RMSE in the maximum likelihood optimization process of the Kriging parameters in the embodiment;

[0051] Figure 11 It is the schematic diagram of the visualization result after the fusion of part of the slope radar data and part of the GNSS three-dimensional deformation data in the embodiment;

[0052] Figure 12 It is the comparison result graph of the detection accuracy of the slope radar data, the GNSS three-dimensional deformation data and the data after the combined observation and fusion. Detailed Description of the Preferred Embodiment

[0053] The present invention will be further described in detail below with reference to the embodiments:

[0054] Embodiment

[0055] AsFigure 1 , Figure 2 As shown, a joint observation and fusion method for slope radar and GNSS three-dimensional deformation data, the method comprising:

[0056] S1. Obtain slope radar data and perform discrete value rejection processing. In some embodiments, the discrete value rejection processing method for the slope radar data of the present invention is as follows:

[0057] Adopt proximity value analysis to preliminarily detect discrete values, and then use the Euclidean distance to calculate the neighborhood set of each data point and the mean and standard deviation of the line-of-sight direction displacement within the neighborhood. The suspected discrete value judgment is performed according to the following formula (based on the information of the radar center, illumination point, and line-of-sight direction displacement in the slope radar data, and according to the principle that the displacement values of adjacent illumination points are spatially consistent, the proximity value analysis method is used for the preliminary judgment of discrete values):

[0058] |d LOS，j -u N(j) |>k ld ·σ N(j) , where u N(j) , σ N(j) are the mean and variance within the neighborhood respectively, d LOS，j is the displacement of data point j. k ld is a parameter preset based on the characteristics of the slope radar historical data (in the example of a certain mining area in Xinjiang, k ld takes the value of 3).

[0059] The present invention constructs a random forest decision tree model to select the optimal splitting feature and threshold through the Gini index at the splitting node, so as to maximize the reduction of the Gini index after splitting. t L , t R are the normal value node and discrete value node after splitting of the slope radar data respectively, N L , N R are the sample numbers of the normal value and discrete value nodes respectively, N is the total number of samples of the slope radar data, G(t L ) is the Gini index of the normal value, and G(t R ) is the Gini index of the discrete value. The Gini index expression is as follows: pi is the sample proportion of the normal value or discrete value in the current classification of the normal value and discrete value of the slope radar, c is the number of classification categories, and the radar data is classified into normal points and discrete values. In this embodiment, taking the grid point data of the displacement data in the 100x100 grid of the slope radar data in a certain mining area in Xinjiang as Figure 5 shown, calculate the neighborhood mean and variance:

[0060]

[0061] Normalize the slope radar data according to the neighborhood mean and variance, as shown in Figure 6 the following figure. Then, perform proximity value analysis on each data point in the slope radar data (calculate the neighborhood set of each data point and the mean and standard deviation of the line-of-sight direction displacement within the neighborhood using Euclidean distance), and judge whether the line-of-sight direction displacement deviates from the neighborhood statistical value through the mean and variance conditions, initially separating normal values and discrete values, as shown in Figure 7 the following figure. Based on the proximity value analysis, merge the initially marked gross error points with other normal data points into a dataset, add a marking column to the dataset to mark whether it is a gross error marked by the proximity value analysis, 0 for normal data, and 1 for initial gross error. According to the average proximity distance and local distance mean of the radar data, construct a decision tree, and select the optimal splitting feature and threshold through the Gini index.

[0062] The random forest decision tree model determines the discrete values in the slope radar data by majority voting of all decision trees and eliminates them. The majority voting expression is as follows: T = mod(t1, t2,..., t i ), i = 1, 2,..., n, where T is the total judgment of whether the tree in the random forest is a discrete value, and t i is each decision tree. After constructing the optimal separation tree, eliminate the discrete values in the slope radar data according to the criteria of the separation tree. Since the random forest strengthens the identification of discrete values in the slope radar through multiple feature selections and splits, it can accurately extract the discrete values that meet the conditions in the slope radar data, providing an accurate data basis for subsequent function construction. In the example of slope radar data in a certain mining area in Xinjiang, the random forest decision tree model constructs a decision tree, selects the optimal splitting feature and threshold through the Gini index for discrete value separation, and the random forest separation result of the slope radar data is as shown in Figure 8 the following figure.

[0063] Obtain GNSS three-dimensional deformation data for discrete value elimination processing and spatial interpolation processing. In some embodiments, the method for discrete value elimination processing of GNSS three-dimensional deformation data is as follows:

[0064] Normalize the GNSS three-dimensional deformation data, calculate the mean of the GNSS three-dimensional deformation data within the sliding window using the sliding window, and initially determine the data points that deviate from the mean of the sliding window as discrete values. A simple example is as follows: The dataset of GNSS three-dimensional deformation data is represented as X = {(x 1，gnss , y 1，gnss , z l，gnss ), (x 2，gnss , y 2，gnss , z 2，gnss ), (x 3，gnss , y 3，gnss , z 3，gnss)......, perform standardization processing on the data, and then use the sliding window technique to detect outliers in GNSS within consecutive time periods. When setting the sliding window and anomaly detection, considering the overall change of GNSS data, set an appropriate window size and calculate the mean value of GNSS data within the window; judge each GNSS data point, and preliminarily determine the data points that deviate from the mean value of the sliding window as discrete values:

[0065] |X′ i，gnss -u ck |>k gnss σ ck , i = 1, 2...n; X′ i.gnss is the displacement data of GNSS over time after standardization, u ck is the mean value calculated by the window, σ ck is the standard deviation calculated by the window, k gnss is a variable that controls the discrete value threshold (determined according to the characteristics of GNSS historical data). If the above formula conditions are met, it is regarded as a discrete value. In this embodiment, taking the displacement data (selecting displacement among displacement, velocity, and acceleration) in the GNSS three-dimensional deformation data of a certain mining area in Xinjiang as an example, the time range is intercepted from January 1st to April 1st, and Figure 3 the displacement data change graph over time in the GNSS three-dimensional deformation data shown in the figure is obtained. By calculation, u N(j) and σ N(j) are obtained. In the example of a certain mining area in Xinjiang, set the sliding window size to 20 data points, calculate the mean value and variance of the sliding window, and screen out discrete value points.

[0066] Classify the preliminarily separated discrete values through a support vector machine with the goal of maximizing the classification interval. Use the label +1 to represent normal data points and -1 to represent discrete data points; in some embodiments, select the displacement, velocity, and acceleration of the separated data as input values, and then perform support vector machine classification. Map the data to a high-dimensional feature space through a radial basis function kernel, and find a hyperplane in this space to separate normal data points from discrete values. The support vector machine SVM can maximize the interval between the two types of samples while ensuring that the sample points are correctly classified. The classification optimization is expressed as:

[0067] The constraint condition expression of the support vector machine is: d n，svm (w svm T g n，svm +b svm )>1, g n，svm represents the vector combination of displacement, velocity, and acceleration input for each data point, d n，svm represents the category of the samples after rough separation, w svmis the normal vector of the hyperplane and is perpendicular to the decision hyperplane. b svm is the bias term and is used to adjust the position of the decision hyperplane; the constraint conditions of the support vector machine ensure that each data point is classified on one side of the hyperplane and is at least 1 away from the decision hyperplane. In the example of a certain mining area in Xinjiang, the classification effect of GNSS three-dimensional deformation data through the support vector machine is as Figure 4 shown. The black solid line represents the classification hyperplane, the dotted line represents the boundary interval, the shaded area is the interval range, and the points close to the interval boundary are the support vectors.

[0068] Through the combination of the hyperplane normal vector w svm and the bias term b svm , calculate the distance from each data point to the hyperplane. If the distance is less than or equal to -1, then determine that the data point is a discrete value and perform elimination processing.

[0069] For the discrete values that cannot be completely separated initially, introduce a soft margin by introducing a slack variable ξ n，svm , with as the goal, and the constraint condition expression is as follows:

[0070] d n，svm (w svm T g i，svm +b svm )>1 - ξ n，svm , where C svm is a hyperparameter used to control the maximization of the soft margin and the error penalty, and further separate and eliminate the discrete values that cannot be completely separated; in the example of a certain mining area in Xinjiang, C is taken as 10 and ξ is taken as 0.3. In some embodiments, after determining the optimal classification criterion through different combinations of w svm and b svm , perform discrete value detection and elimination. Calculate the distance from each sample point to the hyperplane through the decision function f(x), and determine whether it is a discrete value. If f(x) = -1, then determine that the data point is a discrete value; if f(x) = +1, then it is a normal data point; after the detection is completed, eliminate all the samples that are discrete values and only retain the normal value data to obtain the new dataset after elimination.

[0071] In some embodiments, the spatial interpolation processing of GNSS three-dimensional deformation data in method S1 is performed by the Kriging method using spatial correlation for interpolation processing; the Kriging interpolation formula is as follows:

[0072] where is the predicted value of the unknown position x0, is the value on the known data point x i numbered i, and λ i is the Kriging weight numbered i.

[0073] In a further embodiment, the maximum likelihood method is introduced in the Kriging interpolation process of the Kriging method for parameter optimization. By taking partial derivatives and optimizing the parameters in the negative gradient direction, the optimized interpolation model can effectively process the characteristics of GNSS data, thereby providing more accurate results for GNSS spatial interpolation. When using the Kriging method to interpolate using spatial correlation, it is necessary to model the spatial variability of GNSS data; this variability is quantified by calculating the variogram. For two different station positions (the station positions are the GNSS data station positions) x and x+h, the variogram γ(h) is: where Z (x) and Z (x+h) are the observed values at two stations x and x+h, h is the spatial distance between the positions, and a Gaussian model is used to fit the variogram. The Gaussian model expression is as follows: In the formula, a is the range, c0 is the sill, c1 is the structural part. The three parameters of the range, sill, and structural part are the key parameters affecting the fitting effect of Kriging interpolation. In order to reduce uncertainty and increase the spatial interpolation accuracy of the model, it is necessary to fit the variogram parameters. Set the initial parameter range of θ=(c0, c1, a), and then calculate the covariance under the initial parameters, which describes the correlation between each pair of data points; for two data points x i and x j , the covariance σ ij between them is: where h ij is the spatial distance between data points x i and x j .

[0074] Calculate the likelihood function: where θ=(c0, c1, a) are the Gaussian model parameters, ∑(θ) is the covariance matrix of θ, and z is the vector of GNSS station observed values.

[0075] The simplified calculation of the log-likelihood function is as follows:

[0076]

[0077] Calculate the gradient of the log-likelihood function with respect to the parameters and update the parameter θ k

[0078]

[0079] where θ k is the parameter value at the k-th iteration, a is the learning rate, which controls the update amplitude, is the gradient of the log-likelihood function at θ k .

[0080] In each iteration, by calculating the gradient of the current parameter combination θ k , the gradient of the log-likelihood function with respect to each parameter is obtained, and then the parameters are updated along the direction of the negative gradient, gradually approaching the maximum value of the log-likelihood function to obtain the optimal parameter combination of θ = (c0, c1, a) for Kriging interpolation. is the value of the unknown position x0 to be predicted. is the known data point x i The values on it. The Kriging equations are as follows: In the formula, γ(h i0 ) is the variogram value between the GNSS station x i and the target point x0, γ(h ij ) is the variogram value between the GNSS station x i and the GNSS station x j , and u is the mean value of the interpolation area of the target point. Finally, the Kriging equations are solved to obtain the target three-dimensional deformation result.

[0081] An example of Kriging interpolation of the target point with the coordinates of 10 GNSS stations is as follows:

[0082] Table 1 GNSS Station Coordinates

[0083] Station X(m) Y(m) Z(Displacement) Station_1 3.745401188 0.205844943 3.059264474 Station_2 9.507143064 9.699098522 0.697469303 Station_3 7.319939418 8.324426408 1.460723243 Station_4 5.986584842 2.123391107 1.831809216 Station_5 1.560186404 1.818249672 2.280349921 Station_6 1.559945203 1.834045099 3.925879807 Station_7 0.580836122 3.04242243 0.998368911 Station_8 8.661761458 5.247564316 2.571172192 Station_9 6.011150117 4.319450186 2.962072844 Station_10 7.080725778 2.912291402 0.232252064

[0084] To minimize the error of the Kriging interpolation model, the optimization objective is the variogram parameter θ = (c0, c1, a). The initial parameter ranges are set as nugget value c0: (0, 0.1), sill value c1: (0.5, 2.0), range a: (10, 100). Calculate the covariance matrix and take the partial derivative, and optimize the parameters along the direction of the negative gradient. After the MSE no longer decreases significantly, the optimized Kriging interpolation result is finally obtained.

[0085] Table 2 Iteration Process of (c0, c1, a) and Its MSE Results

[0086] <![CDATA[c0]]> <![CDATA[c1]]> a MSE 0.074 1.4979 24.2182 0.147344 0.0155 1.3784 66.8489 0.09742 0.0971 1.0215 38.7676 0.06729 0.0802 1.486 42.1057 0.05173 0.0151 1.8012 91.1919 0.04319 0.0043 1.3741 77.1295 0.03576 0.0213 1.6152 61.2535 0.01546 0.0649 0.8341 88.5739 0.0168 0.0615 1.6733 12.8962 0.0146 0.0764 1.9569 45.7422 0.0137

[0087] Use the optimal variogram parameter values to perform Kriging interpolation on 6 points.

[0088] Table 3 Kriging Interpolation Results of Target Points

[0089] GNSS Point Number x y z 1 2.324 -4.172 7.847 2 -3.933 1.237 -5.634 3 0.834 -2.534 3.114 4 -7.243 4.617 -8.347 5 5.573 -3.334 1.473 6 -0.113 2.734 6.348

[0090] This embodiment adopts the maximum likelihood optimization interpolation iteration of Kriging parameters. Figure 10The figure shows the results of optimizing the variogram parameters of the Kriging interpolation method using the maximum likelihood method. It can be seen from the figure that as the number of optimization iterations increases, the RMSE gradually decreases, demonstrating that the maximum likelihood optimization can effectively improve the performance of the interpolation model.

[0091] S2. Calculate the displacement data of the slope radar data and convert it to the BLH coordinate system. In some embodiments, the displacement data of the slope radar data of the present invention can be processed by the following three-dimensional displacement conversion relationship. The expression of the three-dimensional displacement conversion relationship is as follows: Where S e , S n , S u are the projection vectors in the longitude, latitude, and elevation directions of the slope radar data respectively, and (B i , L i , H i ) are the coordinates in the geocentric rectangular coordinate system respectively, and (B0, L0, H0) is the BLH coordinate of the radar equipment center point of the slope radar data.

[0092] In some embodiments, the method for converting the displacement data of the slope radar data is as follows:

[0093] First, convert the BLH coordinate (B0, L0, H0) of the radar equipment center point to the ECEF coordinate (X0, Y0, Z0). The expression is as follows:

[0094] Where a is the semi-major axis of the WGS84 ellipsoid, and e 2 is the square of the first eccentricity.

[0095] According to the polar coordinates (r i , θ i , a i ) of the slope radar pixel points, calculate the displacement (x′ i , y′ i , z′ i ) of the pixel points in the ENU coordinate system. The expressions of (x′ i , y′ i , z′ i ) are as follows:

[0096] x′ i = r i cos(a i )sin(θ i )

[0097] y′ i = r i cos(a i )cos(θ i )

[0098] z′i = r i sin(a i )

[0099] Convert ENU coordinates to ECEF coordinates

[0100]

[0101] Convert local coordinates (x′ i , y′ i , z′ i ) to geocentric rectangular coordinates (X i , Y i , Z i ).

[0102] X i = X0 + E i

[0103] Y i = Y0 + N i

[0104] Z i = Z0 + U i

[0105] Based on the geocentric rectangular coordinates (X i , Y i , Z i ), calculate the (B i , L i , H i ) of the slope radar pixel points.

[0106]

[0107] Complete the process of converting slope radar pixel points to the BLH coordinate system. The present invention can construct the following three-dimensional displacement conversion relationship for conversion: Where S e , S n , S u are the projection vectors in the longitude, latitude, and elevation directions respectively. Convert polar coordinates to radar pixel BLH coordinates through polar coordinate parameters, convert the position of slope radar pixel points to the GNSS coordinate system, match the position of GNSS monitoring points with radar pixel points, and establish the spatial correspondence between GNSS and radar images, providing a conversion coefficient for converting the displacement in the line-of-sight direction of slope radar to three-dimensional direction displacement in the subsequent fusion of slope radar and GNSS data. In the example of slope radar data in a certain mining area in Xinjiang, according to Figure 9 Calibration and decomposition conversion principle (by converting the radial distance and angle into X, Y, Z coordinates in three-dimensional space, the spatial distribution and height change of radar signals can be analyzed). The example is as follows:

[0108] Table 4 Example Data of Slope Radar

[0109]

[0110]

[0111] The calculation results through three-dimensional space calibration are as follows in the table:

[0112] Table 5 Three-dimensional Calibration Results of Pixel Points

[0113] Target X Y Z Target 1 -2.50 9.01 4.63 Target 2 1.97 -6.879627191 -6.88 Target 3 -8.83 7.323522915 2.022 Target 4 4.16 -9.588310114 9.39 Target 5 6.64 -5.753217786 -6.36

[0114] S3. Construct a joint observation model for slope radar data and GNSS three-dimensional deformation data in the BLH coordinate system. The weight iteration optimization expressions for the two types of observation data, slope radar data and GNSS three-dimensional deformation data, in the joint observation model are as follows: i represents the type of observation data, and W i ′ represents the weight corresponding to the type of observation data i after iteration, and W i represents the weight corresponding to the type of observation data i before iteration, until the weight variances of the two types of observation data are approximately equal, that is The weight iteration optimization process in this embodiment is as follows in the table:

[0115] Table 6 Iterative Change Process of w1 and w2

[0116]

[0117] Calculate the convergence condition. After reaching the convergence condition, convergence is completed.

[0118] In some embodiments, the joint observation model is solved by minimizing the residual based on the least squares principle, and its result expression is as follows:

[0119] X = (A T PA) -1 A T PL, where A is the coefficient matrix, A T is the transpose of the coefficient matrix, P is the weight matrix, and L is the observation matrix of the joint observation model.

[0120] The joint observation model of the present invention is composed of a slope radar observation equation and a GNSS observation equation.

[0121] The expression of the slope radar observation equation is as follows:

[0122] Among them, V los is the deformation value in the radar line-of-sight direction, S e , S n , S uThey are the projection vectors of the slope radar data in the longitude, latitude, and elevation directions, respectively, D e , D n , D u They are the deformation values in the longitude, latitude, and elevation directions that need to be calculated, respectively, D′ los is the initial displacement in the radar line-of-sight direction, D′ e , D′ n , D′ u They are the initial deformation values in the longitude, latitude, and elevation directions, respectively.

[0123] The GNSS observation equation expression is as follows:

[0124] where L e，gnss , L n，gnss , L u，gnss They are the deformation values in the longitude, latitude, and elevation directions obtained by decomposing the GNSS three-dimensional deformation data, respectively, D e , D n , D u They are the deformation values in the longitude, latitude, and elevation directions that need to be calculated, respectively, V e，gnss , V n，gnss , V u，gnss They are the residuals in the east-west, north-south, and vertical directions, respectively.

[0125] S4. The joint observation model takes the weights after iterative optimization and performs data fusion on the slope radar data and the GNSS three-dimensional deformation data according to the weights. The observation matrix L expression of the joint observation model combining the slope radar observation equation and the GNSS observation equation is as follows:

[0126]

[0127] is the fusion result X of the joint observation model.

[0128] In this embodiment, taking a certain mining area in Xinjiang as an example, the visualization results after fusing some slope radar data and some GNSS three-dimensional deformation data are as Figure 11 shown, Figure 11 which shows the visualization results of the three-dimensional displacement field after completing the fusion of GNSS and slope radar data. In the example of a certain mining area in Xinjiang, the comparison results of the detection accuracies of the slope radar data, the GNSS three-dimensional deformation data, and the data after joint observation fusion are as Figure 12 shown. It can be seen from the figure that the RMSE of the data after the joint observation fusion of the present invention is significantly better than that of using only the GNSS three-dimensional deformation data or the slope radar data, indicating that compared with the GNSS three-dimensional deformation data or the slope radar data, the data after the joint observation fusion can more accurately reflect the three-dimensional deformation field of the slope and improve the monitoring quality.

[0129] A joint observation and fusion system for slope radar and GNSS three-dimensional deformation data, comprising a slope radar data processing module, a GNSS three-dimensional deformation data processing module, and a joint observation model. The slope radar data processing module is used to obtain slope radar data and perform discrete value elimination processing. The GNSS three-dimensional deformation data processing module is used to obtain GNSS three-dimensional deformation data and perform discrete value elimination processing and spatial interpolation processing. The joint observation model includes a coordinate transformation calculation module. The coordinate transformation calculation module is used to calculate the displacement data of the slope radar data and transform it to the BLH coordinate system. The coordinate transformation calculation module calibrates the slope radar data and the GNSS three-dimensional deformation data to the same BLH coordinate system. The weight iterative optimization expression of the two types of observation data, namely the slope radar data and the GNSS three-dimensional deformation data, in the joint observation model is as follows: i represents the category of observation data, W i ′ represents the weight corresponding to the observation data category i after iteration, W i represents the weight corresponding to the observation data category i before iteration, until the weight variances of the two types of observation data are approximately equal, that is The joint observation model takes the weights after iterative optimization and performs data fusion on the slope radar data and the GNSS three-dimensional deformation data according to the weights.

[0130] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, and improvements made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A combined observation and fusion method for 3D deformation data of slope radar and GNSS, characterized in that: The method includes: S1. Obtain slope radar data and perform discrete value rejection processing, and obtain GNSS three-dimensional deformation data and perform discrete value rejection processing and spatial interpolation processing; S2. Calculate the displacement data of the slope radar data and convert it to the BLH coordinate system; S3. Construct a joint observation model for slope radar data and GNSS three-dimensional deformation data in the same BLH coordinate system. The weight iteration optimization expressions for the two types of observation data, i.e., slope radar data and GNSS three-dimensional deformation data, in the joint observation model are as follows: i represents the category of observation data, and W′ i represents the weight corresponding to the category i of observation data after iteration, and W i represents the weight corresponding to the category i of observation data before iteration, until the weight variances of the two types of observation data are approximately equal, i.e., S4. The combined observation model takes the iteratively optimized weight and performs data fusion of the slope radar data and the GNSS three-dimensional deformation data according to the weight.

2. The joint observation and fusion method of slope radar and GNSS three-dimensional deformation data according to claim 1, wherein: In method S3, the combined observation model solves the residual minimization processing based on the least squares principle, and its result expression is as follows: X = (A T PA) -1 A T PL, where A is the coefficient matrix, A T is the transpose of the coefficient matrix, P is the weight matrix, and L is the observation matrix of the joint observation model.

3. The joint observation and fusion method of slope radar and GNSS three-dimensional deformation data according to claim 2, characterized in that: The combined observation model is composed of a slope radar observation equation and a GNSS observation equation. The expression of the slope radar observation equation is as follows: Among them, V los is the deformation value in the radar line-of-sight direction, S e , S n , S u are the projection vectors in the longitude, latitude, and elevation directions of the slope radar data respectively, D d , D n , D u are the deformation values in the longitude, latitude, and elevation directions to be calculated respectively, D′ los is the initial displacement in the radar line-of-sight direction, D′ e , D′ n , D′ u are the initial deformation values in the longitude, latitude, and elevation directions respectively; The expression of the GNSS observation equation is as follows: Among them, L e，gnss , L n，gnss , L u，gnss are the deformation values in the longitude, latitude, and elevation directions of the GNSS three-dimensional deformation data decomposition respectively. D e , D n , D u are the deformation values in the longitude, latitude, and elevation directions that need to be calculated respectively. v e，gnss , v n，gnss , v u，gnss are the residuals in the east-west, north-south, and vertical directions respectively; The expression of the observation matrix L of the combined observation model combining the slope radar observation equation and the GNSS observation equation is as follows: It is the fusion result X of the joint observation model.

4. The joint observation and fusion method of slope radar and GNSS three-dimensional deformation data according to claim 1, wherein: In method S1, the method for rejecting discrete values of GNSS three-dimensional deformation data is as follows: Perform standardization processing on the GNSS three-dimensional deformation data, calculate the mean value of the GNSS three-dimensional deformation data within the sliding window using the sliding window, and initially determine the data points that deviate from the mean value of the sliding window as discrete values; Classify the initially separated discrete values through a support vector machine with the goal of maximizing the classification interval, use the label +1 to represent normal data points, and -1 to represent discrete data points; The constraint expression of the support vector machine is: d n，svm (w svm T g n，svm +b svm )>1, g n，svm represents the combination of displacement, velocity, and acceleration vectors of each data point input, d n，svm represents the category of the samples after rough separation, w svm is the normal vector of the hyperplane, perpendicular to the decision hyperplane; b svm is the bias term and is used to adjust the position of the decision hyperplane; the constraint condition of the support vector machine classifies each data point on one side of the hyperplane and at least a distance of 1 from the decision hyperplane; Through the hyperplane normal vector w svm and the bias term b svm in combination, calculate the distance from each data point to the hyperplane. If the distance is less than or equal to -1, determine that the data point is a discrete value and perform rejection processing.

5. The combined observation and fusion method for slope radar and GNSS three-dimensional deformation data according to claim 4, characterized in that: For the discrete values obtained by the initial separation that cannot be completely separated, a soft margin is introduced by introducing slack variables ξ n，svm , with as the objective, and the constraint expression is as follows: d n，svm (w svm Tg i，svm +b svm ) > 1 - ξ n，svm where C svm is a hyperparameter used to control the maximization of the soft margin and the error penalty, and further separate and eliminate discrete values that cannot be completely separated.

6. The combined observation and fusion method of slope radar and GNSS three-dimensional deformation data according to claim 1, characterized in that: The method for rejecting discrete values of slope radar data is as follows: Use proximity value analysis to initially detect discrete values, then calculate the neighborhood set of each data point and the mean and standard deviation of the line-of-sight direction displacement within the neighborhood using the Euclidean distance, and judge the suspected discrete values according to the following formula: |d LOS,j -u N(j) | > k ld ·σ N(j) , where u N(j) , σ N(j) are the mean and variance within the field respectively, and d LOS，j is the displacement of data point j; k ld is a parameter preset based on the characteristics of the historical data of the slope radar. The random forest decision tree model is constructed by splitting nodes to select the optimal splitting feature and threshold through the Gini index, achieving the maximization of the reduction of the Gini index after splitting. t L and t R are the normal value node and the discrete value node after splitting the slope radar data respectively. N L and N R are the sample numbers of the normal value node and the discrete value node respectively. N is the total number of samples of the slope radar data. G(t L ) is the Gini index of the normal value, and G(t R ) is the Gini index of the discrete value. The random forest decision tree model determines and rejects the discrete values in the slope radar data by the majority vote of all decision trees.

7. The combined observation and fusion method for 3D deformation data of slope radar and GNSS according to claim 1, characterized in that: In method S1, the spatial interpolation processing of the GNSS three-dimensional deformation data is performed by the Kriging method using spatial correlation for interpolation processing, and the Kriging interpolation formula is as follows: wherein is the predicted value of the unknown bit x0, is the value of the known data point x with the number i i at, and λ i is the Kriging weight with the number i.

8. The joint observation and fusion method of slope radar and GNSS three-dimensional deformation data according to claim 1, characterized in that: In method S2, the expression of the three-dimensional displacement conversion relationship is as follows: where S e , S n , S u are the projection vectors in the longitude, latitude, and elevation directions of the slope radar data respectively. (B i , L i , H i ) are the coordinates in the geocentric rectangular coordinate system respectively, and (B0, L0, H0) is the BLH coordinate of the center point of the radar device of the slope radar data.

9. A combined observation and fusion system for 3D deformation data of slope radar and GNSS, characterized in that: It includes a slope radar data processing module, a GNSS three-dimensional deformation data processing module, and a joint observation model. The slope radar data processing module is used to obtain slope radar data and perform discrete value rejection processing. The GNSS three-dimensional deformation data processing module is used to obtain GNSS three-dimensional deformation data and perform discrete value rejection processing and spatial interpolation processing. The joint observation model includes a coordinate transformation calculation module. The coordinate transformation calculation module is used to calculate the displacement data of the slope radar data and transform it to the BLH coordinate system. The coordinate transformation calculation module calibrates the slope radar data and the GNSS three-dimensional deformation data to the same BLH coordinate system. The weight iterative optimization expression of the two types of observation data, namely the slope radar data and the GNSS three-dimensional deformation data, in the joint observation model is as follows: i represents the category of observation data, W′ i represents the weight corresponding to the i-th iteration of the observation data category, W i represents the weight corresponding to the i-th iteration of the observation data category before iteration, until the weight variances of the two types of observation data are approximately equal, that is The joint observation model takes the weights after iterative optimization and performs data fusion on the slope radar data and the GNSS three-dimensional deformation data according to the weights.

Images