Parameter Estimation for Precession Cone-Shaped Targets Based on Range–Frequency–Time Radar Data Cube

Abstract

A radar echo signal received from a cone-shaped target with precession contains micro-Doppler (m-D) information from different effective scattering centers. By taking full advantage of the m-D information, this paper proposes a parameter estimation algorithm for precession cone-shaped targets based on the range–frequency–time radar data cube (RDC). We build scattering center models of precession cone-shaped targets with the occlusion effect. The Binary Mask method is first utilized to obtain a high-resolution range-Doppler (RD) sequence. On this basis, the range–frequency–time RDC can be extracted from the RD sequence. In order to approach the actual case, we discuss the parameter estimation algorithm under different radar lines-of-sight (LOS). The most attractive attribute of this algorithm is that it can conduct in-depth research on m-D parameter estimation from a three-dimensional (3D) domain. Finally, the experimental results illustrate the effectiveness of the proposed method.

1. Introduction

Generally, micro-motion is defined as the tiny movement of an object in addition to the movement of the main body, which can more accurately describe the movement characteristics of the object [1]. For example, the flutter of aircraft wings, the clap of human hands, the precession of the space target, and so on. This additional motion of targets can give rise to the modulation effect of echo signals, which is referred to as the micro-Doppler (m-D) effect [2]. M-D effect is widespread and difficult to imitate. Therefore, it is given extensive attention from researchers in the field of target recognition [3,4,5].

In the early study stage, Chen et al. realized the unification of the cone-shaped target micro-motion mode [6]. Usually, affected by spin stabilizers, the micro-motion of a cone-shaped target will be expressed in the form of precession [7]. Then, Hongwei Gao et al. used the engineering approximation method to evaluate the effective scattering centers [8]. Furthermore, XiaoFeng Ai et al. demonstrated that the scattering center of a cone-shaped target is divided into localized scattering centers (LSCs) and sliding-type scattering centers generated by edge diffraction (SSCE) using electromagnetic simulation [9]. In addition, A. R. Persico et al. confirmed the occlusion effect at certain radar line-of-sight (LOS) angles [10]. Therefore, the radar echo signal of the precession cone-shaped target contains the m-D information from different effective scattering centers with the occlusion effect, thus making accurate parameter estimation difficult.

In radar applications, target feature extraction and parameter estimation are often achieved through different radar images. Therefore, opportune radar images are inevitably used in order to efficiently represent the m-D phenomena. Micro-motion target parameter estimation based on two-dimensional (2D) radar images is widely considered in the literature. Xueru Bai et al. used high-resolution imaging to reconstruct precession cone-shaped target scattering centers; additionally, micro-motion parameters were obtained at the same time [11]. In-Oh Choi et al. established an efficient framework for cone-shaped target parameters estimation while saving computational resources [12]. Yu Zhou et al. investigated a time–frequency (TF) curve extraction novel method called CSRDI-MGPTF and then discussed a parameter estimation method under the occlusion effect [13]. Nannan Zhu et al. used the phase-derived range (PDR) method based on a high-resolution range profile (HRRP) to realize parameter estimation under low SNR [14]. However, the inverse synthetic aperture radar (ISAR) image ignores the time-varying characteristics of scattering centers on micro-motion targets, the TF spectrogram ignores the range information, and the HRRP ignores the Doppler information. In short, using these methods, it is difficult to make full use of target echo information. Consequently, we must construct a new form of radar signal expression to represent target m-D characteristics overall.

Compared to 2D radar images, range–frequency–time joint-variable representations are more effective in target parameter estimation because they offer the possibility of using comprehensive m-D information. Several multidimensional processing techniques have been developed over the past few years [15,16,17,18,19]. For example, He Y et al. proposed a novel radar signal concept called range-Doppler surface (RDS), which can contain range, frequency, and time information [17]. S. Z. Gurbuz et al. summarized the current method of human micro-motion recognition and presented the range–time–frequency radar data cube (RDC) recognition method, but the research lacked a qualitative analysis of the model [18]. Baris Erol developed a boosting scheme using RDC-based processing to increase human classification performance [19]. However, the literature [17,18,19] only exploits the spatial shape of the range–frequency–time RDC and does not delve into its intrinsic mechanism and the rich signal characteristics it contains.

The key contribution of this paper is to suggest a new parameter estimation method based on range–frequency–time RDC. Compared with available methods, the proposed method is capable of: (a) making full use of the range, frequency, and time information; (b) realizing the distinction between LSC and SSCE in range–frequency–time RDC; (c) estimating the m-D parameters and structure parameters even under the occlusion effect.

The paper is organized as follows. In Section 2, precession cone-shaped target model and m-D characteristics analysis are introduced. In Section 3, the range–frequency–time RDC construct method is depicted. Section 4 introduces the proposed parameter estimation method for precession cone-shaped targets. Simulation results via electromagnetic computation verify the effectiveness of the analyses in Section 5. Conclusions and future discussions are drawn in the last section.

2. Cone-Shaped Target Micro-Doppler Characteristic Analysis

2.1. Echo Signal and Cone-Shaped Target Precession Model

Assuming the target’s translation motion is precisely compensated, the target echo can be described as:

$$s\left(t\right)={\displaystyle \sum }_k{\sigma }_k{e}^{j2\pi {\displaystyle \int }{f}_k\left(t\right)dt}$$

(1)

where ${\sigma }_k$ and $f_k\left(t\right)$ denote the scattering coefficient and instantaneous Doppler frequency of the $k\mathrm{th}$ scattering point, respectively.

Affected by spin and lateral disturbance, the cone-shaped target finally moves in the form of precession. The instantaneous Doppler frequency $f_k\left(t\right)$ in Equation (1) can be denoted as the m-D modulation of the scattering point.

As shown in Figure 1, the system $\mathit{O-XYZ}$ is the right-hand coordinate system. To facilitate analysis, the mass of the cone-shaped target is set at the origin of coordinate system. ${\omega }_c$ represents the cone rotation angular velocity, ${\omega }_s$ represents the spin angular velocity, and the elevation and azimuth angles of the target symmetry axis are $\gamma$ and $\alpha$, respectively. The angle between the LOS and the $O-Z$ axis is ${\theta }_0$.

According to the theory of scattering center, Point A is fixed on the top of the cone-shaped target, which can be assumed as a localized scattering center (LSC).

Thus, the micro-motion range and m-D frequency of Point A can be rewritten as follows [20]:

$$\{\begin{array}{l}{f}_A=\frac{2{\omega }_c{h}_1\mathit{sin}\gamma \mathit{sin}{\theta }_0\mathit{cos}\left({\omega }_{c}t+\alpha \right)}{\lambda }\\ r_A=h_{1}F\left(t\right)\end{array}$$

(2)

where $\lambda$ is the wavelength; $h_1$ is the distance from the top of the cone to the center of mass; $t$ is the time variable;$F\left(t\right)=\mathit{cos}\gamma \mathit{cos}{\theta }_0+\mathit{sin}\gamma \mathit{sin}{\theta }_0\mathit{sin}\left({\omega }_{c}t+\alpha \right)$.

In Equation (2), the m-D of the localized scattering center in the cone-shaped target conforms to the law of sinusoidal modulation.

Points B and C are located at the edge of the cone-shaped target, which can be considered as sliding-type scattering centers generated by edge diffraction (SSCE). Thus, the micro-motion range and m-D frequency of Point B and Point C can be rewritten as follows [20]:

$$\{\begin{array}{l}{f}_B=-\frac{2{\omega }_c\mathit{sin}\gamma \mathit{sin}{\theta }_0\mathit{cos}\left({\omega }_{c}t+\alpha \right)}{\lambda }\left\{h_2+d\frac{F\left(t\right)}{\sqrt{1-F^2\left(t\right)}}\right\}\\ r_B=-h_{2}F\left(t\right)+r\sqrt{1-F^2\left(t\right)}\end{array}$$

(3)

$$\{\begin{array}{l}{f}_C=-\frac{2{\omega }_c\mathit{sin}\gamma \mathit{sin}{\theta }_0\mathit{cos}\left({\omega }_{c}t+\alpha \right)}{\lambda }\left\{h_2-d\frac{F\left(t\right)}{\sqrt{1-F^2\left(t\right)}}\right\}\\ r_C=-h_{2}F\left(t\right)-r\sqrt{1-F^2\left(t\right)}\end{array}$$

(4)

where $h_2$ is the distance from the center of mass to the center of the bottom, and $d$ is the radius of the target bottom.

According to Equations (3) and (4), the m-D of sliding-type scattering centers in cone-shaped targets does not conform to the sinusoidal modulation law. In addition, we find that the m-D frequency of Point B and C is only different in the sign of $d\frac{F\left(t\right)}{\sqrt{1-F^2\left(t\right)}}$, and the micro-motion range of Point B and C is only different in the sign of $d\sqrt{1-F^2\left(t\right)}$. Based on this, we implement subsequent parameter estimation.

2.2. Occlusion Effect

Since not all the scattering points appear in radar observation time, the signal in Equation (1) turns into:

$$s_{occlusion}\left(t\right)={\displaystyle \sum }_k{\sigma }_k{h}_k\left(t\right)e^{j2\pi {\displaystyle \int }{f}_k\left(t\right)dt}$$

(5)

where $h_k$ is the occlusion function of the $k\mathrm{th}$ scattering point.

The occlusion function $h_k\left(t\right)$ depends on the ${\theta }_0$ and semi-angle ${\phi }_0$ in Figure 1. As shown in

Table 1. The value of occlusion function for cone-shaped target scattering points.

	${\mathit{h}}_{\mathit{A}}\left(\mathit{t}\right)$	${\mathit{h}}_{\mathit{B}}\left(\mathit{t}\right)$	${\mathit{h}}_{\mathit{C}}\left(\mathit{t}\right)$
${\theta }_0<{\phi }_0-\gamma$	1	1	1
${\phi }_0-\gamma \le {\theta }_0<\frac{\pi }{2}-{\phi }_0-\gamma$	1	1	0
$\frac{\pi }{2}-{\phi }_0-\gamma \le {\theta }_0<\frac{\pi }{2}-\gamma$	1	1	0
$\frac{\pi }{2}-\gamma \le {\theta }_0<\pi -{\phi }_0-\gamma$	1	1	1
$\pi -{\phi }_0-\gamma \le {\theta }_0<\pi -\gamma$	0	1	1

, the value of occlusion function is 0 or 1. For the localized scattering center A, the occlusion function $h_A\left(t\right)=0$ for $\pi -{\phi }_0-\gamma \le {\theta }_0<\pi -\gamma$, which means that occlusion occurs. For sliding-type scattering centers B, these scattering centers induce modulation in radar echo signal all the time, so the function $h_B\left(t\right)=1$ for all possibilities of $\theta$. When ${\phi }_0-\gamma \le {\theta }_0<\frac{\pi }{2}-\gamma$, sliding-type scattering center C is invisible, so the function $h_C\left(t\right)=0$. In addition, according to electromagnetic theory, the RCS of LSC is larger than those of SSCE [12].

3. Range–Frequency–Time Radar Data Cube Construction

In this paper, the first step of range–frequency–time RDC construction is to obtain the RD sequence. As shown in Figure 2, the RD sequence can be obtained by using the received signal.

Let us assume the radar emits a linear frequency modulated (LFM) signal; the received radar echo signal from the $k$ scattering point on the $n\mathrm{th}$ frame is found as follows:

$$s_{nk}(\widehat{t},t_n)={\displaystyle \sum }_k{\sigma }_{k}rect\left(\frac{\widehat{t}-\frac{2R_k\left(t_n\right)}{c}}{\tau }\right)\mathit{exp}\left(j2\pi f_c\left(t-\frac{2R_k\left(t_n\right)}{c}\right)\right)\mathit{exp}\left(j2\pi \frac{1}{2}\mu {\left(\widehat{t}-\frac{2R_k\left(t_n\right)}{c}\right)}^2\right)$$

(6)

where $\widehat{t}$ is fast time; $t_n$ is slow time; $\mu$ is frequency modulation slope; $f_c$ is the carrier frequency; $c$ is the velocity of light; $\tau$ is time duration. $R_k\left(t_n\right)$ is the distance between the $k$ scattering point and the radar at a slow time $t_n$, which can be solved using Equations (2)–(4).

It is worth noting that since the rotation speed of the space target is much faster than the ordinary turntable model, using traditional Fourier transform and range-Doppler algorithms cannot obtain high-resolution RD images [21]. Generally, RD images can be seen as a sparse 2D matrix [22]. Therefore, this paper uses the sparse reconstruction method to reconstruct the target signal. The goal of the sparse reconstruction is to improve the image resolution with a limited pulse number. In our work, we utilize one of the sparse reconstruction algorithms, called the two-dimensional gradient projection sequential order one negative exponential (2D-GP-SOONE) algorithm [23]. It is shown that the 2D-GP-SOONE algorithm based on the sparse reconstruction approach outperforms the conventional range-Doppler algorithm. We provide a brief introduction to the 2D-GP-SOONE algorithm.

Assume $\mathit{X}\in {ℝ}^{\overline{M}\times \overline{N}}$ as an unknown 2D sparse RD matrix and $\mathit{S}\in {ℝ}^{M\times N}$ as a known echo signal matrix. $\mathit{A}\in {ℝ}^{M\times \overline{M}}$ and $\mathit{B}\in {ℝ}^{N\times \overline{N}}$ represent the partial Fourier dictionary matrixes for range and azimuth compression, respectively. The model of radar echo signal can be described as:

$$\mathit{S}=\mathit{A}\mathit{X}{\mathit{B}}^T$$

(7)

where ${(•)}^T$ denotes the transpose of a matrix.

The two-dimensional non-convex functions in sequential order one negative exponential function (SOONE) can be expressed as

$$G_{\delta }\left(\mathit{X}\right)={\displaystyle \sum }_{i,j}\mathit{exp}\left(\left(-\frac{\left|{\mathit{X}}_{ij}\right|}{\delta }\right)\right)$$

(8)

where $\delta$ is an auxiliary variable.

Then, the gradient projection (GP) method is used to reconstruct RD sparse matrices. At last, a precise estimation of the RD sparse signal can be defined as follows:

$$\widehat{\mathit{X}}=\mathit{X}+{\mathit{A}}^{+}\left(\mathit{S}-\mathit{A}\mathit{X}{\mathit{B}}^T\right){\left(B^{+}\right)}^T$$

(9)

where ${(•)}^{+}$ denotes the pseudo-inverse of a matrix.

The 2D-GP-SOONE method is not discussed further, since it is not the focus of this paper. Interested readers may refer to [23] for more details.

After that, RD images are superimposed along the slow time axis to form an RD sequence.

To realize accurate scattering center extraction from the RD sequence, we propose a feature enhancement method based on Binary Mask. Furthermore, the feature enhancement method can also perform LSC data association automatically. It is worth noting that for the parameter estimation method proposed in this paper, the association for SSCE is non-essential.

We assume that $\mathit{X}\in {ℝ}^{\overline{M}\times \overline{N}}$ represents the restructuring RD image matrix, where $\overline{M}$ and $\overline{N}$ denotes range cells and frequency cells, respectively. The local statistical information around each pixel can be estimated as follows [24]:

$$E=\frac{1}{UV}{\displaystyle \sum }_{u,v\in \eta }\mathit{X}\left(u,v\right)$$

(10)

$${\sigma }^2=\frac{1}{UV}{\displaystyle \sum }_{u,v\in \eta }{\mathit{X}}^2\left(u,v\right)-{\mu }^2$$

(11)

where $\eta$ is the $U-\mathrm{by}-V$ local neighborhood in image matrix $\mathit{X}$, $E$ represent the local mean; ${\sigma }^2$ represents local variance.

After that, the pixel-wise adaptive Wiener filter can be expressed as:

$$w\left(u,v\right)=E+\frac{{\sigma }^2-{\sigma }_n{}^2}{{\sigma }^2}\left(\mathit{X}\left(u,v\right)-E\right)$$

(12)

where ${\sigma }_n{}^2$ denotes the variance of noise.

The pixel-wise adaptive Wiener filter can retain useful features of the image matrix while suppressing noise. We define the RD image matrix after the pixel-wise adaptive Wiener filter is ${\mathit{X}}_{w1}$. As described in Section 2, the RCS of LSC is larger than those of SSCE. Therefore, Binary Mask is then created to segment different intensities sections. The Binary Mask is a matrix of the same size as the matrix ${\mathit{X}}_{w1}$.

The first step of constructing a Binary Mask is to calculate thresholds of matrix ${\mathit{X}}_{w1}$ by using Otsu’s method [25]. After solving the threshold values $\left({\lambda }_{11},{\lambda }_{12},{\lambda }_{13}\right)$ from ${\mathit{X}}_{w1}$, the first Binary Mask can be expressed as:

$$\mathit{B}{\mathit{M}}_1\left(\overline{m},\overline{n}\right)=\{\begin{array}{c}1,{\mathit{X}}_{w1}\left(\overline{m},\overline{n}\right)>{\lambda }_{13}\\ 0,{\mathit{X}}_{w1}\left(\overline{m},\overline{n}\right)\le {\lambda }_{13}\end{array}$$

(13)

where $\overline{m}\in \overline{M}$, $\overline{n}\in \overline{N}$.

This implies that the elements in $\mathit{B}{\mathit{M}}_1(\overline{m},\overline{n})$ satisfying ${\mathit{X}}_{w1}(\overline{m},\overline{n})>{\lambda }_{13}$ are retained, and the rest are set to 0. The elements in $\mathit{B}{\mathit{M}}_1(\overline{m},\overline{n})$ satisfying ${\mathit{X}}_{w1}(\overline{m},\overline{n})\le {\lambda }_{12}$ indicate that the pixel point belongs to the background or noise of the image. It does not contain specific information and should be discarded. The elements in $\mathit{B}{\mathit{M}}_1(\overline{m},\overline{n})$ satisfying ${\lambda }_{12}<{\mathit{X}}_{w1}(\overline{m},\overline{n})\le {\lambda }_{13}$ indicate that the pixel point belongs to weak scattering points or strong noise. The elements in $\mathit{B}{\mathit{M}}_1(\overline{m},\overline{n})$ satisfying ${\mathit{X}}_{w1}(\overline{m},\overline{n})>{\lambda }_{13}$ indicate that the pixel point belongs to strong scattering points.

Then, the image matrix ${\mathit{X}}_{LSC}$, which only contains strong scattering point information, can be extracted by the following formula:

$${\mathit{X}}_{LSC}=\mathit{B}{\mathit{M}}_1\odot \mathit{X}$$

(14)

where $\odot$ denotes the Hadamard product.

We assume that the number of frames in the RD sequence is $\overline{K}$. Then, repeating the above operation, we can obtain the range and frequency information of strong scattering points in each frame of RD images. On the other hand, the association processing of strong scattering points is also realized.

Next, the image matrix ${\mathit{X}}_{\xi }$ can be obtained by the following equation:

$${\mathit{X}}_{\xi }={\left(\mathit{B}{\mathit{M}}_1\right)}^{-}\odot \mathit{X}$$

(15)

where ${(•)}^{-}$ denotes the logical not operation.

In Equation (15), the image matrix ${\mathit{X}}_{\xi }$ does not consist of strong scattering point signal components. We define ${\mathit{X}}_{\xi }$ after the pixel-wise adaptive Wiener filter is ${\mathit{X}}_{w2}$. After that, thresholds $\left({\lambda }_{21},{\lambda }_{22},{\lambda }_{23}\right)$ of matrix ${\mathit{X}}_{w2}$ can be obtained by using Otsu’s method. Similarly, the second Binary Mask can be expressed as $\mathit{B}{\mathit{M}}_2$.

Then, the image matrix ${\mathit{X}}_{SSCE}$, which only contains weak scattering point information, can be extracted using the following formula:

$${\mathit{X}}_{SSCE}=\mathit{B}{\mathit{M}}_2\odot \mathit{X}$$

(16)

To sum up, strong scattering point information can be extracted from image matrix ${\mathit{X}}_{LSC}$, and weak scattering point information can be extracted from image matrix ${\mathit{X}}_{SSCE}$.

Then, the coordinate position of scattering points in the RD sequence constitutes the range–frequency–time radar data cube.

4. Parameter Estimation Method

Firstly, the micro-motion period is estimated by the RD sequence. We choose $t_0$ as the initial time point and select the RD image ${\mathit{I}}_{ij}$ at $t_0$ as the reference image. The correlation coefficient $sim$ between ${\mathit{I}}_{ij}$ and other sequence images ${\mathit{S}}_{ij}$ can be expressed as:

$$sim=\frac{{{\displaystyle \sum }}_i{{\displaystyle \sum }}_i\left({\mathit{I}}_{ij}-\overline{I}\right)\left({\mathit{S}}_{ij}-\overline{S}\right)}{\sqrt{\left({{\displaystyle \sum }}_i{{\displaystyle \sum }}_j{\left({\mathit{I}}_{ij}-\overline{I}\right)}^2\right)\left({{\displaystyle \sum }}_i{{\displaystyle \sum }}_j{\left({\mathit{S}}_{ij}-\overline{S}\right)}^2\right)}}$$

(17)

where $\overline{I}$ and $\overline{S}$ represent the average of ${\mathit{I}}_{ij}$ and ${\mathit{S}}_{ij}$, respectively.

When reference image ${\mathit{I}}_{ij}$ and the contrast image ${\mathit{S}}_{ij}$ achieve the best matching, the correlation coefficient $sim$ takes the maximum value. Suppose the time of the best-matched image is $t_c$, then the estimation period of cone rotation is $\widehat{T}=t_c-t_0$ and the value of estimation cone rotation angular velocity is ${\widehat{\omega }}_c=\frac{2\pi }{\widehat{T}}$.

Then, the parameters including $h_1$, $h_2$, $d$, and precession angle $\gamma$ can be estimated using the range–frequency–time RDC. Affected by the occlusion effect, the parameter estimation method is discussed in three cases. What is more, in any case, it is only necessary to distinguish the strong scattering point from the weak scattering point; this means that the weak scattering point does not need to be associated.

When ${\theta }_0<{\phi }_0-\gamma$ or $\frac{\pi }{2}-\gamma \le {\theta }_0<\pi -{\phi }_0-\gamma$, effective scattering points $A$, $B$, and $C$ are visible. We assume that the number of frames in the RD sequence is $\overline{K}$. The range and frequency information of scattering point $A$ can be defined as $r_{LSC}\in {ℝ}^{\overline{K}\times 1}$ and $f_{LSC}\in {ℝ}^{\overline{K}\times 1}$, respectively. The range and frequency information of SSCE can be defined as $r_{SSCE}\in {ℝ}^{\overline{K}\times 2}$ and $f_{SSCE}\in {ℝ}^{\overline{K}\times 2}$, respectively.

We define $r^{*}=r_{LSC}+r_{SSCE\_1}+r_{SSCE\_2}$, $f^{*}=f_{LSC}-\left(f_{SSCE\_1}+f_{SSCE\_2}\right)$, where $r_{SSCE\_1}$ and $r_{SSCE\_2}$ represent the first and second columns of matrix $r_{SSCE}$, respectively. $f_{SSCE\_1}$ and $f_{SSCE\_2}$ represent the first and second columns of matrix $f_{SSCE}$, respectively.

We can define the 3D characteristic curve $s^{*}$ of the cone-shaped target, which can be written as

$$\{\begin{array}{l}{r}^{*}=r_{LSC}+r_{SSCE\_1}+r_{SSCE\_2}\\ =r_A+r_B+r_C=a+b\mathit{sin}\left({\omega }_c{t}^{*}+{\phi }_1\right)\\ f^{*}=f_{LSC}-\left(f_{SSCE\_1}+f_{SSCE\_2}\right)\\ =f_A-\left(f_B+f_C\right)=c\mathit{cos}\left({\omega }_c{t}^{*}+{\phi }_2\right)\\ t^{*}=\mathsf{\Delta }t\overline{k},\overline{k}=1,2,\cdots ,\overline{K}\end{array}$$

(18)

where $\mathsf{\Delta }t$ is interval time between each RD image, $a=\left(h_1-2h_2\right)\mathit{cos}\gamma \mathit{cos}{\theta }_0$, $b=\left(h_1-2h_2\right)\mathit{sin}\gamma \mathit{sin}{\theta }_0$, $c=\frac{2{\omega }_c\left(h_1+2h_2\right)\mathit{sin}\gamma \mathit{sin}{\theta }_0}{\lambda }$.

From Equation (18), it can be seen that parameters $a$, $b$, and $c$ are only affected by $h_1$, $h_2$, and precession angles $\gamma$ (${\theta }_0$ is pre-determined). At the same time, we find that $a$ represents the median value of the sinusoidal curve, and $b$ and $c$ represent the amplitude of the sinusoidal curve. Thus, the complex problem of extracting cone-shaped target parameters is transformed into a simple problem of estimating the mean value and amplitude of the sinusoidal curve.

As shown in Section 2, when precession angle $\gamma$ is estimated, the value of $\widehat{F}\left(t\right)$ can be described. The estimation of the bottom radius is as follows:

$$\widehat{d}=\frac{\left|r_{SSCE\_1}-r_{SSCE\_2}\right|}{2\sqrt{1-{\widehat{F}}^2\left(t\right)}}$$

(19)

When ${\phi }_0-\gamma \le {\theta }_0<\frac{\pi }{2}-\gamma$, which means only effective scattering point $A$ and $B$ are visible. By analyzing Equation (3), it can be seen that both the micro-motion and the m-D representation of scattering point $B$ are composed of two parts; the first part is the sinusoidal term, and the second part is the non-sinusoidal term. The sinusoidal part satisfies the following formula:

$$\{\begin{array}{l}{r}_{B-s}=-h_{2}F\left(t\right)=-\frac{h_2}{h_1}{r}_A\\ f_{B-s}=-\frac{2{\omega }_c\mathit{sin}\gamma \mathit{sin}{\theta }_0\mathit{cos}\left({\omega }_{c}t+\alpha \right)}{\lambda }{h}_2=-\frac{h_2}{h_1}{f}_A\end{array}$$

(20)

The non-sinusoidal part satisfies the following formula:

$$\{\begin{array}{l}{r}_{B-n}=r\sqrt{1-F^2\left(t\right)}\\ f_{B-n}=-\frac{2{\omega }_c\mathit{sin}\gamma \mathit{sin}{\theta }_0\mathit{cos}\left({\omega }_{c}t+\alpha \right)dF\left(t\right)}{\lambda \sqrt{1-F^2\left(t\right)}}\end{array}$$

(21)

To realize the parameter estimation under the occlusion effect, the compensation coefficient $\eta$ is defined as:

$$\eta =\frac{h_2}{h_1},\eta \in \left(0,1\right)$$

(22)

Under this case, the range and frequency information of scattering point $A$ can be defined as $r_{LSC}\in {ℝ}^{\overline{K}\times 1}$ and $f_{LSC}\in {ℝ}^{\overline{K}\times 1}$, respectively. The range and frequency information of SSCE can be defined as $r_{SSCE}\in {ℝ}^{\overline{K}\times 1}$ and $f_{SSCE}\in {ℝ}^{\overline{K}\times 1}$, respectively. The 3D characteristic curve $s^{*}$ can be rewritten as

$$\{\begin{array}{l}{r}^{*}=r_A+r_B+r_C=\left(1-2\eta \right)r_{LSC}\\ f^{*}=f_A-\left(f_B+f_C\right)=\left(1+2\eta \right)f_{LSC}\\ t^{*}=\Delta t\overline{k},\overline{k}=1,2,\cdots ,\overline{K}\end{array}$$

(23)

Because the compensation coefficient $\eta$ is unknown, we choose a suitable increasing sequence for $\widehat{\eta }=\left[{\widehat{\eta }}_1,\dots ,{\widehat{\eta }}_J\right]$. Each estimation compensation coefficient ${\widehat{\eta }}_j$ corresponds to a 3D characteristic curve $s_j^{*}$, which is given by:

$$\{\begin{array}{l}{r}_j{}^{*}=\left(1-2{\widehat{\eta }}_j\right)r_A\\ f_j{}^{*}=\left(1+2{\widehat{\eta }}_j\right)f_A\\ t_j{}^{*}=\Delta t\overline{k},\overline{k}=1,2,\cdots ,\overline{K}\end{array}$$

(24)

Then, we can obtain the estimated values of $\left({\widehat{h}}_{1j},{\widehat{h}}_{2j},{\widehat{\gamma }}_j\right)$ under different ${\widehat{\eta }}_j$. Similarly, the value of ${\widehat{F}}_j\left(t\right)$ under different ${\widehat{\eta }}_j$ can be described. The estimation of the bottom radius under different ${\widehat{\eta }}_j$ is as follows:

$${\widehat{d}}_j=\frac{\left|{\widehat{\eta }}_j{r}_{LSC}+r_{SSCE}\right|}{\sqrt{1-{\widehat{F}}_j{}^2\left(t\right)}}$$

(25)

Finally, the set of parameter $\left({\widehat{h}}_{1j},{\widehat{h}}_{2j},{\widehat{d}}_j,{\widehat{\gamma }}_j\right)$ is calculated from (24) and (25) under the compensation coefficient ${\widehat{\eta }}_j$.

By this means, the value of target micro-motion range and m-D frequency under different ${\widehat{\eta }}_j$ can be rewritten $({\widehat{r}}_{LSCj},{\widehat{f}}_{LSCj},{\widehat{r}}_{SSCEj},{\widehat{f}}_{SSCEj})$. To find the approximation of the compensation coefficient $\eta$, the normalized error can be expressed as:

$${\delta }_j=\frac{\left|{\widehat{r}}_{LSCj}-r_{LSC}\right|}{\mathit{max}\left|{\widehat{r}}_{LSCj}-r_{LSC}\right|}+\frac{\left|{\widehat{r}}_{SSCEj}-r_{SSCE}\right|}{\mathit{max}\left|{\widehat{r}}_{SSCEj}-r_{SSCE}\right|}+\frac{\left|{\widehat{f}}_{LSCj}-f_{LSC}\right|}{\mathit{max}\left|{\widehat{f}}_{LSCj}-f_{LSC}\right|}+\frac{\left|{\widehat{f}}_{SSCEj}-f_{SSCE}\right|}{\mathit{max}\left|{\widehat{f}}_{SSCEj}-f_{SSCE}\right|}$$

(26)

According to the above analysis, the value of ${\delta }_j$ corresponds to compensation coefficient ${\widehat{\eta }}_j$. When the value of ${\delta }_j$ reaches the minimum, it means that the sinusoidal part in the m-D information of scattering point B is compensated most completely. This implies that the estimation compensation coefficient ${\widehat{\eta }}_j$ is closer to the real value $\eta$. In other words, the set of parameters $({\widehat{h}}_{1j},{\widehat{h}}_{2j},{\widehat{d}}_j,{\widehat{\gamma }}_j)$ corresponding to the best matching estimation compensation coefficient ${\widehat{\eta }}_j$ is the final parameters’ estimated value.

When ${\phi }_0-\gamma \le {\theta }_0<\frac{\pi }{2}-\gamma$, which means only effective scattering points $B$ and $C$ are visible. The parameter estimation procedure under this condition is similar to case 2.

In this case, the range and frequency information of SSCE can be defined as $r_{SSCE}\in {ℝ}^{\overline{K}\times 2}$ and $f_{SSCE}\in {ℝ}^{\overline{K}\times 2}$, respectively. Assume $r_{SSCE\_1}$ and $r_{SSCE\_2}$ represent the first and second columns of matrix $r_{SSCE}$, respectively. $f_{SSCE\_1}$ and $f_{SSCE\_2}$ represent the first and second columns of matrix $f_{SSCE}$, respectively. The 3D characteristic curve $s^{*}$ can be rewritten as:

$$\{\begin{array}{l}{r}^{*}=r_A+r_B+r_C=\left(1-\frac{1}{2\eta }\right)\left(r_{SSCE\_1}+r_{SSCE\_2}\right)\\ f^{*}=f_A-\left(f_B+f_C\right)=-\left(1+\frac{1}{2\eta }\right)\left(f_{SSCE\_1}+f_{SSCE\_2}\right)\\ t^{*}=\Delta t\overline{k},\overline{k}=1,2,\cdots ,\overline{K}\end{array}$$

(27)

We choose a suitable increasing sequence for $\widehat{\eta }=\left[{\widehat{\eta }}_1,\dots ,{\widehat{\eta }}_J\right]$.

Similarly, the set of parameters $\left({\widehat{h}}_{1j},{\widehat{h}}_{2j},{\widehat{\gamma }}_j\right)$ can be estimated.

The estimation of the bottom radius under different ${\widehat{\eta }}_j$ is as follows:

$${\widehat{d}}_j=\frac{\left|r_{SSCE\_1}-r_{SSCE\_2}\right|}{2\sqrt{1-F_j{}^2\left(t\right)}}$$

(28)

Figure 3 shows the overall architecture for precession cone-shaped parameter estimation.

5. Simulation Results and Discussion

5.1. Validation of Algorithm

Below, as shown in Figure 4, simulated data using a smoothed cone-shaped CAD model is generated to prove the validity of the proposed parameter estimation method. Then, to obtain electromagnetic data, FEKO software is utilized with the Physical Optics (PO) method. In the experiments below, the parameters are set as shown in

Table 2. Parameters of EM model.

Start Frequency	10 GHz	Carrier Frequency	12 GHz
End frequency	12 GHz	Polarization	VV
Range of LOS	$0^{\circ }\sim {180}^{\circ }$	Conic node radius	0.075 m
Bottom radius	0.6 m	Height	2.8 m

.

In actual application, the received radar echo signal often contains white Gaussian noise; therefore, we add white Gaussian noise $n\left(t\right)$ into EM data $s\left(t\right)$. The intensity of noise is defined as:

$${\mathrm{SNR}=10\log }_{10}\left(\frac{E\left\{{\left|s\left(t\right)\right|}^2\right\}}{E\left\{{\left|n\left(t\right)\right|}^2\right\}}\right)$$

(29)

Next, we perform simulations under two cases with different LOS.

Case 1: ${\theta }_0={100}^{\circ }$, satisfied $\pi -{\phi }_0-\gamma \le {\theta }_0<\pi -\gamma$. The imaging result of the 2D-GP-SOONE algorithm is shown in Figure 5a.

Case 2: ${\theta }_0={30}^{\circ }$, satisfied ${\phi }_0-\gamma \le {\theta }_0<\frac{\pi }{2}-{\phi }_0-\gamma$. The imaging result of the 2D-GP-SOONE algorithm is shown in Figure 5b.

As shown in Figure 5, the RCS of LSC is larger than those of SSCE, which provides a basis for LSC data association.

In this simulation, the number of RD sequence frames is 100. We perform RD images every 0.02 s.

Figure 6 presents RD sequence similarity estimation result. As shown in Figure 6a, frame 7 is selected as the reference image, and the value of the correlation coefficient reaches the peak when the number of sequence frames is 29, 51, 74, and 96. Therefore, the estimated periods are 0.44 s, 0.44 s, 0.46 s, and 0.44 s, respectively. The average period is $\widehat{T}=0.445$ s. Accordingly, the estimation value of cone rotation angular velocities in case 1 is ${\overline{\omega }}_c$. Similarly, as shown in Figure 6b, the estimation value of cone rotation angular velocities in case 2 is ${\overline{\omega }}_c=\frac{2\pi }{\widehat{T}}={4.4944}^{\pi }\mathrm{rad}/\mathrm{s}$.

Parameter estimation method under case 1 can be expounded in Figure 7. As depicted in Figure 7a, the range–frequency–time radar data cube can be obtained from the RD sequence, while automatic association with LSC data is implemented by the Binary Mask algorithm.

As exhibited in Figure 7b, when the 3D characteristic curve $s^{*}$ is projected into the time–range domain and the time–frequency domain, the projections all take the form of Sine. The median and amplitude of the projection curve can be obtained by using the fitting tool. Therefore, the parameter estimation values of $s^{*}$ are $a=-0.0679$, $b=0.0819$, and $c=694.5732$. Accordingly, the parameter estimation values of the cone-shaped target are ${\widehat{h}}_1=1.9995 \mathrm{m}$, ${\widehat{h}}_2=0.7998 \mathrm{m}$, $\widehat{\gamma }={11.9977}^{\circ }$, and $\widehat{d}=0.5934 \mathrm{m}$.

Parameter estimation method under case 2 can be expounded in Figure 8. As shown in Figure 8a, affected by the occlusion effect, the range–frequency–time radar data cube provides LSC and one part of the SSCE information.

To achieve parameter estimation, the compensation coefficient is set as: $J=100$, ${\widehat{\eta }}_J=1$, and ${\widehat{\eta }}_1=0.01$. As shown in Figure 8b, when the compensation coefficient $\widehat{\eta }=0.40$, the value of normalized error reaches the minimum. This implies that the set of parameters at $\widehat{\eta }=0.40$ is the final parameter estimated value. Accordingly, the parameter estimation values of the cone-shaped target are ${\widehat{h}}_1=2.0395 \mathrm{m}$, ${\widehat{h}}_2=0.8066 \mathrm{m}$, $\widehat{\gamma }={11.7172}^{\circ }$, and $\widehat{d}=0.6110 \mathrm{m}$.

In both cases, the 3D reconstructed image of the targets is shown in Figure 9. Target geometric parameters under case 1 can be found in Figure 9a. Similarly, target geometric parameters under case 1 can be found in Figure 9b. This shows that the parameter estimation algorithm in this paper can achieve the high-precision 3D reconstruction of the target under the occlusion effect.

5.2. Algorithm Performance Analysis

To verify the effectiveness of the proposed compensation coefficient, the simulation under the different values of $\frac{h_2}{h_1}$ is implemented.

It can be seen from Figure 10 that the value of normalized error ${\delta }_j$ reaches the minimum when the estimation compensation coefficient satisfies ${\eta }_j=\frac{h_2}{h_1}$.

Next, to evaluate the performance of the proposed method under different SNR scenarios, we set the SNR ranges from −8 dB to 8 dB by using Equation (29).

Here, the Mean Absolute Percentage Error (MAPE) is introduced to the quantitative analysis of the accuracy of parameters estimation. MAPE can be defined as $\frac{100\%}{\dot{N}}{\displaystyle \sum }_{\dot{n}=1}^{\dot{N}}\left|\frac{\widehat{y}\left(\dot{n}\right)-y}{y}\right|$, where $\dot{N}$ is the number of Monte Carto experiments; $\widehat{y}\left(\dot{n}\right)$ represents the estimate of the parameter $y$ in the $\dot{n}-th$ experiment. Monte Carlo experiments are carried out 100 times under different SNR conditions, which means $\dot{N}=100$. Figure 11 describes the MAPE of five target parameters under different SNRs.

As shown in Figure 11, with the increase in SNR, the estimation accuracy of parameters $h_1$, $h_2$, and $d$ and precession angle $\gamma$ still increase, and the estimation accuracy of cone rotation angular velocity ${\omega }_c$ remains constant. The cone rotation angular velocity ${\omega }_c$ is estimated by the correlation coefficient $sim$ between different RD sequences. Since each RD sequence frame is affected by noise, the increase in noise will only reduce the similarity coefficient between the matched frame and the reference frame and will not affect the final estimation result. In fact, when the total time is constant, the higher the number of RD sequence frames, the higher the accuracy of cone rotation angular velocity estimation. At the same time, the single most striking observation to emerge from the data comparison occurs when the SNR < −6 dB, the MAPE of $h_1$, $h_2$, $d$, and $\gamma$ increases sharply. The accuracy of the parameter estimation method in this paper is affected by the RD imaging performance. When the SNR < −6 dB, part of the noise intensity is greater than the weak scattering point intensity, causing weak scattering point information to vanish.

Below, we hope to evaluate the performance of the proposed method under different LOS.

Figure 12 presents the variation of MAPE with LOS. We can see that the proposed method has better performance under any LOS, and the MARE of all estimated parameters is kept below 2%.

In order to verify the effectiveness of the method, the proposed parameter estimation method is compared with the methods in [11,12,13,14].

The results shown in

Table 3. MARE of estimated parameters with different methods.

SNR (dB)	Method	H	d	$\mathit{\gamma }$
8	Method in [11]	1.1456	2.3117	1.5514
	Method in [12]	0.6811	2.3064	1.1740
	Method in [13]	0.2341	1.3117	0.7542
	Method in [14]	1.4418	1.8452	1.8924
	Proposed method	0.1297	1.2335	0.7281
4	Method in [11]	2.2778	2.2447	1.9515
	Method in [12]	1.5088	2.6920	1.4463
	Method in [13]	1.2345	1.8764	2.3749
	Method in [14]	1.8864	2.1523	2.1291
	Proposed method	1.2549	1.7454	1.6949
0	Method in [11]	2.4135	3.9482	2.5326
	Method in [12]	2.6396	3.2837	2.6079
	Method in [13]	3.5632	5.5478	5.4750
	Method in [14]	2.7256	1.8448	5.3267
	Proposed method	1.6531	1.8346	2.3568
−4	Method in [11]	8.1545	7.3670	8.1754
	Method in [12]	8.0829	6.9260	7.1740
	Method in [13]	12.3467	13.5623	9.7423
	Method in [14]	7.3691	8.3445	5.8345
	Proposed method	2.9851	2.6969	3.9784
−8	Method in [11]	17.645	16.487	13.015
	Method in [12]	15.5556	13.9438	11.9201
	Method in [13]	21.4789	21.4587	23.1426
	Method in [14]	19.5678	19.3748	20.7918
	Proposed method	11.4079	7.2017	8.2215

indicate that the MAPE under different SNRs of the parameter estimation method proposed in this paper is obviously lower than methods in [11,12,13,14].

Next, we compare average calculation time to evaluate the computational complexity of different methods, as shown in

Table 4. Calculation time of different methods.

	Method in [11]	Method in [12]	Method in [13]	Method in [14]	Proposed Method
Time	21.483 s	0.864 s	42.952 s	2.317 s	1.098 s

. Our experiments are executed on an Intel i7 processor. The software platforms are MATLAB R2020a. It can be seen from Table 4 that the computing time of the proposed method in this paper is shorter than [11,12,13,14] and slightly higher than [12].

6. Conclusions

This paper proposed a precession cone-shaped target parameter estimation method, which uses a new form of comprehensive tool, called the range–frequency–time radar data cube. The method in this paper can realize parameter estimation under any LOS. Additionally, compared with traditional parameter estimation methods, it makes full use of the target range, frequency, and time information, has higher parameter estimation accuracy, and strong robustness. Finally, EM simulation verifies the advantages of our novel method.

The research lays the foundation for micro-Doppler-based fine target recognition. It is worth noting that the parameter estimation method in this paper is suitable for cone-shaped targets. Practically, the target’s structure is complex, such as streamlined structure targets, cone-cylinder targets, etc. Therefore, future work will be focused on parameter estimation with different target structures.