Removing Striping Noise from Cloudy Level 2 Sea Surface Temperature and Ocean Color Datasets

Abstract: This paper introduces a new destriping algorithm for remote sensing data. The method is based on combined Haar Stationary Wavelet transform and Fourier filtering. State-of-the-Art methods based on the discrete wavelet transform (DWT) may not always be effective and may cause different artifacts. Our contribution is three-fold: i) we propose to use the Undecimated Wavelet transform (UWT) to avoid as much as possible shortcomings of the classical DWT; ii) we combine a spectral filtering and UWT using the simplest possible wavelet, the Haar basis, for a computational efficiency; iii) we handle 2D fields with missing data, as commonly observed in ocean remote sensing data due to atmospheric conditions (e.g., cloud contamination). The performances of the proposed filter are tested and validated on the suppression of horizontal strip artifacts in cloudy L2 Sea Surface Temperature (SST) and ocean color snapshots.


Introduction
Passive sensors on board remote sensing platforms use several scanning methods to generate land and sea surface imagery.Images provided by the Sun-synchronous Earth orbit satellites are achieved by a combination of progressive scanning lines in the cross-track direction while the sensor platform motion is along the in-track direction.The resulting images are often contaminated by several types of noise.These undesired artifacts have a direct impact on the visual quality of the provided images.If not corrected, these noises will then further cause processing errors.In this paper, we will deal with striping noise patterns.This type of noise are often present in sea surface temperature and ocean color images provided by infrared and optical imaging spectrometers (e.g.MODIS, VIIRS...).It consists in sharp repetitive patterns which take the form of stripes over the entire image [1] (See Fig. 1).
The reduction of these stripe artifacts is an important research topic.A large number of destriping algorithms have been recently suggested.All the scene-based methods of the destriping literature exploits geometrical priors about the noise.These priors are related to the regular periodicity of the noise.One may cite a variety of approaches based on low-pass filtering implemented in the spatial or frequency domain [2][3][4][5][6].A common feature shared by these methods is that they give rise to blur artifacts.More sophisticated filtering approaches have been proposed.Multiresolution analyses using wavelet decompositions were investigated in [7,8].More recently, variational methods were introduced and explored in [9,10].These methods may however be prohibitively expensive for large datasets.
Reducing striping artifacts in an effective manner without blurring the images still remains challenging.Moreover, the considered case-study applications, infrared sea surface temperature and ocean color observations [11], do not involve cloud-penetrating sensors, What may result in a very high rate of missing data (gaps) in the provided images.The traditional Fourier or wavelet filtering techniques cannot handle such gaps in images and require modifications.In this paper, we address the removal of striping noise in ocean remote sensing images, possibly involving missing data as illustrated in Fig. 1.We develop a destriping algorithm based on a combined wavelet-Fourier filtering.Our algorithm can be regarded as an extension of [7].We evaluate our contribution for real ocean satellite-derived images with a focus on both SST and ocean color imagery.
This paper is organized as follows.In section 2, we provide a short review of the assumptions required by the wavelet decomposition and Fourier transform.In section 3, a detailed technical description of the proposed destriping algorithm is given.We report numerical experiments with real ocean remote sensing data, including a comparison to state-of-the-art approaches in Section 4. Section 5 concludes this paper., (e) and (f) using respectively our method, [7] with Haar wavelet, [7] with Daubechies-4 and using a variationnal approach proposed by [10].

Problem statement
Let us consider an observed image u sn (i) defined in a rectangular domain i ∈ Ω, affected by an additive stripe-type noise.The image degradation model can be expressed by the following equation where u(i) would be the true value at pixel i and sn(i) is the striping noise perturbation.The analysis of satellite images shows that the striping noise can be considered as a structured noise in which the large variability is along the y axis of the image, as illustrated in Fig. 1.By exploiting this prior on the spatial structure of the undesired noise, the filtering problem consists in removing the striping noise of the images without introducing any blurring effects.Following [7], the proposed approach relies on an appropriate decomposition of the image u sn (i) so that the striping noise effect can be isolated from the original hidden image.Notice that we will not deal with other stationary noise, which may be present in the images and removed using appropriate methods.Several filtering approaches have been developed for the removal of striping noise in satellite images.Following the idea that striping noise is a superposition of quasi-perfect periodic signals and can be easily identified in the 2D Fourier spectrum, one can construct a filter for removing it at a given frequency in the Fourier domain [6].The weak point of this method is the fact that this filter does not only remove part of the spatial frequency component related to the undesired stripe noise, but also eliminates and reduces the part of the same component present in the real (physical) signal.In order to avoid this over-denoising effect, [7] proposes to perform this spectral filtering method after a first finite-level discrete dyadic wavelet transform.In this algorithm only the wavelet coefficients (details) are assumed to contain the undesired striping noise and are filtered in the Fourier domain.All the approximation coefficients are kept and the resulting image is obtained by the inverse wavelet transform using the denoised coefficients.The decimated Wavelet analysis (DWT) takes advantage of scaling and directional properties to detect and remove striping patterns in the wavelet domain.The DWT is a non-redundant decomposition [12].This is particularly interesting for storage and computational efficiency purposes.Nevertheless, for reconstruction-related applications, which is our use case here, the DWT does not fulfill the translation-invariance property, what may lead to a large number of artifacts when modifying its wavelet coefficients.

Proposed destriping approach: THE UWT-Fourier based destriping scheme
Following [7], we propose to tackle the problem of removing striping noise through a combined wavelet-Fourier approach.As previously mentioned, destriping with the traditional (orthogonal) discrete wavelet transform sometimes exhibits visual artifacts.These artifacts are caused by the sensibility of these algorithms to translation.The Undecimated wavelet transform (also called stationary wavelet transform) was designed to overcome the lack of translation-invariance of the DWT [12].This property is achieved by removing the decimation step in the orthogonal wavelet transform.
Haar-based UWT decomposition: In the proposed Desptriping scheme, the 2D Haar wavelet transform is the proposed analysis technique.The Haar basis function is well known as the first and the simplest wavelet analysis.The associated scaling and wavelet functions (denoted respectively by φ(x) and ρ(x)) are illustrated in Fig. 2. The major advantages with the use of the Haar analysis are the following: 1. Interpretability: the form of Haar filter is simple and easy to implement; 2. Computational efficiency: unlike the continuous wavelets, fast calculations are obtained, which is important for large satellite derived data products; 3. The inverse transform is performed without any edge effect artifacts.This a key feature in our case as we deal with images involving missing data.
Notice that for applications where reconstruction is needed, the Haar transform also has limitations.
Images reconstructed with the Haar filter may exhibit block-like artifacts when the decimation is involved.The considered UWT approach resolve this problem.The original image u sn is represented in the UWT domain by a sequence of details at different scales and orientations along with an approximation image at a predefined coarsest scale.
where U J 0 represents the approximation image at the lowest scale J and The original image can be obtained using its coefficients by the inverse UWT.

Fourier filtering:
We assume that the noise is periodic, invariant along the x-axis and distributed over several spectral component.Given the Haar UWT decomposition of the noisy image u sn , we further assume that the striping noise is only present in the vertical and diagonal bands of the UWT decomposition and perform a Fourier filtering independently for each detail image of these two where F and F −1 refer, respectively, to the Fourier and the inverse Fourier transform.The denoised image ũ then resorts to where W −1 is the inverse UWT transform, j and k are, respectively, the scale level and the orientation of the wavelet subband.

Fourier filter design:
In the 2D wavenumber domain k = (k x , k y ) the ideal horizontal (resp.vertical) stripes are almost located near the high frequency part in the vertical (resp.horizontal) direction, i.e (0, k y ) [see Fig. 3 for an illustration].The destriping Fourier filter is designed to remove this wavelengths from the Fourier transform of the noisy UWT coefficients.For this purpose, we apply a band-pass filter around k y = 0.This can be achieved by the pointwise multiplication of the FFT coefficients with an inverted Gaussian function where σ controls the width of the filter in the k x -direction.Fig. 3 shows an example of such a function.Since the observed striping artifacts are almost horizontal (resp.vertical), the value of σ must be small so that the Fourier coefficients are set to zero only at k x = 0 (resp.at k y = 0).Thereby, the filtering process is expected to eliminate striping artifacts without producing blur effects.We define the method noise of u as the image difference This method noise (or noise residue) should be as similar as possible to an image composed only of striping patterns.
Destriping in presence of missing data: Destriping in the presence of missing data is a very challenging task, especially when considering a Fourier and wavelet analysis.These decompositions cannot handle images with missing data and require the images to be interpolated prior to the computation of the decompositions.A classical zero-padding strategy may result in a poor estimation of the decomposition and may produce severe visual artifacts due to the introduced discontinuities.It may be noted that we do not address the joint removal of striping noise and interpolation of missing data areas.The goal is rather to apply as preprocessing step an appropriate interpolation scheme, which will result in no noticeable discontinuities of the denoised image at the boundaries of missing data areas.For this purpose, we consider the harmonic image inpainting as described in [13].The method smoothly interpolates inward from the pixel values on the outer boundary of the missing regions.In the following, we will briefly explain the method.Let us consider an image denoted by f and defined in a rectangular domain denoted by Ω. Suppose that this image is only known at a subset Ω k ⊂ Ω.The harmonic inpainting method consists in filling in the missing region by solving the following Dirichlet boundary value problem where ∂ n denotes the derivative operator normal to the boundaries.We can also consider higher-order differential operators as interpolant (e.g., the biharmonic operator ∆ 2 ).

Experimental Results
Experimental setting: Several BT, SST and ocean color snapshots acquired by MODIS Aqua/Terra [14], VIIRS NPP [15] and TIRS Landsat8 [16] were selected to illustrate and validate the performance of the proposed algorithm.We choose scenes heavily affected with striping noise.
Examples of destriping The visual improvement of Modis-Terra ocean color snapshots is illustrated in Fig. 4 and Fig. 5. Smaller images (≈ 400 × 400 pixels) compared to the entire received granule was selected.The reason is that the striping effect is visually hard to observe on significantly large images.Visually, the proposed preprocesing result in a clear improvement of the visual quality of the snapshots.They do not contain any residual horizontal stripes and do not reveal any blur  effect.Fronts and sharp gradients areas are degraded by stripe noise in the original images.This occurs because stripes lead to larger vertical gradients.These geometrical features are significantly enhanced in the resulting destriped snapshots.In panels Fig. 4(f) and Fig. 5(f), we plotted the averaged Fourier power spectrum.The stripes are revealed in the Fourier spectrum by several impulses (or pics) located at different wavelengths, often at Mid and High frequency components.We can observe from the analysis of the the spectral densities before and after the destriping process that striping noise components are no more observed in the spectra of the processed images.
We also applied our algorithm to SST snapshots provided by Modis sensor onboad Aqua platform.As shown in Fig. 6, we reported similar to those obtained for ocean colour snapshots.transform do have an effect on the quality of the obtained images.Fig. 9 shows the joint influence of these key parameters.It suggests that the typical variance parameter must be small.Large values produce a blurring effect.In the various results illustrated in this paper, the default value for the σ parameter is set to 5. Notice that even with smaller values the algorithm gives good performances.Regarding the UWT, a suitable number of decomposition shall not exceed 5 levels.
To deal with the missing data issue while using the Fourier filtering, our method uses a preprossing step which consists in filling in the missing areas using the Laplacian inpainting method.Fig. 8 stresses the benefit of such an interpolant comparing with the classical zero-padding method.This method has been chosen for a practical considerations, since it is parameter-free and inexpensive in computer storage space (relatively sparse matrices to invert) compared to inpainting method based on high-order differential operators (e.g., biharomic inpainting).This method is especially suitable for images involving high rate of missing data where the discretisation of the Laplcian operator gives rise to large matrices.From the reported experiment here using a bi-harmonic interpolant, we expect other inpainting methods based in diffential operator (e.g., AMLE [13]) to lead to similar destriping performances, at the expense of an increased computational complexity.

Comparison to state-of-the-art algorithms:
We performed a comparisons to two recent state-of-the-art algorithms, namely the Fourier-wavelet scheme proposed by [7] and the non-local variational approach proposed by [10].
Compared to our approach, [7] uses the decimated wavelet transform (DWT) which is less effective, especially in the heavy stripped image.In addition, the use of the Haar filter is not suitable and the number of taps (i.e., nonzero coefficients) for the chosen wavelet form must be large.This is illustrated by the results reported in Fig. 9.The quality of the image generated by our method is superior to the images obtained using [7].Unlike panels Fig. 9(d)-(e), all the stripes was removed and the image contains no artifacts.
Fig. 9 also shows the improved destriping performance of our algorithm compared to [10].We further compare in Fig. 10 the method noise resulting from our algorithm and [10].The visual inspection of the associated Fourier power spectra shown in panels Fig. 10(c)-(d) suggests that the energy of our method noise is almost distributed in the narrow horizontal wavenumbers band near k y = 0.By contrast, the energy related to the method noise of the destriped image using [10] is distributed over a broad spectral band and does not conform to the prior assumption about the geometrical nature of the undesired noise.A direct impact of this observation can be seen in the averaged Fourier spectrum reported in panel Fig. 10(e).In fact, The signal spectral magnitude is attenuated for frequencies located near the mid-wavenumber region.
To achieve a quantitative comparison with the considered sate-of-the-art methods for the ocean color maps we compute the mean of the cross-track profiles of each image.This quantitative metric measurement gives a good indication of the strength of the striping noise present in a given image.It consists in calculating the average value along each scan line.The presence of stripes translates to the mean cross-track profile by a strong periodicity.Using this metric, the goodness of a given destriping algorithm is described in terms of its ability to remove this periodic component and obtain a smooth signal.We plot in Fig. 11 the cross-track profiles of images displayed in Fig. 4 and Fig. 5.As expected, the original destriped image exhibits a strong periodicity.The profiles of the destriped images using our approach as well as [7] and [10] are also reported.While our approach and [7] resort to slow-varying profiles with no periodic pattern, [10] does not sufficiently remove these periodic components which are still visible at the end of the scan.
Impact of the destriping for the analysis of SST and Ocean Color snapshots: We illustrate the potential of the proposed destriping algorithm for an operational use of the resulting sea surface fields.The first application deals with the application of our algorithm as a preprocessing step in the automatic detection of SST fronts [17].For this purpose, we apply a Sobel filter to a subset of an original Modis SST maps at full sensor resolution of ≈ 1km and using a downsampled resolution of ≈ 4km (Notice that in [17] the edge detection algorithm is performed using data at ≈ 4km).We    perform the same operation using the destriped version of the considered image.Fig. 12 illustrates the results obtained for this experiment.As expected, in both the full and downsampled resolution of the real image, strong gradients are exhibited for real SST fronts and also caused by the vertical stripes.As such, an automatic detection algorithm would hardly be able to discriminate these two classes of gradient patterns.By contrast, the effect of the destriping is clear in the processed images.It reveals more clearly the geometry of the SST fronts, which were occluded in some cases by the striping artifacts.This example illustrates that our destriping scheme can enhance the detection of thermal fronts using simple automatic front detection algorithms.The retrieval of Level-2 SST products from a Level-1 brightness temperature data may be considered as a second potential use of the proposed destriping approach.The linear [18] (resp.nonlinear [19]) SST retrieval algorithms are typically based on linear (resp.nonlinear) combination of brightness temperature extracted from several channels.Brightness temperature datasets are also involved with stripe artifacts.Thus, it is more convenient to reduce stripe noise before performing the retrieval.Here we report such an application using VIIRS data.Brightness temperatures from the 375 m resolution Imagery Bands (I-Band) are used with 750 m resolution SST fields obtained from the VIIRS Moderate Resolution Bands (M-Band) to obtain 375 m SST fields.The algorithm consists in computing regression coefficients by a rolling-window analysis.We report two examples in Fig. 13.They illustrate the benefits of the destriping prior to the application of SST retrieval algorithm.

Conclusion
As developed, an efficient methodology can be applied to help removing striping noise artifacts in ocean remote sensing images.Such artifacts are common due to the scanning process underlying the formation of satellite-derived sea surface observations.As proposed and tested, the novel scheme  combines a UWT decomposition of the image to a Fourier filter.Contrary to most state-of-the-art techniques, this proposed scheme can very efficiently deals with missing data.On different real satellite-derived images, we demonstrated the relevance of the proposed approach compared to previous work.The use of the UWT is regarded as a key component, which brings clear benefit compared to the DWT [7] and variational prior [10].As illustrated, we anticipate the impact of the proposed destriping scheme for the further analysis of satellite-derived sea surface fields, especially to analyze signatures of sub-mesoscale upper ocean dynamics.

Figure 1 .
Figure 1.Destriping of an ocean brightness temperature snapshot obtained by TIRS Lanndsat 8 on September 11, 2014.The original data are shown in panel (b) while the destriped data are shown in panels (b), (d), (e) and (f) using respectively our method,[7] with Haar wavelet,[7] with Daubechies-4 and using a variationnal approach proposed by[10].

Figure 3 .
Figure 3. (a) A simulated perfect vertical Gaussian striped sheet.(b) The associated 2D power spectrum.All non-zero values are located near the high frequency part in the horizontal direction (c) The inverted Gaussian function considered as the Fourier filter in our algorithm

Figure 4 .
Figure 4. Illustration of the destriping of ocean color products: (a) chlorophyll-a concentration obtained by Modis Terra on December 22, 2015 around 06:20 UCT in the Arabian sea region.The destriped data are shown in panel (b).The gradient magnitude of the original image and the destriped image are shown respectively in panel (d) and (e).Comparison of the averaged Fourier power spectrum is illustrated in panel (f).

Figure 5 .
Figure 5. Illustration of the Destriping of ocean color products: (a) chlorophyll-a concentration obtained by Modis Terra on December 22, 2015 around 06:20 UCT in the Arabian sea region.The destriped data are shown in panel (b).The gradient magnitude of the original image and the destriped image are shown respectively in panel (d) and (e).Comparison of the averaged Fourier power spectrum is illustrated in panel (f).

Figure 6 .
Figure 6.Illustration of the destriped results for SST snapshots derived from Aqua Modis sensor.(a)Original images (b) Destriped images with the proposed algorithm.(c) Gradient magnitude of the original fields.(d) Gradient magnitude of the Destriped images

Figure 7 .Figure 8 .
Figure 7. Destriping results of the image shown in Fig.5.1(a) for different values of the variance related to the Fourier filter g α (Eq.3) and for different wavelet decomposition levels.

Figure 11 .
Figure 11.Mean cross-track profiles of the original and destriped images shown in Fig.5 compared with the state-of-the-art algorithms.

Figure 12 .
Figure 12.Illustration of the impact of the striping noise in the detection of thermal fronts from SST snapshots.The magnitude gradient field of a Modis snapshot computed using the Sobel operator using (a) the full image resolution (b) downsampled version at 4km (c) destriped image at the full sensor resolution (d) the destriped image at 4km resolution.

PreprintsFigure 13 .
Figure 13.Illustration of the potential use of the proposed destriping algorithm.The destriped VIIRS SST data are used along with the I-Band BT at 375 m to produce SST product at 375 m.

Figure 14 .
Figure 14.Power spectral density (PSD) plots for SST snapshots in Fig 5.6.The frequency components related to the periodic stripes in VIIRS 750 m (Black curves) are extremely attenuated.The graphs have been artificially translated for sake of clarity .