A blade secondary icing de-icing method and system based on icing identification

Abstract

This invention relates to the field of blade de-icing technology, and discloses a method and system for de-icing blades under secondary icing based on icing recognition. The method includes acquiring initial surface images of the wind turbine blades, spatially clustering liquid water pixels to obtain water distribution areas; acquiring real-time blade temperature and real-time wind speed in the nacelle, and performing time-series recursion of ice-water states based on the water distribution areas to obtain dynamic evolution characteristics of the ice layer; estimating heat flux and calculating the heat dissipation gradient, and predicting ice layer thickness through icing rate integration; determining secondary icing risk areas or safe areas based on the predicted ice layer thickness and water distribution areas; when secondary icing risk areas exist, planning de-icing action parameters, generating de-icing action commands, and performing closed-loop position and temperature control on the de-icing actuator. This method can prevent equipment failures caused by secondary icing of blades.

Description

Technical Field

[0001] This invention relates to the field of blade de-icing technology, and in particular to a method and system for secondary de-icing of blades based on icing identification. Background Technology

[0002] Currently, during the construction and operation maintenance of offshore wind turbines in cold sea areas, blade icing exhibits different melting, flow, and recondensation states depending on temperature, wind speed, and the location of the blade's curvature. Moisture generated during de-icing tends to remain at the blade's leading edge, root, and recessed areas, forming secondary icing at low temperatures. This means that after the ice layer on the blade surface melts during de-icing, the residual moisture re-condenses into ice at low temperatures. To accurately adjust the de-icing range and intensity, monitoring and feedback control technologies capable of continuously sensing changes in ice thickness, spatial distribution, and residual moisture are required.

[0003] In one existing technology, a resistance heating film, heating cable, or conductive composite material is laid inside the leading edge of the blade, and a controller provides energization to different zones according to preset temperatures and heating durations. The heating elements convert electrical energy into heat, which is transferred through the blade shell to the icing interface, melting the ice or reducing the adhesion between the ice and the blade. Then, the centrifugal force and airflow generated by the blade rotation remove the loosened ice and meltwater, achieving the purpose of de-icing. However, this existing technology controls the internal heating elements of the blade to perform timed, zoned heating through a preset program, using heat conduction to melt the ice and centrifugal force to remove the meltwater. This results in the inability to continuously identify residual moisture and dynamic changes in the ice layer in the recessed areas, and the de-icing intensity cannot be adjusted according to the actual ice conditions.

[0004] Therefore, existing technologies are insufficient to effectively prevent equipment failures caused by secondary icing of blades. Summary of the Invention

[0005] This invention provides a method and system for de-icing blades after secondary icing based on icing identification, in order to prevent equipment failure caused by secondary icing of blades.

[0006] Firstly, in order to solve the above-mentioned technical problems, the present invention provides a method for de-icing secondary icing of blades based on icing identification, comprising: An initial surface image of the wind turbine blades is acquired, and the liquid water pixel spatial clustering is performed on the initial surface image to obtain the water distribution area; Real-time temperature data of the blade surface and real-time wind speed data of the nacelle anemometer are collected, and the ice-water state time sequence is recursively extrapolated by combining the moisture distribution area to obtain the dynamic evolution characteristics of the ice layer. Thermodynamic flux is estimated based on the dynamic evolution characteristics of the ice layer to obtain the initial heat flux. The heat conduction gradient is obtained based on the initial heat flux. The ice formation rate is then predicted based on the integral of the heat conduction gradient to obtain the predicted ice layer thickness. Based on the predicted ice thickness and the moisture distribution area, the risk area of ​​secondary icing is determined, and the risk area or safe area of ​​secondary icing is obtained. When the secondary icing risk area is obtained, the de-icing action parameters are planned according to the secondary icing risk area and the predicted ice thickness to obtain the de-icing action command, and the position and temperature closed-loop control of the de-icing actuator is performed according to the de-icing action command.

[0007] Secondly, the present invention provides a blade secondary icing and de-icing system based on icing identification, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the method described above.

[0008] Thirdly, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method described above.

[0009] Compared with the prior art, the present invention has the following beneficial effects: (1) This invention identifies liquid water pixels through grayscale mapping, local contrast calculation, and Sobel gradient extraction, and obtains the water distribution area by density spatial clustering based on Euclidean distance. Since the gradient magnitude can characterize the water trace boundary and the local contrast variance can reflect the grayscale consistency within the region, density spatial clustering can aggregate spatially adjacent liquid water pixels into a continuous region, solving the problem of difficulty in identifying the location of residual water in the leaf indentation and leading edge, realizing the quantitative characterization of the residual water distribution range, and providing a spatial basis for secondary icing monitoring.

[0010] (2) This invention queries the phase change rate of ice and water based on real-time temperature and wind speed, and uses the Markov state transition method to recursively deduce the change process between liquid water, ice-water mixture, and solid ice. Furthermore, it predicts the ice layer thickness through thermodynamic flux estimation, heat dissipation conduction gradient calculation, and icing rate integral. Since temperature and wind speed can reflect the phase change environment of the blade, the state transition matrix can continuously describe the change of water morphology, and the icing rate integral can accumulate the ice layer growth within the predicted time period, this invention solves the problem of the difficulty in continuously obtaining the dynamic evolution process of the ice layer after de-icing, and realizes the time-series prediction of secondary icing state and ice layer thickness.

[0011] (3) This invention determines the risk area of ​​secondary icing based on the predicted ice thickness and moisture distribution area, and generates de-icing action commands by planning energy density, action time, and motion path through surface reconstruction, ice adhesion strength calculation, and PID control, and implements position and temperature closed-loop control on the de-icing actuator. Since the risk grid connectivity aggregation can limit the de-icing range, the ice adhesion strength and predicted thickness can provide the basis for adjusting the output power and action time, and the closed-loop control can correct the position and temperature according to the execution state, it solves the problem that fixed range and fixed intensity de-icing cannot adapt to the actual ice conditions, realizes directional de-icing of the secondary icing area, and reduces the risk of equipment failure caused by secondary icing of blades. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the process of the blade secondary icing and de-icing method based on icing identification provided in an embodiment of the present invention. Detailed Implementation

[0013] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0014] Reference Figure 1 This invention provides a method for de-icing secondary icing of blades based on icing identification, comprising the following steps: S11, Acquire an initial surface image of the wind turbine blades, and perform liquid water pixel spatial clustering on the initial surface image to obtain the water distribution area; S12, collect real-time temperature data of the blade surface and real-time wind speed data of the nacelle anemometer, and combine them with the water distribution area to perform time-series deduction of ice and water state to obtain the dynamic evolution characteristics of the ice layer. S13. Thermodynamic flux is estimated based on the dynamic evolution characteristics of the ice layer to obtain the initial heat flux. The heat conduction gradient is obtained based on the initial heat flux. The freezing rate is predicted based on the heat conduction gradient to obtain the predicted ice layer thickness. S14, Based on the predicted ice thickness and the moisture distribution area, determine the risk area of ​​secondary icing to obtain the risk area or safe area of ​​secondary icing. S15, when the secondary icing risk area is obtained, the de-icing action parameters are planned according to the secondary icing risk area and the predicted ice thickness to obtain the de-icing action command, and the position and temperature closed-loop control of the de-icing actuator is performed according to the de-icing action command.

[0015] In step S11, an initial surface image of the wind turbine blades is acquired, and the initial surface image is subjected to spatial clustering of liquid water pixels to obtain the water distribution area, including: A weighted average method is used to perform pixel-by-pixel grayscale mapping on the red channel pixel values, green channel pixel values, and blue channel pixel values ​​of each pixel in the initial surface image to obtain the target grayscale matrix. The absolute value of the difference between each center pixel value and the average value of its neighboring pixels in the target grayscale matrix is ​​calculated, and the variance of the absolute value of the difference around each center pixel is calculated to obtain the pixel contrast matrix and the local contrast variance matrix. The horizontal and vertical partial derivatives of the pixel contrast matrix are obtained by using the Sobel operator. The square root of the sum of the squares of the horizontal and vertical partial derivatives is then taken to obtain the gradient magnitude matrix. Pixels whose gradient magnitude in the gradient magnitude matrix is ​​greater than a preset gradient magnitude threshold and whose corresponding variance in the local contrast variance matrix is ​​less than a preset contrast variance threshold are marked as liquid water pixels, thus obtaining a set of liquid water pixels. The set of liquid water pixels is clustered using a density space clustering method based on Euclidean distance to obtain the water distribution region.

[0016] Specifically, the initial surface image originates from a visible light image acquisition device installed at the wind turbine blade detection location, and the acquisition object is the outer surface of the wind turbine blade. The data content of the initial surface image consists of pixels arranged in rows and columns, with each pixel containing red, green, and blue channel pixel values. When performing pixel-by-pixel grayscale mapping for each pixel, the red channel pixel value is multiplied by 0.299, the green channel pixel value by 0.587, and the blue channel pixel value by 0.114. These three products are then added together to obtain the corresponding grayscale value for that pixel. The above three weights are determined based on the fixed weights for visible light image brightness grayscale conversion, and the sum of the three weights is 1. All grayscale values ​​are written into a matrix according to the original row and column positions of the initial surface image to obtain the target grayscale matrix.

[0017] After obtaining the target grayscale matrix, each pixel in the target grayscale matrix is ​​taken as the center pixel. The eight neighboring pixels within a 3×3 window surrounding the center pixel are read. The grayscale values ​​of the eight neighboring pixels are added together and divided by 8 to obtain the average value of the neighboring pixels. Then, the center pixel value is subtracted from the average value of the neighboring pixels, and the absolute value of the difference is taken to obtain the pixel contrast value corresponding to the center pixel. Boundary pixels are supplemented by copying the boundary grayscale values ​​to complete the 3×3 window. For example, if the first row of pixels lacks the pixel above it, the grayscale value of the corresponding position in the first row is copied to the pixel above it, ensuring that each center pixel has a complete 3×3 calculation window. All pixel contrast values ​​are arranged according to their original row and column positions to obtain the pixel contrast matrix. Next, using the absolute values ​​of the differences within a 3×3 window surrounding each center pixel as the calculation object, the absolute values ​​of the differences within the window are added together and divided by the number to obtain the average difference of the window. The average difference of the window is subtracted from the absolute value of each difference within the window and squared. The sum of the squared results is then divided by the number of absolute values ​​of the differences within the window to obtain the local contrast variance corresponding to that center pixel. The local contrast variance matrix is ​​obtained by arranging all the local contrast variances in their original row and column positions.

[0018] When performing Sobel calculations on the pixel contrast matrix, within a 3×3 window corresponding to each center pixel, the contrast values ​​of the three pixels in the left column are multiplied by -1, -2, and -1 respectively; the contrast values ​​of the three pixels in the right column are multiplied by 1, 2, and 1 respectively; and the middle column is multiplied by 0. These nine products are then summed to obtain the horizontal partial derivative of the center pixel. Similarly, the contrast values ​​of the three pixels in the upper row are multiplied by -1, -2, and -1 respectively; the contrast values ​​of the three pixels in the lower row are multiplied by 1, 2, and 1 respectively; and the middle row is multiplied by 0. These nine products are then summed to obtain the vertical partial derivative of the center pixel. Finally, the square root of the sum of the squared horizontal and vertical partial derivatives at the same pixel location is taken to obtain the gradient magnitude at that pixel location. All gradient magnitudes are then arranged in their original row and column positions to form a gradient magnitude matrix.

[0019] The gradient amplitude threshold is determined by historical blade surface image samples, which include blade surface images confirmed by maintenance records to contain liquid water and blade surface images confirmed by maintenance records to contain no liquid water. The maintenance records are obtained through blade surface calibration. During calibration, a calibration template with coordinate markers and a closed opening is fixed to the blade surface. Water is injected into the closed opening using a graduated titrator. Pixels covered inside the opening of the calibration template are recorded as pixels inside the liquid water layer, pixels covered one pixel width inward and one pixel width outward from the opening edge are recorded as pixels at the liquid water boundary, and pixels outside the template are recorded as pixels without liquid water. This provides a unified standard for determining the pixel states in the historical blade surface image samples.

[0020] The gradient magnitude matrix is ​​calculated for historical blade surface image samples in the same manner. Gradient magnitudes of pixels at the liquid water boundary in images confirming the presence of liquid water are extracted. These extracted gradient magnitudes are sorted in ascending order, and the value at the 25th percentile after sorting is used as the gradient magnitude threshold. The contrast variance threshold is determined by the local contrast variance of pixels within the liquid water region in images confirming the presence of liquid water. This local contrast variance is sorted in ascending order, and the value at the 75th percentile after sorting is used as the contrast variance threshold. The corresponding gradient magnitude and local contrast variance are read for each pixel in the target grayscale matrix. When the gradient magnitude of the same pixel is greater than the gradient magnitude threshold and the corresponding local contrast variance is less than the contrast variance threshold, the pixel is marked as a liquid water pixel. All liquid water pixels are aggregated according to their original row and column positions to obtain a liquid water pixel set. When performing density space clustering on the set of liquid water pixels, the row and column coordinates of each liquid water pixel are used as two-dimensional coordinates. The Euclidean distance between any two liquid water pixels is obtained by taking the square root of the sum of the square of the difference between their row coordinates and the square of the difference between their column coordinates.

[0021] The cluster radius is determined by the 90th percentile of the distance between adjacent liquid water pixels in historical blade surface image samples, and the minimum number of pixels is determined by the 10th percentile of the number of pixels in the minimum effective water region in historical blade surface image samples. The minimum effective water region refers to a water region whose actual coverage area is not less than 10% of the single-action area of ​​the de-icing actuator. First, the actual coverage area is converted into the number of pixels according to the calibration ratio of the image acquisition device. Then, the number of pixels in the historical minimum effective water region is sorted from smallest to largest, and the 10th percentile is taken as the minimum number of pixels. If the calculation result has a decimal, it is rounded up. Taking any liquid water pixel as the center, the number of liquid water pixels whose Euclidean distance to it is not greater than the cluster radius is counted. When this number is not less than the minimum number of pixels, the liquid water pixel is taken as the core pixel, and its neighboring liquid water pixels are grouped into the same cluster. Continue to perform neighborhood expansion on newly included core pixels until no new liquid water pixels are added, resulting in a water distribution region; repeat the above process for unclassified liquid water pixels in the liquid water pixel set to obtain the complete water distribution region.

[0022] In step S12, real-time temperature data of the blade surface and real-time wind speed data of the nacelle anemometer are collected. Combined with the moisture distribution area, a time-series deduction of the ice-water state is performed to obtain the dynamic evolution characteristics of the ice layer, including: Subtract the preset freezing point temperature from the real-time temperature data, and arrange the resulting temperature difference with the real-time wind speed data at the same acquisition time to obtain a phase change query sequence. The ice crystal growth rate sequence is obtained by querying the pre-constructed ice-water phase transition rate table according to the phase transition query sequence. The ice crystal growth rate sequence is converted into a state transition matrix sequence between liquid water state, ice-water mixed state and solid ice state using the Markov state transition method. An initial water state vector is constructed based on the liquid water pixels in the water distribution area, and the initial water state vector is multiplied sequentially by each state transition matrix in the state transition matrix sequence to obtain the dynamic evolution characteristics of the ice layer.

[0023] Specifically, the real-time temperature data comes from the blade surface temperature sensor, and the real-time wind speed data comes from the nacelle anemometer; the two are matched one-to-one according to the acquisition time. The real-time temperature data consists of the temperature values ​​of each measuring point on the blade surface at continuous acquisition times, and the real-time wind speed data consists of the wind speed values ​​output by the nacelle anemometer at the same acquisition time. A preset freezing point temperature is subtracted from the real-time temperature data at each acquisition time to obtain the temperature difference; the preset freezing point temperature is determined based on the water source and measured salinity of the water body covering the current moisture distribution area. For rainwater or residual water whose salinity measurement results are within the instrument's detection error range and can be approximately treated as freshwater, the freezing point temperature is set to 0℃. For sea fog deposits, saline meltwater, or other residual water with detected salinity, samples from the same source are collected, and the temperature at which the ice-water phase transition plateau begins to appear is determined using the cooling curve method, or the corresponding freezing point temperature is read from a pre-calibrated salinity-freezing point correspondence table. Freezing point parameters are established separately for different water sources; saline samples are not uniformly treated as 0℃. By arranging the temperature difference values ​​at the same acquisition time with the corresponding real-time wind speed data, a phase change query sequence is obtained.

[0024] The ice-water phase change rate table is constructed from historical blade surface operating data. This historical data includes historical real-time temperature data, historical real-time wind speed data, historical blade surface images, and liquid water, ice-water mixed, and solid ice states confirmed by maintenance records. The state corresponding to the historical blade surface images is determined by the initial water addition mass, the unfrozen water mass, and the attached ice mass. The initial water addition mass is obtained from the mass difference before and after water addition; the unfrozen water mass is obtained from the mass increase resulting from the absorption of water from the blade surface by pre-weighed absorbent material; and the attached ice mass is obtained by subtracting the unfrozen water mass from the initial water addition mass. A liquid water state is defined as having 0 attached ice mass; an ice-water mixed state is defined as having an attached ice mass greater than 0 and an unfrozen water mass greater than 0; and a solid ice state is defined as having 0 unfrozen water mass.

[0025] When constructing the ice-water phase change rate table, historical blade surface operation data were grouped according to water source and salinity range. The historical real-time temperature data was then subtracted from the freezing point temperature corresponding to the measured salinity of each group to obtain the historical temperature difference. Historical temperature differences were divided into intervals of 1℃, and historical real-time wind speed data were divided into intervals of 1 m / s. For historical samples falling within the same temperature difference and wind speed intervals and with the same sampling time interval, the number of samples transitioning from liquid water to an ice-water mixture state, from an ice-water mixture state to a solid ice state, from an ice-water mixture state back to a liquid water state, and from a solid ice state back to an ice-water mixture state were counted. Simultaneously, based on mass reduction records or outflow records from water distribution area boundaries at adjacent sampling times, the number of samples experiencing water outflow or evaporation in each current state was counted. Divide the number of each transition sample by the number of samples belonging to the corresponding current state at the current sampling time to obtain the state transition probability and mass loss probability within adjacent sampling time intervals. The probabilities are dimensionless proportions. Then divide each state transition probability by the corresponding sampling time interval to obtain the phase transition rate in s⁻¹, and write it into the ice-water phase transition rate table.

[0026] When using the Markov state transition method, liquid water, ice-water mixture, and solid ice are considered as three phase states, with a 3x3 phase transition matrix corresponding to each acquisition time. A separate water loss ratio vector, corresponding one-to-one with each of the three current phase states, is also set to record the mass loss caused by water film flowing out of the leaf detection area or water evaporation. Each row in the matrix represents the current phase state, and each column represents the phase state at the next time step, with the row and column order being liquid water, ice-water mixture, and solid ice states, respectively. The transition probabilities of liquid water transferring to an ice-water mixture, ice-water mixture transferring to solid ice, ice-water mixture returning to liquid water, and solid ice returning to an ice-water mixture are filled into the corresponding matrix elements. For the same current phase, the sum of all transition probabilities to other phases and the corresponding water loss probability is first calculated. When the sum is not greater than 1, the retention ratio of each phase is 1 minus the sum. When the sum is greater than 1, the all transition probabilities and water loss probabilities of the phase are normalized by the same ratio so that the normalized sum is 1, and then the retention ratio of the phase is set to 0. Matrix elements that have not undergone direct transition are set to 0. The sum of the phase transition matrix element in each row and the corresponding water loss probability in that row is 1, so the solid ice retention ratio is not fixed at 1. When the leaf is heated and the ambient temperature rises, solid ice can return to an ice-water mixture state, ice-water mixture state can return to liquid water state, and the total mass of each phase can decrease when the water film flows out or evaporates. An initial water state vector is constructed based on the liquid water pixels in the water distribution area, with the initial proportions of liquid water, ice-water mixture, and solid ice being 1, 0, and 0, respectively. Using each water distribution area or its spatial grid node as an independent recursive unit, a state transition matrix is ​​generated based on the real-time temperature at the corresponding location and the real-time wind speed at the same moment. The initial water state vector of each recursive unit is multiplied sequentially by the corresponding state transition matrix, and the water loss ratio is simultaneously subtracted to obtain the three-phase remaining mass ratio and cumulative loss ratio with spatial coordinates and time identifiers at each acquisition time. The remaining freezeable water volume is then calculated in conjunction with the initial water volume of each spatial grid node. After renormalizing the three-phase ratios in the remaining water, they are arranged in chronological order to form the dynamic evolution characteristics of the ice layer, while the cumulative loss ratio is retained for mass correction in ice layer thickness prediction.

[0027] In step S13, thermodynamic flux is estimated based on the dynamic evolution characteristics of the ice layer to obtain the initial heat flux. Then, the heat conduction gradient is calculated based on the initial heat flux to obtain the heat conduction gradient. Finally, the freezing rate is integrally predicted based on the heat conduction gradient to obtain the predicted ice layer thickness, including: The proportions of liquid water and solid ice states are extracted from the dynamic evolution characteristics of the ice layer according to the predicted time to obtain the liquid-solid phase transition ratio sequence. The liquid-solid phase change ratio sequence and the preset water salinity are arranged according to the corresponding prediction time to obtain the thermodynamic state input vector; The thermodynamic state input vector is input into a pre-constructed thermodynamic flux estimation model to perform heat flux mapping, thereby obtaining the initial heat flux.

[0028] Determine the blade surface grid nodes corresponding to the initial heat flux and the outward normal direction of the blade surface; According to Fourier's law of heat conduction and the outward normal direction of the blade surface, the negative value of the initial heat flux is divided by the preset equivalent thermal conductivity of the blade to obtain the normal temperature gradient of each grid node at each prediction time. The normal temperature gradient is discretized according to the grid nodes and the prediction time to obtain the heat conduction gradient.

[0029] Extract the gradient values ​​of the heat conduction gradient at each prediction time within a preset prediction period to obtain a heat conduction gradient sequence; Linear interpolation is performed on each gradient value in the heat conduction gradient sequence in the pre-constructed icing kinetic curve to obtain the icing rate sequence; The predicted ice thickness is obtained by numerically integrating the icing rate sequence over the predicted time period using the trapezoidal integral method.

[0030] Specifically, the proportions of liquid water and solid ice states are extracted from the dynamic evolution characteristics of the ice layer according to the predicted time, resulting in a liquid-solid phase transition ratio sequence. Preset water salinity is determined separately for rainwater, sea fog sediment, meltwater, or other actual water sources; water samples from different sources are not mixed and averaged. Within the same ice-covering season, at least 10 water samples are collected from the same source covering the current water distribution area, and the conductivity of each sample is measured to obtain the water sample conductivity. A water conductivity-salinity conversion table is pre-established using standard saline samples. Standard saline samples with salinities of 0, 0.5, 1, 2, 3, and 5 are prepared, and the conductivity of each standard saline sample is measured. The conductivity and corresponding salinity of each standard saline sample are then entered into the water conductivity-salinity conversion table. During the conversion, the conductivity of the water sample and the conductivity of the water quality are compared with the conductivity in the salinity conversion table. When the conductivity of the water sample equals a certain conductivity in the table, the corresponding salinity is directly read as the salinity value. When the conductivity of the water sample is between two adjacent conductivity values ​​in the table, the difference between the water sample conductivity and the smaller conductivity is calculated first, and then the difference between the larger conductivity and the smaller conductivity is calculated. The position ratio is obtained by dividing the first difference by the second difference. Then, the smaller salinity is added to the position ratio and multiplied by the difference between the two salinities to obtain the salinity value. The above conversion is performed on no less than 10 water samples, and the arithmetic mean of all salinity values ​​from the same water source is taken as the preset water salinity of that source. This salinity is also used in step S12 to determine the corresponding freezing point temperature and in step S13 to construct the thermodynamic state input vector. The liquid water state ratio and solid ice state ratio at each prediction time in the liquid-solid phase change ratio sequence are arranged in correspondence with the preset water salinity to obtain the thermodynamic state input vector.

[0031] The thermodynamic flux estimation model is a linear mapping model. The input variables are the proportion of liquid water state, the proportion of solid ice state, and a preset water salinity. The output variable is the initial heat flux. When constructing the thermodynamic flux estimation model, historical samples are grouped according to the blade surface temperature difference interval and the real-time wind speed interval, and the corresponding flux coefficients are solved for each group of historical samples. During online estimation, the flux coefficients for the corresponding intervals are retrieved based on the current temperature difference and real-time wind speed. The thermodynamic state input vector is then substituted into the flux coefficients to obtain the initial heat flux. The model's internal parameters include the liquid water state proportion flux coefficient, the solid ice state proportion flux coefficient, the salinity flux coefficient, and a constant term. When constructing the thermodynamic flux estimation model, each historical sample is written as a single line of data, which sequentially includes the liquid water state proportion, the solid ice state proportion, the salinity value, and a constant 1. These lines of data form an input data table; the measured heat fluxes corresponding to the same historical sample form an output data column.

[0032] The measured heat flux is the heat transferred per unit time through a unit blade surface area, i.e., heat flux density, with units of W / m². Taking the direction from the blade surface towards the external environment as the positive direction of heat transfer, the sensible heat flux generated by temperature change during the k-th sampling time interval is calculated using the following formula:

[0033] In the formula, is the sensible heat flux generated by temperature change during the kth sampling time interval, in W / m². The equivalent heat capacity of the blade surface, residual moisture, and attached ice layer participating in the temperature change within the heat exchange detection area is expressed in J / K. and These are the blade surface temperatures at two adjacent sampling times, in K or °C, and the difference between them is in K. The blade surface area of ​​the heat exchange detection zone is expressed in m². The time interval between two adjacent sampling moments is expressed in seconds. The equivalent heat capacity is determined by the blade material parameters and the mass and specific heat capacity of each component within the heat transfer detection area, or by a temperature rise calibration test under known heating power. The parameter calculation and calibration process further includes, according to… Calculations were performed, where mb, mw, and mi represent the masses of the blade surface layer, residual moisture, and attached ice layer involved in temperature changes within the heat transfer detection area, respectively, in kg; cb, cw, and ci represent the specific heat capacities of the corresponding materials, in J / (kg·K). During the calibration test, the blade detection area was placed under stable wind speed and ambient temperature conditions. A constant power Pcal was applied by a calibrated heater for a calibration time Δtcal, and the average surface temperature change ΔTcal before and after calibration was recorded. Simultaneously, the ambient heat dissipation Qloss during the same period was obtained from the no-load test, and then... Calculate the equivalent heat capacity. Repeat the calibration at least three times with different heating powers, take the average of the calculated results, and use the deviation of each calculated value from the average value not exceeding the preset calibration error as a validity condition. hour, This indicates that heat is transferred from the blade surface to the external environment; when hour, This indicates the heat absorbed by the blade surface. When there is a change in the ice-water phase ratio between adjacent sampling times, the latent heat of phase change is calculated based on the change in the mass of solid ice and the latent heat of phase change of water. This latent heat of phase change is then divided by the area of ​​the heat transfer detection region and the sampling time interval to obtain the latent heat flux of phase change, where the freezing process corresponds to a positive value and the melting process corresponds to a negative value. The latent heat flux of phase change is added to the sensible heat flux generated by the temperature change to obtain the measured heat flux of the corresponding historical sample. The measured heat fluxes corresponding to each historical sample are then combined to form an output data column.

[0034] When solving for coefficients in the input data table, the input data table is first transposed to obtain a transposed input data table. Then, the transposed input data table is multiplied by the input data table to obtain a 4x4 matrix. The transposed input data table is then multiplied by the output data column to obtain a 4x1 matrix. The inverse of the first matrix is ​​calculated and multiplied by the second matrix to obtain a 4x1 coefficient column. The four values ​​in the coefficient column are used as the liquid water state proportional flux coefficient, the solid ice state proportional flux coefficient, the salinity flux coefficient, and a constant term, respectively. When the first matrix does not have an inverse matrix, 0.001 is added to each diagonal element of the first matrix before calculating the inverse matrix. After construction, the liquid water state proportionality in the thermodynamic state input vector is multiplied by the liquid water state proportional flux coefficient, the solid ice state proportionality is multiplied by the solid ice state proportional flux coefficient, and the preset water salinity is multiplied by the salinity flux coefficient. These are then added to the constant term to obtain the initial heat flux.

[0035] The preset equivalent thermal conductivity of the blade is determined by the factory testing data of the blade material, or by equivalent conversion based on the thermal conductivity and thickness of the blade surface composite material, residual moisture, and attached ice layer. Taking the outward normal direction of the blade surface as the positive direction, according to Fourier's law of heat conduction and the outward normal direction of the blade surface, the normal temperature gradient of the i-th spatial grid node at the k-th prediction time is calculated according to the following formula:

[0036] In the formula, The normal temperature gradient of the i-th spatial grid node at the k-th prediction time is expressed in K / m. This represents the initial heat flux at the corresponding spatial grid node and prediction time, in W / m³. 2 ; The equivalent thermal conductivity of the blade surface structure corresponding to the spatial grid node is expressed in W / (m·K); n is the normal coordinate pointing from the blade surface to the external environment. A negative sign indicates that the direction of heat transfer is opposite to the direction of temperature increase. The normal temperature gradients of each spatial grid node at each prediction time are arranged according to the spatial grid node number and the prediction time to obtain a discrete heat conduction gradient sequence. When querying the freezing rate, the sign of the normal temperature gradient is retained to distinguish the direction of heat transfer; when heat is transferred from the blade surface to the external environment and the temperature of the corresponding spatial grid node is not higher than the freezing point temperature of the water source, the magnitude of the normal temperature gradient is used to query the freezing kinetic curve; when heat is transferred from the external environment to the blade surface, the freezing rate is determined to be 0, or the corresponding melting rate is queried when a melting kinetic curve is pre-constructed and used as the negative growth amount of the ice layer thickness.

[0037] The preset prediction period is determined by the de-icing control cycle, taking the time length between two adjacent de-icing control judgments; the prediction time interval is taken as the sampling interval of real-time temperature data. The icing kinetic curve is constructed from historical blade surface operation data, with the horizontal data being the heat conduction gradient and the vertical data being the icing rate. During construction, the heat conduction gradients in the historical samples are sorted from smallest to largest, and the increment of the measured ice layer thickness at adjacent sampling times is divided by the sampling time interval to obtain the corresponding icing rate; when multiple historical samples have the same heat conduction gradient, the corresponding icing rates are added together and divided by the number of samples to obtain the average icing rate corresponding to that heat conduction gradient. For any gradient value in the heat conduction gradient sequence, adjacent nodes less than or equal to the gradient value and adjacent nodes greater than or equal to the gradient value are found in the icing kinetic curve. Based on the positional ratio of the gradient value between the two adjacent nodes, the icing rates of the two adjacent nodes are linearly interpolated to obtain the corresponding icing rate; when the gradient value is less than the minimum heat conduction gradient of the curve, the icing rate corresponding to the minimum heat conduction gradient is taken; when the gradient value is greater than the maximum heat conduction gradient of the curve, the icing rate corresponding to the maximum heat conduction gradient is taken. The above processing is performed on the gradient values ​​of each spatial grid node at all prediction times to obtain an icing rate sequence corresponding to each spatial grid node. When using the trapezoidal integral method, the icing rates of two adjacent prediction times are added together, divided by 2, and then multiplied by the time interval between the two prediction times to obtain the ice thickness increment for that time period. For each spatial grid node, all ice thickness increments within the prediction time period are added together to obtain the integral ice thickness of that node. The remaining freezeable water volume is calculated based on the initial water volume and cumulative loss ratio of that node, and the remaining freezeable water volume is divided by the ice density and the blade surface area corresponding to that node to obtain the maximum formable ice thickness of that node. The smaller value between the integral ice thickness and the maximum formable ice thickness is taken as the predicted ice thickness of that node, and the predicted ice thickness is associated with and saved with the spatial coordinates of the corresponding grid node.

[0038] In step S14, a secondary icing risk area is determined based on the predicted ice thickness and the moisture distribution area, resulting in a secondary icing risk area or a safe area, including: When the predicted ice layer thickness is greater than a preset thickness threshold, the spatial grid nodes in the moisture distribution area are marked as risk grid nodes, thus obtaining a risk grid set; The risk grid set is connected and aggregated according to the spatial distance between adjacent risk grid nodes to obtain the secondary icing risk area; When the predicted ice layer thickness is less than or equal to the thickness threshold, the spatial grid nodes in the moisture distribution area are marked as safe grid nodes, thus obtaining a safe grid set; The safe region is obtained by fitting the convex hull boundary of the outer grid nodes in the safe grid set.

[0039] Specifically, step S13 outputs the predicted ice thickness for each spatial grid node and associates and saves the predicted ice thickness with the spatial coordinates of the corresponding grid node, thus ensuring a one-to-one correspondence between the predicted ice thickness and the water distribution area according to spatial location. When converting liquid water pixels in the water distribution area into spatial grid nodes, the intrinsic and extrinsic parameters of the image acquisition device and the coordinates of the blade surface calibration points are first read. Then, the row and column coordinates of the liquid water pixels are converted into imaging rays in the coordinate system of the image acquisition device based on the intrinsic parameters. Based on the extrinsic parameters, the imaging rays are converted to the coordinate system of the blade surface, and the intersection of the imaging ray and the blade surface calibration point is used as the spatial grid node. The thickness threshold is determined by the historical blade icing treatment record, which includes the measured ice thickness at each icing event and whether a de-icing action was triggered. The measured ice thicknesses that have triggered de-icing actions are sorted from smallest to largest, and the value at the 10th percentile after sorting is taken as the thickness threshold. For each spatial grid node in the moisture distribution area, if the predicted ice thickness corresponding to that node is greater than a preset thickness threshold, the node is marked as a risk grid node, resulting in a risk grid set. When performing connectivity aggregation on the risk grid set, the spatial distance between adjacent risk grid nodes is calculated. The spatial distance is obtained by taking the square root of the sum of the squared differences in the lateral, longitudinal, and normal coordinates of two risk grid nodes in the blade surface coordinate system. The connectivity distance threshold is determined by the 95th percentile of the distance between adjacent spatial grid nodes in the moisture distribution area. The 95th percentile is chosen to cover the vast majority of normal adjacent spacing caused by image calibration errors or local defects, while excluding the largest 5% of abnormal cross-region spacing, thus avoiding the incorrect merging of separate icing areas. When determining the percentile, historically manually labeled secondary icing areas are used as the verification benchmark. Connected regions are generated using the 90th, 92.5th, 95th, and 97.5th percentiles, respectively. The intersection-union ratio (IUR) and false merging rate of the generated regions and manually labeled regions are calculated. In this embodiment, the 95th percentile has the highest average IUR and the false merging rate meets the preset upper limit, so it is selected as the connected distance threshold. When the spatial distance between two risk grid nodes is not greater than the connected distance threshold, they are assigned to the same connected region. Adjacency judgment is performed on newly assigned risk grid nodes until no new risk grid nodes are added, thus obtaining the secondary icing risk region. For each spatial grid node in the moisture distribution region, when the predicted ice thickness corresponding to the node is less than or equal to the thickness threshold, the node is marked as a safe grid node, thus obtaining a safe grid set. When fitting the convex hull boundary of the outer grid nodes in the safe grid set, the safe grid node with the smallest lateral coordinate is selected as the starting node. Boundary nodes are selected in order of increasing polar angle. When a new node is added and the boundary bends inward, the previous boundary node is deleted until the smallest convex boundary that can surround all safe grid nodes is formed. The area enclosed by the convex boundary is taken as the safe area.

[0040] In step S15, when the secondary icing risk area is obtained, de-icing action parameters are planned based on the secondary icing risk area and the predicted ice thickness to obtain de-icing action commands. Then, the de-icing actuator is subjected to pose and temperature closed-loop control based on the de-icing action commands, including: The spatial coordinate data of the secondary icing risk area are reconstructed using Delaunay surfaces to obtain a three-dimensional spatial mesh. Based on the three-dimensional spatial mesh, the corresponding dynamic evolution features of the ice layer are extracted to obtain the ice crystal structure distribution features. The ice crystal structure distribution characteristics and the predicted ice layer thickness are input into a pre-constructed linear regression model to obtain the ice layer adhesion strength. The ice adhesion strength is input into a PID controller for control quantity mapping to obtain energy density and action time; Based on the energy density and the action time, the adjacent risk grids in the three-dimensional spatial grid are continuously arranged to obtain the de-icing action command containing the movement trajectory, movement speed and output power.

[0041] Specifically, after obtaining the secondary icing risk area, the spatial coordinate data of each risk grid node in the secondary icing risk area is read. The spatial coordinate data includes lateral coordinates, longitudinal coordinates, and normal coordinates. When reconstructing the spatial coordinate data of the secondary icing risk area using the Delaunay surface, the risk grid nodes are first projected onto the blade surface parameter plane. Adjacent nodes are connected in the parameter plane to form triangular elements, and it is checked whether there are other risk grid nodes inside the circumcircle of each triangular element. If there are other risk grid nodes inside the circumcircle, the current connecting edge of the triangular element is deleted, and the connecting edge between the relative nodes is used to re-form the triangular element until no other risk grid nodes are contained inside the circumcircle of any triangular element. Then, the triangular elements are restored to the blade surface coordinate system according to the three-dimensional spatial coordinates of the corresponding nodes to obtain a three-dimensional spatial grid. Based on the position of each grid node in the three-dimensional spatial grid, the proportion of liquid water state, ice-water mixed state, and solid ice state corresponding to the same spatial position is extracted from the dynamic evolution characteristics of the ice layer to obtain the ice crystal structure distribution characteristics.

[0042] The linear regression model is used to map the distribution characteristics of ice crystal structure and predicted ice thickness to ice adhesion strength. Since the sum of the proportions of liquid water, ice-water mixture, and solid ice states is 1, to avoid the three-phase proportions being completely collinear with the constant term, the liquid water state proportion is removed from the model input, retaining only the ice-water mixture proportion, solid ice state proportion, and predicted ice thickness, and a constant term is set. The liquid water state proportion can be obtained by subtracting the ice-water mixture proportion and solid ice state proportion from 1. When constructing the model, historical samples are divided into training and validation sets in chronological order. First, the mean and standard deviation of the ice-water mixture proportion, solid ice state proportion, and predicted ice thickness in the training set are calculated. Then, the corresponding mean is subtracted from each variable and divided by the corresponding standard deviation to obtain dimensionless standardized variables. Validation samples and online samples use the same mean and standard deviation fixed in the training set. The three standardized variables and the constant 1 for each training sample are written sequentially into a row, forming a four-column adhesion strength input data table, and the corresponding measured ice adhesion strength forms the output data column. When using ridge regression, only the regression coefficients of the three standardized variables are regularized using the L2 norm, without regularizing the constant term. The model is trained one by one from the candidate regularization coefficient set, and the mean squared error of the validation set is calculated. The regularization coefficient with the smallest validation error is selected, and 0.001 is not fixed. The four coefficients obtained correspond to the proportion of ice-water mixture, the proportion of solid ice, the predicted ice thickness, and the constant term, respectively. During online calculation, the three input variables are first standardized according to the training parameters, and then substituted into the regression coefficients to obtain the ice adhesion strength.

[0043] When inputting the ice adhesion strength into the PID controller, the control deviation is first obtained by subtracting the preset target adhesion strength from the ice adhesion strength. The preset target adhesion strength is determined by the blade de-icing test, taking the maximum residual adhesion strength of the completely de-iced sample. The internal parameters of the PID controller include proportional coefficient, integral coefficient, and derivative coefficient. The proportional coefficient, integral coefficient, and derivative coefficient are determined by step response tuning. Initially, the integral coefficient and derivative coefficient are set to 0, and the proportional coefficient is gradually increased until the de-icing actuator output reaches a stable response. 60% of this proportional coefficient is taken as the final proportional coefficient, the final proportional coefficient is divided by the stable response time to obtain the integral coefficient, and the final proportional coefficient is multiplied by one-eighth of the stable response time to obtain the derivative coefficient. The current control deviation is multiplied by the proportional coefficient, the result of accumulating the control deviation over sampling time is multiplied by the integral coefficient, the difference between the current control deviation and the control deviation at the previous sampling time is divided by the sampling time interval and multiplied by the derivative coefficient, and the three results are added together to obtain the control quantity. When mapping the control quantity to energy density and action time, the control quantity is first limited according to the rated output range of the de-icing actuator. The lower limit of energy density is the minimum energy density required for the de-icing actuator to generate continuous ice melting on the blade surface. This lower limit is determined by the energy density at which continuous ice melting first occurs during the de-icing test. The upper limit of energy density is the energy density corresponding to the allowable surface temperature rise of the blade material. This upper limit is determined by the temperature resistance test data of the blade material. The lower limit of action time is the shortest response time required for the de-icing actuator to complete one stable output. The upper limit of action time is the maximum continuous output time before the blade surface temperature reaches the allowable surface temperature. The normalized control quantity is obtained by subtracting the lower limit of the control quantity from the limited control quantity and then dividing by the difference between the upper and lower limits. The energy density is obtained by subtracting the lower limit of energy density from the upper limit of energy density, multiplying by the normalized control quantity, and then adding the lower limit of energy density. The action time is obtained by subtracting the lower limit of action time from the upper limit of action time, multiplying by the normalized control quantity, and then adding the lower limit of action time. Based on energy density and action time, adjacent risk grids in the 3D spatial grid are arranged into continuous paths according to spatial distance from near to far to obtain the movement trajectory. The movement speed is obtained by dividing the movement trajectory length by the action time, and the output power is obtained by multiplying the energy density by the action area of ​​the de-icing actuator and then dividing by the action time. Finally, a de-icing action command containing the movement trajectory, movement speed, and output power is obtained, and the de-icing actuator is subjected to pose and temperature closed-loop control according to the de-icing action command. When the de-icing action command is executed, the actual pose of the actuator and the blade surface temperature of the current action area are collected according to the feedback cycle of the de-icing actuator, and compared with the target pose and target temperature corresponding to the movement trajectory to obtain pose deviation and temperature deviation. The movement trajectory or movement speed is corrected according to the pose deviation, and the output power and action time are corrected according to the temperature deviation. When the blade surface temperature reaches the target temperature or the allowable surface temperature of the blade material, the energy output of the corresponding risk grid is stopped.

[0044] In summary, this invention discloses a blade secondary icing removal method based on icing recognition, comprising: acquiring initial surface images of wind turbine blades and performing spatial clustering of liquid water pixels to obtain the moisture distribution area; acquiring real-time temperature data of the blade surface and real-time wind speed data of the nacelle anemometer, and combining the moisture distribution area to perform time-series recursion of ice-water state to obtain the dynamic evolution characteristics of the ice layer; estimating thermodynamic flux, calculating the heat dissipation conduction gradient, and predicting the icing rate integral based on the dynamic evolution characteristics of the ice layer to obtain the predicted ice layer thickness; determining the secondary icing risk area based on the predicted ice layer thickness and the moisture distribution area to obtain the secondary icing risk area or safe area; when the secondary icing risk area is obtained, planning the de-icing action parameters based on the secondary icing risk area and the predicted ice layer thickness to obtain the de-icing action command, and performing posture and temperature closed-loop control of the de-icing actuator based on the de-icing action command.

[0045] This invention identifies and spatially clusters liquid water pixels in images of wind turbine blade surfaces, using the water distribution area as the spatial basis for predicting secondary icing. It obtains the dynamic evolution characteristics of the ice layer through temporal recursion between real-time temperature data, real-time wind speed data, and the water distribution area. The predicted ice thickness is obtained through continuous calculation of heat flux, heat conduction gradient, and icing rate. Furthermore, by determining the risk area for secondary icing and planning the de-icing action parameters, the icing identification results are transformed into the movement trajectory, speed, and output power of the de-icing actuator, achieving early identification of secondary icing areas on the blades, predictive assessment of ice thickness, and closed-loop control of the de-icing action.

[0046] The second embodiment of the present invention provides a blade secondary icing and de-icing system based on icing identification, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the method described above.

[0047] It should be noted that the blade secondary icing and de-icing system based on icing identification provided in this embodiment of the invention is used to execute all the process steps of the blade secondary icing and de-icing method based on icing identification in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0048] It should be noted that the system embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0049] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A method for de-icing secondary icing of blades based on icing recognition, characterized in that, include: An initial surface image of the wind turbine blades is acquired, and the liquid water pixel spatial clustering is performed on the initial surface image to obtain the water distribution area; Real-time temperature data of the blade surface and real-time wind speed data of the nacelle anemometer are collected, and the ice-water state time sequence is recursively extrapolated by combining the moisture distribution area to obtain the dynamic evolution characteristics of the ice layer. Thermodynamic flux is estimated based on the dynamic evolution characteristics of the ice layer to obtain the initial heat flux. The heat conduction gradient is obtained based on the initial heat flux. The ice formation rate is then predicted based on the integral of the heat conduction gradient to obtain the predicted ice layer thickness. Based on the predicted ice thickness and the moisture distribution area, the risk area of ​​secondary icing is determined, and the risk area or safe area of ​​secondary icing is obtained. When the secondary icing risk area is obtained, the de-icing action parameters are planned according to the secondary icing risk area and the predicted ice thickness to obtain the de-icing action command, and the position and temperature closed-loop control of the de-icing actuator is performed according to the de-icing action command.

2. The method for secondary icing and de-icing of blades based on icing identification according to claim 1, characterized in that, The step of performing liquid water pixel spatial clustering on the initial surface image to obtain the water distribution region includes: A weighted average method is used to perform pixel-by-pixel grayscale mapping on the red channel pixel values, green channel pixel values, and blue channel pixel values ​​of each pixel in the initial surface image to obtain the target grayscale matrix. The absolute value of the difference between each center pixel value and the average value of its neighboring pixels in the target grayscale matrix is ​​calculated, and the variance of the absolute value of the difference around each center pixel is calculated to obtain the pixel contrast matrix and the local contrast variance matrix. The horizontal and vertical partial derivatives of the pixel contrast matrix are obtained by using the Sobel operator. The square root of the sum of the squares of the horizontal and vertical partial derivatives is then taken to obtain the gradient magnitude matrix. Pixels whose gradient magnitude in the gradient magnitude matrix is ​​greater than a preset gradient magnitude threshold and whose corresponding variance in the local contrast variance matrix is ​​less than a preset contrast variance threshold are marked as liquid water pixels, thus obtaining a set of liquid water pixels. The set of liquid water pixels is clustered using a density space clustering method based on Euclidean distance to obtain the water distribution region.

3. The method for secondary icing and de-icing of blades based on icing identification according to claim 1, characterized in that, The process of performing time-series recursion on the ice-water state to obtain the dynamic evolution characteristics of the ice layer includes: Subtract the preset freezing point temperature from the real-time temperature data, and arrange the resulting temperature difference with the real-time wind speed data at the same acquisition time to obtain a phase change query sequence. The ice crystal growth rate sequence is obtained by querying the pre-constructed ice-water phase transition rate table according to the phase transition query sequence. The ice crystal growth rate sequence is converted into a state transition matrix sequence between liquid water state, ice-water mixed state and solid ice state using the Markov state transition method. An initial water state vector is constructed based on the liquid water pixels in the water distribution area, and the initial water state vector is multiplied sequentially by each state transition matrix in the state transition matrix sequence to obtain the dynamic evolution characteristics of the ice layer.

4. The method for secondary icing and de-icing of blades based on icing identification according to claim 1, characterized in that, The step of estimating the thermodynamic flux based on the dynamic evolution characteristics of the ice layer to obtain the initial heat flux includes: The proportions of liquid water and solid ice states are extracted from the dynamic evolution characteristics of the ice layer according to the predicted time to obtain the liquid-solid phase transition ratio sequence. The liquid-solid phase change ratio sequence and the preset water salinity are arranged according to the corresponding prediction time to obtain the thermodynamic state input vector; The thermodynamic state input vector is input into a pre-constructed thermodynamic flux estimation model to perform heat flux mapping, thereby obtaining the initial heat flux.

5. The method for secondary icing and de-icing of blades based on icing identification according to claim 1, characterized in that, The step of calculating the heat conduction gradient based on the initial heat flux to obtain the heat conduction gradient includes: Determine the blade surface grid nodes corresponding to the initial heat flux and the outward normal direction of the blade surface; According to Fourier's law of heat conduction and the outward normal direction of the blade surface, the negative value of the initial heat flux is divided by the preset equivalent thermal conductivity of the blade to obtain the normal temperature gradient of each grid node at each prediction time. The normal temperature gradient is discretized according to the grid nodes and the prediction time to obtain the heat conduction gradient.

6. The method for secondary icing and de-icing of blades based on icing identification according to claim 1, characterized in that, The step of integrally predicting the freezing rate based on the heat conduction gradient to obtain the predicted ice thickness includes: Extract the gradient values ​​of the heat conduction gradient at each prediction time within a preset prediction period to obtain a heat conduction gradient sequence; Linear interpolation is performed on each gradient value in the heat conduction gradient sequence in the pre-constructed icing kinetic curve to obtain the icing rate sequence; The predicted ice thickness is obtained by numerically integrating the icing rate sequence over the predicted time period using the trapezoidal integral method.

7. The method for secondary icing and de-icing of blades based on icing identification according to claim 1, characterized in that, The step of determining the secondary icing risk area based on the predicted ice thickness and the moisture distribution area to obtain the secondary icing risk area or safe area includes: When the predicted ice layer thickness is greater than a preset thickness threshold, the spatial grid nodes in the moisture distribution area are marked as risk grid nodes, thus obtaining a risk grid set; The risk grid set is connected and aggregated according to the spatial distance between adjacent risk grid nodes to obtain the secondary icing risk area; When the predicted ice thickness is less than or equal to the thickness threshold, the spatial grid nodes in the moisture distribution area are marked as safe grid nodes, thus obtaining a safe grid set; The safe region is obtained by fitting the convex hull boundary of the outer grid nodes in the safe grid set.

8. The method for secondary icing and de-icing of blades based on icing identification according to claim 1, characterized in that, The step of planning de-icing action parameters based on the secondary icing risk area and the predicted ice thickness to obtain de-icing action instructions includes: The spatial coordinate data of the secondary icing risk area are reconstructed using Delaunay surfaces to obtain a three-dimensional spatial mesh. Based on the three-dimensional spatial mesh, the corresponding dynamic evolution features of the ice layer are extracted to obtain the ice crystal structure distribution features. The ice crystal structure distribution characteristics and the predicted ice layer thickness are input into a pre-constructed linear regression model to obtain the ice layer adhesion strength. The ice adhesion strength is input into a PID controller for control quantity mapping to obtain energy density and action time; Based on the energy density and the action time, the adjacent risk grids in the three-dimensional spatial grid are continuously arranged to obtain the de-icing action command containing the movement trajectory, movement speed and output power.

9. A blade secondary icing and de-icing system based on icing recognition, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the method according to any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the program implements the method described in any one of claims 1 to 8.

Images