Urban sensing using existing fiber-optic networks

Data description and pre-processing

This study uses DAS data collected over six days, from September 20th to 25th, 2023, in San Jose, California. Data acquisition was conducted using a QuantX interrogator provided by Luna – OptaSense⁵², which functioned at a sampling frequency of 200 Hz and was configured with a gauge length of 10 meters. The interrogated telecommunication fiber-optic cable spans approximately 50 kilometers in length. This setup facilitated 50,000 distributed dynamic-strain sensors along the whole cable, with each sensor spaced 1 meter apart. For seismic source mapping, we identified three specific regions, as shown in Fig. 1a. These regions were selected based on the layout of the telecommunication cables, which form a loop within each area, thereby providing a geometric constraint for beamforming. Specifically, the lengths of the telecommunication cables in regions 1, 2, and 3 measure 9.5 km, 4 km, and 5.6 km, respectively, which all fall within the operational sensing range of the DAS system.

The preprocessing of DAS signals involves two primary steps: bandpass filtering and excluding outliers. After signal detrending, a bandpass filter (1 to 20 Hz) is applied. This band is selected to include urban seismic signals from sources of interest above 1 Hz while filtering out low-frequency quasi-static signals, high-frequency noise, and near-field signals. Subsequently, the data is segmented into 10-minute windows for analyzing the temporal variations in SSP. To mitigate the effects of large amplitude outliers, such as those caused by direct impacts on the fiber, signals in each window that exceed the 99th percentile in amplitude are replaced with the median value obtained from a spatial window of 50 DAS channels. This window is centered on the channel being replaced at the corresponding time. Furthermore, to reduce computational costs while preserving sufficient spatial resolution, DAS signals are subsampled by a factor of 10 in the spatial domain.

Seismic source power estimation with DAS

We employ the beamforming technique^31,53 to estimate the power of urban seismic sources at different locations by shifting and stacking the DAS recorded seismic traces according to a wave propagation model. The beamforming approach aims to estimate the coherent wave energy traversing the DAS array and to characterize its propagation attributes⁵⁴. Our analysis focuses on analyzing the propagation of surface waves to locate urban seismic sources. The beamforming is performed in the frequency domain because of the dispersion characteristics of surface waves, where the wave velocity varies with frequency⁵⁵. Our beamforming implementation is based on the software package Acoular⁵⁶.

Initially, we consider the analysis of a single seismic source located at x_s using a DAS array with N channels. The seismic signal recorded by the r-th DAS channel at position x_r and frequency ω is:

where s(ω) denotes the seismic amplitude characterized at the source location at frequency ω, and a(ω, x_r, x_s) is the transfer function that accounts for signal attenuation and phase delays at frequency ω. Assuming a monopole source with a frequency-independent attenuation, the transfer function is defined as:

where r_s,r = ∣x_s − x_r∣ denotes the Euclidean distance between the source and the DAS channel at x_r, v_ω is the wave propagation velocity at frequency ω, and Q is the attenuation quality factor. Phase distortion due to the attenuation is not considered here. The vector of seismic signals at the DAS array at frequency ω due to a source at x_s is given by p(ω) = a(ω; x_s)s(ω), where vector a(ω; x_s) includes all individual transfer functions and compensates for the time delays and attenuation of the seismic wave traveling from the source to the N DAS channels.

The estimated seismic power at a position x_t and frequency ω is estimated using the real-value auto-power spectrum:

where the superscript H denotes the Hermitian transpose, and \({\mathbb{E}}[\cdot ]\) is the expectation operator. The cross-spectral matrix \({\mathbb{E}}\left[{\bf{p}}(\omega ){{\bf{p}}}^{H}(\omega )\right]\) has dimensions N × N, with each element at the i-th row and j-th column corresponding to the cross-spectrum density function of the signals between the i-th and j-th DAS channels at frequency ω. The cross-spectrum function is the Fourier transform of the cross-covariance function of the time-domain signals between these channels. The steering vector h(ω, x_t) is to calculate the weighted sum of the DAS signals using complex-valued weight factors at the location x_t and frequency ω. Applying this steering vector maximizes output power when the assumed and actual source positions match: S(x_t = x_s) > S(x_t ≠ x_s). The calculated output power also should approximate the source power, i.e., \(S(\omega ;{{\bf{x}}}_{t}={{\bf{x}}}_{s})\approx {\mathbb{E}}[s(\omega ){s}^{H}(\omega )]\). Various steering vector formulations exist; we adopt the formulation as suggested in refs. ^57,58:

In our analysis, we make two additional assumptions. First, when performing beamforming, the directionality effect of DAS is not considered. DAS records the combined horizontal components of the wavefield, influenced by incident angle and cable orientation. Characterizing directional sensitivity is challenging due to unknown source polarization. Moreover, seismic signals from multiple sources are considered to be uncorrelated, allowing S(ω; x_t) to sum their contributions. This algorithm remains valid even in scenarios involving multiple sources.

Prior to beamforming, it is necessary to estimate two key parameters for each studied region: wave propagation velocity and the attenuation quality factor.

Surface wave velocity estimation

To estimate surface wave velocities, we construct Virtual Shot Gathers (VSGs) from vehicle-induced surface waves recorded by DAS. As vehicles travel along a roadside DAS array, they produce two types of signals: the quasi-static deformation due to the vehicles’ weight and the vehicle-induced surface waves, predominantly Rayleigh waves. To utilize these surface waves, we first separate the quasi-static deformation (below 1 Hz) and surface waves (1-20 Hz) using low-pass and band-pass filtering of the DAS data, respectively. The quasi-static deformation signals are used to track the locations of moving vehicles on the DAS record through a specialized Kalman filter²³. This algorithm calculates the arrival times of vehicles at each DAS channel using a prominence-based peak detection method and recursively determines the posterior probability of vehicle arrival times across space (in the direction of vehicle motion) by integrating spatial-dependent vehicle detection results from multiple DAS channels. Vehicle tracking results are obtained by converting the estimated arrival times into vehicle locations.

The tracked vehicle trajectories on the DAS records further enable us to select corresponding surface-wave windows (Supplementary Fig. 1a and b). To avoid interference from nearby vehicle-induced wavefields, we select isolated vehicles with a minimum separation of 25 seconds between other vehicles. We construct VSGs by performing cross-correlations of the pivot trace with other traces within these surface-wave windows. Vehicles passing by generate both forward- and backward-propagating waves on either side of the vehicle’s trajectory (Supplementary Fig 1c). We distinctly handle the positive and negative offset sections of wave interferometry¹⁴. For the backward-propagating waves, i.e., waves traveling in the direction opposite to the vehicle movement, we define the cross-correlation function \(C\left({{\bf{x}}}_{s},{{\bf{x}}}_{r},\tau \right)\) as

where \(u\left(t;{\bf{x}}\right)\) is the recorded DAS trace at time t and location x, \(C\left({{\bf{x}}}_{s},{{\bf{x}}}_{r},\tau \right)\) denotes the inter-channel correlation between recorded DAS strain wavefields at two channel pairs x_s and x_r with τ denoting the time lag. Here, \(C\left({{\bf{x}}}_{s},{{\bf{x}}}_{r},\tau \right)\) approximates the wavefield as if we place a source at x_s and a receiver at x_r. For the negative-offset section in VSGs, channels traversed by the vehicle (x_r < x_s) are cross-correlated with the pivot trace u(t, x_s) within the time window [t_s + ϵ, t_s + ϵ + Δt]. Here, t_s denotes the vehicle’s arrival time at the virtual source location x_s, Δt denotes the selected time window length for cross-correlation, and ϵ is a time lag introduced to avoid near-field effects. For the positive-offset section of VSGs, the cross-correlation is performed in the time window of [t_r + ϵ, t_r + ϵ + Δt], where t_r represents the arrival time of the traveling vehicle at virtual receiver location x_r. In the case of forward-propagating waves, we employ the time window [t_r − ϵ − Δt, t_r − ϵ] for the negative-offset section (x_r < x_s), and [t_s − ϵ − Δt, t_s − ϵ] for the positive-offset segment (x_r ≥ x_s). Then, we stack VSGs from many individual vehicles at the same location to form the final VSGs to enhance the signal-to-noise ratio.

The approach we used here to build the VSGs is more computationally efficient than the conventional ambient noise approach, as only windowed data is required for the cross-correlation. Moreover, this method yields more robust VSGs construction without relying on the assumption of uniform source distribution in the ambient noise interferometry⁵⁹.

We perform the aforementioned algorithm to construct VSGs along the fiber cable across all three regions of interest. We obtain VSGs at 2, 5, and 8 locations by stacking 485, 2520, and 3492 vehicles along the fiber in Region 1, 2, and 3, respectively. For all three regions, each VSG utilized 300 DAS channels with channel spacing of 1 m, corresponding to an offset range of -150 m to 150 m. The variation in the number of vehicles isolated for VSGs is due to differences in fiber layouts and traffic volumes across the studied regions. In particular, Region 1 experiences less traffic compared to the other two regions. Supplementary Fig. 2a–c show the stacked VSGs in the three regions, and Supplementary Fig. 2d–f show the phase velocity dispersion images via the phase-shift method⁶⁰. The anti-causal component of the VSGs (negative time lag) typically originates from vehicles in distant opposing traffic lanes¹⁴ and is excluded from the dispersion analysis.

For our beamforming algorithm, we use the fundamental mode dispersion curve in each region, assuming a uniform velocity model within each region. These dispersion curves are estimated by averaging the results across all VSG locations within each region, resulting in a single dispersion curve per region. The fundamental mode dispersion curves for surface waves derived from the VSGs are shown in Supplementary Fig. 2g–i. The study area exhibits geological diversity, as illustrated in the Supplementary Fig. 3. In the computed frequency range (1.5 Hz to 8 Hz), Region 1 exhibits an S-wave velocity structure ranging from around 400 m/s to 650 m/s, whereas Region 2 displays a broader velocity profile, spanning from 320 m/s to 750 m/s. Region 3 has lower velocity structures compared to the other two regions. Based on these dispersion curves, we do not observe significant variations in surface wave velocity across different locations within each region (gray lines), supporting the adequacy of using a single surface wave velocity profile for our current approach. Also, we focus on low frequencies (1.5 Hz to 8 Hz), which propagate through the less variable parts of the subsurface compared to the top few meters.

Attenuation quality factor estimation

The attenuation quality factor, Q, is another essential parameter for computing the transfer function in Eq. (2) before estimating the seismic source power. We assume a constant attenuation quality factor⁶¹ across the subsurface in each region. The logarithm of the spectral ratio between two seismic signals at different locations can be expressed as:

where A(ω; x_i) and A(ω; x_j) represent the energy of seismic waves at frequency ω at locations x_i and x_j, respectively. The constant C is an intercept term, and r_i,j = ∣x_i − x_j∣ denotes the Euclidean distance between the two locations.

In our analysis, seismic signals induced by trucks are utilized to estimate the attenuation factor because trucks are prevalent and strong seismic sources in all three studied regions. Also, the truck-induced seismic signals exhibit a broader peak in the dominant frequency response compared to other urban seismic sources, such as trains or construction activities (Supplementary Fig. 4). The first step in estimating the attenuation factor involves transforming the truck-induced seismic signals to the frequency domain using a Fourier transform for each DAS signal. Based on Eq. (6), Q is then estimated by solving the following least square problem:

where A₀(ω) is the reference power spectral density (PSD) of DAS signals at frequency ω. This reference energy is computed as the average PSD from DAS channels located at a distance of 50 to 90 meters from the truck source to avoid signal clipping near the source that could underestimate the energy. A_r(ω) is the average PSD of signals from DAS channels within a distance of r − d/2 to r + d/2 from the reference signals’ center, with d set at 10 meters. This calculation also aids in assessing the sensing coverage of our DAS system for various urban seismic sources (Supplementary Fig. 5). Supplementary Fig. 6a–c show the attenuation of truck-induced seismic power in the 2.8-3.4 Hz range across the three studied regions. The plots in Supplementary Fig. 6d–f show the logarithmic spectral ratios at varying distances from the truck source, where the slope of the blue dashed line is proportional to Q⁻¹. The estimated attenuation quality factors for regions 1, 2, and 3 are 6, 8, and 14, respectively.

It is important to note that DAS signal amplitude is influenced not only by subsurface structures but also by the ground coupling of the fiber-optic cable and the directional sensitivity of the wavefield. The lack of coupling data across the telecommunication network limits our ability to isolate attenuation solely attributed to subsurface structures. Furthermore, accounting for the directional sensitivity of the wavefield requires detailed information on source polarization⁶², which is often challenging to obtain for passive sources. Therefore, the attenuation determined reflects a composite effect of subsurface characteristics, cable coupling, and directionality. In this context, the derived Q is an approximation for these multifaceted influences.

Seismic source mapping

To construct spatial maps of SSP across the three studied regions (Fig. 4a), we have partitioned each area into a grid layout. Each grid cell within these grids measures 50 meters by 50 meters, a size chosen to balance the granularity of the spatial map. This ensures that each cell adequately represents local variations while maintaining manageable data volumes for processing. The power spectrum of the seismic source in the frequency range of 1.5 to 8 Hz for each grid cell is estimated using Eq. (3). This specific frequency range is selected to effectively capture the primary frequencies associated with urban seismic activities while excluding high-frequency near-field noise. This choice is further supported by observations that many common urban activities display dominant frequencies within this spectrum (Supplementary Fig. 4). To minimize the influence of near-field sources, such as vehicles traversing directly above the fiber, the seismic source power for each grid cell is estimated using DAS channels located more than 100 meters away from that grid cell. Seismic source mapping can be conducted for a specific frequency of interest or across the entire frequency band by aggregating the maps of each individual frequency. Supplementary Fig. 7 visualizes the seismic source intensity (seismic source power per unit of area) of each census block, with the hatch pattern indicating the land use zoning categories.

The accuracy and effectiveness of our seismic source mapping are verified through a series of validations. These tests compare the estimated source locations with known coordinates of various urban activities. The validation scenarios include tracking moving trains (e.g., Supplementary Fig. 8 and Fig. 2a, d), locating trucks (e.g., Fig. 2b, e), identifying a construction site (Fig. 2c, f), and monitoring activities around a school district (Supplementary Fig. 9). We observe high SSP values aligning with the shape of the fiber array, primarily because the fiber routes follow major roads, such as arterial streets and highways, in each region. These roads experience higher traffic volumes, generating more traffic-induced seismic power compared to smaller roads or areas farther from the fiber array. Consequently, the observed high-power regions largely reflect the alignment of seismic source hotspots along these major roads. Furthermore, SSP maps are calculated for each 10-minute segment of the DAS data. These maps illustrate the daily rhythms and spatial variations of urban dynamics (Supplementary Movies 1 and 2).

There are several limitations and future directions of our methodology that should be noted. First, accurate seismic source localization using the DAS system requires coverage from telecommunication optical fibers across sufficient back azimuths for effective beamforming. A better result is achieved with telecommunication optical fibers encircling the studied regions, creating a geometric constraint essential for effective beamforming. For instance, it would be challenging to determine the source locations using only a straight fiber cable. Also, our analysis reveals that the sensing range of DAS data varies with the energy and frequency of different sources (Supplementary Fig. 5). It is difficult to accurately detect weak and high-frequency seismic sources distant from the optical fiber network. Nonetheless, the extensive urban telecommunication network offers the potential for multi-back azimuth coverage and dense sensing infrastructure. Future research could explore the optimal configuration of optical fibers to maximize the accuracy of seismic source mapping and quantify the sensing coverage for different sources.

Furthermore, although our methodology utilizes distant channels to reduce the influence of near-field effects, the attenuation of far-field sources can still bias the spatial distribution of energy toward near-field activity. Future work should aim to address these limitations by developing techniques to differentiate between near-field and far-field contributions more effectively. While regional heterogeneity of subsurface properties, including the surface wave velocity and attenuation, are taken into account, improving the spatial resolution of these properties could enhance the accuracy of seismic source mapping. Future work should also consider the complexities associated with existing dark fibers, including variations in coupling conditions and conduit materials. Addressing these challenges could involve developing efficient calibration methods. For example, controlled driving tests using vehicles with known weights, combined with precise information about the fiber’s location and depth, could provide a practical and scalable solution. Looking ahead, as our urban sensing approach leveraging existing fiber-optic networks gains traction, the expanded spatial coverage and increased density of these networks will provide richer datasets and finer spatial resolution, enabling the development of more accurate and detailed three-dimensional surface wave velocity models. It is also worthwhile to explore alternative beamforming algorithms, which could potentially improve spatial resolution and enhance side-lobe suppression³¹. Given the challenges associated with the unknown characterization of urban sources, we do not include the effects of fiber and source directionality in our model. Future studies, particularly those involving controlled experiments, could provide valuable insights into the impact of source characteristics and DAS directional sensitivities on source mapping results. The recorded seismic signals demonstrate distinct patterns for various types of urban events (Supplementary Fig. 4). This suggests that pattern recognition and machine learning algorithms could be applied to analyze massive amounts of DAS data, enabling efficient and automated monitoring of diverse urban activities.