A Hyperspectral Image Anomaly Detection Method Based on a New Multivariate Skewed t-Distribution Model

Abstract

The present invention relates to a hyperspectral image anomaly detection method based on a new multivariate skewed t-distribution model, and realizes the anomaly detection of hyperspectral images by using a method of modeling with an autoencoder and a new multivariate skewed t-distribution model. The specific method includes: First, use the RX algorithm to preprocess the original hyperspectral image, and select the background part of the result and input it into the stacked denoising autoencoder network for training; Then, input the original hyperspectral image into the trained stacked denoising autoencoder, and the output is the reconstructed hyperspectral image; Finally, subtract the reconstructed hyperspectral image from the original hyperspectral image to obtain the reconstruction error, and then use the multivariate skewed t-distribution model to model the reconstruction error, calculate the anomaly score to obtain the anomaly detection result.

Description

Technical Field

[0001] The present invention belongs to the field of remote sensing image processing technology, and particularly relates to a hyperspectral image anomaly detection method based on an autoencoder and a new multivariate skew t-distribution model. Background Technique

[0002] Hyperspectral images add a spectral dimension on the basis of traditional remote sensing images and are three-dimensional data structures containing the spatial features and spectral information of targets. Hyperspectral image detection techniques can generally be divided into target detection and anomaly detection. Due to the characteristic that abnormal targets can be detected without prior knowledge of the spectral information of the target object, hyperspectral anomaly detection is widely used in many fields such as agriculture, meteorology, and military. The research on hyperspectral anomaly detection has become increasingly important.

[0003] Hyperspectral image anomaly detection based on autoencoders mainly includes two steps: image reconstruction and anomaly detection. Image reconstruction rebuilds the image by finding the background signal in the hyperspectral image, and anomaly detection evaluates whether each pixel belongs to the background based on an anomaly scoring mechanism. Pixels that do not belong to the background are abnormal pixels. The quality of the anomaly scoring mechanism directly affects the final anomaly detection result. At the same time, since most of the current anomaly detection algorithms based on autoencoders are deviation-based methods, the reconstruction error is considered to follow the Laplace distribution and normal distribution in the form of L1 and L2 norms, and the effect is not good in actual application and lacks robustness. Summary of the Invention

[0004] To solve the defects existing in the above-mentioned prior art, the purpose of the present invention is to provide a hyperspectral image anomaly detection method based on a new multivariate skew t-distribution model, which realizes hyperspectral image anomaly detection based on an autoencoder and a new multivariate skew t-distribution model.

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

[0006] A hyperspectral image-based anomaly detection method includes the following steps:

[0007] Step 1, for the original hyperspectral image data X ∈ R MN×B , use the RX algorithm for preprocessing to obtain its background image part; where M and N represent the spatial dimensions of the image, MN is the number of pixels, and B represents the spectral dimension;

[0008] Step 2, use the background image part as the input of the stack denoising autoencoder network based on the spectral angle cosine, and train the network to obtain a trained stack denoising autoencoder;

[0009] Step 3: Take the original hyperspectral image as the input of the trained stacked denoising autoencoder, select the output of the output layer, and obtain the reconstructed hyperspectral image data \(X'\in R\) MN×B ;

[0010] Step 4: Calculate the difference \(d\) i ' between the reconstructed spectral vector \(x\) i ' of the \(i\)-th pixel in the reconstructed hyperspectral image and the original spectral vector \(x\) i , so as to obtain a spectral vector difference image, and use a local standard deviation filter for the spectral vector difference image, so as to obtain a more prominent reconstruction error \(r\) i \(\in R\) B , \(i = 1, 2, \ldots, MN\);

[0011] Step 5: Input the reconstruction error into the MVSkt distribution model, calculate the anomaly score \(AD(r\) i ), and judge the category of the pixel to be measured according to the anomaly score \(AD(r\) i ), and output the anomaly detection result.

[0012] In the said Step 1, the acquisition process of the background image part is as follows:

[0013] Step 1.1: Let the spectral vector of the \(i\)-th pixel in the original hyperspectral image be \(x\) i \(\in R\) B ;

[0014] Step 1.2: Assume that the background pixels in the original hyperspectral image follow a Gaussian distribution. By constructing a target and background window, calculate the difference value between the pixel to be measured and the background pixel in the local area of the image using the Mahalanobis distance, and judge the type of the pixel to be measured according to the size of the difference value. The judgment result is calculated according to the following formula:

[0015]

[0016] In the formula, \(\mu\) b represents the mean of the background, \(C\) b represents the covariance matrix of the background, \(\eta\) represents the RX algorithm decision threshold. If \(RX(x\) i ) \(\lt \eta\), it is initially judged that the pixel to be measured is the background part, otherwise, it is initially judged that the pixel to be measured is the abnormal part;

[0017] Step 1.3: For large abnormal target images containing more than 1% of the total number of pixels in the image, set the RX algorithm threshold to 500 and obtain the background image using the fixed-value replacement method; for small abnormal target images containing less than 1% of the total number of pixels in the image, set the RX algorithm threshold to 350 and obtain the background image using the surrounding pixel replacement method; finally, divide the original hyperspectral image into a background image part and an abnormal image part according to the detection results of the RX algorithm.

[0018] In the said Step 2, the training process of the stacked denoising autoencoder is as follows:

[0019] Step 2.1: Use the stacked denoising autoencoder network and, during training, use a reconstruction loss function based on the spectral angle cosine cos(θ SA ), where the calculation formula of cos(θ SA ) is as follows:

[0020]

[0021] In the formula, A and B represent the output reconstructed spectral vectors and the input spectral vectors of the autoencoder network, and O represents the number of spectral bands;

[0022] Incorporate the spectral angle into the reconstruction loss function E SA to obtain the following formula:

[0023]

[0024] In the formula, f(z (L) ) represents the output reconstructed spectral vector A of the autoencoder, y represents the input spectral vector B, L represents the number of layers of the autoencoder, and z (L) in f(z (L) ) is expressed as:

[0025] z (l) = W (l-1) a (l-1) + b (l-1)

[0026] In the formula, a (l) = f(z (l) ) represents the spectral vector after being processed by the l-th layer of the autoencoder, a (1) = x is the original spectral vector, l = L, L - 1, L - 2, …, 2, L represents the number of layers of the autoencoder, and W (l-1) , b (l-1) represent the weight matrix and the bias matrix of the (l - 1)-th layer of the autoencoder respectively;

[0027] Incorporate the regularization term to obtain the final reconstruction loss function E SA as follows:

[0028]

[0029] In the formula, W and b represent the weight matrix and the bias matrix respectively, and y (m) represents the input spectral vector of the m-th pixel, λ is the regularization parameter, and I and J represent the number of pixels in the l-th layer and the l+1-th layer respectively;

[0030] Step 2.2: Use the Sigmoid function as the activation function in the stacked denoising autoencoder network. In its hidden layer, it is set that the number of neurons in hidden layers 1 to 5 is symmetrically distributed based on the number of neurons in hidden layer 3, and the number of neurons in the first to third autoencoders decreases layer by layer;

[0031] Step 2.3: After inputting the background image, first perform pre-training. In the first step of training, the parameters from the input layer to hidden layer 1 and from hidden layer 5 to the output layer are trained to realize the reconstruction of the image with Gaussian noise added to the input layer. After the training is completed, the obtained parameters are fixed; in the second step of training, the parameters from hidden layer 1 to hidden layer 2 and from hidden layer 4 to hidden layer 5 are trained to realize the reconstruction of the output image of hidden layer 1. After the training is completed, the obtained parameters are fixed; in the third step of training, the input and output parameters of the middlemost hidden layer 3 are trained to realize the reconstruction of the output image of hidden layer 2. After the training is completed, the obtained parameters are fixed;

[0032] Step 2.4: After pre-training, the values of each parameter in the network are initialized to appropriate values. Then, the hidden layers of the trained parameters are integrated together for training to realize the reconstruction of the image with Gaussian noise added to the input layer, and the parameters obtained from the previous training of each layer are adjusted within a small range. Finally, a trained stacked denoising autoencoder is obtained.

[0033] In step 2.2, the network learning rate is set to 10 -3 , the epoch is set to 100, and the batch size is set to 400.

[0034] In step 4, the reconstruction error r i is obtained as follows:

[0035] Step 4.1: For the reconstructed hyperspectral image X′∈R MN×B , calculate the difference d i between the reconstructed spectral vector x i ′ of each pixel and the original spectral vector x i according to the following formula, so as to obtain the spectral vector difference image;

[0036] d i =x i -x i ′

[0037] Step 4.2, preprocess the obtained spectral vector difference image using a local standard deviation filter. After filtering by the filter, the local standard deviation of the abnormal pixels is higher than that of the background pixels, so that the reconstruction error r with more prominent anomalies is obtained through filtering. i 。

[0038] In the said step 5, the calculation process of the anomaly score AD(r i ) is as follows:

[0039] Step 5.1, use the Multivariate Skewed t-distribution (MVSkt) model to simulate the reconstruction error r i , in this model, r i is expressed as the following formula:

[0040]

[0041] In the formula, w i represents a Gaussian random vector, z represents a random variable of the generalized inverse Gaussian distribution, m is the mean vector, b is the bias vector, T is the inverse covariance matrix, it is set that b > 1, α < 0, γ → +∞, and m and T satisfy the normal-Wishart prior distribution, as shown in the following formula:

[0042]

[0043] In the formula, m0, λ, v0 and Ψ0 are hyperparameters;

[0044] Then the probability distribution function of the multivariate skewed t-distribution model is as shown in the following formula:

[0045]

[0046] In the formula R i =(r i -m) T T(r i -m), K α (·) is the modified Bessel function of the second kind of order α, is the modified Bessel function of the second kind order;

[0047] Step 5.2, when the parameters are determined, calculate the negative logarithm of the probability distribution function as the anomaly score to achieve anomaly detection, then its anomaly score is as shown in the following formula,

[0048]

[0049] Step 5.3, input the reconstruction error r obtained in Step 4 i into the anomaly scoring formula for calculation, and the final anomaly score AD(r i ) is obtained;

[0050] Step 5.4, use the threshold segmentation detection method to judge the category of the pixel to be measured according to the anomaly score AD(r i ), and obtain the final anomaly detection result.

[0051] Compared with the prior art, the present invention has lower requirements for parameters, better adaptability of the algorithm, and realizes good anomaly detection on this basis, and has good robustness. Description of the Drawings

[0052] Figure 1 is the pseudo-color map of the original hyperspectral image of the present invention.

[0053] Figure 2 is the flow chart of the present invention.

[0054] Figure 3 is the partial diagram of the background image of the RX algorithm preprocessing result.

[0055] Figure 4 is the structural diagram of the stacked denoising autoencoder network.

[0056] Figure 5 is the AUC result diagram of the stacked denoising autoencoder training process.

[0057] Figure 6 is the comparison diagram of the anomaly detection results in the example of the present invention, where (a) is the detection result of the MVC method, (b) is the detection result of the MVL method, and (c) is the detection result of the MVSkt method.

[0058] Figure 7 is the three-dimensional curve comparison diagram of the anomaly detection results in the example of the present invention, where (a) is the detection result of the MVC method, (b) is the detection result of the MVL method, and (c) is the detection result of the MVSkt method.

[0059] Figure 8 is the ROC diagram of the comparison of the anomaly detection results in the example of the present invention.

[0060] Table 1 is the AUC diagram of the comparison of the anomaly detection results in the example of the present invention. Detailed Embodiments

[0061] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present 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 only used to explain the present invention and are not used to limit the present invention.

[0062] As mentioned above, in the existing hyperspectral image anomaly detection methods, due to reasons such as imperfect scoring mechanisms and inaccurate reconstruction errors, the detection effect is poor, especially the lack of robustness. Therefore, establishing a new scoring mechanism and ensuring high accuracy and high robustness is a key point in the field of hyperspectral anomaly detection. The present invention is a series of improvements made to the above problems.

[0063] To demonstrate the effectiveness of the method of the present invention, real hyperspectral images are used for anomaly detection. The experiment uses the Abu-airport image, as Figure 1 shown. The data is from the Airport-Beach-Urban (ABU) hyperspectral dataset collected by an airborne visible / infrared imaging spectrometer (AVIRIS). The spatial resolution of the image is 17.2 m per pixel, and the spectral bands range from 365 to 2500 nm. The size of the image used in the experiment is 100×100, the number of bands is 205, and two airplanes in the image are used as anomaly targets. Referring to Figure 2 , the basic process of anomaly detection for the Abu-airport image data is as follows:

[0064] Step 1, for the original hyperspectral image X∈R MN×B use the RX algorithm for preprocessing to obtain the hyperspectral image belonging to the background part, that is, its background image part. Where M and N represent the spatial dimensions of the image, MN is the number of pixels, and B represents the spectral dimension. The specific steps can be described as follows:

[0065] Step 1.1, the original hyperspectral image is a three-dimensional matrix of 100×100×205, with a total of 100×100 spatial pixel points, containing the spectral data of the ground target within 205 bands. The spectral information corresponding to each spatial pixel is represented by a one-dimensional vector of size 1×1×205. Let the spectral vector of the i-th pixel in the original hyperspectral image be x i ∈R B .

[0066] Step 1.2, assuming that the background pixels in the original hyperspectral image follow a Gaussian distribution, by constructing a target and background window, calculate the difference value between the pixel to be measured and the background pixel using the Mahalanobis distance in the local area of the image, and judge the type of the pixel to be measured according to the size of the difference value. The judgment result is calculated according to the following formula.

[0067]

[0068] In the formula, μ b represents the mean of the background, C b represents the covariance matrix of the background, η represents the RX algorithm determination threshold. The meaning of this formula is: if RX(x i) If η is less than a certain value, it is preliminarily determined that the pixel to be measured is a background part; otherwise, it is preliminarily determined that the pixel to be measured is an abnormal part.

[0069] Step 1.3: Perform anomaly detection on the aircraft image using the RX algorithm.

[0070] The anomaly detection situation of the RX algorithm is crucial for the subsequent training of the autoencoder network. Since the RX algorithm determines anomalies based on a set threshold, it is necessary to conduct experimental research on the selection of the threshold. For easy comparison and observation, the detected anomaly pixel values are replaced with a fixed value of 0, which is represented as pure black pixels in the image. In the present invention, an anomaly target image containing more than 1% of the total number of pixels in the image is called a large anomaly target image, and an anomaly target image containing less than 1% of the total number of pixels in the image is called a small anomaly target image. For the large anomaly target image, the RX algorithm threshold is set to 500, and the background image is obtained using the fixed value replacement method; for the small anomaly target image, the RX algorithm threshold is set to 350, and the background image is obtained using the surrounding pixel replacement method. Exemplarily, in this embodiment, the aircraft target is relatively large and obviously belongs to a large target. Therefore, the RX algorithm threshold is set to 500, and the background image is obtained using the fixed value replacement method. Finally, the original hyperspectral image is divided into a background image part and an anomaly image part according to the detection results of the RX algorithm, a total of two parts. The processed background image part is as Figure 3 shown.

[0071] Step 2: Use the background image part as the input of the stacked denoising autoencoder network based on the spectral angle cosine, and train the network to obtain a trained stacked denoising autoencoder. The specific training process is as follows:

[0072] Step 2.1: The present invention uses a stacked denoising autoencoder network, which combines a stacked autoencoder with five hidden layers and a denoising autoencoder that adds random Gaussian noise. The specific structure of the network is as Figure 4 shown.

[0073] During training, a reconstruction loss function based on the spectral angle cosine cos(θ SA ) is used, where the calculation formula of cos(θ SA ) is as follows:

[0074]

[0075] In the formula, A and B represent the output reconstruction spectral vector and the input spectral vector of the autoencoder network, and O represents the number of spectral bands;

[0076] Incorporating the spectral angle into the reconstruction loss function E SA results in the following formula:

[0077]

[0078] wherein, f(z (L) ) represents the output reconstructed spectral vector A of the autoencoder, y represents the input spectral vector B, L represents the number of layers of the autoencoder, and z in f(z (L) ) is expressed as: (L) is expressed as:

[0079] z (l) = W (l-1) a (l-1) + b (l-1)

[0080] wherein, a (l) = f(z (l) ) represents the spectral vector after being processed by the l-th layer of the autoencoder, a (1) = x is the original spectral vector, l = L, L - 1, L - 2, …, 2, L represents the number of layers of the autoencoder, and W (l-1) , b (l-1) respectively represent the weight matrix and the bias matrix of the (l - 1)-th layer of the autoencoder;

[0081] Including the regularization term, the final reconstruction loss function E SA is as follows:[[]]

[0082]

[0083] wherein, W and b respectively represent the weight matrix and the bias matrix, y (m) represents the input spectral vector of the m-th pixel, λ is the regularization parameter, and I and J respectively represent the number of pixels in the l-th layer and the (l + 1)-th layer.

[0084] Step 2.2, the Sigmoid function is used as the activation function in the stacked denoising autoencoder network, and the network learning rate is set to 10 -3, the epoch is set to 100 and the batch size is set to 400. In its hidden layer, it is set that hidden layers 1 to 5 are symmetrically distributed based on the number of neurons in hidden layer 3. Considering the purpose of network data feature extraction, the number of neurons in the first to the third autoencoders decreases layer by layer. And in order to avoid too many variables, in this embodiment, the number of encoding units in the first and second encoders is set to 100 and 50 respectively, that is, the number of neurons in hidden layer 1 and hidden layer 5 is 100, the number of neurons in hidden layer 2 and hidden layer 4 is 50, and for the number of neurons in the middlemost hidden layer 3, since it is directly related to the spectral dimension of the final output image, a total of seven values of 5, 10, 15, 20, 25, 30, and 35 are set as the number of neurons in hidden layer 3 for experiments. Three autoencoder networks, MSE - SDAE, CSA - SDAE, and SID - SDAE, are used to train the Pavia bridge image and the Urban road image respectively, and combined with the collaborative representation algorithm for anomaly detection. Judging according to the AUC result, 30 is selected as the number of neurons in hidden layer 3, and the AUC result is as Figure 5 shown.

[0085] Step 2.3, after inputting the background image, first perform pre - training. The first - step training is to train the parameters from the input layer to hidden layer 1 and from hidden layer 5 to the output layer to realize the reconstruction of the image with Gaussian noise added to the input layer. After the training is completed, fix the parameters obtained from the training. The second - step training is to train the parameters from hidden layer 1 to hidden layer 2 and from hidden layer 4 to hidden layer 5 to realize the reconstruction of the output image of hidden layer 1. After the training is completed, fix the parameters obtained from the training. The third - step training is to train the input and output parameters of the middlemost hidden layer 3 to realize the reconstruction of the output image of hidden layer 2. After the training is completed, fix the parameters obtained from the training.

[0086] Step 2.4, after the pre - training of the previous three steps, the parameter values in the network are all initialized to appropriate values. Then integrate the hidden layers of the parameters obtained from the previous three - step training together for training to realize the reconstruction of the image with Gaussian noise added to the input layer, and make a small - range adjustment to the parameters obtained from the previous training of each layer. Finally, a trained stacked denoising autoencoder is obtained.

[0087] Step 3, for the trained stacked denoising autoencoder, input the original hyperspectral image for image reconstruction. The specific steps include: taking the original hyperspectral image data X ∈ R MN×B as the input and passing it into the stacked denoising autoencoder based on the spectral angle cosine, and selecting the output result of the output layer to obtain the reconstructed hyperspectral image data X′ ∈ R MN×B .

[0088] Step 4, calculate the reconstructed spectral vector x i ′ of the i - th pixel in the reconstructed hyperspectral image and the original spectral vector xi Difference d i , thus obtaining a spectral vector difference image, and then using a local standard deviation filter on the spectral vector difference image, so as to obtain a more prominent reconstruction error r through filtering i ∈R B (i = 1, 2, …, MN). The specific steps can be described as follows:

[0089] Step 4.1, for the reconstructed hyperspectral image X′ ∈ R output by the autoencoder MN×B , calculate the reconstructed spectral vector x of each pixel according to the following formula i ′ and the original spectral vector x i Difference d i , thus obtaining a spectral vector difference image.

[0090] d i = x i - x i ′

[0091] Step 4.2, preprocess the obtained spectral vector difference image using a local standard deviation filter. After filtering by the filter, the local standard deviation of abnormal pixels is higher than that of background pixels, so as to obtain a more prominent reconstruction error r through filtering i ∈R B (i = 1, 2, …, MN).

[0092] Step 5, input the reconstruction error r i into the MVSkt distribution model, calculate the anomaly score AD(r i ) at different thresholds, judge the category of the pixel to be measured according to the anomaly score AD(r i ), and output the anomaly detection result. And draw the corresponding ROC curve and AUC surface plot according to the false alarm probability and detection probability in the detection result, and analyze the performance of the MVSkt distribution model by comparing with other detection models. The specific steps can be described as follows:

[0093] Step 5.1, the present invention uses a Multivariate Skewed t-distribution (MVSkt) model to simulate the reconstruction error r i , in this model r i can be expressed as the following formula

[0094]

[0095] In the formula, w i represents a Gaussian random vector z represents a random variable of the generalized inverse Gaussian distribution Let \(m\) be the mean vector, \(b\) be the bias vector, and \(T\) be the inverse covariance matrix. Assume \(b > 1\), \(\alpha < 0\), \(\gamma\rightarrow+\infty\), and \(m\) and \(T\) satisfy the normal-Wishart prior distribution, as shown in the following equation.

[0096]

[0097] where \(m_0\), \(\lambda\), \(\nu_0\), and \(\Psi_0\) are hyperparameters.

[0098] Then the probability density function of the multivariate skew t-distribution model is shown in the following equation.

[0099]

[0100] where \(R\) i \(=(r\) i - m)^T(r\) T - m)\), \(K\) i \((\cdot)\) is the modified Bessel function of the second kind of order \(\alpha\), α and \(K\) is the modified Bessel function of the second kind of order \(\alpha\).

[0101] Step 5.2, calculate each parameter in the model. Using the VB method, first, we need to consider the maximum marginal likelihood estimation problem, as shown in the following equation.

[0102]

[0103] where \(\Theta\) represents the set of parameters of the MVSkt distribution model. Integrate over the latent variables \(\Phi=\{m, T, z\}\), and use the approximate posterior distribution \(q(\Phi)\). The marginal log-likelihood in the equation can be expressed as

[0104] \(\log f(r\) 1:N |\Theta)=F(q,\Theta)+D\) KL \((f(\Phi|r\) 1:N ,\Theta)||q)\)

[0105] where the solutions of \(F(q,\Theta)\) and \(D\) KL \((f||q)\) are shown in the following equations respectively.

[0106] \(F(q,\Theta)=\int\log((f(r\) 1:N ,\Phi|\Theta)) / q(\Phi))q(\Phi)d\Phi\)

[0107] \(D\) KL \((f||q)=-\int\log((f(\Phi r\) 1:N ,\Theta)) / q(\Phi))q(\Phi)d\Phi\)

[0108] where \(D\) KL \((f||q)\) represents \(f(\Phi|r\)1:N , Θ) and the KL divergence between q(Φ).

[0109] The VB method uses the factorization of the exact posterior value to approximately represent the posterior distribution q(Φ)), as shown in the following equation

[0110] f(Φ∣r 1:N , Θ) ≈ q(Φ) = q(m)q(T)q(z)

[0111] This approximation method facilitates the calculation of the expected integral. Therefore, in the VB method, it is necessary to first find q(m), q(T), and q(z), and then update the parameter set iteratively.

[0112] For q(m), its distribution should satisfy the following equation

[0113]

[0114] where const. represents a constant.

[0115] The obtained q(m) should follow a multivariate normal distribution, that is where

[0116]

[0117] C m = ((E[z -1 N + λ)E[T]) -1

[0118]

[0119] For q(T), its distribution should satisfy the following equation

[0120]

[0121] The obtained q(T) should follow a Wishart distribution, that is where

[0122]

[0123] v N = (τ + N + 1)

[0124] For q(z), its distribution should satisfy the following equation

[0125]

[0126] The obtained q(z) should follow a generalized inverse Gaussian distribution, that is where

[0127] c = Ntr(E(T)bbT )

[0128]

[0129] The calculation formula for the expected value is as follows:

[0130] E[m]=μ m

[0131]

[0132]

[0133]

[0134] After the initial parameters are set, the next step is to iteratively update the parameters. By maximizing the marginal log-likelihood representation, that is, by maximizing the lower bound of F(q,Θ), the parameters are found. The lower bounds of m0, Ψ0, β, and λ are given by the following equations respectively.

[0135]

[0136]

[0137]

[0138]

[0139] Calculate the maximized estimator, where the variables and parameters are

[0140] m0 = E[m]

[0141] Ψ0=(τ / (τ - L - 1))(E[T]) -1

[0142] β=-(2α / (E[z -1 ))

[0143]

[0144] Some of the parameters are as follows:

[0145] α = L / 2

[0146] τ = (((N - L - 1)p + L + 1) / (1 - p))

[0147] p = ((τ - L - 1) / (N + τ - L - 1))

[0148] 0 ≤ p ≤ 1

[0149] Step 5.3, when the parameters are determined, calculate the negative logarithm of the probability distribution function as the anomaly score to achieve anomaly detection. Then, the anomaly score AD is shown as follows:

[0150]

[0151] Step 5.4, input the reconstruction error r i ∈R B (i = 1, 2, …, MN) obtained in Step 4 into the anomaly score formula for calculation, that is, the final anomaly score AD(r i ) is obtained. Use the threshold segmentation detection method to judge the category of the pixel to be measured according to the anomaly score AD(r i ) to obtain the final anomaly detection result.

[0152] Step 5.5, draw the false color map of the corresponding final anomaly score AD(r i ), as shown in Figure 6 . Draw the three-dimensional curve map of the corresponding final anomaly score AD(r i ), as shown in Figure 7 . Compared with other templates, it highlights the anomaly more accurately.

[0153] Step 5.6, estimate the false alarm rate and accuracy according to the template and draw the ROC curve graph and AUC area graph. The ROC curve graph is as shown in Figure 8 . The AUC statistical results are compared with other models on the basis of considering the running time, and the results are shown in the following table.

[0154] Algorithm Name AUC Time (s) GRX 0.8403 0.0693 LRX 0.9052 41.1269 CR 0.9069 6.6130 MVC 0.9115 6.9031 MVJ 0.9105 4.7637 MVL 0.9101 4.7676 MVSt 0.9105 4.8084 MVN 0.8458 0.1245 MVSkt in This Chapter 0.9353 6.7372

[0155] According to the comparison results, it can be seen that the AUC value of the MVSkt algorithm proposed by the present invention for anomaly detection on the Abu-airport image is significantly higher than that of other algorithms. The AUC values of the GRX and MVN algorithms are the lowest. In addition, the effects of several other anomaly detection algorithms based on autoencoders are better than those of traditional algorithms. It can be seen that for anomaly detection in complex scenes, the detection effect of the normal mixture mean-variance distribution method based on autoencoders is more applicable than that of traditional methods. Compared with several other multivariate normal mixture mean-variance distributions, the MVSkt algorithm further improves the detection performance of the autoencoder method. In terms of the running time of the algorithm, the time-consuming of the algorithm proposed by the present invention is at the intermediate level and is relatively close to the time-consuming of several multivariate distribution model algorithms. Compared with the traditional algorithms used, the multivariate distribution algorithm based on autoencoders has lower requirements for parameters and better adaptability of the algorithm.

Claims

1. A hyperspectral image anomaly detection method based on a new type of multivariate skewed t-distribution model, characterized in that, Including the following steps: Step 1. For the original hyperspectral image data \(X\in R\) MN×B , preprocess it using the RX algorithm to obtain its background image part; where \(M\) and \(N\) represent the spatial dimensions of the image, \(MN\) is the number of pixels, and \(B\) represents the spectral dimension. Step 2: Use the background image part as the input of the stacked denoising autoencoder network based on the spectral angle cosine, and train the network to obtain a trained stacked denoising autoencoder; Step 3: Use the original hyperspectral image as the input of the trained stacked denoising autoencoder, and select the output of the output layer to obtain the reconstructed hyperspectral image data X′∈R MN×B ; Step 4, calculate the reconstructed spectral vector x of the i-th pixel in the reconstructed hyperspectral image i ′ and the original spectral vector x i to obtain the difference d i of the spectral vectors, thereby obtaining a spectral vector difference image, and applying a local standard deviation filter to the spectral vector difference image, so as to obtain a more prominent reconstruction error r through filtering i ∈R B , where i = 1, 2, …, MN; Step 5: Input the reconstruction error into the MVSkt distribution model to calculate the anomaly score AD(r i ), and determine the class of the pixel to be measured according to the anomaly score AD(r i ), and output the anomaly detection result.

2. The hyperspectral image anomaly detection method based on the novel multivariate skew t-distribution model according to claim 1, characterized in that, In the above Step 1, the process of obtaining the background image part is as follows: Step 1.1, let the spectral vector of the \(i\)-th pixel in the original hyperspectral image be \(\mathbf{x}\) i \(\in\mathbb{R}\) B ; Step 1.2: Assume that the background pixels in the original hyperspectral image follow a Gaussian distribution. By constructing the target and background windows, calculate the difference value between the pixel to be measured and the background pixels in the local area of the image using the Mahalanobis distance, and judge the type of the pixel to be measured according to the size of the difference value. The judgment result is calculated according to the following formula: where μ b represents the mean of the background, C b represents the covariance matrix of the background, η represents the decision threshold of the RX algorithm. If RX(x i ) < η, it is preliminarily determined that the pixel to be measured is the background part; otherwise, it is preliminarily determined that the pixel to be measured is the abnormal part; Step 1.3: For a large abnormal target image containing more than 1% of the total number of pixels in the image, set the RX algorithm threshold to 500, and use the fixed-value replacement method to obtain the background image; for a small abnormal target image containing less than 1% of the total number of pixels in the image, set the RX algorithm threshold to 350, and use the surrounding pixel replacement method to obtain the background image; finally, divide the original hyperspectral image into the background image part and the abnormal image part according to the detection result of the RX algorithm.

3. The hyperspectral image anomaly detection method based on the novel multivariate skewed t-distribution model according to claim 1, wherein, In the above Step 2, the training process of the stacked denoising autoencoder is as follows: Step 2.1: Use a stacked denoising autoencoder network. During training, use a reconstruction loss function based on the spectral angle cosine cos(θ SA ), where the calculation formula of cos(θ SA ) is as follows: In the formula, A and B represent the output reconstructed spectral vector and the input spectral vector of the autoencoder network, and O represents the number of spectral bands; Incorporate the spectral angle into the reconstruction loss function E SA There is the following formula: where, f(z (L) ) represents the output reconstructed spectral vector A of the autoencoder, y represents the input spectral vector B, L represents the number of layers of the autoencoder, and z (L) in f(z (L) is expressed as: z (l) = W (l-1) a (l-1) + b (l-1) where a (l) = f(z (l) ) represents the spectral vector after being processed by the l-th layer of the autoencoder, a (1) = x is the original spectral vector, l = L, L - 1, L - 2, …, 2, L represents the number of layers of the autoencoder, W (l-1) , b (l-1) respectively represent the weight matrix and the bias matrix of the (l - 1)-th layer of the autoencoder; Including the regularization term, the final reconstruction loss function E is obtained SA As shown in the following equation where W and b represent the weight matrix and the bias matrix respectively, and y (m) represents the input spectral vector of the m-th pixel, λ is the regularization parameter, and I and J represent the number of pixels in the l-th layer and the l+1-th layer respectively; Step 2.2: Use the Sigmoid function as the activation function in the stacked denoising autoencoder network. In its hidden layer, set the number of neurons in hidden layers 1 to 5 to be symmetrically distributed based on the number of neurons in hidden layer 3, and the number of neurons in the first to third autoencoders decreases layer by layer; Step 2.3: After inputting the background image, first perform pre-training. In the first step of training, train the parameters from the input layer to hidden layer 1 and from hidden layer 5 to the output layer to realize the reconstruction of the image with Gaussian noise added to the input layer. After the training is completed, fix the obtained parameters; in the second step of training, train the parameters from hidden layer 1 to hidden layer 2 and from hidden layer 4 to hidden layer 5 to realize the reconstruction of the output image of hidden layer 1. After the training is completed, fix the obtained parameters; in the third step of training, train the input and output parameters of the middle hidden layer 3 to realize the reconstruction of the output image of hidden layer 2. After the training is completed, fix the obtained parameters; Step 2.4: After pre-training, the parameter values in the network are all initialized to appropriate values. Then, integrate the hidden layers of the obtained parameters for training to realize the reconstruction of the image with Gaussian noise added to the input layer, and make a small adjustment to the parameters obtained from the previous training of each layer. Finally, obtain a trained stacked denoising autoencoder.

4. The hyperspectral image anomaly detection method based on the novel multivariate skewed t-distribution model according to claim 3, wherein, In step 2.2, the network learning rate is set to 10 -3 , the epoch is set to 100, and the batch size is set to 400.

5. The hyperspectral image anomaly detection method based on the new multivariate skewed t-distribution model according to claim 1, wherein In the said step 4, the reconstruction error r i is obtained as follows: Step 4.1, for the reconstructed hyperspectral image X′∈R MN×B , calculate the difference d i ′ between the reconstructed spectral vector x i of each pixel and the original spectral vector x i according to the following formula, so as to obtain the spectral vector difference image; d i = x i - x i ′ Step 4.2, preprocess the obtained spectral vector difference image using a local standard deviation filter. After filtering by the filter, the local standard deviation of abnormal pixels is higher than that of background pixels, so that a more prominent reconstruction error r is obtained through filtering i .

6. The hyperspectral image anomaly detection method based on the novel multivariate skewed t-distribution model according to claim 1, wherein The abnormal score AD(r i ) is calculated as follows: Step 5.1, use the Multivariate Skewed t-distribution (MVSkt) model to simulate the reconstruction error r i , in this model, r i is expressed as the following formula: where \(w\) i denotes a Gaussian random vector, \(z\) denotes a random variable of the generalized inverse Gaussian distribution, \(m\) is the mean vector, \(b\) is the bias vector, \(T\) is the inverse covariance matrix. It is set that \(b > 1\), \(\alpha < 0\), \(\gamma\rightarrow+\infty\), and \(m\) and \(T\) satisfy the normal-Wishart prior distribution as shown in the following formula: In the formula, m0, λ, v0, and Ψ0 are hyperparameters; Then the probability distribution function of the multivariate skew t-distribution model is shown as the following formula: where R i =(r i -m) T T(r i -m), K α (·) is the modified Bessel function of the second kind of order α, is the modified Bessel function of the second kind of order Step 5.2: When the parameters are determined, calculate the negative logarithm of the probability distribution function as the anomaly score to realize anomaly detection. Then its anomaly score is shown as the following formula, Step 5.3, input the reconstruction error r obtained in Step 4 i into the anomaly scoring formula for calculation, and the final anomaly score AD(r i ) can be obtained; Step 5.4, using the threshold segmentation detection method, determine the category of the pixel to be measured according to the anomaly score AD(r i ), and obtain the final anomaly detection result.

Images