Help on developing the H/V method using Obspy and Numpy

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 Client

from obspy.signal import seisSim
from obspy.xseed import Parser

import 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

Hi,

without having looked to much in the details.. maybe this here might be
of interest: https://github.com/krischer/HtoV-Toolbox ?

best,
Tobias

Hey Guilherme,

for the spectral smoothing you will probably want to use the so called Konno-Ohmachi smoothing introduced in their paper:

Konno, K., & Ohmachi, T. (1998). Ground-Motion Characteristics Estimated from Spectral Ratio between Horizontal and Vertical Vomponents of Microtremor. Bulletin of the Seismological Society of America, 88(1), 228–241.

It essentially it a smoothing window that looks to have a constant width on a logarithmic axis. ObsPy actually has an implementation of it:
http://docs.obspy.org/packages/autogen/obspy.signal.konnoohmachismoothing.html

Quite a while ago I by chance wrote something very similar utilising ObsPy:
https://github.com/krischer/HtoV-Toolbox

It was written while ObsPy was still in its infancy so you most likely will have to adjust some things to make it work. It does pretty much all the steps you described so it might be a good starting point. Additionally it performs the spectral calculations using a multi taper estimation.

Cheers!

Lion