Terahertz signal aliasing elimination method and system based on deconvolution

Abstract

The invention provides a terahertz signal aliasing elimination method and system based on deconvolution, and relates to the technical field of detection, and the method comprises the steps: collecting time domain waveforms of a reference region and a detection region, and obtaining geometric and material parameters; calculating a point spread function based on the reference waveform and geometric prior; calculating a time axis standardization mapping function based on geometric prior, and carrying out standardization processing on the observation waveform; performing deconvolution processing on the normalized time axis to obtain a restored pulse train; and generating a flight time chart based on the recovered pulse train, and executing curvature correction in combination with geometric priori to obtain a defect distribution diagram. According to the method, by introducing reference region prior and time standardization, point spread function drift caused by curvature and refractive index gradual change is effectively overcome, noise amplification in a high-absorption medium is inhibited, signal aliasing can be stably removed in a polyethylene hot melting joint scene, and accurate defect positioning can be achieved.

Description

Technical Field

[0001] This application relates to the field of detection technology, and more specifically, to a method and system for eliminating aliasing of terahertz signals based on deconvolution. Background Technology

[0002] Terahertz time-domain spectroscopy and imaging are widely used for non-destructive testing of pipe fusion joints due to their strong penetrability to non-metallic materials such as polyethylene. During the testing process, terahertz pulses undergo multiple reflections and scatterings as they propagate through multilayer polyethylene fusion joints, forming overlapping echo signals. To mitigate the overlap effect, existing research often employs common signal processing methods such as window function truncation, frequency domain filtering, and wave packet decomposition. These methods can separate some overlapping information when detecting flat or parallel layered structures, but their basic assumption is that the system response remains spatially stable. However, fusion joints have curvature and molten bead geometry characteristics, accompanied by a gradual change in radial refractive index, causing the point spread function to drift significantly with position. Traditional frequency domain or time domain processing methods struggle to establish a physically consistent model.

[0003] Existing technologies still have significant limitations when dealing with curved weld joints and highly absorbent polyethylene media: on the one hand, there is a lack of signal characterization models that can be adaptively updated with geometric changes; on the other hand, under strong absorption conditions, conventional methods are prone to introducing noise or geometric distortion while eliminating aliasing, making it difficult to obtain a defect distribution that corresponds to the real structure.

[0004] Therefore, how to overcome the adverse effects of point spread function drift and time axis nonlinearity in curved hot-melt joints and strongly absorbing polyethylene media, so as to robustly eliminate echo aliasing and obtain a defect distribution that matches the geometry, has become an urgent technical problem to be solved. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this application provides a method and system for terahertz signal aliasing elimination based on deconvolution.

[0006] In a first aspect, this application provides a terahertz signal aliasing cancellation method based on deconvolution, comprising:

[0007] The time-domain waveforms of the reference region and the detection region are acquired to obtain the first reference waveform and the first observed waveform; the geometric and material parameters are obtained to obtain the first geometric prior parameters;

[0008] The first point spread function is obtained by calculating the point spread function based on the first reference waveform and the first geometric prior parameters;

[0009] The time axis normalized mapping function is calculated based on the first geometric prior parameters to obtain the first time normalized mapping function;

[0010] Applying the first time normalization mapping function to the first observed waveform yields the first normalized observed waveform;

[0011] Deconvolution processing is performed on the first normalized observed waveform to obtain the restored pulse train;

[0012] A time-of-flight map is generated based on the restored pulse train, and curvature correction is performed based on the first geometric prior parameters to obtain a defect distribution map.

[0013] Optional, also includes:

[0014] The first observed waveform is deconvolved on the unnormalized original time axis to obtain the second restored pulse train;

[0015] Applying the first time normalization mapping function to the second restored pulse train yields the second normalized restored pulse train;

[0016] Based on the commutability difference value between the restored pulse train and the second normalized restored pulse train, the first coupling constraint value is obtained;

[0017] Based on the first coupling constraint value, the convolution residual between the first normalized observation waveform and the first point spread function, and the sparsity index of the restored pulse train, a joint iterative update is performed to obtain the updated first time normalized mapping function, the updated first point spread function, and the updated restored pulse train.

[0018] Optional, also includes:

[0019] Based on the updated restored pulse train, the first stability index and the second stability index are calculated to obtain the first stopping condition;

[0020] The joint iterative update is terminated when the first stopping condition is met, and the target recovery pulse train is output as the recovery pulse train.

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

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

[0023] Based on the first reference waveform and the first direction set, calculate the direction-related point diffusion function set to obtain the first direction kernel cluster;

[0024] Based on the first normalized observation waveform and the first direction set, directional filtering is performed to obtain the first direction up-dimensional data;

[0025] Based on the first direction up-dimensional data and the first direction kernel cluster, perform sparse deconvolution processing to obtain the restored pulse tensor of the direction index;

[0026] Perform directional neighborhood smoothing on the restored pulse tensor of the directional index to obtain a directional continuous restored pulse tensor;

[0027] The continuous recovery pulse tensor of the direction is projected onto the time axis and the scan axis according to the first direction set to obtain the recovery pulse train.

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

[0029] An initial direction set is generated based on the first geometric prior parameters and the normal information to obtain the first initial direction set;

[0030] The first direction indication is obtained by calculating the direction indication within a preset time window and a preset scanning window based on the first normalized observation waveform.

[0031] The first direction consistency score is obtained by calculating the direction consistency score based on the first direction indication and the first initial direction set;

[0032] Within the angular spacing range defined by the first geometric prior parameter, the first initial direction set is shrunk and resampled based on the first direction consistency score to obtain the second direction set.

[0033] Optional, also includes:

[0034] Based on the first initial direction set and the second direction set, the change in the set is calculated to obtain the stability condition of the first direction;

[0035] 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 used as a new initial direction set and the direction set generation steps are repeated until the first direction stability condition is met.

[0036] Optionally, the step of performing group sparse deconvolution processing based on the first direction up-dimensional data and the first direction kernel cluster to obtain the restored pulse tensor of the direction index includes:

[0037] Based on the first geometric prior parameters and the normal information, the directional adjacency relationship is determined, and a first directional compatibility matrix is ​​generated;

[0038] In the group sparse solution, a shared constraint is applied to the non-zero support of the direction dimension based on the first direction compatibility matrix, so that the non-zero direction corresponding to the sampling at the same time is limited to the adjacent direction specified by the first direction compatibility matrix, and the restored pulse tensor of the direction index is output.

[0039] Optionally, the step of determining the directional adjacency relationship based on the first geometric prior parameters and the normal information, and generating the first directional compatibility matrix, includes:

[0040] Based on the first geometric prior parameters and the normal information, the equivalent incident angle corresponding to each direction is determined, and the phase characteristics of the reflection coefficient corresponding to the equivalent incident angle are evaluated within the terahertz working frequency band. Based on the in-band phase sign consistency, mutually exclusive direction pairs are marked to obtain the first incompatible set.

[0041] After excluding the first incompatible set, the group delay order of the main echoes at adjacent scanning positions is calculated based on the first geometric prior parameters, and the direction pairs are filtered based on the criterion of maintaining the group delay order to obtain the first compatible set.

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

[0043] Optionally, obtaining the defect distribution map includes:

[0044] The equivalent incident angle of each scanning position is determined based on the first geometric prior parameters, and the representative frequency of each restored pulse is determined based on the spectral energy centroid of the restored pulse train, thus obtaining the first representative frequency set.

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

[0046] Phase consistency verification is performed based on the first phase tag set and the second phase tag set. When there is inconsistency, polarity correction and micro-delay correction are performed on the corresponding recovery pulse to obtain a phase-aligned recovery pulse train.

[0047] A time-of-flight map is generated based on the phase-aligned restored pulse train, and an angle-resolved group refractive index is determined based on the equivalent incident angle and the first representative frequency set. Curvature correction is then performed to obtain the defect distribution map.

[0048] Secondly, this application provides a terahertz signal aliasing cancellation system based on deconvolution, comprising:

[0049] The acquisition module is used to acquire time-domain waveforms of the reference area and the detection area to obtain the first reference waveform and the first observed waveform; and to acquire geometric and material parameters to obtain the first geometric prior parameters.

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

[0051] The restoration module is used to perform deconvolution processing based on the first normalized observed waveform to obtain the restored pulse train;

[0052] The output module is used to generate a time-of-flight map based on the restored pulse train and perform curvature correction based on the first geometric prior parameters to obtain a defect distribution map.

[0053] Compared with existing technologies, this application obtains waveforms in the reference region and calculates the point spread function by combining geometric and material parameters. Then, it introduces time-axis normalization mapping, ensuring that the system response based on deconvolution remains consistent with actual propagation characteristics under complex conditions such as curved hot-melt joints and gradual radial refractive index changes. This approach avoids the model mismatch problem caused by point spread function drift with position in existing methods, effectively improving the stability of echo separation and reducing artifacts caused by signal overlap. Simultaneously, the introduction of time-axis normalization transforms the originally non-stationary convolution process into an approximately stationary process, reducing the risk of noise amplification and providing a robust foundation for signal processing in high-absorption polyethylene materials.

[0054] Furthermore, after obtaining the restored pulse train, this invention performs curvature correction by combining geometric prior parameters, accurately mapping the signal processing results into a time-of-flight map and generating a defect distribution map. This process maintains consistency between temporal data and spatial structure in complex curved surfaces, enabling reliable characterization of defect location and morphology. Compared to existing general processing methods based on flat or simplified layered models, this invention can output detection results consistent with the actual geometry in real-world hot-melt joint scenarios, significantly improving defect location accuracy and imaging reliability, and possessing outstanding industrial application value. Attached Figure Description

[0055] Figure 1 A flowchart of a terahertz signal aliasing cancellation method based on deconvolution provided in this application embodiment;

[0056] Figure 2 A flowchart illustrating a joint update method provided in this application embodiment;

[0057] Figure 3 A flowchart illustrating a method for performing deconvolution processing is provided in this application embodiment;

[0058] Figure 4 This is a schematic diagram of a terahertz signal aliasing cancellation system based on deconvolution, provided as an embodiment of this application. Detailed Implementation

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

[0060] See Figure 1 The diagram shows a flowchart of a terahertz signal aliasing cancellation method based on deconvolution provided in this application embodiment, including steps S101 to S106, wherein:

[0061] S101: Acquire time-domain waveforms of the reference area and the detection area to obtain the first reference waveform and the first observed waveform; obtain geometric and material parameters to obtain the first geometric prior parameters;

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

[0063] S103: Calculate the time axis normalized mapping function based on the first geometric prior parameters to obtain the first time normalized mapping function;

[0064] S104: Apply the first time normalization mapping function to the first observed waveform to obtain the first normalized observed waveform;

[0065] S105: Perform deconvolution processing based on the first normalized observation waveform to obtain the restored pulse train;

[0066] S106: Generate a time-of-flight map based on the restored pulse train, and perform curvature correction based on the first geometric prior parameters to obtain a defect distribution map.

[0067] Regarding the above S101:

[0068] In this embodiment, a terahertz time-domain spectroscopy detection system is used to inspect the thermofusion joint of a polyethylene pipe. First, a reference signal is acquired in a uniform, defect-free area on the outer wall of the joint to obtain the reference time-domain waveform under stable conditions, denoted as the first reference waveform. Simultaneously, a corresponding time-domain waveform is acquired in the detection area of ​​the joint under test, denoted as the first observed waveform. All of these waveforms are emitted by a terahertz pulse source and received by a detector, outputting as time-series data.

[0069] Based on this, geometric and material parameters describing the physical properties of the joint are obtained as the first geometric prior parameters. The geometric parameters may include the radius of curvature of the molten bead, the wall thickness, the pipe diameter, and the detection path angle; the material parameters may include the refractive index, dispersion characteristics, and absorption coefficient of polyethylene. These parameters can be obtained through external profile measurement, mechanical measurement, or a combination of material database lookup and actual measurement. The obtained geometric and material parameters serve 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] Regarding S102 above:

[0071] In one implementation, the point spread function is calculated based on a first reference waveform acquired in the reference region and the aforementioned first geometric prior parameters. Specifically, the original response signal of the system under uniform and defect-free conditions is first characterized using the time-domain waveform of the reference region. Since the reference region does not contain structural defects, its echo signal 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. These geometric prior parameters include the radius of curvature of the molten bead, pipe diameter, wall thickness information, and incident path angle. The material parameters include the refractive index curve, dispersion relation, and absorption coefficient of polyethylene within the operating frequency band. The optical path lengths for different propagation paths are determined using the geometric parameters and converted into time delays to correct the time delay distribution of the reference waveform. Phase expansion is performed on the frequency domain signal using the material parameters to compensate for group delay distortion caused by dispersion, while an exponential attenuation correction is applied to the amplitude component to eliminate energy loss caused by medium absorption.

[0073] After the above corrections are completed, the obtained response waveform is normalized, and a minimum phase constraint is introduced to ensure that the generated system pulses meet the requirements of causality and energy concentration. The waveform obtained through this process serves as the first point diffusion function, which can reflect the system response characteristics of the thermofusion joint under actual geometric and material conditions, and serves as the kernel for subsequent deconvolution steps.

[0074] Regarding the above S103:

[0075] In one implementation, the time axis normalized mapping function is calculated using a first geometric prior parameter as input. Because the thermofusion joint has gradually changing curvature and radial refractive index, the propagation path length and group delay at different positions on the same scan line are no longer linearly related. Directly using the original time axis would introduce position-dependent distortion. Therefore, this embodiment establishes a geometric model that combines the radius of curvature, wall thickness, and incident angle information from the geometric prior parameters to calculate the equivalent optical path length corresponding to each sampling point. The equivalent optical path length is then converted into group delay via light speed and refractive index, thus obtaining the time delay distribution of each sampling point.

[0076] Based on the above, the time delay distribution is compared with an ideal linear time axis to generate a one-to-one mapping relationship. This mapping relationship can reparameterize the non-uniform sampling points on the original time axis into uniform sampling points on the normalized time axis, ensuring that 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 standard sampling points in the normalized domain, and signal reconstruction can be achieved using interpolation or resampling methods. The resulting normalized mapping function is defined as the first-time normalized mapping function, which eliminates the time-axis nonlinearity caused by the gradual changes in curvature and refractive index, enabling the deconvolution process to be performed under approximately stationary time-domain conditions.

[0078] Regarding S104 above:

[0079] In one implementation, the first observed waveform obtained in the detection area is mapped to a normalized time axis to obtain a first normalized observed waveform.

[0080] First, the aforementioned first-time normalization mapping function is invoked to establish a one-to-one mapping table between the original time coordinates and the normalized time coordinates. This mapping table uses the sampling point number or the original sampling time as an index to output the normalized target time. Preferably, it can be stored as an array or matrix using a lookup table method for quick subsequent retrieval.

[0081] Since the original sampling points often fall in non-uniform locations within the normalized domain, numerical interpolation methods are needed to avoid distortion. Specifically, linear interpolation, spline interpolation, or resampling algorithms based on Fourier interpolation can be used. In this embodiment, cubic spline interpolation is preferred to smoothly reconstruct the waveform. The reconstruction of the sampling points can be achieved using common numerical computation tools, such as MATLAB, Python's SciPy.interpolate, or the LabVIEW interpolation library.

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

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

[0084] Regarding the above S105:

[0085] In one implementation, deconvolution processing is performed on the first normalized observation waveform to obtain the corresponding restored pulse train. Specifically, the aforementioned first point spread function is first loaded into the deconvolution operation module as the system convolution kernel. The first point spread function has been established by combining the reference region waveform, the radius of curvature of the molten bead, the pipe wall thickness, and the refractive index characteristics of the material, and can reflect the propagation characteristics of the detection system under the geometric conditions of the fusion joint.

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

[0087] Considering that the high absorption characteristics of polyethylene material can easily cause noise amplification during the inversion process, this embodiment introduces sparsity constraints and convolution residual constraints during the inversion process. The sparsity constraint is used to limit the restored pulse train to exhibit an isolated pulse structure in the time domain, to conform to the physical characteristics of the echo event. For example, L1 norm regularization can be applied to the restored result, or soft thresholding can be performed during iteration to suppress low-amplitude pseudo-pulses. The convolution residual constraint is used to ensure that the difference between the restored pulse train after convolution with the first point spread function and the first normalized observed waveform remains 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 passing condition.

[0088] In terms of solution methods, an iterative inversion strategy can be adopted. For example, the Iterative Shrink Thresholding Algorithm (ISTA) or the Alternating Direction Multiplier Method (ADMM) can be used as optimization tools. In each iteration, the convolution residual is calculated based on the current restored pulse train, and a sparsity operator is applied to the pulse train to update and obtain a new estimate. If the residual change is lower than a set threshold in several consecutive iterations, or the number of iterations reaches the upper limit, such as 100, the calculation is terminated, and the final restored pulse train is output.

[0089] The resulting restored pulse train appears in the time domain as a series of independent pulses after dealiasing, and its pulse position and amplitude can correspond to the actual path of multiple reflections in the actual thermofusion joint. For example, when there is an incomplete fusion defect in the detection area, an additional pulse with a shorter delay than the normal wall thickness will appear in the restored pulse train, thus providing basic data for the subsequent generation of time-of-flight maps and curvature correction.

[0090] Regarding S106 above:

[0091] First, each pulse event in the restored 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 the time of flight is calculated using the sampling interval. 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. For example, grayscale values ​​or pseudo-color mapping can be used to represent the magnitude of the time of flight, allowing echoes at different depth positions to be visually presented in the image.

[0092] Subsequently, curvature correction is performed on the time-of-flight map using the first geometric prior parameters. Because the molten beads in the thermofused joint have curved geometric features, directly using the uncorrected time-of-flight map will cause the same physical depth to appear as a curved trajectory on the image, affecting the accurate interpretation of the defect location. Therefore, this embodiment calculates the equivalent propagation path length at each scanning position based on the molten bead's radius of curvature and wall thickness information, and remaps the original time-of-flight data to the geometrically corrected coordinate system. For example, a polar coordinate model can be established to unfold the time-of-flight data affected by curvature into equal-depth layers in Cartesian coordinates, ensuring that the defect's position in the image matches its actual geometric position.

[0093] The time-of-flight map after curvature correction is defined as a defect distribution map, which can realistically reflect the defect situation in the detection area spatially. For example, when there is an incomplete fusion defect, the corrected defect distribution map will show local abnormal depth values ​​or energy voids in the corresponding area; when there is porosity, it will appear as point-like or patchy abnormal signals in the image. In this way, this embodiment not only completes the dealiasing in the signal domain, but also accurately projects the processing results onto the actual geometric structure, providing a reliable basis for defect detection and localization.

[0094] Optional See also Figure 2 This is a flowchart of a joint update method provided in an embodiment of this application. In another embodiment, in order to examine and constrain the non-commutativity that may arise between "time axis normalization processing" and "deconvolution processing" in scenarios with curved surfaces and gradually changing refractive indices, this embodiment... Figure 1 Based on the illustrated process, a dual-path consistency check and joint update mechanism is introduced, specifically including steps S201 to S204, wherein:

[0095] S201: Perform deconvolution processing on the first observed waveform on the unnormalized original time axis to obtain the second restored pulse train;

[0096] S202: Apply the first time normalization mapping function to the second restored pulse train to obtain the second normalized restored pulse train;

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

[0098] S204: Based on the first coupling constraint value, the convolution residual between the first normalized observation waveform and the first point spread function, and the sparsity index of the restored pulse train, perform a joint iterative update to obtain the updated first time normalized mapping function, the updated first point spread function, and the updated restored pulse train.

[0099] Regarding the above S201:

[0100] Without time axis normalization, a deconvolution is performed on the first observed waveform to obtain the second restored pulse train. Specifically, the aforementioned first-point spread function is used as the convolution kernel to establish an inversion model on the original time axis with the goal of "sparseness constraint + convolution residual constraint". For example, the iterative shrinkage threshold method or the alternating direction multiplier method can be used to solve the problem. In each iteration, the convolution residual is calculated and a sparsity operator is applied to suppress low-amplitude pseudo-pulses, while keeping the causal and minimum phase constraints of the first-point spread function unchanged.

[0101] Preferably, before proceeding to the solution, the first observed waveform can be bandpass filtered and amplitude normalized to ensure that the input and convolution kernel have the same bandwidth and dimensions. The result of this step preserves the time delay details before curvature unfolds in the original time coordinate system, serving as a benchmark for subsequent path consistency verification.

[0102] Regarding the above S202:

[0103] The second restored pulse train obtained in S201 is input into the aforementioned first time normalization mapping function and mapped to the normalized time axis to obtain the second normalized restored pulse train. To avoid distortion caused by sampling misalignment, spline interpolation resampling is preferred. To eliminate the deviation between the overall amplitude and micro-delay, amplitude normalization and cross-correlation alignment can be performed on the main peak windows of the two pulse trains (correcting the whole-segment offset by sample-level micro-shifting). The mapping table can adopt a lookup table structure (array or matrix) for quick access; interpolation and cross-correlation can be implemented using common numerical libraries such as MATLAB, Python (SciPy.interpolate, NumPy.fft), or LabVIEW.

[0104] Regarding the above S203:

[0105] Under a unified normalized time axis and the same scan position index, the difference between the restored pulse train and the second normalized restored pulse train obtained in S202 is measured, and the first coupling constraint value is formed accordingly. To enhance physical relevance, the difference measurement may include the following complementary indicators:

[0106] Amplitude envelope difference: Within a preset window, the mean square error or correlation coefficient of the envelopes of two signals is measured to characterize the difference in energy distribution;

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

[0108] Morphological consistency: Compare the statistics of peak number, peak width and peak spacing to avoid false consistency caused by "peak splitting / merging".

[0109] The aforementioned multidimensional differences are aggregated into a single scalar according to preset weights, defined as an commutative difference value. For subsequent location-based updates, the difference value, along with its corresponding scan location and time window index, is archived together to form an indexed representation of the first coupling constraint value. For example, the weights can be tuned offline using a validation dataset, and the window can be a local interval of the main peak ± a number of sampling points.

[0110] Regarding the above S204:

[0111] Using the first coupling constraint value, convolutional residual, and sparsity index as joint objectives, the time-normalized mapping function, point spread function, and restored pulse train are alternately updated, outputting the updated first time-normalized mapping function, the updated first point spread function, and the updated restored pulse train. For example, to ensure physical feasibility, this embodiment applies the following boundary conditions to each updated object:

[0112] For the time-normalized mapping function, only small adjustments are allowed at key nodes, and the mapping function is kept monotonic and smooth. For example, monotonic splines or piecewise linear functions are used, and the node displacement in each round is limited to no more than a preset number of pixels.

[0113] For the point spread function, only fine-tuning of the amplitude attenuation parameter and phase is allowed to maintain causality and minimum phase constraint, preventing the generation of non-physical long tails or pre-ringing.

[0114] For the restored pulse train, the iterative solution with sparsification and residual dual constraints is continued, and the decrease in commutativity difference is used as an additional convergence signal.

[0115] In practical implementation, the update order can be either "pulse train - mapping - point spread function" or "pulse train - point spread function - mapping", as long as the joint decrease of the first coupling constraint value and the convolution residual is used as the iterative advancement condition after each round of iteration. 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 stopping condition will be further explained in subsequent implementations.

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

[0117] Through the dual-path consistency and joint update mechanism of S201 to S204, without changing the hardware acquisition process, the model mismatch between time normalization mapping and point spread function is dynamically corrected by using data-driven coupling constraints. This suppresses the systematic error caused by the non-commutativity of geometry-convolution in curved surface scenarios, thereby improving the consistency and interpretability of the restored pulse train in the normalization domain.

[0118] Optional, also includes:

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

[0120] The joint iterative update is terminated when the first stopping condition is met, and the target recovery pulse train is output as the recovery pulse train.

[0121] In another embodiment, to avoid information depletion-type "false convergence" or artifacts caused by over-iteration in the joint iterative update, this embodiment... Figure 1 Based on the process and S201 to S204, a shutdown criterion is set.

[0122] In practical implementation, based on the updated restored pulse train output by S204, and combined with the updated time normalization mapping function and point spread function for the corresponding round, the amplitude of the restored pulse train is normalized and aligned with microsecond delays on the normalized time axis. For example, cross-correlation alignment can be performed within a local window of the main peak ± several sampling points, and the high absorption frequency band can be weighted by sub-band or bandpass filtering to output a standardized dataset for index calculation.

[0123] Within the local window, for the restored pulse train and its convolutional residual in two adjacent iterations (the kth iteration and the (k-1th iteration), the envelope gradient change rate and residual energy reduction are calculated respectively, and the change of the non-zero support set is statistically analyzed, for example, by the set Hamming distance or intersection-union ratio.

[0124] For example, when the envelope gradient change rate is less than 1% to 3% and the residual energy decrease is less than a preset threshold, the first stability of the scan position is recorded as met; all scan positions are aggregated to obtain a global score of the first stability index.

[0125] Under the same normalized time axis and scan position index, the consistency with the geometric prior is evaluated based on the updated recovered pulse train. Specifically, this may include: performing a monotonicity test on the delay order of the main echo groups of adjacent scan positions, and performing a consistency test on the polarity and phase continuity of the recovered pulses, for example, determining whether a phase flip occurs as a non-physical jump by analyzing the consistency of the signal phase or sign.

[0126] For example, a higher weight is assigned to the phase consistency term in the high absorption frequency band to mitigate the impact of amplitude attenuation on the judgment; when the monotonicity test passes and the phase / polarity consistency rate reaches a preset ratio, the second stability of that scan position is recorded as having met the standard. All scan positions are aggregated to obtain a global score for the second stability index.

[0127] The first stability index and the second stability index are combined to form the first stopping condition.

[0128] For example, a "dual threshold + hysteresis" rule can be adopted: when the first stability index is higher than threshold A and the second stability index is higher than threshold B, and this condition is met for K consecutive iterations (e.g., K is 3 to 5), and the proportion of qualified scan positions is not lower than a preset proportion (e.g., not less than 90%), then the first stopping condition is considered met. To avoid jitter caused by local noise, local masks can be set for individual scan positions that do not meet the standards, and their inclusion in the global judgment can be delayed.

[0129] When the first stopping condition is met, the joint iterative update described in S204 is terminated, the updated restored pulse train of the current round is output as the target restored pulse train, and the corresponding time normalization mapping function and point spread function are frozen for subsequent imaging use; if the first stopping condition is not met, S204 is returned to enter the next round of joint update. For example, a minimum and a maximum number of iteration rounds can be set simultaneously, for example, no less than 10 rounds and no more than 200 rounds. If the first stopping condition is not met even after reaching the maximum number of rounds, the final freeze is performed based on the principle of prioritizing the second stability index to ensure that physical consistency is not sacrificed.

[0130] In this way, determining the timing of stopping the joint iteration from the two complementary dimensions of numerical stability and physical consistency avoids "false convergence" caused by residual convergence alone, and suppresses non-physical artifacts caused by over-iteration, thereby ensuring that the output target restoration pulse train is interpretable and reproducible in the terahertz detection scenario with curved surfaces and gradually changing refractive indices.

[0131] Optional, see Figure 3 A flowchart of a deconvolution processing method provided in this application embodiment includes steps S301 to S306, wherein:

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

[0133] S302: Calculate the set of diffusion functions for direction-related points based on the first reference waveform and the first direction set to obtain the first direction kernel cluster;

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

[0135] S304: Based on the first direction up-dimensional data and the first direction kernel cluster, perform group sparse deconvolution processing to obtain the restored pulse tensor of the direction index;

[0136] S305: Perform directional neighborhood smoothing on the restored pulse tensor of the directional index to obtain a directional continuous restored pulse tensor;

[0137] S306: Project the continuous recovery pulse tensor of the direction onto the time axis and the scan axis according to the first direction set to obtain the recovery pulse train.

[0138] First, in the description of this embodiment, Indicates the scan position index or the horizontal coordinate along the scan path, used to identify spatial sampling points within the same scan line or between adjacent scan 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, it is preferable to apply neighborhood median or spline smoothing to the normal field, and obviously abnormal abrupt changes are removed or backfilled according to threshold rules.

[0142] Regarding S302:

[0143] The first set of directions is determined based on the first geometric prior parameters and the normal vector table. Specifically, discrete sampling is performed within the effective angular domain related to incident / reflection, adaptively increasing the direction sampling in sections with large curvature changes, and appropriately sparsely sampling in regions with gentle curvature, to form a direction list with sufficient coverage and controllable redundancy.

[0144] For example, the number of directions can be 12 to 24, and the direction step is adaptively adjusted according to the local radius of curvature. The set of directions is stored in the form of an ordered list or lookup table index and associated with the scan position index for sharing and use in subsequent steps.

[0145] Regarding S303:

[0146] Based on the first reference waveform, the first direction set, and the first geometric prior parameters, the set of diffusion functions of direction-related points is calculated to obtain the first direction kernel cluster.

[0147] Specifically, based on the first-point diffusion function, the reference response is subjected to directional consistency transformation according to the possible propagation path differences and phase delays in each direction, and amplitude normalization and bandwidth consistency correction are performed to make kernels in different directions comparable. For example, this can be pre-calculated and cached in a MATLAB or Python environment. This is to support subsequent fast convolution / deconvolution operations.

[0148] Regarding S304:

[0149] Based on the first normalized observation waveform and the first direction set, directional filtering is performed to construct the first direction up-dimensional data. Specifically, the first normalized observation waveform is applied to various... The corresponding steerable filters or directional time-frequency atoms (e.g., orientation-aligned Gabor atoms) extract the directional response within a frequency band matching the bandwidth of the first directional cluster, outputting a form such as The data tensor is upgraded in directional dimensions. To avoid boundary effects, mirroring or zero-padding boundary processing is preferred; to maintain comparability between channels, energy normalization is performed on each directional channel.

[0150] Regarding S305:

[0151] Based on the first-direction up-dimensional data and the first-direction kernel cluster, sparse deconvolution processing is performed to obtain the restored pulse tensor of the direction index, and neighborhood smoothing is performed in the direction dimension. Specifically, with the same Location spans all The amplitude vectors of the channels are treated as a "group," and a group sparse solution strategy is employed to activate them together in a small number of adjacent directions, in order to conform to the physical property that the direction gradually changes with position. For example, group sparse solution can be implemented within the ADMM or Iterative Shrinking Threshold (ISTA) framework; to suppress non-physical direction jumps, further... Applying neighborhood smoothing or anisotropic regularization (e.g., weighted averaging of adjacent directional channels) to the output direction continuously restores the pulse tensor In practice, directional channels below the energy threshold can be zeroed out to improve robustness and computational efficiency.

[0152] Regarding S306:

[0153] The directional continuous recovery pulse tensor is projected onto the time axis and scan axis according to the first direction set to obtain the recovery pulse train. Specifically, based on the directional response intensity, consistency with the normal, or preset directional weights, the pulse train is... The dimensions are weighted and converged, and the energy of each directional kernel is normalized before projection to avoid amplitude bias caused by directional bias.

[0154] For example, it is possible in each The position is weighted and summed over adjacent directional channels and then thresholded to obtain a sparse pulse representation in the time domain; the restored pulse train obtained through this projection step is then returned to... This domain seamlessly integrates with the subsequent imaging and curvature correction steps in the above process.

[0155] Through the processing chain of S301 to S306 described above, this embodiment explicitly incorporates the directional element into the signal characterization and deconvolution solution process. In the directional dimension, the directional support of the echo is constrained by group sparsity and neighborhood smoothing, which suppresses the temporal artifacts caused by directional mismatch and multi-angle mixing within pixels, thereby improving the physical interpretability and spatial consistency of the restored pulse train.

[0156] Optionally, determining the direction set based on the first geometric prior parameters and the normal information includes:

[0157] An initial direction set is generated based on the first geometric prior parameters and the normal information to obtain the first initial direction set;

[0158] The first direction indication is obtained by calculating the direction indication within a preset time window and a preset scanning window based on the first normalized observation waveform.

[0159] The first direction consistency score is obtained by calculating the direction consistency score based on the first direction indication and the first initial direction set;

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

[0161] In another embodiment, the process of determining the direction set based on the first geometric prior parameters and normal information includes four steps: adaptive candidate generation, data verification, score aggregation, and constrained update, so as to form a direction set consistent with the local propagation direction under the conditions of surface and gradual change of refractive index.

[0162] First, based on the first geometric prior parameters such as the radius of curvature of the molten bead, the wall thickness, and the incident path angle, and combined with the registered normal vector field, a feasible incident / reflection angle domain is determined for each scanning position, and a first initial direction set is discretized within this domain. To adapt to different curvature regions, this embodiment preferably sets the discretization density adaptively according to the local curvature: smaller angle steps are taken where the curvature is large, and the steps are appropriately widened where the curvature is gentle. The direction set is stored in an ordered list bound to the scanning position index for subsequent retrieval.

[0163] Subsequently, a time window related to the main echo and a scanning window centered on the current scanning position are selected from the first normalized observation waveform to calculate the first direction indication quantity used to reflect the "data evidence" of each candidate direction.

[0164] For example, a steerable filter or directional time-frequency atom aligned with each candidate direction can be applied to obtain the directional response energy, and the energy can be averaged or the peak value can be taken as the energy component within a window; at the same time, the consistency between the predicted phase symbol and the measured phase symbol corresponding to the candidate direction can be compared based on the phase of the analytical signal to obtain the phase consistency component; if necessary, the direction-guided cross-correlation peak value can be introduced as a morphological component to measure the reachability of echo alignment under the assumption of that direction.

[0165] The above components are normalized in magnitude and scale, forming a set of first direction indicator vectors indexed by direction.

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

[0167] After obtaining the directional indication, a directional consistency scoring model combining geometric priors and data evidence is established to score each of the first initial directional sets.

[0168] For example, the energy component, phase consistency component and morphology component can be weighted and aggregated into a first directional consistency score of 0 to 1 according to a preset weight; the weight can be empirically adjusted through offline validation set or production sample, 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 scoring, neighborhood smoothing or robust aggregation can be performed on scores in the same direction within the scanning window, thereby reducing the impact of local isolated noise.

[0170] After the score is obtained, the first initial direction set is shrunk and resampled according to the angular spacing constraint defined by the first geometric prior parameter to obtain the second direction set.

[0171] Specifically, when the distance between adjacent candidate directions is less than the minimum distance, they are weighted and merged according to their consistency scores to avoid redundancy and overfitting caused by excessively dense dispersion; when the distance between adjacent candidate directions is greater than the maximum distance and both ends have high scores, a new direction is inserted at the midpoint or the point where the score gradient is maximum to avoid missed detections caused by excessively sparse dispersion.

[0172] To prevent non-physical jumps in the orientation set between adjacent scan positions, this embodiment limits the addition / deletion ratio of each update to no more than a preset upper limit of the initial set size, such as no more than 30%, and sets an upper limit constraint on the set differences between adjacent positions, such as the total angular offset of the set differences not exceeding a certain number of discrete levels.

[0173] For example, the minimum angular spacing can be 2° to 5°, and the maximum angular spacing can be 8° to 15°; the energy threshold can be set to 10% to 30% of the global maximum value, and the phase consistency rate threshold can be set to 70% to 90%; the above values ​​are exemplary ranges, and can be adjusted according to the bandwidth of the detection system and the absorption characteristics of the material.

[0174] At the implementation level, the direction indication calculation and score aggregation can be completed using steerable filtering, Hilbert analytical phase and cross-correlation library functions in MATLAB or Python environments, and a lookup table structure is used to maintain a bidirectional index of the direction set and scan position. The second direction set obtained through the above-mentioned constrained update serves as the adaptive direction basis for subsequent direction kernel cluster construction, direction dimensionality increase and group sparse deconvolution, thereby simultaneously satisfying the requirements of coverage and smoothness in the direction dimension, and suppressing the direction bias caused by normal measurement error or absorption.

[0175] Optional, also includes:

[0176] Based on the first initial direction set and the second direction set, the change in the set is calculated to obtain the stability condition of the first direction;

[0177] 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 used as a new initial direction set and the direction set generation steps are repeated until the first direction stability condition is met.

[0178] In another embodiment, to suppress round-trip jitter and over-updates in regions with low signal-to-noise ratio or curvature abrupt changes, and to determine the convergence state of the direction set while satisfying physical consistency, this embodiment calculates the set change based on the first initial direction set of the previous round and the second direction set of the current round after completing the adaptive generation of the direction set, thereby forming the first direction stability condition; when the stability condition is satisfied, the current direction set is frozen, otherwise the second direction set is used as the new initial set and the generation is repeated until the stability condition is satisfied.

[0179] Specifically, the two sets of directions are first paired under the same scan position index. To avoid mismatches, bidirectional nearest neighbor matching in ascending order of angle is used, and the allowed angle difference is set to 1–2°. Elements exceeding this range are no longer considered to be in the same direction. Thus, elements that were not matched in the previous set are recorded as deletion candidates, and newly appearing elements in the current set are recorded as new candidates. A source mapping is established for the new directions obtained by merging in the previous round, so as to statistically analyze angle drift and convergence trends later.

[0180] After pairing is completed, the set change at the scan position is calculated. The set change is a multi-component weighted aggregation indicator, which includes at least: the proportion of the number of additions and deletions to the size of the previous set, used to measure the size change; the average and maximum angular drift between matched directions, used to reflect the overall translation of directions; the overall separation of the two discrete sets on the angular axis, for which the maximum nearest neighbor difference process is used, that is, the angular difference from each direction in set one to the nearest neighbor in set two is calculated and the maximum value is taken, and then the larger of the result and the reverse calculation result is taken, and normalized to the full angular domain; the change in the mean and dispersion (such as variance or quantile difference) of the directional consistency score, used to reflect the stability of the data evidence; and the spatial smoothness across positions, for which the 90th quantile of the angular offset after pairing is statistically analyzed in the neighborhood centered on the current scan position (e.g., ±2–5 positions) and normalized to suppress local jumps.

[0181] Considering the low amplitude reliability of terahertz frequencies in the high absorption band, this embodiment increases the weight of components related to phase tag stability and group delay monotonicity during aggregation, for example, to about 0.3–0.5, and applies subband weighting or out-of-band suppression to energy-related components to enhance the ability to identify physical consistency.

[0182] To avoid interference from low signal-to-noise ratio (SNR) locations, a mask is constructed that is not included in the stability assessment. Specifically, the noise envelope is estimated in the baseline region outside the main peak, and the ratio of the main peak energy to the noise energy is calculated. When this ratio is below 6–10 dB, the scan location is marked as a mask point. The ensemble change of the mask point is only used for logging and is not included in global aggregation and stability assessment. Subsequently, a joint global and local assessment is performed to form the first directional stability condition: when the global ensemble change is below the global change threshold, and the proportion of locations whose ensemble changes at each scan location are below the local change threshold is not lower than the lower limit of spatial proportion (e.g., 80%–90%), and this condition is met for 2–4 consecutive updates, stability is considered achieved.

[0183] To avoid jitter caused by repeated crossings near the threshold, a consistent hysteresis is set for entry / exit: after entering stability, if slight fluctuations occur, a stricter set of exit thresholds is used for exit judgment, that is, the global and local thresholds are each increased by about 10%, the lower limit of the space proportion is reduced by about 5 percentage points, and exit is performed if the conditions are not met in one round, thereby improving the robustness of the decision.

[0184] When the stability condition of the first direction is met, the second direction set of the current round is frozen as the new first direction set 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 used as the new initial direction set, and the direction set generation process is returned to continue iterating. To suppress update oscillations and control computational overhead, a step size limit is imposed on each round of updates: the net increase / decrease ratio is no higher than 20%–30% of the size of the previous round set, the maximum allowable angle drift in a single iteration is no more than 2–4°, and the total number of directions is kept within a preset range, for example, no less than 8–12 and no more than 1.5 times the initial set. In regions with significant curvature abrupt changes, the above upper limits can be tightened accordingly to prevent non-physical jumps.

[0185] At the implementation level, set pairing, angle difference quantiles, and set distance can be calculated in common numerical environments, and the analytical phase can be obtained through Hilbert transform. Each round records an update log, including the net number of additions and deletions, average and maximum angle drift, set distance, score change, spatial consistency, mask ratio, and whether the set is frozen and the threshold group used, for subsequent auditing and parameter tuning.

[0186] Through the above control, the adaptive direction set can achieve good convergence and reproducibility in the terahertz detection scenario with curved surfaces and gradually changing refractive indices, thereby improving the reliability of subsequent deconvolution and imaging results.

[0187] Optionally, the step of performing group sparse deconvolution processing based on the first direction up-dimensional data and the first direction kernel cluster to obtain the restored pulse tensor of the direction index includes:

[0188] Based on the first geometric prior parameters and the normal information, the directional adjacency relationship is determined, and a first directional compatibility matrix is ​​generated;

[0189] In the group sparse solution, a shared constraint is applied to the non-zero support of the direction dimension based on the first direction compatibility matrix, so that the non-zero direction corresponding to the sampling at the same time is limited to the adjacent direction specified by the first direction compatibility matrix, and the restored pulse tensor of the direction index is output.

[0190] In one implementation, in order to suppress non-physical discrete jumps in the direction dimension and ensure that the directional activation at the sampling point at the same time is consistent with the geometric continuity of the detection surface, this embodiment first determines the directional adjacency relationship based on the first geometric prior parameter and normal information, and generates a first direction compatibility matrix for constraint solving.

[0191] Specifically, for each spatial sampling point on each scan line, using the previously registered normal vector and the corresponding candidate direction set, it is determined whether any two directions are allowed to be considered "adjacent" according to the principles of angular adjacency and geometric continuity.

[0192] In implementation, an adjacency radius matching the directional discrete step size can be set. Directional pairs whose angular difference is within this radius are considered adjacent. Directional pairs that maintain a smooth change across adjacent scan positions are also marked as "adjacent." Directional pairs that do not meet these conditions are marked as "non-adjacent." These adjacency relationships are recorded in a sparse binary manner, forming a first directional compatibility matrix. Each row and column corresponds to a discrete directional element. Allowed values ​​in the matrix indicate that the direction and another direction can co-occur in the solution, while disallowed values ​​indicate that they should not be activated simultaneously at the same sampling point.

[0193] To reduce the impact of noise on adjacency relationships, it is preferable to smooth the normal field before constructing the matrix and to use mirroring or extension methods at the boundaries of the direction set to avoid incorrectly eliminating boundary directions due to insufficient neighbors.

[0194] After obtaining the first direction compatibility matrix, a shared constraint is introduced in the group sparse deconvolution solution process: the non-zero directions allowed to appear at the same sampling point should form a connected subset connected by the "adjacent" relationship.

[0195] This embodiment employs an iterative strategy of "estimation, projection, and refinement": First, without imposing orientation compatibility constraints, a preliminary group sparse deconvolution is performed on the orientation-upgraded data and orientation kernel clusters to obtain an amplitude estimate of the orientation index. Then, based on the first orientation compatibility matrix, a "feasibility projection" is performed on this amplitude estimate, i.e., the orientation activations sampled at the same time are filtered, retaining only adjacent orientations connected to the current main activation orientation, while suppressing or merging sporadic activations falling on non-adjacent relationships into their nearest adjacent connected components. Subsequently, the convolution residual is recalculated on the new feasible support, and the amplitude is updated. The above "estimation-projection-refinement" steps are performed alternately until the change in residual is lower than a preset threshold or the number of iterations reaches its upper limit.

[0196] In practical implementation, "feasible projection" can be specifically implemented as a combination of two types of operations: one is support pruning based on the compatibility matrix, which directly sets isolated directions that do not meet the adjacency condition to zero; the other is amplitude redistribution based on adjacent directions, which smooths the amplitude within a small set of adjacent directions according to energy and phase consistency to reduce directional jaggedness caused by noise. To ensure numerical stability, it is preferable to perform a light neighborhood smoothing on the direction dimension at the end of each iteration, while maintaining the sparsity of the time dimension.

[0197] Regarding parameter selection, to ensure that the constraints are neither too strict and swallow up the real multi-path directions, nor too loose and allow non-physical activation, this embodiment sets the directional adjacency radius to a range of the same order of magnitude as the local angular discrete interval, preferably 1 to 2 times the local angular discrete interval; the angular discrete interval is the angle between two adjacent discrete directions; in the adaptive direction set, the angular discrete interval is dynamically determined according to the scan position. The allowed adjacent depths at the same sampling point can be one or two discrete levels; the minimum amplitude threshold and iteration stop threshold triggered by "feasible projection" can be tuned according to the system noise level, for example, terminating when the convolution residual decreases by less than a few percent in several consecutive iterations.

[0198] For example, the orientation compatibility matrix can be pre-generated and cached as a sparse adjacency structure; group sparse solution and feasible projection can be implemented in common numerical computing environments; to avoid the boundary orientation being systematically suppressed, energy normalization can be performed on each directional channel before projection to balance the response differences of different directional kernels.

[0199] In this way, the final output direction index restoration pulse tensor presents a connected and smooth support distribution in the direction dimension and maintains a sparse pulse shape in the time dimension. The result is used for subsequent direction projection and imaging processes, which can effectively suppress time domain artifacts caused by direction mismatch in curved surfaces and scenarios with gradually changing refractive indices.

[0200] Optionally, the step of determining the directional adjacency relationship based on the first geometric prior parameters and the normal information, and generating the first directional compatibility matrix, includes:

[0201] Based on the first geometric prior parameters and the normal information, the equivalent incident angle corresponding to each direction is determined, and the phase characteristics of the reflection coefficient corresponding to the equivalent incident angle are evaluated in the terahertz working frequency band. Based on the in-band phase sign consistency, mutually exclusive direction pairs are marked to obtain the first incompatible set.

[0202] After excluding the first incompatible set, the group delay order of the main echoes at adjacent scanning positions is calculated based on the first geometric prior parameters, and the direction pairs are filtered based on the criterion of maintaining the group delay order to obtain the first compatible set.

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

[0204] In one implementation, in order to match the adjacency relationship of the directional dimension with the terahertz electromagnetic propagation characteristics, surface geometry features and refractive index gradient scene, this embodiment determines the directional adjacency relationship based on the first geometric prior parameter and normal information, and generates a directional compatibility matrix accordingly to constrain the sparse deconvolution solution process.

[0205] Specifically, for the scanning area of ​​the polyethylene hot-melt joint, an initial set of directions is determined based on the surface normal and geometric parameters at the scanning location, such as the radius of curvature of the molten bead and the pipe wall thickness. For example, at a typical scanning location, the initial set of directions calculated from the normal can be discretized into angles such as 0°, 5°, 10°, 15°, and 20°. These discrete angles represent the angles between the incident and reflected directions and the local normal, and serve as the initial set for subsequent analysis of directional compatibility.

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

[0207] For example, for angles of 0°, 5°, and 10°, the reflected phase signs calculated near the 0.5 THz frequency point may be positive, positive, and negative, respectively. This indicates that there is a stable phase opposition between angle 10° and 0° and 5°, that is, 10° and the former two are in a "phase-exclusive" relationship. This phase-exclusive relationship is recorded as the first incompatible set, which is used as a constraint for the subsequent construction of the direction compatibility matrix.

[0208] After excluding the first incompatible set, this embodiment further evaluates the order of the main echo group delays at adjacent scan positions. Specifically, for example, at scan positions x=100 and x=101, if the order of the main echo group delays corresponding to angles 5° and 10° in the direction set is "5° earlier than 10°", then the direction pair is considered to conform to the group delay monotonicity and belong to the direction compatible case; conversely, if the order is reversed at x=101, such as 10° earlier than 5°, it is recorded as not conforming to geometric continuity and is not compatible. Through this judgment, the direction pairs that conform to the group delay monotonicity are recorded as the first compatible set, which is used to construct the direction compatibility matrix.

[0209] Based on the first incompatible set and the first compatible set, a first direction compatibility matrix is ​​constructed. This matrix uses the initial direction set as an index, and each matrix element represents whether two directions are allowed to be activated simultaneously at the same sampling time. During matrix initialization, all elements are set to prohibit adjacency; subsequently, direction pairs belonging to the first compatible set are marked as allowed adjacency. For example, if angles 0° and 5° are allowed to be adjacent, the corresponding position in the matrix is ​​marked as allowed. For direction pairs belonging to the first incompatible set, such as 5° and 10° having a phase-exclusive relationship, they are marked as prohibited adjacency.

[0210] To ensure the robustness of the numerical solution, this embodiment further performs sparsification and symmetry processing after matrix generation to eliminate isolated allowed terms and ensure the consistency of the adjacency relationships of direction pairs. At the boundaries of the direction set, mirroring or extension methods can be used to supplement the minimum adjacency degree to avoid the boundary directions being systematically suppressed during projection due to a lack of adjacent channels.

[0211] During the group sparse deconvolution solution stage, non-zero directional activations at the same time sampling point are restricted to the allowed adjacency relationships specified by the first directional compatibility matrix. Specifically, the solution process begins with an initial group sparse solution without adjacency constraints, yielding initial directional activation estimates. Subsequently, the initial activations at each time sampling point are feasiblely projected according to the directional compatibility matrix, suppressing or merging isolated activations that do not satisfy the adjacency relationship into the nearest neighbor compatible direction. This process is repeated iteratively over multiple rounds until the convolution residual change falls below a set threshold (e.g., 1–3%) or reaches the maximum number of iterations, e.g., 50–100, ultimately yielding the restored pulse tensor of the directional index.

[0212] Through the above process of constructing and solving the directional compatibility matrix, in the terahertz detection scenario with curved surfaces and gradually changing refractive indices, the non-zero activation of the directional dimension is always confined to the set of directions that conform to geometric continuity and phase consistency during the deconvolution process, thereby effectively suppressing non-physical directional jumps and temporal artifacts, and improving the physical consistency and spatial positioning accuracy of the final restored pulse train.

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

[0214] The equivalent incident angle of each scanning position is determined based on the first geometric prior parameters, and the representative frequency of each restored pulse is determined based on the spectral energy centroid of the restored pulse train, thus obtaining the first representative frequency set.

[0215] Based on the equivalent incident angle and the first representative frequency set, the interface reflection phase label is calculated to obtain the first phase label set; based on the analytical signal phase of the restored pulse train at the pulse peak, the phase label is calculated and measured to obtain the second phase label set.

[0216] Phase consistency verification is performed based on the first phase tag set and the second phase tag set. When there is inconsistency, polarity correction and micro-delay correction are performed on the corresponding recovery pulse to obtain a phase-aligned recovery pulse train.

[0217] A time-of-flight map is generated based on the phase-aligned restored pulse train, and an angle-resolved group refractive index is determined based on the equivalent incident angle and the first representative frequency set. Curvature correction is then performed to obtain the defect distribution map.

[0218] In one embodiment, to avoid false depth and geometric distortion caused by terahertz phase reversal and dispersion, this embodiment introduces phase tag calibration and angle-resolved group refractive index mapping after obtaining the restored pulse train, and completes the construction of time-of-flight map and curvature correction accordingly, and outputs a defect distribution map.

[0219] In practice, firstly, the equivalent incident angle is determined for each scan line and each scan position based on the first geometric prior parameter. The equivalent incident angle is calculated based on the local normal of the detection surface and the geometric relationship between transmission and reception, and is associated with the scan coordinates point by point. Based on this, representative frequencies are extracted from each pulse of the recovery pulse train.

[0220] The specific approach is as follows: Set a time window near the pulse peak, for example, ±5 to 15 sampling points from the peak. Perform windowing and fast spectrum analysis on the signal within the window, and calculate the spectral energy centroid as the representative frequency of the pulse. When the system bandwidth is wide, several sub-bands can be divided within the effective frequency band and the energy centroid can be calculated separately. Finally, the representative frequency set is obtained by majority vote or weighted average.

[0221] Secondly, by combining the equivalent incident angle and the representative frequency set, the expected interface reflection phase labels are calculated to form the first phase label set. For this purpose, a Fresnel reflection model can be used, and the refractive index and absorption parameters of polyethylene from the material library can be called upon to evaluate the phase sign and transition characteristics of the reflection coefficient under different incident angles and polarization conditions within the working frequency band. When a sub-band strategy is adopted, the phase signs within each sub-band are statistically analyzed to achieve majority consensus. If necessary, the weight of frequency bands close to the phase transition interval is reduced to decrease edge uncertainty.

[0222] Simultaneously, measured phase information is obtained from the restored pulse train: the analytical signal phase is calculated at the peak of each pulse, and phase unwrapping is performed by combining neighborhood expansion to obtain the measured phase label corresponding to the time window, forming a second phase label set.

[0223] In addition, to enhance robustness, neighborhood consistency smoothing can be performed on the measured phase tags between adjacent scan positions, and the pulse can be set not to participate in tagging when the signal-to-noise ratio of the main peak is too low.

[0224] Next, consistency checks and corrections are performed on the two sets of phase labels. When the first and second phase labels are inconsistent, polarity correction is first performed on the corresponding pulses according to the polarity rules, i.e., the entire pulse is flipped; if a small phase residual still exists, micro-delay correction is performed while maintaining polarity. Micro-delay correction can be implemented in two ways: one is to perform sub-sampling precision interpolation repositioning on the pulse peak neighborhood in the time domain, such as parabolic interpolation or fractional delay interpolation based on window functions, to slightly shift 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 restored pulse train is obtained.

[0225] To avoid overcorrection, this 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 on the allowable time shift in a single instance, for example, no more than half a sampling interval.

[0226] Subsequently, a time-of-flight map is generated based on the phase-aligned reconstructed pulse train: the arrival times of the main echo and key secondary echoes are extracted at each scan position, stacked in scan order to form a set of two-dimensional time-of-flight profiles, and the time-of-flight magnitude is represented by grayscale or pseudo-color mapping. To stably convert the time of flight into geometric depth, this embodiment does not use a fixed refractive index, but instead uses an angle-resolved group refractive index for path length reconstruction.

[0227] The specific method is as follows: based on the equivalent incident 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 flight time is converted into the propagation distance in the medium point by point; when the sub-band strategy is adopted, 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] For example, in the commonly used frequency band of polyethylene, the value of the group refractive index may vary slowly with frequency in the range of about 1.47 to 1.55, and there will also be slight differences at different incident angles. These differences are compensated by the aforementioned angle-resolved mapping.

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

[0230] Polar coordinate expansion can be used to expand the uniformly deep layers along the surface to equidistant layers in Cartesian coordinates. For regions with rapid curvature changes, it is preferable to reduce the local expansion window and introduce neighborhood smoothing to reduce numerical instability. The image obtained after curvature correction is defined as a defect distribution map: unfused defects usually appear as local abrupt changes in depth or energy voids; pores or inclusions can appear as point-like or patchy anomalous regions; interlayer separation usually presents as a continuous shallow anomalous band after phase alignment.

[0231] In practical use, contour lines or isoenergy lines can be overlaid on the defect distribution map to assist in interpretation, and size, connectivity and morphological thresholds can be set to output the region labels and spatial coordinates of candidate defects.

[0232] In this way, the mapping from time of flight to spatial depth avoids false depths caused by phase reversal, while maintaining consistency between the time and geometric domains under the conditions of curved surfaces and gradual changes in refractive index. This makes the final defect distribution map more consistent with the actual structure in terms of location and shape, providing a reliable basis for subsequent quantitative assessment and engineering treatment.

[0233] Based on the same inventive concept, this application also provides a terahertz signal aliasing cancellation system based on deconvolution, which corresponds to a terahertz signal aliasing cancellation method based on deconvolution. Since the principle of the system in this application is similar to the aforementioned terahertz signal aliasing cancellation method based on deconvolution in this application, the implementation of the system can refer to the implementation of the method, and the repeated parts will not be described again.

[0234] Reference Figure 4 The diagram shown is a schematic of a terahertz signal aliasing cancellation system based on deconvolution according to an embodiment of this application. The system includes:

[0235] The acquisition module 10 is used to acquire time-domain waveforms of the reference area and the detection area to obtain the first reference waveform and the first observed waveform; and to acquire geometric and material parameters to obtain the first geometric prior parameters.

[0236] Processing module 20 is configured to calculate a point spread function based on the first reference waveform and the first geometric prior parameters to obtain a first point spread function; calculate a time axis normalized mapping function based on the first geometric prior parameters to obtain a first time normalized mapping function; and apply the first time normalized mapping function to the first observed waveform to obtain a first normalized observed waveform.

[0237] The restoration module 30 is used to perform deconvolution processing based on the first normalized observation waveform to obtain the restored pulse train;

[0238] Output module 40 is used to generate a time-of-flight map based on the restored pulse train and perform curvature correction based on the first geometric prior parameters 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 this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application.

Claims

1. A method for terahertz signal aliasing cancellation based on deconvolution, characterized in that, include: The time-domain waveforms of the reference region and the detection region are acquired to obtain the first reference waveform and the first observed waveform; Obtain the geometric and material parameters to obtain the first geometric prior parameters; The first point spread function is obtained by calculating the point spread function based on the first reference waveform and the first geometric prior parameters; The time axis normalized mapping function is calculated based on the first geometric prior parameters to obtain the first time normalized mapping function; Applying the first time normalization mapping function to the first observed waveform yields the first normalized observed waveform; Deconvolution processing is performed on the first normalized observed waveform to obtain the restored pulse train; A time-of-flight map is generated based on the restored pulse train, and curvature correction is performed based on the first geometric prior parameters to obtain a defect distribution map.

2. The terahertz signal aliasing cancellation method based on deconvolution according to claim 1, characterized in that, Also includes: The first observed waveform is deconvolved on the unnormalized original time axis to obtain the second restored pulse train; Applying the first time normalization mapping function to the second restored pulse train yields the second normalized restored pulse train; Based on the commutability difference value between the restored pulse train and the second normalized restored pulse train, the first coupling constraint value is obtained; Based on the first coupling constraint value, the convolution residual between the first normalized observation waveform and the first point spread function, and the sparsity index of the restored pulse train, a joint iterative update is performed to obtain the updated first time normalized mapping function, the updated first point spread function, and the updated restored pulse train.

3. The terahertz signal aliasing cancellation method based on deconvolution according to claim 2, characterized in that, Also includes: Based on the updated restored pulse train, the first stability index and the second stability index are calculated to obtain the first stopping condition; The joint iterative update is terminated when the first stopping condition is met, and the target recovery pulse train is output as the recovery pulse train.

4. The terahertz signal aliasing cancellation method based on deconvolution according to claim 1, characterized in that, The deconvolution process based on the first normalized observed waveform includes: Collect the normal information of the detection surface, determine the direction set based on the first geometric prior parameter and the normal information, and obtain the first direction set; Based on the first reference waveform and the first direction set, calculate the direction-related point diffusion function set to obtain the first direction kernel cluster; Based on the first normalized observation waveform and the first direction set, directional filtering is performed to obtain the first direction up-dimensional data; Based on the first direction up-dimensional data and the first direction kernel cluster, perform sparse deconvolution processing to obtain the restored pulse tensor of the direction index; Perform directional neighborhood smoothing on the restored pulse tensor of the directional index to obtain a directional continuous restored pulse tensor; The continuous recovery pulse tensor of the direction is projected onto the time axis and the scan axis according to the first direction set to obtain the recovery pulse train.

5. The terahertz signal aliasing cancellation method based on deconvolution according to claim 4, characterized in that, The determination of the direction set based on the first geometric prior parameters and the normal information includes: An initial direction set is generated based on the first geometric prior parameters and the normal information to obtain the first initial direction set; The first direction indication is obtained by calculating the direction indication within a preset time window and a preset scanning window based on the first normalized observation waveform. The first direction consistency score is obtained by calculating the direction consistency score based on the first direction indication and the first initial direction set; Within the angular spacing range defined by the first geometric prior parameter, the first initial direction set is shrunk and resampled based on the first direction consistency score to obtain the second direction set.

6. The terahertz signal aliasing cancellation method based on deconvolution according to claim 5, characterized in that, Also includes: Based on the first initial direction set and the second direction set, the change in the set is calculated to obtain the stability condition of the first direction; 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 used as a new initial direction set and the direction set generation steps are repeated until the first direction stability condition is met.

7. The terahertz signal aliasing cancellation method based on deconvolution according to claim 4, characterized in that, The step of performing sparse deconvolution processing on the first direction-upgraded data and the first direction kernel cluster to obtain the restored pulse tensor of the direction index includes: Based on the first geometric prior parameters and the normal information, the directional adjacency relationship is determined, and a first directional compatibility matrix is ​​generated; In the group sparse solution, a shared constraint is applied to the non-zero support of the direction dimension based on the first direction compatibility matrix, so that the non-zero direction corresponding to the sampling at the same time is limited to the adjacent direction specified by the first direction compatibility matrix, and the restored pulse tensor of the direction index is output.

8. The terahertz signal aliasing cancellation method based on deconvolution according to claim 7, characterized in that, The step of determining the directional adjacency relationship based on the first geometric prior parameters and the normal information, and generating the first directional compatibility matrix, includes: Based on the first geometric prior parameters and the normal information, the equivalent incident angle corresponding to each direction is determined, and the phase characteristics of the reflection coefficient corresponding to the equivalent incident angle are evaluated within the terahertz working frequency band. Based on the in-band phase sign consistency, mutually exclusive direction pairs are marked to obtain the first incompatible set. After excluding the first incompatible set, the group delay order of the main echoes at adjacent scanning positions is calculated based on the first geometric prior parameters, and the direction pairs are filtered based on the criterion of maintaining the group delay order to obtain the first compatible set. The direction pairs belonging to the first compatible set are set to allowed adjacency in the first direction compatibility matrix, and the direction pairs belonging to the first incompatible set are set to prohibited adjacency in the first direction compatibility matrix.

9. A terahertz signal aliasing cancellation method based on deconvolution according to claim 1, characterized in that, The obtained defect distribution map includes: The equivalent incident angle of each scanning position is determined based on the first geometric prior parameters, and the representative frequency of each restored pulse is determined based on the spectral energy centroid of the restored pulse train, thus obtaining the first representative frequency set. Based on the equivalent incident angle and the first representative frequency set, the interface reflection phase label is calculated to obtain the first phase label set; based on the analytical signal phase of the restored pulse train at the pulse peak, the phase label is calculated and measured to obtain the second phase label set. Phase consistency verification is performed based on the first phase tag set and the second phase tag set. When there is inconsistency, polarity correction and micro-delay correction are performed on the corresponding recovery pulse to obtain a phase-aligned recovery pulse train. A time-of-flight map is generated based on the phase-aligned restored pulse train, and an angle-resolved group refractive index is determined based on the equivalent incident angle and the first representative frequency set. Curvature correction is then performed to obtain the defect distribution map.

10. A terahertz signal aliasing cancellation system based on deconvolution, characterized in that, include: The acquisition module is used to acquire time-domain waveforms in the reference area and the detection area to obtain the first reference waveform and the first observed waveform. Obtain the geometric and material parameters to obtain the first geometric prior parameters; The processing module is used to calculate the point spread function based on the first reference waveform and the first geometric prior parameters to obtain the first point spread function; The time axis normalization mapping function is calculated based on the first geometric prior parameters to obtain the first time normalization mapping function; the first time normalization mapping function is applied to the first observed waveform to obtain the first normalized observed waveform; The restoration module is used to perform deconvolution processing based on the first normalized observed waveform to obtain the restored pulse train; The output module is used to generate a time-of-flight map based on the restored pulse train and perform curvature correction based on the first geometric prior parameters to obtain a defect distribution map.

Images