Weighing Geophysical Data With Trans‐Dimensional Algorithms: An Earthquake Location Case Study

In geophysical inverse problems, the distribution of physical properties in an Earth model is inferred from a set of measured data. A necessary step is to select data that are best suited to the problem at hand. This step is performed ahead of solving the inverse problem, generally on the basis of expert knowledge. However, expert‐opinion can introduce bias based on pre‐conceptions. Here we apply a trans‐dimensional algorithm to automatically weigh data on the basis of how consistent they are with the fundamental hypotheses made to solve the inverse problem. We demonstrate this approach by inverting arrival times for the location of a seismic source in an elastic half‐space, assuming a point‐source and uniform weights in concentric shells. The key advantage is that the data do no longer need to be selected by an expert, but they are assigned varying weights during the inversion procedure.


10.1029/2023GL102983
2 of 10 Here we propose a novel approach to incorporate the choice of weights for the data in the inversion process (or, more precisely, the variance of data noise).Our approach is based on trans-dimensional Markov chain Monte Carlo (McMC) sampling (Piana Agostinetti et al., 2021;Piana Agostinetti & Sgattoni, 2021) and works by proposing and accepting/rejecting data weighing schemes following the Metropolis algorithm (Sambridge & Mosegaard, 2002) where the data weighing schemes have a variable number of parameters (Malinverno, 2002;Sambridge et al., 2006).The complexity of the weighing scheme is dictated by the data themselves, rather than by user-defined choices made during pre-processing.The assigned weights depend on how closely different data match the fundamental assumptions made in solving the inverse problem.
The novel approach is developed in a fully consistent Bayesian framework, which means that Posterior Probability Distributions (PPD) for the weights are retrieved at the end of the McMC sampling.From a Bayesian point of view, what we call "weights" are data variances (i.e., data uncertainties).In Hierarchical Bayes, estimated data uncertainties is known to incorporate both measurement errors, modeling errors and other sources of uncertainty (e.g., Bodin et al., 2012).Along these lines, our approach can be considered a step further in the context of the algorithms for error covariance estimation (Bodin et al., 2012;Dettmer & Dosso, 2012;Galetti et al., 2016;Kolb & Lekić, 2014;Mustać & Tkalčić, 2015;Piana Agostinetti et al., 2021;Piana Agostinetti & Sgattoni, 2021;Steininger et al., 2013;Xiang et al., 2018).In particular, Ghalenoei et al. (2022) developed an approach, where two separate trans-D samplings are applied to both error and geophysical models.
We test our approach in the geophysical inverse problem of locating a seismic point source using P-and S-wave arrival times recorded by sensors in a seismic network.In this inverse problem, data are generally downweighted with distance of the sensor from the seismic source or are removed in a pre-processing step if the sensors are farther than a chosen distance from the source.In our novel approach, we define a set of spherical shells centered on the source (Figure 1a).All sensors within a shell are assigned the same weight (Figure 1b), but more complex weight assignments can be made (e.g., weights that vary linearly with distance from the source within each shell; see Figure 1c).The number of shells, their radii and weights are unknown, and will be defined by the McMC sampling.The stations that receive the largest weights will be those that measure arrival times consistent with the fundamental assumptions made in the inverse problem (namely, a point-wise seismic source and constant P-and S-wave velocities in the rock volume).
Our natural laboratory is Kiirunavaara mine (Sweden), a 6 km-long active mine with more than 200 seismic sensors in a 3D configuration that spans along the exploited rock volume (Dineva et al., 2022).Given such an extensive seismic network, events can be well located in three dimensions.We selected two seismic events.The first is a man-made blast, used to calibrate the seismic network (Figure 1d).The actual location of this seismic source is known within <1 m and can be immediately used to evaluate our results.The second is a natural M w 4.2 multi-phase seismic event that occurred on 18 May 2020 (Dineva et al., 2022) and it was recorded on all working sensors in the mine (Figure 1e).Our experiment is structured as follows.We first compute a reference solution for the calibration blast by applying a standard McMC algorithm (Riva et al., 2023).In this reference solution, we do not use our novel approach, but we solve for the source location by removing data from sensors at a range of distances from a preliminary location of the seismic source, as done in standard seismological workflows in mines.This is intended to simulate a range of possible expert opinions on the distance threshold for data selection (here we assume that the hypocentral distance is of such utmost importance that observational quality differences may be neglected, which is not the case in crustal studies).We then apply our novel approach to the complete data set for the calibration blast and compare the results with those in the reference solution.Finally, we apply our methodology to the natural seismic event, to observe its performance for a larger and more complex event.

Data and Methods
In our study, we performed three different experiments based on two data sets of seismic measurements.Raw seismic waveforms were recorded in the Kiirunavaara mine (Sweden) by a seismic network (253 seismic sensors spanning the 3D volume around the ore body (Dineva et al., 2022)).From the raw seismic waveforms, P-and S-wave arrival times were automatically retrieved and manually revised by the local seismic system provider (Institute of Mine Seismology, IMS).These arrival times are the data used to solve the geophysical inverse problem of locating the seismic source.In the first data set, the seismic source is a man-made blast used to calibrate the seismic system.The position of the seismic source is known with an accuracy less than 1 m.This data set can be used as the input of a robust field test, as the source location is known with an accuracy higher than the half-width of the characteristic 10.1029/2023GL102983 3 of 10

10.1029/2023GL102983
4 of 10 wavelength of the P-wave (given a blast corner frequency close to 1,000 Hz, if we assume Vp = 6,600 m/s, e.g., the half-wavelength is 3.3 m).Due to the limited amount of explosive used (5 kg), seismic waves have been clearly recorded at 57 seismic sensors only (with 57 P-and 9 S-wave arrival times), which translates in a minimum and maximum sensor distance from the source of 83 and 710 m, respectively.The second data set is based on P-and S-waves generated during a M w 4.2 event that occurred on 18 May 2020.This large event partially destroyed the mine infrastructure in a section about 1,300 m wide (Dineva et al., 2022).The resulting P-and S-waves were recorded at 151 seismic sensors (151 P-wave and 81 S-wave arrival times) as far as about 2,130 m from the preliminary seismic source location.Our preliminary seismic source is the official location (Dineva et al., 2022), even though the damaged area may be as large as 700 × 250 m.For the present study, we specifically revised all P-and S-wave arrival times associated to this event to include as many seismic sensors as possible.It is worth noting that some seismic sensors close to the rupture area were destroyed during the event, limiting the availability of data near the source.The closest sensor to the preliminary seismic source is about 140 m away.

The Reference Solution: McMC Location of a Seismic Event
Our fist experiment obtains a reference solution that consists of estimated source locations obtained using only sensors within a maximum distance from the calibration blast (which is equivalent to have equal weight for all stations within the max distance and ignore stations that are beyond the maximum distance).By considering a range of possible maximum sensor distances, we aim to reproduce the results that would be obtained by different expert opinions.
The observed arrival time data in a vector d are   obs  , where x = P or S, i = 1, …, N, and N = N P + N S = 66, with N P and N S being number of P-and S-wave arrival times, respectively.We apply a standard Markov chain Monte Carlo approach to solve the inverse problem of locating the source.We make four simplifying assumptions: (a) The rock volume is a homogeneous half-space; (b) The seismic event is a point source; (c) The arrival times have associated uncertainties equal to the sampling rate (i.e., σ 0 = 1/6,000 s); and (d) The covariance matrix of the data errors , with I being the identity matrix (see Riva et al., 2023) for details on the methodology).
The model vector m in the inverse problem contains eight parameters.Six parameters are related to the physical model: the coordinates of the seismic source (X s , Y s , Z s ), the origin time measured with respect to the first P-wave arrival time (OT), the P-wave velocity of the half-space (V P ), and the ratio between P-wave velocity and S-wave velocity in the half-space (V P /V S ).Two additional "hyperparameters" (Malinverno & Briggs, 2004) π P and π S , defined below, control the data uncertainties.Thus: In our experiment, all prior probability distributions are uniform within the minimum and maximum values given in Table S1 in Supporting Information S1.
For the given assumptions, the seismic ray paths for model m are straight lines from the seismic source in (X s , Y s , Z s ) to the known position of each seismic sensor.Thus, predicted arrival times   pred  (with x = P, S and i = 1, …, N) can be easily computed from the source-sensor distance, V P , and V P /V S .The vector   =  obs  −  pred  contains the residual differences between predicted and observed arrival times.The hyperparameters π P and π S multiply the error variances that define the error covariance matrix as follows:   2  () =  2 0 ⋅ 10 2  , where x = P or S. Therefore, the error covariance C e (m) is a diagonal matrix that contains the values of   2  () for each P-or S-wave arrival time.In this context, the likelihood function can be written as  where the determinant of the diagonal covariance matrix is the product of its diagonal entries, that is, () .The complete results for one choice of the maximum distance are reported in Figure S1 in Supporting Information S1 for references.

Assigning Data Weights With a Trans-Dimensional Algorithm
The novel approach applied in our second and third experiments (using arrival times from the calibration blast and the natural event) follows the algorithm presented in Piana Agostinetti and Sgattoni (2021), except that here we use spherical shells in space rather than change points in time.In our case, the spherical shells are centered at an approximate preliminary location of the source (see in Supporting Information S1 for a discussion on determining a preliminary source location).The models sampled by McMC are composed by six physical parameters, as done earlier, plus a variable number of parameters related to the spherical shells: the number of shells k, the k-vector of shell radii r k , and the k + 1-vectors of the weights for P-and S-wave arrivals in each shell w k,P and w k,S .Thus:  = (, , , ,  ,  ∕ , , , , , , ) As for the previous reference solution, all Priors for the spherical shell parameters are uniform probability distributions (see Table S2 in Supporting Information S1).Given the model vector m, we compute a modified version of the covariance matrix of the data errors C e (m) that accounts for the data weights as follows: where W(m) is a diagonal matrix that contains the weight assigned to each arrival time; as   *  is diagonal, C e (m) is also diagonal.We remember here, as stated in Section 1, that an equivalent view is that the algorithm samples variances, rather than weights, in each of the spherical shells as in a hierarchical Bayes strategy (Malinverno & Briggs, 2004).In fact, in a Bayesian framework, "data noise" includes observational and theoretical errors.For example, here, as the source-station distance increases, the theoretical "noise," due to a wrong velocity model hypothesis, increases.The resulting variances in the diagonal matrix C e (m) will be consistent with the size of the residual differences between predicted and observed arrival times in the vector   =  obs  −  pred  .Writing δ i as the distance from the source of the sensor recording the i-th arrival time, the entries of W are where w i (m) is where x = P or S depending on whether the i-th arrival time is for a P-or S-wave and  r is the ordered version of the vector r k (i.e.,   1 < . . .<   ).The likelihood function is as in Equation 1 with the covariance matrix C e (m) of Equation 2.
Composing the covariance matrix C e (m) needs a preliminary source location to compute the source-sensor distances δ i .In general, an approximate preliminary location is available soon after the event takes place and can be safely used.If a preliminary event location was not available, the position of the sensor that receives the earliest P-wave arrival can be used as well.Here we use as preliminary locations the actual location of the calibration blast and the location of the natural event published in Dineva et al. (2022); the coordinates of these preliminary locations are in Table S3 in Supporting Information S1.

Results
The reference solution results are shown in Figure 2. Starting with all the available data (all 57 seismic sensors that recorded the blast to a maximum distance of 800 m from the source), we get a posterior mean event location 10.1029/2023GL102983 6 of 10 which is about 12 m away from the blast, with estimated uncertainties as large as 7 m.We then start removing data from sensors farther than 700 m, 600 m, etc., in steps of 100 m (see Figure 2 and "Materials and Methods").The event location uncertainties and the differences with the actual blast position reach a minimum for a maximum sensor distance of 300 m (19 sensors).Considering sensors closer to the source (200 m, 5 sensors) results in an increase in uncertainties and location error.
Our novel approach applied to the blast data gives results that are consistent with those obtained in the reference solution (Figure 3 and Figure S2 in Supporting Information S1).The variation of weights with distance for both P-and S-wave arrival times follows a simple pattern, with a single step decrease at about 380 ± 30 m from the source (Figures 3b and 3c).The weights for S-wave arrival times decrease much more sharply than those for P-waves.This main step is well defined, as seen from the histogram of the sampled shell radii (Figure 3d), although the histogram of the number of shells has a maximum between 5 and 7 (Figure 3a).The sampled weights result in a cloud of event locations that closely reproduces what was found in the reference solution for a maximum distance of 400 m (red vs. black dots in Figure 3e).
In crustal studies, it has been observed that event location uncertainties depend on the azimuthal coverage (Husen et al., 2012).Here we computed the azimuthal coverage of the 3D distribution of seismic sensors (see "Materials and Methods").Azimuthal coverage reaches a nearly stable value at a distance of ca.300 m from the source, and it does not change substantially at greater distances (Figure 3d).The best reference solution was found when selecting stations only within 300 m from the source, which is also close to the distance where the weights obtained in our new method decrease substantially.We apply our data-space exploration algorithm to the arrival times of the natural event (Figure 4).This event has a magnitude M w 4.2, it is composed of several subsequent processes, where the extent of the very first sub-event S1 is likely ca.100-200 m (Dineva et al., 2022)).The final posterior distribution of the source location is close to that initially estimated (Figures 4a and 4b).The pattern of weights with distance is more complex compared to that obtained for the blast.There is a main step at about 1,230 ± 70 m, but also three other maxima in the histogram of shell radii (marked with colored arrows in Figure 4c).The weights for the P-wave arrival times slightly increase from the origin to 150 ± 50 m (gray arrow) and remain near a maximum value between 150 ± 50 and 500 ± 60 m (red dashed arrow).At greater distances, the weight decrease slightly to a nearly constant value out to 1,230 ± 70 m (red arrow), where there is a sharp decrease of almost one order of magnitude.The weights increase again at about 1,900 ± 60 m (blue arrow).We also conducted a test to check whether the overall pattern of weights with distance is significantly affected by the simple parameterization of constant weights in each spherical shell.To this end, we implemented an alternative parameterization where weights are defined at the shell boundaries and vary linearly within each shell (Figure 1c).The pattern of weights with distance obtained with linearly varying weights is very similar to that obtained with constant weights (see yellow contours in Figure 4c, Figures S3 and S4 in Supporting Information S1).The choice of parameterization does not seem to strongly control the variation of weights with distance.

Discussion
In our first test with a controlled blast, the reference solution found with all seismic stations within 300 m from the source could seem to outperform our novel approach, because the mean PPD event location is slightly closer to the blast true position with respect to the mean PPD solution found with the novel approach.However, this difference is small (less than 2 m) and the pattern of the weights found with the novel approach closely mimics a step function (with a step at about 380 m).In case another maximum distance is picked (200, 400, 500 …meters), that is, in case of a different expert opinion, the solution found by our novel approach is better than the reference solution.Thus, we can affirm that our approach outperforms and is more robust than standard approaches in case of limited information about the data.
Comparing the results obtained in the two tests carried out with our novel approach, we note that the pattern of weights with distance seems to be event-dependent and is not a constant in a particular sensor network.While  2b).The yellow square shows the posterior mean source location obtained with the data-weighing method.Colored circles are the posterior mean source locations obtained in the reference solution (same as in Figure 2b).The yellow sun indicates the true position of the calibration blast.
10.1029/2023GL102983 8 of 10 further research would be necessary to determine which event parameters (e.g., magnitude, location) affect the weight pattern, the results indicate that a static workflow for all events would probably introduce artifacts and underestimate the actual uncertainties.In contrast, our approach is adapted to each single event, giving a solution that is statistically consistent and parsimonious (in terms of complexity of the weight pattern parameterization).Alternatively, a more detailed approach to data weighting (e.g., considering a more general 2D/3D spatial distribution of weights, instead of concentric shells, as done in Ghalenoei et al., 2022) could help in reducing the difference in the results obtained for different events.
The relationship between azimuthal coverage of the event and our results is not straightforward.In the controlled blast, the main decrease in the weights we obtains is near the distance where the azimuthal coverage increases substantially (Figures 3d and 4c).On the other hand, there is no clear correspondence between weight patterns and azimuthal coverage in the test with a natural event.This suggests that azimuthal coverage is only one of the factors affecting the reliability of the inverted source location.A workflow based on this parameter (e.g., where distant seismic sensors are removed once the azimuthal gap decreases below a certain threshold) may not give optimal results.In fact, if the gap is larger than 180° with stations in the epicenter near vicinity, a moderately distant station closing this gap may be very useful if the real subsurface velocities are not perfectly well known (which is almost never the case).On the other hand, closing a gap to significantly less than 180° with a single very distant station is at least questionable (if not useless) when considering the uncertainties of phase identification and frequency difference in first arriving/visible wavelets.
The pattern of weights allows us to interpret the results in terms of specific properties of the rock volumes at different distances from the source.We suggest that the seismic sensors closest to the natural event (at distances <150 m, first gray circle in Figure 4a), very likely are in the source area, where the assumption of a point-wise seismic source is not realistic for such a large event.Between the gray and the dashed red circle (distances of 150-500 m) the weights reach their highest values, indicating where the inverse problem assumptions should be valid.Indeed, all sensors within the red dashed circle in Figure 4a are located on the same side of the ore body, where the rock volume is expected to be comparatively homogeneous.Between the dashed and solid red circles in Figure 4a (distances of 500-1,200 m) the weights are still high, but less than in the previous interval.This is likely due to some ray-paths partially crossing the ore body and thus violating the homogeneous rock assumption.Farther than 1,200 m from the source (red circle in Figure 4a), the seismic rays start to densely sample the ore body and the surrounding rocks on both sides of the ore body itself.Here we can expect that the assumption of a homogeneous rock finally breaks down, and the weights decrease significantly.Further investigations are needed to confirm our hypothesis and to check how complex pattern in weights could be related to a less circularity in the data distribution around the seismic source in the case of the natural event than in the case of the blast.
In a more general context, our novel approach can be applied to most of the scientific inference problems, where huge amount of data need to be pre-processed in some way, without introducing bias related to preconceptions of the data-analysts.We mention that our approach only works if data can be ordered or clusterized in some way.Here for example, they are "ordered" with regards to the source-sensor distance.In this case, "ordering" is necessary, but it is not the only way of performing the trans-dimensional data-space exploration.To apply our approach, we need either a metric to be used to "measure" some kind of data-point distance in the data-space, or, equivalently, some kind of data characterization which enables data clustering, where the trans-dimensional approach is used to define the number of data cluster from the data themselves.

Figure 1 .
Figure 1.The arrival time weights w k are associated with a set of k concentric 3D spherical shells with radii r k , centered on a preliminary event location.(a) 2D representation of the spherical shells.(b) Constant weights within each shell.(c) Alternative parameterization with linearly varying weights within each shell.(d)-(e) Seismic data used in this study: sensor locations projected on a horizontal plane (circles) and arrival times for P-waves relative to the earliest recorded arrival time (circle colors).Dashed circles show the distance in meters from the preliminary event location.Important geological features are depicted with colored dots: gray = ore body; blue = clay zones; red/purple = diapir/diabase.Panel (d) shows seismic arrival times for a calibration blast (yellow Sun symbol), and (e) for a M w 4.2 event (yellow star).The red box in (e) shows the smaller area plotted in (d).

Figure 2 .
Figure 2. (a) Blast locations obtained in the reference solution with data for sensors that are within different maximum distances from the actual source location.The last column reports the distance D between the posterior mean and true location of the source.The best location (where D is minimum) is obtained with data from sensors up to 300 m from the blast.The X (South), Y (East), and Z (depth) columns list the posterior mean value for the location coordinates.(b) Source locations sampled by the McMC algorithm for different maximum source-sensor distances projected on the X-Z vertical plane (dots).The maximum source-sensor distances are 800 m (black dots), 700 m (dark blue), 600 m (pink), 500 m (yellow), 400 m (red), 300 m (light blue), and 200 m (green).Colored circles are posterior mean locations.The yellow sun indicates the true position of the calibration blast.

Figure 3 .
Figure 3. Application of the novel data-weighing method to blast recordings.(a) Posterior probability density function (PDF) of the number of spherical shells, approximated by the histogram obtained by McMC sampling.(b) Posterior PDF of the weights assigned to P-wave arrival times as a function of source-sensor distance.Green crosses indicate the distance of each sensor from the source.(c) As in (b) for S-wave arrival times.(d) Posterior PDF of shell radii.The blue dashed line indicates the azimuthal gap as a function of distance from the source (see "Materials and Methods" for a definition).The red dashed line indicate the prior probability distribution for the shell distance.(e) Sampled source locations projected onto a X-Z vertical plane (black dots) compared to source locations in the reference solution for a maximum source-sensor distance of 400 m (red dots; see Figure2b).The yellow square shows the posterior mean source location obtained with the data-weighing method.Colored circles are the posterior mean source locations obtained in the reference solution (same as in Figure2b).The yellow sun indicates the true position of the calibration blast.

Figure 4 .
Figure 4. Application of the novel data-weighing method to recordings of the M w 4.2 natural event.(a) Seismic network geometry (same as in Figure 1e).Colored circles report the position of the main modes in the histogram of sampled shell radii, indicated with colored arrows in panel (c).The inset (b) plots the sampled source locations projected onto a X-Z vertical plane (blue dots) compared to the preliminary location of the event (yellow star).(c) Posterior probability density function (PDF) of shell radii, approximated by the histogram obtained by McMC sampling.The colored arrows indicate the main modes in the posterior PDF, corresponding to the boundaries of the source area (gray arrow), homogeneous rock volume with all sensors on the same side of the ore body (dashed red), homogeneous rock volume (red), heterogeneous rock volume (blue).(d) Posterior PDF of the weights assigned to P-wave arrival times as a function of source-sensor distance.Green crosses indicate the distance of each sensor from the source.The yellow contours display the posterior PDF of the data weights obtained with the alternative parameterization in Figure 1c (see also Figures in Supporting InformationS1).