Analysis and Denoising of Hyperspectral Remote Sensing Image in the Curvelet Domain

Abstract

A new denoising algorithm is proposed according to the characteristics of hyperspectral remote sensing image (HRSI) in the curvelet domain. Firstly, each band of HRSI is transformed into the curvelet domain, and the sets of subband images are obtained from different wavelength of HRSI. And then the detail subband images in the same scale and same direction from different wavelengths of HRSI are stacked to obtain new 3-D datacubes of the curvelet domain. Again, the characteristics analysis of these 3-D datacubes is performed. The analysis result shows that each new 3-D datacube has the strong spectral correlation. At last, due to the strong spectral correlation of new 3-D datacubes, the multiple linear regression is introduced to deal with these new 3-D datacubes in the curvelet domain. The simulated and the real data experiments are performed. The simulated data experimental results show that the proposed algorithm is superior to the compared algorithms in the references in terms of SNR. Furthermore, MSE and MSSIM in each band are utilized to show that the proposed algorithm is superior. The real data experimental results show that the proposed algorithm effectively removes the common spotty noise and the strip noise and simultaneously maintains more fine features during the denoising process.

1. Introduction

Hyperspectral remote sensing image (HRSI) can be viewed as three-dimensional data consisting of one-dimensional spectral information and two-dimensional spatial information. With the fast development of hyperspectral remote sensing technology, HRSI can describe the characteristics of Earth objects more comprehensively and explicitly therefore, it is widely applied in many fields including agriculture, forestry, geological surveys, environmental monitoring, and military recon. Although over the last decades the development of imaging spectrometers is rapid, HRSI is still affected by many complex factors during the processing of acquisition and transmission, which will produce a mass of noises. The data that are contaminated with noise can cause a failure to extract valuable information and hamper further interpretation. In presence of noise in the image, extraction of all the useful information becomes difficult and noise can lead to artefacts and loss of spatial resolution [1]. The noise also affects the target detection, classification, and segmentation, so it is necessary to study the characteristic of HRSI for denoising [2].

Though the spatial correlation of HRSI is weaker than the nature image, the two-dimensional spatial information of HRSI is similar to the nature image [3]. When some traditional linear denoising methods are used to denoise the HRSI, they cannot get the satisfactory results for application such as Wiener filter. Compared with the nature image, the one-dimensional spectral information of HRSI is its particular characteristic. Minimum noise fraction (MNF), Savitzky-Golay filter, and wavelet denoising are the most popular among the existing denoising methods of hyperspectral imagery [4–7]. These methods are able to smooth the spectral while they have a negative impact on the sharp signal features.

In recent years, many denoising methods for HRSIs are constantly introduced with the development of hyperspectral remote sensing technology. Most of denoising methods combine correlation of spatial and spectral domain. Currently, wavelet denoising methods are widely used. The seminal work on signal denoising via wavelet thresholding or shrinkage proposed by Donoho and Johnstone [6, 7] shows that various wavelet thresholding schemes for denoising have near-optimal properties in the minimax sense. The observed noisy data in a local neighborhood is utilized to perform an approximate maximum a posteriori estimation of the variance for each coefficient [8]. Then, the predicted image coefficients are obtained by an approximate minimum mean squared error estimation procedure. Şendur and Selesnick developed a bivariate shrinkage function for image denoising [9, 10]. Their results show that the estimated wavelet coefficients depend on the parent coefficients. The smaller parent coefficients mean the greater shrinkage. In [11], Buades et al. proposed a new filter called the nonlocal means (NLM) filter, which takes account of the two most important attributes of a denoising algorithm: detail preservation and noise removal. Unfortunately, as most of filtering methods, NLM filter is also based on the additive Gaussian signal-independent (SI) noise assumption, so it cannot be directly applied for HRSI with signal-dependent (SD) noise [12]. In order to fully utilize the spatial information and the spectral information, Atkinson et al. proposed a denoising method that uses discrete Fourier transforms, 2-D discrete wavelet transforms, and soft thresholding of wavelet coefficients to denoise HRSI [13]. In [14], Othman and Qian developed a noise reduction algorithm—hybrid spatial-spectral noise-reduction algorithm (HSSNR) for hyperspectral datacube. The algorithm resorts to the spectral derivative domain, where the noise level is elevated, and benefits from the dissimilarity of the signal regularity in the spatial and spectral domains. Chen and Qian [15] combined wavelet denoising with dimensionality reduction for HRSI by using bivariate wavelet thresholding, wavelet packets, and principal component analysis (PCA). In [16], Chen and Qian developed a new denoising method using PCA and wavelet shrinkage to avoid removing the fine features of the datacubes during denoising process. PCA is used to decorrelate the fine features of the datacube from noise for detail preservation. In this paper the method in [16] is named as PCABS.

Since wavelet has good time-frequency-localization property and multiresolution analysis property, it is widely and successfully applied in several fields [17]. However, wavelet is not perfect. Wavelet is mainly applied to representation of isotropic singularity object, while for anisotropic singularity, such as boundary and linear features of an image, it is not a good representation tool. In other words wavelet is good representation for point singularity, but for the curves it is relatively weak [18–20]. The complex wavelet [18], curvelet transform [19], and contourlet transform [20] are proposed to overcome the drawbacks of the wavelet transform. This paper focuses on the curvelet transform applied for HRSI.

Compared with normal three-dimensional data cube of fixed variance of additive noise, the noise level of HRSI may vary dramatically from band to band. The noise standard deviation in each band of HRSI is not constant; in particular, there exist some bands at which the atmosphere absorbs so much light that the signal received from the surface is unreliable [21]. And SD noise may no longer be neglected with respect to the SI noise in the new-generation hyperspectral sensors [22, 23]. In this paper, the proposed denoising method is available for the both SD noise and SI noise. The representation of HRSI in another domain is proposed for detail preservation. By analyzing spectral correlation of HRSI in the new representation domain, we know the fact that the spectral correlation of HRSI in the new representation domain is strong. So the multiple linear regression (MLR) method is proposed to predict the representation of pure signal in the curvelet domain, in order to remove the noise of HRSI.

In the rest of this paper, the mathematical tools curvelet transform and MLR model are introduced in Section 2. Section 3 gives the analysis of HRSI in the curvelet domain. In Section 4 the denoising process is proposed in the curvelet domain. The simulated data experiment and the real data experiment are performed, and the experimental results are presented in Section 5. Finally, Section 6 draws the conclusion.

2. Mathematical Tools

In the section, the mathematical tools curvelet transform and MLR model are introduced. The curvelet transform, pioneered by Candès and Donoho, is shown to be optimal in a certain sense for functions in the domain with curved singularities [19]. Due to the strong spectral correlation of HRSI, MLR model has been widely applied to the HRSI [24, 25].

2.1. The Curvelet Transform

The curvelet transform is a new multiresolution analysis framework and widely applied in various image processing problems. The curvelet decomposition can be equivalently stated in the following four steps: (1) subband decomposition, (2) smooth partitioning, (3) renormalization, and (4) ridgelet transform. In short, the curvelet obtained by bandpass filtering of multiscale ridgelets with passband is rigidly linked to the scale of spatial localization. The discrete curvelet transform [26] is given as follows.

The digital image is used as an example to introduce the discrete curvelet transform. 𝐼 stands for the original image, whose size is 𝑁×𝑁. And then 𝐼 is decomposed by binary wavelet transform with 𝐽 scales, and we have 𝐼=𝐶_𝐽+^𝐽_𝑗=1𝐷_𝑗, (1) where 𝐶_𝐽 is low-frequency component in the coarsest scale 𝐽, while {𝐷_𝑗}_{𝑗=1,2,…𝐽} is high-frequency component in different scale. Here, 𝑗=1 is the finest scale. We now employ a sketch of the discrete curvelet transform algorithm.(1)Perform the binary wavelet transform with 𝐽 scales on the original image 𝐼, and obtain subband sequence 𝐶_𝐽 and {𝐷_𝑗}_{𝑗=1,2,…𝐽}.(2)Let the size of initial block (the subband in the finest scale) be 𝐵_min, and set 𝐵₁=𝐵_min.(3)For 𝑗=1,2,…,𝐽(a)partition the subband 𝐷_𝑗 with a block size 𝐵_𝑗;(b)apply the digital ridgelet transform to each block;(c)If 𝑗 module 2=1𝐵_𝑗+1=2𝐵_𝑗, and else 𝐵_𝑗+1=𝐵_𝑗.

Since each step of the previous decomposition process is invertible, the inverse curvelet transform is an invertible process.

The curvelet transform overcomes the major drawback that wavelets cannot really represent two-dimensional objects with edges sparsely and captures more directional information besides the horizontal, vertical, and diagonal directions. The system approximately obeys the scale relationship width≈length². Therefore, it is optimal in a certain sense for functions in the domain with curved singularities.

It is known that the curvelet transform competes surprisingly well with the ideal adaptive rate. The approximation error is obtained as 𝑓−𝑓^𝐶_𝑚²₂≤𝐶⋅𝑚⁻²log𝑚³,𝑚⟶∞, (2) where 𝑓^𝐶_𝑚 is the 𝑚 biggest terms in the curvelet frame expansion to approximate 𝑓.

Curvelet is optimal in the sense that no other representation can yield a smaller asymptotic error with the same number of terms. Because of its surprising properties for image processing, a fast and accurate discrete curvelet transform operating on digital data is necessary. Candès et al. [27, 28] presented two 2-D discrete curvelet transforms for the second generation curvelet, which is curvelet via transforms and curvelet via wrapping of specially selected Fourier samples. Compared with the first generation curvelet, they are conceptually simpler, faster, and far less redundant. The curvelet transform used in this work is based on unequally-spaced fast Fourier transform.

2.2. Multiple Linear Regression (MLR)

Due to the strong spectral correlation of new 3-D datacubes in the curvelet domain, the MLR model is introduced to predict the representation of pure HRSI in the curvelet domain. It is assumed that the HRSI has 𝐵 spectral bands and each band of HRSI has 𝑀×𝑁 pixels. Let 𝑋 denote a 𝑃×𝐵 matrix of the 𝐵 spectral observed vectors of size 𝑃(𝑃=𝑀×𝑁). The 𝑃×1 vector 𝑋_𝑘(𝑘=1,2,…,𝐵) is the 𝑘th column vector of the matrix 𝑋, namely, lexicographically arranging 𝑘th band image into a column vector 𝑋_𝑘. In this paper, 𝑋_𝑘 is the vector predicted for the signal 𝑋_𝑘 of band 𝑘. The 𝐿 adjacent bands (not including itself) are utilized to perform MLR, where 𝐿 is even; that is, 𝑋_𝑘=𝑋_𝜆𝑘𝛽_𝑘(𝑘=1,2,…,𝐵), (3) where the 𝑃×𝐿 matrix 𝑋_𝜆𝑘 consists of the adjacent column vectors of 𝑋_𝑘 (not including the 𝑘th column vector). 𝑋_𝜆𝑘 has the form 𝑋_𝜆𝑘=𝑋₁,…,𝑋_𝑘−1,𝑋_𝑘+1,…,𝑋_𝐿𝐿1≤𝑘≤2𝑋_{𝑘−𝐿/2},…,X_𝑘−1,𝑋_𝑘+1,…,𝑋_𝑘+𝐿/2𝐿2𝐿<𝑘<𝐵−2𝑋_{𝐵−𝐿+1},…,𝑋_𝑘−1,𝑋_𝑘+1,…,𝑋_𝐵𝐿𝐵−2≤𝑘≤𝐵. (4)

𝛽_𝑘 is the regression vector of size 𝐿×1. For 𝑘=1,2,…,𝐵, the least squares estimator of the regression vector 𝛽_𝑘 is given by 𝛽_𝑘=𝑋^𝑇_𝜆𝑘𝑋_𝜆𝑘⁻¹𝑋^𝑇_𝜆𝑘𝑋_𝑘. (5)

3. Analysis of HRSI in the Curvelet Domain

HRSI is a datacube, having two spatial dimensions and a third spectral dimension. Fixing the wavelength band yields a 2-D image of the scene at a particular wavelength. So it may also be visualized as a stack of 2-D band images, each corresponding to a certain wavelength [3]. Since the bands are so closely spaced in wavelength, images in adjacent bands are highly correlated. In the rest of the section, HRSI is transformed into the curvelet domain band by band and the spectral correlation of the HRSI in the curvelet domain will be mainly discussed to denoise by performing the MLR.

Let 𝐶(⋅) and 𝐶⁻¹(⋅) denote the 2-D curvelet transform and 2-D inverse curvelet transform. Let {𝑋_𝑖}^𝐵_𝑖=1 be HRSI, where 𝑋_𝑖 is the 𝑖th band image and 𝐵 is the number of bands. 𝐶_𝑖=𝐶(𝑋_𝑖) stands for the coefficients of band 𝑋_𝑖 in the curvelet domain. A 𝐽-level curvelet transform is performed on each 𝑋_𝑖. Setting aside the coarsest coefficients 𝐶^𝐿_𝑖 of 𝐶_𝑖, the same scale and same directional coefficients from each high frequency coefficient 𝐶^𝐻_𝑖(𝐶^𝐻_𝑖=𝐶_𝑖/𝐶^𝐿_𝑖) are stacked as a new 3-D datacube. The previous process is shown in Figure 1. In this way, several new 3-D datacubes {𝐹^𝑡}^𝑆_𝑡=1 are obtained where 𝑆 is the number of data. However, transforming {𝐹^𝑡}^𝑆_𝑡=1 into 𝐶_𝑖(𝑖=1,2,…𝐵) is an inverse transform of the previous process. The spatial size of each datacube 𝐹^𝑡 varies with the scale and the direction, but each datacube has the same 𝐵 spectral bands. The new three-dimensional datacubes are new representation of HRSI in the curvelet domain. The spectral correlation of 𝐹^𝑡 will be discussed. In this paper, each band 𝑋_𝑖 is decomposed into 5 levels. The first scale is coarsest scale with only one direction. From the second scale to the fourth scale, there are eight directions. At the finest scale, one direction is obtained. Setting aside the first scale, the total number of the other directions from different scale is 25; that is, 𝑆=25.

Figure 1 

The process of obtaining new 3-D datacubes {𝐹𝑡}𝑆𝑡=1.

Let 𝐹^𝑡_𝑙(𝑚,𝑛) denote the value of the position (𝑚,𝑛,𝑙) of the datacube. Let 𝐻 denote the set of spatial coordinates of 𝐹^𝑡, and let |𝐻| denote the element number of set 𝐻.

The spectral correlation factor of 𝐹^𝑡 between band 𝑘 and band 𝑙 is defined as 𝜌_𝑡()=∑𝑘,𝑙_{(𝑚,𝑛)∈𝐻}𝐹^𝑡_𝑘(𝑚,𝑛)−𝑢_𝑘𝐹^𝑡_𝑙(𝑚,𝑛)−𝑢_𝑙∑_{(𝑚,𝑛)∈𝐻}𝐹^𝑡_𝑘(𝑚,𝑛)−𝑢_𝑘²∑_{(𝑚,𝑛)∈𝐻}𝐹^𝑡_𝑙(𝑚,𝑛)−𝑢_𝑙², (6) where 𝑢_𝑘=(1/|𝐻|)∑_{(𝑚,𝑛)∈𝐻}𝐹^𝑡_𝑘(𝑚,𝑛) and 𝑢_𝑙=(1/|𝐻|)∑_{(𝑚,𝑛)∈𝐻}𝐹^𝑡_𝑙(𝑚,𝑛). In this paper, let 𝑘=𝑙+1; thus 𝜌_𝑡(𝑘,𝑘+1) stands for the correlation factor of neighboring bands for datacube 𝐹^𝑡. In order to analyze the spectral correlation of HRSI in the curvelet domain, AVIRIS hyperspectral data Jasper Ridge whose size is 256 × 256 × 224 (width × height × band) is used as an experimental example. Figure 2 shows the spectral correlation of each datacube 𝐹^𝑡 in the curvelet domain and the spectral correlation of the datacube Jasper Ridge.

Figure 2 

The spectral correlation factor of {𝐹𝑡}𝑆𝑡=1 (a) eight directions of the second scale, (b) eight directions of the third scale, (c) eight directions of the fourth scale, (d) the finest scale, and (e) the spectral correlation factor of the Jasper Ridge.

From Figure 2, it is clear that the absolute value of linear correlation factor between neighboring bands is close to 1 from the second scale to the fourth scale except for the bands 1–5, 105–115, and 150–170, which is due to water absorption and low signal-to-noise ratio (SNR) for the HRSI. These bands are often viewed as junk band. From the second scale to the finest scale, the correlation factor of junk bands in the curvelet domain becomes much lower, but the other bands are still keeping close to 1. The result in the curvelet domain is consistent with the original HRSI. It also indicates that the curvelet coefficients of hyperspectral data between neighboring bands have significant linear correlation from the same scale and the same direction. The curvelet transform keeps the strong spectral correlation of HRSI, and specially at some directions from the second scale to the forth scale the spectral correlation in the curvelet domain becomes stronger. The reason for this phenomenon is the strong spectral correlation of HRSI, which means that the pixels in the same spatial location of each band image are similar. Images at different wavelengths capture the same scene, which has the same physical structure; thus the profiles of different wavelength images are similar [29]. After the curvelet transform, the coefficients of datacube 𝐹^𝑡 at close wavelengths maintain this similarity as well and have significant linear correlation. Since the most of subtle noise of HRSI is maintained in the finest scale, the correlation factor of the finest scale is lower than the other scales. But in the finest scale, the most of the correlation factors are still greater than 0.85. Due to the strong spectral correlation of HRSI in the curvelet domain, a new denoising algorithm based on MLR is proposed in the next section.

4. Proposed Denoising Algorithm

In this section, we summarize our denoising algorithm. According to the correlation factor of HRSI in the curvelet domain in Section 3, the MLR model is performed on the datacubes {𝐹^𝑡}^𝑆_𝑡=1 for denoising.

The denoising process is as follows (Figure 3). (1)Input the noisy datacube 𝑋.(2)A 𝐽-level curvelet transform (CT) is performed on 𝑋_𝑖, and 𝐶_𝑖=𝐶(𝑋_𝑖)(𝑖=1,2,…𝐵) is obtained.(3)Set aside the coarsest scale coefficients 𝐶^𝐿_𝑖 of 𝐶_𝑖, and then the same scale and same directional coefficients from each high frequency coefficients 𝐶^𝐻_𝑖(𝐶^𝐻_𝑖=𝐶_𝑖/𝐶^𝐿_𝑖) are stacked as a 3-D datacube. The 3-D datacubes {𝐹^𝑡}^𝑆_𝑡=1 are obtained.(4)MLR is performed on 𝐹^𝑡, 𝑡=1,2,…𝑆, and 𝐹^𝑡(𝑡=1,2,…𝑆) is obtained. (5)For 𝑖=1,2,…𝐵, combine the new high frequency coefficients 𝐶^𝐻_𝑖 extracted from each 𝐹^𝑡_𝑖(𝑡=1,2,…𝑆) with the corresponding coarsest scale coefficients 𝐶^𝐿_𝑖 to obtain the prediction coefficients 𝐶_𝑖 of the 𝐶_𝑖.(6)The inverse curvelet transform (ICT) is performed to obtain the denoised image 𝑋_𝑖; that is, 𝑋_𝑖=𝐶⁻¹(𝐶_𝑖).(7)Output the denoised datacube 𝑋.

Figure 3 

Block diagram of the proposed method in this paper.

5. Experiment

Signal-to-noise ratio (SNR) is a key parameter on measuring the HRSI quality. So in this paper, we utilize SNR to evaluate the proposed algorithm. Here the SNR is defined as 𝑃SNR=^𝑋𝑃_𝑁, (7) where 𝑃_𝑋 is the power of the pure signals 𝑃𝑋(𝑖,𝑗,𝑘), and 𝑃_𝑁 is the noise power in the noisy signals 𝑋(𝑖,𝑗,𝑘), while (𝑖,𝑗,𝑘) stands for the position of the pixel in the HRSI; that is, ∑SNR=^{𝑀,𝑁,𝐵𝑖=1,𝑗=1,𝑘=1}𝑃𝑋𝑖,𝑗,𝑘²∑^{𝑀,𝑁,𝐵}_{𝑖=1,𝑗=1,𝑘=1}𝑃𝑋𝑖,𝑗,𝑘−𝑋𝑖,𝑗,𝑘². (8)

5.1. Simulated Data Experiment

The simulated experiment of the noise reduction is carried out on AVIRIS images, Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field provided by JPL, NASA. The size of datacube we extracted from the Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field for testing is 256 × 256 × 224 (width × height × band). Figure 4 shows the band no.40 of these HRSIs. For the AVIRIS images a 28 m × 28 m ground sample distance (GSD) datacube is derived by spatially averaging the 4 m × 4 m GSD datacube elevating the nominal SNR to 7000 : 1. Having such high SNR, this datacube is viewed as a pure datacube [14], which is used as a reference to measure the SNR before and after denoising. The image is corrupted by Gaussian white noise. It is different from the simple stationary additive noise model that is simulated by adding noise with a fixed standard deviation to the datacube; the noise variance is proportional to the average signal amplitude of each band, but the noise in each band is still additive noise. The SNR of the simulated noisy data is 600 : 1, which is chosen by comprehensive requirement of the users and the machine design parameters.

(a) Cuprite

(b) Jasper Ridge

(c) Low Altitude

(d) Lunar Lake

(e) Moffett Field

Figure 4 

Band no.80 of AVIRIS images.

In order to indicate that the nature image denoising methods cannot be immediately used for HRSI noise reduction, the 2-D complex wavelet with bivariate shrinkage (CWBS) [10] and the curvelet denoising (CD) [26] are used for comparison. From Table 1, it is shown that the denoised results in terms of SNR are obtained by the two nature image denoising methods; even the lower SNR is obtained. In order to illustrate the superiority of the proposed algorithm in our paper, the HSSNR [14] and the PCABS [16] are also performed to compare with the proposed method in this paper. The parameters chosen for the PCABS are consistent with [16]. In this paper, through utilizing [0.2×𝐵] bands to get the best result through the lots of simulated experiments, the proposed algorithm utilizes [0.2×𝐵] bands to perform the MLR, where [𝑥] means getting the nearest even integer of 𝑥. Table 1 shows the SNR of the hyperspectral datacubes Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field after denoising by the HSSNR, the PCABS, and the proposed algorithm. The results indicate that the proposed method is best method for denoising HRSI in terms of SNR.

Table 1 

SNR of AVIRIS data Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field.

In order to deeply analyze the proposed algorithm, the mean square of errors (MSE) in each band and the mean structural similarity (MSSIM) [30] between before and after denoising are utilized in each band. The lower MSE means that the denoised image is more similar to the original pure image, while MSSIM is close to 1 means that the denoised image is more similar to the original pure image. The MSE of band 𝑘 denoised image is defined as MSE_𝑘=1𝑀𝑁^𝑀,𝑁_{𝑖=1,𝑗=1}𝑋^𝑘_𝑖,𝑗−𝑃𝑋^𝑘_𝑖,𝑗², (9) and the MSSIM of band 𝑘 denoised image is defined as MSSIM_𝑘𝑋^𝑘,𝑃𝑋^𝑘=1𝑄^𝑄_𝑞=1SSIM𝑋^𝑘_𝑞,𝑃𝑋^𝑘_𝑞, (10) where 𝑋^𝑘_𝑡 and 𝑃𝑋^𝑘_𝑡 are the image contents at the 𝑞th local window, and 𝑄 is the number of local windows of the image. According to [30], the SSIM is defined as SSIM(𝐱,𝐲)=2𝜇_𝐱𝜇_𝐲+𝐶₁2𝜎_𝐱𝐲+𝐶₂𝜇²_𝐱+𝜇²_𝐲+𝐶₁𝜎²_𝐱+𝜎²_𝐲+𝐶₂, (11) where 𝐶₁=(𝐾₁𝑉)² and 𝐶₂=(𝐾₂𝑉)² (𝐾₁≪1 and 𝐾₂≪1 are two small constants, 𝐾₁=0.01 and 𝐾₂=0.03) and 𝑉 is the dynamic range of the pixel values. In this paper, we use 11 × 11 circular-symmetric Gaussian weighting function 𝐰={𝑤_𝑖∣𝑖=1,2,…,𝑁}, with standard deviation of 1.5 samples, normalized to unit sum (∑^𝑁_𝑖=1𝑤_𝑖=1). The estimates of local statistics 𝜇_𝐱, 𝜎_𝐱 and 𝜎_𝐱𝐲 are then modified accordingly as 𝜇_𝐱=^𝑁_𝑖=1𝑤_𝑖𝑥_𝑖,𝜎_𝐱=^𝑁_𝑖=1𝑤_𝑖(𝑥_𝑖−𝜇_𝐱)^21/2,𝜎_𝐱𝐲=^𝑁_𝑖=1𝑤_𝑖𝑥_𝑖−𝜇_𝐱𝑦_𝑖−𝜇_𝐲, (12)(13)(14) where 𝜇_𝐲 and 𝜎_𝐲 are consistent with the form 𝜇_𝐱 and 𝜎_𝐱.

From Table 1, it is shown that the PCABS method and the proposed method obtain higher SNR. As a result, Figure 5 only shows MSE and MSSIM in each band of the AVIRIS datacubes denoised by the proposed algorithm and the PCABS. The MSE of each band obtained by the proposed algorithm is lower than the PCABS in nearly all bands, and the MSSIM of each band obtained by the proposed algorithm is closer to 1 than the PCABS, which indicates that the proposed algorithm is superior to PCABS.

(a) Cuprite  

(b) Jasper Ridge

(c) Lunar Lake

(d) Low Altitude

(e) Moffett Field

Figure 5 

The MSE and the MSSIM in each band obtained by PCABS and the proposed method.

The computational complexity, that is, the number of floating-point operations (flops), of the proposed method can be analyzed as follows. The complexity of the curvelet transform is in order of O(𝑃𝐵⋅log(√𝑃)) flops, where 𝑃 is the number of pixels in the spatial domain and 𝐵 is the number of bands. The complexity of the MLR is approximately the number of floating-point operations 4𝑃𝐿𝐵+2𝑃𝐿²𝐵+𝐿³𝐵 flops where 𝑃 and 𝐵 are defined as before and 𝐿 is the number of bands used to perform MLR. Therefore, the computational complexity of the proposed method is mainly contributed by the MLR. The computational complexity of the proposed method is greater than that of the HSSNR and PCABS. The computational times of the HSSNR, the PCABS, and the proposed method are given in Table 2. The previous denoising approaches are implemented using the Matlab programming language and run on a PC with a Pentium 2.70 GHz Dual-Core CPU and a 1024 MB RAM.

Table 2 

The computational time of AVIRIS data Cuprite, Jasper Ridge, Low Altitude, Lunar Lake, and Moffett Field (units: s).

5.2. Real Data Experiment

In this paper the OMIS (operational modular imaging spectrometer) data that is developed by the Shanghai Institute of Technical Physics of the Chinese Academy of Sciences is used for real data experiment to verify the correctness and performance of algorithm. It has 128 spectral bands ranging from visible to thermal infrared wavelength. The size of datacube we extracted from OMIS data for testing is 256 × 256 × 128 (width × height × band). We perform denoising for the original OMIS data. Figure 6 shows band nos. 20, 40, 60, 80, and 100 of the original image and the denoised image obtained by PCABS and the proposed method, as well as the difference between the original image and the denoised image obtained by PCABS and the proposed method. In Figure 6, the denoised images (b) and (c) show that the proposed method is feasible; especially the band no. 80 can be seen obviously. The image (d) is the difference between the original image and the denoised image obtained by PCABS, and the image (e) is the difference between the original image and the denoised image obtained by the proposed algorithm. The difference shows that the two algorithms can effectively remove the common spotty noise and the strip noise from the noisy OMIS. But the features such as the edge of the image can be observed in the difference. The reason for this phenomenon is that the fine features are removed as noise during the denoising process. So the less fine features in the difference mean the better denoised results. From Figures 6(d) and 6(e), the image (d) contains more fine features than the image (e), which is easy to distinguish by human eyes. The result shows that the PCABS method removes more fine features during the denoising process. Therefore the proposed algorithm is superior to the PCABS in terms of removing noise and simultaneously maintaining fine features during the denoising process.

Figure 6 

Band nos. 20, 40, 60, 80, and 100 of (a) the original image and (b) the denoised image obtained by the PCABS, (c) the proposed algorithm, and (d) the difference between the original image and the denoised image obtained by PCABS, as well as (e) the difference between the original image and the denoised image obtained by the proposed algorithm.

6. Conclusions

In this paper, the spectral correlation of HRSI in the curvelet domain is discussed. By the analysis, in the curvelet domain, the strong spectral correlation of the hyperspectral remote sensing image is kept; even in some directions and scales it becomes stronger. So a new denoising algorithm is proposed; the MLR is performed in the curvelet domain to denoise the HRSI.

Simulated experimental result shows that the proposed method improves the quality of HRSI significantly in terms of SNR, MSE of each band, and MSSIM of each band. It is also seen that the denoised results obtained by the two algorithms are not content in the bands 1–5. This is a problem that will be studied in the future. For the real OMIS data, the results show that the proposed method is valid. The proposed method obtains better results in terms of detail preservation and noise removal during the denoising process.

Acknowledgment

This work was supported by the National Natural Science Foundation of China under Project 61101183 and Project 41201363.