A filter and low-rank decomposition based spatial-spectral joint hyperspectral image anomaly detection method

Abstract

The application relates to an abnormality detection method based on a hyperspectral image. The main body is based on a space-spectrum combined feature extraction method of filtering and low-rank decomposition to perform abnormality detection on the hyperspectral image. The specific method comprises the following steps: firstly, in the spatial dimension, a reduced dimension image is obtained through a data dimension reduction and eigenvalue weighted fusion method, and then an improved spatial filtering method is used to extract the spatial features of the image to obtain an initial spatial feature image. In the spectral dimension, a background reconstruction image of the approximate background is obtained by using a Tucker decomposition method on the original hyperspectral image, and a background dictionary of the image is obtained by using an improved k-means clustering method, then the background dictionary is input into a low-rank decomposition model to obtain a sparse matrix, and an initial spectral feature image is obtained, finally, the initial spectral feature image is fused with the spatial feature image to realize abnormality detection.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing image processing technology, and specifically relates to an anomaly detection method for space-spectral joint hyperspectral images based on filtering and low-rank decomposition. Background Technology

[0002] Hyperspectral remote sensing images integrate spatial and spectral features, enabling them to better reflect the differences between various ground features, thus possessing advantages unmatched by traditional optical images. Hyperspectral image anomaly detection, as an unsupervised target detection algorithm, can detect anomalous targets by mining differences in the spectral curves of ground features without requiring any prior spectral information. Therefore, it has high application value and has become one of the research hotspots in the field of remote sensing technology, widely applied in various sectors of the national economy. Traditional anomaly detection algorithms only extract single or a small number of features from the image to complete anomaly detection, making it difficult to comprehensively describe the image's features, thus limiting the improvement of algorithm detection performance. Summary of the Invention

[0003] To overcome the shortcomings of the prior art, the present invention aims to provide an anomaly detection method for spatial-spectral joint hyperspectral images based on filtering and low-rank decomposition, so as to solve the problem that existing numerical algorithms only utilize spatial or spectral information in hyperspectral images, making it difficult to fully and comprehensively describe the true information of the image, and that the background dictionary is easily contaminated during low-rank decomposition.

[0004] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0005] A method for anomaly detection in space-spectrum joint hyperspectral images based on filtering and low-rank decomposition includes the following steps:

[0006] Step 1, process the original hyperspectral image H∈R M×N×B Dimensionality reduction and weighted fusion are used to obtain the dimensionality-reduced image R. pca ∈R M×N Where M and N represent the spatial dimensions of the image, and B represents the spectral dimension;

[0007] Step 2, extract R pca From the spatial information in the image, obtain the initial spatial feature map R. spa ∈R M×N ;

[0008] Step 3: Obtain the hyperspectral background reconstruction image R from the original hyperspectral image H through Tucker decomposition. bg ∈R M×N×B ;

[0009] Step 4, based on R bgRepresentative pixels are selected from the separated background information to form the background dictionary D in the low-rank decomposition model;

[0010] Step 5: Perform low-rank decomposition on the original hyperspectral image H using the background dictionary D to obtain the initial spectral feature map R. spe ∈R M×N ;

[0011] Step 6, R spa and R spe Perform dot product fusion to obtain the final detection result image R. result ∈R M×N .

[0012] Compared with existing technologies, this invention fully explores the spatial features, spectral features, low rank and sparsity of hyperspectral images, providing a more comprehensive description of the images. This allows for the integrated use of the advantages of each feature image to achieve anomaly detection and improve image anomaly detection performance. Attached Figure Description

[0013] Figure 1 This is a pseudo-color image of the original hyperspectral image of this invention.

[0014] Figure 2 This is the ground truth map for the abnormal target.

[0015] Figure 3 This is a flowchart of an example of the present invention.

[0016] Figure 4 The following is a comparison of anomaly detection results in an example of the present invention, wherein (a) is a pseudo-color image, (b) is the true distribution map of the target, (c) is RX, (d) is LRX, (e) is CRD, (f) is WSCF, (g) is RPCA-RX, (h) is LRASR, and (i) is the algorithm of the present invention.

[0017] Figure 5 This is a three-dimensional curve comparing the anomaly detection results in an example of the present invention.

[0018] Wherein (a) is a three-dimensional diagram of the true distribution of abnormal targets, (b) is a three-dimensional network diagram of the RX algorithm, (c) is the LRX algorithm, (d) is the CRD algorithm, (e) is the WSCF algorithm, (f) is the RPCA-RX algorithm, (g) is the LRASR algorithm, and (h) is the algorithm of this invention.

[0019] Figure 6 This is a ROC curve of the results of each anomaly detection algorithm in the examples of this invention.

[0020] (a) shows the ROC curves of the RX algorithm, LRX algorithm, CRD algorithm, WSCF algorithm, RPCA-RX algorithm, LRASR algorithm and the algorithm of this invention, and (b) shows a magnified view of the ROC curves of the above algorithms in the range of 0-0.2 false alarm rate. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0022] To demonstrate the effectiveness of the method of this invention, real hyperspectral images were used for anomaly detection. The experiment used real hyperspectral image data acquired by an airborne visible / infrared imaging spectrometer (AVIRIS) sensor, and the experimental scenario was a real airport, such as... Figure 1 As shown. The spatial resolution of the image is 3.5 m / pixel, and the spectral range is 0.37–2.5 μm. The image size used in the experiment was 100 × 100, with 189 bands, and three aircraft parked on the tarmac were considered anomalous targets. The ground truth image of the anomalous targets is shown below. Figure 2 As shown, the abnormal targets are distributed in the upper half of the image, while the remaining area is identified as the background area.

[0023] For experimental image data, refer to Figure 3 The basic process of the space-spectrum joint hyperspectral image anomaly detection method based on filtering and low-rank decomposition of this invention is as follows:

[0024] Step 1, process the original hyperspectral image H∈R M×N×B The dimensionality is reduced using algorithms such as PCA, and then weighted and fused to obtain the dimensionality-reduced image R. pca ∈R M×N Where M and N represent the spatial dimensions of the image, and B represents the spectral dimension.

[0025] In an embodiment of the present invention, the PCA algorithm is used to transform the original hyperspectral image H∈R M×N×B Reducing the dimension to p, the principal components are arranged according to their energy content, with the smallest principal components retaining the majority of image information. Therefore, the first few principal components can be selected to replace the original data for subsequent spatial feature extraction. The specific method is as follows:

[0026] Based on the covariance matrix, the p largest eigenvalues ​​are obtained as {λ1, λ2, ..., λ...} p The obtained principal components are then sorted in descending order of information content, where the p principal components can be represented as {H1, H2, ..., H...}. p}; In subsequent spatial feature extraction operations, only p images are subjected to a weighted average fusion operation based on the feature value magnitude. The mathematical expression is:

[0027]

[0028] In the formula, λ i Let R represent the i-th eigenvalue. pca The weighted fused image can replace the original image for subsequent processing, reducing the complexity of data processing and improving the efficiency of the algorithm.

[0029] In this embodiment, the images corresponding to the top three principal components containing the most information can be selected for weighted fusion.

[0030] Step 2, extract R pca From the spatial information in the image, obtain the initial spatial feature map R. spa ∈R M×N .

[0031] In this step, an improved spatial filtering algorithm is used to extract image R from the dimensionality-reduced image. pca The improved spatial filter incorporates spatial information and improves the filter weights based on the traditional bilateral filter. The traditional bilateral filter introduces weights based on the gray-level difference between pixels into the Gaussian filter, combining the characteristics of Gaussian and mean filtering, resulting in good edge preservation and noise reduction. Therefore, the improved spatial filter is based on the idea of ​​the traditional bilateral filter, but with improvements to the filter weights. Specifically, the improved spatial filtering algorithm makes the following improvements to the bilateral filter:

[0032] Step 2.1: The first distance weight of the improved spatial filter is determined by the Euclidean distance between the reconstruction center and its neighboring pixels. Specifically, the weight formula based on the Euclidean distance between pixels can be expressed as:

[0033]

[0034] In the formula, x is the center pixel of the template window, i.e., the pixel to be reconstructed, with coordinates (k, l), ε is the pixels surrounding the center pixel to be reconstructed, with coordinates (u, v), and d d σ represents the Euclidean distance between two points. d The distance weight parameter indicates that the closer a pixel is to the pixel to be reconstructed, the more similar it is to the pixel to be reconstructed, and the greater its contribution to the reconstruction of the center pixel.

[0035] Step 2.2: The gradient in an image represents the magnitude and direction of changes in pixel grayscale values. The smaller the gradient value, the smoother the image. Therefore, the second gradient weight of the improved spatial filter is determined by the gradient difference between pixel grayscale values. Specifically, the gradient weight formula of the improved spatial filter is:

[0036]

[0037] In the formula, σ r Let d be the gradient difference weight parameter. r Let G(x) and G(ε) represent the gradient difference between points x and ε. The mathematical formula for G(x) can be expressed as:

[0038]

[0039] In the formula, X represents the gray value of the pixel x to be reconstructed. The greater the difference in pixel gradient values ​​between two points, the smaller the weight coefficient.

[0040] Step 2.3: The improved spatial filter combines the distance weights and gradient weights in the spatial domain. The kernel of the improved spatial filter is obtained by multiplying the distance kernel by the gradient kernel. Specifically, the kernel of the improved spatial filter is represented as:

[0041] w bf (x,ε)=w d (x,ε)w r (x,ε)

[0042] In summary, the improved formula for calculating the spatial filter is:

[0043]

[0044] In the formula, ε i Let represent the other pixels in the convolution kernel whose center pixel is the pixel to be reconstructed, n be the number of pixels in the improved spatial filter convolution kernel, and k(x) represent the regularization parameter, expressed as:

[0045]

[0046] R spa (x) represents the pixel value at position x after spatial filtering. After calculating all pixels, the initial spatial feature map R is obtained. spa .

[0047] Table 1 below shows the parameter analysis of the experimental image data. From Table 1, we can see that when the window size r is 8, the weight parameter σ of the spatial filtering algorithm... d and σ r The algorithm performs best when the values ​​are 0.6 and 0.4 respectively.

[0048] Table 1

[0049]

[0050]

[0051] Step 3: Reconstruct the background of the original hyperspectral image H using Tucker decomposition to obtain the hyperspectral background reconstructed image R. bg ∈R M×N×B .

[0052] The process of this step can be described in detail as follows:

[0053] Step 3.1, for an Nth-order tensor X, its Tucker decomposition can be expressed as:

[0054] X = g × 1U (1) ×2U (2) …× N U (N)

[0055] In the formula, g represents the kernel tensor. Factor matrices representing different modes, The kernel tensor represents the degree of interaction between different dimensions.

[0056] In this embodiment, the hyperspectral data is represented as a third-order tensor, i.e., H X ∈R M×N×B If the image has M rows and N columns, and B bands, then the original hyperspectral image H∈R M×N×B Tucker decomposition, represented as:

[0057]

[0058] In the formula, A∈R M×P , B∈R N×U , C∈R B×R In the formula, g X ×1A represents the 1-mode tensor product of g and A, where g X ×2B represents the 2-mode tensor product, g X ×3C represents the 3-mode tensor product, where the size of the factor matrices A, B, and C under different modes are M×P, N×U, and B×R, respectively. X Let g represent the kernel tensor in this embodiment, where g∈R P×U×R g pur For g X The (p,u,r)th element, the p-th vector of matrix A is a p The u-th vector of matrix B is b u The r-th vector of matrix C is c r .

[0059] Step 3.2: In the low-rank representation model, the hyperspectral image is divided into background based on the low-rank matrix and anomalies based on the sparse matrix. Tucker decomposition can be used to reconstruct the background image. The third-order tensor Tucker decomposition can be solved using the iterative alternating least squares method, transforming the Tucker decomposition into an optimization problem:

[0060]

[0061] In the formula, I is the identity matrix, and the kernel tensor must satisfy the following condition:

[0062] g X =X×1A T ×2B T ×3C T

[0063] Therefore, the objective function of the Tucker decomposition can be optimized as follows:

[0064]

[0065] Step 3.3: The Tucker decomposition is optimized using a higher-order singular value decomposition algorithm. Factor matrices A, B, and C of the solution matrix are solved in each mode, and the kernel tensor g of the tensor projected in each mode is calculated. X The result {g} is obtained X The method is as follows: Sort A, B, and C in descending order based on their eigenvalues, select the first few eigenvalues ​​for each eigenvalue, and substitute their corresponding eigenvectors into g. X =X×1A T ×2B T ×3C T The kernel tensor g was calculated. X The reconstructed background image R is obtained through matrix multiplication. bg .

[0066] Step 4, based on R bg Representative pixels are selected from the separated background information to form the background dictionary D in the low-rank decomposition model.

[0067] Specifically, this step can employ an improved k-means clustering algorithm, which refines the strategy for determining cluster centers to separate the reconstructed background image R. bg To select representative pixels from n samples, the k most representative cluster centers are chosen. The specific steps are as follows:

[0068] Step 4.1: Randomly select the first cluster center. For the remaining n-1 samples, assign each sample x... i The minimum spatial Euclidean distance from all current cluster centers is denoted as D(x). i ), D(x) for all samplesi The maximum value in ) is D max (x i ).

[0069] Step 4.2: To achieve a more ideal clustering effect, the cluster centers should be relatively dispersed. That is, the next cluster center should be chosen from points farther away from all current cluster centers. Therefore, for the next cluster center, D(x) should be... i (less than D) max (x i The selection probability of a sample with a probability of 1 / 2 is reduced to zero, and the D(x) of the remaining samples in the dataset with a non-zero probability is reduced to zero. i The sum of these two terms is denoted as sum(D(x)). i Then D(x) i ) / sum(D(x i That is, sample x. i The probability of being selected as the next cluster center. After calculating the probabilities of all points, the roulette wheel method is used to determine the next cluster center.

[0070] Step 4.3: Repeat the minimum distance calculation and cluster center determination operations until the kth cluster center, i.e., k pixels, is selected to form the background dictionary D.

[0071] In this embodiment, 15 categories are selected in the simulation experiment, and 25 pixels are taken for each category.

[0072] Step 5: Perform low-rank decomposition on the original hyperspectral image using the obtained background dictionary D to obtain the initial spectral feature map R. spe ∈R M×N .

[0073] In this step, the process of performing low-rank decomposition on the original hyperspectral image H is as follows:

[0074] Step 5.1, convert the original hyperspectral image H∈R M×N×B Substituting into the low-rank representation model, we get

[0075] H = DZ + E

[0076] In the formula, Z represents the coefficient matrix and E represents the sparse matrix.

[0077] Step 5.2: Anomaly targets exhibit sparse properties; therefore, the anomaly target detection result map can be obtained by solving the sparse matrix E. In the low-rank representation model, to maximize the sparsity of the matrix, i.e., to minimize the number of non-zero elements in the sparse matrix, l can be used. 2,1 The norm constrains the sparse matrix E, and the objective function is expressed as:

[0078]

[0079] In the formula, rank(·) represents the rank function, λ represents the coefficients of the balancing terms, and ||·|| 2,1 Represent l 2,1 The norm is defined as the sum of the l2 norms of all columns of a matrix, E. 2,1 The norm can be expressed as:

[0080]

[0081] Step 5.3: Since the rank of Z in the above equation is non-convex, the optimal solution cannot be found. Therefore, an auxiliary variable W is introduced to transform the objective function into...

[0082]

[0083] In the formula ||·|| * This represents the calculation of the nuclear norm.

[0084] Step 5.4: Construct a Lagrange function to solve the optimization problem of the objective function, expressed as follows:

[0085]

[0086] In the formula, Y 1 Y 2 These are Lagrange multipliers, and μ > 0 represents the penalty parameter.

[0087] Step 5.5, solve for the Lagrange function, the process is as follows:

[0088] (1) Fix Z, E, Y 1 Y 2 Update variable W, and solve using the formula:

[0089]

[0090] In the formula, Θ represents the singular value thresholding operation;

[0091] (2) Fixed W, E, Y 1 Y 2 Update variable Z, and solve using the formula:

[0092] Z t+1 =(D T D+I) -1 (D T XD T E t +W t+1 +(D T Y t 1 +Y t 2 ) / μ t )

[0093] In the formula, k represents the number of iterations.

[0094] (3) Fix W, Z, Y 1 Y 2 Update variable E, and solve using the formula:

[0095]

[0096] In the formula, S represents soft threshold contraction.

[0097] (4) Fix W, Z, and E, and update the Lagrange multiplier (Y). 1 ,Y 2 The solution formula is:

[0098]

[0099]

[0100] Step 5.6: After solving the low-rank representation model of the original image H, the l2 norm of the sparse matrix E is calculated to obtain the spectral feature image R. spe ∈R M×N The formula is as follows:

[0101]

[0102] Where R spe (p) represents the spectral feature value at each pixel location in the image. Calculating all spectral feature values ​​yields the spectral feature image R. spe .

[0103] Step 6, process the spatial feature image R spa ∈R M×N and spectral feature image R spe ∈R M×N Perform dot product fusion to obtain the final detection result image R. result ∈R M×N The anomaly detection result is output, expressed by the formula:

[0104] R result =R spa ·R spe

[0105] In the formula, R result The final detection result image of the algorithm of this invention is R. spa With R spe These represent the initial spatial feature map and the initial spectral feature map, respectively.

[0106] Depend on Figure 4It can be seen that the RX, LRX, WSCF, and RPCA-RX algorithms can effectively detect the location of aircraft targets, but the background suppression effects of these four algorithms are generally poor, and false detections occur in the lower right corner of the image. The CRD algorithm produces the worst results, suppressing both the abnormal target and the background, making it impossible to detect the target's location visually. The LRASR algorithm shows good background suppression, making the target outline mostly visible, but some background remains. The algorithm of this invention has the best background suppression effect, with a clear target outline and virtually no missed or false detections. Subjective visual analysis of the detection results of each algorithm shows that the algorithm of this invention has the best detection effect, followed by the LRASR algorithm, and the CRD algorithm has the worst.

[0107] Depend on Figure 5 The 3D image of the detection results clearly shows that the RX, CRD, WSCF, RPCA-RX, and LRASR algorithms exhibited false detections in the upper left corner, with relatively large amplitudes. Among them, the RPCA-RX and LRASR algorithms showed better background suppression. The LRX algorithm's background suppression effect was average. The algorithm of this invention showed the best background suppression effect, highlighting abnormal targets with larger target peaks. Based on the above subjective evaluation analysis, for the experimental image data, the algorithm of this invention showed the best detection performance compared to the comparison algorithms.

[0108] Depend on Figure 6 It can be seen that when the false alarm rate is less than 0.05, the detection probability of algorithms other than the CRD algorithm increases rapidly. When the false alarm rate is constant, the ROC curve of the algorithm of this invention is relatively higher than that of the LRASR algorithm. Further analysis using the enlarged view of the upper left corner of the ROC curve... Figure 6 (b) It can be seen that the algorithm of the present invention first achieves the detection probability of 1, followed by LRASR, and the CRD algorithm has the worst detection performance on the experimental data, which shows the advantage of the algorithm of the present invention.

[0109] Table 2 shows the area under the ROC curve (AUC) values ​​of each algorithm in the examples of this invention.

[0110] Table 2

[0111]

[0112] As shown in Table 2, the CRD algorithm has the lowest AUC value, significantly lower than the other algorithms. The AUC values ​​of the RX and WSCF algorithms are almost equal. The AUC values ​​of LRX, LRASR, and the algorithm of this invention are higher, outperforming the other four comparison algorithms. The algorithm of this invention has the highest AUC value, reaching 0.9994. Combined with subjective evaluation analysis, it can be concluded that the algorithm of this invention has the best detection performance on the experimental image data.

Claims

1. A method for anomaly detection in space-spectrum joint hyperspectral images based on filtering and low-rank decomposition, characterized by: Includes the following steps: Step 1, process the original hyperspectral image H∈R M×N×B Dimensionality reduction and weighted fusion are used to obtain the dimensionality-reduced image R. pca ∈R M×N ;in M , N Represents the spatial dimension of an image. B Indicates spectral dimension; Step 2, extract R pca From the spatial information in the image, obtain the initial spatial feature map R. spa ∈R M×N ; Step 3: Obtain the hyperspectral background reconstruction image R from the original hyperspectral image H through Tucker decomposition. bg ∈R M×N×B ; Step 4, based on R bg Representative pixels are selected from the separated background information to form the background dictionary D in the low-rank decomposition model; Step 5: Perform low-rank decomposition on the original hyperspectral image H using the background dictionary D to obtain the initial spectral feature map R. spe ∈R M×N ; Step 6, R spa and R spe Perform dot product fusion to obtain the final detection result image R. result ∈R M×N .

2. The spatial-spectral joint hyperspectral image anomaly detection method based on filtering and low-rank decomposition according to claim 1, characterized in that, In step 1, the PCA algorithm is used to transform the original hyperspectral image H∈R M×N×B Down to p The method is as follows: Calculated from the covariance matrix p The largest eigenvalues ​​are { λ 1, λ 2, …, λ p }, and sort the obtained principal component components in descending order of information content, where p The principal components are represented as {H1, H2, ..., H...} p }; When performing subsequent spatial feature extraction operations, only for p The images are fused using a weighted average method based on their feature values. The mathematical expression for this method is: In the formula, Indicates the first i Each feature value.

3. The spatial-spectral joint hyperspectral image anomaly detection method based on filtering and low-rank decomposition according to claim 1, characterized in that, In step 2, an improved spatial filtering algorithm is used to extract R. pca Regarding the spatial information in the image, the improved spatial filtering algorithm makes the following improvements to the bilateral filter: Step 2.1, the first distance weight of the improved spatial filter is determined by the Euclidean distance between the reconstruction center and the neighboring pixels; Step 2.2, the second gradient weight of the improved spatial filter is determined by the gradient difference between pixel gray values; Step 2.3: The improved spatial filter combines the distance weights and gradient weights in the spatial domain. The kernel of the improved spatial filter is obtained by multiplying the distance kernel and the gradient kernel.

4. The spatial-spectral joint hyperspectral image anomaly detection method based on filtering and low-rank decomposition according to claim 3, characterized in that, Step 2.1, the weighting formula based on the Euclidean distance between pixels, is expressed as: In the formula, x The center pixel of the template window, i.e., the pixel to be reconstructed, has coordinates ( k , l ), ε The pixels surrounding the central pixel to be reconstructed have coordinates of... , d d This represents the Euclidean distance between two points. σ d The distance weight parameter indicates that the closer a pixel is to the pixel to be reconstructed, the more similar it is to the pixel to be reconstructed, and the greater its contribution when reconstructing the center pixel. In step 2.2, the improved gradient weight formula for the spatial filter is as follows: In the formula, σ r These are the gradient difference weight parameters. d r Indicates gradient difference, G ( x )and G ( ε )express x Point and ε gradient value at the point G ( x The mathematical formula for ) is: In the formula, X Represents the pixel to be reconstructed x The greater the difference in pixel gradient values ​​between two points, the smaller the weight coefficient. In step 2.3, the kernel representation of the improved spatial filter is as follows: The improved formula for calculating the spatial filter is as follows: In the formula, ε i The center pixel is the pixel to be reconstructed. x Other pixels in the convolution kernel n To improve the number of pixels in the spatial filter convolution kernel, k ( x ) represents the regularization parameter, expressed as: After spatial filtering calculation x The pixel values ​​of the location pixels are calculated, and after calculating all pixels, the initial spatial feature map R is obtained. spa .

5. The spatial-spectral joint hyperspectral image anomaly detection method based on filtering and low-rank decomposition according to claim 1, characterized in that, In step 3, the hyperspectral background reconstruction image R is obtained. bg The process is as follows: Step 3.1: Represent the hyperspectral data as a third-order tensor, i.e., H X ∈R M×N×B The rows and columns of the image are respectively M and N The number of bands in the image is B Then the original hyperspectral image H∈R M×N×B Tucker decomposition, represented as: In the formula, This represents the 1-mode tensor product of g and A. Represents the 2-mode tensor product. This represents the 3-mode tensor product, where the size of the factor matrices A, B, and C under different modes are M×P, N×U, and B×R, respectively. Let g represent the kernel tensor, g∈R P×U×R , g pur for The The nth element, the nth element of matrix A p The vectors are a p The first matrix B u The vectors are b u The first matrix C r The vectors are c r ; Step 3.2: The third-order tensor Tucker decomposition is solved using the iterative alternating least squares method, and the objective function of the Tucker decomposition is optimized as follows: In the formula, I is the identity matrix; Step 3.3: The Tucker decomposition is optimized using a higher-order singular value decomposition algorithm. Factor matrices A, B, and C of the solution matrix are solved in each mode, and the kernel tensor of the tensor projected in each mode is calculated. The result is { The method is as follows: Sort A, B, and C in descending order based on their eigenvalues, select the first few eigenvalues ​​for each eigenvalue, and substitute their corresponding eigenvectors into the eigenvectors. The kernel tensor was calculated. The reconstructed background image R is obtained through matrix multiplication. bg .

6. The spatial-spectral joint hyperspectral image anomaly detection method based on filtering and low-rank decomposition according to claim 1, characterized in that, Step 4 employs an improved method. k Mean clustering algorithm to separate and reconstruct the background image R bg As a sample, select representative pixels, that is... n Select from 1 sample k The most representative cluster centers are identified through the following steps: Step 4.1: Randomly select the first cluster center, and for the remaining cluster centers... n -1 samples, each sample x i The minimum spatial Euclidean distance from all current cluster centers is denoted as . D ( x i ), all samples D ( x i The maximum value in ) is D max ( x i ); Step 4.2, for the next cluster center, D ( x i Less than D max ( x i The selection probability of a sample with a probability of 1 / 2 is reduced to zero, and the remaining samples in the dataset with a non-zero probability are... D ( x i The sum of these two elements is denoted as . sum ( D ( x i )),but D ( x i ) / sum ( D ( x i That is, the sample x i The probability of being selected as the next cluster center; After calculating the probability of all points, the next cluster center is determined using the roulette wheel method; Step 4.3: Repeat the minimum distance calculation and cluster center determination operations until the first cluster center is selected. k Cluster centers, i.e. k The background dictionary D consists of 10 pixels.

7. The spatial-spectral joint hyperspectral image anomaly detection method based on filtering and low-rank decomposition according to claim 1, characterized in that, Step 5, the process of performing low-rank decomposition on the original hyperspectral image H, is as follows: Step 5.1, convert the original hyperspectral image H∈R M×N×B Substituting into the low-rank representation model, we get: In the formula, Z represents the coefficient matrix and E represents the sparse matrix; Step 5.2: Obtain the anomaly target detection result map by solving the sparse matrix E; use l 2,1 The norm constrains the sparse matrix E to maximize its sparsity; the objective function is expressed as: In the formula, rank( ) represents the rank function. λ The coefficients of the balancing terms are represented. represent l 2,1 Norm, defined as the norm of all columns of a matrix l The sum of the 2-norms; Step 5.3: Introduce the auxiliary variable W, and use the nuclear norm of matrix Z to approximate the rank of matrix Z, thus transforming the objective function into... In the formula Represents the calculation of the nuclear norm; Step 5.4: Construct a Lagrange function to solve the optimization problem of the objective function; Step 5.5: Solve for the Lagrange function; Step 5.6: After solving the low-rank representation model of the original image H, calculate the sparse matrix E. l The 2-norm can be used to obtain the spectral feature image R. spe ∈R M×N This represents the spectral feature value of each pixel in the image. Calculating all spectral feature values ​​yields the spectral feature image R. spe .

8. The spatial-spectral joint hyperspectral image anomaly detection method based on filtering and low-rank decomposition according to claim 7, characterized in that, In step 5.2, E l 2,1 Norm is represented as: In step 5.4, the Lagrange function is expressed as: In the formula, Y 1 , Y 2 They are Lagrange multipliers, This represents the penalty parameter.

9. The method for anomaly detection in space-spectrum joint hyperspectral images based on filtering and low-rank decomposition according to claim 7 or 8, characterized in that, Step 5.5 involves solving the Lagrange function, as follows: (1) Fix Z, E, Y 1 , Y 2 Update variable W, and solve using the following formula: In the formula, This represents the singular value thresholding operation; (2) Fix W, E, Y 1 , Y 2 Update variable Z, and solve using the following formula: In the formula, Indicates the number of iterations; (3) Fix W, Z, Y 1 , Y 2 Update variable E, and solve using the following formula: In the formula, S This indicates soft threshold contraction; (4) Fix W, Z, E, and update the Lagrange multipliers. Y 1 , Y 2 The solution formula is: 。 10. The spatial-spectral joint hyperspectral image anomaly detection method based on filtering and low-rank decomposition according to claim 1, characterized in that, In step 6, R spa and R spe Perform dot product fusion, that is: R result =R spa ·R spe In the formula, R result This is the final detection result image.

Images