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

This paper proposes a direct position determination (DPD) method for a digital modulation signal based on time di ﬀ erence 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 veri ﬁ ed using comprehensive simulations.


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][2][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 concen-trate the raw data from each sensor in a data processing center and then estimate the source location in one step.
DPD methods [17][18][19][20][21][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 mixedinteger 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 case The rest of this paper is organized as follows. In Section 2, 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, O N×n is a N × n matrix with all-zero entries and I N is a N dimensional identity matrix; A T , A H , A -1 , and A † are the transpose, the conjugate transpose, the inverse, and the Moore-Penrose inverse of matrix A; kak denotes the Euclidean norm of vector a; <·> i,j is the element at i th row and j th column of its argument and <·> i denotes the i th element; vecf·g, Ef·g, and b·c are the vectorization, 2 Wireless Communications and Mobile Computing expectation, and rounded down operation; ⊗ is the Kronecker product.

Signal Model.
We consider a two-dimensional (D = 2) scenario shown in Figure 1 including an unknown emitter at p = ½x y T ∈ ℝ D×1 and a location system equipped with M known sensors, which are located at u m = ½x m y m T ∈ ℝ D×1 . Suppose the emitter transmits a digital signal with a known modulation scheme, but the symbol sequence is unknown to the location system. 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 the m th sensor is given by for m = 1, 2, ⋯, M, where τ m = kp − u m k/c represents the signal propagation time between the emitter and the m th sensor, and t 0 is transmitted time. c is the signal propagation speed. b m is the m th path loss coefficient and w m ðtÞ is a complex additive white Gaussian noise with zero mean and the power density σ 2 n . By denoting with T the symbol period and the baseband chip waveform hðtÞ, the baseband complex envelope of the signal sðtÞ can be written as where x l ∈ f0, 1g represents the l th symbol of the unknown transmitted sequence. In this paper, the chip waveform hðtÞ 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 = ½x 1 x 2 ⋯ x L T can be regarded as a functional mapping relation sðxÞ. 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 the m th sensor and the first sensor is given by the following formula for m = 1, 2, ⋯, M. Then, the received signal can be rewritten as Using the properties of the Fourier transform that and defining the function sinc ðνÞ = sin ðπνÞ/ðπνÞ, the delayed signal in (4) is given bỹ After being sampled at t = nT s in each interval, where T s denotes the sampling period, the received signal of sensor m can be rewritten as and then the delayed signal waveform in (6) is transformed intõ In practice, the coefficient K cannot be infinite, so we refer K to a limited non-infinite integer K which is larger than bτ m,1 /T s c, so we havẽ To minimize the truncation errors, using a higher order K will lead to higher modeling accuracy. We further composite all the N observation samples into a vector and yield in which

Wireless Communications and Mobile Computing
With the definition SincðkÞ = sinc ðk − τ m,1 /T s Þ, we have the matrix for m = 2, ⋯, M and A 1 = I N . Multiplied with the special matrix A m , the transmitted waveforms is delayed with τ m,1 . Therefore, the waveformss m ðtÞ is produced with the s 1 ðtÞ and matrix A m with corresponding time delay. Different from the conventional methods in the frequency domain, A m 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 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 sequence x 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 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 vector x whose CRB is CRBðxÞ, if we define a new real vector y = Jx, where J is an invertible matrix, the CRB of vector y can be obtained by CRBðyÞ = J ⋅ CRBðxÞ ⋅ J T [35].  Let b φ 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 b φ satisfies the following information inequality where J φ is the FIM of φ defined as The log-likelihood ratio of the corresponding signal model (10) is given by 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 where 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 Then, the CRB for the unknown waveform is given by In the case of known waveforms, the emitter position p and path loss coefficient b are retained in ω o = ½p T Re ðb T Þ Im ðb T Þ T . Based on the log-likelihood function (17), the FIM of unknown parameter ω o can be partitioned into sub-matrices, which reduces to Then, we can compute the p corner of the CRB matrix as follow The diagonals of CRB

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 that w is white complex Gaussian noise with zero mean and variance σ 2 n , the parameters ρ = ½x T p T b T T need to be estimated. To decouple the uninterested parameters b, the vector r of observation samples can be rewritten as where : After omitting the constant terms, we can minimize the log-likelihood function (17) Submitting (26) into the objective function (17), we have where Note that the optimization model is a mixed-integer nonlinear programming problem, which can be formulated as follows It can be seen that ML sequence estimation (MLSE) has been used to obtain the data sequence of the transmitted signal. If the estimated parameter p 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 ρ = ½x T p T b T T into two groups, fxg and fp, bg. 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 estimatep 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 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 contains 2 L 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 aŝ 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 f0, 1g 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 = ½x 1 , x 2 , ⋯, x L T and y = ½y 1 , y 2 , ⋯, y L T , which are called parents here. The following Laplace distribution is used 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 = ½w 1 , w 2 , ⋯, w L T and z = ½z 1 , z 2 , ⋯, z L T can be obtained as 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 element x l is equal to either 0 or 1 with probability 0.5 if x l 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.

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 estimationp by the following estimator Because the traditional grid search requires repeated computations, the estimator (34) can be minimized using Newton's iteration method. The updated vector ω o = ½p T Re ðb T Þ Im ðb T Þ T is given by where μ q ∈ ð0, 1Þ is the step size of the q th 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 s ðp,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 f 2 ðp,xÞ with respect to ω o , the i th element of the gradient vector hðω o Þ is Following some algebraic manipulations, the ði, jÞ element in the Hessian matrix Hðω o Þ is given by With the matrix transformation vecðXYZÞ = ðZ T ⊗ XÞ ⋅ vecðYÞ, the gradient vector hðω o Þ can be expressed as 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 with 0 < c 1 < 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 7 Wireless Communications and Mobile Computing DPD method is shown in Table 1, and for clarity, the flow diagram can be seen in Figure 2. Table 2 summarizes the numerical complexity, where K step1 , K step2 , and K iter are the amount of looping in the procedure, respectively.

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 Section 4 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 where N exp is the number of the Monte Carlo trials andp ðiÞ is the i th 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.

Near-Field Source.
In this subsection, the performance of the proposed DPD method is examined for both known and unknown transmitted waveforms. Figure 3 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.
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. Figure 4 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].
Further, we compare the distributions of the source position estimation when SNR =5 [dB]. In Figure 5, 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.

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 Figure 7. The other conditions are the same as those in the first set. As shown in Figure 8, 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. Require: Obtain the initial guess of emitter positionp 0 ; the received signal of each sensor r m , 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 to X i elite and generate P random discrete data sequences 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 Figure 7 and QPSK is used as the modulation scheme for 10 symbols. The symbol rate varies from 5 [kHz] to 35 [kHz]. The rolloff factor of the root-raised cosine filter is 0.35, and 20 samples are taken of each symbol. Figure 11 shows the performance of the proposed method and the existing DPD methods for unknown and known waveforms, respectively, at SNR =10 [dB]. The

Items in first stage The number of multiplications Items in second stage
The number of multiplications      13 Wireless Communications and Mobile Computing 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.

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 Figure 7. 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. Figure 12 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.

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   Table 3 are averaged of the total number of Monte Carlo experiments in each case. As 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 Section 4. It can be seen from the statistical data that sequence estimations take the greatest computational power in the procedure of the proposed algorithm.

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.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.