Change detection method for multi-temporal hyperspectral images based on neural network hierarchical random walk

Abstract

The present invention discloses a method for detecting changes in multi-temporal hyperspectral images based on convolutional neural networks and hierarchical random walks, and belongs to the field of remote sensing image processing. The algorithm architecture of feature extraction-prior recognition-result optimization is adopted, which mainly includes three branches, namely, a 2DCNN branch for feature extraction of multi-temporal hyperspectral images, an affinity branch for characterizing the local similarity of multi-temporal hyperspectral images, and a hierarchical random walk link for optimizing the final change detection results. Experimental results show that the overall accuracy of the present invention on the three data sets of Irrigated Agricultural Area, Wetland Agricultural Area and River Area reached 0.950, 0.947 and 0.962 respectively, and the Kappa coefficients were 0.852, 0.878 and 0.755 respectively, which effectively improved the accuracy of change detection.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing image processing, and in particular to a multi-temporal hyperspectral image change detection method based on convolutional neural networks and hierarchical random walks, which can effectively improve the classification consistency between adjacent pixels, alleviate the problems of weakened boundary information and fragmented spatial information, and achieve improved measurement accuracy under a smaller training set. Background Art

[0002] The surface environment is constantly changing due to natural processes and human factors, and these changes are closely related to human living environment and safety. Therefore, rapid and accurate assessment and monitoring of changes in surface use and cover is not only crucial for the utilization and management of land resources, but also provides an important guarantee for better understanding the relationship and interaction between man and nature and achieving sustainable development. Since the nanoscale spectral resolution of hyperspectral remote sensing image data can provide fine spectral features, it brings great opportunities for the recognition and detection of fine changes in the surface using multi-temporal hyperspectral remote sensing images. At the same time, the rapidly expanding amount of data and excessive redundant features also bring challenges to the research of multi-temporal hyperspectral remote sensing image change detection algorithms.

[0003] Based on the rapid development of computer technology and remote sensing technology, researchers have proposed many methods for change detection in multi-temporal remote sensing images. The most common method is Change Vector Analysis (CVA), which is a typical unsupervised method that generates the amplitude and direction of changes by subtracting spectral vectors. Later, some scholars proposed some improved CVA methods. Adar et al. proposed a spectral overlap threshold method in 2012 to obtain a threshold that can distinguish between "changed" and "unchanged" areas. Although it can effectively suppress noise, its complexity is extremely high. Yuan et al. proposed a new semi-supervised distance metric learning method in 2015, which uses spectral features to detect change areas, but lacks adaptability to specific problems, which may have a great impact on the results of change detection. Zhang et al. constructed a mixed vector using change vectors and spectral angles in 2017, and then used an adaptive fusion strategy to generate the final change detection result map, but ignored the consideration of spatial information and texture information of hyperspectral images, which may have a great impact on the results of change detection.

[0004] In recent years, the computing power and data acquisition capabilities of computing devices have increased rapidly. A substantial increase in computing power can alleviate the inefficiency of training, and a substantial increase in training data can reduce the risk of "overfitting". As a result, complex models represented by deep learning technologies such as convolutional neural networks (CNNs) have been gradually applied to the task of change detection in hyperspectral remote sensing images, and have achieved results that are superior to traditional methods. In 2017, Gong et al. proposed a discriminative feature learning method for change detection based on restricted Boltzmann machines (RBMs), which converts dual-phase remote sensing images into feature space and then compares paired features to generate the final change detection map. However, the network does not have the powerful feature extraction and learning capabilities of CNN, which weakens the representation of features. In 2019, Wang et al. proposed a general end-to-end two-dimensional convolutional neural network framework (GETNET), which aims to learn features from higher-level multi-source data. However, the inherent noise of the pseudo training set will lead to a decrease in the performance of the algorithm. Kevin et al. (2019) used a pre-trained CNN model for semantic segmentation to detect changed regions in an unsupervised manner. However, this method heavily relied on the ability of the pre-trained model to perform semantic segmentation.

[0005] In general, although the hyperspectral image change detection technology has developed to a certain extent in the past two decades, most of the existing methods have limitations such as poor result accuracy, high algorithm redundancy, and complex calculations. Hyperspectral image data itself has the characteristics of high dimensionality and high redundancy, which makes the change information of multi-temporal hyperspectral images more implicit, mixed, and difficult to identify. It cannot be processed in some change detection algorithms, which brings difficulties to the change detection task. In addition, the methods applied to multispectral image change detection may not be able to give a correct change representation after compressing high-dimensional data. Although feature extraction or band selection can be used to reduce the dimensionality of hyperspectral images, fine change information may be lost to some extent, thereby reducing the accuracy of change detection. For the hyperspectral image change detection task, the main difficulty is to effectively extract change information from the high-dimensional feature space. Summary of the invention

[0006] The present invention aims to solve the above-mentioned technical problems existing in the prior art and to provide a multi-temporal hyperspectral image change detection method based on convolutional neural network and hierarchical random walk, which can effectively improve the classification consistency between adjacent pixels, alleviate the problems of weakening boundary information and fragmentation of spatial information, and achieve improved measurement accuracy under a smaller training set.

[0007] The technical solution of the present invention is: a multi-temporal hyperspectral image change detection method based on convolutional neural network and hierarchical random walk, which is characterized by the following steps:

[0008] Step 1. Establish and initialize the convolutional neural network N for hyperspectral image change detection CD , the N CD Contains 1 subnetwork N for feature extraction feature and 1 subnetwork N for classification cls ;

[0009] Step 2. Input the hyperspectral image I1∈R of the first phase M×N×D , the hyperspectral image of the second phase I2∈R M ×N×D , the manually labeled pixel point coordinate set and label set, for N CD Training is performed, where M and N are the spatial scales of the hyperspectral image, and D is the spectral dimension;

[0010] Step 3. Input all pixels of the multi-phase hyperspectral difference image DI for data preprocessing, and use the trained convolutional neural network N CD Complete pixel classification;

[0011] Step 4. Model the multi-temporal hyperspectral difference image as a weighted undirected connected graph, calculate the Euclidean distance between each pixel and its adjacent pixels, and obtain the weight matrix and affinity matrix;

[0012] Step 5. Calculate the transition probability and arrival probability, optimize the random walk and obtain the final labeling result.

[0013] The step 1 is performed as follows:

[0014] Step 1.1 Create and initialize subnetwork N feature , contains 6 groups of convolutional layers, namely Conv2_0, Conv2_1, Conv2_2, Conv2_3, Conv2_4, Conv2_5;

[0015] Conv2_0 includes one convolution operation layer, wherein the convolution layer contains 128 convolution kernels of size 1×1, and each convolution kernel performs convolution operation with a step size of 1 pixel;

[0016] The Conv2_1 includes 1 convolution operation, 1 BatchNorm normalization operation and 1 activation operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation;

[0017] The Conv2_2 includes 1 convolution operation, 1 BatchNorm normalization operation and 1 activation operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation;

[0018] The Conv2_3 includes 1 convolution operation, 1 activation operation and 1 maximum pooling operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation, and the pooling layer performs the maximum pooling operation with a pooling kernel of size 2×2;

[0019] The Conv2_4 includes 1 convolution operation, 1 BatchNorm normalization operation and 1 activation operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation;

[0020] The Conv2_5 includes 1 convolution operation, 1 BatchNorm normalization operation, 1 activation operation, 1 maximum pooling operation and 1 Flatten operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation, and the pooling layer performs the maximum pooling operation with a pooling kernel of size 2×2;

[0021] Step 1.2 Create and initialize subnetwork N cls , contains 1 set of fully connected layers, Dense1;

[0022] The Dense1 has num classification units and adopts Softmax as the activation function, where num represents the total number of object categories to be classified.

[0023] The step 2 is performed as follows:

[0024] Step 2.1 Obtain multi-temporal hyperspectral difference image DI according to the Log-Ratio method of formula (1);

[0025] DI=|log(I2-I1)| (1)

[0026] Step 2.2: Extract all pixel points with labels from the multi-temporal hyperspectral difference image DI based on the manually labeled pixel point coordinate set. in, Represents X D The i-thTrain Pixel points, M Train Indicates the total number of pixel points with labels;

[0027] Step 2.3 According to the definition of formula (2), we get in, represents the standardized set of hyperspectral image pixel points with labels. express The i Train Pixel points;

[0028]

[0029] Step 2.4 Each pixel point is taken as the center to divide DI into a series of 9×9 hyperspectral pixel block sets X H1 ;

[0030] Step 2.5 X H1 Each pixel block in is flipped upside down to obtain the hyperspectral pixel block set X H2 ;

[0031] Step 2.6: X H1 Add Gaussian noise with a variance of 0.01 to each pixel block in the image to obtain the hyperspectral pixel block set X H3 ;

[0032] Step 2.7 X H1 Each pixel block in is randomly rotated clockwise by n×90 degrees with its center point as the rotation center to obtain the hyperspectral pixel block set X H4 , where n represents a value randomly selected from the set {1,2,3};

[0033] Step 2.8: Will As the training set of the convolutional neural network, let the number of iterations iter←1 and execute step 2.9;

[0034] Step 2.9: Use network N CD Extract the features of the training set and classify the test results;

[0035] Step 2.9.1 Using subnetwork N feature Training set for hyperspectral images Perform feature extraction to obtain the feature F of the hyperspectral image;

[0036] Step 2.9.2 Use subnetwork N cls Classify the features and calculate the classification prediction result TR pred ;

[0037] Step 2.10: According to the definitions of formula (3) and formula (4), the weighted cross entropy is used as the loss function;

[0038]

[0039]

[0040] in, Indicates the jth cls The weight of the class, Indicates that the pixel belongs to the jth cls The probability of a ground-like object, Indicates the jth ground-truth training sample cls The number of features;

[0041] Step 2.11: If all pixel blocks in the training set have been processed, go to step 2.12; otherwise, take a group of unprocessed pixel blocks from the training set and return to step 2.9;

[0042] Step 2.12 Let iter←iter+1. If the number of iterations iter>Total_iter, the trained convolutional neural network N is obtained. CD , go to step 3, otherwise, use the reverse error propagation algorithm based on stochastic gradient descent and the prediction loss L ω-C Update N ahd , go to step 2.9 to reprocess all pixel blocks in the training set, and Total_iter represents the preset number of iterations.

[0043] The step 3 is performed as follows:

[0044] Step 3.1 Extract all pixel points in DI to form a set in, Indicates T D The i Test Pixels, M Test Represents the total number of all pixels;

[0045] Step 3.2 According to the definition of formula (6), D After standardization, we get in, represents the standardized set of hyperspectral image pixel points with labels. express The i Test Pixel points;

[0046]

[0047] Step 3.3 The DI is divided into a series of 9×9 hyperspectral pixel blocks with each pixel point as the center to form a hyperspectral image test set

[0048] Step 3.4 Use network N CD Extract the features of the test set and classify the test results;

[0049] Step 3.4.1 Using subnetwork N feature Test set for hyperspectral images Perform feature extraction to obtain the feature T of the hyperspectral image;

[0050] Step 3.4.2 Use subnetwork N cls Classify the features and calculate the classification prediction result TR pred ;

[0051] Step 3.5: Use subnetwork N cls Classify the feature T to calculate the classification prediction result TE pred .

[0052] The step 4 is performed as follows:

[0053] Step 4.1 Model the hyperspectral difference image DI as a weighted undirected connected graph G = {V, E}, where V = {v i |i=1,…,M×N} is the vertex, i.e., the set of hyperspectral image pixel points, v i represents the i-th pixel point in V; E = {e ij |i=1,…,M×N;j=1,…,M×N;i<j} is the set of edges connecting each pixel point, e ij Represents the element at position (i, j) in E; define Eu{eu ij |i=1,…,M×N;j=1,…,M×N} is a sparse matrix that stores the Euclidean distances between each pixel and all its neighboring pixels. Its size is a sparse matrix of (M×N)×(M×N), as shown in formula (7), where eu ij Represents the element at position (i,j) in Eu:

[0054] eu ij =||I D (i)-I D (j)||2 (7)

[0055] Step 4.2 optimizes Eu to obtain the weight matrix W of the hyperspectral difference image DI = {ω ij |i=1,…,M×N;j=1,…,M×N}, where ω ijRepresents the element at position (i, j) in W. Its calculation method is shown in formula (8), where σ is the control parameter and ε is a very small constant:

[0056]

[0057] Step 4.3: A=D by matrix operation -1 The affinity matrix A of W to the pixel points of the hyperspectral difference image is A = {a ij |i=1,…,M×N;j=1,…,M×N;i<j} is used for prediction, and its calculation method is shown in formula (9), where D is a diagonal matrix, d i =∑ i≠j ω ij represents the element at position (i,i) in D, a ij Represents the element at position (i,j) in A:

[0058]

[0059] The step 5 is performed as follows:

[0060] Step 5: Calculate the transition probability and arrival probability, optimize the random walk and obtain the final marking result;

[0061] Step 5.1: The classification prediction result TE calculated in step 3.5 pred It is expressed as P = {p i |i=1,…,M×N} is used as a global prior, where p i is the i-th element of P, and defines the set of seeds of the initial training point location information S = {s i |i S =1,…,M×N} is used as local guidance, where s i is the i-th element of S;

[0062] Step 5.2 constructs the Laplacian matrix L of the non-regularized graph, and the calculation method is shown in formula (10);

[0063] L=DW (9)

[0064] Step 5.3 Use prior information to construct a graph G with prior information e ; Assume that the random walk starts from each unlabeled pixel and calculate the probability r k Make it reach the marked pixel first; for the standard random walk, the probability r k By calculating the energy function The minimum value of is obtained, and the calculation method is shown in formula (11), where the energy function The analytical solution to is calculated by solving the linear system of equations:

[0065]

[0066] Step 5.4 Construct the prior probability diagonal matrix The calculation method is shown in formula (12), where diag(P) = diag(p i |i=1,…,M×N) is a diagonal matrix with all elements in P as main diagonal elements;

[0067] P k =Adiag(P) (12)

[0068] Step 5.5 Construct non-spatial energy function The calculation method is shown in formula (13), where S k =diag(s i |i=1,…,M×N) is the seed matrix, C is the weight matrix controlling the seed trade-off, (r k -S k ) T C(r k -S k ) is the energy of random walk to reach the training seed, whose minimum value can aggregate nodes to seed nodes in their category; (r k -J) T P k (r k -J) and Respectively represent the arrival probability of the node in the prior distribution of the current category label and other category labels;

[0069]

[0070] Step 5.6: Combine the spatial energy function and the non-spatial energy function according to formula (14) to express it as E k (r k ), where λ is the weight number that controls the trade-off between spatial energy and non-spatial energy. Each prior distribution node is connected to all nodes in V;

[0071]

[0072] Step 5.7: The transition probability matrix Q on the set V∪Δ∪S∪P is expressed according to formula (15), where: θ is the weight parameter of the prior distribution, c is the weight parameter of the seed, and Δ represents the set of unlabeled pixels;

[0073]

[0074] Step 5.8 Based on the transition probability Q of the prior graph, the probability of a random walk from a node to a seed or a prior node is r ik The formula is shown in formula (16):

[0075]

[0076] Step 5.9 The calculation method of the final marking result R(i) is shown in formula (17);

[0077]

[0078] Compared with the prior art, the present invention improves the accuracy of change detection at the two levels of feature extraction and decision prediction, and includes three branches: a 2D CNN branch for multi-phase hyperspectral image feature extraction and classification prediction results; an affinity branch for characterizing the local similarity of multi-phase hyperspectral images and obtaining an affinity matrix; and a hierarchical random walk model for optimizing the final change detection results. The model architecture of feature extraction, prior recognition and result optimization proposed in the present invention can not only avoid the highly redundant information and highly complex data structure of the hyperspectral remote sensing image data itself, but also effectively improve the classification consistency between adjacent pixels, alleviate the problem of weakening boundary information and fragmenting spatial information due to the presence of a large receptive field and pooling layer in the convolution layer of the deep convolutional neural network, and thus improve the accuracy of the multi-phase hyperspectral image change detection algorithm under a simple network structure with small training samples. Therefore, the present invention has the characteristics of good feature extraction quality, high detection accuracy, and few required training samples. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 This is a diagram of the detection network structure of an embodiment of the present invention.

[0080] Figure 2 This is a comparison chart of the change detection results of the Irrigated Agricultural Area dataset using the present invention, the CVA method, the PCDA method, the ED method, the Image Regr method, the CNN method, and the DSFA method.

[0081] Figure 3 This is a comparison chart of the change detection results of the Wetland Agricultural Area data set between the present invention and the CVA method, PCDA method, ED method, Image Regr method, CNN method and DSFA method.

[0082] Figure 4 This is a comparison chart of the change detection results of the River Area data set between the present invention and the CVA method, PCDA method, ED method, Image Regr method, CNN method and DSFA method. DETAILED DESCRIPTION

[0083] The present invention provides a multi-temporal hyperspectral image change detection method based on convolutional neural network and hierarchical random walk, such as Figure 1 As shown, follow the steps below:

[0084] Step 1. Establish and initialize the convolutional neural network N for hyperspectral image change detection CD , the N CD Contains 1 subnetwork N for feature extraction feature and 1 subnetwork N for classification cls ;

[0085] Step 1.1 Create and initialize subnetwork N feature , contains 6 groups of convolutional layers, namely Conv2_0, Conv2_1, Conv2_2, Conv2_3, Conv2_4, Conv2_5;

[0086] Conv2_0 includes one convolution operation layer, wherein the convolution layer contains 128 convolution kernels of size 1×1, and each convolution kernel performs convolution operation with a step size of 1 pixel;

[0087] The Conv2_1 includes 1 convolution operation, 1 BatchNorm normalization operation and 1 activation operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation;

[0088] The Conv2_2 includes 1 convolution operation, 1 BatchNorm normalization operation and 1 activation operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation;

[0089] The Conv2_3 includes 1 convolution operation, 1 activation operation and 1 maximum pooling operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation, and the pooling layer performs the maximum pooling operation with a pooling kernel of size 2×2;

[0090] The Conv2_4 includes 1 convolution operation, 1 BatchNorm normalization operation and 1 activation operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation;

[0091] The Conv2_5 includes 1 convolution operation, 1 BatchNorm normalization operation, 1 activation operation, 1 maximum pooling operation and 1 Flatten operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation, and the pooling layer performs the maximum pooling operation with a pooling kernel of size 2×2;

[0092] Step 1.2 Create and initialize subnetwork N cls , contains 1 set of fully connected layers, Dense1;

[0093] The Dense1 has num classification units and uses Softmax as the activation function, where num represents the total number of object categories to be classified;

[0094] Step 2. Input the hyperspectral image I1∈R of the first phase M×N×D , the hyperspectral image of the second phase I2∈R M ×N×D , the manually labeled pixel point coordinate set and label set, for N CD Training is performed, where M and N are the spatial scales of the hyperspectral image, and D is the spectral dimension;

[0095] Step 2.1 Obtain multi-temporal hyperspectral difference image DI according to the Log-Ratio method of formula (1);

[0096] DI=|log(I2-I1)| (1)

[0097] Step 2.2: Extract all pixel points with labels from the multi-temporal hyperspectral difference image DI based on the manually labeled pixel point coordinate set. in, Represents X D The i-th Train Pixel points, M Train Indicates the total number of pixel points with labels;

[0098] Step 2.3 According to the definition of formula (2), we get in, represents the standardized set of hyperspectral image pixel points with labels. express The i Train Pixel points;

[0099]

[0100] Step 2.4 Each pixel point is taken as the center to divide DI into a series of 9×9 hyperspectral pixel block sets X H1 ;

[0101] Step 2.5 X H1 Each pixel block in is flipped upside down to obtain the hyperspectral pixel block set X H2 ;

[0102] Step 2.6: X H1 Add Gaussian noise with a variance of 0.01 to each pixel block in the image to obtain the hyperspectral pixel block set X H3 ;

[0103] Step 2.7 X H1 Each pixel block in is randomly rotated clockwise by n×90 degrees with its center point as the rotation center to obtain the hyperspectral pixel block set X H4 , where n represents a value randomly selected from the set {1,2,3};

[0104] Step 2.8: Will As the training set of the convolutional neural network, let the number of iterations iter←1 and execute step 2.9;

[0105] Step 2.9: Use network N CD Extract the features of the training set and classify the test results;

[0106] Step 2.9.1 Using subnetwork N feature Training set for hyperspectral images Perform feature extraction to obtain the feature F of the hyperspectral image;

[0107] Step 2.9.2 Use subnetwork N cls Classify the features and calculate the classification prediction result TR pred ;

[0108] Step 2.10: According to the definitions of formula (3) and formula (4), the weighted cross entropy is used as the loss function;

[0109]

[0110]

[0111] in, Indicates the jth cls The weight of the class, Indicates that the pixel belongs to the jth cls The probability of a ground-like object, Indicates the jth ground-truth training sample cls The number of features;

[0112] Step 2.11: If all pixel blocks in the training set have been processed, go to step 2.12; otherwise, take a group of unprocessed pixel blocks from the training set and return to step 2.9;

[0113] Step 2.12 Let iter←iter+1. If the number of iterations iter>Total_iter, the trained convolutional neural network N is obtained. CD , go to step 3, otherwise, use the reverse error propagation algorithm based on stochastic gradient descent and the prediction loss L ω-C Update N ahd , go to step 2.9 to reprocess all pixel blocks in the training set, and the Total_iter represents the preset number of iterations;

[0114] Step 3. Input all pixels of the multi-phase hyperspectral difference image DI for data preprocessing, and use the trained convolutional neural network N CD Complete pixel classification;

[0115] Step 3.1 Extract all pixel points in DI to form a set in, Indicates T D The i Test Pixels, M Test Represents the total number of all pixels;

[0116] Step 3.2 According to the definition of formula (6), D After standardization, we get in, represents the standardized set of hyperspectral image pixel points with labels. express The i Test Pixel points;

[0117]

[0118] Step 3.3 The DI is divided into a series of 9×9 hyperspectral pixel blocks with each pixel point as the center to form a hyperspectral image test set

[0119] Step 3.4 Use network N CD Extract the features of the test set and classify the test results;

[0120] Step 3.4.1 Using subnetwork N feature Test set for hyperspectral images Perform feature extraction to obtain the feature T of the hyperspectral image;

[0121] Step 3.4.2 Use subnetwork N cls Classify the features and calculate the classification prediction result TR pred ;

[0122] Step 3.5: Use subnetwork N cls Classify the feature T to calculate the classification prediction result TE pred ;

[0123] Step 4. Model the multi-temporal hyperspectral difference image as a weighted undirected connected graph, calculate the Euclidean distance between each pixel and its adjacent pixels, and obtain the weight matrix and affinity matrix (affinity branch);

[0124] Step 4.1 Model the hyperspectral difference image DI as a weighted undirected connected graph G = {V, E}, where V = {v i |i=1,…,M×N} is the vertex, i.e., the set of hyperspectral image pixel points, v i represents the i-th pixel point in V; E = {e ij |i=1,…,M×N;j=1,…,M×N;i<j} is the set of edges connecting each pixel point, e ij Represents the element at position (i, j) in E; define Eu{eu ij |i=1,…,M×N;j=1,…,M×N} is a sparse matrix that stores the Euclidean distances between each pixel and all its neighboring pixels. Its size is a sparse matrix of (M×N)×(M×N), as shown in formula (7), where eu ij Represents the element at position (i,j) in Eu:

[0125] eu ij =||I D (i)-I D (j)||2 (7)

[0126] Step 4.2 optimizes Eu to obtain the weight matrix W of the hyperspectral difference image DI = {ω ij |i=1,…,M×N;j=1,…,M×N}, where ω ij Represents the element at position (i, j) in W. Its calculation method is shown in formula (8), where σ is the control parameter and ε is a very small constant:

[0127]

[0128] Step 4.3: A=D by matrix operation -1 The affinity matrix A of W to the pixel points of the hyperspectral difference image is A = {a ij|i=1,…,M×N;j=1,…,M×N;i<j} is used for prediction, and its calculation method is shown in formula (9), where D is a diagonal matrix, d i =∑ i≠j ω ij represents the element at position (i,i) in D, a ij Represents the element at position (i,j) in A:

[0129]

[0130] Step 5. Calculate the transition probability and arrival probability, optimize the random walk and obtain the final marking result;

[0131] Step 5.1: The classification prediction result TE calculated in step 3.5 pred It is expressed as P = {p i |i=1,…,M×N} is used as a global prior, where p i is the i-th element of P, and defines the set of seeds of the initial training point location information S = {s i |i S =1,…,M×N} is used as local guidance, where s i is the i-th element of S;

[0132] Step 5.2 constructs the Laplacian matrix L of the non-regularized graph, and the calculation method is shown in formula (10);

[0133] L=DW (9)

[0134] Step 5.3 Use prior information to construct a graph G with prior information e ; Assume that the random walk starts from each unlabeled pixel and calculate the probability r k Make it reach the marked pixel first; for the standard random walk, the probability r k By calculating the energy function The minimum value of is obtained, and the calculation method is shown in formula (11), where the energy function The analytical solution to is calculated by solving the linear system of equations:

[0135]

[0136] Step 5.4 Construct the prior probability diagonal matrix The calculation method is shown in formula (12), where diag(P) = diag(p i |i=1,…,M×N) is a diagonal matrix with all elements in P as main diagonal elements;

[0137] P k =Adiag(P) (12)

[0138] Step 5.5 Construct non-spatial energy function The calculation method is shown in formula (13), where S k =diag(s i |i=1,…,M×N) is the seed matrix, C is the weight matrix controlling the seed trade-off, (r k -S k ) T C(r k -S k ) is the energy of random walk to reach the training seed, whose minimum value can aggregate nodes to seed nodes in their category; (r k -J) T P k (r k -J) and Respectively represent the arrival probability of the node in the prior distribution of the current category label and other category labels;

[0139]

[0140] Step 5.6: Combine the spatial energy function and the non-spatial energy function according to formula (14) to express it as E k (r k ), where λ is the weight number that controls the trade-off between spatial energy and non-spatial energy. Each prior distribution node is connected to all nodes in V;

[0141]

[0142] Step 5.7: The transition probability matrix Q on the set V∪Δ∪S∪P is expressed according to formula (15), where: θ is the weight parameter of the prior distribution, c is the weight parameter of the seed, and Δ represents the set of unlabeled pixels;

[0143]

[0144] Step 5.8 Based on the transition probability Q of the prior graph, the probability of a random walk from a node to a seed or a prior node is r i k The formula is shown in formula (16):

[0145]

[0146] Step 5.9 The calculation method of the final marking result R(i) is shown in formula (17);

[0147]

[0148] To verify the effectiveness of the present invention, experiments were carried out using the publicly available Irrigated Agricultural Area dataset, Wetland Agricultural Area dataset and River Area dataset as examples. The overall accuracy (OA) and Kappa coefficient (KC) were used as objective indicators to evaluate the transformation detection results. The evaluation results of the present invention were compared with those of the CVA method, PCDA method, ED method, Image Regr method, CNN method and DSFA method.

[0149] The present invention has achieved the result closest to the Ground Truth. In the change detection task, KC is considered to be the most convincing objective evaluation index. The KC value of the present invention is 1.8%, 2.6% and 2.8% higher than the second best algorithm on the three data sets of Irrigated Agricultural Area, Wetland Agricultural Area and Wetland Agricultural Area, respectively. Therefore, it can be shown that the present invention is superior to other algorithms in terms of accuracy and has better performance.

[0150] from Figure 2 , Figure 3 and Figure 4It can be seen that since the changes in the Irrigated Agricultural Area dataset itself are mostly the expansion process of farmland, its spatiotemporal evolution process is relatively simple. Therefore, CVA, PCDA, ED, Image Regr and the method of the present invention have achieved good results. However, the result graph of the CNN method has more noise points, and although the result graph of the DSFA method has fewer noise points than the CNN method, there is a serious misdetection phenomenon in the upper left corner of the image, that is, it detects pixels that have not changed as changed pixels. The present invention combines convolutional neural network with hierarchical random walk. Compared with the results of CNN algorithm, the introduction of hierarchical random walk effectively solves the noise problem, and still achieves better results even in the case of a smaller training sample set. Although the Wetland Agricultural Area dataset itself changes regularly, there is a serious missed detection phenomenon in CVA, ED, and Image Regr methods. These three methods cause many false negative pixels in the result map, so that many changed pixels are not detected. Even in PCDA and DSFA methods, the visual effect on this dataset is worse, especially the DSFA method, which does not distinguish between changed and unchanged pixels well, and the false detection and missed detection phenomena are relatively serious. The result map of CNN method and the present invention on this dataset is better, and it can be seen that the present invention is more detailed in edge processing than the CNN algorithm. River The changes in the Area dataset itself are very complex, and the change results in all algorithms are two-sided. No good results were obtained under the CVA, PCDA and ED algorithms, while the results obtained by the remaining four algorithms are relatively considerable. Among them, the result graphs of the ImageRegr and DSFA methods have more false positive pixels, that is, they tend to detect unchanged pixels as changed pixels, while the present invention and the CNN algorithm have some false negative pixels. Although there are some missed detection phenomena, they retain more detailed information than other algorithms and can capture more subtle changes.

[0151] Comprehensive Table 1, Table 2, Table 3, Figure 2 , Figure 3 , Figure 4The comparison results show that the present invention effectively improves the classification consistency between adjacent pixels and alleviates the problems of weakened boundary information and fragmented spatial information by constructing a multi-phase hyperspectral image change detection model based on a deep convolutional neural network and a hierarchical random walk model, thereby improving the accuracy of the change detection algorithm under a simple network structure with a small training sample. The overall accuracy of the present invention on the three data sets of Irrigated Agricultural Area, Wetland Agricultural Area and River Area reached 0.950, 0.947 and 0.962 respectively, and the Kappa coefficients were 0.852, 0.878 and 0.755 respectively, which effectively improved the accuracy of change detection.

[0152] Table 1 Objective evaluation indicators of the Irrigated Agricultural Area dataset under different algorithms

[0153]

[0154] Table 2 Objective evaluation indicators of the Wetland Agricultural Area dataset under different algorithms

[0155]

[0156] Table 3 Objective evaluation indicators of the River Area dataset under different algorithms

[0157]

Claims

1. A multi-temporal hyperspectral image change detection method based on convolutional neural network and hierarchical random walk, characterized in that Proceed as follows: Step 1. Establish and initialize the convolutional neural network N for hyperspectral image change detection CD , the N CD Contains 1 subnetwork N for feature extraction feature and 1 subnetwork N for classification cls ; Step 2. Input the hyperspectral image I1∈R of the first phase M×N×D , the second phase of the hyperspectral image I2∈R M×N×D , manually annotated pixel point coordinate set and label set, for N CD Training is performed, where M and N are the spatial scales of the hyperspectral image, and D is the spectral dimension; Step 3. Input all pixels of the multi-phase hyperspectral difference image DI for data preprocessing, and use the trained convolutional neural network N CD Complete pixel classification; Step 4. Model the multi-temporal hyperspectral difference image as a weighted undirected connected graph, calculate the Euclidean distance between each pixel and its adjacent pixels, and obtain the weight matrix and affinity matrix; Step 5. Calculate the transition probability and arrival probability, optimize the random walk and obtain the final labeling result.

2. The multi-temporal hyperspectral image change detection method based on convolutional neural network and hierarchical random walk according to claim 1 is characterized in that The step 1 is carried out as follows: Step 1.1 Create and initialize subnetwork N feature , contains 6 groups of convolutional layers, namely Conv2_0, Conv2_1, Conv2_2, Conv2_3, Conv2_4, Conv2_5; Conv2_0 includes one convolution operation layer, wherein the convolution layer contains 128 convolution kernels of size 1×1, and each convolution kernel performs convolution operation with a step size of 1 pixel; The Conv2_1 includes 1 convolution operation, 1 BatchNorm normalization operation and 1 activation operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation; The Conv2_2 includes 1 convolution operation, 1 BatchNorm normalization operation and 1 activation operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation; The Conv2_3 includes 1 convolution operation, 1 activation operation and 1 maximum pooling operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation, and the pooling layer performs the maximum pooling operation with a pooling kernel of size 2×2; The Conv2_4 includes 1 convolution operation, 1 BatchNorm normalization operation and 1 activation operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation; The Conv2_5 includes 1 convolution operation, 1 BatchNorm normalization operation, 1 activation operation, 1 maximum pooling operation and 1 Flatten operation, wherein the convolution layer contains 256 convolution kernels of size 3×3, each convolution kernel performs convolution operation with a step size of 1 pixel, and uses the nonlinear activation function ReLU as the activation function for operation, and the pooling layer performs the maximum pooling operation with a pooling kernel of size 2×2; Step 1.2 Create and initialize subnetwork N cls , contains 1 set of fully connected layers, Dense1; The Dense1 has num classification units and adopts Softmax as the activation function, where num represents the total number of object categories to be classified.

3. The multi-temporal hyperspectral image change detection method based on convolutional neural network and hierarchical random walk according to claim 2 is characterized in that The step 2 is performed as follows: Step 2.1 Obtain multi-temporal hyperspectral difference image DI according to the Log-Ratio method of formula (1); DI=|log(I2-I1)| (1) Step 2.2: Extract all pixel points with labels from the multi-temporal hyperspectral difference image DI based on the manually labeled pixel point coordinate set. in, Represents X D The i-th Train Pixel points, M Train Represents the total number of pixel points with labels; Step 2.3 According to the definition of formula (2), we get in, represents the standardized set of hyperspectral image pixel points with labels. express The i Train Pixel points; Step 2.4 Each pixel point is taken as the center to divide DI into a series of 9×9 hyperspectral pixel block sets X H1 ; Step 2.5 X H1 Each pixel block in is flipped upside down to obtain the hyperspectral pixel block set X H2 ; Step 2.6: X H1 Add Gaussian noise with a variance of 0.01 to each pixel block in the image to obtain the hyperspectral pixel block set X H3 ; Step 2.7 X H1 Each pixel block in is randomly rotated clockwise by n×90 degrees with its center point as the rotation center to obtain the hyperspectral pixel block set X H4 , where n represents a value randomly selected from the set {1,2,3}; Step 2.8: Will As the training set of the convolutional neural network, let the number of iterations iter←1 and execute step 2.9; Step 2.9: Use network N CD Extract the features of the training set and classify the test results; Step 2.9.1 Using subnetwork N feature Training set for hyperspectral images Perform feature extraction to obtain the feature F of the hyperspectral image; Step 2.9.2 Use subnetwork N cls Classify the features and calculate the classification prediction result TR pred ; Step 2.10: According to the definitions of formula (3) and formula (4), the weighted cross entropy is used as the loss function; in, Indicates the jth cls The weight of the class, Indicates that the pixel belongs to the jth cls The probability of a ground-like object, Indicates the jth ground-truth training sample cls The number of features; Step 2.11: If all pixel blocks in the training set have been processed, go to step 2.12; otherwise, take a group of unprocessed pixel blocks from the training set and return to step 2.9; Step 2.12 Let iter←iter+1. If the number of iterations iter>Total_iter, the trained convolutional neural network N is obtained. CD , go to step 3, otherwise, use the reverse error propagation algorithm based on stochastic gradient descent and the prediction loss L ω-C Update N ahd , go to step 2.9 to reprocess all pixel blocks in the training set, and Total_iter represents the preset number of iterations.

4. The multi-temporal hyperspectral image change detection method based on convolutional neural network and hierarchical random walk according to claim 3 is characterized in that The step 3 is performed as follows: Step 3.1 Extract all pixel points in DI to form a set in, Indicates T D The i Test Pixels, M Test Represents the total number of all pixels; Step 3.2 According to the definition of formula (6), D After standardization, we get in, represents the standardized set of hyperspectral image pixel points with labels. express The i Test Pixel points; Step 3.3 The DI is divided into a series of 9×9 hyperspectral pixel blocks with each pixel point as the center to form a hyperspectral image test set Step 3.4 Use network N CD Extract the features of the test set and classify the test results; Step 3.4.1 Using subnetwork N feature Test set for hyperspectral images Perform feature extraction to obtain the feature T of the hyperspectral image; Step 3.4.2 Use subnetwork N cls Classify the features and calculate the classification prediction result TR pred ; Step 3.5: Use subnetwork N cls Classify the feature T to calculate the classification prediction result TE pred .

5. The multi-temporal hyperspectral image change detection method based on convolutional neural network and hierarchical random walk according to claim 4 is characterized in that The step 4 is performed as follows: Step 4.1 Model the hyperspectral difference image DI as a weighted undirected connected graph G = {V, E}, where V = {v i |i=1,…,M×N} is the vertex, i.e., the set of hyperspectral image pixel points, v i represents the i-th pixel point in V; E = {e ij |i=1,…,M×N;j=1,…,M×N;i<j} is the set of edges connecting each pixel point, e ij Represents the element at position (i, j) in E; define Eu{eu ij |i=1,…,M×N;j=1,…,M×N} is a sparse matrix that stores the Euclidean distances between each pixel and all its neighboring pixels. Its size is a sparse matrix of (M×N)×(M×N), as shown in formula (7), where eu ij Represents the element at position (i,j) in Eu: I ij =||I D (i)-I D (j)||2 (7) Step 4.2 optimizes Eu to obtain the weight matrix W of the hyperspectral difference image DI = {ω ij |i=1,…,M×N;j=1,…,M×N}, where ω ij represents the element at position (i, j) in W. Its calculation method is shown in formula (8), where σ is the control parameter and ε is a very small constant: Step 4.3: A=D by matrix operation -1 The affinity matrix A of W to the pixel points of the hyperspectral difference image is A = {a ij |i=1,…,M×N;j=1,…,M×N;i<j} is used for prediction, and its calculation method is shown in formula (9), where D is a diagonal matrix, d i =Σ i≠j ω ij represents the element at position (i,i) in D, a ij Represents the element at position (i,j) in A:

6. The multi-temporal hyperspectral image change detection method based on convolutional neural network and hierarchical random walk according to claim 5 is characterized in that The step 5 is performed as follows: Step 5: Calculate the transition probability and arrival probability, optimize the random walk and obtain the final marking result; Step 5.1: The classification prediction result TE calculated in step 3.5 pred It is expressed as P = {p i |i=1,…,M×N} is used as a global prior, where p i is the i-th element of P, and defines the set of seeds of the initial training point location information S = {s i |i S =1,…,M×N} is used as local guidance, where s i is the i-th element of S; Step 5.2 constructs the Laplacian matrix L of the non-regularized graph, and the calculation method is shown in formula (10); L=DW (9) Step 5.3 Use prior information to construct a graph G with prior information e ; Assume that the random walk starts from each unlabeled pixel and calculate the probability r k Make it reach the marked pixel first; for the standard random walk, the probability r k By calculating the energy function The minimum value of is obtained, and the calculation method is shown in formula (11), where the energy function The analytical solution to is calculated by solving the linear system of equations: Step 5.4 Construct the prior probability diagonal matrix The calculation method is shown in formula (12), where diag(P) = diag(p i |i=1,…,M×N) is a diagonal matrix with all elements in P as main diagonal elements; P k =Adiag(P) (12) Step 5.5 Construct non-spatial energy function The calculation method is shown in formula (13), where S k =diag(s i |i=1,…,M×N) is the seed matrix, C is the weight matrix controlling the seed trade-off, (r k -S k ) T C(r k -S k ) is the energy of random walk to reach the training seed, whose minimum value can aggregate nodes to seed nodes in their category; (r k -J) T P k (r k -J) and Respectively represent the arrival probability of the node in the prior distribution of the current category label and other category labels; Step 5.6: Combine the spatial energy function and the non-spatial energy function according to formula (14) to express it as E k (r k ), where λ is the weight number that controls the trade-off between spatial energy and non-spatial energy. Each prior distribution node is connected to all nodes in V; Step 5.7: The transition probability matrix Q on the set V∪Δ∪S∪P is expressed according to formula (15), where: θ is the weight parameter of the prior distribution, c is the weight parameter of the seed, and Δ represents the set of unlabeled pixels; Step 5.8 Based on the transition probability Q of the prior graph, the probability of a random walk from a node to a seed or a prior node is r i k The formula is shown in formula (16): Step 5.9 The calculation method of the final marking result R(i) is shown in formula (17);

Images