Abstract
Each pixel in the hyperspectral unmixing process is modeled as a linear combination of endmembers, which can be expressed in the form of linear combinations of a number of pure spectral signatures that are known in advance. However, the limitation of Gaussian random variables on its computational complexity or sparsity affects the efficiency and accuracy. This paper proposes a novel approach for the optimization of measurement matrix in compressive sensing (CS) theory for hyperspectral unmixing. Firstly, a new Toeplitz-structured chaotic measurement matrix (TSCMM) is formed by pseudo-random chaotic elements, which can be implemented by a simple hardware; secondly, rank revealing QR factorization with eigenvalue decomposition is presented to speed up the measurement time; finally, orthogonal gradient descent method for measurement matrix optimization is used to achieve optimal incoherence. Experimental results demonstrate that the proposed approach can lead to better CS reconstruction performance with low extra computational cost in hyperspectral unmixing.
1. Introduction
Compressive sensing (CS) theory [1, 2] is a new developed theoretical framework on signal sampling and data compression, which indicates that if a signal is sparse or compressible in a certain transform domain, the transformed higher-dimensional signal can be projected onto a lower dimensional space by a measurement matrix. It leads to nonadaptive measurement encoding on the signal at a rate far below the Nyquist sampling rate, converting from sampling the signal itself to sampling the information contained in the signal. Therefore it has recently gained more and more attention in various areas of applied mathematics, computer science, and electrical engineering.
Design of measurement matrix is a research hotspot in CS, and measurement matrix optimization has become an inevitable trend to construct a new measurement matrix system. In recent years, scholars have yielded many optimization methods [3–22] to design measurement matrix to reduce the minimum coherence of Gram matrix. These are typically fallen into three categories: iterative thresholding method [3–14], gradient iteration process [15–19], and Tensor product [20–22].
Zhang et al. [6] proposed that the Kronecker product measurement matrix based on orthogonal basis can maintain nonlinear correlation between columns’ vector from high-dimensional data. But Kronecker product [5, 6] leads to low sampling efficiency and poor computational complexity, which limit the deep study for measurement matrix.
By using iterative thresholding method, Elad [7] iteratively reduces the average mutual coherence using a shrinkage operation followed by singular value decomposition (SVD) step and shrinks the elements of Gram matrix. Lustig et al. [8] defined an incoherence criterion and proposed a Monte Carlo scheme for random incoherent sampling. Abolghasemi et al. [9] attempted a kind of nonuniform sampling by segmenting the input signal and taking samples with different rates from each segment. Duarte-Carvajalino and Sapiro [10] take advantage of an eigenvalue decomposition process followed by a KSVD-based algorithm to optimize measurement matrix and learn dictionary basis, respectively. Wang et al. [11] propose to generate colored random projections using an adaptive scheme. Next they [12] propose a variable density sampling strategy by exploiting the prior information about the statistical distributions of natural images in the wavelet domain. Although the algorithm is simple and easy to understand, the research [13, 14] found that iterative thresholding method has slow convergence speed and is easy to fall into local minimum problems in practical application. Meanwhile it damages Restricted Isometric Property (RIP) and eventually may cause the collapse of the BP algorithm.
By using gradient iteration process, Xu et al.’s algorithm [15] first shrinks and updates elements in Gram matrix with Equiangular Tight Frame (ETF). Li et al. [18] constructs the dimensional orthogonal matrix in SVD; Abolghasemi et al.’s algorithm [16] lies in the innovation of the gradient iteration process to obtain measurement matrix. Zhang et al.’s algorithm [17] adopts the spherical search steepest descent method. Tian et al.’s algorithm [19] shrinks Gram matrix with the orthogonal gradient factor matrix to reduce the maximum and average mutual coherence of measurement matrix. The appearance of different pursuit rates brings about the saddle-point steady-state solution, which only guarantees a local minimum solution.
Unfortunately, most algorithms neglect intrasensor correlations between the samples of high-dimensional data like 3D color image, video, and hyperspectral image; it may affect the multispectral features excessively and destroy the original structure of high-dimensional data and it ultimately affects the precision of hyperspectral unmixing.
Toeplitz-structured chaotic measurement matrix is one of deterministic measurement matrices in CS, which requires only independent variables and operations. But it still has three attractive weaknesses: optimal incoherence being unachieved, reconstruction precision being insufficient, and measurement time being unbearable.
To eliminate the weaknesses, we are inspired to work on optimizing the Toeplitz-structured chaotic measurement matrix to obtain better results from high-dimensional signals. Meanwhile, the possibility of fusing these two attractive optimized methods is certain: rank revealing QR factorization with eigenvalue decomposition from Xu et al.’s algorithm [15] and orthogonal gradient descent approach from Tian et al.’s algorithm [19] to obtain a new optimized method to overcome the computational complexity. This is my intuitive idea of this paper.
The key contribution of this work can be elaborated as follows:(1)In the domain of designing the measurement matrix, it is crucial to achieve high quality with implementing effective hardware. In this work, some pseudo-random chaotic elements can approximate the random structure component, which satisfies RIP property with overwhelming probability. So these pseudo-random chaotic elements are applied to Toeplitz-structured chaotic measurement matrix, which forms a new measurement matrix (TSCMM). Compared with the others, we attempt to prove that it satisfies RIP conditions.(2)According to the properties of the separated contents, Gram matrix is improved by a rank revealing QR factorization with eigenvalue decomposition to speed up convergence rate. The results of experiments show that the proposed methods can greatly reduce computation complexity building on the convergence and robustness.(3)From a practical point of view, high computational complexity imposes restrictions on achieving optimal incoherence; orthogonal gradient descent approach is proposed to acquire measurement matrix optimization. The improved scheme can effectively reduce the reconstruction error and acquire satisfied image quality compared to other conventional methods.
The rest of this paper is organized as follows. The main part of this paper starts with review of some measuring coherence criteria and develops to RIP in Section 2. Based on the previous analysis, the proposed approach for measurement matrix optimization in Section 3 is developed. Particularly motivated by TSCMM, QR factorization with eigenvalue decomposition and improved gradient descent approach is then suggested in the following three subsections, starting with Toeplitz-structured chaotic measurement matrix in Section 3.1, followed by QR factorization algorithm in Section 3.2 and an orthogonal gradient descent approach for measurement matrix optimization in Section 3.3. Experimental results in Section 4 and conclusions in Section 5 are presented.
2. Problem Formulation and Analysis
From the viewpoint of mathematics, CS sample procession is approximated by recovering from far incomplete measurements: where in some basis is sparse or compressible representation, while is so-called CS measurement matrix. Because of , the signal is measured through the projection by the measurement matrix which leads to being highly underdetermined.
If the D-dimensional sample signal and independent measurements result are unknown, the Kronecker product measurement matrix [5] can be expressed as . When each sensor obtaining its independent measurements is the same measurement matrix , the joint measurement matrix can be expressed as , where denotes the identity matrix, as shown in Figure 1.
The constants for the matrix are intrinsically tied to the singular values of all column submatrices of a certain size. If are matrices with restricted isometry constants (RIP) , the structure of Kronecker product matrices yields simple bounds for their RIP that can be expressed as
Considering the D-dimensional Kronecker sparsifying basis and a global measurement basis or frames obtained through a Kronecker product of individual measurement bases, the definition of mutual coherence is presented as
High-dimensional Kronecker compressive sensing (HKCS) [6] proposed the optimal synthetic sensing matrix by taking Kronecker products of individual optimal sensing matrix in each dimension. The optimal sensing matrix that minimizes the mutual coherence of the projection matrix can be expressed as
With the same sampling rate, matrices of HKCS have relatively smaller mutual coherence. It can be written as
It also indicates that the optimization process is dividable, which preserves the block feature of Kronecker product matrix and enables fast low-scale matrix computation. The overall video acquisition is decomposed as shown in Figure 2.
The high-dimensional Kronecker products measurement matrix is our optimization goal, as shown in Figure 3.
If is defined as Gram matrix and minimum square error cost function is defined as , the optimization problem can be written as
The goal of eliminating the correlation is to minimize the difference between Gram matrix and identity matrix in the form of Frobenius norm. Considering Kronecker product properties,
If the value can be replaced with its corresponding eigenvalue decomposition and , (7) can then become as follows:
Because is real diagonal matrices, . If , then .
Here, (8) can be reduced to
Supposing that is the elements in , gradient decrease iteration method is used to minimize mean square error (MSE). , where is step size and . is gradient value of ; then . According to Gerschgorin theorem, the column coherence of can be deduced as follows:
If , the infimum of the column coherence is called Welch bound:
Equiangular Tight Frame (ETF) is derived from (11) if the constraints are equality.
3. The Proposed Approach
Based on previous conclusions, the proposed algorithm aims to optimize the Toeplitz-structured chaotic measurement matrix to obtain better results from hyperspectral unmixing. Therefore, the research content in this paper mainly consists of three parts: designed TSCMM, optimized Gram matrix, and orthogonal gradient descent approach, as shown in Figure 4.
The study is given by taking the following methods: firstly, to obtain easy hardware implemented, pseudo-random chaotic elements are used to form a new Toeplitz-structured chaotic measurement matrix (TSCMM), as discussed in Section 3.1; to overcome unbearable Cost Time, Gram matrix is improved by a rank revealing QR factorization with eigenvalue decomposition as discussed in Section 3.2; to achieve optimal incoherence, orthogonal gradient descent method for measurement matrix optimization is presented in Section 3.3. Finally, the improved scheme is presented through explicit analysis and discussion.
3.1. Xu-TSCMM
Recently, [23] is written by myself completely and probed into its initial theory and research. Discrete chaotic system function is proposed to generate a series of pseudo-random numbers. Based on those elements, Toeplitz-structured chaotic measurement matrix (TSCMM) is produced to guarantee the incoherence criterion. To reduce the building time of TSCMM, Circulant/block-diagonal splitting structure is attached on TSCMM. Although about one-third of matrix values are eliminated, the measurement matrix is proved to satisfy Johnson-Lindenstrauss (J-L) lemma and achieves the goal of satisfying RIP. Logistic map [24], as the simplest dynamic systems evidencing chaotic behavior, is described as follows:where is the discrete state. While parameter is 4, the sequence satisfies beta distribution with and , and the next probability density function has been used for simulation.
Set as the output sequence generated by (12) with initial condition , and let denote the regularization of as the following form: .
Approximately, can be considered as random variable, and it satisfies the following distribution: . Then, by selecting different initial conditions , one can obtain vectors with dimension , which enables us to construct the following matrix scaled by:
Here, (13) is called the beta-like matrix.
According to (13), set one initial condition , and generate a sequence in the chaotic system. Then, Toeplitz-structured matrix is constructed in the following form:
Here, is for normalization and is called Toeplitz-structured chaotic measurement matrix (TSCMM), which meets J-L theorem.
3.2. Duarte-ETF Method
Minimum coherence property of ETF has been the main target to find a feasible solution. It is impossible to solve the problem exactly because of the complexity, while the structure of Gram matrix has changed so much that the selection of new units in the following step is very difficult.
To minimize (14), an optimization approach is adopted to reduce the maximum and average mutual coherence of measurement matrix. It shrinks Gram matrix based on ETF theory. The method can minimize the global mutual coherent coefficient of TSCMM by adjusting the eigenvalues above zero to the average value of the sum of these eigenvalues without changing the sum. After performing alternating minimization, the optimized measurement matrix can be constructed from the output Gram matrix with a rank revealing QR factorization with eigenvalue decomposition.
Theorem 1. Given a measurement matrix and representing matrix , there exists a matrix and Gram matrix , where is the column normalization from . If the real positive definite matrix has , the following equality and , where is the column of which is obtained.
Based on Theorem 1, to minimize the largest absolute values of the off-diagonals in the corresponding Gram matrix, we can determine the eigenvalues of Gram matrix by solving the following optimization problems:where is the element of Gram matrix. The optimization problem (15) is to minimize the square sum of the element of Gram matrix if the sum of the characteristic value remains constant. Though the eigenvalue decomposition of Gram matrix has , the value has been gradually approaching to minimize the square sum of . Because real symmetric matrix eigenvalue decomposition is orthogonal, the square sum of nondiagonal elements from Gram matrix gradually decreases and further it attains the effect as follows:
Finally, is obtained after several iterations, which is also the optimal Gram matrix .
3.3. The Orthogonal Gradient Descent Approach
According to the definition of Gram matrix, Gram matrix is the product of measurement matrix and sparse matrix. Therefore, the optimal measurement matrix can be directly derived from the optimal Gram matrix. However, as for overcomplete sparse representation based on redundant dictionary, it is very hard to design an effective algorithm to construct measurement matrix. In order to solve this problem, orthogonal gradient descent method is employed to get the optimal measurement matrix .
If the optimal Gram matrix is obtained, the optimal measurement matrix is as follows:
And then the complex problem from optimal Gram matrix is transformed into simple minimum of . By determining the derivative of , a skew-symmetric matrix with the measurement matrix and the gradient matrix is used to obtain revision factor:
Next, if is identity matrix, the orthogonal matrix can be expressed through the Cayley transform to ensure the positive definiteness of the revision factor:
The orthogonal gradient factor matrix can be obtained to update the gradient direction:
Combining (19) and (20), (21) can be rewritten as follows:
Finally, update measurement matrix with the orthogonal gradient factor matrix :
Because is updating ratio, gradually converges at the minimum value. Then the optimal measurement matrix is the goal of our pursuit. According to the linearity property of Toeplitz measurement matrix [25], satisfies J-L property with overwhelming probability. From Theorem 1, it has been proven that J-L condition can replace RIP condition. So also satisfies RIP property with overwhelming probability.
The flow diagram of the proposed method will be given as shown in Algorithm 1.
| ||||||||||||||||||||||||||||||||||
4. Experiments and Result Analysis
To illustrate the effectiveness of the proposed approach, the most widely used hyperspectral images in unmixing, such as Cuprite, Urban, and Jasper Ridge, were selected in the spectral range from 380 nm to 2500 nm; each channel band width is up to 9.46 nm. All high-dimensional data is provided by the standard hyperspectral library of 224 bands which comes from [26]. To reduce complexity, there are only 128 × 128 pixel blocks of original image which starts from the th pixel. Further, only 8 channels (from 11 to 81 bands every 10 bands) were remained due to dense water vapor and atmospheric effects.
In the course of the experiment, the signal sparsity method is Fourier basis and reconstruction algorithm is Stagewise Orthogonal Matching Pursuit (StOMP) [27]. Various kinds of measurement matrix (Circulant, Toeplitz, Toeplitz-structured chaotic measurement matrix (TSCMM) [23], Elad-Optimization Method (EOM) [7], and Duarte-Carvajalino and Sapiro’s Method (Duarte-ETF) [10]) are employed to illustrate the effectiveness of proposed approach for hyperspectral unmixing.
To demonstrate the efficiency of these methods, traditional evaluation methods can generally be divided into two categories: subjective assessment and objective evaluation. Mean Squared Error (MSE) and Peak-Signal-to-Noise Ratio (PSNR) as one of the most important indices from objective evaluation determine the quality of recovery image, while Cost Time (CT) verifies the efficiency of the proposed approach. They all testify experiment results of recovery signals, built on laptop with AthlonTM Processor 1.60 G HZ, 1 GB RAM, Matlab 7.0, and Windows XP operation platform.
4.1. Cuprite
To illustrate the use of the hyperspectral analysis process, a sample scene covers the Cuprite mining district in western Nevada, USA, from NASA’s Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) is provided. The data provided here is one of the most widely used hyperspectral images in unmixing study. There are 210 wavelengths ranging from 400 nm to 2500 nm, resulting in a spectral resolution of 10 nm.
In Figure 5, the first image (Top) was taken in blue light, the second image (middle) was taken in red light, and the third image (bottom) was taken in near infrared light centered at a wavelength of 750 nanometers.
Figure 5 presents the subjective evaluation by Circulant, Toeplitz, TSCMM, EOM, Duarte-ETF, and proposed method. Compared with other results, the performance from Figure 5(b) is the worst. The reason is that the elements from Circulant measurement matrix follow periodic repetition permutation, which does not satisfy RIP property with overwhelming probability. Figure 5(c) clearly demonstrates that Toeplitz measurement matrix can avoid the problem well. Because of the property of pseudo-random of chaotic sequence, the performance from TSCMM has further improved, as shown in Figure 5(d). The result from Figures 5(e)–5(g) shows that there is a significant impact on different bands using different optimization methods. Near infrared image is less affected by dust and gas. Visible blue channel has strong capability to penetrate water and visible red channel can reflect the health status of plants. Therefore, the sorted off-diagonal entries of the measurement matrix from EOM are likely more sparse and diagonal entries are more concentrated. The results are clearly shown in Figure 5(g) that the proposed approach had a significant performance compared to any others and closely resembles original image. Furthermore, the objective evaluations, which include Mean Squared Error (MSE), Peak-Signal-to-Noise Ratio (PSNR), and Cost Time (CT), can avoid artificial error and draw compelling conclusion. The results can be clearly seen from different methods on recovery Cuprite, as shown in Table 1.
Figure 5 and Table 1 report the recovery quality of the proposed method on recovery Cuprite. The following observations are summarized: of all evaluating indicators considered here, traditional Circulant had the worst performance in both subjective and objective evaluations; since the introduction of Toeplitz, the performance gets major improvement on image quality while improved Toeplitz-structured matrix method (TSCMM) is slightly better than classical Toeplitz matrix method; EOM has significant performance in image quality; however, the optimization process is usually an iterative process, which is also a very complicated and time-consuming process; Duarte-ETF has better contrast and lower computational complexity; the proposed method takes advantage of improved Toeplitz-structured matrix to speed up the convergence speed and improve traditional optimization method to better recovery high-dimensional image. Experimental results show that the proposed method has a better overall performance.
4.2. Urban
University of California Santa Barbara (UCSB) built an urban spectral library for the Goleta/Santa Barbara area. The hyperspectral data of Urban were acquired between late May and early June, 2001, using an ASD full range instrument on loan from the Jet Propulsion Laboratory. These spectra of Urban are characterized by 499 roofs, 179 roads, 66 sidewalks, 56 parking lots, 40 road paints, 37 types of vegetation, 47 types of nonphotosynthetic vegetation, 88 bare soil and beach spectra, 27 acquired from tennis courts, and 50 more from miscellaneous surfaces.
Experiments on the hyperspectral data of Urban demonstrate that the proposed scheme substantially improves the reconstruction accuracy. Clearly it comes to the same conclusion from Figure 6 that, compared with the previous five methods, the effect of proposed approach is obviously superior to any other methods and is the most similar to original image.
To evaluate and compare the proposed method, the following performance indices, such as Average Gradient (AG), Edge-Intensity (EI), Figure Definition (FD), Gray Mean (GM), Standard Deviation (SD), Space Frequency (SF), Variance (VAR), and Structural Similarity (SSIM), were used. It shows that the objective evaluation indices enhance the experiment rigor and convincing. The results are shown in Table 2.
The results of Table 2 show that the method has a higher performance than traditional Toeplitz or Circulant matrix method. Although improved Toeplitz-structured matrix method (TSCMM) is slightly better than classical Toeplitz, both classical optimization measurement matrix method (EOM) and proposed method have significant performance in image quality. Furthermore, the proposed method takes advantage of improved Toeplitz-structured matrix to speed up the convergence speed and improve traditional optimization method to recover better high-dimensional image. Experimental results show that the proposed method has a better overall performance.
4.3. Jasper Ridge
The hyperspectral image of Jasper Ridge was obtained on June 2, September 4, and October 6, 1992, which was calibrated to surface reflectance. The image was the most popular source to analyze with spectral mixture analysis using library endmembers representing green foliage, nonphotosynthetic vegetation, and soils characteristic of the site. Field-based vegetation was obtained from US Geological Service.
From Figure 7, it is obvious that the worst effects from traditional Circulant have reached being almost intolerable. While the results from TSCMM and Toeplitz are almost similar, the former has only slight improvement compared to the latter. On the other hand, the performance of the proposed approach was significantly improved compared to that of EOM or Duarte-ETF. The study concluded that the proposed approach had a significant performance compared to that of others. Furthermore, the objective evaluations are shown in Table 3. The results can be clearly seen from different methods on recovery Jasper Ridge.
The results of Table 3 show that the proposed method has a higher performance than traditional Toeplitz or Circulant matrix method. TSCMM takes advantage of improved Toeplitz-structured matrix to speed up the convergence speed and improve traditional Toeplitz or Circulant matrix method to recover better high-dimensional data. Although EOM has the lower column coherence and faster convergence, it weakens RIP condition and causes recovery performance degradation. While the absolute values from Duarte-ETF concentrate around mutual coherence, this can make the equivalent dictionary as close as possible to an ETF. But this algorithm has high computational complexity. Furthermore, experimental results from the proposed method show that the proposed method has a better overall performance.
4.4. Hyperspectral Unmixing
The fourth experiment performs an experimental evaluation of the accuracy of the standard hyperspectral unmixing districts known as Cuprite [28]. To reduce complexity, there are only 188 channels (3–103, 114–147, and 168–220 bands) that were remained due to dense water vapor and atmospheric effects. In the course of the experiment, the signal sparsity method is Fourier basis and reconstruction algorithm is StOMP.
The results are shown in Figure 8 for hyperspectral unmixing and closely resemble those obtained from hyperspectral data. Figure 8 compares the unmixing performance of the proposed method with different endmembers and the abundance from different endmembers is totally dissimilar. These features from the proposed method ensure clear and accurate the abundances.
(a) Alunite
(b) Andradite
(c) Pyrope
(d) Nontronite
(e) Dumortierite
(f) Kaolinite
(g) Chalcedony
(h) Kaolinite
(i) Buddingtonite
By accessing information from USGS Digital Spectral Library [29], the unmixing performance has been almost correct. Furthermore, different endmembers from USGS_1995_Library [30] are used to verify that the prediction model of hyperspectral unmixing scheme is accurate in Figure 8.
Apparently, the proposed unmixing results (blue thin thread) have a strong correlation with these endmembers (red thin thread) from USGS_1995_Library. From Figure 9(a), a steep sloping line from Alunite suggests that the unmixing endmember (blue thin thread) has remarkable similarity. This same conclusion has been made in studies from Figures 9(d)–9(h). On the other hand, the unsatisfied results from Figures 9(b) and 9(c) have been caused by smooth curve. The proposed method has good accuracy and is robust to traditional filtering, compression, cutting, and noise attack.
(a) Alunite
(b) Andradite
(c) Pyrope
(d) Nontronite
(e) Dumortierite
(f) Kaolinite
(g) Chalcedony
(h) Kaolinite
(i) Buddingtonite
5. Conclusions
In this paper, to overcome the limitation of Toeplitz-structured chaotic measurement matrix, an improved measurement matrix has been carried out in the hyperspectral unmixing process to achieve multiple endmembers of hyperspectral image. And in theory, it proves that this matrix has retained the RIP property with overwhelming probability. Experimental results demonstrate that the proposed method to design of measurement matrix leads to better CS reconstruction performance with low extra computational cost. Compared with some traditional measurement matrix, an improved method has highest technical feasibility, lowest computational complexity, and least computation time consumption in the same recovery quality. The proposed method can take the special advantage in hyperspectral unmixing process and explore the practical satellite system to remote sensing.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research is supported by Chongqing Engineering Laboratory for Detection, Control and Integrated System. The project is also funded by Key Technology Research and Industrialization of Fire Monitoring and Early Warning Sensor Network for High Voltage Transmission Line (KJZH17124). This research is funded by Chongqing Education Commission Foundation (KJ1400612). This project is also granted financial support from a Cooperative Project of Chongqing Technology and Business University (990516001).