A survey on multiscale mollifier decorrelation of seismic data

In this survey paper, we present a multiscale post-processing method in exploration. Based on a physically relevant mollifier technique involving the elasto-oscillatory Cauchy–Navier equation, we mathematically describe the extractable information within 3D geological models obtained by migration as is commonly used for geophysical exploration purposes. More explicitly, the developed multiscale approach extracts and visualizes structural features inherently available in signature bands of certain geological formations such as aquifers, salt domes etc. by specifying suitable wavelet bands.


Introduction
The success of any geophysical exploration project is determined by the quality of the available exploration data. These data have to be analyzed and interpreted as accurate as possible. During the realization of exploration projects, project managers are confronted with different kinds of data. The available data consist of density data sets from gravimetric surveys and/or data from magnetic surveys, but mainly from migration results from conducted seismic exploration.
In Freeden and Blick (2013), a novel mollifier technique for post-processing of exploration data as an improved interpretation technique based on the concept of a "geophysically relevant" wavelet construction was proposed. The concept for this technique goes back to an idea discussed in Freeden and Schreiner (2006) to obtain mollifier approximations of the deflections of the vertical in gravitational theory. Freeden and Blick (2013) was the point of departure for a series of mollifier approaches in gravimetry, magnetometry, and seismics in Freeden's group. In the case of gravimetry, the mollifier method has been worked out and was realized by Freeden (2021), Freeden and Bauer (2020), Freeden and Nashed (2020), Freeden and Sansò (2020), Berg et al. (2020), Blick et al. (2018a), Freeden and Nashed (2018c), Freeden and Nashed (2018d), Blick et al. (2017), Möhringer (2014), and Freeden and Gerhards (2013) (see also the references therein). A similar technique applied to migration results via the Helmholtz equation is discussed in Blick et al. (2018b) and Augustin et al. (2014). For a summary of both methods, the reader is referred to Freeden et al. (2019), Freeden and Nutz (2015), and Freeden (2010). Further approaches were discussed in Blick (2015) for the acoustic wave equation and in Blick and Eberle (2019) for the static Cauchy-Navier equation. As a recent contribution, Freeden (2021) and Blick et al. (2021) also applied the mollifier approach to magnetometry, in addition to gravimetry. The presented paper gives a survey about perspectives in the field of seismics.
It should be noted that the elastic potential as discussed in this paper is not used for inversion purposes, instead, we apply the methodological framework only for post-processing of already existing geological models to establish a better interpretability. However, an inversion process can be applied as discussed e.g., in Freeden (2021), Freeden and Nashed (2020), Freeden and Nashed (2018d), Freeden and Nashed (2018e), as well as in the monograph (Freeden and Nashed 2018a) and more specifically, in the articles (Freeden and Nashed 2018b) and (Freeden and Nashed 2018c) contained therein.
Following some considerations in Blick and Eberle (2019) for density data modeling, we develop a method which enables us to specify also particular directional characteristics of a migration result. In order to establish such a decorrelation technique, we mathematically make the transition from the Helmholtz equation to the elasto-oscillatory Cauchy-Navier equation. In doing so, on the one hand, we leave the classically motivated Helmholtz approach, on the other hand, we are able to detect specific directional features by the elastic integral due to the tensorial nature of the fundamental solution.

Signature decorrelation based on the Helmholtz equation
In what follows, we summarize the results of the Helmholtz approach as given by Blick (2015) in order to present the idea of a physically relevant wavelet construction using Helmholtz wavelets. For the proofs of those results, the reader is referred to the cited literature.
The start of the wavelet development is the representation of the solution U of the Helmholtz equation in a regular region B, i.e., an open, bounded, and connected set in R 3 that contains the origin and whose closed, compact and locally c (2) -smooth boundary ∂B is free of double points and divides the space R 3 into the inner space B and the outer space B ext = R 3 \B. Here, α(x) denotes the solid angle in x subtended by the boundary ∂B.
If B is a regular region, we have Further, let F : B → R be a Hölder continuous function. For a fixed k 0 ≥ 0, the Helmholtz integral equation relating a given contrast function F (i.e., acoustically given seismic data or a migration result) by convolution against the Helmholtz fundamental solution to the Helmholtz potential U is a solution of Eq. (1). By mollification of G(Δ + k 2 0 ; ·), we obtain a "mollifier potential scaling function" G τ (Δ + k 2 0 ; ·) with scale parameter τ > 0. Helmholtz differentiation results in the Haar-like "mollifier source scaling function" Φ τ (k 0 ; ·) = (Δ + k 2 0 )G τ (Δ + k 2 0 ; ·). Hence, even if our main focus lies on the application of the source scaling function, their theoretical construction always starts on the potential level.
Our considerations are summarized in the following theorem (cf. Müller 1969).
Theorem 1 If F is Hölder continuous in B, then the Helmholtz differential equation holds true for all x ∈ B. Here, α(x) denotes the solid angle in x subtended by the boundary ∂B.
Summarizing our considerations, we are led to a result, which builds the theoretical basis for our approach of geological feature extraction in a post-processing process (cf. Blick 2015).

Theorem 2
The "τ -mollifier Helmholtz potential functions" and the "τ -mollifier Helmholtz contrast functions" provided that F is Hölder continuous in B.
In fact, we have that both the real and imaginary part of U − U l τ converges with order O(τ 2 ) as τ tends to zero.
A similar approach for the elasto-static Cauchy-Navier potential for the decorrelation of density data is discussed in Blick and Eberle (2019). For those cases, it could be shown that the source scaling functions are normalized, i.e., the respective integral in R 3 is equal to one for all τ > 0. Note that this is a helpful feature used in the theory of singular integrals, since the integral kernels are Dirac sequences (for more information, see e.g., Stein 1971;Hörmander 1998 andWienholtz et al. 2009). This is not the case here. Nevertheless, we are led to the next theorem (cf. Blick 2015).

Theorem 3 The volume integral of
The theorem is of special importance, since Eq. (5) is used to normalize the wavelets developed later on in this paper.

Signature decorrelation based on the elasto-oscillatory Cauchy-Navier equation
Our task now is to adopt the Helmholtz scheme for the elasto-oscillatory Cauchy-Navier equation. We again start with the potential equation, and more specifically, its solution via fundamental solutions. Introducing C (n) (B) and c (n) (B) as the spaces of n-times continuously differentiable tensor and vector functions in B, respectively, we adopt a scaled version of the setup introduced by Kupradze (1979) and present the classical theory of elasticity for homogeneous and isotropic media. A homogeneous and isotropic elastic medium is composed of a regular region B of the three-dimensional Euclidean space and a set of the constant quantities ρ 0 , λ 0 , and μ 0 satisfying the conditions where ρ 0 is the constant density of the medium and λ 0 , μ 0 are the constant Lamé parameters. Using scaled Lamé parameters to enhance readability, we start with the following definition.

Definition 1
The oscillatory state of the medium B(ρ 0 , λ 0 , μ 0 ) corresponding to the acceleration f is the pair [u, σ ], which satisfies the conditions iii) with f ∈ R 3 , u ∈ C 3 , σ ∈ C 3×3 and the oscillation frequency ω ∈ R.
Substituting (7) into (6), we obtain the (scaled) oscillation equation (of classical elasticity) or the (scaled) oscillatory Cauchy-Navier equation of the medium B(ρ 0 , λ 0 , μ 0 ) corresponding to the acceleration f in terms of the displacement components where By use of the matrix differential operator we can reformulate Eq. (8) as We adopt the (tensor) fundamental solution G( A(∂, ω), x) of Eq. (9) from Kupradze 1979. G( A(∂, ω), x) is given by where It is easy to calculate that All in all, we are confronted with a potential that represents a solution of the Cauchy-Navier equation (9).

Remark 1
Since we only talk about the elasto-oscillatory Cauchy-Navier equation in this paper, we will just use the term Cauchy-Navier equation and drop the elastooscillatory part in order to keep the nomenclature short.

Remark 2
The Cauchy-Navier equation and the mollification of its fundamental solution can be applied in a number of other applications such as wave inversion (e.g., seismic imaging, Aki and Richards 2002) and specialized wave propagation (e.g., coupling of interior/exterior wave propagation problems, Eberle 2018Eberle , 2019Eberle , 2020. Contrary to the construction of scaling functions for the scalar fundamental solutions discussed so far, we are faced with the problem, that we cannot simply obtain a mollification in a ball by application of the Taylor expansion as already presented in Blick and Eberle (2019). In order to mollify the fundamental solution, we rewrite Eq. (10) in terms of G(Δ + k 2 0 ; ·). Hence, where I is the 3 × 3 identity tensor, denotes the fundamental solution w.r.t. the Laplace operator and is formally included by setting k 0 = 0 in G(Δ + k 2 0 ; ·). Now, we are able to introduce the mollification G l τ ( A(∂, ω); ·), l = 1, 2, by substituting G(Δ + k 2 m ; |x|) by G l τ (Δ + k 2 m ; |x|), l = 1, 2, as well as G(Δ; |x|) by G 1 τ (Δ; |x|). Note that we exchange G τ (Δ; |x|) only with the partial taylorization G 1 τ (Δ; |x|) instead of G l τ (Δ; |x|), l = 1, 2. The reason is twofold. First, by only taking the case l = 1, we keep the number of terms small especially regarding that we have to take the fourth power of G l τ (Δ; |x|). Secondly, G 1 τ (Δ; |x|) is the mollified Laplace fundamental solution discussed in Blick (2015) and Blick and Eberle (2019) and hence should be applied to mollify G(Δ; ·).

Remark 3
Here, we lose the radial symmetric property of G(Δ + k 2 0 ; ·) due to the application of the Hessian matrix in the definition of G( A(∂, ω); ·). However this is advantageous for decorrelation purposes, since we are now able to better highlight the spreading direction of the layers contained in the data.
Due to the continuous differentiability of G l τ (Δ + k 2 0 ; ·), we directly obtain the following lemma.
As in the Helmholtz case, the "τ -mollifier Cauchy-Navier potential functions" converge to the Cauchy-Navier potential u.
Theorem 4 Suppose that B is a regular region in R 3 . Further on, let x ∈ B be arbitrary and let f : B → R 3 be continuous. Then, Hence, with and assuming 0 < τ < 1, we cut out a small ball B τ 2 (x) of radius τ 2 around the center x from the integral, so that we get Next, we split f 0 into its positive and negative parts f + 0 and f − 0 , respectively and split (G( A(∂, ω); ·) − G l τ ( A(∂, ω); ·) into its real ( ) and imaginary part ( ). Hence, we have as τ → 0+. The same steps lead to equivalent results for f − 0 , so that we obtain Now, we take a look at deriving the Cauchy-Navier source scaling function corresponding to the potential scaling function G l τ ( A(∂, ω); ·) τ >0 by applying the operator A(∂, ω) to G l τ ( A(∂, ω); ·).
Theorem 5 For x, y ∈ R 3 , we find is called the mollifier Cauchy-Navier source scaling function and is given by with T l 1 , T l 2 , T l 3 and T l 4 as in (15) -(16), respectively, and O denotes the zero tensor.
The terms T l 1 , T l 2 , T l 3 and T l 4 where computed by use of the MATLAB Symbolic Math Toolbox and are given by Similar to the Helmholtz source wavelets, Φ l τ (ω; ·) is singular in the imaginary part for l = 1, and Φ l τ (ω; ·) is singular in both real and imaginary parts for l = 2. Nonetheless, both scaling functions are integrable and the analytical disadvantage can be easily overcome by a numerical integration scheme as discussed in Blick (2015). Figures 2 and 3 depict profiles of the scaling functions G l τ ( A(∂, ω); ·) and Φ l τ (ω; ·).
Remark 4 Taking a look at Fig. 1 and comparing it to the cuts depicted on the right side of Figs. 2 and 3, it is easy to observe the close relationship of the Cauchy-Navier and Helmholtz scaling functions.
Lemma 2 For all τ > 0, the so-called volume integral of the mollifier scaling function Φ l τ (ω; ·) given by Proof By use of polar coordinates and by virtue of elementary calculation, we find Substituting β m , α m and k m as given in Eqs. (11)-(12) and taking the limit τ → 0+ via l'Hospital's rule provides the desired result. The scaling function Φ l τ (ω; ·) allows us to define the "τ -mollifier Cauchy-Navier contrast functions" f l τ given by This leads to the following theorem.

Theorem 6
The function Φ l τ (ω; ·) is piecewise continuously differentiable in R 3 \{0}. In addition, assume that B is a regular region in R 3 and that f : B → R 3 is continuous holds true for all x ∈ B, where α denotes the solid angle in x subtended by the boundary ∂B.
Proof We split Φ l τ (ω; ·) into its real and imaginary parts so that We split Φ l τ (ω; ·) i j into its positive and negative parts, so that where (·) + and (·) − denote the positive and negative parts of the corresponding functions, respectively. Since f is continuous and Φ l τ (ω; ·) + i j as well as Φ l τ (ω; ·) − i j are integrable and do not change the sign, the mean value theorem of integration guarantees the existence of ξ 1 , ξ 2 ∈ B ∩ B τ (x), so that According to Lemma 2 and observing that Φ l τ (ω; ·) is point symmetric and that ∂B as the boundary of a regular region is locally c (2) -smooth, we have Substituting the last equation into Eq. (17), we get Now we need to estimate the integral on the right hand side. Exemplary taking a look at the case l = 1 (cf. Eq. (14)), we split the exponential function into cosinus and sinus and see that Φ l τ (ω; x − y) + i j is composed of the following functions: 1. Terms of the form 0 ≤ cos(rk m )r n τ m−15 for n + m ≥ 12 and 0 ≤ r ≤ τ ≤ π 2k m contribute to both the real and imaginary part. With polar coordinates, we get τ 0 cos(rk m )r n+2 τ m−15 dr = τ m+n−12 where p H q (a, b, z) denotes the generalized hypergeometric function (Magnus et al. 1966;Luke 1969) of order p, q. Here, p denotes the length of the vector a and q the length of the vector b. We have with Γ as the Gamma function denotes the Pochhammer symbol. It should be noted that p H q (a, b, z) is convergent for |z| < ∞ if p ≤ q. Hence for n + m ≥ 13, we have lim τ →0+ τ 0 cos(rk m )r n+2 τ m−15 dr = 0 and for n + m = 12, the integral is bounded.
2. Observing that β 1 = −β 2 in Eq. (14) and k 2 ≥ k 1 , we find terms of the form 0 ≤ (cos(rk 1 ) − cos(rk 2 )) r n τ m−15 for n + m = 10 and 0 ≤ r ≤ τ sufficiently small. Integration leads to lim τ →0+ τ 0 (cos(rk 1 ) − cos(rk 2 )) r n+2 τ m−15 dr = τ m+n−12 where C is a constant. These terms contribute to the real part. 3. For the imaginary part, we have terms of the form 0 ≤ sin(rk m )r n τ m−15 for n + m ≥ 11 and 0 ≤ r ≤ τ ≤ π k m . Integration yields which for τ → 0+ converges to zero in the case of m + n ≥ 12 and is bounded for m + n = 11. 4. Terms of the form 0 ≤ C k m cos(rk m )r n τ 11−m + sin(rk m )r n−1 τ 11−m for τ ≤ π/(2k m ) and C a positive constant contribute to the imaginary part. Observing (19) and (20) All in all using the triangle inequality, we can now state that for l = 1, both the real and imaginary part of the integral of Φ l τ (ω; x − ·) + i j are bounded for τ → 0+. Similar observations for the case l = 2 show that the real part of the integral of Φ l τ (ω; x − ·) + i j is bounded and the integral of the imaginary part is of order O τ −1 as τ → 0+. Hence, we have for l = 1, 2, since f is continuous. Following the same ideas, it is easy to show that for l = 1, we find if f is continuous. Hence, it remains to prove the imaginary case for l = 2. However, if we want to realize that the imaginary part for l = 2 also converges to zero, we have to assume that f ∈ c (1) (B).
In order to continue the proof, we present estimates only used in the case l = 2. Here, we observe, that any function f ∈ c (1) (B) can be written as It is easy to explicitly calculate that the theorem holds true for all linear functions due to the point symmetry of the scaling functions. Hence, we can assume without loss of generality, that D α f (x) = 0 for |α| ≤ 1 and x ∈ B fixed. In that case, we apply Taylor's theorem to approximate f j (ξ i ) by expansion around x, so that as τ tends to zero. All in all, combining these estimates with the observation that as τ → 0+, we follow the previous considerations up to Eq. (18), so that which proves the theorem.

Wavelet representation of potential and contrast functions
Next, we deal with the mathematical mechanisms for interpretation and understanding of available seismic information inside a regular region B. In order to do that, we again take a look at the Helmholtz case in order to make the reader more familiar with the idea of seismic decorrelation and introduce the associated notation. Our purpose is to demonstrate, how the multiscale procedure for the potential canonically transfers to seismic data by use of "Helmholtz derivatives" as shown in Blick (2015). As already noted, the source wavelets in this paper do not have zero mean (neither in the real part, nor the entire function). Since the wavelets are constructed by taking the difference of scaling functions with different scales, this results in the following two cases which we want to avoid.
(k 0 ;·) is large. Calculating the difference, i.e., taking a look at the band-pass filtered signal, the result will be dominated by the low-pass filtered signal generated by filtering with Φ l τ 2 (k 0 ; ·) and hence, we do not get any new information from the decorrelation.
(k 0 ;·) positive, then the band-pass filtered version is in fact the summation instead of the subtraction of the information contained in the two lowpass filtered signals.
Hence, we normalize by the real part of the volume integral of the source scaling functions in order to construct wavelets that have zero mean in the real part. It would also be possible to construct the wavelets Ψ l τ j (k 0 ; ·) in such a way that the entire wavelet has zero mean. However, this leads to a change in phase, which we want to avoid.
Seen from a numerical point of view, it is remarkable that both wavelet functions y → W l τ j (Δ + k 2 0 ; |x − y|) and y → Ψ l τ j (k 0 ; |x − y|) vanish outside a ball around the center x due to their construction. These functions are space-limited showing a ball as local support. Furthermore, the ball becomes smaller and smaller with increasing scale parameter j, so that more and more high frequency phenomena can be highlighted without changing the features outside the ball.
The associated "τ j -mollifier potential wavelet functions" and the "τ j -mollifier contrast wavelet functions" are given by The τ j -potential wavelet functions and the τ j -contrast wavelet functions, respectively, characterize the successive detail information contained in U l In other words, we are able to recover the potential U and the contrast function, i.e., the "geological signatures" F, respectively, in form of "band structures" and As a consequence, the essential problem to be solved in multiscale extraction of geological features is to identify those detail information, i.e., band structures, which specifically contain desired geological characteristics, for example, aquifers and salt domes.
Thus for x ∈ B and F ∈ C (0) (B), we finally end up with the following multiscale reconstruction and In addition, if F is Hölder continuous, we have All in all, the potential U as well as the contrast function, i.e., the "geological signature" F can be expressed in an additive way as a low-pass filtered signal U l τ 0 and F l τ 0 and successive band-pass filtered signals (WU ) l τ j and (W F) l τ j , j = 1, 2, . . . , respectively. It should be mentioned that this multiscale approach is constructed such that, within the spectrum of all wavebands (cf. Eqs. (21) and (22)), certain rock formations may be associated to a specific band within the multiscale reconstruction. Each scale parameter in the decorrelation is assigned to a data distribution which corresponds to the associated waveband and thus, leads to a low-pass approximation of the data at a particular resolution. The wavelet contributions are obtained as part of a multiscale approximation by calculating the difference between two consecutive scaling functions. In other words, the wavelet transformation (filtering) of a signal constitutes the difference of two low-pass filters, thus it may be regarded as a band-pass filter. Due to our construction, the wavelets show an increasing space localization as the scale increases. In this way, the characteristic signatures of a signal can be detected in certain frequency bands.
In the same way as for the Helmholtz equation, we obtain a multiscale procedure for the potential u as well as the contrast function f ∈ c (l−1) . Again suppose that {τ j } j∈N 0 is a positive, monotonically decreasing sequence with lim j→∞ τ j = 0 and with V Φ l τ j (ω;·) = 0 for all j ∈ N 0 . We consider the differences . W l τ j ( A(∂, ω); ·) and Ψ l τ j (ω; ·) are called "mollifier Cauchy-Navier potential wavelet function" and "mollifier Cauchy-Navier source wavelet function".
The associated "τ j -mollifier Cauchy-Navier potential wavelet functions" and the "τ j -mollifier Cauchy-Navier contrast wavelet functions" are given by and (a) 2D cut of the 3D Marmousi migration result.
(b) Interpretation of the Marmousi model.  Martin et al. 2006) the decorrelation ability of the wavelets. We specifically accept the downside of the model that it is constant in the x 2 direction in favor of the fact that we have a complete interpretation of the model as depicted in Fig. 4. The essential ingredients of the decorrelation applied to the Helmholtz equation are illustrated for 2D cuts of a 3D version of the Marmousi migration result by the scheme depicted in Fig. 5. The figure shows that the Helmholtz potential does not yield usable information about the geology because of its smoothness (Fig. 5, top left), since the differences may be understood as a discrete version of Eq. (1). However, if we go over to the multiscale decorrelation in terms of band-pass filtered data U τ j − U τ j−1 of the potential U at scale j, structural information of the Marmousi model become visible (Fig. 5, bottom left). The key idea of the method is that the potential wavelets generating the potential decorrelation (Fig. 5, left) can be correlated via the application of the Helmholtz operator to Haar-type wavelets in the migration level (right). In fact, the decorrelation of the migration result in band signals (Fig. 5, bottom right) clearly shows that information about the geological interfaces becomes available.
The aforementioned wavelet construction is particularly powerful because of its "geophysical relevance", i.e., it forms a compromise reflecting the underlying physics (in accordance with the underlying differential equation) while still delivering an adequate multiscale decorrelation of geological signatures. Transitions of geological  Fig. 5 Schematic visualization of the multiscale decorrelation mechanism as depicted in Blick et al. (2018b) strata can be detected, but of course not the specific geological formations themselves. Nevertheless, the structure of the geological configuration and even the faults become visible by mollifier decorrelation. The wavelets employed for establishing the Helmholtz scheme constitute radial basis functions, so that they are only dependent on the mutual distance of two points of the area under investigation. This means that no specific directionally-reflected information can be verified by the model described above.

Decorrelation using Helmholtz wavelets
We are more interested in the decorrelation of the source data, even if our mathematical setup allows us to also decorrelate the potential. We borrow the results of the Helmholtz case from (Blick 2015) and present them here as a reference for the comparison to the Cauchy-Navier decorrelation scheme.
We concentrate on the real part of the decorrelation of the 3D synthetic Marmousi migration result for k 0 = 0.036 rad m based on a conducted parameter study in Blick (2015). The range of k 0 for the study is given by 0.0018 rad m , 0.1 rad m and results from the choice of an excitation frequency between 10 and 150 rad s for occuring rock velocities between 1500 and 5500 m s . For a comparison to the decorrelation for other values of k 0 , the reader is also referred to (Blick 2015). The sequence τ j is adopted from (Blick 2015) as well. The sequence is constructed in such a way that for scales j = 1, . . . , 5, the real part of V Φ 2 τ (k 0 ;·) is close to one. Since a smaller spacing of τ j is also required for large scale parameters j, a dyadic sequence for scales j = 6, . . . , 9 is taken.

Remark 5
The depicted figures in this section denote 2D cuts of the complete 3D convolution result of the 3D Marmousi migration result as used in Blick (2015).
A decorrelation of the 3D Marmousi migration result F via the partially and fully taylorized normalized Helmholtz source scaling functions and wavelets is illustrated in Figs. 6, 7, 8 and 9 for k 0 = 0.036 rad m . We observe that the convolution with the partially taylorized normalized Helmholtz source scaling function and wavelet tends to oscillate heavily and hence is blurred due to ghost images for low scales. This phenomenon gets worse if k 0 is increased. Comparing with Fig. 4b, we notice that the salt dome in the lower left part of the migration result and some coarser structures are clearly represented in the band-pass filtered signal at scales j = 3.
In comparison, the decorrelation via the fully taylorized normalized Helmholtz source scaling function and wavelet does not oscillate as heavily. As such, the characteristics of the salt dome can be seen at scales j = 4 in the band-pass filtered data. The structure changes for increasing scales j for both types of scaling functions. The dominating layer at scale j = 6 is the shale layer directly above the salt dome. Moreover, we find finer layers in the band-pass filtered signal such as sand layers for higher scales.
Another useful property is that the noise in the upper left corner can be filtered out with both methods. Here, the band-pass filtered signal of the partially taylorized (h) Band-pass filtered migration result at scale = 9. j Fig. 9 Real part of the multiscale approximation of the Marmousi migration result F by convolution with the normalized and fully taylorized (l = 2) Helmholtz source scaling function and wavelet for k 0 = 0.036 rad m source scaling function does not contain the noisy data for all depicted scales, whereas the band-pass filtered signal of the fully taylorized source scaling function is noiseless starting with scales j = 4.

Decorrelation using Cauchy-Navier wavelets
For the Cauchy-Navier wavelet decorrelation, we have to prepare the data set beforehand. The reason is, that the wavelets are based on the Cauchy-Navier equation, which models a homogeneous medium with constant Lamé parameters λ 0 and μ 0 , as well as a constant density ρ 0 . Hence, we choose a reference medium, in this case sandstone, with the parameters ρ 0 = 2066.38 kg m 3 , λ 0 = 1.9 · 10 9 Pa and μ 0 = 6.3 · 10 9 Pa (Gopalakrishnan 2016) to decorrelate the migration result. Hence, we have λ = 9.1948 × 10 5 m 2 s 2 and μ = 3.0488 × 10 6 m 2 s 2 . We use the frequency ω = 95.3 rad s so that k 1 equates to k 0 = 0.036 rad m used in the Helmholtz case. Further, we need vectorial input data f . Taking f = Fe i , i = 1, 2, 3, where {e 1 , e 2 , e 3 } denotes the canonical orthonormal system in R 3 , the decorrelation will be concentrated along the x i direction in the i − th component of the resulting vector. Moreover the remaining two components of the solution vector contain the decorrelation in the diagonal direction of x i and x j , j = i. Hence, for a complete decorrelation, we choose f = (1, 1, 1) T F = (e 1 + e 2 + e 3 )F (cf. Fig. 10) and present the resulting decorrelation in the form of convolutions * by We are now prepared to take a look at the decorrelation depicted in Figs. 11, 12, 13 and 14 using Cauchy-Navier wavelets and concentrate on its real part. We only present selected scales of low-pass and band-pass data since a full decorrelation would be to extensive for the scope of this paper. The Cauchy-Navier wavelets and scaling functions are able to highlight the same structural information as the corresponding decorrelation via Helmholtz wavelets. More specifically, both, the partially and the fully taylorized Cauchy-Navier wavelets and scaling functions can filter out the noise in the upper left corner of the dataset as can be seen by taking a look at the low-pass filtered data in Fig. 11 for the partially and Fig. 13 for the fully taylorized case.
Both methods present noiseless band-pass filtered data starting at scale j = 4 (cf. Figs. 12 and 14) which also contain the structural information of the salt dome on the lower left side in the pictures on the diagonal. The shale layers discussed in the Helmholtz case can also be found in Figs. 12 and 14.
As pointed out numerous times throughout this paper, we developed the decorrelation based on the Cauchy-Navier equation in order to highlight directional characteristics of migration results.
It turns out, that the partially taylorized Cauchy-Navier source wavelets and scaling functions oscillate heavily inside its support especially for large ω. As such, they behave similarly to the partially taylorized Helmholtz source wavelets and scaling functions. In both the Helmholtz and the Cauchy-Navier case, this may result in ghost images in the decorrelation for small scales j.
Further, directional characteristics in the data can clearly be highlighted via the Cauchy-Navier wavelets and source scaling functions. For example, Figs. 11 and 13 show structures extending along the directions of x 1 = x 2 and x 1 = −x 2 in the off-diagonal entries in both the low-pass and band-pass filtered signal. This behavior continues along all depicted scales, but is more pronounced for lower scales due to the size of the support of the wavelets and scaling functions.
It is clear that the mixed directions involving the direction x 3 are close to zero with an amplitude of 10 −12 and contain only numeric noise. This is due to the symmetry of Φ l τ j (ω; ·) and that f is constant in the x 3 direction. The major improvement compared to the decorrelation using Helmholtz wavelets is, that we can further decorrelate the model so that we can precisely highlight structures stretching in the x i direction, as well as diagonally in the data set. Therefore, an interpretation of a given migration result using Cauchy-Navier wavelets will be more precise then by applying Helmholtz wavelets. We discussed the scalar wavelet decorrelation using mollifier Helmholtz wavelets and presented an approach for the directional interpretability by tensorial mollifier Cauchy-Navier wavelet and scaling functions. Decorrelation results were presented and compared for both types of wavelets. The mollifier Helmholtz wavelets lead to a good decorrelation of the data, which aids the interpretation effort in exploration. The mollifier Cauchy-Navier wavelets however, allow a more precise decorrelation since the data set is even further disassembled so that structures stretching in either of the three Cartesian directions as well as all diagonal directions can be highlighted. In fact, the decorrelation via the mollifier Cauchy-Navier source wavelet yields further information on the migration result by highlighting signatures spreading in the three Cartesian coordinate directions on the diagonal and the mixed directions on the offdiagonal of the decorrelation.