Change Detection Method for Multi-temporal SAR Imagery

Abstract

The present invention provides a multi-phase SAR image change detection method, comprising: constructing a similarity matrix between a first covariance matrix of a first-phase dual-polarization SAR image and a second covariance matrix of a second-phase dual-polarization SAR image; calculating a difference map between the first-phase dual-polarization SAR image and the second-phase dual-polarization SAR image according to the similarity matrix; performing initial segmentation on the difference map to obtain an initial label of each pixel in the difference map, wherein the initial label is used to distinguish a change class and an unchanged class in the difference map; according to the initial label, optimizing the initial labels of the change class and the unchanged class in the difference map based on a generalized Gaussian distribution and a contextual relationship between each pixel and the nearest pixel in an eight-neighborhood direction to obtain a change binary map.

Description

Technical Field

[0001] The present disclosure relates to the technical field of radar image processing, and in particular to a multi-temporal SAR image change detection method. Background Art

[0002] Synthetic Aperture Radar (SAR) image change detection uses SAR images covering the same area at different times to compare and analyze to obtain information on surface cover changes. The basic process of change detection mainly includes three technical links, namely: preprocessing → difference map generation → difference map analysis (segmentation). Difference map generation and difference map analysis are the top priorities of the change detection process. However, in the existing technology, change detection is basically carried out around single-polarization or full-polarization radar image data, and detection of dual-polarization radar image data is very rare. Summary of the invention

[0003] In view of this, the present disclosure provides, on one hand, a multi-phase SAR image change detection method, including: constructing a similarity matrix between a first covariance matrix of a first-phase dual-polarization SAR image and a second covariance matrix of a second-phase dual-polarization SAR image; calculating a difference map between the first-phase dual-polarization SAR image and the second-phase dual-polarization SAR image according to the similarity matrix; performing initial segmentation on the difference map to obtain an initial labeling of each pixel in the difference map, wherein the initial labeling is used to distinguish between a changed class and an unchanged class in the difference map; according to the initial labeling, optimizing the initial labels of the changed class and the unchanged class in the difference map based on a generalized Gaussian distribution and a contextual relationship between each pixel and the nearest pixel in an eight-neighborhood direction to obtain a binary change map.

[0004] According to an embodiment of the present disclosure, the construction of a similarity matrix between the covariance matrix of a first phase dual-polarization SAR image and the covariance matrix of a second phase dual-polarization SAR image includes: constructing the first covariance matrix according to backscattering information of a co-polarization channel and backscattering information of a cross-polarization channel corresponding to the first phase dual-polarization SAR image; constructing the second covariance matrix according to backscattering information of a co-polarization channel and backscattering information of a cross-polarization channel corresponding to the second phase dual-polarization SAR image; and constructing the similarity matrix according to the first covariance matrix and the second covariance matrix.

[0005] According to an embodiment of the present disclosure, the similarity matrix R is:

[0006]

[0007]

[0008] Among them, S xx is the backscattering information of the co-polarization channel, S xy is the backscatter information of the cross-polarization channel, * indicates complex conjugate calculation; <·> indicates window averaging calculation, |·| indicates complex modulus calculation, represents the conjugate transpose of the matrix, C1 is the first covariance matrix, and C2 is the second covariance matrix.

[0009] According to an embodiment of the present disclosure, the calculating of the difference map between the first phase dual-polarization SAR image and the second phase dual-polarization SAR image according to the similarity matrix includes: solving the maximum eigenvalue and the minimum eigenvalue of the similarity matrix; calculating the similarity measure between the first phase dual-polarization SAR image and the second phase dual-polarization SAR image according to the maximum eigenvalue and the minimum eigenvalue; and obtaining the difference map according to the similarity measure.

[0010] According to an embodiment of the present disclosure, the similarity measure SIM is:

[0011]

[0012] Among them, λ max , min are the maximum eigenvalue and minimum eigenvalue of the similarity matrix respectively.

[0013] According to an embodiment of the present disclosure, the initial segmentation of the difference map to obtain the initial labeling of each pixel in the difference map includes: statistical modeling of the changed class and the unchanged class in the difference map based on a generalized Gaussian distribution model to obtain statistical characteristics of the changed class and the unchanged class; based on the statistical characteristics, using the Kittler-Illingworth threshold selection criterion to determine the threshold to segment the difference map to obtain the initial labeling of each pixel.

[0014] According to an embodiment of the present disclosure, according to

[0015]

[0016]

[0017]

[0018] Construct the generalized Gaussian distribution model, where X l is the grayscale of the grayscale image of the difference map, X l =0, 1, 2, ..., L-1, L is the maximum pixel value of the grayscale image, m iis the mean value of the pixels in the changed or unchanged class in the difference map, β i is the shape parameter of the generalized Gaussian distribution model, σ i is the standard deviation of the changed or unchanged class in the difference map, Γ(·) is the gamma function, ω i is a change class or a change class, i represents the number of the change class or the unchanged class, c represents the change class, and u represents the unchanged class; wherein the Kittler-Illingworth threshold selection criterion is to meet the minimum cost function as the criterion, determine the threshold in the grayscale range [0, L-1], and obtain the m according to the determined threshold i , σ i and β i .

[0019] According to an embodiment of the present disclosure, wherein, according to the initial labeling, the initial labeling of the changed class and the unchanged class in the difference map is optimized based on the generalized Gaussian distribution and the contextual relationship between each pixel and the nearest pixel in the eight-neighborhood direction, including: according to the initial labeling, based on the generalized Gaussian distribution and the contextual relationship between each pixel and the nearest pixel in the eight-neighborhood direction, constructing the Markov random field model energy function of the difference map; solving the minimized Markov random field model energy function to obtain the optimal labeling of each pixel in the difference map; outputting the change binary map according to the optimal labeling.

[0020] According to an embodiment of the present disclosure, the energy function of the Markov random field model is:

[0021]

[0022] D p (f p )=ln(G(X l |f p )),

[0023]

[0024] Wherein, P represents the set of all pixels in the difference map, p represents the pixel number, q represents the number of the nearest pixel of the p-th pixel in the eight-neighborhood direction, N is the set consisting of the p-th pixel and the nearest pixel, and f p represents the initial label of the p-th pixel, f q represents the initial label of the qth pixel, φ is the smoothing parameter, and λ is the balancing factor.

[0025] According to an embodiment of the present disclosure, the energy function of the Markov random field model for solving the minimization includes: randomly selecting a pixel from the difference map, and when the category mark of the pixel is a changed class, obtaining an unchanged class based on the initial mark; calculating the difference between the energies corresponding to the initial mark of the changed class and the initial mark of the unchanged class based on the energy function of the Markov random field model; if the difference meets a preset condition, updating the category mark of the current pixel from the changed class to the unchanged class; iterating in this way, traversing other pixels in the difference map, and obtaining the sum of the energy differences corresponding to the difference map, until the sum of the energy differences is less than a preset threshold or the number of iterations reaches a preset number, and obtaining the optimal labeling of each pixel in the difference map. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] The above and other objects, features and advantages of the present disclosure will become more apparent through the following description of the embodiments of the present disclosure with reference to the accompanying drawings, in which:

[0027] Figure 1 The flowchart of the multi-temporal SAR image change detection method provided by the embodiment of the present disclosure is schematically shown. DETAILED DESCRIPTION

[0028] In order to make the purpose, technical solutions and advantages of the present disclosure more clearly understood, the present disclosure is further described in detail below in combination with specific embodiments and with reference to the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present disclosure, rather than all of the embodiments. Based on the embodiments in the present disclosure, all other embodiments obtained by ordinary technicians in the field without creative work are within the scope of protection of the present disclosure.

[0029] The terms used herein are only for describing specific embodiments and are not intended to limit the present disclosure. The terms "comprise", "include", etc. used herein indicate the existence of the features, steps, operations and / or components, but do not exclude the existence or addition of one or more other features, steps, operations or components.

[0030] In the present disclosure, unless otherwise clearly specified and limited, the terms "installed", "connected", "connected", "fixed" and the like should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection, an electrical connection, or can communicate with each other; it can be a direct connection, or it can be indirectly connected through an intermediate medium, it can be the internal connection of two elements or the interaction relationship between two elements. For ordinary technicians in this field, the specific meanings of the above terms in the present disclosure can be understood according to specific circumstances.

[0031] In the description of the present disclosure, it is necessary to understand that the terms "longitudinal", "length", "circumferential", "front", "rear", "left", "right", "top", "bottom", "inside", "outside", etc., indicating the orientation or positional relationship, are based on the orientation or positional relationship shown in the accompanying drawings, and are only for the convenience of describing the present disclosure and simplifying the description, and do not indicate or imply that the subsystem or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be understood as a limitation on the present disclosure.

[0032] Throughout the drawings, the same elements are represented by the same or similar reference numerals. Conventional structures or configurations will be omitted when they may cause confusion in the understanding of the present disclosure. The shapes, sizes, and positional relationships of the components in the drawings do not reflect the actual size, proportion, and actual positional relationship. In addition, in the claims, any reference symbol between brackets shall not be construed as a limitation to the claims.

[0033] Similarly, in order to simplify the present disclosure and help understand one or more of the various disclosed aspects, in the above description of the exemplary embodiments of the present disclosure, the various features of the present disclosure are sometimes grouped together into a single embodiment, figure, or description thereof. The description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" and the like means that the specific features, structures, materials, or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present disclosure. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described may be combined in any one or more embodiments or examples in a suitable manner.

[0034] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more of the feature. In the description of the present disclosure, the meaning of "plurality" is at least two, such as two, three, etc., unless otherwise clearly and specifically defined.

[0035] The purpose of the present invention is to propose a multi-phase dual-polarization SAR automatic change detection method based on similarity matrix eigenvalue and Markov random field model for the characteristics of dual-polarization SAR data. The method first constructs a similarity matrix between the covariance matrices of the dual-polarization SAR data of the previous and next phases, and constructs a similarity measure by calculating the maximum and minimum eigenvalues ​​of the similarity matrix to calculate the difference map between the previous and next phase images. Then, the Kittler-Illingworth automatic threshold method based on generalized Gaussian distribution and the image segmentation method based on Markov random field model are combined to realize the automatic and accurate segmentation of the change class and the unchanged class in the difference map. The generalized Gaussian distribution can be used to better fit the statistical characteristics of the change class and the unchanged class in the difference map, and the threshold is automatically determined in combination with the Kittler-Illingworth threshold selection criterion, so as to realize the initial marking of the change area and the unchanged area in the difference map. However, SAR uses radar echo signal coherent imaging, and the obtained image often has more serious "granular" speckle noise, which is easy to cause false detection. Therefore, the energy function of the Markov random field model of the difference map is constructed by comprehensively utilizing the generalized Gaussian distribution and the eight-neighborhood spatial structure. The eight-neighborhood spatial structure can well describe the advantages of the context information of each pixel in the difference map, so as to effectively suppress the false detection caused by speckle noise. At the same time, the simulated annealing global optimal optimization algorithm based on the Metropolis criterion is adopted to continuously iterate and obtain the optimal label of each pixel in the difference map, so as to achieve accurate segmentation of the changed and unchanged classes.

[0036] Figure 1 The flowchart of the multi-temporal SAR image change detection method provided by the embodiment of the present disclosure is schematically shown.

[0037] like Figure 1 As shown, the method may include, for example, operations S101 to S104.

[0038] In operation S101, a similarity matrix is ​​constructed between a first covariance matrix of a dual-polarization SAR image of a first phase and a second covariance matrix of a dual-polarization SAR image of a second phase.

[0039] In the embodiments of the present disclosure, multi-temporal phases generally refer to the characteristics of a set of remote sensing images in a time series, that is, a set of images of the same region acquired at different times. The first temporal phase may be before the second temporal phase.

[0040] In the disclosed embodiment, the first covariance matrix can be constructed based on the backscattering information of the co-polarization channel and the backscattering information of the cross-polarization channel corresponding to the first phase dual-polarization SAR image. The second covariance matrix is ​​constructed based on the backscattering information of the co-polarization channel and the backscattering information of the cross-polarization channel corresponding to the second phase dual-polarization SAR image. The backscattering information refers to the collection of the waves reflected from the target after the radar sends a detection wave to the target (such as the ground, the ocean) to obtain the relevant information of the target contained in the radar image.

[0041] For example, the covariance matrix C can be defined as:

[0042]

[0043] Among them, S xx is the backscatter information of the common polarization channel in dual-polarization SAR (generally complex data), S xy is the backscatter information of the cross-polarization channel in the dual-polarization SAR. The subscript x represents horizontal polarization or vertical polarization, the subscript y represents vertical polarization or horizontal polarization, * represents complex conjugate calculation, <·> represents window averaging calculation, and |·| represents complex modulus calculation.

[0044] Then construct a similarity matrix based on the first covariance matrix and the second covariance matrix. The similarity matrix R can be defined as:

[0045]

[0046] in, represents the conjugate transpose of the matrix, C1 is the first covariance matrix, and C2 is the second covariance matrix.

[0047] In operation S102, a difference map between the dual-polarization SAR image at the first time phase and the dual-polarization SAR image at the second time phase is calculated according to the similarity matrix.

[0048] In the disclosed embodiment, the maximum eigenvalue and the minimum eigenvalue of the similarity matrix can be solved, and then the similarity measure between the first phase dual-polarization SAR image and the second phase dual-polarization SAR image can be calculated according to the maximum eigenvalue and the minimum eigenvalue, and finally a difference map can be obtained according to the similarity measure.

[0049] Exemplarily, the similarity measure SIM may be:

[0050]

[0051] Among them, λ max , min are the maximum eigenvalue and minimum eigenvalue of the similarity matrix respectively.

[0052] In operation S103, the difference map is initially segmented to obtain an initial label of each pixel in the difference map, wherein the initial label is used to distinguish a changed class and an unchanged class in the difference map.

[0053] In the embodiment of the present disclosure, the Kittler-Illingworth automatic threshold method based on the generalized Gaussian distribution can be used to perform initial segmentation on the difference map. Specifically, the following can be done: first, based on the generalized Gaussian distribution model, the change class ω in the difference map is segmented. c and unchanged class ω u Perform statistical modeling to obtain the change class ω c and unchanged class ω u Based on the statistical characteristics, the Kittler-Illingworth threshold selection criterion is used to determine the threshold to segment the difference map and obtain the initial label of each pixel p. Since the generalized Gaussian distribution has only one more shape parameter β than the Gaussian distribution, by adjusting the shape parameter β, the generalized Gaussian distribution can approximate a variety of distributions, such as when β = 2, it is approximately Gaussian distribution; when β = 1, it is approximately Laplace distribution. Therefore, the generalized Gaussian distribution can better fit the statistical characteristics of the changed and unchanged classes in the difference map, and the threshold is automatically determined by combining the Kittler-Illingworth threshold selection criterion, which can achieve a better initial labeling of the changed and unchanged regions in the difference map.

[0054] The model of the generalized Gaussian distribution can be defined as:

[0055]

[0056]

[0057]

[0058] Among them, X l is the grayscale of the grayscale image of the difference map, X l =0, 1, 2, ..., L-1, L is the maximum pixel value of the grayscale image, m i is the mean value of pixels in the changed or unchanged classes in the difference map, β i is the shape parameter of the generalized Gaussian distribution model, β i ≥0,σ i is the standard deviation of the changed or unchanged class in the difference map, Γ(·) is the gamma function, ω i is the changing class or is the changing class, i represents the number of the changing class or the unchanged class, c represents the changing class, and u represents the unchanged class.

[0059] The Kittler-Illingworth threshold selection criterion is to minimize the cost function J(T), determine the threshold T in the grayscale range [0, L-1], and obtain the corresponding m of the changed class or unchanged class according to the determined threshold. i , σ i and β i Specific examples include:

[0060]

[0061]

[0062]

[0063] H c (T) = 1-H u (T)

[0064]

[0065]

[0066]

[0067]

[0068] Where h(·) represents the grayscale histogram of the difference map.

[0069] In operation S104, according to the initial labels, the initial labels of the changed class and the unchanged class in the difference map are optimized based on the generalized Gaussian distribution and the contextual relationship between each pixel and the nearest pixel in the eight-neighborhood direction to obtain a changed binary map.

[0070] In the disclosed embodiment, the energy function of the Markov random field model of the difference map can be constructed based on the initial label, the generalized Gaussian distribution and the contextual relationship between each pixel and the nearest pixel in the eight-neighborhood direction. The minimized energy function of the Markov random field model is solved to obtain the optimal label of each pixel in the difference map. The change binary map is output according to the optimal label.

[0071] Specifically, the Markov random field model E is defined as the data energy E d and smoothing energy E s The linear combination of data energy E d The generalized Gaussian distribution is used to characterize the statistical characteristics of the changed and unchanged classes in the difference map. The smoothing energy Es uses an eight-neighborhood spatial structure to describe the contextual dependency of each pixel p in the difference map with the nearest pixel q in the eight-neighborhood direction, namely:

[0072]

[0073] D p (f p )=ln(G(X l |f p )),

[0074]

[0075] Where P represents the set of all pixels in the difference map, p represents the pixel number, q represents the number of the nearest pixel of the p-th pixel in the eight-neighborhood direction, N is the set of the p-th pixel and the nearest pixel, and f p represents the initial label of the p-th pixel, f q represents the initial label of the qth pixel, φ is the smoothing parameter, which is generally 0.9 by default, and λ is the data energy E d and smoothing energy E s The default balance factor is 1.

[0076] The simulated annealing optimization algorithm of the Metropolis criterion can be used to iteratively minimize the energy function of the Markov random field model to obtain the optimal label f′ of each pixel in the difference map p (f′ p ∈{ω c ,ω u}). The specific process can be:

[0077] First, the simulated annealing control parameters are preset, including the initial annealing temperature T0, the iteration termination threshold t and the total number of iterations K.

[0078] Next, a pixel is randomly selected from the difference map and labeled according to the initial difference map. Determine the category label of the current pixel, and the category label of the pixel is the change class ω c In the case of, get an unchanged class ω from the initial tag combination u , based on the energy function of the Markov random field model, using the generalized Gaussian distribution model parameter m i , σ i , β i , and the smoothing parameter Φ, the difference ΔE between the corresponding energies of the two category labels before and after is calculated.

[0079] Next, the Metropolis criterion is used to determine whether the difference meets the preset conditions. If exp(ΔE / T k )≥α, α∈(0,1) is a random number, then the category of the current pixel is marked from the change category ω c Update to unchanged classω u , and calculate the absolute value of the energy difference before and after the update |ΔE|,

[0080] In this way, the above method is used to iterate over other pixels in the difference map, update the category label of each pixel, and calculate the total energy ∑|ΔE| of the entire difference map that satisfies the Metropolis criterion.

[0081] Finally, it is determined whether the iteration meets the termination condition, including whether the total energy ∑|ΔE| is less than the iteration termination threshold t or whether the number of iterations reaches the preset number K. If the current iteration k does not meet the termination condition, the temperature is lowered T k+1 =λT k , k∈[0, K-1], λ defaults to 0.98. Repeat the above process until the iteration termination condition is met, obtain the optimal label of each pixel in the difference map, achieve accurate segmentation of the changed class and the unchanged class, and output the changed binary map.

[0082] In summary, the multi-phase SAR image change detection method provided by the embodiment of the present disclosure constructs a similarity matrix between the covariance matrices of the dual-polarization SAR data of the previous and next phases, and uses the maximum and minimum eigenvalues ​​of the similarity matrix to construct a similarity measure to calculate the difference map between the previous and next phase images, which can effectively highlight the change area between the previous and next phase images; uses the hierarchical image segmentation strategy of "Kittler-Illingworth automatic threshold initial segmentation based on generalized Gaussian distribution → precise segmentation based on Markov random field model", uses the generalized Gaussian distribution to better fit the statistical characteristics of the changed class and the unchanged class in the difference map, and uses the eight-neighborhood spatial structure to well describe the context information of each pixel in the difference map, thereby effectively suppressing the false detection caused by speckle noise, and realizing the automatic and precise segmentation of the changed class and the unchanged class in the difference map, which plays an important role in SAR change detection and is of great significance for promoting the practical application of dual-polarization SAR change detection. It is of great significance for promoting the practical application of dual-polarization SAR change detection.

[0083] The specific embodiments described above further illustrate the purpose, technical solutions and beneficial effects of the present disclosure. It should be understood that the above description is only a specific embodiment of the present disclosure and is not intended to limit the present disclosure. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present disclosure should be included in the protection scope of the present disclosure.

Claims

1. A multi-temporal SAR image change detection method, comprising: Constructing a similarity matrix between a first covariance matrix of a dual-polarization SAR image of a first time phase and a second covariance matrix of a dual-polarization SAR image of a second time phase; Calculating a difference map between the first phase dual-polarization SAR image and the second phase dual-polarization SAR image according to the similarity matrix; Performing an initial segmentation on the difference map to obtain an initial label of each pixel in the difference map, wherein the initial label is used to distinguish a changed class and an unchanged class in the difference map; According to the initial labeling, the initial labeling of the changed class and the unchanged class in the difference map is optimized based on the generalized Gaussian distribution and the contextual relationship between each pixel and the nearest pixel in the eight-neighborhood direction to obtain a changed binary map, including: according to the initial labeling, based on the generalized Gaussian distribution and the contextual relationship between each pixel and the nearest pixel in the eight-neighborhood direction, constructing a Markov random field model energy function of the difference map; solving the minimized Markov random field model energy function to obtain the optimal labeling of each pixel in the difference map, and outputting a changed binary map according to the optimal labeling; The energy function of the Markov random field model is: Wherein, P represents the set of all pixels in the difference map, p represents the pixel number, q represents the number of the nearest pixel of the p-th pixel in the eight-neighborhood direction, N is the set consisting of the p-th pixel and the nearest pixel, represents the initial label of the p-th pixel, represents the initial label of the qth pixel, is the smoothing parameter, is the balance factor; The energy function of the Markov random field model to be minimized includes: Randomly select a pixel from the difference map, and when the class of the pixel is marked as a changed class, obtain an unchanged class according to the initial mark; Calculating the difference between the energy corresponding to the initial mark of the changed class and the energy corresponding to the initial mark of the unchanged class based on the Markov random field model energy function; If the difference satisfies a preset condition, the category label of the current pixel is updated from the changed category to the unchanged category; In this way, it is iterated to traverse other pixels in the difference map to obtain the sum of the energy differences corresponding to the difference map, until the sum of the energy differences is less than a preset threshold or the number of iterations reaches a preset number, and the optimal labeling of each pixel in the difference map is obtained.

2. The multi-temporal SAR image change detection method according to claim 1, wherein: The process of constructing a similarity matrix between the covariance matrix of the dual-polarization SAR image of the first phase and the covariance matrix of the dual-polarization SAR image of the second phase includes: Constructing the first covariance matrix according to the backscattering information of the co-polarization channel and the backscattering information of the cross-polarization channel corresponding to the first phase dual-polarization SAR image; Constructing the second covariance matrix according to the backscattering information of the co-polarization channel and the backscattering information of the cross-polarization channel corresponding to the second phase dual-polarization SAR image; The similarity matrix is ​​constructed according to the first covariance matrix and the second covariance matrix.

3. The multi-temporal SAR image change detection method according to claim 2, wherein: The similarity matrix for: ， in, is the backscattering information of the co-polarization channel, is the backscattering information of the cross-polarization channel, * indicates the complex conjugate calculation; Indicates that the window is averaged. Represents complex modulus calculation, represents the conjugate transpose of a matrix, is the first covariance matrix, is the second covariance matrix.

4. The multi-temporal SAR image change detection method according to claim 1, wherein: The calculating the difference map between the first phase dual-polarization SAR image and the second phase dual-polarization SAR image according to the similarity matrix comprises: Solving the maximum eigenvalue and the minimum eigenvalue of the similarity matrix; Calculating a similarity measure between the first phase dual-polarization SAR image and the second phase dual-polarization SAR image according to the maximum eigenvalue and the minimum eigenvalue; The difference map is obtained according to the similarity measure.

5. The multi-temporal SAR image change detection method according to claim 4, wherein: The similarity measure SIM is: in, , are the maximum eigenvalue and minimum eigenvalue of the similarity matrix respectively.

6. The multi-temporal SAR image change detection method according to claim 1, wherein: The initial segmentation of the difference map to obtain the initial label of each pixel in the difference map includes: Performing statistical modeling on the changed class and the unchanged class in the difference graph based on a generalized Gaussian distribution model to obtain statistical characteristics of the changed class and the unchanged class; Based on the statistical characteristics, the Kittler–Illingworth threshold selection criterion is used to determine the threshold to segment the difference map, thereby obtaining the initial label of each pixel.

7. The multi-temporal SAR image change detection method according to claim 6, wherein: according to The generalized Gaussian distribution model is constructed, wherein, is the grayscale of the grayscale image of the difference map, =0,1,2,…,L-1, L is the maximum pixel value of the grayscale image, is the mean value of pixels in the changed or unchanged class in the difference map, is the shape parameter of the generalized Gaussian distribution model, is the standard deviation of the changed or unchanged class in the difference map, is the gamma function, is a changed class or is a changed class, i represents the number of the changed class or unchanged class, c represents the changed class, and u represents the unchanged class; The Kittler–Illingworth threshold selection criterion is to minimize the cost function. Determine a threshold value, and obtain the , and .

Images