Free Access
Issue
A&A
Volume 575, March 2015
Article Number A122
Number of page(s) 8
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201425203
Published online 09 March 2015

© ESO, 2015

1. Motivation

Helioseismology is considered to be a state-of-the-art method of solar research to learn about the solar interior (see review by Gizon et al. 2010). Among various helioseismic methods, the local helioseismology plays an important role in studying the spatially localised details. The time–distance helioseismology (Duvall et al. 1993) is one of the approaches. It consists of tools for measuring and interpretation of travel times of solar waves. Since its formulation, it has been used extensively, especially to infer information about flows and sound-speed anomalies in the upper layers of solar convection zone.

The method assumes that the travel time τa is bound to the plasma anomalies δqα with respect to the background by a linear equation (1)where is the sensitivity kernel coming from the forward modelling, r = (x,y) is a horizontal position vector in a Cartesian coordinate system (the remaining vertical component will be denoted z), and na is the travel-time random noise. Index a denotes the selection of the waves (the combination of the spatio-temporal filter, spatial averaging, and the distance between the travel-time measurement points) and indices α or β select the perturbation δq (flow components, sound speed perturbation, etc.).

One of the main goals of time–distance helioseismology is to infer the structure of convective flows in the near-surface layers of solar convection zone from surface measurements of wave travel times. Formally it means that we replace δqβ in (1) by vβ, where β = (x,y,z), and search for it. This can only be achieved by inverse modelling.

To solve the inverse problem, one has to set up the cost function with various terms, which can be cast to a system of linear equations. How to perform such task has been described in detail and discussed elsewhere. The resulting flow estimate v is given by (2)where are the inversion weights to be determined and indicates the component of the inversion averaging kernel (a linear combination of the sensitivity kernels with inversion weights). The averaging kernel describes the smoothing of the real solar convective flows and is constructed to peak around a target depth z0. The second term on the right-hand side of Eq. (2) indicates the realisation of the random noise.

For the inversion method which is the basis of this study, the reader is referred to Jackiewicz et al. (2012) and Švanda et al. (2011). The cost function contains four terms balanced by three trade-off parameters: the quality of the fit (the misfit) of the averaging kernel to a user-defined target function localised in the Sun, the level of the random noise, the level of the pollution of the inverted flow component by other components (the cross-talk), and the term that ensures that the resulting inversion weights are spatially localised (which is a requirement in order to fulfil the mathematical assumptions).

One of the biggest issues in time–distance helioseismology is the presence of the realisation noise. Pressure and surface gravity waves are randomly excited by the vigorous convection and hence their statistics inherently contain a large random component. This acts as a random noise in any of the helioseismic observables and travel times of the waves are no exception. Various helioseismic methods deal with the random noise differently. Sometimes the noise is ignored (e.g. Kosovichev & Duvall 1997), sometimes the variance of the travel times is considered (e.g. Zhao et al. 2001), and in the ideal case the full travel-time noise covariance matrix is considered (e.g. Couvidat et al. 2005; Jackiewicz et al. 2008; Švanda et al. 2011).

The precise knowledge of the noise covariance matrix is essential, as this error propagates through the inversion procedure and translates into the realisation of the random error in the inverted maps. The determination of the full noise-covariance matrix is not straightforward. There are two approaches being used. Either the covariance matrix is read out directly from the data by measuring a large set of travel-time maps, or it is derived from the model using a Monte Carlo-like approach. Both approaches have their own advantages and disadvantages. The data-driven approach would probably be ideal; however, there might be time-varying systematic errors in the travel-time measurements, which will also affect the estimated error level of the inverted quantity. The model-driven methods rely on how precisely the used model agrees with the observed power spectrum of the waves. It would be ideal if one could minimise the effect of imprecise knowledge of travel-time noise on the knowledge of solar flows. Moreover, several studies (Woodard 2007; Jackiewicz et al. 2008) showed that the inversions for snapshot of the solar flows averaged over a few hours starts to be dominated by the random-noise component already at very shallow depths.

In Švanda et al. (2011) we suggested the use of the ensemble averaging approach. The ensemble averaging uses a strong assumption that the random-noise realisations in the individual representatives of the ensemble are independent, and hence in the average the random-noise level scales as , where N is the number of representatives. When N is large enough, it allows the noise term in the inversion cost function to be relaxed and thus decreases the effect of the possible inaccurate knowledge of the noise covariance matrix. The ensemble averaging approach is certainly not useful for studying the snapshots of the flows; however, when investigating the set of the representatives of the same phenomenon (supergranules, sunspots, etc.), it seems extremely powerful. In recent years it seems to have been a standard method of helioseismic research (e.g. Duvall et al. 2006; Duvall & Birch 2010; Švanda 2012; Birch et al. 2013; Švanda et al. 2014).

2. Inversion

The results presented in this note were obtained using the real data coming from the Helioseismic and Magnetic Imager (HMI) archive. Using the standard tracking and mapping pipeline (created and maintained within the German Science Center for SDO by H. Schunker & R. Burston), 24-h Dopplergram datacubes were tracked on a daily basis, from 8 May 2010 to 12 July 2010. These datacubes covered the central part of solar disc roughly having 60 degrees on a side in the Postel projection with a cadence of 45 s. Each datacube was processed in a standard way in a travel-time measurement pipeline. From each frame of the datacube, the mean image that captures most of the pattern of supergranulation was removed. The datacubes then underwent the spatio-temporal filtering for a set of filters (f to p4 ridge filters and eleven standard phase-speed filters) to retain only the waves with desired properties. The travel times were then measured from the filtered datacubes using the linearised Gizon & Birch (2004) approach for a set of distances for each of the spatio-temporal filters using the centre-to-annulus and centre-to-quadrant averaging schemes. Travel-time maps of the waves sensitive towards the surface (the f mode and first four phase-speed filters) in a centre-to-annulus geometry show clearly the pattern of divergence centres corresponding to supergranulation.

These travel-time maps were inverted in an inversion pipeline. The inversion code is written in MATLAB language and described in details in Švanda et al. (2011). It utilises the Born-approximation sensitivity kernels (Birch & Gizon 2007)1 and full travel-time covariance matrix measured directly from a large set of travel-time maps. The inversion implements the multichannel subtractive-OLA approach (Jackiewicz et al. 2012) with additional terms of the cross-talk minimisation and inversion weights localisation. The cost function is regularised strongly about these terms as they are both possible sources of biases. Both the travel-time and inversion pipelines were validated against the surface measurements (Švanda et al. 2013).

Three different combinations of sensitivity kernels were used. The first set combined all ridge-filtered kernels (f to p4 with a set of annulus radii in the range of 5–20 pixels with a step of 1 pixel, hence 240 independent kernels), the second combined all phase-speed-filtered kernels (eleven standard phase-speed filters with five distances for each of them as presented in Table 1 of Couvidat et al. 2006, in total 165 kernels), and the third is the inversion that combined all these kernels at once (405 kernels in total). The inversion using the combined filtering scheme was used only recently (Švanda 2013; DeGrave et al. 2014b,a) and the tests by Švanda (2013) showed that for a testing depth they provided flow estimates that were highly comparable. These tests, however, were performed in a different regime, when the cost function was regularised strongly about the random-noise term.

A different regime was used in this study by running the inversion suitable for the ensemble averaging approach. In this case, the requirement on the random-noise term of the inversion cost function can be relaxed, which allows the inversion to find the solution with a smaller misfit. In other words, the inversion produces a better localised averaging kernel. This kernel has fewer sidelobes. The sidelobes make the interpretations difficult.

The testing sample consisted of 17 474 individual supergranular cells, identified by a watershed algorithm in the centre-to-annulus travel-time maps. The segmentation algorithm and its implementation is described elsewhere (Švanda et al. 2014). The resulting flow maps of all three flow components were aligned about the centres of all detected supergranules and averaged; however, for the following discussion only the vx component (the component in the direction of solar rotation) was used. The claims are exactly the same for the remaining horizontal component vy and the problems become even worse for a weak vertical vz component, for which the random-noise constraint was relaxed beyond the reasonable signal-to-noise ratio.

thumbnail Fig. 1

Inversions for the horizontal flows in an average supergranule at depths 1.9 and 2.9 Mm, together with the display of the corresponding averaging kernels. The left column represents the inversion utilising combined ridge+phase-speed filtering, the middle column is for the ridge filters only, and the right column for phase-speed filters only.

Open with DEXTER

thumbnail Fig. 2

Same as Fig. 1, but for the depths of 4.3 and 6.2 Mm.

Open with DEXTER

3. Interpretation issues

The aim of this study was to test the credibility of the time–distance inversions for convective flows. In the misfit-dominated regime the main burden of the inversion quality is carried by the sensitivity kernels that come from forward modelling. Only a slight regularisation is performed by the random-noise term. A strong regularisation about both the cross-talk term and the term ensuring the localisation of the weights should guarantee the minimisation of biases.

All three different travel-time filtering schemes in the inversion provide similar averaging kernels and similar random-noise levels. Equation (2) implies that the estimates of the inverted flow should also be similar. This claim was verified by varying the trade-off parameters for the inversions utilising the same set of measurements.

Table 1

Tests of the inversion using the synthetic data.

The inverted flow maps together with the corresponding averaging kernels and noise levels are presented in Figs. 1 and 2. It seems that the above-mentioned implication does not hold when different sets of measurements are used in inversions. For each investigated depth (1.9, 2.9, 4.3, and 6.2 Mm) one can obtain different answers (even the sign) just by selecting a different set of measurements, when the random-noise levels are comparable (within a factor of two, given the magnitudes of the inverted flows with very large signal-to-noise ratios), and when the averaging kernels are very similar. It has to be noted that because of the strong regularisation of the cost function about the cross-talk term, the components of the averaging kernels that are not in the direction of the inversion (Eq. (2), first term on the right-hand side for αβ) are negligible (not shown). The second apparent issue is that within the inversions based on the same set of measurement, the sign of the flow reverses repeatedly as one goes deeper. This does not seem to be supported from the side of the theory of solar convection.

The inversions for statistical ensembles are extreme. In this extreme regime, however, we see clearly that they are not robust. There is no reason not to believe that in the case when the noise-regularisation term is stronger and the resulting flow maps are seemingly realistic, the problems persist to an unknown (possibly smaller) extent. For instance, the number of reversals can be manipulated by the strength of the noise-regularisation term. Generally speaking, the stronger the regularisation, the fewer the reversals. It is interesting that in the past the depth of the supergranulation was estimated from the depth where the inverted flow reversed its sign. Duvall et al. (1997), Duvall (1998), Zhao & Kosovichev (2003), and Švanda et al. (2009) reported the reversals depth within a few Mm below the surface. Recently, DeGrave et al. (2014b) used a similar time–distance inversion to this study to validate it by using the state-of-the-art simulation of solar convection and saw the supergranular flow reversal, even when it was not present in the simulation. The authors claimed that it was due to the shortcomings of the helioseismic inversion methods.

3.1. Verification of the inversion using synthetic data

A natural explanation of the issues described above would be an error (a programmer’s bug) in the inversion code, or a mathematical problem in the inversion itself (e.g. a degeneracy of the matrix to be inverted). To eliminate such possibility, a test involving the synthetic data was applied, similarly to the tests performed by Švanda et al. (2011).

A snapshot from the simulation of the convective Sun-like flows (Ustyugov 2006) was convolved with the appropriate set of the sensitivity kernels, in accordance with Eq. (1), to create a set of synthetic travel-time maps. A random realisation of the travel-time noise having the covariance matrix of the inversion was added to these maps to mimic the random excitation of solar waves. In total, three different sets of travel-time maps were computed, one for the ridge filters, one for the phase-speed filters, and one for all filters combined. These travel-time maps were inverted by using exactly the same inversion weights as used above.

The resulting flow maps were compared on a pixel-to-pixel basis to the flow maps obtained by the smoothing of the simulation with the inversion averaging kernels. Essentially, the left-hand side of Eq. (2) was directly compared to the first term of the right-hand side of the same equation, henceforth denoted as . The root-mean-squared (rms) value of the random-noise term (the second term on the right-hand side of Eq. (2)) was then compared to the expected level of the random noise, another output of the inversion. The results are summarised in Table 1.

One can see that the inversion behaves as expected. The inverted flow maps vinv are highly correlated to the expected flow maps vakern. The determined level of the random noise is consistent with the predicted level from the inversion. The worsts case (i.e. depth 6.2 Mm for the phase-speed kernels, which has the lowest correlation coefficient and signal-to-noise ratio) can be seen in Fig. 3. Even this worst case scenario demonstrates that the inversion behaves well and no sign of an obvious bug is found.

3.2. Does the inverted flow fit the observed travel times?

In SOLA inversions it is not guaranteed that the inverted flow model provides a reasonable fit to the observed travel times. It is obvious already from Eq. (1): in order to obtain forward-modelled travel times, the knowledge of the continuous flow model is needed, which is not returned by SOLA. A proper reconstruction is probably possible (a work in progress); however, for a simple illustration a much rougher estimate can be made, (3)where is a horizontally averaged averaging kernel. Such an estimate approximately fulfils Eq. (2) with the random-noise term excluded. The components of three different models from inversions based on different filtering schemes are shown in Fig. 4MP + R for the combined ridge+phase-speed filtering scheme, MR for ridge and MP for phase-speed filters based inversions, respectively. The suspicious alterations of the sign of the horizontal flow with depth are prominently visible. Models MP + R and MP are similar at depths less than 5 Mm and differ by sign deeper down. Model MR is significantly different from the other two.

From these models, forward-modelled travel-time maps were computed following Eq. (1) with the noise term neglected. This was fully justified by the ensemble averaging technique, which increased the signal-to-noise of the measured travel-time maps by more than two orders. A misfit between the forward-modelled travel times and the observed ones τa was estimated from (4)where M was the total number of travel-time measurements indexed by a. When is larger the forward-modelled travel times do not fit the observed travel times as closely.

The fit was evaluated for each considered vector-flow model separately for all three sets of travel-time measurements. The results are summarised in Table 2. One would naively expect the main-diagonal terms to be much smaller than the off-diagonal terms. From a visual inspection of both the observed and forward-modelled travel-time maps it becomes clear that all the fits are far from being satisfactory. Although they are very different (even by a sign), all three models fit the observed travel times equally badly.

The simple reconstruction of the continuous flow field probably also affects the results, but given that the alteration of the sign is visible already in the tomographic maps and that the averaging kernels are very confined around the target depth and have negligible sidelobes, it can be expected that a more realistic estimate of the continuous flow field will not dramatically change the conclusions.

thumbnail Fig. 3

A direct comparison of the expected horizontal flow values to the inverted ones on a pixel-to-pixel basis (only every tenth pixel is plotted to make the plot simpler). The solid line represents the linear fit to the crosses, the thin dotted lines the intervals of the predicted random-noise (1σ interval), and the thick dashed line is the line with the unity slope.

Open with DEXTER

4. Possible causes

There are several suspicious contributors possibly responsible for the problems.

Incompatibility of the travel time maps with the inversion weights. As has been pointed out many times in the past by various authors, it is crucial that any step in the data processing (mapping, filtering, the way the travel times are measured) must be taken into account when computing the corresponding sensitivity kernels. In our case, this was enforced as much as possible. The pixel size was exactly the same in both the data processing and the computation sensitivity kernels, and the spatio-temporal filtering was done not only using the same code, but also using the exact same files, in which the filters were stored. The Born sensitivity kernels (Birch & Gizon 2007) are consistent with the Gizon & Birch (2004) definition of the travel times. The issue with unknown impact are the possible non-linearities in both the forward and inverse problems.

Table 2

Match of the forward-modelled travel times for three considered models to the observed ones.

Difference of the power spectrum of the real data and the model. The sensitivity kernels are computed using the power spectrum, which comes from the model and which is considered Sun-like (Christensen-Dalsgaard et al. 1996). However, the eigenfrequencies deviate from the data eigenfrequencies for higher order modes (already p6 mode frequencies are considerably off), the eigenfrequencies are not available for very small wave numbers (k< 0.1 Mm-1). To investigate the influence of this issue, I recomputed sensitivity kernels, travel times, and inversions with an additional spatio-temporal filter, which limited the power spectrum to the region where the model power spectrum closely matched the power spectrum of the used datacubes. Hence, low wave numbers were filtered out, as were all the signals of oscillations beyond the p6 ridge. The introduction of this additional filter had only a minor impact on the results.

thumbnail Fig. 4

Estimated continuous horizontal velocities obtained from four tomographic maps at depths of 1.9, 2.9, 4.3, and 6.4 Mm. The magnitudes are saturated at levels ± 0.5 km s-1.

Open with DEXTER

thumbnail Fig. 5

Contributions of various filtering schemes to the inversion averaging kernel in case of the combined filtering scheme.

Open with DEXTER

Inaccuracy of sensitivity kernels. Should the sensitivity kernels forward-modelled from a reference solar model be different from those of the real Sun, we must expect this deviation to affect the results. The difference is expected to be small if the reference model is Sun-like. Unfortunately, there does not seem to be a direct way to measure the sensitivity kernels from solar data, perhaps except for the iterative inversions (Hanasoge & Tromp 2014; Hanasoge 2014), introduced to helioseismology only recently. There might be a way to at least verify the total integral of the kernel, which is the work in progress and will be reported on in a later paper. This point embodies both the effects of the differences between the reference model and the real Sun and the possible issues with the kernel computation, both theoretical and numerical. This point is also supported by DeGrave et al. (2014b), where the authors (using a realistic convection simulation) noticed a significantly decreasing correlation between the measured and forward-modelled travel times for ridge filters with higher orders (p3 and beyond) and also for larger phase speeds.

Mathematically posed problem. This possibility is closely related to the previous point. The assumption is that the sensitivity kernels are accurate (again, by accurate I mean those of the real Sun) to within 1%. The combined inversion combines 405 such kernels. These 1% errors translate through the inversion and may be amplified to an unknown extent. There are some hints that this is likely the most probably cause. The contributions of the averaging kernel from two different filtering schemes are shown in Fig. 5. One can see that in all four cases discussed in this research note the two components from different filtering schemes largely subtract from each other. Should one filtering scheme have an unknown systematical bias, the resulting real averaging kernel (i.e. the averaging kernels obtained by convolving inversion weights with the sensitivity kernels of the real Sun, which are unknown) is different and that would easily explain even the change of sign of the flow. A curious reader may object that the way out is not to combine the two filtering schemes. It is shown in Fig. 1 that the inverted horizontal flow based solely on the phase-speed filtering scheme reverses its sign between the depths of 1.9 Mm and 2.9 Mm, which is not realistic. One has to bear in mind that the subtractions occur in every inversion, e.g. the contribution of point-to-annulus measurements is subtracted from the contribution of the point-to-quadrant measurements, etc. Another indication for this explanation is that the magnitude of the inverted flow that uses the phase-speed filters seems to be a bit unrealistic (around 2 km s-1 for the depth of 1.9 Mm).

Both last points imply that the “real” inversion averaging kernels may be very different from those predicted by the inversion, perhaps even having a negative total integral. The cross-talk contribution is also not constrained, which is estimated from the cross-talk components of the averaging kernel.

5. Lessons learned

This work is based on many tens of thousands of CPU-hours of trial-and-error runs. In the case of the inversions suitable for the tomography of the flow snapshot (see e.g. Švanda 2013) the issues are not clearly visible; however, the question is whether the stronger regularisation of the solution (about the random-noise term in this case) removes the issues or hides them instead. By studying the literature on flow inversions in supergranules one has to conclude that all the inversions except for the inversion of Woodard (2007) indicated the reversal of the horizontal flow at various depths. It may be seen as suspicious, because the state-of-the-art simulations (Stein et al. 2009; Rempel et al. 2009) do not indicate such reversals. It may also easily be that the large-amplitude flows in supergranulation recently reported by time–distance inversions (Švanda 2012) are artefacts of similar issues, as in that case f mode and p1 and p2 ridges were used, necessarily leading to subtractions in the inversion.

From comparisons of the inversion results with the direct surface measurements (e.g. Ambrož 2005; Georgobiani et al. 2007; Švanda et al. 2007, 2013) it seems that the very near-surface inversions (hence involving f modes or acoustic waves with small phase speeds) do not suffer from the discussed problems. As we showed recently (Švanda et al. 2013), the flow inversion using the f mode ridge is not only highly correlated with the inferences from the surface granule tracking, but it also provides the properly scaled magnitude of the flow. Such validations against the independently obtained measurements fully justify the scientific results obtained on surface flow fields. It has to be pointed out that no larger subtractions occurred in those inversions.


1

All point-to-point sensitivity kernels were computed by a code kc3 written by Aaron Birch, which uses the normal mode summation approach with considered modes up to a radial order of 8, quadrupole source at the depth of 100 km, and the observation height of 300 km. The correlation time was chosen to be 48 s. A very fine grid in both the wave number (6 × 10-5 Mm-1) and frequency (43 μHz) space was chosen for a initial computation; however, the results do not seem to be extremely sensitive to the grid selection. The point-to-annulus and point-to-quadrant kernels were computed subsequently utilising spatial averaging.

Acknowledgments

This work was supported by the Czech Science Foundation (grant P209/12/P568). All computations were performed using the Sunquake compute cluster at Astronomical Institute of Academy of Sciences in Ondřejov, the tracked and mapped datacubes were obtained from data processing pipelines at Max-Planck-Institut für Sonnensystemforschung (MPS), Göttingen, Germany, which is funded by the German Aerospace Center (DLR). The solar measurements were kindly provided by the HMI consortium. The point-to-point travel-time sensitivity kernels were obtained using the code of Aaron Birch deployed at the MPS. This work was done within the institute research project RVO:67985815 to Astronomical Institute of Czech Academy of Sciences. This work benefits from couloir discussions following the talk Švanda: Current issues with time–distance inversions for flows at AsÚ at LWS Helioseismology Workshop #4: Solar Subsurface Flows from Helioseismology: Problems and Prospects which was held at Stanford University, July 21–23, 2014.

References

All Tables

Table 1

Tests of the inversion using the synthetic data.

Table 2

Match of the forward-modelled travel times for three considered models to the observed ones.

All Figures

thumbnail Fig. 1

Inversions for the horizontal flows in an average supergranule at depths 1.9 and 2.9 Mm, together with the display of the corresponding averaging kernels. The left column represents the inversion utilising combined ridge+phase-speed filtering, the middle column is for the ridge filters only, and the right column for phase-speed filters only.

Open with DEXTER
In the text
thumbnail Fig. 2

Same as Fig. 1, but for the depths of 4.3 and 6.2 Mm.

Open with DEXTER
In the text
thumbnail Fig. 3

A direct comparison of the expected horizontal flow values to the inverted ones on a pixel-to-pixel basis (only every tenth pixel is plotted to make the plot simpler). The solid line represents the linear fit to the crosses, the thin dotted lines the intervals of the predicted random-noise (1σ interval), and the thick dashed line is the line with the unity slope.

Open with DEXTER
In the text
thumbnail Fig. 4

Estimated continuous horizontal velocities obtained from four tomographic maps at depths of 1.9, 2.9, 4.3, and 6.4 Mm. The magnitudes are saturated at levels ± 0.5 km s-1.

Open with DEXTER
In the text
thumbnail Fig. 5

Contributions of various filtering schemes to the inversion averaging kernel in case of the combined filtering scheme.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.