Wearable Sensor Localization Considering Mixed Distributed Sources in Health Monitoring Systems

In health monitoring systems, the base station (BS) and the wearable sensors communicate with each other to construct a virtual multiple input and multiple output (VMIMO) system. In real applications, the signal that the BS received is a distributed source because of the scattering, reflection, diffraction and refraction in the propagation path. In this paper, a 2D direction-of-arrival (DOA) estimation algorithm for incoherently-distributed (ID) and coherently-distributed (CD) sources is proposed based on multiple VMIMO systems. ID and CD sources are separated through the second-order blind identification (SOBI) algorithm. The traditional estimating signal parameters via the rotational invariance technique (ESPRIT)-based algorithm is valid only for one-dimensional (1D) DOA estimation for the ID source. By constructing the signal subspace, two rotational invariant relationships are constructed. Then, we extend the ESPRIT to estimate 2D DOAs for ID sources. For DOA estimation of CD sources, two rational invariance relationships are constructed based on the application of generalized steering vectors (GSVs). Then, the ESPRIT-based algorithm is used for estimating the eigenvalues of two rational invariance matrices, which contain the angular parameters. The expressions of azimuth and elevation for ID and CD sources have closed forms, which means that the spectrum peak searching is avoided. Therefore, compared to the traditional 2D DOA estimation algorithms, the proposed algorithm imposes significantly low computational complexity. The intersecting point of two rays, which come from two different directions measured by two uniform rectangle arrays (URA), can be regarded as the location of the biosensor (wearable sensor). Three BSs adopting the smart antenna (SA) technique cooperate with each other to locate the wearable sensors using the angulation positioning method. Simulation results demonstrate the effectiveness of the proposed algorithm.


Introduction
Wearable health-monitoring systems (WHMS) have emerged as an effective way of improving the performance of remote diagnoses and patients' monitoring [1]. As the world population is aging, the healthcare costs will increase, as well. There has been a need to monitor a patient's status while he or she is out of the hospital in his or her personal environment [2]. Motivated by recent technological advances in microelectronics and wireless communication, the wearable sensors of body area networks (BANs), which are parts of the wireless sensor networks (WSN) [3,4], can be used to locate and monitor patient's status. Then, the feedback information about patient's health condition can be provided to caregivers, i.e., medical center, supervising professional physician or even the user himself or herself [5]. The localization of the patient (wearable sensor) is an important parameter for the medical center or caregiver, especially for patients with cardiovascular diseases.
For the 5G link, the peak data rate will likely be in the range of tens of Gbps, which is suitable for real-time communication between wearable sensors and caregivers via a base station (BS), etc. [6]. A general WHMS architecture consisting of the system's functionality and components is depicted in Figure 1. The patients in the hospital or geracomium are relatively intensive. The wearable sensors of patients can be grouped together, and they can communicate with a BS simultaneously on the same resource block. Thus, the BS and the patients construct a virtual multiple input and multiple output (VMIMO) system. At present, the BS mostly adopts a smart antenna (SA) to reduce interference, increase coverage and provide geographic information [7]. The fixed multi-beam antenna and the adaptive array of antennas are two basic types of SA. The former turns the beams on or off while the patient is moving; the latter processes and combines received signals to maximize the signal-to-interference and noise ratio (SINR) [8,9]. The intersecting point of two rays, which come from two different directions measured by two uniform rectangle arrays (URA), can be regarded as the patient's location. Thus, the accuracy of direction-of-arrival (DOA) estimation [10,11] is crucial in both types of SA. In real applications, the effect of angular spread cannot be ignored due to the scattering, reflection, diffraction and refraction of the transmitted signal. If the signal model is regarded as a point source, the DOA estimation performance of the VMIMO system [12,13] will degrade significantly. Generally, a spatially-distributed source consists of hundreds of point sources [14]. The distributed source can be categorized into incoherently-distributed (ID) and coherently-distributed (CD) sources corresponding to rapidly and slowly time-varying channels, respectively. For ID source estimation, several DOA estimation algorithms have been proposed in [15][16][17][18]. A new subspace-based algorithm without eigendecomposition of the covariance matrix has been proposed based on the Capon estimator [15] and the property of the inverse of the covariance matrix [16]. However, the 2D nominal DOA of the user terminal (UT) has to be estimated based on spectrum peak searching, and the computational complexity is relatively high [15,16]. The maximum likelihood (ML) estimator, which has the optimal estimation performance, has been proposed in [17,18]. However, only one distributed source can be estimated in [18]. For CD source estimation, many DOA estimation algorithms have been proposed, as well [19][20][21]. The 2D DOA estimation for two distributed source models (parametric and nonparametric) has been solved based on the MUSIC-based method [19]. A two-step procedure enabling decoupling the estimation of DOA from that of the angular spread has been proposed in [20]. For 2D DOA estimation, the sequential one-dimensional searching (SOS) method has been proposed based on uniform circular arrays (UCA) [21]. However, 1D searching is needed, as well.
In this paper, we adopt the angulation positioning method for patients' (wearable sensors) localizations based on three BSs with the cooperation of multiple VMIMO systems. To the best of our knowledge, there have been few reports about DOA estimation for ID and CD sources jointly. Thus, we propose a 2D DOA estimation algorithm under the coexistence of ID and CD sources. To be more specific, the main contributions of this paper are listed as follows.
(1) Based on the spatiotemporal separation technique, including the second-order blind identification (SOBI) algorithm and the minimum description length (MDL) criterion, the mixed array manifold matrix, including ID and CD sources, is obtained. ID and CD sources can be separated from the amplitude information of the elements of the mixed array manifold matrix. (2) Based on first-order Taylor series expansion of the steering vector and signal subspace algorithm, we extend the ESPRIT to 2D DOA estimation. Three sub-arrays are constructed to form two rotational invariant relationships, and then, the ESPRIT algorithm is used for DOA estimation of ID sources. (3) Based on the application of generalized steering vectors (GSVs), two rotational invariant relationships are constructed, as well. We use ESPRIT to estimate eigenvalues of the two rotational invariant matrices. The DOAs of CD sources contained in the eigenvalues can be obtained finally. (4) Compared to the MUSIC-based method, the proposed algorithm bypasses the spectrum peak searching because of the closed expression of the proposed estimator. Thus, the proposed algorithm has much lower computational complexity, which is suitable for real-time application. (5) The Cramér-Rao bound (CRB) containing ID and CD sources is derived, whereas the known CRB is only valid for the estimation of CD or ID sources.
This paper is organized as follows. The localization scheme and problem formulation are introduced in Section 2. The source separation algorithm is given in Section 3. The proposed DOA estimation algorithm under the coexistence of ID and CD sources is proposed in Section 4. The simulation results are shown and analyzed in Section 5. The conclusions are drawn in Section 6.
Notation: In this paper, the operator (·) T , (·) H and E {·} denote the transpose, conjugate transpose and expectation, respectively. The boldface uppercase letters and boldface lowercase letters denote matrices and column vectors, respectively. The symbol diag{z 1 , z 2 } stands for a diagonal matrix whose diagonal entries are z 1 and z 2 . The symbol blkdiag{Z 1 , Z 2 } stands for a block diagonal matrix, whose diagonal entries are matrices Z 1 and Z 2 .

The Localization Framework
The localization scheme of wearable sensors is shown in Figure 2. Three URAs are installed at the top of three BSs, respectively. There are three resource blocks in Figure 2. The green circles (biosensors) and the BS1 construct VMIMO1; the yellow circles (biosensors) and the BS2 construct VMIMO2; the blue circles (biosensors) and the BS3 construct VMIMO3. These biosensors (wearable sensors) can measure significant physiological parameters, like oxygen saturation, heart rate, body and skin temperature, respiration rate, electrocardiogram, etc. Generally, the obtained data measured by biosensors are delivered to the BS based on a wireless link via a central node, such as a mobile phone in the same resource block; then, the BS transmits the aggregated vital signs to a medical center or a caregiver. The doctor and caregiver can realize the patient's physical condition via a user interface, such as a PC, mobile phone, etc. In theory, the intersecting point of two rays, which come from two different directions measured by two uniform rectangle arrays (URA), can be regarded as the patient's (wearable sensor) location. However, for one biosensor, the BS can only give one estimated result of DOA in one resource block. We need the help of another BS to give another estimated result of DOA. Thus, multiple VMIMO systems (more than two) should cooperate with each other to locate the positions of biosensors. Generally, we use three BSs to locate the positions of biosensors. The third estimated result of BS3 is to guarantee the estimated accuracy and stability. Each biosensor has a unique ID, and the information of all IDs of biosensors has already been stored in the data processing center (DPC). We can match two different rays estimated by two different BSs according to the unique ID. Thus, the position of the biosensor can be located. The DOA estimation algorithm can estimate multiple DOAs simultaneously. Thus, multiple positions of biosensors can be obtained simultaneously. The details can be found in Appendix A1.

The Problem Formulation of DOA Estimation
In VMIMO systems, the rapidly time-varying and slowly time-varying channels corresponding to ID and CD sources both exist. Thus, the mixed distributed sources, which are a more general case, should be considered. In this subsection, the received data model of mixed distributed sources is constructed, and it is adopted in the separation and 2D DOA estimation of ID and CD sources.
Assume that K narrowband uncorrelated distributed sources impinge on a URA with the number of elements M = M x M y , where M x and M y are the number of elements respectively placed along x and y directions, as shown in Figure 3, i.e., a BS of a VMIMO system with M elements communicates with K UTs simultaneously using SA technique [8]. The inter-element spacing along the x and y direction is half the wavelength. In a cellular mobile communication system, due to the reflection and scattering in the multipath propagation, the angular spread cannot be ignored. The point source model is replaced by a distributed source model. Thus, the M × 1 received data of a BS can be expressed as [14]: (1) where t = 1, . . . , T is a discrete series corresponding to the sampling points from one to T, and T is the snapshot number. The matrix A = [A ID , A CD ] = h 1 , . . . , h K ID , h K ID +1 , . . . , h K ∈ C M×K is the array manifold matrix with the source (biosensor) number K = K ID + K CD , where K ID is the source number of ID sources and K CD is the number of CD sources; h k is the steering vector of the k-th signal source.
, . . . , and: are the array manifold matrices corresponding to ID and CD sources, respectively. The signals vector s (t) = s T ID (t) , s T CD (t) T ∈ C K×1 is the signal transmitted by biosensors.
s ID (t) = s 1 (t) , . . . , s K ID (t) ∈ C K ID ×1 and s CD (t) = s K ID +1 (t) , . . . , s K (t) ∈ C K CD ×1 are the signal vectors corresponding to ID and CD sources, respectively. γ k,j (t), θ k,j (t) and φ k,j (t) are the complex gain, azimuth and elevation of the j-th path for the k-th ID source, respectively, which satisfy 0 ≤ θ k,j (t) ≤ π, 0 ≤ φ k,j (t) ≤ π/2, k = 1, . . . , K ID , j = 1, . . . , L k . L k is the number of multipaths of the k-th ID source. The steering vector a θ k,j (t) , φ k,j (t) ∈ C M×1 is the array response corresponding to θ k,j (t) and φ k,j (t). The m-th element of the array is given by: where ω = 2πd λ, d and λ are the inter-sensor spacing and the wavelength, respectively. The azimuth and elevation can be respectively written in another form as: whereθ k andφ k are nominal azimuth and elevation of the k-th ID source, respectively,θ k,j (t) and φ k,j (t) are referred to as the angular spreads, which correspond to random angular deviations with zero mean and standard deviations σ θ k and σ φ k , respectively.
In addition, some assumptions are considered throughout this paper.
(1) ID sources are uncorrelated with CD sources, and they are uncorrelated with the noise.
(2) The complex gains γ k,j (t), k = 1, . . . , K ID , j = 1, . . . , L k , t = 1, . . . , T, are independent and identically distributed (i.i.d.) complex-valued zero-mean random variables: The noise n (t), t = 1, . . . , T is an i.i.d. random variable both in temporally-and spatially-complex-valued circularly symmetric zero-mean Gaussian variables, whose covariance matrix is given by: where E n (t) n H (t) is the covariance matrix. (4) The number of multipaths L k is sufficiently large, k = 1, . . . , K ID . (5) The BS equipped with M elements is much larger than the number of biosensors K [22]. (6) The signal powers S IDk = |s IDk (t)| 2 , k = 1, . . . , K ID of ID sources are known as a prior; the cellular area of the same resource block is not large; thus the power loss can be ignored.
Finally, we emphasize that the task is to estimate 2D nominal DOAθ k andφ k for ID and CD sources, respectively, k = 1, . . . , K, based on the received snapshot data x (t), t = 1, . . . , T.

The Distributed Source Classification Based on SOBI
In this section, we will separate the ID and CD sources from the received signals with the help of the SOBI algorithm [23,24].
The covariance matrix of the k-th ID source can be modeled as [25]: where B (Φ k ) is the real-valued symmetric Toeplitz matrix (the derivation is given in Appendix A2). The channel of the k-th CD source is defined as [26]: where g (µ k ) is a real-valued term. The SOBI algorithm is a blind signal processing technique based on a second-order statistic characteristic. For the k-th source s k (t), the covariance matrix with nonzero time lag e can be expressed as [23]: Then, the correlation covariance of x (t) with time lag e is given by: where . Note that R (e l ), l = 1, . . . , L contain the identical array manifold with respect to different e ∈ {e 1 , . . . , e L }. It can be proven that if the difference among the diagonal entries of D e is obvious and A is full column rank, then multiple R (e l ) can be joint diagonalized using A or its permutation matrix. Then, the array manifold matrix A can be accurately estimated using R (e l ) with different time lag e. Assume that there are two correlation covariance matrices R (e 1 ) and R (e 2 ); we have: Because D e 1 D −1 e 2 is a diagonal matrix, the unique estimate of array manifold matrix A can be obtained by taking the eigenvalue decomposition (EVD) of D e 1 D −1 e 2 . In real applications, R (e) can be estimated by the limited sample, which is given by: Then, the SOBI algorithm is summarized as Algorithm 1 [24].
Algorithm 1 Implementation of the SOBI Algorithm.
As shown in Equation (9), for the ID source, the amplitude of the entries of the covariance matrix are getting smaller as angular spreads increase; as shown in Equation (10), the amplitude of entries of generalized steering vectors are getting smaller as angular spreads increase.
For distributed sources, the amplitude of entries of array manifold matrix A is usually less than one. However, if n (source number (SN)) is larger than the actual source number, after the signals' separation, the amplitude of entries of array manifold matrix A is usually lager than one based on experiment. When SN is equal to the actual source number, the amplitude of entries of array manifold matrix A is usually not larger than one. Usually, only the ID source may cause this phenomenon. If the ID sources exist, SN is larger than the actual source number based on the minimum description length (MDL) criterion [27]. Thus, the amplitude of any entries of array manifold matrix A is larger than one, and SN > 1; we can distinguish it as an ID source. Then, SN is reduced by one, and the SOBI is applied again. This procedure can be done iteratively until all ID sources are separated. Based on the discussion above and Algorithm 1 , the separation algorithm for ID and CD sources is summarized as Algorithm 2.

Algorithm 2
The Separation Algorithm for ID and CD Sources.

DOA Estimation for ID Source
For the VMIMO system, the existing subspace-based and covariance matching-based algorithms are sophisticated for 2D DOA estimation, because of their tremendous computational complexity of multidimensional searching. In this section, based on the array manifold matrix estimated in the previous section, an ESPRIT-based algorithm is proposed to deal with the problem of 2D DOA estimation with low computational complexity.
Based on first-order Taylor series expansion of a θ k,j (t) , φ k,j (t) , the steering vector in Equation (4) can be approximated as: where the terms after the second term are ignored. If standard deviations σ θ k and σ φ k are sufficiently small, the first-order Taylor series expansion is almost equal to a θ k,j (t) , φ k,j (t) . Then, the received data of ID source can be rewritten as: where c k, It can be seen from Equation (16) that the relationship between the received data of the ID source and a θ k,j (t) , φ k,j (t) is linear; so is its partial derivatives. Thus, the received data of ID sources can be rewritten as: where: Then, the new array manifold matrix is expressed as: the elements of c(t) are functions of the incident signal, the path gains and angular deviations. It can be known that A is only determined by the nominal DOAs,θ k andφ k , k = 1, . . . , K ID . Thus, we can obtain the DOA of ID sources from A, and the covariance matrix of c (t) can be expressed as: It is a diagonal matrix with S IDk = |s Based on Algorithm 2, the steering vectors of ID sources are separated. The estimator ofÂ ID is obtained. Then, the estimator of snapshot datax ID (t) for ID sources can be obtained by: Then, the covariance matrix ofÂ ID is expressed as: whereR IDs = diag S ID1 , . . . , S IDK ID . It can be known thatR IDs is a normal and positive diagonal matrix. In general, B ID is a full column rank matrix; the EVD ofR ID is given by: where Q s ∈ C 3K ID ×3K ID is a diagonal matrix with the entries of 3K ID large eigenvalues ofR ID . The remaining M − 3K ID small eigenvalues ofR ID equal to σ 2 n . Their corresponding subspace are signal subspace U s ∈ C M×3K ID and noise subspace U n ∈ C M×(M−3K ID ) , respectively.Q s = Q s − σ 2 n I 3K ∈ R 3K ID ×3K ID . Based on Equations (22) and (23), we can know that: SinceQ s is a diagonal matrix, as well, there exists a nonsingular matrix T ∈ C 3K ID ×3K ID , which satisfies: The linear relationship between B ID and U s will be used to estimate the nominal DOAs.
In order to obtain 2D angle estimation, azimuth and elevation, two rotation-invariant relationships have to be constructed [28]. As shown in Figure 4, the whole URA is divided into three sub-arrays. Although more sub-arrays can be divided, the computation complexity increases rapidly in the VMIMO system. In order to reduce the computational complexity, only two rotation-invariant relationships are adopted. This is different from the conventional ESPRIT algorithm; this rotational-invariant relationship contains a θ k ,φ k and its partial derivatives. The array manifold matrices of sub-arrays have the same form as that of URA. The steering vector of the l-th sub-array is defined as a l θ k ,φ k ∈ CM ×1 , l = 1, 2, 3 andM = (M x − 1) M y − 1 . In order to estimate the subspaces U ID1 , U ID2 and U ID3 of three sub-arrays, we need to construct a selection matrix as follows: where d 1 = 0, d 2 = 1 and d 3 = M x . In (26), the floor operator makes m − 1/M x − 1 = n, ∀n = 0, . . . , M y − 2 when m = n (M x − 1) + 1, . . . , (n + 1) (M x − 1). It can be seen that for the m-th row of J l , only the (m + m − 1/M x − 1 + d l )-th entry is one, and the other entries are zeros. Thus, J l assigns the (m + m − 1/M x − 1 + d l )-th of U s into the subspaces U ID1 , U ID2 and U ID3 belonging to three sub-arrays, respectively. This coincides with the relationships between the sub-arrays and the URA. Then, the estimators of subspaces U ID1 , U ID2 and U ID3 can be expressed as: Proposition 1. The subspaces U ID1 , U ID2 and U ID3 belonging to three sub-arrays have the relationships as follows: where: Proof. See Appendix A3.
Thus, eigenvalues of V 1 and V 2 are diagonal entries of W 2,1 and W 3,1 , respectively. However, W 2,1 and W 3,1 are not diagonal matrices, but upper triangular matrices. Then, the conventional ESPRIT algorithm cannot be used directly. Therefore, in order to estimate the diagonal elements of W 2,1 and W 3,1 , V 1 and V 2 need to be estimated from the three subspaces U ID1 , U ID2 and U ID3 . According to Equation (28), the estimators of V 1 and V 2 can be obtained by employing the well-known total least-squares (TLS) method. Referring to the similar derivation in [29], we construct a new matrix: and the rank ofŨ ID1 is 6K ID . The EVD ofŨ ID1 is expressed as: where U x ∈ C 6K ID ×6K ID and P x ∈ C 6K ID ×6K ID are eigenvectors and eigenvalues ofŨ ID1 , respectively. Then, U x ∈ C 6K ID ×6K ID can be partitioned into four block matrices as: where U xab ∈ C 3K ID ×3K ID , a, b = 1, 2. The estimator of V 1 is expressed as: In order to estimate nominal DOAs, the EVD ofV 1 is expressed as: where T 1 ∈ C 3K ID ×3K ID and P 1 ∈ C 3K ID ×3K ID are eigenvectors and eigenvalues ofV 1 , respectively. Similar to the process of estimatingV 1 , the estimator of V 2 is expressed as: where T 2 ∈ C 3K ID ×3K ID and P 2 ∈ C 3K ID ×3K ID are the eigenvectors and eigenvalues ofV 2 , respectively. However, the EVDs ofV 1 andV 2 are completed separately. The eigenvalues between P 1 and P 2 have to be matched. Denote the k-th diagonal entry of P d as P dk , d = 1, 2. The eigenvectors of the identical source are strongly correlated; thus, we can construct the sequencing matrix G = T H 1 T 2 to match P 1k and P 2i . According to the coordinate of the maximal entry in the columns (or rows) of matrix G, we adjust the order of eigenvectors. Then, the parameter pairing is completed. The estimators of [W 2,1 ] k+(l−1)K,k+(l−1)K and [W 3,1 ] k+(l−1)K,k+(l−1)K , which are the eigenvalues P 1,3(k−1)+l and P 2,3(k−1)+l , l = 1, 2, 3 ofV 1 andV 2 , respectively, are obtained. According to Equations (34) and (35), we have: The estimators of the nominal 2D DOAsθ k andφ k for ID sources are given by: where k = 1, . . . , K ID . Thus, the 2D DOA estimation for ID sources is completed. The pseudo-code is summarized as Algorithm 3.

Algorithm 3
The DOA Estimation Algorithm for ID Sources.
1: Estimate the covariance matrixR ID according to Equation (25); 2: Take the EVD ofR ID ; Q s is a diagonal matrix with the entries of 3K ID ; large eigenvalues ofR ID , U s ∈ C M×3K ID are the corresponding signal subspace; 3: Divide the URA into three sub-arrays; calculate the selection matrix according to Equations (A5)-(A7); construct two different rotational invariant relationships according to Equation (46); 4: Construct new matrixŨ ID1 ; perform the EVD of Equation (31), partitioning the matrix U x into four blocks; the eigenvalues ofV 1 andV 2 can be obtained according to Equation (33); 5: Estimate 2D DOA for ID sources according to Equations (38) and (39).

DOA Estimation for the CD Source
We first define threeM × M selection matrices: The received data of the three sub-arrays for CD sources can be expressed as: where n 1 , n 2 , n 3 are additive noise vectors and s k (t) is a random process of the k-th signal source.
The generalized steering vectors [30] of the k-th CD source belonging to three sub-arrays can be respectively expressed as: where k = 1, . . . , K CD .

Proposition 2. For a small angular extension, there is an approximate rotational invariant relationships between
where k = 1, . . . , K CD .
Then, we can rewrite them in matrix form as: where: Based on Algorithm 2, the array manifold matrix for CD sources, which is denoted asÂ CD , can be estimated. By multiplying the selection matrices in Equation (40), the array manifold matrices of three sub-arrays can be respectively expressed as: According to Equation (50), two rotational invariant relationships can be respectively expressed as: The k-th diagonal entries ofŴ x andŴ y are defined as ε xk and ε yk , respectively. The estimators of the nominal 2D DOA θ k and φ k for CD sources are respectively given by: where k = 1, . . . , K CD . Thus, the 2D DOA estimation for CD sources are completed. The pseudo-code is summarized as Algorithm 4.

Algorithm 4
The DOA Estimation Algorithm for CD Sources.
1: Define threeM × M selection matrices according to Equation (40); 2: Based on Algorithm 2, the array manifold matrixÂ CD for CD sources is estimated; 3: Divide the URA into three sub-arrays; the array manifold matrix of these can be obtained by multiplying the corresponding selection matrix based on Equation (50); 4: Calculate the two rotational invariant relationships,Ŵ x andŴ y , according to Equation (51); 5: Estimate 2D DOA for CD sources according to Equations (52) and (53).
Thus, the 2D DOA estimation for ID and CD sources can be summarized as follows.

Algorithm 5
The DOA Estimation Algorithm for Mixed Sources.
1: Separate the ID and CD sources based on SOBI algorithms Algorithm 1 and Algorithm 2; 2: Estimate the 2D DOA for ID sources according to Algorithm 3; 3: Estimate the 2D DOA for CD sources according to Algorithm 4; Remark: The proposed algorithm can estimate 2D DOAs of ID or CD sources or joint ID and CD sources. If only ID or CD sources exist in the incident signals, the separation of ID and CD sources can be avoided. It is worth noting that the part of the proposed approach for 2D DOA estimation of ID sources can also be applied to other scenarios by exploiting the rotational invariance property of the array's structure, such as uniform linear arrays (ULAs) and uniform cylindrical arrays (UCyA) [13]. The part of the proposed approach for 2D-DOA estimation of the CD source has a similar characteristic, i.e., the rotational invariance property of the antenna array's structure is exploited. Thus, the proposed approach can be applied to other scenarios by exploiting the rotational invariance property of the array's structure, i.e., there exist three sub-arrays (the array's structure can be arbitrary), which can construct two different rotational invariance relationships; the 2D DOA estimation for ID and CD sources can be achieved with little modification of the proposed algorithm.

Simulation Results
In this section, we will illustrate the performance of the proposed algorithm. The MUSIC-based method is used after source separation. The PM [30], MUSIC-based algorithm [31] and the CRB are considered for the performance comparison. In the simulation, the proposed algorithm is called the ESPRIT-based algorithm.
The simulation condition is introduced as follows. The snapshot number is T = 1000. The group number of the sampling covariance matrix is L = 10 in the SOBI algorithm. The numbers of the elements are 100, i.e., there are M x = 10 and M y = 10 in the x-axis and y-axis, respectively. The number of biosensors is K = 5, i.e., one ID source and one CD source. The number of multipaths for the ID source is N = 50. The nominal azimuth and elevation for the ID source are 20 • , 30 • , 30 • , 50 • , respectively. The azimuth and elevation angular spreads for ID source are both equal to For the CD source, the nominal azimuth and elevation for the CD source are 60 • , 40 • , 70 • , 50 • and 80 • , 40 • , respectively. The azimuth and elevation angular spreads for the CD source are both equal to σ ϑ 2 = σ φ 2 = 2 • . The ID source and the CD source both satisfy the Gaussian-shaped distribution. The noise is the additive Gaussian white noise, and it is not correlated with signals. Five hundred Monte Carlo trials are taken in the simulation. The average root mean square error (RMSE) is defined as: where κ stands forθ k (m) andφ k (m), which are the estimates of θ k and φ k of the m-th Monte Carlo trials, respectively.
As shown in Figure 5, the RMSE of the azimuth estimation for the ID source versus SNR is depicted. It can be seen that the RMSE of the ESPRIT algorithm is larger than that of the MUSIC algorithm, and the CRB has a low SNR. When the SNR increases, the RMSE of the ESPRIT algorithm decrease rapidly. However, the MUSIC algorithm varies slowly as the SNR changes. The PM has the largest RMSE of all of the algorithms. This is mainly caused by the signal and noise subspaces not being orthogonal. The subspace method has the lowest RMSE. However, the computational complexity is tremendous. The RMSE of the elevation estimation for the ID source versus SNR is depicted in Figure 6. The RMSE of the azimuth and elevation estimation for the CD source versus SNR is depicted in Figures 7  and 8, respectively. The curve trend of the PM, the ESPRIT algorithm, the MUSIC and subspace algorithm and the CRB in Figures 6-8 are similar. It can be seen that when SNR is low, the RMSE of the ESPRIT algorithm is larger than that of MUSIC, the subspace algorithms and the CRB. The RMSE of the MUSIC and subspace algorithms become smaller as the SNR increases, but the reduced amplitude is not obvious.   As shown in Figure 9, the RMSE of the azimuth estimation for ID source versus snapshot number is depicted. It can be seen that the RMSE of the ESPRIT algorithm is larger than that of the MUSIC and subspace algorithms.
The RMSE of the elevation estimation for the ID source versus snapshot number from UTs is depicted in Figure 10. The RMSE of the azimuth and elevation estimation for the CD source versus snapshot number is depicted in Figures 11 and 12, respectively. The curve trends of the PM, ESPRIT, MUSIC and subspace algorithms and the CRB in Figures 10-12 are similar, as well. It can be seen that the RMSE of the ESPRIT algorithm is larger than that of the MUSIC and subspace algorithms and the CRB.  We evaluate the averaged CPU times of PM, MUSIC, the subspace algorithm and ESPRIT in the following experiment. The simulation condition is the same as Figure 1. The snapshot number is fixed at 2000, and the SNR is fixed at 3 dB. The experiment is carried out in MATLAB v.8.3.0 on a PC with a Windows 7 system and a 3-GHz CPU. Table 1 presents the the averaged CPU times of PM, MUSIC, the subspace algorithm and ESPRIT. The MUSIC and subspace method is the most time consuming, since the spectrum peak searching is needed. PM costs the least time of all, since the spectrum peak searching is not needed. Thus, ESPRIT is chosen as the DOA estimation algorithm for the biosensor's localization, since it performs well in the estimated accuracy and resolution probability, and it does not cost much time to execute the algorithm. It should be noted that the curve of the MUSIC algorithm does not change much in the above figure. This phenomenon is mainly caused by the MUSIC algorithm being directly used by multiplying the estimated array manifold matrix without any form transformation.
The proposed algorithm does not perform well when the SNR is low. However, as we all know, the BS of VMIMO is equipped with a large number of elements. We can deal with this problem by increasing the number of elements equipped in the BS as a trade-off. The proposed algorithm has a larger RMSE than that of the MUSIC algorithm, but the proposed algorithm has much lower computational complexity compared to the MUSIC algorithm.

Conclusions
For the location of the biosensor (wearable sensor) in health monitoring systems, first, a localization scheme based on a multiple VMIMO system is constructed. The intersecting point of two rays, which come from two different directions measured by two uniform rectangle arrays (URA), can be regarded as the location of the biosensor (wearable sensor). Then, a 2D DOA estimation algorithm under the coexistence of ID and CD sources is proposed. The computational complexity of the proposed algorithm is much lower than that of other multidimensional parameter searching algorithms, such as the MUSIC algorithm. The ID and CD sources are processed separately. Three sub-arrays are selected to construct rotational-invariant matrices. The simulation result confirms that the proposed algorithm outperforms the MUSIC-based algorithm when the SNR is larger than a certain threshold value. In future work, we will focus on the effective separation method between the ID and CD sources with low computational complexity. In addition, the angular spreads' estimation should be considered in our future work.

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

A1. The Localization Method Based on Multiple VMIMO
A similar scheme has been proposed in our previous work [32,33]. It is used for the localization of patients in the emergency healthcare system and partial-discharge source diagnosis and localization in an industrial high-voltage insulation system. The effectiveness of the proposed scheme in [32,33] has been verified. The real test has been done for partial-discharge source localization [33]. In order to have an intuitive impression, we give the localization method of the proposed scheme as follows.
As shown in Figure A1, three URAs cooperate with each other to locate the biosensors' positions. The azimuths θ 1 , φ 1 and θ 2 , φ 2 can be estimated by using the algorithm proposed in Section 4. The coordinates of reference Points A, B and C corresponding to VMIMO1(URA1), VMIMO2(URA2) and VMIMO3(URA3) are (x A , 0, z), (0, 0, z) and (0, y C , 0), respectively. Based on the measurements of URA1 and URA2, we have: Figure A1. The localization of the biosensor.
Based on the measurements of URA2 and URA3, we would have another estimated result. These two results can be used to improve the estimated stability of the proposed scheme.

A2. The Derivation of Equation (9)
Assume the k-th ID source impinges on the array; the covariance matrix of received data can be modeled as: where q k (θ i , φ i ) , θ j , φ j ; σ θ , σ φ is the kernel function of angular auto-correlation. p (α) is the von Mises angular function, whose definition is: where J 0 (α) is first kind zero-order Bessel function, whose definition is: If σ θ > 0 and σ φ > 0, we have: where p I θ, φ, σ θ , σ φ is the angular power density of the ID source and B (Φ k ) is a real-valued symmetric Toeplitz matrix.

A4. Cramér-Rao Bound
In this section, the CRB of 2D DOA estimation, including ID and CD sources, is derived. Generally, the received data of URA can be expressed as: where x (t), a complex Gaussian vector with zero mean, h k (t) is the k-th channel parameter, H ∈ C M×K is the channel parameter matrix and s (t) ∈ C K×1 is the signal vector. The covariance matrix of x (t) is expressed as: For T statistically-independent observations of x (t), the unknown DOAs vector of R x is defined as: where θ = θ 1 , . . . , θ K ID , θ K ID +1 , . . . , θ K T and φ = φ 1 , . . . , φ K ID , φ K ID +1 , . . . , φ K T . The general expression of the (m, n)-th element in the Fisher information matrix (FIM) can be expressed as: Since only the m-th column ofḢ ϑ m is nonzero, then it can be represented asḢ where the m-th column of the identity matrix is defined as γ m K . H θ is the derivative matrix of the channel parameter matrix, which is expressed as: Then, Equation (A30) can be rewritten as: The FIM corresponding to θ can be expressed as: