DAS (distributed acoustic sensing)

Distributed acoustic sensing (DAS) turns a fibre-optic cable into a dense line of seismic sensors, and has become one of the fastest-growing acquisition modalities (downhole, along roads and railways, on dark telecom fibre, on the seafloor). An interrogator unit sends laser pulses down the fibre and measures the phase of the light Rayleigh-backscattered from each section of it. That phase is sensitive to the axial strain of the fibre, so the instrument effectively records strain (or, more commonly, strain rate) along the fibre direction at thousands of positions (“channels”) simultaneously.

Deepwave models DAS data with the deepwave.das module, entirely as receiver-side post-processing of the particle-velocity receivers that the elastic (and acoustic) propagators already provide – no new propagator or kernel is involved, and everything is differentiable, so DAS data can be used directly as an FWI objective.

What DAS measures

For a fibre with unit tangent \(\hat t\) at position \(x\), the axial strain rate is the tangential-tangential component of the strain-rate tensor of the particle-velocity field \(v\):

\[\dot\varepsilon_{tt}(x, t) = \hat t \cdot \bigl(\nabla v(x, t)\bigr) \cdot \hat t = \sum_{ij} t_i \, t_j \, \frac{\partial v_i}{\partial x_j}.\]

A real interrogator does not measure this pointwise quantity. It compares the phase of backscattered light from two sections of fibre a gauge length \(L_g\) apart (typically metres to tens of metres), which is the strain rate averaged over the gauge – exactly a finite difference of the tangential particle velocity between the two gauge ends:

\[\dot\varepsilon(x, t) \approx \frac{v_t(x + \tfrac12 L_g \hat t,\, t) - v_t(x - \tfrac12 L_g \hat t,\, t)}{L_g}, \qquad v_t = \hat t \cdot v.\]

This finite-difference-of-two-points form is both the physical truth (the interferometric gauge) and the natural implementation on top of Deepwave: place pairs of particle-velocity receivers at the two gauge ends, project onto the tangent, difference, and divide by \(L_g\). The gauge ends are generally off-grid; Deepwave’s floating-point receiver coordinates place them exactly (band-limited Hicks interpolation, with the staggered velocity components handled per-component – see Usage). As \(L_g \to 0\) the recipe converges (quadratically, like any centred difference) to the pointwise \(\dot\varepsilon_{tt}\). Keeping \(L_g\) finite instead reproduces what an instrument with that gauge length actually reports, including the loss of sensitivity when the gauge spans an appreciable fraction of a wavelength.

Broadside insensitivity (\(\cos^2\theta\) directivity)

The double projection onto \(\hat t\) gives DAS its characteristic directivity. For a plane P wave travelling in direction \(\hat n\) with particle velocity \(v = A\, \hat n\, g(t - \hat n \cdot x / c)\),

\[\frac{\partial v_i}{\partial x_j} = -\frac{A}{c}\, n_i\, n_j\, g'(\cdot) \quad\Longrightarrow\quad \dot\varepsilon_{tt} = -\frac{A}{c}\, g'(\cdot)\, (\hat t \cdot \hat n)^2 = -\frac{A}{c}\, g'(\cdot) \cos^2\theta,\]

where \(\theta\) is the angle between the fibre and the propagation direction. A fibre is therefore most sensitive to P waves travelling along it and blind to P waves crossing it at right angles (“broadside”), with the amplitude scale set by \((\omega / c)\, |v|\). This falls out of the implementation automatically and is verified numerically in tests/test_das.py: P waves crossing a straight fibre at \(0^\circ/45^\circ/90^\circ\) give measured amplitude ratios of \(1 : 0.507 : 0.019\) against the theoretical \(1 : 0.5 : 0\). (S waves have a different pattern, \(\sin\theta\cos\theta\), which the same machinery also reproduces automatically – it is all just \(\hat t \cdot \nabla v \cdot \hat t\).)

Recipe: DAS with the elastic propagator

The fibre is described by an ordered [n, ndim] tensor of points in grid cells (a polyline; one DAS channel per point). deepwave.das.fiber_gauge_locations() turns it into per-channel gauge endpoint pairs and unit tangents; the endpoints become float-coordinate velocity receiver locations (each spatial component recorded at the same physical points); deepwave.das.strain_rate_from_velocities() combines the recorded velocities into strain rate:

import torch
import deepwave
from deepwave import das

# Fibre: a straight cable at depth, along x (points in grid cells).
fiber_points = torch.stack(
    [
        torch.full((100,), 60.0),          # y (cells)
        40.0 + torch.arange(100.0),        # x (cells): one channel/cell
    ],
    dim=1,
)
gauge_length = 20.0                        # physical units (e.g. m)
plus, minus, tangents = das.fiber_gauge_locations(
    fiber_points, gauge_length, grid_spacing=5.0,
)

# Record both velocity components at both gauge ends of every channel.
locations = torch.cat([plus, minus]).unsqueeze(0)  # [1, 2n, 2], float
out = deepwave.elastic(
    lamb, mu, buoyancy, 5.0, dt,
    source_amplitudes_y=source_amplitudes,
    source_locations_y=source_locations,
    receiver_locations_y=locations,
    receiver_locations_x=locations,
    pml_freq=freq,
)

n = fiber_points.shape[0]
vy, vx = out.receiver_amplitudes_y, out.receiver_amplitudes_x
strain_rate = das.strain_rate_from_velocities(
    {"y": vy[:, :n], "x": vx[:, :n]},   # + gauge ends
    {"y": vy[:, n:], "x": vx[:, n:]},   # - gauge ends
    tangents,
    gauge_length,
)                                       # [shot, channel, time]

Each velocity component is passed under the channel name the propagators use for it ("y" for receiver_amplitudes_y, the same key that a receivers= mapping takes), so nothing depends on the order you write them in: the function matches the keys to the columns of tangents itself, using the library’s slowest-to-fastest axis order. In 3-D, also record receiver_locations_z and pass {"z": vz, "y": vy, "x": vx}. Because the velocity components live at half time steps on the staggered grid, the strain-rate samples correspond to times \((t - 0.5)\,\Delta t\), like every velocity receiver output.

Which propagators? The elastic propagator is the natural home for DAS (fibres respond to elastic ground motion, including S waves). The acoustic propagator also provides particle-velocity receivers (receiver_locations_z/_y/_x), so the identical recipe applies to it. The scalar propagator records only its scalar wavefield – there are no particle-velocity receivers to pair up – so it is not supported.

Gauge averaging as post-processing

strain_rate_from_velocities() already models the gauge physically, as the two-point difference of the interferometric measurement itself. A second, independent form of gauge averaging shows up whenever the fibre is sampled far more densely than the gauge length: interrogators commonly report each channel as a moving average over the gauge window of an underlying finer-grained per-channel signal, rather than (or in addition to) a single two-point difference spanning the whole gauge. gauge_average() implements exactly that – a windowed mean along the channel dimension of any per-channel tensor:

smoothed = das.gauge_average(
    strain_rate, gauge_length=20.0, channel_spacing=1.0, dim=-2,
)

For a window of n = round(gauge_length / channel_spacing) channels, an odd n gives an exactly centred moving average; an even n (no window of channels is exactly centred on any one channel) uses the convention left = n // 2, right = n - 1 - left, i.e. the reported value at channel i is the average actually centred half a channel toward the lower-index side (see the function’s docstring for the exact statement and the boundary/replicate-padding behaviour).

Helically wound fibres

Real DAS cables often wind the fibre helically around a central strength member rather than laying it straight. Following Kuvshinov (2016, Geophysical Prospecting 64(3), 671-688, “Interaction of helically wound fibre-optic cables with plane seismic waves”, doi:10.1111/1365-2478.12303): at wrap angle \(\alpha\) to the cable axis, the fibre’s effective axial strain rate is a weighted combination of the cable-frame normal strain rates (the keyword-only wrap_angle, \(\alpha\) in radians),

\[\dot\varepsilon_\text{helix} = \cos^2\alpha \; \dot\varepsilon_{zz} + \frac{\sin^2\alpha}{2} \left( \dot\varepsilon_{xx} + \dot\varepsilon_{yy} \right),\]

where \(z\) is the cable axis (the fibre tangent used everywhere else in this module) and \(x, y\) are any orthonormal pair spanning the transverse plane (their strain-rate sum is invariant to that choice). \(\alpha = 0\) reduces this exactly to the straight-fibre axial response above. At the magic angle \(\alpha = \arccos(1/\sqrt3) \approx 0.9553\) radians (\(54.7356^\circ\)), the three weights are each exactly 1/3, so the fibre becomes a pure volumetric-strain (dilatation) sensor – the winding angle real cables commonly use because it removes shear-wave/orientation ambiguity.

helical_fiber_response() computes \(\dot\varepsilon_\text{helix}\) from the axial strain rate (exactly strain_rate_from_velocities()’s output using the fibre tangent) and one (2-D, plane strain) or two (3-D) transverse strain rates, which are obtained with the same two-point gauge recipe, oriented normal to the fibre instead of along it – see the function’s docstring for the exact endpoint arithmetic and the derivation, and for the attribution note on SWEEP’s helical_das_response, consulted as an API cross-reference (not a code source; it hardcodes only two angles with unnormalised integer weights, whereas this implements Kuvshinov’s continuous formula).

Units

• fiber_points (and the returned endpoint locations) are in grid cells, like every location argument in Deepwave.

• gauge_length is physical (the same units as grid_spacing, e.g. metres). fiber_gauge_locations() converts it to a cell offset internally; strain_rate_from_velocities() divides the physical velocity difference by it, so with velocities in m/s and gauge_length in metres the strain rate comes out in 1/s.

• The returned tangents are physical direction cosines (anisotropic grid_spacing is accounted for).

• channel_spacing (gauge_average()) is physical, the same units as gauge_length – it is only ever used as their ratio, so any consistent physical unit works.

• wrap_angle (helical_fiber_response()) is in radians, measured from the cable axis, like every other angle in Deepwave. It is keyword-only. Convert a wrap angle you have in degrees with math.radians: the magic angle is math.acos(1 / math.sqrt(3)) (0.9553 radians, 54.7356 degrees).

Differentiability

strain_rate_from_velocities() is a composition of plain differentiable torch operations, and the velocity receivers themselves are differentiable outputs of the propagator, so the gradient of any loss on DAS strain-rate channels flows back to the material models (lamb/mu/buoyancy), the source amplitudes, and any other differentiable input – FWI with DAS data requires no additional setup, exactly like FWI on conventional receivers. This is checked against a finite difference in tests/test_das.py. gauge_average() (a windowed mean) and helical_fiber_response() (a weighted sum) are likewise plain differentiable torch operations, checked directly against their known constant gradients in tests/test_das.py.

Limitations

• Straight-segment tangents. Tangents are estimated by centred finite differences of the fiber_points polyline (one-sided at the fibre’s two ends), which assumes the fibre is locally straight between neighbouring points. Sample a curved fibre finely relative to its radius of curvature, or the tangents (and the placement of the gauge ends, which lie on the local tangent line rather than the curved fibre) become approximate.

• Gauge ends must stay in-model. No bounds checking is performed here; keep the fibre, plus half a gauge length at each end (plus the Hicks interpolation halfwidth), inside the model, or the propagator will raise.

• Strain (as opposed to strain-rate) output is not provided as a separate helper: Deepwave’s elastic receivers record particle velocity, whose gauge difference is strain rate. If your instrument delivers strain, integrate the strain rate over time (e.g. torch.cumsum(strain_rate, dim=-1) * dt), or differentiate your observed strain data instead, whichever better suits your processing.

• No 1-D support (a fibre needs a tangent direction, so at least two spatial dimensions and two fibre points).

• helical_fiber_response() neglects the cable-frame shear-strain terms, following Kuvshinov (2016): they average to zero over one wrap for a symmetric periodic helix. It also does not build the transverse gauge endpoints for you (unlike fiber_gauge_locations() for the axial case) – see its docstring for the short recipe.

API

fiber_gauge_locations Function

deepwave.das.fiber_gauge_locations(fiber_points: Tensor, gauge_length: float, grid_spacing: float | Sequence[float]) → tuple[Tensor, Tensor, Tensor]

Compute per-channel DAS gauge endpoints and tangents from fibre points.

fiber_points describes the fibre as an ordered polyline: point i and point i + 1 are adjacent along the fibre (the order matters, the absolute spacing does not need to be uniform). One DAS channel is placed at each input point (so n points give n channels, matching a common real-world convention where the channel spacing is the point spacing of fiber_points, independent of the gauge length). The channel’s tangent is estimated from its neighbours by a centred finite difference ((fiber_points[i+1] - fiber_points[i-1]) / 2); the two endpoints of the fibre use a one-sided (forward/backward) difference instead. The tangent’s sign (which neighbour is “forward”) does not affect the strain rate computed by strain_rate_from_velocities() (the projection is quadratic in the tangent), only its axis.

The two gauge endpoints for each channel are then placed at fiber_points[i] +/- (gauge_length / 2) * tangent[i], converting the physical half-gauge-length offset into a grid-cell displacement with grid_spacing.

Parameters:

• fiber_points – Ordered fibre sample points, shape [n, ndim] (n >= 2), in grid cells (float). Consecutive points must be distinct.

• gauge_length – The physical gauge length, in the same physical units as grid_spacing (e.g. metres). Must be positive.

• grid_spacing – The physical size of a grid cell: a single number for isotropic spacing, or one number per dimension.

Returns:

• plus_locations: Shape [n, ndim], grid cells: the + tangent gauge endpoint of each channel.

• minus_locations: Shape [n, ndim], grid cells: the - tangent gauge endpoint of each channel.

• tangents: Shape [n, ndim]: the unit tangent of each channel, in physical direction cosines (accounts for anisotropic grid_spacing).

Return type:

Tuple[torch.Tensor, torch.Tensor, torch.Tensor]

strain_rate_from_velocities Function

deepwave.das.strain_rate_from_velocities(v_receivers_plus: Mapping[str, Tensor], v_receivers_minus: Mapping[str, Tensor], tangents: Tensor, gauge_length: float) → Tensor

Compute DAS strain rate from paired particle-velocity receivers.

Combines particle-velocity receivers recorded at the two gauge endpoints (typically the plus_locations/minus_locations returned by fiber_gauge_locations(), fed as receiver locations to deepwave.elastic() or deepwave.acoustic(), one call per spatial component) into the tangential gauge-averaged strain rate

\[\dot\varepsilon = \frac{v_t^+ - v_t^-}{L_g}, \qquad v_t^\pm = \hat t \cdot v^\pm.\]

Each velocity component is passed under the channel name the propagators use for it ("y" for receiver_amplitudes_y, the same key a receivers= mapping takes), so which component is which is stated rather than inferred from position: the columns of tangents are the spatial dimensions in the library’s slowest-to-fastest order, so the keys are "y" and "x" in 2-D and "z", "y" and "x" in 3-D, and this function matches them up itself.

Parameters:

• v_receivers_plus – The particle-velocity receiver amplitudes recorded at the + tangent gauge endpoints, as a mapping from velocity channel name to a tensor of shape [shot, channel, time]: one entry per spatial dimension, {"y": ..., "x": ...} in 2-D and {"z": ..., "y": ..., "x": ...} in 3-D.

• v_receivers_minus – The matching receivers at the - tangent gauge endpoints, with the same keys and shapes as v_receivers_plus.

• tangents – Unit tangents, shape [channel, ndim] (as returned by fiber_gauge_locations()).

• gauge_length – The physical gauge length used to place the endpoints (same value passed to fiber_gauge_locations()). Must be positive.

Returns:

Strain rate, shape [shot, channel, time].

Return type:

torch.Tensor

Raises:

• DeepwaveTypeError – If either velocity argument is not a mapping (for example a sequence of tensors).

• DeepwaveValueError – If the keys are not exactly the velocity channels of a run with tangents.shape[-1] spatial dimensions, if gauge_length is not positive, or if the receiver tensors do not all have the same [shot, channel, time] shape with channel matching tangents.shape[0].

gauge_average Function

deepwave.das.gauge_average(data: Tensor, gauge_length: float, channel_spacing: float, dim: int = -2) → Tensor

Apply an interrogator’s finite-gauge moving average to per-channel data.

strain_rate_from_velocities() already models the gauge physically, as a two-point difference (the interferometric measurement itself). This function models a second, independent source of gauge averaging that shows up whenever a fibre is sampled far more densely than the gauge length: many interrogators report each channel as a moving average over the gauge window of an underlying finer-grained per-channel signal (e.g. many closely spaced two-point strain-rate estimates, or any other per-channel quantity), rather than (or in addition to) a single two-point difference spanning the whole gauge. Applying this to the output of strain_rate_from_velocities(), computed with a small internal two-point gauge on densely spaced channels, approximates that measurement.

Given a window of n = round(gauge_length / channel_spacing) channels (clamped to at least 1; ties round to even, Python’s round()), the output at channel i is the mean of the n channels [i - left, i + right] (channels beyond the array ends are replicated from the nearest in-bounds channel, so the array length and channel positions are unchanged):

\[\bar d_i = \frac{1}{n} \sum_{k=-\text{left}}^{\text{right}} d_{i+k}.\]

Even/odd window convention. For odd n (left == right == (n - 1) / 2), the window is exactly centred on channel i: an unbiased centred moving average. For even n, no window of channels is exactly centred on any one channel (the true window centre falls half a channel spacing away, between two channels); this function uses left = n // 2, right = n - 1 - left, i.e. the window includes one more channel on the lower-index side than the higher-index side. Equivalently, the value reported at channel i is the average actually centred at position i - 0.5 channels. If this half-channel offset matters for your application (e.g. exact two-way time picks), shift your channel coordinates by -0.5 * channel_spacing when n is even, or choose gauge_length/channel_spacing so that n is odd.

Parameters:

• data – Per-channel data of any shape, e.g. [shot, channel, time] (the default dim=-2 matches the channel axis of strain_rate_from_velocities()’s output).

• gauge_length – The physical gauge length (same physical units as channel_spacing, e.g. metres). Must be positive.

• channel_spacing – The physical spacing between adjacent channels (same physical units as gauge_length). Must be positive.

• dim – The channel dimension of data to average over. Defaults to -2.

Returns:

Gauge-averaged data, the same shape and dtype as data.

Return type:

torch.Tensor

helical_fiber_response Function

deepwave.das.helical_fiber_response(axial_strain_rate: Tensor, transverse_strain_rates: Sequence[Tensor], *, wrap_angle: float) → Tensor

Combine cable-frame normal strain rates into a helical-fibre response.

Real DAS cables often wind the fibre helically around a central strength member rather than laying it straight; the fibre’s local strain then responds to a weighted combination of the cable’s normal strain-rate components, not just the axial one. Following Kuvshinov (2016, Geophysical Prospecting 64(3), 671-688, “Interaction of helically wound fibre-optic cables with plane seismic waves”, doi:10.1111/1365-2478.12303), for a fibre wound at wrap angle \(\alpha\) to the cable axis, in radians (\(\alpha = 0\) is a straight fibre along the cable; \(\alpha = \pi/2\) is a fibre wound purely circumferentially), and neglecting the shear-strain terms that time-average to zero for a symmetric periodic helix, the effective axial strain rate the fibre reports is

\[\dot\varepsilon_\text{helix} = \cos^2\alpha \; \dot\varepsilon_{zz} + \frac{\sin^2\alpha}{2} \left( \dot\varepsilon_{xx} + \dot\varepsilon_{yy} \right),\]

where \(z\) is the cable’s local axial direction (the deployment tangent \(\hat t\) used throughout this module) and \(x, y\) are any orthonormal pair spanning the plane normal to it – the transverse-plane sum \(\dot\varepsilon_{xx} + \dot\varepsilon_{yy}\) is invariant to that choice (it equals \(\operatorname{tr}(\dot\varepsilon) - \dot\varepsilon_{zz}\), the strain-rate tensor’s trace restricted to the transverse plane, which does not depend on the transverse basis), so any orthonormal transverse pair gives the same answer. The three weights sum to 1 (\(\cos^2\alpha + \sin^2\alpha/2 + \sin^2\alpha/2 = 1\)): this is always a weighted average of the three cable-frame normal strain rates.

\(\alpha = 0\) reduces this exactly to strain_rate_from_velocities()’s straight-fibre axial response (weight 1 on the axial term, 0 on the transverse terms). At the “magic angle” \(\alpha = \arccos(1/\sqrt3) \approx 0.9553\) radians (\(54.7356^\circ\)), \(\cos^2\alpha = \sin^2\alpha / 2 = 1/3\): the three weights are exactly equal, so the fibre responds only to the trace of the strain-rate tensor, i.e. it becomes a pure volumetric-strain (dilatation) sensor, insensitive to shear – the winding angle real interrogator cables commonly use because it removes shear-wave/ orientation ambiguity.

How to obtain the inputs from Deepwave’s velocity receivers:

axial_strain_rate is exactly the output of strain_rate_from_velocities(), using the real fibre’s tangent (fiber_gauge_locations()’s tangents), as in this module’s main recipe.

Each entry of transverse_strain_rates uses the same two-point gauge construction, but oriented along a direction \(\hat n\) orthogonal to the tangent instead of along it: given the channel points and a unit \(\hat n\) in physical direction-cosine space (orthogonal to tangents),

half_offset_cells = gauge_length / 2 * n_hat / grid_spacing
n_plus = channel_points + half_offset_cells
n_minus = channel_points - half_offset_cells

(the same arithmetic fiber_gauge_locations() uses internally for the tangent) – record velocities at n_plus/n_minus and call strain_rate_from_velocities() with \(\hat n\) in place of tangents. Because \(\dot\varepsilon_{nn}\) is quadratic in \(\hat n\), its sign does not matter.

• 2-D (elastic/acoustic with ndim=2): there is only one in-plane transverse direction (rotate the physical tangent (t0, t1) -> (-t1, t0), a 90-degree rotation that preserves unit length regardless of grid anisotropy since it is done in physical space). The 2-D propagators enforce plane strain (no out-of-plane motion), so the missing out-of-plane transverse component is exactly 0 – pass a length-1 sequence containing just the in-plane transverse strain rate.

• 3-D: pass both orthogonal transverse strain rates as a length-2 sequence (any orthonormal pair spanning the plane normal to the tangent; see the trace argument above for why the specific choice does not matter).

Note

This function implements Kuvshinov’s continuous formula for general \(\alpha\) from first principles. SWEEP’s helical_das_response (github.com/DeepWave-KAUST/sweep, MIT licence, src/sweep/equations/das.py, commit 9eff76deadabcc7c84de761f35ad42637e5a959e) was consulted as a cross-reference for the overall API shape (normal strain-rate components + wrap angle in, one fibre-response tensor out), but it hardcodes only its two named angles with unnormalised integer weights (1:1:1 at its “35.3 deg”, 4:1:1 at its “54.7 deg”, i.e. not divided by the sum of weights); no code from it is used here.

Parameters:

• axial_strain_rate – The cable-axial normal strain rate \(\dot\varepsilon_{zz}\), any shape (typically [shot, channel, time]).

• transverse_strain_rates – A sequence of 1 (2-D, plane strain) or 2 (3-D) transverse normal strain-rate tensors, each the same shape as axial_strain_rate.

• wrap_angle – The helical wrap angle \(\alpha\), in radians (like every other angle in Deepwave), measured from the cable axis: 0 is a straight fibre, math.pi / 2 is purely circumferential winding, and the magic angle is math.acos(1 / math.sqrt(3)). Keyword-only. An angle in degrees is converted with math.radians.

Returns:

The helical-fibre effective strain rate \(\dot\varepsilon_\text{helix}\), the same shape as axial_strain_rate.

Return type:

torch.Tensor

Raises:

DeepwaveValueError – If transverse_strain_rates does not contain 1 or 2 tensors, or if one of them does not have the same shape as axial_strain_rate.