Online Sparse DOA Estimation Based on Sub–Aperture Recursive LASSO for TDM–MIMO Radar

: The least absolute shrinkage and selection operator (LASSO) algorithm is a promising method for sparse source location in time–division multiplexing (TDM) multiple–input, multiple– output (MIMO) radar systems, with notable performance gains in regard to resolution enhance-ment and side lobe suppression. However, the current batch LASSO algorithm suffers from high– computational complexity when dealing with massive TDM–MIMO observations, due to high– dimensional matrix operations and the large number of iterations. In this paper, an online LASSO method is proposed for efﬁcient direction–of–arrival (DOA) estimation of the TDM–MIMO radar based on the receiving features of the sub–aperture data blocks. This method recursively reﬁnes the location parameters for each receive (RX) block observation that becomes available sequentially in time. Compared with the conventional batch LASSO method, the proposed online DOA method makes full use of the TDM–MIMO reception time to improve the real–time performance. Addition-ally, it allows for much less iterations, avoiding high–dimensional matrix operations, allowing the computational complexity to be reduced from O (cid:0) K 3 (cid:1) to O (cid:0) K 2 (cid:1) . Simulated and real–data results demonstrate the superiority and effectiveness of the proposed method.


Introduction
In recent years, colocated multiple-input, multiple-output (MIMO) technology has received increasing amounts of attention and has been extensively used for assisted driving, remote sensing [1], geological exploration [2], and free space optics communication (FSOC) [3]. In radar remote sensing applications, the MIMO technology improves the angular resolution while reducing the physical channel count and antenna aperture, compared with conventional array antennas [4][5][6]. In FSOC applications, MIMO technology is also utilized to overcome attenuation and turbulence effects in atmospheric channels, thereby obtaining channel gain and improving the reliability of FSOC systems [7,8]. In terms of hardware complexity, the frequency-modulated continuous wave (FMCW) sequence of time division multiplexing (TDM) is the simplest option for colocated MIMO. Therefore, more practical applications and mechanistic studies are achieved in this way [9][10][11].
Targeted direction-of-arrival (DOA) estimation is an important research hotspot in MIMO radar signal processing [12]. The most fundamental method is the classical delayand-sum (DAS) approach, which recovers a target source by weighting and delaying the echoes. However, the DAS approach suffers from a low angular resolution and a high side lobe. Moreover, the resolution of the DAS approach can only be improved by increasing the number of channels, which is very uneconomical. Many methods have been proposed to effectively improve the DOA resolution of MIMO radar [13][14][15][16][17][18][19][20]. Based on the rotational invariance property of spatially correlated matrix signal subspaces, the estimation of signal parameters via the invariance technique (ESPRIT) algorithm was proposed to improve the DOA estimation accuracy of MIMO radar [13,14]. The parallel factor analysis (PARAFAC) algorithm was proposed to improve the multitarget resolution [15]. A dimensionalityreducing Capon method has been proposed for the DOA estimation of monostatic MIMO radars [17], and this method reduces the dimensions of the MIMO radar sub-aperture through a dimensionality-reducing transformation, thereby accelerating the Capon algorithm. In [1,21], the iterative adaptive approach (IAA) was introduced for MIMO radar imaging, and it was shown that the IAA can work with few or even a single snapshot, and can provide a higher azimuthal signal source localization accuracy. However, this method suffers from a high computational complexity due to its high dimensional matrix inverse computation per iteration. The sparse representation (SR) has been used in an estimation method that contains DOA estimation for colocated TDM-MIMO radars [19]. The least absolute shrinkage and selection operator (LASSO) method [22] is an efficient tool in the fields of application of sparse linear regression [23,24]. This method consists of a linear least square optimization regularized by an additional penalty term of the L 1 -norm as a measure of sparseness. The LASSO method is demonstrated to improve the DOA estimation performance of the array signal, with higher resolution and a lower side lobe [25]. Due to its superiority, many studies have been carried out on LASSO, including computational efficiency and theoretical analysis, etc. [26][27][28].
In spite of the performance gain, the current batch LASSO algorithm can only handle the complete aperture data already acquired by the MIMO radar. Since the LASSO algorithm requires operations, such as matrix inversion and matrix multiplication and a large number of iterations; the computational complexity of batch implementation is notably high, especially for massive TDM-MIMO data sets. There are several advances in algorithmic acceleration that might be combined with the LASSO algorithm to mitigate the computational burden. Some fast matrix inversion methods, such as the Newman series [29] and Jacobi method [30,31], were proposed to handle the problem of high computational complexity. In [32][33][34], the Neumann series (NS) is considered to perform matrix inversion approximation (MIA). The main contribution of the Neumann method is to convert the matrix inverse calculation into matrix multiplication, which is more suitable for hardware platforms, but the complexity is equal to or higher than that of the exact calculation method, such as the QR-based method [33]. The Jacobi method reduces the complexity from O K 3 to O LK 2 , where L is the number of iterations. However, the Jacobi method converges slowly and, thus, implies higher latency [35]. Although the above fast matrix inversion methods can reduce the computational complexity of the matrix inversion, the approximation error is introduced, which should be carefully controlled. Furthermore, they only reduce the complexity of the matrix inversion, whereas the burden caused by high-dimensional matrix multiplication and a large number of iterations in the batch LASSO algorithm cannot be well handled.
To solve the high-complexity problem of the LASSO method, this paper presents a subaperture recursive LASSO method to allow for online DOA estimation. Online processing refers to the way of processing comparing with the traditional batch processing [28,36,37]. In the way of batch processing, the processor must wait for collecting the measurement completely, and process the measurement taking additional time cost. In contrast, online processing processes the stream of MIMO sub-aperture measurements accordingly; that is, updating the result when a new sub-aperture measurement is available. The benefit of online processing is that it can save the time cost of measurement collection. Meanwhile, online processing also allows for much less iterations, avoiding high-dimensional matrix inversion and multiplication. From the following analysis, it can be seen that the computational complexity can be reduced from O K 3 to O K 2 . This paper is structured as follows. In Section 2, the MIMO radar signal model is reviewed and the TDM-MIMO sub-aperture data reception scheme is derived. Then, the proposed online TDM-MIMO radar DOA estimation method is derived in detail. In Section 3, simulated and real-data results demonstrate the effectiveness of the proposed method. In Section 4, the quantitative comparison and future work are discussed. Section 5 presents the conclusions.

Materials and Methods
A colocated MIMO radar that contains isotropic uniform linear receiving and transmitting arrays (ULA) is considered. The receiving array element spacing is d r = λ/2, where λ denotes the carrier wavelength. The location distributions of M r receiving antennas and M t transmitting antennas are shown in Figure 1 (with four transmitters and four receivers as an example). The transmitting array element spacing is d t = M t · d r . According to the MIMO theory, the equivalent virtual array of this MIMO can be considered a single-input, multiple-output (SIMO) antenna with 16 elements.

Signal Model
Suppose that there are K impinging signals with angular parameter, θ k , k = 1, . . . , K. According to Figure 1, the steering expression of the m-th virtual array element is written as where m t = 1, 2, . . . , M t and m r = 1, 2, . . . , M r represent the sequence number of transmitting and receiving array elements, respectively. In addition, m, m t and m r satisfy the relationship m = (m t − 1)M r + m r . For all virtual elements, the steering vector can be denoted as where M = M t · M r . We assume that there are L TDM periods and each period contains N pulses. Hence, the -th complex baseband receiver signal y( ) ∈ C M can be denoted as [19] where represents the Hadamard product, s( ) the unknown and deterministic target signals, e( ) the additive zero mean Gaussian noise, and σ 2 the variance. ω d denotes the Doppler frequency of the target. γ = t ⊗ 1 M r represents the TDM-induced Doppler frequency shift (DFS), where t = [t 1 , t 2 , t 3 , t 4 ] represents the time-division scheme of the signal transmission, and 1 M r denotes the full one-column vector with a length of M r . To simplify the model, it is assumed that ω d = 0, and N = 1. For the -th snapshot of a given TDM period, model (3) can be simplified to where the steering matrix is given by The receiver signals y of the M channels are arranged as 1, 2, · · · , M t according to the sub-aperture order.
The TDM-MIMO radar system transmits the signals in a time-division manner per the transmitting scheme, as shown by Figure 2a. Only one transmitting antenna is available for each time slot. In the case of the four receivers, the received signals are shown in Figure 2b. The blocks with order 1-4 denote the sub-aperture data received in the first period, and so on. Any set of sub-apertures is equivalent to a SIMO subarray with M r receivers. In the traditional batch processing method, all sub-aperture data must be received completely to formulate the full-aperture data before processing. When dealing with a massive TDM-MIMO data set, the dimension of aperture data will become quite high. As a result, the high dimensionality of the full-channel data increases the temporal and spatial complexity. Moreover, the period of reception will introduce additional delay. Both of these reasons lead to poor real-time performance of batch processing. To solve this problem, we propose an online processing method that processes sub-aperture data along receiving. The proposed method allows beamforming to begin directly after the first sub-aperture data are received. Doing so introduces two benefits, one is that it can save the time cost of measurement collection, the other is that online processing allows for much less iterations, avoiding high-dimensional matrix inversion and multiplication, resulting in lower temporal and spatial complexity.

Proposed Method
Based on model (4), the LASSO method attempts to recover s by solving the following convex optimization problem [22] arg min where λ is a parameter that weighs the sparsity of s. Equation (6) can be solved directly by convex optimization, such as the CVX tool box in MATLAB. However, the source code is not open, which cannot be ported to a hardware platform. In recent years, many iterative methods, such as split Bregman [38], IRLS [39], etc., have been proposed. For example, the split Bregman solution can be written as [38] where d and b denote the introduced temporary vectors, µ an additional regularization parameter, q the iteration number, I the identity matrix. The function shrink(x, η) = sign(x) max(|x| − η, 0), and sign(·) denotes the sign function. It is seen that matrix inversion, matrix multiplication, and a considerable number of iterations are required, so the computational complexity is quite high, which is not suitable for engineering applications. To solve this problem, starting from the cyclic minimization [28,40,41], we propose a subaperture recursive online processing method.

Cyclic Minimization
By the minimization of Equation (6), only one component s i in s is considered at a time. Letỹ i y − ∑ k =i a k s k . Then, the cost function can be written as [28] To solve the scalar minimization problem, the i-th variable to be solved is transformed into the polar form s i = r i e jθ i , where r i 0 and θ i ∈ [−π, π). Then, the quadratic term in Equation (6) can be rewritten as By substituting (9) into (8), a cost function with r i and θ i denoting independent variables is obtained According to (10), the minimizedθ i can be expressed aŝ Then, let where * denotes the conjugate transpose of a matrix. Hence, Equation (10) can be rewritten as The first-order and second-order derivatives of Equation (13) are Therefore, (13) is convex, and the minimization of r should be achieved at dJ/dr = 0 and denoted asr Ultimately, the element-level minimization can be expressed as [28] According to Equations (16) and (11), when updating each elementŝ i , the convex optimization cost function (8) decreases monotonically. Therefore, the resultŝ of the target source distribution can be achieved by traversing i = 1, 2, · · · , K for Equation (17).

Proposed Online Strategy
For TDM-MIMO radars, the sub-aperture data arrive sequentially over time. The p-th sub-aperture data and the corresponding steering matrix can be expressed as Vector y p ∈ C M r contains data for the p-th sub-aperture, and the length is M r , i.e., data blocks 1-4 in Figure 2b (if p = 1). Similarly, the steering matrix A p ∈ C M r ×K contains the steering vectors of the corresponding subarray elements. The optimization problem (6) can then be rewritten as arg min Let s ∈ R K denote the estimated result from the p − 1-th sub-aperture data. Based on the cyclic minimization, we first introduce an intermediate variable where s i denotes the current estimate of the i-th direction, and a p,i ∈ C M r the steering vector of the i-th direction for the p-th sub-aperture, namely, the i-th column of A p . The intermediate variable and recursive calculation variables are introduced: In minimizing r i , only the calculation of β i and γ i is needed. Considering Equation (22), the expression (20) can be simplified to where ζ i denotes the i-th element of ζ p , and Γ p ii the i-th diagonal element of matrix Γ p . Similarly, Equation (11) can also be represented by recursive variables aŝ Hence, the i-th element of the DOA is obtained: By traversing i = 1, 2, · · · , K for Equations (24)-(26), the target source distribution resultŝ obtained by solving the p-th sub-aperture data are obtained. Variables z p can be updated through z p = z p + a p,i s i −ŝ i . Then, according to (21), the update of ζ p can be written as where Γ p i denotes the i-th column of Γ p . For each new sub-aperture data block, ζ p is initialized to ζ p = ρ p − Γ p s. For the first sub-aperture (p = 1), the estimate is initialized as s = 0. Furthermore, the cyclic calculation (24)-(27) requires Q iterations. According to experience, good results are usually obtained after 8 to 10 iterations.
With the arrival of the subsequent MIMO sub-aperture data, the DOA resultŝ can be obtained by reiterating (21)-(27) using Equation (21) (which is achieved through the previous result s by the previous p − 1 sub-aperture data) and the new sub-aperture data y p . The online processing flowchart for the sub-aperture update is shown in Figure 3.
Calculate XX by (22) Figure 3. Flowchart of the online sub-aperture update.

Computational Complexity Analysis
The computational complexity of the DAS method, the IAA method [1,21,42], the traditional LASSO method, and the proposed method are compared, as shown in Table 1.

Method Calculation Times of Multiplication and Division
Computational Complexity Table 1, it is seen that the proposed methods reduce the computational complexity by an order of magnitude compared with the IAA and LASSO method. In addition, the complexity of the proposed method can reach the order of magnitude of the RLS method (O K 2 ). Table 2 shows the computational complexity for the methods with K = 4M r M t and varying M r × M t . Typically, when M r × M t = 1000, the proposed method reduces the required computational time by 194 times and 37 times, compared with the LASSO and IAA, respectively.

Results
In the simulation and experimental verification, both the target and the platform are considered stationary. All of our simulations and experiments were conducted using a 64-bit MATLAB R2018a on a PC workstation with an Intel Core i5-9500 CPU, 3.0 GHz and 16 GB RAM. The proposed method is compared with the conventional DAS, IAA, and LASSO methods in terms of computational time and DOA accuracy. Because this study only concerns the angular estimation performance, the targets were set at the same distance. The root mean square error (RMSE) of the angular estimates is used for evaluating the performance of simulation results and is defined as where θ p denotes the true target value of the p-th target grid, andθ p represents the estimate of the p-th target grid. The signal-to-noise ratio (SNR) is defined as where P s stands for the signal power, and σ 2 is the variance in the additive Gaussian white noise.

Simulation Results
The azimuth of two unrelated sources with equal power is assumed to be [−5 • , 0 • ]. The scan range is between −90 • and 90 • , and the number of points is K=512. The SNR of the signal is SNR = 10 dB. A MIMO radar system with two transmitters and four receivers (carrier frequency f c = 77 GHz; receiving array element spacing: d r = λ/2 = 1.9 mm; transmitting array element spacing: d t = M t · d r = 7.6 mm) is investigated. Figure 4a illustrates the results of 10 Monte Carlo experiments with the traditional DAS method. Although the average processing time is only 0.0002 s, effectively distinguishing between the two adjacent targets is impossible. Figure 4b shows the results of the 10 Monte Carlo experiments with the IAA method. The average processing time is 0.0879 s. The results of the 10 IAA experiments show lower side lobes, and the adjacent targets are basically distinguished and located, which is better than that shown in Figure 4a.
However, the IAA method also has drawbacks, such as an unsatisfactory resolution, a high side lobe, and high computational complexity. Figure 4c shows the results of the 10 LASSO experiments. It is seen that they can well distinguish the adjacent target directions, and the angular estimation error is highly non-significant. Nonetheless, the LASSO method has a high algorithmic complexity and requires the support of complete aperture data, thus exhibiting poor real-time performance. Furthermore, its average processing time is 0.2403 s. As shown in Figure 4d, the 10 results of the proposed method are nearly consistent with the LASSO results. Compared with the results of the traditional DAS and IAA methods, the results of the proposed method provide better resolving of adjacent targets and side lobe suppression.  The average processing time for each sub-aperture data update takes only 0.0215 s. Furthermore, Figure 5 illustrates the DOA result at the arrival of each sub-aperture data block under simulation conditions by the proposed method. The real-time DOA results are progressively improved as the MIMO sub-aperture data are received, and the best performance is achieved at the receipt of the last data block (the DOA performance of the LASSO method). Figure 6 shows the RMSE and CRB of angle estimation by the four methods for SNRs ranging from −5 to 30 dB [44].
Moreover, we can reduce the estimation error of DOA by compensation in two ways. One is the compensation for mutual coupling of array elements [45,46]. The other is the motion compensation: the Doppler phase needs to be considered for the moving platform [47,48]. Good research has been carried out on the above two compensation ways. Because this article mainly focuses on real-time signal processing algorithms, it will not be further discussed here.

One-Dimensional Point Target Experiment
In this section, two adjacent corner reflectors and a MIMO radar with two transmitters and four receivers are used to verify the DOA performance of the proposed method. The layout of the radar system and corner reflectors is shown in Figure 7a; the corner reflectors are illustrated in Figure 7b. The radar and corner reflectors are at the same height, and the carrier frequency of the former is 77 GHz.  Figure 8a shows the DOA result of the DAS method for which the processing time is 0.0024 s. Due to the short distance between the two corner reflectors, their azimuths are indistinguishable. Figure 8b shows the DOA result of the IAA method for which the processing time is 0.1304 s. Figure 8a,b are consistent with the corresponding simulation results. This method effectively improves the DOA in resolution. The corner reflector is clearly separable, but the side lobes are still not well suppressed. Figure 8c shows the DOA result of the LASSO method for which the processing time is 0.2139 s. By taking advantage of the sparseness of the target sources, the LASSO method distinguishes the targets very well, and the side lobe is also quite well suppressed. Compared with traditional methods, the resolution is greatly improved. However, the LASSO method is not applicable to the real-time processing of TDM-MIMO systems. Figure 8d shows the DOA result of the proposed method. This method maintains the same angular resolution as the LASSO method and can update the DOA results of the targets in real time as the MIMO sub-aperture data are received in real time. In addition, the proposed method offers excellent online processing and a super-resolution capability. The average processing time for each sub-aperture data update is 0.0318 s. According to the peaks, the angular estimation errors of the LASSO method and the proposed method are 0.15 • and 0.17 • , respectively, after comparing with the calibrated average angle of the targets.

Two-Dimensional Surface Target Experiment
In this section, MIMO radar data are used to verify the performance of the proposed method. The original optical scene photographed by the drone is shown in Figure 9. The experimental parameters are listed in Table 3. The carrier frequency, the number of array elements, and the beam width have a determined relationship, which mainly affects the azimuth resolution. The pulse width and pulse repetition interval affect the signal to noise ratio, and bandwidth affects the range resolution.   Figure 10a shows the 2D results of the DAS method for which the processing time is 0.0973 s. It is seen that a row of vehicles in the parking lot is vaguely distinguished. However, due to the low resolution and poor side lobe suppression ability of the DAS method, the imaging results with regard to the vehicle in the center of the scene suffer from high side lobes, and the vehicles on the right side of the scene are blurred. Figure 10b shows the 2D results of the IAA method [1] for which the processing time is 124.8980 s. It is seen that the resolution of this image is effectively improved. The side lobes of the strong scattering point are well suppressed. Figure 10c shows the 2D result of the LASSO method for which the processing time is 294.1526 s. It is seen that the resolution of this image is further improved. The side lobes of the car in the center of the scene are significantly suppressed, and the car contour becomes much clearer. The number and position of vehicles on the right side of the image can be easily determined. Unfortunately, both the IAA and the LASSO methods have enormous computational complexity, and cannot implement the online processing of the MIMO data. Figure 10d shows the 2D results of the proposed method for which the processing time is 4.8980 s. It is seen that this image is much better than that shown in Figure 10a,b. The proposed method maintains the same angular resolution as the LASSO method and can update the imaging results as the MIMO sub-aperture data are received in real time. The image entropy (IE) [49] is employed to quantitatively evaluate the performance of the 2D surface target experiment, where low values of IE commonly indicate that the image is well recovered. The image entropy of the above processed results is shown in Table 4. It can be seen that the entropy of LASSO and the proposed method are smaller than that of the DAS and IAA methods. This proves that the results of LASSO and the proposed method are clearer than those of other methods. More importantly, compared with Figure 8c,d, it can be seen that the proposed method has almost no performance loss compared with the traditional LASSO.

Results Analysis
In this paper, a LASSO-based sparse DOA estimation method for online processing of TDM-MIMO radar sub-aperture data is proposed. In this paper, a detailed comparison of the computational complexity, simulation results, measured data of point targets, and measured data of surface targets by various methods was carried out.
The results obtained in this paper show that the proposed online DOA method works well under both simulated and experimental conditions. It should be noted that the computational complexity of the IAA method in Table 1 does not consider the number of iterations, and its total computational complexity should be expressed as O Q(M r M t ) 3 , where Q represents the number of iterations of IAA. Referring to the article [21], the iteration number of the IAA method is set to 15 in all simulations and measurements in this paper.
Since the proposed method can output DOA results from sub-aperture data in real time, we only give the complexity of each sub-aperture recursion in Table 1. The total calculation times of multiplication and division is expressed as (Q + 1)M t K 2 + 9Q + M r 2 + M r M t K, and the computational complexity is O (Q + 1)M t K 2 , where Q represents the number of iterations of the proposed method. M t represents the number of transmitting array elements, namely, the number of recursions required for full aperture. M r represents the number of receiving array elements, namely, the sub-aperture length. In addition, according to the configuration of the transmitting and receiving array elements, M t and M r can also interchange roles. Comparing (c) and (d) of Figures 4, 8 and 10, it is obvious that the proposed method can maintain the recovery performance of the traditional LASSO method while performing lowcomplexity online DOA estimation. It should be mentioned that the selection of parameter λ in Equation (7) refers to the article [28]. In the proposed method, λ = m t M r log 2 K, where m t = 1, 2, . . . , M t is the current recursive sub-aperture ordinal.
In the follow-up research, on the one hand, under the conditions of the moving platform and moving target, Doppler phase and model error compensation need to be considered to further improve the DOA accuracy of the proposed method. On the other hand, after the motion compensation, the relevant hardware platform should be built to implement the TDM-MIMO data's real-time processing on FPGA or DSP using the proposed method.

Extension to Optical Communication
In free-space optical communication, high-speed transfer of effective information can be achieved by quantum cascade laser (QCL) [50,51]. The article [52], discusses the possibility of realizing high-power broadband QCL. The article [53], demonstrates that optical radiation with a wavelength of about 10 µm in limited visibility is characterized by better transmission properties than near-infrared waves. The article [54] analyzes the performance of optical communication channels in the marine environment. The article [55] analyzes the performance of free space optics communication under atmospheric turbulence. This paper mainly discusses the TDM-MIMO radar in the microwave band. Moreover, in the field of FSOC, MIMO technology can be used to solve the problem in attenuation and atmospheric turbulence, which contributes to the atmospheric channel gain and the reliability of the FSOC system. In the article [8,56], a maximum likelihood parameter estimation method is used to improve the bit error rate (BER) performance of the MIMO FSOC system. Due to the large number of massive FSOC-MIMO antennas and the complex atmospheric turbulent propagation environment [3,57], the computational complexity of matrix inversion in the channel estimation is quite high. Therefore, traditional channel estimation methods cannot be directly applied in massive MIMO. A low-complexity massive MIMO channel estimation method must be employed. The proposed online subaperture recursive LASSO estimation method may also be employed to lower the estimation error and improve the running speed. This paper mainly discusses radar applications and the above articles belong to the field of communications. We should note that although the mathematical principle is the same, when this method is applied to MIMO optical communication, minor changes are required, and we will not go into much detail here.

Conclusions
In this paper, an sub-aperture recursive LASSO method is proposed for the online DOA estimation of TDM-MIMO radar using sub-aperture data blocks. The proposed method has two advantages. First, it makes full use of the TDM-MIMO reception time to improve the real-time performance of the DOA algorithm (i.e., online DOA estimation), and the complexity and memory usage per update are remarkably reduced when compared with traditional methods. Second, the proposed method is superior to the traditional DAS and IAA methods with respect to the DOA resolution because it can reach the DOA estimation accuracy of traditional LASSO methods. Simulations and experimental data verify the effectiveness of the proposed method. The online DOA estimation of actual TDM-MIMO radar will be realized through hardware programming in the future.  Data Availability Statement: This study did not report any data.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: