Mixed-flow water turbine draft tube vortex strip characterization method based on pressure pulsation signals

Abstract

The invention discloses a mixed-flow water turbine draft tube vortex strip characterization method based on a pressure pulsation signal, and the method comprises the following steps: collecting the pressure pulsation signal of a draft tube during normal operation, preprocessing the signal, and predicting the pressure pulsation through a modal decomposition method; a sample pressure pulsation prediction error is calculated through a signal actual value and a prediction value, a threshold value is calculated in a self-adaptive mode through an SPOT algorithm, and early warning is conducted once when a monitoring index value exceeds the threshold value; chaos characteristic parameters of the primary early warning pressure pulsation signals are calculated; and a vortex strip early warning coefficient is defined by combining the chaos characteristic parameters of the pressure pulsation signals in the normal operation state, and when the early warning coefficient is larger than a threshold value, it is judged that an unstable vortex strip exists in the draft tube. The method is simple to operate, does not depend on a complete fault data set, and solves the problem that the flow state of the draft tube of the mixed-flow water turbine is difficult to accurately and quickly recognize in the prior art.

Description

Technical Field

[0001] The present invention belongs to the technical field of hydraulic turbines, and in particular, relates to a method for characterizing the draft tube vortex band of a Francis turbine based on pressure pulsation signals. Background Art

[0002] In the global energy pattern, hydropower, as a clean and renewable energy form, occupies a crucial position. The Francis turbine, with its high efficiency and wide applicability, has become the core power generation equipment in many hydropower stations. However, during the operation of the Francis turbine, the draft tube vortex band phenomenon is an extremely critical and complex problem, which has a profound impact on the stable operation, efficiency, and service life of the unit. The draft tube vortex band is an unstable flow structure caused by the special flow conditions of the water flow in the draft tube. When the vortex band is formed, it will oscillate at a certain frequency and amplitude, and this oscillation will generate strong pressure pulsations. These pressure pulsations will propagate along the water flow to the entire turbine flow passage, causing periodic forces on each component of the turbine. It may cause the runner blades to bear alternating stresses, accelerate the fatigue damage of the blade materials, shorten the service life of the blades, and even cause serious faults such as blade cracks. For the unit shafting, the pressure pulsations caused by the vortex band will generate fluctuations in axial force and radial force, causing the shafting to vibrate, destroying the stability of the shafting, increasing the wear between the shafting and the bearings, reducing the accuracy and reliability of the bearings, and thus affecting the operation smoothness of the entire unit.

[0003] With the continuous development of the hydropower industry, the requirements for the operation stability and reliability of Francis turbines are increasing day by day. In order to effectively solve the many problems brought by the draft tube vortex band, it has become the key to accurately monitor the state of the draft tube vortex band. By monitoring the state of the vortex band in real time and accurately, the internal flow conditions of the turbine can be grasped in a timely manner, the fault risks that may be caused by the vortex band can be warned in advance, and a strong basis can be provided for taking corresponding control measures. This is of great significance for ensuring the safe and stable operation of the turbine, improving the power generation efficiency, reducing the operation and maintenance costs, and promoting the sustainable development of the hydropower industry.

[0004] During the actual operation of a power station, the pressure pulsation at the draft tube part of the unit is usually monitored. However, only parameters such as its amplitude, waveform, and frequency domain can be obtained, and it is difficult to directly judge the internal flow state of the unit. Therefore, extracting the unit's flow state information through pressure pulsation signals is an ideal way for turbine state assessment. Currently, the commonly used methods for extracting the dynamic characteristics of the turbine draft tube mainly include extracting the characteristic entropy of the pressure pulsation signal. For example, the Chinese invention patent with the publication number CN103955601B and the publication date of February 15, 2017, judges the operation state of the turbine by analyzing the approximate entropy of the pressure pulsation signal. However, the inventor noticed that the single-scale characteristic entropy is calculated based on the logarithmic conditional probability mean of the new mode appearance when the signal sequence changes in dimension, which focuses more on the statistical analysis of the time series. During our research process, it was found that this method is often not ideal for characterizing the complexity of the chaotic system, especially in terms of its fractal structure and self-similarity. Moreover, this method still relies on the experience of certain operators and has a low degree of intelligence. Summary of the Invention

[0005] The purpose of the present invention is to provide a method for characterizing the vortex band in the draft tube of a Francis turbine based on pressure pulsation signals, which is simple to operate and has high discrimination accuracy, and solves the problem that it is difficult to accurately and quickly identify the flow state of the draft tube in the prior art.

[0006] To achieve this purpose, the method collects the pressure pulsation signals in the Francis turbine through high-precision pressure pulsation sensors installed at the cross-section of the draft cone tube of the Francis turbine, and accurately discriminates the development state of the vortex band in the draft tube of the unit.

[0007] The technical solution adopted by the present invention is a method for characterizing the vortex band in the draft tube of a Francis turbine based on pressure pulsation signals, which is specifically implemented according to the following steps:

[0008] Step 1: Collect the pressure pulsation signal P(t) of the draft tube during normal operation through high-precision pressure pulsation sensors installed at the cross-section of the draft cone tube of the Francis turbine, and preprocess the pressure pulsation signal.

[0009] Step 2: Use the modal decomposition method to decompose the pressure pulsation signal collected in Step 1 into multiple IMF components.

[0010] Step 3: Use the BiLSTM model to predict each IMF component. Superimpose the prediction results of each IMF component to obtain the prediction result, and retain the prediction model.

[0011] Step 4: Predict the pressure pulsation signal under unknown working conditions through the model. Calculate the prediction error of each sample pressure pulsation based on the actual value and the predicted value of the pressure pulsation signal, and use this error as the monitoring index.

[0012] Step 5: Adaptively calculate the threshold using the preset SPOT algorithm. If the monitored metric value exceeds the threshold, an alarm will be generated once, and it is considered that the unit is operating under unstable conditions at this time;

[0013] Step 6: Calculate the chaotic characteristic parameters of the early warning pressure pulsation signal, the attractor diffusion radius C R,a and the correlation dimension C D,a . Combine the chaotic characteristic parameters C R,0 and C D,0 of the pressure pulsation signal under normal operating conditions to define the vortex band early warning coefficient: When the early warning coefficients C R and C D are greater than the set threshold, it can be judged that there is an unstable vortex band in the draft tube.

[0014] The present invention further includes the following steps:

[0015] The data preprocessing process in Step 1 includes removing abnormal pressure pulsation data according to the Pauta criterion and normalizing the data.

[0016] The normalization process is as follows:

[0017]

[0018] where X i is the pressure pulsation signal data value, X min is the data minimum value, and X max is the data maximum value.

[0019] The specific implementation steps of Step 3 are as follows: Take the first 400 data points of the pressure pulsation data as input data, and the 401st data value as output data, and divide them into a sample. Subsequent samples are processed in the same way. After normalizing the training samples, input them into the BiLSTM deep learning model, and use ten-fold cross-validation and grid search to optimize the hyperparameters, the number of BiLSTM layers, the number of neurons, and the learning rate of the prediction model. When the number of training times reaches the set maximum number, save the model as the trained pressure pulsation prediction model.

[0020] The prediction error in Step 4 is calculated by the following formula:

[0021] D error = |P p - P t |

[0022] where D error represents the sample pressure pulsation prediction error; P p represents the sample predicted value; P t represents the sample actual value.

[0023] The specific steps of Step 5 are as follows: Let D be any number in the error sequence {D1, D2, …, D K}; Calculate the anomaly threshold of the first n anomaly scores, where n < K; Initialize the threshold T, obtain the peak point, and fit the parameters according to the Generalized Pareto Distribution (GPD):

[0024]

[0025] where F t (d) is the GPD fitting distribution function; t is the initial threshold; γ is the shape parameter obtained by GPD fitting, σ is the scale parameter obtained by GPD fitting, and d represents the data points of the data set;

[0026] Subsequently, adaptively calculate the anomaly threshold Z q :

[0027]

[0028] where t represents the initial threshold, i.e., the 98% quantile of the initialized data, is the estimated value of the scale parameter; γ is the estimated value of the shape parameter, q represents the probability that the sample exceeds the threshold; n is the sample size, and N t is the number of points where the peak is greater than t; For the error D n+i , i ≤ K - n: If d n+i > Z q , it is determined as an anomaly; If T < D n+i < Z q , refit the parameters according to GPD to obtain a new Z q ; If D n+i < t, no processing is performed.

[0029] Step 6, the attractor diffusion radius C R,a and the correlation dimension C D,a of the primary warning pressure pulsation signal; Specifically: Adaptively calculate the optimal delay time τ and the optimal embedding dimension m through the mutual information coefficient method and the false nearest neighbor method criterion; Map the pressure pulsation signal to a high-dimensional phase space through the optimal delay time τ and the optimal embedding dimension m, and calculate the high-order matrix Y i :

[0030] Y i =(Y i , Y i+τ …, Y i+(m-1)τ ), i = 1, 2, …, N - (m - 1)τ (4);

[0031] where N is the length of the signal pressure pulsation signal, and Y i is the i-th phase point in the matrix. Assume the high-order matrix Y iThe points in the phase space distribution each correspond to a coordinate (x1, y1, z1), (x2, y2, z2)…(x n , y n , z n ); then the high-order matrix Y of the pressure pulsation signal i The diffusion radius of the attractor is defined as the maximum value of the distances among all points:

[0032]

[0033] where is the center of all points.

[0034] The correlation dimension C D,a is defined as the slope of the correlation integral C(R) curve:

[0035]

[0036] R = exp(linspace(log(r min ), log(r max ), M))

[0037] where R is the similarity radius; M i (R) is the number of phase points within the similarity radius R; x i and x k are the i-th and k-th phase points; r min is the minimum value of the similarity radius, r max is the maximum value of the similarity radius, and M is the number of phase points in the linear region.

[0038] Combining the chaos characteristic parameters C R,0 and C D,0 of the pressure pulsation signal under the normal operation state, define the vortex band warning coefficient: When the warning coefficients C R and C D are greater than the set threshold values, it is considered that there is a significant spiral vortex band in the draft tube. Generally, the warning threshold values C R and C D are taken as 80 and 0.3 respectively.

[0039]

[0040] where C R,a and C D,a are the diffusion radius and correlation dimension of the pressure pulsation signal under the warning working condition; C R,0 and C D,0 are the diffusion radius and correlation dimension of the pressure pulsation signal under the normal operation state.

[0041] The method for characterizing the draft tube vortex band of a Francis turbine based on the pressure pulsation signal provided by the present invention has the following beneficial effects:

[0042] 1. The method is simple to operate and has high discrimination accuracy, solving the problem that it is difficult to accurately and quickly identify the flow state of the draft tube in the prior art. By continuously collecting the pressure pulsation signals of the draft tube during normal operation through high-precision pressure pulsation sensors pre-buried at the cross-section of the draft cone tube of the Francis turbine, there is no need to modify the unit. Then, the operating conditions of the unit are discriminated and a primary warning is given through the trained pressure pulsation prediction model based on deep learning. Since there is a lack of flow condition labels and incomplete fault types in actual engineering, the traditional classifier-based method cannot achieve the monitoring of abnormal flow conditions. However, the method proposed in the present invention adaptively discriminates abnormal flow states through the pressure pulsation prediction deviation and does not rely on a complete fault data set, and can be directly extended to the actual machine.

[0043] 2. When the unit is in an abnormal flow condition, the vortex band state in the draft tube is discriminated through the characteristics of the pressure pulsation signal attractor. When the internal flow of the Francis turbine is stable, the eddy current intensity in the draft tube is low and the pressure pulsation is small. At this time, the pressure pulsation in the draft tube is mainly dominated by the dynamic and static interference between the runner and the stay vanes, the diffusion radius of the pressure pulsation signal attractor is small, and the correlation dimension of the pressure pulsation is relatively large. As the turbine deviates from the optimal condition, that is, during the process of load reduction, the circumferential velocity of the water flow at the runner outlet begins to appear, and a vortex band begins to form in the draft tube. At this time, the diffusion radius of the pressure pulsation signal attractor begins to slowly increase, but the formation of the low-frequency vortex band causes the correlation dimension of the pressure pulsation signal to drop suddenly. As the degree of deviation of the turbine from the optimal condition increases, the intensity of the vortex band in the draft tube first increases and then gradually decreases, causing the diffusion radius of the pressure pulsation signal attractor to first increase and then decrease, while the correlation dimension shows a trend of first decreasing and then rising. The present invention quantifies the severity of the vortex band after discriminately identifying abnormal flow states through the evolution law of the characteristics of the pressure pulsation signal attractor of the turbine with the development of the vortex band. Description of the Drawings

[0044] The present invention will be further described below in conjunction with the drawings and embodiments:

[0045] Figure 1 is the flow schematic diagram of the method of the present invention;

[0046] Figure 2 is the trend diagram of the diffusion radius of the pressure pulsation signal attractor changing with the load in the method of the present invention;

[0047] Figure 3 is the trend diagram of the correlation dimension of the pressure pulsation signal attractor changing with the load in the method of the present invention;

[0048] Figure 4 is the modal decomposition result of the pressure pulsation signal under the normal operation condition (100% load) in Embodiment 2 of the present invention;

[0049] Figure 5 It is the prediction result of the training set of the pressure pulsation signal prediction model under the normal operation condition (100% load) in Embodiment 2 of the present invention;

[0050] Figure 6 It is the prediction result of the test set of the pressure pulsation signal prediction model under the normal operation condition (100% load) in Embodiment 2 of the present invention;

[0051] Figure 7 It is the primary warning result of the vortex band condition (70% load) in Embodiment 2 of the present invention.

[0052] Figure 8 It is the curve graph of the pressure pulsation signal delay time under the vortex band condition (70% load) in Embodiment 2 of the present invention;

[0053] Figure 9 It is the curve graph of the embedding dimension of the pressure pulsation signal under the vortex band condition (70% load) in Embodiment 2 of the present invention;

[0054] Figure 10 It is the prediction result of the pressure pulsation signal under the vortex band condition (50% load) in Embodiment 3 of the present invention;

[0055] Figure 11 It is the primary warning result of the vortex band condition (30% load) in Embodiment 4 of the present invention. Specific embodiments

[0056] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0057] Embodiment 1

[0058] A method for characterizing the vortex band in the draft tube of a Francis turbine based on pressure pulsation signals, based on the evolution law of the attractor characteristics of the turbine pressure pulsation signals with the development of the vortex band, quantifies the severity of the vortex band by collecting the pressure pulsation signals of the turbine unit after discriminating the abnormal flow regime in the primary warning. As Figure 1 shown, it is implemented by the following steps:

[0059] Step 1: Collect the pressure pulsation signals during the operation of the turbine through the pressure pulsation sensors pre-buried on the draft cone tube of the turbine. The pressure pulsation measurement point is located at 0.6 times the diameter from the rotation center of the runner on the wall of the draft tube. Remove the abnormal pressure pulsation data according to the Pauta criterion and normalize the data.

[0060] The normalization process is as follows:

[0061]

[0062] Where X i is the data value of the pressure pulsation signal, X min is the minimum data value, X max is the maximum data value.

[0063] Step 2: Use the modal decomposition method to decompose the pressure pulsation signal collected in Step 1 into multiple IMF components;

[0064] Step 2.1: Add white noise of a certain magnitude to the time series P(t) of the pressure pulsation signal after preprocessing to form a signal sequence to be decomposed:

[0065] P i (t) = P(t) + β o E1(w i (t))

[0066] In the formula, P i (t) is the signal to be decomposed with added noise; i = 1, 2, 3…, i represents the number of times of adding white noise; β0 is the scale coefficient of adding white noise, w i (t) is white noise with zero mean and unit variance; E1(w i (t) is the first EMD component of w i (t);

[0067] Step 2.2: Perform EMD decomposition on the signal P i (t) to be decomposed with added noise to obtain the first residual component r1(t) and the IMF component IMF1(t)

[0068] r1(t) = M(P i (t))

[0069] IMF1(t) = P i (t) - r1(t)

[0070] In the formula, the operator M is the local mean of the signal, that is, the average after summing the signal amplitudes.

[0071] Step 2.3: Calculate the kth residual component r k (t) and the kth IMF component IMF k (t)

[0072] r k (t) = M(r k-1 (t) + β k-1 E k (w i (t)))

[0073] IMFk r(t) = r k-1 r(t) - r k r(t)

[0074] where k = 2, 3, 4…, N; β k-1 represents the scale coefficient of the noise added in the k-th decomposition;

[0075] Step 2.4: Repeat the above two steps until the residual component is r(t), a monotonic function, and calculate k IMF components of the pressure pulsation signal. k

[0076] Step 3: Use the BiLSTM model to predict each IMF component. Take the first 400 data points of the pressure pulsation data as the input data and the 401st data value as the output data, and divide them into a sample. Subsequent samples are processed in the same way. After normalizing the training samples, input them into the BiLSTM deep learning model. Use ten-fold cross-validation and grid search to optimize the hyperparameters, the number of BiLSTM layers, the number of neurons, and the learning rate of the prediction model. When the number of training times reaches the set maximum number, save the model as the trained pressure pulsation prediction model. Superimpose the prediction results of each IMF component to obtain the prediction result and retain the prediction model.

[0077] Step 4: Predict the pressure pulsation signal under unknown working conditions through the model. Calculate the prediction error of the pressure pulsation for each sample based on the actual value and the predicted value of the pressure pulsation signal, and use this error as the monitoring index. Calculate through the following formula:

[0078] D error = |P p - P t |

[0079] where D error represents the prediction error of the pressure pulsation for the sample; P p represents the predicted value of the sample; P t represents the actual value of the sample.

[0080] Step 5: Use the preset SPOT algorithm to adaptively calculate the threshold. When the value of the monitoring index exceeds the threshold, an alarm will be generated once, and it is considered that the unit is operating under unstable working conditions at this time.

[0081] Step 5.1: Let D be any number in the error sequence {D1, D2, …, D K}; Calculate the anomaly threshold of the first n anomaly scores, where n < K; Initialize the threshold T, obtain the peak point, and fit the parameters according to the generalized Pareto distribution (GPD) of

[0082]

[0083] ​Among them, F t (d) is the GPD fitting distribution function; t is the initial threshold; γ is the shape parameter obtained by GPD fitting, σ is the scale parameter obtained by GPD fitting, d represents the data points of the data set;

[0084] Step 5.2. Subsequently, adaptively calculate the anomaly threshold Z q :

[0085]

[0086] Among them, t represents the initial threshold, that is, the 98% quantile of the initialization data, is the estimated value of the scale parameter; γ is the estimated value of the shape parameter. q represents the probability that the sample exceeds the threshold; n is the number of samples, and N t is the number of points whose peak value is greater than t. For the error D n+i , i ≤ K - n: If d n+i > Z q , it is determined as an anomaly; if T < D n+i < Z q , re - fit the parameters according to GPD to obtain a new Z q ; if D n+i < t, no processing is performed.

[0087] Step 6. Calculate the chaotic characteristic parameters of the primary warning pressure pulsation signal, the attractor diffusion radius C R,a and the correlation dimension C D,a . Combine the chaotic characteristic parameters C R,0 and C D,0 of the pressure pulsation signal under the normal operation state to define the vortex - band warning coefficient: When the warning coefficients C R and C D are greater than the set threshold, it can be judged that there is an unstable vortex - band in the draft tube.

[0088] Step 6.1. Adopt the mutual - information coefficient method to calculate the mutual information between the pressure pulsation signal and the corresponding delayed signal in the pressure pulsation signal matrix, and adaptively select the optimal delay time τ according to the first local minimum of the mutual information;

[0089]

[0090] In the formula, I is the mutual information, N is the length of the pressure pulsation signal, T is the set delay time between [1, 50]; x i is the pressure pulsation signal at the i - th point; p(x i ) is the probability density of the pressure pulsation signal; p(x i + T) is the probability density of the delayed signal; p(x i , x ip(x,y; +T) is the joint probability density of the two signals; starting from the minimum value of the delay time T, calculate the mutual information, increase the value of T, until the first minimum value of I(T) appears, at this time T is the optimal delay time τ of a certain measurement point;

[0091] Step 6.2, after determining the optimal delay time, adaptively calculate the optimal embedding dimension according to the false nearest neighbor method criterion;

[0092]

[0093] In the formula, and represent the distances between two points in the (d + 1)-dimensional and d-dimensional phase spaces respectively, and Rtol is the threshold value in the range of [10, 50]; is the j-th phase point vector in the reconstructed phase space, and the superscript r represents the phase point coordinates in the reconstructed phase space; is the nearest neighboring phase point of the phase point ; L is the number of phase points after the reconstruction of the pressure pulsation signal; starting from the minimum value of the dimension d, calculate the proportion of false nearest neighbor points, gradually increase the value of d until all false nearest neighbor points disappear, at this time d is the optimal embedding dimension m of a certain measurement point;

[0094] Step 6.3, map the pressure pulsation signal to a high-dimensional phase space through the optimal delay time τ and the optimal embedding dimension m, and calculate the high-order matrix Y i :

[0095] Y i =(Y i , Y i+τ …, Y i+(m-1)τ ), i = 1, 2, …, N - (m - 1)τ (4);

[0096] where N is the length of the signal pressure pulsation signal, and Y i is the i-th phase point in the matrix. Assume that each point in the phase space distribution of the high-order matrix Y i corresponds to a coordinate (x1, y1, z1), (x2, y2, z2)…(x n , y n , z n ). Then the diffusion radius of the attractor of the pressure pulsation signal's high-order matrix Y i is defined as the maximum value of the distances among all points:

[0097]

[0098] where is the center of all points.

[0099] Step 6.4, the correlation dimension C D,a is defined as the slope of the correlation integral C(R) curve:

[0100]

[0101] R = exp(linspace(log(r min ), log(r max ), M))

[0102] where R is the similarity radius; M i (R) is the number of phase points within the similarity radius R; x i and x k are the i-th and k-th phase points; r min is the minimum value of the similarity radius, r max is the maximum value of the similarity radius, and M is the number of phase points in the linear region.

[0103] Combined with the chaotic characteristic parameters C R,0 and C D,0 of the pressure pulsation signal under normal operating conditions, a vortex band warning coefficient is defined: when the warning coefficients C R and C D are greater than the set threshold, it is considered that there is a significant spiral vortex band in the draft tube.

[0104]

[0105] where C R,a and C D,a are the diffusion radius and correlation dimension of the pressure pulsation signal under the warning condition; C R,0 and C D,0 are the diffusion radius and correlation dimension of the pressure pulsation signal under normal operating conditions. Generally, the warning thresholds C R and C D are taken as 80 and 0.3 respectively.

[0106] In this embodiment, the abnormal flow regime is adaptively discriminated through the pressure pulsation prediction deviation, without relying on a complete fault data set. Based on the formation principle of the turbine vortex band, after the abnormal flow regime is adaptively identified through the pressure pulsation prediction model, by analyzing the correlation characteristics between the evolution state of the turbine vortex band and the pressure pulsation signal, more effective cavitation characteristic indexes, the attractor diffusion radius and the correlation dimension, are proposed based on the chaos theory; the disadvantages of being vulnerable to noise in the traditional time-frequency domain feature extraction are effectively avoided, the accuracy of quantifying the turbine vortex band is improved, the calculation method is simple and the accuracy of vortex band characterization is improved.

[0107] The vortex band characterization method of the present invention is based on the formation principle of the vortex band in the water turbine and its influence law on the pressure pulsation signal. Therefore, the pressure pulsation signal in the draft tube is collected by monitoring the pressure pulsation sensor, and after the adverse flow regime is adaptively identified by the pressure pulsation prediction model and the SPOT algorithm, the variation law of the diffusion radius and correlation dimension of the pressure pulsation signal attractor with the evolution of the vortex band is used to quantify the degree of the vortex band in the draft tube of the water turbine. When the unit is in an adverse flow condition, the state of the vortex band in the draft tube is judged by the characteristics of the pressure pulsation signal attractor. When the internal flow of the Francis turbine is stable, the eddy current intensity in the draft tube is low and the pressure pulsation is small. At this time, the pressure pulsation in the draft tube is mainly dominated by the dynamic and static interference between the runner and the stay vanes, the diffusion radius of the pressure pulsation signal attractor is small, and the correlation dimension of the pressure pulsation is relatively large, as shown in Figure 2 and Figure 3 in the stable operation area shown. As the water turbine deviates from the optimal operating condition, that is, during the process of reducing the load, the circumferential velocity of the water flow at the runner outlet begins to appear, and a vortex band begins to form in the draft tube. At this time, the diffusion radius of the pressure pulsation signal attractor begins to increase slowly, but the formation of the low-frequency vortex band causes the correlation dimension of the pressure pulsation signal to drop suddenly; as the degree of deviation of the water turbine from the optimal operating condition increases, the vortex band intensity in the draft tube first increases and then gradually decreases, causing the diffusion radius of the pressure pulsation signal attractor to first increase and then decrease, while the correlation dimension develops in a trend of first decreasing and then increasing (as shown in Figure 2 、 3 ). The present invention quantifies the severity of the vortex band by the evolution law of the characteristics of the pressure pulsation signal attractor of the water turbine with the development of the vortex band after the adverse flow regime is discriminated in the first warning.

[0108] Example 2

[0109] A method for characterizing the vortex band in the draft tube of a Francis turbine based on the pressure pulsation signal is implemented by the following steps:

[0110] Step 1: The pressure pulsation signal during the operation of the water turbine is measured by the pressure pulsation sensor pre-buried on the draft cone tube of the water turbine. The pressure pulsation measurement point is located at 0.6 times the diameter from the rotation center of the runner on the draft tube wall. For subsequent application standardization, the pressure pulsation value is converted into a dimensionless pressure pulsation C p , and the method is as follows:

[0111]

[0112] where p represents the collected pressure pulsation data, represents the average value of the pressure pulsation data, ρ represents the density of water, g represents the acceleration due to gravity, and H represents the head value of this condition.

[0113] Taking a Francis turbine prototype with 9 runner blades, 20 stay vanes and 20 guide vanes as an example, the pressure pulsation sensor used is the M112A22 series high-precision pressure pulsation sensor. The sampling frequency is set to 10.24 kHz, the sensitivity is 14.5 mV / kPa, and the measurement error is less than 0.5% FS. Start the Francis turbine test bench to make the unit operate at full load under normal conditions, gradually reduce the load, and collect the pressure pulsation signals at the draft tube.

[0114] In this example, the operating parameters of the turbine under the experimental conditions are shown in Table 1:

[0115] Table 1 shows the operating parameters of the turbine under the conditions of Example 2

[0116] Experimental conditions Rotational speed / (r / min) <![CDATA[Flow rate / (m 3 / s)]]> 100% load condition 1188.11 0.314 70% load condition 1188.11 0.22

[0117] After the data acquisition is completed, remove the abnormal pressure pulsation data and normalize the data according to the 3σ criterion. The normalization process is as follows:

[0118]

[0119] where X i is the data value of the pressure pulsation signal, X min is the minimum value of the data, and X max is the maximum value of the data.

[0120] Step 2: Use the modal decomposition method to decompose the pressure pulsation signal collected in Step 1 into multiple IMF components;

[0121] Step 2.1: Add white noise of a certain magnitude to the time series P(t) of the preprocessed pressure pulsation signal to form a signal series to be decomposed:

[0122] P i (t) = P(t) + β o E1(w i (t))

[0123] where P i (t) is the signal to be decomposed with added noise; i = 1, 2, 3…, i represents the number of times of adding white noise; β0 is the scale coefficient of the added white noise, w i (t) is white noise with zero mean and unit variance; E1(w i (t) is the first EMD component of w i (t);

[0124] Step 2.2: Perform EMD decomposition on the signal P i (t) with added noise to obtain the first residual component r1(t) and the IMF component IMF1(t)

[0125] r1(t) = M(P i (t))

[0126] IMF1(t) = P i (t) - r1(t)

[0127] In the formula, the operator M is the local mean of the signal, that is, the average is obtained after summing the signal amplitudes.

[0128] Step 2.3: Calculate the k-th residual component r k (t) and the k-th IMF component IMF k (t)

[0129] r k (t) = M(r k-1 (t) + β k-1 E k (w i (t)))

[0130] IMF k (t) = r k-1 (t) - r k (t)

[0131] In the formula, k = 2, 3, 4…, N; β k-1 represents the scale coefficient of the noise added in the k-th decomposition;

[0132] Step 2.4: Repeat the above two steps until the residual component is r k (t) is a monotonic function, and k IMF components of the pressure pulsation signal are calculated.

[0133] The modal decomposition result of the pressure pulsation signal under the normal operating condition (100% load) is as Figure 4 shown. It can be found that 9 IMF components are obtained through adaptive modal decomposition. The frequencies of the components from IMF1 to IMF9 gradually decrease, indicating that the characteristics in the original signal are decoupled after modal decomposition, avoiding the possible problem of feature autocorrelation in the process of directly predicting the original signal. Targeted prediction of each component can improve the accuracy of pressure pulsation prediction.

[0134] Step 3: Use the BiLSTM model to predict each IMF component. Take the first 400 data points of the pressure pulsation data as input data and the 401st data value as output data, and divide them into a sample. Subsequent samples are processed in the same way. After normalizing the training samples, input them into the BiLSTM deep learning model. Use ten-fold cross-validation and grid search to optimize the hyperparameters of the prediction model, the number of BiLSTM layers, the number of neurons, and the learning rate. When the number of training times reaches the set maximum number, save the model as the trained pressure pulsation prediction model. Superimpose the prediction results of each IMF component to obtain the prediction result.

[0135] In this embodiment, a total of 2740 groups of pressure pulsation samples under normal operating conditions are obtained, among which 2500 groups are divided into the training set and 240 groups are the test set. In this example, the initial learning rate is 0.005, and the adam gradient descent algorithm is used. The prediction results of the training set and the test set obtained by training the pressure pulsation signal prediction model 400 times are as Figure 5 and Figure 6 shown. It can be found that whether it is the training set or the prediction set, the change trend of the predicted value of the pressure pulsation is consistent with the actual value, the numerical values are basically equal, and the prediction error is small, indicating that the trained pressure pulsation prediction model has good prediction performance and strong generalization performance. The evaluation indexes of the prediction results are shown in Table 2. It can be found that the determination coefficients are all greater than 0.99, and the mean absolute errors are also small, further indicating the effectiveness and accuracy of the prediction model.

[0136] Table 2 shows the prediction result indexes of the water turbine pressure pulsation signal under the working conditions of Embodiment 2

[0137] Sample set <![CDATA[Coefficient of determination R 2 > Mean absolute error Mean bias error Training set 0.997 3.569e-07 -2.747e-07 Test set 0.992 5.012e-07 -3.343e-07

[0138] Step 4: Retain the prediction model. In this embodiment, a total of 3108 pressure pulsation samples are collected under the condition of reducing the load to 70% for subsequent prediction. Predict the pressure pulsation signal during the load reduction process through the model, and calculate the pressure pulsation prediction error of each sample based on the actual value and the predicted value of the pressure pulsation signal. Take this error as the monitoring index. Calculate through the following formula:

[0139] D error =|P p -P t |

[0140] where D error represents the pressure pulsation prediction error of the sample; P [ represents the predicted value of the sample; P t represents the actual value of the sample.

[0141] Step 5: Use the preset SPOT algorithm to adaptively calculate the threshold. When the monitored index value exceeds the threshold, an alarm will be generated once, and it is considered that the unit is operating under unstable conditions at this time.

[0142] Step 5.1: Let D be any number in the error sequence {D1, D2, …, D K}; Calculate the anomaly threshold of the first n anomaly scores, where n < K; Initialize the threshold T, obtain the peak point, and obtain the parameters according to the Generalized Pareto Distribution (GPD) fitting:

[0143]

[0144] where F t (d) is the GPD fitting distribution function; t is the initial threshold; γ is the shape parameter obtained by GPD fitting, σ is the scale parameter obtained by GPD fitting, and d represents the data points of the data set;

[0145] Step 5.2: Subsequently, adaptively calculate the anomaly threshold Z q :

[0146]

[0147] where t represents the initial threshold, that is, the 98% quantile of the initialization data, is the estimated value of the scale parameter; γ is the estimated value of the shape parameter, q represents the probability that the sample exceeds the threshold; n is the number of samples, and N t is the number of points where the peak value is greater than t. For the error D n+i , i ≤ K - n: If d n+i > Z q , it is determined as an anomaly; if T < D n+i < Z q , re - fit the parameters according to GPD to obtain a new Z q ; if D n+i < t, no processing is performed.

[0148] In this embodiment, the results of the first warning under the condition of reducing the load to 70% are as Figure 7 shown. It can be found that the prediction error is almost 0 in the initial process of load reduction. As the load further decreases, after the vortex band is generated, the prediction error increases rapidly by nearly five times, exceeding the calculated adaptive threshold, and a warning is issued, indicating that the unit is operating under poor conditions at this time.

[0149] Step 6: Calculate the chaotic characteristic parameters of the first warning pressure pulsation signal, the attractor diffusion radius C R,a and the correlation dimension C D,a . Combine the chaotic characteristic parameters C R,0 and CD,0 Define the vortex band warning coefficient: When the warning coefficients C R and C D are greater than the set threshold value, it can be judged that there is an unstable vortex band in the draft tube.

[0150] Step 6.1: Adopt the mutual information coefficient method to calculate the mutual information between the pressure pulsation signal and the corresponding delayed signal in the pressure pulsation signal matrix, and adaptively select the optimal delay time τ according to the first local minimum of the mutual information;

[0151]

[0152] In the formula, I is the mutual information, N is the length of the pressure pulsation signal, T is the set delay time between [1, 50]; x i is the pressure pulsation signal at the i-th point; p(x i ) is the probability density of the pressure pulsation signal; p(x i +T) is the probability density of the delayed signal; p(x i , x i +T) is the joint probability density of the two signals; Calculate the mutual information starting from the minimum value of the delay time T, increase the value of T until the first minimum value of I(T) appears, and at this time T is the optimal delay time τ of a certain measurement point;

[0153] Step 6.2: After determining the optimal delay time, adaptively calculate the optimal embedding dimension according to the false nearest neighbor method criterion;

[0154]

[0155] In the formula, and represent the distances between two points in the (d + 1)-th and d-th dimensional phase spaces respectively, Rtol is the threshold value of [10, 50], is the j-th phase point vector in the reconstructed phase space, the superscript r represents the phase point coordinates in the reconstructed phase space, is the nearest neighbor phase point of the phase point ; L is the number of phase points after the pressure pulsation signal is reconstructed; Calculate the proportion of false nearest neighbor points starting from the minimum value of the dimension d, gradually increase the value of d until all false nearest neighbor points disappear, and at this time d is the optimal embedding dimension m of a certain measurement point;

[0156] Table 3 shows the optimal phase space reconstruction parameters of the turbine pressure pulsation signal under different working conditions in Example 2

[0157] Experimental conditions Embedding dimension m Delay time τ 100% load condition 3 14 70% load condition 5 60

[0158] The delay time curve and the embedding dimension curve of the pressure pulsation signal under the 70% load condition are as shown in Figure 8 and Figure 9as shown

[0159] As Figure 8 shown, in this example, the first minimum point of the mutual information of the draft tube pressure pulsation signal in the vortex band state appears when the delay time reaches 60. Therefore, the optimal delay time is 60. As Figure 9 shown, in this example, the false nearest neighbor rate of the mutual information of the draft tube pressure pulsation signal in the vortex band state decreases to 0 when the embedding dimension reaches 5, and the false nearest neighbor rate no longer changes when the embedding dimension continues to increase. Therefore, the optimal embedding dimension is 5. The optimal phase space reconstruction parameters adaptively determined for the pressure pulsation sequences under different working conditions are shown in Table 3.

[0160] Step 6.3: Map the pressure pulsation signal to a high-dimensional phase space through the optimal delay time τ and the optimal embedding dimension m, and calculate the high-order matrix Y i :

[0161] Y i =(Y i , Y i+τ …, Y i+(m-1)τ ), i = 1, 2, …, N-(m-1)τ (4);

[0162] where N is the length of the signal pressure pulsation signal, and Y i is the i-th phase point in the matrix. Assume that the points distributed in the high-order matrix Y i phase space all correspond to a coordinate (x1, y1, z1), (x2, y2, z2)…(x n , y n , z n ). Then the diffusion radius of the attractor of the pressure pulsation signal's high-order matrix Y i is defined as the maximum value of the distances among all points:

[0163]

[0164] where is the center of all points. In this embodiment, the calculated diffusion radii of the attractor under the normal operating condition and the early warning condition are 2.84e-07 and 3.21e-05 respectively.

[0165] Step 6.4: The correlation dimension C D,a is defined as the slope of the correlation integral C(R) curve:

[0166]

[0167] R = exp(linspace(log(r min ), log(r max ), M))

[0168] where R is the similarity radius; Mi (R) is the number of phase points within the similarity radius R; x i and x k are the i-th and k-th phase points; r min is the minimum value of the similarity radius, r max is the maximum value of the similarity radius, and M is the number of phase points in the linear region. In this embodiment, the calculated correlation dimensions of the attractors under normal operating conditions and warning conditions are 2.633 and 1.254 respectively.

[0169] Combined with the chaotic characteristic parameters C R,0 and C D,0 of the pressure pulsation signal under normal operating conditions, define the vortex band warning coefficient: when the warning coefficients C R and C D are greater than the set threshold, it is considered that there is a significant spiral vortex band in the draft tube.

[0170]

[0171] where C R,a and C D,a are the diffusion radius and correlation dimension of the pressure pulsation signal under warning conditions; C R,0 and C D,0 are the diffusion radius and correlation dimension of the pressure pulsation signal under normal operating conditions. In this embodiment, the calculated warning coefficients C R and C D are 112 and 0.524 respectively, both exceeding the set warning threshold, indicating that it is in the vortex band condition at this time. Due to the lack of flow condition labels and incomplete fault types in actual engineering, the traditional classifier-based method cannot achieve the monitoring of poor flow conditions. The above embodiment discriminates poor flow states through the adaptive prediction deviation of pressure pulsation and does not rely on a complete fault data set. This illustrates the effectiveness of the present invention in characterizing the vortex band of a Francis turbine in this case.

[0172] Embodiment 3

[0173] A method for characterizing the vortex band in the draft tube of a Francis turbine based on pressure pulsation signals is implemented by the following steps:

[0174] Step 1, collect the pressure pulsation signal during the operation of the turbine through the pressure pulsation sensor pre-buried on the draft cone of the turbine. The pressure pulsation measurement point is located at 0.6 times the diameter from the center of rotation of the runner on the draft tube wall. Convert the pressure pulsation value to a dimensionless pressure pulsation C for subsequent application standardization p , the method is as follows:

[0175]

[0176] where p represents the collected pressure pulsation data, represents the average value of the pressure pulsation data, ρ represents the density of water, g represents the acceleration due to gravity, and H represents the head value under this operating condition.

[0177] In this example, the operating parameters of the water turbine under the experimental conditions are shown in Table 4:

[0178] Table 4 shows the operating parameters of the water turbine under the conditions of Example 3

[0179] Experimental conditions Rotational speed / (r / min) <![CDATA[Flow rate / (m 3 / s)]]> 50% load condition 1188.11 0.157

[0180] After the data acquisition is completed, the abnormal pressure pulsation data is removed according to the 3σ criterion and the data is normalized. The normalization process is as follows:

[0181]

[0182] where X i is the data value of the pressure pulsation signal, X min is the minimum value of the data, and X max is the maximum value of the data.

[0183] Step 2: Use the modal decomposition method to decompose the pressure pulsation signal collected in Step 1 into multiple IMF components;

[0184] Step 2.1: Add white noise of a certain magnitude to the time series P(t) of the pressure pulsation signal after preprocessing to form a signal sequence to be decomposed:

[0185] P i (t) = P(t) + β o E1(w i (t))

[0186] In the formula, P i (t) is the signal to be decomposed with added noise; i = 1, 2, 3…, i represents the number of times of adding white noise; β0 is the scale coefficient of adding white noise, w i (t) is white noise with zero mean and unit variance; E1(w i (t) is the first EMD component of w i (t);

[0187] Step 2.2: Perform EMD decomposition on the signal P i (t) to be decomposed with added noise to obtain the first residual component r1(t) and the IMF component IMF1(t)

[0188] r1(t) = M(P i (t))

[0189] IMF1(t) = P i (t) - r1(t)

[0190] The operator M in the formula is the local mean of the signal, that is, the signal amplitude is summed and then averaged.

[0191] Step 2.3: Calculate the kth residual component r k (t) and the kth IMF component IMF k (t)

[0192] r k (t) = M (r k-1 (t)+β k-1 E k (w i (t)))

[0193] IMF k (t) = r k-1 (t)-r k (t)

[0194] Where k = 2, 3, 4 +, N; β k-1 represents the scale coefficient of the noise added to the k-th decomposition;

[0195] Step 2.4: Repeat the above two steps until the residual component is r k (t) The monotonic function is calculated until k IMF components of the pressure pulsation signal are obtained.

[0196] Step 3: Use the BiLSTM model to predict each IMF component. The first 400 data points of the pressure pulsation data are used as input data, and the 401st data value is used as output data. The two are divided into one sample. The subsequent samples are analogous. After normalizing the training samples, input them into the BiLSTM deep learning model. Ten-fold cross validation and grid search are used to optimize the hyperparameters, number of BiLSTM layers, number of neurons, and learning rate of the prediction model. When the number of training times reaches the set maximum number, save the changed model as the trained pressure pulsation prediction model. Superimpose the prediction results of each IMF component to obtain the prediction result.

[0197] Step 4: retain the prediction model. In this embodiment, a total of 3171 pressure pulsation samples are collected under the load condition of 50% for subsequent prediction. The pressure pulsation signal in the process of load reduction is predicted by the model. The prediction result of the pressure pulsation signal under the load condition of 50% is as follows: Figure 10 As shown, it can be clearly found that the pressure pulsation prediction error is large under this working condition. Based on the actual value and predicted value of the pressure pulsation signal, the pressure pulsation prediction error of each sample is calculated and used as a monitoring indicator. Calculated by the following formula:

[0198] D error =|P p -P t |

[0199] Among them, D error represents the prediction error of the sample pressure pulsation; P p represents the sample predicted value; P t represents the actual value of the sample.

[0200] Step 5: Use the preset SPOT algorithm to adaptively calculate the threshold. When the monitored index value exceeds the threshold, an alarm will be generated, and it is considered that the unit is operating under unstable conditions at this time.

[0201] Step 5.1: Let D be any number in the error sequence {D1, D2, …, D K}; Calculate the anomaly threshold of the first n anomaly scores, where n < K; Initialize the threshold T, obtain the peak point, and obtain the parameters according to the generalized Pareto distribution (GPD) fitting:

[0202]

[0203] Among them, F t (d) is the GPD fitting distribution function; t is the initial threshold; γ is the shape parameter obtained by GPD fitting, σ is the scale parameter obtained by GPD fitting, and d represents the data points of the data set;

[0204] Step 5.2: Subsequently, adaptively calculate the anomaly threshold Z q :

[0205]

[0206] Among them, t represents the initial threshold, that is, the 98% quantile of the initialized data, is the estimated value of the scale parameter; γ is the estimated value of the shape parameter, q represents the probability that the sample exceeds the threshold; n is the number of samples, and N t is the number of points whose peak value is greater than t. For the error D n+i , i ≤ K - n: If d n+i > Z q , it is determined as an anomaly; if T < D n+i < Z q , re-fit the parameters according to GPD to obtain a new Z q ; if D n+i < t, no processing is performed.

[0207] Step 6: Calculate the chaotic characteristic parameters of the early warning pressure pulsation signal, the attractor diffusion radius C R,a and the correlation dimension C D,a . Combine the chaotic characteristic parameters C R,0 and C D,0 of the pressure pulsation signal under normal operating conditions to define the vortex band early warning coefficient: When the early warning coefficient C R and C DWhen it is greater than the set threshold value, it can be determined that there is an unstable vortex band in the draft tube.

[0208] Step 6.1: Use the mutual information coefficient method to calculate the mutual information between the pressure pulsation signal and the corresponding delayed signal in the pressure pulsation signal matrix, and adaptively select the optimal delay time τ according to the first local minimum of the mutual information.

[0209]

[0210] In the formula, I is the mutual information, N is the length of the pressure pulsation signal, T is the set delay time between [1, 50]; x i is the pressure pulsation signal at the i-th point; p(x i ) is the probability density of the pressure pulsation signal; p(x i +T) is the probability density of the delayed signal; p(x i , x i +T) is the joint probability density of the two signals; Calculate the mutual information starting from the minimum value of the delay time T, increase the value of T until the first minimum value of I(T) appears, and at this time T is the optimal delay time τ of a certain measurement point.

[0211] Step 6.2: After determining the optimal delay time, adaptively calculate the optimal embedding dimension according to the false nearest neighbor method criterion.

[0212]

[0213] In the formula, and represent the distances between two points in the (d + 1)-th and d-th dimensional phase spaces respectively, Rtol is the threshold value of [10, 50], is the j-th phase point vector in the reconstructed phase space, the superscript r represents the phase point coordinates in the reconstructed phase space, is the nearest neighboring phase point of the phase point ; L is the number of phase points after the reconstruction of the pressure pulsation signal; Calculate the proportion of false nearest neighbor points starting from the minimum value of the dimension d, gradually increase the value of d until all false nearest neighbor points disappear, and at this time d is the optimal embedding dimension m of a certain measurement point.

[0214] Table 5 shows the optimal phase space reconstruction parameters of the turbine pressure pulsation signal under different working conditions in Embodiment 3

[0215] Experimental conditions Embedding dimension m Delay time τ 50% load condition 3 91

[0216] The optimal phase space reconstruction parameters adaptively determined for the pressure pulsation sequence under the 50% load condition are shown in Table 5.

[0217] Step 6.3: Map the pressure pulsation signal to the high-dimensional phase space through the optimal delay time τ and the optimal embedding dimension m, and calculate the high-order matrix Yi :

[0218] Y i = (Y i , Y i+τ …, Y i+(m-1)τ ), i = 1, 2, …, N-(m-1)τ(4);

[0219] where N is the length of the signal pressure pulsation signal, and Y i is the i-th phase point in the matrix. Assume that all the points in the high-order matrix Y i phase space distribution correspond to a coordinate (x1, y1, z1), (x2, y2, z2)…(x n , y n , z n ). Then the high-order matrix Y i of the pressure pulsation signal has the attractor diffusion radius defined as the maximum value of the distances among all the points:

[0220]

[0221] where is the center of all the points. In this embodiment, the calculated attractor diffusion radius under the warning condition is 0.0073.

[0222] Step 6.4, Correlation dimension C D,a is defined as the slope of the correlation integral C(R) curve:

[0223]

[0224] R = exp(linspace(log(r min ), log(r max ), M))

[0225] where R is the similarity radius; M i (R) is the number of phase points within the similarity radius R; x i and x k are the i-th and k-th phase points; r min is the minimum value of the similarity radius, r max is the maximum value of the similarity radius, and M is the number of phase points in the linear region. In this embodiment, the calculated attractor correlation dimension under the warning condition is 1.107.

[0226] Combining the chaotic characteristic parameters C R,0 and C D,0 of the pressure pulsation signal in the normal operation state, define the vortex band warning coefficient: When the warning coefficients C R and C D are greater than the set threshold, it is considered that there is a significant spiral vortex band in the draft tube.

[0227]

[0228] where C R,a and C D,a are the diffusion radius and correlation dimension of the pressure pulsation signal under the warning condition; C R,0 and C D,0 are the diffusion radius and correlation dimension of the pressure pulsation signal under the normal operation state. In this embodiment, the calculated warning coefficients C R and C D are 2.57e4 and 0.579 respectively, both exceeding the set warning threshold, indicating that it is in the vortex band condition at this time. The above embodiment discriminates the abnormal flow pattern through the adaptive discrimination of the pressure pulsation prediction deviation and does not rely on a complete fault data set. It illustrates the effectiveness of the present invention for the vortex band characterization of a Francis turbine.

[0229] Embodiment 4

[0230] A method for characterizing the vortex band in the draft tube of a Francis turbine based on pressure pulsation signals is implemented according to the following steps:

[0231] Step 1: Measure the pressure pulsation signal during the operation of the turbine through the pressure pulsation sensor pre-buried on the draft cone tube of the turbine. The pressure pulsation measurement point is located at 0.6 times the diameter from the rotation center of the runner on the draft tube wall. To convert the pressure pulsation value into a dimensionless pressure pulsation C p , the method is as follows:

[0232]

[0233] where p represents the collected pressure pulsation data, represents the average value of the pressure pulsation data, ρ represents the density of water, g represents the acceleration due to gravity, and H represents the head value of this condition.

[0234] In this example, the operating parameters of the turbine under the experimental conditions are shown in Table 6:

[0235] Table 6 shows the operating parameters of the turbine under the conditions of Embodiment 4

[0236] Experimental conditions Rotational speed / (r / min) <![CDATA[Flow rate / (m 3 / s)]]> 30% load condition 1188.11 0.094

[0237] After the data acquisition is completed, the abnormal pressure pulsation data is removed according to the Pauta criterion and the data is normalized. The normalization process is as follows:

[0238]

[0239] where X i is the data value of the pressure pulsation signal, X min is the minimum value of the data, X maxis the data maximum value.

[0240] Step 2: Use the modal decomposition method to decompose the pressure pulsation signal collected in Step 1 into multiple IMF components;

[0241] Step 2.1: Add white noise of a certain magnitude to the time series P(t) of the preprocessed pressure pulsation signal to form a signal sequence to be decomposed:

[0242] P i (t) = P(t) + β o E1(w i (t))

[0243] Wherein, P i (t) is the signal to be decomposed with added noise; i = 1, 2, 3…, i represents the number of times of adding white noise; β0 is the scale coefficient of adding white noise, w i (t) is white noise with zero mean and unit variance; E1(w i (t) is the first EMD component of w i (t);

[0244] Step 2.2: Perform EMD decomposition on the signal P i (t) to be decomposed with added noise to obtain the first residual component r1(t) and the IMF component IMF1(t)

[0245] r1(t) = M(P i (t))

[0246] IMF1(t) = P i (t) - r1(t)

[0247] Wherein, the operator M is the local mean of the signal, that is, the average after summing the signal amplitudes.

[0248] Step 2.3: Calculate the kth residual component r k (t) and the kth IMF component IMF k (t)

[0249] r k (t) = M(r k-1 (t) + β k-1 E k (w i (t)))

[0250] IMF k (t) = r k-1 (t) - r k (t)

[0251] Wherein k = 2, 3, 4…, N; β k-1Denote the scale coefficient of the noise added in the k-th decomposition;

[0252] Step 2.4: Repeat the above two steps until the residual component is a monotonic function r k (t), and calculate k IMF components of the pressure pulsation signal.

[0253] Step 3: Use the BiLSTM model to predict each IMF component. Take the first 400 data points of the pressure pulsation data as the input data and the 401st data value as the output data, and divide them into a sample. Subsequent samples are processed in the same way. After normalizing the training samples, input them into the BiLSTM deep learning model, and use ten-fold cross-validation and grid search to optimize the hyperparameters of the prediction model, the number of BiLSTM layers, the number of neurons, and the learning rate. When the number of training reaches the set maximum number, save the model as the trained pressure pulsation prediction model. Superimpose the prediction results of each IMF component to obtain the prediction result.

[0254] Step 4: Retain the prediction model. In this embodiment, a total of 2,870 groups of pressure pulsation samples are collected under the condition of reducing the load to 30% for subsequent prediction. Predict the pressure pulsation signal during the process of reducing the load through the model. The prediction result of the pressure pulsation signal under the 50% load condition is as Figure 10 shown. It can be clearly found that the prediction error of the pressure pulsation under this condition is relatively large. Calculate the pressure pulsation prediction error of each sample based on the actual value and the prediction value of the pressure pulsation signal, and use this error as the monitoring index. Calculate through the following formula:

[0255] D error = |P p - P t |

[0256] where D error represents the pressure pulsation prediction error of the sample; P p represents the predicted value of the sample; P t represents the actual value of the sample.

[0257] Step 5: Use the preset SPOT algorithm to adaptively calculate the threshold. When the value of the monitoring index exceeds the threshold, an alarm will be generated once, and it is considered that the unit is operating under an unstable condition at this time.

[0258] Step 5.1: Let D be any number in the error sequence {D1, D2, …, D K}; calculate the anomaly threshold of the first n anomaly scores, where n < K; initialize the threshold T, obtain the peak point, and fit the parameters according to the generalized Pareto distribution (GPD) of

[0259]

[0260] where, Ft (d) is the GPD fitting distribution function; t is the initial threshold; γ is the shape parameter obtained by GPD fitting, σ is the scale parameter obtained by GPD fitting, and d represents the data points of the data set;

[0261] Step 5.2. Subsequently, adaptively calculate the anomaly threshold Z q :

[0262]

[0263] where t represents the initial threshold, i.e., the 98% quantile of the initialized data, is the estimated value of the scale parameter; γ is the estimated value of the shape parameter. q represents the probability that the sample exceeds the threshold; n is the number of samples, and N t is the number of points with peak value greater than t. For the error D n+i , i ≤ K - n: If d n+i > Z q , it is determined as an anomaly; if T < D n+i < Z q , re - fit the parameters according to GPD to obtain a new Z q ; if D n+i < t, no processing is performed.

[0264] In this embodiment, the primary warning result at 30% load condition is reduced as Figure 11 shown. It can be found that the prediction error is almost 0 during the initial reduction process. As the load further decreases, the error increases rapidly after the generation of the vortex band, and a primary warning is issued, indicating that the unit is operating under an abnormal condition at this time.

[0265] Step 6. Calculate the chaotic characteristic parameters of the primary warning pressure pulsation signal, the attractor diffusion radius C R,a and the correlation dimension C D,a . Combine the chaotic characteristic parameters C R,0 and C D,0 of the pressure pulsation signal under normal operating conditions to define the vortex band warning coefficient: When the warning coefficients C R and C D are greater than the set threshold, it can be judged that there is an unstable vortex band in the draft tube.

[0266] Step 6.1. Adopt the mutual information coefficient method to calculate the mutual information between the pressure pulsation signal and the corresponding delayed signal in the pressure pulsation signal matrix, and adaptively select the optimal delay time τ according to the first local minimum of the mutual information;

[0267]

[0268] Wherein, I is the mutual information, N is the length of the pressure pulsation signal, T is the set delay time between [1, 50]; x i is the pressure pulsation signal at the i-th point; p(x i ) is the probability density of the pressure pulsation signal; p(x i +T) is the probability density of the delayed signal; p(x i , x i +T) is the joint probability density of the two signals; Calculate the mutual information starting from the minimum value of the delay time T, increase the value of T until the first minimum value of I(T) appears, and at this time T is the optimal delay time τ of a certain measurement point;

[0269] Step 6.2: After determining the optimal delay time, adaptively calculate the optimal embedding dimension according to the false nearest neighbor method criterion;

[0270]

[0271]

[0272] Wherein, and respectively represent the distances between two points in the (d + 1)-th and d-th dimensional phase spaces, Rtol is the threshold value between [10, 50], is the j-th phase point vector in the reconstructed phase space, the superscript r represents the phase point coordinates in the reconstructed phase space, is the nearest neighboring phase point of the phase point ; L is the number of phase points after the reconstruction of the pressure pulsation signal; Calculate the proportion of false nearest neighbor points starting from the minimum value of the dimension d, gradually increase the value of d until all false nearest neighbor points disappear, and at this time d is the optimal embedding dimension m of a certain measurement point;

[0273] Table 7 shows the optimal phase space reconstruction parameters of the water turbine pressure pulsation signal under different working conditions in Example 4

[0274] Experimental conditions Embedding dimension m Delay time τ 30% load condition 8 88

[0275] The optimal phase space reconstruction parameters adaptively determined for the pressure pulsation sequence under the 30% load condition are shown in Table 7.

[0276] Step 6.3: Map the pressure pulsation signal to the high-dimensional phase space through the optimal delay time τ and the optimal embedding dimension m, and calculate the high-order matrix Y i :

[0277] Y i =(Y i , Y i+τ …, Y i+(m-1)τ ), i = 1, 2, …, N-(m - 1)τ(4);

[0278] where N is the length of the signal pressure pulsation signal, and Y i is the i-th phase point in the matrix. Assume that the high-order matrix Y i points in the phase space distribution all correspond to a coordinate (x1, y1, z1), (x2, y2, z2)…(x n , y n , z n ). Then the high-order matrix Y of the pressure pulsation signal i The attractor diffusion radius is defined as the maximum value of the distances among all points:

[0279]

[0280] where is the center of all points. In this embodiment, the calculated attractor diffusion radius under the warning condition is 0.0058.

[0281] Step 6.4, Correlation dimension C D,a is defined as the slope of the correlation integral C(R) curve:

[0282]

[0283] R = exp(linspace(log(r min ), log(r max ), M))

[0284] where R is the similarity radius; M i (R) is the number of phase points within the similarity radius R; x i and x k are the i-th and k-th phase points; r min is the minimum value of the similarity radius, r max is the maximum value of the similarity radius, and M is the number of phase points in the linear region. In this embodiment, the calculated attractor correlation dimension under the warning condition is 1.564.

[0285] Combining the chaotic characteristic parameters C R,0 and C D,0 of the pressure pulsation signal in the normal operation state, define the vortex band warning coefficient: When the warning coefficients C R and C D are greater than the set threshold, it is considered that there is a significant spiral vortex band in the draft tube.

[0286]

[0287] where C R,a and C D,a are the diffusion radius and correlation dimension of the pressure pulsation signal under the warning condition; C R,0 and C D,0The diffusion radius and correlation dimension of the pressure pulsation signal under normal operating conditions. In this embodiment, the calculated early warning coefficient C R and C D are 2.04e4 and 0.41 respectively, both exceeding the set early warning threshold, indicating that it is in the vortex band condition at this time. The above embodiment discriminates the abnormal flow state through the adaptive discrimination of the pressure pulsation prediction deviation and does not rely on a complete fault data set. It illustrates the effectiveness of the present invention for the characterization of the vortex band of a Francis turbine.

[0288] The above embodiments are only the preferred technical solutions of the present invention and should not be regarded as limitations on the present invention. The embodiments and features in the present application can be arbitrarily combined with each other without conflict. The protection scope of the present invention should be the technical solutions recorded in the claims, including the equivalent replacement solutions of the technical features in the technical solutions recorded in the claims. That is, the equivalent replacement improvements within this scope are also within the protection scope of the present invention.

Claims

1. A method for characterizing vortex bands in a Francis turbine draft tube based on pressure pulsation signals, characterized in that: The method comprises the following steps: Step 1, collecting a pressure pulsation signal P(t) of the tailwater pipe during normal operation by means of a high-precision pressure pulsation sensor installed at the cross section of the tailwater cone pipe of the Francis turbine, and preprocessing the pressure pulsation signal; Step 2, using a modal decomposition method to decompose the pressure pulsation signal collected in step 1 into multiple IMF components; Step 3: Use the BiLSTM model to predict each IMF component, superimpose the prediction results of each IMF component, obtain the prediction result, and retain the prediction model; Step 4: predict the pressure pulsation signal of the unknown working condition through the model, calculate the pressure pulsation prediction error of each sample based on the actual value and the predicted value of the pressure pulsation signal, and use the error as a monitoring indicator; Step 5: The preset SPOT algorithm is used to adaptively calculate the threshold. When the value of the monitoring indicator exceeds the threshold, an alarm will be generated, indicating that the unit is operating in an unstable condition. Step 6: Calculate the chaotic characteristic parameters of the warning pressure pulsation signal, the attractor diffusion radius C R,a And the correlation dimension C D,a , combined with the chaotic characteristic parameter C of the pressure pulsation signal under normal operating conditions R,0 and C D,0 Definition of vortex belt warning coefficient: When the warning coefficient C R and C D When it is greater than the set threshold, it can be determined that an unstable vortex exists in the tailwater pipe.

2. The method for characterizing the vortex band of the Francis turbine draft tube based on the pressure pulsation signal according to claim 1 is characterized in that: The data preprocessing process in step 1 includes removing abnormal pressure pulsation data according to the Laida criterion and normalizing the data. The normalization process is as follows: Where X i is the pressure pulsation signal data value, X min is the minimum value of the data, X max is the maximum value of the data.

3. The method for characterizing the vortex band of the Francis turbine draft tube based on the pressure pulsation signal according to claim 2 is characterized in that The specific process of step 3 is: the first 400 data points of the pressure pulsation data are used as input data, and the 401st data value is used as output data, and the two are divided into one sample. Subsequent samples are analogous to this. After the training samples are normalized, they are input into the BiLSTM deep learning model. Ten-fold cross validation and grid search are used to optimize the hyperparameters, number of BiLSTM layers, number of neurons, and learning rate of the prediction model. When the number of training times reaches the set maximum number, the modified model is saved as the trained pressure pulsation prediction model.

4. The method for characterizing the vortex band of the Francis turbine draft tube based on the pressure pulsation signal according to claim 3 is characterized in that The specific calculation process of the sample pressure pulsation prediction error in step 4 is: D error =|P p -P t | Where D error represents the sample pressure pulsation prediction error; P p represents the sample prediction value; P t Represents the actual value of the sample.

5. The method for characterizing the vortex band of the Francis turbine draft tube based on the pressure pulsation signal according to claim 4 is characterized in that The specific calculation process of step 5 is: Let D be any number in the error sequence {D1, D2, …, D K}; calculate the anomaly threshold of the first n anomaly scores, where n < K; initialize the threshold T, obtain the peak points, and fit the parameters according to the Generalized Pareto Distribution (GPD): Among them, F t (d) is the GPD fitting distribution function; t is the initial threshold; γ is the shape parameter obtained by GPD fitting, σ is the size parameter obtained by GPD fitting, and d represents the data point of the data set; Then adaptively calculate the abnormal threshold Z q : where t represents the initial threshold, i.e., the 98% quantile of the initialization data, is the estimated value of the size parameter; γ is the estimated value of the shape parameter, q represents the probability that the sample exceeds the threshold; n is the number of samples, and N t is the number of points with peak greater than t; for the error D n+i , i ≤ K - n: if d n+i > Z q , it is determined as abnormal; if T < D n+i < Z q , re - fit the parameters according to GPD to obtain a new Z q ; if D n+i < t, no processing is performed.

6. The method for characterizing the vortex band of the Francis turbine draft tube based on the pressure pulsation signal according to claim 5 is characterized in that The attractor diffusion radius C of the primary warning pressure pulsation signal in step 6 R,a And the correlation dimension C D,a The specific calculation process is: The optimal delay time τ and the optimal embedding dimension m are adaptively calculated by the mutual information coefficient method and the false nearest neighbor method. The pressure pulsation signal is mapped to the high-dimensional phase space by the optimal delay time τ and the optimal embedding dimension m, and the high-order matrix Y is calculated. i : AND i (And i ,AND i+τ …,AND i+(m-1)τ ),i=1,2,…,N-(m-1)τ Where N is the length of the signal pressure pulsation signal, Y i is the i-th phase point in the matrix; assuming that the high-order matrix Y i Each point in the phase space distribution corresponds to a coordinate (x1, y1, z1), (x2, y2, z2)…(x n ,y n ,z n ). Then the high-order matrix Y of the pressure pulsation signal is i The attractor diffusion radius is defined as the maximum distance among all points: in is the center of all points; Correlation Dimension C D,a Defined as the slope of the correlation integral C(R) curve: R=exp(linspace(log(r min ),log(r max ),M)) Where R is the similarity radius; M i (R) is the number of phase points within the similarity radius R; x i and x k are the i-th and k-th phase points; r min is the minimum similarity radius, r max is the maximum value of similarity radius, and M is the number of phase points in the linear region.

7. The method for characterizing the vortex band of the Francis turbine draft tube based on the pressure pulsation signal according to claim 6 is characterized in that In step 6, the vortex belt warning coefficient C R and C D The specific calculation process is: Combined with the chaotic characteristic parameter C of the pressure pulsation signal under normal operation R,0 and C D,0 Definition of vortex belt warning coefficient: When the warning coefficient C R and C D When the value is greater than the set threshold, it is considered that there is a significant spiral vortex in the tailwater pipe. Among them C R,a and C D,a is the diffusion radius and correlation dimension of the pressure pulsation signal under the early warning condition; C R,0 and C D,0 is the diffusion radius and correlation dimension of the pressure pulsation signal under normal operating conditions.

8. The method for characterizing vortex bands in the draft tube of a Francis turbine based on pressure pulsation signals according to claim 7, characterized in that: Warning threshold C R and C D Take 80 and 0.3 respectively.

Images