CN114325706A - Distributed scatterer filtering method - Google Patents
Distributed scatterer filtering method Download PDFInfo
- Publication number
- CN114325706A CN114325706A CN202111675777.4A CN202111675777A CN114325706A CN 114325706 A CN114325706 A CN 114325706A CN 202111675777 A CN202111675777 A CN 202111675777A CN 114325706 A CN114325706 A CN 114325706A
- Authority
- CN
- China
- Prior art keywords
- pixel
- filtering
- phase
- image
- distributed
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Abstract
The invention discloses a distributed scatterer filtering method, which takes a long-time sequence radar complex image set as input. And (4) pre-estimating a space-time baseline and coherence between the images, eliminating the images with poor imaging quality and the image pairs with poor interference effect, and selecting the most appropriate 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 sequence intensity value of the target point in the neighborhood of the target point and serve as homogeneous areas. And 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 the expression of the scattering characteristics of the target point. And then, image filtering is carried out based on a sample coherent matrix, intensity filtering and phase filtering are respectively carried out, and triangular phase optimization is required after the phase filtering. And finally, taking the filtered intensity information and the complex information corresponding to the optimized phase information as a result of filtering the distributed scatterers of the long-time sequence radar complex image set.
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 (D-InSAR) technique achieves the purpose of extracting surface deformation by removing a terrain phase component in an interference phase, and is commonly used for monitoring surface deformation caused by sudden disasters such as earthquakes, landslides and the like. However, the quality of the interference fringes in the differential interferogram is affected by a number of decoherence factors, including: time-space decorrelation, atmospheric decorrelation, Doppler decorrelation, volume scattering decorrelation, system noise decorrelation and the like; meanwhile, various parameter errors exist in the satellite system, such as: the introduced phase trend also has certain influence on the resolving result of the earth surface deformation, such as the slope distance error, the base line length error, the base line inclination angle error, the orbit error and the like. For a single interferogram, the phase information of the interferogram is greatly influenced by accidental factors, and the difficulty of separating terrain residual errors, atmospheric effect phases and deformation phases in the phase information is also high, 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 terrain residual error, atmospheric effect phase and deformation phase in phase information by extracting point targets continuously keeping high coherence on a long time sequence and establishing an equation, and the terrain residual error and linear deformation are obtained through inversion. The method obviously reduces the influence of space-time incoherent and atmospheric effects on the monitoring precision of the surface deformation, and is widely applied to the fields of surface sedimentation in urban areas and mining areas, surface deformation caused by occurrence or development of geological disasters and the like.
The existing typical time series InSAR technology based on high coherence points includes a permanent Scatterer technology (PS-InSAR) method, a Coherent Target Analysis (CTA) method, a space-time unwrapting Network (STUN) method, and an Interferometric Point Target Analysis (IPTA) method. The method is used for establishing a model on a long-time sequence interference diagram by taking high coherence points as targets, and is widely applied to the fields related to deformation monitoring, including urban ground subsidence caused by exploitation of underground water and construction of public facilities such as traffic pipelines, mining area environment geological problems caused by exploitation of underground mineral resources, and surface deformation caused by geological disasters and surface deformation caused by geological structures.
Techniques based on permanent scatterer object models rely on high coherence point densities and have limited application in non-urban areas. In non-urban areas, buildings and single large bare rocks are few, and bare land and low vegetation cover more. Most of resolution units corresponding to the surface feature types do not have dominant point-like scatterers, and backscattering coefficients of all scatterers are approximately the same, so that the scatterers are called distributed scatterers. Researchers at home and abroad combine scattering characteristics of ground objects and distributed scatterers on the basis of permanent scatterers, and the application of the time series InSAR technology in non-urban areas is expanded.
The method divides SAR data into a plurality of sets according to the principle that a space baseline in a set is shortest and a space baseline in a set is longest, increases the quality and quantity of interferograms, and improves the space density and time sampling rate of deformation points; the method comprises the steps of identifying homogeneous pixels through hypothesis testing, preliminarily extracting distributed scatterers, then performing feature decomposition on phases of the distributed scatterers, selecting optimal features as values of the distributed targets, and processing the optimal features and the point targets together, so that the space density of deformation points is greatly improved; the QPS technology is used for setting an interferogram subset combination in a time sequence interferogram by using the principle of time coherence maximization so as to invert the terrain residual phase and deformation phase of a target.
Most of current research or application of distributed scatterer filtering technology is to filter and optimize a center point based on a homogeneous region in a single estimation window, and most of the research or application focuses on a method for extracting and decomposing sample coherence matrix information. However, in a complex ground object type distribution scene, a single estimation window is not favorable for maintaining image details, and a sample coherence matrix obtained on the premise is a result under an average principle, and cannot accurately describe the scattering characteristics of the ground object corresponding to the central pixel.
Disclosure of Invention
In order to solve the technical problems, on the basis of the existing theory, the single estimation window is improved, the probability weight and the distance weight are introduced to obtain a weighted sample coherence matrix, and on the basis, advanced computer technology is utilized to realize the processing of homogeneous region identification, complex image filtering, triangular phase optimization and the like to filter the distributed scatterers.
The invention designs a distributed scatterer filtering method consisting 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:
step 2, pre-estimating a space-time baseline and coherence between the images, preprocessing and screening data, selecting a main image, and registering the data by taking the main image as a reference;
step 4, setting a determined sliding estimation window, and identifying pixel-by-pixel homogeneous regions of the radar image;
step 6, in the estimation window, carrying out adaptive spatial intensity filtering and phase filtering processing on the central pixel based on the homogeneous region and the PG filter to obtain the filtered intensity and phase;
step 7, triangular phase optimization is required after filtering;
step 8, setting a homogeneous pixel quantity threshold and a goodness-of-fit threshold, taking pixels meeting conditions as distributed scatterers, and updating the distributed scatterers by using optimized phases and filtering strength to obtain an updated radar image data set;
and 9, subsequently, performing joint calculation on the permanent scatterer and the distributed scatterer as high-coherence point targets on the updated radar image data set according to a PS-InSAR process.
Has the advantages that:
according to the invention, probability weight and distance weight are introduced into an adaptive spatial filtering method based on an estimation window, a weighted sample coherent matrix is calculated by combining scattering characteristics of distributed scatterers, and image filtering, triangular phase optimization and other processing are carried out based on a homogeneous region on the basis of the weighted sample coherent matrix, so that a technical process of filtering the distributed scatterers is formed, the filtered weighted sample coherent matrix can more accurately represent the scattering characteristics of a central pixel corresponding to a target ground object, the problem of loss of image details after filtering in a complex ground object type distribution scene is solved, and the information of PS points is kept as much as possible while identifying the distributed scatterers and carrying out filtering and noise reduction. 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 surface deformation is inverted in a non-urban area.
Drawings
FIG. 1(a) is a flow chart of a time series InSAR method for jointly resolving a permanent scatterer and a distributed scatterer according to the present invention;
fig. 1(b) is a schematic flow chart of a distributed scatterer filtering method of the present invention;
FIG. 2 is a pvalue neighborhood matrix and statistical homography mask: (a) supposing a probability matrix, (b) counting the homographic masks, and (c) communicating the homogeneous regions of the regions;
FIG. 3 is a schematic view of an investigation region;
FIG. 4(a) is a time series interferogram;
FIG. 4(b) is a time series correlation coefficient plot;
fig. 5 is a comparison between an original radar complex image and an interference fringe pattern and a coherence coefficient pattern generated by a model filtering result: (a) an original interference fringe pattern, (b) an original coherence coefficient pattern, (c) a filtered interference fringe pattern, and (d) a filtered coherence coefficient pattern;
fig. 6 is a comparison of measurement points of the original image dataset and the filtered image dataset inversion: (a) measuring points of inversion of the original image data set; (b) filtering the measurement points of the image dataset inversion;
fig. 7 is a comparison of residual terrain phase errors inverted from an original image dataset and a filtered image dataset: (a) residual terrain phase errors inverted from the original image data set, (b) residual terrain phase errors inverted from the filtered image data set;
FIG. 8 is a comparison of atmospheric and orbital phase errors for the inversion of an original image dataset and a filtered image dataset: (a) atmospheric and orbital phase errors inverted from the original image dataset, and (b) atmospheric and orbital phase errors inverted from the filtered image dataset.
Detailed Description
The technical solutions in the embodiments of the present invention will be described clearly and completely with reference to the accompanying drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, rather than all embodiments, and all other embodiments obtained by a person skilled in the art based on the embodiments of the present invention belong to the protection scope of the present invention without creative efforts.
On a long-time radar complex image set, a point target which continuously keeps high coherence is a stable point target. Based on the stable point target, a deformation model is established, the interference of space-time incoherent and atmospheric effects can be overcome, and more accurate terrain residual error and deformation phase are obtained through inversion.
As shown in fig. 1(a), distributed scatterer filtering is typically performed as a pre-processing part of a time-series InSAR process that jointly solves for persistent scatterers and distributed scatterers. In the area with low distribution density of the permanent scatterers, after the distributed scatterers are filtered based on the radar complex image sequence, the image noise is reduced, the coherence is improved, and then the permanent scatterers and the distributed scatterers are used as high-coherence point targets to be jointly solved.
As shown in fig. 1(b), a distributed scatterer filtering method of the present invention includes the following steps:
And 2, pre-estimating a space-time baseline and coherence between the images, eliminating the images with poor imaging quality and the image pairs with poor interference effect, selecting the most appropriate main image, and carrying out data registration by taking the main image as a reference.
And 3, sequencing the registered radar complex images according to imaging time, so as to facilitate subsequent processing.
And 4, setting a determined sliding estimation window, and identifying pixel-by-pixel homogeneous regions of the radar image. The specific method comprises the following steps: and in each sliding window, judging whether the time sequence amplitude vector corresponding to the time sequence complex vector corresponding to each pixel and the central pixel in the sliding window is from a certain same distribution population or not by a double-sample hypothesis test method, and if the hypothesis cannot be overturned, taking the pixel as the statistical same distribution pixel of the central pixel. After a statistical homographic pixel mask of a central pixel of a sliding window is obtained by using a double-sample detection method, connected region identification needs to be carried out on a binary mask image. The identified Pixel set consistent with the central Pixel label is the homogeneous region of the central Pixel, and the number of pixels consistent with the central Pixel label is the number of homogeneous pixels (SHP) of the central Pixel.
wherein x isijRepresents the pixel amplitude/phase values 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 formula for P (i, j) is as follows:
wherein p isijThe hypothesis that represents the pixel of (i, j) versus the center pixel examines the p-value, a represents a given level of significance, β is an exponential coefficient, the larger the coefficient, the greater the effect of the hypothesis that examines the p-value on the filter, and vice versa. G (i, j) is a two-dimensional Gaussian kernel, a computing publicThe formula is as follows:
wherein, σ represents a Gaussian standard deviation which reflects the dispersion degree of the weight data in the template, and the smaller σ is, the more concentrated the distribution is, the larger the template center coefficient is, the smaller the peripheral coefficient is, and the smaller the smoothness degree is; the larger σ, the closer the template center coefficient is to the perimeter coefficient, and the gaussian kernel approximates the mean template when σ → + ∞.
And 6, in the estimation window, carrying out self-adaptive spatial intensity filtering processing on the central pixel based on the homogeneous region and the PG filter to obtain the intensity after filtering.
And 7, in the estimation window, performing phase optimization on the central 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; based on the sample covariance matrix, normalization processing is carried out on all pixel amplitude values in the homogeneous region of each time phase, and a sample coherence matrix of the central pixel is estimated. And optimizing the time sequence phase vector of the central pixel according to the coherence coefficient and the interference phase in the sample coherence matrix to obtain an optimized phase. And simultaneously comparing the optimized interference phase vector with the original interference phase vector to evaluate the quality of the goodness of fit (GOF).
And 8, according to the set homogeneous pixel number threshold and the fitting goodness threshold, taking the pixels meeting the conditions as the distributed scatterers, and updating the distributed scatterers by using the optimized phases and the filtering strength. And obtaining an updated radar image data set.
And 9, subsequently, performing joint calculation on the permanent scatterer and the distributed scatterer as high-coherence point targets on the updated radar image data set according to a PS-InSAR process.
According to the embodiment of the invention, the filtering part of the distributed scatterer in the invention is realized by adopting a python programming language, and the realization process can be summarized as the following steps:
1. setting the size of an estimation window, assuming parameters such as a test probability value threshold, a homogeneous pixel quantity threshold, a goodness-of-fit index threshold, a significance level alpha, a test 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 the intensity filtering and the phase optimization pixel by pixel;
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 in the current sliding window, sequentially carrying out Anderson-Darling test on the central pixel and the adjacent pixels to obtain a sliding window hypothesis probability (p-value) matrix, as shown in figure 2(a), and obtaining a statistical homodistribution mask of the central pixel according to a given significance level alpha, as shown in figure 2 (b). And identifying a connected region of the statistical same-distribution pixel mask image. The identified pixel set consistent with the central pixel label is the homogeneous region of the central pixel, as shown in fig. 2(c), and the number of pixels consistent with the central pixel label is the homogeneous pixel number of the central pixel. The specific method of the Anderson-Darling test comprises the following steps: and judging whether the time sequence amplitude vectors X and Y corresponding to the time sequence complex number vector corresponding to each pixel and the central pixel in the sliding window come from a certain same distribution population, if the assumption cannot be overturned, the pixel is a statistical same distribution pixel of the central pixel. The inspection steps are as follows:
a) combining X and Y and sequencing from small to large to obtain a combined sample Z ═ Z1，z2，...，zk]，k＝2N；
b) With the double weighted sum of the squared differences of the single sample and the combined sample as the statistic, the calculation formula is as follows:
wherein M isiIs the set { x: x belongs to X and X is less than or equal to ziThe number of elements of.
c) Under the given significance level alpha, looking up the AD double-sample critical value table corresponding to the small sample, if AD > ADα(N, N), then H is rejected0I.e. the two samples are considered to be from the ensemble of different distributions.
5. Combining the hypothesis probability, the gaussian weight, and the weight coefficient and the significance level α corresponding to the gaussian weight and the gaussian weight, calculating a weight mask for filtering at each position of the sliding window, normalizing the weights, and defining a gaussian filter (hereinafter referred to as a PG filter) based on a hypothesis test p value, wherein the specific method is as follows: for a sliding estimation window of (2m +1) × (2n +1), the filter calculation formula for the center pixel (coordinates set to (0, 0)) is as follows:
wherein x isijRepresents the pixel amplitude/phase values 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 formula for P (i, j) is as follows:
wherein p isijThe hypothesis that represents the pixel of (i, j) versus the center pixel examines the p-value, a represents a given level of significance, β is an exponential coefficient, the larger the coefficient, the greater the effect of the hypothesis that examines the 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, σ represents a Gaussian standard deviation which reflects the dispersion degree of the weight data in the template, and the smaller σ is, the more concentrated the distribution is, the larger the template center coefficient is, the smaller the peripheral coefficient is, and the smaller the smoothness degree is; the larger σ, the closer the template center coefficient is to the perimeter coefficient, and the gaussian kernel approximates the mean template when σ → + ∞.
6. In the estimation window, carrying out self-adaptive spatial intensity filtering processing on the central pixel based on the homogeneous region and the PG filter to obtain the intensity after filtering;
7. within the estimation window, the center pixel is phase optimized based on the homogenous 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; based on the sample covariance matrix, normalization processing is carried out on all pixel amplitude values in the homogeneous region of each time phase, and a sample coherence matrix of the central pixel is estimated. According to a coherence coefficient and an interference phase in a sample coherence matrix, setting a radar image phase value with the first imaging time to be 0, and according to N (N-1)/2 interference phases in a weighted sample coherence matrix, solving a maximum likelihood estimation value of the N-1 time sequence interference phases by adopting a quasi-Newton method BFGS algorithm in an optimization algorithm so as to enable the maximum likelihood estimation value to meet a time sequence triangle principle. Meanwhile, because the winding phase value is [ - Π, Π ], in the optimization process, boundary setting is required; obtaining the optimized phase of the time sequence phase vector of the central pixel;
8. calculating optimal goodness of fit according to the interference phases before and after optimization; goodness of fit is expressed as the average of the real part of the interfering phase difference:
wherein N represents the number of time-series images,
representing the interference phase, θ, of the nth and the kth imagesn-θkRepresenting the interference phase of the optimized nth image and the optimized kth image, and Re representing the real part of the complex number. As can be seen from the above formula, the goodness-of-fit results are [ -1, 1 [)]The closer the result is to 1, the better the fitting effect, and the closer to-1, the worse the fitting effect.
9. And replacing the original complex image value of the position corresponding to the pixel meeting the requirement by the filtered intensity value and phase value according to the set homogeneous region pixel quantity threshold and the optimal fitting goodness threshold to obtain the updated long-time sequence radar complex image set.
According to the embodiment of the invention, a certain mining area in north Hebei Chengdu area of China is used as an application and verification case area of the invention. 115 degrees 54 '-119 degrees 15' E, 40 degrees 12 '-42 degrees 37' N at Chengdu city of Hebei province, east-west span of about 280km, south-north span of about 270km, total area of about 39500km2. The topography of Chengde city is reduced from northwest to southeast, the north part is inner Mongolian plateau, the middle part is distributed in shallow mountains, the south part belongs to Yanshan mountain range, and the middle-low mountain, inter-mountain basin and river valley are taken as main parts, and the forest coverage rate reaches 55.8%. The Chengdi city is bounded by a Kangbao-Binfield fracture zone, and spans two I-level structural units of a middle court ground platform and an inner Mongolian-Daxing Anlingland trough fold zone. The stratum in the area is developed completely, the mineralization condition is good, and the mineral resources are rich. The elevation ranges of the selected study area are 398-596m, average elevation 468m, maximum slope 88 degrees, and average slope 25 degrees. As shown in FIG. 3, the long-time-series radar image set for the study area schematic case area is derived from high-resolution Sentinel-1 radar image data (https:// scihub. copernius. eu/dhus/#/home) of European Space Agency (ESA) and its precise orbit data (https:// qc. Sentinel1.eo. ESA. int/aux _ poeorb /). The initial data set is 42 scenes in the Sentinel-1A ascending rail IW mode image with the imaging time between 2018 and 10 months 2020, and 25 scenes in the image which meets the requirements are reserved after space-time baseline estimation. An auxiliary DEM for removing terrain phase is derived from the space Shuttle Radar terrain mapping Mission (SRTM, https:// earth x plorer. usgs. gov /) resolution of about 30 m.
The distributed scatterer filtering is usually used as a preprocessing part in the process of time sequence InSAR for jointly resolving a permanent scatterer and a 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 measurement points of the time sequence InSAR technology, establishes a deformation model, separates terrain residual errors, atmospheric effect phases and deformation phases in phase information, and obtains terrain residual errors and linear deformation through inversion.
The distributed scatterer filtering technology is operated to obtain an updated long-time sequence radar image data set, and in order to better show the effects of the technology in reducing noise, improving coherence and keeping details, the section mainly shows the effects of a filtered radar complex image intensity image and an interferogram and a coherence coefficient image generated by interference of the filtered radar complex image intensity image in time and space angles. In a time angle, a time series interferogram and a time series coherence coefficient map are shown; in a spatial perspective, a time-series mean intensity map and a time-series mean coherence coefficient map are shown. The model operation results show a time-series interferogram as shown in fig. 4(a), and a time-series coherence coefficient map of fig. 4 (b);
the invention compares and verifies an interference fringe pattern and a coherence coefficient pattern which are generated by an original radar complex image and a filtering radar complex image, and discusses the effectiveness and stability of a distributed scatterer filtering model for extracting the distributed scatterers as shown in fig. 5(a) - (d) and table 1. In fig. 5(a) - (d), the interference fringes reflect the distribution of the interference phase without unwrapping, the coherence factor reflects the quality of the interference phase, and the visual effect can be seen as the effect of noise reduction and coherence improvement of the filtering. In order to evaluate the distortion condition of the filtered image, 5 image quality evaluation indexes of a mean value, a standard deviation, a signal-to-noise ratio, an equivalent visual number and an edge retention coefficient are selected to compare the image quality of the interference fringe pattern before and after filtering, as shown in table 1, it can be seen that the difference between the mean value and the variance of the original interference fringe pattern and the filtered interference fringe pattern is 0.0138 and-0.0305 respectively, which indicates that the statistical characteristics of the images before and after filtering are not changed greatly, the equivalent visual number is 4.65E-04 and is 9.07E-04 respectively, which indicates that a certain noise reduction effect exists after filtering, and the signal-to-noise ratio and the edge retention coefficient of the filtered interference fringe pattern compared with the original interference fringe pattern are 0.9871 and 0.8165 respectively, which indicates that the distortion of the image is small and the edge retention is good. Fig. 5 is a comparison between the original radar complex image and an interference fringe pattern and a coherence coefficient pattern generated by a model filtering result.
TABLE 1 comparison of interference fringe pattern image quality evaluation indexes generated by original radar complex image and model filtering result
In the case, surface deformation inversion is respectively carried out on the long-time sequence radar image data sets before and after the distributed scatterer filtering, the density of a measuring point, the residual topographic phase error, the atmospheric error phase and the orbit error phase are compared and verified, and the improvement condition of the distributed scatterer technology on the density and the quality of the measuring point in the time sequence InSAR technology is discussed.
Fig. 6 shows the measurement point distribution of the inversion of the original image dataset and the filtered image dataset, and fig. 6(a) - (b) show the measurement point distribution of the inversion of the earth surface deformation of the long-time sequence radar image dataset before and after the distributed scatterer filtering, and it can be seen that the distribution of the earth surface deformation rates obtained by the inversion is more consistent and is-14.92-12.23 mm/y and-12.36-11.87 mm/y, and the distribution also has stronger consistency in the 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 filtering can significantly improve the high coherent point density.
TABLE 2 measurement point comparison of original image dataset and filtered image dataset inversion
FIG. 7 is a comparison of residual terrain phase errors inverted from an original image dataset and a filtered image dataset; fig. 7(a) - (b) are residual terrain phase error distribution conditions of ground surface deformation inversion performed on long-time sequence radar image datasets before and after distributed scatterer filtering, respectively, and it can be seen that residual terrain phase errors obtained by inversion have strong consistency in spatial distribution. Table 3 shows that, when the density of the measurement points obtained from the filtered image dataset is compared with the original image dataset, the maximum value and the mean value of the residual topographic phase error are both reduced, and it can be seen that the inversion accuracy can be improved by the distributed scatterer filtering.
TABLE 3 residual topographic phase error contrast for original image dataset and filtered image dataset inversion
Fig. 8(a) - (b) show the atmospheric and orbital phase error distributions of the long-time sequence radar image datasets before and after the distributed scatterer filtering, respectively, where earth surface deformation inversion is performed, and it can be seen that the atmospheric and orbital phase errors obtained by inversion have strong consistency in both time distribution and spatial distribution. Compared with the original image data set, the error range of the atmospheric and orbit phase errors obtained by the filtering image data set is reduced to some extent, and the inversion accuracy of the atmospheric and orbit phase errors by the filtering of the distributed scatterers is improved to some extent.
Although illustrative embodiments of the present invention have been described above 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 various changes may be apparent to those skilled in the art, and it is intended that all inventive concepts utilizing the inventive concepts set forth herein be protected without departing from the spirit and scope of the present invention as defined and limited by the appended claims.
Claims (4)
1. A method of filtering distributed scatterers, comprising the steps of:
step 1, taking a long-time sequence radar complex image set, orbit data and a research area reference DEM as input;
step 2, pre-estimating a space-time baseline and coherence between the images, preprocessing and screening data, selecting a main image, and registering the data by taking the main image as a reference;
step 3, sequencing the registered radar complex images according to imaging time, so as to facilitate subsequent processing;
step 4, setting a determined sliding estimation window, and identifying pixel-by-pixel homogeneous regions of the radar image;
step 5, respectively defining probability and distance weight according to the same distribution probability and distance, and calculating a weighted sample coherent matrix in a homogeneous region as an expression of scattering characteristics of a target point;
step 6, in the estimation window, carrying out adaptive spatial intensity filtering and phase filtering processing on the central pixel based on the homogeneous region and the PG filter to obtain the filtered intensity and phase;
step 7, triangular phase optimization is required after filtering;
step 8, setting a homogeneous pixel quantity threshold and a goodness-of-fit threshold, taking pixels meeting conditions as distributed scatterers, and updating the distributed scatterers by using optimized phases and filtering strength to obtain an updated radar image data set;
and 9, subsequently, performing joint calculation on the permanent scatterer and the distributed scatterer as high-coherence point targets on the updated radar image data set according to a PS-InSAR process.
2. The method for filtering distributed scatterers according to claim 1, wherein in step 4, a determined sliding estimation window is set, and pixel-by-pixel homogeneous region identification is performed on the radar image by the specific method:
in each sliding window, judging whether a time sequence amplitude vector corresponding to a time sequence complex vector corresponding to each pixel and a central pixel in the sliding window is from a certain same distribution population or not by a double-sample hypothesis testing method, if the hypothesis cannot be overturned, taking the pixel as a statistical same distribution pixel of the central pixel, and after obtaining a statistical same distribution pixel mask of the central pixel in the sliding window by using the double-sample hypothesis testing method, identifying a connected region of a binary mask image; the identified pixel set consistent with the central pixel label is the homogeneous area of the central pixel, and the number of the pixels consistent with the central pixel label is the number of the homogeneous pixels of the central pixel.
3. The method for filtering distributed scatterers according to claim 1, wherein in the step 5, probability and distance weights are defined by the same distribution probability and distance, respectively, and a weighted sample coherence matrix in a homogeneous region is calculated as an expression of scattering characteristics of a target point, specifically as follows:
defining a gaussian filter based on hypothesis testing p-values, for a sliding estimation window of (2m +1) × (2n +1), m, n being a natural number with its central pixel coordinate set to (0, 0), the filter calculation formula for the central pixel is as follows:
wherein x isijThe pixel amplitude/phase values representing the (i, j) locations, P (i, j) representing the weights corresponding to the hypothetical test P values, G (i, j) representing the Gaussian kernel weights, and P (i, j) being calculated as follows:
wherein p isijA hypothetical test p-value representing the pixel at the (i, j) position and the center pixel, a representing a given level of significance, β being an exponential coefficient, the larger the coefficient, the greater the hypothetical test p-value's impact on the filter, and vice versa, G (i, j) being a two-dimensional gaussian kernel, the calculation being as follows:
wherein, σ represents a Gaussian standard deviation which reflects the dispersion degree of the weight data in the template, and the smaller σ is, the more concentrated the distribution is, the larger the template center coefficient is, the smaller the peripheral coefficient is, and the smaller the smoothness degree is; the larger σ, the closer the template center coefficient is to the perimeter coefficient, and the gaussian kernel approximates the mean template when σ → + ∞.
4. The method for filtering distributed scatterers according to claim 1, wherein after the filtering in step 7, triangular phase optimization is required, specifically comprising:
in the estimation window, phase optimization is performed on the central pixel based on the homogeneous region and the PG filter, and 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; based on the sample covariance matrix, normalizing all pixel amplitude values in the homogeneous region of each time phase, estimating a sample coherent matrix of a central pixel, optimizing a time sequence phase vector of the central pixel according to a coherent coefficient and an interference phase in the sample coherent matrix to obtain an optimized phase, and simultaneously comparing the optimized interference phase vector with an original interference phase vector to evaluate the quality of the fitting goodness.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111675777.4A CN114325706A (en) | 2021-12-31 | 2021-12-31 | Distributed scatterer filtering method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111675777.4A CN114325706A (en) | 2021-12-31 | 2021-12-31 | Distributed scatterer filtering method |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114325706A true CN114325706A (en) | 2022-04-12 |
Family
ID=81023352
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111675777.4A Pending CN114325706A (en) | 2021-12-31 | 2021-12-31 | Distributed scatterer filtering method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114325706A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115143877A (en) * | 2022-05-30 | 2022-10-04 | 中南大学 | SAR offset tracking method, device, equipment and medium based on strong scattering amplitude suppression and abnormal value identification |
| CN116047519A (en) * | 2023-03-30 | 2023-05-02 | 山东建筑大学 | Point selection method based on synthetic aperture radar interferometry technology |
-
2021
- 2021-12-31 CN CN202111675777.4A patent/CN114325706A/en active Pending
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115143877A (en) * | 2022-05-30 | 2022-10-04 | 中南大学 | SAR offset tracking method, device, equipment and medium based on strong scattering amplitude suppression and abnormal value identification |
| 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 |
| CN116047519A (en) * | 2023-03-30 | 2023-05-02 | 山东建筑大学 | Point selection method based on synthetic aperture radar interferometry technology |
| CN116047519B (en) * | 2023-03-30 | 2023-06-16 | 山东建筑大学 | Point selection method based on synthetic aperture radar interferometry technology |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Dai et al. | Coastline extraction from repeat high resolution satellite imagery | |
| CN114325706A (en) | Distributed scatterer filtering method | |
| CN109388887B (en) | Quantitative analysis method and system for ground settlement influence factors | |
| CN113960595A (en) | Surface deformation monitoring method and system | |
| CN109212522B (en) | High-precision DEM inversion method and device based on double-base satellite-borne SAR | |
| Wang et al. | Construction and optimization method of the open-pit mine DEM based on the oblique photogrammetry generated DSM | |
| CN113866764A (en) | Landslide susceptibility improvement evaluation method based on InSAR and LR-IOE model | |
| Ma et al. | Minimum spanning tree co-registration approach for time-series Sentinel-1 TOPS data | |
| Jiang et al. | Application of multitemporal InSAR covariance and information fusion to robust road extraction | |
| CN111856459B (en) | Improved DEM maximum likelihood constraint multi-baseline InSAR phase unwrapping method | |
| CN103871039A (en) | Generation method for difference chart in SAR (Synthetic Aperture Radar) image change detection | |
| 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 | |
| CN116047519A (en) | Point selection method based on synthetic aperture radar interferometry technology | |
| Khoshboresh-Masouleh et al. | A deep learning method for near-real-time cloud and cloud shadow segmentation from gaofen-1 images | |
| Tao et al. | Simulation-based building change detection from multiangle SAR images and digital surface models | |
| CN115079172A (en) | MTInSAR landslide monitoring method, equipment and storage medium | |
| Thissen | Automating Surface Water Detection for Rivers | |
| CN113687353A (en) | DS target phase optimization method based on homogeneous pixel time sequence phase matrix decomposition | |
| Sui et al. | Processing of multitemporal data and change detection | |
| Shen et al. | PolInSAR complex coherence nonlocal estimation using shape-adaptive patches matching and trace-moment-based NLRB estimator | |
| Luo et al. | Ice flow velocity mapping in East Antarctica using historical images from 1960s to 1980s: recent progress | |
| CN115540908A (en) | InSAR interference fringe matching method based on wavelet transformation | |
| Wang et al. | An improved SAR radiometric terrain correction method and its application in polarimetric SAR terrain effect reduction | |
| Luo et al. | Despeckling multi-temporal polarimetric SAR data based on tensor decomposition |
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 |