1. A hyperspectral anomaly detection method based on spatio-spectral feature joint constraint is characterized in that,

step 1: establishment based on S_1/2The norm and denoised spectral information low-rank model comprises the following sub-steps:

step 1.1: let any one of the original hyperspectral images be represented asH and W respectively represent the height and width of each wave band in the hyperspectral image, and d represents the number of the wave bands of the hyperspectral image;

step 1.2: carrying out tensor matrixing on omega according to wave band dimensionality to obtain a two-dimensional matrixWherein d is the number of wave bands, n represents the number of pixels, and n is H multiplied by W in number;

step 1.3: modeling the matrixed hyperspectral D to obtain

D＝Z+S+G

Wherein the content of the first and second substances,a matrix representing the background is generated,a matrix of anomalies is represented that is,representing a noise matrix;

step 1.4: establishing S_1/2Low rank representation model:

s.t.D＝Z+S+G

wherein, minimizeCharacterized by the low rank of the background matrix Z, andthe sparsity of the anomaly matrix is represented;

step 1.5: for S in step 1.4_1/2The low rank representation model imposes a denoising constraint term,

step 1.6: establishing an optimization objective function:

step 2: carrying out variation constraint on the background matrix Z to obtain a spatial information constraint model based on total variation regularization:

SSTV(Z)＝||D_hZD_s||₁+||D_vZD_s||₁

wherein the operatorIs a spatial domain variation operator and is a spatial domain variation operator,is a spectral domain variation operator; the 1 norm of the matrix is the sum of the absolute values of the elements in the matrix.

And step 3: establishing a final hyperspectral anomaly detection optimization objective function according to the spectral information low-rank model obtained in the step 1 and the spatial information constraint model obtained in the step 2 as follows:

wherein λ is₁,λ₂,λ₃Is a coefficient used to balance the importance of items;

and 4, step 4: establishing a Lagrange function according to the hyperspectral anomaly detection optimization target function obtained in the step 3, fixing all other variables to solve when one variable is optimized based on an alternating direction multiplier method, obtaining an optimized background matrix Z and an optimized anomaly matrix S, and calculating l of each row of the anomaly matrix S₂The norm, when greater than a given threshold, may be determined to be abnormal.

2. The method for hyperspectral anomaly detection based on spatio-spectral feature joint constraint according to claim 1, characterized in that for the spatial domain variation operator, substantially for each band (D)ⁱH × W on i ═ 1, 2., d), i.e., n pixels, are subjected to a variation operation. The spatial domain variation operator D when calculated under a two-dimensional matrix data model expanded along the spectral dimension_hAnd D_vAll-to-two-dimensional hyperspectral matrixIs differentiated by building a different operator matrix, D_hThe operators differentiating horizontally adjacent pixels, D, in each band_vThe operators differentiate vertically adjacent pixels on each band. Aiming at the variation operator of the spectral domain, the variation operator is essentially to perform variation operation on different wave bands of the same pixel, and when the calculation is performed under a two-dimensional matrix data model stretched and expanded along the spectral dimension, the spectral domain variation operator D_sIs to a two-dimensional hyperspectral matrixThe rows are differentiated.

3. The hyperspectral anomaly detection method based on spatio-spectral feature joint constraint according to claim 1, wherein in the step 3, the coefficient is [0.0001,0.001,0.01,0.1,1,10,100,1000 ].

4. The hyperspectral anomaly detection method based on spatio-spectral feature joint constraint according to claim 1, wherein the method is qualitatively and quantitatively analyzed by using ROC curve and AUC value.

Background

In the field of hyperspectral image processing, target detection without acquiring target spectrum prior is called abnormal detection. In view of the fact that the target prior acquisition has certain difficulty and limitation, the unsupervised anomaly detection algorithm without prior has wider application scenes and more distinct practical value. An "anomaly" in a hyperspectral image is a small portion of other pixels defined as having significant spectral feature differences from the vast majority of background pixels. For example, on a large area of clear apron obtained by hyperspectral camera photography, one or two separate airplanes may be defined as anomalous targets in the hyperspectral image, while other parts such as runways and the like may be considered as backgrounds.

Xuchao, Zhantianming (hyper-spectral anomaly detection method based on low-rank total variation regularization, computer science and exploration, 2020, 14 (12): 2140-2149) enhances anomaly detection results by establishing a low-rank decomposition model and applying total variation regularization constraint on spectral dimensions. However, the method adopts the kernel norm approximate matrix rank, has scaling error, and adopts S with higher approximation degree_1/2The norm approximates the matrix rank, so that the precision is improved; in the paper, only the spectral dimension is subjected to variation constraint, and the utilization of the spatial domain information in the hyperspectral image is ignoredSpatial information in the hyperspectral image is obtained; meanwhile, in the modeling process, the noise in the image is considered, and the noise is restrained in the solving process, so that the anomaly detection effect is further enhanced.

Disclosure of Invention

The technical problem solved by the invention is as follows: in order to reduce the influence of noise on hyperspectral anomaly detection and overcome the defect that spatial domain information is not considered in the prior art, the invention provides a hyperspectral anomaly detection method based on spatial-spectral feature joint constraint_1/2The norm replaces the nuclear norm to approximate the matrix rank, so that the approximation precision is improved; in order to reduce the influence of noise on hyperspectral anomaly detection, a denoising constraint term is further added, the anomaly detection precision is improved, and full mining and high utilization of spectral information are realized; the method further adopts a total variation regularization constraint term, and introduces variation operators on 2 spatial dimensions and 1 spectral dimension respectively, so that the utilization rate of the spatial information of the hyperspectral image is improved while the image is smooth.

The technical scheme of the invention is as follows: a hyperspectral anomaly detection method based on spatial-spectral feature joint constraint includes the following steps: establishment based on S_1/2The norm and denoised spectral information low-rank model comprises the following sub-steps:

step 1.1: let any one of the original hyperspectral images be represented asH and W respectively represent the height and width of each wave band in the hyperspectral image, and d represents the number of the wave bands of the hyperspectral image;

step 1.2: carrying out tensor matrixing on omega according to wave band dimensionality to obtain a two-dimensional matrixWherein d is the number of wave bands, n represents the number of pixels, and n is H multiplied by W in number;

step 1.3: modeling the matrixed hyperspectral D to obtain

D＝Z+S+G

Wherein the content of the first and second substances,a matrix representing the background is generated,a matrix of anomalies is represented that is,representing a noise matrix;

step 1.4: establishing S_1/2Low rank representation model:

s.t.D＝Z+S+G

wherein, minimizeCharacterized by the low rank of the background matrix Z, andthe sparsity of the anomaly matrix is represented;

step 1.5: for S in step 1.4_1/2The low rank representation model imposes a denoising constraint term,

step 1.6: establishing an optimization objective function:

step 2: carrying out variation constraint on the background matrix Z to obtain a spatial information constraint model based on total variation regularization:

SSTV(Z)＝||D_hZD_s||₁+||D_vZD_s||₁

wherein the operatorIs a spatial domain variation operator and is a spatial domain variation operator,is a spectral domain variation operator; the 1 norm of the matrix is the sum of the absolute values of the elements in the matrix.

And step 3: establishing a final hyperspectral anomaly detection optimization objective function according to the spectral information low-rank model obtained in the step 1 and the spatial information constraint model obtained in the step 2 as follows:

wherein λ is₁,λ₂,λ₃Is a coefficient used to balance the importance of items;

and 4, step 4: establishing a Lagrange function according to the hyperspectral anomaly detection optimization target function obtained in the step 3, fixing all other variables to solve when one variable is optimized based on an alternating direction multiplier method, obtaining an optimized background matrix Z and an optimized anomaly matrix S, and calculating l of each row of the anomaly matrix S₂The norm, when greater than a given threshold, may be determined to be abnormal.

The further technical scheme of the invention is as follows: for the spatial domain variation operator, substantially for each band (D)ⁱH × W on i ═ 1, 2., d), i.e., n pixels, are subjected to a variation operation. The spatial domain variation operator D when calculated under a two-dimensional matrix data model expanded along the spectral dimension_hAnd D_vAll-to-two-dimensional hyperspectral matrixIs differentiated by building a different operator matrix, D_hOperator differential is on each bandHorizontally adjacent picture elements, D_vThe operators differentiate vertically adjacent pixels on each band. Aiming at the variation operator of the spectral domain, the variation operator is essentially to perform variation operation on different wave bands of the same pixel, and when the calculation is performed under a two-dimensional matrix data model stretched and expanded along the spectral dimension, the spectral domain variation operator D_sIs to a two-dimensional hyperspectral matrixThe rows are differentiated.

The further technical scheme of the invention is as follows: the coefficients take the values of [0.0001,0.001,0.01,0.1,1,10,100,1000 ].

The further technical scheme of the invention is as follows: and carrying out qualitative and quantitative analysis on the method by using the ROC curve and the AUC value.

Effects of the invention

The invention has the technical effects that: the beneficial effects of the invention mainly comprise:

(1) in the traditional low-rank representation model, the nuclear norm of the matrix is consistently adopted, namely the sum of matrix singular values approximates the rank of the matrix, and the invention proposes that S of the matrix is used_1/2The norm approximates the rank of the matrix, namely the sum of the root of the singular value can approximate the rank of the matrix with higher precision relative to the nuclear norm;

(2) the method is used for modeling the noise in a target function aiming at the problem that the influence of the noise on the hyperspectral image anomaly detection is not considered in the traditional hyperspectral anomaly detection algorithm, and only modeling is carried out on the background and the anomaly;

(3) most of the existing hyperspectral anomaly detection algorithms only consider spectral features, neglect of utilization of hyperspectral image spatial features and constraint of spatial domain information, and aiming at the characteristic, the invention applies total variation regularization constraint to spatial dimensions and spectral dimensions, thereby realizing utilization of the hyperspectral image spatial features.

Drawings

FIG. 1 is a method flow diagram:

FIG. 2 is a diagram showing the comparison between the true value graph of the test data set and the detection effect graph of the present invention and other comparison methods, wherein (a) is the true value graph, (b) is the RX method graph, (c) is the LRX method graph, (d) is the CRD method graph, (e) is the LRASR method graph, (f) is the LSMAD method graph, (g) is the GTVLRR method graph, and (h) is the invented method graph

FIG. 3 is a graph of the ROC test results of several comparative methods of the present invention under a test data set, plotted as a threshold cutoff for the effect of models at different thresholds, with the upper curve model being superior to the lower curve enclosed by it.

FIG. 4 is an abnormal object of a building group concentrated in the middle of an image

Detailed Description

In the description of the present invention, it is to be understood that the terms "center", "longitudinal", "lateral", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", "clockwise", "counterclockwise", and the like, indicate orientations and positional relationships based on those shown in the drawings, and are used only for convenience of description and simplicity of description, and do not indicate or imply that the device or element being referred to must have a particular orientation, be constructed and operated in a particular orientation, and thus, should not be considered as limiting the present invention.

Referring to fig. 1-3, the technical scheme of the invention is as follows:

(1) based on Schatten 1/2quasi-norm (S)_1/2Norm) and denoised spectral information low-rank model

For any one matrix of m rows and n columnsIts rank function, S_1/2The norms and nuclear norms are defined as follows:

wherein σ_i(i 1.., min (m, n)) is a singular value of the matrix.

Let any one of the original hyperspectral images be represented asH and W respectively represent the height and width of each wave band in the hyperspectral image, and d represents the number of the wave bands of the hyperspectral image. The method carries out modeling processing on the hyperspectral image on a second-order matrix space, so that tensor matrixing is carried out on the original image according to the wave band dimension to obtain a two-dimensional matrixWhere d is the number of bands, n represents the number of pixels, and n is equal in number to the product of H and W, i.e., n is H × W. Modeling the stretched hyperspectral data matrix D to obtain:

D＝Z+S+G

wherein the content of the first and second substances,a matrix representing the background is generated,a matrix of anomalies is represented that is,representing a noise matrix. The low rank representation model can then adopt S respectively_1/2L of norm and matrix_2,1The norm is a constrained representation of low rank and sparsity. So that S as shown below can be established_1/2Low rank representation model:

s.t.D＝Z+S+G

wherein, minimizeCharacterized by the low rank of the background matrix Z, andi.e. l representing each row of the anomaly matrix S₂And the sum of the norms represents the sparsity of the abnormal matrix S. The low-rank representation model realizes the depiction of the whole structure of the hyperspectral data matrix.

Meanwhile, in order to further enhance the effect of anomaly detection, a denoising constraint term is applied to the model, namely:

the optimization objective function is thus established as follows:

the target function realizes modeling and constraint of hyperspectrum based on spectral information.

(2) Spatial information constraint model based on total variation regularization

The model makes full use of the hyperspectral spectral information, but ignores the use of the hyperspectral image spatial information. The invention carries out variation constraint on a matrix Z aiming at the prior of background data smoothing, and respectively multiplies a variation operator of a spatial domain by a left side and multiplies a variation operator of a spectral domain by a right side, namely:

SSTV(Z)＝||D_hZD_s||₁+||D_vZD_s||₁

wherein the operatorIs a spatial domain variation operator and is a spatial domain variation operator,is a spectral domain variation operator. The concrete description is as follows:

for the spatial domain variation operator, substantially for each band (D)ⁱH × W on i ═ 1, 2., d), i.e., n pixels, are subjected to a variation operation. The spatial domain variation operator D when calculated under a two-dimensional matrix data model expanded along the spectral dimension_hAnd D_vAll-to-two-dimensional hyperspectral matrixIs differentiated by building a different operator matrix, D_hThe operators differentiating horizontally adjacent pixels, D, in each band_vThe operators differentiate vertically adjacent pixels on each band.

Aiming at the variation operator of the spectral domain, the variation operator is essentially to perform variation operation on different wave bands of the same pixel, and when the calculation is performed under a two-dimensional matrix data model stretched and expanded along the spectral dimension, the spectral domain variation operator D_sIs to a two-dimensional hyperspectral matrixThe rows are differentiated.

(3) Hyperspectral anomaly detection model based on spatial-spectral feature joint constraint

Comprehensively considering the spatial domain information and the spectral domain information, establishing a final hyperspectral anomaly detection optimization objective function as follows:

wherein λ is₁,λ₂,λ₃Is a parameter for balancing the importance of each item, i.e. a coefficient, which is determined in an experiment by means of artificial setting, inWithin certain limits to achieve optimal results. The regulation range in the present invention is [0.0001,0.001,0.01,0.1,1,10,100,1000]。

(4) Model optimization and solution

And establishing a Lagrange function according to the objective function and the constraint condition, fixing all the other variables when one variable is optimized based on an alternating direction multiplier method, and specifically solving the process as described in the following specific example. After the optimized background matrix Z and the abnormal matrix S are obtained, calculating l of each row₂Norm, greater than a given threshold, can be determined to be abnormal.

For the hyperspectral anomaly detection method, the detection performance is qualitatively and quantitatively analyzed by adopting a classical Receiver Operating characteristic Curve (ROC Curve) and an Area Under the Curve value (AUC). The ROC curve is formed by drawing a series of target detection rate and false alarm rate point pairs, the relation between two detection probabilities can be directly reflected, a discrimination threshold value is given, a corresponding group of detection rate and false alarm rate values can be accurately calculated, and different threshold values are respectively set to obtain a series of point pairs. In order to carry out accurate quantitative analysis, an area under an ROC curve is integrated to obtain an AUC value, so that the quality of the method is measured, and the performance of the abnormality detection method is intuitively evaluated.

The following takes a certain hyperspectral image of the coastal city background as an example, and the technical scheme of the invention is explained in detail by combining the accompanying drawings. In this data, the background is an image of a city along the coast, and what is considered as an abnormal object is a group of buildings concentrated in the middle area of the image, as shown in the following figure. The algorithm can detect the position of the building group in the urban image under the condition of no spectral prior information, and can be used for further urban planning and construction. H-W-100 and d-207, that is, 100 × 100 images per band, 207 bands in total, and after tensor matrixing expansion, tensor matrixing is obtained as a resultWherein n is 10000 and d is 207, i.e. eachThe wave band comprises 10000 pixels and 207 wave bands in total. The final model is as follows:

namely:

step 1: the model is solved by an alternating direction multiplier method. For solving the objective function, auxiliary variables P, Q, V are introduced₁The model deformation is:

s.t.P＝D_hZD_s,Q＝D_vZD_s,Z＝V₁

the following augmented lagrange function was constructed:

wherein, B₁,B₂,B₃Is an introduced lagrange multiplier.

According to an optimization principle, the variables in the function are respectively optimized in sequence as follows:

1) optimizing variable V₁：The relevant literature proves that the result of the optimization equation is as follows:

V₁ ^*＝H_λ(Z+B₃)＝Udiag(H_λ(σ))V^T

wherein H_λ(σ)＝(h_λ(σ₁),h_λ(σ₂),...,h_λ(σ_r))^TU and V are the pair matrix Z + B₃And (3) performing singular value decomposition to obtain a left singular matrix and a right singular matrix, namely:

Z+B₃＝Udiag(σ)V^T

wherein λ is 1/ν.

2) Optimizing a variable P:

3) optimizing a variable Q:the relevant literature proves that the optimization equation of the above form can be solved by using a soft threshold operator, namely:

4) optimizing variable Z:

according to the relevant literature and the nature of the matrix kronecker product: that is, if the matrix A ═ BCD is known, there areLower case letters respectively represent the vector form of the matrix to which the upper case letters correspond, i.e. the vertical stack of matrix columns. In the above process of optimizing the matrix Z, can orderThe function of the optimization matrix Z is rewritten as:

after derivation, we can obtain:

wherein the content of the first and second substances,and (4) carrying out simulation solution by using an LSQR method.

5) Optimizing a variable S:the relevant literature proves thatThe optimization result is as follows:

step 2: finally obtaining a sparse abnormal matrix S, and calculating l of each row₂And the norm, wherein the pixels corresponding to the corresponding columns with the norm greater than the set threshold eta can be identified as abnormal pixels. And (4) drawing an ROC curve according to different threshold values, so that the performance of the method can be judged.

Table 1 demonstrates the superiority of the model of the invention in quantitative terms and in the form of specific values. The numerical result is an AUC value, which refers to the area of an ROC curve, and the larger the numerical value, the better the model effect

TABLE 1 AUC values