SUPER-RESOLUTION OF HYPERSPECTRAL IMAGES USING COMPRESSIVE SENSING BASED APPROACH

Over the past decade hyper spectral (HS) image analysis has turned into one of the most powerful and growing technologies in the field of remote sensing. While HS images cover large area at fine spectral resolution, their spatial resolutions are often too coarse for the use in various applications. Hence improving their resolution has a high payoff. This paper presents a novel approach for super-resolution (SR) of HS images using compressive sensing (CS). Besides ill-posedness of SR problem, the main challenge in HS super-resolution is to preserve spectral contents among all bands while increasing their spatial resolutions. In this work, we first obtain an initial estimate of the super-resolution on a reduced dimension HS data. The HS observations of different wavelengths are represented as linear combination of smaller number of basis image planes (BIPs) using principal component analysis (PCA). The novelty of our approach lies in using CS based approach to super-resolve the most informative PCA transformed image representing highest spectral variance (i.e. the first principal component). Our approach uses low and high spatial resolution dictionaries of patches generated by random sampling of raw patches of PCA transformed images that are generated using the training images having similar statistical properties. Using the sparsity constraint, low resolution test patch is represented as a sparse linear combination of relevant dictionary elements adaptively, that leads to initial estimate of super-resolved PCA image having maximum spectral variability. Since SR is an ill-posed problem, we obtain the final solution using a regularization framework considering the sparse coefficients obtained by the CS approach and the autoregressive (AR) parameters obtained from the initial estimate. The remaining PCA images are up-scaled using regularization, considering the same AR parameters which were obtained from super-resolved PCA image having maximum spectral variability. Experiments are conducted on real HS images collected by Airborne Visible/Infrared Imaging Spectrometer (AVIRIS). Visual inspections and quantitative comparison confirm that our method enhances spatial information without introducing significant spectral distortion.


INTRODUCTION
Hyperspectral sensors collect information in the form of reflectance spectra in very narrow (generally 10 nm wide) contiguous bands simultaneously in the visible to mid infrared portion of the spectrum 0.4-2.5µm.They represent an evolution in technology from multispectral sensors, which typically collect spectral information in only a few discrete, noncontiguous bands (Shaw and Manolakis, 2002).HS images generally contain hundreds of spectral bands (Qian et al., 2009).Being spectrally overdetermined, they provide ample spectral information to identify and differentiate spectrally unique materials (Landgrebe, 2002).The high spectral resolution of hyperspectral sensors preserves important properties of the spectrum (e.g., shape of the narrow absorption bands) and makes possible better discrimination of different materials on the ground (Shaw and Manolakis, 2002).Hence they have been proven to be a powerful source for the monitoring of the Earth surface and the atmosphere on global as well as local scales (Tralli et al., 2005).Nowadays, HS images are widely used in wide range of military and defence applications which include target detection and tracking of objects (Manolakis et al., 2003), agriculture planning (Jusoff and Yusoff, 2008), forest inventory (Goodenough et al., 2004), and urban monitoring (Heiden et al., 2012) to mention a few.These applications require high spectral and high spatial resolution data for an accurate determination of object properties.
In practice, few satellites have high spatial resolution ( < 5m ) sensors available onboard.For example Quickbird satel-lite collects panchromatic (PAN) image of spatial resolution 0.7m and four multispectral bands of 2.8m resolution.Recently launched WorldView−1 and WorldView−2 satellites carry an imaging instrument specifically designed to meet the requirements of very high spatial resolution and more number of spectral bands.World View−1 provides a single PAN image of half a meter resolution, while WorldView−2 provides a high spatial resolution (0.46m) PAN image and eight spectral bands having resolution of 1.84m.Although they have high spatial resolution data, they do not provide better discrimination of many surface materials and environmentally relevant information due to their limited spectral resolution.An improved space borne hyperspectral sensor (e.g.Hyperion and Chiris) provides very large number of narrow spectral bands, but their spatial resolution is very less (17-34m).For airborne hyperspectral sensors (e.g.AVIRIS, HYDICE, HyMap) the spatial resolution is dictated largely by the height of the aircraft (Vane et al., 1993).As the height increases the extent of coverage also increases but this decreases the spatial resolution.At low spatial resolutions, the averaging effect on the pixels degrades the performance of spectral detection algorithm (Gu et al., 2008).
For high spatial resolution, the sensor must have a small Instantaneous Field of View (IFOV).However, this reduces the signal to noise (SNR) ratio.In order to improve the SNR, one has to widen the bandwidth allowing more light to enter into the sensor while acquiring the individual bands.Unfortunately, this reduces the spectral resolution of the sensor.Since there remains a trade-off between the spectral and the spatial resolutions in the hyperspectral sensor (Shaw and hua K. Burke, 2003), we need to perform post processing of the HS images.Super-resolution enhancement refers to an algorithmic approach for increasing the spatial details (Chan et al., 2010).Recently, many researchers have attempted to increase the spatial resolution of the HS images by fusing the PAN image and the hyperspectral data (Winter and Winter, 2001, Hardie et al., 2004, Eismann and Hardie, 2004, Capobianco et al., 2007).However, the limitation of these algorithms is that they require the images to be registered.In this paper we present a novel approach for super-resolution of HS images from the perspective of CS which is not dependent on the registration of the different HS images.Our method needs only HR and LR registered HS images of any compatible scenery to create training dictionaries.This is one time offline procedure.Since it is not always possible to get high spatial resolution data of the same scene, our method provides solution in such circumstances.

THEORETICAL BACKGROUND
Presently, compressive sensing theory is drawing major attention in various image processing applications (Candes et al., 2006)..According to this theory, certain class of images with M pixels need only N nonadaptive samples (N < M ) for faithful recovery (Donoho, 2006).Intuitively this specifies recovery of larger number of pixel values from less number of observations.Thus CS provides solution to inverse problem of recovery of more number of values from less number of observations (samples).Super-resolution is also an inverse problem of similar kind.This motivates us to use CS based approach for super-resolution of HS images.

Compressive Sensing
Compressive sensing presents a new approach in signal processing for sparse signal recovery (Donoho, 2006) which directly measures a compressed representation of the signal.The key assumption is that the most signals that arise in nature are sparse whose digital representation requires few nonzero coefficients.CS exploits this sparsity, by allowing a digital signal to be reconstructed using only few linear measurements when compared to the size of the original signal.As long as the measurement matrix satisfies a Restricted Isometry Property (RIP), the exact signal recovery is possible from these measurements (Cand'es, 2008).
Suppose x is an unknown digital image vector in R N that we want to acquire.Normally this should require N pixels.But if we know apriori that x is compressible in certain transform domain (e.g.wavelet, Fourier), then as per CS theory we can acquire x by measuring only M linear projections rather than all N pixels.If these projections are properly chosen (projection matrix or the measurement matrix satisfy the RIP), the size of M can be smaller than N , the size of image.If the image is K sparse, random projection works if M = O(Klog(N/K)) (Candes and Tao, 2007).Mathematically, under the sparsity assumption, a signal x ∈ R N , can be recovered by solving the l1 minimization using standard optimization tool such as linear programming (Bruckstein et al., 2009) i.e., the problem can be posed as where A fundamental ingredient to deploy (1) in applications is the dictionary φ.There are three different ways to construct the dictionary (Bruckstein et al., 2009): (1) Preconstructed dictionaries, like wavelets (Starck et al., 2007), contourlets (Do and Vetterli, 2005) etc.They are generally used for "cartoon-like" images, assumed to be piecewise smooth having smooth boundaries (Elad et al., 2005, Buades et al., 2010).
(2) Tunable dictionaries, in which a basis or frame is generated under the control of particular parameter (discrete or continuous): wavelet packets (Meyer et al., 2000) (parameter is time-frequency subdivision) or bandelettes (Pennec and Mallat, 2005)(parameter is spatial position).
(3) A training database of signal instances similar to those anticipated in the application, and build an empirically learned dictionary (Elad and Aharon, 2006).Here the entries in dictionary are chosen from the empirical data rather than from some theoretical model.Such a dictionary can then be used in the application as a fixed and redundant dictionary.In our application we explore the third option.

Principal Component Analysis
Hyperspectral images are composed of large number of spectral bands (e.g., AVIRIS acquires 224 bands).Hence, applying superresolution technique to each band separately is prohibitive because of time complexity.In addition individual band SR does not make use of the information present across the bands (Akgun et al., 2005).Information is present across these bands in the form of spectral signatures and the identification of ground materials of interest is based on their unique spectral signatures (Shaw and Manolakis, 2002).Besides the ill-posed nature of SR problem, the HS image SR task becomes difficult since preservation of spectral correlation combined with the SR is more challenging.The spectral content of HS images are inherently low dimensional, hence this must be exploited.PCA plays a central role in the analysis of multivariate data (Jolliffe, 1986).Suppose we have a dataset of B hyperspectral bands with size of each band as M × M pixels.Assume that F represents the pixel vector of size B × 1 along spectral dimension of HS image.A set of principal components of dimension M 2 × K are computed from the first K eigen vectors E=[e1, e2, ...., e K ] T B×K , which in turn is computed from the covariance matrix Σ of size B × B from the given data set where Here m F is the vector representing average band intensity.We retain only the top K eigen vectors based on the magnitude of the eigenvalues.The most obvious feature of the principal components is that the maximum spectral variability of hyperspectral bands is contained in the first principal component and each subsequent component has variability decreasing exponentially (Holden and LeDrew, 1999).Percentage of variability included in i th principal component is given by, Since K is less than B, one cannot reconstruct hyperspectral image exactly.However, the hyperspectral bands are highly correlated, and hence only a small value of K suffices to reconstruct HS image to obtain the required details.

Close Approximation to Super-resolution using CS
Here, a dictionary based approach is used on the first LR principal component to learn the HR details.Training database of the same class is used to create a set of the HR and LR patches, so that BIPs of PCA better represent the materials of interest.Thus the choice of BIPs depends on the application.Note that dictionaries are not created directly from all bands of HS image, instead first PCA transformed image is used.Suppose LR HS image of size M × M need to be superresolved to a size of qM × qM , where q is the super-resolution factor.Due to similar statistical properties, HS test image patches can be represented as a sparse linear combination of LR dictionary elements using (2).The same sparsity holds good for the corresponding HR unknown image (Divekar and Ersoy, 2009).Hence we can recover the HR image using the HR dictionary.
The proposed algorithm to obtain initial estimation of SR PCA-I image is described below.Here xp gives sparse representation of test patch yp in terms of LR dictionary (DL) elements adaptively .9. Obtain SR patch using ẑp = DH xp.This is SR patch in transform domain.10.Repeat steps (7) to ( 9) for all patches of Y to obtain transform domain SR PCA-I image Ẑ.Note that the spatial dependency still exists within the pixels of this image.However, different PCA components of HS image are uncorrelated.Hence there is no spatial dependency of pixels among different PCA components.The above steps gives us the initial estimate of the SR PCA-I image Z having a size of qM × qM .In our experiments, we considered patch size of b=2 and b=4 for q=2 and q=4, respectively.We used a dictionary of 10000 raw patches.We mention here that one may use a dictionary training algorithm such as K-SVD (Aharon et al., 2006).We observed that, the improvement in result is not very significant compared to the time complexity involved in training the dictionaries.Hence raw patches belonging to the same class were used in our experiment.

Final Solution using Regularization
Since we are estimating the initial SR image using patch based approach, the spatial correlation is not considered whenever there are discontinuities at the patch boundaries.Hence we need to regularize it further to obtain a better solution by restricting the solution space by using AR and sparsity priors.For regularization purpose one needs to have data fitting term and regularization term.Hence we first model the image formation in order to get the data fitting term.
Let Y be the PCA transformed LR HS image of size M × M and and Z be the corresponding HR HS image of size qM × qM , then the model of image formation is represented as: where y and z represent the lexicographically ordered vectors of size M 2 × 1 and q 2 M 2 × 1, respectively with z representing the SR vector.D is the downsampling matrix taking care of aliasing caused as a result of downsampling.For an integer downsampling factor of q, matrix D consists of q 2 non-zero elements along each row at appropriate locations.Here n is the independent and identically distributed (i.i.d.) noise vector with zero mean and variance σ 2 n and has same size as that of y.In this model the LR intensity is the average of the HR intensities over a neighborhood of q 2 pixels corrupted with additive noise.

Estimation of Autoregressive Parameters
We consider learned super-resolved image as the initial estimate and regularize it further to obtain the final solution.We characterize the statistical dependence of pixel values on its neighbors, by using an AR model, where the pixel value at a location is expressed as a linear combination of its neighborhood pixel values and an additive noise (Joshi et al., 2006).We use initially estimated SR PCA-I image using CS approach as an AR model.Note that, PCA components are mutually uncorrelated across bands but there exist spatial dependency among the pixels in each of the PCA components.We estimate the AR parameters from the PCA-I component and use them for all the other bands.We use a homogeneous AR model and derive a set of parameters for entire SR image.Suppose z(s) is the gray level value of the image pixel at location s = (i, j) in an N × N image, where i = 1, 2, ..., N and j = 1, 2, ..., N .The AR model for z(i, j) can be expressed as (Kashyap and Chellappa, 1983) where Ns is the neighborhood of pixel s, r being a neighborhood index with r ∈ Ns, and ρ are unknown parameters; n(.) is an independent and identically distributed noise sequence with zero mean and unit variance; ρ is the variance of the white noise that generates the spatial data for the given AR parameters.Here we use fifth order neighborhood as a compromise between local and global texture representation, that requires to estimate a total of eight AR model parameters Br using the iteration scheme given in (Kashyap and Chellappa, 1983).We are considering same neighborhood size around each pixel.The extracted AR parameters can also be used to regularize other PCA images, as texture properties of other PCA images are partially correlated to the texture properties of the main PCA image.

Regularization with Sparsity Coefficients and AR Parameters
Super-resolution is an ill-posed inverse problem.Prior information can enhance the quality of the solution considerably.Hence we obtain the final solution using regularization framework.Sparse coefficients obtained during the CS framework and AR parameters obtained from the initial super-resolved image are considered as prior informations.The prior knowledge about sparsity coefficients is used in determining weightage of dictionary atoms to represent HR patches.It provides the constraint of sparsity in final solution.The AR parameters plays the role of maintaining the contextual constraint used to regularize the solution.Using a data-fitting term, sparsity and AR prior terms , the final cost function is written as Here β represents the penalty for departure from sparseness in z.
The above cost function is convex.Hence it can be minimized by using a simple optimization technique such as gradient descent.
In order to provide a good initial guess and to speedup the convergence, the result obtained by CS based technique is used as the initial estimate for z.

EXPERIMENTS AND RESULTS
In order to evaluate the performance of the proposed SR technique, experiments are conducted on two different hyperspectral  Figure 5: Plot showing detailed performance of Correlation coefficients for q = 2 of AVIRIS HS bands 1-141 ated LR HS images and compared the results against the original HS image.We limit AVIRIS HSI data in the wavelength range of 0.4µm and 1.79µm to conduct the simulation to reduce the time complexity, as this will not make any significant difference in the results of our proposed algorithm.After removing few bands having low signal to noise ratio (SNR), 141 bands from the original HS image of AVIRIS were super-resolved by a factor of 2 and 4, respectively.The band removal was based on visual inspection of the images.
Effectiveness of the proposed algorithm on the natural scene is presented in fig. 2. Note that, here results are presented only on the PCA-I image instead of HS bands.Visual inspection of fig.2c shows that initial estimate of SR PCA-I using CS based approach gives a quality, comparable to the groundtruth.Since the neighborhood relations are not considered in patches while obtaining the initial estimate, we observe shading of letters written on a ball and also observe little blockiness.Applying regularization to initially estimated SR PCA-I helps to reduce this effect considerably as shown in fig.2d, which is closer to the groundtruth.This shows effectiveness of regularization in proposed algorithm.
In fig. 3  angle mapper (SAM) (Kruse et al., 1993).These are generally used in the multiresolution fusion techniques in order to measure the spatial and spectral fidelity of the fused multispectral images.Since SAM measure is insensitive to variable gain that results from the topographic illumination effects (Kruse et al., 1993) we choose it for measuring the spectral fidelity.

DISCUSSION
The proposed method achieves higher spatial correlations indicated by the CC and ERGAS, which are more closer to reference value.Our approach is also compared with iterative backprojection (IBP) algorithm (Chan et al., 2010).As shown in fig. 3 visual quality achieved by proposed method is better than the IBP method.We can see that pure white patches in LR image is converted to grayish patches in IBP SR image.In IBP the choice of back-projection filter is arbitrary and incorporation of prior information is difficult, which results in higher spatial as well as spectral distortion in SR image compared to proposed approach.
In the proposed approach the use of regularization helps us to achieve better spectral fidelity in terms of lesser value of SAM.
As we increase the variability retained in PCA components, spatial and spectral distortions are reduced considerably.Here we can extend the CS based approach to all significant PCA components, which further improves the performance.

CONCLUSIONS
We have presented a novel approach to recover the high spatial resolution and high spectral resolution HS image using CS based learning and global AR prior model.The advantages of the proposed technique are: 1) no need of supplementary spatial information in registered form, 2) it shows high spatial fidelity and low spectral distortions, and 3) once the low resolution and high resolution dictionaries are created from training dataset, the HS images captured by a low resolution sensor can be super-resolved using our approach.Quantitative comparison of score indices show that our method enhances spatial information without introducing significant spectral distortion.

Figure 1 :
Figure 1: Block diagram of HS image super-resolution.ofBIPs using PCA transform described in section 2.2.The superresolution is applied on the reduced set of PCA transformed images.To reduce computational burden of l1-minimization, we first apply the CS based approach only to first PCA component to obtain initial SR estimate of PCA-I.Regularization based on AR parameters and sparsity priors is performed on this image to obtain the final SR PCA-I image.Remaining transformed images PCA-II-K are interpolated and regularized using AR parameters obtained from the same initial SR image.Here our assumption is that the AR parameters learned from the initial SR PCA-I are valid for all the PCA images, as we are learning the spatial dependence and not their explicit values.Applying inverse PCA on this set of SR PCA images yields super-resolved HS image.
1. Generate mean subtracted LR and HR training HS images of B bands with Lm = [L1, L2, .., LB] and Hm = [H1, H2, .., HB], respectively.Here all bands of LR and HR HS images are arranged lexicographically to convert them in vectors of size M 2 × 1 and q 2 M 2 × 1, respectively.2. Determine basis eigen vectors corresponding to LR and HR training HS images.Retain basis eigen vectors e h and e l corresponding to maximum spectral variability of data (i.e., highest variance), each of size 1 × B. 3. Create HR and LR transformed images BH and BL by projecting Hm and Lm on their corresponding basis eigen vectors generated in step (2).BH = e h × H T m ; BL = e l × L T m ; Convert BH and BL into matrices of size qM × qM and M × M respectively.4. For each b × b patch of BL choose the corresponding patch from BH .Arrange them in lexicographic order to form the vectors V l and V h of size b 2 × 1 and q 2 b 2 × 1, respectively.5. Repeat step (4) to create dictionaries DL and DH of size b 2 × M 2 /b 2 and q 2 b 2 × M 2 /b 2 , respectively.6. Project the mean subtracted LR hyperspectral test image on basis eigen vector e l and obtain the transformed LR test image Y of size M × M .7. Consider a patch of size b × b from the transformed LR test image Y, convert it in lexicographic order to obtain a vector yp.8. Solve CS based l1-minimization optimization problem min||xp||1 such that yp = DL xp + e.
Figure 6: Surface plots showing detailed performance of SAM for q = 2 on AVIRIS data: (a) Bicubic interpolation, and (b) Proposed approach.
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume I-7, 2012 XXII ISPRS Congress, 25 August -01 September 2012, Melbourne, Australia Block diagram of the proposed approach is shown in fig.1.In our approach, we first represent the HS observations from different wavelengths, as a weighted linear combination of small number

Table 1 :
Table1shows quantitative comparison between bicubic interpolation, IBP and proposed approach.Results are listed for different amount of spectral variability (I) retained after transformation.CC is averaged over all bands of HS image to obtain a global measurement of spatial distortion and SAM is averaged over all pixels to yield a global measurement of spectral distortion.As seen from the table 1, our method provides scores that are more closer to reference values compared to bicubic interpolation and IBP.Lower value of ERGAS in the proposed method indicates less global distortion in super-resolved HS image.CC for all the HS bands plotted in fig.5shows that proposed method gives better spatial fidelity.The surface plots in fig.6aand 6b represent spectral fidelity of each specific pixel in super-resolved HS images for q = 2. Proposed method gives a maximum value of 116.22 degree while bicubic interpolation gives 144.89 degree of SAM.Lower values of maximum as well as average SAM indicate that proposed method provides better spectral fidelity.Quantitative comparison of SR results on AVIRIS data for q=2 and q=4