Spatial-Spectral Joint Compressed Sensing for Hyperspectral Images

Compressed sensing is one of the key technologies to reduce the volume of hyperspectral image for real-time storage and transmission. Reconstruction based on spectral unmixing show tremendous potential in hyperspectral compressed sensing compared with other conventional algorithms that directly reconstruct images. In this article, a joint spatial-spectral joint compressed sensing scheme is proposed. In this scheme, compressed hyperspectral data are collected by spatial-spectral hybrid compressed sampling. As for the reconstruction, an objective function is developed by introducing the fidelity constraints of spatial and spectral measurements and the row sparsity constraint of abundance that guarantee the precise reconstruction of hyperspectral images and obtain endmembers and abundances as by-products accurately. An augmented Lagrangian type algorithm is meticulously elaborated to solve the above optimization problem. Extensive experimental results on several real hyperspectral datasets indicate that the proposed approach can achieve a reconstruction accuracy higher than those of other state-of-the-art methods. The efficiency and feasibility of the proposed scheme give it great potential in hyperspectral compressed sensing.


I. INTRODUCTION
Hyperspectral remote sensing has been extensively applied in fields such as the military, agriculture, and exploration fields, because of its ability to provide rich, detailed information on ground objects [1], [2]. However, with the continuous increase in resolution, massive volumes of hyperspectral images (HSIs) put very high stresses on airborne/satelliteborne imaging systems. Compressed sensing or compressive sampling (CS) is an effective technique for addressing this issue.
CS theory [3], [4] integrates signal acquisition and sensing, and signal is compressed as they are acquired. According to CS theory, the sampling rates (srs) for sparse or compressible signals can be far lower than Nyquist srs. Sparse representation have been widely used in image and video processing The associate editor coordinating the review of this manuscript and approving it for publication was Hengyong Yu .
[5]- [7]. The application of CS in the hyperspectral imaging field has also received wide attention. A series of hyperspectral compressive sensing (HCS) schemes have been proposed [8]- [20]. The special three-dimensional (3D) data structure allows diverse CS methods for HSIs. Note that HSIs can be either spatially sampled, similar to ordinary images, or spectrally sampled [9], [20]. However, spatial and spectral joint is a common method for HSIs processing [21]. Therefore, spatial-spectral hybrid CS method [10], [11], [14], [17], [18] is also proposed for hyperspectral compressive sampling. In this article, the spatial-spectral hybrid CS method proposed in our previous work [14] is used in the spatial-spectral joint CS for hyperspectral data acquisition.
Reconstruction algorithms have always been an important topic of research in the CS field, with HSIs being no exception. Conventional reconstruction algorithms generally directly reconstruct images based on prior information contained in the original data, such as the sparsity and total VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ variation (TV). For HSIs, this one-or two-dimensional prior information is generally extended to 3D information using algorithms, such as 3D sparsity, TV, and low-rank in 3DCS [10]. The structural sparsity of HSIs is also used as prior information for reconstruction [15], [22]. This algorithms are characterized by direct reconstruction of original data and can yield relatively good reconstructed data. However, due to the massive volumes of HSIs, the computational complexity of this algorithm is relatively high.
Since the special 3D structure of HSIs, tensor-based methods have been developed for HCS reconstruction in recent years [23]- [28]. In [23], the tensor Tucker decomposition was employed to describe the global spatial-spectral correlation among all HSI bands, and incorporated into weighted 3D TV regularization to enhance the performance of reconstruction. Xue et al. [25] proposed a nonlocal tensor sparse and low-rank regularization (NTSRLR) approach, which can encode essential structured sparsity of an HSI and explore its advantages for HCS reconstruction task. In order to exploit the spectral consistency and nonlocal self-similarity in the spatial dimension, the framework combining nonlocal lowrank regularization and hyper-Laplacian regularization was proposed in [24]. Chen et al. [28] assumed that the HSIs lies in a low-dimensional subspace instead of designing the regularization of the low-rank approximation, and proposed a subspace-based nonlocal tensor ring decomposition method (SNLTR) for HCS reconstruction. Although these tensorbased works can obtain satisfactory reconstruction results, the terrible computational complexity hinders the universality of the algorithm.
Compressive-projection principal component analysis (CPPCA) [9] and the relevant improved algorithms [29]- [32] reconstruct original data by estimating eigenvectors using principal component analysis (PCA) instead of directly recovering the original data. CPPCA exhibits relatively high reconstruction performance while effectively reducing the computational load. However, there is a critical defect in CPPCA. This reconstruction algorithm fails at low srs, and its reconstruction performance decreases rapidly when sr is lower than 0.1.
As a simple, effective prior assumption, linear mixed models (LMMs) are extensively applied in spectral unmixing. Under an LMM assumption, the spectrum of any pixel can be considered a mixture of the characteristic signals of the endmembers and the corresponding abundances. Thus, in recent years, researchers have used LMMs to reconstruct and decompose HSIs into endmember and abundance estimates and thereby have indirectly reconstructed original HSIs [12], [14], [16]- [18], [33]- [35]. The decrease in the volume of objective data for reconstruction not only significantly increases the reconstruction speed but also considerably improves the reconstruction performance. Wang et al. [14] proposed a spatial-spectral hybrid CS method. In the reconstruction process, this method extracts endmembers by the spatially compressed data and subsequently estimates abundances based on spectral compressed measurements. How-ever, the estimated abundances and extracted endmembers do not necessarily optimally match the original data. In [18], an improved version was proposed. In this method, two objective functions are established, one for spatial measurements and one for spectral measurements. Then, the endmembers and abundances are estimated by continuous iterations. However, it is difficult to optimize the two objective functions at the same time. In other words, when one objective function converges to the optimal point, the other objective function is not necessarily near the optimal point and may even deviate significantly from it. In addition, it is very difficult to measure the contribution of the optimal values of two independent objective functions. In practical application, it is very difficult to adjust the degree of contribution of the optimal values of two independent objective functions.
In view of the abovementioned non-uniformity issue with two objective functions, this article presents a single objective function to optimize both the spatial measurements and compressed spectral data. This algorithm is referred to as spatialspectral joint CS (SSJCS). In addition, inspired by sparse unmixing [36]- [38], a row sparsity constraint is introduced for the abundances in CS reconstruction to further improve the reconstruction accuracy. Another contribution of this article is the design of an effective algorithm for addressing the aforementioned optimization problem. This algorithm allows rapid and high-accuracy HCS reconstruction based on spectral unmixing.
The remainder of this article is organized as follows. Section II introduces the SSJCS algorithm for HSIs under an LMM and its mathematical expression. Section III details the process of the proposed reconstruction algorithm. Section IV presents experimental comparisons of SSJCS and several state-of-the-art reconstruction algorithms based on several widely used HSIs and demonstrates the excellent performance of SSJCS. Finally, Section V provides the conclusions of the paper.

II. SPATIAL-SPECTRAL HYBRID COMPRESSED SAMPLING
X ∈ R L×N is used to represent an HSI, where L is the number of bands and N is the number of pixels. Here, the spatial image of each band is pulled into one row. In other words, each row of X represents the greyscale image of any band, and each column of X represents the spectral curve of any pixel. Thus, the famous LMM can be represented by the following equation: where E ∈ R L×p represents the endmember matrix and S ∈ R p×N represents the corresponding abundance matrix. p is the number of types of ground objects contained in the scene and is also referred to as the number of endmembers. In the spectral unmixing field, two important physical constraints on an LMM, namely, the abundance non-negativity constraint and abundance sum-to-one constraint, generally cannot be ignored. However, in CS reconstruction, it is found that these two constraints have no significant impact on the reconstruction performance but increase the computational complexity of the reconstruction. Therefore, in this article, these two physical constraints are not taken into consideration. SSJCS divides hyperspectral data acquisition into two parts, namely, spectral and spatial CS. These two parts can be realized in parallel using two lenses or in sequence using a single lens [17]. The mathematical expression of spectral CS is as follows: where A L ∈ R l×L is the spectral CS matrix (a uniformly distributed random matrix can be used as A L in experiments). l L means that the number of observation bands and is far smaller than that of original bands, which implies that spectral compressive is achieved. Y L ∈ R l×N is the spectral CS measurements. Then, the spectral sr can be defined as The realized form of spatial CS is as follows: where A N ∈ R N ×n is the spatial CS matrix. n N means that the number of observation pixels and is far smaller than that of original pixels, which implies that spatial compressive is achieved. Y N ∈ R L×n represents the spatial CS measurements. Thus, the spatial sr can be defined as sr N = Ln/LN = n/N . It is worth noting that to facilitate the extraction of initial endmembers, the measurement matrix designed by Wang et al. [14] is still used as A N . In other words, each column of A N is a one-hot vector. This vector has only one element with a value of 1, and each of its other elements has a value of 0. In addition, to ensure the randomness of the measurement matrix, an element with a value of 1 can be randomly placed at any position. Y L and Y N are the actual data acquired by SSJCS. Thus, the total equivalent sr is the sum of sr L and sr N , i.e., sr = sr L + sr N . Considering that the number of bands is often far smaller than that of pixels in an HSIs, we let sr L sr N when setting the srs to ensure that a sufficient number of wavenumber dimensions can be acquired. In the subsequent design, sr N is set to a fixed value of 0.01. The total sr is varied by altering sr L .

III. SPATIAL-SPECTRAL JOINT RECONSTRUCTION
The goal of SSJCS reconstruction is to recover the original data X at high accuracy based on the acquired compressed data Y L and Y N under the given A L and A N . To achieve this, the following objective function is established to solve the reconstruction problem: where C F ≡ trace CC T is referred to as the Frobenius norm of C, the subscript T in C T denotes the transpose of the matrix, S 2,1 ≡ p i=1 s i 2 is defined as the 2,1 norm of S [38], [39] to describe the row sparsity of the abundance, and α and β are regular-term parameters that are used to weight the second and third terms in the objective functions, respectively. In Equation (4), the spectral and spatial observations are combined under one objective function, and their weights are adjusted through α. In addition, the row sparsity constraint on the abundances is increased, which can help E and S approach their respective optimal values. The optimization of Equation (4) requires the simultaneous minimization of E and S and is a nonconvex optimization problem, which is difficult to solve directly. Thus, the alternating direction method of multipliers (ADMM) [40] is used to optimize Equation (4). First, auxiliary matrices Thus, the unconstrained optimization problem in Equation (4) can be converted to the following constrained optimization problem: Equation (5) can be simplified to the following compact form: where I is the identity matrix. Equation (6) is a more compact constrained optimization problem than equation (5). The matrices G and B are constant matrices created by the given measurement matrices A L , A N and identity matrix I. The variable matrix R is composed of matrices R 1 , R 2 , R 3 , R 4 VOLUME 8, 2020 Algorithm 1 ADMM pseudocode for solving Equation (6) 1.
Until some stopping criterion is satisfied. and R 5 , and T is composed of endmember matrix, abundance matrix and zero matrix. g (R) corresponds to the objective function 1 2 (6) is to minimize the objective function g (R) by the variables R and T under the constraint of GT + BR = 0.
Thus, the augmented Lagrangian in Equation (6) can be expressed as follows: where µ > 0 is a positive constant used to control the rate of convergence of the algorithm and D represents Lagrange multipliers. In each iteration of the ADMM, Algorithm 1 sequentially executes the optimization of T by L (Step 3), the optimization of R by L (Step 4), and the update of the Lagrange multipliers (Step 5).
After the k th iteration, the following residuals are defined: where resp and resd are primal and dual residuals, respectively. Boyd et al. [41] suggest that a reasonable termination criterion is that the primal and dual residuals must be small, i.e., resp ≤ ε p and resd ≤ ε d (10) where ε p > 0 and ε d > 0 are feasibility tolerances, respectively. The same tolerance is used for ε p and ε d .
where ε 0 = 10 −4 is an absolute tolerance. The update equation of the ADMM shows that a relatively high µ value generates a relatively large penalty for violations of primal feasibility and therefore tends to generate relatively small primal residuals. Conversely, the definition of resd suggests that a relatively low µ value tends to reduce dual residuals, albeit at the expense of reducing the penalty for primal feasibility. This may result in relatively large primal residuals. Therefore, the initial value of µ (µ 0 ) is set to 0.5, and µ is continuously adjusted based on resp and resd during each iteration.
Considering that the ADMM may converge to a local optimal solution, the selected µ 0 is vitally important for the final result. The special structure of A N renders it easy to estimate the endmember E 0 from Y N . Research on linear spectral unmixing shows that E 0 is very close to the final endmember. E 0 is extracted by the following procedure. The hyperspectral signal identification by minimum error (HySime) algorithm [42] is first employed to estimate the number of endmembers p from Y N . Then, the vertex component analysis (VCA) algorithm [43] is used to extract E 0 from Y N , i.e., where vca signifies the VCA algorithm. The initial abundance can be estimated based on Y L using the method of least squares: where C −1 is the inverse matrix of C. In Algorithm 1, T 0 and R 0 can be initialized using E 0 and S 0 based on equation (13) and (14), and D 0 is initialized to an all-zero matrix. The details of the optimizations with respect to T and R of ADMM Algorithm 1 (SSJCS) are shown in the Appendix at the end of the paper.

IV. EXPERIMENTAL RESULTS AND DISCUSSION
In this section, the effectiveness of the proposed SSJCS algorithm is examined through comparison experiments based on four real-world HSIs. The SSJCS algorithm is compared with CPPCA (a spectral CS algorithm) [9], spectral compressive acquisition (SpeCA) (a spatial-spectral dual sampling algorithm) [17], and spatial-spectral compressed reconstruction based on spectral unmixing (SSCR_SU) (a spatial-spectral dual sampling algorithm) [18]. For CPPCA, it is only necessary to set one parameter, namely, the number of eigenvectors. Default settings can be used for all the other parameters. According to Fowler's recommendation [9], a value in the range of 3-30 that can realize optimal reconstruction performance is selected as the number of eigenvectors. SpeCA and SSCR_SU reconstruct data by spectral unmixing. SSCR_SU uses the HySime algorithm to estimate the number of endmembers, whereas SpeCA is unable to estimate the number of endmembers in advance. Therefore, the number of endmembers is set based on the real number of types of ground objects during the experiments.

A. PARAMETER SETTINGS
The proposed SSJCS algorithm is evidently sensitive to the settings of the regular-term parameters α and β. Unavoidably, in the experiments, it is necessary to give a suitable range for parameter settings. Thus, an attempt is first made to find stable and effective parameter settings through a set of experiments. The Samson dataset used in the experiments is a simple dataset 1 available online. This dataset contains 952 × 952 pixels of 156 channels with wavelengths between 401 and 889 nm. In this experiment, a small area (95 × 95) containing three types of ground objects (soil, trees, and water) is extracted for testing. Figure 4(a) shows a red-greenblue (RGB) pseudo-colour map of the Samson dataset (the red, green, and blue channels correspond to the 105, 50, and 18 bands, respectively). The variation in the reconstruction performance of SSJCS with α and β is measured using the average signal-to-noise ratio (SNR), which is defined as follows: where X i andX i are the original and reconstructed images of the band i, respectively. Figure 1 shows the effects of α and β on the average SNR of the reconstructed image at various srs. As demonstrated in Figure 1, the variation in the SNR with each of α and β exhibits basically similar trends at various srs. As α gradually increases, the SNR first increases and then decreases. As β decreases, the SNR similarly first increases and then decreases. This is because the parameters α and β in the objective function (4) weigh the fidelity term of spatial compressed data and the sparsity term of abundance respectively. When α increases, the proportion of fidelity terms in the objective function increases, and when β decreases, the proportion of sparse terms decreases. The change trend of SNR in Figure 1 shows that the appropriate values of these two parameters are conducive to further improve the reconstruction accuracy. However, at a relatively large β (e.g., β = 1), the SNR is more sensitive to α. As α increases, the SNR increases to its peak value and subsequently decreases. However, as shown in Figure 1, at a relatively small β, the SNR decreases slowly as α increases. This results in a relatively large flat area on the curved surface in Figure 1. The SNR in this area is close to the peak value. This facilitates the setting of α and β. Therefore, α and β are set to 100 and 10 −4 , respectively, in the subsequent experiments. While these settings do not necessarily ensure the maximum SNR at various srs, they can ensure that the SNR is near its maximum value. Figure 2 shows the variation in the residuals of the iteration. As demonstrated in Figure 2, as the number of iterations increases, the primal and dual residuals decrease rapidly to relatively small values. After 30 iterations, as the number of iterations increases, there are relatively insignificant changes in the primal and dual residuals. This suggests that the algorithm has converged to a relatively optimal solution and that further increasing the number of iterations will relatively insignificantly affect the reconstruction performance. Therefore, the maximum number of iterations maxiters is set to 30 in the subsequent experiments.

B. SNR COMPARISON
In this section, the proposed SSJCS algorithm is compared with CPPCA, SpeCA, and SSCR_SU in terms of the SNR of reconstructions based on the Samson, Cuprite, PaviaU, and Urban datasets. In the experiments, sr is increased from 0.1 to 0.5 at a step size of 0.1.
The Cuprite dataset consists of data acquired by the Airborne Visible/Infrared Imaging Spectrometer on June 19, 1997 over a mining area in southern Florida [43] covered primarily by some minerals and a small amount of vegetation. A sub-image containing 250 × 191 pixels is used in the experiment. This image contains a total of 224 bands between 400 and 2,500 nm. After removing the noise and water absorption bands, the image of the remaining 188 bands is used in the experiment. Due to the minor spectral difference between some similar minerals, it is generally considered that this area contains 12 different types of ground objects. Figure 4(a) shows an RGB pseudo-colour image of the Cuprite dataset (the red, green, and blue channels correspond to 40, 20, and 10 bands, respectively).
The Urban dataset, often used in HS image processing, was acquired by the HS Digital Imagery Collection Experiment sensor over an urban district of Copperas Cove, Texas, USA. This dataset contains 210 bands with wavelengths between 400 and 2,500 nm and has a spatial resolution of 307 × 307. After removing the water and noise absorption bands, the remaining 162 bands are used in the comparison experiment. This area contains six types of ground objects. Figure 4(a) shows an RGB pseudo-colour image of the Urban dataset (the red, green, and blue channels correspond to 28, 11, and 2 bands, respectively).
The PaviaU dataset was acquired by the Reflective Optics System Imaging Spectrometer over Pavia, northern Italy. After removing 12 noise bands, this dataset contains 610 × 340 imagery data points in 103 bands and has a geometric ground resolution of 1.3 m. In addition, this dataset contains nine types of ground objects. Figure 4(a) shows an RGB pseudo-colour image of the PaviaU dataset (the red, green, and blue channels correspond to 50, 30, and 5 bands, respectively). Figure 3 shows the reconstruction performance of the algorithms on the four datasets at various srs. The figure shows that the SNR curve of the image reconstructed by each algorithm on each dataset increasing with increasing sr. The SNR curves corresponding to the spectral unmixingbased reconstruction algorithms are flatter, whereas the SNR curve corresponding to CPPCA increases more rapidly as the sr increases. This means that the spectral unmixing-based reconstruction algorithms exhibit higher performance at low srs. This is corroborated by the results in Figure 3(a) and (c). The SNRs corresponding to SpeCA and SSCR_SU are higher than those corresponding to CPPCA at low srs but lower than those corresponding to CPPCA at higher srs. However, SSJCS demonstrates the optimal performance on all of the datasets. For some datasets, the SNRs corresponding to SSJCS are far higher than those corresponding to the other algorithms. For example, at an sr of 0.2, the SNR corresponding to SSJCS is higher than those corresponding to the other algorithms by nearly 10 dB for the Samson dataset. At an  sr of 0.5, the SNR corresponding to SSJCS is higher than those corresponding to the other algorithms by 8 dB. This sufficiently demonstrates the effectiveness of the proposed SSJCS algorithm. SSCR_SU and SpeCA exhibit similar performance on the four datasets. In most cases, SSCR_SU slightly outperforms SpeCA. It is worth noting that CPPCA exhibits poor reconstruction performance on all the datasets at srs below 0.1. This suggests that CPPCA fails at a srs below 0.1.
The images reconstructed by SSJCS based on the PaviaU dataset are compared with those reconstructed by CPPCA reported by Chen et al. [32] and the relevant improved algorithms (class-dependent CPPCA (C-CPPCA) [31], multiple hypotheses (MH)(uniform grouping (U))-CPPCA, MH(non-uniform grouping (NU))-CPPCA, MH(U)-C-CPPCA, MH(NU)-C-CPPCA [32], and multi-task Bayesian CS (MT-BCS) [44]). Table 1 summarizes the results. As demonstrated in Table 1, SSJCS far outperforms CPPCA and the relevant improved algorithms in reconstruction. At an sr of 0.2, the SNR of the image reconstructed by SSJCS is higher than that of the image reconstructed by MH(NU)-C-CPPCA by more than 14 dB. This adequately demonstrates the high efficiency of the proposed algorithm. In addition, by comparing the reconstruction by CPPCA in Figure 3(d) and Table 1, it is found that the reconstruction performance of CPPCA is higher. This is because when running the CPPCA code, the key parameter-the assumed number of eigenvectors-was adjusted to its optimum.
To further confirm the effectiveness of the proposed algorithm, we compared the mean peak signal-to-noise ratio (MPSNR) index on the PaviaU and Urban datasets with several state-of-the-art reconstruction algorithms, which is defined in [25]. The comparison algorithm includes our previous work, hyperspectral compressed sensing based on a multitype mixing model (HCS_MMM) and distributed compressed hyperspectral sensing (DCHS), as well as joint tensor decomposition and reweighted TV regularization method (JTenRe3-DTV) [23], nonlocal tensor sparse and low-rank regularization (NTSRLR) [25], and subspace-based nonlocal tensor ring decomposition method (SNLTR) [28]. It should be noted that the PaviaU and Urban datasets have been cropped into subimages of 300×300×103 and 200×200×184 in references [25] and [28], respectively. The experimental results of all algorithms in Table 2 and Table 3 are for the cropped subimages. The proposed algorithm additionally adds the reconstructed MPSNR of the whole area without cropping. The cropped and cropped results are marked in Table 2 and Table 3 with the dimensions of the dataset, respectively. VOLUME 8, 2020  For PaviaU dataset, all algorithms have poor reconstruction performance at 0.02 sampling rate, although our algorithm achieves the best reconstruction results on uncropped dataset. The sampling rate is increased from 0.02 to 0.05, and the reconstruction accuracy of all algorithms has been greatly improved. For example, the MPSNR of HCS_MMM has even increased by nearly 13dB. At lower sampling rates, DCHS performs better than SSJCS since DCHS is specially designed for low sampling rate environments. In fact, SSJCS can also perform better in low sampling rate environment if the appropriate prior constraint of abundance is added, but this is beyond the scope of this article. At a sampling rate above 0.15, SSJCS can achieve the best reconstruction performance better than other algorithms. Comparing the performance of SSJCS on cropped and uncropped datasets, the cropping of the data seems to have no effect on the algorithm. The reconstruction accuracy of the two is basically close, except for the 0.02 sampling rate. This is because the PaviaU dataset only crops spatial pixels without changing the spectrum. Moreover, the spatial resolution is still high after cropping, which has little effect on the algorithm, and even the reconstruction performance is improved at some specific sampling rates.
For the Urban dataset, SNLTR obtains the best reconstruction performance under low sampling. At higher sampling rates, SSJCS performs better on the uncropped dataset. Here, the performance of SSJCS in cropped and uncropped datasets is quite different. The reconstruction result of the uncropped dataset is obviously better than the cropped dataset. This is mainly caused by two reasons. One is that the cropped image is too small, which affects the reconstruction performance. The other is that the cropped dataset retains 184 bands from the original 210 bands. A small amount of noise is still included in some bands, which seriously degrades the accuracy of LMM. Figure 4 shows the original pseudo-colour images of the four datasets as well as pseudo-colour images reconstructed by various algorithms at a sr of 0.2. It is very difficult to determine the differences between the details of the images reconstructed by the algorithms. This is because CPPCA and the other three spectral unmixing-based reconstruction algorithms use spectral CS measurements. Reconstruction of spectral CS measurements allows relatively complete preservation of spatial information. Therefore, the images reconstructed by these algorithms are very close to the original image in terms of spatial details.
However, for the PaviaU dataset, it can still be found that the image reconstructed by CPPCA is relatively dark, whereas the images reconstructed by the other three spectral unmixing-based reconstruction algorithms are also slightly darker than the original image but closer to the original image than the image reconstructed by CPPCA. This suggests that the performance of these algorithms in recovering spectral information is slightly lower than that in recovering spatial information. CPPCA performs even more poorly than the other three algorithms. SpeCA, SSCR_SU, and SSJCS are similar in terms of reconstruction performance. The minors differences between the images reconstructed by these algorithms cannot be distinguished visually.
In order to clearly distinguish the reconstruction results of different algorithms visually, Figure 5 provides the residual images between the original images and the reconstructed images achieved by different algorithm with 0.5 sampling rate. When the sampling rate is 0.5, the performance of several algorithms on the Cuprite dataset is similar, so the brightness of their residual images is also similar. However, on the other three data sets, the residual images of the proposed SSJCS algorithm are darker than the other three algorithms. This means that SSJCS has the lowest residual and the best reconstruction performance. We can also see from the residual images that CPPCA performs better than SpeCA and SSCR_SU at 0.5 sampling rate on Samson and Urban datasets.

C. SAM COMPARISON
In this section, the abovementioned four reconstruction algorithms are compared using the average spectral angle mapper (SAM), which is defined below: where X j andX j are the original and reconstructed spectral curves of the j th pixel, respectively. The lower the SAM is, the closer the reconstructed spectral curve is to the original spectral curve. Table 4 summarizes the average SAMs obtained by the algorithms when they are used to reconstruct data based on the selected datasets at various srs. The results in Table 4 are similar to those in Figure 3. The proposed SSJCS algorithm exhibits no notable advantages on the Cuprite dataset but is significantly advantageous to the other three algorithms on the other three datasets. In some cases, the average SAM of reconstruction by SSJCS is even less than half that of reconstruction by the best of the other three algorithms. At a sr of 0.5, the SAM achieved by CPPCA is lower than that achieved by SpeCA on each of the four datasets, is far lower than that achieved by SSCR_SU on the Samson and Urban datasets and is close to that achieved by SSCR_SU on Cuprite and PaviaU datasets. This suggests that at high srs, CPPCA outperforms SpeCA and SSCR_SU in recovering spectral information. However, the performance of CPPCA is still slightly lower than that of the proposed algorithm, particularly on the PaviaU dataset. This further demonstrates that the proposed algorithm is effective in reconstructing spectral information. In addition, it is noted that SpeCA exhibits the highest reconstruction performance on the PaviaU dataset at VOLUME 8, 2020 a sr of 0.1. However, this does not undermine the performance of SSJCS. Under this condition, the SAM achieved by SSJCS is higher than the optimal value by only 0.02. This suggests that SSJCS still exhibits excellent reconstruction performance.
To more visually observe the recovery of spectral information, specific pixels in the four datasets at a sr of 0.5 are selected, and the original spectral curves and the spectral curves reconstructed by various algorithms are plotted, as shown in Figure 6. Specifically, the pixel at (50, 50) in the Samson dataset is selected, while the pixel at (100, 100) in each of the other three datasets is selected. In addition, to facilitate observation, the original spectral curves are enlarged locally. A comparison of the locally enlarged spectral curves shows that the spectral curve reconstructed by SSJCS based on each of the four datasets relatively satisfactorily approaches the original spectral curve. While SSJCS performs somewhat poorly on the PaviaU dataset, it still significantly outperforms the other three algorithms. All the algorithms exhibit relatively excellent reconstruction performance on the Cuprite dataset. This is consistent with the results in Table 4. The spectral curves reconstructed by SpeCA and SSCR_SU based on the Samson and Urban datasets deviate relatively significantly from the original curve, whereas SSJCS and CPPCA are able to recover the original spectral curve with relatively low distortion.

D. SSIM COMPARISON
The final evaluation index used in the experiments is the structural similarity index (SSIM), which is an important index for measuring the structural consistency between the original and reconstructed images. The average SSIM is defined as follows: where µ o and µ r are the average values of the original and reconstructed images X i andX i of the i th band, respectively, σ o and σ r are the variances of X i andX i , respectively, σ or is the covariance of X i andX i , and c 1 = (k 1 ) 2 and c 2 = (k 2 ) 2 are constants used to maintain stability (k 1 = 0.01, k 2 = 0.03, and is the dynamic range of the pixel value) [45]. The higher the SSIM is, the more structurally similar the reconstructed image is to the original image. Table 5 summarizes the average SSIMs of reconstruction by various algorithms. As demonstrated in Table 5, the proposed SSJCS algorithm exhibits the best performance on all three datasets except for the Urban dataset. SSCR_SU outperforms SSJCS on the Urban dataset. For the Urban dataset, the average SSIM achieved by SSJCS is again the highest only at a sr of 0.5. This suggests that while SSCR_SU performs moderately in conventional indices (SNR and SAM), this algorithm exhibits excellent performance in structural similarity. Nevertheless, SSJCS still outperforms SSCR_SU on the other datasets.
In summary, the proposed SSJCS algorithm outperforms the other state-of-the-art algorithms in both spatial information recovery and spectral curve reconstruction. SSJCS also exhibits excellent reconstruction performance even at the VOLUME 8, 2020  structural similarity. All the experiments demonstrate that the proposed SSJCS algorithm is an effective HCS reconstruction algorithm.

V. CONCLUSION
This article presented an effective HCS reconstruction algorithm based on spatial-spectral hybrid CS data of HSIs, which reconstructed the original HSIs by establishing a uniform objective function from the spatial-spectral joint measurements. In addition, a row sparsity constraint is introduced for abundances to improve reconstruction accuracy. An approach is developed based on the augmented Lagrangian method and alternating minimization. This approach effectively addresses the reconstruction problem of the objective function. Experimental comparisons in multiple indices show that the proposed algorithm notably outperforms other algorithms in SNR and SAM and exhibits excellent performance in SSIM. In the future, a suitable endmember constraint will be developed for the objective function to further improve reconstruction accuracy.

VI. APPENDIX
In this Appendix, the implementation of the SSJCS algorithms described in Section 3 is explained in detail. First, let us begin with the varying augmented Lagrangian in Equation (8).
While keeping the other variables unchanged, the minimization of E by the augmented Lagrangian can be expressed as follows: arg min which has a closed-form solution of where I L is an L × L identity matrix. Similarly, to calculate S k+1 , S is optimized by the augmented Lagrangian as follows: Its solution is as follows: where I N is an N × N identity matrix. Under normal circumstances, the spatial resolution of hyperspectral (HS) images is relatively high, i.e., N is a relatively large number. As a result, the dimensionality of A N A T N and I N is very large. Fortunately, each column of the designed A N corresponds to a one-hot vector. Thus, A N A T N and I N are both sparse matrices. These matrices do not occupy a large amount of memory. In addition, inverting A N A T N and I N does not involve a high computational load.
The R 1 optimization subproblem can be converted to arg min whose closed-form solution can be obtained as where I p is a p × p identity matrix. The objective function of the R 2 minimization subproblem is as follows: which can be easily deduced as The optimization problem for R 3 can be specified as arg min which has the analytical solution of To compute R 4 , we solve the optimization problem arg min which has a closed-form solution Finally, R 5 is computed by solving the following optimization problem arg min R 5 β R 5  She is currently a Professor with the College of Electronic Science, National University of Defense Technology. Her current research interests include intelligent information processing and information system technology. He is currently an Associate Professor with the School of Biomedical Engineering and Imaging Medicine, Army Medical University. His research interests include hyperspectral images processing, pattern recognition, and machine learning.