Wind power blade structure health monitoring method based on multi-field decoupling optical fiber sensing

Abstract

The invention discloses a wind power blade structure health monitoring method based on multi-field decoupling optical fiber sensing, particularly relates to the technical field of wind power equipment state monitoring, and is used for solving the problem of health state misdiagnosis caused by coherent interference of environmental noise and damage signals. The method comprises the following steps: collecting vibration, strain and temperature signals and generating a background noise vector; calculating a frequency domain coherence function of the vibration signal and the noise vector to identify an interference danger frequency band; analyzing the polarization angle distribution of the vibration signal in the frequency band, and inverting the physical path topology of the noise transmitted into the blade; when the path passes through a predefined weak area, improving a damage weight coefficient of a strain signal of a corresponding frequency band and carrying out Hilbert phase shift reconstruction to generate an anti-interference strain feature vector; calculating the energy dissipation rate of the vector in a dangerous frequency band, and judging an energy leakage event caused by damage by combining a viscoelastic constitutive threshold value of the composite material; the damage type and the spatial position are judged according to the time domain correlation between the energy leakage event and the temperature signal, and the hidden damage detection reliability is improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of wind power equipment state monitoring, and particularly relates to a wind power blade structure health monitoring method based on multi-field decoupling optical fiber sensing. BACKGROUND

[0002] In the wind power blade structure health monitoring, the optical fiber sensing technology gradually replaces the traditional electrical sensor due to the advantages of anti-electromagnetic interference, corrosion resistance, etc. The prior art synchronously collects strain, vibration, temperature and other multi-physical field signals by deploying an optical fiber sensor network on the surface of the blade, separates the cross-sensitivity interference by combining a decoupling algorithm, and then realizes crack, delamination and other damage identification. Especially in the complex sound and vibration environment of the wind farm, multi-parameter fusion diagnosis is considered as an effective means to improve the fault detection rate.

[0003] However, when the environmental noise and damage signal meet certain coherence conditions, the background interference such as gearbox vibration transmission and raindrop impact noise inherent in the blade operation may interfere constructively or destructively with the real damage signal in the frequency domain, resulting in unpredictable energy distortion of the decoupled multi-physical quantity feature vector. This coherence camouflage effect makes the system unable to distinguish between noise amplification artifacts and real damage features, causing misdiagnosis of the health state - including false alarm of external vibration as blade crack, or covering up early bearing damage and other hidden faults. SUMMARY

[0004] The present application provides a wind power blade structure health monitoring method based on multi-field decoupling optical fiber sensing to solve the technical problems in the prior art.

[0005] The technical solution of the present application to solve the above technical problems is as follows: The present application provides the following technical solution: The wind power blade structure health monitoring method based on multi-field decoupling optical fiber sensing comprises: S1, synchronously collecting vibration signals, strain signals and temperature signals through an optical fiber sensor network; S2, extracting gearbox vibration transmission path features and raindrop impact noise spectrum features to generate a background noise feature vector; S3, calculating the coherence function value of the vibration signal and the background noise feature vector in the frequency domain, and identifying the frequency band with a coherence function value greater than a preset threshold as an interference danger band; S4, analyzing the polarization angle distribution of the vibration signal in the interference danger band, and inversely calculating the physical path topology of the background noise entering the blade according to the polarization principal axis deviation; S5, when the physical path topology passes through a pre-defined structure weak area, the damage weight coefficient of the strain signal in the corresponding interference danger band is increased and a phase shift reconstruction is performed to generate an anti-interference strain feature vector; S6, calculate the energy dissipation rate of the anti-interference strain feature vector in the interference danger frequency band, and mark it as a damage-induced energy leakage event when the energy dissipation rate exceeds the theoretical threshold of the material; S7, determine the damage type and location of the blade structure according to the time domain correlation between the damage-induced energy leakage event and the temperature signal.

[0006] Further, the vibration signal, the strain signal and the temperature signal are synchronously collected by the optical fiber sensor network, including: The fiber grating sensor array is arranged on the surface of the wind turbine blade, and the vibration signal, the strain signal and the temperature signal are collected in a synchronous sampling manner; The vibration signals are spatially grouped according to the blade spanwise position, and the amplitude variation characteristics of each group of vibration signals in the main vibration frequency band of the blade are extracted; The collected strain signal is subjected to temperature drift compensation operation, and the compensated strain signal is stored according to the blade cross-sectional area; The synchronously collected temperature signal is mapped to a three-dimensional coordinate grid of the blade to form a temperature gradient distribution field.

[0007] Further, the gear box vibration transmission path characteristics and the raindrop impact noise spectrum characteristics are extracted to generate a background noise feature vector, including: The path coupling component transmitted to the blade is separated from the gear box base vibration signal, and the frequency response function amplitude characteristics of the corresponding path coupling component are extracted; The short-time spectrum envelope of the raindrop impact noise signal is analyzed, and the envelope peak frequency and its harmonic distribution characteristics are extracted; The frequency response function amplitude characteristics and the envelope peak frequency and its harmonic distribution characteristics are aligned and fused according to the frequency band; The fused features are normalized to generate a background noise feature vector.

[0008] Further, the coherence function value between the vibration signal and the background noise feature vector in the frequency domain is calculated, and the frequency band with a coherence function value greater than a preset threshold is identified as an interference danger frequency band, including: The vibration signal is subjected to fast Fourier transform to obtain a vibration signal spectrum; The background noise feature vector is divided into a reference spectrum according to the frequency band; The coherence function value between the vibration signal spectrum and the reference spectrum is calculated to generate a coherence function spectrum; The frequency interval with a coherence function value greater than a preset threshold is located in the coherence function spectrum; Adjacent frequency intervals with a coherence function value continuously greater than a preset threshold are merged to form an interference danger frequency band.

[0009] Further, the coherence function value between the vibration signal spectrum and the reference spectrum is calculated to generate a coherence function spectrum, including: convert the background noise feature vector into a reference spectrum consistent with the frequency spectrum resolution of the vibration signal according to the frequency band division; calculate the coherence function value based on the ratio of cross-power spectral density to self-power spectral density, and generate a coherence function spectrum.

[0010] Further, the polarization angle distribution of the vibration signal in the interference risk frequency band is analyzed, and the physical path topology of the background noise entering the blade is inverted according to the polarization principal axis deviation, including: Perform complex analytic signal synthesis on the vibration signal in the interference risk frequency band to separate the orthogonal vibration components; Calculate the instantaneous polarization angle according to the orthogonal vibration components to form a polarization angle time sequence; Statistically analyze the probability density distribution of the polarization angle time sequence to determine the polarization principal axis direction angle; Calculate the standard deviation of the polarization angle time sequence relative to the polarization principal axis direction angle as the polarization principal axis deviation; When the polarization principal axis deviation exceeds the health reference value, the azimuth and pitch angle of the background noise entering path are inverted according to the mapping relationship between the polarization principal axis direction angle and the blade geometric coordinate system; Combine the azimuth and pitch angle to generate the physical path topology.

[0011] Further, the vibration signal in the interference risk frequency band is subjected to complex analytic signal synthesis to separate the orthogonal vibration components, including: Perform Hilbert transform on the vibration signal to generate orthogonal phase components; Combine the original signal and the orthogonal phase component to form a complex analytic signal, and separate the real and imaginary parts as orthogonal vibration components.

[0012] Further, when the physical path topology passes through a pre-defined structural weak area, the damage weight coefficient of the strain signal in the corresponding interference risk frequency band is increased and phase shift reconstruction is performed to generate an anti-interference strain feature vector, including: Match the physical path topology with the spatial coordinate range of the pre-defined structural weak area for judgment; When the physical path topology falls within the range of the pre-defined structural weak area, perform time domain weighting operation on the strain signal in the interference risk frequency band to increase the damage weight coefficient; Perform Hilbert transform phase shift processing on the strain signal with the increased damage weight coefficient to construct orthogonal strain signal components; Synthesize the orthogonal strain signal components and the original strain signal to generate an anti-interference strain feature vector.

[0013] Further, the energy dissipation rate of the anti-interference strain feature vector in the interference risk frequency band is calculated, and when the energy dissipation rate exceeds the material theoretical threshold, it is marked as a damage-induced energy leakage event, including: Extracting time domain signal segments corresponding to the interference danger frequency band from the anti-interference strain feature vector; Calculating the time domain attenuation rate of the corresponding time domain signal segment as the energy dissipation rate; Determining the material theoretical threshold according to the viscoelastic constitutive relation of the blade composite material; Comparing the energy dissipation rate with the material theoretical threshold; When the energy dissipation rate exceeds the material theoretical threshold, mark the current time domain signal segment as a damage-induced energy leakage event.

[0014] Further, according to the time domain correlation between the damage-induced energy leakage event and the temperature signal, the type and position of the blade structure damage are distinguished, comprising: Extracting the temperature signal change amount in the time window corresponding to the damage-induced energy leakage event; Calculating the time domain lag amount of the temperature signal change amount and the starting time of the damage-induced energy leakage event; Determining the heat conduction characteristic category based on the time domain lag amount, and associating the type of blade structure damage; Locating the temperature abnormal area according to the temperature gradient distribution field of the damage-induced energy leakage event occurrence period; Mapping the center coordinates of the temperature abnormal area to the blade three-dimensional coordinate grid, and outputting the position of the blade structure damage.

[0015] The beneficial effects of the present application are: 1. By identifying the coherent danger frequency band of the vibration signal and the background noise, and combining the polarization angle distribution to analyze the noise physical path topology, the coherent camouflage effect is actively dissociated. Unlike traditional static decoupling, the spatial relationship between the noise transmission path and the weak area of the blade is adaptively improved to increase the strain signal weight, and the interference noise attached to the real damage signal is stripped through the phase shift reconstruction technology, which eliminates the energy distortion caused by environmental interference such as gearbox vibration transmission and raindrop impact, so that the decoupled multi-physical field signal has the real damage characterization ability.

[0016] 2. An energy leakage and material characteristic correlation mechanism is established, an energy dissipation threshold is set based on the viscoelastic constitutive relation of the composite material, a damage event is determined by the energy abnormal leakage of the anti-interference strain feature vector, and damage classification and positioning are realized by combining the time domain response characteristics of the temperature signal, so that the system can accurately capture the early energy dissipation characteristics of hidden damage such as blade matrix cracking and delamination, avoid the false alarm of external vibration as a crack caused by the coherent effect of traditional methods, and solve the problem that hidden faults such as bearing damage are covered by noise, significantly improving the accuracy and reliability of damage identification in complex acoustic vibration environment. BRIEF DESCRIPTION OF DRAWINGS

[0017] Figure 1A flowchart of a wind turbine blade structural health monitoring method based on multi-field decoupling optical fiber sensing of the present application is shown in the figure. Figure 2 A flowchart of generating an anti-interference strain feature vector of the present application is shown in the figure. DETAILED DESCRIPTION

[0018] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative work fall within the scope of protection of the present application.

[0019] Embodiment: Figure 1 The wind turbine blade structural health monitoring method based on multi-field decoupling optical fiber sensing of the present application is shown in the figure, which comprises: S1, synchronously collecting vibration signals, strain signals and temperature signals through an optical fiber sensor network; S2, extracting gear box vibration transmission path features and raindrop impact noise spectrum features to generate background noise feature vectors; S3, calculating the coherence function values of the vibration signals and the background noise feature vectors in the frequency domain, and identifying the frequency bands with coherence function values greater than a preset threshold as interference danger frequency bands; S4, analyzing the polarization angle distribution of the vibration signals in the interference danger frequency bands, and inversely calculating the physical path topology of the background noise entering the blade according to the polarization principal axis deviation degree; S5, when the physical path topology passes through a pre-defined structural weak area, the damage weight coefficient of the strain signals in the corresponding interference danger frequency band is increased and a phase shift reconstruction is performed to generate an anti-interference strain feature vector; S6, calculating the energy dissipation rate of the anti-interference strain feature vector in the interference danger frequency band, and marking as a damage-induced energy leakage event when the energy dissipation rate exceeds a material theoretical threshold; S7, judging the blade structure damage type and position according to the time domain correlation between the damage-induced energy leakage event and the temperature signals.

[0020] In the process of deploying FBG sensor array on the surface of wind turbine blade, the sensor nodes are fixed on the leading edge, trailing edge and spar locations of the blade with equal spacing along the spanwise direction. Each sensor node contains a vibration sensing unit, a strain sensing unit and a temperature sensing unit, and all sensing units transmit signals through the same optical fiber link. A synchronous trigger acquisition mode is adopted, in which the central processing unit sends a synchronization pulse signal to all sensor nodes to ensure that the timestamps of vibration signals, strain signals and temperature signals are aligned to the level of microseconds. During the synchronous acquisition process, the vibration signal sampling frequency is set to 2000 Hz, which is based on the consideration that it covers the main vibration frequency band of 0-1000 Hz of the blade and meets the Nyquist sampling theorem; the strain signal sampling frequency is set to 500 Hz, which is based on the consideration that the static strain change rate is less than 1 Hz; the temperature signal sampling frequency is set to 1 Hz, which is based on the consideration that the thermal inertia characteristic of the blade temperature change time constant is greater than 60 seconds.

[0021] After the acquisition is completed, the vibration signals are spatially grouped according to the spanwise position of the blade: the blade is divided into ten equal-length spanwise sections from the root to the tip, and each section contains three adjacent sensor nodes. For each vibration signal group of the spanwise section, the amplitude variation characteristics in the main vibration frequency band of the blade are extracted: the band-pass filter is used to retain the blade first-order bending vibration frequency band signal of 1-100 Hz, and then the root mean square value of the peak-to-peak value of the signal in this frequency band is calculated as the vibration amplitude characteristic of the spanwise section. For example, the vibration amplitude characteristic acquisition process of the root section includes: intercepting the vibration signals of the three nodes of the root section for 300 seconds, performing 1-100 Hz band-pass filtering on each node signal, calculating the peak-to-peak value of the filtered signal in each 10-second time window, and finally taking the root mean square value of all time window peak-to-peak values as the section characteristic.

[0022] Temperature drift compensation is performed on the acquired strain signals: the compensation is realized through a pre-established temperature-strain coupling coefficient lookup table, which is obtained through the following calibration process: in the no-load state, the blade sample is placed in a temperature-controlled box, and the temperature is raised from -20 degrees Celsius to 60 degrees Celsius in steps of 5 degrees Celsius, and the original readings of the strain sensor at each temperature point are recorded. The rate of change of the readings with temperature is defined as the coupling coefficient. During actual compensation, the coupling coefficient at the corresponding position is queried according to the current temperature signal value, and the compensation formula "compensated strain value = original strain value minus the difference between current temperature and reference temperature multiplied by coupling coefficient" is used for calculation, where the reference temperature is set to 25 degrees Celsius. The compensated strain signals are stored according to the blade cross-sectional area: the blade cross section is divided into four regions: leading edge region, trailing edge region, pressure surface region and suction surface region, and the strain signals of the same region are stored in independent data channels, and each channel data is organized in the form of a two-dimensional matrix corresponding to time sequence and spatial position.

[0023] Mapping the synchronously collected temperature signals to a three-dimensional coordinate grid of the blade: a three-dimensional coordinate grid containing one hundred thousand grid nodes was established based on the blade computer-aided design model, with a grid node spacing of 0.01 meters. When mapping the temperature signals, first determine the nearest neighbor grid node of each temperature sensor in the three-dimensional grid, and assign the measured temperature value of the sensor to the node; for grid nodes without sensors, the inverse distance weighted interpolation method is used to calculate the temperature value, specifically, the distances of the node to the adjacent three sensors are calculated, and the weighted average value of each sensor temperature value divided by the square of the distance is taken. After completing the assignment of all nodes, calculate the temperature difference between adjacent grid nodes and divide by the node spacing to generate the temperature gradient distribution field. The distribution field is stored in the form of a three-dimensional array, and each element of the array contains three components, representing the temperature change rate of the grid node along the spanwise direction, chordwise direction and thickness direction of the blade, with a unit of degrees Celsius per meter. For example, the temperature gradient calculation process of a grid node in the tip region is as follows: obtain the temperature values of the node and its six adjacent nodes along the three coordinate axes, calculate the temperature difference between the adjacent nodes in the spanwise direction divided by 0.01 meters, the temperature difference between the adjacent nodes in the chordwise direction divided by 0.01 meters, and the temperature difference between the adjacent nodes in the thickness direction divided by 0.01 meters, to obtain the component values in the three directions.

[0024] When separating the path coupling component transmitted to the blade from the gearbox base vibration signal, the vibration signal collected by the three-axis acceleration sensor mounted on the gearbox base is obtained, and the signal sampling frequency is set to 5000 Hz to meet the high-frequency vibration analysis requirement. The adaptive noise cancellation technology is used to separate the path coupling component: the vibration signal at the root of the blade is taken as the main input signal, and the vibration signal of the gearbox base is taken as the reference input signal, and the least mean square algorithm is used to iteratively adjust the finite impulse response filter coefficients, and the filter order is set to 128, which is determined by pre-experiment according to the signal correlation time. When the filter coefficient variation is less than the convergence threshold of 0.001, the iteration is stopped, and the enhanced path coupling component is output. Extract the frequency response function amplitude feature of the path coupling component: calculate the frequency response function of the separated path coupling component, and record the amplitude value corresponding to each frequency point in the range of 1 to 2000 Hz with a resolution interval of 1 Hz, to form a frequency-amplitude feature sequence. For example, the amplitude feature acquisition process at 100 Hz is as follows: a 100 Hz sinusoidal excitation is applied to the gearbox base, the response amplitude at the root of the blade is measured, and the ratio of the response amplitude to the excitation amplitude is calculated as the feature value at this frequency point, which is dimensionless.

[0025] When analyzing the short-time spectral envelope of raindrop impact noise signals, the raindrop impact noise signals are collected using the sound pressure sensors installed on the blade surface, which are arranged and synchronously collected in S1. The signal sampling frequency is set to 10000 Hz to cover the broadband characteristics of raindrop impact. The noise signal is processed using short-time Fourier transform: the signal is divided into time windows of 50 milliseconds, and adjacent time windows overlap by 25 milliseconds, which is based on the average duration of raindrop impact events. The signal in each time window is subjected to a 1024-point fast Fourier transform to obtain a frequency spectrum. The envelope peak frequency and its harmonic distribution characteristics are extracted: in the frequency spectrum of each time window, the first three peak points with the largest amplitude are identified as the envelope peak frequency, and it is detected whether there are harmonic components with an amplitude exceeding 30% of the fundamental wave amplitude at integer multiples of these peak frequencies. The threshold is determined by statistical analysis of laboratory raindrop impact experiments. Record all detected fundamental and harmonic frequency values to form a feature set. For example, in a certain time window, the fundamental frequency of 120 Hz is detected, and the amplitude of its second harmonic of 240 Hz is 35% of the fundamental wave, so the feature pair (120 Hz, 240 Hz) is recorded.

[0026] The frequency response function amplitude characteristics and the envelope peak frequency and its harmonic distribution characteristics are aligned and fused by frequency band: the frequency range of 1 to 2000 Hz is divided into 20 equal-width frequency bands, each with a width of 100 Hz. In each frequency band, the fusion operation is performed: first, the arithmetic mean of the frequency response function amplitude corresponding to all frequency points in the frequency band is calculated as the first characteristic value; second, the number of envelope peak frequencies and their harmonic distribution characteristics appearing in the frequency band is counted as the second characteristic value; finally, the two characteristic values are combined according to the formula "the fusion characteristic value is equal to the first characteristic value multiplied by the second characteristic value". The frequency band boundary processing rule is: if the envelope peak frequency is equal to the frequency band boundary value, it is counted into the adjacent two frequency bands at the same time. For example, in the frequency band of 100 Hz to 200 Hz, the frequency response function amplitude average is 0.25, and the envelope peak frequency appears 3 times, so the fusion characteristic value is 0.25x3=0.75.

[0027] Normalization of the fused features: first, traverse all the 20 frequency band fused feature values to find the maximum and minimum values; then map each frequency band fused feature value to the range of 0 to 1 according to the formula "normalized feature value equals current feature value minus minimum value divided by the difference between maximum value and minimum value". In special cases, if the maximum value equals the minimum value, set all feature values to 0.5; if the feature value exceeds the original range, truncate it to the nearest boundary value. The normalized 20 frequency band feature values are arranged in order from low to high frequency to form a 20-dimensional background noise feature vector. The feature vector is stored in a floating-point number array format, and each element corresponds to the normalized feature value of a frequency band. For example, when the minimum value of all frequency band fused feature values is 0.1 and the maximum value is 0.9, and the original feature value of a certain frequency band is 0.5, its normalized value is (0.5-0.1) ÷ (0.9-0.1) = 0.5.

[0028] When performing fast Fourier transform on the vibration signal, the time domain signal of the target spanwise section is extracted from the vibration signal packet data stored in S1. The signal is windowed using a Hanning window function, and the window length is set to 1024 sampling points, for example, which is determined to be 0.5 Hz according to the frequency resolution requirement. Fast Fourier transform is performed on the windowed signal to obtain the complex form of the vibration signal spectrum, and the spectrum range covers 0 to 1000 Hz, with each frequency point interval of 0.5 Hz. The spectrum data is stored in a two-dimensional array with real and imaginary parts, and the first dimension corresponds to the frequency order, and the second dimension stores the real and imaginary values. For example, when processing a certain 10-second vibration signal in the midspan section, 20 segments of 1024-point signals are intercepted (sampling rate 2000 Hz), and the modulus value of each segment after transformation is averaged as the final spectrum.

[0029] When converting the background noise feature vector into a reference spectrum according to the frequency band division, first confirm that the background noise feature vector is a 20-dimensional normalized array generated by S2, and each dimension corresponds to a 100 Hz frequency band (1-100 Hz to 1901-2000 Hz). To match the 0.5 Hz resolution of the vibration signal spectrum, perform linear interpolation on each 100 Hz frequency band: assign the feature value of the frequency band to the center frequency point of the frequency band, and connect the adjacent center points with a straight line to generate a continuous spectrum. After interpolation, the frequency axis of the reference spectrum is completely aligned with the vibration signal spectrum, and each 0.5 Hz frequency point has a corresponding amplitude. For example, the 5th frequency band (401-500 Hz) has a feature value of 0.8, so the center frequency 450 Hz has an amplitude of 0.8, and each 0.5 Hz decreases to 0 at the frequency band boundaries 401 Hz and 500 Hz.

[0030] When calculating the coherence function value between the vibration signal spectrum and the reference spectrum, the following procedures are processed: first, the cross-power spectral density is calculated, and the conjugate product of the vibration signal spectrum and the reference spectrum is taken, and the average of the windowing results is taken; second, the vibration signal self-power spectral density (the average of the square of the vibration signal spectrum modulus) and the reference signal self-power spectral density (the average of the square of the reference spectrum modulus) are calculated respectively; finally, the coherence function value is calculated point by point according to the formula "coherence function value = |cross-power spectral density| vibration signal self-power spectral density x reference signal self-power spectral density)". The calculation result forms a coherence function spectrum and is stored as a one-dimensional array. For example, at 300 Hz, if the square of the modulus of the cross-power spectral density is 0.25, the vibration self-power spectral density is 1.0, and the reference self-power spectral density is 0.5, then the coherence value is 0.25 / (1.0x0.5)=0.5. 2 When locating the frequency interval with a coherence function value greater than a preset threshold in the coherence function spectrum, the preset threshold is set to 0.6, and the threshold is determined by historical data statistics of healthy blades: 100 groups of blade vibration signals and background noise signals in a non-damaged state are collected, the coherence value distribution at each frequency point is calculated, and the 95% quantile value is taken as the threshold. Scan the entire coherence function spectrum, mark all frequency points with a coherence value greater than 0.6, and aggregate the continuously marked frequency points into independent intervals. For example, all three frequency points (150.0 Hz, 150.5 Hz, 151.0 Hz) in the range of 150.0 Hz to 152.5 Hz have a coherence value greater than 0.6, forming an interval [150.0 Hz, 151.0 Hz].

[0031] When merging adjacent frequency intervals with a coherence function value continuously greater than a preset threshold, the maximum merging interval is set to 5 Hz, which is determined by experimental statistics according to the noise correlation bandwidth. If the frequency boundary difference of two intervals is less than 5 Hz, the two intervals are merged into a single continuous interval. The final interval formed after merging is the interference danger frequency band. For example, there are an interval [150.0 Hz, 151.0 Hz] and an interval [151.5 Hz, 153.0 Hz], and the boundary difference 151.5 Hz-151.0 Hz=0.5 Hz is less than 5 Hz, which is merged into the [150.0 Hz, 153.0 Hz] frequency band.

[0032] ​When synthesizing the complex analytic signal from the vibration signal in the interference risk band, the vibration signal time domain data in the corresponding frequency range is extracted from the interference risk band identified by S3. First, the Hilbert transform is performed on the vibration signal: it is realized by a finite impulse response filter, and the filter coefficients are designed according to the 90-degree phase shift requirement, and the order is set to 256, for example, to meet the phase accuracy requirement. The signal output by the Hilbert transform is the quadrature phase component, which maintains the same sampling rate and length as the original vibration signal. The original vibration signal and the quadrature phase component are combined to form a complex analytic signal: the original signal is taken as the real part, and the quadrature phase component is taken as the imaginary part to form a complex sequence. The real part and the imaginary part are separated as orthogonal vibration components and stored as two independent one-dimensional arrays. For example, for the interference risk band signal with a center frequency of 150 Hz and a bandwidth of 5 Hz, after band-pass filtering, the above processing is performed to obtain the real part vibration component array and the imaginary part vibration component array.

[0033] When calculating the instantaneous polarization angle according to the orthogonal vibration component, the instantaneous polarization angle is calculated point by point according to the formula "instantaneous polarization angle = arctan (imaginary part value / real part value)", and the calculation result is in radians. During the calculation process, when the real part is zero, the polarization angle is directly determined as +π / 2 or -π / 2 according to the sign of the imaginary part. The instantaneous polarization angles of all calculation points are arranged in time sequence to form a polarization angle time sequence, which is stored as a one-dimensional floating point array. For example, at a certain time, the real part is 0.8 and the imaginary part is 0.6, then the polarization angle is arctan (0.6 / 0.8) ≈ 0.6435 radians.

[0034] When calculating the probability density distribution of the polarization angle time sequence, the range of -π / 2 to π / 2 radians is divided into 36 equal intervals, and each interval has a width of π / 36 radians (5 degrees). The frequency of each angle value in the polarization angle time sequence falling into each interval is calculated, and the proportion of the frequency in the total number of points is taken as the probability density value. The principal axis direction angle is determined as the center angle value of the interval corresponding to the maximum probability density. If multiple intervals have the same maximum probability density, the arithmetic mean of the center angles of these intervals is taken. For example, it is found that the probability density of interval 10-15 degrees (center angle 12.5 degrees) is the highest at 0.12, and the principal axis direction angle is 12.5 degrees.

[0035] When calculating the standard deviation of the polarization angle time sequence relative to the principal axis direction angle, first, all polarization angle values are subtracted from the principal axis direction angle to obtain a deviation angle sequence, and then the standard deviation is calculated according to the formula "standard deviation = sqrt(sum(deviation angle 2 The number of data points -1). The calculation result is taken as the polarization principal axis deviation, and the unit is radian. For example, when the deviation angle sequence is [0.1, -0.2, 0.05] radians, the standard deviation = sqrt((0.01+0.04+0.0025) / 2) ≈ sqrt(0.02625) ≈ 0.162 radians. When the polarization principal axis deviation exceeds the health benchmark value, the path inversion process is triggered. The health benchmark value is determined by statistical analysis of the health blade historical data: 100 sets of polarization data in undamaged state are collected, the polarization principal axis deviation value of each interference band is calculated, and the 99th maximum value is taken as the benchmark value after ascending arrangement. If the current polarization principal axis deviation is greater than the benchmark value, the inversion is performed according to the mapping relationship between the polarization principal axis direction angle and the blade geometric coordinate system: first, the blade coordinate system is established, taking the blade root center as the origin, the spanwise direction as the positive direction of the Z axis, the chordwise direction as the positive direction of the X axis, and the thickness direction as the positive direction of the Y axis; then the projection angle of the incoming path in the XY plane is calculated according to the formula "azimuth angle=polarization principal axis direction angle×proportionality coefficient K1", and the proportionality coefficient K1 is determined by the sound wave propagation calibration experiment, which measures the polarization principal axis direction angle corresponding to the sound wave at different incident angles; the angle between the path and the XY plane is calculated according to the formula "pitch angle=arctan(polarization principal axis direction angle×proportionality coefficient K2)", and the proportionality coefficient K2 is determined by the material acoustic test experiment according to the acoustic impedance characteristics of the blade material. For example, when the polarization principal axis direction angle is 30 degrees, K1=1.2, and K2=0.8, the azimuth angle=30×1.2=36 degrees, and the pitch angle=arctan(30×0.8)=arctan(24)≈87.4 degrees.

[0036] When the azimuth angle and the pitch angle are combined to generate the physical path topology, the two angle values are combined into a binary tuple of azimuth angle value and pitch angle value, and the corresponding interference dangerous band information is associated. The physical path topology is stored in a structured data format, including a start time field, a duration field, an azimuth angle field, a pitch angle field, and a frequency band range field. For example, a topology record is generated: start time 102.3 seconds, duration 5.6 seconds, azimuth angle 36.0 degrees, pitch angle 87.4 degrees, and frequency band range [148.5, 151.5] Hz.

[0037] Figure 2 The flowchart of generating the anti-interference strain feature vector of the present application is given, and the specific implementation of the generation logic of the anti-interference strain feature vector is: When matching the physical path topology with the spatial coordinate range of the predefined structural weak zone, the data of the predefined structural weak zone in the three-dimensional model of the blade is first loaded. The predefined structural weak zone is determined in advance through finite element stress analysis of the blade. The specific method is to calculate the principal stress distribution of each region of the blade under the rated load, and mark the region with a principal stress value exceeding 80% of the material yield strength as a weak zone. Each weak zone is defined by a minimum bounding box in space, and the bounding box includes six parameters of the minimum X coordinate, the maximum X coordinate, the minimum Y coordinate, the maximum Y coordinate, the minimum Z coordinate, and the maximum Z coordinate. The matching judgment process is as follows: the azimuth angle parameter and the pitch angle parameter in the physical path topology are extracted, and the path ray parameter equation is established in combination with the blade position corresponding to the path starting time point; it is detected whether the path ray intersects with the minimum bounding box of the predefined structural weak zone, and the intersection determination is based on the calculation of the intersection point coordinates of the ray and the six planes of the bounding box. If there is an intersection point and the intersection point coordinates satisfy Xmin≤Xp≤Xmax and Ymin≤Yp≤Ymax and Zmin≤Zp≤Zmax, it is determined that the intersection exists. Among them, Xmin and Xmax represent the minimum boundary value and the maximum boundary value of the bounding box in the X axis direction respectively; Ymin and Ymax represent the minimum boundary value and the maximum boundary value of the bounding box in the Y axis direction respectively; Zmin and Zmax represent the minimum boundary value and the maximum boundary value of the bounding box in the Z axis direction respectively; Xp, Yp, and Zp represent the three-dimensional coordinate values of the intersection point of the ray and the bounding box plane. When at least one intersection point is detected, it is determined that the physical path topology falls within the range of the predefined structural weak zone. For example, the bounding box range of a certain weak zone is X[1.2, 1.5]m, Y[0.3, 0.6]m, and Z[15.2, 15.8]m. The path ray parameter equation is calculated to intersect with the Y=0.3 plane at the coordinate (1.35, 0.3, 15.5). The coordinate satisfies 1.2≤1.35≤1.5, 0.3≤0.3≤0.6, and 15.2≤15.5≤15.8, so it is determined that the matching is successful.

[0038] When the physical path topology falls into the range of the predefined structure weak zone, a time domain weighting operation is performed on the strain signal in the interference dangerous frequency band. First, the time domain signal corresponding to the interference dangerous frequency band is extracted from the compensated strain signal stored in S1, and the sampling frequency is 500 Hz. The initial value of the damage weight coefficient is 1.0, and the promotion rule is: when the path passes through the weak zone, the weight coefficient is multiplied by the promotion factor K, and the value of K is determined by the historical damage data statistics. The specific operation is: multiply the value of each sampling point of the strain signal by the current damage weight coefficient to form a weighted strain signal sequence. The promotion factor K calculation process is: analyze 100 groups of historical damage event data, and count the proportion P of the path passing through the weak zone at the time of damage occurrence. According to the formula K = 1 + 2 × P, the promotion factor is calculated. For example, if P = 0.3 is counted, then K = 1 + 2 × 0.3 = 1.6, and the original weight 1.0 is promoted to 1.6. When the path does not pass through the weak zone, the damage weight coefficient remains 1.0 unchanged.

[0039] When the strain signal after promoting the damage weight coefficient is subjected to Hilbert transform phase shift processing, a Hilbert transform filter is first designed. The filter adopts a finite impulse response structure, and the order is set to 128, for example, to meet the phase accuracy requirement, and the passband range is set to the frequency boundary of the current interference dangerous frequency band. The weighted strain signal is subjected to Hilbert transform to output the orthogonal phase component. The orthogonal strain signal component is constructed as: the original weighted strain signal is used as the real part component, and the Hilbert transform output is used as the imaginary part component. For example, for a strain signal with a center frequency of 150 Hz and a bandwidth of 5 Hz, after being processed by a 128-order Hilbert filter, an orthogonal imaginary part component with the same length as the original signal is obtained.

[0040] When the original strain signal is synthesized with the orthogonal strain signal component, the real part component and the imaginary part component are combined into a complex sequence to form an anti-interference strain feature vector. The feature vector is a complex array, and the array length is equal to the number of strain signal sampling points. Each element contains a real part value and an imaginary part value. The feature vector is stored with additional metadata, including timestamp, spatial position, damage weight coefficient, interference frequency band range, and other information.

[0041] When extracting the time-domain signal segment corresponding to the interference risk frequency band from the anti-interference strain feature vector, first load the anti-interference strain feature vector data generated by S5. The feature vector is in the form of a complex array, containing real and imaginary parts. The extraction process is as follows: according to the interference risk frequency band range determined by S3, for example [148.5, 151.5] Hz, perform band-pass filtering on the feature vector data. The band-pass filter uses a finite impulse response structure, and the order is set to 128 to meet the passband accuracy requirement, the passband boundary is set to the upper and lower limit frequencies of the interference risk frequency band, and the transition bandwidth is set to 2 Hz. After filtering, the time-domain signal segment corresponding to the frequency band is retained, which contains the start time point, end time point, and complex strain data sequence. For example, from the 5.6-second segment from 102.3 seconds to 107.9 seconds, a signal segment containing 1000 complex data points is extracted, with a sampling interval of 0.002 seconds.

[0042] When calculating the time-domain decay rate of the corresponding time-domain signal segment as the energy dissipation rate, the logarithmic decay rate calculation method is used. First, calculate the envelope of the signal segment: take the modulus value sequence of the complex signal as the envelope data, and the modulus value calculation formula is |z| = sqrt(real part 2 imaginary part 2 ). Then calculate the decay rate according to the formula "decay rate = [ln(An / A{n+m})] / (m x Δt)", where An represents the nth envelope peak value, A{n+m} represents the nth+m envelope peak value, Δt is the sampling interval, and m is the interval period number. In actual operation, the average value of the continuous 5 envelope peak values is taken as the energy dissipation rate, with a unit of per second (s. For example, the envelope peak value sequence of a certain segment is [1.2, 0.8, 0.6, 0.4, 0.3] microstrain, and the sampling interval is 0.002 seconds, then the decay rate = [ln(1.2 / 0.8) + ln(0.8 / 0.6) + ln(0.6 / 0.4) + ln(0.4 / 0.3)] / (4 x 0.002) ≈ 125 s.+ -1 ) -1 When determining the material theoretical threshold value according to the viscoelastic constitutive relation of the blade composite material, the generalized Maxwell model is used for theoretical derivation. First, obtain the blade material parameters: obtain the storage modulus E', loss modulus E'', loss factor tan δ, etc. through dynamic mechanical analysis experiment. Then calculate according to the formula "theoretical threshold value = π x fc x tan δ", where fc is the center frequency of the interference risk frequency band, tan δ is the material loss factor, and π is the constant of circular ratio. Material parameters are tested under standard environmental conditions: temperature 20±1℃, relative humidity 50±5%, test frequency range covers 10 to 2000 Hz. For example, when the center frequency is 150 Hz and the material loss factor is 0.02, the theoretical threshold value = π x 150 x 0.02 ≈ 9.42 s. -1 When comparing the size of the energy dissipation rate with the material theoretical threshold value, a numerical comparison operation is directly performed. The calculated energy dissipation rate and the material theoretical threshold value of the corresponding frequency band are compared in size, and the judgment condition is "energy dissipation rate > material theoretical threshold value". The comparison result is output as a Boolean value, true indicating that it exceeds the threshold value, and false indicating that it does not exceed. For example, the energy dissipation rate 125s is greater than the theoretical threshold value 9.42s, so the judgment result is true. When the energy dissipation rate is less than or equal to the material theoretical threshold value, no operation is performed, and the original state is maintained. -1-1 When the energy dissipation rate exceeds the material theoretical threshold value, the current time domain signal segment is marked as a damage-induced energy leakage event. The marking operation is: adding a "damage mark" field in the metadata of the anti-interference strain feature vector, and setting the field value to "energy leakage"; at the same time, an independent event record is generated, including the start time field, the duration field, the spatial position field, the energy dissipation rate value field, the theoretical threshold value field, and the frequency band range field. For example, a record is generated: start time 102.3 seconds, duration 5.6 seconds, spatial position "front edge area", energy dissipation rate 125s, theoretical threshold value 9.42s, frequency band range [148.5, 151.5] Hz, damage mark "energy leakage". -1-1 When extracting the temperature signal change amount in the time window corresponding to the damage-induced energy leakage event, the S6 marked damage-induced energy leakage event record is first loaded. The temperature signal segment in the corresponding time window is extracted from the temperature signal data stored in S1, and the time window is set to the start time to the end time of the damage-induced energy leakage event. The temperature signal change amount is calculated as the difference between the maximum and minimum values of the temperature signal in the current time window, in units of degrees Celsius. For example, the time window of an event is from 102.3 seconds to 107.9 seconds, and the temperature signal rises from 25.6 degrees Celsius to 28.3 degrees Celsius, so the temperature signal change amount is 28.3-25.6=2.7°C.

[0043] When calculating the time domain lag amount of the temperature signal change amount and the start time of the damage-induced energy leakage event, the time point at which the temperature signal change amount reaches the preset change threshold is located as the temperature significant change time. The preset change threshold is set to 0.5 degrees Celsius, which is determined according to the measurement accuracy of the temperature sensor and the statistical characteristics of the background temperature fluctuation: 100 groups of healthy blade temperature data are analyzed, and the temperature fluctuation standard deviation σ is calculated, and the value of σ is taken as the change threshold. The time domain lag amount is calculated as the temperature significant change time minus the start time of the damage-induced energy leakage event, in units of seconds. For example, the event start time is 102.3 seconds, and the temperature reaches 26.1 degrees Celsius (change amount 0.5 degrees Celsius) at 103.5 seconds, so the time domain lag amount is 103.5-102.3=1.2 seconds.

[0044] When determining the heat conduction characteristic category based on the time domain lag quantity, a mapping relationship between the time domain lag quantity and the heat conduction characteristic is established. The mapping rule is: when the time domain lag quantity is less than 0.5 seconds, it is determined as fast heat conduction characteristic; when the time domain lag quantity is between 0.5 seconds and 2 seconds, it is determined as medium-speed heat conduction characteristic; when the time domain lag quantity is greater than 2 seconds, it is determined as slow-speed heat conduction characteristic. The heat conduction characteristic category is related to the blade structure damage type: fast heat conduction corresponds to matrix cracking damage, medium-speed heat conduction corresponds to delamination damage, and slow-speed heat conduction corresponds to fiber fracture damage. The mapping relationship is established by material heat conduction experiment and damage sample verification: three types of damage samples are prepared, the heat wave propagation speed is measured, and the conversion relationship between the speed and the time domain lag quantity is established. For example, the time domain lag quantity 1.2 seconds belongs to the range of 0.5 seconds to 2 seconds, corresponding to the medium-speed heat conduction characteristic, and is related to the delamination damage type.

[0045] When locating the temperature abnormal area according to the temperature gradient distribution field of the damage enabled energy leakage event occurrence period, first, the blade surface temperature sensor array data is obtained. The temperature gradient distribution field is generated by a bicubic spline interpolation algorithm, and the spatial resolution is 5 millimeters. The positioning method is: calculating the temperature change rate (temperature change per unit time) of each grid point, marking the position where the temperature change rate is greater than 2 times the standard deviation of the health benchmark value as a temperature abnormal point; clustering adjacent temperature abnormal points to form a temperature abnormal area, and the area boundary is determined by a convex hull algorithm. The health benchmark value is obtained by statistical analysis of 100 sets of healthy blade data: collecting temperature data under the same working condition, calculating the mean and standard deviation of the temperature change rate at each position. For example, the mean of the health temperature change rate at a specific position is 0.3 degrees Celsius per second, and the standard deviation is 0.05 degrees Celsius per second, then the abnormal threshold is 0.3+2×0.05=0.4℃ / s.

[0046] When mapping the temperature abnormal area center coordinates to the blade three-dimensional coordinate grid, first, the geometric center coordinates of the temperature abnormal area are calculated. The geometric center calculation method is: taking the arithmetic average of all sensor position coordinates in the area. Then, the center coordinates are converted to the blade three-dimensional coordinate system: according to the blade installation azimuth angle parameter and the pitch angle parameter, the sensor local coordinates are converted to the blade global coordinates through the coordinate rotation matrix. The output blade structure damage position is a three-dimensional coordinate point, the format is X coordinate value, Y coordinate value, Z coordinate value three-tuple, and the precision is 0.001 meters. For example, the temperature abnormal area center in the sensor coordinate system is (0.35, 1.42), and after coordinate conversion, the blade global coordinates (15.235, 0.782, 3.418) are obtained.

[0047] The calculations involved in the embodiments are all de-dimensioned to calculate their numerical values, and the preset parameters and threshold values in the calculations are set by those skilled in the art according to the actual situation.

[0048] It should be noted that the application can be deployed in the device itself to realize embedded application, or run on PC terminal or other terminal with user interface, so as to meet various hardware environment and use requirements.

[0049] The above-described embodiments can be implemented in whole or in part by software, hardware, firmware, or any combination thereof. When implemented by software, the above-described embodiments can be implemented in whole or in part in the form of a computer program product. The computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, the processes or functions described in the embodiments of the present application are generated in whole or in part. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable devices. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another, for example, the computer instructions can be transmitted from one website, computer, server or data center to another website, computer, server or data center through wireless or wired direction. The wired transmission mode includes optical fiber, twisted pair, coaxial cable, etc. The wireless transmission includes infrared, microwave, etc. The computer-readable storage medium can be any available medium that can be accessed by a computer or a data storage device such as a server, data center, etc. containing one or more available medium collections. The available medium can be a magnetic medium (for example, floppy disk, hard disk, magnetic tape), optical medium (for example, DVD), or semiconductor medium. The semiconductor medium can be a solid state disk.

[0050] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the above-described system, device and module can refer to the corresponding process in the foregoing method embodiments, which will not be described here.

[0051] In several embodiments provided in the present application, it should be understood that the disclosed system, device and method can be implemented in other ways. For example, the above-described device embodiments are only schematic, for example, the division of the modules is only a logical function division, and actual implementation can have another division manner, for example, a plurality of modules or components can be combined or integrated into another system, or some features can be ignored or not executed. In addition, the coupling or direct coupling or communication connection between the displayed or discussed each other can be indirect coupling or communication connection through some interfaces, devices or modules, which can be electrical, mechanical or other forms.

[0052] The modules described as separate components may or may not be physically separate, and the components shown as modules may or may not be physical modules, and may be located in one place or distributed on multiple network modules. Part or all of the modules can be selected according to actual needs to achieve the purpose of the embodiment of the present application.

[0053] In addition, the functional modules in each embodiment of the present application can be integrated into one processing module, or each module can exist physically alone, or two or more modules can be integrated into one module.

[0054] The functions, if realized in the form of software function modules and sold or used as independent products, can be stored in a computer readable storage medium. Based on this understanding, the technical solutions of the present application essentially or the part that contributes to the prior art or part of the technical solutions can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes a number of instructions for causing a computer device (which can be a personal computer, a server or a network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present application. The aforementioned storage medium includes: a U disk, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk, and various program code storage media.

[0055] The above is only a specific embodiment of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art can easily think of changes or replacements within the technical scope disclosed in the present application, which should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

[0056] Finally: the above is only the preferred embodiment of the present application, and is not used to limit the present application, any modification, equivalent replacement, improvement, etc. within the spirit and principles of the present application, should be included in the protection scope of the present application.

Claims

1. A wind turbine blade structural health monitoring method based on multi-field decoupled fiber sensing, characterized in that, Comprising: S1, synchronously collecting vibration signals, strain signals and temperature signals through an optical fiber sensor network; S2, extracting gear box vibration transmission path characteristics and raindrop impact noise spectrum characteristics to generate background noise feature vectors; S3, calculating the coherence function value of the vibration signal and the background noise feature vector in the frequency domain, identifying the frequency band with a coherence function value greater than a preset threshold as an interference dangerous frequency band; S4, analyzing the polarization angle distribution of the vibration signal in the interference dangerous frequency band, and inversely calculating the physical path topology of the background noise entering the blade according to the polarization principal axis deviation; S5, when the physical path topology passes through a pre-defined structure weak area, the damage weight coefficient of the strain signal in the corresponding interference dangerous frequency band is increased and a phase shift reconstruction is performed to generate an anti-interference strain feature vector; S6, calculating the energy dissipation rate of the anti-interference strain feature vector in the interference dangerous frequency band, and marking as a damage energy leakage event when the energy dissipation rate exceeds the theoretical threshold of the material; S7, judging the blade structure damage type and position according to the time domain correlation between the damage energy leakage event and the temperature signal.

2. The method for wind turbine blade structural health monitoring based on multi-field decoupling fiber-optic sensing according to claim 1, characterized in that, Synchronously collecting vibration signals, strain signals and temperature signals through an optical fiber sensor network, comprising: Arranging an optical fiber grating sensor array on the surface of the wind turbine blade to synchronously sample vibration signals, strain signals and temperature signals; Grouping the vibration signals according to the blade spanwise position, and extracting the amplitude variation characteristics of each group of vibration signals in the main vibration frequency band of the blade; Performing temperature drift compensation on the collected strain signals, and storing the compensated strain signals according to the blade cross-sectional area; Mapping the synchronously collected temperature signals to a three-dimensional coordinate grid of the blade to form a temperature gradient distribution field.

3. The method of claim 2, wherein the method further comprises: Extracting gear box vibration transmission path characteristics and raindrop impact noise spectrum characteristics to generate background noise feature vectors, comprising: Separating the path coupling components transmitted to the blade from the gear box base vibration signals, and extracting the frequency response function amplitude characteristics of the corresponding path coupling components; Analyzing the short-time spectrum envelope of the raindrop impact noise signal, and extracting the envelope peak frequency and its harmonic distribution characteristics; Aligning and fusing the frequency response function amplitude characteristics and the envelope peak frequency and its harmonic distribution characteristics according to the frequency band; Normalizing the fused characteristics to generate background noise feature vectors.

4. The method of claim 3, wherein the method further comprises: Calculating the coherence function value of the vibration signal and the background noise feature vector in the frequency domain, identifying the frequency band with a coherence function value greater than a preset threshold as an interference dangerous frequency band, comprising: Performing fast Fourier transform on the vibration signal to obtain the vibration signal spectrum; Converting the background noise feature vector into a reference spectrum according to the frequency band division; Calculating the coherence function value between the vibration signal spectrum and the reference spectrum to generate a coherence function spectrum; Locating the frequency interval with a coherence function value greater than a preset threshold in the coherence function spectrum; Merging adjacent frequency intervals with a coherence function value continuously greater than a preset threshold to form an interference dangerous frequency band.

Images