A blind spectral unmixing in wavelet domain

—In this paper, a wavelet-based energy minimization framework is developed for joint estimation of endmembers and abundances without assuming pure pixels while considering noisy scenario. Spectrally dense and overlapped hyperspectral data is represented using biorthogonal wavelet bases that yield a compact linear mixing model in wavelet (sparse) domain. It acts as the data term and helps reducing solution space of the unmixed components. Three prior terms are added to better handle ill-posedness, i.e., logarithm of determinant volume regularizer enforces minimum endmember simplex, smoothness (spatial) prior on individual abundance maps while spectral constraint through learning dictionary of abundances. Alternating nonnegative least-squares is employed to optimize the regularized matrix factorization functional in wavelet domain. We conduct theoretical analysis, discuss convergence and algorithmic details. Experiments are conducted on synthetic and three real benchmark hyperspectral data AVIRIS Cuprite, HYDICE Urban, and AVIRIS Jasper Ridge. The efﬁcacy of proposed algorithm is evaluated by comparing results with state-of-the-art.


I. INTRODUCTION
Hyperspectral imaging [1] offers a set of co-registered images with a greater spectral resolution typically 10nm at the cost of poor spatial (ground) resolution typically 1m 2 to 30m 2 . Due to the poor ground resolution, reflectance from multiple materials (objects) fall within a single pixel location and present so-called mixed-pixels accross the hyperspectral data. This needs to be unmixed in order to better understand and quantitatively analyze the remotely acquired scene.
Spectral unmixing determines number of spectrally distinct materials signatures, extracts their signature values called endmembers, and estimates corresponding fractional contributions called abundances [2]. Considering passive remote sensing, endmembers are nonnegative while abundances are often constrained to nonnegative and sum-to-one at each pixel location in a scene. It is highly active research area since the resultant material (abundance) maps allow quantitative analysis of the scene especially of physically inaccessible sites. The signatures in hyperspectral data are statistically dependent, and the observed mixtures in the data are either linear or non-linear in nature [3]. Note that the mixing characteristics, high dimensionality of data and the physical constraints on the unmixed components, place the hyperspectral unmixing beyond the reach of most source separation algorithms [4].
Joint estimation of endmembers and abundances is a challenging problem without assuming pure pixels and under noisy scenario. To this end, the wavelet transform with its inherent multiresolution and discriminating capability towards signal analysis can be sought. In this paper, we employ wavelet based features to jointly unmix the data. A novel blind hyperspectral unmixing method is developed by proposing a compact LMM in the sparse wavelet space (data term) along with total of three constraints (prior terms) on both endmembers and abundances in the wavelet domain. The compact LMM using biorthogonal wavelet yields a most compressible wavelet bases that help improving spectral unmixing process. To our knowledge, this is the first time that constraints for endmembers as well as for the abundances are jointly considered in the wavelet domain to tackle the spectral unmixing without pure pixel assumption in noisy scenario.
This paper is arranged as per the following sections. We first summarize the related works in section II. In section III, we begin by presenting compact LMM and formulating the problem. The proposed approach is then discussed in section IV including optimization, convergence details and theoretical analysis. Section V further validates our approach by conducting experiments on a synthetic and three real hyperspectral data sets. Finally, section VI concludes the paper with possible future directions.

II. RELATED WORK
Over the years, several approaches for estimating number of endmembers, endmember signatures and abundances have been proposed [2], [5], [6], [7]. Blind spectral unmixing algorithms without pure pixel assumption are essentially needed for better adaptability in many real cases [5]. To this end, structured matrix factorization model such as nonnegative matrix factorization (NMF) based unmixing is found suitable with its non-negativity and unsupervised nature [6]. However, note that even after applying NMF, the joint estimation problem remains ill-posed and it is necessary to introduce additional constraints to make it a better-posed [8].
In [9], authors presented an idea providing different algorithms for minimizing the objective function that calculates minimum-volume enclosing simplex. It is a cyclic minimization procedure in which a sequence of linear programs is solved. Notably such algorithms are generally developed under a noiseless setting. The algorithms [10] and [11] are recent variants of such algorithm. However, these are computationally prohibitive and may not converge. Besides, such algorithms implement a robust version of the minimum volume that violates the nonnegativity constraint on abundances. In [12], volume minimization regularizer is presented, however, it is sensitive to outliers and also comes with heavy computation burden. To address such issues, a noise-robust version of the volume regularizer [13] is recently derived, which adopts a modified log-determinant (log-det) loss function. The works in [14] and [15] use volume minimization NMF with successive convex optimization and alternating optimization, respectively. The major drawback in these algorithms is that they are developed under noiseless setting, and thus only work well for datasets with a high signal-to-noise ratio (SNR). More recently, a comparison of the three forms of volume regularizers, namely determinant, log-det, and nuclear, is discussed in [16]. It is found that for highly mixed data sets, the logdet regularizer exhibited better results. Recently an approach is proposed that simultaneously counts and extracts endmembers [17] using a new volume calculation formula. Authors have introduced a new way to calculate a volume for developing a real-time implementation of maximum simplex volume algorithm to extract endmembers [18].
The algorithms to estimate abundances are mainly based on the least-squares, geometry-based and sparse regression approaches [19], [20], [21]. The ill-posedness is handled by use of data-driven regularization that helps making it a better posed [22]. We recently proposed wavelet based approach for estimating abundances [23]. Similarly, curvelet domain-based NMF for abundance estimation is also recently proposed in [24].

III. A COMPACT LMM FOR HYPERSPECTRAL DATA MODEL
Hyperspectral sensors acquire hundreds of contiguous electromagnetic bands of narrow bandwidth (10nm) [1]. To this end, we observe that the spectrally dense and overlapped hyperspectral data is compressible in the wavelet space. We found that biorthogonal family of wavelets best characterize the hyperspectral data by capturing the variations compactly and completely. Fig. 1 (a) shows the biorthogonal wavelets and scaling function approximations for both decomposition and reconstruction. Fig. 1 (b) illustrates three benchmark real hyperspectral data AVIRIS Cuprite [25], HYDICE Urban [26], and AVIRIS Jasper Ridge [27] represented using the biorthogonal discrete wavelet transform (DWT). One can see from Fig. 1 (b) that magnitudes of the resultant DWT coefficients for all the data sets are decaying as per the power law [28]. This validates the fact that hyperspectral data is compressible in wavelet space and hence it can be better represented using biorthogonal wavelet. Since DWT is an orthonormal transform, it conveniently preserves the structure of LMM. Therefore, we resort to represent the hyperspectral data in the wavelet domain using biorthogonal bases. The DWT of hyperspectral data vector r is an inner product defined as, where, <> is inner product operator, ψ s,b are referred to as wavelet bases, and r w denotes resultant DWT coefficients of the data r. Hence, using properties of DWT, we can conveniently re-write the well adopted LMM as, Here, size of data matrix R w is K × N, E w represents the noise matrix of size K × N, for e number of endmembers, size of endmembf endmeer matrix in wavelet domain is M w is K × e, and A w is the abundance matrix of size e × N.
Hence, it can be seen that Eq. (3) is a compact LMM for hyperspectral data in wavelet domain. For hyperspectral unmixing, first step is to find number of endmembers which is either known a priori or to be estimated from available data using virtual dimensionality (VD) [29], [30]. In this paper, we find the number of endmembers using our recently developed multiple hypothesis based approach [31]. This controls the incurrent of Type I error while estimating the number of endmembers using the VD.
Problem statement: The hyperspectral data in wavelet domain is now considered as a data cube, x × y × K, where x and y represent spatial dimensions, and K represents the wavelet spectral dimension of the remotely acquired scene. Given N such hyperspectral data vectors in wavelet space Eq.
(3), our objective is joint estimation of endmember signatures and corresponding abundance maps of the scene.

ABUNDANCES
In this section, we propose blind hyperspectral unmixing in wavelet domain for the defined problem in section III. Given theê using our recent work [31] and referring to the data model (compact LMM) defined in wavelet domain Eq. (3), we now propose to solve the simultaneous endmember extraction and abundance estimation problem in an optimization framework by minimizing the generalized Euclidean distance between R w and M w A w , is a e N -dimensional column vector of all ones. Physical constraints are involved in estimation of endmembers (nonnegative) and abundances (nonnegative and sum-to-one). Here, · F denotes the Frobenius norm. It can be seen from Eq. (4) that the joint estimation of endmembers and abundances is severely ill-posed problem due to the correlated signatures, variations among abundances as well as presence of noise and outliers in the data. Though representing the data in wavelet space helps reducing solution space of unmixed components (Eq. (3)). Further use of data-driven regularization with appropriate priors is useful in better handling the illposedness [32]. In this work, we incorporate constraints on  [25], HYDICE Urban [26], and AVIRIS Jasper Ridge [27] represented using bior wavelet at optimum nodes of decomposition. It can be observed that the wavelet coeffcients for the three real hyperspectral data sets decay as per the power law [28].
endmembers and abundances, and Eq. (4) is now expanded to, Here J 1 (M w ) is penalty function calculating the simplex volume determined by the estimated endmembers, J 2 (A w ) is the smoothness constraint (spatial) on estimated abundance maps, and J 3 (A w ) is the sparseness constraint (spectral) on resultant abundance vectors. The λ 1 , λ 2 and λ 3 are the regularization parameters to be estimated.
We employ log-det function for calculating the simplex volume determined by the estimated endmembers. It effectively restricts the data within a minimum volume simplex whose vertices represent the compact wavelet bases.
where M w is the endmember matrix in wavelet domain, I e is the identity matrix of size e × e, and δ > 0 is a constant to lower bound V logdet . We include Eq. (6) as J 1 (A w ) prior term in Eq. (5). Now, for better abundance estimation, we consider two additional constraints on the abundances, apart from nonngativity and sum-to-one. First, we employ smoothness constraint on abundances, where L is matix representation of the derivative operator. The abundances are often (spatially) smooth in their respective abundance maps, i.e., distribution of fractions of an individual material is found locally homogenous in the scene. Therefore, we impose smoothness constraint over individual abundance map Eq. (7), and add it as the prior term J 2 (A w ) in Eq. (5) [23].
It is found that number of endmembers in a scene is far less than number of bands in the data. Further, a hyperspec-tral image pixel is generally composed of relatively lesser number of endmembers than the available endmembers in the scene and/or dimensionality of data [33]. Hence considering unmixing as a combinatorial problem, efficient techniques based on sparsity-induced regularizers have been investigated. Motivated with this, we consider the sparseness constraint J 3 (A w ) spectrally on abundance vectors. This is done through a learning dictionary that leads to a best representation under strict sparsity constraints [34]. The update of the dictionary columns is combined with an update of the sparse representations thereby accelerating the convergence. In this work, the dictionary of abundances is learnt as follows, s.t.: A w(x,y,e) ≤ T ∀ (x, y), where T is number of non-zeros in A w , i.e., sparseness. Here, A init w is the abundances calculated using FCLS [35] for the purpose of initialization.
Substituting Eq. (6), Eq. (7), and Eq. (8) into Eq. (5), the final cost function of the proposed joint estimation becomes, Finally, estimated endmember matrix M can be obtained from M w by applying the inverse DWT using the biorthogonal bases. A sufficient condition for the reconstruction of M w is, where, <> is inner product operator, · is the dot product operator, ψ s,b form an orthonormal bases, and M denotes resultant endmember matrix. Further, A w is the abundance (coefficient) matrix of size e×N that represents the abundances of the scene. It can be seen that minimizing the final cost function Eq. (9) with respect to both the unknown matrices M w and A w is a combinatorial optimization problem. To this end, we employ alternating non-negative least squares method [36], [37]. Accordingly, Eq. (9) is divided into two sub-problems with following update rule, and where the parameters γ 1 and γ 2 are the step sizes selected based on Armijo rule [38] with k as iteration number. We update endmember matrix M k w , while keeping the abundance matrix A k w fixed, and viceversa. Projected gradient learning [37] imposes the non-negative constraint on both M k w and A k w . If the new estimate does not satisfy the constraint, a projective function is used to project the point, back to the feasible set. The function max() in Eq. (11) and Eq. (12) is used to set the negative components to zero and keep the nonnegative components unchanged. This is consistent with the passive remote sensing used in the hyperspectral imaging.
The first order derivative of Eq. (9) with respect to M w for gradient learning in Eq. (11) is, . (13) The problem is tailored to compute each element of Substituting Eq. (14) into Eq.(9) and we get, Similarly, the first order derivative of Eq. (9) with respect to abundance matrix A w for gradient learning in Eq. (12) is calculated as follows, The convergence of proposed joint estimation is considered in the mean-squared error (MSE) rate sense [39] as follows, where R = M A. Observe that Eq. (17) tracks the MSE norms as the iteration n increases, and has a constant value once it reaches a very small value of γ, i.e., asymptotically constant.

A. Analysis of the proposed approach
In this subsection, we present theoretical analysis of the proposed approach as well as compare it with other approaches. Referring to the LMM, given the data r, estimating e, M, and α, is a challenging inverse ill-posed problem because: 1) multiple solutions exist, 2) ill-posedness due to noise and outliers, 3) endmember signatures are strongly correlated result in ill-conditioned endmember matrix, and 4) abundances are physically constrained (nonnegative and sum-to-one) and, in addition, varied mixing found among the abundances over a scene. Hence, the blind inversion (spectral unmixing) becomes severely ill-posed.
We observe that data is compressible in wavelet space (Fig.  1) and hence we formulate the problem in sparse wavelet domain that effectively reduces the solution space of the unmixed components by the proposed compact LMM (Eq. (3). It changes topology of the enhanced solution space in order to address ill-posedness. The multiple solutions are also due to variations in abundances and noise in the data. Nevertheless, the abundances are often spatially smooth within an abundance map while found sparsely varying location to location across an abundance vector. From this prior knowledge, it is expected that employing smoothness constraint spatially and considering the sparsity constraint spectrally on abundance vectors improve abundance estimation performance. Therefore to achieve a stable solution, we incorporate smoothness J 2 (A w ) and sparseness J 3 (A w ) constraints on abundances (Eq. (9)).
It is found that endmember matrix M is generally illconditioned, i.e., though the signatures are linearly independent but highly correlated, i.e., statistically dependent. To address this, we seek maximally linearly independent bases in biorthogonal DWT space (Eq. (3)). This essentially better aligns the eigen vectors of resulting M w for a convenient projection for improved estimations of endmembers and abundances. To this end, the proposed compact LMM Eq.
(3) promotes linearly more independent and statistically less correlated set of basis vectors. Note that we consider an optimum node of DWT decomposition. For better reconstruction of endmember signatures, we implement proposed algorithm separately on approximation and detail DWT coefficients. This gives two basis matrices with respect to approximation and detail wavelet coefficients of the data. We combine both the basis matrices and apply inverse DWT to reconstruct the endmembers back in the data domain.
We propose to use biorthogonal wavelets for representing the hyperspectral measurements. The most significant feature of biorthogonality is perfect reconstruction [40]. The biorthogonal wavelet gives invertible matrix with perfect reconstruction. We found that biorthogonal wavelets offer locality and compact support while preserving the advantages of orthonormality. The analysis and synthesis wavelets have different number of vanishing moments and regularities. We use the biorthogonal family of wavelets which has greater number of vanishing moments for analysis. Hence, the proposed algorithm in wavelet domain better handles the ill-posedness by representing given linearly independent and highly correlated signatures into a linearly more independent while statistically less correlated set of vectors. Hence, the proposed method handles the ill-conditioned endmember matrix.
The proposed method do not require the pure pixel assumption. MVCNMF [12] uses a volume regularizer which is the squared volume of the convex hull of columns of endmember matrix without the origin. The proposed approach considers the log-det function J 1 (M w ) for calculating simplex volume aims to minimize the rank of endmember matrix considering also the origin in the volume (Eq. (9)). Calculating volume in the wavelet space converts the highly correlated endmember signatures to less correlated signatures making the better conditioned endmember matrix. Rvolmin [13] is a volume minimization technique developed under a noiseless setting, and thus works well for data sets having higher SNR. VRNMF [16] algorithms with different volume regularizers are competitive with volume minimization works like MVCNMF [12] and RVolMin [13]. The proposed approach, unlike VRNMF, considers a more suitable volume regularizer and considered it in the wavelet space so as to improve the stability of the solution as well as show robustness against noise and variations (Eq. (9)). The proposed method uses the optimum regularization parameters by carrying out the sensitivity analysis.

V. EXPERIMENTAL RESULTS
In this section, we discuss the results obtained for the simultaneous estimation of endmembers and abundances using proposed blind algorithm on different hyperspectral data sets. We compare our results with state-of-the-art algorithms.
We first conduct experiments on a synthesized hyperspectral data cube. To assess the noise sensitivity of the proposed algorithm, we have added different levels of white Gaussian noises in the data set and compare the results with state-of-theart approaches. We fix the values of parameters by conducting sensitivity analysis, as well as present average execution times. Next, we show experiments on three well-known real hyperspectral data sets, i.e., the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) at the cuprite mining site area, USA [27], [26], [25], the HYperspectral Digital Imagery Collection Experiment (HYDICE) sensor at the location of Copperas Cove, near Fort Hood, Texas, USA in October 1995 [27], [26], [25]; and Jasper Ridge acquired via the AVIRIS sensor by the Jet Propulsion Laboratory [27], [26], [25]. Initialization of endmembers is done through VCA [41] and abundances are initialized by FCLS [35] in implementing our algorithm.
A. Experiments on synthetic hyperspectral data set 1) Data generation: The LMM is used to synthesize hyperspectral data. We then add white Gaussian noise of different variances to simulate hyperspectral data set with varying SNR. The number of endmembers in synthetic data set is set to five (5). The selected endmembers from the USGS [42] are Asphalt (gds317), Brick (gds350), Fiberglass (gds374), Sheetmetal (gds352) and Vinylplastic (gds372) spread in 431 bands ranging from 400nm − 2500nm. Hence, they form the endmember matrix M of size 431 × 5. The spheric Gaussian field with different parameters are considered for generating abundance maps in the data set. We resort to them since they belong to the standard Gaussian family which are positive definite as well as sufficiently flexible to generate variedly mixed pixels for data. Nonnegative and sum-to-one abundance vectors of size 5 × 1 are randomly generated for every location using the chosen spheric Gaussian field. Note that they contain both smooth regions and edges with significantly high heterogeneous regions making blind estimation process challenging. Using a spheric Gaussian abundance vector α of size 5×1 and the endmember matrix M of size 431×5, a pixel vector r is constructed as r = Mα, i.e., generating a 431 × 1 reflectance data vector. Considering five different abundance maps each of size 128 × 128, a 431-band ground truth data cube is generated, i.e., 128×128×431. The mean image of the synthetic data set along with marked endmembers is shown in Fig. 2.  2) Parameters setting: All the algorithms are implemented in MATLAB (R2019a) platform and run on a laptop with Intel(R) Core(TM) i3-7100 CPU of 2.40 GHz with 4 GB RAM. Independent and identically distributed (IID) Gaussian noises have been added to both the synthetic data sets at different proportions of SNR. Referring to proposed approach as in Eq. (9), it involves three regularization parameters to be estimated. The parameters are fixed based on conducting detail sensitivity analysis at different noise levels. Fig. 3, and Fig. 4 show sample sensitivity analyses of parameter λ 1 for the synthetic data at 20dB SNR for endmember and abundances performance metrics, respectively. Similarly sensitivity analyses of parameters λ 2 , and λ 3 are carried out for fixing optimum values in our experiments. We have chosen optimum values and fixed λ 1 = 0.0001, λ 2 = 0.001, and λ 3 = 0.001 for implementing our algorithm. Similarly optimum values for subproblems Eq. (12) and Eq. (13) are considered as γ 1 = 0.001, and γ 2 = 0.001. One may note that we have used optimum values of various parameters while implementing different state-of-the-art algorithms as mentioned in the respective papers.
3) Performance metrics: We use three metrics to evaluate the performance of our algorithm for estimated endmembers and abundances from the synthetic data sets. The first metric root-mean-squared error (RMSE) [43] measures the meansquared error between true endmember and its estimate. The second metric spectral angle distance (SAD) [44] measures the signature shape similarity between the true endmember and its estimate. See that the SAD [44] is robust against intensity variations at each band unlike the RMSE [43]. Unlike both RMSE and SAD, the third metric spectral information divergence (SID) [45] uses information theoretic measure and considers each pixel vector as a random variable. It measures the probabilistic discrepancy between true and estimated spectra [46].
The ground truth is not available for the real data sets. Hence, in order to perform the quantitative analysis on all the three real data sets, we first find data reconstruction error (DRE) maps in terms of RMSE as shown in [22]. We calculate RMSE value at each pixel location between the available data and reconstructed data using estimated endmembers and abundances by different approaches. For better estimates, we expect lower values of DRE at every location that is quantified by calculating the mean and standard deviation of a DRE map. We also quantify the error with another measure called qualitywith-no-reference (QNR) [47] that generalizes the notion of universal quality index for multiband images. The QNR value is calculated as (1 − d bands ).(1 − d image ), where d bands amounts for the inter-band distortions across the spectral space of data, and d image considers the average distortions over the spatial locations. One can see that the higher values of QNR close to unity indicates better estimates of underlying unmixed components.
4) Results discussion: First, number of endmembers from the synthetic data sets is estimated using our recent MH-NWHFC [31], and the result is listed in Table I. It is clear from Table I that our proposed MH-NWHFC [31] is able to estimate five (5) endmembers from the data at different SNR values. Next, signatures of these five endmembers of the synthetic data set are estimated and shown in Fig. 5. One can notice that the endmembers estimated by the proposed approach (Fig. 5 (e)) are almost the same as the ground truth ( Fig. 5 (a)). Fig. 5 also shows the comparison of the estimated endmembers by different algorithms at 20 dB SNR. The proposed approach allows us to extract features related to scale (or frequency) from the hyperspectral measurements. Further, when individual coefficients are used, it allows the user to extract features related to localized activity within a hyperspectral data vector. The proposed compact LMM and log-det volume regularizer on endmembers in wavelet domain help achieving better performance. Finally, estimated abundance maps are displayed in Fig. 6. One can seen from figure that the estimated maps are visually consistent and the patterns in our estimated abundances are found very close to that of the ground truth abundances and consistent with the state-ofthe-art approaches. Proposed abundances have simultaneously maintained spatial patterns with the smoothness constraint while maintained global (spectral) heterogeneity with sparsity constraint in the wavelet space.
We now quantify the results using the three performance metrics by considering noise sensitivity analysis. Table II and  Table III provides the average error scores of estimated endmembers and estimated abundances, respectively, by proposed approach and compare the results with state-of-the-art methods. The three performance metrics RMSE, SAD, and SID for quantitative assessment are tabulated. We calculate metric values between the estimated and the reference ground truth endmember signatures as well as for abundances listed in Table  II and Table III. It can be seen that the proposed solution better handles the ill-posedness and provides improved estimates of both the endmembers and the abundances.

B. Experiment on real AVIRIS Cuprite data
Here in this subsection, we demonstrate the efficacy of the proposed joint estimation algorithm on a benchmark real hyperspectral data set acquired by the AVIRIS. The data is acquired at the location Cuprite mining site, USA [25]. The data is of 224 contiguous wavelength bands with a range of 370 nm to 2480 nm. The spectral resolution is 10 nm. After removal of noisy and water absorption bands, 188 bands are retained, and a region of 250 × 191 pixels is considered [25]. The mean image of the portion of AVIRIS Cuprite data is shown in Fig. 7 (a). We use bior3.1 that yields a sparse representation of size 250×191×49 whereas the given Cuprite data cube is 250×191×188. We first use our recent work MH-VD [31] to estimate the number of endmembers. As shown in Table IV, we fix it to 12 in our experiment. Fig. 8 shows the estimated endmembers. The endmembers are consistent with respect to state-of-the-art algorithms [13], [16], [41]. Abundance maps estimated using the proposed approach are  [41], (c) MVCNMF [12], (d) VRNMF [16], and (e) proposed.
shown in Fig. 9. These estimates of abundances is consistent when compared to state-of-the-art approach like [48]. Table  V shows the quantitative analysis on all the three datasets. In the table, µ represents mean, and σ represents standard deviation of a DRE map. One can observe from the Table V that better estimate of both endmembers and abundances using the proposed approach decreases the reconstruction error in terms of RMSE for AVIRIS Cuprite data. Further, it is evident from the QNR values in Table V that the data reconstruction using the estimated unmixed component of proposed approach better matches with both spectral and spatial characteristics of given data.

C. Experiment on real HYDICE Urban data
In this subsection, the proposed approach is evaluated on another real data collected by HYDICE. Referring to mean image of HYDICE Urban data in Fig. 7 (b), there are 307×307 pixels. Every pixel is of 2m×2m area. This data has 210 bands with a spectral resolution of 10nm. 162 bands are retained by removing Some bands of dense water vapor and atmospheric effects. First, we use our recent work MH-NWHFC [31] to estimate number of endmembers in this data. Referring to Table IV, number of endmembers is found to be 4. This is consistent to the four land cover types which are generally found in this data set viz, Asphalt, Grass, Tree, Road. Fig.  10 shows that the endmembers estimated using proposed approach and we found them consistent with respect to stateof-the-art approach [16]. Abundance maps corresponding to each endmember are then estimated using our approach and displayed in Fig. 11 which also found visually consistent when compared to the maps obtained using state-of-the-art approach. Table V shows the data reconstruction accuracy in terms of mean and standard deviation of DRE map. QNR values in Table V indicates that the data reconstruction using the proposed approach is better for HYDICE Urban data when compared to results with respect to state-of-the-art methods.

D. Experiment on real AVIRIS Jasper Ridge data
In this subsection, the proposed approach is evaluated on another real data collected by AVIRIS at Jasper Ridge. Each image in the dataset is of 512×614 pixels. The data is remotely acquired in 224 contiguous channels starting from 380 to 2500 nm with up to 9.46 nm spectral resolution within the wavelength range. Referring to mean image of AVIRIS Jasper Ridge data in Fig. 7 (c) we consider a subimage of 100 × 100 pixels. Four land cover types are generally estimated in this data set, viz, Tree, Water, Dirt, and Road. It is found to be 4 referring to Table IV. Endmembers estimated using proposed approach is shown in Fig. 12. Endmembers are consistent when compared to state-of-the-art [16]. Fig. 13 shows the abundance maps corresponding to each endmember which are estimated using our approach and displayed, which are found visually consistent when compared to the maps obtained using state-of-the-art approach [16]. Table V shows the better estimate of both endmembers and abundances for AVIRIS Jasper Ridge using the proposed approach mentioning the mean and standard deviation of the DRE map. Higher QNR values in Table V indicates that the data reconstruction of AVIRIS Jasper Ridge using the proposed approach is better when compared to state-of-the-art methods.

VI. CONCLUSION AND FUTURE WORK
In this paper, we have proposed a novel blind spectral unmixing approach in the wavelet domain without assuming pure pixel and under noisy data setting. Use of biorthogonal wavelet bases is shown for obtaining a compact linear mixing model (LMM) in the wavelet space for the spectrally dense and overlapped hyperspectral data. The compact LMM in wavelet domain has helped reducing solution space of unmixed components. We found that the log-det volume constraint on endmembers in wavelet domain resulted in more accurate endmember estimates. The smoothness prior spatially on the abundance maps while sparsity constraint on abundance vectors (spectrally) help achieving a numerically stable solution for abundances. Theoretical analysis, convergence details, and algorithmic aspects are presented. Qualitative and quantitative result analysis as well as comparison with the state-of-the art algorithms are demonstrated. The experiments validate the efficacy of the proposed algorithm on jointly estimating endmembers and abundances in wavelet domain. We further validate the fact that inclusion of both the spectral and spatial contextual information especially in the wavelet domain while unmixing a scene is better for handling the severe ill-posedness and making the problem better posed.
Future work can be carried as follows: Further optimization techniques can be investigated for more efficient implementation. One could also choose to extend the proposed approach to account for spectral variability in the endmember signatures. Finally, deep learning architectures primarily autoencoders, convolutional neural networks (CNN) and generative adversial networks (GAN) have started showing success in hyperspectral unmixing [49]. It is investigated that the wavelet scattering transform has interesting parallelism with the structure of CNN  Figure 8. Estimated endmembers using proposed approach for AVIRIS Cuprite. Signatures are consistent with state-of-the-art [12], [13], [16], [41]. Asphalt Grass Tree Roof Figure 10. Estimated endmembers using proposed approach for HYDICE Urban. Signatures are consistent with state-of-the-art [16]. [50]. Hence the proposed work could be further explored in development of better CNNs for hyperspectral unmixing. Figure 11. Estimated abundance maps using proposed approach for HYDICE Urban. The visual patterns are consistent with state-of-the-art [16].
Tree Water Dirt Road Figure 12. Estimated endmembers using proposed approach for AVIRIS Jasper Ridge. Signatures are consistent with state-of-the-art [16]. Figure 13. Estimated abundance maps using proposed approach for AVIRIS Jasper Ridge. The visual patterns are consistent with state-of-the-art [16].