CN113343173B - Brillouin frequency shift extraction method - Google Patents

Brillouin frequency shift extraction method Download PDF

Info

Publication number
CN113343173B
CN113343173B CN202110608992.6A CN202110608992A CN113343173B CN 113343173 B CN113343173 B CN 113343173B CN 202110608992 A CN202110608992 A CN 202110608992A CN 113343173 B CN113343173 B CN 113343173B
Authority
CN
China
Prior art keywords
cross
brillouin
correlation
frequency
frequency shift
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110608992.6A
Other languages
Chinese (zh)
Other versions
CN113343173A (en
Inventor
王曰海
苏梁灏
陈蓓
刘笑之
杨建义
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN202110608992.6A priority Critical patent/CN113343173B/en
Publication of CN113343173A publication Critical patent/CN113343173A/en
Application granted granted Critical
Publication of CN113343173B publication Critical patent/CN113343173B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/15Correlation function computation including computation of convolution operations
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/18Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Monitoring And Testing Of Transmission In General (AREA)
  • Complex Calculations (AREA)

Abstract

The invention discloses a brillouin frequency shift extraction method, which comprises the following steps: constructing a cross-correlation feature matrix according to the Brillouin gain spectrum, and constructing a corresponding Brillouin frequency shift dictionary set; the method comprises the steps of carrying out Brillouin frequency shift retrieval, carrying out dot product on a Brillouin gain spectrum and a cross-correlation feature matrix to obtain a cross-correlation degree matrix, and retrieving a value of Brillouin frequency shift from a Brillouin frequency shift dictionary set according to a position with the largest cross-correlation degree in the cross-correlation degree matrix; the cross-correlation method does not need to interpolate the Brillouin gain spectrum, can process a plurality of spectral lines at the same time, has high processing speed, and provides a feasible way for realizing the rapid measurement of optical fiber sensing.

Description

Brillouin frequency shift extraction method
Technical Field
The invention belongs to the technical field of distributed optical fiber sensing, and particularly relates to a Brillouin frequency shift extraction method.
Background
The Brillouin optical time domain analyzer has the capability of simultaneously monitoring temperature and strain, and can be applied to scenes such as fire alarm, oil and gas pipeline leakage detection, power cables and the like which need to monitor temperature and stress.
The working principle of the Brillouin Optical Time Domain Analyzer (BOTDA) is that pumping light is injected into two ends of a sensing optical fiber at one end, continuous detection light is injected into the other end, and when two beams of light meet in the optical fiber, if the difference of frequencies of the two beams of light is in a brillouin spectrum, brillouin scattering is generated, and energy transfer from the pumping light to the detection light occurs. The Brillouin Gain Spectrum (BGS) at each position in the optical fiber can be obtained by continuously adjusting the frequency of the probe light. The brillouin gain spectrum is a lorentz curve shape:
wherein g 0 For peak gain factor, v B For Brillouin Frequency Shift (BFS), Δv is the full-width half-maximum of the brillouin gain spectrum.
Since the brillouin frequency shift has a linear relationship with temperature and strain, a change in temperature or strain can be reflected by a change in brillouin frequency shift.
The BFS is obtained from BGS usually by using a lorentz fitting method, such as gaussian newton method and lm algorithm, but in this way, the accuracy of the result is often related to the set initial value, when the initial value is far away from the optimal value, the error obtained by calculation may fall into the local optimal value, so that the accuracy is poor, and a large number of iterations are required, and the processing time is long. In addition, a cross-correlation method can be adopted, and the existing research shows that the cross-correlation function of two Lorentz-shaped curves is also a Lorentz function, the cross-correlation of an ideal Lorentz-shaped curve and a Lorentz curve with noise is mainly determined by the correlation curve, but not by noise, after the cross-correlation operation is carried out, the frequency position with the largest cross-correlation degree is taken, and the value of the Brillouin frequency shift is reversely deduced; this approach results in correlation of the accuracy of the cross correlation method with the frequency scan interval, so in order to improve the accuracy of the cross correlation, interpolation of the brillouin gain spectrum is commonly used, but the interpolation operation requires a lot of processing time, which makes real-time temperature monitoring a troublesome problem.
Disclosure of Invention
The invention aims at overcoming the defects of the prior art and provides a rapid brillouin frequency shift extraction method. The invention has the characteristic of faster processing speed when the precision is equal to that of the traditional cross-correlation method.
The technical scheme adopted by the invention comprises the following steps:
(1) Constructing a cross-correlation feature matrix: the input data D consists of M Brillouin gain spectrums, the dimension is Mx N, wherein N is the data length of a single Brillouin gain spectrum, and the corresponding cross-correlation characteristic matrix P is constructed as follows:
wherein F is ,……,F L Is L curves with different peak positions and data length N; the construction method of the curve with different peak positions is as follows: the method comprises the steps of equally dividing frequencies in a selected sampling frequency range according to required precision, and constructing by adopting a Gaussian function or a Lorentz function;
(2) Constructing a brillouin frequency shift dictionary set: the brillouin frequency shift dictionary set S constructed according to the cross-correlation feature matrix is as follows:
wherein f 1 ,……,f L Is the curve F in the feature matrix P 1 ,……,F L The frequency of the corresponding peak position; (3) brillouin frequency shift search: by solving a cross-correlation degree matrix R of a cross-correlation characteristic matrix P and input data D, a curve position with the largest cross-correlation degree in R is searched out from S, and the frequency of a peak position corresponding to the curve position with the largest cross-correlation degree is taken as an output result bfs:
R=P·D T
bfs=S[argmax(R)]
wherein argmax is the curve position where the cross correlation is the greatest.
In the above technical solution, further, the lorentz curve has a functional form:
where g is the gain factor, v B For brillouin shift, v is the frequency and Δv is the full-width half maximum of the brillouin gain spectrum.
Further, the cross-correlation feature matrix P in step (1) has a dimension lxn, where L is the frequency range [ start, finish ] of the input data]And the frequency step size deltaf,
further, the frequency f of the brillouin shift dictionary set in step (2) 1 ,……,f L Corresponding to start, start+Δf, … …, start+lΔf, respectively.
In general, compared with the traditional cross-correlation brillouin frequency shift extraction method, the method has the following innovation points:
(1) Only the multiplication of the matrix is involved, the method is simple, and the processing time is fast.
(2) The constructed cross-correlation feature matrix can be simultaneously applied to the acquired M Brillouin gain spectrums, and interpolation of each spectrum is not needed.
(3) Because the traditional cross-correlation method uses Lorentz curves to carry out cross-correlation operation on input data, the extraction precision of Brillouin frequency shift is limited by frequency scanning intervals, the common precision improving mode is interpolation on the input data, the processing time is long, the capability of parallel processing data is poor, the invention creatively uses the Lorentz curves to construct a cross-correlation feature matrix, the matrix only needs to be initialized once, the subsequent extraction algorithm only involves simple matrix multiplication, the parallel processing data has good capability, and the processing time is fast.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the embodiments or the description of the prior art will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention.
Fig. 1 is a flow chart of a brillouin frequency shift extracting method provided by the invention;
FIG. 2 is a comparison of mean square errors of processing results of the method and the interpolation cross-correlation method of the present invention and the Levenberg-Marquardt algorithm under different line widths;
FIG. 3 is a comparison of processing times for an interpolated cross-correlation method and a method of the present invention, with the ordinate being the interpolated cross-correlation processing time/the processing time for the method of the present invention;
FIG. 4 is a mean square error comparison of processing results of an interpolation cross-correlation method and the method of the invention at different frequency scanning steps of 2MHz-18MHz and at intervals of 4MHz under different line widths;
FIG. 5 is a comparison of the interpolated cross-correlation and the processing time of the method of the present invention in the case of FIG. 5, with the interpolated cross-correlation processing time on the ordinate versus the processing time of the method of the present invention;
FIG. 6 is a graph of the data to be processed with a signal-to-noise ratio of 5dB versus a reference curve in a cross-correlation feature matrix;
fig. 7 is a graph of the degree of cross-correlation of data with a signal to noise ratio of 5dB to be processed after feature extraction.
Detailed Description
The present invention will be described in further detail with reference to the drawings and examples, in order to make the objects, technical solutions and advantages of the present invention more apparent. It should be understood that the detailed description is presented by way of example only and is not intended to limit the scope of the invention.
The accuracy of the cross correlation method in extracting the brillouin frequency shift is limited by the frequency scanning interval, and a great deal of time is required for interpolating the data to improve the accuracy. The method does not need interpolation, and can remarkably improve the data processing speed under the condition of being compared with the precision of the traditional cross-correlation method.
The method for deriving the brillouin shift according to the present embodiment is constructed by using the flowchart shown in fig. 1. Specifically, the whole method mainly comprises a data preprocessing part, a cross-correlation feature matrix is constructed, and a Brillouin frequency shift dictionary set and Brillouin frequency shift retrieval are constructed.
Data preprocessing
The data preprocessing module constructed in the embodiment is based on the fact that the BGS obtained through actual testing contain direct current bias, and under different environmental conditions, the BGS have different amplitudes, so that normalization processing is required to be carried out on data, and comparability exists among various data. For BGS data y= [ y ] 1 ,y 2 ,y 3 ,…,y n ]Corresponding sampling frequency range f= [ f 1 , 2 ,…,f n ]Normalized dataThe method comprises the following steps:
where mean (y) represents the average value of data y and std (y) represents the standard deviation of data.
Constructing a cross-correlation feature matrix and a Brillouin frequency shift dictionary set
The cross-correlation feature matrix P constructed in this embodiment is composed of a plurality of sets of lorentz curves with different brillouin frequency shifts, the sampling frequency range of the brillouin gain spectrum in this embodiment is [10630,10828] MHz, the frequency scanning interval is 2MHz, and according to the conventional cross-correlation method, the brillouin frequency shift extraction precision is limited to be more than 1MHz, and a frequency precision, for example, 0.2MHz, is set, so that the whole frequency band can be divided into 991 frequency points in total, and these frequency points form a brillouin frequency shift dictionary set, a brillouin gain spectrum can be generated according to each frequency point, wherein the line width is set to a given minimum line width, and the BGS line width range in the data set is [20,70 MHz ], and then the line width of the characteristic curve (i.e., the full-height half-width of the brillouin gain spectrum) is selected to be 20MHz, so that the 991×100 feature matrix can be constructed, and the 991 curves respectively correspond to the 991 frequency points, and the brillouin frequency shift dictionary needs to be constructed in a simplified process only once.
Brillouin frequency shift retrieval
The embodiment is mainly based on the fact that the greater the cross-correlation degree of the two curves is, the more similar the two curves are, and the closer the brillouin frequency shift of the two curves is. After the cross-correlation feature matrix is constructed, for the input data D, a cross-correlation degree matrix R is calculated for it:
R=P·D T
the corresponding frequency S [ pos ] is obtained as the brillouin shift of the brillouin gain spectrum from the position pos=argmax (R) where the maximum value appears in the cross correlation degree matrix.
It can be seen from fig. 2 that the cross-correlation method has better robustness to noise than lm (Levenberg-Marquardt) algorithm, and at the same time, under the condition of ensuring accuracy (the accuracy of the method is equivalent to that of the interpolation method), and the result of fig. 3 shows that the method has three orders of magnitude improvement in processing time compared with the interpolation method. Fig. 4 shows the results of the method and the interpolation cross-correlation method of the present invention under different frequency scanning intervals (2 MHz, 6MHz, 10MHz, 16 MHz), and it can be found that the two are of equal accuracy under the condition of full-width half-maximum of the fixed brillouin gain spectrum, while fig. 5 shows that the processing time of the method of the present invention is more improved relative to the processing time of the interpolation cross-correlation method as the frequency scanning interval increases. Fig. 6 shows one curve of a brillouin gain spectrum and a cross-correlation feature matrix with a signal-to-noise ratio of 5dB, fig. 7 shows a cross-correlation degree curve of the brillouin gain spectrum and the cross-correlation feature matrix with a signal-to-noise ratio of 5dB after operation, and the abscissa corresponds to frequencies in a brillouin frequency shift dictionary set, and the brillouin frequency shift corresponding to the input brillouin gain spectrum can be deduced to be about 10730MHz through the peak value of the cross-correlation degree curve.
The foregoing detailed description of the preferred embodiments and advantages of the invention will be appreciated that the foregoing description is merely illustrative of the presently preferred embodiments of the invention, and that no changes, additions, substitutions and equivalents of those embodiments are intended to be included within the scope of the invention.

Claims (4)

1. A brillouin shift extraction method, comprising the steps of:
(1) Constructing a cross-correlation feature matrix: the input data D consists of M Brillouin gain spectrums, the dimension is M x N, wherein N is the data length of a single Brillouin gain spectrum, and the corresponding cross-correlation characteristic matrix P is constructed as follows:
wherein F is 1 ,……,F L Is L curves with different peak positions and data length N; the construction method of the curve with different peak positions is as follows: the method comprises the steps of equally dividing frequencies in a selected sampling frequency range according to required precision, and constructing by adopting a Gaussian function or a Lorentz function;
(2) Constructing a brillouin frequency shift dictionary set: the brillouin frequency shift dictionary set S constructed according to the cross-correlation feature matrix is as follows:
wherein f 1 ,……,f L Is the curve F in the feature matrix P 1 ,……,F L The frequency of the corresponding peak position;
(3) Brillouin frequency shift search: by solving a cross-correlation degree matrix R of the cross-correlation feature matrix P and the input data D, a curve position with the greatest cross-correlation degree in R is retrieved from S, and a frequency at which a peak position corresponding to the curve position with the greatest cross-correlation degree is located is used as an output result bfs:
R=P·D T
bfs=S[argmax(R)]
wherein argmax is the curve position where the cross correlation is the greatest.
2. The brillouin shift extraction method according to claim 1, wherein the lorentz function has the form:
wherein g is the gain factor, v B For brillouin shift, v is the frequency and Δv is the full-width half maximum of the brillouin gain spectrum.
3. The brillouin shift extraction method according to claim 1, wherein the dimension of the cross correlation feature matrix P in the step (1) is lxn, where L is a frequency range [ start, finish ] of the input data]And the frequency step size deltaf,
4. a brillouin shift extraction method according to claim 3, wherein the brillouin shift dictionary set in step (2) has a frequency f 1 ,……,f L Corresponding to start, start+Δf, … …, start+lΔf, respectively.
CN202110608992.6A 2021-06-01 2021-06-01 Brillouin frequency shift extraction method Active CN113343173B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110608992.6A CN113343173B (en) 2021-06-01 2021-06-01 Brillouin frequency shift extraction method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110608992.6A CN113343173B (en) 2021-06-01 2021-06-01 Brillouin frequency shift extraction method

Publications (2)

Publication Number Publication Date
CN113343173A CN113343173A (en) 2021-09-03
CN113343173B true CN113343173B (en) 2023-07-21

Family

ID=77474102

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110608992.6A Active CN113343173B (en) 2021-06-01 2021-06-01 Brillouin frequency shift extraction method

Country Status (1)

Country Link
CN (1) CN113343173B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113819932B (en) * 2021-09-28 2023-05-02 北京卫星环境工程研究所 Brillouin frequency shift extraction method based on deep learning and mathematical fitting
CN117490738A (en) * 2023-10-31 2024-02-02 中国南方电网有限责任公司超高压输电公司广州局 Optical fiber sensing method, optical fiber sensing device, computer equipment, medium and computer product

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2825104A1 (en) * 2011-01-27 2012-08-02 Ramot At Tel Aviv University Ltd. Distributed and dynamical brillouin sensing in optical fibers
CN104296673A (en) * 2014-10-22 2015-01-21 中国电子科技集团公司第四十一研究所 Brillouin spectrum signal quality improving method
CN110260897A (en) * 2019-06-18 2019-09-20 华中科技大学 A kind of Brillouin optical time domain analysis instrument denoising method dictionary-based learning and system

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130018633A1 (en) * 2011-07-12 2013-01-17 University Of New Brunswick Method and apparatus for central frequency estimation

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2825104A1 (en) * 2011-01-27 2012-08-02 Ramot At Tel Aviv University Ltd. Distributed and dynamical brillouin sensing in optical fibers
CN104296673A (en) * 2014-10-22 2015-01-21 中国电子科技集团公司第四十一研究所 Brillouin spectrum signal quality improving method
CN110260897A (en) * 2019-06-18 2019-09-20 华中科技大学 A kind of Brillouin optical time domain analysis instrument denoising method dictionary-based learning and system

Also Published As

Publication number Publication date
CN113343173A (en) 2021-09-03

Similar Documents

Publication Publication Date Title
CN113343173B (en) Brillouin frequency shift extraction method
Pereda et al. Analyzing the stability of the FDTD technique by combining the von Neumann method with the Routh-Hurwitz criterion
US9046000B2 (en) Method for detecting foreign object damage in turbomachinery
Kuypers et al. Maximum accuracy evaluation scheme for wireless SAW delay-line sensors
CN116825243B (en) Multi-source data-based thermal barrier coating service life prediction method and system
CN114674352B (en) Distributed disturbance sensing and demodulation method based on Rayleigh scattering spectrum dissimilarity
Abbasnejad et al. FPGA-based implementation of a novel method for estimating the Brillouin frequency shift in BOTDA and BOTDR sensors
Wei et al. Inversion probability enhancement of all-fiber CDWL by noise modeling and robust fitting
CN106768895A (en) A kind of detection method of list, multimode fibre range self-adapting
Zhu et al. Improve accuracy and measurement range of sensing in km-level OFDR using spectral splicing method
CN112050942B (en) Optical fiber interference spectrum cavity length correction method based on phase compensation
CN107014409B (en) A kind of long range optical frequency domain reflection-based optical fiber Distributed Multi destabilization sensing method
KR20030045068A (en) Improved structure identification using scattering signatures
CN112082498A (en) Noise suppression sensing method based on phase measurement method OFDR strain and temperature
Chen et al. Fast algorithm for parameter estimation of LFM signals under low SNR
Abbasi et al. Estimation of time-of-flight based on threshold and peak analysis method for microwaves signals reflected from the crack
Osman et al. Method to improve fault location accuracy against cables dispersion effect
CN110274620B (en) Brillouin scattering signal denoising method based on spectrum center alignment
CN113155267A (en) OFDR system vibration detection method and system based on secondary correlation, storage medium and terminal
CN111707304A (en) Method for rapidly demodulating discrete cavity length of variable-step-length optical fiber F-P sensor
Wang et al. Study on wavelet denoising algorithm based on acoustic emission leakage of heaters
Mehdiyeva et al. Analysis of the mathematical model of AC signals and methods of digital measurements of their integral parameters
CN115265613A (en) Multi-frequency-interval Brillouin frequency shift extraction method and device
CN110823530A (en) Method for obtaining quality factor of micro-resonant cavity
Gao et al. DAVAR method for random noise signal process of FOG based on optimal window

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant