Remote sensing image change detection method based on double non-negative matrix factorization combined with sparse constraints

Abstract

The present invention discloses a remote sensing image change detection method based on double non-negative matrix decomposition with joint sparse constraints, which belongs to the field of image change detection technology and aims to solve the problems that the existing image change detection methods lack consideration of semantic information of ground objects, are sensitive to noise, and have low processing efficiency. It comprises: S1: based on the first phase remote sensing image and the second phase remote sensing image taken at different times in the same geographical area, using L 2,1 DNMF model, generating a differentially enhanced differential feature matrix; S2: clustering the differentially enhanced differential feature matrix to obtain a binary change result; S3: outputting the binary change result as a remote sensing image change detection result of a joint sparse constrained dual non-negative matrix decomposition. The present invention can detect subtle changes at the sub-pixel level, thereby improving the change detection performance and change detection efficiency.

Description

Technical Field

[0001] The present invention relates to the technical field of image change detection, and in particular to a remote sensing image change detection method based on dual non-negative matrix decomposition with joint sparse constraints. Background Art

[0002] With the rapid development of remote sensing technology, remote sensing images are more conducive to our acquisition of richer information about land objects. Accurate analysis of this information is of great significance for understanding the impact of rapid population growth and urbanization on resources and the environment. Change detection is an effective tool for observing surface changes. It identifies regional surface differences by comparing a set of images acquired at different times under the same scene. Hyperspectral images contain rich spectral information and can be used to distinguish spectrally similar objects, making them the main data source for change detection. Hyperspectral image change detection has become one of the most active research topics in remote sensing image processing and has been widely used in many application areas, such as land cover change, environmental monitoring, and natural disaster assessment.

[0003] In the past few decades, many excellent ideas have been applied to change detection and achieved good results. However, there are still some defects in practical applications. Some of these methods directly use differential calculation to obtain change vectors. On the differential results, algebraic calculations or feature extraction methods are used to extract difference features to represent the intensity of the change. Although these methods can quickly obtain results, they usually only focus on the differences between bands and lack consideration of the semantic information of the ground objects. In addition, they are sensitive to noise and require high-quality input data. These defects will limit the effect of change detection. There are also some methods that process the original dual-phase images separately and use the processing results to generate change vectors for change detection. For example, classification-based methods and unmixing-based methods. These change detection methods have been shown to perform well in terms of detection effects. However, hyperspectral remote sensing images have high resolution and contain rich spectral information. This type of method will take longer to process and has higher requirements on computers, requiring computers to be able to process more data.

[0004] Therefore, the change detection method for hyperspectral remote sensing images needs to consider how to improve the efficiency of processing and analysis while maintaining good detection results. Summary of the invention

[0005] The purpose of the present invention is to provide a remote sensing image change detection method based on dual non-negative matrix decomposition with joint sparse constraints, so as to solve the problems that the existing image change detection methods lack consideration of semantic information of ground objects, are sensitive to noise, and have low processing efficiency.

[0006] The technical solution of the present invention to solve the above technical problems is as follows:

[0007] The present invention provides a remote sensing image change detection method using dual non-negative matrix decomposition with joint sparse constraints, comprising:

[0008] S1: Based on the first phase remote sensing image and the second phase remote sensing image taken at different times in the same geographical area, L 2,1 -DNMF model, generating differentially enhanced difference feature matrix;

[0009] S2: clustering the difference feature matrix of the difference enhancement to obtain a binary change result;

[0010] S3: Output the binary change result as a remote sensing image change detection result of a joint sparse constrained dual non-negative matrix decomposition.

[0011] Optionally, the S1 includes:

[0012] S11: Expanding the first temporal phase remote sensing image and the second temporal phase remote sensing image into a first temporal phase observation matrix and a second temporal phase observation matrix respectively;

[0013] S12: Generate L according to the first phase measurement matrix and the second phase measurement matrix 2,1 -Objective function and constraints of DNMF model;

[0014] S13: Using the L 2,1 -The objective function and constraint conditions of the DNMF model decompose the first phase observation matrix and the second phase observation matrix into a shared basis matrix, a first multi-feature weight matrix and a second multi-feature weight matrix, and obtain the original differential enhanced difference feature matrix according to the first multi-feature weight matrix and the second multi-feature weight matrix;

[0015] S14: according to the objective function and the constraint conditions, using the augmented Lagrange multiplier method to solve the objective function and the constraint conditions, and obtain an update rule for the shared basis matrix, the first multi-feature weight matrix, the second multi-feature weight matrix and the original difference enhanced difference feature matrix;

[0016] S15: Initializing variables using the NNDSVD method to obtain an initial shared basis matrix, an initial first multi-feature weight matrix, an initial second multi-feature weight matrix, and an initial differentially enhanced difference feature matrix;

[0017] S16: using the update rule of the shared base matrix, the update rule of the first multi-feature weight matrix, the update rule of the second multi-feature weight matrix and the update rule of the differentially enhanced difference feature matrix to update the current shared base matrix, the current first multi-feature weight matrix, the current second multi-feature weight matrix and the current differentially enhanced difference feature matrix respectively, to obtain a new shared base matrix, a new first multi-feature weight matrix, a new second multi-feature weight matrix and a new differentially enhanced difference feature matrix;

[0018] S17: Use the new shared base matrix, the new first multi-feature weight matrix, the new second multi-feature weight matrix and the new differentially enhanced difference feature matrix as the current shared base matrix, the current first multi-feature weight matrix, the current second multi-feature weight matrix and the current differentially enhanced difference feature matrix and return to S16 for iterative update until the cutoff condition is met to obtain the differentially enhanced difference feature matrix.

[0019] Optionally, in step S11, the first time phase observation matrix and / or the second time phase observation matrix is ​​a two-dimensional matrix, and the two dimensions of the two-dimensional matrix are the number of bands and the number of pixels of the image, respectively.

[0020] Optionally, in S12, the L 2,1 The objective function of the DNMF model is:

[0021]

[0022] The constraints are:

[0023]

[0024] Where X is the first phase observation matrix and Y is the second phase observation matrix and N = I × J represents the number of pixels in the image; I, J and L correspond to the number of rows, columns and bands of the image respectively; A is the shared basis matrix of X and Y and H1 is the first multi-feature weight matrix; H2 is the second multi-feature weight matrix, C is the difference feature matrix of difference enhancement; ‖·‖ F is the Frobenius norm, ‖·‖ 2,1 For L 2,1 -norm; λ is the regularization parameter; (·) T Represents the transpose of the matrix; 1 R represents a column vector of all 1s with dimension R, 1 N represents a column vector of all ones with dimension N.

[0025] The S14 includes:

[0026] S141: Generate L according to the objective function and the constraint conditions 2,1 -Augmented Lagrangian function of DNMF model;

[0027] S142: Fix the L 2,1 The other variables in the augmented Lagrangian function of the -DNMF model are solved for the remaining single variables to obtain the update rules of the shared basis matrix, the update rules of the first multi-feature weight matrix, the update rules of the second multi-feature weight matrix and the update rules of the differential enhanced difference feature matrix.

[0028] Optionally, the L 2,1 The augmented Lagrangian function of the -DNMF model is:

[0029]

[0030] Where L represents the augmented Lagrangian function, A represents the shared basis matrix, H1 represents the first multi-feature weight matrix, H2 represents the second multi-feature weight matrix, M1, M2, M3 and M4 are Lagrangian multipliers; μ is the penalty parameter and μ ≥ 0; X is the first phase observation matrix and Y is the second phase observation matrix and N = I × J represents the number of pixels in the image; I, J and L correspond to the number of rows, columns and bands of the image respectively; C is the difference feature matrix of differential enhancement; ‖·‖ F is the Frobenius norm, ‖·‖ 2,1 For L 2,1 -norm; λ is the regularization parameter; tr(·) represents the trace of the matrix.

[0031] Optionally, the update rule of the shared basis matrix includes:

[0032]

[0033] Among them, A i+1 represents the shared basis matrix of the i+1th iteration, A i represents the shared basis matrix of the i-th iteration, * represents the multiplication symbol, X is the first phase observation matrix and Y is the second phase observation matrix and N=I×J represents the number of pixels in the image; I, J and L correspond to the number of rows, columns and bands of the image respectively; H1 represents the first multi-feature weight matrix, and H2 represents the second multi-feature weight matrix;

[0034] The updating rules of the first multi-feature weight matrix H1 include:

[0035]

[0036] in, represents the first multi-feature weight matrix of the i+1th iteration, represents the first multi-feature weight matrix of the i-th iteration, represents the augmented matrix of the first phase observation matrix X and represents the augmented matrix of the shared basis matrix A and C is the difference feature matrix of differential enhancement; μ is the penalty parameter, H2 represents the second multi-feature weight matrix, and δ is the parameter used to control the influence of normalization constraints; 1 R represents a column vector of all 1s with dimension R, 1 N represents a column vector of all 1s with dimension N;

[0037] The update rules of the second multi-feature weight matrix H2 include:

[0038]

[0039] in, represents the augmented matrix of the second observation matrix and

[0040] Optionally, the S142 further includes:

[0041] Using the L 2,1 -The augmented Lagrangian function of the DNMF model solves the objective function to obtain the update rule of the difference feature matrix enhanced by the difference, the update rule of the Lagrangian multiplier and the update rule of the penalty parameter;

[0042] The fast update rule of the difference feature matrix of the difference enhancement includes:

[0043]

[0044] Among them, C i+1 represents the difference feature matrix of the difference enhancement of the i+1th iteration, Ω represents the shrinkage operator, λ is the regularization parameter; μ is the penalty parameter; H1 represents the first multi-feature weight matrix, H2 represents the second multi-feature weight matrix; M4 represents the Lagrange multiplier;

[0045] The update rule of Lagrange multiplier M4 is:

[0046] M4=M4+μ(C-H1+H2)

[0047] The update rule of penalty parameter μ is:

[0048] μ=min(ρμ,maxμ)

[0049] Here, ρ is the learning rate used for parameter updating.

[0050] Optionally, the S15 includes:

[0051] S151: horizontally splicing the first temporal phase remote sensing image and the second temporal phase remote sensing image to obtain a spliced ​​three-dimensional image;

[0052] S152: Expanding the spliced ​​three-dimensional image into a two-dimensional matrix;

[0053] S153: performing singular value decomposition on the two-dimensional matrix to obtain a left singular value matrix, a singular value matrix and a right singular value matrix of the matrix;

[0054] S154: Take the largest target number of singular values ​​and the corresponding left singular value matrix and right singular value matrix;

[0055] S155: determining an initial shared basis matrix according to the product of the square root of the singular value and the corresponding column in the corresponding left singular value matrix, and determining the initial first multi-feature weight matrix and the initial second multi-feature weight matrix according to the product of the square root of the singular value and the corresponding column in the corresponding right singular value matrix;

[0056] S156: Subtract the initial first multi-feature weight matrix from the initial second multi-feature weight matrix to obtain an initial differentially enhanced difference feature matrix.

[0057] Optionally, in the step S2, a multi-layer k-means clustering method is used to cluster the difference feature matrix of the difference enhancement;

[0058] The multi-layer k-means clustering algorithm includes the following sub-steps:

[0059] S21: Select 3 pixels as the initial cluster centers of the change class, the unchanged class, and the uncertain class;

[0060] S22: Calculate the Euclidean distance between the pixel point and each of the initial cluster centers, and assign the pixel point to the cluster center closest to it;

[0061] S23: Repeat step S22 until all pixels are allocated;

[0062] S24: updating the position of one of the cluster centers, and moving the cluster center to the position of the data mean belonging to the cluster;

[0063] S25: Repeat step S24 until the positions of all the cluster centers are updated;

[0064] S26: Repeat steps S22-S25 until the position of the cluster center no longer changes or reaches a predetermined number of iterations.

[0065] S27: Repeat steps S21-S26 for the pixels in the uncertain class until the uncertain class cannot be divided any further, so as to obtain a binary change result.

[0066] The present invention has the following beneficial effects:

[0067] 1. The NMF used in the present invention extracts feature weights in units of pixels and can detect subtle changes at the sub-pixel level;

[0068] 2. H1 and H2 in the present invention are multi-feature weight matrices learned based on the common features of the two images, and C is an enhanced representation of the difference calculation of H1 and H2. This process makes full use of the relevant feature information between the two-phase images, which can further improve the detection performance;

[0069] 3. The L having a row thinning promoting effect in the present invention 2,1 The norm constraint makes the value of the unchanged feature in C as close to zero as possible, which improves the separability of the changed and unchanged features and enables clustering to obtain more accurate binary change results;

[0070] 4. The multi-layer k-means clustering used by the present invention to generate binary change results is an unsupervised classification method that does not require labeled training data and is applicable to unlabeled data sets. In addition, the time complexity of the k-means clustering algorithm is linear, so it performs well on large-scale data sets and can help improve the efficiency of change detection;

[0071] 5. The objective function of the present invention can simultaneously extract relevant information of the dual-phase image to obtain a difference feature map, and uses a fast solution algorithm with a fast convergence speed, which greatly improves the detection efficiency and has great advantages in processing large-scale data. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 It is a flow chart of the remote sensing image change detection method of the present invention with dual non-negative matrix decomposition combined with sparse constraints;

[0073] Figure 2 This is a flowchart of the joint sparse constrained dual non-negative matrix factorization remote sensing image change detection method. DETAILED DESCRIPTION

[0074] The principles and features of the present invention are described below in conjunction with the accompanying drawings. The examples given are only used to explain the present invention and are not used to limit the scope of the present invention.

[0075] The present invention provides a joint sparse constrained dual non-negative matrix decomposition remote sensing image change detection method; the present invention uses a non-negative matrix (NMF) to learn a non-negative basis matrix of a dual-phase hyperspectral image, and enables them to share the basis matrix. When the basis matrices are consistent, two non-negative multi-feature weight matrices can be learned based on the basis matrix, and their changes can be obtained by differential calculation. Therefore, the difference information can be directly extracted from the two images without processing the original images separately, thereby effectively improving the detection efficiency. In the differential matrix, the difference of the unchanged features is not equal to zero, and the separability of the changed and unchanged features is not particularly good. In order to obtain a more meaningful sparse differential feature map, the present invention uses L with row sparsity promotion effect. 2,1 Norm regularized double non-negative matrix factorization (L 2,1 -DNMF) model can efficiently learn the difference information of dual-phase hyperspectral images and complete the change detection task.

[0076] In practical applications, hyperspectral remote sensing image data usually has a large number of spectral bands, contains rich information, and has a large amount of data. Change detection requires processing multiple data, takes a long time to calculate, and consumes a lot of computing resources. 2,1 -DNMF model is used for change detection, which can complete the detection task in a shorter time, saving computing resources and time costs.

[0077] The present invention provides a remote sensing image change detection method based on double non-negative matrix decomposition with joint sparse constraints. Figure 1 As shown, including:

[0078] S1: Based on the first phase remote sensing image and the second phase remote sensing image taken at different times in the same geographical area, L 2,1 -DNMF model, generating differentially enhanced difference feature matrix;

[0079] The first phase remote sensing image and the second phase remote sensing image are three-dimensional dual-phase data. The three dimensions of the three-dimensional dual-phase data provided by the present invention are the number of rows, the number of columns and the number of bands of the image. The two images are taken in the same geographical area but at different times. Here, X represents the image before the change, and Y represents the image after the change.

[0080] Optionally, the S1 includes:

[0081] S11: Expand the first temporal phase remote sensing image and the second temporal phase remote sensing image into a first temporal phase observation matrix and a second temporal phase observation matrix respectively; that is, expand the I×J×L three-dimensional dual-temporal phase hyperspectral data into a first temporal phase observation matrix and the second phase observation matrix Where I, J and L correspond to the number of rows, columns and spectral dimension (representing the number of bands of the hyperspectral image) of the image respectively, and N = I × J represents the number of pixels in the image.

[0082] S12: Generate L according to the first phase measurement matrix and the second phase measurement matrix 2,1 -Objective function and constraints of DNMF model;

[0083] L 2,1 The objective function of the DNMF model is:

[0084]

[0085] The constraints are:

[0086]

[0087] Where X is the first phase observation matrix and Y is the second phase observation matrix and N = I × J represents the number of pixels in the image; I, J and L correspond to the number of rows, columns and bands of the image respectively; A is the shared basis matrix of X and Y and H1 is the first multi-feature weight matrix; H2 is the second multi-feature weight matrix, C is the difference feature matrix of difference enhancement; ‖·‖ F is the Frobenius norm, ‖·‖ 2,1 For L 2,1 -norm; λ is the regularization parameter; (·) T Represents the transpose of the matrix; 1 R represents a column vector of all 1s with dimension R, 1 N represents a column vector of all ones with dimension N.

[0088] S13: Using the L 2,1 -The objective function and constraint conditions of the DNMF model decompose the first phase observation matrix and the second phase observation matrix into a shared basis matrix, a first multi-feature weight matrix and a second multi-feature weight matrix, and obtain the original differential enhanced difference feature matrix according to the first multi-feature weight matrix and the second multi-feature weight matrix;

[0089] S14: according to the objective function and the constraint conditions, using the augmented Lagrange multiplier method to solve the objective function and the constraint conditions, and obtain an update rule for the shared basis matrix, the first multi-feature weight matrix, the second multi-feature weight matrix and the original difference enhanced difference feature matrix;

[0090] Optionally, S14 includes:

[0091] S141: Generate L according to the objective function and the constraint conditions 2,1 -Augmented Lagrangian function of DNMF model;

[0092] L 2,1 The augmented Lagrangian function of the -DNMF model is:

[0093]

[0094] Where L represents the augmented Lagrangian function, A represents the shared basis matrix, H1 represents the first multi-feature weight matrix, H2 represents the second multi-feature weight matrix, M1, M2, M3 and M4 are Lagrangian multipliers; μ is the penalty parameter and μ ≥ 0; X is the first phase observation matrix and Y is the second phase observation matrix and N = I × J represents the number of pixels in the image; I, J and L correspond to the number of rows, columns and bands of the image respectively; C is the difference feature matrix of differential enhancement; ‖·‖ F is the Frobenius norm, ‖·‖ 2,1 For L 2,1 -norm; λ is the regularization parameter; tr(·) represents the trace of the matrix.

[0095] S142: Fix the L 21 The other variables in the augmented Lagrangian function of the -DNMF model are solved for the remaining single variables to obtain the update rules of the shared basis matrix, the update rules of the first multi-feature weight matrix, the update rules of the second multi-feature weight matrix and the update rules of the differential enhanced difference feature matrix.

[0096] The update rules for the shared basis matrix include:

[0097] When other variables are fixed, the sub-optimization problem of A is:

[0098]

[0099] can be rewritten as:

[0100]

[0101] The partial derivative of F with respect to A is:

[0102]

[0103] According to the KKT condition, that is, and M1.*A=0 is applied to Generate update rules for A;

[0104] That is, update using the following formula:

[0105]

[0106] Among them, A i+1 represents the shared basis matrix of the i+1th iteration, A i represents the shared basis matrix of the i-th iteration, * represents the multiplication symbol, X is the first phase observation matrix and Y is the second phase observation matrix and N=I×J represents the number of pixels in the image; I, J and L correspond to the number of rows, columns and bands of the image respectively; H1 represents the first multi-feature weight matrix, and H2 represents the second multi-feature weight matrix;

[0107] Similar to the solution of A, the sub-optimization problems of H1 and H2 are:

[0108]

[0109]

[0110] In order to satisfy the normalization constraint and Change the above two formulas into:

[0111]

[0112]

[0113] in δ is a parameter used to control the influence of the normalization constraint. Then, the update rules for H1 and H2 are:

[0114] The update rules of the first multi-feature weight matrix H1 include:

[0115]

[0116] in, represents the first multi-feature weight matrix of the i+1th iteration, represents the first multi-feature weight matrix of the i-th iteration, represents the augmented matrix of the first phase observation matrix X and represents the augmented matrix of the shared basis matrix A and C is the difference feature matrix of differential enhancement; μ is the penalty parameter, H2 represents the second multi-feature weight matrix, and δ is the parameter used to control the influence of normalization constraints; 1 R represents a column vector of all 1s with dimension R, 1 N represents a column vector of all ones with dimension N.

[0117] The update rules of the second multi-feature weight matrix H2 include:

[0118]

[0119] in, represents the augmented matrix of the second observation matrix and

[0120] The A obtained in this step is the shared base matrix of the dual-phase image extracted by double non-negative matrix decomposition. This base matrix combines the common feature information of the two images, and the decomposition dimension R is the number of extracted features. Each column of H1 and H2 represents the weight coefficient of R features at a pixel. By performing differential calculations on the two matrices H1 and H2, a feature-based change vector can be obtained.

[0121] In addition, S142 of the present invention further includes:

[0122] Using the L 2,1 The augmented Lagrangian function of the -DNMF model obtains the subproblem of solving the difference feature matrix of differential enhancement:

[0123]

[0124] C is the difference feature matrix of difference enhancement; ‖·‖ F is the Frobenius norm, ‖·‖ 2,1 For L 2,1 -norm; λ is the regularization parameter; H1 is the first multi-feature weight matrix; H2 is the second multi-feature weight matrix, μ is the penalty parameter and μ ≥ 0; M4 is the Lagrange multiplier;

[0125] Then we obtain the fast update rule of the difference feature matrix of the difference enhancement, the update rule of the Lagrange multiplier and the update rule of the penalty parameter;

[0126] The updating rules of the difference feature matrix of the difference enhancement include:

[0127]

[0128] Among them, C i+1 represents the difference feature matrix of the difference enhancement of the i+1th iteration, Ω represents the shrinkage operator, λ is the regularization parameter; μ is the penalty parameter; H1 represents the first multi-feature weight matrix, H2 represents the second multi-feature weight matrix; M4 represents the Lagrange multiplier;

[0129] The update rule of Lagrange multiplier M4 is:

[0130] M4=M4+μ(C-H1+H2)

[0131] The update rule of penalty parameter μ is:

[0132] μ=min(ρμ,maxμ)

[0133] Where ρ is the learning rate for parameter update. The differential feature matrix C of differential enhancement can represent the change intensity of each feature at each pixel. A value greater than 0 indicates a positive change, and a value less than 0 indicates a negative change.

[0134] S15: Initializing variables using the NNDSVD method to obtain an initial shared basis matrix, an initial first multi-feature weight matrix, an initial second multi-feature weight matrix, and an initial differentially enhanced difference feature matrix;

[0135] The non-negative matrix factorization algorithm is a non-convex optimization process in the iterative process, and its result depends on the initial value and is not unique. The selection of the initial values ​​A, H1 and A, H2 will directly affect the iterative results of the decomposition algorithm.

[0136] There are many methods for initializing NMF, among which the most popular ones are based on singular value decomposition (SVD) and clustering. Common non-negative matrix factorization initialization methods include: non-negative double singular value decomposition initialization (NNDSVD), K-means clustering initialization (K-means), fuzzy C-means clustering initialization (FCM), etc. The present invention uses the NNDSVD initialization method to determine the initial values ​​of A, H1 and H2, which effectively solves the problem that the NMF algorithm converges slowly and is easy to fall into local minima.

[0137] Based on this, reference Figure 2 As shown, optionally, the S15 includes:

[0138] S151: horizontally splicing the first temporal phase remote sensing image and the second temporal phase remote sensing image to obtain a spliced ​​three-dimensional matrix;

[0139] S152: Expand the spliced ​​three-dimensional matrix into a two-dimensional matrix

[0140] S153: Perform singular value decomposition (SVD) on the two-dimensional matrix, which can be expressed as Z = U∑V T ; Get the left singular value matrix, singular value matrix and right singular value matrix of the matrix; where the left singular value matrix and the right singular value matrix They are all orthogonal matrices. The singular value matrix ∑ is an L×2*N dimensional matrix, whose r×r dimensional diagonal elements are the singular values ​​of the matrix Z, usually arranged in descending order, and r is the rank of the matrix.

[0141] S154: Take the R largest singular values and the corresponding left singular value matrix and the right singular value matrix where R < min(L, N) is the dimension of non - negative matrix factorization.

[0142] S155: Determine the shared basis matrix according to the product of the square root of the singular value and the corresponding column in the U R and determine the first multi - feature weight matrix and the second multi - feature weight matrix according to the product of the square root of the singular value and the corresponding column in the V R Set the first column of W as:

[0143] Similarly, set the first row of H as:

[0144]

[0145] For the next singular vectors: w2 and h2, due to the orthogonality of singular value decomposition, w2 and h2 may contain negative elements. For the matrix M = σ2w2h2

[0146]

[0147] , replace the negative elements in M with 0, and keep the rest unchanged, that is, take the positive part M T of the matrix M. M + is used as an approximation of M. Perform singular value decomposition on M + , and use the first pair of singular values of M + as W(:, 2) and H(2, :), and perform similar processing for the initialization of other columns and rows of W and H later. + W is the initial value of the present invention A. Take the first N columns of H as the initial value of H1, and the last N columns of H as the initial value of H2.

[0148] S156: Subtract the initial first multi - feature weight matrix from the initial second multi - feature weight matrix to obtain the initial differential enhanced difference feature matrix.

[0149] The initial value of the matrix C is defined as the difference between the initial value of H1 and the initial value of H2.

[0150]

[0151] ​S16: using the update rule of the shared base matrix, the update rule of the first multi-feature weight matrix, the update rule of the second multi-feature weight matrix and the update rule of the differentially enhanced difference feature matrix to update the current shared base matrix, the current first multi-feature weight matrix, the current second multi-feature weight matrix and the current differentially enhanced difference feature matrix respectively, to obtain a new shared base matrix, a new first multi-feature weight matrix, a new second multi-feature weight matrix and a new differentially enhanced difference feature matrix;

[0152] S17: Use the new shared base matrix, the new first multi-feature weight matrix, the new second multi-feature weight matrix and the new differentially enhanced difference feature matrix as the current shared base matrix, the current first multi-feature weight matrix, the current second multi-feature weight matrix and the current differentially enhanced difference feature matrix and return to S16 for iterative update until the cutoff condition is met to obtain the differentially enhanced difference feature matrix.

[0153] S2: clustering the difference feature matrix of the difference enhancement to obtain a binary change result;

[0154] That is, in order to produce the final binary change result, the multi-layer clustering method is used to divide the pixels into two classes on the difference feature matrix obtained by differential enhancement. This process is unsupervised and does not require any manual labeling. More accurate change detection results can be obtained through multiple clustering.

[0155] The present invention utilizes a multi-layer k-means clustering method to cluster the difference feature matrix of the difference enhancement;

[0156] The multi-layer k-means clustering algorithm consists of the following sub-steps:

[0157] S21: Select 3 pixels as the initial cluster centers of the change class, the unchanged class, and the uncertain class;

[0158] S22: Calculate the Euclidean distance between the pixel point and each of the initial cluster centers, and assign the pixel point to the cluster center closest to it;

[0159] The calculation formula of Euclidean distance is:

[0160]

[0161] Wherein, d is the distance, c is the pixel, μ is the cluster center, n is the number of features in the pixel, and i1 is each feature in the pixel;

[0162] S23: Repeat step S22 until all pixels are allocated;

[0163] S24: updating the position of one of the cluster centers, and moving the cluster center to the position of the data mean belonging to the cluster;

[0164] S25: Repeat step S24 until the positions of all the cluster centers are updated;

[0165] S26: Repeat steps S22-S25 until the position of the cluster center no longer changes or reaches a predetermined number of iterations, at which time C is clustered into three categories, representing the unchanged class (represented by "0"), the changed class (represented by "1") and the uncertain class (represented by "2").

[0166] S27: Repeat steps S21-S26 for the pixels in the uncertain class, mark the unchanged class as 0, the changed class as 1, and the uncertain class as 2. This is done until the uncertain class cannot be divided any further, so as to obtain a binary change result.

[0167] S3: Output the binary change result as a remote sensing image change detection result of a joint sparse constrained dual non-negative matrix decomposition.

[0168] The present invention has the following advantages:

[0169] 1. The NMF used in the present invention extracts feature weights in units of pixels and can detect subtle changes at the sub-pixel level.

[0170] 2. H1 and H2 in the present invention are multi-feature weight matrices learned based on the common features of the two images, and C is an enhanced representation of the difference calculation of H1 and H2. This process makes full use of the relevant feature information between the two-phase images to further improve the detection performance.

[0171] 3. The L having a row thinning promoting effect in the present invention 2,1 The -norm constraint makes the value of the unchanged feature in C as close to zero as possible, which improves the separability of the changed and unchanged features and enables clustering to obtain more accurate binary change results.

[0172] 4. The multi-layer k-means clustering used by the present invention to generate binary change results is an unsupervised classification method that does not require labeled training data and is applicable to unlabeled data sets. In addition, the time complexity of the k-means clustering algorithm is linear, so it performs well on large-scale data sets and can help improve the efficiency of change detection.

[0173] 5. The objective function of the present invention can simultaneously extract relevant information of the dual-phase image to obtain a difference feature map, and uses a fast solution algorithm with a fast convergence speed, which greatly improves the detection efficiency and has great advantages in processing large-scale data.

[0174] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A remote sensing image change detection method based on dual non-negative matrix factorization with joint sparse constraints, characterized in that: include: S1: Based on the first phase remote sensing image and the second phase remote sensing image taken at different times in the same geographical area, L 2,1 -DNMF model, generating differentially enhanced difference feature matrix; S2: clustering the difference feature matrix of the difference enhancement to obtain a binary change result; S3: Outputting the binary change result as a remote sensing image change detection result of a joint sparse constrained dual non-negative matrix decomposition; The S1 includes: S11: Expanding the first temporal phase remote sensing image and the second temporal phase remote sensing image into a first temporal phase observation matrix and a second temporal phase observation matrix respectively; S12: Generate L according to the first phase measurement matrix and the second phase measurement matrix 2,1 -Objective function and constraints of DNMF model; S13: Using the L 2,1 -The objective function and constraint conditions of the DNMF model decompose the first phase observation matrix and the second phase observation matrix into a shared basis matrix, a first multi-feature weight matrix and a second multi-feature weight matrix, and obtain the original differential enhanced difference feature matrix according to the first multi-feature weight matrix and the second multi-feature weight matrix; S14: according to the objective function and the constraint conditions, using the augmented Lagrange multiplier method to solve the objective function and the constraint conditions, and obtain an update rule for the shared basis matrix, the first multi-feature weight matrix, the second multi-feature weight matrix and the original difference enhanced difference feature matrix; S15: Initializing variables using the NNDSVD method to obtain an initial shared basis matrix, an initial first multi-feature weight matrix, an initial second multi-feature weight matrix, and an initial differentially enhanced difference feature matrix; S16: using the update rule of the shared base matrix, the update rule of the first multi-feature weight matrix, the update rule of the second multi-feature weight matrix and the update rule of the differentially enhanced difference feature matrix to update the current shared base matrix, the current first multi-feature weight matrix, the current second multi-feature weight matrix and the current differentially enhanced difference feature matrix respectively, to obtain a new shared base matrix, a new first multi-feature weight matrix, a new second multi-feature weight matrix and a new differentially enhanced difference feature matrix; S17: Using the new shared base matrix, the new first multi-feature weight matrix, the new second multi-feature weight matrix and the new differentially enhanced difference feature matrix as the current shared base matrix, the current first multi-feature weight matrix, the current second multi-feature weight matrix and the current differentially enhanced difference feature matrix, and returning to S16 for iterative update until the cutoff condition is met, to obtain the differentially enhanced difference feature matrix; In the step S11, the first time phase measurement matrix and / or the second time phase measurement matrix is ​​a two-dimensional matrix, and the two dimensions of the two-dimensional matrix are the number of bands and the number of pixels of the image respectively; In S12, the L 2,1 The objective function of the DNMF model is: The constraints are: Where X is the first phase observation matrix and Y is the second phase observation matrix and N = I × J represents the number of pixels in the image; I, J and L correspond to the number of rows, columns and bands of the image respectively; A is the shared basis matrix of X and Y and H1 is the first multi-feature weight matrix; H2 is the second multi-feature weight matrix, C is the difference feature matrix of difference enhancement; ‖·‖ F is the Frobenius norm, ‖·‖ 2,1 For L 2,1 -norm; λ is the regularization parameter; (·) T Represents the transpose of the matrix; 1 R represents a column vector of all 1s with dimension R, 1 N represents a column vector of all 1s with dimension N; The S15 includes: S151: horizontally splicing the first temporal phase remote sensing image and the second temporal phase remote sensing image to obtain a spliced ​​three-dimensional image; S152: Expanding the spliced ​​three-dimensional image into a two-dimensional matrix; S153: performing singular value decomposition on the two-dimensional matrix to obtain a left singular value matrix, a singular value matrix and a right singular value matrix of the matrix; S154: Take the largest target number of singular values ​​and the corresponding left singular value matrix and right singular value matrix; S155: determining an initial shared basis matrix according to the product of the square root of the singular value and the corresponding column in the corresponding left singular value matrix, and determining the initial first multi-feature weight matrix and the initial second multi-feature weight matrix according to the product of the square root of the singular value and the corresponding column in the corresponding right singular value matrix; S156: Subtract the initial first multi-feature weight matrix from the initial second multi-feature weight matrix to obtain an initial differentially enhanced difference feature matrix.

2. The method for remote sensing image change detection based on dual non-negative matrix decomposition with joint sparse constraints according to claim 1, characterized in that: The S14 includes: S141: Generate L according to the objective function and the constraint conditions 2,1 -Augmented Lagrangian function of DNMF model; S142: Fix the L 2,1 The other variables in the augmented Lagrangian function of the -DNMF model are solved for the remaining single variables to obtain the update rules of the shared basis matrix, the update rules of the first multi-feature weight matrix, the update rules of the second multi-feature weight matrix and the update rules of the differential enhanced difference feature matrix.

3. The method for remote sensing image change detection based on dual non-negative matrix decomposition with joint sparse constraints according to claim 2, characterized in that: The L 2,1 The augmented Lagrangian function of the -DNMF model is: Where L represents the augmented Lagrangian function, A represents the shared basis matrix, H1 represents the first multi-feature weight matrix, H2 represents the second multi-feature weight matrix, M1, M2, M3 and M4 are Lagrangian multipliers; μ is the penalty parameter and μ ≥ 0; X is the first phase observation matrix and Y is the second phase observation matrix and N = I × J represents the number of pixels in the image; I, J and L correspond to the number of rows, columns and bands of the image respectively; C is the difference feature matrix of differential enhancement; ‖·‖ F is the Frobenius norm, ‖·‖ 2,1 For L 2,1 -norm; λ is the regularization parameter; tr(·) represents the trace of the matrix.

4. The method for remote sensing image change detection based on dual non-negative matrix decomposition with joint sparse constraints according to claim 2 or 3, characterized in that: The update rules of the shared basis matrix include: Among them, A i+1 represents the shared basis matrix of the i+1th iteration, A i represents the shared basis matrix of the i-th iteration, * represents the multiplication symbol, X is the first phase observation matrix and Y is the second phase observation matrix and N=I×J represents the number of pixels in the image; I, J and L correspond to the number of rows, columns and bands of the image respectively; H1 represents the first multi-feature weight matrix, and H2 represents the second multi-feature weight matrix; The updating rules of the first multi-feature weight matrix H1 include: in, represents the first multi-feature weight matrix of the i+1th iteration, represents the first multi-feature weight matrix of the i-th iteration, represents the augmented matrix of the first phase observation matrix X and represents the augmented matrix of the shared basis matrix A and C is the difference feature matrix of differential enhancement; μ is the penalty parameter, H2 represents the second multi-feature weight matrix, and δ is the parameter used to control the influence of normalization constraints; 1 R represents a column vector of all 1s with dimension R, 1 N represents a column vector of all 1s with dimension N; The update rules of the second multi-feature weight matrix H2 include: in, represents the augmented matrix of the second observation matrix and 5. The method for remote sensing image change detection based on dual non-negative matrix decomposition with joint sparse constraints according to claim 2 or 3, characterized in that: The S142 further includes: Using the L 2,1 -The augmented Lagrangian function of the DNMF model solves the objective function and the constraint conditions to obtain the update rule of the difference feature matrix of the difference enhancement, the update rule of the Lagrangian multiplier and the update rule of the penalty parameter; The updating rules of the difference feature matrix of the difference enhancement include: Among them, C i+1 represents the difference feature matrix of the difference enhancement of the i+1th iteration, Ω represents the shrinkage operator, λ is the regularization parameter; μ is the penalty parameter; H1 represents the first multi-feature weight matrix, H2 represents the second multi-feature weight matrix; M4 represents the Lagrange multiplier; The update rule of Lagrange multiplier M4 is: M4=M4+μ(C-H1+H2) The update rule of penalty parameter μ is: μ=min(ρμ,maxμ) Here, ρ is the learning rate used for parameter updating.

6. The method for remote sensing image change detection based on dual non-negative matrix decomposition with joint sparse constraints according to claim 1, characterized in that: In the step S2, clustering the difference feature matrix of the difference enhancement is performed using a multi-layer k-means clustering method; The multi-layer k-means clustering algorithm includes the following sub-steps: S21: Select 3 pixels as the initial cluster centers of the change class, the unchanged class, and the uncertain class; S22: Calculate the Euclidean distance between the pixel point and each of the initial cluster centers, and assign the pixel point to the cluster center closest to it; S23: Repeat step S22 until all pixels are allocated; S24: updating the position of one of the cluster centers, and moving the cluster center to the position of the data mean belonging to the cluster; S25: Repeat step S24 until the positions of all the cluster centers are updated; S26: Repeat steps S22-S25 until the position of the cluster center no longer changes or reaches a predetermined number of iterations. S27: Repeat steps S21-S26 for the pixels in the uncertain class until the uncertain class cannot be divided any further, so as to obtain a binary change result.

Images