CN113281749A - Time sequence InSAR high-coherence point selection method considering homogeneity - Google Patents

Time sequence InSAR high-coherence point selection method considering homogeneity Download PDF

Info

Publication number
CN113281749A
CN113281749A CN202110614249.1A CN202110614249A CN113281749A CN 113281749 A CN113281749 A CN 113281749A CN 202110614249 A CN202110614249 A CN 202110614249A CN 113281749 A CN113281749 A CN 113281749A
Authority
CN
China
Prior art keywords
coherence
coherence coefficient
coefficient
image
time sequence
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.)
Granted
Application number
CN202110614249.1A
Other languages
Chinese (zh)
Other versions
CN113281749B (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

Abstract

The invention discloses a time sequence InSAR high-coherence point selection method considering homogeneity, which comprises the following steps: s1, performing interference processing on the time sequence SAR image to obtain a time sequence coherence coefficient image; s2, acquiring clustering areas of low, middle 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, superposing the PS points and the clustered coherence coefficient images, calculating the mean value of the coherence coefficient of all the PS points and the mean value of the residual medium and high coherence coefficient pixel points in the coherence coefficient images, and obtaining a coherence coefficient threshold; s4, secondary screening is carried out on the residual medium and high coherence coefficient pixel points through a coherence coefficient threshold, interpolation is carried out on the missing pixels in the result image by using the average value of the coherence coefficient in the size of the filtering window, and the final result is used as the high coherence point of the SBAS-InSAR. The invention furthest reserves the integrity of a high-coherence region and high coherence points of the high-coherence region which are partially lost due to interference factors on the basis of considering the pixels with the same property coherence coefficient.

Description

Time 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
China is a country with rapid development of urbanization process, the area of cities and towns is continuously increased, and ground settlement is taken as a geological disaster seriously threatening the safety of cities, so that great damage is caused to the production and the living of people and the sustainable development of cities and towns. Therefore, the method carries out long-term effective deformation monitoring on the urban area where the settlement has or is occurring, finds the settlement area and the deformation development trend thereof in time, and is a necessary condition for realizing the early maintenance of the settlement area and further changing 'passive disaster avoidance and relief' into active 'disaster prevention and control'. The traditional earth surface deformation monitoring method mostly focuses on ground monitoring technology and underground monitoring technology, such as geodetic surveying method, GPS method, drilling point location monitoring and geophysical exploration technology. The traditional deformation monitoring method is difficult to realize the arrangement of monitoring points in a large range, so that the monitoring space is limited in a complex way. In recent years, with the rapid development of communication technology and sensor technology, Interferometric Synthetic Aperture Radar (InSAR) has taken great advantage in urban ground subsidence monitoring by the characteristics that it obtains large-scale ground deformation information in all weather, all day time, cloud penetration and fog penetration, periodicity and high resolution, and particularly, the occurrence of Differential-Interferometric SAR (D-InSAR) greatly improves the ground surface micro deformation detection capability, and D-InSAR can realize centimeter-to-millimeter-level ground deformation measurement without ground control points. However, interference loss correlation, inaccurate atmospheric delay phase estimation, DEM (digital elevation model) error and other interference factors reduce the measurement precision of the D-InSAR, and meanwhile, the D-InSAR cannot acquire time sequence deformation information of a monitoring area, so that certain limitation is generated in the aspect of disaster prediction and evaluation of a ground subsidence area. Therefore, on the basis of D-InSAR, a Permanent Scatterer (PS), a small baseline set (SBAS), stamps (stanford method permanent scatterer technology), and the like are proposed in succession, collectively referred to as time series InSAR.
The SBAS-InSAR technology is proposed by the Berardino team of the Italian scholars, and the influence related to space-time loss is effectively reduced on the basis of controlling a space-time base line. The technology is combined with the atmospheric interference phase in the stable scatterer separation interference phase, and then the time sequence deformation result of the research area is obtained. The core of the technology for realizing the earth surface deformation monitoring is that statistical analysis and modeling are carried out on interference phases of scatterers with stable scattering characteristics in a time sequence image, effective estimation and separation of each component in the interference phases are realized, on the basis, fitting modeling is carried out on each component, a corresponding resolving model is obtained, and finally, the absolute deformation phase in the interference phases in a research area is obtained under the condition of minimizing other phase interference. A representative method comprises: the high coherence point is obtained by adopting a coherence coefficient threshold method, however, the method generally only considers the selection of the high coherence point according to the coherence coefficient, and the number and the reliability of the coherence points obtained according to the coherence coefficient threshold reduce the inversion accuracy of each component in the interference phase to a certain extent under the condition that the interference factors such as terrain residual phase, atmospheric phase, elevation error phase and the like are not removed. In addition, in the calculation method of the coherence coefficient, the size of a calculation window has a large influence on the selection of the coherent points, and the problem of missing selection of the isolated and effective high-coherence points is easily caused by the overlarge window; conversely, a window that is too small is likely to cause false selection of an unstable target near a high coherence point. The amplitude dispersion index threshold value method is used for obtaining the high-coherence point, the correlation between amplitude information and interference is basically irrelevant, and the interference of the SAR images with the loss correlation 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 source; the method is time-consuming in calculation, the selected high coherence point usually cannot take into account the edge and internal characteristics of the high coherence region, and the loss of the coherence point is easily caused; in areas with fewer artificial buildings, the reliability of the high coherence point obtained by using amplitude information is not high. By adopting a double-threshold detection technology, the method has a large requirement on the density of the stable scatterer, improves the reliability of high coherence points, reduces the number of coherence points and has a 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 purpose, the technical scheme adopted by the invention is as follows:
a time sequence InSAR high coherence point selecting method considering homogeneity comprises the following steps:
s1, performing interference processing on the time sequence SAR image to obtain a time sequence coherence coefficient image;
s2, taking the time sequence coherence coefficient image as sample input of a GMM-EM algorithm, obtaining clustering areas of low, middle and high pixel points of the time sequence coherence coefficient image and Gaussian parameters corresponding to each pixel point, and removing the pixel points of the low coherence area according to the Gaussian parameters;
s3, superposing the PS points and the clustered coherence coefficient images, calculating the mean value of the coherence coefficient of all the PS points and the mean value of the residual medium and high coherence coefficient pixel points in the coherence coefficient images, and obtaining a coherence coefficient threshold value according to the mean value of the coherence coefficient pixel points and the residual medium and high coherence coefficient pixel points;
s4, secondary screening is carried out on the residual medium and high coherence coefficient pixel points through a coherence coefficient threshold, interpolation is carried out on the missing pixels in the result image by using the average value of the coherence coefficient in the size of the filtering window, and the final result is used as the high coherence point of the SBAS-InSAR.
Further, the time sequence SAR image in step S1 is a C-band high resolution time sequence SAR image acquired through a satellite; the interference processing includes at least deplature, filtering, and unwrapping.
Further, the specific steps of acquiring the clustering regions of the low, medium and high types of pixel points of the time sequence coherence coefficient image in step S2 are as follows:
the coherence coefficient image is used as a difference image gamma of the same geographic position in SAR images with different time phasesdThe expression is as follows:
Figure BDA0003097385400000031
the time series of coherence coefficients for any resolution cell in the N interferometric pairs can be found according to equation (1) above: gamma ray1、γ2、γ3、...、γNOn the basis, the average value of the coherence coefficient of each pixel is calculated as follows:
Figure BDA0003097385400000032
will be provided with
Figure BDA0003097385400000033
As a difference image in which pixels are associated with probability density distributions
Figure BDA0003097385400000034
Can be distributed by high-coherence coefficient pixel-like elements
Figure BDA0003097385400000035
Medium-coherence coefficient pixel-like 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; omega0、ω1And ω2Respectively representing the prior probabilities of high, medium and low coherence classes of picture elements,
Figure BDA00030973854000000310
linear combination of N components in the image of the coherence coefficient, namely:
Figure BDA00030973854000000311
wherein:
Figure BDA00030973854000000312
is the probability density distribution of the nth component, p (n) is the prior probability of a certain mixture component in the GMM, the probability density function of each component obeys a gaussian distribution, i.e.:
Figure BDA00030973854000000313
wherein: mu.snAnd
Figure BDA00030973854000000314
respectively representing the mean and the variance of the nth component; the Gaussian distribution function satisfies:
Figure BDA00030973854000000315
a priori probability of the mixture component being satisfied
Figure BDA00030973854000000316
And p is more than or equal to 0 and less than or equal to 1 (n);
as can be seen from the formulas (4) and (5), GMM has N components in total, including each component
Figure BDA00030973854000000317
And the prior probability p (N) (1, 2, 3, N) is the solution needed, i.e. at GMM can be determined by the following parameters:
Figure BDA0003097385400000041
further, the step S2 of obtaining gaussian parameters corresponding to the clustering areas of the low, medium and high pixel points of the time sequence coherence coefficient image specifically comprises the following steps:
parameter initialization step: clustering samples by using a K-means clustering algorithm, and taking the mean value of each type as mu0And calculate
Figure BDA0003097385400000042
Taking the proportion of various samples in the total number of the samples as an initial posterior probability;
an estimation step: data of
Figure BDA0003097385400000043
As sample data, complete data
Figure BDA0003097385400000044
By the variable Z ═ { Z1，z2，z3，...，zNIs determined by an estimate of z, where z isnAnd difference image
Figure BDA0003097385400000045
Having the same dimension, finish data YdThe log-likelihood function of (a) is:
Figure BDA0003097385400000046
wherein z isn(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;
a maximization step: parameters used by next iteration in M steps
Figure BDA0003097385400000049
By a variable zn(i, j) estimating, and so on, 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 middle and high coherence coefficient pixel points in the coherence coefficient image in step S4 is
Figure BDA00030973854000000413
The mean value of all the PS point coherence coefficients is
Figure BDA00030973854000000414
And calculating the ratio of each pixel point to the total pixel point number by the following formula, wherein the ratio is greater than the average value:
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: the method uses the GMM algorithm to perform cluster analysis on the coherence coefficients of the ground objects in the researched area, obtains the distribution conditions of high, medium and low coherence coefficients, better aggregates pixel points belonging to the same property together, eliminates the pixel points with low coherence coefficients according to three types of Gaussian parameters, keeps the integrity of the high coherence area and high coherence points partially lost due to interference factors on the basis of considering the pixel with the same property to the maximum extent, and can reduce the influence of missed selection and wrong selection of the high coherence points caused by overlarge or undersize of a calculation window to a certain extent.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to the drawings without creative efforts.
Fig. 1 is a flow chart of the time-series InSAR high coherence point selection method considering homogeneity according to the present invention.
FIG. 2 is a diagram of the combination of interference pairs in the present invention.
Fig. 3 is a gaussian mixture distribution histogram of the coherence coefficient in the present invention.
FIG. 4 is a graph of the coherence factor and PS dot overlap distribution in the present invention.
FIG. 5 is a graph of statistics of coherence points for a conventional method of the present invention and a 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, and the scope of the present invention will be more clearly and clearly defined.
Referring to fig. 1, in order to accurately obtain a high coherence point in a coherent image and overcome the defect that the conventional coherent coefficient threshold based on manual determination is used to obtain the high coherence point and cannot retain edge information of a high coherence region, the embodiment mainly considers the problems of selection precision and number of high coherence scatterers with the same property, and uses a GMM algorithm to cluster homogeneous points first and then sets a threshold in combination with the coherence coefficient of a PS point, so as to maximally retain the spatial distribution of the high coherence region. The embodiment discloses a time sequence InSAR high coherence point selection method considering homogeneity, which specifically comprises the following steps:
and step S1, performing interference processing on the time sequence SAR image to acquire a time sequence coherence coefficient image.
In the embodiment, a high-resolution time sequence SAR image of a C wave band is obtained through a European Sentinel-1A (Sentinel-1A) satellite, optimal interference is obtained according to a preset space-time baseline threshold value to perform interference processing on a combination, and a time sequence coherence coefficient image is obtained through interference processing such as flattening, filtering and unwrapping. The interference pair combination result is shown in fig. 2.
And step S2, taking the time sequence coherence coefficient image as sample input of the GMM-EM algorithm, obtaining clustering areas of low, middle and high pixel points of the time sequence coherence coefficient image and Gaussian parameters corresponding to each pixel point, and removing the pixel points of the low coherence area according to the Gaussian parameters.
GMM is independent based on each pixel value in each interference image, and the coherence coefficient image is used as a difference image gamma of the same geographic position in SAR images in different time phasesdThe expression is:
Figure BDA0003097385400000061
the time series of coherence coefficients for any resolution cell in the N interferometric pairs can be found according to equation (1): gamma ray1、γ2、γ3、...、γN. On the basis, the average value of the coherence coefficient of each pixel is calculated as follows:
Figure BDA0003097385400000062
will be provided with
Figure BDA0003097385400000063
As a difference image in which pixels are associated with probability density distributions
Figure BDA0003097385400000064
Can be distributed by high-coherence coefficient pixel-like elements
Figure BDA0003097385400000065
Medium-coherence coefficient pixel-like 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; omega0、ω1And ω2Respectively representing the prior probabilities of high, medium and low coherence class pixels.
Figure BDA00030973854000000610
Linear combination of N components in the image of the coherence coefficient, namely:
Figure BDA00030973854000000611
wherein:
Figure BDA00030973854000000612
is the probability density distribution of the nth component, p (n) is the prior probability of a certain mixture component in the GMM, the probability density function of each component obeys a gaussian distribution, i.e.:
Figure BDA0003097385400000071
wherein: mu.snAnd
Figure BDA0003097385400000072
respectively representing the mean and the variance of the nth component; the Gaussian distribution function satisfies:
Figure BDA0003097385400000073
a priori probability of the mixture component being satisfied
Figure BDA0003097385400000074
And p is more than or equal to 0 and less than or equal to 1 (n).
As can be seen from the formulas (4) and (5), GMM has N components in total, including each component
Figure BDA0003097385400000075
And the prior probability p (N) (1, 2, 3, N) is the solution needed, i.e. at GMM can be determined by the following parameters:
Figure BDA0003097385400000076
the present embodiment adopts an EM algorithm to solve parameters in a gaussian mixture model, which is an iterative approximation method. Typically a set of parameters is randomly initialized
Figure BDA0003097385400000077
And calculate
Figure BDA0003097385400000078
The posterior probability p (n) generated by a certain component in the GMM is used for iterative optimal value solution, and each iteration consists of two processes of E step (solving the maximum expectation of an observed value) and M step (solving the maximization). Two steps
Figure BDA0003097385400000079
Alternating until a local optimum is converged, wherein the steps are independent of each other, and the specific steps are as follows:
a parameter initialization step, in which this embodiment clusters samples using a K-Means (K-Means) clustering algorithm, and uses the mean of each type as μ0And calculate
Figure BDA00030973854000000710
And taking the proportion of each type of sample in the total number of the samples as the initial posterior probability.
Estimation step, data
Figure BDA00030973854000000711
As sample data, complete data
Figure BDA00030973854000000712
By the variable Z ═ { Z1，z2，z3，...，zNIs determined by an estimate of z, where z isnAnd difference image
Figure BDA00030973854000000713
Having the same dimension, finish data YdThe log-likelihood function of (a) is:
Figure BDA00030973854000000714
wherein z isn(i, j) is a posterior probability, i.e.:
Figure BDA0003097385400000081
wherein the content of the first and second substances,t represents the number of iterations, parameter
Figure BDA0003097385400000082
Is the estimated parameter for the t-th iteration.
A maximization step: parameters used by next iteration in M steps
Figure BDA0003097385400000083
Can be represented by the variable zn(i, j) estimating, and so on, and finally obtaining various parameter values of the Gaussian mixture model, namely:
Figure BDA0003097385400000084
Figure BDA0003097385400000085
Figure BDA0003097385400000086
presetting a component n in GMM to be 3 in an EM algorithm; according to the number of sample points obtained in 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 the three coherence coefficient intervals are mixed is shown in fig. 3 (a). Determining the mean value of each Gaussian distribution in the distribution range of each component coherence coefficient 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 mixture components of the GMM based on the mean value of the coherence coefficient 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 area is, further, the more reliable the deformation phase obtained by inverting the stable area is, however, in the research area, vegetation, forest land, water body and the like are all low coherence areas, and the coherence coefficient is generally distributed between 0 and 0.3; the soil quality of the bare leaked soil or the doped rock is a medium coherent region, and the coherence coefficient is distributed between 0.3 and 0.55; the artificial buildings such as houses, roads and the like are high-coherence areas, and the coherence coefficient is distributed between 0.55 and 1. This is the basis for dividing the original image with coherent coefficients into 3 types of pixels, low, medium and high. In the embodiment, a GMM-EM algorithm is adopted to perform cluster analysis on pixel points in a research area for removing low-coherence ground objects such as vegetation and water bodies, the pixel points with the same property (namely, the pixel points with coherence coefficients 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 with the coherence coefficients is respectively low, medium and high; the method is the basis for dividing the clustered pixels of the coherence coefficient into three categories of low, medium and high.
And S3, superposing the PS points and the clustered coherence coefficient images, calculating the mean value of the coherence coefficient of all the PS points and the mean value of the residual medium and high coherence coefficient pixel points in the coherence coefficient images, and obtaining a coherence coefficient threshold value 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 the inversion of the deformation information of the target area completely depending on a high coherence point in the coherence coefficient image is unreliable to a certain extent. The PS points can show similar and larger intensity information in all time sequence images, the phase dispersion is smaller, the points do not need to analyze the phase information, and the points also show the characteristic of high coherence coefficient after interference processing. Therefore, the PS point is selected under the condition that whether the scattering characteristic of the pixel is stable or not is considered, and the PS point is more reliable than the coherent point directly selected through the coherence coefficient threshold value.
In the present embodiment, the coherence coefficient threshold is calculated in consideration of the PS point and the homogeneous coherence point acquired by the GMM. Firstly, dividing pixels in a coherence coefficient graph into three types according to the size of a coherence coefficient: and low, medium and high, acquiring the distribution condition of pixel points with the same property in a research area through a GMM-EM algorithm, and removing the low-coherence coefficient component according to the Gaussian parameter of the homogeneous point. Respectively calculating the average value of the coherence coefficients of the residual image element points
Figure BDA0003097385400000091
And the mean value of the correlation coefficient of the PS point
Figure BDA0003097385400000092
And the proportion of each pixel point larger than the average value in the total pixel point number.
Figure BDA0003097385400000093
Figure BDA0003097385400000094
Finally, selecting a threshold D capable of retaining the number of the high-coherence points and the accuracy thereof to the maximum extent according to the parameters, wherein the calculation formula is as follows:
Figure BDA0003097385400000095
and step S4, secondary screening is carried out on the residual medium and high coherence coefficient pixel points through a coherence coefficient threshold value D, interpolation is carried out on the missing pixels in the resulting image by using the average value of the coherence coefficient in the size of a filtering window, so as to recover the missing high coherence pixels caused by other interference factors, and the final result is used as the high coherence point of the SBAS-InSAR. And solving linear deformation and elevation errors by using SVD (singular value decomposition), and finally obtaining the time sequence deformation phase of the high coherence point by estimating and eliminating each component in the interference phase of the high coherence region.
In this embodiment, PS points are superimposed on a coherence coefficient map with low coherence regions removed, the superimposed situation of sub-regions within a study range is shown in fig. 4, and the distribution situation of artificial buildings within a region can be seen from the amplitude mean value map in fig. 4(a), in view of that there are many high coherence points (artificial buildings) in the study range and the distribution is dense, the threshold value of the coherence coefficient when PS points are selected is set to 0.55, and then the coherence coefficient of PS points is distributed between 0.55 and 0.9. Fig. 4(b) shows that the PS dot positions are consistent with the distribution of the high-brightness pixel points in the amplitude mean graph. 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 area selected according to the position keeps more reliable. The coherence coefficient mean distribution graph is shown in fig. 4(c), and the superposition condition of PS points and coherence coefficients is shown in fig. 4(d), which shows that the pixel points with high coherence coefficients have the same consistency with the distribution of PS points and high-brightness pixel points in the amplitude mean graph. And then setting a coherence coefficient threshold value to select the final high-quality coherence point according to the distribution range of the coherence coefficient of the coincidence region in the PS point and coherence coefficient mean value graph. The mean value of the coherence coefficients of the low coherence distribution in the three components of the coherence coefficient Gaussian mixture model is 0.4, the lowest coherence coefficient of the PS points is 0.55, the calculation result according to the formula is taken as a threshold value, namely the coherence coefficient threshold value is set to be 0.4, and finally 518576 high coherence points are determined. Finally, it can be seen from the statistical conditions of the method of the present embodiment and the conventional method for selecting the coherence points that the high coherence points selected by the two methods are all concentrated between 0.4 and 0.85, but the pixel point ratio of the coherence coefficient greater than 0.6 selected by the method of the present embodiment is significantly better than the pixel point selected by the conventional method, and the statistical results of the two methods are shown in fig. 5.
The method uses the coherence coefficient as an input sample of a GMM clustering algorithm to perform clustering analysis on pixels bearing coherence coefficient information; obtaining Gaussian parameters of low coherence areas in coherence system images, and removing the class of image element points by combining the parameters; the PS points are superposed on the polymerized coherence coefficient graph, the coherence coefficient threshold is estimated by adopting the coherence coefficient mean of the PS points and the coherence coefficient mean of the two remaining GMM components, and then the pixel points in the two remaining GMM components are selected again; and finally, interpolating the finally obtained high-coherence area to be used as a measuring and calculating target in the SBAS-InSAR, solving linear deformation and elevation errors by using SVD (singular value decomposition), and finally obtaining a deformation time sequence of the high-coherence area through estimation and elimination of each component in an interference phase. Compared with high coherence points obtained by directly using a coherence coefficient threshold and an amplitude dispersion index, the method uses a GMM algorithm to perform cluster analysis on the coherence coefficients of the ground objects in a research area to obtain the distribution conditions of high, medium and low three types of coherence coefficients, then better aggregates the image element points belonging to the same property together, removes the image element points with low coherence coefficients according to three types of Gaussian parameters, furthest retains the integrity of the high coherence area and high coherence points partially lost due to interference factors on the basis of considering the image elements with the same property, and can reduce the influence of missing and wrong selection of the high coherence points due to overlarge or undersize of a calculation window to a certain extent.
Although the embodiments of the present invention have been described with reference to the accompanying drawings, various changes or modifications may be made by the patentees within the scope of the appended claims, and within the scope of the invention, as long as they do not exceed the scope of the invention described in the claims.

Claims (5)

1. A time sequence InSAR high coherence point selection method considering homogeneity is characterized by comprising the following steps:
s1, performing interference processing on the time sequence SAR image to obtain a time sequence coherence coefficient image;
s2, taking the time sequence coherence coefficient image as sample input of a GMM-EM algorithm, obtaining clustering areas of low, middle and high pixel points of the time sequence coherence coefficient image and Gaussian parameters corresponding to each pixel point, and removing the pixel points of the low coherence area according to the Gaussian parameters;
s3, superposing the PS points and the clustered coherence coefficient images, calculating the mean value of the coherence coefficient of all the PS points and the mean value of the residual medium and high coherence coefficient pixel points in the coherence coefficient images, and obtaining a coherence coefficient threshold value according to the mean value of the coherence coefficient pixel points and the residual medium and high coherence coefficient pixel points;
s4, secondary screening is carried out on the residual medium and high coherence coefficient pixel points through a coherence coefficient threshold, interpolation is carried out on the missing pixels in the result image by using the average value of the coherence coefficient in the size of the filtering window, and the final result is used as the high coherence point of the SBAS-InSAR.
2. The method for selecting a time sequence InSAR high-coherence point considering homogeneity according to claim 1, wherein 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 processing includes at least deplature, filtering, and unwrapping.
3. The method for selecting time series InSAR high-coherence point considering homogeneity according to claim 1, wherein the specific steps of obtaining the clustering areas of the low, middle and high types of pixel points of the time series coherence coefficient image in the step S2 are as follows:
the coherence coefficient image is used as a difference image gamma of the same geographic position in SAR images with different time phasesdThe expression is as follows:
Figure FDA0003097385390000011
the time series of coherence coefficients for any resolution cell in the N interferometric pairs can be found according to equation (1) above: gamma ray1、γ2、γ3、...、γNOn the basis, the average value of the coherence coefficient of each pixel is calculated as follows:
Figure FDA0003097385390000012
will be provided with
Figure FDA0003097385390000019
As a difference image in which pixels are associated with probability density distributions
Figure FDA0003097385390000013
Can be distributed by high-coherence coefficient pixel-like elements
Figure FDA0003097385390000014
Medium-coherence coefficient pixel-like distribution
Figure FDA0003097385390000015
And low coherence coefficient pixel-like distribution
Figure FDA0003097385390000016
The composition is as follows:
Figure FDA0003097385390000017
wherein:
Figure FDA0003097385390000018
a coherence coefficient value representing a pixel in the difference image; omega0、ω1And ω2Respectively representing the prior probabilities of high, medium and low coherence classes of picture elements,
Figure FDA0003097385390000021
linear combination of N components in the image of the coherence coefficient, namely:
Figure FDA0003097385390000022
wherein:
Figure FDA0003097385390000023
is the probability density distribution of the nth component, p (n) is the prior probability of a certain mixture component in the GMM, the probability density function of each component obeys a gaussian distribution, i.e.:
Figure FDA0003097385390000024
wherein: mu.snAnd
Figure FDA0003097385390000025
respectively representing the mean and the variance of the nth component; the Gaussian distribution function satisfies:
Figure FDA0003097385390000026
a priori probability of the mixture component being satisfied
Figure FDA0003097385390000027
And p is more than or equal to 0 and less than or equal to 1 (n);
as can be seen from the formulas (4) and (5), GMM has N components in total, including each component
Figure FDA0003097385390000028
And the prior probability p (N) (1, 2, 3, N) is the solution needed, i.e. at GMM can be determined by the following parameters:
Figure FDA0003097385390000029
4. the method for selecting the time sequence InSAR high coherence point considering the homogeneity according to claim 1, wherein the step S2 of obtaining Gaussian parameters corresponding to the clustering regions of the low, medium and high pixel points of the time sequence coherence coefficient image comprises the following specific steps:
parameter initialization step: clustering samples by using a K-means clustering algorithm, and taking the mean value of each type as mu0And calculate
Figure FDA00030973853900000210
Taking the proportion of various samples in the total number of the samples as an initial posterior probability;
an estimation step: data of
Figure FDA00030973853900000211
As sample data, complete data
Figure FDA00030973853900000212
By the variable Z ═ { Z1，z2，z3，...，zNIs determined by the estimation of the frequency band,wherein z isnAnd difference image
Figure FDA00030973853900000213
Having the same dimension, finish data YdThe log-likelihood function of (a) is:
Figure FDA00030973853900000214
wherein z isn(i, j) is a posterior probability,
Figure FDA0003097385390000031
t represents the number of iterations, parameter
Figure FDA0003097385390000032
Is the estimated parameter of the t-th iteration;
a maximization step: parameters used by next iteration in M steps
Figure FDA0003097385390000033
By a variable zn(i, j) estimating, and so on, and finally obtaining various parameter values of the Gaussian mixture model, namely:
Figure FDA0003097385390000034
Figure FDA0003097385390000035
Figure FDA0003097385390000036
5. the homogeneity-aware sequential InSAR high coherence point selection method according to claim 1Characterized in that the average value of the residual middle and high coherence coefficient image element points in the coherence coefficient image in the step S4 is
Figure FDA0003097385390000037
The mean value of all the PS point coherence coefficients is
Figure FDA0003097385390000038
And calculating the ratio of each pixel point to the total pixel point number by the following formula, wherein the ratio is greater than the average value:
Figure FDA0003097385390000039
Figure FDA00030973853900000310
the calculation formula of the coherence coefficient threshold value D is as follows:
Figure FDA00030973853900000311
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 true CN113281749A (en) 2021-08-20
CN113281749B 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)

Cited By (3)

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

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102608584A (en) * 2012-03-19 2012-07-25 中国测绘科学研究院 Time sequence InSAR (Interferometric Synthetic Aperture Radar) deformation monitoring method and device based on polynomial inversion model
CN108627834A (en) * 2018-06-07 2018-10-09 北京城建勘测设计研究院有限责任公司 A kind of subway road structure monitoring method and device based on ground InSAR
CN108802727A (en) * 2018-04-13 2018-11-13 长沙理工大学 A kind of sequential InSAR highway deformation monitoring models and calculation method for taking rheological parameter into account
CN108957456A (en) * 2018-08-13 2018-12-07 伟志股份公司 Landslide monitoring and EARLY RECOGNITION method based on multi-data source SBAS technology
US20190004171A1 (en) * 2015-12-28 2019-01-03 Shenzhen Institute Of Terahertz Technology And Innovation Millimeter wave holographic three-dimensional imaging detection system and method
US20190227190A1 (en) * 2016-07-06 2019-07-25 China Communication Technology Co., Ltd System and method for debugging millimeter wave security inspection instrument
CN110058237A (en) * 2019-05-22 2019-07-26 中南大学 InSAR point Yun Ronghe and three-dimensional deformation monitoring method towards High-resolution SAR Images
CN111812645A (en) * 2020-06-10 2020-10-23 西南交通大学 Satellite interferometry method for deformation of frozen soil in season

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102608584A (en) * 2012-03-19 2012-07-25 中国测绘科学研究院 Time sequence InSAR (Interferometric Synthetic Aperture Radar) deformation monitoring method and device based on polynomial inversion model
US20190004171A1 (en) * 2015-12-28 2019-01-03 Shenzhen Institute Of Terahertz Technology And Innovation Millimeter wave holographic three-dimensional imaging detection system and method
US20190227190A1 (en) * 2016-07-06 2019-07-25 China Communication Technology Co., Ltd System and method for debugging millimeter wave security inspection instrument
CN108802727A (en) * 2018-04-13 2018-11-13 长沙理工大学 A kind of sequential InSAR highway deformation monitoring models and calculation method for taking rheological parameter into account
CN108627834A (en) * 2018-06-07 2018-10-09 北京城建勘测设计研究院有限责任公司 A kind of subway road structure monitoring method and device based on ground InSAR
CN108957456A (en) * 2018-08-13 2018-12-07 伟志股份公司 Landslide monitoring and EARLY RECOGNITION method based on multi-data source SBAS technology
CN110058237A (en) * 2019-05-22 2019-07-26 中南大学 InSAR point Yun Ronghe and three-dimensional deformation monitoring method towards High-resolution SAR Images
CN111812645A (en) * 2020-06-10 2020-10-23 西南交通大学 Satellite interferometry method for deformation of frozen soil in season

Cited By (6)

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

Also Published As

Publication number Publication date
CN113281749B (en) 2023-05-23

Similar Documents

Publication Publication Date Title
CN113281749A (en) Time sequence InSAR high-coherence point selection method considering homogeneity
Zhao et al. Generation of long-term InSAR ground displacement time-series through a novel multi-sensor data merging technique: The case study of the Shanghai coastal area
CN109388887B (en) Quantitative analysis method and system for ground settlement influence factors
CN106950556A (en) Heritage area deformation monitoring method based on distributed diffusion body sequential interference SAR technology
CN113960595A (en) Surface deformation monitoring method and system
CN102144174A (en) Identification and analysis of persistent scatterers in series of sar images
CN112444188B (en) Multi-view InSAR sea wall high-precision three-dimensional deformation measurement method
Wang et al. A new likelihood function for consistent phase series estimation in distributed scatterer interferometry
Jiang et al. Delineation of built-up land change from SAR stack by analysing the coefficient of variation
CN112329265A (en) Satellite remote sensing rainfall refinement space estimation method and system
Hakim et al. InSAR time-series analysis and susceptibility mapping for land subsidence in Semarang, Indonesia using convolutional neural network and support vector regression
CN115079172A (en) MTInSAR landslide monitoring method, equipment and storage medium
Vajedian et al. Extracting sinkhole features from time-series of TerraSAR-X/TanDEM-X data
Guan et al. Fusion of public DEMs based on sparse representation and adaptive regularization variation model
CN113238228B (en) Three-dimensional earth surface deformation obtaining method, system and device based on level constraint
CN114325706A (en) Distributed scatterer filtering method
CN113687353A (en) DS target phase optimization method based on homogeneous pixel time sequence phase matrix decomposition
Yaseen et al. Local interpolation of coseismic displacements measured by InSAR
CN112051572A (en) Method for monitoring three-dimensional surface deformation by fusing multi-source SAR data
Refice et al. On the use of anisotropic covariance models in estimating atmospheric DInSAR contributions
CN116299455A (en) Facility deformation analysis method based on PSInSAR and SqueseAR
Feng et al. Automatic selection of permanent scatterers-based GCPs for refinement and reflattening in InSAR DEM generation
CN113008202B (en) Ground settlement monitoring method integrating different synthetic aperture radar interferometry
CN105551021B (en) The building method of estimation of falling loss rate based on multidate full-polarization SAR
CN113625241A (en) Differential settlement monitoring and early warning method

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