[obspy-users] Very long period data PPSD elevated noise

Tobias Megies megies at geophysik.uni-muenchen.de
Fri Jan 16 11:39:03 CET 2015


Hi Adam,

good to hear that it's cleared up.
If you think you have a way of processing that's better for the non-flat
parts or maybe if you even can help improve testing / cross checking the
PPSD routine, please don't hesitate to open an issue or pull request on
github.

best,
Tobias


On 01/12/2015 08:25 PM, Ringler, Adam wrote:
> Tobais and Lion,
> 
> Thanks for the help.  I just confirmed the elevated vlp noise is coming
> from the water level.  I switched it to 30 and it looked slightly
> elevated (however a definite improvement) and at 10 it went below the
> NLNM at about 1000 seconds period.  To me this suggests it is probably
> better to remove the response from the power spectra and not on the time
> series for estimating pdfs at very long periods.  Using paztofreqresp
> seems to solve this issue.
> 
> Thanks again,
> Ada,
> 
> 
> 
> 
> 
> On Thu, Jan 8, 2015 at 12:56 PM, Tobias Megies
> <megies at geophysik.uni-muenchen.de
> <mailto:megies at geophysik.uni-muenchen.de>> wrote:
> 
>     Hi Adam,
> 
>     looking at the response
>     (http://service.iris.edu/irisws/evalresp/1/query?net=IU&sta=TUC&loc=00&cha=LHZ&time=2013-04-21T10:00:00&output=plot),
>     I suspect that the "Trace.simulate()" call inside PPSD is responsible.
>     It looks like it's using obspy.signal.invsim.seisSim's default values in
>     the PAZ/overall sensitivity/waterlevel deconvolution. seisSim is one of
>     signal's oldest routines and back then a default value of 600 for the
>     waterlevel got implemented, which is much too high from my experience
>     and might cause noise to blow up in your example. I'm using 10-30
>     usually.
> 
>     Maybe you can try to check if hardcoding a more conservative value in
>     spectral_estimation
>     (https://github.com/obspy/obspy/blob/master/obspy/signal/spectral_estimation.py#L640-L641)
>     helps, or you could try to exchange that call with the "seedresp"-option
>     in simulate, providing a RESP file or Parser object (from dateless) and
>     going via evalresp (specifying a pre filtering).
> 
>     hope it helps,
>     Tobias
> 
>     On 12/30/2014 10:31 PM, Ringler, Adam wrote:
>     > Hello Users,
>     >
>     > My apology if this is the wrong place to send this.
>     >
>     > I noticed a potential artifact in using the PPSD function at very long
>     > periods.  It appears that when you deconvolve the data, convert
>     the data
>     > to acceleration and then estimate the psd you might be elevating the
>     > long period spectra.
>     >
>     > Attached is an example of using the PPSD function from
>     > obspy.spectral_estimation (example.png).  I have also included an
>     > example of first estimating the power spectra and then converting to
>     > acceleration by dividing out the response and using a factor of
>     omega to
>     > switch to acceleration.
>     >
>     > As you can see the second plot rides a few dB above the NLNM while
>     > (which I would expect) the first example shows a "hump" with elevated
>     > noise levels past a couple 100 seconds period (which I would not
>     expect).
>     >
>     > Below are more details.
>     >
>     > Best,
>     > Adam
>     >
>     > I added the following code to a local version of
>     spectral_estimation.py
>     > (line 663) to make this change
>     >
>     >
>     ##########################################################################
>     >
>     >         resp =
>     >
>     pazToFreqResp(paz['poles'],paz['zeros'],paz['gain']*paz['sensitivity'],self.sampling_rate,
>     > \
>     >    nfft=self.nfft)
>     >         resp = resp[1:]
>     >         resp = resp[::-1]
>     >         respval = np.absolute(resp*np.conjugate(resp))
>     >         w = 2.0*math.pi*_freq[1:]
>     >         w = w[::-1]
>     >         spec =  (w**2) * spec / respval
>     >         # go to dB
>     >         spec = np.log10(spec)
>     >         spec *= 10
>     >
>     >
>     > I also commented out the tr.simulate and the gradient (lines 634 to
>     > 648).  I am using w to convert to acceleration and am using respval to
>     > remove the response. I named this PPSD2 and then ran the following
>     > comparison:
>     >
>     > #!/usr/bin/env python
>     >
>     > import sys
>     > import os
>     > sys.path.append(os.getcwd())
>     >
>     > import spectral_estimation
>     >
>     > from obspy.signal import PPSD, psd
>     > from obspy.core import read, UTCDateTime, Stream
>     > from obspy.xseed import Parser
>     >
>     > #Read in test data
>     > st = read('tucCHECK.seed')
>     >
>     > #Read in dataless and get paz for test data
>     > parser = Parser("/APPS/metadata/SEED/IU.dataless")
>     > paz = parser.getPAZ(st[0].id,datetime=st[0].stats.starttime)
>     >
>     >
>     > #Make a PPSD with a really long window
>     > ppsd = PPSD(st[0].stats, paz,ppsd_length=50000,overlap=0.0)
>     > ppsd.add(st)
>     > ppsd.plot(filename='example.png',period_lim=(2.0,3000))
>     >
>     > #Compare with the frequency domain correction
>     > ppsd2 = spectral_estimation.PPSD2(st[0].stats,
>     > paz,ppsd_length=50000,overlap=0.0)
>     > ppsd2.add(st)
>     > ppsd2.plot(filename='example2.png',period_lim=(2.0,3000))
>     >
>     >
>     > _______________________________________________
>     > obspy-users mailing list
>     > obspy-users at lists.swapbytes.de <mailto:obspy-users at lists.swapbytes.de>
>     > http://lists.swapbytes.de/mailman/listinfo/obspy-users
>     >
> 
>     --
>     Dipl.-Geophys. Tobias Megies
> 
>     Geophysikalisches Observatorium
>     Ludwigshöhe 8
>     82256 Fürstenfeldbruck
> 
>     Ludwig-Maximilians-Universität
>     Department für Geo- und Umweltwissenschaften
>     Sektion Geophysik
>     Theresienstrasse 41/IV
>     80333 München
> 
>     Tel: +49 (0) 89 2180-73981
>          +49 (0) 89 2180-4326
>     Mail: tobias.megies at geophysik.uni-muenchen.de
>     <mailto:tobias.megies at geophysik.uni-muenchen.de>
>     _______________________________________________
>     obspy-users mailing list
>     obspy-users at lists.swapbytes.de <mailto:obspy-users at lists.swapbytes.de>
>     http://lists.swapbytes.de/mailman/listinfo/obspy-users
> 
> 
> 
> 
> _______________________________________________
> obspy-users mailing list
> obspy-users at lists.swapbytes.de
> http://lists.swapbytes.de/mailman/listinfo/obspy-users
> 

-- 
Dipl.-Geophys. Tobias Megies

Geophysikalisches Observatorium
Ludwigshöhe 8
82256 Fürstenfeldbruck

Ludwig-Maximilians-Universität
Department für Geo- und Umweltwissenschaften
Sektion Geophysik
Theresienstrasse 41/IV
80333 München

Tel: +49 (0) 89 2180-73981
     +49 (0) 89 2180-4326
Mail: tobias.megies at geophysik.uni-muenchen.de


More information about the obspy-users mailing list