OTHR Spectrum Reconstruction of Maneuvering Target with Compressive Sensing

High-frequency (HF) over-the-horizon radar (OTHR) works in a very complicated electromagnetic environment. It usually suffers performance degradation caused by transient interference. In this paper, we study the transient interference excision and full spectrum reconstruction ofmaneuvering targets.The segmental subspace projection (SP) approach is applied to suppress the clutter and locate the transient interference. After interference excision, the spectrum is reconstructed from incomplete measurements via compressive sensing (CS) by using a redundant Fourier-chirp dictionary. An improved orthogonal matching pursuit (IOMP) algorithm is developed to solve the sparse decomposition optimization. Experimental results demonstrate the effectiveness of the proposed methods.


Introduction
High-frequency (HF) sky wave over-the-horizon radar (OTHR) can provide long-range detection of targets in large surveillance areas [1][2][3][4].However, OTHR signal is usually contaminated by impulsive interference, such as radio frequency interference, industrial interference, impulsive noise, lightning impulsive, and meteor trail echoes [5,6].Transient interference usually only exists in several pulse repetition periods, but its intensity is great.It will submerge the target signal in the Doppler frequency domain and make target detection fail.Filtering out these interferences is available [7,8], while it destroys the integrity of signal in the time domain.What is more, modern OTHR systems are multimode and expected to detect multitargets, simultaneously.To satisfy these requirements, discontinuous sampling mode is usually applied which results in parts of signal missing for a single target.Interference excision and data interpolation have been proposed to mitigate the interference and reconstruct OTHR signals [9,10].Prediction model based interpolation and extrapolation algorithms can solve the data missing problem in some situations, where the available data is fitted into a linear prediction model.And the missed data is able to be recovered by the model coefficients and available observations [11][12][13].
Recently, a new framework, referred to as compressive sensing (CS), for simultaneous sensing and compression has received considerable attention and has been successfully applied in many fields, such as signal/image processing, radar imaging, communication, geophysics, and remote sensing [14][15][16][17][18][19].In CS, it has been shown that a signal which is sparse or has a sparse representation in some bases can be recovered from a small number of random nonadaptive linear measurements.The original signal can be reconstructed with incomplete measurement data through convex optimization that uses the sparseness as a priori information.The capability of CS to reconstruct a sparse signal from far fewer nonadaptive measurements provides a new perspective for data reduction in spectrum recovery of targets without signal distortion.There have been several approaches that apply CS for radar signal processing [20][21][22][23][24][25][26].Based on this idea, the spectrum recovery of OTHR signal with incomplete measurements can be converted into an optimization of sparse signal reconstruction, where the sparse priority of targets in Doppler domain is exploited to reconstruct the missing data [27].Current literatures show that OTHR signal of stably moving targets can be reconstructed based on Fourier dictionary via sparse optimization in an optimal manner. 2 International Journal of Antennas and Propagation Generally speaking, most of the currently existing methods for signal reconstruction deal with target echoes as multicomponent complex sinusoidal signals.For signal reconstruction, the sinusoidal signal assumption is only suitable to express the echoed signal from stably moving target within a relatively short coherent integral time (CIT).However, in OTHR, the pulse repetition frequency (PRF) is usually very low; in terms of achieving high Doppler resolution, a long CIT is necessary.For example, when the PRF is 100 Hz and CIT contains 1024 pulses, then the time period of CIT is 10.24 s.In this case, for noncooperative maneuvering targets, such as flying aircrafts and vessels on the sea, serious maneuver is involved in their motion and in yielding time-variant Doppler modulation in the echoed signal.In terms of maneuvering target imaging and detection, the signal components should be modeled as chirps or high-order frequency modulated signals [28,29].In this case, multiple-sinusoidal model is maybe insufficient to characterize the maneuvering moving target only by Doppler frequency.
In this paper, we extend the application of CS on OTHR signal reconstruction, which was introduced in [27] in dealing with spectrum recovery for stably moving targets.Specifically, we consider the signal reconstruction of maneuvering targets which was contaminated by transient interferences.First, the segmental subspace projection (SP) and adaptive threshold are adopted for transient interference excision.Then, a novel algorithm is proposed to reconstruct the complete OTHR signal of maneuvering targets with partial measurements.The integrated Doppler spectrum is reconstructed based on CS by using particular Fourier-chirp dictionary.In general, the number of Doppler cells of targets is much smaller than the number of entire Doppler cells in the spectrum domain.Thus, the OTHR signal can be regarded as sparse in Doppler domain.Under a redundant Fourier-chirp basis, the signal of maneuvering targets can be sparsely represented.Then by sparse optimization, we can reconstruct the signal of targets.The spectrum reconstruction problem is shifted into a problem of a  1 -norm constrained optimal problem.By solving this problem efficiently, the orthogonal matched pursuit (OMP) algorithm is refined to reconstruct the integrate frequency spectrum of OTHR targets.
This paper is organized as follows.In Section 2, the signal model of OTHR is briefly introduced.Then, we analyze the characteristics of transient interference and present the segmental SP method for interference excision.In Section 3, we formulate the sparse optimization for full spectrum reconstruction in the framework of CS and extensions to Fourier-chirp dictionary for maneuvering targets are discussed in detail.What is more, the improved OMP (IOMP) algorithm for signal reconstruction is given.Section 4 gives the real-data results to show the effectiveness and advantages of the proposed method.Finally, some conclusions are summarized in the last section.

Transient Interference Excision
Suppose that the radar transmits a burst of  coherent pulses in a CIT.After range compression, the echoed signal in a certain azimuth-range cell within the CIT can be modeled as where (  ) is the ground and ocean clutter and (  ) is the target signal.(  ) and (  ) denote the transient interference and the additive noise, respectively.The target motion will produce Doppler modulation on radar signal.Assuming a certain range cell contains  stably moving targets with different velocities, we have the signal within the CIT as follows: where   and   = 2V  / are the th target's backward scattering amplitude and Doppler frequency, respectively.V  denotes the velocity of th target and  is wavelength.By applying the conventional Doppler analysis to (  ), we have   = Φx, where Φ is the Fourier transform matrix and   is the Doppler spectrum of targets, clutter, transient interference, and noise.However, for maneuvering targets, serious maneuver is involved in their motion and in yielding time-variant Doppler modulation in the echoed signal; then the traditional spectral analysis method based on Fourier transform will no longer apply.In terms of maneuvering target detection, the signal components should be modeled as chirps or high-order frequency modulated signals.Considering the acceleration of moving targets, the signal is rewritten as follows: where   and   = 2  / are Doppler frequency and chip rate of th target, respectively.  represents the radial acceleration of th target.
In OTHR, the dominant contribution of echo energy is produced by ground and ocean clutter.The transient interference () may only occur in a shorter portion of entire CIT, whose energy is usually much stronger than that of target signals but weaker than the energy of clutter.Thus, the clutter should be removed prior to the detection and excision of transient interference because of its strong intensity.Considering these characteristics of transient interference, the signal of interference can be defined as where  is the number of transient interference and   and   denote the duration and position of th interference.  is the complex envelope of th interference.
In general, ground clutter has a very narrow spectrum and the frequency domain excision method can work well in this case.However, sea clutter has a relatively wide International Journal of Antennas and Propagation spectrum, which may overlap partly with the spectrum of slow moving targets (e.g., ships) or transient interference.But sea clutter exhibits very strong correlation among adjacent range and azimuth cells, while it is not true for transient interference.This feature can be exploited for adaptive sea clutter suppression.Assume x denotes the interested azimuth-range cell (referred to as primary data).We can take the data from  range and azimuth cells in the vicinity of the primary data as secondary data to estimate the sea clutter covariance matrix of the interested cell.First, we divide CIT into  shorter nonoverlapping segments and each segment contains  ( = /) samples for segmental SP processing.The choice of  depends on the following factors: (1) the degree of freedom needed to suppress the clutter and (2) the affordable computational cost.Define the primary data x = [x  1 , x  2 , . . ., x  ℎ , . . ., x   ]  and x ℎ = [((ℎ − 1) + 1), . . ., (ℎ)]  , where (⋅)  denotes transpose operation.The data from adjacent  azimuth-range cells are collected as reference to estimate the clutter covariance matrix.Let x ℎ, = [  ((ℎ − 1) + 1), . . .,   (ℎ)]  denote the data vector of the th ( = 1, 2, . . ., ) reference cell.Then the clutter covariance matrix Rℎ can be estimated by where (⋅)  denotes conjugate transpose operation.By making eigenvalue decomposition (EVD) on Rℎ , we have where   and s  denote the th eigenvalue and eigenvector, respectively.Assume  1 >  2 > ⋅ ⋅ ⋅ >    and there are  principal eigenvalues corresponding to the clutter.Then the  principal eigenvectors span the clutter subspace and the rest eigenvectors span the signal (including targets, transient interference, and noise) subspace.The clutter subspace can be expressed as where S ℎ is composed of  principal eigenvectors corresponding to the  dominating eigenvalues.Define the ℎth signal SP matrix as By projecting the primary data x into the signal subspace we can get the signal after clutter suppression as x ×1 = Px.Generally, transient interference is much stronger than the noise level.After clutter removal, transient interference emerges in x; then by comparing to an adaptive amplitude threshold, we can locate the position and determine the duration of the transient interference.The threshold  is given by  = mean an ⋅ , where mean an represents the mean amplitude of the samples within CIT and  is a constant scale factor used to achieve a desired constant false alarm probability for transient interference detection.It is set as 10− 15 dB experientially.After clutter mitigation and interference excision, the signal becomes where x  ×1 ( < ) is the signal sequence after discarding the samples which are contaminated by the interference.Filtering out the clutter and interference will cause signal missing, while it destroys the integrity of signal in time domain.The process can be explained in Figure 1.The target signal suffers from sidelobe leakage due to data missing by interference excision, which will degrade the moving target detection performance.Data interpolation can be used to remedy this effect.AR model is typically used to predict the missing data and the AR coefficients can be calculated efficiently using the well-known Burg algorithm [30].Also, we can reconstruct the integrate frequency spectrum of OTHR signal by solving  1 -norm constrained optimal problem.

Full Spectrum Reconstruction
The developing theory of compressed sensing indicates that an unknown sparse signal can be exactly recovered from a very limited number of measurements with high probability by solving an optimization problem when some special conditions are met [14][15][16].Considering y ∈ C  , suppose there exists a basis Φ = { 1 ,  2 ,  3 , . . .,   } satisfying y = Φ, where  is a sparse vector.Making a linear measurement process to y, we have z = Fy +  = FΦ + , where F is a  ×  ( < ) random measurement matrix, FΦ is defined as dictionary, and  is additive noise.The CS theory indicates that if the matrix FΦ has optimal restricted isometric property (RIP), accurate recovery of  with high probability can be achieved by solving a convex problem as min (      θ    1 ) , where ‖ ⋅ ‖  denotes   norm and min(⋅) denotes minimization.
After clutter mitigation and interference excision, the spectrum cells of targets usually take up only a fraction of whole Doppler cells.Thus, the OTHR signal is sparse in the Doppler domain.Redefine the signal in (10) as follows: where Φ, F, and Ψ are defined as the Fourier basis of size ×, sensing matrix of size × ( < ), and redundant time-frequency dictionary, respectively. ×1 represents the integrated Doppler spectrum of targets.x  is a  × 1 vector by removing the samples contaminated by interference from x. F is the interference removing matrix of size  × , which is constructed by selecting  rows from an  ×  identity matrix with the indices of missing rows corresponding to those of the samples excised in x.First, we discuss the spectrum recovery of stably moving target.Suppose the PRF is   , which is large enough to ensure that the target's Doppler is not ambiguous.The time sequence is t = [1 : ]  /  , the Doppler resolution is Δ  =   /, and the Doppler sequence is f  : [1 : ]  ⋅ Δ  .The Fourier basis can be expressed as The redundant time-frequency dictionary can be formulated as The observation time sequence t  is a subset of the integrated time sequence.t  is a  × 1 ( < ) vector by removing the time samples contaminated by interference from t.The schematic diagram of sparse signal geometry and dictionary construction is shown in Figure 2.
Then (10) can be written as where  is redefined as the vector whose larger components correspond to the complex amplitudes of the target.It is desired that the measured data x  is not integrated; namely, within incomplete CIT, we expect to extract considerably strong targets exactly.Taking the additional noise into account, the modified convex problem with constraint to estimate the complex amplitudes is min (      θ    1 ) , where  denotes noise level.θ denotes the estimated vector of .In the case of low signal to noise ratio (SNR), we should set a high  for good estimation.In the following realdata experiment, we approximately estimate it from adjacent azimuth-range cells.
For signal reconstruction previously mentioned in ( 12)-( 16), we adopt the sinusoidal signal assumption.It is only suitable to express the echoed signal from stably moving target within a relatively short CIT.In OTHR, the PRF is very low; to achieve high Doppler resolution, the observation time is relatively long.The CIT always lasts for a few seconds or even ten seconds.At this time, the signal model in ( 13) and ( 14) no longer meets the characteristics of maneuvering moving targets.
Next, we consider the spectrum recovery of maneuvering moving target.High resolution target detection of maneuvering targets is challenging due to its motion complexity during long CIT.After range compression, the signal s for an azimuth-range cell in (3) is rewritten as follows: where  targets are assumed in the interested cell and  is the noise.The th component is given by Δ is the pulse repetition interval and Δ = 1/  .Let   = (1/( ⋅ Δ)) ⋅   and   = (1/( ⋅ Δ) 2 )  .Then the signal can be expressed as After transient interference excision, the incomplete signal in (10) can be rewritten as To clarify the expression, we define the index of pulses consisting of the valid data set as I  ×1 and the index of pulses corresponding to full data of CIT as I = [1 2 ⋅ ⋅ ⋅ ]  ×1 .I  is a subset of I.The interference excision scheme determines the index selection of I  .
Clearly, each component is a frequency modulated signal with unknown Doppler frequency, chirp rate, and complexvalued amplitude.For our convenient derivation, let  and  stand for the Doppler frequency and chirp rate, respectively.And then a partial Fourier-chirp dictionary is able to be constructed The Fourier-chirp basis can be formulated as e (, ) =   ⊙   , where "⊙" denotes Hadamard product,  = exp(−(2/)),  = exp(−(/ 2 )Δ), and Δ is the grid step of chirp rate.
Setting the Doppler series corresponds to k = [1 : ]  ×1 and the chirp rate extends to the vector l = [−/2 + 1 : /2]  ×1 ⋅ Δ, where  is supposed to be an even integer and set great enough that l includes the chirp rates of all signal components in (20).Then, the signal of a certain range cell can be rewritten as Our goal is to reconstruct the unknown vector  based on the definite partial Fourier-chirp dictionary E and the incomplete data x  .The partial Fourier-chirp dictionary E is deterministic but not parametric.So we can construct the dictionary in advance.The key step is to estimate the optimal parameters   ,   , and   for each signal component, by which the full spectrum of targets can be reconstructed straightforwardly.In OTHR, the spectrum cells of targets usually take up only a fraction of whole Doppler cells.Thus, the signal can be regarded as sparse in the Doppler domain.Such sparseness is utilized to recover the high-dimensional vector  by a  1 -norm optimization The noise term  = ‖‖ 2 can be estimated from adjacent azimuth-range cells.One can see that the optimization is composed of two different terms: the  2 -norm constraint that preserves the data fidelity of the solution and the  1 -norm optimization that imposes most elements to be small with a few large ones in accordance with the sparse feature.The optimization problem can be efficiently solved by some approaches [31][32][33][34].However, these reconstruction algorithms are complex and computationally intensive.Thus, the software implementation of these algorithms is extremely slow and power consuming.Orthogonal matching pursuit (OMP) was proposed by Tropp and Gilbert to be an efficient and reliable reconstruction algorithm [35].OMP is less complex and it is a greedy algorithm that finds the closely correlated values in each iteration.The very large scale integration (VLSI) and other hardware implementations can be found in [36][37][38].OMP reconstruction can be divided into two main stages: optimization which finds the closely correlated values and least square problem.For large scale dictionary, the implementation of correlation is time consuming since it often requires a large number of matrix multiplications.
Based on the procedure of standard orthogonal matching pursuit (OMP) algorithm, an improved solver named IOMP International Journal of Antennas and Propagation  is developed to the optimization problem (24), which implements the correlation by fast Fourier transform (FFT).The signal component estimation is detailed as described below.
Suppose we need to estimate the first signal component of x  , which is denoted as e 1 .Its Doppler frequency  1 and chirp rate  1 are achieved when the inner product of x  and the basis in E reach their maximum Intuitively, decreasing the Doppler and chirp rate grid steps will increase the resolution of recovered velocity and acceleration of target.However, if we decrease the grid steps, the dimension of E increases and leads to computational inefficiency.On the other hand, it may increase the correlation between the columns of dictionary E and degrade the recovery performance.Lowering the correlation between adjacent columns of dictionary is necessary for accurate signal recovery in CS.For a discrete Fourier-chirp dictionary, the Doppler and chirp rate steps should be selected appropriately [39,40].Then, we apply FFT instead of matrixvector multiplication directly to enhance efficiency in inner product computation in (25).Denote the inner product matrix corresponding to the dictionary and the measurement x  as Apparently, D is a  ×  matrix.Instead of calculating its element one by one by vector dot product, we may obtain a column through the following fast Fourier transform computation: where x  and   are  × 1 vectors based on x  and   , with the missing data zero-padded.By seeking the maximum element in D, we can get the estimation of ⟨ 1 ,  1 ⟩.Then, we may achieve the signal estimation as follows: Then, the residual is calculated by r 1 = x  − x  e 1 .
The procedure for IOMP reconstruction can be summarized as follows.
(2) Find the ⟨  ,   ⟩ by FFT calculation and maximization.Consider where max 2 represents finding the maximum value of a matrix.From the index   , we can determine the signal parameter ⟨  ,   ⟩.
( (4) Perform modified Gram-Schmidt orthogonalization [35,38] using the e(  ,   ) and E −1 in order to determine q.Calculate the new residual according to (5) Increment  and return to step (2) until the energy of residual signal is below the noise level.
Our goal is to recover the Doppler frequency, chirp rate, and complex-valued amplitude of each target.In (31), the reconstructed θ represents the complex-valued amplitude of targets.From index set Λ and signal set E, we can deduce each target's signal with certain Doppler frequency and chirp rate.Assume the procedure is implemented by  iterations.
After the signal decomposition into the partial Fourier-chirp dictionary, with the above estimated parameters the signal can be reconstructed through the following formulations: By the sparse decomposition under the Fourier-chirp dictionary, the amplitude, Doppler frequency, and chirp rate of each scattering center are achieved simultaneously.
International Journal of Antennas and Propagation

Results and Discussion
The proposed algorithm was evaluated using the experimental data from a trial sky wave OTHR.The PRF is   = 100 Hz and a CIT contains  = 1024 pulses.Herein, 1024 pulses with coherent time [−5.12,5.12] s are included in the data set.There are transient interference and a stably moving target in the interested azimuth-range cell.In order to validate the reconstruction performance of maneuvering moving target, we add a simulated maneuvering target into the real measured experimental data.The maneuvering target is constructed according to (3).The Doppler frequency, chirp rate, and SNR of it are set as  2 = 10.3Hz,  2 = 0.52 Hz/s, and 0 dB (before Doppler coherent integration), respectively.The real and imaginary part of the raw echo data are shown in Figure 3(a).In time domain, we can see that the targets are drowned in clutter and there exists a transient interference in the pulse indexes from 650 to 750. Figure 3 The number of reference cells (the secondary data) is  = 32 (4 azimuth cells × 8 range cells).The degree of freedom for segmental SP processing is  = 16 and the CIT is divided into  = / = 64 segments for clutter covariance matrix estimation.Figure 4(a) depicts the distribution of the eigenvalues of  covariance matrices after sorting; we can see that the two maximums of eigenvalues correspond to the clutter energy and  is set as 2. Figure 4(b) shows the comparison of the Doppler spectrum of the testing azimuth-range cell before and after segmental SP.It can be noted that the clutter is suppressed but most transient interference is retained by segmental SP.The time domain waveform of the processed signal and the excision threshold is shown in Figure 4(c).The transient interference shows up within the time interval 6.5-7.5 s.The scale factor  is set as 10 dB.After interference excision, we get the incomplete data x  ×1 ,  = 900.The dictionary E is constructed according to (21)- (22).The Doppler resolution is Δ  =   / = 0.098 Hz and the chirp rate grid is set as Δ  = 0.25 Hz/s.The length of chirp rate sequence is set as  = 8.Thus, the dimension of partial Fourier-chirp dictionary is  × ( ⋅ ) = 900 × 8192.Then, we can use OMP algorithm to solve the optimization problem in (24).The implementation of correlation is time consuming and resource demanding since it often requires matrix-vector multiplications, especially for such a large scale dictionary.According to (26)-( 27), we apply FFT instead of matrix-vector multiplication directly to enhance efficiency in inner product computation.In each iteration of the IOMP, we get the index   and determine the signal parameter ⟨  ,   ⟩.Then, we can obtain the Doppler frequency and chirp rate exactly and recover each signal according to (32).In this experiment, the Doppler frequency and chirp rate of two targets we estimated are ⟨−210, 0⟩ and ⟨105, 2⟩, respectively.Figure 5(a) shows the Doppler spectrum of FFT results by setting the contaminated samples to zero. Figure 5(b) depicts the comparison of the Doppler spectrum by convex optimization and IOMP optimization.It is shown that the spectrum of target 2 is spreading because of target's maneuverability.What is more, the granting sidelobes are relatively high because of data missing when the contaminated samples are set to zero.CS reconstruction can significantly focus the target's spectrum and reduce the granting sidelobes, which is helpful for further analysis of detection of small moving targets.More importantly, using CS reconstruction, we can estimate the Doppler frequency and chirp rate simultaneously.In this way, we can obtain the focused spectrum for maneuvering target, as shown in Figure 5(b) clearly.Both convex optimization and IOMP can reconstruct the spectrum of stably and maneuvering moving targets exactly.Compared to convex optimization, the IOMP significantly improves the computational efficiency, with slight error in reconstruction of sidelobes and background noise.

Conclusion
In this paper, we focus on the transient interference mitigation and full spectrum reconstruction of maneuvering targets for OTHR.The segmental subspace projection and adaptive threshold are used for transient interference excision.And integrated Doppler spectrum of maneuvering targets is reconstructed based on compressive sensing.Instead of directly recovering the stably moving target from incomplete measurements by using a partial Fourier dictionary, in the proposed scheme, a redundant Fourier-chirp dictionary is constructed for maneuvering moving targets.According to the characteristics of the dictionary, an improved orthogonal matching pursuit algorithm is proposed to solve the sparse optimization efficiently.The data experiments validate the effectiveness of the proposed approaches for the OTHR signal reconstruction of maneuvering targets.

Figure 2 :
Figure 2: OTHR sparse signal geometry and dictionary construction.

Figure 3 :
Figure 3: The raw echoed data of OTHR.
Transient interference excision in time domain

Figure 5 :
Figure 5: CS reconstruction of stably and maneuvering targets.
(b) shows the distribution of targets, clutter, and transient interference in frequency domain with all pulses by hamming windowed fast Fourier transform.It is shown that the clutter spectrum band is about [−2, 3] Hz.And the transient interference spans a Doppler frequency band from −8 Hz to −2 Hz.The energy of interference is stronger than that of target signals but weaker than that of clutter.There are 2 targets in Figure 3(b): target 1 is a stably moving target with Doppler frequency of −20.5 Hz and target 2 is a maneuvering target whose Doppler spectrum is broaden and shifted.Due to the maneuverability of target, Fourier transform did not properly account for the motion and it results in energy divergence and the spectrum defocus for maneuvering moving target.