For a square integrable function f ( t ) ∈ L 2 ( R ) , ψ ( t ) is the continuous wavelet generator function, ψ a , b ( t ) is the continuous wavelet basis function, and the continuous wavelet transform of f ( t ) is defined as H f ( a , b ) = 〈 f ( t ) , ψ a , b ( t ) 〉 = a − 1 / 2 ∫ − ∞ + ∞ f ( t ) ψ ( t − b a ¯ ) d t   ,   a > 0 (3) Where H f ( a , b ) is the result of f ( t ) continuous wavelet transform, a is the scale factor, b is the time shift factor, and ψ ( ⋅ ) ¯ is the conjugate of the wavelet generating function. Taking a = a 0 j , b = a 0 j k ( j is the discretization index, k is the discretization coefficient, j , k ∈ Z ), the two are discretized, and the discrete wavelet transform can be defined as H f ( j , k ) = 〈 f ( t ) ,   ψ j , k ( t ) 〉 = a 0 − j 2 ∫ − ∞ + ∞ f ( t ) ψ ( a 0 j − k ) d t ¯ (4) Where H f ( j , k ) is the result of discrete wavelet transform of f ( t ) [15, 16].

In the detection of the sample under test, the detected signal usually contains certain noise, which can be considered as Gaussian white noise. The signal containing noise can be simply expressed as follows y j = f j + δ ⋅ z j   ,   j = 1 , ... , n (5) Where y j is the signal containing noise, f j is the ideal noise-free signal, z j is the doped noise, δ is the noise level, and n is the sampling length. The essence of signal denoising is to use the different characteristics of signal and noise to extract noise δ ⋅ z j from noise containing signal y j and restore the ideal signal f j . For the ideal effective signal, it is continuous in the time domain. Therefore, after the discrete wavelet transform, the energy of the effective signal is mainly concentrated in the low frequency subband, and the modulus of the generated wavelet coefficient would be relatively large. For Gaussian white noise, there is no continuity in the time domain, showing a strong randomness. Meanwhile, after conducting discrete wavelet transform, it also has a strong randomness, thus the energy of noise in wavelet domain is mainly concentrated in the high frequency subband, often still think of the Gaussian white noise, produced by the modulus of wavelet coefficient will be relatively small. Based on the above characteristics, the wavelet coefficients corresponding to noise are still satisfy the Gaussian white noise distribution. According to the definition of standard deviation and variance of random signals, the standard deviation of signals reflects the dispersion degree from each point in the discrete signal to the mean value of signals. Therefore, this paper proposes to take the standard deviation of wavelet coefficients of noise containing signals under the decomposition of each scale of discrete wavelet transform as the threshold to reflect the distribution range of noise wavelet coefficients to a certain extent.

Define the scale function φ ( t ) ∈ L 2 ( R ) by an integer multiple translation k to get φ 0 , k ( t ) , then φ 0 , k ( t ) = φ ( t − k )   ,   k ∈ Z (6) If space V 0 = S p a n k { φ 0 , k ( t ) } ¯ is defined by a signal linearly expressed by signal φ 0 , k ( t ) , then V 0 is called a closed signal space of φ 0 , K ( t ) .By translating and stretching the scale function, a new scale function φ j , k ( t ) is obtained, which is expressed as φ j , k ( t ) = 2 j / 2 φ ( 2 j t − k ) ,   j ∈ Z   ,   k ∈ Z (7) Similarly, the newly closed signal space spanned by φ j , k ( t ) is φ j , k ( t ) = V j = S p a n k { φ j , k ( t ) } ¯ = S p a n k { 2 j / 2 φ ( 2 j t − K } ¯   ,   j ∈ Z   ,   k ∈ Z (8) After further derivation of the scale function φ j , k ( t ) , the multi-resolution analysis equation corresponding to the scale function coefficient is obtained, and the expression is φ ( t ) = ∑ k h 0 [ k ] ⋅ 2 φ ( 2 t − k ) (9) The equation shows that different scale functions correspond to different h 0 [ k ] [17, 18]. When the noise containing signal y i ∈ V j + 1 , the noise containing signal can be expressed as y i = ∑ k c j [ k ] ⋅ φ j , k ( t ) + ∑ k d j [ k ] ⋅ ψ j , k ( t ) (10)

Among them, there is c j [ k ] = 〈 y i , φ j , k ( t ) 〉 ,   d j [ k ] = 〈 y i , ψ j , k ( t ) 〉 (11)

According to the multi-resolution analysis equation φ t = ∑ k h 0 [ n ] ⋅ 2 φ ( 2 t − n ) of the scale function, we can obtain φ j , k ( t ) = ∑ n h 0 [ n − 2 k ] ⋅ φ j + 1 , n ( t ) (12)

Then there is c j [ k ] = 〈 y i , φ j , k ( t ) 〉 = ∑ n h 0 [ n − 2 k ] ⋅ 〈 y i , φ j + 1 , n ( t ) 〉 = ∑ n h 0 [ n − 2 k ] ⋅ c j + 1 [ n ] (13)

Similarly, according to ψ ( t ) = ∑ n h 1 [ n ] ⋅ 2 φ ( 2 t − n ) , the wavelet coefficient d j [ k ] at each scale and each decomposition layer can be obtained as d j [ k ] = ∑ n h 1 [ n − 2 k ] ⋅ c j + 1 [ n ] (14) Set the standard deviation of the wavelet coefficient d j [ k ] as threshold U , and the expression is U = 1 M − 1 ∑ q = 1 M { ( ∑ n h 1 [ n − 2 k ] ⋅ c j + 1 [ n ] ) q − ∑ q = 1 M ( ∑ n h 1 [ n − 2 k ] ⋅ c j + 1 [ n ] ) q M } 2 (15) where M is the number of wavelet coefficients of the noise containing signal under the decomposition of each scale of the discrete wavelet transform.

According to the characteristics of Gaussian distribution, setting the wavelet coefficient within the threshold range to zero can suppress the interference of noise to the maximum extent, namely, the effect of hard threshold function. Therefore, this paper selects the hard threshold function to carry out the hard threshold shrinkage denoising process for the wavelet coefficient of the noise containing signal. When the absolute value of the wavelet coefficient of the noise containing signal is less than the set threshold, it is made to be zero; otherwise, it keeps unchanged. After shrinkage denoising, the wavelet coefficients are reconstructed and restored to the terahertz detection signal, and the amplitude imaging method is used for imaging processing to obtain the preliminary terahertz detection image.

Guided filtering [9–11] is a filtering method with good performance in image denoising, defogging and detail enhancement. Its core idea is to guide the input image through a guidance image based on the local linear model, so that the overall contour feature of the output image is similar to that of the input image, and the texture detail region is similar to that of the guidance image. The output image q can be expressed as follows q = f g u i d e r ( p , I , r , ε ) (16) Where p is the input image, i.e., the image after wavelet denoising, I is the guidance image, r is the size of the filtering window, ε is the regularization parameter, r and ε can be determined by empirical value. In this paper, r is 3 and ε is 0.64. The assumption between output image q and guidance image I is as follows : there is a linear relationship in local window w v centered on pixel v , and the output expression of a pixel can be expressed as q i = g v I i + o v   ,   ∀ i ∈ w v (17) Where i corresponds to the pixel index in the local window, g v and o v are linear coefficients in the window. By calculating the derivative gradient on both sides of the above equation, ∇ q = g ∇ I can be obtained. It can be seen that when the input guidance image has gradient in a certain region, the output image will also maintain the corresponding gradient. Therefore, the guided filter has good edge preservation performance while smoothing the background. In order to obtain the optimal solution of g v and o v , it is necessary to retain the effective information of the input image p as much as possible in the output image q . Even if the difference between the two is the smallest, the method is usually to introduce the minimization cost function and minimize it to solve the optimization problem. E ( g v , o v ) is the minimization cost function, and the expression is as follows E ( g v , o v ) = ∑ i ∈ w v [ ( g v I i + o v − p i ) 2 + ε g v 2 ] (18) Where ε g v 2 is the penalty term, ε is used to prevent g v too large regularization parameters. The Eq. 18 is solved by linear regression model to minimize it, and the optimal solution of g v and o v can be obtained as g v = 1 | w | ∑ i ∈ w v I i p i − μ v p v ¯ σ v 2 + ε (19) o v = p v ¯ − g v μ v (20) Where μ v and σ v 2 are the mean value and variance of pixels in local window w v of the guidance image, | w | is the number of pixels contained in w v , and p v ¯ = 1 | w | ∑ i ∈ w v p i is the mean value of image p to be filtered in the window. When using Eq. 19 to calculate the linear coefficient, the g v and o v calculated in different window w v are obviously different, which can be solved by averaging the values obtained from g v and o v , so the output image can be expressed as q i = 1 | w | ∑ i ∈ w k ( g v I i + o v ) = g i ¯ I i + o i ¯ (21) Where g i ¯ = 1 | w | ∑ v ∈ w v g v and o i ¯ = 1 | w | ∑ v ∈ w v o v are the average of linear coefficients of all local windows under the same pixel index.

Unsharp masking [19, 20] is a common sharpening enhancement technique. Its basic principle is that the high frequency part of the image is obtained by subtracting the original image from the blurred image that is obtained by low-pass filtering, and then the high frequency part of the image is gained and superimposed with the original image to obtain the resulting image with enhanced details and edges. The unsharp masking process can be expressed as follows l ( i , j ) = f ( i , j ) + ζ × { f ( i , j ) − L P [ f ( i , j ) ] } (22) where f ( i , j ) is the original input image, L P is the low-pass filter, ζ is the gain coefficient of high-frequency details, and l ( i , j ) is the enhanced output image.

In this paper, the mean filter is selected as the low-pass filter smoothing filter, and the mean filter template with a size of 3 × 3 is generated to conduct convolution filtering on the image after wavelet denoising, so as to obtain its blurred image. According to the principle of unsharp masking, the difference operation between the original image and the blurred image is used to obtain the mask image, that is, the image reflecting the detail information of the image. After ζ factor multiplication of the mask image, it is superimposed with the original image to obtain the sharpened image with enhanced detail and edge information in the image after wavelet denoising, so as to achieve the effect of image sharpening. Using this linear unsharp masking method to process the terahertz detection image can get better sharpening effect. While the sharpness of the details is improved, the smooth region in the original image is not affected and preserved. Finally, the sharpened image of unsharp masking is used as the guidance image to conduct the guided filtering on the image after wavelet denoising, so that the texture details in the image after wavelet denoising are similar to those in the guidance image.

Image P m after wavelet denoising is processed by guided filtering technology with its own image P m and sharpened image P s as the guidance image, the detail layer P w reflecting the image detail information and the base layer P e reflecting the image contour are obtained. The obtained detail layer P w and the base layer P e are superposed with the gain of β and γ coefficients respectively, and the terahertz detection resulting image P o u t with improved contrast, enhanced detail and edge information is obtained. The process can be expressed as follows P o u t = β × P w + γ × P e (23) Where P o u t is the output terahertz detection resulting image, β and γ are the gain coefficients, P w is the detail layer, and P e is the base layer.

Polymethacrylimide foam composites are heat-resistant composites with the highest strength and stiffness to weight ratio. As the core material of sandwich structure, polymethacrylimide foam composites are widely used in aerospace, shipbuilding, military and other fields [21]. As one of the three major polymer materials, rubber has high elasticity than metal, high mechanical strength and good bending resistance. Because of its unique properties, rubber is widely used in medical and health, power communication, civil engineering and other fields [22]. Glass fiber reinforced composites are widely used in the bearing and electromagnetic wave transmission structures of unmanned aerial vehicles due to their excellent strength-weight ratio, aerodynamic performance and microwave permeability [23]. With the increasingly wide application of the above composites, the demand for its quality control, nondestructive testing and evaluation is also increasing. In the process of preparation and application, the above composites will be affected by external forces and temperature variations, resulting in internal defects such as cracks, voids and adhesive debonding, affecting the normal use of composite materials, and when the damage reaches a certain level, it will even cause a catastrophe [24], so the accurate detection of composite materials has important practical significance.

After the imaging detection of the above composite materials, the amplitude imaging method was used to image the obtained terahertz detection signal, and Gaussian filtering and the method proposed in this paper were used to denoise and enhance the terahertz detection image. When the hard threshold shrinkage denoising method based on discrete wavelet transform was used to denoise the original detection signal, for the adhesive debonding defect with a diameter of 10.0 mm in sample 1, “sym4” was selected as the wavelet denoising combination with 3 layers of decomposition layer, and “sym5” was selected as the wavelet denoising combination with four layers of decomposition layer for the adhesive debonding defect with a diameter of 6.0 mm; For the delamination defects in sample 2, “sym5” was selected as the wavelet denoising combination with two layers of decomposition layer; For the adhesive debonding defect with a diameter of 10.0 mm and the adhesive debonding defect with a diameter of 3.0 mm in sample 3, the combination of wavelet denoising with “sym5” as the wavelet basis and the decomposition layer of 3 layers was selected for hard threshold shrinkage denoising.

Standard deviation refers to the degree of dispersion of image pixel gray value relative to the mean value. If the standard deviation is larger, it indicates that the gray scale in the image is more dispersed and the image quality is better. Its calculation formula is S D = 1 G R ∑ i = 1 G ∑ j = 1 R ( F ( i , j ) − u ) 2 (24) where G and R represent the sizes of the image, F ( i , j ) represents the gray value of pixel F ( i , j ) corresponding to image F , and u represents the average gray value of all pixels in the image.

The mean gradient represents the ability of image to express details, which is a measure of image sharpness. Usually, the greater the mean gradient is, the more the details of the image are, and the higher the sharpness is. The calculation formula of the mean gradient is as follows M G = 1 G R ∑ i = 1 G ∑ j = 1 R [ ∂ f ( x , y ) ∂ x ] 2 + [ ∂ f ( x , y ) ∂ y ] 2 (25) where ∂ f ( x , y ) ∂ x represents the gradient in the horizontal direction, and ∂ f ( x , y ) ∂ y represents the gradient in the vertical direction.

Information entropy is an effective method to evaluate information content in images. For an image with gray scale in the range of [0,L-1], the expression of information entropy is I E = − ∑ i = 0 L − 1 p ( s i ) log 2 ⁡ p ( s i ) (26) Where p ( s i ) is the probability of gray level s i , the greater the information entropy indicates that the image contains more details.

Image sharpness is an important index to measure the quality of the image, which can better correspond to the subjective feelings of people. The energy gradient can evaluate the image sharpness in real time, and the greater the energy gradient value of the image, the better the image sharpness. The definition of the function is as follows E G = ∑ x = 1 G ∑ y = 1 R ( | f ( x + 1 , y ) − f ( x , y ) | 2 + | f ( x , y + 1 ) − f ( x , y ) | 2 ) (27) Where f ( x , y ) represents the gray value of image f at the corresponding pixel ( x , y ) .

The local contrast of the image can reflect the dynamic range of the gray distribution of the image pixels. The larger the dynamic range is, the higher the contrast is. The formula for calculating the local contrast is L C = 1 G R ∑ x = 1 G ∑ y = 1 R I x , y max − I x , y min I x , y max + I x , y min + C (28) Where I x , y max and I x , y min represent the maximum and minimum pixel gray values in a sub-block centered on pixel ( x , y ) , respectively. In this paper, the size of the sub-block is 5 × 5. In order to avoid the case of denominator 0, C is a smaller constant with a value of 0.1.

The calculation results of objective evaluation indexes of different methods are shown in Table 1. It can be seen from the table that the standard deviation, mean gradient, information entropy, energy gradient and local contrast of the five objective evaluation index data are higher than the corresponding values of the original detection image and the detection image after Gaussian filtering. It shows that the method in this paper has obvious advantages in contrast enhancement and level of detail and sharpness improvement.