Method for extracting perturbation characteristics and pattern recognition based on sensing signal of microstructured optical fiber

Abstract

This invention provides a disturbance feature extraction and pattern recognition method based on microstructured optical fiber sensing signals, belonging to the field of signal processing technology. By fusing the curvature and torsion statistical ratio features of the polarization-state Bonga sphere trajectory, the fractal dimension and gray-level co-occurrence matrix contrast-dissimilarity difference features of the time-frequency graph, and the correlation coefficient and mutual information features between the polarization and intensity domains, this invention significantly improves the distinguishability and recognition accuracy of different disturbance types, especially exhibiting stronger robustness against easily confused events such as periodic vibrations and sudden impacts. Simultaneously, the classification mechanism based on Mahalanobis distance and adaptive threshold discrimination enables the system to effectively reject unknown disturbance types, avoiding the high false alarm rate problem of traditional closed-set classifiers in open environments. Therefore, in practical applications, it balances recognition accuracy and generalization ability, and the entire feature extraction process does not involve complex iterations, meeting real-time processing requirements.

Description

Technical Field

[0001] This invention belongs to the field of signal processing technology, specifically relating to a method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals. Background Technology

[0002] Distributed fiber optic sensing technology, with its advantages of resistance to electromagnetic interference, corrosion resistance, and long-distance continuous monitoring, has been widely used in perimeter security, structural health monitoring, and oil and gas pipeline safety early warning. Optical time-domain reflectometry (OTDR) based on backscattering Rayleigh scattering and its polarization-sensitive variant (polarized optical time-domain reflectometer, POTDR) senses external disturbances by detecting the evolution of the polarization state of backscattered light in the optical fiber, exhibiting high sensitivity and fast response speed. In recent years, the introduction of microstructured optical fibers (such as photonic crystal fibers and Bragg grating fibers) has further enhanced the fiber's response to weak disturbances. Their unique mode coupling characteristics make the polarization state more sensitive to changes in external parameters such as stress, vibration, and temperature. At the signal processing level, existing methods mainly extract statistical features (such as changes in polarization degree and the elliptic parameters of the Stokes parametric trajectory) from the polarization state sequence or time-frequency domain features (such as the short-time Fourier transform spectrum and wavelet packet energy spectrum) from the light intensity sequence, and use classifiers such as support vector machines and neural networks to identify disturbance events. However, the above feature extraction is often limited to a single physical domain (polarization domain or intensity domain) and does not fully explore the geometric properties of polarization trajectory, thus failing to make full use of the high-dimensional information contained in microstructured optical fibers.

[0003] The main shortcomings of existing technologies are as follows: First, most methods only utilize the instantaneous values ​​of polarization state sequences or their simple statistics (mean, variance), ignoring the local geometric features of the polarization state trajectory on the Bonga sphere, especially the degree of trajectory bending and twisting characterized by curvature and torsion, which are directly physically related to the mechanical modes of external disturbances (such as shearing, torsion, and compression); Second, time-frequency domain feature extraction mostly stays at the level of one-dimensional statistics such as energy spectrum or spectral peaks, lacking systematic modeling of the texture structure of the time-frequency matrix (such as fractal self-similarity and the local gray-level change patterns characterized by the gray-level co-occurrence matrix), resulting in insufficient ability to distinguish between sudden impact disturbances and periodic vibration disturbances; Third, the cross-domain coupling relationship between the polarization domain and the intensity domain is often ignored. In fact, the same external disturbance can simultaneously cause polarization state rotation and light intensity fluctuations, and the linear correlation and nonlinear dependence between the two contain rich information for disturbance type discrimination. Summary of the Invention

[0004] The purpose of this section is to outline some aspects of the embodiments of the present invention and to briefly introduce some preferred embodiments. Some simplifications or omissions may be made in this section, as well as in the abstract and title of the present application, to avoid obscuring the purpose of this section, the abstract and title of the invention. Such simplifications or omissions shall not be used to limit the scope of the present invention.

[0005] In view of the aforementioned existing problems, the present invention is proposed.

[0006] Therefore, the technical problem solved by this invention is: how to achieve accurate identification of known perturbation types and effective rejection of unknown perturbation types.

[0007] To address the aforementioned technical problems, the present invention provides the following technical solution: A perturbation feature extraction and pattern recognition method based on microstructured optical fiber sensing signals includes: acquiring polarization state sequences and light intensity sequences; generating a sequence of trajectory points on a Bonga sphere based on the polarization state sequence, and calculating the discrete curvature sequence and discrete torsion sequence of the trajectory points; performing time-frequency transformation on the light intensity sequence to obtain a time-frequency matrix, calculating the fractal dimension of the time-frequency matrix, and extracting the contrast and dissimilarity of the gray-level co-occurrence matrix of the time-frequency matrix; calculating the first statistical ratio feature and the second statistical ratio feature based on the curvature sequence and the torsion sequence, and calculating the third difference feature based on the contrast and dissimilarity; calculating the correlation coefficient between the curvature sequence and the total energy sequence obtained from the total energy of each frame of the time-frequency matrix, and calculating the mutual information between the short-time energy sequence obtained from the light intensity sequence and the arc length change rate sequence obtained from the polarization state sequence; combining the first statistical ratio feature, the second statistical ratio feature, the third difference feature, the fractal dimension, the contrast, the correlation coefficient, and the mutual information into a feature vector, normalizing it, and inputting it into a classifier to obtain the classification result.

[0008] As a preferred embodiment of the present invention, the calculation of the discrete curvature sequence and discrete torsion sequence of trajectory points includes: mapping each polarization state to trajectory points on a Bonga sphere, with each trajectory point corresponding to an equal time interval; using spherical linear interpolation for the Bonga sphere trajectory between every two adjacent original polarization states, generating a fixed number of interpolated trajectory points at equal arc length intervals to obtain a total interpolated trajectory point sequence; using a fixed number of consecutive points in the total interpolated trajectory point sequence as a sliding window, calculating the discrete curvature and discrete torsion of the center point within the window using the discrete Frenet formula to obtain the discrete curvature sequence and discrete torsion sequence.

[0009] In a preferred embodiment of the present invention, the calculation of the first statistical ratio feature and the second statistical ratio feature includes: performing median filtering on the discrete curvature sequence and the discrete torsion sequence respectively to obtain a denoised curvature sequence and a denoised torsion sequence; calculating the arithmetic mean of the denoised curvature sequence as the curvature mean, calculating the standard deviation of the denoised curvature sequence as the curvature standard deviation, and using the ratio of the curvature mean to the curvature standard deviation as the first statistical ratio feature; calculating the arithmetic mean of the denoised torsion sequence as the torsion mean, calculating the standard deviation of the denoised torsion sequence as the torsion standard deviation, and using the ratio of the torsion mean to the torsion standard deviation as the second statistical ratio feature.

[0010] In a preferred embodiment of the present invention, the calculation of the fractal dimension of the time-frequency matrix and the extraction of the contrast and dissimilarity of the gray-level co-occurrence matrix of the time-frequency matrix include: performing a short-time Fourier transform on the light intensity sequence using a Hanning window and a specified overlap length to obtain a time-frequency matrix, wherein the rows of the time-frequency matrix correspond to frequencies and the columns correspond to time frames; calculating the fractal dimension of the time-frequency matrix using the box counting method to obtain a fractal dimension scalar; quantizing the pixel values ​​of the time-frequency matrix into gray levels of a specified level to generate a gray-level co-occurrence matrix; calculating the contrast scalar of the gray-level co-occurrence matrix and simultaneously calculating the dissimilarity scalar; and using the difference between the contrast scalar and the dissimilarity scalar as a third difference feature.

[0011] As a preferred embodiment of the present invention, the calculation of correlation coefficient and mutual information includes: calculating the sum of squares of all frequency amplitudes in each time frame of the time-frequency matrix to obtain the total energy sequence; interpolating the curvature sequence to make the number of sampling points equal to the length of the total energy sequence to obtain the resampled curvature sequence, and calculating the correlation coefficient between the resampled curvature sequence and the total energy sequence; calculating the short-time energy sequence of the light intensity sequence using a Hanning window and a specified overlap length; calculating the Bonga sphere great circle arc length between adjacent points based on the polarization state sequence to obtain the original arc length change rate sequence; performing linear interpolation on the short-time energy sequence and the original arc length change rate sequence respectively to make them have the same number of sampling points on the same time axis, and discretizing them into a specified number of equally spaced intervals, and calculating the mutual information between them.

[0012] As a preferred embodiment of the present invention, the combination into a feature vector and normalization includes: sequentially concatenating the first statistical ratio feature, the second statistical ratio feature, the third difference feature, the fractal dimension scalar, the contrast scalar, the correlation coefficient, and the mutual information into an original feature vector; reading the global mean vector and the global standard deviation vector of each component from a pre-stored training set statistical file; subtracting the corresponding global mean from each component of the original feature vector and dividing by the corresponding global standard deviation to obtain the normalized feature vector.

[0013] As a preferred embodiment of the present invention, the following steps are taken: The input classifier to obtain the classification result includes: extracting the sample covariance matrix and feature template vector for each perturbation category from the training set, wherein the feature template vector is the mean of the normalized feature vectors of all samples in the corresponding category; adding a fixed regularization coefficient to the covariance matrix of each category and multiplying it by the identity matrix to obtain a regularized covariance matrix; calculating the Mahalanobis distance between the normalized feature vector to be classified and the feature template vector of each category, wherein the covariance matrix of the Mahalanobis distance uses the regularized covariance matrix of the corresponding category; selecting the category with the smallest Mahalanobis distance as the candidate category, and simultaneously calculating the difference between the smallest Mahalanobis distance and the second smallest Mahalanobis distance; if the smallest Mahalanobis distance is less than the mean of the Mahalanobis distances of the corresponding category in the training set plus a specified multiple of the standard deviation, and the difference is greater than a preset difference threshold, then the corresponding candidate category is output; otherwise, an unknown perturbation category is output.

[0014] As a preferred embodiment of the present invention, the spherical linear interpolation generates a fixed number of interpolation trajectory points at equal intervals of circular arc length, and the arc lengths between adjacent interpolation points are equal; the discrete Frenet formula uses a fixed number of sliding windows to calculate the curvature of a discrete point sequence, and the curvature and torsion of the window center point are estimated by the first and second derivatives of the discrete curve, respectively.

[0015] As a preferred embodiment of the present invention, before calculating the correlation coefficient between the resampled curvature sequence and the total energy sequence, the resampled curvature sequence and the total energy sequence are respectively subjected to zero-mean processing; the linear interpolation of the short-time energy sequence and the original arc length change rate sequence is time-axis aligned so that the two have the same number of sampling points on the same time axis, and is discretized into a specified number of intervals in the form of an isofrequency histogram, and the mutual information is calculated using the classical entropy formula.

[0016] In a preferred embodiment of the present invention, the pixel spacing of the gray-level co-occurrence matrix is ​​one and the direction is zero degrees, and the gray-level quantization is a specified number of levels; the contrast scalar is calculated by weighted sum of the squares of the differences between adjacent gray levels, and the dissimilarity scalar is calculated by weighted sum of the absolute values ​​of the differences between adjacent gray levels.

[0017] The beneficial effects of this invention are as follows: Compared with the prior art, the technical effects of this invention are as follows: By fusing the curvature and torsion statistical ratio features of the polarization state Bonga sphere trajectory, the fractal dimension and gray-level co-occurrence matrix contrast-dissimilarity difference features of the time-frequency graph, and the correlation coefficient and mutual information features between the polarization domain and the intensity domain, this invention significantly improves the distinguishability and recognition accuracy of different disturbance types, especially exhibiting stronger robustness to easily confused events such as periodic vibrations and sudden impacts. At the same time, the classification mechanism based on Mahalanobis distance and adaptive threshold discrimination enables the system to effectively reject unknown disturbance types, avoiding the problem of high false alarm rate of traditional closed-set classifiers in open environments. Thus, in practical applications, it balances recognition accuracy and generalization ability, and the entire feature extraction process does not involve complex iterations, meeting the requirements of real-time processing. Attached Figure Description

[0018] Figure 1 This is a flowchart of the disturbance feature extraction and pattern recognition method based on microstructured optical fiber sensing signals described in this invention.

[0019] Figure 2 This is a flowchart of the sub-process for extracting the polarization trajectory curvature and torsion features in this invention.

[0020] Figure 3 This is a flowchart of the time-frequency domain feature and cross-domain interaction feature extraction sub-process in this invention. Detailed Implementation

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

[0022] like Figures 1-3 As shown, the disturbance feature extraction and pattern recognition method based on microstructured optical fiber sensing signals of the present invention includes: S1: Collect the polarization state sequence and intensity sequence of backscattered light from the microstructured optical fiber.

[0023] In some embodiments, step S1 specifically includes: The laser beam enters through the first port of the fiber optic circulator and exits through the second port to the microstructured optical fiber. The microstructured optical fiber is made of photonic crystal fiber or Bragg grating fiber, and its length is set according to the monitoring range, typically ranging from 1 km to 10 km. As the laser beam propagates through the microstructured optical fiber, it generates backscattered Rayleigh light along its path. This backscattered light carries information about the polarization state and intensity of external disturbances (such as vibration, temperature changes, and pressure) experienced by the fiber. The backscattered light returns along the fiber to the second port of the circulator and is then output from the third port to the subsequent detection module.

[0024] The backscattered light output from the third port of the circulator is split into two paths by a 1×2 fiber beam splitter: the first path enters the polarization state measurement module, and the second path enters the light intensity measurement module. The splitting ratio of the beam splitter is 50:50 to maintain the energy balance of the two signals.

[0025] The first backscattered light passes sequentially through a polarization state analyzer consisting of a polarization beamsplitter and four photodetectors. The polarization beamsplitter decomposes the incident light into horizontal and vertical polarization components, which are received by two high-speed photodetectors respectively. Simultaneously, the 45° linear and circular polarization components are measured using a combination of a quarter-wave plate and the polarization beamsplitter. The four detectors synchronously acquire data at the same sampling rate (typically 1MHz to 10MHz, set according to the perturbation frequency range), obtaining the time-series signals of four Stokes parameters (S0, S1, S2, S3). For each sampling time i, normalization is performed to obtain the Bonga spherical coordinates.

[0026] The coordinates of the Bonga sphere at all sampling times are arranged in chronological order to form a polarization state sequence. The observation time window is set to 1ms to 1s based on the expected duration of the disturbance event.

[0027] The second backscattered light is directly incident on a high-speed photodetector. This detector, along with the four detectors in the polarization state measurement module, is triggered by the same clock source to ensure strict time synchronization. The detector outputs an analog voltage signal, which is converted into a digital signal by an analog-to-digital converter at the same sampling rate, yielding a light intensity sequence. The light intensity sequence and the polarization state sequence correspond one-to-one in time, with a time deviation of less than one sampling period.

[0028] The acquired polarization state and light intensity sequences are stored as discrete digital signals in a buffer or memory. Each record contains a timestamp, a normalized Stokes vector, and a light intensity value. The sequence length is... This corresponds to a complete observation time window.

[0029] In the technical solution of this disclosure, low-noise, high-sensitivity acquisition of perturbation signals along the fiber is achieved by injecting narrow-linewidth continuous laser into a microstructured optical fiber and using a circulator to separate the backscattered light. The polarization state measurement module adopts a four-channel synchronous detection method, which can completely acquire the polarization evolution trajectory on the Bonga sphere, providing a physical basis for subsequent curvature and torsion analysis.

[0030] S2: Generate a sequence of trajectory points on the Bonga sphere based on the polarization state sequence, and calculate the discrete curvature sequence and discrete torsion sequence of the trajectory points.

[0031] In some embodiments, step S2 specifically includes the following sub-steps: S2.1: Map the polarization state sequence to trajectory points on a Bonga sphere, with each point corresponding to an equal time interval.

[0032] Read the polarization state sequence stored in step S1 ,in, And they satisfy the condition that their sum is 1. Each It is directly considered as a point on a 3D Bonga sphere. Since the sampling interval is constant, each trajectory point corresponds to an equal time interval, and the time axis is aligned. The sampling interval is... ,in, The sampling frequency.

[0033] S2.2: Spherical linear interpolation is used for the Bonga sphere trajectory between every two adjacent original polarization states to generate a fixed number of interpolated trajectory points at equal arc length intervals.

[0034] For any adjacent original polarization states and Calculate the length of the great circle arc on the Bonga sphere: ; Interpolation points are generated using the spherical linear interpolation formula. Let the interpolation parameters... The coordinates of the interpolation point are: ; in, For a sine function, to ensure that the arc lengths between adjacent interpolation points are equal, the interval is... Divided into equal parts Section, take , ,get Interpolation points (including endpoints). correspond , correspond ,middle Interpolation points. For the sake of brevity, in this embodiment of the invention, each interval is referred to as generating [interpolation points]. An interpolation sub-segment or a new one Each point, used in... As the standard.

[0035] The arc length between adjacent interpolation points is .

[0036] right arrive Repeat the above interpolation, and concatenate the beginning and end of each subsequence (adjacent intervals share endpoints). ), thus obtaining the total interpolation trajectory point sequence ,in, .

[0037] Optionally, when When the value is close to 0 (i.e., adjacent polarization states are almost identical), to avoid numerical instability, linear interpolation is directly used instead of spherical linear interpolation, and the corresponding processing logic is given in the specification.

[0038] S2.3: Using a fixed number of consecutive points in the total interpolation trajectory point sequence as a sliding window, the discrete curvature and discrete torsion of the center point within the window are calculated using the discrete Frenet formula.

[0039] Set the length of the sliding window (Odd number), the index of the point within the window is ,in The center point of the window. From 3 to Traverse and calculate the discrete curvature at the center point. and discrete torsion .

[0040] Discrete curvature is calculated using the three-point method, utilizing points : ; in, denoted by , represents the three-dimensional Euclidean distance, and × represents the cross product of vectors. Geometrically, the above formula represents the reciprocal of the circumradius determined by the three points.

[0041] Discrete torsion requires four points The calculation formula is: ; in, The torsion is a scalar triple product of three vectors. The numerator represents the directed volume of the parallelepiped spanned by the three vectors, and the denominator is the square of the curvature-related quantity. The sign of the torsion indicates the direction in which the curve deviates from the plane.

[0042] It should be noted that the above torsion formula requires discrete points to be parameterized with approximately equal arc lengths. This step generates interpolation points with equal arc length intervals through spherical linear interpolation, satisfying this prerequisite. If the denominator approaches zero due to numerical errors, the torsion of that point is set to zero and marked as an outlier, which can be removed in subsequent median filtering.

[0043] For each valid center point calculate and After traversal, a discrete curvature sequence is obtained. and discrete torsion sequence Both sequences have a length of 1. .

[0044] Finally, the calculated curvature and torsion sequences are used as the output of step S2 and fed into subsequent steps for median filtering and feature extraction.

[0045] It should be noted that in the technical solution of this disclosure, spherical linear interpolation is used to densify the sparse original polarization state sequence, generating interpolation trajectory points with equal arc length intervals, thus ensuring the numerical stability of subsequent curvature and torsion calculations. Compared with linear interpolation, spherical linear interpolation strictly constrains the interpolation points to the great arc of the Bonga sphere, conforming to the physical laws of polarization state evolution and avoiding geometric distortion introduced by the interpolation points deviating from the spherical surface. Furthermore, by employing a five-point sliding window combined with the discrete Frenet formula, the curvature and torsion of the curve can be estimated simultaneously within a local range, exhibiting better noise resistance than a three-point window.

[0046] S3: Perform time-frequency transformation on the light intensity sequence to obtain the time-frequency matrix, calculate the fractal dimension of the time-frequency matrix, and extract the contrast and dissimilarity of the gray-level co-occurrence matrix of the time-frequency matrix.

[0047] S3.1: Perform a short-time Fourier transform on the light intensity sequence using a Hanning window and a specified overlap length to obtain the time-frequency matrix.

[0048] Read the light intensity sequence acquired in step S1, and use Short Time Fourier Transform (STFT) to convert the one-dimensional time-domain signal into a two-dimensional time-frequency representation: Set the window length of the Hanning window to be There are sampling points, and the number of overlapping points between adjacent windows is . That is, the overlap rate is 50%. The actual time span corresponding to the window length is... Typical values ​​range from 0.256ms to 2.56ms (depending on...). The Hanning window function is defined as follows: ; For the first Extract the light intensity sequence starting from the index of each time frame. The beginning Let there be 1 point, denoted as 1. After windowing, perform a Discrete Fourier Transform (DFT): ; in, For frequency index, For time frame indexing, The total number of frames is calculated using the following formula: ; Because the Fourier transform of a real signal has conjugate symmetry, it typically only retains the non-negative frequency components, i.e., it takes... (like (The number is even). Therefore, the time-frequency matrix The dimension is ,in, The matrix represents the number of frequency points, with each column corresponding to the frequency distribution of a time frame, and each row corresponding to the change of a fixed frequency over time. Matrix elements. That is, the amplitude spectrum (or energy spectrum, taken as square). In this embodiment, the amplitude spectrum is used.

[0049] It should be noted that converting the light intensity sequence into a time-frequency matrix using short-time Fourier transform preserves the time-varying characteristics of the signal and reveals the evolution of the frequency distribution over time, providing rich two-dimensional information for subsequent texture feature extraction. Furthermore, the use of the Hanning window effectively suppresses spectral leakage, and the overlapping window design ensures a balance between time and frequency resolution.

[0050] S3.2: The fractal dimension of the time-frequency matrix is ​​calculated using the box counting method to obtain the fractal dimension scalar.

[0051] time-frequency matrix Consider it as a two-dimensional grayscale image with dimensions of Line (frequency) × Column (time). The fractal dimension is calculated using box counting and is used to characterize the texture self-similarity and complexity of the matrix. The specific steps are as follows:

[0052] (1) Set the matrix elements Normalization to The interval is used to obtain the normalized matrix. .

[0053] (2) Set the range of values ​​for the box size S, usually taking the range of values ​​for S. ,in, No more than Half of it. For a typical time-frequency matrix (such as... ), where S can be 2, 4, 8, 16, 32, 64.

[0054] (3) Normalize the matrix Consider it as a surface in three-dimensional space, where Coordinates represent planar positions and grayscale values. Let the height be . For a side length of , Given a square box, find the maximum gray value within the (i,j)th grid cell. and minimum value The number of boxes required for this grid for: ; in, The maximum height after degree normalization (usually taken as 1). The sum of the number of boxes in all grids is... Insufficient at the boundary The grid is processed according to the actual size.

[0055] It should be noted that, since the grayscale value range of the time-frequency matrix (normalized to 0-1) differs from the spatial dimension (number of pixels), the scale in the height direction of the box is uniformly represented by S (pixel unit) for both the physical and spatial dimensions. This is an empirical simplification. If the dynamic range of the time-frequency matrix is ​​too small, the grayscale can be linearly stretched to [0, ... ],in Can be taken as .

[0056] (4) Change the box size S and repeat step 3 to obtain a series of point pairs. .

[0057] (5) Fit the line using the least squares method: Among them, slope That is, the fractal dimension scalar. , This is the intercept of the fitted line.

[0058] Among them, fractal dimension scalar The value of fractal dimension typically ranges from 2 to 3 (the theoretical maximum value is 3 for two-dimensional images). This value reflects the roughness of the texture of the time-frequency plot: the more complex the disturbance signal and the richer the frequency components, the larger the fractal dimension.

[0059] It should be noted that fractal dimension, as a measure of the overall complexity of time-frequency graphs, can effectively distinguish between periodic vibrations (low fractal dimension, close to 2) and random shock disturbances (high fractal dimension, close to 3). The embodiment of this invention uses box counting method for calculation, which is simple to implement and robust to scale changes.

[0060] S3.3: Quantize the pixel values ​​of the time-frequency matrix into gray levels of a specified number to generate a gray-level co-occurrence matrix.

[0061] First, the time-frequency matrix Each element is quantized to grayscale levels. Let the number of quantization levels be . (Typical value range: 8–32). Calculate the maximum value of the matrix. and minimum value Map each element to Integer gray levels: ; Obtain the gray quantization matrix Its size is Elements are integers from 0 to 15.

[0062] Then, a gray-level co-occurrence matrix is ​​constructed to describe the spatial dependencies between pixel gray levels. This embodiment uses a pixel spacing (step size) of 1 and a direction of 0° (horizontal direction, i.e., along the time axis). For each gray level pair... Statistics in the matrix The number of adjacent pixel pairs that satisfy the following condition: the gray level of the left pixel is The right-hand pixel (shifted 1 pixel horizontally to the right) has a grayscale level of [missing value]. Border handling: Pixels in the last column that have no right neighbor are not included in the calculation.

[0063] After the statistics are completed, the counting matrix is ​​divided by the total number of adjacent pixel pairs (i.e., ), to obtain the normalized gray-level co-occurrence matrix , dimension Each element Represents grayscale level and The probability of them occurring simultaneously under a given spatial relationship.

[0064] S3.4: Calculate the contrast scalar and dissimilarity scalar of the gray-level co-occurrence matrix, and calculate the difference between the two as the third difference feature.

[0065] Specifically, based on the normalized gray-level co-occurrence matrix Define the following two texture features: Contrast ratio reflects the degree of drastic change in local grayscale values, and its calculation formula is as follows: ;

[0066] The larger the value, the greater the difference in grayscale between adjacent pixels in the image, and the clearer the texture.

[0067] Dissimilarity is similar to contrast, but uses absolute values ​​instead of squares, making it less sensitive to noise. The calculation formula is: ; Dissimilarity also reflects the degree of grayscale difference, but the weight increases linearly, while the weight of contrast increases quadratically.

[0068] The difference between contrast and dissimilarity is used as the third difference feature. This difference feature can highlight areas in the texture with strong and uneven gray-level changes, and has higher sensitivity to certain types of disturbances (such as the breakage of time-frequency pattern stripes caused by sudden impacts).

[0069] Finally, the fractal dimension calculated above is... Contrast, dissimilarity, and third difference features Package and output.

[0070] It should be noted that, in this embodiment of the invention, the contrast and dissimilarity of the gray-level co-occurrence matrix characterize the gray-level variation features of the time-frequency map from the perspective of local texture: contrast is sensitive to abrupt changes, while dissimilarity is sensitive to slow changes, and the difference between the two can further amplify the differences between the two types of perturbations. By combining fractal dimension with GLCM features, a multi-scale time-frequency domain feature set from global to local is formed, thereby making up for the inadequacy of the expressive power of a single feature.

[0071] S4: Calculate the first and second statistical ratio features based on the curvature and torsion sequences, and calculate the third difference feature based on contrast and dissimilarity.

[0072] S4.1: Perform median filtering on the discrete curvature sequence and the discrete torsion sequence respectively to obtain the denoised curvature sequence and the denoised torsion sequence.

[0073] Read the output discrete curvature and discrete torsion sequences. Since noise may be introduced during the original curvature and torsion calculation process (such as shot noise from the photodetector and environmental vibration interference), median filtering is used for smoothing.

[0074] Set the filter window length to (Odd number). For the curvature sequence, for each index Take three values ​​within the window (from the second element to the second-to-last element of the sequence). After sorting them, the median is taken as the filtered value. Boundary points (the first and last elements of the sequence) are either preserved in their original value or their neighboring values ​​are copied. The same operation is performed on the torsion sequence to obtain the denoised torsion sequence. and The length remains unchanged.

[0075] It should be noted that median filtering (with a window length of 3) effectively suppresses the influence of isolated noise points on statistical features when denoising discrete curvature and torsion sequences, while preserving the local trend of the sequence. Compared with mean filtering, median filtering is better at preserving edge information and is suitable for handling transient jumps that may exist in curvature sequences.

[0076] S4.2: Calculate the arithmetic mean of the denoised curvature sequence as the curvature mean, and calculate the standard deviation of the denoised curvature sequence as the curvature standard deviation.

[0077] The arithmetic mean (first sample moment) of the denoised curvature sequence is defined as: ; Furthermore, the standard deviation of curvature (the square root of the second central moment) is defined as: ; The summation iterates through all... (common (4 items), denominator adopts (Unbiased estimation). This reflects the average level of curvature of the polarization trajectory on the Bonga sphere. This reflects the intensity of curvature fluctuations.

[0078] S4.3: The ratio of the mean curvature to the standard deviation of curvature is used as the first statistical ratio feature.

[0079] Among them, the first statistical ratio characteristic for: ; This ratio characterizes the relative dispersion of the curvature sequence. When the curvature sequence has small fluctuations (small standard deviation) and a large mean, A larger value indicates that the trajectory exhibits a stable, large curvature (such as continuous oscillation); when the curvature sequence fluctuates wildly (large standard deviation), A smaller value indicates that the trajectory contains abrupt changes or noise. This ratio is dimensionless and has a certain degree of robustness to changes in signal amplitude.

[0080] S4.4: Calculate the arithmetic mean of the denoised torsion sequence as the torsion mean, and calculate the standard deviation of the denoised torsion sequence as the torsion standard deviation.

[0081] Similarly, the arithmetic mean of the denoised torsion sequence is: ; The standard deviation of the torsion is: ; in, This reflects the average degree of torsion of the polarization trajectory deviating from the plane curve. This reflects the volatility of the torsion.

[0082] S4.5: Use the ratio of the mean to the standard deviation of torsion as the second statistical ratio feature.

[0083] Among them, the second statistical ratio characteristic for: ; and similar, Characterizes the relative stability of the torsion sequence. For periodic torsion (such as rotational perturbations) Non-zero and Smaller Larger; for trajectories dominated by random noise Approaching zero and Larger Smaller.

[0084] It should be noted that using the ratio of the mean to the standard deviation as a statistical feature, rather than directly using the mean or standard deviation, has the following advantages: the ratio is dimensionless, eliminating the amplitude influence caused by different sampling rates or trajectory scales; this ratio can simultaneously reflect the central tendency and dispersion of the sequence, containing more information than a single statistic; for certain types of perturbations (such as periodic vibrations and random shocks), the ratio feature has significantly higher discriminative power than the mean or standard deviation alone. Experiments show that the curvature ratio feature is sensitive to polarization state changes caused by fiber bending, while the torsion ratio feature is sensitive to torsional perturbations; the two are complementary.

[0085] S4.6: Calculate the third difference feature based on contrast and dissimilarity.

[0086] Read the contrast scalar and dissimilarity scalar of the gray-level co-occurrence matrix calculated in step S3. Define the third difference feature. This is the difference between the two. It should be noted that this difference feature can highlight areas with strong and uneven gray-level jumps in the time-frequency image. Contrast is weighted by squares, assigning higher weights to larger gray-level differences; dissimilarity is weighted by absolute values, with the weights increasing linearly. After subtracting the two, a larger difference indicates the presence of a small number of but extremely large gray-level jumps in the texture (such as the breakage of stripes in the time-frequency image caused by a sudden impact); a difference close to zero or a negative value indicates that the gray-level changes in the texture are relatively gentle or uniform.

[0087] Furthermore, the third difference feature utilizes the difference between contrast and dissimilarity to further amplify the non-uniform texture features in the time-frequency plot. Since contrast is more sensitive to large gray-level differences, while dissimilarity treats all gray-level differences equally, the difference between the two can reveal whether there are a few but significant gray-level abrupt changes. This is of unique value for identifying sudden disturbances (such as transient changes caused by external impacts or local overheating).

[0088] S5: Calculate the correlation coefficient between the curvature sequence and the total energy sequence obtained from the total energy of each frame of the time-frequency matrix, and calculate the mutual information between the short-time energy sequence obtained from the light intensity sequence and the arc length change rate sequence obtained from the polarization state sequence.

[0089] In some embodiments, step S5 specifically includes the following sub-steps: S5.1: Calculate the sum of squares of all frequency amplitudes in each time frame of the time-frequency matrix to obtain the total energy sequence.

[0090] For each time frame Calculate the sum of squared amplitudes at all frequency points within the frame, which is defined as the total energy: ; Obtain the total energy sequence , length is This sequence reflects the energy fluctuations of the optical signal in the time-frequency domain over time: when a disturbance occurs, the energy is redistributed within a specific frequency band, and the total energy sequence shows corresponding fluctuations.

[0091] S5.2: Interpolate the curvature sequence to make the number of sampling points equal to the length of the total energy sequence, thus obtaining a resampled curvature sequence.

[0092] Since the sampling points of the curvature sequence correspond to the positions of equal arc lengths of the Bonga sphere trajectory points, while the sampling points of the total energy sequence correspond to equal time frames, their sampling mechanisms are different, and the number of sampling points is generally not equal. Usually much larger To calculate the correlation between the two, the curvature sequence needs to be interpolated to the same time axis as the total energy sequence.

[0093] First, establish the time coordinates of the curvature sequence: the first... Each point corresponds to an index in the total interpolation trajectory point sequence. (Because curvature calculation uses a five-point sliding window, with the window center point being the first...) (One interpolation point). This interpolation point is located on a spherical linear interpolation arc segment between two original polarization states, corresponding to the time coordinate. for: ; in, That is, the index of the original sampling interval to which the interpolation point belongs; That is, the sub-segment number of the interpolation point within the current interval ( (corresponding to the starting point of the interval); the above calculations ensure that the time coordinates of each curvature point are accurate to the interpolation sub-interval level. Based on this, a time coordinate array of the curvature sequence can be constructed. .

[0094] Since the number of sampling points in the curvature sequence is usually much larger than the length of the total energy sequence, to align them to the same time axis, cubic spline interpolation is first performed on the curvature sequence to make its time coordinates match those of the total energy sequence. Completely identical. The length obtained after interpolation is... resampled curvature sequence If non-monotonic or overshoot phenomena occur during the interpolation process, piecewise linear interpolation can be used instead (this embodiment of the invention will not be elaborated on, and the technical solution can be selected according to actual needs).

[0095] Cubic spline interpolation can ensure the continuity of the first and second derivatives of the interpolation curve, and is suitable for smooth curvature changes.

[0096] S5.3: Perform zero-mean processing on the resampled curvature sequence and the total energy sequence respectively, and calculate the Pearson correlation coefficient between the two.

[0097] Specifically, to eliminate the influence of the sequence mean on the correlation coefficient, the two sequences are first zero-mean normalized: ; ; Then calculate the Pearson correlation coefficient. : ; Correlation coefficient The range of values ​​is Positive values ​​indicate a positive correlation between curvature change and total energy change (e.g., a disturbance simultaneously increases polarization curvature and signal energy), while negative values ​​indicate a negative correlation. This coefficient, as the first component of the cross-domain interaction characteristic, is denoted as... .

[0098] S5.4: Calculate short-time energy sequences for light intensity sequences using the Hanning window and a specified overlap length.

[0099] The original light intensity sequence is reread using the same window parameters as in step S3.1: a Hanning window length of 256 and an overlap number of 128. For each time frame... Extract light intensity sampling points within the window Calculate the short-time energy of this frame: ; in, Using the Hanning window function, a short-time energy sequence is obtained. The length is also The time axis is completely consistent with the total energy sequence. Short-time energy reflects the instantaneous intensity fluctuation of the optical signal in the time domain and is sensitive to amplitude changes caused by disturbances.

[0100] S5.5: Calculate the Bonga sphere great circle arc length between adjacent points based on the polarization state sequence to obtain the original arc length change rate sequence.

[0101] First, read the original polarization state sequence from step S1.3 (before interpolation). For each pair of adjacent points, calculate the great circle arc length of both on the Bonga sphere: ; Because the original polarization state sequence corresponds to equal time intervals Therefore, the rate of change of arc length is defined as: ;

[0102] Obtain the original arc length change rate sequence , length is The time axis is aligned with the sampling time of the light intensity sequence (each Corresponding time interval The average rate of change. This sequence characterizes the angular velocity of the polarization state rotating on the Bonga sphere, reflecting the instantaneous change in the modulation rate of the polarization state due to external disturbances.

[0103] S5.6: Perform linear interpolation on the short-time energy sequence and the original arc length change rate sequence respectively, so that they have the same number of sampling points on the same time axis, and discretize them into a specified number of equally spaced intervals, and calculate their mutual information.

[0104] Because the short-time energy sequence length is (Number of time frames), length of the original arc length change rate sequence is (Number of sampling points, usually) Much larger Both need to be interpolated to the same time axis to achieve alignment. The time axis of the short-time energy series is chosen as the reference, where... Defined as the first The center moment of each time frame: ; The rate of change of the original arc length The midpoint of the corresponding time interval is: ; The two may have a fixed offset on the timeline. : ; To ensure physical synchronization, we must address Perform correction: ; Linear interpolation is used to calculate in Value at: ; in, (Boundary) or At times, only one-sided linear extrapolation is used, or one frame at the beginning and one at the end are discarded. After interpolation, a sequence of the same length as the short-time energy sequence is obtained. .

[0105] Thus, we have obtained two sequences of equal length: and .

[0106] Furthermore, to calculate mutual information, the two continuous sequences need to be discretized. The number of discrete intervals needs to be determined. Calculate the value ranges of the two sequences respectively. and and divide it into equal parts An equal-frequency histogram (equal-interval method is also acceptable). For each time frame ,Sure The interval number that falls into ,as well as The interval number that falls into Joint statistical distribution frequency: ; Its edge distribution is as follows: ; Mutual information Defined as: ; The base of the logarithm is usually set to 2, in which case the unit of mutual information is bits. This value quantifies the statistical dependence between two sequences: the greater the mutual information, the more information is shared between the short-time energy sequence and the arc length change rate sequence, meaning that there is a strong correlation between the instantaneous fluctuations in light intensity and the polarization state rotation rate.

[0107] Finally, the calculated correlation coefficients and mutual information as output , recorded as .

[0108] It should be noted that the technical solution of this disclosure embodiment is designed with two cross-domain interaction features, which characterize the physical coupling relationship between the polarization domain and the time-frequency domain from the two levels of linear correlation and nonlinear dependence, respectively.

[0109] The first characteristic is the Pearson correlation coefficient between the curvature sequence and the total energy sequence. The curvature sequence reflects the degree of bending of the polarization state on the Bonga sphere, while the total energy sequence reflects the intensity distribution of the signal in the time-frequency domain. For the same external disturbance (such as vibration or shock), changes in polarization state and energy often exhibit a certain degree of synchronicity: as the disturbance intensifies, the curvature of the polarization trajectory increases, and simultaneously, the signal energy diffuses towards higher frequencies, causing fluctuations in the total energy. The correlation coefficient quantifies the degree of this linear correlation, providing a basis for distinguishing different types of disturbances. For example, low-frequency periodic vibrations may produce a strong positive correlation, while sudden shocks may produce a weak negative correlation. Zero-mean processing eliminates the influence of the DC component, ensuring that the correlation coefficient only reflects the consistency of the fluctuation pattern.

[0110] The second feature is the mutual information between the short-time energy sequence and the arc length change rate sequence. Short-time energy describes the instantaneous intensity of light from a time-domain perspective, while the arc length change rate describes the angular velocity of polarization state rotation. Although their physical meanings differ, both are driven by the same perturbation source. Mutual information, as a nonlinear measure, can capture the nonlinear statistical dependence between the two, overcoming the limitations of linear correlation coefficients. For example, when the perturbation signal exhibits nonlinear modulation or saturation effects, mutual information can still effectively reflect its correlation. Interpolating the two sequences to the same time axis and discretizing them into equally spaced intervals ensures computational consistency and numerical stability.

[0111] The two features mentioned above together construct a cross-domain interactive feature set, which complements the aforementioned polarization domain statistical ratio feature and time-frequency domain texture feature. This allows the final feature vector to comprehensively characterize the fiber optic disturbance signal from multiple dimensions, significantly improving the classifier's discriminative ability. The feature extraction process does not involve complex iterations, meeting the real-time processing requirements of automotive or embedded systems.

[0112] S6: Combine the first statistical ratio feature, the second statistical ratio feature, the third difference feature, the fractal dimension, the contrast, the correlation coefficient, and the mutual information into a feature vector, normalize it, and input it into the classifier to obtain the classification result.

[0113] S6.1: Concatenate them into a single original feature vector.

[0114] These seven features are concatenated in a fixed order to form a seven-dimensional original feature vector: ; The superscript T indicates transpose. The components of this vector have different dimensions and ranges of values. Directly using the original values ​​for classification will lead to features with larger dimensions dominating the distance metric, so normalization is required.

[0115] S6.2: During the offline training phase, collect at least 1000 labeled samples. Each sample contains known perturbation types (such as no perturbation, periodic vibration, impact, temperature drift, etc.) or is explicitly marked as unknown (for training the rejection threshold). Extract seven-dimensional features for each sample, calculate the global mean and global standard deviation of each component of all training samples (excluding the unknown class), and store them as a JSON format file. Subtract the corresponding global mean from each component of the original feature vector and divide by the corresponding global standard deviation to obtain the normalized feature vector.

[0116] Specifically, during the offline training phase, a large number of labeled samples (containing various perturbation events, such as undisturbed events, periodic vibrations, shocks, temperature drifts, etc.) are used to extract the aforementioned seven-dimensional features, and the global mean of each feature component is calculated across all training samples. and global standard deviation These statistics are stored in a configuration file (such as a JSON or binary file).

[0117] For the current sample to be identified, read the pre-stored global mean vector. and global standard deviation vector Z-score normalization is applied to the original feature vectors: ; in, The first feature vector is the first feature vector. Each component (unnormalized eigenvalue); For the standardized first Each component.

[0118] Obtain the normalized eigenvectors After standardization, the mean of each component is 0 and the standard deviation is 1, which eliminates the differences in dimensions and scales, and makes each feature contribute equally to the classifier.

[0119] S6.3: Extract the sample covariance matrix and feature template vector for each perturbation category from the training set.

[0120] Assume there is a total There are 1 known disturbance categories, denoted as _____. During the offline training phase, for each category... Collect the normalized feature vectors of all training samples of this class, and calculate the mean vector as the feature template for this class: ; in, For category The training sample set, This represents the number of samples in this class. Simultaneously, calculate the sample covariance matrix for this class: ; in, It is a 7×7 symmetric positive definite matrix that characterizes the internal scattering structure of the eigenvectors of this class. The template vector and covariance matrix are pre-computed and stored for online classification.

[0121] Known disturbance category The criteria need to be defined according to the actual application scenario, such as undisturbed vibration, periodic vibration (frequency range 10-200Hz), impact disturbance (duration less than 10ms), temperature drift (rate of change less than 0.1°C / s), etc. The number of training samples for each category should be no less than 30 to ensure the reliability of the covariance matrix estimation.

[0122] S6.4: Add a fixed regularization coefficient to the covariance matrix of each category and multiply it by the identity matrix to obtain the regularized covariance matrix.

[0123] It should be noted that, due to the limited number of training samples or the linear correlation between features, the covariance matrix... The covariance matrix may be singular or ill-conditioned, causing Mahalanobis distance calculation to fail. To address this issue, a ridge regression regularization method is used, adding a diagonal perturbation to each covariance matrix: .in, It is a 7×7 identity matrix. is the regularization coefficient. The value of is determined on the training set through five-fold cross-validation. The value that maximizes the average classification accuracy across all categories is selected from the range provided. In this embodiment, the optimal value on the training set is [value to be specified in this embodiment]. After regularization It is strictly positive definite and reversible.

[0124] S6.5: Calculate the Mahalanobis distance between the normalized feature vector to be classified and the feature template vector of each class. The covariance matrix of the Mahalanobis distance is the regularized covariance matrix of the corresponding class.

[0125] Specifically, for the normalized feature vector of the sample to be classified Calculate its value to each category separately. Mahalanobis distance : ;

[0126] It should be noted that Mahalanobis distance takes into account the variance of each feature component and the correlation between features, and has a higher discriminative power than Euclidean distance. The smaller the distance, the closer the sample is to the template class.

[0127] Furthermore, using Mahalanobis distance as a similarity measure has the following core advantages: it automatically weights the variances of each feature through the covariance matrix, automatically assigning smaller weights to features with larger variances; it considers the correlation between features, avoiding the distortion of distance by redundant information; and Mahalanobis distance is dimensionless and independent of feature scale. By estimating the covariance matrix separately for each category, the classifier can adapt to the internal distribution differences of different disturbance types (e.g., the feature variance of the periodic vibration category is smaller, while the variance of the shock disturbance category is larger).

[0128] Furthermore, we first find the minimum Mahalanobis distance. and their corresponding categories Then find the second smallest Mahalanobis distance. And calculate the difference between the two. This difference reflects the classifier's confidence in the current sample's classification: the larger the difference, the more significantly the candidate class is superior to other classes, and the more reliable the classification result.

[0129] S6.6: Determine the classification result or unknown perturbation category based on the threshold.

[0130] Specifically, during the offline training phase, for each category Calculate the mean Mahalanobis distance of all training samples in this class to their own templates. and standard deviation .

[0131] Specifically, for each training sample ,calculate: ;

[0132] Then calculate the mean of these distances. and standard deviation These statistics are also stored in advance.

[0133] For a sample to be classified, if both of the following conditions are met, a candidate class is output; otherwise, an unknown perturbation class is output: Condition 1: + ,in, For a specified multiple, joint optimization is performed using cross-validation to maximize the F1 score on the validation set. This embodiment takes... =2.0. This condition ensures that the distance between a sample and a candidate class does not exceed the normal fluctuation range of the distance between the sample and the training samples of that class. For the training phase, for the category The arithmetic mean of the Mahalanobis distances of all training samples to their own class templates reflects the central tendency of the distances within that class of samples. For the training phase, for the category The sample standard deviation of the Mahalanobis distance from all training samples to their own class template reflects the dispersion of the distances within that class.

[0134] Condition 2: If the difference exceeds a preset threshold, joint optimization is performed using cross-validation to maximize the F1 score on the validation set; in this embodiment, this threshold is 0.8. This condition ensures that the candidate class is significantly better than other classes, avoiding ambiguous decisions. If both conditions are met, the candidate class is output. As the classification result of the perturbation event; otherwise, output the unknown perturbation category (indicating that the sample does not belong to any known perturbation type, or the classification confidence is insufficient).

[0135] Finally, the classification results (known category labels or unknown identifiers) are output in a standard format for use by upper-layer applications (such as alarm systems, data loggers, or operation and maintenance platforms). Optional auxiliary information such as Mahalanobis distance values ​​and differences can also be output for debugging or confidence analysis.

[0136] It should be noted that the two conditions in the threshold determination mechanism control the output quality from the perspectives of absolute distance and relative discriminative power: Condition 1 prevents abnormal samples that are far from all known categories from being misclassified as a known class; Condition 2 prevents arbitrary selection between two close categories. When either condition is not met, the unknown disturbance category is output, giving the system the ability to reject such disturbances, making it suitable for disturbance monitoring tasks in open environments.

[0137] The method also includes one or more processors and memory.

[0138] The memory is used to store operable instructions that, when executed by the one or more processors, cause the one or more processors to perform operations, including the flow of the disturbance feature extraction and pattern recognition method based on microstructured fiber optic sensing signals described in the foregoing embodiments, especially... Figure 1 The flowchart of the method is shown.

[0139] Other aspects disclosed in the embodiments of the present invention also propose a computer-readable medium for storing software including instructions executable by one or more computers, which, upon execution, cause the one or more computers to perform operations including the flow of the disturbance feature extraction and pattern recognition method based on microstructured optical fiber sensing signals of the foregoing embodiments, particularly... Figure 1 The flowchart of the method is shown.

[0140] It should be recognized that embodiments of the present invention may be implemented or carried out by computer hardware, a combination of hardware and software, or by computer instructions stored in a non-transitory computer-readable storage medium.

[0141] The method can be implemented using standard programming techniques, including a non-transitory computer-readable storage medium configured with a computer program in the computer program, wherein the storage medium is configured such that the computer operates in a specific and predefined manner.

[0142] Each program can be implemented in a high-level procedural or object-oriented programming language to communicate with the computer system; however, if required, the program can be implemented in assembly or machine language.

[0143] In any case, the language can be either compiled or interpreted.

[0144] Furthermore, for this purpose, the program can run on programmed application-specific integrated circuits.

[0145] The processes described herein (or variations and / or combinations thereof) can be executed under the control of one or more computer systems configured with executable instructions, and can be implemented by hardware or a combination thereof as code (e.g., executable instructions, one or more computer programs, or one or more applications) that commonly executes on one or more processors. The computer program includes a plurality of instructions executable by one or more processors.

[0146] Furthermore, the method can be implemented in any suitable computing platform, including but not limited to personal computers, minicomputers, mainframes, workstations, networked or distributed computing environments, standalone or integrated computer platforms, or in communication with charged particle tools or other imaging devices.

[0147] Various aspects of the present invention can be implemented in machine-readable code stored on a non-transitory storage medium or device, whether portable or integrated into a computing platform, such as a hard disk, optical read and / or write storage medium, RAM, ROM, etc., such that it can be read by a programmable computer, and when the storage medium or device is read by the computer, it can be used to configure and operate the computer to perform the processes described herein.

[0148] Furthermore, machine-readable code, or parts thereof, can be transmitted via wired or wireless networks.

[0149] When such media includes instructions or programs that combine with a microprocessor or other data processor to implement the steps described above, the invention described herein includes these and other different types of non-transitory computer-readable storage media.

[0150] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals, characterized in that, include: Obtain the polarization state sequence and light intensity sequence; Generate a sequence of trajectory points on the Bonga sphere based on the polarization state sequence, and calculate the discrete curvature sequence and discrete torsion sequence of the trajectory points; The time-frequency matrix is ​​obtained by performing time-frequency transformation on the light intensity sequence. The fractal dimension of the time-frequency matrix is ​​calculated, and the contrast and dissimilarity of the gray-level co-occurrence matrix of the time-frequency matrix are extracted. The first and second statistical ratio features are calculated based on curvature and torsion sequences, and the third difference feature is calculated based on contrast and dissimilarity. Calculate the correlation coefficient between the curvature sequence and the total energy sequence obtained from the total energy of each frame of the time-frequency matrix, and calculate the mutual information between the short-time energy sequence obtained from the light intensity sequence and the arc length change rate sequence obtained from the polarization state sequence; The first statistical ratio feature, the second statistical ratio feature, the third difference feature, the fractal dimension, the contrast, the correlation coefficient, and the mutual information are combined into a feature vector, which is then normalized and input into the classifier to obtain the classification result.

2. The method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals according to claim 1, characterized in that, The discrete curvature sequence and discrete torsion sequence of the calculated trajectory points include: Each polarization state is mapped to a trajectory point on a Bonga sphere, with each trajectory point corresponding to an equal time interval; Spherical linear interpolation is used for the Bonga sphere trajectory between every two adjacent original polarization states. A fixed number of interpolated trajectory points are generated at intervals of equal-sized circular arc lengths to obtain the total interpolated trajectory point sequence. Using a fixed number of points in the total interpolation trajectory point sequence as a sliding window, the discrete curvature and discrete torsion of the center point within the window are calculated using the discrete Frenet formula, resulting in a discrete curvature sequence and a discrete torsion sequence.

3. The method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals according to claim 2, characterized in that: The calculation of the first and second statistical ratio characteristics includes: Median filtering is applied to the discrete curvature sequence and the discrete torsion sequence respectively to obtain the denoised curvature sequence and the denoised torsion sequence. The arithmetic mean of the denoised curvature sequence is calculated as the curvature mean, and the standard deviation of the denoised curvature sequence is calculated as the curvature standard deviation. The ratio of the curvature mean to the curvature standard deviation is used as the first statistical ratio feature. The arithmetic mean of the denoised torsion sequence is calculated as the torsion mean, and the standard deviation of the denoised torsion sequence is calculated as the torsion standard deviation. The ratio of the torsion mean to the torsion standard deviation is used as the second statistical ratio feature.

4. The method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals according to claim 1, characterized in that, Calculating the fractal dimension of the time-frequency matrix and extracting the contrast and dissimilarity of the gray-level co-occurrence matrix of the time-frequency matrix includes: A short-time Fourier transform is performed on the light intensity sequence using a Hanning window and a specified overlap length to obtain a time-frequency matrix, wherein the rows of the time-frequency matrix correspond to frequencies and the columns correspond to time frames. The fractal dimension of the time-frequency matrix is ​​calculated using the box counting method, and the fractal dimension scalar is obtained. The pixel values ​​of the time-frequency matrix are quantized into gray levels of a specified number to generate a gray-level co-occurrence matrix; Calculate the contrast scalar of the gray-level co-occurrence matrix, and simultaneously calculate the dissimilarity scalar; The difference between the contrast scalar and the dissimilarity scalar is used as the third difference feature.

5. The method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals according to claim 1, characterized in that, Calculating correlation coefficients and mutual information includes: The total energy sequence is obtained by calculating the sum of squares of the amplitudes of all frequencies in each time frame of the time-frequency matrix. Interpolate the curvature sequence to make the number of sampling points equal to the length of the total energy sequence to obtain the resampled curvature sequence, and calculate the correlation coefficient between the resampled curvature sequence and the total energy sequence; Calculate short-time energy sequences for light intensity sequences using a Hanning window and a specified overlap length; The length of the Bonga sphere great circle between adjacent points is calculated based on the polarization state sequence, and the original arc length change rate sequence is obtained. Linear interpolation is performed on the short-time energy sequence and the original arc length change rate sequence to ensure that they have the same number of sampling points on the same time axis. The sequences are then discretized into a specified number of equally spaced intervals, and their mutual information is calculated.

6. The method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals according to any one of claims 3, 4, and 5, characterized in that, Combining and normalizing the feature vectors includes: The first statistical ratio feature, the second statistical ratio feature, the third difference feature, the fractal dimension scalar, the contrast scalar, the correlation coefficient, and the mutual information are concatenated in sequence into an original feature vector; Read the global mean vector and global standard deviation vector of each component from the pre-stored training set statistics file. Subtract the corresponding global mean from each component of the original feature vector and divide by the corresponding global standard deviation to obtain the normalized feature vector.

7. The method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals according to claim 6, characterized in that, The classification results obtained by inputting the classifier include: Extract the sample covariance matrix and feature template vector for each perturbation category from the training set. The feature template vector is the mean of the normalized feature vectors of all samples in the corresponding category. Add a fixed regularization coefficient to the covariance matrix of each category and multiply it by the identity matrix to obtain the regularized covariance matrix; The Mahalanobis distance is calculated between the normalized feature vector to be classified and the feature template vector of each category. The covariance matrix of the Mahalanobis distance is the regularized covariance matrix of the corresponding category. The class with the smallest Mahalanobis distance is selected as the candidate class. At the same time, the difference between the smallest Mahalanobis distance and the second smallest Mahalanobis distance is calculated. If the smallest Mahalanobis distance is less than the mean of the Mahalanobis distances of the corresponding class in the training set plus a specified multiple of the standard deviation, and the difference is greater than a preset difference threshold, then the corresponding candidate class is output; otherwise, the unknown perturbation class is output.

8. The method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals according to claim 2, characterized in that, The spherical linear interpolation uses equal-sized circular arc lengths to generate a fixed number of interpolation trajectory points, with equal arc lengths between adjacent interpolation points; the discrete Frenet formula uses a fixed-point sliding window to calculate the curvature of a discrete point sequence, and the curvature and torsion of the window center point are estimated by the first and second derivatives of the discrete curve, respectively.

9. The method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals according to claim 5, characterized in that, Before calculating the correlation coefficient between the resampled curvature sequence and the total energy sequence, both the resampled curvature sequence and the total energy sequence are zero-mean processed. The linear interpolation of the short-time energy sequence and the original arc length change rate sequence is time-axis aligned so that they have the same number of sampling points on the same time axis, and is discretized into a specified number of intervals using an isofrequency histogram. Mutual information is calculated using the classical entropy formula.

10. The method for perturbation feature extraction and pattern recognition based on microstructured optical fiber sensing signals according to claim 4, characterized in that, The pixel spacing of the gray-level co-occurrence matrix is ​​one, the direction is zero, and the gray level is quantized to a specified number of levels; The contrast scalar is calculated by weighting the sum of the squares of the differences between adjacent gray levels, while the dissimilarity scalar is calculated by weighting the sum of the absolute values ​​of the differences between adjacent gray levels.

Images