**Directional** ` _{0} **Sparse Modeling for Image Stripe** **Noise Removal**

**Hong-Xia Dou, Ting-Zhu Huang *, Liang-Jian Deng, Xi-Le Zhao, and Jie Huang**

School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu, Sichuan, 611731, P. R. China

***** Corresponding author: tingzhuhuang@126.com
Version December 26, 2017 submitted to MDPI

**Abstract:**Remote sensing images are often polluted by stripe noise, which leads to negative impact

1

on visual performance. Thus, it is necessary to remove stripe noise for the subsequent applications,

2

e.g.,classification, recognition,etc. This paper commits to remove the stripe noise to enhance the

3

visual quality of images, in the meanwhile preserves image details of stripe-free regions. Instead

4

of solving the underlying image by various algorithms, we first estimate the stripe noise from the

5

degraded images, then computing the final destriping image by the difference of the known stripe

6

image and the estimated stripe noise. In this paper, we propose a non-convex`_{0}sparse model for

7

remote sensing image destriping by taking full consideration of the intrinsically directional and

8

structural priors of stripe noise, as well as the locally continuous property of underlying image.

9

Moreover, the proposed non-convex model is solved by a proximal alternating direction method of

10

multipliers (PADMM) based algorithm and we also give the corresponding theoretical analysis of the

11

proposed algorithm. Extensive experimental results on simulated and real data demonstrate that

12

the proposed method outperforms recently state-of-the-art destriping methods, both visually and

13

quantitatively.

14

**Keywords:** Non-convex`_{0}sparse model; PADMM based algorithm; Mathematical program with

15

equilibrium constraints (MPEC); Stripe noise removal.

16

**1. Introduction**

17

Stripe noise (all denoted as “stripes” in this paper), which is generally caused by the inconsistence

18

of the detecting element scanning or the influence of the detector moving and temperature changes,etc.,

19

are an universal phenomenon in remote sensing images. They will result in a bad influence not only on

20

visual quality but also on subsequent applications in remote sensing images. Therefore, it is necessary

21

to remove stripes and simultaneously maintain the healthy pixels from the degraded images. In

22

general, the stripes have strongly directional and structural information,e.g.,pixels normally damaged

23

on row by row or column by column.

24

Recently, many approaches for destriping problems have been proposed, which may be roughly

25

divided into three categories, mainly including filtering-based methods, statistics-based methods

26

and optimization-based methods. Note that the proposed method belongs to the category of

27

optimization-based methods.

28

The filtering-based methods, which are easy to obtain the results with various filters, have been

29

widely utilized for remote sensing image destriping, see [1–4]. In [1], Chenet al.propose an approach

30

for remote sensing image destriping tasks based on a finite-impulse response filter (FIR) in frequency

31

domain, as well as exhibit the results on the experimental CMODIS data. However, the given method

32

unavoidably leads to ringing and ripple artifacts. In [3], the wavelet analysis and adaptive fourier

33

Submitted toMDPI, pages 1 – 25 www.mdpi.com/journal/notspecified

zero-frequency amplitude normalization are used for hyperspectral image destriping problems, and

34

this wavelet-based method shows promising ability for both stripes and random noise.

35

The statistics-based methods are mainly to analyze the distribution of stripes. These approaches

36

hold strong directional characters, to formulate excellent priors for the remote sensing image destriping,

37

e.g.,[5–11]. In [7], Weinrebet al.introduce a method based on matching empirical distribution functions

38

(EDFs) for GOES-7 data, while the limitations and unstable property are caused by assuming the

39

similarity and regularity among the stripes. To conquer the instability when the stripes are irregular or

40

nonlinear, Rakwatinet al.[9] introduce a method, using both histogram-matching algorithm and local

41

least squares fitting, to remove the stripes of Aqua MODIS band 6. In [10], spectral moment matching

42

(SpcMM) method, which can remove various frequencies stripes in a specific band automatically, is

43

proposed for Hyperion image destriping. In addition, Shenet al.[11] employ a piece-wise destriping

44

method, which uses correction coefficients of each portion by considering neighbouring normal row,

45

for nonlinear and irregular stripes, but it can not automatically select a threshold to divide the image

46

into different parts.

47

Recently, the optimization-based methods show superiorities for remote sensing image destriping

48

problems,e.g.,[12–23]. The image destriping generally results in an ill-posed problem which fails

49

to obtain a meaningful, stable and unique solution. Therefore, a common strategy for ill-posed

50

problem is to construct a regularization model via investigating the underlying image priors. For

51

the optimization-based methods, they focus on searching and discovering the intrinsically prior

52

knowledge to generate reasonable regularization models. In [17], the authors present a unidirectional

53

total variational (UTV) model for MODIS image stripes removal by fully considering the directional

54

information of stripes. The UTV model is motivated by the classical TV model and the analysis of

55

directional stripes. Changet al.[21] propose an optimization model combining the UTV with sparse

56

priors of stripes applying to denoising and destriping simultaneously. In [22], the authors utilize the

57

split Bregman iteration method with an anisotropic spectral-spatial total variation regularization to

58

remove multispectral image stripes.

59

In summary, although these optimization-based methods can yield excellent results of removing

60

stripes, there still exists much room to improve. Most of them are implemented only from the

61

perspective of noise removal, but without considering the typical properties of stripes,e.g.,directional

62

and structural properties. Even though considering these properties, the formulated sparse destriping

63

models fail to accurately depict the typical properties of stripes, see [24], [25]. Moreover, the designed

64

algorithms for non-convex models,e.g.,`_{0}sparse model, can not obtain the most precise solution.

65

These motivate us to develop a more reasonable model and effectively design the corresponding

66

algorithm, which theoretically guarantees the convergence, to solve the remote sensing destriping

67

problems.

68

In this paper, to remove the stripes of remote sensing images, we propose a non-convex sparse

69

model which mainly consists of three sparse priors, including an`_{0}sparse prior by fully considering

70

the directional property of stripes (y-axis), an`_{1}sparse prior by considering the discontinuity of

71

underlying image (x-axis), and the sparsity of stripes by considering the structural property of stripes.

72

Moreover, we design a PADMM based algorithm to solve the proposed non-convex sparse model. In

73

particular, the convergence to the KKT point of the optimization problem is theoretically proven in

74

the work. Results of several simulated and real images show that the proposed method is superior to

75

recently state-of-the-art destriping methods.

76

The contributions of this work are summarized as follows

77

1) Fully considering the latent priors of stripes, we formulate an`_{0}sparse model which depicts

78

the intrinsically sparse character more accurately than`_{1}sparse model.

79

2) We solve the non-convex model by a designed PADMM based algorithm which we have given

80

the corresponding theoretical analysis of the proposed algorithm by this paper (see Appendix A).

81

3) The proposed method, which is less sensitive to related parameters, outperforms recently

82

several state-of-the-art image destriping methods.

83

The outline of this paper is organized as follows. In Section 2, we will briefly introduce the related

84

work. The proposed model and detailed solving algorithm will be shown in Section 3. In section 4, we

85

compare the proposed method with some state-of-the-art remote sensing image destriping methods,

86

and discuss the results with different stripes. Finally, conclusions are drawn in Section 5.

87

**2. Related work**

88

2.1. Destriping problem formulation

89

The striping effects in remote sensing images mainly make up of additive and multiplicative components [15]. However, the multiplicative stripes can be described as additive case by the logarithm [26]. Thus, many researches more focus on the additive stripes model

**b**(x,y) =**u**(x,y) +**s**(x,y) (1)
where **b**(x,y)_{,} _{u}(x,y) _{and} _{s}(x,y) separately denote the components of the observe image, the
underlying image and stripes at the location(x,y). For convenience, a matrix-vector form can be
written as follows

**b**=**u**+**s**, (2)

where **b**,**u** and**s** ∈ R^{n} represent the lexicographical order vectors of **b**(x,y), **u**(x,y) and**s**(x,y),

90

respectively. The purpose of our work is to estimate the stripes**s**, then the underlying image will be

91

recovered by the formula of**u**=**b**−**s**.

92

2.2. UTV for remote sensing image destriping

93

The total variation (TV) model, which is first proposed by Rudin, Oshaer and Fatemi (ROF) [27], has shown powerful ability in many image applications,e.g.,image unmixing [28], image deblurring [29], image inpainting [30],etc.It has the following form

E(**u**) = ^{1}
2
Z

Ω||**u**−**b**||^{2}+*λ*TV(**u**), (3)
where*λ*is a positive regularization parameter, andTV(**u**)represents the regularization expressed as

TV(u) = Z

Ω|∇**u**|=
Z

Ω

s
*∂***u**

*∂*x
2

+
*∂***u**

*∂*y
2

dxdy. (4)

In many approaches, **s**(x,y) is usually regarded as constant in a given line. Although this

94

assumption has shown stability in MOS-B, it fails in MODIS. Not only predominant nonlinear effects,

95

but also the data quality of random stripes have been obtained in many emissive bands. Thus, more

96

realistic assumptions are introduced to design an efficient destriping method.

97

Without loss the generality, we can assume that the stripes are along the vertical direction (y-axis).

Fully considering the directional property of stripes, the authors in [17] consider the following relation

*∂***s**(x,y)

*∂*y

*∂***s**(x,y)

*∂*x

, (5) where we denote y-axis is along stripes direction, and x-axis is across stripes direction. By the relation in Eq. (5), we have

Z

Ω

*∂***s**(x,y)

*∂*y

dxdy Z

Ω

*∂***s**(x,y)

*∂*x

dxdy, (6)

which means

TVy(**s**)TVx(**s**) (7)

whereTVxandTVyare horizontal and vertical variations, respectively. The authors in [17] encourage the robustness of stripes removal by minimizing the unidirectional total variation (UTV) model as follows

E(**u**) =TVy(**u**−**b**) +*λ*TVx(**u**), (8)

which can be solved by Euler-Lagrange equation based algorithm.

98

In [17], the UTV model can effectively deal with remote sensing image destriping problems,

99

which has been demonstrated holding promising ability on Aqua and Terra MODIS data. Although

100

TV model preserves image edges well, it can not accurately depict the specifically directional property

101

of stripes, and leads to undesired results. The UTV model that involves unidirectional constraint can

102

remove stripes excellently in the meanwhile not destroy the underlying image details. Inspired by

103

the UTV model, we fully consider the intrinsically directional and structural priors of stripes and the

104

continuous property of the underlying image. Finally, we form a unidirectional and sparse based

105

optimization model.

106

**3. The proposed method**

107

Combining the stripes model (2), we will give the proposed optimal model with unidirectional

108

prior motivated by the extension of the UTV model. In what follows, the detailed explanations of the

109

proposed model and the corresponding solving algorithm will be exhibited.

110

3.1. The proposed model

111

3.1.1. Local smoothness along stripe direction

112

The stripes of remote sensing images are generally appeared with column-by-column (y-axis)
or row-by-row (x-axis), without loss of generality, we view all stripes as column-by-column case to
formulate the finally directional model^{1}. Considering the smoothness within the stripes, the difference
between adjacent pixels is quite small, or even close to zero, thus we generally use sparse prior for
this character along the stripe direction (y-axis). The first regularization for the difference within the
stripes is given as follows

R_{1}=||∇_{y}**s**||_{0}, (9)

where∇_{y}is a partial difference operator along stripe direction^{2}. Comparing with some popular

113

sparse measures,e.g.,`_{1}-norm and`_{p}-norm (0 < p<1), the`_{0}-norm that stands for the number of

114

non-zero elements of a vector is the most accurate measure to depict sparse property, thus here we

115

employ`_{0}-norm to describe the sparsity of∇ys. Although this term will lead to the non-convexity

116

of the proposed model, we utilize the designed PADMM based algorithm to guarantee the solution

117

converging to the KKT point.

118

3.1.2. Local continuity of the underlying image

119

In general, the underlying image**u**along x-axis is viewed as being continuous. When adding
column-by-column stripes**s**to the underlying image, the local continuity of**u**is broken, which means
that we should force∇x**u**being small to keep the continuity of**u**. By this assumption and the relation
**u**=**b**−**s**, we utilize the following`_{1}-norm regularization to describe the local continuity of the
underlying image

R_{2}=||∇_{x}(_{b}−_{s})||_{1}, (10)

1 The row-by-row stripes can be easily rotated to column-by-column stripes to fit in the proposed model.

2 ∇_{y}urepresents the vector form of∇_{y}UwhereUis a 2D image. The similar meaning is∇_{x}u.

where∇_{x} represents the difference operator in the across-stripe direction. Note that this term is

120

actually the second term of the UTV model (8).

121

3.1.3. Global sparsity of stripes

122

In many destriping approaches,e.g.,[24,25,31,32], the stripes can be naturally viewed as being
sparse when the stripes are not heavy. Inspired by their excellent works, here we take the`_{1}-norm to
depict the sparsity of stripes, see as follows

R3=||**s**||_{1}. (11)

Even though the stripes are heavy, this sparse term (11) is still necessary to retain, since it can

123

effectively avoid the undesired effect and keep the robustness of the proposed method (see more

124

discussion from the results section).

125

Combining the above three regularization terms, we finally formulate the`_{0}sparse model for
remote sensing image destriping,

min**s** ||∇_{y}_{s}||_{0}+*µ*||_{s}||_{1}+*λ*||∇_{x}(_{b}−_{s})||_{1}, (12)
where*µ*and*λ*are two positive parameters.

126

(a) (b)

**Figure 1.**The number of nonzero of**s**(a) and∇_{y}**s**(b), where**s**is estimated from a real image example
(see Fig. 4) by the method [24]. It is clear that∇_{y}**s**is more sparse than s.

Note that, the proposed model (12) is similar as the model in [24], since they both employ the

127

directional property of stripes. However, there still exists an important difference that the model in [24]

128

enforces`_{1}norm to∇_{y}**s**and`_{0}norm to**s**whereas our model enforces`_{1}norm to s and`_{0}norm to∇_{y}**s**.

129

It can be seen that our model is more reasonable than the model in [24], because∇_{y}_{s}is significantly

130

more sparse than**s**. For instance, Fig.1shows the number of non-zeros of**s**(Fig.1(a)) and∇y**s**(Fig.

131

1(b)), where**s**is estimated from a real image example by the method [24], it is clear that∇_{y}**s**is almost

132

all around 0, whereas**s**is not. The`_{0}norm is the best way to depict sparsity, thus our model which

133

enforces`_{0}_{norm to}∇y**s**.

134

In what follows, we will exhibit how to solve the proposed non-convex sparse model by

135

introducing the PADMM based algorithm, as well as give the theoretical analysis of the convergence.

136

3.2. The solution

137

Before solving the proposed model (12), we first present an excellent work,i.e.,Mathematical

138

program with equilibrium constraints (MPEC) [31], to transfer the non-convex`_{0}regularization term

139

to the other equivalent one.

140 141

**Equivalent MPEC reformulation:** For the non-convex `_{0} regularization term, there exist many
approaches to approximate it, e.g., `_{1}-norm [33], the logarithm function [34] or the penalty

decomposition algorithm (PDA) [35]. In this work, we are inspired by a recently elegant work,
i.e.,MPEC, to transfer the`_{0}regularization term to an equivalent problem, so that we can design a
PADMM based algorithm to efficiently solve the equivalent model, in the meanwhile theoretically
guarantee the convergence.**Lemma:** [MPEC equation [31]] For any given**w**∈R^{n}, it holds that

||**w**||_{0}= min

**0**≤**v**≤**1**h**1**,**1**−**v**i, s.t.**v** |**w**|=0, (13)
and**v**^{∗}=**1**−sign(|**w**|)is the unique optimal solution of the minimization problem (13).

142

**Proof:**See details in [31].

143

From Lemma3.2, the`_{0}-norm sparse optimization model in Eq. (12) is equivalent to

**0**≤**v**≤**1**,**s**min h**1**,**1**−**v**i+*µ*||**s**||_{1}+*λ*||∇_{x}(**b**−**s**)||_{1},
s.t. **v** |∇_{y}**s**|=0,

(14)

where denotes the elementwise product. According to the analysis of [31], if**s**^{∗} is the globally

144

optimal solution of Eq. (12), then(_{s}^{∗}_{,}_{1}−sign(|∇y**s**^{∗}|))is the unique global minimizer of Eq. (14).

145

Note that the Eq. (14) is still a non-convex problem, and the non-convexity is only caused by the

146

constraint**v** |∇_{y}**s**|=0. However, this problem (14) is similar to the main problem in [31], which

147

is efficiently solved by a PADMM^{3}based algorithm that theoretically guarantees the convergence.

148

Therefore, we employ the designed PADMM based algorithm to solve the resulted problem (14), as

149

well as give the theoretical analysis of the convergence.

150

In the following, we will use the PADMM based algorithm to solve the optimization problem (14).

151

3.3. PADMM based Algorithm

152

Considering the non-smooth`_{1}terms in problem (14), we take the following variable substitutions
to get the new optimization problem,

**0**≤**v**≤**1**,smin h**1**,**1**−**v**i+*µ*||**z**||_{1}+*λ*||**w**||_{1},

s.t **v** |**h**|=0,∇_{y}**s**=**h**,**s**=**z**,∇_{x}(**b**−**s**) =**w**,

(15)

with the auxiliary variables**h**,**z**,**w** ∈ R^{n}. The augmented Lagrangian functionLof Eq. (15) is as
follows

L(**h**,**z**,**w**,**v**,**s**,**π**_{1},**π**_{2},**π**_{3},**π**_{4},*β*_{1},*β*_{2},*β*_{3},*β*_{4})

=h**1**,**1**−**v**i+*µ*||**z**||_{1}+*λ*||**w**||_{1}+h∇_{y}**s**−**h**,**π**_{1}i

+^{β}^{1}

2 ||∇_{y}**s**−**h**||^{2}_{2}+h**s**−**z**,* π*2i+

^{β}

^{2}

2 ||**s**−**z**||^{2}_{2}
+h∇_{x}(**b**−**s**)−**w**,* π*3i+

^{β}

^{3}

2 ||∇_{x}(**b**−**s**)−**w**||^{2}_{2}
+h**v** |**h**|,**π**_{4}i+^{β}^{4}

2 ||**v** |**h**|||^{2}_{2},

(16)

where* π*1,

*2,*

**π***3and*

**π***4are Lagrange multipliers, and*

**π***β*

_{1},

*β*

_{2},

*β*

_{3}and

*β*

_{4}are positive parameters. The

153

minimization problem (16) can be solved by the PADMM based algorithm. Next, we discuss the

154

solution of each subproblem.

155

3 Actually, PADMM method is an extended version of ADMM method, which has been applicated to many image applications, e,g.,image deblurring [36], image denoising [37], tensor completion [38],etc.

1) The**h**-subproblem can be written to the minimized problem as follows
min**h** h∇_{y}_{s}^{k}−**h**,**π**_{1}^{k}i+ ^{β}^{1}

2 ||∇_{y}_{s}^{k}−_{h}||^{2}_{2}
+h_{v}^{k} |_{h}|,**π**_{4}^{k}i+ ^{β}^{4}

2 ||_{v}^{k} |_{h}|||^{2}_{2}.

(17)

Now, let hiis the i-th pixel of**h**and we discuss two situations when the element hi 6=0,
if h_{i}>0,

hi = (*β*_{1}(∇_{y}_{s})_{i}+ (**π**_{1}^{k})_{i})−(**π**_{4}^{k})_{i}(_{v}^{k})_{i}

*β*_{1}+*β*_{4}(**v**^{k})_{i}(**v**^{k})_{i} ^{,} ^{(18)}
if h_{i}<0,

hi= (−1)−(*β*_{1}(∇_{y}_{s}^{k})_{i}+ (**π**_{1}^{k})_{i})−(**π**_{4}^{k})_{i}(_{v}^{k})_{i}

*β*_{1}+*β*_{4}(_{v}^{k})_{i}(_{v}^{k})_{i} ^{.} ^{(19)}
In summary, the**h**-subproblem has the closed-form solution as follows

**h**^{k+1}=sign(**q**^{k})∗ |_{q}^{k}| −**π**_{4}^{k}_{v}^{k}

*β*_{1}+*β*_{4}**v**^{k}**v**^{k}, (20)

where**q**^{k} =*β*_{1}∇_{y}_{s}^{k}+**π**_{1}^{k}.

156

2) The**z**-subproblem is given as follows

min**z** *µ*||**z**||_{1}+h**s**^{k}−**z**,**π**_{2}^{k}i+ ^{β}^{2}

2 ||**s**^{k}−**z**||^{2}_{2}, (21)
which has the closed-form solution by soft-thresholding strategy [39]

**z**^{k+1}=**Shrink**(**s**^{k}+^{π}

k 2

*β*2, *µ*
*β*2

), (22)

where**Shrink**(**a**, T) =sign(**a**)∗max(|**a**−T|, 0).

157

3) Similar to**z**-subproblem,**w**-subproblem is written as follows
min**w** *λ*k**w**k_{1}+^{β}^{3}

2 ||∇x(**b**−**s**^{k})−**w**+^{π}

3k

*β*_{3}||^{2}_{2}. (23)

The problem (23) has the following closed-form solution by the soft-shrinkage formulation,
**w**^{k+1}=**Shrink**(**q**^{k}, *λ*

*β*_{3}), (24)

where**q**^{k} =∇x(**b**−**s**^{k}) + ^{π}_{β}^{3}^{k}

3.

158

4) The**v**-subproblem can be written as follows

**0**≤**v**≤**1**min h**v**,**c**^{k}i+^{β}^{4}

2 ||**v** |**h**^{k+1}|||^{2}_{2}, (25)

where**c**^{k} =**1**−**π**_{4}^{k} |**h**^{k+1}|. Combining with the constraint**0**≤**v**≤**1**, it has the closed-form solution,
**v**^{k+1}=min(**1**, max(**0**, −**c**^{k}

*β*_{4}|**h**^{k+1}| |**h**^{k+1}|)). (26)

**Algorithm 1:**The algorithm for model (12)

**Input:**The observed image**b**(with stripes), the parameters*λ*,*µ*,*β*_{i},i=1, 2, 3, 4,
the constant*κ*∈(0, ^{1}

*β*1||∇^{T}_{y}||^{2}+*β*2+*β*3||∇^{T}_{x}||^{2}), the maximum number of iterationsM_{iter},
and the calculation accuracytol.

**Output:**The stripes**s**
**Initialize:**

1)k←0,**v**^{0}←**1**,**s**^{0}←**b**,rho←1
**While**rho>tol and k<M_{iter}

2)k←k+1

3) Solve**h**^{k}by Eq. (20)
4) Solve**z**^{k}by Eq. (22)
5) Solve**w**^{k}by Eq. (24)
6) Solve**v**^{k}by Eq. (26)
7) Solve**s**^{k}by Eq. (28)

8) Update the multipliers* π*i,i=1, 2, 3, 4, by Eq. (29)
9) Calculate the error

rho=||∇_{y}**s**^{k+1}−**h**^{k+1}||_{2}+||**s**^{k+1}−**z**^{k+1}||_{2}+||∇_{x}(**b**−**s**^{k+1})−**w**^{k+1}||_{2}+||**v**^{k+1} |**h**^{k+1}|||_{2}.
**Endwhile**

159

5) Here, PADMM based algorithm needs to introduce an extra convex proximal term^{1}_{2}||**s**−**s**^{k}||^{2}_{D},
which is defined as||**x**||^{2}_{D} =**x**^{T}**Dx**, andDis a symmetric positive definite matrix. The**s**-subproblem
becomes a strong convex optimization problem as

min**s** h∇_{y}_{s}−_{h}^{k+1},**π**_{1}^{k}i+ ^{β}^{1}

2 ||∇_{y}_{s}−_{h}^{k+1}||^{2}_{2}
+h**s**−**z**^{k+1},**π**_{2}^{k}i+^{β}^{2}

2 ||**s**−**z**^{k+1}||^{2}_{2}
+h∇x(**b**−**s**)−**w**^{k+1},**π**_{3}^{k}i
+ ^{β}^{3}

2 ||∇_{x}(_{b}−_{s})−_{w}^{k+1}||^{2}_{2}+^{1}

2||_{s}−_{s}^{k}||^{2}_{D},

(27)

where

**D**= ^{1}

*κ***I**−(*β*_{1}∇_{y}^{T}∇_{y}+*β*_{2}+*β*_{3}∇^{T}_{x}∇_{x}),

*κ*∈ 0, 1

*β*_{1}||∇_{y}||^{2}_{2}+*β*_{2}+*β*_{3}||∇_{x}||^{2}_{2}

! . Then, Eq. (27) will be equivalent to:

**s**^{k+1}=argmin

**s**

1

2||**s**−**g**^{k}||^{2}_{2}, (28)

where**g**^{k} =**s**^{k}−*κ*[*β*_{1}(∇y**s**^{k}−**h**^{k+1}) +*β*_{2}(**s**^{k}−**z**^{k+1})−*β*_{3}∇^{T}_{x}(∇x**b**− ∇x**s**^{k}−**w**^{k+1})].

160

6) Finally, we update the Lagrangian multipliers by
**π**_{1}^{k+1}=**π**_{1}^{k}+*β*_{1}(∇_{y}**s**^{k+1}−**h**^{k+1}),
**π**_{2}^{k+1}=**π**_{2}^{k}+*β*_{2}(_{s}^{k+1}−_{z}^{k+1}),

**π**_{3}^{k+1}=**π**_{3}^{k}+*β*_{3}(∇_{x}(**b**−**s**^{k+1})−**w**^{k+1}),
**π**_{4}^{k+1}=**π**_{4}^{k}+*β*_{4}(_{v}^{k+1} |_{h}^{k+1}|).

(29)

Combining steps 1) to 6), we formulate the final algorithm to iteratively solve the proposed`_{0}

161

sparse model (12). In particular, the subproblems all have the closed-form solutions to ensure the

162

accuracy of the algorithm. Finally, the solving process has been summarized in Algorithm 1.

163

In Algorithm 1,*λ*,*µ*,*β*_{1},*β*_{2},*β*_{3},*β*_{4}are some pre-defined parameters,tolandMiterrepresent the

164

positive tolerance value and the maximum iterations, respectively. In this work, we settol=1/255

165

andM_{iter}=10^{3}. In the following, we discuss the convergence of the Algorithm 1.

166

**4. Experiment results**

167

In this section, we compare the proposed method with several state-of-the-art destriping methods,

168

including the wavelet Fourier adaptive filter (WFAF) [3], the statistical linear destriping (SLD) [26],

169

the unidirectional total variation model (UTV) [17], the global sparsity and local variational (GSLV)

170

[24], and the Low-Rank Single-Image Decomposition (LRSID) [25], on both simulated and real remote

171

sensing data. The codes of these methods, except the GSLV method, are available^{4}. As suggested in

172

[25], we utilize the same periodic/nonperiodic stripes function adding stripes intensity [0, 255] to the

173

underlying images. By the similar measure as in [25], the degraded images were normalized between

174

[0, 1]. All experiments are conducted in MATLAB (R2016a) on a desktop with 16Gb RAM and Inter(R)

175

Core(TM) CPU i5-4590: @3.30GHz.

176

To evaluate the effects of different destriping methods, we will compare several qualitative and quantitative assessments. On the qualitative aspect, we show the visual results, the mean cross-track profile and the power spectrum of different methods. We also employ some acknowledged indexes, i.e.,peak signal-to-noise ratio (PSNR)[40], structural similarity index (SSIM) [40] and the relative error (ReErr), to evaluate the performance of different approaches. The ReErr formula is as follows,

ReErr= ||**s**_{added}−**s**_{restored}||_{2}

||**s**added||_{2} ^{,}

where the**s**addedand**s**restoredrepresent the added stripes and restored stripes by different methods,

177

respectively. Then, we will discuss how to select parameters. We note that we test the comparing

178

methods according to the default or suggested parameters in their papers and codes.

179

4.1. Simulated experiments

180

In simulated experiments, the stripes with periodic (Per) and nonperiodic (NonPer) noise are

181

mainly determined by “Intensity” andr. Here, the “Intensity” means the added absolution value of

182

the stripe scope, and therrepresents the stripes ratio level within the remote sensing images. For

183

convenience to compare, different stripes added to remote sensing images will be denoted as a vector

184

with three elements,e.g.,(Per, 10, 0.2) which represents the periodic stripes, the “Intensity” 10 and

185

stripes ratio 0.2.

186

We take six experimental images, which the first, second, third and sixth examples are available

187

on the website^{5}, and the forth and fifth examples are available on the website^{6}, to test the performance

188

of different methods. To compare these methods clearly, we zoom in destriping details on the bottom

189

left or bottom right of the image.

190

**1) Periodic Stripes.**For the periodic stripes case, we only take one example,i.e.,the first column

191

of Fig.2with added stripes (Per, 10, 0.2), to compare the performance. Almost of all existing methods

192

performs quite excellent due to the simple structures of periodic stripes. The first column of Fig.2also

193

demonstrates the consistent conclusion that all comparing approaches remove the periodic stripes and

194

well preserve the image details of stripe-free regions.

195

4 http://www.escience.cn/people/changyi/codes.html.

5 DigitalGlobe withhttp://www.digitalglobe.com/product-samples.

6 MODIS data withhttps://ladsweb.nascom.nasa.gov/

**Table 1.**The ReErr results between**s**_{added}and**s**_{restored}for different methods

images (a) (b) (c) (d) (e) (f)

WFAF 0.1588 0.2828 0.2519 0.2468 0.2386 0.2574
SLD 0.0874 0.1670 0.1723 0.1664 0.1330 0.1346
UTV 0.0831 0.1542 0.2371 0.1818 0.1314 0.1375
GSLV 0.0867 0.1030 0.2385 0.1926 0.0912 0.1654
LRSID 0.0917 0.1884 0.2731 0.2125 0.1450 0.1897
Ours **0.0193** **0.0693** **0.0365** **0.0892** **0.0304** **0.0813**

**2) Nonperiodic Stripes.**For the non periodic stripes case, we test five remote sensing images from

196

the second column to the end column in Fig.2with added stripes (NonPer, 100, 0.6), (NonPer, 50, 0.2),

197

(NonPer, 60, 0.4), (NonPer, 100, 0.4) and (NonPer, 50, 0.6), respectively. Then, we display the destriping

198

results of WFAF, SLD, UTV, GSLV, LRSID and the proposed method for different simulated remote

199

sensing images starting from third row to the end row. See the visual results of the second column, the

200

WFAF method has a obvious black line and changes the intensity contrast of the underlying image

201

significantly. Although the other comparing methods can remove stripes, some regions change the

202

intensity contrast of the underlying image on the left and the right parts, and the proposed method

203

shows a good performance. Then, from the third to sixth examples, we can clearly observe the residual

204

stripes and blurring effects resulted by the others comparing methods. Moreover, our method not only

205

removes stripes completely but also preserves image details well. From Fig.3, we display the smaller

206

patches of Fig.2for visual quality comparisons, and ours results have a better performance than the

207

others.

208

Fig.4shows the estimated stripes based on Fig.2. From Fig.4, we know that the other comparing

209

methods may generate blurring effect and change intensity contrast. Meanwhile, the estimated stripes

210

of the proposed method neither eliminate image structures nor bring in blurring effects for both

211

periodic and nonperiodic stripes cases.

212

In Fig.5, we show the difference/residuals between the added stripes and restored ones. Although

213

ours results have some residuals, the proposed method shows a better performance than the others

214

compared methods. Moreover, we utilize the ReErr results to show the differences/residuals of Fig.5

215

in quantitative aspect. The ReErr results have shown in Table1. From Table1, our results outperform

216

than the other compared methods.

217

**Figure 2.**The visual results of different simulated images. From top to bottom: underlying images,
degraded images, the destriping results of WFAF, SLD, UTV, GSLV, LRSID and Ours. The degraded
images in the second row are respectively added the stripes (from left to right): (Per, 10, 0.2), (NonPer,
100, 0.4), (NonPer, 50,0.2), (NonPer, 60, 0.4), (NonPer, 100, 0.4) and (NonPer, 50, 0.6). Readers are
recommended to zoom in all figures for better visibility.

**Figure 3.**The zoom results of different simulated images in Fig.2. From top to bottom: zoom of the
underlying images, the degrsded images, the destriping results of WFAF, SLD, UTV, GSLV, LRSID and
Ours. Note that the levels of stripes are same as Fig.2.

**Figure 4.**The stripes**s**of different simulated images in Fig.2. From top to bottom: the added stripes
on the underlying image, the extracted stripe components of WFAF, SLD, UTV, GSLV, LRSID and Ours.

Note that the levels of stripes are same as Fig.2.

**Figure 5.**The difference of the added stripes and restored ones. From top to bottom: the difference
results of WFAF, SLD, UTV, GSLV, LRSID and Ours. Note that the levels of stripes are same as Fig.2.

**2) Averagely quantitative performance on 32 test images.** To quantitatively test robustness

218

and effectiveness of the proposed method, Table2and Table3report the averagely quantitative

219

comparisons on 32 remote sensing images, which are randomly selected from three websites^{7}. In the

220

tables, the best PSNR and SSIM results have been identified in bold. Especially, we compare these

221

methods on 32 remote sensing images with fixed parameters for each method.

222

Table2shows the PSNR and SSIM results on periodic stripes with different stripe levels. Although

223

variance of PSNR is not the smallest, the SSIM of the proposed method holds the best performance,

224

and SSIM is an important index to indicate stability on structural similarity of one method. Moreover,

225

our method has the best mean value results of PSNR and SSIM which show the significant advantages

226

than the other comparing methods.

227

7 1) “DigitalGlobe” withhttp://www.digitalglobe.com/product-samples. 2) some subimages of “hyperspectral image of Washington DC Mall” withhttps://engineering.purdue.edu/~biehl/MultiSpec/. 3) “MODIS” data withhttps://ladsweb.

nascom.nasa.gov/

**Table 2.**The mean value of PSNR and SSIM of 32 images with periodic noise

Intensity Intensity=10 Intensity=50 Intensity=100

Ratio r=0.2 r=0.6 r=0.2 r=0.6 r=0.2 r=0.6

WFAF 41.400±3.601 41.702±3.870 37.160±1.975 37.553±1.975 32.196±1.457 32.501±1.732
SLD 42.037±2.927 41.048±2.909 41.710±2.930 41.957±2.928 40.614±2.549 41.644±2.836
PSNR UTV 42.030±3.229 41.032±2.886 40.920±2.773 43.086±2.298 41.470±3.385 41.058±3.299
GSLV 42.552±2.955 42.630±2.886 42.202±3.058 43.533±2.856 43.431±3.091 43.801±2.705
LRSID 43.948±2.104 42.775±2.010 42.308±2.169 44.548±1.976 43.779±2.500 44.035±2.014
Ours **52.918**±**4.074** **49.497**±**3.956** **52.853**±**4.910** **49.212**±**4.390** **52.854**±**4.902** **49.182**±**4.368**
WFAF 0.9934±0.0058 0.9936±0.0062 0.9887±0.0084 0.9905±0.0078 0.9818±0.0103 0.9847±0.0085

SLD 0.9966±0.0029 0.9966±0.0029 0.9965±0.0031 0.9965±0.0032 0.9962±0.0033 0.9964±0.0037
SSIM UTV 0.9959±0.0027 0.9959±0.0027 0.9911±0.0025 0.9928±0.0023 0.9954±0.0024 0.9937±0.0076
GSLV 0.9991±0.0077 0.9968±0.0076 0.9916±0.0079 0.9903±0.0082 0.9966±0.0085 0.9969±0.0053
LRSID 0.9990±0.0107 0.9945±0.0056 0.9932±0.0044 0.9947±0.0032 0.9936±0.0047 0.9957±0.0031
Ours **0.9994**±**0.0007** **0.9987**±**0.0011** **0.9994**±**0.0013** **0.9986**±**0.0016** **0.9994**±**0.0062** **0.9986**±**0.0019**

**Table 3.**The mean value of PSNR and SSIM of 32 images with nonperiodic noise

Intensity Intensity=10 Intensity=50 Intensity=100

Ratio r=0.2 r=0.6 r=0.2 r=0.6 r=0.2 r=0.6

WFAF 40.971±2.523 39.372±2.249 30.536±1.508 37.609±2.263 24.849±1.573 22.594±1.541
SLD 41.476±2.592 40.935±2.201 35.964±1.510 42.007±3.020 30.963±1.414 28.403±1.729
PSNR UTV 41.153±2.880 38.615±2.041 35.648±1.527 42.505±3.010 31.055±4.687 31.599±2.578
GSLV 42.282±2.359 39.018±1.654 41.985±1.239 39.838±2.903 36.184±1.399 35.408±2.472
LRSID 42.672±1.418 39.034±1.302 42.814±1.349 40.497±2.024 37.779±1.212 33.559±1.132
Ours **48.801**±**3.985** **44.700**±**3.784** **49.057**±**4.791** **49.057**±**4.492** **44.365**±**5.106** **39.452**±**4.494**
WFAF 0.9925±0.0056 0.9903±0.0069 0.9744±0.0104 0.9905±0.0081 0.9364±0.0207 0.9029±0.0565

SLD 0.9965±0.0031 0.9952±0.0031 0.9950±0.0041 0.9964±0.0032 0.9907±0.0060 0.9823±0.0142
SSIM UTV 0.9958±0.0029 0.9934±0.0052 0.9937±0.0042 0.9914±0.0056 0.9886±0.0193 0.9851±0.0122
GSLV 0.9982±0.0016 0.9917±0.0042 0.9962±0.0101 0.9967±0.0088 0.9956±0.0091 0.9933±0.0152
LRSID 0.9983±0.0032 0.9934±0.0113 0.9891±0.0070 0.9962±0.0042 0.9975±0.0091 0.9924±0.0402
Ours **0.9991**±**0.0006** **0.9956**±**0.0035** **0.9990**±**0.0010** **0.9986**±**0.0016** **0.9979**±**0.0012** **0.9942**±**0.0042**

For the nonperiodic stripes, we show the mean value results in Table3. The WFAF method shows

228

the instability, and the PSNR and SSIM of LRSID method are consistent with the results in [25]. From

229

the two tables, our method always shows a good performance significantly.

230

In Fig. 6, we take two examples of Table2to show the PSNR and SSIM performance of all

231

comparing methods on each image. The y-axis stands for the value of PSNR or SSIM and the x-axis

232

represents thei-th image of 32 examples. Fig.6(I) and Fig.6(II) are the PSNR and SSIM performance

233

of stripes (Per, 100, 0.6), and Fig.6(III) and Fig.6(IV) are that of stripes (NonPer, 50, 0.2), respectively.

234

Although the PSNR results fluctuate with respect to different images, our method holds the best

235

PSNR results on almost of all images. Moreover, the SSIM results show the best performance with the

236

smallest variance, which is consistent with the results of Table2and Table3. From Fig.6, our method

237

is superior to the other comparing methods.

238

0 5 10 15 20 25 30 35

i-th image

25 30 35 40 45 50 55 60

PSNR

WFAF SLD UTV GSLV LRSID Ours

0 5 10 15 20 25 30 35

i-th image

0.94 0.95 0.96 0.97 0.98 0.99 1

SSIM

WFAF SLD UTV GSLV LRSID Ours

0 5 10 15 20 25 30 35

i-th image

25 30 35 40 45 50 55 60

PSNR

WFAF SLD UTV GSLV LRSID Ours

0 5 10 15 20 25 30 35

i-th image

0.94 0.95 0.96 0.97 0.98 0.99 1

SSIM

WFAF SLD UTV GSLV LRSID Ours

(I) (II) (III) (IV)

**Figure 6.**The PSNR and SSIM performance on 32 images for the stripes (Per, 100, 0.6) and (NonPer,
50, 0.2). The x-axis represents each image and the quantitative results are shown in y-axis. (I) and (II)
are the PSNR and SSIM results for the stripes (Per, 100, 0.6), respectively. (III) and (IV) respectively
represent the PSNR and SSIM performance of the stripes (NonPer, 50, 0.2).

4.2. Real experiments

239

We also display the destriping results of six methods for six real remote sensing images, which

240

are also available on the website^{8}, see Fig.7. Similar to Fig.2, the six real images with different stripes

241

are shown in the first row, and the destriping results of all comparing methods are presented from the

242

second row to the end row.

243

In Fig.7, for the first, fifth and last real images, the proposed method not only removes the stripes

244

completely, but also preserves image details on stripe-free regions well. Note that the methods GSLV

245

and LRSID fail to obtain excellent results for the first image as the mentioned in their papers. For the

246

forth column, there are also several stripe residuals with WFAF and SLD, and the wide black shadow

247

areas appear by the UTV, GSLV and LRSID methods. Moreover, the destriping results of the WFAF

248

and SLD leave significant stripes for the second image, and still exist the wispy stripes for the third

249

example. According to several real experiments, the results demonstrate the universal effectiveness

250

and stability of the proposed method.

251

4.3. More discussion

252

**1) Qualitative Analysis.** For the further comparisons of different destriping methods for

253

simulated and real remote sensing images, we show the following two assessments. One is the

254

mean cross-track profile that the x-axis stands for the column number of an image and the y-axis

255

represents the mean value of each column, see Fig.8and Fig.10. The other is the power spectrum that

256

the x-axis is the normalized frequencies of an image, and the y-axis shows the spectral magnitude with

257

a logarithmic scale, see Fig.9and Fig.11.

258

In simulated experiments, the mean cross-track profile of the first image of Fig.2has been shown

259

in Fig.8. Note that Fig.8(a) shows the mean cross-track profile of the underlying image, and Fig.8(b)

260

is the result of the degraded image. Moreover, Fig.8(c)-(f) are the mean cross-track profile results of

261

the six destriping methods, respectively. From the overall perspective, Fig.8(d) and Fig.8(e) have

262

obvious change of the intensity contrast. Seeing the details, Fig.8(c)-(g) have some mild fluctuations

263

which are different with the underlying image in Fig. 8(a). The proposed method shows the best

264

performance, since it is almost same as the original one.

265

In addition, the power spectrum results of the second image of Fig.2has been shown in Fig.9.

266

We denote the power spectrum results as Fig.9(a)-(h) which represent the power spectrum results of

267

the underlying image, the degraded image and the destriping results of six methods, respectively. Fig.

268

9(c)-(g) have more fluctuations which indicate these methods may have the stripe residuals or bring a

269

8 https://ladsweb.nascom.nasa.gov/

**Figure 7.**The visual results of different real images. From top to bottom: the real images, the destriping
results of WFAF, SLD, UTV, GSLV, LRSID and Ours. Readers are recommended to zoom in all figures
for better visibility.

little new noise in their destriping processes. For our method,i.e.,Fig. 9(f), it not only removes all

270

stripes, but also preserves almost the essential details such as edges.

271

0 50 100 150 200 250 300 350 400

Line number

0.35 0.4 0.45 0.5

Mean value

0 50 100 150 200 250 300 350 400

Line number

0.2 0.3 0.4 0.5 0.6 0.7

Mean value

0 50 100 150 200 250 300 350 400

Line number

0.35 0.4 0.45 0.5

Mean value

0 50 100 150 200 250 300 350 400

Line number

0.35 0.4 0.45 0.5

Mean value

(a) (b) (c) (d)

0 50 100 150 200 250 300 350 400

Line number

0.35 0.4 0.45 0.5

Mean value

0 50 100 150 200 250 300 350 400

Line number

0.35 0.4 0.45 0.5

Mean value

0 50 100 150 200 250 300 350 400

Line number

0.35 0.4 0.45 0.5

Mean value

0 50 100 150 200 250 300 350 400

Line number

0.35 0.4 0.45 0.5

Mean value

(e) (f) (g) (h)

**Figure 8.**Spatial mean cross-track profiles for simulated image of the first simulated example of Fig.2.

(a) Underlying image. (b) Degraded image. Destriping results by (c) WFAF, (d) SLD, (e) UTV, (f) GSLV, (g) LRSID, (h) Ours.

0 0.1 0.2 0.3 0.4

Normalized frequency

0 5 10 15

Power spectrum

0 0.1 0.2 0.3 0.4

Normalized frequency

0 5 10 15

Power spectrum

0 0.1 0.2 0.3 0.4

Normalized frequency

0 5 10 15

Power spectrum

0 0.1 0.2 0.3 0.4

Normalized frequency

0 5 10 15

Power spectrum

(a) (b) (c) (d)

0 0.1 0.2 0.3 0.4

Normalized frequency

0 5 10 15

Power spectrum

0 0.1 0.2 0.3 0.4

Normalized frequency

0 5 10 15

Power spectrum

0 0.1 0.2 0.3 0.4

Normalized frequency

0 5 10 15

Power spectrum

0 0.1 0.2 0.3 0.4 0.5

Normalized frequency

0 5 10 15

Power spectrum

(e) (f) (g) (h)

**Figure 9.**Power spectrum for simulated image of the second example of Fig.7. (a) Underlying image.

(b) Degraded image. Destriping results by (c) WFAF, (d)SLD, (e) UTV, (f) GSLV, (g) LRSID, (h) Ours.

In real experiments, we also show the mean cross-track profile and the power spectrum in Fig.

272

10and Fig.11, respectively. Fig.10shows the mean cross-track profile results of the first column of

273

Fig.7. Note that Fig.10(a) is the mean cross-track profile result of the first real remote sensing image,

274

and Fig. 10(b)-(g) show the profile results of six destriping methods, respectively. In general, the

275

profiles of the destriping method should smoothen huge fluctuates and maintain primary structure

276

information. However, the profiles of WFAF and LRSID have obvious fluctuations where the stripes

277

still exist, and that of SLD is over-smooth missing a lot of underlying image details. In Fig.10(d) and

278

Fig. 10(e), although stripes are mostly removed, the destriping profiles have some mild burrs and

279

too much smoothness because of the unidirectional property of UTV and the global sparsity of GSLV,

280

respectively. In addition, the profile of the proposed method,i.e.,Fig.10(g), can realize the desired

281

result both on removing stripes and keeping underlying image details.

282

In Fig.11, the power spectrum results of the forth example of Fig. 7are plotted. Fig. 11(a)-(h)

283

represent the power spectrum results of the forth real remote sensing image and six destriping methods,

284

respectively. We observe that the real remote sensing image in Fig.11(a) has much fluctuates where

285

stand for stripes. According to the power spectrum results of the six methods in Fig.11(b)-(f), although

286

the stripes are almost removed well, there are still some slight blurring regions, while the proposed

287

method shows the best performance in Fig.11(g).

288

0 50 100 150 200 250 300 350 400 450 500

Line number

0.1 0.2 0.3 0.4 0.5

Mean value

0 50 100 150 200 250 300 350 400 450 500

Line number

0.1 0.2 0.3 0.4 0.5

Mean value

0 50 100 150 200 250 300 350 400 450 500

Line number

0.1 0.2 0.3 0.4 0.5

Mean value

0 50 100 150 200 250 300 350 400 450 500

Line number

0.1 0.2 0.3 0.4 0.5

Mean value

(a) (b) (c) (d)

0 50 100 150 200 250 300 350 400 450 500

Line number

0.1 0.2 0.3 0.4 0.5

Mean value

0 50 100 150 200 250 300 350 400 450 500

Line number

0.1 0.2 0.3 0.4 0.5

Mean value

0 50 100 150 200 250 300 350 400 450 500

Line number

0.1 0.2 0.3 0.4 0.5

Mean value

(e) (f) (g)

**Figure 10.** Spatial mean cross-track profiles for the first real example of Fig. 7. (a) Real image.

Destriping results by (b) WFAF, (c)SLD, (d) UTV, (e) GSLV, (f) LRSID, (g) Ours.

**2) The influence of different regularization terms in the proposed model.** Fully considering

289

the destriping problem (2) and the optimization model (12), we assume that R_{2}is a necessary term,

290

since R2is the only term to describe the property of the underlying image**u**. To confirm whether both

291

R1and R2are necessary priors as well as have significant contribution for destriping performance,

292

in Fig.12, we give the mean value of PSNR and SSIM for 32 images as before. Here, R_{12}represents

293

R_{1}+R_{2}, R_{23} stands for R_{2}+R_{3}and R_{123}represents R_{1}+R_{2}+R_{3}(i.e., the proposed model). Please

294

find the definitions of R1, R2, R3from Eq. (9), Eq. (10) and Eq. (11), respectively.

295

Fig.12(I) and Fig.12(II) show the mean value of PSNR and the mean value of SSIM on 32 images

296

same as before for periodic stripes. The periodic stripe levels (a)-(f) are (Per, 10, 0.2), (Per, 10, 0.6), (Per,

297

50, 0.2), (Per, 50, 0.6), (Per, 100, 0.2) and (Per, 100, 0.6), respectively. Moreover, Fig.12(III) and Fig.12

298

(IV) display the mean value of PSNR and the mean value of SSIM on 32 images for nonperiodic stripes.

299

The nonperiodic stripe levels (a)-(f) stand for (NonPer, 10, 0.2), (NonPer, 10, 0.6), (NonPer, 50, 0.2),

300

(NonPer, 50, 0.6), (NonPer, 100, 0.2) and (NonPer, 100, 0.6), respectively.

301

From the results in Fig.12, we can conclude three points. 1) The results both PSNR and SSIM of

302

the proposed model (i.e.,R123) perform the best than those of the other two models. 2) For R12and R23,

303

R23shows more stability than R12as the green bars do not significantly change with different stripes.

304