Wavelength deviation compensation method and system for rapid wavelength analysis of far-end fiber bragg grating

Abstract

The invention provides a wavelength deviation compensation method and system for rapid wavelength analysis of a far-end fiber bragg grating, and relates to the technical field of fiber bragg grating sensing, and the method comprises the steps: collecting grating reflection spectrum signals of different transmission distances, carrying out frequency domain transformation to separate phase components, and building a mapping relation between the transmission distance and dispersion cumulant; calculating an initial wavelength offset estimation value; identifying a refractive index section, and establishing a segmentation compensation coefficient; constructing a wavelength-strain-temperature three-dimensional space, and optimizing a compensation coefficient according to the trajectory curvature; and updating the compensation coefficient according to the residual gradient. The wavelength analysis precision of the far-end fiber bragg grating can be improved, and the influence of dispersion on measurement is reduced.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of fiber grating sensing, and in particular to a wavelength deviation compensation method and system for remote fiber grating wavelength rapid analysis. BACKGROUND

[0002] Fiber grating sensing technology is a technology for detecting physical parameters by using the characteristic that the center wavelength of a fiber grating changes with external physical quantities (such as strain, temperature, etc.). A remote fiber grating sensing system allows the sensing element to be placed in an area away from the demodulation device, and transmits the reflected spectrum signal through a long-distance optical fiber to realize real-time monitoring of the remote environmental parameters.

[0003] During long-distance transmission of the optical fiber, due to the existence of the fiber dispersion effect, the optical signal will have a wavelength-dependent delay during transmission, resulting in distortion of the received grating reflection spectrum. This distortion causes the traditional wavelength analysis algorithm (such as the centroid method, Gaussian fitting method, etc.) to have a large error in the remote grating wavelength analysis, affecting the measurement accuracy. Especially in a multi-grating cascade system, gratings at different positions are affected by different degrees of dispersion, and their wavelength analysis errors are also different.

[0004] The current remote fiber grating wavelength analysis technology still has some defects and deficiencies. The traditional wavelength analysis method cannot effectively distinguish between the intrinsic wavelength shift of the grating and the wavelength shift caused by transmission dispersion. In a long-distance transmission scenario, the cumulative effect of dispersion causes the analysis accuracy to decrease significantly with the increase of the transmission distance. The existing dispersion compensation technology is mainly based on a uniform dispersion model and cannot adapt to the existence of refractive index sudden change sections in the optical fiber, such as local dispersion changes caused by non-uniform factors such as optical fiber joints and bending. This makes the compensation effect unsatisfactory. The existing technology lacks an adaptive wavelength deviation compensation mechanism. In the case of a drastic change in environmental conditions, it is difficult for fixed compensation parameters to track real-time dispersion changes, especially the fluctuation of the fiber dispersion parameter caused by temperature changes, which makes the measurement accuracy unstable in the long-term monitoring process. SUMMARY

[0005] The embodiment of the present application provides a wavelength deviation compensation method and system for remote fiber grating wavelength rapid analysis, which can solve the problems in the prior art.

[0006] In a first aspect of the embodiment of the present application, a wavelength deviation compensation method for remote fiber grating wavelength rapid analysis is provided, comprising: Collecting the reflection spectrum signals of a plurality of fiber gratings that have experienced different transmission distances; performing frequency domain transformation on the reflection spectrum signal, separating a linear phase component and a nonlinear phase component, establishing a linear mapping relationship between a transmission distance and a dispersion accumulation, calculating a transmission dispersion offset and an intrinsic wavelength offset, and obtaining an initial wavelength offset estimation value of gratings at different spatial positions; According to the initial wavelength offset estimation value, identify the refractive index mutation section and the refractive index stable section, and establish a segmented compensation coefficient. By tensor decomposition, calculate the correction, correct the initial wavelength offset estimation value, and obtain the compensated grating center wavelength. 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 a preset curvature threshold, reconstruct the segmented compensation coefficient according to the historical initial wavelength offset estimation value sequence, and output the physical quantity demodulation data. Map the wavelength reconstruction residual in the physical quantity demodulation data to the three-dimensional space, and update the segmented compensation coefficient according to the residual gradient direction.

[0007] In an alternative embodiment, performing frequency domain transformation on the reflection spectrum signal, separating a linear phase component and a nonlinear phase component, establishing a linear mapping relationship between a transmission distance and a dispersion accumulation, calculating a transmission dispersion offset and an intrinsic wavelength offset, and obtaining an initial wavelength offset estimation value of gratings at different spatial positions includes: Transform the reflection spectrum signal from time domain to frequency domain to obtain a frequency domain spectrum distribution; Extract a continuous phase response curve in the frequency domain spectrum distribution, perform polynomial fitting on the phase response curve according to a fitting weight, separate a linear phase component composed of linear fitting parameters and a nonlinear phase component composed of residual phase distribution; According to the linear phase component, calculate the group delay difference corresponding to different transmission distances, and establish a linear mapping relationship between the transmission distance and the dispersion accumulation; Convert the nonlinear phase component back to the wavelength domain, and extract the intrinsic center wavelength of each grating position; According to the linear mapping relationship, convert the group delay difference to the dispersion accumulation of each spatial position, and combine the intrinsic center wavelength to calculate the transmission dispersion offset and the intrinsic wavelength offset of each grating; Vector composition of transmission dispersion offset and intrinsic wavelength offset, to obtain the initial wavelength offset estimation value of gratings at different spatial positions.

[0008] In an alternative embodiment, extracting a continuous phase response curve in the frequency domain spectrum distribution, performing polynomial fitting on the phase response curve according to a fitting weight, and separating a linear phase component composed of linear fitting parameters and a nonlinear phase component composed of residual phase distribution includes: Phase unwrapping processing is performed on the complex amplitudes of the frequency points in the frequency domain spectrum distribution to generate a continuous phase response curve; Fitting weights are assigned to the sampling points in the continuous phase response curve according to the amplitude information, and polynomial fitting is performed on the continuous phase response curve according to the fitting weights to separate linear fitting parameters representing a frequency first power relationship and nonlinear fitting parameters representing a frequency high power relationship; The linear fitting parameters are multiplied by the frequency points to obtain a linear phase component corresponding to the transmission distance, and the linear phase component is subtracted from the continuous phase response curve to obtain a residual phase distribution; The ratio of the phase difference value of adjacent frequency points in the residual phase distribution to the frequency interval is calculated to generate a phase change rate sequence, the positive and negative value distribution of the phase change rate sequence is statistically determined to determine the grating period change direction, the adjacent difference value of the phase change rate sequence is extracted to obtain the grating period gradual change rate, and the residual phase distribution is labeled with structural characteristics according to the grating period change direction and the grating period gradual change rate. The residual phase distribution labeled with structural characteristics is determined as a nonlinear phase component corresponding to the grating structure.

[0009] In an optional embodiment, the refractive index mutation section and the refractive index stable section are identified according to the initial wavelength shift estimation value, and a segmented compensation coefficient is established, the initial wavelength shift estimation value is corrected by calculating a correction amount through tensor decomposition to obtain a compensated grating center wavelength. The initial wavelength shift estimation value is segmented and counted according to the transmission distance, the gradient change rate of the initial wavelength shift estimation value in each transmission distance section is calculated, the transmission distance section with a gradient change rate exceeding a preset change rate threshold is marked as a refractive index mutation section, and the transmission distance section with a gradient change rate within the change rate threshold is marked as a refractive index stable section. The boundary positions of the refractive index mutation section and the refractive index stable section are extracted, and a segmented compensation coefficient is established at the boundary positions. Principal component analysis is performed on the temperature field component and the strain field component to determine the tensor dimension size, a third-order tensor structure is constructed and filled, a constrained alternating least squares iterative operation is performed, and 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 are separated. Element-wise multiplication is performed on the first eigenvector and the segmented compensation coefficient to obtain a transmission dispersion correction amount, and tensor outer product operation is performed on the second eigenvector and the third eigenvector to obtain a multi-physical field coupling correction amount. The transmission dispersion correction amount and the multi-physical field coupling correction amount are subtracted from the initial wavelength shift estimation value to obtain a compensated grating center wavelength.

[0010] In an alternative embodiment, principal component analysis is performed on the temperature field component and the strain field component to determine the tensor dimension size, a third-order tensor structure is constructed and filled, a constrained alternating least squares iterative operation is performed, and 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 are separated. The principal component analysis is performed on the temperature field component and the strain field component, respectively, the corresponding principal component variance contribution rates are extracted, the cumulative sum of the principal component variance contribution rates is calculated, the cumulative variance contribution rate is determined, and when the cumulative variance contribution rate reaches a preset variance threshold, the number of principal components corresponding to the cumulative variance contribution rate is determined as the temperature field dimension size and the strain field dimension size. A third-order tensor structure is constructed according to the number of sampling points of the initial wavelength shift estimation value, the temperature field dimension size, and the strain field dimension size. The three-dimensional coordinate index of the initial wavelength shift estimation value in the third-order tensor structure is traversed, when the tensor element corresponding to the three-dimensional coordinate index does not exist, the filled tensor element with a Euclidean distance less than a preset distance threshold from the three-dimensional coordinate index is searched, the Euclidean distance reciprocal 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 a filled third-order tensor structure. The filled third-order tensor structure is subjected to an alternating least squares iterative operation with non-negativity constraint and orthogonality constraint, and when the relative change rate of the third-order tensor reconstruction error is less than a preset convergence threshold continuously, the iteration is terminated, and the factor matrix of each dimension is extracted to determine the first eigenvector, the second eigenvector, and the third eigenvector.

[0011] In an alternative embodiment, a wavelength-strain-temperature three-dimensional space is constructed, the compensated grating center wavelength is mapped as a trajectory in the three-dimensional space, the trajectory curvature is calculated, when the trajectory curvature exceeds a preset curvature threshold, the segmented compensation coefficient is reconstructed according to the historical initial wavelength shift estimation value sequence, and the physical quantity demodulation data is outputted, including: A three-dimensional space with wavelength, strain, and temperature as coordinate axes is established, the compensated grating center wavelength and the strain value and temperature value at the corresponding time are combined to form a three-dimensional coordinate point sequence, and a trajectory in the three-dimensional space is formed. The trajectory tangent vector is obtained by differential operation on the three-dimensional coordinate point sequence, the angle change rate of adjacent trajectory tangent vectors is calculated to determine the trajectory curvature. When the trajectory curvature exceeds the preset curvature threshold, the wavelength shift mode corresponding to the curvature abnormal period in the historical initial wavelength shift estimation value sequence is identified, and the frequency domain feature component of the wavelength shift mode is extracted. Matching the current initial wavelength offset estimation value with the frequency domain feature component, assigning compensation weights of different refractive index sections according to the matching degree, and reconstructing the segmented compensation coefficient; Applying the reconstructed segmented compensation coefficient to the compensated grating center wavelength, separating the wavelength shift caused by strain from the wavelength shift caused by temperature, and outputting physical quantity demodulation data.

[0012] In an optional embodiment, when the trajectory curvature exceeds the preset curvature threshold, a wavelength offset mode corresponding to a curvature abnormal period in the sequence of historical initial wavelength offset estimation values is identified, and extracting the frequency domain feature component of the wavelength offset mode comprises: When the trajectory curvature exceeds the preset curvature threshold, a time window sliding mechanism is established in the sequence of historical initial wavelength offset estimation values, and a time window length is set; The historical trajectory curvature of the historical initial wavelength offset estimation value in each time window is calculated, and the time window whose historical trajectory curvature exceeds the preset curvature threshold is marked as a curvature abnormal period; The initial wavelength offset estimation value in the curvature abnormal period is extracted to form an abnormal wavelength offset sub-sequence; The correlation coefficient of the abnormal wavelength offset sub-sequence and the current initial wavelength offset estimation value is calculated, the abnormal wavelength offset sub-sequence with a correlation coefficient greater than a preset correlation coefficient threshold is screened, and a wavelength offset mode is determined; The wavelength offset mode is subjected to frequency domain transformation to obtain an amplitude spectrum and a phase spectrum, the frequency component corresponding to the amplitude peak value in the amplitude spectrum is extracted as a main frequency feature, and the phase gradient change rate exceeding a preset phase gradient threshold is extracted from the phase spectrum as a phase mutation feature; The main frequency feature and the phase mutation feature are combined to form the frequency domain feature component.

[0013] In a second aspect of the embodiment of the application, a wavelength deviation compensation system for fast analysis of a remote fiber grating wavelength is provided, comprising: A first unit is configured to collect reflection spectrum signals of a plurality of fiber gratings experiencing different transmission distances; A second unit is configured to perform frequency domain transformation on the reflection spectrum signals, separate linear phase components and nonlinear phase components, establish a linear mapping relationship between transmission distance and dispersion accumulation, calculate transmission dispersion offset and intrinsic wavelength offset, and obtain initial wavelength offset estimation values of gratings at different spatial positions; A third unit is configured to identify refractive index mutation sections and refractive index stable sections according to the initial wavelength offset estimation values, establish segmented compensation coefficients, calculate correction amounts through tensor decomposition, and correct the initial wavelength offset estimation values to obtain compensated grating center wavelengths; The fourth unit is configured to construct a wavelength-strain-temperature three-dimensional space, map the compensated grating center wavelength to a track in the three-dimensional space, calculate a track curvature, and when the track curvature exceeds a preset curvature threshold, reconstruct a segmented compensation coefficient according to a historical initial wavelength offset estimation value sequence, and output physical quantity demodulation data. The fifth unit is configured to map a wavelength reconstruction residual in the physical quantity demodulation data to the three-dimensional space, and update the segmented compensation coefficient according to a residual gradient direction.

[0014] In a third aspect, the present application provides an electronic device, comprising: a processor; a memory for storing processor-executable instructions; The processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0015] In a fourth aspect, the present application provides a computer-readable storage medium having computer program instructions stored thereon, wherein the computer program instructions are executed by a processor to implement the method described above.

[0016] In the present application, the linear phase component and the nonlinear phase component are separated by frequency domain transformation, a linear mapping relationship between the transmission distance and the dispersion accumulation is established, the accurate demodulation of the remote grating center wavelength is realized, the influence of the dispersion caused by the long-distance transmission on the grating measurement is effectively eliminated, and the monitoring accuracy of the fiber grating sensing network is improved. BRIEF DESCRIPTION OF DRAWINGS

[0017] Figure 1 A flowchart of the wavelength deviation compensation method for the remote fiber grating wavelength rapid analysis of the present application embodiment; Figure 2 A fiber sensing signal compensation and separation flowchart. DETAILED DESCRIPTION

[0018] In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the protection scope of the present application.

[0019] The technical solutions of the present application will be described in detail below with specific embodiments. The following specific embodiments can be combined with each other, and the same or similar concepts or processes can not be described in some embodiments.

[0020] Figure 1 The flowchart of the wavelength deviation compensation method for the remote fiber grating wavelength rapid analysis of the embodiments of the present application is shown as Figure 1 The method comprises: Collecting reflection spectrum signals of a plurality of fiber gratings with different transmission distances; Performing frequency domain transformation on the reflection spectrum signals, separating linear phase components and nonlinear phase components, establishing a linear mapping relationship between transmission distance and dispersion accumulation, calculating transmission dispersion offset and intrinsic wavelength offset, and obtaining initial wavelength offset estimation values of gratings at different spatial positions; According to the initial wavelength offset estimation values, identifying refractive index sudden change sections and refractive index stable sections, establishing segmented compensation coefficients, calculating correction values through tensor decomposition, correcting the initial wavelength offset estimation values, and obtaining compensated grating center wavelengths; Constructing a wavelength-strain-temperature three-dimensional space, mapping the compensated grating center wavelengths into a trajectory in the three-dimensional space, calculating the trajectory curvature, when the trajectory curvature exceeds a preset curvature threshold, reconstructing the segmented compensation coefficients according to a historical initial wavelength offset estimation value sequence, and outputting physical quantity demodulation data; Mapping the wavelength reconstruction residual in the physical quantity demodulation data to the three-dimensional space, and updating the segmented compensation coefficients according to the residual gradient direction.

[0021] In an optional implementation, performing frequency domain transformation on the reflection spectrum signals, separating linear phase components and nonlinear phase components, establishing a linear mapping relationship between transmission distance and dispersion accumulation, calculating transmission dispersion offset and intrinsic wavelength offset, and obtaining initial wavelength offset estimation values of gratings at different spatial positions comprises: Transforming the reflection spectrum signals from time domain to frequency domain to obtain frequency domain spectrum distribution; A continuous phase response curve is extracted from the frequency domain spectrum distribution, and a polynomial fitting is performed on the phase response curve according to a fitting weight, so as to separate a linear phase component composed of linear fitting parameters and a nonlinear phase component composed of a residual phase distribution; According to the linear phase component, a group delay difference corresponding to different transmission distances is calculated, and a linear mapping relationship between the transmission distance and the dispersion accumulation is established. The nonlinear phase component is converted back to the wavelength domain, and the intrinsic center wavelength of each grating position is extracted. According to the linear mapping relationship, the group delay difference is converted into the dispersion accumulation of each spatial position, and the transmission dispersion offset and the intrinsic wavelength offset of each grating are calculated in combination with the intrinsic center wavelength. The transmission dispersion offset and the intrinsic wavelength offset are vector synthesized to obtain the initial wavelength offset estimation value of the grating at different spatial positions.

[0022] In a specific embodiment, the reflection spectrum signal of the fiber grating sensor is obtained, which is usually the corresponding relationship data of wavelength and intensity obtained by the grating reflected light through a spectrum analyzer or other equipment. The collected reflection spectrum signal is transformed from time domain to frequency domain through fast Fourier transform algorithm. For example, the Fourier transform is performed on the reflection spectrum signal with a sampling frequency of 10 GHz and a data point number of 1024, and the corresponding frequency domain spectrum distribution is obtained.

[0023] After obtaining the frequency domain spectrum distribution, it is necessary to extract the phase response information therefrom. By calculating the amplitude and phase of the frequency domain signal, a group of phase values distributed on the frequency axis can be obtained. Since the phase value usually jumps between -π and π, phase unwrapping processing is required to make the phase curve change continuously. In experimental tests, the grating reflection signal in the wavelength range of 1530 nm to 1560 nm is processed, and the continuous phase curve is obtained after unwrapping.

[0024] The continuous phase response curve after unwrapping is polynomial fitted, and the signal intensity of different frequency points is considered as the fitting weight during the fitting, and the frequency point with higher signal intensity has higher weight. The fitting adopts the least square method, and a 3-order polynomial is used for fitting. In the fitting result, the first-order term coefficient corresponds to the linear phase component, representing the group delay characteristic; the remaining high-order terms and the fitting residual constitute the nonlinear phase component, reflecting the intrinsic characteristics of the grating and the wavelength shift caused by dispersion. In actual tests, the linear phase component coefficient is about -2.5×10 -10 s, representing the basic delay characteristic of the optical signal in the transmission medium.

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

[0026] The nonlinear phase component is converted back to the wavelength domain through inverse Fourier transform to extract the intrinsic center wavelength of each grating position. In the wavelength domain, the nonlinear phase component appears as distortion of the wavelength spectrum, and the peak position corresponds to the intrinsic center wavelength of the grating. Through the peak search algorithm, the center wavelength of each grating can be accurately determined. For example, in the test, the intrinsic center wavelengths of the four gratings at different positions are 1541.23 nm, 1545.67 nm, 1550.12 nm and 1555.48 nm, respectively.

[0027] The group delay difference is converted to the dispersion accumulation of each spatial position using the established linear mapping relationship. For example, for the grating with a distance of 20 km, the dispersion accumulation corresponding to the 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 the dispersion effect, and the difference between the center wavelength and the reference wavelength and the accumulated dispersion amount need to be considered in the calculation. For the grating near 1550 nm, the wavelength offset caused by 20 km of transmission distance is about 0.05 nm.

[0028] While determining the transmission dispersion offset, the intrinsic wavelength offset is extracted by analyzing the nonlinear phase component. The intrinsic wavelength offset represents the wavelength change caused by the characteristics of the grating itself or external stress, temperature and other factors. In the experiment, a temperature change of 1°C causes an intrinsic wavelength offset of about 0.01 nm.

[0029] The transmission dispersion offset and the intrinsic wavelength offset are vectorially combined to obtain the initial wavelength offset estimate of the grating at different spatial positions. The vectorial combination considers the directionality of the two offsets, and the offset for increasing the wavelength is positive, and the offset for decreasing the wavelength is negative. Through experimental verification, the error between the wavelength offset value estimated by the method and the actual measurement value is within ±0.003 nm, which meets the requirements of high-precision distributed sensing.

[0030] In an alternative embodiment, a continuous phase response curve is extracted from the frequency domain spectral distribution, a polynomial fitting is performed on the phase response curve according to fitting weights, and a linear phase component composed of linear fitting parameters and a nonlinear phase component composed of residual phase distribution are separated, including: A phase unwrapping process is performed on the complex amplitudes of the frequency points in the frequency domain spectral distribution to generate a continuous phase response curve; Fitting weights are assigned to the sampling points in the continuous phase response curve according to amplitude information, a polynomial fitting is performed on the continuous phase response curve according to the fitting weights, and linear fitting parameters representing a frequency first power relationship and nonlinear fitting parameters representing a frequency high power relationship are separated; The linear fitting parameters are multiplied by the frequency points to obtain a linear phase component corresponding to the transmission distance, and the linear phase component is subtracted from the continuous phase response curve to obtain a residual phase distribution; The ratio of the phase difference value of adjacent frequency points to the frequency interval in the residual phase distribution is calculated to generate a phase change rate sequence, the direction of grating period change is determined through the positive and negative value distribution statistics of the phase change rate sequence, the grating period gradual change rate is extracted through the adjacent difference value of the phase change rate sequence, and the residual phase distribution is structure feature labeled according to the direction of grating period change and the grating period gradual change rate; The residual phase distribution labeled with structure features is determined as the nonlinear phase component corresponding to the grating structure.

[0031] In a specific embodiment, when performing phase unwrapping processing on the complex amplitudes of the frequency points in the frequency domain spectral distribution, frequency domain spectral distribution data is obtained, which contains frequency point information and corresponding complex amplitude information. The complex amplitude can be represented in the form of amplitude and phase, where the initial phase value is usually limited to the range of -π to π. Since the phase value has a periodic jump, phase unwrapping processing is needed to obtain a continuous phase response curve. The phase unwrapping processing uses a phase accumulation method, which detects the phase difference value between adjacent frequency points. When the phase difference value is detected to exceed π, it is determined as a phase jump point, and all phase values after this point are adjusted by an appropriate 2π integer multiple to make the overall phase curve show continuous change characteristics. For example, when the phase value at frequency point f1 is 2.9π and the phase value at adjacent frequency point f2 is -2.95π, the actual phase difference is about 0.05π instead of 5.85π. After phase unwrapping processing, the phase of f2 point is adjusted to 3.05π to ensure the continuity of the phase curve.

[0032] When assigning fitting weights to each sampling point in the continuous phase response curve according to the amplitude information, considering that the signal-to-noise ratios of different frequency points in the frequency domain spectrum distribution are different, the frequency points with larger signal amplitudes usually have higher phase accuracy. Therefore, a weight assignment strategy positively correlated with the amplitudes of the frequency points is adopted, and the weight coefficient can be set to be proportional to the square of the amplitude. For example, for a frequency point with an amplitude of A, the weight corresponding to the frequency point can be set as A 2 / A_max 2 , where A_max is the maximum amplitude of all frequency points. In specific operations, a threshold value of the amplitude can be set, and when the amplitude of a frequency point is lower than the threshold value, the weight of the frequency point is set to zero, so as to exclude the interference of noise on the fitting result.

[0033] When performing polynomial fitting on the continuous phase response curve according to the fitting weights, the weighted least squares method is used for polynomial fitting, and the polynomial order can be selected as 3-5. In the fitting process, the frequency is taken as the independent variable, and the phase value is taken as the dependent variable, to construct a polynomial fitting model. By minimizing the weighted residual sum of squares, the polynomial coefficients are solved. Among the solved polynomial coefficients, the coefficient of the first order term corresponds to the linear phase characteristic and is directly related to the transmission distance; the coefficients of the high order terms represent the nonlinear relationship between the frequency and the phase, and are related to the characteristics of the grating structure. For example, for certain test data, after 4-order polynomial fitting, the coefficients a0=-0.12, a1=9.83×10 -5 , a2=2.15×10 -10 , a3=-4.37×10 -15 , and a4=1.92×10 -20 are obtained, where a1 is the linear fitting parameter.

[0034] When multiplying the linear fitting parameter with each frequency point, the linear parameter a1 is directly multiplied by the corresponding frequency value to obtain the linear phase component. For example, for a point with a frequency of f, the linear phase component is a1×f. The linear phase component is subtracted from the continuous phase response curve, and the remaining part is the residual phase distribution, which contains the nonlinear phase characteristics related to the grating structure.

[0035] When calculating the ratio of the phase difference between adjacent frequency points in the residual phase distribution and the frequency interval, 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 the calculated phase change rates are arranged in frequency order to form a sequence of phase change rates. By counting the number distribution of positive and negative values in the sequence, the change direction of the grating period is determined. When the number of positive values is significantly more than the number of negative values, it indicates that the grating period is increasing; otherwise, it indicates that the grating period is decreasing. The difference between adjacent points in the sequence of phase change rates is calculated to obtain the second-order phase change rate, and the numerical distribution characteristics of the second-order phase change rate reflect the grating period gradient rate. For example, when the mean value of the second-order phase change rate is positive and shows an increasing trend, it indicates that the grating period gradient rate is accelerating; when the mean value of the second-order phase change rate is close to zero and the fluctuation is small, it indicates that the grating period changes at an approximately constant speed.

[0036] When the change direction of the grating period and the grating period gradient rate are used to label the structural characteristics of the residual phase distribution, the residual phase distribution is compared with the theoretical model to determine the grating structure type. For example, when the residual phase distribution shows a quadratic function characteristic and the grating period changes linearly, it can be labeled as a chirped grating structure; when the residual phase distribution shows a high-order polynomial characteristic and the grating period gradient rate changes nonlinearly, it can be labeled as a non-uniform chirped grating structure. Through structural characteristic labeling, a clear correspondence between the residual phase distribution and the grating structure is established, and the nonlinear phase component corresponding to the grating structure is obtained.

[0037] In an optional implementation, the refractive index abrupt change section and the refractive index stable section are identified according to the initial wavelength shift estimation value, and a segmented compensation coefficient is established. The initial wavelength shift estimation value is corrected by tensor decomposition to obtain the compensated grating center wavelength, including: The initial wavelength shift estimation value is segmented and counted according to the transmission distance, and the gradient change rate of the initial wavelength shift estimation value in each transmission distance section is calculated. The transmission distance section with a gradient change rate exceeding a preset change rate threshold is marked as a refractive index abrupt change section, and the transmission distance section with a gradient change rate within the change rate threshold is marked as a refractive index stable section. The boundary positions of the refractive index abrupt change section and the refractive index stable section are extracted, and a segmented compensation coefficient is established at the boundary positions. Principal component analysis is performed on the temperature field component and the strain field component 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. The first characteristic vector is multiplied with the segment compensation coefficient element by element to obtain a transmission dispersion correction quantity; and the second characteristic vector is tensor outer product with the third characteristic vector to obtain a multi-physical field coupling correction quantity. The transmission dispersion correction quantity and the multi-physical field coupling correction quantity are subtracted from the initial wavelength shift estimation value to obtain a compensated grating center wavelength.

[0038] In a specific embodiment, the collected wavelength shift data is affected by dispersion effect and environmental disturbance during 100 km transmission, resulting in deviation of the wavelength shift data. After obtaining the initial wavelength shift estimation value, it needs to be segmented. The initial wavelength shift estimation value is segmented according to every 5 km of transmission distance, and the gradient change rate of the initial wavelength shift estimation value in each segment is calculated. The gradient change rate calculation method is to calculate the ratio of the wavelength shift value difference of adjacent measurement points in the segment to the distance difference, to obtain the gradient change rate sequence in each segment, and then calculate the standard deviation. When the standard deviation exceeds the preset change rate threshold of 0.01 nm / km, the segment is marked as a refractive index mutation segment; when the standard deviation is within 0.01 nm / km, the segment is marked as a refractive index stable segment. For example, in actual application, the 0-5 km segment gradient change rate standard deviation is 0.008 nm / km, which is marked as a refractive index stable segment; the 5-10 km segment gradient change rate standard deviation is 0.025 nm / km, which is marked as a refractive index mutation segment.

[0039] The boundary positions of the refractive index mutation segment and the refractive index stable segment are extracted, and the segment compensation coefficient is established at the boundary position. For example, there is a boundary at 5 km, and the compensation coefficient is established at this position. The compensation coefficient of the refractive index stable segment is set to 1.0, and the compensation coefficient of the refractive index mutation segment is determined by the ratio of the local wavelength shift to the ideal wavelength shift. For example, the compensation coefficient of the 5-10 km segment is calculated as the ratio of the average measured wavelength shift value of the segment to the ideal wavelength shift value under the same temperature strain condition, and the compensation coefficient is 1.23. After all the boundary positions establish the corresponding segment compensation coefficient, a complete segment compensation coefficient table is formed.

[0040] The principal component analysis is performed on the temperature field component and the strain field component to determine the tensor dimension size. The collected temperature data forms a temperature matrix, and the strain data forms a strain matrix, and the eigenvalues are calculated respectively. The number of characteristic vectors corresponding to the eigenvalues with a cumulative contribution rate of 95% is selected as the respective dimension size. For example, the temperature field dimension is 4, and the strain field dimension is 3. According to the determined dimension size, a three-order tensor structure is constructed, the first dimension is the transmission distance (set as 20 measurement points), the second dimension is the temperature field dimension (4), and the third dimension is the strain field dimension (3), forming a 20×4×3 three-order tensor. The initial wavelength shift estimation value is filled into the tensor to form a complete data tensor.

[0041] An alternating least squares iteration operation with constraints is performed to optimize the three eigenvectors. The upper limit of iteration is set to 100 times and the convergence threshold is set to 0.0001. In each iteration, the eigenvectors of two dimensions are fixed, and the eigenvector of the third dimension is optimized. The optimization is performed in turn until the convergence condition is reached or the upper limit of iteration is reached. Through this process, the first eigenvector of the transmission distance dimension (20-dimensional vector), the second eigenvector of the temperature field dimension (4-dimensional vector), and the third eigenvector of the strain field dimension (3-dimensional vector) are obtained.

[0042] The first eigenvector is multiplied element by element with the segmented compensation coefficient to obtain the transmission dispersion correction. Specifically, each element in the 20-dimensional first eigenvector is multiplied by the segmented compensation coefficient at the corresponding transmission distance. For example, the 6th measurement point is located in the refractive index mutation section, and the first eigenvector element value is 0.15 nm. The corresponding segmented compensation coefficient is 1.23, and the calculated transmission dispersion correction is 0.1845 nm. The second eigenvector and the third eigenvector are subjected to tensor outer product operation to obtain the multi-physical field coupling correction. The tensor outer product operation is to multiply the 4-dimensional second eigenvector and the 3-dimensional third eigenvector to generate a 4x3 matrix. The matrix represents the influence of the coupling effect of the temperature field and the strain field on the wavelength shift. The sum of the matrix elements is the total multi-physical field coupling correction, for example, the calculation result is 0.0722 nm.

[0043] The transmission dispersion correction and the multi-physical field coupling correction are subtracted from the initial wavelength shift estimation value to obtain the compensated grating center wavelength. For example, the initial wavelength shift estimation value of a certain measurement point is 1.2563 nm, the transmission dispersion correction is 0.1845 nm, and the multi-physical field coupling correction is 0.0722 nm. The compensated grating center wavelength is 1.2563 nm-0.1845 nm-0.0722 nm=0.9996 nm. The same compensation calculation is performed for 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, and the relative error is reduced by 88%, effectively improving the measurement accuracy in long-distance applications.

[0044] As shown in FIG. 1, a flowchart of the optical fiber sensing signal compensation and separation is shown. Figure 2

[0045] ​In an alternative embodiment, the principal component analysis is performed on the temperature field component and the strain field component to determine the tensor dimension size, a third-order tensor structure is constructed and filled, a constrained alternating least squares iterative operation is performed, and 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 are separated and obtained, including: The principal component analysis is performed on the temperature field component and the strain field component, respectively, the corresponding principal component variance contribution rates are extracted, the cumulative sum of the principal component variance contribution rates is calculated, the cumulative variance contribution rate is determined, and when the cumulative variance contribution rate reaches a preset variance threshold, the number of principal components corresponding to the cumulative variance contribution rate is determined as the temperature field dimension size and the strain field dimension size. A third-order tensor structure is constructed according to the number of sampling points of the initial wavelength shift estimation value, the temperature field dimension size, and the strain field dimension size. The three-dimensional coordinate index of the initial wavelength shift estimation value in the third-order tensor structure is traversed, when the tensor element corresponding to the three-dimensional coordinate index does not exist, the filled tensor element with a Euclidean distance less than a preset distance threshold from the three-dimensional coordinate index is searched, the Euclidean distance reciprocal of each filled tensor element is calculated to determine the interpolation weight, each filled tensor element is weighted and summed according to the interpolation weight, and the three-dimensional coordinate index position is filled to obtain a filled third-order tensor structure. The alternating least squares iterative operation with non-negativity constraint and orthogonality constraint is performed on the filled third-order tensor structure, the iteration is terminated when the relative change rate of the third-order tensor reconstruction error is less than a preset convergence threshold, and the factor matrix of each dimension is extracted to determine the first eigenvector, the second eigenvector, and the third eigenvector.

[0046] In a specific embodiment, the principal component analysis is performed on the temperature field component and the strain field component to determine the appropriate dimension size. For the temperature field component, data samples including multiple temperature points such as 20℃, 30℃, 40℃, 50℃, and 60℃ are collected. By calculating the covariance matrix and performing eigenvalue decomposition, the corresponding eigenvalues are obtained as [3.75, 0.82, 0.31, 0.09, 0.03]. The variance contribution rates of each principal component are calculated as [0.75, 0.16, 0.06, 0.02, 0.01], and the cumulative variance contribution rates are calculated as [0.75, 0.91, 0.97, 0.99, 1.00]. The preset variance threshold is set to 0.95, and when the cumulative contribution rate reaches 0.97, the number of principal components corresponding to the cumulative variance contribution rate is 3, and therefore the temperature field dimension size is determined to be 3.

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

[0048] Based on the above analysis results, a three-order tensor structure is constructed. The first dimension of the tensor represents the transmission distance sampling point, and the size is the number of sampling points of the initial wavelength offset estimation value 10; the second dimension represents the temperature field, and the size is 3; and the third dimension represents the strain field, and the size is 4. Therefore, a 10×3×4 three-order tensor structure is constructed.

[0049] For this tensor structure, missing data elements need to be filled. For example, there is no actual measurement data at the coordinate (5, 2, 3). The preset distance threshold is set to 1.5, and the filled tensor elements with a Euclidean distance less than 1.5 from the coordinate are searched, and four points (4, 2, 3), (5, 1, 3), (5, 2, 2) and (6, 2, 3) are found, whose Euclidean distances are 1.0, 1.0, 1.0 and 1.0, the distance reciprocals are 1.0, 1.0, 1.0 and 1.0, and the normalized interpolation weights are all 0.25. Assuming that the values of the four filled elements are 0.82, 0.75, 0.78 and 0.85 respectively, the filling value at (5, 2, 3) is obtained by weighted summation according to the interpolation weights as 0.8. Through a similar method, the filling of the entire tensor is completed.

[0050] The three-order tensor after filling is subjected to alternating least squares iteration operation with non-negativity constraint and orthogonality constraint. The preset convergence threshold is set to 10 -4 , and three factor matrices are initialized, wherein 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, and all elements are initialized to random non-negative values.

[0051] In the first iteration, the second and third factor matrices are fixed, and the first factor matrix is updated. The obtained 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 to obtain [0.51, 0.63, 0.58]. The first and second factor matrices are fixed, and the third factor matrix is updated to obtain [0.45, 0.52, 0.56, 0.47]. The tensor reconstruction error is calculated to be 0.085.

[0052] In the second iteration, the above process is repeated to obtain 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 tensor reconstruction error is calculated to be 0.078, and the relative change rate is (0.085-0.078) / 0.085≈0.082.

[0053] The iteration continues until the 15th iteration, and the three factor matrices converge to the first factor matrix [0.30, 0.31, 0.32, 0.33, 0.34, 0.33, 0.32, 0.31, 0.30, 0.28], the second factor matrix [0.54, 0.60, 0.59], and the third factor matrix [0.43, 0.53, 0.55, 0.49]. At this time, the tensor reconstruction error is 0.056, and the relative change rate is 0.0009, which is less than the preset convergence threshold 10 -4 Therefore, the iteration is terminated.

[0054] 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 sensing data analysis and processing.

[0055] In an optional embodiment, a wavelength-strain-temperature three-dimensional space is constructed, the compensated grating center wavelength is mapped as a trajectory in the three-dimensional space, the trajectory curvature is calculated, when the trajectory curvature exceeds a preset curvature threshold, the segmented compensation coefficient is reconstructed according to a historical initial wavelength offset estimation value sequence, and the physical quantity demodulation data is output, which includes: A three-dimensional space is established with wavelength, strain, and temperature as coordinate axes, and the compensated grating center wavelength, strain value, and temperature value at the corresponding time are combined to form a three-dimensional coordinate point sequence, and a trajectory in the three-dimensional space is formed; The three-dimensional coordinate point sequence is subjected to differential operation to obtain a trajectory tangent vector, and the angle change rate of adjacent trajectory tangent vectors is calculated to determine the trajectory curvature; When the trajectory curvature exceeds a preset curvature threshold, a wavelength shift mode corresponding to a curvature abnormal period in the historical initial wavelength shift estimation value sequence is identified, and a frequency domain feature component of the wavelength shift mode is extracted; The current initial wavelength shift estimation value is matched with the frequency domain feature component, a compensation weight of different refractive index sections is assigned according to the matching degree, and the segmented compensation coefficient is reconstructed; The reconstructed segmented compensation coefficient is applied to the compensated grating center wavelength, the wavelength shift caused by strain and the wavelength shift caused by temperature are separated, and physical quantity demodulation data is output.

[0056] In a specific embodiment, on the basis of obtaining the compensated grating center wavelength, a three-dimensional space is established with wavelength, strain, and temperature as coordinate axes, for further analyzing the wavelength shift characteristics. The compensated grating center wavelength, strain value, and temperature value at the corresponding time are combined to form a three-dimensional coordinate point sequence, and a trajectory in the three-dimensional space is formed. For example, at a certain time, if the compensated grating center wavelength is 1550.214 nanometers, 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. At a frequency of 10 samples per second, 10 minutes of continuous sampling can obtain 6000 such three-dimensional coordinate points, forming a complete trajectory.

[0057] The three-dimensional coordinate point sequence is subjected to differential operation to obtain a trajectory tangent vector. The differential operation is realized by using the central difference method. For the i-th point in the sequence, the tangent vector is equal to the coordinate difference between the i+1-th point and the i-1-th point divided by the corresponding time interval. An appropriate time window size is selected for smoothing processing, and a typical value is 5 to 10 sampling points. For example, for the three-dimensional coordinate point (1550.214, 120, 25.3) at a certain time, if the previous coordinate is (1550.212, 119.5, 25.2) and the next coordinate is (1550.217, 120.8, 25.4), the sampling interval is 0.1 seconds, then the tangent vector of the point is ((1550.217-1550.212) / 0.2, (120.8-119.5) / 0.2, (25.4-25.2) / 0.2), i.e. (0.025, 6.5, 1).

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

[0059] When the trajectory curvature exceeds a preset curvature threshold, it indicates that the relationship between the grating center wavelength and strain or temperature has changed abnormally, which can be caused by external interference or changes in the state of the optical fiber. The preset curvature threshold is set according to the actual application environment, and is usually between 0.03 and 0.1 rad / s. When the curvature exceeds the threshold, the wavelength shift pattern corresponding to the period of abnormal curvature in the historical initial wavelength shift estimation value sequence needs to be identified. For example, if in a 10-minute sampling sequence, the trajectory curvature from the 245th second to the 268th second continuously exceeds the threshold of 0.05 rad / s, reaching 0.07 rad / s, then the wavelength shift estimation value sequence in this period is extracted as the abnormal pattern.

[0060] The frequency domain feature components of the wavelength shift pattern are extracted, and the time domain wavelength shift sequence is converted to the frequency domain using the Fast Fourier Transform. To obtain more accurate spectral characteristics, a Hanning window function can be applied to reduce the spectral leakage effect. For the above-mentioned 23-second abnormal wavelength shift pattern, assuming a sampling rate of 10 Hz, there are 230 sampling points. After Fourier transform, it can be found that there are obvious amplitude peaks at 2 Hz and 5 Hz, which are 0.015 nm and 0.008 nm respectively. These frequencies and amplitudes form the frequency domain feature components of the wavelength shift pattern.

[0061] The current initial wavelength shift estimation value is matched with the identified frequency domain feature components, and the matching degree is calculated. The matching degree is calculated based on the frequency domain correlation of the current wavelength shift sequence and the historical abnormal pattern. Taking the last 10 seconds of wavelength shift data, the frequency domain representation is also obtained by Fourier transform, and the amplitudes at 2 Hz and 5 Hz are calculated, which are assumed to be 0.014 nm and 0.007 nm. Compared with the frequency domain feature components of the historical abnormal pattern, the matching degree can be represented as the normalized value of the frequency point amplitude difference, such as (1-|0.015-0.014| / 0.015) x (1-|0.008-0.007| / 0.008) = 0.94.

[0062] The compensation weight of different refractive index sections is assigned according to the matching degree, and the segmented compensation coefficient is reconstructed. The higher the matching degree is, the more similar the interference mode of the current wavelength shift estimation value is to the historical abnormal mode, and the stronger the compensation adjustment needs to be. Assuming that the original segmented compensation coefficient is 0.92, 1.08 and 0.95 at 2500 meters, 5500 meters and 8000 meters respectively, and the matching degree is 0.94, the compensation coefficient can be adjusted according to the matching degree, and the new compensation coefficient 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 coefficient is 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.

[0063] The reconstructed segmented compensation coefficient is applied to the compensated grating center wavelength, further improving the accuracy of wavelength shift compensation. Taking the point with a transmission distance of 5000 meters as an example, if the compensated grating center wavelength is 1550.214 nanometers, the point is located in the second refractive index section, and the corresponding reconstructed compensation coefficient is 1.133, then the final compensation result is 1550.214×1.133=1756.393 nanometers.

[0064] Based on the accurately compensated grating center wavelength, the wavelength drift caused by strain and the wavelength drift caused by temperature are separated. This separation process takes advantage of the different characteristics of strain and temperature on wavelength. Generally, the wavelength drift caused by temperature change has a lower frequency characteristic, while the wavelength drift caused by strain change has a relatively higher frequency characteristic. Through a band-pass filter set, wavelength drifts of different frequency bands can be separated. The low-frequency part (0-0.5 Hz) mainly corresponds to the drift caused by temperature, and the high-frequency part (0.5-5 Hz) mainly corresponds to the drift caused by strain. For example, for the final compensation result 1756.393 nanometers, after frequency separation, the wavelength drift caused by temperature can be 1.245 nanometers, and the wavelength drift caused by strain can be 0.523 nanometers.

[0065] According to the wavelength sensitivity parameters of the fiber Bragg grating, the separated wavelength drifts are converted into actual physical quantities. Assuming that the temperature sensitivity is 10 picometers per degree Celsius and the strain sensitivity is 1.2 picometers per micro-strain, then the wavelength drift of 1.245 nanometers caused by temperature corresponds to a temperature change of 124.5 degrees Celsius, and the wavelength drift of 0.523 nanometers caused by strain corresponds to a strain change of 435.8 micro-strains. These converted physical quantities are the final demodulation data, which can be used for subsequent monitoring analysis. Through this wavelength deviation compensation and physical quantity separation method, the accuracy of remote fiber grating wavelength analysis can be improved to within 0.005 nanometers, and the corresponding temperature and strain analysis accuracy can reach 0.5 degrees Celsius and 5 micro-strains respectively.

[0066] In an optional embodiment, when the trajectory curvature exceeds the preset curvature threshold, identifying the wavelength shift pattern corresponding to the curvature abnormal period in the sequence of historical initial wavelength shift estimation values comprises: When the trajectory curvature exceeds the preset curvature threshold, a time window sliding mechanism is established in the sequence of historical initial wavelength shift estimation values, and a time window length is set; The historical trajectory curvature corresponding to the historical initial wavelength shift estimation values in each time window is calculated, and the time window in which the historical trajectory curvature exceeds the preset curvature threshold is marked as a curvature abnormal period; The initial wavelength shift estimation values in the curvature abnormal period are extracted to form an abnormal wavelength shift sub-sequence; The correlation coefficient of the abnormal wavelength shift sub-sequence and the current initial wavelength shift estimation value is calculated, the abnormal wavelength shift sub-sequence with a correlation coefficient greater than a preset correlation coefficient threshold is screened, and a wavelength shift pattern is determined; The wavelength shift pattern is subjected to frequency domain transformation to obtain an amplitude spectrum and a phase spectrum, the frequency component corresponding to the amplitude peak value in the amplitude spectrum is extracted as a main frequency feature, and the phase gradient change rate exceeding a preset phase gradient threshold is extracted from the phase spectrum as a phase mutation feature; The main frequency feature and the phase mutation feature are combined to form the frequency domain feature component.

[0067] In a specific embodiment, a time window sliding mechanism is established to locate the curvature abnormal period in the sequence of historical initial wavelength shift estimation values, then the wavelength shift pattern is determined through correlation analysis, and finally the feature component is extracted through frequency domain transformation.

[0068] The method of the present application triggers the process of identifying the wavelength shift pattern by monitoring the trajectory curvature in real time. When the trajectory curvature exceeding the preset curvature threshold (for example, 0.85 rad / s 2 ) is detected, the time window sliding mechanism is started to be performed on the sequence of historical initial wavelength shift estimation values. The sequence of historical initial wavelength shift estimation values contains wavelength shift estimation value data collected at a sampling frequency of 5 Hz in the past 30 minutes. The system sets the time window length to 60 seconds, and the window sliding step length to 5 seconds.

[0069] For each time window, the trajectory curvature corresponding to the historical initial wavelength shift estimation values in the window is calculated. Specifically, for the wavelength shift estimation values λt and λt+1 of each two adjacent time points t and t+1 in the window, the difference Δλ = λt+1 - λt is calculated, and then the trajectory curvature is calculated according to the change rate of the adjacent difference. If the maximum value of the trajectory curvature in a certain time window exceeds the preset curvature threshold 0.85 rad / s2 Then, the time window is marked as a curvature anomaly period.

[0070] For example, when the current trajectory curvature is detected as 0.92 rad / s 2 , which exceeds the preset curvature threshold. By the time window sliding mechanism, the system identifies 3 curvature anomaly periods in the historical data, which are located at [5min, 6min], [12min, 13min] and [22min, 23min] time periods within the past 30 minutes.

[0071] From these curvature anomaly periods, initial wavelength shift estimation values are extracted, which form anomaly wavelength shift sub-sequences A, B and C, respectively. Sub-sequence A contains 60 seconds x 5Hz = 300 sampling points, with wavelength shift values ranging from [1532.25nm, 1532.46nm]; sub-sequence B also contains 300 sampling points, with wavelength shift values ranging from [1532.31nm, 1532.52nm]; sub-sequence C also contains 300 sampling points, with wavelength shift values ranging from [1532.28nm, 1532.49nm].

[0072] The correlation coefficients of these anomaly wavelength shift sub-sequences with the current initial wavelength shift estimation value sequence (300 sampling points within the last 60 seconds, with wavelength shift values ranging from [1532.29nm, 1532.50nm]) are calculated. The correlation coefficient of sub-sequence A with the current sequence is 0.72, the correlation coefficient of sub-sequence B is 0.91, and the correlation coefficient of sub-sequence C is 0.88. If the preset correlation coefficient threshold is set as 0.80, then sub-sequences B and C are selected as the wavelength shift pattern.

[0073] For the determined wavelength shift pattern, frequency domain transformation is performed to obtain the amplitude spectrum and phase spectrum. Preferably, the Fast Fourier Transform algorithm is used to convert 300 time domain sampling points into frequency domain representation. For sub-sequence B, frequency domain transformation is performed, and the obtained amplitude spectrum shows significant peaks at 0.25Hz, 0.75Hz and 1.2Hz, with peak amplitudes of 0.056, 0.038 and 0.027, respectively; the phase spectrum shows phase values at these frequency points as 0.45π, 1.2π and 1.7π, respectively. For sub-sequence C, the same operation is performed, and the obtained amplitude spectrum shows significant peaks at 0.28Hz, 0.72Hz and 1.25Hz, with peak amplitudes of 0.052, 0.041 and 0.023, respectively; the phase spectrum shows phase values at these frequency points as 0.42π, 1.25π and 1.75π, respectively.

[0074] The frequency components corresponding to the amplitude peak values are extracted from the amplitude spectrum as the main frequency features. In the sub-sequence B, the main frequency features are [0.25 Hz, 0.75 Hz, 1.2 Hz]; in the sub-sequence C, the main frequency features are [0.28 Hz, 0.72 Hz, 1.25 Hz].

[0075] The phase gradient of the phase spectrum, i.e., the phase change rate between adjacent frequency points, is calculated. In the sub-sequence B, the phase gradient in the frequency range [0.70 Hz, 0.80 Hz] is 2.5π / Hz, which exceeds the preset phase gradient threshold value 2.0π / Hz; in the sub-sequence C, the phase gradient in the frequency range [0.68 Hz, 0.78 Hz] is 2.6π / Hz, which also exceeds the preset threshold value. The frequency components with the phase gradient change rate exceeding the threshold value are extracted as the phase mutation features, which are [0.75 Hz] and [0.72 Hz], respectively.

[0076] The main frequency features and the phase mutation features are combined to form the frequency domain feature components. For the sub-sequence B, the frequency domain feature components are {main frequency: [0.25 Hz, 0.75 Hz, 1.2 Hz], phase mutation: [0.75 Hz]}; for the sub-sequence C, the frequency domain feature components are {main frequency: [0.28 Hz, 0.72 Hz, 1.25 Hz], phase mutation: [0.72 Hz]}. Through these frequency domain feature components, the wavelength shift abnormal pattern can be accurately identified, which provides an important basis for subsequent wavelength correction.

[0077] In this embodiment, the extracted frequency domain feature components are used for adjustment of the wavelength compensation model, so that the system can perform targeted correction for specific types of wavelength shift patterns, and the measurement accuracy and stability of the optical fiber sensing system in a complex environment are significantly improved.

[0078] The wavelength deviation compensation system for remote fiber grating wavelength rapid analysis in the embodiment of the application comprises: A first unit is configured to collect reflection spectrum signals of a plurality of fiber gratings that have experienced different transmission distances; A second unit is configured to perform frequency domain transformation on the reflection spectrum signals, separate linear phase components and nonlinear phase components, establish a linear mapping relationship between transmission distance and dispersion accumulation, calculate transmission dispersion offset and intrinsic wavelength offset, and obtain initial wavelength offset estimation values of gratings at different spatial positions; A third unit is configured to identify refractive index mutation sections and refractive index stable sections according to the initial wavelength offset estimation values, establish segmented compensation coefficients, calculate correction values through tensor decomposition, correct the initial wavelength offset estimation values, and obtain compensated grating center wavelengths. The fourth unit is configured to construct a wavelength-strain-temperature three-dimensional space, map the compensated grating center wavelength to a track in the three-dimensional space, calculate a track curvature, when the track curvature exceeds a preset curvature threshold, reconstruct a segmented compensation coefficient according to a historical initial wavelength offset estimation value sequence, and output physical quantity demodulation data. The fifth unit is configured to map a wavelength reconstruction residual in the physical quantity demodulation data to the three-dimensional space, and update the segmented compensation coefficient according to a residual gradient direction.

[0079] In a third aspect, the present application provides an electronic device, comprising: a processor; a memory for storing processor-executable instructions; wherein the processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0080] In a fourth aspect, the present application provides a computer-readable storage medium having stored thereon computer program instructions, which, when executed by a processor, implement the method described above.

[0081] The present application can be a method, apparatus, system and / or computer program product. The computer program product can include a computer-readable storage medium having stored thereon computer-readable program instructions that, when executed by a processor, implement various aspects of the present application.

[0082] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

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, identify the refractive index abrupt change segment and the refractive index stable segment, establish segmented compensation coefficients, calculate the correction amount through tensor decomposition, 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 center wavelength of the compensated grating 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 then performed on the phase response curve according to the fitting weights to separate the linear phase component formed by the linear fitting parameters and the nonlinear phase component formed by the residual phase distribution. 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 optic 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