# A Mixed-Integer Solution for TDOA Based Direct Positioning of Digital Modulation Signal

Authors: Zeyu Yang; Ding Wang; Xiukun Ren; Fuquan Nie; Bin Yang
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

This paper proposes a direct position determination (DPD) method for a digital modulation signal based on time difference of arrival (TDOA) measurements. Unlike the two-step positioning process, the measurements are used to directly estimate the source position. Fully utilizing information in a transmitted waveform can improve the accuracy of source localization in a single-step method. Based on maximum likelihood (ML) estimates, when the digital modulation scheme is known to the location system, an alternating iterative DPD method is developed to locate the emitter. The objective function of the proposed algorithm takes into account the source position and transmitted symbol sequence, which makes it a mixed-integer optimization problem. In particular, ML sequence estimation is adopted and the complex envelope is restored using a genetic algorithm. Then, based on Newton’s method, an alternating iterative algorithm is proposed to update the position and symbol sequence results. Compared with the existing DPD method, the proposed algorithm gives more accurate location results for unknown waveforms. In addition, the proposed algorithm can reach the Cramér-Rao bound (CRB) for known signal waveforms, as verified using comprehensive simulations.


---

## Body

## 1. Introduction

High-precision passive localization is very important in many applications, such as signal processing, underwater acoustic, radar, navigation, traffic, and smart home technology. The conventional decentralized positioning methods require two-steps procedure [1–3]. Under the decentralized framework, firstly, the measurements, such as the angle of arrival (AOA) [4, 5], the time of arrival (TOA) [6, 7], the time difference of arrival (TDOA) [8, 9], frequency of arrival (FOA) [10, 11], the frequency difference of arrival (FDOA) [12, 13], the received signal strength (RSS) [14], and the gain ratios of arrival (GROA) [15, 16], are extracted from the received signals, and then the data collected by sensors are stored in data processing center to estimate the source position using various location methods. However, in order to achieve the optimal performance of the two-step methods, measurements must correspond to a single emitter. In contrast, direct position determination (DPD) methods concentrate the raw data from each sensor in a data processing center and then estimate the source location in one step.DPD methods [17–22] can provide the optimal source location results in various localization scenarios, because they can make full use of the signal observations [23]. Weiss [17] proposed a maximum likelihood- (ML-) based DPD method, which was applied to a single source of unknown and known waveforms, and then extended to multiple sources of known waveforms [18]. It indicates that the accurate signal information is helpful to optimizing the model and improving the positioning accuracy. In addition, the methods based on the maximum likelihood have been proven to reach the corresponding Cramér-Rao bound (CRB). It is generally believed that if the signal waveform is known, the DPD methods can accurately estimate the target position. Considering that it is difficult to obtain the original signal waveform in practical application, a series of DPD methods have been developed with specific signal properties to improve the estimation accuracy. Reuven and Weiss [24] proposed a cyclic DPD method that combines centralized processing and cyclostationary exploitation for narrowband cyclostationary signals. Furthermore, Wang et al. [25] constructed an ML-based estimator and devised an alternating iteration algorithm by using the constant modulus of the phase-modulated source. On the basis of subspace data fusion (SDF), the characteristics of orthogonal frequency division multiplexing signals were used to optimize TOA and AOA measurements [26]. Some conventional decentralized methods, such as the NC-MUSIC algorithm, use the noncircular property of complex signals to provide high-resolution DOA measurements [27]. Yin et al. [28] proposed the NC-SDF DPD method, which combines the SDF criterion and signal noncircularity. To reduce the computational complexity, a decoupled ML-based DPD algorithm for noncircular sources is developed based on Doppler shifts and AOA measurements [29]. Nevertheless, DPD methods can apply the product of various signal properties of the transmitted signals and improve the localization accuracy in the associated scenarios.In the harsh electromagnetic environment, the waveform of the received signal is usually difficult to determine. Unfortunately, the abovementioned DPD methods did not consider the digital modulation signal, especially the modulation scheme known to the location system. Nowadays, digital modulation signals (e.g., binary phase shift keying (BPSK) and quadrature phase shift keying (QPSK), and quadrature amplitude modulation (QAM) signals) are widely used in modern communication systems. The Global Positioning System (GPS) also provides high-precision position perception through satellite constellations [30] using specific waveforms, so typical examples include communication transmitters, such as airplanes, drones, or transport ships. The extraction of the received signal features can greatly improve positioning accuracy. If the digital modulation mode is known, the baseband complex envelope of the received signal can be reconstructed from the estimated data sequence. Then, the positioning problem of digital modulation signal can be transformed into the combination of position estimation and data sequence estimation. Because the estimated parameters include the source position and the data sequence of the transmitted waveforms, the cost function of the proposed algorithm is actually a mixed-integer optimization problem [31]. Note that the existing DPD methods usually use exhaustive search methods, which is impractical when the search dimension is determined by the length of the data sequence. The previous optimization methods are inevitable for integer optimization problems. Therefore, in this case, it is important to develop the mixed-integer DPD method, and the improvement of the prior information of digital modulation schemes is still unclear.In this paper, we propose a mixed-integer DPD method for the digital modulation signal source based on original signal observations, which contain the TDOA measurements. First, the received signals are modeled with a known digital modulation scheme. Different from the existing DPD methods of processing observations in the frequency domain, the proposed method employs the sinc function [32] to reconstruct the time-domain model based on TDOA measurements. Then, a nonlinear cost function that considers the source position and data sequence is established by using the ML criterion. In addition, an iterative technique is applied to update the data sequence and source position alternatively to reduce computational complexity. A modified genetic algorithm (GA) [33] is adopted in the mixed-integer solution. In order to evaluate the performance of the proposed method, the CRBs with the known and unknown waveforms are derived. However, as described in detail later, the proposed DPD method for a digital modulation signal cannot be regarded as a native fusion of the data sequence and source position estimation.The main contributions are summarized as follows:(1)
Firstly, the time-domain signal model based on TDOA measurements is established in this paper. In addition, an ML-based cost function is formulated for digital modulation signals. In order to improve the positioning accuracy, the digital modulation scheme is known to the positioning systems(2)
To solve the mixed-integer optimization problem, a computationally sufficient DPD method was developed to locate the emitter source based on a modified genetic algorithm and Newton’s method. The proposed solution can guarantee the optimal position estimation, rather than simply combining the process of estimating the source position and the transmission sequence(3)
The CRB on the estimated variance of unknown and a known signal waveform is derived. They are used to evaluate the performance of the proposed method. Simulation results show that the proposed method outperforms the existing DPD method for unknown waveforms and its performance can reach the CRB for the known waveform caseThe rest of this paper is organized as follows. In Section2, a digital modulation signal model is described to formulate the localization problem based on TDOA measurements, and the CRBs with unknown and known signal waveform are derived for evaluation. In Section 3, an alternating iterative DPD method is developed for data sequence estimation and source position estimation. Numerical simulations are conducted to compare the performance of the proposed algorithm with existing methods in Section 4 and conclusions are drawn in Section 5.


### 1.1. Notations

We use uppercase and lowercase boldface to represent the matrices and vectors, respectively. In this paper,ON×n is a N×n matrix with all-zero entries and IN is a N dimensional identity matrix; AT, AH, A‐1, and A† are the transpose, the conjugate transpose, the inverse, and the Moore-Penrose inverse of matrix A; a denotes the Euclidean norm of vector a; <·>i,j is the element at ith row and jth column of its argument and <·>i denotes the ith element; vec·, E·, and · are the vectorization, expectation, and rounded down operation; ⊗ is the Kronecker product.


## 1.1. Notations

We use uppercase and lowercase boldface to represent the matrices and vectors, respectively. In this paper,ON×n is a N×n matrix with all-zero entries and IN is a N dimensional identity matrix; AT, AH, A‐1, and A† are the transpose, the conjugate transpose, the inverse, and the Moore-Penrose inverse of matrix A; a denotes the Euclidean norm of vector a; <·>i,j is the element at ith row and jth column of its argument and <·>i denotes the ith element; vec·, E·, and · are the vectorization, expectation, and rounded down operation; ⊗ is the Kronecker product.


## 2. Problem Formulation

### 2.1. Signal Model

We consider a two-dimensional (D=2) scenario shown in Figure 1 including an unknown emitter at p=xyT∈ℝD×1 and a location system equipped with M known sensors, which are located at um=xmymT∈ℝD×1. Suppose the emitter transmits a digital signal with a known modulation scheme, but the symbol sequence is unknown to the location system.Figure 1
Location model: the emitter source is located byM stations equipped with sensors.Without loss of generality, all sensors are synchronized in time and the intercepted signals only propagate in the line of sight (LOS). Thus, the received waveform of themth sensor is given by
(1)rmt=bmst−τm−t0+wmt,for m=1,2,⋯,M, where τm=p−um/c represents the signal propagation time between the emitter and the mth sensor, and t0 is transmitted time. c is the signal propagation speed. bm is the mth path loss coefficient and wmt is a complex additive white Gaussian noise with zero mean and the power density σn2. By denoting with T the symbol period and the baseband chip waveform ht, the baseband complex envelope of the signal st can be written as
(2)st=∑l=1Lxl⋅ht−lT,where xl∈0,1 represents the lth symbol of the unknown transmitted sequence. In this paper, the chip waveform ht is a known root-raised cosine filter. If the modulation scheme is known to the location system, the relationship between the transmitted waveform and the data sequence x=x1x2⋯xLT can be regarded as a functional mapping relation sx. Based on a certain type of digital modulation scheme, the discrete sequence is converted into a signal waveform using a root-raised cosine filter.In passive localization problem, the TDOA measurements contain the source position information but it is difficult to obtain the transmitted time. By setting the first sensor as the reference, the TDOA between themth sensor and the first sensor is given by the following formula
(3)τm,1=τm−τ1=1cp−um−p−u1,for m=1,2,⋯,M. Then, the received signal can be rewritten as
(4)rmt=bms~t−τm,1+wmt,where s~t=st−τ1−t0. Using the properties of the Fourier transform that
(5)Fe−jωτm,1=12π∫−ππe−jωτm,1ejωtdω=sinct−τm,1,and defining the function sincν=sinπν/πν, the delayed signal in (4) is given by
(6)s~t−τm,1=∫−∞∞sincν−τm,1s~t−νdν.After being sampled att=nTs in each interval, where Ts denotes the sampling period, the received signal of sensor m can be rewritten as
(7)rmn=bms~nTs−τm,1+wmn,and then the delayed signal waveform in (6) is transformed into
(8)s~nTs−τm,1≜∑k=−∞∞sinck−τm,1/TssnTs−kTs.In practice, the coefficientK cannot be infinite, so we refer K to a limited non-infinite integer K which is larger than τm,1/Ts, so we have
(9)s~nTs−τm,1≈∑k=−KKsinck−τm,1/TssnTs−kTs.To minimize the truncation errors, using a higher orderK will lead to higher modeling accuracy. We further composite all the N observation samples into a vector and yield
(10)r¯=As~x+w¯in which
(11)r¯=r1Tr2T⋯rMTT∈ℝMN×1A=b1A1Tb2A2T⋯bMAMTT∈ℝMN×Ns~=sTs−τ1−t0s2Ts−τ1−t0⋯sNTs−τ1−t0T∈ℝN×1w¯=w1Tw2T⋯wMTT∈ℝMN×1..With the definitionSinck=sinck−τm,1/Ts, we have the matrix
(12)Am≜Sinc0⋯Sinc−k⋯Sinc−K0⋯⋯0⋮⋮⋮⋱0⋮Sinck⋯Sinc0⋯Sinc−kSinc−K⋱⋮⋮⋮⋱⋮⋱0SincK⋯Sinck⋯Sinc0Sinc−kSinc−K⋱⋮0⋱⋮⋱⋮⋱⋮⋱⋮⋱0⋮SincKSinckSinc0⋯Sinc−k⋯Sinc−K0⋱⋮⋱⋮⋮⋮SincKSinck⋯Sinc0⋯Sinc−k⋮0⋱⋮⋮⋮0⋯⋯0SincK⋯Sinck⋯Sinc0,for m=2,⋯,M and A1=IN. Multiplied with the special matrix Am, the transmitted waveform s~ is delayed with τm,1. Therefore, the waveformssmt is produced with the s1t and matrix Am with corresponding time delay. Different from the conventional methods in the frequency domain, Am can realize the fractional time delay, which is much more precise and direct to model the received signals.


### 2.2. Comparison with the CRBs

CRB is the lowest possible variance that can be achieved by an unbiased estimator [34]. Here, the CRB is needed to be deduced to evaluate the performance of the proposed method. Let us compose all the unknown parameters into a real vector. The unknown parameters are contained in the vector
(13)ρ=xTpTRebTImbTT.It can be seen that in the proposed method, the transmitted digital modulation scheme is called prior information, and it is inappropriate to differentiate the cost function with respect to the integer data sequencex together with p and b. Therefore, the CRBs are considered to be derived for both known and unknown waveforms, which means the vector ρ can be transformed to
(14)φ=Res~xTIms~xTpTRebTImbTT.And then the CRB of unknown waveform can be devised first. In order to reduce the calculation of the inverse of high dimensional Fisher information matrix (FIM), a novel proposition is introduced to derive the CRB, which is based on the following proposition.Proposition 1.
For a real vectorx whose CRB is CRBx, if we define a new real vector y=Jx, where J is an invertible matrix, the CRB of vector y can be obtained by CRBy=J⋅CRBx⋅JT [35].Letφ^ be an unbiased estimate of the parameter vector φ. The corresponding CRB bounds the mean square error of any unbiased estimator of φ and the mean square error (MSE) matrix of φ^ satisfies the following information inequality
(15)Eφ−φ^φ−φ∧T≻¯Jφ−1,where Jφ is the FIM of φ defined as
(16)Jφ=−E∂2lnfr¯φ∂φ∂φT.The log-likelihood ratio of the corresponding signal model (10) is given by
(17)lnfr¯φ∝1σn2∑m=1Mrm−bmAms~x2=1σn2r¯−As~x2,up to an additive constant, and thus we can derive the CRBs of the unknown waveforms. After performing differentiation, the associated CRB matrix for φ can be computed by
(18)CRBφu−1=Jφu=2σn2Re∂As~x∂φTH∂As~x∂φT=2σn2ReαβHαβ,where
(19)α=∂As~∂Res~xT∂As~∂Ims~xT,β=∂As~/∂ωoT=∂As~/∂pT∂As~/∂RebT∂As~/∂ImbT,in which ωo=pTRebTImbTT. It can be seen that the parameter vector determines the dimension of the CRB matrix. Therefore, we can reduce computational complexity of matrix inversion by using Proposition 1. Using Proposition 1 and Schur-Complementary theorem [36], the FIM of ωo can be reformulated as
(20)Jωou=JppuJpbuJbpuJbbu=2σn2ReβHΠα⊥β.Then, the CRB for the unknown waveform is given by(21)CRBpu=Jppu−JpbuJbbu−1Jbpu−1.In the case of known waveforms, the emitter positionp and path loss coefficient b are retained in ωo=pTRebTImbTT. Based on the log-likelihood function (17), the FIM of unknown parameter ωo can be partitioned into sub-matrices, which reduces to
(22)CRBωok−1=Jωok=2σn2Re∂A¯spb∂ωoTH∂A¯spb∂ωoT=JppkJpbkJbpkJbbk.Then, we can compute thep corner of the CRB matrix as follow
(23)CRBpk=Jppk−JpbkJbbk−1Jbpk−1.The diagonals ofCRBpu and CRBpk both contain D values that are the variance bounds for the parameters in p. The trace of CRBpu and CRBpk are compared with the localization mean square errors of the proposed estimator in Section 4.


## 2.1. Signal Model

We consider a two-dimensional (D=2) scenario shown in Figure 1 including an unknown emitter at p=xyT∈ℝD×1 and a location system equipped with M known sensors, which are located at um=xmymT∈ℝD×1. Suppose the emitter transmits a digital signal with a known modulation scheme, but the symbol sequence is unknown to the location system.Figure 1
Location model: the emitter source is located byM stations equipped with sensors.Without loss of generality, all sensors are synchronized in time and the intercepted signals only propagate in the line of sight (LOS). Thus, the received waveform of themth sensor is given by
(1)rmt=bmst−τm−t0+wmt,for m=1,2,⋯,M, where τm=p−um/c represents the signal propagation time between the emitter and the mth sensor, and t0 is transmitted time. c is the signal propagation speed. bm is the mth path loss coefficient and wmt is a complex additive white Gaussian noise with zero mean and the power density σn2. By denoting with T the symbol period and the baseband chip waveform ht, the baseband complex envelope of the signal st can be written as
(2)st=∑l=1Lxl⋅ht−lT,where xl∈0,1 represents the lth symbol of the unknown transmitted sequence. In this paper, the chip waveform ht is a known root-raised cosine filter. If the modulation scheme is known to the location system, the relationship between the transmitted waveform and the data sequence x=x1x2⋯xLT can be regarded as a functional mapping relation sx. Based on a certain type of digital modulation scheme, the discrete sequence is converted into a signal waveform using a root-raised cosine filter.In passive localization problem, the TDOA measurements contain the source position information but it is difficult to obtain the transmitted time. By setting the first sensor as the reference, the TDOA between themth sensor and the first sensor is given by the following formula
(3)τm,1=τm−τ1=1cp−um−p−u1,for m=1,2,⋯,M. Then, the received signal can be rewritten as
(4)rmt=bms~t−τm,1+wmt,where s~t=st−τ1−t0. Using the properties of the Fourier transform that
(5)Fe−jωτm,1=12π∫−ππe−jωτm,1ejωtdω=sinct−τm,1,and defining the function sincν=sinπν/πν, the delayed signal in (4) is given by
(6)s~t−τm,1=∫−∞∞sincν−τm,1s~t−νdν.After being sampled att=nTs in each interval, where Ts denotes the sampling period, the received signal of sensor m can be rewritten as
(7)rmn=bms~nTs−τm,1+wmn,and then the delayed signal waveform in (6) is transformed into
(8)s~nTs−τm,1≜∑k=−∞∞sinck−τm,1/TssnTs−kTs.In practice, the coefficientK cannot be infinite, so we refer K to a limited non-infinite integer K which is larger than τm,1/Ts, so we have
(9)s~nTs−τm,1≈∑k=−KKsinck−τm,1/TssnTs−kTs.To minimize the truncation errors, using a higher orderK will lead to higher modeling accuracy. We further composite all the N observation samples into a vector and yield
(10)r¯=As~x+w¯in which
(11)r¯=r1Tr2T⋯rMTT∈ℝMN×1A=b1A1Tb2A2T⋯bMAMTT∈ℝMN×Ns~=sTs−τ1−t0s2Ts−τ1−t0⋯sNTs−τ1−t0T∈ℝN×1w¯=w1Tw2T⋯wMTT∈ℝMN×1..With the definitionSinck=sinck−τm,1/Ts, we have the matrix
(12)Am≜Sinc0⋯Sinc−k⋯Sinc−K0⋯⋯0⋮⋮⋮⋱0⋮Sinck⋯Sinc0⋯Sinc−kSinc−K⋱⋮⋮⋮⋱⋮⋱0SincK⋯Sinck⋯Sinc0Sinc−kSinc−K⋱⋮0⋱⋮⋱⋮⋱⋮⋱⋮⋱0⋮SincKSinckSinc0⋯Sinc−k⋯Sinc−K0⋱⋮⋱⋮⋮⋮SincKSinck⋯Sinc0⋯Sinc−k⋮0⋱⋮⋮⋮0⋯⋯0SincK⋯Sinck⋯Sinc0,for m=2,⋯,M and A1=IN. Multiplied with the special matrix Am, the transmitted waveform s~ is delayed with τm,1. Therefore, the waveformssmt is produced with the s1t and matrix Am with corresponding time delay. Different from the conventional methods in the frequency domain, Am can realize the fractional time delay, which is much more precise and direct to model the received signals.


## 2.2. Comparison with the CRBs

CRB is the lowest possible variance that can be achieved by an unbiased estimator [34]. Here, the CRB is needed to be deduced to evaluate the performance of the proposed method. Let us compose all the unknown parameters into a real vector. The unknown parameters are contained in the vector
(13)ρ=xTpTRebTImbTT.It can be seen that in the proposed method, the transmitted digital modulation scheme is called prior information, and it is inappropriate to differentiate the cost function with respect to the integer data sequencex together with p and b. Therefore, the CRBs are considered to be derived for both known and unknown waveforms, which means the vector ρ can be transformed to
(14)φ=Res~xTIms~xTpTRebTImbTT.And then the CRB of unknown waveform can be devised first. In order to reduce the calculation of the inverse of high dimensional Fisher information matrix (FIM), a novel proposition is introduced to derive the CRB, which is based on the following proposition.Proposition 1.
For a real vectorx whose CRB is CRBx, if we define a new real vector y=Jx, where J is an invertible matrix, the CRB of vector y can be obtained by CRBy=J⋅CRBx⋅JT [35].Letφ^ be an unbiased estimate of the parameter vector φ. The corresponding CRB bounds the mean square error of any unbiased estimator of φ and the mean square error (MSE) matrix of φ^ satisfies the following information inequality
(15)Eφ−φ^φ−φ∧T≻¯Jφ−1,where Jφ is the FIM of φ defined as
(16)Jφ=−E∂2lnfr¯φ∂φ∂φT.The log-likelihood ratio of the corresponding signal model (10) is given by
(17)lnfr¯φ∝1σn2∑m=1Mrm−bmAms~x2=1σn2r¯−As~x2,up to an additive constant, and thus we can derive the CRBs of the unknown waveforms. After performing differentiation, the associated CRB matrix for φ can be computed by
(18)CRBφu−1=Jφu=2σn2Re∂As~x∂φTH∂As~x∂φT=2σn2ReαβHαβ,where
(19)α=∂As~∂Res~xT∂As~∂Ims~xT,β=∂As~/∂ωoT=∂As~/∂pT∂As~/∂RebT∂As~/∂ImbT,in which ωo=pTRebTImbTT. It can be seen that the parameter vector determines the dimension of the CRB matrix. Therefore, we can reduce computational complexity of matrix inversion by using Proposition 1. Using Proposition 1 and Schur-Complementary theorem [36], the FIM of ωo can be reformulated as
(20)Jωou=JppuJpbuJbpuJbbu=2σn2ReβHΠα⊥β.Then, the CRB for the unknown waveform is given by(21)CRBpu=Jppu−JpbuJbbu−1Jbpu−1.In the case of known waveforms, the emitter positionp and path loss coefficient b are retained in ωo=pTRebTImbTT. Based on the log-likelihood function (17), the FIM of unknown parameter ωo can be partitioned into sub-matrices, which reduces to
(22)CRBωok−1=Jωok=2σn2Re∂A¯spb∂ωoTH∂A¯spb∂ωoT=JppkJpbkJbpkJbbk.Then, we can compute thep corner of the CRB matrix as follow
(23)CRBpk=Jppk−JpbkJbbk−1Jbpk−1.The diagonals ofCRBpu and CRBpk both contain D values that are the variance bounds for the parameters in p. The trace of CRBpu and CRBpk are compared with the localization mean square errors of the proposed estimator in Section 4.


## 3. Location Algorithm: Digital Signal with Known Modulation and Unknown Symbol Sequences

In this section, a mixed-integer DPD optimization model is developed for the digitally modulated signal. For previous works, the estimators for the cases of known and unknown waveforms are investigated in [18]. When the transmitted signals are known or the waveforms can be restored based on the symbols and the known modulation scheme, the transmitted waveforms can be regarded as prior information, which can be used to improve the accuracy of position estimation. Based on the concept of restoring transmitted waveforms with a known modulation scheme, an efficient iterative algorithm is proposed for estimating the source position.Assume thatw¯ is white complex Gaussian noise with zero mean and variance σn2, the parameters ρ=xTpTbTT need to be estimated. To decouple the uninterested parameters b, the vector r¯ of observation samples can be rewritten as
(24)r¯=A¯sb+w¯where
(25)A¯s=A1s~ON×1ON×1ON×1ON×1A2s~ON×1ON×1⋮⋮⋮ON×1⋯ON×1AMs~.After omitting the constant terms, we can minimize the log-likelihood function (17) with respect to the unknown path loss coefficient b, which is given by
(26)b^=A¯sHA¯s−1A¯sHr¯.Submitting (26) into the objective function (17), we have
(27)f2p,x=2r¯H⋅ΠA¯s⊥⋅r¯,where
(28)ΠA¯s⊥=IMN−A¯sA¯sHA¯s−1A¯sH.Note that the optimization model is a mixed-integer nonlinear programming problem, which can be formulated as follows(29)minp,xf2p,xs.t.p∈ℝD×1,xl∈0,1,l=1,⋯,L.It can be seen that ML sequence estimation (MLSE) has been used to obtain the data sequence of the transmitted signal. If the estimated parameterp acts as a constant vector, the abovementioned cost function (29) can be regarded as the ML estimator of the data sequence, and it is a standard integer nonlinear optimization problem. Based on the idea of a decoupled iteration solution, we can update the source position and the data sequence alternatively. We divide the unknown parameters ρ=xTpTbTT into two groups, x and p,b. The two groups are optimized in turn, until the cost function finally converges to the minimum. The optimization model is divided into the following two stages.


### 3.1. First Stage: Sequence Estimation

First of all, by minimizing the cost function (29) with an appropriate initial estimate p^ using a method for unknown waveforms, MLSE can be simply derived. Then, the previous mixed-integer nonlinear programming problem can be transformed into the following pure-integer nonlinear optimization problem
(30)minxf2x,p^s.t.xl∈0,1,l=1,⋯,L.The above sequence estimator is the first step in the process of alternating iteration. Because the length of the transmitted data sequence is finite, the finite feasible set contains2L members, and the number of feasible solutions increases nonlinearly with the length of data sequence. Here, a feasible set can be regarded as a population and a feasible solution as an individual. An improved real-coded genetic algorithm utilized in the GA toolbox in MATLAB [33] can be used to solve the integer optimization problem. Without loss of generality, the fitness functions in the GA are derived based on the original objective function. Through mathematical deduction, the objective function (30) can be rewritten as
(31)x^MLSE=argmaxxl∈0,1,l=1,⋯,Lffitnessx,p^=argmaxxl∈0,1,l=1,⋯,Lr¯HA¯sx,p^A¯sHx,p∧A¯sx,p∧−1A¯sHx,p^r¯,which is a standard integer nonlinear optimization problem. Different from the conventional algorithm for mixed-integer optimization problems, the elements in each sequence are limited to 0,1 without equality and inequality constraints in the proposed estimator, but the fitness function is equal to the original objective function.The extended Laplace crossover operator and the power mutation operator are applied, and the elite individuals in the last generation are used to make a mating pool. With the Laplace crossover operator, two offspring are generated from two elite members,x=x1,x2,⋯,xLT and y=y1,y2,⋯,yLT, which are called parents here. The following Laplace distribution is used
(32)ℓl=u−vlogλl,γl≤0.5u+vlogλl,γl>0.5,where λl,γl∈0,1 are uniform random numbers. Because of the restriction of a discrete sequence, these parameters u and v are set to zero and unity, respectively. Then, the elements of the offspring w=w1,w2,⋯,wLT and z=z1,z2,⋯,zLT can be obtained as
(33)wl=xl+ℓlxl−yl,zl=yl+ℓlxl−yl.Then, after applying crossover and mutation to the old population, the truncation procedure is carried out to ensure that the integer restriction is satisfied. The elementxl is equal to either 0 or 1 with probability 0.5 if xl is not an integer.The complexity of calculation depends on the signal duration and the oversampling method. Considering the running time in real-world application, the number of samples may be limited.


### 3.2. Second Stage: Position Estimation

When the updated sequence converges, we move to the second stage. After estimating the sequencex^MLSE with an imprecise source position estimation, it is necessary to update the position estimation p^ by the following estimator
(34)p^=argminp∈ℝD×1,bf2p,x^=argminp∈ℝD×1,br¯H⋅ΠA¯sp,x^⊥⋅r¯,where
(35)ΠA¯sp,x^⊥=IMN−A¯sp,x^A¯sHp,x∧A¯sp,x∧−1A¯sHp,x^.Because the traditional grid search requires repeated computations, the estimator (34) can be minimized using Newton’s iteration method. The updated vector ωo=pTRebTImbTT is given by
(36)ωq+1o=ωqo−μqH−1ωqohωqo,where μq∈0,1 is the step size of the qth iteration. hωo and Hωo are the gradient vector and the Hessian matrix of (34), respectively.Applying the first-order derivation operator of the orthogonal projection matrixΠA¯sp,x^⊥, we can obtain the gradient vector hωo. Therefore, Hessian matrix Hωo which is essential to Newton’s-type method can be derived based on equation (38). As a result, the approximation has no effect on the numerical value of the estimation. Differentiating f2p,x^ with respect to ωo, the ith element of the gradient vector hωo is
(37)∂f2p,x^∂<ωo>i=Rer¯H⋅−ΠA¯sp,x^⊥⋅∂A¯sp,x^∂<ωo>i⋅A¯s†p,x^−ΠA¯sp,x∧⊥⋅∂A¯sp,x∧∂<ωo>i⋅A¯s†p,x∧H⋅r¯.Following some algebraic manipulations, thei,j element in the Hessian matrix Hωo is given by
(38)∂2f2p,x^∂<ωo>i∂<ωo>j≈2Rer¯H⋅A¯s†p,x∧H⋅∂A¯sHp,x^∂<ωo>i⋅ΠA¯sp,x^⊥⋅∂A¯sp,x^∂<ωo>j⋅A¯s†p,x^⋅r¯.With the matrix transformationvecXYZ=ZT⊗X⋅vecY, the gradient vector hωo can be expressed as
(39)hωo=−ReVA¯sp,x^⋅A¯s†p,x^r¯⊗r¯HΠA¯sp,x∧⊥T+ΠA¯sp,x∧⊥Hr¯⊗r¯HA¯s†p,x∧HT,where
(40)VA¯sp,x^=vec∂A¯sp,x∧∂pvec∂A¯sp,x∧∂Rebvec∂A¯sp,x∧∂ImbT.Based on the property thatΠA¯sp,x^⊥=ΠA¯sp,x^⊥⋅ΠA¯sp,x^⊥, the Hessian matrix Hωo can be formulated as follow
(41)Hωo≈2ReVA¯sp,x^⋅A¯s†p,x^r¯⊗ΠA¯sp,x^⊥A¯s†p,x∧r¯T⊗ΠA¯sp,x^⊥⋅VA¯sp,x^T,which is essential to the Newton’s-type method. Then, we can compute the gradient vector hωo and Hωo for Newton’s iteration algorithm in (36).The main formulas are given above. With appropriate step-size factors [37], the updated result can converge to the corresponding optimal result in several iterations, which is verified by the experiments in Section 4. The step-size factor μq can be adaptively determined to follow the strong Wolfe conditions [38], which means
(42)f2ωqo+μq⋅H−1ωqo⋅hωqo,x^≤f2ωqo,x^+μq⋅hTωqo⋅H−1ωqo⋅h1ωqo,h1Tωqo+μq⋅H−1ωqo⋅hωqo⋅H−1ωqo⋅hωqo≤c1hTωqo⋅H−1ωqo⋅hωqo,with 0<c1<1.Note that the data sequence is initially estimated using a relatively inaccurate position estimation. A good initial estimate is very important for subsequent iteration. In this paper, the position estimate for the first MLSE in (31) is set as the initial estimate, which is accurate enough for the proposed convergence. The procedure of the proposed DPD method is shown in Table 1, and for clarity, the flow diagram can be seen in Figure 2.Table 1
The direct localization algorithm for digital modulation signals.
Require: Obtain the initial guess of emitter positionp^0; the received signal of each sensor rm,m=1,2,⋯,M; thresholds ε and δ (sufficiently small and positive values, which are 10-6 in this paper).
A. First-stage processing (sequence estimation):
(i) Apply crossover and mutation toXelitei and generate P random discrete data sequences Xi=x1i,x1i,⋯,xP−pi,Xelitei as the ith population. The subset Xelitei contains the elite individuals in the last generation
(ii) Update the maximum of fitness value and the corresponding sequencex^=xbesti in Xi
(iii) Update the searching generationi⟵i+1
B. Second-stage processing (position estimation):
(i) Compute the gradient vectorhω and the hessian matrix Hω using (39) and (41)
(ii) Update the iterative resultωq+1 via (36), then q⟵q+1
(iii) Ifhωq≤ε, obtain p^=ωq1:D
(iv) Ifp^−p^0≤δ, output p^ as the estimated position, otherwise, reset p^0=p^ and return to the first-stage
Output: The estimated positionp^.Figure 2
Flow diagram of the proposed DPD method.


### 3.3. Computational Complexity Analysis

Table2 summarizes the numerical complexity, where Kstep1, Kstep2, and Kiter are the amount of looping in the procedure, respectively.Table 2
Computational complexity of each iteration.
Items infirst stage
The number of multiplications
Items insecond stage
The number of multiplications
f
fitness
x
,
p
^
M
2
N
2
+
1
M
N
ω
q
+
1
o
D
+
2
M
D
+
2
M
+
1
f
2
p
,
x
^
M
2
N
2
+
1
M
N
H
ω
o
4
M
3
N
2
+
2
M
2
N
2
h
ω
o
2
M
N
D
+
2
M
2
M
N
+
4
M
+
1
Inverse ofHωo
D
+
2
M
3
−
2
D
+
2
M
2
+
6
D
+
2
M
−
1
Total offirst stage
S
1
=
K
step
1
M
2
N
2
+
1
M
N
Total ofsecond stage
S
2
=
K
step
2
M
3
N
3
+
4
M
3
N
2
+
2
M
2
N
2
+
M
N
+
2
D
+
2
M
2
M
2
N
2
+
4
M
2
N
+
M
N
+
D
+
2
M
3
−
D
+
2
M
2
+
7
D
+
2
M
−
1
Total
K
iter
S
1
+
S
2


## 3.1. First Stage: Sequence Estimation

First of all, by minimizing the cost function (29) with an appropriate initial estimate p^ using a method for unknown waveforms, MLSE can be simply derived. Then, the previous mixed-integer nonlinear programming problem can be transformed into the following pure-integer nonlinear optimization problem
(30)minxf2x,p^s.t.xl∈0,1,l=1,⋯,L.The above sequence estimator is the first step in the process of alternating iteration. Because the length of the transmitted data sequence is finite, the finite feasible set contains2L members, and the number of feasible solutions increases nonlinearly with the length of data sequence. Here, a feasible set can be regarded as a population and a feasible solution as an individual. An improved real-coded genetic algorithm utilized in the GA toolbox in MATLAB [33] can be used to solve the integer optimization problem. Without loss of generality, the fitness functions in the GA are derived based on the original objective function. Through mathematical deduction, the objective function (30) can be rewritten as
(31)x^MLSE=argmaxxl∈0,1,l=1,⋯,Lffitnessx,p^=argmaxxl∈0,1,l=1,⋯,Lr¯HA¯sx,p^A¯sHx,p∧A¯sx,p∧−1A¯sHx,p^r¯,which is a standard integer nonlinear optimization problem. Different from the conventional algorithm for mixed-integer optimization problems, the elements in each sequence are limited to 0,1 without equality and inequality constraints in the proposed estimator, but the fitness function is equal to the original objective function.The extended Laplace crossover operator and the power mutation operator are applied, and the elite individuals in the last generation are used to make a mating pool. With the Laplace crossover operator, two offspring are generated from two elite members,x=x1,x2,⋯,xLT and y=y1,y2,⋯,yLT, which are called parents here. The following Laplace distribution is used
(32)ℓl=u−vlogλl,γl≤0.5u+vlogλl,γl>0.5,where λl,γl∈0,1 are uniform random numbers. Because of the restriction of a discrete sequence, these parameters u and v are set to zero and unity, respectively. Then, the elements of the offspring w=w1,w2,⋯,wLT and z=z1,z2,⋯,zLT can be obtained as
(33)wl=xl+ℓlxl−yl,zl=yl+ℓlxl−yl.Then, after applying crossover and mutation to the old population, the truncation procedure is carried out to ensure that the integer restriction is satisfied. The elementxl is equal to either 0 or 1 with probability 0.5 if xl is not an integer.The complexity of calculation depends on the signal duration and the oversampling method. Considering the running time in real-world application, the number of samples may be limited.


## 3.2. Second Stage: Position Estimation

When the updated sequence converges, we move to the second stage. After estimating the sequencex^MLSE with an imprecise source position estimation, it is necessary to update the position estimation p^ by the following estimator
(34)p^=argminp∈ℝD×1,bf2p,x^=argminp∈ℝD×1,br¯H⋅ΠA¯sp,x^⊥⋅r¯,where
(35)ΠA¯sp,x^⊥=IMN−A¯sp,x^A¯sHp,x∧A¯sp,x∧−1A¯sHp,x^.Because the traditional grid search requires repeated computations, the estimator (34) can be minimized using Newton’s iteration method. The updated vector ωo=pTRebTImbTT is given by
(36)ωq+1o=ωqo−μqH−1ωqohωqo,where μq∈0,1 is the step size of the qth iteration. hωo and Hωo are the gradient vector and the Hessian matrix of (34), respectively.Applying the first-order derivation operator of the orthogonal projection matrixΠA¯sp,x^⊥, we can obtain the gradient vector hωo. Therefore, Hessian matrix Hωo which is essential to Newton’s-type method can be derived based on equation (38). As a result, the approximation has no effect on the numerical value of the estimation. Differentiating f2p,x^ with respect to ωo, the ith element of the gradient vector hωo is
(37)∂f2p,x^∂<ωo>i=Rer¯H⋅−ΠA¯sp,x^⊥⋅∂A¯sp,x^∂<ωo>i⋅A¯s†p,x^−ΠA¯sp,x∧⊥⋅∂A¯sp,x∧∂<ωo>i⋅A¯s†p,x∧H⋅r¯.Following some algebraic manipulations, thei,j element in the Hessian matrix Hωo is given by
(38)∂2f2p,x^∂<ωo>i∂<ωo>j≈2Rer¯H⋅A¯s†p,x∧H⋅∂A¯sHp,x^∂<ωo>i⋅ΠA¯sp,x^⊥⋅∂A¯sp,x^∂<ωo>j⋅A¯s†p,x^⋅r¯.With the matrix transformationvecXYZ=ZT⊗X⋅vecY, the gradient vector hωo can be expressed as
(39)hωo=−ReVA¯sp,x^⋅A¯s†p,x^r¯⊗r¯HΠA¯sp,x∧⊥T+ΠA¯sp,x∧⊥Hr¯⊗r¯HA¯s†p,x∧HT,where
(40)VA¯sp,x^=vec∂A¯sp,x∧∂pvec∂A¯sp,x∧∂Rebvec∂A¯sp,x∧∂ImbT.Based on the property thatΠA¯sp,x^⊥=ΠA¯sp,x^⊥⋅ΠA¯sp,x^⊥, the Hessian matrix Hωo can be formulated as follow
(41)Hωo≈2ReVA¯sp,x^⋅A¯s†p,x^r¯⊗ΠA¯sp,x^⊥A¯s†p,x∧r¯T⊗ΠA¯sp,x^⊥⋅VA¯sp,x^T,which is essential to the Newton’s-type method. Then, we can compute the gradient vector hωo and Hωo for Newton’s iteration algorithm in (36).The main formulas are given above. With appropriate step-size factors [37], the updated result can converge to the corresponding optimal result in several iterations, which is verified by the experiments in Section 4. The step-size factor μq can be adaptively determined to follow the strong Wolfe conditions [38], which means
(42)f2ωqo+μq⋅H−1ωqo⋅hωqo,x^≤f2ωqo,x^+μq⋅hTωqo⋅H−1ωqo⋅h1ωqo,h1Tωqo+μq⋅H−1ωqo⋅hωqo⋅H−1ωqo⋅hωqo≤c1hTωqo⋅H−1ωqo⋅hωqo,with 0<c1<1.Note that the data sequence is initially estimated using a relatively inaccurate position estimation. A good initial estimate is very important for subsequent iteration. In this paper, the position estimate for the first MLSE in (31) is set as the initial estimate, which is accurate enough for the proposed convergence. The procedure of the proposed DPD method is shown in Table 1, and for clarity, the flow diagram can be seen in Figure 2.Table 1
The direct localization algorithm for digital modulation signals.
Require: Obtain the initial guess of emitter positionp^0; the received signal of each sensor rm,m=1,2,⋯,M; thresholds ε and δ (sufficiently small and positive values, which are 10-6 in this paper).
A. First-stage processing (sequence estimation):
(i) Apply crossover and mutation toXelitei and generate P random discrete data sequences Xi=x1i,x1i,⋯,xP−pi,Xelitei as the ith population. The subset Xelitei contains the elite individuals in the last generation
(ii) Update the maximum of fitness value and the corresponding sequencex^=xbesti in Xi
(iii) Update the searching generationi⟵i+1
B. Second-stage processing (position estimation):
(i) Compute the gradient vectorhω and the hessian matrix Hω using (39) and (41)
(ii) Update the iterative resultωq+1 via (36), then q⟵q+1
(iii) Ifhωq≤ε, obtain p^=ωq1:D
(iv) Ifp^−p^0≤δ, output p^ as the estimated position, otherwise, reset p^0=p^ and return to the first-stage
Output: The estimated positionp^.Figure 2
Flow diagram of the proposed DPD method.


## 3.3. Computational Complexity Analysis

Table2 summarizes the numerical complexity, where Kstep1, Kstep2, and Kiter are the amount of looping in the procedure, respectively.Table 2
Computational complexity of each iteration.
Items infirst stage
The number of multiplications
Items insecond stage
The number of multiplications
f
fitness
x
,
p
^
M
2
N
2
+
1
M
N
ω
q
+
1
o
D
+
2
M
D
+
2
M
+
1
f
2
p
,
x
^
M
2
N
2
+
1
M
N
H
ω
o
4
M
3
N
2
+
2
M
2
N
2
h
ω
o
2
M
N
D
+
2
M
2
M
N
+
4
M
+
1
Inverse ofHωo
D
+
2
M
3
−
2
D
+
2
M
2
+
6
D
+
2
M
−
1
Total offirst stage
S
1
=
K
step
1
M
2
N
2
+
1
M
N
Total ofsecond stage
S
2
=
K
step
2
M
3
N
3
+
4
M
3
N
2
+
2
M
2
N
2
+
M
N
+
2
D
+
2
M
2
M
2
N
2
+
4
M
2
N
+
M
N
+
D
+
2
M
3
−
D
+
2
M
2
+
7
D
+
2
M
−
1
Total
K
iter
S
1
+
S
2


## 4. Simulation Results

In this section, several Monte Carlo simulations are introduced to verify the effectiveness of the proposed DPD methods. The performance of the proposed method is compared with that of the DPD methods with unknown and known waveforms, and the CRB given in Section4 is also compared. The simulation results exceeded 200 Monte Carlo tests on average. We focus on the root mean square error (RMSE) to evaluate the performance of our methods, which is defined by
(43)RMSE=1Nexp∑i=1Nexpp−p∧i22where Nexp is the number of the Monte Carlo trials and p∧i is the ith position estimate. Besides, two different localization scenarios are considered. The following simulations are carried out for these two cases. The source and sensor positions are indicated by random symbols.


### 4.1. Near-Field Source

In this subsection, the performance of the proposed DPD method is examined for both known and unknown transmitted waveforms. Figure3 shows the geometric position of the near-field scenario, where the sensors are located at (-3, 3) [km], (-3, -3) [km], (3, -3) [km], and (3, 3) [km] and the emitter source is located at (1.4, 1.88) [km]. Then, a QPSK signal is transmitted from the emitter, and 20 samples are taken of each symbol. The roll-off factor of the root-raised cosine filter is set as 0.35. The data sequence of the QPSK signal is random, the symbol rate is 10 [kHz], the length of the symbol sequence is 10, and the order K is equal to 401 in this case.Figure 3
Geometry for near-field source localization.Here, the proposed method and existing DPD methods for both unknown and known signals are compared. The signal-to-noise (SNR) varied from -10 to 35 [dB] at 5 [dB] intervals. Figure4 shows the RMSE of the proposed method and existing DPD methods based on the above conditions. It can be seen that for known waveforms, the performance of the proposed method can reach that of the DPD method. By making use of the modulation scheme, when the waveform is unknown, the proposed method outperforms the existing DPD method. Because the data sequence of the transmitted signals is unknown to the location system, the proposed method has a lower threshold than that of the DPD method for known waveforms. In addition, the estimation accuracy of the proposed DPD method can attain the corresponding CRB when the SNR is 0 [dB].Figure 4
RMSE of the localization methods and associated CRBs versus SNR.Further, we compare the distributions of the source position estimation when SNR =5 [dB]. In Figure5, the estimation results of these methods are compared together. In order to evaluate the accuracy of these methods intuitively, we introduce the uncertainty ellipse of circular error probability (CEP) [39] as an approximate measure to evaluate their performance. Figure 6 shows that 50% of the estimates of the two methods are encapsulated in an ellipse. The two ellipses of the proposed method and the DPD for known waveforms are about the same size, which is in good agreement with the RMSE results of the source position estimation.Figure 5
True position of near-field emitter and distributions of position estimates for proposed method and existing methods at SNR =5 [dB].Figure 6
Uncertainty ellipses and partial distributions of Figure5.


### 4.2. Far-Field Source

In this section, the performance of locating a far-field source is studied. The geometric structure of the position of the far-field scenario is seen in Figure7. The sensors are located at (0, 0) [km], (-43.76, 25.4) [km], (43.76, 25.4) [km], and (0, -42.5) [km] and the emitter source is located at (70, 94) [km].Figure 7
Geometry for far-field source localization.The other conditions are the same as those in the first set. As shown in Figure8, the proposed DPD method can reach the corresponding CRB. Due to the model errors, the estimation accuracy of the DPD method for unknown waveforms cannot achieve CRB even at high SNR. Figure 9 shows the distributions of these methods at SNR =5 [dB]. In Figure 10, the ellipse of the proposed method is slightly larger than that of the existing DPD method for known waveforms, but much smaller than that of unknown case. In general, our method performs significantly better and has a higher threshold compared with the DPD method for unknown waveforms. Besides, this method can also achieve the same accuracy as that for the DPD method for known waveforms under a sufficiently high SNR.Figure 8
RMSE of the localization methods and associated CRBs versus SNR.Figure 9
True position of far-field emitter and distributions of position estimates for proposed method and existing methods at SNR =5 [dB].Figure 10
Uncertainty ellipses and partial distributions of Figure9.


### 4.3. Different Symbol Rates

In this subsection, the performance of the proposed DPD method is evaluated for different symbol rates. Because the data sequence is actually random, the symbol rate is one of the most important attributes of a digital modulation signal as it determines the bandwidth of the received signal. In the following simulations, the position geometry is the same as that in Figure7 and QPSK is used as the modulation scheme for 10 symbols. The symbol rate varies from 5 [kHz] to 35 [kHz]. The roll-off factor of the root-raised cosine filter is 0.35, and 20 samples are taken of each symbol.Figure11 shows the performance of the proposed method and the existing DPD methods for unknown and known waveforms, respectively, at SNR =10 [dB]. The performance of the proposed method is close to that of the DPD method for known waveforms, which is consistent with the previous two experiments. As shown in Figure 11, the estimation accuracy of the proposed method and the existing DPD methods decreases with the increasing symbol rate.Figure 11
RMSE of the localization methods and associated CRBs versus symbol rate.


### 4.4. Different Sequence Length

We examined the performance of the proposed DPD method for another digital modulation scheme and various data sequence lengths. BPSK modulation scheme is used, and the location geometry of the emitting source and the sensors is the same as that in Figure7. The roll-off coefficient of the root-raised cosine filter is 0.35, the symbol rates are both 10 kHz, and 20 samples are taken of each symbol. In addition, the symbol sequence is also random in this case.Figure12 shows the RMSE of the source position estimation and the CRBs for various number of symbols. The results demonstrate that the proposed method can reach the relevant CRBs with high SNR. For a given number of symbols in a data sequence, the proposed method has a comparable performance. In addition, the CRBs and the RMSEs for the proposed DPD method show that the accuracy of the estimation results increases with increasing number of symbols in the transmitted data sequence.Figure 12
RMSEs and associated CRBs versus SNR for BPSK.


### 4.5. Running-Time Comparison

In this section, the running time of the proposed DPD method is compared with the previous DPD methods. All the simulations are conducted using MATLAB R2017a, which is equipped with an Intel Core i5-8400 CPU operating @ 2.8GHz and 16 GB of random access memory (RAM). The running-time results in Table3 are averaged of the total number of Monte Carlo experiments in each case.Table 3
Comparison of run time (SNR =5 [dB]).
Localization algorithm
Running time in second
Sequence estimation(20 samples/symbol)
Position estimation
Total
Proposed
36.99
2.37
39.36
DPD with unknown waveform (grid search)
0
145.08
145.08
DPD with known waveform (grid search)
0
148.68
148.68As a result, the proposed methods can converge to the optimal result with a good initial guess in several iterations. Compared with the previous DPD methods based on the grid search, the proposed method has a significant superiority in the position estimation stage. The following conditions for the simulation are the same as those in Section4. It can be seen from the statistical data that sequence estimations take the greatest computational power in the procedure of the proposed algorithm.


## 4.1. Near-Field Source

In this subsection, the performance of the proposed DPD method is examined for both known and unknown transmitted waveforms. Figure3 shows the geometric position of the near-field scenario, where the sensors are located at (-3, 3) [km], (-3, -3) [km], (3, -3) [km], and (3, 3) [km] and the emitter source is located at (1.4, 1.88) [km]. Then, a QPSK signal is transmitted from the emitter, and 20 samples are taken of each symbol. The roll-off factor of the root-raised cosine filter is set as 0.35. The data sequence of the QPSK signal is random, the symbol rate is 10 [kHz], the length of the symbol sequence is 10, and the order K is equal to 401 in this case.Figure 3
Geometry for near-field source localization.Here, the proposed method and existing DPD methods for both unknown and known signals are compared. The signal-to-noise (SNR) varied from -10 to 35 [dB] at 5 [dB] intervals. Figure4 shows the RMSE of the proposed method and existing DPD methods based on the above conditions. It can be seen that for known waveforms, the performance of the proposed method can reach that of the DPD method. By making use of the modulation scheme, when the waveform is unknown, the proposed method outperforms the existing DPD method. Because the data sequence of the transmitted signals is unknown to the location system, the proposed method has a lower threshold than that of the DPD method for known waveforms. In addition, the estimation accuracy of the proposed DPD method can attain the corresponding CRB when the SNR is 0 [dB].Figure 4
RMSE of the localization methods and associated CRBs versus SNR.Further, we compare the distributions of the source position estimation when SNR =5 [dB]. In Figure5, the estimation results of these methods are compared together. In order to evaluate the accuracy of these methods intuitively, we introduce the uncertainty ellipse of circular error probability (CEP) [39] as an approximate measure to evaluate their performance. Figure 6 shows that 50% of the estimates of the two methods are encapsulated in an ellipse. The two ellipses of the proposed method and the DPD for known waveforms are about the same size, which is in good agreement with the RMSE results of the source position estimation.Figure 5
True position of near-field emitter and distributions of position estimates for proposed method and existing methods at SNR =5 [dB].Figure 6
Uncertainty ellipses and partial distributions of Figure5.


## 4.2. Far-Field Source

In this section, the performance of locating a far-field source is studied. The geometric structure of the position of the far-field scenario is seen in Figure7. The sensors are located at (0, 0) [km], (-43.76, 25.4) [km], (43.76, 25.4) [km], and (0, -42.5) [km] and the emitter source is located at (70, 94) [km].Figure 7
Geometry for far-field source localization.The other conditions are the same as those in the first set. As shown in Figure8, the proposed DPD method can reach the corresponding CRB. Due to the model errors, the estimation accuracy of the DPD method for unknown waveforms cannot achieve CRB even at high SNR. Figure 9 shows the distributions of these methods at SNR =5 [dB]. In Figure 10, the ellipse of the proposed method is slightly larger than that of the existing DPD method for known waveforms, but much smaller than that of unknown case. In general, our method performs significantly better and has a higher threshold compared with the DPD method for unknown waveforms. Besides, this method can also achieve the same accuracy as that for the DPD method for known waveforms under a sufficiently high SNR.Figure 8
RMSE of the localization methods and associated CRBs versus SNR.Figure 9
True position of far-field emitter and distributions of position estimates for proposed method and existing methods at SNR =5 [dB].Figure 10
Uncertainty ellipses and partial distributions of Figure9.


## 4.3. Different Symbol Rates

In this subsection, the performance of the proposed DPD method is evaluated for different symbol rates. Because the data sequence is actually random, the symbol rate is one of the most important attributes of a digital modulation signal as it determines the bandwidth of the received signal. In the following simulations, the position geometry is the same as that in Figure7 and QPSK is used as the modulation scheme for 10 symbols. The symbol rate varies from 5 [kHz] to 35 [kHz]. The roll-off factor of the root-raised cosine filter is 0.35, and 20 samples are taken of each symbol.Figure11 shows the performance of the proposed method and the existing DPD methods for unknown and known waveforms, respectively, at SNR =10 [dB]. The performance of the proposed method is close to that of the DPD method for known waveforms, which is consistent with the previous two experiments. As shown in Figure 11, the estimation accuracy of the proposed method and the existing DPD methods decreases with the increasing symbol rate.Figure 11
RMSE of the localization methods and associated CRBs versus symbol rate.


## 4.4. Different Sequence Length

We examined the performance of the proposed DPD method for another digital modulation scheme and various data sequence lengths. BPSK modulation scheme is used, and the location geometry of the emitting source and the sensors is the same as that in Figure7. The roll-off coefficient of the root-raised cosine filter is 0.35, the symbol rates are both 10 kHz, and 20 samples are taken of each symbol. In addition, the symbol sequence is also random in this case.Figure12 shows the RMSE of the source position estimation and the CRBs for various number of symbols. The results demonstrate that the proposed method can reach the relevant CRBs with high SNR. For a given number of symbols in a data sequence, the proposed method has a comparable performance. In addition, the CRBs and the RMSEs for the proposed DPD method show that the accuracy of the estimation results increases with increasing number of symbols in the transmitted data sequence.Figure 12
RMSEs and associated CRBs versus SNR for BPSK.


## 4.5. Running-Time Comparison

In this section, the running time of the proposed DPD method is compared with the previous DPD methods. All the simulations are conducted using MATLAB R2017a, which is equipped with an Intel Core i5-8400 CPU operating @ 2.8GHz and 16 GB of random access memory (RAM). The running-time results in Table3 are averaged of the total number of Monte Carlo experiments in each case.Table 3
Comparison of run time (SNR =5 [dB]).
Localization algorithm
Running time in second
Sequence estimation(20 samples/symbol)
Position estimation
Total
Proposed
36.99
2.37
39.36
DPD with unknown waveform (grid search)
0
145.08
145.08
DPD with known waveform (grid search)
0
148.68
148.68As a result, the proposed methods can converge to the optimal result with a good initial guess in several iterations. Compared with the previous DPD methods based on the grid search, the proposed method has a significant superiority in the position estimation stage. The following conditions for the simulation are the same as those in Section4. It can be seen from the statistical data that sequence estimations take the greatest computational power in the procedure of the proposed algorithm.


## 5. Conclusion

This paper proposed a DPD method for positioning the digital modulation signal based on TDOA measurements. By utilizing the digital modulation scheme of the transmitted signals, we devise an ML-based cost function for the given digital modulation signals. Considering that the estimator is a mixed-integer optimization problem, a modified genetic algorithm and Newton method are used to estimate the data sequence and source location alternately. Numerical simulations have proven that the proposed algorithm can converge to the optimal solution in several iterations. In addition, the CRB for unknown waveform and unknown waveform are derived, and the performance of the proposed method is compared. The results of existing DPD methods without modulation schemes are investigated to show the superior performance of the proposed algorithm. Numerical results and the distribution of estimates reveal the benefits of the modulation information and show that the proposed algorithm can approach that of the existing method with known waveforms.Currently, the proposed algorithm only uses the TDOA information to locate single emitter. In the future work, we will extend the proposed algorithm to hybrid TDOA/FDOA localization in multiple-target localization scenarios.



---

*Source: 5104172-2022-09-27.xml*