A statistical approach to estimating station biases and error levels

Yael Radzyner, Meirav Galun and Boaz Nadler April 2023

Statistical Model for Station Biases and Error Levels

Consider a seismic network of N monitoring stations.
Suppose a seismic event j of unknown (true) magnitude Y_j was detected by a subset of stations S_j. Denote by X_i,j the reported magnitudes by these stations, i ∈ S_j.
The currently employed standard method to compute a network estimate of the event magnitude is by a simple averaging (excluding outliers) of the reported magnitudes

Y_net(event j) =

\frac{1}{| S j |}

∑_{i∈ S_j} X_i,j
The above formula implicitly assumes that stations have very small (if any) systematic errors (biases) and that all stations have the same (or very similar) error levels, i.e., the random errors in their estimated event magnitudes all have the same standard deviations. Namely the above averaging is consistent with the following model,

X_i,j = Y_j + σ ξ_i,j
where σ is the error level at all stations and ξ_i,j are random variables with zero mean and unit variance.
As we illustrate in our manuscript, using a large dataset from the REB, empirical data is inconsistent with the above model.

In our work we instead consider the following model. Each station i has two unknown parameters: a station bias b_i and a station error level σ_i.
We assume that for an event j of unknown magnitude Y_j, for each station i that reported a magnitude for the event, the follows holds

X_i,j = Y_j + b_i + σ_i ξ_i,j
Given a large collection of reported station magnitudes X_i,j over a large set of events j the goal is to estimate the vector of station biases b = (b₁,···,b_N) and the vector of station error levels σ = (σ₁,···,σ_N)

Inference Method

The key idea in our approach to estimate the vectors b and σ is to construct quantitites that do not depend on the unknown event magnitudes Y_j, which may be viewed as nuisance parameters.
To this end we consider for any event j and for any pair of stations i,k that detected this event, the following quantity
X_i,j - X_k,j = b_i - b_k + σ_iξ_i,j - σ_kξ_k,j
Conveniently in the above quantity the unknown magnitude Y_j has been cancelled out.
Error Level Estimation:
Let v = (σ₁²,···, σ_N²). Based on the above, we propose to estimate v by minimizing the following functional
T(v) = ∑_i≠k |M_i,k | · ( V_i,k - v_i - v_k )² where V_i,k = Var( X_i,j - X_k,j ) is the empirical variance of the above differences, over all events j jointly detected by stations i,k and and M_i,k is the total number of such events jointly detected by stations i,k In the above the summation is only over pairs of stations for which |M_i,k| ≥ 2. The above is a quadratic objective, and its solution is given by the solution of an N×N system of linear equations.

Stations' Biases Estimation:

Similarly, to estimate the station biases we first compute the following empirical means D_i,k =

\frac{1}{| M i,k |}

· ∑_{j ∈ M_i,k} ( X_i,j - X_k,j )
We then propose to estimate the station biases by minimizing the following objective G (b) =

\frac{|M i,k |}{σ i 2 + σ k 2}

· ∑_i≠k ( D_i,k - b_i + b_k )²
Since the true error levels are unknown, in the above we replace σ_i² by v_i, estimated by the solution of the linear system described above. Estimating the station biases also amounts to solving an N× N system of linear equations. Importantly, this system is rank defficient. Assuming the graph of stations having at least 2 jointly detected events is connected the rank defficiently is one. We thus find the solution of minimal norm, and then remove the remaining single degree of freedom by the normalization ∑_i b_i = 0.
Python Code and Demo:
The following Python script generates random data from a network of N monitoring stations and then estimates the station error levels and biases as described above.
demo_seismic_station_bias_error_levels.py
Download Python source