GNSS Station Layout Method Based on InSAR Deformation and Its Application

Abstract

The present invention provides a GNSS station layout method based on InSAR deformation and its application. Based on InSAR historical deformation, through spatial sampling and the Kriging interpolation method in geostatistics, under the premise of taking into account cost and accuracy, the reference number, spatial distribution and corresponding deformation prediction accuracy of the GNSS station layout in the deformation area are given through multi-scale iteration. This can not only enable the GNSS station to capture the spatial characteristics of the deformation to the greatest extent possible, but also improve the accuracy of the subsequent InSAR and GNSS fusion. The present invention fully utilizes the advantages of InSAR in deformation monitoring, such as large scale, high precision and high spatial resolution, provides guidance and reliable basis for the layout of stations during GNSS deformation monitoring, and helps to better capture the spatial characteristics of regional deformation. In the current context of frequent geological disasters, it is very beneficial to the monitoring and understanding of geological disasters.

Description

Technical Field

[0001] The present invention relates to synthetic aperture radar interferometry (InSAR) technology and global navigation satellite system (GNSS) technology in the field of geodesy based on remote sensing images, and particularly relates to a GNSS station layout method based on InSAR deformation and its application. Background Art

[0002] The global navigation satellite system (GNSS), as one of the most mature precision positioning technologies at present, has the advantages of all-weather, automation and long-term continuous observation. Therefore, it has been widely used in many engineering surveying and natural disaster monitoring fields. However, the cost of GNSS is relatively high, and it is impossible to achieve large-area and high-density layout. Therefore, the spatial resolution of its monitoring results is relatively low, and it cannot reflect the spatial variation of deformation in the monitoring area. The synthetic aperture radar interferometry (InSAR) technology has achieved great breakthroughs and achievements in various fields such as volcanoes, earthquakes, mining areas, glaciers, atmosphere, permafrost, etc. with its unique advantages after decades of development, and it is one of the effective means for current geological disaster monitoring. And it has the advantages of high spatial resolution, high precision and large-scale rapid monitoring. However, the temporal resolution of InSAR monitoring results is limited by the revisit period of the sensor. Even when using the combined observation of constellations A and B of Sentinel-1 satellite, its revisit period is only six days, which still cannot fully and correctly reflect the deformation characteristics of sudden, non-linear and large deformation gradient such as mining areas, landslides, and groundwater.

[0003] At present, many studies have been dedicated to the integration of InSAR and GNSS technologies, such as using GNSS to correct the atmospheric delay in InSAR, jointly solving three-dimensional deformation with GNSS and InSAR, and fusing the deformation results of InSAR and GNSS to obtain surface deformation results with high spatio-temporal resolution. However, these studies have hardly considered how to deploy GNSS points to achieve the best effect of InSAR and GNSS integration, and how to balance the relationship between cost and monitoring accuracy. In particular, deformation information is the most intuitive precursor reflection of geological disasters. However, the station layout for traditional GNSS deformation monitoring mainly relies on empirical guidance and lacks theoretical basis. Therefore, the obtained GNSS layout scheme is only a feasible one, resulting in certain cost waste and low deformation monitoring accuracy. At the same time, its deformation monitoring results are also difficult to reflect the spatial distribution characteristics of deformation in the entire study area. Therefore, how to consider the rationality of the spatial distribution and the number of GNSS deployments to maximize the integration accuracy of InSAR and GNSS deformation while minimizing the cost is an application problem that urgently needs to be solved at present. Summary of the Invention

[0004] The purpose of the present invention is to overcome the defect that the station layout lacks guidance and reliable basis in current GNSS deformation monitoring, and to propose a GNSS station layout method and its application based on InSAR deformation. This method is based on InSAR historical deformation, and through spatial sampling and Kriging interpolation method in geostatistics, on the premise of considering cost and accuracy, the reference number, spatial distribution of GNSS points in the deformation area and their corresponding deformation prediction accuracy are given iteratively, which can not only make the GNSS points capture deformation characteristics as much as possible, but also improve the accuracy of the subsequent integration of InSAR and GNSS.

[0005] The present invention first provides the following technical solutions:

[0006] A GNSS station layout method based on InSAR deformation, which includes the following steps:

[0007] S1. InSAR data processing: Obtain the differential interferogram atlas of radar images of a single or multiple orbits, and calculate the one-dimensional / two-dimensional / three-dimensional average surface deformation rate within a period of time through InSAR technology.

[0008] S2. Determine the candidate point set: Use all InSAR deformation points as the initial point set for GNSS station layout. Through quadtree downsampling, select GNSS station candidate points that can reflect the deformation spatial characteristics in all deformation directions to reduce the quantity and density of the initial point set. At the same time, obtain the mask area where stations cannot be deployed according to external data such as the geological structure and topography of the study area, and remove the candidate points in the mask area to obtain the final candidate point set.

[0009] S3. Establish the objective function model: First, estimate the variogram model of the InSAR deformation rate field. Meanwhile, with the minimum prediction variance as the goal, interpolate the entire deformation field through the Kriging method according to the selected variogram and measurement stations, and establish a model that minimizes the root mean square error (RMSE) between the interpolated deformation and the InSAR-monitored deformation.

[0010] S4. Determine the initial measurement station distribution: The initial measurement station distribution includes two types of measurement stations. One is the measurement stations that must be arranged according to the needs of the project and the item, and these measurement stations will not change subsequently; the other is the measurement stations that need to be iteratively optimized using the objective function, and the initial positions of these measurement stations are determined by sequentially deleting candidate points according to the change of the root mean square.

[0011] S5. Solve the objective function model by multi-scale iterative optimization: Solve the objective function model using the multi-scale iterative optimization method according to the InSAR deformation, the initial measurement stations, the candidate point set, and the variogram of the deformation field until the accuracy of the model solution meets the requirements.

[0012] According to some preferred embodiments of the present invention, the S1 includes:

[0013] S11. For the obtained radar images, after registration, selection of small baseline networking, interference, removal of flat ground and terrain, phase filtering, phase unwrapping, atmospheric correction, and geocoding, an interference map set is formed.

[0014] S12. Assume that there is only one spatio-temporal baseline set, then the time series deformation and the average deformation rate can be solved by the least squares method through the Small BAseline Subsets (SBAS) method.

[0015] S13. According to the data situation of the obtained research area, one-dimensional, two-dimensional, or three-dimensional deformations can be further obtained respectively:

[0016] One-dimensional deformation: If there is only one orbit, the Line of Sight (LOS) deformation monitored by SBAS-InSAR can be directly used, or it can be converted into the vertical deformation:

[0017] D U =D LOS / cosθ inc

[0018] Where D U represents the vertical deformation, D LOS is the LOS deformation monitored by SBAS-InSAR, and θ inc represents the local radar incident angle.

[0019] Two-dimensional deformation: If ascending and descending orbit data can be obtained, the vertical and east-west two-dimensional deformations can be obtained by fusing the ascending and descending orbits and ignoring the north-south deformation:

[0020]

[0021]

[0022] where D U and D E are the vertical and east-west deformations to be solved respectively; represents the deformation monitored by the ascending orbit data, is the local incidence angle of the radar of the ascending orbit satellite, is the radar azimuth angle of the ascending orbit satellite; represents the deformation monitored by the descending orbit data, is the local incidence angle of the radar of the descending orbit satellite, is the radar azimuth angle of the descending orbit satellite.

[0023] The parameters in the above formula can be directly solved by the least squares method.

[0024] Three-dimensional deformation: If ascending and descending orbits of different sensors can be obtained, classical methods such as D-InSAR, Offset-Tracking, Multi-Aperture InSAR, etc. are used to obtain the deformations in multiple directions, and finally the east-west, north-south, and vertical deformations are solved by the least squares method.

[0025] According to some preferred embodiments of the present invention, the S2 includes:

[0026] S21. Take all InSAR deformation points as initial candidate points, respectively obtain candidate point sets in different directions through quadtree downsampling, and then take the union of the candidate point sets in different directions to obtain the candidate point set of the GNSS station.

[0027] S211. Take the InSAR deformation point set as the initial candidate point (this InSAR deformation point set can be either in raster data format or vector data format);

[0028] S212. Considering that the variance can reflect the spatial variation characteristics of the deformation, a variance threshold T of the deformation is set based on the InSAR deformation; [[ID=4G]]

[0029] S213. Use the quadtree recursive algorithm to continuously subdivide the InSAR point set into squares, and calculate the variance of the InSAR deformation in each square:

[0030]

[0031] Among them, σ 2 represents the variance of InSAR points within the square, n represents the number of InSAR points within the square, and d i represents the deformation of the i-th point within the square, is the mean value of n InSAR deformation points within the square.

[0032] If the variance is less than or equal to T, stop the segmentation; if the variance is greater than T, continue the segmentation until the variance of the deformation in each square is less than or equal to T. Use the point set after quadtree downsampling as the candidate point set in this direction.

[0033] S214. Find the union of the candidate point sets in multiple directions to obtain the initial candidate point set for the entire region.

[0034] S22. Obtain external data such as available geological structures and topographic features in the study area, and mask the areas where it is not suitable to set up measurement stations, such as areas with too large DEM gradients and overly lush vegetation areas;

[0035] S23. Remove the candidate points in the S24 masked area from the initial candidate point set to obtain the final candidate point set.

[0036] According to some preferred embodiments of the present invention, the S3 includes:

[0037] S31. Calculate the variation value of the two-dimensional deformation field using all InSAR deformation points (when the data volume is large, random sampling can be used to reduce the data volume):

[0038]

[0039] In the formula, h is the distance between deformation points, also known as the potential difference; N(h) is the number of pairs of all observation points with a distance of h; z(x i ) and z(x i +h) respectively represent the InSAR deformation observation values with a relative distance of h; is the variance at a distance of h, that is, the variation value, which increases with the increase of h within a certain range. When the distance of the measured points is greater than the maximum correlation distance, this value tends to be stable. Through this step, the variance of the InSAR deformation field and the sample values of the distance h can be obtained.

[0040] The goal is to establish the relationship between the variance and the spatial distance h between point pairs, that is, the variogram. Since there are many classical models in the stable process of the variogram, such as the exponential model, Gaussian model, spherical model, etc. According to the above-obtained sample values, select the model with the best fitting degree as the variogram and fit and calculate the model parameters.

[0041] S32. Assume that a total of n candidate points are obtained in a certain area D in step S2, and the set D n is used to represent them. Now, it is necessary to select N points from these n candidate points to deploy GNSS stations. Then, the following sets and symbols are defined:

[0042] Use the set P N to represent the set where these GNSS points are distributed; use (x i , y i , v i ) to represent the position and its deformation rate of the i-th point in the current point set, where i = 1, 2, … N, which can be simplified to

[0043] Use S M to represent a series of P N with the same number of GNSS stations but different distributions, where M represents that there are M sets of P N in this set. Then Then where and both represent the k-th point set.

[0044] S33. Considering that Kriging interpolation can perform unbiased optimal estimation on regionalized variables, therefore, Kriging interpolation is respectively performed on the above M sets of P N to obtain the complete regional deformation field after interpolation. Then, it is compared with the deformation field monitored by InSAR, and the root mean square error RMSE of each GNSS distribution set is calculated respectively:

[0045]

[0046] where d i represents the deformation of the i-th point in the deformation field monitored by InSAR; represents the deformation of the i-th point of Kriging interpolation; N(D) represents the number of InSAR points in the monitoring area D; RMSE k represents the root mean square error of the k-th point set in the set S M .

[0047] S34. The deployment of GNSS stations involves the coordinate positions of each station. Therefore, the objective function should be a function of the position of each station. Taking the RMSE calculated in step S33 as an index, the problem of solving the position of the GNSS station deployment can be transformed into the following problem of minimizing RMSE:

[0048]

[0049] where W is the weight or scale factor, and R represents the root mean square error of different GNSS station deployments.

[0050] More generally, when the two-dimensional or three-dimensional deformation of the study area is conditionally obtained, the deformation in each dimension should be taken into account simultaneously. At this time, the above formula is converted into the following form:

[0051]

[0052] where Ndim = 1, 2, 3, representing the one-dimensional, two-dimensional, and three-dimensional deformations of the study area obtained respectively; R i represents the root mean square error in the i-th direction; W i represents the weight in the i-th direction.

[0053] The weight can adjust the proportion of deformations in different dimensions. For example, when focusing on settlement, the weight of settlement can be increased, and the distribution result of GNSS station layout will more reflect the spatial distribution of settlement. At the same time, in order to reflect the binding force of the weight, it should meet the following requirements:

[0054]

[0055] According to some preferred embodiments of the present invention, the S4 includes:

[0056] S41. Considering the trade-off between the model solution efficiency and accuracy in step S34, first determine an initial point distribution as close as possible to the optimal point distribution, and then iteratively solve the initial points. Therefore, this step is to determine the initial point distribution.

[0057] S411. First, make the following definitions:

[0058] Delete one point in turn from n candidate points, and the remaining all points form a point set Then

[0059]

[0060] where i = 1, 2,..., n, indicating that there are n such point sets; n - 1 indicates that there are n - 1 remaining points in the point set; k ≠ i indicates that the i-th point among the original n candidate points is deleted. Then use these n Q point sets to form a new set F n , then

[0061]

[0062] where n represents the number of elements in the set F n There are.

[0063] S412. First, perform Kriging interpolation on these n candidate points using the variogram obtained in step S31 to obtain the interpolated deformation field. Then, calculate its root mean square error (RMSE) with the InSAR deformation monitoring according to the formula for calculating RMSE in S33, denoted as r′0.

[0064] S413. For each point set in the set F n perform Kriging interpolation using the variogram obtained in step S31 to obtain the interpolated deformation field. Then, calculate its RMSE with the InSAR deformation monitoring respectively according to the formula for calculating RMSE in S33, denoted as r′ n , where n represents that there are n RMSE results.

[0065] S414. Calculate

[0066] r′ = r′ n - r′0

[0067] Record the point set corresponding to the minimum value in r′ and the candidate point (x i′ , y i′ , v i′ ) that is deleted. This minimum value indicates that deleting this point has the least impact on the overall interpolation result among all candidate points. Therefore, this point can be removed, and the point set is used as the new candidate point set.

[0068] S415. Let n = n - 1, and repeat steps S412 - S414 until n = 3 (because Kriging interpolation requires at least 3 points).

[0069] S416. Finally, obtain the order of deleting each candidate point from the initial n candidate points. When the number of GNSS stations N in the given study area is known, delete the first n - N points, and the point set Q′ N composed of the remaining N points is the initial station layout for this area.

[0070] S42. Determine the point set Q M composed of the station positions that must be set up due to engineering and project requirements. The positions of these stations will not change subsequently.

[0071] S43. To ensure that the number of stations is always N, it is necessary to find the point closest to each point in Q N in the initial point set Q′ M respectively, and replace it with the position of the point in Q M . Thus, the initial station Q N is obtained.

[0072] According to some preferred embodiments of the present invention, the S5 includes:

[0073] S51. The initial measurement station Q determined through steps S2 to S4 N It still needs to be further iteratively optimized to fit a higher deformation field as much as possible. To improve the computational efficiency in a large-scale scenario, we propose a multi-scale iterative optimization method for the final solution. First, downscale the initial candidate points to reduce the number of candidate points: increase the variance threshold in step S2 to a certain value, and at this time, a new candidate point set D' is obtained n , and at this time D' n has a much lower density and number than D n ;

[0074] S52. Iteratively solve the initial measurement station Q N The better solution X' under the small-scale candidate point set D' n . N .

[0075] S521. For each point in the initial measurement station Q N except the fixed point Q M (assuming the current point is ), search for points in the small-scale candidate point set D' n whose distance is less than the large distance threshold T1 (such as twice the range of the variogram), and form a point set C n1 , where n1 represents that there are n1 points in this point set.

[0076] S522. First, perform Kriging interpolation on these n candidate points in the initial measurement station Q N using the variogram obtained in step S31 to obtain its interpolated deformation field, and then calculate its root mean square error with InSAR deformation monitoring according to the calculation formula of the root mean square error in S33, denoted as r″0.

[0077] S523. Replace the current point n1 with each point in C and perform Kriging interpolation using the variogram obtained in step S31 respectively to obtain its interpolated deformation field, and then calculate its root mean square error r″ n with InSAR deformation monitoring according to the calculation formula of the root mean square error in S33 respectively.

[0078] S524. Calculate

[0079] r″ = r″ n - r″0

[0080] Record the minimum value r″ in r″ min and its corresponding j represents the jth point in C corresponding to this minimum value. When r″ n1 min ​< 0, indicating that the point is better than the current point so the point is used to replace the previous point Otherwise, the current point remains unchanged, thus completing this iteration.

[0081] S525. Let n = n + 1 (skip when encountering the fixed point Q M ), and repeat steps S523 - S524 until all points are iterated and the initial station is updated to obtain the updated station X' N , and let Q N = X' N .

[0082] S526. Repeat steps S521 - S525 until the distance between the corresponding positions of each point in the previous and current stations Q N is less than a given threshold. At this time, the final X' N is the better solution of the initial station Q N in the small - scale candidate point set D' n .

[0083] S53. Iteratively solve the final solution X N of X' n under the large - scale candidate point set D N . Similar to the process of solving the station X' N , but the initial station Q N needs to be replaced by X' N , the small - scale candidate point set D' n needs to be replaced by the large - scale candidate point set D n , and the larger distance threshold T1 is replaced by the smaller distance threshold T2 (such as half of the range of the variogram). Then repeat steps S521 - S526 to solve the final optimized GNSS station layout result X N .

[0084] According to the above - mentioned method for determining GNSS stations based on InSAR deformation with multi - scale iteration optimization, a system for determining GNSS station layout based on InSAR deformation can be obtained.

[0085] The above - mentioned GNSS station layout method or system can be used in the layout of stations for GNSS deformation monitoring.

[0086] The GNSS station layout method based on InSAR deformation of the present invention makes full use of the advantages of InSAR in large-scale, high-precision, and high-spatial-resolution deformation monitoring. For a certain deformation field monitored by InSAR, it can quickly give the reference number and location of GNSS stations to be arranged, can capture deformation characteristics to the greatest extent and save costs, and can also be used to improve the accuracy of the later InSAR and GNSS fusion. In the current context of frequent geological disasters, it is very beneficial to the monitoring and understanding of geological disasters. Description of the Drawings

[0087] Figure 1 It is a schematic diagram of a specific implementation process of the method of the present invention;

[0088] Figure 2 It is the deformation rate map obtained by SBAS-InSAR in Example 1;

[0089] Figure 3 It is the distribution map of candidate stations in Example 1;

[0090] Figure 4a It is the initial station distribution in Example 1;

[0091] Figure 4b It is the deformation result map predicted by using the initial stations in Example 1;

[0092] Figure 5a It is the final station distribution map in Example 1;

[0093] Figure 5b It is the deformation result map predicted by using the final stations in Example 1. Detailed Implementation Modes

[0094] The present invention will be described in detail below in conjunction with the embodiments and the drawings. However, it should be understood that the embodiments and the drawings are only used for exemplary description of the present invention, and cannot constitute any limitation to the protection scope of the present invention. All reasonable transformations and combinations within the scope of the inventive concept of the present invention fall within the protection scope of the present invention.

[0095] As Figure 1 shown, a specific GNSS station layout method based on InSAR deformation includes the following steps:

[0096] S1. InSAR data processing: Obtain the differential interferogram set of radar images of a single or multiple orbits, and calculate the one-dimensional / two-dimensional / three-dimensional average surface deformation rate over a period of time through InSAR technology.

[0097] In some specific embodiments, step S1 may further include:

[0098] S11. For the acquired radar images, first use orbit information for coarse alignment, then use intensity cross-correlation for fine alignment; set a small baseline threshold to select interferometric pairs that meet the temporal and spatial baseline thresholds (the temporal baseline is 200 days and the spatial baseline is 800 meters) for interferometry, flatland and terrain removal, and phase filtering; select unwrapping reference points for phase unwrapping, and then perform atmospheric correction and geocoding.

[0099] S12. Assuming there is only one spatiotemporal baseline set, establish a deformation rate and terrain residual model, and obtain the surface deformation through least squares solution.

[0100] S13. Based on the data of the study area, one-dimensional, two-dimensional or three-dimensional deformation can be further obtained:

[0101] One-dimensional deformation: If there is only one track, the Line of Sight (LOS) deformation monitored by SBAS-InSAR can be used directly, or it can be converted into vertical deformation:

[0102] D U =D LOS / cosθ inc

[0103] Among them D U Indicates vertical deformation, D LOS is the LOS deformation monitored by SBAS-InSAR, θ inc represents the radar local incidence angle.

[0104] Two-dimensional deformation: If ascending and descending orbit data are available, vertical and east-west two-dimensional deformation can be obtained by fusing the ascending and descending orbits and ignoring the north-south deformation:

[0105]

[0106]

[0107] Among them D U and D E are the vertical and east-west deformations to be determined respectively; represents the deformation monitored by orbit raising data, is the radar local incidence angle of the orbit-raising satellite, is the radar azimuth of the orbit-raising satellite; represents the deformation monitored by the orbit-dropping data, is the radar local incidence angle of the descending satellite, is the radar azimuth of the descending satellite.

[0108] The parameters in the above formula can be directly solved by least squares.

[0109] Three-dimensional deformation: If the ascending and descending orbits of different sensors can be obtained, classical methods such as D-InSAR, Offset-Tracking, and Multi-Aperture InSAR are used to obtain the deformations in multiple directions, and finally the deformations in the east-west, north-south, and vertical directions are solved by least squares.

[0110] S2. Determine the candidate point set: Use all InSAR deformation points as the initial point set for GNSS station layout. Through quadtree downsampling, select GNSS station candidate points that can reflect the deformation space characteristics in all deformation directions to reduce the quantity and density of the initial point set. At the same time, obtain the mask area where stations cannot be laid out based on external data such as the geological structure and topography of the study area, and remove the candidate points in the mask area to obtain the final candidate point set.

[0111] In some specific embodiments, step S2 may further include:

[0112] S21. Use all InSAR deformation points as the initial candidates, and respectively obtain candidate point sets in different directions through quadtree downsampling. Then, take the union of the candidate point sets in different directions to obtain the candidate point set for GNSS stations.

[0113] S211. Use the InSAR deformation point set as the initial candidate (this InSAR deformation point set can be in either raster data format or vector data format);

[0114] S212. Considering that variance can reflect the spatial variation characteristics of deformation, a variance threshold T = 3.5 mm for deformation is set based on InSAR deformation;

[0115] S213. Use the quadtree recursive algorithm to continuously divide the InSAR point set into squares, and calculate the variance of the InSAR deformation in each square:

[0116]

[0117] where σ 2 represents the variance of the InSAR points within the square, n represents the number of InSAR points within the square, d i represents the deformation of the i-th point within the square, is the mean of the n InSAR deformation points within the square.

[0118] If the variance is less than or equal to T, stop the division; if the variance is greater than T, continue the division until the variance of the deformation in each square is less than or equal to T. Use the point set obtained by quadtree downsampling as the candidate point set for this direction.

[0119] S214. Obtain the candidate point set for the entire region by taking the union of the candidate point sets in multiple directions.

[0120] S22. Obtain external data such as available geological structures and topographic features in the study area, and mask out areas that are not suitable for setting up survey stations, such as areas with too large DEM gradients or overly lush vegetation.

[0121] S23. Remove the candidate points in the masked area of S24 from the initial candidate point set to obtain the final candidate point set.

[0122] S3. Establish an objective function model: First, estimate the variogram model of the InSAR deformation rate field. At the same time, with the goal of minimizing the prediction variance, interpolate the entire deformation field through the Kriging method according to the selected variogram and survey stations, and establish a model that minimizes the root mean square error (RMSE) between the interpolated deformation and the InSAR monitored deformation.

[0123] In some specific embodiments, step S3 may further include:

[0124] S31. Calculate the variation value of the two-dimensional deformation field using all InSAR deformation points (when the data volume is large, random sampling can be used to reduce the data volume):

[0125]

[0126] In the formula, h is the spacing between deformation points, also known as the displacement difference; N(h) is the number of pairs of all observation points with a spacing of h; z(x i ) and z(x i + h) respectively represent the InSAR deformation observation values with a relative distance of h; is the variance at a spacing of h, that is, the variation value, which increases with the increase of h within a certain range. When the distance between the measured points is greater than the maximum correlation distance, this value tends to be stable. Through this step, the variance of the InSAR deformation field and the sample values of the spacing h can be obtained.

[0127] This goal is to establish the relationship between the variance and the spatial distance h between point pairs, that is, the variogram. Since there are many classical models in the stable process of the variogram, such as the exponential model, Gaussian model, spherical model, etc. According to the above obtained sample values, select the model with the best fitting degree as the variogram and fit and calculate the model parameters.

[0128] S32. Assume that a total of n candidate points are obtained in a certain area D in step S2, and the set D n is used to represent it. Now it is necessary to select N points from these n candidate points to set up GNSS survey stations, then define the following sets and symbols:

[0129] Use the set \(P\) N to represent the set of the distributions of these GNSS points; use \((x\) i , y\) i , v\) i ) to represent the position and its deformation rate of the \(i\)-th point in the current point set, where \(i = 1, 2, \cdots, N\), which can be simplified to

[0130] Use \(S\) M to represent a series of \(P\)'s with the same number of GNSS stations but different distributions N , where \(M\) represents that there are \(M\) sets of \(P\) in this set N , then Then where and both represent the \(k\)-th point set.

[0131] S33. Considering that Kriging interpolation can perform unbiased optimal estimation on regionalized variables, therefore, perform Kriging interpolation on the above \(M\) sets of \(P\) N respectively to obtain the complete regional deformation field after interpolation. Then compare it with the deformation field monitored by InSAR, and calculate the root mean square error RMSE of each GNSS distribution set respectively:

[0132]

[0133] where \(d\) i represents the deformation of the \(i\)-th point in the deformation field monitored by InSAR; represents the deformation of the \(i\)-th point by Kriging interpolation; \(N(D)\) represents the number of InSAR points in the monitoring area \(D\); RMSE k represents the root mean square error of the \(k\)-th point set in the set \(S\) M .

[0134] S34. The layout of GNSS stations involves the coordinate positions of each station. Therefore, the objective function should be a function of the position of each station. Taking the RMSE calculated in step S33 as an index, the problem of solving the position of the GNSS station layout can be transformed into the following problem of minimizing RMSE:

[0135]

[0136] where \(W\) is the weight or scale factor, and \(R\) represents the root mean square error of different GNSS station layouts.

[0137] More generally, when the two-dimensional or three-dimensional deformation of the study area can be obtained conditionally, the deformation of each dimension should be considered simultaneously. At this time, the above formula is transformed into the following form:

[0138]

[0139] where Ndim = 1, 2, 3, respectively representing the one-dimensional, two-dimensional, and three-dimensional deformations of the research area obtained; R i represents the root mean square error in the i-th direction; W i represents the weight in the i-th direction.

[0140] The weight can adjust the proportion of deformations in different dimensions. For example, when focusing on settlement, the weight of settlement can be increased, and the GNSS station layout results will more reflect the spatial distribution of settlement. At the same time, to reflect the binding force of the weight, it should meet the following requirements:

[0141]

[0142] S4. Determine the initial station distribution: The initial station distribution includes two types of stations. One is the stations that must be laid out according to the needs of the project and the project, and these stations will not change subsequently; the other is the stations that need to be iteratively optimized using the objective function, and these stations should be roughly close to the final distribution to improve the efficiency and accuracy of subsequent iterative solutions.

[0143] In some specific embodiments, step S4 may further include:

[0144] S41. Considering the trade-off between the model solution efficiency and accuracy in step S34, first determine an initial point distribution that is as close as possible to the optimal point distribution, and then perform iterative solution on the initial points. Therefore, this step is to determine the initial point distribution.

[0145] S411. First, make the following definitions:

[0146] Delete one point in turn from n candidate points, and the remaining all points form a point set Then

[0147]

[0148] where i = 1, 2,..., n, indicating that there are n such point sets; n - 1 indicates that there are n - 1 remaining points in the point set; k ≠ i indicates that the i-th point among the original n candidate points is deleted. Then use these n Q point sets to form a new set F n , then

[0149]

[0150] where n represents the number of elements in set F n There are elements in it.

[0151] S412. First, perform Kriging interpolation on these n candidate points using the variogram obtained in step S31 to obtain the interpolated deformation field. Then, calculate its root mean square error (RMSE) with the InSAR deformation monitoring according to the RMSE calculation formula in S33, denoted as r′0.

[0152] S413. For each point set in the set F n , perform Kriging interpolation using the variogram obtained in step S31 to obtain the interpolated deformation field. Then, calculate its RMSE with the InSAR deformation monitoring respectively according to the RMSE calculation formula in S33, denoted as r′ n , where n represents that there are n RMSE results.

[0153] S414. Calculate

[0154] r′ = r′ n - r′0

[0155] Record the point set corresponding to the minimum value in r′ and the candidate point (x i′ , y i′ , v i′ ) that is deleted. This minimum value indicates that deleting this point has the least impact on the overall interpolation result among all candidate points. Therefore, this point can be removed, and the point set is used as the new candidate point set.

[0156] S415. Let n = n - 1, and repeat steps S412 - S414 until n = 3 (because Kriging interpolation requires at least 3 points).

[0157] S416. Finally, obtain the order of deleting each candidate point from the initial n candidate points. When the number of GNSS stations N in the given study area is known, delete the first n - N points, and the point set Q′ N composed of the remaining N points is the initial station layout in this area.

[0158] S42. Determine the point set Q M composed of the station positions that must be set up due to engineering and project requirements. The positions of these stations will not change subsequently.

[0159] S43. To ensure that the number of stations is always N, it is necessary to find the point closest to each point in Q T in the initial point set Q′ M respectively, and replace it with the position of the point in Q M . Thus, the initial station Q N is obtained.

[0160] S5. Multi-scale iterative optimization to solve the objective function model: According to the InSAR deformation, the initial measurement stations, the candidate point set, and the variogram of the deformation field, a multi-scale iterative optimization method is used to solve the objective function model until the accuracy of the model solution meets the requirements.

[0161] In some specific embodiments, step S5 may further include:

[0162] S51. The initial measurement station Q determined through steps S2 to S4 N still needs to be further iteratively optimized to make it fit a higher deformation field as much as possible. To improve the computational efficiency in a large-scale scenario, we propose a multi-scale iterative optimization method for the final solution. First, downscale the initial candidate points to reduce the number of candidate points: increase the variance threshold in step S2 from 3.5 mm to 5 mm, and at this time, a new candidate point set D' is obtained n , and at this time, D' n has a density and number much smaller than D n ;

[0163] S52. Iteratively solve the initial measurement station Q N for a better solution X' n under the small-scale candidate point set D' N .

[0164] S521. For each point in the initial measurement station Q N except the fixed point Q M (assuming the current point is ), search for points in the small-scale candidate point set D' n with a distance less than the large distance threshold T1 = 200 m and form a point set C n1 , where n1 represents that there are n1 points in this point set.

[0165] S522. First, perform Kriging interpolation on these n candidate points in the initial measurement station Q N using the variogram obtained in step S31 to obtain its interpolated deformation field, and then calculate its root mean square error with the InSAR deformation monitoring according to the root mean square error calculation formula in S33, denoted as r″0.

[0166] S523. Replace the current point n1 with each point in C and perform Kriging interpolation using the variogram obtained in step S31 respectively to obtain its interpolated deformation field, and then calculate its root mean square error r″ n with the InSAR deformation monitoring according to the root mean square error calculation formula in S33 respectively.

[0167] S524. Calculate

[0168] r″ = r″ n -r″0

[0169] Record the minimum value r″ in r″ min and its corresponding Let j represent the j-th point in C corresponding to this minimum value. When r″ n1 is less than 0, it means min the point is better than the current point so use the point to replace the previous point Otherwise, the current point remains unchanged, and thus this iteration is completed.

[0170] S525. Let n = n + 1 (skip when encountering the fixed point Q M ), and repeat steps S523 - S524 until all points are iterated, and update the initial survey station to obtain the updated survey station X′'' N , and let Q N = X′ N .

[0171] S526. Repeat steps S521 - S525 until the distance between the corresponding positions of each point in the previous and current survey stations Q N is less than 10m. At this time, the final X′ N is the solution of the initial survey station Q N in the small-scale candidate point set D′ n with a better solution.

[0172] S53. Iteratively solve the final solution X N of X′ n under the large-scale candidate point set D N . Similar to the process of solving the survey station X′ N , but the initial survey station Q N needs to be replaced by X′ N , the small-scale candidate point set D′ n needs to be replaced by the large-scale candidate point set D n , and the larger distance threshold T1 = 200m needs to be replaced by the smaller distance threshold T2 = 100m. Then repeat steps S521 - S526 to solve the final GNSS survey station layout optimization result X N .

[0173] Through steps S1 to S5 above, we leverage the advantages of InSAR in deformation monitoring, including large scale, high precision, and high spatial resolution, to determine the placement of GNSS stations. For a given deformation field monitored by InSAR, we can quickly determine the reference number and locations of GNSS stations, maximizing the capture of deformation characteristics and reducing costs. This can also be used to improve the accuracy of subsequent InSAR and GNSS fusion.

[0174] Example 1

[0175] This example uses 29 scenes of Sentinel-1SAR data from January to December 2018 in a certain area.

[0176] Figure 2 This is the deformation rate map obtained by the SBAS-InSAR method (unit: mm). The blank areas around it are the masked-out areas that are irrelevant and unsuitable for station deployment. Figure 3 is the candidate point result obtained by this method; Figure 4a and Figure 4b The initial station distribution and deformation interpolation results obtained by this method are shown in Figure 2. Compared with the InSAR monitoring results, the RMSE is 3.33 mm. Figure 5a and Figure 5b This is the final station distribution and deformation interpolation result after iterative optimization of this method. Compared with the InSAR monitoring results, its RMSE is 2.97mm;

[0177] The results show that the GNSS station layout obtained by this method can indeed well reflect the spatial characteristics of deformation. Relative to the maximum settlement of -50 mm, the final RMSE of this method is only 2.97, indicating that the station layout method proposed by this method is practical and effective.

[0178] The above embodiments are merely preferred embodiments of the present invention, and the scope of protection of the present invention is not limited thereto. All technical solutions within the scope of protection of the present invention are within the scope of protection of the present invention. It should be noted that improvements and modifications that can be made by persons of ordinary skill in the art without departing from the principles of the present invention are also considered within the scope of protection of the present invention.

Claims

1. A GNSS station layout method based on InSAR deformation, characterized in that Including the following steps: S1. InSAR data processing: Obtain the differential interferogram set of radar images of a single or multiple orbits, and calculate the one-dimensional / two-dimensional / three-dimensional average surface deformation rate within a certain period of time through InSAR technology; S2. Determine the candidate point set: Use all InSAR deformation points as the initial point set for GNSS station layout. Through quadtree downsampling, select GNSS station candidate points that can reflect the deformation space characteristics in all deformation directions to reduce the quantity and density of the initial point set. At the same time, obtain the mask area where stations cannot be laid out based on external data such as the geological structure and topography of the research area, and remove the candidate points in the mask area to obtain the final candidate point set; S3. Establish the objective function model: First, estimate the variogram model of the InSAR deformation rate field. At the same time, with the minimum prediction variance as the goal, interpolate the entire deformation field through the Kriging method according to the selected variogram and measurement stations, and establish a model that minimizes the root mean square error between the interpolated deformation and the InSAR monitored deformation; S4. Determine the initial station distribution: The initial station distribution includes two types of stations. One is the stations that must be laid out according to the needs of the project and will not change subsequently. The other is the stations that need to be iteratively optimized using the objective function. The initial positions of these stations are determined by sequentially deleting candidate points according to the change in the root mean square; S5. Solve the objective function model through multi-scale iterative optimization: Solve the objective function model using the multi-scale iterative optimization method according to the InSAR deformation, initial stations, candidate point set, and variogram of the deformation field until the accuracy of the model solution meets the requirements.

2. The GNSS station layout method based on InSAR deformation according to claim 1, wherein The S1 includes the following sub-steps: S11. For the obtained radar images, after registration, selection of small baseline networking, interference, removal of flat ground and topography, phase filtering, phase unwrapping, atmospheric correction, and geocoding, form the interferogram set; S12. Assume there is only one spatio-temporal baseline set, then solve the time series deformation and average deformation rate by the least squares method through the short baseline set SBAS method; S13. According to the data situation of the obtained research area, further obtain one-dimensional, two-dimensional, or three-dimensional deformations respectively: One-dimensional deformation: If there is only one orbit, the line-of-sight LOS deformation monitored by SBAS-InSAR can be directly used and converted into vertical deformation: D U = D LOS / cosθ inc Among them, D U represents the vertical deformation, and D LOS is the LOS deformation monitored by SBAS-InSAR, and θ inc represents the local incident angle of the radar; Two-dimensional deformation: If ascending and descending orbit data can be obtained, obtain the two-dimensional deformation in the vertical and east-west directions by ascending and descending orbit fusion and ignoring the north-south deformation: Among them, D U and D E are the vertical and east-west deformations to be determined, respectively; represents the deformation monitored by ascending-track data, is the local radar incidence angle of the ascending-track satellite, is the radar azimuth angle of the ascending-track satellite; represents the deformation monitored by descending-track data, is the local radar incidence angle of the descending-track satellite, is the radar azimuth angle of the descending-track satellite; The parameters in the above formula are directly solved by the least squares method; Three-dimensional deformation: If ascending and descending orbits of different sensors can be obtained, use the classical method to obtain the deformations in multiple directions, and finally solve the deformations in the east-west, north-south, and vertical directions by the least squares method.

3. The GNSS station layout method based on InSAR deformation according to claim 1, wherein The S2 includes the following sub-steps: S21. Use all InSAR deformation points as the initial candidates, respectively obtain the candidate point sets in different directions through quadtree downsampling, and then take the union of the candidate point sets in different directions to obtain the candidate point set of GNSS stations; S211. Use the InSAR deformation point set as the initial candidate; S212. Set the variance threshold T of the deformation based on the InSAR deformation; S213. Use the quadtree recursive algorithm to continuously subdivide the InSAR point set into squares, and calculate the variance of the InSAR deformation in each square: Among them, σ 2 represents the variance of the InSAR points within the square, n represents the number of InSAR points within the square, and d i represents the deformation of the i-th point within the square, is the mean value of the n InSAR deformation points within the square; If the variance is less than or equal to T, stop dividing; if the variance is greater than T, continue dividing until the variance of the deformation in each square is less than or equal to T; Use the point set obtained by downsampling the quadtree as the candidate point set in this direction; S214. Take the union of the candidate point sets in multiple directions to obtain the initial candidate point set for the entire region; S22. Obtain the available geological structure and topographic and geomorphic external data of the study area, and mask the areas that are not suitable for setting up survey stations; S23. Remove the candidate points in the masked area of S24 from the initial candidate point set to obtain the final candidate point set.

4. The GNSS station layout method based on InSAR deformation according to claim 1, wherein The said S3 includes the following sub-steps: S31. Calculate the variation value of the two-dimensional deformation field using the InSAR deformation points: Where h is the spacing between deformation points, also known as the potential difference; N(h) is the number of pairs of all observation points with a spacing of h; z(x i ) and z(x i +h) respectively represent the InSAR deformation observation values with a relative distance of h; is the variance at a spacing of h, that is, the variation value; the variance of the InSAR deformation field and the sample value of the spacing h are obtained through this step; The goal is to establish the relationship between the variance and the spatial distance h between point pairs, that is, the variogram; according to the sample values obtained above, select the model with the best goodness of fit as the variogram and fit and calculate the model parameters; S32. Assume that a total of n candidate points are obtained in a certain area D in step S2, and the set D n is used to represent it. Now, it is necessary to select N points from these n candidate points to set up GNSS stations. Then, the following sets and symbols are defined: Use the set P N to represent the set of the GNSS point distributions; use (x i , y i , v i ) to represent the position and its deformation rate of the i-th point in the current point set, where i = 1, 2, … N, simplified as Denoted by S M represents a series of sets P N with the same number of GNSS stations but different distributions, where M represents that there are M sets P in this set N , then then where and both represent the k-th point set; S33. For the above M Ps N perform Kriging interpolation respectively to obtain the complete regional deformation field after interpolation; then compare it with the deformation field monitored by InSAR, and calculate the root mean square error RMSE of each GNSS distribution set respectively: Among them, d i represents the deformation of the i-th point in the deformation field monitored by InSAR; represents the deformation of the i-th point of Kriging interpolation; N(D) represents the number of InSAR points in the monitoring area D; RMSE k represents the set S M the root mean square error of the k-th point set in; S34. Taking the RMSE calculated in step S33 as an index, the problem of solving the layout position of the GNSS survey station is transformed into the following problem of minimizing RMSE: Where W is the weight or scaling factor, and R represents the root mean square error of different survey station layouts; When the two-dimensional or three-dimensional deformation of the study area can be obtained under certain conditions, taking into account the deformation in each dimension, the above formula is transformed into the following form: where Ndim = 1, 2, 3, representing the one-dimensional, two-dimensional, and three-dimensional deformations of the research area obtained respectively; R i represents the root mean square error in the i-th direction; W i represents the weight in the i-th direction; The weight can adjust the proportion of deformations in different dimensions. For example, when focusing on subsidence, increase the weight of subsidence, then the layout result of the GNSS survey station will more reflect the spatial distribution of subsidence; At the same time, in order to reflect the binding force of the weight, it should meet the following requirements:

5. The GNSS station layout method based on InSAR deformation according to claim 4, wherein The said S4 includes the following sub-steps: S41. Considering the trade-off between the model solving efficiency and accuracy in step S34, first determine an initial point distribution as close as possible to the optimal point distribution, and then iteratively solve the initial points; S411. First, make the following definitions: Delete one point from n candidate points in sequence, and all the remaining points form a point set Then: where \(i = 1, 2, \ldots, n\), indicating that there are \(n\) such point sets; \(n - 1\) represents the remaining \(n - 1\) points in the point set; \(k\neq i\) means that the \(i\)-th point among the original \(n\) candidate points is deleted; then use these \(n\) \(Q\) point sets to form a new set \(F\) n , then: where n represents the number of elements in set F n ; S412. First, perform Kriging interpolation on these n candidate points using the variogram obtained in step S31 to obtain the interpolated deformation field, and then calculate its root mean square error with the InSAR deformation monitoring according to the formula for calculating the root mean square error in S33, denoted as r′0; S413. Perform Kriging interpolation on each point set in the set F n using the variogram obtained in step S31 to obtain its interpolated deformation field, and then calculate the root mean square error between it and the InSAR deformation monitoring according to the root mean square error calculation formula in S33, denoted as r′ n , where n represents that there are n root mean square error results; S414. Calculate r′＝r′ n -r′0 The point set corresponding to the minimum value in the record r' and the candidate point (x i′ , y i′ , v i′ ) to be deleted; this minimum value indicates that among all candidate points, deleting this point has the least impact on the overall interpolation result. Remove this point and use the point set as the new candidate point set; S415. Let n = n - 1, and repeat steps S412 - S414 until n = 3; S416. Finally, obtain the order of deleting each candidate point from the initial n candidate points; when the number of GNSS stations N in the given research area is known, delete the first n - N points, and the remaining N points form a point set Q'. N This is the initial station layout for this area. S42. Determine the point set Q composed of the survey station positions that must be arranged due to engineering and project requirements. M The positions of these survey stations will not change subsequently. S43. To ensure that the number of measurement stations is always N, it is necessary to find the points closest to each point in Q in the initial point set Q' T respectively, and replace them with the positions in Q M . Thus, the initial measurement station Q M is obtained N .

6. The GNSS station layout method based on InSAR deformation according to claim 4, wherein The said S5 includes the following sub-steps: S51. First, downscale the initial candidate points to reduce the number of candidate points: increase the variance threshold in step S2 to a certain value, and at this time, a new candidate point set D' is obtained n , and at this time D' n has a density and number much smaller than D n ; S52. Iteratively solve for the initial survey station Q N in the small-scale candidate point set D′ n for a better solution X′ N ; S521. For the initial measurement station Q N except for the fixed point Q M in each point, assume the current point is search for points with a distance less than the distance threshold T1 in the small-scale candidate point set D′ n and form a point set C n1 , where n1 represents that there are n1 points in this point set; S522. First, for these n candidate points in the initial measurement station Q N perform Kriging interpolation using the variogram obtained in step S31 to obtain the interpolated deformation field, and then calculate its root mean square error with the InSAR deformation monitoring according to the calculation formula of the root mean square error in S33, denoted as r″0; S523, use C in sequence n1 Each point in replaces the current point The variogram obtained in step S31 is used for Kriging interpolation to obtain the interpolated deformation field, and then the root mean square error r" between the interpolated deformation field and the InSAR deformation monitoring is calculated according to the root mean square error calculation formula in S33. n ; S524 Calculate r″＝r″ n -r″0 Record the minimum value r″ in r″ min and its corresponding j indicates that the minimum value corresponds to C n1 The jth point in min <0, description Point than the current point Better, so use Replace the front point Otherwise the current point Keep still and complete this iteration; S525. Let n = n + 1. When the fixed point Q is encountered M skip it and repeat steps S523 - S524 until all points are iterated and the initial station is updated to obtain the updated station X'. N and let Q N = X'. N ; S526. Repeat steps S521 - S525 until the distance between the corresponding positions of each point in the previous and current survey stations Q N is less than a given threshold; at this time, the final X' N is the solution of the initial survey station Q N which is a better solution under the small-scale candidate point set D' n ; S53. Iteratively solve for X' N For the final solution X n under the large-scale candidate point set D N , the process is the same as solving for the survey station X' N , but the initial survey station Q N needs to be replaced by X' N , the small-scale candidate point set D' n needs to be replaced by the large-scale candidate point set D n , the distance threshold T1 needs to be replaced by a new distance threshold T2, and the new distance threshold T2 is less than the distance threshold T1; then repeat steps S521 - S526 to solve for the final optimized GNSS survey station layout result X N .

7. Application of the GNSS survey station layout method based on InSAR deformation according to any one of claims 1 - 6 in synthetic aperture interferometric radar measurement and global navigation satellite system.

Images