CN116052009A - GNSS station measurement layout method based on InSAR deformation and application thereof - Google Patents

GNSS station measurement layout method based on InSAR deformation and application thereof Download PDF

Info

Publication number
CN116052009A
CN116052009A CN202211356074.XA CN202211356074A CN116052009A CN 116052009 A CN116052009 A CN 116052009A CN 202211356074 A CN202211356074 A CN 202211356074A CN 116052009 A CN116052009 A CN 116052009A
Authority
CN
China
Prior art keywords
deformation
insar
station
point
points
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202211356074.XA
Other languages
Chinese (zh)
Inventor
毛文祥
李志伟
许兵
朱焱
侯景鑫
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202211356074.XA priority Critical patent/CN116052009A/en
Publication of CN116052009A publication Critical patent/CN116052009A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/13—Satellite images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/94—Hardware or software architectures specially adapted for image or video understanding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/188—Vegetation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Multimedia (AREA)
  • Theoretical Computer Science (AREA)
  • Astronomy & Astrophysics (AREA)
  • Remote Sensing (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Software Systems (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention provides a GNSS station layout method based on InSAR deformation and an application thereof, which take InSAR historical deformation as a basis, give out the reference number, the spatial distribution and the corresponding deformation prediction precision of the GNSS station layout of a deformation area through multi-scale iteration on the premise of considering cost and precision by using a Kriging interpolation method in space sampling and geostatistics, so that the GNSS station can capture the spatial characteristics of deformation to the greatest extent, and the precision of post-InSAR and GNSS fusion can be improved. The invention fully utilizes the advantages of InSAR in large scale, high precision and high spatial resolution on deformation monitoring, provides guidance and reliable basis for station arrangement during GNSS deformation monitoring, and is beneficial to better capturing the spatial characteristics of regional deformation. Under the background of frequent geological disasters at present, the monitoring and recognition of the geological disasters are very facilitated.

Description

GNSS station measurement layout method based on InSAR deformation and application thereof
Technical Field
The invention relates to a synthetic aperture radar interferometry (InSAR, synthetic Aperture Radar Interferometry) technology and a global navigation satellite system (GNSS, global Navigation Satellite System) technology in the field of geodetic measurement based on remote sensing images, in particular to a GNSS station layout method based on InSAR deformation and application thereof.
Background
Global navigation satellite system (GNSS, global Navigation Satellite System) is one of the most mature precision positioning technologies at present, which has the advantages of all-weather, automation and long-time continuous observation, and is therefore widely used in many engineering measurement and natural disaster monitoring fields. However, GNSS has a high cost, and cannot be deployed in a large area and in a high density, so that the spatial resolution of the monitoring result is low, and the spatial change of the deformation of the monitoring area cannot be reflected. The technology of synthetic aperture radar interferometry (InSAR, interferometric Synthetic Aperture Radar) has been developed for decades, and has made great breakthroughs and achievements in various fields such as volcanic, earthquake, mining areas, glaciers, atmosphere, frozen soil and the like rapidly with the unique advantages, so that the technology is one of the effective means for monitoring geological disasters at present. And it has the advantages of high spatial resolution, high accuracy and large-scale rapid monitoring. However, the time resolution of the InSAR monitoring result is limited by the revisit period of the sensor, and even if the A constellation and the B constellation of the Sentinel-1 satellite are used for joint observation, the revisit period is six days, so that the deformation characteristics of the sensor still cannot be completely and correctly reflected when the deformation of the sensor, such as mining areas, landslide, groundwater and the like, is strong in burstiness, obvious in nonlinear characteristics and large in deformation gradient.
There have been many studies on fusion of InSAR and GNSS technologies, such as correcting the atmospheric delay in InSAR using GNSS, combining GNSS and InSAR to solve three-dimensional deformation, and fusing InSAR and GNSS deformation results to obtain a surface deformation result with high space-time resolution, however, these studies have hardly considered how to arrange GNSS points to make the fusion of InSAR and GNSS best, and how to balance the relationship between cost and monitoring accuracy. In particular, the deformation information is used as the most intuitive precursor reflection of geological disasters, the conventional station layout for GNSS deformation monitoring mainly depends on experience guidance and lacks theoretical basis, so that the obtained GNSS layout scheme is only one feasible scheme, a certain cost waste and low deformation monitoring precision are caused, and the deformation monitoring result is difficult to reflect the spatial distribution characteristics of deformation of the whole research area. Therefore, how to consider the space distribution and the rationality of the number of the layout of the GNSS, so that the fusion accuracy of the InSAR and the GNSS deformation is the highest and the cost is the application problem to be solved.
Disclosure of Invention
The invention aims to overcome the defect that the current GNSS deformation monitoring station layout lacks guidance and reliable basis, and provides an InSAR deformation-based GNSS station layout method and application thereof.
The invention firstly provides the following technical scheme:
a GNSS station layout method based on InSAR deformation comprises the following steps:
s1, inSAR data processing: and obtaining a differential interference atlas of radar images of a single track or a plurality of tracks, and calculating the one-dimensional/two-dimensional/three-dimensional average deformation rate of the earth surface in a period of time by using an InSAR technology.
S2, determining a candidate point set: and taking all InSAR deformation points as an initial point set for GNSS station measurement layout, and selecting GNSS station measurement candidate points capable of showing deformation space characteristics of all deformation directions through quadtree downsampling so as to reduce the number and density of the initial point set. And meanwhile, obtaining a mask region in which a measuring station cannot be arranged according to external data such as geological structures, topography and the like of the research region, removing candidate points in the mask region, and obtaining a final candidate point set.
S3, establishing an objective function model: firstly, estimating a variation function model of an InSAR deformation rate field, simultaneously taking a minimum prediction variance as a target, interpolating the whole deformation field through a kriging method according to a selected variation function and a measuring station, and establishing a model capable of minimizing Root Mean Square Error (RMSE) of interpolation deformation and InSAR monitoring deformation.
S4, determining initial station measurement distribution: the initial station distribution comprises two types of stations, one type of station is a station which is determined to be arranged according to engineering and project requirements, and the station is not changed later; the other is a station requiring iterative optimization using an objective function, and the initial position of the station is determined according to the root mean square variation by sequentially deleting candidate points.
S5, solving an objective function model through multi-scale iterative optimization: and solving an objective function model by adopting a multi-scale iterative optimization method according to InSAR deformation, an initial station, a candidate point set and a variation function of a deformation field until the precision of model solving reaches the requirement.
According to some preferred embodiments of the invention, the S1 comprises:
s11, registering, selecting a small baseline network, interfering, removing land and topography, phase filtering, phase unwrapping, atmosphere correction and geocoding are carried out on the obtained radar image to form an interference atlas.
S12, assuming that only one space-time baseline set exists, the time sequence deformation and average deformation rate can be obtained through a short baseline set (Small BAseline Subsets, SBAS) method through least square.
S13, according to the acquired data conditions of the research area, one-dimensional, two-dimensional or three-dimensional deformation can be further acquired respectively:
one-dimensional deformation: if there is only one track, the Line of Sight (LOS) deformation of SBAS-InSAR monitoring can be directly used, or can be converted into a vertical deformation:
D U ＝D LOS /cosθ inc
wherein D is U Represents the deformation in the vertical direction, D LOS LOS directional deformation, θ for SBAS-InSAR monitoring inc Representing the radar local angle of incidence.
Two-dimensional deformation: if the lifting rail data can be obtained, the two-dimensional deformation in the vertical direction and the east-west direction can be obtained by fusing the lifting rails and neglecting the deformation in the north-south direction:
Figure BDA0003919932250000031
Figure BDA0003919932250000032
wherein D is U And D E Respectively the vertical deformation and the east-west deformation to be solved;
Figure BDA0003919932250000033
representing the deformation of the monitoring of the lift rail data,
Figure BDA0003919932250000034
radar local angle of incidence for an orbiting satellite, +.>
Figure BDA0003919932250000035
Radar azimuth for an orbiting satellite; />
Figure BDA0003919932250000036
Representing the deformation of the derailment data monitoring, +.>
Figure BDA0003919932250000037
Radar local angle of incidence for a down-track satellite, +.>
Figure BDA0003919932250000038
Is the radar azimuth of the down-track satellite.
The parameters in the above equation can be solved directly by least squares.
Three-dimensional deformation: if the lifting rail of different sensors can be obtained, the deformation of the lifting rail in multiple directions is obtained by classical methods such as D-InSAR, offset Tracking (Offset-Tracking), multi-Aperture InSAR (Multi-Aperture InSAR) and the like, and finally the deformation in the east-west, south-north and vertical directions is solved by least square.
According to some preferred embodiments of the invention, the S2 comprises:
s21, taking all InSAR deformation points as initial candidate points, respectively obtaining candidate point sets in different directions through quadtree downsampling, and then obtaining a union set of the candidate point sets in different directions to obtain a candidate point set of the GNSS station.
S211, taking the InSAR deformation point set as an initial candidate point (the InSAR deformation point set can be in a grid data format or a vector data format);
s212, considering that the variance can reflect the spatial variation characteristic of deformation, setting a variance threshold T of deformation based on InSAR deformation;
s213, continuously subdividing the InSAR point set into squares by using a quadtree recursion algorithm, and calculating the variance of InSAR deformation in each square:
Figure BDA0003919932250000039
wherein sigma 2 Representing the variance of the InSAR points in the square, n represents the number of InSAR points in the square, and d i Representing the deformation of the i-th point in the square,
Figure BDA00039199322500000310
is the mean of n InSAR deformation points in the square.
If the variance is less than or equal to T, stopping segmentation; if the variance is greater than T, continuing to divide until the variance of the deformation in each square is less than or equal to T. The set of points downsampled by the quadtree is taken as the candidate set of points for that direction.
S214, obtaining a union set of candidate point sets in multiple directions, and obtaining an initial candidate point set of the whole area.
S22, obtaining external data such as available geological structures, topography and the like of a research area, and masking off areas unsuitable for laying measuring stations, such as areas with overlarge DEM gradient and overfull vegetation;
S23, removing the candidate points in the mask area from the initial candidate point set to obtain a final candidate point set.
According to some preferred embodiments of the invention, the S3 comprises:
s31, calculating variation values of the two-dimensional deformation field by using all InSAR deformation points (when the data volume is large, the data volume can be reduced by adopting random sampling):
Figure BDA0003919932250000041
wherein h is the distance between deformation points, also called the potential difference; n (h) is the paired number of all observation points with h as a distance; z (x) i ) And z (x) i +h) represent InSAR deformation observations at a relative distance h, respectively;
Figure BDA0003919932250000042
the variance, i.e., the variance value, of the spacing h increases with increasing h over a range. When the actual point distance is greater than the maximum correlation distance, the value tends to stabilize. By this step the variance of the InSAR deformation field can be obtained +.>
Figure BDA0003919932250000043
Sample value from the distance h.
The aim is to establish the variance
Figure BDA0003919932250000044
And the spatial distance h between pairs of points, i.e. the variation function. Since the variation function has many classical models in the stabilization process, such as an exponential model, a gaussian model, a spherical model, etc. And selecting a model with the best fitting degree as a variation function according to the obtained sample value, and fitting and calculating the model parameters.
S32, assuming that n candidate points are obtained in total in a certain region D in the step S2, using a set D n Indicating that it is now necessary to choose N points from the N candidate points to lay GNSS stations, the following set and symbol are defined:
with set P N Representing a set of these GNSS point distributions; by (x) i ,y i ,v i ) Representing the ith in the current point setThe position of the point and its deformation rate, where i=1, 2, … N, can be simplified to
Figure BDA0003919932250000045
By S M Representing a series P with the same number of GNSS stations but different distribution N Wherein M represents a total of M sets P in the set N Then
Figure BDA0003919932250000046
Then->
Figure BDA0003919932250000047
Wherein->
Figure BDA0003919932250000048
And->
Figure BDA0003919932250000049
All represent the kth point set.
S33, taking the Kriging interpolation into consideration, performing unbiased optimal estimation on the regional variable, thereby performing M P' S on the regional variable N And (5) respectively performing Kerling interpolation to obtain an interpolated complete area deformation field. And then comparing the deformation field monitored by the InSAR with the deformation field monitored by the InSAR, and respectively calculating Root Mean Square Errors (RMSE) of each GNSS distribution set:
Figure BDA00039199322500000410
wherein d i Representing deformation of an ith point in the InSAR monitored deformation field;
Figure BDA00039199322500000411
representing the deformation of the ith point of the kriging interpolation; n (D) represents the number of InSAR points in the monitoring area D; RMSE k Representation set S M Root mean square error for the kth point set.
The layout of the GNSS stations involves the coordinate position of each station, and thus the objective function should be a function of each station position. Taking the RMSE calculated in step S33 as an index, the position solution problem of the GNSS station arrangement may be converted into the following RMSE minimization solution problem:
Figure BDA0003919932250000051
Wherein W is a weight or a scale factor, and R represents root mean square error of different station arrangements.
More generally, when a two-dimensional or three-dimensional deformation of the investigation region is conditionally acquired, the deformation of each dimension should be taken into account at the same time, and the above equation is converted into the following form:
Figure BDA0003919932250000052
where ndim=1, 2,3, respectively representing the one-, two-and three-dimensional deformations of the investigation region obtained; r is R i Root mean square error representing the i-th direction; w (W) i Indicating the weight of the i-th direction.
The weight can adjust the duty ratio of deformation in different dimensions, for example, when the weight is focused on sedimentation, the sedimentation weight can be increased, and the GNSS station layout result can reflect more sedimentation space distribution. Meanwhile, in order to embody the constraint force of the weight, the following requirements should be satisfied:
Figure BDA0003919932250000053
according to some preferred embodiments of the invention, the S4 comprises:
s41, considering the problem of model solving efficiency and accuracy compromise in the step S34, firstly determining an initial point location distribution which is as close to the optimal point location distribution as possible, and then carrying out iterative solving on the initial point location. The step therefore consists in determining an initial distribution of points.
S411, first, the following definition is performed:
from n candidate pointsSequentially deleting one point, and forming a point set by all the rest points
Figure BDA0003919932250000054
Then
Figure BDA0003919932250000055
Where i=1, 2, …, n, represents n such point sets; n-1 represents the remaining n-1 points in the point set; k+.i indicates that the ith point of the original n candidate points is deleted. Then the n Q point sets are used to form a new set F n Then
Figure BDA0003919932250000056
Wherein n represents the set F n With one element.
S412, performing Kerling interpolation on the n candidate points by adopting the variation function obtained in the step S31 to obtain an interpolated deformation field, and calculating the root mean square error of the interpolation field and InSAR deformation monitoring according to the calculation formula of the root mean square error in the step S33, wherein the root mean square error is denoted as r '' 0 。
S413, pair set F n Each point set in (2) is subjected to Kerling interpolation by adopting the variation function obtained in the step S31 to obtain an interpolation deformation field, and then the root mean square error of the deformation field and InSAR deformation monitoring is calculated according to a calculation formula of the root mean square error in the step S33 and is recorded as r '' n N represents the root mean square error result, where n is the number.
S414, calculating
r′＝r′ n -r′ 0
Recording the point set corresponding to the minimum value in r
Figure BDA0003919932250000061
And deleted candidate points (x i′ ,y i′ ,v i′ ). The minimum value represents the influence of the point deleted from all candidate points on the overall interpolation resultMinimum, the point can therefore be removed and the set of points
Figure BDA0003919932250000062
As a new candidate point set.
S415, let n=n-1, and repeat steps S412 to S414 until n=3 (since at least 3 points are needed for kriging).
S416, finally obtaining the order of deleting each candidate point from the initial n candidate points. When the number N of GNSS stations in a research area is given, deleting the first N-N points and forming a point set Q 'consisting of the rest N points' N I.e. the initial start station for the area.
S42, determining a point set Q consisting of station positions which are required to be laid out according to engineering and project requirements M This part of the station position will not subsequently change.
S43, in order to ensure that the number of measuring stations is always N, a set Q 'is required to be collected at the initial point' N Find the distance Q respectively M Closest to each point in (a) and using Q M The point positions in the test station are replaced, so as to obtain an initial measuring station Q N 。
According to some preferred embodiments of the invention, the S5 comprises:
s51, determining an initial measuring station Q through the steps S2-S4 N Further iterative optimization is required to fit as much higher deformation fields as possible. In order to improve the calculation efficiency in a large-scale scene, we propose a multi-scale iterative optimization method for final solution. Firstly, the initial candidate points are downscaled to reduce the number of the candidate points: increasing the variance threshold in step S2 to a certain value, and obtaining a new candidate point set D' n At this time D' n The density and number of (C) are far less than D n ；
S52, iteratively solving initial measuring station Q N At a small scale candidate point set D' n The following more optimal solution X' N 。
S521 for initial measuring station Q N Middle-divided fixed point Q M Each point outside (assuming the current point is
Figure BDA0003919932250000063
) At a larger distance threshold T 1 (e.g., 2-fold change of the variation function), at a small scale candidate point set D' n The search distance is less than T 1 And form a point set C n1 Where n1 represents that there are n1 points in the set of points.
S522, first, for initial station Q N Performing Kerling interpolation on the n candidate points in the step S31 by adopting a variation function obtained in the step S31 to obtain an interpolated deformation field, and then calculating the root mean square error of the deformation field and InSAR deformation monitoring according to a calculation formula of the root mean square error in the step S33, and marking the root mean square error as r' 0 。
S523, sequentially using C n1 Each of which replaces the current point
Figure BDA0003919932250000064
And performing Kriging interpolation by using the variation function obtained in the step S31 to obtain an interpolated deformation field, and then respectively calculating the root mean square error r' of the deformation field and InSAR deformation monitoring according to the calculation formula of the root mean square error in the step S33 n 。/>
S524, calculate
r″＝r″ n -r″ 0
Record the minimum value r 'in r' min And corresponding to
Figure BDA0003919932250000065
j represents the minimum value corresponding to C n1 Is the j-th point in (a). When r' min <0, description->
Figure BDA0003919932250000071
Point is +. >
Figure BDA0003919932250000072
Better, thus use the point +.>
Figure BDA0003919932250000073
Replace the anterior point->
Figure BDA0003919932250000074
Otherwise, the current point->
Figure BDA0003919932250000075
And keeping the iteration still, thus completing the iteration.
S525, let n=n+1 (encountering fixed point Q M Skip) and repeat steps S523-S524 until all points are iterated, and update the initial station to obtain updated station X '' N And let Q N ＝X′ N 。
S526, repeating the steps S521-S525 until the station Q is detected twice N The distance between the corresponding positions of each of the points in (a) is less than a given threshold. At this time last X' N I.e. solving the initial measuring station Q N At a small scale candidate point set D' n The following is a better solution.
S53, iteratively solving X' N At a large scale candidate point set D n Final solution X below N . With solving station X' N But requires the initial station Q N By X' N Instead, a small scale candidate point set D' n With a large set of candidate points D n Instead, a larger distance threshold T 1 With a smaller distance threshold T 2 (e.g., half the variation of the variation function). And repeating the steps S521-S526 to obtain the final GNSS station layout optimization result X N 。
According to the method for determining the GNSS station by multi-scale iterative optimization based on InSAR deformation, a system for determining the arrangement of the GNSS station based on InSAR deformation can be obtained.
The GNSS station measurement layout method or system can be used for the station measurement layout of GNSS deformation monitoring.
The GNSS station layout method based on InSAR deformation fully utilizes the advantages of large scale, high precision and high spatial resolution of InSAR in deformation monitoring. The method can be used for rapidly giving out the reference number and the reference positions of the GNSS station arrangement aiming at a certain deformation field of InSAR monitoring, capturing deformation characteristics to the greatest extent, saving cost and improving the precision of the late InSAR and GNSS fusion. Under the background of frequent geological disasters at present, the monitoring and recognition of the geological disasters are very facilitated.
Drawings
FIG. 1 is a schematic flow chart showing one embodiment of the method of the present invention;
FIG. 2 is a graph of deformation rate obtained by SBAS-InSAR in example 1;
FIG. 3 is a plot of candidate station profiles of example 1;
FIG. 4a is an initial station distribution of example 1;
FIG. 4b is a graph of deformation results predicted using an initial station in example 1;
FIG. 5a is a final station diagram of example 1;
FIG. 5b is a graph of the deformation results predicted using the final station in example 1.
Detailed Description
The present invention will be described in detail with reference to the following examples and drawings, but it should be understood that the examples and drawings are only for illustrative purposes and are not intended to limit the scope of the present invention in any way. All reasonable variations and combinations that are included within the scope of the inventive concept fall within the scope of the present invention.
A specific GNSS station layout method based on InSAR deformation as shown in fig. 1 includes the following steps:
s1, inSAR data processing: and obtaining a differential interference atlas of radar images of a single track or a plurality of tracks, and calculating the one-dimensional/two-dimensional/three-dimensional average deformation rate of the earth surface in a period of time by using an InSAR technology.
In some embodiments, step S1 may further include:
s11, roughly registering the acquired radar image by adopting track information, and then finely registering by adopting intensity cross-correlation; setting a small baseline threshold value, selecting interference pairs meeting a space-time baseline threshold value (the time baseline is 200 days and the space baseline is 800 meters), and carrying out interference, land leveling, terrain removal and phase filtering; the unwrapping reference point is selected for phase unwrapping, followed by atmospheric correction and geocoding.
S12, assuming that only one space-time baseline set exists, establishing a deformation rate and terrain residual error model, and obtaining the surface deformation through least square calculation.
S13, according to the acquired data conditions of the research area, one-dimensional, two-dimensional or three-dimensional deformation can be further acquired respectively:
one-dimensional deformation: if there is only one track, the Line of Sight (LOS) deformation of SBAS-InSAR monitoring can be directly used, or can be converted into a vertical deformation:
D U ＝D LOS /cosθ inc
Wherein D is U Represents the deformation in the vertical direction, D LOS LOS directional deformation, θ for SBAS-InSAR monitoring inc Representing the radar local angle of incidence.
Two-dimensional deformation: if the lifting rail data can be obtained, the two-dimensional deformation in the vertical direction and the east-west direction can be obtained by fusing the lifting rails and neglecting the deformation in the north-south direction:
Figure BDA0003919932250000081
Figure BDA0003919932250000082
wherein D is U And D E Respectively the vertical deformation and the east-west deformation to be solved;
Figure BDA0003919932250000083
representing the deformation of the monitoring of the lift rail data,
Figure BDA0003919932250000084
radar local angle of incidence for an orbiting satellite, +.>
Figure BDA0003919932250000085
Radar azimuth for an orbiting satellite; />
Figure BDA0003919932250000086
Representing the deformation of the derailment data monitoring, +.>
Figure BDA0003919932250000087
Radar local angle of incidence for a down-track satellite, +.>
Figure BDA0003919932250000088
Is the radar azimuth of the down-track satellite.
The parameters in the above equation can be solved directly by least squares.
Three-dimensional deformation: if the lifting rail of different sensors can be obtained, the deformation of the lifting rail in multiple directions is obtained by classical methods such as D-InSAR, offset Tracking (Offset-Tracking), multi-Aperture InSAR (Multi-Aperture InSAR) and the like, and finally the deformation in the east-west, south-north and vertical directions is solved by least square.
S2, determining a candidate point set: and taking all InSAR deformation points as an initial point set for GNSS station measurement layout, and selecting GNSS station measurement candidate points capable of showing deformation space characteristics of all deformation directions through quadtree downsampling so as to reduce the number and density of the initial point set. And meanwhile, obtaining a mask region in which a measuring station cannot be arranged according to external data such as geological structures, topography and the like of the research region, removing candidate points in the mask region, and obtaining a final candidate point set.
In some embodiments, step S2 may further include:
s21, taking all InSAR deformation points as initial candidate points, respectively obtaining candidate point sets in different directions through quadtree downsampling, and then obtaining a union set of the candidate point sets in different directions to obtain a candidate point set of the GNSS station.
S211, taking the InSAR deformation point set as an initial candidate point (the InSAR deformation point set can be in a grid data format or a vector data format);
s212, considering that the variance can reflect the spatial variation characteristic of deformation, setting a variance threshold T=3.5mm of deformation based on InSAR deformation;
s213, continuously subdividing the InSAR point set into squares by using a quadtree recursion algorithm, and calculating the variance of InSAR deformation in each square:
Figure BDA0003919932250000091
wherein sigma 2 Representing the variance of the InSAR points in the square, n represents the number of InSAR points in the square, and d i Representing the deformation of the i-th point in the square,
Figure BDA0003919932250000092
is the mean of n InSAR deformation points in the square.
If the variance is less than or equal to T, stopping segmentation; if the variance is greater than T, continuing to divide until the variance of the deformation in each square is less than or equal to T. The set of points downsampled by the quadtree is taken as the candidate set of points for that direction. .
S214, obtaining a union set of candidate point sets in multiple directions to obtain a candidate point set of the whole area.
S22, obtaining external data such as available geological structures, topography and the like of a research area, and masking off areas unsuitable for laying measuring stations, such as areas with overlarge DEM gradient and overfull vegetation;
s23, removing the candidate points in the mask area from the initial candidate point set to obtain a final candidate point set.
S3, establishing an objective function model: firstly, estimating a variation function model of an InSAR deformation rate field, simultaneously taking a minimum prediction variance as a target, interpolating the whole deformation field through a kriging method according to a selected variation function and a measuring station, and establishing a model capable of minimizing Root Mean Square Error (RMSE) of interpolation deformation and InSAR monitoring deformation.
In some embodiments, step S3 may further include:
s31, calculating variation values of the two-dimensional deformation field by using all InSAR deformation points (when the data volume is large, the data volume can be reduced by adopting random sampling):
Figure BDA0003919932250000093
wherein h is the distance between deformation points, also called the potential difference; n (h) is the paired number of all observation points with h as a distance; z (x) i ) And z (x) i +h) represent InSAR deformation observations at a relative distance h, respectively;
Figure BDA0003919932250000094
The variance, i.e., the variance value, of the spacing h increases with increasing h over a range. When the actual point distance is greater than the maximum correlation distance, the value tends to stabilize. By this step the variance of the InSAR deformation field can be obtained +.>
Figure BDA0003919932250000095
Sample value from the distance h.
The aim is to establish the variance
Figure BDA0003919932250000101
And the spatial distance h between pairs of points, i.e. the variation function. Since the variation function has many classical models in the stabilization process, such as an exponential model, a gaussian model, a spherical model, etc. And selecting a model with the best fitting degree as a variation function according to the obtained sample value, and fitting and calculating the model parameters.
S32, assuming that n candidate points are obtained in total in a certain region D in the step S2, using a set D n Indicating that it is now necessary to choose N points from the N candidate points to lay GNSS stations, the following set and symbol are defined:
with set P N Representing a set of these GNSS point distributions; by (x) i ,y i ,v i ) Representing the position of the i-th point in the current point set and its deformation rate, where i=1, 2, … N, can be simplified to
Figure BDA0003919932250000102
By S M Representing a series P with the same number of GNSS stations but different distribution N Wherein M represents a total of M sets P in the set N Then
Figure BDA0003919932250000103
Then->
Figure BDA0003919932250000104
Wherein->
Figure BDA0003919932250000105
And->
Figure BDA0003919932250000106
All represent the kth point set.
S33, taking the Kriging interpolation into consideration, performing unbiased optimal estimation on the regional variable, thereby performing M P' S on the regional variable N And (5) respectively performing Kerling interpolation to obtain an interpolated complete area deformation field. And then comparing the deformation field monitored by the InSAR with the deformation field monitored by the InSAR, and respectively calculating Root Mean Square Errors (RMSE) of each GNSS distribution set:
Figure BDA0003919932250000107
wherein d i Representing deformation of an ith point in the InSAR monitored deformation field;
Figure BDA0003919932250000108
representing the deformation of the ith point of the kriging interpolation; n (D) represents the number of InSAR points in the monitoring area D; RMSE k Representation set S M Root mean square error for the kth point set.
The layout of the GNSS stations involves the coordinate position of each station, and thus the objective function should be a function of each station position. Taking the RMSE calculated in step S33 as an index, the position solution problem of the GNSS station arrangement may be converted into the following RMSE minimization solution problem:
Figure BDA0003919932250000109
wherein W is a weight or a scale factor, and R represents root mean square error of different station arrangements.
More generally, when a two-dimensional or three-dimensional deformation of the investigation region is conditionally acquired, the deformation of each dimension should be taken into account at the same time, and the above equation is converted into the following form:
Figure BDA00039199322500001010
Where ndim=1, 2,3, respectively representing the one-, two-and three-dimensional deformations of the investigation region obtained; r is R i Root mean square error representing the i-th direction; w (W) i Indicating the weight of the i-th direction.
The weight can adjust the duty ratio of deformation in different dimensions, for example, when the weight is focused on sedimentation, the sedimentation weight can be increased, and the GNSS station layout result can reflect more sedimentation space distribution. Meanwhile, in order to embody the constraint force of the weight, the following requirements should be satisfied:
Figure BDA0003919932250000111
s4, determining initial station measurement distribution: the initial station distribution comprises two types of stations, one type of station is a station which is determined to be arranged according to engineering and project requirements, and the station is not changed later; the other is a station requiring iterative optimization with an objective function, which should be approximately close to the final distribution to improve the efficiency and accuracy of the subsequent iterative solutions.
In some embodiments, step S4 may further include:
s41, considering the problem of model solving efficiency and accuracy compromise in the step S34, firstly determining an initial point location distribution which is as close to the optimal point location distribution as possible, and then carrying out iterative solving on the initial point location. The step therefore consists in determining an initial distribution of points.
S411, first, the following definition is performed:
from n candidate pointsDeleting one point in turn, and forming a point set by all the rest points
Figure BDA0003919932250000112
Then
Figure BDA0003919932250000113
Where i=1, 2, …, n, represents n such point sets; n-1 represents the remaining n-1 points in the point set; k+.i indicates that the ith point of the original n candidate points is deleted. Then the n Q point sets are used to form a new set F n Then
Figure BDA0003919932250000114
Wherein n represents the set F n With one element.
S412, performing Kerling interpolation on the n candidate points by adopting the variation function obtained in the step S31 to obtain an interpolated deformation field, and calculating the root mean square error of the interpolation field and InSAR deformation monitoring according to the calculation formula of the root mean square error in the step S33, wherein the root mean square error is denoted as r '' 0 。
S413, pair set F n Each point set in (2) is subjected to Kerling interpolation by adopting the variation function obtained in the step S31 to obtain an interpolation deformation field, and then the root mean square error of the deformation field and InSAR deformation monitoring is calculated according to a calculation formula of the root mean square error in the step S33 and is recorded as r '' n N represents the root mean square error result, where n is the number.
S414, calculating
r′＝r′ n -r′ 0
Recording the point set corresponding to the minimum value in r
Figure BDA0003919932250000115
And deleted candidate points (x i′ ,y i′ ,v i′ ). The minimum value represents that after deleting the point in all candidate points, the point has the least influence on the whole interpolation result Small, so that the point can be removed and the point set
Figure BDA0003919932250000116
As a new candidate point set.
S415, let n=n-1, and repeat steps S412 to S414 until n=3 (since at least 3 points are needed for kriging).
S416, finally obtaining the order of deleting each candidate point from the initial n candidate points. When the number N of GNSS stations in a research area is given, deleting the first N-N points and forming a point set Q 'consisting of the rest N points' N I.e. the initial start station for the area.
S42, determining a point set Q consisting of station positions which are required to be laid out according to engineering and project requirements M This part of the station position will not subsequently change.
S43, in order to ensure that the number of measuring stations is always N, a set Q 'is required to be collected at the initial point' T Find the distance Q respectively M Closest to each point in (a) and using Q M The point positions in the test station are replaced, so as to obtain an initial measuring station Q N 。
S5, solving an objective function model through multi-scale iterative optimization: and solving an objective function model by adopting a multi-scale iterative optimization method according to InSAR deformation, an initial station, a candidate point set and a variation function of a deformation field until the precision of model solving reaches the requirement.
In some embodiments, step S5 may further include:
s51, determining an initial measuring station Q through the steps S2-S4 N Further iterative optimization is required to fit as much higher deformation fields as possible. In order to improve the calculation efficiency in a large-scale scene, we propose a multi-scale iterative optimization method for final solution. Firstly, the initial candidate points are downscaled to reduce the number of the candidate points: increasing the variance threshold in step S2 to 3.5mm to 5mm, a new candidate point set D 'is obtained' n At this time D' n The density and number of (C) are far less than D n ；
S52, iteratively solving initial measuring station Q N At a small scalePoint selection set D' n The following more optimal solution X' N 。
S521 for initial measuring station Q N Middle-divided fixed point Q M Each point outside (assuming the current point is
Figure BDA0003919932250000121
) At a larger distance threshold T 1 =200m, at a small scale candidate point set D' n The search distance is less than T 1 And form a point set C n1 Where n1 represents that there are n1 points in the set of points.
S522, first, for initial station Q N Performing Kerling interpolation on the n candidate points in the step S31 by adopting a variation function obtained in the step S31 to obtain an interpolated deformation field, and then calculating the root mean square error of the deformation field and InSAR deformation monitoring according to a calculation formula of the root mean square error in the step S33, and marking the root mean square error as r' 0 。
S523, sequentially using C n1 Each of which replaces the current point
Figure BDA0003919932250000122
And performing Kriging interpolation by using the variation function obtained in the step S31 to obtain an interpolated deformation field, and then respectively calculating the root mean square error r' of the deformation field and InSAR deformation monitoring according to the calculation formula of the root mean square error in the step S33 n 。
S524, calculate
r″＝r″ n -r″ 0
Record the minimum value r 'in r' min And corresponding to
Figure BDA0003919932250000123
j represents the minimum value corresponding to C n1 Is the j-th point in (a). When r' min <0, description->
Figure BDA0003919932250000124
Point is +.>
Figure BDA0003919932250000125
Better, thus use the point +.>
Figure BDA0003919932250000126
Replace the anterior point->
Figure BDA0003919932250000127
Otherwise, the current point->
Figure BDA0003919932250000128
And keeping the iteration still, thus completing the iteration.
S525, let n=n+1 (encountering fixed point Q M Skip) and repeat steps S523-S524 until all points are iterated, and update the initial station to obtain updated station X '' N And let Q N ＝X′ N 。
S526, repeating the steps S521-S525 until the station Q is detected twice N The distance between the corresponding positions of each of the points in (a) is less than 10m. At this time last X' N I.e. solving the initial measuring station Q N At a small scale candidate point set D' n The following is a better solution.
S53, iteratively solving X' N At a large scale candidate point set D n Final solution X below N . With solving station X' N But requires the initial station Q N By X' N Instead, a small scale candidate point set D' n With a large set of candidate points D n Instead, a larger distance threshold T 1 A smaller distance threshold T for =200m 2 =100m instead. And repeating the steps S521-S526 to obtain the final GNSS station layout optimization result X N 。
Through the S1-S5, the advantages of large scale, high precision and high spatial resolution of InSAR in deformation monitoring are fully utilized to determine the layout of the GNSS station. The method can be used for rapidly giving out the reference number and the reference positions of the GNSS station arrangement aiming at a certain deformation field of InSAR monitoring, capturing deformation characteristics to the greatest extent, saving cost and improving the precision of the late InSAR and GNSS fusion.
Example 1
This example was performed using 29 scene Sentinel-1SAR data from 1 month in 2018 to 12 months in 2018 in a certain area.
FIG. 2 is a graph of deformation rate (in mm) obtained by the SBAS-InSAR method, with blank areas around the graph being areas of uncorrelated mask fall and unsuitable for station placement; FIG. 3 is a candidate point result obtained by the present method; FIGS. 4a and 4b show the initial station distribution and its deformation interpolation results obtained by the method, compared with InSAR monitoring results, the RMSE is 3.33mm; FIGS. 5a and 5b are the final station distribution and its deformation interpolation results after iterative optimization of the method, compared with InSAR monitoring results, the RMSE is 2.97mm;
from the results, the GNSS station layout obtained by the method can indeed reflect the spatial characteristics of deformation well, and compared with the maximum settlement of-50 mm, the final RMSE of the method is only 2.97, which proves that the station layout method provided by the method is practical and effective.
The above embodiments are merely preferred embodiments of the present invention, and the scope of the present invention is not limited to the above embodiments. All technical schemes belonging to the concept of the invention belong to the protection scope of the invention. It should be noted that modifications and adaptations to the present invention may occur to one skilled in the art without departing from the principles of the present invention and are intended to be within the scope of the present invention.

Claims (7)

1. The GNSS station measurement layout method based on InSAR deformation is characterized by comprising the following steps of:
s1, inSAR data processing: acquiring a differential interference atlas of radar images of a single track or a plurality of tracks, and calculating the average deformation rate of the earth surface in one dimension/two dimension/three dimension in a period of time by an InSAR technology;
s2, determining a candidate point set: taking all InSAR deformation points as an initial point set for GNSS station measurement layout, and selecting GNSS station measurement candidate points capable of showing deformation space characteristics of all deformation directions through quadtree downsampling so as to reduce the number and density of the initial point set; meanwhile, a mask area in which a measuring station cannot be arranged is obtained according to geological structure and topography external data of a research area, candidate points in the mask area are removed, and a final candidate point set is obtained;
s3, establishing an objective function model: firstly, estimating a variation function model of an InSAR deformation rate field, simultaneously taking a minimum prediction variance as a target, interpolating the whole deformation field through a kriging method according to a selected variation function and a measuring station, and establishing a model capable of minimizing root mean square error of interpolation deformation and InSAR monitoring deformation;
s4, determining initial station measurement distribution: the initial station distribution comprises two types of stations, one type of station is a station which is determined to be arranged according to engineering and project requirements, and the station is not changed later; the other is a station which needs to be iteratively optimized by adopting an objective function, and the initial position of the station is determined according to the change of root mean square by sequentially deleting candidate points;
S5, solving an objective function model through multi-scale iterative optimization: and solving an objective function model by adopting a multi-scale iterative optimization method according to InSAR deformation, an initial station, a candidate point set and a variation function of a deformation field until the precision of model solving reaches the requirement.
2. The GNSS station layout method based on InSAR deformation of claim 1, wherein S1 includes the sub-steps of:
s11, registering an acquired radar image, selecting a small baseline network, interfering, removing land and topography, filtering phases, unwrapping phases, correcting atmosphere and geocoding to form an interference atlas;
s12, assuming that only one space-time baseline set exists, obtaining time sequence deformation and average deformation rate by a short baseline set SBAS method through least square solution;
s13, further acquiring one-dimensional, two-dimensional or three-dimensional deformation according to the acquired data conditions of the research area:
one-dimensional deformation: if only one track exists, the line-of-sight LOS deformation monitored by SBAS-InSAR can be directly used to convert the line-of-sight LOS deformation into vertical deformation:
D U ＝D LOS /cosθ inc
wherein D is U Representing a vertical shapeVariable, D LOS LOS directional deformation, θ for SBAS-InSAR monitoring inc Representing the radar local angle of incidence;
Two-dimensional deformation: if the lifting rail data can be obtained, the two-dimensional deformation in the vertical direction and the east-west direction is obtained by fusing the lifting rails and neglecting the deformation in the north-south direction:
Figure FDA0003919932240000021
Figure FDA0003919932240000022
wherein D is U And D E Respectively the vertical deformation and the east-west deformation to be solved;
Figure FDA0003919932240000023
deformation indicative of monitoring of the lift rail data, +.>
Figure FDA0003919932240000024
Radar local angle of incidence for an orbiting satellite, +.>
Figure FDA0003919932240000025
Radar azimuth for an orbiting satellite; />
Figure FDA0003919932240000026
Representing the deformation of the derailment data monitoring, +.>
Figure FDA0003919932240000027
Radar local angle of incidence for a down-track satellite, +.>
Figure FDA0003919932240000028
Radar azimuth for the orbiting satellite;
the parameters in the above formula are directly solved by least square;
three-dimensional deformation: if the lifting rail with different sensors can be obtained, the deformation of the lifting rail in multiple directions is obtained by adopting a classical method, and finally the deformation in the east-west, south-north and vertical directions is solved by least square.
3. The GNSS station layout method based on InSAR deformation of claim 1, wherein S2 includes the sub-steps of:
s21, taking all InSAR deformation points as initial candidate points, respectively obtaining candidate point sets in different directions through quadtree downsampling, and then obtaining a union set of the candidate point sets in different directions to obtain a candidate point set of a GNSS measuring station;
s211, taking the InSAR deformation point set as an initial candidate point;
S212, setting a variance threshold T of deformation based on InSAR deformation;
s213, continuously subdividing the InSAR point set into squares by using a quadtree recursion algorithm, and calculating the variance of InSAR deformation in each square:
Figure FDA0003919932240000029
wherein sigma 2 Representing the variance of the InSAR points in the square, n represents the number of InSAR points in the square, and d i Representing the deformation of the i-th point in the square,
Figure FDA00039199322400000210
is the average value of n InSAR deformation points in the square;
if the variance is less than or equal to T, stopping segmentation; if the variance is greater than T, continuing to divide until the variance of deformation in each square is less than or equal to T; taking the point set of the quadtree downsampling as a candidate point set of the direction;
s214, obtaining a union set of candidate point sets in multiple directions to obtain an initial candidate point set of the whole area;
s22, obtaining available geological structure and topography external data of a research area, and masking off an area unsuitable for laying a measuring station;
s23, removing the candidate points in the mask area from the initial candidate point set to obtain a final candidate point set.
4. The GNSS station layout method based on InSAR deformation of claim 1, wherein S3 includes the sub-steps of: s31, calculating a variation value of the two-dimensional deformation field by using InSAR deformation points:
Figure FDA0003919932240000031
Wherein h is the distance between deformation points, also called the potential difference; n (h) is the paired number of all observation points with h as a distance; z (x) i ) And z (x) i +h) represent InSAR deformation observations at a relative distance h, respectively;
Figure FDA0003919932240000032
the variance of the interval h is a variation value; by this step the variance of InSAR deformation field is obtained>
Figure FDA0003919932240000033
Sample value from the distance h;
the aim is to establish the variance
Figure FDA0003919932240000034
And the spatial distance h between pairs of points, i.e., the variation function; according to the obtained sample value, selecting a model with the best fitting degree as a variation function, and fitting and calculating the model parameters;
s32, assuming that n candidate points are obtained in total in a certain region D in the step S2, using a set D n Indicating that it is now necessary to choose N points from the N candidate points to lay GNSS stations, the following set and symbol are defined:
with set P N Representing a set of these GNSS point distributions; by (x) i ,y i ,v i ) Indicating whenThe position of the i-th point in the front set and its deformation rate, where i=1, 2, … N, is simplified to
Figure FDA0003919932240000035
By S M Representing a series P with the same number of GNSS stations but different distribution N Wherein M represents a total of M sets P in the set N Then
Figure FDA0003919932240000036
Then->
Figure FDA0003919932240000037
Wherein->
Figure FDA0003919932240000038
And->
Figure FDA0003919932240000039
All represent the kth point set;
s33, for the M above P N Respectively performing Kerling interpolation to obtain an interpolated complete area deformation field; and then comparing the deformation field monitored by the InSAR with the deformation field monitored by the InSAR, and respectively calculating Root Mean Square Errors (RMSE) of each GNSS distribution set:
Figure FDA00039199322400000310
wherein d i Representing deformation of an ith point in the InSAR monitored deformation field;
Figure FDA00039199322400000311
representing the deformation of the ith point of the kriging interpolation; n (D) represents the number of InSAR points in the monitoring area D; RMSE k Representation set S M Root mean square error of the kth point set;
s34, taking the RMSE calculated in the step S33 as an index, converting the position solving problem of the GNSS station arrangement into the following RMSE minimization solving problem:
Figure FDA00039199322400000312
wherein W is a weight or a scale factor, and R represents root mean square errors of different station arrangements;
when the two-dimensional or three-dimensional deformation of the research area is obtained under the condition, the deformation of each dimension is considered, and the above expression is converted into the following form:
Figure FDA00039199322400000313
where ndim=1, 2,3, respectively representing the one-, two-and three-dimensional deformations of the investigation region obtained; r is R i Root mean square error representing the i-th direction; w (W) i A weight representing the i-th direction;
the weight can adjust the duty ratio of deformation in different dimensions, if the weight of sedimentation is increased when the weight is focused on sedimentation, the GNSS station arrangement result can reflect more sedimentation space distribution; meanwhile, in order to embody the constraint force of the weight, the following requirements should be satisfied:
Figure FDA0003919932240000041
5. The InSAR deformation-based GNSS station layout method of claim 4 wherein S4 includes the sub-steps of:
s41, considering the problem of model solving efficiency and accuracy compromise in the step S34, firstly determining an initial point location distribution which is as close to the optimal point location distribution as possible, and then carrying out iterative solving on the initial point location;
s411, first, the following definition is performed:
one point is deleted in turn from the n candidate points,all the remaining points form a point set
Figure FDA0003919932240000042
Then:
Figure FDA0003919932240000043
where i=1, 2, …, n, represents n such point sets; n-1 represents the remaining n-1 points in the point set; k+.i indicates that the i-th point of the original n candidate points is deleted; then the n Q point sets are used to form a new set F n Then:
Figure FDA0003919932240000044
wherein n represents the set F n Has one element;
s412, performing Kerling interpolation on the n candidate points by adopting the variation function obtained in the step S31 to obtain an interpolated deformation field, and calculating the root mean square error of the interpolation field and InSAR deformation monitoring according to the calculation formula of the root mean square error in the step S33, wherein the root mean square error is denoted as r '' 0 ；
S413, pair set F n Each point set in (2) is subjected to Kerling interpolation by adopting the variation function obtained in the step S31 to obtain an interpolation deformation field, and then the root mean square error of the deformation field and InSAR deformation monitoring is calculated according to a calculation formula of the root mean square error in the step S33 and is recorded as r '' n N represents the root mean square error result with n;
s414, calculating
r′＝r′ n -r′ 0
Recording the point set corresponding to the minimum value in r
Figure FDA0003919932240000045
And deleted candidate points (x i′ ,y i′ ,v i′ ) The method comprises the steps of carrying out a first treatment on the surface of the The minimum value represents the pair of all candidate points after deleting the pointThe overall interpolation result has minimal impact, remove the point, and set the point +.>
Figure FDA0003919932240000046
As a new candidate point set;
s415, let n=n-1, and repeat steps S412 to S414 until n=3;
s416, finally obtaining the sequence of deleting each candidate point from the initial n candidate points; when the number N of GNSS stations in a research area is given, deleting the first N-N points and forming a point set Q 'consisting of the rest N points' N I.e. the initial start-up station for the area;
s42, determining a point set Q consisting of station positions which are required to be laid out according to engineering and project requirements M The position of the station is not changed later;
s43, in order to ensure that the number of measuring stations is always N, a set Q 'is required to be collected at the initial point' T Find the distance Q respectively M Closest to each point in (a) and using Q M The point positions in the test station are replaced, so as to obtain an initial measuring station Q N 。
6. The GNSS station layout method based on InSAR deformation of claim 4, wherein S5 includes the sub-steps of: :
s51, firstly, downscaling the initial candidate points to reduce the number of the candidate points: increasing the variance threshold in step S2 to a certain value, and obtaining a new candidate point set D' n At this time D' n The density and number of (C) are far less than D n ；
S52, iteratively solving initial measuring station Q N At a small scale candidate point set D' n The following more optimal solution X' N ；
S521 for initial measuring station Q N Middle-divided fixed point Q M Each point other than the current point is assumed to be
Figure FDA0003919932240000051
With a distance threshold T 1 At the small scale candidate point set D' n Searching inThe distance is less than T 1 And form a point set C n1 Wherein n1 represents that there are n1 points in the point set;
s522, first, for initial station Q N Performing Kerling interpolation on the n candidate points in the step S31 by adopting a variation function obtained in the step S31 to obtain an interpolated deformation field, and then calculating the root mean square error of the deformation field and InSAR deformation monitoring according to a calculation formula of the root mean square error in the step S33, and marking the root mean square error as r' 0 ；
S523, sequentially using C n1 Each of which replaces the current point
Figure FDA0003919932240000058
And performing Kriging interpolation by using the variation function obtained in the step S31 to obtain an interpolated deformation field, and then respectively calculating the root mean square error r' of the deformation field and InSAR deformation monitoring according to the calculation formula of the root mean square error in the step S33 n ；
S524 calculation
r″＝r″ n -r″ 0
Record the minimum value r 'in r' min And corresponding to
Figure FDA0003919932240000052
j represents the minimum value corresponding to C n1 The j-th point in (a); when r' min <0, description->
Figure FDA0003919932240000053
Point is +.>
Figure FDA0003919932240000054
Better, thus use the point +.>
Figure FDA0003919932240000055
Replace the anterior point- >
Figure FDA0003919932240000056
Otherwise, the current point->
Figure FDA0003919932240000057
The iteration is completed in such a way that the iteration is kept still;
s525, let n=n+1, encounters the fixed point Q M Skipping and repeating the steps S523-S524 until all points are iterated, and updating the initial station to obtain an updated station X'. N And let Q N ＝X′ N ；
S526, repeating the steps S521-S525 until the station Q is detected twice N The distance between the corresponding positions of each of the points in (a) is less than a given threshold; at this time last X' N I.e. solving the initial measuring station Q N At a small scale candidate point set D' n The following better solution;
s53, iteratively solving X' N At a large scale candidate point set D n Final solution X below N With solving the measuring station X' N But requires the initial station Q N By X' N Instead, a small scale candidate point set D' n With a large set of candidate points D n Instead, the distance threshold T 1 With a new distance threshold T 2 Instead, a new distance threshold T 2 Less than the distance threshold T 1 The method comprises the steps of carrying out a first treatment on the surface of the And then repeating the steps S521-S526 to solve the final GNSS station layout optimization result X N 。
7. Use of the InSAR deformation based GNSS station layout method of any of claims 1 to 6 in synthetic aperture interferometry and global navigation satellite systems.
CN202211356074.XA 2022-11-01 2022-11-01 GNSS station measurement layout method based on InSAR deformation and application thereof Pending CN116052009A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211356074.XA CN116052009A (en) 2022-11-01 2022-11-01 GNSS station measurement layout method based on InSAR deformation and application thereof

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211356074.XA CN116052009A (en) 2022-11-01 2022-11-01 GNSS station measurement layout method based on InSAR deformation and application thereof

Publications (1)

Publication Number Publication Date
CN116052009A true CN116052009A (en) 2023-05-02

Family

ID=86122459

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211356074.XA Pending CN116052009A (en) 2022-11-01 2022-11-01 GNSS station measurement layout method based on InSAR deformation and application thereof

Country Status (1)

Country Link
CN (1) CN116052009A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116659429A (en) * 2023-08-01 2023-08-29 齐鲁空天信息研究院 Multi-source data high-precision time sequence earth surface three-dimensional deformation resolving method and system

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116659429A (en) * 2023-08-01 2023-08-29 齐鲁空天信息研究院 Multi-source data high-precision time sequence earth surface three-dimensional deformation resolving method and system

Similar Documents

Publication Publication Date Title
CN110058236B (en) InSAR and GNSS weighting method oriented to three-dimensional surface deformation estimation
CN110111274B (en) Method for calibrating exterior orientation elements of satellite-borne push-broom optical sensor
CN114791273B (en) InSAR deformation monitoring result interpretation method for landslide
Takaku et al. Quality improvements of ‘AW3D’global DSM derived from ALOS prism
CN112284332B (en) High-rise building settlement monitoring result three-dimensional positioning method based on high-resolution INSAR
CN114689015B (en) Method for improving elevation precision of optical satellite stereoscopic image DSM
CN112882032B (en) Method and device for dynamically monitoring geological disaster SAR in key area of gas pipeline
CN115201825B (en) Atmospheric delay correction method in InSAR (interferometric synthetic aperture radar) inter-seismic deformation monitoring
CN116052009A (en) GNSS station measurement layout method based on InSAR deformation and application thereof
CN111964667B (en) geomagnetic-INS (inertial navigation System) integrated navigation method based on particle filter algorithm
CN102436679B (en) Medium-resolution remote sensing image discrete point DEM (Digital Elevation Model) construction method based on medium value filtering
CN116305902A (en) Flood maximum submerged depth space simulation method based on multi-mode remote sensing
CN113238228B (en) Three-dimensional earth surface deformation obtaining method, system and device based on level constraint
CN115457022B (en) Three-dimensional deformation detection method based on live-action three-dimensional model front-view image
CN116485857A (en) High-time-resolution glacier thickness inversion method based on multi-source remote sensing data
CN110927765A (en) Laser radar and satellite navigation fused target online positioning method
CN116310901A (en) Debris flow material source dynamic migration identification method based on low-altitude remote sensing
KR101941132B1 (en) Apparatus and method for extending available area of regional ionosphere map
CN105093222A (en) Automatic extraction method for block adjustment connection points of SAR image
CN115563582A (en) Multi-scale evaluation method, device and system for geological risk along oil and gas pipeline
Magnússon et al. The bedrock and tephra layer topography within the glacier filled Katla caldera, Iceland, deduced from dense RES-survey
CN109387872A (en) Surface-related multiple prediction technique
He et al. Robust Estimation of Landslide Displacement from Multi-temporal UAV Photogrammetry-Derived Point Clouds
Wu et al. Exploring Convolutional Neural Networks for Predicting Sentinel-C Backscatter between Image Acquisitions
Prakash et al. Assessment of Urban Built-Up Volume Using Geospatial Methods: A Case Study of Bangalore

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination