# Direct Position Determination with Single Sensor Based on Signal Periodicity

Authors: Cheng Wang; Ding Wang; Lu Gao; Bin Yang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Due to practical limitations on size and cost, aerial vehicles generally cannot equip complicated sensors to form sensor array for target localization. In this paper, we investigate the direct position determination (DPD) of stationary source via single moving sensor. First, we analyze artificial signal structure and construct the DPD model with the frame periodicity of artificial signal. The model incorporates Doppler information extracted from both transformation frames and adjacent samples into target localization. Secondly, we consider the effect of oscillator instability and present an iterative solution for joint estimation of target location and phase noise caused by oscillator imperfection. The proposed technique fully exploits periodic structure of artificial wireless signal, which leads to significant enhancement in localization performance. Both theoretical analysis and simulations are presented to confirm its effectiveness.


---

## Body

## 1. Introduction

Passive localization of emitters is an important issue in many scientific and engineering applications. Traditional two-step localization methods extract intermediate parameter such as angle or time delay of incident signals to determine source positions [1]. In past decades, direct position determination (DPD) technique has become the focus of intensive research. DPD directly exploits sensor outputs to localize target without estimating intermediate parameter. It outperforms traditional two-step methods especially for low signal-to-noise ratio (SNR) and/or limited signal samples [2, 3].Among current DPD study, there has been considerable interest for geolocation based on single moving station. This is mainly due to the fact that both synchronization and data transmission between stations are not required when all processing is done within a single station [4]. Demissie first developed and analyzed a mathematical DPD model with moving station [5]. Based on array data, Capon iterative optimization technique and characteristic of noncircular source are introduced into DPD, respectively, for improving the positioning accuracy [6, 7]. These single-station methods are designed based on the angle information and array sensors. For moving station such as unmanned aerial vehicle, restrictions on payload dimensions and weight pose various challenges on the deployment of multichannel receiver or sensor array.To address such issue, a variety of techniques such as reflector-aided geolocation have been proposed [4]. Among existing solutions, single-station DPD with Doppler shift is regarded as the most promising one. This approach utilizes Doppler shift caused by sensor motion to determine target location and enables single sensor geolocation [8, 9]. Though effective and simple, single-station DPD based on Doppler would encounter two problems in practice. First, complex-modulated symbols would bring phase fluctuation into received signal, so it demands approaches that could separate unknown symbols and recover Doppler shift [10]. For instance, methods proposed in [11, 12] utilize perfect temporal coherence of strictly noncircular signals to remove phase ambiguity caused by complex modulation. Single-channel technique used in [13, 14] adopts phase curve fitting for joint estimation of both modulated symbols and Doppler information. However, these methods are limited to specific modulations, which makes them pale in practical application. Another inevitable challenge is the effect of oscillator instability. Slight phase noise of oscillator will lead to error in location result and must be considered [15].In this paper, we consider exploiting signal periodicity to cope with above problem and investigate the DPD algorithm based on single moving sensor. It is worth noting that most artificial signals have periodic structure such as periodic radar pulse and communication signals that are coded in frame for transformation. This periodicity allows for time difference of arrival (TDOA) localization with single moving sensor [16, 17]. Similarly, single-channel system presented in [18] virtually rotates receiving antenna with the period of target signal to provide angle information. Motivated by previous work, we analyze the signal structure and exploit signal frame periodicity to abstract Doppler information for target localization. Moreover, practical environment is considered, where phase noise induced by oscillator imperfections is not known exactly but obeys Gaussian distributions with known correlation coefficients. We establish the DPD model with oscillator phase noise. Then, an iterative solution is provided for joint estimation of target location and phase noise. Simulation results show the proposed DPD algorithm can achieve better location accuracy compared with existing methods.The main contributions of this paper are listed as follows:(1)
With the analysis of wireless signal structure, we divide each signal frame into periodic preamble and data block that contains random information symbols. Then, a novel DPD model is provided for dealing with such partially periodic signal. It incorporates the Doppler information extracted from both transformation frames and adjacent samples into target localization.(2)
We combine the established signal model with the probability distribution of oscillator phase noise. To avoid multidimensional search, an iterative optimizing scheme is then presented for updating target position and phase noise alternately. The algorithm fully exploits the location information contained in Doppler shifts and minimizes the effect of oscillator imperfections.(3)
The Cramér–Rao bound (CRB) of the source position estimation based on the received signal model is derived.The rest of the paper is organized as follows. Section2 introduces the signal model and formulates the problem. Section 3 presents the DPD method with Doppler information of received signal. Performance analysis is given via CRB in Section 4. In Section 5, simulation experiments are provided for validating the proposed estimator. Finally, conclusion is drawn in Section 6.Notation. Throughout the paper, bold lowercase letters are used for vectors, and bold upper case letters are for matrices. Superscripts T, H, and ∗ represent the transpose, conjugate transpose, and complex conjugate, respectively. Small letter j denotes the imaginary number, i.e., −1.  symbolizes absolute value and ⊗ represents the Kronecker product. Re {·} and Im {·} signify the real and imaginary parts, respectively ⊙ represents the Hadamard product.  represents the Euclidean norm. Im denotes m×m identity matrix. 1n denotes n-dimensional unit vector, and 0m is m-dimensional zero vector.


## 2. Signal Model and Problem Formulation

### 2.1. Data Model for Single Receiving Sensor

Consider the scenario as pictured in Figure1, where a flying sensor receives the direct path signal from stationary emitter. The sensor intercepts signals at K different points along the trajectory and collects each batch of data in a short interval. After downsampling, the baseband signal collected in the kth interception interval can be written as(1)rkt=βkstexpj2πfkpt+jϕkt+nkt,where p is the target coordinate vector; s(t) represents signal waveform; nkt is the additive noise; βk stands for complex fading factor; ϕkt denotes random phase noise which is mainly introduced by oscillator drift [15]; and fkp denotes the Doppler shift caused by platform motion, that is,(2)fkp=fccvkTp−qkp−qk,where fc represents signal frequency; c is the speed of light; and qk and vk denote coordinate vector and velocity vector of the receiving sensor at the kth interval, respectively.Figure 1
The scenario of source geolocation with single moving sensor.Equation (2) describes the relationship between fkp and source position, which implies that information of source position is involved in the residual frequency of rkt. Generally, sensor position qk and velocity vk are considered to be known a priori. However, in contrast with active systems where transmitted signal is under control by the receiver, st is not available for passive receivers. The complex-modulated symbols contained in st would bring phase ambiguity into rkt and restrict the Doppler shift recovery [10]. Meanwhile, the existence of ϕkt also poses challenges on the abstraction of Doppler shift.


### 2.2. Structure of Artificial Wireless Signal

To get rid of complex-modulated sequencest, we first analyze the structure of artificial wireless signal in the sequel. As shown in Figure 2, it is well known that most artificial signals are characterized by frame transmission and each frame consists of preamble block and data block [19]. Data blocks of signals generally contain audio or video data while the preamble blocks carry advanced information like data coding parameters, service labels, current date, synchronization sequence, and so on. Note that both frame structure and preamble data remain almost unchanged once communication service begins [20]. Hence, we can draw the conclusion that there exists periodicity in artificial signals and frame length TF is exactly the length of period. Such frame periodicity phenomenon also occurs in signals of periodic radar pulse or the radar rotating with constant revolution speed.Figure 2
Typical structure of wireless signals.Moreover, since pulse-shaping filters in digital modulation are real-valued, the phase ambiguity is mainly bought by information symbols ofst. As such, there exists no phase change among samples within each symbol. As demonstrated in Figure 2, adjacent samples of an oversampled signal belong to the same symbol with a high probability, so phases of these samples generally remain the same. For narrow-band signal with low baud rate, this phenomenon is frequently encountered.Both frame periodicity and the phase invariance between continuous samples has been successfully utilized in non-cooperative signal processing, such as direction of arrival estimation [19] and blind frequency recovery [21]. Based on above analysis, we list the assumptions in this study:(A1) Signal sourcestcan be divided into preamble sequence sat and data sequence sbt, that is, st=sat+sbt. sat is independent with sbt and they satisfy(3)Esatsa∗t−mTF≈σa2,m=0,1,…,M−1,Esbtsb∗t−mTF≈0,ifm≠0,σb2,ifm=0,whereTF denotes frame length and M is the number of involved frames. σa2,σb2 are the power of sat, sbt, respectively. In practice, TF can be obtained from signal standards or directly estimated from received data via TDOA measurements [15], so it is assumed to be known to the receiver.(A2) Narrow-band signal and a slow-fading environment are considered, that is,st remains unchanged in a short period of signal samples and fades independently among observation intervals. Hence, stsatisfies(4)st≈st−nTs,n=0,1,…,N−1,where integerN is the number of involved samples and Ts stands for sampling interval.(A3) The oscillator phase noiseϕkt,k=1,…,K is statistically independent from interval to interval. Within each observation interval, ejϕkt is a Gaussian stationary process that satisfies(5)Eejϕkt=1,Eejϕkt−1ejϕkt−τ−1∗=gτ,Eejϕkt−1ejϕkt−τ−1=g¯τ,whereτ denotes time delay and correlation coefficients gτ and g¯τ are assumed to be known to the receiver.Above assumptions describe the temporal coherence of artificial wireless signals. Though above discussions focus on the frame structure of artificial signals, the repetition of spreading code or cyclic prefix in orthogonal frequency division multiplexing system also meets Assumption A1. Interested readers may refer to [22] for details. Relying on above assumptions, the problem that we address now is to determine target position directly with rkt expressed in (1) and (2).


## 2.1. Data Model for Single Receiving Sensor

Consider the scenario as pictured in Figure1, where a flying sensor receives the direct path signal from stationary emitter. The sensor intercepts signals at K different points along the trajectory and collects each batch of data in a short interval. After downsampling, the baseband signal collected in the kth interception interval can be written as(1)rkt=βkstexpj2πfkpt+jϕkt+nkt,where p is the target coordinate vector; s(t) represents signal waveform; nkt is the additive noise; βk stands for complex fading factor; ϕkt denotes random phase noise which is mainly introduced by oscillator drift [15]; and fkp denotes the Doppler shift caused by platform motion, that is,(2)fkp=fccvkTp−qkp−qk,where fc represents signal frequency; c is the speed of light; and qk and vk denote coordinate vector and velocity vector of the receiving sensor at the kth interval, respectively.Figure 1
The scenario of source geolocation with single moving sensor.Equation (2) describes the relationship between fkp and source position, which implies that information of source position is involved in the residual frequency of rkt. Generally, sensor position qk and velocity vk are considered to be known a priori. However, in contrast with active systems where transmitted signal is under control by the receiver, st is not available for passive receivers. The complex-modulated symbols contained in st would bring phase ambiguity into rkt and restrict the Doppler shift recovery [10]. Meanwhile, the existence of ϕkt also poses challenges on the abstraction of Doppler shift.


## 2.2. Structure of Artificial Wireless Signal

To get rid of complex-modulated sequencest, we first analyze the structure of artificial wireless signal in the sequel. As shown in Figure 2, it is well known that most artificial signals are characterized by frame transmission and each frame consists of preamble block and data block [19]. Data blocks of signals generally contain audio or video data while the preamble blocks carry advanced information like data coding parameters, service labels, current date, synchronization sequence, and so on. Note that both frame structure and preamble data remain almost unchanged once communication service begins [20]. Hence, we can draw the conclusion that there exists periodicity in artificial signals and frame length TF is exactly the length of period. Such frame periodicity phenomenon also occurs in signals of periodic radar pulse or the radar rotating with constant revolution speed.Figure 2
Typical structure of wireless signals.Moreover, since pulse-shaping filters in digital modulation are real-valued, the phase ambiguity is mainly bought by information symbols ofst. As such, there exists no phase change among samples within each symbol. As demonstrated in Figure 2, adjacent samples of an oversampled signal belong to the same symbol with a high probability, so phases of these samples generally remain the same. For narrow-band signal with low baud rate, this phenomenon is frequently encountered.Both frame periodicity and the phase invariance between continuous samples has been successfully utilized in non-cooperative signal processing, such as direction of arrival estimation [19] and blind frequency recovery [21]. Based on above analysis, we list the assumptions in this study:(A1) Signal sourcestcan be divided into preamble sequence sat and data sequence sbt, that is, st=sat+sbt. sat is independent with sbt and they satisfy(3)Esatsa∗t−mTF≈σa2,m=0,1,…,M−1,Esbtsb∗t−mTF≈0,ifm≠0,σb2,ifm=0,whereTF denotes frame length and M is the number of involved frames. σa2,σb2 are the power of sat, sbt, respectively. In practice, TF can be obtained from signal standards or directly estimated from received data via TDOA measurements [15], so it is assumed to be known to the receiver.(A2) Narrow-band signal and a slow-fading environment are considered, that is,st remains unchanged in a short period of signal samples and fades independently among observation intervals. Hence, stsatisfies(4)st≈st−nTs,n=0,1,…,N−1,where integerN is the number of involved samples and Ts stands for sampling interval.(A3) The oscillator phase noiseϕkt,k=1,…,K is statistically independent from interval to interval. Within each observation interval, ejϕkt is a Gaussian stationary process that satisfies(5)Eejϕkt=1,Eejϕkt−1ejϕkt−τ−1∗=gτ,Eejϕkt−1ejϕkt−τ−1=g¯τ,whereτ denotes time delay and correlation coefficients gτ and g¯τ are assumed to be known to the receiver.Above assumptions describe the temporal coherence of artificial wireless signals. Though above discussions focus on the frame structure of artificial signals, the repetition of spreading code or cyclic prefix in orthogonal frequency division multiplexing system also meets Assumption A1. Interested readers may refer to [22] for details. Relying on above assumptions, the problem that we address now is to determine target position directly with rkt expressed in (1) and (2).


## 3. Proposed DPD Method

### 3.1. DPD Model with Signal Periodicity

Based on Assumption A2, a new observation vectorrt∈ℂN×1 can be constructed as follows:(6)rkt=rkt,rkt−Ts,…,rkt−N−1TsT≈βkstej2πfkpt+jϕkt,…,stej2πfkpt−N−1Ts+jϕkt−N−1TsT+nkt=βkΓ1,kakpstej2πfkpt+nkt,where N×1 noise vector nkt=nkt,nkt−Ts,…,nkt−N−1TsT. η1,k=ejϕkt,…,ejϕkt−N−1TsT describes the effect of oscillator phase noise and diagonal matrix Γ1,k=diagη1,k. The N×1 steering vector akp has the following form:(7)akp=1,e−j2πfkpTs,…,e−j2πfkpN−1TsT.According to Assumption A1,rk′t∈ℂM×1 can be similarly constructed as(8)rk′t=rkt,rkt−TF,…,rkt−M−1TFT=βksat+sbtej2πfkpt+jϕkt,…,sat−M−1TF+sbt−M−1TFej2πfkpt−M−1TF+jϕkt−M−1TFT+nk′t=βkΓ2,kbkpsat+bkp⊙sbtej2πfkpt+nk′t,where sbt=sbt,…,sbt−M−1TFT. η2,k=ejϕkt,…,ejϕkt−M−1TFT and Γ2,k=diagη2,k. The M×1 steering vector bkp has the following form:(9)bkp=1,e−j2πfkpTF,…,e−j2πfkpM−1TFT.With (6)–(9), we can construct array vector xkt∈ℂMN×1 as(10)xkt=rkTt,rkTt−TF,…,rkTt−M−1TFT=βkΓ1,kakp⊗stej2πfkpt+jϕkt,…,st−M−1TFej2πfkpt−M−1TF+jϕkt−M−1TFT+wkt=βkΓ1,kakp⊗Γ2,kbkpsat+Γ2,kbkp⊙sbtej2πfkpt+wkt,where noise vector wkt=nkTt,…,nkTt−M−1TFT. Using properties of Hadamard product and Assumption A1, we find that(11)Rb,k=EΓ2,kbkp⊙sbtΓ2,kbkp⊙sbtH=Γ2,kbkpbkHpΓ2,kH⊙EsbtsbHt=σb2Γ2,kbkpbkHpΓ2,kH⊙IM.Since the diagonal entries ofΓ2,kbkpbkHpΓ2,kH are all 1, Rb,kbecomes the diagonal matrix with σb2 as its diagonal entry, i.e., Rb,k=σb2IM. Then, with properties of Kronecker product and the independence between sat and sbt, covariance matrix of xkt leads to(12)Rk=ExktxkHt=βk2σa2Γ1,kakp⊗Γ2,kbkpΓ1,kakp⊗Γ2,kbkpH+Γ1,kakpakHpΓ1,kH⊗Rb,k+Wk=βk2σa2Γ1,kakp⊗Γ2,kbkpΓ1,kakp⊗Γ2,kbkpH+σb2Γ1,kakp⊗IMΓ1,kakp⊗IMH+Wk=βk2σa2c1,kp,ηkc1,kHp,ηk+σb2c2,kp,ηkc2,kHp,ηk+Wk,where c1,kp,ηk=Γ1,kakp⊗Γ2,kbkp. c2,kp,ηk=Γ1,kakp⊗IM. ηk=η1,kT,η2,kTT is the column vector with dimensions M+N. The noise covariance Wk=EwktwkHt is determined by temporal correlation of noise, so its elements are allowed to differ from each other. In the following, Wk is supposed to satisfy Wk=Σσk2, where σk denotes unknown noise power. MN×MN matrix Σ represents noise covariance matrix. Σ is generally known to the receiver, since it can be determined, for example, using sample statistics from a number of independent, identical experiments.With above discussions, the Markov-like estimates ofp with ηk,k=1,…,K can be obtained by fitting covariance matrix of received data to (12) in a weighted least squares sense. However, a multidimensional search is essential for estimating target position p and phase noise vector simultaneously, which has dominant computational complexity [23]. To solve the problem, an iteration optimizing scheme is adopted here. Specifically, the optimization of target position p is performed with fixed ηk,k=1,…,K in each iteration. Then, the refined estimate of ηk,k=1,…,K could be obtained with respect to the latest estimated p. The updating procedure continues iteratively until converged. In the next subsection, the above iterative solution will be derived.


### 3.2. Updation of Target Position

According to (12), Rk can be expressed in vector form as(13)γk=vecRk=σa2Φ1,kp,ηkβk2+σb2Φ2,kp,ηkβk2+vecΣσk2=Φkp,ηkςk,where vec(·) is a vector obtained by stacking columns of the argument on top of each other. Φ1,kp,ηk=vecc1,kp,ηkc1,kHp,ηk and Φ2,kp,ηk=vecc2,kp,ηkc2,kHp,ηk. Φkp,ηk=Φ1,kp,ηk,Φ2,kp,ηk,vecΣ and ςk=σa2βk2,σb2βk2,σk2T.Denoteηkb=η1,kbT,η2,kbTT as the estimation of ηk. Then, the Markov-like estimate of p with ηkbis obtained as(14)p^=argminp∑k=1Kγk−γ^kHCkγk−γ^k,where Ck denotes weighting matrix, γ^k=vecR^k, and R^k is the covariance matrix constructed by data samples of xkt. Generally, Ck is chosen to be the inverse of the asymptotic covariance of the residuals [23]. Substitution of (13) into (14) yields(15)p^=argminp∑k=1KCk0.5Φkp,ηkbςk−γ^k2,where Ck0.5 is a Hermitian square root factor of Ck. According to (15), the estimation of ςk that minimizes the cost function yields(16)ς^k=ΦkHp,ηkbCkΦkp,ηkb−1ΦkHp,ηkbCkγ^k.Substituting (16) into (15) leads to(17)p^=argminp∑k=1KCk0.5Φkp,ηkbΦkHp,ηkbCkΦkp,ηkb−1ΦkHp,ηkbCkγ^k−γ^k2=argminp∑k=1Kγ^kHCkγ^k−γ^kHCkΦkp,ηKbΦkHp,ηkbCkΦkp,ηkb−1ΦkHp,ηkbCkγ^k.As∑k=1Kγ^kHCkγ^ is constant during optimization, equation (17) becomes(18)p^=argmaxp∑k=1Kγ^kHCkΦkp,ηkbΦkHp,ηkbCkΦkp,ηkb−1ΦkHp,ηkbCkγ^k.Target positionp can be directly updated via grid search with (18). Next section investigates the procedure for phase noise vector updating.


### 3.3. Updation of Phase Noise Vector

Equation (6) can be rewritten as(19)rkt=βkΓ1,kakpstej2πfkpat+nkt=βkAkpη1,kstej2πfkpt+nkt,where Akp=diagakp. Then, the eigenvalue decomposition of its covariance matrix has the following form:(20)R1,k=ErktrkHt=λk1uk1uk1H+Uk1Λk1Uk1H,where vector uk1 denotes the eigenvector related to the largest eigenvalue λk1. Uk1 and Λk1 contain noise eigenvectors and eigenvalues of R1,k, respectively. Based on subspace theory, Akpη1,k are orthogonal to the subspace spanned by Uk1 [8]. Moreover, based on Assumption A3, the covariance matrix of η1,k and η2,k is assumed to satisfy(21)Eη1,k−1Nη1,k−1NH=G1,Eη2,k−1Mη2,k−1MH=G2,where the m, nth entry of G1 and G2 is gm−nTs and gm−nTF, respectively.Denotep^ as the vector containing estimated source location. With (20) and (21), η1,k,k=1,…,K can be updated by minimizing the following cost function:(22)fη1,1,…,η1,K=α∑k=1Kη1,kHAkHp^Uk1Uk1HAkp^η1,k+∑k=1KG1−0.5η1,k−1N2,where real constant α is generally selected according to the noise level. Based on A3, the presumed oscillator phase noise is statistically independent from interval to interval. Hence, the cost function in (22) can be decoupled as(23)fη1,k=αη1,kHAkHp^Uk1Uk1HAkp^η1,k+η1,k−1NHG1−1η1,k−1N.Taking the complex gradient with respect toη1,k and solving gives(24)η^1,k=αAkHp^Uk1Uk1HAkp^+G1−1−1G1−11N.Similarly, equation (8) can be expressed as(25)rk′t=βkΓ2,kbkpsat+Γ2,kbkp⊙sbtej2πfkpt+nk′t=βkBkpη2,ksat+Γ2,kbkp⊙sbtej2πfkpt+nk′t,where Bkp=diagbkp. Since Rb,k=σb2IM, the eigenvalue decomposition of covariance matrix of rk′t leads to(26)R2,k=Erk′trk′Ht=λk2uk2uk2H+Uk2Λk2Uk2H,where vector uk2 denotes the eigenvector related to the largest eigenvalue λk2. Uk2 and Λk2 contain noise eigenvectors and eigenvalues of R2,k, respectively. Similar to (23), η2,k,k=1,…,K can be updated by minimizing the following cost function:(27)fη2,k=αη2,kHBkHp^Uk2Uk2HBkp^η2,k+η2,k−1MHG2−1η2,k−1M.Taking the complex gradient with respect toη2,k and solving gives(28)η^2,k=αBkHp^Uk2Uk2HBkp^+G2−1−1G2−11M.The steps of above target localization scheme are listed in Table1, and related flowchart is shown in Figure 3. As shown in the table, step 3 is engaged to update target location based on received signal, while step 4 is applied for determining phase noise vector η1,k,η2,k,k=1,…,K. The iteration may be stopped by comparing estimation results between adjacent iteration. The main complexity of proposed method is involved in step 3. The computational complexity for estimating p^ via grid search is about o8JKM2N2, where J represents the number of grid points in the search. When there exists no oscillator phase noise, it merely requires steps 1∼3 to determine target location.Table 1
Implementation of the proposed method.
(1) Choose a small constante and set η1,kb=1N, η2,kb=1M.(2) Stack the received data vectorrkt,rk′t,xkt and then compute related covariance matrices and subspace.(3) Substituteη1,kb, η2,kb into (18) to estimate p^ via grid search.(4) Insertp^ into (24) and (28) for estimating η1,ka, η2,ka, k=1,…,K, respectively.(5) If∑k=1Kη1,ka−η1,kb2+η2,ka−η2,kb2 is smaller than e, regard p^ as the position estimation. Otherwise, set η1,kb=η1,ka, η2,kb=η2,ka and return to step 3.Figure 3
Flowchart of the proposed method.


## 3.1. DPD Model with Signal Periodicity

Based on Assumption A2, a new observation vectorrt∈ℂN×1 can be constructed as follows:(6)rkt=rkt,rkt−Ts,…,rkt−N−1TsT≈βkstej2πfkpt+jϕkt,…,stej2πfkpt−N−1Ts+jϕkt−N−1TsT+nkt=βkΓ1,kakpstej2πfkpt+nkt,where N×1 noise vector nkt=nkt,nkt−Ts,…,nkt−N−1TsT. η1,k=ejϕkt,…,ejϕkt−N−1TsT describes the effect of oscillator phase noise and diagonal matrix Γ1,k=diagη1,k. The N×1 steering vector akp has the following form:(7)akp=1,e−j2πfkpTs,…,e−j2πfkpN−1TsT.According to Assumption A1,rk′t∈ℂM×1 can be similarly constructed as(8)rk′t=rkt,rkt−TF,…,rkt−M−1TFT=βksat+sbtej2πfkpt+jϕkt,…,sat−M−1TF+sbt−M−1TFej2πfkpt−M−1TF+jϕkt−M−1TFT+nk′t=βkΓ2,kbkpsat+bkp⊙sbtej2πfkpt+nk′t,where sbt=sbt,…,sbt−M−1TFT. η2,k=ejϕkt,…,ejϕkt−M−1TFT and Γ2,k=diagη2,k. The M×1 steering vector bkp has the following form:(9)bkp=1,e−j2πfkpTF,…,e−j2πfkpM−1TFT.With (6)–(9), we can construct array vector xkt∈ℂMN×1 as(10)xkt=rkTt,rkTt−TF,…,rkTt−M−1TFT=βkΓ1,kakp⊗stej2πfkpt+jϕkt,…,st−M−1TFej2πfkpt−M−1TF+jϕkt−M−1TFT+wkt=βkΓ1,kakp⊗Γ2,kbkpsat+Γ2,kbkp⊙sbtej2πfkpt+wkt,where noise vector wkt=nkTt,…,nkTt−M−1TFT. Using properties of Hadamard product and Assumption A1, we find that(11)Rb,k=EΓ2,kbkp⊙sbtΓ2,kbkp⊙sbtH=Γ2,kbkpbkHpΓ2,kH⊙EsbtsbHt=σb2Γ2,kbkpbkHpΓ2,kH⊙IM.Since the diagonal entries ofΓ2,kbkpbkHpΓ2,kH are all 1, Rb,kbecomes the diagonal matrix with σb2 as its diagonal entry, i.e., Rb,k=σb2IM. Then, with properties of Kronecker product and the independence between sat and sbt, covariance matrix of xkt leads to(12)Rk=ExktxkHt=βk2σa2Γ1,kakp⊗Γ2,kbkpΓ1,kakp⊗Γ2,kbkpH+Γ1,kakpakHpΓ1,kH⊗Rb,k+Wk=βk2σa2Γ1,kakp⊗Γ2,kbkpΓ1,kakp⊗Γ2,kbkpH+σb2Γ1,kakp⊗IMΓ1,kakp⊗IMH+Wk=βk2σa2c1,kp,ηkc1,kHp,ηk+σb2c2,kp,ηkc2,kHp,ηk+Wk,where c1,kp,ηk=Γ1,kakp⊗Γ2,kbkp. c2,kp,ηk=Γ1,kakp⊗IM. ηk=η1,kT,η2,kTT is the column vector with dimensions M+N. The noise covariance Wk=EwktwkHt is determined by temporal correlation of noise, so its elements are allowed to differ from each other. In the following, Wk is supposed to satisfy Wk=Σσk2, where σk denotes unknown noise power. MN×MN matrix Σ represents noise covariance matrix. Σ is generally known to the receiver, since it can be determined, for example, using sample statistics from a number of independent, identical experiments.With above discussions, the Markov-like estimates ofp with ηk,k=1,…,K can be obtained by fitting covariance matrix of received data to (12) in a weighted least squares sense. However, a multidimensional search is essential for estimating target position p and phase noise vector simultaneously, which has dominant computational complexity [23]. To solve the problem, an iteration optimizing scheme is adopted here. Specifically, the optimization of target position p is performed with fixed ηk,k=1,…,K in each iteration. Then, the refined estimate of ηk,k=1,…,K could be obtained with respect to the latest estimated p. The updating procedure continues iteratively until converged. In the next subsection, the above iterative solution will be derived.


## 3.2. Updation of Target Position

According to (12), Rk can be expressed in vector form as(13)γk=vecRk=σa2Φ1,kp,ηkβk2+σb2Φ2,kp,ηkβk2+vecΣσk2=Φkp,ηkςk,where vec(·) is a vector obtained by stacking columns of the argument on top of each other. Φ1,kp,ηk=vecc1,kp,ηkc1,kHp,ηk and Φ2,kp,ηk=vecc2,kp,ηkc2,kHp,ηk. Φkp,ηk=Φ1,kp,ηk,Φ2,kp,ηk,vecΣ and ςk=σa2βk2,σb2βk2,σk2T.Denoteηkb=η1,kbT,η2,kbTT as the estimation of ηk. Then, the Markov-like estimate of p with ηkbis obtained as(14)p^=argminp∑k=1Kγk−γ^kHCkγk−γ^k,where Ck denotes weighting matrix, γ^k=vecR^k, and R^k is the covariance matrix constructed by data samples of xkt. Generally, Ck is chosen to be the inverse of the asymptotic covariance of the residuals [23]. Substitution of (13) into (14) yields(15)p^=argminp∑k=1KCk0.5Φkp,ηkbςk−γ^k2,where Ck0.5 is a Hermitian square root factor of Ck. According to (15), the estimation of ςk that minimizes the cost function yields(16)ς^k=ΦkHp,ηkbCkΦkp,ηkb−1ΦkHp,ηkbCkγ^k.Substituting (16) into (15) leads to(17)p^=argminp∑k=1KCk0.5Φkp,ηkbΦkHp,ηkbCkΦkp,ηkb−1ΦkHp,ηkbCkγ^k−γ^k2=argminp∑k=1Kγ^kHCkγ^k−γ^kHCkΦkp,ηKbΦkHp,ηkbCkΦkp,ηkb−1ΦkHp,ηkbCkγ^k.As∑k=1Kγ^kHCkγ^ is constant during optimization, equation (17) becomes(18)p^=argmaxp∑k=1Kγ^kHCkΦkp,ηkbΦkHp,ηkbCkΦkp,ηkb−1ΦkHp,ηkbCkγ^k.Target positionp can be directly updated via grid search with (18). Next section investigates the procedure for phase noise vector updating.


## 3.3. Updation of Phase Noise Vector

Equation (6) can be rewritten as(19)rkt=βkΓ1,kakpstej2πfkpat+nkt=βkAkpη1,kstej2πfkpt+nkt,where Akp=diagakp. Then, the eigenvalue decomposition of its covariance matrix has the following form:(20)R1,k=ErktrkHt=λk1uk1uk1H+Uk1Λk1Uk1H,where vector uk1 denotes the eigenvector related to the largest eigenvalue λk1. Uk1 and Λk1 contain noise eigenvectors and eigenvalues of R1,k, respectively. Based on subspace theory, Akpη1,k are orthogonal to the subspace spanned by Uk1 [8]. Moreover, based on Assumption A3, the covariance matrix of η1,k and η2,k is assumed to satisfy(21)Eη1,k−1Nη1,k−1NH=G1,Eη2,k−1Mη2,k−1MH=G2,where the m, nth entry of G1 and G2 is gm−nTs and gm−nTF, respectively.Denotep^ as the vector containing estimated source location. With (20) and (21), η1,k,k=1,…,K can be updated by minimizing the following cost function:(22)fη1,1,…,η1,K=α∑k=1Kη1,kHAkHp^Uk1Uk1HAkp^η1,k+∑k=1KG1−0.5η1,k−1N2,where real constant α is generally selected according to the noise level. Based on A3, the presumed oscillator phase noise is statistically independent from interval to interval. Hence, the cost function in (22) can be decoupled as(23)fη1,k=αη1,kHAkHp^Uk1Uk1HAkp^η1,k+η1,k−1NHG1−1η1,k−1N.Taking the complex gradient with respect toη1,k and solving gives(24)η^1,k=αAkHp^Uk1Uk1HAkp^+G1−1−1G1−11N.Similarly, equation (8) can be expressed as(25)rk′t=βkΓ2,kbkpsat+Γ2,kbkp⊙sbtej2πfkpt+nk′t=βkBkpη2,ksat+Γ2,kbkp⊙sbtej2πfkpt+nk′t,where Bkp=diagbkp. Since Rb,k=σb2IM, the eigenvalue decomposition of covariance matrix of rk′t leads to(26)R2,k=Erk′trk′Ht=λk2uk2uk2H+Uk2Λk2Uk2H,where vector uk2 denotes the eigenvector related to the largest eigenvalue λk2. Uk2 and Λk2 contain noise eigenvectors and eigenvalues of R2,k, respectively. Similar to (23), η2,k,k=1,…,K can be updated by minimizing the following cost function:(27)fη2,k=αη2,kHBkHp^Uk2Uk2HBkp^η2,k+η2,k−1MHG2−1η2,k−1M.Taking the complex gradient with respect toη2,k and solving gives(28)η^2,k=αBkHp^Uk2Uk2HBkp^+G2−1−1G2−11M.The steps of above target localization scheme are listed in Table1, and related flowchart is shown in Figure 3. As shown in the table, step 3 is engaged to update target location based on received signal, while step 4 is applied for determining phase noise vector η1,k,η2,k,k=1,…,K. The iteration may be stopped by comparing estimation results between adjacent iteration. The main complexity of proposed method is involved in step 3. The computational complexity for estimating p^ via grid search is about o8JKM2N2, where J represents the number of grid points in the search. When there exists no oscillator phase noise, it merely requires steps 1∼3 to determine target location.Table 1
Implementation of the proposed method.
(1) Choose a small constante and set η1,kb=1N, η2,kb=1M.(2) Stack the received data vectorrkt,rk′t,xkt and then compute related covariance matrices and subspace.(3) Substituteη1,kb, η2,kb into (18) to estimate p^ via grid search.(4) Insertp^ into (24) and (28) for estimating η1,ka, η2,ka, k=1,…,K, respectively.(5) If∑k=1Kη1,ka−η1,kb2+η2,ka−η2,kb2 is smaller than e, regard p^ as the position estimation. Otherwise, set η1,kb=η1,ka, η2,kb=η2,ka and return to step 3.Figure 3
Flowchart of the proposed method.


## 4. CRB on Position Estimation

In this section, we use the Cramér–Rao bound (CRB) on position estimation as a performance benchmark. By the properties of maximum-likelihood estimators, the asymptotic distribution ofp^ is Gaussian with covariance matrix that equals the CRB.Based on Assumption A3,η1,k, η2,k satisfy(29)Eη1,k−1Nη1,k−1NT=G¯1,Eη2,k−1Mη2,k−1MT=G¯2,where the m, nth entry of G¯1 and G¯2 is g¯m−nTs and g¯m−nTF, respectively. Substituting (21) into (29) and solving gives(30)EReη1,k−1NReTη1,k−1N=ReG1+ReG¯12,EReη2,k−1MReTη2,k−1M=ReG2+ReG¯22,EImη1,k−1NImTη1,k−1N=ReG1−ReG¯12,EImη2,k−1MImTη2,k−1M=ReG2−ReG¯22,EReη1,k−1NImTη1,k−1N=ImG¯1−ImG12,EReη2,k−1MImTη2,k−1M=ImG¯2−ImG22,EImη1,k−1NReTη1,k−1N=ImG1+ImG¯12,EImη2,k−1MReTη2,k−1M=ImG2+ImG¯22.Denote2KM+N×1 column vector η=Reη1,1T,Imη1,1T,Reη2,1T,Imη2,1T,…,Reη1,KT,Imη1,KT,Reη2,KT,Imη2,KTT. With Assumption A3 and (30), the covariance matrix of η can be rewritten as(31)Eη−η¯η−η¯H=blkdiagG1′,G2′,…,G1′,G2′,where G1′=ReG1+ReG¯1/2ImG¯1−ImG1/2ImG¯1+ImG1/2ReG1−ReG¯1/2, G2′=ReG2+ReG¯2/2ImG¯2−ImG2/2ImG¯2+ImG2/2ReG2−ReG¯2/2, and 2KM+N×1 vector η¯=1NT,0NT,1MT,0MT,…T.According to investigated DPD model and the probability distribution of the presumed phase noise vector, the log-likelihood function yields(32)fγk|ζ=−∑k=1K1σk2Ck0.5Φkp,ηkςk−γ^k2−G′−0.5η−η¯2,where ζ=p,ς,ηT, ς=ς1,ς2,…,ςKT and G′=blkdiagG1′,G2′,…,G1′,G2′.With conclusions in [24], the Fisher information matrix for all the real unknowns ζ is computed by(33)Jζ=2∑k=1K∂Φkp,ηkςk∂ζHCk−0.5Πk⊥Ck−0.5∂Φkp,ηkςk∂ζ+OOOG′−1,where Πk⊥=IMN−C−0.5Φkp,ηkΦkHp,ηkC−1Φkp,ηk−1ΦkHp,ηkC−0.5/σk2 and Odenotes zero matrix. Without loss of generality, we suppose weighting matrix C=IMN and denote(34)Ypp=∑k=1K∂Φkp,ηkςk∂pHΠk⊥∂Φkp,ηkςk∂p,Ypς=∑k=1K∂Φkp,ηkςk∂pHΠk⊥∂Φkp,ηkςk∂ς,Ypη=∑k=1K∂Φkp,ηkςk∂pHΠk⊥∂Φkp,ηkςk∂η,Yηη=∑k=1K∂Φkp,ηkςk∂ηHΠk⊥∂Φkp,ηkςk∂η,Yςς=∑k=1K∂Φkp,ηkςk∂ςHΠk⊥∂Φkp,ηkςk∂ς,Yςη=∑k=1K∂Φkp,ηkςk∂ςHΠk⊥∂Φkp,ηkςk∂η.Then, Jζ can be rewritten as(35)Jζ=2YppYpςYpηYpςTYςςYςηYpηTYςηTYηη+0.5G′−1.With matrix inversion lemma, CRB in the presence of oscillator noise is given as(36)CRBp=12Ypp−YpςYςς−1YpςT−Ζ−1,where Ζ=Ypη−YpςYςς−1YςηYηη+0.5G′−1−YςηTYςς−1Yςη−1Ypη−YpςYςς−1YςηT.


## 5. Numerical Results

In this section, simulation results are provided to demonstrate the effectiveness of proposed algorithm. Without loss of generality, we assume a planar geometry and depict both sensor trajectory and target position in Figure4. The sensor moves from [−20, 20]T (km) to [0, 20]T (km) with constant velocity v = [250, 0]T (m/s). It intercepts the transmitted signal once per 4 kilometer. The target is fixed at [20,0]T (km).Figure 4
Scenario of the target and receiving sensor.To enable DPD based on frame periodicity, target signal is assumed to have the transmission frame of 3.5 ms and each frame contains 12% periodic preamble symbols. Symbols are transmitted using QPSK modulation with 6 kBd. Signal is received with the sampling rate of 2.5 MHz and its carrier frequencyfc = 1.5 GHz. Parameters related to target signal are summarized in Table 2. When adopting DPD processing, data related to 3 signal frames (M = 3) are collected for localization during each observation interval and signal waveform is assumed to remain unchanged within 3 adjacent samples (N = 3).Table 2
Parameters of target signal.
Carrier frequency (GHz)Modulation typeBaud rate (kBd)Frame length (ms)Percentage of preamble symbols (%)Sampling rate (MHz)1.5QPSK6003.5122.5Fading factorβk,k=1,…,K is independently and randomly drawn from a complex Gaussian random generator with zero mean and the variance of 1. The correlation coefficient is supposed to decrease linearly with τ, i.e., gτ=σg2/τ, where σg2 denotes the oscillator noise power. Unless otherwise stated, it is fixed at 0 dB. Root mean square error (RMSE) is adopted to evaluate the localization estimation, which is computed by(37)RMSE=11000∑i=11000p−p^i2,where p^i is the estimation result in the ith Monte Carlo trial.Firstly, we contrast proposed method with the two-step localization algorithm (designated as two-step) presented in [10] and the array-based DPD algorithm (designated as DPD-array) introduced in [9]. DPD-array determines target localization with both Doppler information abstracted from adjacent samples and the angle estimated with sensor array. Here DPD-array is applied with a uniform linear array with 3 antennas, each spacing half-wavelength apart. The two-step method utilizes nonlinear transformation to recover Doppler information from QPSK modulated signal samples.Figure5 demonstrates position estimation results, and it shows that RMSE of the proposed method coincides with CRB when input SNR is larger than 3 dB. The result in Figure 5 indicates that the proposed method has the ability to achieve better estimation accuracy than both DPD-array and two-step methods. The improvement mainly results from the employment of signal frame periodicity. Since nonlinear transformations may amplify additive noise or produce interfering components, the achievable performance of two-step method cannot be satisfied at low SNR region.Figure 5
Estimated RMSEs versus SNR.Secondly, we evaluate the location RMSEs of tested methods with different source positions when SNR is 5 dB. As Figure6 shows, the performance of each algorithm would degrade with the increase of X coordinate, since the transmitter and observer are spaced further apart with larger X. The result implies that the proposed method is able to attain better estimation accuracy with worse location geometry.Figure 6
Estimated RMSEs versus different source positions.Thirdly, we discuss some important factors with respect to the proposed algorithm. Specifically, we change the number of preamble symbols to alter preamble percentage in Figure7 and alter the involved signal frames M in Figure 8. Other parameters including frame length remain unchanged. Figure 7 shows the influence of percentage of preamble symbols on proposed estimator. Since nonperiodic data symbols cannot produce bkp defined in (9) for localization, the proposed method achieves lower RMSE with less preamble symbols involved in Figure 7. When SNR is 5 dB, the achievable RMSE with 5% preamble symbols approaches 0.5 km. RMSEs versus the number of involved signal frames M are shown in Figure 8. The results demonstrate that the performance of proposed method degrades with less signal frames. The phenomenon is due to the fact that we exploit Doppler information for target localization. The accuracy of Doppler shift estimation is mainly determined by signal collection time, which corresponds to MTF in the proposed method.Figure 7
Estimated RMSEs versus percentage of preamble symbols.Figure 8
Estimated RMSEs with different number of involved frames.Moreover, we evaluate the proposed method when signal period is not exactly known to the receiver, i.e.,T^F=TF+Te, where Te is the error caused by oscillator or other imperfections. Figure 9 shows location estimates with distinct period errors. From the figure, we find that the performance of the proposed method degrades significantly with error larger than 1.5 μs. This phenomenon can be predicated. When period error and symbol duration are of the level, the frame periodicity assumption given in (A1) would become invalid and the Doppler extracted from transformation frames cannot be utilized for DPD. Note that the proposed method still works well when Te < 1 μs, which indicates that it provides certain robustness to such imperfections. Figure 10 shows RMSEs versus the power of oscillator phase noise when SNR = 5 dB. As phase noise brings ambiguity for localization, it can be seen from Figure 10 that RMSEs would become larger with the increase of noise power. For phase noise smaller than 5 dB, its influence can be ignored.Figure 9
Estimated RMSEs with period error.Figure 10
Estimated RMSEs versus phase noise power.


## 6. Conclusion

In this paper, a novel DPD technique has been introduced to achieve single sensor geolocation. First, we analyze the typical structure of artificial wireless signals and establish the DPD model including oscillator instability. Then, an iterative optimization scheme is presented for the joint estimation of target position and oscillator phase noise. The proposed method incorporates Doppler information caused by sensor movement into target localization and achieves good localization performance. Though analysis of the paper mainly considers single source, the proposed technique can also be applied to nonperiodic source with known sequence or multiple sources with distinct periods. DPD with frame periodicity has a number of associated research challenges, including low-complexity implementation, robust algorithm for wideband or multipath source, investigation of interference effects, and determination of the value of signal period. These will be the focus of our future work.



---

*Source: 6649287-2021-05-13.xml*