Nondestructive flaw detection method and system

Abstract

The invention relates to a nondestructive flaw detection method and system, and the method comprises the steps: carrying out the dynamic coupling wave mode decomposition of a response signal, and obtaining a basic wave mode component and an energy ratio; aiming at a plurality of discrete time-frequency scale units, constructing a multi-scale time-frequency energy flow graph cluster based on the basic wave mode component and the energy ratio; calculating the direction gradient of each energy flow graph in the energy flow graph cluster by taking the energy amplitude as the weight to obtain a gradient vector field; determining a suspicious region according to topological structure characteristics of the gradient vector field and an abnormal threshold value of a wave mode energy ratio in the energy flow graph cluster; determining the boundary of the suspicious region based on the trajectory of the gradient vector field; performing depth distribution inversion calculation on the suspicious region according to the energy attenuation gradient of the multi-scale energy flow diagram and the attenuation of each basic wave mode energy in the depth direction; and performing regularization tomography by using the energy ratio and spatial distribution of each fundamental wave mode as multi-channel input, and reconstructing the three-dimensional geometrical morphology of the suspicious region.

Description

Technical Field

[0001] The present invention relates to the field of flaw detection, and in particular to a non-destructive flaw detection method and system. Background Art

[0002] Non-destructive testing (NDT) plays an irreplaceable role in ensuring safe equipment operation, improving product quality, and preventing catastrophic accidents. It inspects and tests the internal and surface structure, properties, and condition of test pieces without damaging or affecting their performance. Traditional NDT techniques, including ultrasonic testing, radiographic testing, magnetic particle testing, penetrant testing, and eddy current testing, constitute a crucial technical support for industrial safety and quality control. Ultrasonic testing, due to its strong penetration, sensitivity to area-type defects, and relatively low cost, has been widely used in the inspection of welds, forgings, and composite materials. While advanced methods such as phased array ultrasound, time-of-flight diffraction (TOFD), and synthetic aperture focusing (SAFT) have been developed to improve inspection accuracy and image quality, elastic waves propagating through complex structures produce multiple wave modes, resulting in conversions and superpositions, as well as complex interactions with defects. This makes the received signals extremely complex, making it difficult to accurately extract characteristic information related to the defects. Furthermore, the precise positioning, sizing, depth assessment, and three-dimensional topography of defects have always been challenging research areas in ultrasonic testing. Traditional evaluation methods based on signal amplitude or time of flight are susceptible to interference from factors such as material inhomogeneity, structural geometry, and detection coupling conditions, resulting in insufficient reliability and accuracy in test results. Furthermore, the key to improving detection capabilities lies in effectively utilizing multimodal ultrasonic information to fully explore the defect information contained in wavefield energy distribution, propagation paths, and mode conversions, and presenting this information to inspectors in an intuitive manner to assist in comprehensive defect assessment. Summary of the Invention

[0003] In order to solve the above problems, the present invention provides a non-destructive testing method, which includes the following steps: An excitation source generates elastic waves in a predetermined area of ​​the measured object, and a sensor array collects response signals; the response signals are subjected to dynamic coupled wave mode decomposition to obtain fundamental wave mode components and energy proportions; and a multi-scale time-frequency energy flow graph cluster based on the fundamental wave mode components and energy proportions is constructed for multiple discrete time-frequency scale units; Calculating the directional gradient of each energy flow graph in the energy flow graph cluster with the energy amplitude as the weight to obtain a gradient vector field; determining a suspicious area based on the topological structural characteristics of the gradient vector field and the abnormal threshold of the wave mode energy ratio in the energy flow graph cluster; and determining the boundary of the suspicious area based on the trajectory of the gradient vector field; The depth distribution inversion calculation of the suspicious area is performed based on the energy attenuation gradient of the multi-scale energy flow diagram and the depth-direction attenuation of the energy of each fundamental wave mode; and the energy proportion and spatial distribution of each fundamental wave mode are used as multi-channel input to perform regularized tomography and reconstruct the three-dimensional geometric shape of the suspicious area.

[0004] Optionally, the dynamic coupled wave mode decomposition is performed on the response signal to obtain the fundamental wave mode component and energy proportion, specifically: The response signal collected by each sensor is converted to the frequency-wavenumber domain using a two-dimensional fast Fourier transform. In the frequency-wavenumber domain, the isolation band width is determined based on the theoretical dispersion curves of different wave modes. The fundamental wave mode components are separated using multiple predefined two-dimensional bandpass filter banks that match the theoretical dispersion curve shapes. The amplitude spectrum of the fundamental wave mode component is integrated in the frequency-wavenumber domain to obtain the total energy of the fundamental wave mode component, and the proportion of the fundamental wave mode component in the total energy of all fundamental wave modes at the current acquisition point is calculated to obtain the energy proportion.

[0005] Optionally, the determining of the suspicious area based on the topological structural characteristics of the gradient vector field and the abnormal threshold of the wave mode energy ratio in the energy flow graph cluster is specifically as follows: Calculate the divergence field and curl field of each weighted directional gradient field; The area with a value lower than the negative threshold in the divergence field is regarded as the energy convergence area, and the area with a value higher than the positive threshold is regarded as the energy divergence area; the area with a value higher than the amplitude threshold in the curl field is regarded as the energy rotation area; The spatial union of the energy convergence region, the energy divergence region and the energy rotation region is regarded as the topological singular region; Calculating the normalized deviation of the energy proportion of each fundamental wave mode at each spatial position relative to the mean energy proportion of the mode in the entire measured area in the energy flow graph cluster; When the energy ratio normalized deviations of at least two fundamental wave modes are greater than a preset value, the spatial position is marked as an energy ratio abnormal point; The spatial intersection of the topological singular region and all energy proportion abnormal point sets is taken as the suspicious region.

[0006] Optionally, the determining the boundary of the suspicious area based on the trajectory of the gradient vector field is specifically: Inside each suspicious area, the point with the largest modulus of the gradient vector field is selected as the seed point; Starting from the seed point, the fourth-order Runge-Kutta method is used to bidirectionally trace gradient streamlines and equipotential lines along the direction of the gradient vector field and the direction perpendicular to the gradient vector field until the gradient modulus of the streamlines or equipotential lines is lower than the background gradient level predetermined based on the statistical characteristics of the signal of the entire object under test; the outermost envelope of all closed gradient streamlines and equipotential line segments or those truncated at the edge of the suspicious area is used as the boundary of the suspicious area.

[0007] Optionally, the depth distribution inversion calculation of the suspicious area is performed based on the energy attenuation gradient of the multi-scale energy flow diagram and the depth-direction attenuation of the energy of each fundamental wave mode, specifically: Extract the energy value sequence corresponding to different center frequencies in the multi-scale time-frequency energy flow graph cluster for each point in each suspicious area, and calculate the attenuation gradient of the energy value with frequency; Establishing a depth-attenuation lookup table for each fundamental wave mode when there are no defects in the measured object material, the lookup table containing theoretical attenuation characteristics of wave mode energy as propagation depth increases at different frequencies; constructing a cost function based on the weighted mean square error between a theoretical energy attenuation gradient calculated from the theoretical attenuation characteristics and the actually observed multi-scale energy attenuation gradient; The conjugate gradient method is used to minimize the cost function to obtain the optimal depth value of the point; the above calculation is performed on all points in the suspicious area to obtain the depth distribution.

[0008] Optionally, the energy proportion and spatial distribution of each fundamental wave mode are used as multi-channel input to perform regularized tomography to reconstruct the three-dimensional geometric shape of the suspicious area, specifically: The energy proportion of the N fundamental wave modes decomposed at each sensor collection point is used as the N-dimensional feature vector of the spatial position; The N-dimensional feature vectors of all the collected points in the suspicious area and the adjacent areas are used to form N independent two-dimensional data channels; An algebraic reconstruction technique is used as the basic tomography algorithm, and a total variation regularization term is utilized to preserve edges and suppress noise. The tomographic imaging voxel grid size is set to 1 / 5 of the minimum spacing of the sensor array. The following optimization problem is solved iteratively: minimizing the data fidelity term plus the regularization coefficient multiplied by the total variation regularization term, where the data fidelity term is the L2 norm difference between the theoretical value and the actual observed value of the two-dimensional data of each channel calculated from the reconstructed three-dimensional voxels through a preset physical acoustic model, and λ is the regularization coefficient. Iteration is performed until the change in the reconstructed voxel model between two consecutive iterations is less than a preset convergence value, thereby obtaining the three-dimensional geometric shape of the suspicious area.

[0009] The present invention also provides a non-destructive testing system, which includes the following modules: An acquisition and decomposition module is configured to generate elastic waves in a predetermined area of ​​the measured object by an excitation source, and to collect response signals from a sensor array; perform dynamic coupled wave mode decomposition on the response signals to obtain fundamental wave mode components and energy proportions; and construct a multi-scale time-frequency energy flow graph cluster based on the fundamental wave mode components and energy proportions for multiple discrete time-frequency scale units; a preliminary identification module for calculating the directional gradient of each energy flow graph in the energy flow graph cluster, weighted by energy amplitude, to obtain a gradient vector field; determining a suspicious area based on the topological structural characteristics of the gradient vector field and an abnormal threshold of the wave mode energy ratio in the energy flow graph cluster; and determining the boundary of the suspicious area based on the trajectory of the gradient vector field; The morphology extraction module is used to perform depth distribution inversion calculation of the suspicious area based on the energy attenuation gradient of the multi-scale energy flow diagram and the depth-direction attenuation of the energy of each fundamental wave mode; and use the energy proportion and spatial distribution of each fundamental wave mode as multi-channel input to perform regularized tomography and reconstruct the three-dimensional geometric morphology of the suspicious area.

[0010] Optionally, the dynamic coupled wave mode decomposition is performed on the response signal to obtain the fundamental wave mode component and energy proportion, specifically: The response signal collected by each sensor is converted to the frequency-wavenumber domain using a two-dimensional fast Fourier transform. In the frequency-wavenumber domain, the isolation band width is determined based on the theoretical dispersion curves of different wave modes. The fundamental wave mode components are separated using multiple predefined two-dimensional bandpass filter banks that match the theoretical dispersion curve shapes. The amplitude spectrum of the fundamental wave mode component is integrated in the frequency-wavenumber domain to obtain the total energy of the fundamental wave mode component, and the proportion of the fundamental wave mode component in the total energy of all fundamental wave modes at the current acquisition point is calculated to obtain the energy proportion.

[0011] Optionally, the determining of the suspicious area based on the topological structural characteristics of the gradient vector field and the abnormal threshold of the wave mode energy ratio in the energy flow graph cluster is specifically as follows: Calculate the divergence field and curl field of each weighted directional gradient field; The area with a value lower than the negative threshold in the divergence field is regarded as the energy convergence area, and the area with a value higher than the positive threshold is regarded as the energy divergence area; the area with a value higher than the amplitude threshold in the curl field is regarded as the energy rotation area; The spatial union of the energy convergence region, the energy divergence region and the energy rotation region is regarded as the topological singular region; Calculating the normalized deviation of the energy proportion of each fundamental wave mode at each spatial position relative to the mean energy proportion of the mode in the entire measured area in the energy flow graph cluster; When the energy ratio normalized deviations of at least two fundamental wave modes are greater than a preset value, the spatial position is marked as an energy ratio abnormal point; The spatial intersection of the topological singular region and all energy proportion abnormal point sets is taken as the suspicious region.

[0012] Optionally, the determining the boundary of the suspicious area based on the trajectory of the gradient vector field is specifically: Inside each suspicious area, the point with the largest modulus of the gradient vector field is selected as the seed point; Starting from the seed point, the fourth-order Runge-Kutta method is used to bidirectionally trace gradient streamlines and equipotential lines along the direction of the gradient vector field and the direction perpendicular to the gradient vector field until the gradient modulus of the streamlines or equipotential lines is lower than the background gradient level predetermined based on the statistical characteristics of the signal of the entire object under test; the outermost envelope of all closed gradient streamlines and equipotential line segments or those truncated at the edge of the suspicious area is used as the boundary of the suspicious area.

[0013] Optionally, the depth distribution inversion calculation of the suspicious area is performed based on the energy attenuation gradient of the multi-scale energy flow diagram and the depth-direction attenuation of the energy of each fundamental wave mode, specifically: Extract the energy value sequence corresponding to different center frequencies in the multi-scale time-frequency energy flow graph cluster for each point in each suspicious area, and calculate the attenuation gradient of the energy value with frequency; Establishing a depth-attenuation lookup table for each fundamental wave mode when there are no defects in the measured object material, the lookup table containing theoretical attenuation characteristics of wave mode energy as propagation depth increases at different frequencies; constructing a cost function based on the weighted mean square error between a theoretical energy attenuation gradient calculated from the theoretical attenuation characteristics and the actually observed multi-scale energy attenuation gradient; The conjugate gradient method is used to minimize the cost function to obtain the optimal depth value of the point; the above calculation is performed on all points in the suspicious area to obtain the depth distribution.

[0014] Optionally, the energy proportion and spatial distribution of each fundamental wave mode are used as multi-channel input to perform regularized tomography to reconstruct the three-dimensional geometric shape of the suspicious area, specifically: The energy proportion of the N fundamental wave modes decomposed at each sensor collection point is used as the N-dimensional feature vector of the spatial position; The N-dimensional feature vectors of all the collected points in the suspicious area and the adjacent areas are used to form N independent two-dimensional data channels; An algebraic reconstruction technique is used as the basic tomography algorithm, and a total variation regularization term is utilized to preserve edges and suppress noise. The tomographic imaging voxel grid size is set to 1 / 5 of the minimum spacing of the sensor array. The following optimization problem is solved iteratively: minimizing the data fidelity term plus the regularization coefficient multiplied by the total variation regularization term, where the data fidelity term is the L2 norm difference between the theoretical value and the actual observed value of the two-dimensional data of each channel calculated from the reconstructed three-dimensional voxels through a preset physical acoustic model, and λ is the regularization coefficient. Iteration is performed until the change in the reconstructed voxel model between two consecutive iterations is less than a preset convergence value, thereby obtaining the three-dimensional geometric shape of the suspicious area.

[0015] The present invention performs dynamic coupled wave mode decomposition on the response signals collected by the sensor array to obtain the fundamental wave mode components and energy proportions, and constructs a multi-scale time-frequency energy flow graph cluster for multiple discrete time-frequency scale units, which can more carefully reveal the propagation characteristics of elastic waves in the measured object and the interaction mechanism with defects. In addition, the directional gradient of each energy flow graph in the energy flow graph cluster is calculated with the energy amplitude as the weight to obtain the gradient vector field, and the suspicious area is determined by combining its topological structure characteristics with the abnormal threshold of the wave mode energy proportion, and the boundary of the suspicious area is determined based on the trajectory of the gradient vector field, thereby improving the accuracy of defect positioning and the clarity of boundary identification. By performing depth distribution inversion calculation on the suspicious area based on the energy attenuation gradient of the multi-scale energy flow graph and the depth-direction attenuation of the energy of each fundamental wave mode, an effective assessment of the defect depth is achieved, and the three-dimensional geometric shape of the suspicious area can be reconstructed, providing an intuitive basis for the nature judgment and risk assessment of the defect. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 This is a flow chart of Example 1; Figure 2 is a schematic diagram of the frequency-wavenumber spectrum; Figure 3 is a schematic diagram of a bandpass filter; Figure 4 Schematic diagram of a topologically singular region; Figure 5 Schematic diagram of abnormal energy ratio points; Figure 6 Schematic diagram of the final suspicious area. DETAILED DESCRIPTION

[0017] In order to enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions in the embodiments of this specification will be clearly and completely described below in conjunction with the drawings in the embodiments of this specification. Obviously, the embodiments described are only part of the embodiments of this application, not all of the embodiments. Based on the embodiments of this specification, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of this application.

[0018] Those skilled in the art will appreciate that the embodiments of this specification may be provided as methods, systems, or computer program products. Therefore, the embodiments of this specification may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Furthermore, the embodiments of this specification may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0019] The present invention is described below with reference to a specific embodiment.

[0020] like Figure 1 The non-destructive testing method shown includes the following steps: S1, an excitation source generates elastic waves in a predetermined area of ​​the measured object, and a sensor array collects response signals; the response signals are subjected to dynamic coupled wave mode decomposition to obtain fundamental wave mode components and energy proportions; and a multi-scale time-frequency energy flow graph cluster based on the fundamental wave mode components and energy proportions is constructed for multiple discrete time-frequency scale units; By applying pulse or continuous wave excitation at specific locations on the surface of the object being measured, the elastic waves required for detection are generated, and multiple sensors arranged in an optimal manner synchronously record the full wavefield response signal containing rich propagation information. In order to separate the multiple fundamental wave modes aliased in the signal and quantify their respective energy contributions, an adaptive time-frequency analysis method based on wavelet packet transform is used in one embodiment to automatically select the optimal decomposition basis based on the local characteristics of the signal, thereby effectively analyzing the energy distribution of different wave mode components as they change with time and frequency. Through techniques such as variational mode decomposition or Hilbert-Huang transform, the instantaneous frequency and amplitude characteristics of the signal are extracted, and then a multi-scale energy distribution map that can dynamically reflect the continuous flow of energy in time and space dimensions is constructed to obtain an energy flow graph cluster.

[0021] S2, calculating the directional gradient of each energy flow graph in the energy flow graph cluster with the energy amplitude as the weight to obtain a gradient vector field; determining a suspicious area based on the topological structure characteristics of the gradient vector field and the abnormal threshold of the wave mode energy ratio in the energy flow graph cluster; and determining the boundary of the suspicious area based on the trajectory of the gradient vector field; When calculating the directional gradient of the energy flow graph, in one embodiment, a gradient calculation method based on a structural tensor is employed. Energy amplitude is introduced as an adaptive weighting factor to enhance the significance of gradient features in high-energy regions, thereby generating a gradient vector field that better reflects the primary energy propagation path. In another embodiment, the basic gradient components are calculated at each location in each energy flow graph, for example, using the Sobel operator to obtain the gradients in the x and y directions. Then, at each location, the energy amplitude of that point is used as a weight and multiplied by the calculated basic gradient components to obtain the energy amplitude-weighted directional gradient at that point. These weighted directional gradients calculated for all points in the energy flow graph are aggregated to obtain the gradient vector field. When identifying suspicious areas, in one embodiment, machine learning-based pattern recognition techniques are employed. For example, a deep neural network model is trained, taking into account various topological properties of the gradient vector field, such as divergence, curl, and the curvature and density of the vector lines. Combined with the statistical characteristics of wave mode energy proportions extracted from the energy flow graph cluster, the model automatically learns and identifies abnormal patterns, thereby identifying suspicious damage areas. The precise boundaries of the suspicious area can be determined through image segmentation algorithms, such as the watershed algorithm or the Markov random field model, which uses the intensity and direction information of the gradient vector field to guide the generation and optimization of the boundary so that it accurately fits the actual contour of the defect.

[0022] S3, based on the energy attenuation gradient of the multi-scale energy flow diagram and the depth-direction attenuation of the energy of each fundamental wave mode, the depth distribution inversion calculation of the suspicious area is performed; and the energy proportion and spatial distribution of each fundamental wave mode are used as multi-channel input to perform regularized tomography to reconstruct the three-dimensional geometric shape of the suspicious area.

[0023] When performing depth distribution inversion of a suspicious area, in one embodiment, a three-dimensional finite element simulation model is established that incorporates the acoustic properties of the material and the geometric parameters of the structure. The depth parameters of the defects in the model are iteratively adjusted, and the simulated surface wave field energy attenuation characteristics are compared with the experimentally observed multi-scale energy flow diagram attenuation gradient. A global optimization algorithm, such as a genetic algorithm or particle swarm optimization, is then used to search for the best-matching depth value. To reconstruct the three-dimensional geometric morphology of the suspicious area, in one embodiment, a tomographic imaging inversion framework based on Bayesian theory is employed. The energy proportions and spatial distribution information of different fundamental wave modes are incorporated into the inversion process as multivariate prior information. By solving the posterior probability distribution that takes into account the data noise model and model parameter uncertainty, a probabilistically meaningful three-dimensional morphology of the defect is reconstructed. Methods such as Gaussian process regression can be used to interpolate and fuse sparse or incomplete data to improve imaging quality.

[0024] In an optional embodiment, the dynamic coupled wave mode decomposition is performed on the response signal to obtain the fundamental wave mode component and energy proportion, specifically: The response signal collected by each sensor is converted to the frequency-wavenumber domain using a two-dimensional fast Fourier transform. In the frequency-wavenumber domain, the isolation band width is determined based on the theoretical dispersion curves of different wave modes. The fundamental wave mode components are separated using multiple predefined two-dimensional bandpass filter banks that match the theoretical dispersion curve shapes. The amplitude spectrum of the fundamental wave mode component is integrated in the frequency-wavenumber domain to obtain the total energy of the fundamental wave mode component, and the proportion of the fundamental wave mode component in the total energy of all fundamental wave modes at the current acquisition point is calculated to obtain the energy proportion.

[0025] Specifically, suppose a composite laminate with a thickness of 5 mm is being tested. The excitation source generates a broadband ultrasonic pulse somewhere on the plate, and a group of sensors arranged in a linear array, for example, 64, receives the complex response signal after the elastic wave propagates in the material. For each time series signal received by the sensor, such as a data containing 2048 sampling points, it is combined with the spatial arrangement information of the entire sensor array and a two-dimensional fast Fourier transform is performed together. The transformation converts the original time-space domain signal into the frequency wavenumber domain, and the energy of the signal is distributed according to its frequency component and the spatial frequency, i.e., the wavenumber, along the direction of the sensor array. For example, at a frequency of 300 kHz, energy may be observed to be concentrated in areas with wavenumbers of 0.5 radians per millimeter and 1.2 radians per millimeter. Different wavenumbers correspond to different wave propagation modes.

[0026] In the frequency wavenumber domain, based on the dispersion curves of Lamb waves and other fundamental wave modes known to composite materials or obtained through theoretical calculations, such as the theoretical wavenumbers corresponding to the S0 mode and the A0 mode at different frequencies. For example, theoretically, the wavenumber of the A0 mode at 300 kHz should be 1.2 radians per millimeter, and the S0 mode corresponds to 0.5 radians per millimeter, such as Figure 2 As shown. Using the theoretical dispersion curve, a certain isolation band width is set, such as selecting plus or minus 10% of the theoretical wave value as the boundary, to obtain a series of two-dimensional bandpass filters. Each filter matches the dispersion curve shape of a specific fundamental wave mode, such as Figure 3As shown in Figure 3 , applying these filters to data in the frequency-wavenumber domain effectively separates the aliased fundamental wave mode components, such as the S0 and A0 mode signal components, from the total signal. After separation, for each fundamental wave mode, such as the S0 mode, its amplitude spectrum in the frequency-wavenumber domain is integrated to obtain a numerical value. For example, the total energy of the S0 mode is 80. Similarly, the total energy of the A0 mode might be 120. At the sensor acquisition point currently under consideration, the total energy of all fundamental wave modes is 200. The S0 mode contributes 0.4% of the energy, and the A0 mode contributes 0.6%.

[0027] In an optional embodiment, for a plurality of discrete time-frequency scale units, a multi-scale time-frequency energy flow graph cluster based on the fundamental wave mode components and energy proportions is constructed, specifically: Set the time window length set T = {t1, t2, ..., tM} and the center frequency set F = {f1, f2, ..., fN} as the parameters of the discrete time-frequency scaling unit, where ti+1 = ti*Qt, fj+1 = fj*Qf, and Qt and Qf are constant coefficients greater than 1; For each sensor acquisition point and each time-frequency scale unit defined by a specific time window length tk and a center frequency fl, the fundamental wave mode component signals decomposed at the acquisition point are bandpass filtered and time-windowed, where the center frequency of the bandpass filter is fl, the bandwidth is fl / 4, and the time-windowing window function is a Hanning window with a length of tk; Calculate the energy value of each fundamental wave mode within the time-frequency scale unit and its proportion in the total energy of the unit; map the energy proportion of the fundamental wave mode of each sensor acquisition point under the time-frequency scale unit to the spatial coordinates of the acquisition point to obtain a two-dimensional energy flow diagram corresponding to the time-frequency scale unit and the fundamental wave mode; The energy flow graphs corresponding to all time-frequency scale units and all fundamental wave modes constitute the multi-scale time-frequency energy flow graph cluster.

[0028] It is necessary to first determine a series of discrete time-frequency observation scales. For the time window length, start with a shorter length. For example, set the first time window t1 to 10 microseconds. The subsequent time window lengths are equal to tk multiplied by a constant coefficient Qt greater than one, such as Qt is 2, and then successively obtain 20 microseconds, 40 microseconds, and so on. Assume that a total of five such time window lengths are set. Similarly, for the center frequency, a frequency set must also be determined. For example, starting with the first center frequency f1 of 50 kHz, the subsequent center frequencies are equal to fl multiplied by a constant coefficient Qf greater than one, such as Qf is 1.2, thus obtaining 60 kHz, 72 kHz, and so on. Assume that a total of four such center frequencies are set. Then, each specific time window length tk and a specific center frequency fl constitute an independent discrete time-frequency scale unit, serving as an observation window of different resolutions.

[0029] For each sensor acquisition point and each defined time-frequency scale unit, the fundamental wave mode component signals decomposed at that acquisition point, such as the S zero-mode Lamb wave and A zero-mode Lamb wave components, are processed. For a time-frequency scale unit with a 10-microsecond time window and a center frequency of 50 kHz, the S zero-mode signal component is extracted and bandpass filtered with a center frequency of 50 kHz and a bandwidth set to one-quarter of the center frequency, i.e., 12.5 kHz. A Hanning window is then used to intercept a 10-microsecond portion of the filtered signal. The energy of the S zero-mode within this 10-microsecond time window and near the 50 kHz frequency band is calculated. Similarly, the energy of the A zero-mode and other possible fundamental wave modes within this time-frequency scale unit is calculated. The energy contribution of each fundamental wave mode to the total energy of all fundamental wave modes within the current time-frequency scale unit is then calculated. The energy contribution of a specific fundamental wave mode, such as the S-zero mode, at the current sensor location and within the current 10 microsecond, 50 kHz time-frequency scale unit is assigned to the spatial coordinates of the sensor. The same operation is performed for all sensor acquisition points, resulting in a two-dimensional energy flow diagram corresponding to the S-zero mode at the current time-frequency scale. Repeating this process for all M time window lengths, N center frequencies, and all identified fundamental wave modes results in a cluster of multi-scale time-frequency energy flow diagrams. For example, if there are five time scales, four frequency scales, and two fundamental wave modes, the resulting cluster will be 40 different two-dimensional energy flow diagrams.

[0030] In an optional embodiment, the determining of the suspicious area based on the topological structural characteristics of the gradient vector field and the abnormal threshold of the wave mode energy ratio in the energy flow graph cluster is specifically as follows: Calculate the divergence field and curl field of each weighted directional gradient field; The area with a value lower than the negative threshold in the divergence field is regarded as the energy convergence area, and the area with a value higher than the positive threshold is regarded as the energy divergence area; the area with a value higher than the amplitude threshold in the curl field is regarded as the energy rotation area; The spatial union of the energy convergence region, the energy divergence region and the energy rotation region is regarded as the topological singular region; Calculating the normalized deviation of the energy proportion of each fundamental wave mode at each spatial position relative to the mean energy proportion of the mode in the entire measured area in the energy flow graph cluster; When the energy ratio normalized deviations of at least two fundamental wave modes are greater than a preset value, the spatial position is marked as an energy ratio abnormal point; The spatial intersection of the topological singular region and all energy proportion abnormal point sets is taken as the suspicious region.

[0031] After obtaining the detailed direction and intensity map of the surface wave energy flow of the measured component, that is, the gradient vector field, the divergence of this field is calculated. For example, in a small area of ​​a steel plate, if the calculated divergence value is particularly low, for example, lower than a set negative threshold of -0.8, the energy is concentrated in this area, which is recorded as an energy convergence area. This may be because there is a defect that absorbs or scatters energy. On the contrary, if the divergence value is particularly high, such as higher than the positive threshold of 0.8, the energy diffuses outward from here to form an energy divergence area. At the same time, the curl of the gradient field is also calculated, which reflects the degree of rotation of the energy flow. If the curl value at a certain place exceeds an amplitude threshold, such as 1.2, it means that the energy has undergone obvious vortex or rotation here, which is recorded as an energy rotation area, which is common at the edge of a defect or a complex interaction point. All identified energy convergence areas, energy divergence areas, and energy rotation areas are spatially merged to form a preliminary set of topologically singular regions. These areas are areas where the energy flow behavior is abnormal, such as Figure 4 shown.

[0032] In order to improve the accuracy, the abnormal situation of the energy ratio of the wave mode in the energy flow diagram cluster is introduced as another basis for judgment. Assume that two Lamb wave modes are mainly analyzed in the detection, and they are called mode A and mode B. For the entire measured area, the energy ratio of mode A may be 0.6 on average, and mode B is 0.4. The normalized deviation of the energy ratio of mode A and mode B at each point relative to their respective global average ratios is calculated point by point. For example, at a specific point, the energy ratio of mode A reaches 0.9, and its normalized deviation may be 4.5, while the energy ratio of mode B drops to 0.1, and its normalized deviation is -5.0. If it is set that the energy ratio normalized deviation of at least two basic wave modes is greater than a preset value, such as 3.0, then this point is marked as an energy ratio abnormal point, such as Figure 5As shown. That is, the wave energy distribution at this point is significantly different from the normal situation, and it is likely to be affected by a local defect. The obtained topological singular region is spatially compared with all these energy ratio abnormal point sets. Only those regions that belong to the topological singular region and are marked as energy ratio abnormal points at the same time are finally determined to be the real suspicious regions. For example, if a point is both an energy convergence area and shows that the energy ratio deviation of mode A and mode C exceeds 3, then it is confirmed as a suspicious defect location, such as Figure 6 shown.

[0033] In an optional embodiment, the determining the boundary of the suspicious area based on the trajectory of the gradient vector field is specifically: Inside each suspicious area, the point with the largest modulus of the gradient vector field is selected as the seed point; Starting from the seed point, the fourth-order Runge-Kutta method is used to bidirectionally trace gradient streamlines and equipotential lines along the direction of the gradient vector field and the direction perpendicular to the gradient vector field until the gradient modulus of the streamlines or equipotential lines is lower than the background gradient level predetermined based on the statistical characteristics of the signal of the entire object under test; the outermost envelope of all closed gradient streamlines and equipotential line segments or those truncated at the edge of the suspicious area is used as the boundary of the suspicious area.

[0034] Specifically, suppose a suspected area is identified on a composite material panel under inspection, perhaps encompassing a 10 cm x 10 cm square. Within this square, the wave energy flow gradient vector field has been calculated. The direction of each point in this vector field represents the direction of energy flow, while the magnitude, or modulus, of the vector represents its intensity. To precisely determine the boundaries of the suspected area, the point with the largest gradient vector modulus within this 100 cm² area is first searched. For example, at coordinates X = 5 cm, Y = 6 cm, the gradient modulus reaches 0.9, the highest value in the entire area. This point is then selected as the starting seed point for tracking. From this seed point, two primary tracking paths are initiated: one along the direction of the gradient vector at that point, for example, toward X = 5.1 cm, Y = 6.05 cm, and the other perpendicular to the gradient vector, for example, toward X = 4.95 cm, Y = 6.1 cm. Using the fourth-order Runge-Kutta numerical method, the computer gradually extends these two initial directions and their opposite directions in small steps, plotting gradient streamlines and equipotential lines. A background gradient level (assuming it's 0.15) is pre-set based on the signal statistics of the entire composite material being tested, such as the average or percentile gradient modulus values ​​across all regions. When the gradient modulus of a tracing gradient streamline or equipotential line falls below the background level of 0.15 at its current location, tracing of that line is terminated. For example, if a gradient streamline originating from a seed point gradually decreases from an initial 0.9 unit to 0.14 at its tip after extending 3 cm, the trace of that streamline is stopped. Tracing in multiple directions from the seed point yields a series of line segments. The outermost edges of all closed curves traced completely, or line segments interrupted due to gradient modulus values ​​falling below the background value or reaching the initially defined edge of the suspect area, are connected to form an envelope. This envelope forms the boundary of the suspect area, which may no longer be the original 10 cm × 10 cm square but a more irregular shape that more closely resembles the actual defect shape, such as an ellipse with a major axis of 7 cm and a minor axis of 4 cm.

[0035] In an optional embodiment, the depth distribution inversion calculation of the suspicious area is performed based on the energy attenuation gradient of the multi-scale energy flow diagram and the depth-direction attenuation of the energy of each fundamental wave mode, specifically: Extract the energy value sequence corresponding to different center frequencies in the multi-scale time-frequency energy flow graph cluster for each point in each suspicious area, and calculate the attenuation gradient of the energy value with frequency; Establishing a depth-attenuation lookup table for each fundamental wave mode when there are no defects in the measured object material, the lookup table containing theoretical attenuation characteristics of wave mode energy as propagation depth increases at different frequencies; constructing a cost function based on the weighted mean square error between a theoretical energy attenuation gradient calculated from the theoretical attenuation characteristics and the actually observed multi-scale energy attenuation gradient; The conjugate gradient method is used to minimize the cost function to obtain the optimal depth value of the point; the above calculation is performed on all points in the suspicious area to obtain the depth distribution.

[0036] For any point within the suspicious region, let's assume it's detection point A. Through the previous steps, a multi-scale time-frequency energy flow diagram cluster has been constructed for this detection point A. This means we have information about the energy at different center frequencies. For example, for the symmetric S0 mode of a Lamb wave, the energy value is 0.9 at a center frequency of 1 MHz, 0.7 at 1.5 MHz, 0.4 at 2 MHz, and 0.2 at 2.5 MHz. This energy value sequence of 0.9, 0.7, 0.4, and 0.2 indicates that the energy of the S0 mode at this point decays as the frequency increases. Based on this energy value sequence over frequency, the rate of energy decay is calculated, which is the attenuation gradient of the energy value over frequency. For example, from 1 MHz to 2.5 MHz, the energy decreases by 0.7. This overall downward trend and the rate of change within each frequency range constitute information about the attenuation gradient. This calculation is repeated for each point within the suspicious region, resulting in an energy attenuation gradient distribution map for the entire region.

[0037] A depth attenuation lookup table for the material in a defect-free state is pre-established. The lookup table records the theoretical attenuation of the energy of Lamb waves, such as the S0 mode, at different frequencies as the depth of propagation in the material increases. For example, a 1.5 MHz S0 wave will attenuate 5% of its energy for every 1 mm of depth it propagates in a sound aluminum plate, while a 2 MHz S0 wave will attenuate 8% of its energy for every 1 mm of depth it propagates. A cost function is constructed to measure the difference between the energy attenuation gradient actually observed at detection point A and the energy attenuation gradient that should theoretically be generated if the defect is at a certain depth. For example, assuming the defect is at a depth of 2 mm, the theoretical attenuation gradient at this time is calculated according to the lookup table, and then compared with the actually observed attenuation gradient to obtain a difference; then, assuming the depth is 2.1 mm, the theoretical attenuation gradient is recalculated and compared again to obtain another difference. In one embodiment, the cost function is the weighted square sum of these differences. Numerical optimization algorithms, such as the conjugate gradient method, are used to continuously adjust the assumed depth values. The goal is to find the optimal depth value that minimizes the cost function. In other words, at this depth, the theoretically calculated attenuation characteristics best match the observed attenuation characteristics. This process is repeated for all points within the suspicious area, and the optimal depth values ​​for each point are obtained, thus depicting the depth distribution of the entire suspicious area.

[0038] In an optional embodiment, the energy proportion and spatial distribution of each fundamental wave mode are used as multi-channel input to perform regularized tomography to reconstruct the three-dimensional geometric shape of the suspicious area, specifically: The energy proportion of the N fundamental wave modes decomposed at each sensor collection point is used as the N-dimensional feature vector of the spatial position; The N-dimensional feature vectors of all the collected points in the suspicious area and the adjacent areas are used to form N independent two-dimensional data channels; An algebraic reconstruction technique is used as the basic tomography algorithm, and a total variation regularization term is utilized to preserve edges and suppress noise. The tomographic imaging voxel grid size is set to 1 / 5 of the minimum spacing of the sensor array. The following optimization problem is solved iteratively: minimizing the data fidelity term plus the regularization coefficient multiplied by the total variation regularization term, where the data fidelity term is the L2 norm difference between the theoretical value and the actual observed value of the two-dimensional data of each channel calculated from the reconstructed three-dimensional voxels through a preset physical acoustic model, and λ is the regularization coefficient. Iteration is performed until the change in the reconstructed voxel model between two consecutive iterations is less than a preset convergence value, thereby obtaining the three-dimensional geometric shape of the suspicious area.

[0039] Specifically, for example, through mode decomposition, three primary guided wave modes are identified at each sensor acquisition point: the symmetric mode S0, the antisymmetric mode A0, and the shear-horizontal mode SH0, where N is 3. At each sensor location on the material surface, such as the point with coordinates x1 and y1, the percentage of each of these three modes' energy relative to the total wave energy is calculated. Specifically, at that point, the S0 mode might account for 50% of the energy, the A0 mode for 35%, and the SH0 mode for 15%. The resulting three-dimensional eigenvector is [0.5, 0.35, 0.15], representing the energy distribution characteristics of the wave modes at that point in space. If 100 sensors are deployed in the inspection area encompassing the suspect area and its surroundings, 100 such three-dimensional eigenvectors will be generated. These data are organized into three independent two-dimensional data channels: the first channel is a two-dimensional plot of the S0 mode energy contribution at all sensor locations; the second channel is a two-dimensional plot of the A0 mode energy contribution; and the third channel is a two-dimensional plot of the SH0 mode energy contribution.

[0040] These three 2D data channels—the spatial energy distribution maps of the S0, A0, and SH0 modes—serve as multi-channel inputs to drive the subsequent 3D tomography process, aiming to reconstruct the 3D structure of the defect within the suspect region. To this end, the 3D space beneath the suspect region is divided into numerous tiny cubic units, or voxels. For example, if the minimum spacing of the sensor array is 10 mm, the side length of a voxel can be set to 2 mm. Algebraic reconstruction techniques are employed as the underlying algorithm, adjusting the acoustic properties of each voxel, such as the speed of sound or attenuation coefficient, through iterative optimization. The optimization goal is to ensure that the theoretical 2D data for each wave mode in the sensor plane, calculated based on the current 3D voxel model and the pre-defined physical acoustic model, is as consistent as possible with the three observed 2D data channels. In one embodiment, this consistency is measured using the L2 norm difference, a data fidelity term. Furthermore, a total variation regularization term is introduced to ensure that the reconstructed defect morphology has sharp edges and a smooth interior, avoiding excessive artifacts. By minimizing the sum of the data fidelity term and the regularization coefficient multiplied by the total variation regularization term, the algorithm continuously adjusts the voxel properties until the change in the three-dimensional voxel model calculated after two consecutive iterations is less than a preset small convergence threshold; the calculation formula for minimizing the sum of the data fidelity term and the regularization coefficient multiplied by the total variation regularization term is , where U represents the three-dimensional geometric shape of the area to be reconstructed, and P(U) is the theoretical value of the two-dimensional data of each channel calculated based on the currently reconstructed three-dimensional voxel model U. is the two-dimensional data of each channel actually observed, is the regularization term, and TV(U) is the total variation regularization term. This results in a reconstruction that reflects the true 3D geometry of the internal defects in the suspected area. For example, if there is an internal delamination defect, the S0 mode may be insensitive to it and transmit, while the A0 mode may be strongly scattered or undergo energy conversion on its surface. These different responses are captured by their respective energy channels and jointly inverted through multi-channel tomography to obtain the size, shape, and spatial posture of the defect.

[0041] It should be noted that the terms "comprises," "includes," or any other variations thereof are intended to encompass non-exclusive inclusion, such that a process, method, commodity, or apparatus that includes a series of elements includes not only those elements but also other elements not explicitly listed, or includes elements inherent to such process, method, commodity, or apparatus. In the absence of further limitations, an element defined by the phrase "comprises a..." does not exclude the presence of other identical elements in the process, method, commodity, or apparatus that includes the element.

[0042] This specification may be described in the general context of computer-executable instructions, such as program modules, executed by a computer. Generally, program modules include routines, programs, objects, components, data structures, and the like that perform specific tasks or implement specific abstract data types. This specification may also be practiced in distributed computing environments where tasks are performed by remote processing devices connected through a communications network. In a distributed computing environment, program modules may be located in both local and remote computer storage media, including storage devices.

[0043] It will be apparent to those skilled in the art that various modifications and variations may be made to one or more embodiments of this specification. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of one or more embodiments of this specification shall be included within the scope of the claims of this specification.

Claims

1. A non-destructive testing method, characterized in that: The method comprises the following steps: An excitation source generates elastic waves in a predetermined area of ​​the measured object, and a sensor array collects response signals; the response signals are subjected to dynamic coupled wave mode decomposition to obtain fundamental wave mode components and energy proportions; and a multi-scale time-frequency energy flow graph cluster based on the fundamental wave mode components and energy proportions is constructed for multiple discrete time-frequency scale units; Calculating the directional gradient of each energy flow graph in the energy flow graph cluster with the energy amplitude as the weight to obtain a gradient vector field; determining a suspicious area based on the topological structural characteristics of the gradient vector field and the abnormal threshold of the wave mode energy ratio in the energy flow graph cluster; and determining the boundary of the suspicious area based on the trajectory of the gradient vector field; The depth distribution inversion calculation of the suspicious area is performed based on the energy attenuation gradient of the multi-scale energy flow diagram and the depth-direction attenuation of the energy of each fundamental wave mode; and the energy proportion and spatial distribution of each fundamental wave mode are used as multi-channel input to perform regularized tomography and reconstruct the three-dimensional geometric shape of the suspicious area.

2. The method according to claim 1, characterized in that The dynamic coupled wave mode decomposition of the response signal is performed to obtain the fundamental wave mode component and energy proportion, specifically: The response signal collected by each sensor is converted to the frequency-wavenumber domain using a two-dimensional fast Fourier transform. In the frequency-wavenumber domain, the isolation band width is determined based on the theoretical dispersion curves of different wave modes. The fundamental wave mode components are separated using multiple predefined two-dimensional bandpass filter banks that match the theoretical dispersion curve shapes. The amplitude spectrum of the fundamental wave mode component is integrated in the frequency-wavenumber domain to obtain the total energy of the fundamental wave mode component, and the proportion of the fundamental wave mode component in the total energy of all fundamental wave modes at the current acquisition point is calculated to obtain the energy proportion.

3. The method according to claim 1, characterized in that The determining of the suspicious area based on the topological structural characteristics of the gradient vector field and the abnormal threshold of the wave mode energy ratio in the energy flow graph cluster is specifically as follows: Calculate the divergence field and curl field of each weighted directional gradient field; The area where the value of the divergence field is lower than the negative threshold is regarded as the energy convergence area, and the area where the value is higher than the positive threshold is regarded as the energy divergence area; The area in the curl field with a value higher than the amplitude threshold is regarded as the energy rotation area; The spatial union of the energy convergence region, the energy divergence region and the energy rotation region is regarded as the topological singular region; Calculating the normalized deviation of the energy proportion of each fundamental wave mode at each spatial position relative to the mean energy proportion of the mode in the entire measured area in the energy flow graph cluster; When the energy ratio normalized deviations of at least two fundamental wave modes are greater than a preset value, the spatial position is marked as an energy ratio abnormal point; The spatial intersection of the topological singular region and all energy proportion abnormal point sets is taken as the suspicious region.

4. The method according to claim 1, wherein The determining of the boundary of the suspicious area based on the trajectory of the gradient vector field is specifically as follows: Inside each suspicious area, the point with the largest modulus of the gradient vector field is selected as the seed point; Starting from the seed point, the fourth-order Runge-Kutta method is used to bidirectionally trace gradient streamlines and equipotential lines along the direction of the gradient vector field and the direction perpendicular to the gradient vector field until the gradient modulus of the streamlines or equipotential lines is lower than the background gradient level predetermined based on the statistical characteristics of the signal of the entire object under test; the outermost envelope of all closed gradient streamlines and equipotential line segments or those truncated at the edge of the suspicious area is used as the boundary of the suspicious area.

5. The method according to claim 1, characterized in that The depth distribution inversion calculation of the suspicious area is performed based on the energy attenuation gradient of the multi-scale energy flow diagram and the depth direction attenuation of the energy of each fundamental wave mode, specifically: Extract the energy value sequence corresponding to different center frequencies in the multi-scale time-frequency energy flow graph cluster for each point in each suspicious area, and calculate the attenuation gradient of the energy value with frequency; Establishing a depth-attenuation lookup table for each fundamental wave mode when there are no defects in the measured object material, the lookup table containing theoretical attenuation characteristics of wave mode energy as propagation depth increases at different frequencies; constructing a cost function based on the weighted mean square error between a theoretical energy attenuation gradient calculated from the theoretical attenuation characteristics and the actually observed multi-scale energy attenuation gradient; The conjugate gradient method is used to minimize the cost function to obtain the optimal depth value of the point; the above calculation is performed on all points in the suspicious area to obtain the depth distribution.

6. The method according to claim 1, characterized in that The energy proportion and spatial distribution of each fundamental wave mode are used as multi-channel input to perform regularized tomography to reconstruct the three-dimensional geometric shape of the suspicious area, specifically: The energy proportion of the N fundamental wave modes decomposed at each sensor collection point is used as the N-dimensional feature vector of the spatial position; The N-dimensional feature vectors of all the collected points in the suspicious area and the adjacent areas are used to form N independent two-dimensional data channels; An algebraic reconstruction technique is used as the basic tomography algorithm, and a total variation regularization term is utilized to preserve edges and suppress noise. The tomographic imaging voxel grid size is set to 1 / 5 of the minimum spacing of the sensor array. The following optimization problem is solved iteratively: minimizing the data fidelity term plus the regularization coefficient multiplied by the total variation regularization term, where the data fidelity term is the L2 norm difference between the theoretical value and the actual observed value of the two-dimensional data of each channel calculated from the reconstructed three-dimensional voxels through a preset physical acoustic model, and λ is the regularization coefficient. Iteration is performed until the change in the reconstructed voxel model between two consecutive iterations is less than a preset convergence value, thereby obtaining the three-dimensional geometric shape of the suspicious area.

7. A non-destructive testing system, characterized in that: The system includes the following modules: An acquisition and decomposition module is configured to generate elastic waves in a predetermined area of ​​the measured object by an excitation source, and to collect response signals from a sensor array; perform dynamic coupled wave mode decomposition on the response signals to obtain fundamental wave mode components and energy proportions; and construct a multi-scale time-frequency energy flow graph cluster based on the fundamental wave mode components and energy proportions for multiple discrete time-frequency scale units; a preliminary identification module for calculating the directional gradient of each energy flow graph in the energy flow graph cluster with the energy amplitude as the weight, and obtaining a gradient vector field; Determining a suspicious area based on the topological structural characteristics of the gradient vector field and the abnormal threshold of the wave mode energy ratio in the energy flow graph cluster; and determining the boundary of the suspicious area based on the trajectory of the gradient vector field; The morphology extraction module is used to perform depth distribution inversion calculation of the suspicious area based on the energy attenuation gradient of the multi-scale energy flow diagram and the depth-direction attenuation of the energy of each fundamental wave mode; and use the energy proportion and spatial distribution of each fundamental wave mode as multi-channel input to perform regularized tomography and reconstruct the three-dimensional geometric morphology of the suspicious area.

8. The system according to claim 7, characterized in that The dynamic coupled wave mode decomposition of the response signal is performed to obtain the fundamental wave mode component and energy proportion, specifically: The response signal collected by each sensor is converted to the frequency-wavenumber domain using a two-dimensional fast Fourier transform. In the frequency-wavenumber domain, the isolation band width is determined based on the theoretical dispersion curves of different wave modes. The fundamental wave mode components are separated using multiple predefined two-dimensional bandpass filter banks that match the theoretical dispersion curve shapes. The amplitude spectrum of the fundamental wave mode component is integrated in the frequency-wavenumber domain to obtain the total energy of the fundamental wave mode component, and the proportion of the fundamental wave mode component in the total energy of all fundamental wave modes at the current acquisition point is calculated to obtain the energy proportion.

9. The system according to claim 7, wherein: The determining of the suspicious area based on the topological structural characteristics of the gradient vector field and the abnormal threshold of the wave mode energy ratio in the energy flow graph cluster is specifically as follows: Calculate the divergence field and curl field of each weighted directional gradient field; The area where the value of the divergence field is lower than the negative threshold is regarded as the energy convergence area, and the area where the value is higher than the positive threshold is regarded as the energy divergence area; The area in the curl field with a value higher than the amplitude threshold is regarded as the energy rotation area; The spatial union of the energy convergence region, the energy divergence region and the energy rotation region is regarded as the topological singular region; Calculating the normalized deviation of the energy proportion of each fundamental wave mode at each spatial position relative to the mean energy proportion of the mode in the entire measured area in the energy flow graph cluster; When the energy ratio normalized deviations of at least two fundamental wave modes are greater than a preset value, the spatial position is marked as an energy ratio abnormal point; The spatial intersection of the topological singular region and all energy proportion abnormal point sets is taken as the suspicious region.

10. The system according to claim 7, wherein: The determining of the boundary of the suspicious area based on the trajectory of the gradient vector field is specifically as follows: Inside each suspicious area, the point with the largest modulus of the gradient vector field is selected as the seed point; Starting from the seed point, the fourth-order Runge-Kutta method is used to bidirectionally trace gradient streamlines and equipotential lines along the direction of the gradient vector field and the direction perpendicular to the gradient vector field until the gradient modulus of the streamlines or equipotential lines is lower than the background gradient level predetermined based on the statistical characteristics of the signal of the entire object under test; the outermost envelope of all closed gradient streamlines and equipotential line segments or those truncated at the edge of the suspicious area is used as the boundary of the suspicious area.

Images