A Unidirectional Total Variation and Second-Order Total Variation Model for Destriping of Remote Sensing Images

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.


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 [6][7][8].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 [10][11][12][13][14][15].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 [20][21][22][23][24][25] 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 [28][29][30] 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.

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 as where  , and  , are the gain and offset parameters, respectively, and  , is the sum of linear error and sensor noise.According to (1), for healthy pixel (, ),  , = 1 and  , = 0.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 f and u are two vectors of the observed image and true image, respectively.A is a diagonal matrix whose diagonal elements are the gains of each pixel.B is the offset vector, and  is the vector containing errors and noise.

The Shen-MAP Model.
To get the true image u by the given degraded image f, Shen et al. give the below estimation: According to Bayes' rule, we get Because (u | f) is independent of f, (f) can be regarded as a constant.Thus (4) can be rewritten as Equation ( 5) can be rewritten as the following equation due to the monotonic logarithm function () = ln(): 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 (f | u) = ().The errors and noises are assumed to be Gaussian independent with mean zero, by which we can get the PDF of first term where  1 ∈ R is a constant and K is a diagonal matrix with the diagonal elements being variances of each pixels.Rewriting (7) as where Q is also a diagonal matrix with where   and   are the diagonal elements of Q and K, 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 (u) to preserve the reconstructed image edges and other detailed information [33].The Huber-Markov prior is formulated as where  2 ∈ R is a constant,  is a clique within the set of all image cliques , and (⋅) is the Huber function defined as where  ∈ R 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: where The 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.

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 ∇ 2  = ( +  ,  +  ,  +  ,  +  ) 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 R × as H.The discrete gradient operator is a mapping ∇ : 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 div : V → H 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 ( ) , with It is easy to know that  +  =  +  ( +  ) =  +  ( +  ) =  +  .The discrete second-order divergence operator is denoted as div 2 : W → H with the adjointness property We define where  = ( We will present a split Bregman method to solve the above minimization problem in the next subsection.

The Split Bregman
The corresponding constrained minimization problem is min ,u We can solve ( 27) by using the split Bregman method.
We use the split Bregman method to solve problem (28).The matrices A, B, and Q in (15) are obtained via the same method mentioned in [12].

𝑔 Subproblem. From (28), we get
whose optimality condition is By choosing periodic boundary conditions, (31) can be solved by fast Fourier transform (FFT denoted by F(⋅)) and get where and  right 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 matrixvector 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.

u Subproblem.
According to (28), we have It leads to This is a closed form solution, so (35) can be solved simply.
and we obtain the result as follows: ( where 0/0 = 0 is assumed.We define the stopping criteria as follows: where  ∈ R is a predetermined positive value.A direct implementation of the USTV model is presented in Algorithm 1.

− (
The optimization problem is well structured since the four variables , u, V, and  can be separated into two groups  and (u, V, )  and the convergence of the method can be guaranteed by the existing ADMM theory [36].Moreover, the variables u, V, 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.

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 by where   ,   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 by where std is the standard deviation and max and min are the std thresholds corresponding to 1 and 0, 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 matrixvector 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  = 1.0 × 10 −6 ,  = 5 × 10 4 ,  = 1,  = 0.8,  1 = 10 max(, ),  2 = 20, and  3 = 10 in all experiments.From Figures 2-4, 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 0.1, 0.2, 0.3, 0.4, and 0.5 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 as  where   and  sd are the mean and standard deviation of DNs, respectively.In particular, the ICV index is computed on homogeneous regions within a window of 10 × 10 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  de for the destriped image.  and  de can be computed where Dis is the distance from the origin in the Fourier space and BW  is the stripe noise region of the averaged power spectrum down the columns.Three 10 × 10 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.

Conclusion
In this paper, we proposed unidirectional TV and secondorder 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.
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 f − B = C, and (15) can be written as min u  2 ‖Q (Au − C)‖ 2 +       ∇  u     1 +       ∇ 2 u     1 .

Table 1 :
ICVs and NRs of the comparison experiments.