**Adaptive Multiclass Mahalanobis Taguchi System for Bearing Fault Diagnosis under Variable Conditions**

**Ning Wang 1,2,3, Zhipeng Wang 1,2,3,\* , Limin Jia 1,2,3,\*, Yong Qin 1,2,3, Xinan Chen 1,2,3 and Yakun Zuo 1,2,3**


Received: 16 November 2018; Accepted: 19 December 2018; Published: 21 December 2018

**Abstract:** Bearings are vital components in industrial machines. Diagnosing the fault of rolling element bearings and ensuring normal operation is essential. However, the faults of rolling element bearings under variable conditions and the adaptive feature selection has rarely been discussed until now. Thus, it is essential to develop a practicable method to put forward the disposal of the fault under variable conditions. Considering these issues, this paper uses the method based on the Mahalanobis Taguchi System (MTS), and overcomes two shortcomings of MTS: (1) MTS is an effective tool to classify faults and has strong robustness to operating conditions, but it can only handle binary classification problems, and this paper constructs the multiclass measurement scale to deal with multi-classification problems. (2) MTS can determine important features, but uses the hard threshold to select the features, and this paper selects the proper feature sequence instead of the threshold to overcome the lesser adaptivity of the threshold configuration for signal-to-noise gain. Hence, this method proposes a novel method named adaptive Multiclass Mahalanobis Taguchi system (aMMTS), in conjunction with variational mode decomposition (VMD) and singular value decomposition (SVD), and is employed to diagnose the faults under the variable conditions. Finally, this method is verified by using the signal data collected from Case Western Reserve University Bearing Data Center. The result shows that it is accurate for bearings fault diagnosis under variable conditions.

**Keywords:** fault diagnosis; bearing; SVD; VMD; adaptive Multiclass Mahalanobis Taguchi System

### **1. Introduction**

Rolling element bearings have wide applications in industrial machines and are one of the most critical components. If faults occur in bearings, equipment could be damaged and disasters might happen consequently. Therefore, it is essential to monitor the health conditions of bearings. The analysis of vibration signals has been a hot research pot and used to detect faults of bearings. It is crucial to recognize faults occurring in bearings and avoid fatal breakdowns as early as possible. For decades, many researchers have conducted extensive research on fault diagnosis. At present, the fault diagnosis methods are divided into model-based methods and data-driven methods. The model-based methods generally build on the physics of the process, generating the residuals between the measure process variables and estimates [1], such as Hidden Markov Modeling (HMM) which is successfully applied to bearing fault detection and diagnosis [2]. Autoregressive modelling [3]

also has had excellent performance in bearing fault diagnosis. The accelerated degradation testing (ADT) [4] method is useful in fault diagnosis and lifespan prediction, and reference [5] presents a new approach using observer-based residual generation with no complicated design constraints to establish the relationship between the state estimation error and the fault signal. For the data-driven method, which is based on the historical data and does not need accurate mathematic and priori knowledge, it has wide applications in fault diagnosis. For example, Bayesian network is an excellent data-driven diagnosis method [6–8]. Machine learning algorithms, such as K-nearest Neighbor (KNN), Deep Convolutional Neural Networks (DCNN), and auto-encoders are effective in the fault diagnosis of various bearings [9–12]. There is also decision tree [13], Support Vector Machine (SVM) [14], and wavelet transform [15]. Artificial neural network (ANN) and Trace Ratio Linear Discriminant Analysis are other methods of diagnosing the bearings fault [16,17]. In addition, the vibration signal under the various operating conditions (especially in low rotational speed) is non-stationary and non-linear, the characteristic defect frequencies move continuously with the change of rotating speed, and is the same as the bearing that is going to break down. If the bearing progresses toward failure, the nonlinear features also start to be dominated by stochastic signal. With the vibration signal measured for diagnosing the bearing fault under variable condition, it will be difficult to diagnose the fault of bearing by using the traditional methods. Reference [18] proposed a signal selection scheme based upon two order tracking techniques from complicated non-stationary operational measured vibrations. Reference [19] reviewed features extraction methods and its application on bearing vibration signal and presents an empirical study of feature extraction methods in low rotational speed. Reference [20] proposes the estimation of instantaneous speed relative fluctuation (ISRF) in a vibration speed. Reference [21] proposes Stacked Convolutional Autoencoders (SCAE) together with DCNN in stationary and non-stationary speed operation. Graph-based rebalance semi-supervised learning (GRSSL) [22], weighted self-adaptive evolutionary extreme learning machine (WSaE-ELM) [23] and Singular Spectrum Analysis [24] are effective in diagnosing the fault under variable conditions.

However, although the aforementioned methods are effective for bearing faults diagnosis, the feature selection part is often non-adaptive or unexplainable. The conventional methods mainly extract and select features manually, which relies heavily on the experts' knowledge and experience. Since the signals acquired in the real world might be various in many different aspects, the features selected manually might be sensitive in the variations of operation conditions and import inevitable errors for fault diagnosis. The deep learning algorithms can extract features automatically and overcome this drawback. However, deep learnings are data-hungry and require plenty of training data which are hardly acquired in practice, especially the faulty data under different conditions. Besides, the features acquired by deep learnings are unexplainable. Therefore, the aforementioned methods can hardly diagnose the faults under various operation conditions in practice. In contrast, the Mahalanobis Taguchi System (MTS) is more robust than the other methods in various operating conditions [25].

MTS offers a tool to determine important features and optimize the system. It is a different form compared with the other classification methods, because this classification model of measurement scale is constructed by using the class samples. It is useful to diagnose bearing faults under various conditions because the different pattern could be identified by using the Mahalanobis distance (MD) and Taguchi method. In this paper, MD is used to calculate the distance of the correlations between the benchmark and others, and the distance could be measured without the volatility of data. The advantage of MD is that it takes into consideration the correlations between the features and this consideration is very important in pattern analysis, which is why MTS is suitable for bearing fault diagnosis under various conditions [26]. On the other hand, the Taguchi method is used to select features without manual intervention, which could improve the robustness of the algorithm. MTS also offers an effective tool for multivariate analysis [27], considering that the bearing faults can be classified according to locations, such as inner race, outer race and rolling element [28]. However, when the conventional MTS is used for bearing fault diagnosis, misclassifications might occur due to less adaptivity of the threshold configuration for signal-to-noise gain. During the feature selection, the

threshold is normally set as a constant, which might result in the overfitting problem. If the threshold value is too large, some critical features might be eliminated. On the contrary, if the threshold value is too low, some useless or harmful feature might be selected. Therefore, if the threshold value does not march the training data sufficiently, misclassifications emerge.

To overcome this drawback, this paper presents a novel method named adaptive Multiclass Mahalanobis Taguchi system (aMMTS) for bearing fault diagnosis. This method employs the MTS for multi-classification by considering different conditions as different benchmarks respectively. The results are based on the minimum MDs between the data and each benchmark data, and the label of the data is determined to be consistent with the label of the benchmark data whose MD is minimum. The method can be described briefly as follows: Firstly, after the Mahalanobis space (MS) is constructed by the two-level orthogonal array of the Taguchi method, aMMTS calculates the MDs from the data to the benchmark data and obtains features' signal-to-noise ratios (SNRs) and gain values by using two-level orthogonal array. Secondly, the features are selected adaptively by recalculating the MDs via rearranging the order of features' gain values by ascending and descending. Therefore, the proposed method is able to select the best classification result according to the adaptive chosen sequence of features' SNRs instead of a hard threshold. Here, the sequence of SNR is determined by a function, which selects several maximum or minimum features to calculated MDs. Finally, a set of features with the best results is selected as the final feature vector. In this method, two different sets of training samples are employed to calculate the SNRs respectively and obtain the final feature vectors respectively. By the aforementioned improvement, the proposed aMMTS is capable to overcome the drawback of the conventional MTS and prevent the over-fitting problem. Therefore, the aMMTS is insensitive to the operation conditions and can be employed for bearing fault diagnosis.

Moreover, this method is combined with variational mode decomposition (VMD) [29] and singular value decomposition (SVD) to diagnose the faults. VMD is an entirely non-recursive algorithm, and is used to decompose the signal. It has been proven that due to the characteristics of nonlinear vibration in the bearings, VMD is more efficient than empirical mode decomposition (EMD) and Fourier transform (FT) under variable condition. SVD is used to extract the features. Therefore, VMD and SVD are employed in this paper.

The rest of this article is organized as follows: Section 2 introduces the algorithms involved in this paper. Section 3 illustrates the experiments to validate the proposed method. Section 4 is the conclusions.

### **2. Methodology**

In this paper, the main steps of fault diagnosis are signal decomposition, feature extraction and fault detection. The detailed process and method are as follows:

Step1: The VMD in conjunction with wavelet denoising is employed to eliminate the noises and decompose the raw signals;

Step2: Extracting features from the decomposed signals by SVD;

Step3: The proposed aMMTS is employed for the fault diagnosis.

The steps of this method are shown in Figure 1.

### *2.1. VMD*

After the wavelet denoising is used to estimate the noise of the raw signal, this paper employs VMD to decompose non-stationary signals. VMD can decompose a signal into different simple intrinsic mode functions, whose frequency center and bandwidth are band-limited and determined by iterative searching for the optimal solution of the variational model. The constrained formula is given as:

$$\min\_{\mu\_k \downarrow \omega\_k} \{ \sum\_k ||\partial\_t[(\delta(t) + \frac{j}{\pi t}) \* \mu\_k()t]| e^{-j\omega\_k t} ||\_2^2 \} \qquad \text{s.t.} \sum \mu\_k = f \tag{1}$$

*Sensors* **2019**, *19*, 26

where *μ<sup>k</sup>* is the sub-signals, *ω<sup>k</sup>* represents the center frequency of sub-modes. The optimal solution can be solved as the minimization problem, which could be addressed by introducing a quadratic penalty and Lagrangian multipliers [30]:

$$L(\{\mu\_k\}\{\omega\_k\}, \lambda) = a \sum\_k |\langle \partial\_l(\delta(t) + \frac{j}{\pi t}) \* \mu\_k(t) \rangle| e^{-j\omega\_k t} |\langle \frac{\pi}{2} + ||f(t) - \sum\_k \mu\_k(t)||\_2^2 + \langle \lambda(t), f(t) - \sum\_k \mu\_k(t) \rangle \tag{2}$$

*α* denotes the balancing parameter of the constraint.

**Figure 1.** The scheme of the proposed fault diagnosis method.

All sub-signals *μ<sup>k</sup>* are updated for all *ω* ≥ 0 as follow:

$$\mu\_k^{n+1} \gets \frac{\hat{f} - \sum\_{ik} \mu\_i^n + \frac{\hat{\lambda}^n}{2}}{1 + 2a(\omega - \omega\_k^n)^2} \tag{3}$$

*μ*ˆ1 *<sup>k</sup>*, *<sup>ω</sup>*<sup>ˆ</sup> <sup>1</sup> *<sup>k</sup>* and *<sup>λ</sup>*<sup>ˆ</sup> <sup>1</sup> are initialized to all zeroes.

All center frequency of sub-modes *ω<sup>k</sup>* are updated as follow:

$$
\omega\_k \leftarrow \frac{\int\_0^\infty \omega \left| \boldsymbol{\mu}\_k^{n+1}(\omega) \right|^2 d\omega}{\int\_0^\infty \left| \boldsymbol{\mu}\_k^{n+1}(\omega) \right|^2 d\omega} \tag{4}
$$

*Sensors* **2019**, *19*, 26

End for

$$
\hat{\lambda}^{n+1} \gets \hat{\lambda}^n + \tau (\hat{f} - \sum\_{k} \hat{\mu}\_k^{n+1}) \tag{5}
$$

Until:

$$\frac{\sum\_{k} \left\| \left\| \hat{\mu}\_{k}^{n+1} - \hat{\mu}\_{k}^{n} \right\| \right\|\_{2}^{2}}{\left\| \hat{\mu}\_{k}^{n} \right\|\_{2}^{2}} < \varepsilon \tag{6}$$

### *2.2. SVD*

After the signals are decomposed into several modes by VMD, the features are extracted from the modes by SVD, and can be constructed as a feature matrix. SVD is a powerful tool for feature extraction in linear algebra. According to SVD, the matrix could be decomposed as follow:

$$X = \mathsf{U}\mathsf{l}\omega\mathsf{V}^T\tag{7}$$

where *X* represents a *m* × *n* matrix. There are two orthogonal matrices: matrix *U*(*m* × *m*) and *V*(*n* × *n*), and a singular diagonal matrix *ω*(*ωij* = 0, *i* = *j* and *ω*<sup>11</sup> ≥ *ω*<sup>22</sup> ≥ ··· ≥ 0), the diagonal element *ω*11, *ω*22, ··· , *ωmm* is the singular value of *X*. *U* is called left singular vector, and the columns of *U*. *V* is called right singular vector. The columns of *U* and *V* especially are orthogonal to each other, and are base vector. To obtain more intrinsic information in the matrix, the singular vectors are selected. As a consequence, SVD is employed to decompose the eigenmatrix, and obtained the singular value vectors (*ω*11, *ω*22, ··· , *ωmm*).

### *2.3. Mahalanobis–Taguchi System*

Mahalanobis–Taguchi System is a pattern recognition method integrated by the MD, orthogonal table and other tools such as SNR that are proposed by Taguchi [31], who introduces the experimental design of the field SNR to pattern recognition, which can reduce the dimensions of data, and use the orthogonal table to construct the MS. The MSs are used to calculate the MDs of the experimental data, and the valid features are distinguished by the SNR. Then, the MSs are recalculated by using the valid features. Finally, the results are obtained. The calculation of the MD is described as follows:

### 2.3.1. Mahalanobis Distance

The MD is a method of using normal data to normalize the fault data to compute the average distance between points and groups using normal data, the calculation formula of MD is as follows:

$$\mathbf{^jMD\_j} = \frac{1}{k} \mathbf{Z\_{i\bar{j}}^T \mathbf{C}^{-1} Z\_{i\bar{j}}} \mathbf{C}^{-1} \mathbf{Z\_{i\bar{j}}} \tag{8}$$

$$Z\_{i\bar{j}} = \left\{ z\_{1\circ}, z\_{2\circ}, z\_{3\circ}, \dots, z\_{k\circ} \right\}, z\_{i\bar{j}} = \frac{x\_{i\bar{j}} - \overline{x}\_{\bar{i}}}{s\_{\bar{i}}} \tag{9}$$

*MDj* represents the MD of the *j*th sample, *k* represents the number of the feature, *xij* represents the *i*th feature's value of the of the *j*th sample *xi* represents the mean of *i*th feature, *si* represents the Standard deviation of *i*th feature, *C*−<sup>1</sup> represents the inverse matrix of the correlation coefficient matrix.

### 2.3.2. Taguchi Method

In the Mahalanobis–Taguchi System, the MD measures the deviation of the test value from the normal value. The Taguchi method is able to select the features which have a larger contribution to identifying bearing faults, and then use the selected valid features to calculate the MD; the method is as below:

### Orthogonal Array

Selecting the appropriate two-level orthogonal array, and then the *k*-original features obtained by VMD and SVD, are arranged into each column of the orthogonal array. In the orthogonal array, "1" indicates that the feature is selected, "2" indicates that the feature is not selected, and a MS is generated according to each row of the orthogonal array.

### SNR and Its Gain

The main function of signal noise ratio is to select a valid feature, the calculation formula of generating SNR *η<sup>i</sup>* in the *i* line based on the orthogonal array.

$$\eta\_i = -10 \lg \frac{1}{N} \sum\_{j=1}^{N} MD\_{ij} \tag{10}$$

*j* ∈ [1, *N*] represents the number of training samples.

*η<sup>i</sup>* represents the recognition effects of the characteristic feature, the valid feature is selected by comparing the mean of SNR of each characteristic feature at two levels. The formula is as follow:

$$\overline{\eta\_{\dot{j}}} = \frac{\sum \eta\_{\dot{i}}}{m} \dot{j} = 1,2 \tag{11}$$

$$\overline{\Delta \eta\_{\bar{j}}} = \sum \left( \frac{\overline{\eta\_{\bar{j}}}}{m} \right) \tag{12}$$

*j* = 1, 2 represents two levels, *i* represents the number of rows in the MS, *η*<sup>1</sup> represents that in the level of '1', the recognition effects to identify abnormal conditions of using this feature. *η*<sup>2</sup> represents that in the level of '2', the recognition effects to identify abnormal conditions of not using this feature. If Δ*η<sup>j</sup>* represents the SNR gain, R = {Δ*η<sup>j</sup>* ' ' ' *η<sup>i</sup>* > 0 that indicates that this feature is a valid feature, if Δ*η<sup>j</sup>* < 0 that indicates that this feature is not a valid feature.

### 2.3.3. Adaptive Multiclass MTS

After SNR gain is calculated, to overcome the drawback of the conventional MTS during the threshold selection, this paper presents the adaptive multiclass MTS. There are several following improvements:

(1) Solving the multiple-classification problem. Selecting samples from each kind as the benchmark data. Then, the distances between other data and each benchmark data are calculated by MTS. Therefore, the label of the benchmark data with the minimum distance is selected as the label of the training data.

(2) Selecting the features adaptively. The feature sequence is selected several times, and the best fault recognition is the one which is minus MDs between it and benchmark. It solves the error problem caused by hard threshold selection.

(3) Avoiding the overfitting problem. Since the SNR gains are calculated by training samples, different training samples are employed to calculate the MDs and identify bearing faults, and the difference validation samples are set to validate the identified result and recalculate MDs.

This method is shown in Figure 2.

This adaptive multiclass MTS can be described as follow:

First, the data are labeled, the multi-class MSs are constructed, MDs between MSs are calculated and SNR gains are obtained. This step is shown as Figure 3; *m* is the number of samples, *n* is the number of features, and *j* is the number of the label. The training data are divided into three parts, and one of them is named as Benchmark. The data *Aj* and *Cj* represents one kind of data respectively (such as normal, fault of inner race, outer race, rolling element), and data B includes all kinds of data.

**Figure 2.** The step of Multiclass–Mahalanobis–Taguchi system (aMMTS).

Second, new sequences of feature parameters are generated by the sequence of features' SNR gains in ascending and descending order, then two collections are obtained from the above sequences based on the ascending or descending order, with the ascending and descending collection as follows:

$$\sigma\_k = (\overline{\triangle \eta\_{ki}}, \dots, \overline{\triangle \eta\_{kj}} \Big| \overline{\triangle \eta\_{ki}} > \overline{\triangle \eta\_{kj}}, i < j, i \in [1, N]) \tag{13}$$

$$q\_k = \left(\overline{\triangle \eta\_{ki}}, \dots, \overline{\triangle \eta\_{kj}}\Big|\overline{\triangle \eta\_{ki}} < \overline{\triangle \eta\_{kj}}, i > j, i \in [1, N]\right) \tag{14}$$

$$\mathcal{R} = \{r\_k | k \in [1, N]\}\tag{15}$$

$$\mathbf{Q} = \{q\_k | k \in [1, N]\}\tag{16}$$

This step is shown in Figure 4:

 *n* η η η § · ¨ ¸ ¨ ¸ ¨ Δ Δ Δ ¸ ¨ ¸ ¨ ¸ © ¹ U*k k kj ki kj* \_ > @ + ++ + η ηη η L <sup>&</sup>gt; *i ji N* < ∈ *qk ki kj ki kj* \_ > @ + ++ + η ηη η <sup>&</sup>lt; *i ji N* > ∈ *eg* Δ >Δ >Δ > >Δ >Δ >Δ η η η ηη η *n i* ( ) ( ) ( ) ( ) *m m R* η η η ηηη ηηη η ½ <sup>Δ</sup> ° ° Δ Δ =Δ Δ Δ ® ¾ ° ° ΔΔΔ Δ ¯ ¿ ( ) ( ) ( ) ( ) *i i Q* η η η ηηη ηηη η ½ <sup>Δ</sup> ° ° Δ Δ =Δ Δ Δ ® <sup>¾</sup> ° ° ΔΔΔ Δ ¯ ¿

Third, and the positions of feature's SNR gain are the same as the positions of corresponding features in the sequences, the features that are used to recalculate the MDs are selected by the corresponding feature's SNR gain. This step is shown in Figure 5. The MDs that are between each kind of *A* and *Ci* are calculated.

( ) ( ) ( ) ( ) *m m R* η η η ηηη ηηη η ½ <sup>Δ</sup> ° ° Δ Δ =Δ Δ Δ ® ¾ ° ° ΔΔΔ Δ ¯ ¿ *n n m m mn CC C CC C CC C* § · ¨ ¸ ¨ ¸ ¨ ¸ ¨ ¸ © ¹ ( ) Δ Δ η η *CC C <sup>n</sup> A C C m m A A A A A A* § · ¨ ¸ ¨ ¸ ¨ ¸ ¨ ¸ © ¹ *MD* ( ) SRVLWLRQ FDOFXODWH 0'V *n z z zn AA A AA A* § · ¨ ¸ © ¹ %HQFKPDUN 7UDLQLQJ *n n m m mn BB B BB B BB B* § · ¨ ¸ ¨ ¸ ¨ ¸ ¨ ¸ © ¹ \$ *B C* \$ \$ \$M *n z z zn Aj Aj Aj Aj Aj Aj* § · ¨ ¸ © ¹ *n z z zn AA A AA A* § · ¨ ¸ ¨ ¸ © ¹ *n z z zn CC C CC C* § · ¨ ¸ © ¹ & & &M *n z z zn Cj Cj Cj Cj Cj Cj* § · ¨ ¸ © ¹ *n z z zn CC C CC C* § · ¨ ¸ © ¹ \$ *j MD MD MD* § · ¨ ¸ ¨ ¸ ¨ ¸ ¨ ¸ © ¹ UHSHDWWR *A Aj*

**Figure 5.** The third step of aMMTS.

*Sensors* **2019**, *19*, 26

Forth, the proper sequence of SNR gain is chosen by a function, which is according to the minimum MD. If two labels of data corresponding to minimum MD are the same, the recognition result is right, and accumulate the number of right recognition results. *S* is the recognition result. This step is shown in Figure 6.

**Figure 6.** The forth step of aMMTS.

Finally, the recognition result is verified if there is a unique optimal recognition result, and the SNR gain's position in the sequence is the feature's order. If there are several optimal recognition results, repeat step 3 to recalculate the result. Afterwards, a set of sequences with the best recognition effect is determined as the feature sequence, which solves the self-adaptation problem of thresholds.

### **3. Results**

In this paper, the experimental data are from Case Western Reserve University Bearing Data Center. This experiment involved three different faults that occurred on three components: inner race, outer race and rolling element. The vibration signals were acquired under four different speeds: 1797 r/min, 1772 r/min, 1750 r/min, and 1730 r/min, and the sampling frequency was set to 12 kHz. To demonstrate the aMMTS, this study randomly selected the data in the dataset under the defect of 0.07 inches. The number of samples are shown in Table 1.


**Table 1.** The number of samples.

There were 2192 samples; 548 for inner race, 548 for outer race, 548 for rolling element and 548 for normal. The data are divided into three parts: training data, validation data and test data. In order to avoid the overfitting caused by the training data, training samples were used to construct MS, generate the SNR gain and calculate MDs by using the sequences of SNR gains, and were divided into three parts, with one of the parts set as benchmark group. To avoid the over-fitting problem, group A was used to construct MS and generate the SNR gain, and group B were used to calculate the MDs with the sequences of SNR gain and identify faults.

Validation samples were used to verify the recognition result if there exists the same minimum MDs, and the sequence was selected according to the best result.

Test samples were used to validate the proposed method.

### *3.1. Signal Decomposition by Using VMD and Wavelet Denosing*

Above all, this study employed wavelet denoising to remove the noise from the raw signals. First, the Daubechies 5 (db5) was used to decompose the signal, and obtained the wavelet decomposition vector and the bookkeeping vector. Second, thresholds wavelet coefficient was calculated by setting the detail vector which would be compressed as [1–3] and the vector which is the corresponding percentages of lower coefficients as [100,90,80], and using the wavelet decomposition vector and the bookkeeping vector. Lastly, the thresholds, Daubechies 5 (db5) and decomposed signals were used to reconstruct the denoising signal. Then the VMD was used to decompose the signal, and was needed to give the preset IMF component number *K* and penalty parameter *α* which constrained the moderate bandwidth. The value of *α* toke the default value 1024, the value of *K* was 8. An example is shown in Figure 7.

**Figure 7.** The intercepted signal.

IMFs are as shown in Figure 8.

**Figure 8.** The IMFs.

### *3.2. Feature Extraction by Using SVD*

SVD was used to analyze the IMFs. After the signal decomposition, the IMF matrix was decomposed by SVD, and obtained singular value vectors. The singular value vectors were considered as features and formed the feature matrix. Then, the feature matrix was used to diagnose the fault by aMMTS. To avoid the over-fitting problem, the features were divided into training samples, validation samples and test samples. The features of the above IMFs of those were shown in Table 2.


**Table 2.** The feature of IMFs.

The features obtained by SVD are shown in Table 3.


### *3.3. Fault Diagnosis Using aMMTS*

After the feature extraction, the aMMTS was used to identify and diagnose fault modes. The steps of aMMTS are as follow:

Firstly, the MS of training and benchmark were constructed, the eight-factor and two-level orthogonal array is shown in Table 4, and the MS based on Table 2 is shown in Table 5;


**Table 4.** The eight-factor and two-level orthogonal array.


**Table 5.** The Mahalanobis space (MS) based on the inner race.

Secondly, the MD was calculated, and SNR gain was also obtained by benchmark samples and training samples. The SNR gain is shown in Table 6;


**Table 6.** The signal-to-noise ratio (SNR) gain of features.

Thirdly, the MDs between the benchmark samples and validation samples were calculated by using the ascending and descending order of SNR;

Fourthly, the validation samples were used to verify the correctness of feature selection which existed more than one smallest MD;

Fifthly, the best sequence was chosen and set as the sequence of features.

Lastly, the best sequence was used to identify the test samples. Took the benchmark is outer race as the example shown in Figure 9.

**Figure 9.** The classification result of outer race.

Finally, the test sample was used to test the result of the method, and the benchmarks were inner race, rolling element, outer race and normal. The results are shown in Table 7 and the MDs between benchmark and test sample are shown in Figure 10.

**(c) Figure 10.** *Cont*.

**(d)** 

**Figure 10.** The Mahalanobis distances (MDs) between benchmark and testing data: (**a**) The MDs between benchmark (Inner race) and testing data; (**b**) The MDs between benchmark (Outer race) and testing data; (**c**) The MDs between benchmark (Rolling element) and testing data; (**d**) The MDs between benchmark (Normal) and testing data.

**Table 7.** The recognition result.


As shown in Table 7 and Figure 10, this method accurately classified and diagnosed the fault of the bearing by using the different benchmarks. The recognition results of normal and outer race reached 100%. However, it is not accurate enough to diagnose the fault of inner race and rolling element. However, in the normal and the fault of inner race, it is effective in industrial application.

### **4. Discussion**

Rolling element bearings are one of the most frequently used components in rotating machineries. This paper presents the method based on the wavelet denoising VMD-SVD-aMMTS to diagnose the fault of bearings under the variable conditions. Firstly, VMD is used to decompose the signal. Secondly, SVD is used to extract the feature. The adaptive aMMTS uses the feature sequences and multi-benchmarks to overcome the drawback of MTS for adaptive feature selection, multi-classification and over-fitting. The experimental result shows that the method could accurately diagnose faults effectively.

However, in the actual situation, there is an imbalance between fault data and normal data. In this method, aMMTS lacks research on the imbalanced study. The absence of faulty data may create a new problem, such over-fitting. Therefore, additional experiments under imbalanced data should be done to improve the method.

**Author Contributions:** N.W. collected and analyzed the data, made charts and diagrams, conceived and performed the experiments and wrote the paper; Z.W. and L.J. conceived the structure, provided guidance and modified the manuscript; X.C. analyzed the data and contributed analysis tools. Y.Q. provided guidance. Y.Z. revised the reviews.

**Funding:** This research was funded by National Key R&D Program of China grant number (No.2016YFB1200402).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Variational Bayesian Based Adaptive Shifted Rayleigh Filter for Bearings-Only Tracking in Clutters**

### **Jing Hou , Yan Yang \* and Tian Gao**

School of Electronic and Information, Northwestern Polytechnical University, Xi'an 710072, China; jhou0825@nwpu.edu.cn (J.H.); tiangao@nwpu.edu.cn (T.G.)

**\*** Correspondence: yangyan7003@nwpu.edu.cn

Received: 11 February 2019; Accepted: 26 March 2019; Published: 28 March 2019

**Abstract:** This paper considers bearings-only target tracking in clutters with uncertain clutter probability. The traditional shifted Rayleigh filter (SRF), which assumes known clutter probability, may have degraded performance in challenging scenarios. To improve the tracking performance, a variational Bayesian-based adaptive shifted Rayleigh filter (VB-SRF) is proposed in this paper. The target state and the clutter probability are jointly estimated to account for the uncertainty in clutter probability. Performance of the proposed filter is evaluated by comparing with SRF and the probability data association (PDA)-based filters in two scenarios. Simulation results show that the proposed VB-SRF algorithm outperforms the traditional SRF and PDA-based filters especially in complex adverse scenarios in terms of track continuity, track accuracy and robustness with a little higher computation complexity.

**Keywords:** bearings-only tracking; clutter; variational Bayesian; Shifted Rayleigh Filter

### **1. Introduction**

Bearings-only target tracking is to estimate the current position and velocity of a target using only the noise-corrupted bearing measurements from one or multiple observer platforms. It is an important tracking problem that arises in both military and civilian applications, such as underwater sonar tracking, bistatic radar, air traffic control and computer vision [1].

Because of the intrinsic nonlinearities in the measurement models, it is difficult to acquire an optimal solution of this problem. Several suboptimal algorithms have been developed for bearings-only tracking in the literature. The extended Kalman filter (EKF) [2] in the Cartesian coordinate system is an early attempt. However, it is easy to diverge. To improve the stability of the EKF, the bearings-only tracking problem was formulated in modified polar coordinates, resulting in the modified polar coordinate EKF (MPEKF) [3]. However, it requires good initialization to guarantee convergence. The well-known pseudo-linear estimator (PLE) [4] was also developed to solve the bearings-only tracking problem. However, it gives a biased estimate at long ranges. In recent years, some bias compensation techniques were developed to improve the performance of PLE [5–7]. In addition, more sophisticated nonlinear Kalman filtering algorithms, such as unscented Kalman filter (UKF) [8,9], cubature Kalman filter (CKF) [1] and particle filter (PF) [10,11] were applied for bearings-only target tracking. PF can provide good performance but at the price of heavy computation load. Noteworthily, Clark et al. [12,13] proposed a novel shifted Rayleigh filter (SRF) for bearings-only tracking, which is still based on the approximation of conditional expectations but with novel feature which is performing a calculation to exploit the essential structure of the nonlinearities in a new way. It is shown to exhibit similar performance to the PF in certain challenging scenarios with much lower computational complexity [14].

However, these algorithms do not consider the impact of clutter which makes the bearings-only tracking problem more intractable. Apparently, the classical treatment of clutter in target tracking

problem can be extended to the bearings-only tracking problem. For example, Reference [15] integrated the maximum entropy fuzzy probabilistic data association (MEFPDA) with the square-root cubature Kalman filter (SCKF) to deal with the clutter in bearings-only tracking. Clark et al. also included the effect of clutter measurements into the SRF algorithm in [12]. However, in the classical algorithms, the clutter probability is usually assumed known and constant, which maybe time-varying or hard to determine in advance especially in adverse scenarios. Use of incorrect value of the clutter probability may lead to track accuracy degradation even track loss. A straight-forward idea is to account for the unknown clutter probability in the process of estimation of the state. That is, we need to solve the problem of bearings-only target tracking with uncertain parameter of clutter probability.

As we know, the Bayesian approach is the most general approach of solving the problem with uncertain parameters. However, it is not trivial to get the analytical solution for most Bayesian approaches due to complex nonlinear probability density function or high dimension of integration. Recently, the variational Bayesian (VB)-based adaptive filters [16–18] have drawn extensive attention, which utilize a new simpler, analytically tractable distribution to approximate the true posterior distribution so that avoiding direct calculations of complex integrals. Its adaptive strategy has a strong ability of tracking time-varying parameters. Therefore, we adopt the VB method to jointly estimate the target state and the clutter probability within the framework of SRF in this paper for bearings-only tracking in clutters. By establishing a conjugate exponential model for clutter probability and data association indicator, the proposed filter is derived using the iterative filtering framework. The tracking performance of the proposed VB-SRF is evaluated by comparing with SRF, PDA-SCKF [15] and MEFPDA-SCKF [15] via two simulation examples. It shows that the proposed filter outperforms the traditional SRF and two PDA-SCKF-based filters in complex mismatched scenarios in terms of track continuity and track accuracy but at the cost of higher computation complexity.

The remainder of the paper is organized as follows. Section 2 gives the problem formulation. In Section 3, the variational Bayesian filtering is described. Section 4 derives the VB-based adaptive SRF. Section 5 provides simulation results and performance evaluation of the proposed approach, followed by the conclusions in Section 6.

### **2. The Shifted Rayleigh Filter Algorithm**

### *2.1. The Bearing Model*

Considering the shifted Rayleigh filter (SRF) for bearings-only tracking in *R*2, the state equation and the measurement equation are described as [13]:

$$\mathbf{x}\_{k} \quad = \mathbf{F}\_{k-1}\mathbf{x}\_{k-1} + \mathbf{u}\_{k-1}^{s} + \mathbf{v}\_{k-1} \tag{1}$$

$$\begin{array}{lclcl}\mathbf{y}\_k &= \mathbf{H}\_k \mathbf{x}\_k + \mathbf{u}\_k^{\mathrm{m}} + \mathbf{w}\_k\\ \mathbf{b}\_k &= \boldsymbol{\Pi}(\mathbf{y}\_k) \end{array} \tag{2}$$

where, **x***<sup>k</sup>* is the state vector which describes the position and velocity of the target; **b***<sup>k</sup>* is the noisy bearing measurement. Π denotes the projection of the plane onto the unit circle. That is taking a 2-vector **<sup>y</sup>***<sup>k</sup>* into its normalized form **<sup>y</sup>***k*/||**y***k*||. Then **<sup>b</sup>***<sup>k</sup>* = (*sinθk*, *cosθk*)*T*, where *<sup>θ</sup><sup>k</sup>* is the bearing of the target position relative to the sensor platform. **F***k*−<sup>1</sup> and **H***<sup>k</sup>* are the state transition matrix and measurement matrix, respectively. **u***<sup>s</sup> <sup>k</sup>*−<sup>1</sup> and **<sup>u</sup>***<sup>m</sup> <sup>k</sup>* are the inputs to the system to increase the versatility of the model, for example, to reflect known perturbations to the dynamics and changes in sensor location; **v***k*−<sup>1</sup> is the Gaussian process noise with zero mean and covariance **Q***v*, and **w***<sup>k</sup>* is the Gaussian measurement noise with zero mean and covariance **Q***<sup>w</sup>* and independent of **v***k*.

The unusual point in the measurement model is that the noise **w***<sup>k</sup>* is modelled as additive noise present in an "augmented" measurement, **y***k*, of the Cartesian coordinates of target relative to the sensor platform, which is projected onto the plane to generate the actual bearing measurement **b***k*. It is different with the traditional "angle plus white noise" model, expressed as

$$\theta\_k = \arctan(d\_1/d\_2) + \epsilon\_k \tag{3}$$

where, *d*1, *d*<sup>2</sup> are the components of the displacement vector **d***<sup>k</sup>* = **H***k***x***<sup>k</sup>* + **u***<sup>m</sup> <sup>k</sup>* .  *<sup>k</sup>* is the sensor noise with Gaussian distribution N (0, *<sup>σ</sup>*2) and independent of the displacement vector **<sup>d</sup>***k*.

Actually, as explained in [13], the shifted Rayleigh bearing model (2) can be related with the traditional model (3) by making a variant on the shifted Rayleigh noise model

$$\begin{array}{ll} \mathbf{y}'\_k &= \mathbf{d}\_k + ||\mathbf{d}\_k|| \mathbf{e}' \\ \mathbf{b}'\_k &= \Pi(\mathbf{y}'\_k) \end{array} \tag{4}$$

where, **<sup>e</sup>** is an <sup>N</sup> (**0**, *<sup>σ</sup>*2**I**2×2) distributed random variable. The only difference with model (2) is the noise term **w***<sup>k</sup>* used to construct the augmented measurement **y***<sup>k</sup>* is replaced by ||**d***k*||**e** . ||**d***k*||**e** differs from **w***k*, but has identical first and second moments, and is uncorrelated with **d***k*.

The angular *θ <sup>k</sup>* of the modified vector bearing **b** *<sup>k</sup>* can be represented as

$$\theta\_k' = \arctan(d\_1/d\_2) + \epsilon\_k' \tag{5}$$

where  *<sup>k</sup>* is a zero mean random variable, restricted to [−*π*, *π*], independent of **d***k*, with density *ασ*(·):

$$a\_{\sigma}(\theta) = \frac{e^{-1/2\sigma^2}}{2\pi} \left( 1 + \sqrt{2\pi} \frac{\cos\theta}{\sigma} F\_{\text{normal}} \left( \frac{\cos\theta}{\sigma} \right) e^{1/2(\cos\theta/\sigma)^2} \right) \tag{6}$$

where, *Fnormal*(·) is the cumulative distribution function of a standard N (0, 1) variable.

Note that *θ* given by (5), is very close to the bearing *θ* in the standard model (3). The only difference is that the densities of the noise terms used in their construction are *ασ*(*θ*) and the normal N (0, *<sup>σ</sup>*2) density, respectively. Reference [13] plots the two densities for *<sup>σ</sup>*<sup>2</sup> ≤ 1. It shows the two density functions are virtually indistinguishable.

Given the bearing model, the SRF is to calculate the estimates of the conditional mean and covariance of the target state **x***k*, given measurements up to time *k*, **b**1:*k*. That is,

$$\mathbf{\dot{x}}\_{k|k} = E[\mathbf{x}\_k | \mathbf{b}\_{1:k}]\_\prime \quad \mathbf{P}\_{k|k} = cov[\mathbf{x}\_k | \mathbf{b}\_{1:k}] \tag{7}$$

The formulas for the SRF algorithm can be seen in [13].

### *2.2. The Treatment of Clutter*

Accounting for the effects of clutter on the bearing measurements, we represent a cluttered bearing measurement as

$$z\_k = (1 - r\_k)\theta\_k + r\_k l I\_k \tag{8}$$

where *Uk* is the bearing measurement of the clutter, which is assumed to be an uniform random variable on [−*π*, *π*]; *θ<sup>k</sup>* is the bearing measurement of the actual target from the sensor, *rk* is defined as data association indicator at time *k*.

$$r\_k = \begin{cases} 1 & \text{if the measurement is from cluster} \\ 0 & \text{if the measurement is from target} \end{cases} \tag{9}$$

The prior of *rk* is assumed independent of the previous data associations and can be described as

$$p(r\_k) = \begin{cases} \ \emptyset & \text{if} \quad r\_k = 1\\ 1 - \emptyset & \text{if} \quad r\_k = 0 \end{cases} \tag{10}$$

Then the likelihood of the measurement is

$$p(z\_k|\mathbf{x}\_k, r\_k) = \begin{cases} 1/2\pi & \text{if } \quad r\_k = 1\\ f(\theta\_k|\mathbf{x}\_k) & \text{if } \quad r\_k = 0 \end{cases} \tag{11}$$

It can be represented in two forms:

$$p(z\_k|\mathbf{x}\_k, r\_k) = \frac{1}{2\pi\mathfrak{r}}r\_k + f(\theta\_k|\mathbf{x}\_k)(1 - r\_k) \tag{12}$$

$$\mathbf{x} = (\frac{1}{2\pi})^{r\_k} f(\theta\_k | \mathbf{x}\_k)^{1 - r\_k} \tag{13}$$

Note that the two forms are used for different purpose in the later section. *f*(*θk*|**x***k*) is the likelihood of the target bearing measurement, whose expression is given in Appendix A.

Here, we need to explain the reason for using this representation (8) of the cluttered measurement. The information carried in (8) is identical with the projection measurement **z***<sup>b</sup> <sup>k</sup>* = [*sinzk*, *coszk*]. So the mean and covariance of the target state **x***<sup>k</sup>* conditioned on cluttered measurement **z***<sup>b</sup> <sup>k</sup>*, which are the aims needs to be achieved in the SRF, are equivalent with these conditioned on *zk*. Furthermore, this representation (8) is more convenient for calculation. Thus, *zk* instead of **z***<sup>b</sup> <sup>k</sup>* is adopted to represent the clutter measurement.

Based on the cluttered measurement model, given clutter probability *ξ*, the state estimate and covariance can be calculated as follows.

Suppose the density of **<sup>x</sup>***k*−<sup>1</sup> given measurements up to *<sup>k</sup>* − 1, *pk*−1|*k*−1(**x***k*−1), is normal with mean **<sup>x</sup>**¯ *<sup>k</sup>*−1|*k*−<sup>1</sup> and covariance **P¯** *<sup>k</sup>*−1|*k*−1. The posterior density of state **<sup>x</sup>***<sup>k</sup>* conditioned on cluttered bearing measurements **z**1:*<sup>k</sup>* is given as

$$p(\mathbf{x}\_k|\mathbf{z}\_{1:k}) = p(\mathbf{x}\_k|z\_k, \mathbf{x}\_{k-1} \sim p\_{k-1|k-1}(\mathbf{x}\_{k-1}))$$

$$= q\_k(0)p\_{k|k}(\mathbf{x}\_k) + q\_k(1)p\_{k|k-1}(\mathbf{x}\_k) \tag{14}$$

where *qk*(*i*) = *<sup>p</sup>*(*rk* = *<sup>i</sup>*|*zk*, **<sup>z</sup>**1:*k*−1), *pk*|*k*(**x***k*) is the non-normal density of **<sup>x</sup>***<sup>k</sup>* conditioned on *rk* = 0 and *zk*, or, equivalently, on *<sup>θ</sup>k*, and *pk*|*k*−1(**x***k*) is the density of **<sup>x</sup>***<sup>k</sup>* when there is no target measurement at time *k*. It is normal with mean **x**ˆ *<sup>k</sup>*|*k*−<sup>1</sup> and covariance **P***k*|*k*−1. Thus, the state estimate and covariance at time *k* can be obtained as

$$\dot{\mathbf{x}}\_{k|k} = E\left[\mathbf{x}\_k | z\_k, \mathbf{x}\_{k-1} \sim p\_{k-1|k-1}(\mathbf{x}\_{k-1})\right]$$

$$= q\_k(0)\hat{\mathbf{x}}\_{k|k} + q\_k(1)\hat{\mathbf{x}}\_{k|k-1} \tag{15}$$

$$\begin{split} \mathbb{P}\_{k|k} &= \text{cov}[\mathbf{x}\_{k}|z\_{k}, \mathbf{x}\_{k-1} \sim p\_{k-1|k-1}(\mathbf{x}\_{k-1})] \\ &= q\_{k}(0) \left( \mathbf{P}\_{k|k} + \left( \hat{\mathbf{x}}\_{k|k} - \bar{\mathbf{x}}\_{k|k} \right) \left( \hat{\mathbf{x}}\_{k|k} - \bar{\mathbf{x}}\_{k|k} \right)^{T} \right) \\ &+ q\_{k}(1) \left( \mathbf{P}\_{k|k-1} + \left( \hat{\mathbf{x}}\_{k|k-1} - \bar{\mathbf{x}}\_{k|k} \right) \left( \hat{\mathbf{x}}\_{k|k-1} - \bar{\mathbf{x}}\_{k|k} \right)^{T} \right) \end{split} \tag{16}$$

where **x**ˆ *<sup>k</sup>*|*<sup>k</sup>* and **P***k*|*<sup>k</sup>* are the state estimate and its covariance based on the actual target measurement. **x**ˆ *<sup>k</sup>*|*k*−<sup>1</sup> and **P***k*|*k*−<sup>1</sup> are the predicted target state estimate and covariance. All of these can be obtained using the basic formulas of SRF.

The conditional densities *qk*(0) and *qk*(1) are given by the equations

$$q\_k(0) = c\_k (1 - \zeta) f\_k(\theta\_k | \mathbf{z}\_{1:k-1}) \tag{17}$$

$$q\_k(1) = \frac{c\_k \zeta}{2\pi} \tag{18}$$

where *ck* is the normalizing constant and *fk*(*θk*|**z**1:*k*−1) is the density of the actual target bearing *<sup>θ</sup><sup>k</sup>* conditioned on measurements **z**1:*k*−1. The expression is given in Appendix B.

However, in a complex environment, the probability of clutter is time-varying or hard to determine in advance. In this case, the probability *ξ* is unknown, so the above formulas are not applicable. In this paper, we resort the VB method to find the joint posterior density of **x***<sup>k</sup>* and *ξ* so as to account for the uncertainty in clutter probability.

### **3. Variational Bayesian Filtering**

In this section, we first review the conjugate exponential (CE) model, and then gives the variational Bayesian solution of the CE model.

### *3.1. Conjugate Exponential Model*

Given measurements **<sup>z</sup>**1:*k*−1, the posterior of the system state *<sup>p</sup>*(**x***k*−1|**z**1:*k*−1) and the posterior of the parameter *<sup>p</sup>*(**r***k*−1|**z**1:*k*−1), we assume the complete-data likelihood in the exponential family:

$$p(\mathbf{x}\_k, \mathbf{z}\_k | \mathbf{r}\_k, \mathbf{z}\_{1:k-1}) = g(\mathbf{r}\_k) f(\mathbf{x}\_k, \mathbf{z}\_k) e^{\oint (\mathbf{r}\_k)^T u(\mathbf{x}\_k, \mathbf{z}\_k)} \tag{19}$$

where, *φ*(**r***k*) is the vector of natural parameters **r***k*, *u* and *f* are known functions, and *g* is a normalization constant:

$$g(\mathbf{r}\_k)^{-1} = \int f(\mathbf{x}\_k, \mathbf{z}\_k) e^{\phi(\mathbf{r}\_k)^T \boldsymbol{u}(\mathbf{x}\_k, \mathbf{z}\_k)} d\mathbf{x}\_k d\mathbf{z}\_k$$

The parameter prior is conjugate to the complete-data likelihood:

$$p(\mathbf{r}\_k | a\_k^-, \boldsymbol{\beta}\_k^-) = h(a\_k^-, \boldsymbol{\beta}\_k^-) g(\mathbf{r}\_k)^{\boldsymbol{\beta}\_k^-} e^{\boldsymbol{\phi}(\mathbf{r}\_k)^T a\_k^-} \tag{20}$$

where *α*− *<sup>k</sup>* and *β*<sup>−</sup> *<sup>k</sup>* are hyperparameters of the prior, and *h* is a normalization constant. Note the prior *p*(**r***k*|*α*<sup>−</sup> *<sup>k</sup>* , *β*<sup>−</sup> *<sup>k</sup>* ) is said to be conjugate to the likelihood *p*(**x***k*, **z***k*|**r***k*) if and only if the posterior

$$p(\mathbf{r}\_k | \boldsymbol{\alpha}\_k, \boldsymbol{\beta}\_k) \propto p(\mathbf{r}\_k | \boldsymbol{\alpha}\_k^-, \boldsymbol{\beta}\_k^-) p(\mathbf{x}\_k, \mathbf{z}\_k | \mathbf{r}\_k)$$

is of the same parametric form as the prior. Then we call models that satisfy Equations (19) and (20) conjugate-exponential.

### *3.2. VB Approximation Method*

Applying Bayes' rule, we have the joint posterior of **x***<sup>k</sup>* and *rk* as

$$p(\mathbf{x}\_k, \mathbf{r}\_k | \mathbf{z}\_{1:k}) \propto p(\mathbf{x}\_k, \mathbf{z}\_k | \mathbf{r}\_k, \mathbf{z}\_{1:k}) p(\mathbf{r}\_k | \boldsymbol{\alpha}\_k^-, \boldsymbol{\beta}\_k^-) \tag{21}$$

The analytic solution to (21) would be difficult to calculate. Here, we use the VB method to approximate the true posterior distribution with a product of tractable marginal posteriors [17].

$$p(\mathbf{x}\_{k\prime}, \mathbf{r}\_k | \mathbf{z}\_{1:k}) \approx Q(\mathbf{x}\_k, \mathbf{r}\_k) = Q\_x(\mathbf{x}\_k) Q\_r(\mathbf{r}\_k) \tag{22}$$

where *Qx*(**x***k*) and *Qr*(**r***k*) are unknown approximating marginal densities of **x***<sup>k</sup>* and **r***k*.

The basic idea of VB approximation is to minimize the Kullback- Leibler (KL) divergence between the approximating posterior and the true posterior:

$$KL\left[Q\_{\mathbf{x}}(\mathbf{x}\_{k})Q\_{\mathbf{r}}(\mathbf{r}\_{k})||p(\mathbf{x}\_{k},\mathbf{r}\_{k}|\mathbf{z}\_{1:k})\right] = \int Q\_{\mathbf{x}}(\mathbf{x}\_{k})Q\_{\mathbf{r}}(\mathbf{r}\_{k}) \times \log\left(\frac{Q\_{\mathbf{x}}(\mathbf{x}\_{k})Q\_{\mathbf{r}}(\mathbf{r}\_{k})}{p(\mathbf{x}\_{k},\mathbf{r}\_{k}|\mathbf{z}\_{1:k})}\right) d\mathbf{x}\_{k} d\mathbf{r}\_{k}\tag{23}$$

Given the measurements **z**1:*k*, we can minimize the KL divergence with respect to the probability densities *Qx*(**x***k*) and *Qr*(**r***k*) in turn, while keeping the other fixed. Then, the following equations can be given as:

$$Q\_{\mathbf{x}}(\mathbf{x}\_{k}) \propto \exp\left(\langle \ln p(\mathbf{x}\_{k}, \mathbf{r}\_{k}, \mathbf{z}\_{k} | \mathbf{z}\_{1:k-1}) \rangle\_{\mathbf{r}\_{k}}\right) \tag{24}$$

$$Q\_r(\mathbf{r}\_k) \propto \exp\left( \langle \ln p(\mathbf{x}\_{k\prime}, \mathbf{r}\_{k\prime} | \mathbf{z}\_{1:k-1}) \rangle\_{\mathbf{x}\_k} \right) \tag{25}$$

where ·**x***<sup>k</sup>* and ·**r***<sup>k</sup>* denote the expectations with respect to *Qx*(**x***k*) and *Qr*(**r***k*), respectively. Obviously, it is not an explicit solution since the distribution of each parameter is dependent on the other and neither distributions is known. The mechanism of VB method is to firstly give the initial values of the parameters and then use expectation-maximum (EM) algorithm to iteratively calculate *Qx*(**x***k*) and *Qr*(**r***k*) until convergence. For the above CE models, *Qx*(**x***k*) and *Qr*(**r***k*) can be obtained from the following procedure [19]:

(1) The VB expectation step yields:

$$Q\_x(\mathbf{x}\_k) \propto f(\mathbf{x}\_k, \mathbf{z}\_k) e^{\langle \phi(\mathbf{r}\_k) \rangle\_{\mathbf{r}\_k}^T \mu(\mathbf{x}\_k, \mathbf{z}\_k)} = p(\mathbf{x}\_k | \mathbf{z}\_k, \langle \phi(\mathbf{r}\_k) \rangle\_{\mathbf{r}\_k}) \tag{26}$$

(2) The VB maximization step yields that *Qr*(**r***k*) is conjugate and of the form

$$Q\_{\mathbf{r}}(\mathbf{r}\_k) = h(\boldsymbol{a}\_k, \boldsymbol{\beta}\_k^- \mathbf{g}(\mathbf{r}\_k)^{\beta\_k} e^{\boldsymbol{\Phi}(\mathbf{r}\_k)^T \boldsymbol{a}\_k} \tag{27}$$

where, *α<sup>k</sup>* and *β<sup>k</sup>* are the hyper-parameters, and

$$
\mu\_k = \alpha\_k^- + \langle \mu(\mathbf{x}\_{k'}, \mathbf{z}\_k) \rangle\_{\mathbf{x}\_k} \tag{28}
$$

$$
\beta\_k = \beta\_k^- + n
$$

where *n* is the dimension of the measurement.

### **4. VB Based Adaptive Shifted Rayleigh Filter with Unknown Clutter Probability**

Considering the system model (1) and the measurement model (8) described in Section 2, we adopt the VB method within the SRF framework to get the joint estimation of the target state and the clutter probability.

The core is to determine the posterior approximation *Q*(**x***k*,*rk*, *ξ*). Assume factorization *Q*(**x***k*,*rk*, *ξ*) ≈ *Qx*(**x***k*)*Q*(*rk*, *ξ*), we can obtain *Qx*(**x***k*), *Qr*(*rk*) and *Q<sup>ξ</sup>* (*ξk*) at each time *k* through the following procedure.

(1) Optimization of *Qx*(**x***k*) for fixed *Q*(*rk*, *ξ*).

First, by using the first form of the measurement likelihood (12), the complete-data likelihood is presented as

$$\begin{split} p(\mathbf{x}\_{k}, z\_{k} | r\_{k}, \mathbf{z}\_{1:k-1}) &= p(z\_{k} | \mathbf{x}\_{k}, r\_{k}) p(\mathbf{x}\_{k} | r\_{k}, \mathbf{z}\_{1:k-1}) \\ &= [\frac{1}{2\pi} r\_{k} + f(\theta\_{k} | \mathbf{x}\_{k})(1 - r\_{k})] N(\mathbf{x}\_{k}; \mathbf{z}\_{k|k-1}, \mathbf{P}\_{k|k-1}) \end{split} \tag{30}$$

Then according to (24), we can get

$$\begin{split} Q\_x(\mathbf{x}\_k) &\approx \exp\{ \langle \text{lnp}(z\_k, \mathbf{x}\_k | r\_k, \mathbf{z}\_{1:k-1}) \rangle\_{r\_k \bar{r}} \} \\ &= [\frac{1}{2\pi} \langle r\_k \rangle\_{r\_k} + f(\theta\_k | \mathbf{x}\_k) (1 - \langle r\_k \rangle\_{l^q\_k}) ] \mathcal{N}(\mathbf{x}\_k; \mathbf{\hat{x}}\_{k|k-1}, \mathbf{P}\_{k|k-1}) \\ &= \frac{1}{2\pi} \mathcal{N}(\mathbf{x}\_k; \mathbf{\hat{x}}\_{k|k-1}, \mathbf{P}\_{k|k-1}) \langle r\_k \rangle\_{r\_k} + f(\theta\_k | \mathbf{x}\_k) \mathcal{N}(\mathbf{x}\_k; \mathbf{\hat{x}}\_{k|k-1}, \mathbf{P}\_{k|k-1}) (1 - \langle r\_k \rangle\_{r\_k}) \\ &\approx \frac{1}{2\pi} \mathcal{N}(\mathbf{x}\_k; \mathbf{\hat{x}}\_{k|k-1}, \mathbf{P}\_{k|k-1}) \langle r\_k \rangle\_{r\_k} + \mathcal{N}(\mathbf{x}\_k; \mathbf{\hat{x}}\_{k|k}, \mathbf{P}\_{k|k}) f(\theta\_k | \mathbf{z}\_{1:k-1}) (1 - \langle r\_k \rangle\_{r\_k}) \end{split} \tag{31}$$

The approximation sign in (32) is because the following:

$$\begin{split} f\left(\theta\_{k}|\mathbf{x}\_{k}\right) \mathcal{N}(\mathbf{x}\_{k}; \mathbf{\hat{x}}\_{k|k-1}, \mathbf{P}\_{k|k-1}) &= p(\theta\_{k}|\mathbf{x}\_{k}, \mathbf{z}\_{1:k-1}) p(\mathbf{x}\_{k}|\mathbf{z}\_{1:k-1}) \\ &= p(\mathbf{x}\_{k}|\theta\_{k'}, \mathbf{z}\_{1:k-1}) p(\theta\_{k}|\mathbf{z}\_{1:k-1}) \\ &\approx \mathcal{N}(\mathbf{x}\_{k}; \mathbf{\hat{x}}\_{k|k'}, \mathbf{P}\_{k|k}) f(\theta\_{k}|\mathbf{z}\_{1:k-1}) \end{split} \tag{33}$$

where *<sup>f</sup>*(*θk*|**z**1:*k*−1) is derived in Appendix B.

Comparing (32) with the posterior density (14) of SRF, the difference lies in the weights. Except for a normalization constant, the clutter probability *ξ* used before in (14) has been replaced by *rkrk* in (32), which is updated online.

(2) Optimization of *Q*(*rk*, *ξ*) for fixed *Qx*(**x***k*).

We use the VB method again by factorizing *Q*(*rk*, *ξ*) ≈ *Qr*(*rk*)*Q<sup>ξ</sup>* (*ξ*). Assume the conjugate prior of *rk* is binomial distributed with parameter *ξ*. That is,

$$p(r\_k | \xi) = \xi^{r\_k} (1 - \xi)^{1 - r\_k} \tag{34}$$

and *p*(*ξ*) follows beta distribution with parameters *α*<sup>1</sup> and *α*2.

$$p(\mathfrak{J}; \mathfrak{a}\_1, \mathfrak{a}\_2) = \frac{1}{B(\mathfrak{a}\_1, \mathfrak{a}\_2)} \mathfrak{J}^{\mathfrak{a}\_1 - 1} (1 - \mathfrak{J})^{\mathfrak{a}\_2 - 1} \tag{35}$$

where *B*(*α*1, *α*2) = Γ(*α*1)Γ(*α*2)/Γ(*α*<sup>1</sup> + *α*2).

To have the form of (19), the complete-data likelihood *<sup>p</sup>*(**x***k*, *zk*|*rk*, **<sup>z</sup>**1:*k*−1) is re-derived using the second form (13) of *p*(*zk*|**x***k*,*rk*) as

$$\begin{split} p(\mathbf{x}\_k, z\_k | r\_k, \mathbf{z}\_{1:k-1}) &= (1/2\pi)^{r\_k} [f(\theta\_k | \mathbf{x}\_k)]^{1-r\_k} \mathcal{N}(\mathbf{x}\_k; \mathbf{\hat{x}}\_{k|k-1}, \mathbf{P}\_{k|k-1}) \\ &= f(\theta\_k | \mathbf{x}\_k) \mathcal{N}(\mathbf{x}\_k; \mathbf{\hat{x}}\_{k|k-1}, \mathbf{P}\_{k|k-1}) \exp\{r\_k [\ln(1/2\pi) - \ln[f(\theta\_k | \mathbf{x}\_k)]]\} \end{split} \tag{36}$$

Rewriting the conjugate prior of *rk* in the form of (20), we can get

$$\begin{split} p(r\_k|\underline{\xi}) &= \underline{\xi}^{r\_k} (1 - \underline{\xi})^{1 - r\_k} \\ &= (1 - \underline{\xi}) \exp[r\_k \ln(\frac{\underline{\xi}}{1 - \underline{\xi}})] \end{split} \tag{37}$$

Then applying (27), *Qr*(*rk*) can be obtained as

$$Q\_r(r\_k) = (1 - \eta\_k) \exp[r\_k \ln(\frac{\eta\_k}{1 - \eta\_k})] \tag{38}$$

where, *η<sup>k</sup>* is the hyper-parameter and updated as

$$\ln(\frac{\eta\_k}{1-\eta\_k}) = \langle \ln(\frac{\mathcal{J}}{1-\mathcal{J}}) \rangle\_{\tilde{\mathcal{J}}} + \langle [\ln(1/2\pi) - \ln[f(\theta\_k|\mathbf{x}\_k)]] \rangle\_{\mathbf{x}\_k} \tag{39}$$

Likewise, rewriting the conjugate prior of *ξ* in the form of (20), we can get

$$\begin{split} p(\xi) &= \frac{1}{B(\mathfrak{a}\_1, \mathfrak{a}\_2)} \xi^{\mathfrak{a}\_1 - 1} (1 - \xi)^{\mathfrak{a}\_2 - 1} \\ &= \frac{1}{B(\mathfrak{a}\_1, \mathfrak{a}\_2)} (1 - \xi)^{\mathfrak{a}\_1 + \mathfrak{a}\_2 - 2} (1 - \xi)^{-(\mathfrak{a}\_1 - 1)} \xi^{\mathfrak{a}\_1 - 1} \\ &= \frac{1}{B(\mathfrak{a}\_1, \mathfrak{a}\_2)} (1 - \xi)^{\mathfrak{a}\_1 + \mathfrak{a}\_2 - 2} \exp[(\mathfrak{a}\_1 - 1) \ln(\frac{\xi}{1 - \frac{\mathfrak{a}}{\xi}})] \end{split} \tag{40}$$

Then applying (27), *Q<sup>ξ</sup>* (*ξ*) can be obtained as

$$Q\_{\tilde{\xi}}(\tilde{\xi}) = \frac{1}{B(a\_1', a\_2')} (1 - \tilde{\xi})^{a\_1' + a\_2' - 2} \exp[(a\_1' - 1) \ln(\frac{\tilde{\xi}}{1 - \tilde{\xi}})] \tag{41}$$

with hyper-parameters

$$
\alpha\_1' = \alpha\_1 + \langle r\_k \rangle\_{r\_k} \tag{42}
$$

$$
\pi\_2' = \pi\_2 + n - \langle r\_k \rangle\_{r\_k} \tag{43}
$$

where *n* is the dimension of the measurement.

According to the approximated posteriors of *Qr*(*rk*) and *<sup>Q</sup><sup>ξ</sup>* (*ξ*), *rkrk* and *ln*( *<sup>ξ</sup>* <sup>1</sup>−*<sup>ξ</sup>* )*<sup>ξ</sup>* can be obtained as:

$$
\langle r\_k \rangle\_{r\_k} = \eta\_k \tag{44}
$$

$$\langle \ln(\frac{\tilde{\xi}}{1-\tilde{\xi}}) \rangle\_{\tilde{\xi}} = \psi(a\_1') - \psi(a\_2') \tag{45}$$

where *ψ*(·) is the digamma function.

Taking expectation and covariance on the posterior *Qx*(**x***k*), the conditional mean and covariance of the target state can then be obtained. We summarize the entire filtering procedure of the VB-based SRF (VB-SRF) in Algorithm 1.

**Algorithm 1** : VB-SRF. **(1) Initialization**: **<sup>x</sup>**¯0|0, **<sup>P</sup>**¯ <sup>0</sup>|0, **<sup>Q</sup>***v*, **<sup>Q</sup>***w*, *<sup>η</sup>*0, *<sup>α</sup>*1,0, *<sup>α</sup>*2,0 **(2) Prediction: <sup>x</sup>**<sup>ˆ</sup> *<sup>k</sup>*|*k*−<sup>1</sup> <sup>=</sup> **<sup>F</sup>***k*−1**x**¯ *<sup>k</sup>*−1|*k*−<sup>1</sup> <sup>+</sup> **<sup>u</sup>***<sup>s</sup> k*−1 **<sup>P</sup>***k*|*k*−<sup>1</sup> <sup>=</sup> **<sup>F</sup>***k*−1**P**¯ *<sup>k</sup>*−1|*k*−1**F***<sup>T</sup> <sup>k</sup>*−<sup>1</sup> <sup>+</sup> **<sup>Q</sup>***<sup>v</sup>* **<sup>S</sup>***<sup>k</sup>* <sup>=</sup> **<sup>H</sup>***k***P***k*|*k*−1**H***<sup>T</sup> <sup>k</sup>* <sup>+</sup> **<sup>Q</sup>***<sup>m</sup> k <sup>η</sup>k*|*k*−<sup>1</sup> <sup>=</sup> *ρηk*−1, *<sup>α</sup>*1,*k*|*k*−<sup>1</sup> <sup>=</sup> *ρα*1,*k*−1, *<sup>α</sup>*2,*k*|*k*−<sup>1</sup> <sup>=</sup> *ρα*2,*k*−<sup>1</sup> where *ρ* is the scale factor and 0 < *ρ* ≤ 1. **(3) Update:** the update of VB-SRF utilizes iterate filtering framework. **(3.a) First set**: **x**¯ (0) *<sup>k</sup>*|*<sup>k</sup>* <sup>=</sup> **<sup>x</sup>**<sup>ˆ</sup> *<sup>k</sup>*|*k*−1, **<sup>P</sup>**¯ (0) *<sup>k</sup>*|*<sup>k</sup>* <sup>=</sup> **<sup>P</sup>***k*|*k*−1, *<sup>η</sup>*(0) *<sup>k</sup>* <sup>=</sup> *<sup>η</sup>k*|*k*−1, *<sup>α</sup>*(0) 1,*<sup>k</sup>* <sup>=</sup> *<sup>α</sup>*1,*k*|*k*−1, *<sup>α</sup>*(0) 2,*<sup>k</sup>* <sup>=</sup> *<sup>α</sup>*2,*k*|*k*−<sup>1</sup> **(3.b) Calculate state estimation and its covariance using SRF when the measurement is from the target**: **<sup>K</sup>***<sup>k</sup>* <sup>=</sup> **<sup>P</sup>***k*|*k*−1**H***<sup>T</sup> <sup>k</sup>* **<sup>S</sup>**−<sup>1</sup> *k ε<sup>k</sup>* = (**b***<sup>T</sup> <sup>k</sup>* **<sup>S</sup>**−<sup>1</sup> *<sup>k</sup>* **<sup>b</sup>***k*)−1/2**b***<sup>T</sup> <sup>k</sup>* **<sup>S</sup>**−<sup>1</sup> *<sup>k</sup>* (**H***k***X**<sup>ˆ</sup> *<sup>k</sup>*|*k*−<sup>1</sup> <sup>+</sup> **<sup>u</sup>***<sup>m</sup> k* ) *γ<sup>k</sup>* = (**b***<sup>T</sup> <sup>k</sup>* **<sup>S</sup>**−<sup>1</sup> *<sup>k</sup>* **<sup>b</sup>***k*)−1/2*ρn*(*εk*) *δ<sup>k</sup>* = (**b***<sup>T</sup> <sup>k</sup>* **<sup>S</sup>**−<sup>1</sup> *<sup>k</sup>* **<sup>b</sup>***k*)−1/2[<sup>2</sup> <sup>+</sup> *<sup>ε</sup>kρ*2(*εk*) <sup>−</sup> *<sup>ρ</sup>*<sup>2</sup> <sup>2</sup>*εk*] *<sup>ρ</sup>*2(*εk*) = *<sup>ε</sup><sup>k</sup> <sup>e</sup>* <sup>−</sup>*ε*<sup>2</sup> *<sup>k</sup>*/2+ <sup>√</sup>2*π*(*ε*<sup>2</sup> *<sup>k</sup>*+1)*Fnormal*(*ε<sup>k</sup>* ) *e* <sup>−</sup>*ε*<sup>2</sup> *<sup>k</sup>*/2+ <sup>√</sup>2*π*(*ε<sup>k</sup>* )*Fnormal*(*ε<sup>k</sup>* ) **<sup>x</sup>**<sup>ˆ</sup> *<sup>k</sup>*|*<sup>k</sup>* = (**<sup>I</sup>** <sup>−</sup> **<sup>K</sup>***k***H***k*)**x**<sup>ˆ</sup> *<sup>k</sup>*|*k*−<sup>1</sup> <sup>−</sup> **<sup>K</sup>***k***u***<sup>m</sup> <sup>k</sup>* + *γk***K***k***b***<sup>k</sup>* **<sup>P</sup>***k*|*<sup>k</sup>* = (**<sup>I</sup>** <sup>−</sup> **<sup>K</sup>***k***H***k*)**P***k*|*k*−<sup>1</sup> <sup>+</sup> *<sup>δ</sup>k***K***k***b***k***b***<sup>T</sup> <sup>k</sup>* **<sup>K</sup>***<sup>T</sup> k* **(3.c) For** *j* = 1 : *N*, iterate the following *N* (*N* denotes iterated times) steps: • **Calculate the fused state estimation and its covariance**: **x**¯ (*j*) *<sup>k</sup>*|*<sup>k</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup>*π<sup>c</sup> <sup>η</sup>*(*j*−1) *<sup>k</sup>* **<sup>x</sup>**<sup>ˆ</sup> *<sup>k</sup>*|*k*−<sup>1</sup> <sup>+</sup> <sup>1</sup> *<sup>c</sup>* (<sup>1</sup> <sup>−</sup> *<sup>η</sup>*(*j*−1) *<sup>k</sup>* )*f*(*θk*|**z**1:*k*−1)**x**<sup>ˆ</sup> *<sup>k</sup>*|*<sup>k</sup>* **<sup>P</sup>**¯ (*j*) *<sup>k</sup>*|*<sup>k</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup>*π<sup>c</sup> <sup>η</sup>*(*j*−1) *<sup>k</sup>* (**P***k*|*<sup>k</sup>* + (**x**<sup>ˆ</sup> *<sup>k</sup>*|*<sup>k</sup>* <sup>−</sup> **<sup>x</sup>**¯ *<sup>k</sup>*|*k*)(**x**<sup>ˆ</sup> *<sup>k</sup>*|*<sup>k</sup>* <sup>−</sup> **<sup>x</sup>**¯ *<sup>k</sup>*|*k*)*T*) +<sup>1</sup> *<sup>c</sup>* (<sup>1</sup> <sup>−</sup> *<sup>η</sup>*(*j*−1) *<sup>k</sup>* )*f*(*θk*|**z**1:*k*−1)(**P***k*|*k*−<sup>1</sup> + (**x**<sup>ˆ</sup> *<sup>k</sup>*|*k*−<sup>1</sup> <sup>−</sup> **<sup>x</sup>**¯ (*j*) *<sup>k</sup>*|*k*)(**x**<sup>ˆ</sup> *<sup>k</sup>*|*k*−<sup>1</sup> <sup>−</sup> **<sup>x</sup>**¯ (*j*) *<sup>k</sup>*|*k*)*T*) where *c* = <sup>1</sup> <sup>2</sup>*<sup>π</sup> <sup>η</sup>*(*j*−1) *<sup>k</sup>* <sup>+</sup> *<sup>f</sup>*(*θk*)(<sup>1</sup> <sup>−</sup> *<sup>η</sup>*(*j*−1) *<sup>k</sup>* ) is a normalization term, and *<sup>f</sup>*(*θk*|**z**1:*k*−1) can be obtained using (A6). • **Update parameters**: *ln*( *<sup>η</sup>*(*j*) *k* <sup>1</sup>−*η*(*j*) *k* ) = *ψ*(*α* (*j*−1) 1,*<sup>k</sup>* ) − *ψ*(*α* (*j*−1) 2,*<sup>k</sup>* ) + *ln*(1/2*π*) − *lnf*(*θk*|**x**¯ (*j*) *k*|*k*) *α* (*j*) 1,*<sup>k</sup>* = *α* (*j*−1) 1,*<sup>k</sup>* <sup>+</sup> *<sup>η</sup>*(*j*) *k α* (*j*) 2,*<sup>k</sup>* = *α* (*j*−1) 2,*<sup>k</sup>* <sup>−</sup> *<sup>η</sup>*(*j*) *<sup>k</sup>* + 1 • **End for and set <sup>x</sup>**¯ *<sup>k</sup>*|*<sup>k</sup>* <sup>=</sup> **<sup>x</sup>**¯ (*N*) *<sup>k</sup>*|*<sup>k</sup>* , **<sup>P</sup>**¯ *<sup>k</sup>*|*<sup>k</sup>* <sup>=</sup> **<sup>P</sup>**¯ (*N*) *<sup>k</sup>*|*<sup>k</sup>* , *<sup>η</sup><sup>k</sup>* <sup>=</sup> *<sup>η</sup>*(*N*) *<sup>k</sup>* , *<sup>α</sup>*1,*<sup>k</sup>* <sup>=</sup> *<sup>α</sup>*(*N*) 1,*<sup>k</sup>* , *<sup>α</sup>*2,*<sup>k</sup>* <sup>=</sup> *<sup>α</sup>*(*N*) 2,*<sup>k</sup>* .

### **5. Simulation Results**

To evaluate the performance of the VB-SRF algorithm, two scenarios which are almost the same with these in [12,13] are utilized. The differences lie in the clutter probability in scenario 1 and the sensor tracks in scenario 2 which were not detailed in [12]. The two scenarios are very representative. In scenario 1, a maneuvering sensor is used in order to satisfy the condition of observability in bearings-only tracking. In scenario 2, multiple distributed sensors with large noise variance are utilized, which make the problem more challenging. The tracking performance of the VB-SRF algorithm was compared with the SRF algorithm, the MEFPDA-SCKF algorithm and PDA-SCKF algorithm in terms of track loss, track accuracy and computation complexity. The track loss is declared when the track error is large enough that making the filter diverge. Root mean square (RMS) error is used to show the track accuracy. In addition, the computation complexity is reflected by the computation time of each filter. The simulation codes can be downloaded through Github [20].

### *5.1. Scenario 1*

In scenario 1, a target moves along a horizontal track, with zero vertical displacement, according to a white noise acceleration model. The state of the target is represented as **x***<sup>k</sup>* = [*x*1,*k*, *x*˙1,*k*], where *x*1,*<sup>k</sup>* and *x*˙1,*<sup>k</sup>* are the horizontal distance and velocity at time *k*, respectively. The observer platform follows an approximately parallel track at a constant average speed. The horizontal and vertical displacements of the platform *x<sup>p</sup> <sup>k</sup>* = [*x<sup>p</sup>* 1,*k*, *<sup>x</sup><sup>p</sup>* 2,*k*] *<sup>T</sup>* are governed by the equations:

$$\mathfrak{x}\_{1,k}^p = 4k + \mathfrak{x}\_{1,k}^p \tag{46}$$

$$\mathfrak{x}\_{2,k}^p = \mathfrak{A}0 + \mathfrak{x}\_{2,k}^p \tag{47}$$

in which *x*˜ *p* 1,*<sup>k</sup>* and *x*˜ *p* 2,*<sup>k</sup>* are zero mean Gaussian white noise processes, both with variance *q* = 1.

The measurement is the angle (in radians) of the line-of-sight of this target from the platform. The sensor noise is Gaussian white noise with variance *σ*<sup>2</sup> = (0.05)2rad2 = 2.862deg2 . The true clutter probability is 0.8. Other parameters are detailed in [13]. The configuration of the observer platform and target is illustrated in Figure 1. The bearing measurement of the target is presented in Figure 2. It can be seen that there is no abrupt change.

**Figure 1.** Target-observer geometry in Scenario 1.

The clutter probability used in SRF is set as *pc* = 0.8, which is the same with the true clutter probability. The initial values of VB-SRF parameters are *η*<sup>0</sup> = 0.8, *α*1,0 = 2, *α*2,0 = 10. The clutter density *λ* is calculated using −*kln*(1 − *pc*)/2*π* in MEFPDA-SCKF and PDA-SCKF. Figure 3 presents the RMS target position errors of the four filters using 1000 Monte Carlo runs. It can be seen that SRF and VB-SRF have comparable performance under the correct clutter probability. MEFPDA-SCKF and PDA-SCKF have a little better tracking accuracy than SRF and VB-SRF.

**Figure 2.** The measurement of the target in Scenario 1.

**Figure 3.** RMS target position errors with correct clutter probability in Scenario 1.

In challenging scenarios, the clutter probability maybe unknown or time varying. The pre-set parameters are probably inaccurate. So here we set mistuned clutter probabilities for the four filters: (1) *pc* = *η*<sup>0</sup> = 0.7; (2) *pc* = *η*<sup>0</sup> = 0.5; (3) *pc* = *η*<sup>0</sup> = 0.3. The percentages of track losses in 1000 Monte Carlo runs are given in Table 1. Clearly, the VB-SRF algorithm outperforms the SRF algorithm due to fewer track losses. Moreover the proportion of track losses increases as the mistuning aggravates. Especially, there are 2.9% tracks are lost for SRF while only 0.1% tracks are lost for VB-SRF in the worst case. In addition, we can see that MEFPDA-SCKF and PDA-SCKF have no track loss in all the three mistuned cases.

The RMS position errors of the four filters with different mistuned clutter probabilities are shown in Figure 4. We only consider the runs without track loss. From the figure, we can see that the tracking accuracy of VB-SRF is slightly better than SRF when *η*<sup>0</sup> ≥ 0.5. However, the performance difference is not obvious. When *pc* = *η*<sup>0</sup> = 0.3, VB-SRF exhibits distinct superiority over SRF. It implies that VB-SRF is more robust than SRF. Meanwhile, MEFPDA-SCKF and PDA-SCKF show better tracking accuracy than SRF and VB-SRF under all the three mistuned cases. They are almost not affected by the mistuning. It shows MEFPDA-SCKF and PDA-SCKF are more accurate and robust than SRF and VB-SRF under this simple scenario.

Table 2 shows the computation time of the four filters with 100 Monte Carlo runs. It is clear that VB-SRF has the maximum time of computation, which is about 2 times of that of SRF. The computation time of MEFPDA-SCKF and PDA-SCKF are comparable and both smaller than SRF and VB-SRF.

On the whole, the VB-SRF algorithm outperforms the SRF in terms of track continuity, track accuracy and robustness especially in severely mismatched scenarios but with higher computation complexity. MEFPDA-SCKF and PDA-SCKF perform better than SRF and VB-SRF in all aspects. It illustrates that the PDA-SCKF-based strategy has superiority over the SRF-based strategy in handing the clutters in simple scenarios.


**Table 1.** The percentages of track losses of four filters in two scenarios.

 7LPHV 506SRVLWLRQHUURUVP 9%í65) 65) 3'\$í65&.) 0()3'\$í65&.) (**a**) *pc* = *η*<sup>0</sup> = 0.7 9%í65) 65) 3'\$í65&.) 0()3'\$í65&.)

(**b**) *pc* = *η*<sup>0</sup> = 0.5

**Figure 4.** *Cont.*

**Figure 4.** RMS target position errors with different mistuned clutter probabilities in Scenario 1.



### *5.2. Scenario 2*

For scenario 2, the aim is to track a single target from several drifting sonobuoys. A monitoring aircraft estimates the positions of the drifting sonobuoys by observing the direction of arrival of sensor transmissions. The sonobuoys track the position of the target by means of noisy bearings measurements.

The state is 12-dimensional:

$$\mathbf{x}\_k = [\mathbf{x}\_{0,k'} \dot{\mathbf{x}}\_{0,k'} \dot{\mathbf{y}}\_{0,k'} \dot{\mathbf{y}}\_{0,k'} \mathbf{x}\_{1,k'} \boldsymbol{y}\_{1,k'} \mathbf{x}\_{2,k'} \mathbf{y}\_{2,k'} \mathbf{x}\_{3,k'} \mathbf{y}\_{3,k'} \boldsymbol{u}\_{1,k'} \boldsymbol{u}\_{1,k}]^T \tag{48}$$

the first four components represent the (*x*, *y*) coordinates of the position and velocity of the target, the next six, the coordinates of the positions of the three sonobuoys, and the last two, those of the drift current effecting all three sonobuoys.

Six simultaneous measurements are made at each time step. Three of these are measurements of the bearing angles of the sonobuoys from the monitoring platform and they are uncluttered. Three are the bearing angles of the target from the sonobuoys, which are subject to clutter. The standard deviation of monitoring sensor bearing noise and sonobuoy sensor bearing noise are 0.8◦ and 16◦, respectively. The true probability of clutter is set as 0.667. In addition, the bearing of the clutter is uniformly distributed over [−*π*, *π*]. Other simulation parameters can be referred to [12]. 200 Monte Carlo runs are performed to evaluate the performance of the proposed filter. Figure 5 shows the behaviour of the estimates of target and sonobuoy positions provided by both the SRF and VB-SRF, for a typical simulation.

**Figure 5.** Typical tracks of target, drifting sonobuoy sensors, together with the estimated tracks.

The RMS target position errors of the four filters with correct clutter probability assumption are given in Figure 6. Compared with scenario 1, the differences between the four filters are more dramatic in scenario 2. This is probably because scenario 2 is more complex in which multiple sensors are used to observe the target and give abruptly changing and severely noise-corrupted bearing measurements of the target, shown in Figure 7. Thus, even a minor change in filtering strategy could result in large variations in performance. Meanwhile, seen from Figure 7, an abrupt change (almost from +180◦ to −180◦) occurs in the target bearing measurement from sonobuoy sensor 3 at *k* = 62 s, which leads to several track losses shown in Table 1 and much larger position errors of MEFPDA-SCKF and PDA-SCKF. Whereas, SRF and VB-SRF are less affected by the abrupt bearing variation since the value of the projected measurement *bk* = (*sinθk*, *cosθk*)*<sup>T</sup>* is invariant when there is a 360◦ change in bearing *θk*. They have comparable tracking accuracy. It is hard to decide which one is better.

**Figure 6.** RMS target position errors with correct clutter probability in Scenario 2.

**Figure 7.** The target measurements from three sonobuoy sensors in Scenario 2.

To compare the performance under adverse scenarios, we set mistuned clutter probabilities as: (1) *pc* = *η*<sup>0</sup> = 0.5; (2) *pc* = *η*<sup>0</sup> = 0.3. In case (1), there is no track loss in SRF and VB-SRF and the RMS position errors of the two filters are shown in Figure 8a. We can see that the RMS position errors of SRF are slightly increased compared with the case with no mistuning, while the RMS position errors of VB-SRF remain almost unchanged. In case (2), as shown in Table 1, 2.7% tracks are lost for SRF while no track is lost for VB-SRF. Meanwhile, as can be seen in Figure 8b, VB-SRF has much smaller RMS errors than the SRF. For MEFPDA-SCKF and PDA-SCKF, unlike with scenario 1, they have higher percentage of track losses and larger RMS position errors than VB-SRF in both mistuned cases. It shows that VB-SRF is superior to PDA-SCKF-based algorithms in more challenging scenarios. In addition, from Table 2, we can see that the computation time of VB-SRF is twice of the SRF and three times of MEFPDA-SCKF and PDA-SCKF. Overall, all these reveal that the proposed VB-SRF algorithm has significant performance superiority in severely mismatched and complex cases at the cost of a little higher computation complexity.

**Figure 8.** *Cont*.

**Figure 8.** RMS target position errors with mistuned clutter probability in Scenario 2.

### **6. Conclusions**

Bearings-only target tracking in the presence of clutter is a difficult problem because of the nonlinearity of the measurement model, the measurement origin uncertainty and the observability of the target. The Shifted Rayleigh filter (SRF) is shown to exhibit good performance for bearings-only target tracking in certain challenging scenarios through exploiting the essential structure of the nonlinearities in a new way. However, the clutter probability is assumed known and constant in SRF, which may not match with the truth especially in adverse scenarios. Therefore, to handle the bearings-only target tracking in clutters with uncertain clutter probability, a variational Bayesian-based adaptive shifted Rayleigh filter (VB-SRF) is proposed in this paper. By establishing a conjugate exponential model of the clutter probability and the data association indicator, the approximated posterior probability densities of the target state and the clutter parameters are iteratively calculated using the VB expectation and maximization steps. Finally, joint estimation of the target state and the clutter probability are achieved in the framework of SRF. The tracking performance of the proposed filter is compared with SRF, PDA-SCKF and MEFPDA-SCKF via two simulation examples. It shows that the proposed filter outperforms the other three filters in terms of track continuity and track accuracy with a little higher computation complexity in complex adverse scenarios. In addition, it also reveals that the proposed VB-SRF exhibits better robustness than the traditional SRF.

**Author Contributions:** J.H. conceived this paper, derived the method, and wrote the original draft ; Y.Y. made the investigation; and T.G. revised the paper and provided some valuable suggestions.

**Funding:** This research was funded by Aeronautical Science Foundation of China (No. 2016ZC53033), and Shaanxi Provincial Key Research and Development Programs (No. 2017ZDXM-GY-06, No. 2017GY-057), and Xi'an Science and Technology Planning Project–Scientific and Technological Innovation Guidance Project (No. 201805042YD20CG26 (8)).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**

The following abbreviations are used in this manuscript:



### **Appendix A. Derivation of** *<sup>f</sup>* **(***θk|***x***k***)**

Given the true state **x***k*, the "augmented" measurement **y***<sup>k</sup>* is assumed to be N (**m***k*, **Q***w*) variable, where **m***<sup>k</sup>* = **H***k***x***<sup>k</sup>* + **u***<sup>m</sup> <sup>k</sup>* . *θ<sup>k</sup>* is the bearing of **y***k*.

Let **y** *<sup>k</sup>* <sup>=</sup> **<sup>Q</sup>**−1/2 *<sup>w</sup>* **<sup>y</sup>***k*, then **<sup>y</sup>** *<sup>k</sup>* ∼ N (**Q**−1/2 *<sup>w</sup>* **<sup>m</sup>***k*,**I**). Its bearing has a density about its mean of the form *ασ*(*θ*). More precisely, let

$$\mathbf{g}\_k(\theta) = [(\sin \theta, \cos \theta)\mathbf{Q}\_w^{-1}(\sin \theta, \cos \theta)^T]^{-1/2} \tag{A1}$$

and define

$$h\_k(\theta) = \arctan\left(\frac{a\_{11}\sin\theta + a\_{12}\cos\theta}{a\_{21}\sin\theta + a\_{22}\cos\theta}\right) \tag{A2}$$

where [*aij*] = **Q**−1/2 *<sup>w</sup>* .

Actually, *gk*(*θ*) and *hk*(*θ*) are the reciprocal length and bearing of the transformed unit vector **<sup>Q</sup>**−1/2 *<sup>w</sup>* **<sup>y</sup>***<sup>k</sup>* ||**y***<sup>k</sup>* || . The bearing *hk*(*θk*) then has *<sup>α</sup>gk* (*θ<sup>m</sup> <sup>k</sup>* )(*hk*(*θk*) <sup>−</sup> *hk*(*θ<sup>m</sup> <sup>k</sup>* )) as its density, where *<sup>θ</sup><sup>m</sup> <sup>k</sup>* is the bearing of **m***k*. Inserting a Jacobian term, we can obtain the likelihood of measurement *θk*:

$$f(\theta\_k | \mathbf{x}\_k) = \frac{g\_k^2(\theta\_k)}{(\det \mathbf{Q}\_w)^{1/2}} \mathbf{a}\_{\S \boldsymbol{k}}(\theta\_k^m) \left(\boldsymbol{h}\_k(\theta\_k) - \boldsymbol{h}\_k(\theta\_k^m)\right) \tag{A3}$$

### **Appendix B. Derivation of** *<sup>f</sup>* **(***θk|***z1:***<sup>k</sup>−***1)**

Given the previous measurements **y**1:*k*−1, the "augmented" measurement **y***<sup>k</sup>* is assumed to be <sup>N</sup> (**y**<sup>ˆ</sup> *<sup>k</sup>*, **<sup>S</sup>***k*) variable, where ˆ**y***<sup>k</sup>* <sup>=</sup> **<sup>H</sup>***k***x**<sup>ˆ</sup> *<sup>k</sup>*|*k*−<sup>1</sup> <sup>+</sup> **<sup>u</sup>***<sup>m</sup> <sup>k</sup>* , **<sup>S</sup>***<sup>k</sup>* <sup>=</sup> **<sup>H</sup>***k***P***k*|*k*−1**H***<sup>T</sup> <sup>k</sup>* + **Q***w*.

Let **y** *<sup>k</sup>* <sup>=</sup> **<sup>S</sup>**−1/2 *<sup>k</sup>* **y***k*, then **y** *<sup>k</sup>* ∼ N (**S**−1/2 *<sup>k</sup>* **y**ˆ *<sup>k</sup>*,**I**). The bearing of **y** *<sup>k</sup>* has a density about its mean of the form *ασ*(*θ*).

Let

$$\mathbf{g}\_k'(\theta) = [(\sin \theta, \cos \theta)\mathbf{S}\_k^{-1}(\sin \theta, \cos \theta)^T]^{-1/2} \tag{A4}$$

$$h\_k'(\theta) = \arctan\left(\frac{s\_{11}\sin\theta + s\_{12}\cos\theta}{s\_{21}\sin\theta + s\_{22}\cos\theta}\right) \tag{A5}$$

where [*sij*] = **S**−1/2 *<sup>k</sup>* .

*g <sup>k</sup>*(*θ*) and *h <sup>k</sup>*(*θ*) are the reciprocal length and bearing of the transformed unit vector **<sup>S</sup>**−1/2 *<sup>k</sup>* **y***<sup>k</sup>* ||**y***<sup>k</sup>* || . The bearing *h <sup>k</sup>*(*θk*) then has *αg k* (ˆ *<sup>θ</sup><sup>k</sup>* )(*h <sup>k</sup>*(*θk*) − *h k*(ˆ *θk*)) as its density, where ˆ *θ<sup>k</sup>* is the bearing of **y**ˆ *<sup>k</sup>*. Inserting a Jacobian term, we can obtain the probability density function of measurement *θ<sup>k</sup>* given previous measurements **z**1:*k*−1:

$$f(\theta\_k | \mathbf{z}\_{1:k-1}) = \frac{g\_k'^2(\theta\_k)}{(det \mathbf{S}\_k)^{1/2}} a\_{\underline{\mathcal{S}}\_k'(\theta\_k)} (h\_k'(\theta\_k) - h\_k(\theta\_k)) \tag{A6}$$

### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **DOA Estimation and Self-calibration under Unknown Mutual Coupling**

### **Dong Qi \* , Min Tang, Shiwen Chen, Zhixin Liu and Yongjun Zhao**

National Digital Switching System Engineering and Technological Research Center (NDSC), Zhengzhou 86-450001, China; tangminmvp@126.com (M.T.); ndsccsw@126.com (S.C.); liuzhixin54@sina.com (Z.L.); zhaoyjzz@163.com (Y.Z.)

**\*** Correspondence: qidong1027@126.com; Tel.: +86-17513261827

Received: 22 January 2019; Accepted: 21 February 2019; Published: 25 February 2019

**Abstract:** In practical applications, the assumption of omnidirectional elements is not effective in general, which leads to the direction-dependent mutual coupling (MC). Under this condition, the performance of traditional calibration algorithms suffers. This paper proposes a new self-calibration method based on the time-frequency distributions (TFDs) in the presence of direction-dependent MC. Firstly, the time-frequency (TF) transformation is used to calculate the space-time-frequency distributions (STFDs) matrix of received signals. After that, the estimated steering vector and corresponding noise subspace are estimated by the steps of noise removing, single-source TF points extracting and clustering. Then according to the transformation relationship between the MC coefficients, steering vector and MC matrix, we deduce a set of linear equations. Finally, with two-step alternating iteration, the equations are solved by least square method in order to estimate DOA and MC coefficients. Simulations results show that the proposed algorithm can achieve direction-dependent MC self-calibration and outperforms the existing algorithms.

**Keywords:** DOA estimation; direction-dependent mutual coupling; time-frequency distribution; self-calibration

### **1. Introduction**

DOA estimation, as an important branch of array signal processing, is widely used in radar, sonar, radio astronomy and other fields [1]. In the past decades, many classical algorithms have been proposed which perform well in ideal situations, such as MUSIC, ESPRIT and other subspace-based algorithms. However, in practical applications, the performance of above algorithms suffers due to the effect of gain/phase uncertainties [2], sensor position perturbation [3] and mutual coupling (MC) [4], among which the mutual coupling caused by mutual excitation of array elements are common in engineering and makes the estimation accuracy deteriorate seriously.

A number of methods have been proposed for DOA estimation in the presence of mutual coupling, which can be classified into two types: active-calibration [5] and self-calibration [6]. The active-calibration method which makes use of the calibration sources whose DOAs are exactly known, can achieve high DOA accuracy with low computation complexity. However, the existence of the calibration sources increases the additional cost of the system, and the performance of the algorithm deteriorates rapidly in the presence of DOA errors of the calibration sources, which is inevitable in practice.

On the contrary, self-calibration method is preferable since it does not require any prior knowledge of source locations and accomplishes the DOA estimation and error calibration online. In [7], an iterative algorithm is proposed for the estimation of DOA and MC coefficients, however, the result will converge to the local optimum if the initial values deviate far from the real ones. For uniform circular arrays, a self-calibration method is proposed based on rank-reduction estimator by using

the complex symmetric Toeplitz property in [8]. The algorithm only needs one-dimensional search which lowers the computational complexity, but its parameter estimation is prone to be ambiguous. Moreover, reference [9,10] utilize alternating iteration and recursive estimation to solve this problem, respectively. Recently the methods which make use of instrumental sensors for array calibration has also been developed [11,12]. They exploit the fact that only part of the new array has mutual coupling or other errors after adding instrumental sensors into the original array. But in practice, it is impossible to obtain the ideal instrumental sensors.

The above algorithms are only adopted to direction-independent mutual coupling which is modelled with a single matrix. As it is established under the assumption of the omnidirectional antenna array, the model becomes ineffective when the array elements are not omnidirectional antennas. However, in practical engineering, due to the limitations of the manufacture and the working environment of the antenna, the array elements have directional beam pattern in general [13]. As a result, the mutual coupling is direction-dependent, leading to performance degradation of the existing algorithms. Thus, it is of great practical significance to study the DOA estimation in the presence of direction-dependent mutual coupling. Few papers are proposed to solve this problem. Based on rank-rare theory, a method of 2D-DOA estimation for direction-dependent MC is proposed in [14]. However, in order to estimate the direction-dependent MC coefficients, the algorithm adopts the idea of receiving mutual-impedance method proposed in [15], which is an off-line measurement algorithm. Therefore, the algorithm fails in the time-varying systems. In [16] Ahmet proposed a method to calibrate the direction-dependent mutual coupling and estimate the DOAs. The algorithm divides the angle search range into several sectors by means of angle sector. By comparing the spectral peaks in each sector, the angle interval of initial angle and the corresponding MC coefficients are estimated, and then self-calibration is completed by iterations. However, this method is easy to fail when the angular spacing between the DOAs of the incident signals is small or the initial value deviation is large.

Motivated by these facts, in this paper a new algorithm for estimation of DOAs and MC coefficients is proposed based on the idea of time-frequency distributions (TFDs) which has been widely applied in blind source separation [17,18]. Firstly, the space-time-frequency distributions (STFDs) matrix of the received signal is solved. Then, the single-source time-frequency (TF) points are extracted by denoising and removing the cross-terms at each TF point, and the optimal STFDs matrix of each signal is estimated by clustering the single-source TF points. Finally, the STFDs matrix is decomposed into eigenvectors and noise subspaces of each signal, and an alternating iteration method based on least squares is proposed for estimation of DOA and MC coefficients. The simulation results show that the algorithm can provide satisfactory performance in case of direction-dependent MC. The main contributions are as follows:


The rest of the paper is organized as follows: Section 2.1 is devoted to the problem formulation. Then an approach based on the TFDs is proposed to obtain the steering vector and noise subspace of each signal by separating the mixed signals, and DOAs and MC coefficients are estimated by a two-step alternate iteration in Section 2.2. Next the algorithm analysis is provided in Section 2.3. Simulations are conducted to illustrate the effectiveness of the proposed methods in Section 3 and conclusions are finally drawn in Section 4.

### **2. Models and Methods**

### *2.1. Array Signal Model*

Considering *<sup>K</sup>* far-field narrow-band signals *<sup>s</sup>k*(*t*)(*<sup>k</sup>* <sup>=</sup> 1, 2, ··· , *<sup>K</sup>*) impinging on the array which is composed of *M* elements, and the directions of arrival are {*θ*1, *θ*2, ··· , *θK*}, respectively. Then the received signals at the *t*-th sample can be expressed as:

$$\mathbf{X}(t) = \mathbf{A}(\theta)\mathbf{s}(t) + \mathbf{N}(t)t = 1, 2, \cdots, T \tag{1}$$

where *<sup>X</sup>*(*t*)=[*x*1(*t*), *<sup>x</sup>*2(*t*),..., *xM*(*t*)]*<sup>T</sup>* are the array outputs. *<sup>N</sup>*(*t*)=[*n*1(*t*), *<sup>n</sup>*2(*t*),..., *nM*(*t*)]*<sup>T</sup>* denotes zero-mean additive white Gaussian noise. *<sup>A</sup>*(*θ*) = [*a*(*θ*1), *<sup>a</sup>*(*θ*2), ··· , *<sup>a</sup>*(*θK*)] is the ideal manifold matrix, *a*(*θi*) is the steering vector of the *i*-th signal.

In the presence of direction-dependent MC, the real steering vector is written as:

$$b(\theta) = \mathcal{C}(\theta)a(\theta) \tag{2}$$

where *<sup>C</sup>*(*θ*) <sup>∈</sup> <sup>C</sup>*M*×*<sup>M</sup>* is the MC matrix.

For a uniform linear array or a uniform circular array model, the MC matrix can be expressed as a band complex symmetric Toeplitz matrix or a three-band complex cyclic matrix, respectively. The following transformation form is used to represent the MC matrix uniformly:

$$\mathcal{C}(\theta\_k) = \text{CMC}(\mathfrak{c}(\theta\_k)) \tag{3}$$

where *CMC*(·) is the operation of constructing the MC matrix using MC coefficients. *c*(*θk*)=[*c*1*k*, *c*2*k*, ··· , *cpk*, **0***M*−*p*] *<sup>T</sup>* <sup>∈</sup> <sup>C</sup>*M*×<sup>1</sup> are the MC coefficients, **<sup>0</sup>***M–p* is a 1 <sup>×</sup> (*M–p*) zero vector and *<sup>p</sup>* is the degree of freedom of MC, which equals the number of non-zero elements of *c*(*θk*), so in the calculation process, only the first *p* elements are solved. Thus, the simplified vector [*c*1*k*, *c*2*k*, ··· , *cpk*] *T* is used in the formula derivation and simulation conditions. Taking the uniform linear array as an example, the MC matrix of *θ<sup>k</sup>* can be expressed as:

$$\mathbf{C}(\theta\_k) = \begin{bmatrix} c\_{1k} & \cdots & c\_{pk} \\ \vdots & \ddots & & \ddots \\ c\_{pk} & & c\_{1k} & & c\_{pk} \\ & & \ddots & & \ddots & \vdots \\ & & & c\_{pk} & \cdots & c\_{1k} \end{bmatrix}\_{M \times M} \tag{4}$$

Then the array receiving model is expressed as

$$\mathbf{X}(t) = \mathbf{B}(\theta)\mathbf{s}(t) + \mathbf{N}(t) \tag{5}$$

where *<sup>B</sup>*(*θ*)=[*b*(*θ*1), *<sup>b</sup>*(*θ*2), ··· , *<sup>b</sup>*(*θK*)] = [*C*(*θ*1)*a*(*θ*1),*C*(*θ*2)*a*(*θ*2), ··· ,*C*(*θK*)*a*(*θK*)] is the real array manifold matrix containing direction-dependent MC.

### *2.2. Method of DOA and MC Coefficients Estimation*

In this section, the TFDs are introduced to obtain the true steering vector and noise-subspace of each signal, with which DOA and MC coefficients are estimated. The algorithm flow is shown in Figure 1.

**Figure 1.** The flow diagram of the proposed algorithm.

The first step is obtaining the STFD matrix by quadric time-frequency transform of *X*(*t*). Secondly the single-source self-term TF points are extracted by denoising and removing the cross-term and common-term generated by the signals. After Step 2, the optional STFD matrix of each signal is obtained by clustering TF points of the same signal. Then the steering vector and noise-subspace of each signal are estimated by the eigen-decomposition of STFD matrix. Finally, the two-step iteration is performed to estimate the DOAs and corresponding MC coefficients. The following is the method of DOA and MC coefficient estimation. Firstly, the basic concept of the TFDs and the steps of steering vector estimation are introduced.

### 2.2.1. Steering Vector Estimation Based on TFDs

For a single signal *x*(*t*), the TFDs in discrete time form can be expressed as:

$$\rho\_{\text{xx}}(t,f) = \sum\_{k=-\infty}^{+\infty} \sum\_{l=-\infty}^{+\infty} \mathbf{x}(t+k+l)\mathbf{x}^\*(t+k-l)\boldsymbol{\varrho}(k,l)e^{-j4\pi f l} \tag{6}$$

where *ϕ*(*k*, *l*) is the kernel function, (·) <sup>∗</sup> denotes the conjugate operator, *ρxx*(*t*, *f*) is the self-term TF point of *x*(*t*).

For two signals *x*1(*t*) and *x*2(*t*), the cross-time-frequency distributions in discrete time form can be expressed as:

$$\rho\_{X1\overline{x}2}(t,f) = \sum\_{k=-\infty}^{+\infty} \sum\_{l=-\infty}^{+\infty} x\_1(t+k+l)x\_2^\*(t+k-l)\,\rho(k,l)e^{-j4\pi fl} \tag{7}$$

*ρx*1*x*<sup>2</sup> (*t*, *f*) is the cross-term TF point of *x*1(*t*) and *x*2(*t*).

Therefore, the STFDs matrix of *X*(*t*) is defined as

$$D\_{\mathbf{X}\mathbf{X}}(t,f) = \sum\_{k=-\infty}^{+\infty} \sum\_{l=-\infty}^{+\infty} \mathbf{X}(t+k+l)\mathbf{X}^{H}(t+k-l)\,\boldsymbol{\varrho}(k,l)e^{-j4\pi f l} \tag{8}$$

According to the definition in paper [17], TF points in time-frequency domains can be divided into three classes including self-term TF points, single-source self-term TF points and cross-term TF points which are expressed as (*ta*, *fa*),(*tas*, *fas*) and (*tc*, *fc*), respectively.

The energy of self-term TF points is generated by one or more sources whose cross-term energy is approximately zero. Single-source self-term TF points are generated by only the self-term of single signal, and the cross-term TF points are generated by the cross-term of the signals whose self-term energy is nearly zero.

Under the noiseless condition, substituting Equation (5) into Equation (8), we obtain:

$$\begin{array}{rcl} \mathbf{D}\_{\mathbf{X}\mathbf{X}}(t,f) &=& \sum\_{k=-\infty}^{+\infty} \sum\_{l=-\infty}^{+\infty} \mathbf{B} \mathbf{s}(t+k+l) \mathbf{s}^{H}(t+k-l) f(k,l) e^{-j4\pi f l} \mathbf{B}^{H} \\ &= \mathbf{B} \left( \sum\_{k=-\infty}^{+\infty} \sum\_{l=-\infty}^{+\infty} \mathbf{S}(t+k+l) \mathbf{S}^{H}(t+k-l) f(k,l) e^{-j4\pi f l} \right) \mathbf{B}^{H} \\ &= \mathbf{B} \mathbf{D}\_{\mathbf{s}\mathbf{s}}(t,f) \mathbf{B}^{H} \end{array} \tag{9}$$

where *DXX* is the STFDs matrix of array outputs. *Dss* is the STFDs matrix of the incident signals whose principal diagonal elements are generated by the self-term of incident signals and the non-principal diagonal elements are corresponding to the cross-terms between the signals.

*Sensors* **2019**, *19*, 978

For the self-term TF points, since the signal energy is mainly generated by the self-term of sources, the non-diagonal element of *Dss* is approximately zero, so it can be expressed as a diagonal matrix, so the STFDs matrix can be expressed as:

$$\mathbf{D} \mathbf{X} \mathbf{x} (t\_{\mathbf{d}\prime} f\_{\mathbf{d}}) = \mathbf{B} \begin{bmatrix} \rho\_{s\_1 s\_1} (t\_{\mathbf{d}\prime} f\_{\mathbf{d}}) & 0 & \dots & 0 \\ 0 & \rho\_{s\_2 s\_2} (t\_{\mathbf{d}\prime} f\_{\mathbf{d}}) & & \vdots \\ \vdots & & \ddots & 0 \\ 0 & \dots & 0 & \rho\_{\delta^{\prime} \mathbf{x}\_{\mathbf{k}}} (t\_{\mathbf{d}\prime} f\_{\mathbf{d}}) \end{bmatrix} \mathbf{B}^H \tag{10}$$

On this basis, when the energy at the TF points is generated only by the single source *si* (*i* = 1, 2 ··· *<sup>K</sup>*), *DXX*(*tas*, *fas*) can be expressed as:

$$D\_{\mathbf{X}\mathbf{X}}(t\_{\text{as}\star}f\_{\text{as}}) := \rho\_{s\_i s\_i}(t\_{\text{as}\star}f\_{\text{as}}) \mathfrak{a}\_i \mathfrak{a}\_i^H \tag{11}$$

It is deduced from Equation (11) that the steering vector of each signal can be obtained by eigen-decomposition of the STFDs matrix at the TF points of single-source self-term. The procedure of extracting TF points from single source term is as follows.

*Step 1*: Remove the noise: In order to reduce the computation burden and improve the estimation accuracy, it is necessary to set an appropriate threshold to remove the TF points with low energy which may be generated by noise. For each time slice (*tp*, *f*), apply Equation (12) for all frequency *fq* points in this slice, and then the TF points with large energy remain:

$$\frac{||D\_{\text{XX}}(t\_{p\prime}f\_q)||}{\max\_f \{||D\text{xx}(t\_{p\prime}f)||\}} > \varepsilon\_1 \tag{12}$$

where *ε*<sup>1</sup> is a small positive real number and • represents the operator of 2-norm.

*Step 2*: Extract the TF points of self-term: After removing the noise points, the retaining TF points mainly include self-term TF points and cross-term TF points. At the self-term TF points, the STFDs matrix is approximately diagonal, and the values of the principal diagonal elements are much larger than those of the other elements, so the STFDs matrix at self-term TF points yields: *trace*{*DXX*(*t*, *<sup>f</sup>*)}

$$\frac{\text{trace}\{\mathbf{D}\_{\text{XX}}(t,f)\}}{\|\mathbf{D}\_{\text{XX}}(t,f)\|} > \varepsilon\_2 \tag{13}$$

where *ε*<sup>2</sup> is a positive real number close to but less than 1.

*Step 3*: Extract the TF points of single-source self-term: For signals which overlap in time-frequency domain, the self-term TF points may be composed of multiple signals. Therefore, it is necessary to extract the single-source self-term TF points from the self-term TF points, which can be accomplished by Equation (14):

$$\left| \frac{\lambda\_{\text{max}} \{ \mathbf{D} \mathbf{x} (t, f) \}}{\text{trace} \{ \mathbf{D}\_{\text{XX}} (t, f) \}} - 1 \right| \le \varepsilon\_3 \tag{14}$$

where *<sup>ε</sup>*<sup>3</sup> is a small positive real threshold and *<sup>λ</sup>*max{*DXX*(*t*, *<sup>f</sup>*)} is the largest eigenvalue of *DXX*(*t*, *f*).

With the three steps above, the single-source self-term TF points are obtained. Then the steering vectors of each signal as well as the noise subspace could be estimated by the eigen-decomposition of their STFDs matrix. However, the above derivation is completed without considering the noise. In the presence of noise, the steering vectors estimated by only a few TF points are biased. Therefore, it is necessary to obtain more accurate information by clustering the multiple single-source self-term TF points of the same signal. A time-frequency clustering method is provided in Step 4.

*Step 4*: Cluster the TF points of the same signal: The steering vector of each signal can be estimated as the principal eigenvector of STFDs matrix at each TF point. Regarding the steering vector which contains the DOA information as a feature, all self-term TF points can be classified into *Q*(*Q*≥*K*) categories by the classification algorithm. That is to say, if the following conditions are satisfied, TF points of (*t*1, *f*1),(*t*2, *f*2) belong to the same category:

$$d(a(t\_1, f\_1), a(t\_2, f\_2)) < \varepsilon\_4 \tag{15}$$

where *d*(*x,y*) is the Euclidean distance between *x* and *y*, and *ε*<sup>4</sup> is a small positive threshold.

After clustering, we extract the first *K* categories which contain the most TF points, and obtain the TFDs of the *K* signals, respectively. Then the STFDs matrices at those TF points belonging to the first *K* categories are summed and averaged. Finally, the eigen-decomposition of the average STFDs matrix is performed to estimate the steering vector of each signal and the corresponding noise subspace, which are denoted as &*b*(*θk*) and *<sup>E</sup>n*(*θk*), respectively, *<sup>k</sup>* <sup>=</sup> 1, 2, ··· , *<sup>K</sup>*.

The above steps of removing the noise, extracting self-term TF points, and clustering TF points of same signal aimed to essentially select the TF points that satisfy the specific conditions, and then the matrix at the TF point is processed.

### 2.2.2. DOA and MC Coefficients Estimation

For a uniform linear array or uniform circular array model, the transformation relationship between the MC coefficients vector and MC matrix can be expressed by Equation (16):

$$\mathfrak{b}(\theta\_k) = \mathfrak{C}(\theta\_k)\mathfrak{a}(\theta\_k) = T(\theta\_k)\mathfrak{c}(\theta\_k)k = 1, 2, \dots, \, ^\prime \text{,} K \tag{16}$$

where *T*(*θk*) is transformation matrix which contains the direction information.

For a uniform linear array, the transformation matrix can be expressed as:

$$T(\theta\_k) = \mathbf{Q}\_1(\theta\_k) + \mathbf{Q}\_2(\theta\_k) \tag{17}$$

where:

$$\begin{aligned} [Q\_1(\theta\_k)]\_{ij} &= \begin{cases} [a(\theta\_k)]\_{i+j-1} & i+j \le M+1\\ 0 & i+j > M+1 \end{cases} \\\ [Q\_2(\theta\_k)]\_{ij} &= \begin{cases} [a(\theta\_k)]\_{i-j+1} & i \ge j \ge 2\\ 0 & \text{otherwise} \end{cases} \end{aligned} \tag{18}$$

Similarly, for a uniform circular array the transformation matrix can be expressed as:

$$T(\theta\_k) = \mathbf{Q}\_1(\theta\_k) + \mathbf{Q}\_2(\theta\_k) + \mathbf{Q}\_3(\theta\_k) + \mathbf{Q}\_4(\theta\_k) \tag{19}$$

where *<sup>Q</sup>*1(*θk*), *<sup>Q</sup>*2(*θk*), *<sup>Q</sup>*3(*θk*), *<sup>Q</sup>*4(*θk*) yield:

$$\begin{aligned} [Q\_1(\theta\_k)]\_{ij} &= \begin{cases} [a(\theta\_k)]\_{i+j-1} & i+j \le M+1\\ 0 & i+j > M+1 \end{cases} \\ [Q\_2(\theta\_k)]\_{ij} &= \begin{cases} [a(\theta\_k)]\_{i-j+1} & i \ge j \ge 2\\ 0 & otherwise \end{cases} \\ [Q\_3(\theta\_k)]\_{ij} &= \begin{cases} [a(\theta\_i)]\_{M+1+i-j} & i < j \le l \\ 0 & otherwise \end{cases} \\ [Q\_4(\theta\_k)]\_{ij} &= \begin{cases} [a(\theta\_k)]\_{i+j-M-1} & 2 \le i \le l, i+j \ge M+2\\ 0 & otherwise \end{cases} \end{aligned} \tag{20}$$

As there is a multiple relation between the estimated steering vector and the actual steering vector, we deduce that:

$$\mathfrak{b}(\theta\_k) = \rho\_k \check{\mathfrak{b}}(\theta\_k) = \mathfrak{C}(\theta\_k)\mathfrak{a}(\theta\_k) = \mathfrak{T}(\theta\_k)\mathfrak{c}(\theta\_k) \tag{21}$$

where *ρ<sup>k</sup>* is the multiplier. And the MC coefficients *c*(*θk*) can be calculated by the least squares method:

$$\widetilde{\mathfrak{C}}(\theta\_k) = \rho\_k \left( T^H(\theta\_k) T(\theta\_k) \right)^{-1} T^H(\theta\_k) \widetilde{\mathfrak{b}}(\theta\_k) \tag{22}$$

According to Equation (22), we find that there also exists a corresponding coefficient relationship between the estimated MC coefficients &*c*(*θk*) and the true MC coefficients *<sup>c</sup>*(*θk*). Because the first element of &*c*(*θk*) should be 1, the actual MC coefficients can be obtained by normalizing the first element of &*c*(*θk*) with Equation (23):

$$
\hat{\mathfrak{c}}(\theta\_k) = \frac{\breve{\mathfrak{c}}(\theta\_k)}{[\breve{\mathfrak{c}}(\theta\_k)]\_1} \tag{23}
$$

where [&*c*(*θk*)] <sup>1</sup> represents the first element of &*c*(*θk*).

Then DOA of the *k*-th signal is estimated according to Equation (24) (root MUSIC algorithm could be applied to solve the functions) with the new MC matrix *C*ˆ(*θk*)= CMC(*c*ˆ(*θk*)):

$$\min\_{\theta\_k} \left\| E\_n^H(\theta\_k) \mathcal{C}(\theta\_k) \mathcal{a}(\theta\_k) \right\|^2 k = 1, 2, \cdots, \le K \tag{24}$$

However, as the direction *θ<sup>k</sup>* and the MC coefficients are unknown, the function cannot be solved directly. Based on the above steps, a new method for DOA estimation and MC coefficients is proposed based on the two-step alternating iteration.

Firstly, initialize the MC coefficients and compute the DOAs ˆ *θ*1, ˆ *<sup>θ</sup>*2, ··· , <sup>ˆ</sup> *θ<sup>K</sup>* through the conventional algorithms. Then construct the transformation matrix *T*(ˆ *θ*1), *T*(ˆ *<sup>θ</sup>*2), ··· , *<sup>T</sup>*(<sup>ˆ</sup> *θK*) according to ˆ *θ*1, ˆ *<sup>θ</sup>*2, ··· , <sup>ˆ</sup> *θK*. On this basis, the MC coefficients is estimated by Equations (22) and (23) with which the new MC matrix is constructed. Finally the DOA of each signal is obtained using Equation (24), and transformation matrix with the new angle information is constructed. In this way, iterations are carried out alternately until the estimation bias is less than the threshold or the iterations have been repeated for certain times. The proposed algorithm is summarized in Table 1.

**Table 1.** Algorithm for estimation of DOA and direction-dependent MC based on TFDs.

Step 7. Construct the MC matrix using the new MC coefficients and estimate the DOAs using (24).

Step 8. Repeat Step 5 to Step 7 until the estimation errors is less than the threshold or the iterations have been repeated for certain times.

### *2.3. Algorithmic Analysis*

1) For the selection of the TFD kernel function, it is known that different TFD kernel functions correspond to different TFD transforms, which can be divided into linear time-frequency transforms and quadratic time-frequency transforms. The Wigner-Ville distribution (WVD) is one of the quadratic TFDs that has better performance in TF focusing and resolution. These two elements are significant in the extraction and clustering of TF points. However, due to the interaction between different signals, cross-terms are generated, which cause false time-frequency information and degrade the estimation accuracy. Therefore, the smoothed pseudo-Wigner-Ville distribution (SPWVD) is utilized to suppress the cross-terms between signals by windowing

Step 1. Collect N snapshots and calculate the STFDs of received signal with (9).

Step 2. Extract the single source TF points based on (12)–(14).

Step 3. Cluster the single source TF points into *Q*(*Q*≥*K*) categories using (15).

Step 4. Estimate &*b*(*θk*) and *<sup>E</sup>n*(*θk*) of each signal with the *<sup>K</sup>* largest classes.

Step 5. Construct the transformation matrix with &*b*(*θk*).

Step 6. Estimate the MC coefficients based on Equations (22) and (23).

method in this paper. Similarly, the short-time Fourier transform (STFT) can also provide the accurate time-frequency distribution to solve the problem and we should choose the kernel function flexibly for the different conditions.


### **3. Numerical Simulation**

This section illustrates the effectiveness of the algorithm through simulations with other methods as comparison. Without loss of generality, a uniform linear array composed of 7 elements is selected in this paper, and *p* = 3. Assuming that there are three far-field narrowband LFM (linear frequency modulation) signals in the space, and the duration of three signals is the same. The sampling frequency of system is 50 MHz. The parameters of signals are listed in Table 2.


**Table 2.** Parameters of Linear Frequency Modulation Signals.

### *3.1. Calibration Results of Proposed Algorithm*

In this section, we provide the results of each step shown in Section 3.1 to verify the effectiveness of the algorithm. Figure 2a shows a three-dimensional time-frequency diagram of the received signal. From which, we find that cross-terms exist in the time-frequency domain of the three signals. Although the smooth pseudo-Wigner distribution is used to remove most of the cross-terms, there still remain some cross-terms in the overlapping part. Figure 2b is a two-dimensional time-frequency distribution which has been binarily processed. It can be seen that except for the time-frequency points of the signal itself, there are also noise points and the time-frequency points of cross-terms between different signals.

**Figure 2.** Time-Frequency Distribution of received signals. (**a**) Time-Frequency Distribution (3-dimensional); (**b**) Time-Frequency Distribution (2-dimensional).

Figure 3 indicates the processing results introduced in Section 3.1. Figure 3a shows the result of removing the noise points, in which we find that nearly all the noise points can be filtered, while only the TF points with high energy remain. Figure 3b shows the result of the self-term TF points extraction. It is seen that the cross-term TF points are removed but the self-term TF points in the overlapping area are removed as well. Figure 3c shows the result of the single-source self-term TF points extraction. We find that there are only single-source TF points remaining while the overlapping self-term TF points of multiple signals are removed. Figure 3d–f are the TF distribution diagrams of Signals 1–3 after clustering, respectively. It is found that the TF points of each signal can be obtained by clustering with spatial information as feature. And they are matched with the modulation frequency of the real signals.

**Figure 3.** Processing steps of TFDs. (**a**)The TFDs after removing noise points; (**b**) The TFDs of self-term points; (**c**) The TFDs of single-source self-term points; (**d**)The TFDs of signal 1 after clustering; (**e**) The TFDs of signal 2 after clustering; (**f**) The TFDs of signal 3 after clustering.

Figure 4a–c show the variation of the estimated DOAs versus the iteration times of three signals respectively. We find that after a few times of iteration, the estimated DOAs gradually approach the real value and finally stabilize. Figures 5 and 6 indicate variation of the estimated DOAs and MC coefficients versus the iteration times respectively. It can be seen from the figures that the algorithm has high estimation accuracy in DOA and MC coefficients under the condition of *SNR* = 10 *dB* after a few times of iteration.

**Figure 4.** Variation of the estimated DOAs versus the iterations of three signals. (**a**) Variation of the estimated DOAs versus the iterations of signal 1; (**b**) Variation of the estimated DOAs versus the iterations of signal 2;(**c**) Variation of the estimated DOAs versus the iterations of signal 3.

**Figure 5.** Variation of the estimated DOAs versus the iteration times.

**Figure 6.** Variation of the estimated MC coefficients versus the iteration times.

*3.2. RMSE Comparison versus Input Signal-Noise-Ratio (SNR)*

In this section, performance of the proposed algorithm is investigated with different input SNR from −10 dB to 14 dB, which is compared with PEDDMC algorithm [16], conventional MUSIC and the Cramer-Rao lower bound (CRLB) with unknown mutual coupling [7]. The RMSE are obtained through 200 Monte-Carlo simulations. The calculation formula of RMSE is as follows:

$$RMSE = \sqrt{\sum\_{n=1}^{N\_s} \sum\_{i=1}^{K} \left(\theta\_i - \hat{\theta}\_{i,n}\right)^2 / (K\text{N}\_s)}\tag{25}$$

As shown in the Figure 7, when the SNR is low, estimation errors of the three algorithms are large. With the increase of SNR, performance of the proposed algorithm and PEDDMC algorithm are improved. However, the proposed algorithm has lower RMSE because the TFDs algorithm has better anti-noise performance, and in addition the DOA estimation and error calibration for each signal are carried out separately. As a result, the proposed method has better performance and its RMSE follows the CRLB closely. As for the MUSIC algorithm, it is a classical super-resolution algorithm with superior estimation performance and robustness under ideal conditions, but it fails in the presence of mutual coupling conditions, since it has no ability to achieve calibration. As a result its performance does not improve as SNR increases.

**Figure 7.** Comparison of DOA estimation performance versus SNR.

### *3.3. RMSE Comparison versus Input Snapshots*

This section compares the performance of the algorithms with different number of snapshots from 100 to 1000. We fix the SNR at 15 dB and calculate the output RMSE through 200 Monte Carlo simulations. From Figure 8, it is known that the estimation performance of proposed algorithm is poor under the condition of small snapshots. This is because small snapshots will lead to the TF resolution deterioration of time-frequency distribution, which affects the extracting and clustering of TF points from single-source self-term. With the increase of snapshots, the RMSE of the proposed algorithm gradually becomes lower than that of the PEDDMC algorithm and becomes very close to the CRLB.

**Figure 8.** Comparison of DOA estimation performance versus snapshots.

### **4. Conclusions**

Considering the problem that the performance of traditional calibration algorithms degrades in the presence of direction-dependent mutual coupling, this paper introduces the time-frequency analysis method into array signal processing and proposes a new self-calibration algorithm based on alternating iteration. The simulation results show that the proposed algorithm is effective in the presence of direction-dependent mutual coupling and outperforms the existing algorithms. Though the computational complexity of the proposed algorithm is higher because it requires more snapshots to ensure the estimation accuracy of time-frequency transform, compared with the other existing algorithms, the proposed algorithm can perform DOA estimation under multipath or underdetermined conditions.

### **5. Patents**

The authors would like to thank the editor and anonymous reviewers for their careful reading and constructive comments which provide an important guidance for our paper writing and research work. This work was supported by the National Natural Science Foundation of China under Grant 61703433.

**Author Contributions:** All authors contributed extensively to the study presented in this manuscript. D.Q. designed the main idea, methods and experiments, interpreted the results and wrote the original paper. M.T. and Z.L supervised the main idea, edited the manuscript and provided many valuable suggestions on visualization. S.C. and Y.Z. offered help with the funding.

**Funding:** This research was funded by National Natural Science Foundation of China, grant number 61703433.

**Acknowledgments:** The authors would like to thank anonymous reviewers for their valuable comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Time Di**ff**erence of Arrival (TDoA) Localization Combining Weighted Least Squares and Firefly Algorithm**

### **Peng Wu, Shaojing Su, Zhen Zuo \*, Xiaojun Guo, Bei Sun and Xudong Wen**

College of Intelligence Science and Technology, National University of Defense Technology, Changsha 410073, China; pengwu9510@163.com (P.W.); susj-5@163.com (S.S.); jeanakin@nudt.edu.cn (X.G.); beys1990@163.com (B.S.); wenxudong13@163.com (X.W.)

**\*** Correspondence: z.zuo@nudt.edu.cn

Received: 11 April 2019; Accepted: 2 June 2019; Published: 4 June 2019

**Abstract:** Time difference of arrival (TDoA) based on a group of sensor nodes with known locations has been widely used to locate targets. Two-step weighted least squares (TSWLS), constrained weighted least squares (CWLS), and Newton–Raphson (NR) iteration are commonly used passive location methods, among which the initial position is needed and the complexity is high. This paper proposes a hybrid firefly algorithm (hybrid-FA) method, combining the weighted least squares (WLS) algorithm and FA, which can reduce computation as well as achieve high accuracy. The WLS algorithm is performed first, the result of which is used to restrict the search region for the FA method. Simulations showed that the hybrid-FA method required far fewer iterations than the FA method alone to achieve the same accuracy. Additionally, two experiments were conducted to compare the results of hybrid-FA with other methods. The findings indicated that the root-mean-square error (RMSE) and mean distance error of the hybrid-FA method were lower than that of the NR, TSWLS, and genetic algorithm (GA). On the whole, the hybrid-FA outperformed the NR, TSWLS, and GA for TDoA measurement.

**Keywords:** TDoA; weighted least squares; firefly algorithm; hybrid-FA

### **1. Introduction**

Target localization based on a group of sensor nodes whose positions are known has been extensively studied in research on signal processing [1–3]. It has been applied widely in military and civil fields, including sensor networks [4], wireless communication [2], radar [5], navigation, and so forth [6–8]. Commonly adopted positioning methods include the signal's time of arrival (ToA) [9], time difference of arrival (TDoA), frequency difference of arrival (FDoA) [10–13], or doppler shift [14]. Compared with FDoA, the TDoA and ToA methods can achieve higher positioning accuracy and require only one channel for each sensor node to perform the measurement, which can minimize the load requirement for a single-sensor node.

For a passive location system based on TDoA, once the measured data are obtained, the range difference between the target and two different sensor nodes can be calculated. In this connection, a set of hyperbolic equations or hyperboloids can be obtained and the solution of the equations is the coordinate of the target. Generally, the solving algorithms commonly adopted include iterative, analytical, and search methods.

The procedure for solving equations from the TDoA method is complex and difficult because the equations are nonlinear, and many studies have been carried out on how to solve this issue. The main idea of the Taylor series method is to expand the first Taylor series of the nonlinear positioning equation at the initial estimation of the target position and then solve the equations by iteration [15]. The advantage of this method is that it can fuse multiple observation data. Yang et al. transformed the equation into a constrained weighted least squares (CWLS) estimation problem by introducing auxiliary variables, and then the Newton iteration method was adopted to solve the problem [16]. In [17], nonconvex TDoA localization was transformed into a convex semidefinite programming (SDP) problem, and the approximate result was taken as the initial value for the Newton iteration method. All of these methods are iterative. Compared with iterative methods, the closed-form method does not need the initial estimation of the target's location and iterative solving is not necessary either. For example, Chan and Ho [18] transformed nonlinear equations into pseudolinear equations by introducing auxiliary variables. Then, the equations were solved by two-step weighted least squares (TSWLS). One downside of this algorithm is that the result is substantially different than the actual position when the signal-to-noise ratio (SNR) is low. Considering this problem, the constrained total least squares (CTLS) method was proposed in [19–21]. While it is not a closed-form method and Newton iterations are needed, the complexity of the CTLS method is much higher than that of the TSWLS method. The approximate maximum likelihood (AML) method [22] was proposed, which can obtain a linear equation from the maximum likelihood function and then the target location can be calculated. The AML method has better positioning performance than the TSWLS method.

In addition to the traditional TSWLS, iteration methods, and so on, many scholars have investigated new methods to enhance positioning accuracy. Two new shrinking-circle methods were proposed (SC-1 and SC-2) to solve a TDoA-based localization problem in a 2-D space [23]. Additionally, a weighted least squares (WLS) algorithm with the cone tangent plane constraint for hyperbolic positioning was proposed, which added the distance between the target and the reference sensor as a new dimension [24]. The theoretical bias of maximum likelihood estimation (MLE) is derived when sensor location errors and positioning measurement noise both exist [25]. Using a rough estimated result by MLE to subtract the theoretical bias can deliver a more accurate source location estimation. Apart from this, research based on certain typical algorithms has been carried out to extract and calculate the TDoA of ultrahigh frequency (UHF) signals [26]. The AML algorithm was proposed for determining a moving target's position and velocity by utilizing TDoA and FDoA measurements [27].

It is also efficient to use a search algorithm to calculate the position of a target. A hybrid genetic algorithm (GA) was proposed to enhance solution accuracy [18]. Nature-inspired algorithms are powerful algorithms for optimization. The firefly algorithm (FA) is one such nature-inspired algorithm, which was proposed in 2008. Using the FA for multimodal optimization applications with high efficiency has been proposed [28,29].

In general, among the methods for solving TDoA equations, analytical and iterative methods both have limitations. Research on algorithms that are robust and have low computational complexity is still worthy of study. Search algorithms for TDoA measurements can provide accurate results, although the efficiency will inevitably decrease when there are many estimated parameters [29]. Therefore, it is necessary to develop a highly efficient search algorithm for TDoA.

This paper is organized as follows: Section 2 introduces the basic model of TDoA measurement. Section 3 formulates the basic principle of WLS. Section 4 provides the main steps of the FA. Section 5 details the hybrid-FA methods proposed in this paper. Section 6 presents the results of simulations and experiments to support the theoretical analysis.

### **2. Problem Description**

In this section, 2-D target localization based on TDoA measurement is presented in the line-of-sight environment. Assume that there are *N* (*N* ≥ 3) sensor nodes, which can also be called basic sensors (*BSs*), to determine the position of the target. The coordinates of the sensor nodes are known, which are *si* = (*ai*, *bi*) *<sup>T</sup>*, *<sup>i</sup>* <sup>∈</sup> {1, 2,..., *<sup>N</sup>*}, where [·] *<sup>T</sup>* denotes the matrix transpose. Assume that the target's coordinate is *p* =(*x*, *y*) *T*.

As shown in Figure 1a, there are three basic sensors in the 2-D plane to determine the position of the target which form two groups of hyperbolas [24]. The hyperbola has two intersections in the absence of noise and there is one ambiguous position in them. When noise exists, the other two groups of hyperbolas have the other two intersections and both intersections have errors. In order to avoid the ambiguous position, it is advisable to increase the number of the sensors. As demonstrated in Figure 1b, the four basic sensors form three groups of hyperbolas and there is only one intersection without noise, which is the estimated position of the target. When noise exists, it is necessary to follow certain principles to obtain the optimal results.

**Figure 1.** The principle of time difference of arrival (TDoA) measurement. (**a**) The diagram when there are three basic sensors of TDoA. (**b**) Multiple hyperbolas for the optimal position.

Take the first basic sensor *BS*<sup>1</sup> as a reference sensor and assume that the signal propagates in a straight line between the target and each basic sensor without considering the influence of non-line-of-sight propagation. Assume that the times when the signal arrives at basic sensors *BS*<sup>1</sup> and *BSi* are *t*<sup>1</sup> and *ti*, respectively, and the propagation speed of the signal is *c*. The range of difference between the target and two basic sensors *BS*<sup>1</sup> and *BSi* is " *ri*,1# . This paper assumes that range difference errors {*ni*} are independent Gaussian random variables with zero mean and known variance <sup>σ</sup><sup>2</sup> *<sup>i</sup>* , i.e., N 0, σ<sup>2</sup> *i* . We can obtain

$$r\_{i,1} = \mathfrak{c}|t\_1 - t\_i|\tag{1}$$

$$r\_{i,1} = d\_{i,1} + n\_{i,1}, i \in \{2, \ldots, N\}. \tag{2}$$

Thus,

$$|c|t\_1 - t\_i| = |d\_{i,1} + n\_{i,1}| \tag{3}$$

where *di*,1= *di* − *d*1. Here, distances between the target and the receiver pair *BS*<sup>1</sup> and *BSi* can be expressed as follows:

$$d\_1 = \sqrt{\left(x - a\_1\right)^2 + \left(y - b\_1\right)^2} \tag{4}$$

$$d\_i = \sqrt{(x - a\_i)^2 + \left(y - b\_i\right)^2}, \ i \in \{2, \ldots, N\}. \tag{5}$$

Actually, the process of obtaining results based on TDoA measurements is the process of solving the *N* − 1 equations as shown in Equation (3) and obtaining the optimal solution.

### **3. WLS Method**

Usually, there are iterative methods, such as those mentioned in Section 1, to solve the equations, for which the computational burden is heavy. In this section, the WLS method is introduced based on TDoA measurements [29]. The sum of squares of residuals is defined as *JNLS*(%*x*):

$$f\_{\rm NLS}(\widetilde{\boldsymbol{\alpha}}) = \min \sum\_{i=1}^{N} R\_i^2(\widetilde{\boldsymbol{\alpha}}) \tag{6}$$

where %*<sup>x</sup>* represents the optimization variable, and residual *Ri* (%*x*) can be expressed as

$$R\_j(\overleftarrow{\bf x}) = \overleftarrow{r\_{i,1}} - r\_{i,1} \tag{7}$$

where %*ri*,1 is the measured value. Therefore, the optimal solution *<sup>p</sup>*<sup>ˆ</sup> according to the principle of minimum variance is

$$\mathfrak{p} = \arg\min\_{\mathbf{x} \in \mathbb{R}^2} I\_{\text{NLS}}(\mathbf{\overline{x}}). \tag{8}$$

Nonlinear hyperbolic equations can be transformed as follows:

$$
\sigma\_{i,1} + \sqrt{\left(\mathbf{x} - a\_1\right)^2 + \left(y - b\_1\right)^2} = \sqrt{\left(\mathbf{x} - a\_i\right)^2 + \left(y - b\_i\right)^2} + n\_{i,1}, \ i \in \{2, \ldots, N\}. \tag{9}
$$

After mathematical transformation, we can obtain

$$a\_i(x - a\_1)(a\_i - a\_1) + (y - b\_1)(b\_i - b\_1) + r\_{i,1}d\_1 = \frac{1}{2} [\left(a\_i - a\_1\right)^2 + \left(b\_i - b\_1\right)^2 - r\_{i,1}^2] + d\_i n\_{i,1}, i \in \{2, \dots, N\} \tag{10}$$

where the second-order term *ni*,1<sup>2</sup> is ignored and *ei*,1= *dini*,1. We can obtain

$$AX = \Theta + e$$

in which

$$A = \begin{bmatrix} a\_2 - a\_1 & b\_2 - b\_1 & r\_{2,1} \\ a\_3 - a\_1 & b\_3 - b\_1 & r\_{3,1} \\ \vdots & \vdots & \vdots \\ a\_N - a\_1 & b\_N - b\_1 & r\_{N,1} \end{bmatrix} \tag{12}$$

$$X = \begin{bmatrix} x - a\_1 \ y - b\_1 \ d\_1 \end{bmatrix} \tag{13}$$

$$\Theta = \frac{1}{2} \begin{bmatrix} \left(a\_2 - a\_1\right)^2 + \left(b\_2 - b\_1\right)^2 - r\_{2,1}^2\\ \left(a\_3 - a\_1\right)^2 + \left(b\_3 - b\_1\right)^2 - r\_{3,1}^2\\ \vdots\\ \left(a\_N - a\_1\right)^2 + \left(b\_N - b\_1\right)^2 - r\_{N,1}^2 \end{bmatrix} \tag{14}$$

$$\begin{vmatrix} \mathbf{i} \\ (a\_N - a\_1)^2 + (b\_N - b\_1)^2 - r\_{N,1}^2 \end{vmatrix}$$

$$\boldsymbol{\varepsilon} = \begin{bmatrix} \boldsymbol{\varepsilon}\_{2,1} \ \boldsymbol{\varepsilon}\_{3,1} \cdots \ \boldsymbol{\varepsilon}\_{N,1} \end{bmatrix}^T. \tag{15}$$

Then, the WLS objective function can be expressed as

$$J\_{WLS}(X) = \left(AX - \Theta\right)^{\mathrm{T}} \mathcal{W}(AX - \Theta) \tag{16}$$

where the weighting matrix is *<sup>W</sup>* = *E* " *eeT* #−<sup>1</sup> .

Thus,

$$\text{At}\_{WLS} = \left(A^T W A\right)^{-1} A^T W \Theta. \tag{17}$$

The WLS method is often adopted because of its simplicity and lower computational burden. As the equations are approximated in the process of simplification, it has low accuracy.

### **4. Firefly Algorithm**

In this section, the principle of the FA is introduced, which is a kind of heuristic algorithm inspired by the flickering behavior of fireflies. The three following idealized rules are needed for the model of this algorithm [30]:


Through the attraction between the brighter and less bright fireflies, the fireflies will eventually gather around the brightest firefly, which can realize the optimization of the objective function. We can use FA search methods to obtain the optimal result that satisfies Formula (18):

$$J\_{NLS}(\overleftarrow{\mathfrak{x}}) = \min \sum\_{i=1}^{N} \mathcal{R}\_i^2(\overleftarrow{\mathfrak{x}}). \tag{18}$$

In the search space, fireflies move towards brighter fireflies continuously to complete optimization until the preset termination condition of the algorithm is reached.

Assuming that the number of fireflies is *N* and the dimension is *D*, the positions of the *i th* and *j th* fireflies are *xi* <sup>=</sup> (*xi*1, *xi*2, ··· , *xiD*), *<sup>i</sup>* <sup>=</sup> 1, 2, ··· , *<sup>N</sup>* and *xj* <sup>=</sup> *xj*1, *xj*2, ··· , *xjD* , *j* = 1, 2, ··· , *N*. *rij* is the distance between the *i th* and *j th* fireflies, which can be calculated as follows:

$$r\_{i\bar{j}} = \|\mathbf{x}\_{i\bar{i}} - \mathbf{x}\_{\bar{j}}\| = \sqrt{\sum\_{d=1}^{D} \left(\mathbf{x}\_{id} - \mathbf{x}\_{\bar{j}d}\right)^{2}}.\tag{19}$$

Among them, *xid* and *xjd* represent the positions of the *i th* and *j th* fireflies, respectively. The relative brightness of the fireflies is defined as

$$I = I\_0 \mathfrak{e}^{-\gamma r\_{ij}\mathfrak{z}} \tag{20}$$

where *I*<sup>0</sup> represents the brightness of the firefly, which is proportional to the value of the objective function. γ is the coefficient of absorbing light intensity, which is usually defined as a constant. *rij* denotes the distance of fireflies *i* and *j*.

The attractiveness of the fireflies is defined as follows:

$$
\beta = \beta\_0 e^{-\gamma r\_{ij}^{-2}} \tag{21}
$$

where β<sup>0</sup> denotes the factor of maximum attraction degree, indicting the attractiveness of the position with the maximum brightness. From Formula (21), we understand that attractiveness decreases with the increase of the distance and the coefficient of absorbing light intensity.

The update of the location is

$$\mathbf{x}\_{\rm id}(t+1) = \mathbf{x}\_{\rm id}(t) + \beta \left(\mathbf{x}\_{\rm id}(t) - \mathbf{x}\_{\rm id}(t)\right) + \boldsymbol{\varepsilon} \cdot \boldsymbol{\alpha}\_{\rm i}(t) \tag{22}$$

where *xid*(*t*) and *xjd*(*t*) are the positions of the *i th* and *j th* fireflies after the *t th* generation. α*i*(*t*) denotes the step factor of the *t th* generation. The range of the value ε is [−0.5, 0.5], which sequences with uniform distribution.

The process of optimization is as follows. Fireflies with varying degrees of brightness are randomly dispersed in the solution space. The brightness and attractiveness of the fireflies can be calculated according to Equations (20) and (21), respectively. The less bright fireflies will move towards the brighter one. In order to avoid falling into the local optimum, the perturbation term ε · α*i*(*t*) is added to the process of location updating. Finally, the fireflies will gather around the firefly with highest brightness. The optimal result can thus be obtained. The flowchart for this can be found below.

**Step 1.** Initialize the parameters in the algorithm. Set the number of fireflies *N*, the factor of maximum attraction degree β0, and maximum iteration number or convergence criterion.

**Step 2.** Initialize the location of the fireflies randomly and calculate the value of the objective function as the original brightness.

**Step 3.** Calculate the brightness and attractiveness of fireflies referring to Equations (20) and (21), respectively, and determine the moving direction of the fireflies according to their relative brightness.

**Step 4.** Update the location of the fireflies according to Equation (22) and add the perturbation terms.

**Step 5.** Recalculate the brightness of the fireflies after updating the location of the fireflies.

**Step 6.** When the convergence criterion is satisfied or the maximum number of iterations reached, go to the next step; otherwise, go to Step 3.

**Step 7.** Output the global extremum and optimal value.

### **5. Hybrid-FA Method**

While it usually takes more time to use a search algorithm than iterative methods, this method provides higher accuracy and thus has great potential in practical applications. There are certain reasons for the longer solution time. One significant cause is that it will search the optimal result in the global scope. In this context, if some reasonable regional restrictions are given, a search algorithm, including the FA method, will reduce the computation amount as well as ensure the accuracy of the result.

In this study, the WLS and FA methods were combined for optimal implementation. Since it is easy to obtain the initial result by the WLS method, the initial result can be used to provide the limited area for the FA method.

Assume that the search area is square and the length of a side is *l*. If the result obtained by the WLS method is (*xwls*, *ywls*), which is the initial value, then the constrained region can be given for the FA method as [*xwls* <sup>±</sup> <sup>1</sup> <sup>4</sup> *<sup>l</sup>*] <sup>×</sup> [*ywls* <sup>±</sup> <sup>1</sup> <sup>4</sup> *l*] to ensure the target falls into the restricted region as much as possible. Additionally, new firefly positions out of the restricted region are ignored in this method.

As shown in Figure 2, the WLS and FA methods are combined for optimal implementation. The initial value is obtained by the WLS method, then the initial result can be used to provide the restrained area for the FA method to search for the optimal result. There are two conditions when the algorithm ends, fulfilling the convergence criterion or implementing maximum iterative times set in advance. As for the former condition, the values of the objective function obtained in the *<sup>i</sup>* <sup>−</sup> <sup>1</sup>*th* and *<sup>i</sup> th* iteration are compared. Assuming that the objective function is F, then the ending condition can be expressed as follows:

$$\|X(\mathbf{i}^{th}) - X(\mathbf{i} - \mathbf{1}^{th})\| \le \varepsilon\_1 \tag{23}$$

$$\|\|\mathbb{F}(\dot{\mathbf{r}}^{th}) - \mathbb{F}(\dot{\mathbf{r}} - \mathbf{1}^{th})\| \le \varepsilon\_2 \tag{24}$$

where ε<sup>1</sup> and ε<sup>2</sup> are positive predetermined numbers.

**Figure 2.** The diagram of hybrid firefly algorithm (hybrid-FA) method.

### **6. Results**

### *6.1. Preprocessing*

Simulations and indoor experiments were conducted and the results of the proposed method and other commonly used methods are presented and compared here. Before the simulations and experiments, the definition of the SNR and the evaluation index are given.

In this paper, the SNR of the signal is defined as

$$\text{SNR} = 10 \log \frac{d\_{i,1}}{\sigma\_i^2} \text{ dB} \tag{25}$$

where *di*,1 denotes the distance difference between the target and the basic sensors *BS*<sup>1</sup> and *BSi*, and σ*<sup>i</sup>* represents the standard deviation of the noise.

Therefore, when used in practice, if the SNR is known in advance, the variance of noise can be obtained:

$$
\sigma\_i^2 = \frac{d\_{i,1}^{-2}}{10^{\text{SNR}/10}}.\tag{26}
$$

Otherwise, when the SNR is not known in advance, the variance of noise can be obtained using the approximate method. The performance index of the receiver is usually known according to the specifications or can be measured by testing. Here, the signal's time of arrival was measured by the receiver. Assume that the measurement error targeting at the time of arrival of different receivers is no more than Δˆ*ti*(*i* = 1, 2, ..., *M*), where M denotes the number of the receivers. Assume that the time difference of arrival between *BS*<sup>1</sup> and *BS*<sup>i</sup> is Δ*t*1,*i*.

Then, we can get

$$
\Delta t\_{1,i} = |t\_1 - t\_i| \tag{27}
$$

where *t*<sup>1</sup> and *ti* denote the measured values when the signal arrives at the basic sensors *BS*<sup>1</sup> and *BSi*, respectively.

The variance of the noise of Δ*t*1,*<sup>i</sup>* can be approximated to

$$
\sigma\_i^2 = \varepsilon^2 \left(\Delta \hat{t}\_1^{'2} + \Delta \hat{t}\_i^{'2}\right) \tag{28}
$$

where *c* represents the propagation velocity of the signal.

The localization performance was evaluated referring to the root-mean-square error (RMSE), which is defined as

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left[ \left( \pounds\_i - \mathfrak{x} \right)^2 + \left( \pounds\_i - \mathfrak{y} \right)^2 \right]} \tag{29}$$

where *n* represents the number of simulation times, (*x*, *y*) are the real position coordinates of the target, and (*x*ˆ*i*, *y*ˆ*i*) are the estimated positions based on the *i th* calculation. RMSE was used to measure the average coordinate distance between the estimated target position and the actual target position. The lower it is, the higher the accuracy is.

### *6.2. Simulation Conditions*

The simulation was performed using Matlab 2014a and all the results were obtained on the same computer with a 1.8 GHz CPU and 8 G RAM. Assume that the coordinates of the four basic sensors were *BS1*(0,0), *BS2*(0,10), *BS3*(10,10), and *BS4*(10,0), where *BS1* is the reference basic sensor. The layout of the four sensor nodes is shown in Figure 3.

**Figure 3.** The layout of four sensor nodes in the simulation experiment.

### *6.3. The Robustness*

The results of the simulations were as follows when the SNR = 30 dB. In order to have a better display, the vertical coordinate represents -RMSE. Therefore, the optimal result was the coordinate when the -RMSE reached the maximum.

As illustrated in Figure 4, the FA method searched for the optimal result in the global region ([0,10] × [0,10]). From the figure, we can see that the result obtained by the FA method was even closer to the actual target location than that obtained by the WLS method, while the hybrid-FA method simply searched in the square marked in red. In this connection, it reduced the computational burden and maintained high accuracy.

restricted region by Hybrid-FA method

**Figure 4.** The diagram of hybrid-FA and FA.

Table 1 shows the comparison between FA and hybrid-FA under the same condition when SNR = 30 dB. When the RMSE reached 0.03741 m, the hybrid-FA method only needed to iterate 30 times, while 100 times was required for the FA method by itself. The simulation results demonstrate the efficiency of the hybrid-FA method.

**Table 1.** The root-mean-square error (RMSE) of the hybrid-FA and FA methods with different numbers of iterations.


In this section, the results of the commonly used algorithms CWLS [31], Newton–Raphson (NR) [29], TSWLS [18], and GA [30], were compared with the hybrid-FA. The GA method is a search algorithm that is commonly used for optimization. For these methods, the number of iterations was 30 and the coordinate of the target was set as (2,3). Figure 5 illustrates the compared results.

**Figure 5.** The comparison of the four algorithms for TDoA measurement.

As demonstrated in Figure 5, with the increase of SNR, the RMSE of each algorithm tended to decrease, which means the accuracy of the position had improved. When the SNR = 10 and 15 dB, the RMSE of TSWLS was the lowest, while the RMSE of the other five methods was typically higher than 2 m. It was hard to achieve high accuracy when the SNR was very low, as the acquired data was limited. The data of FDoA or TOA should be combined to get higher accuracy. When the SNR = 20, 25, 35, and 40 dB, the RMSE of the hybrid-FA was the lowest, which was 0.46077, 0.03788, 0.02427, and 0.0177 m, respectively. When SNR = 30 dB, the RMSE of the CWLS method was 0.0001 m lower than that of hybrid-FA method. On the whole, the hybrid-FA method was better than the NR, TSWLS, and GA methods when the SNR ranged from 20 to 40 dB in the simulations. The result of GA is also shown in Figure 5. The RMSE of GA was higher than that of the proposed scheme, which illustrates that the proposed scheme is more appropriate for optimization than GA.

### *6.4. Experiment*

Experiments were carried out to verify the rationality of the algorithms [23,32]. The coordinates of four speakers were *BS1*(0,0)m, *BS2*(0,10)m, *BS3*(10,10)m, and *BS4*(10,0)m, which were used to generate signals. The speakers emitted chirp signals, which were continuous impulse signals of 2.5 kHz. A phone was placed at the same height to receive the sound signal as well as record the receiving time. In this experiment, time-division multiplexing was adopted. The emitting cycle was 1 s, and the speakers emitted 100 ms long signals one after another. The speaker of *BS1* emitted 100 ms long signals at the beginning of every emitting cycle. The speakers of *BS2*, *BS3*, and *BS4* emitted 100 ms long signals at the 250th, 500th, and 750th ms, respectively. Take the speaker of *BS1* as a reference speaker. The receiving time data were saved to a text file and exported to the computer to be processed in Matlab 2014a.

Assume that *t j <sup>i</sup>* represents the time when the phone receives the signal from the *BSi* speaker at the jth emitting cycle. The time difference of arrival between *BS1* and *BSi* can be expressed as

$$
\Delta t = \left| t\_i^j - t\_1^j - \frac{i-1}{4}T \right| \tag{30}
$$

where *T* denotes the emitting cycle. Thus, the range of difference between the receiver and speakers *BS1* and *BSi* can be expressed as

$$r\_{i\_r1}^{\rho} = \mathcal{c}^{\rho} \Delta t$$

where *c<sup>o</sup>* denotes the propagation velocity of the signal. Then, the equations can be obtained and solved by the localization algorithms.

Firstly, the experiment was conducted to verify the performance of the five localization methods. There were 19 test positions of the receiver and the distribution of receiver occupied the search region as much as possible. At each test position, 50 trials were conducted under the same conditions. In total, 950 trials were carried out in this experiment. The sketch of the experiment is shown in Figure 6.

In Figure 6, the red points represent the position of the receiver and the distance of the two adjacent test points is 2.5 m. For the results, the trials with RMSE greater than 2.5 m were considered as bad results. The result of RMSE is the mean of the results from all trials for each method.

The RMSE of the five methods in this experiment and the amount of bad results are shown in Figure 7. The RMSE of the hybrid-FA method was 0.6778 m, which was lower than that of NR, TSWLS, and GA and 0.0031 m higher than that of CWLS. As for the amount of bad results, for the hybrid-FA method, it was 90, which was less than that of the NR, TSWLS, and GA methods and 3 more than that of CWLS. It can be concluded that the performance of the hybrid-FA method was superior to that of the NR, TSWLS, and GA methods for TDoA measurement.

**Figure 6.** The sketch of the first experiment.

**Figure 7.** The results of the first experiment. (**a**) The RMSE of the five methods. (**b**) The number of bad results from the five methods.

In the experiment, the phone moved along the red path slowly, as shown in Figure 8. The A series of data was recorded and the results were calculated according to the different methods. The amount of the test position was 19 in the experiment. The discrete position sequence was obtained, then the Kalman filter with the same parameters was used to smooth the motion trail. The final results are shown below.

Figure 9 illustrates the trajectory tracking of the CWLS, hybrid-FA, NR, TSWLS, and GA methods. In Figure 9, the blue point is the discrete position solved by localization algorithms, the black line is the actual trajectory of the target, the red point is the estimated position obtained by smoothing the discrete position sequence using the Kalman filter, and the red line is the smoothed trajectory of the target. It was difficult to arrive at a conclusion solely through observation. For this reason, the mean distance error was introduced to compare the performances of the different methods. Assume that the coordinate of the smoothed position is *xo i* , *yo i* for each method and *d<sup>o</sup> <sup>i</sup>* is the distance between the smoothed position and the line y = x in the coordinate system, which is the actual moving path of the receiver. Thus, the mean distance error is defined as follows.

$$\text{mean distance error} = \frac{1}{\text{N}} \sum\_{i=1}^{\text{N}} d\_i^o \tag{32}$$

**Figure 8.** The setting of second experiment.

**Figure 9.** The trajectory tracking of the five methods based on TDoA. (**a**) The trajectory tracking of constrained weighted least squares (CWLS). (**b**) The trajectory tracking of hybrid-FA. (**c**) The trajectory tracking of Newton–Raphson (NR). (**d**) The trajectory tracking of two-step weighted least squares (TSWLS). (**e**) The trajectory tracking of the genetic algorithm (GA).

Table 2 shows the mean distance error of the CWLS, hybrid-FA, NR, TSWLS, and GA methods. The mean distance error of the hybrid-FA method in this experiment was 0.03419 m, which was less than that of the NR, TSWLS, and GA methods and 0.000985 m more than that of CWLS. It can be concluded that the hybrid-FA method outperformed the NR, TSWLS, and GA methods for TDoA measurement.

**Table 2.** The mean distance error of different methods.


### **7. Conclusions**

For TDoA measurement, a good algorithm should balance calculation and precision. In this paper, a hybrid-FA method was proposed that combined the WLS and FA methods, which used the result from WLS with low computational burden to provide a reasonable limit to the search region for the FA method. The results of the proposed method were compared with the CWLS, NR, TSWLS, and GA methods using simulations and two experiments, which demonstrated the validity and limitations of the proposed method.

As expected, the hybrid-FA method could cut down the computation of the algorithm with high accuracy compared with using the FA only. Additionally, the hybrid-FA method was compared with the CWLS, NR, TSWLS, and GA methods using simulations and experiments. The RMSE of the hybrid-FA method was lower than that of the NR, TSWLS, and GA methods when the SNR ranged from 20 to 40 dB in the simulations. The result of the first experiment showed that the RMSE of the hybrid-FA method was 0.6778 m, which was lower than that of NR, TSWLS, and GA. The results of the second experiment illustrated that the mean distance error of the hybrid-FA method was 0.03419 m, which was lower than that of NR, TSWLS, and GA. On the whole, the hybrid-FA method outperformed the NR, TSWLS, and GA methods for TDoA measurement.

**Author Contributions:** Funding acquisition, S.S.; Methodology, P.W.; Resources, X.W.; Software, X.G.; Supervision, B.S.; Validation, Z.Z.

**Funding:** This work was supported by National Natural Science Foundation of China under Grant No. 51575517, and the Natural Science Foundation of Hunan Province under Grant No. 2018JJ3607.

**Conflicts of Interest:** The authors declare no conflict of interests.

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Letter*
