Method and device for detecting defects of a curved weld based on topological constraints of phase space

Abstract

The present application relates to a kind of based on phase space topological constraint's curved surface weld defect detection method and device.The method includes to the measured curved surface weld test piece is applied single point transient pulse excitation, collection excitation force signal and vibration response signal, response signal is reconstructed and extracts measured recursion quantization analysis index vector in phase space;Candidate function library is constructed, and multiple candidate hysteresis dynamic equation is generated by sparse regression;The output signal of each candidate equation is reconstructed again in phase space, extracts predicted recursion quantization analysis index vector and calculates the phase space topological deviation between it and measured recursion quantization analysis index vector;With time domain fitting error as threshold, with topological deviation minimum as criterion to determine final equation;According to measured recursion quantization analysis index vector, the structural characteristics and coefficient characteristics of final equation are constructed defect discriminant feature vector, and the probability of existence of defect is output by probability classification model.The method improves the sensitivity and explainability of curved surface weld defect detection.

Description

Technical Field

[0001] This invention relates to the field of nondestructive testing technology, and in particular to a method and apparatus for detecting defects in curved surface welds based on phase space topological constraints. Background Technology

[0002] In recent years, curved surface components (such as pressure vessels, engine casings, and composite material skins) have been widely used in aerospace and high-end equipment. These components are mainly connected by welds, and the early detection of microcracks in the welds is directly related to property safety and equipment operational safety. Closed fatigue cracks repeatedly come into contact and open and close under stress, exhibiting a typical "breathing effect," which leads to strong nonlinear characteristics and hysteresis behavior in the structural vibration response. Traditional linear detection methods (natural frequency shift, modal shape changes, etc.) have extremely low sensitivity to changes in interface contact state. In the closed state, defects do not significantly change the linear stiffness and damping parameters of the structure, resulting in a high rate of missed detection. Summary of the Invention

[0003] This invention provides a method and apparatus for detecting defects in curved surface welds based on phase space topological constraints, in order to overcome at least one of the above-mentioned technical problems existing in the prior art.

[0004] To achieve the above objectives, the embodiments of the present invention adopt the following technical solutions: In a first aspect, the present invention provides a method for detecting defects in curved surface welds based on phase space topological constraints, comprising: A single-point transient pulse excitation was applied to the weld specimen of the curved surface of the pipeline under test, and the excitation force signal and vibration response signal were acquired simultaneously. The acquired excitation force signal and vibration response signal were preprocessed by windowing, high-pass filtering and amplitude normalization to obtain preprocessed excitation force signal and preprocessed vibration response signal. The delay time and embedding dimension are determined based on the preprocessed vibration response signal, and the phase space trajectory matrix of the measured vibration response signal is constructed. A recursive matrix is ​​established based on the phase space trajectory matrix, and the measured recursive quantitative analysis index vector is extracted. A candidate function library containing polynomial terms, hysteresis basis function terms, differential operator terms, and cross-coupling terms is constructed. Using the preprocessed excitation force signal as input, sparse regression is performed on multiple candidate sparse regularization parameters to obtain multiple candidate hysteresis dynamic equations. Based on each candidate hysteresis dynamic equation, a predicted vibration response signal is generated, and the time-domain fitting error between the predicted vibration response signal and the preprocessed vibration response signal is calculated. The predicted vibration response signal is then reconstructed in phase space, and a predicted recursive quantization analysis index vector is extracted. The phase space topology deviation between the predicted recursive quantization analysis index vector and the measured recursive quantization analysis index vector is calculated. Among the candidate hysteresis dynamic equations whose time-domain fitting error is less than or equal to a preset error threshold, the final sparse regularization parameter and the final hysteresis dynamic equation are determined based on the criterion of minimizing the phase space topological deviation. Based on the measured recursive quantitative analysis index vector and the structural and coefficient features of the final hysteresis dynamic equation, a defect discrimination feature vector is constructed, and the probability of defect existence is obtained through a pre-trained probability classification model.

[0005] In one possible implementation of the first aspect, the delay time is determined by the first local minimum of the average mutual information function; the embedding dimension is determined by the dimension corresponding to the first time the proportion of false nearest neighbors is less than a preset threshold; the neighborhood radius of the recursion matrix is ​​determined according to the quantile of the distance between the phase space trajectory points, or determined in reverse according to the target recursion rate; the recursive quantification analysis index vector includes at least one or more of the following: recursion rate, determinism, average diagonal length, maximum diagonal length, divergence, diagonal length distribution entropy, laminar flow rate, and dwell time.

[0006] In one possible implementation of the first aspect, the highest order of the polynomial term does not exceed 5; the hysteresis basis function term includes Bouc-Wen type hysteresis basis function form and Duhem type hysteresis basis function form; the cross-coupling term is the product term of the polynomial term and the hysteresis basis function term.

[0007] In one possible implementation of the first aspect, the preprocessed excitation force signal is used as input to perform sparse regression on multiple candidate sparse regularization parameters to obtain multiple candidate hysteresis dynamic equations, including: The NSGA-III multi-objective optimization algorithm is used to filter the candidate function items in the candidate function library; The chromosome of the NSGA-III multi-objective optimization algorithm adopts a binary encoding method, and each gene bit indicates whether the corresponding candidate function item in the candidate function library is selected into the equation structure; The initial population is generated by a combination of random initialization and prior initialization. The random initialization sets each candidate function term to 1 with a preset probability, and the prior initialization sets the constant term, linear term, excitation input term and at least one hysteresis basis function term to 1. The crossover operation uses either uniform crossover or single-point crossover, with a crossover probability of 0.9. The mutation operation uses bit-flip mutation, with a mutation probability of 1 / d to 0.1, where d is the number of candidate function terms; With minimizing temporal fitting error and phase space topological deviation as dual optimization objectives, non-dominated sorting and crowding distance calculation are performed, and the next generation population is selected according to non-dominated level and crowding distance. When the termination condition is met, output the Pareto front solution set; The candidate function terms corresponding to each solution in the Pareto front solution set are combined as the structure of the candidate hysteresis dynamic equation, and sparse regression is performed on the structure of the candidate hysteresis dynamic equation to obtain the candidate hysteresis dynamic equation.

[0008] In one possible implementation of the first aspect, performing sparse regression on the structure of the candidate hysteresis dynamics equation includes: The sparse regression solution is obtained using either the LASSO algorithm or the sequential threshold least squares algorithm; among them... When using the sequential threshold least squares algorithm, the least squares solution is used as the initial coefficients. Threshold truncation and partial least squares reconstruction are iteratively performed until the set of activated candidate options no longer changes or the preset maximum number of iterations is reached. The sparse truncation threshold τ is used as the initial coefficient. s The relationship between the candidate sparse regularization parameter λ and the parameter λ satisfies , This is the initial least squares coefficient matrix. This represents the maximum absolute value of the initial coefficient, and α is the threshold proportionality coefficient.

[0009] In one possible implementation of the first aspect, generating the predicted vibration response signal based on each of the candidate hysteresis dynamic equations includes: Using the initial state of the preprocessed vibration response signal as the initial condition, the fourth-order Runge-Kutta method or a variable-step ordinary differential equation solver is used to numerically integrate each of the candidate hysteresis dynamic equations to generate the predicted vibration response time series corresponding to each candidate hysteresis dynamic equation.

[0010] In one possible implementation of the first aspect, the defect discrimination feature vector includes one or more of the following: recursive quantization analysis index, the number of activated terms in the final hysteresis dynamic equation, the proportion of nonlinear terms, the proportion of hysteresis terms, the norm of equation coefficients, and the equivalent hysteresis dissipation index.

[0011] In one possible implementation of the first aspect, the probabilistic classification model is a Gaussian process classification model, a support vector machine probabilistic model, or a probabilistic neural network model; When the probabilistic classification model is a Gaussian process classification model, a radial basis function kernel function is used; The hyperparameters of the Gaussian process classification model are determined by maximizing the log-marginal likelihood function under the Laplace approximation.

[0012] In one possible implementation of the first aspect, after obtaining the probability of defect presence through a pre-trained probabilistic classification model, the method further includes: The state of the curved weld is determined based on the comparison between the probability of the defect and the preset probability threshold. Output the probability of the defect's existence and the state of the curved weld.

[0013] Compared with the prior art, the present invention has at least the following beneficial effects: The present invention provides a method for detecting defects in curved welds based on phase space topological constraints. This method automatically discovers the hysteresis dynamic equation through sparse regression, overcoming the limitations of traditional methods that rely on prior models. Furthermore, it determines the final equation based on minimizing the phase space topological deviation between the predicted vibration response and the measured vibration response, solving the problem that existing methods may cause phase space structure distortion due to selecting solutions solely based on time-domain fitting errors. Finally, it integrates the structural and coefficient characteristics of the equation with recursive quantitative analysis indicators, and outputs the probability of defect existence through a probabilistic classification model, combining the interpretability of the dynamic model with the quantitative reliability of defect determination.

[0014] Secondly, the present invention provides a surface weld defect detection device based on phase space topological constraints, comprising: The signal acquisition and preprocessing unit is used to apply single-point transient pulse excitation to the weld specimen of the curved surface of the pipe under test, simultaneously acquire the excitation force signal and the vibration response signal, and perform windowing truncation, high-pass filtering and amplitude normalization preprocessing on the acquired excitation force signal and the vibration response signal respectively to obtain the preprocessed excitation force signal and the preprocessed vibration response signal. The phase space reconstruction unit is used to determine the delay time and embedding dimension based on the preprocessed vibration response signal, construct the phase space trajectory matrix of the measured vibration response signal, establish a recursive matrix based on the phase space trajectory matrix, and extract the measured recursive quantitative analysis index vector. The equation identification unit is used to construct a candidate function library containing polynomial terms, hysteresis basis function terms, differential operator terms and cross-coupling terms. Using the preprocessed excitation force signal as input, it performs sparse regression on multiple candidate sparse regularization parameters to obtain multiple candidate hysteresis dynamic equations. The topology selection unit is used to generate predicted vibration response signals based on each of the candidate hysteresis dynamic equations, calculate the time-domain fitting error between the predicted vibration response signal and the preprocessed vibration response signal, reconstruct the phase space of the predicted vibration response signal again, extract the predicted recursive quantization analysis index vector, calculate the phase space topology deviation between the predicted recursive quantization analysis index vector and the measured recursive quantization analysis index vector, and determine the final sparse regularization parameter and the final hysteresis dynamic equation based on the minimum phase space topology deviation among the candidate hysteresis dynamic equations whose time-domain fitting error is less than or equal to a preset error threshold. The defect discrimination unit is used to construct a defect discrimination feature vector based on the measured recursive quantitative analysis index vector and the structural features and equation coefficient features of the final hysteresis dynamic equation, and to obtain the probability of defect existence through a pre-trained probability classification model.

[0015] Thirdly, the present invention provides an electronic device comprising: at least one processor and at least one memory, wherein the memory stores computer-readable instructions; the computer-readable instructions are executed by one or more of the processors, causing the electronic device to implement the method for detecting surface weld defects based on phase space topological constraints as in any implementation of the first aspect.

[0016] Fourthly, the present invention provides a storage medium having a computer-executable program stored thereon, the computer-executable program being used to cause a computer to execute a method for detecting defects in curved welds based on phase space topological constraints as described in any implementation of the first aspect.

[0017] Understandably, the beneficial effects achieved by the apparatus of the second aspect, the electronic device of the third aspect, and the storage medium of the fourth aspect provided above can be referred to the beneficial effects of the first aspect and any of its possible design embodiments, which will not be repeated here. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention; Figure 2 A flowchart illustrating a method for detecting defects in curved welds based on phase space topological constraints, provided in an embodiment of the present invention; Figure 3 This is a structural block diagram of a curved surface weld defect detection device based on phase space topological constraints, provided in an embodiment of the present invention. Detailed Implementation

[0020] The technical solutions of the embodiments of the present invention will be described below with reference to the accompanying drawings. In the description of the present invention, unless otherwise stated, " / " indicates that the objects before and after are in an "or" relationship. For example, A / B can represent A or B. The "or" in the present invention is merely a description of the relationship between the related objects, indicating that three relationships can exist. For example, A or B can represent: A alone, A and B simultaneously, and B alone. A and B can be singular or plural. Furthermore, in the description of the present invention, unless otherwise stated, "multiple" refers to two or more. "At least one of the following" or similar expressions refer to any combination of these items, including any combination of single or plural items.

[0021] Furthermore, to facilitate a clear description of the technical solutions of the embodiments of the present invention, the terms "first" and "second" are used in the embodiments of the present invention to distinguish identical or similar items with essentially the same function and effect. Those skilled in the art will understand that the terms "first" and "second" do not limit the quantity or execution order, and the terms "first" and "second" are not necessarily different.

[0022] In this embodiment of the invention, the terms "exemplary" or "for example" are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" or "for example" in this embodiment of the invention should not be construed as superior or more advantageous than other embodiments or designs. Specifically, the use of terms such as "exemplary" or "for example" is intended to present the relevant concepts in a concrete manner for ease of understanding.

[0023] In recent years, curved surface components (such as pressure vessels, engine casings, and composite material skins) have been widely used in aerospace and high-end equipment. These components are mainly connected by welds, and the early detection of microcracks in the welds is directly related to property safety and equipment operational safety. Closed fatigue cracks repeatedly come into contact and open and close under stress, exhibiting a typical "breathing effect," which leads to strong nonlinear characteristics and hysteresis behavior in the structural vibration response. Traditional linear detection methods (natural frequency shift, modal shape changes, etc.) have extremely low sensitivity to changes in interface contact state. In the closed state, defects do not significantly change the linear stiffness and damping parameters of the structure, resulting in a high rate of missed detection.

[0024] In view of this, on the one hand, embodiments of the present invention provide a method for detecting defects in curved welds based on phase space topological constraints, comprising: applying a single-point transient pulse excitation to the tested curved weld specimen of a pipeline, simultaneously acquiring excitation force signals and vibration response signals, and performing preprocessing on the acquired excitation force signals and vibration response signals by windowing truncation, high-pass filtering, and amplitude normalization to obtain preprocessed excitation force signals and preprocessed vibration response signals; determining the delay time and embedding dimension based on the preprocessed vibration response signals to construct a phase space trajectory matrix of the measured vibration response signals; establishing a recursive matrix based on the phase space trajectory matrix and extracting the measured recursive quantitative analysis index vector; constructing a candidate function library containing polynomial terms, hysteresis basis function terms, differential operator terms, and cross-coupling terms, and using the preprocessed excitation force signals as input, performing sparse regression on multiple candidate sparse regularization parameters to obtain... Multiple candidate hysteresis dynamic equations are obtained; a predicted vibration response signal is generated based on each candidate hysteresis dynamic equation, and the time-domain fitting error between the predicted vibration response signal and the preprocessed vibration response signal is calculated; the phase space of the predicted vibration response signal is reconstructed again and a predicted recursive quantization analysis index vector is extracted, and the phase space topological deviation between the predicted recursive quantization analysis index vector and the measured recursive quantization analysis index vector is calculated; among the candidate hysteresis dynamic equations whose time-domain fitting error is less than or equal to a preset error threshold, the final sparse regularization parameter and the final hysteresis dynamic equation are determined based on the minimum phase space topological deviation; a defect discrimination feature vector is constructed based on the structural features and equation coefficient features of the measured recursive quantization analysis index vector and the final hysteresis dynamic equation, and the probability of defect existence is obtained through a pre-trained probabilistic classification model.

[0025] This invention provides a method for detecting defects in curved welds based on phase space topological constraints. It automatically discovers hysteresis dynamic equations through sparse regression, overcoming the limitations of traditional methods that rely on prior models. Furthermore, it determines the final equation based on minimizing the phase space topological deviation between the predicted and measured vibration responses, solving the problem that existing methods may distort the phase space structure by relying solely on time-domain fitting errors. Finally, it integrates the structural and coefficient characteristics of the equations with recursive quantitative analysis indicators, outputting the probability of defect presence through a probabilistic classification model, thus combining the interpretability of the dynamic model with the quantitative reliability of defect determination.

[0026] In some embodiments, the method for detecting defects in curved welds based on phase space topological constraints provided by the present invention can be executed by any electronic device 20 with data processing capabilities, such as a general-purpose computer, personal computer, laptop computer, switch, or tablet computer. The specific implementation of the electronic device 20 is not limited here.

[0027] Figure 1 A schematic diagram of the hardware structure of an electronic device provided in an embodiment of the present invention is shown. The electronic device 20 includes a processor 210, a memory 220, and a communication interface 230.

[0028] Processor 210 may include one or more processing cores. Processor 210 connects to various parts within electronic device 20 using various interfaces and lines, and performs various functions and processes data of electronic device 20 by running or executing instructions, programs, code sets, or instruction sets stored in memory 220, and by calling data stored in memory 220. Optionally, processor 210 may be implemented using at least one of the following hardware forms: Central Processing Unit (CPU), Graphics Processing Unit (GPU), Digital Signal Processing (DSP), Field-Programmable Gate Array (FPGA), and Programmable Logic Array (PLA).

[0029] The memory 220 may include random access memory (RAM) or read-only memory (ROM). Optionally, the memory 220 may include a non-transitory computer-readable storage medium. The memory 220 may be used to store instructions, programs, code, code sets, or instruction sets. The memory 220 may include a program storage area. This program storage area may store instructions for implementing an operating system, instructions for implementing at least one function, instructions for implementing the various method embodiments described above, etc.

[0030] Communication interface 230 is used to communicate with other devices, equipment or communication networks, such as data storage devices, image processing devices or Ethernet, wireless access network (RAN), wireless local area network (WLAN), etc.

[0031] In terms of physical implementation, the aforementioned devices (such as processor 210, memory 220, and communication interface 230) can each be devices within the same device (such as a laptop computer). Alternatively, at least two of these devices can be located within the same device, i.e., as different devices within the same device, similar to the deployment of devices or components in a distributed system.

[0032] It is understood that the structure illustrated in this embodiment does not constitute a specific limitation on the electronic device 20. In other embodiments of the present invention, the electronic device 20 may include more or fewer components than illustrated, or combine some components, or split some components, or have different component arrangements. The illustrated components may be implemented in hardware, software, or a combination of software and hardware.

[0033] The following description, in conjunction with the accompanying drawings, illustrates a method for detecting defects in curved welds based on phase space topological constraints, as provided in an embodiment of the present invention.

[0034] like Figure 2 As shown in the figure, the present invention provides a method for detecting defects in curved surface welds based on phase space topological constraints. This method may include, but is not limited to, the following: S1: Apply a single-point transient pulse excitation to the weld seam specimen of the pipe under test, and simultaneously acquire the excitation force signal and vibration response signal. Then, perform windowing truncation, high-pass filtering and amplitude normalization preprocessing on the acquired excitation force signal and vibration response signal to obtain the preprocessed excitation force signal and preprocessed vibration response signal.

[0035] The specimens in the embodiments of the present invention may be, but are not limited to, curved steel plate specimens of 300mm×60mm×2mm.

[0036] The embodiments of the present invention may, but are not limited to, use a hammer to apply a single-point transient pulse excitation, with the excitation direction in the negative Z-axis direction and the response acquisition direction in the positive Z-axis direction. This single-point transient pulse excitation can excite the vibration response of the structure over a wide frequency range, and is particularly suitable for capturing the nonlinear hysteresis behaviors such as contact, separation, and frictional dissipation exhibited by closed cracks under vibrational loads.

[0037] According to the parameter information of the defective specimen, the laser focal length can be set by the VibroRemoteConsole laser controller, and the vibration response can be recorded by the laser vibration meter acquisition system. Channel 1 is set as the input channel (force signal) with a sensitivity of 2.42mV / N; Channel 2 is set as the output channel (response signal).

[0038] The key parameters of the laser vibrometer acquisition system can be set as follows: The sampling frame rate is 5.12 kHz; the spectral line count is 1600; and the bandwidth is 2 kHz.

[0039] The linear averaging mode excitation was used, with a trigger delay of -20ms, a trigger threshold of 0.4%, and a trigger hysteresis of 0.2%. The linear averaging mode excitation effectively suppressed random noise interference and improved the signal-to-noise ratio by averaging multiple excitations. The trigger delay was set to a negative value to ensure that the signal prior to the excitation moment could be acquired as a baseline reference.

[0040] After obtaining the excitation force signal and vibration response signal, they are preprocessed respectively.

[0041] Specifically, embodiments of the present invention may, but are not limited to, perform the following preprocessing on the excitation force signal and the vibration response signal respectively: (1) Adding a window to cut off Since the excitation force signal and the vibration response signal have different data characteristics, different types of window functions are applied in this embodiment of the invention, as follows: The force signal is a transient impact signal. In this embodiment of the invention, a rectangular window function is applied to the excitation force signal, which can completely preserve the original characteristics of the impact waveform.

[0042] The response signal is a damped oscillation signal. In this embodiment of the invention, an exponential window function (main window offset of 5ms, damping constant of 100ms) is applied to the vibration response signal, which can effectively suppress noise amplification at the tail of the signal, reduce spectral leakage, and improve the signal-to-noise ratio.

[0043] (2) High-pass filtering In this embodiment of the invention, a high-pass filter is applied to the two types of windowed signals, with 3dB frequency 1 at 5Hz and 3dB frequency 2 at 1.2kHz. The high-pass filter removes low-frequency environmental vibrations and DC offset components, while retaining mid-to-high frequency components related to the weld structure's dynamic characteristics. The lower cutoff frequency of 5Hz is used to eliminate extremely low-frequency interference such as baseline drift and temperature drift; the upper cutoff frequency of 1.2kHz is used to filter out high-frequency electronic noise while retaining the main frequency band where the nonlinear response of weld defects is located.

[0044] (3) Amplitude normalization The embodiments of the present invention perform amplitude normalization processing on the filtered signal to eliminate amplitude differences caused by different excitation intensities and different sensor sensitivities, providing a unified dimensional benchmark for subsequent data processing.

[0045] After the above preprocessing, the preprocessed excitation force signal and the preprocessed vibration response signal are obtained.

[0046] S2: Determine the delay time and embedding dimension based on the preprocessed vibration response signal, and construct the phase space trajectory matrix of the measured vibration response signal; establish a recursive matrix based on the phase space trajectory matrix, and extract the measured recursive quantitative analysis index vector.

[0047] Specifically, according to Takens' embedding theorem, for a one-dimensional time series y(t), when the embedding dimension m satisfies m≥2d+1 (d is the fractal dimension of the original dynamic system), there exists a differential homeomorphic mapping from the original dynamic system to the reconstructed phase space. That is, the reconstructed phase space trajectory has the same dynamic topology as the original system in the sense of differential homeomorphism. Based on this, this embodiment of the invention maps the one-dimensional response signal to a high-dimensional phase space, reconstructing a phase space trajectory topologically equivalent to the original dynamic system, which serves as the input for subsequent recursive quantization analysis.

[0048] Among them, delay time and embedding dimension are two core parameters for phase space reconstruction. Embodiments of the present invention can determine the delay time and the embedding dimension in the following ways, but are not limited to: (1) Determination of delay time This invention employs the average mutual information method to determine the optimal delay time. The average mutual information function measures the statistical dependence between the signal y(t) and its delayed signal y(t+τ), and is defined as: .

[0049] Among them, p(y(t) i )) and p(y(t) j +τ)) are the univariate probability density functions (obtained through histogram method or kernel density estimation method, usually by dividing the signal value range into M equal intervals and counting the proportion of sampling points in each interval to the total number of points), p(y(t) i ),y(t j +τ)) is a two-dimensional joint probability density function (the joint frequency ratio of statistical signal values ​​falling in the i-th interval and delayed signal values ​​falling in the j-th interval), both of which can be directly estimated from the measured time series.

[0050] The optimal delay time is defined as the time when the AMI(τ) curve first reaches a local minimum, with the search range limited to [1, 100] sampling points. The physical meaning of this criterion is as follows: when τ is too small, y(t) and y(t+τ) are highly correlated, and the phase space trajectory is concentrated near the main diagonal, failing to fully develop the system's dynamic structure; when τ is too large, the two tend to be statistically independent, the phase space trajectory is too diffuse, and noise amplification is introduced. The first local minimum of AMI(τ) corresponds to the balance point between the two, representing "maximum decorrelation and minimum redundancy."

[0051] (2) Determining the embedding dimension This invention employs a pseudo-nearest neighbor method to determine the optimal embedding dimension. When the embedding dimension m is too low, points that are adjacent in the real phase space but not adjacent in the low-dimensional projection will be incorrectly identified as nearest neighbors; these points are called pseudo-nearest neighbors. As the embedding dimension increases, the proportion of pseudo-nearest neighbors gradually decreases and tends to stabilize. The dimension that first falls below a preset threshold is the optimal embedding dimension.

[0052] Specifically, after determining the delay time τ, a phase space state vector is constructed for the candidate embedding dimension m: .

[0053] In the formula, i=1,2,…,N-mτ.

[0054] For any state vector, search for its nearest neighbor among all other state vectors except itself. : .

[0055] To avoid the influence of false nearest neighbors caused by temporally adjacent sampling points, this embodiment of the invention introduces the Theiler window constraint |ji|>w. The window width w is determined based on the normalized autocorrelation function: the normalized autocorrelation function of the preprocessed response signal is calculated, and the sampling interval corresponding to the first decay to below 1 / e or the first crossing of zero is denoted as w. a Let w = max(τ, w a The sampling point is limited to the range of 1τ to 3τ. This eliminates sampling points that are temporally adjacent but whose dynamic states are not truly similar, reducing the impact of temporal correlation in the pseudo-nearest neighbor judgment.

[0056] Calculate the distance between the state vector and its nearest neighbor in the m-dimensional phase space: .

[0057] When the embedding dimension increases from m to m+1, calculate the distance increment in the direction of the new dimension: .

[0058] If the following condition is met, then the nearest neighbor is determined to be a pseudo-nearest neighbor: .

[0059] In the formula, R tol This is a pre-set relative distance threshold based on the noise level of the response signal and reconstruction stability. It is typically set to 10~20. The higher the noise level, the greater the relative distance threshold. tol The value should be appropriately increased to avoid noise being misjudged as an effective dynamic structure; the cleaner the signal, the better. tol It can be reduced appropriately to retain more detailed information.

[0060] The proportion of pseudo nearest neighbors under different candidate embedding dimensions is statistically analyzed. .

[0061] In the formula, Let m be the number of points that are determined to be pseudo-nearest neighbors in the embedding dimension m. The number of points in the phase space involved in the calculation.

[0062] The minimum embedding dimension corresponding to the first time FNN(m) falls below a preset proportion threshold (usually 1%~5%) and tends to stabilize is determined as the final embedding dimension.

[0063] After determining the delay time τ and the embedding dimension m, this embodiment of the invention maps the one-dimensional preprocessed vibration response signal y(t)={y1,y2,…,yN} to an m-dimensional phase space according to Takens' embedding theorem, constructing a high-dimensional phase space trajectory matrix X: .

[0064] In the formula, y i The amplitude of the preprocessed vibration response signal at the i-th sampling time is represented by τ, where i = 1, 2, ..., N; τ is the delay time; m is the embedding dimension; and L = N - (m - 1)τ is the total number of valid state vectors after reconstruction. For the reconstructed phase space trajectory matrix, its i-th row [y i ,y i+τ ,…,y i+(m−1)τ ] represents a state point in phase space.

[0065] The L state points arranged in chronological order in this matrix constitute a high-dimensional phase space trajectory that is homeomorphic to the original dynamic system, preserving the implicit nonlinear dynamic characteristics in the original time series, and providing a data foundation for the subsequent establishment of recursive matrices and the extraction of recursive quantitative analysis indicators.

[0066] After constructing the phase space trajectory matrix, a recursive matrix is ​​established based on the phase space trajectory matrix, and the measured recursive quantitative analysis index vector is extracted.

[0067] Specifically, the recursion matrix is ​​a binary representation of the recursive relationships between all pairs of state points in the phase space trajectory. In this embodiment of the invention, the Euclidean distance between any two state points in the phase space is calculated and compared with a preset neighborhood radius threshold to establish the recursion matrix: .

[0068] In the formula, X i and X j X represents the i-th and j-th state point vectors in the phase space trajectory matrix, respectively. i =[yi ,y i+τ ,…,y i+(m-1)τ ], X j =[y j ,y j+τ ,…,y j+(m-1)τ ]; The Euclidean distance between two state points is calculated using the following formula: ; ε is a preset neighborhood radius threshold. If the Euclidean distance between two state points is less than or equal to ε, then the two state points are considered to have a recursive relationship in the phase space (i.e., the system returns to a neighboring state during the evolution process); otherwise, there is no recursive relationship. H(.) is the Heaviside step function, specifically: .

[0069] This function acts as a binarization operator, mapping continuous distance metrics to discrete elements of a recursive matrix: it assigns a value of 1 when the distance is less than or equal to a threshold (recursion occurs), and assigns a value of 0 when the distance is greater than the threshold (recursion does not occur).

[0070] Specifically, the neighborhood radius ε is a core parameter for establishing the recursive matrix, and its value directly affects the sensitivity and robustness of the recursive quantification analysis index. In this embodiment of the invention, it can be determined in the following two ways: Method 1: Quantile method. Calculate the Euclidean distance between all pairs of state points in the phase space trajectory matrix. (i,j=1,2,…,L,i≠j), forming the distance set D={d ij} i≠j Let the p-quantile of this set be the neighborhood radius ε, i.e., ε = Q. p (D), where Q p (·) is the quantile function, and p is preferably taken as 10%~15%. The advantage of this method is that the neighborhood radius and the natural scale of the phase space trajectory are adaptively matched, avoiding the subjectivity of manually setting the threshold; it is suitable for scenarios where the scale of the phase space trajectory changes greatly with the working conditions or sensor position.

[0071] Method 2: Target Recursion Rate Inverse Method. Set the target recursion rate RR. target (Typically 1%~5%), use binary search or other numerical optimization methods to find the recurrence rate RR(ε) that makes the actual recurrence rate RR(ε) equal to RR. targetThe neighborhood radius ε. The advantages of this method are: it provides a unified benchmark for the recurrence rate among different detection signals, ensuring the comparability of recurrence quantification analysis indicators under different specimens and detection conditions; it is suitable for scenarios that require horizontal comparison of multiple detection signals, such as the comparative analysis of RQA (Recurrence Quantification Analysis) indicators between healthy samples and defective samples.

[0072] It should be noted that either of the above two methods can be used, or the choice can be made flexibly according to the actual signal characteristics and detection requirements, without any restrictions.

[0073] The recursion matrix establishes the connection between the geometric structure of the phase space trajectory (the proximity relationship between state points) and the dynamic characteristics of the system. Based on the recursion matrix, this embodiment of the invention extracts a series of quantitative indicators to quantitatively describe the statistical characteristics of the recursion structure. These quantitative indicators may include, but are not limited to: the recursion rate (RR) extracted based on the global recursion point density; and deterministic diagonal length (DET), average diagonal length (L), and maximum diagonal length (L') extracted based on the distribution of the diagonal structure. max Divergence DIV and diagonal length distribution entropy ENTR; laminar flow rate LAM and dwell time TT extracted based on the distribution of vertical line structure.

[0074] The above quantitative indicators are combined into a vector Q. meas =[RR,DET,L,L max ,DIV,ENTR,LAM,TT] T , which serves as a description of the phase space topological features of the measured vibration response signal.

[0075] S3: Construct a candidate function library containing polynomial terms, hysteresis basis function terms, differential operator terms, and cross-coupling terms. Using the preprocessed excitation force signal as input, perform sparse regression on multiple candidate sparse regularization parameters to obtain multiple candidate hysteresis dynamic equations.

[0076] Specifically, in this embodiment of the invention, the hysteresis dynamic equation to be discovered is expressed in explicit state-space form as follows: .

[0077] In the formula, X(t) is the system state variable vector; This is a matrix of candidate function libraries composed of state variables; Let be the sparse coefficient matrix to be identified.

[0078] The system state variable vector is defined as follows: .

[0079] In the formula, y(t) is the vibration response signal (displacement response) after preprocessing; y(t) is the first time derivative of the vibration response signal (velocity response), which is numerically calculated from y(t) using the central difference method or the Savitzky-Golay filtering method; z(t) is the internal state variable (i.e., hysteresis variable) used to characterize the contact, opening and closing and hysteresis dissipation behavior of the weld defect interface, and its initial value can be 0.

[0080] Because the hysteretic behavior of a closed crack under vibration load cannot be determined solely by y(t) and Explicit description requires an additional internal variable to remember the system's loading history. This embodiment of the invention introduces z(t), enabling the candidate function library to include Bouc-Wen type hysteresis basis functions (such as...). , (etc.), thereby capturing the nonlinear contact and frictional dissipation behavior of the weld interface.

[0081] The candidate function library is the foundation of sparse regression, and its design directly determines the expressive power of the final discovery equation. The candidate function library constructed in this embodiment of the invention includes the following four categories: (1) Polynomial terms The polynomial terms are used to describe the linear and nonlinear algebraic relationships between system state variables, including constant terms, first-order and higher-order terms of state variables, and interaction terms between state variables: .

[0082] Since the nonlinearity of the weld structure mainly originates from contact nonlinearity, its nonlinearity intensity is typically no higher than order 5. Furthermore, excessively high polynomial orders can lead to an explosion in the size of the candidate function library, increasing the risk of overfitting and computational burden. Therefore, in this embodiment of the invention, the highest order of the polynomial terms does not exceed order 5, approximately 40 to 50 terms.

[0083] (2) Delayed basis function term The hysteresis basis function term is used to describe the hysteresis behavior during the contact and opening / closing processes of weld defect interfaces. Embodiments of this invention incorporate multiple possible forms of hysteresis basis functions into the candidate library simultaneously, including, but not limited to: Bouc-Wen type hysteresis basis function form: , Where n is the hysteresis evolution parameter (which can be 1 to 3), controlling the shape of the hysteresis loop.

[0084] Duhem-type hysteresis basis function form: , . Among them, sgn() is the symbolic function.

[0085] Generalized hysteresis basis function form: , , .

[0086] The evolution of the hysteresis internal state variable z(t) is related to velocity. For the driver, its update form can be: .

[0087] In the formula, β and γ are hysteresis shape parameters that control the amplitude and fullness of the hysteresis loop.

[0088] Bouc-Wen type hysteresis basis functions and These correspond to the hysteresis dissipation term (related to velocity direction and hysteresis amplitude) and hysteresis stiffness term (related to velocity direction) in the hysteresis restoring force, respectively. These two terms together determine the width, fullness, and symmetry of the hysteresis loop. Duhem-type hysteresis basis functions. and It focuses on describing the non-smooth mapping relationship between displacement and hysteresis variables, and is suitable for modeling the alternating "viscous-slip" behavior in interface contact.

[0089] (3) Differential operator terms The differential operator terms are used to describe the derivative relationships between system state variables, including , Etc. In this embodiment of the invention, the differential operator term is introduced, enabling the candidate function library to include damping and inertial terms.

[0090] (4) Cross-coupling terms The cross-coupling term, a product of a polynomial term and a hysteresis basis function term, describes the coupling between nonlinear effects and hysteresis effects during the opening and closing of the weld defect interface. This embodiment of the invention introduces this cross-coupling term to capture the interaction between the nonlinear contact stiffness and hysteresis dissipation of the defect interface.

[0091] In this embodiment of the invention, the initial candidate library size is approximately 100-150 items. To reduce computational burden and improve the stability of sparse regression, this embodiment of the invention first pre-screens the candidate library using stepwise regression.

[0092] Specifically, for each candidate basis function, using the state derivative component as the target variable, the partial F-statistic for univariate regression is calculated: ; ; .

[0093] In the formula, q is the target derivative vector to be fitted; This is the predicted value obtained by regression using only the j-th candidate basis function; The mean of the objective derivative; SSR j The sum of squares is the regression sum, reflecting the variance explained by the j-th candidate basis function; SSE j F represents the residual sum of squares, reflecting the variance that the j-th candidate basis function failed to explain; N is the number of sampling points; j This is a partial F-statistic; the larger the value, the stronger the explanatory power of the candidate basis function for the target derivative.

[0094] Press F j Sort all candidates in descending order, retain the 1 / 3 to 1 / 2 of items with the largest F-statistic, and reduce the candidate pool size to [missing value]. item.

[0095] Get the filtered Then, using the preprocessed excitation force signal as input, sparse regression is performed on multiple candidate sparse regularization parameters to obtain multiple candidate hysteresis dynamic equations.

[0096] In one feasible implementation, embodiments of the present invention may, but are not limited to, employ the NSGA-III multi-objective optimization algorithm to screen candidate function terms in the candidate function library in order to obtain the structure of candidate hysteresis dynamic equations.

[0097] Specifically, the filtered The candidate options are denoted as The optimal decision variable is a binary choice vector: .

[0098] Where δ j =1 indicates that the j-th candidate is selected into the governing equation structure, δ j =0 indicates that it is not selected.

[0099] The binary selection vector is optimized using the NSGA-III multi-objective optimization algorithm. The optimization objectives include: Time-domain fitting error target: ; .

[0100] In the formula, To predict the coefficient of determination between the vibration response and the measured vibration response, Let be the predicted vibration response value at the i-th sampling time. Let be the measured vibration response value at the i-th sampling point. This represents the average measured vibration response.

[0101] Phase space topological fidelity target: ; .

[0102] In the formula, The RQA index vector for the predicted signal. This is the RQA index vector of the measured signal. Let ε be the Euclidean norm, and let ε be a small positive number to prevent the denominator from being zero.

[0103] In this embodiment of the invention, the chromosomes of the NSGA-III multi-objective optimization algorithm can be encoded in binary form, with each chromosome representing a candidate equation structure selection vector b=[b1,b2,…,b…]. d ], where d is the number of candidate function terms after pre-screening, b j =1 indicates that the j-th candidate function term is selected into the candidate hysteresis dynamic equation, b j =0 indicates that the candidate function item was not selected.

[0104] The initial population of the NSGA-III multi-objective optimization algorithm can be generated by combining random initialization and prior initialization, as follows: Random initialization: Each candidate function term is set to 1 with a preset probability p0, preferably p0 is 0.1~0.3, to ensure that the initial equation structure is sparsity.

[0105] Prior initialization: Set the constant term, linear term, excitation input term and at least one hysteresis basis function term to 1 to ensure that the candidate equation has basic dynamic expression capability.

[0106] The crossover operation of the NSGA-III multi-objective optimization algorithm can employ uniform crossover or single-point crossover, with a crossover probability P. c The preferred value is 0.9.

[0107] The mutation operation of the NSGA-III multi-objective optimization algorithm uses bit-flip mutation. For each gene locus in the chromosome, with a mutation probability P... m Perform a toggle between 0 and 1, P m The optimal value is 1 / d to 0.11, where d is the number of candidate function terms. To avoid an empty equation structure, if all gene positions in the mutated individual are 0, at least one linear or lag term is randomly selected and set to 1.

[0108] In this embodiment of the invention, the NSGA-III multi-objective optimization algorithm first performs non-dominated sorting based on the two objective function values: individuals that are not simultaneously superior to other individuals are divided into the first non-dominated layer, and the remaining individuals are divided into subsequent non-dominated layers in turn.

[0109] For individuals in the same non-dominated layer, a reference point association mechanism is used to maintain population distribution diversity. In the case of a two-dimensional objective, crowding distance can also be calculated, and the degree of crowding of individuals can be determined according to the sum of the distances between adjacent individuals after normalization of the objective function, and individuals with larger crowding distances are given priority to be retained.

[0110] When the termination condition is met, output the Pareto front solution set.

[0111] Specifically, the termination condition in the embodiments of the present invention may include, but is not limited to, any of the following: (1) Reach the maximum number of generations of evolution (200~500 generations); (2) The change in the optimal phase space topological deviation of the Pareto front for several consecutive generations (20 to 50 generations) is less than the preset threshold δ (δ can be taken as 10). -4 ~10 -3 ); (3) The non-dominated solution set does not change for several consecutive generations (20 to 50 generations).

[0112] After the NSGA-III outputs a set of non-dominated Pareto front solutions, the candidate function terms corresponding to each solution in the Pareto front solution set are combined to form the structure of the candidate hysteresis dynamic equations. Then, sparse regression is performed on the structure of the candidate hysteresis dynamic equations to obtain the candidate hysteresis dynamic equations.

[0113] Specifically, in this embodiment of the invention, sparse regression is performed on each candidate individual on the Pareto front to obtain specific equation coefficients: .

[0114] In the formula, It is a state data matrix composed of measured vibration response signal, response velocity signal, excitation force signal, and hysteresis internal state variables. This is the corresponding state derivative data matrix. This is the candidate function library matrix. Let be the sparse coefficient matrix to be identified, and λ be the candidate sparse regularization parameter.

[0115] The larger the value of λ, the heavier the L1 (L1 norm, the sum of the absolute values ​​of all elements in a vector) penalty, the more coefficients are compressed to zero, and the sparser the equation structure; the smaller the value of λ, the more non-zero coefficients are retained, and the more complex the equation structure.

[0116] In one feasible implementation, embodiments of the present invention may, but are not limited to, employ the LASSO algorithm or the sequential threshold least squares algorithm for sparse regression solution.

[0117] Specifically, when using the LASSO algorithm, the above least squares problem with L1 regularization is solved directly.

[0118] When using the sequential threshold least squares algorithm, the least squares solution is first used as the initial coefficients: .

[0119] In the formula, The Moore-Penrose pseudoinverse of the candidate function library matrix.

[0120] Then, threshold truncation is performed: for the coefficients obtained in the k-th iteration... The absolute value is less than the preset threshold τ s The coefficients are set to zero. The set of indices of the retained non-zero coefficients is denoted as: .

[0121] Then, unregularized least squares fitting is re-performed on the retained candidates: .

[0122] In the formula, This is a candidate function library submatrix composed of the columns corresponding to the retained candidate options.

[0123] Repeat the threshold truncation and partial least squares reconstruction until the set of activated candidates no longer changes, or until the preset maximum number of iterations is reached (e.g., 200 times).

[0124] Sparse cutoff threshold τ s The following relationship exists between the candidate sparse regularization parameter λ and the parameter: .

[0125] In the formula, This is the initial least squares coefficient matrix. The initial coefficient represents the maximum absolute value, and α is the threshold scaling factor, ranging from 0.05 to 0.2. For detection signals with high noise levels, α can be appropriately increased; for signals with high signal-to-noise ratios and where weak nonlinear terms need to be retained, α can be appropriately decreased.

[0126] In one feasible implementation, the sparse truncation threshold τ in the embodiments of the present invention s It can also be determined through cross-validation.

[0127] Specifically, the sampled data is divided into a training segment and a validation segment according to time sequence. The training segment is used to identify candidate hysteresis dynamic equations, and the validation segment is used to calculate the time-domain fitting error and phase-space topological deviation of the predicted vibration response; in the candidate threshold set The performance was verified one by one, and the τ value that satisfies the time-domain fitting error constraint and minimizes the phase space topological deviation was selected. s As the final threshold.

[0128] After performing the above sparse regression on each candidate structure on the Pareto front, the corresponding equation coefficients are obtained, thus generating multiple candidate hysteresis dynamic equations. Each candidate equation simultaneously contains: Equation structure: Which candidate function terms are activated (with non-zero coefficients); Equation coefficients: The specific coefficient values ​​corresponding to each activation term.

[0129] S4: Generate predicted vibration response signals based on each of the candidate hysteresis dynamic equations, calculate the time-domain fitting error between the predicted vibration response signals and the preprocessed vibration response signals; and reconstruct the phase space of the predicted vibration response signals again and extract the predicted recursive quantization analysis index vector, calculate the phase space topological deviation between the predicted recursive quantization analysis index vector and the measured recursive quantization analysis index vector.

[0130] In this embodiment of the invention, the generation of predicted vibration response signals based on each of the candidate hysteresis dynamic equations can be understood as follows: For each candidate hysteresis dynamic equation, its corresponding coefficient matrix Substituting the structure into the state-space equations: .

[0131] In the formula, u is the preprocessed excitation signal (system input), and z is the hysteresis internal state variable.

[0132] The predicted vibration response time series is generated using a numerical integration method. The specific steps are as follows: Using the initial state of the measured signal as the initial condition: .

[0133] Where y(t0) is the initial displacement value of the preprocessed vibration response signal; v(t0) is the initial velocity value (i.e., The value at time t0); z(t0) is the initial hysteresis internal state variable (usually taken as 0, indicating that the weld interface is in a natural state without residual deformation at the initial time).

[0134] The candidate equations are numerically integrated using the fourth-order Runge-Kutta method (RK4) or a variable-step ordinary differential equation solver (such as MATLAB's ode45). The integration step size can be 1 / 10 of the sampling interval, or adaptively determined based on the numerical integration error.

[0135] After integration, the predicted vibration response time series corresponding to the candidate hysteresis dynamic equation is obtained. Its length is related to the measured vibration response time series. They must be identical so that subsequent point-by-point comparisons can be performed.

[0136] In this embodiment of the invention, the time-domain fitting error is used to measure the predicted vibration response time series. Compared with the measured vibration response time series Consistency at the waveform level.

[0137] Specifically, in this embodiment of the invention, the time-domain fitting error can be represented by the normalized root mean square error (NRMSE): .

[0138] in, The Euclidean norm (L2 norm) of a vector is defined as follows: .

[0139] It should be noted that the time-domain fitting error in the embodiments of the present invention can also be the root mean square error (RMSE), the normalized mean square error (NMSE), or determined by the coefficient of determination R. 2 The error index obtained from the conversion (E) t =1−R 2 Equivalent substitutions are not limited here.

[0140] While calculating the time-domain fitting error between the predicted vibration response signal and the preprocessed vibration response signal, this embodiment of the invention also applies the predicted vibration response time series. Perform phase space reconstruction and recursive quantization analysis identical to the measured vibration response in S2 to extract the predictive recursive quantization analysis index vector. Then calculate the predictive recursive quantitative analysis index vector. With the measured recursive quantitative analysis index vector The phase space topological deviation between them.

[0141] Specifically, the phase space topology deviation is a quantitative indicator that measures the difference between the candidate equation and the real system at the phase space level. In this embodiment of the invention, the phase space topology deviation is defined as the predictive recursive quantitative analysis index vector. Distance between the measured recursive quantitative analysis index vector and: .

[0142] In the formula, d(.) is the distance metric between recursive quantitative analysis index vectors. The distance metric can be standardized Euclidean distance, weighted Euclidean distance, or Mahalanobis distance.

[0143] Taking Mahalanobis distance as an example, the phase space topological deviation is expressed as: .

[0144] In the formula, The regularized covariance matrix is ​​estimated from the RQA index vector of the benchmark or training samples: ; ; .

[0145] In the formula, K is the number of baseline samples or training samples, preferably K≥30; Let K be the recursive quantitative analysis index vector for the Kth sample. η is the sample mean of the recursive quantitative analysis index vector; I is the regularization coefficient; and I is the identity matrix.

[0146] The introduction of the regularization term η·I guarantees the covariance matrix It is always invertible (positive definite), avoiding the problem of singular or ill-conditioned covariance matrix caused by insufficient training samples or high correlation between indicators.

[0147] Ultimately, each candidate hysteresis dynamic equation (corresponding to a λ value) was evaluated using two quantitative indicators: time-domain fitting error and phase-space topological deviation.

[0148] S5: Among the candidate hysteresis dynamic equations whose time-domain fitting error is less than or equal to a preset error threshold, the final sparse regularization parameter and the final hysteresis dynamic equation are determined based on the criterion of minimizing the phase space topological deviation.

[0149] Specifically, in this embodiment of the invention, the time-domain fitting error is used as the first screening threshold to eliminate candidate equations with severe time-domain distortion; the phase-space topological deviation is used as the second selection criterion to select the candidate equation with the smallest phase-space topological deviation from the qualified candidates in the time domain. The screening process can be understood as follows: First screening: Time-domain fitting error threshold filtering.

[0150] In this embodiment of the invention, the preset time-domain error threshold E th It can be determined based on, but is not limited to, the following factors: Noise level of the measured signal: the higher the noise level, the higher the E th The restrictions should be relaxed appropriately to avoid mistakenly eliminating valid candidate equations. Engineering experience: In weld defect detection scenarios, E th The value is typically between 0.05 and 0.20, depending on the stability of the specimen material and the testing conditions. Cross-validation or benchmarking: Pre-testing can be performed on healthy specimens to statistically analyze the distribution of time-domain fitting error, using E... th =μ E +2σ E (Mean plus twice the standard deviation) is used as the threshold.

[0151] Among all candidate equations that satisfy Candidate equations are directly eliminated.

[0152] The first screening establishes a basic physical credibility threshold, ensuring that the equations do not violate basic physical consistency at the time domain level.

[0153] The second screening: minimizing the phase space topological deviation for optimal selection.

[0154] After the first round of screening, the remaining candidate equations satisfy... That is, they all have basic time-domain fitting accuracy.

[0155] Among the time-domain qualified candidate equations, the phase space topological deviation D is selected. topo The candidate equation with the smallest (λ) is taken as the final equation. Mathematically, this is equivalent to searching from the subset of candidate equations that satisfy the time-domain constraints for the equation that minimizes D. topo The λ value that reaches its minimum value is the candidate sparse regularization parameter corresponding to that λ value, which is the final sparse regularization parameter, denoted as λ. ,and The corresponding candidate hysteresis dynamic equation was determined as the final hysteresis dynamic equation.

[0156] Because of D topo (λ) measures the difference between the candidate equation and the real system at the recursive structure level in phase space. Choose D. topo The equation with the smallest (λ) means that, while ensuring the basic time-domain fitting accuracy, we choose the equation that is closest to the real system at the level of dynamics.

[0157] After determining the final hysteresis dynamic equation, in order to improve the accuracy of coefficient estimation, embodiments of the present invention may, but are not limited to, further refine the selected set of activation terms. Perform unregularized least squares fitting: .

[0158] in It is a candidate function library submatrix containing only the activated candidates in the final equation.

[0159] By selecting the set of activation items Performing unregularized least squares fitting can eliminate the coefficient shrinkage bias caused by L1 regularization and obtain unbiased final coefficient estimates.

[0160] The final hysteresis dynamic equation outputs the following information: Equation structure: a combination of activated candidate function terms; Equation coefficients: The coefficient values ​​corresponding to each activation term; Final regularization parameters: .

[0161] S6: Construct a defect discrimination feature vector based on the measured recursive quantitative analysis index vector and the structural features and equation coefficient features of the final hysteresis dynamic equation, and obtain the probability of defect existence through a pre-trained probability classification model.

[0162] In this embodiment of the invention, the eight indicators in the measured recursive quantization analysis index vector describe the recursive behavior characteristics of the system in phase space from different dimensions, as follows: Recurrence rate (RR): Reflects the overall density of recursive points in the phase space trajectory, characterizing the degree of periodicity or regularity of the system's motion; Determinism (DET): Reflects the determinism and predictability of the system's evolutionary behavior; the higher the value, the more deterministic the system's behavior tends to be. Average diagonal length (L): reflects the average persistence of the diagonal structure in the recursive structure and is related to the predictable time scale of the system. Maximum diagonal length (L) max ): Reflects the length of the longest continuous evolution path during the system's evolution process and is an important indicator for measuring system stability; Divergence (DIV): The reciprocal of the maximum diagonal length, used to estimate the maximum Lyapunov exponent of the system, characterizing the degree of chaos in the system; Diagonal length distribution entropy (ENTR): reflects the complexity of the diagonal length distribution in a recursive structure; the higher the value, the more complex the system's dynamic behavior. Laminar flow rate (LAM): reflects the proportion of time a system's phase space trajectory stays near certain states, characterizing the intermittent nature of the system; Dwell time (TT): The average length of the vertical line, which measures the average duration for which the system remains in a laminar flow state.

[0163] When defects exist in the weld, the nonlinear stiffness and damping characteristics of the structure change, and the recursive structure of the phase space trajectory also changes accordingly. The above 8 RQA indicators will show a systematic shift, and therefore can be used as effective features for defect identification.

[0164] In this embodiment of the invention, the structural features of the equation reflect the "structure" of the equation itself. When defects exist in the weld, the hysteresis nonlinearity of the system is enhanced, and the number and proportion of nonlinear and hysteresis terms in the equation will increase accordingly. Specifically, this includes: (1) Number of activated items (N) active ) The total number of candidate function terms with non-zero coefficients in the final hysteresis dynamic equation is counted. Under healthy conditions, the dynamic behavior of the weld structure is close to linear, and the number of activated terms in the equation is relatively small (usually mainly linear terms and a small number of nonlinear terms). Under defective conditions, hysteresis nonlinearity and contact nonlinearity are enhanced, and more nonlinear terms and hysteresis basis function terms are activated to describe the interface opening and closing behavior, resulting in a significant increase in the number of activated terms.

[0165] (2) Proportion of nonlinear terms (R) nonlinear ) Calculate the proportion of nonlinear terms (polynomial terms of order ≥ 2) out of all activated terms in the final equation. This proportion reflects the overall level of nonlinearity in the system. The proportion of nonlinear terms is typically higher in defective states than in healthy states.

[0166] (3) Proportion of delayed terms (R) hysteresis ) The proportion of hysteresis basis function terms (including Bouc-Wen type, Duhem type, and generalized hysteresis basis function terms) to all activated terms in the final equation is calculated. This proportion directly reflects the relative importance of hysteresis behavior in system dynamics. The more significant the "breathing effect" of the weld closure crack, the higher the proportion of hysteresis terms.

[0167] In this embodiment of the invention, the coefficient characteristics of the equation reflect the "weight strength" of each term in the equation. When defects exist in the weld, the nonlinear stiffness and hysteresis dissipation of the system are enhanced, and the corresponding coefficient norm and dissipation index will increase accordingly. Specifically, this includes: (1) Equation coefficient norm Calculate all non-zero coefficients of the final hysteresis dynamic equation. L2 norm: .

[0168] The norm of the equation coefficients reflects the amplitude level of the overall dynamic response of the system. The presence of cracks will change the equivalent stiffness and damping of the system, and the norm of the coefficients will change accordingly.

[0169] (2) Norm of lag term coefficient Calculate the L2 norm of the coefficients of all hysteresis basis function terms in the final equation: .

[0170] in This is the index set of hysteresis basis function terms. This index specifically quantifies the strength of hysteresis behavior and is highly sensitive to the "breathing effect" of closed cracks. The more severe the defect, the stronger the hysteresis dissipation, and the larger the norm of the hysteresis term coefficients.

[0171] (3) Equivalent hysteresis dissipation index The energy dissipation capability of a system due to hysteresis nonlinearity can be comprehensively characterized by an approximate estimate over the time domain by the integral of the product of the hysteresis term coefficients and the magnitudes of the corresponding hysteresis basis function terms. .

[0172] in This represents the average amplitude of the j-th hysteresis basis function term over the sampling time. This index directly measures the frictional energy dissipation and contact nonlinearity at the weld interface; a larger value indicates a more significant hysteresis dissipation caused by defects.

[0173] The above three types of features are combined into the original defect discrimination feature vector: Where P is the dimension of the defect discrimination feature. Let i be the i-th original feature.

[0174] To eliminate differences in units and numerical ranges among different features, this embodiment of the invention performs Z-score standardization on all features: .

[0175] In the formula, and denoted as the mean and standard deviation of feature i on the benchmark training samples (healthy specimens and specimens with known defects), respectively.

[0176] After standardization, the dimensions of each feature are unified and eliminated, with a mean of 0 and a variance of 1, preventing features with large numerical ranges from dominating the classification model's decisions. The standardized defect discrimination feature vector is denoted as: .

[0177] In this embodiment of the invention, the probabilistic classification model may be, but is not limited to, a Gaussian process classification model, a support vector machine probabilistic model, or a probabilistic neural network model. The following explanation uses the Gaussian process classification model as an example.

[0178] In this embodiment of the invention, the Gaussian Process Classification (GPC) model uses a Radial Basis Function (RBF) function, defined as: .

[0179] In the formula, and Let be the i-th and j-th standardized defect discrimination feature vectors, respectively, and l be the length scaling parameter used to control the rate at which sample similarity decays with increasing feature distance. Let V be the signal variance.

[0180] Hyperparameters of Gaussian process classification model The posterior distribution of GPC is determined by maximizing the log-marginal likelihood function under the Laplace approximation. Since the posterior distribution of GPC typically lacks a closed-form analytical solution, the Laplace approximation is used for the solution. .

[0181] In the formula, is the maximum a posteriori estimate of the latent function; K is the kernel matrix; W is the second derivative matrix of the negative log-likelihood with respect to the latent function value.

[0182] The conjugate gradient method or quasi-Newton method (such as BFGS) is used to... Regarding the maximization of hyperparameters to obtain the optimal hyperparameters .

[0183] The Gaussian model is trained in advance, and the training process can be understood as follows: Step 1: Constructing the training set This invention employs supervised learning to train a probabilistic classification model. The training set contains two classes of samples: Healthy samples: collected from defect-free curved surface weld specimens under different excitation positions, different measuring point positions, and multiple repeated excitation conditions, and labeled as y=0; Defect samples: collected from curved weld specimens with different crack lengths (e.g., 5mm, 10mm, 15mm, 20mm, 25mm) and crack depths, and labeled as y=1.

[0184] The total number of training samples shall be no less than 60 sets, of which no less than 20 sets are healthy samples and no less than 40 sets are defective samples (at least a number of repeated response signals shall be collected for each defect level) to ensure that the model learns the features of each category of samples fully and avoids underfitting or overfitting due to insufficient samples.

[0185] Step 2: Divide the dataset into training and test sets. In this embodiment of the invention, all samples are divided into a training set and a test set using stratified random sampling: Training set: Used for learning model parameters and hyperparameters, accounting for 70%~80%; Test set: Used to evaluate the generalization performance of the model, accounting for 30%~20%.

[0186] Stratified random sampling can ensure that the proportion of healthy samples and defective samples in the training set and test set is basically consistent with the overall distribution (i.e., the training set and test set are divided according to the same healthy / defective ratio).

[0187] Step 3: Model Training and Testing The Gaussian process classification model is trained using the training set. After training, the classification accuracy, recall, specificity, and area under the receiver operating characteristic (ROC) curve (AUC) are calculated using the test set to verify the effectiveness of the Gaussian process classification model.

[0188] After training, the standardized defect discrimination feature vector corresponding to the welded surface component under test is... Input a Gaussian process classification model, and the model outputs the posterior probability that the sample belongs to the defect category: .

[0189] In the formula, The feature matrix of the training set; The training set label vector; . The value reflects the confidence level that the sample under test has defects. The closer the value is to 1, the higher the confidence level of the model in the tested curved welded component being in a defective state; the closer the value is to 0, the higher the confidence level of the model in the tested curved welded component being in a healthy state.

[0190] In one feasible implementation, embodiments of the present invention obtain the probability of defect presence. Then, the state can be determined based on, but is not limited to, a preset probability threshold.

[0191] Specifically, in this embodiment of the invention, the state of the weld on the curved surface of the tested pipeline can be divided into three states, and the judgment criteria are as follows: When P defect ≥P H At this point: the tested curved surface weld component is determined to be in a defective state. The model then concludes with a high confidence level that a closed crack exists in the weld, recommending repair or replacement.

[0192] When P defect ≤P L At this point, the tested curved surface weld component is determined to be in a healthy state. The model assumes the weld structure is intact with high confidence and requires no further processing.

[0193] When P L <P defect <P HIf the measured curved surface weld component is in an uncertain state, the model confidence level is insufficient, prompting the testing personnel to conduct supplementary testing (such as adjusting the excitation position, adding measuring points, or using other non-destructive testing methods for verification).

[0194] In this embodiment of the invention, the upper probability threshold P H A value of 0.95 can be used, with a lower probability threshold P. L A value of 0.05 is acceptable. It should be noted that the above threshold can be adjusted based on the training sample size, false positive cost, false negative cost, and engineering safety level, and is not limited here.

[0195] The final output includes the probability of defect presence and the state of the curved weld.

[0196] Based on the aforementioned method for detecting surface weld defects based on phase space topological constraints provided in the first aspect, embodiments of the present invention provide a device for detecting surface weld defects based on phase space topological constraints, such as... Figure 3 As shown, the surface weld defect detection device based on phase space topological constraints includes: The signal acquisition and preprocessing unit 110 is used to apply a single-point transient pulse excitation to the weld seam specimen of the pipe under test, simultaneously acquire the excitation force signal and the vibration response signal, and perform windowing truncation, high-pass filtering and amplitude normalization preprocessing on the acquired excitation force signal and the vibration response signal respectively to obtain the preprocessed excitation force signal and the preprocessed vibration response signal. The phase space reconstruction unit 120 is used to determine the delay time and embedding dimension based on the preprocessed vibration response signal, construct the phase space trajectory matrix of the measured vibration response signal, establish a recursive matrix based on the phase space trajectory matrix, and extract the measured recursive quantitative analysis index vector. The equation identification unit 130 is used to construct a candidate function library containing polynomial terms, hysteresis basis function terms, differential operator terms and cross-coupling terms. Using the preprocessed excitation force signal as input, it performs sparse regression on multiple candidate sparse regularization parameters to obtain multiple candidate hysteresis dynamic equations. The topology selection unit 140 is used to generate predicted vibration response signals according to each of the candidate hysteresis dynamic equations, calculate the time-domain fitting error between the predicted vibration response signal and the preprocessed vibration response signal, and reconstruct the phase space of the predicted vibration response signal again and extract the predicted recursive quantization analysis index vector, calculate the phase space topology deviation between the predicted recursive quantization analysis index vector and the measured recursive quantization analysis index vector; and determine the final sparse regularization parameter and the final hysteresis dynamic equation based on the minimum phase space topology deviation among the candidate hysteresis dynamic equations whose time-domain fitting error is less than or equal to a preset error threshold. The defect discrimination unit 150 is used to construct a defect discrimination feature vector based on the measured recursive quantitative analysis index vector and the structural features and equation coefficient features of the final hysteresis dynamic equation, and to obtain the probability of defect existence through a pre-trained probability classification model.

[0197] Based on the aforementioned method for detecting surface weld defects based on phase space topological constraints provided in the first aspect, this embodiment of the invention also provides a storage medium storing a computer-executable program. The computer-executable program is used to cause a computer to execute the method for detecting surface weld defects based on phase space topological constraints as described in any implementation of the first aspect. Explanations of the relevant content and descriptions of the beneficial effects of any of the computer-readable storage media provided above can be found in the corresponding embodiments described above, and will not be repeated here.

[0198] Those skilled in the art will understand that the program for implementing all or part of the steps of the above embodiments, which can be executed by a program instructing related hardware, can be stored in a computer-readable storage medium. The storage medium mentioned above can be a read-only memory, a random access memory, etc. The processing unit or processor mentioned above can be a central processing unit, a general-purpose processor, an application-specific integrated circuit (ASIC), a microprocessor (DSP), a field-programmable gate array (FPGA), or other programmable logic devices, transistor logic devices, hardware components, or any combination thereof.

[0199] This invention also provides a computer program product containing instructions that, when executed on a computer, cause the computer to perform any of the methods described in the above embodiments. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the flow or function according to the embodiments of this invention is generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, computer instructions can be transmitted from one website, computer, server, or data center to another via wired (e.g., coaxial cable, fiber optic, digital subscriber line (DSL)) or wireless (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium accessible to a computer or a data storage device such as a server or data center that integrates one or more available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., DVD), or a semiconductor medium (e.g., SSD), etc.

[0200] It should be noted that the devices for storing computer instructions or computer programs provided in the embodiments of the present invention, such as, but not limited to, the aforementioned memory, computer-readable storage medium, and communication chip, are all non-transitory. Those skilled in the art should recognize that the functions described in the embodiments of the present invention in one or more of the above examples can be implemented using hardware, software, firmware, or any combination thereof. When implemented using software, these functions can be stored in a computer-readable storage medium or transmitted as one or more instructions or code on a computer-readable storage medium. Computer-readable storage media include computer storage media and communication media, wherein communication media include any medium that facilitates the transmission of computer programs from one place to another. Storage media can be any available medium accessible to general-purpose or special-purpose computers.

[0201] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

Claims

1. A method for detecting defects in curved surface welds based on phase space topological constraints, characterized in that, include: A single-point transient pulse excitation was applied to the weld specimen of the curved surface of the pipeline under test, and the excitation force signal and vibration response signal were acquired simultaneously. The acquired excitation force signal and vibration response signal were preprocessed by windowing, high-pass filtering and amplitude normalization to obtain preprocessed excitation force signal and preprocessed vibration response signal. The delay time and embedding dimension are determined based on the preprocessed vibration response signal, and the phase space trajectory matrix of the measured vibration response signal is constructed. A recursive matrix is ​​established based on the phase space trajectory matrix, and the measured recursive quantitative analysis index vector is extracted. A candidate function library containing polynomial terms, hysteresis basis function terms, differential operator terms, and cross-coupling terms is constructed. Using the preprocessed excitation force signal as input, sparse regression is performed on multiple candidate sparse regularization parameters to obtain multiple candidate hysteresis dynamic equations. Based on each candidate hysteresis dynamic equation, a predicted vibration response signal is generated, and the time-domain fitting error between the predicted vibration response signal and the preprocessed vibration response signal is calculated. The predicted vibration response signal is then reconstructed in phase space, and a predicted recursive quantization analysis index vector is extracted. The phase space topology deviation between the predicted recursive quantization analysis index vector and the measured recursive quantization analysis index vector is calculated. Among the candidate hysteresis dynamic equations whose time-domain fitting error is less than or equal to a preset error threshold, the final sparse regularization parameter and the final hysteresis dynamic equation are determined based on the criterion of minimizing the phase space topological deviation. Based on the measured recursive quantitative analysis index vector and the structural and coefficient features of the final hysteresis dynamic equation, a defect discrimination feature vector is constructed, and the probability of defect existence is obtained through a pre-trained probability classification model.

2. The method for detecting defects in curved surface welds based on phase space topological constraints according to claim 1, characterized in that, The delay time is determined by the first local minimum of the average mutual information function; the embedding dimension is determined by the dimension corresponding to the first time the proportion of false nearest neighbors is less than a preset threshold; the neighborhood radius of the recursion matrix is ​​determined according to the quantile of the distance between phase space trajectory points, or determined in reverse according to the target recursion rate; the recursive quantification analysis index vector includes at least one or more of the following: recursion rate, determinism, average diagonal length, maximum diagonal length, divergence, diagonal length distribution entropy, laminar flow rate, and dwell time.

3. The method for detecting defects in curved surface welds based on phase space topological constraints according to claim 1, characterized in that, The highest order of the polynomial term does not exceed 5; the lag basis function term includes Bouc-Wen type lag basis function form and Duhem type lag basis function form; the cross-coupling term is the product term of the polynomial term and the lag basis function term.

4. The method for detecting defects in curved surface welds based on phase space topological constraints according to claim 1, characterized in that, The process involves using the preprocessed excitation force signal as input and performing sparse regression on multiple candidate sparse regularization parameters to obtain multiple candidate hysteresis dynamic equations, including: The NSGA-III multi-objective optimization algorithm is used to filter the candidate function items in the candidate function library; The chromosome of the NSGA-III multi-objective optimization algorithm adopts a binary encoding method, and each gene bit indicates whether the corresponding candidate function item in the candidate function library is selected into the equation structure; The initial population is generated by a combination of random initialization and prior initialization. The random initialization sets each candidate function term to 1 with a preset probability, and the prior initialization sets the constant term, linear term, excitation input term and at least one of the hysteresis basis function terms to 1. The crossover operation uses either uniform crossover or single-point crossover, with a crossover probability of 0.

9. The mutation operation uses bit-flip mutation, with a mutation probability of 1 / d to 0.1, where d is the number of candidate function terms; With minimizing temporal fitting error and phase space topological deviation as dual optimization objectives, non-dominated sorting and crowding distance calculation are performed, and the next generation population is selected according to non-dominated level and crowding distance. When the termination condition is met, output the Pareto front solution set; The candidate function terms corresponding to each solution in the Pareto front solution set are combined as the structure of the candidate hysteresis dynamic equation, and sparse regression is performed on the structure of the candidate hysteresis dynamic equation to obtain the candidate hysteresis dynamic equation.

5. The method for detecting defects in curved surface welds based on phase space topological constraints according to claim 4, characterized in that, The sparse regression operation on the structure of the candidate hysteresis dynamics equations includes: Sparse regression is solved using the LASSO algorithm or the sequential threshold least squares algorithm. When using the sequential threshold least squares algorithm, the least squares solution is used as the initial coefficients. Threshold truncation and partial least squares reconstruction are performed iteratively until the set of activated candidate options no longer changes or the preset maximum number of iterations is reached. The sparse truncation threshold τ... s The relationship between the candidate sparse regularization parameter λ and the parameter λ satisfies , This is the initial least squares coefficient matrix. This represents the maximum absolute value of the initial coefficient, and α is the threshold proportionality coefficient.

6. The method for detecting defects in curved surface welds based on phase space topological constraints according to claim 1, characterized in that, The step of generating predicted vibration response signals based on each of the candidate hysteresis dynamic equations includes: Using the initial state of the preprocessed vibration response signal as the initial condition, the fourth-order Runge-Kutta method or a variable-step ordinary differential equation solver is used to numerically integrate each of the candidate hysteresis dynamic equations to generate the predicted vibration response time series corresponding to each candidate hysteresis dynamic equation.

7. The method for detecting defects in curved surface welds based on phase space topological constraints according to claim 1, characterized in that, The defect discrimination feature vector includes one or more of the following: recursive quantitative analysis index, the number of activated terms in the final hysteresis dynamic equation, the proportion of nonlinear terms, the proportion of hysteresis terms, the norm of the equation coefficients, and the equivalent hysteresis dissipation index.

8. The method for detecting defects in curved surface welds based on phase space topological constraints according to claim 1, characterized in that, The probabilistic classification model is a Gaussian process classification model, a support vector machine probabilistic model, or a probabilistic neural network model. When the probabilistic classification model is a Gaussian process classification model, a radial basis function kernel function is used; The hyperparameters of the Gaussian process classification model are determined by maximizing the log-marginal likelihood function under the Laplace approximation.

9. The method for detecting defects in curved surface welds based on phase space topological constraints according to claim 1, characterized in that, After obtaining the probability of defect presence through the pre-trained probabilistic classification model, the method further includes: The state of the curved weld is determined based on the comparison between the probability of the defect and the preset probability threshold. Output the probability of the defect's existence and the state of the curved weld.

10. A device for detecting defects in curved surface welds based on phase space topological constraints, characterized in that, include: The signal acquisition and preprocessing unit is used to apply single-point transient pulse excitation to the weld specimen of the curved surface of the pipe under test, simultaneously acquire the excitation force signal and the vibration response signal, and perform windowing truncation, high-pass filtering and amplitude normalization preprocessing on the acquired excitation force signal and the vibration response signal respectively to obtain the preprocessed excitation force signal and the preprocessed vibration response signal. The phase space reconstruction unit is used to determine the delay time and embedding dimension based on the preprocessed vibration response signal, construct the phase space trajectory matrix of the measured vibration response signal, establish a recursive matrix based on the phase space trajectory matrix, and extract the measured recursive quantitative analysis index vector. The equation identification unit is used to construct a candidate function library containing polynomial terms, hysteresis basis function terms, differential operator terms and cross-coupling terms. Using the preprocessed excitation force signal as input, it performs sparse regression on multiple candidate sparse regularization parameters to obtain multiple candidate hysteresis dynamic equations. The topology selection unit is used to generate predicted vibration response signals based on each of the candidate hysteresis dynamic equations, calculate the time-domain fitting error between the predicted vibration response signal and the preprocessed vibration response signal, reconstruct the phase space of the predicted vibration response signal again, extract the predicted recursive quantization analysis index vector, calculate the phase space topology deviation between the predicted recursive quantization analysis index vector and the measured recursive quantization analysis index vector, and determine the final sparse regularization parameter and the final hysteresis dynamic equation based on the minimum phase space topology deviation among the candidate hysteresis dynamic equations whose time-domain fitting error is less than or equal to a preset error threshold. The defect discrimination unit is used to construct a defect discrimination feature vector based on the measured recursive quantitative analysis index vector and the structural features and equation coefficient features of the final hysteresis dynamic equation, and to obtain the probability of defect existence through a pre-trained probability classification model.

Images