Sensibility to noise of new multifractal fusion methods for ocean variables

The repeated observation of the same signatures of mesoscale and submesoscale features in different ocean variables indicates that some common, non-linear processes affect them to a significant extent. A new method to exploit these common signatures to improve the quality of a noisy variable (i.e. increasing the signal-to-noise ratio) using another variable as template has recently been introduced. The method is based on superimposing the multifractal structure of singularity exponents from the template variable to the variable to be enhanced. In this paper, we will discuss the sensitivity of this method to the presence of noise of different types and amplitude. Our results indicate that multifractal methods can be a key to enhancing the existing databases of remote sensing images and give hints about non-linear dynamics of the ocean.


Introduction
The introduction of new remote sensing platforms for the observation of oceans has posed more questions than the answers they have provided.Previously unobserved processes can now be seen by repeated observations (as for instance with the SMOS -Soil Moisture and Ocean Salinity -mission (Font et al., 2012) and the large-scale deployment of the Argo floater system) and improved numerical models.Making sense of the large amount of new data and at the same time understanding the associated processes is a pressing challenge.
On the other hand, ocean structures associated with many different processes (frontogenesis and eddy generation Hallberg andGnanadesikan, 2006, eddy propagation Isern-Fontanet et al., 2003;Chelton et al., 2007, variations in sea roughness, upwelling Mason et al., 2011, etc.) have been shown to leave a recognisable footprint in many different variables.This common footprint is the result of similar nonlinear terms within the equations defining the evolution of ocean scalars (such as sea surface temperature -SST, sea surface salinity -SSS, chlorophyll concentration -CC, etc.), even if those non-linear terms are not the dominant contribution to the dynamics of the particular scalars.This common signature in ocean scalars indicates that a part of the information conveyed by the different ocean variables is in fact redundant.Understanding this redundancy is useful in two different ways.On one hand it can serve to assess and even to measure the magnitude of the associated terms in the evolution equations.On the other hand, that redundancy can be exploited to code information in a more compact way and even to enhance data quality when one of the data sources is corrupted or otherwise contaminated.
The introduction of the microcanonical multifractal formalism (Isern-Fontanet et al., 2007;Turiel et al., 2008b) provides a unified formalism to model and quantify this redundancy in a consistent mathematical model.The common signature in all ocean variables is the multifractal structure underlying all of them, which is dimensionless and independent of the signal amplitude of any particular scalar.Evidence is now extensive about the validity of this approach to deal with ocean variables (Turiel et al., 2005(Turiel et al., , 2008a(Turiel et al., , 2009)).Recently, a non-parametric, non-linear method to fuse the information coming from different scalars using their common multifractal structure has been proposed (Umbert et al., 2014) for dealing with a specific application (improving the quality of SMOS maps).
The goal of this paper is to analyse the performance of this novel approach to fuse the information from different scalars as a data quality enhancer.To that goal, we have started from controlled situations, in which the amplitude and characteristics of noise are known and a ground truth is available.The consequences of this study go beyond a simple denoising algorithm, as this method settles the basis for the implementation of new dynamic equations based on a new dynamic, more robust entities coming from multifractal analysis.

Data
For the analysis developed in this paper the output of a global circulation model has been used.Our data set consists of output from the Ocean general circulation model For the Earth Simulator (OFES) (Masumoto et al., 2004;Masumoto, 2010).The simulation has been carried out with a model developed starting from Princeton GFDL MOM-3 (a z coordinate, eddy resolving, explicit free-surface model) by adapting it to the parallel architecture of the Earth Simulator.The simulation has been run on the near-global ocean, spun up for 50 yr under climatological forcing taken from monthly mean NCEP (United States National Centers for Environmental Prediction) atmospheric data.After that period the OFES is forced by the daily mean NCEP reanalysis for 48 yr from 1950 to 1998.See Masumoto et al. (2004) for additional details on the forcings.
The output of the model corresponds to daily data for the last 8 yr of simulation.Horizontal angular resolution is the same in both the zonal, φ, and meridional, θ, directions, with values of θ = φ = 1/10 • .The calculations are restricted to variables in surface layer, namely temperature and salinity, to emulate remote sensing data.The results shown in this paper refer to 1 January 1990 (see Fig. 1), although similar ones are obtained for any other day.

Noise sources
All noise sources n β (x) in this paper are Gaussianly distributed with zero mean and marginal standard deviation σ n which is the same at all the points.The spatial correlation matrix of the noise is given by the exponent β of its power spectrum, namely: The value of the exponent can be taken between β = 0 (purely decorrelated noise) and β = 2 (strongly correlated noise).In this study we take β as β = 0, 1 and 2.
3 Non-linear methods used in this study: theoretical review

Singularity analysis
The existence of a multifractal structure in synoptic maps of ocean scalars has been studied for several years now and can be given a geometrical meaning by means of the appropriate analysis.A number of studies have shown that singularity fronts forming a multifractal hierarchy can be extracted from SST maps of different types and resolutions (Turiel et al., 2005(Turiel et al., , 2008a)).It has also been observed that other scalars, such as chlorophyll concentration maps (Nieves et al., 2007) and even brightness temperature maps (Isern-Fontanet et al., 2007), have exactly the same singularity fronts that can be put into correspondence among the different scalars.The emergence of such a singular structure is the result of the advective forces acting on a quasi-2-D turbulence regime, and can be thus observed for scales ranging from kilometres to the planetary scale.Notice that 2-D turbulence is not a requirement for having singularity fronts: 3-D turbulence also generates singularity fronts, although they form compiles 2-D manifolds in general (Kestener and Arnéodo, 2003).Quasi-2-D turbulence is required to be able to identify singularity fronts derived from maps of surface variables with dynamic structures.Singularity analysis is the keystone of the so-called microcanonical multifractal formalism (MMF); in essence, the fluid can be understood as a hierarchical arrangement of different fractal components, each one with its own fractal dimensional and characterized by a particular value of the singularity exponent.The existence of scaling properties of scalars submitted to specific types of turbulent flow can be tracked back to Kolmogorov's seminal works (Kolmogorov, 1941a, b).Later works showed that structure functions in fully developed turbulence (2-and 3-D) present multiscale dependence on the scale parameter (Dubrulle, 1994;She and Leveque, 1994;She and Waymire, 1995).The intuition that such multiscale behaviour is the result of the hierarchy of multiple fractal components is due to Parisi and Frisch (Parisi and Frisch, 1985), although the explicit extraction of the fractal components of an ocean variable (just for the surface layer) was not done until Turiel et al. (2005).With the explicit calculation of the singularity exponents of an ocean scalar MMF can be applied to provide a dynamic interpretation of the resulting singular fronts.For a comprehensive presentation of MMF the reader is referred to Turiel et al. (2008b).
Obtaining the singularity exponents associated with a given point x of a scalar map s should be easy by just studying the neighbourhood of the point, as for simple scalars: where d is the dimension of the domain of s (d = 2 in our case) and r is any displacement vector, with r denoting its modulus.In practical terms, real signals s are affected by many artifacts, as noise, discretization errors and spurious long-range correlations.Hence, the calculation of singularity exponents requires the use of an appropriate interpolation scheme, the most usual one being wavelet projections.Let s be an arbitrary 2-D scalar signal.It will be said that s has a singularity exponent h(x) at the point x if, for any wavelet function (Mallat, 1999), the following relation holds: where r stands for an arbitrary scale parameter (normalized by the integral scale so it is dimensionless and smaller than 1) and o r h(x) is a term that becomes negligible compared to r h(x) when r goes to zero.The amplitude function α(x) does not depend on the particular scale at which the wavelet projection is calculated and has the same units as the gradient of the scalar.
The left-hand side, i.e.T |∇s|(x, r), is called the wavelet projection of the gradient modulus of s over the wavelet , and represents a local zoom of variable size around the point x.What is important in Eq. ( 3) is the function h(x), which is called the singularity exponent of the function s at the point x.The singularity exponent is, by construction, a dimensionless measure of the regularity or irregularity of the function s around the point x, extending the usual concepts of differentiability (integer values of h(x) have the regular meaning of the amount of continuous derivatives accepted by the function) and of Hölder exponents (Rudin, 1987;Mallat and Huang, 1992).
If a singularity exponent can be assigned at each point of a signal s, the signal is said to be multifractal in the mi-crocanonical sense (although strictly speaking some additional conditions should hold, see Turiel et al., 2008b).All the scalar variables studied in this paper are multifractal in the microcanonical sense, which can be verified in a standard way (Turiel et al., 2006).
Singularity exponents are related to ocean circulation and thus they are not specific to any particular scalar under study (as has been shown in the case of remote sensed chlorophyll concentration and SST in Nieves et al., 2007).In fact, it has been verified that to a good approximation singularity lines coincide with the streamlines of the flow (Turiel et al., 2009), confirming that singularity exponents are characterized by the flow and are not scalar-specific, as shown in Figs. 2 and 3. Singularity analysis is useful to track fronts or streamlines similarly to the methods described in Cayula and Cornillon (1995); Ullman and Cornillon (2000).
However, the coincidence of the isolines of singularity exponents associated with two different scalars does not imply that the values of singularity exponents of those two scalars coincide.It turns out that, at least with our output of OFES, the singularities from SST and from SSS coincide within the numerical accuracy of the singularity analysis methods and the uncertainties associated with the model dynamics.In Fig. 4 we show the conditioned histogram of the singularity exponents of SSS conditioned by the value of the singularity exponents of SST.The maximum probability line is given by h SSS = h SST − 0.08, while more than 95 % of the probability is concentrated within ±0.15 of this line, so empirically h SSS = h SST − 0.08 ± 0.15 with a good approximation.

Fusion method
The complete theoretical discussion of the fusion method can be found in Umbert et al. (2014); we summarize the main results here.
Let us suppose that we have two scalar variables, s and θ, that are multifractal according to MMF.It is additionally assumed that both scalars have the same singularity exponents at all points.We have hence: By using some theorems from functional analysis (extensions of Riesz's representation theorem) we conclude from Eq. ( 4) that the gradients of θ and s must be related by a 2 × 2 smooth matrix ρ, namely: (5) The matrix ρ (x) must be slowly varying, meaning that it has a small gradient; otherwise it will introduce new singularity exponents in s not present in θ, which contradicts our first hypothesis, Eq. (4).
Equation ( 5), although appealing, does not provide a constructive rule for calculating the matrix ρ from the data.In 3. Singularity analysis is useful to track fronts or streamlines similarly to the methods described in Cayula and Cornillon (1995); Ullman and Cornillon (2000).2014) a simplifying hypothesis is discussed: if this matrix is almost everywhere proportional to the identity matrix we can conclude that: where, as for ρ, the functions a and b must have small gradients.
Let us now suppose that the function s is contaminated by some source of noise n, so in fact our data are θ(x) (assumed to be noiseless) and s (x) = s (x)+n (x).We can construct a filtered version of s, denoted by s f , by applying the following relation: where â(x) and b(x) are estimates of the actual parameters a(x) and b(x) in Eq. ( 6).These estimates must be obtained from the values of s and θ by an appropriate method which takes into account the value of both variables and which enforces that both â and b must have small gradients almost everywhere.In Umbert et al. (2014) one such method was proposed, in which â and b are calculated by linear regressions weighted around each point (Nieves et al., 2007).The total weight of a point x, N (x), is defined as follows: For any function f we define its local average around the point x, f x , as: We have thus: and then Eq. ( 7) is applied to obtain the filtered signal s f .This algorithm was shown in Umbert et al. (2014) to lead to a large increase in the quality of SMOS SSS maps using OSTIA SST maps as templates, and a partial restoration of the structure of singularity exponents.The method, however, is far from perfect (the singularity exponents of the filtered signal s f and those of the template θ differ significantly), so further research must be conducted to improve the method.
The method was introduced directly to process remote sensing data, without any prior study on its sensitivity to noise and on the expected quality of the results.This paper intends to fill this gap by analysing in detail the performance of the method in front of several sources of noise of controlled spatial correlation and amplitude.

No noise
In the absence of noise, the fusion algorithm creates a map which slightly deviates from the original, as can be observed in Fig. 5.The map of singularity exponents, Fig. 6 shows that some structures have become a bit blurred.This is probably a consequence of the approximation implied by Eq. ( 7).The error map at the bottom of Fig. 5 has a standard deviation of 0.3.This is an indication that the error correction by our multifractal fused method cannot go below this threshold within the present formulation.

Decorrelated noise (β = 0)
Decorrelated noise is the easiest to filter away by our fusion algorithm.In Fig. 7 we show an example of white noise added to the original SSS map and of the result of applying the fusion algorithm to it, as well as the difference map.Significantly, the average error is negligible (0.02) and the standard deviation of the error map (obtained by averaging over all spatial locations) is 0.33 for input noise of unity standard deviation, a considerable reduction in the error, to very nearly the same level occurring when no noise is introduced.Nevertheless, the restoration of singularity exponents, Fig. 8, is very incomplete: just the singularity exponents associated with the main currents can be recognised.Notice however that the level of input noise is rather high.

Correlated noise (β = 1)
Correlated noise is still very well filtered away by our fusion algorithm.In Fig. 9 we show an example of noise added to the original SSS map with the spatial correlation of the noise proportional to k −1 and of the result of applying the fusion algorithm to it, as well as the difference map.
The error map has a spatial average of 0.02 and a standard deviation of 0.36 for input noise of unity standard deviation, a considerable error reduction slightly worse than in the previous case.The error map is still associated with significant geophysical structures.The restoration of singularity exponents, Fig. 10, is worse than in the previous case and centered on the main currents.

Strongly correlated noise (β = 2)
The correction on strongly correlated noise is very challenging for our fusion algorithm, with an important deformation of large-scale oceanic features.In Fig. 11 we show an example of noise added to the original SSS map with the spatial correlation of the noise proportional to k −2 and of the result of applying the fusion algorithm to it, as well as the difference map.
Strongly correlated noise introduces large-scale structures, which are not of geophysical origin, but even so the error map has a negligible spatial average of 0.014.Submitted to such a strong perturbation, the standard deviation of the error is larger in this case, arriving at a value 0.66 for input noise of unity standard deviation; significant but not quite satisfactory.The error map has a fuzzy structure, and the main resulting structures are void of geophysical meaning.The restoration of singularity exponents, Fig. 10, is even worse than in the previous case and even the main currents are hard to recognise.

Dependence of the reconstruction quality on noise intensity
To quantify the quality of our fusion algorithm, we define the output noise level, σ o , as the standard deviation of the error image, calculated over the whole image (in fact the average error is always lower than 0.02 and negligible by comparison, so the standard deviation and the root squared error are almost identical).We define the input noise, σ i , as the standard deviation of the injected noise.
In Fig. 13 we show the evolution of the ratios of output noise to input noise as the size of input noise increases for the three types of noise studied.We observe that there is a crossover value σ c below which fusion is not useful: for any input noise σ i < σ c the output noise σ o is equal to σ c (and then the ratio σ o /σ i diverges).On the other hand, when the level of input noise is large enough the capability of our fusion algorithm to remove noise stagnates at a fixed percentage of input noise (the asymptotic maximum percentages of noise removal are 91 % (β = 0), 83 % (β = 1), and 43 % (β = 2).Those ratios are really very good, implying that fusion is always a good strategy even for very noisy images.
Fig. 9. Top: SSS after the introduction of correlated noise with input noise of standard deviation σ i = 1 psu.Middle: SSS map after filtering the map on top using the SST map as a template.Bottom: Differences between the original (no noise) SSS map and the fused SSS map.

Comparison with median filters
Fusion methodologies have the appeal of invoking the use of synergistic approaches to combine the information coming from different variables; however, it is a relatively sophisticated methodology to filter noise which does not depend on signal amplitude, as the one used in this paper.Many conventional filters, as for instance median filters or average filters, could lead to better or comparable results and are easier to implement.We have thus compared the results of fusion filters with median filters (average filters always have worse performance than median filters).
To define our median filters, we will consider square neighbourhoods of (2N + 1) × (2N + 1) points (at the given resolution) centered at each point of the map.For each square neighbourhood we compute the median of the points in the neighbourhood (the value of the signal which has exactly 50 % probability of occurrence, that is, the one exactly in the middle position if the list of (2N +1)×(2N +1) values were ordered), and we assign this local median value to that point to obtain the median filtered map.Obviously the quality of this map depends on the size of the neighbourhood, which is parametrized by N. Larger Ns tend to reduce the impact of the overimposed Gaussian noise, but at the cost of increasing the error associated with smoothening of the original signal.In Fig. 14 we show an example of the application of median filters of different sizes to a map of SSS corrupted with spatially decorrelated Gaussian noise.
As in Sect.4.5 we have calculated the evolution of the ratio of output to input noise as the amplitude of input noise increases for median filters of different sizes, compared to the performance of fusion.As shown in Figs. 15, 16, and 17, median filters achieve good reconstruction only when the amplitude of noise is small enough.In general, median filters have better performance when noise is less correlated and when the filter size increases, at the cost of sacrificing resolution.At large sizes (5 × 5 or greater) median and average filters become very similar for input noises greater than 0.5 (figures not shown) because the median and the average of the Gaussian noise coincide.Taking this into account we can predict that the output noise for a median filter applied to a map contaminated with decorrelated noise with N ≥ 2 and σ i > 0.5 behaves as: where σ g (N ) is the inner geophysical variability of SSS in the (2N + 1) × (2N + 1) neighbourhoods.SSS is a variable with relatively low variability to scales of less than 1 degree (Jorda et al., 2011), so up to N = 5 the first term in Eq. ( 15) can be neglected as compared to σ i for σ i > 0.5.We obtain thus This explains why the curves in Fig. 15 asymptotically tend to straight lines, and why the asymptotical value is proportional to 1 2N+1 .When noise is correlated, it is possible to infer that the formula above can be generalized to: where γ → 0 as β → ∞.This explains why median filters have lower performance with correlated noise.In contrast, fusion offers a much greater reduction of noise, probably because it is a non-linear method.

Summary and future work
In this paper we have shown that ocean scalars of different types possess a common multifractal structure.This correspondence has been highlighted in previous works (Nieves et al., 2007;Isern-Fontanet et al., 2007), but it has only been used recently (Umbert et al., 2014) to derive a new algorithm for data fusion.
Our multifractal fusion method can be used to enhance the signal quality of a variable contaminated with noise by using another variable as a template.The method is constructed to reintroduce the multifractal structure of the noised variable based on that of the template.This guarantees the reestablishment of the structure of short-range spatial correlations in the signal, but it seems to impact also long-range correlations.
We have shown the effect of the fusion algorithm when a source of uncorrelated noise of given amplitude is added.Our results indicate that until a threshold in noise amplitude of about 0.3 is reached the application of our denoising algorithm does not improve the quality of the signal; even a noiseless image degrades in quality to 0.3 after the application of the fusion algorithm.The results of the application of The estimated asymptotic values for the input/output ratios are 0.09, 0.17, and 0.57, respectively.The crossovers (passings by 1) occur for input noise levels of 0.33, 0.33, and 0.40, respectively.the fusion algorithm to remote sensing data (SMOS SSS maps using OSTIA SST templates) discussed in Umbert et al. (2014) show that in fact the quality could be even better, with post-fusion errors of around 0.20 and even smaller.This difference may indicate that the correspondence in the multifractal structure of the maps of SSS and SST resulting from OFES is not as good as in real data from the oceans.The result indicates that reinforcing the correspondence between ocean scalars in numerical models of the ocean could potentially improve their assessment quality about the actual ocean state, which is important both for operational applications and for the study of ocean processes.The introduction of an assimilation scheme accounting for the singularity structure of the scalars should be studied.Fig. 17.Evolution of the ratio output noise/input noise depending on the value of input noise for strongly correlated noise (β = 2).Red curve: fusion; green curve: median filter with 3×3 neighbourhoods; blue curve: median filter with 5 × 5 neighbourhoods; purple curve: median filter with 11 × 11 neighbourhoods.For large input noises fusion performs better than the others; the crossover points occur for input noise levels of about 0.5, 0.55, and 0.65, respectively.
The existence of a common singularity structure is a central hypothesis of our algorithm, which seems to be confirmed for more and more examples from remote sensing and from numerical models.The development of an appropriate theoretical framework to explain the creation of singularities is one of the major challenges in this field of research, but will help to unify the evolution equations for different scalars and most probably control their numerical behaviour and improve their quality.
The method presented in Umbert et al. (2014) and analysed here is just a first step in the theory of multifractal methods for data fusion.From the experiences shown in this paper it is clear that the method is not fully able to restore the multifractal structure using that of the template.Most likely this limitation comes from the approximation used (the matrix issued from Riesz's representation theorem considered to be proportional to the identity matrix).The extension for a representation with a non-trivial matrix and a better inference strategy for the smooth matrix elements should be actively sought.This matrix will not only serve to improve the quality of existing remote sensing maps, but also will serve to extrapolate maps in missing areas and will provide geophysical insight into the dynamic relations governing the relation between the different variables (in particular, providing canonical classification of regions according to the physical and chemical properties of their water masses -mostly reminiscent of the biogeochemical provinces).

Fig. 4 .
Fig. 4.Histograms of h SSS conditioned by h SST as derived from OFES simulation for the whole globe on 1 January 1990.For each column of the histogram the value of the SST-related variable is fixed, so the values in that column represent the distribution of probability of values of the SSS-related variable for that fixed value of the SST-related variable.To enhance visualization, in this figure the histogram is normalized by its maximum probability in each row (so the brightest colour in each row represents the mode or maximum of probability along that row).

Fig. 5 .
Fig. 5. Top: Original SSS.Middle: SSS map after filtering using the SST map as a template.Bottom: Map of differences.

Fig. 7 .
Fig. 7. Top: SSS after the introduction of decorrelated noise with input noise of standard deviation σ i = 1 psu.Middle: SSS map after filtering the map on top using the SST map as a template.Bottom: Differences between the original (no noise) SSS map and the fused SSS map.

Fig. 8 .
Fig. 8. Top: Singularity exponents associated with SSS after noise of unity amplitude and β = 0. Bottom: Singularity exponents derived from the fused SSS map.

Fig. 10 .
Fig. 10.Top: Singularity exponents associated with SSS after noise of unity amplitude and β = 1.Bottom: Singularity exponents derived from the fused SSS map.

Fig. 11 .
Fig. 11.Top: SSS after the introduction of strongly correlated noise with input noise of standard deviation σ i = 1 psu.Middle: SSS map after filtering the map on top using the SST map as a template.Bottom: Differences between the original (no noise) SSS map and the fused SSS map.

Fig. 12 .
Fig. 12. Top: Singularity exponents associated with SSS after noise of unity amplitude and β = 2. Bottom: Singularity exponents derived from the fused SSS map.

Fig. 13 .
Fig.13.Evolution of the ratio output noise/input noise depending on the value of input noise.Red curve: decorrelated noise (β = 0); green curve: correlated noise (β = 1); blue curve: strongly correlated noise (β = 2); a constant level of 1 is also drawn (in black) for reference.The estimated asymptotic values for the input/output ratios are 0.09, 0.17, and 0.57, respectively.The crossovers (passings by 1) occur for input noise levels of 0.33, 0.33, and 0.40, respectively.

Fig. 14 .
Fig. 14. Results of applying median filters to a noisy salinity with no spatial correlation and amplitude 1 (same as in Fig. 7, top).The filters are calculated over boxes of variable sizes centered around each point.Results are for 3 × 3 boxes (top), 5 × 5 boxes (middle), and 11 × 11 (bottom).