# Micromotion Feature Extraction of Space Target Based on Track-Before-Detect

Authors: Yijun Chen; Qun Zhang; Ying Luo; Tat Soon Yeo
Journal: Journal of Sensors
 (2017)

---

## Abstract

The micromotion feature of space target provides an effective approach for target recognition. The existing micromotion feature extraction is implemented after target detection and tracking; thus the radar resources need to be allocated for target detection, tracking, and feature extraction, successively. If the feature extraction can be implemented by utilizing the target detecting and tracking pulses, the radar efficiency can be improved. In this paper, by establishing a feedback loop between micromotion feature extraction and track-before-detect (TBD) of target, a novel feature extraction method for space target is proposed. The TBD technology is utilized to obtain the range-slow-time curves of target scatterers. Then, micromotion feature parameters are estimated from the acquired curve information. In return, the state transition set of TBD is updated adaptively according to these extracted feature parameters. As a result, the micromotion feature parameters of space target can be extracted concurrently with implementing the target detecting and tracking. Simulation results show the effectiveness of the proposed method.


---

## Body

## 1. Introduction

The monitoring and recognition of space targets have become more and more important since these space targets such as fragments of satellites, rocket bodies, and other space debris are hazardous to aerospace activities with the increasing outer space explorations [1]. Generally, a space target has complex micromotion such as spinning, precession, and rolling, in addition to the body translation [2]. These micromotions will induce the micro-Doppler (m-D) effect, which can be regarded as a unique signature, providing additional information for classification, recognition, and identification of the target [3]. Therefore, since the concept of micro-Doppler was introduced to the radar signal processing field [4, 5], the micromotion feature extraction of space target has drawn extensive attention of scholars [6–16].Based on the fixed scatterer model [7], the micromotion feature of space target was usually analyzed with the time-frequency distribution technique, and the micromotion feature of target was extracted with parameter transformation method [7]. However, the time-frequency methods have some limitations due to the fact that the change of micro-Doppler frequency with time is highly nonlinear [8]. A feature extraction method based on empirical-mode decomposition (EMD) was proposed to estimate the oscillation frequency of the truck’s surface in [9]. To get the higher potential payoffs from the exploitation of m-D effect in wideband radar system, the micromotion feature extraction for wideband radar based on complex image orthogonal matching pursuit decomposition is proposed in [10].However, many space targets are of conical and cone-cylinder shapes, these rotationally symmetric targets cannot be described exactly with the fixed scatterer model. Their scattering centers of the edge of the conic bottom slide with the radar line of sight (LOS) that means the scattering centers local at the two intersections of the conic bottom edge and the radar wave incident plane [11]. Based on the slide scatterer model, the mathematical expression of micro-Doppler is derived, which is validated by the darkroom data in [12, 13]. However, these works did not provide effective micromotion feature extraction methods for the slide scatterer model. Unfortunately, most of the existing micromotion feature extraction methods for the fixed scatterer model will not function properly because the micro-Doppler curve of a slide scatterer is much different from that of a fixed scatterer. Although some methods based on the complex generalized Radon transform and extended Hough transform are proposed for micromotion feature extraction based on the slide scatterer model [14, 15], the computation load of these methods is heavy, and the methods need high pulse repetition frequency (PRF) to avoid the frequency domain aliasing phenomenon. In this case, a lot of radar resources needs to be allocated for micromotion feature extraction.In fact, the existing micromotion feature extraction methods of space target are implemented after target detection and tracking. It is necessary to allocate the limited radar resources for target detecting, tracking, and feature extraction, successively. In the case of multitarget monitoring, the allocation contradiction of radar resources will be serious. To overcome this problem, we intend to seek a novel approach which can implement target detecting and tracking simultaneously when extracting the micromotion features of a space target, which can be available for both the fixed scatterer model and the slide scatterer model.In recent years, the track-before-detect (TBD) technology has shown good performance in weak target detecting and tracking. Unlike traditional techniques that declare the presence of a target at each scan, the TBD technology processes more consecutive scans jointly, and then it declares the presence of a target and its track [16–18]. Through the interscan accumulation, the TBD technology can improve the probability of target detection. The idea of establishing a feedback loop between feature extraction and TBD of target has been provided in our previous preliminary work [19]. In this paper, a micromotion feature extraction method of space target based on TBD technology is further proposed. In the method, the tracks of target scatterers can be obtained during the TBD of target, and the tracks are the range-slow-time curves of target scatterers indeed. We call the range-slow-time curve as “range trajectory” in this paper. On this basis, we attempt to add the micromotion feature extraction into the process of target detecting and tracking by establishing a feedback loop, namely, that the micromotion feature parameters are extracted via fitting the obtained range trajectory according to the mathematical expression of m-D effect, and in turn the extracted micromotion feature parameters are utilized to update the parameters of TBD adaptively. As a result, the micromotion feature extraction, detecting, and tracking of space target can be implemented simultaneously with the information feedback, which can provide real-time and effective information for target recognition. By changing the mathematical expressions of different kinds of m-D effects for range trajectory fitting, the fixed scatterer model, the slide scatterer model, and any kinds of micromotion forms can be processed with the proposed method.This paper is organized as follows. Taking the conical target contains a fixed scatterer and two slide scatterers as an example, the micromotion feature of the target precession is analyzed in Section2. Combining with the TBD technology, the micromotion feature extraction method is presented in detail in Section 3. Simulations are presented in Section 4 and some conclusions are made in the last section.


## 2. Micromotion Feature of Conical Target with Precession

When the space target is flying outside the atmosphere, it needs to spin around its own axis of symmetry to maintain stability. Meanwhile, it needs to cone around a spatial directional axis due to the lateral force caused by projectile separation and bait release. These two types of rotation formed the precession of target. For a conical target, it usually consists of three dominant scatterers (i.e., one tip-cone scatterer and two cone-base scatterers), where the tip-cone scatterer can be treated as a fixed tip-cone scatterer, and the cone-base scatterers are slide cone-base scatterers. The geometry of conical target with precession is shown in Figure1, where O is the target center, and the target is rotating around z-axis with precession angle ϑ and precession frequency ω rad/s. The distance between the radar and the target center is denoted as Rc. The radar LOS is denoted as n⃑, y-axis is in the plane determined by n⃑ and z-axis, and x-axis is established according to the right-hand rule. The LOS can be represented as n⃑=sin⁡βy→+cos⁡βz→, where β is the angle between LOS and z-axis. The target consists of three scatterers denoted as a, c, and d, respectively. The radius of the base circle is r0, the distance between the target centroid and scatterer a is denoted as oa, and the vertical distance between the target centroid and the bottom circle is represented as ob. According to the rigid body dynamics, the velocity of target should be along z-axis and the velocity is given as v.Figure 1
The geometry of conical target with precession.Assume that the angle betweenx-axis and the projection of vector oa→ on the xoy plane is ϕ0 at time t=0, then the angle between LOS and oa→ is denoted as φ, and it can be calculated by(1)cos⁡φ=sin⁡βsin⁡ϑsin⁡ωt+ϕ0+cos⁡βcos⁡ϑ.The projection ofoa→, oc→, and od→ in the line-of-sight direction at time t can be represented as ra(t), rc(t), and rd(t), respectively; that is,(2)rat=-oa·cos⁡φ-cos⁡β·v·t+Rc=-oasin⁡βsin⁡ϑsin⁡ωt+ϕ0+cos⁡βcos⁡ϑ-cos⁡β·v·t+Rc(3)rct=r0·sin⁡φ+ob·cos⁡φ-cos⁡β·v·t+Rc=r01-cos2⁡βcos2⁡ϑ-sin2⁡βsin2⁡ϑsin2⁡ωt+ϕ0-2cos⁡βcos⁡ϑsin⁡βsin⁡ϑsin⁡ωt+ϕ0+obsin⁡βsin⁡ϑsin⁡ωt+ϕ0+cos⁡βcos⁡ϑ-cos⁡β·v·t+Rc(4)rdt=-r0·sin⁡φ+ob·cos⁡φ-cos⁡β·v·t+Rc=-r01-cos2⁡βcos2⁡ϑ-sin2⁡βsin2⁡ϑsin2⁡ωt+ϕ0-2cos⁡βcos⁡ϑsin⁡βsin⁡ϑsin⁡ωt+ϕ0+obsin⁡βsin⁡ϑsin⁡ωt+ϕ0+cos⁡βcos⁡ϑ-cos⁡β·v·t+Rc.Reference [1] points out that, in the wideband radar system, the micromotion feature of target can be described by range-slow-time image, where the peak of range profile appears to be a range-slow-time curve (range trajectory) which is determined by the scatterer range rt. And the range trajectory reflects the micromotion feature of target scatterer. Therefore, from (2)-(4), the micromotion feature of cone-tip scatterer a is sinusoid, while that of the cone-base scatterers c and d is quasisinusoid, which deviates from the standard sinusoid.During the target observation, the shielded effect should be considered. Assume that the cone half angle of target isΥ. Obviously, φ is changing with time t, and φ∈β-ϑ,β+ϑ. The observable scatterers are different when φ locates in different interval range of [0,π]. Therefore, we divide [0,π] into four observation areas as follows:Area 1: 0≤φ<Υ, the observable scatterers are a, c, and d.Area 2: Υ≤φ<π/2, the observable scatterers are a and d.Area 3: π/2≤φ<π-Υ, the observable scatterers are a, c, and d.Area 4: π-Υ≤φ<π, the observable scatterers are c and d.The impact of the shielded effect on the proposed micromotion feature extraction method will be discussed in Section3.Next, a simple simulation is given to validate the correctness of theoretical analysis. The pulse duration isTp=1 μs, carrier frequency is fc=10GHz, signal bandwidth is B=3GHz, pulse repetition frequencies is PRF=60Hz, and the coherent processing time is Tc=1 s. The radius of the base circle is r0=1m, the distance between the target center and cone-tip scatterer a is oa=3m, and the vertical distance between the target center and the bottom circle is ob=0.3m.To observe the micromotion feature of target more intuitively, we assume thatRc=0m and v=0m/s. The precession angle is ϑ=15°, precession frequency is ω=8π rad/s, the angle between LOS and z-axis is β=145°, and the initial angle is ϕ0=π/100 rad. We can obtain the range-slow-time image as shown in Figure 2.Figure 2
Range-slow-time image of precession target.From Figure2, we can see that the range profile peak of scatterer a changes with slow-time following the sinusoidal form, while those of the sliding scatterers c and d are quasisinusoids (deviating from the standard sinusoid form), which confirms the theoretical analysis.


## 3. Micromotion Feature Extraction Based on TBD

Due to the requirement of additional observation pulses after target detection and tracking, the existing micromotion feature extraction methods usually occupy much radar resources. Also, the real-time performance and radar efficiency are not satisfied. So we intend to combine the micromotion feature extraction with target detecting and tracking in the way of information feedback and then implement the micromotion feature extraction, detecting and tracking of space target simultaneously. It will save the radar resources and improve the radar efficiency and real-time performance of micromotion feature extraction.In this section, we establish a feedback loop between micromotion feature extraction and TBD of target. Firstly, based on TBD technology, the target scatterer range trajectory information is backtracked along with the energy accumulation process, and then the micromotion feature parameters of target can be fitted with these obtained trajectory information. In return, the extracted micromotion feature parameters are used to update the state transition set of TBD adaptively, and the results of micromotion feature extraction are considered into the declaration of the target presence. As a result, micromotion feature extraction, target detecting, and tracking can be completed simultaneously.


### 3.1. Range Trajectory Backtracking and Micromotion Feature Extraction

For simpleness, the range trajectory backtracking for single scatterer based on TBD is described firstly. The monitoring area is divided intoNr×Nθ grids according to the range and azimuthal angle, and each grid is denoted as a state (i,j), i=1,…,Nr, j=1,…,Nθ, which represents the position i-Nr/2·Δr+R0,j-Nθ/2Δθ+θ0, where Δr and Δθ are the stepped increasement of the range and azimuthal angle, and R0,θ0 is the center of the monitoring area. Assume the radar transmits wideband signal pt^. At each scan, the beams towards all the azimuthal angles (Nθ angles) are formed by the radar. The beam width is denoted as ϕB, each state (i,j) will be hit by M+1 successive beams, where(5)M=ϕBΔθ,ϕBΔθis evenϕBΔθ-1,ϕBΔθis odd.Obviously, the sequence number of these M+1 beams should be j-M/2,j-M/2+1,…,j+M/2.At thekth scan, the echo signal of the state xk=ik,jk of the mth beam can be represented as(6)sk,m,xkt^=σkik,jkGmjk·pt^-2ik-Nr/2·Δr+2R0c,t^∈-Tp2,Tp2,m∈1,…,M+1,where Gm· represents the beampattern of the mth beam and Gmjk represents the obtained transmit gain at the jkth azimuthal angle with the mth beam and σk(ik,jk) is the backscattered amplitude of the state (ik,jk) at the kth scan.After performing range compression, the high-resolution range profile (HRRP) can be obtained as(7)Sk,m,xkFf=Tp·Gmjk·σkik,jk·sincFf+2ik-Nr/2·Δr+2R0c.Thus, at thekth scan, the measured value of each state xk=(ik,jk) can be defined as(8)Zkxk=∑m=1M+1Sk,m,xkFfFf=-2ik-Nr/2·Δr+2R0/c2,ik=1,…,Nr,jk=1,…,Nθ.At thekth scan, the cumulative energy of state xk is denoted as Ixk, and it can be calculated as(9)Ixk=Zkxk+maxxk-1∈ΓxkIxk-1,where Γxk is the state transition set. Γxk contains all the possible state xk-1 which can transit to state xk. Set an appropriate detection threshold for the cumulative energy function I(xk) after K scans accumulation, and the state sequence whose cumulative energy is larger than the threshold can be backtracked according to(10)Bkxk=arg⁡maxxk-1∈Γxk⁡Ixk-1,where Bkxk is the backtracking function, which is used to record the state corresponding to the maximum cumulative energy at each scan. Assume that the recorded state of the kth scan is xk=ik,jk, the estimated scatterer range trajectory can be denoted as R(k)=ik-Nr/2·Δr+R0, k=1,…,K, and the angle trajectory is Θ(k)=jk-Nθ/2Δθ+θ0, k=1,…,K.For the cone-tip scatterera, the estimated range trajectory(11)Rk=rak·Δtshould be equal to (2), where Δt is the time interval between the two adjacent scans. Without loss of generality, it holds Δt=1/PRF. Similarly, for the cone-base scatterers c and d, the estimated range trajectories Rk should be equal to (3) and (4), respectively.However, the TBD method described above can obtain only one scatterer range trajectory, which cannot meet the requirement of getting range trajectories of each scatterers. Although some TBD methods for multitarget detecting and tracking have been proposed [20, 21], they require that the state of different targets can not be the same. However, in this paper, the different scatterers are usually at the same azimuthal angle and may be of the same range at some scans (i.e., the range trajectories may be intersected). That means the different scatterers may have the same state, which leads to the existing method which can not be used. To resolve this problem, we improve the TBD method as follows.Assume the target consists ofP observable scatterers. At the kth scan, for each azimuthal angle jk, select any P states to form an expanded state ik,1,jk,ik,2,jk,…,ik,P,jk, where the P states can be the same; that is, ik,p1 can be equal to ik,p2 when p1≠p2. The measured value of the expanded state yk=ik,1,jk,ik,2,jk,…,ik,P,jk is defined as(12)Zkyk=∑p=1PZkik,p,jk.On the basis, the cumulative energy shown as (9) is conducted in terms of the expanded state, where the state correlation is necessary. For example, there are 6 kinds of states correlation approaches when P=3, shown as Figure 3.Figure 3
State correlation approaches.Among these 6 kinds of states correlation approaches, only one is consistent with the actual scatterers trajectories. In [22], we have proposed a state correlation approach which can be used for states correlation.AfterK scans, the cumulative energy Iyk should be compared with a detection threshold Tα:(13)Tα=γKP,where γ is a constant which will affect the target detecting performance, in this paper we call γ as “detection threshold coefficient.” The state sequences whose cumulative energy is larger than the threshold can be backtracked according to (10).It should be pointed out that the states correlation leads to the dependence of the energy accumulation of each scatterer. A specific expression of the detection threshold coefficient is difficult to be derived with a given false alarm probability. Just as [20], the detection threshold coefficient can be chosen from Monte-Carlo experiments.Assume that the number of the cumulative energyIyk which is larger than the threshold is Q. Thus, Q state sequences will be obtained, and each state sequence contains P range trajectories Rp(k)=ik-Nr/2·Δr+R0   p=1,…,P,k=1,…,K and P angle trajectories Θp(k)=jk-Nθ/2Δθ+θ0   p=1,…,P,k=1,…,K. Due to the fact that intersections number of any two range trajectories will be small, a condition for selecting the reasonable state sequence from the Q state sequences is defined as(14)CRp1k=Rp2k<ςK,∀p1≠p2,where CRp1k=Rp2k represents the intersections number of Rp1(k) and Rp2(k) and ς is a constant used to control the ratio of intersections number to scans number K. The state sequences which satisfy the condition shown as (14) are selected, and the corresponding range trajectories and angle trajectories are obtained from the target scatterers.For each possible value ofP, the TBD procedure proposed above can be conducted to obtain the range trajectories of P scatterers. The shielded effect has been discussed in Section 2: when φ falls with inArea 1 orArea 3, three scatterers are observable (corresponding to P=3); when φ falls withinArea 2 orArea 4, two scatterers are observable (corresponding to P=2). Therefore, the possible values of P are P=3 and P=2. Thus, the TBD procedure of P=3 and P=2 should be conducted.In the case ofP=3, if there is state sequence which satisfies (14) which can be selected out, three range trajectories (i.e., R1k, R2k, and R3k) can be obtained. Firstly, we assume that R1k is the range trajectory of the cone-tip scatterer a and R2k and R3k are the range trajectories of the cone-base scatterers c and d, respectively. R1(k), R2(k), and R3(k) can be represented as(15)R1k=rak·Δt+εk,R2k=rck·Δt+εk,R3k=rdk·Δt+εk,where ε(k) represents the error induced from the noise and TBD procedure. We can estimate the micromotion feature parameter vector PA=oa,ϕ0,β,ϑ,ω,v,ob,r0,Rc by fitting R1k, R2k, and R3k according to the curve form of standard sinusoid shown as (2) and quasisinusoid shown as (3)-(4) with the least squares method:(16)minPAR1k+oasin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0+cos⁡βcos⁡ϑ+cos⁡β·vk·Δt-Rc+R2k-r01-cos2⁡βcos2⁡ϑ-sin2⁡βsin2⁡ϑsin2⁡ωk·Δt+ϕ0-2cos⁡βcos⁡ϑsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0-obsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0+cos⁡βcos⁡ϑ+cos⁡β·vk·Δt-Rc+R3k+r01-cos2⁡βcos2⁡ϑ-sin2⁡βsin2⁡ϑsin2⁡ωk·Δt+ϕ0-2cos⁡βcos⁡ϑsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0-obsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0+cos⁡βcos⁡ϑ+cos⁡β·vk·Δt-Rc.Equation (16) can be solved with the Levenberg–Marquardt method [23], which is sensitive to the initial values. Therefore, how to set the appropriate initial values is proposed as follows. The EMD method [24] can be used to separate R1(k) into a set of intrinsic mode functions (IMF) which is descended by frequency. Thus, we can obtain(17)I1,R1k=Ba+Cak,I2,R1k=Aasin⁡ωak+ϕ0aAa=-oasin⁡βsin⁡ϑ+ξ,Ba=-oacos⁡βcos⁡ϑ+Rc+ξ,Ca=-cos⁡β·v+ξ,ωa=ω+ξ,ϕ0a=ϕ0+ξ,where ξ is the error from ε(k) and EMD method. According to (3) and (4), we can get(18)R+k=R2k+R3k=2obsin⁡βsin⁡ϑsin⁡ωk+ϕ0+cos⁡βcos⁡ϑ-2cos⁡β·vk+2Rc+εkR-k=R2k-R3k=2r01-cos2⁡βcos2⁡ϑ-sin2⁡βsin2⁡ϑsin2⁡ωk·Δt+ϕ0-2cos⁡βcos⁡ϑsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0+εk.Similarly, separatingR+(k) with the EMD method, we can obtain(19)I1,R+k=Bcd+Ccdk,I2,R+k=Acdsin⁡ωcdk+ϕ0cdAcd=2obsin⁡βsin⁡ϑ+ξ,Bcd=2obcos⁡βcos⁡ϑ+2Rc+ξ,Ccd=-2cos⁡β·v+ξ,ωcd=ω+ξ,ϕ0cd=ϕ0+ξ.The initial values|oa|ini, ϕ0ini, βini, ϑini, ωini, vini, |ob|ini, r0ini, and Rcini for (19) can be obtained according the equation set:(20)-oainicos⁡βinicos⁡ϑini+Rcini=I1,R102obinicos⁡βinicos⁡ϑini+2Rcini=I1,R+0cos⁡βini·vini=-I1,R1k2-I1,R1k1/k2-k1+I1,R+k2-I1,R+k1/k2-k12oainisin⁡βinisin⁡ϑini=max⁡I2,R1-min⁡I2,R122obinisin⁡βinisin⁡ϑini=max⁡I1,R+-min⁡I1,R+2ϕ0ini=arsin⁡I2,R10/-oainisin⁡βinisin⁡ϑini+arsin⁡I2,R+0/2obinisin⁡βinisin⁡ϑini2,sin⁡βinisin⁡ϑini>0arsin⁡I2,R10/oainisin⁡βinisin⁡ϑini+arsin⁡I2,R+0/-2obinisin⁡βinisin⁡ϑini2,sin⁡βinisin⁡ϑini<0ωini=arg⁡maxω⁡FFTI2,R1+arg⁡maxω⁡FFTI2,R+k22r0ini1-cos2⁡βinicos2⁡ϑini-sin2⁡βinisin2⁡ϑinisin2⁡ϕ0ini-2cos⁡βinicos⁡ϑinisin⁡βinisin⁡ϑinisin⁡ϕ0ini=R-0,where arsin⁡· is the inverse function of sin⁡·. In (20), the equations’ number is one less than unknown parameters number. Thus, we defined the search interval and the search stepped increasement of r0 as rmin,rmax and Δr0, respectively. For each initial values r0inisn=rmin+sn-1·Δr0, sn=1,2,…,SN, and SN=rmax-rmin/Δr0, a group of initial values of |oa|ini(sn), ϕ0ini(sn), βini(sn), ϑini(sn), ωini(sn), vini(sn), obinisn, and Rcinisn can be obtained according to (20). Based on the initial values, the corresponding micromotion feature parameter vector PAsn=[oasn,ϕ0sn,βsn,ϑsn, ωsn,vsn,obsn, r0sn,Rcsn] can be obtained by solving (16) with the Levenberg–Marquardt method. The fitting error with the snth group initial values is calculated as(21)Esn=13KR1k+oasnF1sn+R2k-r0snF2sn-obsnF1sn+R3k+r0snF2sn-obsnF1sn,F1sn=sin⁡βsnsin⁡ϑsnsin⁡ωsnk+ϕ0sn+cos⁡βsncos⁡ϑsn+cos⁡βsn·vsnk-Rcsn,F2sn=1-cos2⁡βsncos2⁡ϑsn-sin2⁡βsnsin2⁡ϑsnsin2⁡ωsnk+ϕ0sn-2cos⁡βsncos⁡ϑsnsin⁡βsnsin⁡ϑsnsin⁡ωsnk+ϕ0sn.AllE(sn) is compared with each other, and the estimated micromotion feature parameter vector is defined as(22)PA^≜PAarg⁡minsn⁡Esn.The fitting error is denoted as(23)E^≜minsn⁡Esn.It should be pointed thatPA^ and E^ are obtained under the assumption that R1(k) is the range trajectory of the cone-tip scatterer a. Therefore, we rewrite them as PA^(R1) and E^R1.Next, we assumeR2(k) and R3k are the range trajectory of the cone-tip scatterer, respectively. The micromotion feature parameter vectors PA^R2 and PA^R3 can be obtained, and the corresponding fitting errors are denoted as E^R2 and E^R3. The micromotion feature parameter vector corresponding to the minimum fitting error is selected out as the final micromotion feature parameter vector(24)PA^f≜PA^arg⁡minR1,R2,R3⁡E^.The fitting error is denoted as(25)E^f≜minR1,R2,R3⁡E^.Similarly, in the case ofP=2, the micromotion feature parameter vector can be obtained. What should be pointed out is that if the two range trajectories are corresponding to the cone-base scatterers c and d, respectively, the parameter |oa| can not be obtained.


### 3.2. Adaptive Update of State Transition Set

Just as mentioned in (18), the state transition set contains all the possible state yk-1 which can transit to state yk, and the construction of state transition set will significantly affect the efficiency and performance of TBD and micromotion feature extraction. In this paper, the state transition set is updated adaptively, by forecasting the state of scatterers in next scan according to the extracted micromotion feature parameters.Assume that, in thekth scan, the cumulative energy Iyk is larger than the threshold; the extracted micromotion feature parameter vector can be obtained with the method proposed in Section 3.1, which is represented as(26)PA^f,yk=oayk,ϕ0yk,βyk,ϑyk,ωyk,vyk,obyk,r0yk,Rcyk.The corresponding fitting error is denoted asE^f,yk.The state transition set can be determined in light of the extracted micromotion feature parameters. The stateyk will belong to the state transition set Γyk+1 (i.e., yk∈Γ(yk+1)) when it satisfies(27)ik+1,1-Nr2·Δr+R0-Δra<fRCs,yk·Δr,jk+1-jk<2,for scattereraik+1,2-Nr2·Δr+R0-Δrc<fRCs,yk·Δr,jk+1-jk<2,for scatterercik+1,3-Nr2·Δr+R0-Δrd<fRCs,yk·Δr,jk+1-jk<2,for scattererd,where(28)Δra=-oayksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+cos⁡βykcos⁡ϑyk-cos⁡βykvykk+1+RcykΔrc=r0yk1-cos2⁡βykcos2⁡ϑyk-sin2⁡βyksin2⁡ϑyksin2⁡ωykk+1+ϕ0yk-2cos⁡βykcos⁡ϑyksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+obyksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+cos⁡βykcos⁡ϑyk-cos⁡βyk·vykk+1+Rcyk.Δrd=-r0yk1-cos2⁡βykcos2⁡ϑyk-sin2⁡βyksin2⁡ϑyksin2⁡ωykk+1+ϕ0yk-2cos⁡βykcos⁡ϑyksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+obyksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+cos⁡βykcos⁡ϑyk-cos⁡βyk·vykk+1+Rc,ykCs,yk represents the consistency of the micromotion feature parameter vectors between k-1th scan and kth scan, and fR· is an adaptive adjustment function which is used to control the size of searching range gate during TBD according to Cs,yk. In this paper, Cs,yk is defined as(29)Cs,yk=meanPA^f,yk-PA^f,yk-1PA^f,yk-1,where mean· represents the mean value of a vector. Obviously, the smaller Cs,yk is, the higher consistency is.Obviously, ifCs,yk is relatively small, it shows that the consistency of the extracted micromotion feature parameter vectors is high, we can consider that the precision of the micromotion feature extraction is well, and the forecast of state information in the next scan is accurate. In this case, the searching range gate can be decreased to reduce the computation load. On the contrary, if Cs,yk is relatively large, the forecasted precision is low and the searching range gate should be increased appropriately. Therefore, fR· should be an increasing function, which is defined as(30)fRCs,yk=rm+λ·Cs,yk,where λ is a constant coefficient, rm is the minimum of searching range gate, which is set to be 5 in this paper, and the maximum value of fR(·) is limited to 20.


### 3.3. Adaptive Starting and Ending of Precession Feature Extraction

In TBD algorithms, the detection threshold affects the target detecting performance directly. However, how to set an appropriate detection threshold is difficult, especially in the multiscatterers environment.As the micromotion feature extraction has been included into the process of target detecting and tracking by establishing a feedback loop, the pertinent issue now is how to declare the presence of a target and when to start and end the micromotion feature extraction algorithm.In the proposed method, the energy accumulation value and the result of micromotion feature extraction are both taken into full consideration to address the problem above, which is different from the traditional signal processing method. Just as we know, in traditional method, only when the energy accumulation value is larger than the threshold, it declares the presence of a target, and then the extraction of its micromotion feature will be implemented by transmitting additional observation pulses. In comparison, in the proposed method, by adding the micromotion feature extraction into the process of target detecting and tracking, on condition that the precision of micromotion feature extraction is high enough, it can declare the presence of a target and the micromotion feature is extracted successfully at the same time, although the energy accumulation value dose not achieve the required level.Assume the minimum and maximum total number of scans that are jointly processed in TBD areKN and KM, respectively. Firstly, set two detection threshold coefficients: the lower detection threshold coefficient γ1 and the higher detection threshold coefficient γ2. In the kth (k≥KN) scan if the cumulative energy Iyk is larger than Tα2=γ2·k·P, we declare the presence of a target and the micromotion feature parameters can be obtained with the proposed method described in Section 3.1. On the other hand, if the cumulative energy I(yk) is larger than Tα1=γ1·k·P and smaller than Tα2, backtracking the target scatterers trajectories according to (10) and (14). On this basis, the extracted micromotion feature parameter vector at kth scan can be obtained. Go on to cumulate energy for the data of k+1th scan, and extract the micromotion feature parameters at k+1th scan. The consistency of the extracted micromotion feature parameter vectors Cs,yk can be calculated according to (29). If Cs,yk and fitting error E^f,yk are both relatively small (satisfies Cs,yk<TC and E^f,yk<TE), we can declare the presence of a target and get the micromotion feature parameters. The energy accumulation is no longer needed. Otherwise, update the state transition set according to (27), and the energy accumulation of the data of k+2th scan is needed. Repeat the steps described above until it satisfies Cs,yk<TC and E^f,yk<TE or Iyk>Tα2, or it reached the KMth scan. In conclusion, the flowchart of micromotion feature extraction of space target based on TBD is shown in Figure 4.Figure 4
Micromotion feature extraction of space target based on TBD.There are three points should be noticed for the proposed method.(1) The values ofγ1 and γ2 are chosen from Monte-Carlo experiments with the given false alarm probability PFA,γ1 and PFA,γ2 in the case the presence of a target is declared when the energy accumulation value is larger than the threshold. The corresponding detection probabilities are denoted as PD,γ1 and PD,γ2. However, the result of micromotion feature extraction is taken into consideration to declare the presence of a target in this paper; thus the final false alarm probability PFA and detection probability PD obtained from the proposed method are hard to be calculated from a specific expression with the variables of PFA,γ1, PD,γ1, PFA,γ2, and PD,γ2. Numerous experiments have shown that PFA will be a little higher than PFA,γ2. Therefore the value of PFA,γ2 can be chosen according to the desired false alarm probability, which is set as PFA,γ2=0.005 in this paper. The value of PFA,γ1 will affect the detection probability and computation load, and it is set as PFA,γ1=0.5 in this paper from numerous experiments with the consideration of maximizing the detection probability and minimizing the computation load.(2) When the value ofP is larger than the observable scatterers number, the cumulative energy and the trajectories which are corresponding to the target will not be different from that corresponding to the noise. Thus, the target will not be detected. On the contrary, if the value of P is smaller than the observable scatterers number, only P trajectories can be obtained. Thus, some scatterers trajectories will be lost and the corresponding micromotion feature parameters can not be extracted. Therefore, the value of P should be set as the maximum possible value firstly, which can be preset by experience. Then the value of P decreases gradually until a target is detected and the micromotion feature is extracted or P=1.(3) The computational load of the proposed method is proportional to the search range. Therefore, to reduce the computational load, we assume that target detecting and tracking with narrow-band radar is conducted firstly, and a lower threshold is used to detect the possible targets; then the obtained coarse position and velocity information of targets are utilized to determine the search range of the proposed method.


## 3.1. Range Trajectory Backtracking and Micromotion Feature Extraction

For simpleness, the range trajectory backtracking for single scatterer based on TBD is described firstly. The monitoring area is divided intoNr×Nθ grids according to the range and azimuthal angle, and each grid is denoted as a state (i,j), i=1,…,Nr, j=1,…,Nθ, which represents the position i-Nr/2·Δr+R0,j-Nθ/2Δθ+θ0, where Δr and Δθ are the stepped increasement of the range and azimuthal angle, and R0,θ0 is the center of the monitoring area. Assume the radar transmits wideband signal pt^. At each scan, the beams towards all the azimuthal angles (Nθ angles) are formed by the radar. The beam width is denoted as ϕB, each state (i,j) will be hit by M+1 successive beams, where(5)M=ϕBΔθ,ϕBΔθis evenϕBΔθ-1,ϕBΔθis odd.Obviously, the sequence number of these M+1 beams should be j-M/2,j-M/2+1,…,j+M/2.At thekth scan, the echo signal of the state xk=ik,jk of the mth beam can be represented as(6)sk,m,xkt^=σkik,jkGmjk·pt^-2ik-Nr/2·Δr+2R0c,t^∈-Tp2,Tp2,m∈1,…,M+1,where Gm· represents the beampattern of the mth beam and Gmjk represents the obtained transmit gain at the jkth azimuthal angle with the mth beam and σk(ik,jk) is the backscattered amplitude of the state (ik,jk) at the kth scan.After performing range compression, the high-resolution range profile (HRRP) can be obtained as(7)Sk,m,xkFf=Tp·Gmjk·σkik,jk·sincFf+2ik-Nr/2·Δr+2R0c.Thus, at thekth scan, the measured value of each state xk=(ik,jk) can be defined as(8)Zkxk=∑m=1M+1Sk,m,xkFfFf=-2ik-Nr/2·Δr+2R0/c2,ik=1,…,Nr,jk=1,…,Nθ.At thekth scan, the cumulative energy of state xk is denoted as Ixk, and it can be calculated as(9)Ixk=Zkxk+maxxk-1∈ΓxkIxk-1,where Γxk is the state transition set. Γxk contains all the possible state xk-1 which can transit to state xk. Set an appropriate detection threshold for the cumulative energy function I(xk) after K scans accumulation, and the state sequence whose cumulative energy is larger than the threshold can be backtracked according to(10)Bkxk=arg⁡maxxk-1∈Γxk⁡Ixk-1,where Bkxk is the backtracking function, which is used to record the state corresponding to the maximum cumulative energy at each scan. Assume that the recorded state of the kth scan is xk=ik,jk, the estimated scatterer range trajectory can be denoted as R(k)=ik-Nr/2·Δr+R0, k=1,…,K, and the angle trajectory is Θ(k)=jk-Nθ/2Δθ+θ0, k=1,…,K.For the cone-tip scatterera, the estimated range trajectory(11)Rk=rak·Δtshould be equal to (2), where Δt is the time interval between the two adjacent scans. Without loss of generality, it holds Δt=1/PRF. Similarly, for the cone-base scatterers c and d, the estimated range trajectories Rk should be equal to (3) and (4), respectively.However, the TBD method described above can obtain only one scatterer range trajectory, which cannot meet the requirement of getting range trajectories of each scatterers. Although some TBD methods for multitarget detecting and tracking have been proposed [20, 21], they require that the state of different targets can not be the same. However, in this paper, the different scatterers are usually at the same azimuthal angle and may be of the same range at some scans (i.e., the range trajectories may be intersected). That means the different scatterers may have the same state, which leads to the existing method which can not be used. To resolve this problem, we improve the TBD method as follows.Assume the target consists ofP observable scatterers. At the kth scan, for each azimuthal angle jk, select any P states to form an expanded state ik,1,jk,ik,2,jk,…,ik,P,jk, where the P states can be the same; that is, ik,p1 can be equal to ik,p2 when p1≠p2. The measured value of the expanded state yk=ik,1,jk,ik,2,jk,…,ik,P,jk is defined as(12)Zkyk=∑p=1PZkik,p,jk.On the basis, the cumulative energy shown as (9) is conducted in terms of the expanded state, where the state correlation is necessary. For example, there are 6 kinds of states correlation approaches when P=3, shown as Figure 3.Figure 3
State correlation approaches.Among these 6 kinds of states correlation approaches, only one is consistent with the actual scatterers trajectories. In [22], we have proposed a state correlation approach which can be used for states correlation.AfterK scans, the cumulative energy Iyk should be compared with a detection threshold Tα:(13)Tα=γKP,where γ is a constant which will affect the target detecting performance, in this paper we call γ as “detection threshold coefficient.” The state sequences whose cumulative energy is larger than the threshold can be backtracked according to (10).It should be pointed out that the states correlation leads to the dependence of the energy accumulation of each scatterer. A specific expression of the detection threshold coefficient is difficult to be derived with a given false alarm probability. Just as [20], the detection threshold coefficient can be chosen from Monte-Carlo experiments.Assume that the number of the cumulative energyIyk which is larger than the threshold is Q. Thus, Q state sequences will be obtained, and each state sequence contains P range trajectories Rp(k)=ik-Nr/2·Δr+R0   p=1,…,P,k=1,…,K and P angle trajectories Θp(k)=jk-Nθ/2Δθ+θ0   p=1,…,P,k=1,…,K. Due to the fact that intersections number of any two range trajectories will be small, a condition for selecting the reasonable state sequence from the Q state sequences is defined as(14)CRp1k=Rp2k<ςK,∀p1≠p2,where CRp1k=Rp2k represents the intersections number of Rp1(k) and Rp2(k) and ς is a constant used to control the ratio of intersections number to scans number K. The state sequences which satisfy the condition shown as (14) are selected, and the corresponding range trajectories and angle trajectories are obtained from the target scatterers.For each possible value ofP, the TBD procedure proposed above can be conducted to obtain the range trajectories of P scatterers. The shielded effect has been discussed in Section 2: when φ falls with inArea 1 orArea 3, three scatterers are observable (corresponding to P=3); when φ falls withinArea 2 orArea 4, two scatterers are observable (corresponding to P=2). Therefore, the possible values of P are P=3 and P=2. Thus, the TBD procedure of P=3 and P=2 should be conducted.In the case ofP=3, if there is state sequence which satisfies (14) which can be selected out, three range trajectories (i.e., R1k, R2k, and R3k) can be obtained. Firstly, we assume that R1k is the range trajectory of the cone-tip scatterer a and R2k and R3k are the range trajectories of the cone-base scatterers c and d, respectively. R1(k), R2(k), and R3(k) can be represented as(15)R1k=rak·Δt+εk,R2k=rck·Δt+εk,R3k=rdk·Δt+εk,where ε(k) represents the error induced from the noise and TBD procedure. We can estimate the micromotion feature parameter vector PA=oa,ϕ0,β,ϑ,ω,v,ob,r0,Rc by fitting R1k, R2k, and R3k according to the curve form of standard sinusoid shown as (2) and quasisinusoid shown as (3)-(4) with the least squares method:(16)minPAR1k+oasin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0+cos⁡βcos⁡ϑ+cos⁡β·vk·Δt-Rc+R2k-r01-cos2⁡βcos2⁡ϑ-sin2⁡βsin2⁡ϑsin2⁡ωk·Δt+ϕ0-2cos⁡βcos⁡ϑsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0-obsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0+cos⁡βcos⁡ϑ+cos⁡β·vk·Δt-Rc+R3k+r01-cos2⁡βcos2⁡ϑ-sin2⁡βsin2⁡ϑsin2⁡ωk·Δt+ϕ0-2cos⁡βcos⁡ϑsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0-obsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0+cos⁡βcos⁡ϑ+cos⁡β·vk·Δt-Rc.Equation (16) can be solved with the Levenberg–Marquardt method [23], which is sensitive to the initial values. Therefore, how to set the appropriate initial values is proposed as follows. The EMD method [24] can be used to separate R1(k) into a set of intrinsic mode functions (IMF) which is descended by frequency. Thus, we can obtain(17)I1,R1k=Ba+Cak,I2,R1k=Aasin⁡ωak+ϕ0aAa=-oasin⁡βsin⁡ϑ+ξ,Ba=-oacos⁡βcos⁡ϑ+Rc+ξ,Ca=-cos⁡β·v+ξ,ωa=ω+ξ,ϕ0a=ϕ0+ξ,where ξ is the error from ε(k) and EMD method. According to (3) and (4), we can get(18)R+k=R2k+R3k=2obsin⁡βsin⁡ϑsin⁡ωk+ϕ0+cos⁡βcos⁡ϑ-2cos⁡β·vk+2Rc+εkR-k=R2k-R3k=2r01-cos2⁡βcos2⁡ϑ-sin2⁡βsin2⁡ϑsin2⁡ωk·Δt+ϕ0-2cos⁡βcos⁡ϑsin⁡βsin⁡ϑsin⁡ωk·Δt+ϕ0+εk.Similarly, separatingR+(k) with the EMD method, we can obtain(19)I1,R+k=Bcd+Ccdk,I2,R+k=Acdsin⁡ωcdk+ϕ0cdAcd=2obsin⁡βsin⁡ϑ+ξ,Bcd=2obcos⁡βcos⁡ϑ+2Rc+ξ,Ccd=-2cos⁡β·v+ξ,ωcd=ω+ξ,ϕ0cd=ϕ0+ξ.The initial values|oa|ini, ϕ0ini, βini, ϑini, ωini, vini, |ob|ini, r0ini, and Rcini for (19) can be obtained according the equation set:(20)-oainicos⁡βinicos⁡ϑini+Rcini=I1,R102obinicos⁡βinicos⁡ϑini+2Rcini=I1,R+0cos⁡βini·vini=-I1,R1k2-I1,R1k1/k2-k1+I1,R+k2-I1,R+k1/k2-k12oainisin⁡βinisin⁡ϑini=max⁡I2,R1-min⁡I2,R122obinisin⁡βinisin⁡ϑini=max⁡I1,R+-min⁡I1,R+2ϕ0ini=arsin⁡I2,R10/-oainisin⁡βinisin⁡ϑini+arsin⁡I2,R+0/2obinisin⁡βinisin⁡ϑini2,sin⁡βinisin⁡ϑini>0arsin⁡I2,R10/oainisin⁡βinisin⁡ϑini+arsin⁡I2,R+0/-2obinisin⁡βinisin⁡ϑini2,sin⁡βinisin⁡ϑini<0ωini=arg⁡maxω⁡FFTI2,R1+arg⁡maxω⁡FFTI2,R+k22r0ini1-cos2⁡βinicos2⁡ϑini-sin2⁡βinisin2⁡ϑinisin2⁡ϕ0ini-2cos⁡βinicos⁡ϑinisin⁡βinisin⁡ϑinisin⁡ϕ0ini=R-0,where arsin⁡· is the inverse function of sin⁡·. In (20), the equations’ number is one less than unknown parameters number. Thus, we defined the search interval and the search stepped increasement of r0 as rmin,rmax and Δr0, respectively. For each initial values r0inisn=rmin+sn-1·Δr0, sn=1,2,…,SN, and SN=rmax-rmin/Δr0, a group of initial values of |oa|ini(sn), ϕ0ini(sn), βini(sn), ϑini(sn), ωini(sn), vini(sn), obinisn, and Rcinisn can be obtained according to (20). Based on the initial values, the corresponding micromotion feature parameter vector PAsn=[oasn,ϕ0sn,βsn,ϑsn, ωsn,vsn,obsn, r0sn,Rcsn] can be obtained by solving (16) with the Levenberg–Marquardt method. The fitting error with the snth group initial values is calculated as(21)Esn=13KR1k+oasnF1sn+R2k-r0snF2sn-obsnF1sn+R3k+r0snF2sn-obsnF1sn,F1sn=sin⁡βsnsin⁡ϑsnsin⁡ωsnk+ϕ0sn+cos⁡βsncos⁡ϑsn+cos⁡βsn·vsnk-Rcsn,F2sn=1-cos2⁡βsncos2⁡ϑsn-sin2⁡βsnsin2⁡ϑsnsin2⁡ωsnk+ϕ0sn-2cos⁡βsncos⁡ϑsnsin⁡βsnsin⁡ϑsnsin⁡ωsnk+ϕ0sn.AllE(sn) is compared with each other, and the estimated micromotion feature parameter vector is defined as(22)PA^≜PAarg⁡minsn⁡Esn.The fitting error is denoted as(23)E^≜minsn⁡Esn.It should be pointed thatPA^ and E^ are obtained under the assumption that R1(k) is the range trajectory of the cone-tip scatterer a. Therefore, we rewrite them as PA^(R1) and E^R1.Next, we assumeR2(k) and R3k are the range trajectory of the cone-tip scatterer, respectively. The micromotion feature parameter vectors PA^R2 and PA^R3 can be obtained, and the corresponding fitting errors are denoted as E^R2 and E^R3. The micromotion feature parameter vector corresponding to the minimum fitting error is selected out as the final micromotion feature parameter vector(24)PA^f≜PA^arg⁡minR1,R2,R3⁡E^.The fitting error is denoted as(25)E^f≜minR1,R2,R3⁡E^.Similarly, in the case ofP=2, the micromotion feature parameter vector can be obtained. What should be pointed out is that if the two range trajectories are corresponding to the cone-base scatterers c and d, respectively, the parameter |oa| can not be obtained.


## 3.2. Adaptive Update of State Transition Set

Just as mentioned in (18), the state transition set contains all the possible state yk-1 which can transit to state yk, and the construction of state transition set will significantly affect the efficiency and performance of TBD and micromotion feature extraction. In this paper, the state transition set is updated adaptively, by forecasting the state of scatterers in next scan according to the extracted micromotion feature parameters.Assume that, in thekth scan, the cumulative energy Iyk is larger than the threshold; the extracted micromotion feature parameter vector can be obtained with the method proposed in Section 3.1, which is represented as(26)PA^f,yk=oayk,ϕ0yk,βyk,ϑyk,ωyk,vyk,obyk,r0yk,Rcyk.The corresponding fitting error is denoted asE^f,yk.The state transition set can be determined in light of the extracted micromotion feature parameters. The stateyk will belong to the state transition set Γyk+1 (i.e., yk∈Γ(yk+1)) when it satisfies(27)ik+1,1-Nr2·Δr+R0-Δra<fRCs,yk·Δr,jk+1-jk<2,for scattereraik+1,2-Nr2·Δr+R0-Δrc<fRCs,yk·Δr,jk+1-jk<2,for scatterercik+1,3-Nr2·Δr+R0-Δrd<fRCs,yk·Δr,jk+1-jk<2,for scattererd,where(28)Δra=-oayksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+cos⁡βykcos⁡ϑyk-cos⁡βykvykk+1+RcykΔrc=r0yk1-cos2⁡βykcos2⁡ϑyk-sin2⁡βyksin2⁡ϑyksin2⁡ωykk+1+ϕ0yk-2cos⁡βykcos⁡ϑyksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+obyksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+cos⁡βykcos⁡ϑyk-cos⁡βyk·vykk+1+Rcyk.Δrd=-r0yk1-cos2⁡βykcos2⁡ϑyk-sin2⁡βyksin2⁡ϑyksin2⁡ωykk+1+ϕ0yk-2cos⁡βykcos⁡ϑyksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+obyksin⁡βyksin⁡ϑyksin⁡ωykk+1+ϕ0yk+cos⁡βykcos⁡ϑyk-cos⁡βyk·vykk+1+Rc,ykCs,yk represents the consistency of the micromotion feature parameter vectors between k-1th scan and kth scan, and fR· is an adaptive adjustment function which is used to control the size of searching range gate during TBD according to Cs,yk. In this paper, Cs,yk is defined as(29)Cs,yk=meanPA^f,yk-PA^f,yk-1PA^f,yk-1,where mean· represents the mean value of a vector. Obviously, the smaller Cs,yk is, the higher consistency is.Obviously, ifCs,yk is relatively small, it shows that the consistency of the extracted micromotion feature parameter vectors is high, we can consider that the precision of the micromotion feature extraction is well, and the forecast of state information in the next scan is accurate. In this case, the searching range gate can be decreased to reduce the computation load. On the contrary, if Cs,yk is relatively large, the forecasted precision is low and the searching range gate should be increased appropriately. Therefore, fR· should be an increasing function, which is defined as(30)fRCs,yk=rm+λ·Cs,yk,where λ is a constant coefficient, rm is the minimum of searching range gate, which is set to be 5 in this paper, and the maximum value of fR(·) is limited to 20.


## 3.3. Adaptive Starting and Ending of Precession Feature Extraction

In TBD algorithms, the detection threshold affects the target detecting performance directly. However, how to set an appropriate detection threshold is difficult, especially in the multiscatterers environment.As the micromotion feature extraction has been included into the process of target detecting and tracking by establishing a feedback loop, the pertinent issue now is how to declare the presence of a target and when to start and end the micromotion feature extraction algorithm.In the proposed method, the energy accumulation value and the result of micromotion feature extraction are both taken into full consideration to address the problem above, which is different from the traditional signal processing method. Just as we know, in traditional method, only when the energy accumulation value is larger than the threshold, it declares the presence of a target, and then the extraction of its micromotion feature will be implemented by transmitting additional observation pulses. In comparison, in the proposed method, by adding the micromotion feature extraction into the process of target detecting and tracking, on condition that the precision of micromotion feature extraction is high enough, it can declare the presence of a target and the micromotion feature is extracted successfully at the same time, although the energy accumulation value dose not achieve the required level.Assume the minimum and maximum total number of scans that are jointly processed in TBD areKN and KM, respectively. Firstly, set two detection threshold coefficients: the lower detection threshold coefficient γ1 and the higher detection threshold coefficient γ2. In the kth (k≥KN) scan if the cumulative energy Iyk is larger than Tα2=γ2·k·P, we declare the presence of a target and the micromotion feature parameters can be obtained with the proposed method described in Section 3.1. On the other hand, if the cumulative energy I(yk) is larger than Tα1=γ1·k·P and smaller than Tα2, backtracking the target scatterers trajectories according to (10) and (14). On this basis, the extracted micromotion feature parameter vector at kth scan can be obtained. Go on to cumulate energy for the data of k+1th scan, and extract the micromotion feature parameters at k+1th scan. The consistency of the extracted micromotion feature parameter vectors Cs,yk can be calculated according to (29). If Cs,yk and fitting error E^f,yk are both relatively small (satisfies Cs,yk<TC and E^f,yk<TE), we can declare the presence of a target and get the micromotion feature parameters. The energy accumulation is no longer needed. Otherwise, update the state transition set according to (27), and the energy accumulation of the data of k+2th scan is needed. Repeat the steps described above until it satisfies Cs,yk<TC and E^f,yk<TE or Iyk>Tα2, or it reached the KMth scan. In conclusion, the flowchart of micromotion feature extraction of space target based on TBD is shown in Figure 4.Figure 4
Micromotion feature extraction of space target based on TBD.There are three points should be noticed for the proposed method.(1) The values ofγ1 and γ2 are chosen from Monte-Carlo experiments with the given false alarm probability PFA,γ1 and PFA,γ2 in the case the presence of a target is declared when the energy accumulation value is larger than the threshold. The corresponding detection probabilities are denoted as PD,γ1 and PD,γ2. However, the result of micromotion feature extraction is taken into consideration to declare the presence of a target in this paper; thus the final false alarm probability PFA and detection probability PD obtained from the proposed method are hard to be calculated from a specific expression with the variables of PFA,γ1, PD,γ1, PFA,γ2, and PD,γ2. Numerous experiments have shown that PFA will be a little higher than PFA,γ2. Therefore the value of PFA,γ2 can be chosen according to the desired false alarm probability, which is set as PFA,γ2=0.005 in this paper. The value of PFA,γ1 will affect the detection probability and computation load, and it is set as PFA,γ1=0.5 in this paper from numerous experiments with the consideration of maximizing the detection probability and minimizing the computation load.(2) When the value ofP is larger than the observable scatterers number, the cumulative energy and the trajectories which are corresponding to the target will not be different from that corresponding to the noise. Thus, the target will not be detected. On the contrary, if the value of P is smaller than the observable scatterers number, only P trajectories can be obtained. Thus, some scatterers trajectories will be lost and the corresponding micromotion feature parameters can not be extracted. Therefore, the value of P should be set as the maximum possible value firstly, which can be preset by experience. Then the value of P decreases gradually until a target is detected and the micromotion feature is extracted or P=1.(3) The computational load of the proposed method is proportional to the search range. Therefore, to reduce the computational load, we assume that target detecting and tracking with narrow-band radar is conducted firstly, and a lower threshold is used to detect the possible targets; then the obtained coarse position and velocity information of targets are utilized to determine the search range of the proposed method.


## 4. Simulations

In this section, some simulations are carried out to verify the efficiency of the proposed algorithm.Assume that the target center is located at0,0,1000 km at the beginning of observation (i.e., Rc=1000km) and the velocity of target is 500 m/s. The other simulation parameters are the same with those described in Section 2. The geometry of radar and target is illustrated as shown in Figure 5.Figure 5
Geometry of radar and target.The parameters of the TBD procedure are set as follows:KN=15, KM=30, Nr=Nθ=1000, Δr=0.05m, Δθ=0.005°, R0=1000 m, θ0=0°, ϕB=0.15°, ς=0.17,  TC=0.1, TE=2·Δr=0.1 m, PFA,γ1=0.5, and PFA,γ2=0.005. First, we consider the problem of the detection threshold coefficient setting. When the energy accumulation value is larger than the threshold, it declares the presence of a target, and then the false alarm probability versus detection threshold coefficient is reported in Figure 6. Therefore, for the given PFA,γ1=0.5 and PFA,γ2=0.005, the detection threshold coefficients are set as γ1=1 and γ2=1.3.Figure 6
False alarm probability versus detection threshold coefficient.The micromotion feature extraction based on TBD ofP=3 is conducted with SNR=8 dB which is added to the HRRP of the target. Due to the fact that Fourier transform with Nr (Nr=1000) points is taken to obtain the HRRP, the average received SNR of the original echo is obtained by subtracting 10log⁡1000=30 dB from the reported value. After KN=15 scans, among all the state sequences whose cumulative energy is larger than Tα1 and satisfying the condition shown as (14), one state sequence is corresponding to the target, as shown in Figure 7(a), and the other state sequences are induced from noise; one of them is shown as Figure 7(b). For convenience, the ordinate value has reduced by a constant of 1×106.Figure 7
Backtracked range trajectories: (a) induced by target; (b) induced by noise.
(a)
(b)In the following text, for the state ofyk, the cumulative energy, range trajectories, angle trajectories, estimated micromotion feature parameter vector, fitting error, consistency of the extracted micromotion feature parameter vector, state transition set of target, and noise in next scan are denoted as IT(yk), IN(yk), Ryk,T, Ryk,N, Θyk,T, Θyk,N, PA^f,yk,T, PA^f,yk,N, E^f,yk,T, E^f,yk,N, Cs,yk,T, Cs,yk,N, ΓTyk+1, and ΓNyk+1, respectively.Using any two values of the range trajectory, the probable slope can be obtained, which can be utilized for the coarse slope compensation of trajectories, as shown in Figure8.Figure 8
Backtracked range trajectories with coarse slope compensation; (a) induced by target; (b) induced by noise.
(a)
(b)Based on the range trajectories of target and noise, respectively, the estimated micromotion feature parameter vectorsPA^f,yk,T and PA^f,yk,N can be obtained according to (13)–(22). At this time, the energy accumulation values IT(yk) and INyk are both smaller than threshold Tα2, and the consistencies Cs,yk,T and Cs,yk,N can not be obtained due to the inexistence of the extracted micromotion feature parameter vector at the former scan. Therefore, it is necessary to continue to carry out energy accumulation. For ITyk and IN(yk), respectively, update state transition sets ΓT(yk+1) and ΓN(yk+1) according to (27) by utilizing the obtained precession feature parameter vectors PA^f,yk,T and PA^f,yk,N.After 24 scans, forIN(yk), we can obtain the consistency Cs,yk,N=3.92 and fitting error E^f,yk,N=0.26, which does not satisfy Cs,yk<TC and E^f,yk<TE. On the contrary, for IT(yk), the consistency Cs,yk,T=0.08 and fitting error E^f,yk,T=0.05, which satisfies Cs,yk<TC and E^f,yk<TE. Therefore, we declare the presence of a target, and the micromotion feature parameter vector PA^f,yk,T can be obtained, as shown in Table 1. And the update process of the consistency, fitting error, and the size of searching range gate are shown in Figure 9.Table 1
Extracted micromotion feature parameter vector.
Parameter
|
o
a
| (m)
|
o
b
| (m)
r
0 (m)
ϕ
0 (rad)
β (°)
ϑ (°)
ω (rad/s)
v (m/s)
R
c (m)
Ture value
3.000
0.300
1.000
0.031
2.530
0.261
25.132
500.000
10
6
Estimated value
2.89
2
0.319
0.911
0
.
028
2.625
0
.
252
26.01
3
492.114
1.000023×106
Error (%)
3.60
6.33
8.90
9.67
3.75
3.45
3.52
1.58
0.00Figure 9
Update process. (a) Consistency of the extracted micromotion feature parameters. (b) Fitting error. (c) Size of searching range gate.
(a)
(b)
(c)From Figure8, we can see that the consistency of the extracted micromotion feature parameter vector is increasing (i.e., Cs,yk is decreasing) with energy accumulation, and Cs,yk,T can reduce to the value smaller than TC=0.1, while Cs,yk,N is always higher than TC=0.1. Similarly, the size of searching range gate has the same change trend with the consistency. At the same time, the fitting error increases with energy accumulation, and E^f,yk,T can stabilize at the value which is smaller than TE, while E^f,yk,N will be higher than TE. Therefore, the state sequence corresponding to the target can be selected out, and the other state sequences induced from noise can be eliminated. Furthermore, the obtained micromotion feature parameter vector shown in Table 1 is closed to the theoretical value, which illustrates the effectiveness of the proposed method.Next, the detection performance and micromotion feature extraction performance of the proposed method are compared with the traditional TBD method and the extended Hough transform method, respectively. With the proposed method, the false alarm probability and the detection probability arePFA=0.0052 and PD=0.8823. With the traditional TBD method, fixing the false alarm probability as 0.0052, the detection probability is 0.4827. Obviously, the detection performance of the proposed method is better than that of the traditional TBD method. Further, the micromotion feature extraction performance is defined as the estimated error of each micromotion feature parameter. With the extended Hough transform method, the micromotion feature parameter vector is obtained shown in Table 2. Due to the low SNR, the micromotion feature parameter vector can not be extracted effectively with the extended Hough transform method, while the high micromotion feature extraction performance can be obtained with the proposed method, shown as Table 1.Table 2
Extracted micromotion feature parameter vector in with extended Hough transform method.
Parameter
o
a (m)
o
b (m)
r
0 (m)
ϕ
0 (rad)
β (°)
ϑ (°)
ω (rad/s)
v (m/s)
R
c (m)
Ture value
3.000
0.300
1.000
0.031
2.530
0.261
25.132
500.000
106
Estimated value
6.200
0.500
0.600
0.100
1.700
0.100
13.800
653.000
1.000043
×
10
6
Error (%)
106.67
66.67
40.00
222.58
32.81
61.69
45.09
30.60
0.00Finally, the detection performance and micromotion feature extraction performance with different SNRs are considered. Fixing the false alarm probability as 0.0052, the detection probability of the proposed method is better than that of the traditional TBD method, shown in Figure10. Without loss of generality, we think the micromotion feature extraction is successful when the estimated error of each micromotion feature parameter is less than 10%. The success of micromotion feature extraction with different SNRs is shown in Figure 11. We can see that the success of micromotion feature extraction with the proposed method is much higher than that with the extended Hough transform method when SNR is lower than 11 dB. What is more important is that the micromotion feature parameters can be extracted without transmitting additional pulses, which can save the radar resources and provide real-time information for target recognition.Figure 10
Detection probability versus SNR.Figure 11
Success of micromotion feature extraction versus SNR.


## 5. Conclusions

A micromotion feature extraction method for space target based on TBD is proposed in this paper. This method establishes a feedback loop between micromotion feature extraction and TBD of target. As a result, the micromotion feature can be extracted concurrent with target detecting and tracking, and the detection performance can be improved. The steps of the method are described in detail, and some simulations are given to illustrate its effectiveness. It is noted that the cone-shaped target is taken as an example for describing and verifying the proposed feature extracting method. In fact, the proposed method is not influenced by the shape of the target. The proposed feature extracting method is suitable for the fixed scatterer model, slide scatterer model, and any other kinds of micromotion forms by using the corresponding fitting curves.However, the radar resource is limited; the resource saturation will be an important problem when there are many targets coexisting in the radar monitoring area. Thus, the reasonable and effective resources scheduling algorithms are important for exploiting the benefits of the proposed method. The related study will be presented in another independent manuscript.



---

*Source: 8723042-2017-08-22.xml*