obspy.signal.array_analysis.array_rotation_strain

array_rotation_strain(subarray, ts1, ts2, ts3, vp, vs, array_coords, sigmau)

This routine calculates the best-fitting rigid body rotation and uniform strain as functions of time, and their formal errors, given three-component ground motion time series recorded on a seismic array. The theory implemented herein is presented in the papers [Spudich1995], (abbreviated S95 herein) [Spudich2008] (SF08) and [Spudich2009] (SF09).

This is a translation of the Matlab Code presented in (SF09) with small changes in details only. Output has been checked to be the same as the original Matlab Code.

Note

ts_ below means “time series”

Parameters:

• vp (float) – P wave speed in the soil under the array (km/s)

• vs (float) – S wave speed in the soil under the array Note - vp and vs may be any unit (e.g. miles/week), and this unit need not be related to the units of the station coordinates or ground motions, but the units of vp and vs must be the SAME because only their ratio is used.

• array_coords (numpy.ndarray) – array of dimension na x 3, where na is the number of stations in the array. array_coords[i,j], i in arange(na), j in arange(3) is j coordinate of station i. units of array_coords may be anything, but see the “Discussion of input and output units” above. The origin of coordinates is arbitrary and does not affect the calculated strains and rotations. Stations may be entered in any order.

• ts1 (numpy.ndarray) – array of x1-component seismograms, dimension nt x na. ts1[j,k], j in arange(nt), k in arange(na) contains the k’th time sample of the x1 component ground motion at station k. NOTE that the seismogram in column k must correspond to the station whose coordinates are in row k of in.array_coords. nt is the number of time samples in the seismograms. Seismograms may be displacement, velocity, acceleration, jerk, etc. See the “Discussion of input and output units” below.

• ts2 (numpy.ndarray) – same as ts1, but for the x2 component of motion.

• ts3 (numpy.ndarray) – same as ts1, but for the x3 (UP or DOWN) component of motion.

• sigmau (float or numpy.ndarray) –

  standard deviation (NOT VARIANCE) of ground noise, corresponds to sigma-sub-u in S95 lines above eqn (A5). NOTE: This may be entered as a scalar, vector, or matrix!

  • If sigmau is a scalar, it will be used for all components of all stations.

  • If sigmau is a 1D array of length na, sigmau[i] will be the noise assigned to all components of the station corresponding to array_coords[i,:]

  • If sigmau is a 2D array of dimension na x 3, then sigmau[i,j] is used as the noise of station i, component j.

  In all cases, this routine assumes that the noise covariance between different stations and/or components is zero.

• subarray (numpy.ndarray) – NumPy array of subarray stations to use. I.e. if subarray = array([1, 4, 10]), then only rows 1, 4, and 10 of array_coords will be used, and only ground motion time series in the first, fourth, and tenth columns of ts1 will be used. Nplus1 is the number of elements in the subarray vector, and N is set to Nplus1 - 1. To use all stations in the array, set in.subarray = arange(na), where na is the total number of stations in the array (equal to the number of rows of in.array_coords. Sequence of stations in the subarray vector is unimportant; i.e. subarray = array([1, 4, 10]) will yield essentially the same rotations and strains as subarray = array([10, 4, 1]). “Essentially” because permuting subarray sequence changes the d vector, yielding a slightly different numerical result.

Returns:

Dictionary with fields:

A: (array, dimension 3N x 6)

data mapping matrix ‘A’ of S95(A4)

g: (array, dimension 6 x 3N)

generalized inverse matrix relating ptilde and data vector, in S95(A5)

ce: (4 x 4)

covariance matrix of the 4 independent strain tensor elements e11, e21, e22, e33

ts_d: (array, length nt)

dilatation (trace of the 3x3 strain tensor) as a function of time

sigmad: (scalar)

standard deviation of dilatation

ts_dh: (array, length nt)

horizontal dilatation (also known as areal strain) (eEE+eNN) as a function of time

sigmadh: (scalar)

standard deviation of horizontal dilatation (areal strain)

ts_e: (array, dimension nt x 3 x 3)

strain tensor

ts_s: (array, length nt)

maximum strain ( .5*(max eigval of e - min eigval of e) as a function of time, where e is the 3x3 strain tensor

cgamma: (4 x 4)

covariance matrix of the 4 independent shear strain tensor elements g11, g12, g22, g33 (includes full covariance effects). gamma is traceless part of e.

ts_sh: (array, length nt)

maximum horizontal strain ( .5*(max eigval of eh - min eigval of eh) as a function of time, where eh is e(1:2,1:2)

cgammah: (3 x 3)

covariance matrix of the 3 independent horizontal shear strain tensor elements gamma11, gamma12, gamma22 gamma is traceless part of e.

ts_wmag: (array, length nt)

total rotation angle (radians) as a function of time. I.e. if the rotation vector at the j’th time step is w = array([w1, w2, w3]), then ts_wmag[j] = sqrt(sum(w**2)) positive for right-handed rotation

cw: (3 x 3)

covariance matrix of the 3 independent rotation tensor elements w21, w31, w32

ts_w1: (array, length nt)

rotation (rad) about the x1 axis, positive for right-handed rotation

sigmaw1: (scalar)

standard deviation of the ts_w1 (sigma-omega-1 in SF08)

ts_w2: (array, length nt)

rotation (rad) about the x2 axis, positive for right-handed rotation

sigmaw2: (scalar)

standard deviation of ts_w2 (sigma-omega-2 in SF08)

ts_w3: (array, length nt)

”torsion”, rotation (rad) about a vertical up or down axis, i.e. x3, positive for right-handed rotation

sigmaw3: (scalar)

standard deviation of the torsion (sigma-omega-3 in SF08)

ts_tilt: (array, length nt)

tilt (rad) (rotation about a horizontal axis, positive for right handed rotation) as a function of time tilt = sqrt( w1^2 + w2^2)

sigmat: (scalar)

standard deviation of the tilt (not defined in SF08, From Papoulis (1965, p. 195, example 7.8))

ts_data: (array, shape (nt x 3N))

time series of the observed displacement differences, which are the di in S95 eqn A1

ts_pred: (array, shape (nt x 3N))

time series of the fitted model’s predicted displacement difference Note that the fitted model displacement differences correspond to linalg.dot(A, ptilde), where A is the big matrix in S95 eqn A4 and ptilde is S95 eqn A5

ts_misfit: (array, shape (nt x 3N))

time series of the residuals (fitted model displacement differences minus observed displacement differences). Note that the fitted model displacement differences correspond to linalg.dot(A, ptilde), where A is the big matrix in S95 eqn A4 and ptilde is S95 eqn A5

ts_m: (array, length nt)

Time series of M, misfit ratio of S95, p. 688

ts_ptilde: (array, shape (nt x 6))

solution vector p-tilde (from S95 eqn A5) as a function of time

cp: (6 x 6)

solution covariance matrix defined in SF08

Warnings

This routine does not check to verify that your array is small enough to conform to the assumption that the array aperture is less than 1/4 of the shortest seismic wavelength in the data. See SF08 for a discussion of this assumption.

This code assumes that ts1[j,:], ts2[j,:], and ts3[j,:] are all sampled SIMULTANEOUSLY.

Notes

1. Note On Specifying Input Array And Selecting Subarrays

   This routine allows the user to input the coordinates and ground motion time series of all stations in a seismic array having na stations and the user may select for analysis a subarray of Nplus1 <= na stations.

2. Discussion Of Physical Units Of Input And Output

   If the input seismograms are in units of displacement, the output strains and rotations will be in units of strain (unitless) and angle (radians). If the input seismograms are in units of velocity, the output will be strain rate (units = 1/s) and rotation rate (rad/s). Higher temporal derivative inputs yield higher temporal derivative outputs.

   Input units of the array station coordinates must match the spatial units of the seismograms. For example, if the input seismograms are in units of m/s^2, array coordinates must be entered in m.

3. Note On Coordinate System

   This routine assumes x1-x2-x3 is a RIGHT handed orthogonal coordinate system. x3 must point either UP or DOWN.