CN114325706B - Distributed scatterer filtering method - Google Patents
Distributed scatterer filtering method Download PDFInfo
- Publication number
- CN114325706B CN114325706B CN202111675777.4A CN202111675777A CN114325706B CN 114325706 B CN114325706 B CN 114325706B CN 202111675777 A CN202111675777 A CN 202111675777A CN 114325706 B CN114325706 B CN 114325706B
- Authority
- CN
- China
- Prior art keywords
- pixel
- phase
- filtering
- coherence
- image
- 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.)
- Active
Links
- 238000001914 filtration Methods 0.000 title claims abstract description 56
- 238000000034 method Methods 0.000 title claims abstract description 40
- 238000009826 distribution Methods 0.000 claims abstract description 33
- 239000011159 matrix material Substances 0.000 claims abstract description 27
- 230000000694 effects Effects 0.000 claims abstract description 17
- 238000005457 optimization Methods 0.000 claims abstract description 14
- 238000003384 imaging method Methods 0.000 claims abstract description 8
- 238000012360 testing method Methods 0.000 claims description 16
- 238000004364 calculation method Methods 0.000 claims description 14
- 239000013598 vector Substances 0.000 claims description 13
- 238000012545 processing Methods 0.000 claims description 12
- 238000011160 research Methods 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 3
- 238000013441 quality evaluation Methods 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000010998 test method Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 2
- 238000005516 engineering process Methods 0.000 abstract description 14
- 230000001427 coherent effect Effects 0.000 abstract description 5
- 238000012876 topography Methods 0.000 description 15
- 238000010586 diagram Methods 0.000 description 9
- 238000005259 measurement Methods 0.000 description 5
- 230000014759 maintenance of location Effects 0.000 description 3
- 238000005065 mining Methods 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 238000007476 Maximum Likelihood Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000000354 decomposition reaction Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 229910052500 inorganic mineral Inorganic materials 0.000 description 2
- 238000012423 maintenance Methods 0.000 description 2
- 239000011707 mineral Substances 0.000 description 2
- 238000011282 treatment Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 238000000767 Anderson–Darling test Methods 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000005305 interferometry Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 230000002085 persistent effect Effects 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
- 238000004804 winding Methods 0.000 description 1
Landscapes
- Radar Systems Or Details Thereof (AREA)
Abstract
The invention relates to a distributed scatterer filtering method, which takes a long-time sequence radar complex image set as input. Estimating the space-time base line and coherence between images, removing images with poor imaging quality and image pairs with poor interference effect, and selecting the most suitable main image. And carrying out data registration based on the main image. In the space-time dimension, the distributed scatterer filtering technology identifies pixels which are statistically distributed with the time-series intensity values of the target point in the neighborhood of the target point as a homogeneous region. And respectively defining probability and distance weight by using the probability and the distance of the same distribution, and calculating a weighted sample coherence matrix in the homogeneous region as the expression of the scattering characteristics of the target point. Then image filtering is carried out based on the sample coherent matrix, intensity filtering and phase filtering are respectively carried out, and triangular phase optimization is needed after the phase filtering. And finally, taking complex information corresponding to the filtered intensity information and the optimized phase information as a result of the long-time sequence radar complex image set distributed scatterer filtering.
Description
Technical Field
The invention relates to the technical field of radar imaging, in particular to a distributed scatterer filtering method combining hypothesis testing and Gaussian distribution.
Background
The differential synthetic aperture radar interferometry (DIFFERENTIAL INSAR, D-InSAR) technology achieves the purpose of extracting the surface deformation by removing the topographic phase component in the interference phase, and is commonly used for monitoring the surface deformation caused by sudden disasters such as earthquake, landslide and the like. However, the quality of interference fringes in a differential interferogram is affected by a number of decoherence factors, including: space-time decoherence, atmospheric decoherence, doppler decoherence, bulk scattering decoherence, system noise decoherence and other reasons; while various parameter errors exist in satellite systems, such as: the introduced phase trend has certain influence on the calculation result of the surface deformation, such as the skew error, the baseline length error, the baseline inclination angle error, the orbit error and the like. For a single interferogram, the phase information is greatly influenced by accidental factors, and the difficulty of separating the topography residual error, the atmosphere effect phase and the deformation phase in the phase information is also great, so that the accuracy and the reliability of the D-InSAR in the surface deformation measurement application are limited.
The time sequence InSAR technology separates the topography residual error, the atmosphere effect phase and the deformation phase in the phase information by extracting point targets which continuously keep high coherence on a long time sequence and establishing an equation, and inversion is carried out to obtain topography residual errors and linear deformation. The influence of space-time incoherence and atmospheric effect on the ground surface deformation monitoring precision is obviously reduced, and the method is widely applied to the fields of ground surface subsidence in urban areas and mining areas, ground surface deformation caused by geological disasters or development and the like.
Typical high coherence point based time series InSAR techniques exist as the permanent scatterer technique (PERSISTENT SCATTERER, PS-InSAR) method, the Coherent target analysis (Coherent TARGETANALYSIS, CTA) method, the spatio-temporal unwrapping network (Spatio-Temporal Unwrapping Network, STUN) method, and the interference point target analysis (Interferometric Point TARGETANALYSIS, IPTA) method. The methods aim at high coherence points, build a model on a long-time sequence interferogram, and are widely applied to the fields related to deformation monitoring, including urban ground subsidence caused by exploitation of underground water, construction of public facilities such as traffic pipelines and the like, mining area environment geological problems caused by exploitation of underground mineral resources, and ground surface deformation caused by geological disasters and ground surface deformation caused by geological structures.
Techniques based on a model of the permanent scatterer target rely on high coherence point densities and are limited in application in non-urban areas. In non-urban areas, the buildings and single large bare rock are less, and bare land and low vegetation cover more. The resolution units corresponding to the ground object types are mostly free of dominant point-shaped scatterers, and the backscattering coefficients of all the scatterers are approximately the same and are called distributed scatterers. The research scholars at home and abroad combine the scattering characteristics of the ground feature, and combine the distributed scatterer on the basis of the permanent scatterer, so that the application of the time sequence InSAR technology in non-urban areas is expanded.
The method divides SAR data into a plurality of sets according to the principle that the space base lines in the sets are shortest and the space base lines among the sets are longest, so that the quality and the number of interferograms are increased, and the space density and the time sampling rate of deformation points are improved; squeeSAR, identifying homogeneous pixels through hypothesis testing, primarily extracting a distributed scatterer, then performing feature decomposition on the phase of the distributed scatterer, selecting an optimal feature as a value of the distributed target, and processing the optimal feature and the point target together, so that the space density of deformation points is greatly improved; in the QPS technology, in a time sequence interference diagram, the principle of maximizing the time coherence is adopted, and an interference diagram subset combination is set, so that the topography residual phase and the deformation phase of a target are inverted.
The research or application of the current distributed scatterer filtering technology mostly filters and optimizes the central point based on the homogeneous region in a single estimation window, and mostly focuses on the extraction and decomposition method of the sample coherence matrix information. However, in a complex ground object type distribution scene, a single estimation window is not beneficial to the maintenance of image details, and on the premise that the obtained sample coherence matrix is the result under the average principle, the scattering characteristics of the ground object corresponding to the center pixel cannot be accurately described.
Disclosure of Invention
In order to solve the technical problems, the invention improves a single estimation window on the basis of the prior theory, introduces probability weight and distance weight to obtain a weighted sample coherent matrix, and realizes the treatments of homogeneous region identification, complex image filtering, triangular phase optimization and the like on the basis of the weighted sample coherent matrix by utilizing an advanced computer technology to filter a distributed scatterer.
The invention designs a distributed scatterer filtering method composed of homogeneous region identification, complex image filtering and triangular phase optimization based on the scattering characteristics of the distributed scatterer and by combining hypothesis test and Gaussian distribution.
The technical scheme of the invention is a distributed scatterer filtering method, which comprises the following steps:
step1, taking a long-time sequence radar complex image set, track data and a research area reference DEM as inputs;
Step 2, estimating space-time base lines and coherence between images, preprocessing and screening data, selecting main images, and registering the data with the main images as the reference;
step 3, sequencing the registered radar complex images according to imaging time, so that subsequent processing is facilitated;
step 4, setting a determined sliding estimation window, and identifying homogeneous areas of the radar image pixel by pixel;
Step 5, respectively defining probability and distance weight according to the probability and the distance of the same distribution, and calculating a weighted sample coherence matrix in the homogeneous region as expression of scattering characteristics of the target point;
Step 6, in the estimation window, performing self-adaptive spatial intensity filtering and phase filtering processing on the central pixel based on the homogeneous region and the PG filter to obtain filtered intensity and phase;
step 7, triangular phase optimization is needed after filtering;
step 8, setting a homogeneous pixel quantity threshold and a fitting goodness threshold, using pixels meeting the conditions as distributed scatterers, and updating the distributed scatterers by using the optimized phases and the optimized filtering intensities to obtain an updated radar image data set;
and 9, carrying out joint calculation on the permanent scatterer and the distributed scatterer serving as a high-coherence point target on the updated radar image data set according to the PS-InSAR flow.
The beneficial effects are that:
According to the invention, the probability weight and the distance weight are introduced into the adaptive spatial filtering method based on the estimation window, the scattering characteristics of the distributed scatterers are combined, the weighted sample coherence matrix is calculated, and based on the same-quality area, the image filtering, the triangular phase optimization and other treatments are carried out, so that a technical flow of the distributed scatterer filtering is formed, the filtered weighted sample coherence matrix can more accurately represent the scattering characteristics of the central pixel corresponding to the target ground object, the problem of loss of the filtered image details under the complex ground object type distribution scene is solved, and the PS point information is reserved as much as possible while the distributed scatterers are identified, filtered and noise is reduced. The filtered long-time sequence radar complex image set can be used for jointly solving the time sequence InSAR technology of the distributed scatterer and the permanent scatterer, and the time sequence InSAR technology based on the permanent scatterer cannot extract enough measuring points when the ground surface deformation is inverted in a non-urban area.
Drawings
FIG. 1 (a) is a flow chart of a method for jointly resolving a time series InSAR of a permanent scatterer and a distributed scatterer in accordance with the present invention;
FIG. 1 (b) is a schematic flow chart of a method for filtering a distributed scatterer according to the present invention;
FIG. 2 is pvalue neighborhood matrix and statistically co-distributed mask: (a) assuming a probability matrix, (b) counting a uniformly distributed mask, (c) communicating regions of homogeneity;
FIG. 3 is a schematic illustration of a study area;
FIG. 4 (a) is a time series interferogram;
FIG. 4 (b) is a time series coherence coefficient map;
FIG. 5 is a diagram of interference fringes and a diagram of coherence coefficient comparing the original radar complex image with the model filtering result: (a) an original interference fringe pattern, (b) an original coherence coefficient pattern, (c) a filtered interference fringe pattern, (d) a filtered coherence coefficient pattern;
FIG. 6 is a measurement point contrast of an original image dataset versus a filtered image dataset inversion: (a) measuring points inverted from the raw image dataset; (b) filtering the inverted measurement points of the image dataset;
FIG. 7 is a comparison of residual terrain phase errors for an inversion of an original image dataset and a filtered image dataset: (a) A residual topography phase error inverted from the original image dataset, (b) a residual topography phase error inverted from the filtered image dataset;
FIG. 8 is an atmospheric and orbital phase error contrast for an original image dataset inverted with a filtered image dataset: (a) An atmospheric and orbital phase error of inversion of the original image dataset, (b) an atmospheric and orbital phase error of inversion of the filtered image dataset.
Detailed Description
The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention, and it is apparent that the described embodiments are only some embodiments of the present invention, but not all embodiments, and all other embodiments obtained by those skilled in the art without the inventive effort based on the embodiments of the present invention are within the scope of protection of the present invention.
On a long-time sequence radar complex image set, a point target which continuously keeps high coherence is a stable point target. Based on a stable point target, a deformation model is established, so that the interference of space-time incoherence and atmospheric effect can be overcome, and more accurate terrain residual errors and deformation phases can be obtained through inversion.
As shown in fig. 1 (a), distributed scatterer filtering is typically used as a pre-processing part in a time series InSAR process that jointly solves for a permanent scatterer and a distributed scatterer. In the area with low distribution density of the permanent scatterers, after the distributed scatterers are filtered based on the radar complex image sequence, image noise is reduced, coherence is improved, and then the permanent scatterers and the distributed scatterers are used as high-coherence point targets for joint calculation.
As shown in fig. 1 (b), a distributed scatterer filtering method of the present invention includes the steps of:
and step1, taking a long-time sequence radar complex image set, track data and a research area reference DEM as inputs.
And 2, estimating space-time baselines and coherence among the images, removing the image with poor imaging quality and the image pair with poor interference effect, selecting the most suitable main image, and carrying out data registration based on the main image.
And step 3, sequencing the registered radar complex images according to imaging time, so that subsequent processing is facilitated.
And 4, setting a determined sliding estimation window, and identifying homogeneous areas of the radar image pixel by pixel. The specific method comprises the following steps: in each sliding window, judging whether a time sequence amplitude vector corresponding to a time sequence complex vector corresponding to each pixel in the sliding window and a central pixel of the pixel is from a certain same distribution overall by a double-sample hypothesis test method, and if the hypothesis cannot be overturned, the pixel is a statistical same-distribution pixel of the central pixel. After the statistical uniform distribution pixel mask of the central pixel of the sliding window is obtained by using the double-sample detection method, the connected region identification is required to be carried out on the binary mask image. The identified pixel set consistent with the center pixel label is the homogeneous region of the center pixel, and the number of pixels consistent with the center pixel label is the homogeneous pixel number (StatisticallyHomogeneous Pixel, SHP) of the center pixel.
Step 5, defining a Gaussian filter based on the hypothesis test p value, (p-value gaussian filter is hereinafter referred to as PG filter), and for a sliding estimation window of (2m+1) x (2n+1), the filter calculation formula of the center pixel (coordinates are set as (0, 0)) is as follows:
Where x ij represents the pixel amplitude/phase value of (i, j), P (i, j) represents the weight corresponding to the hypothesis test P value, and G (i, j) represents the gaussian kernel weight. The calculation formula of P (i, j) is as follows:
Where p ij represents the hypothetical test p value for the pixel of (i, j) and the center pixel, α represents a given level of saliency, β is an exponential coefficient, the larger the coefficient, the greater the effect of the hypothetical test p value on the filter and vice versa. G (i, j) is a two-dimensional Gaussian kernel, and the calculation formula is as follows:
Wherein sigma represents a Gaussian standard deviation, which reflects the discrete degree of weight data in the template, the smaller sigma is, the more concentrated the distribution is, the larger the center coefficient of the template is, the smaller the surrounding coefficient is, and the smaller the smoothness is; the larger the σ, the closer the template center coefficient is to the perimeter coefficient, and the gaussian kernel approximates the mean template when σ→+ -infinity.
And 6, in the estimation window, performing self-adaptive spatial intensity filtering processing on the central pixel based on the homogeneous region and the PG filter to obtain the filtered intensity.
And 7, in the estimation window, carrying out phase optimization on the center pixel based on the homogeneous region and the PG filter. The specific method comprises the following steps: estimating a sample covariance matrix of the center pixel based on the homogeneous region and the PG filter; and carrying out normalization processing on all pixel amplitude values in the homogeneous region of each same phase based on the sample covariance matrix, and estimating a sample coherence matrix of the central pixel. And optimizing the time sequence phase vector of the center pixel according to the coherence coefficient and the interference phase in the sample coherence matrix to obtain an optimized phase. And meanwhile, comparing the optimized interference phase vector with the original interference phase vector to evaluate the quality of goodness of fit (GOF).
And 8, according to the set homogeneous pixel quantity threshold and the set goodness-of-fit threshold, taking the pixels meeting the conditions as a distributed scatterer, and updating the distributed scatterer by using the optimized phase and the optimized filtering intensity. And obtaining an updated radar image data set.
And 9, carrying out joint calculation on the permanent scatterer and the distributed scatterer serving as a high-coherence point target on the updated radar image data set according to the PS-InSAR flow.
According to the embodiment of the invention, the distributed scatterer filtering part is realized by adopting a python programming language, and the realization process can be summarized into the following steps:
1. Setting parameters such as an estimation window size, an assumption checking probability value threshold, a homogeneous pixel quantity threshold, a fitting goodness index threshold, a significance level alpha, an assumption checking probability value weight coefficient, a Gaussian distribution standard deviation and the like;
2. Reading a radar complex image time sequence registered to the main image, and setting a processing range;
3. Starting to calculate the complex value of the central pixel in each estimation window after intensity filtering and phase optimization;
4. And taking the current pixel as a center, and identifying a homogeneous region of the center pixel in the estimation window according to the set size of the estimation window. The specific method comprises the following steps: and (3) in the current sliding window, carrying out Anderson-Darling inspection on the central pixel and the neighborhood pixels in sequence to obtain a sliding window hypothesis probability (p-value) matrix, as shown in fig. 2 (a), and obtaining a statistical equal-distribution mask of the central pixel according to a given significance level alpha, as shown in fig. 2 (b). And carrying out connected region identification on the pixel mask images with the same statistical distribution. The identified pixel set consistent with the center pixel label is the homogeneous region of the center pixel, as shown in fig. 2 (c), and the number of pixels consistent with the center pixel label is the homogeneous number of pixels of the center pixel. The specific method for Anderson-Darling test comprises the following steps: judging whether the time sequence amplitude vectors X and Y corresponding to the time sequence complex vectors corresponding to each pixel and the central pixel in the sliding window come from a certain same distribution overall, and if the assumption can not be overturned, the pixel is the statistical same-distribution pixel of the central pixel. The test steps are as follows:
a) Combining X and Y and sorting from small to large to obtain a combined sample z= [ Z 1,z2,...,zk ], k=2n;
b) Taking the double weighted sum of the square differences of the single sample and the combined sample as a statistic, the calculation formula is as follows:
Where M i is the set { x: x is less than or equal to z i }, which is X and X is a number of elements of (a).
C) At a given significance level α, look at the AD double sample threshold table for small samples, if AD > AD α (N, N), then H 0 is rejected, i.e., the two samples are considered to be from a population of different distributions.
5. The method comprises the steps of combining hypothesis probability, gaussian weight, weight coefficient and significance level alpha corresponding to the Gaussian weight, calculating weight masks used for filtering at each position of a sliding window, normalizing the weights, and defining a Gaussian filter (p-value gaussian filter, hereinafter referred to as PG filter) based on a hypothesis test p value, wherein the specific method comprises the following steps: for a sliding estimation window of (2m+1) × (2n+1), the filter calculation formula for its center pixel (coordinates set to (0, 0)) is as follows:
Where x ij represents the pixel amplitude/phase value of (i, j), P (i, j) represents the weight corresponding to the hypothesis test P value, and G (i, j) represents the gaussian kernel weight. The calculation formula of P (i, j) is as follows:
Where p ij represents the hypothetical test p value for the pixel of (i, j) and the center pixel, α represents a given level of saliency, β is an exponential coefficient, the larger the coefficient, the greater the effect of the hypothetical test p value on the filter and vice versa. G (i, j) is a two-dimensional Gaussian kernel, and the calculation formula is as follows:
Wherein sigma represents a Gaussian standard deviation, which reflects the discrete degree of weight data in the template, the smaller sigma is, the more concentrated the distribution is, the larger the center coefficient of the template is, the smaller the surrounding coefficient is, and the smaller the smoothness is; the larger the σ, the closer the template center coefficient is to the perimeter coefficient, and the gaussian kernel approximates the mean template when σ→+ -infinity.
6. In the estimation window, performing self-adaptive spatial intensity filtering processing on the central pixel based on the homogeneous region and the PG filter to obtain filtered intensity;
7. Within the estimation window, the center pixel is phase optimized based on the homogeneous region and the PG filter. The specific method comprises the following steps: estimating a sample covariance matrix of the center pixel based on the homogeneous region and the PG filter; and carrying out normalization processing on all pixel amplitude values in the homogeneous region of each same phase based on the sample covariance matrix, and estimating a sample coherence matrix of the central pixel. According to the coherence coefficient and the interference phase in the sample coherence matrix, the radar image phase value with the first imaging time is set to be 0, according to N (N-1)/2 interference phases in the weighted sample coherence matrix, a quasi-Newton method BFGS algorithm in an optimization algorithm is adopted, and the maximum likelihood estimation value of the N-1 time sequence interference phases is solved, so that the maximum likelihood estimation value meets the time sequence triangle principle. Meanwhile, as the winding phase value is [ -II, II ], the boundary setting is needed in the optimization process; obtaining the phase after the time sequence phase vector of the central pixel is optimized;
8. Calculating the best fitting goodness according to the interference phases before and after optimization; the goodness of fit is expressed as the average of the real part of the interference phase difference as follows:
wherein N represents the number of time-series images, Represents the interference phase of the nth image and the kth image, θ n-θk represents the interference phase of the nth image and the kth image after optimization, and Re represents the real part of complex number. From the above equation, the result of the goodness of fit is [ -1,1], the closer the result is to 1, the better the fitting effect is, the closer to-1, the worse the fitting effect is.
9. And replacing the original complex image values of the positions corresponding to the pixels meeting the requirements by the filtered intensity values and the filtered phase values according to the set pixel quantity threshold value of the homogeneous region and the best fitting goodness threshold value, so as to obtain an updated long-time sequence radar complex image set.
According to the embodiment of the invention, a certain mining area in the northwest of China is taken as an application and verification case area of the invention. 115 DEG 54'-119 DEG 15' E,40 DEG 12'-42 DEG 37' N, east-west span about 280km, north-south span about 270km and total area about 39500km 2 at the Dewar city of Hebei province. The urban topography of the Deck is reduced from northwest to southeast, the north part is an inner Mongolia plateau, the middle part is shallow mountain distribution, the south part belongs to the Yanshan mountain, the low and medium mountains, the basin between mountains and the valley of the river are the main, and the forest coverage rate reaches 55.8 percent. The underway city is bounded by a Kangbao-enclosure breaking zone with two class I building blocks spanning the medium-oriented base station and the inner Mongolian-great Khingan geosyncline fold. The formation in the region is complete, the ore forming condition is good, and the mineral resources are rich. The elevation range of the selected study area is 398-596m, the average elevation 468m, the maximum gradient 88 degrees, and the average gradient 25 degrees. As shown in FIG. 3, the long-time sequence radar image set for the case area of the schematic diagram of the investigation region is derived from the high-resolution Sentinel-1 radar image data (https:// scihub. Crunicus. Eu/dhus/#/home) of the European space agency (European SPACEAGENCY, ESA) and its precise orbit data (https:// qc. Sentinel1.Eo. Esa. Int/aux_ poeorb /). The initial data set is 42 scenes of Sentinel-1A track-lifting IW mode images with imaging time between 10 months in 2018 and 10 months in 2020, and 25 scenes of images meeting the requirements are reserved after space-time base line estimation. The auxiliary DEM for removing the terrain phase is derived from the spacecraft radar terrain mapping mission (Shuttle Radar Topography Mission, SRTM, https:// earthoxplorer. Usgs. Gov /) resolution of about 30m.
The distributed scatterer filtering is usually used as a preprocessing part in the process of the time sequence InSAR for jointly resolving the permanent scatterer and the distributed scatterer, has the filtering effects of reducing noise and improving coherence, can extract more high-coherence points meeting requirements on a long time sequence compared with an original data set, improves the density and quality of measuring points of the time sequence InSAR technology, establishes a deformation model, separates the topography residual error, the atmospheric effect phase and the deformation phase in phase information, and performs inversion to obtain topography residual errors and linear deformation.
The operation of the distributed scatterer filtering technology is to obtain an updated long-time sequence radar image dataset, and in order to better show the effects of the technology on noise reduction, coherence improvement and detail maintenance, the method is mainly shown by the effects of a filtered radar complex image intensity diagram and interference patterns and coherence coefficient diagrams generated by interference thereof on time and space angles. Showing a time sequence interference pattern and a time sequence coherence coefficient pattern in a time angle; at a spatial angle, a time-series average intensity map and a time-series average coherence coefficient map are shown. A time series interference diagram shown in fig. 4 (a) and a time series coherence coefficient diagram shown in fig. 4 (b) shown in the model operation result;
the invention compares and verifies the interference fringe pattern and the coherence coefficient pattern obtained by generating the original radar complex image and the filter radar complex image, as shown in fig. 5 (a) - (d) and table 1, and discusses the effectiveness and stability of the distributed scatterer filter model on the extraction of the distributed scatterer. The interference fringes in fig. 5 (a) - (d) reflect the distribution of the unreeled interference phases, the coherence coefficient reflects the quality of the interference phases, and the effect of noise reduction and coherence improvement of the filtering can be seen by visual effect. In order to evaluate the distortion condition of the filtered image, 5 image quality evaluation indexes of mean value, standard deviation, signal-to-noise ratio, equivalent vision number and edge retention coefficient are selected to compare the image quality of the interference fringe image before and after filtering, as shown in table 1, the mean value and variance of the original interference fringe image and the average value and variance of the filtered interference fringe image are respectively different by 0.0138 and-0.0305, which indicates that the statistical characteristic of the images before and after filtering is not greatly changed, the equivalent vision number is respectively 4.65E-04 and is 9.07E-04, which indicates that the filtered interference fringe image has a certain noise reduction effect after filtering, and compared with the signal-to-noise ratio and the edge retention coefficient of the original interference fringe image, which indicates that the distortion of the image is smaller and the edge retention is better. FIG. 5 is a comparison of the interference fringe pattern and the coherence coefficient pattern generated by the filtering result of the original radar complex image and the model.
Table 1 interference fringe pattern image quality evaluation index contrast generated by original radar complex image and model filtering result
In the case, the long-time sequence radar image data sets before and after the distributed scatterer filtration are subjected to earth surface deformation inversion respectively, the density of measuring points, the residual topography phase error, the atmosphere and the orbit error phase are subjected to contrast verification, and the improvement condition of the distributed scatterer technology on the density and the quality of the measuring points in the time sequence InSAR technology is discussed.
Comparing the original image dataset with the inverted measuring points of the filtered image dataset in fig. 6, and respectively carrying out the surface deformation inversion on the long-time sequence radar image datasets before and after the filtering of the distributed scatterer in fig. 6 (a) - (b), wherein the obtained surface deformation rate distribution obtained by inversion is more consistent and is respectively-14.92-12.23 mm/y and-12.36-11.87 mm/y, and the obtained surface deformation rate distribution also has stronger consistency in spatial distribution. Table 2 shows that the measurement point density obtained from the filtered image dataset is increased by about 6 times compared to the original image dataset, and it can be seen that the distributed scatterer filter can significantly improve the high coherence point density.
Table 2 measurement point contrast for inversion of raw image dataset and filtered image dataset
FIG. 7 is a comparison of the phase error of the residual terrain inverted from the original image dataset and the filtered image dataset; fig. 7 (a) - (b) are respectively the residual topography phase error distribution situations of the surface deformation inversion of the long-time sequence radar image data sets before and after the distributed scatterer filtration, and it can be seen that the residual topography phase errors obtained by inversion have stronger consistency in spatial distribution. Table 3 shows that the density of the measuring points obtained by filtering the image data set is reduced compared with the original image data set, and the maximum value and the average value of the residual topography phase error are reduced, so that the inversion accuracy of the distributed scatterer filtering can be improved.
TABLE 3 contrast of residual topography phase errors for inversion of raw image dataset and filtered image dataset
Fig. 8 (a) - (b) are respectively the atmospheric and orbit phase error distribution conditions of the surface deformation inversion of the long-time sequence radar image data sets before and after the filtering of the distributed scatterer, and it can be seen that the atmospheric and orbit phase errors obtained by inversion have stronger consistency in time distribution and space distribution. Compared with the original image data set, the error range of the atmospheric and orbit phase errors obtained by filtering the image data set is reduced, and the inversion accuracy of the atmospheric and orbit phases is improved by the distributed scatterer filtering.
While the foregoing has been described in relation to illustrative embodiments thereof, so as to facilitate the understanding of the present invention by those skilled in the art, it should be understood that the present invention is not limited to the scope of the embodiments, but is to be construed as limited to the spirit and scope of the invention as defined and defined by the appended claims, as long as various changes are apparent to those skilled in the art, all within the scope of which the invention is defined by the appended claims.
Claims (4)
1. A method for filtering a distributed scatterer, comprising the steps of:
step1, taking a long-time sequence radar complex image set, track data and a research area reference DEM as inputs;
Step 2, estimating space-time base lines and coherence between images, preprocessing and screening data, selecting main images, and registering the data with the main images as the reference;
step 3, sequencing the registered radar complex images according to imaging time, so that subsequent processing is facilitated;
step 4, setting a determined sliding estimation window, and identifying homogeneous areas of the radar image pixel by pixel;
Step 5, respectively defining probability and distance weight according to the probability and the distance of the same distribution, and calculating a weighted sample coherence matrix in the homogeneous region as expression of scattering characteristics of the target point;
Step 6, in the estimation window, performing self-adaptive spatial intensity filtering and phase filtering processing on the central pixel based on the homogeneous region and the PG filter to obtain filtered intensity and phase;
step 7, triangular phase optimization is needed after filtering;
step 8, setting a homogeneous pixel quantity threshold and a fitting goodness threshold, using pixels meeting the conditions as distributed scatterers, and updating the distributed scatterers by using the optimized phases and the optimized filtering intensities to obtain an updated radar image data set;
and 9, carrying out joint calculation on the permanent scatterer and the distributed scatterer serving as a high-coherence point target on the updated radar image data set according to the PS-InSAR flow.
2. The method of claim 1, wherein the step 4 of setting a determined sliding estimation window, and performing pixel-by-pixel homogeneous region identification on the radar image comprises the following steps:
Judging whether each pixel in the sliding window is from a certain same distribution ensemble or not according to a double-sample hypothesis test method and a time sequence amplitude vector corresponding to a time sequence complex vector corresponding to a central pixel of the pixel in the sliding window, if the assumption cannot be overturned, the pixel is a statistical same-distribution pixel of the central pixel, and after obtaining a statistical same-distribution pixel mask of the central pixel of the sliding window by using the double-sample hypothesis test method, carrying out connected region identification on a binary mask image; the identified pixel set consistent with the center pixel label is the homogeneous region of the center pixel, and the number of pixels consistent with the center pixel label is the homogeneous number of the center pixel.
3. The method according to claim 1, wherein the step 5 is to define probability and distance weight with probability and distance of the same distribution, respectively, and calculate weighted sample coherence matrix in homogeneous region as expression of scattering property of target point, specifically as follows:
Defining a Gaussian filter based on an assumption test p value, and for a sliding estimation window of (2m+1) × (2n+1), m, n is a natural number, the coordinates of a central pixel are set to be (0, 0), and a filter calculation formula of the central pixel is as follows:
Wherein x ij represents the pixel amplitude value/phase value of the (i, j) position, P (i, j) represents the weight corresponding to the hypothesis test P value, G (i, j) represents the gaussian kernel weight, and the calculation formula of P (i, j) is as follows:
Wherein p ij represents the hypothetical test p value for the pixel at the (i, j) position and the center pixel, α represents a given level of saliency, β is an exponential coefficient, the larger the coefficient, the greater the effect of the hypothetical test p value on the filter and vice versa, G (i, j) is a two-dimensional gaussian kernel, the following formula is calculated:
Wherein sigma represents a Gaussian standard deviation, which reflects the discrete degree of weight data in the template, the smaller sigma is, the more concentrated the distribution is, the larger the center coefficient of the template is, the smaller the surrounding coefficient is, and the smaller the smoothness is; the larger the σ, the closer the template center coefficient is to the perimeter coefficient, and the gaussian kernel approximates the mean template when σ→+ -infinity.
4. The method of claim 1, wherein the filtering in step 7 requires triangular phase optimization, and the method specifically comprises:
In the estimation window, the phase optimization is carried out on the center pixel based on the homogeneous region and the PG filter, and the specific method is as follows: estimating a sample covariance matrix of the center pixel based on the homogeneous region and the PG filter; based on the sample covariance matrix, carrying out normalization processing on all pixel amplitude values in the same homogeneous region of each same phase, estimating a sample coherence matrix of a central pixel, optimizing a time sequence phase vector of the central pixel according to a coherence coefficient and an interference phase in the sample coherence matrix to obtain an optimized phase, and simultaneously comparing the optimized interference phase vector with an original interference phase vector to carry out fitting goodness-of-quality evaluation.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111675777.4A CN114325706B (en) | 2021-12-31 | 2021-12-31 | Distributed scatterer filtering method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111675777.4A CN114325706B (en) | 2021-12-31 | 2021-12-31 | Distributed scatterer filtering method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114325706A CN114325706A (en) | 2022-04-12 |
| CN114325706B true CN114325706B (en) | 2024-10-25 |
Family
ID=81023352
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111675777.4A Active CN114325706B (en) | 2021-12-31 | 2021-12-31 | Distributed scatterer filtering method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114325706B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114779249B (en) * | 2022-04-19 | 2024-08-20 | 东北大学 | GB-SAR (Global positioning System-specific absorption Rate) interference image filtering method and system |
| CN115143877B (en) * | 2022-05-30 | 2023-04-25 | 中南大学 | SAR offset tracking method, device, equipment and medium based on strong scattering amplitude suppression and outlier identification |
| CN116047519B (en) * | 2023-03-30 | 2023-06-16 | 山东建筑大学 | Point selection method based on synthetic aperture radar interferometry technology |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106950556A (en) * | 2017-05-03 | 2017-07-14 | 三亚中科遥感研究所 | Heritage area deformation monitoring method based on distributed diffusion body sequential interference SAR technology |
| CN107392861A (en) * | 2017-06-29 | 2017-11-24 | 南京航空航天大学 | A kind of rarefaction representation SAR image method for reducing speckle based on Gauss ratio mixed model |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| IT1394733B1 (en) * | 2009-07-08 | 2012-07-13 | Milano Politecnico | PROCEDURE FOR FILTERING INTERFEROGRAMS GENERATED BY IMAGES ACQUIRED ON THE SAME AREA. |
| CN108051810B (en) * | 2017-12-01 | 2020-06-09 | 南京市测绘勘察研究院股份有限公司 | InSAR distributed scatterer phase optimization method |
-
2021
- 2021-12-31 CN CN202111675777.4A patent/CN114325706B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106950556A (en) * | 2017-05-03 | 2017-07-14 | 三亚中科遥感研究所 | Heritage area deformation monitoring method based on distributed diffusion body sequential interference SAR technology |
| CN107392861A (en) * | 2017-06-29 | 2017-11-24 | 南京航空航天大学 | A kind of rarefaction representation SAR image method for reducing speckle based on Gauss ratio mixed model |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114325706A (en) | 2022-04-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114325706B (en) | Distributed scatterer filtering method | |
| CN109558859B (en) | Mining area distribution information extraction method and system based on DInSAR and DCNN | |
| Modava et al. | Hierarchical coastline detection in SAR images based on spectral‐textural features and global–local information | |
| CN113960595A (en) | Surface deformation monitoring method and system | |
| Wang et al. | Construction and optimization method of the open-pit mine DEM based on the oblique photogrammetry generated DSM | |
| CN113919226B (en) | Mining vegetation ecological cumulative effect disturbance range identification method based on weight | |
| Jiang et al. | Delineation of built-up land change from SAR stack by analysing the coefficient of variation | |
| CN113281749A (en) | Time sequence InSAR high-coherence point selection method considering homogeneity | |
| CN114091274A (en) | Landslide susceptibility evaluation method and system | |
| Kwak et al. | A new approach for rapid urban flood mapping using ALOS-2/PALSAR-2 in 2015 Kinu River Flood, Japan | |
| CN115079172A (en) | MTInSAR landslide monitoring method, equipment and storage medium | |
| CN111623749B (en) | Novel railway goaf boundary extraction method based on D-InSAR | |
| CN117333468A (en) | Flood disaster monitoring method for multi-mode time sequence PolSAR image | |
| Feng et al. | A hierarchical network densification approach for reconstruction of historical ice velocity fields in East Antarctica | |
| CN113238228A (en) | Level constraint-based InSAR three-dimensional surface deformation acquisition method, system and device | |
| Zhang et al. | Using a phase-congruency-based detector for glacial lake segmentation in high-temporal resolution sentinel-1a/1b data | |
| Xu et al. | Mapping and analyzing the annual dynamics of tidal flats in the conterminous United States from 1984 to 2020 using Google Earth Engine | |
| Du et al. | Extracting water body data based on SDWI and threshold segmentation: A case study in permafrost area surrounding Salt Lake in Hoh Xil, Qinghai-Tibet Plateau, China | |
| Wu et al. | ESRGAN-based DEM super-resolution for enhanced slope deformation monitoring in lantau island of Hong Kong | |
| CN118196637A (en) | Potential landslide identification method based on InSAR deformation and influence factor coupling | |
| CN117310705B (en) | Flood disaster rapid detection method based on dual-polarized SAR image | |
| Sui et al. | Processing of multitemporal data and change detection | |
| CN112529837A (en) | Remote sensing image change detection algorithm based on coupling discrimination feature self-learning network | |
| CN117274821A (en) | Multi-polarization SAR farmland flood detection method and system considering rainfall influence | |
| Liu et al. | Research on automatic recognition of active landslides using InSAR deformation under digital morphology: A case study of the Baihetan reservoir, China |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |