Sensitivity kernels for transmission fibre optics

SUMMARY

Fibre-optic sensing based on transmission offer an alternative to scattering-based distributed acoustic sensing (DAS). The ability to interrogate fibres that are thousands of kilometres long opens opportunities for studies of remote regions, including ocean basins. However, by averaging deformation along the fibre, transmission systems produce integrated instead of distributed measurements. They defy traditional interpretations in terms of simple seismic phases, thereby inherently requiring a full-waveform approach. For this, we develop a formalism to calculate sensitivity kernels of transmitted optical phase changes with respect to (Earth) structure using optical phase delay measurements. We demonstrate that transmission-based sensing can effectively provide distributed measurements when optical phase delays are analysed in different time windows. The extent to which a potentially useful sensitivity coverage can be achieved depends on the fibre geometry, and specifically on its local curvature. This work establishes a theoretical foundation for tomographic inversions and experimental design using transmission-based optical sensing.

Computational seismology, Seismic instruments, Seismic tomography, Theoretical seismology, Wave propagation

1 INTRODUCTION

Distributed acoustic sensing (DAS) is an emerging family of technologies that permits seismic data acquisition with high spatio-temporal resolution and large bandwidth using fibre-optic cables (e.g. Hartog 2017; Lindsey et al. 2020; Paitz et al. 2021). The relative ease of deploying cables in challenging terrain, and opportunities to co-use existing telecommunication infrastructure, have enabled seismological applications that would have been substantially more difficult or expensive using conventional seismic instruments (e.g. Lindsey et al. 2017; Martin et al. 2017; Walter et al. 2020; Klaasen et al. 2021). While Rayleigh back-scattering allows DAS to achieve distributed measurements with an effective channel spacing in the centimetre range, it also limits the length of the fibre that can be interrogated to several tens of kilometres.

Alternative systems based on optical transmission overcome this limitation by measuring deformation-induced changes of phase (Marra et al. 2018; Bogris et al. 2021, 2022; Bowden et al. 2022) or polarization (Mecozzi et al. 2021). Reaching interrogation distances of thousands of kilometres, opens new opportunities to study seismic activity and Earth structure in remote regions, including ocean basins. The main drawback of transmission-based systems lies in the averaging of deformation along the fibre. Hence, in contrast to DAS, the measurement is not distributed but integrated. It therefore remains unclear if transmission measurements could be used at all to infer Earth structure.

In the following sections, we develop a formalism for the calculation of sensitivity kernels for measurements of deformation-induced optical phase delays. Using simple conceptual examples, we demonstrate that dissecting the phase delay time-series can effectively enable distributed measurements, thereby providing a theoretical foundation for tomographic inversions and optimal experimental design.

2 THEORETICAL DEVELOPMENTS

We consider a fibre with position |$\mathbf {x}=\hat{\mathbf {x}}(s)$| parametrized in terms of the arc length s ∈ [0, L], where L is the total length, as illustrated in Fig. 1. In terms of the effective refractive index r(s), the vacuum speed of light c, and the circular frequency ω of the optical signal, its time-dependent phase change |$\theta (t)=\Delta \dot{\varphi }(t)$| can be represented in two equivalent forms (Bowden et al. 2022; Fichtner et al. 2022),

$$\begin{equation*} \theta (t) = \frac{\omega }{c} \int _{s=0}^L r(s) \, \dot{\varepsilon }[\hat{\mathbf {x}}(s),t]\, {\rm d}s = - \frac{\omega }{c} \int _{s=0}^L \frac{{\rm d}}{{\rm d}s} \left[ r(s)\, \mathbf {e}(s) \right] \cdot \dot{\mathbf {u}}[\hat{\mathbf {x}}(s),t]\, {\rm d}s\, , \end{equation*}$$

(1)

where an overdot denotes a time derivative. Eq. (1) relates the optical phase change θ(t) to the deformation of the medium, expressed in terms of the displacement field |$\mathbf {u}(\mathbf {x},t)$| or the axial component |$\varepsilon =\mathbf {e}\cdot \mathbf {E} \cdot \mathbf {e}$| of the strain tensor |$\mathbf {E}=(\nabla \mathbf {u}+\nabla \mathbf {u}^T)/2$|⁠. The first variant of eq. (1) allows us to compare transmission measurements θ(t) to axial strain rate measurements |$\dot{\varepsilon }[\hat{\mathbf {x}}(s),t]$| from DAS. The second variant of eq. (1), which we will use in the following paragraphs, exhibits the dependence of θ(t) on fibre curvature and spatial variations of the effective refractive index.

Figure 1.

Schematic illustration of fibre deformation. The undeformed fibre, shown as black curve, is represented by the position vector |$\hat{\mathbf {x}}(s)$|⁠, which is parametrized in terms of the arc length s ∈ [0, L]. The displacement field |$\mathbf {u}[\hat{\mathbf {x}}(s),t]$| in blue moves |$\hat{\mathbf {x}}(s)$| to |$\hat{\mathbf {x}}(s)+\mathbf {u}[\hat{\mathbf {x}}(s)]$|⁠. The result is the deformed fibre in grey. The local tangent vector |$\mathbf {e}(s)$| is displayed as a thick black arrow. Upon deformation, an optical signal, shown in orange, acquires a phase change Δϕ(t).

The dynamic fields |$\mathbf {u}$| and ε depend on the distribution of medium properties, such as wave speeds and density, collected in the model vector |$\mathbf {m}(\mathbf {x})$| and omitted in the notation to avoid clutter. To infer a plausible |$\mathbf {m}(\mathbf {x})$| from the comparison of observed and calculated phase changes, θ^obs and θ, we require a suitable measurement functional. Among numerous available options (e.g. Gee & Jordan 1992; Fichtner et al. 2008; Bozdağ et al. 2011) we choose the cross-correlation time-shift within a windowed section of the phase change time-series (Luo & Schuster 1991). We define the windowed cross-correlation as

$$\begin{equation*} \mathcal {C}(\tau )=\int _{t=-\infty }^{\infty } w(t) \, \theta ^\text{obs}(t)\, \theta (t+\tau )\, {\rm d}t\, , \end{equation*}$$

(2)

where w(t) denotes a window function that isolates a selected part of the time-series. The time-shift |$\mathcal {T}$| between the windowed versions of θ^obs and θ is defined as the value of τ where |$\mathcal {C}(\tau )$| reaches its global maximum. Hence, by differentiating (2), we obtain

$$\begin{equation*} 0 = \dot{\mathcal {C}}(\mathcal {T}) = \int _{t=-\infty }^{\infty } w(t) \, \theta ^\text{obs}(t)\, \dot{\theta }(t+\mathcal {T})\, {\rm d}t = - \int _{t=-\infty }^{\infty } \dot{\theta }_w^\text{obs}(t-\mathcal {T})\, \theta (t)\, {\rm d}t\, , \end{equation*}$$

(3)

with the convenient definition of the windowed phase change |$\theta _w^\text{obs}=w\theta ^\text{obs}$|⁠. Eq. (3) defines the measurement functional |$\mathcal {T}$| implicitly. Variations |$\delta \mathbf {m}(\mathbf {x})$| of medium parameters induce variations |$\delta \mathcal {T}$| of the measurement. Invoking implicit function differentiation, we obtain

$$\begin{equation*} \delta \mathcal {T} = - \int _{t=-\infty }^{\infty } \dot{\theta }_w^\text{obs}(t-\mathcal {T})\, \delta \theta (t)\, {\rm d}t \, / \int _{t=-\infty }^{\infty } \ddot{\theta }_w^\text{obs}(t-\mathcal {T})\, \theta (t)\, {\rm d}t \, . \end{equation*}$$

(4)

Making the common assumption that |$\theta _w^\text{obs}$| is approximately a time-shifted version of θ_w, that is |$\theta _w^\text{obs}(t-\mathcal {T})\approx \theta _w(t)$| (e.g. Luo & Schuster 1991; Dahlen et al. 2000), we can simplify eq. (4) to

$$\begin{equation*} \delta \mathcal {T} = ||\dot{\theta }||_w^{-2} \, \int _{t=-\infty }^{\infty } \dot{\theta }_w(t)\, \delta \theta (t)\, {\rm d}t\, , \end{equation*}$$

(5)

with the squared norm |$||\dot{\theta }||_w^{2}=\int _{t=-\infty }^{\infty } w(t)\, \dot{\theta }^2(t)\, {\rm d}t$|⁠. To relate variations of the measurement functional |$\delta \mathcal {T}$| to variations in medium parameters |$\delta \mathbf {m}(\mathbf {x})$|⁠, we substitute the variation of eq. (1) into (5),

$$\begin{equation*} \delta \mathcal {T} = - ||\dot{\theta }||_w^{-2} \, \int _{t=-\infty }^{\infty } \int _{s=0}^L \ddot{\theta }_w(t) \, \mathbf {a}(s) \cdot \delta \mathbf {u}[\hat{\mathbf {x}}(s),t]\, {\rm d}s\, , \end{equation*}$$

(6)

where we defined the vector |$\mathbf {a}(s)=(\omega /c)\, d\left[ r(s)\, \mathbf {e}(s) \right]/{\rm d}s$| for notational convenience. Since |$\delta \mathbf {u}$| can in practice not be computed efficiently for large numbers of arbitrary perturbations |$\delta \mathbf {m}(\mathbf {x})$|⁠, we apply the adjoint method (e.g. Tromp et al. 2005; Fichtner et al. 2006). For this, we write the (seismic) wave equation symbolically in terms of the wave equation operator |$\mathbf {L}$| as |$\mathbf {L}[\mathbf {u}(\mathbf {m}),\mathbf {m}]=\mathbf {f}$|⁠, where |$\mathbf {f}$| is the source term. Multiplying the variation of the wave equation, |$\delta \mathbf {L}[\mathbf {u}(\mathbf {m}),\mathbf {m}] + \mathbf {L}[\delta \mathbf {u}(\mathbf {m}),\mathbf {m}]=\mathbf {0}$| with an arbitrary field |$\mathbf {u}^\dagger$| and integrating over time and the spatial domain |$\oplus \subset \mathbb {R}^3$|⁠, yields

$$\begin{equation*} \int _{t=-\infty }^{\infty } \int _\oplus \mathbf {u}^\dagger \cdot \delta \mathbf {L}[\mathbf {u}(\mathbf {m}),\mathbf {m}] \, {\rm d}\mathbf {x}\, {\rm d}t + \int _{t=-\infty }^{\infty } \int _\oplus \mathbf {u}^\dagger \cdot \mathbf {L}[\delta \mathbf {u}(\mathbf {m}),\mathbf {m}]\, {\rm d}\mathbf {x}\, {\rm d}t = 0\, . \end{equation*}$$

(7)

Invoking the adjoint |$\mathbf {L}^\dagger$| of |$\mathbf {L}$| and adding eq. (6) to eq. (7), we find

$$\begin{equation*} \delta \mathcal {T} = \int _{t=-\infty }^{\infty } \int _\oplus \delta \mathbf {u} \cdot \left[ \mathbf {L}^\dagger (\mathbf {u}^\dagger ,\mathbf {m}) - ||\dot{\theta }||_w^{-2} \ddot{\theta }_w(t) \, \mathbf {a}(s) \, \delta _L(\mathbf {x}) \right]\, {\rm d}\mathbf {x}\, {\rm d}t + \int _{t=-\infty }^{\infty } \int _\oplus \mathbf {u}^\dagger \cdot \delta \mathbf {L}[\mathbf {u}(\mathbf {m}),\mathbf {m}]\, {\rm d}\mathbf {x}\, {\rm d}t\, . \end{equation*}$$

(8)

In (8), the curve δ-distribution δ_L is defined as |$\int _{s=0}^L f[\hat{\mathbf {x}}(s)]\, {\rm d}s = \int _\oplus f(\mathbf {x})\, \delta _L(\mathbf {x})\, {\rm d}\mathbf {x}$| for any function f in ⊕. We can now eliminate |$\delta \mathbf {u}$| from (8) by forcing the first term on the right-hand side to zero and thereby defining the adjoint field |$\mathbf {u}^\dagger$| as the solution of the adjoint equation

$$\begin{equation*} \mathbf {L}^\dagger (\mathbf {u}^\dagger ,\mathbf {m}) = ||\dot{\theta }||_w^{-2} \ddot{\theta }_w(t) \, \mathbf {a}(s) \, \delta _L(\mathbf {x})\, . \end{equation*}$$

(9)

The right-hand side of eq. (9) is the adjoint source, which δ_L localizes along the fibre. The time evolution of the adjoint source is controlled by the second derivative of the windowed computed phase changes, |$\ddot{\theta }_w$|⁠, and its position-dependent orientation is given by the vector |$\mathbf {a}(s)$|⁠. What is left of eq. (8) can be conveniently written in terms of the structural sensitivity kernel |$\mathbf {K}$|⁠, that is the volumetric sensitivity of the time-shift |$\mathcal {T}$|⁠,

$$\begin{equation*} \delta \mathcal {T} = \int _{t=-\infty }^{\infty } \int _\oplus \mathbf {u}^\dagger \cdot \delta \mathbf {L}[\mathbf {u}(\mathbf {m}),\mathbf {m}]\, {\rm d}\mathbf {x}\, {\rm d}t\, = \int _\oplus \mathbf {K}(\mathbf {x})\cdot \delta \mathbf {m}\, . \end{equation*}$$

(10)

Explicit expressions of |$\mathbf {K}$| for various medium parameters can be found, for example in Tromp et al. (2005) or Fichtner (2010), and will not be repeated here.

3 CONCEPTUAL EXAMPLES

The following examples are intended to illustrate the generation of phase change signals θ(t) and their time-dependent sensitivity to medium parameters. We make the plausible assumption that the effective refractive index r(s) is constant over a seismic wavelength, that is |$\mathcal {O}(10)$| km in our examples, thereby allowing us to ignore its derivative. Not trying to mimic a specific acquisition system, we set ωr/c = 1 m s^–2 for simplicity. The normalization factor |$||\dot{\theta }||_w^{-2}$| in the adjoint source (9) ensures that traveltime kernels are unaffected by this choice, and by any other amplitude scaling, for example via the seismic moment of the wavefield source. With this setting, the vector |$\mathbf {a}(s)$| equals the non-normalized normal vector |${\rm d}\mathbf {e}(s)/{\rm d}s$|⁠. To avoid any complications, the elastic medium for our calculations is unbounded, isotropic and perfectly elastic, with P velocity α = 8000 m s^–1, S velocity β = 5000 m s^–1 and density ρ = 3000 kg m^–3. Well-known analytical solutions for moment tensor and single force sources may be found, for example in Aki & Richards (2002).

Starting with the simplest possible case, we consider an explosive source that only radiates a P wave, as shown in Fig. 2. The displacement field |$\mathbf {u}(\mathbf {x},t)$| deforms two different fibres, plotted as black and red curves, respectively. Since |$\mathbf {u}(\mathbf {x},t)$| interacts with different parts of the fibres at different times, the resulting phase changes θ(t) contain multiple oscillations that are more complex than the seismic P wavelet. Even though the red fibre is only half as long as the black fibre, it produces phase changes of nearly the same amplitude, because θ(t) is proportional to |$\mathbf {a}(s)={\rm d}\mathbf {e}(s)/{\rm d}s$|⁠, that is to curvature. Hence, the example in Fig. 2 illustrates that larger curvature may compensate shorter length, and vice versa.

Figure 2.

Forward modelling example using two different fibre geometries with lengths of L = 1989 km (black) and L = 1065 km (red). The wavefield radiates from an explosive source at the location of the star. A snapshot of the x-component is shown in the left-hand panel, together with the source–time function; a Heaviside function filtered between 0.1 and 0.5 Hz. The resulting time-series θ(t) are shown to the right-hand panel, in the colour of the corresponding fibre geometry.

Applying the formalism developed in Section 2, allows us to compute sensitivity kernels K_α for relative perturbations of P velocity, δln α, some examples of which are displayed in Fig. 3. Using different windows w(t) for the time-shift measurements |$\mathcal {T}$|⁠, produces kernels with different spatial coverage. They exhibit a Fresnel zone structure with vanishing sensitivity along the ray path, known from finite-frequency traveltime measurements on seismometer recordings (e.g. Dahlen et al. 2000). In accord with eq. (1), the kernels connect the source to segments of the fibre with large curvature.

Figure 3.

Sensitivity kernels K_α for relative P velocity perturbations δln α. Time-shift measurements |$\mathcal {T}$| in different windows, shown in the lower panel, produce kernels with different spatial coverage that indicate the section of the fibre where the signal has been primarily generated.

Modifying the source from an explosion to a double-couple moment tensor generates S waves in addition to P waves. As illustrated in Fig. 4, they also appear in the phase change time-series θ(t) as a sequence of oscillations produced by different fibre segments. The fibre segment responsible for a specific oscillation can again be determined by sensitivity kernel analysis, now providing sensitivity with respect to both P and S velocity. Fig. 4 corroborates that segments with large curvature behave similar to discrete, that is distributed measurement points, thereby providing window-specific sensitivity coverage, despite the integrated nature of the phase change time-series.

Figure 4.

Sensitivity kernels K_α and K_β for relative velocity perturbations δln α and δln β, respectively. The phase change time-series θ(t) and the different measurement windows used to calculate kernels are shown in the lower panel.

4 DISCUSSION AND CONCLUSIONS

Using adjoint techniques for structural sensitivity analysis, we have shown with conceptual examples that optical sensing systems based on the transmission of deformation-induced phase changes can effectively be used to make space-distributed measurements. For this, the fibre must contain curved segments that behave similar to localized sensors because the sensitivity of a fibre segment to deformation is proportional to its local curvature |${\rm d}\mathbf {e}(s)/{\rm d}s$|⁠. Provided that the spacing of curved segments is larger than a seismic wavelength, they can produce a sequence of distinguishable wavelets in the phase change time-series θ(t), each representing local, that is distributed, deformation. According to eq. (1), the fibre normal vector |${\rm d}\mathbf {e}(s)/{\rm d}s$| should be roughly parallel to the displacement field polarization |$\mathbf {u}$|⁠.

In the presence of a suitably shaped fibre, the time-dependent analysis of θ(t) may provide a set of traveltime (or other) measurements and sensitivity kernels that are useful for the solution of tomographic inverse problems. In this context, strongly curved fibre segments roughly mimic a network of conventional seismic instruments. This may, indeed, be beneficial for imaging remote regions such as ocean basins, where telecommunication cables already exist. For this, however, the fibre geometry, that is, the instrument response of the system, needs to be known with sufficient accuracy, where the meaning of ‘sufficient’ is application-dependent.

It must, unfortunately, be suspected that developers of future telecommunications infrastructure might ignore the wishes of structural seismologists and deploy cables in accord with economic and political boundary conditions, instead of maximizing curvature. Nevertheless, the methods presented in Section 2 provide the theoretical foundation for optimal experimental cable design, as well as tools for the sensitivity analysis of existing fibre-optic cable installations.

ACKNOWLEDGEMENTS

The authors would like to thank the editor in chief Jörg Renner and the reviewer Yann Capdeville for their constructive comments. Andreas Fichtner gratefully acknowledges United Airlines for a 7-hr delay at New Orleans airport, which provided ample time to develop the forward modelling theory, summarized in eq. (1). This work was partially funded by the RISE project under the European Union’s Horizon 2020 program (grant 821115).

DATA AVAILABILITY STATEMENT

No data have been used for this research.