A Hybrid Global Minimization Scheme for Accurate Source Localization in Sensor Networks

We consider the localization problem of multiple wideband sources in a multi-path environment by coherently taking into account the attenuation characteristics and the time delays in the reception of the signal. Our proposed method leaves the space for unavailability of an accurate signal attenuation model in the environment by considering the model as an unknown function with reasonable prior assumptions about its functional space. Such approach is capable of enhancing the localization performance compared to only utilizing the signal attenuation information or the time delays. In this paper, the localization problem is modeled as a cost function in terms of the source locations, attenuation model parameters and the multi-path parameters. To globally perform the minimization, we propose a hybrid algorithm combining the differential evolution algorithm with the Levenberg-Marquardt algorithm. Besides the proposed combination of optimization schemes, supporting the technical details such as closed forms of cost function sensitivity matrices are provided. Finally, the validity of the proposed method is examined in several localization scenarios, taking into account the noise in the environment, the multi-path phenomenon and considering the sensors not being synchronized.


Introduction
A challenging and highly demanding signal-processing application is the localization of signal sources using the physical measurements at some sensors in the environment. Source localization has become an important task in various applications such as mobile communications, global positioning system (GPS), radar, sonar, navigation, seismology, and geophysics [1][2][3][4][5].
During the recent decades various algorithms have been proposed to estimate the location of the signal sources. These methods utilize different signal characteristics at different sensors and generally can be classified in three main categories: using the time difference of arrival (TDOA); analyzing the signal direction of arrival (DOA) at distinct arrays; and using the differences in the signal amplitude or received energy level. For a constant propagation speed, the TDOA among different sensors is proportional to the source-sensor range differences and may be estimated through methods such as cross-correlation (CC) [6] or its generalized version (GCC) [7]. The source locations can then be estimated using geometric methods such as linear, spherical or hyperbolic intersections [8][9][10]. To estimate the DOA, for narrowband signals, high resolution algorithms such as multiple signal classification (MUSIC) [11] and maximum likelihood (ML) [12] are proposed. In [13], those authors propose an approximate maximum likelihood method (AML) for wideband signals using spectral properties of the signal when rather long sample streams are available. In this method, the corresponding cost function can be directly expressed in terms of the source locations or in a far-field case, may be expressed in terms of the relative time delays followed by a postprocessing step to find the source locations from the corresponding DOAs. The post-processing step may be carried out through geometric methods such as cross bearing or a machine learning approach such as the support vector machine (SVM) method [14]. Using the differences in the signal intensity or energy level for the purpose of localization is a more recent technique [15,16]. Theoretically, this class of localization can be considered for both narrowband and wideband signals by only taking into account the attenuation information and usually neglecting the time delay information. For these methods, a precise attenuation model in the environment is inevitable for an accurate localization.
Moreover, from an optimization perspective the resulting cost functions in these kind of approaches usually undergo many local optima and saddle points which require considering specific optimization schemes [17].
In this article, we manage the problem of localization of multiple wideband sources by coherently taking into account the TDOA and the amplitude attenuation pattern. Our method generalizes the AML approach to utilize the signal attenuation characteristics. We provide a more robust algorithm in which the targets should simultaneously satisfy the correct time delays among the sensors and provide sensible level of attenuation at each sensor. Unlike the aforementioned energy and intensity based methods which not only ignore the time delay stamps but also require knowing the signal attenuation model, we benefit using the delay information and as a generalization to our recent study in [18], leave the space for not knowing an exact signal attenuation model in the environment by suggesting an appropriate functional space for it. We minimize a cost function which is obtained through maximum likelihood approach from which the locations, attenuation model parameters, and the multi-path parameters are obtained. To apply the minimization, we propose a hybrid approach combining the differential evolution algorithm [19] with the Levenberg-Marquardt algorithm [20]. This combination provides a minimization scheme which is likely to globally search for the optima and rather quickly converges to the accurate results. Through simulations and Cramér-Rao bound, we verify the effectiveness of the novel method introduced in this paper.
This article is organized as follows. In Section 2, we propose a general form for the received signal at every sensor and later provide an adaptive model for the signal attenuation based on Laurent polynomials. In Section 3, a maximum likelihood estimation of the source location and attenuation parameters is proposed. We also provide the Cramér-Rao bound for this estimation problem. For the purpose of minimization in Section 4, a hybrid approach combining the differential evolution algorithm with the Levenberg-Marquardt is proposed for which the combination algorithm and closed form equations for calculation of the Jacobian are provided. In Section 5, we examine the efficiency of proposed method through some examples, and finally there are some concluding remarks in Section 6.

Signal model
Although the general approach proposed in this article is in theory independent of signal nature and the type of sensors used, in order to make reasonable simulations we consider acoustic source localization. Consider N acoustic sources having unknown locations r S n . Each source is omni-directionally emitting a signal s n (t), n = 1, ..., N at the time frame t. We also consider M acoustic microphones that are placed in known positions r M m , m = 1, ..., M, in the same environment. For every source in the environment, the function that describes the signal attenuation at a specific point is a(r), where r is the distance from the point to the source. In general, the signal attenuation function may be a function of various parameters such as signal frequency, medium inhomogeneity, etc. To simplify the problem, in this article we consider this function to be an identical form for all sources and solely function of the distance to the source. However, unlike some previous energy-based localizations [16,17] in which the attenuation is known to be proportional to r -1 , the actual form of a(·) is considered unknown function here. This type of modeling provides an additional flexibility to the problem for more realistic scenarios where the inverse proportionality of a(·) to r is violated because of other parameters, such as signal bandwidth and medium inhomogeneity. Considering s n (t -1 × N s /v) to be the signal measured one length unit away from every source, with N s being the samples per second and v being the propagation speed, ideally the overall received signal samples from the acoustic sources at every microphone is modeled as where ρ m,n = ||r M m − r S n || is the distance from n th source to m th microphone, and τ m,n = r m,n N s /v is the corresponding time samples delay in receiving the signal. The received signal in (1) is normalized to each microphone gain level to decrease the number of notations. A more realistic model takes into account the noise in the environment and also the signals going through a multi-path channel before arriving at every sensor; hence, we rewrite the received signal as Pm,n p=1 γ m,n,p s n (t −τ m,n,p ) + w m (t). (2) The term w m (t) represents the background noise which is considered to be a zero-mean white Gaussian with variance s 2 for the sake of this article: Gaussianity is not a limiting assumption in this article. Between the n th source and m th microphone, we consider Pm , n indirect paths each causing g m,n , p loss in the signal amplitude andτm,n,p delay in the signal reception, modeling the multi-path phenomenon. Beside the positions r S n , which are the main unknowns of the localization problem, the signals s n (t), the multi-path parameters g m,n,p andτm,n,p , and the propagation loss function a(·) are also unknown and should be determined based on the received signals at the sensors. The appearance of τ m,n (which is related to the unknown quantities r S n ) andτm,n,p as the argument of an unknown signal s n (t) causes an extra complexity for any optimization scheme performed to solve the localization problem. However, this problem may be remedied by applying the discrete Fourier transform to (2) to extract the delays and form an equivalent equation in which the unknown parameters are separated in individual terms: For where, X m (f), S n (f) and ξ m (f) are the data, signal, and noise spectra respectively. As stated in [13], we emphasize on the fact that, for (3) to be a valid equivalent form of (2), we need n t to be large enough to avoid edge effects and accordingly n f >n t . In general, having more samples from the signal better poses the problem.

A low-order representation of signal attenuation model
As discussed earlier, our assumption about the attenuation model in the environment in this article is an identical model for all sources, which is only dependent on the distance of the point to the acoustic source. In an ideal environment, a(r) can be well modeled as a multiple of r -1 . Since there are different parameters involved in the attenuation modeling, a(r) is being considered as an unknown here. However, in order to keep the wellposedness of the problem, we choose it to be an element of a low-dimensional function space. For this sake, we consider a(r) to be a Laurent polynomial of limited order as In this model, only negative powers of r are considered, which is because for an attenuation model, we are physically required to have In many applications, the low-order representation of a(r) in (4) is acceptable enough to model the attenuation, and usually considering only few terms in the series (i.e., L rather small), would suffice for the localization problem.

Derivation
Based on the general attenuation model proposed, matching of the data spectra with those of the model can be expressed by using (4) in (3) as The central limit theorem states that ξ m (f), which is a transformed zero mean Gaussian random variable to the frequency domain, should be a complex zero mean Gaussian with variance n t s 2 . For every frequency bin f having (6) can be written in a matrix form as where for which ⊗ represents the Kronecker product, and I N × N the identity matrix of size N × N, and where The matrix H(f) is related to the multipath parameters and its elements are Rewriting the negative log-likelihood function to estimate the unknown parameters θ including the source positions, source signal spectrums, multi-path parameters and quantities b ℓ , we have and Similar to the idea in [13], for a real valued signal, we can only consider up to n f /2 frequency bins and form Q with blocks of Q(f) for f = 0, 1, ..., n f /2. We would like to highlight the fact that in [13], the zero frequency bin is ignored due to producing a constant term in the likelihood function; however, in our approach, the matrices K(0) and H(0) are still dependent on r m,n and the multi-path parameters g m,n,p and hence are worth being considered.
Clearly, the minimization in (11) is equivalent to minimizing Q H (f) Q(f) for every f. Considering the unknown signal spectrum S(f), the minima should satisfy for f = 0, ..., n f /2, and therefore the unknowns of the minimization reduce to the source positions, multi-path parameters, and the attenuation coefficients. Considering a 2D localization problem, as the case in the example section, neglecting the multi-path the vector of unknowns would be where x S n and y S n are, respectively, the x and y components of the position vector r S n . In case of multi-path, the parameters g m,n,p andτm,n,p are also included in θ. The approach is clearly not only limited to 2D Cartesian systems and 3D Scenarios, and other coordinate systems may also be considered.

Cramér-Rao lower bounds for the estimated parameters
For an unbiased parameter estimation problem, the Cramér-Rao lower bound (CLRB) is a theoretical lower bound on the variance of the problem estimates. Based on (7), the total model relating the parameters of interest to the complete dataset is where .., S(n f /2) T ] T contains the signal spectra of all the sources, and ξ = [ξ(0) T , ..., ξ(n f /2) T ] T is the corresponding noise vector. Moreover, G(θ; S) =KS for which the matrixK explicitly dependent on θ is The Cramér-Rao lower bound is defined as the diagonal elements of the inverse Fisher matrix F, which for the model in (17) is representable as [21] where and R ξ is the noise covariance matrix which for our problem is simply n t s 2 I. The matrix ∂G/∂ϑ is composed of the sub-blocks ∂G/∂S , ∂G/∂r S n , ∂G/∂β , ∂G/∂γ m,n,p , and ∂G/∂τ m,n,p Clearly For the θ parameters, sinceK is composed ofK(f ) , we only discuss the sensitivity ofK(f ) to every class of parameters. Based on the fact that where To calculate ∂K(f )/∂x S n , we have where The matrix ∂R (f )/∂x S n is a matrix with the same size as R ℓ (f), with all columns being zero except the n th column. Simply applying the derivative shows that the (m, Specifying the elements of the Fisher matrix F yields the CRLB values for all the estimations.

Minimization strategy
The minimization in (11) may be performed through various optimization schemes, most generally categorized as global and local optimizations. For a global optimization different approaches such as deterministic, stochastic, or evolutionary and metaheuristic methods may be considered [22][23][24]. Clearly for an accurate localization, global minimizers of (11) are required. However, in general, using global methods to optimize an arbitrary function may be iteratively or computationally expensive. As an alternative to this and specifically for a least squares cost function as (11), local search methods such as gradient descent and quasi-Newton methods may be considered [20]. Although these methods can be relatively faster than the global ones, there is always a chance of getting trapped into a local minima. In the context of localization, although for good initial estimates of the source relatively fast methods such as the gradient descent and alternating projection are proposed, to increase the chances of finding a global minima, the process usually involves exhaustive search methods such as the grid search and multiresolution search [13,16].
For the purpose of this article, we consider a hybrid approach combining a global search method with a fast local search method [25,26]. Hybrid methods have received considerable interests in different areas in the recent years [26][27][28][29]. More specifically, we consider a hybrid combination of the differential evolution (DE) algorithm [19] as successful evolutionary search with the Levenberg-Marquardt algorithm (LMA) [20,30] as a rather fast and robust local search method. Before getting to the combination scheme, we provide a brief description of each method highlighting the main technical issues specifically in the context of our localization problem.

Differential evolution algorithm
DE is among the metaheuristic and evolutionary global optimization schemes. Simplicity and successful performance are the main advantages of this algorithm. Considering θ = [θ 1 , θ 2 , ..., θ D ] to be the vector of problem unknowns with size D, at every generation G of the algorithm, N P parameter vectors θ i,G = [θ 1,i,G , θ 2,i,G , ..., θ D,i,G ], i = 1, 2, ..., N P , are generated. The initial population is randomly chosen with a uniform distribution in the search region. For this study, we consider the DE/ rand/1/bin, which is a general and widely used strategy of this algorithm [19,31]. For every generation, three main operations are performed as follows:

Mutation
In this phase, a mutant vector μ i,G is generated as where r 1 , r 2 , and r 3 are randomly selected indices among 1, 2, ..., N P , and F ε 0[2] is a constant real scalar controlling the difference vector amplification.

Crossover
A mixing with the mutant vector to increase the diversity of the population is performed by generating new trial vectors υ i,G of length D, defined as with d = 1, 2, ..., D. Here, C R ε [0,1] is the crossover constant, r(d) [0,1] is the d th evaluation of a uniform random number generator in [0,1] and k(i) ε {1, 2, ..., D} is a randomly chosen index ensuring that υ i,G takes at least one of the elements of μ i,G .

Selection
At this step a next generation population member θ i,G+1 is produced by a selection among θ i,G and υ i,G . This selection is based on the fitness, and basically, the vector with the lower cost proceeds to the next generation.

A Levenberg-Marquardt algorithm for the local minimization
As the local minimization scheme, we suggest using the LMA. Our attention toward this algorithm is based on several advantages. LMA is basically considered as a Newton type method and provides a rather quadratic convergence. Meanwhile, this algorithm benefits from stability and uses a trust region approach [30]. The other feature of this method, considered as an advantage over other methods such as the gradient descent, is its suitability for cases where there are different variables of different types as the cost function arguments. In fact, LMA is almost independent of variable scaling, while for methods such as the gradient descent, minimizing a cost function dependent on a set of variables with different natures and scales requires appropriate parameter scaling to guarantee a proper convergence [30]. This is a demanding feature for our problem where the θ vector in general consists of the source locations, attenuations coefficients, and the multi-path parameters.
In the LMA, which is an iterative algorithm, we start with a θ (0) as the starting point. At every iteration, having θ (i) already in hand, θ (i+1) can be obtained by solving where Q is the vertical vector of length Mn f /2 shown in (12) and obtained for values θ (i) at that iteration. The parameter l (i) is the damping factor, obtained at every iteration based on the trust region approach [20,30]. The Jacobian matrix J θ contains the sensitivities of Q to every element of θ. In order to run the algorithm, we need to know the Jacobian matrix at every iteration, obtaining which is discussed in the Appendix.

The hybrid combination scheme
For the purpose of combining the DE with the LMA, we propose using a sequential hybridization approach [26]. In this approach, the DE initially starts the minimization by generating consecutively more fitting generations. After a certain number of generations or after getting relatively slow in decreasing the fitness, the best θ in the last generation is passed to the LMA algorithm as an initialization. The minimization continues until convergence. An illustration of this algorithm is provided in Figure 1.

Simulation results
To examine the method developed in the previous section, we consider some localization examples in this section. In the first example, we consider a reverberation-free environment to show the efficiency of the method for such cases and provide a comparative study for this scenario. The second example brings more realistic issues such as the multi-path, and sensor synchronization error into the problem, and examines the performance of the proposed method for such cases.
Before proceeding with the examples, we would like to highlight a fact regarding the relationship between the cost function and the matrixK . Referring to (15), it can be easily verified that scalingK(f ) by a scalar does not change the cost function. In other words, if the b ℓ and g m,n,p values are simultaneously scaled by a scalar value, then the cost function remains the same. Therefore, we rewrite the attenuation model in (4) as which somehow normalizes a(·) and unifies the representation. Clearly, since the desired unknowns of the problem are the acoustic source coordinates, obtaining a multiple of the attenuation and multi-path coefficients is non-problematic. The true attenuation model to be used in this article is a(r) = r -1.25 [15].

Example 1
For the purpose of this example, we consider the sensors to be placed in the first quarter of the x-y plane as a spiral array of M = 40 microphones. The spiral is represented in a parametric form as where the angles s are equally spaced in [2π, 4π]. Our purpose of choosing such sensor arrangement was to provide a non-symmetric and still reproducible arrangement. The sensor locations are shown in Figure 2. The sources used in this example are wideband sources with center frequency 500 Hz and 200 Hz bandwidths. The sampling frequency is 4 KHz. The number of samples available from the sources at every sensor is n t = 4000 (i.e., the signal duration is 1 s), and the number of frequency bins is taken to be n f = 4100. The signal to noise ratio (SNR) at every sensor is 20 dB, and the speed of propagation is considered to be the speed of sound as v = 345 m/s. In the proposed minimization scheme and more specifically the DE part, we take G max = 5. Moreover, we set F = 0.8, C R = 1, and N P = 40. This parameter setting was selected as a general DE setting; however, more discussions on determining the DE parameters are available in [19]. The general attenuation model is considered to be a(r) = r -1 + b 1 r -2 + b 2 r -3 , for which the values b 1 and b 2 are in charge of tuning the unknown model. There is no reverberation in the environment (i.e.,K = K ), and all sensors are synchronized in receiving the signal.
To provide a better understanding of the problem, in Figure 2, the cost function behavior for a known attenuation model is shown. In Figure 2a, the cost is shown when the source is located at point (4, 3) within the sensors convex hull. All positions are in meters. Figure 2b shows the cost when the source is located at (12,10) outside the sensors region. In both cases the cost functions are rather well-behaved functions away from the sensors. Intuitively, for two neighboring points in the domain, sudden variation of the cost function with respect to both time delay criteria and attenuation model constraints is unlikely, and hence the resulting cost functions are usually expected to be rather slow varying and well behaved, away from the sensors.
In Figure 3, we have shown the iterative procedure of finding a single source, once located at (4, 3) and once at (12,10). For the first case, the source location is estimated to be at (4.002, 3.008) and the attenuation coefficients are estimated to be b 1 = −23.85 and b 2 = 27.93. In the second case the source estimation is (11.999, 10.002) and the attenuation coefficients are found to be b 1 = 4.19 and b 2 = 1.79. We observe that both localization results accurately match the exact source positions. The attenuation coefficients obtained in both cases are only in charge of fitting the low order model to the true model for the source-sensor ranges in each problem and due to different problems they do not necessarily need to be in the same ranges. By providing this low-order The cost function corresponding to a source located at (12,10).
attenuation model, we provide the flexibility to the problem for accurately estimating the sources. As a more challenging problem, we consider concurrent localization of the two sources located at (4, 3) and (12,10). Figure 4 shows the iterative procedure of finding the sources. The estimated source locations are (4.000, 3.001) and (11.989, 9.993), and the attenuation parameters are estimated to be b 1 = 4.77 and b 2 = −3.94. Again, an accurate match between the exact source locations and the estimated ones is observable.
We further examine our proposed method through a comparison with the AML method developed in [13]. For this purpose, we start reducing the signal samples by reducing the signal duration from 1 to 0.1 seconds and observing the error caused in the source estimation. Here, we consider the single source localization for the source being located at (12,10). Figure 5 shows the resulting error, as the signal duration decreases in both methods. As it is clearly observed, using both time delay and attenuation information helps our method provide better estimates of  the source locations even with less available data, compared to the AML method which only uses the time delay information. We further examine the performance of both methods for various SNR values. In Figure 6, the CRLB is calculated for the same single source scenario with the source located at (12,10). The RMS errors in estimating the x and y components of the source are obtained through 50 independent noise realizations for every SNR value shown in the figure. Again, the proposed method shows an acceptable performance regarding the closeness to the CRLB and the superior performance compared to the AML method.

Example 2
In a more realistic scenario, we examine the performance of the proposed method in a noisy environment where sensor synchronization error and reverberation are likely to happen. The sensor network configuration is shown in Figure 7, where three circular arrays each composed of 25 sensors centered at points (15,5), (2,15), and (5,28) are considered. The acoustic source is located at (35, 25), and the signal specifications are the same as the previous example. For this example G max is taken to be 20 to benefit more from a global search of a cost function which may not be as well behaved as the previous example because of bringing more unknown parameters into the problem. The low-order attenuation model considered in this example is a(r) = r -1 + br -2 , with b as the tuning parameter. Again, an SNR of 20 dB is considered at all sensors for all the experiments.
We first examine the case that the sensors are not exactly synchronized to receive the data. For this where ζ is a random variable uniformly distributed around zero. Equation (34) basically models the asynchronous measurements of the sensor data. In Table 1, we have provided the localization results for three different synchronization error variances 0.5, 1, and 2 ms. Clearly the phase error is a destructive phenomenon in TDOA localization algorithms; however, considering the localization errors in Table 1, one would observe that exploiting the attenuation information beside the phase information enables our algorithm to perform a rather accurate localization task in case of sensors being out of synchronism.
Furthermore, a more challenging problem is when the reverberation is also taken into account. In theory, for the emitted signal to arrive at every measuring sensor, an individual multi-path filter should be considered.
Although the formulation in this article is general, for the purpose of this example, we have made a reasonable and practical assumption that for all the sensors within each array, the filter representing the multi-path is identical. In general, the sensor network may be represented as a collection of several clusters with each being composed of sensors closely placed and each cluster treated as a single receiving node. This assumption prevents dealing with a large collection of unknowns (g m,n,p and τ m,n,p ) for every source-sensor pair and aggregates them into fewer parameters each assigned to the clusters.
To generate a reverberated signal, we use the multipath FIR filters shown in Figure 8 where three or four shifted scales of the signal are added to it. For the localization purpose, however, we only consider finding the main indirect path. In other words, for every array shown in Figure 7, only one multi-path coefficient g and one multi-path delayτ is to be estimated which totally brings six unknowns associated with the multi-path phenomenon into the minimization problem. The remaining minimization unknowns are the source coordinates and the attenuation coefficients as before. The fourth  row of Table 1 shows the localization result for this problem. Although the number of unknowns were relatively higher than the previous examples and the cost function is clearly not as well behaved as before, using DE as the initial minimization scheme provides a suitable starting state for the LMA, and this sequential technique helps the algorithm make a rather accurate localization in a noisy and reverberated environment. The fifth row of Table 1 corresponds to the case of having both the multi-path and the synchronization issues, for which the results are still promising. The progressive estimates of the target throughout the minimization are shown in Figure 9.

Conclusion
In this article, we proposed an efficient method for localization of multiple wideband sources based on both signal attenuation and time delay information. The method developed in this article models the localization problem as a minimization problem and provides an additional flexibility of not being exactly aware of the signal attenuation model. We propose a certain function space for the unknown model, and tune it iteratively along to estimate the signal source locations. The minimization scheme used here is a hybrid algorithm, combining the differential evolution with the LMA. This combination increases the chances of finding a global minima while benefits from the speed and computational advantages of Newton methods. The accuracy and performance of the method is examined through several simulations depicting a noisy environment, a multi-path environment, and lack of synchronization among sensors. In the simulations, we compared our approach with the approximate maximum likelihood method which show the superiority of the proposed method.   Figure 9 The progressive estimates of a single source located at (35, 25). The first jumps and good initials correspond to applying the DE.