A φ-OTDR system data processing method to reduce the phase recovery error of the DCM algorithm

Abstract

The present invention belongs to the technical field of distributed optical fiber sensing of phase-sensitive optical time-domain reflectometry, and specifically relates to a φ-OTDR system data processing method for reducing the phase recovery error of a DCM algorithm. The method comprises the following steps: S1, reconstructing a two-dimensional space-time matrix of three-channel Rayleigh scattering signals output by the system; S2, performing time-domain shifted difference and accumulation processing on the two-dimensional space-time matrix respectively; S3, performing weighted averaging to obtain a multi-channel fused differential trace, judging whether there is vibration based on the differential trace, and if so, obtaining the vibration intrusion position; S4, interpolating the space-time matrix data at the extracted intrusion source position using a resampling method; and S5, performing phase demodulation processing on the resampled data using a DCM algorithm. The present invention reduces the phase recovery error of the DCM algorithm, expands the amplitude-frequency response range, and improves the anti-noise performance. Based on the present invention, high-speed linear and precise demodulation of long-distance and full vibration information can be achieved, thereby ensuring the stability and accuracy of the sensing system.

Description

Technical Field

[0001] The present invention belongs to phase sensitive optical time domain reflectometry In the field of distributed optical fiber sensing technology, specifically a method for reducing the phase recovery error of the DCM algorithm is provided. System data processing method. Background Art

[0002] In recent years, The system has sparked a research boom in the entire fiber optic sensing community due to its successful monitoring of ocean dynamics, volcanic events, and seismic waves. The systems can be divided into two categories: one is based on amplitude detection system, also known as distributed vibration sensing (DVS); the other is based on phase detection The system is also called distributed acoustic sensing (DAS). Compared with the former, the latter can not only locate the vibration signal, but also completely restore the waveform of the vibration signal. At present, many phase demodulation schemes have been proposed, such as digital coherent demodulation [Z. Pan, et al. ACP 2011, 8311 (8), 1–6], in-phase quadrature (I / Q) demodulation [Z. Wang, et al. Opt. Express, 24 (2), 853–858, 2016], 3×3 coupler demodulation [A. Masoudi, et al. Meas. Sci. Technol., 24 (8), 085204, 2013] and phase carrier generation (PGC) demodulation [Y. Shang, et al. Meas. J. Int. Meas. Confed., 79, 222–227, 2016]. Since the 3×3 coupler phase demodulation scheme adopts a direct detection structure and is relatively insensitive to polarization, it is favored by more and more researchers.

[0003] The differential cross multiplication (DCM) algorithm is an algorithm for realizing vibration information detection and phase recovery in a direct detection phase demodulation scheme, and has been widely used in many cases [Patent No.: CN202110571279.9; Heng Qian, et al. IEEE Photonics technol. Lett., 32(8), 473-476, 2020]. However, in In the system, signal processing is achieved by reconstructing the collected one-dimensional data into a two-dimensional space-time matrix and performing phase demodulation in the time domain. Since the time domain sampling rate of the discrete time system depends on the pulse repetition frequency, and the differential operation accuracy of the algorithm depends on the system sampling rate, when the system sampling rate is low, there is a large error between the demodulation result and the phase of the loaded vibration signal, the demodulated phase waveform is distorted, and the spectrum energy has harmonic dispersion and transfer. In addition, the demodulation result of the algorithm will be affected by factors such as light intensity fluctuations and signal distortion, which seriously reduces the accuracy. The phase demodulation performance of the system makes it difficult to meet the requirements of high-precision, low-cost, and large-scale applications in the field of optical fiber sensing. Summary of the Invention

[0004] The present invention overcomes the shortcomings of the prior art and aims to solve the following technical problems: providing a method for reducing the phase recovery error of the DCM algorithm. The system data processing method can increase the sampling rate and operation speed without increasing the cost, thereby improving the phase recovery quality, while also reducing the effects of noise and signal distortion.

[0005] In order to solve the above technical problems, the technical solution adopted by the present invention is: a method for reducing the phase recovery error of the DCM algorithm The system data processing method comprises the following steps:

[0006] S1. Obtain the Rayleigh scattering signals of the first, second, and third channels, and reconstruct the two-dimensional space-time matrices respectively;

[0007] S2. Perform time domain moving difference and accumulation processing on each channel data in the spatial domain to obtain three difference curves;

[0008] S3. Perform a weighted average of the three differential curves to obtain a multi-channel fused differential trace. Determine whether there is vibration based on the differential trace. If so, obtain the vibration intrusion location.

[0009] S4. Extract the signal strength of the three channels at the vibration intrusion location and perform interpolation operations on them respectively to obtain their respective interpolation functions to implement resampling. The calculation formula of the interpolation function is:

[0010] S i (t) = A i +B i t+C i t 2 +D i t 3 ;

[0011] Among them, t∈[t i ,t i+1 ], A i , Bi , C i , D i They are the interpolation function in the time interval [t i ,t i+1 ] interpolation coefficient, t i and t i+1 They represent the time corresponding to the i-th data point and the i+1-th data point respectively. The calculation formula of the interpolation coefficient is:

[0012]

[0013] Among them, I m (i) and I m (i+1) represents the intensity value of the i-th and i+1-th time domain sampling points corresponding to the vibration intrusion position of the m-th channel, h represents the time domain sampling time interval, α i and α i+1 Respectively represent the intermediate coefficients corresponding to the i-th and i+1-th time domain sampling points;

[0014] S5. Perform phase demodulation processing using DCM based on the resampled data obtained in S4 to obtain vibration phase information.

[0015] Preferably, the calculation formula of the intermediate coefficient is:

[0016]

[0017] Among them, I m (1)~I m (N time ) represent the 1st to Nth positions corresponding to the vibration intrusion position of the mth channel time The intensity value of the time domain sampling point, N time Indicates the number of sampling points in the time domain, α1~α Ntime Indicates the 1st to Nth time The intermediate coefficients corresponding to the sampling points.

[0018] Preferably, in S2, the formula for performing time domain moving difference processing on each channel data is:

[0019] Diff m (i,j)=Trace m (i+x,j)-Trace m (i,j) where 1≤i≤N time -x,1≤j≤N space ;

[0020] The formula for accumulation processing is:

[0021] where 1≤i≤Ntime -x,1≤j≤N space ;

[0022] Among them, Sum m (j) represents the cumulative intensity difference corresponding to the jth spatial domain sampling point of the mth channel, Diff m (i, j) represents the difference between the jth spatial domain sampling point and the ith time domain sampling point. m (i+x,j) and Trace m (i, j) represents the intensity value of the j-th spatial domain in the m-th channel at the i+x-th and i-th time domain sampling points, x represents the differential interval, N time Indicates the number of time domain sampling points.

[0023] Preferably, the value of x is greater than 5.

[0024] Preferably, in S3, the three differential curves are Gaussian filtered and normalized before weighted averaging.

[0025] Preferably, in S3, the method for normalizing the three difference curves is:

[0026] where 1≤j≤N space ;

[0027] Among them, Norm m (j) represents the normalized intensity value of the difference curve corresponding to the mth channel at the jth spatial domain sampling point, Sum m (j) represents the cumulative intensity difference corresponding to the jth spatial domain sampling point of the mth channel, min(Sum m ) and max(Sum m ) represent the minimum intensity value and the maximum intensity value of the mth difference curve respectively.

[0028] Preferably, in S3, the formula for performing weighted averaging on the three difference curves is:

[0029]

[0030] Among them, DS(j) represents the intensity value of the weighted average difference trace at the jth spatial domain sampling point, Norm m (j) represents the normalized intensity value of the difference curve corresponding to the mth channel at the jth spatial domain sampling point. w1, w2 and w3 are the weights of the three difference curves, which are calculated as follows:

[0031]

[0032] Among them, w mIndicates the weight of the differential curve corresponding to the mth channel, N space Indicates the number of spatial domain sampling points.

[0033] Preferably, in S5, according to the resampling result obtained in S4, phase demodulation processing is performed using DCM to obtain vibration phase information.

[0034] Compared with the prior art, the present invention has the following beneficial effects:

[0035] (1) The present invention uses the resampling method to interpolate the collected data, solving the problem The DCM phase demodulation algorithm in the system significantly reduces the phase recovery error and broadens the amplitude-frequency response range due to the discontinuous waveform of the acquired signal caused by the limited pulse repetition frequency.

[0036] (2) The present invention improves The DCM phase demodulation algorithm in the system reduces the demodulation phase low-frequency noise introduced by light intensity fluctuations, thereby improving the demodulation quality and the detection effect of the system;

[0037] (3) The present invention can complement the spatial domain Rayleigh scattering intensity of the three-way acquisition signal of the 3×3 coupler, suppress the interference fading effect, and significantly improve the positioning accuracy;

[0038] (4) The subsequent processing algorithm adopted by the present invention only needs to The system only needs to process the collected signals without increasing hardware costs;

[0039] (5) The signal processing method proposed in the present invention has the advantages of fast running speed, high demodulation accuracy, and applicability to different types of vibration scenarios.

[0040] In summary, the present invention proposes a method to reduce the phase recovery error of the DCM algorithm. The system data processing method has the beneficial effects of improving the demodulation accuracy of the DCM algorithm, reducing noise and signal distortion, reducing costs, and having a wide range of applications. The development of phase demodulation technology in the field of distributed optical fiber sensing provides new ideas and methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 A method for reducing the phase recovery error of the DCM algorithm provided by the embodiment of the present invention Flowchart of the system data processing method;

[0042] Figure 2 The embodiment of the present invention uses a 3×3 Michelson interferometer. Detection optical path device;

[0043] Figure 3 This is the phase demodulation flow chart of the DCM algorithm;

[0044] Figure 4 Comparison of the effects before and after multi-channel fusion positioning;

[0045] Figure 5 To compare the demodulated phase waveforms before and after the resampling method is implemented;

[0046] Figure 6 Comparison of the amplitude-frequency response range before and after the resampling method is implemented. DETAILED DESCRIPTION

[0047] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below. Obviously, the described embodiments are part of the embodiments of the present invention, not all the embodiments; based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0048] like Figure 1 As shown, the embodiment of the present invention provides a method for reducing the phase recovery error of the DCM algorithm. System data processing method, which The system generates a distributed Rayleigh scattering signal containing the intrusion source location and vibration information. It then reconstructs a two-dimensional space-time matrix from the first and second channel Rayleigh backscattering signals output by the 3×3 Michelson interferometer and the third channel Rayleigh scattering signal output by the third output port of the connected circulator. Each channel's data undergoes moving difference and superposition processing, and the resulting three-way differential traces are used to locate the vibration intrusion location using a channel fusion algorithm, achieving precise positioning for interference fading suppression. The space-time matrix data at the vibration location is extracted and resampled using an interpolation algorithm. This process effectively increases the number of sampling points. The resampled data is phase-demodulated using the DCM algorithm, expanding the DCM algorithm's amplitude-frequency response range and improving its noise immunity. This method includes two-dimensional space-time matrix reconstruction, multi-channel fusion positioning, vibration source identification and extraction, resampling, and DCM phase demodulation.

[0049] Specifically, the data processing flow of the present invention is as follows.

[0050] S1. Obtain the Rayleigh scattering signals of the first, second, and third channels, and reconstruct the two-dimensional space-time matrices respectively.

[0051] Specifically, if Figure 2 As shown, the present invention adopts The system is based on a 3×3 Michelson interferometer to achieve distributed acoustic / vibration sensing. The system consists of four parts: a narrow linewidth pulse light source modulation module 20, a Rayleigh scattering echo detection module 21, a 3×3 interferometric demodulation module 22, and a signal acquisition and processing module 23.

[0052] (1) The narrow-linewidth pulse light source modulation module 20 is used to achieve modulation and shaping of highly coherent pulses. This module includes a narrow-linewidth laser 201, a signal generator 202, an acousto-optic modulator 203, and a pulse light amplifier 204. The narrow-linewidth laser 201 emits highly coherent laser light, which is modulated by the acousto-optic modulator 203 into pulse light with a pulse width T and an optical frequency f. After being amplified by the pulse light amplifier 204, it is injected into the Rayleigh scattering echo detection module 21.

[0053] (2) The Rayleigh scattering echo detection module 21 is used to achieve backscattered Rayleigh light transmission and high-quality amplification. The module includes an optical circulator 211, a sensing fiber 212, an erbium-doped fiber amplifier 214, and a filter 215. The sensing fiber 212 adopts a single-mode fiber G652D. The backscattered Rayleigh light signal excited by the pulse in the fiber is further amplified by the erbium-doped fiber amplifier 214 to an average power of 10μW, and then enters the 3×3 interferometric demodulation module 22 through the 0.4nm bandwidth filter 214. In this embodiment, a piezoelectric ceramic tube 213 is set in the sensing fiber 212 to simulate the external intrusion signal, and various vibration modes such as triangle wave, square wave, and sine wave can be loaded.

[0054] (3) The 3×3 interference demodulation module 22 uses a 3×3 coupler 222 and a Michelson interferometer to achieve zero-difference detection and polarization-independent interference output, and the phase difference of the three output signals meets 120°. After the output signal of the Rayleigh scattering echo detection module 21 enters port 1 of the 3×3 coupler 222 through the circulator 221, it is output from ports 4 and 6 respectively and enters the two interference arms of the Michelson interferometer. The first Faraday rotator mirror 223 and the second Faraday rotator mirror 224 are set on the two interference arms of the Michelson interferometer. When reflecting the input light, its polarization principal axis is rotated 90°, and the polarization states of the two reflected light beams are orthogonal to the input light. After returning, the two reflected light beams interfere inside the 3×3 coupler 222 and are output from ports 1, 2, and 3, with the modulation phases being 120°. Among them, ports 1, 2, and 3 are located on one side of the 3×3 coupler 222, and ports 4, 5, and 6 are located on the other side of the 3×3 coupler 222. The signal output from port 1 is outputted through the circulator 221 , and then, together with the signals output from ports 2 and 3 , are injected into the signal acquisition and processing module 23 .

[0055] (4) The signal acquisition and processing module 23 includes a first avalanche photodetector 331, a second avalanche photodetector 232, a third avalanche photodetector 233, a four-channel high-speed data acquisition card 234, and a computer control mainboard 235. The input ends of the first avalanche photodetector 231, the second avalanche photodetector 232, and the third avalanche photodetector 233 are respectively connected to port 3 of the circulator 221 and port 2 and port 3 of the 3×3 coupler 222, and the output ends are connected to the four-channel high-speed data acquisition card 234. The communication buckle of the four-channel high-speed data acquisition card 234 is connected to the computer control mainboard 235 for data upload and download, and its trigger output port is connected to the signal generator 202 for controlling the start and end of detection of the system.

[0056] In this embodiment, The three-way data points output by the high-speed data acquisition card 234 in the system are used to construct a time-space matrix. The data points are reconstructed in two dimensions and output three Rayleigh scattering matrices (referred to as Trace1, Trace2, and Trace3) with a phase difference of 120°. The size of the Rayleigh scattering matrix is ​​N. time ×N space , N space Indicates the number of spatial domain sampling points, which is determined by the oscilloscope sampling rate f osc and the pulse repetition frequency f pulse The ratio of N space =f osc / f pulse ; N time Represents the number of time domain sampling points, corresponding to the number of light pulses emitted.

[0057] S2. Perform time domain moving difference processing on each channel data in the spatial domain to obtain three difference curves.

[0058] In this embodiment, positioning is performed through multi-channel fusion. The specific method is: performing time-domain moving difference and accumulation operations on the three-channel Rayleigh scattering matrices Trace1, Trace2, and Trace3 output by the time-space matrix operation, and then performing spatial domain differential trace fusion operations to obtain reliable positioning of the intrusion source. The specific calculation principle is as follows:

[0059] Since the Rayleigh scattering intensity corresponding to any point in the spatial domain is the superposition of the Rayleigh scattered light waves of the scattering points within the pulse width, the change in Rayleigh scattering intensity is random and independent of each other. Single-mode optical fiber can be regarded as a sensor array composed of a large number of point sensors, which can realize distributed multi-point monitoring. When an intrusion event occurs near a certain position on the optical fiber, the scattering points in the optical fiber are affected by the strain and an elastic-photoelectric effect occurs. The phase of the Rayleigh scattered light at this location is modulated and the intensity changes, which is reflected as time domain intensity fluctuations. When performing moving differences on the Rayleigh scattering matrices Trace1, Trace2 and Trace3 of the three channels, the operation will be performed in the time domain, with a difference interval of x, and three dimensions of (N time -x)×N space The differential time-space matrix Diff1, Diff2 and Diff3. Taking the mth channel data as an example (m=1,2,3), the differential intensity value Diff obtained by the jth spatial domain sampling point and the ith time domain sampling point m (i,j) is:

[0060] Diff m (i, j)=Trace m (i+x, j)-Trace m (i, j) where 1≤i≤N time -x, 1≤j≤N space ;(1)

[0061] Where, Trace m (i+x,j) and Trace m (i, j) represents the intensity value of the jth spatial domain sampling point at the i+xth and ith time domain sampling points in the mth channel. x represents the differential interval. The choice of x can be determined according to the difference of the actual signal. It is generally greater than 5 to ensure a good differential intensity. At this time, the size of the differential matrix is ​​(N time -x)×N space The three differential time-space matrices Diff1, Diff2 and Diff3 are respectively accumulated and averaged in the time domain to obtain three intensity difference curves (denoted as Sum1, Sum2 and Sum3). Then the cumulative intensity difference Sum corresponding to the jth spatial domain sampling point of the mth channel is m (j) is written as:

[0062] where 1≤i≤N time -x,1≤j≤N space ;(2)

[0063] Among them, Sum m Represents the differential curve of the mth channel.

[0064] S3. Perform weighted averaging on the three differential curves to obtain a multi-channel fused differential trace. Determine whether there is vibration based on the differential trace. If so, obtain the vibration intrusion location.

[0065] Due to the phase transmission characteristics of the 3×3 coupler, the intensities of the three differential curves at different locations are complementary, so they can be fused to improve positioning accuracy. Gaussian filtering and normalization are performed on the three differential curves Sum1, Sum2, and Sum3, respectively, to map their intensity values ​​to the range of [0, 1]. The specific normalization method is:

[0066] where 1≤j≤N space ;(3)

[0067] Among them, Norm m (i) represents the normalized intensity value of the difference curve of the mth channel at the jth spatial domain sampling point, min(Sum m ) and max(Sum m ) represent the minimum intensity value and the maximum intensity value of the difference curve of the mth channel respectively.

[0068] By performing weighted averaging on the differential curves Norm1, Norm2, and Norm3 of the three normalized channels, the multi-channel fused differential trace DS is obtained to improve the accuracy of the fused differential trace. The intensity value DS(j) of the weighted averaged differential trace at the jth spatial domain sampling point is:

[0069]

[0070] Where w1, w2 and w3 are the weights of the differential curves of the three channels, and their sum is 1, which represents the importance ratio of the three differential traces in calculating the weighted average. The weight w of the mth channel is m The calculation method is:

[0071]

[0072] The presence of an intrusion source can be determined and its location information extracted through a multi-threshold segmented search of the differential trace. This operation first performs a multi-threshold segmented search on the multi-channel fused differential trace DS. When the intensity of the j-th spatial domain sampling point exceeds the threshold, vibration is determined to be occurring at that point. Otherwise, no vibration is detected, and subsequent operations are terminated.

[0073] like Figure 4 As shown, (a) is the locally enlarged time domain waveform of the Rayleigh scattering traces (Trace1, Trace2 and Trace3) of the three channels; (b) is the comparison between the differential traces after multi-channel fusion and the differential traces before fusion. Figure 4 As shown in (a), the Rayleigh scattering traces of the three channels are complementary in intensity, which can suppress the influence of interference fading noise. Figure 4 The differential trace results in (b) verify the noise suppression effect. The differential trace after multi-channel fusion can eliminate false alarm points and improve positioning accuracy.

[0074] S4. Extract the signal strengths of the three channels at the vibration intrusion location, perform interpolation operations on them respectively, and obtain their respective interpolation functions to achieve resampling.

[0075] In this embodiment, the time domain intensity variation curve Trace at the extracted intrusion position (corresponding to the Jth spatial domain sampling point) is obtained. m (N time ,J), use the interpolation algorithm to interpolate data to improve the phase demodulation accuracy. The specific implementation process is as follows:

[0076] Step 1: Extract Trace based on the intrusion source location information provided after vibration source judgment and extraction (assuming that there is an intrusion event at the Jth spatial domain sampling point). m The intensity sequence signals corresponding to the J-th spatial domain sampling point are recorded as I1, I2 and I3 for subsequent processing.

[0077] Step 2: Perform Gaussian filtering on the I1, I2, and I3 signals to remove the effects of noise and signal distortion.

[0078] Step 3: Perform interpolation operations on the filtered values ​​I1, I2, and I3 respectively. The number of interpolation data points is selected according to user requirements. The precision, accuracy, smoothness, and computational complexity of the interpolated SI1, SI2, and SI3 meet the requirements of phase demodulation. The calculation principle of the interpolation function is as follows:

[0079] Taking the mth channel data as an example (m=1,2,3), for signal I m , in the time domain axis, the number of sampling points is N time , the corresponding time interval is [t1,t Ntime ]. Perform cubic polynomial interpolation between every two time nodes, and the interpolation function corresponding to the i-th data point to be interpolated is:

[0080] SI i (t) = A i +B i t+C i t 2 +D i t 3 ; (6)

[0081] Among them, t∈[t i ,t i+1], A i , B i , C i , D i They are the interpolation functions in the interval [t i ,t i+1 ], the value range of i is {1≤i≤N time -2}. To ensure that the interpolated waveform meets the requirements of smoothness and continuity, three calculation conditions must be met: ①SI i (t i )=I m (i), SI i (t i+1 )=I m (i+1);②SI i '(t i+1 )=SI i+1 '(t i+1 );③SI i ”(t i+1 )=SI i+1 ”(t i+1 ). Due to t i+1 and t i The intervals between them are equal and are recorded as h. The following expression holds:

[0082] A i =I m (i) (7)

[0083] A i +B i h+C i h 2 +D i h 3 =I m (i+1) (8)

[0084] B i +2C i h+3D i h 2 =B i+1 (9)

[0085] 2C i +6D i h=2C i+1 (10)

[0086] Among them, I m (i) and I m (i+1) represents the intensity values ​​of the i-th and i+1-th time domain sampling points of the m-th channel.

[0087] Let the intermediate coefficient α i =2C i, substituting into formula (11), we can get:

[0088]

[0089] A i , D i Substituting into formula (9), we can get:

[0090]

[0091] A i , B i , C i , D i Substituting into formula (10), we can see that:

[0092]

[0093] Under natural boundary conditions, α1=0, Therefore, according to formula (14), the following linear equations can be constructed, and the intermediate coefficient vector α has a unique solution:

[0094]

[0095] Among them, I m (1)~I m (N time ) represent the 1st to Nth channels of the mth channel respectively. time The intensity value of the time sampling point, N time Indicates the number of time domain sampling points, α1~α Ntime Indicates the 1st to Nth time The intermediate coefficients corresponding to the sampling points.

[0096] By solving the above linear equations, we can get the intermediate coefficient vectors α1~α related to the m-th channel acquisition signal. Ntime Furthermore, the interpolation function SI of the mth channel collected signal i The coefficient (A) of (t) i , B i , C i , D i ) can be calculated by the following formula:

[0097]

[0098] Among them, I m (i) and I m (i+1) represents the intensity values ​​corresponding to the i-th and i+1-th time domain sampling points extracted from the spatial domain sampling point J at the intrusion position of the m-th channel, that is, I m (i) = Trace m (i,J). h represents ti+1 and t i The time interval between them satisfies the relationship h=1 / f pulse , α i and α i+1 Respectively represent the intermediate coefficients corresponding to the i-th and i+1-th time domain sampling points;

[0099] S5. Based on the resampled data obtained in step S4, vibration phase information is obtained by DCM phase demodulation.

[0100] In this embodiment, the DCM phase demodulation method used is the existing DCM phase demodulation method. The phase demodulation process is shown in Figure 3 The specific implementation process is as follows:

[0101] Step 1: Remove DC influence through DC removal module 30. The specific calculation method is to obtain DC signals by adding the resampled signals SI1, SI2 and SI3 output by the resampling module through adder 301. The DC signals are added to SI1, SI2 and SI3 through adders 302, 303 and 304 respectively to remove DC influence, and obtain the DC-removed signals a, b and c;

[0102] Step 2: Calculate the vibration phase differential information through the differential cross-multiplication module 31. The specific calculation method is as follows: a, b, and c enter differentiators 311, 312, and 313 for forward difference quotient processing to obtain d, e, and f; a, b, c, d, e, and f enter the cross-multiplier 314 to obtain the signal X containing the vibration phase differential information;

[0103] Step 3: Calculate the phase differential value by removing the light intensity fluctuation module 32. The specific calculation method is to send the DC-removed signals a, b and c to the squarers 321, 322 and 323 and the adder 324 to obtain the light intensity fluctuation term M, which is then sent to the multiplier 325 and multiplied by the coefficient As the denominator of the divider 326, an X / M operation 326 is performed to obtain a phase differential value;

[0104] Step 4: Run the integration module 33. Perform an integration operation 331 on the value obtained by X / M to obtain the vibration phase change φ to be measured.

[0105] In this embodiment, the resampling method is used to reduce The phase recovery error of the DCM algorithm in the system is due to the fact that when the DCM algorithm demodulates the phase, the calculation results of the differentiators 311, 312 and 313 based on Taylor expansion are highly dependent on the continuity of the two signals in the time domain. For example, at the measurement time t i The differential value d of the DC-free signal a of the first channel is approximately calculated as d(t i )=[a(t i+1)-a(t i )] / h, the corresponding error δ d (t i )for:

[0106]

[0107] Among them, a (n) is the nth-order derivative function of a, and the mathematical form of a is:

[0108]

[0109] Where D and f are the intrusion signal amplitude and frequency respectively, t is the time domain sampling time, φ envir is the phase introduced by the environmental disturbance, h represents t i+1 and t i The time interval between them satisfies the relationship: h = 1 / f pulse δ e (t i ) and δ f (t i ) is calculated using the same method. d , δ e and δ f The size of will be transmitted in the following calculation process, and the final demodulation phase error is:

[0110]

[0111] From equations (16) and (18), we can see that the demodulation phase error is determined by the intrusion signal amplitude, frequency and f pulse Ok. Due to the traditional System f pulse It is not adjustable and is limited to the kHz level due to the km-level optical fiber length. Therefore, the demodulation phase error is often large and cannot accurately restore the vibration phase change information. Therefore, the embodiment of the present invention uses the resampling method to interpolate the number of data points to achieve adjustable time domain sampling rate, and the demodulation phase error is no longer limited to f pulse , the error will be greatly reduced. Figure 5 (a) is the comparison of the demodulated phase time domain waveform before and after resampling, and (b) is the comparison of the spectrum signal-to-noise ratio. It can be seen from the figure that the recovery quality of the demodulated phase obtained by the data processing method provided by the embodiment of the present invention is effectively improved, the light intensity fluctuation noise is suppressed, and the spectrum signal-to-noise ratio is improved by 52.5dB. Figure 6 The root mean square error (RMSE) distribution of the demodulated phase before (a) and after (b) resampling is compared. As can be seen from the figure, the amplitude-frequency response range (Response area) of the demodulated phase is greatly widened by the resampling method in the embodiment of the present invention.

[0112] In summary, the present invention discloses a method for reducing the phase recovery error of the DCM algorithm. This system data processing method preprocesses the three-channel sampling signals of the #imgpt46# system based on a 3×3 coupler using a multi-channel fusion positioning algorithm and resampling methods. This significantly reduces the phase recovery error of the DCM algorithm and eliminates the influence of light intensity waveform noise. This method can achieve high-speed, linear, and precise phase demodulation of full vibration information over long distances, ensuring the stability and accuracy of the sensing system. This method has great application prospects and market demand.

[0113] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method to reduce the phase recovery error of DCM algorithm φ -OTDR system data processing method, characterized in that, The following steps are involved: S1. Obtain the Rayleigh scattering signals of the first, second, and third channels, and reconstruct the two-dimensional space-time matrices respectively; S2. Perform time-domain moving difference and accumulation processing on each channel data in the spatial domain to obtain three difference curves; S3. Perform a weighted average of the three differential curves to obtain a multi-channel fused differential trace. Determine whether there is vibration based on the differential trace. If so, obtain the vibration intrusion location. S4. Extract the signal strength of the three channels at the vibration intrusion location and perform interpolation operations on them respectively to obtain their respective interpolation functions to implement resampling. The calculation formula of the interpolation function is: ； in, t∈ , A i , B i , C i , D i They are the interpolation function in the time interval The interpolation coefficients on and Respectively represent i data points and i +1 data point corresponding to the time, the interpolation coefficient is calculated as follows: ； in, and They represent the mth channel corresponding to the vibration intrusion position. i and i +1 intensity value of the time domain sampling point, h represents the sampling time interval in the time domain, and Respectively represent i and i +The intermediate coefficient corresponding to 1 time domain sampling point; S5. Perform phase demodulation processing using DCM based on the resampled data obtained in S4 to obtain vibration phase information; The calculation formula of the intermediate coefficient is: ； in, Respectively represent m The channel is located at the first to The intensity value of the time domain sampling point, Indicates the number of sampling points in the time domain, Indicates the 1st to The intermediate coefficient corresponding to the sampling points; In S2, the formula for performing time domain moving difference processing on each channel data is: ,in ; The formula for accumulation processing is: ,in ; in, represents the number of spatial domain sampling points, Indicates the m The channel in j The cumulative intensity difference corresponding to the spatial domain sampling points, Indicates the j The spatial domain sampling point is i The difference value of the time domain sampling points, and Indicates the m In the channel j The spatial domain in i + x and i The intensity value of the time domain sampling point, x represents the difference interval, Indicates the number of time domain sampling points.

2. A method for reducing the phase recovery error of the DCM algorithm according to claim 1 φ -OTDR system data processing method, characterized in that, described x The value of is greater than 5.

3. A method for reducing the phase recovery error of the DCM algorithm according to claim 1 φ -OTDR system data processing method, characterized in that, In S3, the three difference curves are Gaussian filtered and normalized before weighted averaging.

4. A method for reducing the phase recovery error of the DCM algorithm according to claim 1 φ -OTDR system data processing method, characterized in that, In S3, the method for normalizing the three difference curves is: , in ; in, Indicates the m The differential curve corresponding to each channel is j The normalized intensity value of the spatial domain sampling point, Indicates the m The channel in j The cumulative intensity difference corresponding to the spatial domain sampling points, and Respectively represent m The minimum and maximum intensity values ​​of the difference curve.

5. A method for reducing the phase recovery error of the DCM algorithm according to claim 1 φ -OTDR system data processing method, characterized in that, In S3, the formula for weighted averaging the three difference curves is: ； in, DS ( j ) indicates the weighted average difference trace at the j The intensity value of the spatial domain sampling point, Indicates the m The differential curve corresponding to each channel is j The normalized intensity value of the spatial domain sampling point, are the weights of the three difference curves, which are calculated as follows: ； in, Indicates the m The weight of the difference curve corresponding to each channel.

6. A method for reducing the phase recovery error of the DCM algorithm according to claim 1 φ -OTDR system data processing method, characterized in that, In S5, based on the resampling result obtained in S4, phase demodulation processing is performed using DCM to obtain vibration phase information.

Images