CN101770027B - Ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion - Google Patents
Ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion Download PDFInfo
- Publication number
- CN101770027B CN101770027B CN2010101067941A CN201010106794A CN101770027B CN 101770027 B CN101770027 B CN 101770027B CN 2010101067941 A CN2010101067941 A CN 2010101067941A CN 201010106794 A CN201010106794 A CN 201010106794A CN 101770027 B CN101770027 B CN 101770027B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msubsup
- msup
- gps
- mtd
- 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.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 89
- 238000012544 monitoring process Methods 0.000 title claims abstract description 46
- 230000004927 fusion Effects 0.000 title claims abstract description 12
- 238000006243 chemical reaction Methods 0.000 claims abstract description 9
- 238000004422 calculation algorithm Methods 0.000 claims description 46
- 238000012937 correction Methods 0.000 claims description 42
- 239000013598 vector Substances 0.000 claims description 21
- 238000005259 measurement Methods 0.000 claims description 20
- 239000011159 matrix material Substances 0.000 claims description 18
- 230000008569 process Effects 0.000 claims description 16
- 230000000007 visual effect Effects 0.000 claims description 8
- 238000007500 overflow downdraw method Methods 0.000 claims description 7
- 230000001629 suppression Effects 0.000 claims description 7
- 238000013507 mapping Methods 0.000 claims description 6
- 238000005457 optimization Methods 0.000 claims description 5
- 238000010587 phase diagram Methods 0.000 claims description 5
- 238000002790 cross-validation Methods 0.000 claims description 4
- 238000001914 filtration Methods 0.000 claims description 4
- 239000005436 troposphere Substances 0.000 claims description 4
- 238000004458 analytical method Methods 0.000 claims description 3
- 230000001427 coherent effect Effects 0.000 claims description 3
- 230000001360 synchronised effect Effects 0.000 claims description 3
- 238000005516 engineering process Methods 0.000 abstract description 22
- 238000004364 calculation method Methods 0.000 description 12
- 238000012545 processing Methods 0.000 description 10
- 230000008901 benefit Effects 0.000 description 9
- 238000011156 evaluation Methods 0.000 description 6
- 238000004804 winding Methods 0.000 description 6
- 230000008859 change Effects 0.000 description 5
- 230000009466 transformation Effects 0.000 description 5
- 238000000137 annealing Methods 0.000 description 4
- 238000001514 detection method Methods 0.000 description 4
- 239000000047 product Substances 0.000 description 4
- 230000002123 temporal effect Effects 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 3
- 230000001934 delay Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000012952 Resampling Methods 0.000 description 2
- CLOMYZFHNHFSIQ-UHFFFAOYSA-N clonixin Chemical compound CC1=C(Cl)C=CC=C1NC1=NC=CC=C1C(O)=O CLOMYZFHNHFSIQ-UHFFFAOYSA-N 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000006073 displacement reaction Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000012886 linear function Methods 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 238000004062 sedimentation Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 238000002592 echocardiography Methods 0.000 description 1
- 239000012467 final product Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000005305 interferometry Methods 0.000 description 1
- 230000033001 locomotion Effects 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000010845 search algorithm Methods 0.000 description 1
- 238000002922 simulated annealing Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
Images
Landscapes
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
The invention discloses a ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion, which comprises the following steps: 1, distributing GPS points; 2, collecting GPS signals and collecting SAR images; 3, resolving GPS points; 4, inhabiting SAR image spot noise; 5, establishing the coordinate conversion relationship between the GPS and SAR data; 6, utilizing GPS data for modifying SAR satellite trajectory errors; 7, forming an InSAR differential interference map; 8, utilizing GPS points for inverting the atmospheric delay; 9, modifying the InSAR atmospheric delay errors; 10, unwrapping the InSAR differential interference map; 11, converting the unwrapping phase into deformation information; 12, carrying out geographical coding on the InSAR deformation information; 13, fusing InSAR and GPS deformation information to obtain high-precision ground surface three-dimensional deformation results; and 14, utilizing the Kalman wave filter estimation and obtaining the ground surface three-dimensional deformation results with high spatial resolution. Because the invention adopts the ground surface three-dimensional deformation monitoring technology with the fused InSAR and GPS data, the application limitation of the single monitoring technology can be broken through, and the spatial resolution of the ground surface three-dimensional deformation monitoring results is greatly improved. The monitoring precision, particularly the precision in the vertical direction is high, and at the same time, the integral reliability of the system is enhanced.
Description
Technical Field
The invention relates to a method for monitoring three-dimensional deformation of a ground surface in a space-to-ground mode, in particular to a method for monitoring three-dimensional deformation of the ground surface based on InSAR and GPS data fusion, which is suitable for monitoring three-dimensional deformation of the ground surface in regions such as urban environment, landslide and earthquake with high spatial-temporal resolution and high precision.
Background
In recent years, with the increasingly mature satellite positioning theory and the improvement of the performance of receiving devices, Global Positioning System (GPS) has been widely used for deformation monitoring. The GPS can provide a three-dimensional deformation measurement value with high time resolution and high precision, is not influenced by weather, does not need to look through between point locations, is flexible in station arrangement, and has strong adaptability to earth surface conditions. However, due to the limitation of the number of receivers and the network arrangement array, the GPS has low spatial resolution and small coverage area, and some serious potential safety hazards may be missed. Specifically, the existing GPS deformation monitoring method mainly has the following disadvantages: (1) the method of net arrangement observation is adopted, the labor intensity is high, and the cost of a monitoring system is high; (2) by adopting a point observation method, the spatial resolution of the monitoring data is low, and the requirements of large-area deformation analysis and geological disaster prediction are difficult to meet; (3) in a deep mountain canyon or an urban building dense area, the GPS receiver antenna is shielded, so that the number of received GPS satellites is reduced, the GPS positioning precision is greatly reduced, and the requirement on deformation observation precision cannot be met.
Synthetic aperture radar interferometry (InSAR) is a new surface detection technology which is rapidly developed in recent years, and more countries and regions utilize the InSAR differential technology to detect ground deformation phenomena caused by mining, earthquakes, volcanic motions and the like. Compared with the GPS technology, the InSAR technology has the advantages of high vertical monitoring precision, large monitoring range, near space continuity and the like. However, the InSAR technology is low in time resolution and very sensitive to factors such as atmospheric delay errors, satellite orbit errors, earth surface conditions and temporal incoherence, which causes difficulty in application of the InSAR technology to earth surface deformation detection.
Disclosure of Invention
The invention aims to overcome the defects in the prior art, comprehensively utilize InSAR and GPS technologies, break through the application limitation of a single technology, and provide a surface three-dimensional deformation monitoring method based on the fusion of InSAR and GPS data.
The basic principle of the invention is as follows:
the InSAR technology is a new earth surface detection technology which is rapidly developed in recent years, and the principle of the InSAR technology is to monitor the deformation of the earth surface by processing received earth surface echoes with different time phases. Compared with the GPS technology, the InSAR technology has the advantages of high vertical monitoring precision, large monitoring range, near space continuity and the like. However, the InSAR technology is low in time resolution and very sensitive to factors such as atmospheric delay errors, satellite orbit errors, earth surface conditions and temporal incoherence, which causes difficulty in application of the InSAR technology to earth surface deformation detection.
The invention considers the fusion of GPS and InSAR data, breaks through the limitation of single technology application, exerts the respective advantages thereof and greatly improves the resolution capability of a space domain and a time domain. However, the data fusion of the GPS and the InSAR has the following key technologies to be solved: (1) correcting the orbit error of the InSAR satellite by using the GPS; (2) correcting the atmospheric delay influence of InSAR data by using a GPS; (3) GPS-assisted InSAR phase unwrapping; (4) and fusing InSAR and GPS data in a time-space domain.
In consideration of the fact that a high-precision monitoring method needs to be invented, errors of observed data need to be corrected, GPS data are fused according to an InSAR data processing flow, the InSAR errors corrected by the GPS are utilized, then a GPS and InSAR fusion method based on an improved Markov random field is adopted, fusion of the InSAR and the GPS data in a time-space domain is achieved, and surface three-dimensional deformation information with high space-time resolution is obtained.
1. GPS correction of orbital errors in InSAR data
The orbit information of the SAR satellite is used in the steps of InSAR image registration, phase-to-elevation conversion, geocoding and the like. Therefore, obtaining accurate orbit information is important for improving the precision of InSAR measurement. Because the precision of the orbit information of the InSAR data is not very high, for example, the precision of the basic orbit data of ERS cannot meet the requirement of measurement precision in the flight direction of 2-4m, the vertical track direction of 1-2m and the radial direction of 5m, the invention introduces GPS observation data, accurately measures the ground coordinates of characteristic ground objects (such as intersections of rivers and roads, water towers and the like) or artificial corner reflectors on the InSAR image by using GPS, and then obtains the accurate orbit coordinates of SAR satellites by using the reverse process of geocoding, thereby achieving the purpose of correcting the orbit error of the InSAR data.
2. GPS inversion atmospheric delay correction InSAR atmospheric delay error
When the electromagnetic wave emitted by the SAR sensor passes through the atmosphere, the propagation speed changes, and the propagation path is also bent, so that a delay phenomenon is generated. The atmospheric delays of SAR images observed in different time in the same region are not completely the same, so that the InSAR interferometric phase diagram is inevitably polluted. Therefore, the atmospheric delay is inverted by the GPS, the troposphere estimation of the GPS is firstly obtained by adopting a random process method, then the inter-station time domain bidyne method is adopted to estimate the atmospheric delay correction, the strong correlation between the atmospheric delay and the elevation is considered, an improved inverse distance weighted interpolation method which integrates the GPS elevation information is proposed to encrypt the atmospheric delay correction according to the theory of the geographic space correlation, and then the encrypted atmospheric delay correction is eliminated from the InSAR interference phase diagram, so that the InSAR atmospheric delay error is corrected.
3. GPS assisted InSAR phase unwrapping
In order to convert the phase into elevation or deformation information, the winding phase in the range of main values in the interferogram must be recovered to obtain the absolute phase, i.e. the phase needs to be unwrapped. The existence of residual points makes the unwrapping phase discontinuous and even an un-unwrapped 'island' appears. For this purpose, the invention introduces a small number of GPS control points, and takes the GPS control points as new unwrapping calculation points, thereby eliminating 'island' and making the unwrapping phase of the whole area continuous.
The invention provides two different GPS-assisted InSAR phase unwrapping methods according to the characteristics of error propagation of an unwrapping algorithm, namely a GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm and a Markov random field-based GPS-assisted InSAR unwrapping algorithm introducing residual points.
4. InSAR and GPS data fusion method based on improved Markov random field
For geological disasters such as landslide and earthquake, large-range deformation can be caused in a short time, and the repeated observation period of InSAR about one month is not enough, so that interpolation estimation in a time domain is needed. In addition, the GPS is based on observation of finite points, and the spatial density of the points cannot meet the requirement of monitoring, so that interpolation estimation is required in the spatial domain. The invention adopts an InSAR and GPS data fusion method based on an improved Markov random field, decomposes a one-dimensional deformation result of the InSAR visual line into three-dimensional deformation information through GPS constraint, and obtains a three-dimensional deformation field with high spatial resolution; then, performing kriging interpolation and encryption on the three-dimensional deformation field with high spatial resolution on high time frequency GPS data (the acquired ground settlement amount) in a time domain, so that a quasi-GPS result is encrypted into a quasi-GPS time sequence in the time domain to obtain the three-dimensional deformation field with high time resolution; and finally, estimating all points in the grid by using Kalman filtering to obtain a surface three-dimensional deformation result with high space-time resolution, and monitoring the surface deformation with high space-time resolution.
The invention relates to a surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion, which comprises the following steps:
1. laying GPS points
The distributed GPS points are used as control points for converting GPS and InSAR data coordinates and correcting InSAR errors, and also used as monitoring points for monitoring key parts of a measuring area. Thus, the deployed GPS points are divided into two categories: the first type is a stable characteristic point (called as a first type GPS point) in a monitoring area (the stability of the stable characteristic point is not influenced by deformation), such as the intersection of a river and a road, and the type of point is mainly used for unifying GPS data and InSAR data under the same reference coordinate system, correcting SAR satellite orbit errors and inverting atmospheric delay; the second type is the key monitoring part of the monitoring area (called the second type GPS point).
When the first type of GPS points are distributed, attention needs to be paid to the distribution of the points, and the points are uniformly distributed in a measuring area as much as possible.
2. Collecting and recording GPS signals and collecting SAR images
Because the orbit error and the atmospheric delay error of InSAR data are corrected by using GPS data, the acquired GPS signal is synchronized with the SAR image acquisition time, and the acquisition time period of the GPS information can be preferably extended forward and backward by about 10 minutes.
3. Three-dimensional coordinate information of two types of GPS points is solved according to GPS signals
The two types of GPS points are separately resolved, and the three-dimensional coordinate information of the resolved first type of GPS points comprises: the plane position information and the elevation information, and the three-dimensional coordinate information of the second type of GPS points comprises: planar displacement and vertical sedimentation.
4. Speckle noise suppression for SAR images
The speckle noise of the SAR image is effectively inhibited, and the method is beneficial to improving the accuracy of the InSAR monitoring result and ensuring the reliability of the monitoring result. The present invention employs a modified Lee algorithm (which is an existing algorithm).
5. Searching the characteristic points on the SAR image, and establishing the coordinate conversion relation between the searched characteristic points and the first type GPS points
In order to fuse the data of InSAR and GPS, the data of InSAR and GPS must be unified under the same reference coordinate system. And searching characteristic ground objects on the InSAR image, and measuring the ground coordinates of the characteristic ground objects by using the GPS. Hypothesis characterizationThe image coordinate of the point target is (I)row，Icol) With ground coordinates of (G)lat，Glan). Under a unified projection relationship, the transformation of the two sets of coordinate systems can be expressed as:
or
L＝BX (2)
Wherein, L is a GPS coordinate matrix of the characteristic point target, B is a corresponding image coordinate matrix, and x is a conversion matrix of 3 multiplied by 2 order. Thus, at least 3 feature point data are required to determine 6 parameters in the transformation matrix. When there are more feature point data, 6 transformation parameters can be obtained by least square adjustment. Assuming an observation error of
Then the observation equation is:
L＝BX+□ (3)
the estimated value of the transformation matrix is:
wherein,，D□is that
The variance matrix of (2). Finally, the atmosphere obtained by the GPS is extended by utilizing a conversion matrixAnd data translation and rotation such as late correction and track correction are converted into data under an InSAR image coordinate system.
6. Correction of orbital errors in SAR satellites using GPS data
A method combining a mapping equation method and least square estimation is adopted, and specifically the method comprises the following steps:
in the correction of the orbit error, the SAR image after the noise suppression is used as input data, a stable characteristic feature on the InSAR image is searched, and the ground coordinates thereof are measured by the GPS. According to the conversion relation between the radar coordinates (namely the slant range and the visual angle from the characteristic ground object to the InSAR image) of the characteristic ground object and the ground coordinates thereof, a terrain mapping equation is constructed, namely
Wherein X, Y and Z are ground coordinates of the characteristic points, a1，a2，…，c3For SAR satellite attitude angle
Elements of the constructed rotation matrix, λ0Is a scale factor, r is the slant distance from the feature to the SAR image, theta is the view angle, Xs，Ys，ZsIs the orbital position of the SAR satellite.
The orbit parameters and baseline parameters of the SAR satellite are then expressed as a linear function of time, namely:
Xs＝Xs0+Xs(t-t0)
. (6)
B＝B0+B(t-t0)
and (3) linearizing the topographic mapping equation to obtain an observation error equation:
V＝AX-L (7)
in the formula,
V＝[VX，VY，VZ]T
L＝[LX，LY，LZ]
wherein, Δ Xs0，ΔYs0，ΔZs0And
are each t0Track position and attitude correction number at time;
and
respectively corresponding rate of change correction; delta lambda0，ΔB0Are each t0Scale and baseline correction at time;
respectively corresponding rate of change correction; delta phi0Is a phase constant correction number; a is11，a12，…a317For partial derivatives of unknown parameters；VX，VY，VZThe difference between the known values of X, Y, Z and their approximate values, respectively.
Since the error equation contains 17 unknowns, at least more than 6 first-class GPS points are required to solve the unknown parameters. The method adopts a least square method to estimate 17 unknown parameters in an observation error equation, thereby obtaining the correction value of the SAR satellite orbit data error.
7. Forming InSAR differential interferograms
Firstly, carrying out fine registration on a main SAR image and a sub SAR image, wherein the registration precision is 1/8 pixels; secondly, resampling the auxiliary image into a coordinate system of the main image; finally, complex conjugate multiplication is carried out on the main and auxiliary images after registration to generate an interference pattern; then eliminating the flat land effect phase; and finally, carrying out differential interference processing by adopting a two-rail method or a three-rail method to obtain a differential interference pattern.
8. Utilizing the first type GPS points to invert the atmospheric delay and unify the atmospheric delay into the coordinate system of InSAR
And inverting the atmospheric delay by using the GPS, firstly obtaining troposphere estimation of the GPS by adopting a random process method, secondly estimating atmospheric delay correction by adopting an inter-station time domain double difference method, and unifying the atmospheric delay correction into a coordinate system of the InSAR. Wherein, the time domain double difference method between stations is as follows:
if point A is a GPS point on the SAR image and point B is another GPS point on the same SAR image, the tropospheric delay of A, B points on epoch j obtained by GPS solution is respectively marked as DA j、DB jThe difference between the tropospheric delays between the stations is
By taking A as a reference station, the equation (8) can be used for deducing the single difference between tropospheric delay stations of other GPS points.
Assuming that the two sites A, B, the two single differences established by equation (8) at the two epochs j (main image acquisition time) and k (sub-image acquisition time) are:
from the two single differences of equation (9), the epoch double difference can be found as:
as can be seen from equation (10), the estimated atmospheric correction double difference has two forms: one is the time domain double difference between stations, i.e. the difference between stations is first and then the difference between epochs, and the other is the time domain double difference between stations, i.e. the difference between epochs is first and then the difference between stations. Because the inter-station difference is only related to a single SAR image, the inter-epoch difference can be formed by combining with other SAR images with the inter-station difference more freely and flexibly, and therefore, the inter-station time domain double difference which is more beneficial to practical application is adopted.
9. Atmospheric delay correction InSAR differential interferogram atmospheric delay error by utilizing GPS inversion
Considering that the atmospheric delay has strong correlation with the elevation, according to the theory of geographic spatial correlation, introducing an elevation influence factor, constructing a spatial different correlation model, and providing an improved inverse distance weighted interpolation method for fusing GPS elevation information to encrypt the model. And encrypting the atmospheric delay inverted by the GPS, and then eliminating the encrypted atmospheric delay from the InSAR interferometric phase diagram. The basic principle of the improved inverse distance weighted interpolation method for fusing GPS elevation information is as follows:
inverse distance weighted interpolation is a commonly used spatial interpolation method that assumes that the degree of similarity of two things increases as the distance between them decreases. Therefore, the Inverse Distance Weighting (IDW) interpolation performs weighted averaging using the distance between the estimated value point and the observation point as a weight, and the closer the estimated value point is, the higher the weight is given to the observation point. Inverse Distance Weighted (IDW) interpolation can be expressed as:
Z(x0，y0，z0) Is a point (x)0，y0，z0) An estimated value of (c), Z (x)i，yi，zi) Is a point (x)i，yi，zi) Where n is the number of observations used for estimation, λiIs a point (x)i，yi，zi) The weight corresponding to the observed value, di0(i 1 … N) is the distance from observation point i to the evaluation point, p (p > 0) is the power of the distance, and
in practice, p is usually 2.
As can be seen from equation (11), the inverse distance weighted interpolation method considers the influence of the horizontal distance and the elevation difference between the evaluation point and the observation point on the evaluation point to be equivalent. However, this is not the case, for example, in the interpolation problem of atmospheric delay, the horizontal distance and the height difference between the estimation point and the observation point have different effects on the estimation point. In order to reflect the difference of the horizontal distance and the elevation difference between an estimated value point and an observed point on the estimation value, according to the geographic space correlation theory, an elevation influence factor alpha is introduced, a space autocorrelation model is constructed, and an improved inverse distance weighting interpolation method fusing elevation information is provided, wherein the specific interpolation formula is as follows:
wherein, Z (x)0，y0，z0) Is a point (x)0，y0，z0) An estimated value of (c), Z (x)i，yi，zi) Is a point (x)i，yi，zi) Where n is the number of observations used for estimation, λiIs a point (x)i，yi，zi) The weight corresponding to the observed value, hiIs a point (x)i，yi，zi) Elevation of (c), h0Is a point (x)0，y0，z0) Elevation of (b)0、b1For the fitting coefficient to be found, di0 sFor horizontal distance from observation point i to evaluation point, di0 hIs the difference in elevation, k, of the observation point i and the estimation point1、k2The slope of a linear fit curve of the difference in atmospheric delay as a function of horizontal distance and elevation difference, respectively, abs (-) represents the absolute value.
The improved inverse distance weighted interpolation method for integrating elevation information is based on the theory of geographic spatial correlation, further assumes the spatial distribution of various geographic elements has intrinsic correlation on the basis of the assumption of spatial autocorrelation, and infers a spatial autocorrelation model of the spatial distribution of other elements according to the spatial distribution of one element. Therefore, the improved inverse distance weighted interpolation method for fusing elevation information is more beneficial to improving the interpolation precision when the elevation of the station to be estimated has larger change than that of the surrounding sampling station.
10. Unwrapping the InSAR differential interferogram after atmospheric correction
The invention designs two GPS-assisted InSAR unwrapping algorithms: the system comprises a GPS assisted InSAR Goldstein branch tangent unwrapping algorithm and a Markov random field-based GPS assisted InSAR unwrapping algorithm introducing residual points. The two methods have respective advantages and disadvantages, the GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm adopts an accurate integral algorithm, the unwrapping precision is high, the calculation speed is high, the partial unwrapping is only realized, and the unwrapping range and precision are greatly influenced by the noise level; the Markov random field GPS-assisted InSAR unwrapping algorithm introduced with the residual points is a global unwrapping algorithm and has strong noise suppression capability, but because the simulated annealing optimization search algorithm is adopted, the unwrapping speed is slow, and the possibility of falling into local optimization inevitably occurs, so that the unwrapping error is caused. Therefore, the user can select different unwrapping algorithms according to different needs. The basic principles of the two unwrapping algorithms are as follows:
(1) GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm
When phase unwrapping is performed using the Goldstein branch-and-cut method, such a situation occurs: a small area surrounded by branch cuts is completely isolated from other areas of the interferogram, such a small area is visually referred to as an "unwrapped island". The formation reason of the 'disentangled islands' is many, such as small islands in water, surrounded open spaces and the like. For the 'unwinding island', the traditional processing method comprises the following steps: (a) when the island area is small, the island area can be abandoned to be unwound during unwinding, only other areas can be unwound, and after the final products (such as DEM or deformation) of other areas are obtained, the island area is interpolated by using a proper interpolation algorithm, so that the final products of the island area are obtained. The disadvantage of this method is that the phase information on the island is ignored and the interpolation process introduces large errors. (b) And when the island area is large, reselecting a de-winding calculation point in the island area to perform independent phase de-winding on the island area. Although the method uses the phase information on the island, the island area and other areas of the interferogram are not unified in the unwinding reference due to the adoption of the new unwinding calculation point, so that the interpretation of the final product is difficult. Aiming at the problems existing in the traditional processing method, the InSAR deformation measurement precision is usually centimeter level, and the GPS deformation measurement precision can reach millimeter level.
The principle of the GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm is as follows:
introducing discrete GPS points, selecting one GPS point as a starting point, solving the relative elevations or deformations of other GPS points relative to the starting point, and converting the relative elevations or deformations into unwrapped phase values by using the following formula:
wherein phi istopo uAnd phidefo uRespectively, unwrapped terrain interference phase and unwrapped deformation interference phase, and deltah and deltar are respectively the relative elevation and deformation relative to the total unwrapped calculation point.
And then, the island points are used as starting points to be unwound again until all the island points are solved.
The method is a local unwrapping algorithm and has the advantages of high unwrapping speed, high accuracy, unified reference and capability of evaluating unwrapping accuracy.
(2) Markov random field-based GPS-assisted InSAR unwrapping algorithm introducing residual points
The identification and connection of residual points are the biggest characteristics of the Goldstein branch tangent unwrapping algorithm, and can effectively isolate error points and avoid the global propagation of errors. However, island can be generated by unwinding the Goldstein branch tangent, and the pixel points on the branch tangent cannot be unwound, so that the effective range of unwinding is reduced. The GPS-assisted unwrapping algorithm based on the Markov random field proposed by Gudmundsson does not consider the influence of residual points, and can cause the global propagation of errors. Therefore, the invention introduces the residual error point identification technology into the GPS auxiliary unwrapping algorithm based on the Markov random field, provides the GPS auxiliary InSAR unwrapping algorithm based on the Markov random field with the introduced residual error points, and realizes the optimal integration of the unwrapping precision and the unwrapping range. The method converts the InSAR unwrapping process into an energy solving function U (K is K, Y is I)w) The course of the minimum, namely:
obtaining a global optimal solution at a high rate and probability by adopting a repeated rapid annealing strategy, namely, giving a large temperature decay rate cool and a termination temperature TsWhen T is less than TsWhen T is equal to T0And repeating the rapid annealing until a global optimal solution is reached. Initial integer number matrix k used in iteration, initial unwrapping interferogram IuUnwrapped interferogram I'uAnd the proportional coefficient r is calculated as follows:
where k is the initial integer matrix, round (□) is the operator, Iu GPSFor unwrapped interferograms inverted from GPS interpolation results, IwFor the winding interferogram, λ is the wavelength, IuIs an initial unwrapped interferogram, I'uFor unwrapping interferograms, r is the scale factor.
Different from a method that the Goldstein branch tangent unwrapping algorithm firstly identifies residual points and then connects branch tangents, the residual point identification and marking strategy adopted by the method is as follows: and calculating a 2 multiplied by 2 closed path unit accumulated value S, and marking all 4 pixel points as residual points if S is not zero, thereby completing the identification and marking of the residual points. The advantage of such marking the residual points is: (1) the error of marking the position of the residual error point, which is possibly caused by artificially appointing the position of the residual error point (the upper left corner point) in the traditional residual error point identification method, is avoided, so that the propagation of the error is more effectively inhibited; (2) complex branch tangent connection is not performed any more, so that the calculation efficiency is improved.
11. Converting InSAR unwrapping phase into deformation information
Wherein, □ RdFor the amount of deformation of the line of sight, phidTo the unwrapped phase.
12. InSAR deformation information geocoding
The step is a process of converting the result obtained by the previous processing under the radar coordinate system into a geographic coordinate system, and is carried out by adopting a range-Doppler (R-D) positioning model. In the R-D positioning model, the position of any pixel on the SAR image is determined by three equations: a distance equation, a doppler centroid plane equation, and an earth model equation.
Distance equation: <math>
<mrow>
<mi>R</mi>
<mo>=</mo>
<msup>
<mrow>
<mo>[</mo>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>R</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>R</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>·</mo>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>R</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>R</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>]</mo>
</mrow>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
</msup>
</mrow>
</math>
doppler equation: <math>
<mrow>
<msub>
<mi>f</mi>
<mi>D</mi>
</msub>
<mo>=</mo>
<mfrac>
<mn>2</mn>
<mi>λR</mi>
</mfrac>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>V</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>V</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>·</mo>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>R</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>R</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
</mrow>
</math>
an earth model equation:
in the formula, R is the distance between the SAR sensor and a ground point,
position vectors for satellite and ground points, respectively, ". is a vector point product, fDFor doppler shift, λ is the radar wavelength,
in the form of a vector of the velocity of the satellite,
for the ground point velocity vector, when using the WGS-84 coordinate system,
is zero, xt，ytAnd ztAs a ground point position vector
The projection on the x, y and z axes, a is the earth's major half axis, a is 6378.137km, b is the earth's minor half axis, b is 6356.752km, and h is the average normal height. The geocoding process of the measurement result can be completed by solving the three equations.
13. A Global Positioning System (GPS) and interferometric synthetic aperture radar (InSAR) fusion method based on an improved Markov random field is adopted to fuse deformation information of InSAR and GPS to obtain a ground surface three-dimensional deformation monitoring result with high spatial and temporal resolution
Firstly, decomposing a one-dimensional deformation result of the InSAR visual line into three-dimensional deformation information through GPS constraint to obtain a three-dimensional deformation field with high spatial resolution. And obtaining an optimal result by adopting a decomposition optimization method, wherein a cross validation weighting method is adopted for determining the weight of the GPS constraint condition. The specific principle is as follows:
constructing a relation between a one-dimensional deformation result and a three-dimensional deformation result of the InSAR visual line:
dlos＝sinθcos(αh-3π/2)dn+sinθsin(αh-3π/2)de-cosθdu
(17)
＝[un，ue，uu][dn，de，du]T
wherein d islosFor InSAR line-of-sight distortion measurement, [ u [ u ] ]n，ue，uu]Is a line of sight unit vector, [ u ]n，ue，uu]＝[0.34-0.095，0.935]，dn、deAnd duThe optimal value is the northeast direction, theta is the radar view angle, and the value of theta changes along with the distance from the pixel point to the satellite subsatellite point; alpha is alphahThe angle between the flight direction of the satellite and the north direction can be read from SAR data.
Determining a total energy function for modeling calculation by using a Markov random field:
wherein,
in the formula (18), U represents the total energy, dlos iFor InSAR line-of-sight distortion measurement, [ u [ u ] ]n i，ue i，uu i]Is a line of sight unit vector, dn i、de iAnd du iOptimized value for northeast heaven direction, σins iFor InSAR observation error, dkrig ni、dkrig eiAnd dkrig uiFor northeast direction kriging interpolation, sigmakrig ni、σkrig eiAnd σkrig uiThe kriging interpolation standard deviation in the northeast direction.
Since 4 × M non-negative terms are included in equation (18), equation (18) reaches the global minimum when all non-negative terms reach the minimum. Thus, pair dn i、de iAnd du iCalculating the partial derivative and making the partial derivative zero to obtain a calculation formula of the analytical optimization method, namely:
the weight determination method in the formula (18) has two problems of how to evaluate InSAR measurement precision and Kriging interpolation precision, and for effectively solving the two problems, the invention determines the weight by adopting a cross validation weighting method, and the specific process is as follows:
●, registering the GPS observation value and the InSAR observation value, converting the GPS three-dimensional observation value to the InSAR sight direction, and estimating the measurement variance of the InSAR pixel corresponding to the GPS observation station according to the formula (20):
wherein n is the number of GPS observation stations, dGPS，i losConverting the three-dimensional observation value of the ith GPS observation station into the InSAR sight line direction result, dins，i losAnd obtaining an InSAR ith pixel sight line direction observation value.
● calculate the measurement variance for each pixel point of InSAR:
in the formula, σins，i 2And gammaiMeasure variance estimates and coherence coefficients, E (γ), for InSAR ith pixel points, respectivelyGPS) And the average value of the coherent coefficients of the InSAR pixels corresponding to the GPS observation station.
● dividing the GPS observation values into A, B groups, performing ordinary Kriging interpolation on the group A, and solving the variance between the B groups of GPS three-dimensional observation values and the corresponding interpolation:
wherein n' is the number of GPS observation stations in the B group, dGPS，i n、dGPS，i eAnd dGPS，i uIs the observation value of the ith GPS observation station of the B group in the northeast direction, dkrig，i n、dkrig，i eAnd dkrig，i uAnd interpolating for the ith GPS observation station Crigger of the B group in the northeast direction.
● calculate a kriging interpolation variance correction factor:
in the formula, E [ (sigma)krig，n 2)B]、E[(σkrig，e 2)B]And E [ (sigma)krig，u 2)B]And respectively the average value of the Kriging interpolation variances in the northeast direction of the GPS observation stations in the B groups.
● correction of kriging variance:
wherein (sigma)krig，i n)2、(σkrig，i e)2And (σ)krig，i u)2And the kriging interpolation variance in the northeast direction of the ith pixel point is obtained.
● determining the weight:
and then, performing Kriging interpolation and encryption on the three-dimensional deformation field with high spatial resolution on the high-time frequency GPS data (the acquired ground settlement amount thereof) in a time domain, so that the quasi-GPS result is encrypted into a quasi-GPS time sequence in the time domain to obtain the three-dimensional deformation field with high time resolution.
And finally, estimating all points in the grid by using Kalman filtering to obtain a surface three-dimensional deformation result with high space-time resolution.
The invention discloses a surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion, which breaks through the limitation of single technology application by integrating the advantages of InSAR and GPS, provides a new technology with high space-time resolution and high precision for surface monitoring, and is suitable for large-range surface deformation monitoring.
Drawings
FIG. 1: a flow chart of data processing of the method of the invention;
FIG. 2: a flow chart of correcting track errors of InSAR data by the GPS;
FIG. 3: a schematic diagram of correcting InSAR atmospheric delay errors by GPS inversion atmospheric delay;
FIG. 4: the schematic diagram of the unwinding island is generated;
FIG. 5: considering residual points, and based on a Markov random field GPS auxiliary InSAR phase unwrapping data processing flow chart;
Detailed Description
The present invention will be described in further detail with reference to examples.
Example (b):
1. laying GPS points
The method lays two types of GPS points: the first type is a stable characteristic point in a monitoring area, such as a river and road intersection, and the point is mainly used for unifying GPS data and InSAR data to the same reference coordinate system, correcting SAR satellite orbit errors and inverting atmospheric delay; the second type is the key monitoring part of the monitoring area.
When the first type of GPS points are distributed, attention needs to be paid to the distribution of the points, and the points are uniformly distributed in a measuring area as much as possible.
2. Collecting and recording GPS signals and collecting SAR images
And acquiring a GPS signal and an SAR image which are synchronous in time, wherein the acquired GPS information is forwards and backwards delayed for about 10 minutes relative to the acquisition time of the SAR image.
3. Three-dimensional coordinate information of two types of GPS points is solved according to GPS signals
The two types of GPS are separately solved, and the three-dimensional coordinate information of the first type of GPS point comprises the following steps: the plane position information and the elevation information, and the three-dimensional coordinate information of the second type of GPS points comprises: planar displacement and vertical sedimentation.
4. Speckle noise suppression for SAR images
The method adopts an improved Lee algorithm to inhibit speckle noise of the SAR image.
5. Searching the characteristic points on the SAR image, and establishing the coordinate conversion relation between the searched characteristic points and the first type GPS points
Firstly, searching more than 3 InSAR images for feature point target ground objects, and measuring the ground coordinates by using a GPS. Then, the estimated value of the transformation matrix between the two sets of coordinate systems is calculated by the following formula:
wherein, L is the ground coordinate matrix of the characteristic point target, B is the corresponding image coordinate matrix,D□is an observation error
The variance matrix of (2).
And finally, translating and rotating data obtained by the GPS, such as atmospheric delay correction, orbit correction and the like, into data under an InSAR image coordinate system by using a conversion matrix.
6. Correction of orbital errors in SAR satellites using GPS data
In the correction of the orbit error, first, the SAR image after the noise suppression is used as input data, and the target ground object of the feature point on the 6 or more InSAR images is searched and the ground coordinates thereof are measured by the GPS.
Secondly, a topographic mapping equation is constructed, i.e.
Wherein X, Y and Z are ground coordinates of the feature point target ground object, a1，a2，…，c3For SAR satellite attitude angle
The formed rotation matrix elements include the radar wavelength lambda, the light speed c, the slant distance from the feature point target ground object to the SAR image r, the view angle theta, and the Xs，Ys，ZsPhi is the unwrapping phase for the orbital position of the SAR satellite.
The orbit parameters and baseline parameters of the SAR satellite are then expressed as a linear function of time, namely:
Xs＝Xs0+Xs(t-t0)
B＝B0+B(t-t0)
and (3) linearizing the topographic mapping equation to obtain an observation error equation:
V＝AX-L
in the formula,
V＝[VX，VY，VZ]T
L＝[LX，LY，LZ]
wherein, Δ Xs0，ΔYs0，ΔZs0Andare each t0Track position and attitude correction number at time;
and
respectively corresponding rate of change correction; delta lambda0，ΔB0Are each t0Scale and baseline correction at time;
respectively corresponding rate of change correction; delta phi0Is a phase constant correction number; a is11，a12，…a317Partial derivatives of each unknown parameter; vx，VY，VZThe difference between the known values of X, Y, Z and their approximate values, respectively.
And (3) estimating 17 unknown parameters in the observation error equation by adopting a least square method, thereby obtaining the correction value of the SAR satellite orbit data error.
7. Forming InSAR differential interferograms
Firstly, carrying out fine registration on a main SAR image and a sub SAR image, wherein the registration precision is 1/8 pixels; secondly, resampling the auxiliary image into a coordinate system of the main image; finally, complex conjugate multiplication is carried out on the main and auxiliary images after registration to generate an interference pattern; then eliminating the flat land effect phase; and finally, carrying out differential interference processing by adopting a two-rail method or a three-rail method to obtain a differential interference pattern.
8. Utilizing the first type GPS points to invert the atmospheric delay and unify the atmospheric delay into the coordinate system of InSAR
Firstly, acquiring troposphere estimation of a GPS (global positioning system) by adopting a random process method, then estimating atmospheric delay correction by adopting an inter-station time domain double-difference method, and unifying the atmospheric delay correction into a coordinate system of an InSAR (interferometric synthetic aperture radar). Wherein, the time domain double difference method between stations is as follows:
wherein, the point A is a GPS reference point on the SAR image, the point B is another GPS point on the same SAR image, j and k are respectively the acquisition time of the main and auxiliary images, DA jAnd AB jThe tropospheric delays at epoch j for the two points A, B resolved by the GPS are shown.
9. Atmospheric delay correction InSAR differential interferogram atmospheric delay error by utilizing GPS inversion
Firstly, an improved inverse distance weighted interpolation method which integrates GPS elevation information is adopted to encrypt the atmospheric delay of GPS inversion, and then the encrypted atmospheric delay is eliminated from an InSAR interferometric phase diagram. The improved inverse distance weighted interpolation method for fusing the GPS elevation information comprises the following steps:
wherein, Z (x)0，y0，z0) Is a point (x)0，y0，z0) An estimated value of (c), Z (x)i，yi，zi) Is a point (x)i，yi，zi) Where n is the number of observations used for estimation, λiIs a point (x)i，yi，zi) The weight corresponding to the observed value, hiIs a point (x)i，yi，zi) Elevation of (b)0、b1For the fitting coefficient to be found, di0 sFor horizontal distance from observation point i to evaluation point, di0 hIs the difference in elevation, k, of the observation point i and the estimation point1、k2The slope of a linear fit curve of the difference in atmospheric delay as a function of horizontal distance and elevation difference, respectively, abs (-) represents the absolute value.
10. Unwrapping the InSAR differential interferogram after atmospheric correction
The method designs two GPS-assisted InSAR unwrapping algorithms: the system comprises a GPS assisted InSAR Goldstein branch tangent unwrapping algorithm and a Markov random field-based GPS assisted InSAR unwrapping algorithm introducing residual points.
The two methods have respective advantages and disadvantages, and the GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm is a local unwrapping algorithm, has the advantages of high unwrapping speed, high precision, uniform reference and capability of evaluating unwrapping precision, but has larger influence on unwrapping range and precision by noise level; the Markov random field GPS-assisted InSAR unwrapping algorithm introduced with the residual points is a global unwrapping algorithm, and is strong in noise suppression capability and low in unwrapping speed. The user can select different unwrapping algorithms according to different needs.
1) GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm
In the GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm, the data processing steps are as follows:
(1) the GPS observers are registered into the wrapped interferogram.
(2) And solving a main value of the winding phase gradient field, identifying residual points and connecting branch tangents.
(3) A total unwrapping calculation point, such as point 1, is selected from GPS observation stations located outside the island, and the relative elevations or deformations of other GPS points with respect to point 1 are determined and converted to unwrapped phase values using the following equation:
in the above formula, phitopo uAnd phidefo uRespectively, unwrapped terrain interference phase and unwrapped deformation interference phase, and deltah and deltar are respectively the relative elevation and deformation relative to the total unwrapped calculation point.
(4) And (4) unwrapping the interferogram by using a Goldstein branch tangent method by taking the point 1 as an unwrapping calculation point until all unwrappable points are unwrapped. Taking the case of fig. 4 as an example, in this process all points outside the island will be successfully unwound, while points inside the island cannot be unwound.
(5) The interferogram is scanned to find the GPS point which is not unwound yet, namely point 2, and a new round of phase unwinding is carried out by taking the point as a new unwinding calculation point.
(6) And (5) repeating the step until the whole interference pattern is wound.
(7) And finding out GPS points which are unwound but have not been calculated, namely the points 3 and 4, and carrying out statistical analysis on the known value of the unwinding phase obtained by inversion in the step 3 and the estimation of the unwinding phase obtained in the step 4-5, thereby realizing the precision evaluation of the unwinding algorithm.
(8) And (4) giving the known value of the unwrapping phase inverted by the GPS to a corresponding point on the unwrapping interferogram, and ending the whole unwrapping process.
2) Markov random field-based GPS-assisted InSAR unwrapping algorithm introducing residual points
The method converts the InSAR unwrapping process into an energy solving function U (K is K, Y is I)w) Minimum valueThe process of (a), namely:
obtaining a global optimal solution at a high rate and probability by adopting a repeated rapid annealing strategy, namely, giving a large temperature decay rate cool and a termination temperature TsWhen T is less than TsWhen T is equal to T0And repeating the rapid annealing until a global optimal solution is reached. Initial integer number matrix k used in iteration, initial unwrapping interferogram IuUnwrapped interferogram I'uAnd the proportional coefficient r is calculated as follows:
where k is the initial integer matrix, round (□) is the operator, Iu GPSFor unwrapped interferograms inverted from GPS interpolation results, IwFor the winding interferogram, λ is the wavelength, IuIs an initial unwrapped interferogram, I'uFor unwrapping interferograms, r is the scale factor.
Wherein the residual point identification and marking strategy is: and calculating a 2 multiplied by 2 closed path unit accumulated value S, and marking all 4 pixel points as residual points if S is not zero, thereby completing the identification and marking of the residual points. And after the residual error point marking is finished, the non-residual error points are unwound by using a GPS auxiliary InSAR unwinding algorithm based on a Markov random field. After the step of unwrapping, marking all unwrapped points as fixed domains, and performing an InSAR unwrapping algorithm assisted by a Markov random field GPS only by adopting smooth constraint once again, thereby completing the unwrapping of the whole interferogram. Therefore, the Markov random field GPS-assisted InSAR unwrapping algorithm considering residual points firstly unwrapps the high-quality primitive points, and then further unwrapps the low-quality points by taking the unwrapping values of the high-quality primitive points as constraints, so that error propagation is restrained to the maximum extent. The Markov random field GPS-assisted InSAR unwrapping algorithm considering the residual points realizes better integration of the unwrapping precision and range. The data processing flow of the Markov random field GPS assisted InSAR unwrapping algorithm considering the residual points is shown in FIG. 5.
11. Converting InSAR unwrapping phase into deformation information
Wherein, □ RdFor the amount of deformation of the line of sight, phidTo the unwrapped phase.
12. InSAR deformation information geocoding
The method adopts a distance Doppler (R-D) positioning model, and converts the result of the radar coordinate system obtained by the previous processing into a process of the geographic coordinate system. In the R-D positioning model, the position of any pixel on the SAR image is determined by three equations: a distance equation, a doppler centroid plane equation, and an earth model equation.
Distance equation: <math>
<mrow>
<mi>R</mi>
<mo>=</mo>
<msup>
<mrow>
<mo>[</mo>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>R</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>R</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>·</mo>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>R</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>R</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>]</mo>
</mrow>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
</msup>
</mrow>
</math>
doppler equation: <math>
<mrow>
<msub>
<mi>f</mi>
<mi>D</mi>
</msub>
<mo>=</mo>
<mfrac>
<mn>2</mn>
<mi>λR</mi>
</mfrac>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>V</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>V</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>·</mo>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>R</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>R</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
</mrow>
</math>
an earth model equation:
in the formula, R is the distance between the SAR sensor and a ground point,
position vectors for satellite and ground points, respectively, ". is a vector point product, fDFor doppler shift, λ is the radar wavelength,
in the form of a vector of the velocity of the satellite,
for the ground point velocity vector, when using the WGS-84 coordinate system,
is zero, xt，ytAnd ztAs a ground point position vector
Projection onto x, y and z axesA is a long half shaft, b is a short half shaft, h is an average normal height, a is 6378.137km, and b is 6356.752 km. The geocoding process of the measurement result can be completed by solving the three equations.
13. A Global Positioning System (GPS) and interferometric synthetic aperture radar (InSAR) fusion method based on an improved Markov random field is adopted to fuse deformation information of InSAR and GPS to obtain a ground surface three-dimensional deformation monitoring result with high spatial and temporal resolution
Firstly, decomposing a one-dimensional deformation result of the InSAR visual line into three-dimensional deformation information through GPS constraint. The specific algorithm is as follows:
(1) establishing the relation between the one-dimensional deformation result and the three-dimensional deformation result of the visual line
dlos＝sinθcos(αh-3π/2)dn+sinθsin(αh-3π/2)de-cosθdu
(28)
＝[un，ue，uu][dn，de，du]T
Wherein d islosFor InSAR line-of-sight distortion measurement, [ u [ u ] ]n，ue，uu]Is a line of sight unit vector, [ u ]n，ue，uu]＝[0.34，-0.095，0.935]，dn、deAnd duThe optimal value is the northeast direction, theta is the radar view angle, and the value of theta changes along with the distance from the pixel point to the satellite subsatellite point; alpha is alphahThe angle between the flight direction of the satellite and the north direction can be read from SAR data.
(2) Determining a total energy function
Wherein,
wherein, Wins iIs InSAR constraint weight, Wkrig ni，Wkrig ei，Wkrig uiIs GPS constraint weight, dlos iFor InSAR line-of-sight distortion measurement, [ u [ u ] ]n i，ue i，uu i]Is a line of sight unit vector, dn i、de iAnd du iOptimized value for northeast heaven direction, σins iFor InSAR observation error, dkrig ni、dkrig eiAnd dkrig uiFor the kriging interpolation in the northeast direction,σkrig ni、σkrig eiand σkrig uiThe kriging interpolation standard deviation in the northeast direction.
(3) Determining a constraint weight by adopting a cross validation weighting method, which comprises the following specific steps:
●, the GPS and InSAR observations are registered, and the GPS three-dimensional observations are converted to InSAR sight direction.
● estimating the measurement variance of InSAR image elements corresponding to the GPS observation station according to the following formula:
wherein n is the number of GPS observation stations, dGPS，i losConverting the three-dimensional observation value of the ith GPS observation station into the InSAR sight line direction result, dins，i losAnd obtaining an InSAR ith pixel sight line direction observation value.
The measurement variance of ● InSAR per pixel point can be calculated by:
in the formula, σins，i 2And gammaiMeasure variance estimates and coherence coefficients, E (γ), for InSAR ith pixel points, respectivelyGPS) And the average value of the coherent coefficients of the InSAR pixels corresponding to the GPS observation station.
● the GPS observations were divided into A, B groups and the A group was interpolated by ordinary kriging.
●, calculating the variance between the B groups of GPS three-dimensional observed values and the corresponding interpolation:
wherein n' is the number of GPS observation stations in the B group, dGPS，i n、dGPS，i eAnd dGPS，i uIs the ith G in B group in northeast directionPS Observation station observed value, dkrig，i n、dkrig，i eAnd dkrig，i uAnd interpolating for the ith GPS observation station Crigger of the B group in the northeast direction.
● calculate a kriging interpolation variance correction factor:
in the formula, E [ (sigma)krig，n 2)B]、E[(σkrig，e 2)B]And E [ (sigma)krig，u 2)B]And respectively the average value of the Kriging interpolation variances in the northeast direction of the GPS observation stations in the B groups.
● correction of kriging variance:
wherein (sigma)krig，i n)2、(σkrig，i e)2And (σ)krig，i u)2And the kriging interpolation variance in the northeast direction of the ith pixel point is obtained.
● weight:
(4) obtaining maximum posterior estimation of three-dimensional deformation field by direct analytic method
And then, performing Kriging interpolation and encryption on the three-dimensional deformation field with high spatial resolution on the high-time frequency GPS data (the acquired ground settlement amount thereof) in a time domain, so that the quasi-GPS result is encrypted into a quasi-GPS time sequence in the time domain to obtain the three-dimensional deformation field with high time resolution.
And finally, estimating all points in the grid by using Kalman filtering to obtain a surface three-dimensional deformation result with high space-time resolution.
Claims (2)
1. A surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion comprises the following steps:
step one, laying GPS points
There are two types of deployed GPS points: the first type is a stable characteristic point in a monitoring area and is used for unifying GPS data and InSAR data under the same reference coordinate system, correcting SAR satellite orbit errors and inverting atmospheric delay; the second type is monitoring points of a monitoring area;
step two, collecting and recording GPS signals and collecting SAR images
The acquired GPS signal is synchronized with the SAR image acquisition time;
step three, calculating three-dimensional coordinate information of two types of GPS points according to GPS signals
The three-dimensional coordinate information of the first type of GPS monitoring points comprises: the position information and the elevation information, and the three-dimensional coordinate information of the second type of GPS monitoring points comprises: position information and deformation information;
step four, speckle noise suppression of SAR image
Step five, searching the characteristic points on the SAR image, and establishing the coordinate conversion relation between the searched characteristic points and the first type GPS points
Step six, correcting the orbit error of the SAR satellite by using the GPS data
A method combining a mapping equation method and least square estimation is adopted;
step seven, forming an InSAR differential interference pattern
Step eight, inverting the atmospheric delay by utilizing the first type of GPS points and unifying the atmospheric delay into the coordinate system of the InSAR
Acquiring troposphere estimation of a GPS (global positioning system) by adopting a random process method, and estimating atmospheric delay correction by adopting an inter-station time domain double-difference method;
step nine, correcting the atmospheric delay error of the InSAR differential interferogram by utilizing the atmospheric delay of GPS inversion
Encrypting the atmospheric delay inverted by the GPS by adopting an improved inverse distance weighted interpolation method for fusing GPS elevation information, and then eliminating the encrypted atmospheric delay from the InSAR differential interference phase diagram;
step ten, unwrapping the InSAR differential interferogram after atmospheric correction
The InSAR differential interferogram after atmosphere correction is unwound by adopting a GPS-assisted InSAR Goldstein branch tangent unwinding algorithm or a Markov random field-based GPS-assisted InSAR unwinding algorithm introducing residual points;
eleven, converting the InSAR unwrapping phase into deformation information
Wherein, □ RdFor the amount of deformation of the line of sight, phidFor unwrapping phase, λ is the radar wavelength;
step twelve, InSAR deformation information geocoding
Distance equation: <math>
<mrow>
<mi>R</mi>
<mo>=</mo>
<msup>
<mrow>
<mo>[</mo>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>R</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>R</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>·</mo>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>R</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>R</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>]</mo>
</mrow>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
</msup>
</mrow>
</math>
doppler equation: <math>
<mrow>
<msub>
<mi>f</mi>
<mi>D</mi>
</msub>
<mo>=</mo>
<mfrac>
<mn>2</mn>
<mi>λR</mi>
</mfrac>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>V</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>V</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>·</mo>
<mrow>
<mo>(</mo>
<mover>
<msub>
<mi>R</mi>
<mi>s</mi>
</msub>
<mo>→</mo>
</mover>
<mo>-</mo>
<mover>
<msub>
<mi>R</mi>
<mi>t</mi>
</msub>
<mo>→</mo>
</mover>
<mo>)</mo>
</mrow>
</mrow>
</math>
an earth model equation:
in the formula, R is the distance between the SAR sensor and a ground point,
position vectors for satellite and ground points, respectively, ". is a vector point product, fDIn order to be the doppler shift frequency,
in the form of a vector of the velocity of the satellite,
is the ground point velocity vector, xt，ytAnd ztAs a ground point position vector
Projection on x, y and z axes, wherein a is a half axis of the earth's major axis, b is a half axis of the earth's minor axis, and h is an average normal height;
solving the three equations to obtain the InSAR deformation information geocode;
step thirteen, adopting a GPS and InSAR fusion method based on the improved Markov random field, fusing deformation information of the InSAR and the GPS to obtain a high-precision surface three-dimensional deformation monitoring result;
and finally, estimating all points in the grid by using Kalman filtering to obtain a surface three-dimensional deformation result with high space-time resolution.
2. The earth's surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion of claim 1, characterized in that: the GPS and InSAR fusion method based on the improved Markov random field in the step thirteen comprises the following steps:
firstly, decomposing a one-dimensional deformation result of the InSAR visual line into three-dimensional deformation information through GPS constraint; the specific algorithm is as follows:
(1) establishing the relation between the one-dimensional deformation result and the three-dimensional deformation result of the visual line
dlos＝sinθcos(αh-3π/2)dn+sinθsin(αh-3π/2)de-cosθdu
＝[un，ue，uu][dn，de，du]T
Wherein d islosFor InSAR line-of-sight distortion measurement, [ u [ u ] ]n，ue，uu]Is a line of sight unit vector, [ u ]n，ue，uu]＝[0.34，-0.095，0.935]，dn、deAnd duThe optimal value is the northeast direction, theta is the radar view angle, and the value of theta changes along with the distance from the pixel point to the satellite subsatellite point; alpha is alphahThe included angle between the satellite flight direction and the north direction is read out from SAR data;
(2) determining a total energy function
Wherein,
wherein,
for the InSAR constraint weight, the weight is,
in order to constrain the weight values for the GPS,
for InSAR line-of-sight distortion measurement,is a line-of-sight unit vector,
and
for the optimization of the values in the north east direction,
in order to observe the error of the InSAR,
and
for the kriging interpolation in the northeast direction,
and
the standard deviation of the Kriging interpolation in the northeast direction;
(3) determining a constraint weight by adopting a cross validation weighting method, which comprises the following specific steps:
registering the GPS and InSAR observation values, and converting the GPS three-dimensional observation value to the InSAR sight direction;
estimating the measurement variance of the InSAR image element corresponding to the GPS observation station according to the following formula:
wherein n is the number of GPS observation stations,
converting the three-dimensional observation value of the ith GPS observation station into the InSAR sight line direction result,
the observation value of the ith pixel sight line direction of the InSAR is obtained;
the measurement variance of each pixel point of InSAR is calculated by the following formula:
in the formula,
and gammaiMeasure variance estimates and coherence coefficients, E (γ), for InSAR ith pixel points, respectivelyGPS) The coherent coefficient mean value of the InSAR pixels corresponding to the GPS observation station;
dividing the GPS observed values into A, B two groups, and performing common Kriging interpolation on the group A;
solving the variance between the three-dimensional observation values of the group B GPS and the corresponding interpolation values:
wherein n' is the number of the GPS observation stations in the B group,
and
for the north east direction of the day B group of i GPS observatory observations,andinterpolating for the ith GPS observation station kriging in the B group in the northeast direction;
calculating a kriging interpolation variance correction factor:
in the formula,
and
respectively is a Kriging interpolation variance mean value in the northeast direction of the GPS observation station in the group B;
correcting the kriging variance:
in the formula,
andinterpolating variance for the pixel point of the ith pixel point in the direction of northeast sky by kriging; and (3) weighting:
(4) obtaining maximum posterior estimation of three-dimensional deformation field by direct analytic method
When the equation coefficient matrix is not zero, solving a group of unique solutions, wherein the solutions are the maximum posterior estimation of the three-dimensional deformation field with high spatial resolution;
and then, performing Kriging interpolation and encryption on the three-dimensional deformation field on the time domain for the high-time frequency GPS data, so as to encrypt the quasi-GPS result into a quasi-GPS time sequence in the time domain, and obtain the three-dimensional deformation with high time resolution.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2010101067941A CN101770027B (en) | 2010-02-05 | 2010-02-05 | Ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2010101067941A CN101770027B (en) | 2010-02-05 | 2010-02-05 | Ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN101770027A CN101770027A (en) | 2010-07-07 |
| CN101770027B true CN101770027B (en) | 2012-05-16 |
Family
ID=42503006
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN2010101067941A Expired - Fee Related CN101770027B (en) | 2010-02-05 | 2010-02-05 | Ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN101770027B (en) |
Families Citing this family (64)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101964009B (en) * | 2010-09-17 | 2012-09-05 | 中煤地航测遥感局有限公司 | System and method for manufacturing 3D products based on interferometric synthetic aperture radar (INSAR) |
| RU2480788C2 (en) * | 2010-12-27 | 2013-04-27 | Общество с ограниченной ответственностью "Интеллектуальные радиооптические системы" | Radar system of remote earth sensing |
| CN102927934B (en) * | 2012-11-07 | 2015-01-28 | 中南大学 | Method for obtaining mining area earth surface three-dimensional deformation fields through single interferometric synthetic aperture radar (InSAR) interference pair |
| CN103091675B (en) * | 2013-01-11 | 2014-07-30 | 中南大学 | Mining lot exploiting and monitoring method based on interferometric synthetic aperature radar (InSAR) technology |
| CN103091676A (en) * | 2013-01-22 | 2013-05-08 | 中国矿业大学 | Mining area surface subsidence synthetic aperture radar interferometry monitoring and calculating method |
| CN103279945B (en) * | 2013-04-26 | 2015-11-25 | 北京理工大学 | A kind of interferometric phase image unwrapping method based on Quality Map daoyin technique and Branch cut |
| CN103576149B (en) * | 2013-06-05 | 2015-11-18 | 河海大学 | A kind of foundation interference radar three-dimensional deformation extraction method based on amplitude information |
| CN103454636B (en) * | 2013-09-08 | 2015-05-20 | 西安电子科技大学 | Differential interferometric phase estimation method based on multi-pixel covariance matrixes |
| CN104391296A (en) * | 2014-10-15 | 2015-03-04 | 淮海工学院 | Radar three-dimensional deformation field reconstruction technology based on general least squares adjustment |
| CN104732513A (en) * | 2014-11-27 | 2015-06-24 | 国网青海省电力公司西宁供电公司 | Landslide analytical method |
| CN104730551B (en) * | 2015-03-12 | 2017-03-22 | 北京理工大学 | Space-ground bistatic differential interferometry baseline coordinate and deformation quantity measurement method |
| CN105158760B (en) * | 2015-08-10 | 2017-04-26 | 中南大学 | Method for inverting underground fluid volume change and three dimension surface deformation using InSAR |
| CN105866777B (en) * | 2016-03-29 | 2018-10-16 | 北京理工大学 | The bistatic PS-InSAR three-dimensional deformations inversion method of the multi-period navigation satellite of multi-angle |
| CN106407714A (en) * | 2016-10-14 | 2017-02-15 | 珠海富鸿科技有限公司 | Air pollution assessment method and device based on CALPUFF system |
| CN106932777A (en) * | 2017-04-12 | 2017-07-07 | 西南石油大学 | Interfering synthetic aperture radar based on temperature baselines is to optimum option method |
| CN107123134A (en) * | 2017-05-03 | 2017-09-01 | 长安大学 | A kind of Dangerous Rock Body landslide monitoring method of feature based |
| CN107102332B (en) * | 2017-05-11 | 2019-12-06 | 中南大学 | InSAR three-dimensional earth surface deformation monitoring method based on variance component estimation and stress strain model |
| CN107329140B (en) * | 2017-07-28 | 2019-06-25 | 安徽威德萨科技有限公司 | A kind of road and bridge holistic health monitoring method |
| CN107389029B (en) * | 2017-08-24 | 2019-10-29 | 北京市水文地质工程地质大队 | A kind of surface subsidence integrated monitor method based on the fusion of multi-source monitoring technology |
| CN107918127A (en) * | 2017-11-20 | 2018-04-17 | 武汉大学 | A kind of road slope deformation detecting system and method based on vehicle-mounted InSAR |
| CN107991676B (en) * | 2017-12-01 | 2019-09-13 | 中国人民解放军国防科技大学 | Troposphere error correction method of satellite-borne single-navigation-pass InSAR system |
| CN108303735A (en) * | 2018-01-30 | 2018-07-20 | 单新建 | The earthquake deformation acquisition methods of subband interferometry based on optimized parameter setting |
| CN108507454B (en) * | 2018-03-09 | 2019-12-03 | 北京理工大学 | One kind being based on navigation satellite Bi-InSAR deformation inverted image extracting method |
| CN108594224B (en) * | 2018-03-30 | 2021-09-03 | 中国电力工程顾问集团中南电力设计院有限公司 | Three-dimensional time sequence deformation monitoring method fusing different platforms and orbit SAR data |
| CN109001731B (en) * | 2018-06-06 | 2022-11-01 | 中国电子科技集团公司第二十九研究所 | SAR interferometric phase unwrapping reference point selection method, equipment and storage medium |
| CN108983232B (en) * | 2018-06-07 | 2020-06-09 | 中南大学 | InSAR two-dimensional surface deformation monitoring method based on adjacent rail data |
| CN112285752B (en) * | 2018-10-08 | 2023-12-15 | 闽江学院 | Single-point positioning method and device with high positioning precision |
| CN109522384A (en) * | 2018-11-20 | 2019-03-26 | 广州方舆科技有限公司 | The online 3-D scanning service system and terminal laid for website |
| CN109738892B (en) * | 2019-01-24 | 2020-06-30 | 中南大学 | Mining area earth surface high-space-time resolution three-dimensional deformation estimation method |
| CN109839635B (en) * | 2019-03-13 | 2022-12-27 | 武汉大学 | Method for extracting elevation of height measurement foot points through Cryosat-2 SARIn mode L1 b-level waveform data |
| CN109917382B (en) * | 2019-03-19 | 2020-09-22 | 中国海洋大学 | Sea tide load displacement influence evaluation and correction method in coastal zone InSAR interferogram |
| CN110109106A (en) * | 2019-04-23 | 2019-08-09 | 中国电力科学研究院有限公司 | A kind of InSAR interferometric phase unwrapping method in region with a varied topography |
| CN110044327B (en) * | 2019-04-29 | 2021-10-12 | 上海颖川佳固信息工程股份有限公司 | Infrastructure settlement monitoring method and system based on SAR data and GNSS data |
| CN110055945B (en) * | 2019-05-22 | 2021-05-25 | 马培峰 | Method and device for monitoring soil consolidation settlement and related equipment |
| CN110208879B (en) * | 2019-05-28 | 2021-03-02 | 中国科学院、水利部成都山地灾害与环境研究所 | Method for forecasting easiness of rock mass landslide in freeze-thaw zone in high-altitude mountain area |
| CN110333508B (en) * | 2019-07-19 | 2021-02-19 | 中南大学 | Multisource SAR data-based joint inversion method for time-space sliding distribution after same-seismic |
| CN112749470B (en) * | 2019-10-31 | 2023-07-14 | 北京华航无线电测量研究所 | Layout optimization fitting method for structural deformation sensor |
| CN111142119B (en) * | 2020-01-10 | 2021-08-17 | 中国地质大学（北京） | Mine geological disaster dynamic identification and monitoring method based on multi-source remote sensing data |
| CN111721241A (en) * | 2020-05-08 | 2020-09-29 | 北京理工大学 | GNSS-InBSAR and GB-InSAR cross-system fusion three-dimensional deformation measurement method |
| CN113405447A (en) * | 2020-05-19 | 2021-09-17 | 湖南北斗微芯产业发展有限公司 | Track traffic deformation monitoring method, device and equipment integrating InSAR and GNSS |
| CN111522006B (en) * | 2020-06-29 | 2020-10-09 | 航天宏图信息技术股份有限公司 | Earth surface settlement monitoring method and device fusing Beidou and InSAR data |
| CN112097733A (en) * | 2020-07-28 | 2020-12-18 | 兰州交通大学 | Surface deformation monitoring method combining InSAR technology and geographic detector |
| CN111998766B (en) * | 2020-08-31 | 2021-10-15 | 同济大学 | Surface deformation inversion method based on time sequence InSAR technology |
| CN112050725A (en) * | 2020-09-14 | 2020-12-08 | 广东省核工业地质局测绘院 | Surface deformation monitoring method fusing InSAR and GPS |
| CN112540369A (en) * | 2020-11-27 | 2021-03-23 | 武汉大学 | Landslide three-dimensional deformation resolving method and system integrating GNSS and lifting rail time sequence InSAR |
| CN112540370A (en) * | 2020-12-08 | 2021-03-23 | 鞍钢集团矿业有限公司 | InSAR and GNSS fused strip mine slope deformation measurement method |
| CN112835043B (en) * | 2021-01-06 | 2023-03-21 | 中南大学 | Method for monitoring two-dimensional deformation in any direction |
| CN113176544B (en) * | 2021-03-05 | 2022-11-11 | 河海大学 | Mismatching correction method for slope radar image and terrain point cloud |
| CN113091600B (en) * | 2021-04-06 | 2022-12-16 | 长沙理工大学 | Monitoring method for monitoring deformation of soft soil foundation by utilizing time sequence InSAR technology |
| CN113343504A (en) * | 2021-08-02 | 2021-09-03 | 成都理工大学 | Three-dimensional landslide motion risk probability evaluation method |
| CN113624122B (en) * | 2021-08-10 | 2022-09-20 | 中咨数据有限公司 | Bridge deformation monitoring method fusing GNSS data and InSAR technology |
| CN113589286B (en) * | 2021-09-28 | 2021-12-14 | 中国矿业大学 | Unscented Kalman filtering phase unwrapping method based on D-LinkNet |
| CN113960596B (en) * | 2021-10-20 | 2023-05-05 | 苏州深蓝空间遥感技术有限公司 | Landslide three-dimensional deformation monitoring method based on Beidou and PS-InSAR |
| CN113740852B (en) * | 2021-11-04 | 2022-04-12 | 国网湖北省电力有限公司超高压公司 | Method for monitoring deformation of power transmission tower by ground-based SAR (synthetic aperture radar) of artificial corner reflector |
| CN114170763B (en) * | 2021-12-06 | 2022-09-23 | 山东省地质矿产勘查开发局八〇一水文地质工程地质大队 | Landslide monitoring system and method |
| CN114877798B (en) * | 2022-03-31 | 2022-12-02 | 北京建筑大学 | Vortex wave/IMU (inertial measurement unit) fused building deformation monitoring method and system |
| CN114966685B (en) * | 2022-05-24 | 2023-04-07 | 中国水利水电科学研究院 | Dam deformation monitoring and predicting method based on InSAR and deep learning |
| CN114757238B (en) * | 2022-06-15 | 2022-09-20 | 武汉地铁集团有限公司 | Method and system for monitoring deformation of subway protection area, electronic equipment and storage medium |
| CN114791273B (en) * | 2022-06-24 | 2022-09-13 | 中国铁道科学研究院集团有限公司铁道建筑研究所 | InSAR deformation monitoring result interpretation method for landslide |
| CN115358311B (en) * | 2022-08-16 | 2023-06-16 | 南昌大学 | Multi-source data fusion processing method for ground surface deformation monitoring |
| CN117109426B (en) * | 2023-08-28 | 2024-03-22 | 兰州交通大学 | Three-dimensional deformation field modeling method fusing GNSS/InSAR observation data |
| CN117331078B (en) * | 2023-10-26 | 2024-10-01 | 内蒙古至远创新科技有限公司 | Method and system for calculating differential deformation rate based on InSAR data |
| CN117213443B (en) * | 2023-11-07 | 2024-03-19 | 江苏省地质调查研究院 | Construction and updating method of ground settlement monitoring network with integration of heaves, earth and depth |
| CN118097897B (en) * | 2024-04-29 | 2024-08-13 | 泉州中科星桥空天技术有限公司 | Wide-area important geological disaster early-stage identification method based on SAR technology |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101126811A (en) * | 2007-09-29 | 2008-02-20 | 北京交通大学 | Method for detecting lakeshore and extracting lake profile from SAR image |
| US7446705B1 (en) * | 2007-10-24 | 2008-11-04 | Wisconsin Alumni Research Foundation | Method and apparatus for determining parameters for a parametric expression characterizing the phase of an acquired signal |
| CN101493520A (en) * | 2009-01-16 | 2009-07-29 | 北京航空航天大学 | SAR image variation detecting method based on two-dimension gamma distribution |
| CA2652639A1 (en) * | 2008-02-06 | 2009-08-06 | Halliburton Energy Services, Inc. | Geodesy via gps and insar integration |
| CN101634705A (en) * | 2009-08-19 | 2010-01-27 | 西安电子科技大学 | Method for detecting target changes of SAR images based on direction information measure |
| CN101634709A (en) * | 2009-08-19 | 2010-01-27 | 西安电子科技大学 | Method for detecting changes of SAR images based on multi-scale product and principal component analysis |
-
2010
- 2010-02-05 CN CN2010101067941A patent/CN101770027B/en not_active Expired - Fee Related
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101126811A (en) * | 2007-09-29 | 2008-02-20 | 北京交通大学 | Method for detecting lakeshore and extracting lake profile from SAR image |
| US7446705B1 (en) * | 2007-10-24 | 2008-11-04 | Wisconsin Alumni Research Foundation | Method and apparatus for determining parameters for a parametric expression characterizing the phase of an acquired signal |
| CA2652639A1 (en) * | 2008-02-06 | 2009-08-06 | Halliburton Energy Services, Inc. | Geodesy via gps and insar integration |
| CN101493520A (en) * | 2009-01-16 | 2009-07-29 | 北京航空航天大学 | SAR image variation detecting method based on two-dimension gamma distribution |
| CN101634705A (en) * | 2009-08-19 | 2010-01-27 | 西安电子科技大学 | Method for detecting target changes of SAR images based on direction information measure |
| CN101634709A (en) * | 2009-08-19 | 2010-01-27 | 西安电子科技大学 | Method for detecting changes of SAR images based on multi-scale product and principal component analysis |
Non-Patent Citations (2)
| Title |
|---|
| 罗海滨等.基于DInSAR方法监测南京地表沉降的结果与分析.《高技术通讯》.2008,第18卷(第04期),418-421. * |
| 罗海滨等.应用INSAR与GPS集成技术检测地表形变探讨.《遥感技术与应用》.2006,第21卷(第06期),第493-496页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN101770027A (en) | 2010-07-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN101770027B (en) | Ground surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion | |
| Catalão et al. | Merging GPS and atmospherically corrected InSAR data to map 3-D terrain displacement velocity | |
| CN112986993B (en) | InSAR deformation monitoring method based on space constraint | |
| Sumantyo et al. | Long-term consecutive DInSAR for volume change estimation of land deformation | |
| CN105487065A (en) | Time sequence satellite borne radar data processing method and device | |
| Jiang | Sentinel-1 TOPS co-registration over low-coherence areas and its application to velocity estimation using the all pairs shortest path algorithm | |
| CN109471104B (en) | Method for acquiring three-dimensional movement amount of earth surface from SAR data of two parallel tracks | |
| Mao et al. | Estimation and compensation of ionospheric phase delay for multi-aperture InSAR: An azimuth split-spectrum interferometry approach | |
| Tang et al. | High-spatial-resolution mapping of precipitable water vapour using SAR interferograms, GPS observations and ERA-Interim reanalysis | |
| Liu et al. | Correction of positional errors and geometric distortions in topographic maps and DEMs using a rigorous SAR simulation technique | |
| CN117541929A (en) | Deformation risk assessment method for large-area power transmission channel of InSAR in complex environment | |
| Liu et al. | A comparative study of DEM reconstruction using the single-baseline and multibaseline InSAR techniques | |
| Marghany | DEM reconstruction of coastal geomorphology from DINSAR | |
| Chang et al. | InSAR atmospheric distortions mitigation: GPS observations and NCEP FNL data | |
| Wang et al. | A high precision DEM extraction method based on InSAR data | |
| Deo et al. | Evaluation of interferometric SAR DEMs generated using TanDEM-X data | |
| Shiramizu et al. | Generation of a high-accuracy regional DEM based on ALOS/PRISM imagery of East Antarctica | |
| Nitti et al. | Evaluation of DEM-assisted SAR coregistration | |
| Wieczorek | EVALUATION OF DEFORMATIONS IN THE URBAN AREA OF OLSZTYN USING SENTINEL-1 SAR INTERFEROMETRY. | |
| Wu et al. | Regression analysis of errors of sar-based dems and controlling factors | |
| Krynski et al. | Estimation of height changes of GNSS stations from the solutions of short vectors and PSI measurements | |
| Belgued et al. | Application of radar space triangulation to the calibration of interferometric DEM | |
| Ng | Advanced satellite radar interferometry for ground surface displacement detection | |
| Lin | Ice surface topography digital elevation model by interferometric SAR method | |
| Shimada et al. | Repeat pass SAR interferometry of the Pi-SAR (L) for DEM generation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20120516 Termination date: 20150205 |
|
| EXPY | Termination of patent right or utility model |