Active terahertz (THz) imaging systems have been shown to provide excellent capabilities for detecting and imaging hidden threats, contraband, and internal structures^[1,2] in terms of nondestructive testing. A THz system that has the advantage of standoff observation is generally a multireflector continuous-wave system that collects signals of objects at different locations through rotating mirror scanning^[3]. The research on this kind of system mainly focuses on imaging theory^[4,5], the detection method^[6], or an image-processing algorithm, while the imaging response of different targets is mostly based on actual measurement.

The propagation and scattering of THz waves should be considered to predict the detection results of the imaging system. Although the point spread function model has been used to convolve the objective function to obtain simulated THz images^[7], this method generally does not consider the scattering caused by the object itself. In the THz band, the scattering characteristics of rough surfaces have been measured and calculated^[8,9]. Based on ray tracing, the bidirectional reflectance distribution function has been used to simulate THz image of human body with rough surface^[10]. A physical optics-based THz imaging simulator simulates the THz images of metal plates and mannequins^[11]. It is not enough to only consider the rough surfaces scattering for detecting the internal structure of objects. The scattering characteristics of void structures inside the material are predicted by a point source synthesis model^[12] or a ray-tracing model^[13] in the case of vertical specular reflection. There is still a lack of THz imaging simulation that considers both the internal structure and the incidence angle.

In the visible band, the imaging system collects light from different reflection areas, causing the detected characteristics of the measured object to change, which is the same as the imaging results in the THz band^[14]. The visible light scattering field from surface defects of optical components has been simulated by using an electromagnetic model^[15–17]. The scattering field is then converted into a far-field distribution known as a visible image by using the angular spectrum theory. Although the electromagnetic simulation models have been applied to THz band scattering^[18,19], the calculation time required for electrically large problems is very long.

In this paper, the THz images of foams with embedded stripe defects measured by an active bifocal THz imaging system^[20] are rapidly simulated based on geometric optics (GO) and compared with the measurements. Due to the electrically large reflectors and imaging at long standoff distances, angular spectrum theory and ray tracing are used to obtain the incident field of the foam surface and the electric field distribution at the detector. The electric field distribution, which is modulated by ideal embedded stripe defects on the foam surface, is calculated by GO and the finite-element method (FEM), respectively. Section 2 establishes the numerical model of the imaging system. In Section 3, the images of defective foams at 0.22 THz are presented and compared with simulations. Section 4 provides some brief conclusions.

Figure 1(a) shows a photograph of the quasi-optical system, which includes the main mirror (MR), submirror (SR), tertiary feed reflector mirror (TFR), rotating mirror (RM), source (S), detector (D), and beam splitter (BS). The building process and principle of the system that only changes the imaging distance to achieve bifocal imaging are given in the previous work^[20]. Figure 1(b) shows the parameters that describe the system. The standoff distance of 0.22 THz imaging is 8 m, and the field of view (FoV) is approximately $50\mathrm{cm}\times 90\mathrm{cm}$^[14]. The 0.22 THz sources are model AMC666. The homemade 0.22 THz detector has a noise equivalent power (NEP) of ${10}^{-10}\mathrm{W}/{\mathrm{Hz}}^{1/2}$.

The reflector surfaces of the system have an electrically large size, and pixel-by-pixel acquisition needs to be simulated, while the THz wave is a high-frequency electromagnetic wave. Considering the limited computational storage space and simulation efficiency, we simulate THz images of defective foams in two-dimensional space ($x–z$ plane). Since the system has a larger scanning range in the $y_t$ direction, we selected strip defects with a symmetrical distribution in the $y_t$ axis for imaging in Section 3.

The term $k$ is the working wavenumber, and $r(x)$ is the distance from the point on the front plane of the TFR to the radiation source. The incident wave on the front plane of the TFR can be expressed as $$U_{tf}(x)=e^{-ikr(x)}.$$(1)

Considering the actual aperture size, the aperture function is introduced, $$P(x)=\{\begin{array}{cc}1& \mathrm{inside}\mathrm{the}\mathrm{aperture}\\ 0& \mathrm{other}\end{array}.$$(2)

Equation (1) is rewritten as $$U_{tf}(x)=e^{-ikr(x)}·P(x).$$(3)

According to the theory of angular spectrum diffraction, the monochromatic light-field distribution on a certain plane can be regarded as the superposition of monochromatic plane waves propagating in different directions. The amplitudes and phases of each plane wave depend on the argument and amplitude of the angular spectrum. The angular spectrum of the TRF incident wave is $$A(k_x)={\int }_{-\infty }^{+\infty }{U}_{tf}(x)·e^{-i{k}_{x}x}\mathrm{d}x,$$(4)where $k_x$ is the component of the wave vector in the $x$ direction. We trace each plane wave component of the TRF incident wave by GO. The plane wave in the $k_x$ direction starts from points $x_i$ ($i=1,2,3,\dots$) on the incident plane of the TRF and is reflected by the RM, the SR, and the MR in turn before reaching the imaging plane. The rays intersect the imaging plane ($x_t$ axis) at the points $x_{ti}$ with the total optical path of $l(x_{ti})$ and the phase lag of $e^{ikl(x_{ti})}$ ($i=1,2,3,\dots$). For focusing imaging, the intersections $x_{ti}$ are scattered. By interpolating the total optical path $l(x_{ti})$ and initial coordinate $x_i(x_{ti})$ of discretely distributed incident points $x_{ti}$, the total optical path $l(x_t)$ and initial coordinate $x(x_t)$ of continuously distributed incident points on the $x_t$ axis are obtained. At this point, the incident wave in the imaging plane can be expressed as $$U(x_t)={\int }_{-\infty }^{+\infty }A(k_x)·e^{i{k}_{x}x(x_t)}·e^{ikl(x_t)}\mathrm{d}{k}_x.$$(5)

Next, we establish the scattering model of the defective foams, as shown in Fig. 2. Defects are set as ideal rectangular defects, regardless of edge roughness, changes in refractive index, etc. Since an iron plate with strong reflective properties for THz waves was used in the experiment, the reflective surface is set as an ideal electrical conductor. GO calculations of the scattering field are similar to those of the incident field. In Fig. 2(a), the incident point ($x_{\text{foam}}$) and total optical path ($l_{\text{foam}}$) of the scattering field on the foam surface are obtained by geometric optical tracing of monochromatic light in different propagation directions. The scattering field is expressed as $$U_{\text{foam}}(x_{\text{foam}})={\int }_{-\infty }^{+\infty }{A}_i(k_x)e^{i{k}_{x}x(x_{\text{foam}})}{e}^{ik{l}_{\text{foam}}(x_{\text{foam}})}\mathrm{d}{k}_x,$$(6)where $A_i$ is the angular spectrum of the incident fields. Using the FEM to calculate the scattering field distribution of the foam surface, it is necessary to select the perfectly matched layer boundary conditions to eliminate the influence of the reflected waves caused by the truncated boundary, as shown in Fig. 2(b). The area of FEM simulation is $15\mathrm{cm}\times 2.2\mathrm{cm}$. Equation (5) is used to obtain the incident field distribution and the coordinate point $x_{t\max }$ of maximum incident intensity. The horizontal range of the simulation area ranges from $x_{t\max }-7.5\mathrm{cm}$ to $x_{t\max }+7.5\mathrm{cm}$. The distance between the ideal reflecting surface and the lower surface of the foam, and that between the incident field surface and the upper surface of the foam are both 0.1 cm. THz imaging simulation is achieved on a minimum scale by constantly moving the simulation area to mimic point-by-point scanning operations.

Finally, we calculate the electric field distribution at the detector [$U_{\text{detector}}(x)$] analogously to the way that we calculated the incident field. The received light intensity is expressed as $$I\propto {\int }_s{|U_{\text{detector}}(x)|}^2\mathrm{d}s.$$(7)

The foam board is ethylene-vinyl acetate copolymer (EVA) with a size of $60\mathrm{cm}\times 120\mathrm{cm}$ and a thickness of 2 cm. We machined the embedded strip defects on one side of the EVA boards. Defects with both lengths and intervals of 10 cm are symmetrically distributed in the central axis of the EVA boards. Figures 3(a) and 3(b) show EVA boards with embedded horizontal strip defects. The widths of the strip defect are 0.6 cm in Fig. 3(a) and 3 cm in Fig. 3(b). We tested defects with depths of 0.5, 1, and 1.5 cm. The edge of the machined defects is relatively smooth, so we treat it as an ideal rectangular defect. Figure 3(c) shows the defect-free side of the EVA board. The defective side of the EVA board is adhered to the iron plate, which is fixed on the blackboard, as shown in Fig. 3(d). We have adjusted the system magnification to obtain THz images of different measured objects. To compare THz images, we use the maximum value per image for normalization and then display it on a logarithmic scale.

We measured the refractive index ($n_{\mathrm{EVA}}=n_r-i·n_i$) of the EVA boards by the fiber-based THz TDS system^[21,22]. At 0.22 THz, $n_r=1.0144$ and $n_i=1.16\times {10}^{-4}$. The images of defect-free foam boards at 0.22 THz are shown in Fig. 4(a). We compare the intensity at the center axis of the defect THz images with the simulations. Figure 4(b) shows the imaging results and simulations of the defect-free EVA at 0.22 THz. The simulations obtained by GO and FEM are almost the same. There are obvious wavy intensity changes in the imaging light spots, while the intensity of the simulated spot gradually decreases with increasing incident angle. There are two main reasons for the difference between imaging and simulation. One is the position of the imaging plane; the other is the adjustment accuracy of the imaging system. In the experiment, manually adjusting the system and the target position of the standoff distance is very complicated. Small errors in the real position of the target and the adjustment cause changes in the measurement intensity.

Figure 5 shows the THz images of the defective EVA board. The incident angle increases with the distance from the center of the FoV. As shown in Figs. 6(a1) and 6(a2), under the constant aperture, the receiving intensity gradually decreases, resulting in light spots [Fig. 4(b)]. The existence of defects changes the distribution of the light field. Figures 6(b1), 6(c1), 6(d1), and 6(e1) show a partial deviation from the original direction of specular reflection. Reflection deviation mainly occurs near the defect edge, as shown by the black dotted line in Figs. 6(b1) and 6(c1). Therefore, the specular reflection component decreases, and the diffuse reflection component increases. In Figs. 6(b2), 6(c2), 6(d2), and 6(e2), the black box shows the light-field distribution of the 60 cm aperture. The far-field distribution of the defect scattering field is basically not different, except for the coordinate migration. The greater the incidence angle, the greater the diffuse reflection proportion of the collected reflected light. Imaging features change from specular reflection to diffuse reflection.

For the embedded defects, specular reflection imaging with reduced intensity at the defect is inside the light spot, and diffuse reflection imaging with increased intensity at the defect is outside the light spot shown in Fig. 5. When the defect width increases to 3 cm, regardless of the imaging features, the intensity changes at the defect appear as a saddle shape in the $y_t$ direction, which is consistent with the change trend of the intensity peak caused by edge scattering in dark-field imaging^[16]. The intensity value variation, which is affected by the scattering effect at the edge of the defect, is not only related to the incident angle, but also related to the defect size.

Figure 7 compares the 0.22 THz imaging results and the simulations of the horizontal strip defects. The similarity between the simulation and imaging results is very high in the whole FoV. We simulated the resolvable distance of the two 0.5-cm-width defects. Figure 8 shows the GO results of the resolvable distance. In most areas, the discriminability of two defects decreases with the decrease of the distance between defects. Especially in the case of $y_t=\pm 40\mathrm{cm}$ and $y_t=\pm 50\mathrm{cm}$, the measurement results are very similar to Ref. [14]. However, at $y_t=\pm 30\mathrm{cm}$, the discriminability of two defects is the lowest. Previous work on the discriminability measurements of defects ignored the angular variation in diffuse reflection imaging. This simulation method can help us analyze imaging differences and obtain more THz image information. Because $n_{\mathrm{EVA}}$ is very close to 1, the multiple reflections and refractions of THz waves are not considered, and only one transmission is considered in the GO method. The difference between the two simulation methods is small, which shows the feasibility of GO approximate calculation. There are 745 simulated sampling points with an interval of about 0.15 cm. Using FEM requires grid generation, which takes about 6 s for each sampling point, while using GO method requires only one ray tracing, which takes less than 10 min for all sampling points.

In this paper, the diffuse and specular reflection imaging behavior of defective foam is measured and simulated at 0.22 THz. Due to the electrically large size of the reflectors, the theory of angular spectrum diffraction and ray tracing is used to calculate the incident and receiving electric field distribution of the system with a large FoV and long standoff distance. FEM and GO models of defective foams are established to calculate the near-field scattering light distribution. The proposed GO modeling method can accurately and quickly calculate the intensity information of defective foams. More complex properties of the measured object can be simulated by multiple geometric optical tracing. The simulation model theoretically reveals the imaging process, helping to explain the THz images. Further work will build more THz image data sets of measured objects.