Directly Estimating Endmembers for Compressive Hyperspectral Images

The large volume of hyperspectral images (HSI) generated creates huge challenges for transmission and storage, making data compression more and more important. Compressive Sensing (CS) is an effective data compression technology that shows that when a signal is sparse in some basis, only a small number of measurements are needed for exact signal recovery. Distributed CS (DCS) takes advantage of both intra- and inter- signal correlations to reduce the number of measurements needed for multichannel-signal recovery. HSI can be observed by the DCS framework to reduce the volume of data significantly. The traditional method for estimating endmembers (spectral information) first recovers the images from the compressive HSI and then estimates endmembers via the recovered images. The recovery step takes considerable time and introduces errors into the estimation step. In this paper, we propose a novel method, by designing a type of coherent measurement matrix, to estimate endmembers directly from the compressively observed HSI data via convex geometry (CG) approaches without recovering the images. Numerical simulations show that the proposed method outperforms the traditional method with better estimation speed and better (or comparable) accuracy in both noisy and noiseless cases.


Introduction
Hyperspectral images (HSI) are collections of hundreds of images that have been acquired simultaneously in narrow and adjacent spectral bands, typically by airborne sensors [1,2]. Through the continued development of sensing technology, the spectral and spatial resolution for HSI has increased significantly. For example, the NASA Jet Propulsion Laboratory's Airborne Visible Infra-Red Imaging Spectrometer (AVIRSI) covers the wavelength region from 0.4~2.5 microns using 244 spectral channels at the nominal spectral resolution of 10 nm [3]; the spatial resolution of the hyperspectral imager in the Tiangong 1 aircraft is 5 m [4]. High spectral and spatial resolution results in HSI providing a wealth of information for accurate target detection and identification, leading to many applications including environmental monitoring, agriculture planning, and mineral exploration [5]. However, it also makes the volume of data very large, which introduces a significant challenge to data transmission, storage and analysis. Due to the extremely large volume of HSI data, compression technology has received considerable interest in recent years.
In conventional HSI sensing systems, the full data are acquired and are then compressed before transmission. This paradigm has several disadvantages: first, all the data should be stored; second, the computationally costly implementation of the compression is required to reside on board, housed within the sensing modality. Typically, the sensor platform is a severely resource-constrained environment such as a plane or satellite. As an alternative to the conventional sensing systems, compressive sensing (CS) [6,7] is an effective approach to acquire and compress the data in only one step. CS theory shows that only a small collection of a sparse or compressible signal contains enough information for stable signal recovery. Distributed CS (DCS) extends the single signal CS to multiple signals [8][9][10]. By exploiting both intra-and inter-signal correlation structures, DCS can reduce the number of measurements of each signal effectively, saving on the costs of data storage, communication and processing. DCS is very suitable for multi-channel applications, such as HSI.
Blind hyperspectral unmixing (HU) is one of the most prominent research topics in signal processing (SP) for hyperspectral remote sensing [11,12]. Blind HU aims to identify endmembers present in a captured scene, as well as their proportions [13]. There are many methods for blind HU such as pixel purity index (PPI) [14], N-FINDR [15], vertex component analysis (VCA) [16], SSCBSS [17], hypGMCA [18], and modified VCA (MVCA) [19], which are all based on the Nyquist sampling theorem. There are also some HU methods based on the CS theory, such as CSU [20] and the method proposed in [5], but they all assume that the endmembers are known as side information. Endmember estimation is a key step to identify the materials in HSI, and in many applications, the endmembers are unknown.
The traditional method for endmember estimation under the CS/DCS framework consists of 2 steps: (1) recovering the HSI data by CS/DCS methods and (2) estimating the endmembers from the recovered data by HU methods. The recovery step takes considerable time and also introduces errors into the estimation step, which will degrade the speed and accuracy of the endmember estimation.
In this paper, by designing a type of coherent measurement matrix, we propose a novel method that estimates the endmembers directly from the compressive HSI with convex geometry (CG) approaches, which outperforms the traditional method with better estimation speed and better (or comparable) accuracy.
The paper is structured as follows. The necessary theoretical background and notations are provided in Section 2. In Section 3, we describe the proposed method in detail. The performance of the proposed method is demonstrated in Section 4 in comparison to the traditional method. We conclude the paper in Section 5. Important acronyms used in this paper are listed in Table 1. HU refers to any process that separates the pixel spectral from a hyperspectral image into a collection of constituent spectral or spectral signatures, called endmembers and a set of fractional abundances, one set per pixel [12]. Mixing models can be characterized as either linear or nonlinear [12,13]. The linear mixing model (LMM) is a very representative model for HSI, and it is an acceptable approximation of the light scattering mechanisms in many real scenarios. In this paper, we only focus on the LMM.  Figure 1.
Owing to physical constraints, the proportions are non-negative and satisfy the full additivity constraint: Obviously, row of is [ ], and row of is [ ]. From the signal processing aspect, can be seen as the source matrix whose th column contains the proportion of source at each pixel; is the mixing matrix, and is the mixture matrix.

Convex Geometry Approaches for HU
There are several key approaches for HU, including convex geometry (CG) approaches, statistical approaches, sparse regression approaches and nonnegative matrix factorization [12,13]. CG approaches are very popular and effective in HU. A vast majority of HU developments, if not all, are directly or intuitively related to concepts introduced in CG studies [13]. In this paper, we only focus on the CG approaches.
We introduce some mathematical notations in convex analysis: spanned space, affine hull and convex hull.
Note that the vertices of conv{ 1 , ⋯ , } are , ⋯ , ; thus, in CG approaches, the inference of the endmember matrix is equivalent to identifying the vertices of the simplex.

Compressive Sensing
CS theory indicates that if a signal is sparse or compressive in some basis, it can be exactly recovered by a small number of measurements, much less than the number required by the Nyquist sampling theory.
be a sparse signal ( ≪ ). The sparse basis is ϵ × with a sparse coefficient vector ϵ . The signal can be denoted as with ‖ ‖ 0 = , where ‖ ‖ 0 denotes the number of non-zero entries in . is the × measurement matrix, where < . The observation vector consisits of linear projections of : where = is called the sensing matrix. The design of is critical for CS. A sufficient condition for stable signal recovery is that satisfies the RIP (Restricted Isometry Property) [21]. For each integer = 1, 2, ⋯, define the isometry constant δ of the matrix as the smallest number such that (1 − δ )‖ ‖ 2 2 ≤ ‖ ‖ 2 2 ≤ (1 + δ )‖ ‖ 2 2 holds for all -sparse vectors . We will loosely say that the matrix obeys the RIP of order if δ is not too close to one [7]. The correlation between and is shown in Figure 3.
The length of is much smaller than the length of , and thus, Equation (9) is underdetermined and the solution is ill-posed [23], for the matrix has more columns than rows. The most original approach for solving this problem is to find the sparsest vector , which seeks a solution to the 0 minimization problem The 0 minimization is NP-hard and computationally intractable [23]. Fortunately, it has been proven that the 1 minimization method can also exactly recover the signal under some conditions [7,21]. The 1 minimization is given by: The 1 minimization is also called the basis pursuit (BP), whose computational complexity is ( 2 3/2 ) [24]. The recovery speed is very slow, especially when is very large in the HSI application.

Distributed Compressive Sensing
DCS [8][9][10] is a combination of distributed source coding (DSC) and CS. In the DCS framework, multichannel sensors measure signals that are each individually sparse in some domain and also correlated from sensor to sensor. The DCS theory rests on a concept called the joint sparsity of a signal ensemble. There are three joint sparse models (JSM): JSM-1, JSM-2 and JSM-3. In JSM-1, all signals are sparse and have common sparse components, while each signal contains sparse innovation parts. In JSM-2, all signals share the same support set with different amplitudes. In JSM-3, all signals have non-sparse common parts and sparse innovations.
JSM-2 is the most concise model, as shown in Figure 4, and has been applied in compressive HSI [27]. A key prior that will be essential for compressive HSI is that each source image/mixture image is piecewise smooth in the spatial domain, implying a sparse representation in the DCT (discrete cosine transform) domain or the wavelet domain.
According to the assumption above, the image in each spectral band has a sparse representation, whose observations measured by the CS method can be written as From our assumption; each source image is also sparse in some domain; and thus; the sparse representation of each mixture (spectral image) has the same sparse location with different coefficients due to the different mixing parameters of each source.
The SOMP (Simultaneous Orthogonal Matching Pursuit) method is proposed in [28] to reconstruct all of the signals that fall into the JSM-2 simultaneously, and this algorithm outperforms the OMP algorithm when dealing with multiple signals [8] and in compressively sensed HSI reconstruction applications [27].
Hence, the compressive observation of HSI can be denoted as To summarize, the task of endmember estimation from compressive HSI can be described as follows: given the compressive observation matrix , as well as the measurement matrix and the sparse basis , estimate the endmember matrix .

Traditional Method for Endmember Estimation from Compressive HSI
The framework of the traditional method for solving the problem mentioned above is shown in Figure 5. It contains two steps: in the first step, it recovers the hyperspectral image ( = 1, 2, ⋯ , ) by solving the DCS problem described in Equation (12) and then estimates the endmember matrix (from the recovered mixtures) by solving the endmember estimation problem, as shown in Equation (3). In Figure 5, ̂ is the recovered value of , and ̂ is the estimated value of . The aim of endmember estimation is to estimate the endmember matrix , not the image matrix , which is only an intermediate representation used to calculate . However, the recovery of the images is a necessary step in the traditional method, as shown in Figure 5. It consumes a great deal of time and may also introduce errors to the estimation step.
In the next section, we will propose a new method that can estimate directly from the compressive HSI data without recovering the images, which leads to a better estimation speed.

Framework Description of the Proposed Method
As discussed in Section 2, if we can estimate endmembers directly from the compressive HSI data without recovering the images, omitting the recovery step will greatly reduce the complexity of computation; we will obtain a much better estimation speed and possibly also a better estimation accuracy.
The compressive observation of HSI is denoted as follows: where = can be regarded as the compressive measurement of the source . The value = = (∑ =1 ) is the compressive measurement of mixture , and it is also the mixture of all of the ( = 1,2, ⋯ , ). is the element of matrix at row column .
Equation (14) can be considered as LMM, as shown in Equation (3). Thus, we wish to estimate the endmember matrix directly via CG approaches such as PPI, N-FINDR, VCA, and MVCA.
Unfortunately, the properties of the measurement matrix make it impossible to directly use the CG approaches on Equation (14). Incoherence is a critical property indicating that the structures of the measurement matrix used in CS that, unlike the signals of interest, have a dense representation in the basis , and random matrices are largely incoherent with any fixed basis , making the sensing matrix hold RIP with overwhelming probability [7]. Hence, = is not sparse in basis , and [ ] (similar to the definition of [ ] in Section 2.1) cannot hold the non-negative and full additivity constraints, as shown in Equation (2), due to its dense and random character. We can see that: | [ ] ∈ }, = 1, 2, ⋯ , [ ] is the element in matrix at column row . From Equation (15), we can see that is not a convex hull, and [ ] is not a convex combination of endmember signatures. Actually, is the space spanned by { 1 , ⋯ , } (see Figure 2 for the relationship between the spanned space and the convex hull), and [ ] is a point in the space. In this condition, the endmember signatures { 1 , ⋯ , } are no longer vertices of a simplex, which means that we cannot use CG approaches to estimate the endmembers directly.
To the knowledge of the authors and the referenced materials, we have not found a measurement matrix that not only satisfies the incoherence (or RIP), but also makes = sparse in basis or [ ] hold the non-negative and full additivity constraints. It seems impossible to estimate endmembers from compressive HSI data directly. We propose a novel method, as shown in Figure 6. A coherent measurement matrix is used to compressively measure the HSI, and the endmembers can be directly estimated from the observations. First, we design a coherent measurement matrix that makes [ ] hold the non-negative and full additivity constraints.
The coherent matrix that does not hold the RIP cannot be used for exact and robust signal recovery, as mentioned above. In this paper, the aim is to estimate endmembers directly from the compressive HSI data, not to recover the HSI data. Therefore, we do not have to use an incoherent matrix (meaning, we do not have to use a random matrix). We give an example of such a coherent matrix. Let ϵ × be an identity matrix. We construct the measurement matrix using parts of . For example, we select one row in every ( = 1, 2, … ) rows from to compose ϵ R × (the compression ratio is , = round ( ), where round() is the operation of rounding towards the nearest integer). It is clear that [ ] holds the non-negative and full additivity constraints. With this kind of measurement matrix, it is possible to use CG approaches, such as PPI or VCA, to estimate the endmembers directly by solving the problem shown in Equation (14).
Other measurement matrices that make [ ] in = hold the non-negative and full additivity constraints can be used in this proposed method. This type of measurement matrix achieves undersampling of HSI. It will lose some proportion information. Therefore, the undersampled data cannot be used to recover the proportion information. If this information is required, one can recover it by the data captured by an incoherent measurement matrix, as shown in Figure 6.
In this paper, we use the VCA algorithm to estimate the endmembers directly in the proposed method. We assume the presence of pure pixels in the undersampled data as required by VCA, and we also assume the presence of all of the materials in the data. These assumptions are easy to realize for the increasing spatial resolution of hyperspectral sensors.
In the proposed method, we first use a coherent measurement matrix to compressively sense HSI and then use the VCA method to estimate the endmember directly from the HSI observations without recovering the images, which is a necessary step in the traditional method. As shown in the dashed parts in Figure 6, if the proportion information is required, one can use another incoherent measurement matrix to capture the global information of HSI, which can be used along with the estimated endmembers ̂ to recover the proportions [5,20].

Analysis of the Computational and Memory Complexity
Under the linear observation model, HSI data are in a subspace of dimension . Typically, is far less than . The dimensionality ranges from around 100 to 250, whereas ranges from about 3 to 20 [29]. In the VCA algorithm, the HSI data dimensionality is reduced by PCA or SVD. The computational complexity of the proposed method is ( 2 ) [16].
The traditional methods used in this paper are SOMP-VCA (SOMP for HSI data recovery and VCA for endmember estimation) and OMP-MVCA. The computational complexity of the SOMP is ( ), where is the sparsity of the signals. Typically, and are both larger than , so > 2 ≫ 2 . The computational complexity of VCA used in the traditional method is ( 2 ) (each spectral image length is , while the length of the compressive spectral image in the proposed method is ). The total computational complexity of the traditional method is ( + 2 ). In most cases, > 2 ≫ 2 , so the complexity can be simply denoted as ( ). The computational complexity of the SOMP-VCA is much larger than that of the proposed method. 2 > 2 , and thus, we can see that the computational complexity of the traditional method is larger than that of the proposed method, even if we use other CS recovery methods besides the SOMP algorithm. The computational complexity of OMP is ( ) when ≪ , otherwise, the computational complexity is the larger one of ( 3 ) and ( ). The computational complexity of MVCA is smaller than that of VCA. So the complexity of OMP-MVCA is dominated by OMP.
In the proposed method, the memory required is ( ) due to dimensionality reduction and data compression. In the SOMP-VCA method, the memory required by SOMP is ( ) and the memory required by VCA is ( ). ≫ > ; therefore, the memory required by the SOMP-VCA method can simply be ( ), which is much larger than that of the proposed method. Similar to the analysis of SOMP-VCA, the memory required by OMP-MVCA is simply ( ).

Evaluations of the Proposed Method
In this section, we first evaluate the practicability of the proposed method. Then, we compare the performance of the proposed method with the performance of the SOMP-VCA method.
The data used in this paper are (1) the semi-synthetic HSI of a rural suburb of Geneva and (2) real-work urban HSI, both of which are also used in [5] (This dataset is available at [30]. We acknowledge Mohammad Golbabaee, Simon Arberet and Vandergheynst for providing the dataset.). In the rural HSI, = 64 × 64, = 3, = 64, all the pixels are pure (which means that each pixel only contains one material), as shown in Figure 7. In the urban HSI, = 256 × 256, = 6, = 171 (the first 16 bands are used in this paper), and parts of the pixels are pure, as shown in Figure 8.
In this paper, the size of HSI data is reduced by SVD from × to × in VCA algorithm.
where θ is the angle between vector and ̂(the estimated value of ). Based on , we estimate the rms error distance: where rmsSAE measures the distance between and ̂ for = 1, 2, ⋯ , . (SAE stands for the Signature Angle Error).
We use ̃ to denote the coherent measurement matrix described in Section 3.1. In noisy cases, white Gaussian noise is added to the observation, i.e., ̃=̃+ . The SNR (signal-to-noise ratio) is from 20 to 40 dB with a step length of 10 dB. is the number of pixels of a hyperspectral image, and is the number of compressive measurements of an image. The compression ratio ( = round( )) is from 2 to 20 with a step length of 1 in the proposed method. We also consider the non-compression (Nyquist sampling) data, i.e., = 1. Each experiment is repeated 50 times to calculate the mean result. The platform used in this paper is with the following hardware and software: (1)   From Figures 9 and 10, we can see that as the compression ratio t grows, the rmsSAE does not change a lot, especially when the SNR is large, compared with the non-compression case t = 1 (this is due to the performance and the assumption of the VCA algorithm). However, the runtime decreases greatly as t increases compared with the non-compression case. We can see that the proposed method can estimate the endmember with comparable accuracy to the Nyquist-based method (t = 1) with much faster estimation speed as t increases, as long as the pure pixel and all material presence assumptions hold. We note that the value of t is user-defined as long as the presence of each material and pure pixel assumption holds.
To assess the performance of the proposed method for a large number of endmembers, we vary the number from = 10 to = 20. The endmember data are mineral signatures extracted from the U.S. Geological Survey (USGS) spectral library [16]. Each endmember signature consists of = 224 spectral bands. Each synthetic hyperspectral image has 4096 pixels.
From Figure 11 we can see that the rmsSAE values increase roughly with the increase of the number of endmembers under the same noise level. And the rmsSAE values also decrease roughly with the increase of SNR, as expected. Comparing Figure 11a,b, we can also see that the performance of the proposed is comparabe when = 10 and = 20, as long as the pure pixel and all material presence assumptions hold.

Comparison of the Proposed Method and the SOMP-VCA Method
Next, we compare the performance of the proposed method with the traditional SOMP-VCA and OMP-MVCA methods.
In SOMP-VCA and OMP-MVCA methods, the compression ratio ranges from 2 to 10 with a step length of 2. The incoherent measurement matrix is a random Gaussian matrix [21]. Each of the six signals (endmembers) in the urban HSI data has = 65,536 points. They require too much memory, as described in Section 3.2, when processed by a computer (the computer used in the paper with 4 GB RAM cannot afford enough memory to run the SOMP-VCA/OMP-MVCA algorithm when = 2 and 4). In the experiments below, we divide each signal into several segments, each of which consists of = 4906 points (The length of each signal in the rural HSI data is 4096). Therefore, is a × matrix, where = round( ). To recover the HSI data by the traditional method, the sparsity of each segment is required as a priori condition. We set = round( 100 ) in the experiments (Generally, it is very hard for people to get the exact value of in practical applications. This is also a disadvantage of traditional methods). In noisy cases, the SNR is from to 20 to 40 dB with a step length of 10 dB. Each experiment was repeated 50 times. From Tables 2 and 3, we can see that at the same compression ratio, the rmsSAEs of the proposed method are much smaller than the rmsSAEs of both SOMP-VCA and OMP-MVCA methods when it is used for both the synthetic and real data. The reason is that both SOMP-VCA and OMP-MVCA methods first recover the images and then estimate endmembers from the recovered images, and the recovery step will introduce error to the estimation step. The proposed method directly estimates the endmembers by the VCA method without the recovery step. We can see that in the noiseless case, the proposed method can estimate the endmembers exactly, as some values of rmsSAE in Tables 2 and 3 are 0.0000.   Tables 4 and 5. We can see that under the same condition, the standard deviation values of the proposed method are also smaller than the values of SOMP-VCA or OMP-MVCA methods. The results indicate that the convergence of the proposed method is the best among the three methods. From Tables 6 and 7, we can see that the time consumed by the proposed method is much smaller than the time consumed by SOMP-VCA or OMP-MVCA methods. As analyzed in Section 3, the proposed method estimates the endmembers in one step, while both SOMP-VCA and OMP-MVCA use two steps. Therefore, from Tables 2 to 7, the estimation accuracy and speed of the proposed method are both better than those of SOMP-VCA and OMP-MVCA methods. Different CS recovery algorithms affect the performance and the runtime for estimating the endmembers. To eliminate the influence of the choice of the CS method, we suppose that the sparsity of each signal is known and that the CS method can recover the HSI data accurately, and the recovered data are as accurate as the Nyquist-based data at the same noise level. From Figures 9 and 10, we can see that the performance of the proposed method ( > 1) is comparable to the performance of the method with Nyquist-based data ( = 1). Therefore, the performance of the proposed method is comparable to or better than the performance of the traditional CS base method. The same is true of the runtime time of the proposed method: it is less than that of the traditional CS base method with other CS recovery methods, as analyzed in Section 3.2.

Conclusions
In this paper, we proposed a new method to directly estimate the endmembers from the compressive observations of the HSI data, while traditional methods first have to recover the HSI data from the compressive observations and then estimate the endmembers. Simulation results demonstrated that the proposed method outperforms the traditional method with better estimation speed and better (or comparable) accuracy.