Abnormal Rayleigh wave ellipticity points with some velocity models

[chengsong2333]

Dear Prof. Luu, I am a PhD student at Seoul National University, I really appreciate your great code. I enconter some abnormal ellipticity points compared with Geopsy(gpell) results:

You can load and plot the model by

RWE_disba_geopsy_model.txt
ell_data_freq.txt

import numpy as np
from disba import Ellipticity, EllipticitySensitivity
np.loadtxt('RWE_disba_geopsy_model.txt')
f_hob, _, _ = np.loadtxt('ell_data_freq.txt').T
t = 1/f_hob[::-1]
ell = Ellipticity(*velocity_model_i.T, algorithm='dunkin',dc=0.0001)
rel = ell(t, mode=0)
plt.semilogx(1/t, np.abs(rel[1]), label='disba')
plt.title('RWE')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.legend()

I believe in your code, it is just some abnormal points under some extreme stiuations, I would appreciate if you could fix the code to smooth these abnormal points

[keurfonluu]

Hi @chengsong2333,

Are you using the same frequency axis for disba and geopsy? If you compute the ellipticity between 0.49181134 and 0.50655995 seconds where you observe "strange" ellipticity values, it doesn't look like abnormal points.

Rather, looking at CPS's documentation about ellipticity, it seems to be due to sharp contrasts in velocities in your model.

However, I am not sure why geopsy doesn't exhibit such behavior.

[chengsong2333]

Dear Prof. Luu, Thanks so much for your timely reply, I read carefull about the description of elltpticity, I am actually also aware of the meaning of peaks of ellipticity because I did inversion of ellipticity using your code :). Following your period between 0.49181134 and 0.50655995, I generate 1000 points using Geopsy, the output file is mars6.txt, and I load the frequency/period to compute ellipticity using your code:

mars6.txt

import numpy as np
from disba import Ellipticity, EllipticitySensitivity
velocity_model_i = np.loadtxt('RWE_disba_geopsy_model.txt')
freq, ell_gpell = np.loadtxt('mars6.txt', skiprows = 6).T
t = 1/freq[::-1]
ell = Ellipticity(*velocity_model_i.T)
rel = ell(t, mode=0)
plt.semilogx(1/t, np.abs(rel[1]), label='disba')
plt.semilogx(freq, np.abs(ell_gpell), '--', label='gpell')
plt.title('RWE')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.legend()

But it does not mean Geopsy(gpell) is 100% correct. So I try to compute forward calculation of ellipticity using CPS, but I fail to predict lower frequency parts, I have attached the zipped file, I calculate ellipticity with these command:

Luu_file.zip

sprep96 -HR 0 -HS 10 -M mars2.d -d dfile_mars -PARR period_file_hob -NMOD 1 -L -R
sdisp96
sdpsrf96 -R -XMIN 0 -XMAX 4 -YMIN 3.0 -YMAX 5.0
sregn96
sdpegn96 -R -C -ASC

And plot:

import numpy as np
from disba import Ellipticity, EllipticitySensitivity
velocity_model_i = np.loadtxt('RWE_disba_geopsy_model.txt')
f_hob, _, _ = np.loadtxt('ell_data_freq.txt').T
cps_data = np.loadtxt('SREGN.ASC', skiprows=1)
ell_cps = cps_data[:,-1]
f_cps = cps_data[:,3]
t = 1/f_cps [::-1]

ell = Ellipticity(*velocity_model_i.T, algorithm='dunkin',dc=0.0001)
rel = ell(t, mode=0)

plt.semilogx(1/t, np.abs(rel[1]), label='disba')
plt.semilogx(f_cps, np.abs(ell_cps), label='cps')
plt.title('RWE')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.legend()

That is the information I can provide, I have no idea what is going on, and I would really appreciate if you could figure it out