https://orcid.org/0000-0003-3309-6785
1 INTRODUCTION
We would like to thank Prof. Helffrich for this opportunity to further elaborate on our work on the determination of crustal thickness from ambient noise autocorrelations (Becker & Knapmeyer-Endrun 2018, BK18 hereafter). In his comment, Prof. Helffrich states that insufficient methodological details are given in BK18 as his attempt to reproduce fig. 7(b) failed. The increase in reflectivity seen at later lag times, between 16 and19 s, was absent in Prof. Helffrich’s result. He suggests that this difference is caused by differences in the filter design of the bandpass filter applied to the correlations caused by an insufficient description of the filter parameters in BK18 and that the signal between 16 and 19 s is due to beating caused by the narrow passband of the filter. Furthermore, Prof. Helffrich states that the importance of the length of the weighting window when identifying peaks in the second derivative of the envelope of the correlation is not discussed sufficiently in BK18.
In the following, we comment on both criticisms and provide further details. In particular, we show that the applied filtering does not lead to beating and is not the cause of the discrepancy in the results and that the influence of the weighting window on peak identification is minor.
2 FILTERING AND PROCESSING OF AUTOCORRELATIONS
The overall processing procedure summarized in the comment is correct, although some pre-processing steps, which are also mentioned in BK18, such as removing the mean and linear trend and eliminating time windows containing data gaps or spikes, are omitted. As described in BK18, we applied a zero-phase fourth-order digital Butterworth highpass with a cut-off frequency of 0.5 Hz to the raw data, and a zero-phase eighth-order Butterworth bandpass filter with a passband of 1–2 Hz to filter the correlations. In BK18, we did not consider that the filter implementation for bandpass filters applies two passes, so that the selected filter order is multiplied by two; this means we incorrectly stated a filter order of four for the Butterworth bandpass filter. We cannot give any further details on these filters, as asked for in the comment, as the passband and order completely describe a Butterworth filter, which by definition has a roll-off of 20 dB per decade for every pole (e.g. Rorabaugh 1993). We did not perform the data processing in Python (i.e. SciPy and NumPy; Millman & Aivazis 2011), as is assumed in the comment, but in Matlab®, where no further parameters can be provided in the filter design. As the Matlab® documentation points out that transfer function design for the Butterworth filters can get numerically unstable even at low orders, we also checked the filter responses against the zero-pole-gain design type to ensure filter stability.
Prof. Helffrich hypothesizes that the details of the applied filter have an important influence on the result and that the signal between 16 and 19 s in the autocorrelations for station PBUR is probably due to beating as the passband of the filter is quite narrow. However, this is not the case as even the unfiltered data show the discussed features between 16 and 19 s (Fig. 1a), in contrast to fig. 1(b) of the comment. These features are also persistent for different filter orders, that is, appear for smaller and larger roll-offs than in the filter we used originally (Figs 1b and c). We also increased the width of the passband to 1–4 Hz without any striking differences regarding the existence of the signal between 16 and 19 s lag time (Fig. 1d), though amplitudes of the autocorrelation generally appear lower for lag times larger than 5 s. Results shown here are for 1 hr time windows used to calculate the autocorrelation to be directly comparable to the results shown in the comment. However, similar results are found for 3 hr time windows, as originally used in BK18. We also note that it is not uncommon to use narrow passbands when analysing auto- and cross-correlations in terms of Moho-reflected body waves, for example, 0.5–1 Hz (Poli et al.2012; Tibuleac & von Seggern 2012) or 2–4 Hz (Gorbatov et al.2013; Kennett et al.2015; Saygin et al.2017).
Investigation of the influence of filtering on the autocorrelations. Shown are daily linear stacks of phase autocorrelation calculated over 1 hr time windows for the year 2013 for station PBUR, similar to fig. 7(b) in BK18 and fig. 1 of the comment by Prof. Helffrich. a) Raw stacks without any filtering applied to the autocorrelations. b) Autocorrelations filtered with a zero-phase Butterworth bandpass filter with a passband of 1–2 Hz and an order of 4. c) Autocorrelations filtered with a zero-phase Butterworth bandpass filter with a passband of 1-2 Hz and an order of 16. d) Autocorrelations filtered with a zero-phase Butterworth bandpass filter with a passband of 1-4 Hz and an order of 8.
As a conclusion, the narrow passband and filter order are not responsible for the increased reflectivity at later lag times and the differences must originate from other processing steps. We note that the main effect of filtering the data, with either passband and independent of filter order, is an amplitude reduction around 15 s lag time. This amplitude reduction was interpreted as reflectivity change caused by the Moho in BK18, whereas the amplitude increase between 16 and 19 s was not interpreted in terms of crustal structure.
One possible origin could be the implementation order of the pre-processing steps of downsampling, removing the mean and linear trend and tapering, as some of them are omitted from the summary provided by Prof. Helffrich. Theoretical considerations indicate that the order in which these operations are applied should have no influence on the result as they are all linear. We confirmed this by comparing results for different processing orders. The downsampling algorithm as implemented in Matlab® by default uses an eighth-order lowpass Chebyshev Type 1 IIR filter. Using an Finite impulse respone (FIR) filter of order 30 instead leads to slight differences in the pre-processed seismograms (Fig. 2a), but does not have any clear effect on the autocorrelation result especially after stacking (Figs 2b and c). This is also in agreement with the observations by Prof. Helffrich that different downsampling algorithms cannot explain the difference in the autocorrelations. Regarding his further tests of steps 1 and 2, as stated in BK18, the removal of the instrument response did not affect our results, and neither did the use of 1 or 3 hr time windows for the calculation of the correlations.
Investigation of the influence of the filter used in decimating the seismograms. (a) One arbitrary hour of data from PBUR on 2013 January 1 after pre-processing. Black line corresponds to decimation with IIR Chebychev lowpass of order 8, grey line to decimation with FIR lowpass filter of order 30. (b) Phase autocorrelations for the same hour of data as shown in (a). (c) Linear stack of phase autocorrelations calculated over 1 hr windows (including the one shown in (b)) for one day of data.
Another possible source of the different results could be the spectral smoothing algorithm. In the following we will briefly summarize the method used in BK18. After cutting the data into 3 hr time windows, applying the above pre-processing steps and a zero-phase fourth-order digital Butterworth highpass filter with a cut-off frequency of 0.5 Hz, a spectral smoothing algorithm is applied to the amplitude spectrum. A fast Fourier transform is used to calculate the spectrum. For each frequency, the average amplitude within a long- and a short-frequency window is calculated. These windows are centred on the frequency under consideration. For the 3 hr time window of data, the short-frequency window consists of 10 samples and the long-frequency window consists of 10 000 samples, with each sample corresponding to 4.63 × 10−5 Hz. Next the ratio between the average amplitudes in the short window and the long window is determined for each frequency. In addition, the mean of all short-term long-term average ratios over the whole spectrum is calculated and compared to the individual ratios for each frequency. If an individual single ratio is larger than the mean ratio, indicating significant amplification, the spectral amplitude at the considered frequency is replaced by the long-term average. The description in BK18 is based on the inverse quantity, that is, the amplitude ratio between the long and the short window, with spectral amplitudes being replaced if an individual ratio was smaller than the mean. The description here reflects the actual implementation. Both descriptions are not mathematically equivalent, but we investigated the influence of the two different implementations on the resulting phase autocorrelations for station PBUR and do not see a difference in the result (Fig. 3).
Investigation of the influence of the different implementations of spectral smoothing on the autocorrelation results. Shown are daily linear stacks of phase autocorrelation calculated over 1 hr time windows for the year 2013 for station PBUR, similar to fig. 1 of the comment by Prof. Helffrich. (a) Stacks of phase autocorrelations calculated with the spectral smoothing algorithm as described in the text. (b) Stacks of phase autocorrelations calculated with the spectral smoothing algorithm as originally described in BK18.
Additionally, in the case of PBUR, the application of spectral smoothing does not have a strong influence on the resulting stacked autocorrelations and the increase in reflectivity at later lag times is also apparent when omitting the spectral smoothing from the processing flow, again also before applying any filtering to the autocorrelations (Fig. 4). This means that the spectral smoothing does not cause the signals under discussion. The minimum in amplitude at around 15 s lag time can also already be detected in the phase-weighted stacks without filtering, and with or without spectral smoothing (Figs 4c and d).
Investigation of the influence of spectral smoothing on the autocorrelations. Shown are stacked phase autocorrelation using 3 hr time windows for the year 2013 for station PBUR. The left-hand column shows stacks of autocorrelations calculated without spectral smoothing, the right-hand column shows stacks of autocorrelations calculated with spectral smoothing. (a,b) Linear stack without filtering: (c,d) phase-weighted stack (Schimmel & Paulssen 1997) without filtering; (e,f) linear stack after application of zero-phase digital Butterworth bandpass with corner frequencies 1 and 4 Hz and order 8; (g,h) phase-weighted stack after application of zero-phase digital Butterworth bandpass with corner frequencies 1 and 4 Hz and order 8.
3 WEIGHTING OF SECOND DERIVATIVE
In order to determine the reflectivity change associated with the Moho, we propose in BK18 to use a weighted second derivative of the envelope of the autocorrelation and to identify local maxima of this function within a lag time window determined from prior information. The weight is defined as the average value of the envelope over a 5 s time window centred on the considered time sample. Prof. Helffrich argues that the length of the weighting time window plays a significant role in the peak selection process. It is however equally possible to determine local maxima without weighting of the second derivative. An example for this is shown in Figs 5(a)–(c).
Investigation of the effect of weighting of the second derivative of the envelope on peak detection for station PBUR. (a) Time–frequency domain phase weighted stack of autocorrelations calculated over 3 yr of 3 hr time windows and filtered between 1 and 2 Hz with a zero-phase digital Butterworth bandpass of order 8. The dashed line shows the envelope, the grey region shows the target window and the red line shows the picked lag time. (b) Normalized second derivative of the envelope in (a) without weighting. (c) Weighted second derivative of the envelope in (a), using the average amplitude of the envelope within a 5 s time window as weight. Red triangle marks maximum of the weighted derivative within the a priori time window. (d) Time-frequency domain phase weighted stack of autocorrelations calculated over 3 yr of 3 hr time windows and filtered between 2 and 4 Hz with a zero-phase digital Butterworth bandpass of order 8, zoomed on the first 10 s of lag time related to shallow structure. (e) Normalized second derivative of the envelope in (d) without weighting. (f) Weighted second derivative of the envelope in (d), using the average amplitude of the envelope within a 5 s time window for weighting. Grey dashed lines mark the three largest local maxima identified in the weighted derivative that correspond to amplitude minima of the envelope.
It can be seen that the presence of the local maximum at around 15.6 s in the weighted second derivative does not arise from the weighting; rather, the maximum is already present in the second derivative without weighting. The weighting simply enhances its amplitude in comparison to later peaks starting at 20 s. Without a priori information, for example, to define a target window, this could help in the identification of the relevant peak related to a reflectivity change at the Moho. However, as relative amplitudes before 15 s are also amplified by the weighting, its effect is very minor in the case of PBUR. Things are different for other stations, for example, for GEC2 (Fig. 6).
Same as Fig. 5 for station GEC2.
Here, later peaks in the second derivative with higher amplitudes are effectively suppressed by the weighting and earlier peaks not significantly enhanced, so that the maximum around 12 s becomes more ‘global’ by the weighting. Without the weighting, it is however in both examples still possible to identify the local maximum related to the reflectivity change at the Moho.
We agree that the availability of prior information is a necessary prerequisite for the automatic determination of the Moho depth, first of all because the measured lag times need to be converted to depth based on some velocity model, but also to define a target window. It is also true that when trying to identify shallow structural details, different filtering, weighting or modelling options might be necessary. For example, filtering between 1 and 2 Hz is not suitable to detect shallow crustal features. Indeed, Saygin et al. (2017) filtered at higher frequencies, between 2 and 4 Hz, to identify P-wave reflections within the Jakarta Basin. As the application of the automatic picking based on second derivatives of the envelopes to determine the Moho reflection proved difficult when considering this frequency range in BK18 (fig. A4), this might also be the case when this automated processing is applied to signals at shorter lag times, that is, from shallower depth above 10 km which were not considered in BK18. However, variations in autocorrelation amplitude within the first 5 s are clearly visible in fig. A4 of BK18, showing autocorrelations for 25 broad-band stations filtered between 2 and 4 Hz. So while the automatic picking algorithm might fail, information on shallow structure is still contained in the autocorrelations. Figs 5 and 6(d)–(f) show that changes in the autocorrelation function can still be properly identified based on the second derivatives of the envelopes. Their relation to shallow structure needs further investigation, though.
4 CONCLUSIONS
Clearly, we would expect the method and results presented in BK18 to be reproducible and we are concerned that Prof. Helffrich failed in his attempt to reproduce fig. 7(b) of BK18, showing the temporal consistency of autocorrelation results for one out of 57 stations analysed in the above paper. Bandpass filtering of the autocorrelations with insufficient information on the used filters, which was proposed as a source of the discrepancy by Prof. Helffrich, can be excluded as a cause, though, as the high amplitudes between 16 and 19 s lag time absent in his results are already present in our data set before any filtering is applied. Besides, the data analysis software we are using contains the standard roll-off for Butterworth filters, which are solely defined by order and frequency band, as given in BK18, so there is no additional information that can be provided regarding the filters. We discuss some other processing steps that might be responsible for the discrepancy in results, that is, pre-processing of the seismograms, the spectral smoothing algorithm or the phase autocorrelation itself. However, we find no evidence of any influence of the order of the pre-processing steps of the filter used in decimating the data on the results, and even without spectral smoothing, our stacked autocorrelations still show an amplitude increase between 16 and 19 s lag time. The weighting of the second derivative of the envelope helps in identifying local maxima related to reflectivity changes in the autocorrelation, but is not a necessary step, as these local maxima can also be identified in the second derivative without weighting. As such, the significance of the window length used for weighting is only minor, and it would still be possible to detect changes in the correlation function at lag times shorter than 5 s based on peaks in the second derivative of the envelope of the autocorrelation.
REFERENCES
Author notes
Now at: Ingenieurgesellschaft Nordwest mbH, Oldenburg, Germany
Now at: University of Cologne, Earthquake Observatory Bensberg, Bergisch Gladbach, Germany
