A method for extracting vibration perception information features based on tunnel spatial correlation

Abstract

A vibration perception information feature extraction method based on tunnel spatial correlation relates to a vibration perception information extraction method. First, the collected vibration signal is processed and denoised to improve signal quality. Second, feature extraction is performed, and dimension reduction is achieved using difference analysis and correlation analysis. Then, a signal fusion feature correction algorithm based on the spatial correlation of adjacent measuring points is proposed. A distance-related credibility weight calculation method based on quadratic function weight distribution is constructed. The credibility weight calculation is combined with the probability distribution correlation of kernel density estimation to achieve fusion feature value correction of the spatial correlation of adjacent measuring points. Finally, the recognition accuracy is calculated by combining the SVM classification model. By processing and denoising the vibration signal, the quality of the collected signal is improved. Feature extraction and dimensionality reduction are achieved based on difference analysis and correlation analysis. A signal fusion feature correction algorithm based on spatial correlation is proposed to ensure high-precision feature extraction.

Description

Technical Field

[0001] The present invention relates to a vibration perception information extraction method, in particular to a vibration perception information feature extraction method based on tunnel spatial correlation, and belongs to the technical field of subway operation monitoring. Background Art

[0002] Accurately identifying abnormal disturbances in subway tunnel operations and extracting their features can effectively ensure the safety of subway structures. When distributed optical fiber is used to monitor subway cross-section strain, the system's high sensitivity often causes the collected monitoring signals to contain considerable noise, which seriously affects the extraction and discovery of effective information and makes it more difficult to identify the vibration time of abnormal disturbances during subway operations.

[0003] Therefore, in order to solve the problem of how to accurately extract the most effective features of vibration perception information, a vibration perception information feature extraction method based on the spatial correlation of subway tunnels is urgently needed to accurately extract the abnormal disturbance information features, thereby eliminating the safety hazards caused by abnormal disturbances in the subway protection zone and ensuring the safe operation of the subway. Summary of the Invention

[0004] To address the shortcomings of the background technology, the present invention provides a vibration perception information feature extraction method based on tunnel spatial correlation. It improves the quality of the collected signal by processing and denoising the vibration signal, then realizes feature extraction and dimensionality reduction based on difference analysis and correlation analysis, and proposes a signal fusion feature correction algorithm based on spatial correlation to ensure high-precision feature extraction.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a method for extracting vibration perception information features based on tunnel spatial correlation, comprising the following steps:

[0006] Step 1: Lay out in subway tunnels Distributed optical fiber collects vibration signals, processes outliers and normalizes them, and uses wavelet thresholding to remove noise, achieving data preprocessing and noise reduction of vibration signals.

[0007] Step 2: Extract features from the processed vibration signal based on the time and frequency domains. Use difference analysis and correlation analysis to process features, thus achieving multi-dimensional feature extraction and dimensionality reduction of the vibration signal.

[0008] Step 3: Construct a distance-related credibility weight calculation method based on quadratic function weight distribution, combine it with the probability distribution correlation of kernel density estimation to calculate the credibility weight, complete the information fusion of the adjacent measurement points, and realize the algorithm correction of the signal fusion feature of spatial correlation. Specifically:

[0009] 3.1 Construct a signal fusion feature correction algorithm based on spatial correlation. The process is as follows:

[0010] (1) Calculate the characteristic value vector T of the signal collected at the target measurement point p (p) :

[0011]

[0012] Where, Indicates the nth characteristic calculation value of the characteristic vector of the signal collected at the target measurement point p;

[0013] According to the characteristic calculation values ​​of the target measuring point p and the 2m+1 measuring points before and after it, the characteristic calculation value matrix of the target measuring point p range is formed:

[0014]

[0015] Among them, record:

[0016]

[0017] Then, the target measurement point p range characteristic calculation value matrix is ​​expressed as:

[0018]

[0019] Where, Represents the target measurement point p range characteristic calculation value matrix, with a dimension of (2m+1)×n, Represents the nth feature calculation value vector within the target measurement point p, with a dimension of (2m+1)×1;

[0020] (2) Establish a distance-related credibility weight vector based on the distance between the target measuring point p and the 2m+1 measuring points before and after it.

[0021]

[0022]

[0023] Where, It represents the reliability weight coefficient of the distance between the mth measuring point and the target measuring point p. Represents the reliability weight coefficient of the distance between the mth measuring point and the target measuring point p;

[0024] (3) Calculate the characteristic values ​​of the target measurement point p and the 2m+1 measurement points before and after it, and establish the probability distribution related credibility weight coefficient for each feature through probability density estimation Then establish the probability distribution related credibility weight matrix of the feature population first:

[0025]

[0026]

[0027] Where, Represents the probability distribution related credibility weight vector of the nth feature of the total 2m+1 measurement points before and after the target measurement point p, with a dimension of 1×(2m+1)×1. Represents the probability distribution related credibility weight coefficient of the nth feature of the mth measurement point before the target measurement point p, Represents the credibility weight coefficient of the probability distribution of the nth feature of the mth measurement point after the target measurement point p;

[0028] but:

[0029]

[0030] Where, Represents the probability distribution of the feature population and the associated credibility weight matrix, with dimension n×(2m+1);

[0031] (4) Taking into account the distance-related credibility weight and the probability distribution-related credibility weight, the comprehensive weight of each feature before and after the measurement point is calculated and normalized as the basis for correcting the feature calculation value. The process is as follows:

[0032] For the nth feature before and after the target measurement point p, the Hadamard product of the distance-related credibility weight vector and the probability distribution-related credibility weight vector is calculated and normalized to obtain the normalized comprehensive weight vector of the nth feature before and after the target measurement point p. First:

[0033]

[0034] The calculation method of each element is as follows:

[0035]

[0036] Where, Represents the comprehensive weight vector of the nth feature before and after the target measurement point p, Represents the Hadamard product calculation symbol;

[0037] Normalize the comprehensive weight vector to obtain the normalized comprehensive weight vector:

[0038]

[0039] The normalization calculation method for each element is as follows:

[0040]

[0041] Where, Represents the normalized comprehensive weight vector of the nth feature before and after the target measurement point p;

[0042] (5) Take the characteristic calculation value of the target measuring point p and the signals collected from a total of 2m+1 measuring points, and correct the characteristic calculation value of the target measuring point p. The process is as follows:

[0043] According to the previous steps (1)-(4), the normalized comprehensive weight vector of each feature before and after the target measurement point p is calculated in sequence to form the normalized comprehensive weight matrix of the target measurement point p feature:

[0044]

[0045] Then, the nth characteristic calculation value T of the target measurement point p is n (p) Correction value It is expressed as follows:

[0046]

[0047] Then, the nth characteristic calculation value vector T of the target measurement point p is (p) Correction vector It is expressed as follows:

[0048]

[0049] Where E represents the identity matrix of dimension n;

[0050] 3.2 Calculate the distance-related credibility weight for quadratic function weight distribution. The process is as follows:

[0051] (1) Since the distance relationship between different measuring points increases linearly, it is assumed that as the distance increases, the credibility of the measuring point in correcting the target measuring point p will gradually decay in the trend of a quadratic function. This decay is expressed by a quadratic function:

[0052]

[0053] Where, It represents the credibility of the qth measuring point when correcting the target measuring point p, a represents the credibility attenuation coefficient, and a is taken as -1 / (2m) 2 , w max Indicates the reliability of the measuring point when correcting its own information;

[0054] (2) Calculate the credibility of each measuring point when correcting the target measuring point p information, and normalize the weights so that the sum of all weights is 1, forming a credibility weight vector related to the distance between the measuring points before and after the target measuring point p. First:

[0055]

[0056] Where, represents the credibility weight of the qth measuring point when correcting the target measuring point p;

[0057] but:

[0058]

[0059] Where, Represents the distance-related credibility weight vector of the target measurement point p;

[0060] 3.3 The credibility weight is calculated based on the probability distribution of kernel density estimation. The process is as follows:

[0061] (1) Calculate the characteristic calculation values ​​of all the features of the target measuring point p and the 2m+1 measuring points before and after it, and obtain For the nth feature, the feature calculation value of the target measurement point p and the 2m+1 measurement points before and after is

[0062] (2) Perform kernel density estimation by calculating the kernel to obtain the probability density function f of the nth feature calculation value n (T) describes the probability distribution of the nth characteristic calculation value before and after the target measurement point p, and is calculated using the following formula:

[0063]

[0064] Where K(x) represents the kernel function of the Epanechnikov kernel;

[0065] (3) f n (T) Interpolation is performed to calculate the probability of occurrence of the characteristic calculation value of the target measuring point p and the total 2m+1 measuring points before and after it, and this value is used as the credibility weight coefficient related to the probability distribution. The credibility weight coefficient related to the distance of the qth measuring point near the target measuring point p is as follows:

[0066]

[0067] Where, Represents the confidence weight coefficient associated with the probability distribution of the nth feature calculation value of the qth measurement point near the target measurement point p;

[0068] (4) Establishing the credibility weight of the probability distribution of the feature population as follows:

[0069]

[0070] Thus, the probability distribution of the feature population is established with a confidence weight Characterize the probability distribution of the calculated characteristic values ​​of the target measuring point before and after the target measuring point p, so as to correct the calculated characteristic value of the collected signal of the target measuring point p;

[0071] Step 4: According to steps 1, 2, and 3, the signal fusion feature correction algorithm based on spatial correlation is used to process the feature calculation value, and the SVM classification model is used to calculate the model to obtain the recognition accuracy of the signal feature.

[0072] Compared with the existing technology, the beneficial effects of the present invention are: the present invention is aimed at distributed optical fiber collection of vibration signals, and studies the noise reduction algorithm of subway tunnel structure vibration signals based on wavelet decomposition technology. By processing and denoising the vibration signals, the quality of the collected signals is improved. For vibration event pattern recognition, vibration signals can be extracted, and feature extraction and dimensionality reduction are realized based on correlation and difference analysis. The spatial correlation characteristics between adjacent measuring points of distributed optical fiber collection vibration signals are utilized to propose a signal fusion feature correction algorithm based on spatial correlation. Compared with the existing method, the accuracy and efficiency of feature extraction of abnormal disturbance vibration information in tunnels are improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 is a flow chart of the present invention;

[0074] Figure 2 is a signal waterfall diagram collected in the embodiment;

[0075] Figure 3 is a comparative diagram of the distribution of vibration characteristics on different vibration events in the embodiment;

[0076] Figure 4 is a Pearson correlation coefficient graph of consecutive vibration events in the embodiment;

[0077] Figure 5 is a Pearson correlation coefficient graph of shock vibration events in the embodiment;

[0078] Figure 6 is a Pearson correlation coefficient graph of continuous + shock vibration events in the embodiment;

[0079] Figure 7 is a Pearson correlation coefficient graph of the embodiment without interference events;

[0080] Figure 8 is a graph of the SVM recognition accuracy of the uncorrected features in the embodiment;

[0081] Figure 9 is a graph of the SVM recognition accuracy of the corrected features in the embodiment. DETAILED DESCRIPTION

[0082] The technical solutions of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0083] like Figure 1 As shown, a method for extracting vibration perception information features based on tunnel spatial correlation includes the following steps:

[0084] Step 1: Lay out in subway tunnels Distributed optical fiber collects vibration signals, processes outliers and normalizes them, and uses the wavelet threshold method to remove noise, thus achieving data preprocessing and noise reduction of vibration signals. Specifically:

[0085] 1.1 First, process the abnormal value of the collected vibration signal. The process is as follows:

[0086] (1) Select the calculation quartile to establish the recognition threshold of the outlier, sort the vibration signals according to the amplitude, and calculate the first quartile Q1 and the third quartile Q3;

[0087] (2) Calculate the recognition threshold of outliers (y min ,y max )as follows:

[0088] y max =Q3+3×(Q3-Q1)

[0089] y min =Q1-3×(Q3-Q1)

[0090] will exceed (y min ,y max )’s vibration signal is marked as an outlier;

[0091] (3) Delete outlier points.

[0092] 1.2 The vibration signal after outlier processing is normalized using the centering and standardization method. The process is as follows:

[0093] (1) Centralization is to subtract the average value from the vibration signal after outlier processing, using the following formula:

[0094]

[0095]

[0096] Where y urepresents the u-th vibration signal after outlier processing, y represents the average value of the vibration signal after outlier processing, N represents the number of sampling points of the vibration signal after outlier processing, y' u represents the u-th vibration signal after centralization;

[0097] (2) The formula used for standardization is as follows:

[0098]

[0099] Where y u,sta represents the u-th vibration signal after normalization, and σ represents the standard deviation of the vibration signal.

[0100] 1.3 Use wavelet threshold method to remove noise. The process is as follows:

[0101] First, the coif wavelet-based noise reduction process is used. This process can significantly reduce the noise signal, improve the signal-to-noise ratio, and improve the signal quality without losing much information, which is convenient for subsequent processing.

[0102] Secondly, a threshold is set for the obtained wavelet coefficients, and the wavelet coefficients with absolute values ​​less than the threshold are assigned a smaller value, preferably 0, to achieve the effect of noise reduction. The effect is better when the vibration signal is processed. When the vibration signal is detected, the Rigrsure threshold method is used to reduce the noise of the vibration signal;

[0103] Finally, a suitable threshold function is selected to further optimize the noise reduction effect. The improved threshold function used in this invention is as follows:

[0104]

[0105] Where TH represents the set threshold, f represents the wavelet coefficient, and Sign represents the Sign function.

[0106] The vibration signal after outlier processing and normalization is subjected to the above steps in sequence to remove noise.

[0107] Step 2: Extract features from the processed vibration signal based on the time and frequency domains. Use difference analysis and correlation analysis to process features, and achieve multi-dimensional feature extraction and dimensionality reduction of the vibration signal. Specifically:

[0108] 2.1 Extract the multi-dimensional eigenvalues ​​of the vibration signal from the time domain and frequency domain perspectives, and summarize the common T features.

[0109] 2.2 One-way ANOVA was used to analyze the differences. The process is as follows:

[0110] (1) First, the features are organized into a feature matrix, where each row represents a feature and each column represents a sample, and a label vector with a length equal to the number of samples is created;

[0111] (2) Perform a one-way analysis of variance for each feature, calculate the differences of each feature between different vibration event categories, and output a W value. The W value represents the probability of observing the actual difference under the null hypothesis;

[0112] (3) Select the differential features according to the magnitude of the W value for subsequent classifier training.

[0113] 2.3 Perform correlation analysis and evaluation based on the Pearson correlation coefficient method, and the process is as follows:

[0114] (1) Use the extracted features to construct a feature matrix;

[0115] (2) Calculate the Pearson correlation coefficient between each pair of features in the feature matrix to form a correlation coefficient matrix;

[0116] (3) Since it is difficult to express the correlation between features, the present invention uses the method of generating a heat map for display. The Pearson correlation coefficient between features is represented by the data interval [-1, 1]. The closer the Pearson coefficient is to -1, the stronger the negative correlation relationship; the closer it is to +1, the stronger the positive correlation relationship; and the closer it is to 0, the weaker the correlation relationship.

[0117] Step 3: Construct a method for calculating weights based on the distance-related credibility of quadratic function weight allocation, combine the probability distribution related to kernel density estimation for calculating credibility weights, complete the information fusion of adjacent front and rear measurement points, and realize the algorithm correction of the signal fusion feature of spatial correlation. Specifically:

[0118] 3.1 Construct an algorithm for correcting the signal fusion feature based on spatial correlation, and the process is as follows:

[0119] (1) Calculate the feature calculation value vector T of the signal collected by the target measurement point p (p) :

[0120]

[0121] In the formula, represents the nth feature calculation value of the feature vector of the signal collected by the target measurement point p.

[0122] According to the feature calculation values of the target measurement point p and a total of 2m + 1 (m < p) measurement points before and after, form the target measurement point p range feature calculation value matrix:

[0123]

[0124] ​

[0125]

[0126] Then, the target measurement point p range characteristic calculation value matrix is ​​expressed as:

[0127]

[0128] Where, Represents the target measurement point p range characteristic calculation value matrix, with a dimension of (2m+1)×n, Represents the nth feature calculation value vector within the target measurement point p, with a dimension of (2m+1)×1;

[0129] (2) Establish a distance-related credibility weight vector based on the distance between the target measuring point p and the 2m+1 measuring points before and after it.

[0130]

[0131] Where, It represents the reliability weight coefficient of the distance between the mth measuring point and the target measuring point p. Represents the reliability weight coefficient of the distance between the mth measuring point and the target measuring point p;

[0132] (3) Calculate the characteristic values ​​of the target measurement point p and the 2m+1 measurement points before and after it, and establish the probability distribution related credibility weight coefficient for each feature through probability density estimation Then establish the probability distribution related credibility weight matrix of the feature population first:

[0133]

[0134] Where, Represents the probability distribution related credibility weight vector of the nth feature of the total 2m+1 measurement points before and after the target measurement point p, with a dimension of 1×(2m+1)×1. Represents the probability distribution related credibility weight coefficient of the nth feature of the mth measurement point before the target measurement point p, Represents the credibility weight coefficient of the probability distribution of the nth feature of the mth measurement point after the target measurement point p.

[0135] but:

[0136]

[0137] Where, Represents the probability distribution of the feature population and the associated credibility weight matrix, with dimension n×(2m+1);

[0138] (4) Taking into account the distance-related credibility weight and the probability distribution-related credibility weight, the comprehensive weight of each feature before and after the measurement point is calculated and normalized as the basis for correcting the feature calculation value. The process is as follows:

[0139] For the nth feature before and after the target measurement point p, the Hadamard product of the distance-related credibility weight vector and the probability distribution-related credibility weight vector is calculated and normalized to obtain the normalized comprehensive weight vector of the nth feature before and after the target measurement point p. First:

[0140]

[0141] The calculation method for each element is as follows:

[0142]

[0143] Where, Represents the comprehensive weight vector of the nth feature before and after the target measurement point p, Indicates the Hadamard product calculation symbol.

[0144] Normalize the comprehensive weight vector to obtain the normalized comprehensive weight vector:

[0145]

[0146] The normalization calculation method for each element is as follows:

[0147]

[0148] Where, Represents the normalized comprehensive weight vector of the nth feature before and after the target measurement point p;

[0149] (5) Take the characteristic calculation value of the target measuring point p and the signals collected from a total of 2m+1 measuring points, and correct the characteristic calculation value of the target measuring point p. The process is as follows:

[0150] According to the previous steps (1)-(4), the normalized comprehensive weight vector of each feature before and after the target measurement point p is calculated in sequence to form the normalized comprehensive weight matrix of the target measurement point p feature:

[0151]

[0152] Then, the calculated value of the nth characteristic of the target measuring point p is Correction value It is expressed as follows:

[0153]

[0154] Then, the nth characteristic calculation value vector T of the target measurement point p is (p) Correction vector It is expressed as follows:

[0155]

[0156] Where E represents the identity matrix of dimension n;

[0157] 3.2 Calculate the distance-related credibility weight for quadratic function weight distribution. The process is as follows:

[0158] (1) Since the distance relationship between different measuring points increases linearly, it is assumed that as the distance increases, the credibility of the measuring point in correcting the target measuring point p will gradually decay in the trend of a quadratic function. This decay is expressed by a quadratic function:

[0159]

[0160] Where, It represents the credibility of the qth measuring point when correcting the target measuring point p, a represents the credibility attenuation coefficient, and a is taken as -1 / (2m) 2 , w max Indicates the reliability of the measuring point when correcting its own information;

[0161] (2) Calculate the credibility of each measuring point when correcting the target measuring point p information, and normalize the weights so that the sum of all weights is 1, forming a credibility weight vector related to the distance between the measuring points before and after the target measuring point p. First:

[0162]

[0163] Where, It represents the credibility weight of the qth measuring point when correcting the target measuring point p.

[0164] but:

[0165]

[0166] Where, Represents the distance-related credibility weight vector of the measurement points before and after the target measurement point p.

[0167] 3.3 The credibility weight is calculated based on the probability distribution of kernel density estimation. The process is as follows:

[0168] (1) Calculate the characteristic calculation values ​​of all the features of the target measuring point p and the 2m+1 measuring points before and after it, and obtain For the nth feature, the feature calculation value of the target measurement point p and the 2m+1 measurement points before and after is

[0169] (2) Perform kernel density estimation by calculating the kernel to obtain the probability density function f of the nth feature calculation value n (T) describes the probability distribution of the nth characteristic calculation value before and after the target measurement point p, and is calculated using the following formula:

[0170]

[0171] Where K(x) represents the kernel function of the Epanechnikov kernel;

[0172] (3) f n (T) Interpolation is performed to calculate the probability of occurrence of the characteristic calculation value of the target measuring point p and the total 2m+1 measuring points before and after it, and this value is used as the credibility weight coefficient related to the probability distribution. The credibility weight coefficient related to the distance of the qth measuring point near the target measuring point p is as follows:

[0173]

[0174] Where, Represents the confidence weight coefficient associated with the probability distribution of the nth feature calculation value of the qth measurement point near the target measurement point p;

[0175] (4) Establishing the credibility weight of the probability distribution of the feature population as follows:

[0176]

[0177] Thus, the probability distribution of the feature population is established with a confidence weight Characterize the probability distribution related credibility weight of the characteristic calculation values ​​of the measuring points before and after the target measuring point p, so as to correct the characteristic calculation value of the collected signal of the target measuring point p.

[0178] Step 4: Extract the signal fusion feature correction algorithm based on spatial correlation according to steps 1, 2, and 3, and process the feature calculation value. Use the SVM classification model to calculate the model to obtain the recognition accuracy of the signal feature. Specifically:

[0179] 4.1 Take the T features extracted in step 2 as the feature parameters to be corrected, analyze these T features based on the one-way analysis of variance method in step 2.2 and summarize them into a feature matrix;

[0180] 4.2 Based on the correlation analysis evaluation in step 2.3, extract the calculated values ​​of the features and construct a feature matrix. Perform correlation analysis on each feature in each vibration event and draw a Pearson correlation coefficient graph.

[0181] 4.3 Based on the Pearson correlation coefficient graph, highly correlated features and similar strongly correlated features are eliminated, and the remaining Y features are retained for classifier training to complete feature fusion;

[0182] 4.4 Based on the signal fusion feature correction algorithm based on spatial correlation proposed in step 3, the Y features retained in step 4.3 are fused and corrected;

[0183] 4.5 Using the fused Y features as key features, use SVM to build a recognition model;

[0184] 4.6 Use the fused Y key features as the sample set, randomly select a part of them as the training set, and train the remaining part as the test set. At the same time, establish a comparison group. The comparison group uses the unfused and corrected values ​​of these Y features. Train the SVM classifiers for the test set and the comparison group separately. Record the events recognized by the two trained SVM classifiers, and calculate the recognition accuracy of each event and the overall event for comparison.

[0185] Example

[0186] Built in the laboratory The distributed fiber optic disturbance sensing system uses a 1km long single-mode optical fiber as the propagation medium and a distributed fiber acoustic wave detector (DAS) as the light source and data acquisition device. The laid optical fiber is a single-mode single-core optical fiber with a total length of 1km. Vibration signals are not collected for about 500m after the jumper access point to ensure stability during the signal acquisition process.

[0187] Four vibration signals are simulated: environmental disturbance, impact vibration, continuous vibration, and impact + continuous vibration. When collecting environmental disturbance signals, the input vibration signal can be directly collected without assuming there is no interference. Impact vibration is simulated by striking the tabletop with a hammer. After a single strike, the hammer is immediately moved away from the tabletop, and a signal of the gradually decaying vibration is collected. Continuous vibration signals are simulated by manually controlling the hammer to strike the tabletop continuously at a certain frequency. Impact + continuous vibration signals are simulated by using a hammer to randomly combine both impacts and continuous strikes.

[0188] The theoretical methods from step 1.1 to step 1.2 are used to perform outlier processing and normalization on the collected vibration signals respectively.

[0189] After completing the above data processing, the collected signal is denoised using the wavelet threshold method in step 1.3 to obtain a high-quality Vibration signal data, combined with Figure 2 The following is a waterfall diagram of a signal collected after the above processing. After the above data collection and processing, a total of 620 sets of data are obtained. The specific situation is shown in Table 1:

[0190] Table 1 Experimental data collection

[0191]

[0192] According to the multi-dimensional characteristic values ​​of the disturbance signal extracted from multiple aspects of time domain and frequency domain, a variety of common signal characteristic parameters are sorted out and summarized. In this implementation case, 20 basic features are selected as the features to be corrected, namely: maximum value (Max), minimum value (Min), average value (Me), mean peak value (PM), rectified average value (Arv), variance (Var), standard deviation (Std), kurtosis (Kurt), skewness (Skew), root mean square (RMS), form factor (Sh), peak factor (Cr), impulse factor (Im), margin factor (Cl), center of gravity frequency (CF), mean square frequency (MSF), frequency variance (RMSF), band energy (BE), relative power spectrum entropy (RPSD) and peak to average ratio (Pars).

[0193] The peak-to-average ratio in this embodiment is improved to more accurately reveal the dynamic changes of the optical fiber sensing signal and improve the accuracy and reliability of signal processing. The improvement process is as follows:

[0194] The improved peak-to-average ratio can select the 85% upper quantile of the signal absolute amplitude when calculating, so as to reduce the influence of the local large noise signal on the data feature and improve the discrimination of the feature. The improved peak-to-average ratio can better meet the requirements of the proposed method. Processing and feature analysis of collected data.

[0195] Step 1: Calculate the absolute value of the signal |y(t)|;

[0196] Step 2: Select the part above the 85% quantile in |y(t)| and calculate the average value as the peak value y of the signal mean,peak ;

[0197] Step 3: Calculate the peak-to-average ratio:

[0198]

[0199] Calculate the characteristic values ​​of four vibration events, among which Table 2 is the feature extraction result of continuous vibration event, Table 3 is the feature extraction result of impact vibration, Table 4 is the feature extraction result of impact + continuous vibration, and Table 5 is the feature extraction result of non-interference event. Figure 3 As shown in Figure 3, a comparative distribution diagram of vibration features across different vibration events is plotted. Direct observation of the calculated feature values ​​extracted from the four vibration events reveals significant differences in the different features across different vibration events, which to some extent reflects the overall ability of the features studied in this paper to distinguish different signals.

[0200] Table 2 Continuous vibration event signal feature extraction results

[0201]

[0202] Table 3. Results of feature extraction of shock vibration event signals

[0203]

[0204]

[0205] Table 4. Feature extraction results of shock + continuous vibration event signals

[0206]

[0207] Table 5. Results of feature extraction of non-interference event signals

[0208]

[0209]

[0210] The 20 extracted features were used as the feature parameters to be corrected. Based on the algorithm in step 2.2, this example refers to the general W value limit of 0.05. That is, when the W value is less than 0.05, the feature is considered to be of great help in the classification task. These 20 features were analyzed using a one-way analysis of variance method and summarized into a feature matrix. The calculated W values ​​are shown in Table 6 below. The results show that the selected 20 features can well describe the differences in events.

[0211] Table 6 W-value of feature difference analysis

[0212]

[0213] Based on the method proposed in step 2.3, the characteristic calculation values ​​of the four vibration events are extracted to construct a characteristic matrix. The correlation analysis of each feature is performed in each vibration event. The results are calculated based on the Pearson correlation coefficient. Figures 4 to 7 ,in Figure 4 is the Pearson correlation coefficient diagram of continuous vibration events, Figure 5 is the Pearson correlation coefficient diagram of shock vibration events, Figure 6 is the Pearson correlation coefficient diagram of continuous + shock vibration events, Figure 7 The Pearson correlation coefficient diagram of the non-interference event. Comprehensively analyze the values ​​of the correlation coefficients of different features in the four types of events, make an overall evaluation of the correlation between the features, and combine Figures 4 to 7 It can be seen that the characteristics of the vibration signal have the following correlations:

[0214] (1) Among the four types of vibration events, there is a strong correlation between Var, Std, and RMS, and one of the three can be selected;

[0215] (2) Among the four types of vibration events, there is a strong correlation between Sh, Cr, Im, CI, and Kurt, and some features can be considered to be discarded;

[0216] (3) Among the four types of vibration events, there is a strong correlation between CF, MSF, RMSF, and BE. It is possible to consider discarding some features to improve the efficiency of subsequent calculations;

[0217] (4) Some of the other features not mentioned above have strong correlations with other features when analyzing a certain vibration event, but no strong correlations are found when analyzing other vibration events. It can be inferred that this phenomenon is caused by the nature of the event itself, rather than defects in the features themselves. Therefore, highly correlated features and similar strongly correlated features are eliminated, that is, the five features Var, Std, Sh, Kurt and CI are discarded, and the remaining 15 features are retained for subsequent classifier training.

[0218] Based on the signal fusion feature correction algorithm based on spatial correlation proposed in step 3, the 15 features retained by each vibration event data sample are fused and corrected. In this implementation case, each measurement point uses the information of the 10 measurement points before and after for correction during calculation, that is, the value of m is 10. When calculating the distance-related credibility weight, the credibility attenuation coefficient a is -0.0025, and the credibility w of the measurement point when correcting its own information is max When calculating the credibility weight associated with the probability distribution based on kernel density estimation, the Epanechnikov kernel is selected as the calculation kernel for kernel density estimation. Similarly, the characteristic calculation values ​​of the 10 measurement points before and after the measurement point are used to estimate the probability distribution of the characteristic values.

[0219] The 15 fused features were used as key features, combined with an SVM-based model for recognition. The extracted features served as the sample set, with 75% randomly selected as the training set and the remaining 25% as the test set. This trained model was designated SVM-I. A comparison group was also established, using the unfused, corrected values ​​of the 15 features to train the SVM classifier. This trained model was designated SVM-II. Events accurately recognized by the SVM classifiers from both training runs were recorded, and the individual and overall recognition accuracy rates were calculated for comparison. Table 7 shows the SVM recognition accuracy comparison.

[0220] Table 7 Comparison of SVM recognition accuracy using features trained before and after correction

[0221]

[0222]

[0223] SVM recognition effect comparison chart combined Figures 8 and 9 As shown, Figure 8 is the SVM recognition accuracy graph without corrected features, Figure 9 Figure 2 shows the SVM recognition accuracy of features corrected by the method of the present invention. By comparison, it was found that the signal fusion feature correction algorithm based on spatial correlation proposed in the present invention can effectively identify the calculated signal feature values. The recognition accuracy of the trained model has increased by about 5%, and the recognition accuracy of some events has increased by 10%. The overall recognition accuracy has increased from 82.5% to 87.74%, a significant improvement. Therefore, the signal fusion feature correction algorithm based on spatial correlation proposed in the present invention can effectively improve the quality of the calculated signal feature values ​​and can directly improve the recognition rate of the trained classifier.

[0224] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above and that the invention can be implemented in other configurations without departing from the spirit or essential characteristics of the invention. Therefore, the embodiments should be considered in all respects as illustrative and non-restrictive, and the scope of the invention is defined by the appended claims, not the foregoing description, and all variations coming within the meaning and range of equivalents of the claims are intended to be embraced therein. Any reference sign in a claim should not be construed as limiting the claim to which it relates.

[0225] In addition, it should be understood that although this specification is described in terms of implementation methods, not every implementation method contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other implementation methods that can be understood by those skilled in the art.

Claims

1. A method for extracting vibration perception information features based on tunnel spatial correlation, characterized by: The following steps are involved: Step 1: Lay out in subway tunnels -OTDR distributed optical fiber collects vibration signals, processes outliers and normalizes them, and uses wavelet threshold method to remove noise, thus achieving data preprocessing and noise reduction of vibration signals; Step 2: Extract features from the processed vibration signal based on the time and frequency domains. Use difference analysis and correlation analysis to process features, thus achieving multi-dimensional feature extraction and dimensionality reduction of the vibration signal. Step 3: Construct a distance-related credibility weight calculation method based on quadratic function weight distribution, combine it with the probability distribution correlation of kernel density estimation to calculate the credibility weight, complete the information fusion of the adjacent measurement points, and realize the algorithm correction of the signal fusion feature of spatial correlation. Specifically: 3.1 Construct a signal fusion feature correction algorithm based on spatial correlation. The process is as follows: (1) Calculate the characteristic value vector T of the signal collected at the target measurement point p (p) : Where, Indicates the nth characteristic calculation value of the characteristic vector of the signal collected at the target measurement point p; According to the characteristic calculation values ​​of the target measuring point p and the 2m+1 measuring points before and after it, the characteristic calculation value matrix of the target measuring point p range is formed: Among them, record: Then, the target measurement point p range characteristic calculation value matrix is ​​expressed as: Where, Represents the target measurement point p range characteristic calculation value matrix, with a dimension of (2m+1)×n, Represents the nth feature calculation value vector within the target measurement point p, with a dimension of (2m+1)×1; (2) Establish a distance-related credibility weight vector based on the distance between the target measuring point p and the 2m+1 measuring points before and after it. Where, It represents the reliability weight coefficient of the distance between the mth measuring point and the target measuring point p. Represents the reliability weight coefficient of the distance between the mth measuring point and the target measuring point p; (3) Calculate the characteristic values ​​of the target measurement point p and the 2m+1 measurement points before and after it, and establish the probability distribution related credibility weight coefficient for each feature through probability density estimation Then establish the probability distribution related credibility weight matrix of the feature population first: Where, Represents the probability distribution related credibility weight vector of the nth feature of the total 2m+1 measurement points before and after the target measurement point p, with a dimension of 1×(2m+1)×1. Represents the probability distribution related credibility weight coefficient of the nth feature of the mth measurement point before the target measurement point p, Represents the credibility weight coefficient of the probability distribution of the nth feature of the mth measurement point after the target measurement point p; but: Where, Represents the probability distribution of the feature population related credibility weight matrix, dimension is n × (2m + 1); (4) Taking into account the distance-related credibility weight and the probability distribution-related credibility weight, the comprehensive weight of each feature before and after the measurement point is calculated and normalized as the basis for correcting the feature calculation value. The process is as follows: For the nth feature before and after the target measurement point p, the Hadamard product of the distance-related credibility weight vector and the probability distribution-related credibility weight vector is calculated and normalized to obtain the normalized comprehensive weight vector of the nth feature before and after the target measurement point p. First: The calculation method for each element is as follows: Where, Represents the comprehensive weight vector of the nth feature before and after the target measurement point p, Represents the Hadamard product calculation symbol; Normalize the comprehensive weight vector to obtain the normalized comprehensive weight vector: The normalization calculation method for each element is as follows: Where, Represents the normalized comprehensive weight vector of the nth feature before and after the target measurement point p; (5) Take the characteristic calculation value of the target measuring point p and the signals collected from a total of 2m+1 measuring points, and correct the characteristic calculation value of the target measuring point p. The process is as follows: According to the previous steps (1)-(4), the normalized comprehensive weight vector of each feature before and after the target measurement point p is calculated in sequence to form the normalized comprehensive weight matrix of the target measurement point p feature: Then, the calculated value of the nth characteristic of the target measuring point p is Correction value It is expressed as follows: Then, the nth characteristic calculation value vector T of the target measurement point p is (p) Correction vector It is expressed as follows: Where E represents the identity matrix of dimension n; 3.2 Calculate the distance-related credibility weight for quadratic function weight distribution. The process is as follows: (1) Since the distance relationship between different measuring points increases linearly, it is assumed that as the distance increases, the credibility of the measuring point in correcting the target measuring point p will gradually decay in the trend of a quadratic function. This decay is expressed by a quadratic function: Where, It represents the credibility of the qth measuring point when correcting the target measuring point p, a represents the credibility attenuation coefficient, and a is taken as -1 / (2m) 2 , w max Indicates the reliability of the measuring point when correcting its own information; (2) Calculate the credibility of each measuring point when correcting the target measuring point p information, and normalize the weights so that the sum of all weights is 1, forming a credibility weight vector related to the distance between the measuring points before and after the target measuring point p. First: Where, represents the credibility weight of the qth measuring point when correcting the target measuring point p; but: Where, Represents the distance-related credibility weight vector of the target measurement point p; 3.3 The credibility weight is calculated based on the probability distribution of kernel density estimation. The process is as follows: (1) Calculate the characteristic calculation values ​​of all the features of the target measuring point p and the 2m+1 measuring points before and after it, and obtain For the nth feature, the feature calculation value of the target measurement point p and the 2m+1 measurement points before and after is (2) Perform kernel density estimation by calculating the kernel to obtain the probability density function f of the nth feature calculation value n (T) describes the probability distribution of the nth characteristic calculation value before and after the target measurement point p, and is calculated using the following formula: Where K(x) represents the kernel function of the Epanechnikov kernel; (3) f n (T) Interpolation is performed to calculate the probability of occurrence of the characteristic calculation value of the target measuring point p and the total 2m+1 measuring points before and after it, and this value is used as the credibility weight coefficient related to the probability distribution. The credibility weight coefficient related to the distance of the qth measuring point near the target measuring point p is as follows: Where, Represents the confidence weight coefficient associated with the probability distribution of the nth feature calculation value of the qth measurement point near the target measurement point p; (4) Establishing the credibility weight of the probability distribution of the feature population as follows: Thus, the probability distribution of the feature population is established with a confidence weight Characterize the probability distribution of the calculated characteristic values ​​of the target measuring point before and after the target measuring point p, so as to correct the calculated characteristic value of the collected signal of the target measuring point p; Step 4: According to steps 1, 2, and 3, the signal fusion feature correction algorithm based on spatial correlation is used to process the feature calculation value, and the SVM classification model is used to calculate the model to obtain the recognition accuracy of the signal feature.

2. The method for extracting vibration perception information features based on tunnel spatial correlation according to claim 1, characterized in that: The outlier processing process in step 1 is as follows: (1) Select the calculation quartile to establish the recognition threshold of the outlier, sort the vibration signals according to the amplitude, and calculate the first quartile Q1 and the third quartile Q3; (2) Calculate the recognition threshold of outliers (y min ,y max )as follows: y max =Q3+3×(Q3-Q1) y min =Q1-3×(Q3-Q1) will exceed (y min ,y max )’s vibration signal is marked as an outlier; (3) Delete outlier points.

3. The method for extracting vibration perception information features based on tunnel spatial correlation according to claim 2, characterized in that: The normalization process in step 1 adopts a centralized and standardized approach, and the process is as follows: (1) Centralization is to subtract the average value from the vibration signal after outlier processing, using the following formula: Where y u represents the u-th vibration signal after outlier processing, represents the average value of the vibration signal after outlier processing, N represents the number of sampling points of the vibration signal after outlier processing, y' u represents the u-th vibration signal after centralization; (2) The formula used for standardization is as follows: Where y u,sta represents the u-th vibration signal after normalization, and σ represents the standard deviation of the vibration signal.

4. The method for extracting vibration perception information features based on tunnel spatial correlation according to claim 3, characterized in that: The process of removing noise using the wavelet threshold method in step 1 is as follows: First, the coif wavelet basis is used for noise reduction; Secondly, a threshold is set for the obtained wavelet coefficients, and the wavelet coefficients with absolute values ​​less than the threshold are assigned a smaller value to achieve the effect of noise reduction. -When detecting OTDR vibration signals, select the Rigrsure threshold method to reduce the noise of the vibration signal; Finally, the threshold function is selected to optimize the noise reduction effect. The improved threshold function is as follows: Where TH represents the set threshold, f represents the wavelet coefficient, and Sign represents the Sign function.

5. The method for extracting vibration perception information features based on tunnel spatial correlation according to claim 1, characterized in that: The step 2 is specifically as follows: 2.1 Extract the multi-dimensional eigenvalues ​​of the vibration signal from the time domain and frequency domain perspectives, and summarize T features; 2.2 One-way ANOVA was used to analyze the differences. The process is as follows: (1) First, the features are organized into a feature matrix, where each row represents a feature and each column represents a sample, and a label vector with a length equal to the number of samples is created; (2) Perform one-way ANOVA on each feature, calculate the differences between different vibration event categories and output a W value, which represents the probability of observing actual differences under the null hypothesis; (3) Select the features with differences according to the value of W; 2.3 Correlation analysis and evaluation based on the Pearson correlation coefficient method, the process is as follows: (1) Construct a feature matrix using the extracted features; (2) Calculate the Pearson correlation coefficient between each pair of features in the feature matrix to form a correlation coefficient matrix; (3) The data interval [-1, 1] is used to represent the Pearson correlation coefficient between each feature. The closer the Pearson coefficient is to -1, the stronger the negative correlation is. The closer it is to +1, the stronger the positive correlation is. The closer it is to 0, the weaker the correlation is.

6. The method for extracting vibration perception information features based on tunnel spatial correlation according to claim 1, characterized in that: The step 4 of obtaining the recognition accuracy of the signal features includes the following steps: According to steps one, two, and three, the features extracted and fused and corrected are used as key features. The recognition model is built using SVM. The key features are used as the sample set, and a part of them is randomly selected as the training set. The remaining part is trained as the test set. At the same time, a comparison group is established. The comparison group uses the unfused and corrected values ​​of the features. The SVM classifiers are trained for the test set and the comparison group respectively. The events recognized by the two trained SVM classifiers are recorded, and the recognition accuracy of each event and the overall event is calculated for comparison.

Images