A method and system for determining vibration propagation direction by phase sensitive optical time domain reflectometry

Abstract

The present application relates to the technical field of optical fiber vibration monitoring, in particular to a phase-sensitive optical time domain reflectometry vibration propagation direction determination method and system, comprising obtaining full-line echo in-phase component and quadrature component data; constructing complex echo sequence of each distance gate, and performing self-calibration compensation on the preprocessed phase sequence to obtain the compensated phase sequence; based on the quality index, screening the candidate distance gate, estimating the arrival time delay of each gate pair in the forward gate pair set and constructing gate pair weight, robustly fusing the arrival time delay of each gate pair based on the gate pair weight, constructing the reverse gate pair set corresponding to the forward gate pair set and obtaining the reverse fusion time delay, performing symmetric difference on the forward fusion time delay and the reverse fusion time delay, comparing the unbiased time delay with the preset direction determination threshold, and outputting the consistency index reflecting the consistency degree of multi-gate pair estimation, the present application improves the accuracy and long-term operation consistency of the direction determination result.

Description

Technical Field

[0001] This invention relates to the field of fiber optic vibration monitoring technology, and in particular to a method and system for determining the vibration propagation direction using phase-sensitive optical time-domain reflection. Background Technology

[0002] Distributed fiber optic vibration monitoring technology has been widely used in perimeter security, pipeline inspection, rail transit and other fields. The distributed vibration monitoring system based on phase-sensitive optical time-domain reflectometer (Φ-OTDR) can realize vibration detection distributed along the optical fiber. When external vibration or strain acts on a certain section of optical fiber, it causes a phase change in the corresponding backscattered signal, that is, to obtain the vibration response information at each position along the line. The spatial sampling unit discretely divided along the optical fiber is usually called the distance gate.

[0003] In existing technologies, vibration propagation direction determination based on single-end access often employs the adjacent distance gate cross-correlation method. This method selects two adjacent distance gates near the target event location, extracts the vibration response sequences corresponding to the two distance gates, and uses the delay corresponding to the maximum value of the cross-correlation function as the arrival time difference between the two distance gates. The vibration propagation direction is then determined based on the sign of this time difference.

[0004] However, the aforementioned direction determination method based on the cross-correlation of adjacent distance gates has significant shortcomings when used in engineering fields. First, this method is susceptible to phase drift caused by laser frequency shift and ambient temperature changes. When the phase baselines of different distance gates drift inconsistently, the main peak position of the cross-correlation curve is prone to shift, leading to fluctuations in the estimated time difference of arrival. Second, under Rayleigh scattering fading and polarization state changes, the echo quality of individual distance gates decreases, weakening the ability of the cross-correlation results to characterize the actual propagation delay. For vibrations with small amplitude or low frequency, the misjudgment rate increases significantly. In addition, there are slight differences in the equivalent bandwidth and group delay response characteristics of different distance gates, resulting in a fixed bias in the estimation results. This bias is difficult to eliminate in the long term using a simple threshold, leading to regular deviations in the direction conclusions and poor long-term repeatability. At present, a vibration propagation direction determination method and system based on phase-sensitive optical time-domain reflection is needed. Summary of the Invention

[0005] To address the issues of insufficient stability and high misjudgment rate in direction determination under the constraints of single-end access and a single sensing fiber, and under conditions of laser frequency shift, slow phase drift, and Rayleigh scattering fading, this invention provides a method and system for determining the vibration propagation direction using phase-sensitive optical time-domain reflection.

[0006] In a first aspect, the present invention provides a method for determining the vibration propagation direction of phase-sensitive optical time-domain reflection, which adopts the following technical solution: A method for determining the vibration propagation direction of phase-sensitive optical time-domain reflectometry, comprising: Acquire the in-phase and quadrature components of the full-line echo collected by a single-ended single-fiber phase-sensitive optical time-domain reflectometry system; Construct complex echo sequences for each range gate, extract the phase sequences for each range gate in the slow time dimension, and perform preprocessing operations on the phase sequences to obtain preprocessed phase sequences; A set of reference range gates is selected outside the target event region. The common drift term is estimated based on the set of reference range gates, and the preprocessed phase sequence is self-calibrated and compensated to obtain the compensated phase sequence. Within the target event area, candidate distance gates are screened based on quality indicators to obtain an effective distance gate set, and a set of forward gate pairs that satisfy the preset gate spacing constraint is constructed within the effective distance gate set; Arrival delay is estimated for each gate pair in the forward gate pair set and gate pair weights are constructed. Based on the gate pair weights, the arrival delays of each gate pair are robustly fused to obtain the forward fused delay. Construct a set of reverse gate pairs corresponding to the set of forward gate pairs and obtain the reverse fusion delay. Perform symmetric difference between the forward fusion delay and the reverse fusion delay to obtain the unbiased delay. The unbiased time delay is compared with a preset direction determination threshold. Based on the comparison result, the propagation direction of vibration along the fiber axis is output, and a consistency index reflecting the degree of consistency of the multi-gate pair estimation is also output.

[0007] Furthermore, the extraction of the phase sequence of each distance gate in the slow time dimension includes defining an index identifier, defining the spatial sampling unit after discretization along the optical fiber as a distance gate according to the index identifier, and defining the pulse sequence dimension as slow time. The in-phase and quadrature components corresponding to each distance gate and each slow time are synthesized into a complex echo sequence. Based on the complex echo sequence, the phase principal value sequence is extracted through a four-quadrant arctangent operation. The expression for the phase principal value sequence is: , in, The phase principal value of the i-th distance gate at the m-th slow time sampling point. For orthogonal components, For in-phase components, This is the arctangent operation for the four quadrants.

[0008] Furthermore, the preprocessing operation of the phase sequence includes phase expansion of the phase principal value sequence, detrending of the expanded continuous phase sequence to obtain a detrended phase sequence, and finally bandpass filtering of the detrended phase sequence to obtain a preprocessed phase sequence. The phase unrolling includes: calculating the principal phase difference of the i-th distance gate between adjacent slow-time sampling points, and then flipping the principal phase difference back to the interval. The corrected principal value difference is obtained internally, and then accumulated point by point in slow time to obtain the expanded continuous phase sequence: , in, Let be the continuous phase sequence after expansion of the i-th distance gate at the m-th slow time sampling point. To reflect the principal phase difference back to the interval The corrected principal value difference obtained within the inner, It is the continuous phase sequence after the expansion of the i-th distance gate at the (m-1)-th slow time sampling point.

[0009] Furthermore, the self-calibration compensation of the preprocessed phase sequence includes setting a protection distance outside the target event region and determining the distance gates located outside the protection distance as the reference distance gate set; For each slow time, the median of the preprocessed phase sequence of each range gate in the reference range gate set is taken to estimate the common drift term; The common drift term is removed from the preprocessed phase sequence of each range gate to obtain the compensated phase sequence: , in, The phase sequence after compensation for the distance gate z at slow time t. This is the preprocessed phase sequence of the distance gate z at slow time t. This is the estimate of the common drift term at slow time t.

[0010] Furthermore, the process of filtering candidate distance gates based on quality indicators includes determining a set of candidate distance gates within the target event region; Calculate the quality index of each candidate distance gate, wherein the quality index includes at least one of amplitude quality index, amplitude stability index, phase availability index and related availability index; The quality indicators are compared with the corresponding preset thresholds. Candidate distance gates that meet the screening conditions are included in the effective distance gate set, while candidate distance gates that do not meet the screening conditions are removed from the effective distance gate set. The relevant availability metrics are calculated using the following formula: , in, For the i-th candidate distance gate, the relevant availability metric is... The amplitude of the main peak of cross-correlation, This is the amplitude of the second peak. To prevent extremely small constants with a denominator of zero.

[0011] Furthermore, the construction of the forward door pair set satisfying the preset door spacing constraint includes selecting a forward door pair consisting of two effective distance doors from the effective distance door set, wherein each forward door pair is denoted as (i,j), and The door spacing between the two doors meets the preset constraints. ; All forward door pairs that satisfy the door spacing constraint are defined as the set of forward door pairs: , in, For the forward gate pair set, For the effective distance gate set, and For effective distance gate index, To preset the minimum door spacing, This is the preset maximum door spacing.

[0012] Furthermore, the symmetrical difference between the forward fusion delay and the reverse fusion delay includes calculating the discrete cross-correlation function for each gate pair in the forward gate pair set based on the compensated preprocessed phase sequence, and taking the delay corresponding to the maximum value of the cross-correlation function as the estimated arrival delay value of the gate pair. Gate pair weights are constructed based on the ratio of the main peak amplitude to the secondary peak amplitude of the cross-correlation of each gate pair, and the arrival delay estimates of each gate pair are weighted and fused by median to obtain the forward fusion delay. Construct a set of reverse door pairs corresponding to the set of forward door pairs, wherein each door pair in the set of reverse door pairs satisfies the same door spacing constraint as the set of forward door pairs; The forward fusion delay is obtained by differentiating the backward fusion delay and taking half of the result. The expression for the forward fusion delay is as follows: , in, For forward fusion latency, Let be the estimated arrival delay of gate pair (i,j). Let (i,j) be the weight of the gate pair. For the forward gate pair set, This is the weighted median operator.

[0013] Furthermore, the output reflects a consistency index of the degree of consistency of the multi-gate pair estimation, including comparing the absolute value of the unbiased delay with a preset direction determination threshold. When the absolute value of the unbiased delay is less than or equal to the preset direction determination threshold, an uncertain state is output. When the unbiased delay is greater than the preset direction determination threshold, it is determined that the vibration propagates from the near end to the far end. When the unbiased delay is less than the negative preset direction determination threshold, it is determined that the vibration propagates from the far end to the near end. Based on the deviation between the estimated arrival delay of each gate pair and the fusion delay, and combined with the gate pair weights, a consistency index is calculated: , in, As a consistency indicator, For a set of gate pairs, Let (a, b) be the weight of the gate pair. Let be the estimated arrival time delay of gate pair (a,b). This is the estimated delay value after fusion. This is the consistency tolerance threshold. This is an indicator function.

[0014] Secondly, a vibration propagation direction determination system based on phase-sensitive optical time-domain reflectometry includes: Narrow linewidth laser, beam splitting unit, local oscillator branch conditioning module, coherent receiving module, high-speed data acquisition module and digital signal processing module; The narrow linewidth laser is connected to the beam splitting unit, which divides the laser into an emission branch and a local oscillator branch. The emission branch is connected to the sensing fiber via a pulse modulation module, an optical amplification module, and the first port of an optical circulator in sequence. The local oscillator branch is connected to the local oscillator input terminal of the coherent receiving module via the local oscillator branch conditioning module. The third port of the optical circulator is connected to the signal input terminal of the coherent receiving module. The output terminal of the coherent receiving module is connected to the input terminal of the high-speed data acquisition module. The output terminal of the high-speed data acquisition module is connected to the input terminal of the digital signal processing module.

[0015] Furthermore, the digital signal processing module is used to process the digital signal output by the high-speed data acquisition module to determine the propagation direction of the vibration along the axis of the sensing fiber.

[0016] In summary, the present invention has the following beneficial technical effects: 1. This invention estimates the common drift term by selecting a set of reference distance gates outside the target event region and performs self-calibration compensation on the preprocessed phase sequence, thereby effectively suppressing the slow phase drift caused by laser frequency shift and ambient temperature change, and improving the long-term stability and repeatability of direction determination.

[0017] 2. This invention obtains an effective range gate set by screening candidate range gates based on quality indicators within the target event region, and performs weighted robust fusion of the arrival delays of multiple forward gate pairs, thereby achieving adaptive elimination of abnormal range gates under Rayleigh scattering fading and polarization state change conditions, and reducing the direction misjudgment rate.

[0018] 3. By constructing a set of reverse gate pairs corresponding to the forward gate pair set and performing symmetrical difference, this invention achieves self-elimination of the systematic bias introduced by the difference in equivalent response of distance gates and fixed link delay, thereby improving the accuracy of direction determination results and long-term operational consistency.

[0019] 4. This invention achieves synchronous output of vibration propagation direction and confidence level by comparing unbiased time delay with a preset direction determination threshold and outputting a consistency index that reflects the consistency of multi-gate pair estimation, thus providing a quantitative basis for secondary decision-making in engineering alarms. Attached Figure Description

[0020] Figure 1 This is a schematic diagram of the overall process of a method for determining the vibration propagation direction of phase-sensitive optical time-domain reflection according to an embodiment of the present invention.

[0021] Figure 2 This is a schematic diagram of the phase processing flow according to an embodiment of the present invention.

[0022] Figure 3 This is a schematic diagram of the spatial relationship between the reference gate and the target gate in an embodiment of the present invention.

[0023] Figure 4 This is a schematic diagram of multi-gate pairs and symmetric differential debiasing according to an embodiment of the present invention.

[0024] Figure 5 This is a flowchart of the propagation direction output based on unbiased time delay and direction determination threshold in an embodiment of the present invention. Detailed Implementation

[0025] The present invention will be further described in detail below with reference to the accompanying drawings.

[0026] Example 1 Reference Figure 1 This embodiment of a method for determining the vibration propagation direction of phase-sensitive optical temporal reflection includes: S1. Acquire the in-phase and quadrature components of the full-line echo collected by the single-ended single-fiber phase-sensitive optical time-domain reflectometry system; S2. Construct complex echo sequences for each range gate, extract the phase sequence of each range gate in the slow time dimension, and perform preprocessing operations on the phase sequence to obtain the preprocessed phase sequence. S3. Select a set of reference range gates outside the target event region, estimate the common drift term based on the set of reference range gates, and perform self-calibration compensation on the preprocessed phase sequence to obtain the compensated phase sequence; S4. Within the target event area, candidate distance gates are screened based on quality indicators to obtain a set of effective distance gates, and a set of forward gate pairs that satisfy the preset gate spacing constraint is constructed within the set of effective distance gates. S5. Estimate the arrival delay of each gate pair in the forward gate pair set and construct the gate pair weights. Based on the gate pair weights, robustly fuse the arrival delays of each gate pair to obtain the forward fusion delay. Construct a set of reverse gate pairs corresponding to the set of forward gate pairs and obtain the reverse fusion delay. Perform symmetric difference between the forward fusion delay and the reverse fusion delay to obtain the unbiased delay. S6. Compare the unbiased delay with the preset direction determination threshold, output the propagation direction of vibration along the fiber axis based on the comparison result, and output a consistency index reflecting the degree of consistency of the multi-gate pair estimation.

[0027] Specifically, a method for determining the vibration propagation direction of phase-sensitive optical time-domain reflectometry includes the following: like Figure 1 As shown in Figure S1, this embodiment adopts a single-ended coherent detection system. The system includes a narrow-linewidth laser, a beam splitting unit, a pulse modulation module, an optical amplification module, an optical circulator, a sensing fiber, a local oscillator branch conditioning module, a coherent receiving module, a high-speed data acquisition module, and a digital signal processing module. Through a single-ended single-fiber phase-sensitive optical time-domain reflectometry system, a detection pulse is emitted into the sensing fiber at a fixed pulse repetition frequency. The Rayleigh backscattered echo in the fiber is received, and the in-phase component electrical signal and the quadrature component electrical signal are obtained through coherent reception. The in-phase component electrical signal and the quadrature component electrical signal are synchronously sampled and digitized to obtain the in-phase component data and the quadrature component data at each distance gate and each slow time sampling moment along the entire line.

[0028] like Figure 2 As shown in Figure S2, the spatial sampling unit after discretization along the optical fiber is defined as a distance gate, identified by index z; the pulse sequence dimension is defined as slow time, identified by index t, and the spatial position corresponding to the distance gate z is: , in, Given the system spatial sampling interval, based on the in-phase and quadrature component data obtained in step S1, the in-phase and quadrature components corresponding to each range gate and each slow time are synthesized into a complex echo sequence. For the t-th pulse and range gate z, the complex echo is represented as: , in, For in-phase components, For orthogonal components, Using the imaginary unit, based on the complex echo sequence, the phase principal value sequence of each range gate in the slow time dimension is extracted through four-quadrant arctangent operation to avoid ambiguity in quadrant determination using ordinary arctangent. For the i-th range gate and the m-th slow time sampling point, the complex echo is denoted as: , in, and These are the in-phase and quadrature components, respectively, thus yielding the phase principal value sequence of the range gate in the slow time dimension: , in, The phase principal value of the i-th distance gate at the m-th slow time sampling point. This is the arctangent operation for the four quadrants.

[0029] The phase is obtained directly from the arctangent of the four quadrants. They are usually distributed within the main value interval, and adjacent sampling points are separated by a certain distance. Abrupt transitions may occur at boundary times, therefore, it is necessary to expand the phase along the slow-time direction to obtain a continuous phase sequence. This embodiment adopts a point-by-point expansion method based on the principal value difference of the phase. Specifically, the principal value phase difference between adjacent slow-time sampling points of the i-th distance gate is calculated: , The principal phase difference is returned to the interval. Within this, the corrected principal value difference is obtained: when hour, ; when hour, ; when hour, ; Then, by accumulating the data point by point at a slower pace, we obtain the expanded continuous phase sequence: , in, Let be the continuous phase sequence after expansion of the i-th distance gate at the m-th slow time sampling point. To reflect the principal phase difference back to the interval The corrected principal value difference obtained within the inner, The continuous phase sequence after the expansion of the i-th distance gate at the (m-1)-th slow time sampling point can be eliminated by the above processing, thus making the phase change more realistically reflect the optical path change caused by vibration.

[0030] The unfolded continuous phase sequence typically contains both the true vibrational components and a low-frequency, slowly varying trend introduced by laser frequency shift, temperature changes, and receiver link drift. To suppress this type of slowly varying baseline, detrending processing is performed on the unfolded continuous phase sequence. Trend estimation is then performed on the unfolded continuous phase sequence to obtain the trend term. The trend term is removed from the expanded continuous phase sequence to obtain the detrended phase sequence: , in, Let be the detrended phase sequence of the i-th distance gate at the m-th slow time sampling point. This is a trend item. When balancing engineering implementation complexity and real-time performance, the trend item... A trend estimate can be obtained by using a moving average with a fixed window length; alternatively, it can be obtained by low-pass filtering, polynomial fitting, or local regression methods.

[0031] To further highlight the effective response within the target vibration frequency band and suppress high-frequency noise, bandpass filtering is performed on the detrended phase sequence to obtain a preprocessed phase sequence: , in, This is the preprocessed phase sequence of the i-th distance gate at the m-th slow time sampling point. It is a bandpass filter operator, and its passband range is set according to the frequency range of the target vibration.

[0032] After the above-mentioned phase extraction, phase expansion, detrending and bandpass filtering, the preprocessed phase sequence not only retains the main dynamic information caused by vibration propagation, but also reduces the impact of phase entanglement, low-frequency baseline fluctuation and high-frequency noise on subsequent propagation time delay estimation, thus providing a stable input for subsequent common drift compensation, multi-gate pair time delay fusion and direction determination.

[0033] S3. Under single-ended coherent detection conditions, laser frequency shift, ambient temperature changes, and receiver link phase drift introduce a slowly varying phase term that is approximately common along the optical fiber, causing synchronous shifts in the phase sequences of each range gate. This results in drift or spurious peaks in subsequent arrival delay estimations based on inter-gate correlation. For ease of explanation, the preprocessed phase sequence obtained in step S2 is represented as follows: , in, This is the preprocessed phase sequence of the distance gate z at slow time t. The effective component caused by the propagation of real vibrations. For common drift items, It is noise.

[0034] like Figure 3 As shown, a protection distance is set outside the target event area. The range gates located outside the protection distance and unaffected by events are defined as the reference range gate set, denoted as . The protection distance This is used to ensure that the reference distance gate set does not include the event-affected area, thereby avoiding interference from the actual vibration signal on the common drift estimation. The protection distance is determined comprehensively based on the system spatial resolution, the event area positioning error, the differential distance used for phase demodulation, and the maximum spatial range that vibration may affect within the decision time window. Specifically, the protection distance is not less than the maximum value among the following four: the length corresponding to the system spatial resolution, the differential distance or gauge length used for phase demodulation, the upper limit of the positioning error of the target event area, and the maximum estimated distance of vibration propagation within the decision time window, with a safety margin reserved on this basis to avoid contamination of the reference distance gate by the vibration response, pulse broadening, and positioning error of the target event. In an optional embodiment, the protection distance... Determined according to the following formula: , in, For system spatial resolution, The differential distance or gauge length used for phase demodulation. This represents the upper limit of the positioning error for the target event area. The maximum estimated velocity of the vibration to be detected propagating along the sensing optical fiber. The current determination time window length. For safety margin.

[0035] The protection distance can also be verified by the stability of the reference distance gate. Specifically, the protection distance is gradually increased, and the energy, median absolute deviation, and correlation with the target event area signal of the phase sequence within the reference distance gate set are statistically analyzed. When the vibration energy of the reference distance gate is lower than a preset background energy threshold, and its correlation with the target event area signal is lower than a preset pollution threshold, the corresponding distance is determined as the protection distance. This ensures that the reference distance gate set mainly includes common drift components and background noise, but not target vibration propagation components.

[0036] For each slow time t, the reference distance gate set Robust statistics are taken from the preprocessed phase sequences of each range gate to estimate the common drift term. Preferably, the median is used as the robust statistics. , in, This is an estimate of the common drift term at slow time t. To perform median calculations on all range gates within the reference range gate set, the estimated common drift term is finally removed from the preprocessed phase sequence of each range gate, resulting in the compensated phase sequence: , in, The compensated phase sequence of the distance gate z at slow time t is obtained by suppressing the common drift term through the above self-calibration compensation, so that the subsequent arrival delay estimation is mainly determined by the propagation-related component, thereby improving the long-term stability and repeatability of the direction determination.

[0037] S4. Based on the event detection results, energy mutation distribution, or preset judgment region, determine a candidate range gate set, denoted as C, within the target event region. Due to factors such as Rayleigh scattering fading, local low signal-to-noise ratio, polarization state changes, and phase expansion instability, the candidate set C may contain anomalous range gates unsuitable for direction determination. Directly using all candidate range gates for gate pair construction and time delay estimation can easily lead to enhanced cross-correlation spurious peaks, time delay estimation shifts, and final direction misjudgment. Therefore, a quality index is calculated for each range gate within the candidate range gate set C, and an effective range gate set V is constructed according to the screening rules.

[0038] In this embodiment, the quality indicators include at least one or more of the following categories: amplitude quality indicators, amplitude stability indicators, phase availability indicators, and related availability indicators. For the i-th candidate distance gate, within the decision time window W, its echo amplitude sequence is: , in, Let be the echo amplitude of the i-th candidate distance gate at the m-th sampling point. For in-phase components, These are orthogonal components. The median is taken as the amplitude level indicator. , in, Let be the amplitude quality index for the i-th candidate distance gate. To calculate the median of all sampling points within the decision time window W, when Below the preset lower limit At that time, it is considered that the range gate is in a Rayleigh scattering fading state or the reception quality is insufficient.

[0039] The amplitude stability index is: , in, Let be the amplitude stability index for the i-th candidate distance gate. For standard deviation calculation, For mean calculation, It is a minimal constant, when Greater than the preset stability threshold This indicates that the distance gate fluctuates significantly within the current time window, resulting in unstable reception.

[0040] For the expanded phase sequence, we can count the number of abnormal jumps, the number of expansion failures, or the proportion of excessive phase abrupt changes between adjacent sampling points within the decision time window. For example, we can define the abnormal jump count as... And define the phase availability index as: , in, Let be the phase availability index for the i-th candidate distance gate. Indicates the decision time window The total number of sampling points within, For the i-th candidate distance gate in the decision time window The abnormal jump count within, when Below the preset threshold If the phase unfolding reliability of the distance gate is insufficient, it can be considered that the gate does not participate in the subsequent gate pair construction.

[0041] Cross-correlation analysis is performed on the candidate distance gate with its adjacent gates, reference template, or other distance gates that passed the initial screening, and the significance of the main peak is calculated. Preferably, the significance ratio of the main peak can be defined as: , in, For the i-th candidate distance gate, the relevant availability metric is... Indicates the amplitude of the main cross-correlation peak. This represents the amplitude of the secondary peak or the representative value of the sidelobe. To prevent constants with a denominator of zero. When Less than the preset threshold If the peak corresponding to the distance gate is not prominent, it indicates that the ability to represent the actual propagation delay is weak, and it can be removed.

[0042] Based on the above quality indicators, this embodiment adopts the following screening logic: when the candidate distance gates simultaneously meet the following conditions... , , , If the candidate distance gate is selected, it is included in the effective distance gate set V; otherwise, it is removed from the effective distance gate set. The quality index threshold is determined as follows: The preset thresholds for each quality indicator can be determined during the system installation and commissioning phase based on background data under targetless vibration conditions, calibration data under known vibration conditions, and preset misjudgment rate requirements. These thresholds can be adaptively updated during system operation based on recent background data. For amplitude quality indicators, the lower amplitude threshold is used to exclude range gates in Rayleigh scattering fading states. This threshold can be determined based on the statistical distribution of echo amplitudes from each range gate under targetless vibration conditions. For example, the lower quantile of the amplitude quality indicators for reference range gates or historical normal range gates can be used as the lower amplitude threshold; or it can be determined according to the following formula: , in, and These are the mean and standard deviation of the quality index for the normal distance gate amplitude, respectively. This is a coefficient set according to the allowed rejection ratio. When the amplitude quality index of the candidate distance gate is lower than... When this occurs, it is classified as a low-amplitude or fading distance gate.

[0043] For the amplitude stability index, its upper threshold is used to exclude distance gates that exhibit excessive amplitude fluctuations within the judgment time window. This threshold can be determined based on the higher quantile of the amplitude stability index under targetless vibration conditions, or according to the following formula: , in, and These represent the mean and standard deviation of the amplitude stability index for a normal distance gate, respectively. When the amplitude stability index of a candidate distance gate is higher than... At that time, it was considered that its receiving status was unstable.

[0044] For the phase availability index, its lower limit threshold is determined based on the allowable proportion of abnormal phase jumps. Let the maximum allowable proportion of abnormal jumps within the judgment time window be... Then the phase availability threshold can be expressed as: , When the phase availability index of the candidate distance gate is lower than If the number of phase expansion failures or abnormal transitions exceeds the allowable range, it will not participate in the subsequent gate pair construction.

[0045] For relevant availability metrics, the lower limit threshold is determined based on calibration data with known propagation direction and known arrival delay. The availability metrics for correct and incorrect gate pairs are calculated separately, and the relevant availability threshold is determined based on the preset maximum incorrect gate pair reception rate. In one alternative implementation, the relevant availability threshold is determined according to the following formula: , in, is any candidate relevant availability threshold from the candidate threshold set. This indicates that the relevant availability threshold will be set to... At that time, the proportion of incorrect gate pairs that are incorrectly included in the set of valid gate pairs. The preset maximum error gate pair reception rate. When the relevant availability metrics of the candidate distance gate... Below At that time, it was believed that the main peak of its cross-correlation was not prominent enough and did not participate in the subsequent gate pair construction.

[0046] During long-term system operation, background statistics for each quality indicator can be continuously updated during periods when no target event occurs. The corresponding thresholds are updated using a sliding time window or exponential weighting method to adapt to long-term state changes caused by variations in optical power, ambient temperature, and device aging. To prevent sudden threshold changes due to short-term anomalies, upper and lower limits can be set for the threshold update magnitude.

[0047] For different application scenarios, a combination of two, three, or all of the quality indicators can be used for screening as needed. Through the above quality assessment and screening, range gates with sufficient amplitude, stable state, reliable phase, and significant correlation peaks can be preferentially retained within the event region. This ensures that subsequent propagation delay estimation and direction determination no longer rely on a few failed or fading gates, reducing the probability of misjudgment under Rayleigh scattering fading and polarization perturbation conditions, and improving robustness and repeatability in engineering applications.

[0048] After obtaining the effective range gate set V, construct a set of forward gate pairs within the effective range gate set V. Each forward gate pair consists of two effective distance gates, denoted as . ,in ,and The door spacing between the two doors meets the preset constraints: , In the formula, To preset the minimum door spacing, The maximum gate spacing is preset. The purpose of setting gate spacing constraints is to avoid excessively small propagation delays when gates are too close, which would be detrimental to stable estimation, and to avoid excessive waveform differences when gates are too far apart, which would reduce the reliability of cross-correlation. This yields the set of forward gate pairs: , in, For the forward gate pair set, For the effective distance gate set, and For effective distance gate index, To preset the minimum door spacing, This is the preset maximum door spacing.

[0049] like Figure 4 As shown, S5, for the set of forward gate pairs Any forward gate pair in Within the decision time window W, the preprocessed phase sequence after compensation obtained in step S3 is used. and Calculate the discrete cross-correlation function: , in, For the door Discrete cross-correlation function, For discrete delay quantities, and These are the compensated preprocessed phase sequences for range gate i and range gate j, respectively. The delay corresponding to the maximum value of the cross-correlation function is used as the arrival delay estimate for this gate pair. , in, For the door The discrete delay estimate, if the sampling time interval is The corresponding actual time delay is: , in, For the door The actual time delay is significant because the correlation quality of different gate pairs may vary, and a few failed or low-quality gate pairs can easily affect the overall result. Therefore, this embodiment further constructs weights based on the significance of the correlation peaks of each gate pair. Preferably, the weights can be defined as follows: The weights of each door pair are: , in, For the door The weight, This represents the amplitude of the main peak of the cross-correlation of the gate. This is the amplitude of the secondary peak or the representative value of the sidelobe. This weight reflects the prominence of the main peak relative to the secondary peak or sidelobe. The clearer the main peak and the more reliable the time delay estimate, the greater the corresponding weight.

[0050] After obtaining the delay estimates and corresponding weights for all forward gate pairs, a robust fusion method is used to obtain the forward fusion delay. Preferably, the weighted median is used to fuse the actual time delays of all gate pairs in the forward gate pair set. Specifically, the forward fusion delay is obtained by fusing the delays of each forward gate pair according to their weights. , in, For forward fusion latency, As a weighted median operator, it can be similarly obtained that the reverse fusion delay is obtained by fusing the delays of each reverse gate pair according to the weight. Compared with the simple weighted average, the weighted median is less sensitive to outliers and abnormal gate pairs, and is more suitable for situations where there are fading gates, spurious peak gates and short-term abnormal disturbances in the engineering field.

[0051] Considering the potential systematic biases in real-world systems, such as differences in equivalent bandwidth, group delay, and fixed processing link delay between gates of different distances, a set of forward gate pairs is constructed. Corresponding reverse door pair set Each gate pair in the reverse gate pair set consists of distance gates from the effective distance gate set V, and satisfies the same gate spacing constraint as the forward gate pair. The reverse fusion delay is obtained by calculating the discrete cross-correlation function based on the compensated preprocessed phase sequence, constructing gate pair weights, and using weighted median fusion, in the same manner as the forward gate pair. .

[0052] The unbiased delay is obtained by taking the difference between the forward fusion delay and the backward fusion delay and then taking half of the difference: , in, To eliminate the unbiased time delay after the system is fixed biased, this symmetrical differential method can reduce the impact of inter-gate response inconsistency and fixed link delay on the final direction determination without additional external calibration, so that the direction determination results have higher long-term stability and repeatability.

[0053] S6. Obtain the unbiased delay obtained in step S5. After that, the direction is not output immediately, but its absolute value is first compared with a preset direction determination threshold. Comparison; the preset direction determination threshold Based on the system's time resolution, gate spacing, and expected propagation speed range, it is used to distinguish between effective propagation delay and minor deviations caused by noise or estimation fluctuations.

[0054] like Figure 5 As shown, when the absolute value of the unbiased delay is less than or equal to the preset direction determination threshold, i.e. If the current event direction information is insufficient or the multi-gate estimation consistency is insufficient, an uncertain state is output to avoid giving incorrect conclusions when the direction information is unreliable; when the unbiased delay is greater than the preset direction determination threshold, i.e. Based on The sign determines the direction of propagation: when When, it is determined that the vibration propagates from the near end to the far end; when When the vibration is detected, it is determined that the vibration propagates from the far end to the near end.

[0055] To characterize the reliability of direction determination, a consistency index is output synchronously. This consistency index is formed by the degree of consistency between the gate pair estimates and the fused value, calculated based on the deviation between the arrival delay estimates and the fused delay of each gate pair, combined with the gate pair weights. , in, As a consistency indicator, For a set of gate pairs, Let (a, b) be the weight of the gate pair. Let be the estimated arrival time delay of gate pair (a,b). This is the estimated delay value after fusion. This is the consistency tolerance threshold. This is an indicator function; it takes a value of 1 when the condition within the parentheses is true, and a value of 0 otherwise. The larger the consistency index C, the more consistent the estimation results of multiple pairs, and the higher the reliability of the corresponding direction judgment. This consistency index is used for subsequent alarm decisions, result screening, or manual review, providing a quantitative basis for reliability in the engineering field.

[0056] Example 2 The difference between this embodiment and Embodiment 1 is that this embodiment provides a vibration propagation direction determination system based on phase-sensitive optical time-domain reflection; The system in this embodiment adopts a single-end access coherent detection system, including a narrow linewidth laser, a beam splitting unit, a pulse modulation module, an optical amplification module, an optical circulator, a sensing fiber, a local oscillator branch conditioning module, a coherent receiving module, a high-speed data acquisition module, and a digital signal processing module.

[0057] The output of the narrow-linewidth laser is connected to the input of a beam splitter, which divides the laser output into a transmitter branch and a local oscillator branch. The transmitter branch generates the probe pulse, while the local oscillator branch provides the local oscillator light required for coherent reception. The transmitter branch is sequentially connected to a pulse modulation module, an optical amplification module, and the first port of an optical circulator. The second port of the optical circulator is connected to a sensing fiber to inject the probe pulse. The local oscillator branch is connected to the local oscillator input of the coherent receiver module via a local oscillator conditioning module. The third port of the optical circulator is connected to the signal input of the coherent receiver module, and the output of the coherent receiver module is connected to the input of a high-speed data acquisition module. The output of the high-speed data acquisition module is connected to the input of a digital signal processing module.

[0058] The pulse modulation module modulates the continuous laser into a probe pulse with a predetermined pulse width and repetition frequency. The optical amplification module increases the probe pulse power to meet the power requirements for propagation along the sensing fiber and for exciting Rayleigh backscattering. The local oscillator branch conditioning module includes one or more of a polarization adjuster, a variable optical attenuator, or a delay adjustment unit to adjust the polarization state, power, and time delay matching of the local oscillator light, ensuring optimal interference conditions between the local oscillator light and the echo signal light in the coherent receiving module, thereby improving the interference reception quality. The coherent receiving module causes the echo signal light to interfere with the local oscillator light and outputs corresponding in-phase and quadrature component electrical signals. The high-speed data acquisition module synchronously samples and digitizes the in-phase and quadrature component electrical signals to obtain digital values ​​at each distance gate and each slow-time sampling moment along the entire line.

[0059] When the probe pulse enters the sensing fiber, the Rayleigh backscattered echoes generated at various locations within the fiber return via an optical circulator and are output from the third port of the circulator to the signal input of the coherent receiver module. Simultaneously, the local oscillator branch, split by the beam splitter, is input to the local oscillator input of the coherent receiver module after passing through a local oscillator branch conditioning module. The echo signal light interferes with the local oscillator light in the coherent receiver module, producing in-phase and quadrature component electrical signals. These in-phase and quadrature component electrical signals are synchronously sampled and digitized by the high-speed data acquisition module before being sent to the digital signal processing module.

[0060] The digital signal processing module is used to process the digital signals output by the high-speed data acquisition module to determine the direction of vibration propagation along the axis of the sensing fiber. Specifically, the digital signal processing module constructs complex echo sequences for each range gate and extracts phase sequences. The phase sequences are then subjected to phase expansion, detrending processing, and bandpass filtering to obtain a preprocessed phase sequence. Outside the target event region, a reference range gate set is selected. Based on this set, a common drift term is estimated, and the preprocessed phase sequence is self-calibrated to obtain a compensated phase sequence. Within the target event region, candidate range gates are screened based on quality indicators to obtain an effective range gate set. Within this effective range gate set, a forward gate pair set satisfying a preset gate spacing constraint is constructed. Arrival delays for each gate pair in the forward gate pair set are estimated, and gate pair weights are constructed. Based on these weights, robust fusion of the arrival delays of each gate pair is performed to obtain the forward fusion delay. A reverse gate pair set corresponding to the forward gate pair set is constructed, and the reverse fusion delay is obtained. Symmetrical difference is performed between the forward and reverse fusion delays to obtain the unbiased delay. The unbiased delay is compared with a preset direction determination threshold. Based on the comparison result, the propagation direction of vibration along the sensing fiber axis is output, along with a consistency index reflecting the degree of consistency in the estimation of multiple gate pairs.

[0061] Thus, the system forms the following signal link: the output light from the narrow-linewidth laser is split by a beam splitter; one path is modulated by a pulse modulation module, amplified by an optical amplification module, and then injected into the sensing fiber through an optical circulator; the other path, as the local oscillator light, is conditioned by a local oscillator branch conditioning module and then sent to a coherent receiving module; the Rayleigh backscattered echo in the sensing fiber is extracted by the optical circulator and interferes with the local oscillator light in the coherent receiving module, obtaining in-phase and quadrature component electrical signals; these in-phase and quadrature component electrical signals are synchronously sampled and digitized by a high-speed data acquisition module and then sent to a digital signal processing module to complete subsequent direction determination processing. The above system connection relationship and signal link ensure stable determination of vibration propagation direction under single-end access and a single sensing fiber, without the need for double-end deployment or additional sensors, and improves the stability and repeatability of determination under adverse engineering conditions such as laser frequency shift, temperature drift, Rayleigh scattering fading, and polarization state changes.

[0062] The above are all preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Therefore, all equivalent changes made in accordance with the structure, shape and principle of the present invention should be covered within the scope of protection of the present invention.

Claims

1. A method for determining the vibration propagation direction of phase-sensitive optical time-domain reflectometry, characterized in that, include: Acquire the in-phase and quadrature components of the full-line echo collected by a single-ended single-fiber phase-sensitive optical time-domain reflectometry system; Construct complex echo sequences for each range gate, extract the phase sequences for each range gate in the slow time dimension, and perform preprocessing operations on the phase sequences to obtain preprocessed phase sequences; A set of reference range gates is selected outside the target event region. The common drift term is estimated based on the set of reference range gates, and the preprocessed phase sequence is self-calibrated and compensated to obtain the compensated phase sequence. Within the target event area, candidate distance gates are screened based on quality indicators to obtain an effective distance gate set, and a set of forward gate pairs that satisfy the preset gate spacing constraint is constructed within the effective distance gate set; Arrival delay is estimated for each gate pair in the forward gate pair set and gate pair weights are constructed. Based on the gate pair weights, the arrival delays of each gate pair are robustly fused to obtain the forward fused delay. Construct a set of reverse gate pairs corresponding to the set of forward gate pairs and obtain the reverse fusion delay. Perform symmetric difference between the forward fusion delay and the reverse fusion delay to obtain the unbiased delay. The unbiased time delay is compared with a preset direction determination threshold. Based on the comparison result, the propagation direction of vibration along the fiber axis is output, and a consistency index reflecting the degree of consistency of the multi-gate pair estimation is also output.

2. The method for determining the vibration propagation direction of phase-sensitive optical time-domain reflection according to claim 1, characterized in that, The step of extracting the phase sequence of each distance gate in the slow time dimension includes defining an index identifier, defining the spatial sampling unit after discretization along the optical fiber as a distance gate according to the index identifier, and defining the pulse sequence dimension as slow time. The in-phase and quadrature components corresponding to each distance gate and each slow time are synthesized into a complex echo sequence. Based on the complex echo sequence, the phase principal value sequence is extracted through a four-quadrant arctangent operation. The expression for the phase principal value sequence is: ， in, The phase principal value of the i-th distance gate at the m-th slow time sampling point. For orthogonal components, For in-phase components, This is the arctangent operation for the four quadrants.

3. The method for determining the vibration propagation direction of phase-sensitive optical time-domain reflection according to claim 2, characterized in that, The preprocessing operation of the phase sequence includes phase expansion of the phase principal value sequence, detrending of the expanded continuous phase sequence to obtain a detrended phase sequence, and finally bandpass filtering of the detrended phase sequence to obtain a preprocessed phase sequence. The phase unrolling includes: calculating the principal phase difference of the i-th distance gate between adjacent slow-time sampling points, and then flipping the principal phase difference back to the interval. The corrected principal value difference is obtained internally, and then accumulated point by point in slow time to obtain the expanded continuous phase sequence: ， in, Let be the continuous phase sequence after expansion of the i-th distance gate at the m-th slow time sampling point. To reflect the principal phase difference back to the interval The corrected principal value difference obtained within the inner, It is the continuous phase sequence after the expansion of the i-th distance gate at the (m-1)-th slow time sampling point.

4. The method for determining the vibration propagation direction of phase-sensitive optical time-domain reflection according to claim 1, characterized in that, The self-calibration compensation of the preprocessed phase sequence includes setting a protection distance outside the target event region and determining the distance gates located outside the protection distance as the reference distance gate set; For each slow time, the median of the preprocessed phase sequence of each range gate in the reference range gate set is taken to estimate the common drift term; The common drift term is removed from the preprocessed phase sequence of each range gate to obtain the compensated phase sequence: ， in, The phase sequence after compensation for the distance gate z at slow time t. This is the preprocessed phase sequence of the distance gate z at slow time t. This is the estimate of the common drift term at slow time t.

5. The method for determining the vibration propagation direction of phase-sensitive optical time-domain reflection according to claim 1, characterized in that, The process of filtering candidate distance gates based on quality indicators includes determining a set of candidate distance gates within the target event region; Calculate the quality index of each candidate distance gate, wherein the quality index includes at least one of amplitude quality index, amplitude stability index, phase availability index and related availability index; The quality indicators are compared with the corresponding preset thresholds. Candidate distance gates that meet the screening conditions are included in the effective distance gate set, while candidate distance gates that do not meet the screening conditions are removed from the effective distance gate set. The relevant availability metrics are calculated using the following formula: ， in, For the i-th candidate distance gate, the relevant availability metric is... The amplitude of the main peak of cross-correlation, This is the amplitude of the secondary peak. To prevent extremely small constants with a denominator of zero.

6. The method for determining the vibration propagation direction of phase-sensitive optical time-domain reflection according to claim 5, characterized in that, The construction of the forward door pair set that satisfies the preset door spacing constraint includes selecting forward door pairs consisting of two effective distance doors from the effective distance door set, wherein each forward door pair is denoted as (i,j), and The door spacing between the two doors meets the preset constraints. ; All forward door pairs that satisfy the door spacing constraint are defined as the set of forward door pairs: ， in, For the forward gate pair set, For the effective distance gate set, and For effective distance gate index, To preset the minimum door spacing, This is the preset maximum door spacing.

7. The method for determining the vibration propagation direction of phase-sensitive optical time-domain reflection according to claim 1, characterized in that, The symmetrical difference between the forward fusion delay and the reverse fusion delay includes calculating the discrete cross-correlation function for each gate pair in the forward gate pair set based on the compensated preprocessed phase sequence, and taking the delay corresponding to the maximum value of the cross-correlation function as the estimated arrival delay value of the gate pair. Gate pair weights are constructed based on the ratio of the main peak amplitude to the secondary peak amplitude of the cross-correlation of each gate pair, and the arrival delay estimates of each gate pair are weighted and fused by median to obtain the forward fusion delay. Construct a set of reverse door pairs corresponding to the set of forward door pairs, wherein each door pair in the set of reverse door pairs satisfies the same door spacing constraint as the set of forward door pairs; The forward fusion delay is obtained by differentiating the backward fusion delay and taking half of the result. The expression for the forward fusion delay is as follows: ， in, For forward fusion latency, Let be the estimated arrival delay of gate pair (i,j). Let (i,j) be the weight of the gate pair. For the forward gate pair set, This is the weighted median operator.

8. The method for determining the vibration propagation direction of phase-sensitive optical time-domain reflection according to claim 1, characterized in that, The output reflects the consistency index of the multi-gate pair estimation consistency, including comparing the absolute value of the unbiased delay with a preset direction determination threshold. When the absolute value of the unbiased delay is less than or equal to the preset direction determination threshold, an uncertain state is output. When the unbiased delay is greater than the preset direction determination threshold, it is determined that the vibration propagates from the near end to the far end. When the unbiased delay is less than the negative preset direction determination threshold, it is determined that the vibration propagates from the far end to the near end. Based on the deviation between the estimated arrival delay of each gate pair and the fusion delay, and combined with the gate pair weights, a consistency index is calculated: ， in, As a consistency indicator, For a set of gate pairs, Let (a, b) be the weight of the gate pair. Let be the estimated arrival time delay for gate pair (a,b). This is the estimated latency value after fusion. This is the consistency tolerance threshold. This is an indicator function.

9. A vibration propagation direction determination system based on phase-sensitive optical time-domain reflectometry, executed according to claim 1, characterized in that, include: Narrow linewidth laser, beam splitting unit, local oscillator branch conditioning module, coherent receiving module, high-speed data acquisition module and digital signal processing module; The narrow linewidth laser is connected to the beam splitting unit, which divides the laser into an emission branch and a local oscillator branch. The emission branch is connected to the sensing fiber via a pulse modulation module, an optical amplification module, and the first port of an optical circulator in sequence. The local oscillator branch is connected to the local oscillator input terminal of the coherent receiving module via the local oscillator branch conditioning module. The third port of the optical circulator is connected to the signal input terminal of the coherent receiving module. The output terminal of the coherent receiving module is connected to the input terminal of the high-speed data acquisition module. The output terminal of the high-speed data acquisition module is connected to the input terminal of the digital signal processing module.

10. The vibration propagation direction determination system for phase-sensitive optical time-domain reflectometry according to claim 9, characterized in that, The digital signal processing module is used to process the digital signal output by the high-speed data acquisition module to determine the propagation direction of vibration along the axis of the sensing fiber.

Images