Abstract

The analysis of urban seismic signals offers valuable insights into urban environments and society. Yet, accurate detection and localization of seismic sources on a city-wide scale with conventional seismographic network is unavailable due to the prohibitive costs of ultra-dense seismic arrays required for imaging high-frequency anthropogenic sources. Here, we leverage existing fiber-optic networks as a distributed acoustic sensing system to accurately locate urban seismic sources and estimate how their intensity varies over time. By repurposing a 50-kilometer telecommunication fiber into an ultra-dense seismic array, we generate spatiotemporal maps of seismic source power (SSP) across San Jose, California. Our approach overcomes the proximity limitations of urban seismic sensing, enabling accurate localization of remote seismic sources generated by urban activities, such as traffic, construction, and school operations. We also show strong correlations between SSP values and environmental noise levels, as well as various persistent urban features, including land use patterns and demographics.

Introduction

Cities, as epicenters of human activity and major contributors to global emissions, have drawn increased research interests aimed at improving urban policies, enhancing the quality of life, and promoting sustainability¹. Urban science employs quantitative and modeling approaches from various disciplines to achieve these goals^2,3. Central to this field is the availability of large-scale datasets capturing various urban “signals” with precise spatial and temporal resolution^4,5. Despite advances in urban sensing and big data approaches^6,7,8, comprehensive urban data access remains limited due to issues, including data ownership, privacy concerns, and high costs^9,10. Recently, urban seismic signals have been used to monitor urban environments cost-effectively and pervasively. These signals carry valuable information for characterizing urban dynamic features such as environmental conditions, traffic patterns, and cultural and social activities^{11,12,13,14,15}. However, seismic signals related to urban activities, particularly those at frequencies above 1 Hz, are subject to scattering and attenuation over short distances¹⁶. Consequently, accurately estimating the spatiotemporal distribution and intensity of urban seismic sources, a process known as seismic source mapping, requires ultra-dense seismic arrays, which are not feasible to install in urban areas due to high costs and disruptions. Hence, it is unclear how these urban seismic sources are distributed within the city and how they are related to dynamic and persistent urban features.

Here, we explore the use of existing and ubiquitous urban infrastructure – telecommunication optical fibers – to map urban seismic sources. By repurposing existing fiber-optic cables into a sensing network using distributed acoustic sensing (DAS)^17,18, we create an ultra-dense and cost-effective seismic array capable of covering extensive urban areas. Integrating seamlessly with the existing telecommunication infrastructure, our approach provides a scalable solution for continuously recording urban seismic signals. While DAS has been used to detect and locate distant natural seismic sources like earthquakes^19,20 and ocean waves^21,22, its applications in urban science^23,24,25,26, such as traffic and footstep detection, have been limited by proximity sensing - they have been used to monitor ground motions and quasi-static deformations due to activities that occur close to the sensors, such as traffic along the road on which the fiber is laid. Moreover, existing DAS applications for urban sensing have been developed for specific purposes, like traffic estimation or bridge health monitoring^27,28, and are not adaptable for general-purpose urban sensing.

In this study, we integrate seismic interferometry and beamforming algorithms to overcome the proximity limitations of previous urban seismic sensing approaches. Our method models the propagation of seismic surface waves to estimate the spatiotemporal distribution of seismic source power (SSP), defined as the average seismic energy per unit of time within an area. Our approach achieves precise localization and estimation of various seismic sources occurring remotely from optical fibers on a two-dimensional spatial urban scale. We validate our seismic source mapping results against various urban seismic activities with known coordinates and timings, including city-wide traffic movements and human activities at construction sites and schools. Additionally, we find that SSP intensity is influenced by and significantly correlated with several persistent urban features, such as land use patterns, average daily traffic, points of interest density, and demographics.

Results

Scalable seismic source mapping

The continuous seismic monitoring capabilities of DAS using existing telecommunication fiber networks present a unique Large-N array²⁹ to capture seismic signals from various urban seismic sources. The recorded seismic signals are analyzed to estimate spatiotemporal maps of SSP (Fig. 1). This analysis utilizes seismic data generated from passive sources and collected through existing, unused fiber-optic cables, known as dark fibers. Therefore, it is highly scalable and integrates seamlessly with the existing telecommunication infrastructure. By connecting an interrogator unit to one end of a 50-km long dark fiber in San Jose, California, we created 50,000 virtual dynamic-strain sensors (DAS channels) continuously sampled at 200 Hz and distributed in one-meter spacing along the fiber, mainly along the roadways. These virtual sensors, each costing the equivalent of a few dollars when splitting the cost of the interrogator unit, operate without interrupting telecommunication services or necessitating on-site sensor installation and maintenance^17,18. We continuously recorded seismic signals from the fiber-optic cables encircling three key regions over six days from September 20th to 25th, 2023. The dataset is detailed in Methods section.

Fig. 1: Seismic source mapping using existing telecommunication fiber-optic cables as seismic sensing arrays using distributed acoustic sensing (DAS).

Methods

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:

$$p(\omega ;{{\bf{x}}}_{r})=a(\omega ;{{\bf{x}}}_{r},{{\bf{x}}}_{s})s(\omega ),$$  (1)

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:

$$a(\omega ;{{\bf{x}}}_{r},{{\bf{x}}}_{s})=\frac{1}{\sqrt{{r}_{s,r}}}\exp \left(-\frac{\omega {r}_{s,r}}{2Q{v}_{\omega }}\right)\exp \left(-\frac{i\omega {r}_{s,r}}{{v}_{\omega }}\right),$$  (2)

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:

$$S(\omega,{{\bf{x}}}_{t})={{\bf{h}}}^{H}(\omega ;{{\bf{x}}}_{t}){\mathbb{E}}\left[{\bf{p}}(\omega ){{\bf{p}}}^{H}(\omega )\right]{\bf{h}}(\omega ;{{\bf{x}}}_{t}),$$  (3)

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:

$${\bf{h}}(\omega ;{{\bf{x}}}_{t})=\frac{1}{\sqrt{N}}\frac{{\bf{a}}(\omega,{{\bf{x}}}_{t})}{\sqrt{{{\bf{a}}}^{H}(\omega ;{{\bf{x}}}_{t}){\bf{a}}(\omega ;{{\bf{x}}}_{t})}}.$$  (4)

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

$$C\left({{\bf{x}}}_{s},{{\bf{x}}}_{r},\tau \right)=\left\{\begin{array}{l}\mathop{\int}\nolimits_{{t}_{s}+\epsilon }^{{t}_{s}+\epsilon+\Delta t}u\left(t+\tau ;{{\bf{x}}}_{r}\right)\cdot u\left(t;{{\bf{x}}}_{s}\right)dt,\quad {{\bf{x}}}_{r} < {{\bf{x}}}_{s}\hfill\\ \mathop{\int}\nolimits_{{t}_{r}+\epsilon }^{{t}_{r}+\epsilon+\Delta t}u\left(t-\tau ;{{\bf{x}}}_{r}\right)\cdot u\left(t;{{\bf{x}}}_{s}\right)dt,\quad {{\bf{x}}}_{r}\ge {{\bf{x}}}_{s},\end{array}\right.$$  (5)

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:

$$\ln \left(\frac{A(\omega ;{{\bf{x}}}_{i})}{A(\omega ;{{\bf{x}}}_{j})}\right)=-\frac{\omega {r}_{i,j}}{2Q{v}_{\omega }}+C,$$ (6)

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:

$$\hat{Q},\hat{C}=\arg \mathop{\min }\limits_{C,Q}\sum _{\omega }\sum _{r}{\left[\ln \left(\frac{{A}_{0}(\omega )}{{A}_{r}(\omega )}\right)+\frac{\omega r}{2Q{v}_{\omega }}-C\right]}^{2},$$ (7)

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.