CN115082337A - Hyperspectral mixed noise removing method based on double total variation - Google Patents

Abstract

The present invention provides a method for removing hyperspectral mixed noise based on double total variation, obtains a hyperspectral image containing mixed noise, selects an SSTV regular term to describe the global sparse structure of the hyperspectral image, and at the same time designs a weighted norm to constrain the hyperspectral image The local sparse structure of the spatial direction difference image of the principal component map of the image, introduces double constraints, selects more sparse global features and local features, the sparse regular term is used to isolate sparse noise, and the norm describes the Gaussian noise of the image. Separate hyperspectral image denoising model, introduce auxiliary variables to optimize the hyperspectral image denoising model based on double total variation, use the ALM algorithm framework to optimize the proposed model, and convert complex problems into multiple simple sub-problems to alternate Iterative solution is performed to obtain a denoised hyperspectral image; the invention can effectively improve the applicability and accuracy of hyperspectral image denoising.

Description

A method for removing hyperspectral hybrid noise based on double total variation

technical field

The patent of the present invention belongs to the field of remote sensing image processing, and more specifically, relates to a method for removing hyperspectral hybrid noise based on double total variation.

Background technique

Hyperspectral remote sensing technology effectively combines advanced spectral technology and traditional two-dimensional image imaging technology. Hyperspectral images have rich spectral information, and such images with very high spectral resolution have wide application value in practice. During the imaging process, due to electromagnetic interference, imaging environment and transmission process, the collected hyperspectral images will be polluted by various mixed noises, including Gaussian noise, impulse noise, stripe noise and dead line noise. Mixing noise in images reduces the image quality and affects the accuracy and reliability of further applications, such as classification, unmixing, and object detection of hyperspectral images. Therefore, hybrid noise removal has important research value in hyperspectral image processing.

At present, many denoising methods have been proposed at home and abroad, and good results have been achieved. The current hyperspectral denoising methods mainly include three types: spatial domain denoising method, spectral domain denoising method and model-based optimization denoising method. In recent years, sparse and low-rank models have been widely used in hyperspectral image denoising, decomposing and estimating the tensor rank from different angles, and the effects are not the same. More appropriate regularization constraints should be selected to retain more useful and detailed information to further improve the effect of hyperspectral image denoising.

The remote sensing image is regarded as a tensor data for modeling and expression, and then by constraining the variable to be determined, it can satisfy a certain tensor property and realize the effective estimation of information restoration. The two core problems in solving remote sensing image information restoration are establishing an effective tensor optimization model and designing an efficient algorithm for solving the model. Therefore, designing an efficient solution algorithm to find a satisfactory solution to the model is an important task at present.

SUMMARY OF THE INVENTION

In order to overcome the above problems existing in the prior art, the patent of the present invention proposes a method for removing hyperspectral hybrid noise based on double total variation.

The present invention is achieved through the following technical solutions:

A hyperspectral hybrid noise removal method based on double total variation. The hyperspectral image is a three-dimensional tensor data, that is, two-dimensional spatial dimension and one-dimensional spectral dimension. Tucker decomposition is used to describe the similarity of the overall structure of high-dimensional images. The visual saliency map acts as a weighting factor on the original image, and then is added to the original image to enhance the saliency features. The spatial segmentation smoothness prior of the image is another effective regular term for mixed noise removal. The SSTV regular term is selected to describe the high The global sparse structure of spectral images, while designing the reweighted L ₁ norm to constrain the local sparse structure of the spatial direction difference image of the principal component map of hyperspectral images, introducing double constraints, selecting more sparse global features and local features, L ₁ sparse regular term It is used to isolate sparse noise, and the F-norm characterizes the Gaussian noise of the image. The established model adopts the ALM algorithm framework to optimize the proposed model, and converts the complex problem into multiple simple sub-problems for alternate iterative solutions. Spectral image, the steps of this method are as follows:

Step 1: Obtain a hyperspectral image Y containing mixed noise;

Step 2: Select the SSTV regular term to describe the global sparse structure of the hyperspectral image, and design the reweighted L ₁ norm to constrain the local sparse structure of the spatial direction difference image of the principal component map of the hyperspectral image, introduce double constraints, and select a more sparse global Features and local features, L ₁ sparse regular term is used to isolate sparse noise, F-norm describes the Gaussian noise of the image, introduces double constraints of sparse structure, obtains more accurate local features and global features, better remove noise, and restore hyperspectral Image, combined with tensor decomposition, build a hyperspectral image denoising model based on double total variation;

Step 3: Introduce auxiliary variables Z, R ₁ and R ₂ to optimize the step 2 hyperspectral image denoising model based on double total variation. The introduction of auxiliary variables makes the problem of step 2 based on the hyperspectral image denoising model based on double total variation Separable, the separation sub-problem is easier to solve, where X=Z, _Dw X=R ₁ , DH=R ₂ , Z is the hyperspectral image equal to X, R ₁ is the differential image of X, and R ₂ is the principal component difference image of the graph;

Step 4: Use the ALM algorithm to solve the optimized hyperspectral image denoising model based on double total variation, and obtain a denoised hyperspectral image.

Further technical solutions include:

The specific process of step 2 is as follows:

The mathematical model for establishing a hyperspectral image denoising model based on prior information is:

Y=X+S+N;

是实数空间，I₁、I₂分别是高光谱图像的长和宽，I₃为光谱的数量，X表示干净的高光谱图像，S表示稀疏噪声，N表示高斯噪声，X、S、N与Y大小相同；高光谱遥感影像混合噪声去除的目标是从观测的Y中复原出真实的影像X，这是一个病态反问题；where Y represents the noisy hyperspectral image,

is a real number space, I ₁ , I ₂ are the length and width of the hyperspectral image, respectively, I ₃ is the number of spectra, X represents a clean hyperspectral image, S represents sparse noise, N represents Gaussian noise, X, S, N and The size of Y is the same; the goal of removing the mixed noise of hyperspectral remote sensing images is to restore the real image X from the observed Y, which is an ill-conditioned inverse problem;

Construct a hyperspectral image denoising model based on double total variation:

Among them, D _w (·)=[w ₁ ×D _x (·); w ₂ ×D _y (·); w ₃ D _z (·)] is an anisotropic space spectral difference operator, D(·)= [D _x ( ); _Dy ( )] is the principal component spatial difference operator, D _x is the difference operator in the horizontal direction of space, _{Dy is the difference operator in the vertical direction of space, and D z is} the difference operator in the spectral direction Difference operator, G is the weighting factor, λ ₁ , λ ₂ , λ ₃ are non-negative regularization parameters used to balance the weights between items, ε is the variance of Gaussian noise, || || ₁ represents the tensor L ₁ norm, || || _F represents the Frobenius norm of the tensor;

Set the sparse coefficient λ ₁ =λ ₂ =1, λ ₃ =10, and the isotropic difference coefficient w=[w ₁ ,w ₂ ,w ₃ ] takes the value [1,1,0.6];

The x, y and z directions respectively represent the spatial horizontal direction, the spatial vertical direction and the spectral direction in the high-dimensional remote sensing data, and D _x , _Dy , and D _z represent the difference operators of the three-dimensional tensor in three different directions, respectively. The value of any position (i, j, p) is defined as follows:

The PCA dimension reduction is performed on the hyperspectral image to obtain the principal component map H, and the edge information map W of the principal component map is obtained through the horizontal and vertical differences, and G is used as a weighting factor, which is updated according to the following formula:

Among them, eps is a very small number. In the implementation process, the matlab built-in function eps() is used to obtain this very small number to avoid the denominator being zero;

PCA dimensionality reduction is realized by matlab embedded function, and the first three principal components are selected;

The specific process of step three is as follows:

Introducing auxiliary variables Z, R ₁ and R ₂ to optimize the model, where X=Z, _Dw X=R ₁ , DH=R ₂ , the proposed hyperspectral image denoising model based on double total variation is equivalently expressed as the following Model:

The specific process of step 4 is as follows:

Based on the ALM solution framework, the augmented Lagrangian function of the hyperspectral image denoising model based on double total variation is written in the following form:

where μ is the penalty coefficient, Λ ₁ , Λ ₂ , Λ ₃ , Λ ₄ are Lagrange multipliers;

Initialize X=Y, Z=R1 ₌ R2=S= _{Λ1 = Λ2} = _Λ3 = _Λ4 =0 tensor, μ=0;

Under the framework of ALM, we can alternately iteratively solve each variable while fixing other variables, and initialize k=0;

Step 4.1, fix other variables, update X ^k+1 :

The sub-problem is solved by the HOOI algorithm. When G ^k+1 and U _i ^k+1 are obtained by the HOOI algorithm, the update of X is substituted into the Tucker decomposition to obtain, namely:

X ^k+1 =C ^k+1 × ₁ U ₁ ^k+1 × ₂ U ₂ ^k+1 × ₃ U ₃ ^k+1 ;

The tensor rank parameter in the HOOI algorithm is set to [120, 120, 10];

Step 4.2, fix other variables, update Z ^k+1 :

The subproblem is a least squares problem, which is obtained by solving the linear equation as follows:

μ(I+D _w ^T D _w )Z=μX+Λ ₂ +D _w ^T (μR ₁ −Λ ₃ );

The above formula is solved using the Fast Fourier Transform:

After obtaining Z ^k+1 of the current iteration process through the above formula, continue to further compensate and update Z ^k+1 ;

Perform SLIC superpixel segmentation on the principal component map H obtained in step 2 to obtain a visual saliency map M, and use the visual saliency map M as a weighting factor to perform element multiplication with each band of the hyperspectral image Z ^k+1 in the current iterative process. , and then added to the hyperspectral image Z ^k+1 of the current iterative process, thereby enhancing the salient features of the hyperspectral image, and obtaining a hyperspectral smooth region with better restoration effect;

Bring Z ^k+1 into Z ^k+1 (:,:,i)+a×M×Z ^k+1 (:,:,i) The calculated result is Z ^{k +1} after compensation update;

Among them, M is the visual saliency map, a is the coefficient of saliency feature enhancement, the SLIC superpixel segmentation is implemented by matlab embedded function, the value range of the visual saliency map coefficient a is set to [0, 1], and the superpixel segmentation parameters are set as The default value, the expected number of superpixels is in the range [300,600];

Step 4.3, fix other variables, update R ₁ ^k+1 :

Using the soft-threshold contraction operator, the closed-form solution of R ₁ is obtained as follows:

Among them, Shrinkage is the soft threshold shrinkage operator, Shrinkage(m,n)=sign(m)max(|m|-n,0), sign(m) is the sign function, when m>0, the value is 1, m = 0, the value is 0, otherwise the value is -1;

Step 4.4, fix other variables, update R ₂ ^k+1 :

Using the soft-threshold contraction operator, the closed _- form solution of R2 is obtained as follows:

Step 4.5, fix other variables, update S ^k+1 :

Using the soft-threshold contraction operator, the closed-form solution of S is obtained as follows:

Step 4.6, fix other variables, update Lagrange multipliers:

当满足Rel≤10^-6，则终止迭代，输出无噪图像；若Rel＞10^-6，则重复步骤四中步骤4.1至步骤4.6，继续交替迭代更新，直至满足迭代终止条件即Rel≤10^-6或迭代次数k达到预设最大迭代次数40时，输出无噪图像X。Step 4.7, judge whether the iteration termination condition is met, the relative change error

When Rel≤10 ^-6 is satisfied, the iteration is terminated, and a noise-free image is output; if Rel>10 ^-6 , steps 4.1 to 4.6 in step 4 are repeated, and the alternate iterative update is continued until the iteration termination condition is satisfied, that is, Rel≤10 ^{- 6} or when the number of iterations k reaches a preset maximum number of iterations of 40, a noise-free image X is output.

The method for removing hyperspectral mixed noise based on double total variation proposed by the present invention selects SSTV regular term to describe the global sparse structure of hyperspectral image, and re-weights L ₁ norm to constrain the local difference image of hyperspectral image principal component map space. Sparse structure, introducing double constraints, selects more sparse global features and local features. The detail information of hyperspectral images plays an important role in practical applications. The traditional SSTV regularization method has poor recovery effect on local details with serious mixed noise pollution in hyperspectral images. The reweighted L ₁ norm constraint is added to select sparser local features. , improve the peak signal-to-noise ratio and structural similarity of the restored image, and then obtain a hyperspectral restored image with better mixed noise removal effect.

Description of drawings

1 is an overall block diagram of a method for removing hyperspectral hybrid noise based on double total variation provided by the present invention;

Figure 2 is the real image of the 160th band of the IndianPines hyperspectral image;

Figure 3 shows the noise image of the 160th band of the IndianPines hyperspectral image;

Figure 4 is the restored image of the 160th band of the IndianPines hyperspectral image;

Detailed ways

The present invention will be further explained below in conjunction with the accompanying drawings.

The invention proposes a hyperspectral hybrid noise removal method based on double total variation. The test image of this method is the Indian Pines hyperspectral image recognized at home and abroad, and the size of the tested Indian Pines hyperspectral image is 145×145×224. The Indian Pines hyperspectral image with mixed noise is the model input Y, which is tested according to the following technical solutions:

Step 1: Obtain a hyperspectral image Y containing mixed noise;

Step 2: Select the SSTV regular term to describe the global sparse structure of the hyperspectral image, and design the reweighted L ₁ norm to constrain the local sparse structure of the spatial direction difference image of the principal component map of the hyperspectral image, introduce double constraints, and select a more sparse global Features and local features, L ₁ sparse regular term is used to isolate sparse noise, F-norm describes the Gaussian noise of the image, introduces double constraints of sparse structure, obtains more accurate local features and global features, better remove noise, and restore hyperspectral Image, combined with tensor decomposition, build a hyperspectral image denoising model based on double total variation; the specific process is:

The mathematical model for establishing a hyperspectral image denoising model based on prior information is:

Y=X+S+N;

是实数空间，I₁、I₂分别是高光谱图像的长和宽，I₃为光谱的数量，X表示干净的高光谱图像，S表示稀疏噪声，N表示高斯噪声，X、S、N与Y大小相同；高光谱遥感影像混合噪声去除的目标是从观测的Y中复原出真实的影像X，这是一个病态反问题；where Y represents the noisy hyperspectral image,

is a real number space, I ₁ , I ₂ are the length and width of the hyperspectral image, respectively, I ₃ is the number of spectra, X represents a clean hyperspectral image, S represents sparse noise, N represents Gaussian noise, X, S, N and The size of Y is the same; the goal of removing the mixed noise of hyperspectral remote sensing images is to restore the real image X from the observed Y, which is an ill-conditioned inverse problem;

Construct a hyperspectral image denoising model based on double total variation:

Among them, D _w (·)=[w ₁ ×D _x (·); w ₂ ×D _y (·); w ₃ D _z (·)] is an anisotropic space spectral difference operator, D(·)= [D _x ( ); _Dy ( )] is the principal component spatial difference operator, D _x is the difference operator in the horizontal direction of space, _{Dy is the difference operator in the vertical direction of space, and D z is} the difference operator in the spectral direction Difference operator, G is the weighting factor, λ ₁ , λ ₂ , λ ₃ are non-negative regularization parameters used to balance the weights between items, ε is the variance of Gaussian noise, || || ₁ represents the tensor L ₁ norm, || || _F represents the Frobenius norm of the tensor;

Set the sparse coefficient λ ₁ =λ ₂ =1, λ ₃ =10, and the isotropic difference coefficient w=[w ₁ ,w ₂ ,w ₃ ] takes the value [1,1,0.6];

The x, y and z directions respectively represent the spatial horizontal, spatial vertical and spectral directions in high-dimensional remote sensing data, and D _x , _Dy , and D _z represent the difference operators of three-dimensional tensors in three different directions and at any position. The value of (i,j,p) is defined as follows:

The PCA dimension reduction is performed on the hyperspectral image to obtain the principal component map H, and the edge information map W of the principal component map is obtained through the horizontal and vertical differences, and G is used as a weighting factor, which is updated according to the following formula:

Among them, eps is a very small number. In the implementation process, the matlab built-in function eps() is used to obtain this very small number to avoid the denominator being zero;

PCA dimensionality reduction is realized by matlab embedded function, and the first three principal components are selected;

Step 3: Introduce auxiliary variables Z, R ₁ and R ₂ to optimize the step 2 hyperspectral image denoising model based on double total variation. The introduction of auxiliary variables makes the problem of step 2 based on the hyperspectral image denoising model based on double total variation Separable, the separation sub-problem is easier to solve, where X=Z, _Dw X=R ₁ , DH=R ₂ , Z is the hyperspectral image equal to X, R ₁ is the differential image of X, and R ₂ is the principal component The differential image of the graph; the specific process is:

Introducing auxiliary variables Z, R ₁ and R ₂ to optimize the model, where X=Z, _Dw X=R ₁ , DH=R ₂ , the proposed hyperspectral image denoising model based on double total variation is equivalently expressed as the following Model:

Step 4: Use the ALM algorithm to solve the optimized hyperspectral image denoising model based on double total variation, and obtain the denoised hyperspectral image; the specific process is as follows:

Based on the ALM solution framework, the augmented Lagrangian function of the hyperspectral image denoising model based on double total variation is written in the following form:

where μ is the penalty coefficient, Λ ₁ , Λ ₂ , Λ ₃ , Λ ₄ are Lagrange multipliers;

Initialize X=Y, Z=R1 ₌ R2=S= _{Λ1 = Λ2} = _Λ3 = _Λ4 =0 tensor, μ=0;

Under the framework of ALM, we can alternately iteratively solve each variable while fixing other variables, and initialize k=0;

Step 4.1, fix other variables, update X ^k+1 :

The sub-problem is solved by the HOOI algorithm. When G ^k+1 and U _i ^k+1 are obtained by the HOOI algorithm, the update of X is substituted into the Tucker decomposition to obtain, namely:

X ^k+1 =C ^k+1 × ₁ U ₁ ^k+1 × ₂ U ₂ ^k+1 × ₃ U ₃ ^k+1 ;

The tensor rank parameter in the HOOI algorithm is set to [120, 120, 10];

Step 4.2, fix other variables, update Z ^k+1 :

The subproblem is a least squares problem, which is obtained by solving the linear equation as follows:

μ(I+D _w ^T D _w )Z=μX+Λ ₂ +D _w ^T (μR ₁ −Λ ₃ );

The above formula is solved using the Fast Fourier Transform:

After obtaining Z ^k+1 of the current iteration process through the above formula, continue to further compensate and update Z ^k+1 ;

Perform SLIC superpixel segmentation on the principal component map H obtained in step 2 to obtain a visual saliency map M, and use the visual saliency map M as a weighting factor to perform element multiplication with each band of the hyperspectral image Z ^k+1 in the current iterative process. , and then added to the hyperspectral image Z ^k+1 of the current iterative process, thereby enhancing the salient features of the hyperspectral image, and obtaining a hyperspectral smooth region with better restoration effect;

Bring Z ^k+1 into Z ^k+1 (:,:,i)+a×M×Z ^k+1 (:,:,i) The calculated result is Z ^{k +1} after compensation update;

Among them, M is the visual saliency map, a is the coefficient of saliency feature enhancement, SLIC superpixel segmentation is implemented by matlab embedded function, the value range of the visual saliency map coefficient a is set as [0, 1], and the superpixel segmentation parameters are set is the default value, and the expected number of superpixels ranges from [300,600];

Step 4.3, fix other variables, update R ₁ ^k+1 :

The sub-problem can be operated with the soft threshold shrinkage operator, and the closed-form solution of R ₁ can be obtained as follows:

Among them, Shrinkage is the soft threshold shrinkage operator, Shrinkage(m,n)=sign(m)max(|m|-n,0), sign(m) is the sign function, when m>0, the value is 1, m = 0, the value is 0, otherwise the value is -1;

Step 4.4, fix other variables, update R ₂ ^k+1 :

The subproblem can be operated with the soft threshold shrinkage operator, and the closed-form solution of R ₂ is obtained as follows:

Step 4.5, fix other variables, update S ^k+1 :

The subproblem can be operated with the soft threshold shrinkage operator, and the closed-form solution of S is obtained as follows:

Step 4.6, fix other variables, update Lagrange multipliers:

当满足Rel≤10^-6，则终止迭代，输出无噪图像；若Rel＞10^-6，则重复步骤四中步骤4.1至步骤4.6，继续交替迭代更新，直至满足迭代终止条件即Rel≤10^-6或迭代次数k达到预设最大迭代次数40时，输出无噪图像X。Step 4.7, judge whether the iteration termination condition is met, the relative change error

When Rel≤10 ^-6 is satisfied, the iteration is terminated, and a noise-free image is output; if Rel>10 ^-6 , steps 4.1 to 4.6 in step 4 are repeated, and the alternate iterative update is continued until the iteration termination condition is satisfied, that is, Rel≤10 ^{- 6} or when the number of iterations k reaches a preset maximum number of iterations of 40, a noise-free image X is output.

In particular, each sub-problem in step 4 will be solved alternately and iteratively based on the ALM framework, until Rel≤10 ^-6 or the iteration number reaches the preset maximum number of iterations 40, the iteration is terminated, and the output X is the noise-free IndianPines hyperspectral image, To achieve the purpose of removing the mixed noise of hyperspectral images.

The results of noise removal of the present invention are evaluated from different perspectives, including: visual comparison and quantitative comparison. The indicators of quantitative comparison, the average peak signal-to-noise ratio (MPSNR) and the average structural similarity (MSSIM), the larger the value, the closer the denoised image is to the reference image, and the better the denoising effect.

A method for removing hyperspectral mixed noise based on double total variation proposed by the present invention, the experimental platform is realized based on MATLAB software, and the PCA principal component analysis and SLIC algorithm both use MATLAB embedded functions. The results of the case simulation experiment are shown in the figure. Figure 2 is the real image of the 160th band of the Indian Pines hyperspectral image, Figure 3 is the image of the 160th band of the Indian Pines hyperspectral image contaminated by noise, and Figure 4 is the Indian Pines hyperspectral image. The 160th band of the image is restored by this method. By analyzing Figure 2, Figure 3 and Figure 4, it can be seen that the image details restored by this method are very clear. The MPSNR of the Indian Pines hyperspectral noise image is 13.7389dB, and the MSSIM is 0.2025. The MPSNR of the Indian Pines hyperspectral restored image restored by this method is increased to 38.8024dB, and the MSSIM is increased to 0.9848, which proves the effectiveness and feasibility of this method.

Claims (4)

1. A hyperspectral hybrid noise removal method based on double total variation. Hyperspectral image is a three-dimensional tensor data, that is, two-dimensional spatial dimension and one-dimensional spectral dimension. Tucker decomposition is used to describe the similarity of the overall structure of high-dimensional images. , the visual saliency map is used as a weighting factor on the original image, and then added to the original image to enhance the saliency features. The spatial slice smoothness prior of the image is another effective regular term for mixed noise removal. Select the SSTV regular term to Describe the global sparse structure of hyperspectral images, and design the reweighted L ₁ norm to constrain the local sparse structure of spatial direction difference images of the principal component map of hyperspectral images, introduce double constraints, select more sparse global features and local features, and L ₁ is sparse The regular term is used to isolate the sparse noise, and the F-norm characterizes the Gaussian noise of the image. The established model adopts the ALM algorithm framework to optimize the proposed model, and converts the complex problem into multiple simple sub-problems to solve alternately and iteratively to obtain denoising. The post-hyperspectral image is characterized in that the steps of the method are as follows: Step 1: Obtain a hyperspectral image Y containing mixed noise; Step 2: Select the SSTV regular term to describe the global sparse structure of the hyperspectral image, and design the reweighted L ₁ norm to constrain the local sparse structure of the spatial direction difference image of the principal component map of the hyperspectral image, introduce double constraints, and select a more sparse global Features and local features, L ₁ sparse regular term is used to isolate sparse noise, F-norm describes the Gaussian noise of the image, introduces double constraints of sparse structure, obtains more accurate local features and global features, better remove noise, and restore hyperspectral Image, combined with tensor decomposition, build a hyperspectral image denoising model based on double total variation; Step 3: Introduce auxiliary variables Z, R ₁ and R ₂ to optimize the step 2 hyperspectral image denoising model based on double total variation. The introduction of auxiliary variables makes the problem of step 2 based on the hyperspectral image denoising model based on double total variation Separable, the separation sub-problem is easier to solve, where X=Z, _Dw X=R ₁ , DH=R ₂ , Z is the hyperspectral image equal to X, R ₁ is the differential image of X, and R ₂ is the principal component difference image of the graph; Step 4: Use the ALM algorithm to solve the optimized hyperspectral image denoising model based on double total variation, and obtain a denoised hyperspectral image. 2. a kind of hyperspectral hybrid noise removal method based on double total variation according to claim 1, is characterized in that, the concrete process of step 2 is: The mathematical model for establishing a hyperspectral image denoising model based on prior information is: Y=X+S+N; where Y represents the noisy hyperspectral image,

is a real number space, I ₁ , I ₂ are the length and width of the hyperspectral image, respectively, I ₃ is the number of spectra, X represents a clean hyperspectral image, S represents sparse noise, N represents Gaussian noise, X, S, N and The size of Y is the same; the goal of removing the mixed noise of hyperspectral remote sensing images is to restore the real image X from the observed Y, which is an ill-conditioned inverse problem; Construct a hyperspectral image denoising model based on double total variation:

Among them, D _w (·)=[w ₁ ×D _x (·); w ₂ ×D _y (·); w ₃ D _z (·)] is an anisotropic space spectral difference operator, D(·)= [D _x ( ); _Dy ( )] is the principal component spatial difference operator, D _x is the difference operator in the horizontal direction of space, _{Dy is the difference operator in the vertical direction of space, and D z is} the difference operator in the spectral direction Difference operator, G is the weighting factor, λ ₁ , λ ₂ , λ ₃ are non-negative regularization parameters used to balance the weights between items, ε is the variance of Gaussian noise, || || ₁ represents the tensor L ₁ norm, || || _F represents the Frobenius norm of the tensor; Set the sparse coefficient λ ₁ =λ ₂ =1, λ ₃ =10, and the isotropic difference coefficient w=[w ₁ ,w ₂ ,w ₃ ] takes the value [1,1,0.6]; The x, y and z directions respectively represent the spatial horizontal direction, the spatial vertical direction and the spectral direction in the high-dimensional remote sensing data, and D _x , _Dy , and D _z represent the difference operators of the three-dimensional tensor in three different directions, respectively. The value of any position (i, j, p) is defined as follows:

The PCA dimension reduction is performed on the hyperspectral image to obtain the principal component map H, and the edge information map W of the principal component map is obtained through the horizontal and vertical differences, and G is used as a weighting factor, which is updated according to the following formula:

Among them, eps is a very small number. In the implementation process, the matlab built-in function eps() is used to obtain this very small number to avoid the denominator being zero. PCA dimensionality reduction is realized by matlab embedded function, and the first three principal components are selected. 3. a kind of hyperspectral hybrid noise removal method based on double total variation as described in claim 2, is characterized in that, the concrete process of step 3 is: Introducing auxiliary variables Z, R ₁ and R ₂ to optimize the model, where X=Z, _Dw X=R ₁ , DH=R ₂ , the proposed hyperspectral image denoising model based on double total variation is equivalently expressed as the following Model:

. 4. a kind of hyperspectral hybrid noise removal method based on double total variation as described in claim 3, is characterized in that, the concrete process of step 4 is: Based on the ALM solution framework, the augmented Lagrangian function of the hyperspectral image denoising model based on double total variation is written in the following form:

where μ is the penalty coefficient, Λ ₁ , Λ ₂ , Λ ₃ , Λ ₄ are Lagrange multipliers; Initialize X=Y, Z=R1 ₌ R2=S= _{Λ1 = Λ2} = _Λ3 = _Λ4 =0 tensor, μ=0; Under the framework of ALM, we can alternately iteratively solve each variable while fixing other variables, and initialize k=0; Step 4.1, fix other variables, update X ^k+1 :

The sub-problem is solved by the HOOI algorithm. When G ^k+1 and U _i ^k+1 are obtained by the HOOI algorithm, the update of X is substituted into the Tucker decomposition to obtain, namely: X ^k+1 =C ^k+1 × ₁ U ₁ ^k+1 × ₂ U ₂ ^k+1 × ₃ U ₃ ^k+1 ; The tensor rank parameter in the HOOI algorithm is set to [120, 120, 10]; Step 4.2, fix other variables, update Z ^k+1 :

The subproblem is a least squares problem, which is obtained by solving the linear equation as follows: μ(I+D _w ^T D _w )Z=μX+Λ ₂ +D _w ^T (μR ₁ −Λ ₃ ); The above formula is solved using the Fast Fourier Transform:

After obtaining Z ^k+1 of the current iteration process through the above formula, continue to further compensate and update Z ^k+1 ; Perform SLIC superpixel segmentation on the principal component map H obtained in step 2 to obtain a visual saliency map M, and use the visual saliency map M as a weighting factor to perform element multiplication with each band of the hyperspectral image Z ^k+1 in the current iterative process. , and then added to the hyperspectral image Z ^k+1 of the current iterative process, thereby enhancing the salient features of the hyperspectral image, and obtaining a hyperspectral smooth region with better restoration effect; Bring Z ^k+1 into Z ^k+1 (:,:,i)+a×M×Z ^k+1 (:,:,i) The calculated result is Z ^k+1 after compensation update; Among them, M is the visual saliency map, a is the coefficient of saliency feature enhancement, the SLIC superpixel segmentation is implemented by matlab embedded function, the value range of the visual saliency map coefficient a is set to [0, 1], and the superpixel segmentation parameters are set as The default value, the expected number of superpixels is in the range [300,600]; Step 4.3, fix other variables, update R ₁ ^k+1 :

Using the soft-threshold contraction operator, the closed-form solution of R ₁ is obtained as follows:

Among them, Shrinkage is the soft threshold shrinkage operator, Shrinkage(m,n)=sign(m)max(|m|-n,0), sign(m) is the sign function, when m>0, the value is 1, m = 0, the value is 0, otherwise the value is -1; Step 4.4, fix other variables, update R ₂ ^k+1 :

Using the soft-threshold contraction operator, the closed _- form solution of R2 is obtained as follows:

Step 4.5, fix other variables, update S ^k+1 :

Using the soft-threshold contraction operator, the closed-form solution of S is obtained as follows:

Step 4.6, fix other variables, update Lagrange multipliers:

Step 4.7, judge whether the iteration termination condition is met, the relative change error

When Rel≤10 ^-6 is satisfied, the iteration is terminated, and a noise-free image is output; if Rel>10 ^-6 , steps 4.1 to 4.6 in step 4 are repeated, and the alternate iterative update is continued until the iteration termination condition is satisfied, that is, Rel≤10 ^{- 6} or when the number of iterations k reaches a preset maximum number of iterations of 40, a noise-free image X is output.

Images