CN106651778A - Spectral imaging method based on self-adaptive coupling observation and non-linear compressed learning - Google Patents

Abstract

The invention discloses a spectral imaging method based on self-adaptive coupling observation and non-linear compressed learning to mainly solve the problem that the prior art cannot ensure that matrix observation is not related to a dictionary obtained by learning. The method is realized by the following steps: 1, projecting signals in an original space to a feature space; 2, conducting dictionary learning in the feature space through a KPCA method and solving a sparse dictionary; 3, conducting random initialization on an observation matrix, conducting iteration training on a perception matrix formed by a product of the sparse dictionary and the initialized observation matrix, and obtaining an observation matrix that has undergone coupling optimization; 4, realizing non-linear compressed perception spectral imaging through a nucleus compressed perception method; and 5, restoring original signals through a pre-image method. The result of an experiment shows that in a condition of different sampling rates, the spectral imaging method, compared with a conventional method taking a Gaussian random matrix as an observation matrix, has a great reconstitution effect and can be used for obtaining low-cost and high-quality high spectral images.

Description

Spectral imaging method based on adaptive coupling observation and nonlinear compression learning

Technical Field

The invention belongs to the technical field of image processing, and particularly relates to a compression spectrum imaging method which can be used for low-cost and high-quality acquisition of a hyperspectral image.

Background

Compressed sensing is a new sampling theory developed in the signal processing technology field in recent years, and by utilizing the sparsity of signals, accurate recovery of information can be realized under the condition that the sampling rate is far less than the traditional Nyquist sampling rate. The existing methods for compressing sensing are all based on an explicit linear sparse representation model. The linear sparse representation model has the advantages of simplicity, intuition, easiness in understanding, easiness in operation and the like. However, the actual scene information is complex, and it is difficult to obtain sufficiently sparse representation under a linear sparse representation model. If a linear sparse representation model is used, the sparsity of obtaining sparse coefficients is low, which results in more measurements needed to recover the original signal. On the other hand, for the application of the compressive sensing technology in compressive imaging, a better observation system is designed to improve the imaging quality. Most of the existing compression imaging methods are based on random observation, and an observation matrix is selected as a Gaussian random matrix or a Bernoulli random matrix. However, such random observation matrix can be uncorrelated with most orthogonal dictionaries, and meets the condition of limited equidistant characteristics, but the random observation matrix still belongs to a non-adaptive observation matrix, which has universality, but for different signals, because the information of the signals is not fully mined, and the self-adaptability and optimality are not provided, a high-quality compressed imaging result cannot be obtained under a small number of sensors.

Disclosure of Invention

The invention aims to provide a spectral imaging method based on adaptive coupling observation and nonlinear compression learning, aiming at overcoming the defects of the prior art, so that an observation matrix and a learned dictionary have better non-correlation, and the quality of compression imaging is improved under a small amount of observation.

The technical scheme of the invention is that through a nonlinear kernel function, a signal in an original space is projected onto a feature space, and dictionary learning is carried out in the feature space; performing iterative training on a perception matrix formed by the product of a sparse dictionary and an initial observation matrix by a method similar to K-SVD (K-singular value decomposition) to obtain an observation matrix subjected to coupling optimization; the nonlinear compressed sensing spectral imaging is realized by a nuclear compressed sensing method. The method comprises the following implementation steps:

(1) selecting three groups n₁×n₂×n₃Randomly selecting n spectral bands except the 10 th spectral band to construct a training sample matrix Y ═ Y₁,y₂,…,y_j,…,y_n]Wherein n is₁×n₂Is the size of the hyperspectral image, n₃Denotes the number of total bands, y_jRepresenting a column vector drawn by the jth spectral segment, wherein j is 1,2, …, n and n are the number of training samples;

(2) training a dictionary by utilizing a training sample matrix Y, solving a sparse dictionary of the training samples by adopting a Kernel Principal Component Analysis (KPCA) method, and recording the sparse dictionary asn is the number of samples, and K is the number of atoms of the obtained sparse dictionary;

(3) setting the number of measured values as m, and randomly generating a Gaussian random matrix as an initial observation matrixi＝1,2,…,m，φ_iIs to initialize an observation matrix phi₀A row vector of (a);

(4) combining a training sample matrix Y and an initial observation matrix phi according to a sparse dictionary A₀Obtaining a sensing matrix G ═ f (Φ)₀ ^T)^T(Y) A, training a perception matrix G to obtain a final observation matrix phi, wherein f is a nonlinear mapping function, and T represents transposition;

(5) drawing three groups of 10 th spectral band images into column vectors as the measurementTest samples, respectively denoted as e₁,e₂,e₃；

(6) According to the observation matrix phi obtained in the step (4), three test samples e in the step (5) are subjected to₁,e₂,e₃Carrying out nonlinear compression imaging to obtain a measured value m ═ f_k(<Φ,e_i>) Wherein e is_iDenotes the ith group of test specimens, i from 1 to 3, f_kIs a selected kernel function;

(7) computing sparse coefficients using least squaresWhereinRepresents a pseudo-inverse;

(8) reconstructing the original image by a pre-image method according to the sparse coefficient β and the sparse dictionary AWhereinThe restored image for the ith group.

Compared with the prior art, the invention has the following advantages:

1. the adaptive coupling optimization algorithm in the kernel space provided by the invention enables the observation matrix and the sparse dictionary obtained through dictionary learning to have higher non-correlation, and can better meet the limited equidistant characteristic condition of a compressed sensing theory.

2. According to the method, the observation matrix is obtained by performing coupling optimization according to the sparse dictionary obtained by dictionary learning, the information of the signal can be fully mined, and the method has self-adaptability and better reconstruction effect.

Drawings

FIG. 1 is a flow chart of an implementation of the present invention;

FIG. 2 is a comparison graph of the reconstruction effect of the InianPines test image at a sampling rate of 10% using the method of the present invention and a Gaussian random matrix as the observation matrix;

FIG. 3 is a plot of mean square error of a recovered image reconstructed from a test image IndianPines at a sampling rate of 10% using the method of the present invention and a Gaussian random matrix as the observation matrix;

FIG. 4 is a comparison graph of the reconstruction effect of the Moffet of the test image at a sampling rate of 10% using the method of the present invention and a Gaussian random matrix as the observation matrix;

FIG. 5 is a plot of the mean square error of the recovered image reconstructed from the test image Moffet at a sampling rate of 10% using the method of the present invention and a Gaussian random matrix as the observation matrix;

FIG. 6 is a comparison graph of the reconstruction effect of the test image WashtonDC at a sampling rate of 10% using the method of the present invention and a Gaussian random matrix as an observation matrix;

fig. 7 is a mean square error curve of an image restored by reconstructing a test image WashtonDC at a sampling rate of 10% using the method of the present invention and a gaussian random matrix as an observation matrix.

Detailed description of the invention

Referring to fig. 1, the specific implementation steps of the present invention are as follows:

step 1, constructing a training sample matrix.

Obtaining three groups n from an image library₁×n₂×n₃Taking the 10 th spectral band image of each group of hyperspectral images as a test sample, taking other n spectral bands as training samples, and drawing each spectral band image of the training samples intoColumn vectors forming a training sample matrix of size w × n, Y ═ Y₁,y₂,…,y_j,…,y_n]Wherein w ═ n₁×n₂y_jRepresents the column vector into which each spectral fragment is pulled, j ═ 1,2, …, n.

And 2, training the dictionary by using the training sample Y.

In the nonlinear compression imaging, the dictionary training method comprises methods of kernel principal component analysis, kernel K-SVD, kernel independent component analysis and the like, and the time for training the dictionary by the kernel principal component analysis method is relatively short, so that the method adopts a kernel principal component analysis KPCA method to train a training sample to obtain a sparse dictionary, and comprises the following steps: .

2a) Selecting the kernel function as a polynomial kernel function k (x, y) ═ f_k(<x,y>)＝(<x,y>+c)^dWherein x and y are input variables of the kernel function, c is an intercept parameter of the kernel function, and the value of c is 0.5, d is an exponential parameter, and the value of d is 5;

2b) using the kernel function selected in 2a) to solve the kernel matrix of the sample matrix Y:

C＝k(Y,Y)＝(<Y,Y>+c)^d

whereinn is the number of samples;

2c) performing eigenvalue decomposition on the kernel matrix C, i.e. C ═ V Σ V^TWherein Σ is a diagonal matrix, V is an eigenvector matrix corresponding to the eigenvalue, and T represents transposition;

2d) setting the contribution degree Z to be 0.9, calculating the contribution rate of each positive eigenvalue in the diagonal matrix sigma in the sum of all positive eigenvalues, sequentially accumulating the contribution rates of the positive eigenvalues according to the sequence of the eigenvalues in descending order, stopping when the accumulated sum is greater than the contribution rate Z, and recording the number of the accumulated eigenvalues as p;

2e) according to the first p mostLarge eigenvector of eigenvalue, and corresponding eigenvector matrix V_pThen the sparse dictionary A ═ V_p。

Step 3, generating an initial observation matrix phi₀。

And (3) randomly generating a Gaussian random matrix with the size of m × w as an initial observation matrix by setting the number of the measured values as m:φ_iis to initialize an observation matrix phi₀I-th row vector of (1, 2, …, m).

Step 4, combining a training sample matrix Y and an initial observation matrix phi according to the sparse dictionary A₀And obtaining a perception matrix G.

A sparse dictionary A, a training sample matrix Y and an initial observation matrix phi₀Multiplying the three to obtain a perception matrix: g ═ f (Φ)₀ ^T)^Tf (Y) A, where f represents the nonlinear mapping function and T represents the transposition.

And 5, training the perception matrix G to obtain a final observation matrix phi.

5a) Selecting the kernel function as a polynomial kernel function k (x, y) ═ f_k(<x,y>)＝(<x,y>+c)^dC is an intercept parameter of the kernel function, and the value of c is 0.5, d is an exponential parameter, and the value of d is 5;

5b) graham matrix G for computing perception matrix G^TG, solving the following optimization problem to obtain an observation matrix with minimum correlation with the sparse dictionary in the feature space:

wherein I_KIs an identity matrix with the size of K × K, and F represents the Frobenius norm of the matrix;

5c) will optimize the problemConversion to G^TG≈I_KLet the sensing matrix G be f (phi)₀ ^T)^Tf (Y) substitution of A into G^TG≈I_KAnd to the two sides, respectively, the sparse dictionary A is multiplied on the left, and the transposition A is multiplied on the right^TThe following relationship is obtained:

AA^Tf(Y)^Tf(Φ₀ ^T)f(Φ₀ ^T)^Tf(Y)AA^T≈AA^T；

5d) for the above AA^TPerforming singular value decomposition, i.e. AA^T＝VΛV^TWhereinis an orthogonal matrix obtained after singular value decomposition,is a diagonal matrix, n is the number of training samples;

5e) mixing AA^T＝VΛV^TSubstituting the formula in 5c) to obtain the following relation:

VΛV^Tf(Y)^Tf(Φ₀ ^T)f(Φ₀ ^T)^Tf(Y)VΛV^T≈VΛV^T

5f) according to the sample matrix Y and the initialized observation matrix phi₀Computing a transition kernel matrix K in combination with a kernel function_nm：

5g) Reducing the formula of 5e) to:

wherein the intermediate variable

5h) Writing the diagonal matrix Λ as diag (λ)₁,…,λ_n) Introduction of intermediate variables_mnWritten in the form of row vectors_mn＝[τ₁,…,τ_m]^TAnd combining the objective functionsIs converted intoWherein v is_i＝[λ₁τ_i,1,…,λ_nτ_i,n]；

5i) Will be provided withThe expansion is as follows:

meanwhile, an error matrix is defined as:

5j) according to a transition kernel matrix K_nmAnd an orthogonal matrix V initialized by the following equation_mn：

5k) Setting the maximum value of the cycle times t as m, wherein t starts to cycle from 1;

5l) in the t-th cycle, an error matrix E is calculated_tAnd decomposing the characteristic value to obtain the maximum characteristic value ξ_tAnd the corresponding characteristic vector u_t；

5m) update according to the following formulaEach row vector of (a):

5n) cycle m times, an intermediate variable with minimal correlation is obtained

5o) computing a transition kernel matrix

5p) solving the following formula to obtain an observation matrix phi subjected to coupling optimization:

wherein k (x, y) ═ f_k(<x,y>)。

And 6, acquiring a test image.

Taking the 10 th spectral band image of the three groups of hyperspectral images as a test image, and drawing the test image into column vectors which are respectively marked as e₁,e₂,e₃。

Step 7, according to the dictionary A obtained in the step 2 and the observation matrix phi obtained in the step 4, the three test images e in the step 5 are subjected to nuclear compression sensing₁,e₂,e₃And carrying out nonlinear compression imaging.

7a) According to the compression observation equation M-G beta form, calculating a measurement vector M and a perception matrix G:

7b) obtaining a sparse coefficient by adopting a least square algorithm according to the measured value vector M and the perception matrix G obtained by calculation

Step 8, reconstructing the original image according to the sparse coefficient β and the sparse dictionary A

The pre-image method is used for reconstructing an original image through the following formulai＝1,2,3：

Wherein u is_pThe p-th column of the unit orthogonal basis is represented, wherein p is 1,2, …, w is the number of pixel points of the hyperspectral image, and A β is [ c ]₁,c₂,…,c_j,…,c_n]^T，c_jJ element, f, representing A β_kFor the previously selected polynomial kernel function,is f_kIs the inverse function of (c).

The effects of the present invention can be further illustrated by the following experiments:

1) conditions of the experiment

Three groups of hyperspectral images used in the experiment are typical AVIRIS hyperspectral data: IndianPines, Moffet and WashtonDC. IndianPines data were obtained in 1992 from AVIRIS sensors imaging the northwest agricultural area of indiana, and Moffet images were obtained in 1992 from AVIRIS sensors imaging the Moffet area of california at 8 months; the wavelength range of the two groups of images is 0.4 um-2.5 um, 224 spectral bands are totally obtained, 200 spectral bands are obtained after all pixels are removed to be 0 and opaque wave bands are removed, and the spatial resolution is 20 m. The WashtonDC image is obtained by imaging a Washtinton DC Mall area by a HYDICE spectrometer, the wavelength range is 0.4-2.5 um, 210 spectral bands are totally obtained, 191 bands are selected after preprocessing, and the spatial resolution is 2.8 m. The size of IndianPines images is 145 × 145 × 200, the size of Moffet images is 145 × 145 × 200, and the size of WashtonDC images is 145 × 145 × 191.

Experiment simulation environment: the software MATLAB 2012R is used as a simulation tool, the CPU is AMD A8-5550M, the dominant frequency is 2.10GHz, the memory is 16G, and the operating system is Windows 7 flagship edition.

And randomly extracting the 10 th spectral band from each group of hyperspectral images as a test image, and taking images of other spectral bands of each group of hyperspectral images as training samples.

2) Emulated content

Simulation 1: under different sampling rates of 0.1-20%, the method of the invention and the existing method using the Gaussian random matrix as the observation matrix are respectively adopted to carry out nonlinear compressed sensing simulation experiments on the test images, and the experimental results are shown in Table 1.

TABLE 1 comparison of experiments for two methods at different sampling rates

As can be seen from Table 1, with the continuous increase of the sampling rate, the PSNR of the two methods continuously rises, which shows that the reconstruction effect is steadily increased, but the method of the present invention has the largest improvement amplitude. Under the same sampling rate, the PSNR of the method is highest, and the reconstruction effect is best.

Simulation 2: under the sampling rate of 10%, the method of the present invention and the existing method using the gaussian random matrix as the observation matrix are respectively adopted to perform the compressed sensing simulation experiment on the test image IndianPines, and the experimental results are shown in fig. 2 and fig. 3, wherein:

FIG. 2(a) is an original image of the 10 th spectral band of the test image IndianPines;

FIG. 2(b) is a reconstructed image using a conventional method using a Gaussian random matrix as an observation matrix, wherein the PSNR is 34.8322 dB;

FIG. 2(c) is a reconstructed image with a PSNR of 37.3694dB using the method of the present invention;

figure 3 is a plot of the mean square error for two methods at different sampling rates.

Simulation 3: under the sampling rate of 10%, the method of the present invention and the existing method using a gaussian random matrix as an observation matrix are respectively adopted to perform a compressed sensing simulation experiment on the test image Moffet, and the experimental results are shown in fig. 4 and 5, wherein:

FIG. 4(a) is an original image of a 10 th spectral band of a test image Moffet;

FIG. 4(b) is a reconstructed image using a conventional method using a Gaussian random matrix as an observation matrix, which has a PSNR of 47.2954 dB;

FIG. 4(c) is a reconstructed image with a PSNR of 49.9599dB using the method of the present invention;

figure 5 is a plot of the mean square error for two methods at different sampling rates.

And (4) simulation: under the sampling rate of 10%, the method of the present invention and the existing method using the gaussian random matrix as the observation matrix are respectively adopted to perform the compressed sensing simulation experiment on the test image WashtonDC, and the experimental results are shown in fig. 6 and 7, wherein:

FIG. 6(a) is an original image of the 10 th spectral band of the test image WashtonDC;

FIG. 6(b) is a reconstructed image using a conventional method using a Gaussian random matrix as an observation matrix, which has a PSNR of 44.1374 dB;

FIG. 6(c) is a reconstructed image with a PSNR of 46.4503dB using the method of the present invention;

figure 7 is a plot of the mean square error for two methods at different sampling rates.

As can be seen from the experimental results of fig. 2, 4 and 6, under the same sampling rate, compared with the existing method using a gaussian random matrix as an observation matrix, the PSNR of the method is higher, and the reconstruction effect is the best.

As can be seen from the graphs of fig. 3, 5 and 7, as the sampling rate increases, the mean square error of both methods shows a downward trend; but under the same sampling rate, the mean square error of the recovered image is smaller, which shows the superiority of the method.

Claims (3)

1. A spectral imaging method based on adaptive coupling observation and nonlinear compression learning comprises the following steps:

(1) selecting three groups n₁×n₂×n₃The hyperspectral image except the 10 th spectrum segment randomly selects n spectrum segments to construct a sample matrix: y ═ Y₁,y₂,…,y_j,…,y_n]Wherein n is₁×n₂Representing the size of the hyperspectral image, n₃Is the total number of spectral bands, y, of the hyperspectral image_jRepresents the column vector drawn from the jth spectral slice, j is 1,2, …, n, n isThe number of training samples;

(2) training a dictionary by utilizing a training sample matrix, and calculating a sparse dictionary of the training samples by adopting a Kernel Principal Component Analysis (KPCA) method, and recording the sparse dictionary asn is the number of samples, and K is the number of atoms of the obtained sparse dictionary;

(3) setting the number of measured values as m, and randomly generating a Gaussian random matrix as an initial observation matrixi＝1,2,…,m，φ_iIs to initialize an observation matrix phi₀A row vector of (a);

(4) combining a training sample matrix Y and an initial observation matrix phi according to a sparse dictionary A₀Obtaining a sensing matrix G ═ f (Φ)₀ ^T)^T(Y) A, training a perception matrix G to obtain a final observation matrix phi, wherein f is a nonlinear mapping function, and T represents transposition;

(5) drawing three groups of 10 th spectral band images into column vectors as test samples, and respectively marking as e₁,e₂,e₃；

(6) According to the observation matrix phi obtained in the step (4), three test samples e in the step (5) are subjected to₁,e₂,e₃Carrying out nonlinear compression imaging to obtain a measured value M-f_k(<Φ,e_i>) Wherein e is_iDenotes the ith group of test specimens, i from 1 to 3, f_kIs a selected kernel function;

(7) computing sparse coefficients using least squaresWhereinRepresents a pseudo-inverse;

(8) reconstructing the original image by a pre-image method according to the sparse coefficient β and the sparse dictionary AWhereinThe restored image for the ith group.

2. The method of claim 1, wherein the perception matrix G is trained in step (4) by:

2a) selecting the kernel function as a polynomial kernel function k (x, y) ═ f_k(<x,y>)＝(<x,y>+c)^dC is an intercept parameter of the kernel function, and the value of c is 0.5, d is an exponential parameter, and the value of d is 5;

2b) computing the gram matrix of the perception matrix G, i.e. G^TAnd G, solving the following optimization problem to obtain an observation matrix with minimum correlation with the sparse dictionary in the feature space:

$min\left|\right|G^{T}G-I_K|{|}_F^2,$

wherein I_KIs an identity matrix with the size of K × K, and F represents the Frobenius norm of the matrix;

2c) converting the above optimization problem to G^TG≈I_KLet the sensing matrix G be f (phi)₀ ^T)^Tf (Y) A in (G)^TG≈I_KAnd to the two sides, respectively, the sparse dictionary A is multiplied on the left, and the transposition A is multiplied on the right^TThe following relationship is obtained:

AA^Tf(Y)^Tf(Φ₀ ^T)f(Φ₀ ^T)^Tf(Y)AA^T≈AA^T

2d) for the above AA^TPerforming singular value decomposition, i.e. AA^T＝VΛV^TWhereinis an orthogonal matrix obtained after singular value decomposition,is a diagonal matrix and n is the number of training samples.

2e) Mixing AA^T＝VΛV^TSubstituting the relation in 3c) to obtain the following relation:

VΛV^Tf(Y)^Tf(Φ₀ ^T)f(Φ₀ ^T)^Tf(Y)VΛV^T≈VΛV^T

2f) according to the sample matrix Y and the initialized observation matrix phi₀Computing a transition kernel matrix K in combination with a kernel function_nm：

2g) Simplifying the relation of 3e) to:

$\begin{array}{c}{\mathrm{V\&Lambda;V}}^T{K}_{nm}{K}_{nm}^K{\mathrm{V\&Lambda;V}}^T\&ap;{\mathrm{V\&Lambda;V}}^T\\ {\mathrm{\&Lambda;\&Gamma;}}_{mn}^T{\mathrm{\&Gamma;}}_{mn}\mathrm{\&Lambda;}\&ap;\mathrm{\&Lambda;}\end{array},$

wherein,

2h) writing the diagonal matrix Λ as diag (λ)₁,…,λ_n) Will be_mnWritten in the form of row vectors_mn＝[τ₁,…,τ_m]^TAn objective functionConversion to:wherein v is_i＝[λ₁τ_i,1,…,λ_nτ_i,n]^T；

2i) Mixing the aboveThe expansion is as follows:meanwhile, an error matrix is defined as:

$E_t=\mathrm{\&Lambda;}-\underset{i=1,i\&NotEqual;t}{\overset{m}{\&Sigma;}}{v}_i{v}_i^T;$

2j) according to a transition kernel matrix K_nmAnd an orthogonal matrix V initialized by the following equation_mn：

${\mathrm{\&Gamma;}}_{mn}=K_{nm}^{T}V$

2k) Setting the maximum value of the cycle times t as m, wherein t starts to cycle from 1;

2l) in the t-th cycle, an error matrix E is calculated_tAnd decomposing the characteristic value to obtain the maximum characteristic value ξ_tAnd the corresponding characteristic vector u_t，

2m) update by the following equationEach row vector of (a):

$\sqrt{{\mathrm{\&xi;}}_t}{u}_t={\&lsqb;{\mathrm{\&lambda;}}_1{\mathrm{\&tau;}}_{t,1}^{\&prime;},...,{\mathrm{\&lambda;}}_n{\mathrm{\&tau;}}_{t,n}^{\&prime;}\&rsqb;}^T$

2n) are cycled m times to obtain a signal with minimal correlation

2o) computing a transition kernel matrix

2p) solving the following formula to obtain an observation matrix phi subjected to coupling optimization:

wherein k (x, y) ═ f_k(<x,y>)。

3. The method of claim 1, wherein said step (8) is performed by reconstructing the original image by a pre-image method according to the sparse coefficient β and the sparse dictionary AThe method is carried out by the following formula:

${\hat{e}}_i=\underset{p=1}{\overset{w}{\&Sigma;}}<e_i,u_p>u_p=\underset{p=1}{\overset{w}{\&Sigma;}}{f_k}^{-1}\left(\underset{j=1}{\overset{n}{\&Sigma;}}{c}_{j}k\right(y_j,u_p\left)\right)u_p$

wherein u is_pThe p-th column of the unit orthogonal base is expressed, p is 1,2, …, w is the number of the hyperspectral image pixel points, A β is [ c ]₁,c₂,…,c_j,…,c_n]^T，c_jJ element, f, representing A β_kFor a previously selected polynomial kernel, f_k ^-1Is f_kIs the inverse function of (c).