A High Precision Terahertz Wave Image Reconstruction Algorithm

With the development of terahertz (THz) technology, the applications of this spectrum have become increasingly wide-ranging, in areas such as non-destructive testing, security applications and medical scanning, in which one of the most important methods is imaging. Unlike remote sensing applications, THz imaging features sources of array elements that are almost always supposed to be spherical wave radiators, including single antennae. As such, well-developed methodologies such as Range-Doppler Algorithm (RDA) are not directly applicable in such near-range situations. The Back Projection Algorithm (BPA) can provide products of high precision at the the cost of a high computational burden, while the Range Migration Algorithm (RMA) sacrifices the quality of images for efficiency. The Phase-shift Migration Algorithm (PMA) is a good alternative, the features of which combine both of the classical algorithms mentioned above. In this research, it is used for mechanical scanning, and is extended to array imaging for the first time. In addition, the performances of PMA are studied in detail in contrast to BPA and RMA. It is demonstrated in our simulations and experiments described herein that the algorithm can reconstruct images with high precision.


Introduction
THz imaging technology is widely used in nondestructive testing (NDT) [1][2][3], security screening [4][5][6][7][8][9][10], and medical scanning [11][12][13][14]. One typical THz imaging mode is based on the time-domain spectroscopic (TDS) approach, which is predicated on simultaneous determination of the field amplitude, phase and polarization [15]. It has been widely used, and its performance has been highly enhanced by researchers over the years [16,17]. In order to improve the spatial resolution to that beyond the realm of the diffraction limit, THz near-field technology with detecting probes of subwavelength aperture size should be considered [15,18,19]. However it may not be practical in some situations, such as in human screening scenarios, where, a detector array would be a workable alternative. In particular, charge-coupled devices (CCDs) in the THz range (5)(6)(7)(8)(9)(10)(11)(12)(13) [20], and bolometers [21] with much higher responsivities in the THz frequency range are available and can be configured to suit specific imaging applications. Reading out online will be allowed in the form of single-shot recording at a pixel size of a few micrometers. Thus it will be preferred for real-time systems. However, the budget and complexity of a large-sized system are invariably determining factors in some applications, e.g., in a human THz imager with a size of 1 mˆ0.5 m, such a detector array might not be practical. Another mode pertaining to conventional optics and microwaves [22] is readily deployable without drastic modifications to realize image reconstruction: an antenna array could be utilized whose elements can be reduced substantially by the employment while being free of interpolation errors. It has been applied in image reconstruction of targets under the illumination of THz Gaussian beams [43,44], as such, it will be reasonable for quasi-optical systems. However when it comes to the scenario of the multistatic case, the angular divergence of the beam, which is 120 degrees in our work, will be much larger than 30 degrees considered previously [45]. Thus the precision of Gaussian beam approximation will no longer obtain. In the present study, PMA is applied to THz imaging reconstruction under the monostatic as well as the spherical wave assumption. PMA is then extended to the multistatic case for the first time, and made available in the simulation. Finally, the performance of PMA is studied in comparison with BPA and RMA in detail in the simulations and experiments.

Monostatic Case
In order to reconstruct a 3D image, one could follow a number of pathways, one of which being to apply a 2D scanning aperture. The simplest scenario can be a planar scanning aperture as depicted in Figure 1. The dashed grid forms the antenna track during scanning, while each node indicates a sample position, which means equal interval sampling along y-axis and z-axis. Targets to be scanned shall be located in the dashed rectangular box area, which will be the volume to be focused.
Sensors 2016, 16,1139 3 of 18 reconstruction of targets under the illumination of THz Gaussian beams [43,44], as such, it will be reasonable for quasi-optical systems. However when it comes to the scenario of the multistatic case, the angular divergence of the beam, which is 120 degrees in our work, will be much larger than 30 degrees considered previously [45]. Thus the precision of Gaussian beam approximation will no longer obtain. In the present study, PMA is applied to THz imaging reconstruction under the monostatic as well as the spherical wave assumption. PMA is then extended to the multistatic case for the first time, and made available in the simulation. Finally, the performance of PMA is studied in comparison with BPA and RMA in detail in the simulations and experiments.

Monostatic Case
In order to reconstruct a 3D image, one could follow a number of pathways, one of which being to apply a 2D scanning aperture. The simplest scenario can be a planar scanning aperture as depicted in Figure 1. The dashed grid forms the antenna track during scanning, while each node indicates a sample position, which means equal interval sampling along y-axis and z-axis. Targets to be scanned shall be located in the dashed rectangular box area, which will be the volume to be focused. As mentioned above, the image reconstruction process is naturally an inverse problem of electromagnetic field scattering. Therefore Helmholtz equation solved with the Born approximation is a primary part in matched filtering algorithms. In addition, it is found that if attenuation in free space is ignored, it will not impact too much in the image reconstructed. This approximation could be reasoned as follows: if the distance between the source and the observation points do not vary too much (z is a constant), only variations on the x-y plane in the object volume are taken into account. Additionally, the phase factors of the echoes play a major role when coherent summations are performed in our imaging scenario. Of course, the final product will just be a reference image to some extent, and each voxel presents the complex reflection coefficient of a position. Suppose the antenna is positioned at ( , , ), while a scatter point of the target is positioned at ( , , ). Thus the objective of the imaging task is to reconstruct the reflection coefficient of ( , , ), which is presented as ( , , ). Here the antenna positions form a planar array, and it is always assumed to be parallel to one of the coordinate planes, whose third coordinate component is supposed to be here. According to the assumptions made above, the system response can be expressed as [6]: where k presents the wavenumber of the signal used, and X, Y and Z span the domain to be imaged, which is represented by the dashed rectangular box in As mentioned above, the image reconstruction process is naturally an inverse problem of electromagnetic field scattering. Therefore Helmholtz equation solved with the Born approximation is a primary part in matched filtering algorithms. In addition, it is found that if attenuation in free space is ignored, it will not impact too much in the image reconstructed. This approximation could be reasoned as follows: if the distance between the source and the observation points do not vary too much (z is a constant), only variations on the x-y plane in the object volume are taken into account. Additionally, the phase factors of the echoes play a major role when coherent summations are performed in our imaging scenario. Of course, the final product will just be a reference image to some extent, and each voxel presents the complex reflection coefficient of a position. Suppose the antenna is positioned at px, y, z L q, while a scatter point of the target is positioned at px, y, zq. Thus the objective of the imaging task is to reconstruct the reflection coefficient of px, y, zq, which is presented as p px, y, zq. Here the antenna positions form a planar array, and it is always assumed to be parallel to one of the coordinate planes, whose third coordinate component is supposed to be z L here. According to the assumptions made above, the system response can be expressed as [6]: s`x 1 , y 1 , k˘" x,y,z p px, y, zq e´j where k presents the wavenumber of the signal used, and X, Y and Z span the domain to be imaged, which is represented by the dashed rectangular box in Figure 1. Here the antenna pattern has not been considered. This is because antenna pattern and attenuation in free space will affect the amplitude of each voxel in approximately the same manner, whereas the phase of each voxel should be the main differentiator in this algorithm, just like in the case of synthetic aperture technology. A 2D Fourier transformation is performed to both sides of Equation (1), after which we get: x,y p px, y, zq e´j k x 1 x e´j k y 1 y e´j k z pz´z L q dxdy where S´k x 1 , k y 1 , k¯is the 2D Fourier transformation(2DFT) of s`x 1 , y, k˘relative to the position variables x and y, and P´k x 1 , k y 1 , z¯is the 2DFT of p px, y, zq relative to x and y. k z is the wavenumber which is introduced to simplify the formula, and it is expressed as: According to the theory of matched filter, Equation (2) can be transformed to: So if the signal and the echo of frequency domain is known, the reflection coefficient function p px, y, zq will be obtained from: p px, y, zq " k x 1 ,k y 1 P´k x 1 , k y 1 , z¯e jk x 1 x e jk y 1 y dk x 1 dk y 1 Perhaps it has been noticed that S´k x 1 , k y 1 , k¯cannot be obtained from the echoes directly. It will be natural to conclude that S`x 1 , y, k˘can be expressed as: wherein E`x 1 , y 1 , k˘refers to the received signal at the location`x 1 , y 1˘i n wavenumber domain, and o pkq refers to the transmitted signal in wavenumber domain. Equation (4) is a critical step to complete the pulse compression along the range direction. If discrete signals are obtained, Equation (4) will be realized in the form: In Equation (7), k i is the sample in wavenumber domain. ∆k is the sample interval of wavenumber k.
In summary, the whole algorithm can be realized following the steps as shown in Figure 2. Firstly, the signals received should be transformed to frequency (wavenumber) domain. The echoes are often baseband complex signals. If only real signals are available, Hilbert Transformation should be performed. Then according to Equation (2), 2DFT will be performed along both cross range directions. If the excitation source is not a Dirac delta function, which is always the case, the system response must be calculated from Equation (6). The pulse compression along the range direction is realized by Equation (7). In fact, ∆k is an irrelevant factor in the process that can be ignored. Finally, 2DFT and 2DIFT will be realized by FFT.

Multistatic Case
In the foregoing consideration, only one antenna is applied to scan targets mechanically. It will take much longer time than electronic scanning, which can be realized by Multiple Input Multiple Output (MIMO) or Phased Array (PA) technology. When MIMO technology is used, the corresponding scenario called the multistatic case here, the amount of antenna elements are greatly reduced than the monostatic case. The array geometry, which is called a plus array [46], and the setup of the experiment, are demonstrated schematically in Figure 3.  Suppose the locations of an arbitrary transmitting element and its corresponding receiving element are ( , , ), and ( , , ), respectively, as shown in Figure 3b, similar to Equation (1), the system response in the multistatic case of a scattering problem can be described by:

Multistatic Case
In the foregoing consideration, only one antenna is applied to scan targets mechanically. It will take much longer time than electronic scanning, which can be realized by Multiple Input Multiple Output (MIMO) or Phased Array (PA) technology. When MIMO technology is used, the corresponding scenario called the multistatic case here, the amount of antenna elements are greatly reduced than the monostatic case. The array geometry, which is called a plus array [46], and the setup of the experiment, are demonstrated schematically in Figure 3.

Multistatic Case
In the foregoing consideration, only one antenna is applied to scan targets mechanically. It will take much longer time than electronic scanning, which can be realized by Multiple Input Multiple Output (MIMO) or Phased Array (PA) technology. When MIMO technology is used, the corresponding scenario called the multistatic case here, the amount of antenna elements are greatly reduced than the monostatic case. The array geometry, which is called a plus array [46], and the setup of the experiment, are demonstrated schematically in Figure 3.  Suppose the locations of an arbitrary transmitting element and its corresponding receiving element are ( , , ), and ( , , ), respectively, as shown in Figure 3b, similar to Equation (1), the system response in the multistatic case of a scattering problem can be described by:  Suppose the locations of an arbitrary transmitting element and its corresponding receiving element are px t , y t , z t q, and px r , y r , z r q, respectively, as shown in Figure 3b, similar to Equation (1), the system response in the multistatic case of a scattering problem can be described by: In Equation (8), the reflection coefficient function p px, y, zq is the objective to be imaged. Compared to Equation (1), it is found that the system response is a 5D matrix. Different groups of transmitting and receiving antennas increase the dimensionality of data, which can be interpreted as the increase of degrees of freedom. That of course is one of the salient features of MIMO technology.
Then 4DFT will be performed on both sides of Equation (8) along the four cross range directions, thereafter we obtain: where the operator FT 4D r. . .s means 4DFT of the first four dimensions of the data, with the functions: As g t px t , y t q and g r px r , y r q have no coupling between their integral variables, it follows that: The functions G t`kxt , k yt , k˘and G r`kxr , k yr , k˘can be expressed as: Substitute Equations (12)-(14) into Equation (9), it is recast in the form: x,y,z p px, y, zq e´j k x x e´j k y y e´j k z z dxdydz For convenience, a new set of variables k x , k y and k z are introduced, defined as: If P`k x , k y , z˘expresses the 2DFT of p px, y, zq, it will be given as: x,y p px, y, zq e´j k x x e´j k y y dxdy Substituting Equation (19) into Equation (15), and we will get: According to the theory of matched filter, the formulas: and: p px, y, zq " k x ,k y P`k x , k y , z˘e jk x x e jk y y dk x dk y (22) together reconstruct the final 3D image.
The corresponding block diagram of the process can be summarized as shown in Figure 4. Compared to Figure 2, two main issues should be noticed:

1.
Phase compensations will be taken according to Equation (21), and when the array is a planar one, which is the case in Figure 3b it is always assumed z t " 0 and z r " 0 for convenience; 2.
Dimensionality reduction means S`k xt , k yt , k xr , k yr , k˘will be transformed into S`k x , k y , k˘, where a 5D data matrix is reduced to a 3D one. This is performed according to Equations (16) and (17).  (22) together reconstruct the final 3D image.
The corresponding block diagram of the process can be summarized as shown in Figure 4. Compared to Figure 2, two main issues should be noticed: 1. Phase compensations will be taken according to Equation (21), and when the array is a planar one, which is the case in Figure 3b it is always assumed = 0 and = 0 for convenience; 2. Dimensionality reduction means ( , , , , ) will be transformed into ( , , ), where a 5D data matrix is reduced to a 3D one. This is performed according to Equations (16) and (17).

Discussion on Sampling Criteria and Spatial Resolution
According to the Nyquist sampling theorem, the sensors must sample "close" enough to cover the support region of the wavenumber components. Otherwise, it will suffer from spectrum aliasing in frequency domain.
The derivation of the sampling criteria is presented in Appendix A. According to Equations (A5)-(A8), the spatial sampling interval must satisfy:

Discussion on Sampling Criteria and Spatial Resolution
According to the Nyquist sampling theorem, the sensors must sample "close" enough to cover the support region of the wavenumber components. Otherwise, it will suffer from spectrum aliasing in frequency domain. The derivation of the sampling criteria is presented in Appendix A. According to Equations (A5)-(A8), the spatial sampling interval must satisfy: where u can be replaced by x or y and s can be replaced by t or r. In Equation (9), L us refers to the array's aperture size of the transmitting/receiving elements in the u-direction, and R 0 is the closest distance between the antenna array and the target's front surface. Of course, if equivalent elements are considered, the sampling requirement in the same direction may become less stringent. The spatial resolution for the multistatic case is discussed in Appendix A, and can be summarized as follows: where u represents the cross range direction, which can be replaced by x or y. λ c is the center wavelength. Here the sampling criteria and resolution are discussed only for the multistatic case. As the monostatic case can be deduced easily from the equations above, or be found in other papers [6,47], it will not be discussed here.

Simulation Results
In this section, two sample simulations will be undertaken to prove the effectiveness of our algorithm. To describe its performance objectively, two classical algorithms, RMA and BPA, will be introduced to provide contrasts. In the following, the coordinate system adopted in the previous section will be modified, such that range direction is set to be parallel to the x-axis to adapt to the settings in the electromagnetic field simulation software.

Monostatic Case
In this case, electromagnetic field simulation data of specific scenarios, calculated with the method of moment (MoM), will be utilized to test the performance of the three algorithms, including their applicability and efficiency. The parameters used in these two simulations are listed in Table 1. Two kinds of targets have been used to reconstruct images as demonstrated by Figure 5. The first one aims to test the performance of positioning and focusing, and the second tests the ability of the algorithms to present objects with complex shapes.
In Figure 5a, seven metal points that are extremely small (the radii are about a thirtieth of the center wavelength) are used as ideal scatter point targets. Two kinds of targets have been used to reconstruct images as demonstrated by Figure 5. The first one aims to test the performance of positioning and focusing, and the second tests the ability of the algorithms to present objects with complex shapes.  In Figure 6, the results are demonstrated in the form of slice figures. The first and the third row are obtained from BPA and PMA, respectively. According to their actual coordinates, the slice positions are just where they should be. However, four points appear in error slice positions in Figure 6d In Figure 6, the results are demonstrated in the form of slice figures. The first and the third row are obtained from BPA and PMA, respectively. According to their actual coordinates, the slice positions are just where they should be. However, four points appear in error slice positions in Figure 6d-f. That means RMA's performance of positioning in range direction is poor.
To demonstrate the final results in the form of a 2D image, a function is applied to present 3D complex matrices. The brightness of each pixel represents the maximum modulus of the complex voxels along the range direction, while which ranges from red to blue with increasing distance, of each pixel indicates the position where the maximum modulus is situated. Additionally, noise suppression is introduced in the function with a threshold to filter out low amplitude voxel values, which are supposed to be the noise. In the simulations, the threshold is set to be 0.1 as the noise comes mainly from sidelobes, whereas a higher threshold shall be adopted in treating experimental data. In Figure 7, the results of the algorithms are presented with the plot function. It can be seen that the left two points are nearest to the scanning plane, whose z-coordinate component is 0.03. While the two points on the right side whose color tends to be blue must be farther than the central ones.  To demonstrate the final results in the form of a 2D image, a function is applied to present 3D complex matrices. The brightness of each pixel represents the maximum modulus of the complex voxels along the range direction, while which ranges from red to blue with increasing distance, of each pixel indicates the position where the maximum modulus is situated. Additionally, noise suppression is introduced in the function with a threshold to filter out low amplitude voxel values, which are supposed to be the noise. In the simulations, the threshold is set to be 0.1 as the noise comes mainly from sidelobes, whereas a higher threshold shall be adopted in treating experimental data. In Figure 7, the results of the algorithms are presented with the plot function. It can be seen that the left two points are nearest to the scanning plane, whose z-coordinate component is 0.03. While the two points on the right side whose color tends to be blue must be farther than the central ones. The brightness of each pixel represents the maximum modulus of the complex voxels along the range direction, while the color, which ranges from red to blue with increasing distance, of each pixel indicates the position where the maximum modulus is situated.
Next a metal fan is used to test the ability of each algorithm to present objects with complex shapes. The profile of the fan is shown in Figure 5b. The nearer to the central part of the model, the better resolution will be required to depict it. Therefore it reflects the imaging resolution to some extent. Figures 8 and 9 demonstrate the 3D images reconstructed by the three algorithms. It seems that they have nearly the same resolution across range direction as well as the sidelobe-performance. Finally, the efficiency of each algorithm will be discussed. Here it is presented as the processing time consumed, listed in Table 2. It can be seen that RMA achieves the highest efficiency amongst the three algorithms, while BPA takes too long to be used in a real-time system; PMA performs much The brightness of each pixel represents the maximum modulus of the complex voxels along the range direction, while the color, which ranges from red to blue with increasing distance, of each pixel indicates the position where the maximum modulus is situated.
Next a metal fan is used to test the ability of each algorithm to present objects with complex shapes. The profile of the fan is shown in Figure 5b. The nearer to the central part of the model, the better resolution will be required to depict it. Therefore it reflects the imaging resolution to some extent. Figures 8 and 9 demonstrate the 3D images reconstructed by the three algorithms. It seems that they have nearly the same resolution across range direction as well as the sidelobe-performance. The brightness of each pixel represents the maximum modulus of the complex voxels along the range direction, while the color, which ranges from red to blue with increasing distance, of each pixel indicates the position where the maximum modulus is situated.
Next a metal fan is used to test the ability of each algorithm to present objects with complex shapes. The profile of the fan is shown in Figure 5b. The nearer to the central part of the model, the better resolution will be required to depict it. Therefore it reflects the imaging resolution to some extent. Figures 8 and 9 demonstrate the 3D images reconstructed by the three algorithms. It seems that they have nearly the same resolution across range direction as well as the sidelobe-performance. Finally, the efficiency of each algorithm will be discussed. Here it is presented as the processing time consumed, listed in Table 2. It can be seen that RMA achieves the highest efficiency amongst the three algorithms, while BPA takes too long to be used in a real-time system; PMA performs much Next a metal fan is used to test the ability of each algorithm to present objects with complex shapes. The profile of the fan is shown in Figure 5b. The nearer to the central part of the model, the better resolution will be required to depict it. Therefore it reflects the imaging resolution to some extent. Figures 8 and 9 demonstrate the 3D images reconstructed by the three algorithms. It seems that they have nearly the same resolution across range direction as well as the sidelobe-performance. Finally, the efficiency of each algorithm will be discussed. Here it is presented as the processing time consumed, listed in Table 2. It can be seen that RMA achieves the highest efficiency amongst the three algorithms, while BPA takes too long to be used in a real-time system; PMA performs much Finally, the efficiency of each algorithm will be discussed. Here it is presented as the processing time consumed, listed in Table 2. It can be seen that RMA achieves the highest efficiency amongst the three algorithms, while BPA takes too long to be used in a real-time system; PMA performs much faster than BPA, retaining high precision of BPA, which makes compromise between BPA and RMA.

Multistatic Case
In this section, the performance of PMA for the multistatic case, the procedure of which has been described in detail in Section 2.2, will be tested with electromagnetic field simulation briefly here. The metal fan in Figure 7b is used as the target, and the antenna array is arranged in the geometry shown in Figure 3a. One hundred and one transmitting and receiving elements are used in this scenario, and the operation procedure will be as follows: each transmitting element radiates signals in sequence while all the receiving elements receive signals at the same time. In fact, the plus array imaging is simpler than the case mentioned in Section 2.2, as x t = 0 and y r = 0. However, the principle of the algorithm can be proved without loss of generality.
The results of multistatic array imaging are demonstrated in Figure 10. Compared to Figures 8 and 9, it is found that the performance is on par with that of the monostatic case.

Multistatic Case
In this section, the performance of PMA for the multistatic case, the procedure of which has been described in detail in Section 2.2, will be tested with electromagnetic field simulation briefly here. The metal fan in Figure 7b is used as the target, and the antenna array is arranged in the geometry shown in Figure 3a. One hundred and one transmitting and receiving elements are used in this scenario, and the operation procedure will be as follows: each transmitting element radiates signals in sequence while all the receiving elements receive signals at the same time. In fact, the plus array imaging is simpler than the case mentioned in Section 2.2, as xt = 0 and yr = 0. However, the principle of the algorithm can be proved without loss of generality.
The results of multistatic array imaging are demonstrated in Figure 10. Compared to Figures 8 and 9, it is found that the performance is on par with that of the monostatic case.

Experimental Results
In this section, we present results of the monostatic experiment. Figure 11 demonstrates the setup of the experiment and the target to be imaged. A model N5247A network analyzer (Agilent, Santa Clara, CA, USA) is used to transmit, receive and process the signals. The signal used is a step-frequency continuous wave, and will be received after down-conversion. The signals received are analyzed in frequency domain, where it will be convenient to process. Mechanical raster scanning is realized by a 2D scanning platform. The parameters applied in this experiment are listed in Table 3.

Experimental Results
In this section, we present results of the monostatic experiment. Figure 11 demonstrates the setup of the experiment and the target to be imaged. A model N5247A network analyzer (Agilent, Santa Clara, CA, USA) is used to transmit, receive and process the signals. The signal used is a step-frequency continuous wave, and will be received after down-conversion. The signals received are analyzed in frequency domain, where it will be convenient to process. Mechanical raster scanning is realized by a 2D scanning platform. The parameters applied in this experiment are listed in Table 3.

Experimental Results
In this section, we present results of the monostatic experiment. Figure 11 demonstrates the setup of the experiment and the target to be imaged. A model N5247A network analyzer (Agilent, Santa Clara, CA, USA) is used to transmit, receive and process the signals. The signal used is a step-frequency continuous wave, and will be received after down-conversion. The signals received are analyzed in frequency domain, where it will be convenient to process. Mechanical raster scanning is realized by a 2D scanning platform. The parameters applied in this experiment are listed in Table 3.   The slice figure results in Figure 12 prove the poor positioning performance of RMA again. The slice position agrees well between Figure 12a,c, while an error of 0.08 m happens in RMA. In Figure 13, some notable cases happen in Figure 13b. Firstly, some red artifacts appear around the image. It results from the errors introduced by interpolation. In fact, the poor positioning performance must be relevant to defocus in the range direction from this point of view.  The slice figure results in Figure 12 prove the poor positioning performance of RMA again. The slice position agrees well between Figure 12a,c, while an error of 0.08 m happens in RMA. In Figure 13, some notable cases happen in Figure 13b. Firstly, some red artifacts appear around the image. It results from the errors introduced by interpolation. In fact, the poor positioning performance must be relevant to defocus in the range direction from this point of view. Secondly, the color of the fan in Figure 13b tends to be blue while the other two images appear red. This is also due to the appearance of artifacts. The plot function of the 3D reflectivity coefficients will be scaled to the maximum datum after noises-removing operation. The color in an image does not refer to an absolute range but a relative quantity. Therefore the three images cannot be contrasted  The slice figure results in Figure 12 prove the poor positioning performance of RMA again. The slice position agrees well between Figure 12a,c, while an error of 0.08 m happens in RMA. In Figure 13, some notable cases happen in Figure 13b. Firstly, some red artifacts appear around the image. It results from the errors introduced by interpolation. In fact, the poor positioning performance must be relevant to defocus in the range direction from this point of view. Secondly, the color of the fan in Figure 13b tends to be blue while the other two images appear red. This is also due to the appearance of artifacts. The plot function of the 3D reflectivity coefficients will be scaled to the maximum datum after noises-removing operation. The color in an image does Secondly, the color of the fan in Figure 13b tends to be blue while the other two images appear red. This is also due to the appearance of artifacts. The plot function of the 3D reflectivity coefficients will be scaled to the maximum datum after noises-removing operation. The color in an image does not refer to an absolute range but a relative quantity. Therefore the three images cannot be contrasted directly in color. It simply means that the artifacts appear nearer to the scanning aperture than the target.

Discussion
In this section, some phenomena in the simulations and experiments will be discussed further. It is known that RMA takes the advantage of FFT in all three dimensions and achieves the highest efficiency. However as the samples distribute unequally in the range direction, interpolation called "Stolt interpolation [40]" should be applied to rearrange the spectrum. Thus a sinc interpolation is introduced and the error caused by the sinc interpolation is smaller than the resolution in the range direction in principle. However, as the interpolation samples are limited, slight defocus in the range direction will be introduced after all. And one of the ramifications may be its poor positioning performance. From another point of view, the reason may be explained in the following way: the sinc interpolation is conducted by moving data in wavenumber domain along the range direction with a sinc weighting function, thus the position tends to shift. Another problem is the artifacts appearing in Figure 13b. However, it is noticed that in the results of the simulations ( Figure 9b) the artifacts do not appear. No noise is included in the simulations, while it will certainly appear in the experimental data. As interpolation causes defocus, the signals will be added incoherently in some voxels. Thus peaks may appear in some unexpected positions that are considered as artifacts.
The algorithm of BP is very classical. Its principle is simple and easy to understand. Recently some modified BPAs have been proposed, which aim to accelerate. However, they are not yet suitable for real-time systems. Whereas BPA processes the echoes directly without any approximation, it is supposed to be closest to ideality in resolution. Thus in this article the result of BPA is introduced as a criterion of image quality. Furthermore, it is preferable to mesh the grid more densely which will introduce more computation.
In fact, Phase-shift Migration Algorithm relates to RMA and BPA internally. It applies Equation (7) directly without carrying out interpolation. Thus it retains the performance of BPA, but operates a little slower than RMA. However, PMA can be implemented in parallel, and has the potential of real-time processing.

Conclusions
In this research, the Phase-shift Migration Algorithm (PMA) is proposed in terahertz imaging, the principle of which is derived in both monostatic and multistatic cases. The sampling criteria and spatial resolution evaluation are obtained in multistatic case. Electromagnetic field simulation and experimental results verify the applicability and performance of PMA. Compared with RMA and BPA, the features and drawbacks of PMA are summarized below:

‚
High precision of positioning relative to RMA, proved in the simulations and experiments; ‚ Capability to operate without knowing the distance a priori. In RMA, the distance from the scanning array to the objective scene center must be known before image reconstruction, or the final images may deteriorate; however this is not the case in PMA, as long as the frequency sampling interval satisfies the requirement of the scope of the objective; ‚ Conserve the focusing and resolution performance of BPA in principle without any approximation; ‚ More computation intensive than RMA, though can be implemented in parallel; ‚ Must operate in equal-interval sampling case. This is a fatal restriction, especially in multistatic cases. But BPA will not suffer from the same.
This algorithm has the potential of working in multiple arrays. The principle of the idea can be found in reference [14], and our future work may focus on the algorithm used in the scenario of multistatic case, especially in sparse arrays.
In this appendix we outline the derivation of sampling criteria and spatial resolution. In the array imaging case, the path of the electromagnetic wave is depicted in Figure 3b, and expressed as: px´x t q 2`p y´y t q 2`p z´z t q 2`b px´x r q 2`p y´y r q 2`p z´z r q 2 " R t px, y, zq`R r px, y, zq (A1) Taking the partial derivative of R px, y, zq with respect to x t results in: px´x t q 2`p y´y t q 2`p z´z t q 2 " x t´x R t px, y, zq If a step-frequency continuous wave is used as the excitation source, the phase of the received signal will be: ϕ px t q "´2 π λ R px, y, zq Supposing ∆x is the sample interval, the Nyquist sampling theorem dictates:ˇˇˇB ϕ px r q Bx tˇ∆ x t ď π (A4) Now suppose the planar antenna array's aperture size in the x-direction is L x , and the objective domain's scope in the x-direction is D x . When the transmitting and receiving elements are chosen to be on the edge of the object box in the x-direction, and the scatter point is exactly on the other side of the edge, the left side of Equation (A4) reaches its maximum. Thus the sampling criteria in the x-direction can be obtained as: where L xt refers to the array's aperture size of the transmitting elements in the x-direction, and R 0 is the nearest distance between the antenna array and the target's front surface. In the same manner the other sampling criteria in the x-direction and the y-direction can be concluded: Derivation of spatial resolution: From the spectral point of view, the resolution ρ u in the u-direction (u can be x, y or z) relates to the support domain of wavenumber ∆k u : In this regard, the resolution in the range direction and the lateral directions can be obtained as: where ρ r and ρ a (the subscript a can be x or y) mean the resolution along and across the range direction, while B is the bandwidth of the signal used. c represents the travelling speed of electromagnetic wave in free space. As k x is calculated from: x´x t R t px,y,zq`x´x r R r px,y,zq ı «´2 π λ c 2x´x t´xr R 0 (A12) ∆k x can be obtained as: Thus we will get the lateral resolution along the xand y-direction: as well as: ρ y " λ c R 0 L yt`Lyr (A15)