A Hyperspectral Image Denoising Method Based on Gaussian-Wishart Prior

Abstract

The present invention discloses a hyperspectral image denoising method based on Gaussian-Wishart prior, which includes the following steps: 1) Input the noisy hyperspectral image. W and H are the length and width of the spatial dimension of the hyperspectral image respectively, and S is the number of spectral bands; 2) Extract overlapping full-band image patches with a size of dw×dH×S, and use K-means++ to cluster the image patches into groups and input them into the non-parametric Bayesian CP decomposition model; 3) In the Bayesian model, derive the posterior probability formula of the parameters through the relationship between the conjugate prior and the likelihood function, and update the parameters to be estimated sequentially by Gibbs sampling; 4) Reconstruct the image group from the weight coefficients and factor matrices of CP decomposition; 5) Restore the image patches, and perform mean smoothing at the overlapping parts of the image patches to reconstruct the image of this iteration; 6) After balancing the results of this iteration and the original input image through the regularization coefficient, use them as the input image for the next iteration. After multiple iterations of steps 1-6, output the reconstruction result of the last iteration. The present invention does not require the noise variance as input and has good denoising effects and applicability.

Description

Technical Field

[0001] The present invention relates to the field of image processing, and in particular to a hyperspectral image denoising method based on Gaussian-Wishart prior, which can be used to remove various types of noises such as Gaussian noise, salt-and-pepper noise, and mixed noise brought during the acquisition process of hyperspectral images. Background Art

[0002] Hyperspectral images have been continuously developed and widely used in many monitoring tasks due to their rich available spectral information. Typical applications include target detection, environmental monitoring, ground object recognition, medical tomography analysis, etc. However, due to the large amount of hyperspectral image data and the influence of factors such as the environment and acquisition equipment during imaging, the image is inevitably damaged by noise. Denoising is a research hotspot in the field of hyperspectral images. As a key preprocessing process, it provides quality assurance for subsequent applications that need to extract spatial-spectral information, such as spectral unmixing, classification, and target detection.

[0003] Regarding the problem of hyperspectral image denoising, many methods have been proposed by domestic and foreign researchers. Classic denoising algorithms include: wavelet transform algorithm, singular value decomposition algorithm, and optimized reconstruction algorithm, etc. These algorithms vectorize the high-dimensional data and adopt the processing methods of two-dimensional images, which will not only destroy the internal structure of the high-dimensional data, lose the original information of the data, but also lose the advantage of rich spectral information of the high-dimensional data. In order to better utilize information such as spectral correlation, sparse and low-rank models are applied to hyperspectral image denoising, especially the method combining tensor decomposition and Bayesian. Through sparse and low-rank priors, it not only improves the representation ability of the three-dimensional structure of hyperspectral images, but also takes into account the spatial local self-similarity and spectral correlation of hyperspectral images.

[0004] Tensor decomposition mainly includes two decomposition forms: Tucker decomposition and CANDECOMP / PARAFAC (CP) decomposition. CP decomposition decomposes a tensor into the sum of multiple rank-one tensors, and the number of rank-one tensors is the CP rank. Tucker decomposition is a high-order SVD method that decomposes a tensor into a core tensor and factor matrices in each mode. For an N-dimensional tensor, there are N factor matrices, and the tensor rank includes all the ranks in these N modes. Solving the tensor rank for both methods is an NP-Hard problem. The application of non-parametric Bayesian methods can estimate the tensor rank while learning the factor matrices. However, since a large input needs to be artificially estimated as the rank, there is a correlation between the factors. Traditional non-parametric Bayesian methods assume that the factors are uncorrelated, and the flexibility and denoising performance of the algorithms need to be improved. In addition, in actual situations, the noise is often unknown, and existing tensor low-rank methods often require the noise variance as an input, resulting in poor adaptability in practical applications. Summary of the Invention

[0005] In order to overcome the deficiencies of the prior art and improve the denoising performance, the present invention proposes a hyperspectral image denoising method based on Gaussian-Wishart prior. First, the image is segmented into blocks, and the K-means++ clustering method is used to divide the blocks into multiple groups, effectively utilizing the non-local self-similarity in the spatial dimension of the image. Then, the non-parametric Bayesian CP decomposition model is used to decompose the tensors of each group, the Gaussian-Wishart prior is adopted to explore the correlation between factors, and the Gibbs sampling method is used for parameter inference. Finally, the regularization parameter is used to balance the input at each time, and the denoising result is obtained through multiple iterations.

[0006] The technical solution adopted by the present invention to solve its technical problems is: a hyperspectral image denoising method based on Gaussian-Wishart prior, and the specific operation steps are as follows:

[0007] Step 1) Obtain a hyperspectral image contaminated by noise where represents the noisy image, represents the noise, which can be Gaussian noise, salt-and-pepper noise, and mixed noise, etc., is the clean image to be restored. W and H are the length and width of the spatial dimension of the hyperspectral image respectively, and S is the number of spectral bands.

[0008] Step 2) Extract overlapping full-band image blocks of size d w ×d H ×S, where d w , d H is the size of the image block. The K-means++ clustering is used to divide the image blocks into C groups, and then the image blocks in the groups are transformed into two-dimensional data of (d w d H )×S, and the groups are input into the non-parametric Bayesian CP decomposition model in turn.

[0009] Step 3) In the non-parametric Bayesian CP decomposition model, each input tensor group of size (d w d H )×S×L is represented as Its CP decomposition can be expressed as:

[0010]

[0011] represents the clean tensor group to be reconstructed, λ r represents the weight coefficient of the factor, R represents the CP rank, and 1≤r≤R. represents the factor vector, and the factor matrices

[0012] In Bayesian analysis, it is necessary to specify the prior distribution and discuss the probability distribution of parameters. Assume that the noise ε follows a Gaussian distribution with precision τ ε , and τ ε follows a gamma prior; the prior of the weight coefficient λ r of CP decomposition is set to the MultiplicativeGamma Process (MGP) prior; assume that the factor row vectors of CP decomposition follow a multivariate Gaussian distribution, and the mean and variance follow a Gaussian-Wishart prior. The posterior probability formula of the parameters is derived through the relationship between the conjugate prior and the likelihood function, and the parameters to be estimated, namely the above CP decomposition coefficients, factor matrices, and noise precision, etc., are updated sequentially by Gibbs sampling.

[0013] Step 4) Reconstruct the image group from the CP decomposition weight coefficients and factor matrices obtained by Gibbs sampling: This is the reconstructed clean image group.

[0014] Step 5) After all image groups are reconstructed, restore the image patches and mean-smooth the overlapping parts of the image patches to reconstruct the image recovered this time.

[0015] Step 6) After balancing the result reconstructed this time and the original input image with the regularization coefficient α n , use it as the input image for the next iteration.

[0016]

[0017] Among them represents the original noisy hyperspectral image, represents the result this time, f represents the non-parametric Bayesian CP decomposition model, represents the result of the next processing. After N iterations of steps 1-6, the reconstruction result of the Nth time is output as the final hyperspectral image denoising result of this method.

[0018] Compared with the prior art, the present invention adopts the above technical solutions and has the following beneficial effects:

[0019] This technical solution is mainly used to solve the problem of noise interference in hyperspectral images, and can remove various noises such as Gaussian noise, salt-and-pepper noise, and mixed noise contained in hyperspectral images. K-means++ is used to cluster image patches, effectively utilizing the non-local similarity of hyperspectral images in the spatial dimension; CP decomposition is adopted to better represent the internal structure of high-dimensional data; a non-parametric Bayesian method based on Gaussian-Wishart prior is used, which can not only adaptively infer the CP rank and noise intensity, but also flexibly explore the correlation between factors; the input image is balanced through a regularization coefficient, and the texture and edge information in the image are deeply extracted, further improving the denoising ability for noise. The present invention does not require the noise variance as input and has good denoising effects and applicability. Description of the Drawings

[0020] Figure 1 is a flowchart of the method of the present invention.

[0021] Figure 2 is a hyperspectral image containing Gaussian noise and a hyperspectral image after denoising using the present invention.

[0022] Figure 3 is a hyperspectral image containing salt-and-pepper noise and a hyperspectral image after denoising using the present invention.

[0023] Figure 4 is a hyperspectral image containing Gaussian-salt-and-pepper mixed noise and a hyperspectral image after denoising using the present invention. Detailed Embodiments

[0024] The technical solution of the present invention will be further described below with reference to the drawings.

[0025] A hyperspectral image denoising method based on Gaussian-Wishart prior of the present invention includes the following steps:

[0026] Step 1) Noise Observation Model

[0027] Obtain a hyperspectral image contaminated by noise Where represents the noise image, represents the noise, which can be Gaussian noise, salt-and-pepper noise, and mixed noise, etc., is the clean image to be restored. W and H are the length and width of the spatial dimension of the hyperspectral image respectively, and S is the number of spectral bands.

[0028] Step 2) Clustering Method

[0029] Due to the non-local similarity in the spatial dimension of hyperspectral images, K-means++ clustering can be adopted. K-means++ optimizes the distribution of the initial clustering centers based on the original K-means clustering, improving the clustering effect. First, extract overlapping full-band image patches of size d w ×d H ×S from the hyperspectral image, where d w , d H is the size of the image patch. Set the scanning step size as p, 1 ≤ p < min(d w , d H ), and segment the hyperspectral image into overlapping image patches. The number of patches Blocks = ((W - d w ) / p + 1) * ((H - d h ) / p + 1). In this method, d w = d H = 6, p = 2. Use K-means++ clustering to divide the image patches into C (C = Blocks / 100) groups. To reduce the data dimension and computational complexity, transform the image patches in the group into two-dimensional data of (d w d H ) × S. Assume that a group contains L image patches, then the size of the group is (d w d H ) × S × L. Input all the groups into the non-parametric Bayesian CP decomposition model in sequence.

[0030] Step 3) Core model - Non-parametric Bayesian CP decomposition model based on Gaussian-Wishart prior

[0031] Input all the groups in step 2) in a loop. Each input three-dimensional tensor group is represented as Its noise model can be expressed as:

[0032]

[0033] represents the clean tensor group to be reconstructed, and ε represents the noise contained in the tensor group .

[0034] The CP decomposition for reconstructing the clean tensor group can be expressed as:

[0035]

[0036] Where λ r represents the weight coefficient of the factor, R represents the CP rank, initialized as S, 1 ≤ r ≤ R. represents the factor vector, and the factor matrix That is, for a three-dimensional tensor group, there are three factor matrices, and each factor matrix contains R factor column vectors.

[0037] In Bayesian analysis, it is necessary to specify the prior distribution and then discuss the probability distribution of the parameters.

[0038] Step 3-1) Assume that the model noise ε follows a Gaussian distribution with precision τ ε and τ ε follows a gamma prior. That is According to formulas (1) and (2), the likelihood function of the model can be expressed as:

[0039]

[0040] where F represents the Frobenius norm, which is defined as the square root of the sum of the squares of the absolute values of the elements in the tensor.

[0041] For each element y i in it and its noise model, i = (i1 i2...i K ) 1 ≤ i k ≤ I K , formula 2 can be rewritten as:

[0042]

[0043] is the element of the factor matrix U k at (i k , r).

[0044] Correspondingly, formula (3) can be rewritten as:

[0045]

[0046] Step 3-2) Set the Multiplicative Gamma Process (MGP) multiplicative gamma process prior for the weight coefficient λ r in the model:

[0047]

[0048] λ r follows a Gaussian distribution with mean 0 and precision τ r , δ l follows a gamma distribution, and a c is the hyperparameter of the gamma distribution, which is set to 3 in the experiment. Use Gibbs sampling to update the parameter λ r to be estimated, and its update process is as follows:

[0049]

[0050]

[0051] The update process of variance and mean is as follows:

[0052] Step 3-3) Considering the correlation between factor matrices, for the row vectors of the factor matrix Set a multivariate Gaussian prior, and for the unknown parameters μ k , Λ k Set a Gaussian-Wishart prior:

[0053]

[0054] The most important role of the Wishart distribution is to describe the covariance of a multivariate normal distribution sample, as the conjugate prior of a multivariate normal distribution with an unknown mean and precision matrix (the inverse of the covariance matrix):

[0055]

[0056] Derive the posterior probability formulas of the factor row vector sum and parameters such as mean and variance through the relationship between the conjugate prior and the likelihood function. Gibbs sampling updates the parameters to be estimated according to the posterior formula in turn, and the process is as follows:

[0057] Row vector Posterior probability distribution:

[0058]

[0059]

[0060]

[0061] Parameter μ k , Λ k Posterior probability distribution:

[0062]

[0063]

[0064]

[0065]

[0066]

[0067] Among them,

[0068] Step 3-4) To infer the noise precision τ ε , set a flexible gamma prior τ ε~Gamma(a0, b0), and usually the initial values a0 = b0 = 1 are set.

[0069] Parameter τ ε Posterior probability distribution of:

[0070]

[0071]

[0072] Step 4) includes the following steps:

[0073] Reconstruct the image group from the CP decomposition weight coefficients and factor matrices obtained by Gibbs sampling in step 3: That is the reconstructed clean image group.

[0074] Step 5) includes the following steps:

[0075] After all the image groups are reconstructed, each two-dimensional data block in the group is restored to a three-dimensional tensor and restored to the position of the original image. For the overlapping parts of the blocks, averaging is used for smoothing processing, and finally the hyperspectral image of this restoration is reconstructed. Step 6) includes the following steps:

[0076] After balancing the result of this reconstruction and the original input image with the regularization coefficient α n as the input image for the next iteration, which is expressed as follows:

[0077]

[0078] where represents the original noisy hyperspectral image, represents the result of this time, f represents the non-parametric Bayesian CP decomposition model, represents the result of the next processing. After N iterations of steps 1 - 6, the reconstruction result of the Nth time is output as the final hyperspectral image denoising result of this method. In the method experiment, α n = 0.2, N = 3.

[0079] Some of the images listed in the attached drawings are the denoising results under different noises. The content described in the embodiments of this specification is only an enumeration of the implementation forms of the inventive concept. The protection scope of the present invention should not be regarded as limited to the specific forms stated in the embodiments. The protection scope of the present invention also extends to equivalent technical means that can be thought of by those skilled in the art according to the inventive concept of the present invention.

Claims

1. A hyperspectral image denoising method based on non-parametric Bayesian CP decomposition with Gaussian-Wishart prior, comprising the following steps: Step 1) Obtain the original noisy hyperspectral image where represents the noisy image, represents the noise, including Gaussian noise, salt-and-pepper noise, and mixed noise, is the clean image to be restored, W and H are the length and width of the spatial dimension of the hyperspectral image respectively, and S is the number of spectral bands; Step 2) Extract overlapping full-band image patches of size d w ×d h ×S, where d w and d h respectively represent the width and height of the extracted image patches. Use K-means++ clustering to divide the image patches into C groups, and then transform the image patches in the groups into (d w d h )×S two-dimensional data, and input the groups into the non-parametric Bayesian CP decomposition model in sequence; Step 3) In the non-parametric Bayesian CP decomposition model, assume that the noise ε follows a Gaussian distribution with precision τ ε , τ ε follows a Gamma prior, set a multiplicative gamma process prior for the weight coefficient λ of CP decomposition; assume that the factor row vectors of CP decomposition follow a multivariate Gaussian distribution, and the mean and variance follow a Gaussian-Wishart prior. Derive the posterior probability formula of the parameters through the relationship between the conjugate prior and the likelihood function, and update the parameters to be estimated in turn by Gibbs sampling. The parameters to be estimated are the weight coefficient λ r of CP decomposition, the factor matrix U r , and the noise precision τ k ; ε ; Step 4) Reconstruct the image group from the CP decomposition weight coefficients and factor matrices obtained by Gibbs sampling: where R represents the CP rank, 1 ≤ r ≤ R, λ r represents the weight coefficient of the factor, represents the factor vector, and the factor matrix Step 5) After all the image groups are reconstructed, restore the image patches and smooth the overlapping parts of the image patches by averaging to reconstruct the image restored this time; Step 6) After balancing the result of this iteration and the original noisy hyperspectral image with the regularization coefficient a n use it as the input for the next iteration; After multiple iterations of steps 1-6, output the reconstruction result of the last iteration.

2. A hyperspectral image denoising method based on nonparametric Bayesian CP decomposition with Gaussian-Wishart prior as claimed in claim 1, characterized in that: Step 2) comprises the following steps: Input the original noisy hyperspectral image Initialize the size of the sub-blocks as d w , d h , where d w , d h represent the width and height of the extracted image blocks respectively, the scanning step is p, 1 ≤ p < min(d w , d h ), segment the hyperspectral image into overlapping image blocks, the number of blocks is ((W - d w ) / p + 1) × (H - d h ), through the K-means++ clustering method, set the initial number of clusters C, cluster the similar blocks into a group with a four-dimensional tensor of dimension d w× d h ×S×L, where L is the number of blocks in a group, unfold each image block in the group along the spatial dimension into two-dimensional data, and each group becomes a three-dimensional tensor of dimension (d w d h )×S×L, input these C groups of three-dimensional tensors into the next step.

3. A hyperspectral image denoising method based on non-parametric Bayesian CP decomposition with Gaussian-Wishart prior as claimed in claim 1, characterized in that: Step 5) comprises the following steps: After all the image groups are reconstructed, each two-dimensional data block in the group is restored to a three-dimensional tensor and restored to the position of the original image. For the overlapping part of the blocks, smoothing processing is performed by taking the average value, and finally the hyperspectral image restored this time is reconstructed.

4. A hyperspectral image denoising method based on nonparametric Bayesian CP decomposition with Gaussian-Wishart prior as claimed in claim 1, characterized in that: Step 6) comprises the following steps: After balancing the results of this reconstruction and the original noisy hyperspectral image with the regularization coefficient a n it is used as the input image for the next iteration, expressed as follows: Among them represents the result of this reconstruction, and f represents the non-parametric Bayesian CP decomposition model. represents the result of the next processing; after N iterations of steps 1-6, the reconstruction result of the Nth time is output as the final hyperspectral image denoising result.

Images