Terahertz time-domain spectral imaging method based on femtosecond laser pumping-detection principle

Abstract

The invention discloses a terahertz time-domain spectral imaging method based on a femtosecond laser pumping-detection principle. The terahertz time-domain spectral imaging method comprises the following steps: acquiring a terahertz wave beam signal of a to-be-detected sample through a terahertz time-domain spectral system; gaussian fitting is carried out on the terahertz wave beam signal by using the terahertz wave beam model; and performing imaging reconstruction on the Gaussian fitting result of the terahertz wave beam signal by adopting a maximum likelihood and maximum expected value strategy to obtain an imaging reconstruction result. According to the terahertz time-domain spectral imaging method based on the femtosecond laser pumping-detection principle, the reconstruction technology combining the Gaussian beam model and the MLEM algorithm is introduced, the wave path calculation error is effectively reduced, the scanning imaging quality is improved, and high-resolution two-dimensional terahertz imaging is achieved; according to the method, the structural characteristics of fiber distribution, agglomeration, interruption and the like in the cigarette can be clearly presented, the defects of image blurring and multiple artifacts of a traditional back projection algorithm are overcome, and an efficient and reliable solution is provided for nondestructive testing and quality evaluation of the internal structure of a cigarette product.

Description

Terahertz time-domain spectral imaging method based on femtosecond laser pump-probe principle Technical Field

[0001] This invention relates to the field of terahertz time-domain spectroscopy, and more specifically, to a terahertz time-domain spectroscopy imaging method based on the femtosecond laser pump-probe principle. Background Technology

[0002] Terahertz waves (THz) generally refer to electromagnetic waves with frequencies ranging from 0.1 to 10 THz, located between the microwave and infrared bands. Terahertz waves combine the excellent penetrability of microwaves with the high resolution of light waves, exhibiting characteristics such as strong penetration, good security, and high spectral resolution. They have broad application potential in non-destructive testing of materials, identification of pharmaceutical components, security inspection, biomedical imaging, and communications. Since the vibrational and rotational energy levels of many materials and molecules fall within the terahertz band, terahertz time-domain spectroscopy (THz-TDS), as an emerging spectroscopic analysis and imaging technique, can effectively extract the optical parameters of samples, enabling qualitative and quantitative analysis.

[0003] THz-TDS systems are typically built on the femtosecond laser pump-probe principle: after the femtosecond laser pulse is split, one pump beam excites a photoconductive antenna or nonlinear crystal to generate a broadband terahertz pulse, while the other probe beam overlaps with the terahertz signal in an electro-optic crystal. Through the electro-optic effect, the terahertz electric field information is converted into a change in the probe beam's polarization. Finally, a balanced detector quantizes the signal and reconstructs the time-domain waveform. After Fourier transform, frequency-domain information can be obtained, allowing the extraction of optical parameters such as the sample's absorption coefficient, refractive index, and dispersion characteristics. Compared to Fourier transform infrared spectroscopy (FTIR), THz-TDS not only offers higher temporal resolution and non-destructive testing capabilities but also directly acquires amplitude and phase information, thus providing greater advantages in quantitative analysis and imaging.

[0004] However, existing THz-TDS systems still have some key problems in application: 1) System errors: Most parameter calculations assume that THz waves are plane waves, while the actual propagation process is more in line with the Gaussian beam model. Although the difference between the two calculation results is small, it is still not negligible in high-precision imaging and parameter measurement; 2) Random errors: mainly from sample thickness measurement errors and incident angle errors. In experiments, the thickness measurement error can reach 1%, and the incident angle deviation will further amplify the error in the path calculation, thus affecting the accuracy of optical parameter extraction; 3) Beam diffusion effect: In scanning imaging, the Gaussian distribution of terahertz waves leads to a large point spread function, which in turn causes a decrease in resolution and image blurring; 4) Limitations of data acquisition and imaging algorithms: Traditional filtered back projection (FBP) algorithms are sensitive to noise and it is difficult to obtain high-quality reconstruction results; conventional inversion methods have low recognition rates when baseline drift, astigmatism and environmental noise are present, generally around 90%, which is difficult to meet the high-precision requirements of practical applications; 5) Insufficient real-time performance: THz scanning imaging process is often time-consuming, which is not conducive to rapid detection and online monitoring scenarios.

[0005] Therefore, there is an urgent need for a terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle. Summary of the Invention

[0006] The purpose of this invention is to provide a terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle to solve the problems in the prior art and achieve high-resolution two-dimensional terahertz imaging.

[0007] This invention provides a terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle, comprising:

[0008] The terahertz beam signal of the sample under test is obtained by a terahertz time-domain spectroscopy system.

[0009] The terahertz beam signal is fitted with Gaussian using a pre-constructed terahertz beam model;

[0010] The Gaussian fitting results of the terahertz beam signal are used to perform imaging reconstruction by employing the maximum likelihood maximum expectation strategy, and the imaging reconstruction results are obtained.

[0011] The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle described above, preferably, involves acquiring the terahertz beam signal of the sample under test using a terahertz time-domain spectroscopy system, including:

[0012] The output pulse of the femtosecond laser source is split into pump light and probe light by a beam splitter;

[0013] The pump light emitted by the terahertz transmitting unit is focused and then irradiates the low-temperature grown gallium arsenide photoconductive antenna to generate broadband terahertz pulses.

[0014] The terahertz pulse generated by the terahertz transmission optical path is collimated by a lens and transmitted to the sample under test to form a terahertz light spot signal.

[0015] The terahertz detection unit, wherein the detection light and the transmitted terahertz wave are collinearly incident on the ZnTe electro-optic crystal, and the terahertz electric field intensity is converted into a polarization state change through the Pockels effect, and the signal is obtained by the Wollaston prism beam splitting and the balanced detector.

[0016] Full waveform sampling is achieved by adjusting the optical path of the detection optical path through the delay control module;

[0017] The time-domain signal is subjected to Fourier transform by the data processing module to obtain frequency-domain information and thus the terahertz beam signal of the sample under test.

[0018] The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle described above, preferably, involves Gaussian fitting of the terahertz beam signal using a pre-constructed terahertz beam model, comprising:

[0019] During the scanning imaging process, the received signal is represented by the convolution of the transmitted waveform and the sample function;

[0020] The beam waist radius represents the beam diameter;

[0021] Energy distribution is represented by Rayleigh distance.

[0022] The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle described above preferably includes, during the scanning imaging process, representing the received signal by the convolution of the emitted waveform and the sample function, comprising:

[0023] The received signal during the scanning process can be represented by the following formula:

[0024]

[0025] in, This indicates the result of the scanning imaging. Let n represent the function of the sample to be tested, and n represent the additive noise. The beam model represents the point spread function, and * represents the convolution operation;

[0026] In terahertz transmission imaging, the beam passes through the sample under test along the z-axis. Since the sample is three-dimensional, the precise three-dimensional model of the imaging process can be represented by the following formula:

[0027]

[0028] in, Indicates the location of the terahertz transmitter. The position of the terahertz detection end is indicated, and the accurate model of the PSF is expressed by the following formula. :

[0029]

[0030] Where NA represents the numerical aperture and k represents the cutoff factor. This represents the distance from the sample to the z-axis. Indicates the attenuation coefficient. The frequency at the focal plane is . The spot function,

[0031] The representation of the beam waist radius to indicate the beam spot diameter includes:

[0032] The waist radius is expressed by the following formula: The distance from the focal point is Gaussian beam spot diameter :

[0033]

[0034] in, Indicates Rayleigh distance;

[0035] The representation of energy distribution using Rayleigh distance includes:

[0036] The Gaussian beam energy distribution is expressed by the following formula. :

[0037]

[0038] in, This represents the distance from the Gaussian beam to the optical axis. This represents the energy at the center point of the beam.

[0039] In the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle described above, preferably, the beam model of the point spread function... This describes the relationship between numerical aperture, attenuation coefficient, spot function, and frequency, wherein the beam model of the point spread function... The construction process includes:

[0040] In a terahertz two-dimensional scanning imaging system, the imaging response can be considered as a linear space-invariant model, and the point spread function (PSF) can be approximated by a Gaussian beam, with the expression:

[0041]

[0042] in: Indicates terahertz frequency, The frequency-dependent beam radius is used to reflect the scattering and dispersion characteristics of the beam in space.

[0043] The minimum spot radius of the Gaussian focal spot formed by a terahertz beam through a focusing lens or reflector is expressed by the following formula. The physical relationship between the numerical aperture (NA) of the system and the system's numerical aperture (NA):

[0044]

[0045] in: This indicates the terahertz wavelength, which decreases as the frequency increases. The larger the NA, the smaller the light spot, the sharper the PSF, and the higher the resolution.

[0046] The minimum spot radius of the Gaussian focal spot formed by a terahertz beam passing through a focusing lens or reflector. The physical relationship between the PSF and the numerical aperture (NA) of the system is expressed by the following formula, which represents the frequency dependence of the PSF:

[0047] ;

[0048] In terahertz systems, propagation distance This leads to spot expansion, and the propagation distance dependence of PSF can be expressed by the following formula:

[0049] ;

[0050] Rayleigh length is expressed by the following formula:

[0051] ;

[0052] The attenuation coefficient of terahertz waves in a medium due to absorption loss is expressed by the following formula:

[0053]

[0054] in, This represents the frequency-dependent absorption coefficient; the higher the frequency, the stronger the attenuation.

[0055] The propagated PSF expression, expressed in its complete form with an attenuation factor, is given by the following formula:

[0056]

[0057] in, This indicates that PSF is determined by frequency. Numerical aperture (NA), spot function and attenuation coefficient A joint decision.

[0058] The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle described above preferably includes the following: The imaging reconstruction result obtained by using a maximum likelihood maximum expectation strategy to perform Gaussian fitting of the terahertz beam signal includes:

[0059] Construct a projection model to express the imaging problem as a matrix equation;

[0060] MLEM iterative reconstruction based on projection model to iteratively solve the image;

[0061] Deconvolution operations are introduced during the iteration process to perform super-resolution optimization;

[0062] The reconstructed terahertz image is output, resulting in a high-resolution two-dimensional image.

[0063] The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle described above preferably includes the following: Constructing a projection model to express the imaging problem as a matrix equation includes:

[0064] The scan data can be represented by a matrix equation as follows:

[0065]

[0066] middle, Represents the column vector of the original image. Represents the weight matrix, This represents the collected data, where n represents the noise level. Let y represent the projection weight matrix of size M×N, where M represents the total number of rays, N represents the total number of image pixels, each element value represents the length of the ray passing through that pixel, and y represents the data collected, which is of size M×1.

[0067] The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle described above, preferably, involves iterative image reconstruction using MLEM based on a projection model, which includes:

[0068] In the imaging model, the response of the terahertz system is expressed by the following formula:

[0069]

[0070] in, This represents the measured signal at the scan point. This represents the spatial reflectance distribution to be restored. To represent the system matrix containing the terahertz beam spread function;

[0071] Based on the EM derivation using Poisson likelihood, the update formula for terahertz MLEM is obtained as follows:

[0072]

[0073] Among them, forward projection It is used to calculate the predicted measurement, and the back projection is used to realize the error backpropagation, and through the normalization term Ensure the stability and non-negativity of updates;

[0074] MLEM iteration is performed based on the update formula of terahertz MLEM.

[0075] The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle described above, preferably, includes the following step: MLEM iteration based on the terahertz MLEM update formula.

[0076] Expected steps: Estimate the current image Next, calculate the expected value of the observed data with respect to the complete likelihood function, and obtain the result for each measurement. Belongs to each pixel The probability contribution;

[0077] Maximization step: Using the statistical weights obtained from the expectation step, maximize the expected likelihood function to update the image estimate.

[0078] The desired steps include:

[0079] The predicted measurement is calculated using the following formula:

[0080]

[0081] The attribution coefficient is calculated using the following formula:

[0082]

[0083] The attribution coefficient represents the ratio of measurement to prediction. This ratio reflects the deviation between the current model's prediction and the actual measurement. It is a weighting factor for back projection and is a simplified result of the following expectation:

[0084] ,

[0085] The maximization step includes:

[0086] The update formula for MLEM is expressed by the following formula:

[0087]

[0088] in, This indicates that the measurement error is back-projected back to pixel j according to the system response matrix; the denominator is... It is a normalization term to ensure that the updated value is stable and remains non-negative.

[0089] The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle described above, preferably, includes the introduction of deconvolution operations during the iteration process for super-resolution optimization, comprising:

[0090] Convolutional modeling: When using the current image estimate to generate the forward projection, it is expressed by the following formula:

[0091]

[0092] This is represented as a convolution of the image with a terahertz PSF to simulate the spread blurring of a real beam;

[0093] Calculate the correction factor: Based on the ratio of the measured value to the convolutional predicted value, the following formula is obtained:

[0094] ;

[0095] Backprojection and deconvolution updates are performed based on the correction factor:

[0096]

[0097] The back-projection part is equivalent to deconvolving the correction factor back to the pixel position.

[0098] This invention provides a terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle. By introducing a reconstruction technique combining a Gaussian beam model and the MLEM algorithm, it effectively reduces path length calculation errors and improves scanning imaging quality, achieving high-resolution two-dimensional terahertz imaging. By introducing a Gaussian beam model into the imaging reconstruction process and combining iterative optimization with the maximum likelihood maximum expectation algorithm, it effectively compensates for beam dispersion effects in actual terahertz wave propagation. A projection weight matrix containing prior beam knowledge is constructed, and deconvolution operations are performed during the MLEM iteration process, significantly improving image resolution and signal-to-noise ratio. The method of this invention can clearly present the structural features of internal fibers, agglomeration, and discontinuity in cigarettes, overcoming the shortcomings of traditional back-projection algorithms such as image blurring and numerous artifacts. It provides an efficient and reliable solution for non-destructive testing and quality assessment of the internal structure of cigarette products. Attached Figure Description

[0099] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described below with reference to the accompanying drawings, wherein:

[0100] Figure 1 is a flowchart of an embodiment of the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle provided by the present invention;

[0101] Figure 2 is a schematic diagram of the terahertz beam model;

[0102] Figure 3 shows the measured terahertz light spot signal;

[0103] Figure 4 shows the Gaussian fitting results of Figure 3;

[0104] Figure 5 shows the direct scanning imaging results of fresh leaves;

[0105] Figure 6 shows the imaging results of fresh leaves using the maximum likelihood algorithm.

[0106] Figure 7 shows the terahertz amplitude and phase signal distribution of the finished cigarette through direct scanning imaging;

[0107] Figure 8 shows the amplitude and phase signal distribution of the finished cigarette sticks as captured by the maximum likelihood algorithm.

[0108] Figure 9 is a schematic diagram of the discontinuous imaging analysis results of the finished cigarette sticks. Detailed Implementation

[0109] Various exemplary embodiments of the present disclosure will now be described in detail with reference to the accompanying drawings. The descriptions of the exemplary embodiments are merely illustrative and are in no way intended to limit the present disclosure or its application or use. The present disclosure may be implemented in many different forms and is not limited to the embodiments described herein. These embodiments are provided so that the present disclosure will be thorough and complete, and will fully express the scope of the disclosure to those skilled in the art. It should be noted that, unless specifically stated otherwise, the relative arrangement of components and steps, the composition of materials, numerical expressions, and values ​​set forth in these embodiments should be interpreted as exemplary only and not as limiting.

[0110] The terms “first,” “second,” and similar terms used in this disclosure do not indicate any order, quantity, or importance, but are merely used to distinguish different parts. Terms such as “including” or “contains” mean that the element preceding the term encompasses the element listed after it, and do not exclude the possibility of encompassing other elements as well. Terms such as “above” and “below” are used only to indicate relative positional relationships; when the absolute position of the described object changes, this relative positional relationship may also change accordingly.

[0111] In this disclosure, when a specific component is described as being located between a first component and a second component, an intermediary component may or may not be present between the specific component and the first or second component. When a specific component is described as connecting to other components, the specific component may be directly connected to the other components without having an intermediary component, or it may not be directly connected to the other components but may have an intermediary component.

[0112] All terms used in this disclosure (including technical or scientific terms) have the same meaning as understood by one of ordinary skill in the art to which this disclosure pertains, unless otherwise specifically defined. It should also be understood that terms defined in a general dictionary, such as a dictionary, should be interpreted as having a meaning consistent with their meaning in the context of the relevant art, and not as having an idealized or highly formalized meaning, unless expressly defined herein.

[0113] Techniques, methods, and equipment known to those skilled in the art may not be discussed in detail, but where appropriate, they should be considered part of the specification.

[0114] To address the problems existing in current THz-TDS systems during application, such as systematic errors, random errors, beam blurring effects, limitations in data acquisition and imaging algorithms, and insufficient real-time performance, several improvement methods have been proposed by academia and engineering applications. These include reducing errors by optimizing thickness measurement methods, eliminating baseline drift through signal preprocessing, and improving recognition rates using statistical learning methods. However, these methods often rely on complex calculations, accurate initial values, or are sensitive to noise, making it difficult to balance system stability and imaging resolution.

[0115] Therefore, there is an urgent need for an imaging system and algorithm that can simultaneously address systematic error correction, random error suppression, and beam dispersion compensation. This invention introduces a Gaussian beam model into the existing THz-TDS system and combines it with the Maximum Likelihood Expectation Maximization (MLEM) algorithm for image reconstruction. This is expected to improve resolution while enhancing the system's robustness to noise and errors, providing a better solution for the application of terahertz waves in high-resolution imaging and material detection.

[0116] As shown in Figure 1, the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle provided in this embodiment includes the following steps in actual execution:

[0117] Step S1: Obtain the terahertz beam signal of the sample under test using a terahertz time-domain spectroscopy (THz-TDS) system.

[0118] Specifically, a terahertz time-domain spectroscopy experimental system based on the femtosecond laser pump-probe principle acquires the terahertz beam signal of the sample under test. In one embodiment of the terahertz time-domain spectroscopy imaging method based on the femtosecond laser pump-probe principle of the present invention, step S1 may specifically include:

[0119] Step S11: The output pulse of the femtosecond laser source is split into pump light and probe light by a beam splitter.

[0120] In this invention, the femtosecond laser source is a femtosecond laser with a center wavelength of 800 nm and a pulse width of less than 100 fs. The laser pulse emitted by the femtosecond laser source is split into two paths by a beam splitter: one path is used as pump light, and the other path is used as probe light.

[0121] Step S12: The pump light emitted by the terahertz transmitting unit is focused and then irradiates the low-temperature grown gallium arsenide (LT-GaAs) photoconductive antenna to generate a broadband terahertz pulse.

[0122] The generated broadband terahertz pulses have a spectral range of 0.1–5 THz.

[0123] Step S13: The generated terahertz pulse is collimated by a lens and transmitted to the sample under test through the terahertz transmission optical path to form a terahertz light spot signal.

[0124] The generated terahertz pulses are collimated by a high-density polyethylene lens before penetrating the sample under test for terahertz transmission and interaction with the sample. In practice, the sample is fixed on a three-dimensional scanning frame, which can be precisely adjusted in the X, Y, and Z directions and can rotate 360°.

[0125] Step S14: Terahertz detection unit, where the detection light and the transmitted terahertz wave are collinearly incident on the ZnTe electro-optic crystal, and the terahertz electric field intensity is converted into a polarization state change through the Pockels effect, and the signal is obtained by the Wollaston prism beam splitting and the balanced detector.

[0126] Among them, the balanced detector and the balanced photodetector, the polarization state change signal is split by the Wollaston prism and the balanced photodetector to obtain the orthogonal polarization intensity difference signal.

[0127] Step S15: Adjust the optical path of the detection optical path through the delay control module to achieve full waveform sampling.

[0128] Specifically, a linear displacement stage is installed along the probe optical path, with a scanning range of 0–300 ps and a stepping accuracy better than 0.1 fs. The full waveform of the terahertz electric field is acquired by controlling the delay.

[0129] Step S16: Perform Fourier transform on the time domain signal through the data processing module to obtain frequency domain information and obtain the terahertz beam signal of the sample under test.

[0130] Specifically, the time-domain electric field signal ETHz(t) is recorded and subjected to Fast Fourier Transform (FFT) to obtain the corresponding frequency domain information. The system operates at a temperature of 21°C and a humidity of 4%, with a spectral resolution of 12.5 GHz and an average of 1024 scans.

[0131] Step S2: Perform Gaussian fitting on the terahertz beam signal using a pre-constructed terahertz beam model.

[0132] The terahertz beam is focused by a lens in the system to form a spot that conforms to a Gaussian beam model. In one embodiment of the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle of the present invention, step S2 may specifically include:

[0133] Step S21: During the scanning imaging process, the received signal is represented by the convolution of the transmitted waveform and the sample function.

[0134] In one embodiment of the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle of the present invention, step S21 may specifically include:

[0135] Step S211: The received signal during the scanning process is expressed by the following formula:

[0136]

[0137] in, This indicates the result of the scanning imaging. Let n represent the function of the sample to be tested, and n represent the additive noise. The beam model represents the point spread function, and * represents the convolution operation.

[0138] In terahertz imaging systems, the THz beam is focused into a spot by a set of lenses. Studies have shown that the terahertz beam is a Gaussian beam, and the received signal during scanning can be considered as the convolution of the transmitted waveform and the target function. The beam model of the point spread function is described below. This describes the relationship between numerical aperture, attenuation coefficient, beam spread function, and frequency; specifically, it refers to the beam model of the point spread function. The construction process includes:

[0139] Step SA1: In a terahertz two-dimensional scanning imaging system, the system's imaging response can be considered as a linear spatially invariant (LSI) model. The point spread function (PSF) is approximated by a Gaussian beam, and its expression is:

[0140]

[0141] in: Indicates terahertz frequency, This represents the frequency-dependent beam radius, used to reflect the scattering and diffusion characteristics of the beam in space.

[0142] Step SA2: The minimum spot radius of the Gaussian focal spot formed by a terahertz beam through a focusing lens or reflector is expressed by the following formula. The physical relationship between the numerical aperture (NA) of the system and the system's numerical aperture (NA):

[0143]

[0144] in: This indicates the terahertz wavelength, which decreases as the frequency increases. The larger the NA, the smaller the light spot, the sharper the PSF, and the higher the resolution.

[0145] Step SA3: Determine the minimum spot radius of the Gaussian focal spot formed by the terahertz beam through the focusing lens or reflector. The physical relationship between the PSF and the numerical aperture (NA) of the system is expressed by the following formula, which represents the frequency dependence of the PSF:

[0146] .

[0147] Step SA4: In a terahertz system, the propagation distance This leads to spot expansion, and the propagation distance dependence of PSF can be expressed by the following formula:

[0148] .

[0149] Step SA5: Express the Rayleigh length using the following formula:

[0150] .

[0151] Step SA6: The attenuation coefficient of terahertz waves in the medium due to absorption loss is expressed by the following formula:

[0152]

[0153] in, This represents the absorption coefficient related to frequency; the higher the frequency, the stronger the attenuation.

[0154] Step SA7: Express the propagated PSF expression in its complete form with attenuation factor, as shown by the following formula:

[0155]

[0156] in, This indicates that PSF is determined by frequency. Numerical aperture (NA), spot function and attenuation coefficient Together, these factors determine that high frequency and high NA conditions can reduce the beam size and improve system resolution, while propagation distance and absorption loss lead to beam dispersion and reduce reconstruction accuracy.

[0157] Step S212: In terahertz transmission imaging, the beam passes through the sample under test along the z-axis. Since the sample under test is three-dimensional, the accurate three-dimensional model of the imaging process can be represented by the following formula:

[0158]

[0159] in, Indicates the location of the terahertz transmitter. The position of the terahertz detection end is indicated, and the accurate model of the PSF is expressed by the following formula. :

[0160]

[0161] Where NA represents the numerical aperture and k represents the cutoff factor. This represents the distance from the sample to the z-axis. Indicates the attenuation coefficient. The frequency at the focal plane is . The light spot function.

[0162] Step S22: The beam waist radius represents the beam diameter.

[0163] Specifically, the waist radius is expressed by the following formula: The distance from the focal point is Gaussian beam spot diameter :

[0164]

[0165] in, This represents the Rayleigh distance.

[0166] Step S23: Represent the energy distribution using Rayleigh distance.

[0167] Specifically, the Gaussian beam energy distribution is represented by the following formula. :

[0168]

[0169] in, This represents the distance from the Gaussian beam to the optical axis. The energy at the beam center point is so small that it can be considered constant within the Rayleigh distance range.

[0170] Within the Rayleigh distance range, beam variation is negligible and can be considered constant, thus simplifying model calculations.

[0171] Step S3: Use the maximum likelihood maximum expectation strategy to perform imaging reconstruction on the Gaussian fitting result of the terahertz beam signal to obtain the imaging reconstruction result.

[0172] In step S3, two-dimensional scanning imaging is performed based on the MLEM algorithm. The terahertz MLEM algorithm is an iterative reconstruction method based on maximum likelihood estimation, used to recover high-quality two-dimensional terahertz scanning images under conditions of strong noise, finite aperture, and beam dispersion. In one embodiment of the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle of the present invention, step S3 may specifically include:

[0173] Step S31: Construct a projection model to express the imaging problem as a matrix equation.

[0174] Specifically, the scan data is represented by a matrix equation as follows:

[0175]

[0176] middle, Represents the column vector of the original image. Represents the weight matrix, This represents the collected data, where n represents the noise level. Let R represent a projection weight matrix of size M×N, where M represents the total number of rays, N represents the total number of pixels in the image, and each element represents the length of the ray passing through that pixel. Let y represent the recorded data, which is of size M×1. Since each ray can only pass through a few pixels, R is a sparse matrix. In practical implementations, this equation is difficult to solve directly due to ill-conditioned and noise issues, and requires an iterative algorithm.

[0177] Step S32: Perform MLEM iterative reconstruction based on the projection model to iteratively solve the image.

[0178] In one embodiment of the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle of the present invention, step S32 may specifically include:

[0179] Step S321: In the imaging model, the response of the terahertz system is expressed by the following formula:

[0180]

[0181] in, This represents the measured signal at the scan point. This represents the spatial reflectance distribution to be restored. Let represent the system matrix containing the terahertz beam point spread function (PSF).

[0182] The energy or amplitude measured during terahertz detection typically follows a Poisson distribution. MLEM achieves a progressive approximation of the true reflection / transmission distribution by constructing a forward physical model of the terahertz scanning system and maximizing the probability of the measured data under this model. PSF is usually represented by a Gaussian beam model to characterize the effects of numerical aperture, frequency, and propagation dispersion on the beam spot, thereby accurately describing the energy diffusion behavior of terahertz waves on the sample surface.

[0183] Step S322: Based on the EM derivation using Poisson likelihood, the update formula for terahertz MLEM is obtained as follows:

[0184]

[0185] Among them, forward projection It is used to calculate the predicted measurement, and the back projection is used to realize the error backpropagation, and through the normalization term Ensure the stability and non-negativity of updates.

[0186] Due to the system matrix By incorporating the frequency dependence characteristics, spot diffusion, and propagation attenuation models of terahertz Gaussian beams, MLEM can dynamically compensate for the diffusion effect of terahertz waves in each iteration, allowing the reconstructed image to gradually converge to the true reflection distribution. Compared to traditional deconvolution or interpolation methods, terahertz MLEM has stronger noise robustness and spatial resolution recovery capabilities, making it particularly suitable for terahertz two-dimensional scanning imaging under conditions of large spot size, small target, or low signal-to-noise ratio.

[0187] Step S323: Perform MLEM iteration according to the update formula of terahertz MLEM.

[0188] In one embodiment of the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle of the present invention, step S323 may specifically include:

[0189] Step S3231, Expectation Step (E-step): In the current image estimation... Next, calculate the expected value of the observed data with respect to the complete likelihood function, and obtain the result for each measurement. Belongs to each pixel The probability contribution.

[0190] In one embodiment of the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle of the present invention, step S3231 may specifically include:

[0191] Step S32311: Calculate the predicted measurement using the following formula:

[0192]

[0193] Step S32312: Calculate the attribution coefficient using the following formula:

[0194]

[0195] The attribution coefficient represents the ratio of measurement to prediction. This ratio reflects the deviation between the current model prediction and the actual measurement. It is a weighting factor for the back projection and is a simplified result of the following expectation:

[0196] .

[0197] In this invention, the role of the E-step is to estimate the contribution of each measurement value to each pixel, providing statistical weights for the M-step.

[0198] Step S3232, Maximization Step (M step): Using the statistical weights obtained from the expectation step, maximize the expected likelihood function to update the image estimate.

[0199] Specifically, the update formula for MLEM is expressed by the following formula:

[0200]

[0201] in, This indicates that the measurement error is back-projected back to pixel j according to the system response matrix; the denominator is... It is a normalization term to ensure that the updated value is stable and remains non-negative.

[0202] This M-step ensures that the likelihood function increases monotonically in each iteration, causing the image estimate to gradually approach the maximum likelihood solution.

[0203] Step S33: Introduce deconvolution operation during the iteration process to perform super-resolution optimization.

[0204] Step S33 embodies the deconvolution processing mechanism based on Gaussian beam prior, and the resolution can be further improved through super-resolution optimization. In one embodiment of the terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle of the present invention, step S33 may specifically include:

[0205] Step S331, Convolutional Modeling: When using the current image estimate to generate the forward projection, it is expressed by the following formula:

[0206]

[0207] This is represented as a convolution of the image with a terahertz PSF to simulate the spread blurring of a real beam;

[0208] Step S332, Calculate the correction factor: Based on the ratio of the measured value to the convolutional predicted value, the following formula is obtained:

[0209] .

[0210] Step S333: Perform backprojection and deconvolution update based on the correction factor:

[0211]

[0212] The back-projection part is equivalent to deconvolving the correction factor back to the pixel position. Therefore, the MLEM update process is essentially an adaptive deconvolution.

[0213] Gaussian beam prior information introduces the system's true PSF into matrix A, achieving physically consistent convolution / deconvolution; it compensates for the large spot size and diffusion effect of terahertz waves; it gradually sharpens the image in each iteration, restoring the high-frequency structure; and it effectively improves the spatial resolution of two-dimensional imaging.

[0214] Step S34: Output the reconstructed terahertz image to obtain a high-resolution two-dimensional image.

[0215] The output high-resolution two-dimensional images can clearly reflect the internal structure of the sample.

[0216] In this invention, a hierarchical dependency exists between the beam model, projection model, and MLEM reconstruction in terahertz two-dimensional imaging. First, a terahertz beam model is established using Gaussian beam propagation theory to obtain the system's point spread function (PSF) and its energy distribution on the pixel plane; this distribution is used as a priori to construct the projection matrix. This allows it to simultaneously incorporate geometric scanning relationships and beam intensity weighting, thereby accurately describing the beam intensity from the actual object. To measurement data The forward imaging process is then performed. Subsequently, MLEM iteratively utilizes this projection model for statistical inversion, continuously updating the image estimate through data consistency constraints and PSF prior constraints. The introduction of PSF incorporates convolution / deconvolution compensation during the iteration process, which can suppress terahertz beam dispersion to a certain extent and achieve super-resolution reconstruction. Ultimately, the imaging result is jointly determined by the physical beam model, the accuracy of the projection matrix, and the iterative optimization of MLEM.

[0217] In one embodiment of the present invention, the terahertz beam signal was actually measured as shown in Figure 3, and Gaussian fitting was performed, the result of which is shown in Figure 4. In Figure 4, the blue dashed line represents the actual measured beam signal. Resolution is generally expressed using the full width at half maximum (FWHM), and analysis shows that the beam's FWHM is approximately 1 mm. The red solid line represents the Gaussian fitting result, and it can be observed that the fitting result is very close to the measured value.

[0218] Furthermore, to verify the effectiveness of the maximum likelihood algorithm (MLEM), transmission imaging was performed on fresh leaves as targets in some embodiments of the present invention. The leaves were approximately 0.5 mm thick and placed at the focal point on a plane perpendicular to the terahertz wave optical path. Data was collected from the leaves at 1 mm sampling intervals, and the direct scanning imaging results are shown in Figure 5. High-resolution imaging of the leaves was then performed using the method of the present invention, and the maximum likelihood algorithm imaging results are shown in Figure 6. It can be observed that the original scanning imaging results shown in Figure 5 can only roughly distinguish the main vein structure of the leaf, while the branch veins are difficult to distinguish. The maximum likelihood algorithm imaging results shown in Figure 6 demonstrate that the imaging results using the maximum likelihood algorithm can clearly distinguish all vein structures and effectively suppress artifacts and ringing phenomena.

[0219] Furthermore, to verify the effectiveness of the maximum likelihood algorithm (MLEM) in tobacco product detection, in some embodiments of this invention, transmission imaging was performed on finished cigarettes as targets. The agglomeration imaging analysis results of the finished cigarettes are shown in Figures 7 and 8. In terahertz imaging data, high-density areas have low signal intensity, while low-density areas have high signal intensity. To facilitate density distribution analysis, this invention sets high-density areas in red and low-density areas in blue. Figure 7 shows the terahertz amplitude and phase signal distribution of the direct scanning imaging of the finished cigarette. Analysis of Figure 7 reveals that the original scanning imaging results can only roughly distinguish the overall outline structure inside the cigarette. The fine fiber distribution and material agglomeration inside are blurred and artifacts to a certain extent, making accurate identification impossible. Figure 8 shows the amplitude and phase signal distribution of the maximum likelihood algorithm imaging of the finished cigarette. Analysis of Figure 8 reveals that the imaging results using the maximum likelihood algorithm can clearly distinguish the distribution of fibers and dense areas inside the cigarette, highlighting detailed structures. At the same time, artifacts and ringing effects are effectively suppressed, and the imaging quality is significantly better than that of the direct scanning imaging method. Furthermore, in Figure 8, denser agglomerated areas (such as tobacco clumps) are shown as red high-density blocks, which contrast sharply with the surrounding blue low-density areas.

[0220] After the sample to be tested was prepared, the cigarette was subjected to terahertz scanning to obtain its transmission imaging data, and the results were analyzed, as shown in Figure 9. By performing terahertz scanning on a cigarette with artificial discontinuity defects, the imaging results shown in Figure 9 clearly reveal the discontinuous regions inside the cigarette and their spatial distribution characteristics. Compared with a complete cigarette, the discontinuous regions shown in Figure 9 exhibit significant changes in transmission intensity, verifying the sensitivity and feasibility of terahertz imaging in detecting internal structural anomalies in cigarettes. Experimental results show that the terahertz imaging analysis method used in this invention can effectively reveal the discontinuous characteristics of cigarettes, providing a theoretical basis for the non-destructive testing of the internal structure of finished cigarettes, and verifying the feasibility and practicality of the proposed analysis method.

[0221] The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle provided in this invention effectively reduces path length calculation errors and improves scanning imaging quality by introducing a reconstruction technique combining a Gaussian beam model and the MLEM algorithm, achieving high-resolution two-dimensional terahertz imaging. By introducing a Gaussian beam model into the imaging reconstruction process and combining iterative optimization with the maximum likelihood maximum expectation algorithm, the beam dispersion effect in actual terahertz wave propagation is effectively compensated. A projection weight matrix containing prior beam knowledge is constructed, and deconvolution operation is performed during the MLEM iteration process, significantly improving image resolution and signal-to-noise ratio. The method of this invention can clearly present the structural features of internal fiber distribution, agglomeration, and discontinuity in cigarettes, overcoming the shortcomings of traditional back-projection algorithms such as image blurring and numerous artifacts, providing an efficient and reliable solution for non-destructive testing and quality assessment of the internal structure of cigarette products.

[0222] The embodiments of this disclosure have now been described in detail. To avoid obscuring the concept of this disclosure, some details known in the art have not been described. Those skilled in the art can fully understand how to implement the technical solutions disclosed herein based on the above description.

[0223] While specific embodiments of this disclosure have been described in detail by way of examples, those skilled in the art should understand that the examples are for illustrative purposes only and not intended to limit the scope of this disclosure. Those skilled in the art should understand that modifications can be made to the above embodiments or equivalent substitutions can be made to some technical features without departing from the scope and spirit of this disclosure. The scope of this disclosure is defined by the appended claims.

Claims

1. A terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle, characterized in that, include: The terahertz beam signal of the sample under test is obtained by a terahertz time-domain spectroscopy system; the terahertz beam signal is fitted with Gaussian using a pre-constructed terahertz beam model; the Gaussian fitting result of the terahertz beam signal is reconstructed by using the maximum likelihood maximum expectation strategy to obtain the imaging reconstruction result.

2. The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle according to claim 1, characterized in that, The acquisition of the terahertz beam signal of the sample under test using a terahertz time-domain spectroscopy system includes: the output pulse of a femtosecond laser source is split into pump light and probe light by a beam splitter; the pump light emitted by the terahertz emission unit is focused and then irradiates the low-temperature grown gallium arsenide photoconductive antenna to generate a broadband terahertz pulse; the generated terahertz pulse is collimated by a lens and transmitted to the sample under test through the terahertz transmission optical path to form a terahertz spot signal; the terahertz detection unit, in which the probe light and the transmitted terahertz wave are collinearly incident on a ZnTe electro-optic crystal, converts the terahertz electric field intensity into a polarization state change through the Pockels effect, and obtains the signal through beam splitting by a Wollaston prism and detection by a balanced detector; the optical path of the probe optical path is adjusted by a delay control module to achieve full waveform sampling; and the time-domain signal is Fourier transformed by a data processing module to obtain frequency domain information, thus obtaining the terahertz beam signal of the sample under test.

3. The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle according to claim 1, characterized in that, The method of using a pre-constructed terahertz beam model to perform Gaussian fitting on the terahertz beam signal includes: during the scanning imaging process, representing the received signal by the convolution of the transmitted waveform and the sample function; representing the spot diameter by the beam waist radius; and representing the energy distribution by the Rayleigh distance.

4. The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle according to claim 3, characterized in that, The method of representing the received signal by the convolution of the transmitted waveform and the sample function during the scanning imaging process includes: representing the received signal during the scanning process by the following formula: in, This indicates the result of the scanning imaging. Let n represent the function of the sample to be tested, and n represent the additive noise. The beam model represents the point spread function, where * denotes the convolution operation. In terahertz transmission imaging, the beam passes through the sample under test along the z-axis. Since the sample is three-dimensional, the precise three-dimensional model of the imaging process is represented by the following formula: in, Indicates the location of the terahertz transmitter. The position of the terahertz detection end is indicated, and the accurate model of the PSF is expressed by the following formula. ： Where NA represents the numerical aperture and k represents the cutoff factor. This represents the distance from the sample to the z-axis. Indicates the attenuation coefficient. The frequency at the focal plane is . The beam spot function, wherein the beam waist radius represents the beam spot diameter, includes the following formula: the waist radius is expressed as... The distance from the focal point is Gaussian beam spot diameter ： in, Representing the Rayleigh distance; the representation of energy distribution using the Rayleigh distance includes: representing the Gaussian beam energy distribution using the following formula. ： in, This represents the distance from the Gaussian beam to the optical axis. This represents the energy at the center point of the beam.

5. The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle according to claim 4, characterized in that, The beam model of the point spread function This describes the relationship between numerical aperture, attenuation coefficient, spot function, and frequency, wherein the beam model of the point spread function... The construction process includes: In a terahertz two-dimensional scanning imaging system, the system's imaging response can be considered as a linear space-invariant model, and the point spread function (PSF) is approximated by a Gaussian beam, expressed as: in: Indicates terahertz frequency, The frequency-dependent spot radius reflects the scattering and dispersion characteristics of the beam in space; the minimum spot radius of the Gaussian focal spot formed by a terahertz beam through a focusing lens or reflector is expressed by the following formula. The physical relationship between the numerical aperture (NA) of the system and the system's numerical aperture (NA): in: The wavelength of a terahertz beam decreases as the frequency increases. A larger NA (nanometer) results in a smaller spot size, a sharper PSF (photosensitive element), and higher resolution. The minimum spot radius is determined by the Gaussian focal spot formed by the terahertz beam passing through a focusing lens or reflector. The physical relationship between the PSF and the numerical aperture (NA) of the system is expressed by the following formula, which represents the frequency dependence of the PSF: In terahertz systems, propagation distance This leads to spot expansion, and the propagation distance dependence of PSF can be expressed by the following formula: Rayleigh length is expressed by the following formula: The attenuation coefficient of terahertz waves in a medium due to absorption loss is expressed by the following formula: in, The absorption coefficient is related to frequency; the higher the frequency, the stronger the attenuation. The PSF expression after propagation, expressed in its complete form with an attenuation factor, is given by the following formula: in, This indicates that PSF is determined by frequency. Numerical aperture (NA), spot function and attenuation coefficient A joint decision.

6. The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle according to claim 1, characterized in that, The method of using the maximum likelihood maximum expectation strategy to perform imaging reconstruction on the Gaussian fitting results of the terahertz beam signal to obtain the imaging reconstruction results includes: constructing a projection model to express the imaging problem as a matrix equation; performing MLEM iterative reconstruction based on the projection model to iteratively solve the image; introducing deconvolution operation during the iteration process to perform super-resolution optimization; and outputting the reconstructed terahertz image to obtain a high-resolution two-dimensional image.

7. The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle according to claim 6, characterized in that, The construction of the projection model to express the imaging problem as a matrix equation includes: the scan data is represented by the matrix equation as follows: middle, Represents the column vector of the original image. Represents the weight matrix, This represents the collected data, where n represents the noise level. Let y represent the projection weight matrix of size M×N, where M represents the total number of rays, N represents the total number of image pixels, each element value represents the length of the ray passing through that pixel, and y represents the data collected, which is of size M×1.

8. The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle according to claim 6, characterized in that, The MLEM iterative reconstruction based on the projection model, to iteratively solve for the image, includes: in the imaging model, the response of the terahertz system is expressed by the following formula: in, This represents the measured signal at the scan point. This represents the spatial reflectance distribution to be restored. To represent the system matrix containing the terahertz beam spread function, the update formula for the terahertz MLEM is derived based on the Poisson likelihood EM derivation as follows: Among them, forward projection It is used to calculate the predicted measurement, and the back projection is used to realize the error backpropagation, and through the normalization term Ensure the stability and non-negativity of the updates; perform MLEM iterations according to the update formula of terahertz MLEM.

9. The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle according to claim 8, characterized in that, The MLEM iteration based on the terahertz MLEM update formula includes: Expectation step: in the current image estimation... Next, calculate the expected value of the observed data with respect to the complete likelihood function, and obtain the result for each measurement. Belongs to each pixel The probability contribution; maximization step: using the statistical weights obtained from the expectation step, maximize the expected likelihood function to update the image estimate, wherein the expectation step includes: calculating the predicted measurement using the following formula: The attribution coefficient is calculated using the following formula: The attribution coefficient represents the ratio of measurement to prediction. This ratio reflects the deviation between the current model's prediction and the actual measurement. It is a weighting factor for back projection and is a simplified result of the following expectation: The maximization step includes: expressing the update formula for MLEM using the following formula: in, This indicates that the measurement error is back-projected back to pixel j according to the system response matrix; the denominator is... It is a normalization term to ensure that the updated value is stable and remains non-negative.

10. The terahertz time-domain spectral imaging method based on the femtosecond laser pump-probe principle according to claim 6, characterized in that, The introduction of deconvolution operations during the iteration process for super-resolution optimization includes: Convolution modeling: When using the current image estimate to generate the forward projection, it is expressed by the following formula: This is represented as the convolution of the image with a terahertz PSF to simulate the diffusion blur of a real beam; the correction factor is calculated using the following formula based on the ratio of the measured value to the convolutional predicted value: Backprojection and deconvolution updates are performed based on the correction factor: The back-projection part is equivalent to deconvolving the correction factor back to the pixel position.

Images