Phaseless Terahertz Coded-Aperture Imaging for Sparse Target Based on Phase Retrieval Algorithm

Phaseless terahertz coded-aperture imaging (PL-TCAI) is a novel radar computational imaging method that utilizes the coded aperture and the incoherent detector array to achieve forward-looking and high-resolution imaging without relying on relative motion. In this paper, we propose a more reasonable and compact architecture for the PL-TCAI system and derive the imaging model of PL-TCAI based on the random frequency-hopping signal. Since most phase retrieval algorithms for PL-TCAI utilize only the intensity of echo signals to accurately reconstruct the target, excessive measurement samples are usually required. In order to reduce the number of measurement samples required for imaging, this paper proposes a sparse Wirtinger flow algorithm with optimal stepsize (SWFOS) by using the sparse prior of the target. The specific procedures of the SWFOS algorithm include the support recovery, initialization by truncated spectral method, iteration via gradient descent scheme, hard threshold operation, and stepsize optimization of iteration. Numerical simulations are performed, and the results show that the SWFOS algorithm not only has good performance for the PR problem, but can also sharply reduce the number of measurement samples required for imaging in the PL-TCAI system.


Introduction
Terahertz coded-aperture imaging (TCAI) [1][2][3][4] is a novel high-resolution imaging method that does not rely on any relative motion between the radar and the target. In TCAI, the computer-controlled electronic coded aperture [5][6][7][8] can modulate the amplitude or phase of the transmitting signal to form a space-time two-dimensional random orthogonal radiation field, thereby randomly generating various illumination modes in the imaging area. Random radiation fields can help to improve radar imaging resolution, even though they are accompanied by a significant reduction in detection distance [9]. The azimuth resolution of TCAI is not related to the Doppler frequency, so it does not depend on the relative motion with the target, which can achieve forward-looking and staring imaging [9,10]. For the TCAI system, the reference signal corresponding to the imaging area can be derived and then combined with the receiving signal, and then the matrix imaging equation of TCAI can be established. In recent years, with the increasing demand for forward imaging and high-resolution imaging, TCAI has received more and more attention. For example, since the frequency range of terahertz waves is 0.1-10 THz, TCAI can be used for human security screening at airports, stations, or other important descent method, which directly affects the performance of the algorithms. If the stepsize is too large, the algorithms will diverge easily, while if the stepsize is too small, the algorithms will converge slowly. For the SPARTA, SWF, and Thresholding Wirtinger Flow algorithms, the stepsizes are selected empirically, and the default stepsize functions or stepsize constants are proposed for the measurement vectors that conform to the standard Gaussian random distribution [26][27][28][29][30][31][32][33][34][35]. However, the measurement vectors of PL-TCAI conform to sub-Gaussian random distribution, and the stepsize selection criteria or functions based on the standard Gaussian random distribution are no longer applicable to PL-TCAI. Different PL-TCAI systems may require different stepsize selection criteria, and it is undoubtedly difficult to choose a new stepsize for PL-TCAI based on experience or repeated experiments, which is almost impossible to achieve in the practical PL-TCAI systems. Therefore, the optimal stepsize operation is necessary in the different PL-TCAI systems, which avoids the empirical selection of stepsize by calculating the stepsize of each iteration.
This paper is organized as follows. In Section 2, a new PL-TCAI configuration, which combines the incoherent detector array with the electronic coded aperture, is introduced in detail, and the imaging model of PL-TCAI based on the random frequency-hopping (R-FH) signal is derived. Section 3 analyzes the target reconstruction principle of PL-TCAI and proposes a sparse PR algorithm for PL-TCAI, in which the specific procedure of the SWFOS algorithm is given. In Section 4, numerical simulations are implemented to demonstrate the general performance of the SWFOS algorithm in the PR problem and the feasibility of the SWFOS algorithm to reconstruct targets in the proposed PL-TCAI. Finally, we summarize the main work and contributions of this paper in Section 5.

Proposed PL-TCAI Architecture
The configuration of our proposed PL-TCAI system is shown in Figure 1. It mainly includes a computer, a transmitter, a reflectarray coded aperture with an incoherent detector array. The computer is used to control different devices and to process echo data for imaging. During the transmitting process, the computer controls the transmitter to generate a specified form of terahertz signal and simultaneously loads the random modulation scheme to the coded aperture to modulate the terahertz wave. In the receiving process, the computer is used for computational imaging by processing the intensity data of the echo signal collected by the incoherent detector array. The transmitter is designed to emit terahertz continuous waves with a certain beamwidth of the main lobe, ensuring that it can sufficiently illuminate the entire coded aperture over a distance, while the side lobes need to be low enough to not affect the transmitted signal. Due to the different materials and modulation principles, the coded aperture usually has two modulation modes for terahertz waves to control their spatiotemporal independence, namely random phase modulation and random amplitude modulation. In this PL-TCAI system, a reflectarray coded aperture made of metamaterials [36][37][38] is selected for random phase modulation, which has proven to be practical in the terahertz band. The reflectarray coded aperture has the advantages of scalability, simple fabrication, and low cost [39]. It can achieve a wide range of phase shifts, and is also effective for large bandwidth signals. Here, we adopt a 1-bit digitally controlled reflectarray coded aperture [37], which is more convenient to load coding factors repeatedly and quickly. The reflectarray coded aperture is shown in the red block in Figure 1, which modulates terahertz waves into the spatiotemporal-independent radiation field in the imaging area. The incoherent detector array, which can only detect the intensity of the echo signal without phase information, is designed to be combined with the reflectarray coded aperture, as shown in the blue block in Figure 1. This design makes the PL-TCAI system more compact and integrated. Each incoherent detector element can independently transfer the intensity information of the echo signal to the computer for processing, which is inspired by multiple-output technology. The acquisition speed of the echo signal increases with the number of incoherent detector elements, but it should be noted that, in order to ensure the independence between different rows of the reference-signal matrix, the distance between the elements must be large enough. It is assumed that the target is located in an imaging area divided into W grid cells according to system resolution, and the strong scatterers of the target are exactly at the center of the grid cells. When the target appears in the imaging area, it changes the scattering coefficients corresponding to the grid cells. The purpose of the PL-TCAI system is to calculate the scattering coefficients of all grid cells in the imaging area so that the specific position, shape, and size of the target in the imaging area can be obtained from the distribution of the scattering coefficients. Obviously, the size and number of grid cells are the representation of imaging resolution and imaging area size, respectively. grid cells. When the target appears in the imaging area, it changes the scattering coefficients corresponding to the grid cells. The purpose of the PL-TCAI system is to calculate the scattering coefficients of all grid cells in the imaging area so that the specific position, shape, and size of the target in the imaging area can be obtained from the distribution of the scattering coefficients. Obviously, the size and number of grid cells are the representation of imaging resolution and imaging area size, respectively.  Figure 1. The architecture of the phaseless terahertz coded-aperture imaging (PL-TCAI) system based on combination of incoherent detector array and reflectarray coded aperture.

Imaging Model
According to Figure 1, the imaging model of the proposed PL-TCAI can be deduced in detail. In the transmitting process, the signal transmitted by the transmitter may be a single-frequency signal, a step-frequency signal, or a random frequency-hopping (R-FH) signal. Here, assuming the PL-TCAI system adopts a terahertz R-FH signal, it can be written as Equation (1): where () n st is the transmitting signal at time n t , A is the amplitude of the signal, j is the imaginary unit, c f is the carrier frequency, and n t f is the random frequency hopping at time n t .
The reflectarray coded aperture includes M phase modulation elements, and the position vector corresponding to the m-th element is m r . The transmitting signal ( , ) nm st r arriving at the m-th coded aperture element can be expressed as Equation (2): where 0 R is the position vector of transmitter, and c is the propagation speed of the wave. At the coded aperture, the transmitting signal is randomly phase modulated by the array elements. Controlled by the computer, the random phase modulation factor ( )  mn t corresponding to the mth element at time n t is loaded into the coded aperture. Then, the imaging area is divided into W grid cells, and the center position vector of the w-th grid cell is w r . The incoherent detector array has a total of Q elements, and the q-th element position vector is q r . Each detector element can detect the intensity of the signal and transmit it back to the computer for processing, according to the multiple-output technology. Hence, at time n t , for the echo signal detected by the q-th detector array element, the reference signal ( ) ( , ) q nw St r of the w-th imaging grid cell can be expressed as Equation (3):

Imaging Model
According to Figure 1, the imaging model of the proposed PL-TCAI can be deduced in detail. In the transmitting process, the signal transmitted by the transmitter may be a single-frequency signal, a step-frequency signal, or a random frequency-hopping (R-FH) signal. Here, assuming the PL-TCAI system adopts a terahertz R-FH signal, it can be written as Equation (1): where s(t n ) is the transmitting signal at time t n , A is the amplitude of the signal, j is the imaginary unit, f c is the carrier frequency, and f t n is the random frequency hopping at time t n . The reflectarray coded aperture includes M phase modulation elements, and the position vector corresponding to the m-th element is r m . The transmitting signal s(t n , r m ) arriving at the m-th coded aperture element can be expressed as Equation (2): where R 0 is the position vector of transmitter, and c is the propagation speed of the wave. At the coded aperture, the transmitting signal is randomly phase modulated by the array elements. Controlled by the computer, the random phase modulation factor ϕ m (t n ) corresponding to the m-th element at time t n is loaded into the coded aperture. Then, the imaging area is divided into W grid cells, and the center position vector of the w-th grid cell is r w . The incoherent detector array has a total of Q elements, and the q-th element position vector is r q . Each detector element can detect the intensity of the signal and transmit it back to the computer for processing, according to the multiple-output technology. Hence, at time t n , for the echo signal detected by the q-th detector array element, the reference signal S (q) (t n , r w ) of the w-th imaging grid cell can be expressed as Equation (3): Then, the intensity of the echo signal received by the q-th detector element at time t n can be expressed as Equation (4): where β w is the scattering coefficient corresponding to the w-th grid cell. According to Equation (6), we can derive the intensity of the echo signal received by the entire incoherent detector array at time t n , which can be written as Equation (5): where Here, [·] T and | · | denote the matrix transposition and the absolute value, respectively. Aside from that, the reference signal matrix S(t n ) at time t n is written as Equation (6): After the frequency of the transmitting signal is randomly hopped N times and the coded aperture is loaded with the coding scheme N times, the intensity detector array samples the echo signal N times and obtains L = Q × N samples. Finally, in combination with Equation (5), the imaging model of PL-TCAI is deduced as where Sr = [Sr(t 1 ), Sr(t 2 ), . . . , Sr(t N )] T and S = [S(t 1 ), S(t 2 ), . . . , S(t N )] T are the final received signal vector and the final reference signal matrix, respectively. Furthermore, ω = [ω 1 , ω 2 , . . . , ω L ] T is the additive measurement noise vector at the receiving terminal. As can be seen from Equations (6) and (7), for samples of the same size, sampling and coding with a single detector will be L = Q × N times, but when using a detector array with Q elements, there will be less sampling and coding, which will only take N times. Therefore, the proposed PL-TCAI with an intensity detector array can greatly reduce the number of sampling and coding, and directly improve the imaging efficiency.
Imaging processing is used to solve Equation (7) in order to accurately estimate β, which represents the target-scattering characteristics. Obviously, the imaging capability of PL-TCAI is affected by the structure of the reference signal matrix S. For PL-TCAI, an ideal S should be that all rows and columns are independent of each other. Fortunately, under the random phase modulation of the coded aperture, the columns of S can be independent of each other. In addition, by using the R-FH signal and increasing the distance between the detector elements, the independence of the rows of S can be increased.

Target Reconstruction Principle
Reconstructing the target from Equation (7), where only the intensity of the echo signal is measured, is a typical PR problem [26][27][28][29][30][31][32][33][34][35]. Classic algorithms for solving PR problems include the GS algorithm, the PhaseLift algorithm, the WF algorithm, etc. The GS algorithm is highly dependent on a priori information about the target, and often gets stuck into local minimums without converging during the solution process. The PhaseLift algorithm solves a PR problem by establishing a higher-dimensional linear equation, which converts a non-convex problem into a convex problem. Although the PhaseLift algorithm can provide an accurate estimate, as the dimension of β increases, the memory required for the calculation will far exceed the capabilities of ordinary computers because of matrix lifting. However, for a non-convex problem of PR, the WF algorithm directly adopts the gradient descent method to iteratively update the estimated β, searching for the global minimum. In order to successfully and accurately estimate β, the WF algorithm often needs excessive samples, where the number of samples L may satisfy L/W ≥ 6 [20,26].
As is well-known, obtaining a large number of samples is very unfriendly for a fast-imaging system. For the proposed PL-TCAI, although a detector array is adopted to speed up the detection of echoes, with the increase of the number of grid cells W, the number of samples L required by the WF algorithm becomes extremely large, thereby increasing the number of coding and sampling. In order to reduce the number of coding and sampling while increasing the speed of imaging, it is necessary to explore an algorithm that requires fewer samples. Recently, it was confirmed that utilizing the sparse prior information of the target can effectively reduce the samples required by the algorithm [33,34].
When the detected target is k sparse, it means that there are k strong scattering points in the imaging area. Therefore, for solving Equation (7), we can express the problem as the following: whereβ is an estimated scattering coefficient vector, and · 0 denotes the zero norm. S i is the vector of i-th row vector of S, and Sr i is the i-th element of Sr. Here, ω i is the i-th measurement noise, and according to the central limit theorem (CLT) [40], the noise ω i conforms to the Gaussian distribution, which can be written as ω i ∼ N (0, σ). Therefore, the probability density function of ω i can be written as the following: Here, the problem of solving the probability density function of ω i in Equation (9) can be transformed into the problem of solving the estimatedβ, which can be solved by the following maximum likelihood function: In fact, Equation (10) is a non-convex quadratic system, which has many local minimums, and driving it to converge to a local minimum will be very difficult. This problem is customarily called non-deterministic polynomial hard. However, as mentioned in the previous section, the rows and columns of the reference signal matrix S obtained in PL-TCAI are independent of each other. In particular, each row of S is affected by the random phase modulation, resulting in a random distribution of the row vector S i , i = 1, 2, . . . , L, which has a benign geometrical structure. Hence, it is quite possible to find a global minimum from Equation (10) to obtain the estimatedβ.

SWFOS Algorithm
To search for the global minimum of Equation (10), we propose a sparse phase retrieval algorithm with the optimal stepsize (SWFOS) for fast iterative update. Recovering the support of β is arranged in the first step of the algorithm, and then we initialize the estimatedβ 0 via a truncated spectral method under the recovered supports. Finally, the estimatedβ is iteratively updated using the gradient descent method, where each iteration adopts the optimal stepsize and adds a hard thresholding operation. The following is a detailed introduction to the SWFOS algorithm.

Recovering the Support
It is assumed that the scattering coefficient vector conforms to β ∈ R W and β 0 = k, which represents that there are only k strong scattering points in the imaging area with W grid cells. Using the methods mentioned in References [33,34], we assume that β has a support β sup that complies with β sup ∈ R W and β sup 0 = k. Here, the support β sup is an index set, which is used to constrain the estimatedβ during the initialization procedure. Besides, we define the random variables Here, S i is the i-th row vector of the reference signal matrix S, and we assume that S i ∼ N (0, I W ). Thus, based on the moment of Gaussian variables, it can be calculated that E S 4 i,j = 3 and E S 2 i,j = 1. Finally, according to the rotational invariance property of Gaussian distributions, the following equation can be obtained: where β / j is obtained by removing the j-th entry from β, while S / j is also obtained in the same way.
Note that E Y i,j = 2β 2 j + β 2 2 will be valid at j ∈ β sup and β j 0. Otherwise, if j β sup or β j = 0, there will be E Y i,j = β 2 2 . With the help of the strong law of large numbers, we know that, when L is large enough, the sample average will approach the ensemble one. With the increase of L, it is easy to get the following: It can be seen from Equations (11) and (12) that Y i, j will be affected by β j . For example, when β 2 j is very large, the corresponding Y i,j will also be very large. It should be noted that E Y i,j = 2β 2 j + β 2 2 will be valid at j ∈ β sup and β j 0. Therefore, we can sort out the largest k values from Y i,j can be equal to the support β sup with a probability of at least 1 − 6/L when L ≥ C 0 k 2 log(LW) for some constant C 0 and minimum nonzero entries β min := min j∈β sup β j on the order of 1 + √ k / β 2 .

Initialization via Truncated Spectral Method
Next, we need to initialize the estimatedβ 0 under the recovered support β sup 0 . It can be noticed from Reference [26] that the spectral method is a very attractive initialization procedure, which initializesβ 0 by calculating the leading eigenvector of Y SM = 1/L L i=1 Sr i S i T S i . If the initializedβ 0 is sufficiently accurate, the estimatedβ can be guaranteed to converge in the iterative update. However, this method is sensitive to the complexity of the sample, which needs to exceed W log W. The smaller the sample complexity, the lower the success rate of initialization and the larger the relative error of initialization. Fortunately, in response to this shortcoming, an improved algorithm called the truncated spectral method is proposed in Reference [27], where the initialization is completed by calculating the leading eigenvector of Moreover, since the support β sup 0 of β is recovered in the previous section, we need to constrain S i by deleting some elements that are not in the support β sup 0 . Thus, we can get a special matrix, which is written as the following: where α y is a truncation threshold.

Iteration via Gradient Descent Scheme
One of the most important steps in solving Equation (10) is to iteratively update the estimated β by using a gradient descent scheme, which ultimately searches for the global minimum in each iteration. Based on the Wirtinger derivative [26], the gradient of f β in Equation (10) can be defined as the following: Therefore, in the τ-th iteration, the estimatedβ can be updated once in the direction of the negative gradient, which is expressed as Equation (15): where µ τ is the stepsize of the τ-th iteration. When τ = 0,β 0 is obtained from the previous section.
In each iteration of Equation (14), a hard threshold operation can be performed on the calculatedβ τ+1 , which can be written as Equation (16) Here, T k βτ+1 represents an operation that keeps k entries of the largest absolute values ofβ τ+1 and sets the other entries ofβ τ+1 to zero. The hard threshold operation can reduce the freedom dimension and constrain the searching domain [33,34].

Optimizing the Iteration Stepsize
It is worth noting that the iterative update in Equation (15) is performed along the negative gradient direction, and each step of the iteration needs to set a stepsize µ τ . In References [26][27][28][29][30][31][32][33][34], the stepsize µ τ in each iteration is chosen by experience. In the early iterations, a small µ τ should be taken due to the large relative error of the estimatedβ, but, as the iterations count increases, a larger µ τ can be adopted for accelerating the convergence rate. Actually, if µ τ is too small, it will cause the algorithm to converge slowly, and if µ τ is too large, it will cause the algorithm to diverge. Thus, picking out a reasonable µ τ in each iteration is challenging but quite necessary for the algorithm. In order to optimize µ τ in each iteration, we can formulate the problem of selecting µ τ of the τ-th iteration as the following: Equation (17) is an optimization problem about the stepsize µ τ , and many efficient methods can be used to solve it. Here, we adopt a novel method detailed in Reference [35]. By combining Equations (10) and (17), it is easy to get the following equation: Obviously, Equation (18) is a univariate quartic. Hence, in order to obtain an optimal µ τ , the function f βτ − µ τ ∇ f βτ / 1/L L l=1 Sr l must satisfy the following optimization condition: After fully expanding Equation (19), it will be written as Equation (20): where Sr i τ = S iβ τ 2 , and Q = S i ∇ f βτ / 1/L L l=1 Sr l . Moreover, (·) * and Re(·) denote the complex conjugate and the real part of the complex number, respectively. It can be seen from Equation (20) that it is a univariate cubic equation about µ τ , which can be simply written as the following: where the coefficients can be calculated separately by According to the relevant mathematical knowledge, when the coefficients of the univariate cubic equation are all real values, the form of the root has two possible cases, and we will only adopt its minimum real value [35].
Finally, the procedure of the SWFOS algorithm is detailed in Algorithm 1.
Set the optimal stepsize µ τ by computing the root of the following:

Numerical Simulation
In this section, the numerical simulation results are presented to illustrate the performance of the SWFOS algorithm and the feasibility of the PR algorithm in PL-TCAI. First, the comparison results between the SWFOS algorithm and other algorithms are given, showing the advantages of SWFOS in some aspects. Second, the reason why the SWFOS algorithm can be adopted in the PL-TCAI system is given, which proves the practicability of SWFOS. Finally, the performance results of the SWFOS algorithm in the PL-TCAI system are given, which demonstrate that the SWFOS algorithm can effectively reduce the number of measurement samples and shorten imaging time at low SNRs. All simulations are performed on a computer with Intel Xeon CPU Gold 5188 at 2.3 GHz and 128 GB of memory. In the simulation, the ratio between the number of samples L and the number of grid cells W can be written as L/W. The imaging performance is evaluated by relative imaging error (RIE) which can be calculated as the following: whereβ is the estimated target scattering coefficient vector and · 2 denotes the Euclidean norm of a vector. Here, we adopt the empirical success rate of 100 Monte Carlo trials to evaluate the recovery performance of various algorithms, where a success of a trial is declared if the RIE is below -100 dB.
In order to estimate the correlation between different columns of S, we adopt the autocorrelation coefficients of columns, which is also called the space independence function in Reference [9] and is defined as the following: where ρ column is the autocorrelation coefficients of columns. S · j and S i· are the j-th column and the i-th row of the matrix S, respectively. Moreover, the operation cov(A, B) refers to find the covariance for two vectors A and B, and the operation σ(A) refers to find the standard deviation of vector A.

Performance of SWFOS Algorithm
To demonstrate the performance of the SWFOS algorithm, the SWF, WFOS, SPARTA, TAF, and RAF algorithms are implemented as comparisons. Among them, the sparse PR algorithms include the SWFOS, SWF, and SPARTA algorithms, and the others are non-sparse algorithms. In this simulation, considering only the complex case, we assume that all reference signal vectors {S i } 1≤i≤L satisfy the independent standard complex Gaussian distribution, which can be expressed as {S i } 1≤i≤L ∼ N (0, I W /2) + jN (0, I W /2). We set the target scattering coefficient vector β as a real Gaussian random vector, satisfying β ∼ N (0, I) and β ∈ R 1000 . For fairness, the initialization of all algorithms is obtained by 100 power iterations, and a fixed number of gradient iterations T = 1000 are run for all algorithms.
From Figure 2a, we can find that the SWFOS, SWF, and SPARTA algorithms can successfully recover β accurately at a low ratio L/W due to the consideration of sparsity priori, while conventional algorithms, such as TAF, RAF, and WFOS algorithms, require a high ratio L/W. This intuitively proves that when the target is sparse, the number of samples will be greatly reduced if the target imaging is reconstructed accurately by the sparse algorithm. In addition, we can also carefully observe that the recovery performance of the SWFOS algorithm is slightly better than the SWF and SPARTA algorithms in the complex case. As can be seen from the curves, the SWFOS algorithm has a stable 100% empirical success rate at L/W ≥ 1.1, while the SWF and SPARTA algorithms have a 100% empirical success rate at L/W ≥ 2.6 and L/W ≥ 1.6, respectively. Figure 2b is a result of the SWFOS, SWF, and SPARTA algorithms, respectively, recovering β of various sparsity k. From Figure 2b, we can see that in the complex case, as the sparsity k increases, the empirical success rate of all algorithms decreases significantly. However, we can still find that the SWFOS algorithm is more robust to sparsity k than the other two algorithms. The curve visually shows that the empirical success rate of the SWFOS algorithm fluctuates slightly and keeps in the high level when k ≤ 20 and fixed L/W = 1. In combination with Figures 2a,b, we can also conclude that the robustness of the SWFOS algorithm to the sparsity k can be improved by increasing the ratio L/W. In order to compare the convergence of these three sparse algorithms, we analyze the RIEs of each iteration of algorithms during operation. As can be seen from Figure 2c, the convergence rate of the SWFOS algorithm, which adopts the optimal stepsize, is significantly faster than that of the SWF and SPARTA algorithms. When k = 10 and L/W = 1 are fixed in the complex case, the RIE of the SWFOS algorithm decreases to -700 dB after 350 iterations, while that of the SWF algorithm and the SPARTA algorithm falls to the same level after 650 iterations and 500 iterations, respectively. This result demonstrates our analysis in the previous section that the SWFOS algorithm can determine the optimal stepsize and direction in each gradient descent iteration to quickly search for the global minimum.
All of Figure 2a-c is obtained in the absence of noise, while the effects of noise on these sparse algorithms are discussed in Figure 2d. We still fix k = 10 and L/W = 1, and add Gaussian white noise to the imaging model, which is described as Equation (7). The SNR varies from -10 to 50 dB, and all sparse algorithms perform T = 1000 iterations, respectively. From Figure 2d, we can see that all algorithms are sensitive to noise, and the RIEs of all algorithms increase gradually as the SNR decreases. However, even if the SNR is reduced to 10 dB, the three sparse algorithms can still recover β. In addition, it can be seen from the curves that the SWFOS algorithm performs slightly better than the SWF and SPARTA algorithms, especially in high SNR conditions.

Application of SWFOS Algorithm in PL-TCAI
We have mentioned above that the SWFOS algorithm performs well at solving the sparse PR problem when all reference signal vectors {S i· } 1≤i≤L satisfy the independent standard complex Gaussian distribution. However, considering the actual situation, it is difficult for {S i· } 1≤i≤L to meet this condition in the proposed PL-TCAI system. In order to prove the feasibility of the SWFOS algorithm applied to the PL-TCAI system, we analyze the characteristics of the reference signal matrix S derived from the PL-TCAI system. Here, we assume that imaging area in the PL-TCAI system is divided into W = 32 × 32 grid cells, and L = 1W samples are obtained by the incoherent detector array, which only receives the intensity information of the echo signal. In the simulations, a THz R-FH signal is used as the transmitting signal, and the basic parameter settings are given in Table 1. Table 1. Basic parameter settings used in the simulations.

Parameter Value
Center frequency 345 GHz Pulse width 100 µs  Figure 3a is a three-dimensional frequency histogram, in which can be the data distribution characteristics of each row of S can be visually observed. From Figure 3a, we see that the real values of each row of S is mainly between -50 and 50, and the distribution trend of the amount is gradually decreasing from 0 to both sides. For further analysis, the probability distribution function (PDF) of each row of S is plotted in Figure 3b. As shown by the curves of the different colors in Figure 3b, we can conclude that the real values of each row of S conform to the sub-Gaussian distribution, and the same results also appear in the imaginary values of each row of S. The reason for these results is that we adopt the coded aperture in the PL-TCAI system to randomly modulate the phase of the transmitting signal, making the reference signal corresponding to the imaging area random. According to the central limit theorem (CLT) [40], when directly affected by the random phase modulation, each S i· becomes a random vector and the elements in S i· approximately obey the Gaussian distribution. In fact, the vector S i· is the mathematical expression of the instantaneous radiation field of the reference signal. In order to vividly display the modulated reference signal, we plot the radiation field of the reference signal at a certain moment in Figure 3c. From Figure 3c, we can see that the distribution of the instantaneous radiation field is obviously random, which is consistent with our previous analysis. Furthermore, Figure 3d is the result of the space independence function of the PL-TCAI system, which is calculated from Equation (23). As can be seen from Figure 3d, the columns of S are random and independent, which corresponds to the result that the elements in S i· conform to the random sub-Gaussian distribution. It is worth mentioning that Figure 3c,d can also describe the imaging ability of the PL-TCAI system, where the more random the radiation field is and the better the space independence function is, the stronger the imaging ability is. Next, to verify the validity of the proposed algorithm in the PL-TCAI system, we adopt the SWFOS, SWF, and SPARTA algorithms to solve the imaging equation of the PL-TCAI system, where S is derived from the PL-TCAI system but does not conform to the independent standard complex Gaussian distribution. As in Section 4.1, the initialization of all algorithms is obtained by 100 power iterations, and a fixed number of gradient iterations T = 1000 are run for all algorithms. Since the priori sparsity k of the target is unknown in actual radar imaging, in all sparse algorithms, we set k ≈ √ W = 32 as an estimate to ensure a high empirical success rate. W is the number of grid cells in the imaging area and is also the dimension of the scattering coefficient vector β. Figure 4a is the original sparse target, and Figure 4b-d shows the imaging results for the sparse target when L/W = 1.5 in the PL-TCAI system with the setting SNR = 20 dB. From Figure 4b-d, we can see that only the SWFOS algorithm accurately reconstructs the target, and the other two algorithms do not converge and cannot reconstruct the target. In order to analyze the performance of three sparse algorithms in the PL-TCAI system, we calculate the empirical success rates with L/W increasing. As can be seen from Figure 5, only the SWFOS algorithm can work properly, and the other two algorithms cannot accurately solve the imaging equation. The SWFOS algorithm has a stable 100% empirical success rate at L/W ≥ 1.7, while the SWF and SPARTA algorithms have an 0% empirical success rate for all these ratios. Compared with the ideal S, the performance of the SWFOS algorithm in PL-TCAI is degraded because the vector S i· conforms to the random sub-Gaussian distribution. Because the stepsize of the SWF and SPARTA algorithms is chosen empirically, for an actual PL-TCAI system, when S i· does not satisfy the independent standard complex Gaussian distribution, it is necessary to reselect a new stepsize. However, for the SWF and SPARTA algorithms, it is extremely difficult to select the appropriate stepsize for the PL-TCAI system based solely on experience. Such a dilemma will never appear in the SWFOS algorithm, which can calculate an optimal stepsize in each iteration to accurately solve the imaging equation.  Furthermore, in order to verify that the proposed SWFOS algorithm can significantly reduce the number of samples needed for sparse target imaging in the PL-TCAI system, we compare the performance of the SWFOS and WFOS algorithms with different numbers of samples when setting SNR = 20 dB. Figure 6 shows the empirical success rates of SWFOS and WFOS in different L/W ratios. From Figure 6, we can see that, when reconstructing the sparse target, the SWFOS algorithm has a high empirical success rate at L/W = 1.5, while, for the WFOS algorithm, without using sparse prior of targets, its empirical success rate is still very low at L/W = 6, which is only about 30%. Figure 7 shows the imaging results of the SWFOS algorithm and the WFOS algorithm for the sparse target in the different number of samples, and we can see that the SWFOS algorithm can accurately reconstruct the target even if L/W = 1.5. Under the same conditions, the WFOS algorithm will fail to reconstruct the target obviously. Of course, with the increase of L/W, the WFOS algorithm can gradually reconstruct the target. As shown in Figure 7b, most strong scatterers of the target can be reconstructed when L/W = 6, but there are still many weak spurious scatterers. To further evaluate the imaging performance of the proposed algorithm in the PL-TCAI system, we record the RIEs and runtime of two algorithms, and the result is shown in Table 2. From Table 2, we can see that, although the SWFOS algorithm which utilizes the sparse prior needs fewer samples to reconstruct the sparse target, it consumes less time and makes the RIE smaller.   Finally, Figure 8a-c illustrates the imaging results of the SWFOS algorithm for the sparse target when L/W = 1.5 in the PL-TCAI system with different SNRs. From Figure 8, we can find that the sparse target can be fully identified at SNR = 15 dB, while, when SNR = 10 dB, the sparse target can only be partially reconstructed. To analyze the influence of noise on the SWFOS algorithm in the PL-TCAI system, we calculate the RIEs with the SNR increasing when setting L/W = 1.5 and k = 10, and the result is plotted in Figure 9. From Figure 9, we can see that the SWFOS algorithm can perform well in the PL-TCAI system with low SNRs, such as RIE = -70dB when SNR = 20 dB. With the improvement of SNR condition, the imaging quality of SWFOS algorithm is greatly improved, and it can reconstruct the sparse target accurately.

Discussion
Numerical simulation results show that the SWFOS algorithm can solve the Phase Retrieval problem well and can also be effectively applied in PL-TCAI. However, the derivation method of the reference signal S (q) (t n , r w ) in the PL-TCAI proposed in this paper is based on the distance and time delay, and its accuracy is not completely applicable to the far field. Since the accuracy of reference signal S (q) (t n , r w ) has a great influence on imaging results of PL-TCAI, ensuring the accuracy of reference signal S (q) (t n , r w ) will be the main work we are facing in the actual PL-TCAI system. In addition, since the atmospheric water vapor has a high absorption for terahertz waves, it is difficult for PL-TCAI to achieve high-resolution imaging at a long distance in the atmosphere, which is also a huge challenge for all terahertz imaging.

Conclusions
In this paper, we propose a compact PL-TCAI system configuration based on the combination of the incoherent detector array and the coded aperture, which can achieve fast imaging, high-resolution imaging, and forward-looking imaging without any relative motion. We also propose a sparse phase retrieval algorithm named SWFOS for the sparse target in PL-TCAI, which aims to reduce the number of measurement samples required for imaging. Firstly, the proposed configuration of PL-TCAI is introduced in detail, in which the multiple-output technology adopted by the incoherent detector array can significantly reduce the number of coding and sampling compared with the single detector. Then, according to the proposed system configuration, the imaging equation of PL-TCAI which adopts the R-FH signal as the transmitting signal is derived in detail. Besides, we analyze the target reconstruction principle of PL-TCAI and propose the SWFOS algorithm for sparse target imaging. The specific procedures of the SWFOS algorithm include the support recovery, initialization by truncated spectral method, iteration via gradient descent scheme, hard threshold operation, and iteration stepsize optimization. Finally, numerical simulation results show the basic performance of the SWFOS algorithm and demonstrate the feasibility of applying the SWFOS algorithm to sparse target imaging in PL-TCAI. In conclusion, the proposed structure and algorithm can significantly reduce the amount of coding and sampling required for sparse target imaging, shorten the imaging time of the PL-TCAI system, and provide a new method for high-speed PTCAI, at nearly real-time frame rates.