A method and system for terahertz signal aliasing elimination based on deconvolution

Abstract

The application provides a terahertz signal aliasing elimination method and system based on deconvolution, and relates to the technical field of detection. The method comprises the following steps: collecting time-domain waveforms of a reference area and a detection area, and obtaining geometric and material parameters; calculating a point spread function based on the reference waveform and geometric priori; calculating a time axis normalization mapping function based on the geometric priori, and performing normalization processing on the observation waveform; performing deconvolution processing on the normalized time axis to obtain a recovered pulse train; generating a time-of-flight map based on the recovered pulse train, and performing curvature correction combined with the geometric priori to obtain a defect distribution map. By introducing the reference area priori and time normalization, the application effectively overcomes the point spread function drift caused by the curvature and the refractive index gradient, suppresses the noise amplification in the high-absorption medium, and can stably remove the signal aliasing and realize the accurate positioning of defects in the polyethylene hot melt joint scene.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of detection, in particular to a terahertz signal aliasing elimination method and system based on deconvolution. BACKGROUND

[0002] Terahertz time-domain spectroscopy and imaging are widely used in the nondestructive testing of pipe fusion joints due to their strong penetration of non-metallic materials such as polyethylene. During the detection process, the terahertz pulse propagates in the multi-layer polyethylene fusion joint and experiences multiple reflections and scattering, forming overlapping echo signals. In order to weaken the overlapping effect, existing researches mostly use window function truncation, frequency domain filtering, wave packet decomposition and other general signal processing methods. These methods can separate part of the overlapping information when detecting flat plates or parallel layered structures, but their basic assumption is that the system response remains stable in space. However, the fusion joint has curvature and bead geometry characteristics, accompanied by a gradual change in the radial refractive index, making the point spread function drift significantly with position, and it is difficult for traditional frequency domain or time domain processing methods to establish a physically consistent model.

[0003] The existing technology still has significant limitations when facing curved fusion joints and strongly absorbing polyethylene media: on the one hand, there is a lack of signal representation model that can adaptively update with geometric changes; on the other hand, under strong absorption conditions, conventional methods are prone to introduce noise or geometric distortion while eliminating aliasing, making it difficult to obtain a defect distribution corresponding to the real structure.

[0004] Therefore, how to overcome the adverse effects of point spread function drift and time axis nonlinearity in curved fusion joints and strongly absorbing polyethylene media, and robustly eliminate echo aliasing while obtaining a defect distribution that matches the geometric structure, has become a technical problem that needs to be solved. SUMMARY

[0005] In view of the deficiencies of the prior art, the present application provides a terahertz signal aliasing elimination method and system based on deconvolution.

[0006] In a first aspect, the present application provides a terahertz signal aliasing elimination method based on deconvolution, comprising:

[0007] Collecting the time-domain waveforms of the reference region and the detection region to obtain a first reference waveform and a first observation waveform; obtaining geometric and material parameters to obtain a first geometric prior parameter;

[0008] Calculating a point spread function based on the first reference waveform and the first geometric prior parameter to obtain a first point spread function;

[0009] Calculating a time axis normalization mapping function based on the first geometric prior parameter to obtain a first time normalization mapping function;

[0010] applying the first time normalization mapping function to the first observed waveform to obtain a first normalized observed waveform;

[0011] performing deconvolution processing based on the first normalized observed waveform to obtain a recovered pulse train;

[0012] generating a time-of-flight map based on the recovered pulse train and performing curvature correction based on the first geometric prior parameter to obtain a defect distribution map.

[0013] Optionally, the method further comprises:

[0014] performing deconvolution processing on the first observed waveform on an un-normalized original time axis to obtain a second recovered pulse train;

[0015] applying the first time normalization mapping function to the second recovered pulse train to obtain a second normalized recovered pulse train;

[0016] calculating an exchangeability difference value based on the recovered pulse train and the second normalized recovered pulse train to obtain a first coupling constraint value;

[0017] performing joint iterative updating based on the first coupling constraint value, a convolution residual between the first normalized observed waveform and the first point spread function, and a sparsity index of the recovered pulse train to obtain an updated first time normalization mapping function, an updated first point spread function, and an updated recovered pulse train.

[0018] Optionally, the method further comprises:

[0019] calculating a first stability index and a second stability index based on the updated recovered pulse train to obtain a first stop condition;

[0020] terminating the joint iterative updating when the first stop condition is satisfied and outputting a target recovered pulse train as the recovered pulse train.

[0021] Optionally, the performing deconvolution processing based on the first normalized observed waveform comprises:

[0022] acquiring normal information of a detection surface, determining a direction set based on the first geometric prior parameter and the normal information to obtain a first direction set;

[0023] calculating a direction-dependent point spread function set based on the first reference waveform and the first direction set to obtain a first direction kernel cluster;

[0024] performing direction filtering based on the first normalized observed waveform and the first direction set to obtain first direction up-dimensional data;

[0025] performing group sparse deconvolution processing on the first direction kernel cluster and the first direction elevated dimension data to obtain a direction indexed deconvolved pulse tensor;

[0026] performing direction neighborhood smoothing processing on the direction indexed deconvolved pulse tensor to obtain a direction continuous deconvolved pulse tensor;

[0027] projecting the direction continuous deconvolved pulse tensor to a time axis and a scanning axis according to the first direction set to obtain the deconvolved pulse train.

[0028] Optionally, the determining the direction set based on the first geometric prior parameter and the normal information comprises:

[0029] generating an initial direction set based on the first geometric prior parameter and the normal information to obtain a first initial direction set;

[0030] calculating a direction indicator based on the first normalized observation waveform within a preset time window and a preset scanning window to obtain a first direction indicator;

[0031] calculating a direction consistency score based on the first direction indicator and the first initial direction set to obtain a first direction consistency score;

[0032] performing contraction and resampling on the first initial direction set based on the first direction consistency score within an angle interval range defined based on the first geometric prior parameter to obtain a second direction set.

[0033] Optionally, the method further comprises:

[0034] calculating a set change quantity based on the first initial direction set and the second direction set to obtain a first direction stability condition;

[0035] outputting the second direction set as the first direction set when the first direction stability condition is satisfied, or repeating the generating of the direction set with the second direction set as a new initial direction set when the first direction stability condition is not satisfied until the first direction stability condition is satisfied.

[0036] Optionally, the performing group sparse deconvolution processing on the first direction kernel cluster and the first direction elevated dimension data to obtain a direction indexed deconvolved pulse tensor comprises:

[0037] determining a direction adjacency relationship based on the first geometric prior parameter and the normal information to generate a first direction compatible matrix;

[0038] In the group sparse solution, a shared constraint is imposed on the non-zero support of the direction dimension according to the first direction compatibility matrix, so that the non-zero directions corresponding to the same time sample are limited to the adjacent directions specified by the first direction compatibility matrix, and a recovered pulse tensor of the direction index is output.

[0039] Optionally, the first direction compatibility matrix is generated based on the first geometric prior parameter and the normal information.

[0040] The equivalent incidence angle corresponding to each direction is determined based on the first geometric prior parameter and the normal information, and the reflection coefficient phase characteristics corresponding to the equivalent incidence angle are evaluated in the terahertz working frequency band, and based on the in-band phase sign consistency, the direction pairs with mutually exclusive phases are marked to obtain a first incompatible set;

[0041] After excluding the first incompatible set, the main echo group time delay order of the adjacent scanning positions is calculated based on the first geometric prior parameter, and the direction pairs are screened based on the group time delay order as a criterion to obtain a first compatible set;

[0042] The direction pairs belonging to the first compatible set are set as allowed adjacency in the first direction compatibility matrix, and the direction pairs belonging to the first incompatible set are set as prohibited adjacency in the first direction compatibility matrix.

[0043] Optionally, the defect distribution map is obtained by:

[0044] The equivalent incidence angle of each scanning position is determined based on the first geometric prior parameter, and the representative frequency of each recovered pulse is determined based on the spectral energy center of the recovered pulse train to obtain a first representative frequency set;

[0045] The interface reflection phase label is calculated based on the equivalent incidence angle and the first representative frequency set to obtain a first phase label set, and the measured phase label is calculated based on the analytical signal phase of the recovered pulse train at the pulse peak to obtain a second phase label set;

[0046] The phase consistency verification is performed based on the first phase label set and the second phase label set, and when the consistency is inconsistent, the polarity correction and the micro-time delay correction are performed on the corresponding recovered pulse to obtain a phase-aligned recovered pulse train;

[0047] The time-of-flight map is generated based on the phase-aligned recovered pulse train, and the angle-resolved group refractive index is determined based on the equivalent incidence angle and the first representative frequency set, and the curvature correction is performed to obtain the defect distribution map.

[0048] In a second aspect, the present application provides a terahertz signal aliasing elimination system based on deconvolution, comprising:

[0049] The acquisition module is configured to acquire time-domain waveforms of the reference region and the detection region to obtain a first reference waveform and a first observation waveform; and obtain geometric and material parameters to obtain a first geometric prior parameter;

[0050] The processing module is configured to calculate a point spread function based on the first reference waveform and the first geometric prior parameter to obtain a first point spread function; calculate a time-axis normalization mapping function based on the first geometric prior parameter to obtain a first time normalization mapping function; and apply the first time normalization mapping function to the first observation waveform to obtain a first normalized observation waveform;

[0051] The recovery module is configured to perform deconvolution processing based on the first normalized observation waveform to obtain a recovered pulse train;

[0052] The output module is configured to generate a time-of-flight map based on the recovered pulse train, and perform curvature correction based on the first geometric prior parameter to obtain a defect distribution map

[0053] Compared with the prior art, the present application obtains waveforms in the reference region and calculates a point spread function in combination with geometric and material parameters, and then introduces a time-axis normalization mapping, so that the system response based on deconvolution can be consistent with the actual propagation characteristics under the complex conditions of a curved hot melt joint and a radially refractive index gradient. This scheme avoids the model mismatch problem caused by the drift of the point spread function with position in the existing method, effectively improves the stability of echo separation, and reduces the artifacts caused by signal overlap. At the same time, the introduction of the time-axis normalization converts the originally non-stationary convolution process into an approximately stationary processing process, reducing the risk of noise amplification and providing a robust foundation for signal processing in high-absorption polyethylene materials.

[0054] Further, after obtaining the recovered pulse train, the present application performs curvature correction in combination with the geometric prior parameter, accurately maps the signal processing result to a time-of-flight map, and generates a defect distribution map. Through this process, the consistency of time-domain data and spatial structure can be maintained in a complex curved structure, and reliable characterization of defect position and morphology can be achieved. Compared with the existing general processing method based on a flat plate or a simplified layered model, the present application can output detection results consistent with the real geometry in the actual hot melt joint scene, significantly improve the defect positioning accuracy and imaging reliability, and has outstanding industrial application value. BRIEF DESCRIPTION OF DRAWINGS

[0055] Figure 1 A flowchart of a terahertz signal aliasing elimination method based on deconvolution provided by an embodiment of the present application;

[0056] Figure 2 A flowchart of a joint updating method provided by an embodiment of the present application;

[0057] Figure 3 A flowchart of a method for performing deconvolution processing is provided for the embodiments of the present application.

[0058] Figure 4 A schematic diagram of a terahertz signal aliasing elimination system based on deconvolution is provided for the embodiments of the present application. DETAILED DESCRIPTION

[0059] The technical solutions in the embodiments of the present application will be clearly and completely described with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application, not all the embodiments.

[0060] Referring to Figure 1 A flowchart of a method for performing deconvolution processing is provided for the embodiments of the present application.

[0061] S101: Collecting time-domain waveforms of a reference region and a detection region to obtain a first reference waveform and a first observation waveform; obtaining geometric and material parameters to obtain a first geometric prior parameter;

[0062] S102: Calculating a point spread function based on the first reference waveform and the first geometric prior parameter to obtain a first point spread function;

[0063] S103: Calculating a time-axis normalization mapping function based on the first geometric prior parameter to obtain a first time normalization mapping function;

[0064] S104: Applying the first time normalization mapping function to the first observation waveform to obtain a first normalized observation waveform;

[0065] S105: Performing deconvolution processing based on the first normalized observation waveform to obtain a recovered pulse train;

[0066] S106: Generating a time-of-flight map based on the recovered pulse train, and performing curvature correction based on the first geometric prior parameter to obtain a defect distribution map.

[0067] For the above S101:

[0068] In the present embodiment, a terahertz time-domain spectroscopy detection system is used to detect a polyethylene pipe hot melt joint. First, a reference signal is collected at a uniform defect-free area on the outer wall of the joint to obtain a reference time-domain waveform of the system under stable conditions, denoted as a first reference waveform. At the same time, a corresponding time-domain waveform is collected at the detection region of the joint to be tested, denoted as a first observation waveform. The above waveforms are obtained by a terahertz pulse source emitting and a detector receiving, and are output in the form of time series data.

[0069] On this basis, the geometric and material parameters for describing the physical characteristics of the joint are obtained as the first geometric prior parameters. The geometric parameters can include the curvature radius of the bead, the pipe wall thickness, the pipe diameter size, and the probe path angle, etc. The material parameters can include the refractive index, the dispersion characteristics, and the absorption coefficient of the polyethylene, etc. These parameters can be obtained by external contour measurement, mechanical measurement, or a combination of material database table lookup and actual measurement. The obtained geometric and material parameters are used as inputs for subsequent point spread function modeling and time axis normalization, ensuring that the signal processing model conforms to the actual joint structure and propagation environment.

[0070] For the above S102:

[0071] In an embodiment, the calculation of the point spread function is based on the first reference waveform obtained in the reference area and the first geometric prior parameters described above. Specifically, the time-domain waveform of the reference area is first used to represent the original response signal of the system under uniform and defect-free conditions. Since the reference area does not contain structural defects, the echo signal thereof can accurately reflect the transmission and reception characteristics of the detection system itself.

[0072] Subsequently, the first geometric prior parameters are introduced into the modeling and correction process of the reference waveform. The geometric prior parameters include the curvature radius of the bead, the pipe diameter size, the wall thickness information, and the incident path angle, etc. The material parameters include the refractive index curve, the dispersion relationship, and the absorption coefficient of the polyethylene in the working frequency band. The optical path length of different propagation paths is determined by the geometric parameters, and is converted into a time delay to correct the time delay distribution of the reference waveform. The phase unwrapping is performed on the frequency domain signal by the material parameters to compensate for the group delay distortion caused by dispersion, and an exponential type attenuation correction is applied to the amplitude component to eliminate the energy loss caused by medium absorption.

[0073] After the above correction is completed, the obtained response waveform is normalized and a minimum phase constraint is introduced, so that the generated system impulse satisfies the causality and energy concentration requirements. The waveform obtained through this process is the first point spread function, which can reflect the system response characteristics of the hot melt joint under actual geometric and material conditions, and is used as the operation kernel for the subsequent deconvolution step.

[0074] For the above S103:

[0075] In an embodiment, the calculation of the time axis normalization mapping function takes the first geometric prior parameter as input. Due to the characteristics of the hot melt joint having a bead curvature and a radial refractive index gradient, the propagation path length and the group delay at different positions on the same scan line are no longer in a linear relationship, and directly using the original time axis will introduce position-dependent distortion. Therefore, the embodiment establishes a geometric model, combines the curvature radius, wall thickness and incidence angle information in the geometric prior parameter, and calculates the equivalent optical path corresponding to each sampling point. The equivalent optical path is converted into group delay through the speed of light and the refractive index, thereby obtaining the time delay distribution of each sampling point.

[0076] On the basis of the above, the time delay distribution is compared with the ideal linear time axis to generate a one-to-one mapping relationship. The mapping relationship can re-parameterize the non-uniform sampling points of the original time axis into uniform sampling points on the normalized time axis, ensuring that the waveforms at different positions can be processed in a unified time coordinate system during subsequent deconvolution processing.

[0077] Specifically, the original sampling time can be mapped to the standard sampling point in the normalized domain, and interpolation or resampling methods can be used to realize signal reconstruction. The normalized mapping function generated in this way is defined as the first time normalization mapping function, which eliminates the time axis nonlinearity caused by curvature and refractive index gradient, so that the deconvolution process can be performed in the condition of approximately stationary time domain.

[0078] For S104 described above:

[0079] In an embodiment, the first observed waveform obtained by the detection area is mapped to the normalized time axis to obtain a first normalized observed waveform.

[0080] First, the aforementioned first time normalization mapping function is called to establish a one-to-one mapping table between the original time coordinate and the normalized time coordinate. The mapping table takes the sampling point serial number or the original sampling time as the index and outputs the target time after normalization. Preferably, it can be stored as an array or a matrix in a lookup table manner for subsequent fast calling.

[0081] Since the original sampling points often fall in non-uniform positions in the normalized domain, numerical interpolation methods are needed to avoid distortion. Specifically, linear interpolation, spline interpolation or resampling algorithm based on Fourier interpolation can be used. In this embodiment, cubic spline interpolation is preferred to smooth the reconstructed waveform. With commonly used numerical operation tools, such as MATLAB, SciPy.interpolate in Python or LabVIEW interpolation library, the reconstruction of sampling points can be realized.

[0082] After the interpolation is completed, the waveform amplitude obtained by resampling is normalized to ensure that the waveform energy does not shift with the change of coordinates. Subsequently, the cross-correlation value of the normalized waveform and the original waveform at the main pulse position is calculated as a consistency check index; when the cross-correlation coefficient is greater than a preset threshold, it is considered that the normalized mapping is correct.

[0083] Through the above steps, the original observation waveform is resampled to the normalized time axis to obtain a first normalized observation waveform. The waveform is uniformly distributed in the time domain, and the main pulse position is consistent with the geometrically corrected optical path, providing a unified coordinate reference for subsequent deconvolution.

[0084] For the above S105:

[0085] In an embodiment, a deconvolution process is performed on the first normalized observation waveform to obtain a corresponding recovered pulse train. Specifically, the first point spread function described above is first loaded into the deconvolution operation module as a system convolution kernel. The first point spread function has been established in combination with the reference region waveform, the bead curvature radius, the pipe wall thickness, and the material refractive index characteristics, and can reflect the propagation characteristics of the detection system under the geometric conditions of the hot melt joint.

[0086] Subsequently, the first normalized observation waveform is input into the deconvolution operation module. During the operation, the input signal is first subjected to frequency domain transformation to obtain a complex frequency spectrum representation. For example, the time domain waveform can be converted into a frequency domain spectrum line using fast Fourier transform (FFT), and then subjected to constrained inversion using the frequency domain transfer function of the first point spread function to preliminarily recover the time domain pulse train.

[0087] Considering that the high absorption characteristics of polyethylene material are prone to cause noise amplification during inversion, the present embodiment introduces sparsity constraint and convolution residual constraint during inversion. The sparsity constraint is used to limit the recovered pulse train to exhibit an isolated pulse structure in the time domain, so as to conform to the physical characteristics of echo events. For example, L1 norm regularization can be applied to the recovered result, or soft thresholding operation can be performed during iteration to suppress low-amplitude pseudo-pulses. The convolution residual constraint is used to ensure that the difference between the recovered pulse train and the first normalized observation waveform after convolution by the first point spread function is kept within a preset range. For example, the mean square error of the residual can be set to be less than 5% of the energy of the original waveform as a qualified condition.

[0088] On the solving method, an iterative inversion strategy can be adopted. Exemplarily, an iterative shrinkage threshold algorithm (ISTA) or an alternating direction method of multipliers (ADMM) can be selected as an optimization solving tool. In each iteration, the convolution residual is calculated based on the current recovered pulse train, and a sparse operator is applied to the pulse train to update the new estimated result. If the residual change amplitude is lower than the set threshold in continuous several iterations, or the number of iterations reaches the upper limit, for example, 100 times, the calculation is terminated, and the final recovered pulse train is output.

[0089] The obtained recovered pulse train exhibits a series of independent pulses after dealiasing in the time domain, and the pulse position and amplitude can correspond to the actual path of multiple reflections in the hot melt joint. For example, when there is an unfused defect in the detection area, an additional pulse with a shorter delay than the normal wall thickness will appear in the recovered pulse train, thereby providing basic data for the generation and curvature correction of the subsequent time-of-flight map.

[0090] For the above S106:

[0091] First, each pulse event in the recovered pulse train is mapped to time-of-flight data. Specifically, the peak position of each pulse is extracted as the arrival time of the echo, and converted to time-of-flight with a sampling interval. In order to facilitate two-dimensional imaging, multiple time-of-flight data obtained from the same scan line are arranged in order of scan position to form a time-of-flight profile. By combining the profiles of all scan lines, a time-of-flight map is constructed. Exemplarily, gray value or pseudo-color mapping can be used to represent the size of time-of-flight, so that echoes at different depth positions can be intuitively presented in the image.

[0092] Subsequently, the time-of-flight map is corrected for curvature in combination with the first geometric prior parameter. Since the hot melt joint has a curved surface geometry, directly using the uncorrected time-of-flight map will cause the same physical depth to appear as a curved trajectory in the image, affecting the accurate interpretation of the defect position. Therefore, according to the bead curvature radius and wall thickness information, the equivalent propagation path length of each scan position is calculated, and the original time-of-flight data is remapped to the geometric corrected coordinate system. Exemplarily, by establishing a polar coordinate model, the time-of-flight data under the influence of curvature can be unfolded into an isobath layer in the Cartesian coordinate system, so that the position of the defect in the image is consistent with the actual geometric position.

[0093] The time-of-flight image after curvature correction is defined as a defect distribution map, which can truly reflect the defect condition of the detection area in space. For example, when there is an unfused defect, a local abnormal depth value or an energy void will appear in the corresponding area in the corrected defect distribution map; when there is a gas hole, it will appear as a point or patch-shaped abnormal signal in the image. In this way, the embodiment not only completes the signal domain deconvolution, but also accurately projects the processing result into the actual geometric structure, providing a reliable basis for defect detection and positioning.

[0094] Optionally referring to Figure 2 For the flowchart of the joint updating method provided by the embodiment of the application, in another implementation, in order to verify and constrain the non-commutativity that may be generated by the "time axis normalization processing" and the "deconvolution processing" in the curved surface and refractive index gradient scene, the embodiment introduces a double-path consistency verification and joint updating mechanism on the basis of the flowchart shown in Figure 1 , and specifically includes steps S201-S204, wherein:

[0095] S201: performing deconvolution processing on the first observation waveform on the original time axis without normalization to obtain a second recovered pulse train;

[0096] S202: applying the first time normalization mapping function to the second recovered pulse train to obtain a second normalized recovered pulse train;

[0097] S203: calculating a commutativity difference value based on the recovered pulse train and the second normalized recovered pulse train to obtain a first coupling constraint value;

[0098] S204: performing joint iterative updating based on the first coupling constraint value, a convolution residual between the first normalized observation waveform and the first point spread function, and a sparsity index of the recovered pulse train to obtain an updated first time normalization mapping function, an updated first point spread function, and an updated recovered pulse train.

[0099] For the above S201:

[0100] Under the condition that the time axis is not normalized, the first observation waveform is deconvolved once to obtain a second recovered pulse train. Specifically, the first point spread function is called as a convolution kernel to establish an inversion model with "sparsity constraint + convolution residual constraint" as the target on the original time axis; for example, the iterative shrinkage threshold method or the alternating direction multiplier method can be used for solving, and the convolution residual is calculated and the sparsification operator is applied to suppress low-amplitude pseudo-pulses in each iteration, while keeping the causality and minimum phase constraint of the first point spread function unchanged.

[0101] Preferably, band-pass filtering and amplitude normalization are performed on the first observed waveform before entering the solver, so that the input and the convolution kernel are in the same bandwidth and dimension. The result of this step preserves the time delay details before the curvature is unfolded in the original time coordinate system, as a contrast benchmark for the subsequent path consistency check.

[0102] For the above S202:

[0103] The second recovered pulse train obtained in S201 is input into the aforementioned first time normalization mapping function, mapped to the normalized time axis, to obtain a second normalized recovered pulse train. To avoid distortion caused by sampling misalignment, spline interpolation resampling is preferably used; to eliminate overall amplitude and micro time delay deviations, the main peak window of the two pulse trains can be exemplarily subjected to amplitude normalization and cross-correlation alignment (to correct the entire segment offset in a sample-level micro-shift manner). The mapping table can adopt a lookup table structure (array or matrix) for fast calling; interpolation and cross-correlation can be realized through common numerical libraries such as MATLAB, Python (SciPy.interpolate, NumPy.fft), or LabVIEW.

[0104] For the above S203:

[0105] On the unified normalized time axis and the same scanning position index, the recovered pulse train and the second normalized recovered pulse train obtained in S202 are subjected to difference measurement, and a first coupling constraint value is formed accordingly. To enhance physical pertinence, the difference measurement can include the following complementary indicators:

[0106] Amplitude envelope difference: the envelope of the two columns of signals is subjected to mean square difference or correlation coefficient measurement within a preset window, for describing energy distribution difference;

[0107] Phase / polarity consistency: phase flipping and polarity mismatch are distinguished based on the analytical signal phase or sign consistency; for example, in the high absorption segment of the terahertz operating frequency band, the phase term weight can be increased to reduce the interference of amplitude attenuation on the measurement;

[0108] Morphology consistency: the statistics of peak number, peak width, and peak spacing are compared to avoid false consistency caused by “peak shape splitting / merging”.

[0109] The above multi-dimensional differences are aggregated into a single scalar according to preset weights, defined as the commutativity difference value; for subsequent positioning-type updating, the difference value and its corresponding scanning position, time window index are archived together, forming an indexed representation of the first coupling constraint value. Exemplarily, the weights can be offline set by a verification data set, and the window can take a local interval of the main peak ± several sampling points.

[0110] For the above S204:

[0111] With the first coupling constraint value, the convolution residual and the sparsity index as the joint target, the time normalization mapping function, the point spread function and the recovered pulse train are updated alternately, and the updated first time normalization mapping function, the updated first point spread function and the updated recovered pulse train are output. To ensure physical feasibility, the embodiment imposes the following boundary conditions on each update object:

[0112] For the time normalization mapping function, only small adjustments are allowed at key nodes, and the mapping function is kept monotonic and smooth, for example, a monotonic spline or piecewise linear function is used, and the node displacement is limited to no more than a preset pixel per round;

[0113] For the point spread function, only the amplitude attenuation parameter and the phase are allowed to be fine-tuned, and the causality and minimum phase constraints are continuously maintained to prevent the generation of non-physical long tails or pre-ringing;

[0114] For the recovered pulse train, continue to use the iterative solution with sparsity and residual double constraints, and use the decrease of commutativity difference as an additional convergence signal.

[0115] In a specific implementation, the update order of "pulse train-mapping-point spread function" can be used, or the order of "pulse train-point spread function-mapping" can be used, as long as the joint decrease of the first coupling constraint value and the convolution residual is used as the iterative promotion condition at the end of each round. When the preset convergence criterion is met, the current round of joint update is terminated, and the updated three objects are obtained; the determination of the stop condition can be further described in subsequent embodiments.

[0116] For example, the upper limit of the iteration can be set to 50-200 rounds, the upper limit of the single mapping node displacement can be set to 0.1-0.5 sampling intervals, and the upper limit of the single fine-tuning of the point spread function phase can be set to a small angle interval to avoid excessive correction.

[0117] Through the double-path consistency and joint update mechanism of S201-S204, without changing the hardware acquisition process, the data-driven coupling constraint is used to dynamically correct the model mismatch of the time normalization mapping and the point spread function, and to suppress the systematic error caused by the geometric-convolution non-commutativity in the curved surface scene, thereby improving the consistency and interpretability of the recovered pulse train in the normalization domain.

[0118] Optionally, it further comprises:

[0119] Based on the updated recovered pulse train, a first stability index and a second stability index are calculated to obtain a first stop condition;

[0120] When the first stop condition is met, the joint iterative update is terminated, and the target recovered pulse train is output as the recovered pulse train.

[0121] In another embodiment, to avoid information exhaustion type "false convergence" or artifacts caused by over iteration in joint iterative update, the embodiment sets a stop criterion based on S201-S204 Figure 1 The process and S201-S204 are based on setting a stop criterion.

[0122] In a specific implementation, based on the updated recovered pulse train output by S204, combined with the updated time-normalized mapping function and point spread function of the corresponding round, the amplitude of the recovered pulse train is normalized and the micro-delay is aligned on the normalized time axis. For example, cross-correlation alignment can be performed in a local window of ± several sampling points of the main peak, and weighting is performed on the high absorption frequency band by sub-band or band-pass filtering, and a standardized data set for index calculation is output.

[0123] In the local window, the envelope gradient change rate and residual error energy drop are calculated for the recovered pulse train and its convolution residual error of the adjacent two rounds (round k and round k-1), and the change amount of the non-zero support set is counted, for example, using the set Hamming distance or intersection ratio.

[0124] For example, when the envelope gradient change rate is less than 1%-3% and the residual error energy drop is lower than the preset threshold, it is recorded as the first stability of the scanning position meeting the standard; aggregate all scanning positions to obtain the global score of the first stability index.

[0125] Under the same normalized time axis and scanning position index, the consistency of the updated recovered pulse train with the geometric prior is evaluated. Specifically, it can include: monotonicity test on the main echo group time delay order of adjacent scanning positions, and consistency test on the polarity and phase continuity of the recovered pulse, for example, by analyzing the signal phase or symbol consistency to determine whether the phase flip appears non-physical jump.

[0126] For example, the phase consistency term is given a higher weight in the high absorption frequency band to reduce the influence of amplitude attenuation on the judgment; when the monotonicity test passes and the phase / polarity consistency rate reaches the preset proportion, it is recorded as the second stability of the scanning position meeting the standard. Aggregate all scanning positions to obtain the global score of the second stability index.

[0127] The first stability index and the second stability index are jointly determined to form the first stop condition.

[0128] For example, the "double threshold + hysteresis" rule can be used: when the first stability index is higher than threshold A and the second stability index is higher than threshold B, and it is satisfied for K consecutive iterations, for example, K is 3-5, and the proportion of scanning positions meeting the standard is not less than the preset proportion, for example, not less than 90%, it is considered that the first stop condition is reached. To avoid jitter caused by local noise, local masks can be set for individual scanning positions that do not meet the standard and delayed to be included in global judgment.

[0129] When the first stop condition is met, the joint iterative update in S204 is terminated, the updated pulse train in the current round is output as the target pulse train, and the corresponding time-normalized mapping function and the point spread function are frozen for subsequent imaging; if the first stop condition is not met, return to S204 for the next round of joint update. Exemplarily, a minimum and a maximum iteration round can be set, for example, no less than 10 rounds and no more than 200 rounds, and when the maximum round is reached and the first stop condition is still not met, the principle of giving priority to the second stability index is used for final freezing to ensure that the physical consistency is not sacrificed.

[0130] In this way, the stopping time of the joint iteration is determined from the two complementary dimensions of numerical stability and physical consistency, which avoids "false convergence" caused by only residual convergence, and suppresses non-physical artifacts caused by over-iteration, thereby ensuring that the output target pulse train has interpretability and reproducibility in the curved and refractive index gradient terahertz detection scene.

[0131] Optionally, referring to Figure 3 A flowchart of a method for performing deconvolution processing provided by an embodiment of the present application includes steps S301-S306, wherein:

[0132] S301: Collect normal information of a detection surface, determine a direction set based on the first geometric prior parameter and the normal information, and obtain a first direction set;

[0133] S302: Calculate a direction-dependent point spread function set based on the first reference waveform and the first direction set, and obtain a first direction kernel cluster;

[0134] S303: Perform direction filtering based on the first normalized observation waveform and the first direction set, and obtain first direction up-dimensional data;

[0135] S304: Perform group sparse deconvolution processing based on the first direction up-dimensional data and the first direction kernel cluster, and obtain a direction-indexed pulse recovery tensor;

[0136] S305: Perform direction neighborhood smoothing processing on the direction-indexed pulse recovery tensor, and obtain a direction-continuous pulse recovery tensor;

[0137] S306: Project the direction-continuous pulse recovery tensor to the time axis and the scanning axis according to the first direction set, and obtain the pulse train.

[0138] First, in the description of the present embodiment, represents a scanning position index or a transverse coordinate along a scanning path, and is used to identify spatial sampling points within the same scanning line or between adjacent scanning lines; This represents the time sampling point on the normalized time axis (the time after reparameterization by the first time normalization mapping function), used to unify the time domain scale at different locations; The direction variable related to the local normal of the detected surface is defined as the angle between the incident / reflected direction and the normal. Represents the elements of the discretized direction set, { } is an ordered finite set obtained by discretizing the effective angular domain; Indication and direction The corresponding direction-related point diffusion function (direction kernel) is obtained by combining the first reference waveform with the first geometric prior parameters through direction consistency transformation and normalization. The tensor represents the direction of the up-dimensional data, which is the first normalized observation waveform after being compared with each other. The response of the corresponding directional filter bank in the three-dimensional domain of scan position × direction × time; The restored pulse tensor representing the direction index is given by... With directional nuclei { Under the constraints of}, the result is obtained after grouping sparse deconvolution and performing neighborhood smoothing in the directional dimension; the above "directional projection" refers to the same Position Pair along The dimensions are weighted and aggregated according to preset weights, and then normalized to their amplitudes, thus returning to... The domain forms a restored pulse train.

[0139] In another embodiment, in order to explicitly characterize the effect of the incident / reflection direction on echo aliasing under curved surface and gradually changing refractive index conditions, and to suppress artifacts caused by orientation mismatch in the directional dimension, this embodiment... Figure 1 Based on the overall process shown, a processing chain is introduced, which includes directional dimensionality increase, a set of sparse deconvolutions, directional neighborhood smoothing, and directional projection. Specifically, it includes steps S301 to S306.

[0140] Regarding S301:

[0141] The normal information of the surface being tested is collected and a one-to-one correspondence is established with the scanning coordinate system. For example, the surface point cloud of the molten bead region can be obtained through a profilometer or structured light 3D scanning, and the unit normal vector of each point is estimated according to a preset grid. The externally measured normal field is registered to the system scanning coordinate system via calibration points, and the corresponding normal vector table is obtained at each scanning line and each sampling position using linear or spline interpolation. To reduce the influence of external measurement noise, neighborhood median or spline smoothing is preferably applied to the normal field, and obvious abrupt changes are removed or backfilled according to threshold rules.

[0142] Regarding S302:

[0143] A first set of directions is determined based on the first geometric prior parameter and the normal vector table. Specifically, discrete sampling is performed in the effective angle domain related to incidence / reflection, and the direction sampling is adaptively encrypted in a section with large curvature variation and is appropriately sparse in a flat curvature area to form a direction list that covers sufficiently and has controllable redundancy.

[0144] Exemplarily, the number of directions can be 12-24, and the direction step is adaptively adjusted according to the local curvature radius. The direction set is stored in the form of an ordered list or a lookup table index, and is associated with the scanning position index for subsequent steps to share and call.

[0145] For S303:

[0146] A set of direction-dependent point spread functions is calculated based on the first reference waveform, the first set of directions, and the first geometric prior parameter, to obtain a first direction kernel cluster.

[0147] Specifically, taking the first point spread function as a reference, the reference response is directionally consistent transformed according to the possible propagation path difference and phase delay of each direction, and amplitude normalization and bandwidth consistency correction are performed to make the kernels of different directions comparable. Exemplarily, the can be pre-calculated and cached in the MATLAB or Python environment to support subsequent fast convolution / deconvolution operations.

[0148] For S304:

[0149] Direction filtering is performed based on the first normalized observation waveform and the first set of directions to construct first direction high-dimensional data. Specifically, the first normalized observation waveform is applied to each corresponding steerable filter or directional time-frequency atom (such as a direction-aligned Gabor atom) to extract a direction response in a frequency band matched with the bandwidth of the first direction kernel cluster, and output a direction high-dimensional data tensor in the form of To avoid boundary effects, mirror or zero padding boundary processing is preferably used; to maintain comparability between channels, energy normalization is performed on each direction channel.

[0150] For S305:

[0151] Group sparse deconvolution processing is performed based on the first direction high-dimensional data and the first direction kernel cluster to obtain a restored pulse tensor indexed by directions, and neighborhood smoothing is performed in the direction dimension. Specifically, the same position is used to cross all The amplitude vectors of the channels are collectively activated in a few neighboring directions as a "group" to comply with the physical characteristic that the direction varies slowly with position. Exemplarily, the group sparse solution can be implemented under the ADMM or iterative shrinkage threshold (ISTA) framework; to suppress the non-physical direction jump, further neighborhood smoothing or anisotropic regularization (e.g., weighted average of the neighboring direction channels) is applied to the output direction continuous deconvolution pulse tensor . In implementation, the direction channels below the energy threshold can be set to zero to improve robustness and computational efficiency. . In implementation, the direction channels below the energy threshold can be set to zero to improve robustness and computational efficiency.

[0152] For S306:

[0153] The direction continuous deconvolution pulse tensor is projected onto the time axis and the scanning axis according to the first direction set to obtain a deconvolution pulse train. Specifically, according to the direction response intensity, the normal consistency or the preset direction weight, the dimension is weighted and converged, and the energy of each direction kernel is normalized before projection to avoid amplitude bias caused by direction preference.

[0154] Exemplarily, the neighboring direction channels can be summed and thresholded at each position to obtain a sparse pulse representation in the time domain; the deconvolution pulse train obtained through the projection step returns to the domain, seamlessly connecting with the subsequent imaging and curvature correction steps in the above process.

[0155] Through the processing chain of S301-S306 described above, the present embodiment explicitly incorporates the direction element into the signal representation and deconvolution solving process, and collectively constrains the direction support of the echo in the direction dimension with group sparsity and neighborhood smoothing, thereby suppressing the time domain artifacts caused by direction mismatch and intra-pixel multi-angle mixing, and further improving the physical interpretability and spatial consistency of the deconvolution pulse train.

[0156] Optionally, the determination of the direction set based on the first geometric prior parameter and the normal information comprises:

[0157] generating an initial direction set based on the first geometric prior parameter and the normal information to obtain a first initial direction set;

[0158] calculating a direction indicator quantity within a preset time window and a preset scanning window based on the first normalized observation waveform to obtain a first direction indicator quantity;

[0159] calculating a direction consistency score based on the first direction indicator quantity and the first initial direction set to obtain a first direction consistency score;

[0160] Within the angular interval range defined based on the first geometric prior parameter, the first initial direction set is contracted and resampled based on the first direction consistency score to obtain a second direction set.

[0161] In another embodiment, the process of determining the direction set based on the first geometric prior parameter and the normal information comprises four steps of adaptive candidate generation, data validation, score aggregation, and limited update to form a direction set consistent with the local propagation direction under the conditions of curved surface and refractive index gradient.

[0162] First, according to the first geometric prior parameters such as the curvature radius of the bead, the thickness of the tube wall, the incident path angle, and in combination with the registered normal vector field, the feasible incident / reflection angle domain is determined for each scanning position, and the first initial direction set is obtained by discretizing the angle domain. In order to adapt to different curvature regions, the discretization density is preferably set adaptively according to the local curvature: smaller angle step is taken in the region with larger curvature, and the angle step is appropriately relaxed in the region with flat curvature. The direction set is stored in the form of an ordered list and bound with the scanning position index for subsequent calling.

[0163] Subsequently, the time window related to the main echo and the scanning window centered on the current scanning position are selected in the first normalized observation waveform, and the first direction indicator is calculated to reflect the "data evidence" of each candidate direction.

[0164] Exemplarily, a steerable filter or a directional time-frequency atom aligned with each candidate direction can be applied to obtain the directional response energy, and the energy is averaged or peaked in the window to obtain the energy component; at the same time, the consistency of the predicted phase sign corresponding to the candidate direction with the measured phase sign is obtained based on the analytical signal phase to obtain the phase consistency component; if necessary, the cross-correlation peak guided by the direction can be introduced as a morphology component to measure the accessibility of echo alignment under the direction assumption.

[0165] The above components are respectively amplitude-normalized and scale-unified to form a set of first direction indicator vectors indexed by direction.

[0166] Considering that the amplitude of terahertz is easily affected by attenuation in the high absorption band, the phase consistency component is given a higher weight, and the energy component is subjected to out-of-band suppression or sub-band weighting to reduce the interference of absorption on the judgment.

[0167] After obtaining the direction indicator, a direction consistency score model combining geometric prior and data evidence is established to score the first initial direction set one by one.

[0168] Exemplarily, the energy component, the phase consistency component and the morphology component can be weighted and aggregated into a first direction consistency score normalized to 0-1 with preset weights; the weights can be set through an offline verification set or a production sample experience, and can be adaptively fine-tuned according to the signal-to-noise ratio of the scanning position.

[0169] To improve the spatial stability of the score, the scores of the same direction within the scanning window can be smoothed or robustly aggregated in the neighborhood, so as to reduce the influence of local isolated noise.

[0170] After the score is obtained, the first initial direction set is contracted and resampled according to the angle interval constraint defined by the first geometric prior parameter to obtain a second direction set.

[0171] Specifically, when the interval of adjacent candidate directions is less than the minimum interval, they are merged by weighting according to their consistency scores, so as to avoid redundancy and overfitting caused by too dense discretization; when the interval of adjacent candidate directions is greater than the maximum interval and the scores at both ends are high, a new direction is inserted at the midpoint or the place with the maximum score gradient, so as to avoid missing caused by too sparse discretization.

[0172] To prevent the direction set from producing non-physical jumps between adjacent scanning positions, the embodiment limits the addition / deletion ratio of each update to be no more than a preset upper limit of the initial set size, for example, no more than 30%, and sets an upper limit constraint on the set difference of adjacent positions, for example, the total angular offset of the set difference is no more than several discrete levels.

[0173] Exemplarily, the minimum angle interval can be 2°-5°, and the maximum angle interval can be 8°-15°; the energy threshold can be set at 10%-30% of the global maximum value, and the phase consistency rate threshold can be set at 70%-90%; the above numerical values are exemplary ranges, and the actual setting can be made according to the bandwidth of the detection system and the absorption characteristics of the material.

[0174] In terms of implementation, the direction indicator calculation and score aggregation can be completed in MATLAB or Python environment by using steerable filtering, Hilbert analytic phase and cross-correlation library functions, and the direction set and the bidirectional index of the scanning position are maintained in a lookup table structure. The second direction set obtained through the above limited update serves as the adaptive direction basis for subsequent direction kernel cluster construction, direction dimension lifting and group sparse deconvolution, so as to meet the coverage and smoothness requirements in the direction dimension, and suppress the direction bias caused by normal measurement error or absorption.

[0175] Optionally, it further comprises:

[0176] calculating a set change quantity based on the first initial direction set and the second direction set to obtain a first direction stability condition;

[0177] output the second direction set as the first direction set when the first direction stability condition is satisfied; and repeat the generating of the direction set as a new initial direction set and the generating of the direction set until the first direction stability condition is satisfied when the first direction stability condition is not satisfied.

[0178] In another embodiment, to suppress the back-and-forth jitter and over-update of the adaptive direction set in low signal-to-noise ratio or curvature mutation area, and determine the convergence state of the direction set under the premise of satisfying the physical consistency, the first direction stability condition is formed based on the first initial direction set of the last round and the second direction set of the current round after the adaptive generation of the direction set is completed. When the stability condition is satisfied, the current direction set is frozen, otherwise the second direction set is taken as a new initial set to repeat the generation until the stability condition is satisfied.

[0179] Specifically, first, the direction sets of two rounds are paired under the same scan position index. To avoid mispairing, a bidirectional nearest neighbor matching in ascending order of angle is adopted, and the allowed angle difference is set to 1-2°, and elements exceeding this range are no longer considered as the same direction. Thus, the elements in the last round set that are not matched are recorded as deletion candidates, and the newly appeared elements in the current round set are recorded as new addition candidates; the new directions obtained by merging in the last round are established to have a source mapping, so as to subsequently count the angle drift and convergence trend.

[0180] After pairing is completed, the set change quantity of the scan position is calculated. The set change quantity is a multi-component weighted aggregation index, at least including: the proportion of the number of new additions and deletions to the size of the last round set, for measuring the size change; the average angle drift and the maximum angle drift between the matched directions, for reflecting the overall translation of the directions; the overall separation degree of the two discrete sets on the angle axis, for which the maximum nearest neighbor difference process is adopted, that is, the angle difference between each direction in set one and the nearest neighbor in set two is calculated and the maximum value is taken, and the larger one is taken with the reverse calculation result, and is normalized in the full angle domain; the change of the mean and dispersion (such as variance or quantile difference) of the direction consistency score, for reflecting the stability of the data evidence; and the spatial smoothness across positions, for which the 90% quantile of the paired angle offset is counted in the neighborhood (for example, ±2-5 positions) centered on the current scan position and normalized, to suppress local jumps.

[0181] Considering that the amplitude reliability of terahertz is low in the high absorption band, the weight of the components related to the stability of the phase label and the monotonicity of the group delay is increased in the aggregation, for example, to about 0.3-0.5, and the energy-related components are weighted in subbands or suppressed out of band, to enhance the recognition ability of the physical consistency.

[0182] To avoid interference determination in low signal-to-noise ratio positions, a mask that does not count in stability determination is constructed. The specific method is: estimate the noise envelope in the baseline interval outside the main peak, calculate the ratio of the main peak energy to the noise energy, and when the ratio is lower than 6-10 dB, mark the scanning position as a mask point; The change of the mask point set is only used for log recording and is not included in the global aggregation and stability determination. Then, the global and local joint determination is carried out to form the first direction stability condition: when the global set change is lower than the global change threshold, and the set change of each scanning position is lower than the local change threshold, the proportion of the position is not lower than the spatial proportion lower limit, for example, 80%-90%, and the above conditions are all met in the last 2-4 rounds of update, that is, it is considered to reach stability.

[0183] To avoid jitter caused by repeated crossing of the threshold value, a consistent hysteresis is set for entering / leaving: after entering stability, if there is slight fluctuation, a stricter exit threshold group is used for determination, that is, the global and local thresholds are each increased by about 10%, and the spatial proportion lower limit is reduced by about 5 percentage points, and it is exited when it does not meet for 1 round, thereby improving the robustness of decision-making.

[0184] When the first direction stability condition is met, the second direction set of the current round is frozen as the new first direction set, which is used for subsequent direction kernel cluster construction, direction dimensionality increase, and group sparse deconvolution; When the stability condition is not met, the second direction set is taken as the new initial direction set, and the direction set generation process is returned to continue iteration. To suppress update shock and control computational overhead, each round of update applies a step limit: the net increase / deletion ratio is not higher than 20%-30% of the last round of set size, the maximum angle drift allowed at a time is not more than 2-4°, the total number of directions is kept within a preset interval, for example, not less than 8-12, and not higher than 1.5 times the initial set; In the area where the curvature mutation is significant, the above upper limit can be tightened accordingly to prevent non-physical jumps.

[0185] In terms of implementation, the set pairing, angle difference bit number, and set distance can be calculated in a common numerical environment, and the analytical phase can be obtained through Hilbert transform; Each round records the update log, including the net increase / deletion number, the average and maximum angle drift, the set distance, the score change, the spatial consistency, the mask proportion, and whether it is frozen and the threshold group adopted, which is used for subsequent audit and parameter setting.

[0186] Through the above control, the adaptive direction set obtains good convergence and reproducibility in the terahertz detection scene of curved surface and refractive index gradient, thereby improving the reliability of the subsequent deconvolution and imaging results.

[0187] Optionally, the group sparse deconvolution processing based on the first direction dimensionality increased data and the first direction kernel cluster to obtain the recovered pulse tensor of the direction index comprises:

[0188] determining a direction adjacency relationship based on the first geometric prior parameter and the normal information, and generating a first direction compatibility matrix;

[0189] imposing a sharing constraint on the non-zero support of the direction dimension according to the first direction compatibility matrix in the group sparse solution, limiting the non-zero directions corresponding to the same time sample to the adjacent directions specified by the first direction compatibility matrix, and outputting the recovered impulse tensor of the direction index.

[0190] In an embodiment, in order to suppress non-physical discrete jumps in the direction dimension and make the direction activation at the same time sample consistent with the detection surface geometry continuity, the embodiment first determines a direction adjacency relationship based on the first geometric prior parameter and the normal information, and generates a first direction compatibility matrix for constraint solving.

[0191] Specifically, for each spatial sampling point on each scan line, the registered normal vector and the corresponding direction candidate set are used to determine whether any two directions are allowed to be considered "adjacent" according to the angle adjacency and geometric continuity principle.

[0192] In implementation, an adjacency radius matching the direction discrete step can be set, so that the directions with an angle difference within the radius are marked as "adjacent", and the directions that remain smooth in adjacent scan positions are also marked as "adjacent". Directions that do not meet the above conditions are marked as "non-adjacent". The above adjacency relationship is recorded in a sparse binary manner, i.e., a first direction compatibility matrix is formed, where each row and each column corresponds to a discrete direction element. The allowed value in the matrix indicates that the direction and another direction can appear together in solving, and the disallowed value indicates that the two directions should not be activated at the same time sample.

[0193] In order to reduce the influence of noise on the adjacency relationship, it is preferred to smooth the normal field before constructing the matrix, and to use mirroring or continuation at the boundary of the direction set to avoid the boundary direction being incorrectly excluded due to insufficient neighbors.

[0194] After obtaining the first direction compatibility matrix, a sharing constraint is introduced in the group sparse deconvolution solving process: the non-zero directions allowed to appear at the same time sample should form a connected subset connected by "adjacent" relationships.

[0195] The embodiment adopts an iterative strategy of "estimation, projection, and refinement": first, a preliminary group sparse deconvolution is performed on the direction elevation data and the direction kernel cluster without applying the direction compatibility constraint to obtain the amplitude estimation of the direction index; then, a "feasible projection" is performed on the amplitude estimation according to the first direction compatibility matrix, that is, the direction activation at the same time sampling is screened, only the adjacent directions connected with the current main active direction are retained, and the sporadic activation falling on the non-adjacent relationship is suppressed or incorporated into the nearest adjacent connected block; then, the convolution residual is recalculated on the new feasible support and the amplitude is updated. The above steps of "estimation, projection, and refinement" are alternately performed until the change of the residual is lower than a preset threshold or the iteration round reaches an upper limit.

[0196] In a specific implementation, the "feasible projection" can be specifically implemented as a combination of two types of operations: one type is support clipping based on the compatibility matrix, that is, the isolated direction not meeting the adjacent condition is directly set to zero; the other type is amplitude redistribution based on the adjacent direction, that is, the amplitude is smoothed within a small range of adjacent set according to the energy and phase consistency, so as to weaken the direction sawtooth caused by noise. In order to ensure numerical stability, it is preferred to perform a light neighborhood smoothing on the direction dimension at the end of each iteration, while keeping the sparsity of the time dimension from being destroyed.

[0197] In terms of parameter selection, in order to make the constraint neither too strict to swallow the real multi-path direction nor too loose to allow non-physical activation, the embodiment sets the direction adjacency radius to a range of the same order of magnitude as the local angular discrete interval, preferably 1-2 times the local angular discrete interval; the angular discrete interval is the included angle between two adjacent discrete directions; in the adaptive direction set, the angular discrete interval is dynamically determined according to the scanning position. The allowed adjacent depth at the same time sampling can be one to two discrete levels; the minimum amplitude threshold triggered by the "feasible projection" and the iteration stop threshold can be adjusted according to the system noise level, for example, the iteration is terminated when the convolution residual decreases by a few percentage in a plurality of consecutive iterations.

[0198] For example, the direction compatibility matrix can be pre-generated and cached as a sparse adjacency structure; the group sparse solution and the feasible projection can be implemented in a common numerical calculation environment; in order to avoid that the boundary direction is systematically suppressed, the energy of each direction channel can be normalized before projection to balance the response difference of different direction kernels.

[0199] In this way, the final output direction index restores the pulse tensor, which presents a connected and smooth support distribution in the direction dimension and keeps a sparse pulse form in the time dimension, and the result is used for subsequent direction projection and imaging process, which can effectively suppress the time domain artifacts caused by direction mismatch in the scene of curved surface and refractive index gradient.

[0200] Optionally, the determining of the direction adjacency relationship based on the first geometric prior parameter and the normal information, and generating a first direction compatibility matrix comprises:

[0201] The equivalent incident angle corresponding to each direction is determined based on the first geometric prior parameter and the normal information, and the reflection coefficient phase characteristics corresponding to the equivalent incident angle are evaluated in the terahertz working frequency band, and phase-incompatible direction pairs are marked based on in-band phase sign consistency to obtain a first incompatible set;

[0202] After excluding the first incompatible set, the main echo group time delay order of adjacent scanning positions is calculated based on the first geometric prior parameter, and the direction pairs are screened based on the group time delay order as a criterion to obtain a first compatible set;

[0203] The direction pairs belonging to the first compatible set are set as allowed adjacency in the first direction compatibility matrix, and the direction pairs belonging to the first incompatible set are set as forbidden adjacency in the first direction compatibility matrix.

[0204] In an embodiment, in order to match the adjacent relationship of the direction dimension with the terahertz electromagnetic propagation characteristics, the surface geometric characteristics, and the refractive index gradient scene, the embodiment determines the direction adjacency relationship based on the first geometric prior parameter and the normal information, and generates a direction compatibility matrix based thereon to constrain the group sparse deconvolution solving process.

[0205] Specifically, first, for the scanning region of the polyethylene hot melt joint, the surface normal and geometric parameters such as the bead curvature radius and the pipe wall thickness at the scanning position are used to determine an initial direction set. For example, at a certain typical scanning position, the initial direction set calculated from the normal can be discretized into angles such as 0°, 5°, 10°, 15°, 20°, and the like. These discrete angles represent the included angles of the incident and reflected directions relative to the local normal, and serve as the initial set for subsequent analysis of direction compatibility.

[0206] Subsequently, the equivalent incident angle is calculated for each initial direction, and the amplitude and phase of the reflection coefficient are determined by using the Fresnel formula according to the reflection characteristics of the polyethylene material in the terahertz working frequency band, such as the 0.1 THz to 3 THz frequency band.

[0207] For example, for angles 0°, 5°, and 10°, the reflection phase signs calculated near the frequency point of 0.5 THz can be positive, positive, and negative, respectively, which indicates that angle 10° has stable phase opposition relative to 0° and 5°, that is, 10° belongs to the “phase-incompatible” relationship with the former two. This phase-incompatible relationship is recorded as the first incompatible set, which is used for subsequent construction of the direction compatibility matrix constraint.

[0208] After excluding the first incompatible set, the embodiment further evaluates the main echo group delay order of adjacent scanning positions. Specifically, for example, at scanning positions x = 100 and x = 101, if the main echo group delay orders corresponding to angles 5° and 10° in the direction set are both “5° earlier than 10°”, it is considered that the direction pair meets the group delay monotonicity and belongs to the direction compatible case; otherwise, if the order is reversed at x = 101, such as 10° earlier than 5°, it is considered that the geometry continuity is not met and is not compatible. Through such judgment, the direction pairs that meet the group delay monotonicity are recorded as the first compatible set, which is used for the construction of the direction compatibility matrix.

[0209] On the basis of obtaining the first incompatible set and the first compatible set, the first direction compatibility matrix is constructed. The matrix takes the initial direction set as the index, and each matrix element represents whether two directions can be activated at the same time sampling. When the matrix is initialized, all elements are set to prohibit adjacency; then the direction pairs belonging to the first compatible set are marked as allowing adjacency, for example, angles 0° and 5° are allowed to be adjacent, and the corresponding position of the matrix is marked as allowed. For the direction pairs belonging to the first incompatible set, for example, there is a phase exclusion relationship between 5° and 10°, which is marked as prohibited adjacency.

[0210] To ensure the robustness of numerical solution, the embodiment further performs sparsification and symmetrization processing after the matrix generation is completed, to eliminate isolated allowed items and ensure the mutual consistency of the adjacent relationship of the direction pairs. At the boundary of the direction set, the minimum adjacency can be supplemented by mirroring or extension to avoid the boundary direction being systematically suppressed in the projection process due to the lack of adjacent channels.

[0211] In the group sparse deconvolution solving stage, the non-zero direction activation at the same time sampling is limited within the allowed adjacent relationship specified by the above-mentioned first direction compatibility matrix. In the specific solving process, first, the initial direction activation estimation is obtained by performing preliminary group sparse solving without applying the adjacency constraint; then, the initial activation at each time sampling is projected according to the direction compatibility matrix, that is, the isolated activation that does not meet the adjacency relationship is suppressed or merged into the nearest neighbor compatible direction. The above-mentioned process is alternately performed through multiple iterations until the convolution residual changes below a set threshold (for example, 1–3%) or reaches the maximum number of iterations, for example, 50–100 times, and finally the recovered pulse tensor of the direction index is obtained.

[0212] Through the above-mentioned direction compatibility matrix construction and solving process, in the terahertz detection scene of curved surface and refractive index gradient, the non-zero activation in the direction dimension is always limited within the direction set that meets the geometric continuity and phase consistency in the deconvolution process, thereby effectively suppressing the non-physical direction jump and time domain artifact, and improving the physical consistency and spatial positioning accuracy of the final recovered pulse train.

[0213] Optionally, the obtaining the defect distribution map comprises:

[0214] Based on the first geometric prior parameter, an equivalent incidence angle of each scanning position is determined, a representative frequency of each recovered pulse is determined based on a spectral energy center of the recovered pulse train, and a first representative frequency set is obtained;

[0215] Based on the equivalent incidence angle and the first representative frequency set, an interface reflection phase label is calculated, and a first phase label set is obtained; based on an analytical signal phase of the recovered pulse train at a pulse peak, a measured phase label is calculated, and a second phase label set is obtained;

[0216] Based on the first phase label set and the second phase label set, phase consistency verification is performed, and when the phase labels are inconsistent, polarity correction and micro-time delay correction are performed on the corresponding recovered pulse, and a phase-aligned recovered pulse train is obtained;

[0217] Based on the phase-aligned recovered pulse train, a time-of-flight map is generated, and an angle-resolved group refractive index is determined based on the equivalent incidence angle and the first representative frequency set, curvature correction is performed, and the defect distribution map is obtained.

[0218] In an embodiment, to avoid false depth and geometric distortion caused by terahertz phase flipping and dispersion, the embodiment introduces phase label calibration and angle-resolved group refractive index mapping after obtaining the recovered pulse train, and completes time-of-flight map construction and curvature correction based on this, and outputs the defect distribution map.

[0219] In a specific implementation, first, based on the first geometric prior parameter, an equivalent incidence angle is determined for each scanning line and each scanning position. The equivalent incidence angle is calculated according to the local normal of the detection surface, the transmission and reception geometry, and is associated with the scanning coordinates point by point. On this basis, the representative frequency of the recovered pulse train is extracted pulse by pulse.

[0220] The specific method is as follows: a time window is set near the pulse peak, for example, 5-15 sampling points around the peak, windowing and fast spectral analysis are performed on the signal in the window, and the spectral energy center is calculated as the representative frequency of the pulse; when the system bandwidth is wide, a plurality of sub-bands can be divided in the effective frequency band, and the energy centers are calculated respectively, and finally the representative frequency set is obtained by majority voting or weighted average.

[0221] Secondly, the expected interface reflection phase label is calculated in combination with the equivalent incidence angle and the representative frequency set to form a first phase label set. To this end, the Fresnel reflection model can be used and the polyethylene refractive index and absorption parameters in the material library are called to evaluate the phase sign and transition characteristics of the reflection coefficient under different incidence angles and polarization conditions within the working frequency band; when the sub-band strategy is used, the phase signs within each sub-band are counted, and the majority consistency is used as the criterion, and if necessary, the frequency band near the phase transition interval is reduced in weight to reduce the edge uncertainty.

[0222] At the same time, the measured phase information is obtained from the recovered pulse train: the analytical signal phase is calculated at the peak of each pulse, and the phase is unwrapped in combination with the neighborhood expansion to obtain the measured phase label corresponding to the time window to form a second phase label set.

[0223] In addition, to enhance robustness, the measured phase label can be smoothed for neighborhood consistency between adjacent scanning positions, and the pulse is set not to participate in marking when the signal-to-noise ratio of the main peak is too low.

[0224] Thirdly, consistency verification and correction are performed on the two sets of phase labels. When the first phase label and the second phase label are inconsistent, polarity correction is first performed on the corresponding pulse according to the polarity rule, i.e., the pulse is flipped as a whole; if there is still a small phase residual, micro-time delay correction is performed under the premise of maintaining the polarity. Micro-time delay correction can adopt two types of implementation paths: one is to perform sub-sampling accuracy interpolation and relocation on the pulse peak neighborhood in the time domain, such as parabolic interpolation or fractional delay interpolation based on window function, to make a small translation of the pulse along the time axis; the other is to apply a linear phase term to the pulse in the frequency domain to achieve equivalent micro-delay compensation. After correction, a phase-aligned recovered pulse train is obtained.

[0225] To avoid over-correction, the embodiment only triggers the above operation when the phase is inconsistent and the signal-to-noise ratio meets the standard, and sets an upper limit to the time translation allowed at a time, for example, not more than half a sampling interval.

[0226] Subsequently, a time-of-flight map is generated based on the phase-aligned recovered pulse train: the arrival time of the main echo and the key secondary echo is extracted at each scanning position, stacked to form a two-dimensional time-of-flight profile set in the scanning order, and the time-of-flight size is represented by gray scale or pseudo-color mapping. To stably convert the time-of-flight into geometric depth, instead of using a fixed refractive index, the embodiment uses an angle-resolved group refractive index to reconstruct the path length.

[0227] The specific method is as follows: according to the equivalent incidence angle and the representative frequency set, the group refractive index under the angle-frequency combination is obtained by querying or interpolating in the material dispersion library, and the time-of-flight is converted into the propagation distance in the medium point by point; when the sub-band strategy is used, the group refractive index of each sub-band can be weighted according to the energy weight to obtain the effective group refractive index of the pulse.

[0228] Exemplarily, in the frequency band commonly used by polyethylene, the value of the group refractive index can slowly change in the range of about 1.47-1.55 with the frequency, and a slight difference can also appear under different incident angles, which is compensated by the above-mentioned angle resolution mapping.

[0229] After the time-distance conversion is completed, curvature correction is performed to eliminate the geometric distortion caused by the bead surface. Specifically, a geometric model of the detection section is established according to the first geometric prior parameter, and the depth value based on the path length is transformed from the scanning coordinate system to the target geometric coordinate system.

[0230] The polar coordinate expansion method can be used to expand the isobath layer along the curved surface to the equidistant layer in the rectangular coordinate; for the area with rapid curvature change, the local expansion window is preferably reduced and the neighborhood smoothing is introduced to reduce numerical instability. The image obtained after curvature correction is defined as a defect distribution map: the unfused defects usually show local depth mutation or energy void; the pores or inclusions can show point-like or patch-like abnormal areas; the interlayer separation usually presents a continuous shallow abnormal band after phase alignment.

[0231] In actual use, the isobath or energy contour line can be superimposed on the defect distribution map to assist in interpretation, and the size, connectivity and morphological threshold are set to output the region label and spatial coordinates of the candidate defects.

[0232] In this way, the time-of-flight to spatial depth mapping avoids the false depth caused by phase reversal, while maintaining the consistency of the time domain and the geometric domain under the conditions of curved surface and refractive index gradient, so that the final defect distribution map is more consistent with the actual structure in position and shape, providing a reliable basis for subsequent quantitative evaluation and engineering disposal.

[0233] Based on the same inventive concept, the embodiments of the present application also provide a terahertz signal aliasing elimination system based on deconvolution corresponding to the terahertz signal aliasing elimination method based on deconvolution. Since the system in the embodiments of the present application solves the problem by similar principles as the above-mentioned terahertz signal aliasing elimination method based on deconvolution, the implementation of the system can be referred to the implementation of the method, and the repeated parts will not be described here.

[0234] Referring to Figure 4 Fig. 1 shows a schematic diagram of a terahertz signal aliasing elimination system based on deconvolution provided by the embodiments of the present application, the system comprises:

[0235] The acquisition module 10 is configured to acquire the time-domain waveform of the reference region and the detection region to obtain a first reference waveform and a first observation waveform; and obtain the geometric and material parameters to obtain a first geometric prior parameter;

[0236] The processing module 20 is configured to calculate a point spread function based on the first reference waveform and the first geometric prior parameter, to obtain a first point spread function; calculate a time-axis normalization mapping function based on the first geometric prior parameter, to obtain a first time normalization mapping function; and apply the first time normalization mapping function to the first observation waveform, to obtain a first normalized observation waveform.

[0237] The recovering module 30 is configured to perform deconvolution processing based on the first normalized observation waveform, to obtain a recovered pulse train.

[0238] The output module 40 is configured to generate a time-of-flight map based on the recovered pulse train, and perform curvature correction based on the first geometric prior parameter, to obtain a defect distribution map.

[0239] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application, rather than limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A deconvolution-based method for terahertz signal aliasing cancellation, characterized in that, The method comprises: collecting time-domain waveforms of a reference region and a detection region to obtain a first reference waveform and a first observation waveform; obtaining geometric and material parameters to obtain a first geometric prior parameter; calculating a point spread function based on the first reference waveform and the first geometric prior parameter to obtain a first point spread function; calculating a time-axis normalization mapping function based on the first geometric prior parameter to obtain a first time normalization mapping function; applying the first time normalization mapping function to the first observation waveform to obtain a first normalized observation waveform; performing deconvolution processing based on the first normalized observation waveform to obtain a recovered pulse train; generating a time-of-flight map based on the recovered pulse train and performing curvature correction based on the first geometric prior parameter to obtain a defect distribution map; The method further comprises: performing deconvolution processing on the first observation waveform on an unnormalized original time axis to obtain a second recovered pulse train; applying the first time normalization mapping function to the second recovered pulse train to obtain a second normalized recovered pulse train; calculating an exchangeability difference value based on the recovered pulse train and the second normalized recovered pulse train to obtain a first coupling constraint value; performing joint iterative updating based on the first coupling constraint value, a convolution residual between the first normalized observation waveform and the first point spread function, and a sparsity index of the recovered pulse train to obtain an updated first time normalization mapping function, an updated first point spread function, and an updated recovered pulse train; The method further comprises: calculating a first stability index and a second stability index based on the updated recovered pulse train to obtain a first stopping condition; terminating the joint iterative updating when the first stopping condition is met, and outputting a target recovered pulse train as the recovered pulse train.

2. The deconvolution-based method for terahertz signal aliasing elimination according to claim 1, wherein, The deconvolution processing based on the first normalized observation waveform comprises: collecting normal information of a detection surface, determining a direction set based on the first geometric prior parameter and the normal information to obtain a first direction set; calculating a set of direction-dependent point spread functions based on the first reference waveform and the first direction set to obtain a first direction kernel cluster; performing direction filtering based on the first normalized observation waveform and the first direction set to obtain first direction-enhanced data; performing group sparse deconvolution processing based on the first direction-enhanced data and the first direction kernel cluster to obtain a direction-indexed recovered pulse tensor; performing direction neighborhood smoothing processing on the direction-indexed recovered pulse tensor to obtain a direction-continuous recovered pulse tensor; projecting the direction-continuous recovered pulse tensor to the time axis and the scanning axis according to the first direction set to obtain the recovered pulse train.

3. The deconvolution-based method for terahertz signal aliasing elimination according to claim 2, wherein, The determination of the direction set based on the first geometric prior parameter and the normal information comprises: generating an initial direction set based on the first geometric prior parameter and the normal information to obtain a first initial direction set; calculating a direction indicator based on the first normalized observation waveform within a preset time window and a preset scanning window to obtain a first direction indicator; calculating a direction consistency score based on the first direction indicator and the first initial direction set to obtain a first direction consistency score; Within an angle interval range defined based on the first geometric prior parameter, a shrinking and resampling is performed on the first initial direction set based on the first direction consistency score, to obtain a second direction set.

4. The deconvolution-based terahertz signal aliasing elimination method of claim 3, wherein, Further comprising: A set change quantity is calculated based on the first initial direction set and the second direction set, to obtain a first direction stability condition; When the first direction stability condition is met, the second direction set is output as the first direction set; when the first direction stability condition is not met, the second direction set is taken as a new initial direction set, and the generation step of the direction set is repeated until the first direction stability condition is met.

5. The deconvolution-based terahertz signal aliasing elimination method of claim 2, wherein, The group sparse deconvolution processing based on the first direction high-dimensional data and the first direction kernel cluster to obtain the recovered pulse tensor of the direction index comprises: A direction adjacency relationship is determined based on the first geometric prior parameter and the normal information, to generate a first direction compatibility matrix; In the group sparse solution, a sharing constraint is imposed on the non-zero support of the direction dimension according to the first direction compatibility matrix, so that the non-zero directions corresponding to the same time sampling are limited to the adjacent directions specified by the first direction compatibility matrix, and the recovered pulse tensor of the direction index is output.

6. The deconvolution-based terahertz signal aliasing elimination method of claim 5, wherein, The first direction compatibility matrix is generated by determining the direction adjacency relationship based on the first geometric prior parameter and the normal information, which comprises: Based on the first geometric prior parameter and the normal information, the equivalent incidence angle corresponding to each direction is determined, and the reflection coefficient phase characteristics corresponding to the equivalent incidence angle are evaluated within the terahertz working frequency band, based on the in-band phase sign consistency, the phase-contradictory direction pairs are marked to obtain a first incompatible set; After excluding the first incompatible set, the main echo group time delay order of the adjacent scanning positions is calculated based on the first geometric prior parameter, and the direction pairs are screened based on the group time delay order as the criterion to obtain a first compatible set; The direction pairs belonging to the first compatible set are set as allowed adjacency in the first direction compatibility matrix, and the direction pairs belonging to the first incompatible set are set as prohibited adjacency in the first direction compatibility matrix.

7. The deconvolution-based terahertz signal aliasing elimination method of claim 1, wherein, The defect distribution map comprises: Based on the first geometric prior parameter, the equivalent incidence angle of each scanning position is determined, and the representative frequency of each recovered pulse is determined based on the spectral energy center of the recovered pulse train, to obtain a first representative frequency set; Based on the equivalent incidence angle and the first representative frequency set, an interface reflection phase label is calculated to obtain a first phase label set; based on the analytical signal phase of the recovered pulse train at the pulse peak, a measured phase label is calculated to obtain a second phase label set; Based on the first phase label set and the second phase label set, a phase consistency check is performed, and when they are inconsistent, polarity correction and micro-time delay correction are performed on the corresponding recovered pulse to obtain a phase-aligned recovered pulse train; Based on the phase-aligned recovered pulse train, a time-of-flight map is generated, and based on the equivalent incidence angle and the first representative frequency set, an angle-resolved group refractive index is determined, and curvature correction is performed to obtain the defect distribution map.

8. A deconvolution-based terahertz signal aliasing cancellation system, characterized by, Comprise: The acquisition module is configured to acquire time-domain waveforms of the reference region and the detection region to obtain a first reference waveform and a first observation waveform; The geometric and material parameters are obtained to obtain a first geometric prior parameter; The processing module is configured to calculate a point spread function based on the first reference waveform and the first geometric prior parameter to obtain a first point spread function; The first geometric prior parameter is used to calculate a time-axis normalization mapping function to obtain a first time normalization mapping function; and the first observation waveform is applied with the first time normalization mapping function to obtain a first normalized observation waveform; The restoration module is configured to perform deconvolution processing on the first observation waveform on an unnormalized original time axis to obtain a second restored pulse train; The first time normalization mapping function is applied to the second restored pulse train to obtain a second normalized restored pulse train; An exchangeability difference value is calculated based on the restored pulse train and the second normalized restored pulse train to obtain a first coupling constraint value; Based on the first coupling constraint value, a convolution residual between the first normalized observation waveform and the first point spread function, and a sparsity index of the restored pulse train, joint iterative updating is performed to obtain an updated first time normalization mapping function, an updated first point spread function, and an updated restored pulse train; First and second stability indexes are calculated based on the updated restored pulse train to obtain a first stop condition; When the first stop condition is met, the joint iterative updating is terminated, and a target restored pulse train is output as the restored pulse train; The output module is configured to generate a time-of-flight map based on the restored pulse train, and perform curvature correction based on the first geometric prior parameter to obtain a defect distribution map.

Images