Acoustic imaging method for wind turbine generator blade fault diagnosis

Abstract

The invention discloses a wind turbine generator blade fault diagnosis-oriented acoustic imaging method, which comprises the following steps of: arranging a uniform circular microphone array to collect the vibration sound of a fan blade, and decomposing a received original sound signal by adopting a CEEMDAN algorithm to obtain an IMF component; calculating energy weighted synthesis kurtosis of the IMF components, screening the IMF components according to the energy weighted synthesis kurtosis, calculating an azimuth angle and a pitch angle of an IMF component array signal, and calculating beam forming output power according to the azimuth angle and the pitch angle; and performing maximum value normalization on the beam forming output power, and mapping the normalized power value into a gray value of the image, thereby obtaining an acoustic image. According to the method, the kurtosis index is synthesized by adopting complete adaptive noise set empirical mode decomposition and energy weighting, fault features are effectively extracted, and the accuracy of fault diagnosis can be improved.

Description

Technical Field

[0001] The present invention relates to the technical fields of wind power generation equipment condition monitoring technology, array signal processing, and acoustic imaging technology, and particularly relates to an acoustic imaging method for blade fault diagnosis of wind turbines. Background Art

[0002] With the growth of wind power production demand, wind turbines are developing towards high power and large scale, and the structural complexity is increasing. They are often deployed in harsh environments and subjected to extreme meteorological conditions and dynamic loads, resulting in an increased risk of blade damage. Given that the operation and maintenance cost accounts for 65%-95% of the total system investment, and the unplanned maintenance caused by faults accounts for 30%-60%, the development of efficient condition monitoring and fault diagnosis technologies has significant economic value.

[0003] The current mainstream fault diagnosis methods such as vibration sensors and acoustic emission have significant shortcomings. For example, vibration sensors need to be embedded inside the structure and cannot monitor external components; the diagnosis dimension is single, and traditional acoustics relies on spectrum analysis and lacks the ability to extract spatio-temporal joint features; the outputs of existing technologies are mostly time-domain curves / spectrum diagrams, making it difficult to achieve fault space mapping.

[0004] In recent years, the acoustic array technology (i.e., microphone array) has shown potential in sound source localization and sound field imaging. For example, an acoustic image is generated through beamforming. However, the existing applications of acoustic arrays are mostly limited to sound source localization, and the visual features of acoustic images are not fully utilized for condition recognition. At the same time, image recognition technologies (such as deep learning) have made significant progress in industrial inspection, but they have not been combined with the images generated by acoustic arrays for wind turbine monitoring.

[0005] Array signal processing technology realizes sound source inversion and sound field reconstruction through the spatial topology of multiple microphones. Its spatial domain filtering characteristics can effectively separate overlapping sound sources and enhance target features, overcoming the positioning ambiguity defect of traditional acoustic diagnosis. Compared with vibration / single-channel acoustic monitoring, the microphone array has the advantage of non-contact measurement and can achieve millimeter-level positioning of spatially distributed faults (such as blade aerodynamic abnormal noises), providing a new paradigm for visual fault diagnosis.

[0006] The inventor proposed a method in the patent application with the application number 201911295468.7 and the title "Method and Equipment for Fault Monitoring of Wind Turbine Blades Based on Acoustic Sensor Array". This method first fixedly installs an acoustic sensor array at the bottom of the support rod of the wind turbine for non-contact acquisition of blade sound signals. The next step is to perform enhancement processing on the acquired sound signals and construct a convolutional neural network (CNN) fault judgment model. Then, fault location is carried out for the faulty blade, including determining the fault area and determining the fault location through the fault azimuth angle. Then, the strong interference direction is suppressed, and then the fault location step is returned to accurately locate the faulty wind blade. Finally, a blade fault type database is established, and incremental training of the CNN fault judgment model is carried out. Although this method has a certain robustness to environmental noise, its adaptability under extreme weather conditions (such as strong wind and heavy rain) may still need to be improved.

[0007] In the article titled "Fault Detection of Wind Turbine Blades Based on Sound Signals", the author proposed a method for fault detection of wind turbine blades based on sound signals. It mainly uses the discrete Fourier transform (DFT) to extract the frequency domain features of the sound signals and uses a multi-layer perceptron (MLP) for classification. First, the acquired sound data is segmented into segments with a length of 4096 for subsequent processing and analysis. Then, the discrete Fourier transform is performed on the sound signals to extract the frequency domain features. Based on the extracted frequency domain features, a multi-layer perceptron is constructed for fault classification. This method only uses the discrete Fourier transform for feature extraction and needs to combine other signal processing techniques and feature extraction methods to improve the accuracy and robustness of fault diagnosis.

[0008] In the article titled "Research on Fault Monitoring Method of Wind Turbine Blades Based on Vibration", the author proposed a method for fault monitoring of wind turbine blades based on vibration signals. It mainly monitors the operating state of the blades by arranging vibration sensors in the inner cavity of the blades. The core steps of this method include data acquisition, data analysis, and fault diagnosis. The data acquisition system consists of a biaxial acceleration sensor, a wireless transmission device, a data acquisition unit, an optical fiber ring network, a field server, and analysis and online monitoring and fault diagnosis software. Its data acquisition strategy is mainly to install a biaxial composite vibration sensor inside the blade to monitor the vibration states such as the flap, pitch, and torsion of the blade, and a built-in temperature sensor is used to collect the temperature of the inner cavity surface of the blade. Data analysis includes two methods: time domain analysis and frequency domain analysis. Time domain analysis compares the trend of characteristic values of different blades, and frequency domain analysis extracts the natural frequency of the blade through Fourier transform and analyzes the frequency change. However, it needs to be customized and adjusted for different models or configurations of wind turbines, which limits its generalization application. Summary of the Invention

[0009] In view of this, the object of the present invention is to provide an acoustic imaging method for wind turbine blade fault diagnosis, so as to solve the technical problem of non-contact monitoring of wind turbine blade faults, realize the positioning and classification of blade faults, and improve the environmental adaptability, real-time performance and accuracy of the detection method.

[0010] The acoustic imaging method for wind turbine blade fault diagnosis of the present invention includes the following steps:

[0011] I) Arrange a uniform circular microphone array to collect the vibration sound of the wind turbine blade;

[0012] II) Decompose the original sound signal received by the microphone array using the CEEMDAN algorithm to obtain IMF components;

[0013] III) Calculate the energy-weighted composite kurtosis of the IMF components, and the calculation formula is as follows:

[0014] Q = ρ k ·K ni ·E k (1)

[0015] In the formula: Q is the energy-weighted composite kurtosis, ρ k is the correlation coefficient between the IMF component and the original sound signal, K ni is the kurtosis after normalization, E k is the energy ratio of each IMF component;

[0016] Among them, the calculation formula of ρ k is as follows:

[0017]

[0018] In the formula: C k (t) is the kth IMF component obtained after CEEMDAN decomposition, x(t) is the original sound signal, δ uk is the standard deviation of C k (t), which is used to measure the degree of dispersion of the data points of C k (t) relative to its mean value; δ x is the standard deviation of the original sound signal x(t), which is used to measure the degree of dispersion of the data points of x(t) relative to its mean value;

[0019] Among them, the calculation formula of K ni is as follows:

[0020]

[0021] In the formula, K(C vk ) is the kurtosis of the IMF component in the frequency domain, and its calculation formula is:

[0022]

[0023] where C vk (n) represents the value of the k-th IMF component of the v-th virtual array element at time point n; N is the total number of sampling points, i.e., the time series length of the signal; δ vk represents the standard deviation of the k-th IMF component of the v-th virtual array element;

[0024] where, E k is calculated as follows:

[0025]

[0026] where T represents the time length of the signal, i.e., the time range of integration;

[0027] IV) Screen the IMF components according to the energy-weighted synthetic kurtosis, retain the IMF components with an energy-weighted synthetic kurtosis greater than the set threshold, and recombine the retained IMF components by channel into an IMF component array signal;

[0028] V) Calculate the arrival azimuth angle θ and elevation angle of the IMF component array signal

[0029] VI) Calculate the beamforming output power according to the azimuth angle and elevation angle

[0030]

[0031] where, w m is the beamforming weight vector, R y is the signal covariance matrix, and H represents the conjugate transpose;

[0032] VII) Normalize the maximum value of the beamforming output power and map the normalized power value to the gray value of the image to obtain the acoustic image.

[0033] Furthermore, in step 2), the CEEMDAN algorithm is used to decompose the received signal of the microphone array, including:

[0034] (1) Add Gaussian white noise to the signal S(t) to be decomposed to obtain a new signal;

[0035] (2) Perform EMD decomposition on the new signal to obtain the first-order IMF component, denoted as IMF1(t), i.e.:

[0036]

[0037] where n is the number of IMF components obtained by decomposition, and N is the total number of decompositions in the whole process;

[0038] (3) After obtaining the first-order IMF component, calculate the residual noise r1(t), that is:

[0039] r1(t) = S(t) - IMF1(t) (8)

[0040] (4) Add a pair of positive and negative Gaussian white noises to r1(t) to obtain a new signal, and perform EMD decomposition on the new signal to obtain the second IMF component, denoted as IMF2(t);

[0041] (5) After obtaining the second-order IMF component, calculate the residual noise r2(t), that is:

[0042] r2(t) = r1(t) - IMF2(t) (9)

[0043] Repeat the above decomposition steps until the obtained residual signal cannot be further decomposed, that is, when the signal is a monotonic function, K IMF components are obtained, and the final residual signal is expressed as:

[0044]

[0045] where IMF k represents the k-th order IMF component obtained by decomposition using the CEEMDAN algorithm.

[0046] Furthermore, the calculation formulas for the azimuth angle θ and the elevation angle are as follows:

[0047]

[0048] where is the phase difference, which is obtained in the following way:

[0049] Map the array response vector of the IMF component array signal obtained in step IV) to the beam space:

[0050]

[0051] where represents the beam space manifold, is the array response vector of the IMF component array signal, and the expression is as follows:

[0052]

[0053] where γ m is the angle between the m-th microphone and the reference microphone. For a uniform circular array M is the number of microphones; k is the wave number, defined as k = 2π / λ, where λ is the wavelength of the signal; R is the radius of the microphone array circle, and j is the imaginary unit;

[0054] The calculation formula of

[0055]

[0056] where

[0057]

[0058] W is a matrix containing multiple beamforming weight vectors, and its row vectors maintain the orthogonality by selecting discrete azimuth angles to maintain of i d is the highest excitation mode order passing through the array aperture, v(α i ) is a complex vector that transforms the change of the azimuth angle α i into a Vandermonde-like structure; the vectors v(α i ) corresponding to multiple azimuth angles α

[0059] Vector transforms the azimuth angle change into a Vandermonde-like structure, expressed as

[0060]

[0061] v(α i ) is the value at the discrete azimuth angle α i , j is the imaginary unit, satisfying j 2 = -1;

[0062] C V V H is used to convert the elevation angle information of the array manifold into a symmetric amplitude taper, where:

[0063] C v = diag{j -d ,..., j -1 , j 0 , j 1 ,..., j d} (18)

[0064]

[0065] In the formula, w m is a weight vector; m is an index variable representing different weight vectors, m ∈ [-d, d]; N is the total number of signal sampling points;

[0066] The array is divided into two overlapping sub - arrays S (0) and S (1) , and the two satisfy:

[0067] S (1) = S (0) Ψ(21)

[0068] Estimate Ψ by the least - squares method or the total least - squares method, and the expression is:

[0069] Ψ=(S (0)H S (0) ) -1 S (0)H S (1) (22)

[0070] Perform eigenvalue decomposition on Ψ, and the phase difference of its eigenvalues is

[0071] Advantages of the present invention:

[0072] 1. The acoustic imaging method for blade fault diagnosis of wind turbine generators according to the present invention utilizes the omnidirectional signal receiving ability of the uniform circular array (UCA) and combines the adaptive beamforming technology to achieve precise positioning of blade faults.

[0073] 2. The acoustic imaging method for blade fault diagnosis of wind turbine generators according to the present invention adopts the complete adaptive noise ensemble empirical mode decomposition (CEEMDAN) and the energy - weighted synthetic kurtosis index to effectively extract fault features and improve the accuracy of fault diagnosis.

[0074] 3. The acoustic imaging method for blade fault diagnosis of wind turbine generators according to the present invention has good robustness to environmental changes and array errors, and can maintain stable imaging performance even in the presence of non - ideal factors such as multipath or mutual coupling between array elements. Description of the Drawings

[0075] Figure 1 is the UCA far - field model.

[0076] Figure 2 is the CEEMDAN algorithm flow chart. Detailed Embodiments

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

[0078] The acoustic imaging method for blade fault diagnosis of wind turbine generators in this embodiment includes the following steps:

[0079] I) Arrange a uniformly circular microphone array to collect the vibration sound of the fan blade.

[0080] When acoustically collecting the sound of the fan blade, the positions of the microphones need to be determined first. Multiple microphones are placed in a certain spatial topological structure, and data can be collected in parallel from different positions. Each microphone is regarded as an element of the array. By applying different phase weights between the microphones, the signal in the desired direction is enhanced, and the signals in other directions are suppressed.

[0081] In this embodiment, a uniformly circular array (UCA) microphone is used, which can achieve 360° omnidirectional signal reception and transmission. The circular symmetric structure makes its signal response consistent in any direction, avoiding the problem of uneven side lobes caused by geometric asymmetry in linear or rectangular arrays. This characteristic is particularly important in complex multipath environments or scenarios with randomly distributed interference sources, and can effectively reduce the direction-dependent error.

[0082] As Figure 1 shown, in this embodiment, M identical microphones are uniformly distributed on the circumference with a radius of R in the x-y plane. A narrowband plane wave with a wave number k = 2π / λ is transmitted from the sound source S on the fan blade. The origin O of the x-y plane coordinate system represents the center of the array, and the point S' is the projection point of S on the x-y plane. The radius of the reference element 0 is used as the reference line. Let θ measured downward from the z-axis be the elevation angle of the sound source arrival direction, and measured counterclockwise between OS and the reference line is the azimuth angle. The array response vector of the uniformly circular microphone array is expressed as

[0083]

[0084] where γ m is the angle between the m-th microphone and the reference microphone (usually the center of the array or the first microphone). For a uniformly circular array M is the number of microphones; k is the wave number, defined as k = 2π / λ, where λ is the wavelength of the signal; R is the radius of the microphone array circle, and j is the imaginary unit;

[0085] II) Use the CEEMDAN algorithm to decompose the original sound signal received by the microphone array to obtain IMF components. Using the CEEMDAN algorithm to decompose the received signal of the microphone array includes:

[0086] (1) Add Gaussian white noise to the signal S(t) to be decomposed to obtain a new signal;

[0087] (2) Perform EMD decomposition on the new signal to obtain the first-order IMF component, denoted as IMF1(t), that is:

[0088]

[0089] Wherein, n is the number of IMF components obtained by decomposition, and N is the total number of decompositions in the whole process;

[0090] (3) After obtaining the first-order IMF component, calculate the residual noise r1(t), that is:

[0091] r1(t) = S(t) - IMF1(t) (3)

[0092] (4) Add a pair of positive and negative Gaussian white noises to r1(t) to obtain a new signal, and perform EMD decomposition on the new signal to obtain the second IMF component, denoted as IMF2(t);

[0093] (5) After obtaining the second-order IMF component, calculate the residual noise r2(t), that is:

[0094] r2(t) = r1(t) - IMF2(t) (4)

[0095] Repeat the above decomposition steps until the obtained residual signal cannot be further decomposed, that is, when the signal is a monotonic function, K IMF components are obtained, and the final residual signal is expressed as:

[0096]

[0097] Wherein, IMF k represents the k-th order IMF component obtained by decomposition using the CEEMDAN algorithm.

[0098] The acoustic signal collected by the microphone array often contains environmental noise, and its amplitude and frequency change randomly with time. Traditional frequency-domain filtering methods are difficult to effectively separate such noise because they rely on global basis functions. In this embodiment, the CEEMDAN algorithm is used to decompose the acoustic signal collected by the microphone array,

[0099] and the signal is disassembled into multi-scale intrinsic mode functions (IMFs) through EMD adaptive decomposition, directly stripping high-frequency random noise and retaining low-frequency target sound source characteristics (such as blade crack acoustic emission signals, etc.).

[0100] In acoustic imaging, multi-source aliasing (such as the superposition of blade damage and wind noise, etc.) will cause beamforming positioning ambiguity. By decomposing the time-scale difference of the signal, EMD preferentially removes noise components that have nothing to do with the spatial distribution of the sound source (such as isotropic wind noise, etc.), enhances the consistency of the space-time characteristics of the target sound source, and improves the accuracy of sound source inversion.

[0101] The CEEMDAN algorithm decomposition solves the limitations of EMD by introducing specific noise at each stage of the decomposition and extracting each Intrinsic Mode Function (IMF), ultimately generating a unique residue. It overcomes problems such as the endpoint effect and mode aliasing. At the same time, this method also introduces signals that can handle positive and negative Gaussian white noise, thereby eliminating residual noise and improving computational efficiency.

[0102] III) Calculate the energy-weighted composite kurtosis of the IMF components, and the calculation formula is as follows:

[0103] Q = ρ k ·K ni ·E k (6)

[0104] In the formula: Q is the energy-weighted composite kurtosis, ρ k is the correlation coefficient between the IMF component and the original sound signal, K ni is the kurtosis after normalization, E k is the energy ratio of each IMF component.

[0105] Among them, the calculation formula of ρ k is as follows:

[0106]

[0107] In the formula: C k (t) is the k-th IMF component obtained after CEEMDAN decomposition, x(t) is the original sound signal, δ uk is the standard deviation of C k (t), which is used to measure the degree of dispersion of the data points of C k (t) relative to its mean value; δ x is the standard deviation of the original sound signal x(t), which is used to measure the degree of dispersion of the data points of x(t) relative to its mean value.

[0108] Among them, the calculation formula of K ni is as follows:

[0109]

[0110] In the formula, K(C vk ) is the kurtosis of the IMF component in the frequency domain, and its calculation formula is:

[0111]

[0112] In the formula, C vk (n) represents the value of the k-th IMF component of the v-th virtual array element at time point n; N is the total number of sampling points, that is, the length of the time series of the signal; δ vkis the standard deviation of the kth IMF component of the vth virtual array element.

[0113] Among them, E k The calculation formula is as follows:

[0114]

[0115] Where T represents the time length of the signal, that is, the time range of integration. T is the total time of signal sampling, which is used to define the time window for energy calculation.

[0116] IV) IMF components are screened according to the energy-weighted synthetic kurtosis, and the IMF components whose energy-weighted synthetic kurtosis is greater than a set threshold are retained, and the retained IMF components are reorganized into IMF component array signals by channel.

[0117] The vibration signal of the fan blade under normal operation can be approximately regarded as a normal distribution. When a fault impact signal is generated, the distribution of the signal amplitude shifts, and the signal no longer satisfies the normal distribution. The larger the energy-weighted composite kurtosis, the more impact components there are in the signal, indicating that there are more impact components in the signal and the more obvious the fault characteristics are.

[0118] The weighted score W of each IMF is calculated by the formula W = Km x Ek x correlation coefficient, which integrates energy, frequency domain characteristics and correlation. IMFs with high W values usually contain more effective information and may correspond to fault characteristics or key vibration modes, so they are retained first.

[0119] In specific screening, an empirical threshold (such as W>0.5) can be set to retain only IMFs above the threshold. By screening, IMFs with high energy, strong correlation with faults and significant frequency domain characteristics are retained, and noise or redundant components are eliminated.

[0120] V) Calculate the arrival azimuth θ and elevation angle of the IMF component array signal Azimuth angle θ and elevation angle The calculation formula is as follows:

[0121]

[0122] in is the phase difference, Obtained by:

[0123] Map the array response vector of the IMF component array signal obtained in step IV) to the beam space:

[0124]

[0125] In the formula, represents the beam space manifold, is the array response vector of the IMF component array signal.

[0126] The calculation formula is as follows:

[0127]

[0128] Where,

[0129]

[0130] W is a matrix containing multiple beamforming weight vectors. Each weight vector is used to weight the signals of different channels to achieve a specific beamforming effect. Its row vectors are orthogonal by selecting discrete azimuth angles to maintain orthogonality. d is the highest excitation mode order passing through the array aperture. v(α i ) is a complex vector that transforms the change of azimuth angle α i into a Vandermonde-like structure. The vectors v(α i ) corresponding to multiple azimuth angles α i are combined into the matrix W.

[0131] The vector transforms the change of azimuth angle into a Vandermonde-like structure, expressed as

[0132]

[0133] v(α i ) is the value at the discrete azimuth angle α i . j is the imaginary unit, satisfying j 2 = -1.

[0134] C V V H is used to convert the elevation information of the array manifold into a symmetric amplitude taper, where:

[0135] C v = diag{j -d ,..., j -1 , j 0 , j 1 ,..., j d} (17)

[0136]

[0137]

[0138] In the formula, w mis a weight vector; m is an index variable representing different weight vectors, m ∈ [-d, d]; N is the total number of signal sampling points.

[0139] For the one obtained by transforming through formula (13) the array is divided into two overlapping sub-arrays S (0) and S (1) , and the two satisfy:

[0140] S (1) = S (0) Ψ (20)

[0141] Estimate Ψ by the least squares method or the total least squares method, and the expression is:

[0142] Ψ = (S (0)H S (0) ) -1 S (0)H S (1) (21)

[0143] Perform eigenvalue decomposition on Ψ, and the phase difference of its obtained eigenvalues

[0144] VI) Calculate the beamforming output power according to the azimuth angle and the elevation angle

[0145]

[0146] where w m is the beamforming weight vector, R y is the signal covariance matrix, and H represents the conjugate transpose;

[0147] VII) Normalize the maximum value of the beamforming output power and map the normalized power value to the gray value of the image, thereby obtaining the acoustic image.

[0148] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that the technical solutions of the present invention can be modified or equivalently replaced without departing from the spirit and scope of the technical solutions of the present invention, and they should all be covered within the scope of the claims of the present invention.

Claims

1. An acoustic imaging method for blade fault diagnosis of wind turbines, characterized in that: Including the following steps: I) Arrange a uniform circular microphone array to collect the vibration sound of the fan blade; II) Decompose the original sound signal received by the microphone array using the CEEMDAN algorithm to obtain IMF components; III) Calculate the energy-weighted composite kurtosis of the IMF components, and the calculation formula is as follows: Q = ρ k ·K ni ·E k (1) Where: Q is the energy-weighted composite kurtosis, ρ k is the correlation coefficient between the IMF component and the original sound signal, K ni is the kurtosis after normalization, E k is the energy ratio of each IMF component; Among them, ρ k is calculated as follows: Where: C k (t) is the k-th IMF component obtained after CEEMDAN decomposition, x(t) is the original sound signal, and δ uk is the standard deviation of C k (t), which is used to measure the degree of dispersion of the data points of C k (t) relative to its mean; δ x is the standard deviation of the original sound signal x(t), which is used to measure the degree of dispersion of the data points of x(t) relative to its mean; Among them, K ni The calculation formula is as follows: where K(C vk ) is the kurtosis of the IMF component in the frequency domain, and its calculation formula is: where C vk (n) represents the value of the k-th IMF component of the v-th virtual array element at time point n; N is the total number of sampling points, that is, the time series length of the signal; δ vk represents the standard deviation of the k-th IMF component of the v-th virtual array element; Among them, E k The calculation formula is as follows: In the formula, T represents the time length of the signal, that is, the time range of integration; IV) Screen the IMF components according to the energy-weighted composite kurtosis, retain the IMF components with the energy-weighted composite kurtosis greater than the set threshold, and recombine the retained IMF components by channel into an IMF component array signal; V) Calculate the azimuth angle θ and elevation angle of the IMF component array signal VI) Calculate the beamforming output power based on the azimuth angle and elevation angle where, w m is the beamforming weight vector, R y is the signal covariance matrix, and H represents the conjugate transpose; VII) Normalize the beamforming output power Perform maximum normalization on it, map the normalized power value to the gray value of the image, and then obtain the acoustic image.

2. The acoustic imaging method for blade fault diagnosis of a wind turbine according to claim 1, characterized in that: In step 2), using the CEEMDAN algorithm to decompose the received signal of the microphone array includes: (1) Add Gaussian white noise to the signal S(t) to be decomposed to obtain a new signal; (2) Perform EMD decomposition on the new signal to obtain the first-order IMF component, denoted as IMF1(t), that is: In the formula, n is the number of IMF components obtained by decomposition, and N is the total number of decompositions in the whole process; (3) After obtaining the first-order IMF component, calculate the residual noise r1(t), that is: r1(t) = S(t) - IMF1(t) (8) (4) Add a pair of positive and negative Gaussian white noises to r1(t) to obtain a new signal, perform EMD decomposition on the new signal, and obtain the second IMF component, denoted as IMF2(t); (5) After obtaining the second-order IMF component, calculate the residual noise r2(t), that is: r2(t) = r1(t) - IMF2(t) (9) Repeat the above decomposition steps until the obtained residual signal cannot be further decomposed, that is, when the signal is a monotonic function, K IMF components are obtained, and the final residual signal is expressed as: where, IMF k represents the k-th IMF component obtained by decomposition using the CEEMDAN algorithm.

3. The acoustic imaging method for blade fault diagnosis of wind turbine according to claim 1, wherein: The azimuth angle θ and the elevation angle described in step V) are calculated as follows: wherein is the phase difference, is obtained by the following method: Map the array response vector of the IMF component array signal obtained in step IV) to the beam space: wherein, represents the beamspace manifold, is the array response vector of the IMF component array signal, and the expression is as follows: where γ m is the angle between the m-th microphone and the reference microphone. For a uniform circular array M is the number of microphones; k is the wave number, defined as k = 2π / λ, where λ is the wavelength of the signal; R is the radius of the circle of the microphone array, and j is the imaginary unit; The calculation formula is as follows: Where, W is a matrix containing multiple beamforming weight vectors, and its row vectors maintain orthogonality by selecting discrete azimuth angles to maintain . d is the highest excitation mode order passing through the array aperture, v(α i ) is a complex vector that transforms the change in azimuth angle α i into a Vandermonde-like structure; the vectors v(α i ) corresponding to multiple azimuth angles α i are combined into the matrix W; vector Convert the azimuth change into a Vandermonde-like structure, expressed as v(α i ) is the value at the discrete azimuth angle α i , where j is the imaginary unit, satisfying j 2 = -1; C V V H For converting the elevation angle information of an array manifold into a symmetric amplitude taper, where: C v = diag{j -d ,...,j -1 ,j 0 ,j 1 ,...,j d} (18) where w m is a weight vector; m is an index variable representing different weight vectors, m ∈ [-d, d]; N is the total number of signal sampling points; After being transformed by formula (13), the array is divided into two overlapping sub-arrays S (0) and S (1) , and the two satisfy: S (1) = S (0) Ψ(21) Estimate Ψ by the least squares method or the total least squares method, and the expression is: Ψ=(S (0)H S (0) ) -1 S (0)H S (1) (22) Perform eigenvalue decomposition on Ψ to obtain its eigenvalues The phase difference of

Images