Radiated energy and max_freq_Er

[krisvanneste]

The computation of radiated energy appears to be very sensitive to the choice of maximum frequency (configuration parameter max_freq_Er). With the default value of None, which implies using the full spectrum, computation of radiated energy fails in some cases, as shown in the plot below for stations GR.GRA1 and GR.GRC1. The log shows that this is because noise energy is larger than signal energy: skipping spectrum. However, the noise spectrum is everywhere below the signal spectrum, except at the highest frequencies (> 9 Hz).

I did some manual computations of the spectral integral of these spectra:

• signal spectrum (fmax = None): spectral_integral = 5725.75
• signal spectrum (fmax = 9): spectral_integral = 5567.65
• noise spectrum (fmax = None): spectral_integral = 16886.83
• noise spectrum (fmax = 9): spectral_integral = 0.75

I was able to avoid this problem by setting fmax to spectral_snratio_fmax from the header of the signal spectrum (8.77 Hz in the case of GR.GRA1) in the radiated_energy_and_apparent_stress function before computing the spectral integral.

[claudiodsf]

Hi Kris, I'm a bit surprised by the value of the noise spectrum integral when fmax=None.

I cannot see how the integral can be that big, given that the noise is consistently two orders of magnitude below the signal for most of the spectrum.

Any thoughts on that? Would it be possible for you to share a working example?

[claudiodsf]

By the way, these spectra are very strange: looks like you're using the response of a downsampled channel (nyquist at ~6 Hz) and what we see above 6 Hz is the antialiasing filter

[krisvanneste]

Claudio,

I couldn't believe it either at first. Looking at the _spectral_integral function though, it is clear that higher frequencies have a higher contribution to the result. The (displacement) spectral amplitudes are multiplied by 2 * np.pi * freq * np.exp(np.pi * t_star * freq)). With t_star=0.25, this gives the following results:

• freq=1: 13.78
• freq=5: 1594.48
• freq=10: 161853.00

Especially when using linear-spaced frequencies for the integration, it is understandable that if the noise spectrum is above the signal spectrum at the highest frequencies, this may result in a larger spectral integral, regardless of the situation at lower frequencies.

You are right about the odd shape of the spectra. These are data for an earthquake in 1983 from the Graeffenberg array in Germany. It's the only digital data we have. I had another look at the instrument response, but there is no anti-alias stage, so if there was indeed one, it is not corrected.

I will see if I can make a standalone example, but as I run everything on the fly without writing waveform data and event data, I need some time to figure out what is the best/easiest way to do this. I could also write the output spectra to miniseed files, but this will not preserve the header information contained inside.

[claudiodsf]

It's definitively the exponential term of the t_star correction that plays a role here!

I understand that it's complicated for you to create a standalone example, don't worry.
I'm just curious if you can produce a spectral plot with the non-attenuation-corrected synthetic spectra (config option plot_spectra_no_attenuation=True), to visually evaluate the effect of attenuation (t_star=0.25 is quite large!)

What's your suggestion to address this problem?

• Is the config option max_freq_Er enough?
• Should we use spectral_snratio_fmax (remember that there is a hard-coded value of 3 here)? Add an option to use it instead of max_freq_Er?
• Should we output in the log messages for which frequencies noise velocity spectrum is larger than the signal velocity spectrum?

Please let me know your point of view 😉

[krisvanneste]

Claudio,

It would still be useful to create a standalone example, but in the meantime here is the plot with the Brune fit without attenuation:

I will think about your suggestions and answer you tomorrow. Note that I am not so familiar with the theory behind radiated energy. I'm wondering if the high susceptibility to high-frequency content has been mentioned in the papers?

[claudiodsf]

A possible solution would be to add an option to estimate the radiated energy in a station-dependent spectral window (fmin, fmax) where the noise is below the signal.

Then, there is already a finite-bandwidth correction term in the code, that can be applied. Currently it is applied only for f>fmax, but it should be possible to correct also for f<fmin, as discussed in Di Bona and Rovelli (1988)

[claudiodsf]

On a side note, your t_star value seems really high to me, leading to quite a low corner frequency for a M4.6 event (stress drop around 0.1 MPa).
You might want to check the t_star - fc trade-off plots that you can obtain using the grid search algorithm.

[claudiodsf]

What I said above is not correct: fc and t_star are positively correlated, so reducing t_star would also reduce the corner frequency.

I think the problem with your fit is the "hole" in the spectrum above 6 Hz: the spectral fit wants to go very low at high frequency to fit this hole and this leads to high t_star values.

I think you should window your spectra to an fmax of 6 Hz

[krisvanneste]

A possible solution would be to add an option to estimate the radiated energy in a station-dependent spectral window (fmin, fmax) where the noise is below the signal.

Then, there is already a finite-bandwidth correction term in the code, that can be applied. Currently it is applied only for f>fmax, but it should be possible to correct also for f<fmin, as discussed in Di Bona and Rovelli (1988)

Claudio,
I think it's evident that the spectral integral should be limited to the range where the signal spectrum is above the noise spectrum. So it should not go beyond spectral_snratio_fmin and spectral_snratio_fmax, which are known when it is calculated. There is already a configuration parameter max_freq_Er. Maybe we should add min_freq_Er as well, and set them to spectral_snratio_fmax and spectral_snratio_fmin when they are not specified (which is what most users will probably do). Optionally, they could also be constrained to these values if they are beyond the frequency range where SNR > 3. Or do you mean to define these separately for each station in the configuration file?

[krisvanneste]

What I said above is not correct: fc and t_star are positively correlated, so reducing t_star would also reduce the corner frequency.

I think the problem with your fit is the "hole" in the spectrum above 6 Hz: the spectral fit wants to go very low at high frequency to fit this hole and this leads to high t_star values.

I think you should window your spectra to an fmax of 6 Hz

Thanks, I will try imposing an fmax of 6 Hz.

[claudiodsf]

I think it's evident that the spectral integral should be limited to the range where the signal spectrum is above the noise spectrum. So it should not go beyond spectral_snratio_fmin and spectral_snratio_fmax, which are known when it is calculated. There is already a configuration parameter max_freq_Er. Maybe we should add min_freq_Er as well, and set them to spectral_snratio_fmax and spectral_snratio_fmin when they are not specified (which is what most users will probably do). Optionally, they could also be constrained to these values if they are beyond the frequency range where SNR > 3. Or do you mean to define these separately for each station in the configuration file?

Ok, so I propose to rename the config parameter max_freq_Er to Er_freq_min_max and use it as follows:

• Er_freq_min_max = noise, noise (default): set min and max frequency to the "noise limits" (i.e. the frequency range where SNR>3)
• Er_freq_min_max = None or Er_freq_min_max = None, None: use the whole spectrum
• Er_freq_min_max = None, noise: use the lowest possible frequency, and set the max frequency to the "noise limit"
• Er_freq_min_max = 1, 10: use frequencies between 1 and 10 Hz
• Er_freq_min_max = 1, noise: use frequencies between 1 and the "noise limit"

What do you think?

[claudiodsf]

Thanks, I will try imposing an fmax of 6 Hz.

Curious to see the results 😉

[krisvanneste]

Ok, so I propose to rename the config parameter max_freq_Er to Er_freq_min_max and use it as follows:

* `Er_freq_min_max = noise, noise` (default): set min and max frequency to the "noise limits" (i.e. the frequency range where SNR>3)

* `Er_freq_min_max = None` or `Er_freq_min_max = None, None`: use the whole spectrum

* `Er_freq_min_max = None, noise`: use the lowest possible frequency, and set the max frequency to the "noise limit"

* `Er_freq_min_max = 1, 10`: use frequencies between 1 and 10 Hz

* `Er_freq_min_max = 1, noise`: use frequencies between 1 and the "noise limit"

What do you think?

Claudio,

This is quite flexible and should address most cases. Maybe use Er_freq_range similar to spectral_sn_freq_range ? On the other hand, there are also a number of _min_max configuration parameters, but these are mostly for constraining computation results...

[krisvanneste]

Curious to see the results 😉

Here is the result with freq2_broadb=6:

t_star is now 0.045, but the corner frequency is now very poorly defined.

When I lower the minimum frequency (bp_freqmin_broadb=0.1, freq1_broadb=0.1), I obtain this:

This looks much better, but the corner frequency remains low (0.7 - 1 Hz). The t_star values remain at 0.045.

I also managed to generate a standalone example with the same function I use to run sourcespec. If you unzip the attached ZIP file in a folder, then you should be able to run it with the following command:
python <path_to_sourcespec>/bin/source_spec.py --configfile <ssp_input_folder>/651.ssp.conf --qmlfile <ssp_input_folder>/event.qml --evname "LIEGE (BE)" --evid 651 --station_metadata <ssp_input_folder>/response.xml --run_id test_variable_winlen_SH --run_id_subdir --outdir <ssp_input_folder>/ssp_output --trace_path <ssp_input_folder>
651_ssp.zip

[claudiodsf]

The last example looks much better, and the t_star value of 0.14-0.17 s is more reasonable.

One can see a bit of the effect of the high pass filter on the left side of the spectrum: maybe you can filter at lower frequencies, e.g. bp_freqmin_broadb = 0.01 while keeping the windowing at freq1_broadb = 0.1

[claudiodsf]

Maybe use Er_freq_range similar to spectral_sn_freq_range ?

ok !

[claudiodsf]

Here is a possible solution: bae2ddf

There is however a problem when using spectral_snratio_fmin and/or spectral_snratio_fmax (noise options in Er_freq_range): those values are computed over a "valid" snratio spectrum (this line), which is defined from the config parameter spectral_sn_frequency_range (see here).

So, if, for example, spectral_sn_frequency_range = 0.1, 2.0, then Er_freq_range will be forced to be between these two values.

Do you remember why we made this choice? Shouldn't we compute spectral_snratio_fmin and spectral_snratio_fmax using the whole spectrum?

[claudiodsf]

P.S. I'm going to add new commits to introduce another feature in ssp_radiated_energy requested by another user: taking into account for the radiation pattern.

This should not affect the current issue.

[krisvanneste]

Here is a possible solution: bae2ddf

There is however a problem when using spectral_snratio_fmin and/or spectral_snratio_fmax (noise options in Er_freq_range): those values are computed over a "valid" snratio spectrum (this line), which is defined from the config parameter spectral_sn_frequency_range (see here).

So, if, for example, spectral_sn_frequency_range = 0.1, 2.0, then Er_freq_range will be forced to be between these two values.

Do you remember why we made this choice? Shouldn't we compute spectral_snratio_fmin and spectral_snratio_fmax using the whole spectrum?

I will look into this tomorrow.

[krisvanneste]

Configuration parameter spectral_sn_freq_range was introduced here: b2b5c77
This was before "my time". I usually set it to None.

[claudiodsf]

I usually set it to None.

I will then change:

snr_valid_freqs = valid_freqs[valid_snratio >= 3]

to

snr_valid_freqs = freqs[weight.snratio >= 3]

[krisvanneste]

This could be more robust, but then option spectral_sn_freq_range will only be used to compute the spectral SNR (spectral_sn_ratio), which is saved in the header. This will then not be in agreement with spectral_snratio_fmin/fmax, which may lead to confusion.

[krisvanneste]

I think the option spectral_sn_freq_range may still be useful for cases where the SNR may have different portions where it is above the threshold. You could have some portions of the spectrum at low or high frequency that are separated from the main portion where SNR is also >= 3, but which you do not want to include. So I guess someone who is using this option knows what [s]he is doing.

[claudiodsf]

I think the option spectral_sn_freq_range may still be useful for cases where the SNR may have different portions where it is above the threshold. You could have some portions of the spectrum at low or high frequency that are separated from the main portion where SNR is also >= 3, but which you do not want to include. So I guess someone who is using this option knows what [s]he is doing.

I think that in this case people should just use (station-specific) spectral windowing, i.e. the options freq1_* and freq2_*.

I would keep my proposal above and maybe rename the following two config options for more clarity:

spectral_sn_min --> average_spectral_sn_min
spectral_sn_freq_range --> average_spectral_sn_freq_range

...or just keep the current names and explain them better in the config file comments.

What do you think?

[krisvanneste]

OK, that's fine with me. Maybe a better explanation is sufficient.

[claudiodsf]

Ok, just pushed the latest commits. I guess we can close this 😉