I’ve been playing around with different ideas to test for cross channel interference or interaction on the RS&B and finally got my head around the Cross_Correlation functions in Obspy which ATM seems by far the best way to do it compared to what I’ve tried.
This is of particular interest to me to identify and eliminate seismic sources of infrasound signal when trying to detect meteors by infrasound.
Attached is the code I’ve used FYI and use if interested. I’m sure it could be polished up and written in a more elegant manner.
from obspy.clients.fdsn import Client from obspy.core import UTCDateTime, Stream from obspy.signal import filter, cross_correlation import matplotlib.pyplot as plt rs = Client('RASPISHAKE') xc = cross_correlation # set the station name and download the response information stn = 'R21C0' # your station name inv = rs.get_stations(network='AM', station=stn, level='RESP') print (inv) # set data start/end times start = UTCDateTime(2022, 6, 6, 23, 17, 16.7) # (YYYY, m, d, H, M, S.ffffff) end = start + 1 # end = start + duration in seconds (S.ffffff) # set the FDSN server location and channel names channels = ['EHZ' , 'HDF'] # ENx = accelerometer channels; EHx or SHZ = geophone channels # get waveforms and put them all into one Stream stream = Stream() for ch in channels: trace = rs.get_waveforms('AM', stn, '00', ch, start, end) stream += trace stream.merge(method=0, fill_value='latest') stream.detrend(type='demean') print(stream) #get instrument response inv = rs.get_stations(network="AM", station=stn, level="RESP") # setup figure and subplots fig = plt.figure(figsize=(12,12)) ax1 = fig.add_subplot(4,1,1) ax2 = fig.add_subplot(4,1,2) ax3 = fig.add_subplot(4,1,3) ax4 = fig.add_subplot(4,1,4) fig.suptitle("AM."+stn+" RS&B Cross Channel Correlation Report: "+start.strftime('%d/%m/%Y %H:%M:%S.%f UTC')) # plot raw traces ax1.plot(stream.times(reftime=start), stream.data, lw=1, color='g') ax2.plot(stream.times(reftime=start), stream.data, lw=1) # calculate cross correlation s = 10 # shift in data points each side to test cc = xc.correlate(stream, stream, s) print(cc) # print correlations for each shift value shift, value = xc.xcorr_max(cc) #find the max correlation and shift print(shift, round(value,3)) # get the limits of the y axis so text can be consistently placed ax2b, ax2t = ax2.get_ylim() #add the correlation to the plot ax2.text(0, ax2t*0.8, 'Raw Shift = '+str(shift)+'/100 s. Correlation Coefficient ='+str(round(value,3))) # set up filter filt = [0.3, 0.5, 10, 12] #edit filter values to suit # remove instrument response and add new traces to the stream stream += stream.remove_response(inventory=inv,pre_filt=filt,output="VEL", water_level=60) stream += stream.remove_response(inventory=inv,pre_filt=filt,output="DEF", water_level=60) # plot filtered traces ax3.plot(stream.times(reftime=start), stream.data, lw=1, color='r') ax4.plot(stream.times(reftime=start), stream.data, lw=1, color='b') #calculate cross correlation of filtered traces cc1 = xc.correlate(stream, stream, s) print(cc1) # print correlations for each shift value shift1, value1 = xc.xcorr_max(cc1) #find the max correlation and shift print(shift1, round(value1,3)) # get ax4 y axis limits so text can be consistently placed ax4b, ax4t = ax4.get_ylim() #add the correlation to the plot ax4.text(0, ax4t*0.8, 'Filtered Shift = '+str(shift1)+'/100 s. Correlation Coefficient ='+str(round(value1,3))) # setup some plot details ax1.set_ylabel("EHZ Raw Counts") ax2.set_ylabel("HDF Raw Counts") ax3.set_ylabel("EHZ Velocity (m/s) ("+str(filt)+" to "+str(filt)+"Hz)") ax4.set_ylabel("HDF Pascals (Pa) ("+str(filt)+" to "+str(filt)+"Hz)") # save the final figure #plt.savefig(stn+' RSandB Corr'+start.strftime('%Y%m%d %H%M%S UTC')) #comment this line out till figure is final # show the final figure plt.show()
This is what the output looks like:
In this case, there’s a 48% correlation between the filtered EHZ and HDF channels and as the HDF signal is leading the EHZ channel by 0.05s it appears to be the HDF channel influencing the EHZ channel.
This obviously works best for short sections of signal. Long sections are likely to reduce the correlation as random noise becomes more significant.