A wind turbine blade unmanned aerial vehicle inspection method and system based on multi-source data

Abstract

The application relates to the technical field of the wind energy industry, and discloses a wind turbine blade unmanned aerial vehicle inspection method and system based on multi-source data, which comprises target radiation images and target echo data of wind turbine blades; initial mapping relationships are obtained by extracting structural features and performing spatial mapping processing according to the target radiation images and the target echo data; deviation correction processing is performed on the initial mapping relationships, and alignment fusion is performed in combination with the structural features to obtain alignment fusion data; the evolution process of damage deterioration is deduced according to the alignment fusion data, the overall severity of wind turbine blade damage evolution is evaluated, and an initial damage degree score is obtained; error compensation correction is performed on the initial damage degree score, a modified damage degree score is obtained, the blade state is evaluated according to the modified damage score, and a blade state evaluation result is obtained. The method can solve the problem that a wind turbine blade state evaluation result is unreliable.

Description

Technical Field

[0001] This invention relates to the field of wind energy technology, and in particular to a method and system for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data. Background Technology

[0002] Currently, wind turbine blades are affected by fatigue loads, environmental erosion, and material aging during long-term operation, resulting in the gradual accumulation of internal cracks, delamination, and surface damage.

[0003] In one existing technology, current wind turbine blade drone inspection technology typically uses a high-definition camera mounted on a drone to capture images of the blade surface. Image processing algorithms are then used to calculate the location coordinates of surface defects, and LiDAR scans the blade profile to infer the location coordinates of internal structural defects. These two methods are then fused according to a preset spatial coordinate mapping rule to obtain an overall picture of the defects, and the condition of the wind turbine blade is then assessed based on these defects. However, in practical applications in the wind energy industry, wind turbine blades are flexible rotating components. This approach can only handle the real-time positioning of a static, single blade, and the fusion rules are designed for a fixed blade geometry model, failing to adapt to adaptive adjustments in dynamic environments. This results in unreliable wind turbine blade condition assessment results.

[0004] In summary, existing technologies suffer from unreliable results in assessing the condition of wind turbine blades. Summary of the Invention

[0005] This invention provides a method and system for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data, in order to solve the problem of unreliable wind turbine blade condition assessment results.

[0006] Firstly, in order to solve the above-mentioned technical problems, the present invention provides a method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data, comprising: Target radiation images and target echo data of wind turbine blades were obtained through drone inspections. Thermal anomaly features on the surface of the wind turbine blades are extracted based on the target radiation image, and defect structure features inside the wind turbine blades are extracted based on the target echo data. The thermal anomaly features and the defect structure features are spatially mapped to obtain an initial mapping relationship; The initial mapping relationship is subjected to deviation correction processing, and the thermal anomaly features and the defect structure features are aligned and fused according to the corrected mapping relationship to obtain aligned and fused data. The aligned and fused data is divided into two periods of data: current and historical data. The spatial morphological features of the two periods of data are extracted, and first and second damage deterioration features are generated based on the spatial morphological features. Based on the first and second damage deterioration characteristics, the evolution process of damage deterioration is deduced and the overall severity of wind turbine blade damage evolution is assessed to obtain an initial damage severity score. The initial damage score is corrected by error compensation to obtain a corrected damage score, and the leaf condition is evaluated based on the corrected damage score to obtain the leaf condition evaluation result.

[0007] Secondly, the present invention provides a wind turbine blade unmanned aerial vehicle (UAV) inspection system based on multi-source data, 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] Compared with the prior art, the present invention has the following beneficial effects: (1) The present invention first filters the initial surface radiation features and internal echo features in the inspection data to obtain the target surface radiation features and target echo data; by cleaning and denoising the original signal, the noise effect caused by environmental interference is eliminated, ensuring the integrity and purity of the data; this standardized processing lays a high-quality data foundation for the subsequent accurate extraction of surface and internal features, so that the subsequent analysis can more accurately reflect the real physical state of the wind turbine blades; (2) This invention extracts thermal anomaly features by calculating the temperature gradient of the target radiation image, and generates a scattering cross section matrix by using a phase mutation sequence and maps it to defect structure features. This parallel extraction and characterization mechanism of multi-source heterogeneous data realizes the dual capture of blade surface damage and internal defects, solves the problem that traditional single data source is difficult to fully characterize complex damage morphology, and thus significantly improves the coverage and identification accuracy of defect detection. (3) The present invention generates a projection transformation matrix based on the extracted thermal anomaly features and defect structure features and performs projection transformation. If the position deviation is greater than the preset threshold, iterative optimization is performed to determine the initial spatial position mapping relationship. On this basis, the data is transformed and spatially overlapped by matching feature point sets to finally obtain aligned and fused data. This registration mechanism of dual verification and iterative optimization effectively overcomes the cumulative spatial registration error caused by the flexible deformation of wind turbine blades and cross-period inspection, and realizes high-precision rigid fusion of multi-source features in a unified coordinate system. (4) This invention performs time segmentation on the aligned and fused data to obtain current and historical data, extracts spatial morphological features to generate damage deterioration features, and then calculates the deterioration increment and spatial expansion values, constructs an evolution trajectory model and determines the quantitative value of the evolution process; this process breaks through the limitation of existing technologies that can only process static single-period data, and uses quantitative indicators to characterize the dynamic evolution trend of damage over time, providing intuitive data support for accurately assessing the long-term degradation process of blade structure; (5) The present invention determines whether an error compensation mechanism needs to be activated based on the quantified value of the evolution process, corrects the fused data and recalculates the quantified value through the deviation distribution matrix, and then predicts the remaining life based on the modal frequency offset and generates a maintenance priority sequence. This feedback closed-loop evaluation system ensures the accuracy and timeliness of the blade condition evaluation results, and can dynamically adjust the evaluation conclusions according to the real-time structural stress changes, thereby effectively improving the scientificity and reliability of wind turbine blade health management. Attached Figure Description

[0009] Figure 1 This is a schematic diagram of a method for inspecting wind turbine blades using unmanned aerial vehicles (UAVs) based on multi-source data, provided in the first embodiment of the present invention. Detailed Implementation

[0010] 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.

[0011] Reference Figure 1 The first embodiment of the present invention provides a method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data, including the following steps: S11, obtains target radiation images and target echo data of wind turbine blades through drone inspection; S12, extract the thermal anomaly features of the wind turbine blade surface based on the target radiation image, and extract the defect structure features inside the wind turbine blade based on the target echo data; S13, Perform spatial mapping processing on the thermal anomaly features and the defect structure features to obtain an initial mapping relationship; S14, perform deviation correction processing on the initial mapping relationship, and align and fuse the thermal anomaly features and the defect structure features according to the corrected mapping relationship to obtain aligned and fused data; S15, perform time segmentation on the aligned and fused data to obtain two periods of data: current and historical data. Extract the spatial morphological features of the two periods of data and generate first and second damage deterioration features based on the spatial morphological features. S16. Based on the first and second damage deterioration characteristics, the evolution process of damage deterioration is deduced and the overall severity of the wind turbine blade damage evolution is assessed to obtain an initial damage degree score. S17, perform error compensation correction on the initial damage degree score to obtain a corrected damage degree score, and evaluate the leaf condition based on the corrected damage score to obtain the leaf condition evaluation result.

[0012] In step S11, it is necessary to obtain target radiation images and target echo data of the wind turbine blades through drone inspection, including: The target's radiation image is acquired using an infrared thermal imager mounted on a drone, and the target's echo data is acquired using a microwave radar.

[0013] It should be noted that when acquiring target radiation images using an infrared thermal imager mounted on a UAV and target echo data using a microwave radar, the system obtains initial surface radiation images and initial internal echo data acquired by the UAV at different times. The initial surface radiation image is subjected to mean filtering to remove isolated noise points, resulting in the target surface radiation image, which is a two-dimensional digital grayscale matrix, with the matrix element values ​​representing the relative radiance of the corresponding pixel. The initial internal echo data is subjected to low-pass filtering to suppress high-frequency clutter, resulting in the target internal echo data, which is a one-dimensional temporal complex sequence, with the sequence element values ​​representing the echo amplitude and phase at the corresponding sampling time.

[0014] In step S12, it is necessary to extract the thermal anomaly features of the wind turbine blade surface based on the target radiation image, and extract the defect structure features inside the wind turbine blade based on the target echo data, including: Calculate the temperature gradient of the target radiation image, extract hot spot regions from the target radiation image based on the temperature gradient and perform morphological analysis to obtain thermal anomaly features; A phase abrupt change sequence is extracted from the target echo data. If the phase abrupt change sequence meets a preset abrupt change condition, the phase abrupt change sequence is inversely synthesized to generate a scattering cross section matrix. The defect structure characteristics are obtained by spatial mapping of the scattering cross-section matrix.

[0015] It should be noted that when calculating the temperature gradient of the target radiation image, the system calculates the horizontal and vertical grayscale change rates for each pixel. A 2×2 symmetric matrix is ​​constructed using these calculated grayscale change rates. The values ​​in the symmetric matrix, from top to bottom and left to right, represent the square of the horizontal grayscale change rate, the product of the horizontal and vertical grayscale change rates, the square of the vertical grayscale change rate, and the product of the horizontal and vertical grayscale change rates, respectively. Two eigenvalues ​​are then obtained by solving the eigenvalue formula of this matrix. The larger eigenvalue is taken as the temperature gradient value for that point. After traversing the entire image, the temperature gradient matrix is ​​output.

[0016] Next, the system compares the temperature gradients in the temperature gradient matrix with a preset anomaly threshold pixel by pixel. Using a 3×3 pixel window as a local area, if the gradient values ​​of all nine pixels within this window are greater than the preset anomaly threshold, the region is considered to have continuously exceeded the threshold. The pixel with the largest gradient in this region is used as a seed point. Each neighboring pixel in the four-neighborhood of the seed point is checked, and the difference between the gradient value of the neighboring pixel and the gradient value of the seed point is calculated. If the difference is less than a preset gradient similarity tolerance, the neighboring pixel is included in the current region, and the newly included pixel is used as a new seed point to continue expanding outwards. This process is repeated until no new neighboring pixels meet the condition. Finally, the connected pixel set is extracted as the hot spot region. When all temperature gradients in the temperature gradient matrix do not exceed the preset anomaly threshold, the system determines that the surface temperature distribution in this region is uniform and there are no signs of abnormal thermal stress concentration. Therefore, the hot spot region extraction and morphological analysis process is not triggered, and the current region is directly marked as having a normal surface condition, skipping the subsequent thermal anomaly feature generation steps. This reduces unnecessary computational overhead and avoids misjudging normal textures as damage.

[0017] The preset gradient similarity tolerance is determined by collecting manually labeled real hot spot areas from multiple sets of historical inspection data, statistically analyzing the distribution of gradient differences between adjacent pixels within each real hot spot, and taking the 85th percentile of all differences as the preset gradient similarity tolerance.

[0018] It is worth noting that the anomaly threshold is based on infrared calibration experimental data from multiple sets of known defective and normal leaves. The distribution of local maximum values ​​in the temperature gradient matrix of the two types of samples is statistically analyzed and sorted from smallest to largest. The 75th percentile value of the local maximum value in the normal leaf sample and the 25th percentile value of the local maximum value in the defective leaf sample are taken, and the average of these two values ​​is used as the initial threshold. Finally, the threshold is fixed in the system.

[0019] Subsequently, a circular structuring element with a radius of 2 pixels is defined as a probe for morphological operations. An opening operation is performed on the binary image of the hotspot region. First, an erosion operation is performed, sliding the structuring element across the image. The center pixel is retained only when the structuring element is completely covered by foreground pixels, thus removing isolated noise pixels and small connecting bridges. Then, a dilation operation is performed, sliding the structuring element across the eroded image. As long as there are foreground pixels within the range of the structuring element, the center pixel is set as foreground, thus restoring the main outline of the hotspot to make it close to the original shape. Then, a closing operation is performed on the opening operation result. First, a dilation operation is performed to fill the small voids inside the hotspot caused by noise or low contrast. Then, an erosion operation is performed to restore the edge expansion caused by dilation, making the hotspot boundary smooth and continuous. After all operations are completed, the aspect ratio and connected area of ​​each connected region are calculated. Regions with excessively large aspect ratios or excessively small areas are identified as pseudo-hotspots and are removed. Finally, the thermal anomaly features reflecting the true damage morphology are output.

[0020] The criteria for defining excessively large aspect ratios or excessively small areas are as follows: For each connected region, the ratio of the longer side to the shorter side of its smallest bounding rectangle is calculated as the aspect ratio. If this ratio is greater than a preset ratio threshold of 3, it is considered that the aspect ratio is too large. The total number of pixels contained in the connected region is calculated as the connected area. If this area is less than a preset area threshold of 50 pixels, it is considered that the area is too small. The preset ratio threshold and preset area threshold are set by collecting multiple sets of historical inspection data, manually confirmed real hot spot areas and pseudo hot spot areas, and statistically analyzing their aspect ratio and area distribution ranges. The 95th percentile of the aspect ratio of real hot spots is taken as the upper limit of the threshold, and the minimum area of ​​real hot spots is taken as the lower limit of the area. After verification and adjustment based on on-site measurement data, the thresholds are fixed in the system.

[0021] Then, when extracting the phase change sequence from the target echo data, the system extracts the echo envelope signal of the target echo data and converts it to the frequency domain through a fast Fourier transform. After setting the negative frequency component to zero, the positive frequency component is doubled, and then a complex analytic signal is obtained through an inverse Fourier transform. The imaginary part of this complex signal is the Hilbert transform result of the original signal. Subsequently, the principal argument value of the complex analytic signal is calculated for each sampling point to obtain the instantaneous phase sequence. Then, the difference operation between adjacent sampling points is performed on the instantaneous phase sequence, that is, the instantaneous phase of the next sampling point is subtracted from the instantaneous phase of the previous sampling point. If the difference exceeds the range of positive and negative π, phase unwrapping is performed by adding or subtracting 2π to eliminate the jump. Finally, a set of difference sequences reflecting the rate of phase change is output as the phase change sequence.

[0022] If the phase jump sequence meets the preset jump condition, the preset jump condition refers to whether there is a phase jump greater than a preset angle threshold in the phase jump sequence. If there is, it indicates that the microwave signal has encountered a physical interface with a large difference in dielectric constant inside the blade. Then, according to the echo signal segment corresponding to the phase jump sequence, the echo signal segment is first pulse compressed along the fast time direction. The broadband linear frequency modulated signal is compressed into a narrow pulse through matched filtering to obtain high-resolution range information. Then, motion compensation is performed on multiple pulse echoes along the slow time direction to eliminate the phase error caused by blade rotation or UAV jitter. Then, a fast Fourier transform is performed on the compensated signal along the slow time direction to convert the Doppler frequency shift into the azimuth position to obtain a range-Doppler two-dimensional image. Finally, constant false alarm rate detection and normalization processing are performed on the image to extract the scattering intensity value of each pixel and generate a scattering cross-section matrix that reflects the scattering characteristics of the target area.

[0023] The process of setting the angle threshold is as follows: microwave radar echo data of blades of the same type under healthy conditions are collected, the normal phase jump distribution caused by layup and adhesive joints is statistically analyzed, the 95th percentile of this distribution is taken as the initial value, and then the lower limit of the jump of known defect samples is combined for engineering correction, and finally the value is fixed in the system. If the phase abrupt change sequence does not meet the preset abrupt change conditions, the system determines that there is no obvious dielectric constant abrupt change interface inside the region, and therefore does not trigger inverse synthetic aperture radar imaging processing, and marks the region as having a normal internal structure, skipping the subsequent scattering cross section matrix generation step.

[0024] Finally, when performing spatial mapping on the scattering cross section matrix, the system first reads the latitude, longitude, altitude, and attitude angle data from the UAV's flight trajectory parameters, where the attitude angles include pitch, roll, and yaw angles. Then, based on these parameters, a transformation relationship from a two-dimensional image coordinate system to a three-dimensional world coordinate system is constructed. The row and column coordinates of each pixel in the scattering cross section matrix are considered as its position on the two-dimensional image plane. The spatial pointing of the pixel relative to the UAV sensor is calculated using the attitude angle data. This pointing is then extended to the three-dimensional space where the blade is located, using the UAV's altitude and latitude and longitude, to obtain the spatial coordinates of the pixel in the blade's three-dimensional physical coordinate system. The scattering intensity value of the pixel is then filled into the corresponding spatial coordinates. After completing the mapping of all pixels, a three-dimensional spatial point cloud dataset is output.

[0025] Subsequently, based on the blade design drawings, a three-dimensional geometric model of the wind turbine blade is established using three-dimensional geometric modeling software. A regular three-dimensional voxel grid is constructed within the spatial range of the blade's three-dimensional geometric model. The scattering intensity value of each point is filled into its respective voxel grid, and finally, a three-dimensional voxel dataset containing scattering intensity values ​​is output as a defect structural feature.

[0026] In step S13, spatial mapping processing is required between the thermal anomaly features and the defect structure features to obtain an initial mapping relationship, including: Extract the set of anomalous centroid coordinates corresponding to the thermal anomaly features and the sequence of structural boundary nodes corresponding to the defect structure features; The set of abnormal centroid coordinates is aligned with the sequence of structural boundary nodes in the main direction to obtain a feature registration vector. The feature registration vector is subjected to projection transformation to establish an initial mapping relationship.

[0027] It should be noted that the thermal anomaly feature is a binary mask image preserved after morphological opening and closing operations, where foreground pixels mark the connected regions of the actual damaged hotspots. The system takes this binary mask of the thermal anomaly feature as input and calculates its zeroth and first moments. The zeroth moment is the sum of the gray values ​​of all pixels within the hotspot region; the first moments are divided into horizontal and vertical first moments, calculated by summing the products of the x-coordinate and gray value of each pixel, and the y-coordinate and gray value of each pixel, respectively. Then, the horizontal first moment is divided by the zeroth moment to obtain the x-coordinate of the centroid, and the vertical first moment is divided by the zeroth moment to obtain the y-coordinate of the centroid, thereby extracting the two-dimensional pixel coordinates representing the geometric center of the anomaly region, forming a set of anomaly centroid coordinates.

[0028] Simultaneously, the system extracts the structural boundary node sequence corresponding to the defect structural features. The defect structural features refer to a complete three-dimensional voxel dataset obtained by spatial mapping of the scattering cross-section matrix. Each voxel in this dataset stores the scattering intensity value at its corresponding spatial location. Next, using this three-dimensional voxel dataset as input, the system calculates the scattering intensity difference between each voxel and its adjacent voxels along the x, y, and z axes, obtaining the first-order gradient components in these three directions. The system then takes the square root of the sum of the squares of these three components to obtain the gradient magnitude of the voxel. After traversing all voxels, a gradient magnitude voxel matrix is ​​generated. A gradient threshold is then set, and voxels with gradient magnitudes greater than this threshold are marked as candidate edge points. Non-maximum suppression is then performed on the candidate edge points along the gradient direction, retaining only voxels with locally maximum magnitudes in the gradient direction as final edge points. Finally, the three-dimensional spatial coordinates of these edge points are output, constructing the structural boundary node sequence.

[0029] The gradient threshold setting process is based on internal echo calibration experimental data of multiple sets of known defective blades and normal blades. The distribution range of voxel gradient amplitude in the two types of samples is statistically analyzed. The mean gradient amplitude of the normal sample plus three times the standard deviation is taken as the initial threshold. Then, the threshold is verified and adjusted through manually labeled defect boundary samples to ensure that the threshold can accurately distinguish between the real structural edge and noise fluctuation.

[0030] Then, all two-dimensional pixel coordinates (x, y) in the anomaly centroid coordinate set are converted to three-dimensional homogeneous coordinates (x, y, 0). The 3×3 covariance matrix of the converted anomaly centroid coordinate set is calculated, and the eigenvalues ​​and eigenvectors of this covariance matrix are solved. The eigenvector corresponding to the largest eigenvalue is taken as the principal direction vector. Similarly, the same 3×3 covariance matrix calculation and eigenvalue decomposition operation is performed on all coordinate points in the structural boundary node sequence to obtain their principal direction vectors. Subsequently, the rotation matrix between the two principal direction vectors is calculated so that one principal direction vector is parallel to the other after rotation. At the same time, the geometric center coordinates of the anomaly centroid coordinate set and the structural boundary node sequence are calculated. The two geometric center coordinates are subtracted to obtain the translation vector, which is a 3×1 column vector. The rotation matrix and the translation vector are combined into a 4×4 homogeneous transformation matrix, i.e., the projection transformation matrix. The 3×3 sub-block in the upper left corner is the rotation matrix, the 3×1 sub-block in the upper right corner is the translation vector, the 1×3 sub-block in the lower left corner is the zero vector, and the element in the lower right corner is 1.

[0031] Next, the projection transformation matrix is ​​transformed. First, each two-dimensional pixel coordinate (x, y) in the set of anomalous centroid coordinates is expanded into a four-dimensional homogeneous coordinate form, that is, the z-coordinate is supplemented with 0 and a dimension with a value of 1 is added, resulting in a four-dimensional homogeneous coordinate (x, y, 0, 1). This four-dimensional homogeneous coordinate is then multiplied by the projection transformation matrix to obtain a new four-dimensional homogeneous coordinate (x', y', z', w). Finally, the first three components of this four-dimensional homogeneous coordinate are divided by the fourth component w to normalize it, resulting in the final three-dimensional spatial coordinate (x' / w, y' / w, z' / w). The set of three-dimensional coordinates output after the above transformation of all anomalous centroid coordinates forms the first mapping matrix.

[0032] If the positional deviation corresponding to the first mapping matrix is ​​greater than the preset mapping optimization threshold, the first mapping matrix is ​​iteratively optimized to obtain the second mapping matrix. Specifically, the root mean square error of the first mapping matrix during feature point alignment is calculated as the positional deviation value. If this value is greater than the preset deviation threshold, the first mapping matrix is ​​iteratively optimized. In each iteration, each point in the abnormal centroid coordinate set is first transformed to the three-dimensional coordinate system of the structural boundary node sequence according to the current transformation parameters. Then, for each transformed point, the node with the smallest Euclidean distance in the structural boundary node sequence is searched as its corresponding point, and the association relationship of all point pairs is established. Then, based on the coordinate data of all point pairs, the optimal rotation matrix and translation vector that minimize the root mean square error of the point pairs are calculated by the singular value decomposition method. The current transformation parameters are updated with the rotation matrix and translation vector to complete one iteration. Then, the root mean square error of the point pairs is recalculated with the updated transformation parameters. If the error is still greater than the preset deviation threshold, the next iteration is entered until the root mean square error is less than the preset deviation threshold or the maximum number of iterations of 50 is reached. Finally, the transformation matrix after iterative convergence is output as the second mapping matrix. If the positional deviation does not exceed the preset deviation threshold, the first mapping matrix is ​​used directly.

[0033] The preset mapping optimization threshold and the preset deviation threshold were both determined through offline calibration experiments. Specifically, multiple sets of pairing data of surface thermal anomaly centroids and internal structural boundary nodes with known spatial correspondences were collected. An affine transformation was performed on each set of data to obtain the first mapping matrix and its root mean square error was calculated. The mean and standard deviation of all group errors were statistically analyzed. The mean plus twice the standard deviation was used as the mapping optimization threshold to determine whether to start iterative optimization. At the same time, the same batch of paired data was input into the algorithm and the root mean square error after each iteration was recorded. The maximum value of the error when all groups of data finally converged was taken as the deviation threshold to control the iteration termination condition and ensure that the optimized registration accuracy meets the engineering requirements. The two thresholds were solidified in the system after cross-validation.

[0034] Finally, when establishing the initial spatial mapping relationship between thermal anomaly features and defect structure features based on the second mapping matrix, the two-dimensional anomaly centroid coordinates are transformed into the three-dimensional structural boundary node coordinate system through the second mapping matrix, and the spatial correspondence between the temperature anomaly region on the blade surface and the internal physical defect boundary is established.

[0035] In step S14, the initial mapping relationship needs to be corrected for deviation. Based on the corrected mapping relationship, the thermal anomaly features and the defect structure features are aligned and fused to obtain aligned and fused data, including: Calculate the positional deviation of the initial mapping relationship. If the positional deviation exceeds a preset spatial deviation threshold, extract the two-dimensional radiation feature point set of the thermal anomaly feature and the three-dimensional echo feature point set of the defect structure feature. Otherwise, use the initial mapping relationship as a corrected mapping relationship. The two-dimensional radiation feature point set and the three-dimensional echo feature point set are matched by nearest neighbor search. The initial mapping relationship is corrected according to the matching result to obtain the corrected mapping relationship. The thermal anomaly features and the defect structure features are aligned and fused according to the modified mapping relationship to obtain aligned and fused data.

[0036] It should be noted that when calculating the positional deviation of the initial mapping relationship, the spatial coordinates of the corresponding feature points in the first mapping matrix and the second mapping matrix in step S13 are subtracted point by point, and the root mean square value of the Euclidean distance of all corresponding points is calculated as the positional deviation value. If the positional deviation value exceeds the preset spatial deviation threshold, a two-dimensional radiation feature point set is extracted from the thermal anomaly features, and a three-dimensional echo feature point set is extracted from the three-dimensional voxel dataset of the defect structure features.

[0037] It should be noted that extracting a two-dimensional radiation feature point set from thermal anomaly features means using the two-dimensional coordinates of each foreground pixel in the binary mask image of thermal anomaly features as feature points to form a two-dimensional radiation feature point set. Extracting a 3D echo feature point set from a 3D voxel dataset of defect structural features involves taking the 3D voxel dataset of defect structural features as input, establishing a 3×3 local neighborhood window centered on each voxel, calculating the covariance matrix of the scattering intensity values ​​of all voxels within this neighborhood, performing eigenvalue decomposition on the covariance matrix to obtain three eigenvalues, defining the ratio of the product of these three eigenvalues ​​to their sum as the Gaussian curvature of the voxel, generating a Gaussian curvature voxel matrix after traversing all voxels, performing Gaussian filtering on this matrix to smooth noise and enhance the curvature extremum region, setting a curvature threshold, marking voxels with Gaussian curvature greater than the threshold as candidate corner points, and then performing non-maximum suppression on the candidate corner points, i.e., retaining only the voxels with the largest curvature in the local neighborhood as the final corner points, and finally outputting the 3D spatial coordinates of these corner points to form a 3D echo feature point set.

[0038] When the positional deviation of the initial spatial location mapping relationship is detected to be less than the preset spatial deviation threshold, the mapping relationship is determined to be accurate and reliable. No fine optimization based on feature point set is required, and the initial mapping relationship is directly used as the corrected mapping relationship.

[0039] The process of setting the preset spatial deviation threshold is as follows: based on multiple sets of calibration experimental data with known precise spatial correspondence, the average Euclidean distance between thermal anomaly features and defect structure features in each set of data after initial mapping is calculated, the mean and standard deviation of all group distances are statistically analyzed, and the mean plus three times the standard deviation is used as the threshold.

[0040] When performing nearest neighbor search matching on the two-dimensional radiation feature point set and the three-dimensional echo feature point set, firstly, the 64-dimensional descriptor of each feature point in the three-dimensional echo feature point set is used as a query library to construct an index structure based on a randomized KD-tree. The branching factor of the tree is set to 4 to balance search speed and accuracy. Then, for each feature point in the two-dimensional radiation feature point set, its descriptor is input as a query vector into the index structure. The nearest and second nearest neighbor candidate nodes are recursively searched in the KD-tree. The Euclidean distance between the query vector and the candidate node descriptor is calculated, and the values ​​of the nearest and second nearest neighbor distances are output. A distance ratio threshold of 0.75 is set. When the ratio of the nearest neighbor distance to the second nearest neighbor distance is less than this threshold, the nearest neighbor match is considered reliable, and the matching pair is retained as a candidate target matching point; otherwise, it is discarded. All retained matching pairs constitute the target matching point pair set, which is the matching result.

[0041] When correcting the initial mapping relationship based on the matching results, the target matching point pair set is input into the thin plate spline interpolation algorithm model. This model uses the coordinates of the two-dimensional radiation feature points in the target matching point pair as the source control points and the coordinates of the corresponding three-dimensional echo feature points as the target control points. By solving a system of equations containing global affine transformation terms and local radial basis function terms, a nonlinear mapping function from the two-dimensional image coordinate system to the three-dimensional spatial coordinate system is calculated. This mapping function is the corrected mapping relationship. Based on this corrected mapping relationship, any two-dimensional pixel coordinates can be mapped to a precise position in three-dimensional space.

[0042] The distance ratio threshold is based on multiple sets of known correct and incorrect matching calibration data. The distribution of the nearest neighbor distance and the second nearest neighbor distance ratio of all correct matching pairs and the ratio distribution of all incorrect matching pairs are statistically analyzed. The median value between the upper quartile of the correct matching ratio distribution and the lower quartile of the incorrect matching ratio distribution is taken as the threshold.

[0043] Then, when aligning and fusing the thermal anomaly features and the defect structure features according to the modified mapping relationship, firstly, the coordinates of each foreground pixel in the binary mask image of the hot spot region corresponding to the thermal anomaly features are input into the modified mapping relationship. The three-dimensional spatial coordinates corresponding to each two-dimensional pixel coordinate are calculated through the nonlinear mapping function to obtain the mapping point set of the hot spot region in three-dimensional space. At the same time, the three-dimensional voxel dataset corresponding to the defect structure features is used as the base space. Subsequently, a voxel mesh fusion algorithm is used to create a three-dimensional voxel mesh with the same spatial range and resolution as the three-dimensional voxel dataset. Each voxel in the mesh contains two channels, one channel for storing the surface thermal radiation intensity value and the other channel for storing the internal echo reflectivity value.

[0044] Next, each 3D coordinate point in the hot spot region mapping point set is traversed, its voxel location is marked, and the corresponding surface thermal radiation intensity value is written into the surface channel of that voxel. Simultaneously, the internal echo reflectivity value of each voxel in the defect structure feature 3D voxel dataset is directly written into the corresponding voxel's internal channel. For voxels with both surface and internal data, the values ​​from the two channels are concatenated and stored. For voxels with only internal data, the surface channel is set to zero; for voxels with only surface data, the internal channel is set to zero. After all voxels are filled, a 3D voxel dataset containing multi-source physical properties is output; this dataset is the aligned and fused data.

[0045] In step S15, the aligned and fused data needs to be time-segmented to obtain current and historical data. Spatial morphological features of the two data periods are extracted, and first and second damage deterioration features are generated based on these spatial morphological features, including: The aligned and fused data is segmented along the time dimension to obtain data from the current period and data from historical periods; Three-dimensional region growing calculation is performed on the current period data to construct a three-dimensional point cloud model. Principal component analysis is performed on the three-dimensional point cloud model to extract the first damage deterioration feature. Voxel connected component analysis was performed on the historical period data to extract the three-dimensional contours of historical defects; Calculate the volume of the three-dimensional contour of the historical defect. If the volume is greater than a preset volume threshold, perform finite element analysis on the three-dimensional contour of the historical defect to obtain the second damage deterioration feature.

[0046] It should be noted that when the aligned and fused data is segmented according to the time dimension to obtain the current period data and the historical period data, the timestamp labels attached to the aligned and fused data are read. Data with timestamps within the last 30 days are classified as current period fused data, and data with timestamps between the last 180 days and 30 days are classified as historical period fused data, outputting two sets of three-dimensional voxel datasets with clear time attributes.

[0047] Next, when constructing a 3D point cloud model using the current period data for 3D region growth calculations, the current period data refers to the aligned and fused current period data obtained from the previous segmentation. This data is a 3D voxel set with three attributes: surface radiation intensity, internal echo reflectivity, and local dielectric constant. First, all voxels are traversed to extract the local dielectric constant, which is then subtracted from the standard dielectric constant of the healthy leaf matrix. Voxels with an absolute difference greater than a preset dielectric threshold are marked as anomalous dielectric constant voxels. These voxels correspond to abnormal internal media regions and serve as input objects for region growth. Subsequently, the top 5 anomalous dielectric constant voxels with the largest absolute difference are selected as initial seed points, each corresponding to an independent defect region to be grown. For each seed point, its 6-neighbor voxels are traversed. If a neighboring voxel simultaneously satisfies two conditions: it is a voxel with an anomalous dielectric constant, and the difference between its surface radiation intensity and the average radiation intensity of the currently grown region is less than a preset radiation tolerance, then the neighboring voxel is included in the growth range of the current defect. The newly included voxel is then used to expand outwards to a 6-neighborhood, and this process is repeated until no new voxels meet the conditions. After growth is complete, the 3D spatial coordinates of all voxels in each defect region are extracted, and the vertex of its smallest circumscribed convex hull is taken as the defect boundary coordinates. The set of all vertices constitutes a 3D point cloud model, which serves as the current spatial morphological feature.

[0048] Then, when extracting the first damage deterioration feature using principal component analysis on the 3D point cloud model, the 3D point cloud model is dimensionality reduced, and the length change rate and maximum depth gradient of the defect region along the principal axis are calculated. These quantitative indicators are output as the first damage deterioration feature. First, the coordinates of all points in the 3D point cloud model are digitized, i.e., the average coordinates of all points are subtracted. Then, the 3D covariance matrix of the point cloud is calculated, and the three eigenvalues ​​and corresponding eigenvectors of the matrix are solved. The eigenvector corresponding to the maximum eigenvalue is taken as the first principal axis direction of the defect. Each point is projected onto this principal axis direction to obtain one-dimensional projected coordinates. The difference between the maximum and minimum values ​​of all projected coordinates is taken as the current length of the defect along the principal axis. Then, the length of the same defect along the same principal axis direction in the previous inspection cycle is retrieved from the historical inspection data stored in the system. The current length is divided by the length of the previous cycle to obtain the length change rate. At the same time, the depth gradient of each point in the point cloud along the direction perpendicular to the principal axis is calculated: for each point, the difference between its depth direction and that of its neighboring points is divided by the distance between the two points, and the maximum value of this gradient among all points is taken as the maximum depth gradient. The first damage deterioration feature output includes two values: the length change rate and the maximum depth gradient.

[0049] Finally, when calculating the volume of the three-dimensional contour of the historical defect, the number of all voxels contained in the contour is counted. Each voxel corresponds to a fixed physical size. For example, if the side length of each voxel is 1 mm, then the volume of a single voxel is 1 cubic millimeter. The total number of voxels is multiplied by the physical volume of a single voxel to obtain the actual volume value of the three-dimensional contour of the historical defect. If the volume is greater than the preset volume threshold of 500 cubic centimeters, finite element analysis is performed on the three-dimensional contour of the historical defect to obtain the second damage deterioration feature. Specifically, the three-dimensional contour of the historical defect is input into the simplified finite element analysis model, and it is discretized into a tetrahedral mesh by Delaunay triangulation and given composite material properties. After applying a unit tensile load at the defect boundary, the stress distribution is solved, and the ratio of the maximum principal stress to the nominal stress in the far field is taken as the stress concentration factor. At the same time, the difference in the number of voxels between the current contour and the earlier historical data is divided by the time interval to obtain the volume expansion rate. The stress concentration factor and the volume expansion rate are packaged and output as the second damage deterioration feature.

[0050] It is worth noting that the process of setting up the simplified finite element analysis model is as follows: During the system software development phase, secondary development is carried out based on the static analysis solver module in the open-source finite element library. By trimming unnecessary dynamic and thermal analysis functions, the core elasticity and mechanics solver components are retained. Then, a material property database is constructed based on the standard mechanical performance parameters of blade composite materials. Parameters such as longitudinal elastic modulus, transverse elastic modulus, shear modulus, and Poisson's ratio are fixed in the system in the form of configuration files. At the same time, the Delaunay triangulation mesh generation algorithm is encapsulated as an independent preprocessing module. The mesh density parameter is determined to be one element per cubic millimeter after optimization through simulation tests of multiple sets of known defective blades. Finally, all modules are compiled and linked into a dynamic link library file. During the initialization of the UAV inspection system, the loader loads the library file into memory and registers the function entry address. When the system runs, the model can be directly called to complete the calculation of stress and volume change rate by passing in the three-dimensional point cloud model and load parameters.

[0051] It should be noted that the preset volume threshold is based on the destructive detection data of multiple sets of known damaged blades. The lower volume limit of different types of damage, such as delamination, cracks, and debonding, can be clearly identified by the ultrasonic phased array in the early stage. The lower volume limit refers to the minimum physical volume value that can be stably and clearly identified as a defect under the ultrasonic phased array detection conditions for a specific type of damage. Damage below this volume cannot be reliably detected due to weak signals or indistinguishability from noise. The median value of the lower volume limit of all types of damage is taken as the threshold.

[0052] In step S16, based on the first and second damage deterioration characteristics, it is necessary to deduce the evolution process of damage deterioration and assess the overall severity of the wind turbine blade damage evolution to obtain an initial damage severity score, including: By comparing the first damage deterioration feature with the second damage deterioration feature, the deterioration increment feature is extracted, and the spatial expansion value of the defect is calculated. The variation features of the edge contour are extracted from the expanded values, and an evolution trajectory model is constructed; Extract the characteristics reflecting the degree of material fatigue degradation from the evolution trajectory model, and calculate the structural stress change characteristics; The stress variation characteristics of the structure are nonlinearly mapped to obtain an initial damage score.

[0053] It should be noted that when extracting the incremental deterioration feature by comparing the first and second damage deterioration features, the first damage deterioration feature containing the current length change rate and the second damage deterioration feature containing the historical volume expansion rate are respectively used as two row vectors. Each row vector also includes auxiliary parameters related to the feature, such as the current depth gradient and the historical stress concentration coefficient, thus forming two feature vectors with the same dimension. These two feature vectors are then input into the feature difference matrix, and the difference vector is obtained by subtracting the corresponding elements. The Euclidean norm of this difference vector is then calculated as the incremental deterioration feature characterizing the accelerated damage state. Subsequently, based on this incremental deterioration feature, the envelope of the inner and outer boundaries of the blade defect is calculated to obtain the volume difference between the defect region in the current period and in the historical period. This difference is then output as the spatial expansion value.

[0054] Next, when extracting edge contour change features from the spatial expansion values, if the spatial expansion value exceeds a preset 100 cubic centimeters, the system first acquires the boundary voxel point cloud of the defect expansion region. These boundary voxels are the set of voxels located at the outermost layer of the defect, extracted from the 3D region growth results. Then, for each boundary voxel, a spherical neighborhood with a radius of 3 voxels is established centered on it. The 3D coordinates of all boundary points within the neighborhood are collected. Principal component analysis is performed on these coordinates to calculate the three eigenvalues ​​of the covariance matrix. The ratios of the three eigenvalues ​​are combined into a curvature tensor. The difference between the maximum and minimum principal curvature is taken as the Gaussian curvature of that point. After traversing all boundary voxels, the curvature distribution field is obtained. Then, the first-order difference of the curvature distribution field is calculated along the expansion direction, and the position with the largest difference is extracted as the maximum gradient of the edge curvature change. This gradient value is the edge contour change feature.

[0055] Next, the evolution trajectory model is constructed using the edge contour change feature. Specifically, the maximum change gradient of edge curvature at different time points is used as the observation sequence input into the Hidden Markov Model. The model predefines K hidden states, each corresponding to a type of spread direction of the defect in three-dimensional space, such as along the blade spanwise, along the chordwise, or along the thickness direction. The state transition probability matrix is ​​a K×K matrix, where the element in the i-th row and j-th column represents the probability of transitioning from the current state i to the next state j. Then, the Baum-Welch algorithm is used to train the observation sequence, calculating the probability of being in each hidden state at each time step. Based on these probabilities, the maximum likelihood estimation of the state transition probability matrix is ​​performed, iteratively updating the value of each element in the matrix until convergence. After training, the most likely hidden state sequence is decoded based on the latest observation values. This sequence is the spread path of the defect in three-dimensional space, and the path is output as the evolution trajectory model representing the direction of damage development.

[0056] The defect expansion region refers to the difference region between the defect contour in the current period and the defect contour in the historical period. That is, when calculating the spatial expansion value by the three-dimensional convex hull algorithm, the spatial region corresponding to the volume increment obtained by subtracting the volume of the convex hull in the historical period from the volume of the convex hull in the current period is composed of voxels in the convex hull in the current period that are not covered by the convex hull in the historical period.

[0057] Then, when extracting the decay characteristics reflecting material fatigue from the evolution trajectory model, the spatial deviation between the current propagation path and the standard safe path is calculated. The three-dimensional spatial coordinate sequence of each time node on the current propagation path is aligned with the corresponding coordinate sequence on the standard safe path. The Euclidean distance between corresponding points in the two sequences after alignment is calculated, and the root mean square value of the distance between all corresponding points is taken as the spatial deviation value, which is the decay characteristic. This deviation value is used as the decay characteristic. Subsequently, the structural stress change characteristics are calculated based on the decay characteristic. The spatial deviation value is used as the geometric defect amplification factor and substituted into the preset stress concentration factor empirical formula. This formula is obtained by fitting the blade material fatigue test data and outputs a dimensionless stress concentration factor as the structural stress change characteristic. Then, the shear stress increment of the blade main beam area under wind load is calculated. The stress concentration factor is multiplied by the rated wind load corresponding to the extreme wind speed recorded in the local meteorological database to obtain the actual shear stress of the main beam section. Then, the theoretical shear stress at the same position under the healthy state of the blade is subtracted. The difference between the two is the shear stress increment.

[0058] When obtaining the initial damage severity score by nonlinearly mapping the structural stress change characteristics, the shear stress increment is nonlinearly mapped and normalized through a multilayer perceptron network, and an integer between 0 and 100 is output as the initial damage severity score. This quantified value is used to characterize the overall severity of the blade damage evolution.

[0059] It is worth noting that the Hidden Markov Model is configured as follows: During the system development phase, based on historical inspection data of multiple sets of known damage evolution paths, the damage evolution process is divided into three implicit states: initial, acceleration, and instability. The observation sign of each state is obtained by discretizing the gradient of the maximum change in edge curvature into five levels. After multiple rounds of expectation maximization iterative optimization, the converged state transition probability matrix and observation probability matrix are solidified in the model configuration file. When the system runs, the model instantiation can be completed by directly loading the configuration file.

[0060] It is worth noting that the setup process of the multilayer perceptron network is as follows: During the system development phase, multiple sets of historical blade inspection data are collected, with the shear stress increment as the input feature and the damage severity level labeled by experts as the label, to construct a training dataset; a fully connected network structure containing one hidden layer with 64 neurons is adopted, and supervised training is performed using the mean squared error loss function and the Adam optimizer. After multiple rounds of iteration until the loss converges, the trained network weights are solidified in the model file; the system loads this weight file to complete the model instantiation during runtime.

[0061] In step S17, the initial damage severity score needs to be corrected for error compensation to obtain a corrected damage severity score. The leaf condition is then evaluated based on the corrected damage score to obtain a leaf condition evaluation result, including: If the initial damage score is greater than the preset quantization error threshold, the deviation distribution matrix corresponding to the calculation error is extracted; The aligned and fused data is subjected to error compensation processing based on the deviation distribution matrix to obtain corrected fused data. The third deterioration loss feature is extracted based on the corrected fusion data, and the initial damage level score is updated based on the deterioration loss feature to obtain the corrected damage level score.

[0062] It should be noted that if the initial damage score is greater than the preset quantification error threshold, the spatial residual between the predicted damage coordinates and the actual observed coordinates in the historical inspection records of the blade is used as input. Cluster analysis is performed on the residual data. First, the blade surface is divided into 5 equal intervals along the spanwise direction from the leaf root to the leaf tip, and 5 equal intervals along the chordwise direction from the leading edge to the trailing edge, thus dividing the blade surface into 25 grid regions of 5×5. Then, all spatial residuals between the predicted damage coordinates and the actual observed coordinates in the historical inspection records are traversed, and the residuals are classified into the corresponding grid regions according to the location of the predicted coordinates corresponding to each residual. In each grid region, the mean of all residual vectors falling into that region is calculated to obtain a three-dimensional average residual vector. Then, the magnitude of this three-dimensional average residual vector is used as the deviation metric of that grid region. The deviation metrics of the 25 grid regions are arranged in the order of spanwise and chordwise to form a 5x5 deviation distribution matrix. The value of each element in the matrix is ​​the deviation metric of the corresponding grid region. If the initial damage score is not greater than the preset quantization error threshold, the current calculation process is determined to be reliable. There is no need to trigger the error compensation mechanism, and the initial damage score is directly used as the corrected damage score.

[0063] The predicted damage coordinates refer to the three-dimensional spatial coordinates of each node on the defect propagation path decoded by the Hidden Markov Model when the system calculates the initial damage severity score based on the evolution trajectory model and the multilayer perceptron network. These coordinates are the theoretical defect positions deduced by the system based on historical data and current features. The actual observed coordinates refer to the true three-dimensional spatial coordinates of the defect obtained by directly acquiring them through infrared thermal imagers and microwave radar and then spatially aligning them during actual inspections at the same or similar time. These coordinates originate from the measured positions of thermal anomaly features and internal structural features in the aligned and fused data after spatial mapping. The spatial residual between the two is the Euclidean distance difference between the corresponding coordinate points in three-dimensional space.

[0064] The preset quantification error threshold is based on historical inspection data with multiple sets of known accurate assessment results. The absolute value of the deviation between the quantification value of the evolution process and the actual severity of damage is calculated for each set of data. The mean and standard deviation of all absolute values ​​of deviation are calculated, and the mean plus three times the standard deviation is used as the threshold.

[0065] Next, when performing error compensation processing on the aligned fused data according to the deviation distribution matrix to obtain the corrected fused data, the deviation distribution matrix and the aligned fused data are first input into the principal component analysis algorithm to extract the top 3 principal components with the largest variance contribution rate. Then, a compensation weight vector of length 10 is generated according to the loading coefficient of the principal components. This compensation weight vector is then input into the Kalman filter algorithm as a state transition parameter to perform iterative smoothing processing on the original aligned fused data, perform high-frequency signal stripping operation, and output the corrected fused data.

[0066] Then, when extracting the third deterioration loss feature based on the modified fusion data, the modified fusion data is decomposed into the energy distribution state within a specific frequency band, the energy mutation peak is input into the pre-trained support vector regression machine model, a new damage evolution score is calculated through kernel function mapping, and an integer is output as the modified damage degree score.

[0067] This Support Vector Regression (SVR) model maps low-dimensional input data to a high-dimensional feature space using a kernel function. In this high-dimensional space, it searches for an optimal hyperplane that minimizes the deviation of all training sample points from this hyperplane. The training process is as follows: First, multiple sets of historical inspection data are collected. Each set contains a known energy mutation peak and its corresponding true damage evolution score. The true score is manually assigned by experts based on blade damage detection results. All energy mutation peaks are used as input feature vectors, and the corresponding true scores are used as output labels, forming the training sample set. During training, a radial basis function (RBF) is chosen as the kernel function. This kernel function measures the similarity between two input feature vectors by calculating the Euclidean distance between them. The closer the distance between two samples, the closer they are in the mapped high-dimensional space. This kernel function implicitly maps nonlinearly separable data in a low-dimensional space to an infinite-dimensional high-dimensional space, making relationships that cannot be fitted by a linear model in the low-dimensional space linearly separable in the high-dimensional space. This adapts to the nonlinear correlation between energy mutation peaks and damage scores in this scenario, and is less prone to overfitting compared to other kernel functions. The width parameter and regularization penalty coefficient of the kernel function are determined through cross-validation, resulting in a regression model that can predict damage evolution scores based on the input energy mutation peaks. In use, the energy mutation peaks extracted from the current corrected fusion data are input into the model. The model calculates the similarity between the peaks and support vectors using the kernel function, then weights and sums the similarities to output the predicted value. The predicted value is rounded to the nearest integer and used as the corrected damage score.

[0068] It should be noted that the assessment of leaf condition based on the modified damage score to obtain the leaf condition assessment result includes: The modal frequency shift of the blade is calculated based on the modified damage severity score, and the modal frequency shift is compared with the preset structural health benchmark to determine the damage range. If the damage range exceeds the preset critical failure threshold, the remaining life prediction value is calculated and a maintenance priority sequence is generated to finally determine the blade condition assessment result. If it does not exceed the threshold, the corresponding blade condition assessment result is output according to the preset decision mapping table.

[0069] It should be noted that when calculating the modal frequency shift of the blade based on the modified damage degree score, the modified damage degree score is used as input. The rainflow counting method is used to statistically analyze the historical load cycles of the blade, decompose the load amplitude and mean into several complete cycles, and combine the Meinell linear fatigue cumulative damage theory to accumulate the damage score caused by each cycle, and output a dimensionless fatigue cumulative degree index.

[0070] Next, the preset stiffness attenuation matrix is ​​updated using the fatigue accumulation level. This matrix is ​​a 6x6 diagonal matrix, where the diagonal elements correspond to the initial bending stiffness of the six key sections of the blade. The fatigue accumulation level is used as an attenuation factor and multiplied by each diagonal element to obtain the reduced stiffness value, thus updating the entire matrix. Then, the modal frequency change is calculated through structural dynamics analysis. The updated stiffness attenuation matrix is ​​substituted into the overall finite element dynamics equation of the blade, and the eigenvalues ​​of the equation are solved to extract the natural frequency in the first flapping direction. This frequency is subtracted from the reference frequency under healthy blade conditions to obtain the modal frequency offset.

[0071] Then, when comparing the modal frequency offset with the preset structural health benchmark to determine the damage range, the modal frequency offset is input into the interval division mapping table. The mapping table pre-divides the offset value range into three intervals: mild, moderate and severe. Each interval corresponds to a specific upper and lower limit of the offset. The corresponding damage range is determined by comparing which interval the offset falls into.

[0072] If the damage interval exceeds a preset critical failure threshold, the predicted remaining life of the blade is calculated based on the parameters of the damage interval. Specifically, the upper and lower limits of the damage interval and the current cumulative operating time of the blade are used as inputs. The scale and shape parameters are solved using the maximum likelihood estimation method. First, it is assumed that the remaining life of the blade follows a two-parameter Weibull distribution, and a likelihood function containing the scale and shape parameters is constructed. This function expresses the probability of observing the existing data (upper and lower limits of the damage interval and cumulative operating time) under the current parameter values. Then, the logarithm of the likelihood function is taken to obtain the log-likelihood function. The partial derivatives with respect to the scale and shape parameters are calculated and set to zero to obtain two simultaneous equations. Since there is no analytical solution to this system of equations, a numerical solution is performed, i.e., starting from the initial guess value, the predicted value is calculated based on the current parameters in each iteration. The gradient and Hessian matrix of the log-likelihood function are calculated. The parameter values ​​are updated using the product of the gradient and the inverse of the Hessian matrix. This process is repeated iteratively until the parameter change between two adjacent iterations is less than a preset convergence tolerance (e.g., 0.001). The parameter values ​​obtained at this point are the scale and shape parameters under maximum likelihood estimation. After obtaining the scale and shape parameters, the current cumulative operating time of the blade is substituted into the cumulative distribution function of the Weibull distribution to calculate the expected time required for the blade to reach the failure threshold from the current moment, outputting a remaining lifetime prediction value in months. Subsequently, this remaining lifetime prediction value and the prediction values ​​of other blades in the same wind farm are arranged in ascending order of time to generate a maintenance priority sequence. The corresponding evaluation rules are triggered according to the position of the blade in this sequence, and the final blade status evaluation result is output. If the critical failure threshold is not exceeded, the corresponding evaluation result is output according to a preset decision mapping table. This decision mapping table associates different damage intervals with standard maintenance recommendations. For example, a mild interval corresponds to extending the inspection cycle, and a moderate interval corresponds to arranging planned maintenance, directly outputting the corresponding text evaluation conclusion.

[0073] The process of setting the preset stiffness attenuation matrix is ​​as follows: During the system initialization phase, the geometric parameters of six key sections are extracted according to the blade design drawings. The six key sections are the blade root section, the maximum chord length section, the section at one-third of the distance from the blade tip, the section at two-thirds of the distance from the blade tip, the blade tip section, and the transition section, including the area, moment of inertia, and torsional constant of each section. Then, combined with the material properties such as the elastic modulus and shear modulus of the blade composite laminate, the initial bending stiffness of each section is calculated using the Euler-Bernoulli beam theory. The six stiffness values ​​are arranged sequentially on the diagonal, and all off-diagonal elements are set to zero, forming a six-by-six diagonal matrix as the initial stiffness attenuation matrix. The values ​​of this matrix are verified by ground static loading tests and then solidified in the system.

[0074] The process of setting up the finite element dynamic equations is as follows: Based on the structural characteristics and modal analysis requirements of the blade, the blade is simplified into a lumped mass model with six key cross-sectional nodes. Each node only considers the lateral displacement degree of freedom, for a total of six degrees of freedom. Fixed constraints are applied to the blade root nodes, and the boundary conditions are reflected in the equations by deleting corresponding rows and columns or setting them to zero. A six-by-six overall mass matrix is ​​constructed using the lumped mass method, where the diagonal elements are the equivalent mass of each cross-section and the off-diagonal elements are zero. The updated stiffness attenuation matrix is ​​directly used as the overall stiffness matrix. Finally, an undamped free vibration equation is established, ignoring the damping term to simplify the eigenvalue solution.

[0075] The process of setting the preset structural health benchmark is as follows: First, under the healthy condition of the blade at the factory or during its first commissioning, the natural frequency of the first-order flapping direction of the blade is obtained as the benchmark frequency through hammer impact modal test or finite element modal analysis. Then, based on a large amount of historical operation and maintenance data of similar blades in the wind farm, the distribution range of the modal frequency offset of the blade from the healthy state to the appearance of detectable damage is statistically analyzed. The mean offset plus one standard deviation is used as the boundary between the mild and moderate ranges, and the mean offset plus three standard deviations is used as the boundary between the moderate and severe ranges. At the same time, the boundary lines are fine-tuned in combination with the blade design safety margin and industry operation and maintenance specifications. Finally, the upper and lower limit values ​​of the three ranges are fixed in the range division mapping table as the structural health benchmark.

[0076] The process of setting up the interval division mapping table is as follows: First, based on a large amount of historical modal test data of the same type of blades in the wind farm, the distribution pattern of the first-order flapping frequency offset of the blades at different service stages is statistically analyzed. The 95th percentile of the offset distribution of all healthy blades is taken as the upper limit of the mild interval, the 50th percentile of the offset distribution of all blades with detectable cracks is taken as the upper limit of the moderate interval, and the 50th percentile of the offset distribution of all blades with severe delamination is taken as the upper limit of the severe interval. At the same time, the boundaries of each interval are modified by engineering in combination with the blade design safety factor and industry operation and maintenance procedures to ensure that the mild interval corresponds to slight degradation that can continue to operate, the moderate interval corresponds to significant changes that require planned maintenance, and the severe interval corresponds to dangerous conditions that require immediate shutdown and maintenance. Finally, the upper and lower limit values ​​of the three intervals are fixed in the system in the form of a lookup table.

[0077] The second embodiment of the present invention provides a wind turbine blade unmanned aerial vehicle (UAV) inspection system based on multi-source data, including: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the method described above.

[0078] It should be noted that the wind turbine blade UAV inspection system based on multi-source data provided in this embodiment of the invention is used to execute all the process steps of the wind turbine blade UAV inspection method based on multi-source data in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0079] It should be noted that the device 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. Furthermore, in the accompanying drawings of the device embodiments provided by this invention, the connection relationships between modules indicate that they have communication connections, which can be specifically implemented as one or more communication buses or signal lines. Those skilled in the art can understand and implement this without any creative effort.

[0080] 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 unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data, characterized in that, include: Target radiation images and target echo data of wind turbine blades were obtained through drone inspections. Thermal anomaly features on the surface of the wind turbine blades are extracted based on the target radiation image, and defect structure features inside the wind turbine blades are extracted based on the target echo data. The thermal anomaly features and the defect structure features are spatially mapped to obtain an initial mapping relationship; The initial mapping relationship is subjected to deviation correction processing, and the thermal anomaly features and the defect structure features are aligned and fused according to the corrected mapping relationship to obtain aligned and fused data. The aligned and fused data is divided into two periods of data: current and historical data. The spatial morphological features of the two periods of data are extracted, and first and second damage deterioration features are generated based on the spatial morphological features. Based on the first and second damage deterioration characteristics, the evolution process of damage deterioration is deduced and the overall severity of wind turbine blade damage evolution is assessed to obtain an initial damage severity score. The initial damage score is corrected by error compensation to obtain a corrected damage score, and the leaf condition is evaluated based on the corrected damage score to obtain the leaf condition evaluation result.

2. The method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data according to claim 1, characterized in that, The acquisition of target radiation images and target echo data of wind turbine blades through drone inspection includes: The target's radiation image is acquired using an infrared thermal imager mounted on a drone, and the target's echo data is acquired using a microwave radar.

3. The method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data according to claim 1, characterized in that, The step of extracting thermal anomaly features from the surface of the wind turbine blades based on the target radiation image and extracting defect structural features from the interior of the wind turbine blades based on the target echo data includes: Calculate the temperature gradient of the target radiation image, extract hot spot regions from the target radiation image based on the temperature gradient and perform morphological analysis to obtain thermal anomaly features; A phase abrupt change sequence is extracted from the target echo data. If the phase abrupt change sequence meets a preset abrupt change condition, the phase abrupt change sequence is inversely synthesized to generate a scattering cross section matrix. The defect structure characteristics are obtained by spatial mapping of the scattering cross-section matrix.

4. The method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data according to claim 1, characterized in that, The step of spatially mapping the thermal anomaly features and the defect structure features to obtain an initial mapping relationship includes: Extract the set of anomalous centroid coordinates corresponding to the thermal anomaly features and the sequence of structural boundary nodes corresponding to the defect structure features; The set of abnormal centroid coordinates is aligned with the sequence of structural boundary nodes in the main direction to obtain a feature registration vector. The feature registration vector is subjected to projection transformation to establish an initial mapping relationship.

5. The method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data according to claim 1, characterized in that, The initial mapping relationship is corrected for deviation, and the thermal anomaly features and defect structure features are aligned and fused according to the corrected mapping relationship to obtain aligned and fused data, including: Calculate the positional deviation of the initial mapping relationship. If the positional deviation exceeds a preset spatial deviation threshold, extract the two-dimensional radiation feature point set of the thermal anomaly feature and the three-dimensional echo feature point set of the defect structure feature. Otherwise, use the initial mapping relationship as a corrected mapping relationship. The two-dimensional radiation feature point set and the three-dimensional echo feature point set are matched by nearest neighbor search to obtain the matching result, and the initial mapping relationship is corrected to obtain the corrected mapping relationship; The thermal anomaly features and the defect structure features are aligned and fused according to the modified mapping relationship to obtain aligned and fused data.

6. The method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data according to claim 1, characterized in that, The process involves time-segmenting the aligned and fused data to obtain current and historical data, extracting spatial morphological features from the two data periods, and generating first and second damage deterioration features based on these spatial morphological features, including: The aligned and fused data is segmented along the time dimension to obtain data from the current period and data from historical periods; Three-dimensional region growing calculation is performed on the current period data to construct a three-dimensional point cloud model. Principal component analysis is performed on the three-dimensional point cloud model to extract the first damage deterioration feature. Voxel connected component analysis was performed on the historical period data to extract the three-dimensional contours of historical defects; Calculate the volume of the three-dimensional contour of the historical defect. If the volume is greater than a preset volume threshold, perform finite element analysis on the three-dimensional contour of the historical defect to obtain the second damage deterioration feature.

7. The method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data according to claim 1, characterized in that, Based on the first and second damage deterioration characteristics, the evolution process of damage deterioration is deduced and the overall severity of wind turbine blade damage evolution is assessed to obtain an initial damage severity score, including: By comparing the first damage deterioration feature with the second damage deterioration feature, the deterioration increment feature is extracted, and the spatial expansion value of the defect is calculated. The variation features of the edge contour are extracted from the expanded values, and an evolution trajectory model is constructed; Extract the characteristics reflecting the degree of material fatigue degradation from the evolution trajectory model, and calculate the structural stress change characteristics; The stress variation characteristics of the structure are nonlinearly mapped to obtain an initial damage score.

8. The method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data according to claim 1, characterized in that, The step of performing error compensation correction on the initial damage severity score to obtain a corrected damage severity score includes: If the initial damage score is greater than the preset quantization error threshold, the deviation distribution matrix corresponding to the calculation error is extracted; The aligned and fused data is subjected to error compensation processing based on the deviation distribution matrix to obtain corrected fused data. The third deterioration loss feature is extracted based on the corrected fusion data, and the initial damage level score is updated based on the deterioration loss feature to obtain the corrected damage level score.

9. The method for unmanned aerial vehicle (UAV) inspection of wind turbine blades based on multi-source data according to claim 1, characterized in that, The step of assessing the leaf condition based on the modified damage score to obtain the leaf condition assessment result includes: The modal frequency shift of the blade is calculated based on the modified damage severity score, and the modal frequency shift is compared with the preset structural health benchmark to determine the damage range. If the damage range exceeds the preset critical failure threshold, the remaining life prediction value is calculated and a maintenance priority sequence is generated to finally determine the blade condition assessment result. If it does not exceed the threshold, the corresponding blade condition assessment result is output according to the preset decision mapping table.

10. A wind turbine blade unmanned aerial vehicle (UAV) inspection system based on multi-source data, 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 9.

Images