Robust low ‐ rank Hankel matrix recovery for skywave radar slow ‐ time samples

In skywave radar, the slow ‐ time samples received in a certain range ‐ azimuth cell are usually processed for signal analysis and target detection. Particularly, to extract the principal components, such as sea clutter and target signal, in slow ‐ time samples, the Hankel matrix structure and subspace ‐ based methods are often applied, which can also reduce the additive Gaussian noise. In the real environment, the slow ‐ time samples sometimes are seriously contaminated by the transient interference, affecting the performance of useful signal reconstruction. The robust low ‐ rank matrix recovery methods have been applied to solve this problem, where the transient interference is assumed to be sparse, following the Laplace distribution. However, the Hankel structures of useful signal and transient interference have not been fully utilized in the conventional methods, which can enhance the reconstruction performance. Here, the robust low ‐ rank Hankel matrix recovery problems is reformulated for skywave radar slow ‐ time samples and solve them with the inexact augmented Lagrange multiplier method. Additionally, the low ‐ rank Hankel matrix completion problem is discussed, where the location of the transient interference is determined. The experimental results have demonstrated the good performance of our proposed methods and also some interesting conclusions are obtained.


| INTRODUCTION
Skywave radar operates at the high-frequency (HF) band , which propagates by the refraction of skywave from the ionosphere. Thereby, it has the capability of detecting targets at thousands of kilometers away and overcoming the line-of-sight limitation of the conventional radar caused by the Earth curvature [1]. Skywave radar has gained much attention in the region of sea surface states remote sensing and overthe-horizon early warning [2,3]. In addition, due to the characteristic of long wavelengths, skywave radar also has the capability to detect stealth targets [4].
The skywave radar slow-time samples, received in a particular range-azimuth cell over the coherent processing interval (CPI), generally contain useful information about targets and sea clutter, and are often used for target detection and signal analysis, such as Doppler processing and spectrum estimation. To reduce the additive Gaussian noise and extract the principal components, for example, the target-plus-clutter echoes, in slow-time samples, the Hankel matrix structure is usually applied. It assumes that the Hankel structure of principal components in slow-time samples holds the low-rank property, and the low-rank matrix is obtained via singular value decomposition (SVD) of the Hankel structured slowtime samples [5]. This method, named as low-rank Hankel matrix recovery, has been applied in several fields, such as speech enhancement [6]. In skywave radar, the Hankel structure can also be applied for clutter suppression and ionospheric decontamination [7,8].
However, the slow-time samples sometimes are seriously corrupted by the transient interference which cannot be efficiently suppressed by the digital beamforming technique [9]. The transient interference may come from regional lighting discharges or local man-made noise sources [10]. The lightingbased transient interference has the potential to degrade radar sensitivity by tens of decibels when no compensation is applied. As a result, it is critical to mitigate the transient interference and restore the corrupted slow-time samples for signal processing and analysis, especially for low-rank Hankel matrix recovery.
Several signal processing techniques for suppressing transient interference in skywave radar have been published [9,[11][12][13]. The key elements of conventional techniques are similar. For the first time, they locate the transient interference in the slow-time domain. The corrupted slow-time samples are then set to zero or censored from the original data. Finally, the missing samples are reconstructed by interpolation methods based on linear prediction to restore clutter coherence over the slow-time samples [1]. An alternative approach for transient interference mitigation was proposed in [14]. It utilizes the adaptive time-frequency analysis technique to parameterize clutter signals and transient interference. The transient interference can then be identified and subtracted from the slowtime samples.
Instead of the conventional methods, the robust low-rank matrix recovery methods for skywave radar slow-time samples are proposed [15][16][17]. These methods do not require the preprocessing, detection and interpolation steps. It assumes that the transient interference is sparse in matrix and reformulates the low-rank matrix recovery problem as a convex optimization problem. The singular value thresholding (SVT) method [18] and the inexact augmented Lagrange multiplier (IALM) method [19] are used to solve the convex optimization problem. Although the study [15] introduces the Hankel structure constraint of the target-plus-clutter echoes in the matrix recovery framework, the above-cited studies have not fully utilized the Hankel structure of slow-time samples. In addition, they suppose that the transient interference is sparse in slow-time domain where a special model has not been given.
Here, we reformulate the robust low-rank Hankel matrix recovery problem for skywave radar slow-time samples. It is assumed that the target-plus-clutter signal is contaminated by the impulsive noise which is modelled as the mixture of Gaussian and Laplace distribution. To solve the low-rank matrix recovery problem, it is relaxed as a convex optimization problem and then the IALM method (named as the IALM-Hankel method) is applied. In addition, we also consider the problem that the location of the transient interference is predetermined, for example, the low-rank Hankel matrix completion problem.
The remainder of the study is organized as follows. A signal model with transient interference received in skywave radar is introduced in Section 2. In Section 3, the robust low-rank Hankel matrix recovery problem and the low-rank Hankel matrix completion problem for signal reconstruction are formulated. Section 4 proposes the solutions for low-rank Hankel matrix recovery and completion by the IALM-Hankel method. Section 5 shows the experimental results of low-rank Hankel matrix recovery and completion. Finally, Section 6 concludes our study.

| SIGNAL MODEL
In skywave radar, the slow-time samples z(m) received in a given range-azimuth cell within the CPI could be denoted as the sum of target-plus-clutter signal e(m), transient interference i(m) and additive Gaussian noise n(m), where m = 1, …, M denotes the slow-time index and M is the length of CPI. The vector form of received signal model is where z, e, n and i ∈ C M�1 . The target-plus-clutter echoes are composed of targets s (m), sea clutter c sea (m) and sometimes ground clutter c ground (m). Assuming Q independent targets are contained in slow-time samples, the target echoes could be expressed as where A s,q , f s,q and φ q are the amplitude, Doppler frequency shift and phase of the qth target respectively, and T is the pulse repetition interval. Sea clutter is the most powerful component in the echoes [16]. According to the Bragg scattering hypothesis, the dominant sea clutter spectrum received by skywave radar can be characterized by two first-order Bragg peaks which are symmetrically placed about 0 Hz in the Doppler domain [20,21]. The Doppler frequency shift of the positive first-order Bragg peak is f b ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi g cos θ=ðπλÞ p , where λ represents the operation wavelength, θ denotes the skywave grazing angle and g is the gravity acceleration. As a result, the sea clutter is modelled as where A r and A a are the receding and advancing amplitude of first-order sea clutter, respectively, and A h (m) is the higher order scatter arising mainly from the second-order interaction. For brevity, we ignore the influence of surface current and ionospheric contamination. Furthermore, the ideal ground clutter with zero Doppler frequency shift is denoted as: where A g and ψ are the amplitude and initial phase of ground clutter, respectively. The Gaussian noise n(m) may come from the receivers, atmosphere or galaxy. It assumes that n ∼ CNð0; σ 2 g IÞ, where I is the identity matrix and σ 2 g represents the Gaussian noise variance.
The transient interference is usually short in slow-time domain, but spread in Doppler domain. A typical transient interference model is constructed by Gaussian basis representation [14]:

-
where A i,k is the complex envelop, σ k is the variance, t k is the location in the time domain and f i,k is the location in the Doppler domain, of the kth transient interference, respectively. K i denotes the number of transient interference.

| PROBLEM FORMULATION
Let H be a linear operator which maps a vector x ∈ C M�1 into the Hankel matrix X ∈ C m 2 �m 1 , namely the Hankel operator [22], denoted as: where ρ m is the number of entries on the mth anti-diagonal of X [23]. Interestingly enough, the inverse Hankel operation is equivalent [24] to finding the least-squares estimate of a Hankel structured matrix, that is, where ‖⋅‖ F is the Frobenius norm of a matrix. Then the received data in Equation (2) can be reshaped in the Hankel matrix form:

HðzÞ ¼ HðeÞ þ HðnÞ þ HðiÞ: ð10Þ
Assume that the target-plus-clutter echoes e are composed of r superimposed sinusoids and r ≪ min(m 1 , m 2 ). Then the Hankel matrix HðeÞ has rank r, namely low-rank (compared with min (m 1 , m 2 )), when the sinusoidal frequencies are constant. However, if the frequencies are time-varying, the rank of Hankel matrix HðeÞ will be more than r. Fortunately, it is also shown that if the frequencies of sinusoids change not rapidly in the CPI, the rank of HðeÞ will be close to r [7].
Here, it is aimed to extract the principal components e from slow-time observations z. Note that the Hankel matrix HðeÞ holds the low-rank property. As a result, our problem could be formulated as a low-rank matrix recovery problem and the signal-plus-clutter signals e are obtained by inverse Hankel operation to the Hankel matrix HðeÞ.
For convenience, we first consider the case that the transient interference i is absent and it is aimed to extract the target-plus-clutter signal e from observation vector z which is corrupted by the additive complex Gaussian noise n, where ‖⋅‖ 2 is the Euclidean norm of a vector and rank(⋅) denotes the matrix rank operator. This procedure is also named as the noise reduction.
Unfortunately, it is difficult to solve the optimization problem (11) directly. A typical non-convex approach is to compute the nearest rank-r approximation via truncated SVD, which has also been extended to the robust low-rank Hankel matrix recovery problem [23,25,26]. However, the nonconvex approaches are not stable and the parameter selection is tricky. Here, we use the nuclear norm minimization as the convex relaxation of the rank minimization problem [27]. Hence, the optimization (11) is relaxed as where ‖⋅‖ * denotes the nuclear norm of a matrix. Furthermore, to eliminate the inequality term, problem (12) can also be reformulated as where α > 0 is a scale parameter that is related with the variance σ g and the size M of the additive Gaussian noise. Next, we consider the presence of transient interference. In skywave radar, the transient interference are mainly due to the lighting discharges and the durations of lighting are typically about 200-400 ms. For ship detection or remote sensing task, the CPI is generally about 10-80 s. Thereby, the vector i of transient interference could be considered as a sparse vector [15]. To describe the sparsity of the transient interference, we assume that the transient interference i(m) follows Laplace distribution, which is often used to describe impulsive noise in robust signal processing [28]. Then, the problem (13) is extended to the robust low-rank Hankel matrix recovery problem, where ‖⋅‖ 1 is denoted as the sum of the absolute values of a vector or matrix, and β > 0 is a scalar parameter, which is determined by the scale parameter σ l and the size M of the Laplace impulsive noise. It is interesting to note that the additive complex Gaussian noise and the transient interference can be assumed as the Laplace impulsive noise together from a statistical perspective. In other words, problem (14) can be simplified as In skywave radar, the location of transient interference sometimes can be predetermined by several methods [29]. As a result, problem (14) can be reshaped as a low-rank Hankel matrix completion problem, where Ω denotes the collection of indexes of the slow-time samples uncorrupted by the transient interference. P Ω ðxÞ is defined as the orthogonal projection which preserves the entries of x in Ω and sets the entries outside Ω to zeros. Compared with the conventional low-rank matrix recovery and completion problems, the study [30] states that under mild incoherence conditions, the Hankel structure allows perfect recovery as soon as the number of observation samples exceeds the order of r log 4 n and is stable against bounded noise. Moreover, it admits perfect signal recovery from the order of r 2 log 3 n even when a constant proportion of the samples are corrupted with arbitrary magnitudes. In the following section, we will discuss how to solve the above convex optimization problems.

| ALGORITHMS
At the beginning, the soft-thresholding (shrinkage) operator is introduced: where x ∈ C and τ > 0. This operator can be extended to vectors and matrices by applying it element-wise. For matrices X ∈ C m 2 �m 1 and W ∈ C m 2 �m 1 , the soft-thresholding operator is often applied to solve the following optimization problems: where UΣV H is the SVD of W, and S τ ðWÞ ¼ arg min In addition, Equation (19) can also be extended to vector form, that is There are a lot of methods proposed to solve the convex optimization problems in Section 3 [18,19,31,32]. Among these methods, the IALM method has good performance and computation efficiency. As a result, we choose the IALM method to solve the above convex optimizations.
Note that problems (13) and (15) are the special case of Equation (14). Hence, it is better to first discuss the derivation process for Equation (14). The augmented Lagrangian function for problem (14) is: Lðe; n; i; y; μÞ ¼ HðeÞ where y ∈ C M�1 is the Lagrange multiplier, R⟨⋅; ⋅⟩ denotes the real part of complex inner product and μ > 0 is a penalty parameter.
Using the idea in the study [19], the iterative scheme of the IALM method for Equation (14) is: where k is the iteration number. Next, we solve the equation in (22) separately. Using the following formula: where a; b ∈ C M�1 , then: Similarly, for n and i, we have: n ¼ arg min n Lðe; n; i; y; μÞ and i ¼ arg min i Lðe; n; i; y; μÞ The above results produce the iterative sequences for problem (14): To be distinguished, we name our proposed method in this section as the IALM-Hankel method. Accordingly, the iterative procedures for Problems (13) and (15) can be obtained by separately ignoring the transient interference term i and the Gaussian noise term n.
Note that the augmented Lagrangian function for problem (16) is: Lðe; n; i; y; μÞ ¼ HðeÞ Hence, Lðe; n; i; y; μÞ where � Ω represents the indexes outside Ω. As a result, the iterations for problem (16) are Here, we discuss the implementation details about the above method. The initial value for n , i and y are set to be 0. For calculation simplicity, we fix the value of μ, that is μ k = μ 0 . In the following section, it will be found that the value of μ will not change the convergence result except the convergence rate. In addition, the value of α and β are selected empirically, which will influence the reconstruction accuracy of e. To terminate the iterations, we use the convergence rate as the stopping criteria, that is: where ϵ > 0 is a scale parameter. Here, a typical value for ϵ, 10 −4 , is chosen. If the iteration number exceeds the preset value K = 20, it will also terminate the iteration procedure.

| EXPERIMENTAL RESULTS
In this section, several experiments are conducted to demonstrate the feasibility of the proposed methods for skywave radar slow-time samples. The CPI contains M = 512 pulses and the pulse repetition interval is T = 40 ms. In addition, m 1 = 257 and m 2 = 256. The simulated target-plus-clutter signal is obtained via truncated SVD for slow-time samples. Signal-to-noise ratio is set as 30 dB. The reconstruction performance is evaluated by the mean square error (MSE): where b x ∈ C M�1 represents the reconstructed signal, and the original target-plus-clutter signal or the corrupted signal is denoted as x ∈ C M�1 . If x is the original target-plus-clutter signal, the MSE describes the difference between the original signal and the reconstructed signal. If x is the corrupted signal, it exhibits the bias of reconstructed signal against the corrupted data. For the first time, we consider the influence of parameter selection on noise reduction. It assumes that the original target-plus-clutter signal is corrupted by the complex Gaussian noise and the IALM-Hankel method for problem (13) is applied. In Figure 1, the scale parameter α is fixed as 10 −5 and we vary the penalty parameter μ, for example, μ = 10 −4 , μ = 10 −5 and μ = 10 −6 . According to Figure 1, it is easy to be concluded that with the increase of iteration number, the reconstructed signal is getting closer to the original signal and is stable with a certain deviation which is independent of the value of μ. However, the penalty parameter μ determines the convergence ratio of the IALM-Hankel method. As a result, it is better to vary the value of μ iteratively for fast algorithms.
On the other side, we fix the value of μ and vary α, for example, α = 10 −4 , α = 10 −5 and α = 10 −6 , as depicted in Figure 2. It is inferred that the value of scale parameter α determines the signal reconstruction performance. Particularly, the scale parameter α is related to the variance of Gaussian noise. As shown in Figure 2(a) and (b), the lower the value of α is than the reasonable value (α = 10 −5 ), the larger the MSE is for the reconstructed signal against the original signal. Interestingly enough, if the value of scale parameter is higher than the reasonable value, the reconstructed signal will get closer to the contaminated samples. As a result, it is important to choose the appropriate scale parameter for good reconstruction. Similar conclusions can also be extended into Problems (14)- (16).
Next, we discuss the case that the transient interference is present in slow-time samples and the IALM-Hankel methods for Gaussian noise and Laplace noise are used separately. For convenience, the slow-time samples corrupted by the transient interference from slow-time index 274-288 are obtained from the real skywave radar data. The parameters for IALM-Hankel methods are displayed in Figure 3. It can be seen that the IALM-Hankel method for Laplace noise can efficiently suppress the transient interference and smooth the samples corrupted by Gaussian noise while the method for Gaussian noise is not. Furthermore, Figure 3(c) shows that the IALM-Hankel method can also reduce the noise level in power spectrum domain which can raise the probability of target detection.
Here, the Hankel structure of signal is utilized to enhance the reconstruction performance. Using the same data in Figure 3, Figure 4 shows the results of MSE against the corrupted samples, which correspond to the IALM method and IALM-Hankel method, respectively. Although the study [30]  has proved that the Hankel structure of signal will enhance the reconstruction performance, it seems that the difference between the IALM methods for real data processing is very small. A possible explanation is that the inverse Hankel operation in Equation (24) will destroy the low-rank property of signal Hankel matrix. As a result, from the perspective of engineering application, there is no need to consider the Hankel structure of signal in iteration procedures.
In Section 3, we respectively model the impulsive noise (transient interference and Gaussian noise) as the Laplace distribution and the mixture of Gaussian and Laplace distribution. To compare the performance of these two kinds of distribution model, the IALM-Hankel methods for problems  (14) and (15) are applied for a section of real data. The real data is contaminated by the transient interference which can be divided into three sections. The slow time indexes of the first and the third section are from 100 to 150 and from 430 to 450, respectively. The second section is from 300 to 320. Particularly, the power of the transient interference in the middle section is stronger than that of the other data sections, as shown in Figure 5(a). The parameter settings for problem (14) are the same as that in Figure 3 while for problem (15) are μ = 5 � 10 −7 and β = 9 � 10 −2 . According to Figure 5(a) and (b), the signal reconstruction performance of mixture model for impulsive noise-contaminated samples is similar to that of the Laplace distribution model. Although the mixture model can describe more complicated impulsive noise distribution, it seems that there is no need to complicate the robust low-rank Hankel matrix recovery problem, due to the additional parameter α. Sometimes, the location of transient interference can be predetermined and the missing data reconstruction method can be used. Instead of the conventional method, we use the matrix completion to reconstruct the missing data, as illustrated in Figure 6. It can be seen that the method used in our study has a good performance in signal reconstruction, similar to the result in Figure 5(a), where the location of transient interference is unknown. For conventional methods, the corrupted slow-time samples are discarded as missing (zeroing the contaminated samples). However, there may be useful information in the corrupted samples as well, particularly when the amplitude of transient interference is small, which can improve the quality of matrix completion compared with the 'zeroing'case. Our method for problem (16) makes full use of the information of the transient interference.
Finally, the Range-Doppler (RD) maps for real data are shown in Figure 7. In skywave radar, the RD map on which the constant false alarm ratio detection is conducted plays a key role. Figure 7(a) displays the result of conventional method, that is, directly applying fast Fourier transform (FFT) with Hanning window for slow-time samples in every range bin. It can be seen that due to the presence of transient interference, the level of noise power is raised which causes the missed detection of targets (targets A and B). Next, we use the subspace-based methods to extract principal components in slow-time samples and apply the FFT-based method to estimate its power spectrum. The result of the method for Gaussian noise is illustrated in Figure 7(b). Although this method can reduce Gaussian noise level, the influence of transient interference still remains. Instead, the methods for Laplace noise, depicted in Figure 7(c) and (d), behave robustly in transient interference environment. Especially, the method used in Figure 7(c) is the IALM-Hankel method while in Figure 7(d) is the SVT-Hankel method. The penalty parameter for the SVT-Hankel method is the same as that of the IALM-Hankel method, that is μ = 5 � 10 −7 , and the step size of it is set as δ = 0.3. It can be seen that although both of the two methods can suppress the transient interference, the IALM-Hankel method has a better performance.

| CONCLUSION
Robust low-rank Hankel matrix recovery for skywave radar slow-time samples has been studied here. To extract the target-plus-clutter echoes in slow-time samples, the Hankel matrix structure is utilized and it assumes that the Hankel structured principal components hold the low-rank property, which deduces the low-rank Hankel matrix recovery problem. Considering that the slow-time samples in skywave radar are usually corrupted by the Gaussian noise and sometimes the transient interference, we reconstruct the robust low-rank Hankel matrix recovery problem, where the transient interference is modelled as the Laplace distribution. However, the robust low-rank Hankel matrix recovery problem cannot be solved directly. To overcome this problem, we replace the rank minimization with the nuclear norm minimization. As a result, the non-convex problem is relaxed as the convex problem which can be solved easily. Noting that it is redundant to model the interference as the mixture of Gaussian and Laplace distribution, we also simplify the noise model as the Laplace distribution. In addition, the case that the location of the transient interference is predetermined is also considered and the low-rank Hankel matrix completion is discussed here. To solve the low-rank Hankel matrix recovery and completion problems, we introduce the augmented Lagrange multiplier In the experimental part, we first study the influence of parameter selection on noise reduction. It is easy to be concluded that the penalty parameter μ only determines the convergence ratio of the IALM-Hankel method while the scale parameter α is related to the signal reconstruction performance. Similar conclusion can also be extended to the scale parameter β for the transient interference. Furthermore, the results of real data processing show the efficiency of the robust low-rank Hankel matrix recovery method for suppressing the transient interference and raising the probability of target detection. In addition, our experimental results show that there is no need to consider the Hankel structure of signal in iteration procedures. Although it is better to model the interference as the mixture of Gaussian and Laplace distribution, the Laplace distribution model still holds good reconstruction performance. Finally, we use the real data corrupted by the transient interference for robust low-rank Hankel matrix recovery. Compared with the conventional method and the SVT-Hankel method, the robust IALM-Hankel method has the best performance.
There are still several problems needed to be solved in the future. For the first time, how to choose the appropriate parameter value for better reconstructing the desired signal. In addition, here, we use the nuclear norm minimization to construct the optimization problem, where a better description for the principal components may be proposed. Particularly, the distribution model of impulsive noise can be extended to the general Gaussian distribution model. All the above are worth carefully studying and deep thoughts.