CN112050725A - Surface deformation monitoring method fusing InSAR and GPS - Google Patents

Surface deformation monitoring method fusing InSAR and GPS Download PDF

Info

Publication number
CN112050725A
CN112050725A CN202010957324.XA CN202010957324A CN112050725A CN 112050725 A CN112050725 A CN 112050725A CN 202010957324 A CN202010957324 A CN 202010957324A CN 112050725 A CN112050725 A CN 112050725A
Authority
CN
China
Prior art keywords
gps
insar
points
point
deformation
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
CN202010957324.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.)
Surveying And Mapping Institute Of Guangdong Nuclear Industry Geological Bureau
Original Assignee
Surveying And Mapping Institute Of Guangdong Nuclear Industry Geological Bureau
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 Surveying And Mapping Institute Of Guangdong Nuclear Industry Geological Bureau filed Critical Surveying And Mapping Institute Of Guangdong Nuclear Industry Geological Bureau
Priority to CN202010957324.XA priority Critical patent/CN112050725A/en
Publication of CN112050725A publication Critical patent/CN112050725A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/004—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring coordinates of points
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/86—Combinations of radar systems with non-radar systems, e.g. sonar, direction finder
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9094—Theoretical aspects
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/14—Receivers specially adapted for specific applications
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/40—Means for monitoring or calibrating

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Electromagnetism (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention relates to a surface deformation monitoring technology based on mutual fusion of InSAR and GPS. The method comprises the following steps: arranging GPS characteristic points and monitoring points; collecting SAR images and GPS signals of corresponding time; resolving three-dimensional coordinate information of the GPS points; speckle noise suppression of the SAR image; searching SAR image characteristic points and calculating the conversion relation between the SAR image characteristic points and GPS characteristic points; track error correction of the SAR satellite; forming an InSAR differential interference pattern; inversion and correction of atmospheric delay; phase unwrapping of the differential interferogram; the unwrapping phase is converted into deformation information; geocoding; and fusing deformation information of the InSAR and the GPS to obtain a deformation monitoring result. The method disclosed by the invention utilizes the mutual combination of InSAR and GPS technologies, breaks through the traditional single monitoring technology, and makes full use of the advantages of the two technologies to make up each other, so that the accuracy of the deformation data is more accurate and reliable.

Description

Surface deformation monitoring method fusing InSAR and GPS
Technical Field
The invention relates to a global positioning system and a synthetic aperture radar interferometry technique. In particular to a mode for improving insufficient monitoring information when the InSAR technology is used for monitoring surface deformation.
Background
The high time resolution and high precision earth surface three-dimensional monitoring is the advantage of the GPS technology, which can obtain high precision three-dimensional deformation information, however, the monitoring of the GPS technology is limited by the environment and cost, and the monitoring result is usually limited in spatial resolution and sparsely distributed, so the GPS technology can be ideally completed for small-range earth surface deformation monitoring, and for large-range monitoring, a large number of points need to be laid, which is not economically feasible.
The InSAR is an important space-to-ground monitoring technology at present, has all-weather all-day-long working capacity, and cannot be influenced by day and night and cloud layers, and meanwhile, the InSAR has the advantages of large coverage, convenience, rapidness, low cost and high spatial resolution, but is influenced by atmospheric delay, satellite orbit errors, surface coverage change and the like. The purpose of monitoring the ground surface deformation is to achieve a preventive effect, and in the InSAR technology, the problem of blurred vision direction exists, so the obtained result is one-dimensional deformation information of the vision direction, and in our life, many catastrophic deformations which occur in two-dimensional and three-dimensional are required to be monitored and prevented.
Disclosure of Invention
Combining the advantages and disadvantages of the two technologies, the invention solves the problems that: (1) correcting the orbit error of the InSAR technology by using GPS data; (2) correcting the atmospheric delay error of the InSAR technology by using GPS data; (3) correcting InSAR interference data by using a GPS; (4) InSAR and GPS fusion method.
(1) Track error of InSAR technology corrected by GPS data
In InSAR data processing, satellite orbit data affects important steps such as image registration and geocoding, so the accuracy of orbit data determines the accuracy of the result. In order to obtain high-precision orbit information, the invention utilizes InSAR visual line displacement to compare with GPS monitoring data, simulates an orbit error correction model and removes satellite orbit errors.
(2) Atmospheric delay error of InSAR technology corrected by GPS data
Deformation information can be obtained by carrying out differential processing on the interferogram, and the obtained deformation information can cause data errors due to relative atmospheric delay between two pixels on the SAR image and between the two SAR images. In order to eliminate such errors, the solution can be made by using the property of phase difference, and the difference of a certain fixed image element on the image reflects the phase difference and the surface deformation. Firstly, a random interpolation model is used for obtaining the estimation of a GPS troposphere, then an atmospheric delay error can be obtained through a double difference algorithm between time domains and stations, and an atmospheric delay error interpolation is obtained through an improved distance inverse ratio weighting interpolation method, so that the atmospheric delay error of InSAR data is removed.
(3) Correcting InSAR interferometric data using GPS
Nbc the unwrapping algorithm is a local unwrapping algorithm, the unwrapping precision is high, but the residual error points cannot be unwrapped; the Gud unwrapping algorithm is a global unwrapping algorithm, and unwrapping can be performed on each pixel, but unwrapping errors are propagated globally, so that unwrapping accuracy is low. In order to improve the unwrapping precision and the unwrapping range, the GPS characteristic point auxiliary InSAR phase unwrapping is carried out by combining Nbc unwrapping algorithm and Gud unwrapping algorithm.
(4) InSAR and GPS fusion method
Because the precision of the deformation field respectively established by the function model constraint and the random model constraint is not as high as that of the deformation field established by the dual-model constraint, the method adopts a mode of fusing InSAR and GPS data of the dual-constraint of the function model and the random model to establish the three-dimensional deformation field.
In order to solve the problems, the invention adopts the following technical scheme: an earth surface deformation monitoring method fusing an InSAR and a GPS is characterized by comprising the following steps:
step one, laying GPS characteristic points and monitoring points;
step two, collecting SAR images and GPS signals of corresponding time periods;
step three, resolving three-dimensional coordinate information of the GPS points;
step four, speckle noise suppression of the SAR image;
searching the SAR image feature points and calculating the conversion relation between the SAR image feature points and the GPS feature points;
sixthly, correcting the orbit error of the SAR satellite;
step seven, forming an InSAR differential interference pattern;
step eight, inverting and correcting atmospheric delay;
step nine, phase unwrapping of the differential interference pattern;
step ten, converting the unwrapping phase into deformation information;
step eleven, geocoding;
and step twelve, fusing deformation information of the InSAR and the GPS to obtain a deformation monitoring result.
Furthermore, in the first step, two types of GPS points are required to be distributed, and the two types of GPS points are divided into feature points and monitoring points; the characteristic points are used for unifying InSAR and GPS data to the same reference coordinate system and correcting a series of errors, so that the selected points have the characteristic of stability in a monitoring area; and the monitoring points are used for acquiring deformation information of the monitoring area.
Further, in the second step, in order to correct a series of errors of the InSAR data by using the GPS, a time period corresponding to the acquired GPS satellite signal needs to be synchronized with a time period of the SAR image.
Further, in the third step, the three-dimensional coordinate information of the distributed GPS characteristic points and the monitoring points is calculated according to the collected GPS satellite signals. Wherein the feature points require position information and elevation information; the monitoring points require position information and deformation information.
Furthermore, in the fourth step, in order to ensure the reliability of the result and need to suppress the speckle noise of the SAR image, the invention adopts wavelet transformation and fuzzy clustering to suppress the speckle noise of the image. The method mainly comprises the steps of firstly carrying out wavelet transform decomposition on an SAR image, dividing decomposed coefficients into two categories by utilizing fuzzy C-means clustering, wherein the first category is energy only containing speckle noise, the second category is energy containing signal energy, the first category is directly subjected to zero setting, and the second category is subjected to noise reduction, so that the purpose of inhibition is achieved.
Further, in the fifth step, according to the distributed GPS characteristic points, points corresponding to the GPS characteristic points on the SAR image are found, and a coordinate conversion relation between the points is established. At least three characteristic points are needed for achieving the purpose of coordinate conversion, and if the number of the characteristic points is more than three, the adjustment can be carried out by using a least square method.
Further, in the sixth step, the monitoring data of the GPS characteristic points and the sight direction data of the InSAR are used for comparative analysis, and a track error correction model is established, so that the track error in the InSAR data is removed. The specific operation is as follows:
converting the three-dimensional displacement of the GPS into the view direction of the InSAR by utilizing a coordinate conversion relation:
Figure BDA0002679027840000031
wherein d isGPSThe displacement of the GPS three-dimensional displacement in the sight line direction after coordinate conversion; n, E, U are components of the GPS three-dimensional displacement in the north-south, east-west and vertical directions respectively; a ishIs the flight direction of the flight platform; θ is the side view of the flight platform.
Assuming that n GPS feature points are distributed, the pixel positions of the n feature points in the image can be obtained through a coordinate conversion formula, which is denoted as I. Then the track error (d) of pixel iorbit(i) ) is:
dorbit(i)＝dInSAR(i)-dGPS(i) (2)
the error of fitting the orbit by the quadric surface model is as follows:
Figure BDA0002679027840000041
Figure BDA0002679027840000042
Figure BDA0002679027840000043
where v is a correction number containing atmospheric delay, terrain error, etc.; the position of a pixel in the image is represented by (i, j); and a is0、a1、a2、a3、a4、a5Are parameters of the fitted model.
Can be expressed in the form of a matrix:
V＝AX-l (6)
wherein the content of the first and second substances,
Figure BDA0002679027840000044
Figure BDA0002679027840000045
Figure BDA0002679027840000046
Figure BDA0002679027840000047
the weight matrix P can be expressed as:
Figure BDA0002679027840000048
wherein gamma is the coherence value of the pixel, and satisfies that gamma is more than 0 and less than 1, and the larger gamma is, the larger the weight is; suppose σ0Is an error in the unit weight, then
Figure BDA0002679027840000049
Pi＝γi。
Solving the least square solution as:
X＝(ATPA)-1APl (7)
the error estimate in the unit weight is:
Figure BDA00026790278400000410
and the allowable error is three times of the error in the unit weight, the value exceeding the allowable error is removed, and the removed value is iterated to obtain the final result, namely the pixel track error, so that the track error of the monitored area is obtained. Thus correcting InSAR data d after orbit errorcorrectionComprises the following steps:
dcorrection＝dInSAR-dorbit (9)
further, in a seventh step, interference flow processing is performed on the SAR image data obtained after the sixth step is completed, so as to obtain an interferogram.
Further, in step eight, deformation information may be obtained by performing difference processing on the interferogram in the previous step, and the obtained deformation information may cause an error of data due to relative atmospheric delay between two pixels on the SAR image and between two SAR images. In order to eliminate such errors, the solution can be made by using the property of phase difference, and the difference of a certain fixed image element on the image reflects the phase difference and the surface deformation. Firstly, a random interpolation model is used for obtaining estimation of a GPS troposphere, and then an atmospheric delay error can be obtained through a double difference algorithm between time domains and stations. After the atmospheric delay error is obtained, the atmospheric delay needs to be corrected. The distributed GPS characteristic points are utilized, an inverse distance ratio weighting interpolation method is adopted, and because the correlation between the atmospheric delay and the elevation is strong, the atmospheric delay correction is reversed by adopting an improved inverse distance ratio weighting interpolation method, and the atmospheric delay is estimated; and correcting atmospheric delay correction of the InSAR differential interferogram using the estimated atmospheric delay. The specific operation is as follows:
obtaining an inverse distance weighted interpolation formula of the atmospheric delay from GPS observation data:
Figure BDA0002679027840000051
wherein the coordinate is east
Figure BDA0002679027840000052
North China
Figure BDA0002679027840000053
Interpolated atmospheric delay of points of
Figure BDA0002679027840000054
Coordinate is east
Figure BDA0002679027840000055
North China
Figure BDA0002679027840000056
GPS delay correction of points to
Figure BDA0002679027840000057
Wherein
Figure BDA0002679027840000058
Representing weights associated with delays of GPS acquisition. N represents the number of GPS points in the vicinity of the interpolated predicted point.
The weight of the inverse distance weighted interpolation method is:
Figure BDA0002679027840000059
Figure BDA00026790278400000510
wherein the distance from the characteristic point i to the point to be interpolated is
Figure BDA00026790278400000511
Modified inverse distance weighted interpolation:
the above principle makes the image proportion of horizontal distance and elevation interpolation identical, however, in practice, the proportion of the two is different, so an influence proportion is added:
Figure BDA00026790278400000512
the horizontal distance of a data point to a point to be interpolated is represented as:
Figure BDA0002679027840000061
the elevation difference from the data point to the point to be interpolated is expressed as:
Figure BDA0002679027840000062
thus, the improved inverse distance weighted interpolation method after improvement is:
Figure BDA0002679027840000063
Figure BDA0002679027840000064
Figure BDA0002679027840000065
further, in step nine, there are many ways to correct InSAR phase unwrapping using the GPS control point. Nbc the unwrapping algorithm is a local unwrapping algorithm, the unwrapping precision is high, but the residual error points cannot be unwrapped; the Gud unwrapping algorithm is a global unwrapping algorithm, and unwrapping can be performed on each pixel, but unwrapping errors are propagated globally, so that unwrapping accuracy is low. In order to improve the unwrapping precision and the unwrapping range, the GPS characteristic point auxiliary InSAR phase unwrapping is carried out by combining the Nbc unwrapping algorithm and the Gud unwrapping algorithm, and the method specifically comprises the following operations:
(1) firstly, phase unwrapping is carried out through an Nbc unwrapping algorithm, and initial values of the whole-cycle unknowns of the island pixels with the non-residual error points and the GPS characteristic points can be obtained.
(2) And removing decimal place of the initial value obtained in the last step to obtain an integer, and taking the integer as the starting value of Gub unwrapping algorithm to correct the unknown number of the whole week. And obtaining the unwrapping values of the non-residual points and the island pixels with the GPS characteristic points.
(3) And interpolating the winding result obtained in the last step, and supplementing residual points and unwrapping values of island pixels without GPS characteristic points.
(4) Marking the non-residual points with the final unwrapping values and the island pixels with the GPS characteristic points obtained in the second step as fixed domains, and marking the residual points with the initial wrapping values and the island pixels without the GPS characteristic points obtained in the third step as non-fixed domains. The winding values from the second and third steps are combined and taken as a starting value to Gub for phase unwrapping in the unwrapping algorithm. Thus yielding our final unwrapped value.
Further, in step ten, the phase of unwrapping in step nine is converted into intuitive deformation information. The specific operation is as follows:
Figure BDA0002679027840000066
wherein Δ RdFor the amount of deformation of the line of sight, phidFor unwrapping phase, λ is the radar wavelength.
Further, in the eleventh step, the deformation result in the SAR coordinate system is projected to the direction of the experimental setup, so as to obtain the deformation amount in the sight line direction. The specific operation is as follows:
distance equation:
Figure BDA0002679027840000071
doppler equation:
Figure BDA0002679027840000072
an earth model equation:
Figure BDA0002679027840000073
wherein R represents SAR satellite sensor and groundThe distance between the points; s represents the satellite sensor t represents a ground point;
Figure BDA0002679027840000074
and
Figure BDA0002679027840000075
respectively representing the position vectors of the SAR satellite and the ground point;
Figure BDA0002679027840000076
represents a doppler shift;
Figure BDA0002679027840000077
representing the satellite velocity vector,
Figure BDA0002679027840000078
Representing a ground speed vector;
Figure BDA0002679027840000079
representing ground point location vectors
Figure BDA00026790278400000710
Projections on the x, y, z axes; a is the earth major semi-axis, b is the earth minor semi-axis, and h is the average normal height.
Further, in the twelfth step, a three-dimensional deformation field is established in a mode that a function model is fused with random model double-constrained InSAR and GPS data. The specific operation is as follows:
the amount of line-of-sight displacement of the synthetic aperture radar at a certain point P is denoted as LInSARAccording to the space position relation of InSAR technology imaging, the following observation equation can be obtained:
Figure BDA00026790278400000711
wherein, VInSARRepresents an observed value residual error; a. theInSAR＝[Sx Sy Sz]Is a design matrix; for three-dimensional displacement of point P
Figure BDA00026790278400000712
And (4) showing. And solving the three-dimensional deformation at the point P by using three or more than three line-of-sight observed values in the SAR interference pair.
Under the principle, the invention utilizes the GPS observation point as prior information to construct random model constraint. Recording the three-dimensional deformation monitoring value of the GPS characteristic point as
Figure BDA00026790278400000713
Will be provided with
Figure BDA00026790278400000714
As a pseudo-observation vector, the pseudo-observation vector is
Figure BDA00026790278400000715
Covariance matrix of
Figure BDA00026790278400000716
Note the book
Figure BDA00026790278400000717
Is a full matrix of
Figure BDA00026790278400000718
The error equation for the additional stochastic model constraint is then:
Figure BDA0002679027840000081
from the least squares it can be known that:
Figure BDA0002679027840000082
while
Figure BDA0002679027840000083
The co-factor matrix of (c) is:
Figure BDA0002679027840000084
wherein the content of the first and second substances,
Figure BDA0002679027840000085
Figure BDA0002679027840000086
in the formula (I), the compound is shown in the specification,
Figure BDA0002679027840000087
and representing the weight array corresponding to the observation value of the radar in the sight line direction.
In the InSAR technology, displacement information in the north-south direction is difficult to collect in the line of sight due to the influence of the flight mode of the flight platform. In order to make up for the defect, the invention utilizes the high-precision GPS characteristic point monitoring result to supplement the displacement information in the north-south direction. The error equation for the following additional functional model constraints can be listed:
Figure BDA0002679027840000088
where C is the coefficient matrix of the constraint equation and W is a constant term. The basic idea of least squares can be found:
Figure BDA0002679027840000089
Figure BDA00026790278400000810
Figure BDA00026790278400000811
the co-factor matrix of (c) is:
Figure BDA00026790278400000812
in the formula
Figure BDA00026790278400000813
Figure BDA00026790278400000814
KSRepresenting the joint coefficient vector.
In summary, the following error equation can be obtained by combining the stochastic model and the functional model dual constraint:
Figure BDA00026790278400000815
the basic idea of using least squares is:
Figure BDA00026790278400000816
Figure BDA00026790278400000817
Figure BDA00026790278400000818
the co-factor matrix of (c) is:
Figure BDA00026790278400000819
wherein the content of the first and second substances,
Figure BDA00026790278400000820
Figure BDA0002679027840000091
Figure BDA0002679027840000092
in summary, there are:
Figure BDA0002679027840000093
therefore, the method comprises the following steps:
Figure BDA0002679027840000094
and:
Figure BDA0002679027840000095
through the analysis, the random model and the function model are more accurate than a single model by using double constraints. According to the basic principle, the three-dimensional deformation of the earth surface is further solved.
Drawings
FIG. 1: data processing flow diagram of the invention
FIG. 2: flow chart for correcting track error of InSAR data by using GPS monitoring point data
FIG. 3: schematic diagram for correcting atmospheric delay error by utilizing GPS (global positioning system) to invert atmospheric delay
FIG. 4: graph of cause of unwinding island formation
FIG. 5: GPS-assisted InSAR phase unwrapping flow chart
DETAILED DESCRIPTION OF EMBODIMENT (S) OF INVENTION
Example 1
As shown in fig. 1-5, the present invention adopts the following technical solutions: an earth surface deformation monitoring method fusing an InSAR and a GPS is characterized by comprising the following steps:
step one, laying GPS characteristic points and monitoring points;
step two, collecting SAR images and GPS signals of corresponding time periods;
step three, resolving three-dimensional coordinate information of the GPS points;
step four, speckle noise suppression of the SAR image;
searching the SAR image feature points and calculating the conversion relation between the SAR image feature points and the GPS feature points;
sixthly, correcting the orbit error of the SAR satellite;
step seven, forming an InSAR differential interference pattern;
step eight, inverting and correcting atmospheric delay;
step nine, phase unwrapping of the differential interference pattern;
step ten, converting the unwrapping phase into deformation information;
step eleven, geocoding;
and step twelve, fusing deformation information of the InSAR and the GPS to obtain a deformation monitoring result.
Furthermore, in the first step, two types of GPS points are required to be distributed, and the two types of GPS points are divided into feature points and monitoring points; the characteristic points are used for unifying InSAR and GPS data to the same reference coordinate system and correcting a series of errors, so that the selected points have the characteristic of stability in a monitoring area; and the monitoring points are used for acquiring deformation information of the monitoring area.
Further, in the second step, in order to correct a series of errors of the InSAR data by using the GPS, a time period corresponding to the acquired GPS satellite signal needs to be synchronized with a time period of the SAR image.
Further, in the third step, the three-dimensional coordinate information of the distributed GPS characteristic points and the monitoring points is calculated according to the collected GPS satellite signals. Wherein the feature points require position information and elevation information; the monitoring points require position information and deformation information.
Furthermore, in the fourth step, in order to ensure the reliability of the result and need to suppress the speckle noise of the SAR image, the invention adopts wavelet transformation and fuzzy clustering to suppress the speckle noise of the image. The method mainly comprises the steps of firstly carrying out wavelet transform decomposition on an SAR image, dividing decomposed coefficients into two categories by utilizing fuzzy C-means clustering, wherein the first category is energy only containing speckle noise, the second category is energy containing signal energy, the first category is directly subjected to zero setting, and the second category is subjected to noise reduction, so that the purpose of inhibition is achieved.
Further, in the fifth step, according to the distributed GPS characteristic points, points corresponding to the GPS characteristic points on the SAR image are found, and a coordinate conversion relation between the points is established. At least three characteristic points are needed for achieving the purpose of coordinate conversion, and if the number of the characteristic points is more than three, the adjustment can be carried out by using a least square method.
Further, in the sixth step, the monitoring data of the GPS characteristic points and the sight direction data of the InSAR are used for comparative analysis, and a track error correction model is established, so that the track error in the InSAR data is removed. The specific operation is as follows:
converting the three-dimensional displacement of the GPS into the view direction of the InSAR by utilizing a coordinate conversion relation:
Figure BDA0002679027840000101
wherein d isGPSThe displacement of the GPS three-dimensional displacement in the sight line direction after coordinate conversion; n, E, U are components of the GPS three-dimensional displacement in the north-south, east-west and vertical directions respectively; a ishIs the flight direction of the flight platform; θ is the side view of the flight platform.
Assuming that n GPS feature points are distributed, the pixel positions of the n feature points in the image can be obtained through a coordinate conversion formula, which is denoted as I. Then the track error (d) of pixel iorbit(i) ) is:
dorbit(i)＝dInSAR(i)-dGPS(i) (2)
the error of fitting the orbit by the quadric surface model is as follows:
Figure BDA0002679027840000111
Figure BDA0002679027840000112
Figure BDA0002679027840000113
wherein v is a group containingCorrection numbers of atmospheric delay, terrain error, and the like; the position of a pixel in the image is represented by (i, j); and a is0、a1、a2、a3、a4、a5Are parameters of the fitted model.
Can be expressed in the form of a matrix:
V＝AX-l (6)
wherein the content of the first and second substances,
Figure BDA0002679027840000114
Figure BDA0002679027840000115
Figure BDA0002679027840000116
Figure BDA0002679027840000117
the weight matrix P can be expressed as:
Figure BDA0002679027840000118
wherein gamma is the coherence value of the pixel, and satisfies that gamma is more than 0 and less than 1, and the larger gamma is, the larger the weight is; suppose σ0Is an error in the unit weight, then
Figure BDA0002679027840000119
Pi＝γi。
Solving the least square solution as:
X＝(ATPA)-1APl (7)
the error estimate in the unit weight is:
Figure BDA00026790278400001110
the tolerance error is three times of the error in the unit weight, the value exceeding the tolerance error is eliminated, and the eliminated value is iterated to obtain the final result, i.e. the pixel trackAnd error, thereby obtaining the track error of the monitoring area. Thus correcting InSAR data d after orbit errorcorrectionComprises the following steps:
dcorrection＝dInSAR-dorbit (9)
further, in a seventh step, interference flow processing is performed on the SAR image data obtained after the sixth step is completed, so as to obtain an interferogram.
Further, in step eight, deformation information may be obtained by performing difference processing on the interferogram in the previous step, and the obtained deformation information may cause an error of data due to relative atmospheric delay between two pixels on the SAR image and between two SAR images. In order to eliminate such errors, the solution can be made by using the property of phase difference, and the difference of a certain fixed image element on the image reflects the phase difference and the surface deformation. Firstly, a random interpolation model is used for obtaining estimation of a GPS troposphere, and then an atmospheric delay error can be obtained through a double difference algorithm between time domains and stations. After the atmospheric delay error is obtained, the atmospheric delay needs to be corrected. The distributed GPS characteristic points are utilized, an inverse distance ratio weighting interpolation method is adopted, and because the correlation between the atmospheric delay and the elevation is strong, the atmospheric delay correction is reversed by adopting an improved inverse distance ratio weighting interpolation method, and the atmospheric delay is estimated; and correcting atmospheric delay correction of the InSAR differential interferogram using the estimated atmospheric delay. The specific operation is as follows:
obtaining an inverse distance weighted interpolation formula of the atmospheric delay from GPS observation data:
Figure BDA0002679027840000121
wherein the coordinate is east
Figure BDA0002679027840000122
North China
Figure BDA0002679027840000123
Interpolated atmospheric delay of points of
Figure BDA0002679027840000124
Coordinate is east
Figure BDA0002679027840000125
North China
Figure BDA0002679027840000126
GPS delay correction of points to
Figure BDA0002679027840000127
Wherein
Figure BDA0002679027840000128
Representing weights associated with delays of GPS acquisition. N represents the number of GPS points in the vicinity of the interpolated predicted point.
The weight of the inverse distance weighted interpolation method is:
Figure BDA0002679027840000129
Figure BDA00026790278400001210
wherein the distance from the characteristic point i to the point to be interpolated is
Figure BDA00026790278400001211
Modified inverse distance weighted interpolation:
the above principle makes the image proportion of horizontal distance and elevation interpolation identical, however, in practice, the proportion of the two is different, so an influence proportion is added:
Figure BDA0002679027840000131
the horizontal distance of a data point to a point to be interpolated is represented as:
Figure BDA0002679027840000132
the elevation difference from the data point to the point to be interpolated is expressed as:
Figure BDA0002679027840000133
thus, the improved inverse distance weighted interpolation method after improvement is:
Figure BDA0002679027840000134
Figure BDA0002679027840000135
Figure BDA0002679027840000136
further, in step nine, there are many ways to correct InSAR phase unwrapping using the GPS control point. Nbc the unwrapping algorithm is a local unwrapping algorithm, the unwrapping precision is high, but the residual error points cannot be unwrapped; the Gud unwrapping algorithm is a global unwrapping algorithm, and unwrapping can be performed on each pixel, but unwrapping errors are propagated globally, so that unwrapping accuracy is low. In order to improve the unwrapping precision and the unwrapping range, the GPS characteristic point auxiliary InSAR phase unwrapping is carried out by combining the Nbc unwrapping algorithm and the Gud unwrapping algorithm, and the method specifically comprises the following operations:
(1) firstly, phase unwrapping is carried out through an Nbc unwrapping algorithm, and initial values of the whole-cycle unknowns of the island pixels with the non-residual error points and the GPS characteristic points can be obtained.
(2) And removing decimal place of the initial value obtained in the last step to obtain an integer, and taking the integer as the starting value of Gub unwrapping algorithm to correct the unknown number of the whole week. And obtaining the unwrapping values of the non-residual points and the island pixels with the GPS characteristic points.
(3) And interpolating the winding result obtained in the last step, and supplementing residual points and unwrapping values of island pixels without GPS characteristic points.
(4) Marking the non-residual points with the final unwrapping values and the island pixels with the GPS characteristic points obtained in the second step as fixed domains, and marking the residual points with the initial wrapping values and the island pixels without the GPS characteristic points obtained in the third step as non-fixed domains. The winding values from the second and third steps are combined and taken as a starting value to Gub for phase unwrapping in the unwrapping algorithm. Thus yielding our final unwrapped value.
Further, in step ten, the phase of unwrapping in step nine is converted into intuitive deformation information. The specific operation is as follows:
Figure BDA0002679027840000141
wherein Δ RdFor the amount of deformation of the line of sight, phidFor unwrapping phase, λ is the radar wavelength.
Further, in the eleventh step, the deformation result in the SAR coordinate system is projected to the direction of the experimental setup, so as to obtain the deformation amount in the sight line direction. The specific operation is as follows:
distance equation:
Figure BDA0002679027840000142
doppler equation:
Figure BDA0002679027840000143
an earth model equation:
Figure BDA0002679027840000144
wherein R represents the distance between the SAR satellite sensor and a ground point; s represents the satellite sensor t represents a ground point;
Figure BDA0002679027840000145
and
Figure BDA0002679027840000146
respectively representing the position vectors of the SAR satellite and the ground point;
Figure BDA0002679027840000147
represents a doppler shift;
Figure BDA0002679027840000148
representing the satellite velocity vector,
Figure BDA0002679027840000149
Representing a ground speed vector;
Figure BDA00026790278400001410
representing ground point location vectors
Figure BDA00026790278400001411
Projections on the x, y, z axes; a is the earth major semi-axis, b is the earth minor semi-axis, and h is the average normal height.
Further, in the twelfth step, a three-dimensional deformation field is established in a mode that a function model is fused with random model double-constrained InSAR and GPS data. The specific operation is as follows:
the amount of line-of-sight displacement of the synthetic aperture radar at a certain point P is denoted as LInSARAccording to the space position relation of InSAR technology imaging, the following observation equation can be obtained:
Figure BDA00026790278400001412
wherein, VInSARRepresents an observed value residual error; a. theInSAR＝[Sx Sy Sz]Is a design matrix; for three-dimensional displacement of point P
Figure BDA00026790278400001413
And (4) showing. To solve the three-dimensional deformation at point P is advantageousAnd solving the sight direction observation value in three or more SAR interference pairs.
Under the principle, the invention utilizes the GPS observation point as prior information to construct random model constraint. Recording the three-dimensional deformation monitoring value of the GPS characteristic point as
Figure BDA00026790278400001414
Will be provided with
Figure BDA00026790278400001415
As a pseudo-observation vector, the pseudo-observation vector is
Figure BDA0002679027840000151
Covariance matrix of
Figure BDA0002679027840000152
Note the book
Figure BDA0002679027840000153
Is a full matrix of
Figure BDA0002679027840000154
The error equation for the additional stochastic model constraint is then:
Figure BDA0002679027840000155
from the least squares it can be known that:
Figure BDA0002679027840000156
while
Figure BDA0002679027840000157
The co-factor matrix of (c) is:
Figure BDA0002679027840000158
wherein the content of the first and second substances,
Figure BDA0002679027840000159
Figure BDA00026790278400001510
in the formula (I), the compound is shown in the specification,
Figure BDA00026790278400001511
and representing the weight array corresponding to the observation value of the radar in the sight line direction.
In the InSAR technology, displacement information in the north-south direction is difficult to collect in the line of sight due to the influence of the flight mode of the flight platform. In order to make up for the defect, the invention utilizes the high-precision GPS characteristic point monitoring result to supplement the displacement information in the north-south direction. The error equation for the following additional functional model constraints can be listed:
Figure BDA00026790278400001512
where C is the coefficient matrix of the constraint equation and W is a constant term. The basic idea of least squares can be found:
Figure BDA00026790278400001513
Figure BDA00026790278400001514
Figure BDA00026790278400001515
the co-factor matrix of (c) is:
Figure BDA00026790278400001516
in the formula
Figure BDA00026790278400001517
Figure BDA00026790278400001518
KSRepresenting the joint coefficient vector.
In summary, the following error equation can be obtained by combining the stochastic model and the functional model dual constraint:
Figure BDA00026790278400001519
the basic idea of using least squares is:
Figure BDA00026790278400001520
Figure BDA00026790278400001521
Figure BDA00026790278400001522
the co-factor matrix of (c) is:
Figure BDA0002679027840000161
wherein the content of the first and second substances,
Figure BDA0002679027840000162
Figure BDA0002679027840000163
Figure BDA0002679027840000164
in summary, there are:
Figure BDA0002679027840000165
therefore, the method comprises the following steps:
Figure BDA0002679027840000166
and:
Figure BDA0002679027840000167
example data
The InSAR sight line deformation and GPS point three-dimensional deformation data are shown in the following table 1:
TABLE 1
Figure BDA0002679027840000168
Figure BDA0002679027840000171
Algorithm scheme
For the simulation example, three schemes are adopted to carry out fusion of GPS and InSAR data to calculate the three-dimensional deformation.
According to the first scheme, a GPS monitoring value is used as a prior condition, an InSAR observation value is used as an observed quantity, and the GPS and InSAR data fusion is restrained by establishing a random model.
And secondly, establishing a function model to restrain the data fusion of the GPS and the InSAR by taking the observed value in the south-north direction of the GPS as a constraint condition
And thirdly, establishing a function model and a random model to perform double-constrained GPS and InSAR data fusion by taking the GPS monitoring value as a prior condition and the observation value in the north-south direction of the GPS as a constraint condition.
Algorithmic comparison
In the case of only InSAR data, the accuracy of the monitor values is estimated from the simulated InSAR data using a 5 x 5 square matrix, which is called direct decomposition. Under the condition of adding GPS data for fusion, the root mean square error of the three-dimensional deformation field established by the three schemes is obtained by using the basic principle of the step twelve, and the result is shown in the following table 2:
TABLE 2
Figure BDA0002679027840000172
It can be seen from the above data that the precision of the direct solution algorithm, the function model constraint, the random model constraint and the dual model constraint are consistent in the N direction. In the direction E, the accuracy of the function model constraint is lower than that of the direct method, because the accuracy of the direction E is greatly reduced by forcibly fixing the direction N under the algorithm. It can be seen from the table that the accuracy of the bi-model constraints, both in the whole and in all directions, is higher than that of the single model constraints and the direct solution. Under the double-model constraint, the precision of the E direction is improved by 18% compared with that of the direct InSAR solution, although the precision of the double-model constraint in the U direction is not greatly improved.
The method combines the advantages of InSAR and GPS technologies, breaks through the limitation of the traditional single surface deformation monitoring technology, and greatly improves the precision in the elevation direction.

Claims (9)

1. An earth surface deformation monitoring method fusing an InSAR and a GPS is characterized by comprising the following steps:
step one, laying GPS characteristic points and monitoring points;
step two, collecting SAR images and GPS signals of corresponding time periods;
step three, resolving three-dimensional coordinate information of the GPS points;
step four, speckle noise suppression of the SAR image;
searching the SAR image feature points and calculating the conversion relation between the SAR image feature points and the GPS feature points;
sixthly, correcting the orbit error of the SAR satellite;
step seven, forming an InSAR differential interference pattern;
step eight, inverting and correcting atmospheric delay;
step nine, phase unwrapping of the differential interference pattern;
step ten, converting the unwrapping phase into deformation information;
step eleven, geocoding;
and step twelve, fusing deformation information of the InSAR and the GPS to obtain a deformation monitoring result.
2. The earth surface deformation monitoring method fusing InSAR and GPS according to claim 1, characterized in that: in the step one, the characteristic points are used for unifying InSAR and GPS data to the same reference coordinate system and correcting a series of errors, so that the selected points have the characteristic of stability in a monitoring area; and the monitoring points are used for acquiring deformation information of the monitoring area.
3. The earth surface deformation monitoring method fusing InSAR and GPS according to claim 1, characterized in that: in the fourth step, in order to ensure the reliability of the result, the speckle noise of the SAR image needs to be suppressed, and the invention adopts wavelet transformation and fuzzy clustering to suppress the speckle noise of the image.
4. The earth surface deformation monitoring method fusing InSAR and GPS according to claim 1, characterized in that: and in the fifth step, according to the distributed GPS characteristic points, finding out points corresponding to the points on the SAR image, and establishing a conversion relation between the points.
5. The earth surface deformation monitoring method fusing InSAR and GPS according to claim 1, characterized in that: and in the sixth step, the monitoring data of the GPS characteristic points and the sight direction data of the InSAR are used for comparative analysis, and a track error correction model is established, so that the track error in the InSAR data is removed.
6. The earth surface deformation monitoring method fusing InSAR and GPS according to claim 1, characterized in that: in step eight, the estimation of the troposphere of the GPS is firstly obtained by using a random interpolation model, and then the atmospheric delay error is calculated by using the double difference between the sites and the time domain.
After the atmospheric delay error is obtained, the atmospheric delay needs to be corrected. The distributed GPS characteristic points are utilized, an inverse distance ratio weighted interpolation method is adopted, and because the correlation between the atmospheric delay and the elevation is strong, the atmospheric delay correction is reversed by adopting the improved inverse distance ratio weighted interpolation method, and the atmospheric delay is estimated; and correcting atmospheric delay correction of the InSAR differential interferogram using the estimated atmospheric delay. The specific operation is as follows:
obtaining an inverse distance weighted interpolation formula of the atmospheric delay from GPS observation data:
Figure FDA0002679027830000021
wherein the coordinate is east
Figure FDA0002679027830000022
North China
Figure FDA0002679027830000023
Interpolated atmospheric delay of points of
Figure FDA0002679027830000024
Coordinate is east
Figure FDA0002679027830000025
North China
Figure FDA0002679027830000026
GPS delay correction of points to
Figure FDA0002679027830000027
Wherein
Figure FDA0002679027830000028
Representing weights associated with delays of GPS acquisition. N represents the number of GPS points in the vicinity of the interpolated predicted point.
The weight of the inverse distance weighted interpolation method is:
Figure FDA0002679027830000029
and
Figure FDA00026790278300000210
Figure FDA00026790278300000211
wherein the distance from the characteristic point i to the point to be interpolated is
Figure FDA00026790278300000212
Modified inverse distance weighted interpolation:
the above principle makes the image proportion of horizontal distance and elevation interpolation identical, however, in practice, the proportion of the two is different, so an influence proportion is added:
Figure FDA00026790278300000213
the horizontal distance of a data point to a point to be interpolated is represented as:
Figure FDA00026790278300000214
the elevation difference from the data point to the point to be interpolated is expressed as:
Figure FDA00026790278300000215
thus, the improved inverse distance weighted interpolation method after improvement is:
Figure FDA00026790278300000216
Figure FDA00026790278300000217
and
Figure FDA00026790278300000218
Figure FDA0002679027830000031
7. the earth surface deformation monitoring method fusing InSAR and GPS according to claim 1, characterized in that: in step nine, the GPS characteristic point auxiliary InSAR phase unwrapping is carried out by combining Nbc unwrapping algorithm and Gud unwrapping algorithm.
8. The earth surface deformation monitoring method fusing InSAR and GPS according to claim 1, characterized in that: in the eleventh step, the deformation result under the SAR coordinate system is projected to the direction set by the experiment to obtain the deformation amount in the sight line direction. The specific operation is as follows:
distance equation:
Figure FDA0002679027830000032
doppler equation:
Figure FDA0002679027830000033
an earth model equation:
Figure FDA0002679027830000034
wherein R represents the distance between the SAR satellite sensor and a ground point; s stands for satelliteSensor t represents a ground point;
Figure FDA0002679027830000035
and
Figure FDA0002679027830000036
respectively representing the position vectors of the SAR satellite and the ground point;
Figure FDA0002679027830000037
represents a doppler shift;
Figure FDA0002679027830000038
representing the satellite velocity vector,
Figure FDA0002679027830000039
Representing a ground speed vector;
Figure FDA00026790278300000310
representing ground point location vectors
Figure FDA00026790278300000311
Projections on the x, y, z axes; a is the earth major semi-axis, b is the earth minor semi-axis, and h is the average normal height.
9. The earth surface deformation monitoring method fusing InSAR and GPS according to claim 1, characterized in that: in the step twelve, the precision of the three-dimensional deformation field established by adopting a mode of fusing the InSAR and GPS data of double constraints of the function model and the random model is higher than that of the deformation field established by single constraints of the function model and the random model. The specific operation is as follows:
the amount of line-of-sight displacement of the synthetic aperture radar at a certain point P is denoted as LInSARAccording to the space position relation of InSAR technology imaging, the following observation equation can be obtained:
Figure FDA00026790278300000312
wherein, VInSARRepresents an observed value residual error; a. theInSAR＝[Sx Sy Sz]Is a design matrix; for three-dimensional displacement of point P
Figure FDA00026790278300000313
And (4) showing. And solving the three-dimensional deformation at the point P by using three or more than three line-of-sight observed values in the SAR interference pair.
Under the principle, the invention utilizes the GPS observation point as prior information to construct random model constraint. Recording the three-dimensional deformation monitoring value of the GPS characteristic point as
Figure FDA0002679027830000041
Will be provided with
Figure FDA0002679027830000042
As a pseudo-observation vector, the pseudo-observation vector is
Figure FDA0002679027830000043
Covariance matrix of
Figure FDA0002679027830000044
Note the book
Figure FDA0002679027830000045
Is a full matrix of
Figure FDA0002679027830000046
The error equation for the additional stochastic model constraint is then:
Figure FDA0002679027830000047
from the least squares it can be known that:
Figure FDA0002679027830000048
while
Figure FDA0002679027830000049
The co-factor matrix of (c) is:
Figure FDA00026790278300000410
wherein the content of the first and second substances,
Figure FDA00026790278300000411
in the formula (I), the compound is shown in the specification,
Figure FDA00026790278300000412
and representing the weight array corresponding to the observation value of the radar in the sight line direction.
In the InSAR technology, displacement information in the north-south direction is difficult to collect in the line of sight due to the influence of the flight mode of the flight platform. In order to make up for the defect, the invention utilizes the high-precision GPS characteristic point monitoring result to supplement the displacement information in the north-south direction. The error equation for the following additional functional model constraints can be listed:
Figure FDA00026790278300000413
where C is the coefficient matrix of the constraint equation and W is a constant term. The basic idea of least squares can be found:
Figure FDA00026790278300000414
Figure FDA00026790278300000415
the co-factor matrix of (c) is:
Figure FDA00026790278300000417
in the formula
Figure FDA00026790278300000418
KSRepresenting the joint coefficient vector.
In summary, the following error equation can be obtained by combining the stochastic model and the functional model dual constraint:
Figure FDA0002679027830000051
the basic idea of using least squares is:
Figure FDA0002679027830000052
Figure FDA0002679027830000053
Figure FDA0002679027830000054
the co-factor matrix of (c) is:
Figure FDA0002679027830000055
wherein the content of the first and second substances,
Figure FDA0002679027830000056
Figure FDA0002679027830000057
Figure FDA0002679027830000058
in summary, there are:
Figure FDA0002679027830000059
therefore, the method comprises the following steps:
Figure FDA00026790278300000510
and:
Figure FDA00026790278300000511
through the analysis, the random model and the function model are more accurate than a single model by using double constraints. According to the basic principle, the three-dimensional deformation of the earth surface is further solved.
CN202010957324.XA 2020-09-14 2020-09-14 Surface deformation monitoring method fusing InSAR and GPS Pending CN112050725A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010957324.XA CN112050725A (en) 2020-09-14 2020-09-14 Surface deformation monitoring method fusing InSAR and GPS

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010957324.XA CN112050725A (en) 2020-09-14 2020-09-14 Surface deformation monitoring method fusing InSAR and GPS

Publications (1)

Publication Number Publication Date
CN112050725A true CN112050725A (en) 2020-12-08

Family

ID=73610759

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010957324.XA Pending CN112050725A (en) 2020-09-14 2020-09-14 Surface deformation monitoring method fusing InSAR and GPS

Country Status (1)

Country Link
CN (1) CN112050725A (en)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112711022A (en) * 2020-12-18 2021-04-27 中国矿业大学 GNSS chromatography-assisted InSAR (interferometric synthetic aperture radar) atmospheric delay correction method
CN112986993A (en) * 2021-02-07 2021-06-18 同济大学 InSAR deformation monitoring method based on space constraint
CN114002719A (en) * 2021-10-12 2022-02-01 广州市城市规划勘测设计研究院 Single-frequency dual-mode multi-antenna fusion RTK positioning method, device, equipment and medium
CN115201825A (en) * 2022-09-16 2022-10-18 眉山环天智慧科技有限公司 Atmospheric delay correction method in InSAR (interferometric synthetic aperture radar) inter-seismic deformation monitoring
WO2022214114A3 (en) * 2021-08-10 2022-11-10 中咨数据有限公司 Bridge deformation monitoring method fusing gnss data and insar technology
CN115358311A (en) * 2022-08-16 2022-11-18 南昌大学 Multi-source data fusion processing method for surface deformation monitoring
CN117109426A (en) * 2023-08-28 2023-11-24 兰州交通大学 Three-dimensional deformation field modeling method fusing GNSS/InSAR observation data

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH07199804A (en) * 1993-12-28 1995-08-04 Nec Corp Topographical map generating device employing three-dimensional information obtained by interference type synthetic aperture radar
US5659318A (en) * 1996-05-31 1997-08-19 California Institute Of Technology Interferometric SAR processor for elevation
CN101770027A (en) * 2010-02-05 2010-07-07 河海大学 Ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion
CN107389029A (en) * 2017-08-24 2017-11-24 北京市水文地质工程地质大队 A kind of surface subsidence integrated monitor method based on the fusion of multi-source monitoring technology
CN110058237A (en) * 2019-05-22 2019-07-26 中南大学 InSAR point Yun Ronghe and three-dimensional deformation monitoring method towards High-resolution SAR Images

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH07199804A (en) * 1993-12-28 1995-08-04 Nec Corp Topographical map generating device employing three-dimensional information obtained by interference type synthetic aperture radar
US5659318A (en) * 1996-05-31 1997-08-19 California Institute Of Technology Interferometric SAR processor for elevation
CN101770027A (en) * 2010-02-05 2010-07-07 河海大学 Ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion
CN107389029A (en) * 2017-08-24 2017-11-24 北京市水文地质工程地质大队 A kind of surface subsidence integrated monitor method based on the fusion of multi-source monitoring technology
CN110058237A (en) * 2019-05-22 2019-07-26 中南大学 InSAR point Yun Ronghe and three-dimensional deformation monitoring method towards High-resolution SAR Images

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
马艳鸽等: "函数模型与随机模型双约束的ＧＰＳ－ＩｎＳＡＲ数据融合方法建立三维形变场", 《大地测量与地球动力学》 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112711022A (en) * 2020-12-18 2021-04-27 中国矿业大学 GNSS chromatography-assisted InSAR (interferometric synthetic aperture radar) atmospheric delay correction method
CN112711022B (en) * 2020-12-18 2022-08-30 中国矿业大学 GNSS chromatography-assisted InSAR (interferometric synthetic aperture radar) atmospheric delay correction method
CN112986993A (en) * 2021-02-07 2021-06-18 同济大学 InSAR deformation monitoring method based on space constraint
WO2022214114A3 (en) * 2021-08-10 2022-11-10 中咨数据有限公司 Bridge deformation monitoring method fusing gnss data and insar technology
CN114002719A (en) * 2021-10-12 2022-02-01 广州市城市规划勘测设计研究院 Single-frequency dual-mode multi-antenna fusion RTK positioning method, device, equipment and medium
CN115358311A (en) * 2022-08-16 2022-11-18 南昌大学 Multi-source data fusion processing method for surface deformation monitoring
CN115358311B (en) * 2022-08-16 2023-06-16 南昌大学 Multi-source data fusion processing method for ground surface deformation monitoring
CN115201825A (en) * 2022-09-16 2022-10-18 眉山环天智慧科技有限公司 Atmospheric delay correction method in InSAR (interferometric synthetic aperture radar) inter-seismic deformation monitoring
CN115201825B (en) * 2022-09-16 2023-01-17 眉山环天智慧科技有限公司 Atmospheric delay correction method in InSAR (interferometric synthetic aperture radar) inter-seismic deformation monitoring
CN117109426A (en) * 2023-08-28 2023-11-24 兰州交通大学 Three-dimensional deformation field modeling method fusing GNSS/InSAR observation data
CN117109426B (en) * 2023-08-28 2024-03-22 兰州交通大学 Three-dimensional deformation field modeling method fusing GNSS/InSAR observation data

Similar Documents

Publication Publication Date Title
CN112050725A (en) Surface deformation monitoring method fusing InSAR and GPS
CN107102333B (en) Satellite-borne InSAR long and short baseline fusion unwrapping method
Rabus et al. The shuttle radar topography mission—a new class of digital elevation models acquired by spaceborne radar
CN113624122B (en) Bridge deformation monitoring method fusing GNSS data and InSAR technology
Crosetto Calibration and validation of SAR interferometry for DEM generation
CN113340191B (en) Time series interference SAR deformation quantity measuring method and SAR system
CN111273293B (en) InSAR residual motion error estimation method and device considering terrain fluctuation
CN107144293A (en) A kind of geometric calibration method of video satellite area array cameras
CN106960174A (en) High score image laser radar vertical control point is extracted and its assisted location method
CN112711021B (en) Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method
CN113514829A (en) InSAR-oriented initial DSM block adjustment method
CN112986949B (en) SAR high-precision time sequence deformation monitoring method and device for diagonal reflector
CN113189551B (en) GB-InSAR heavy rail error compensation method based on scene DEM
KR20120009186A (en) method for manufacturing a digital elevation model using a SAR data
CN115469308A (en) Multi-track InSAR inter-seismic deformation rate field splicing method, device, equipment and medium
CN115097450A (en) Cross-track high-resolution three-SAR (synthetic aperture radar) offset large-gradient landslide deformation estimation method
CN102944308B (en) Attitude error correcting method of time-space joint modulation interference imaging spectrometer
CN112835043B (en) Method for monitoring two-dimensional deformation in any direction
Li et al. Improve the ZY-3 height accuracy using ICESat/GLAS laser altimeter data
CN115343709B (en) Terrain inversion method suitable for distributed spaceborne SAR system
CN114488144A (en) SAR offset three-dimensional deformation estimation method and system with additional DEM constraint
CN114280608A (en) Method and system for removing DInSAR elevation-related atmospheric effect
CN111580101A (en) InSAR baseline error uncontrolled correction method and device based on external DEM
CN115201823B (en) Ground surface deformation monitoring method utilizing BDS-InSAR data fusion
Lombardi et al. Accuracy of high resolution CSK interferometric Digital Elevation Models

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
RJ01 Rejection of invention patent application after publication

Application publication date: 20201208

RJ01 Rejection of invention patent application after publication