CN113281749B - Timing sequence InSAR high coherence point selection method considering homogeneity - Google Patents

Timing sequence InSAR high coherence point selection method considering homogeneity Download PDF

Info

Publication number
CN113281749B
CN113281749B CN202110614249.1A CN202110614249A CN113281749B CN 113281749 B CN113281749 B CN 113281749B CN 202110614249 A CN202110614249 A CN 202110614249A CN 113281749 B CN113281749 B CN 113281749B
Authority
CN
China
Prior art keywords
coherence
coefficient
image
points
pixel
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
Application number
CN202110614249.1A
Other languages
Chinese (zh)
Other versions
CN113281749A (en
Inventor
高贵
张涛
高昇
曹敏
尹文禄
刘伟
刘涛
刘宏立
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest Jiaotong University
Original Assignee
Southwest Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southwest Jiaotong University filed Critical Southwest Jiaotong University
Priority to CN202110614249.1A priority Critical patent/CN113281749B/en
Publication of CN113281749A publication Critical patent/CN113281749A/en
Application granted granted Critical
Publication of CN113281749B publication Critical patent/CN113281749B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/885—Radar or analogous systems specially adapted for specific applications for ground probing

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Image Processing (AREA)

Abstract

The invention discloses a timing sequence InSAR high coherence point selection method considering homogeneity, which comprises the following steps: s1, performing interference processing on a sequential SAR image to obtain a sequential coherence coefficient image; s2, acquiring clustering areas of low, medium and high pixel points of a time sequence coherence coefficient image and Gaussian parameters corresponding to each type, and removing the pixel points of the low coherence area according to the Gaussian parameters; s3, overlapping the PS points and the clustered coherent coefficient images, calculating the coherent coefficient mean value of all the PS points and the mean value of the residual medium and high coherent coefficient pixel points in the coherent coefficient images, and obtaining a coherent coefficient threshold; s4, performing secondary screening on residual pixel points with medium and high coherence coefficients through a coherence coefficient threshold value, interpolating a coherence coefficient average value in a filtering window size for a pixel missing in a result image, and taking a final result as a high coherence point of SBAS-InSAR. The invention reserves the integrity of the high coherence area and the high coherence point of the high coherence area which is partially lost due to interference factors to the maximum degree on the basis of considering the coherence coefficient pixels with the same kind of property.

Description

Timing sequence InSAR high coherence point selection method considering homogeneity
Technical Field
The invention relates to the technical field of geological data processing, in particular to a time sequence InSAR high coherence point selection method considering homogeneity.
Background
The urban land subsidence is a geological disaster which seriously threatens urban safety, and causes huge damage to the production and living of people and the sustainable development of towns. Therefore, the method and the device can be used for effectively monitoring the deformation of the urban area which is already or is in subsidence for a long time, and timely finding out the subsidence area and the deformation development trend thereof, so that the method and the device are necessary conditions for realizing the early maintenance of the subsidence area and further changing the 'passive disaster avoidance and relief' into the active 'disaster prevention and treatment'. The conventional earth surface deformation monitoring method is mostly focused on the ground monitoring technology and the underground monitoring technology, such as a geodetic method, a GPS method, drilling point position monitoring, geophysical detection technology and the like. The traditional deformation monitoring method is difficult to realize large-range monitoring point layout, so that the limit of intersection is generated in the monitoring space. In recent years, along with rapid development of communication technology and sensor technology, the interference synthetic aperture radar (Interferometric Synthetic Aperture Radar, inSAR) acquires large-scale earth surface deformation information in all-weather, all-day, cloud penetrating, mist penetrating, periodic and high-resolution modes, so that great advantages are brought into play in urban ground subsidence monitoring, particularly the capability of earth surface micro deformation detection is greatly improved due to the occurrence of Differential radar interference SAR (D-InSAR), and the D-InSAR can realize earth surface deformation measurement in the centimeter or even millimeter level under the condition that ground control points are not needed. However, interference factors such as interference decorrelation, inaccurate atmospheric delay phase estimation, DEM (DigitalElevationModel) error and the like reduce the measurement accuracy of the D-InSAR, meanwhile, the D-InSAR cannot acquire time sequence deformation information of a monitoring area, and certain limitation is generated in the aspects of disaster prediction and evaluation of a ground subsidence area. Thus, on the basis of D-InSAR, a permanent scatterer (PermanentScatterertechnology, PS), a small baseline set (SBAS), staMPS (Stanfordmethodof permanent scatterer technology), etc. are sequentially proposed, collectively referred to as time series InSAR.
The SBAS-InSAR technology is proposed by the Italy student Berardino team, and the influence of time-space decorrelation is effectively reduced on the basis of controlling a time-space baseline. The technology combines a stable scatterer to separate an atmospheric interference phase in an interference phase, so as to obtain a time sequence deformation result of a research area. The core of the technology for realizing the surface deformation monitoring is to carry out statistical analysis and modeling on the interference phase of a scatterer with stable scattering characteristics in a time sequence image, realize effective estimation and separation of each component in the interference phase, carry out fitting modeling on each component on the basis, acquire a corresponding resolving model, and finally acquire the absolute deformation phase in the interference phase of a research area under the condition of minimizing the interference of other phases. Representative methods include: the method is characterized in that a coherence coefficient threshold value method is adopted to obtain high coherence points, however, the method generally only considers the selection of the high coherence points according to the coherence coefficient, and under the condition that interference factors such as a topography residual error phase, an atmosphere phase, an elevation error phase and the like are not removed, the number and the reliability of the coherence points obtained according to the coherence coefficient threshold value reduce the inversion precision of each component in the interference phase to a certain extent. In addition, from the viewpoint of a calculation method of the coherence coefficient, the size of a calculation window has a larger influence on the selection of the coherence points, and the problem of isolation and effective high coherence point omission is easily caused by the overlarge window; conversely, too small a window is prone to false selection of unstable targets near the high coherence point. The amplitude deviation index threshold method is used for obtaining the high coherence point, the correlation between the amplitude information and the interference relative is basically irrelevant, and the interference of the inter-SAR image uncorrelation or low correlation is overcome to a certain extent. However, the amplitude information is closely related to the image resolution, resulting in a limitation of narrow data sources; the method is time-consuming in calculation, and the selected high coherence points cannot usually consider the edge and internal characteristics of the high coherence region, so that the loss of the coherence points is easy to cause; in areas with fewer artificial buildings, the reliability of high coherence points obtained by adopting amplitude information is not high. By adopting the double-threshold detection technology, the method has larger requirements on the density of the stable scatterer, improves the reliability of high coherence points, reduces the number of the coherence points, and has limited application range.
Disclosure of Invention
The invention aims to provide a time sequence InSAR high coherence point selection method considering homogeneity, so as to overcome the defects in the prior art.
In order to achieve the above purpose, the technical scheme adopted by the invention is as follows:
a time sequence InSAR high coherence point selection method considering homogeneity comprises the following steps:
s1, performing interference processing on a sequential SAR image to obtain a sequential coherence coefficient image;
s2, taking the time sequence coherence coefficient image as a sample input of a GMM-EM algorithm, acquiring clustering areas of low, medium and high pixel points of the time sequence coherence coefficient image and Gaussian parameters corresponding to each type, and removing the pixel points of the low coherence area according to the Gaussian parameters;
s3, overlapping the PS points and the clustered coherent coefficient images, calculating the coherent coefficient mean value of all the PS points and the mean value of the residual medium and high coherent coefficient pixel points in the coherent coefficient images, and acquiring a coherent coefficient threshold according to the mean value and the mean value;
s4, performing secondary screening on residual pixel points with medium and high coherence coefficients through a coherence coefficient threshold value, interpolating a coherence coefficient average value in a filtering window size for a pixel missing in a result image, and taking a final result as a high coherence point of SBAS-InSAR.
Further, the time sequence SAR image in the step S1 is a C-band high-resolution time sequence SAR image obtained by a satellite; the interference treatment comprises at least a de-flattening, a filtering and an unwrapping.
Further, the specific steps of acquiring the clustering areas of the low, medium and high three types of pixel points of the time sequence coherence coefficient image in the step S2 are as follows:
coherence coefficient image as difference image gamma of same geographic position in different time phase SAR images d The expression is:
Figure BDA0003097385400000031
from equation (1) above, the time series of coherence coefficients of any one resolution element of the N interference pairs can be found: gamma ray 1 、γ 2 、γ 3 、...、γ N On the basis, the average value of the coherence coefficient of each pixel is calculated:
Figure BDA0003097385400000032
will be
Figure BDA0003097385400000033
As a difference imageWherein the picture elements are associated with a probability density distribution +.>
Figure BDA0003097385400000034
Can be distributed by high coherence coefficient type pixels>
Figure BDA0003097385400000035
Middle coherence coefficient pixel distribution->
Figure BDA0003097385400000036
And low coherence coefficient pixel-like distribution->
Figure BDA0003097385400000037
The composition is as follows:
Figure BDA0003097385400000038
wherein:
Figure BDA0003097385400000039
a coherence coefficient value representing a pixel in the difference image; omega 0 、ω 1 And omega 2 The prior probabilities of pixels of the high, medium and low coherence class are represented respectively, < >>
Figure BDA00030973854000000310
Linear combination of N components in the coherence coefficient image, namely:
Figure BDA00030973854000000311
wherein:
Figure BDA00030973854000000312
is the probability density distribution of the nth component, p (n) is the prior probability of a certain mixed component in the GMM, and the probability density function of each component obeys the gaussian distribution, namely:
Figure BDA00030973854000000313
wherein: mu (mu) n And
Figure BDA00030973854000000314
mean and variance of the nth component are represented, respectively; the gaussian distribution function satisfies:
Figure BDA00030973854000000315
the prior probability of the mixed components satisfies
Figure BDA00030973854000000316
And p (n) is more than or equal to 0 and less than or equal to 1;
as can be seen from formulas (4) and (5), the GMM has a total of N components, including each component
Figure BDA00030973854000000317
And a priori probabilities p (N) (n=1, 2,3.., N) are needed to be solved, i.e., at GMM can be determined by the following parameters:
Figure BDA0003097385400000041
further, the step S2 of acquiring gaussian parameters corresponding to the clustering areas of the low, medium and high pixel points of the time sequence coherence coefficient image specifically includes the steps of:
parameter initialization: clustering samples by using a K-means clustering algorithm, and taking the mean value of each type as mu 0 And calculate
Figure BDA0003097385400000042
Taking the proportion of various samples to the total number of samples as initial posterior probability;
an estimation step: data
Figure BDA0003097385400000043
As a sampleThis data, complete data->
Figure BDA0003097385400000044
From the variable z= { Z 1 ，z 2 ，z 3 ，...，z N Estimate of }, where z n And difference image->
Figure BDA0003097385400000045
With the same dimensions, finish data Y d Is:
Figure BDA0003097385400000046
wherein z is n (i, j) is a posterior probability,
Figure BDA0003097385400000047
t represents the number of iterations, parameter +.>
Figure BDA0003097385400000048
Is the estimated parameter of the t-th iteration;
maximizing: parameters used in the next iteration in M steps
Figure BDA0003097385400000049
From the variable z n And (i, j) estimating, and the like, and finally obtaining various parameter values of the Gaussian mixture model, namely:
Figure BDA00030973854000000410
/>
Figure BDA00030973854000000411
Figure BDA00030973854000000412
further, the average value of the remaining medium-high coherence coefficient pixel points in the coherence coefficient image in the step S4 is
Figure BDA00030973854000000413
All PS points coherence coefficient mean value is +.>
Figure BDA00030973854000000414
And calculating the ratio of each pel point to the total pel point number greater than the average value by the following formula:
Figure BDA00030973854000000415
Figure BDA0003097385400000051
the calculation formula of the coherence coefficient threshold value D is as follows:
Figure BDA0003097385400000052
compared with the prior art, the invention has the advantages that: according to the invention, the GMM algorithm is used for carrying out cluster analysis on the coherence coefficient of the ground object in the research area, so that the distribution condition of three kinds of coherence coefficients of high, medium and low is obtained, then the pixel points belonging to the same kind of property are better aggregated together, the pixel points with low coherence coefficient are removed according to three kinds of Gaussian parameters, the integrity of the high coherence area and the high coherence point, which is partially lost due to interference factors, are reserved to the maximum extent on the basis of considering the pixels with the coherence coefficient with the same kind of property, and the influence of the missed selection and the wrong selection of the high coherence point, which is caused by the overlarge or the overlarge calculation window, can be reduced to a certain extent.
Drawings
In order to more clearly illustrate the embodiments of the invention or the technical solutions in the prior art, the drawings that are required in the embodiments or the description of the prior art will be briefly described, it being obvious that the drawings in the following description are only some embodiments of the invention, and that other drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of a method for selecting a high coherence point of a time sequence InSAR in consideration of homogeneity.
FIG. 2 is a diagram of interference versus combination in the present invention.
Fig. 3 is a gaussian mixture distribution histogram of the coherence coefficient in the present invention.
Fig. 4 is a plot of coherence coefficient and PS-point superposition profile in the present invention.
FIG. 5 is a graph showing the correlation statistics of the conventional method and the method of the present invention.
Detailed Description
The preferred embodiments of the present invention will be described in detail below with reference to the accompanying drawings so that the advantages and features of the present invention can be more easily understood by those skilled in the art, thereby making clear and defining the scope of the present invention.
Referring to fig. 1, in order to accurately obtain high coherence points in a coherent image, the defect that the edge information of a high coherence area cannot be reserved when a traditional high coherence point is obtained based on a coherence coefficient threshold value determined manually is overcome, the problem of the selection precision and number of high coherence scatterers with the same nature is mainly considered, a GMM algorithm is used for clustering the homogeneous points first, and then the threshold value is set by combining with the coherence coefficient of a PS point, so that the spatial distribution of the high coherence area is reserved to the greatest extent. The embodiment discloses a timing sequence InSAR high coherence point selection method considering homogeneity, which specifically comprises the following steps:
and S1, performing interference processing on the sequential SAR image to obtain a sequential coherence coefficient image.
In this embodiment, a high-resolution sequential SAR image of a C band is obtained through a European Sentinel-1A (Sentinel-1A) satellite, an optimal interference pair combination is obtained according to a preset time-space baseline threshold value, interference processing is performed, and a sequential coherence coefficient image is obtained through interference processing such as de-flattening, filtering and unwrapping. Wherein the interference pair combination result is shown in fig. 2.
S2, taking the time sequence coherence coefficient image as sample input of a GMM-EM algorithm, acquiring clustering areas of low, medium and high pixel points of the time sequence coherence coefficient image and Gaussian parameters corresponding to each type, and removing the pixel points of the low coherence area according to the Gaussian parameters.
The GMM is independent based on the pixel values in each interference image, and the coherence coefficient image is used as a difference image gamma of the same geographic position in different time-phase SAR images d The expression is:
Figure BDA0003097385400000061
the time series of coherence coefficients for any one of the resolution elements of the N interference pairs can be found according to equation (1): gamma ray 1 、γ 2 、γ 3 、...、γ N . On the basis, the average value of the coherence coefficient of each pixel is calculated:
Figure BDA0003097385400000062
will be
Figure BDA0003097385400000063
As a difference image, wherein the pixels are associated with a probability density distribution +.>
Figure BDA0003097385400000064
Can be distributed by high coherence coefficient type pixels>
Figure BDA0003097385400000065
Middle coherence coefficient pixel distribution->
Figure BDA0003097385400000066
And low coherence coefficient pixel-like distribution->
Figure BDA0003097385400000067
The composition is as follows:
Figure BDA0003097385400000068
wherein:
Figure BDA0003097385400000069
a coherence coefficient value representing a pixel in the difference image; omega 0 、ω 1 And omega 2 Representing the prior probabilities of the pixels of the high, medium and low coherence classes, respectively. />
Figure BDA00030973854000000610
Linear combination of N components in the coherence coefficient image, namely:
Figure BDA00030973854000000611
wherein:
Figure BDA00030973854000000612
is the probability density distribution of the nth component, p (n) is the prior probability of a certain mixed component in the GMM, and the probability density function of each component obeys the gaussian distribution, namely:
Figure BDA0003097385400000071
wherein: mu (mu) n And
Figure BDA0003097385400000072
mean and variance of the nth component are represented, respectively; the Gaussian distribution function satisfies:>
Figure BDA0003097385400000073
the prior probability of the mixed components satisfies
Figure BDA0003097385400000074
And p (n) is more than or equal to 0 and less than or equal to 1.
As can be seen from formulas (4) and (5), the GMM has a total of N components, including each component
Figure BDA0003097385400000075
And a priori probabilities p (N) (n=1, 2,3.., N) are needed to be solved, i.e., at GMM can be determined by the following parameters:
Figure BDA0003097385400000076
the embodiment adopts an EM algorithm to solve parameters in a Gaussian mixture model, and is an iterative approximation method. A set of parameters is typically randomly initialized
Figure BDA0003097385400000077
And calculate +.>
Figure BDA0003097385400000078
The posterior probability p (n) generated by a certain component in the GMM is then calculated by iterative optimal value calculation, and each iteration consists of two processes of an E step (obtaining the maximum expected value of the observed value) and an M step (obtaining the maximum value). Two steps->
Figure BDA0003097385400000079
Alternating until the local optimum is converged, wherein the local optimum is independent of each other, and the specific steps are as follows:
in the parameter initialization step, the embodiment uses a K-Means (K-Means) clustering algorithm to cluster samples, and uses various Means as mu 0 And calculate
Figure BDA00030973854000000710
The ratio of each type of sample to the total number of samples is taken as the initial posterior probability.
Estimating step, data
Figure BDA00030973854000000711
As sample data, complete data->
Figure BDA00030973854000000712
From the variable z= { Z 1 ，z 2 ，z 3 ，...，z N Estimate of }, where z n And difference image->
Figure BDA00030973854000000713
With the same dimensions, finish data Y d Is:
Figure BDA00030973854000000714
wherein z is n (i, j) is a posterior probability, i.e.:
Figure BDA0003097385400000081
wherein t represents the number of iterations, the parameter
Figure BDA0003097385400000082
Is the estimated parameter for the t-th iteration.
Maximizing: parameters used in the next iteration in M steps
Figure BDA0003097385400000083
Can be represented by the variable z n And (i, j) estimating, and the like, and finally obtaining various parameter values of the Gaussian mixture model, namely:
Figure BDA0003097385400000084
Figure BDA0003097385400000085
Figure BDA0003097385400000086
/>
presetting a component n=3 in the GMM in an EM algorithm; according to the number of sample points obtained from the average value of the coherence coefficients, 822867 pixel points are distributed in the range of 0.28-0.37, 742701 pixel points are distributed in the range of 0.37-0.61, 394230 pixel points are distributed in the range of 0.61-0.93, and the frequency histogram after mixing the three coherence coefficient intervals is shown in fig. 3 (a). The mean value of each Gaussian distribution is determined by the distribution range of each component coherence coefficient and is respectively as follows: 0.75, 0.49, 0.32; the probability that each sample point belongs to each cluster is set to 0.4, 0.2, respectively. The distribution of the three GMM blend components based on the mean value of the coherence coefficients can be seen from fig. 3 (b).
In the time sequence InSAR technology, the higher the pixel coherence coefficient is, the more stable the corresponding region is, and further, the more reliable the deformation phase obtained by inverting the stable region is, however, in the research region, vegetation, woodland, water body and the like are all low coherence regions, and the coherence coefficient is generally distributed between 0 and 0.3; the soil with bare leakage or the soil doped with rock is a medium coherent area, and the coherence coefficient is distributed between 0.3 and 0.55; the artificial buildings such as houses and roads are high coherence areas, and the coherence coefficient is distributed between 0.55 and 1. This is the basis for dividing the original coherence coefficient image into low, medium and high 3 types of pixels. In the embodiment, a GMM-EM algorithm is adopted to perform cluster analysis on pixel points in a research area where low-coherence ground objects such as vegetation, water bodies and the like are removed, pixel points with the same property (namely, pixel points with the coherence coefficient belonging to the same Gaussian distribution) are clustered, gaussian parameters of 3 types of pixel points are obtained according to the algorithm, and on the basis, the pixel points with the coherence coefficient are respectively low, medium and high; the method is a basis for classifying clustered coherence coefficient pixels into low, medium and high categories.
And S3, superposing the PS points and the clustered coherent coefficient images, calculating the coherent coefficient mean value of all the PS points and the mean value of the residual medium and high coherent coefficient pixel points in the coherent coefficient images, and acquiring a coherent coefficient threshold according to the mean value and the mean value.
The coherence coefficient image is obtained under the condition that the interference image still contains a plurality of interference phases, and deformation information of the inversion target area is unreliable to a certain extent by completely relying on high coherence points in the coherence coefficient image. PS points can show similar, relatively large intensity information in all time-series images, the phase dispersion is relatively small, the points do not need to analyze phase information, and the points also show high coherence coefficient characteristics after interference processing. Therefore, PS points are selected in consideration of whether the scattering characteristics of the pixels are stable, and are more reliable than coherent points selected directly through the coherence coefficient threshold.
In this embodiment, the coherence coefficient threshold is calculated in consideration of the homogeneous coherence points acquired by the PS point and the GMM. Firstly, pixels in a coherent coefficient map are classified into three categories according to the size of the coherent coefficient: the distribution condition of pixel points with the same kind of property in a research area is obtained through a GMM-EM algorithm, and low-coherence coefficient components are removed according to Gaussian parameters of the same mass points. Respectively calculating the average value of the coherence coefficients of the residual pixel points
Figure BDA0003097385400000091
And mean value of coherence coefficients of PS points +.>
Figure BDA0003097385400000092
And the pixel points are larger than the duty ratio of the average value in the total pixel point number.
Figure BDA0003097385400000093
Figure BDA0003097385400000094
Finally, selecting a threshold value D capable of keeping the number of high coherence points and the accuracy of the high coherence point to the maximum extent through the parameters, wherein a calculation formula is as follows:
Figure BDA0003097385400000095
and S4, performing secondary screening on residual medium and high coherence coefficient pixel points through a coherence coefficient threshold D, interpolating a coherence coefficient average value in a filtering window size for a missing pixel in a result image to recover a high coherence pixel missing due to other interference factors, and taking a final result as a high coherence point of the SBAS-InSAR. And solving linear deformation and elevation errors by utilizing SVD, and finally obtaining a time sequence deformation phase of a high coherence point through estimation and elimination of each component in the interference phase of the high coherence region.
In this embodiment, PS points are superimposed in a coherence coefficient map of a low coherence region, the superposition of partial regions in a research range is shown in fig. four, and the distribution of artificial buildings in the region can be seen from an amplitude average map in fig. 4 (a), in view of the fact that the research range is high in coherence points (artificial buildings) and dense in distribution, the coherence coefficient threshold value when PS points are selected is set to be 0.55, and then the coherence coefficient of PS points is distributed between 0.55 and 0.9. Fig. 4 (b) shows that PS dot positions are consistent with the distribution of highlight pels in the amplitude mean map. Compared with the numerical value in the coherence coefficient mean value graph, the coherence coefficient of the PS point is higher, so that the position of the PS point in the whole time sequence image keeps stable scattering characteristics, and the high coherence region selected by the method is more reliable. The average distribution diagram of the coherence coefficient is shown in fig. 4 (c), and the superposition situation of the PS point and the coherence coefficient is shown in fig. 4 (d), so that the pixel point with high coherence coefficient has consistency with the distribution of the high-brightness pixel point in the average distribution diagram of the PS point and the amplitude. And then, setting a coherence coefficient threshold value according to the distribution range of the coherence coefficient of the coincident region in the PS point and coherence coefficient mean value diagram, and selecting the final high-quality coherence point. The average value of the coherence coefficients of low coherence distribution in three components of the coherence coefficient Gaussian mixture model is 0.4, the lowest coherence coefficient of the PS point is 0.55, the calculated result according to the formula is taken as a threshold value, namely the coherence coefficient threshold value is set to be 0.4, and 518576 high coherence points are finally determined. Finally, as can be seen from the statistics of the coherent point selected by the method of the embodiment and the conventional method, the high coherent point selected by the two methods is concentrated between 0.4 and 0.85, but the pixel point with the coherence coefficient greater than 0.6 selected by the method of the embodiment has a significantly better occupancy rate than the pixel point selected by the conventional method, and the statistics of the two methods is shown in fig. 5.
The invention uses the coherence coefficient as an input sample of the GMM clustering algorithm to perform clustering analysis on pixels bearing coherence coefficient information; acquiring Gaussian parameters of a low coherence area in a coherence system image, and eliminating pixel points of the class by combining the parameters; overlapping the PS points on the polymerized coherence coefficient map, estimating a coherence coefficient threshold value by adopting the coherence coefficient mean value of the PS points and the coherence coefficient mean value of the two residual GMM components, and further selecting the pixel points in the two residual GMM components again; and (3) interpolating the finally obtained high-coherence region to serve as a measuring and calculating target in the SBAS-InSAR, solving linear deformation and elevation errors by utilizing SVD, and finally obtaining a deformation time sequence of the high-coherence region through estimation and elimination of each component in the interference phase. Compared with the method for directly using the coherence coefficient threshold value and the amplitude dispersion index to obtain the high coherence points, the method uses the GMM algorithm to perform cluster analysis on the coherence coefficients of ground objects in the research area to obtain the distribution situation of three coherence coefficients of high, medium and low, then well aggregates the pixel points belonging to the same kind of property together, eliminates the pixel points of the low coherence coefficient according to three kinds of Gaussian parameters, reserves the integrity of the high coherence area and the high coherence points of the high coherence area, which are partially lost due to interference factors, to the maximum extent on the basis of considering the pixels of the coherence coefficients of the same kind of property, and can reduce the influence of the missed selection of the high coherence points due to overlarge or overlarge calculation window to a certain extent.
Although the embodiments of the present invention have been described with reference to the accompanying drawings, the patentees may make various modifications or alterations within the scope of the appended claims, and are intended to be within the scope of the invention as described in the claims.

Claims (3)

1. A time sequence InSAR high coherence point selection method considering homogeneity is characterized by comprising the following steps:
s1, performing interference processing on a sequential SAR image to obtain a sequential coherence coefficient image;
s2, taking the time sequence coherence coefficient image as a sample input of a GMM-EM algorithm, acquiring clustering areas of low, medium and high pixel points of the time sequence coherence coefficient image and Gaussian parameters corresponding to each type, and removing the pixel points of the low coherence area according to the Gaussian parameters;
s3, overlapping the PS points and the clustered coherent coefficient images, calculating the coherent coefficient mean value of all the PS points and the mean value of the residual medium and high coherent coefficient pixel points in the coherent coefficient images, and acquiring a coherent coefficient threshold according to the mean value and the mean value;
s4, performing secondary screening on residual pixel points with medium and high coherence coefficients through a coherence coefficient threshold value, interpolating a coherence coefficient average value in a filtering window size for a pixel missing in a result image, and taking a final result as a high coherence point of SBAS-InSAR;
the specific steps of acquiring the clustering areas of the low, medium and high three types of pixel points of the time sequence coherence coefficient image in the step S2 are as follows:
coherence coefficient image as difference image gamma of same geographic position in different time phase SAR images d The expression is:
Figure FDA0004174172670000011
from equation (1) above, the time series of coherence coefficients of any one resolution element of the N interference pairs can be found: gamma ray 1 、γ 2 、γ 3 、...、γ N On the basis, the average value of the coherence coefficient of each pixel is calculated:
Figure FDA0004174172670000012
will be
Figure FDA0004174172670000013
As a difference image, wherein the pixels are associated with a probability density distribution +.>
Figure FDA0004174172670000014
Can be distributed by high-coherence coefficient pixel
Figure FDA0004174172670000015
Middle coherence coefficient pixel distribution->
Figure FDA0004174172670000016
And low coherence coefficient pixel-like distribution->
Figure FDA0004174172670000017
The composition is as follows:
Figure FDA0004174172670000018
wherein:
Figure FDA0004174172670000019
a coherence coefficient value representing a pixel in the difference image; omega 0 、ω 1 And omega 2 The prior probabilities of pixels of the high, medium and low coherence class are represented respectively, < >>
Figure FDA00041741726700000110
Linear combination of N components in the coherence coefficient image, namely:
Figure FDA00041741726700000111
wherein:
Figure FDA0004174172670000021
is the probability density distribution of the nth component, p (n) is the prior probability of a certain mixed component in the GMM, and the probability density function of each component obeys the gaussian distribution, namely:
Figure FDA0004174172670000022
wherein: mu (mu) n And
Figure FDA0004174172670000023
mean and variance of the nth component are represented, respectively; the gaussian distribution function satisfies:
Figure FDA0004174172670000024
/>
the prior probability of the mixed components satisfies
Figure FDA0004174172670000025
And (n) is more than or equal to 0 and less than or equal to 1;
as can be seen from formulas (4) and (5), the GMM has a total of N components, including each component
Figure FDA0004174172670000026
And the prior probability p (n) is the required solution; n=1, 2,3., N, i.e. at GMM may be determined by the following parameters:
Figure FDA0004174172670000027
the step S2 is specifically implemented by acquiring Gaussian parameters corresponding to clustering areas of low, medium and high pixel points of the time sequence coherence coefficient image, wherein the Gaussian parameters comprise the following steps:
parameter initialization: clustering samples by using a K-means clustering algorithm, and taking the mean value of each type as mu 0 And calculate
Figure FDA0004174172670000028
Taking the proportion of various samples to the total number of samples as initial posterior probability;
an estimation step: data
Figure FDA0004174172670000029
As sample data, complete data->
Figure FDA00041741726700000210
From the variable z= { z 1 ，z 2 ，z 3 ，...z N Estimate of }, where z n And difference image->
Figure FDA00041741726700000211
With the same dimensions, finish data Y d Is:
Figure FDA00041741726700000212
wherein z is n (i, j) is a posterior probability,
Figure FDA00041741726700000213
t represents the number of iterations, parameter +.>
Figure FDA00041741726700000214
Is the estimated parameter of the t-th iteration;
maximizing: parameters used in the next iteration in M steps
Figure FDA00041741726700000215
From the variable z n And (i, j) estimating, and the like, and finally obtaining various parameter values of the Gaussian mixture model, namely:
Figure FDA0004174172670000031
Figure FDA0004174172670000032
Figure FDA0004174172670000033
2. the method for selecting high coherence points of time-series InSAR according to claim 1, wherein the time-series SAR image in step S1 is a C-band high resolution time-series SAR image obtained by satellite; the interference treatment comprises at least a de-flattening, a filtering and an unwrapping.
3. The method for selecting high coherence points of time-series InSAR with homogeneity being considered as set forth in claim 1, wherein the average value of the remaining medium and high coherence coefficient pixel points in the coherence coefficient image in the step S4 is
Figure FDA0004174172670000034
All PS points coherence coefficient mean value is +.>
Figure FDA0004174172670000035
And calculating the ratio of each pel point to the total pel point number greater than the average value by the following formula:
Figure FDA0004174172670000036
/>
Figure FDA0004174172670000037
the calculation formula of the coherence coefficient threshold value D is as follows:
Figure FDA0004174172670000038
/>
CN202110614249.1A 2021-06-02 2021-06-02 Timing sequence InSAR high coherence point selection method considering homogeneity Active CN113281749B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110614249.1A CN113281749B (en) 2021-06-02 2021-06-02 Timing sequence InSAR high coherence point selection method considering homogeneity

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110614249.1A CN113281749B (en) 2021-06-02 2021-06-02 Timing sequence InSAR high coherence point selection method considering homogeneity

Publications (2)

Publication Number Publication Date
CN113281749A CN113281749A (en) 2021-08-20
CN113281749B true CN113281749B (en) 2023-05-23

Family

ID=77283259

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110614249.1A Active CN113281749B (en) 2021-06-02 2021-06-02 Timing sequence InSAR high coherence point selection method considering homogeneity

Country Status (1)

Country Link
CN (1) CN113281749B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113643284B (en) * 2021-09-09 2023-08-15 西南交通大学 Polarized synthetic aperture radar image ship detection method based on convolutional neural network
CN114187533B (en) * 2022-02-15 2022-05-03 西南交通大学 GB-InSAR (GB-InSAR) atmospheric correction method based on random forest time sequence classification
CN114488151B (en) * 2022-04-08 2022-06-24 中国科学院空天信息创新研究院 Active and passive combined detection method, device, equipment and medium for observation ship

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102608584B (en) * 2012-03-19 2014-04-16 中国测绘科学研究院 Time sequence InSAR (Interferometric Synthetic Aperture Radar) deformation monitoring method and device based on polynomial inversion model
CN105699494B (en) * 2015-12-28 2018-08-17 深圳市太赫兹科技创新研究院 Millimeter wave hologram three-dimensional image-forming detecting system and method
CN105938206B (en) * 2016-07-06 2017-09-15 华讯方舟科技有限公司 Millimeter wave safety check instrument debugging system and millimeter wave safety check instrument adjustment method
CN108802727B (en) * 2018-04-13 2021-01-15 长沙理工大学 Time sequence InSAR (interferometric synthetic Aperture Radar) highway deformation monitoring model considering rheological parameters and calculating method
CN108627834A (en) * 2018-06-07 2018-10-09 北京城建勘测设计研究院有限责任公司 A kind of subway road structure monitoring method and device based on ground InSAR
CN108957456B (en) * 2018-08-13 2020-09-04 伟志股份公司 Landslide monitoring and early identification method based on multi-data-source SBAS technology
CN110058237B (en) * 2019-05-22 2020-10-09 中南大学 InSAR point cloud fusion and three-dimensional deformation monitoring method for high-resolution SAR image
CN111812645A (en) * 2020-06-10 2020-10-23 西南交通大学 Satellite interferometry method for deformation of frozen soil in season

Also Published As

Publication number Publication date
CN113281749A (en) 2021-08-20

Similar Documents

Publication Publication Date Title
CN113281749B (en) Timing sequence InSAR high coherence point selection method considering homogeneity
US5053778A (en) Generation of topographic terrain models utilizing synthetic aperture radar and surface level data
CN109388887B (en) Quantitative analysis method and system for ground settlement influence factors
CN113960595A (en) Surface deformation monitoring method and system
Liao et al. Urban change detection based on coherence and intensity characteristics of SAR imagery
CN107064933B (en) SAR chromatography building height method based on cyclic spectrum estimation
Zhao et al. A new approach for forest height inversion using X-band single-pass InSAR coherence data
Vajedian et al. Extracting sinkhole features from time-series of TerraSAR-X/TanDEM-X data
Wu et al. A deep learning based method for local subsidence detection and InSAR phase unwrapping: Application to mining deformation monitoring
Berezowski et al. Flooding extent mapping for synthetic aperture radar time series using river gauge observations
CN113687353B (en) DS target phase optimization method based on homogeneous pixel time sequence phase matrix decomposition
CN113238228B (en) Three-dimensional earth surface deformation obtaining method, system and device based on level constraint
González et al. Relative height accuracy estimation method for InSAR-based DEMs
Feng et al. Automatic selection of permanent scatterers-based GCPs for refinement and reflattening in InSAR DEM generation
Soja et al. Model-based estimation of tropical forest biomass from notch-filtered p-band sar backscatter
Refice et al. On the use of anisotropic covariance models in estimating atmospheric DInSAR contributions
Pierce et al. Evaluation of the horizontal resolution of SRTM elevation data
El Hage et al. Effect of image‐matching parameters and local morphology on the geomorphological quality of SPOT DEMs
CN114114258A (en) Power transmission tower deformation monitoring method based on SBAS-InSAR technology
Evers et al. Concept to analyze the displacement time series of individual persistent scatterers
de Haas et al. Optimizing UAV-SfM based topographic change detection with survey co-alignment
de Moraes Frasson Using the Surface Water and Ocean Topography mission data to estimate river bathymetry and channel roughness
Huang Statistical theory for the detection of persistent scatterers in InSAR imagery
Abileah et al. Bathymetry from fusion of multi-temporal Landsat and radar altimetery
Zou et al. New building detection using SAR images with different resolutions

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