roll_stalta: Rolling STA/LTA

In seismicRoll: Fast Rolling Functions for Seismology using 'Rcpp'

Rolling STA/LTA

Description

Fast rolling STA/LTA using C++/Rcpp. Additional performance gains can be achieved by skipping increment values between calculations.

The STA/LTA ratio method is used for automatic detection of seismic signal arrival times.

Usage

roll_stalta(x, n_sta, n_lta, increment = 1)

Arguments

x	an R numeric vector
n_sta	integer STA window size
n_lta	integer LTA window size
increment	integer shift to use when sliding the window to the next location

Details

The roll_stalta function described here does no preprocessing of the incoming data and merely calculates the ratio of the average value in the STA window to the average value in the LTA window. Windows are aligned so that the index is at the left edge of the STA window and at the right edge of the LTA window.

STA(x_i) = \frac{1}{ns}\sum_{j=i}^{i+ns}{x_i}

LTA(x_i) = \frac{1}{nl}\sum_{j=i-nl}^{i}{x_i}

r_i = \frac{STA_i}{LTA_i}

[---------- LTA --------*]........

.......................[*- STA --]

For proper use of this algorithm seismic data should be preprocessed as in the example below with:

• demean, detrend and taper the raw signal

• square the processed signal to get power

With increment=1, this function is equivalent to, eg:

sta <- roll_mean(x,3,1,"left")

lta <- roll_mean(x,30,1,"right")

r <- sta/lta

For increments greater than one, the rolling means above will not align properly, hence the need for a dedicated roll_stalta function.

Values within n_lta-1 of the beginning or n_sta-1 of the end of x are set to NA.

Setting increment to a value greater than one will result in NAs for all skipped-over indices.

Value

A vector of values of the same length as x with each point containing the STA/LTA ratio at that point.

References

First Break Picking

Examples

# Contrived example
x <- rep(c(1,5,3,2,1),each=20)
p <- roll_stalta(x,3,6)
plot(x, pch=17, cex=0.8, ylim=c(0,max(x)),
     main="Test of roll_stalta on artificial data")
points(p,cex=1.5,col='red',type='b')
legend('topleft',
       legend=c('data','STA/LTA'),
       pch=c(17,1),
       col=c('black','red'))

# Real example requiring the 'seismic' package
## Not run: 
require(seismic)
 
# Create a new IrisClient
iris <- new("IrisClient")
  
# Seismic data with a large quake
starttime <- as.POSIXct("2010-02-27 06:30:00", tz="GMT")
endtime <- as.POSIXct("2010-02-27 07:00:00", tz="GMT")
st <- getDataselect(iris,"IU","ANMO","00","BHZ",starttime,endtime)
tr <- st@traces[[1]]
 
# Preprocess the data
x <- DDT(tr)@data
 
# Calculate the first break 'picker'
n_sta <- 3 * tr@stats@sampling_rate
n_lta <- 10 * n_sta
p <- roll_stalta(x^2,n_sta,n_lta)
 
first_break <- which(p == max(p,na.rm=TRUE))

plot(x,type='l',
     main='Test of STA/LTA first break picker on raw seismic data')
abline(v=first_break,col='red')  

## End(Not run)

seismicRoll documentation built on Aug. 22, 2023, 9:13 a.m.