Seismic Source Tracking With Six Degree‐of‐Freedom Ground Motion Observations

Back azimuth (BAz) information can be determined from combined measurements of rotations and translations at a single site. Such six degree‐of‐freedom (6‐DoF) measurements are reasonably stable in delivering similar information compared to a small‐scale array of three‐component seismometers. Here we investigate whether a 6‐DoF approach is applicable to tracking seismic sources. While common approaches determining the timing and location of energy sources generating seismic waves rely on the information of P‐waves, here we use S‐waves. We track back azimuths of directly arriving SH‐waves in the two‐dimensional case, P‐converted SV‐waves, direct SV‐ and direct SH‐waves in the three‐dimensional case. For data analysis, we compare a cross‐correlation approach using a grid‐search optimization algorithm with a polarization analysis method. We successfully recover the rupture path and rupture velocity with only one station, under the assumption of an approximately known fault location. Using more than one station, rupture imaging in space and time is possible without a priori assumptions. We discuss the effects of rupture directivity, supershear rupture velocity, source‐receiver geometry, wavefield interference, and noise. We verify our approach with the analysis of moving traffic noise sources using 6‐DoF observations. The collocated classic seismometer and newly built ring laser gyroscope ROMY near Munich, Germany, allow us to record high‐fidelity, broadband 6‐DoF (particle velocity and rotational rate) ground motions. We successfully track vehicles and estimate their speed while traveling along a nearby highway using the estimated BAz as a function of time of a single station observation.

which shows high robustness with respect to measurement uncertainties. We verify the concept in synthetic 2-D experiments analyzing SH-wave polarity and discuss the applicability and robustness of the developed methodologies. Third, we demonstrate earthquake rupture tracking in 3-D media from the rotation polarization caused by P-converted SV-waves and direct SV-and SH-waves. We analyze the effect of interfering arrivals and non-uniform slip rate distribution. We discuss source-receiver scales and geometry as well as challenges of the method for future global applications. Finally, we show a real data application of the proposed approach, which mimics the rupture imaging process. We track moving traffic seismic sources using a high-resolution 6-DoF point measurement.

Seismic Source Tracking
Most common techniques to image earthquake properties using array data can be divided into two categories which are both based on analyzing the phase information of P-waves. In contrast to finite slip inversions, no detailed knowledge of Green's functions and source properties is necessary. Methods of the first category are based on conventional array measurements. Termed back-projection methods, seismic energy radiation is imaged by applying array beam-forming techniques. Back-projection was for the first time successfully demonstrated for the 2004 Sumatra-Andaman Earthquake (Ishii et al., 2005;Krüger & Ohrnberger, 2005). Directivity effects were utilized to characterize faulting mechanisms (Ammon et al., 2005).
Methods of the second category track earthquake rupture by estimating the back azimuth (BAz) of incoming waves with a single-station. In polarization analysis, the three translational components of standard seismometers can be used to estimate the BAz and incidence angle of incoming waves (Flinn, 1965;Greenhalgh et al., 2005;Montalbetti & Kanasewich, 1970;Vidale, 1986). Bayer et al. (2012) developed a single-station approach to track moving sources by polarization analysis of local and regional P-wave arrivals. They normalize the BAz variation with respect to a known hypocenter. Frohlich and Pulliam (1999) pointed out that, compared to traveltime-based methods, classic single-station approaches suffer from several ambiguities, for example, 180° BAz errors. The joint analysis of translational and rotational motions can help overcome such drawbacks.

6-DoF Ground Motions
The complete wavefield excited by an infinitesimally small deformation can be described by the three components of translation, three components of rotation, and six components of strain (Aki & Richards, 2002). However, until recently, seismology is dominated by translational observations (vertical, N-S, E-W), sometimes combined with strain measurements. Translational motion is the movement of a particle along an axis. In contrast, rotational motion describes the particle movement around an axis. Information on rotations has been widely ignored, mainly, because of measurement difficulties. 6-DoF information is obtained from measuring in addition to three translational components also three components of rotational motion. This increase in information compared to classical observations has the potential to improve existing methods and creates new opportunities for research and industry (e.g., Igel et al., 2015;Schmelzbach et al., 2018). Rotational motions can be derived from arrays of conventional single or multi-component sensors (e.g., Huang, 2003;Spudich et al., 1995;Spudich & Fletcher, 2009;Suryanto et al., 2006). However, these methods are limited by array spacing as well as local heterogeneities and site effects. Classic translational measurements are also sensitive to tilt, i.e., the horizontal components of the rotation vector (Graizer & Kalkan, 2008;van Driel et al., 2015). Recent advances in fiber-optic gyroscopes and ring laser-based sensors show that applicable, single-station measurements for translation and rotation are within reach (Bernauer et al., , 2012Schreiber & Wells, 2013).
The earthquake source process and the interaction of the wavefield with a free surface or heterogeneities of the Earth can excite rotational ground motions. In isotropic media the rotational motion     ( , , ) T x y z ω can be described by a linear combination of spatial derivatives of the translational particle displacement motion  ( , , ) T x y z u u u u (e.g., Cochard et al., 2006): where × denotes the cross product and ∂ k denotes spatial derivatives with respect to x k . The same relation is valid for the time derivatives, the rotation rate  ω and particle velocity v. At the free surface, the stress-free boundary condition is slightly modified (e.g., Schmelzbach et al., 2018). While inside isotropic media the curl operator separates the S-wave field, in anisotropic media even (quasi-)P-waves can have a rotational component (Pham et al., 2010). Local phase velocities and the BAz can be estimated from 6-DoF measurements at a single-station due to the relation between translations and rotations (Edme & Yuan, 2016;Hadziioannou et al., 2012;Igel et al., 2007;Pancha et al., 2000;Sollberger et al., 2018).

BAz Estimation From a Single Station
We first test two different methods to track sources of seismic energy (earthquake rupture) in simple 2-D examples. Both methods are based on a plane wave assumption and analyze the polarity of directly arriving SH-waves at a single station. Since the region of energy radiation moves during the rupture across the fault plane, we utilize sliding windows moving throughout the signal to determine the evolution of the signal source direction. In each time window, the BAz is estimated and a temporal trend can be derived by comparing all windows. When this information is combined with a priori knowledge on the fault or with data from other stations, the rupture propagation and its velocity can be estimated.
The CC (cross-correlation) method is a grid-search optimization algorithm that relies on the interdependence of transverse translational motion and vertical rotation. The CC is a measure for the similarity between two signals and the CC coefficient provides a measure for the degree of similarity (see Figure 1). A CC coefficient of 1 implies perfect similarity, a value of −1 means anti-correlation. Similar to the approach by Igel et al. (2007), we estimate the BAz by rotating the horizontal acceleration components in small steps around all possible BAz (0-360°) and cross-correlating successively with the vertical rotation rate. A zero-lag normalized CC coefficient is used. For a noise-free signal, the CC coefficient in the grid-search is a function without a clear maximum. It is a step-function that jumps from −1 to 1. Therefore, we use the two zero transitions of the step-function instead of the global maximum. We expect that the central position between the zero transitions corresponds to the actual BAz.
The second method was introduced by Sollberger et al. (2018) and we refer to it hereinafter as polarization analysis. In comparison to the CC method, it is more flexible and can be applied to P-, SV-, SH-, Rayleigh-, and Love-waves. Instead of a MUSIC likelihood function (Schmidt, 1986), we use one based on the classical power spectrum, because we expect a more stable result. We assume that the global maximum of the likelihood function is related to the actual BAz. It is necessary to define a parameter space for evaluating the likelihood function in a grid search. While the CC method requires the definition of the BAz increments, the increments for the S-wave velocity and the incident angle must be additionally defined for the polarization analysis. Both methods are illustrated for a plane wave in Figure 1.
The difficulty of retrieving BAz (source directivity) from 3-D observations with solely translational motions is due to two challenges: 1) The 180° ambiguity in BAz estimates if only translational motions are recorded (Langston & Liang, 2008).
Considering that rotational motions are essentially the curl of translational motions, the polarity of rotation will reverse in case of an opposite propagating direction while translation polarity remains unchanged. Thus the joint analysis of rotation and translation will help to remove the 180° ambiguity when locating the sources 2) Translation records suffer from interfering with different types of wavefields at the free surface, i.e., P-and SV/SH-waves, Rayleigh-and Love-waves are generally intermixed in recorded horizontal components. However, rotational motions naturally separate P-and S-waves as P-waves do not generate rotational motions in isotropic media. SV-and SH-waves (the same as Rayleigh-and Love-waves) are also naturally separated despite unknown source locations since SV-waves or Rayleigh-waves only generate rotational motions on horizontal components while SH related (Love-wave related) rotational motions being isolated on vertical components. We can therefore take advantage of the fact that two horizontal rotational components contain exclusively SV-or Rayleigh-waves Based on the plane wave assumption, the ratio between the two horizontal rotational components is directly related to the BAz according to: where   n and   e denote the north-south and east-west components of rotational rate, respectively. This simple relationship is specially useful for estimating source directivity and it is independent of any possible radiation pattern that the source might have. The negative θ BAz value derived from the inverse tangent function is converted to the value within 0 and 180° by adding 180°. To remove the 180° ambiguity, we then compare the rotated transverse component of rotational rate (  t ) based on θ BAz to the vertical component of acceleration (a z ). If they are positively correlated, 180° should be added to θ BAz .

Combining Multiple Stations
It is possible to track the horizontal propagation of a rupture with only one station, assuming that the BAz changes are correctly determined from the seismogram and the fault position is known a priori. An infinitesimal thin ray could be constructed for each estimated BAz in the direction of the directly arriving waves. The intersections of these rays with the fault position would show the temporal evolution of the rupture. In case of an unknown fault, at least two stations are necessary for the tracking process. But for more than two stations the rays will not intersect in exactly one point, since measurement errors and inaccuracies in the methodology can't be excluded completely. We want to take these uncertainties in the BAz into account by using wider beams instead of infinitesimal thin rays. We define the shape function p(x, y, t) of each beam as  where the standard deviation σ i (t) is defined individually for each station i, Φ i (x, y, t) ∈ [0, 180°] denotes the angular distance from an arbitrary point in space to the estimated BAz of a specific station i and N stations denotes the total number of stations. Note that the shape function is time-dependent and is defined for each horizontal position (x, y). The time framework is defined in such a way that t = 0 corresponds to the first arrival at a respective station. For a specific time-step t 0 the amplitudes of all beams are added up and we assume that the most likely source position is close to the maximum value of p(x, y, t 0 ). In Figure 2, we show examples of shape functions for different BAz errors and a 2-D representation of p(x, y, t 0 ) for two stations using Equation 3 with N stations = 2.

Synthetic Case Studies in 2-D
In the following, we use data of elastic wave simulations in 2-D to demonstrate possible applications on a fundamental level and their limitations. 2-D wave propagation simulations are performed using the spectral element package se2wave. The mesh representation and support for MPI parallelism in se2wave is provided via PETSc (Balay et al., 2019(Balay et al., , 1997. We note that while homogeneous 2-D media SH rotation rate ought to be proportional to ground acceleration, the subsequently presented synthetic examples gradually increase in complexity and benefit from consistent and reproducible numerical calculations.
We describe two different test cases. The fault position is known in the first case and we try to track the spatial and temporal evolution with only a single station. In the second case, we assume that the fault position is unknown and the individual results of many stations are combined.

Rupture Tracking With a Single 6-DoF Station
We model a pure strike-slip earthquake embedded in a 2-D homogeneous medium. The unilateral rupture has a constant speed of 80% of the S-wave velocity v s and it is implemented as a line of double-couple point sources. We choose the source time function of each point source to be an ordinary Gaussian. Additionally, we slightly randomize the onset time and the seismic moment to render our synthetic study more realistic. Each station records the horizontal accelerations a x , a y and the vertical rotation rate   z of the directly incoming P-and S-waves (middle panels of Figure 3). These seismograms demonstrate that P-waves do not have a YUAN ET AL.
10.1029/2020JB021112 5 of 22 Figure 2. The BAz uncertainties can be described by a shape function. Left: each curve shows the shape function for a specific BAz error σ. The beam angle describes the broadness around the estimated BAz which corresponds to 0°. Right: Equation 3 is illustrated in 2-D at a certain time-step for two stations (white triangles). An expected error of σ = 1.5° describes the broadness of each beam. Its highest value (yellow area) is expected to be close to the source position. BAz, back azimuth.
amplitude -shape function rotational component in an isotropic and homogeneous medium. Station A only records weak P-wave amplitudes due to the perpendicular position with respect to the rupture. Due to different BAz between source and station, the duration of the SH-arrivals is different for stations A and B. The maximum expected BAz variation for station A is about 11.5 and 5° for station B. We estimate the BAz changes by moving a sliding window of 1.5 s length through the SH-wave signal of the seismograms.
In each window the BAz is estimated by the polarization analysis and the CC method and the results are illustrated in the bottom subplots of Figure 3. Each point in these graphs represents the central position of a time window. The points are color-coded relative to the first and last SH-wave arrival. Both methods provide the same linear trend and the results are nearly perfectly overlapping.  Using SH-waves for rupture tracking in a single station approach. Top: rupture (red line) and receiver positions (blue triangles) are shown in horizontal dimensions x, y. The color-coded rays indicate the estimated BAz variations at each station. The rupture is correctly tracked. Middle: recorded seismograms of horizontal accelerations a x , a y and vertical rotation rate   z . Bottom: the BAz is estimated by two different methods from the direct SH-arrivals. The estimates are divided into three sub-windows in each of which rupture speed is determined. The true rupture velocity is 80% v s . BAz, back azimuth; CC, cross-correlation.

(a)
If the rupture or a certain part of the rupture propagates at a constant rupture speed and the starting and ending point are approximately known, it is possible to estimate the rupture speed by trigonometric considerations. The rupture speed depends on the S-wave velocity, the rupture length, the rupture duration measured at the receiver, and the orientation of rupture direction and receiver. Both bottom plots in Figure 3 are divided in three sub-windows. We note that the signal duration is not necessarily the same in each subplot. We determine the velocity in each sub-window by fitting a straight line through the estimated BAz and express it relative to the known S-wave velocity. In this manner, we estimate rupture velocities close to the real value of 80% v s in all sub-windows.

Direct Estimation of Rupture Velocity
We evaluate the sensitivity of the proposed 6-DoF tracking methods to variations in earthquake rupture propagation speed across the fault. Reliable, far-field estimation of rupture velocity is important to constrain earthquake dynamics, stress drop, and implications for seismic hazard but is inherently difficult because of the intermixing of rupture geometry and rise time in controlling the P-and S-wave pulse shapes (e.g., McGuire & Kaneko, 2018).
We test three different rupture scenarios (Figure 4). The model setup is the same as in the previous 2-D tests, but here only the first half of the rupture has a constant speed of 80% v s . The second half breaks with a constant velocity of 40%, 60% or 150% of v s as shown in Figures 4a-4c, respectively. Earthquake ruptures can propagate at sub-Rayleigh or at intersonic speeds (e.g., Archuleta, 1984;Dunham et al., 2003;Gabriel et al., 2012) and a speed of 150% v s means that the rupture is propagating faster than the radiated SH-waves and close to the P-wave speed. This effect is referred to as supershear rupture speeds.
In Figure 4 we estimate the BAz changes for each case of velocity variation at station A in the same way as in Figure 3. Each column of Figure 4 represents the result for a specific velocity jump. The BAz results are divided into three sub-windows, in which the rupture velocity is calculated, respectively. The final ve-YUAN ET AL.   locity results, which we express relative to the S-wave velocity, are illustrated for each sub-window by a red horizontal line in the lower panels. We indicate the true rupture speed by blue dashed lines. In each test scenario, there is a significant increase or decrease visible from the starting velocity of 80% v s in the first sub-window to the final rupture velocity of 40%, 60%, and 150% in the last sub-windows. While the speed of the second half of the rupture is determined nearly perfectly, the starting velocity is slightly underestimated.

Rupture Tracking in Heterogeneous Media
We expect that 6-DoF rupture tracking is more difficult in heterogeneous materials since reflected and scattered energy will contaminate the directly arriving SH-waves. We rerun the simulation of Figure 3 now perturbing the homogeneous model by adding a normally distributed random material heterogeneity. We add a variation of up to ± 5% to density, P-wave and S-wave velocity in the medium. In the numerical simulations quadrilateral elements are employed, each possessing piece-wise constant material properties and edge lengths ∼ 100 m. Material properties in each element are perturbed independently of each other and no smoothing of the piece-wise constant properties is applied between neighboring elements.
The seismograms recorded at station A are shown in Figure 5. Due to reflections in the material, P-and S-waves are no longer perfectly separated. Reflected phases are visible in all components after the dominant SH-arrival. The lower panel shows the tracking result, in which the true starting and ending BAz for station A are marked by blue dashed lines. The spatial dimensions of the rupture are nearly perfectly estimated.
Although we expect a nearly straight line for the temporal evolution, there are higher deflections than in the homogeneous model. However, on average a rupture speed of 77% v s is determined, which is very close to the true velocity. The BAz deflections are increasing for stations in the higher distance and for stations that are placed in a geometrical orientation, in which the SH-wave amplitudes are less dominant compared to the P-wave amplitudes.

Rupture Tracking for Unknown Simple and Complex Rupture Paths and Directivity Effects
Rupture tracking with only a single station is possible if the fault or more explicitly the rupture path is known a priori. In the following, we assume that the fault position is unknown. Since a single station is not enough to track the rupture in this case, we here combine the BAz estimates of many stations. The BAz is still calculated in a single-station approach at each receiver, but the final tracking results of all stations are combined.
As described in Section 2.2, for a certain time-step we send a virtual beam back from each station in the direction of the rupture. By using a broad beam instead of a thin ray, we here take BAz uncertainties into account. We add up the amplitudes of all beams and the maximum is expected to be the most likely source point. The station coordinates and the BAz changes at each receiver are the input parameters for Equation 3.
However, there is another issue that is referred to as a consistent time-frame. Due to different BAz between stations and rupture, the SH-arrivals at each station have a varying length of time (compare to the seismograms of Figure 3). The signal duration of the target wavefield is mainly affected by source-receiver geometry, rupture velocity, and 3-D subsurface velocities. The time-shift of the sliding window has to be corrected at each station for this effect. Otherwise, an offset of the estimated rupture position from the actual location is expected, even if the BAz is correctly determined. Previous studies neglect such directivity effects expecting only small deviations.
First, we verify that a 6-DoF method provides accurate results for simple and more complex fault geometry then we discuss the importance of directivity effects. We apply a time correction to the BAz estimates by assuming that the start and the end of the directly arriving SH-waves are visible in the signal. This is done in the rotational component of ground motions due to its high sensitivity to shear motions.
In the following, we track a simple unilateral rupture at five stations. The stations are placed in an asymmetrical pattern around the rupture with different distances to the source. The stations are situated in such a way that the resolution is about the same for both spatial dimensions. The medium and rupture parameters are equal to the homogeneous model in the previous section. The SH-waves at each station are picked manually and the BAz change is independently determined of the other stations. In Figure 6, rupture tracking results are shown for three different time-steps (a, b, and c). The estimated starting position is shown in Figure 6a and the ending position in Figure 6c. Animation S1 visualizes the continuous rupture imaging (Movie S1: ms01.avi).
In the first subplot of Figure 6, we show the arrangement of rupture (red line) and stations (white triangles).
The We repeat the same experiment as presented in Figure 7 for a more complex rupture geometry. An animation of the continuous rupture imaging for the complex rupture geometry is provided in the supporting information (Movie S2: ms02.avi). The rupture propagates on three horizontally displaced segments of different lengths. The four subfigures show the rupture tracking at different time steps. Even in this more complex situation, the rupture is correctly tracked and the fault offsets are visible in the final tracking results.
The length of time during which body waves arrive directly varies for different station locations in dependence on the rupture position. Such directivity effects will cause artifacts in the tracking result if the information of many stations is combined. In both previous experiments, we correct for directivity effects by picking the start and end time of the SH-arrivals in the seismograms. We expect that, in real data, it is YUAN ET AL.
10.1029/2020JB021112 9 of 22 Thus, an important question arises: How is rupture tracking affected if only the first arrival is visible in the seismograms? Bayer et al. (2012) neglected directivity effects in a comparable approach with classic 3C data for P-waves and used the same time-shift for the BAz estimation window. While conventional back-projection does not include a time correction for directivity effects (e.g., Ishii et al., 2005), P-wave-based rupture tracking utilizing beam-forming has been shown to require correction for the varying locations of the seismic sources (e.g., Krüger & Ohrnberger, 2005). Figure 8 illustrates the impact of directivity effects in an additional numerical experiment. The rupture is tracked by a small array of four stations which has a relatively small opening angle. The impact of directivity effects is expected to be significantly smaller for an array with a small opening angle. However, if the angle is too small, it is not possible to determine both spatial dimensions of the rupture in good quality. In Figure 8, the station positions are shown in the first map. In the left subplot, we show the final result for BAz estimation in which we corrected for directivity effects. The resolution of the x coordinate is not as good as in the previous results, since a smaller array is used, but the rupture is still tracked correctly. The right subplot shows the result for the same data, but this time the time-shift of the sliding window in the BAz estimation is assumed the same for all stations. The starting position is still correctly determined but later estimated points show a systematic deviation from the true rupture path. Even if the rupture area and its linear trend are roughly tracked, the geometry is not correctly derived. By neglecting directivity effects it is possible to track the beginning of the rupture, but not its complete spatial evolution.
YUAN ET AL.

Rupture Tracking in 3-D Heterogeneous Media
We extend the presented 2-D findings by examining the stability and accuracy of 6-DoF rupture tracking in 3-D heterogeneous media where multi-phases interfere with each other.
The opening angle for a given earthquake, i.e., the detected BAz variation, depends on epicentral distance and station azimuth. The resolving power at a single station will potentially decrease with the increasing epicentral distance while increase with increasing inclination angles. In the following, we first estimate the expected opening angle for 3-D rupture tracking of earthquakes as a function of rupture length, epicentral distances and inclination angles, while avoiding the polarization uncertainty dropping below the noise level. Then we perform 3-D synthetic tests for rupture tracking with 6-DoF measurements.

3-D Single-Station Opening Angles
We define the opening angle of a specific station as the difference between the BAz for the starting position of the unilateral rupture and the BAz for the ending position (Bayer et al., 2012). A large opening angle is desirable to minimize uncertainties which corresponds to a short distance between receivers and earthquakes. But the distance has to be large enough to fulfill the plane wave assumption while keeping in mind that the analysis of the polarization is less efficient for signals in which many different phases are interfering.
The following description is a purely geometrical concept to demonstrate the expected scaling of opening angles. If we define the rupture path as a straight line on a sphere, we can describe the geometry between the receiver and rupture by a large triangle. In Figure 9, we illustrate the opening angle α for fault lengths between 100 and 1,000 km (  windows). The triangle is not necessarily isosceles, which is expressed by the inclination angle δ ∈ [0°, 90°]. The maximum opening angle occurs for δ = 90°, i.e., the station is situated perpendicular to the center of the rupture and the triangle is isosceles. In general, α increases for larger faults as well as for shorter station distances. By applying the spherical law of cosines, the opening angle α can be described by the side lengths s 1 , s 2 , and l, where l denotes the rupture length: In an epicentral distance of 30° the opening angle for a fault of 1,000 km is about 18° (see the upper right window of Figure 9). For example, the mainly unilateral 2004, Great Sumatra-Andaman earthquake ruptured across a length of about 1,200 km. If we increase the inclination angle from 0° to 30° at the same distance, the opening angle decreases to about 9°. This may be still sufficient for an estimation of the spatial and temporal evolution of the rupture with a single station. Significantly shorter ruptures or very small inclination angles, however, will lead to an opening angle of only a few degrees, which is challenging to track.
YUAN ET AL.

Synthetic Case Studies in 3-D
To verify the stability and accuracy of the proposed rupture tracking approach using 6-DoF measurements in 3-D, we calculate synthetic seismograms using Instaseis, an efficient tool for generating synthetic global seismograms using Green's function databases generated in 2.5D axisymmetric spectral element simulations using AxiSEM  and utilizing 1-D axisymmetric velocity models. We here assume the 1-D isotropic Preliminary Reference Earth Model model with attenuation effects and the highest frequency of up to 0.5 Hz to generate the synthetic data. Although Instaseis does not allow the direct analysis of rotational components, we derive them using a densely spaced array, i.e., gradient-based array-derived rotation.
In the following synthetic tests, we place four additional stations surrounding the central station with a spatial interval of 100 m (see the upper left subplot in Figure 10). The array-derived rotation is calculated based on a finite-difference scheme (Langston, 2007;Spudich et al., 1995), as the rotational motions will be simplified to horizontal spatial gradients of translational motions at the free surface where vertical stress equals zero (Robertsson & Curtis, 2002). Instaseis enables us to handle finite ruptures represented by an arbitrary number of point sources. The simulated rupture consists of six subevents and propagates approximately from south-east to north-west (indicated by the black arrow in Figure 10).
All subevents are assigned a uniform faulting mechanism (strike: 336°, rake: 114°, dip: 7°) and are evenly distributed along the fault plane at the same depth (10 km). The total rupture length is about 236 km. Considering that the rupture speed and radiated energy can be largely affected by local structural properties and stress conditions, we slightly randomize the source time functions of each subevent in terms of slip rate and initiation time (see the source time functions in Figures 11b, 11d, 11f, and 11h. White noise has been added YUAN ET AL.

10.1029/2020JB021112
13 of 22 to each synthetic data set, which ends up with a signal-to-noise ratio (SNR) of 25. The SNR value is defined as SNR = 10log 10 (A signal /A noise ) using the logarithmic decibel scale, where A denotes the root mean square amplitude in a certain time window (150 s after the corresponding arrival).
The estimated BAz at stations ST1-4 with the epicentral distances of 25°, 14°, 50°, and 45°, respectively, is expected to continuously increase during the rupture tracking process ( Figure 10). For stations ST1-3, we use the singular value decomposition (SVD) algorithm for robust ratio calculations of Equation 2 (Greenhalgh et al., 2018), with a sliding time window of 25 s. The recorded rotational motions in the two horizontal rotational components are mostly resulting from P-converted SV-waves at the Earth surface at ST1-2 and direct SV-waves at ST3. For station ST4, we apply the CC method to the vertical rotational component and the two horizontal translational components, in order to focus on direct SH-waves, with the same sliding window as the one for ST1-3. We select stations ST3-4 with a larger epicentral distance such that the direct SV-and SH-waves can be separated in time from surface waves. We generate one data set for each station. The BAz estimate for each data set as a function of time is shown in Figures 11a, 11c, 11e, and 11g (dashed black lines denote the theoretical starting and ending BAz of the rupture at the given station). The solid dots represent the theoretical BAz of each subevent. The arrows indicate the onset of direct P-or S-arrivals. The corresponding source time functions are plotted in Figures 11b, 11d, 11f, and 11h.
We show that the estimated BAz during the rupture tracking process is generally accurate and consistent at all four stations ( Figure 11) using either P-converted SV-waves (at ST1-2), direct SV-(at ST3), or direct SHwaves (at ST3). However, the gradient of the estimated BAz variation is not constant.
This can be attributed to two factors: 1) the uneven onset time and slip rate across the finite sources. Theoretically, the slope variation of the estimated BAz is supposed to directly indicate the changes of rupture speed as we have discussed in Section 3.1. However, since we randomize the source time functions of all subevents, the slope of the estimated BAz should not be strictly invariant. In Figure 11d, we notice that there is a stronger subevent between 60 s and 80 s, which may lead to the plateau between 260 s and 280 s in Figure 11c. smaller BAz) and later stronger arrival (with a bigger BAz), which ends up with a bias toward the bigger BAz when performing SVD analysis. 2) the interference of non-direct arrivals generated by different subevents. The earlier free-surface related multiples and the reflections from subsurface discontinuities could both interfere with the later direct arrivals. This is a major issue we may have to deal with in real data analysis.
The intensity of the interference may vary with different factors, such as radiation patterns of earthquakes, epicentral distances, and 3-D velocity structures.
Combining the BAz estimates from multiple 6-DoF stations (ST1-4) and various S-waves (P-converted SV-, direct SV-and SH-waves), we apply the same shape function (Equation 3) to the 3-D example. As shown in Figure 12, the maximum amplitude of the multi-beams focuses on the theoretical positions of each subevent within the corresponding time window.
The influence of varying levels of noise and of the length of the sliding window is evaluated for station ST3 (Figure 10) based on the polarization analysis of direct SV-waves in the same rupture scenario introduced earlier in this section. As shown in Figures 13a and 13b, the error of the estimated BAz increases with the decreasing SNR as expected. When comparing the time-variant BAz and the theoretical BAz of each subevent (solid dots in Figures 13b and 13c), we find that it is quite stable and accurate in cases with an SNR above 20. The overall trend of the BAz variation is still preserved despite the BAz errors from signals with relatively low SNR (≤10).
The length of sliding windows is another factor that may affect the stability and accuracy of the proposed BAz estimation. In Figure 13c, we fix the SNR at 35 and apply three distinct window lengths (10, 20, and 30 s) for analysis with the same highest frequency in the synthetic data (0.5 Hz). A shorter sliding window theoretically provides higher resolution in identifying subevents, but in practice, we tend to apply longer windows to stabilize BAz estimation and mitigate the effect of interfering arrivals and/or noisy data as is shown in Figure 13c.

Real Case Studies -Tracking Moving Traffic Noise Sources Using 6-DoF Observation
Due to the lack of real 6-DoF observations of large earthquakes at suitable epicentral distances, we verify the proposed methodology by tracking moving traffic seismic noise sources, which might be a useful analogy to the real rupture process regardless of the scale difference. The collocated classic triaxial seismometer (STS-2) at station FUR belonging to the (alias?) and newly built ring laser gyroscope ROMY Igel et al., 2020) at the Geophysical Observatory Fürstenfeldbruck near Munich, Germany, allow us to record high-fidelity, broadband 6-DoF (particle velocity and rotational rate) ground motions. In Figure 14a, we show the site map of the 6-DoF station. As the station is not far from the highway (blue line on the map), traffic noise is expected to be dominant in the 6-DoF observation for certain frequencies. To avoid overlapping issues and better identify each passing vehicle, we specially choose a data example around midnight when there is no heavy traffic. After converting the particle velocity record to particle acceleration, the 6-DoF data is then detrended, band-pass filtered to 1-20 Hz and downsampled to 50 Hz.
In the top three panels of Figure 14b, we show a nearly one-hour continuous record of the two horizontal acceleration components and the vertical rotational rate, from which we can see consistent traffic-induced signals. Considering the distance and the relative position of the highway and the station, the BAz variations of inbound vehicles (from southeast to northwest) should roughly decrease from 100° to 0° and then keep decreasing from 360° to 300°. The outbound ones are just the reverse. As is shown in the bottom panel of Figure 14b We only apply the CC method here focusing on the SH/Love waves as the horizontal rotational components of the ring laser gyroscope are still under improvement. From the zoom-in plot of one selected time window (Figure 14c), we can calculate the moving speed of the vehicle through the ratio of distance and time over a certain BAz variation range. The estimated speed when BAz changing from 0° to 70° is ∼90 km/h, which is empirically a reasonable estimate at this time of night.

Discussion
In this study, we explore the potential of using 6-DoF observations to track seismic sources by exploiting 1) the correlation of translational and rotational motion observations of SH-waves and 2) the polarization filtering effect of pure horizontal rotational motions. We demonstrate with synthetic and real data that tracking seismic sources is possible provided that 6-DoF measurements are taken at an appropriate epicentral distance and direction to the source. In synthetic tests, we show that direct estimates of rupture velocity can be derived under sub-Rayleigh and supershear speed variations along the fault. We also show that -as long as the rupture-induced shear waves can be identified, the directivity effect can (and should) be corrected for. In case the rupture path is not known a priori, multiple 6-DoF stations with good azimuth coverage are required to track the rupture process.
The presented synthetic models are all unidirectional propagating. We do expect considerably more complexity if 1) the rupture is bilateral, 2) the finite source is extremely complex, or 3) in presence of complex 3-D velocity structures. For any of the complex rupture scenarios listed above, however, more advanced processing techniques need to be developed in future studies. Despite their simplicity, the presented synthetic examples pave the way for successful observational tracking of moving traffic sources using the highly sensitive ring laser gyroscope ROMY.
The advantage of the proposed method lies in the fast and relatively easily obtained seismic source tracking using point stations and taking S-waves fully into account; potentially complementing P-wave based methods and fully complex data-driven kinematic or dynamic earthquake source inversion (Ide, 2007;Gallovič et al., 2019, and references therein).
Structural discontinuities, e.g. posed by radially stratified media such as included in the 3-D synthetic examples in Section 4.2, as well as 3-D heterogeneity, may introduce interfering P-wave and S-wave arrivals from various epicentral distances. Choosing an optimal sliding window length has to balance trade-offs in resolution and interfering arrivals and noise. In our analysis time windows with a single wave-type arrival are not required. For example, the 3-D synthetic and moving traffic real data examples exhibit direct SV and Rayleigh waves mixed with SH and Love waves, respectively. In difference to conventional array-based techniques, 6-DoF measurements can mostly mitigate the P-coda interference due to the nature of rotational motions allowing to image seismic source processes using direct S-waves. Severe interference of multiple types of waves may hinder the here presented application of rupture tracking methods using S-waves. How-YUAN ET AL. (b) (c) (a) ever, this issue might be mitigated in the case of 6-DoF measurements thanks to the inherent wavefield separation in the rotational components, i.e., only SV-waves or Rayleigh waves are presented in horizontal rotation and only SH-waves or Loves wave are presented in vertical rotation. As is shown in Figures 11e-11h, we are able to capture the rupture process when applying both the polarization and the CC method onto direct SV-waves and direct SH-waves, respectively. This illustrates potential as a useful complement to classical back-projection earthquake rupture imaging which solely relies on P-wave information.
Through the traffic-induced seismic noise real data example, we have further verified the feasibility and effectiveness of tracking various seismic sources with 6-DoF point measurement. The results show that moving sources can be precisely located. The vehicle speed estimates are also reasonable assuming a known highway location. The real data example resembles the rupture imaging process assuming finite earthquake sources and may be readily applied to future 6-DoF observations of earthquake rupture. We note that the relatively short source-receiver distance (of roughly 200 m) in this example is presumably not affecting the proposed analysis technique. Unlike conventional travel-time based methods which are limited in accurately measuring time differences in case of limited source-receiver distance, here the target wave type is always in-phase in translational and rotational components and the back azimuth estimation is relying on their relative amplitude difference.
Compared to the seismic array method, 6-DoF point measurements are superior for alleviating site effects. We note that the complexity in both earthquake rupture and wave propagation within 3-D Earth may challenge BAz estimates. With the advent of the first commercial broadband portable rotational seismometer systems Yuan et al., 2020), direct observations of 6-DoF motions of large earthquakes become feasible. The sensitivity of the instrument allows for recording large earthquakes with high signalto-noise ratios. The presented methods with respect to rotational seismology are applicable to (combined) strain observations. With the increasing accuracy of distributed acoustic sensing (DAS) type measurements application of this method to DAS observations should be further explored (Jousset et al., 2018;Lindsey et al., 2017;Yu et al., 2019).

Conclusions
6-DoF single-station observations allow the extraction of wavefield information comparable to small-scale seismic arrays (e.g., Igel et al., 2015;Schmelzbach et al., 2018;Sollberger et al., 2018). In particular, estimates of phase velocities and sub-receiver physical velocities, incidence and BAz angles are possible. We show that such 6-DoF observations allow in principle to track the location of sources of seismic energy and discuss sensitivity and challenges to methods based on cross-correlation or polarization analysis, respectively. Investigating the potential of emerging 6-DoF observations in the context of earthquake physics, the here developed approaches can be generalized to arbitrary sources of seismic energy such as environmental sources, volcanic sources and atmospheric sources as well as to DAS type strain measurements. While we demonstrate the potential of the proposed point measurement for fast and easy seismic source tracking, future efforts may focus on accounting for complexities such as bilateral earthquake rupture, geometric source complexity and complex 3-D velocity structures.