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

Authors: Min Wang; Ting-Zhu Huang; Xi-Le Zhao; Liang-Jian Deng; Gang Liu
Journal: Mathematical Problems in Engineering
 (2017)

---

## 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.


---

## Body

## 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 [6–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–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–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–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 Section2, 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 fi,j and ui,j denote the DN at (i,j) in the contaminated image and the true image, respectively, the relationship between fi,j and ui,j can be represented as(1)fi,j=Ai,jui,j+Bi,j+εi,j,where Ai,j and Bi,j are the gain and offset parameters, respectively, and εi,j is the sum of linear error and sensor noise. According to (1), for healthy pixel (i,j), Ai,j=1 and Bi,j=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:(2)f=Au+B+ε,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.


### 2.2. The Shen-MAP Model

To get the true imageu by the given degraded image f, Shen et al. give the below estimation: (3)u^=arg⁡maxu⁡pu∣f.According to Bayes’ rule, we get(4)u^=arg⁡maxu⁡pf∣upupf.Because pu∣f is independent of f, pf can be regarded as a constant. Thus (4) can be rewritten as(5)u^=arg⁡maxu⁡pf∣upu.Equation (5) can be rewritten as the following equation due to the monotonic logarithm function fz=ln⁡(z):(6)u^=arg⁡maxu⁡ln⁡pf∣u+ln⁡pu.Two probability density functions (PDF) should be established to get the true imageu^. 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 pf∣u=pε. The errors and noises are assumed to be Gaussian independent with mean zero, by which we can get the PDF of first term(7)pf∣u=1N1exp⁡-12f-Au-BTK-1f-Au-B,where N1∈R is a constant and K is a diagonal matrix with the diagonal elements being variances of each pixels. Rewriting (7) as(8)pf∣u=1N1exp⁡-12Qf-Au-B2,where Q is also a diagonal matrix with (9)qii=kii-1/2,where qii and kii 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 pu to preserve the reconstructed image edges and other detailed information [33]. The Huber-Markov prior is formulated as(10)pu=1N2exp⁡-12λ∑i,j∑c∈C,k=14ρdckui,j,where N2∈R is a constant, c is a clique within the set of all image cliques C, and ρ· is the Huber function defined as (11)ρz=z2,ifz≤μ,2μz-μ2,otherwise,where μ∈R is a threshold parameter used to separate the quadratic and linear regions. dckui,j can be computed by following finite second-order differences to measure the spatial activity for the pixel i,j. (12)dc1ui,j=ui-1,j-2ui,j+ui+1,j,dc2ui,j=ui,j-1-2ui,j+ui,j+1,dc3ui,j=12ui-1,j-1-2ui,j+ui+1,j+1,dc4ui,j=12ui-1,j+1-2ui,j+ui+1,j-1.Substituting (6) using (8) and (10), an equivalent formulation of (6) is obtained as follows: (13)u^=arg⁡minu⁡Eu,where(14)Eu=λ2Qf-Au-B2+∑i,j∑c∈C,k=14ρdckui,j.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.


## 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 fi,j and ui,j denote the DN at (i,j) in the contaminated image and the true image, respectively, the relationship between fi,j and ui,j can be represented as(1)fi,j=Ai,jui,j+Bi,j+εi,j,where Ai,j and Bi,j are the gain and offset parameters, respectively, and εi,j is the sum of linear error and sensor noise. According to (1), for healthy pixel (i,j), Ai,j=1 and Bi,j=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:(2)f=Au+B+ε,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.


## 2.2. The Shen-MAP Model

To get the true imageu by the given degraded image f, Shen et al. give the below estimation: (3)u^=arg⁡maxu⁡pu∣f.According to Bayes’ rule, we get(4)u^=arg⁡maxu⁡pf∣upupf.Because pu∣f is independent of f, pf can be regarded as a constant. Thus (4) can be rewritten as(5)u^=arg⁡maxu⁡pf∣upu.Equation (5) can be rewritten as the following equation due to the monotonic logarithm function fz=ln⁡(z):(6)u^=arg⁡maxu⁡ln⁡pf∣u+ln⁡pu.Two probability density functions (PDF) should be established to get the true imageu^. 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 pf∣u=pε. The errors and noises are assumed to be Gaussian independent with mean zero, by which we can get the PDF of first term(7)pf∣u=1N1exp⁡-12f-Au-BTK-1f-Au-B,where N1∈R is a constant and K is a diagonal matrix with the diagonal elements being variances of each pixels. Rewriting (7) as(8)pf∣u=1N1exp⁡-12Qf-Au-B2,where Q is also a diagonal matrix with (9)qii=kii-1/2,where qii and kii 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 pu to preserve the reconstructed image edges and other detailed information [33]. The Huber-Markov prior is formulated as(10)pu=1N2exp⁡-12λ∑i,j∑c∈C,k=14ρdckui,j,where N2∈R is a constant, c is a clique within the set of all image cliques C, and ρ· is the Huber function defined as (11)ρz=z2,ifz≤μ,2μz-μ2,otherwise,where μ∈R is a threshold parameter used to separate the quadratic and linear regions. dckui,j can be computed by following finite second-order differences to measure the spatial activity for the pixel i,j. (12)dc1ui,j=ui-1,j-2ui,j+ui+1,j,dc2ui,j=ui,j-1-2ui,j+ui,j+1,dc3ui,j=12ui-1,j-1-2ui,j+ui+1,j+1,dc4ui,j=12ui-1,j+1-2ui,j+ui+1,j-1.Substituting (6) using (8) and (10), an equivalent formulation of (6) is obtained as follows: (13)u^=arg⁡minu⁡Eu,where(14)Eu=λ2Qf-Au-B2+∑i,j∑c∈C,k=14ρdckui,j.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.


## 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:(15)u^=argminu⁡λ2Qf-Au-B2+α∇yu1+β∇2u1,where α, β are regularization parameters. ∇u=∇xu,∇yu=Dx+u,Dy+u and ∇2u=Dxx+u,Dyy+u,Dxy+u,Dyx+u 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 asm×n matrix. We denote the space Rm×n as H. The discrete gradient operator is a mapping ∇:H→V, where V=H×H. ∀u∈H, ∇u is given by (16)∇ui,j=Dx+ui,j,Dy+ui,j,with (17)Dx+ui,j=ui,j+1-ui,j,if1≤i≤m,1≤j<n,ui,1-ui,j,if1≤i≤m,j=n,Dy+ui,j=ui+1,j-ui,j,if1≤i<m,1≤j≤n,u1,j-ui,j,ifi=m,1≤j≤n.We useDx+i,j and Dy+i,j to denote forward difference operators with periodic boundary condition (u 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 (18)-div⁡p·u=p·∇u,∀u∈H,p∈V.This property forms the discrete analogue of integration by parts. We define (19)div⁡pi,j=Dx-p1i,j+Dy-p2i,j,p=p1,p2∈V,with (20)Dx-ui,j=ui,j-ui,n,if1≤i≤m,j=1,ui,j-ui,j-1,if1≤i≤m,1<j≤n,Dy-ui,j=ui,j-um,j,ifi=1,1≤j≤n,ui,j-ui-1,j,if1<i≤m,1≤j≤n.We denote backward difference operators with periodic boundary condition asDx-i,j and Dy-i,j. Analogously we define the second-order discrete differential operators as a mapping ∇2:H→W, where W=H×H×H×H. ∀u∈H, ∇2u is given by (21)∇2ui,j=Dxx+ui,j,Dyy+ui,j,Dxy+ui,j,Dyx+ui,j,with (22)Dxx+ui,j=ui,n-2ui,j+ui,j+1,if1≤i≤m,j=1,ui,j-1-2ui,j+ui,j+1,if1≤i≤m,1<j<n,ui,j-1-2ui,j+ui,1,if1≤i≤m,j=n,Dyy+ui,j=um,j-2ui,j+ui+1,j,ifi=1,1≤j≤n,ui-1,j-2ui,j+ui+1,j,if1<i<m,1≤j≤n,ui-1,j-2ui,j+u1,j,ifi=m,1≤j≤n,Dxy+ui,j=ui,j-ui+1,j-ui,j+1+ui+1,j+1,if1≤i<m,1≤j<n,ui,j-u1,j-ui,j+1+u1,j+1,ifi=m,1≤j<n,ui,j-ui+1,j-ui,1+ui+1,1,if1≤i<m,j=n,ui,j-u1,j-ui,1+u1,1,ifi=m,j=n.It is easy to know thatDxy+=Dx+Dy+=Dy+Dx+=Dyx+. The discrete second-order divergence operator is denoted as div2:W→H with the adjointness property (23)div2⁡q·u=q·∇2u,∀u∈H,q∈W.We define (24)div2⁡qi,j=Dxx-q11i,j+Dyy-q22i,j+Dxy-q12i,j+Dyx-q21i,j,where q=q11,q22,q12,q21∈W. We can check that Dxx-=Dxx+, Dyy-=Dyy+, and Dxy-=Dyx- with (25)Dxy-ui,j=ui,j-ui,n-um,j+um,n,ifi=1,j=1,ui,j-ui,j-1-um,j+um,j-1,ifi=1,1<j≤n,ui,j-ui-1,j-ui,n+ui-1,n,if1<i≤m,j=1,ui,j-ui,j-1-ui-1,j+ui-1,j-1,if1<i≤m,1<j≤n.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 f-B=C, and (15) can be written as (26)minuλ2QAu-C2+α∇yu1+β∇2u1.The corresponding constrained minimization problem is(27)ming,uv,wλ2QAu-C2+αv1+βw1,s.t.u=g,v=∇yg,w=∇2g.We can solve (27) by using the split Bregman method.(28)gk+1,uk+1,vk+1,wk+1=arg⁡ming,uv,w⁡λ2QAu-C2+αv1+βw1+λ12b1k+g-u2+λ22b2k+∇yg-v2+λ32b3k+∇2g-w2,(29)b1k+1=b1k+gk+1-uk+1,b2k+1=b2k+∇ygk+1-vk+1,b3k+1=b3k+∇2gk+1-wk+1.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].


#### 3.2.1.g Subproblem

From (28), we get(30)gk+1=arg⁡ming∈H⁡λ12b1k+g-uk2+λ22b2k+∇yg-vk2+λ32b3k+∇2g-wk2,whose optimality condition is(31)λ1I-λ2Dy-Dy++λ3Dxx-Dxx++2Dxy-Dxy++Dyy-Dyy+g=λ1uk-b1k-λ2Dy-vyk-b2yk+λ3Dxx-wxxk-b3xxk+2Dxy-wxyk-b3xyk+Dyy-wyyk-b3yyk.By choosing periodic boundary conditions, (31) can be solved by fast Fourier transform (FFT denoted by F(·)) and get (32)FD·Fg=FEQright,where (33)FD=λ1I-λ2FDy-Dy++λ3FDxx-Dxx++2FDxy-Dxy++FDyy-Dyy+,and EQright is the right-hand side of (31). Here “·” means pointwise multiplication of matrices. Thus, we have (34)gk+1=F-1FEQright/FD,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 solveg subproblem.


#### 3.2.2.u Subproblem

According to (28), we have(35)uk+1=arg⁡minu∈H⁡λ2QAu-C2+λ12b1k+gk+1-u2.It leads to (36)u=λATQTQA+λ1I-1λATQTQC+λ1gk+1+b1k.This is a closed form solution, so (35) can be solved simply.


#### 3.2.3.v Subproblem

Due to (28), we have(37)vk+1=arg⁡minv∈V⁡αv1+λ22b2k+∇ygk+1-v2.It can be solved by shrinkage [35]. Denoting ski,j=b2k+Dygk+1i,j∈R, we can get (38)vk+1i,j=max⁡ski,j-αλ2,0ski,jski,j,where 0/0=0 is assumed.


#### 3.2.4.w Subproblem

Because of (28), we have(39)wk+1=arg⁡minw∈W⁡βw1+λ32b3k+∇2gk+1-w2.We solve (39) adopting the same procedure as in v-subproblem and set (40)tki,j=t1ki,j,t2ki,j,t3ki,j,t3ki,j=b3k+∇2gk+1i,j=b3xxki,j+Dxx+gk+1i,j,b3yyki,j+Dyy+gk+1i,j,b3xyki,j+Dxy+gk+1i,j,b3xyki,j+Dxy+gk+1i,j∈R4,and we obtain the result as follows: (41)wxxk+1i,j=max⁡tki,j-βλ3,0t1ki,jtki,j,wyyk+1i,j=max⁡tki,j-βλ3,0t2ki,jtki,j,wxyk+1i,j=max⁡tki,j-βλ3,0t3ki,jtki,j,where 0/0=0 is assumed. We define the stopping criteria as follows: (42)uk+1-uk2uk2≤η, where η∈R is a predetermined positive value. A direct implementation of the USTV model is presented in Algorithm 1.Algorithm 1:The USTV model solved by the split Bregman method.
Input:u1=f. Initialization:u0=1, g1=u1, v1=∇yg1, w1=∇2g1, b11=1, b21=1, b31=1. k = 0. while Terminate condition is not satisfied do k=k+1. The  g  subproblem: compute gk+1 according to Eq. (30) by given uk, vk, wk. The  u  subproblem: compute uk+1 according to Eq. (35) by given gk+1. The  v  subproblem: compute vk according to Eq. (37) by given gk+1. The  w  subproblem: compute wk according to Eq. (39) by given gk+1. Update on  b1: b1k+1=b1k+gk+1-uk+1. Update on  b2: b2k+1=b2k+∇ygk+1-vk+1. Update on  b3: b3k+1=b3k+∇2gk+1-wk+1. end while Output: Destriping imageuk.


### 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(43)g,u,v,wT=arg⁡ming,uv,w⁡λ2QAu-C2+αv1+βw1+λ1b1λ2b2λ3b3,Ig∇y∇2g-Iu000Iv000Iwuvw+12λ1Igλ2∇yλ3∇2g-λ1Iu000λ2Iv000λ3Iwuvw2.The optimization problem is well structured since the four variablesg, u, v, and w can be separated into two groups g and u,v,wT and the convergence of the method can be guaranteed by the existing ADMM theory [36]. Moreover, the variables u, v, and w 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.


## 3.1. Notations

Without loss of generality, we represent a gray image asm×n matrix. We denote the space Rm×n as H. The discrete gradient operator is a mapping ∇:H→V, where V=H×H. ∀u∈H, ∇u is given by (16)∇ui,j=Dx+ui,j,Dy+ui,j,with (17)Dx+ui,j=ui,j+1-ui,j,if1≤i≤m,1≤j<n,ui,1-ui,j,if1≤i≤m,j=n,Dy+ui,j=ui+1,j-ui,j,if1≤i<m,1≤j≤n,u1,j-ui,j,ifi=m,1≤j≤n.We useDx+i,j and Dy+i,j to denote forward difference operators with periodic boundary condition (u 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 (18)-div⁡p·u=p·∇u,∀u∈H,p∈V.This property forms the discrete analogue of integration by parts. We define (19)div⁡pi,j=Dx-p1i,j+Dy-p2i,j,p=p1,p2∈V,with (20)Dx-ui,j=ui,j-ui,n,if1≤i≤m,j=1,ui,j-ui,j-1,if1≤i≤m,1<j≤n,Dy-ui,j=ui,j-um,j,ifi=1,1≤j≤n,ui,j-ui-1,j,if1<i≤m,1≤j≤n.We denote backward difference operators with periodic boundary condition asDx-i,j and Dy-i,j. Analogously we define the second-order discrete differential operators as a mapping ∇2:H→W, where W=H×H×H×H. ∀u∈H, ∇2u is given by (21)∇2ui,j=Dxx+ui,j,Dyy+ui,j,Dxy+ui,j,Dyx+ui,j,with (22)Dxx+ui,j=ui,n-2ui,j+ui,j+1,if1≤i≤m,j=1,ui,j-1-2ui,j+ui,j+1,if1≤i≤m,1<j<n,ui,j-1-2ui,j+ui,1,if1≤i≤m,j=n,Dyy+ui,j=um,j-2ui,j+ui+1,j,ifi=1,1≤j≤n,ui-1,j-2ui,j+ui+1,j,if1<i<m,1≤j≤n,ui-1,j-2ui,j+u1,j,ifi=m,1≤j≤n,Dxy+ui,j=ui,j-ui+1,j-ui,j+1+ui+1,j+1,if1≤i<m,1≤j<n,ui,j-u1,j-ui,j+1+u1,j+1,ifi=m,1≤j<n,ui,j-ui+1,j-ui,1+ui+1,1,if1≤i<m,j=n,ui,j-u1,j-ui,1+u1,1,ifi=m,j=n.It is easy to know thatDxy+=Dx+Dy+=Dy+Dx+=Dyx+. The discrete second-order divergence operator is denoted as div2:W→H with the adjointness property (23)div2⁡q·u=q·∇2u,∀u∈H,q∈W.We define (24)div2⁡qi,j=Dxx-q11i,j+Dyy-q22i,j+Dxy-q12i,j+Dyx-q21i,j,where q=q11,q22,q12,q21∈W. We can check that Dxx-=Dxx+, Dyy-=Dyy+, and Dxy-=Dyx- with (25)Dxy-ui,j=ui,j-ui,n-um,j+um,n,ifi=1,j=1,ui,j-ui,j-1-um,j+um,j-1,ifi=1,1<j≤n,ui,j-ui-1,j-ui,n+ui-1,n,if1<i≤m,j=1,ui,j-ui,j-1-ui-1,j+ui-1,j-1,if1<i≤m,1<j≤n.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 f-B=C, and (15) can be written as (26)minuλ2QAu-C2+α∇yu1+β∇2u1.The corresponding constrained minimization problem is(27)ming,uv,wλ2QAu-C2+αv1+βw1,s.t.u=g,v=∇yg,w=∇2g.We can solve (27) by using the split Bregman method.(28)gk+1,uk+1,vk+1,wk+1=arg⁡ming,uv,w⁡λ2QAu-C2+αv1+βw1+λ12b1k+g-u2+λ22b2k+∇yg-v2+λ32b3k+∇2g-w2,(29)b1k+1=b1k+gk+1-uk+1,b2k+1=b2k+∇ygk+1-vk+1,b3k+1=b3k+∇2gk+1-wk+1.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].


### 3.2.1.g Subproblem

From (28), we get(30)gk+1=arg⁡ming∈H⁡λ12b1k+g-uk2+λ22b2k+∇yg-vk2+λ32b3k+∇2g-wk2,whose optimality condition is(31)λ1I-λ2Dy-Dy++λ3Dxx-Dxx++2Dxy-Dxy++Dyy-Dyy+g=λ1uk-b1k-λ2Dy-vyk-b2yk+λ3Dxx-wxxk-b3xxk+2Dxy-wxyk-b3xyk+Dyy-wyyk-b3yyk.By choosing periodic boundary conditions, (31) can be solved by fast Fourier transform (FFT denoted by F(·)) and get (32)FD·Fg=FEQright,where (33)FD=λ1I-λ2FDy-Dy++λ3FDxx-Dxx++2FDxy-Dxy++FDyy-Dyy+,and EQright is the right-hand side of (31). Here “·” means pointwise multiplication of matrices. Thus, we have (34)gk+1=F-1FEQright/FD,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 solveg subproblem.


### 3.2.2.u Subproblem

According to (28), we have(35)uk+1=arg⁡minu∈H⁡λ2QAu-C2+λ12b1k+gk+1-u2.It leads to (36)u=λATQTQA+λ1I-1λATQTQC+λ1gk+1+b1k.This is a closed form solution, so (35) can be solved simply.


### 3.2.3.v Subproblem

Due to (28), we have(37)vk+1=arg⁡minv∈V⁡αv1+λ22b2k+∇ygk+1-v2.It can be solved by shrinkage [35]. Denoting ski,j=b2k+Dygk+1i,j∈R, we can get (38)vk+1i,j=max⁡ski,j-αλ2,0ski,jski,j,where 0/0=0 is assumed.


### 3.2.4.w Subproblem

Because of (28), we have(39)wk+1=arg⁡minw∈W⁡βw1+λ32b3k+∇2gk+1-w2.We solve (39) adopting the same procedure as in v-subproblem and set (40)tki,j=t1ki,j,t2ki,j,t3ki,j,t3ki,j=b3k+∇2gk+1i,j=b3xxki,j+Dxx+gk+1i,j,b3yyki,j+Dyy+gk+1i,j,b3xyki,j+Dxy+gk+1i,j,b3xyki,j+Dxy+gk+1i,j∈R4,and we obtain the result as follows: (41)wxxk+1i,j=max⁡tki,j-βλ3,0t1ki,jtki,j,wyyk+1i,j=max⁡tki,j-βλ3,0t2ki,jtki,j,wxyk+1i,j=max⁡tki,j-βλ3,0t3ki,jtki,j,where 0/0=0 is assumed. We define the stopping criteria as follows: (42)uk+1-uk2uk2≤η, where η∈R is a predetermined positive value. A direct implementation of the USTV model is presented in Algorithm 1.Algorithm 1:The USTV model solved by the split Bregman method.
Input:u1=f. Initialization:u0=1, g1=u1, v1=∇yg1, w1=∇2g1, b11=1, b21=1, b31=1. k = 0. while Terminate condition is not satisfied do k=k+1. The  g  subproblem: compute gk+1 according to Eq. (30) by given uk, vk, wk. The  u  subproblem: compute uk+1 according to Eq. (35) by given gk+1. The  v  subproblem: compute vk according to Eq. (37) by given gk+1. The  w  subproblem: compute wk according to Eq. (39) by given gk+1. Update on  b1: b1k+1=b1k+gk+1-uk+1. Update on  b2: b2k+1=b2k+∇ygk+1-vk+1. Update on  b3: b3k+1=b3k+∇2gk+1-wk+1. end while Output: Destriping imageuk.


## 3.2.1.g Subproblem

From (28), we get(30)gk+1=arg⁡ming∈H⁡λ12b1k+g-uk2+λ22b2k+∇yg-vk2+λ32b3k+∇2g-wk2,whose optimality condition is(31)λ1I-λ2Dy-Dy++λ3Dxx-Dxx++2Dxy-Dxy++Dyy-Dyy+g=λ1uk-b1k-λ2Dy-vyk-b2yk+λ3Dxx-wxxk-b3xxk+2Dxy-wxyk-b3xyk+Dyy-wyyk-b3yyk.By choosing periodic boundary conditions, (31) can be solved by fast Fourier transform (FFT denoted by F(·)) and get (32)FD·Fg=FEQright,where (33)FD=λ1I-λ2FDy-Dy++λ3FDxx-Dxx++2FDxy-Dxy++FDyy-Dyy+,and EQright is the right-hand side of (31). Here “·” means pointwise multiplication of matrices. Thus, we have (34)gk+1=F-1FEQright/FD,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 solveg subproblem.


## 3.2.2.u Subproblem

According to (28), we have(35)uk+1=arg⁡minu∈H⁡λ2QAu-C2+λ12b1k+gk+1-u2.It leads to (36)u=λATQTQA+λ1I-1λATQTQC+λ1gk+1+b1k.This is a closed form solution, so (35) can be solved simply.


## 3.2.3.v Subproblem

Due to (28), we have(37)vk+1=arg⁡minv∈V⁡αv1+λ22b2k+∇ygk+1-v2.It can be solved by shrinkage [35]. Denoting ski,j=b2k+Dygk+1i,j∈R, we can get (38)vk+1i,j=max⁡ski,j-αλ2,0ski,jski,j,where 0/0=0 is assumed.


## 3.2.4.w Subproblem

Because of (28), we have(39)wk+1=arg⁡minw∈W⁡βw1+λ32b3k+∇2gk+1-w2.We solve (39) adopting the same procedure as in v-subproblem and set (40)tki,j=t1ki,j,t2ki,j,t3ki,j,t3ki,j=b3k+∇2gk+1i,j=b3xxki,j+Dxx+gk+1i,j,b3yyki,j+Dyy+gk+1i,j,b3xyki,j+Dxy+gk+1i,j,b3xyki,j+Dxy+gk+1i,j∈R4,and we obtain the result as follows: (41)wxxk+1i,j=max⁡tki,j-βλ3,0t1ki,jtki,j,wyyk+1i,j=max⁡tki,j-βλ3,0t2ki,jtki,j,wxyk+1i,j=max⁡tki,j-βλ3,0t3ki,jtki,j,where 0/0=0 is assumed. We define the stopping criteria as follows: (42)uk+1-uk2uk2≤η, where η∈R is a predetermined positive value. A direct implementation of the USTV model is presented in Algorithm 1.Algorithm 1:The USTV model solved by the split Bregman method.
Input:u1=f. Initialization:u0=1, g1=u1, v1=∇yg1, w1=∇2g1, b11=1, b21=1, b31=1. k = 0. while Terminate condition is not satisfied do k=k+1. The  g  subproblem: compute gk+1 according to Eq. (30) by given uk, vk, wk. The  u  subproblem: compute uk+1 according to Eq. (35) by given gk+1. The  v  subproblem: compute vk according to Eq. (37) by given gk+1. The  w  subproblem: compute wk according to Eq. (39) by given gk+1. Update on  b1: b1k+1=b1k+gk+1-uk+1. Update on  b2: b2k+1=b2k+∇ygk+1-vk+1. Update on  b3: b3k+1=b3k+∇2gk+1-wk+1. end while Output: Destriping imageuk.


## 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(43)g,u,v,wT=arg⁡ming,uv,w⁡λ2QAu-C2+αv1+βw1+λ1b1λ2b2λ3b3,Ig∇y∇2g-Iu000Iv000Iwuvw+12λ1Igλ2∇yλ3∇2g-λ1Iu000λ2Iv000λ3Iwuvw2.The optimization problem is well structured since the four variablesg, u, v, and w can be separated into two groups g and u,v,wT and the convergence of the method can be guaranteed by the existing ADMM theory [36]. Moreover, the variables u, v, and w 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 matricesA (gains), B (biases), and Q 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(44)A=σsσr,B=μs-μrσsσr,where μs, μr are the mean values of the striped row and reference row, respectively, and σs, σr are the standard deviations of the striped row and reference row, respectively. The elements in Q represent the reciprocals of the noise standard deviation in different pixel location and they can be determined by(45)qii=ln⁡e-1std-minmax-min+1,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 A and Q 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 η=1.0×10-6, λ=5×104, α=1, β=0.8, λ1=10max⁡(α,β), λ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 1
MODIS images before destriping. (a) Aqua MODIS band 30. (b) Terra MODIS band 28. (c) Terra MODIS band 30.
(a) (b) (c)Figure 2
Destriped results of Aqua MODIS band 30. (a) Original striped image. (b) Moment matching [9]. (c) The Shen-MAP model [12]. (d) The proposed USTV model.
(a) (b) (c) (d)Figure 3
Destriped results of Terra MODIS band 28. (a) Original striped image. (b) Moment matching [9]. (c) The Shen-MAP model [12]. (d) The proposed USTV model.
(a) (b) (c) (d)Figure 4
Destriped results of Terra MODIS band 30. (a) Original striped image. (b) Moment matching [9]. (c) The Shen-MAP model [12]. (d) The proposed USTV model.
(a) (b) (c) (d)Figure5 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.Figure 5
Mean cross-track profiles of MODIS images before (blue line) and after (red line) destriping. (a) Aqua MODIS band 30. (b) Terra MODIS band 28. (c) Terra MODIS band 30.
(a) (b) (c)Figure 6
Mean column power spectrum of the MODIS images before (blue line) and after (red line) destriping. (a) Aqua MODIS band 30. (b) Terra MODIS band 28. (c) Terra MODIS band 30.
(a) (b) (c)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(46)ICV=RmRsd,where Rm and Rsd 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(47)NR=NsNde, where Ns is the power of frequency components produced by stripes in the striped image and Nde for the destriped image. Ns and Nde can be computed by(48)Ns=∑BWNPsDis,Nde=∑BWNPdeDis,Dis∈0.1,0.2,0.3,0.4,0.5,where Dis is the distance from the origin in the Fourier space and BWN 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.Table 1
ICVs and NRs of the comparison experiments.
Original Moment Shen-MAP USTV Aqua band 30 Sample 1 (ICV) 8.17 29.51 37.96 66.23 Sample 2 (ICV) 13.01 23.14 25.29 28.81 Sample 3 (ICV) 8.78 14.84 18.45 23.64 NR 1.00 2.41 2.59 2.64 Terra band 28 Sample 1 (ICV) 10.94 20.59 25.18 39.78 Sample 2 (ICV) 25.31 22.10 38.29 69.35 Sample 3 (ICV) 15.19 14.09 16.77 24.76 NR 1.00 2.92 3.64 4.62 Terra band 30 Sample 1 (ICV) 31.45 75.41 80.88 97.36 Sample 2 (ICV) 39.91 85.55 105.29 180.87 Sample 3 (ICV) 33.26 52.84 67.52 94.05 NR 1.00 1.44 1.60 1.71


## 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 A and B 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.



---

*Source: 4397189-2017-04-09.xml*