Passive Seismic Imaging of Urban Environments Using Distributed Acoustic Sensing: A Case Study from Melbourne, Australia Open Access

High‐resolution near‐surface characterization is essential for assessing seismic hazards (e.g., $VS30$⁠), land subsidence, and groundwater changes. Achieving this level of detail requires sensing capabilities with small interstation spacing to capture short‐wavelength structures. Distributed acoustic sensing (DAS) offers a cost effective and scalable solution for high‐resolution imaging and long‐term monitoring through ultradense seismic arrays (Zhan, 2020). DAS operates by converting Rayleigh back‐scattered light —resulting from intrinsic impurities in the fiber—into measurements of longitudinal strain or strain rates. Existing telecommunication fiber‐optic cables can be repurposed for DAS, turning them into thousands of seismic sensors that cover distances spanning tens of kilometers. This use of dark fiber provides a nonintrusive and economical method for imaging short‐wavelength structures (on the order of meters) across otherwise inaccessible areas or terrains and enables continuous, real‐time monitoring at high frequencies (hundreds of Hertz) over extended periods (months to years).

Recent studies have successfully leveraged dark fiber for high‐resolution near‐surface imaging, ranging from select sites (Dou et al., 2017; Spica et al., 2020) to regional‐scale tomographic model (Yang et al., 2022), which is crucial for improving seismic site characterization, such as obtaining $VS30$ measurements. There are challenges in using dark fibers for passive imaging, and various groups have explored different strategies to improve the process. In seismic noise source selection, Ajo‐Franklin et al. (2019) focused on selecting energetic train signals to enhance cross correlation, and Shragge et al. (2021) targeted storm‐induced microseismic energy to capture low‐frequency (<1.0 Hz) Rayleigh waves for imaging deeper depths. In processing steps, Li et al. (2023) and Smolinski et al. (2024) utilized a frequency–wavenumber filter to eliminate spurious wave directions, and Ehsaninezhad et al. (2024) applied denoising techniques to increase investigation depth. Other groups such as Czarny et al. (2023) explored the effects of weather and Mirzanejad et al. (2024) explored the effects of fiber installation on the phase velocity measurements.

In this study, we present a comprehensive workflow designed to address the challenges of performing passive imaging with DAS in urban environments using dark fiber. These challenges include strong diurnal noise source distribution, nonlinear array geometry, and poor fiber coupling to the ground. Using the workflow, we successfully characterized the top 100 m near‐surface shear‐wave velocity in Melbourne, Australia, in a section with no prior near‐surface model or borehole measurements. Our workflow can be easily adapted to other dark fiber systems with complex array geometry, aiming to maximize the use of all available channels, and requires minimal prior knowledge, including seismic source distribution, cable conditions, and ground coupling for velocity measurements, and initial velocity models for inversion.

The DAS experiment repurposes dark fiber from the Australian Light Infrastructure Research Testbed (ALIRT) in Melbourne, Australia. The 76.6‐km‐long fiber is dedicated to research purposes, connecting two photonics laboratories at Royal Melbourne Institute of Technology (RMIT) and Monash University. A 25 km segment of the ALIRT fiber is transformed into a DAS array with 6252 channels at 4 m spacing. The array traces busy arterial roads and tram lines, starting from the Melbourne Central Business District, where the RMIT City Campus is located, and extending through local roads toward the eastern residential suburbs (Fig. 1a), aligning with a grid‐street plan that includes orthogonal sections and the 105‐m‐long Victoria Bridge over the Yarra River. The physical locations of the channels were determined by tap tests at right‐angle turns and known maintenance pits, whereas the remaining channel locations were interpolated, with an estimated maximum error of 20 m (about five channels).

In the experiment, strain‐rate measurements were obtained using a Silixa iDASv2 interrogator with a fixed 10 m gauge length. The DAS time series, originally recorded at 1000 Hz were down‐sampled to 250 Hz for spectra calculation and to 62.5 Hz for passive imaging to reduce file size for efficient computation. The experiment ran intermittently from December 2021 to October 2023, collecting 178 days of data, some of which were contributed to the global DAS month campaign in February 2023 (Wuestefeld et al., 2024).

The DAS data are strongly influenced by diurnal traffic signals from motor vehicles, buses, and trams (Fig. 1b). The ambient seismic amplitude at 1–60 Hz shows a high power spectral density (PSD) during the day, decreasing by about 8 dB at night, with the quietest period between 1 a.m. and 5 a.m. when public transportation services are suspended (Fig. 1c,d). This PSD change is most pronounced in the initial channels near Melbourne’s Central Business District, where many tram routes are located. There is no significant difference in the PSD pattern between weekdays and weekends in Melbourne, as trams, one of the main seismic noise source, run all night on weekends.

Subsurface imaging using noise interferometry has been widely utilized for large‐scale long‐period seismic data over the last two decades (Campillo et al., 2014). DAS data, with its small station spacing allow the recovery of higher frequency surface waves with shorter wavelengths, making it sensitive to shallow subsurface features. The imaging workflow is similar to many conventional imaging approaches, containing three main steps: (1) cross‐correlation; (2) surface‐wave dispersion analysis; and (3) 1D depth profile inversion. However, there are specific considerations addressing issues unique to passive DAS imaging in urban settings.

The calculation of cross‐correlation functions follows the conventional approach in Bensen et al. (2007) with minor modifications. The 15 min time series are prefiltered between 0.1 and 30 Hz, detrended, and have common‐mode noise minimized by removing the median signal amplitude. Spectral whitening is applied to the data segments, and the cross‐correlation functions (CCFs) are calculated in the frequency domain using a graphic processing units‐based Python tool (Yang, 2025). These CCFs are then transformed back to the time domain and stacked over a two‐week period from 9 to 23 February 2022. A frequency–wavenumber filter targeting phase velocities below 50 m/s and above 5000 m/s is applied to the stacked CCFs to suppress spurious arrivals, including waves traveling in the opposite direction, and to minimize the effects of slow‐moving seismic sources from traffic.

Source distribution is known to affect the characteristics of seismic ambient noise, with surface‐wave cross correlation mainly contributed by sources in the stationary phase (Boschi and Weemstra, 2015). In our case, the peak seismic source is from traffic, which is ideally aligned with the fiber, but requires select time windows to maximize its power. We observe that the diurnal seismic noise pattern significantly impacts the CCFs. Cross correlation using daytime recordings recovers dispersive waves, with amplitudes that decay quickly over short distances. In contrast, nighttime recordings reveal nondispersive waves traveling at a speed of 960 m/s over long distances (more than 1 km) (Fig. 2a). This body wave is likely from man‐made structures, such as metallic water pipes buried near the fiber, which have an acoustic wavespeed of 900–1200 m/s. Therefore, we find that the cross correlations should only be computed when traffic signals are most prominent and stable; in this case, between 07:00 and 20:00 local time (see the supplemental material, available to this article).

Using DAS in noise interferometry is further complicated by channels with poor ground coupling, which exhibit inconsistent amplitudes. These issues often stem from extra loops in maintenance pits or aerial sections like the across the Yarra River on the Victoria Bridge, and must be accounted for in cross‐correlation analysis. The geolocalization approach by Biondi et al. (2023) exploits changes in amplitude in earthquake recordings to identify subpar channels. Here, we directly use the CCFs by normalizing the sum of energy per channel and identifying subpar channels with very low amplitude that exceed a conservative threshold of several standard deviations. We successfully linked these channels to known fiber joints and manholes and identified additional sections not connected to known joints (Fig. 2b). The subpar channels are removed from the dispersion curve analysis.

As the DAS sensors are evenly spaced in both time and space, we can use an array‐beamforming approach to calculate dispersion curves, similar to the approach by Jiang and Denolle (2022). We determine the phase velocity for each frequency band between 1 and 25 Hz by binning the DAS channels into subarrays and analyzing the maximum stacked energy of the shifted CCFs (i.e., slant‐stacks).

Several studies have highlighted the limitations of DAS as a single‐component sensor, primarily sensitive to the axial component, making sensor orientation crucial for accurately recovering Rayleigh or Love waves through cross correlation. Cross‐correlating channels with virtual sources aligned in the same direction mainly recover Rayleigh waves, whereas for oblique or orthogonal channels (up to 90°), a mix of Rayleigh and Love waves are recovered, complicating their separation within the dispersion curves (Song et al., 2021; Fang et al., 2023).

In practice, we can minimize the impact of DAS directivity by selecting subarrays and virtual sources that share the same orientation. For the ALIRT DAS array, we ensure each subarray is geometrically linear by checking the slope between each end. The size of the subarray influences both the largest resolvable wavelength and the sharpness of the dispersion curves. Larger subarrays improve wavelength resolution but can result in more rejected dispersion curves due to nonlinear geometry. We chose a subarray size of 50 channels, which effectively covers wavelengths of 3–25 Hz. To maximize dispersion curve calculations without compromising quality, we performed slant stacks using ±200 channels adjacent to the subarray as virtual sources for each of the subarray channels, and stack the dispersion curves together. Accounting for both subarray geometry and individual channel quality, we achieved a channel utilization of 47.6%, resulting in 2845 dispersion curves spanning a length of 11.38 km (Fig. 3).

Prior to inverting for 1D shear‐wave velocity profiles along the array, the slant‐stacks from the beamforming analysis are averaged at every five channels to improve robustness. The phase velocities of the fundamental mode of Rayleigh waves were manually picked to ensure continuous dispersion curves within 5–25 Hz for the inversion. Some scattering in the picking remains and are considered as part of the noise in the Markov chain Monte Carlo (McMC) fitting (Fig. 4). We focused on the fundamental mode of Rayleigh waves in the 5–25 Hz range, as they are more stable and less prone to contamination from higher modes at higher frequencies. The phase velocity is sensitive to shear‐wave velocity down to 300 m, with higher frequencies corresponding to shallower depths, assuming a velocity of 900 m/s at 1 Hz.

To translate dispersion curves into shear velocities as a function of depth, we use a Bayesian inversion approach based on McMC sampling (Magrini et al., 2024). Bayesian sampling allows a robust exploration of the parameter space (e.g., Sambridge and Mosegaard, 2002) and is therefore suited to addressing the surface‐wave‐velocity inverse problem across a range of scales, from near‐surface (Pasquet and Bodet, 2017) to crustal‐upper mantle depths (Bodin, Sambridge, Tkalcic, et al., 2012; Magrini et al., 2023). Similar Bayesian inversion approach has also been applied to studies using DAS data in Spica et al. (2020) and Smolinski et al. (2024). The global search inherent to this approach is especially effective for regions without prior velocity models, as in our case.

We assume the subsurface is composed of discrete layers overlying a homogeneous half‐space, with each layer being isotropic and uniform in $VS$⁠, $VP$⁠, and density. We further assume a $VP$ to $VS$ ratio at 1.89, which is typical for heavily jointed rocks with low velocities (Barton, 2006), and use the Nafe‐Drake curve (Brocher, 2005) to convert $VP$ to density; this allows us to accurately model density in the presence of low $VP$⁠.

We discretize the subsurface using transdimensional Voronoi tessellation, choosing a uniform prior for the number of layers between 5 and 20, and for the depth of the Voronoi nuclei ranging from 0 to 200 m, which determine the layer interfaces. For $VS$ (Fig. 4), we utilize a uniform prior with depth‐dependent bounds, defined by a deviation of ±40% from the model of Chandler et al. (2005), which was obtained using the spatial autocorrelation method at a site northwest of Melbourne. The Chandler model, and therefore the $VS$ prior, follows a power‐law shape consistent with the Dix‐type relation for surface waves (Haney and Tsai, 2015). Compared to a simple uniform prior, using a discretization‐informed prior—where the range of admissible velocity values varies with depth—improves McMC convergence by narrowing the parameter space and excluding unrealistic models, such as those with higher velocities at the surface compared to 200 m depth. Finally, we adopt a hierarchical approach by assigning uniform priors to both the data noise standard deviation (0–180 m/s) and the noise correlation (0–1), with the correlation function given by a double exponential decay (e.g., Bodin, Sambridge, Rawlinson, and Arroucau, 2012).

We ran 50 Markov chains for 500,000 iterations each, with every chain starting from a different random realization. We saved one model every 200 iterations from the 100,000th iteration onward, resulting in 100,000 posterior samples. The similarity between mean and median velocities indicates a stable inversion process with few outliers. Model uncertainty is shown by the width of the 10th–90th percentile range from posterior models (Fig. 4). The higher quality of the input data (i.e., dispersion curves with defined beams) results in narrower percentile window width. Velocities at depths shallower than 40 m are better constrained due to more measurements and higher frequency sensitivity, whereas deeper velocities are less constrained. The final depth profiles presented are limited to 100 m.

We find that the seismic velocities derived from the inversion compare well with the variations in the geologically mapped lithology of Melbourne. With a span of just 25 km, the ALIRT DAS array covers three major lithological groups, as shown in Figure 5, reflecting the rich Melbourne geological history (Neilson, 2018; Geological Survey of Victoria, 2022). During the Silurian and Lower Devonian epochs, significant deposition of marine siltstone and sandstone formed the Melbourne Formation, potentially extending to 10 km deep. In the Tertiary period, deposition shifted to the Red Bluff Sandstone formation, consisting of clay and sandstone. Outcrops of this formation, ∼15–30 m thick, are found about 20 km south of the DAS array. In the last two million years (late Pliocene to early Holocene), basalt flows from the Newer Volcanic Group covered much of western Victoria, concluding at Melbourne and currently bounded by the Yarra River. The basalt formation is characterized by multiple thin flows, generally estimated to be 0.5–10 m thick, with thicker accumulations observed in ancient valleys, such as at a site 15 km west of the array, where flows reach up to 60 m. After the volcanism in the Quaternary, the study region experienced fluvial erosion, shaping the unconsolidated colluvium on hillsides and alluvium along the river and its tributaries.

The combined 1D profiles in Figure 5 reveal clear variations in shear‐wave velocities and depth, which can be understood in relation to the lithology, marked by color blocks at the top of each profile. Starting from the west of the profile, at sites on Newer Volcanic Group basalt flows (e.g., 6.3 km [Fig. 4], profiles show low velocities of 700–1000 m/s near the surface, increasing to around 2000 m/s at depths of 40–50 m. The velocity profiles at distances of 7.1 and 21.8 km (Fig. 4), associated with the Silurian age Melbourne Formation (Fig. 5), display a flat profile with velocities increasing from 1300 m/s at the surface to 2000 m/s at 100 m. Between distances 10 and 21 km, at sites on Red Bluff Sandstone formation, the velocity profiles show a slow layer of 700 m/s in the top 20–30 m, with velocities increasing to 1800–2300 m/s at 100 m depth.

Sites located on the Newer Volcanic Group and Red Bluff Sandstone formations show a distinct transition to higher velocities at depths of 20–40 m. This transition likely indicates a lithological change from the shallower formations to the older and potentially more consolidated Melbourne Formation siltstone. Near the surface, the shear‐wave speed of basalt, subject to weathering processes, can be as low as 500 m/s, as observed in borehole measurements (Yaede et al., 2015). The Red Bluff Sandstone, being poorly consolidated, exhibits lower velocities compared to the underlying siltstone formation. In contrast, sites on the Melbourne Formation display a less distinct depth transition. Given the consistent lithology, we hypothesize that the reduction in velocity at shallower depths may be due to increased fracturing near the surface, whereas compaction at greater depths leads to an increase in velocity, as seen in velocity measurements along a 65 m borehole (Holbrook et al., 2019).

Sections with Quaternary alluvium and colluvium, such as at 9.2 km exhibit poorer quality dispersion curves. Velocities start at 1200 m/s at the surface, transitioning to 2000 m/s at a depth of 40 m, indicating that the colluvium extends to 40 m and overlies older Melbourne Formation with velocities of ∼2000 m/s, similar to those observed in other profiles. In addition, within the Red Bluff Sandstone formation, a distinct change in velocity profiles between distances of 11 and 12.3 km indicates lithological variations that are not reflected in the published lithology map (Geological Survey of Victoria, 2022). The velocity profiles at these sites are similar to those with colluvium mapped at the surface between the Melbourne Formation (around the 7.7 and 9.1 km distance marks).

A key advantage of passive seismic imaging with dense seismic arrays, in particular DAS is the ability to recover subsurface information that complements and surpasses the depth limitations of other noninvasive geotechnical surveying methods, such as ground penetrating radar (less than 20 m). High‐resolution passive imaging is especially valuable in regions like in Melbourne, where borehole data are absent and the existing lithology map relies on limited outcrops or is completely obscured in urban areas, and is based on interpretations from local geomorphological features. The $VS$ profile obtained from the DAS array not only aligns well with surface lithology but also provides estimates of the sediment layer thickness before transitioning to unweathered bedrock, offering insights into how overlying layers affect bedrock weathering.

One crucial goal in passive imaging is achieving consistent resolution across a wide range of depths, from near surface through the crust, to understand processes like fluid pathways, weathering, and erosion informed by deeper geology (Holbrook et al., 2019). Although the reliance on high‐frequency traffic noise (1–40 Hz) as the primary seismic noise source limits the overall depth resolution, this frequency range focuses on the top 30–100 m and can resolve important parameters such as $VS30$ for seismic hazard characterization. Stacking longer time durations (over two weeks) does not significantly enhance the signal‐to‐noise ratio for lower frequencies (below 5 Hz). The resolvable wavelength also depends on the subarray size, requiring consideration of DAS directivity and array layout, both geometry and channel quality, which are difficult to modify when using dark fiber. Advancements, such as the integration of amorphous arrays‐combining seismometers with DAS (Nayak et al., 2023) and the joint inversion of both Rayleigh and Love waves (Song et al., 2021) hold promise for improving not only the depth resolution but also enhancing spatial resolution, enabling the transition from 1D to 3D high‐resolution velocity models of the near‐surface.

Achieving high‐quality CCFs in urban environments can be challenging. We find that maximizing stationary seismic sources, such as daytime traffic, improves the noise cross‐correlation results. This is important as filtering out acoustic signals from underground pipes using f‐k methods is difficult due to their similar phase velocities and frequencies to those observed in the subsurface. Cross correlations can be contaminated by cross‐line sources, including broadside arrivals, but this can be minimized with f‐k filtering of waves traveling in opposite directions. Despite these difficulties, the stability of CCFs with a short time‐window stacking shows promising potential for long‐term monitoring of subsurface seismic velocity changes, as demonstrated by (Shen et al., 2024).

We demonstrate an effective workflow of strategically selecting traffic noise and channels to improve surface‐wave dispersion curve measurements, directly addressing the challenges of passive imaging using dark fiber in urban environments including strong diurnal noise patterns, array geometry, and poor fiber coupling to the ground. This workflow not only enables the recovery of a high‐resolution near‐surface model for Melbourne, Australia, addressing a crucial gap in the region’s geophysical surveys but also provides a valuable approach for enhancing near‐surface characterization in other seismically active urban areas.

The decimated data from the Australian Light Infrastructure Research Testbed (ALIRT) distributed acoustic sensing (DAS) deployment, along with metadata following Lai et al. (2024), can be accessed at the National Computational Infrastructure (NCI) through the 2030 Geophysics Project (doi: 10.25914/kf96-nb26), with support from NCI, AuScope, and Australian Research Data Commons (ARDC). This study uses the open‐source package BayesBay at https://github.com/fmagrini/bayes-bay (last accessed August 2024) for inversion. The supplemental material contains a concise overview of the passive imaging workflow, along with an h5 file of the top 100 m of the shear‐wave velocity model presented in Figure 5.

The authors acknowledge that there are no conflicts of interest recorded.

The authors sincerely thank the Editor‐in‐Chief Keith Koper, Associate Editor Jonathan Ajo‐Franklin, and two anonymous reviewers for their helpful feedback and suggestions, which improved this article. The team gratefully acknowledges Royal Melbourne Institute of Technology (RMIT) InPac Lab (A. Mitchell, S. Tan, L. Broadley, A. Boes, and B. Corcoran) and Australian Academic and Research Network (AARNet) (T. Rayner) for access to Australian Light Infrastructure Research Testbed (ALIRT) dark fiber; ANU‐RSES technical team (J. Bryne) for computational support; and the National Computational Infrastructure data group (J. Croucher, H. Hollman, and N. Rees) for open data sharing assistance. The authors also thank M. Lambert and W. Zeng at U. Adelaide for the valuable discussions. This research was supported by the Australian Research Council Future Fellowship “Lighting Up Fiber for Seismic Imaging” (project number FT210100440), Discovery Project “Measuring the seismic pulse of the Earth using fiber optics” (project number DP230100277), and Centre of Excellence in Optical Microcombs for Breakthrough Science (project number CE230100006), which are funded by the Australian Government. The distributed acoustic sensing (DAS) experiment was enabled by AuScope via the National Collaborative Research Infrastructure Strategy (NCRIS).