swiss-seismological-service / scdetect

A computationally efficient earthquake detection module for SeisComP
https://scdetect.readthedocs.io
GNU Affero General Public License v3.0
15 stars 6 forks source link

Pick handling for processed/projected streams #3

Closed damb closed 3 years ago

damb commented 3 years ago

Certain phases (e.g. Sg phases) are often picked by seismologists using a projected stream triplet (see also https://service.iris.edu/irisws/rotation/docs/1/help/). A projected stream triplet is identified by a dedicated stream identifier, the projected stream identifier (i.e. in case of a ZRT projection the channel part of the identifier is modified, accordingly). This projected stream identifier is used by SeisComP when storing DataModel::Pick objects at the database. However, the original stream identifier is not stored. Though, the original stream identifier is required in order to generate the template waveforms.

Questions: Is there a certain rule in order to obtain the original stream identifier from a projected stream identifier? If yes, how does it look like?

References:

damb commented 3 years ago

@kaestli

kaestli commented 3 years ago

Not sure whether I get it. You have a pick e.g. on HHT, and you would like to generate a template for that event and that station? And derive the channels to cross-correlate from the pick?. I guess it would be HH* (all three channels)

(but would you cross-correlate only on HHZ, if you start from a P pick labelled with HHZ? in the picker window of seiscomp e.g., all 3 streams are displayed, and you counter-check with all of them when picking. The selection of Z for the pick label is somewhat arbitrary.

damb commented 3 years ago

Not sure whether I get it. You have a pick e.g. on HHT, and you would like to generate a template for that event and that station? And derive the channels to cross-correlate from the pick?. I guess it would be HH* (all three channels)

As we discussed, already, the concept of a station, in fact, does not exist (at least for the time being). The module operates fully stream based. A stream is identified by a stream identifier i.e. NET.STA.LOC.CHA. Currently, wildcards are not supported, neither. Within a template configuration a corresponding stream is specified as follows. E.g.

{
    "initTime": 20,
    "templateWaveformStart": -2,
    "templateWaveformEnd": 10,
    "waveformId": "CH.GRIMS..HHE",
    "templateWaveformId": "CH.GRIMS..HHE",
    "filter": "",
    "templateFilter": "",
    "templatePhase": "Sg"
}

Now, assume that it turns out that there is no pick for phase Sg related with CH.GRIMS..HHE for an origin with originId (the origin identifier is not part of the template stream definition but of the corresponding template detector config), but with CH.GRIMS..HHT and CH.GRIMS..HHR. So, how to determine the original channel(s)? Is this even possible? I.e. which pick time refers to the HHN channel and which one to HHE?

for the rotation HHT, actually information of all 3 channels was used.

That's not correct. If I understand M_{2D} from https://service.iris.edu/irisws/rotation/docs/1/help/ correctly, only the 2 horizontal components are used.

damb commented 3 years ago

Subsource codes are arbitrary and thus need to be determined using a best effort approach. Which one is still not fully clear to me, though.

kaestli commented 3 years ago

When preparing a template event, what you need from a pick at the station is only the time (telling you which time interval of the traces you would have to extract in best hope that they contain waveform energy of the event. However, as all streams of a station (or: sensor location) are recorded at the same place, this energy is also expected to arrive on all channels at the same time. However, which channel was used for picking does not necessarily relate to the question which channel is good for cross correlation. This is basically an operator decision. The waveform arrical time at any stream 1Net.2Sta.3Loc. is good enough to define the relevant time window on any (including any other) stream 1Net.2Sta.3Loc. Ok, you may assume that if channel X was picked, this implies that channel X may contain relevant waveform energy, but not necessarily that it is also good for cross correlation (imagine, e.g. a clipped HHZ signal, which is completely fine for picking, however, it does not describe the event's specific trace shape, so for cross corrlation, one would rather use the colocated strongmotion channel HGZ). For selecting which channel of RT trace to compare to which channel of template trace, just look at the channel pairs defined in the configuration. At the pick time of the event at the station/sensor location you have to look only in order to select the relevant time window.

(side comment: it would actually be possible to define a cross correlation on a HHR (2-dimensional rotation), HHT 2- or (3-dimensional rotation), or similar auxiliary stream. That setram would then need to be continuously computed on the fly following the maths of the IRIS link, and the relative location of the station to the template event. Seismologically it even makes sense (as it describes not only a similar waveform shape, but also a similar origin locaton). However, this is an advanced feature which we should not consider at the time being.)

damb commented 3 years ago

Thanks @kaestli.

For selecting which channel of RT trace to compare to which channel of template trace, just look at the channel pairs defined in the configuration. At the pick time of the event at the station/sensor location you have to look only in order to select the relevant time window.

The only issue still not solved is the example when there are picks for a certain phase on both XXR and XXT channels. One possible solution of course could look like: Take the mean pick time of both picks and include both horizontal components when configuring the Detector e.g. XXN and XXE. Would this proposal be good enough for the time being?

damb commented 3 years ago

(side comment: it would actually be possible to define a cross correlation on a HHR (2-dimensional rotation), HHT 2- or (3-dimensional rotation), or similar auxiliary stream. That setram would then need to be continuously computed on the fly following the maths of the IRIS link, and the relative location of the station to the template event. Seismologically it even makes sense (as it describes not only a similar waveform shape, but also a similar origin locaton). However, this is an advanced feature which we should not consider at the time being.)

I thought about this already. Though, this make only sense if all horizontal components are part of a Detector configuration.

kaestli commented 3 years ago

You will not have the same phase picked on multiple channels of one sensorlocation (it would not make sense; it would mean that e.g. the Sg wave arrived at the coordinates of the sensor location either at time p1 or at time p2. A location algorithm would fail with this). However, you may have different phases picked at the same station, with different arrival times, and thus different time windows. if cross correlation time windows are short, they will contain only or primarily that arrival, and having two templates for intervals of picks corresponding to two phases actually adds useful independent information describing the template event.

(I guess these considerations are beyond what has been done before; @luca-s , how is it in seismatch? Probably they consider only a template snipplet around the first (P) arrival without any further considerations?)

damb commented 3 years ago

You will not have the same phase picked on multiple channels of one sensorlocation (it would not make sense; it would mean that e.g. the Sg wave arrived at the coordinates of the sensor location either at time p1 or at time p2. A location algorithm would fail with this).

This makes sense. (I messed something up with the catalog I've got available, here.)

luca-s commented 3 years ago

@luca-s , how is it in seismatch? Probably they consider only a template snipplet around the first (P) arrival without any further >considerations?

I think we have already covered this in one of the email exchanges. SeisMatch uses only P picks and the template length is long enough to include both P and S. I do not believe we will ever need to consider the S phase and neither the rotation on R T.

luca-s commented 3 years ago

@damb if we ever need to perform a rotation, we can still apply the same logic that SeisComP applies in scolv (I actually used that in scrtdd), but I doubt we will need it

EDIT: this probably doesn't solve your problem. I have understood your question only now. You are just interested in finding the component codes, without projection/rotation. Sorry for the misunderstanding. I believe my next comment can help you with that.

luca-s commented 3 years ago

@damb You might find this function useful.

//! Returns the best matching streams for the vertical and horizontal components
//! having the stream code set to streamCode (without component code, first two letters).
//! Returns true if the resulting 3 components are forming an orthogonal system,
//! false otherwise.
//! The method tries to find 3 orthogonal components and select the first with
//! the lowest dip (largest Z value) as vertical. The remaining two are returned
//! as 1st horizontal (Y axis or North) and 2nd horizontal (X axis or East)
//! respectively.
//! NOTE: Each of the comps entries in res can be nullptr.
SC_SYSTEM_CORE_API bool
getThreeComponents(ThreeComponents &res, const SensorLocation *loc, const char *streamCode, const Core::Time &time);
damb commented 3 years ago

@damb You might find this function useful.

@luca-s, thanks. However, IMO this will be of use, only, in case specifying subsource codes using wildcards are allowed.

damb commented 3 years ago

To conclude, for the time being, the following implementation will be provided:

luca-s commented 3 years ago

@damb so the configuration requires the user to provide both the waveformId and a phase type (Pg, Sg, etc)? If so than it makes sense to me and the issue is solved (for now? :) )

kaestli commented 3 years ago

For full definition, and full transparency on the behaviour, yes. (The [from my point of view :-)] probably better alternative to hardwiring the use of the first phase in the code and hiding the topic away for the time being)

damb commented 3 years ago

If so than it makes sense to me and the issue is solved (for now? :) )

I'll close the issue, when the implementation is done.