WO2023182938A2

WO2023182938A2 - Method for brillouin optical time domain analysis (botda) based dynamic distributed strain sensing

Info

Publication number: WO2023182938A2
Application number: PCT/SG2023/050184
Authority: WO
Inventors: Ying Song; Juanjuan Hu; Hui Dong; Hailiang Zhang
Original assignee: Agency For Science, Technology And Research
Priority date: 2022-03-22
Filing date: 2023-03-21
Publication date: 2023-09-28
Also published as: WO2023182938A3

Abstract

The present disclosure relates to Brillouin Optical Time Domain Analysis (BOTDA) based method and system for improving dynamic strain measurements in distributed fiber optical sensor which involve application of a total variation (TV) filter on a raw image to produce a denoised image which includes a numeral solution. The present disclose also relates to ascertaining two or more evaluation parameters selected from the group consisting of image peak signal to noise ratio (PSNR), similarity structure index metrics (SSIM), and strain vector correlation coefficient (corr_coef) and mean square error (MSE), and also the evaluation performance indicator based on a combination thereof for comparing performance of various denoising approaches.

Description

METHOD FOR BRILLOUIN OPTICAL TIME DOMAIN ANALYSIS (BOTDA) BASED DYNAMIC DISTRIBUTED STRAIN SENSING

Field

[0001] Embodiments of the invention relate to distributed optical fiber sensing, and particularly relate to a denoising method and system for a Brillouin optical time domain analyzer (BOTDA).

Background

[0002] Distributed fiber optical sensing techniques have attracted many research interests due to low cost, capability in long range continuous sensing, chemical inertness, and immunity to external-electromagnetic interference in applications where the optical fiber simultaneously acts as the optical transmitting medium and millions of sensing points.

[0003] Existing distributed fiber optical sensing techniques utilises different physical principles, such as Brillouin scattering, Rayleigh scattering and Raman scattering. Among them, BOTDA sensing is the one of the most attractive techniques that provide distributed temperature and strain measurement with high accuracy and resolution over a long sensing fiber or distance, even in hazardous environment. Thus, BOTDA based distributed fiber optical sensing systems have vast potential applications in civil engineering applications, especially in monitoring the structural health of infrastructures, e.g., tunnels, pipeline, bridge, and in geophysical research.

[0004]A major challenge of distributed BOTDA strain sensing relates to its performance limitations due to signal-to-noise ratio (SNR) degradation of sensing signals. While up to thousands of traces averaging may be used to improve SNR of BOTDA, this approach is very time-consuming and unavailable for dynamic BOTDA strain sensing. To-date, other approaches have been proposed for improving SNR of BOTDA, including: (i) pump pulse coding, (ii) wide-bandwidth frequency modulation, and (iii) image processing techniques. In pump pulse coding technique, pump pulse is modulated to a coded sequence to improve SNR. Simplex coding, bipolar coding and loop coding may be used, however, extra time is required due to frequency scan needed for every coded sequence. Wide- bandwidth frequency modulation technique, such as frequency comb-based sweep-free scheme, slope-assisted scheme, and fast-frequency sweeping scheme, requires complex system, and hardware configuration and modification. Image processing techniques applying existed image denoising approaches have been proposed to cope with dynamic BOTDA measurement challenge without any system configuration and hardware complexity modification. Among them, wavelet denoising (WD) and non-local means filter (NLM) have been proposed to improve the SNR of BOTDA. Furthermore, block-matching 3D filter (BM3D) has been proposed to achieve further BOTDA performance improvement.

Summary

[0005] In a first aspect of the invention, a Brillouin Optical Time Domain Analysis (BOTDA) based method for improving dynamic strain measurements in distributed fiber optical sensor is provided. The method comprises: ascertaining at least one raw image based on a raw measurement being a Brillouin Gain Spectrum (BGS) measurement; and applying a total variation (TV) filter, including: ascertaining a first image processing operator by ascertaining negative partial derivatives of a predetermined previous second image processing operator; ascertaining an iterative denoised image by summing the first image processing operator and the raw image; ascertaining an image gradient of the iterative denoised image; ascertaining a current second image processing operator; ascertaining a current image gradient energy of the iterative denoised image; ascertaining an absolute value difference between a predetermined previous image gradient energy and the current image gradient energy; ascertaining whether the absolute value difference is less than a product of an initial image gradient energy of the raw image and a predetermined iterative parameter; and based on ascertaining the absolute value difference is less than the product, ascertaining the iterative denoised image as a denoised image which includes a numerical solution, wherein the current second image processing operator is ascertained based on a first quotient of a first dividend and a first divisor, wherein the first dividend includes a difference between the previous second image processing operator and a product of a predetermined image dimension related constant and the image gradient, wherein the first divisor is a summation of the following: one; and a product of a second quotient and a square root of a summation of square values of the image gradient, wherein a second dividend of the second quotient is the predetermined image dimension related constant and a second divisor of the second quotient is an update weight, wherein the current image gradient energy is ascertained based on a summation of the following: a square value of the first image processing operator; and a product of the update weight and a square root of a summation of square values of the image gradient.

[0006] In an embodiment of the first aspect, the method further comprises: based on ascertaining the difference is no less than the product: assigning the current second image processing operator as the predetermined previous second image processing operator; assigning the current image gradient energy as the predetermined previous image gradient energy; and repeating the steps of applying the total variation (TV) filter.

[0007] In an embodiment of the first aspect, the method further comprises: curve-fitting the denoised image to produce a curve-fitted denoised image; and based on the curve-fitted denoised image, ascertaining a strain vector for the raw measurement. [0008] In an embodiment of the first aspect, the method further comprises: ascertaining a reference measurement based on a plurality of Brillouin Gain Spectrum (BGS) measurements; ascertaining at least one reference image based on the reference measurement; and curve-fitting the reference image to produce a curve-fitted reference image; and based on the curve-fitted reference image, ascertaining a reference strain vector for the reference measurement.

[0009] In an embodiment of the first aspect, the method further comprises: based on the denoised strain vector and the reference strain vector, evaluating a correlation coefficient (corr_coef), and/or a mean square error (MSE) between the denoised strain vectors and the reference strain vector.

[0010] In an embodiment of the first aspect, the method further comprises: based on the denoised image and the reference image, evaluating an image peak signal to noise ratio (PSNR) and/or a structure similarity index matrix (SSIM) between the denoised image and the reference image.

[0011] In an embodiment of the first aspect, the method further comprises: based on the correlation coefficient (corr_coef), the mean square error (MSE), the image peak signal to noise ratio (PSNR) and/or the structure similarity index matrix (SSIM), deriving a performance indicator of the denoised image; and evaluating performance of the denoised image by comparing the performance indicator of the denoised image against another performance indicator of another denoised image.

[0012] In a second aspect of the invention, a Brillouin optical time domain analyzer (BOTDA) based apparatus for improving dynamic strain measurements in distributed fiber optical sensor is provided. The apparatus comprises: a memory device configured to store storing a plurality of measurements; and a computing processor communicably coupled to the memory device and configured to: ascertain at least one raw image based on a raw measurement being a Brillouin Gain Spectrum (BGS) measurement; and apply a total variation (TV) filter, including: ascertaining a first image processing operator by ascertaining negative partial derivatives of a predetermined previous second image processing operator; ascertain an iterative denoised image by summing the first image processing operator and the raw image; ascertain an image gradient of the iterative denoised image; ascertaining a current second image processing operator; ascertain a current image gradient energy of the iterative denoised image; ascertain an absolute value difference between a predetermined previous image gradient energy and the current image gradient; ascertain whether the absolute value difference is less than a product of an intial image gradient energy of the raw image and a predetermined iterative parameter; and based on ascertaining the absolute value difference is less than the product, ascertain the iterative denoised image as a denoised image which includes a numerical solution, wherein the current second image processing operator is ascertained based on a first quotient of a first dividend and a first divisor, wherein the first dividend includes a difference between the previous second image processing operator and a product of a predetermined image dimension related constant and the image gradient, wherein the first divisor is a summation of the following: one; and a product of a second quotient and a square root of a summation of square values of the image gradient, wherein a second dividend of the second quotient is the predetermined image dimension related constant and a second divisor of the second quotient is an update weight, wherein the current image gradient energy is ascertained based on a summation of the following: a square value of the first image processing operator; and a product of the update weight and a square root of a summation of square values of the image gradient.

[0013] In an embodiment of the second aspect, the computing processor is further configured to: based on ascertaining the difference is no less than the product: assign the current second image processing operator as the predetermined previous second image processing operator; assign the current image gradient energy as the predetermined previous image gradient energy; and repeat the steps of applying the total variation (TV) filter.

[0014] In an embodiment of the second aspect, the computing processor is further configured to: curve-fit the denoised image to produce a curve-fitted denoised image; and based on the curve-fitted denoised image, ascertain a strain vector for the raw measurement.

[0015] In an embodiment of the second aspect, the computing processor is further configured to: ascertain a reference measurement based on a plurality of Brillouin Gain Spectrum (BGS) measurements; ascertain at least one reference image based on the reference measurement; and curve-fit the reference image to produce a curve-fitted reference image; and based on the curve-fitted reference image, ascertain a reference strain vector for the reference measurement.

[0016] In an embodiment of the second aspect, the apparatus further comprises: based on the denoised strain vector and the reference strain vector, evaluating a correlation coefficient (corr_coef), and/or a mean square error (MSE) between the denoised strain vectors and the reference strain vector. [0017] In an embodiment of the second aspect, the apparatus further comprises: based on the denoised image and the reference image, evaluating an image peak signal to noise ratio (PSNR) and/or a structure similarity index matrix (SSIM) between the denoised image and the reference image.

[0018] In an embodiment of the second aspect, the apparatus further comprises: based on the correlation coefficient (corr_coef), the mean square error (MSE), the image peak signal to noise ratio (PSNR) and/or the structure similarity index matrix (SSIM), deriving a performance indicator of the denoised image; and evaluating performance indicator of the denoised image by comparing the performance indicator of the denoised image against another performance indicator of another denoised image.

Brief Description of Drawings

[0019] Figure 1 is a block diagram showing an arrangement for BOTDA strain measurement.

[0020] Figure 2(a) is a block diagram showing an image denoising method for one-shot (single pump pulse) dynamic BOTDA strain measurement.

[0021] Figure 2(b) is a flow sequence for implementing block 23 of Figure 2(b).

[0022] Figure 3(a) shows a ground-truth figure and Figure 3(b) shows the 75th column in the ground-truth figure of Figure 3(a).

[0023] Figure 4 shows an example of strain calculation.

[0024] Figure 5 is a schematic representation of an experimental setup.

[0025] Figures 6(a) to 6(f) show 1 -shot dynamic measurement with multiple image denoising approaches.

[0026] Figures 7(a) to 7(f) show a Lorentzian curve fitting example for a strain located at distance sequence 16 for 1 -shot dynamic measurement with multiple image denoising approaches, respectively.

[0027] Figures 8(a) to 8(d) show the 1 -shot dynamic strain measurement vector obtained after Lorentzian curve fitting and peak detection. [0028] Figures 9(a) to 9(d) show the corresponding zoom-in views for strain sequences 1 to 40.

[0029] Figures 10(a) to 10(d) show violin plots corresponding to four evaluation parameters Corr-coef, MSE, PSNR, and SSIM, respectively, for 1 -shot dynamic strain measurement where the centre dash is the mean.

Detailed Description

[0030] In the following description, numerous specific details are set forth in order to provide a thorough understanding of various illustrative embodiments of the invention. It will be understood, however, to one skilled in the art, that embodiments of the invention may be practiced without some or all of these specific details. In other instances, well known process operations have not been described in detail in order not to unnecessarily obscure pertinent aspects of embodiments being described. In the drawings, like reference numerals refer to same or similar functionalities or features throughout the several views.

[0031] Embodiments described in the context of one of the methods or devices are analogously valid for the other methods or devices. Similarly, embodiments described in the context of a method are analogously valid for a device, and vice versa.

[0032] Features that are described in the context of an embodiment may correspondingly be applicable to the same or similar features in the other embodiments. Features that are described in the context of an embodiment may correspondingly be applicable to the other embodiments, even if not explicitly described in these other embodiments. Furthermore, additions and/or combinations and/or alternatives as described for a feature in the context of an embodiment may correspondingly be applicable to the same or similar feature in the other embodiments.

[0033] In the context of various embodiments, including examples and claims, the articles “a”, “an” and “the” as used with regard to a feature or element include a reference to one or more of the features or elements. The terms “comprising,” “including,” and “having” are intended to be open-ended and mean that there may be additional features or elements other than the listed ones. Identifiers such as “first”, “second” and “third” are used merely as labels, and are not intended to impose numerical requirements on their objects, nor construed in a manner imposing any relative position or time sequence between limitations. The term “and/or” includes any and all combinations of one or more of the associated listed items.

[0034] Figure 1 is a block diagram showing an arrangement for BOTDA strain measurement. The arrangement includes a sensing fiber under test (FUT) which is to be arranged along the asset or structure to be monitored, and a BOTDA interrogator equipment communicably coupled to the FUT.

[0035]A pulsed pump light is injected at a first end of the sensing FUT and propagates along the FUT, while at the opposite end of the FUT, a counter propagating continuous wave (CW) probe light having a specific wavelength is continuously injected. Optical power will transfer between the pump light and the probe light at all positions along the FUT if their frequency difference is close to the local BFS. The BFS has a linear dependance on temperature and tensile strain variation. Thus, by deriving all the local BFS, the temperature or/and strain variation along the FUT can be monitored. The local BFS can be obtained by BOTDA processing through fitting the measured local BGS with Lorentzian or Gaussian profile.

[0036] Embodiments of the invention disclose an image denoising method for improving dynamic distributed strain measurement. The method includes total variation (TV) regularization-based image denoising and noise estimation approach which may be implemented in a BOTDA interrogator. The method may further include ascertaining two or more evaluation parameters selected from the group consisting of image peak signal to noise ratio (PSNR), similarity structure index metrics (SSIM), and strain vector correlation coefficient (corr_coef) and mean square error (MSE) against a reference measurement. [0037] Figure 2(a) is a block diagram showing an image denoising method for one-shot (single pump pulse) dynamic BOTDA strain measurement.

[0038] In block 1 1 , a reference measurement is ascertained. This reference measurement may be ascertained based on a plurality of ground-truth measurements, e.g., static Brillouin gain spectrum (BGS) measurements acquired using conventional BOTDA processing, at multiple positions along the length of the FUT. The reference measurement may be ascertained as an average of the BGS measurements. Count of the BGS measurements may be 5000, or more, or less.

[0039] In block 12, based on the reference measurement, at least one reference image, e.g. a plurality of images, is ascertained. The reference image provides the frequency scan difference between the pump light and the probe light at multiple optical positions along the length of the FUT, e.g. vertical axis of the reference images represents BGS while horizontal axis represents optical fiber distance.

[0040] In block 13, the reference image is curve-fitted to produce a curve-fitted reference image. Lorentzian curve or Gaussian distribution curve may be used. More details on application of Lorentzian curve are provided in later paragraphs.

[0041] In block 14, based on the curve-fitted reference image, a reference strain vector y is ascertained. More details on strain vector are provided in later paragraphs.

[0042] In block 21 , a first raw measurement, being a BGS measurement, is ascertained. The raw measurement may be obtained from one shot measurement or an average measurement, e.g. average of two or five or other number of measurements, from BOTDA interrogator.

[0043] In block 22, based on the raw measurement, at least one raw image, e.g. a plurality of raw images, are ascertained. Raw measurement from the BOTDA interrogator includes two-dimensional (2D) data for which x-axis is distance along the optical cable and y-axis is the frequency which is associated to the distributed strain measured. The raw measurement, e.g. a long 2D image, may be segmented along the distance into a plurality of raw images, e.g. short 2D images.

[0044] In block 23, a total variation (TV) filter is applied to the at least one raw image to ascertain at least one denoised image.

[0045]Generally, image denoising is a procedure to reconstruct a true image from a noisy image. The problem of image denoising can be modeled as follows: z = x + n (Equation 1 ) where z is the observed noisy image, x is the unknown clean image, and n represents additive white Gaussian noise with standard deviation on. The purpose of noise reduction is to decrease the noise while minimizing the loss of original features. In other words, image denoising is a procedure for reducing noise but keeping flat areas smooth, protecting edges from blur, preserving textures characteristics and not generating new artifacts. Owing to solve the clean image x from the Equation 1 may be an ill-posed problem, a unique solution may not be ascertained from an image model with noise. To obtain a good estimation of denoised image x , multiple image denoising methods including spatial domain and transfer domain approaches may be employed.

[0046]The TV denoising used in Figure 2(a) is based on a spatial domain filtering approach utilizing image priors and minimizing an energy function E to calculate the denoised image x modeled as: (Equation 2)

[0047]The motivation for variational denoising methods of Equation 2 is maximum a posterior (MAP) probability estimate. From a Bayesian perspective, the MAP probability estimate P of x is (Equation 3)

[0048] By equivalently formulating the Equation 3 and reformulating the Equation 3 under Gaussian white noise, using the statistical fact that natural images are locally smooth, and the pixel intensity gradually varies in most regions, the TV denoising model x can be presented as: x = arg min- \\z — x\\ + A||V_X||₁ (Equation 4) x 2 where A is the regularization parameter and the V_x is the gradient of x.

[0049] In view of the, block 23 may be implemented based on a flow sequence of Figure 2(b).

[0050] In block 231 of Figure 2(b), the following parameters may be ascertained and/or inputted into a computing processor for performing the remaining blocks:

• At least one original normalized raw image /mo, e.g. noisy image z;

• Initial iterative time value n which may be set at one, i.e. n = 1 ;

• Update weight co which may be ascertained based on a predetermined

/c regularization parameter A, where A = — and k may be a predetermined constant based on image noise variance;

• Iterative parameter £ which may be used to stop iteration when prescribed conditions are met and may be predetermined based on experimental result and requirement;

• Maximum iterative time value N;

• Image dimension related constant C.

[0051] In block 232, at least some of the following parameters may be initialized or their initial values ascertained:

• Initial image gradient Go = diff (imo), where diff refers to first-order partial derivatives; in general, image gradient G may be ascertained based on first-order partial derivatives of the original normalized raw image imo with respect to its horizontal axis x and vertical axis y; image gradient G may refer to change in intensity or brightness of the image; • Initial first image processing operator Fo is set at zero, i.e. Fo = 0;

• Initial second image processing operator

HQ = C and co are predetermined constants;

• Initial image gradient energy

E_o

, where Eo may alternatively be referred to as initial E.

The first image processing operator and the second image processing operator /-/will be further explained in block 234 and block 237 respectively.

[0052] In block 233, a current iterative time value n is ascertained whether it has exceeded a predetermined maximum time value N. If the current iterative time value has not exceeded the predetermined maximum time value N, block 233 proceeds to block 234. If the iterative time value has reached or exceeded the predetermined maximum time value N, block 233 proceeds to block 240.

[0053] In block 234, a first image processing operator is ascertained based on a negative divergence, e.g. partial derivatives, of a previous second image processing operator along both vertical and horizontal axes. This operation may be represented by:

(Equation 5)

is the partial derivative along the horizontal axis the partial

derivative along the vertical axis y of the image. The previous second image processing operator is predetermined. For the first iteration n=1 , the previous second image processing operator is Ho as ascertained in block 232; for subsequent iteration, the previous second image processing operator is ascertained in previous iteration.

[0054] In block 235, based on the first image processing operator F_n ascertained in block 234 and the original normalized raw image imo, an iterative denoised image is ascertained by summing the original normalized raw image imo and the first image processing operator F_n. This operation may be represented by: iterative denoised image = imo + F_n. This iterative denoised image may be subsequently provided as the denoised image output in block 240 subject to block 239.

[0055] In block 236, an image gradient of the iterative denoised image G_n is ascertained. Known techniques may be used to ascertain image gradient. For example, an image gradient may be ascertained based on first-order partial derivatives of the iterative denoised image with respect to its horizontal axis x and vertical axis y.

[0056] In block 237, a current second image processing operator H_n is ascertained, e.g. updated, based on the previous second image processing operator H_n-i, image dimension related constant C, image gradient G_n ascertained in block 236, and update weight co. In particular, the current second image processing operator H_n may be ascertained based on a first quotient of a first dividend and a first divisor, wherein the first dividend includes a difference between the previous second image processing operator H_n-i and a product of a predetermined image dimension related constant C and the image gradient Gn, wherein the first divisor is a summation of the following: one; and a product of a second quotient and a square root of a summation of square values of the image gradient G_n, wherein a second dividend of the second quotient is the predetermined image dimension related constant C and a second divisor of the second quotient is an update weight co. This operation may be represented by:

(Equation 6)

[0057] In block 238, a current image gradient energy E_n is ascertained based on a summation of the following: a square value of the first image processing operator F_n and a product of the update weight co and a square root of a summation of square values of the image gradient G_n. This operation may be represented by: E_n ⁼ -n.² + ⁶⁰ J sum G ²

(Equation 7) [0058] In block 239, an absolute value (modulus) of a difference between a previous image gradient energy E_n-i and a current image gradient energy E_n is ascertained. It is also ascertained whether the absolute value of the absolute value difference is less than a product of a predetermined iterative parameter E and initial image gradient energy of the raw image (initial E or Eo). This operation may be represented by:

\E_{n- 1} — E_n \ < £ (initial E) or

(Equation 8)

For the first iteration n=1 , the previous image gradient energy is Eo as ascertained in block 232; for subsequent iteration, the previous image gradient energy is E_n-i ascertained in previous iteration.

[0059] If the above-mentioned absolute value difference is less than the above- mentioned product, block 239 proceeds to block 240 in which the iterative denoised image is ascertained as the output or final denoised image which includes a numerical solution.

[0060] If the above-mentioned absolute value difference is no less than, e.g. greater than or equal to, the above-mentioned product, block 239 proceeds to block 241 where the iterative time value is incremented, e.g. by one, and thereafter steps of block 234 to block 239 are repeated or iterated. Prior to performing the iteration, the current second image processing operator may be assigned as the predetermined previous second image processing operator; the current image gradient energy may be assigned as the predetermined previous image gradient energy. [0061]Block 23 proceeds to block 13 in which the denoised image is curve-fitted to produce a curve-fitted denoised image. Lorentzian curve or Gaussian distribution curve may be used. Lorentzian function can be described as: (Equation 9)

where A represents amplitude, w is the full width at half maximum (FWHM) of the function peak. x₀ represents the center of Lorentzian Function. Hence, Lorentzian curve fitting (LCF) is configured to ascertain the center frequency x₀.

[0062] LCF is performed separately for each column of the images of block 12 and block 26. The outputs are frequency vectors for all the columns, respectively.

[0063] For example, Figure 3(a) shows a ground-truth figure and Figure 3(b) shows the 75th column in the ground-truth figure of Figure 3(a). The asterisk represents the collected data, the curve represents the LCF results. The frequency was scanned from 10650 MHz to 11349 MHz with a scanning step of 3 MHz. The fitted x₀ is about 10847.5 MHz. It should be noted that 0 of y axis in Figure 3(a) correspond to 10650MHz. Normally, without applying strain on the optical fiber, the Brillouin frequency of a conventional single mode fiber is about 10.85GHz.

[0064] In block 14, based on the curve-fitted denoised image, a denoised strain vector y_t is ascertained.

[0065] It is understood that if strain and/or temperature varies, the Brillouin frequency of the fiber will shift. By monitoring the frequency shift, the strain and/or temperature can be ascertained. The strain and temperature coefficients are known parameters. Different fibers may have slight different strain and temperature coefficients, they can be calibrated. Strain coefficient may be about 19 pe/MHz, temperature coefficient may be about 1 °C/MHz.

[0066]An example of strain calculation is described with reference to Figure 4. It corresponds to one column of the denoised image, which only represents one local position of the optical fiber. The original fitted center is 10849.8MHz, but after strain is applied to the optical fiber, the fitted center is 11017.9MHz. Then, when the frequency shift is 168.1 MHz, strain is a product of frequency shift, e.g. 168.1 MHz, and a predetermined strain coefficient, e.g. 19 pe/MHz. Accordingly, the strain is 168.1 MHz x 19 ps/MHz = 3193.9 pe.

[0067] Blocks 15 and 16 include evaluating measurement performance based on strain vector-based parameters and image-based parameters, respectively.

[0068] Block 15 includes evaluating measurement performance based on strain vectorbased parameters derived from the denoised strain vector and the reference strain vector. Particularly, the vector-based parameters include correlation coefficient (corr_coef), and mean square error (MSE) between the denoised strain vectors and the reference strain vector.

[0069] Given a reference (ground-truth) strain vector y_f, the quantitative measurements for denoised (measured) strain vector y_t corr_coef and MSE with ground-truth strain vector y_t are obtained by: (Equation 10) (Equation 11 )

where u_y and u_y are mean of y_t and y_t , respectively.

[0070] Block 16 includes evaluating measurement performance based on image-based parameters derived from the denoised image and the reference image. Particularly, the image-based parameters include image peak signal to noise ratio (PSNR) and structure similarity index matrix (SSIM) between denoised (measured) image and reference (ground-truth) image.

[0071]Given a reference image x, the quantitative measurements PSNR of a denoised image x is obtained by:

Y1 (Equation 12) (Equation 13)

where ._x , a_x , a_x are the means and variances of x and x, respectively, a_xx is the covariance between x and x , and Ci and C2 are constant values used to avoid instability.

[0072]To illustrate the performance improvement of the above-described methods over conventional methods, an experimental setup and experimental results are provided.

[0073] Figure 5 is a schematic representation of an experimental setup which uses a conventional BOTDA setup to measure BGS distribution along a 100 m FUT. Spatial step of the data and the pulse width were 0.5 m and 10 ns, respectively. The BGS was obtained by scanning the frequency difference with a step of 3 MHz from 10650 MHz to 11349 MHz. The BGS with 5000 times average was collected as reference (ground-truth) measurement. For dynamic measurement, raw measurement with only 1 time measure (1 -shot) was collected for investigating the above-described methods. The sensing fiber is a polyimide coated single mode fiber (SMF). Ends of one segment of the fiber were mounted on two pieces of Lexan Polycarbonate (PC) sheets by using epoxy. The PC sheets were fixed on two stages with a distance about 2.12 m. Static strain was applied to the segment of the fiber between the two stages by fixing one stage and moving the other stage axially. The applied static strain was about 2380 pe for the experiments. The data were collected by using the default settings of a commercial BOTDA interrogator for the conventional SMF. The strain coefficient and temperature coefficient were set to be 18.915 pe/MHz and 0.9765°C/MHz, respectively. A thermocouple temperature sensor was used to roughly compensate the temperature perturbation (resolution is 0.1 °C).

[0074] In the experiment, the same set of raw measurements are denoised using different image denoising approaches, namely, WD, NLM, BM3D and methods of the invention embodiments, followed by curve-fitting using Lorentzian curve to derive the central/peak frequency which corresponds the strain value. Furthermore, performance evaluation of dynamic strain measurement from both denoising images and derived strain vectors are conducted using the above-described four quantitative measurement parameters.

[0075]As an example, Figures 6(a) to 6(f) show 1 -shot dynamic measurement with multiple image denoising approaches, wherein vertical axis is frequency and horizontal axis is distance. Figures 6(a) to 6(f) correspond to images obtained from reference (ground-truth), 1 -shot raw measurement, image denoising approaches WD, NLM, BM3D and the methods of the invention embodiments, respectively. It was observed that after image denoising, there is obvious SNR improvement comparing a denoised image of the invention embodiments as shown in Figures 6(f) with 1 -shot raw measurement image as shown in Figure 6(b). Similarly, a difference was observed comparing a denoised image of the invention embodiments as shown in Figures 6(f) with the reference (ground truth) image shown in Figure 6(a). Furthermore, more distortion is observed from denoised image using WD approach.

[0076] As an example, Figures 7(a) to 7(f) show a Lorentzian curve fitting example for a strain located at distance sequence 16 for 1 -shot dynamic measurement with multiple image denoising approaches, respectively. The horizontal axis represents frequency corresponding to strain magnitude while the vertical axis represents amplitude. It was observed that in Figure 7(f) relating to the denoised image of the invention embodiments, the fitted curve is closer to the reference (ground-truth) fitted curve which leads to the closer peak location obtained. Based on Figures 7(a) to 7(f), in relation to closeness of to the reference (ground-truth) fitted curve, it was observed that order of performance of denoising approaches is in the order of invention embodiments, BM3D, WD, NLM (from best to worst).

[0077]As an example, Figures 8(a) to 8(d) show the 1 -shot dynamic strain measurement vector obtained after Lorentzian curve fitting and peak detection while Figures 9(a) to 9(d) show the corresponding zoom-in views for strain sequences 1 to 40. It is observed that it is more difficult to use quality and visual observation to compare the denoising effect based on the denoising images shown in Figures 7(a) to 7(f) as it is indirectly related to the measured distributed strain vector. However, from Figures 9(a) to 9(d), it is evident that the BM3D approach and invention embodiments provide better dynamic strain measurement as the obtained dynamic strain measurements are closer to the reference (ground-truth).

[0078] To quantitatively compare evaluation parameters from the above experimental results, Table 1 shows PSNR and SSIM ascertained from denoised images, and Corr_coef and MSE between curve-fitted strain vectors of the dynamic raw measurement and reference (ground-truth) measurement. As aforementioned, WD, NLM, BM3D and invention embodiments approaches were used for the collected 1 -shot raw measurement result based on experiment setup shown in Figure 3. The mean of up to 9 times measurement were used to ensure reliable result.

[0079] When using PSNR (mean) as evaluation parameter, which is basis of conventionally used performance evaluation parameter, the measurement results in the order sequence from best to worst (higher value is better) is: invention embodiments (18.62dB), BM3D (16.78dB), NLM (16.56dB), WD (13.60dB) where the raw measurement result is 9.61 dB. The SSIM (mean) provides the same order sequence as the PSNR. However, it is observed that when MSE was used to evaluate the performance, the order sequence from best to worst (lower value is better) is: BM3D (183.54), invention embodiment (197.99), NLM (282.71 ), WD (452.33) which is different from PSNR and SSIM. In addition, it was observed that Corr_coef (near value of 1 is better) provides an order sequence from best to worst: invention embodiments (0.9611 ), BM3D (0.9561 ), WD (0.9490), NLM (0.9370) which is different from PSNR, SSIM and MSE.

Table 1 : Quantitative accuracy evaluation for 1 -shot dynamic strain measurement

[0080] As the four evaluation parameters provide different evaluation result, a unique evaluation parameter for BOTDA strain measurement performance may be formulated as:

Strain measurement index, Ml = SI (PSNR) + SI (SSIM) + SI (MSE) + SI (Corr_coef) where SI is the sequence index in the range of 1 to 4 if five different methods are being compared and a SI of 1 represents best performance while SI of 4 represents worst performance. For the 1 -shot example, the best performance is provided by invention embodiments which Ml (invention embodiments) is 1 +1 +2+1 =5. Ml (BM3D)=7 is ranked second, Ml (NI_M)=13 is ranked third, while Ml (WD)=15 is ranked fourth. Accordingly, performance indicates of denoised images which may be ascertained from various approaches may be derived. Thereafter, performance of the denoised images and/or their respective approaches may be evaluated by comparing the respective performance indicators.

[0081] Aside from accuracy of measurement, measurement uncertainty (robustness) may be another indicator for performance evaluation. Figures 10(a) to 10(d) show violin plots corresponding to four evaluation parameters Corr-coef, MSE, PSNR, and SSIM, respectively, for 1 -shot dynamic strain measurement where the centre dash is the mean. It was observed for PSNR, although invention embodiments result in the highest PSNR, the range of PSNR is almost maintained the same as with other approaches. For SSIM, invention embodiments achieve a more narrow range than BM3D and NLM. For the corr_coef and MSE, although the range of invention embodiments has a slightly wider range as compared with other approaches, the distance is small.

[0082] Table 2 lists the quantitative ranges corresponding to Figures 10(a) to 10(d). In general, the proposed scheme is with best robustness based on the comprehensive evaluation.

[0083] Table 2: Quantitative robustness evaluation for 1 -shot dynamic strain measurement.

[0084] In an aspect of the disclosure, a BOTDA based apparatus, e.g. interrogator, comprises computing processor(s), memory storage device(s) configured to store computer-executable code for performing operations in accordance with the abovedescribed methods, e.g. Figure 2(a) and/or 2(b), display device(s) for displaying outputs generated by performing the operations, communication device(s) configured to receive data and/or signals produced by the FUT and/or transmit data and/or signals produced by the computing processor(s). The computing processor(s) may be communicably coupled to the above-mentioned memory storage device(s), displace device(s), communication device(s) and/or other electronic devices.

[0085] In an aspect of the disclosure, a BOTDA based system for improving dynamic strain measurements in distributed fiber optical sensor comprises a BOTDA interrogator and an optical fiber communicably coupled to the BOTDA interrogator. [0086] In an aspect of the disclosure, a non-transitory computer-readable medium having computer-readable code executable by at least one processor is provided to perform the method/steps as described in the foregoing, e.g. Figures 2(a) and/or 2(b).

[0087] Embodiments of the invention provide at least the following advantages:

- With total variation (TV) regularization-based image denoising and noise estimation approach, dynamic distributed strain measurement is improved.

- Multiple evaluation parameters, e.g. image peak signal to noise ratio (PSNR), similarity structure index metrics (SSIM), and strain vector correlation coefficient (corr_coef) and mean square error (MSE), which are based on image and strain vector, are provided for more comprehensive performance evaluation metrics.

[0088] Other embodiments will be apparent to those skilled in the art from consideration of the specification and practice of the disclosure. Furthermore, certain terminology has been used for the purposes of descriptive clarity, and not to limit the disclosed embodiments. The embodiments and features described above should be considered exemplary.

Claims

1. A Brillouin Optical Time Domain Analysis (BOTDA) based method for improving dynamic strain measurements in distributed fiber optical sensor, the method comprising: ascertaining at least one raw image based on a raw measurement being a Brillouin Gain Spectrum (BGS) measurement; and applying a total variation (TV) filter, including: ascertaining a first image processing operator by ascertaining negative partial derivatives of a predetermined previous second image processing operator; ascertaining an iterative denoised image by summing the first image processing operator and the raw image; ascertaining an image gradient of the iterative denoised image; ascertaining a current second image processing operator; ascertaining a current image gradient energy of the iterative denoised image; ascertaining an absolute value difference between a predetermined previous image gradient energy and the current image gradient energy; ascertaining whether the absolute value difference is less than a product of an initial image gradient energy of the raw image and a predetermined iterative parameter; and based on ascertaining the absolute value difference is less than the product, ascertaining the iterative denoised image as a denoised image which includes a numerical solution, wherein the current second image processing operator is ascertained based on a first quotient of a first dividend and a first divisor, wherein the first dividend includes a difference between the previous second image processing operator and a product of a predetermined image dimension related constant and the image gradient, wherein the first divisor is a summation of the following: one; and a product of a second quotient and a square root of a summation of square values of the image gradient, wherein a second dividend of the second quotient is the predetermined image dimension related constant and a second divisor of the second quotient is an update weight, wherein the current image gradient energy is ascertained based on a summation of the following: a square value of the first image processing operator; and a product of the update weight and a square root of a summation of square values of the image gradient.

2. The method according to claim 1 , further comprising: based on ascertaining the difference is no less than the product: assigning the current second image processing operator as the predetermined previous second image processing operator; assigning the current image gradient energy as the predetermined previous image gradient energy; and repeating the steps of applying the total variation (TV) filter.

3. The method according to claim 1 or claim 2, further comprising: curve-fitting the denoised image to produce a curve-fitted denoised image; and based on the curve-fitted denoised image, ascertaining a strain vector for the raw measurement.

4. The method according to claim 3, further comprising: ascertaining a reference measurement based on a plurality of Brillouin Gain Spectrum (BGS) measurements; ascertaining at least one reference image based on the reference measurement; and curve-fitting the reference image to produce a curve-fitted reference image; and based on the curve-fitted reference image, ascertaining a reference strain vector for the reference measurement.

5. The method according to claim 4, further comprising: based on the denoised strain vector and the reference strain vector, evaluating a correlation coefficient (corr_coef), and/or a mean square error (MSE) between the denoised strain vectors and the reference strain vector.

6. The method according to any one of claim 4 to claim 5, further comprising: based on the denoised image and the reference image, evaluating an image peak signal to noise ratio (PSNR) and/or a structure similarity index matrix (SSIM) between the denoised image and the reference image.

7. The method according to claim 6, further comprising: based on the correlation coefficient (corr_coef), the mean square error (MSE), the image peak signal to noise ratio (PSNR) and/or the structure similarity index matrix (SSIM), deriving a performance indicator of the denoised image; and evaluating performance of the denoised image by comparing the performance indicator of the denoised image against another performance indicator of another denoised image.

8. A Brillouin optical time domain analyzer (BOTDA) based apparatus for improving dynamic strain measurements in distributed fiber optical sensor, the apparatus comprising: a memory device configured to store storing a plurality of measurements; and a computing processor communicably coupled to the memory device and configured to: ascertain at least one raw image based on a raw measurement being a Brillouin Gain Spectrum (BGS) measurement; and apply a total variation (TV) filter, including: ascertaining a first image processing operator by ascertaining negative partial derivatives of a predetermined previous second image processing operator; ascertain an iterative denoised image by summing the first image processing operator and the raw image; ascertain an image gradient of the iterative denoised image; ascertaining a current second image processing operator; ascertain a current image gradient energy of the iterative denoised image; ascertain an absolute value difference between a predetermined previous image gradient energy and the current image gradient; ascertain whether the absolute value difference is less than a product of an intial image gradient energy of the raw image and a predetermined iterative parameter; and based on ascertaining the absolute value difference is less than the product, ascertain the iterative denoised image as a denoised image which includes a numerical solution, wherein the current second image processing operator is ascertained based on a first quotient of a first dividend and a first divisor, wherein the first dividend includes a difference between the previous second image processing operator and a product of a predetermined image dimension related constant and the image gradient, wherein the first divisor is a summation of the following: one; and a product of a second quotient and a square root of a summation of square values of the image gradient, wherein a second dividend of the second quotient is the predetermined image dimension related constant and a second divisor of the second quotient is an update weight, wherein the current image gradient energy is ascertained based on a summation of the following: a square value of the first image processing operator; and a product of the update weight and a square root of a summation of square values of the image gradient.

9. The apparatus according to claim 8, wherein the computing processor is further configured to: based on ascertaining the difference is no less than the product: assign the current second image processing operator as the predetermined previous second image processing operator; assign the current image gradient energy as the predetermined previous image gradient energy; and repeat the steps of applying the total variation (TV) filter.

10. The apparatus according to claim 8 or claim 9, wherein the computing processor is further configured to: curve-fit the denoised image to produce a curve-fitted denoised image; and based on the curve-fitted denoised image, ascertain a strain vector for the raw measurement.

11. The apparatus according to claim 10, wherein the computing processor is further configured to: ascertain a reference measurement based on a plurality of Brillouin Gain Spectrum (BGS) measurements; ascertain at least one reference image based on the reference measurement; and curve-fit the reference image to produce a curve-fitted reference image; and based on the curve-fitted reference image, ascertain a reference strain vector for the reference measurement.

12. The apparatus according to claim 1 1 , further comprising: based on the denoised strain vector and the reference strain vector, evaluating a correlation coefficient (corr_coef), and/or a mean square error (MSE) between the denoised strain vectors and the reference strain vector.

13. The apparatus according to any one of claim 11 to claim 12, further comprising: based on the denoised image and the reference image, evaluating an image peak signal to noise ratio (PSNR) and/or a structure similarity index matrix (SSIM) between the denoised image and the reference image.

14. The apparatus according to claim 13, further comprising: based on the correlation coefficient (corr_coef), the mean square error (MSE), the image peak signal to noise ratio (PSNR) and/or the structure similarity index matrix (SSIM), deriving a performance indicator of the denoised image; and evaluating performance indicator of the denoised image by comparing the performance indicator of the denoised image against another performance indicator of another denoised image.

15. A non-transitory computer-readable medium having computer-readable code executable by at least one processor to perform the method according to any one of claims 1 to 7.