1. Introduction
The echoes of High Frequency Surface Wave Radar contains sea clutter, ionospheric clutter and other kinds of noises. Among them, sea clutter and ionospheric clutter have the most serious impact on targets detection, since the intensity of these two clutters. Often more than tens of dB higher than the signal. The research on HFSWR sea clutter has been very sufficient. For example, Haykin [1] has proved that sea clutter was a chaotic process in both microwave radar and high frequency radar. If the ionospheric clutter of HFSWR can be proved to be a chaotic process, then we can use the properties of the chaotic system to extract the weak target signals.
Target signal extraction in the context of chaotic noise has been a popular research area. Fei [2] studied an image encryption algorithm based on mixed chaotic dynamic systems and external keys. Uchida [3] used two microchip solid-state lasers with external modulation as chaotic signal sources and the simulation results showed that algorithm gived a good performance in the separation and denoising of mixed noisy signals in the presence of a white Gaussian noise or stationary colored noise, when analysis based on the second-order statistics was applicable [4]. Xie [5] studied blind source separation of continuous-time chaotic signals based on fast random search algorithm. A fast random search (FRS) algorithm was, therefore, proposed. Due to the fact that the great amount data of speech communications and real-time communication has been required [6], utilized the intractability of the underdetermined BSS problem to present a dual key speech encryption method. Kim [8] assumed that one had no knowledge about the governing equations of the source signals, and that the mixed signal was simply the sum of the sources. From the perspective of dynamical systems proposed a supervised learning framework that can solve this problem through an intermediate dynamical system. With the mixing fractions being known in advance or unknown, the schemes and theories for the separations of two groups of the mixed optical chaotic signals were proposed in detail, using the VCSEL-based reservoir computing (RC) systems [9]. In particular, in the field of radar signal processing, Haykin [10] was the first to introduce the chaos analysis method into the radar echo signal processing, through the analysis of the IPIX radar sea clutter echo signals, to get the conclusion that sea clutter has the characteristics of chaos, and then use artificial neural networks to study the extraction of small target signals in the background of sea clutter.
Since the stochastic characteristics of the noise will cover the chaotic feature of the ionospheric clutter signal, this paper firstly removed part of the white noise by applying empirical mode decomposition (EMD) method, and then reconstructed the ionospheric clutter through the geometrical properties of the chaotic attractor, and projected the chaotic signal into the tangent space of the attractor manifolds in the local region, whereas the weak harmonic signals (i.e., the target signals) were retained in the complementary space of the local tangent space of the attractor, which in turn achieved the extraction of the weak target signals.
2. Pre-Processing of Experimental Data
The data in this paper were from the High Frequency Surface Wave Radar (HFGWR) established in Weihai, Shandong Province, China. Although the HFSWR was designed to propagate around the sea surface, but there will still be part of the energy irradiated to the ionosphere, which is emitted by the ionosphere and then received by the receiving antenna to form ionospheric clutter eventually.
First, we determined whether the ionospheric clutter signal possessed chaotic characteristics. In the chaos recognition method, it can be comprehensively determined from both qualitative and quantitative aspects. From the viewpoint of causes of ionospheric clutter, the radar echo signal shows randomness and certainty at the same time, and this complex behavior is due to the fact that many kinds of nonlinear motions in the ionosphere are intertwined with each other; From the point of view of quantitative analysis, we estimated the Lyapunov exponent (the degree of dispersion of the attractor), Komogorov entropy (the degree of chaos) and fractal dimension (the number of attractor dimensions) of the clutter signal, and these three parameters are the invariants of attractors and can be used to characterize the chaotic properties.
Another factor that cannot be ignored is noise, which has a serious effect on chaotic signals. Therefore, in order to utilize the chaotic property of ionospheric clutter, this paper adopted an empirical mode decomposition denoising algorithm to remove part of the noise and increased the clutter-to-noise ratio of the ionospheric clutter. The noise reduction method is based on the property that the product of the white noise energy density and the corresponding averaging period tends to be a constant. Figure 1 shows the ionospheric clutter time series before and after the denoising process.
 |
Fig. 1 The time series of ionospheric clutter based on EMD denoising. |
Table 1 showed the estimation results of chaotic invariants of ionospheric clutter time series. According to the calculation results, at least two results can be obtained: firstly, based on the fact that the maximum Lyapunov exponent and the Kolmogorov entropy were larger than zero and the correlation dimension contained decimals with fractal characteristics, it can be concluded that the time series ionospheric clutter was a chaotic process; Secondly, after EMD denoising, the three kinds of chaotic invariants had realized growth, among which the increase of the maximum Lyapunov exponent proved that the chaotic feature of the processed data was enhanced.
Table 1 Basic information about database. |
 |
3. Mathematical Modeling Procedure
Now we demonstrate the chaotic feature of the HFSWR ionospheric clutter signal. Ideally, chaotic signals can be mathematically modeled using one or a set of differential equations. But the ionospheric clutter, as a real physical phenomenon in nature, the exact description of its dynamics is difficult to achieve. Therefore, this paper adopted an improved local tangent space projection algorithm with the following procedure:
Assuming that the ionospheric clutter is \(y \left( n \right)\) and the weak harmonic signal(i.e., the target signal) is \(s \left( n \right)\), the mixed signal is \(r\left(n\right) = y\left(n\right) + s\left(n\right)\), the manifold which the attractor of \(y\left(n\right)\) located in is \(\mathcal{M}\), and the manifold which the attractor of \(r\left(n\right)\) located in is \(\mathcal{M}^{\prime}\). In general, the dimension of the manifold \(\mathcal{M}^{\prime}\) is higher than \(\mathcal{M}\). Since \(s \left( n \right)\) is a deterministic harmonic signal, then the dynamical systems corresponding to \(r\left(n\right)\) and \(y \left(n\right)\) are differential homogeneous, which is the theoretical basis for extracting weak signals. Reconstruct the phase space for \(r\left(n\right)\) and \(y\left(y\right)\), respectively:
| \[\begin{equation*} \begin{aligned} & r(n)=[r(n+\tau(d-1)), r(n+\tau(d-2)), \ldots, r(n)]^T \\ & y(n)=[y(n+\tau(d-1)), y(n+\tau(d-2)), \ldots, y(n)]^T \end{aligned} \tag{1} \end{equation*}\] |
Where, \(\tau\) is the delay time, \(d\) is the reconstruction dimension. Suppose that there are \(Q\) proximity points \(r\left(n_i\right) \left( i=1,2,\dots,Q \right)\) in the neighborhood space \(\Omega_n\). It is required that \(Q\) is greater than \(d\). Then estimate the local tangent space \(T_{y\left(n\right)}M\) of the epidemic in which \(y\left(n\right)\) is located:
-
Estimate the space \(\Omega_n\) center \(\overline{r}\left(n_0\right)=\frac{1}{Q}\sum_{i=1}^{Q}r\left(n_i\right)\), generally \(\overline{r}\left(n_0\right) \neq r\left(n_0\right)\);
-
Construct transition matrix \(A\):
and singular value decomposition of A is:\[\begin{equation*} \begin{split} A= & \left[\frac{r\left(n_1\right)-\bar{r}\left(n_0\right)}{\left\|r\left(n_1\right)-\bar{r}\left(n_0\right)\right\|}, \frac{r\left(n_2\right)-\bar{r}\left(n_0\right)}{\left\|r\left(n_2\right)-\bar{r}\left(n_0\right)\right\|}, \right. \\ &\left.\ldots,\frac{r\left(n_Q\right)-\bar{r}\left(n_0\right)}{\left\|r\left(n_Q\right)-\bar{r}\left(n_0\right)\right\|}\right] \end{split} \tag{2} \end{equation*}\] \[\begin{equation*} \begin{split} A=& \left[u_1, u_2, \ldots, u_d\right] \cdot \\ &\left[\operatorname{diag}\left(\sigma_1^2, \sigma_2^2, \ldots, \sigma_Q^2\right) \mid 0_{d \times(Q-d)}\right] \cdot V^T \end{split} \tag{3} \end{equation*}\] -
Using the singular values of the matrix \(A\) obtained in the previous step can determine the local tangent space dimension \(l\), which is also the dimension of the manifold \(\mathcal{M}\): If there exists \(k\) such that \(\sigma_k^2 \gg \sigma_{k+1}^2\) and \(\sum_{i=1}^k \sigma_i^2 \approx \sum_{i=1}^d \sigma_i^2\), then \(l=k\); otherwise, set a threshold close to 1, denoted as \(v\), \(l\) is equal to the minimum \(k\) that holds the following equation:
\[\begin{equation*} \left[\sum_{i=1}^k \sigma_i^2 / \sum_{i=1}^d \sigma_i^2\right]^{1 / 2} \geq v \tag{4} \end{equation*}\] -
The coordinate point of origin of the local tangent space \(T_{y\left(n\right)}M\) is \(\overline{r}\left(n_0\right)\), and one set of orthogonal bases is \(\left[u_1,u_2,\dots,u_l\right]\)
Denote the nonlinear mappings \(F:\mathcal{M}\rightarrow \mathcal{M}\) and \(F:\mathcal{M}^\prime \rightarrow \mathcal{M}^\prime\) as the dynamical equations for \(y\left(n\right)\) and \(r\left(n\right)\), respectively. If the nonlinear mapping \(F\) corresponding to the reconstructed dynamical equations is restricted to the interior mapping of \(\mathcal{M}\), the reconstructed dynamical equations are guaranteed to come from chaotic noise. In order to satisfy such requirements, it is first necessary to project the neighboring points \(r\left(n_i\right)\left(i=1,2,\dots,Q\right)\) of \(r\left(n\right)\) and the vector after a unit time interval \(r\left(n_i+1\right)\) onto \(T_{y\left(n\right)}M\) and \(T_{y\left(n+1\right)}M\), respectively, and its projection is denoted as:
| \[\begin{equation*} \begin{gathered} y^{\prime}\left(n_i\right)=\bar{r}(n)+U_l U_l^T\left(r\left(n_i\right)-\bar{r}(n)\right) \\ y^{\prime}\left(n_i+1\right)=\bar{r}(n+1)+\tilde{U}_l \tilde{U}_l^T\left(r\left(n_i+1\right)-\bar{r}(n+1)\right) \end{gathered} \tag{5} \end{equation*}\] |
where, \(\overline{r}\left(n\right)\) and \(\overline{r}\left(n+1\right)\) are the centers of \({r}\left(n_i\right)\) and \({r}\left(n_i+1\right)\), respectively; \(U_l=[u_1,u_2,\dots,u_l]\) and \(\tilde{U}_l=[\tilde{u}_1,\tilde{u}_2,\dots,\tilde{u}_l]\) are the unit orthogonal bases of \(T_{y\left(n\right)}M\) and \(T_{y\left(n+1\right)}M\), respectively. The best estimate of the mapping \(f\) is that satisfies the following equation:
| \[\begin{equation*} y^{\prime}\left(n_i+1\right)=\hat{F}\left(y^{\prime}\left(n_i\right)\right), i=1,2, \ldots, Q \tag{6} \end{equation*}\] |
\(\hat{F}\) is the local linear form of the mapping \(F\), which is equivalent to the linear mapping \(L T_{y\left(n\right)} \rightarrow T_{y\left(n+1\right)}\), and make the following equation hold:
| \[\begin{equation*} \tilde{U}_l^T\left(r\left(n_i+1\right)-\bar{r}(n+1)\right)=L U_l^T\left(r\left(n_i\right)-\bar{r}(n)\right) \tag{7} \end{equation*}\] |
The reconstructed dynamical system equations can be obtained using the least squares method as:
| \[\begin{equation*} \begin{split} y^{\prime}(n+1)&=\hat{F}\left(y^{\prime}(n)\right)=\bar{r}(n+1)+\ldots \\ &\tilde{U}_l L U_l^T\left(y^{\prime}(n)-\bar{r}(n)\right)=A_n y^{\prime}(n)+b_n \end{split} \tag{8} \end{equation*}\] |
Where, \(A_n=\tilde{U}_l L U_l^T\) and \(\quad b_n=\bar{r}(n+1)-A_n \bar{r}(n)\). Substituting \(r\left(n\right)\) for \(y^\prime\left(n\right)\) on the right side of the equal sign gives the prediction equation:
| \[\begin{equation*} y^{\prime}(n+1)=A_n r(n)+b_n \tag{9} \end{equation*}\] |
Then take the trajectory from \(n\) to \(n+N\), a set of equations are obtained as follows:
| \[\begin{equation*} \left\{\begin{array}{c} y^{\prime}(n+1)=A_n y^{\prime}(n)+b_n \\ \vdots \\ y^{\prime}(n+N)=A_{n+N-1} y^{\prime}(n+N-1)+b_{n+N-1} \end{array}\right. \tag{10} \end{equation*}\] |
Where, \(y^{\prime}(n)=\left[y^{\prime}(n-\tau(d-1)),\ldots,y^{\prime}(n)\right]\). Equation (10) obtain \(N\times d\) equations and \(N+d\) variables. Since \(y(n)=[y(n-\tau(d-1)), y(n-\tau(d-2)), \ldots, y(n)]^T\) is the state variable of the chaotic system, which contains \(d\) independent variables. According to the deterministic mechanism of chaotic systems, it is known that the value of the chaotic signal can be uniquely determined by these \(d\) variables. Thus there are actually only \(N\) independent equations in the system of Eq. (10). To solve this equation, the generally given solution is to define the distance \(D\) between two chaotic trajectories, then find the minimum value of the distance \(D\). The disadvantage of this solution is that the equations will be in a distorted state due to the instability of the chaotic motion. In order to reduce the scale of variables and increase the computational speed, this paper adopts an improved solving method. The core algorithm is to add a set of transition equations:
| \[\begin{equation*} y^{\prime}(j)=A_{j-1} r(j-1)+b_{j-1}, j=n, \ldots, n+N \tag{11} \end{equation*}\] |
By presetting certain initial values on the desired chaotic trajectory, the solution of the equations can be calculated through iteration. The advantage of this approach is to accelerate the solution speed and ensure that the trajectory of the desired chaotic system is sufficiently close to the trajectory of the observed chaotic system. The more initial values are given, the more accurate the estimation of the chaotic trajectory will be. And the initial values used in the experiment is ten percent of the measured data in our experiment.
In summary, the algorithm on the separation of weak harmonic signals in the background of chaotic clutter in the ionosphere of HFSWR can be summarized as follows:
-
The time series of the HFSWR ionospheric clutter signal is reconstructed in phase space using the time-delay method to obtain the vector \(r(n)=[r(n-\tau(d-1)), r(n-\tau(d-2)), \ldots, r(n)]\).
-
Calculate \(r\left(k\right)\) and its neighborhood points \(\{ r(k_i) \}_{i=1}^Q\), then estimate the local tangent space \(T_{y\left(k\right)}M\) where \(y\left(k\right)\) locates in, and in particular, \(k=n,n+1,\dots,n+N\).
-
Project \(\{r\left(k_i\right)\) and \(\{r\left(k_{i+1}\right)\}\) onto \(T_{y\left(k\right)}M\) and \(T_{y\left(k+1\right)}M\), respectively. Then the least squares method is used to determine the local linear mapping \(L:T_{y\left(k\right)}M \rightarrow T_{y\left(k+1\right)}M\), which can reconstruct the chaotic dynamical system equations.
-
Estimate the chaotic trajectory \(y^\prime(k),k=n,n+1,\dots,n+N\), then the weak harmonic signal \(s\prime\left(k\right)\) can be expressed as \(s^{\prime}(k)=r(k)-y^{\prime}(k)\). Iterate \(n\) and repeat Step1 to Step4 to realize the extraction of weak harmonic signals.
4. Experimental Results
Measured data was from High Frequency Surface Wave Radar Experimental Station, Weihai, Shandong, China and the experimental date was March 20, 2023. The weak harmonic signal was a ship which was picked up on the Automatic Identification System (AIG). The ship was about 200 kilometers from the radar station, which coincided with the distance cell where the ionospheric F1-layer was located. A sequence of ionospheric clutter clutter containing a weak harmonic signal (the ship signal) was shown in Fig. 2, where we can not separate the ship signal from the ionospheric clutter.
 |
Fig. 2 Ionospheric clutter with weak harmonic signals. |
The measured data were processed according to the proposed separation and reconstruction method in Sect. 3, and the weak harmonic signals obtained after the separation from the chaotic ionospheric clutter background were shown in Fig. 3. From the Fig. 3, it can be seen that compared with the ideal situation, the separated harmonic signals were in the same phase, but there was a certain degree of attenuation, which was due to the fact that in the process of local tangent space projection, a part of the energy of harmonic signals was unavoidably projected into the local tangent space, and the energy of signals in the complementary space was reduced.
 |
Fig. 3 Extraction of weak harmonic signals. |
We compare the local tangent space projection (LTSP) method used in this paper with the classical neural network (NN) and blind source separation (BSS), where the experimental data and hardware conditions used are identical, and the comparison results are shown in Table 2. From Table 2, it can be found that the results of BSS for chaotic ionospheric clutter are not satisfactory, and the performance of NN is not stable due to the limitation of the choice of the number of layers and the number of hidden nodes, whereas the methodproposed in this paper utilizes the dynamics of ionospheric clutter and achieves better separation result.
Table 2 Comparison of algorithms performance. |
 |
5. Conclusion
In that paper, we proposed a new method for extracting weak target signals from strong chaotic noise background. This method exploited the property that the trajectory of a chaotic system in the local tangent space was inconsistent with the target signal. We first proved that the ionospheric clutter from HFSWR was a naturally occurring chaotic system by using several characteristic parameters, and then verified the method proposed in that paper using the measured data. The experimental results showed that this separation method could find the weak harmonic signal, i.e., the target signal, from the chaotic ionospheric clutter background. We will continue to optimize the algorithm in future research.
Acknowledgments
This research was funded by the National Natural Science Foundations of China (NSFC) under Grant 62061026 and Doctoral Fund Project of LongDong University under Grant XYBY202001.
References
[1] S. Haykin and H. Leung, “Chaotic model of sea clutter using a neural network,” Proc. SPIE The International Society for Optical Engineering, 1989.
CrossRef
[2] P. Fei, S.-S. Qiu, and L. Min, “An image encryption algorithm based on mixed chaotic dynamic systems and external keys,” Proc. 2005 International Conference on Communications, 2005.
CrossRef
[3] A. Uchida, M. Kuraya, S. Yoshimori, and K. Umeno, “Separation of mixed chaotic signals in microchip lasers by independent component analysis,” 2007 European Conference on Lasers and Electro-Optics, 2007.
CrossRef
[4] K. Pukenas, “Blind separation of noisy pseudoperiodic chaotic signals,” Elektronika Ir Elektrotechnika, vol.91, no.3, pp.31-34, 2009.
URL
[5] Z. Xie and J. Feng, “Blind source separation of continuous-time chaotic signals based on fast random search algorithm,” IEEE Trans. Circuits Syst. II, Exp. Briefs, vol.57, no.6, pp.461-465, 2010.
CrossRef
[6] H. Zhao, S. He, Z. Chen, and X. Zhang, “Dual key speech encryption algorithm based underdetermined BSS,” The Scientific World Journal, vol.2014, pp.1-7, 2014.
CrossRef
[7] J. He, Y. Song, P. Du, and L. Xu, “Analysis of single channel blind source separation algorithm for chaotic signals,” Mathematical Problems in Engineering, vol.2018, 9571510, 2018.
CrossRef
[8] Z. Lu, J. Kim, and D. Bassett, “Chaotic source separation solved by a tank of water through invertible generalized synchronization,” APS, 2020.
[9] D. Zhong, Y. Hu, K. Zhao, W. Deng, P. Hou, and J. Zhang, “Accurate separation of mixed high-dimension optical-chaotic signals using optical reservoir computing based on optically pumped VCSELs,” Opt. Express, vol.30, no.22, pp.39561-39581, 2022.
CrossRef
[10] B.X. Li, and S. Haykin, “Chaotic detection of small target in sea clutter,” IEEE International Conference on Acoustics, 2002.
CrossRef
