A Refinement of Backward Correlation Technique for Precise Brillouin Frequency Shift Extraction

Abstract

A new method for extracting the Brillouin frequency shift (BFS) from the Brillouin gain spectrum (BGS), the modified backward correlation method (MBWC), is presented. The possibilities of using MBWC, and MBWC in combination with the Lorentzian curve fitting (LCF) based on Levenberg–Marquardt (LM) method, are studied. The effectiveness of the new method, and its combination with LM, has been demonstrated for processing spectra with a low signal-to-noise ratio (SNR). The experiments, which were in good agreement with the performed simulation, showed that at SNR = 0 dB, the combined (MBWC + LM) method provided the BFS extraction error of less than 4 MHz, while the state-of-the-art LM algorithm extracted it with the error greater than 4.5 MHz. The advantage of correlation methods becomes more significant with the decreasing SNR: at SNR = −2 dB, the LM’s error is 14.3 MHz, and that of the combined one is 8.1 MHz.

1. Introduction

Currently, there are quite a lot of relevant scientific and technical challenges related to distributed monitoring of various object parameters [1]: buildings, devices made of new materials, vehicles, elements of urban infrastructure, roads, bridges, tunnels, pipelines, and mines. The key parameters to be monitored are temperature, deformations [2], vibration frequency ranges [3], gas concentrations in the medium [4] in which the object under study is located, and some other characteristics. The classic method for observing these quantities is to use traditional pointwise sensors based on piezoelectric elements, thermocouples, laser distance meters, and other well-known technologies. Fiber-optic sensors are also used, both pointwise and distributed [5]. Pointwise fiber-optic sensors, which are mainly fiber Bragg gratings [6], as in the case of the classical approach, do not provide a complete distributed measurement since the study of object characteristics occurs only at those points of the fiber sensor where the necessary reflective structure is formed. If it is necessary to obtain information about each point of the sensor, it is advisable to use distributed fiber-optic sensors [7,8,9,10]. They can be based on Rayleigh [11], Raman [12], or Brillouin [13,14,15] scattering. The intensity of Rayleigh scattering in an optical fiber weakly depends on the impact applied; therefore, special methods of radiation phase extraction are used to obtain information about temperatures and strains [16]. This approach requires rather complicated signal processing and imposes significant requirements on the laser source: both in terms of the coherence length and in the sweeping linearity [17]. Raman scattering is sensitive to temperatures only, so it is possible to observe temperatures and strains in a fiber simultaneously only with duplication of measurements by another method [18]. Brillouin scattering, which is a consequence of the scattering of a photon by an acoustic phonon because of electrostriction, has a well-studied spectral component of the Lorentzian shape, the peak frequency of which is proportional to the temperature and deformation effects exerted on the sensor [19]. Thus, the obtaining of necessary parameters could be carried out by a simple proportional recalculation of the frequency extracted from the measurement data. However, the main obstacle to obtaining the correct temperature or fiber strain values is the spectrum distortion caused by the noise of the registering photodiodes, thermal fluctuations of the material molecules, the interaction between the probe and pump radiation, and other factors [20]. The extremely weak intensities of the Brillouin scattering components make the problem even more difficult.

For the correct extraction of the Brillouin frequency shift (BFS), various methods are used: from Lorentzian curve fitting (LCF) to the use of correlation algorithms, image processing, and neural network solutions. Some of these techniques will be discussed in detail in later parts of this manuscript.

The references describe in detail the comparison of individual methods for detecting BFS at different signal-to-noise ratios (SNR) of Brillouin gain spectra (BGS) [21,22]. There are also papers devoted to the use of various methods in conjunction with each other. Thus, in [23], the positive and negative aspects of the joint use of correlation methods and artificial intelligence are evaluated.

The backward correlation method (BWC) was proposed in [24]. The essence of the method is to inverse the spectrum relative to the center frequency of the scanning range and find such a shift of the reflected spectrum relative to the original one, in which the Lorentzian peaks overlap each other as much as possible. The farther the maximum of the Lorentzian peak is from the center frequency of the range, the greater the BFS is required for superposition. In [24], the degree of peak overlap was proposed to be calculated from the convolution value of the original and reflected shifted spectra. Determining the optimal shift F, at which the convolution maximum was reached, the authors found the position of the maximum of the Brillouin spectrum as f_C + F/2, where f_C is the center frequency in the scanning range (Figure 1).

Theoretically, the dependence of the original and shifted Lorentzian function’s convolution value on the value of the shift between them is also a Lorentzian function with twice the half-width at half-height (for a derivation, see, for example, [23]). Additionally, since the role of the noise terms decreases when convolution is taken due to statistical averaging, this resulting Lorentzian function (hereinafter referred to as the correlogram) should have a greater SNR than the original BGS. Accordingly, it is processing by existing methods can make it possible to extract the position of the Lorentzian peak with greater accuracy than the processing of the original spectrum. While we were conducting experiments for [23], such an attempt was made to combine BWC with the LCF method, but this approach did not bring the improved accuracy of BFS extraction. The reason for this is the fact that the real correlogram differs from the ideal Lorentzian function due to the following factors: the discreteness of the original spectrum, the presence of a noise component, and the limited frequency scanning range.

This work is devoted to studying the possibilities of improving the BWC method for its further joint use in conjunction with the LCF method.

2. Methodology

The most non-obvious point in the BWC method is the question of how to shift the spectrum correctly. Let the original spectrum be a set of pairs of points (f_i, P_i), where f_i is the frequency in the spectrum, P_i is the corresponding power spectral density, i ∈ [0,N], and N + 1 is the number of points in the spectrum. Then, the inversed spectrum could be written as (f_N−i, P_N−i), i ∈ [0,N]. When the spectrum is shifted, say, to the left by k elements, then the k-th element of the reflected spectrum will go to the 0th position, k + 1—to the 1st position, etc., but there are no elements in the unshifted spectrum, corresponding to k last elements in the shifted spectrum. Accordingly, in [24], it was proposed to fill in these missing values with zeros.

If the original spectrum was the sum of an ideal Lorentzian function and a side signal (total power spectral density minus Lorentzian power spectral density), and the average value of the side signal was strictly equal to zero, then such a replacement by zeros would not greatly distort the spectrum during a shift. However, the Brillouin optical time domain analyzer (BOTDA) can only approximately determine what is a useful signal and what is a side signal, especially if the useful signal is small (it is difficult to determine the signal at low SNRs) or, conversely, very large (at high SNRs, it is difficult to determine the noise). Therefore, having a spectrum taken with BOTDA, one cannot be sure that the average side signal equals zero. If the average side signal is not zero, the replacement by zeros distorts the spectrum during the shift.

In this case, the manifestation of the following negative effect is possible [25]—the maximum of the correlogram stops being corresponding to the real position of the Lorentzian peak but shifts to the middle of the BGS. This happens because at zero shift, the spectrum contains the smallest number of artificially introduced zero elements. (Figure 2) Accordingly, the BWC method stops working.

In [25], it was proposed to use Pearson correlation instead of spectra convolution. The search for the correlation maximum made it possible to achieve the same BFS determination accuracy as in [24] for side signals with an average value. However, for the purposes of this work, this approach is not suitable since the Pearson correlogram does not resemble the Lorenz function (Figure 3).

If we knew the exact average value of the side signal, we could subtract it from the spectrum in advance and get the “correct Lorentzian” correlogram (Figure 4a). At the same time, excessive subtraction manifests itself in shifting the correlogram into the negative region (Figure 4b, similar to the dependence in Figure 3).

Thus, to obtain the “correct” correlogram, the authors developed the following iterative algorithm:

(1)

Calculate the average value R of the power spectral density over the entire spectrum. Since it is equal to the sum of the average value of the useful signal S and the average value of the side signal N, it will certainly be higher than N.

(2)

The test value for subtracting Q is set to R.

(3)

The value of Q is subtracted from all points of the spectrum, and a correlogram is constructed. During the construction, the number of the resulting function negative values is counted.

(4)

If there are no negative values (we received a correlogram as in Figure 2), then the subtracted value was insufficient. If there are too many negative values (we have a correlogram as in Figure 4B), then we subtract excessively. In both cases, we proceed to step 5.

Otherwise, we consider that we have subtracted the desired value and exited the loop. Since the real correlogram always has some noise, in practice, the loop was excited when more than 1% and less than 5% of the correlogram points fell into the negative region.

(5)

Q value is corrected by the method of half division of the segment [0,R], and we go back to step 3.

An example of the algorithm operation is shown in Figure 5. In this case, it took four iterations.

Next, the following methods for finding the maximum of BGS were compared:

MBWC—by finding the maximum of the obtained correlogram.

LM, Levenberg–Marquardt—by approximating the original spectrum with a Lorentzian function.

Combined—by approximating the resulting correlogram with a Lorentzian function.

The last two methods used one of the most promising approaches to approximating the Lorentzian function—the Levenberg–Marquardt algorithm. The method is an iterative search for the ideal Lorenzian function, which would approximate the available data with a minimum error. The least squares method is used to calculate the approximation error.

The essence of the LM method, especially as applied to BGS, is considered in detail in [26]. The Lorentzian function model is based on the Cauchy (Cauchy–Lorentz) distribution, the main parameters of which are: amplitude (A), width (w), and maximum position (µ). The expression for building the model can be represented as follows:

$$f(x:A,\mathsf{\mu },w)=\frac{A}{\pi }\left(\frac{w}{{(x-\mathsf{\mu })}^2+w^2}\right)$$

(1)

The LM algorithm is used to select these parameters and minimize the error of the approximating function. This algorithm can be thought of as a combination of the gradient descent algorithm and the Gauss-Newton (GN) algorithm. LM has local convergence GN and the main advantages of the gradient descent algorithm [27]. However, as mentioned earlier, the algorithm has its drawbacks, such as the need to select initial parameters, which can perform a significant role in the accuracy and speed of the entire system.

Comparison of the three methods was carried out both on the generated and experimental data (Figure 6)

3. Simulation

For the given scanning step and SNR, a Lorentzian spectrum was generated in accordance with (1) [26], where Brillouin width w was taken as much as 40 MHz. The scanning ranges from f_s to f_f, and the scanning step was taken. In the central part of this range, the BFS f_b (equals µ in Equation (1)) of the BGS was randomly selected.

Then, the noise component was added. To do this, we analyzed previously measured spectra from [22] and studied the characteristics of signal distribution. The received signal was cleaned from noise by means of multiple averaging. The statistical characteristics of the noise are presented in Figure 7 below.

As can be seen from the graph, it has a normal distribution in the frequency domain. To rigorously confirm the hypothesis about the normal distribution of the data, the Kolmogorov–Smirnov test was used at a significance level of 0.05. The fact that the distribution is normal made it easier to generate a random signal component.

The SNR parameter was defined as ten natural logarithms of the maximum signal $\frac{1}{\pi w}$ to the RMS noise ratio. An error in determining the BFS was defined as the modulus of the difference between the given f_b and found values.

4. Experiment

A series of SMF-28 optical fiber measurements were carried out using a commercial BOTDA analyzer. To reduce measurement instabilities associated with a random change in ambient temperature, the optical fiber was thermostated in a climatic chamber at a temperature of 25 °C (scheme, Figure 8).

Figure 8 is a simplified structure of the BOTDA Omnisens DiTeSt STA-R202 Brillouin analyzer. This instrument has earned a reputation among users for its high accuracy and ease of measurement. A standard analyzer of this type consists of two important parts: a pump radiation module and a probe radiation module. In most schemes of such analyzers, the same narrow-band laser is used, the radiation of which is divided with a coupler and used for probing and pumping. The pump radiation module modulates the radiation in intensity and then amplifies it with an erbium-doped fiber amplifier (EDFA). Thus, a series of rectangular pulses with a duration of the order of several nanoseconds is formed at the output. The time intervals between pulses appeared to be microseconds and were longer than a double passage of the pulse through the fiber sensor. Currently, a pulse length of 10 ns was used, which is equivalent to a 1 m sensor spatial resolution. Standard telecom single mode fiber was spliced with pigtails having E2000 type connectors with APC-polished ferrule by automatic fusion splicing instrument. The fiber was pre-wound on a standard Corning transport spool with minimal tension. Next, the fiber was placed in the Espec MC711 heat chamber. Fiber pinching and twisting, assumed to provoke the signal attenuation and the polarization state variations, were excluded. After that, the heat chamber was set to 25 °C. After heating the coil and the fiber for one hour, the measurement was started. For measurements, the option for the time-optimized maximum BGS search has been disabled. If it is enabled, the frequency resolution in the spectrum is not constant: it is rather low at the side parts of the spectrum and increased near the BGS maximum. This is not suitable for cross-correlation calculations since it requires a constant discrete shift between the functions.

Since the BFS for SMF-28 fiber under normal conditions is about 10.8 GHz, the frequency scanning range was set to 10,700 to 10,900 MHz with a scanning step of 0.5 MHz.

The measurement was carried out at various conditions of data averaging to simulate the signal-to-noise ratio deviations of the BGS. The SNR of the spectra was calculated using the LM algorithm to identify a useful signal, then subtract it from the original spectrum. Hence, obtaining the useful signal amplitude and noise amplitude, we can determine the SNR of the spectrum. The number of averaged BGS varied from 1 to 10,000. To obtain a reference spectrum, the number of averaged spectra was set to 100,000. The reference spectrum was used to calculate the error of the BFS extraction. An optical fiber about 1 km long was used in the experiment. The spatial resolution during measurements was 1 m, which ensured the acquisition of about 1000 spectra at the given settings.

Since the used commercial BOTDA setup could not acquire spectra with extremely low SNR (below 2 dB), the spectra for this range were obtained by adding the random value noise to the measured spectrum with the lowest measured SNR. The procedure was analogous to the one described above.

5. Results and Discussion

Figure 9 shows the results of BFS determination accuracy numerical simulation by three methods for the following parameters: scanning range 10,700–10,900 MHz, scan step 1 MHz (that is, the BGS contains 201 points). f_b, representing BFS, was chosen randomly in the range of 10,750–10,850 MHz. For each given SNR, the error was averaged as an absolute value over 100 BGS realizations.

The combined method began to give advantages only in the area of negative SNRs, where the errors are already so large that, in fact, none of the methods can be considered successful for the practice. In the region of positive SNRs, the best results were obtained using the standard LM algorithm.

Figure 9 shows that for all signal-to-noise ratios less than approximately −3 dB, the backward correlation methods become more efficient than LM. This is due to the fact that at sufficiently high noise levels, it is already difficult to call this curve Lorentzian, and therefore, the fitting of the Lorentzian function occurs with a rather high error. The BWC and MBWC methods do not require strict requirements for a certain signal shape; it is enough that the signal has a clearly localized smooth maximum. The presence of such a maximum in a signal of this type is proved by us in [23], and the smoothness of the line is provided by the mechanism of mutual cancellation of the noise of two signals when calculating the discrete correlation function. We reported on this in more detail in [24]. At the same time, it should be noted that the errors of all methods in this SNR range (<−3 dB) are already too large for any of them to be considered working.

As was previously noted in [24], the greatest gain from the use of BWC (and, consequently, MBWC) should be expected with a large number of points in the spectrum. With an increase in the BGS resolution (that is, a decrease in the frequency scanning step while maintaining the scan range), there is a standard gain in the accuracy of determining the BFS for all algorithms—the error decreases in proportion to the root of the scanning step [28]. However, BWC and MBWC do not work with the original BGS but with the correlogram, having higher SNR, as shown in Figure 2a and Figure 4a. Additionally, the SNR of the correlogram itself grows with increasing fineness of the BGS resolution because the contribution to the convolution from the products of useful signals is proportional to the number of points, while the non-desired contributions from the products of the useful signal with the side signal and from the product of two side signals are proportional only to the square root of the number of points.

Therefore, similar simulations were performed for spectra with a finer resolution: scanning range of 10,700–10,900 MHz, scanning step of 0.5 MHz (that is, the BGS now contains 401 points),and f_b is randomly selected in the range of 10,750–10,850 MHz. The BFS determination error averaged over 100 spectral realizations, which is shown in Figure 10.

In Figure 10, where the spectral resolution is doubled, the effect of error reduction in the area of −5 dB is not seen, which indirectly confirms the earlier hypothesis about the influence of low spectral resolution on the appearance of this effect. The trend and behavior of the curves are in agreement with our previous studies mentioned above. It is noteworthy that in all cases where a chain of two sequentially applied BWC-based methods was used, the gain in accuracy begins already in positive SNRs. This proves that the combined use of different BFS extraction methods enhances their strengths and offsets their weaknesses.

Currently, the correlation algorithms have started to be useful already in the area of positive SNRs, beginning from about 5 dB. Additionally, it is the combined method that demonstrates the best results, the gain of which, in terms of BFS extraction precision in comparison with LM, is about 0.5 MHz at SNR = 5 dB and increases with decreasing SNR (for example, at SNR = 0 dB, the gain is 5 MHz). In terms of more familiar values for sensor users, this is 0.5–5 K more accurate temperature measurement or 10–100 με more accurate strain measurement.

Figure 11 shows the results of processing real spectra obtained in the experiment. There is a good qualitative agreement with the simulation results from Figure 9, namely: starting from an SNR value of about 3 dB, the combined method gives the best results, and the smaller the SNR, the greater the gain. In our opinion, it is the slightly different noise distributions in the experiment and in the model that lead to a small quantitative difference in curve interception positions. Nevertheless, the correspondence between the experimental and simulated data should be recognized as satisfactory. The general trend—the gain of correlation methods and, in particular, the combined one, at low SNRs—is confirmed.

In general, the application of the MBWC method makes it possible to work with a correlogram instead of the original spectrum, where the SNR is higher, and all information about the position of the Brillouin scattering peak must be saved. Therefore, it is not necessary to take LM as a “second step”, as shown in Figure 11 where the results of spectra processing by double sequential application of MBWC are presented for comparison. Although this technique loses to the combination of MBWC + LM, it still outperforms the pure LM method.

The experimental results given in Figure 11 qualitatively confirm the simulation data. The behavior of the curves is practically the same as in Figure 9 and Figure 10. Of course, there are some differences: the intersection points of different curves slightly differ. Consequently, the boundaries of areas in which individual methods have advantages or disadvantages are slightly shifted. We believe that this difference is due to the presence in the experiment of those factors that were not taken into account in the simulation. Since we do not know the exact optical and electronic design of the commercial BOTDA we used, some individual noise components that affect the signal may not have been taken into account. This issue requires further investigation.

It should be noted that BFS extraction error is not the only parameter that has to be taken into consideration when choosing the spectra processing method. For the reader’s convenience, we have summarized the features of the different methods in the following

Table 1. Methods features.

Method/Parameter	BWC (Single)	MBWC (Single)	LM (Single)	MBWC + LM
Case 1 (201 point simulation) error at 0 dB SNR	24.9	6.12	4.48	6.46
Case 2 (401 point simulation) error at 0 dB SNR	24.7	8.27	9.34	4.12
Case 3 (experimental) error at 0 dB SNR	19.5	7.53	4.63	3.95
Efficiency at Low SNRs	Low *	High	Medium	Very High
Efficiency at High SNRs	Medium	Medium	High	Medium
Spectrum Non-symmetry influence	Low	Low	High	Low
Estimated computational costs	Low	Low	High	High
Compatibility with other techniques	No *	Yes	Yes	Yes

* True only for spectra with non-zero averaged side signal.

.

The table shows that when the signal level is equal to the noise level, that is, the signal-to-noise ratio is 0 dB, the combined MBWC + LM method shows the best performance. A slightly lower accuracy is shown by another combined method double MBWC. The advantage of MBWC + LM can be explained by the fact that after calculating the correlogram, which still has little noise, a fitting that eliminates its effect occurs. Double MBWC, as in the classical implementation, simply determines the frequency on the correlogram that has the maximum intensity.

6. Summary

This paper presents a modification of the BWC method named the MBWC, the essence of which is to correctly take into account the average level of the side signal. The effect of SNR and fineness of BGS resolution on the accuracy of the BFS extraction from the BGS by the MBWC method is studied. The possibilities of using MBWC in combination with other methods are also considered.

The applicability area of both MBWC and combined algorithms, based on it, should be the processing of spectra with low SNRs. Examples of such applications are long-distance sensing lines, and measurements under extreme conditions (for example, at elevated temperatures), where SMF is not applicable, and it is necessary to use special optical fibers with obviously high losses. In this case, the joint use of MBWC and LM can give a gain of up to several K or tens of microstrains.

It also should be noted that this work opens a wide front for further research. Uniformization of the input and output parameters of BGS processing methods allows changing the sequence of algorithms to achieve the optimal result both in terms of BFS extraction accuracy and data processing time. In addition, it is possible to include new modules in the calculation sequence. The operation of these modules can be based on completely different principles. Thus, some attention deserves image processing methods, which have recently been increasingly applied to data obtained from distributed fiber optic sensors [29,30,31,32,33]. Distributed sensors based on the principles of Brillouin reflectometry and Brillouin analysis in the time domain are no exception [34,35,36,37]. It is easy to see that the initial data array obtained using such devices is a surface, that is, a three-dimensional graph, where the axes represent the distance, frequency, and signal intensity. If the intensity of the signal is represented as a color, then these data will turn into a bitmap. Processing by such methods will allow one to operate simultaneously in the time and frequency domains. Since many image processing filters are non-linear, there will be an inevitable change in the shape of the spectrum, which will eventually inhibit the use of the Lorentzian function. At the stage after these transformations, the LCF method will not be effective, and the focus will be moved toward correlation methods, including the MBWC algorithm.

Another important area of research is the interaction of the described methods with artificial intelligence and machine learning algorithms [38,39,40,41,42]. We have already made several successful attempts to investigate this, which is presented in [23,43]. However, in both mentioned works, the processing was actually two-stage: first, the analytical algorithms functioned (independently), then the neural network was applied. We believe that it is necessary to add new algorithms at the first stage and place them in a sequence. This will make the work of the neural network more productive, reduce the training time, and increase the accuracy of BFS extraction.