Wind turbine blade operation noise real-time acquisition method

Abstract

The invention discloses a wind turbine blade operation noise real-time acquisition method, which comprises the following steps of: firstly, dividing an airspace according to the rotation characteristic of a wind turbine blade and the characteristics of a beam former so as to ensure that the acquired blade operation noise has enough sound pressure level; secondly, combining angle information of a main shaft of the wind turbine and boundary information of the divided airspace blocks, and carrying out time domain segmentation on the collected sound signals; and finally, inputting an original sound signal obtained by segmenting each block into a corresponding broadband beam former so as to obtain a blade operation noise signal enhanced according to the blocks. The method not only effectively improves the extraction quality of the target signal, but also provides reliable data support for subsequent fault diagnosis, positioning and type identification.

Description

Technical Field

[0001] The invention belongs to the technical field of wind turbine fault detection, and in particular is a method for real-time collection of wind turbine blade operation noise. Technical Background

[0002] The working environment of wind turbines is highly complex, and the acquisition process of their blade operation noise signals is easily affected by a variety of interference noises, including wind noise with fluctuating characteristics, mechanical noise with periodic characteristics, biological noise with strong randomness, and continuous electrical equipment operation noise. Among them, the wind noise pressure level in the wind field environment is significantly higher than the wind turbine blade operation noise, resulting in the blade operation noise often being masked by wind noise, thereby increasing the difficulty of signal acquisition.

[0003] In the prior art, independent microphones are usually used to collect the operating noise of wind turbine blades, and high-pass filters are used to filter out sound signals with frequencies below 2000Hz. The theoretical basis of this method is that a large number of studies have shown that wind noise energy is mainly concentrated in the low-frequency area below 2000Hz. However, this method has obvious limitations: on the one hand, it cannot effectively filter out interference from other types of noise; on the other hand, since independent microphones cannot obtain the azimuth information of the sound source, it is difficult to accurately locate the fault. In recent years, with the popularization of microphone array technology, it has gradually been applied to the collection and analysis of noise signals. Existing methods are mainly divided into two categories: the first method uses a fixed-direction beamformer combined with a high-pass filter to continuously collect acoustic signals along the normal direction of the array geometric center, such as the document "Early warning of damaged wind turbine blades using spatial-temporal spectral analysis of acoustic emissions signals", which introduces a microphone array combined with a fixed main lobe direction beamforming to enhance weak signals and weaken environmental noise signals, and uses short-time Fourier transform (STFT) to create the time-spectrum of the processed signal, and calculates the radiation energy change period and cyclic modulation spectrum from the time-spectrum. Blade fault monitoring is performed based on whether there are periodic or cyclostationary features in the acoustic emission signal. Although this method suppresses non-target acoustic signals to a certain extent, it fails to give full play to the advantages of the beamformer, lacks full use of spatiotemporal information, and cannot effectively remove occasional biological noise, resulting in a large amount of useless information in the output signal, which increases the burden of data storage; the second type of method is based on compressed sensing technology for sound field reconstruction and sound source localization, such as the method proposed in the literature "Damage identification of wind turbine blades using an adaptive method for compressive beamforming based on the generalized minimax-concave penalty function", which uses a generalized minimax concave penalty function to improve sparse recovery capabilities and adaptively adjusts the calculation step size to achieve high spatial resolution acoustic inverse problem solution with less computing time. Experiments have shown that this method can accurately locate the fault position of wind turbine blades. This processing method loses the time-frequency characteristics of the acoustic signal and cannot provide an effective acoustic signal basis for fault type identification based on time-frequency feature extraction and data-driven. Summary of the invention

[0004] The present invention aims to solve the complex environmental challenges faced by wind turbine noise collection, including the problem of multiple noise interferences, especially wind noise masking the blade operation noise, and the limitations of existing technologies in noise removal and sound source localization. To this end, a method for real-time collection of wind turbine blade operation noise is proposed. This method is a new method that combines airspace division and broadband beamforming to improve the signal-to-noise ratio, reduce data redundancy, and provide high-quality data support for subsequent fault diagnosis and type identification, thereby improving the intelligent level of wind turbine operation and maintenance. The method relies on the self-developed 64-element microphone array and data acquisition card in hardware, and the blade rotation airspace division algorithm, signal segmentation algorithm and broadband beamformer in software to achieve the collection of high-quality wind turbine blade operation noise and include spatiotemporal information.

[0005] To achieve the above object, the present invention adopts the following technical solution: a method for real-time collection of wind turbine blade operation noise is designed, characterized in that the method comprises the following steps:

[0006] Step 1): Install a PCB with a uniform rectangular microphone array integrated on the upper surface on the lower section of a wind turbine tower, the wind turbine tower is located in the middle position behind the microphone array, the wind turbine blades and the microphone array are located on the same side of the wind turbine tower, and the wind turbine blades are located directly above the microphone array; the angle between the lower surface of the plane where the microphone array is located and the plumb line of the wind turbine tower is the microphone array installation angle, denoted as α; the cone angle of the wind turbine blade is denoted as β; the horizontal distance between the geometric center of the microphone array and the center of the hub is denoted as D; the distance from the center of the hub to the tip of the wind turbine blade is denoted as L B ; The height of the geometric center of the microphone array from the ground is denoted as H A ; The height from the center of the wheel hub to the ground is denoted as H B ;

[0007] Step 2): Set the coordinate system:

[0008] The center of the microphone array is taken as the origin O of the spatial rectangular coordinate system. The Z axis is perpendicular to the plane where the microphone array is located, and its positive direction points to the wind turbine blades; the X axis is on the plane where the microphone array is located, and is collinear with the symmetry axis of the microphone array that runs through the wind turbine tower, and its positive direction points to the ground; the Y axis passes through point O and is perpendicular to the XOZ plane, and its positive direction conforms to the right-hand coordinate system principle, thereby forming an XYZ coordinate system;

[0009] The hub center is taken as the origin O' of the spatial rectangular coordinate system. The Z' axis is collinear with the main shaft axis of the wind turbine, and its positive direction points to the front. The X' axis is on the plumb line passing through the O' point, and its positive direction points to the ground. The Y' axis passes through the O' point and is perpendicular to the X'O'Z' plane. The positive direction conforms to the right-hand coordinate system principle, thus forming the X'Y'Z' coordinate system.

[0010] Step 3): Divide the blade rotation airspace according to the wind turbine blade rotation characteristics and beamforming main lobe width change characteristics; the specific implementation is as follows:

[0011] Step 3-1) Determine the necessary physical parameters and spatial position parameter values, specifically the blade cone angle β, the microphone array installation angle α, and the height H of the microphone array geometric center from the ground A , Height of wheel hub center from ground H B , the horizontal distance D between the geometric center of the microphone array and the center of the hub and the length L from the blade tip to the center of the hub B ;

[0012] Step 3-2) According to the parameters determined in step 3-1), the airspace of the blade rotation airspace within the range of 45° on the left and right sides of the tower is divided to obtain the block boundary;

[0013] Step 4): Calculate the coordinates of the center of each block in the X'Y'Z' coordinate system;

[0014] Step 5): Calculate the expected main lobe angle value of the beamformer corresponding to each block according to the geometric center coordinates of each block;

[0015] Step 6): Design the FIR broadband beamformer corresponding to each block;

[0016] Step 6-1) Use the DCRCB algorithm or other robust beam design methods to calculate the expected main lobe beam response of each block when the signal frequency is 4KHz; the DCRCB optimization formula can be written as follows:

[0017]

[0018] ||p s || s =M

[0019] Where R is the data covariance matrix, using a simulated noise signal, which is 4KHz white noise; M is the number of array elements, where an array element refers to a microphone in a microphone array; ε0 is the maximum allowable steering vector error; P s represents the direction response vector of the true expected direction obtained according to the covariance matrix; The directional response vector representing the ideal desired direction; (·) H It means to find the conjugate transpose of a matrix.

[0020] Then the weight of the sound signal collected by each microphone in a certain block is:

[0021]

[0022] Then the expected main lobe beam response of a certain block can be obtained as:

[0023]

[0024] Repeat this step until the expected main lobe beam responses of all blocks are calculated;

[0025] Step 6-2) For the frequency range of 2 to 8 kHz of the wind turbine blade operation noise, divide the sub-bands into sub-bands every 100 Hz, and calculate the weights of the constant main lobe response narrowband beamformer at the center frequency of each sub-band; the constant main lobe response narrowband beam design at the center frequency of the sub-band to which a certain block belongs can be written as the following optimization formula:

[0026]

[0027] subject to|ω H p(Ψ i )|≤ξ 0i ,Ψ i ∈Θ SL ,i=1,…,N SL

[0028] ||ω|| 2 ≤ξ0

[0029] Among them, Θ ML Represents the angle value within the main lobe range, Θ SL Indicates the angle value within the sidelobe range; N ML Indicates the total number of main lobe angle intervals that are discretized; N SL Represents the total number of discretized sidelobe angle intervals; Ψ j Represents the jth main lobe angle vector, including the pitch angle and azimuth angle; Ψ i represents the i-th sidelobe angle vector, including the pitch angle and azimuth angle; B d (Ψ j ) indicates that the signal direction is Ψ j The beam response vector at time ; P(Ψ j ) indicates that the signal direction is Ψ j is the steering vector at ; ω represents the weight of the narrowband beamformer; ξ 0i represents the maximum beam response corresponding to the i-th discrete sidelobe angle value; ξ0 represents the maximum total energy of the beamformer weight vector;

[0030] Repeat this step until the narrowband beamformer weights of each frequency sub-band to which each block belongs are calculated;

[0031] Step 6-3) Calculate the expected frequency response of the FIR beamformer to each block;

[0032]

[0033] in, represents the beamformer weight corresponding to the mth array element under the condition that the center frequency of the kth frequency subband and the main lobe points to the center point of the ith block, and takes the conjugate of the weight; T m represents the pre-delay time required for the mth array element; f k represents the center frequency of the kth frequency subband;

[0034] By using the above formula, the expected frequency response of the FIR beamformer corresponding to the mth array element at the center frequency of the kth sub-band is calculated, and this step is repeated until the expected frequency response of the FIR beamformer belonging to each block is calculated;

[0035] Step 6-4) Calculate the FIR broadband beamformer parameters of each block;

[0036] The following optimization formula is proposed to solve the optimal parameter length;

[0037]

[0038] Subject to 0≤L≤100

[0039] Among them, RMSE PB is the FIR filter passband RMS error, max[H c (f p )] represents the maximum frequency response within the stopband frequency range of the FIR filter; L represents the parameter length of the FIR filter; F SB Indicates the passband range of the FIR filter;

[0040] After selecting the optimal parameter length, the parameter solution method of the FIR broadband beamformer belonging to a certain block is expressed as the following optimization formula:

[0041]

[0042] Among them, λ k and λ p are the kth subband weight in the passband range and the pth subband weight in the stopband range, which are generally 1; e(f) = [1, exp(-i2πfT s ),…,exp(-i(L-1)2πfT s )] T is the array steering vector at frequency f, T s Indicates the sampling period; represents the maximum value of the sum of squares of stopband errors; h represents the FIR filter coefficient; H d (f) represents the desired beam response at frequency f;

[0043] Repeat this step until the FIR broadband beamformer parameters of each block are calculated;

[0044] Step 7): synchronously collect the aerodynamic noise signal and angle information of the wind turbine blades in operation; use a microphone array to collect the sound signal for 10 seconds through a data acquisition card to ensure that the wind turbine blades pass through all blocks;

[0045] Step 8): Pre-delay the acoustic signal received by each array element; delay time T m Calculated by the following formula:

[0046] T m = -int[τ m (Ψ0) / T s ]·T s

[0047] Among them, τ m (Ψ0) represents the desired pitch angle γ0 and azimuth angle The signal receiving delay time of the mth array element under the condition is calculated as follows:

[0048]

[0049] Among them, x m ,y m , z m represents the x, y, z coordinates of the mth array element in the XYZ coordinate system, c represents the speed of sound, which is generally 340 m / s;

[0050] Step 9): Segment the acoustic signal collected by each array element by combining the block meridian angle position information and the blade angle information to correspond the acoustic signal to the block one by one; the specific method is: first find the time point when the meridian angle of each block is equal to the blade angle, then extract the acoustic signal segment of the time period between the two time points and attribute it to the corresponding block; finally, input the acoustic signal segment of each block to the FIR broadband beamformer to which the block belongs, and obtain the single-passband acoustic signal belonging to the block.

[0051] Compared with the prior art, the present invention has the following beneficial effects:

[0052] 1) The present invention combines the rotation characteristics of wind turbine blades and the relationship between the main lobe width and azimuth of the beamformer to reasonably divide the blade rotation airspace, ensuring that the collected blade operation noise has a sufficiently high sound pressure level, providing a high-quality data basis for subsequent signal processing. In view of the frequency characteristics of wind turbine blade operation noise, a broadband beamformer with a passband of 2000 to 8000 Hz is designed to effectively filter out wind noise, electrical equipment operation noise and other high-frequency noise interference, significantly improve the target signal-to-noise ratio, and enhance the ability to extract fault features.

[0053] 2) The present invention uses block division and main shaft angle information to perform time domain segmentation on the collected acoustic signals, realizes spatial-temporal joint processing of noise signals, and improves signal processing efficiency. The acoustic signals after each block segmentation are processed using a broadband beamformer optimized for collecting wind turbine blade operation noise to obtain an enhanced blade operation noise signal, providing clear and accurate signal data for subsequent fault diagnosis and type identification.

[0054] 3) The present invention effectively reduces the collection and storage of useless signals and reduces the burden of data storage through spatial domain division and time domain segmentation. The present invention effectively reduces the collection and storage of useless signals and reduces the burden of data storage through spatial domain division and time domain segmentation. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 It is a front view schematic diagram of the coordinate system in step 2) of an embodiment of a method for real-time collection of wind turbine blade operating noise of the present invention.

[0056] Figure 2 It is a side view schematic diagram of the coordinate system in step 2) of an embodiment of a method for real-time collection of wind turbine blade operation noise of the present invention.

[0057] Figure 3 This is a schematic diagram of acoustic signal segmentation in step 9 of an embodiment of a method for real-time collection of wind turbine blade operation noise of the present invention.

[0058] Figure 4 This is a schematic diagram of the location of the block marked as ① in an embodiment of a method for real-time collection of wind turbine blade operating noise of the present invention.

[0059] Figure 5 This is a time domain diagram of the acoustic signal actually collected by element No. 1 of a microphone array in one embodiment of a method for real-time collection of wind turbine blade operation noise of the present invention.

[0060] Figure 6 This is a time-frequency diagram actually collected by element 1 of a microphone array in an embodiment of a method for real-time collection of wind turbine blade operation noise of the present invention.

[0061] Figure 7 This is a time domain diagram of a single-channel acoustic signal obtained from a block marked as ① in an embodiment of a method for real-time collection of wind turbine blade operation noise of the present invention.

[0062] Figure 8 This is a time-frequency diagram of a single-channel acoustic signal obtained from a block marked as ① in an embodiment of a method for real-time acquisition of wind turbine blade operation noise of the present invention. DETAILED DESCRIPTION

[0063] The technical solution of the present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0064] The present invention provides a method for real-time collection of wind turbine blade operation noise, the method comprising the following steps:

[0065] Step 1): Install the PCB with a uniform rectangular microphone array integrated on the upper surface on the lower section of the wind turbine tower. The wind turbine tower is located in the middle position behind the microphone array. The wind turbine blades and the microphone array are located on the same side of the wind turbine tower, and the wind turbine blades are located directly above the microphone array. The angle between the lower surface of the plane where the microphone array is located and the plumb line of the wind turbine tower is the microphone array installation angle, denoted as α. The cone angle of the wind turbine blade is denoted as β. The horizontal distance between the geometric center of the microphone array and the center of the hub is denoted as D. The distance from the center of the hub to the tip of the wind turbine blade is denoted as L B The height of the geometric center of the microphone array from the ground is denoted as H A The height from the center of the wheel to the ground is denoted as H B .

[0066] The microphone array is a uniform rectangular array composed of 64 silicon microphones (silicon microphones mainly use the deformation of silicon chips during vibration to generate voltage and convert sound signals into electrical signals). The 64 silicon microphones are connected to a self-developed 64-channel synchronous data acquisition card to collect signals from each silicon microphone respectively, with a synchronous acquisition frequency of 200KHz.

[0067] Step 2): Set the coordinate system

[0068] The center of the microphone array is taken as the origin O of the spatial rectangular coordinate system, the Z axis is perpendicular to the plane where the microphone array is located, and its positive direction points to the wind turbine blades; the X axis is on the plane where the microphone array is located, and is collinear with the symmetry axis of the microphone array that runs through the wind turbine tower, and its positive direction points to the ground; the Y axis passes through point O and is perpendicular to the XOZ plane, and its positive direction conforms to the right-hand coordinate system principle, thereby forming an XYZ coordinate system.

[0069] The hub center is taken as the origin O' of the spatial rectangular coordinate system. The Z' axis is collinear with the axis of the wind turbine main shaft, and its positive direction points to the front (windward). The X' axis is on the plumb line passing through the O' point, and its positive direction points to the ground. The Y' axis passes through the O' point and is perpendicular to the X'O'Z' plane. The positive direction conforms to the right-hand coordinate system principle, thus forming the X'Y'Z' coordinate system.

[0070] Step 3): Divide the blade rotation airspace according to the wind turbine blade rotation characteristics and beamforming main lobe width change characteristics. The specific implementation is as follows:

[0071] Step 3-1) Determine the necessary physical parameters and spatial position parameter values, specifically the blade cone angle β, the microphone array installation angle α, and the height H of the microphone array geometric center from the ground A , Height of wheel hub center from ground H B , the horizontal distance D between the geometric center of the microphone array and the center of the hub and the length L from the blade tip to the center of the hub B The above parameter values ​​can be obtained through the wind turbine parameters and its installation status and the installation status of the microphone array.

[0072] Step 3-2) According to the parameters determined in step 3-1), the airspace of blade rotation within the range of 45° on the left and right sides of the tower is divided to obtain the block boundary. The specific implementation process is as follows:

[0073] Step 3-2-1) First determine the position of the intersection point A between the normal line of the center point of the plane where the microphone array is located and the blade rotation space, and then calculate the length L of the line connecting point A and the geometric center O' of the hub AO' .

[0074]

[0075] Step 3-2-2) Design a DCRCB beamformer when the target signal frequency is 4KHz and the signal direction is perpendicular to the microphone array plane, and calculate its main lobe width, that is, the maximum difference in pitch angle θ when the gain attenuation is -3dB max The maximum difference between the angle

[0076] Step 3-2-3) The latitude of the block to which point A belongs can be expressed as:

[0077]

[0078] The angle between the meridian of the block to which point A belongs and the positive direction of the X axis is and R3 represents the latitude boundary position of the block to which point A belongs, which is close to the blade root side, and R4 represents the latitude boundary position of the block to which point A belongs, which is close to the blade tip side.

[0079] Step 3-2-4) After obtaining the longitude and latitude positions of the block to which point A belongs, calculate the longitude and latitude positions of the remaining blocks to realize the division of the blade rotation airspace within 45° on the left and right sides of the tower.

[0080] The method for calculating the latitude position is as follows:

[0081]

[0082] Among them, n and m both start from 1 and increase by 1. nIt represents the distance from the nth latitude line to the hub center in the interval from the latitude line on the blade root side of the block to which point A belongs, R m It represents the distance from the mth latitude line to the center of the hub in the interval from the latitude line to the blade tip on the side of the block where point A belongs. Indicates the scaling factor corresponding to the spacing between two adjacent latitude lines. L1 represents the latitude spacing of the block to which point A belongs, and is calculated as follows:

[0083]

[0084] In order to avoid unreasonable latitudes near the blade root and blade tip, the following judgments need to be made during each calculation:

[0085]

[0086] Wherein, μ represents the tolerance coefficient of the spacing between adjacent latitudes, μ∈(0,1]; if the latitude spacing of the n+1th or m+1th block is less than μ times the required spacing, the block will be merged into the nth or mth block.

[0087] The method for calculating the meridian angle is as follows:

[0088]

[0089] in, k represents the number of blocks between the block to which the current latitude belongs and the block to which point A belongs in the X-axis direction. represents the scaling factor corresponding to the angle between two adjacent meridians, then Ω start Indicates the angle between the first meridians on both sides of the X-axis of the block group to which the current latitude belongs. Ω n Indicates the angle between the nth meridian on the side close to -45 degrees and the positive direction of the X axis. Ω m It represents the angle between the mth meridian on the +45 degree side and the positive direction of the X-axis.

[0090] Step 4): Calculate the coordinates of the center of each block in the X'Y'Z' coordinate system. The specific method is to take the middle value of the two latitude positions of each block to obtain R mid ; Then calculate the median value of the angle between the two meridians of each block to get Ω mid . The bold characters here represent row vectors. Considering that the airspace where the wind turbine blades rotate can be approximated as a conical surface, the calculation method for the center coordinates of each block is as follows:

[0091]

[0092] Among them, R mid ∈[0,L B ], x, y, z are the row vectors of the x, y, z coordinate values ​​of the center point of each block in the X'Y'Z' coordinate system.

[0093] Step 5): Calculate the expected main lobe angle value of the beamformer corresponding to each block according to the geometric center coordinates of each block. The specific steps are as follows:

[0094] Step 5-1) Convert the coordinates of the block's geometric center from the X'Y'Z' coordinate system to the XYZ coordinate system. The specific process is as follows:

[0095] First, the coordinate values ​​are translated.

[0096]

[0097] Among them, X1, Y1, and Z1 are row vectors of x, y, and z coordinate values ​​obtained after the center point of each block is translated in the X'Y'Z' coordinate system.

[0098] Secondly, the coordinate values ​​obtained by the translation processing are rotated to obtain their coordinate values ​​in the XYZ coordinate system.

[0099]

[0100] Among them, [·](:,i) represents the coordinate value vector of the center point of the i-th block, Represents the coordinate value vector of the center point of the i-th block obtained after the translation operation, It represents the coordinate value vector obtained by rotating the coordinate value of the center point of the i-th block around the Y axis by α degrees after translation. I represents the total number of blocks.

[0101] Step 5-2) Calculate the expected elevation angle and azimuth angle of the beamformer main lobe corresponding to each block according to the coordinate values ​​x2, y2, and z2 of the center point of each block in the XYZ coordinate system:

[0102]

[0103] Where θ represents the pitch angle, Indicates the direction angle.

[0104] Step 6): Design the FIR wideband beamformer corresponding to each block.

[0105] Step 6-1) Use the DCRCB algorithm or other robust beam design methods to calculate the expected main lobe beam response of each block when the signal frequency is 4KHz. The DCRCB optimization formula can be written as follows:

[0106]

[0107] ||ps || s =M

[0108] Wherein, R is the data covariance matrix. In the present invention, because it is essentially a fixed-direction beamformer application mode, an analog noise signal is used to generate it. The analog noise signal is a 4KHz white noise. M is the number of array elements. In the present invention, an array element refers to a microphone in a microphone array. The number of array elements is 64, that is, 64 silicon microphones. ε0 is the maximum allowable steering vector error. P s Represents the directional response vector of the true expected direction obtained from the covariance matrix. The directional response vector representing the ideal desired direction. (·) H It means to find the conjugate transpose of a matrix.

[0109] Then the weight of the sound signal collected by each microphone in a certain block is:

[0110]

[0111] Then the expected main lobe beam response of a certain block can be obtained as:

[0112]

[0113] Repeat this step until the expected main lobe beam responses of all blocks are calculated.

[0114] Step 6-2) For the wind turbine blade operating noise frequency range of 2 to 8 kHz, divide the sub-bands into sub-bands every 100 Hz, and calculate the constant main lobe response narrowband beamformer weights at the center frequency of each sub-band. The constant main lobe response narrowband beam design at the sub-band center frequency of a certain block can be written as the following optimization formula:

[0115]

[0116] subject to|ω H p(Ψ i )|≤ξ 0i ,Ψ i ∈Θ SL ,i=1,…,N SL

[0117] ||ω|| 2 ≤ξ0

[0118] Among them, Θ ML Represents the angle value within the main lobe range, Θ SL Indicates the angle value within the side lobe range. N ML Indicates the total number of main lobe angle intervals that are discretized. N SL Represents the total number of discretized sidelobe angle intervals.j Represents the jth main lobe angle vector, including the elevation angle and azimuth angle. i Represents the i-th sidelobe angle vector, including the pitch angle and azimuth angle. d (Ψ j ) indicates that the signal direction is Ψ j The beam response vector at time P(Ψ j ) indicates that the signal direction is Ψ j ω is the weight of the narrowband beamformer. ξ 0i Indicates the maximum beam response corresponding to the i-th discrete sidelobe angle value. For example, 0.01 means that the signal gain in this direction must be less than -40 dB. ξ0 indicates the maximum total energy of the beamformer weight vector, which is used to prevent the beamformer from being unstable due to excessive beamformer weights.

[0119] This step is repeated until the narrowband beamformer weights of each frequency sub-band to which each block belongs are calculated.

[0120] Step 6-3) Calculate the expected frequency response of the FIR beamformer to each block.

[0121]

[0122] in, represents the beamformer weight corresponding to the mth array element under the condition that the center frequency of the kth frequency subband and the main lobe points to the center point of the ith block, and takes the conjugate of the weight. m Indicates the pre-delay time required for the mth array element. k represents the center frequency of the kth frequency subband.

[0123] By using the above formula, the expected frequency response of the FIR beamformer corresponding to the mth array element at the kth sub-band center frequency is calculated, and this step is repeated until the expected frequency response of the FIR beamformer belonging to each block is calculated.

[0124] Step 6-4) Calculate the FIR wideband beamformer parameters for each block.

[0125] Since the calculation amount of FIR broadband beamformer parameter determination is large, it is necessary to select a reasonable parameter length to ensure a good beamforming effect while reducing the calculation amount. The following optimization formula is proposed to solve the optimal parameter length.

[0126]

[0127] Subject to 0≤L≤100

[0128] Among them, RMSE PB is the FIR filter passband RMS error, max[Hc (f p )] represents the maximum frequency response within the stopband frequency range of the FIR filter. L represents the parameter length of the FIR filter. SB Indicates the passband range of the FIR filter.

[0129] After selecting the optimal parameter length, the parameter solution method of the FIR broadband beamformer belonging to a certain block can be expressed as the following optimization formula:

[0130]

[0131] Among them, λ k and λ p They are the kth subband weight in the passband and the pth subband weight in the stopband, respectively, and are generally 1. e(f)=[1,exp(-i2πfT s ),…,exp(-i(L-1)2πfT s )] T is the array steering vector at frequency f, T s Indicates the sampling period. represents the maximum value of the sum of squares of the stopband errors. h represents the FIR filter coefficient. d (f) represents the desired beam response at frequency f.

[0132] This step is repeated until the FIR wideband beamformer parameters of each block are calculated.

[0133] Step 7): Synchronously collect the aerodynamic noise signal and angle information of the wind turbine blades in operation. Use a microphone array to collect the acoustic signal for 10 seconds through a data acquisition card to ensure that the wind turbine blades pass through all blocks. The microphone array must be consistent with the array element position when designing the beamformer in the above step.

[0134] The angle information of the wind turbine blades can be obtained by a variety of methods, as long as it is synchronized with the acoustic signal in time. The present invention proposes three angle acquisition methods.

[0135] Method 1: Use angle sensors combined with transmission transposition to obtain real-time angle information of wind turbine blades, and share a synchronous data acquisition card with the microphone array to calculate the real-time angle of each blade;

[0136] Method 2: Pre-delay the 10-second acoustic signal received by a certain array element, then perform short-time Fourier transform, calculate the root mean square energy of each time segment, find the position of the local maximum root mean square energy on the time axis, calculate the time difference between each two local maximum values, record it as ΔT, and record the longest time required for the wind turbine to rotate once as T. max and the shortest duration is T min , first determine whether the time difference is greater than and less than If the time difference is not obtained, it is retained and otherwise eliminated. Then, outliers are eliminated according to the interquartile range (IQR) method. Finally, the average of the retained time differences is calculated as the time required for the blade to rotate 120°. The speed of the wind turbine blade during this period can be further calculated to obtain the angle information. This method cannot achieve the same real-time and authenticity as method 1.

[0137] Method three, combining image information to obtain wind turbine blade angle information. Currently, image processing is widely and maturely used in wind turbine blade identification. The real-time blade angle provided by image recognition can ensure a certain degree of real-time performance.

[0138] Step 8): Pre-delay the acoustic signal received by each array element. Delay time T m Calculated by the following formula:

[0139] T m = -int[τ m (Ψ0) / T s ]·T s

[0140] Among them, τ m (Ψ0) represents the desired pitch angle θ0 and azimuth angle The signal receiving delay time of the mth array element under the condition is calculated as follows:

[0141]

[0142] Among them, x m ,y m , z m represents the x, y, z coordinate values ​​of the mth array element in the XYZ coordinate system, and c represents the speed of sound, which is generally 340m / s.

[0143] Step 9): Segment the acoustic signals collected by each array element by combining the block meridian angle position information and the blade angle information to correspond the acoustic signals to the blocks one by one. The specific method is: first find the time point when the meridian angle of each block is equal to the blade angle, then extract the acoustic signal segment of the time period between the two time points and attribute it to the corresponding block. Finally, input the acoustic signal segment of each block to the FIR broadband beamformer to which the block belongs, and obtain the single-passband acoustic signal that finally belongs to the block, thereby obtaining high-quality wind turbine blade operation noise.

[0144] The acoustic signal obtained through the above steps contains high-quality wind turbine operating noise and greatly suppresses other noise. Through one-to-one correspondence with the blocks, the acoustic signal contains spatial information, providing a solid foundation for subsequent fault diagnosis, positioning and identification.

[0145] The technical effect of the present invention is further explained below in conjunction with a certain wind farm application.

[0146] The audio frequency range that the data acquisition card can collect is 20 to 16000Hz. According to the characteristics of the wind turbine blade operating noise frequency, the passband frequency band of the broadband beamformer should be set between 2000 and 8000Hz. This frequency range includes the noise frequencies caused by three faults: blade root cracks, blade trailing edge cracks, and blade tip damage; the stopband frequency band should be set to 20 to 1800Hz and 8200 to 16000Hz to filter out wind noise, electrical equipment operating noise and other high-frequency noise; the transition band frequency band is 1800 to 2000Hz and 8000Hz to 8200Hz. The passband frequency band is divided into 60 sub-bands according to the 100Hz bandwidth to meet the narrowband signal requirements, the stopband frequency band is divided into 24 sub-bands according to 400Hz, and the transition frequency band is divided into 1 sub-band. The parameter length of the FIR filter is 51.

[0147] The operating noise of a wind turbine with a cracked blade root trailing edge was collected in a wind farm in Zhangjiakou according to the steps of the invention. Figure 5-6 This is the time domain diagram and time-frequency diagram of the acoustic signal actually collected by array element 1. Figure 7-8 The proposed method is followed by Figure 4 The time domain and time-frequency diagrams of the monophonic acoustic signal obtained in the block marked as ① in the figure. It can be seen that the background noise and the wind turbine blade fault acoustic signal within 2300 to 3000 Hz have been suppressed.

[0148] Any matters not described in the present invention are applicable to the prior art.

Claims

1. A method for real-time collection of wind turbine blade operating noise, characterized in that: The method comprises the following steps: Step 1): Install a PCB with a uniform rectangular microphone array integrated on the upper surface on the lower section of a wind turbine tower, the wind turbine tower is located in the middle position behind the microphone array, the wind turbine blades and the microphone array are located on the same side of the wind turbine tower, and the wind turbine blades are located directly above the microphone array; the angle between the lower surface of the plane where the microphone array is located and the plumb line of the wind turbine tower is the microphone array installation angle, denoted as α; the cone angle of the wind turbine blade is denoted as β; the horizontal distance between the geometric center of the microphone array and the center of the hub is denoted as D; the distance from the center of the hub to the tip of the wind turbine blade is denoted as L B ; The height of the geometric center of the microphone array from the ground is denoted as H A ; The height from the center of the wheel hub to the ground is denoted as H B ; Step 2): Set the coordinate system The center of the microphone array is taken as the origin O of the spatial rectangular coordinate system. The Z axis is perpendicular to the plane where the microphone array is located, and its positive direction points to the wind turbine blades; the X axis is on the plane where the microphone array is located, and is collinear with the symmetry axis of the microphone array that runs through the wind turbine tower, and its positive direction points to the ground; the Y axis passes through point O and is perpendicular to the XOZ plane, and its positive direction conforms to the right-hand coordinate system principle, thereby forming an XYZ coordinate system; The hub center is taken as the origin O' of the spatial rectangular coordinate system. The Z' axis is collinear with the main shaft axis of the wind turbine, and its positive direction points to the front. The X' axis is on the plumb line passing through the O' point, and its positive direction points to the ground. The Y' axis passes through the O' point and is perpendicular to the X'O'Z' plane. The positive direction conforms to the right-hand coordinate system principle, thus forming the X'Y'Z' coordinate system. Step 3): Divide the blade rotation airspace according to the wind turbine blade rotation characteristics and beamforming main lobe width change characteristics; the specific implementation is as follows: Step 3-1) Determine the necessary physical parameters and spatial position parameter values, specifically the blade cone angle β, the microphone array installation angle α, and the height H of the microphone array geometric center from the ground A , Height of wheel hub center from ground H B , the horizontal distance D between the geometric center of the microphone array and the center of the hub and the length L from the blade tip to the center of the hub B ; Step 3-2) According to the parameters determined in step 3-1), the airspace of the blade rotation airspace within the range of 45° on the left and right sides of the tower is divided to obtain the block boundary; Step 4): Calculate the coordinates of the center of each block in the X'Y'Z' coordinate system; Step 5): Calculate the expected main lobe angle value of the beamformer corresponding to each block according to the geometric center coordinates of each block; Step 6): Design the FIR broadband beamformer corresponding to each block; Step 6-1) Use the DCRCB algorithm or other robust beam design methods to calculate the expected main lobe beam response of each block when the signal frequency is 4KHz; the DCRCB optimization formula can be written as follows: ||p s || s =M Where R is the data covariance matrix, using a simulated noise signal, which is 4KHz white noise; M is the number of array elements, where an array element refers to a microphone in a microphone array; ε0 is the maximum allowable steering vector error; P s represents the direction response vector of the true expected direction obtained according to the covariance matrix; The directional response vector representing the ideal desired direction; (·) H It means to find the conjugate transpose of a matrix. Then the weight of the sound signal collected by each microphone in a certain block is: Then the expected main lobe beam response of a certain block can be obtained as: Repeat this step until the expected main lobe beam responses of all blocks are calculated; Step 6-2) For the frequency range of 2 to 8 kHz of the wind turbine blade operation noise, divide the sub-bands into sub-bands every 100 Hz, and calculate the weights of the constant main lobe response narrowband beamformer at the center frequency of each sub-band; the constant main lobe response narrowband beam design at the center frequency of the sub-band to which a certain block belongs can be written as the following optimization formula: subject to|ω H p(Ψ i )|≤ξ 0i ,Ψ i ∈Θ SL ,i=1,…,N SL ||ω|| 2 ≤ξ0 Among them, Θ ML Represents the angle value within the main lobe range, Θ SL Indicates the angle value within the sidelobe range; N ML Indicates the total number of main lobe angle intervals that are discretized; N SL Represents the total number of discretized sidelobe angle intervals; Ψ j Represents the jth main lobe angle vector, including the pitch angle and azimuth angle; Ψ i represents the i-th sidelobe angle vector, including the pitch angle and azimuth angle; B d (Ψ j ) indicates that the signal direction is Ψ j The beam response vector at time ; P(Ψ j ) indicates that the signal direction is Ψ j is the steering vector at ; ω represents the weight of the narrowband beamformer; ξ 0i represents the maximum beam response corresponding to the i-th discrete sidelobe angle value; ξ0 represents the maximum total energy of the beamformer weight vector; Repeat this step until the narrowband beamformer weights of each frequency sub-band to which each block belongs are calculated; Step 6-3) Calculate the expected frequency response of the FIR beamformer to each block; in, represents the beamformer weight corresponding to the mth array element under the condition that the center frequency of the kth frequency subband and the main lobe points to the center point of the ith block, and takes the conjugate of the weight; T m represents the pre-delay time required for the mth array element; f k represents the center frequency of the kth frequency subband; By using the above formula, the expected frequency response of the FIR beamformer corresponding to the mth array element at the center frequency of the kth sub-band is calculated, and this step is repeated until the expected frequency response of the FIR beamformer belonging to each block is calculated; Step 6-4) Calculate the FIR broadband beamformer parameters of each block; The following optimization formula is proposed to solve the optimal parameter length; Subject to 0≤L≤100 Among them, RMSE PB is the FIR filter passband RMS error, max[H c (f p )] represents the maximum frequency response within the stopband frequency range of the FIR filter; L represents the parameter length of the FIR filter; F SB Indicates the passband range of the FIR filter; After selecting the optimal parameter length, the parameter solution method of the FIR broadband beamformer belonging to a certain block is expressed as the following optimization formula: Among them, λ k and λ p are the kth subband weight in the passband and the pth subband weight in the stopband respectively; e(f) = [1, exp(-i2πfT s ),…,exp(-i(L-1)2πfT s )] T is the array steering vector at frequency f, T s Indicates the sampling period; represents the maximum value of the sum of squares of stopband errors; h represents the FIR filter coefficient; H d (f) represents the desired beam response at frequency f; Repeat this step until the FIR broadband beamformer parameters of each block are calculated; Step 7): synchronously collect the running noise signal of the wind turbine blades and the angle information of the wind turbine blades; use the microphone array to collect the sound signal for 10 seconds through the data acquisition card to ensure that the wind turbine blades pass through all blocks; Step 8): Pre-delay the acoustic signal received by each array element; delay time T m Calculated by the following formula: T m =-int[τ m (Ψ0) / T s ]·T s Among them, τ m (Ψ0) represents the desired pitch angle θ0 and azimuth angle The signal receiving delay time of the mth array element under the condition is calculated as follows: Among them, x m ,y m , z m represents the x, y, z coordinates of the mth array element in the XYZ coordinate system, and c represents the speed of sound; Step 9): Segment the acoustic signal collected by each array element by combining the block meridian angle position information and the blade angle information to correspond the acoustic signal to the block one by one; the specific method is: first find the time point when the meridian angle of each block is equal to the blade angle, then extract the acoustic signal segment of the time period between the two time points and attribute it to the corresponding block; finally, input the acoustic signal segment of each block to the FIR broadband beamformer to which the block belongs, and obtain the single-passband acoustic signal belonging to the block.

2. A method for real-time collection of wind turbine blade operating noise according to claim 1, characterized in that: In step 1), the microphone array is a uniform rectangular array consisting of 64 silicon microphones.

3. A method for real-time collection of wind turbine blade operating noise according to claim 2, characterized in that: The 64 silicon microphones are connected to a self-developed 64-channel synchronous data acquisition card to collect signals from each silicon microphone respectively, and the synchronous acquisition frequency is 200KHz.

4. A method for real-time collection of wind turbine blade operating noise according to claim 1, characterized in that: The specific implementation process of step 3-2) is as follows: Step 3-2-1) First determine the position of the intersection point A between the normal line of the center point of the plane where the microphone array is located and the blade rotation space, and then calculate the length L of the line connecting point A and the geometric center O' of the hub AO' ; Step 3-2-2) Design a DCRCB beamformer when the target signal frequency is 4KHz and the signal direction is perpendicular to the microphone array plane, and calculate its main lobe width, that is, the maximum difference in pitch angle θ when the gain attenuation is -3dB max The maximum difference between the angle Step 3-2-3) The latitude of the block to which point A belongs can be expressed as: The angle between the meridian of the block to which point A belongs and the positive direction of the X axis is and R3 represents the latitude boundary position of the block to which point A belongs, close to the blade root side, and R4 represents the latitude boundary position of the block to which point A belongs, close to the blade tip side; Step 3-2-4) After obtaining the longitude and latitude positions of the block to which point A belongs, calculate the longitude and latitude positions of the remaining blocks to realize the division of the airspace of blade rotation within the range of 45° on the left and right sides of the tower; The method for calculating the latitude position is as follows: Among them, n and m both start from 1 and increase by 1. n It represents the distance from the nth latitude line to the hub center in the interval from the latitude line on the blade root side of the block to which point A belongs, R m It indicates the distance from the mth latitude line to the hub center in the interval from the latitude line to the blade tip on the side of the block where point A belongs; Indicates the scaling factor corresponding to the spacing between two adjacent latitude lines; L1 indicates the latitude spacing of the block to which point A belongs, and is calculated as follows: In order to avoid unreasonable latitudes near the blade root and blade tip, the following judgments need to be made during each calculation: Wherein, μ represents the tolerance coefficient of the spacing between adjacent latitudes, μ∈(0,1]; if the spacing between the latitudes of the n+1th or m+1th block is less than μ times the required spacing, the block will be merged into the nth or mth block; The method for calculating the meridian angle is as follows: in, k represents the number of blocks between the block to which the current latitude belongs and the block to which point A belongs in the X-axis direction. represents the scaling factor corresponding to the angle between two adjacent meridians, then Ω start Indicates the angle between the first meridians on both sides of the X-axis of the block group to which the current latitude belongs; Ω n Indicates the angle between the nth meridian on the side close to -45 degrees and the positive direction of the X axis; Ω m It represents the angle between the mth meridian on the +45 degree side and the positive direction of the X-axis.

5. The method for real-time collection of wind turbine blade operating noise according to claim 1, characterized in that: The specific method of step 4) is to obtain R by taking the middle value of the two latitude positions of each block. mid ; Then calculate the median value of the angle between the two meridians of each block to get Ω mid ; The bold characters here represent row vectors; considering that the rotating airspace of wind turbine blades can be approximated as a conical surface, the calculation method of the center coordinates of each block is as follows: Among them, R mid ∈[0,L B ], x, y, z are the row vectors of the x, y, z coordinate values ​​of the center point of each block in the X'Y'Z' coordinate system.

6. A method for real-time collection of wind turbine blade operating noise according to claim 1, characterized in that: The specific steps of step 5) are as follows: Step 5-1) Convert the coordinates of the block geometric center from the X'Y'Z' coordinate system to the XYZ coordinate system. The specific process is as follows: First, the coordinate values ​​are translated; Among them, X1, Y1, and Z1 are the row vectors of the x, y, and z coordinate values ​​of the center points of each block after translation in the X'Y'Z' coordinate system; Secondly, the coordinate values ​​obtained by the translation processing are rotated to obtain their coordinate values ​​in the XYZ coordinate system; Among them, [·](:,i) represents the coordinate value vector of the center point of the i-th block, Represents the coordinate value vector of the center point of the i-th block obtained after the translation operation, represents the coordinate value vector obtained by rotating the coordinate value of the center point of the i-th block after translation around the Y axis by α degrees; I represents the total number of blocks; Step 5-2) Calculate the expected elevation angle and azimuth angle of the beamformer main lobe corresponding to each block according to the coordinate values ​​x2, y2, and z2 of the center point of each block in the XYZ coordinate system: Where θ represents the pitch angle, Indicates the direction angle.

7. A method for real-time collection of wind turbine blade operating noise according to claim 1, characterized in that: In step 9), the wind turbine blade angle information can be obtained by any of the following methods: Method 1: Use angle sensors combined with transmission transposition to obtain real-time angle information of wind turbine blades, and share a synchronous data acquisition card with the microphone array to calculate the real-time angle of each blade; Method 2: Pre-delay the 10-second acoustic signal received by a certain array element, then perform short-time Fourier transform, calculate the RMS energy of each time segment, find the position of the local maximum RMS energy on the time axis, calculate the time difference between each two local maximum values ​​and record it as ΔT, and record the longest time required for the wind turbine to rotate once as T. max and the shortest duration is T min , first determine whether the time difference is greater than and less than If the time difference is not reached, it is retained; otherwise, it is eliminated. Then, outliers are eliminated according to the interquartile range method. Finally, the average value of the retained time difference is calculated as the time required for the blade to rotate 120°. The speed of the wind turbine blade during this period can be further calculated to obtain the angle information. Method 3: Combine image information to obtain wind turbine blade angle information.

Images