Hyperspectral Anomaly Detection Method Combining Robust Dictionary and Dual Collaborative Constraint Regularization Terms

Abstract

The present invention proposes a hyperspectral anomaly detection method that integrates a robust dictionary and a dual collaborative constraint regularization term to solve the problems of poor robustness during dictionary construction and insufficient information utilization, resulting in poor detection results. It includes: 1) obtaining a density map using a density estimation model and a local outlier factor, and performing pixel-by-pixel multiplication and pixel-by-pixel addition operations on it to obtain a background suppression map and an anomaly enhancement map; 2) performing binarization processing on the background suppression map and the anomaly enhancement map to obtain their index maps, and constructing a potential anomaly dictionary and a background dictionary; 3) constructing a low-rank and sparse representation model based on the dual collaborative constraint regularization term and performing optimization to obtain an initial detection result; 4) performing a non-linear transformation on the initial detection result using the background suppression map and the anomaly enhancement map to obtain the final detection result. The present invention enhances the robustness of the dictionary and fully considers the global and local characteristics, thereby effectively improving the hyperspectral anomaly detection effect.

Description

Technical Field

[0001] The present invention belongs to the technical field of hyperspectral remote sensing image processing, and further relates to hyperspectral anomaly detection. Specifically, it is a hyperspectral anomaly detection method that combines a robust dictionary and a double collaborative constraint regular term, and can be used for geological exploration, urban planning, and space exploration. Background Technique

[0002] Hyperspectral anomaly detection aims to find pixels in a hyperspectral image that have significant spectral differences from neighboring background pixels. Compared with hyperspectral target detection, since hyperspectral anomaly detection does not require prior knowledge, it is widely used in many fields, such as geological exploration, on-site reconnaissance, urban planning, space exploration, etc.

[0003] Currently, the low-rank and sparse representation model based on the joint dictionary has developed rapidly in the field of hyperspectral anomaly detection. Among them, PAB-DC [Reference: Huyan Ning, Zhang Xiangrong, Zhou Huiyu, Jiao Licheng. Hyperspectral anomaly detection via background and potential anomaly dictionaries construction [J]. IEEE Transactions on Geoscience and Remote Sensing, 2018, 57(4): 2263-2276] is a classic work in this field. This method first decomposes the hyperspectral image into a low-rank component, a sparse component, and a noise component; secondly, it uses the joint sparse representation model to select background pixels to construct a background dictionary. At the same time, it defines the pixel anomaly level using the residual of the joint sparse representation model in the local area to construct a potential anomaly dictionary; finally, it optimizes and solves the constructed model. However, the above method has two problems: 1) Only global information is considered, ignoring local information, resulting in suboptimal detection performance; 2) The constructed potential anomaly dictionary is easily interfered by background pixels, that is, the robustness is not strong enough, which to a certain extent restricts the characterization effect of the anomaly component. Summary of the Invention

[0004] The object of the present invention is to address the deficiencies of the above-mentioned existing technologies by proposing a hyperspectral anomaly detection method that combines a robust dictionary and a double collaborative constraint regularization term, aiming to solve the problems of poor robustness during dictionary construction and suboptimal detection results caused by the failure to utilize global and local information. First, a density estimation model and the local outlier factor are respectively used to construct a robust dictionary, obtaining a potential anomaly dictionary and a background dictionary. Secondly, a low-rank and sparse representation model based on a double collaborative constraint regularization term is constructed and optimized to obtain an initial detection result. Finally, the initial detection result is nonlinearly transformed using a background suppression map and an anomaly enhancement map to obtain the final detection result. The combined dictionary constructed by the present invention has strong robustness and fully considers global and local characteristics, which is beneficial to improving the hyperspectral anomaly detection effect.

[0005] The specific steps for the present invention to achieve the above object are as follows:

[0006] (1) Input the hyperspectral image where H, W, and d respectively represent the height, width, and number of bands of the hyperspectral image;

[0007] (2) Use a density estimation model and the local outlier factor to evaluate all pixels in the hyperspectral image to obtain the first density map D m1 and the second density map D m2 ; and perform pixel-by-pixel multiplication and pixel-by-pixel addition operations on D m1 and D m2 respectively to obtain the background suppression map D1 and the anomaly enhancement map D2;

[0008] (3) Use the density estimation model to cluster the hyperspectral image where each cluster follows a Gaussian distribution, and record the position index of the cluster as where i ∈ [1, n M represents the index of the cluster, and n M represents the total number of clusters; at the positions of the background suppression map D1 and the anomaly enhancement map D2, perform binary processing on them using the 3σ principle to obtain the first index map B1 and the second index map B2; i

[0009] (4) Use the first index map B1 and the second index map B2 to construct the potential anomaly dictionary X (4) Use the first index map B1 and the second index map B2 to construct the potential anomaly dictionary X A and the background dictionary X B ;

[0010] (5) First, convert the hyperspectral image into a two-dimensional image where N = H × W; secondly, according to the potential anomaly dictionary X A and the background dictionary X B, decompose the two-dimensional form image X into a low-rank component X B W, a sparse component X A S, and a noise component E:

[0011] X = X B W + X A S + E,

[0012] where W and S respectively represent the representation coefficient matrices of the low-rank component and the sparse component;

[0013] (6) Take the decomposition formula of X in step (5) as a constraint condition to construct an optimization objective function for hyperspectral anomaly detection based on a double collaborative constraint regular term:

[0014]

[0015] where α, β, λ, and γ respectively represent the first, second, third, and fourth regularization coefficients, s.t. represents the constraint condition, ||·|| * represents the nuclear norm of the matrix, ||·||1 represents the L1 norm of the matrix, ||·|| 2,1 represents the L 2,1 norm of the matrix, represents the square of the Frobenius norm of the matrix;

[0016] (7) Use the alternating direction multiplier method to solve the objective function to obtain W, S, E. By solving the L2 norm of each column vector of X A S and converting it into a two-dimensional form, obtain the initial detection result D3;

[0017] (8) Use the background suppression map D1 and the anomaly enhancement map D2 to perform a non-linear transformation on the initial detection result D3 to obtain the final detection result D F :

[0018]

[0019] where τ1 and τ2 represent the transformation coefficients.

[0020] Compared with the prior art, the present invention has the following advantages:

[0021] First, the present invention uses a density estimation model and a local outlier factor to generate a binary map, and at the same time combines principal component analysis and superpixel segmentation to establish a joint dictionary construction strategy based on a two-stream density map. By constructing a robust joint dictionary through the above strategy, it can better characterize the background component and the anomaly component;

[0022] Second, since the present invention imposes a double collaborative constraint regularization term on the low-rank and sparse representation model, the ability of the model to depict the background component and the abnormal component is further improved by comprehensively utilizing the global information and the local information, thereby effectively improving the detection performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] Figure 1 is a flowchart for implementing the method of the present invention;

[0024] Figure 2 is a schematic diagram of the application of the method of the present invention in a specific example;

[0025] Figure 3 is a schematic diagram of different data sets used in the present invention; where (a) represents the false color image and its label of the Pavia data set; (b) represents the false color image and its label of the HYDICE data set; (c) represents the false color image and its label of the Gulfport data set;

[0026] Figure 4 is a comparison chart of the detection results of the method of the present invention and the mainstream methods on 3 data sets; where (I-a) to (I-h) respectively represent the detection results of RX, CRD, LRASR, LSMAD, PAB-DC, LSDM-MoG, KIFD and the method of the present invention on the Pavia data set, (I-i) represents the label of the Pavia data set; (II-a) to (II-h) respectively represent the detection results of RX, CRD, LRASR, LSMAD, PAB-DC, LSDM-MoG, KIFD and the method of the present invention on the HYDICE data set, (II-i) represents the label of the HYDICE data set; (III-a) to (III-h) respectively represent the detection results of RX, CRD, LRASR, LSMAD, PAB-DC, LSDM-MoG, KIFD and the method of the present invention on the Gulfport data set, (III-i) represents the label of the Gulfport data set;

[0027] Figure 5 is the ROC-(P d , P f ) curve graph of the method of the present invention and the mainstream methods on 3 data sets; where (a), (b), (c) respectively represent the ROC curves of the method of the present invention and the mainstream methods on the Pavia, HYDICE and Gulfport data sets with the abscissa being P f and the ordinate being P d ;

[0028] Figure 6 is the ROC-(P f, τ) curve; where (a), (b), and (c) respectively represent the ROC curves with τ as the abscissa and P as the ordinate of the method of the present invention and the mainstream methods on the Pavia, HYDICE, and Gulfport datasets f ROC curve. Specific implementation manner

[0029] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0030] Embodiment 1:

[0031] Refer to Figure 1 and Figure 2 A hyperspectral anomaly detection method integrating a robust dictionary and a double collaborative constraint regular term proposed by the present invention is specifically implemented according to the following steps:

[0032] Step 1: Input a hyperspectral image where H, W, and d respectively represent the height, width, and number of bands of the hyperspectral image;

[0033] Step 2: Use the density estimation model and the local outlier factor to evaluate the density values of all pixels in the hyperspectral image to obtain the first density map D m1 and the second density map D m2 ; and perform pixel-by-pixel multiplication and pixel-by-pixel addition operations on D m1 and D m2 respectively to obtain the background suppression map D1 and the anomaly enhancement map D2; the density estimation model used in this embodiment is a neural network model composed of 3 fully connected layers, and the corresponding number of nodes is {d, 128, n M}, where n M is set to 8. After the first 2 fully connected layers, the tanh activation function is used for nonlinear transformation, and a dropout layer is added after the activation function to prevent the model from overfitting. After the last fully connected layer, the softmax function is used to normalize the output to between 0 and 1.

[0034] Use the density estimation model to evaluate the density values of all pixels in the hyperspectral image to obtain the first density map D m1 , and the specific steps are as follows:

[0035] (2a1) Use the negative of the maximum likelihood function of the Gaussian mixture model as the loss function of the density estimation model, and use to represent the probability density function of the Gaussian distribution, where μ' and Σ' respectively represent the mean vector and the covariance matrix, and |·| represents the determinant operation, to obtain the probability density function f M (x):

[0036]

[0037] Among them, n M represents the number of Gaussian mixture components, and ω q represents the weighting coefficient of the Gaussian mixture components, and μ' q and Σ' q respectively represent the mean vector and covariance matrix corresponding to the q-th Gaussian component;

[0038] (2a2) The maximum likelihood function of the Gaussian mixture model is expressed as:

[0039]

[0040]

[0041] Among them, N = H·W represents the total number of pixels in the hyperspectral image, and x j represents the j-th pixel in X;

[0042] (2a3) The loss function of the density estimation model is expressed as follows:

[0043]

[0044]

[0045]

[0046]

[0047] Among them, ο jq represents the q-th component output after the j-th sample is sent into the model;

[0048] (2a4) Set hyperparameters, including the optimizer, learning rate, batch size, and maximum number of iterations, and train the density estimation model to obtain a trained density estimation model;

[0049] (2a5) Send the hyperspectral image into the trained density estimation model to obtain the output Use to calculate the density values of all pixels in, and convert it into a two-dimensional form to obtain the first density map D m1 .

[0050] Use the local outlier factor to evaluate the density values of all pixels in the hyperspectral image to obtain the second density map D m2 The steps are as follows:

[0051] (2b1) For the pixel x to be measured t(t ∈ [1, N]), let k-dist(x t ) and N k (x t ) represent the k-th pixel closest to x t and the set of k pixels respectively;

[0052] (2b2) For any pixel in N k (x t ), denote the reachable distance from x t to as: where

[0053]

[0054] represents the Euclidean distance between x and t ;

[0055] (2b3) Use the pixels in N k (x t ) to calculate the local reachability density lrd t (x k ): t where |N

[0056]

[0057] k (x t )| represents the number of pixels in N k t (x t ).

[0058] (2b4) Calculate the local outlier factor LOF(x t ) of x t :

[0059]

[0060] Repeat steps (2b1) to (2b4), traverse all pixels in the hyperspectral image , obtain the density values of all pixels, and convert them into a two-dimensional form to obtain the second density map D m2 .

[0061] Step 3: Use the density estimation model to cluster the hyperspectral image , and each cluster follows a Gaussian distribution. Record the position index of the cluster as where i ∈ [1, n M represents the index of the cluster, and n Mdenotes the total number of clusters; at the positions of the background suppression map D1 and the anomaly enhancement map D2, the 3σ principle is used to perform binarization processing on them respectively, obtaining the first index map B1 and the second index map B2. The above-mentioned use of the density estimation model for the hyperspectral image i to perform clustering, and its implementation is as follows: The hyperspectral image is sent into the trained density estimation model to obtain the output Using to find the category of the j-th sample, where argmax(·) represents the index of the maximum value. Traverse all rows in O in the above manner to obtain the clustering result. Here, the 3σ principle is used for binarization processing. Specifically, the positions where the values are between (μ0 - 3×σ, μ0 + 3×σ) are set to 0, and the remaining positions are set to 1, where μ0 represents the mean of each cluster and σ represents the standard deviation of each cluster. Traverse all clusters in the above manner to obtain the first index map B1 and the second index map B2.

[0062] Step 4: Use the first index map B1 and the second index map B2 to construct the potential anomaly dictionary X A and the background dictionary X B ; The implementation steps are as follows:

[0063] (4.1) Multiply the first index map B1 and the second index map B2 pixel by pixel to obtain the index map product Record the positions where the value in is 1 as P A , and select the pixels in the hyperspectral image with the same positions as P A to construct the potential anomaly dictionary X A ;

[0064] (4.2) Use principal component analysis to perform dimensionality reduction on the hyperspectral image , take the first 3 principal components, and denote them as

[0065] (4.3) Use the simple linear iterative clustering algorithm to perform superpixel segmentation on to obtain the superpixel block set where i' ∈ [1, n s represents the index of the superpixel block, and n s represents the total number of superpixel blocks;

[0066] (4.4) Add the first index map B1 and the second index map B2 pixel by pixel to obtain the index map sum Taking the superpixel block as a unit, calculate the sum of the elements in with the same positions as . If the sum is 0, record The centroid; otherwise, ignore Using P B to represent the set of centroids of the recorded superpixel blocks, and selecting pixels with the same position as P in B to construct the background dictionary X B .

[0067] Step 5: First, convert the hyperspectral image into a two-dimensional image where N = H × W; Second, according to the potential anomaly dictionary X A and the background dictionary X B , decompose the two-dimensional image X into a low-rank component X B W, a sparse component X A S, and a noise component E:

[0068] X = X B W + X A S + E,

[0069] where W and S represent the representation coefficient matrices of the low-rank component and the sparse component, respectively;

[0070] Step 6: Using the decomposition formula of X in Step 5 as a constraint condition, construct an optimization objective function for hyperspectral anomaly detection based on a double collaborative constraint regularization term:

[0071]

[0072] where α, β, λ, and γ represent the first, second, third, and fourth regularization coefficients, s.t. represents the constraint condition, ||·|| * represents the nuclear norm of the matrix, ||·||1 represents the L1 norm of the matrix, ||·|| 2,1 represents the L 2,1 norm of the matrix, represents the square of the Frobenius norm of the matrix;

[0073] Step 7: Use the alternating direction method of multipliers to solve the objective function to obtain W, S, E. By solving the L2 norm of each column vector of X A S and converting it into a two-dimensional form, obtain the initial detection result D3.

[0074] Using the alternating direction method of multipliers to solve the objective function, first introduce auxiliary variables to transform the objective function; then, according to the transformed objective function, construct an augmented Lagrangian function; and then use the alternating direction method of multipliers to optimize and solve to obtain the optimal objective. After W, S, E are obtained, the detection result is obtained according to the following formula:

[0075]

[0076] Among them, b ∈ [1, d] represents the b-th band of the hyperspectral image , j' represents the pixel index; after traversing all pixels in X, the detection result is converted into a two-dimensional form to obtain the initial detection result D3.

[0077] Step 8: Perform a non-linear transformation on the initial detection result D3 using the background suppression map D1 and the anomaly enhancement map D2 to obtain the final detection result D F :

[0078]

[0079] Among them, τ1 and τ2 represent transformation coefficients.

[0080] Example 2: The overall steps of the hyperspectral anomaly detection method in this example are the same as those in Example 1. For the training and testing details of the density estimation model, the following further description is made:

[0081] (a) Training of the density estimation model

[0082] In the training stage, the density estimation model uses the negative of the maximum likelihood function of the Gaussian mixture model as the loss function. Let represent the probability density function of the Gaussian distribution, where μ' and Σ' represent the mean vector and covariance matrix respectively, and |·| represents the determinant operation. Then the probability density function f M (x) of the Gaussian mixture model can be expressed as:

[0083]

[0084] where n M represents the number of Gaussian mixture components, ω q represents the weighting coefficient of the Gaussian mixture components, μ' q and Σ' q represent the mean vector and covariance matrix corresponding to the q-th Gaussian component respectively. Correspondingly, the maximum likelihood function of the Gaussian mixture model can be expressed as:

[0085]

[0086]

[0087] where N = H·W represents the total number of pixels in the hyperspectral image, and x j represents the j-th pixel in X.

[0088] According to the above analysis, the loss function of the density estimation model can be expressed as:

[0089]

[0090]

[0091]

[0092]

[0093] wherein, ο jq represents the q-th component output after the j-th sample is fed into the model. In the training stage, in this embodiment, the optimizer is set to Adagrad, the learning rate is 0.0001, the batch size is the total number of hyperspectral pixels, and the number of iterations is 1000.

[0094] (b) Testing of the density estimation model: After the training is completed, the hyperspectral image X is fed into the trained density estimation model to obtain the output Using calculate the density values of all pixels in X, and convert them into a two-dimensional form to obtain the density map D m1 .

[0095] Example 3: The overall steps of the hyperspectral anomaly detection method in this embodiment are the same as those in Example 1. The specific process of using the alternating direction method of multipliers to solve the objective function, that is, to optimize the objective, is further described as follows:

[0096] (7.1) Introduce auxiliary variables P, Q, R, and V, and transform the objective function into the following form:

[0097]

[0098] s.t. X = X B W + X A S + E, W = P, S = Q, W = R, S = V

[0099] (7.2) Construct the augmented Lagrangian function:

[0100]

[0101] wherein, Y m (m ∈ [1, 5]) represents the Lagrange multiplier, tr[·] represents the trace of the matrix, and η represents the penalty coefficient;

[0102] (7.3) Use the alternating direction method of multipliers for optimization, that is, update one variable each time and fix the remaining variables:

[0103] (7.3a) Update P and fix the remaining variables, then the objective function is transformed into:

[0104]

[0105] Φ1 / η (Σ) i”,i” = max(Σ i”,i” -1 / η, 0)

[0106] where U, Σ, V T represent the left singular matrix, diagonal matrix, and right singular matrix obtained from the singular value decomposition of W + Y2 / η respectively, and Σ i”,i” represents the i-th main diagonal element of the diagonal matrix Σ;

[0107] (7.3b) Update Q, fix the other variables, then the objective function is transformed into:

[0108]

[0109] where soft1(z, u) = sign(z)·max(|z| - u, 0), and sign(·) represents the sign function, i.e., r Q ∈ [1, n] and c Q ∈ [1, N] represent the row index and column index of the variable respectively.

[0110] (7.3c) Update E, fix the other variables, then the objective function is transformed into:

[0111]

[0112] where ||y||2 represents the L2 norm of the vector y, and c E ∈ [1, N] represents the column index of the variable respectively;

[0113] (7.3d) Update R, fix the other variables, then the objective function is transformed into:

[0114]

[0115] (7.3e) Update V, fix the other variables, then the objective function is transformed into:

[0116]

[0117] (7.3f) Update W, fix the other variables, then the objective function is transformed into:

[0118]

[0119] (7.3g) Update S, fix the other variables, then the objective function is transformed into:

[0120]

[0121] (7.3h) Update the Lagrange multiplier Y m (m ∈ [1, 5]):

[0122] Y1 = Y1 + η(X - X B W - X A S - E)

[0123] Y2 = Y2 + η(W - P)

[0124] Y3 = Y3 + η(S - Q)

[0125] Y4 = Y4 + η(W - R)

[0126] Y5 = Y5 + η(S - V)

[0127] (7.3i) Update the penalty coefficient η:

[0128] η = min(ρη, η max )

[0129] Where P, Q, E, R, V, W, S, Y m (m ∈ [1, 5]) have initial values of all zero matrices, η0, ρ, and η max are 10 -4 , 1.1, and 10 10 ;

[0130] Repeat steps (7.3a) to (7.3i) until the iteration termination condition is satisfied:

[0131] max(||X - X B W - X A S - E|| F , ||W - P|| F , ||S - Q|| F , ||W - R|| F , ||S - V|| F ) < ε,

[0132] Where ε represents the iteration convergence value, and its value is 10 -6 ; If the iteration does not converge, terminate the iteration when the number of iterations exceeds 500 times; obtain the optimized objectives W, S, E.

[0133] The effects of the present invention will be further described below in combination with experiments.

[0134] 1. Experimental conditions:

[0135] The experiments of the present invention were carried out in a hardware environment with a CPU main frequency of 2.00 GHz, 16 GB of memory, and a GeForce RTX 2080Ti graphics card, and a software environment of Windows 10 operating system, Python 3.6.13, Pytorch-GPU 1.8.2, and CUDA 10.2.

[0136] 2. Experimental content:

[0137] In this experiment, hyperspectral images were acquired as a dataset using 3 different sensors. This dataset includes the Pavia (acquired under the ROSIS-03 sensor), HYDICE (acquired under the HYDICE sensor), and Gulfport (acquired under the AVIRIS sensor) datasets, as shown in (a), (b), and (c) of Figure 3 respectively. On each dataset, the detection effects of the method of the present invention and seven mainstream methods were compared from two perspectives: qualitative (including detection maps and two ROC curves) and quantitative (including two AUC values).

[0138] Among them, the mainstream methods mainly include RX [Paper: Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution], CRD [Paper: Collaborative representation for hyperspectral anomaly detection], LRASR [Paper: Anomaly detection in hyperspectral images based on low-rank and sparse representation], LSMAD [Paper: A low-rank and sparse matrix decomposition-based mahalanobis distance method for hyperspectral anomaly detection], PAB-DC [Paper: Hyperspectral anomaly detection via background and potential anomaly dictionaries construction], LSDM-MoG [Paper: Low-rank and sparse decomposition with mixture of Gaussian for hyperspectral anomaly detection], and KIFD [Paper: Hyperspectral anomaly detection with kernel isolation forest].

[0139] 3. Simulation results and analysis:

[0140] Figure 4 This is a comparison chart of the detection results of the method of the present invention and the mainstream methods on 3 datasets; among them, (I)-(III) correspond to the Pavia, HYDICE, and Gulfport datasets respectively, (a)-(g) correspond to 7 existing methods including RX, CRD, LRASR, LSMAD, PAB-DC, LSDM-MoG, and KIFD respectively, (h) corresponds to the method of the present invention, and (i) corresponds to the label; specifically: Figure 4In (I-a) to (I-h), they respectively represent the detection results of RX, CRD, LRASR, LSMAD, PAB-DC, LSDM-MoG, KIFD, and the method of the present invention on the Pavia dataset, and (I-i) represents the label of the Pavia dataset; (II-a) to (II-h) respectively represent the detection results of RX, CRD, LRASR, LSMAD, PAB-DC, LSDM-MoG, KIFD, and the method of the present invention on the HYDICE dataset, and (II-i) represents the label of the HYDICE dataset; (III-a) to (III-h) respectively represent the detection results of RX, CRD, LRASR, LSMAD, PAB-DC, LSDM-MoG, KIFD, and the method of the present invention on the Gulfport dataset, and (III-i) represents the label of the Gulfport dataset.

[0141] By comparing the existing methods, the method of the present invention with the labels, it can be seen that the method of the present invention has the best effect and is significantly better than the mainstream methods among the existing 7 kinds of prior arts.

[0142] Figure 5 and Figure 6 respectively represent the ROC-(P d ,P f ) curve and the ROC-(P f ,τ) curve of the method of the present invention and the mainstream methods on 3 datasets; among them, (a) represents the Pavia dataset, (b) represents the HYDICE dataset, and (c) represents the Gulfport dataset.

[0143] The larger the area under the ROC-(P d ,P f ) curve, the better the detection performance. It can be seen from Figure 5 that the detection effect of the method of the present invention is the best; the smaller the area under the ROC-(P f ,τ) curve, the lower the false alarm rate. It can be seen from Figure 6 that the method of the present invention has the lowest false alarm rate and the best performance.

[0144] The following table (Table 1) lists the AUC values (i.e., AUC-(P d ,P f ) and AUC-(P f ,τ)) of the method of the present invention and 7 mainstream methods on 3 datasets. The AUC value is between 0 and 1. For AUC-(P d ,P f ), the closer its value is to 1, the better the detection effect; for AUC-(P f ,τ), the closer its value is to 0, the lower the false alarm rate, that is, the better the performance.

[0145] Table 1 Comparison of AUC values of the method of the present invention and seven mainstream methods on 3 datasets

[0146]

[0147] As can be seen from the above table, on 3 different datasets (Pavia, HYDICE, and Gulfport), the method of the present invention has the largest AUC-(P d ,P f ) and the smallest AUC-(P f ,τ), further indicating that the method of the present invention has the highest detection accuracy and the lowest false alarm rate compared with the existing methods, that is, the best performance.

[0148] The above experimental results prove the correctness and effectiveness of the method proposed by the present invention.

[0149] The parts not described in detail in the present invention belong to the common general knowledge of those skilled in the art.

[0150] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Obviously, for those skilled in the art, after understanding the content and principle of the present invention, various modifications and changes in form and details may be made without departing from the principle and structure of the present invention. However, these modifications and changes based on the idea of the present invention are still within the scope of the claims of the present invention.

Claims

1. A hyperspectral anomaly detection method integrating a robust dictionary and a double collaborative constraint regular term, characterized in that It includes the following steps: (1) Input hyperspectral image where H, W, and d represent the height, width, and number of bands of the hyperspectral image, respectively; (2) Use the density estimation model and the local outlier factor to evaluate all pixels in the hyperspectral image to obtain the density values of all pixels, and get the first density map D m1 and the second density map D m2 ; and perform pixel-by-pixel multiplication and pixel-by-pixel addition operations on D m1 and D m2 respectively to obtain the background suppression map D1 and the anomaly enhancement map D2; (3) Using density estimation model to estimate hyperspectral images Clustering is performed, each cluster obeys Gaussian distribution, and the location index of the cluster is recorded as where i∈[1,n M ] represents the cluster index, n M Indicates the total number of clusters; Po in the background suppression map D1 and the abnormal enhancement map D2 i In terms of position, the 3σ principle is used to perform binarization processing respectively, and the first index map B1 and the second index map B2 are obtained; (4) Construct the potential anomaly dictionary X using the first index map B1 and the second index map B2 A and the background dictionary X B ; (5) First, convert the hyperspectral image into a two-dimensional image where N = H × W; secondly, according to the potential anomaly dictionary X A and the background dictionary X B , decompose the two-dimensional image X into a low-rank component X B W, a sparse component X A S, and a noise component E: X = X B W + X A S + E, Among them, W and S respectively represent the representation coefficient matrices of the low-rank component and the sparse component; (6) Taking the decomposition formula of X in step (5) as a constraint condition, construct an optimization objective function for hyperspectral anomaly detection based on a dual collaborative constraint regular term: s.t. X = X B W + X A S + E where α, β, λ, and γ represent the first, second, third, and fourth regularization coefficients respectively, s.t. represents the constraint condition, ||·|| * represents the nuclear norm of a matrix, ||·||1 represents the L1 norm of a matrix, ||·|| 2,1 represents the L 2,1 norm, represents the square of the Frobenius norm of a matrix; (7) Solve the objective function using the alternating direction method of multipliers to obtain W, S, and E. By solving the L2 norm for each column vector of X A and converting it into a two-dimensional form, the initial detection result D3 is obtained; (8) Perform a non - linear transformation on the initial detection result D3 using the background suppression map D1 and the anomaly enhancement map D2 to obtain the final detection result D F : Among them, τ1 and τ2 represent transformation coefficients.

2. The method according to claim 1, wherein: In step (2), the density estimation model is a neural network model composed of three fully connected layers, with the corresponding number of nodes being {d, 128, n M}, where n M is set to 8. The first two fully connected layers are followed by a tanh activation function for non-linear transformation. A dropout layer is added after the activation function to prevent overfitting of the model. After the last fully connected layer, a softmax function is used to normalize the output between 0 and 1.

3. The method according to claim 1, wherein: The first density map D in step (2) m1 , is obtained according to the following steps: (2a1) Use the negative of the maximum likelihood function of the Gaussian mixture model as the loss function of the density estimation model. With representing the probability density function of the Gaussian distribution, where μ' and Σ' represent the mean vector and covariance matrix respectively, and |·| represents the determinant operation, the probability density function f of the Gaussian mixture model is obtained M (x): where n M denotes the number of Gaussian mixture components, ω q denotes the weighted coefficient of the Gaussian mixture components, μ' q and Σ' q denote the mean vector and covariance matrix corresponding to the q-th Gaussian component, respectively; (2a2) The maximum likelihood function of the Gaussian mixture model is expressed as: where N = H·W represents the total number of pixels in the hyperspectral image, and x j represents the j-th pixel in X; (2a3) The loss function of the density estimation model is expressed as follows: where, ο jq represents the q-th component output after the j-th sample is fed into the model; (2a4) Set hyperparameters, including the optimizer, learning rate, batch size, and maximum number of iterations, and train the density estimation model to obtain a trained density estimation model; (2a5) Feed the hyperspectral image into the trained density estimation model to obtain the output Use to calculate the density values of all pixels in m1 .

4. The method according to claim 1, wherein: The second density map D in step (2) m2 , is obtained according to the following steps: (2b1) For the pixel x to be measured t (t ∈ [1, N]), let k-dist(x t ) and N k (x t ) represent the k-th pixel closest to x t and the set of k pixels respectively; (2b2) for N k (x t ) any pixel Let x t to reachable distance be expressed as: Among them, represents the Euclidean distance t between x and (2b3) Use N k (x t ) to calculate the local reachability density lrd t of x k (x t ): Among them, |N k (x t )| represents the number of pixels in N k (x t ); (2b4) Calculate x t 's Local Outlier Factor LOF(x t ): Repeat steps (2b1) to (2b4) to traverse the hyperspectral image for all pixels therein, obtain the density values of all pixels, and convert them into a two-dimensional form to obtain the second density map D m2 .

5. The method according to claim 1, characterized in that: In step (3), the hyperspectral image is clustered using the density estimation model as follows: The hyperspectral image is fed into the trained density estimation model to obtain the output Using to find the class of the j-th sample, where argmax(·) represents the index of the maximum value. Traverse all rows in O in the above manner to obtain the clustering result.

6. The method according to claim 1, characterized in that: In step (3), for binarization processing, positions where the numerical value is between (μ0 - 3×σ, μ0 + 3×σ) are set to 0, and the remaining positions are set to 1, where μ0 represents the mean of each cluster and σ represents the standard deviation of each cluster. Traverse all clusters in the above manner to obtain the first index map B1 and the second index map B2.

7. The method according to claim 1, characterized in that: In step (4), the first index map B1 and the second index map B2 are used to construct the potential anomaly dictionary X A and the background dictionary X B , and the implementation steps are as follows: (4.1) Multiply the first index map B1 and the second index map B2 pixel by pixel to obtain the index map product Record The positions with a value of 1 in are P A and in the hyperspectral image Select the pixels with the same positions as P A to construct the potential anomaly dictionary X A ; (4.2) Use principal component analysis to reduce the dimensionality of the hyperspectral image and take the first three principal components, which are denoted as (4.3) Use the simple linear iterative clustering algorithm to perform superpixel segmentation to obtain a set of superpixel blocks where i' ∈ [1, n s represents the index of the superpixel block, and n s represents the total number of superpixel blocks; (4.4) Add the first index map B1 and the second index map B2 pixel by pixel to obtain the sum of the index maps Taking superpixel blocks as units, calculate the sum of elements with the same positions as in. If the sum is 0, record the centroid of; otherwise, ignore Denote the set of centroids of superpixel blocks that have been recorded by P B and select pixels with the same positions as P in to construct the background dictionary X B . B .

8. The method according to claim 1, characterized in that: In step (7), the optimization of the objective function using the alternating direction method of multipliers is as follows: (7.1) Introduce auxiliary variables to transform the objective function; (7.2) Construct an augmented Lagrangian function according to the transformed objective function; (7.3) Use the alternating direction method of multipliers to optimize and solve the augmented Lagrangian function to obtain W, S, and E.

9. The method according to claim 1, characterized in that: In step (7), by solving the L2 norm for each column vector of X A S, the formula is as follows: where b ∈ [1, d] represents the b-th band of the hyperspectral image After traversing all pixels in X, the detection results are converted into a two-dimensional form to obtain the initial detection result D3.

Images