Hello,
Im builing a program that analyses the ambient noise to infer the location of water. For this im using the H/V method proposed by Nakumura.
The technique originally proposed by Nogoshi and Igarashi (1971), and wide-spread by Nakamura (1989), consists in estimating the ratio between the Fourier amplitude spectra of the horizontal (H) to vertical (V) components of the ambient noise vibrations recorded at one single station.
The computation of the H/V ratio follows different steps (see Figure 1):
- record a 3-component ambient noise signal
- select of the most stationary time windows (e.g., using an anti-triggering algorithm) in order to avoid transient noise
- compute and smoothing of the Fourier amplitude spectra for each time windows
- Average the two horizontal component (using a quadractic mean)
- compute the H/V ratio for each window
- compute the average H/V ratio
Source: http://www.geopsy.org/documentation/geopsy/hv.html
I leave here the code that I been developing to test this method:
import numpy as np
from obspy.core import *
from obspy.arclink import Clientfrom obspy.signal import seisSim
from obspy.xseed import Parserimport matplotlib.pyplot as plt
st = read('http://examples.obspy.org/COP.BHE.DK.2009.050’)
st += read('http://examples.obspy.org/COP.BHN.DK.2009.050’)
st += read('http://examples.obspy.org/COP.BHZ.DK.2009.050’)
e_trace= st[0].copy()
n_trace = st[1].copy()
z_trace= st[2].copy()
client = Client(“webdc.eu”, 18001, user='test@obspy.org’)
for tr in [e_trace, n_trace, z_trace]:
start=tr.stats.starttime
end=start+3600
tr.trim(start,end) #only want 1hour of data for testing purposes
paz=client.getPAZ(tr.stats.network, tr.stats.station, tr.stats.location, tr.stats.channel,start)
df = tr.stats.sampling_rate
tr.data = seisSim(tr.data, df, paz_remove=paz) #desconvolution, converting the data from counts to m/s
#Get the amplitude values from the resulting traces
e_array=e_trace.data
n_array=n_trace.data
z_array=z_trace.data
#Removing the mean from the trace
e=e-np.mean(e)
n=n-np.mean(n)
z=z-np.mean(z)
#Apply the fourier transformation to the arrays
e_fft=np.fft.fft(e)
n_fft=np.fft.fft(n)
z_fft=np.fft.fft(z)
#Smoothing
e=smooth(e_fft)
n=smooth(n_fft)
V=smooth(z_fft)
ndouble=abs(np.power(n,2))
edouble=abs(np.power(e,2))
#Apply the quadratic mean
H=np.sqrt(np.add(edouble,ndouble))
#Devide the two components
HV=np.divide(H,V)
n=e.size
freq = np.fft.fftfreq(n, d=(1/df))
plt.figure()
plt.plot(freq,H)
plt.xlim(0,2)
plt.figure()
plot(freq,HV)
plt.xlim(0,2)
plt.show()
Where smooth is the following function (source):
def smooth(x,window_len=11,window=‘hanning’):
if x.ndim != 1:
raise ValueError, “smooth only accepts 1 dimension arrays.”
if x.size < window_len:
raise ValueError, “Input vector needs to be bigger than window size.”
if window_len<3:
return x
if not window in [‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’]:
raise ValueError, “Window is on of ‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’”
s=np.r_[2x[0]-x[window_len-1::-1],x,2x[-1]-x[-1:-window_len:-1]]
if window == ‘flat’: #moving average
w=np.ones(window_len,‘d’)
else:
w=eval(‘np.’+window+‘(window_len)’)
y=np.convolve(w/w.sum(),s,mode=‘same’)
return y[window_len:-window_len+1]
My problem right now is finding the best way to smooth the waveform data. I researched quit a lot and tested many different methods and functions from the scipy. Im a newbie in obspy and python for that matter. Very few people in my college can help me.Does anybody have a good smoothing algorithm that can give me?
Thanks so much in advance.
Best regards,
Guilherme