Wavelength deviation compensation method and system for fast wavelength demodulation of a remote fiber grating

Abstract

The application provides a wavelength deviation compensation method and system for remote fiber grating wavelength rapid analysis, relates to the technical field of fiber grating sensing, and comprises the following steps: collecting grating reflection spectrum signals of different transmission distances, performing frequency domain transformation to separate phase components, establishing a mapping relationship between transmission distance and dispersion accumulation, and calculating an initial wavelength offset estimation value; identifying a refractive index section, and establishing a segmented compensation coefficient; constructing a wavelength-strain-temperature three-dimensional space, optimizing the compensation coefficient according to trajectory curvature, and updating the compensation coefficient according to residual gradient. The application can improve the wavelength analysis accuracy of remote fiber grating and reduce the influence of dispersion on measurement.

Description

Technical Field

[0001] This invention relates to the field of fiber optic grating sensing technology, and in particular to a wavelength deviation compensation method and system for rapid wavelength resolution of far-end fiber optic gratings. Background Technology

[0002] Fiber Bragg grating (FBG) sensing technology is a technique that utilizes the characteristic that the center wavelength of a fiber Bragg grating varies with external physical quantities (such as strain and temperature) to detect physical parameters. Remote FBG sensing systems allow sensing elements to be placed far from demodulation equipment, transmitting reflected spectral signals over long distances via optical fibers to achieve real-time monitoring of remote environmental parameters.

[0003] During long-distance fiber optic transmission, the optical signal experiences wavelength-dependent delays due to fiber dispersion, resulting in distortion of the received grating reflection spectrum. This distortion causes significant errors in traditional wavelength resolution algorithms (such as the centroid method and Gaussian fitting method) at the far end of the grating, affecting measurement accuracy. This is especially true in multi-grating cascade systems, where gratings at different locations are affected by dispersion to varying degrees, leading to different wavelength resolution errors.

[0004] Currently, long-range fiber optic grating wavelength resolution technology still has some defects and shortcomings. Traditional wavelength resolution methods fail to effectively distinguish between intrinsic wavelength shift of the grating and wavelength shift caused by transmission dispersion. In long-distance transmission scenarios, the cumulative effect of dispersion causes the resolution accuracy to decrease significantly with increasing transmission distance. Existing dispersion compensation techniques are mainly based on uniform dispersion models, which cannot adapt to the refractive index abrupt changes in optical fibers, such as local dispersion changes caused by non-uniform factors like fiber joints and bends, resulting in poor compensation effects. Existing technologies lack adaptive wavelength deviation compensation mechanisms. Under drastic changes in environmental conditions, fixed compensation parameters are difficult to track real-time dispersion changes, especially the fluctuations in fiber dispersion parameters caused by temperature changes, which makes the measurement accuracy unstable during long-term monitoring. Summary of the Invention

[0005] This invention provides a wavelength deviation compensation method and system for rapid wavelength resolution of far-end fiber Bragg gratings, which can solve the problems in the prior art.

[0006] A first aspect of this invention provides a wavelength deviation compensation method for rapid wavelength resolution of a far-end fiber Bragg grating, comprising:

[0007] The reflectance spectral signals of multiple fiber gratings that have traveled different transmission distances were collected;

[0008] The reflected spectral signal is frequency domain transformed to separate the linear phase component and the nonlinear phase component, a linear mapping relationship between transmission distance and dispersion accumulation is established, the transmission dispersion offset and intrinsic wavelength offset are calculated, and the initial wavelength offset estimate of the grating at different spatial positions is obtained.

[0009] Based on the initial wavelength offset estimate, the refractive index abrupt change segment and the refractive index stable segment are identified, and a segmented compensation coefficient is established. The correction amount is calculated through tensor decomposition to correct the initial wavelength offset estimate and obtain the compensated grating center wavelength.

[0010] A three-dimensional space of wavelength-strain-temperature is constructed, and the compensated grating center wavelength is mapped to a trajectory in the three-dimensional space. The trajectory curvature is calculated. When the trajectory curvature exceeds the preset curvature threshold, the piecewise compensation coefficient is reconstructed based on the historical initial wavelength offset estimation value sequence, and the physical quantity demodulation data is output.

[0011] The wavelength reconstruction residual in the demodulated physical quantity data is mapped to the three-dimensional space, and the piecewise compensation coefficient is updated according to the residual gradient direction.

[0012] In one optional embodiment, the reflected spectral signal is subjected to frequency domain transformation to separate the linear phase component and the nonlinear phase component, a linear mapping relationship between transmission distance and dispersion accumulation is established, and the transmission dispersion offset and intrinsic wavelength offset are calculated to obtain the initial wavelength offset estimate of the grating at different spatial locations, including:

[0013] The reflection spectrum signal is transformed from the time domain to the frequency domain to obtain the frequency domain spectral distribution;

[0014] A continuous phase response curve is extracted from the frequency domain spectral distribution. The phase response curve is then fitted with a polynomial according to the fitting weights to separate the linear phase component composed of the linear fitting parameters and the nonlinear phase component composed of the residual phase distribution.

[0015] Calculate the group delay difference corresponding to different transmission distances based on the linear phase components, and establish a linear mapping relationship between transmission distance and dispersion accumulation.

[0016] The nonlinear phase component is converted back to the wavelength domain to extract the intrinsic center wavelength of each grating position.

[0017] Based on the linear mapping relationship, the group delay difference is converted into the dispersion accumulation at each spatial location, and the transmission dispersion offset and intrinsic wavelength offset of each grating are calculated in combination with the intrinsic center wavelength.

[0018] By vector synthesis of the transmission dispersion offset and the intrinsic wavelength offset, the initial wavelength offset estimate of the grating at different spatial positions is obtained.

[0019] In one optional embodiment, extracting a continuous phase response curve from the frequency domain spectral distribution, performing polynomial fitting on the phase response curve according to the fitting weights, and separating the linear phase component composed of the linear fitting parameters and the nonlinear phase component composed of the residual phase distribution includes:

[0020] Phase unwrapping processing is performed on the complex amplitude at each frequency point in the frequency domain spectral distribution to generate a continuous phase response curve;

[0021] Each sampling point in the continuous phase response curve is assigned a fitting weight according to the amplitude information. The continuous phase response curve is then fitted with a polynomial according to the fitting weight to separate the linear fitting parameters that characterize the first power relationship of frequency and the nonlinear fitting parameters that characterize the higher power relationship of frequency.

[0022] The linear fitting parameters are multiplied by each frequency point to obtain the linear phase component corresponding to the transmission distance. The residual phase distribution is obtained by subtracting the linear phase component from the continuous phase response curve.

[0023] The ratio of the phase difference to the frequency interval of adjacent frequency points in the residual phase distribution is calculated to generate a phase change rate sequence. The direction of grating period change is determined by statistically analyzing the positive and negative values ​​of the phase change rate sequence. The grating period gradient rate is extracted by the adjacent difference values ​​of the phase change rate sequence. The structural features of the residual phase distribution are calibrated based on the direction of grating period change and the grating period gradient rate.

[0024] The residual phase distribution calibrated by structural features is determined as the nonlinear phase component corresponding to the grating structure.

[0025] In one optional embodiment, refractive index abrupt change segments and refractive index stable segments are identified based on the initial wavelength offset estimate, and piecewise compensation coefficients are established. The correction amount is calculated through tensor decomposition to correct the initial wavelength offset estimate, resulting in the compensated grating center wavelength, including:

[0026] The initial wavelength offset estimate is segmented and statistically analyzed according to the transmission distance. The gradient rate of change of the initial wavelength offset estimate in each transmission distance segment is calculated. Transmission distance segments with a gradient rate of change exceeding a preset rate of change threshold are marked as refractive index abrupt change segments, and transmission distance segments with a gradient rate of change within the rate of change threshold are marked as refractive index stable segments.

[0027] Extract the boundary positions between the refractive index abrupt change region and the refractive index stable region, and establish piecewise compensation coefficients at the boundary positions;

[0028] Principal component analysis is performed on the temperature field components and strain field components to determine the tensor dimension size. A third-order tensor structure is constructed and filled. Constrained alternating least squares iterative operation is performed to separate the first eigenvector of the transmission distance dimension, the second eigenvector of the temperature field dimension, and the third eigenvector of the strain field dimension.

[0029] The first eigenvector is multiplied element-wise with the piecewise compensation coefficient to obtain the transmission dispersion correction; the second eigenvector and the third eigenvector are multiplied by tensor cross product to obtain the multiphysics coupling correction.

[0030] The compensated grating center wavelength is obtained by subtracting the transmission dispersion correction and multiphysics coupling correction from the initial wavelength offset estimate.

[0031] In one optional embodiment, principal component analysis is performed on the temperature field components and strain field components to determine the tensor dimension size, a third-order tensor structure is constructed and filled, and a constrained alternating least squares iterative operation is performed to separate the first eigenvector of the transmission distance dimension, the second eigenvector of the temperature field dimension, and the third eigenvector of the strain field dimension, including:

[0032] Principal component analysis is performed on the temperature field component and the strain field component respectively. The variance contribution rate of each principal component is extracted, the cumulative sum of the variance contribution rates of each principal component is calculated, the cumulative variance contribution rate is determined, and the number of principal components corresponding to the cumulative variance contribution rate reaching the preset variance threshold is determined as the temperature field dimension and the strain field dimension.

[0033] A third-order tensor structure is constructed based on the number of sampling points, the dimension of the temperature field, and the dimension of the strain field of the initial wavelength offset estimate.

[0034] The initial wavelength offset estimate is traversed through the three-dimensional coordinate indices of the third-order tensor structure. When there is no measured data for the tensor element corresponding to the three-dimensional coordinate index, the filled tensor element whose Euclidean distance to the three-dimensional coordinate index is less than a preset distance threshold is searched. The inverse of the Euclidean distance of each filled tensor element is calculated to determine the interpolation weight. The filled tensor elements are weighted and summed according to the interpolation weight, and the three-dimensional coordinate index position is filled to obtain the filled third-order tensor structure.

[0035] The filled third-order tensor structure is subjected to alternating least squares iterative operations with nonnegativity and orthogonality constraints. The iteration is terminated when the relative rate of change of the third-order tensor reconstruction error is continuously less than the preset convergence threshold. The factor matrices of each dimension are extracted to determine the first eigenvector, the second eigenvector, and the third eigenvector.

[0036] In one optional embodiment, a wavelength-strain-temperature three-dimensional space is constructed, and the compensated grating center wavelength is mapped to a trajectory in the three-dimensional space. The trajectory curvature is calculated, and when the trajectory curvature exceeds a preset curvature threshold, the piecewise compensation coefficients are reconstructed based on the historical initial wavelength offset estimation sequence. The output physical quantity demodulation data includes:

[0037] A three-dimensional space with wavelength, strain, and temperature as coordinate axes is established. The center wavelength of the compensated grating and the strain and temperature values ​​at the corresponding time are combined to form a three-dimensional coordinate point sequence, thus forming a trajectory in the three-dimensional space.

[0038] The trajectory tangent vector is obtained by performing differential operations on the sequence of three-dimensional coordinate points, and the trajectory curvature is determined by calculating the rate of change of the angle between adjacent trajectory tangent vectors.

[0039] When the trajectory curvature exceeds a preset curvature threshold, identify the wavelength shift pattern corresponding to the period of abnormal curvature in the historical initial wavelength shift estimation value sequence, and extract the frequency domain feature components of the wavelength shift pattern.

[0040] The initial wavelength offset estimate is matched with the frequency domain feature components, and compensation weights for different refractive index segments are assigned according to the matching degree to reconstruct the segmented compensation coefficients.

[0041] The reconstructed segmented compensation coefficients are applied to the center wavelength of the compensated grating to separate the wavelength drift caused by strain from the wavelength drift caused by temperature, and output demodulated physical quantity data.

[0042] In an optional embodiment, when the trajectory curvature exceeds a preset curvature threshold, the wavelength shift pattern corresponding to the period of curvature abnormality in the historical initial wavelength shift estimation sequence is identified, and the frequency domain feature components of the wavelength shift pattern are extracted, including:

[0043] When the trajectory curvature exceeds the preset curvature threshold, a time window sliding mechanism is established in the historical initial wavelength offset estimation value sequence, and the time window length is set.

[0044] For each time window, calculate the corresponding historical trajectory curvature of the historical initial wavelength offset estimate, and mark the time window in which the historical trajectory curvature exceeds the preset curvature threshold as a curvature abnormal period.

[0045] Extract the initial wavelength offset estimate during the curvature anomaly period to form an abnormal wavelength offset subsequence;

[0046] Calculate the correlation coefficient between the abnormal wavelength offset subsequence and the current initial wavelength offset estimate, filter out abnormal wavelength offset subsequences with correlation coefficients greater than a preset correlation coefficient threshold, and determine the wavelength offset pattern;

[0047] The wavelength shift mode is subjected to frequency domain transformation to obtain amplitude spectrum and phase spectrum. The frequency component corresponding to the amplitude peak is extracted from the amplitude spectrum as the main frequency feature. The phase gradient is calculated from the phase spectrum, and the frequency component whose phase gradient change rate exceeds the preset phase gradient threshold is extracted as the phase abrupt change feature.

[0048] The dominant frequency feature is combined with the phase change feature to form the frequency domain feature component.

[0049] A second aspect of the present invention provides a wavelength deviation compensation system for rapid wavelength resolution of a far-end fiber Bragg grating, comprising:

[0050] The first unit is used to collect the reflection spectral signals of multiple fiber gratings that have traveled different transmission distances;

[0051] The second unit is used to perform frequency domain transformation on the reflected spectral signal, separate the linear phase component and the nonlinear phase component, establish a linear mapping relationship between transmission distance and dispersion accumulation, calculate the transmission dispersion offset and intrinsic wavelength offset, and obtain the initial wavelength offset estimate of the grating at different spatial positions.

[0052] The third unit is used to identify the refractive index abrupt change section and the refractive index stable section based on the initial wavelength offset estimate, and to establish a piecewise compensation coefficient. The correction amount is calculated through tensor decomposition to correct the initial wavelength offset estimate and obtain the compensated grating center wavelength.

[0053] The fourth unit is used to construct a wavelength-strain-temperature three-dimensional space, map the compensated grating center wavelength to a trajectory in the three-dimensional space, calculate the trajectory curvature, and when the trajectory curvature exceeds the preset curvature threshold, reconstruct the piecewise compensation coefficients based on the historical initial wavelength offset estimation value sequence and output physical quantity demodulation data.

[0054] The fifth unit is used to map the wavelength reconstruction residual in the demodulated physical quantity data to the three-dimensional space and update the piecewise compensation coefficient according to the residual gradient direction.

[0055] A third aspect of the present invention provides an electronic device, comprising:

[0056] processor;

[0057] Memory used to store processor-executable instructions;

[0058] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0059] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0060] In this embodiment of the invention, linear and nonlinear phase components are separated by frequency domain transformation, establishing a linear mapping relationship between transmission distance and dispersion accumulation. This achieves precise demodulation of the center wavelength of the far-end grating, effectively eliminating the influence of dispersion caused by long-distance transmission on grating measurement and improving the monitoring accuracy of the fiber optic grating sensing network. It enables adaptive compensation for gratings at different spatial locations, solving the problem of poor compensation performance of traditional methods in complex environments. This improves the measurement accuracy of the system in regions with abrupt changes in refractive index and enhances the adaptability of the fiber optic sensing network to environmental changes. Furthermore, a three-dimensional spatial monitoring mechanism for wavelength-strain-temperature and trajectory curvature is constructed. Combined with wavelength reconstruction residual mapping technology, dynamic updates of the compensation coefficients are achieved, enabling the system to have self-correction capabilities. This significantly improves the long-term stability and reliability of the remote fiber optic sensing system, providing strong support for the application of distributed fiber optic sensing technology in fields such as structural health monitoring. Attached Figure Description

[0061] Figure 1 This is a flowchart illustrating the wavelength deviation compensation method for rapid wavelength resolution of far-end fiber Bragg gratings according to an embodiment of the present invention.

[0062] Figure 2 This is a flowchart of fiber optic sensing signal compensation and separation. Detailed Implementation

[0063] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0064] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0065] Figure 1 This is a flowchart illustrating the wavelength deviation compensation method for rapid wavelength resolution of far-end fiber Bragg gratings according to an embodiment of the present invention, as shown below. Figure 1 As shown, the method includes:

[0066] The reflectance spectral signals of multiple fiber gratings that have traveled different transmission distances were collected;

[0067] The reflected spectral signal is frequency domain transformed to separate the linear phase component and the nonlinear phase component, a linear mapping relationship between transmission distance and dispersion accumulation is established, the transmission dispersion offset and intrinsic wavelength offset are calculated, and the initial wavelength offset estimate of the grating at different spatial positions is obtained.

[0068] Based on the initial wavelength offset estimate, the refractive index abrupt change segment and the refractive index stable segment are identified, and a segmented compensation coefficient is established. The correction amount is calculated through tensor decomposition to correct the initial wavelength offset estimate and obtain the compensated grating center wavelength.

[0069] A three-dimensional space of wavelength-strain-temperature is constructed, and the compensated grating center wavelength is mapped to a trajectory in the three-dimensional space. The trajectory curvature is calculated. When the trajectory curvature exceeds the preset curvature threshold, the piecewise compensation coefficient is reconstructed based on the historical initial wavelength offset estimation value sequence, and the physical quantity demodulation data is output.

[0070] The wavelength reconstruction residual in the demodulated physical quantity data is mapped to the three-dimensional space, and the piecewise compensation coefficient is updated according to the residual gradient direction.

[0071] In one optional implementation, the reflected spectral signal is subjected to frequency domain transformation to separate the linear phase component and the nonlinear phase component, a linear mapping relationship between transmission distance and dispersion accumulation is established, and the transmission dispersion offset and intrinsic wavelength offset are calculated to obtain the initial wavelength offset estimate of the grating at different spatial locations, including:

[0072] The reflection spectrum signal is transformed from the time domain to the frequency domain to obtain the frequency domain spectral distribution;

[0073] A continuous phase response curve is extracted from the frequency domain spectral distribution. The phase response curve is then fitted with a polynomial according to the fitting weights to separate the linear phase component composed of the linear fitting parameters and the nonlinear phase component composed of the residual phase distribution.

[0074] Calculate the group delay difference corresponding to different transmission distances based on the linear phase components, and establish a linear mapping relationship between transmission distance and dispersion accumulation.

[0075] The nonlinear phase component is converted back to the wavelength domain to extract the intrinsic center wavelength of each grating position.

[0076] Based on the linear mapping relationship, the group delay difference is converted into the dispersion accumulation at each spatial location, and the transmission dispersion offset and intrinsic wavelength offset of each grating are calculated in combination with the intrinsic center wavelength.

[0077] By vector synthesis of the transmission dispersion offset and the intrinsic wavelength offset, the initial wavelength offset estimate of the grating at different spatial positions is obtained.

[0078] In one specific implementation, the reflection spectrum signal of the fiber Bragg grating sensor is acquired, typically by obtaining wavelength and intensity correspondence data from the reflected light of the grating through a spectrometer or similar device. The acquired reflection spectrum signal is then transformed from the time domain to the frequency domain using a Fast Fourier Transform (FFT) algorithm. For example, performing a FFT on a reflection spectrum signal with a sampling frequency of 10 GHz and 1024 data points yields the corresponding frequency domain spectral distribution.

[0079] After obtaining the frequency domain spectral distribution, phase response information needs to be extracted. By calculating the amplitude and phase of the frequency domain signal, a set of phase values ​​distributed along the frequency axis can be obtained. Since the phase values ​​typically jump between -π and π, phase unwrapping processing is required to make the phase curve exhibit continuous variation. In experimental testing, grating reflection signals in the wavelength range of 1530nm to 1560nm were processed, and a continuously varying phase curve was obtained after unwrapping.

[0080] Polynomial fitting was performed on the continuous phase response curve after unwinding. Signal intensity at different frequencies was considered as fitting weights, with higher weights for frequencies with stronger signal intensity. The least squares method was used, employing a third-order polynomial. In the fitting results, the coefficients of the first-order terms correspond to the linear phase component, representing the group delay characteristics; the remaining higher-order terms and the fitting residuals constitute the nonlinear phase component, reflecting the wavelength shift caused by the grating's intrinsic characteristics and dispersion. In actual testing, the coefficients of the linear phase component were approximately -2.5 × 10⁻⁶. -10 s represents the basic delay characteristic of the optical signal in the transmission medium.

[0081] The group delay difference corresponding to different transmission distances is calculated based on the extracted linear phase components. The group delay difference refers to the difference in transmission speed of different wavelength components of the optical signal due to dispersion effects during transmission. By analyzing the linear phase components of gratings installed at multiple different distances, a linear mapping relationship between transmission distance and cumulative dispersion is established. For example, measurement results at four distance points (0km, 10km, 20km, and 30km) show that the coefficient of correlation between cumulative dispersion and distance is approximately 17 ps / (nm·km), consistent with the dispersion characteristics of standard single-mode fiber.

[0082] The nonlinear phase components are transformed back into the wavelength domain using an inverse Fourier transform, and the intrinsic center wavelengths at each grating position are extracted. In the wavelength domain, the nonlinear phase components manifest as distortions in the wavelength spectrum, with their peak positions corresponding to the intrinsic center wavelengths of the gratings. A peak-finding algorithm can accurately determine the center wavelength of each grating. For example, in the test, the intrinsic center wavelengths of the gratings at four different positions were 1541.23 nm, 1545.67 nm, 1550.12 nm, and 1555.48 nm, respectively.

[0083] Using the established linear mapping relationship, the group delay difference is converted into the cumulative dispersion at each spatial location. For example, for a grating with a distance of 20 km, the cumulative dispersion corresponding to its group delay difference is 340 ps / nm. Combined with the extracted intrinsic center wavelength, the transmission dispersion offset of each grating is calculated. The transmission dispersion offset represents the wavelength measurement deviation caused by dispersion effects; its calculation must consider the difference between the center wavelength and the reference wavelength, as well as the cumulative dispersion. For a grating near 1550 nm, the wavelength offset caused by a 20 km transmission distance is approximately 0.05 nm.

[0084] While determining the transmission dispersion offset, the intrinsic wavelength offset was extracted through analysis of the nonlinear phase component. The intrinsic wavelength offset represents the wavelength change caused by the grating's own characteristics or external factors such as stress and temperature. In the experiment, a 1℃ temperature change caused an intrinsic wavelength offset of approximately 0.01nm.

[0085] The transmitted dispersion offset and the intrinsic wavelength offset are vector-synthesized to obtain the initial wavelength offset estimates for gratings at different spatial locations. The vector synthesis considers the directionality of both offsets, assigning positive values ​​for increasing wavelength and negative values ​​for decreasing wavelength. Experimental verification shows that the wavelength offset estimated by this method has an error within ±0.003 nm compared to the actual measured value, meeting the requirements for high-precision distributed sensing.

[0086] In one optional implementation, extracting a continuous phase response curve from the frequency domain spectral distribution, performing polynomial fitting on the phase response curve according to fitting weights, and separating the linear phase component composed of linear fitting parameters and the nonlinear phase component composed of the residual phase distribution includes:

[0087] Phase unwrapping processing is performed on the complex amplitude at each frequency point in the frequency domain spectral distribution to generate a continuous phase response curve;

[0088] Each sampling point in the continuous phase response curve is assigned a fitting weight according to the amplitude information. The continuous phase response curve is then fitted with a polynomial according to the fitting weight to separate the linear fitting parameters that characterize the first power relationship of frequency and the nonlinear fitting parameters that characterize the higher power relationship of frequency.

[0089] The linear fitting parameters are multiplied by each frequency point to obtain the linear phase component corresponding to the transmission distance. The residual phase distribution is obtained by subtracting the linear phase component from the continuous phase response curve.

[0090] The ratio of the phase difference to the frequency interval of adjacent frequency points in the residual phase distribution is calculated to generate a phase change rate sequence. The direction of grating period change is determined by statistically analyzing the positive and negative values ​​of the phase change rate sequence. The grating period gradient rate is extracted by the adjacent difference values ​​of the phase change rate sequence. The structural features of the residual phase distribution are calibrated based on the direction of grating period change and the grating period gradient rate.

[0091] The residual phase distribution calibrated by structural features is determined as the nonlinear phase component corresponding to the grating structure.

[0092] In one specific implementation, when performing phase unwrapping processing on the complex amplitude at each frequency point in the frequency domain spectral distribution, frequency domain spectral distribution data is obtained, including frequency point information and corresponding complex amplitude information. The complex amplitude can be represented as a combination of amplitude and phase, where the initial phase value is typically limited to the range of -π to π. Since the phase value exhibits periodic jumps, phase unwrapping processing is required to obtain a continuous phase response curve. Phase unwrapping processing employs a phase accumulation method. By detecting the phase difference between adjacent frequency points, when a phase difference exceeding π is detected, it is determined to be a phase jump point. All phase values ​​after this point are adjusted accordingly by integer multiples of 2π, ensuring the overall phase curve exhibits continuous variation characteristics. For example, when the phase value at frequency point f1 is 2.9π, while the phase value at the adjacent frequency point f2 is -2.95π, the actual phase difference is approximately 0.05π rather than 5.85π. After phase unwrapping processing, the phase at point f2 is adjusted to 3.05π, ensuring the continuity of the phase curve.

[0093] When assigning fitting weights to each sampling point in a continuous phase response curve based on amplitude information, considering the differences in signal-to-noise ratio at different frequencies in the frequency domain spectral distribution, frequencies with larger signal amplitudes typically have higher phase accuracy. Therefore, a weighting strategy positively correlated with the amplitude of the frequency point can be adopted, where the weight coefficient is proportional to the square of the amplitude. For example, for a frequency point with an amplitude of A, its corresponding weight can be set to A0. 2 / A_max 2 Where A_max is the maximum amplitude value among all frequency points. In practice, an amplitude threshold can be set. When the amplitude of a frequency point is lower than the threshold, its weight is set to zero, thereby eliminating the interference of noise on the fitting results.

[0094] When performing polynomial fitting on continuous phase response curves according to fitting weights, the weighted least squares method is used for polynomial fitting, and the polynomial order can be selected from 3 to 5. During the fitting process, frequency is used as the independent variable and phase value as the dependent variable to construct a polynomial fitting model. The polynomial coefficients are solved by minimizing the weighted sum of squared residuals. In the obtained polynomial coefficients, the first-order coefficients correspond to linear phase characteristics and are directly related to the transmission distance; the higher-order coefficients characterize the nonlinear relationship between frequency and phase and are related to the grating structure characteristics. For example, for a certain test data, after fitting with a 4th-order polynomial, the coefficients are a0 = -0.12 and a1 = 9.83 × 10⁻⁶. -5 a2 = 2.15 × 10 -10 a3 = -4.37 × 10 -15 a4 = 1.92 × 10 -20 , where a1 is the linear fitting parameter.

[0095] When multiplying the linear fitting parameters by each frequency point, the linear phase component is obtained by directly multiplying the linear parameter a1 by the corresponding frequency value. For example, for a point with frequency f, its linear phase component is a1×f. Subtracting the linear phase component from the continuous phase response curve yields the residual phase distribution, which contains the nonlinear phase characteristics related to the grating structure.

[0096] When calculating the ratio of the phase difference to the frequency interval among adjacent frequency points in the residual phase distribution, for any adjacent frequency points f1 and f2 and their corresponding residual phases φ1 and φ2, the phase change rate (φ2-φ1) / (f2-f1) is calculated. All calculated phase change rates are arranged in frequency order to form a phase change rate sequence. By statistically analyzing the distribution of positive and negative values ​​in this sequence, the direction of change in the grating period is determined. When the number of positive values ​​is significantly greater than that of negative values, it indicates that the grating period is increasing; conversely, it indicates that the grating period is decreasing. Difference between adjacent points in the phase change rate sequence yields the second-order phase change rate, whose numerical distribution characteristics reflect the gradual change rate of the grating period. For example, when the mean of the second-order phase change rate is positive and shows an increasing trend, it indicates that the gradual change rate of the grating period is accelerating; when the mean of the second-order phase change rate is close to zero and the fluctuation is small, it indicates that the grating period is changing at an approximately uniform rate.

[0097] When calibrating the residual phase distribution based on the direction and rate of change of the grating period, the residual phase distribution is fitted and compared with the theoretical model to determine the grating structure type. For example, when the residual phase distribution exhibits quadratic function characteristics and the grating period changes linearly, it can be calibrated as a chirped grating structure; when the residual phase distribution exhibits high-order polynomial characteristics and the rate of change of the grating period changes nonlinearly, it can be calibrated as a non-uniform chirped grating structure. Through structural feature calibration, a clear correspondence is established between the residual phase distribution and the grating structure, obtaining the nonlinear phase component corresponding to the grating structure.

[0098] In one optional implementation, the refractive index abrupt change segment and the refractive index stable segment are identified based on the initial wavelength offset estimate, and a piecewise compensation coefficient is established. The correction amount is calculated through tensor decomposition to correct the initial wavelength offset estimate, resulting in the compensated grating center wavelength, including:

[0099] The initial wavelength offset estimate is segmented and statistically analyzed according to the transmission distance. The gradient rate of change of the initial wavelength offset estimate in each transmission distance segment is calculated. Transmission distance segments with a gradient rate of change exceeding a preset rate of change threshold are marked as refractive index abrupt change segments, and transmission distance segments with a gradient rate of change within the rate of change threshold are marked as refractive index stable segments.

[0100] Extract the boundary positions between the refractive index abrupt change region and the refractive index stable region, and establish piecewise compensation coefficients at the boundary positions;

[0101] Principal component analysis is performed on the temperature field components and strain field components to determine the tensor dimension size. A third-order tensor structure is constructed and filled. Constrained alternating least squares iterative operation is performed to separate the first eigenvector of the transmission distance dimension, the second eigenvector of the temperature field dimension, and the third eigenvector of the strain field dimension.

[0102] The first eigenvector is multiplied element-wise with the piecewise compensation coefficient to obtain the transmission dispersion correction; the second eigenvector and the third eigenvector are multiplied by tensor cross product to obtain the multiphysics coupling correction.

[0103] The compensated grating center wavelength is obtained by subtracting the transmission dispersion correction and multiphysics coupling correction from the initial wavelength offset estimate.

[0104] In one specific implementation, the acquired wavelength offset data is affected by dispersion effects and environmental disturbances during transmission over a 100-kilometer distance, leading to deviations in the wavelength offset data. After obtaining the initial wavelength offset estimate, it needs to be segmented. The initial wavelength offset estimate is segmented and statistically analyzed in 5-kilometer intervals according to the transmission distance, and the gradient rate of change of the initial wavelength offset estimate within each segment is calculated. The gradient rate of change is calculated by dividing the difference in wavelength offset values ​​between adjacent measurement points within a segment by the difference in distance between them, obtaining the gradient rate of change sequence within each segment, and then calculating its standard deviation. When the standard deviation exceeds a preset rate of change threshold of 0.01 nm / km, the segment is marked as a refractive index abrupt change segment; when the standard deviation is within 0.01 nm / km, the segment is marked as a refractive index stable segment. For example, in practical applications, the standard deviation of the gradient change rate in the 0-5 km section is measured to be 0.008 nm / km, which is marked as a stable refractive index section; the standard deviation of the gradient change rate in the 5-10 km section is 0.025 nm / km, which is marked as a sudden change in refractive index section.

[0105] The boundary locations between the refractive index abrupt change sections and the refractive index stable sections are extracted, and piecewise compensation coefficients are established at these boundary locations. For example, if a boundary exists at 5 km, a compensation coefficient is established at that location. The compensation coefficient for the refractive index stable section is set to 1.0, and the compensation coefficient for the refractive index abrupt change section is determined by the ratio of the local wavelength shift to the ideal wavelength shift. For example, the compensation coefficient for the 5-10 km section is calculated as the ratio of the average measured wavelength shift value of this section to the ideal wavelength shift value under the same temperature strain conditions, resulting in a compensation coefficient of 1.23. After establishing corresponding piecewise compensation coefficients at all boundary locations, a complete piecewise compensation coefficient table is formed.

[0106] Principal component analysis was performed on the temperature and strain field components to determine the tensor dimensions. The collected temperature data were used to form a temperature matrix, and the strain data were used to form a strain matrix. Eigenvalues ​​were calculated for each. The number of eigenvectors corresponding to the eigenvalues ​​with a cumulative contribution rate of 95% was selected as their respective dimensions. For example, the temperature field dimension was 4 and the strain field dimension was 3. Based on the determined dimensions, a third-order tensor structure was constructed. The first dimension was the transmission distance (set to 20 measurement points), the second dimension was the temperature field dimension (4), and the third dimension was the strain field dimension (3), forming a 20×4×3 third-order tensor. The initial wavelength offset estimate was filled into this tensor to form a complete data tensor.

[0107] Constrained alternating least squares iterative operations are performed to separate three eigenvectors through iterative optimization. The maximum number of iterations is set to 100, and the convergence threshold is 0.0001. In each iteration, two eigenvector dimensions are fixed, and the third eigenvector dimension is optimized. This process is repeated until the convergence condition or the maximum number of iterations is reached. Through this process, the first eigenvector (20-dimensional vector) representing the transmission distance dimension, the second eigenvector (4-dimensional vector) representing the temperature field dimension, and the third eigenvector (3-dimensional vector) representing the strain field dimension are obtained.

[0108] The transmission dispersion correction is obtained by performing element-wise multiplication of the first eigenvector with the piecewise compensation coefficient. Specifically, each element of the 20-dimensional first eigenvector is multiplied by the piecewise compensation coefficient at the corresponding transmission distance. For example, the 6th measurement point is located in the abrupt refractive index region, with a first eigenvector element value of 0.15 nm and a corresponding piecewise compensation coefficient of 1.23, resulting in a calculated transmission dispersion correction of 0.1845 nm. The multiplication of the second and third eigenvectors by a tensor outer product yields the multiphysics coupling correction. This tensor outer product multiplies the 4-dimensional second eigenvector with the 3-dimensional third eigenvector, generating a 4×3 matrix. This matrix represents the effect of the coupling between the temperature and strain fields on the wavelength shift. The summation of the matrix elements yields the total multiphysics coupling correction, for example, a calculated result of 0.0722 nm.

[0109] The compensated grating center wavelength is obtained by subtracting the transmission dispersion correction and multiphysics coupling correction from the initial wavelength offset estimate. For example, if the initial wavelength offset estimate for a certain measurement point is 1.2563 nm, the transmission dispersion correction is 0.1845 nm, and the multiphysics coupling correction is 0.0722 nm, then the compensated grating center wavelength is 1.2563 nm - 0.1845 nm - 0.0722 nm = 0.9996 nm. The same compensation calculation is performed on all measurement points to obtain a complete sequence of compensated grating center wavelength data. Experimental verification shows that the error between the compensated grating center wavelength and the actual wavelength is reduced from an average of 0.25 nm to 0.03 nm, a relative error reduction of 88%, effectively improving the measurement accuracy in long-distance applications.

[0110] like Figure 2 The diagram shown illustrates the process of fiber optic sensing signal compensation and separation.

[0111] In one optional implementation, principal component analysis is performed on the temperature field components and strain field components to determine the tensor dimension size, a third-order tensor structure is constructed and filled, and a constrained alternating least squares iterative operation is performed to separate the first eigenvector of the transmission distance dimension, the second eigenvector of the temperature field dimension, and the third eigenvector of the strain field dimension, including:

[0112] Principal component analysis is performed on the temperature field component and the strain field component respectively. The variance contribution rate of each principal component is extracted, the cumulative sum of the variance contribution rates of each principal component is calculated, the cumulative variance contribution rate is determined, and the number of principal components corresponding to the cumulative variance contribution rate reaching the preset variance threshold is determined as the temperature field dimension and the strain field dimension.

[0113] A third-order tensor structure is constructed based on the number of sampling points, the dimension of the temperature field, and the dimension of the strain field of the initial wavelength offset estimate.

[0114] The initial wavelength offset estimate is traversed through the three-dimensional coordinate indices of the third-order tensor structure. When there is no measured data for the tensor element corresponding to the three-dimensional coordinate index, the filled tensor element whose Euclidean distance to the three-dimensional coordinate index is less than a preset distance threshold is searched. The inverse of the Euclidean distance of each filled tensor element is calculated to determine the interpolation weight. The filled tensor elements are weighted and summed according to the interpolation weight, and the three-dimensional coordinate index position is filled to obtain the filled third-order tensor structure.

[0115] The filled third-order tensor structure is subjected to alternating least squares iterative operations with nonnegativity and orthogonality constraints. The iteration is terminated when the relative rate of change of the third-order tensor reconstruction error is continuously less than the preset convergence threshold. The factor matrices of each dimension are extracted to determine the first eigenvector, the second eigenvector, and the third eigenvector.

[0116] In one specific implementation, principal component analysis is performed on the temperature field component and the strain field component separately to determine the appropriate dimensional size. For the temperature field component, data samples were collected at multiple temperature points, including 20℃, 30℃, 40℃, 50℃, and 60℃. By calculating the covariance matrix and performing eigenvalue decomposition, the corresponding eigenvalues ​​were obtained as [3.75, 0.82, 0.31, 0.09, 0.03]. The variance contribution rate of each principal component was calculated as [0.75, 0.16, 0.06, 0.02, 0.01], and the cumulative variance contribution rate was [0.75, 0.91, 0.97, 0.99, 1.00]. A preset variance threshold of 0.95 was set. When the cumulative contribution rate reached 0.97, the corresponding number of principal components was 3. Therefore, the dimensional size of the temperature field was determined to be 3.

[0117] Similarly, principal component analysis was performed on the strain field components, and data samples were collected at multiple strain levels including 0 με, 500 με, 1000 με, 1500 με, 2000 με, and 2500 με. The calculated eigenvalues ​​are [4.82, 1.12, 0.44, 0.38, 0.15, 0.09], the variance contribution rate is [0.69, 0.16, 0.06, 0.05, 0.02, 0.01], and the cumulative variance contribution rate is [0.69, 0.85, 0.91, 0.96, 0.99, 1.00]. Based on the same preset variance threshold of 0.95, the number of principal components corresponding to a cumulative contribution rate of 0.96 is 4; therefore, the strain field dimension is determined to be 4.

[0118] Based on the above analysis, a third-order tensor structure was constructed. The first dimension of the tensor represents the transmission distance sampling points, with a size of 10, which is the number of sampling points of the initial wavelength offset estimate; the second dimension represents the temperature field, with a size of 3; and the third dimension represents the strain field, with a size of 4. Therefore, a 10×3×4 third-order tensor structure was constructed.

[0119] For this tensor structure, missing data elements need to be filled. For example, there is no actual measurement data at coordinates (5, 2, 3). Setting a preset distance threshold of 1.5, we search for filled tensor elements with an Euclidean distance less than 1.5 from this coordinate. We find four points: (4, 2, 3), (5, 1, 3), (5, 2, 2), and (6, 2, 3). Their Euclidean distances are 1.0, 1.0, 1.0, and 1.0, respectively, and their reciprocals are also 1.0, 1.0, 1.0, and 1.0. The normalized interpolation weights are all 0.25. Assuming the values ​​of these four filled elements are 0.82, 0.75, 0.78, and 0.85, respectively, a weighted sum based on the interpolation weights yields a filled value of 0.8 at (5, 2, 3). Using a similar method, the entire tensor is filled.

[0120] Perform alternating least squares iterations with nonnegativity and orthogonality constraints on the filled third-order tensor. Set the preset convergence threshold to 10. -4 Initialize three factor matrices, where the first matrix has a dimension of 10×1, the second matrix has a dimension of 3×1, and the third matrix has a dimension of 4×1. All elements are initialized to random non-negative values.

[0121] In the first iteration, the second and third factor matrices are fixed, and the first factor matrix is ​​updated. The calculated first factor matrix is ​​[0.32, 0.33, 0.31, 0.33, 0.34, 0.33, 0.32, 0.31, 0.30, 0.29]. The first and third factor matrices are fixed, and the second factor matrix is ​​updated, resulting in [0.51, 0.63, 0.58]. The first and second factor matrices are fixed, and the third factor matrix is ​​updated, resulting in [0.45, 0.52, 0.56, 0.47]. The tensor reconstruction error is calculated to be 0.085.

[0122] In the second iteration, the above process is repeated to obtain the updated three factor matrices. The first factor matrix is ​​[0.31, 0.32, 0.32, 0.33, 0.34, 0.33, 0.32, 0.31, 0.30, 0.28], the second factor matrix is ​​[0.52, 0.62, 0.59], and the third factor matrix is ​​[0.44, 0.53, 0.55, 0.48]. The calculated tensor reconstruction error is 0.078, and the relative rate of change is (0.085-0.078) / 0.085≈0.082.

[0123] The iteration continued until the 15th iteration, at which point the three factor matrices converged to the following values: first factor matrix [0.30, 0.31, 0.32, 0.33, 0.34, 0.33, 0.32, 0.31, 0.30, 0.28], second factor matrix [0.54, 0.60, 0.59], and third factor matrix [0.43, 0.53, 0.55, 0.49]. At this point, the tensor reconstruction error was 0.056, and the relative rate of change was 0.0009, which was less than the preset convergence threshold of 10 for three consecutive iterations. -4 Therefore, the iteration is terminated.

[0124] The final extracted first feature vector [0.30, 0.31, 0.32, 0.33, 0.34, 0.33, 0.32, 0.31, 0.30, 0.28] represents the feature of the transmission distance dimension; the second feature vector [0.54, 0.60, 0.59] represents the feature of the temperature field dimension; and the third feature vector [0.43, 0.53, 0.55, 0.49] represents the feature of the strain field dimension. These feature vectors can be used for subsequent sensor data analysis and processing.

[0125] In one optional implementation, a wavelength-strain-temperature three-dimensional space is constructed, and the compensated grating center wavelength is mapped to a trajectory in the three-dimensional space. The trajectory curvature is calculated, and when the trajectory curvature exceeds a preset curvature threshold, the piecewise compensation coefficients are reconstructed based on the historical initial wavelength offset estimation sequence. The output physical quantity demodulation data includes:

[0126] A three-dimensional space with wavelength, strain, and temperature as coordinate axes is established. The center wavelength of the compensated grating and the strain and temperature values ​​at the corresponding time are combined to form a three-dimensional coordinate point sequence, thus forming a trajectory in the three-dimensional space.

[0127] The trajectory tangent vector is obtained by performing differential operations on the sequence of three-dimensional coordinate points, and the trajectory curvature is determined by calculating the rate of change of the angle between adjacent trajectory tangent vectors.

[0128] When the trajectory curvature exceeds a preset curvature threshold, identify the wavelength shift pattern corresponding to the period of abnormal curvature in the historical initial wavelength shift estimation value sequence, and extract the frequency domain feature components of the wavelength shift pattern.

[0129] The initial wavelength offset estimate is matched with the frequency domain feature components, and compensation weights for different refractive index segments are assigned according to the matching degree to reconstruct the segmented compensation coefficients.

[0130] The reconstructed segmented compensation coefficients are applied to the center wavelength of the compensated grating to separate the wavelength drift caused by strain from the wavelength drift caused by temperature, and output demodulated physical quantity data.

[0131] In one specific implementation, based on the compensated grating center wavelength, a three-dimensional space is established with wavelength, strain, and temperature as coordinate axes for further analysis of wavelength shift characteristics. The compensated grating center wavelength, along with the corresponding strain and temperature values ​​at different times, forms a sequence of three-dimensional coordinate points, creating a trajectory in this three-dimensional space. For example, at a certain moment, if the compensated grating center wavelength is 1550.214 nm, the corresponding strain value is 120 microstrain, and the temperature value is 25.3 degrees Celsius, then the point (1550.214, 120, 25.3) can be marked in the three-dimensional space. By sampling at a frequency of 10 times per second for 10 minutes, 6000 such three-dimensional coordinate points can be obtained, forming a complete trajectory.

[0132] Differential operations are performed on the sequence of three-dimensional coordinate points to obtain the trajectory tangent vector. The differential operation is implemented using the central difference method. For the i-th point in the sequence, its tangent vector is equal to the coordinate difference between the (i+1)-th and (i-1)-th points divided by the corresponding time interval. An appropriate time window size is selected for smoothing, typically 5 to 10 sampling points. For example, for a three-dimensional coordinate point (1550.214, 120, 25.3) at a certain moment, if the coordinates at the previous moment are (1550.212, 119.5, 25.2) and the coordinates at the next moment are (1550.217, 120.8, 25.4), and the sampling interval is 0.1 seconds, then the tangent vector of this point is ((1550.217-1550.212) / 0.2, (120.8-119.5) / 0.2, (25.4-25.2) / 0.2), that is, (0.025, 6.5, 1).

[0133] Calculate the rate of change of the angle between adjacent tangent vectors to determine the trajectory curvature. The angle calculation is based on the ratio of the dot product of two tangent vectors to their respective magnitudes. The rate of change of the angle is the difference between adjacent angles divided by the corresponding time interval. To reduce the influence of noise, a low-pass filter can be applied to the calculated rate of change of the angle, with the cutoff frequency set to 1 / 10 of the sampling frequency. Assuming that within a certain time period, three consecutive tangent vectors are (0.022, 6.2, 0.9), (0.025, 6.5, 1), and (0.030, 7.1, 1.2), and the calculated adjacent angles are 0.012 radians and 0.017 radians respectively, with a sampling interval of 0.1 seconds, then the rate of change of the angle is (0.017 - 0.012) / 0.1 = 0.05 radians / second.

[0134] When the trajectory curvature exceeds a preset curvature threshold, it indicates an abnormal change in the relationship between the grating center wavelength and strain and temperature, possibly due to external interference or changes in the fiber condition. The preset curvature threshold is set according to the actual application environment, typically between 0.03 and 0.1 radians / second. When the curvature exceeds the threshold, it is necessary to identify the wavelength shift pattern corresponding to the abnormal curvature period in the historical initial wavelength shift estimation value sequence. For example, if in a 10-minute sampling sequence, the trajectory curvature is detected to continuously exceed the threshold of 0.05 radians / second from the 245th second to the 268th second, reaching 0.07 radians / second, then the wavelength shift estimation value sequence within that period is extracted as the abnormal pattern.

[0135] The frequency domain characteristic components of the wavelength shift mode are extracted, and the time-domain wavelength shift sequence is transformed to the frequency domain using a Fast Fourier Transform (FFT). To obtain more accurate spectral characteristics, a Hanning window function can be applied to reduce spectral leakage. For the aforementioned 23-second anomalous wavelength shift mode, assuming a sampling rate of 10 Hz, there are 230 sampling points. After the Fourier transform, significant amplitude peaks may be observed at 2 Hz and 5 Hz, at 0.015 nm and 0.008 nm, respectively. These frequencies and amplitudes constitute the frequency domain characteristic components of this wavelength shift mode.

[0136] The current initial wavelength shift estimate is matched with the identified frequency domain feature components to calculate the matching degree. The matching degree calculation is based on the frequency domain correlation between the current wavelength shift sequence and historical anomalous patterns. Wavelength shift data from the most recent 10 seconds is taken, and its frequency domain representation is obtained through Fourier transform. The amplitudes at 2Hz and 5Hz are calculated, assumed to be 0.014 nm and 0.007 nm, respectively. Compared with the frequency domain feature components of historical anomalous patterns, the matching degree can be expressed as the normalized value of the amplitude difference at the frequency points, such as (1 - |0.015 - 0.014| / 0.015) × (1 - |0.008 - 0.007| / 0.008) = 0.94.

[0137] The compensation weights for different refractive index segments are assigned based on the matching degree, and the segmented compensation coefficients are reconstructed. A higher matching degree indicates that the interference pattern affecting the current wavelength offset estimate is more similar to historical anomaly patterns, requiring stronger compensation adjustments. Assuming the original segmented compensation coefficients are 0.92, 1.08, and 0.95 at 2500m, 5500m, and 8000m respectively, and the matching degree is 0.94, the compensation coefficients can be adjusted based on the matching degree. The new compensation coefficients can be expressed as the product of the original compensation coefficient and (1 + matching degree × adjustment factor). If the adjustment factor is set to 0.05, the reconstructed compensation coefficients are 0.92 × (1 + 0.94 × 0.05) = 0.963, 1.08 × (1 + 0.94 × 0.05) = 1.133, and 0.95 × (1 + 0.94 × 0.05) = 0.995 respectively.

[0138] The reconstructed segmented compensation coefficients are applied to the center wavelength of the compensated grating to further improve the accuracy of wavelength offset compensation. Taking a point with a transmission distance of 5000 meters as an example, if the center wavelength of the compensated grating is 1550.214 nm, and this point is located in the second refractive index segment, the corresponding reconstruction compensation coefficient is 1.133, then the final compensation result is 1550.214 × 1.133 = 1756.393 nm.

[0139] Based on the precisely compensated grating center wavelength, wavelength drift caused by strain is separated from wavelength drift caused by temperature. This separation process utilizes the different characteristics of the effects of strain and temperature on wavelength. Generally, wavelength drift caused by temperature changes has lower frequency characteristics, while wavelength drift caused by strain changes has relatively higher frequency characteristics. Using a bandpass filter bank, wavelength drift in different frequency bands can be separated. The low-frequency portion (0-0.5Hz) mainly corresponds to temperature-induced drift, while the high-frequency portion (0.5-5Hz) mainly corresponds to strain-induced drift. For example, for the final compensation result of 1756.393 nm, after frequency separation, the wavelength drift caused by temperature may be 1.245 nm, and the wavelength drift caused by strain may be 0.523 nm.

[0140] Based on the wavelength sensitivity parameters of the fiber Bragg grating, the separated wavelength drift is converted into actual physical quantities. Assuming a temperature sensitivity of 10 picometers / degree Celsius and a strain sensitivity of 1.2 picometers / microstrain, a wavelength drift caused by a temperature of 1.245 nanometers corresponds to a temperature change of 124.5 degrees Celsius, and a wavelength drift caused by a strain of 0.523 nanometers corresponds to a strain change of 435.8 microstrains. These converted physical quantities are the final demodulated data, which can be used for subsequent monitoring and analysis. Through this wavelength deviation compensation and physical quantity separation method, the accuracy of wavelength resolution of the far-end fiber Bragg grating can be improved to within 0.005 nanometers, with corresponding temperature and strain resolution accuracies reaching 0.5 degrees Celsius and 5 microstrains, respectively.

[0141] In one optional implementation, when the trajectory curvature exceeds a preset curvature threshold, the wavelength shift pattern corresponding to the period of curvature abnormality in the historical initial wavelength shift estimation sequence is identified, and the frequency domain feature components of the wavelength shift pattern are extracted, including:

[0142] When the trajectory curvature exceeds the preset curvature threshold, a time window sliding mechanism is established in the historical initial wavelength offset estimation value sequence, and the time window length is set.

[0143] For each time window, calculate the corresponding historical trajectory curvature of the historical initial wavelength offset estimate, and mark the time window in which the historical trajectory curvature exceeds the preset curvature threshold as a curvature abnormal period.

[0144] Extract the initial wavelength offset estimate during the curvature anomaly period to form an abnormal wavelength offset subsequence;

[0145] Calculate the correlation coefficient between the abnormal wavelength offset subsequence and the current initial wavelength offset estimate, filter out abnormal wavelength offset subsequences with correlation coefficients greater than a preset correlation coefficient threshold, and determine the wavelength offset pattern;

[0146] The wavelength shift mode is subjected to frequency domain transformation to obtain amplitude spectrum and phase spectrum. The frequency component corresponding to the amplitude peak is extracted from the amplitude spectrum as the main frequency feature. The phase gradient is calculated from the phase spectrum, and the frequency component whose phase gradient change rate exceeds the preset phase gradient threshold is extracted as the phase abrupt change feature.

[0147] The dominant frequency feature is combined with the phase change feature to form the frequency domain feature component.

[0148] In one specific implementation, a time window sliding mechanism is established to locate the time period of curvature anomaly in the historical initial wavelength offset estimation value sequence, then the wavelength offset pattern is determined by correlation analysis, and finally the feature components are extracted by frequency domain transformation.

[0149] The method of this invention triggers the process of identifying wavelength shift patterns by real-time monitoring of trajectory curvature. When the trajectory curvature exceeds a preset curvature threshold (e.g., 0.85 rad / s), the method proceeds. 2 When the time window is set to 60 seconds, a sliding time window mechanism is initiated on the historical initial wavelength offset estimation sequence. The historical initial wavelength offset estimation sequence contains wavelength offset estimation data acquired at a sampling frequency of 5 Hz over the past 30 minutes. The system sets the time window length to 60 seconds and the window sliding step size to 5 seconds.

[0150] For each time window, the trajectory curvature corresponding to the historical initial wavelength offset estimate within the window is calculated. Specifically, for the wavelength offset estimates λt and λt+1 of every two adjacent times t and t+1 within the window, their difference Δλ = λt+1 - λt is calculated, and then the trajectory curvature is calculated based on the rate of change of the adjacent differences. If the maximum trajectory curvature within a certain time window exceeds a preset curvature threshold of 0.85 rad / s... 2 If so, then the time window is marked as a period of curvature anomaly.

[0151] For example, when the current trajectory curvature is detected to be 0.92 rad / s 2 The curvature exceeded the preset threshold. Through a time window sliding mechanism, the system identified three periods of abnormal curvature in historical data, which were located in the time intervals of [5 minutes, 6 minutes], [12 minutes, 13 minutes] and [22 minutes, 23 minutes] within the past 30 minutes.

[0152] Initial wavelength shift estimates were extracted from these periods of curvature anomalies, forming anomalous wavelength shift subsequences A, B, and C, respectively. Subsequence A contains 300 sampling points (60 seconds × 5 Hz), with wavelength shift values ​​ranging from [1532.25 nm to 1532.46 nm]; subsequence B also contains 300 sampling points, with wavelength shift values ​​ranging from [1532.31 nm to 1532.52 nm]; and subsequence C also contains 300 sampling points, with wavelength shift values ​​ranging from [1532.28 nm to 1532.49 nm].

[0153] The correlation coefficients of these anomalous wavelength shift subsequences with the current initial wavelength shift estimate sequence (300 sampling points within the last 60 seconds, with wavelength shift values ​​ranging from [1532.29 nm to 1532.50 nm]) are calculated. The correlation coefficients of subsequence A and the current sequence are 0.72, 0.91, and 0.88, respectively. If a preset correlation coefficient threshold of 0.80 is set, subsequences B and C are selected as wavelength shift patterns.

[0154] For a given wavelength shift pattern, a frequency domain transformation is performed to obtain the amplitude and phase spectra. Preferably, a Fast Fourier Transform (FFT) algorithm is used to convert 300 time-domain sampling points into a frequency-domain representation. Performing a frequency domain transformation on subsequence B yields an amplitude spectrum with significant peaks at 0.25 Hz, 0.75 Hz, and 1.2 Hz, with peak amplitudes of 0.056, 0.038, and 0.027, respectively; the phase spectrum shows phase values ​​of 0.45π, 1.2π, and 1.7π at these frequency points. Performing the same operation on subsequence C yields an amplitude spectrum with significant peaks at 0.28 Hz, 0.72 Hz, and 1.25 Hz, with peak amplitudes of 0.052, 0.041, and 0.023, respectively; the phase spectrum shows phase values ​​of 0.42π, 1.25π, and 1.75π at these frequency points.

[0155] The frequency components corresponding to the amplitude peaks are extracted from the amplitude spectrum as the dominant frequency features. In subsequence B, the dominant frequency features are [0.25Hz, 0.75Hz, 1.2Hz]; in subsequence C, the dominant frequency features are [0.28Hz, 0.72Hz, 1.25Hz].

[0156] The phase gradient of the phase spectrum, i.e., the rate of phase change between adjacent frequency points, is calculated. In subsequence B, the phase gradient in the frequency range [0.70Hz, 0.80Hz] is 2.5π / Hz, exceeding the preset phase gradient threshold of 2.0π / Hz; in subsequence C, the phase gradient in the frequency range [0.68Hz, 0.78Hz] is 2.6π / Hz, also exceeding the preset threshold. The frequency components with phase gradient change rates exceeding the threshold are extracted as phase abrupt change features, namely [0.75Hz] and [0.72Hz].

[0157] The dominant frequency characteristics and phase abrupt change characteristics are combined to form frequency domain feature components. For subsequence B, the frequency domain feature components are {dominant frequency: [0.25Hz, 0.75Hz, 1.2Hz], phase abrupt change: [0.75Hz]}; for subsequence C, the frequency domain feature components are {dominant frequency: [0.28Hz, 0.72Hz, 1.25Hz], phase abrupt change: [0.72Hz]}. These frequency domain feature components can accurately identify abnormal wavelength shift patterns, providing an important basis for subsequent wavelength correction.

[0158] In this embodiment, these extracted frequency domain feature components are used to adjust the wavelength compensation model, enabling the system to perform targeted corrections for specific types of wavelength offset modes, which significantly improves the measurement accuracy and stability of the fiber optic sensing system in complex environments.

[0159] The wavelength deviation compensation system for rapid wavelength resolution of far-end fiber Bragg gratings according to embodiments of the present invention includes:

[0160] The first unit is used to collect the reflection spectral signals of multiple fiber gratings that have traveled different transmission distances;

[0161] The second unit is used to perform frequency domain transformation on the reflected spectral signal, separate the linear phase component and the nonlinear phase component, establish a linear mapping relationship between transmission distance and dispersion accumulation, calculate the transmission dispersion offset and intrinsic wavelength offset, and obtain the initial wavelength offset estimate of the grating at different spatial positions.

[0162] The third unit is used to identify the refractive index abrupt change section and the refractive index stable section based on the initial wavelength offset estimate, and to establish a piecewise compensation coefficient. The correction amount is calculated through tensor decomposition to correct the initial wavelength offset estimate and obtain the compensated grating center wavelength.

[0163] The fourth unit is used to construct a wavelength-strain-temperature three-dimensional space, map the compensated grating center wavelength to a trajectory in the three-dimensional space, calculate the trajectory curvature, and when the trajectory curvature exceeds the preset curvature threshold, reconstruct the piecewise compensation coefficients based on the historical initial wavelength offset estimation value sequence and output physical quantity demodulation data.

[0164] The fifth unit is used to map the wavelength reconstruction residual in the demodulated physical quantity data to the three-dimensional space and update the piecewise compensation coefficient according to the residual gradient direction.

[0165] A third aspect of the present invention provides an electronic device, comprising:

[0166] processor;

[0167] Memory used to store processor-executable instructions;

[0168] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0169] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0170] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0171] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions 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 wavelength deviation compensation method for rapid wavelength resolution of far-end fiber optic gratings, characterized in that, include: The reflectance spectral signals of multiple fiber gratings that have traveled different transmission distances were collected; The reflected spectral signal is frequency domain transformed to separate the linear phase component and the nonlinear phase component, a linear mapping relationship between transmission distance and dispersion accumulation is established, the transmission dispersion offset and intrinsic wavelength offset are calculated, and the initial wavelength offset estimate of the grating at different spatial positions is obtained. Based on the initial wavelength offset estimate, the refractive index abrupt change segment and the refractive index stable segment are identified, and a segmented compensation coefficient is established. The correction amount is calculated through tensor decomposition to correct the initial wavelength offset estimate and obtain the compensated grating center wavelength. A three-dimensional space of wavelength-strain-temperature is constructed, and the compensated grating center wavelength is mapped to a trajectory in the three-dimensional space. The trajectory curvature is calculated. When the trajectory curvature exceeds the preset curvature threshold, the piecewise compensation coefficient is reconstructed based on the historical initial wavelength offset estimation value sequence, and the physical quantity demodulation data is output. The wavelength reconstruction residual in the demodulated physical quantity data is mapped to the three-dimensional space, and the piecewise compensation coefficient is updated according to the residual gradient direction.

2. The method according to claim 1, characterized in that, The reflected spectral signal is subjected to frequency domain transformation to separate the linear phase component and the nonlinear phase component, a linear mapping relationship between transmission distance and dispersion accumulation is established, and the transmission dispersion offset and intrinsic wavelength offset are calculated to obtain the initial wavelength offset estimates of the grating at different spatial locations, including: The reflection spectrum signal is transformed from the time domain to the frequency domain to obtain the frequency domain spectral distribution; A continuous phase response curve is extracted from the frequency domain spectral distribution. The phase response curve is then fitted with a polynomial according to the fitting weights to separate the linear phase component composed of the linear fitting parameters and the nonlinear phase component composed of the residual phase distribution. Calculate the group delay difference corresponding to different transmission distances based on the linear phase components, and establish a linear mapping relationship between transmission distance and dispersion accumulation. The nonlinear phase component is converted back to the wavelength domain to extract the intrinsic center wavelength of each grating position. Based on the linear mapping relationship, the group delay difference is converted into the dispersion accumulation at each spatial location, and the transmission dispersion offset and intrinsic wavelength offset of each grating are calculated in combination with the intrinsic center wavelength. By vector synthesis of the transmission dispersion offset and the intrinsic wavelength offset, the initial wavelength offset estimate of the grating at different spatial positions is obtained.

3. The method according to claim 2, characterized in that, A continuous phase response curve is extracted from the frequency domain spectral distribution. A polynomial fit is performed on the phase response curve according to the fitting weights to separate the linear phase component composed of the linear fitting parameters and the nonlinear phase component composed of the residual phase distribution, including: Phase unwrapping processing is performed on the complex amplitude at each frequency point in the frequency domain spectral distribution to generate a continuous phase response curve; Each sampling point in the continuous phase response curve is assigned a fitting weight according to the amplitude information. The continuous phase response curve is then fitted with a polynomial according to the fitting weight to separate the linear fitting parameters that characterize the first power relationship of frequency and the nonlinear fitting parameters that characterize the higher power relationship of frequency. The linear fitting parameters are multiplied by each frequency point to obtain the linear phase component corresponding to the transmission distance. The residual phase distribution is obtained by subtracting the linear phase component from the continuous phase response curve. The ratio of the phase difference to the frequency interval of adjacent frequency points in the residual phase distribution is calculated to generate a phase change rate sequence. The direction of grating period change is determined by statistically analyzing the positive and negative values ​​of the phase change rate sequence. The grating period gradient rate is extracted by the adjacent difference values ​​of the phase change rate sequence. The structural features of the residual phase distribution are calibrated based on the direction of grating period change and the grating period gradient rate. The residual phase distribution calibrated by structural features is determined as the nonlinear phase component corresponding to the grating structure.

4. The method according to claim 1, characterized in that, Based on the initial wavelength shift estimate, abrupt refractive index change segments and stable refractive index segments are identified, and piecewise compensation coefficients are established. The correction amount is calculated through tensor decomposition to correct the initial wavelength shift estimate, resulting in the compensated grating center wavelength, including: The initial wavelength offset estimate is segmented and statistically analyzed according to the transmission distance. The gradient rate of change of the initial wavelength offset estimate in each transmission distance segment is calculated. Transmission distance segments with a gradient rate of change exceeding a preset rate of change threshold are marked as refractive index abrupt change segments, and transmission distance segments with a gradient rate of change within the rate of change threshold are marked as refractive index stable segments. Extract the boundary positions between the refractive index abrupt change region and the refractive index stable region, and establish piecewise compensation coefficients at the boundary positions; Principal component analysis is performed on the temperature field components and strain field components to determine the tensor dimension size. A third-order tensor structure is constructed and filled. Constrained alternating least squares iterative operation is performed to separate the first eigenvector of the transmission distance dimension, the second eigenvector of the temperature field dimension, and the third eigenvector of the strain field dimension. The first eigenvector is multiplied element-wise with the piecewise compensation coefficient to obtain the transmission dispersion correction; the second eigenvector and the third eigenvector are multiplied by tensor cross product to obtain the multiphysics coupling correction. The compensated grating center wavelength is obtained by subtracting the transmission dispersion correction and multiphysics coupling correction from the initial wavelength offset estimate.

5. The method according to claim 4, characterized in that, Principal component analysis is performed on the temperature and strain field components to determine the tensor dimensions. A third-order tensor structure is constructed and filled. Constrained alternating least squares iterative operations are performed to separate the first eigenvector of the transmission distance dimension, the second eigenvector of the temperature field dimension, and the third eigenvector of the strain field dimension, including: Principal component analysis is performed on the temperature field component and the strain field component respectively. The variance contribution rate of each principal component is extracted, the cumulative sum of the variance contribution rates of each principal component is calculated, the cumulative variance contribution rate is determined, and the number of principal components corresponding to the cumulative variance contribution rate reaching the preset variance threshold is determined as the temperature field dimension and the strain field dimension. A third-order tensor structure is constructed based on the number of sampling points, the dimension of the temperature field, and the dimension of the strain field of the initial wavelength offset estimate. The initial wavelength offset estimate is traversed through the three-dimensional coordinate indices of the third-order tensor structure. When there is no measured data for the tensor element corresponding to the three-dimensional coordinate index, the filled tensor element whose Euclidean distance to the three-dimensional coordinate index is less than a preset distance threshold is searched. The inverse of the Euclidean distance of each filled tensor element is calculated to determine the interpolation weight. The filled tensor elements are weighted and summed according to the interpolation weight, and the three-dimensional coordinate index position is filled to obtain the filled third-order tensor structure. The filled third-order tensor structure is subjected to alternating least squares iterative operations with nonnegativity and orthogonality constraints. The iteration is terminated when the relative rate of change of the third-order tensor reconstruction error is continuously less than the preset convergence threshold. The factor matrices of each dimension are extracted to determine the first eigenvector, the second eigenvector, and the third eigenvector.

6. The method according to claim 1, characterized in that, A three-dimensional wavelength-strain-temperature space is constructed, and the compensated grating center wavelength is mapped to a trajectory in the three-dimensional space. The trajectory curvature is calculated. When the trajectory curvature exceeds a preset curvature threshold, the piecewise compensation coefficients are reconstructed based on the historical initial wavelength offset estimation sequence. The output demodulated physical quantity data includes: A three-dimensional space with wavelength, strain, and temperature as coordinate axes is established. The center wavelength of the compensated grating and the strain and temperature values ​​at the corresponding time are combined to form a three-dimensional coordinate point sequence, thus forming a trajectory in the three-dimensional space. The trajectory tangent vector is obtained by performing differential operations on the sequence of three-dimensional coordinate points, and the trajectory curvature is determined by calculating the rate of change of the angle between adjacent trajectory tangent vectors. When the trajectory curvature exceeds a preset curvature threshold, identify the wavelength shift pattern corresponding to the period of abnormal curvature in the historical initial wavelength shift estimation value sequence, and extract the frequency domain feature components of the wavelength shift pattern. The initial wavelength offset estimate is matched with the frequency domain feature components, and compensation weights for different refractive index segments are assigned according to the matching degree to reconstruct the segmented compensation coefficients. The reconstructed segmented compensation coefficients are applied to the center wavelength of the compensated grating to separate the wavelength drift caused by strain from the wavelength drift caused by temperature, and output demodulated physical quantity data.

7. The method according to claim 6, characterized in that, When the trajectory curvature exceeds a preset curvature threshold, the wavelength shift pattern corresponding to the period of abnormal curvature in the historical initial wavelength shift estimation sequence is identified, and the frequency domain feature components of the wavelength shift pattern are extracted, including: When the trajectory curvature exceeds the preset curvature threshold, a time window sliding mechanism is established in the historical initial wavelength offset estimation value sequence, and the time window length is set. For each time window, calculate the corresponding historical trajectory curvature of the historical initial wavelength offset estimate, and mark the time window in which the historical trajectory curvature exceeds the preset curvature threshold as a curvature abnormal period. Extract the initial wavelength offset estimate during the curvature anomaly period to form an abnormal wavelength offset subsequence; Calculate the correlation coefficient between the abnormal wavelength offset subsequence and the current initial wavelength offset estimate, filter out abnormal wavelength offset subsequences with correlation coefficients greater than a preset correlation coefficient threshold, and determine the wavelength offset pattern; The wavelength shift mode is subjected to frequency domain transformation to obtain amplitude spectrum and phase spectrum. The frequency component corresponding to the amplitude peak is extracted from the amplitude spectrum as the main frequency feature. The phase gradient is calculated from the phase spectrum, and the frequency component whose phase gradient change rate exceeds the preset phase gradient threshold is extracted as the phase abrupt change feature. The dominant frequency feature is combined with the phase change feature to form the frequency domain feature component.

8. A wavelength deviation compensation system for rapid wavelength resolution of a far-end fiber Bragg grating, used to implement the method of any one of claims 1-7, characterized in that, include: The first unit is used to collect the reflection spectral signals of multiple fiber gratings that have traveled different transmission distances; The second unit is used to perform frequency domain transformation on the reflected spectral signal, separate the linear phase component and the nonlinear phase component, establish a linear mapping relationship between transmission distance and dispersion accumulation, calculate the transmission dispersion offset and intrinsic wavelength offset, and obtain the initial wavelength offset estimate of the grating at different spatial positions. The third unit is used to identify the refractive index abrupt change section and the refractive index stable section based on the initial wavelength offset estimate, and to establish a piecewise compensation coefficient. The correction amount is calculated through tensor decomposition to correct the initial wavelength offset estimate and obtain the compensated grating center wavelength. The fourth unit is used to construct a wavelength-strain-temperature three-dimensional space, map the compensated grating center wavelength to a trajectory in the three-dimensional space, calculate the trajectory curvature, and when the trajectory curvature exceeds the preset curvature threshold, reconstruct the piecewise compensation coefficients based on the historical initial wavelength offset estimation value sequence and output physical quantity demodulation data. The fifth unit is used to map the wavelength reconstruction residual in the demodulated physical quantity data to the three-dimensional space and update the piecewise compensation coefficient according to the residual gradient direction.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images