A possible issue with synthetic waveforms generated by MTUQ from CPS GF database

[Jinyin-Hu]

Hello,

I'm encountering a problem with the synthetic data generated by MTUQ using CPS Green's function (GF) database. The synthetic waveforms from MTUQ differ from those created by my code. It seems to me there is something wrong there with reading/processing the CPS database in MTUQ. Any suggestions will be helpful. Many thanks.

Here are the details. I used CPS to generate the GFs for six stations, a given 1D model and a given depth, then read them into MTUQ as
model = 'mdj3'
db = open_db('/mypath/grn_2017_2d/mdj3', format='CPS', model=model)
greens = db.get_greens_tensors(stations, origin)
greens_sw = greens.map(process_sw)
where:
process_sw = ProcessData(
filter_type='Bandpass',
freq_min=0.02,
freq_max=0.05,
pick_type='CPS_metadata',
CPS_database='/mypath/grn_2017_2d/',
CPS_model=model,
window_type='surface_wave',
window_length=350,
capuaf_file=path_weights,
apply_scaling = False )

To generate synthetics, I output the greens_sw as:
from mtuq.misfit.waveform import level2 as level2
stations = level2._get_stations(data_sw) #data_sw = data.map(process_sw) is the surface wave data
components = level2._get_components(data_sw)
# collapse structures into NumPy arrays
obs = level2._get_data(data_sw, stations, components) #n_stat x n_comp x npts
greens_mtuq = level2._get_greens(greens_sw, stations, components) #n_stat x n_comp x 6 x npts, and in up-south-east convention

So far, I obtained the observations, array 'obs' and the elementary seismograms corresponding to six MT parameters, greens_mtuq. For comparison, I also used my own code to prepare the elementary seismograms corresponding to six MT parameters, greens_my, as
from CPS import load_CPS_greens_tensor #my CPS.py is attached at the bottom.
greens_my = load_CPS_greens_tensor(origin, stations, path_to_data, dt=0.05, npts=11000,
filter_dict={'fmin':0.02, 'fmax':0.05, 'corners':4, 'zerophase':False})
#the shape is 6 x n_stat x n_comp x npts
Then, I manually cut the 'greens_my' using the windows defined similarly to MTUQ. And here greens_my uses the north-east-down convention. After converting it into the up-south-east conventions by following the same code in MTUQ, I plot the comparison between greens_mtuq and greens_my for one station, as in the figure below. They are quite different.

The CPS GFs here are actually for the velocity model from Dreger et al. (2021) to study the DPRK2017 event. To further confirm this difference, I generated the waveforms for the MT solution from Dreger et al. (2021) by a linear combination of MT and greens_mtuq or greens_my, as
pred_mtuq = numpy.einsum('scet,e->sct', greens_mtuq, mt)
pred_my = numpy.einsum('esct,e->sct', greens_my, mt)
Here are two figures for each predicted dataset:
1)This is for predicted data (red) using greens_mtuq, compared with the observations (black). The variance reduction is 71.1%.

2)This is for predicted data (red) using greens_my, compared with the observations (black). The variance reduction is 90.9%, which is expected by Dreger et al. (2021).

Reference: Dreger et al. (2021). Path Calibration of the Democratic People’s Republic of Korea 3 September 2017 Nuclear Test, SRL, 10.1785/0220210105.
CPS.py

[rmodrak]

Hi @Jinyin-Hu , Many thanks for the report. Stepping back, below are some very general comments that I hope might be relevant. Feedback Robert or others might also be helpful.

-For adding CPS support, our starting point was to compare FK and CPS Green's functions. We noted that they differed only by a sign convention. So, to obtain the CPS machinery, we just applied these sign changes to the already existing FK machinery.

-The above sounds straightforward, but in reality, doing the FK versus CPS comparison was something of an undertaking. Neither solver compiled using modern Fortran compilers. We used legacy gfortran options and made other modifications to get them to compile, but then encountered 'underflow' warnings at runtime. Still though, in the end we were getting CPS versus FK comparisons much closer than the comparisons posted above.

[rmodrak]

Drawing from private functions such as _get_greens may be error prone? Perhaps try the following to generate unprocessed MTUQ synthetics

db = open_db(path_cps, format='CPS')
greens = db.get_greens_tensors(stations, origin)
synthetics = greens.get_synthetics(mt)

then compare with your own unprocessed synthetics?

[Jinyin-Hu]

Hi @rmodrak , thank you for your reply and information. It's good to know how the CPS support is implemented in MTUQ. I just tried greens.get_synthetics(mt), but still got the same waveforms. Actully, I found _get_greens and _get_data functions are used in misfit class. So I assume MTUQ also use them to do the inversion. Aother clue is that I calculated the GFs in displacement and observed data is also displacement. When read the data in MTUQ, I already set [units:m', 'type:displacement'] for function 'read'. Does MTUQ has some default setting about that, e.g., if all waveforms and GFs will be converted to velocity?
Many thank.
Jinyin

[Jinyin-Hu]

Hi @rmodrak , I might found it. Same thought about the type of data and GFs, could you confirm that MTUQ treat the GFs as velocity if not tags provided? If found part of the code in mtuq.mtuq.greens_tensor.CPS.py as

class GreensTensor(GreensTensorBase):
"""
FK Green's tensor object

Overloads base class with machinery for working with CPS-style
Green's functions

"""
def __init__(self, *args, **kwargs):
    super(GreensTensor, self).__init__(*args, **kwargs)
    
    if 'type:greens' not in self.tags:
        self.tags += ['type:greens']

    if 'type:velocity' not in self.tags:
        self.tags += ['type:velocity']

    if 'units:m' not in self.tags:
        self.tags += ['units:m']

As in the 2nd if condition of init, it treat the default type of GFs as velocity (and it is same in FK.py). How can I tell the function this GFs is in displacement? For the data part, I used:
data = read(path_data, format='sac',
event_id=event_id,
station_id_list=station_id_list,
tags=['units:m', 'type:displacement'])
It seems this tags for data was not pass to GFs. And the processing converted the GFs to the 'displacement'. And if I followed the same way to converting my own synthetics into 'displacement' using numpy.cumsum, the comparison in the first figure became. (Please ignore the time shifts which are caused by my own cuting window).

[rmodrak]

Hi @Jinyin-Hu , Perhaps I could defer to others for those questions, and focus for now on a related issue--

From speaking with Herrmann's and Ammon's students, my understanding is there is a single widespread Green's function convention used by CPS itself and in Green's function libraries at several institutions. The goal with MTUQ is just to support this one convention, which I believe involves velocity and CG units.

[rmodrak]

Moving forward, I think we could have an MTUQ example based on CPS Green's functions that exactly reproduces the current examples based on FK. I believe Robert and I independently had gotten quite close, but not 100 percent of the way there. Let me see if we can resume this.

In addition, better documentation would be very helpful, as you were saying, but can be time consuming when CPS does not include this itself.

[rwalkerlewis]

The solution to this might be to change mtuq.greens_tensor.CPS.py:

    def __init__(self, *args, **kwargs):
        super(GreensTensor, self).__init__(*args, **kwargs)

        if 'type:greens' not in self.tags:
            self.tags += ['type:greens']

        if 'type:velocity' not in self.tags:
            self.tags += ['type:velocity']

        if 'units:m' not in self.tags:
            self.tags += ['units:m']

to similar as in mtuq.greens_tensor.SPECFEM3D.py:

    def __init__(self, *args, **kwargs):
        super(GreensTensor, self).__init__(*args, **kwargs)

        if 'type:greens' not in self.tags:
            self.tags += ['type:greens']

        if 'type:velocity' not in self.tags or 'type:displacement' not in self.tags:
            self.tags += ['type:displacement']

        if 'units:m' not in self.tags:
            self.tags += ['units:m']

(just add an or statement that acknowledges displacement is a type)

[rmodrak]

@rwalkerlewis Hi Robert, very helpful thank you.

Do you recall by chance what CPS outputs? I believe it's velocity by default?

[Jinyin-Hu]

Hi @rmodrak @rwalkerlewis , thanks. I think it is velocity by default in CPS from its manual, and it also has an option for displacement. In my case I choose to output displacement. I think this should solve my problem if I output GFs as velocity.

[rwalkerlewis]

Looking at the CPS source, default output for tpulse96, hpulse96, and gpulse96 are all in ground velocity.

The FK reader assumes velocity as well.

[Jinyin-Hu]

Thank you both.