Analysis of distributed optical fibre acoustic sensors through numerical modelling

A distributed optical fibre acoustic sensor is numerically modelled. To increase the flexibility of the model, the building blocks of the sensing system are modelled separately and later combined to form the numerical model. This approach is adopted to facilitate the evaluation of each of the individual building blocks and their effects on the output of the sensor. The numerical model is used to assess the effect of parameters such as the linewidth of the laser source, the width of the probe pulse, and the frequency and amplitude of perturbation on the response of the sensing system. It is shown that the precision and accuracy of the sensing system are affected by the frequency and amplitude of perturbation as well as the pulse width and linewidth of the probe pulse. © 2017 Optical Society of America OCIS codes: (060.2370) Fiber optics sensors; (000.4430) Numerical approximation and analysis; (290.5870) Scattering, Rayleigh. References and links 1. A. E. Alekseev, V. S. Vdovenko, B. G. Gorshkov, V. T. Potapov, and D. E. Simikin, “A phase-sensitive optical time-domain reflectometer with dual-pulse phase modulated probe signal,” Laser Phys. 24(11), 115106 (2014). 2. C. Wang, C. Wang, Y. Shang, X. Liu, and G. Peng, “Distributed acoustic mapping based on interferometry of phase optical time-domain reflectometry,” Opt. Commun. 346, 172–177 (2015). 3. H. Gabai, Y. Botsev, M. Hahami, and A. Eyal, “Optical frequency domain reflectometry at maximum update rate using I/Q detection,” Opt. Lett. 40(8), 1725–1728 (2015). 4. G. Tu, X. Zhang, Y. Zhang, F. Zhu, L. Xia, and B. Nakarmi, “The development of an φ-OTDR system for quantitative vibration measurement,” IEEE Photon. Technol. Lett. 27(12), 2811–2816 (2015). 5. G. Fang, T. Xu, S. Feng, and F. Li, “Phase-sensitive optical time domain reflectometer based on phase-generated carrier algorithm,” J. Lightw. Technol. 33(13), 2811–2816 (2015). 6. Z. Wang, L. Zhang, S. Wang, N. Xue, F. Peng, M. Fan, W. Sun, X. Qian, J. Rao, and Y. Rao, “Coherent φ-OTDR based on I/Q demodulation and homodyne detection,” Opt. express 24(2), 853–858 (2016). 7. G Yang, X. Fan, S. Wang, B. Wang, Q. Liu, and Z. He, “Long-range distributed vibration sensing based on phase extraction from phase-sensitive OTDR,” IEEE Photon. J. 8(3), 6802412 (2016). 8. J. Pastor-Graells, H. F. Martins, A. Garcia-Ruiz, S. Martin-Lopez, and M. Gonzalez-Herraez, “Single-shot distributed temperature and strain tracking using direct detection phase-sensitive OTDR with chirped pulses,” Opt. Express 24(12), 13121–13133 (2016). 9. A. Masoudi and T. P. Newson, “High spatial resolution distributed optical fibre dynamic strain sensor with enhanced frequency and strain resolution,” Opt. Lett. 42(2), 290–293 (2017). 10. X. He, S. Xie, F. Liu, S. Cao, L. Gu, X. Zheng, and M. Zhang, “Multi-event waveform-retrieved distributed optical fiber acoustic sensor using dual-pulse heterodyne phase-sensitive OTDR,” Opt. Lett. 42(3), 442–445 (2017). 11. S. Liehr, Y. S. Muanenda, S. Munzenberger, and K. Krebber, “Relative change measurement of physical quantities using dual-wavelength coherent OTDR,” Opt. Express 25(2), 720-729 (2017). 12. Y. Peled, A. Motil, L. Yaron, and M. Tur, “Slope-assisted fast distributed sensing in optical fibers with arbitrary Brillouin profile,” Opt. Express 19(21), 19845–19854 (2011). 13. J. Urricelqui, A. Zornoza, M. Sagues, and A. Loayssa, “Dynamic BOTDA measurements based on Brillouin phase-shift and RF demodulation,” Opt. Express, 20(24), 26942–26949 (2012). 14. A. Masoudi, M. Belal, and T. P. Newson, “Distributed dynamic large strain optical fiber sensor based on the detection of spontaneous Brillouin scattering,” Opt. Lett. 38(17), 3312–3315 (2013). 15. Y.Mizuno, N. Hayashi, H. Fukuda, K.Y. Song, andK.Nakamura, “Ultrahigh-speed distributed Brillouin reflectometry,” Light Sci. Appl. 5, e16184 (2016). 16. A. Masoudi and T. P. Newson, “Contributed Review: Distributed optical fibre dynamic strain sensing,” Rev. Sci. Instrum. 87(1), 011501 (2016). 17. L.B. Liokumovich, N.A. Ushakov, O.I. Kotov, M.A. Bisyarin, and A.H. Hartog, “Fundamentals of Optical Fiber Sensing Schemes Based on Coherent Optical Time Domain Reflectometry: Signal Model Under Static Fiber Conditions,” J. Lightw. Technol. 33(17), 3660–3671 (2015). Vol. 25, No. 25 | 11 Dec 2017 | OPTICS EXPRESS 32021 #302831 Journal © 2017 https://doi.org/10.1364/OE.25.032021 Received 19 Jul 2017; revised 10 Nov 2017; accepted 17 Nov 2017; published 7 Dec 2017 18. A. H. Hartog and K. Kader, “Distributed fiber optic sensor system with improved linearity,” US Patent 2012/0 067118A1 (2012). 19. A. Masoudi, M. Belal, and T. P. Newson, “Distributed optical fibre audible frequency sensor,” Proc. SPIE 9157, 91573T (2014). 20. A. Masoudi, M. Belal, and T. P. Newson, “A distributed optical fibre dynamic strain sensor based on phase-OTDR,” Meas. Sci. Technol. 24(8), 085204 (2013). 21. S. K. Sheem, T. G. Giallorenzi, and K. Koo, “Optical techniques to solve the signal fading problem in fiber interferometers,” Appl. Opt. 21(4), 689–693 (1982). 22. M. D. Mermelstein, R. Posey, G. A. Johnson, and S. T. Vohra, “Rayleigh scattering optical frequency correlation in a single-mode optical fiber,” Opt. Lett. 26(2), 58–60 (2000). 23. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman & Hall, 1995). 24. G. P. Agrawal, Nonlinear Fiber Optics, 4th ed. (Academic Press, 2007). 25. C. D. Butter and G. B. Hocker, “Fiber optics strain gauge,” Appl. Opt. 17(18), 2867–2869 (1978). 26. G. B. Hocker, “Fiber-optic sensing of pressure and temperature,” Appl. Opt. 18(9), 1445–1448 (1979). 27. P. Healey, “Statistics of Rayleigh backscatter from a single-mode fiber,” Electron. Lett. 21(6), 226–228 (1985). 28. K. Shimizu, T. Horiguchi, and Y. Koyamada, “Characteristics and reduction of coherent fading noise in Rayleigh backscattering measurement for optical fibers and components,” J. Lightw. Technol. 10(7), 982–987 (1992). 29. L. Shi, T. Zhu, Q. He, and S. Huang, “Effect of laser linewidth on phase-OTDR based distributed vibration sensing regime,” Proc. SPIE 9157, 91576H (2014). 30. J. C. Juarez, E. W. Maier, K. N. Choi, and H. F. Taylor, “Distributed fiber-optic intrusion sensor system,” J. Lightw. Technol. 23(6), 2081–2087 (2005).


Introduction
The volume of research on distributed optical fibre dynamic strain sensing has increased substantially during the past five years.The main driving force behind this expansion was the rising demand for this type of sensors in areas such as oil and gas industry, geophysical sciences, and civil engineering.Distributed optical fibre acoustic sensors, also known as distributed vibration sensors (DVS), owe their rising popularity to their capability of mapping vibrations along a long, narrow, and usually inaccessible subterranean spaces such as boreholes, often under harsh environmental conditions.The sensing principle of DVS allows the interrogation unit of such systems to be kept at a safe distance while the sensing fibre connected to the interrogation unit can be encapsulated in several protective layers to withstand harsh conditions.
Up until now, the main focus of the research studies has been on exploring different optical mapping techniques using Rayleigh [1][2][3][4][5][6][7][8][9][10][11] or Brillouin scattering [12][13][14][15].A detailed review of these measurement techniques and their key parameters is documented in [16].To improve on the recent studies, numerical analysis of the sensing techniques investigated so far can be of significant importance.In 2015, Liokumovich et al. have published an article on numerical modelling of Rayleigh backscattering process with an emphasis on phase-sensitive optical time domain reflectometry (ϕ-OTDR) sensing technique [17].In this article, the theory of Rayleigh backscattering was established followed by mathematical modelling of a ϕ-OTDR setup introduced by Hartog et al. [18].The study focused on one of the three ϕ-OTDR sensing techniques addressed in [16] i.e. the sensing technique which directly measures the phase of the backscattered light by converting the backscattered traces from the optical domain to the electrical domain.The study, however, was limited to the analysis of the sensing system under static conditions and did not incorporate the effect of dynamic perturbations.
The aim of this study is to present an alternative approach to numerical analysis of phasesensitive DVS.In this approach, the building blocks of the sensing system such as the sensing fibre, the probe pulse, and the detection unit are modelled separately and later combined to model the behaviour of the sensing system.Such an approach offers two advantages: (i) It allows each element in the system to be modified independently hence facilitating the evaluation of the effect of each element on the output of the system and, (ii) it allows noises and deficiencies inherent to each module to be added separately which, in turn, facilitates refining the numerical model.In this paper, the framework of the numerical analysis is established by modelling the building blocks of the sensing system in their ideal states.In order to avoid complexity and simplify the assessment of the core modules and their interactions, the effect of noises and other deficiencies are not included in this study.Figure 1 shows the building blocks used to model the behaviour of the sensor.The content of the blocks and their role in the numerical analysis are described in the next section.
The new modelling approach is used to assess the behavior of the ϕ-OTDR sensing technique introduced by Masoudi et al. [19,20].This technique uses the phase of the Rayleigh backscattered light to measure dynamic perturbation, but instead of direct phase measurement in the electrical domain, it measures the relative phase of the backscattered light in the optical domain.In this technique, the backscattered light from the sensing fibre is launched into an imbalanced Mach-Zehnder interferometer (IMZI) with two uneven arms to from two similar traces with a temporal shift.By mixing the two traces at the output of the interferometer, the relative phase difference between the adjacent sections of the sensing fibre can be measured.The schematic of the setup used to demonstrate this sensing technique is presented in figure 2. The role of the symmetric 3 × 3 coupler in the setup is to form three interference signals with a relative phase difference of 2π/3 in order to avoid interferometer signal fading [21].In order to avoid coherence fading, a probe pulse with a relatively broad linewidth was used.Furthermore, the sensing system showed no polarization-dependent signal fading [22], a phenomenon which was attributed to the fact that the average beat length of standard single-mode fibre (SMF) was significantly longer than the gauge length of the sensor.
The rest of this paper is structured as follows: section 2 presents the principle of the numerical modelling.In this section, the building blocks of the model and their parameters are described and

Linewidth
Pulse width Spacing it is shown how these blocks interact to simulate the operation of the sensor.Section 3 discusses the simulation procedure and the parameters used to model the sensing system.In section 4, the effect of the parameters such as the linewidth of the laser source and the pulse width on the response of the sensing system are simulated and the results presented.Section 5 is dedicated to the analysis of the simulation results to determine how each parameter affects the output of the sensor and how these parameters can be modified to optimize the performance of the sensing system.The final section summarizes the study and establishes the scope for future work.

Numerical modelling principles
In order to numerically analyse an OTDR system, the backscattered light from the sensing fibre of that system needs to be modelled.In ϕ-OTDR systems where the phase of the Rayleigh backscattering plays an important role, the electric field of the backscattered light has to be modelled in order to conserve the phase information.The backscattered electric field of a ϕ-OTDR system can be modelled by analysing the interaction between the probe pulse and the frozen-in inhomogeneities in the fibre which give rise to Rayleigh scattering.The flowchart of figure 1 shows the elements used to model the ϕ-OTDR sensing system of figure 2. These elements include: • FibGen: The module used to model the sensing fibre; • Probe Pulse: The module used to specify the probe pulse and its parameters; • BEF: The module used to model the backscattered electric field as a result of the interaction between the probe pulse and the sensing fibre; • IMZI: The module used to simulate the transfer function of the imbalanced MZI; • PD: The module used to model the response of the photodetectors to the incident electric field; • Strain: The module used to simulate the effect of dynamic vibrations on the Rayleigh scattering in the sensing fibre.
The remainder of this section describes these modules separately and delineates the way they interact to model the sensing system.Sensing Fibre: As mentioned earlier, Rayleigh scattering occurs as a result of the interaction between light and inhomogeneities in the sensing fibre.The role of the FibGen module is to The first column of the matrix determines the distribution of inhomogeneities as a function of distance while the second column specifies the size of the inhomogeneities.Here, the size of an inhomogeneity determines the strength of the scattering from an interval.To generate this matrix, the following principles were considered: 1.The size of the scatterers in an optical fibre are much smaller than the wavelength of the probe light [23].Due to the limited processing power of computers, modelling the individual inhomogeneities is impractical.The alternative approach adopted is to divide the sensing fibre into n individual intervals and to represent all the scatterers within each interval by a single scattering centre.In this approach, the length of the intervals is kept much smaller than the width of the interrogating pulse.
To model the random distribution of inhomogeneities in the sensing fibre (figure 3), random locations were assigned to the scattering centres using the following formula: where L k is the location of the k-th scattering centre and D is the length of the scattering interval.This randomness determines whether the backscattered waves from each interval are added constructively or destructively which, consequently, forms the peaks and troughs of the backscattered coherent Rayleigh noise (CRN).
2. The size of the inhomogeneities, s k , are independent random quantities with equal distributions that define the magnitude of the scattered light from each scattering centre [17].Unlike the location of the scattering centres that changes as a result of an external perturbation on the fibre, the size of the scatterers remain constant when the fibre undergoes strain.This statement is valid, but only under the assumption that fibre elongation is much smaller than the length of the scattering interval, D. This assumption guarantees that as the fibre is perturbed: (i) the relative spacing of the true scattering elements within interval D remains reasonably fixed and, (ii) the number of the true scattering elements that enter or leave the interval D does not significantly affects the contribution of the remaining scattering elements.
Probe Pulse: To fully characterize the probe pulse, both the spatial and spectral distribution of the pulse needs to be specified as shown in figure 4. For each data point on the probe pulse, ... ...  an amplitude and a phase was assigned.Therefore, the probe pulse can be fully specified by an N × m × 2 matrix called PP where N determines the frequency components forming the probe pulse, m determines the number of elements used to represent the spatial distribution of the probe pulse, and PP(q, k, 1) and PP(q, k, 2) determine the amplitude a(q) k = a(λ q , d k ) and phase ϕ(q) k = ϕ(λ q , d k ) of the data points, respectively, where q ∈ (1, N) and k ∈ (1, m).

k-th interval
Defining the probe pulse in this way provides the flexibility to model the system with different pulse parameters.The linewidth of the probe pulse was modelled by N individual single frequency components with a fixed frequency separation of ∆ν/N.The value of ∆ν should satisfy two conditions: 1.The linewidth of a light pulse has to be greater than the Fourier transform limit of the pulse given by ∆ν × τ p = 0.315 for sech-squared pulse, or where τ p is the probe pulse duration.
2. The value of the frequency separation, ∆ν/N, has to be determined according to the width of the probe pulse.The value of the frequency separation should be sufficiently small to guarantee that the backscattered light from scatterers within the pulse duration interfere.For a light source with a linewidth of ∆ν/N, the coherence length is give by [24] where n is the refractive index.The value of the frequency separation should be small enough to satisfy the inequality L pw L C where L pw is the pulse width.
For each spectral component of the probe pulse, the variation in the phase of the probe pulse depends on the wavelength of that spectral component and the spatial distribution of the pulse.For the spectral component λ q , for instance, a random phase was assigned to the first spatial element ϕ(λ q , d 1 ) and the phase value of the remaining elements was determined through the following equation where ∆d is the spatial separation between two adjacent elements.
Backscattered Electric Field: In an OTDR system, the interaction between the probe pulse and the inhomogeneities in the fibre gives rise to the backscattered light.This interaction was modelled using the fibre (F i b) and probe pulse (PP) matrices introduced earlier.Since the response of some of the modules in the numerical model are wavelength dependent, a separate backscattered electric field (BEF) was modelled for each spectral component in the probe pulse.In other words, for N wavelengths in the probe pulse, λ 1 ∼ λ N , N individual BEF was formed.These BEFs were later combined in the photodetector to model the Rayleigh backscattered trace.
The BEF from any point on the sensing fibre is the superposition of m individual backscattered electric fields corresponding to m spatial elements of the probe pulse.For the spectral component λ = λ q , the BEF from position p on the fibre is given by: where j is the imaginary number and In this equation, s u+p and L u+p denote the size and location of the scattering points in the sensing fibre (Eq.( 1)), and β q is the wavenumber for λ = λ q .This equation shows that the amplitude of the BEF is determined by the product of the size of the scattering points s u+p and the intensity of the probe pulse a(q) u .The phase of the BEF, on the other hand, is determined by the location of the scattering points L u+p , the wavenumber β q , and the phase of the probe pulse ϕ(q) u .Using trigonometric identities, Eq. ( 5) can be simplified to: where and Using Eqs. ( 7)-( 9), the backscattered electric field can be defined by a N × r × 2 matrix EF where N determines the spectral components present in the BEF, r determines the number of the data points in the BEF trace, and EF(p, 1, λ q ) and EF(p, 2, λ q ) determine the amplitude A(λ q , L p ) and phase φ(λ q , L p ) of the BEF, respectively.

Imbalanced MZI:
The role of this module is to model the effect of the imbalanced MZI (IMZI) on the backscattered electric field matrix EF.According to the experimental setup shown in figure 2, the BEF from the sensing fibre is fed into the IMZI via a symmetrical 2 × 2 coupler.After propagating through the two unequal arms with a path different of ∆L, the electric fields are mixed in a symmetrical 3 × 3 coupler.The emergent electric fields from the three arms of the 3 × 3 coupler, O 1 ∼ O 3 , can be represented by three N × r × 2 matrices.Figure 5 illustrates the mathematical procedure used to convert the backscattered electric field EF to the outputs of the IMZI.The amplitudes of the electric fields at points A and B are equal to the amplitude of the EF divided by √ 2. The arm of the coupler connected to A is assumed to be the through port of the coupler and, as a result, the phase of the electric field at point A remains unchanged.The other arm of the coupler connected to B is the coupled port and its phase can be calculated by adding π /2 to the phase of EF.
The delay line is modelled by adding zeros to the beginning of the electric field matrix.The introduction of zero rows to the beginning or the electric field matrix is equivalent to temporally shifting the backscattered field forward in time.The number of zero rows added to the matrix depends on the path imbalance, ∆L, and the spatial separation between data points of the BEF.Therefore, the matrix that represents the electric field at point Ā can be formed by adding zero rows to the beginning of the electric field matrix at point A.
The electric fields at the output of the IMZI can be calculated using the electric field matrices at Ā and B. Each output arm of the 3 × 3 coupler receives a third of the energy of each input arm.Therefore, the amplitude of the electric fields at Ā and B should be multiplied by 1 / √ 3 before being combined in the outputs.The phase relationship between an input and an output port of a 3 × 3 coupler depends on whether the input and output arms form a through port or a coupled port.For instance, the input port Ā and the output port O 1 form a through port which means there are no phase difference between the phase of the electric field at Ā and O 1 .On the other hand, the input port B and the output port O 1 form a coupled port which means there is 2π /3 phase difference between the electric fields at the two ports.Using these principles, the electric fields at the output of the 3 × 3 coupler is given by: 3 ) 3 ) In this equation, E Ā and E B are both sum of sinusoidal waves and can be written as: where α i and β i represent the phase of the through port and coupled port at the output of the 2 × 2 coupler, respectively.Substituting E Ā and E B from Eq. ( 11) into Eq.( 10), gives: where Equation (13) shows that the intensity of the 3 × 3 coupler outputs are 2π/3 out of phase.The phase of the 3 × 3 coupler outputs (Φ 1 (i), Φ 2 (i), and Φ 3 (i) in Eq.( 12)) are a function of the electric fields E Ā and E B .But since they do not affect the detectors, they can be ignored.

Photodetectors:
The photodetector module converts the backscattered electric field to intensity.The intensity of the BEF at the photodetectors can be calculated using the superposition principle: the net electric field at any given time is equal to the sum of all individual fields.Using Eq. ( 12), the intensity of the backscattered electric field of, for instance, the first output arm is where T is determined by the bandwidth of the detector.Expanding this equation gives where ∆ω = (ω i − ω k ) and ∆Φ = (Φ 1 (i) − Φ 1 (k)).In this expansion, the high-frequency terms, sin[( , are ignored since their frequency values are higher than the bandwidth of the detector. Equation ( 15) consists of two terms: the DC component (the first term on the right hand-side of the equation) and the frequency-dependent term (the second term on the right hand-side of the equation).As the BEFs from two separate segments of the optical fibre mix, the backscattered fields with equal frequencies interfere to form the DC component.According to Eq. ( 13), the value of the DC component depends on the difference in the phase of the backscattered electric field.The frequency-dependent term, on the other hand, appears as the by-product of the BEFs mix.The value of the frequency-dependent term depends on both the beat frequencies and the phase differences.Hence, the frequency-dependent term can be viewed as noise, added to the DC component.
To simplify Eq. ( 15), the bandwidth of the detector is initially assumed broad enough to detect the backscattered light from individual scattering centres.Under this assumption, the values of the intensity and phase of the BEF would be constant within the integration interval 0 ∼ T. As a result, the frequency-dependent term in Eq. ( 15) can be modified by reversing the order of integration and summation as follows: The integration period 0 ∼ T is determined by where D is the scattering interval.Since the intensities and phases of the BEFs are constant within the limits of the integral, Eq. ( 16) becomes: This equation was used to convert the BEF to the backscattered intensity trace.The effect of the bandwidth of the detector was modelled using a moving average i.e. the output of a detector with a bandwidth BW is obtained by calculating the moving average of 1 /BW.T data points on the backscattered intensity trace where T is given in Eq. (17).
Perturbation: The role of this module is to simulate perturbations on a single or multiple points along the sensing fibre.The effect of a perturbation on the fibre was assumed to be a periodic elongation of fibre along its main axis.The impact of perturbation on any given section of the fibre was modelled by rearranging the position of the scattering centres with regards to the frequency and amplitude of the perturbation at that section.It is worth mentioning that the rearrangement process does not affect the size of the scattering centres.
The phase of the BEF is a function of both the dimensional and refractive index changes in the fibre.The relationship between the changes in the phase of the BEF, ∆Φ, and the induced strain on the fibre, , is given by [25]: where l is the length of the fibre under strain, β is the propagation constant, n is the refractive index, µ is the Poisson's ratio, and p 11 and p 12 are strain-optic coefficients.Replacing the values of refractive index (n = 1.456),Poisson's ratio (µ = 0.17), and strain-optic coefficients (p 11 = 0.121, p 12 0.27) [26] in Eq. (19) gives Equation (20) indicates that the net change in phase is proportional to 78% of the phase change due to the geometric elongation of the fibre.Therefore, to accurately simulate the effect of induced strain on the phase of the BEF, the effect of strain-optic coefficient (i.e.changes in the RI of a medium as a result of the changes in the strain of the medium) was incorporated in the model.In order to assess the effect of dynamic perturbation on the numerical model, a dynamic perturbation loop was added to the model as shown in the flowchart of figure 1.The role of this loop is to control the periodic changes in the position of the inhomogeneities based on the amplitude and frequency of perturbations.

Numerical modelling procedure
Prior to any analysis, the reliability of the numerical model was assessed through two tests.The first test involved examining the CRN level of the backscattered traces generated by the model as a function of probe pulse linewidth and comparing the results with the mathematical model of Rayleigh backscattering from a single mode fibre [27,28].In this test, the linewidth of 1m probe pulses were increased from 1pm to 100pm in eight steps and the CRN level of the backscattered traces were compared.The second test involved evaluating the changes in the CRN pattern at the output of the IMZI as a result of external perturbation on the fibre and comparing the outcome of the simulation with the experimental results [29].For this test, probe pulses with two different linewidths (3pm and 30pm) were used to generate the backscattered light.The repetition rate of the probe pulse was set to 10µs to trace the changes in the backscattered light as a result of 1kHz sinusoidal strain.
Following the reliability tests, the numerical model was used to evaluate the dynamic behaviour of the model.The main focus of the numerical evaluation was to gain a better understanding of how different parameters of the system affect the sensor output and how they can be modified to improve the accuracy of the sensor.The evaluation included assessing the accuracy of the sensor's response as a function of parameters such as the frequency and amplitude of perturbations, the linewidth of the light source, and the width of the probe pulse.
For the first trial, the sensing system was modelled for a simple scenario in which dynamic strains were imposed on two separate sections of the sensing fibre.The first section was subjected to a 1kHz sinusoidal strain with an amplitude of 2.5µε.The strain was imposed on a meter long section of the fibre from the 7 th meter to the 8 th meter.The second section was subjected to a 2kHz sinusoidal strain with an amplitude of 1.5µε.The strain was imposed on a 2m section of the fibre between the 19 th to the 21 st meter.For this test, the pulse width and the linewidth of the probe pulse was set to 0.5m and 5pm, respectively.
The effect of the frequency and amplitude of perturbations on the accuracy of the sensor was evaluated by assigning a fixed value to all of the parameters in the sensing system while stepping the frequency and amplitude of perturbations through a range of values.To examine the effect of strain, nine different strain levels from 0.1µε to 2.0µε were applied at a fixed frequency of 1kHz.The effect of frequency was studied by stepping the frequency of perturbation from 250Hz to 2000Hz at a fixed strain level of 0.5µε.For both sets of analysis, the pulse width and linewidth of the probe pulse were set to 0.5m and 3pm, respectively.
To evaluate the effect of probe pulse linewidth on the accuracy of the sensor, the behaviour of the system was modelled for a range of linewidths between 0.1pm to 100pm.The same parameters were used for all five linewidths including 1kHz perturbation with strain rate of 0.5µε.For this set of analysis, the pulse width was set to 3m.The path-imbalance of the IMZI was set to 6m to incorporate 3m pulse width.
To assess the effect of pulse width on the accuracy, the sensing system was modelled by stepping the pulse widths from 50cm to 2m in 25cm steps.The frequency and strain of perturbation were set to 1kHz and 1µε, respectively, for all of the simulations while the path-imbalance of the IMZI was fixed at 4m. 3pm linewidth was used for the probe pulse in this analysis.
In order to maintain consistency, fixed values were assigned to a number of parameters of the numerical model such as the detector bandwidth (BW Det = 300M Hz), sampling rate (F Sa = 5GSa/sec), the length of the scattering interval (D = 20mm), and the frequency separation between the wavelength components (∆ν/N = 0.02pm ≡ 2.5M Hz).D = 20mm guaranteed that the scattering interval is much shorter than 50cm pulse width used to analyse the model.It also conforms with the assumption that the size of the scattering elements does not change as a result of perturbation given the maximum elongation assessed in this study was 5µm, equivalent to 0.025% of the length of the interval.

Simulation results
Figure 6 depicts the outcome of the first reliability test.Figure 6(a) shows the backscattered CRN pattern for three different linewidths, namely, 1pm, 10pm, and 100pm.The dashed line in this diagram determines the mean value of the backscattered trace i.e. the backscattered trace for a broadband source with an infinitely broad linewidth.Figure 6(b) represents the probability density function (PDF) of the intensity of the backscattered CRN for the three linewidths.The solid traces in this figure display the simulation results while the dashed traces determine the PDF of the backscattered light obtained through theoretical analysis [27].The dot-dashed line in figure 6(b) indicates the normal line, the line which represents the intensity of the backscattered light with an infinitely short coherence length.Figure 6(c) shows the changes in the CRN level of the backscattered trace as a function of the square root of the probe pulse coherence-length, L coh .
Figure 7 shows the result of the second reliability test of the simulation.In this test, the changes in the pattern of the backscattered CRN for a fibre under dynamic strain was assessed using probe pulses of two different linewidths.Each graph comprises of several backscattered traces, each representing the CRN pattern at a certain time.With the repetition rate of the probe pulse set to 100kHz, each trace represents the backscattered CRN for a 10µs window.The plot on the left shows the simulation result for a narrow linewidth source (3pm) while the plot on the right shows the simulation result for a source with a broader linewidth (30pm).
Figure 8 shows the outcome of the numerically modelled sensing system interrogated by a probe pulse with a linewidth and pulse width of 5pm and 50cm, respectively, while two regions of the sensing fibre were sinusoidally modulated.The first section was subjected to a 1kHz sinusoidal strain with an amplitude of 2.5µε while the second section of the fibre undergone a 2kHz sinusolidal strain with an amplitude of 1.5µε.Figure 8(a) shows the backscattered CRN of the sensing fibre after the IMZI.The regions where dynamic strains were applied are distinguished with two blue boxes.The first box identifies a 2m section between the 7 th and the 9 th meter while the second box identifies a 3m section between the 19 th and the 22 nd meter.Figure 9 shows the behaviour of the numerical model to 1kHz sinusoidal strain at various strain levels.Figure 9(a) exhibits the output of the numerical model as a function of induced strain.The model was used to assess the system response to different strain levels between 100nε and 2µε.The data for each strain level was obtained by averaging the results from twenty separate realisations of the simulation each with a different random collection of scattering centres.The error bar assigned to each data point represents the standard deviation of the simulation results at each strain level.The standard deviation at each point indicates the variation in the outcome of different realisations of the numerical model.Figure 9(b) illustrates the standard deviation of the data points shown in figure 9(a) as a function of induced strain.
Figure 10 shows the behaviour of the sensing system to a fixed strain level at different frequencies.Figure 10(a) exhibits the output of the numerical model to 0.5µε sinusoidal strains at different frequencies ranging from 250Hz to 2000Hz.Like the previous assessment, the data for each frequency value was obtained by averaging the results from twenty separate realisations of the system each with a different random collection of scattering centres.The error bar assigned to each data point represents the standard deviation of those individual realizations.Figure 10(b) illustrates the relationship between the frequency of the induced strain and the standard deviations of the data points in figure 10(a).represents the regression line fitted to the data points.The horizontal axis of this diagram is plotted in logarithmic scale to clarify the relationship between the linewidth of the probe pulse and the standard deviation of the results.Finally, figure 11(b) demonstrates the relationship between the standard deviation of the simulation result as a function of the width of the probe pulse.The system model was used to assess the effect of different pulse widths from 50cm to 2m on the output of the sensor.The standard deviation for each pulse width was obtained by measuring the root mean square of twenty different realization of the model, each with a different random collection of scattering centres.The dashed line in the graph represents the regression line fitted to the data points representing the simulation results.

Discussion
The preliminary test results shown in figures 6 and 7 indicate that the behaviour of the numerical model conforms with the theoretical and experimental studies [27][28][29][30].Figure 6(a) shows an inverse relationship between CRN level and the square root of the linewidth of the probe pulse, (a) Relationship between the standard deviation of the strain level simulated using the numerical model and the linewidth of the probe pulse.For this diagram, the sensing system was modelled for a range of linewidths from 0.1pm to 100pm for a fixed pulse width of 3m.(b) Relationship between the standard deviation of the strain level simulated using the numerical model and the width of the probe pulse.For this diagram, the sensing system was modelled by stepping the pulse widths from 50cm to 2m in 25cm steps while the linewidth of the probe pulse was fixed to 3pm.
∆λ.Here, the CRN level is defined as the RMS of the coherent Rayleigh noise where the mean value used in the calculation is determined by the intensity of the backscattered trace for a broadband source with an infinitely broad linewidth.It can be seen that a probe pulse with 1pm linewidth has a backscattered trace with significantly larger deviation from the normal line compared with the backscattered trace of a probe pulse with 100pm linewidth.The relationship between the coherence length and linewidth of a light source can be calculated by substituting the frequency-domain linewidth ∆ν in Eq. ( 4) with its equivalent spectral-domain linewidth ∆λ: where λ is the wavelength of the light source.A probe pulse with a long coherence length results in a highly coherent backscattered light.For a highly coherent backscattered light, the interaction between the backscattered wavelets from individual scattering centres is more likely to result in constructive or destructive interference, hence forming the peaks and troughs of the green trace shown in figure 6(a).A probe pulse with a short coherence length, on the other hand, results in a backscattered trace which barely deviates from the normal line.That is due to the fact that incoherent waves do not interact with one another to form peaks and troughs.Instead, such waves are added on an intensity basis where the total intensity is equal to the sum of the intensities of uncorrelated wavelengths that form the backscattered light.The solid curves in figure 6(b) which present the PDF of the backscattered intensities have a good agreement with the theoretical analysis of Rayleigh scattering in a one-dimensional optical fibre [27,28] which are shown by the dashed curves.According to the mathematical analysis, a highly coherent backscattered trace has a negative exponential PDF.As the coherence length decreases, the PDF of the intensity distribution starts to converge to the normal line, the line which represents the intensity of the backscattered light with an infinitely short coherence length.This behaviour can be observed in the PDF distribution of figure 6(b).The green curve that represents the PDF of the backscattered light for a source with a long coherence length (∆λ = 1pm ) has a relatively even distribution which spreads well beyond the normal line (dot-dash line in the figure).On the other hand, the blue curve that represents the PDF of the backscattered light with a shorter coherence length (∆λ = 100pm ) has a distribution with a mean value close to the normal line and a relatively small variance.
The discrepancy between the theoretical analysis (the dashed curves) and the numerical results (the solid curves), especially for the data at left hand-side of diagram 6(b), is due to the approximations made in the theoretical analysis.In the theoretical analysis [27], the probe pulse was divided into M subintervals such that the backscattered light at different intervals are uncorrelated.In reality, however, any two adjacent intervals have a degree of coherence depending on their widths, a factor which is incorporated in the numerical analysis.Therefore, the numerical analysis results predict the behaviour of the system more accurately.
For a probe pulse with a coherence length much shorter than the pulse width (L coh L pw ), CRN level of the backscattered light, σ C R N , is proportional to the square root of the coherence length.That is because a light pulse with a spatial width to coherence length ratio of N (L pw /L coh = N) acts as N individual incoherent waves, reducing the CRN level by a factor of 1/

√
N: This proportionality is shown by the diagonal dashed-line in figure 6(c).For a light pulse with a coherence length much longer than the pulse width, on the other hand, the interaction between the light pulse and the scattering elements in the fibre results in a highly coherent backscattered light.For such a pulse, the CRN level of the backscattered light reaches an asymptotic level, independent of the coherence length.This level is marked by a horizontal dashed-line in figure 6(c).The data point on the right hand-side of this figure represents CRN level for a probe pulse between the two extremes.The trend line in this figure shows the evolution of the CRN level from a large L pw /L coh ratio probe pulse (where the pulse width is much longer than the coherence length of the probe pulse) to a small L pw /L coh ratio probe pulse (where the coherence length of the probe pulse is much longer than the pulse width).
The test result depicted in figure 7 shows that the intensity fluctuation of the backscattered light for a narrow linewidth probe pulse is much larger than that of a broad linewidth probe pulse, a behaviour that has been experimentally observed [29].With a narrow linewidth pulse, a small perturbation along the sensing fibre results in a significant variation in the intensity of the backscattered light.This interdependency between the linewidth and intensity fluctuation is why distributed optical fibre intrusion sensors use narrow linewidth or single-frequency laser source to increase their sensitivity [30].
The results of the two reliability tests showed good agreement with the previously reported theoretical and experimental studies hence validating the operation of the numerical model and providing a certain level of confidence in its ability to model OTDR systems.
The results of the preliminary system test, shown in figure 8(a), shows that the backscattered traces deviates from the normal CRN pattern in the two regions of the fibre where dynamic strains are applied.The backscattered trace diagram of figure 8(a) shows that while a 1m and a 2m segments of the fibre underwent dynamic perturbations, a 2m and a 3m sections of the backscattered trace were affected.This is as a result of the interaction between the probe pulse with a length L pw and the section of fibre under perturbation.The backscattered trace starts to deviates from the normal CRN trace as the front-end of the probe pulse enters the perturbed region and continues until the back-end exits that region.
The 3D diagrams of figure 8 exhibits a response similar to what has been previously observed using an experimental setup [20].Figure 9(a) shows a linear relationship between the induced strain and the output of the numerical model.The graph shows that the standard deviation of the output increases with the level of the induced strain.The relationship between the strain level and the standard deviation is shown in figure 9(b).The behaviour of the standard deviation as a function of strain level can be explained with the help of figure 12.This figure shows the distribution of inhomogeneities in two fixed sections of the fibre before and after elongation.In the top diagram which shows the fibre before elongation, the two sections are meters apart.The width of each section, W , represents the length of the fibre covered by the pulse width.The phase of the backscattered light from each section is determined by the distribution of inhomogeneities within that section.For an unperturbed fibre, the phase-difference between the backscattered signals from the two sections is where ϕ 1 and ϕ 2 are two random phases.Any disturbance along the fibre results in redistribution of inhomogeneities within that section.The diagram at the bottom of figure 12 shows a scenario where the fibre is longitudinally stretched by ∆ meters.This change in the length of the fibre results in redistribution of inhomogeneities which leads to variation in the phase-difference between the backscattered light from the two sections.Following the redistribution, the value the phase-difference, ∆ϕ, is given by: where ϕ 1 and ϕ 2 are the phases of the backscattered light from the same stretch of fibre prior elongation.The level of change in the phase-difference can be obtained by subtracting Eq. ( 24) from Eq. ( 23): The first term in Eq. ( 25) represents the linear changes in the phase-difference while the other two terms correspond to non-linear changes.Equation (25) shows that linear changes in the phase-difference is directly proportional to the strain level.The non-linear changes in the phasedifference occur due to redistribution of scatterers within the two regions from which the phase of the backscattered light are measured.From figure 12, it can be seen that as the fibre stretches, the relative distance between the inhomogeneities within each section changes (d → d + ∆d).This relative shift results in disproportionate phase variation between ϕ 2 and ϕ 2 as well as ϕ 1 and ϕ 1 .For a small strain level, the relative shift between scatterers is negligible and, therefore, the deviation in ∆φ is small.The deviation starts to increase as the strain level increases.This correlation can be observed in figure 9(b).Figure 10(a) shows that the output of the differentiate and cross-multiplying (DXM) demodulator is frequency independent.The graph shows that the error bar of each data point covers the expected value of the response.Figure 10(b) shows that low frequency perturbations have higher standard deviation compared with higher frequency perturbations.For any periodic signal, the precision of the measurement is proportional to the number of cycles sampled.For a high frequency perturbation, the signal processing procedure takes advantage of the data acquired over many cycles to evaluate the strain level.For low frequency perturbation, fewer cycles can be sampled during the same period.Therefore, the sensing system demonstrate a better standard deviation for higher frequency perturbations.
The analysis of figure 11(a) indicates that the accuracy of the output of the numerical model is linearly proportional to the log of the probe pulse linewidth.This relationship can be explained with the help of Eq. (15).As mentioned earlier, this equation consists of two terms, the DC component and the frequency-dependent term.From this equation, it can be seen that the contribution of the DC component to the intensity of the backscattered light is proportional to the number of the frequency components present in the probe pulse.Therefore, using a probe pulse with a broader linewidth increases the influence of the DC component (i.e. the desirable component) on the output of the detector.In contrast, the effect of the frequency-dependent term on the output of the detector reduces as the linewidth of the probe pulse increases.The frequency-dependent term in Eq. ( 15) consists of sum of all the beat frequencies present in the probe pulse.This term is a quasi-random number with zero mean.The standard deviation of this term is inversely proportional to the linewidth of the probe pulse which means that as the linewidth of the probe pulse increases, it is less likely that the value of the frequency-dependent term deviates significantly from zero.
Finally, figure 11(b) shows that the standard deviations of the strain levels obtained via numerical modelling have a linear relationship with the duration of the probe pulse.The linear relation between the pulse width and the standard deviation can be explained with the help of figure 12.As mentioned earlier, the non-linear changes in the phase of the backscattered light occurs due to redistribution of scatterers in the fibre.The level of the non-linear changes depends on the extent the inhomogeneities within the width of the probe pulse are displaced relative to one another as the fibre stretches.For a long probe pulse that covers a long section of the fibre, the separation between the scatterers is large, especially for the scatterers that reside at the two ends of that section.Therefore, even a low strain level results in a significant redistribution of the scatterers and, consequently, little correlation between the phase of the backscattered light from that section before and after strain.For a short probe pulse, on the other hand, the relative displacement of the scatterers is small since the pulse covers a shorter section of the fibre.Therefore, the redistribution of the scatterers does not significantly affect the phase of the backscattered light from that section.This relationship between the width of the probe pulse and the standard deviation of the measurement can be seen in the diagram of figure 11(b).

Conclusion and future work
The simulation results presented in this study provide an insight into the operation of DVS systems and the application of numerical modelling in determining the limitations of such systems.In this study, the behaviour of the numerical model in response to changes in various sensing parameters was assessed.The simulation results showed that the precision of the measurement is a function of the frequency and amplitude of perturbations.It was shown that the sensing system has a higher precision when measuring dynamic strains with lower amplitude.The precision of the sensing system as a function of the pulse width and linewidth of the probe pulse was also assessed.It was shown that the precision of the sensor can be enhanced by reducing the width of the probe pulse or by increasing its linewidth.
This work provides a useful basis for an in depth analysis of DVS systems that employ ϕ-OTDR to map vibrations along the sensing fibre.The numerical model used in this study can be further expanded to analyse the effect of parameters such as the polarization of the probe pulse, digitization level of analogue-to-digital converter, non-ideal optical components, and noise on the behaviour of the sensing system.In particular, addition of ASE noise, detector noise, and laser phase noise to the numerical model helps to realize a more realistic model of the sensing system.Future work will also aim to juxtapose the numerical analysis against the experimental results to validate the numerical model.

Fig. 1 .
Fig. 1.The flowchart of the numerical model used to simulate the operation of phase-sensitive distributed vibration sensor.

Fig. 2 .
Fig. 2. The schematic of the experimental setup modelled using the new simulation technique.DFB Laser: Distributed feedback laser, EDFA: Erbium-doped fibre amplifier.

Fig. 3 .
Fig. 3. Schematic representation of sensing fibre model.The red circles in the picture represent the scattering points in the fibre.The diameter of the circles represent the size of the scattering points while the location of the circles indicate the random position of the scatterers.

Fig. 5 .
Fig. 5. Mathematical procedure used to convert the backscattered electric field EF to the outputs of the IMZI.The solid arrows indicate through ports while the dashed arrows indicate coupled ports.

Fig. 6 .Fig. 7 .
Fig. 6.The first reliability test results.(a) Backscattered CRN pattern for a probe pulse with 1pm, 10pm, and 100pm linewidths, (b) Probability density function (PDF) of the backscattered CRN for the three linewidths.The simulation results (the solid curves) are juxtaposed with the theoretical analysis (the dashed curves), (c) Changes in the CRN level of the backscattered trace as a function of the square root of the light-source coherence-length.

Fig. 8 .
Fig. 8. Simulation results of the numerically modelled sensing system interrogated by 50cm probe pulse while two regions of the sensing fibre were sinusoidally modulated at 1kHz and 2kHz with an amplitude of 2.5µε and 1.5µε, respectively.(a) Backscattered CRN from the sensing fibre, (b) 3D representation of the differential phase after differentiate and cross-multiply demodulation in time domain, and (c) 3D representation of the differential phase after differentiate and cross-multiply demodulation in frequency domain.

Figure 8 (
b) represents the output of the numerical model following differential and cross-multiply demodulation.The output of the model represents the differential phase between adjacent sections of the sensing fibre 1m apart.This diagram demonstrates the dynamic changes in the strain along the sensing fibre as a function of time.Figure 8(c) provides the 3D representation of the simulation result in the frequency domain.This diagram was obtained by calculating the FFT of the data shown in figure 8(b) along the time axis.Signal processing procedure used to analyse the output of the numerical model was identical to that used for the experimental setup.The two peaks in the 3D plot correspond to the regions on the fibre subjected to dynamic strain.The first peak identifies a 1001Hz signal with an amplitude of 2.57µε positioned around 8.2m.The second peak in the diagram identifies a 2002Hz signal with an amplitude of 1.31µε positioned around 20.4m.

Figure 11 (Fig. 9 .
Fig. 9. Numerical model response to 1kHz sinusoidal strain to a range of strain levels between 100nε and 2µε.(a) Numerical model output as a function of amplitude of induced strain, (b) Standard deviation of different realization of the system model versus the amplitude of induced strain.

Fig. 10 .
Fig. 10.Numerical model response to 0.5µε perturbations for a frequency range of 250Hz to 2000Hz.(a) Numerical model output as a function of the frequency of the induced strain, (b) Standard deviation of different realization of the model versus the frequency of induced strain.
Fig. 11.(a)  Relationship between the standard deviation of the strain level simulated using the numerical model and the linewidth of the probe pulse.For this diagram, the sensing system was modelled for a range of linewidths from 0.1pm to 100pm for a fixed pulse width of 3m.(b) Relationship between the standard deviation of the strain level simulated using the numerical model and the width of the probe pulse.For this diagram, the sensing system was modelled by stepping the pulse widths from 50cm to 2m in 25cm steps while the linewidth of the probe pulse was fixed to 3pm.

Figure 8 (
b) shows the time domain response of the model for the first 4ms of the simulation.The two sinusoidal waveforms corresponding to 1kHz and 2kHz vibrations can be clearly observed.The position and amplitude of the waveforms indicate that the system model can accurately retrieve the vibration parameters used in the simulation.The two peaks in figure8(c) identify the location, frequency, and amplitude of the perturbations.This result shows that there is no cross-talk between the two regions under strain and the rest of the fibre.In addition, it can be seen that the spatial span of the affected regions are the same in all diagrams of figure 8.

Fig. 12 .
Fig. 12. Distribution of inhomogeneities along two fixed sections of the sensing fibre before and after longitudinal elongation.