# Compressed Sensing-Based Range-Doppler Processing Method for Passive Radar

Authors: Xia Bai; Hejing Guo; Juan Zhao; Tao Shan
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

Passive radar (PR) systems use the existing transmitters of opportunity in the environment to perform tasks such as detection, tracking, and imaging. The classical cross-correlation based methods to obtain the range-Doppler map have the problems of high sidelobe and limited resolution due to the influence of signal bandwidth. In this paper, we propose a novel range-Doppler processing method based on compressed sensing (CS), which performs sparse reconstruction in range and Doppler dimensions to achieve high resolution and reduces sidelobe without excessive computational burden. Results from numerical simulations and experimental measurements recorded with the Chinese standard digital television terrestrial broadcasting (DTTB) based PR show that the proposed method successfully handles the range-Doppler map formatting problem for PR and outperforms the existing CS-based PR processing methods.


---

## Body

## 1. Introduction

Passive radar [1] (PR) is a kind of radar system which uses the existing transmitters of opportunity (such as FM [2, 3], GSM [4], and DVB-T [5, 6]) in space to achieve target detection and tracking and other tasks without special deployment or installation of transmitters. In recent years, PR has been widely concerned in the military and commercial fields because of advantages in terms of low-cost implementation, confidentiality, strong antijamming, and reduced electromagnetic pollution to the environment.In this study, we mainly focus on the generation of a passive radar range-Doppler (RD) map. Based on the matched filtering theory, the classical method uses the cross-ambiguity function (CAF) and fast Fourier transform (FFT) to calculate RD response. However, the classical method faces some challenges. On the one hand, the generated RD map has a high sidelobe level, which may get false target position and Doppler frequency. On the other hand, the signal bandwidth of passive radar is usually narrow compared to that of active radar. It limits range resolution, and this limitation leads to undesirable performance in various applications [7]. In order to improve the range resolution of PR, the multiple broadcast channels from a single transmitter were exploited [8–10], which can be implemented only when the multiband system is used. Besides, a longer integration time can improve the Doppler resolution. But that would cause migration phenomena during integration and require complex compensation [11, 12].In recent years, the application of compressed sensing (CS) has been considered in passive radar. CS is a sparse signal processing technology [13, 14], which can reduce the amount of sampling data and use a small number of measurements to achieve excellent signal reconstruction. CS has been widely used in the field of radar signal processing; it has great potential in improving resolution [15, 16]. With the increasing attention in the field of PR, the RD map of targets for PR can be formulated as a sparse recovery problem, and the works related to CS-based PR processing have been published [17–19]. In [17, 18], a normal CS-based PR processing scheme was proposed to achieve better range and Doppler resolutions. Due to the fact that the dictionary is composed of the template signals with discrete delays and Doppler shifts, this method has high reconstruction accuracy, but it needs a large amount of storage and calculations. [19] proposed a RD map generation algorithm for PR. The extended orthogonal matching pursuit (EOMP) algorithm is employed to obtain the Doppler frequency with an improved resolution and a reduced sidelobe level. But one-dimensional cross-correlation is used to obtain the range compressed profile, which is the same as the classical method.In this paper, we propose a novel range-Doppler processing method for PR. CS is applied to the range domain and Doppler domain, respectively, which no longer requires huge storage space occupied by the dictionary. In addition, our contributions are the following. (1) In the Doppler dimension, considering the sparse characteristics of the signal, we present a modified OMP algorithm based on the multiple measurement vector (MMV) [20] model to improve the processing ability. (2) In the range dimension, we use the Fourier dictionary in the range-frequency domain, which can easily achieve high-resolution range estimation. Furthermore, a global search is to find the most relevant atom of the dictionary matrix for sparse reconstruction of the RD map, which can be treated as another way to use EOMP.The rest of this paper is organised as follows. A brief signal model of PR is introduced in Section2. In Section 3, a novel theoretical derivation of CS-based RD map generation for PR is presented. The comparison with other processing methods is also discussed in this section. Section 4 demonstrates experimental results using the simulated data and the real data. Finally, Section 5 concludes this work.


## 2. Passive Radar Geometry and Signal Model

Figure1 schematically illustrates a typical passive radar geometry, where the system is composed of transmitting station T and receiving station R. It is assumed that there is a moving target P in the scene. RT, RR, and RL represent the transmitter-target, receiver-target, and transmitter-receiver distances. v is the velocity vector of the moving target. β is the bistatic angle. ϕ is the angle between the vector v and the bistatic angle. θT is the emission angle, and θR is the target observation angle.Figure 1
Simplified bistatic geometry for passive radar.The PR receiver consists of two channels, the reference channel and surveillance channel. The reference channel gathers a time-delayed version of the transmitted waveform, and the surveillance channel records the signals scattered from targets. Let the transmitted waveform be represented byxTt, then the signal collected by the reference channel can be written as(1)xt=ArxTt−τL+nrt,where Ar is the complex amplitude, τL=RL/c represents the time delay, and nrt represents the thermal noise in the reference channel.Admittedly, the surveillance channel also contains direct signal and multipath in practice. Here, we assume that disturbance has been removed [21–24], and then, the response of the moving target P can be expressed as(2)yt=AxTt−τpexpj2πfpt+nt,where A is the complex amplitude, τp is the bistatic time delay corresponding to the target location, fp is the Doppler frequency shift related to the target velocity, and nt is the thermal noise in the surveillance channel. According to the geometric relationship between the target and the bistatic radar system, the time delay and instantaneous Doppler frequency of the target can be expressed as [25](3)τp=RR2+RL2+2RRRLsinθR+RRc,fp=2vλcosϕcosβ2.Considering an observation scene consisting ofK scattering points, a generalized expression of the surveillance signal can be written as(4)yt=∑k=1KAkxTt−τkexpj2πfkt+nt,where Ak, fk, and τk are the complex amplitude, the bistatic time delay, and the Doppler frequency shift of the kth scattering point. To simplify the analysis, the thermal fluctuations nrt and nt are neglected, Ar is set to 1, and τL is set to 0. Then, the surveillance signal can be rewritten as(5)yt=∑k=1KAkxt−τkexpj2πfkt.It should be noted that this simplification has no significant impact. For example, the conditionτL=0 can be satisfied as long as the reference signal is added a corresponding time delay.


## 3. CS-Based Processing for PR

In this section, we present our investigation of the following CS-based method to form a range-Doppler map with passive radar data.


### 3.1. Signal Preprocessing

The surveillance signal and the reference signal first need to be divided into multiple short segments [26], as shown in Figure 2. It means that the Doppler frequency change within a segment interval is ignored.Figure 2
Signal segmentation.Lett=τ+mT (where T is the segment interval); t can be called slow time. τ represents time delay, which is called fast time. Now, the surveillance signal is expressed as follows:(6)yτ=∑k=1KAkxτ−τk+mTexpj2πfkτexpj2πfkmT≈∑k=1KAkxτ−τk+mTexpj2πfkmT,where m=0,1,⋯,M−1 is the slow time index and M denotes the number of segments. The two-dimensional discrete form of (6) can be expressed as(7)ym,n=∑k=1KAkxnΔτ−nkΔτ+mTexpj2πfkmT,where Δτ is the sampling time interval and τk=nkΔτ. n=0,1,⋯,N−1 is the fast time index, and N denotes the sample number of each segment. Similarly, the discrete representation of the reference signal is expressed as(8)xm,n=xnΔτ+mT.


### 3.2. Sparse Reconstruction of Doppler Domain

By taking an FFT ofxm,n with respect to n, we have(9)Fxnm,l=∑n=0N−1xm,nexp−j2πlΔfnnΔτ=∑n=0N−1xm,nexp−jl2πNn,whereΔfn=1/NΔτ is the range-frequency bin size and l is the range-frequency index. Similarly, after performing an FFT of ym,n with respect to n, we can obtain(10)Fynm,l=∑n=0N−1ym,nexp−j2πlΔfnnΔτ=∑k=1KAkFxnm,lexp−j2πlΔfnnkΔτexpj2πfkmT=∑k=1KAkFxnm,lexp−jl2πNnkexpj2πfkmT.In order to make a sparse representation of the range-frequency bins, a Doppler dictionaryΨfl∈ℂM×Mgrid is constructed as(11)Ψfl=ψ1l⋯ψql⋯ψMgridl,ψql=Fxnm,lexpj2πMgridqm,m=0,1,⋯,M−1,where Mgrid denotes the number of the Doppler grid. And then, (10) can be expressed as(12)Fynm,l=∑k=1KAkFxnm,lexp−jl2πNnkexpj2πMgridqkm,where qk=fk/Δfq and Δfq=1/MgridT represent the size of the Doppler grid. From this, we may know that Mgrid=M makes the same resolution level as the classical method, Δfq=Δf=1/MT, and the larger Mgrid can generate higher Doppler resolution.For each range-frequency binl, (12) can be rewritten as(13)Fynl=Ψflal,where(14)Fynl=Fyn0,l⋯⋯FynM−1,lT,al=a1⋯aq⋯aMgridT.There areK nonzero elements in al, when the scattering points make different Doppler frequency shifts. And then, the positions and amplitude values are, respectively,(15)qk=fkΔfq,k=1,2,⋯,K,aqk=Akexp−jl2πNnk,k=1,2,⋯,K.If there are the same Doppler frequencies, then the number of nonzero elements becomes smaller. In turn, the complex amplitude is the sum of all coefficients related to the same Doppler.According to the CS theory, the sparsity of the signal ensures the feasibility of reducing the amount of data. The sensing matrix obtained by multiplying the measurement matrixΦ and the Doppler dictionary (sparse basis matrix) Ψfl needs to satisfy the restricted isometric property (RIP). The commonly used measurement matrices are the random Gaussian matrix and partial random unit matrix. In order to facilitate the practical operation, we adopt the method of random extraction of data, which can be expressed as(16)Fyn_subl=ΦmFynl=ΦmΨflal=Ψf_sublal,where Φm∈ℝMsub×M is a partial random unit matrix. In order to obtain the coefficient vector al, (16) needs to be solved. Due to the sparseness of the coefficient vector, it is equivalent to solving the following minimum norm problem:(17)a∧l=argminal0s.t.Fyn_subl=Ψf_sublal.There are many methods to solve (17). The greedy algorithm is widely used in practical application because of its excellent geometric interpretation, good reconstruction effect, and fast reconstruction speed. The most representative greedy algorithm is the OMP algorithm. Considering the consistency of signal models of multiple range-frequency bins, the same operation can be carried out for each range-frequency bin. That is to say, it can solve al through OMP under the single measurement vector (SMV) model according to Fyn_subl and Ψf_subl, respectively.We note that 13 is independent ofl. It is found that the positions of nonzero elements of al in multiple range-frequency bins are the same; that is, the support set of each sparse coefficient vector is the same. This feature means that it can be considered an MMV model (see Figure 3), and we can use the joint sparsity to improve reconstruction performance. However, the existing OMP algorithm under the MMV model cannot be directly applied. Considering that Ψf_subl is different for passive radar data, a modified version of MMV-OMP is proposed here. The pseudocode is shown in Algorithm 1. In order to facilitate the derivation, the matrix a=a∧1,⋯,a∧l,⋯,a∧N which is obtained by the modified MMV-OMP algorithm will be rewritten as Fymq,l.Figure 3
MMV model for passive radar.Algorithm 1: Pseudocode of modified MMV-OMP algorithm for Doppler reconstruction.
Input:  surveillance signal Fyn_sub∈ℂMsub×NdictionaryΨf_sub=Ψfsub1,⋯,Ψfsubl,⋯,ΨfsubN∈ℂMsub×MgridNsparsityKresidual thresholdTh1Initialize: iteration count k=0residual matrixR=Fyn_subestimate support collectionΛ=∅coefficient matrixa=0Mgrid×NWhile  k≤K or RkF>Th1 do k=k+1(Identification)vk⟵argmax∑l=1NΨf_sublHFyn_subl, where vk isthe column index of the largest element inΨf_sublHFyn_subl(Update Index Support)Λk=Λk−1∪vk;(Estimation)w^Λk=Ψf_sublΛk†Fyn_subl, where ·†represents the pseudo-inverse of the matrix,aΛk,n=w^Λk,(Update Residual)Rkl=Fyn_subl−ΨfsublΛkw^Λk.Output  coefficient matrix aEnd


### 3.3. Sparse Reconstruction of Range Domain

After Doppler reconstruction, the two-dimensional data can be expressed as(18)Fymq,l=∑kAkexp−jl2πNnk,q=qk,0,q≠qk.It can be seen from (18) that the signal is independent of Doppler frequency bin q. Therefore, the same time-delay dictionary Ψτ∈ℂN×Ngrid can be constructed for each Doppler bin, which can be expressed as(19)Ψτ=ψ1⋯ψn′⋯ψNgrid,ψn′=exp−j2πNgridn′l,l=0,1,⋯,N−1,whereNgrid denotes the number of time-delay (range) grids. So, the signal of qth Doppler bin can be expressed as(20)Fymq=Ψταq,where(21)Fymq=Fymq,0⋯⋯Fymq,N−1T,αq=α1⋯αn′⋯αNgridT.In whole range-Doppler plane, there areK nonzero elements K<<MgridNgrid, and the positions and amplitude values are, respectively,(22)n′k=τkΔτ′,k=1,2,⋯,K,αn′k=Ak,k=1,2,⋯,K,where Δτ′=ΔτN/Ngrid. It was obvious that Ngrid=N makes the same or better resolution level as the classical method; the larger Ngrid can generate finer range resolution.Similarly, the signal is randomly extracted, and the observation equation can be expressed as(23)Fym_subq=ΦnFymq=ΦnΨταq=Ψτ_subαq,where Φn∈ℝNsub×N is a partial random unit matrix. It is important to note that the randomness has a constraint. For range-frequency bins, the frequency range is determined by the sampling rate fs, which usually satisfies fs≥B (B is the signal bandwidth). This means that some range-frequency bins contain invalid information. Therefore, the random extraction only considers the effective part of the signal bandwidth.The coefficient vectorαq can be solved as follows:(24)α∧q=argminαq0s.t.Fym_subq=Ψτ_subαq.To effectively solve the problem, we utilize a sparse matrix recovery algorithm. Its pseudocode is described in Algorithm2, which can be seen as another way of using EOMP. Similarly, the matrix α=α∧1,⋯,α∧q,⋯,α∧MgridT is rewritten as Iymq,n′.Algorithm 2: Pseudocode of EOMP algorithm for range reconstruction.
Input:  surveillance signal after Doppler processingFym_subT∈ℂNsub×MgriddictionaryΨτ_sub∈ℂNsub×NgridsparsityQresidual thresholdTh2Initialize: iteration count q=0residual matrixr=Fym_subTestimate support collectionL=0Q×Mgridcoefficient matrixα′=0Ngrid×MgridWhile  q≤Q or rkF>Th2 do q=q+1(Identification)uq,bq⟵argmaxΨτ_subHr,where uqandbq are the row index and column index of the largestelement inΨτ_subHr;(Update Index Support)Lq,bq=uq;P=nonzerosLbq,where nonzeros• representsnon-zero elements in the vector;(Estimation)γ=Ψτ_subP†Fym_subbq, where ·†represents the pseudo-inverse of the matrix,α′P,bq=γ;(Update Residual)rbq=Fym_subbq−Ψτ_subPγ.Output  coefficient matrix α=α′TEndThe signal processed by the reconstruction algorithm can be expressed as(25)Iymq,n′=Ak,q=qk,n′=nk′,0,others.Obviously, (25) is the distribution of the scattering coefficient in the time-delay Doppler grid. Therefore, based on the above analysis, we can see that the range-Doppler map for passive radar can be obtained by using the proposed CS-based processing method.


### 3.4. Summary of Method Flow

In order to intuitively show the processing technique, the flowchart of the proposed method is shown in Figure4. The steps are briefly summarized as follows.Step 1.
After segmentation of reference signal and surveillance signal, perform FFT onxm,n and ym,n in the fast time n direction, respectively.Step 2.
According toΦm, construct Doppler dictionary Ψf_sub by using the reference signal of the range-frequency domain and observation vector Fyn_sub by using the surveillance signal of the range-frequency domain, and then perform Algorithm 1 (MMV-OMP) to reconstruct Doppler frequency distribution.Step 3.
According toΦn, construct time-delay dictionary Ψτ_sub and observation matrix Fym_sub, and then perform Algorithm 2 (EOMP) to reconstruct range-Doppler map.Figure 4
Flowchart of the proposed method.There are two points which remain to be explained. In Step2, we actually only need to process Nsub range-frequency bins related to Φn. In addition, the grid sizes (Δfq, Δτ′) and grid numbers (Mgrid, Ngrid) of dictionaries (Ψf_sub, Ψτ_sub) can be set in accordance with the actual condition.


### 3.5. Comparative Analysis

We compare the proposed method with the time-delay/Doppler combination dictionary-based CS-PR method (named as original method here) presented in [17, 18] and Feng et al.’s method presented in [19]. Table 1 demonstrates the signal model, dictionary size of the original method, Feng et al.’s method, and the proposed method, respectively.Table 1
CS-PR method comparison.
MethodsOriginal methodFeng et al.’s methodProposed methodSignal modelrD:yvec=ΨrDαvecD:yr=ΨDαD:Fyn_subl=Ψf_sublalr:Fym_sub=Ψτ_subαDictionary sizeΨrDMsubNsub×MgridNgridΨDMsub×MgridD:Ψf_sublMsub×Mgridr:Ψτ_subNsub×NgridThe major difference among these three methods lies in the signal model and method implementation. In the original method, the reflectivity map matrix is reconstructed by a cascade of 1-D CS reconstruction. All measurements are stacked into a single observation vectoryvec, the state of each time-delay/Doppler combination is stacked into the state vector αvec, and the dictionary ΨrD by discretizing the delay Doppler plane on a grid takes up the most memory. The size is MsubNsub×MgridNgrid, which leads to the memory occupation being too large and the computational burden being huge.In both Feng et al.’s method and the proposed algorithm, the range reconstruction and the Doppler reconstruction are separately completed, which means that the dictionariesΨD, Ψf_subl, Ψτ_sub have much smaller size. The difference between the two methods is that the former only uses CS to estimate Doppler frequency while the latter uses CS in both directions, even though the use of EOMP is sameness. We can obtain high-resolution capability in the range coordinate. This is just what Feng et al.’s method does not have.In addition, we can consider the fact that the input of EOMP is the output of MMV-OMP, which records the support set information. Therefore, the input data size of EOMP can be adjusted according to the size of the support set. As shown in Figure5, there are two execution modes to obtain the range-Doppler map. The second execution mode can further reduce the computational burden due to the small amount of data.Figure 5
The diagram of signal reconstruction.


## 3.1. Signal Preprocessing

The surveillance signal and the reference signal first need to be divided into multiple short segments [26], as shown in Figure 2. It means that the Doppler frequency change within a segment interval is ignored.Figure 2
Signal segmentation.Lett=τ+mT (where T is the segment interval); t can be called slow time. τ represents time delay, which is called fast time. Now, the surveillance signal is expressed as follows:(6)yτ=∑k=1KAkxτ−τk+mTexpj2πfkτexpj2πfkmT≈∑k=1KAkxτ−τk+mTexpj2πfkmT,where m=0,1,⋯,M−1 is the slow time index and M denotes the number of segments. The two-dimensional discrete form of (6) can be expressed as(7)ym,n=∑k=1KAkxnΔτ−nkΔτ+mTexpj2πfkmT,where Δτ is the sampling time interval and τk=nkΔτ. n=0,1,⋯,N−1 is the fast time index, and N denotes the sample number of each segment. Similarly, the discrete representation of the reference signal is expressed as(8)xm,n=xnΔτ+mT.


## 3.2. Sparse Reconstruction of Doppler Domain

By taking an FFT ofxm,n with respect to n, we have(9)Fxnm,l=∑n=0N−1xm,nexp−j2πlΔfnnΔτ=∑n=0N−1xm,nexp−jl2πNn,whereΔfn=1/NΔτ is the range-frequency bin size and l is the range-frequency index. Similarly, after performing an FFT of ym,n with respect to n, we can obtain(10)Fynm,l=∑n=0N−1ym,nexp−j2πlΔfnnΔτ=∑k=1KAkFxnm,lexp−j2πlΔfnnkΔτexpj2πfkmT=∑k=1KAkFxnm,lexp−jl2πNnkexpj2πfkmT.In order to make a sparse representation of the range-frequency bins, a Doppler dictionaryΨfl∈ℂM×Mgrid is constructed as(11)Ψfl=ψ1l⋯ψql⋯ψMgridl,ψql=Fxnm,lexpj2πMgridqm,m=0,1,⋯,M−1,where Mgrid denotes the number of the Doppler grid. And then, (10) can be expressed as(12)Fynm,l=∑k=1KAkFxnm,lexp−jl2πNnkexpj2πMgridqkm,where qk=fk/Δfq and Δfq=1/MgridT represent the size of the Doppler grid. From this, we may know that Mgrid=M makes the same resolution level as the classical method, Δfq=Δf=1/MT, and the larger Mgrid can generate higher Doppler resolution.For each range-frequency binl, (12) can be rewritten as(13)Fynl=Ψflal,where(14)Fynl=Fyn0,l⋯⋯FynM−1,lT,al=a1⋯aq⋯aMgridT.There areK nonzero elements in al, when the scattering points make different Doppler frequency shifts. And then, the positions and amplitude values are, respectively,(15)qk=fkΔfq,k=1,2,⋯,K,aqk=Akexp−jl2πNnk,k=1,2,⋯,K.If there are the same Doppler frequencies, then the number of nonzero elements becomes smaller. In turn, the complex amplitude is the sum of all coefficients related to the same Doppler.According to the CS theory, the sparsity of the signal ensures the feasibility of reducing the amount of data. The sensing matrix obtained by multiplying the measurement matrixΦ and the Doppler dictionary (sparse basis matrix) Ψfl needs to satisfy the restricted isometric property (RIP). The commonly used measurement matrices are the random Gaussian matrix and partial random unit matrix. In order to facilitate the practical operation, we adopt the method of random extraction of data, which can be expressed as(16)Fyn_subl=ΦmFynl=ΦmΨflal=Ψf_sublal,where Φm∈ℝMsub×M is a partial random unit matrix. In order to obtain the coefficient vector al, (16) needs to be solved. Due to the sparseness of the coefficient vector, it is equivalent to solving the following minimum norm problem:(17)a∧l=argminal0s.t.Fyn_subl=Ψf_sublal.There are many methods to solve (17). The greedy algorithm is widely used in practical application because of its excellent geometric interpretation, good reconstruction effect, and fast reconstruction speed. The most representative greedy algorithm is the OMP algorithm. Considering the consistency of signal models of multiple range-frequency bins, the same operation can be carried out for each range-frequency bin. That is to say, it can solve al through OMP under the single measurement vector (SMV) model according to Fyn_subl and Ψf_subl, respectively.We note that 13 is independent ofl. It is found that the positions of nonzero elements of al in multiple range-frequency bins are the same; that is, the support set of each sparse coefficient vector is the same. This feature means that it can be considered an MMV model (see Figure 3), and we can use the joint sparsity to improve reconstruction performance. However, the existing OMP algorithm under the MMV model cannot be directly applied. Considering that Ψf_subl is different for passive radar data, a modified version of MMV-OMP is proposed here. The pseudocode is shown in Algorithm 1. In order to facilitate the derivation, the matrix a=a∧1,⋯,a∧l,⋯,a∧N which is obtained by the modified MMV-OMP algorithm will be rewritten as Fymq,l.Figure 3
MMV model for passive radar.Algorithm 1: Pseudocode of modified MMV-OMP algorithm for Doppler reconstruction.
Input:  surveillance signal Fyn_sub∈ℂMsub×NdictionaryΨf_sub=Ψfsub1,⋯,Ψfsubl,⋯,ΨfsubN∈ℂMsub×MgridNsparsityKresidual thresholdTh1Initialize: iteration count k=0residual matrixR=Fyn_subestimate support collectionΛ=∅coefficient matrixa=0Mgrid×NWhile  k≤K or RkF>Th1 do k=k+1(Identification)vk⟵argmax∑l=1NΨf_sublHFyn_subl, where vk isthe column index of the largest element inΨf_sublHFyn_subl(Update Index Support)Λk=Λk−1∪vk;(Estimation)w^Λk=Ψf_sublΛk†Fyn_subl, where ·†represents the pseudo-inverse of the matrix,aΛk,n=w^Λk,(Update Residual)Rkl=Fyn_subl−ΨfsublΛkw^Λk.Output  coefficient matrix aEnd


## 3.3. Sparse Reconstruction of Range Domain

After Doppler reconstruction, the two-dimensional data can be expressed as(18)Fymq,l=∑kAkexp−jl2πNnk,q=qk,0,q≠qk.It can be seen from (18) that the signal is independent of Doppler frequency bin q. Therefore, the same time-delay dictionary Ψτ∈ℂN×Ngrid can be constructed for each Doppler bin, which can be expressed as(19)Ψτ=ψ1⋯ψn′⋯ψNgrid,ψn′=exp−j2πNgridn′l,l=0,1,⋯,N−1,whereNgrid denotes the number of time-delay (range) grids. So, the signal of qth Doppler bin can be expressed as(20)Fymq=Ψταq,where(21)Fymq=Fymq,0⋯⋯Fymq,N−1T,αq=α1⋯αn′⋯αNgridT.In whole range-Doppler plane, there areK nonzero elements K<<MgridNgrid, and the positions and amplitude values are, respectively,(22)n′k=τkΔτ′,k=1,2,⋯,K,αn′k=Ak,k=1,2,⋯,K,where Δτ′=ΔτN/Ngrid. It was obvious that Ngrid=N makes the same or better resolution level as the classical method; the larger Ngrid can generate finer range resolution.Similarly, the signal is randomly extracted, and the observation equation can be expressed as(23)Fym_subq=ΦnFymq=ΦnΨταq=Ψτ_subαq,where Φn∈ℝNsub×N is a partial random unit matrix. It is important to note that the randomness has a constraint. For range-frequency bins, the frequency range is determined by the sampling rate fs, which usually satisfies fs≥B (B is the signal bandwidth). This means that some range-frequency bins contain invalid information. Therefore, the random extraction only considers the effective part of the signal bandwidth.The coefficient vectorαq can be solved as follows:(24)α∧q=argminαq0s.t.Fym_subq=Ψτ_subαq.To effectively solve the problem, we utilize a sparse matrix recovery algorithm. Its pseudocode is described in Algorithm2, which can be seen as another way of using EOMP. Similarly, the matrix α=α∧1,⋯,α∧q,⋯,α∧MgridT is rewritten as Iymq,n′.Algorithm 2: Pseudocode of EOMP algorithm for range reconstruction.
Input:  surveillance signal after Doppler processingFym_subT∈ℂNsub×MgriddictionaryΨτ_sub∈ℂNsub×NgridsparsityQresidual thresholdTh2Initialize: iteration count q=0residual matrixr=Fym_subTestimate support collectionL=0Q×Mgridcoefficient matrixα′=0Ngrid×MgridWhile  q≤Q or rkF>Th2 do q=q+1(Identification)uq,bq⟵argmaxΨτ_subHr,where uqandbq are the row index and column index of the largestelement inΨτ_subHr;(Update Index Support)Lq,bq=uq;P=nonzerosLbq,where nonzeros• representsnon-zero elements in the vector;(Estimation)γ=Ψτ_subP†Fym_subbq, where ·†represents the pseudo-inverse of the matrix,α′P,bq=γ;(Update Residual)rbq=Fym_subbq−Ψτ_subPγ.Output  coefficient matrix α=α′TEndThe signal processed by the reconstruction algorithm can be expressed as(25)Iymq,n′=Ak,q=qk,n′=nk′,0,others.Obviously, (25) is the distribution of the scattering coefficient in the time-delay Doppler grid. Therefore, based on the above analysis, we can see that the range-Doppler map for passive radar can be obtained by using the proposed CS-based processing method.


## 3.4. Summary of Method Flow

In order to intuitively show the processing technique, the flowchart of the proposed method is shown in Figure4. The steps are briefly summarized as follows.Step 1.
After segmentation of reference signal and surveillance signal, perform FFT onxm,n and ym,n in the fast time n direction, respectively.Step 2.
According toΦm, construct Doppler dictionary Ψf_sub by using the reference signal of the range-frequency domain and observation vector Fyn_sub by using the surveillance signal of the range-frequency domain, and then perform Algorithm 1 (MMV-OMP) to reconstruct Doppler frequency distribution.Step 3.
According toΦn, construct time-delay dictionary Ψτ_sub and observation matrix Fym_sub, and then perform Algorithm 2 (EOMP) to reconstruct range-Doppler map.Figure 4
Flowchart of the proposed method.There are two points which remain to be explained. In Step2, we actually only need to process Nsub range-frequency bins related to Φn. In addition, the grid sizes (Δfq, Δτ′) and grid numbers (Mgrid, Ngrid) of dictionaries (Ψf_sub, Ψτ_sub) can be set in accordance with the actual condition.


## 3.5. Comparative Analysis

We compare the proposed method with the time-delay/Doppler combination dictionary-based CS-PR method (named as original method here) presented in [17, 18] and Feng et al.’s method presented in [19]. Table 1 demonstrates the signal model, dictionary size of the original method, Feng et al.’s method, and the proposed method, respectively.Table 1
CS-PR method comparison.
MethodsOriginal methodFeng et al.’s methodProposed methodSignal modelrD:yvec=ΨrDαvecD:yr=ΨDαD:Fyn_subl=Ψf_sublalr:Fym_sub=Ψτ_subαDictionary sizeΨrDMsubNsub×MgridNgridΨDMsub×MgridD:Ψf_sublMsub×Mgridr:Ψτ_subNsub×NgridThe major difference among these three methods lies in the signal model and method implementation. In the original method, the reflectivity map matrix is reconstructed by a cascade of 1-D CS reconstruction. All measurements are stacked into a single observation vectoryvec, the state of each time-delay/Doppler combination is stacked into the state vector αvec, and the dictionary ΨrD by discretizing the delay Doppler plane on a grid takes up the most memory. The size is MsubNsub×MgridNgrid, which leads to the memory occupation being too large and the computational burden being huge.In both Feng et al.’s method and the proposed algorithm, the range reconstruction and the Doppler reconstruction are separately completed, which means that the dictionariesΨD, Ψf_subl, Ψτ_sub have much smaller size. The difference between the two methods is that the former only uses CS to estimate Doppler frequency while the latter uses CS in both directions, even though the use of EOMP is sameness. We can obtain high-resolution capability in the range coordinate. This is just what Feng et al.’s method does not have.In addition, we can consider the fact that the input of EOMP is the output of MMV-OMP, which records the support set information. Therefore, the input data size of EOMP can be adjusted according to the size of the support set. As shown in Figure5, there are two execution modes to obtain the range-Doppler map. The second execution mode can further reduce the computational burden due to the small amount of data.Figure 5
The diagram of signal reconstruction.


## 4. Experimental Results

In this section, we present experimental results with simulated data and real data. The effectiveness of the proposed processing method is demonstrated.


### 4.1. Simulation Data

We have conducted numerical experiments to investigate the performance of the proposed processing method. A digital television terrestrial broadcasting (DTTB) signal is simulated. The frame structure of signal is shown in Table2, which includes frame header (945 symbols) and frame body (3780 symbols).Table 2
Frame structure of DTTB signal.
Frame headSystem informationDataDTTB945 symbols36 symbols3744 symbolsThe parameters used in the simulation experiment are shown in Tables3 and 4. In the observation scene, the reflection mechanisms are assumed to be a point-like target (target 1) and a line-like target (target 2). Target 2 consists of three scatterer points, which are located in the same Doppler bin and become neighbors in the range direction.Table 3
System parameter.
ParametersSymbolValueCarrier frequencyfc674 MHzBandwidthB7.56 MHzSample frequencyfs10 MHzNumber of segmentsM512Sample numberN256Number of range gridMgrid50Number of Doppler gridNgrid20Table 4
Target scenario parameters.
ParametersTarget 1Target 2Scatterer point 1Scatterer point 2Scatterer point 3Target location index10262830Doppler bin4-4-4-4Amplitude10.70.91At the beginning, the result obtained by the classical CAF is shown in Figure6. As can be seen from Figures 6(a) and 6(b), the CAF can achieve the target scene recovery, but the reconstructed map is out of clarity due to the large sidelobes. On the contrary, the CS-PR methods can be used to remove the sidelobes. The processing results by using the CS-PR methods listed in Table 1 are shown in Figure 7, which are the RD maps obtained by using full samples and partial samples from the same scenario. For the partial sample case, in order to ensure the same amount of data, the original method and the proposed method perform one-half data undersampling processing in the range direction and Doppler direction, respectively, while Feng et al.’s method only performs one-quarter data undersampling processing in the Doppler direction.Figure 6
Result obtained by CAF: (a) the original scene; (b) RD map.
(a)(b)Figure 7
RD maps obtained by (a, b) original CS-PR method,K=4; (c, d) Feng et al.’s method, K=4; (e, f) proposed method, K=Q=4; (a, c, e) with full samples; (b, d, f) with partial samples (1/4).
(a)(b)(c)(d)(e)(f)Figure7 demonstrates that these processing methods based on CS can effectively suppress the sidelobe. The run times of different methods are summarized in Table 5. Apparently, the original CS-PR method uses much longer running time than the other two methods, due to the usage of the dictionary by the discretizing RD plane. In order to more intuitively observe the effect of different methods on the sidelobe suppression, the range profiles obtained by different methods are shown in Figure 8. From the figures, it shows that Feng et al.’s method cannot clearly identify the scatterer points in the range direction. It indicates that the range resolution of the method is limited, even though its run time is the least. Fortunately, the proposed processing method can reconstruct the target scene with a good performance both in accuracy and in computational efficiency.Table 5
Run times for different CS-PR methods.
MethodsFull samplesPartial samplesOriginal method3.365 s0.615 sFeng et al.’s method0.069 s0.054 sProposed method0.161 s0.123 sFigure 8
Range profiles obtained by different methods: (a) target 1; (b) target 2.
(a)(b)It is known that the CS approach offers great potential for better resolution by using a finer dictionary. The grid sizes are set to be the same as CAF in previous experiments (Δfq=Δf, Δτ′=Δτ). Now, the Doppler grid and time-delay grid are only half the original size (Δfq=Δf/2, Δτ′=Δτ/2). As the grid is refined, the challenge is that the dictionary correlation will increase, which may lead to the performance degradation for CS. The Gram matrix is used to verify the dictionary coherence, and results are depicted in Figures 9 and 10. Gf=Ψf_subHΨf_sub and Gτ=Ψτ_subHΨτ_sub represent Gram matrices in Doppler direction and range direction, respectively. As can be seen from Figures 9 and 10, each Gram matrix is close to the unit matrix. Due to this coherence characteristic, CS is able to produce superresolution radar images. Figure 11 shows the RD map obtained by the proposed method when the grid is refined. It is observed that the resolution improvement is achieved by utilizing a fine grid.Figure 9
Gram matrices of original dictionaries (Δfq=Δf, Δτ′=Δτ): (a) Gf; (b) Gτ.
(a)(b)Figure 10
Gram matrices of refined dictionaries (Δfq=Δf/2, Δτ′=Δτ/2): (a) Gf; (b) Gτ.
(a)(b)Figure 11
Result obtained by the proposed method with fine grid: (a) refining grid with2×2; (b) RD map.
(a)(b)In order to demonstrate the advantages of OMP under the MMV model, we compare the Doppler reconstructions by SMV-OMP and MMV-OMP at different signal-to-noise ratios (SNR). The values of SNR are10dB, 0dB, and −10dB. The parameters of target 1 use the following settings. The location index is 10, the Doppler bin is 4, and the amplitude is 0.3. Other parameters remain unchanged. The results are shown in Figures 12–14. For SMV-OMP, there are a number of Doppler reconstruction errors, which will lead to insufficient energy accumulation in the range direction. Figure 15 shows the final RD maps when SNR is −10dB. It can be found that target 1 is not visible in the RD map obtained by SMV-OMP. On the contrary, MMV-OMP has robust performance because it considers the information of multiple observations.Figure 12
Results of sparse reconstruction of Doppler domain (SNR=10dB) by (a) SMV-OMP and (b) MMV-OMP.
(a)(b)Figure 13
Results of sparse reconstruction of Doppler domain (SNR=0dB) by (a) SMV-OMP and (b) MMV-OMP.
(a)(b)Figure 14
Results of sparse reconstruction of Doppler domain (SNR=−10dB) by (a) SMV-OMP and (b) MMV-OMP.
(a)(b)Figure 15
RD maps (SNR=−10dB) obtained by (a) SMV-OMP and (b) MMV-OMP.
(a)(b)


### 4.2. Real Data

We use real data from an DTV-based PR system to further verify the proposed method. The parameters used are shown in Table6. We obtain results for CAF and the proposed method. The RD map obtained by the CAF approach is shown in Figure 16; Figures 17 and 18 show the results of CS reconstruction if 25% of the full data is used. In comparison to the CAF, the proposed method shows good performance owing to CS.Table 6
System parameters.
ParametersSymbolValueCarrier frequencyfc674 MHzBandwidthB7.56 MHzSample frequencyfs10 MHzNumber of segmentsM2048Sample numberN2048Number of range gridMgrid50Number of Doppler gridNgrid200Figure 16
RD map obtained by CAF.Figure 17
RD map obtained by the proposed method.Figure 18
RD map obtained by the proposed method with fine grid.For a more illustrative comparison, two cuts are made along the location of the target in the range and Doppler directions shown in Figures19(a) and 19(b). The proposed method has a considerably lower sidelobe level than the CAF and indeed allows for improved range and Doppler resolutions by the fine gridding.Figure 19
Profiles of target: (a) in range direction; (b) in Doppler direction.
(a)(b)


## 4.1. Simulation Data

We have conducted numerical experiments to investigate the performance of the proposed processing method. A digital television terrestrial broadcasting (DTTB) signal is simulated. The frame structure of signal is shown in Table2, which includes frame header (945 symbols) and frame body (3780 symbols).Table 2
Frame structure of DTTB signal.
Frame headSystem informationDataDTTB945 symbols36 symbols3744 symbolsThe parameters used in the simulation experiment are shown in Tables3 and 4. In the observation scene, the reflection mechanisms are assumed to be a point-like target (target 1) and a line-like target (target 2). Target 2 consists of three scatterer points, which are located in the same Doppler bin and become neighbors in the range direction.Table 3
System parameter.
ParametersSymbolValueCarrier frequencyfc674 MHzBandwidthB7.56 MHzSample frequencyfs10 MHzNumber of segmentsM512Sample numberN256Number of range gridMgrid50Number of Doppler gridNgrid20Table 4
Target scenario parameters.
ParametersTarget 1Target 2Scatterer point 1Scatterer point 2Scatterer point 3Target location index10262830Doppler bin4-4-4-4Amplitude10.70.91At the beginning, the result obtained by the classical CAF is shown in Figure6. As can be seen from Figures 6(a) and 6(b), the CAF can achieve the target scene recovery, but the reconstructed map is out of clarity due to the large sidelobes. On the contrary, the CS-PR methods can be used to remove the sidelobes. The processing results by using the CS-PR methods listed in Table 1 are shown in Figure 7, which are the RD maps obtained by using full samples and partial samples from the same scenario. For the partial sample case, in order to ensure the same amount of data, the original method and the proposed method perform one-half data undersampling processing in the range direction and Doppler direction, respectively, while Feng et al.’s method only performs one-quarter data undersampling processing in the Doppler direction.Figure 6
Result obtained by CAF: (a) the original scene; (b) RD map.
(a)(b)Figure 7
RD maps obtained by (a, b) original CS-PR method,K=4; (c, d) Feng et al.’s method, K=4; (e, f) proposed method, K=Q=4; (a, c, e) with full samples; (b, d, f) with partial samples (1/4).
(a)(b)(c)(d)(e)(f)Figure7 demonstrates that these processing methods based on CS can effectively suppress the sidelobe. The run times of different methods are summarized in Table 5. Apparently, the original CS-PR method uses much longer running time than the other two methods, due to the usage of the dictionary by the discretizing RD plane. In order to more intuitively observe the effect of different methods on the sidelobe suppression, the range profiles obtained by different methods are shown in Figure 8. From the figures, it shows that Feng et al.’s method cannot clearly identify the scatterer points in the range direction. It indicates that the range resolution of the method is limited, even though its run time is the least. Fortunately, the proposed processing method can reconstruct the target scene with a good performance both in accuracy and in computational efficiency.Table 5
Run times for different CS-PR methods.
MethodsFull samplesPartial samplesOriginal method3.365 s0.615 sFeng et al.’s method0.069 s0.054 sProposed method0.161 s0.123 sFigure 8
Range profiles obtained by different methods: (a) target 1; (b) target 2.
(a)(b)It is known that the CS approach offers great potential for better resolution by using a finer dictionary. The grid sizes are set to be the same as CAF in previous experiments (Δfq=Δf, Δτ′=Δτ). Now, the Doppler grid and time-delay grid are only half the original size (Δfq=Δf/2, Δτ′=Δτ/2). As the grid is refined, the challenge is that the dictionary correlation will increase, which may lead to the performance degradation for CS. The Gram matrix is used to verify the dictionary coherence, and results are depicted in Figures 9 and 10. Gf=Ψf_subHΨf_sub and Gτ=Ψτ_subHΨτ_sub represent Gram matrices in Doppler direction and range direction, respectively. As can be seen from Figures 9 and 10, each Gram matrix is close to the unit matrix. Due to this coherence characteristic, CS is able to produce superresolution radar images. Figure 11 shows the RD map obtained by the proposed method when the grid is refined. It is observed that the resolution improvement is achieved by utilizing a fine grid.Figure 9
Gram matrices of original dictionaries (Δfq=Δf, Δτ′=Δτ): (a) Gf; (b) Gτ.
(a)(b)Figure 10
Gram matrices of refined dictionaries (Δfq=Δf/2, Δτ′=Δτ/2): (a) Gf; (b) Gτ.
(a)(b)Figure 11
Result obtained by the proposed method with fine grid: (a) refining grid with2×2; (b) RD map.
(a)(b)In order to demonstrate the advantages of OMP under the MMV model, we compare the Doppler reconstructions by SMV-OMP and MMV-OMP at different signal-to-noise ratios (SNR). The values of SNR are10dB, 0dB, and −10dB. The parameters of target 1 use the following settings. The location index is 10, the Doppler bin is 4, and the amplitude is 0.3. Other parameters remain unchanged. The results are shown in Figures 12–14. For SMV-OMP, there are a number of Doppler reconstruction errors, which will lead to insufficient energy accumulation in the range direction. Figure 15 shows the final RD maps when SNR is −10dB. It can be found that target 1 is not visible in the RD map obtained by SMV-OMP. On the contrary, MMV-OMP has robust performance because it considers the information of multiple observations.Figure 12
Results of sparse reconstruction of Doppler domain (SNR=10dB) by (a) SMV-OMP and (b) MMV-OMP.
(a)(b)Figure 13
Results of sparse reconstruction of Doppler domain (SNR=0dB) by (a) SMV-OMP and (b) MMV-OMP.
(a)(b)Figure 14
Results of sparse reconstruction of Doppler domain (SNR=−10dB) by (a) SMV-OMP and (b) MMV-OMP.
(a)(b)Figure 15
RD maps (SNR=−10dB) obtained by (a) SMV-OMP and (b) MMV-OMP.
(a)(b)


## 4.2. Real Data

We use real data from an DTV-based PR system to further verify the proposed method. The parameters used are shown in Table6. We obtain results for CAF and the proposed method. The RD map obtained by the CAF approach is shown in Figure 16; Figures 17 and 18 show the results of CS reconstruction if 25% of the full data is used. In comparison to the CAF, the proposed method shows good performance owing to CS.Table 6
System parameters.
ParametersSymbolValueCarrier frequencyfc674 MHzBandwidthB7.56 MHzSample frequencyfs10 MHzNumber of segmentsM2048Sample numberN2048Number of range gridMgrid50Number of Doppler gridNgrid200Figure 16
RD map obtained by CAF.Figure 17
RD map obtained by the proposed method.Figure 18
RD map obtained by the proposed method with fine grid.For a more illustrative comparison, two cuts are made along the location of the target in the range and Doppler directions shown in Figures19(a) and 19(b). The proposed method has a considerably lower sidelobe level than the CAF and indeed allows for improved range and Doppler resolutions by the fine gridding.Figure 19
Profiles of target: (a) in range direction; (b) in Doppler direction.
(a)(b)


## 5. Conclusion

In this paper, the problem of compressed sensing-based range-Doppler processing for passive radar is investigated. In order to reduce the sidelobes and improve the resolution, we have proposed a novel CS-PR method, in which the modified MMV-OMP algorithm is used to perform sparse reconstruction of Doppler dimension, and then, the EOMP algorithm is used to perform sparse reconstruction of range dimension. Compared to previous CS-PR methods, we can achieve a high-quality reconstruction of the range-Doppler map of target scenario and do not suffer from the heavy computational burden. The effectiveness of the proposed method is verified by experiments with simulated data and real data. The improved resolution capability will be helpful to widen the extent of application.



---

*Source: 5570498-2021-07-23.xml*