DAS Over Multimode Fibers With Reduced Fading by Coherent Averaging of Spatial Modes

We investigate the performance of distributed acoustic sensing over multi-mode fibers based on heterodyne phase-sensitive optical time-domain reflectometry. We report a mathematical model describing the relation between phase variation and applied strain in the presence of multi-mode propagation that supports the feasibility of distributed acoustic measurements over multi-mode fibers. We also propose a novel coherent averaging method that achieves up to a three-fold reduction of the noise floor compared to state-of-the-art methods.

Abstract-We investigate the performance of distributed acoustic sensing over multi-mode fibers based on heterodyne phase-sensitive optical time-domain reflectometry. We report a mathematical model describing the relation between phase variation and applied strain in the presence of multi-mode propagation that supports the feasibility of distributed acoustic measurements over multi-mode fibers. We also propose a novel coherent averaging method that achieves up to a three-fold reduction of the noise floor compared to state-of-the-art methods.

I. INTRODUCTION
R AYLEIGH-BASED fiber-optic distributed acoustic sensors (DAS) have attracted increasing attention in recent years due to their unique features [1]. While most DAS systems use single-mode fibers (SMFs), there is a growing interest in extending DAS to multi-mode fibers (MMFs). One of the main driving interests is to leverage the existing MMF infrastructure in the oil and gas industries, which is currently used for Raman-based temperature monitoring [2]. A few works started investigating DAS over MMFs, however, primarily focusing on few-mode fibers and not analyzing how the proposed method would scale to MMFs [3], [4], [5]. Regardless of whether the optical fiber is single-or multimode, phase-sensitive optical time domain reflectometry (φ-OTDR) is one of the most commonly used DAS configurations [1]. In recent years, a lot of effort has been put to deal with the signal fading issue of φ-OTDR systems, which is due to both intra-pulse coherent interference and polarization variation along the fiber. Recent techniques have tackled the coherent fading issue combining independent measurements channels (e.g., extracted from different frequency bands) carrying the same disturbance information, which are coherently averaged using techniques such as the so called rotatedvector-sum (RVS) [6], [7]. When DAS is applied to MMFs, mode coupling is an additional issue that needs to be faced; nevertheless, the spatial modes can be used also to reduce signal fading, since they act as independent measurement channels.
In this letter, expanding a preliminary analysis [8], we experimentally and theoretically analyze the feasibility of DAS measurements over MMFs. The proposed setup relies on a three-mode photonic lantern (PhL) as a mode de/multiplexer to exploit the spatial diversity to improve the measurement reliability [9], [10]. Moreover, we also propose a novel coherent averaging method that is optimal in terms of mean-squared error (MSE) and outperforms the RVS method.

II. THEORETICAL MODEL
Fundamentally, DAS in SMFs is made possible by the fact that the phase of the propagating light varies linearly with the strain applied to the fiber. This holds true also for MMFs, but only for each mode individually; differently, the relation between the phase of the total light field and the applied strain may be nonlinear in general. The phase variations measured by a DAS reflect what happens in the perturbed fiber sections, whose length, L, is typically shorter than the DAS gauge length and hence in the order of meters. In a MMF, the electric field transmitted across one of these sections can be written as E out = m a m exp {− j (2π/λ)n m L} E m , where the sum extends over the propagating modes, n m is the effective index of mode m, E m its field distribution (we omitted the transverse coordinates for brevity), and a m its complex amplitude. The perturbed fiber sections are typically no longer than a few meters and in such a short distance the action of the strain is not strong enough to induce a significant transfer of power from one mode to another. In a SMF scenario, this power transfer would correspond to a polarization change, which is indeed largely and safely ignored in DAS systems. For these reasons, the coefficients a m can be assumed constant.
Consistently with the typical assumption made for SMFbased DAS, we assume that the strain is constant over the perturbed fiber section. Therefore, the phase term of E out varies from (2π/λ)n m L to (2π/λ)n m L + φ m , with where ϵ is the applied axial strain, which induces the elongation L → L(1 + ϵ) and the refractive index change n m → (n m + δn m ); moreover, we set n m = n 0 + n m , with n 0 denoting the effective index of the fundamental mode (used here as a reference) and n m is the modal birefringence between the fundamental mode and the mode m [11].  be approximated to the first order in the modal birefringence as δn m ≈ −0.1 n 3 0 ϵ(1 + 3 n m /n 0 ). Inserting this expression into (1), and considering only terms up to the first order in ϵ, where the first term is the usual strain-induced phase shift observed in SMFs [1], whereas the second term is an extra phase shift characteristic of each higher-order mode. It is then clear that the relation between the phase of the transmitted field E out and the applied strain is linear only as long as this extra phase shift is much smaller than both the first term (but this is always verified since n m ≪ n 0 ) and 2π. Due to the round-trip, the phase shift measured by the DAS is actually twice the one calculated before. Therefore, the above condition yields This theoretical result sets an upper bound to the strain applied to the fiber, above which DAS in MMFs would lose linearity. This bound is inversely proportional to the perturbed fiber section. Moreover, it is also inversely proportional to the modal birefringence, which means that MMFs with less modal dispersion, such as graded-index ones, enable measurements on a wider strain range. Considering as an example L = 2 m, n 0 ≈ 1.47, and assuming max{ n m } ≈ 10 −3 , (2) yields ϵ ≪ 1.1 mϵ, which is quite a loose constraint considering the typical strain values of the perturbations encountered in DAS measurements [12].

A. Experimental Setup
To validate the feasibility of performing DAS measurements over MMFs, we implemented the heterodyne φ-OTDR scheme shown in Fig. 1. The output of a narrow-linewidth laser is split by a 90:10 coupler into two branches. One branch acts as the local-oscillator (LO); the other branch is modulated by an acousto-optic modulator (AOM) into rectangular pulses of 100 ns duration with a 110 MHz frequency shift. The pulses are launched into a 3-km-long OM2 MMF fiber (50/125 µm) through one of the ports of a three-mode PhL (Phoenix Photonics 3-PL), which acts also as a spatial demultiplexer for the Rayleigh backscattered light. The PhL is realized by a few mode fiber supporting the LP 01 , LP 11a , and LP 11b modes; the connection between the PhL and the OM2 fiber is realized through a splice, which is adiabatically tapered to minimize its insertion losses and maximize its return losses. In principle, the three spatial modes demultiplexed by the PhL should be sent to three different receivers. For simplicity, in this experiment the three modes have been multiplexed in time by means of fiber delay lines, so to use a single receiver; this comes at the expense of a three-fold reduction of the achievable acoustic bandwidth. The lengths of the fiber delay lines are about 7 km and 14 km, respectively, so that the introduced delays avoid the overlap of the backscattered traces from the 3-km-long MMF link demultiplexed by the PhL. The probe pulse repetition rate is consequently set to T rep = 110 µs. The beating between the LO and the backscattered light is detected by a balancedphotodiode (BPD) with a 400 MHz bandwidth and digitized at 500 MS/s by a 12-bit data acquisition (DAQ) board. The inphase (I) and quadrature (Q) components of the backscattered field are obtained with conventional digital I/Q demodulation.

B. LMS Coherent Averaging
The information obtained from independent measurements of the same event needs to be optimally combined, or averaged, into a single result. Let r m (t n , z k ) denote the complex value of the Rayleigh backscattered trace of the spatial mode m, at an arbitrary position z k , and at time t n . We call a "sequence" the set of values that r m (t n , z k ) takes for specific mode and position over N T consecutive time samples. Sequences from two different modes but at the same position are affected by the same acoustic perturbation; nonetheless, they have an arbitrary, unpredictable and mode-dependent phase offset, θ m , that prevent their direct averaging. Therefore, before the average can be performed, the phases of the sequences must be aligned. The RVS method achieves this alignment by selecting a reference time instant t 1 , and determining the phase offsets of each spatial mode as θ m = arg {r m (t 1 )} [6]. Nevertheless, no indication is given for the selection of this reference time, and, indeed, a poor choice compromises the phase alignment and subsequent averaging. To deal with this issue, the recently proposed double-averaging RVS (DA-RVS) method uses N T consecutive time samples of a randomly preselected sequence, and use them to compute a refined time-domain-averaged estimate of the phase offsets [13]. The method proposed in this letter proceeds in two steps. First, the sequence with the highest energy over the N T time samples is selected as the reference sequence. Then, the phase offsets are computed solving a least-mean-square (LMS) problem that yields the optimum phase offsets values that align each other sequence to the reference one. We refer to the proposed method as the LMS method. In details, assume that r (t n ) = argmax m N T n=1 |r m (t n )| 2 is the reference sequence according to the highest-energy criterion; we determine the phase offset needed to align the generic sequence r m (t n ) tor (t n ) by seeking the angle θ m that minimizes the cost function E = n w(t n )|r (t n ) − r m (t n ) exp( jθ m )| 2 , where n = 1, . . . , N T , and w(t n ) = |r (t n )||r m (t n )| is used to enhance the weight of the higher amplitude samples. Solving the relative LMS problem yields the optimal estimate of the phase θ m as θ m = arg N T n=1 w(t n )r * (t n )r m (t n ) . We remark that this choice is optimal in the MSE sense. After estimating the phase offsets, θ m , with either the RVS, DA-RVS, or the LMS methods, we compute the coherent average as follows: Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.

A. Spatial Modes Averaging
To determine the advantages of using spatial diversity and the performance of the LMS method, we evaluated the root mean square (RMS) value of the local phase variation without mechanical perturbations applied to the MMF. In these settings, any deviation from zero can be attributed to noise. The gauge length to compute the local phase variation is set to 10 m. The RMS value is evaluated over all fiber locations considering 1.1 s of data acquired by the DAQ, which correspond to 1.1 s/T rep = 10, 000 time samples. To investigate the impact of the temporal window size N T on the quality of the estimated phase offsets, the 10, 000 time samples were segmented into chunks of N T time samples, with N T in the range 5 ÷ 100 samples. For each chunk, the RVS, DA-RVS, and LMS methods estimate the phase offsets, which are then used to compute the coherent average chunkwise. Consecutive chunks overlap by 1 time sample, so as to retain the continuity of the phase signal. Figure 2(a) shows the RMS value, estimated over the whole fiber length, of the local phase variation versus the time window size N T . The dashed line refers to the single spatial mode performance, which is independent of N T ; whereas, the cross markers and the triangular markers refer to the LMS average of two and three spatial modes, respectively. As expected, the RMS decreases as the number of averaged modes increases, with a significant drop compared to considering each spatial mode separately. The RMS achieved by the LMS method also decreases as N T increases because the detection of the reference sequence with the highest energy becomes more accurate for higher N T values. Figure 2(b) compares the LMS method with the RVS and DA-RVS methods when all three spatial modes are averaged. LMS outperforms the other methods, due to the more accurate selection of the reference sequence and the MSE-optimal estimate of the phase offsets. The RMS reduction saturates as N T increases; therefore, in the

B. Longitudinal Averaging
So far we analyzed how to reduce signal fading by exploiting the spatial mode diversity. In what follows, we generalize this method by exploiting the information diversity carried by the z-samples of the measured Rayleigh backscattered trace. It is not uncommon that ϕ-OTDR systems oversample the backscattered trace, so that the spacing between consecutive z-samples is typically much smaller than the gauge length [12]. Therefore, including a few sampling points in the analysis rather than just one, causes a marginal degradation of the spatial resolution, while giving the opportunity to exploit also a longitudinal diversity. Actually, a number N Z of sampling points close to each other, can be treated as N Z independent measurements of the same event, to which coherent averaging can be applied to further reduce fading. The performance of this average depends on the spectral content of the probe pulse. In general, the closer to Fourier-limited the probe pulse, the more correlated are the z-samples, and hence less information can be obtained by exploiting longitudinal diversity. In this perspective, the best performance is expected with non-Fourier-limited probe pulses [12]. To some extent, this approach is analogous to the frequency-diversity one proposed in [7]; however, the present one is more general, since it does not require any assumptions on the pulse spectral shape. Previous works reported alternative methods for leveraging longitudinal diversity, yet relying on specific fiber configurations, such as wrapping the fiber around a cylindrical cavity structure [14], or restricting the processing to seismic DAS data [15]. In the experimental setting of Fig. 1, we use pulses of 100 ns, which yield a spatial resolution S R of about 10 m, and we sample the backscattered traces at 500 MS/s, which corresponds to a sampling step δz of about 20 cm. Moreover, the bandwidth of the digital band-pass filter employed for I/Q demodulation is set to 30 MHz, so as to retain the first two side lobes of the sinc-shaped backscattered field spectrum. This parameter setting suggests that even samples closer than S R could be sufficiently uncorrelated to make longitudinal diversity effective. To test this claim, for each spatial mode and for each sample at arbitrary position z k , we consider also the samples at z k ± k · δz. This gives a total of M · N Z = 9 sequences (M = 3 modes and N Z = 3 longitudinal samples) over which to compute the coherent average with the above described methods. The parameter k controls the distance between the averaged samples; the larger k, the lower their statistical correlation, yet the worst the effective spatial resolution S eff R = S R + 2 k · δz. The results are shown in Fig. 3(a), where the RMS of the phase noise floor is plotted versus k (top x-axis) and S eff R (bottom x-axis) for each of the three averaging method considered so far; k = 0 corresponds to the case when no longitudinal diversity is considered. Clearly, all methods benefit from longitudinal-diversity averaging for increasing distance between samples, until the RMS almost saturates for k ≈ 15. Note that this value corresponds to a distance of about 3 m, consistently with the 30 MHz bandwidth of the I/Q bandpass filter. Basically, beyond this distance samples are almost completely uncorrelated and the RMS reduction achieved by the averaging is maximum. We may also test if, for a fixed maximum distance k, there is any benefit in considering more than 3 samples in the range z k ± k · δz. The answer is no, as shown in Fig. 3(b) where the RMS for k = 15 is plotted as a function of N Z : the RMS shows negligible variation as a function of N Z , meaning that the samples within the range are too much statistically correlated to provide any advantage in the averaging. To investigate the response of the system to dynamic strain, we inserted a fiber shaker at about 1.1 km from the input. The spatio-temporal maps of the local phase variation in presence a 80-Hz sinusoidal perturbation and the corresponding spectrograms (computed over the whole acquisition time of 1.1 s) are shown in the left and right columns of Fig. 4, respectively. Figure 4(a) shows the results for a single spatial mode; it can be seen that the local phase variation is strongly compromised by interference fading and polarization fading. Figures 4(b) and 4(c) show the results of the LMS average of two spatial modes and three spatial modes, respectively; when all the three spatial modes are averaged, fading errors are significantly reduced. On the spectrogram maps, it can be observed the presence of a faint third-order harmonic at about 240 Hz, which is attributed to the non-ideal response of the fiber shaker. Figure 4(d) shows the results for the coherent average carried out jointly across the three spatial modes and across the z-axis with parameters N Z = 3 and k = 15; remarkably, the phase variation faithfully reconstructs the applied perturbation without fading errors.
V. CONCLUSION In this letter we investigated DAS over MMFs, which is useful in scenarios like oil and gas plants where MMFs are already installed for temperature monitoring. We reported a mathematical model supporting the feasibility of DAS over MMF. Moreover, we proposed a novel coherent averaging method to reduce fading in backscattered traces; we also extended the method to longitudinal diversity, trading spatial resolution for phase accuracy. The proposed coherent averaging method can reduce the noise floor by up to three times compared to state-of-the-art methods.