Locate Earthquake with UTM projection

Description

Non-linear hypocenter location with UTM geographical projection. Used for locating earthquakes in local or regional settings.

Usage

1
2
3
4
5
XYlocate(Ldat, EQ, vel, maxITER = 10, distwt = 10,
lambdareg = 100, FIXZ
= FALSE, REG = TRUE, WTS = TRUE, STOPPING = TRUE,
RESMAX = c(.4,.5), tolx = 0.005, toly = 0.005,
 tolz = 0.01, PLOT = FALSE)

Arguments

Ldat

list, must inlude: x,y,err, sec, cor (see details)

EQ

list, must inlude: x,y,z, t

vel

list, 1D velocity structure

maxITER

Maximum number of iterations

distwt

distance weighting factor

lambdareg

regularization parameter for damping

FIXZ

logical, TRUE = fix depth, i.e. only calculate x,y,t

REG

logical, TRUE=use regularization

WTS

logical, TRUE==use weighting

STOPPING

logical, TRUE=use stopping criteria

RESMAX

vector, residual max for P and S, default=c(4,5)

tolx

numeric, tolerance in km in x direction

toly

numeric, tolerance in km in y direction

tolz

numeric, tolerance in km in z direction

PLOT

logical, plot results during iterations

Details

Input pick list must have at x,y,z, sec, cor, err elements for each station. If no station correction is available it is set to zero. If no uncertainty (err) is available, it is set to 0.05 sec. Each station must have a finite x-y coordinate and arrival time in seconds. Events are located relative to the minute.

Routine uses the svd in a sequence of linear inversions to estimate the nonlinear location.

Value

List:

EQ

list, Earthquake hypocenter and time

its

number of iterations

rms

rms residual

wrms

wheighted rms residual

used

vector, index of used equations

guesses

list of x,y,z,t intermediate locations when converging

Note

This routine should be called by a wrapper (Vlocate) that applies the algorithm several times and changes parameters based on the quality.

If RESMAX is used and the robust approach yields fewer than 4 equations, the best (smallest) four residuals will be used to determiine the event location.

Author(s)

Jonathan M. Lees<jonathan.lees@unc.edu>

See Also

Vlocate

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
## Not run: 
library(RSEIS)
data(GH)

g1 = GH$pickfile
data(VELMOD1D)
vel= VELMOD1D



 w1 = which(!is.na(g1$STAS$lat))
         sec = g1$STAS$sec[w1]

         N = length(sec)
         Ldat =    list(
           name = g1$STAS$name[w1],
           sec = g1$STAS$sec[w1],
           phase = g1$STAS$phase[w1],
           lat=g1$STAS$lat[w1],
           lon = g1$STAS$lon[w1],
           z = g1$STAS$z[w1],
           err= g1$STAS$err[w1],
           yr = rep(g1$LOC$yr , times=N),
           jd = rep(g1$LOC$jd, times=N),
           mo = rep(g1$LOC$mo, times=N),
           dom = rep(g1$LOC$dom, times=N),
           hr =rep( g1$LOC$hr, times=N),
           mi = rep(g1$LOC$mi, times=N) )

 MLAT = median(Ldat$lat)
    MLON = median(Ldat$lon)
    
    proj = GEOmap::setPROJ(type=2, LAT0=MLAT, LON0=MLON)

####   get station X-Y values in km
    XY = GEOmap::GLOB.XY(Ldat$lat, Ldat$lon, proj)
###   add to Ldat list
    Ldat$x = XY$x
    Ldat$y = XY$y
       wstart = which.min(Ldat$sec)



EQ = list(x=XY$x[wstart], y=XY$y[wstart], z=6, t=Ldat$sec[wstart] )

      
 maxITER = 7
###print(EQ)
    AQ = XYlocate(Ldat,EQ,vel, 
      maxITER = maxITER,
      distwt = 1,
      lambdareg =10 ,
      FIXZ = FALSE,
      REG = TRUE,
      WTS = TRUE,
      STOPPING = TRUE,
      RESMAX = c(0.1,0.1),
      tolx =   0.001,
      toly = 0.001 ,
      tolz = 0.5, PLOT=FALSE)

########  update the new location

AXY = GEOmap::XY.GLOB(AQ$EQ$x, AQ$EQ$y, proj)
AQ$EQ$lat = AXY$lat
AQ$EQ$lon = AXY$lon
if(AQ$EQ$lon>180) { AQ$EQ$lon = AQ$EQ$lon-360 }


plot(c(Ldat$x, AQ$EQ$x) , c(Ldat$y,AQ$EQ$y), type='n' , xlab="km",
ylab="km" )

points(Ldat$x, Ldat$y, pch=6)

points(AQ$EQ$x, AQ$EQ$y, pch=8, col='red')

points(EQ$x, EQ$y, pch=4, col='blue')


legend("topright", pch=c(8,4, 6), col=c("red", "blue", "black"),
 legend=c("Final location", "Initial guess", "Station"))


print(AQ)

EQ$x = 10
EQ$y = 2


## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.