Abstract

Remote sensing images often suffer from stripe noise, which greatly degrades the image quality. Destriping of remote sensing images is to recover a good image from the image containing stripe noise. Since the stripes in remote sensing images have a directional characteristic (horizontal or vertical), the unidirectional total variation has been used to consider the directional information and preserve the edges. The remote sensing image contaminated by heavy stripe noise always has large width stripes and the pixels in the stripes have low correlations with the true pixels. On this occasion, the destriping process can be viewed as inpainting the wide stripe domains. In many works, high-order total variation has been proved to be a powerful tool to inpainting wide domains. Therefore, in this paper, we propose a variational destriping model that combines unidirectional total variation and second-order total variation regularization to employ the directional information and handle the wide stripes. In particular, the split Bregman iteration method is employed to solve the proposed model. Experimental results demonstrate the effectiveness of the proposed method.

1. Introduction

The stripe noise is a common phenomenon that appears in multidetector imaging systems, including the atomic force microscope (AFM) [1], passive millimeter-wave (PMMW) radio [2, 3], moderate resolution imaging spectroradiometer (MODIS) [4], and other systems [5, 6]. The main causes of stripe noise include nonresponse of detectors, relative gain and/or offset variations of detectors, and calibration errors. The above situations severely degrade the quality of the imagery. Therefore, it is critical to recover the pixels contaminated by stripe noise before performing the image interpretation successfully.

A considerable effort has been made to remove stripe noise. The process of removing stripe noise in the striped image is called image destriping. Many methods have been proposed for destriping problems. The first kind of approaches employs digital filtering in the transform domain and suppresses the specific frequency caused by stripes [68]. However, some structural information with the same frequency are also removed along with stripes, in which case this will lead to image blurring and artifacts and can not remove stripes totally. The second kind of methods focuses on the statistical property of digital numbers (DNs) for each detector such as histogram matching and moment matching [4, 9]. To remove stripes totally, we need to adjust the distribution of DNs to a reference one which is based on the similarity assumption of the data. The third kind of methods treats the destriping issue as an inverse problem [1015]. An energy functional is formed by a regularization term and a fidelity term with a regularization parameter to balance the two terms. The minimization of the energy functional is the desired destriping result.

It is known that destriping is an ill-posed problem to reconstruct a high quality image from one degraded image. Recently, the maximum a posteriori (MAP) estimation method becomes popular in image denoising [16], deblurring [17, 18], and superresolution reconstruction [19]. Shen and Zhang [12] first proposed a MAP-based (Shen-MAP) model for the destriping and inpainting problems of remote sensing images. In Shen-MAP model, a Huber-Markov prior which is an alternative between total variation (TV) regularization [2025] and Tikhonov regularization [26] was used and showed a promising performance. However, the model they built did not consider the directional characteristic of stripes in remote sensing images and employed a symmetrical regularization term.

The total variation (TV) model [27] has been attracting more attention in image denoising. It has been proved to be a very efficient denoising approach because of its property of effectively preserving edge information. In [13], the authors proposed a more sophisticated unidirectional TV model and the directional information was used. But when handling heavy irregular stripe noise, it did not work very well. Since TV model will lead to staircase effect, the high-order TV models [2830] were developed to overcome this drawback. Moreover, high-order TV model was also used in image inpainting [31] because of its ability to connect large gaps in inpainting domains. If we consider the heavy stripes in the remote sensing images as dead pixel domains, the destriping problem is a special case of inpainting. For the remote sensing images contaminated with heavy stripe noise, the width of stripes is often large and high-order TV will show its advantage in destriping. In other words, we can employ the high-order TV to handle the heavy stripes with big width. Inspired by this motivation, we propose unidirectional TV and second-order TV model for destriping of remote sensing images which can not only make full use of the stripes’ directional information and preserve edge information but also handle the heavy stripe noise well. The split Bregman iteration method (see details in [32]) is used to solve the proposed model.

The main contributions of this paper include two aspects. First, we propose unidirectional TV and second-order TV (USTV) model for destriping of remote sensing images. Second, we present an efficient algorithm based on the split Bregman method to solve the proposed USTV model. Experiments demonstrate that it can get better visual sense and quantitative results. The rest of this paper is organized as follows. In Section 2, we introduce the image observation model and the Shen-MAP model. In Section 3, we propose our USTV model for image destriping problems and then utilize the split Bregman method to solve the proposed model. In Section 4, the experimental results are presented. We will conclude this paper in Section 5.

2. The MAP Model

2.1. Image Observation Model

A model should be constructed to relate the true image to the contaminated image, by which we can analyze the destriping problem. The model can be linearly described as in [9]. Letting and denote the DN at in the contaminated image and the true image, respectively, the relationship between and can be represented aswhere and are the gain and offset parameters, respectively, and is the sum of linear error and sensor noise. According to (1), for healthy pixel , and . For the destriping problem, the parameters of pixels in a row or a column are often assumed to be the same. In particular, (1) can be rewritten as matrix-vector form:where and are two vectors of the observed image and true image, respectively. is a diagonal matrix whose diagonal elements are the gains of each pixel. is the offset vector, and is the vector containing errors and noise.

2.2. The Shen-MAP Model

To get the true image by the given degraded image , Shen et al. give the below estimation: According to Bayes’ rule, we getBecause is independent of , can be regarded as a constant. Thus (4) can be rewritten asEquation (5) can be rewritten as the following equation due to the monotonic logarithm function :

Two probability density functions (PDF) should be established to get the true image . The first term in (6) is the likelihood density function, which ensures the uniformity of the estimated image to the observed image. According to (2), we have . The errors and noises are assumed to be Gaussian independent with mean zero, by which we can get the PDF of first termwhere is a constant and is a diagonal matrix with the diagonal elements being variances of each pixels. Rewriting (7) aswhere is also a diagonal matrix with where and are the diagonal elements of and , respectively, the second term in (6) is the image prior based on the spatial constraints on the image. A Huber-Markov prior is used for to preserve the reconstructed image edges and other detailed information [33]. The Huber-Markov prior is formulated aswhere is a constant, is a clique within the set of all image cliques , and is the Huber function defined as where is a threshold parameter used to separate the quadratic and linear regions. can be computed by following finite second-order differences to measure the spatial activity for the pixel .

Substituting (6) using (8) and (10), an equivalent formulation of (6) is obtained as follows: whereThe first term of (14) is the fidelity term and the second term is the regularization term, and is the regularization parameter to balance the two terms.

3. The USTV Model

As we can see in Shen-MAP model, the Huber-Markov prior is a symmetrical regularization term which does not consider the directional characteristic of stripes in remote sensing images. Bouali and Ladjal [13] proposed unidirectional TV model to employ the directional information and obtained better results. But when handling heavy irregular stripe noise, it fails to provide satisfying results. When the stripe noise is very heavy, the pixels contaminated by the noise will have less correlation with the true value and the width of stripes in the remote sensing images will be large. In this case, the destriping can be considered as an inpainting problem that fills in the wide stripes of remote sensing images. Recently, high-order TV model was used in image inpainting [31] because of its ability to connect large gaps in inpainting domains. Therefore, to make full use of the directional information and the advantages of high-order TV’s ability to handle heavy stripe noise, we propose unidirectional TV and second-order TV (USTV) model as follows:where , are regularization parameters. and are the first- and second-order discrete differential operators, respectively. We will give the definitions of the notations used in the later subsection. The notations can be found in [34].

3.1. Notations

Without loss of generality, we represent a gray image as matrix. We denote the space as . The discrete gradient operator is a mapping , where . , is given by with

We use and to denote forward difference operators with periodic boundary condition ( is periodically extended). So FFT can be used in our algorithm. We also need the discrete divergence operator that has the adjointness property This property forms the discrete analogue of integration by parts. We define with

We denote backward difference operators with periodic boundary condition as and . Analogously we define the second-order discrete differential operators as a mapping , where . , is given by with

It is easy to know that . The discrete second-order divergence operator is denoted as with the adjointness property We define where . We can check that , , and with

We will present a split Bregman method to solve the above minimization problem in the next subsection.

3.2. The Split Bregman Method for the USTV Model

We first change the unconstrained minimum problem (15) to a constrained problem and then apply the split Bregman method to solve it. Let , and (15) can be written as The corresponding constrained minimization problem isWe can solve (27) by using the split Bregman method.

We use the split Bregman method to solve problem (28). The matrices , , and in (15) are obtained via the same method mentioned in [12].

3.2.1. Subproblem

From (28), we getwhose optimality condition isBy choosing periodic boundary conditions, (31) can be solved by fast Fourier transform (FFT denoted by ) and get where and is the right-hand side of (31). Here “” means pointwise multiplication of matrices. Thus, we have where “/” denotes pointwise division of matrices.

It is worth noting that the periodic boundary condition is a popular boundary assumption used in image processing. Under this assumption, the matrices corresponding to the gradient operator and second-order divergence operator have block circulant with circulant blocks (BCCB) structure, which have a particular spectral decomposition. The matrix-vector multiplications can be calculated by fast Fourier transform more quickly. We can also choose other boundary conditions, for example, zero or reflexive boundary condition, but it may be more difficult or slow to solve subproblem.

3.2.2. Subproblem

According to (28), we haveIt leads to This is a closed form solution, so (35) can be solved simply.

3.2.3. Subproblem

Due to (28), we haveIt can be solved by shrinkage [35]. Denoting , we can get where is assumed.

3.2.4. Subproblem

Because of (28), we haveWe solve (39) adopting the same procedure as in -subproblem and set and we obtain the result as follows: where is assumed. We define the stopping criteria as follows: where is a predetermined positive value. A direct implementation of the USTV model is presented in Algorithm 1.

Input: .
Initialization: , , , , , , .
.
while Terminate condition is not satisfied do
.
The    subproblem: compute according to Eq. (30) by given , , .
The    subproblem: compute according to Eq. (35) by given .
The    subproblem: compute according to Eq. (37) by given .
The    subproblem: compute according to Eq. (39) by given .
Update on  : .
Update on  : .
Update on  : .
end while
Output: Destriping image .
3.3. Convergence Analysis

In this work, we use split Bregman method to solve the optimization problem in (27) and obtained (28). We can reformulate (28) as

The optimization problem is well structured since the four variables , , , and can be separated into two groups and and the convergence of the method can be guaranteed by the existing ADMM theory [36]. Moreover, the variables , , and in the second group are decoupled, and their optimal solutions can be calculated separately. Since all of the subproblems have close form solution, we can solve them easily on their corresponding subproblems in the alternating direction method.

4. Experiments

There are three matrices (gains), (biases), and which should be determined firstly. We calculate them via the same method used in [12]. Given a reference row, the gain and bias of the striped row can be obtained bywhere , are the mean values of the striped row and reference row, respectively, and , are the standard deviations of the striped row and reference row, respectively. The elements in represent the reciprocals of the noise standard deviation in different pixel location and they can be determined bywhere is the standard deviation and and are the std thresholds corresponding to and , respectively. The std was calculated on a neighbor region that does not contain any stripes. Moreover, we do not have to construct the big matrices and and their transpositions during the solution process. Since they are diagonal matrices, the matrix-vector multiplication operations can be obtained by scalar multiplications pixel by pixel.

In our experiments, we test the proposed model on Terra and Aqua MODIS images provided by Shen and Zhang [12] and the test images are shown in Figure 1. We set , , , , , , and in all experiments. From Figures 24, we compare the proposed model with moment matching [9] and Shen-MAP model [12] on Aqua and Terra MODIS images. Visual results of the destriped images clearly show that the proposed USTV model performs destriping more efficiently. The residual stripes of homogeneous regions by the proposed method are fewer than moment matching’s and Shen-MAP model’s. In particular, the stripe noise in Figure 3(a) is much heavier which demonstrates the superiority of the proposed model.

Figure 5 shows the mean cross-track profiles before and after destriping with the proposed USTV model. The horizontal axis represents the line number, and the vertical axis represents the mean DNs of each row. The blue line represents the mean cross-track profile of the striped images, which fluctuates wildly due to the stripes. The red line represents the mean cross-track profile of the destriped images and the fluctuations are suppressed to a large extent. Figure 6 presents the mean column power spectrum that is viewed as a function of normalized frequency before and after destriping with the proposed USTV model. The horizontal axis represents the normalized frequency and the vertical axis represents the averaged power spectrum of all columns. In the case of detector-to-detector stripe noise, the pulses are located at the frequencies of , , , , and cycles. The values of the power spectrum where the pulses are located have been reduced greatly.

We use two quality indexes to evaluate the performance of the proposed USTV model. The first index is inverse coefficient of variation (ICV) [37] defined aswhere and are the mean and standard deviation of DNs, respectively. In particular, the ICV index is computed on homogeneous regions within a window of pixels. The second index is noise reduction (NR) [4, 8] defined as where is the power of frequency components produced by stripes in the striped image and for the destriped image. and can be computed bywhere Dis is the distance from the origin in the Fourier space and is the stripe noise region of the averaged power spectrum down the columns. Three homogeneous regions were selected for the ICV evaluation. The ICV and NR results are presented in Table 1. We compare the proposed USTV model with moment matching [9] and the Shen-MAP model [12]. From Table 1, we can see that the proposed USTV model obtains the highest ICV value and NR value.

5. Conclusion

In this paper, we proposed unidirectional TV and second-order TV model for destriping of the remote sensing images. The split Bregman iteration method was used to solve the proposed model. To the best of our knowledge, this is the first work to combine unidirectional TV and second-order TV for the destriping of remote sensing images. In addition, we employed two image quality indexes ICV and NR to evaluate the performance of different methods. Experiments demonstrated that the proposed model obtained better destriping effects than the moment matching method and the Shen-MAP model in terms of both visual sense and quantitative assessments. However, the degradation process in (2) that we use is a linear model and the results rely on parameters and by moment match. To some extent, it is still a match-based destriping method. A possible future work is to consider the degradation process as an additive process and develop a new unidirectional TV and second-order TV model.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work is supported by 973 Program (2013CB329404) and NSFC (61370147, 61402082).