SLAM-based Joint Calibration of Multiple Asynchronous Microphone Arrays and Sound Source Localization

In [12], based on the time difference of arrival (TDOA) between each pair of microphones, a calibration algorithm was developed to estimate the positions of microphones within a single microphone array. In [13], a bilinear calibration method based on time of flight (TOF) between each sensor source pair was proposed to estimate the microphone and source positions in 3D under the condition that the transmitting time is known. In [14], based on time of arrival (TOA) measurements and assuming knowledge of the distances between the sources and the microphones, a method for joint calibration of the positions of multiple microphones and sound source localization was proposed. In [15], a calibration method using TOA measurements was proposed for the scenario with a planar microphone array and a sound source moving in 3D.

Note that the applicability of the above-mentioned methods is limited in that they all rely on hardware synchronization between microphone channels, which is challenging to implement for robotic platforms in practice due to spatial and cost constraints [11]. Recently, in [16]-[18], a general framework using batch simultaneous localization and mapping has been developed for joint sound source localization and calibration of a single microphone array with asynchronous effects (i.e., clock difference and initial time offset).

Compared to single microphone array-based systems, the calibration of systems with multiple arrays has gained more recent attention. For example, the proposed approach in [19] utilizes direction of arrival (DOA) measurements to determine the sound source location and inter-array TDOA measurements to obtain the microphone array location through exhaustive grid search. The work [20] employs evolutionary algorithms to improve the accuracy and real-time performance of the approach in [19]. Based on DOA and inter-array TDOA measurements, another calibration framework for multiple microphone arrays is proposed in [21] using distributed damped Newton optimization. Note that the above-mentioned methods focus on the 2D case.

For the more general 3D case, there are only a few existing works. In [22], an artificial bee colony algorithm was employed to calibrate the positions and orientation of microphone arrays in 3D. Nevertheless, this method assumes that the sound source position at different moments is partially known and the clocks of the arrays are synchronized using hardware. Simultaneous calibration of positions, orientations, and time offsets of multiple microphone arrays and sound source positions in 3D was explored in [23] and [24].

For spatially distributed microphone arrays, it is necessary to consider both the initial time offsets and the sampling clock differences between the arrays [25], especially in the case of asynchronized scenarios based on the USB protocol and wireless acoustic sensor networks. In the above situations, each microphone array captures acoustic signals through its own microprocessor-controlled analog-to-digital converter and has a unique sampling clock source. Therefore, when launching multiple microphone arrays, differences in initialization result in varying initial time offsets between arrays. Moreover, the microprocessors in these microphone arrays often have limited performance, and the oscillators/crystals used to generate clock signals typically drift around their nominal frequencies. As a result, differences in sampling clocks accumulate over time. Not properly handling the above issue will significantly degrade the performance of sound source localization/tracking algorithms embedded in the arrays [17].

To the best of our knowledge, there is no work that has addressed the simultaneous calibration of positions, orientations, time offsets and sampling clock differences of multiple microphone arrays and sound source positions in 3D. In fact, as for single microphone array, calibration of multiple microphone arrays can be considered as a SLAM problem [26]-[29], where microphone arrays and the moving sound source serve as landmarks in the environment and the robot, respectively. As illustrated in Fig. 1, the acoustic measurements from the microphone arrays and the motion measurements from the robot are utilized in the optimization process, with landmark-robot constraints and robot relative pose constraints enforced, similar to the approach used in full information estimation and batch SLAM [26]-[31]. Then, two important questions arise.

Firstly, it is critical to assess whether the information contained in the measurements is sufficient to estimate the unknown parameters of microphone arrays and sound source locations. This is the so-called observability problem in the SLAM literature [32]-[33]. Although there exist works on observability analysis of SLAM-based calibration of single microphone arrays, in-depth analysis for the case with multiple microphone arrays is lacking.

Secondly, the selection of initial values is crucial because the considered calibration is a nonlinear least squares (NLS) problem, similar to batch SLAM [26], [34]. Many existing algorithms for solving such NLS problems employ the Gauss-Newton method or its variants. These methods typically require reasonable initial guesses; otherwise, the algorithms may converge toward local minima, or in extreme cases, diverge. For some specific problems, novel algorithms with certifiable convergence properties have been proposed in [35]-[37].

Motivated by the above observations, in this paper, we adopt batch SLAM as a general framework for the simultaneous calibration of translations, orientations, time offsets and sampling clock differences of multiple microphone arrays, and sound source positions in 3D. Our contributions are two-fold.

Firstly, we concentrate on the parameter identifiability of the corresponding SLAM problem. As discussed in existing works [32]-[33], SLAM is not observable from a control theoretical perspective. Hence, in the SLAM literature, the observability problem of SLAM has been tackled from an information-theoretic perspective, where all the parameters to be identified are taken to be constant but unknown. From the information-theoretic perspective, Fisher information quantifies the amount of information contained in a set of observations about a set of unknown parameters [33]. Following the above line of argument, when the multiple microphone array calibration problem is formulated as an NLS parameter estimation problem, the full rankness of the associated FIM determines the parameter identifiability or observability of the calibration problem.

Hence, in this paper, by leveraging the FIM approach, we thoroughly investigate the identifiability of the unknown parameters, including translations, orientations, and asynchronous factors between the microphone arrays and the sound source positions. We establish necessary/sufficient conditions under which the FIM and the Jacobian matrix have full column rank, which implies the identifiability of the unknown parameters. Furthermore, we identify several scenarios where the unknown parameters are not uniquely identifiable.

Secondly, we propose an effective framework to initialize the unknown parameters from the measurements, which is used as the initial guess in batch SLAM. Specifically, the initialization procedure is composed of the following major steps: (i) estimation of the sound source position by triangulation; (ii) estimation of distance between the sound source and microphone arrays using 3D geometry; (iii) estimation of microphone array poses using the iterative closest point (ICP) method; (iv) estimation of the asynchronous factors using linear least squares (LLS). As to be explained later in the paper (see Section IV. A), the microphone array pose estimation problem addressed in step (iii) mentioned above is conceptually a point-to-point registration problem, and hence can be tackled effectively using ICP [38]. To validate the effectiveness and robustness of the proposed initialization framework, we have conducted extensive numerical simulations and real experiments. Overall, the proposed pipeline achieves higher accuracy with fast convergence, in comparison to methods that use the noise-corrupted ground truth of the unknown parameters as the initial guess in the optimization, and other state-of-the-art methods in the literature [20], [23].

Compared to existing frameworks, the proposed calibration method requires less prior information. More specifically, the knowledge of the source’s position required in [22]-[23], or the distance between the signal source and the microphones needed in [14] is not required in the proposed framework in this paper. It should also be noted that our previous works documented in [16]-[18] primarily focused on calibrating individual microphones within a single array while in this paper we address the more challenging problem of calibrating multiple microphone arrays.

Finally, we remark that the observability analysis reported in Section III has been previously reported in our conference paper [39]. However, the results of [39] are only applicable for the case where the time interval between consecutive sound source events is fixed. In the current paper, we generalize the results in [39] from the scenario of fixed-interval sound source emissions to arbitrary time intervals (i.e., the interval between every two consecutive sound events can be asynchronous and time-varying). More importantly, we have proposed an effective framework for estimating the initial values of the parameters and conducted extensive simulation studies and real experiments to validate the entire calibration pipeline. All the codes and multimodal dataset used in this paper are publicly available at https://github.com/AISLAB-sustech/Calibration_of_Multi_Mic_Arrays.

Notation: Denote $x$, $\mathbf{x}$, and $\mathbf{X}$ as scalars, vectors, and matrices, respectively. $\mathbf{X}^{\mathrm{T}}$ represents the transpose of matrix $\mathbf{X}$. $\mathbf{I}_{n}$ stands for the identity matrix of $n$ dimensions. $\mathbb{R}^{n}$ denotes the $n$-dimensional Euclidean space. $[a_{1};\cdots;a_{n}]$ denotes $[a_{1}^{\mathrm{T}},\cdots,a_{n}^{\mathrm{T}}]^{\mathrm{T}}$, where $a_{1},\cdots,a_{n}$ are scalars/vectors/matrices with proper dimensions. $diag_{n}(\mathbf{A})$ denotes a block diagonal matrix with $\mathbf{A}$ as block diagonal entries for $n$ times; $diag(\mathbf{A},\mathbf{B})$ denotes a block diagonal matrix with $\mathbf{A}$ and $\mathbf{B}$ as its block diagonal entries; and $\mathbf{0}$ as a matrix of appropriate dimensions with its all entries as 0. $\mathbf{X}>0$ means that $\mathbf{X}$ is a positive definite matrix. We denote $\left\|\mathbf{x}\right\|_{\mathbf{P}}^{2}=\mathbf{x}^{\mathrm{T}}\mathbf{Px}$. Vectors/matrices, with dimensions not explicitly stated, are assumed to be algebraically compatible.

The covariance matrix $\mathbf{C}_{\hat{x}}$ of the estimation error corresponding to the estimated values $\hat{\mathbf{x}}$ and the true values $\check{\mathbf{x}}$ of unknown parameters in the observation model in (13) can be calculated by

For nonrandom vector parameter estimation, the FIM of an unbiased estimator is defined as:

where $\nabla_{x}$ is the gradient operator w.r.t. the vector $\mathbf{x}$, $p(\mathbf{z}|\mathbf{x})$ is the probability distribution function, and the derivatives are calculated at the true value $\check{\mathbf{x}}$ [40, chap. 2]. It can be shown that the covariance matrix of any unbiased estimator $\hat{\mathbf{x}}$ satisfies

i.e., when the $\mathbf{I}_{FIM}$ is singular, the Cram$\acute{\mathrm{e}}$r-Rao lower bound will not exist [40, pp. 165], one or more parameters will be unobservable. As in [32], the Fisher information matrix for the models described in (17) can be formulated as:

where $\mathbf{J}$ is the Jacobian of the observation model $\mathbf{g}(\cdot)$ in (13), and its explicit expressions will be given in (22). When $\mathbf{W}>\mathbf{0},$ one has that

Since the first microphone array is viewed as the reference array, its corresponding parameters are all set to zero. The remaining state vectors contain only $\left(N-1\right)$ microphone arrays parameters $\mathbf{x}_{arr}$ and the sound source position $\mathbf{s}^{k}$ at all $K$ time steps. From the definition of the Jacobian matrix [41, pp. 569], we know that $\mathbf{J}\in\mathbb{R}^{g_{1}\times g_{2}}$, where

From (17)-(20), a necessary and sufficient condition for $\mathbf{I}_{FIM}$ to be nonsingular is that $\mathbf{J}$ has full column rank. For $\mathbf{J}$ to be of full column rank, it is necessary that

where $\left\lceil\cdot\right\rceil$ stands for the ceiling operation generating the least integer not less than the number within the operator. We then have the following results.

Proposition 1: The Jacobian $\mathbf{J}$ can be written as

where $\mathbf{L}^{k}=\frac{\partial\mathbf{y}^{k}(\mathbf{x}_{arr},\mathbf{s}^{k})}{{\partial}\mathbf{x}_{arr}}$, $\mathbf{T}^{k}=\frac{\partial\mathbf{y}^{k}(\mathbf{x}_{arr},\mathbf{s}^{k})}{\partial\mathbf{s}^{k}}$ with $\mathbf{y}^{k}(\mathbf{x}_{arr},\mathbf{s}^{k})$ being the inter-array TDOA and DOA observation model at the $k\raisebox{0.0pt}{-}th$ time instant, $k=1,...,K$ (expression of $\mathbf{y}^{k}(\mathbf{x}_{arr},\mathbf{s}^{k})$ can be found in (10); the detailed expressions of $\mathbf{L}^{k}$ and $\mathbf{T}^{k}$ can be found in (47) and (51) in Appendix A, respectively).

Proof. See Appendix A.   

Given the equivalence of full rankness between the FIM and the Jacobian, in the following, we will focus on investigating conditions under which the Jacobian derived in (22) can or can not be of full column rank.

We firstly have the following results regarding the equivalence of full column rank between the Jacobian (22) and matrix $\mathbf{F}$ in (23) which has a much simpler structure.

Proof. The proof is similar to that of [18, Theorem 1] and is skipped here.   

We next present a necessary condition (Theorem 2) and a sufficient condition (Theorem 3) under which matrix $\mathbf{F}$ in (23) is of full column rank.

Proof. See Appendix A.   

It can be seen from Proposition 1 and Theorems 1-3 that observability of the considered identification question is determined both by the configuration of microphone arrays (i.e., the relative transforms, namely, orientation and translation) and the sound source positions. This raises the question of under what conditions on the microphone array configuration and the sound source trajectory, the necessary conditions in Theorem 2 cannot hold. In this section, we will focus on this question and discover some special cases where observability is impossible. Our major result is stated in Theorems 4-5.

Proof. See Appendix A.   

Proof. See Appendix A.   

The observability analysis presented in the above subsections refers to conditions concerning the ground truth value of the sound source trajectories or the configurations of microphone arrays. Hence, the observability analysis is of theoretical interest as it can serve as guidelines when designing microphone array configurations or the sound source trajectories during the calibration process. One can also rely on the results of Section III.C to avoid the unobservable scenarios from a theoretical point of view.

However, during real calibration processes, the measurements contain noises (i.e., the ground truth is not known a prior). Hence, the observability analysis results obtained above are not directly applicable. It is crucial to develop a reliable algorithmic pipeline that can achieve satisfactory convergence and accuracy. This will be discussed in the next section. One should note that the algorithmic pipeline presented in the sequel can also be applied to the nonobservable cases (but the calibration results will be unreliable). This is because, for these scenarios, the noisy measurements do not contain enough information to estimate the unknown parameters. This is also why the analysis in Sections III.A to III.C is valuable, as it suggests avoiding such unobservable situations when designing the microphone array configurations or the sound source trajectories.

Based on the above arguments, to validate the theoretical analysis, we will discuss both observable and unobservable situations in the numerical simulations in Section V. In the experimental results of Section VI, we will only design experiments that correspond to observable cases.

For notational simplicity, in the sequel, we use $\mathbf{d}_{i}^{k}$ and $T_{i}^{k}$ to denote the Gaussian noise corrupted DOA and inter-array TDOA measurements, respectively. We use $\hat{\cdot}$ to represent the estimates of the unknown scalar/vector/matrix parameters. Our proposed initialization procedure is composed of the following main steps: (i) estimation of the sound source position by triangulation; (ii) estimation of the distance between the sound source and microphone arrays using 3D geometry; (iii) estimation of microphone array poses using ICP; (iv) estimation of the asynchronous factors using LLS.

(i) Estimation of the sound source position by triangulation: Without loss of generality, the initial trajectory of the moving sound source in the global frame $\left\{\mathrm{\mathbf{x}}_{arr\_1}\right\}$ is illustrated in Fig. 2(a). Then, from geometry, the initial position of the sound source can be obtained by triangulation and using the first two consecutive DOA measurements as follows:

where $L=\left\|\mathbf{s}_{\Delta}^{1}\right\|_{2}$, i.e., $L$ is the measured distance that the source moves between the first two consecutive moments, $\left\langle\cdot,\cdot\right\rangle$ is the angle of two vectors, and $\hat{d}_{1}^{1}$ is the distance between the first sound source position and the origin. Note that $\mathbf{s}_{\Delta}^{1},\mathbf{d}_{1}^{1},\mathbf{d}_{1}^{2}$ can be obtained from the relative position and DOA measurements, respectively. Once the initial position of the sound source is obtained as above, the sound source positions at different time steps can be estimated:

(ii) Estimation of the distance between the sound source and microphone arrays using 3D geometry: We calculate the distance between each source node and microphone arrays to provide constraints for estimating microphone array poses. One can construct an over-constrained NLS for estimating the distance $\hat{d}_{i}^{k}$ between each source node and microphone arrays by using the law of cosines constraints. To illustrate, as shown in Fig. 2(b), each microphone array and any four source positions $A$, $B$, $C$, and $D$ at the corresponding time instances form a polyhedron (when the four nodes are coplanar, it is tetrahedral, and when the four nodes are on different planes, it forms a five-vertex hexahedral structure). We construct an NLS problem by enforcing the law of cosines for each face of the polyhedron (including the two inner faces). For the scenario shown in Fig. 2(b), denote the estimated squared distance between any two sound source nodes $a$, $b$ among the four source positions $A$, $B$, $C$, and $D$ as:

where $\mathbf{d}_{i}^{a}$ and $\mathbf{d}_{i}^{b}$ are the unit direction vectors of the corresponding sides with length $\hat{d}_{i}^{a}$ and $\hat{d}_{i}^{b}$, respectively. Denote the difference between $\hat{L}_{ab}^{2}$ and $L_{ab}^{2}$ as:

where $m=1,2,\cdots,6$ and $L_{ab}=\left\|\hat{\mathbf{s}}^{a}-\hat{\mathbf{s}}^{b}\right\|_{2}$. Consider a system of six nonlinear equations, given by $F(d_{i}^{A,B,C,D})=\left[F_{1};F_{2};\cdots;F_{6}\right]$. We use ${d}_{i}^{A,B,C,D}$ to collectively denote the distances between the four sound source positions and the $i\raisebox{0.0pt}{-}th$ microphone array, which can be estimated by solving

Note that the nonlinear optimization problem in (27) features a polynomial cost function with a fixed number of unknown parameter dimensions, namely, four edge lengths. However, the batch optimization problem in (15) has a more intricate objective function, incorporating polynomials, exponentials, and trigonometric functions, with $8(N-1)+3K$ optimization variables, where $N$ and $K$ represents the numbers of microphone arrays and time steps, respectively. Hence, in general, the optimization problem in (27) will be much easier to solve (it can be conveniently solved, for instance, using the trust region reflective method [42]) than the entire batch optimization problem in (15). To improve the estimation accuracy of $\hat{d}_{i}^{k}$ at all time instances $K$, we form combinations by selecting any four sound source positions from all time instances, where the $i\raisebox{0.0pt}{-}th$ microphone array-to-$k\raisebox{0.0pt}{-}th$ source line exists within multiple polyhedra. This implies that we can leverage multiple estimation results to achieve greater accuracy. By solving for the edge lengths of each polyhedron and employing the well-known interquartile range (IQR) method [43, pp. 236], we calculate the average value of these same edges in different polyhedra. This average serves as the estimated distance $\hat{d}_{i}^{k}$ between the $i\raisebox{0.0pt}{-}th$ microphone array and the sound source position at the $k\raisebox{0.0pt}{-}th$ time step.

(iii) Estimation of microphone arrays positions and orientations using ICP: Note that the positions of the sound source in the frame $\left\{\mathrm{\mathbf{x}}_{arr\_i}\right\}$ can be estimated as:

We treat the sound source positions as features in each coordinate frame. To find the transformation that optimally aligns the sound source positions with the reference frame is akin to representing the same features in the reference frame. To tackle this challenge, we formulate an NLS problem to minimize the mapping error of sound source positions between $\left\{\mathrm{\mathbf{x}}_{arr\_i}\right\}$ and $\left\{\mathrm{\mathbf{x}}_{arr\_1}\right\}$:

which is conceptually a point-to-point registration problem that can be tackled effectively using ICP [38]. Hence, as in [38], let $\mathbf{p}$ and $\mathbf{p^{\prime}_{i}}$ be the geometric mean of the source position in $\left\{\mathrm{\mathbf{x}}_{arr\_1}\right\}$ and $\left\{\mathrm{\mathbf{x}}_{arr\_i}\right\}$, and they can be computed based on the estimated sound source positions $\mathbf{\hat{s}}_{1}^{k}$ and $\mathbf{\hat{s}}_{i}^{k}$, respectively. The covariance of the sound source trajectory expressed in the two different frames becomes: