Method for quantitatively evaluating three-dimensional morphology of thinning defect in plate under ultrasonic guided wave field undersampling scanning

Abstract

The invention discloses a method for quantitatively evaluating the three-dimensional shape of a thinning defect in a plate under ultrasonic guided wave field undersampling scanning, belongs to the field of industrial nondestructive testing, and particularly relates to high-precision three-dimensional shape reconstruction and thickness loss quantification of a corrosion defect in a metal plate shell structure (such as aerospace skin, a pressure vessel wall, an oil pipeline and the like). According to the method, sparse sampling guided wave field reconstruction, full waveform inversion (FWI) and priori knowledge constraint are fused, Shannon sampling theorem limitation is broken through, efficient measurement of corrosion depth, contour and distribution under the conditions of a single excitation source and a limited view angle is achieved, and the method is suitable for rapid online detection in a complex industrial environment.

Description

Technical Field

[0001] This invention belongs to the field of industrial nondestructive testing, specifically involving high-precision three-dimensional morphological reconstruction and thickness loss quantification of corrosion defects in metal plate and shell structures (such as aerospace skins, pressure vessel walls, and oil pipelines). By integrating sparse sampling guided wave field reconstruction, full waveform inversion (FWI), and prior knowledge constraints, it overcomes the limitations of Shannon's sampling theorem, achieving efficient measurement of corrosion depth, profile, and distribution under single excitation source and limited viewing angle conditions. It is suitable for rapid online detection in complex industrial environments. Background Technology

[0002] In the field of nondestructive testing, ultrasonic guided waves are widely used for corrosion detection of plate and shell structures due to their advantages such as long-distance propagation and high sensitivity. However, existing technologies have the following bottlenecks:

[0003] 1. Traditional guided wave field imaging technology relies on dense spatial sampling (e.g., the minimum wavelength of A0 mode guided wave at 300kHz is 6.1mm, and the sampling interval needs to be ≤3.05mm), resulting in low data acquisition efficiency. Under sparse sampling, Fourier transform fails due to wavenumber aliasing and cannot separate multimode guided wave components (reference: EB Flynn, SY Chong, GJ Jarmer, J.-R. Lee, Structural imaging through local wavenumber estimation of guided waves, NDT&E International. 59(2013) 1–10. https: / / doi.org / 10.1016 / j.ndteint.2013.04.003.).

[0004] 2. Compressed sensing methods (such as the sparse reconstruction algorithm proposed by Mesnil et al.) can reduce the amount of data through compressed sampling, but their core problem is that they can only locate defects and cannot quantify depth. In addition, these methods require the pre-construction of wavenumber-thickness dictionaries, which are poorly adapted to non-uniform corrosion or complex contour defects (Reference O.Mesnil, M.Ruzzene, Sparse wavefield reconstruction and source detection using compressed sensing, Ultrasonics. 67(2016)94–104. https: / / doi.org / 10.1016 / j.ultras.2015.12.014.).

[0005] 3. Full waveform inversion (FWI) technology approximates measured data by iteratively optimizing model parameters, which can theoretically achieve high-precision quantification. However, traditional FWI requires multi-angle excitation data to cover the entire defect perspective, and the lack of information under a single excitation source leads to ill-conditioned inversion (J. Rao, M. Ratassepp, Z. Fan, Guided wave tomography based on full waveform inversion, IEEE Transactions Ultrasonics, Ferroelectrics, Frequency Control. 63(2016) 737–745. https: / / doi.org / 10.1109 / TUFFC.2016.2536144.). Meanwhile, existing acoustic models neglect guided wave mode conversion effects (such as S0 and A0 mode coupling), resulting in significant depth quantization errors (reference AAEZimmermann, P. Huthwaite, B. Pavlakovic, High-resolution thickness maps of corrosion using SH1 guided wave tomography, Proceedings Royal Society A: Mathematical, Physical Engineering Sciences. 477(2021)20200380. https: / / doi.org / 10.1098 / rspa.2020.0380.).

[0006] The physical limitations of sparse sampling and a single excitation source result in incomplete waveguide field information; furthermore, existing algorithmic frameworks fail to effectively integrate the localization and quantization processes, leading to unstable inversion results. Therefore, a method that can collaboratively optimize sparse reconstruction and the physical model is urgently needed to overcome the limitations of traditional techniques. How can the following objectives be achieved under conditions of a single excitation source and sparse sampling (sampling interval > half a wavelength)?

[0007] 1. High-precision reconstruction of the guided wavefield and extraction of defect scattering characteristics;

[0008] 2. Solve the inverse problem of FWI stably to avoid local optima; Summary of the Invention

[0009] 1. Technical problems to be solved

[0010] Sparse sampling guided wavefield reconstruction: under the condition that the sampling interval is more than half a wavelength, wavenumber aliasing is eliminated and the spatiotemporal signal of the guided wavefield is restored;

[0011] Simultaneous quantization of defect contour and depth: solving the problem of coordinated optimization between compressed sensing localization and FWI depth inversion;

[0012] Single-source FWI stability: By using prior knowledge constraints, noise amplification and divergence during the inversion process are suppressed.

[0013] 2. Core technology innovation

[0014] This invention solves the problem of quantifying thinning defects in thin-walled structures under sparse sampling and single excitation source conditions through two-stage collaborative optimization:

[0015] Phase 1: Defect Localization and Contour Extraction

[0016] Based on the straight-ray model, a waveguide field equation is constructed, representing the measured signal as a linear combination of sparse scattering sources. A sparse optimization problem is solved using basis pursuit denoising (BPDN) to extract the distribution of defect scattering sources. The Quickhull algorithm is used to calculate the convex hull boundary of the scattering source point set, and a characteristic function of the defect region is defined by combining the number of turns integral, providing prior knowledge for stage 2.

[0017] Phase 2: Deep Quantization Based on FWI

[0018] The two-dimensional frequency-domain acoustic equations are discretized using the finite difference method to construct the FWI objective function. The characteristic function of stage 1 is incorporated into the parameter update process to suppress noise in non-defect regions and improve inversion stability. The inverted wave velocity is mapped to the plate thickness distribution through the Rayleigh-Lamb dispersion equation, achieving physical quantification of corrosion depth.

[0019] The core of this invention lies in two-stage collaboration. The localization result of stage 1 provides spatial constraints for stage 2. Through the feature function defined by the convex hull boundary, the FWI parameter update is restricted to the defect region, avoiding noise interference in global optimization. The flowchart of this scheme is as follows: Figure 2 As shown.

[0020] Overall architecture:

[0021] Phase 1: Defect Localization and Contour Extraction

[0022] Phase 1 comprises three steps: First, a sparse model of the ultrasonic guided wave field is constructed by building a ray model; then, basis tracing noise reduction is used to solve for the location and intensity of the wave field scattering source; finally, the defect region characteristic function is defined using the location of the scattering source to complete the solution representation of the defect contour. The technical details of each step are as follows.

[0023] (1) Sparse modeling of ultrasonic guided wave field:

[0024] First, a ray propagation model of the ultrasonic guided wave field is constructed, with subscripts m and s indicating the measurement point and scattering source, respectively. The location of the measurement point is known. Figure 1 The yellow circle in the middle represents the measurement point, and the location of the scattering source is unknown. Figure 1 The gray circle represents the hypothetical location of the scattering source, and (x, y) represents the position of the wave field in a Cartesian coordinate system located on the plane of the plate. Based on the straight-ray model, the measurement point position (x, y) is defined. m ,y m ) and the assumed scattering source location (x) s ,y s Euclidean distance between Calculate the frequency domain guided wavefield equations:

[0025]

[0026] Symbol definition:

[0027] ω: angular frequency;

[0028] i: represents the imaginary unit

[0029] p m (ω): The sound pressure intensity at the measurement point m at angular frequency ω, which is the Fourier transform of the measurement signal in the time domain;

[0030] μ: Ultrasonic guided wave modes, including: antisymmetric mode, symmetric mode and horizontal shear mode;

[0031] N s The total number of assumed scattering sources can be any integer value greater than the number of measurement points;

[0032] The equivalent complex amplitude of the assumed scattering source at position s at angular frequency ω and mode μ is a complex function of the variable ω in the frequency domain;

[0033] The wave number of a wave with mode μ propagating along a straight path between the measurement point m and the scattering source s at an angular frequency ω can be obtained by solving the ultrasonic guided wave dispersion curve.

[0034] In the above modeling, the measurement points can be set to be uniformly distributed in the detection area, with a spacing greater than half the wavelength. The possible locations of the assumed scattering sources can be assumed to be randomly distributed in the detection area.

[0035] (2) Basis Pursuit Denoising (BPDN) Solution:

[0036] Let N m This represents the total number of measurement points. The total number of scattering sources is N. s The total number of measurement points N m Under the given conditions, N can be obtained. m An equation of the same form as (1). This N mThe equations are integrated and represented as matrix equations:

[0037]

[0038] Given p and H, the problem of identifying the ultrasonic guided wave scattering source q caused by defects in the plate is equivalent to solving the inverse problem of (2), which can be expressed as the following sparse optimization problem:

[0039]

[0040] Symbol definition:

[0041] ‖·‖1: L1 norm (the sum of the absolute values ​​of the vector elements);

[0042] ‖·‖2: L2 norm (the square root of the sum of the squares of the vector elements);

[0043] p: Measurement point vector Dimension N m ×1, corresponding to p at each measurement point in equation (1) m A vector composed of (ω), where each component represents the sound pressure intensity value at each measurement point at angular frequency ω, can be obtained by acquiring time-domain waveform signals and then performing Fourier transform.

[0044] q: Scattering source vector Dimension N s ×1, corresponding to different positions in equation (1) The vector formed by these variables is the quantity to be determined and is usually sparse.

[0045] σ: Noise margin, taken as 5% of the measured signal energy;

[0046] H: Ultrasonic guided wave ray propagation matrix with mode μ and dimension N. m ×N s The corresponding elements are as follows, and the physical meaning of each element in the formula is the same as that in (1):

[0047]

[0048] The sparse optimization problem in equation (3) above can be solved using the SPGL1 algorithm.

[0049] Once the actual ultrasonic guided wave scattering source is determined, its location can be obtained from the pre-defined relationship between the assumed scattering source coordinates and the vector q. Since the ultrasonic guided wave scattering source originates from the influence of thinning defects on ultrasonic guided wave propagation, the scattering source is distributed within the boundary or region of the corrosion defect. Therefore, the location and outline of the corrosion defect can be identified based on the location of the actual scattering source; the geometric outline of the corrosion region will be further characterized in the next step.

[0050] (3) Defect profile solution representation:

[0051] First, a scattering source screening process is performed, retaining hypothetical scattering sources with a certain intensity and removing hypothetical scattering sources with weak scattering intensity. A threshold q is set. threshold = 0.1 × max(q), where max(q) represents taking the maximum value of the elements in vector q, and the threshold intensity of the scattering source is taken as 0.1 times the maximum value of the elements in vector q. Using the set threshold q... threshold From the scattering source vector q obtained in the previous step, scattering sources with amplitudes above a preset threshold are selected to form a point set of effective scattering source locations.

[0052] Then, convex hull calculations are used to model the distribution of scattering sources; the location where scattering sources are concentrated is the defect region. The Quickhull algorithm is used to generate the minimum convex hull, denoted as C. s S is the set of scattering source points surrounding the wave field plane. s For all points, the shortest convex closed curve is used. The loop integral is defined using the non-zero winding rule to determine whether the scattering source location is within the defect region. For each point (x, y) located in the Cartesian coordinate system where the wave field is located, the closed curve C is calculated in the polar coordinate system with itself as the origin. s The number of wraps around this point is calculated using the following formula:

[0053]

[0054] Finally, for the defective region, the feature function m(x,y) is defined as follows:

[0055]

[0056] In the setting of the characteristic function m(x,y), if The corresponding point (x, y) is located within the defect area; all other areas are outside the defect. For the area inside the defect, m(x, y) is set to a value of 1. For the area outside the defect, m(x, y) is set to a positive real number close to 0, 0.01.

[0057] Symbol definition:

[0058] ∮: Represents the loop integral;

[0059] θ: Convex hull C s Loop integral variables on;

[0060] Convex hull C s The number of turns around the point (x,y);

[0061] In the next stage, we will use the feature function m(x,y) as prior knowledge and substitute it into equation (12) to improve the performance of thickness loss quantization.

[0062] Phase 2: Deep Quantization Based on Full Waveform Inversion (FWI)

[0063] Phase 2 comprises four steps. First, based on the frequency domain two-dimensional acoustic wave equation, a discrete difference model of wave field propagation is constructed. Next, an FWI objective function is constructed and the parameter update gradient is solved. Then, using the characteristic function obtained in Phase 1, the parameters are updated using mask constraint optimization, iteratively solving for the wave field velocity distribution. Finally, using the dispersion relationship of guided waves in the plate, the wave field velocity distribution is converted into a plate thickness distribution, thereby achieving a quantitative assessment of the three-dimensional morphology of defects in the plate. The technical details of each step are as follows.

[0064] (1) Acoustic wave equation modeling:

[0065] In the aforementioned defined wavefield plane rectangular coordinate system, the frequency domain two-dimensional acoustic wave equation is as follows:

[0066]

[0067] Where v(x,y,ω) represents the wave velocity at position (x,y) in the wave field at angular frequency ω; p(x,y,ω) is the sound pressure at position (x,y) in the wave field at angular frequency ω, which is the frequency domain representation of the disturbance signal generated by the sound source, and is the Fourier transform of the time domain signal at that position; q(x,y,ω) is the sound source at position (x,y) in the wave field, which is the frequency domain representation of the input function driving the sound wave generation, and is a complex function of angular frequency ω.

[0068] The wavefield is divided into a Cartesian coordinate system, and the discrete, densely packed grid points are denoted as {(x... i ,y j )|i∈{1,…,N x},j∈{1,…,N y}}, where N x and N y Let h be the number of discrete grid points in the x and y directions, respectively. x Let x be the discrete dimension in the x-direction, that is: i+1 -x i =h x (i∈{1,…,N x -1}), let h be... y The discrete dimension in the y-direction: y j+1 -y j =h y (j∈{1,…,N y -1}).

[0069] Discretizing equation (7) using the finite difference space yields:

[0070]

[0071] For any given angular frequency, equation (8) can be expressed in matrix form as follows:

[0072]

[0073] Symbol definition:

[0074] ⊙: Hadamard product (element-by-element multiplication);

[0075] Y: The model parameter matrix, composed of the reciprocal of the square of the wave velocity at each location in the wave field, i.e., Y = (1 / v 2 (x i ,y j ,ω),i∈{1,…,N x},j∈{1,…,N y}, where N x and N y These are the number of discrete grid points in the x and y directions, respectively;

[0076] P: Sound pressure field matrix, composed of sound pressure at each location in the wave field, i.e., P = (p(x) / (x)) i ,y j ,ω),i∈{1,…,N x},j∈{1,…,N y} represents the grid (x) i ,y j The sound pressure at angular frequency ω is obtained by the Fourier transform of the amplitude of the time-domain measured waveform at the corresponding measurement point.

[0077] Q: Source term matrix, composed of the positions of sound sources in the wave field, i.e., Q = (q(x)). i ,y j ,ω),i∈{1,…,N x},j∈{1,…,N y}, where the non-zero elements correspond to known wavefield sound sources, and are complex functions of angular frequency ω, (x i ,y j ) indicates its location.

[0078] The second-order difference matrix corresponds to the discretization of the x and y directions on the wave field plane, respectively. It is directly extracted from equation (8). The dimension corresponds to the grid division, and the element value corresponds to the coefficient of equation (8). The matrix has a dimension of N x ×N x , The matrix has a dimension of N y ×Ny The specific form is as follows:

[0079] express Transpose transformation;

[0080] (2) FWI objective function and gradient calculation:

[0081] Based on the acoustic wave equation (9), the node settings in the discrete grid correspond to the wave field measurement point settings, forming a definite second-order difference matrix. and The initial value of the model parameter matrix Y is set to be an identity matrix, and the sound source term Q is set to correspond to the actual ultrasonic guided wave field detection excitation signal. Therefore, by solving the matrix equation (9), the simulation result of the wave field state, i.e., the simulated wave field sound pressure data matrix, denoted as P, can be obtained. sim Furthermore, the simulated wavefield is continuously approximated to the actual wavefield through the following steps, and the model parameter matrix Y continuously approximates the actual situation. The approximation process of the ultrasonic guided wavefield can be equivalently transformed into an optimization problem of solving the objective function at each frequency, which is a typical FWI problem. First, the objective function J is defined as:

[0082]

[0083] Symbol definition:

[0084] ‖·‖2: L2 norm (the square root of the sum of the squares of the vector elements);

[0085] e: Residual vector, representing the difference between the measured and analog signals;

[0086] P meas : Wavefield acoustic pressure measurement data matrix;

[0087] P sim The intermediate quantity obtained by solving (9) is the simulated wave field sound pressure data matrix.

[0088] R: Spatial sampling matrix, which extracts data corresponding to the measurement points from the simulated wave field sound pressure;

[0089] Then solve the optimization problem: The gradient g for parameter update is obtained using the adjoint state method:

[0090]

[0091] Symbol definition:

[0092] p adj The accompanying state variables are solved using the conjugate gradient method.

[0093] diag(p): a diagonal matrix with diagonal elements in the form of a sound pressure matrix p vector;

[0094] diag(p) H : Represents the conjugate transpose of diag(p);

[0095] Re: The operation of taking the real part;

[0096] (3) Mask constraint optimization parameter update:

[0097] Since the wave velocity change caused by plate thickness loss only occurs in the defect region and can be represented by the prior knowledge obtained in Stage 1, a mask-constrained optimization strategy is adopted for parameter updates. The wave field and wave velocity distribution can be solved by iterative iteration. Combining the characteristic function obtained in Stage 1, the gradient update range is constrained as follows:

[0098] vec(Y new ) = vec(Y old )-η·vec(M)⊙g (12)

[0099] Symbol definition:

[0100] Y old The parameter matrix Y in formula (9) is composed of the reciprocals of the squares of the wave velocities at each position in the wave field. As the gradient g is updated, it iterates to the old parameter values ​​before the update.

[0101] Y new As the gradient g is updated, the iteration proceeds to the updated parameter values.

[0102] η: Step size parameter, controls the update speed of the parameter, and takes the value of a real number in (0,1);

[0103] M: Mask matrix, used to refine the update gradient g of the model parameters. It is constructed based on the values ​​of the feature function obtained from formula (6) in stage 1 at each position in the wave field, i.e., M = (m(x i ,y j ),i∈{1,…,N x},j∈{1,…,N y}, where N x and N y These represent the number of discrete points in the x and y directions of the discrete grid points in the wavefield plane rectangular coordinate system, respectively;

[0104] vec(M): Expands matrix M into a vector column by column.

[0105] Corrosion quantification can be achieved by transforming the model parameters. The transformed model parameters Y are converted into phase velocities, and then converted into plate thickness distributions based on the dispersion relationship of the waveguide phase velocities. This dispersion relationship can be calculated by (13) below.

[0106] (4) Wave velocity-thickness conversion:

[0107] By solving the Rayleigh-Lamb dispersion equations, the dispersion relationship between wave velocity and thickness for the wave field in a plate is obtained. The wave velocity distribution is then transformed to reflect the plate thickness distribution. The Rayleigh-Lamb dispersion equations are as follows:

[0108]

[0109] h z The thickness of an isotropic plate;

[0110] k: indicates the wavenumber of the ultrasonic Lamb wave;

[0111] α, β: coefficients, their formulas are:

[0112]

[0113] C L : Represents the velocity of the longitudinal wave, which is an inherent property of the material. Given a specific material, it is a known input quantity, similar to v(x) mentioned earlier. i ,y j ,ω) are irrelevant;

[0114] C T : Represents the velocity of the transverse wave, which is an inherent property of the material. Given a specific material, it is a known input quantity, similar to v(x) mentioned earlier. i ,y j ,ω) are irrelevant.

[0115] Equation (13) provides the relationship between wavenumber and thickness at a given angular frequency ω. Let h... z With a thickness increment of 0.01 mm, numerical calculations yield the wavenumber-dispersion relationship for each discrete thickness within the plate thickness range. Then, through linear interpolation, a complete and continuous wavenumber-dispersion curve containing the discrete thickness values ​​and the regions between them is obtained. Dispersion curves for both A0 and S0 modes are calculated separately, and the combined curves provide three-dimensional wavenumber-frequency-thickness dispersion surface plots for both modes. At each discrete frequency, based on the three-dimensional wavenumber-frequency-thickness dispersion relationship, the plate thickness distribution can be obtained from the wavenumber distribution by looking up a table.

[0116] Comparative experiments with traditional methods show that the present invention has significant advantages in terms of noise tolerance and sensitivity to shallow defects. For example, at the same signal-to-noise ratio (SNR = 20dB), traditional FWI produces a large number of artifacts in the inversion results because it does not introduce spatial constraints; while the present invention effectively suppresses inversion errors in non-defect regions through mask optimization.

[0117] Quantitative accuracy: In a 3mm aluminum plate, the average absolute error of a 2mm deep defect is 0.12mm, with a contour overlap of 92.3%. Traditional FWI requires multi-angle excitation data to cover the entire defect perspective; under a single excitation source, information loss leads to ill-conditioned inversion. This scheme, however, uses defect localization and contour extraction obtained in the first stage to provide prior knowledge for the second-stage FWI solution process, while simultaneously reducing the need to update defect-free areas. This allows for a more effective quantitative assessment of corrosion defects, making the model more stable.

[0118] Efficiency Improvement: The sampling interval has been increased from 3.65 mm to 5 mm, and the single detection time has been reduced from 120 minutes to 75 minutes. By overcoming the limitation of Shannon's sampling theorem on the sampling interval, the number of wavefield sampling points per unit area has been reduced. Undersampling can be used to quantitatively evaluate corrosion defects in spatial waveguide plate structures, thereby improving the efficiency of waveguide field scanning detection technology.

[0119] Robustness: Compared with the traditional FWI method, this method reduces the number of FWI iterations by 40% under a single excitation source (200 times → 120 times), and the iteration speed reaches 10 times / second after GPU acceleration. Since the feature function obtained in stage 1 (6) is used to mask the update of the FWI iterative model in stage 2, the set FWI iteration stopping criterion can be achieved with fewer iterations, that is, the accuracy requirement of the final quantitative evaluation of corrosion defects can be achieved. Figure 5 As shown, the convergence error of this method is significantly better than that of the traditional FWI method. Attached image description:

[0120] Figure 1 Schematic diagram of a guided wavefield scanning system. Left: Layout of the piezoelectric exciter (PZT) and laser scanning probe; Right:

[0121] Signal processing flow (data acquisition → sparse reconstruction → FWI inversion).

[0122] Figure 2 Flowchart of the two-stage method. Stage 1: Sparse scattering source extraction → Convex hull generation → Feature function construction;

[0123] Phase 2: Discretization of acoustic equations → Calculation of adjoint gradients → Update of mask constraints.

[0124] Figure 3 : Snapshot of guided wavefield from finite element simulation (comparison of scattered waves from 0.5mm, 1mm, and 2mm defects). Rows: different defect depths; Columns: wavefield distribution at 20μs, 30μs, and 40μs.

[0125] Figure 4Hardware system block diagram. Excitation module: Signal generator (Agilent 33500B) → Power amplifier (NFHSA4014) → PZT; Receiving module: Laser interferometer (Polytec PSV-500) → Data acquisition card (NI).

[0126] PXIe-5162); Control unit: LabVIEW real-time control of scan path and trigger synchronization.

[0127] Figure 5 Comparison of FWI iteration errors. Blue line: Traditional FWI (without mask constraint); Red line: The method of this invention (with mask constraint).

[0128] Figure 6 : A physical image of an aviation aluminum plate containing rectangular thinning defects.

[0129] Figure 7 : Comparison chart of defect detection output results. The first row of three figures are from the literature (J. Rao, M. Ratassepp, Z.).

[0130] Fan, Guided wave tomography based on fullwaveforminversion, IEEETransactions Ultrasonics, Ferroelectrics, Frequency Control. 63 (2016) 737–745.

[0131] The result of the method in https: / / doi.org / 10.1109 / TUFFC.2016.2536144. is equivalent to ignoring stage 1 and directly using stage 2 in this method; the second line is the output result diagram of this method. Detailed Implementation

[0132] 1. Data Acquisition and Preprocessing

[0133] 1) Instrument hardware configuration and data acquisition:

[0134] Experimental specimen: A 3mm thick aluminum plate with a pre-fabricated rectangular groove measuring 20mm × 30mm and 2mm deep to simulate a clearly defined defect. The measurement area was a 60mm × 70mm rectangular region covering the location of the pre-fabricated defect.

[0135] Excitation source: A single excitation source is located in the lower left corner of the waveguide field scanning area. It is connected to the piezoelectric sensor by a signal generator through a power amplifier for excitation. The excitation source is set to a 5-cycle Hanning window amplitude-modulated signal with a center frequency of 300kHz.

[0136] Sampling interval: The A0 mode wave has a relatively shorter wavelength within the frequency band, with its shortest wavelength being approximately 6.1 mm. The maximum spatial sampling interval satisfying Shannon's sampling theorem is approximately 3.05 mm. The spatial sampling interval for the experimental data was set to 4.5 mm to ensure undersampling of the guided wavefield.

[0137] Data Acquisition: 195 sampling points are arranged in a pre-defined square grid within the guided wavefield scanning area. A laser probe scans each point within the detection area along equidistant vertical lines from left to right. Measurement positions are arranged from bottom to top along the scanning line. Simultaneously, a dual-wave mixing laser interferometer is controlled to acquire the guided wavefield signals at each measurement position, collecting discrete ultrasonic guided wave spatiotemporal signals. The laser interferometer at the receiving end has a sampling rate of 1 GHz and a spatial resolution of 0.1 mm. A servo controller is used to achieve synchronous excitation and reception. The experimental system block diagram is shown below. Figure 4 As shown, the displacement amplitude response signal acquired by the set is Fourier transformed and used as measurement data, which can then be further processed and used for wavefield imaging in the frequency-wavenumber domain.

[0138] 2) Signal preprocessing:

[0139] First, cross-correlation calibration is performed to eliminate signal bias between the elastic model and the acoustic model in FWI, as shown in the following formula:

[0140] △t=argmax τ ∫p ref (t)p meas (t+τ)dt (15)

[0141]

[0142] p ref The amplitude of the simulation reference signal used for matching is the response signal of the flawless plate at each measurement position in the waveguide field calculated based on the acoustic model.

[0143] p meas The measured amplitude of the signal;

[0144] △t: time shift;

[0145] α: Amplitude scaling factor, used to adjust the measured signal to match the reference signal.

[0146] Furthermore, the initial wave envelope is truncated, retaining the initial wave time window (e.g., 0-50 μs) to avoid interference from multiple reflections, and a Fourier transform is performed. A screening criterion based on Frèchet distance is adopted to eliminate measurement points with large waveform differences; the distance screening threshold is set to F. threshold=0.2×max(F), where F is a vector composed of the distance between the measured values ​​of each measurement point on the board and the reference value Frèchet, and max(F) represents taking the maximum value of the elements in vector F. Here, the filtering threshold is set to 0.2 times the maximum distance.

[0147] 2. Algorithm Implementation Process

[0148] 1) Ray matrix construction: corresponding to the sparse modeling of the waveguide field in "Phase 1" of the technical solution, calculate the ray propagation path between each measurement point and the scattering source according to the preset scattering source grid (spacing 1mm) and generate the matrix.

[0149] 2) Sparse optimization solution: Corresponding to the basis pursuit denoising (BPDN) process in "Stage 1", configure the SPGL1 solver and set the noise tolerance σ = 0.05 × 1 / 2 p meas ‖2, maximum number of iterations 1000, output sparse solution q, obtain sparse scattering source distribution.

[0150] 3) Eliminating the influence of excitation sources: The sparse vector q obtained in the previous step contains contributions from both scattering and excitation sources. Each component of vector q... This represents the hypothetical source intensity at the corresponding position on the Cartesian coordinate system of the plate surface. Centered on the known position coordinates (x0, y0) of the excitation source, the area covered by the piezoelectric sensor is defined as the detection blind zone. Within this zone, due to the high degree of overlap between the excitation and scattered signals, it is not considered an effective defect detection area. Components in q located within the aforementioned blind zone index range are removed and assigned a value of zero, resulting in a pure scattering source sparse vector. This makes the subsequently calculated mask matrix more accurate and eliminates the interference of the excitation source on defect detection.

[0151] 4) Defect contour solution: The defect contour solution corresponding to "Stage 1" (3) is represented by first screening the scattering sources, removing weak assumed scattering sources with scattering intensity below the set threshold, and filtering weak false scattering sources caused by noise; then constructing the feature function m(x,y) of the internal region of the defect through the fast convex hull algorithm and the non-zero winding principle. This prior knowledge is used as the mask matrix M in the FWI inversion process in "Stage 2" to participate in the iteration process.

[0152] 5) Gradient calculation: Corresponding to the adjoint state method of FWI in step (2) of “Stage 2”, the acoustic wave equation is discretized into a sparse matrix equation and solved by the conjugate gradient method.

[0153] 6) Mask Constraint Update: In step (3) of “Stage 2”, the inversion noise in non-defect regions is suppressed by combining the defect contour feature function m(x,y). The feature function acts as a spatial filter in the inversion, allowing only parameter updates within the defect region to obtain the phase velocity values ​​of each measurement point in the final iteration. This constraint essentially utilizes the prior information on the defect location provided in Stage 1, significantly reducing the computational load.

[0154] 7) Wave velocity-thickness conversion: The dispersion relationship is obtained by solving the Rayleigh-Lamb dispersion equation. For each discrete thickness value, the wavenumber of each mode at different angular frequencies is calculated, and a three-dimensional dispersion relationship surface diagram is formed by linear interpolation. In the specific conversion implementation, for the Lamb wavenumber distribution at a specific discrete frequency obtained by actual measurement, the wavenumber distribution is directly mapped to the final plate thickness distribution by looking up a table.

[0155] Example: Detection of thinning defects in aerospace aluminum plates

[0156] Input: Material: 3mm thick 2024-T3 aluminum alloy; Defect: Thinning defect measuring 20mm (height) × 30mm (width) × 2mm (depth); Sampling: Horizontal and vertical intervals are both 4.5mm, covering an area of ​​60mm × 70mm, such as... Figure 6 As shown.

[0157] Output: The distribution of thinning defect depth output by this method is as follows: Figure 7 The second row of figures shows the method described in the first row (J. Rao, M. Ratassepp, Z. Fan, Guided wave tomography based on full waveform inversion, IEEE Transactions Ultrasonics, Ferroelectrics, Frequency Control. 63(2016) 737–745. https: / / doi.org / 10.1109 / TUFFC.2016.2536144.). This method fails to produce effective defect detection results. In contrast, the results of this method provide effective quantitative detection of rectangular thinning defects. The detected defect shape and depth closely approximate the actual defect situation, with an average absolute depth error of 0.12 mm and a contour overlap of 92.3%, demonstrating the effectiveness and superiority of this method.

Claims

1. A method for quantitative evaluation of the three-dimensional morphology of thinning defects in a plate under ultrasonic guided wave field undersampling scanning, characterized in that: Phase 1: Defect Localization and Contour Extraction Phase 1 includes three steps: First, a sparse model of the ultrasonic guided wave field is achieved by constructing a ray model; then, the location and intensity of the wave field scattering source are obtained by using basis tracing noise reduction; finally, the feature function of the defect region is defined by using the location of the scattering source to complete the solution representation of the defect contour. Phase 2: Deep Quantization Based on Full Waveform Inversion (FWI) Phase 2 comprises four steps. First, based on the frequency domain two-dimensional acoustic wave equation, a discrete difference model of wave field propagation is constructed. Then, an FWI objective function is constructed and the parameter update gradient is solved. Next, using the characteristic function obtained in Phase 1, the parameters are updated using a mask constraint optimization method, and the wave field velocity distribution is iteratively solved. Finally, using the dispersion relationship of guided waves in the plate, the wave field velocity distribution is converted into a plate thickness distribution, thereby realizing a quantitative assessment of the three-dimensional morphology of defects in the plate.

2. The method according to claim 1, characterized in that: (1) Sparse modeling of ultrasonic guided wave field: First, a ray propagation model of the ultrasonic guided wave field is constructed. The subscripts m and s indicate the measurement point and scattering source, respectively. The location of the measurement point is known, while the location of the scattering source is unknown. (x, y) represents the position of the wave field in a Cartesian coordinate system, which lies on the plane of the plate. Based on the straight ray model, the position of the measurement point (x, y) is defined. m ,y m ) and the assumed scattering source location (x) s ,y s Euclidean distance between Calculate the frequency domain guided wavefield equations: Symbol definition: ω: angular frequency; i: represents the imaginary unit p m (ω): The sound pressure intensity at the measurement point m at angular frequency ω, which is the Fourier transform of the measurement signal in the time domain; μ: Ultrasonic guided wave modes, including: antisymmetric mode, symmetric mode and horizontal shear mode; N s The total number of scattering sources can be set to any integer value greater than the number of measurement points; The equivalent complex amplitude of the set scattering source at position s at angular frequency ω and mode μ is a complex function of the variable ω in the frequency domain; The wave number of a wave with mode μ propagating along a straight path between the measurement point m and the scattering source s at an angular frequency ω can be obtained by solving the ultrasonic guided wave dispersion curve. In the above modeling, the measurement points can be set to be uniformly distributed in the detection area with a spacing greater than half the wavelength; the possible locations of the scattering sources can be set to be randomly distributed in the detection area. (2) Basis Pursuit Denoising (BPDN) Solution: Let N m The total number of measurement points; the total number of scattering sources N s The total number of measurement points N m Under the given conditions, N can be obtained. m N equations of the same form as (1); m The equations are integrated and represented as matrix equations: Given p and H, the problem of identifying the ultrasonic guided wave scattering source q caused by defects in the plate is equivalent to solving the inverse problem of (2), which is expressed as the following sparse optimization problem: Symbol definition: ‖·‖1: L1 norm (the sum of the absolute values ​​of the vector elements); ‖·‖2: L2 norm (the square root of the sum of the squares of the vector elements); p: Measurement point vector Dimension N m ×1, corresponding to p at each measurement point in equation (1) m A vector composed of (ω), where each component represents the sound pressure intensity value at each measurement point at angular frequency ω, can be obtained by acquiring time-domain waveform signals and then performing Fourier transform. q: Scattering source vector Dimension N s ×1, corresponding to different positions in equation (1) The vector formed; σ: Noise margin, taken as 5% of the measured signal energy; H: Ultrasonic guided wave ray propagation matrix with mode μ and dimension N. m ×N s The corresponding elements are as follows, and the physical meaning of each element in the formula is the same as that in (1): The sparse optimization problem in equation (3) above can be solved using the SPGL1 algorithm; Once the actual ultrasonic guided wave scattering source is determined, its position is obtained from the preset relationship between the coordinates of the scattering source and the vector q. Since the ultrasonic guided wave scattering source originates from the influence of the thinning defect on the propagation of the ultrasonic guided wave, the scattering source is distributed within the boundary or region of the corrosion defect. Therefore, the position and outline of the corrosion defect are identified based on the position of the actual scattering source, and the geometric outline of the corrosion region will be further characterized in the next step. (3) Defect profile solution representation: First, scattering sources are screened, retaining those with a certain intensity and removing those with weak intensity; a threshold q is set. threshold = 0.1 × max(q), where max(q) represents taking the maximum value of the elements in vector q, and the threshold intensity of the scattering source is taken as 0.1 times the maximum value of the elements in vector q; using the set threshold q threshold From the scattering source vector q obtained in the previous step, scattering sources with amplitudes above a preset threshold are selected to form a point set of effective scattering source locations. Then, convex hull calculation is used to model the distribution of scattering sources, and the location where the scattering sources are concentrated is the defect region; the Quickhull algorithm is used to generate the minimum convex hull, denoted as C. s S is the set of scattering source points surrounding the wave field plane. s For all points, the shortest convex closed curve is used. The loop integral is defined using the non-zero winding rule to determine whether the location of the scattering source is within the defect region. For each point (x, y) located in the Cartesian coordinate system where the wave field is located, the closed curve C is calculated in the polar coordinate system with itself as the origin. s The number of wraps around this point is calculated using the following formula: Finally, for the defective region, the feature function m(x,y) is defined as follows: In the setting of the characteristic function m(x,y), if W Cs If (x,y)≠0, then the corresponding point (x,y) is located within the defect area, and all other areas are outside the defect. For the internal region of the defect, m(x,y) is set to a value of 1; for the external region of the defect, m(x,y) is set to a positive real number close to 0, 0.

01. Symbol definition: ∮: Represents the loop integral; θ: Convex hull C s Loop integral variables on; Convex hull C s The number of times the windings are around the point (x,y).

3. The method according to claim 1, characterized in that: The feature function m(x,y) is used as prior knowledge and substituted into equation (12) to improve the performance of thickness loss quantization. The technical details for each step are as follows; (1) Acoustic wave equation modeling: In the aforementioned defined wavefield plane rectangular coordinate system, the frequency domain two-dimensional acoustic wave equation is as follows: Where v(x,y,ω) represents the wave velocity at position (x,y) in the wave field at angular frequency ω; p(x,y,ω) is the sound pressure at position (x,y) in the wave field at angular frequency ω, which is the frequency domain representation of the disturbance signal generated by the sound source, and is the Fourier transform of the time domain signal at that position; q(x,y,ω) is the sound source at position (x,y) in the wave field, which is the frequency domain representation of the input function driving the sound wave generation, and is a complex function of angular frequency ω; The wavefield is divided into a Cartesian coordinate system, and the discrete, densely packed grid points are denoted as {(x... i ,y j )|i∈{1,…,N x },j∈{1,…,N y }}, where N x and N y Let h be the number of discrete grid points in the x and y directions, respectively; x Let x be the discrete dimension in the x-direction, that is: i+1 -x i =h x (i∈{1,…,N x -1}), let h be... y The discrete dimension in the y-direction: y j+1 -y j =h y (j∈{1,…,N y -1}); Discretizing equation (7) using the finite difference space yields: For any given angular frequency, equation (8) can be expressed in matrix form as follows: Symbol definition: ⊙: Hadamard product (element-by-element multiplication); Y: The model parameter matrix, composed of the reciprocal of the square of the wave velocity at each location in the wave field, i.e., Y = (1 / v 2 (x i ,y j ,ω),i∈{1,…,N x },j∈{1,…,N y }, where N x and N y These are the number of discrete grid points in the x and y directions, respectively; P: Sound pressure field matrix, composed of sound pressure at each location in the wave field, i.e., P = (p(x) / (x)) i ,y j ,ω),i∈{1,…,N x },j∈{1,…,N y } represents the grid (x) i ,y j The sound pressure at angular frequency ω is obtained by the Fourier transform of the amplitude of the time-domain measured waveform at the corresponding measurement point. Q: Source term matrix, composed of the positions of sound sources in the wave field, i.e., Q = (q(x)). i ,y j ,ω),i∈{1,…,N x },j∈{1,…,N y }, where the non-zero elements correspond to known wavefield sound sources, and are complex functions of angular frequency ω, (x i ,y j ) indicates its location; The second-order difference matrix corresponds to the discretization of the x and y directions on the wave field plane, respectively. It is directly extracted from equation (8). The dimension corresponds to the grid division, and the element value corresponds to the coefficient of equation (8). The matrix has a dimension of N x ×N x , The matrix has a dimension of N y ×N y The specific form is as follows: express Transpose transformation; (2) FWI objective function and gradient calculation: Based on the acoustic wave equation (9), the node settings in the discrete grid correspond to the wave field measurement point settings, forming a definite second-order difference matrix. and The initial value of the model parameter matrix Y is set to be an identity matrix, and the sound source term Q is set to correspond to the actual ultrasonic guided wave field detection excitation signal. Therefore, by solving the matrix equation of (9), the simulation result of the wave field state can be obtained, that is, the simulated wave field sound pressure data matrix, denoted as P. sim Furthermore, the simulated wave field is continuously approximated to the actual wave field through the following steps, and the model parameter matrix Y continuously approximates the actual situation; firstly, the objective function J is defined as: Symbol definition: ‖·‖2: L2 norm; e: Residual vector, representing the difference between the measured and analog signals; P meas : Wavefield acoustic pressure measurement data matrix; P sim The intermediate quantity obtained by solving (9) is the simulated wave field sound pressure data matrix. R: Spatial sampling matrix, which extracts data corresponding to the measurement points from the simulated wave field sound pressure; Then solve the optimization problem: The gradient g for parameter update is obtained using the adjoint state method: Symbol definition: p adj The accompanying state variables are solved using the conjugate gradient method. diag(p): a diagonal matrix with diagonal elements in the form of a sound pressure matrix p vector; diag(p) H : Represents the conjugate transpose of diag(p); Re: The operation of taking the real part; (3) Mask constraint optimization parameter update: Since the wave velocity change caused by plate thickness loss only occurs in the defect region and can be represented by the prior knowledge obtained in Stage 1, a mask-constrained optimization strategy is adopted for parameter updating. The wave field and wave velocity distribution can be solved by iterative iteration. Combining the characteristic function obtained in Stage 1, the gradient update range is constrained as follows: thing(Y new )=thing(Y old )-η·vec(M)⊙g (12) Symbol definition: Y old The parameter matrix Y in formula (9) is composed of the reciprocals of the squares of the wave velocities at each position in the wave field. As the gradient g is updated, it iterates to the old parameter values ​​before the update. Y new As the gradient g is updated, the iteration proceeds to the updated parameter values. η: Step size parameter, controls the update speed of the parameter, and takes the value of a real number in (0,1); M: Wavefield mask matrix, used to refine the update gradient g of the model parameters. It is constructed based on the values ​​of the characteristic functions obtained from formula (6) in stage 1 at each position in the wavefield, i.e., M = (m(x)). i ,y j ),i∈{1,…,N x },j∈{1,…,N y }, where N x and N y These represent the number of discrete points in the x and y directions of the discrete grid points in the wavefield plane rectangular coordinate system, respectively; vec(M): Expands matrix M into a vector column by column; Corrosion quantification is achieved by transforming model parameters. The transformed model parameters Y are converted into phase velocities, and then converted into plate thickness distributions based on the dispersion relationship of the waveguide phase velocities. This dispersion relationship is calculated in (13) below. (4) Wave velocity-thickness conversion: By solving the Rayleigh-Lamb dispersion equations, the dispersion relationship between wave velocity and thickness in the plate wave field is obtained; the wave velocity distribution is then transformed to reflect the plate thickness distribution; the Rayleigh-Lamb dispersion equations are as follows: h z The thickness of an isotropic plate; k: indicates the wavenumber of the ultrasonic Lamb wave; α, β: coefficients, their formulas are: C L : Represents the velocity of longitudinal waves, which is an inherent property of the material; C T : Represents the velocity of a transverse wave, which is an inherent property of the material; Equation (13) yields the relationship between wavenumber and thickness at a certain angular frequency ω. Let h... z The thickness increment step is 0.01 mm. Numerical calculations are performed to obtain the wavenumber-dispersion correspondence of each discrete thickness within the plate thickness range. Then, through linear interpolation, a complete and continuous wavenumber-dispersion relationship curve containing the discrete thickness values ​​and the regions between them is obtained. The dispersion relationship curves under the A0 and S0 modes are calculated separately, and the three-dimensional dispersion relationship surface plots of the "wavenumber-frequency-thickness" of the A0 and S0 modes are summarized. At each discrete frequency, based on the three-dimensional dispersion relationship of "wavenumber-frequency-thickness", the plate thickness distribution can be obtained from the wavenumber distribution by looking up a table.

Images