Anomaly Detection Method for Spaceborne Hyperspectral Images Based on Joint Low-Rank Tensor Approximation

Abstract

The present invention discloses a spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation. The hyperspectral image is represented by a third-order tensor, and the original hyperspectral image tensor is decomposed into the sum of a background tensor and an anomaly tensor. According to the low-rank characteristics of the hyperspectral image in the spatial dimension and the spectral dimension, a joint low-rank constraint regular term and a spatial local smoothness regular term are designed for the background tensor in the low-rank decomposition model, and a per-tube sparse constraint is imposed on the anomaly tensor. Then, the alternating direction multiplier method is used to iteratively solve the above optimization model to obtain the optimal decomposition result. Finally, the final anomaly detection result is obtained according to the optimal anomaly tensor obtained by the decomposition. The method of the present invention can effectively achieve high-precision anomaly detection in the complex background of spaceborne hyperspectral images, and can be effectively applied to anomaly detection in different complex backgrounds such as oceans, land, forests, and gobi, laying an important foundation for subsequent accurate target detection and recognition.

Description

Technical Field

[0001] The present invention belongs to the field of hyperspectral image detection, and in particular relates to a spaceborne hyperspectral image anomaly detection method based on a joint low-rank near-tensor-like. Background Art

[0002] The so-called anomaly detection refers to detecting abnormal targets and suppressing background information without using prior spectral information of the targets, and the detection results can provide regions of interest for subsequent accurate target detection and recognition. Hyperspectral imaging technology simultaneously collects the geometric, radiation, and spectral information of the targets to form an image cube. When imaging the targets, hundreds of continuous and fine radiation intensity data of different wavelengths of the targets are obtained within a wide spectral range, forming the spectral characteristic curve "fingerprint" of the targets, enabling hyperspectral images to distinguish the subtle differences between different objects, thereby realizing the effective recognition of the targets. This unique advantage of hyperspectral images has enabled them to be widely applied in the fields of anomaly detection, target detection, and target recognition.

[0003] In the prior art, the implementation of various hyperspectral image anomaly detection methods depends on two basic characteristics of abnormal targets in hyperspectral images: 1) The abnormal targets are different from their surrounding backgrounds in spectral and spatial characteristics; 2) The occurrence probability of abnormal targets in hyperspectral images is low, occupying fewer pixels.

[0004] The commonly used anomaly detection methods can be roughly divided into five types: methods based on statistical theory, methods based on representation, methods based on decomposition, methods based on joint spatial-spectral, and methods based on deep learning.

[0005] The methods based on statistical theory rely on the construction of statistical models, and the actual images are ever-changing, making it difficult to describe all images with a definite model. The methods based on representation rely on the construction of an over-complete background dictionary, and it is often difficult to construct a dictionary that contains all types of backgrounds but does not contain abnormal targets in practical applications. The methods based on joint spatial-spectral rely more on spatial characteristics and are difficult to detect some targets with unclear spatial characteristics well. The methods based on deep learning require a large amount of data for training, and it is difficult to obtain a large amount of training data for hyperspectral anomaly detection tasks. At the same time, the real-time performance of such methods is difficult to guarantee.

[0006] At present, decomposition-based hyperspectral anomaly detection methods have achieved good detection performance. These methods decompose the original hyperspectral image into a background part and an anomaly part, and constrain the decomposition model according to the inherent characteristics of the background and the anomaly, so as to obtain the optimal decomposition result and the anomaly detection result. However, as mentioned in the literature "Prior-Based Tensor Approximation for Anomaly Detection in Hyperspectral Imagery", in many decomposition-based anomaly detection methods, the hyperspectral image has to be matrixized (rearranging the three-dimensional data into a two-dimensional matrix with spectral vectors as units), which destroys the inherent spatial-spectral characteristics of the hyperspectral data cube and affects the final anomaly detection performance. Summary of the Invention

[0007] The present invention provides a spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation. A spatial and spectral joint low-rank constraint regular term is designed and added to the tensor low-rank decomposition model. The optimal anomaly tensor is decomposed to obtain the final anomaly detection result, which can effectively achieve high-precision detection and recognition of anomalies in different complex backgrounds on the basis of protecting the inherent spatial-spectral characteristics of the hyperspectral data cube.

[0008] A spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation includes the following steps:

[0009] (1) Let the original hyperspectral image to be detected be where represents the real number field, H and W are the spatial dimensions of the hyperspectral image, and B is the number of bands; the original hyperspectral image is decomposed to obtain a background tensor and an anomaly tensor Then the original hyperspectral image is expressed as:

[0010] (2) According to the two-way low-rank characteristics of the spaceborne hyperspectral image, two-dimensional low-rank constraints in the spatial and spectral dimensions and spatial local smoothness constraints are imposed on the background tensor , and tube-wise sparse constraints are imposed on the anomaly tensor to obtain the following optimization model

[0011]

[0012]

[0013] where is the background tensor Measure the joint low-rank property of the spatial and spectral dimensions, which is used to simultaneously constrain the spatial low-rank property and the spectral low-rank property of the background tensor; For the background tensor Spatial local smoothness regularization term of For the anomaly tensor Tube sparsity regularization term of Denote the adjustment variable And Minimize the subsequent objective function, For the constraint condition;

[0014] (3) Use the alternating direction method of multipliers to solve the above optimization model to obtain the optimal background tensor And the optimal anomaly tensor

[0015] (4) Calculate the final anomaly detection result R through the optimal anomaly tensor The mathematical expression is as follows:

[0016]

[0017] In the formula, ||·|| F Is the Frobenius norm of the matrix, Denote the optimal anomaly tensor obtained by decomposition The fiber of modulo 3 corresponding to the spatial position i, j in

[0018] Among them, in step (2), the weighted tensor nuclear norm based on t-SVD is adopted to design the Joint low-rank constraint term of the spatial and spectral dimensions of the background tensor, that is, In the optimization model, the mathematical expression is:

[0019]

[0020] Among them, Denote the weighted tensor nuclear norm, Used to constrain the Spatial dimension low-rank property of the background tensor; Used to constrain the Spectral dimension low-rank property of the background tensor, For the background tensor Permutation tensor of

[0021] Among them, the weighted tensor nuclear norm based on t-SVD is defined as follows:

[0022]

[0023] Among them, represents the weighted tensor nuclear norm, used to characterize the low-rank characteristics of all forward slices, is a third-order tensor, and n1, n2, and n3 are the values of the three dimensions of the tensor respectively; ||·|| w,* represents the weighted nuclear norm, is the tensor obtained by performing Fourier transform on each tube of the nth forward slice of used to characterize the low-rank characteristics of is the mth singular value of m is the corresponding weight, n represents the index of the forward slice of the tensor and m represents the index of the singular value of the forward slice of

[0024] Among them, in step (2), based on the linear total variation norm, a spatial local smoothness constraint is imposed on the background tensor i.e., the mathematical expression in the optimization model is: Mathematical expression:

[0025]

[0026] where β is the regularization parameter; ||·|| F is the F norm of the matrix, is the matrix obtained by performing modulo-1 matricization on the background tensor and the mathematical expression is is the matrix obtained by performing modulo-2 matricization on the background tensor and the mathematical expression is is the finite difference operator in the vertical direction, is the finite difference operator in the horizontal direction; where H and W are the spatial dimensions of the hyperspectral image, and B is the number of bands.

[0027] Among them, is the finite difference operator in the vertical direction and is the finite difference operator in the horizontal direction, and the mathematical expression is:

[0028]

[0029]

[0030] In the formula, i represents the operator matrix DH and D W The row index of, and j represents the operator matrix D H and D W The column index of, and others represent all other positions.

[0031] Among them, in step (2), based on l 1,1,2 norm, apply a per-tube sparse constraint to the abnormal tensor That is, in the optimization model The mathematical representation is:

[0032]

[0033] Among them, λ is the regularization parameter, ||·|| 1,1,2 represents the l 1,1,2 norm of the matrix, is the abnormal tensor, i represents the row index of the abnormal tensor in the two-dimensional space, and j represents the abnormal tensor in the two-dimensional space. The column index, represents the abnormal tensor obtained by decomposition The corresponding spatial position in is the modulo-3 fiber at i, j, and the modulo-3 fiber is the spectral vector in the hyperspectral image.

[0034] Among them, the specific process of step (3) is:

[0035] First, bring and The mathematical representation of into the optimization model, and then we get

[0036]

[0037]

[0038] Secondly, introduce auxiliary variables and respectively constrain the background tensor to make the variables in the objective function in the above optimization model separable. Then the optimization model is rewritten as:

[0039]

[0040]

[0041] Thirdly, according to the rewritten optimization model above, perform augmented Lagrangian function transformation to get the following:

[0042]

[0043] Among them, and is the Lagrange multiplier, α is the balance parameter, and μ is the penalty parameter;

[0044] Finally, solve iteratively according to the solution framework of the alternating direction multiplier method. By updating one variable in the model and fixing the other variables, the above variables are transformed into the solution of independent sub-problems. After the (k + 1)-th iteration, each variable is updated as follows: Update to obtain the sub-problem Update to obtain the sub-problem Update to obtain the sub-problem Update to obtain the sub-problem Update to obtain the sub-problem Update to the sub-problem

[0045] For the Lagrange multiplier and the penalty parameter μ, the update is as follows:

[0046]

[0047]

[0048]

[0049]

[0050]

[0051] μ k+1 = min(ρμ k , μ max )

[0052] Continuously repeat the above variable update process until the iteration stops, to obtain the optimal background tensor and the optimal anomaly tensor

[0053] Among them, the process of updating to obtain the sub-problem is as follows:

[0054]

[0055] Among them, is after the k-th iteration of is after the (k + 1)-th iteration of is after the k-th iteration of is after the k-th iteration of μk is μ after the k-th iteration;

[0056] Update to obtain the sub-problem The process is as follows:

[0057]

[0058] where, is after the k-th iteration is after the k-th iteration k is μ after the k-th iteration;

[0059] Update to obtain the sub-problem The process is as follows:

[0060]

[0061] where, is after the k-th iteration is after the k-th iteration k is μ after the k-th iteration, is the modulo 1 matrixification matrix of;

[0062] Update to obtain the sub-problem The process is as follows:

[0063]

[0064] where, is after the k-th iteration is after the k-th iteration k is μ after the k-th iteration, is the modulo 2 matrixification matrix of;

[0065] Update to obtain the sub-problem The process is as follows:

[0066]

[0067]

[0068] Among them, is after the k-th iteration is after the update of the (k + 1)-th iteration is after the update of the (k + 1)-th iteration is after the update of the (k + 1)-th iteration is after the update of the (k + 1)-th iteration is after the update of the k-th iteration are respectively after the update of the k-th iteration μ k is μ after the k-th iteration.

[0069] Update to the sub-problem The process is as follows:

[0070]

[0071] Among them, is after the (k + 1)-th iteration is after the k-th iteration is after the update of the k-th iteration μ k is μ after the k-th iteration.

[0072] Compared with the prior art, the present invention has the following beneficial effects:

[0073] 1. Represent the hyperspectral image with a tensor, decompose it into the sum of a background tensor and an anomaly tensor, use the weighted tensor nuclear norm based on t-SVD, design a joint low-rank constraint regular term for the background tensor, and at the same time impose constraints on the low-rank characteristics of the spatial dimension and the spectral dimension, which can effectively characterize the background information and complete the estimation of the background tensor; in addition, imposing a spatially piecewise smoothness constraint on the background tensor helps to reduce the inclusion of abnormal pixels in the background tensor, enabling the anomaly tensor to capture as many abnormal pixels as possible;

[0074] 2. The anomaly tensor should only contain abnormal pixels, that is, only the positions of abnormal pixels are non-zero, and this characteristic makes the anomaly tensor sparse. The present invention is based on the l 1,1,2 norm, and imposes a per-tube sparsity constraint on the anomaly tensor, which can make the anomaly tensor contain as few background pixels as possible, thereby improving the detection performance;

[0075] 3. By using the alternating direction method of multipliers for iterative solution, it is possible to decompose and obtain the optimal background tensor and the optimal anomaly tensor results on the basis of protecting the inherent spatial-spectral characteristics of the hyperspectral data cube, and obtain the final anomaly detection result according to the optimal anomaly tensor obtained by the decomposition, effectively realizing the high-precision detection of anomalies;

[0076] 4. The method of the present invention can effectively realize the anomaly detection in different complex backgrounds such as the ocean, land, forest, and gobi, laying an important foundation for subsequent precise target detection and recognition. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1 is a spaceborne pushbroom imaging hyperspectral image containing typical backgrounds and anomaly targets;

[0078] Figure 2 is a flow chart of a spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation according to an embodiment of the present invention;

[0079] Figure 3 is a schematic diagram of the principle of spaceborne pushbroom hyperspectral imaging;

[0080] Figure 4 is a schematic diagram of tensor permutation operation;

[0081] Figure 5 is the hyperspectral image to be detected in an embodiment of the present invention;

[0082] Figure 6 is the anomaly detection result of the spaceborne hyperspectral image in an embodiment of the present invention;

[0083] Figure 7 is the detection results of different methods in the visible and near-infrared bands in an embodiment of the present invention;

[0084] Figure 8 is the detection results of different methods in the short-wave infrared band in an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0085] The present invention will be further described in detail below with reference to the drawings and embodiments. It should be noted that the following embodiments are intended to facilitate the understanding of the present invention and do not limit it in any way.

[0086] The present invention proposes a spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation, which combines the inherent characteristics of the pushbroom hyperspectral imaging payload image and can effectively realize the detection and recognition of anomaly targets in different backgrounds such as the ocean, land, forest, and gobi, including single pixel-level targets, multi pixel-level targets, and regional-level targets. As Figure 1As shown, the spaceborne pushbroom imaging hyperspectral images (from the hyperspectral imagers carried by GF-5 satellite and ZY1-02D satellite) contain different abnormal targets under several typical backgrounds. The anomaly detection method proposed by the present invention represents the hyperspectral image with a third-order tensor, and decomposes the original hyperspectral image tensor into the sum of a background tensor and an anomaly tensor. According to the low-rank characteristics of the spaceborne pushbroom imaging hyperspectral image in both spatial and spectral dimensions, a joint low-rank constraint regular term and a local smoothness constraint are designed for the background tensor in the tensor low-rank decomposition model, and a per-tube sparsity constraint is imposed on the anomaly tensor. Then, the above model is iteratively solved using the alternating direction multiplier method to obtain the optimal decomposition result. Finally, the anomaly detection result is obtained based on the optimal anomaly tensor obtained by the decomposition. This method can effectively achieve high-precision anomaly detection in the complex background of spaceborne pushbroom imaging hyperspectral images, laying an important foundation for subsequent accurate target detection and recognition.

[0087] An anomaly detection method for spaceborne hyperspectral images based on joint low-rank tensor approximation, as Figure 2 shown, includes the following:

[0088] Model construction:

[0089] 1) Let the hyperspectral image to be detected be where H and W are the spatial dimensions of the hyperspectral image, B is the number of bands, denotes the real number field. We assume that the original hyperspectral image can be decomposed into a background tensor and an anomaly tensor the sum of, expressed as In the ideal optimal decomposition result, the background tensor only contains background pixels; the anomaly tensor only contains anomaly pixels, that is, only the tubes corresponding to the anomaly pixels, namely the mode-3 fibers in the third-order tensor, and the spectral vectors corresponding to the hyperspectral image are non-zero. In order to achieve the above decomposition effect, we design a joint low-rank constraint regular term and a local smoothness constraint regular term for the background tensor, and impose a per-tube sparsity constraint on the anomaly tensor to obtain the following optimization model:

[0090]

[0091]

[0092] where, is a measure of the joint low-rank characteristics of the spatial and spectral dimensions of the background tensor , used to simultaneously constrain the spatial low-rank characteristics and spectral low-rank characteristics of the background tensor; is the spatial local smoothness regular term of the background tensor , Is an abnormal tensor of the tube sparsity regularization term; Denotes the adjustment variable and minimizes the subsequent objective function, is the constraint condition;

[0093] 2) The joint low-rank regularization term of the background tensor based on tensor singular value decomposition (t-SVD) in the optimization model For the background tensor Each band corresponds to a forward slice of a third-order tensor. The background of a hyperspectral image usually consists of various types of ground objects, and the spectral characteristics of the same type of ground object are similar. This non-local similarity in the spatial domain results in the background tensor showing low-rank characteristics on the forward slices. For a hyperspectral image in pushbroom imaging mode, only constraining the forward low-rank characteristics of the background tensor is insufficient. As Figure 3 shown, during the pushbroom imaging process, each frame of two-dimensional data acquired by the detector corresponds to the spectral information of a ground strip, corresponding to the horizontal slice in the hyperspectral image tensor. When the satellite flies, multiple consecutive strips form a complete hyperspectral image. Such an imaging process results in the background tensor having low-rank characteristics different from the forward slices on the horizontal slices. Therefore, it is necessary to simultaneously constrain the spatial low-rank characteristics and spectral low-rank characteristics of the background tensor to accurately characterize the characteristics of the background tensor and complete the estimation of the background tensor. The tensor nuclear norm based on t-SVD can be used as a measure of the rank of a third-order tensor. However, considering that the tensor nuclear norm treats eigenvalues of different sizes equally, it will cause larger eigenvalues to decrease more and lose the main information of the background. Therefore, we use the weighted tensor nuclear norm instead of the tensor nuclear norm to avoid the above problems. The weighted tensor nuclear norm based on t-SVD is defined as follows:

[0094]

[0095] where, denotes the weighted tensor nuclear norm, is used to characterize the low-rank characteristics of all forward slices,, is a third-order tensor, and n1, n2, and n3 are the values of the three dimensions of the tensor respectively; ||·|| w,* denotes the weighted nuclear norm, is for the tensor obtained by performing Fourier transform on each tube of the nth forward slice of is used to characterize the low-rank characteristics of is the mth singular value ofm For the corresponding weight, where n represents the tensor index of the forward slice, and m represents the forward slice index of the singular value.

[0096] Based on the above weighted tensor nuclear norm, we impose a bidirectional low-rank constraint on the background tensor, and the optimization model is rewritten as:

[0097]

[0098]

[0099] where is used to constrain the spatial low-rank property of the background tensor ; is the dimension permutation tensor of the background tensor , obtained by rearranging the horizontal slices, and its forward slice corresponds to the horizontal slice, as shown in Figure 4 ; is used to constrain the spectral low-rank property of the background tensor ; α is a balancing parameter used to balance the degree of the two-dimensional low-rank constraints.

[0100] 3) Optimize the local smoothness regularization term of the background tensor based on the linear total variation norm in the optimization model. Background pixels in hyperspectral images are usually continuously distributed, and adjacent pixels exhibit a certain degree of spatial uniformity, which makes the background tensor have the characteristic of spatial piecewise smoothness. Imposing a spatial piecewise smoothness constraint on the background tensor helps to reduce the inclusion of abnormal pixels in the background tensor, enabling the abnormal tensor to capture as many abnormal pixels as possible, thereby improving the detection performance. Based on the linear total variation norm, we impose a spatial local piecewise smoothness constraint on the background tensor, that is, the in the optimization model, which is mathematically expressed as:

[0101]

[0102] where β is the regularization parameter; is the matrix obtained by performing modulo-1 matrixization (arranging column by column with modulo-1 fibers) on the background tensor , and is mathematically expressed as is the matrix obtained by performing modulo-2 matrixization (arranging column by column with modulo-2 fibers) on the background tensor , and is mathematically expressed as ||·|| F is the Frobenius norm of the matrix (usually called the F-norm); and are finite difference operators in the vertical and horizontal directions, and their mathematical expressions are:

[0103]

[0104]

[0105] Substitute into the aforementioned optimization model, and we get:

[0106]

[0107]

[0108] 4) The tube sparsity regular term of the l 1,1,2 -norm in the optimization model. The abnormal tensor should only contain abnormal pixels, that is, only the positions of abnormal pixels are non-zero, and this characteristic makes the abnormal tensor sparse. At the same time, the abnormal targets in the hyperspectral image should have complete spectral vectors, that is, the non-zero pixels in the abnormal tensor should be in the unit of tube fibers. In the present invention, the l 1,1,2 -norm is used to constrain the per-tube sparsity of this abnormal tensor, that is, in the optimization model The mathematical expression is:

[0109]

[0110] where λ is the regularization parameter. Substitute into the aforementioned model, and we get:

[0111]

[0112]

[0113] Solving this model can obtain the optimal decomposition result, and thus obtain the anomaly detection result.

[0114] Model solving:

[0115] 5) Use the alternating direction multiplier method to solve the above optimization model. First, introduce the auxiliary variable to make the variables in the objective function of the optimization model separable, and the optimization model is rewritten as:

[0116]

[0117] Then, according to the above constrained optimization model, the augmented Lagrangian function is obtained as follows:

[0118]

[0119] where and are Lagrange multipliers, and μ is the penalty parameter. According to the solution framework of the alternating direction method of multipliers, in the iterative solution process, by fixing other variables and updating a certain variable, the above problem is transformed into the solution of independent subproblems. Specifically, for the (k + 1)-th iteration, the variables are updated in turn as follows:

[0120] For Fixing other variables gives the following subproblem:

[0121]

[0122] where is after the k-th iteration is after the (k + 1)-th iteration is after the k-th iteration μ k is μ after the k-th iteration. This subproblem is a problem of minimizing the weighted tensor nuclear norm and is solved by the weighted tensor singular value thresholding algorithm. The solution process is as follows: Let First, perform a Fourier transform on each mode-3 fiber of to obtain Then, perform a singular value decomposition on each forward slice of For the i-th forward slice of the singular value decomposition gives Perform a singular value shrinking operation on the singular value matrix S, that is, subtract a parameter from each singular value, and the size of the parameter is where w is the weight corresponding to each singular value. After the singular value shrinking is completed, multiply the updated singular value matrix by the original U and V. After all the forward slices have undergone the above operations, perform an inverse Fourier transform on the mode-3 fibers of the resulting new tensor to obtain the updated where w is the weight corresponding to each singular value. After the singular value shrinking is completed, multiply the updated singular value matrix by the original U and V. After all the forward slices have undergone the above operations, perform an inverse Fourier transform on the mode-3 fibers of the resulting new tensor to obtain the updated

[0123] For Fixing other variables gives the following subproblem:

[0124]

[0125] where is after the k-th iteration is after the k-th iteration μ k is μ after the k-th iteration. This subproblem is a problem of minimizing the weighted tensor nuclear norm and is solved by the weighted tensor singular value thresholding algorithm, similar to the solution of The solution process is as follows:

[0126] Let Perform a tensor permutation operation on it to obtain This permutation operation will rearrange the slices of, and the obtained forward slice of is the original horizontal slice of. Then, perform the same operation as when solving for on to obtain the slice-permuted Finally, perform a permutation operation on to obtain the updated

[0127] For Fixing other variables, the following sub-problem can be obtained:

[0128]

[0129] where is after the k-th iteration is after the k-th iteration μ k is μ after the k-th iteration, is after the k-th iteration mod-1 matrixized matrix of. Let Perform mod-1 matrixization on it (arrange its mod-1 fibers into a matrix by columns) to obtain The closed-form solution of the above sub-problem can be obtained: Rearrange it into a tensor to obtain the updated

[0130] For Fixing other variables, the following sub-problem can be obtained:

[0131]

[0132] where is after the k-th iteration is after the k-th iteration μ k is μ after the k-th iteration, is after the k-th iteration mod-2 matrixized matrix of. Similar to the solution of Let Perform mod-2 matrixization on it (arrange its mod-2 fibers into a matrix by columns) to obtain The closed - form solution of the above sub - problem can be obtained: Rearrange it into a tensor to get the updated

[0133] For Fixing other variables, the following sub - problem can be obtained:

[0134]

[0135] This sub - problem has a closed - form solution as follows:

[0136]

[0137] For Fixing other variables, the following sub - problem can be obtained:

[0138]

[0139] Let The closed - form solution of this sub - problem is as follows:

[0140]

[0141] For the Lagrange multiplier and the penalty parameter μ, the update is as follows:

[0142]

[0143]

[0144]

[0145]

[0146]

[0147] μ k+1 = min(ρμ k , μ max )

[0148] The judgment basis for iterative stopping is to satisfy the following conditions or reach the maximum number of iterations:

[0149]

[0150] where ξ = 10 -6 , and the maximum number of iterations is set to 500.

[0151] Repeat the above update process of each variable until the iteration stops, and the optimal decomposition result can be obtained: the optimal background tensor and the optimal anomaly tensor Through Calculate the final anomaly detection result R, which is mathematically expressed as follows:

[0152]

[0153] To further demonstrate the technical effects of the present invention, a hyperspectral image obtained by an advanced hyperspectral imager (AHSI) carried on the GF-5 satellite is selected as an example to specifically implement the anomaly detection method proposed by the present invention. The hyperspectral imager uses a push-broom imaging method.

[0154] As Figure 5 shown, a hyperspectral image includes a visible and near-infrared band image and a short-wave infrared image. The spatial dimension of the visible and near-infrared band image contains 100×100 pixels, and the number of bands is 150. The short-wave infrared image contains 100×100 pixels, and the number of bands is 180. Perform anomaly detection on the visible and near-infrared band image and the short-wave infrared image respectively with the same process to obtain the visible and near-infrared anomaly detection result and the short-wave infrared anomaly detection result. Assume that the hyperspectral image to be detected is The specific implementation process is as follows:

[0155] 1) Initialize the background tensor and the anomaly tensor: Initialize the auxiliary variable: Initialize the Lagrange multiplier: Initialize each parameter: The initial penalty parameter μ0 of the Lagrange penalty is 10 -2 , the self-adaptive update coefficient ρ of the Lagrange penalty parameter is 1.1, and the iteration stop judgment parameter ξ is 10 -6 , the initial iteration number k = 0, the joint low-rank constraint balance parameter α = 0.6, the smoothness regularization parameter β = 0.25, and the tube-by-tube sparsity regularization parameter λ = 1. For the detection of different images, the parameter α can be appropriately adjusted within the interval [0,1] according to the low-rank characteristics of the spatial and spectral dimensions of the hyperspectral image, and the parameter λ can be appropriately adjusted within the interval [0.1,3] according to the size and quantity of the abnormal targets of interest to obtain the best detection result.

[0156] 2) According to the framework of the alternating direction method of multipliers, start iterative optimization. Assume that k iterations have been performed to obtain

[0157] For the (k + 1)-th iteration, update the auxiliary variable Let First, perform a Fourier transform on each modulo-3 fiber of to obtain Then, perform a singular value decomposition on each forward slice of For the i-th forward slice of , it is expressed as Perform singular value decomposition to obtain Perform a singular value shrinking operation on the singular value matrix S: subtract a parameter from each singular value, where the parameter size is where w m is the weight corresponding to the m-th singular value, and in the present invention, this weight is set to is the m-th eigenvalue of the n-th forward slice , ε = 10 -1 is a fixed constant to avoid the case of a zero denominator. After the singular value shrinking is completed, multiply the updated singular value matrix by the original U and V to obtain a new forward slice. After performing the above operations on all the forward slices, a new tensor is obtained For the obtained new tensor perform an inverse Fourier transform on the modulo 3 fibers to obtain the k + 1 iteratively updated

[0158] 3) Update the auxiliary variable

[0159] Let Perform a tensor permutation operation on it to obtain This permutation operation will rearrange the slices of such that the forward slices of the obtained are the original horizontal slices. Then, perform the same operations as those used in updating on to obtain Finally, perform a permutation operation on to obtain the updated

[0160] 4) Update the auxiliary variable

[0161] Let First, perform a modulo 1 matrixification on it: arrange its modulo 1 fibers in columns to form a matrix, obtaining Then, the matrix form of the closed-form solution of the above sub-problem can be obtained as follows: Finally, rearrange it into a tensor to obtain the updated

[0162] 5) Update the auxiliary variable

[0163] Similar to the update of , let First, perform a modulo 2 matrixification on it: arrange its modulo 2 fibers in columns to form a matrix, obtaining Then, the closed-form solution matrix of the above sub-problem can be obtained: Finally, rearrange it into a tensor to obtain the updated

[0164] 6) Update the background tensor

[0165] According to the above update Update the background tensor at the (k + 1)-th iteration through the following formula

[0166]

[0167] 7) Update the anomaly tensor

[0168] According to the update Let Update the anomaly tensor at the (k + 1)-th iteration through the following formula

[0169]

[0170] 8) Update the Lagrange multiplier tensor and the Lagrange penalty parameter μ k+1 :

[0171] After updating the above variables, update the Lagrange multiplier tensor and the penalty parameter through the following formula:

[0172]

[0173]

[0174]

[0175]

[0176]

[0177] μ k+1 = min(ρμ k , μ max )

[0178] 9) Determine whether the iteration process stops:

[0179] Condition 1: and and

[0180] Condition 2: The number of iterations k > 100

[0181] If any of the above conditions is satisfied, the iterative process is stopped, and the optimal background tensor and the optimal anomaly tensor obtained by decomposition are output. and the optimal anomaly tensor

[0182] 9) Complete anomaly detection:

[0183] According to the optimal anomaly tensor obtained by decomposition The final anomaly detection result is obtained through the following formula:

[0184]

[0185] So far, the anomaly detection method proposed by the present invention has been implemented on a hyperspectral image, and the detection results are as Figure 6 shown.

[0186] Finally, the proposed method of the present invention is compared with five commonly used anomaly detection methods for this hyperspectral image. The five detection methods for comparison are: the RX detection method based on statistics, from the article "Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution"; the LRASR detection method based on representation, from the article "Anomaly Detection in Hyperspectral Images Based on Low-Rank and Sparse Representation"; the LSMAD detection method based on matrix decomposition, from the article "A Low-Rank and Sparse Matrix Decomposition-Based Mahalanobis Distance Method for Hyperspectral Anomaly Detection"; the PTA detection method based on tensor decomposition, from the article "Prior-Based Tensor Approximation for Anomaly Detection in Hyperspectral Imagery"; the RGAE detection method based on deep learning, from the article "Hyperspectral Anomaly Detection With Robust Graph Autoencoders".

[0187] It can be seen that the comparison results in the near-infrared band and the short-wave infrared band are respectively as Figure 7 and Figure 8As can be seen, the proposed method can effectively achieve high-precision detection of anomalies in the complex background of spaceborne pushbroom imaging hyperspectral images, and the detection results are significantly better than other methods.

[0188] The above-described embodiments have elaborated on the technical solutions and beneficial effects of the present invention. It should be understood that the above are only specific embodiments of the present invention and are not used to limit the present invention. Any modifications, supplements, and equivalent replacements made within the scope of the principles of the present invention shall be included within the protection scope of the present invention.

Claims

1. A spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation, characterized in that, It includes the following steps: (1) Let the original hyperspectral image to be detected be where represents the real number field, H and W are the spatial dimensions of the hyperspectral image, and B is the number of bands; the original hyperspectral image is decomposed to obtain a background tensor and an anomaly tensor Then the original hyperspectral image is expressed as: (2) According to the bidirectional low-rank characteristics of spaceborne hyperspectral images, for the background tensor impose low-rank constraints in both the spatial and spectral dimensions and a spatial local smoothness constraint, and for the anomaly tensor impose a per-tube sparsity constraint to obtain the following optimization model Among them, is a measure of the joint low-rank property of the spatial dimension and spectral dimension of the background tensor, and is used to simultaneously constrain the spatial low-rank property and spectral low-rank property of the background tensor; is a measure of the joint low-rank property of the spatial dimension and spectral dimension of the background tensor, and is used to simultaneously constrain the spatial low-rank property and spectral low-rank property of the background tensor; is the background tensor is the spatial local smoothness regularization term of the background tensor is the anomaly tensor is the tube sparsity regularization term of the anomaly tensor; represents adjusting the variables and to minimize the subsequent objective function, is the constraint condition; Optimizing the Tube Sparsity Regularization Term in the Model The mathematical representation is as follows: where λ is the regularization parameter, and ||·|| 1,1,2 represents the l 1,1,2 norm of the matrix, and ||·|| F represents the Frobenius norm of the matrix, is the anomaly tensor, i represents the row index of the anomaly tensor in the two-dimensional space, and j represents the anomaly tensor in the two-dimensional space, and the column index, represents the anomaly tensor obtained by decomposition and the corresponding spatial position in it is the fiber of modulus 3 at i, j. The fiber of modulus 3 is the spectral vector in the hyperspectral image; (3) Solve the above optimization model using the alternating direction multiplier method to obtain the optimal background tensor and the optimal anomaly tensor (4) Through the optimal anomaly tensor Calculate the final anomaly detection result R, which is mathematically expressed as follows: where ||·|| F is the Frobenius norm of the matrix, represents the optimal anomaly tensor obtained by decomposition and the fiber of modulo 3 at the spatial position of i, j in it, and the fiber of modulo 3 is the spectral vector in the hyperspectral image; where i represents the row index of the detection result R, and j represents the column index of the detection result R.

2. The spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation according to claim 1, wherein, In step (2), a weighted tensor nuclear norm based on t-SVD is adopted to design a joint low-rank constraint term for the spatial dimension and spectral dimension of the background tensor , that is, the in the optimization model, and the mathematical representation is as follows: Among them, represents the weighted tensor nuclear norm, which is used to constrain the low-rank property of the spatial dimension of the background tensor ; is used to constrain the low-rank property of the spectral dimension of the background tensor ; is the permutation tensor of the background tensor ; α is the balance parameter, which is used to balance the degree of the two low-rank constraints.

3. The spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation according to claim 2, characterized in that The weighted tensor nuclear norm based on t-SVD is defined as follows: Among them, represents the weighted tensor nuclear norm, which is used to characterize the low-rank characteristics of all forward slices, is a third-order tensor, and n1, n2, and n3 are the values of the three dimensions of the tensor respectively; ‖·‖ w,* represents the weighted nuclear norm, which is the tensor obtained by performing Fourier transform on each tube of the nth forward slice of used to characterize the low-rank characteristics of is the mth singular value of m and w is the corresponding weight, where n represents the index of the forward slice of the tensor and m represents the index of the singular value of the forward slice .

4. The spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation according to claim 1, wherein In step (2), based on the linear total variation norm, for the background tensor apply the spatial local smoothness constraint, that is, in the optimization model The mathematical representation is: where β is the regularization parameter; is the background tensor The matrix obtained by modulo-1 matrixization, mathematically expressed as is the background tensor The matrix obtained by modulo-2 matrixization, mathematically expressed as is the finite difference operator in the vertical direction, is the finite difference operator in the horizontal direction; where H and W are the spatial dimensions of the hyperspectral image, and B is the number of bands.

5. The spaceborne hyperspectral image anomaly detection method based on joint low-rank tensor approximation according to claim 4, characterized in that is a finite difference operator in the vertical direction and is a finite difference operator in the horizontal direction, and the mathematical representation is as follows: wherein, i represents the row index of the operator matrices D H and D W , j represents the column index of the operator matrices D H and D W , and others represent all other positions.

6. The on-orbit hyperspectral image anomaly detection method based on joint low-rank tensor approximation according to claim 1, wherein The specific process of step (3) is as follows: First, after substituting the mathematical representations of and into the optimization model, we obtain Secondly, introduce auxiliary variables Constrain the background tensors respectively, so that the variables in the objective function of the above optimization model can be separated, and then the optimization model is rewritten as: Again, according to the above rewritten optimization model, augmented Lagrangian function transformation is performed to obtain the following: wherein, and are Lagrange multipliers, α is a balance parameter, and μ is a penalty parameter; Finally, iterative solution is performed according to the solution framework of the alternating direction multiplier method. By fixing other variables while updating a certain variable in the model, the above variables are transformed into the solution of independent sub-problems. After the (k + 1)-th iteration, each variable is updated in turn as follows: Update to obtain the sub-problem Update to obtain the sub-problem Update to obtain the sub-problem Update to obtain the sub-problem Update to obtain the sub-problem Update to the sub-problem For the Lagrange multipliers and the penalty parameter μ, the update is as follows: μ k+1 = min(ρμ k , μ max ) Continuously repeat the above variable update processes until the iteration stops to obtain the optimal background tensor and the optimal anomaly tensor μ k is the Lagrange multiplier and the result of the k-th iteration of the penalty parameter μ.

7. The on-orbit hyperspectral image anomaly detection method based on joint low-rank tensor approximation according to claim 6, characterized in that Update Obtain sub-problems The process is as follows: Among them, is after the k-th iteration is after the (k + 1)-th iteration is after the k-th iteration is after the k-th iteration μ k is μ after the k-th iteration; Update Obtain sub-questions The process is as follows: Among them, is after the k-th iteration is after the k-th iteration is after the k-th iteration μ k is μ after the k-th iteration; Update Obtain a sub-problem The process is as follows: Among them, is after the k-th iteration is after the k-th iteration is after the k-th iteration μ k is μ after the k-th iteration, is after the k-th iteration mod 1 matrix of the matrix; Update Obtain a sub-problem The process is as follows: Among them, is after the k-th iteration is after the k-th iteration is after the k-th iteration μ k is μ after the k-th iteration, is after the k-th iteration the modulo 2 matrixification matrix of; Update Obtain sub-problems The process is as follows: Among them, is after the k-th iteration is after the update of the (k + 1)-th iteration is after the update of the (k + 1)-th iteration is after the update of the (k + 1)-th iteration is after the update of the (k + 1)-th iteration is after the update of the k-th iteration are respectively after the update of the k-th iteration μ k is μ after the k-th iteration; Update to sub-problems The process is as follows: Among them, is after the (k + 1)-th iteration is after the k-th iteration is after the k-th iteration update μ k is μ after the k-th iteration.

Images