Patent Publication Number: US-2023146473-A1

Title: Frequency-drift compensation in chirped-pulse-based distributed acoustic sensing

Description:
CROSS REFERENCE 
     This application claims the benefit of U.S. Provisional Pat. Application Serial No. 63/270,199 filed 21OCT. 2021, the entire contents of which being incorporated by reference as if set forth at length herein. 
    
    
     TECHNICAL FIELD 
     This disclosure relates generally to distributed fiber optic sensing (DFOS) systems, methods, and structures and, in particular, to frequency-drift compensation in chirped-pulse-based distributed acoustic sensing (DAS). 
     BACKGROUND 
     Recently, DFOS systems and methods have been employed to provide superior acoustic and / or vibrational monitoring of roadways, bridges, and buildings. The reliability, robustness, and sensitivity of such systems is generally known to be unmatched by existing, legacy systems and methods. Given such characteristics, further improvement of DFOS / DAS techniques - coupled with novel analysis systems and methods would represent a welcome addition to the art. 
     SUMMARY 
     An advance in the art is made according to aspects of the present disclosure directed to frequency drift compensation for coded-DAS systems that use chirped pulses as a probe signal. 
     Frequency drift results in timing jitter in the estimated Rayleigh impulse response, which is obtained by correlating the received Rayleigh backscatter with the originating chirp. Since in any given received frame, neighboring samples will have similar timing jitter, this effect can be compensated by dividing each frame into small blocks where all samples are shifted by the same timing jitter, correlating the amplitude profile from frame to frame, and then re-aligning them in time. 
     Our inventive approach provides an architecture for estimating and correcting timing jitter using an overlap-and-save architecture with low algorithmic complexity and enables coded-DAS systems to employ using “cheaper” lasers with larger frequency drift by using DSP to compensate for the hardware. 
     Viewed from one aspect, our inventive approach estimates timing jitter by correlating the amplitude of the estimated Rayleigh impulse response of every frame with a reference frame, and then re-aligns each frame by the estimated timing jitter. As the amount of timing jitter varies within a frame, every frame is divided into blocks where all samples have similar timing jitter, and perform timing jitter estimation and compensation on a block-by-block, frame-by-frame basis using an overlap-and-save method. Tracking of a slowly changing channel is enabled by allowing the reference frame to be periodically updated. 
     Viewed from a first aspect, the present disclosure describes DFOS systems, methods, and structures for monitoring outdoor cabinets containing fiber optic facilities in which the cabinet / fiber optic cable contained therein are configured to provide superior acoustic sensing. 
     Viewed from a second aspect, the present disclosure describes DFOS systems, methods, and structures for monitoring manhole structures. 
     Finally, viewed from yet another aspect, the present disclosure describes DFOS systems, methods, and structures employing a machine learning-based analysis method that employs a temporal relation network. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWING 
       A more complete understanding of the present disclosure may be realized by reference to the accompanying drawing in which: 
         FIG.  1 (A)  is a schematic diagram illustrating a DFOS system according to aspects of the present disclosure; 
         FIG.  1 (B)  is a schematic diagram illustrating a coded constant-amplitude DFOS system with out-of-band signal generation according to aspects of the present disclosure; 
         FIG.  2 (A)  and  FIG.  2 (B)  are a pair of plots illustrating the impact of frequency shift on the correlation of a chirped pulse for:  FIG.  2 (A)  no frequency shift, and  FIG.  2 (B)  4 MHz frequency shift according to aspects of the present disclosure; 
         FIG.  3    is a schematic diagram illustrating a canonical model of a coded -DAS interrogator where the same laser is used to generate probe signal x(t) and to serve as local oscillator for coherent detection of the Rayleigh backscatter y(t) of the fiber under test (FUT) with Rayleigh impulse response h(t) according to aspects of the present disclosure; 
         FIG.  4    is a schematic diagram showing illustrative architecture for laser frequency drift compensation in coded-DAS according to aspects of the present disclosure; and 
         FIG.  5    is a schematic diagram showing illustrative procedure for compensating the time shift of each block according to aspects of the present disclosure; 
     
    
    
     DESCRIPTION 
     The following merely illustrates the principles of the disclosure. It will thus be appreciated that those skilled in the art will be able to devise various arrangements which, although not explicitly described or shown herein, embody the principles of the disclosure and are included within its spirit and scope. 
     Furthermore, all examples and conditional language recited herein are intended to be only for pedagogical purposes to aid the reader in understanding the principles of the disclosure and the concepts contributed by the inventor(s) to furthering the art and are to be construed as being without limitation to such specifically recited examples and conditions. 
     Moreover, all statements herein reciting principles, aspects, and embodiments of the disclosure, as well as specific examples thereof, are intended to encompass both structural and functional equivalents thereof. Additionally, it is intended that such equivalents include both currently known equivalents as well as equivalents developed in the future, i.e., any elements developed that perform the same function, regardless of structure. 
     Thus, for example, it will be appreciated by those skilled in the art that any block diagrams herein represent conceptual views of illustrative circuitry embodying the principles of the disclosure. 
     Unless otherwise explicitly specified herein, the FIGs comprising the drawing are not drawn to scale. 
     By way of some additional background, we begin by noting that distributed fiber optic sensing (DFOS) is an important and widely used technology to detect environmental conditions (such as temperature, vibration, acoustic excitation vibration, stretch level etc.) anywhere along an optical fiber cable that in turn is connected to an interrogator. As is known, contemporary interrogators are systems that generate an input signal to the fiber and detects / analyzes the reflected/scattered and subsequently received signal(s). The signals are analyzed, and an output is generated which is indicative of the environmental conditions encountered along the length of the fiber. The signal(s) so received may result from reflections in the fiber, such as Raman backscattering, Rayleigh backscattering, and Brillion backscattering. DFOS can also employ a signal of forward direction that uses speed differences of multiple modes. Without losing generality, the following description assumes reflected signal though the same approaches can be applied to forwarded signal as well. 
       FIG.  1 (A)  is a schematic diagram of a generalized, prior-art DFOS system . As will be appreciated, a contemporary DFOS system includes an interrogator that periodically generates optical pulses (or any coded signal) and injects them into an optical fiber. The injected optical pulse signal is conveyed along the optical fiber. 
     At locations along the length of the fiber, a small portion of signal is reflected and conveyed back to the interrogator. The reflected signal carries information the interrogator uses to detect, such as a power level change that indicates - for example -a mechanical vibration. While not shown in detail, the interrogator may include a coded DFOS system that may employ a coherent receiver arrangement known in the art such as that illustrated in  FIG.  1 (B) . 
     The reflected signal is converted to electrical domain and processed inside the interrogator. Based on the pulse injection time and the time signal is detected, the interrogator determines at which location along the fiber the signal is coming from, thus able to sense the activity of each location along the fiber. 
     Those skilled in the art will understand and appreciate that by implementing a signal coding on the interrogation signal enables the sending of more optical power into the fiber which can advantageously improve signal-to-noise ratio (SNR) of Rayleigh-scattering based system (e.g. distributed acoustic sensing or DAS) and Brillouin-scattering based system (e.g. Brillouin optical time domain reflectometry or BOTDR). 
     Operationally, we assume that the DFOS system will be Rayleigh-scattering based system (e.g., distributed acoustic sensing or DAS) and Brillouin-scattering based system (e.g., Brillouin optical time domain reflectometry or BOTDR) with a coding implementation. With such coding designs, these systems will be most likely be integrated with fiber communication systems due to their lower power operation and will also be more affected by the optical amplifier response time. 
     In the arrangement illustratively shown in the block diagram, we assume that the coded interrogation sequence is generated digitally and modulated onto the sensing laser via digital-to-analog-conversion (DAC) and an optical modulator. The modulated interrogation sequence may be amplified to optimal operation power before being directed into the fiber for interrogation. 
     Advantageously, the DFOS operation may also be integrated together with communication channels via WDM in the same fiber. Inside the sensing fiber, the interrogation sequence and the returned sensing signal may be optically amplified - either via discrete (EDFA/SOA) or distributed (Raman) methods. A returned sensing signal is routed to a coherent receiver after amplification and optical band-pass filtering. The coherent receiver detects the optical fields in both polarizations of the signal, downconverting them to 4 baseband lanes for analog-to-digital conversion (ADC) sampling and digital signal processor (DSP) processing. As those skilled in the art will readily understand and appreciate, the decoding operation is done in the DSP to generate the interrogated Rayleigh or Brillouin response of the fiber, and any changes in the response are then identified and interpreted for sensor readouts. 
     With continued reference to the figure, since the coded interrogation sequence is generated digitally, the out-of-band signal is also generated digitally, and then combined with the code sequence before waveforms are created by the DAC. When generated together digitally, the out-of-band signal will only be generated outside the time period of the code sequence, so when added together, the combined waveform will have a constant amplitude. 
     In distributed acoustic sensing (DAS), the interrogator launches a probe signal x(t) into the fiber under test (FUT) to estimate its Rayleigh impulse response h(t) . The received signal is given by the convolution y(t) = x(t)⊗h(t) . Conventionally,  
     
       
         
           
             x 
             
               t 
             
             = 
             
               P 
             
               
             rect 
             
               
                 t 
                 / 
                 T 
               
             
           
         
       
     
      is a pulse with peak power P and duration T. 
     In order to measure time-variation in the Rayleigh impulse response caused by acoustic vibration, the interrogator transmits x(t) periodically at a frame rate of T p . Provided T p  is longer than the round-trip propagation time T rt  of the FUT, the received signal is a sequence of optical time-domain reflectometry (OTDR) traces. 
     The spatial resolution of y(t) is z res  = (c/2n eƒƒ )T and is controlled by the bandwidth of x(t). By measuring the phase relationship between pairs of points at positions z 1  = (c/2n eƒƒ )t 1  and z 2  = (c/2n eƒƒ )t 2  in the received OTDR, and how that phase varies with time, the time-varying longitudinal strain between z 1  and z 2  can be monitored. This technique of phase-OTDR (ɸ-OTDR) is well known in the art. 
     Due to the weak power of Rayleigh backscatter, the reach achievable by DAS is limited. Coded-DAS increases the signal-to-noise (SNR) ratio achieved by DAS to allow for longer FUT. Instead of launching pulses - and as previously noted - the probe signal in coded-DAS is a sequence with autocorrelation function r xx (t) = x(t) ∗ x(t) that is as close to a delta function as possible, with its width constrained only by the bandwidth of x(t). The coded-DAS interrogator performs correlation of the received signal with x(t) to obtain z(t) = x(t) ∗ y(t) = rxx(t)⊗h(t) . Since rxx(0) =  
     
       
         
           
             
               
                 
                   ∫ 
                   
                     − 
                     ∞ 
                   
                   
                     + 
                     ∞ 
                   
                 
                 
                   
                     
                       
                         
                           x 
                           
                             t 
                           
                         
                       
                     
                     2 
                   
                 
               
             
           
         
       
     
      dt = P ▪ T c  is the energy of x(t), the longer the sequence, the larger the received signal. Advantageously, SNR can be increased linearly with T c  without sacrificing spatial resolution. Laser phase noise is the limiting factor for how large T c  can be. 
     One well-known family of sequences with good autocorrelation properties are chirped pulses:  
     
       
         
           
             x 
             
               t 
             
             = 
             
               P 
             
             exp 
             
               
                 j 
                 2 
                 π 
                 α 
                 
                   
                     
                       t 
                       2 
                     
                   
                   2 
                 
               
             
             rect 
             
               
                 
                   t 
                   
                     
                       T 
                       c 
                     
                   
                 
               
             
           
         
       
     
      where T c  is the chirp duration, α is the chirp rate, and  
     
       
         
           
             
               P 
             
           
         
       
     
      is the amplitude of the envelope. The bandwidth of the x(t) is B = αT c , and its autocorrelation is given by:  
     
       
         
           
             
               R 
               
                 x 
                 x 
               
             
             
               t 
             
             = 
             x 
             
               t 
             
             ∗ 
             x 
             
               t 
             
             = 
               
             
               
                 
                   T 
                   c 
                 
                 − 
                 
                   t 
                 
               
             
               
             sin 
             ​ 
             c 
             
               
                 α 
                 t 
                 
                   
                     
                       T 
                       c 
                     
                     − 
                     
                       t 
                     
                   
                 
               
             
             , 
           
         
       
     
     For long chirp duration T c  » T, the width of the main lobe of the sinc in Eq. (2) is T = 1/αT c  = 1/B, which is the same spatial resolution as a rectangular pulse of the same bandwidth used in conventional OTDR. 
     Chirped pulses have the special property that the correlation function between two chirped pulses x 1 (t) and x 2 (t) is only the sinc function in Eq. (2) if their chirp rates α 1  and α 2  match. If this is the case, their correlation peak occurs where their center frequency match. 
       FIG.  2 (A)  and  FIG.  2 (B)  are a pair of plots illustrating the impact of frequency shift on the correlation of a chirped pulse for:  FIG.  2 (A)  no frequency shift, and  FIG.  2 (B)  4 MHz frequency shift according to aspects of the present disclosure. Chirps have the special property that a correlation peak is produced when their instantaneous frequencies align. Frequency modulation of one chirped pulse thus results in a temporal shift of the correlation peak (and also a slight broadening of the main lobe of the correlation function due to reduced bandwidth overlap). 
     An example is shown in  FIG.  2 (A)  and  FIG.  2 (B)  where a chirped pulse of duration T c  = 50 µs and bandwidth B = 10 MHz (α = 2 × 10 11  s -2 ) is correlated with the same chirped pulse that is frequency-shifted by Δv = 4 MHz. It is observed that their correlation x 2 (t) ∗ x 1 (t) is centered at an offset of Δτ = Δν/α = 20 µs . In addition, the width of the correlation is inversely proportional to their bandwidth overlap, which in the example shown is 6 MHz. 
     This property of chirped pulses, where frequency modulation results in temporal shift of the center of the correlation function is important to coded-DAS based on chirped pulses. 
       FIG.  3    is a schematic diagram illustrating a canonical model of a coded -DAS interrogator where the same laser is used to generate probe signal x(t) and to serve as local oscillator for coherent detection of the Rayleigh backscatter y(t) of the fiber under test (FUT) with Rayleigh impulse response h(t) according to aspects of the present disclosure. 
     Consider the canonical model of a DAS system shown in  FIG.  3   , where the same laser is used to generate the probe signal and to serve as local oscillator (LO) for coherent detection of the Rayleigh backscatter. Consider the reflection from the point shown in the FUT in  FIG.  3   . Since there is a propagation delay difference of τ z  = (2n eff /c)z between the reflected signal and the LO, the laser frequency will have drifted over duration τ z . Frequency drift is the manifestation of laser phase noise in the low frequency region and is caused by factors such as temperature fluctuation and mechanical vibration, etc., which causes fluctuation in the optical length of the laser cavity, thus impacting the laser’s center frequency. 
     As per  FIG.  2 (A) , and  FIG.  2 (B) , the impact of frequency drift is uncertainty in time (or position) of the Rayleigh impulse response. Normally, the fiber position z is mapped to time coordinate t = (2n eff /c)z in the Rayleigh impulse response h(t). But due to laser frequency drift, the actual fiber position at time t is z̅ + δz where z̅ = (c/2n eyy )t is the mean position, and δZ is a variable that changes from frame to frame, depending on the laser frequency drift Δv during that frame. When DAS computes a &lt;p-OTDR based on the differential product between two points separated by gauge length z g , instead of monitoring the phase difference between two fixed scatterers at z 1  = z̅ – z g /2 and z 2  = z̅ + z g/ 2, the actual phase difference is computed for two scatterers at z 1  = z̅ + δz – z g /2 and z 2  = z̅ + δz + z g /2. If δZ is larger than the spatial resolution, the resulting DAS will become non-sensical, since scatters separated by larger than the spatial resolution should be independent. 
     The impact of laser frequency drift on the performance of DAS is understood in an experiment in which a 50-km long FUT was probed using chirped pulses of duration T c  = 50 µs and bandwidth B = 10 MHz (spatial resolution z res  = c/2n eff  B ≈ 10 m), at a repetition rate of 1 kHz. At 50 km, a piezo-electric transducer (PZT) was inserted, followed by another 100 m of termination fiber. The PZT was excited with a 67-Hz sine wave with peak-to-peak amplitude of 1.6 rad. 
     When one employs a specific (“bad”) laser, the frequency drift results are so severe that the background noise level is larger than the vibration amplitude of the PZT. With a different (“good”) laser however, the 67-Hz sine wave is easily observed. 
     Although it is possible to build a chirped pulse DAS using the “good” laser, it is more expensive. Since frequency drift only results in linear translation of the time/position axis, and fiber positions close together should experience similar time/position shift as the laser frequency should be stable over short timescales, it should be possible to compensate this time (position) jitter using digital signal processing (DSP), and therefore allow DAS using the “bad” laser with higher phase noise. 
     Notwithstanding, the present disclosure provides a solution for frequency drift compensation for coded-DAS systems that use chirped pulses as the probe signal. Frequency drift results in timing jitter in the estimated Rayleigh impulse response, which is obtained by correlating the received Rayleigh backscatter with the originating chirp. 
     Since in any given received frame, neighboring samples will exhibit similar timing jitter, our inventive disclosure compensates for this effect by dividing each frame into small blocks where all samples are shifted by the same timing jitter, correlating the amplitude profile from frame to frame, and then re-aligning them in time. Our disclosed solution provides an architecture for estimating and correcting timing jitter using an overlap-and-save architecture, and also provides preferred implementations with low algorithmic complexity. Our disclosed solution enables coded-DAS using “cheaper” lasers with larger frequency drift by using DSP to compensate the hardware deficiencies. 
     As will become apparent to those skilled in the art, features of our disclosed solution according to aspects of the present disclosure include: (i) estimation of timing jitter by correlating the amplitude of the estimated Rayleigh impulse response of every frame with a reference frame, and then (ii) re-aligning each frame by the estimated timing jitter. 
     As the amount of timing jitter varies within a frame, our inventive approach divides every frame into blocks where all samples exhibit a similar timing jitter, and then perform timing jitter estimation and compensation on a block-by-block, frame-by-frame basis using an overlap-and-save method. Our inventive approach also enables tracking of a slowly changing channel by allowing the reference frame to be periodically updated. 
       FIG.  4    is a schematic diagram showing illustrative architecture for laser frequency drift compensation in coded-DAS according to aspects of the present disclosure. As illustrated, the architecture shown in  FIG.  4    is for chirped pulses using an overlap-and-save algorithm. Operationally, after correlating the received signal with the known chirp to obtain h j (t) = x(t) ∗ y j (t) in the j-th frame, the time axis is partitioned into N b  overlapping blocks of duration T b  and overlap T ov . Taking frame j = 0 as reference, the square amplitude of each block |h j (t b )| 2  in subsequent frames is compared with |h 0 (t b )| 2  to find their relative offset τ j , b . Each block is then shifted by τ j,b  to obtain the compensated outputs h j (t b ). Normal DAS operations are then performed on h̃ j (t). 
     To further clarify, we let y j (t) be the complex-valued signal vector received due to the probe pulse launched in frame j. Let h j (t) = x(t) ∗ y j (t) be the output of the correlator. The duration of each vector is T rt  equal to the round-trip return time of the FUT. We divide h j (t) into overlapping blocks of duration T B  and overlap T ov . The number of overlapping blocks N b  needs to satisfy (N b  – 1)(T b  – T ov ) + T b  ≥ T rt . Let t b  ε {t: b(Tb – T ov ) ≤ t &lt; b(Tb – T ov ) + Tb} be the duration of block b, where 0 ≤ b &lt; N b , 
     We then take the initial frame h o (t) to be the reference. For every block b, we correlate the amplitude |h j (t b )| 2  of all subsequent frames with the reference |h 0 ( tb )| 2  to estimate the timing jitter τ j,b  = max ∫ t   b  |h j (t  _   τ )| 2  |h 0 (t)| 2  dt during block b of frame j. We then shift h j (t b ) by τ j,b  to obtain the compensated signal h j (t b ) = h j (t b  – τ j,b ). The useful part of the overlap-and-save, b(Tb – T ov ) ≤ t &lt; (b + 1)(T b  – T ov ), is then stored. The frequency-drift compensated signal is h̃(t). 
     The downstream operations needed to estimate vibration at every point in the FUT are then the same as conventional DAS. These operations may include calculating differential beat products at a pre-defined gauge length, diversity-combining the beat products from different polarizations, frequencies, and spatial channels, etc., and then finally taking the unwrapped phase at every location in the FUT. 
     According to our inventive disclosure, the block length T b  should be chosen so that the root mean square (r.m.s.) frequency drift (σ Δv (T b ) over duration T b  is much less than time resolution of the chirp T = 1/B, i.e., the samples in h j (t b ) are all delayed relative to h o (t) by a similar amount down to an accuracy of the time-resolution 1/B. T b  will depend on the laser used by the interrogator. 
     Suppose the laser’s instantaneous frequency is v(t). Let S vv (f ) be the two-sided frequency noise spectrum of is v(t). The frequency drift Δv T   b (t) = v(t) - v(t – T b ) between two time instances separated by T b  has the Fourier transform  
     
       
         
           
             Δ 
             
               v 
               
                 
                   T 
                   b 
                 
               
             
             
               f 
             
             = 
             
               
                 1 
                 − 
                 
                   e 
                   
                     − 
                     j 
                     2 
                     π 
                     f 
                     
                       T 
                       b 
                     
                   
                 
               
             
             v 
             
               f 
             
               
             . 
           
         
       
     
     . Hence, the spectrum of Δv T   b (t) is 
     
       
         
           
             
               S 
               
                 Δ 
                 
                   v 
                   
                     
                       T 
                       b 
                     
                   
                 
               
             
             
                 
               
                 Δ 
                 
                   v 
                   
                     
                       T 
                       b 
                     
                   
                 
               
             
             ( 
             f 
             ) 
             = 
             
               
                 
                   
                     1 
                     − 
                     
                       e 
                       
                         − 
                         j 
                         2 
                         π 
                         f 
                         
                           T 
                           b 
                         
                       
                     
                   
                 
               
               2 
             
             
               S 
               
                 v 
                 v 
               
             
             ( 
             f 
             ) 
             . 
           
         
       
     
     It is well known that the variance  
     
       
         
           
             E 
             
               
                 
                   
                     
                       
                         Δ 
                         
                           v 
                           
                             
                               T 
                               b 
                             
                           
                         
                         
                           t 
                         
                       
                     
                   
                   2 
                 
               
             
           
         
       
     
      is equal to the autocorrelation  
     
       
         
           
             
               R 
               
                 Δ 
                 
                   v 
                   
                     
                       T 
                       b 
                     
                   
                 
                 Δ 
                 
                   v 
                   
                     
                       T 
                       b 
                     
                   
                 
               
             
             
               τ 
             
           
         
       
     
      evaluated at τ = 0, which is equal to the integral of  
     
       
         
           
             
               S 
               
                 Δ 
                 
                   v 
                   
                     
                       T 
                       b 
                     
                   
                 
               
             
             
                 
               
                 Δ 
                 
                   v 
                   
                     
                       T 
                       b 
                     
                   
                 
               
             
             ( 
             f 
             ) 
           
         
       
     
      over all frequencies. Hence,  
     
       
         
           
             E 
             
               
                 
                   
                     
                       
                         Δ 
                         
                           v 
                           
                             
                               T 
                               b 
                             
                           
                         
                         
                           t 
                         
                       
                     
                   
                   2 
                 
               
             
             = 
             
               
                 
                   
                     ∫ 
                     
                       
                           
                         
                           − 
                           ∞ 
                         
                         
                           + 
                           ∞ 
                         
                       
                     
                   
                 
               
             
             1 
             − 
           
         
       
     
     
       
         
           
             
               
                 
                   
                     
                       e 
                       
                         − 
                         j 
                         2 
                         π 
                         f 
                         
                           T 
                           b 
                         
                       
                     
                   
                 
               
               2 
             
             
               S 
               
                 v 
                 v 
               
             
             
               f 
             
               
               
             d 
             f 
             . 
           
         
       
     
      Since Δv T   b (t) is a Gaussian process, T b  should be chosen so that  
     
       
         
           
             
               σ 
               
                 Δ 
                 v 
               
             
             
               
                 
                   T 
                   b 
                 
               
             
             = 
             
               
                 E 
                 
                   
                     
                       
                         
                           
                             
                               Δ 
                               
                                 
                                   V 
                                   
                                     
                                       T 
                                       b 
                                     
                                   
                                 
                               
                             
                           
                         
                       
                       2 
                     
                   
                 
               
             
             = 
             
               1 
               / 
               
                 ς 
                 B 
                 , 
               
             
           
         
       
     
      to ensure the probability  
     
       
         
           
             erfc 
             
               
                 
                   ς 
                   / 
                   
                     
                       2 
                     
                   
                 
               
             
           
         
       
     
      is small that the  samples T b  at the beginning and end of a block have relative timing jitter greater than 1/B. 
     Similarly, the overlap duration T ov  should be larger than twice the maximum value of |τ j,b |that may need to be compensated. As per  FIG.  3   , the r.m.s. value of the timing jitter is σ Δv (t)/a, which increases with distance from the interrogator z = (c/2n eyy )t,. If a constant overlap size is used, then the overlap duration should be chosen as  
     
       
         
           
             
               T 
               
                 o 
                 v 
               
             
             = 
             K 
             
               
                 
                   σ 
                   
                     Δ 
                     v 
                   
                 
                 
                   
                     
                       T 
                       
                         r 
                         t 
                       
                     
                   
                 
               
               α 
             
             , 
           
         
       
     
      to ensure the probability  
     
       
         
           
             erfc 
             
               
                 
                   κ 
                   / 
                   
                     
                       2 
                     
                   
                 
               
             
           
         
       
     
      is small that the time jitter τ j,N   n  of the last block of any frame is larger than the overlap duration. 
     In practice, the Rayleigh impulse response h(t) that is estimated by h j (t) is slowly changing with time due to polarization rotation, temperature fluctuation, etc. If the reference h 0 (t) is held for too long, then eventually, h(t) will have changed so much that |h 0 (t)| 2  is no longer a good match for the incoming frames. Hence, the reference should be updated regularly every N u  frames. i.e., the mN u -th timing-jitter-compensated frame h̃ mn   u (t) should be used as the new reference for compensating frames h mN   u+1 (t) to h(m +1 ) N   u (t). 
       FIG.  5    is a schematic diagram showing illustrative procedure for compensating the time shift of each block according to aspects of the present disclosure. Correlation with the reference is first computed to find the index n j ,b at which r jo , b [m] is maximized over – N ov /2 ≤ m ≤ N ov /2 . The block is then barrel shifted by n j , b  samples. 
     If more accurate compensation is required, the ratios between the two neighboring samples with the peak sample can be computed, and a lookup table can be used to find the fractional shift n&#39; j,b . Fractional shifting of the block is achieved using an interpolation filter. 
     Reduced Complexity Implementation 
     Per the aforementioned description, the computationally expensive operations are: (a) correlating each block |h j (t b )| 2  with |h 0 (t b )| 2 , and (b) shifting each h j (t b ) by τ j,b . In practice, these signals are sampled at a rate of T S  = 1/MB where M is the oversampling ratio relative to the chirp bandwidth B. Similar to digital coherent receivers for telecommunications, M = 2 is often chosen as compromise. Note that τ j,b  does not have to be an integer multiple of T s . 
     The digital correlation  
     
       
         
           
             
               r 
               
                 j 
                 0 
                 , 
                 b 
               
             
             
               m 
             
             = 
             
               
                 ∑ 
                 
                   n 
                   ∈ 
                   
                     n 
                     b 
                   
                 
               
               
                 
                   
                     
                       
                         
                           h 
                           j 
                         
                         
                           
                             n 
                             − 
                             m 
                           
                         
                       
                     
                   
                   2 
                 
                 
                   
                     
                       
                         
                           h 
                           0 
                         
                         
                           n 
                         
                       
                     
                   
                   2 
                 
               
             
           
         
       
     
      is expected to be a sinc 2 ( . ) function sampled at T s . A computationally simple solution is to compute r jo , b [m] only for a range of indices – N ov /2 ≤ m ≤ N ov /2, where N ov  = [T ov /T] is the overlap duration in number of samples, and then find the index n j ,b which maximizes r jo,b [m] over this range. The signals h j  [n b ] are then barrel shifted by n j ,b samples, which does not require any multiplications. This simplified implementation aligns all h j [n b ] to the nearest sample, ensuring the compensated signals have residual timing jitter of at most ±0.5T S . 
     Note that it is possible to reduce timing jitter further with an alternative implementation having slightly higher computational complexity. Since  
     
       
         
           
             
               r 
               
                 j 
                 0 
                 , 
                 b 
               
             
             
               m 
             
             = 
           
         
       
     
     
       
         
           
             
               
                 ∑ 
                 
                   n 
                   ∈ 
                   
                     n 
                     b 
                   
                 
               
               
                 
                   
                     
                       
                         
                           h 
                           j 
                         
                         
                           
                             n 
                             − 
                             m 
                           
                         
                       
                     
                   
                   2 
                 
                 
                   
                     
                       
                         
                           h 
                           0 
                         
                         
                           n 
                         
                       
                     
                   
                   2 
                 
               
             
           
         
       
     
      is a sampled sinc 2 ( . ) function with the width between first nulls ±M samples apart, we can take the index n j ,b which maximizes r jo , b [m], along with its two neighbors n j ,b – 1 and n j ,b + 1. 
     A lookup table indexed by the two ratios r j0,b [n j , b  – 1]/r j,0 , b [n j,b ] and r jo , b [n j,b  + 1]/r j0,b [n j,b]  quantized to a given timing jitter accuracy can then be used to find the fractional sample offset n&#39; j,b&#39;  where  
     
       
         
           
             
               n 
               
                 j 
                 , 
                 b 
               
             
             + 
             
               
                 n 
                 ′ 
               
               
                 j 
                 , 
                 b 
               
             
             = 
             
               
                 
                   τ 
                   
                     j 
                     , 
                     b 
                   
                 
               
               
                 
                   T 
                   s 
                 
               
             
             . 
           
         
       
     
      As with the previous case, h j [n b ] is first barrel-shifted by n j ,b samples. To shift the remaining fractional sample n&#39; j,b&#39;  an interpolation filter h int [n] of N int  taps can be used. The coefficients for h int [n] can once again be stored in a lookup table for different fractional offsets, quantized to a timing jitter accuracy required. This step is analogous to timing recovery in digital coherent receivers for telecommunications, where the received signal is digitally resampled using interpolation filters to keep the samples synchronized with the symbol clock. As N int  can be small (e.g., 3 or 5 taps), resampling using interpolation filters can be of low complexity equal to N int  complex multiplications per sample per polarization. 
     Example Experimental Results 
     We have produced experimental results demonstrating the operation of the frequency-drift compensation according to the present disclosure. The experimental setup is as shown in  FIGS.  2 A  2 B , where chirped pulses of duration T c  = 50 µs and bandwidth B = 10 MHz (spatial resolution z res  = c/2n eff  B ≈ 10 m), at a repetition rate of 1 kHz are generated using a digital-to-analog converter (DAC) which drives a Mach-Zehnder I/Q modulator. 
     The FUT comprises of a 50 km spool of standard single-mode fiber (SSMF), followed by a 12-m long piezo-electric transducer (PZT), whose output is terminated by a 100-m long fiber. The PZT is driven with a 67-Hz sine wave with peak-to-peak amplitude of 1.6 rad. It was observed that distortion caused by frequency drift is so severe that the resulting background noise level is larger than the vibration amplitude of the PZT 
     To evaluate timing jitter associated with laser frequency drift, we evaluated over 200 traces of |h j (t) | 2  recovered after correlating the received signal with the chirp. An exponential decay profile due to fiber attenuation, which has a round-trip loss of ≈ 0.4 dB/km was produced. As expected, there is negligible timing jitter between the beginning of the fiber, but timing jitter of around 1.5 µs peak-to-peak is observed at the end of the FUT, which corresponds to a frequency drift of ~ 300 kHz peak-to-peak. 
     With this jitter characterized, we used our inventive method to compensate frequency drift. Block sizes of T b  = 25 µs and overlap duration of T ov  = 2.5 µs were used. When we applied this technique to data as above, |h j (t)| 2  near the end of the after timing-jitter compensation, the 200 traces become well aligned with each other. 
     When our technique is applied to performance of DAS using the “cheap/bad” laser previously described with frequency drift compensation, the 67-Hz vibration is now clearly observed and the phase spectrum has a noise floor corresponding to a strain sensitivity of 38  
     
       
         
           
             p 
             ∈ 
             / 
             
               
                 Hz 
               
             
             . 
           
         
       
     
      Finally, the performance of DAS using a “good” laser that exhibits much less frequency drift than the “bad” laser produces a strain sensitivity at a position of the PZT also around 38  
     
       
         
           
             p 
             ∈ 
             / 
             
               
                 H 
                 z 
               
             
             . 
           
         
       
     
      Thus, frequency-drift compensation enables the cheaper “bad” laser with larger frequency drift to achieve comparable DAS performance as the expensive “good” laser 
     At this point, while we have presented this disclosure using some specific examples, those skilled in the art will recognize that our teachings are not so limited. Accordingly, this disclosure should only be limited by the scope of the claims attached hereto.