© ESO, 2017

1. SOLA helioseismic inversions in a nutshell

The solar interior is filled with waves travelling between various places. In the convective envelope, the pressure (p) modes largely dominate the spectrum. The propagation of the waves is affected by perturbations in the plasma state parameters, by magnetic fields and last but not least by plasma streaming. The time-distance helioseismology comprises a set of tools used to measure and analyse the wave travel times. In this study we focus on difference travel times, that is, the difference of the measured travel time of waves travelling in the opposite directions. The difference travel times in the quiet Sun regions are mostly sensitive to flows (Burston et al. 2015).

The standard time-distance helioseismic pipeline consists of the following consecutive steps: first the spatio-temporal datacube is prepared using the tracking and mapping pipeline, this datacube is spatio-temporarily filtered to retain only waves of interest and subsequently the travel times are measured from cross-correlations of the filtered signal at two places. These travel times are finally inverted for flows assuming the linear relation between flow vector v^true and the measured difference travel time δτ via the sensitivity kernel K^a (coming usually from forward modelling, e.g. Birch & Gizon 2007; Burston et al. 2015); $\mathit{\delta }{\mathit{\tau }}^{\mathit{a}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{\odot }}{\mathrm{d}}^{\mathrm{2}}{{r}}^{\mathrm{\prime }}\mathrm{d}\mathit{z}\hspace{0.17em}{{K}}^{\mathit{a}}\mathrm{(}{{r}}^{\mathrm{\prime }}\mathrm{-}{r}\mathit{,z}\mathrm{)}\mathrm{·}{{v}}^{\mathrm{true}}\mathrm{\left(}{{r}}^{\mathrm{\prime }}\mathit{,z}\mathrm{\right)}\mathrm{+}{\mathit{n}}^{\mathit{a}}\mathrm{\left(}{r}\mathrm{\right)}\mathit{,}$(1)where the position vector splits into the horizontal r and the vertical z domains. Index a uniquely refers to the selection of the measurement geometry (there is a free choice of spatio-temporal filters, distance of the measurement points, and/or additional spatial averaging). The realisation of the random noise n^a is not known, however its covariance matrix may be measured from the data and used in the estimate of the random-noise level in the inverted flow velocities.

Equation (1) describes the forward problem, which gives the recipe for how to compute the (forward-modelled) travel times when the vector velocity field v is known. The usual requirement is an inverse modelling, hence the derivation of the velocity field in the Sun from the measured set of travel-time maps. To do so, various classes of methods were employed, where two of them are used most often: the regularised least squares (RLS) and optimally localised averaging (OLA). The RLS method (in time-distance helioseismology used for the first time by Kosovichev 1996) seeks to find the models of the solar interior, which provide the best least-squares fit to the measured travel-time maps, while regularising the solution (e.g. by requiring the smooth solution). OLA or its form Subtractive OLA (SOLA; Pijpers & Thompson 1992) is based on explicitly constructed spatially confined averaging kernels by taking a linear combination of sensitivity kernels, while simultaneously having control over the level of random noise in the results. A SOLA-type inversion is the principal method discussed in the current paper using a code validated by Švanda et al. (2011) and a data processing pipeline verified against the direct surface measurements by Švanda et al. (2013).

Both methods result in the estimates of flows at a given target depth. This estimate is a true velocity field smoothed by the averaging kernel and also contains the random-noise component (v^rnd). ${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}\mathrm{\left(}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{\odot }}\sum _{\mathit{\beta }}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathit{v}}_{\mathit{\beta }}^{\mathrm{true}}\mathrm{\left(}{r}\mathit{,z}\mathrm{\right)}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z}\mathrm{+}{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{rnd}}\mathrm{\left(}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathit{,}$(2)where α = (x,y,z) is a α-component of the velocity vector v^true. SOLA method uses a user-given target function, which provides an initial estimate of the resulting averaging kernel. The method balances the spatial misfit of the true averaging kernel, the required target function and the propagation of the random noise into the resulting tomographic maps.

The averaging kernel with components ${\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}$ can be derived from the inversion as a secondary product and requires normalisation so that its spatial integral equals unity. The components contain information about the smoothing of the flow component in the direction of the inversion (β = α) and also about the leakage (the cross-talk) from the other components (for β ≠ α).

The level of the random noise may also be estimated directly from the inversion (in case of OLA method), so at least its root-mean-square (RMS) value ${\mathit{\sigma }}_{\mathit{\alpha }}\mathrm{=}\mathrm{RMS}\mathrm{\left(}{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{rnd}}\mathrm{\right)}$ is known.

In effect, the output of SOLA inversion is a set of tomographic maps at various target depths. By the formulation, the results are therefore regularly gridded (defined on a regularly spaced grid) in the horizontal domain (where the spatial grid is usually regular), whereas in the vertical direction, the coverage of the domain is given by a selection of the target depths, hence usually not regularly gridded. The regularly gridded solutions are required for instance whenever the estimates of spatial derivatives are needed; for example, when computing for all components of divergence, vorticity, helicity or Reynolds stress. The aim of this study is to derive and test a method for estimating the regularly gridded 3D vector flow field from such a set of discrete tomographic maps with known averaging kernels and noise-level estimates. Simplified approaches were used recently, usually using the target depths as true representatives of the location of tomographic maps and interpolating between them in order to get an estimate of the regularly gridded velocity field (e.g. DeGrave & Jackiewicz 2015). We would like to point out that such simplified methods are often not justified, because the typical helioseismic averaging kernels contain sidelobes in their sensitivity and therefore it is difficult to assign only a single depth to such a kernel. The averaging kernels also often peak at a different depth than, for example, where their gravity centre is located, which makes the assignment of a single representative depth even more ambiguous. Our goal is to find a way of reconstructing a reasonable estimate of the vertical variations of the flow that is only sparsely covered by a set of tomographic maps. We would like to use the full structure of the known averaging kernels to properly reconstruct the estimate of the true flow. Our technique aims to combine the information from a set of inversions with different averaging kernels, thereby to some extent resembling the super-resolution (e.g. Hardie et al. 1997) technique, routinely used in the imaging problems.

2. Formulation of the problem

We seek the 3D regularly gridded vector velocity field which, if convolved with a set of given averaging kernels, provides a set of flow maps that is close to the set of tomographic maps from the inversion at given target depths. Inspired by Eq. (2), the sought relation may be written for each target depth z^d as $\frac{\mathrm{1}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}\mathrm{\left(}{r}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}\sum _{\mathit{\beta ,}{{r}}_{\mathit{j}}\mathit{,z}}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{(}{{r}}_{\mathit{j}}\mathrm{-}{r}\mathit{,z}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{)}{\mathit{v}}_{\mathit{\beta }}\mathrm{\left(}{{r}}_{\mathit{j}}\mathit{,z}\mathrm{\right)}\mathit{,}$(3)where we replaced the continuous integrations by their discrete limits: $\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}\hspace{0.17em}\mathrm{\to }\hspace{0.17em}\sum _{{{r}}_{\mathit{j}}}{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}\hspace{0.17em}\mathit{f}\mathrm{\left(}{{r}}_{\mathit{j}}\mathrm{\right)}\hspace{0.5em}\mathrm{and}\hspace{0.5em}\mathrm{\int }\mathrm{d}\mathit{z}\hspace{0.17em}\mathit{f}\mathrm{\left(}{z}\mathrm{\right)}\hspace{0.17em}\mathrm{\to }\hspace{0.17em}\sum _{{\mathit{z}}_{\mathit{j}}}{\mathit{h}}_{\mathit{z}}\hspace{0.17em}\mathit{f}\mathrm{\left(}{{z}}_{\mathit{j}}\mathrm{\right)}\mathit{.}$(4)Note that Eq. (3) is in fact a discretised Eq. (2) without a noise term explicitly given. Constant h_x is the spacing of the horizontal discrete grid (and h_z analogously represents the spacing in the vertical grid). For simplicity, we assume that the grids in both domains are regular and equidistant. For irregular grids our solution must either be modified accordingly, or all vectors, matrices and 3D arrays must be interpolated onto a regular grid. The symbolic summation over r_j represents the summation over all N = n_xn_y discrete values of position vector r_j having dimensions (n_x,n_y). Index d indicates the selection of a target depth. We assume that the velocity vector v is close to the true velocity field v^true. How close the two are is discussed in Sect. 2.2.

For each target depth z^d, there is a separate Eq. (3), however the modelled velocity field v is common to all these equations. The noise level of each of these maps very likely varies, hence we introduced an inverse weighting by the RMS of the random noise, basically meaning that less precise inferences (with larger RMS of random noise) are penalised against more precise inferences (with smaller RMS of random noise). For simplicity, we first discuss the equations to be solved for one target depth and only then combine all available target depths to form the larger problem.

Factors of $\mathrm{1}\mathit{/}\mathrm{(}{\mathit{\sigma }}_{\mathit{\gamma }}^{\mathit{d}}{\mathrm{)}}^{\mathrm{2}}$ in (3) could, in principle, be removed from each side of the equation. In the following, we use a set of Eq. (3) written for a set of target depths z^d, where keeping these factors in the equations simplifies an introduction of additional terms later in the paper, especially when a regularisation of the solution is studied in order to reduce the propagation of the random-noise.

Equation (3) may be written in a Fourier space with a definition of the Fourier transform given by a pair of equations; $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ \mathit{f}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}& \mathrm{=}& \\ \mathit{f}\mathrm{\left(}{r}\mathit{,z}\mathrm{\right)}& \mathrm{=}& \end{array}$where k is a horizontal wave vector and h_k spacing in the Fourier domain: $\frac{\mathrm{1}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\gamma }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}\stackrel{}{{{\mathit{v}}_{\mathit{\gamma }}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\gamma }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathit{h}}_{\mathit{z}}\sum _{\mathit{\beta }}\sum _{\mathit{z}}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\gamma }\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\beta }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}\mathit{,} \mathrm{\forall }{k}\mathit{.}$(7)The star ^∗ denotes the complex conjugate. Equation (7) holds for each wave vector k and target depth z^d. The precise derivation of Eq. (7) is given in Appendix A.

2.1. Solution for the flow

For each target depth z^d, we have three (three velocity components v_γ) Eqs. (7). However, all these equations have a common velocity field $\begin{array}{c}{}_{\mathrm{˜}}\\ {v}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}$ that is to be determined. Not counting the ideal or trivial cases, in reality we expect that it will be impossible to fulfill Eq. (7) for a set of different target depths z^d exactly. By casting the problem written for a set of target depths into the matrix equation, $\mathit{y}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}\mathit{A}\mathrm{\right(}{k}\mathrm{\left)}\mathit{x}\mathrm{\right(}{k}\mathrm{)}\mathit{,} \mathrm{\forall }{k}\mathit{,}$(8)the use of the pseudoinverse solvers allows us to find the solution x in a least-squares sense. Matrices y, A and x constitute of blocks of Eq. (7), namely $\mathit{y}\mathrm{=}\left(\begin{array}{c}\\ \left[\begin{array}{c}\\ \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{\mathit{x}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\stackrel{}{{{\mathit{v}}_{\mathit{x}}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}\\ \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{\mathit{y}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\stackrel{}{{{\mathit{v}}_{\mathit{y}}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}\\ \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\stackrel{}{{{\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}\end{array}\right]\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\mathrm{\downarrow }\mathrm{\left(}\mathit{d}\mathrm{\right)}\end{array}\right)\mathit{,}$(9)$\mathit{A}\mathrm{=}\left[\begin{array}{c}\\ \frac{\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{x}}^{\mathit{x}\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}}{{\left({\mathit{\sigma }}_{\mathit{x}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}& \frac{\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{y}}^{\mathit{x}\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}}{{\left({\mathit{\sigma }}_{\mathit{x}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}& \frac{\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{z}}^{\mathit{x}\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}}{{\left({\mathit{\sigma }}_{\mathit{x}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\\ \frac{\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{x}}^{\mathit{y}\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}}{{\left({\mathit{\sigma }}_{\mathit{y}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}& \frac{\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{y}}^{\mathit{y}\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}}{{\left({\mathit{\sigma }}_{\mathit{y}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}& \frac{\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{z}}^{\mathit{y}\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}}{{\left({\mathit{\sigma }}_{\mathit{y}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\\ \frac{\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{x}}^{\mathit{z}\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}}{{\left({\mathit{\sigma }}_{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}& \frac{\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{y}}^{\mathit{z}\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}}{{\left({\mathit{\sigma }}_{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}& \frac{\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{z}}^{\mathit{z}\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}}{{\left({\mathit{\sigma }}_{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\end{array}\right]\mathit{...}\mathrm{\to }\mathrm{\left(}\mathit{z}\mathrm{\right)}$(10)and $\mathit{x}\mathrm{=}\left(\begin{array}{c}\\ \left[\begin{array}{c}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{x}}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{y}}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{z}}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}\end{array}\right]\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\mathrm{\downarrow }\mathrm{\left(}\mathit{z}\mathrm{\right)}\end{array}\right)\mathit{.}$(11)The arrows in matrices y, A and b indicate the running indices in depth (z) and in target depths (d). The goal is to retrieve vector x, which contains the estimate of the regularly gridded velocity field. Such calculation may be done using standard solvers. We note that in cases where ${\mathrm{𝒦}}_{\mathit{\alpha }}^{\mathit{\beta }}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\mathrm{=}\mathrm{0}$ for α ≠ β, the solution will be simpler as the matrix A will consist of diagonal blocks. Such a solution is equivalent to the case when no cross-talk exists. Usually, there is a (possibly small) contribution of the cross-talk so it is wise to use such information in the reconstruction of the estimate of the true velocity field.

2.2. Averaging kernels

In the previous section, we assumed that the sought velocity field v is close to the real flow v^true. The relation “close” can be quantified by an assumed linear relation ${v}\mathrm{=}\mathrm{ℱ}{{v}}^{\mathrm{true}}\mathrm{+}{{v}}^{\mathrm{noise}}\mathit{,}$(12)where v^noise corresponds to the propagated random noise and ℱ is a linear operator representing the convolution with the new averaging kernels ${\mathrm{ℱ}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{t}}\mathrm{\right)}$. In this case, z_t is the target depth on the vertical grid. For each grid point on the vertical grid, we have a new averaging kernel. Equation (12) is therefore an equivalent of Eq. (2) in the original SOLA problem. Following the notation introduced in Eqs. (12), (2) may also be written in a symbolic operator form ${{v}}^{\mathrm{inv}}\mathrm{=}\mathrm{𝒦}{{v}}^{\mathrm{true}}\mathrm{+}{{v}}^{\mathrm{rnd}}\mathit{,}$(13)and (8) as ${v}\mathrm{=}{\mathrm{𝒜}}^{-1}{{v}}^{\mathrm{inv}}\mathit{,}$(14)where represents the operator form of matrix A (Eq. (10)) that also absorbed the multiplicative factor of (2π)²h_z for simplicity.

We insert (12) into (13) and use (14) to obtain $\mathrm{ℱ}{{v}}^{\mathrm{true}}\mathrm{+}{{v}}^{\mathrm{noise}}\mathrm{=}{\mathrm{𝒜}}^{-1}\mathrm{𝒦}{{v}}^{\mathrm{true}}\mathrm{+}{\mathrm{𝒜}}^{-1}{{v}}^{\mathrm{rnd}}\mathit{.}$(15)By comparing the terms we have $\mathrm{ℱ}\mathrm{=}{\mathrm{𝒜}}^{-1}\mathrm{𝒦} \mathrm{and} {{v}}^{\mathrm{noise}}\mathrm{=}{\mathrm{𝒜}}^{-1}{{v}}^{\mathrm{rnd}}\mathit{.}$(16)The equation above gives a recipe for how to compute the new set of averaging kernels from the old ones and explains how the realisation of the random noise propagates through the reconstruction. After performing the mathematics in detail, the new averaging kernels ${\mathrm{ℱ}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{t}}\mathrm{\right)}$ result from solving the equation ${\mathit{K}}^{\mathit{\alpha }}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}\mathit{A}\mathrm{\right(}{k}\mathrm{\left)}{\mathit{F}}^{\mathit{\alpha }}\mathrm{\right(}{k}\mathrm{)}\mathit{,} \mathrm{\forall }{k}\mathit{,}$(17)for each α, k, and z_t, where ${\mathit{K}}^{\mathit{\alpha }}\mathrm{=}\left(\begin{array}{c}\\ \left[\begin{array}{c}\\ \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{\mathit{x}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{x}}^{\mathit{\alpha }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}\\ \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{\mathit{y}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{y}}^{\mathit{\alpha }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}\\ \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{z}}^{\mathit{\alpha }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}\end{array}\right]\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\mathrm{\downarrow }\mathrm{\left(}\mathit{d}\mathrm{\right)}\end{array}\right)\mathit{,}$(18)

${\mathit{F}}^{\mathit{\alpha }}\mathrm{=}\left(\begin{array}{c}\\ \left[\begin{array}{c}\\ \begin{array}{c}\mathrm{˜}\\ {\mathrm{ℱ}}_{\mathit{x}}^{\mathit{\alpha }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{t}}\mathrm{\right)}\\ \begin{array}{c}\mathrm{˜}\\ {\mathrm{ℱ}}_{\mathit{y}}^{\mathit{\alpha }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{t}}\mathrm{\right)}\\ \begin{array}{c}\mathrm{˜}\\ {\mathrm{ℱ}}_{\mathit{z}}^{\mathit{\alpha }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{t}}\mathrm{\right)}\end{array}\right]\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\mathrm{\downarrow }\mathrm{\left(}\mathit{z}\mathrm{\right)}\end{array}\right)\mathit{,}$(19)and A is given by Eq. (10).

3. Random noise

Tomographic maps from time-distance helioseismology naturally contain some level of random noise, which originates from the random excitation of the waves by convection. The realisation of this noise propagates through our reconstruction procedure into the regularly gridded estimate of the velocity field.

3.1. Random noise-level estimation

Similarly to the determination of the averaging kernel, we could proceed in deriving the propagation of random noise through the procedure, basically leading to the second part of Eq. (16). The tomographic maps as inputs to our reconstruction have known estimates of random noise-levels in the form of RMS of this random component only. As we show later, this might not be enough, as we also need (and use) the estimate of the spatial power spectrum of this random-noise component N(k) (for details see Appendix B). In such a case, the use of the simulated random-noise realisation seems to be a feasible solution to obtain an estimate of the random noise in the reconstructed flow cubes.

The noise part of (16) then written explicitly transforms into $\mathit{\delta y}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}\mathit{A}\mathrm{\right(}{k}\mathrm{\left)}\mathit{\delta x}\mathrm{\right(}{k}\mathrm{)}\mathit{,} \mathrm{\forall }{k}\mathit{,}$(20)where A is given by (10) and for δx and δy explicitly $\mathit{\delta y}\mathrm{=}\left(\begin{array}{c}\\ \left[\begin{array}{c}\\ \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{\mathit{x}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\stackrel{}{{{\mathit{v}}_{\mathit{x}}^{\mathrm{rnd}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}\\ \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{\mathit{y}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\stackrel{}{{{\mathit{v}}_{\mathit{y}}^{\mathrm{rnd}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}\\ \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\right)}^{\mathrm{2}}}\stackrel{}{{{\mathit{v}}_{\mathit{z}}^{\mathrm{rnd}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}^{{\mathit{d}}_{\mathrm{1}}}\mathrm{\right)}\end{array}\right]\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\mathrm{\downarrow }\mathrm{\left(}\mathit{d}\mathrm{\right)}\end{array}\right)\mathit{,} \mathrm{and} \mathit{\delta x}\mathrm{=}\left(\begin{array}{c}\\ \left[\begin{array}{c}\\ \stackrel{}{{{\mathit{v}}_{\mathit{x}}^{\mathrm{noise}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}\\ \stackrel{}{{{\mathit{v}}_{\mathit{y}}^{\mathrm{noise}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}\\ \stackrel{}{{{\mathit{v}}_{\mathit{z}}^{\mathrm{noise}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}\end{array}\right]\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\mathrm{\downarrow }\mathrm{\left(}\mathit{z}\mathrm{\right)}\end{array}\right)\mathit{.}$(21)Values of ${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{rnd}}$ are randomly generated so that their spatial power spectrum is N(k) and their RMS equals ${\mathit{\sigma }}_{\mathit{\alpha }}^{{\mathit{z}}_{\mathrm{0}}}$. The estimate of the random-noise level ${\mathit{\sigma }}_{\mathit{\alpha }}^{{\mathit{z}}_{\mathrm{t}}}$ propagating to the resulting reconstructed velocity field is then computed as ${\mathit{\sigma }}_{\mathit{\alpha }}^{{\mathit{z}}_{\mathrm{t}}}\mathrm{=}\mathrm{RMS}\mathrm{\left[}{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{noise}}\mathrm{\right(}{r}\mathrm{;}{\mathit{z}}_{\mathrm{t}}{\mathrm{\left)}}^{\mathrm{\right]}}\mathit{.}$(22)

3.2. Dealing with random noise

The algorithm described in Sect. 2 actually belongs to the problems of deconvolution in astrophysics, which are known to be unstable in the presence of random noise. Using our synthetic tests, we found that our implementation, as described in the previous Sections in the case of realistic noise power spectrum, boosts the random noise by approximately eleven orders of magnitude, which is unacceptable for any reasonable flow reconstruction.

Thus, we require a method for suppressing the propagation of the random noise. A way out is to avoid the spatial wave numbers, where the noise dominates. This can either be achieved by an introduction of the filter in the k-space domain that, in effect, allows us to reconstruct for only some subset of wave vectors and to replace the contribution of wave vectors with large wave numbers (where the noise is likely do dominate) by zeros. We found that for a small level of noise (less than 1% contribution) when approximately 20% smallest wave numbers were used, the reconstruction led to an acceptable solution. For any larger noise contribution, the solution was not acceptable.

Another option is to introduce an ad-hoc regularisation term to ensure a smooth solution, by assuming that the random noise will cause small-scale oscillations in the solution. Such regularisation is relatively common in helioseismology. We may introduce a smoothing parameter ϵ and, thanks to keeping $\mathrm{1}\mathit{/}\mathrm{(}{\mathit{\sigma }}_{\mathit{\gamma }}^{\mathit{d}}{\mathrm{)}}^{\mathrm{2}}$ factors in Eq. (7), simply modify the matrix A to A′: ${\mathit{A}}^{\mathrm{\prime }}\mathrm{=}\mathit{A}\mathrm{+}\left(\begin{array}{c}\\ \left[\begin{array}{c}\\ \mathit{ϵ}{\mathit{k}}^{\mathrm{2}}\mathrm{0}\mathrm{0}\\ \mathrm{0}\mathit{ϵ}{\mathit{k}}^{\mathrm{2}}\mathrm{0}\\ \mathrm{0}\mathrm{0}\mathit{ϵ}{\mathit{k}}^{\mathrm{2}}\end{array}\right]& \mathit{...}\mathrm{\to }\mathrm{\left(}\mathit{z}\mathrm{\right)}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\mathrm{\downarrow }\mathrm{\left(}\mathit{d}\mathrm{\right)}\end{array}\right)\mathit{.}$(23)Our ad-hoc term resembles that of the RLS inversion, as stated by Couvidat et al. (2006), for example. Then, we invert for x from equation y = (2π)²h_zA′x. The boosting of the noise by eleven orders mentioned above forces us to use a very strong regularisation (with ϵ in the order of 10⁴ or more), which, in effect, allows us to reconstruct only very large-scale components of the spatial spectrum correctly, making the results extremely smooth. All the details of the flow field are lost. We do not find such a result satisfactory. However, the regularisation may become useful, when dealing with numerical stability, for example, when the matrix A is close to singular.

As a final attempt, we adapt solutions from the image restoration problems, namely the Wiener filtering. Wiener filtering introduces something resembling a k-space filter that allows us to weight wave numbers so that only the wave numbers where the signal dominates the random noise contribute to the reconstructed image. The method requires prior knowledge about the expected power spectrum of the real solar flows and random noise, which we have at our disposal (details in Appendix B).

Inspired by the Wiener-filtering algorithm instead of inverting for x from y = (2π)²h_zAx for each k, we obtain an optimal estimate of x as $\mathit{x̂}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}}\mathit{Gy,} \mathrm{where} \mathit{G}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}{\mathit{A}}^{-1}\mathrm{\left(}{k}\mathrm{\right)}\left[\frac{\mathrm{\left|}\mathrm{\right|}\mathit{A}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{\left|}\mathrm{\right|}{\mathit{h}}_{\mathit{k}}^{\mathrm{2}}}{\mathrm{\left|}\mathrm{\right|}\mathit{A}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{\left|}\mathrm{\right|}{\mathit{h}}_{\mathit{k}}^{\mathrm{2}}\mathrm{+}\frac{\mathit{N}\mathrm{\left(}{k}\mathrm{\right)}}{\mathit{S}\mathrm{\left(}{k}\mathrm{\right)}}}\right]\mathrm{·}$(24)Note that A(k) is a matrix, whereas in the original Wiener formulation it was a scalar. We represent the norm of the matrix by the Euclidean norm. N(k) and S(k) represent the estimates of power spectra of the random noise and the signal (for details see Appendix B).

4. Implementation and synthetic tests

We implemented the above given algorithm in Matlab programming language. To assess the performance of the method, we constructed synthetic inversion-like tomographic maps obtained by smoothing the numerical simulation of the solar convection (Fig. 1, source data from Ustyugov 2008) with a set of averaging kernels (Fig. 2). The synthetic averaging kernels were constructed as Gaussians in both directions with a horizontal width of 5 Mm. In the vertical direction, the averaging kernels target depths of 1, 1.9, 2.9, 4.3, 6.2, 9.2 and 12.1 Mm with a full-width-at-half-maximum equal to the target depth of the given kernel.

Fig. 1 Simulated convective flows used for the validation of our reconstruction techniques. The units of the coordinates are in Mm.

Fig. 2 Depth sensitivity of averaging kernels used in construction of seven synthetic tomographic maps.

The reconstruction was performed on a horizontal grid identical to the horizontal grid of the model (with a pixel size of 1.4 Mm) and the vertical grid extending from 10 Mm depth to 400 km above the surface with a vertical sampling of 200 km. For each grid point, we computed the estimate of the reconstructed flow vector, averaging kernels for each flow component, and also the estimate of the noise levels for each flow component. Due to the assumed translation invariance in the horizontal domain, for each flow component, the averaging kernels and the estimates of the noise levels are identical for all horizontal points within the identical depth.

Fig. 3 Examples of original averaging kernels (dashed) and averaging kernels after deconvolution at the same depths (solid). Only horizontally averaged kernels are plotted.

Fig. 4 Comparison of horizontal cuts through averaging kernels (along y = 0) at the target depth 1.9 Mm for the input averaging kernel (black) and the reconstructed one (red).

The averaging kernels ${\mathrm{ℱ}}_{\mathit{\beta }}^{\mathit{\alpha }}$ are, in general, narrower than the averaging kernels ${\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}$ of the input set of tomographic maps (see Fig. 3), which is what we wanted. In Fig. 3, we only compared averaging kernels at corresponding target depth; we remind the reader that due to the sampling of the vertical domain as described above, there is, in fact, 53 averaging kernels in total, one for each grid point in the vertical domain. In our testing case, the overall shape of the new averaging kernels ℱ is not satisfying compared to the original one, as the averaging kernels of the reconstruction contain many sidelobes (the negative ones as well). It must be noted that in a realistic case, the averaging kernels also contain sidelobes (the negative ones as well). Therefore, from this point of view, it is not possible to immediately say which is “worse”.

In Fig. 3, we plotted examples of horizontally integrated kernels for simplicity. The horizontally integrated kernels are independent of noise levels, as by using the Wiener-like filtering, the k = 0 is always reconstructed. In the horizontal domain, the shape of the kernels is influenced by the level of random noise. In the ideal case (no noise present), the deconvolution goes to approximately the level of a pixel size (see Fig. 4); they are close to a numerical Dirac δ function. The described behaviour is exactly the same for all depths as the input averaging kernels had the same full-width-at-half-maximum of 5 Mm. With an increasing level of noise, the kernels get wider and they approach the shape of the input kernels for larger noise levels. Then, in effect, only a little or no deconvolution of the flow in the horizontal domain takes place. We would like to stress that our main goal was to reconstruct the reasonable estimate of the flow in the vertical domain, where the input set of tomographic maps only sparsely covers the vertical domain.

The total integrals of the new averaging kernels are then used to properly scale the magnitude of the flows at each target depth. This re-normalisation is necessary because Wiener-like filtering suppresses some of the Fourier components and therefore naturally introduces a scaling of resulting flow magnitudes. Consequently, we normalise new averaging kernels so that their spatial integral equals unity.

The reconstructed flow compared to the model is shown in Figs. 5–8. One can see that the method works extremely well in cases where no noise is present. The correlation coefficient between Figs. 5a and b is 0.8–0.9 for all components. At first glance, almost all details of the flow field seem to be correctly reconstructed (see also first two rows of Fig. 6). Also, the depth structure is reconstructed very well, which can be seen in Fig. 8 (black and red lines). Small-scale oscillations of the flow are well beyond the resolution limit set by the averaging kernels, but the smoothed trend is captured relatively well. Interestingly, the two cuts (black and red lines of Fig. 8) match one another better in the depth range of 0–6 Mm than deeper down. At depths 0–6 Mm, all the averaging kernels overlap.

We further investigated this issue and found that the best reconstruction is indeed achieved in depth range, which is “included” in several tomographic maps. This gives a recipe for how to select the target function for SOLA inversions that may potentially lead to a reasonable estimate of the regularly gridded velocity field by our method. The target functions should be relatively broad and should overlap. Also, a larger set of different target functions provide a more robust estimate of the true velocity field.

Our reconstruction method also returns reasonable results in the case where the random noise is present in the input data. The fact that the Wiener-like approach out-performs the ad-hoc regularisation is clear from Fig. 7, where we compare example results of reconstruction methods differing by noise treatment. Obviously, without taking care of the noise, the results are completely wrong. The reconstructed flow magnitude oscillates from pixel to pixel and the noise is boosted by almost twelve orders. Adding an ad-hoc smoothing term improves the performance, however very large values of the smoothing parameter are needed, and even in a such case, the results are not satisfactory; they either still contain a large fraction of propagated random noise, or are very smooth with most of the details of the true velocity field lost. Wiener-like approach provides an optimal solution to the problem. Therefore we use this method further and study its properties.

A Wiener-like approach implemented by us allows a reasonable reconstruction of the flows for cases where the signal-to-noise ratio (S/N) is larger than 20 and positive correlation is found even with a S/N of 10 (meaning that the inverted flow maps are polluted with 10% random noise – see Table 1). S/N values may be understood by following the model of inversion for supergranular flows. The typical magnitude of the horizontal flow components is approximately 300 m s^-1, reasonable inversions have the random-noise level of 25 m s^-1. In such a case, S/N equals to 12. A S/N of 100 would approximately represent the flow maps obtained after averaging of 100 supergranule individuals.

The level of random noise is estimated for each depth following the procedure described in Sect. 3.1. As an example, we plotted the estimated RMS of the random noise in the reconstructed horizontal flow in Fig. 9 together with the RMS of the horizontal flow at the same depths for the case where the S/N is 100. At all depths, the S/N ratio of the reconstructed flow is above unity and at both ends of the domain, S/N approaches unity. While at the bottom boundary this is due to the fact that the reconstruction is dominated by noise, at the top boundary, the increase of both the noise RMS and the flow amplitude is an artifact of the normalisation by the new averaging kernel. As seen in Fig. 2, at the top of the domain, there is little power in all original averaging kernels , hence the averaging kernels ℱ at these depths have an arbitrary shape with little localisation. As a consequence, the reconstruction of the flow at these depths is physically meaningless.

It should be noted that the reconstruction behaves reasonably in cases of large signal-to-noise ratio. Interestingly, this requirement conforms well with conclusions described at the beginning of this section, where we showed that rather broad averaging kernels lead to a better reconstruction of the regularly gridded estimate of the flow. In SOLA, localisation in the Sun (essentially the width of the averaging kernel) and the level of random noise are balanced by a trade-off parameter in the inversion cost function. Thus, the choice of a broad target function automatically leads to a lower level of random noise and hence a larger S/N for a given tomographic map. The selection of a set of target functions that overlap prepares a set of tomographic maps suitable for the reconstruction of the regularly gridded estimate of a true velocity field.

Fig. 5 Effect of the random noise to the reconstruction of the regularly gridded velocity field: the case of the near-surface flows. Displayed are results with varying noise level: a) input data; b) no noise. Other maps were reconstructed with Wiener-like filtering to deal with random-noise boosting with varying S/N levels: c) S/N 1000; d) S/N 100; e) S/N 20 and f) S/N 10. In each cut, the parallel components are indicated by arrows, while the perpendicular component is indicated by the colour. Colour-scale of all panels is identical. The units of the coordinates are in Mm. The maps of all three components at the depth of 1 Mm are displayed in Fig. 6.

Fig. 6 Maps of all components of the flow at the depth of 1 Mm for various levels of the random noise. Maps for the x component, the y component and the vertical z component of velocity are found in the left-hand, the middle and the right-hand column, respectively. The rows contain the input velocity field in the first row and reconstructions with increasing level of random noise in order from top to bottom. Panels with given S/N value were reconstructed using a Wiener-like approach. The correlation coefficient for the whole cube between the input velocity field and the reconstruction is given in the titles of the plot panels. Colour-bars have the units of km s-1 and the field-of-view shown here covers 38 Mm × 38 Mm in the horizontal directions.

Fig. 7 Testing approaches to deal with random noise. An input datacube with a S/N of 100 (i.e. 1% random noise contribution – first row) was reconstructed with the basic method ignoring the noise (second row), ad-hoc regularisation with a set of values of regularisation parameter (third to fifth row) and the Wiener-like filtering approach (last row). The correlation coefficient for the whole cube between the input velocity field and the reconstruction is given in the titles of the plot panels. Colour-bars have the units of km s-1 (violet indicates fast variations from pixel to pixel) and the field-of-view shown here covers 38 Mm × 38 Mm in the horizontal directions.

Fig. 8 Cut through the horizontal flow component at one point for the input velocity field (black) and reconstructions with varying noise levels (coloured lines). The shaded regions indicate the respective error bars. Note that the black and red lines do not have errors.

Even the Wiener-like solution is not ideal. The prior knowledge (or an estimate) of a spatial spectrum of the velocities and of the random noise poses a real challenge. It may be estimated reasonably for the surface, for example, in the way outlined in this paper, however, it is expected from the theory that the spatial spectrum of the convective flows changes with depth. Such variability is very difficult to incorporate in the reconstruction.

In the synthetic test described above, we used the model averaging kernels (Fig. 2), which do not have a cross-talk component, whereas matrix A (Eq. (10)) allows us to use the information even from the cross-talk (off-diagonal) terms. The presence of cross-talk is an issue for inversions for weak quantities, such as the vertical component of the flow. In order to estimate the robustness of our method in the presence of cross-talk, we performed additional tests. Firstly, the averaging kernels were computed in the same way as before. Additionally, we added cross-talk components that have the same shape as the components in the direction of the inversion. Such a situation describes the case when the cross-talk is highly correlated with a desired signal and its contribution is twice larger than the contribution of the flow component we invert for. Despite the amount of large cross-talk, the reconstructed flow resembles the model one with a correlation coefficient of ~0.4 for the horizontal components and ~ 0.2 for the vertical component for the case of S/N = 0.05. Secondly, we kept only the cross-talk components and set the averaging kernels in the direction of the inversion to zero. The flow was recovered (with a correlation coefficient of ~ 0.5) with the model for both vertical and horizontal components.

Additional tests showed that the correlation coefficient between the recovered flow and the model is lower when the cross-talk components of the averaging kernels are similar and comparable in peak magnitude to the averaging kernel in the direction of the inversion. A significant difference in peak magnitudes or an opposite sign actually help our method to improve the trustworthiness of the resulting reconstruction. The performed tests point out the difficulties caused by the presence of cross-talk, however they also show, via these poorly representative examples, that our reconstruction method is still able to recover the true structure of the flow reasonably well.

5. Summary

We derived a method for deconvolution of a set of tomographic maps coming from local helioseismology, which allows for construction of a regularly gridded estimate of the 3D vector velocity field. The method was successfully tested using synthetic data. We showed that it provides reasonable estimates of the true velocity field with a lesser smoothing than the original set of tomographic maps. To deal with random noise, we introduced an approach similar to Wiener filtering known from image processing.

The method can naturally be applied to real data coming from SOLA inversions, which will be the topic of future study. It can be applied to any set of tomographic maps with a reasonable signal-to-noise ratio (larger than 50). The best usability can be expected in studies such as the investigation of the flow structure of solar supergranulation or similar features of convective flows.

Fig. 9 Level of random noise in the reconstructed flow (blue) and the RMS of the reconstructed flow (red). Both curves are plotted for the horizontal component of velocity and for the case with a S/N of 100. Obviously, the error increases beyond reasons at the ends of the domain.

Acknowledgments

This work benefits from the student grant awarded in 2014 to M.K. by the Faculty of Mathematics and Physics, Charles University in Prague, which was supervised by M.Š. M.Š. acknowledges the support of the Czech Science Foundation (grant 14-04338S) and of the institute research project RVO:67985815 to Astronomical Institute of Czech Academy of Sciences. The help of Dominik Jančík and Sven Ubik with the visualisation of the flows is greatly acknowledged.

References

Appendix A: Derivation of Eq. (7)

In Sect. 2, we prescribed an ad-hoc set of Eqs. (7) to be solved to obtain an estimate of the reconstructed 3D vector flow defined on a regular grid. In the following, we show that these equations may be properly obtained by a minimisation procedure.

We define a misfit between the inverted flow map at target depth z^d and the modelled tomographic map as $\begin{array}{ccc}{\mathit{\chi }}_{{\mathit{z}}^{\mathit{d}}}^{\mathrm{2}}\mathrm{\left(}{r}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{\alpha }}{\left[{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}\mathrm{\left(}{r}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\mathrm{-}{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}\sum _{\mathit{\beta ,}{{r}}_{\mathit{j}}\mathit{,z}}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{(}{{r}}_{\mathit{j}}\mathrm{-}{r}\mathit{,z}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{)}{\mathit{v}}_{\mathit{\beta }}\mathrm{\left(}{{r}}_{\mathit{j}}\mathit{,z}\mathrm{\right)}\right]}^{\mathrm{2}}\\ & & \end{array}$(A.1)For each target depth z^d and positional vector r, there is a separate Eq. (A.1), however the modelled velocity field v is common to all these equations. Note the inverse weighting by the RMS of the random noise already discussed in Sect. 2.

Equation (A.1) can be seen as a cost function, which we want to minimise with respect to the unknown velocity field v. Using standard methods of calculus we first obtain equation $\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \frac{\mathit{\partial }{\mathit{\chi }}_{{\mathit{z}}^{\mathit{d}}}^{\mathrm{2}}\mathrm{\left(}{r}\mathrm{\right)}}{\mathit{\partial }{\mathit{v}}_{\mathit{\gamma }}\mathrm{\left(}{{r}}_{\mathrm{1}}\mathit{,}{\mathit{z}}_{\mathrm{1}}\mathrm{\right)}}\mathrm{=}\sum _{\mathit{\alpha }}\frac{\mathrm{2}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}\mathrm{(}\mathrm{-}{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}{\mathrm{𝒦}}_{\mathit{\gamma }}^{\mathit{\alpha }}\mathrm{(}{{r}}_{\mathrm{1}}\mathrm{-}{r}\mathit{,}{\mathit{z}}_{\mathrm{1}}\mathit{,}{\mathit{z}}^{\mathit{d}}{\mathrm{)}}^{\mathrm{)}}\\ & & \mathrm{\times }\left[{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}\mathrm{\left(}{r}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\mathrm{-}{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}\sum _{\mathit{\beta }}\sum _{{{r}}_{\mathit{j}}}\sum _{\mathit{z}}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{(}{{r}}_{\mathit{j}}\mathrm{-}{r}\mathit{,z}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{)}{\mathit{v}}_{\mathit{\beta }}\mathrm{\left(}{{r}}_{\mathit{j}}\mathit{,z}\mathrm{\right)}\right]\mathit{,}\end{array}$(A.2)which holds for each z^d and γ. The minimisation is performed with respect to the γ-component of velocity at point (r₁,z₁). Since in each of these points, the equation equals zero, we may sum over all of these points in the computational domain. For each γ, we obtain $\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \sum _{{{r}}_{\mathrm{1}}}\sum _{{\mathit{z}}_{\mathrm{1}}}\sum _{\mathit{\alpha }}\frac{\mathrm{2}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}\mathrm{(}\mathrm{-}{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}{\mathrm{𝒦}}_{\mathit{\gamma }}^{\mathit{\alpha }}\mathrm{(}{{r}}_{\mathrm{1}}\mathrm{-}{r}\mathit{,}{\mathit{z}}_{\mathrm{1}}\mathit{,}{\mathit{z}}^{\mathit{d}}{\mathrm{)}}^{\mathrm{)}}\\ & & \end{array}$(A.3)We note that the first term reiterates the normalisation condition to the averaging kernels $\sum _{{\mathit{z}}_{\mathrm{1}}}\sum _{{{r}}_{\mathrm{1}}}{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}{\mathrm{𝒦}}_{\mathit{\gamma }}^{\mathit{\alpha }}\mathrm{\left(}{{r}}_{\mathrm{1}}\mathit{,}{\mathit{z}}_{\mathrm{1}}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}_{\mathit{\gamma }}^{\mathit{\alpha }}\mathit{,}$(A.4)where ${\mathit{\delta }}_{\mathit{\gamma }}^{\mathit{\alpha }}$ is a Kronecker delta. We solve (A.3) in the Fourier space using a definition given by the pair of Eqs. (5) and (6) to obtain $\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\gamma }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathit{h}}_{\mathit{k}}^{\mathrm{2}}\sum _{{k}}\stackrel{}{{{\mathit{v}}_{\mathit{\gamma }}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\exp \left[\mathrm{i}{k}\mathrm{·}{r}\right]\\ & & \mathrm{+}\frac{\mathrm{1}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\gamma }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}{\mathit{h}}_{\mathit{z}}{\mathit{h}}_{\mathit{k}}^{\mathrm{4}}\sum _{\mathit{\beta }}\sum _{{{r}}_{\mathit{j}}}\sum _{\mathit{z}}\sum _{{{k}}_{\mathrm{1}}}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\gamma }}\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathit{,z}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\\ & & \end{array}$(A.5)Equation (A.5) holds for each z^d, each γ and each r. In the following, we use the orthogonality of the functions exp[ik·r] and solve (A.5) in the k-by-k manner. After some algebra, we obtain $\frac{\mathrm{1}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\gamma }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}\stackrel{}{{{\mathit{v}}_{\mathit{\gamma }}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}{{\mathrm{(}{{\mathit{\sigma }}_{\mathit{\gamma }}^{\mathit{d}}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathit{h}}_{\mathit{z}}\sum _{\mathit{\beta }}\sum _{\mathit{z}}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\gamma }\mathrm{\ast }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}^{\mathit{d}}\mathrm{\right)}\begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\beta }}\end{array}\mathrm{\left(}{k}\mathit{,z}\mathrm{\right)}\mathit{,} \mathrm{\forall }{k}\mathit{,}$(A.6)which is Eq. (7).

Appendix B: Estimates of spatial power spectra

In SOLA, the cost function (Eq. (12) in Švanda et al. 2011) allows us to balance between the localisation (the averaging kernel) and the random-noise level in the resulting flow map. Symbolically written $\mathrm{cost}\mathrm{function}\mathrm{=}\mathrm{avg}\mathit{.}\mathrm{kernel}\mathrm{misfit}\mathrm{+}\mathit{\mu }\hspace{0.17em}\mathrm{random}\mathrm{noise}\mathit{,}$(B.1)where misfit evaluates the localisation of the averaging kernel and μ is a trade-off parameter free of choice. For details, we refer the reader to Pijpers & Thompson (1992). A minimisation of this cost function gives, at the end of the whole procedure, a flow map at a given depth together with an averaging kernel and an estimate for a level of random noise. By choosing values of the trade-off parameter μ, one can easily influence the random-noise level in these maps.

Fig. B.1 Estimates of power spectra of the signal (blue) and the noise (red) obtained from the real data as a function of angular degree l. The power spectrum of the signal (flows) was normalised so that the RMS of the flow is 100 m s-1. The power spectrum of the random noise was also normalised to RMS of 10 m s-1. Such a situation therefore represents S/N = 10.

From the resulting flow map, one can compute a spatial power spectrum P(k) using standard methods. Let us consider two extreme cases. Firstly, let us choose a large μ, so that the random-noise term multiplied by μ dominates the right-hand side of inversion cost function (B.1). This selection causes the random-noise term to be preferentially minimised and the resulting flow maps contain mostly the signal. The resulting spatial power spectrum P(k) approximates the power spectrum of the flows; let us name it S(k). The drawback of this selection is usually a worse localisation in the Sun (a worse-quality averaging kernels with sidelobes).

As a second extreme, let us select a small μ. Then, the misfit term in the cost function (B.1) is preferred and the resulting flow map contains a large fraction of noise. It may even be noise dominated, meaning the level of random noise is much larger (by many orders) than the expected magnitude of the real flows. The power spectrum P(k) approximates the power spectrum of the noise, named N(k). The localisation in the Sun is usually good in this case.

The estimates for the power spectrum of the flows S and the noise N would not be robust enough if only inversion for a single flow map was considered. Also, in the case of a single map, the difference between the averaging kernels for the two kinds of inversions could be large, meaning that the validity of our approach to obtain the power-spectra estimates might be questioned. To improve the robustness, one needs to average both power spectra overa larger set of flow maps. We used 1100 individual flow maps obtained with both signal- and noise-dominated inversions to estimate power spectra S and N. Averaging over a large sample of flow maps also allows us to tune values of needed μ so that the difference (evaluated by eye) of the resulting averaging kernels for noise- and signal-dominated inversion is not too large.

We believe that the methodology described here allows us to derive robust estimates for spatial power spectra of the flows S and the noise N. Both power spectra might then be easily scaled so that their RMS values attain desired values by simply invoking the Parseval theorem and the linearity between the RMS and the scaling of the random sample.

Examples can be seen in Fig. B.1 for a S/N of 10. In such a situation, obviously the noise dominates the signal at large angular degrees; the situation is largely reversed at supergranular scales of l ~ 120 and then again in the low-l part of the spectrum.

All Tables

All Figures

	Fig. 1 Simulated convective flows used for the validation of our reconstruction techniques. The units of the coordinates are in Mm.
In the text

	Fig. 2 Depth sensitivity of averaging kernels used in construction of seven synthetic tomographic maps.
In the text

	Fig. 3 Examples of original averaging kernels (dashed) and averaging kernels after deconvolution at the same depths (solid). Only horizontally averaged kernels are plotted.
In the text

	Fig. 4 Comparison of horizontal cuts through averaging kernels (along y = 0) at the target depth 1.9 Mm for the input averaging kernel (black) and the reconstructed one (red).
In the text

Fig. 5 Effect of the random noise to the reconstruction of the regularly gridded velocity field: the case of the near-surface flows. Displayed are results with varying noise level: a) input data; b) no noise. Other maps were reconstructed with Wiener-like filtering to deal with random-noise boosting with varying S/N levels: c) S/N 1000; d) S/N 100; e) S/N 20 and f) S/N 10. In each cut, the parallel components are indicated by arrows, while the perpendicular component is indicated by the colour. Colour-scale of all panels is identical. The units of the coordinates are in Mm. The maps of all three components at the depth of 1 Mm are displayed in Fig. 6.

In the text

Fig. 6 Maps of all components of the flow at the depth of 1 Mm for various levels of the random noise. Maps for the x component, the y component and the vertical z component of velocity are found in the left-hand, the middle and the right-hand column, respectively. The rows contain the input velocity field in the first row and reconstructions with increasing level of random noise in order from top to bottom. Panels with given S/N value were reconstructed using a Wiener-like approach. The correlation coefficient for the whole cube between the input velocity field and the reconstruction is given in the titles of the plot panels. Colour-bars have the units of km s-1 and the field-of-view shown here covers 38 Mm × 38 Mm in the horizontal directions.

In the text

Fig. 7 Testing approaches to deal with random noise. An input datacube with a S/N of 100 (i.e. 1% random noise contribution – first row) was reconstructed with the basic method ignoring the noise (second row), ad-hoc regularisation with a set of values of regularisation parameter (third to fifth row) and the Wiener-like filtering approach (last row). The correlation coefficient for the whole cube between the input velocity field and the reconstruction is given in the titles of the plot panels. Colour-bars have the units of km s-1 (violet indicates fast variations from pixel to pixel) and the field-of-view shown here covers 38 Mm × 38 Mm in the horizontal directions.

In the text

	Fig. 8 Cut through the horizontal flow component at one point for the input velocity field (black) and reconstructions with varying noise levels (coloured lines). The shaded regions indicate the respective error bars. Note that the black and red lines do not have errors.
In the text

	Fig. 9 Level of random noise in the reconstructed flow (blue) and the RMS of the reconstructed flow (red). Both curves are plotted for the horizontal component of velocity and for the case with a S/N of 100. Obviously, the error increases beyond reasons at the ends of the domain.
In the text

	Fig. B.1 Estimates of power spectra of the signal (blue) and the noise (red) obtained from the real data as a function of angular degree l. The power spectrum of the signal (flows) was normalised so that the RMS of the flow is 100 m s-1. The power spectrum of the random noise was also normalised to RMS of 10 m s-1. Such a situation therefore represents S/N = 10.
In the text