Autoset Mw inversion limits

[claudiodsf]

Hi @krisvanneste,

I just pushed a commit (a15d171) with a new approach for autosetting Mw limits, between 90% of the minimum
of the spectral plateau and 110% of its maximum.

In my tests, this works very well and removes the need for the Mw_0_variability parameter.

However, I just tested it with noise weighting. I would be very glad if you could test it in your case, with frequency weighting.

Thanks!

[krisvanneste]

Claudio,
I will try to test it before the end of the week.
Is that OK with you?
Kris

[claudiodsf]

Sure, no rush!

[claudiodsf]

I also changed the clipping detection algorithm, see d1dbb0f.

[krisvanneste]

I compared the new approach for autosetting the Mw limits with the previous one using the Mw_0_variability parameter on my example with short/missing noise windows, requiring frequency weighting. In this particular case, I obtain better results with the new approach:

• Average Mw (4.1 instead of 3.9) closer to expected value (Mw=4.3 according to Bensberg network)
• Smaller misfits

However, I'm not sure this example is representative. I should do more testing, but unfortunately I'm still struggling to check the instrument responses for our older data.

The change you made looks simple, but there are 2 fundamental differences with the previous approach:

• the bounds are based on the min/max values of the part of the spectrum that is used to estimate the initial inversion parameters (and that part depends on the weighting used) instead of the initial Mw_0 value.
• the additional variability (10%) is defined in linear space, whereas Mw_0_variability corresponds to log space.
  However, in both cases you would expect a reasonable range, so I'm a bit surprised by the relatively large differences in my case.

Below I include the plots and logs obtained with the 2 approaches.

Fitted spectra (old approach):

Fitted spectra (new approach):

Log (old approach):
1306.ssp.log

Log (new approach):
1306.ssp.log

[krisvanneste]

Concerning the frequency weighting, in another program I use inverse frequency weighting. This is something I read in a paper by Kaneko & Shearer (2014) (https://academic.oup.com/gji/article/197/2/1002/617325):

We weight the fit inversely with frequency so that all parts of the spectrum seen in a log–log plot contribute equally (this generally improves the fit to the low-frequency part of the spectrum, which is defined by relatively few points).

If you don't have this paper, I can send it to you.

I was thinking of implementing this in sourcespec as well, but this will thus also require a strategy to determine the initial part of the spectrum (_freq_ranges_for_Mw0_and_tstar0 function)...

But first I will try to test the new clipping detection!

[krisvanneste]

I tested the new clipping detection with some clipped data, but it doesn't seem to be working.
Four stations (BEBN, LCH, MEMB and RCHB) have clipped S phases, but no warnings are given. Here is the trace plot:

I edited the code to turn on the debug plot for these particular stations:

RCHB may be more difficult to detect, but the other stations are quite obviously clipped.

What I do like about the new implementation is that the detection is limited to the wave type, which makes it possible to determine Mw from the P-phase when clipping is limited to the S-phase. However, for this to work properly, I think line 87 should be removed, as it redefines the end of the checked interval:

t2 = (trace.stats.arrivals['S'][1] + config.win_length)

I also implemented clip detection in my code, but it tends to overdetect. Robust clip detection is not so easy as it would seem...

[claudiodsf]

Hi Kris, thanks for your testing.

Mw_0_variablity was also in linear space (magnitude space). The old code is this one:

bounds.Mw_min = Mw_0 - config.Mw_0_variability
bounds.Mw_max = Mw_0 + config.Mw_0_variability

While the new one is:

bounds.Mw_min = np.nanmin(ydata[idx0: idx1]) * 0.9
bounds.Mw_max = np.nanmax(ydata[idx0: idx1]) * 1.1

This new code should provide larger bounds, which make it possible to find a better solution. It is particularly evident when using grid search and looking at misfit plots, as in the attached plots.

Concerning the inversion, it is made in logarithmic frequency spacing (see here) and this should mitigate the problem raised by Kaneko & Shearer.

But, of course, any improvement to the weighting strategy is welcome!

[claudiodsf]

Concerning the clipping detection: I guess there is a problem when the histogram peaks are at the edges of the histogram.

I had to smooth the histogram (orange curve) to avoid secondary peaks, but it obviously doesn't work in your case.

Your third case is hard because the clipping is just on one oscillation, but the other two should work.
Any idea on how to properly smooth the histogram?
I can also try figure out myself, if you send me the traces ;-)

[krisvanneste]

Sorry about the confusion regarding linear/log space. I forgot that you convert the spectra to seismic moment, whereas I was reasoning in amplitude...

[krisvanneste]

I will export the clipped records to miniseed and send them to you.

[krisvanneste]

Here are the miniseed files:
clipped_seismograms.zip
Do you need other data (stationxml, ...)?

[claudiodsf]

Thanks, I'm fine with miniSEED, since clipping is evaluated on uncorrected data.

[claudiodsf]

I moved the discussion on clipping detection to #23.

Let me know if you want to iterate more on Mw limits, or if you think that it is fine.

[krisvanneste]

I think it is fine as the new code indeed provides larger bounds.
The misfit plots you show are interesting. Is there an option in sourcespec to generate them?
I suppose the upper plot corresponds to the old code, and the lower one to the new code?

[claudiodsf]

The misfit plots are generated when using "grid search" (GS) for inv_algorithm.

I'm using more and more grid search (even if slower) because misfit plots give a clear view of parameter interdependence.

And, yes, old code on top plot.

[krisvanneste]

Great, I didn't know that. I will try it too!

[claudiodsf]

I think we can considered this as tested. Closing