# 3D Positioning Algorithm Design for RIS-aided mmWave Systems

Tuo Wu, Cunhua Pan, Yijin Pan, Sheng Hong, Hong Ren, Maged Elkashlan, Feng Shu
and Jiangzhou Wang, _Fellow, IEEE_ (Corresponding author: Cunhua Pan). T. Wu
and M. Elkashlan are with the School of Electronic Engineering and Computer
Science at Queen Mary University of London, London E1 4NS, U.K.
(Email:{tuo.wu, maged.elkashlan}@qmul.ac.uk). C. Pan, Y. Pan and H. Ren are
with the National Mobile Communications Research Laboratory, Southeast
University, Nanjing 210096, China. (Email: {cpan, panyj, hren}@seu.edu.cn). S.
Hong is with Information Engineering School of Nanchang University, Nanchang
330031, China. (Email: shenghong@ncu.edu.cn). F. Shu is with the School of
Electronic and Optical Engineering, Nanjing University of Science and
Technology, Nanjing 210094, China, and also with the School of Information and
Communication Engineering, Hainan University, Haikou 570228, China.(E-mail:
shufeng0101@163.com). J. Wang is with the School of Engineering, University of
Kent, UK. (e-mail: J.Z.Wang@kent.ac.uk).

###### Abstract

In this paper, we investigate a three-dimensional (3D) positioning algorithm
for a millimeter wave (mmWave) system, where the reconfigurable intelligent
surfaces (RIS) are leveraged to enhance the positioning performance of mobile
users (MUs). We propose a two-stage weight least square (TSWLS) algorithm to
obtain the closed-form solution of the MU’s position. In the first stage, we
construct the pseudolinear equations based on the angle of arrival (AOA) and
the time difference of arrival (TDOA) estimation at the RISs, then we obtain a
preliminary estimation by solving these equations using the iterative weight
least square (WLS) method. Based on the preliminary estimation in the first
stage, a new set of pseudolinear equations are obtained, and a finer
estimation is obtained by solving the equations using the WLS method in the
second stage. By combining the estimation of both stages, the final estimation
of the MU’s position is obtained. Further, we study the theoretical bias of
the proposed algorithm by considering the estimation error in both stages.
Simulation results demonstrate the superiority of the proposed positioning
algorithm. Furthermore, it is also shown that the proposed algorithm still
have good positioning performance with low SNR.

###### Index Terms:

Reconfigurable intelligent surface (RIS), intelligent reflecting surface
(IRS), positioning, sensing, millimeter wave (mmWave).

##  I Introduction

In the upcoming sixth generation (6G) wireless network, a high positioning
accuracy is urgently required by various position-related services [1], such
as smart factories [2], autonomous vessels [3], automated vehicles [4] and
mobile user (MU) sensing [5]. Hence, a wealth of wireless positioning
algorithms that are widely deployed in the fifth generation (5G) millimeter
wave (mmWave) network, are now being re-imagined in the 6G network that aims
for a higher positioning acuracy [5].

By deploying positioning algorithms in mmWave systems [6, 7, 8, 9], the
positioning accuracy can be improved as the positioning parameters (e.g., time
of arrival (TOA), time difference of arrival (TDOA) and angle of arrival
(AOA)) can be estimated more accurately with narrow beams. However, since
wireless signals in mmWave band are more vulnerable to blockages [7], the
reliability of high-precision positioning capability faces great challenges
[10]. Reconfigurable intelligent surface (RIS) is a promising technique that
is capable of constructing new line-of-sight (LoS) communication links [11,
12, 13, 14, 15], which can be deployed in the existing radio positioning
systems with various advantages. First, when deploying a wireless positioning
reference, RIS is more preferred than the base station (BS) in terms of
hardware cost [16]. Second, as a passive panel, RIS provides reliable and
high-precision estimation performance with low energy consumption [17].
Finally, the large size of the RIS panel enables high accuracy radio
positioning parameter estimation [18]. Therefore, it is meaningful to develop
positioning algorithms for RIS-aided mmWave systems.

Currently, RIS-aided positioning is extensively studied, including performance
analysis [19, 20, 21, 22, 23] and algorithm design [24, 25, 26, 27]. For
performance analysis, the authors of [19] derived the theoretical positioning
performance bounds, such as Cram$\rm{\acute{e}}$r-Rao lower bound (CRLB), to
provide analytical performance validation. In [20], Liu _et al._ derived the
CRLB of a multiple-RIS-aided mmWave positioning system. In [21], the authors
provided CRLB expressions and the effective achievable data rate (EADR), in
which the phase shifts of the RIS were optimized to maximize the EADR. The
authors of [22] derived the CRLB in terms of the position parameters of the
target node, and provided the theoretical relationship between the optimal
power and the phase shifts of the RIS. To assess the performance of
synchronous and asynchronous signaling schemes, Alouini _et al._ derived the
CRLB in [23]. For algorithm design, [24] proposed a received signal strength
(RSS) based positioning scheme enabled by an RIS and designed an iterative
algorithm to obtain the optimal phase shifts. An optimization method was
proposed in [25] to design the phase shift for RIS-aided positioning and
communication systems. In [26], a reduced-complexity maximum likelihood based
estimation procedure was designed to jointly recover the user position and the
synchronization offset.

However, the above-mentioned contributions mainly focused on performance
analysis, such as CRLB, and did not provide specific positioning algorithms.
For the positioning algorithm design, the existing works mainly focused on
deriving the optimal phase shifts rather than the positioning estimation.
Therefore, practical positioning algorithms based on radio positioning
parameters at the RISs and estimation error analysis are needed. In addition,
to evaluate the positioning performance, the existing contributions mainly
utilized the CRLB and the root mean square error (RMSE). However, the
estimated positioning parameters generally are not biased. As shown in [28],
the bias is considerably high when the signal-noise ratio (SNR) is low. Thus,
it is necessary to conduct bias analysis.

Against the above background, we investigate an mmWave positioning system
aided by multiple RISs in this paper and we design a two-stage weight least
square (TSWLS) algorithm to derive the closed-form solution of the position.
Then, bias analysis is provided for the proposed positioning algorithm. Our
contributions are summarized as follows:

  * 1)

We have proposed a three-dimensional (3D) positioning algorithm for mmWave
systems aided by multiple RISs. As multiple RISs are deployed in the
positioning system, we utilize AOAs and TDOAs at the RISs for positioning
algorithm design. Pseudolinear equations involving the AOAs and TDOAs are
formulated based on the geometry relationship and the estimation error.

  * 2)

The TSWLS algorithm is derived to obtain the closed-form expression of the
MU’s 3D position. To be specific, in the first stage, we construct the
pseudolinear equations using the AOAs and TDOAs parameters at the RISs. Then,
preliminary estimation is derived by solving these pseudolinear equations
using the iterative WLS method. In the second stage, based on the preliminary
estimation in the first stage, a new set of pseudolinear equations are
formulated and the finer estimation is obtained by solving the equations.
Finally, using the estimation of both stages, the final estimation of the MU’s
position is obtained.

  * 3)

To evaluate the performance of the proposed positioning algorithm, we
investigate the theoretical estimation bias. More specifically, by decomposing
the matrices containing the estimation information into true values and
estimation error, the estimation error of both stages can be further derived,
and the theoretical bias of both stages is derived accordingly. Finally, the
bias of the final solution is derived using the bias of both stages.

  * 4)

Simulation results are provided to evaluate the performance of the proposed
positioning algorithm. First, it is shown that the proposed algorithm still
have good positioning performance with low SNR. Then, the estimated 3D
position estimated with the proposed algorithm is more accurate than that
derived by using the geometry relationship directly. Finally, the derived bias
results are in good agreement with the simulation results, which confirms the
accuracy of the derived results.

The remainder of the paper is organized as follows. The system model for the
two-step positioning aided by multiple RISs is described in Section II. The
closed-form solution of the MU’s position is given in Section III. Section IV
derives the bias analysis of the proposed positioning scheme. Simulation
results are given in Section V. Section VI concludes this work.

Figure 1: RIS-aided positioning system model.

##  II System Model

Consider an RIS-aided positioning system, where an MU sends pilot signals to
the BS to locate the MU with the assistance of multiple RISs. The BS is
equipped with a uniform linear array (ULA) of $N_{b}$ antennas, and the MU is
equipped with a single antenna. Moreover, there are $M$ RISs, and each RIS is
a uniform planar array (UPA).

The BS is placed parallel to the x-axis with the center located at
$\mathbf{p}=[x_{p},y_{p},z_{p}]^{T}$. The $i$th RIS is placed parallel to the
y-o-z plane with its center located at
$\mathbf{s}_{i}=[x_{i},y_{i},z_{i}]^{T}$, $i=1,2,\cdots,M$. The UPA-based RIS
has $N_{y,z}=N_{y}\times N_{z}$ reflecting elements, where $N_{y}$ and $N_{z}$
denote the numbers of reflecting elements along the y-axis and z-axis,
respectively. The true position of the MU is
$\mathbf{q}=[x_{q},y_{q},z_{q}]^{T}$ and is assumed to be placed parallel to
the x-o-y plane. The estimated location of the MU is
$\mathbf{\hat{q}}=[\hat{x}_{q},\hat{y}_{q},\hat{z}_{q}]^{T}$. Generally, once
the RISs and BS have been deployed, the coordinates $\mathbf{s}_{i}$ and
$\mathbf{p}$ are known and invariant. In order to locate the MU, we need to
obtain the estimated $\mathbf{\hat{q}}$.

In this paper, multiple estimated radio positioning parameters (AOAs and
TDOAs) are utilized for the wireless positioning. As shown in Fig. 1, there
are two kinds of line-of-sight (LoS) communication links containing the radio
positioning parameters, which are the direct links from the MU to the BS, and
the reflecting links through the RISs.

###  II-A AOA information

For the MU-RIS-BS link, the available angle information includes the AOAs
(azimuth and elevation angles) at the RISs, the AODs (azimuth and elevation
angles) at the RISs, and the AOAs at the BS. For the MU-BS link, the angle
information for positioning is the AOA at the BS. However, since we can only
estimate the azimuths of the AOAs at the BS due to the assumption of ULAs, the
3D coordinates of the MU cannot be estimated without the elevation of the
AOAs. Therefore, we use the AOAs at the RISs for the positioning estimation in
this paper.

In Fig. 1, the azimuth AOA $\theta_{MR,i}$ at the RIS is the angle between the
projection of the wave vector on the x-o-y plane and the y-axis as

$\theta_{MR,i}=\arctan\frac{x_{q}-x_{i}}{y_{q}-y_{i}},$ (1)

where $\theta_{MR,i}\in(0,\pi)$.

The elevation AoA $\phi_{MR,i}$ at the RIS of the MU-RIS link is

$\phi_{MR,i}=\arctan\frac{z_{q}-z_{i}}{\sin\theta_{MR,i}(x_{q}-x_{i})+\cos\theta_{MR,i}(y_{q}-y_{i})},$
(2)

where $\phi_{MR,i}\in(-\pi/2,\pi/2)$.

###  II-B TDOA information

The TDOA information (time delays) can be extracted from the received pilot
signals. Then, the difference of link distances can be readily obtained by
multiplying the light speed with TDOAs. Thus, in the following, we focus on
the difference of distances, as it is equivalent to TDOAs.

Let $R_{BU}$ denote the true distance of the MU-BS link, which can be
calculated as

$R_{BU}=\sqrt{(x_{q}-x_{p})^{2}+(y_{q}-y_{p})^{2}+(z_{q}-z_{p})^{2}}.$ (3)

Additionally, the true distance from the MU to the $i$th RIS can be expressed
as

$R_{RU,i}=\sqrt{(x_{q}-x_{i})^{2}+(y_{q}-y_{i})^{2}+(z_{q}-z_{i})^{2}}.$ (4)

Using $R_{BU}$ as the reference distance, the true distance difference between
$R_{{RU,i}}$ and $R_{{BU}}$ can be written as

$R_{{B,i}}=R_{{RU,i}}-R_{{BU}}.$ (5)

The true distance difference $R_{B,i}$ will be used in the following
positioning algorithm design.

###  II-C Radio positioning parameters estimation

The estimated radio positioning parameters can be modeled as

$\hat{\theta}_{MR,i}={\theta}_{MR,i}+n_{i},\quad\hat{\phi}_{MR,i}={\phi}_{MR,i}+\omega_{i},\quad\hat{R}_{B,i}=R_{{B,i}}+\nu_{i},\vspace{-0.1cm}$
(6)

where $\\{\hat{\theta}_{MR,i},\hat{\phi}_{MR,i},\hat{R}_{B,i}\\}$ denote the
estimation of $\\{{\theta}_{MR,i},{\phi}_{MR,i},R_{{B,i}}\\}$ and
$\\{n_{i},\omega_{i},\nu_{i}\\}$ denote the estimation error. For the sake of
illustration, we collect all the estimation positioning parameters in the
following vectors

$\hat{{\bm{\theta}}}_{\textrm{r}}={\bm{\theta}}_{\textrm{r}}+\mathbf{n},\quad\hat{\bm{\phi}}_{\textrm{r}}=\bm{\phi}_{\textrm{r}}+{\bm{\omega}},\quad\hat{\bf
R}_{B}={\bf R}_{B}+{\bm{\nu}},$ (7)

where

$\displaystyle\hat{\bm{\theta}}_{\textrm{r}}$
$\displaystyle=[\hat{\theta}_{MR,1},\hat{\theta}_{{MR,2}},...,\hat{\theta}_{{MR,M}}]^{T},\bm{\theta}_{\textrm{r}}=[\theta_{{MR,1}},\theta_{{MR,2}},...,\theta_{{MR,M}}]^{T},\hat{\bm{\phi}}_{\textrm{r}}=[\hat{\phi}_{{MR,1}},\hat{\phi}_{{MR,2}},...,\hat{\phi}_{{MR,M}}]^{T},$
$\displaystyle\bm{\phi}_{\textrm{r}}$
$\displaystyle=[\phi_{{MR,1}},\phi_{{MR,2}},...,\phi_{{MR,M}}]^{T},\hat{\bf
R}_{B}=[\hat{R}_{B,1},\hat{R}_{B,2},...,\hat{R}_{B,M}]^{T},{\bf
R}_{B}=[{R}_{B,1},{R}_{B,2},...,{R}_{B,M}]^{T},$ $\displaystyle\mathbf{n}$
$\displaystyle=[n_{1},n_{2},...,n_{M}]^{T},\quad\bm{\omega}=[\omega_{1},\omega_{2},...,\omega_{M}]^{T},\quad{\bm{\nu}}=[\nu_{1},\nu_{2},...,\nu_{M}]^{T}.$
(8)

For tractability, $\mathbf{n}$, $\bm{\omega}$ and ${\bm{\nu}}$ are assumed to
be the additive zero-mean complex Gaussian noise [29]. Furthermore, the
covariance matrices can be written as

$\displaystyle{\bf
Q}_{n}=\left[\begin{array}[]{ccc}\sigma^{2}_{n_{1}}&\cdots&0\\\
\vdots&\ddots&\vdots\\\ 0&\cdots&\sigma^{2}_{n_{M}}\\\
\end{array}\right],\quad{\bf
Q}_{\omega}=\left[\begin{array}[]{ccc}\sigma^{2}_{\omega_{1}}&\cdots&0\\\
\vdots&\ddots&\vdots\\\ 0&\cdots&\sigma^{2}_{\omega_{M}}\\\
\end{array}\right],\quad{\bf
Q}_{\nu}=\left[\begin{array}[]{ccc}\sigma^{2}_{\nu_{1}}&\cdots&0\\\
\vdots&\ddots&\vdots\\\ 0&\cdots&\sigma^{2}_{\nu_{M}}\\\ \end{array}\right],$
(18)

where $\sigma^{2}_{n_{i}}$, $\sigma^{2}_{\omega_{i}}$ and
$\sigma^{2}_{\nu_{i}}$ denote the variance of $n_{i}$, $\omega_{i}$ and
$\nu_{i}$, respectively.

This paper focuses on the position estimation algorithm design, and the phase
shift matrices of the RISs are set to unit matrix for the sake of analysis.

##  III Closed-form Solution

In this section, we propose an algorithm to derive a closed-form solution to
the MU’s position. First, the pseudolinear equations based on the AOA and TDOA
information of the RISs are derived. Then, an algorithm is proposed to
estimate the position of the MU based on the pseudolinear equations.

###  III-A Pseudolinear equations based on the AOA estimation

In this subsection, we derive the pseudolinear equations based on the geometry
relationship of AOAs at the $i$th RIS of the MU-RIS-BS link.

First, from (1), we have

$\cos\theta_{MR,i}(x_{q}-x_{i})-\sin\theta_{MR,i}(y_{q}-y_{i})=0.$ (19)

By replacing the true value $\theta_{MR,i}$ with the estimated
$\hat{\theta}_{MR,i}$, we have

$\eta_{\theta_{MR,i}}=\cos\hat{\theta}_{MR,i}(x_{q}-x_{i})-\sin\hat{\theta}_{MR,i}(y_{q}-y_{i}),$
(20)

where $\eta_{\theta_{MR,i}}$ denotes the residual error of $\theta_{MR,i}$ due
to the estimation error. For the sake of analysis, let
$\breve{\mathbf{g}}_{\theta_{MR,i}}=[-\cos\hat{\theta}_{MR,i}$
$,\sin\hat{\theta}_{MR,i},0]^{T}$. Using the definitions of ${\bf q}$ and
${\bf s}_{i}$, (20) can be rewritten as

$\eta_{\theta_{MR,i}}=\breve{\mathbf{g}}_{\theta_{MR,i}}^{T}{\bf
s}_{i}-\breve{\mathbf{g}}_{\theta_{MR,i}}^{T}\mathbf{q}.$ (21)

Assuming that the angle estimation algorithms have good performance[30], thus
it is reasonable to assume that the estimation error is very small. Hence, we
have $\sin(n_{i})\approx n_{i}$ and $\cos(n_{i})\approx 1$, where $n_{i}$ is
the estimation error of $\theta_{MR,i}$. Consequently, we have the following
approximations

$\displaystyle\sin\hat{\theta}_{MR,i}=\sin(\theta_{MR,i}+n_{i})\approx\sin\theta_{MR,i}+n_{i}\cos\theta_{MR,i},$
$\displaystyle\cos\hat{\theta}_{MR,i}=\cos(\theta_{MR,i}+n_{i})\approx\cos\theta_{MR,i}-n_{i}\sin\theta_{MR,i}.$
(22)

By substituting (III-A) into (20), we have

$\displaystyle\eta_{\theta_{MR,i}}$
$\displaystyle\approx(\cos\theta_{MR,i}-n_{i}\sin\theta_{MR,i})(x_{q}-x_{i})-(\sin\theta_{MR,i}+n_{i}\cos\theta_{MR,i})(y_{q}-y_{i})$
$\displaystyle=-n_{i}R_{{RU,i}}\cos\phi_{MR,i}.$ (23)

In (III-A), (19) and
$(x_{q}-x_{i})\sin\theta_{MR,i}+(y_{q}-y_{i})\cos\theta_{MR,i}=l_{i1}+l_{i2}=R_{{RU,i}}\cos\phi_{MR,i}$
shown in Fig. 2-(a) are used.

(a) The azimuth AOA at the $i$th RIS of the MU-RIS-BS link.

(b) The elevation AOA at the $i$th RIS of the MU-RIS-BS link.

Figure 2: The AOAs at the $i$th RIS of the MU-RIS-BS link.

Next, using (2), we have

$0=\cos{\phi}_{MR,i}(z_{q}-z_{i})-\sin{\phi}_{MR,i}\left(\sin\theta_{MR,i}(x_{q}-x_{i})+\cos\theta_{MR,i}(y_{q}-y_{i})\right).$
(24)

By replacing $\\{\phi_{MR,i},\theta_{MR,i}\\}$ with
$\\{\hat{\phi}_{MR,i},\hat{\theta}_{MR,i}\\}$, we have

$\eta_{\phi_{MR,i}}=\cos\hat{\phi}_{MR,i}(z_{q}-z_{i})-\sin\hat{\phi}_{MR,i}\sin\hat{\theta}_{MR,i}(x_{q}-x_{i})-\sin\hat{\phi}_{MR,i}\cos\hat{\theta}_{MR,i}(y_{q}-y_{i}),$
(25)

where $\eta_{\phi_{MR,i}}$ is the residual error of $\phi_{MR,i}$. For
simplicity, we define
$\breve{\mathbf{g}}_{\phi_{MR,i}}=[\sin\hat{\phi}_{MR,i}\sin\hat{\theta}_{MR,i},$
$\sin\hat{\phi}_{MR,i}\cos\hat{\theta}_{MR,i},-\cos\hat{\phi}_{MR,i}]^{T}$.
Then, by utilizing the definitions of ${\bf q}$ and ${\bf s}_{i}$, (25) can be
expressed as

$\eta_{\phi_{MR,i}}=\breve{\mathbf{g}}_{\phi_{MR,i}}^{T}{\bf
s}_{i}-\breve{\mathbf{g}}_{\phi_{MR,i}}^{T}{\mathbf{q}}.$ (26)

It is assumed that the estimation error is very small [30], we have
$\sin(\omega_{i})\approx\omega_{i}$ and $\cos(\omega_{i})\approx 1$, where
$\omega_{i}$ is the estimation error of $\phi_{MR,i}$. As a result, we have
the following approximations as

$\displaystyle\sin\hat{\phi}_{MR,i}=\sin(\phi_{MR,i}+\omega_{i})\approx\sin\phi_{MR,i}+\omega_{i}\cos\phi_{MR,i},$
$\displaystyle\cos\hat{\phi}_{MR,i}=\cos(\phi_{MR,i}+\omega_{i})\approx\cos\phi_{MR,i}-\omega_{i}\sin\phi_{MR,i}.$
(27)

Then, by substituting (III-A) into (25), and performing some mathematical
manipulations, we have

$\displaystyle\eta_{\phi_{MR,i}}$
$\displaystyle\approx(\cos\phi_{MR,i}-\omega_{i}\sin\phi_{MR,i})(z_{q}-z_{i})$
$\displaystyle\quad-(\sin\phi_{MR,i}+\omega_{i}\cos\phi_{MR,i})((x_{q}-x_{i})\sin\theta_{MR,i}+(y_{q}-y_{i})\cos\theta_{MR,i})$
$\displaystyle=-\omega_{i}R_{{RU,i}}.$ (28)

To derive (III-A), we have used (24) and
$((z_{q}-z_{i})\sin\phi_{MR,i}+((x_{q}-x_{i})\sin{\theta}_{MR,i}+(y_{q}-y_{i})\cos{\theta}_{MR,i})$
$\cos\phi_{MR,i})=d_{i1}+d_{i2}=R_{{RU,i}}$, which can be verified using Fig.
2-(b).

Finally, using (21), (III-A), (26) and (III-A), we have

$\breve{\mathbf{g}}_{\theta_{MR,i}}^{T}\mathbf{s}_{i}-\breve{\mathbf{g}}_{\theta_{MR,i}}^{T}\mathbf{q}\approx-
n_{i}R_{{RU,i}}\cos\phi_{MR,i},\quad\breve{\mathbf{g}}_{\phi_{MR,i}}^{T}\mathbf{s}_{i}-\breve{\mathbf{g}}_{\phi_{MR,i}}^{T}\mathbf{q}\approx-\omega_{i}R_{{RU,i}}.$
(29)

###  III-B Pseudolinear equations based on the TDOA estimation

TDOA estimation is used in this subsection to derive the pseudolinear
equations. First, by taking the square of both sides of (3), we have

$R^{2}_{{BU}}=K_{q}^{2}+K_{p}^{2}-2x_{q}x_{p}-2y_{q}y_{p}-2z_{q}z_{p},$ (30)

where $K_{q}^{2}=x_{q}^{2}+y_{q}^{2}+z_{q}^{2}$ and
$K_{p}^{2}=x_{p}^{2}+y_{p}^{2}+z_{p}^{2}$. Similarly, by taking the square of
both sides of (4), we have

$\displaystyle
R^{2}_{{RU,i}}=K_{q}^{2}+K_{i}^{2}-2x_{q}x_{i}-2y_{q}y_{i}-2z_{q}z_{i},$ (31)

where $K_{i}^{2}=x_{i}^{2}+y_{i}^{2}+z_{i}^{2}$. Then, using (5), we can
obtain $R_{{RU,i}}=R_{{B,i}}+R_{{BU}}$. Thus, we have

$\displaystyle R^{2}_{{RU,i}}=(R_{{B,i}}+R_{{BU}})^{2}.$ (32)

By substituting (31) into (32) and expanding the right hand side of (32), (32)
is rewritten as

$\displaystyle
R^{2}_{{B,i}}+2R_{{B,i}}R_{{BU}}+R^{2}_{{BU}}=K_{q}^{2}+K_{i}^{2}-2x_{q}x_{i}-2y_{q}y_{i}-2z_{q}z_{i}.$
(33)

Moreover, using (30) and (33), we have the following expression

$\displaystyle R^{2}_{{B,i}}+2R_{{B,i}}R_{{BU}}$
$\displaystyle=K_{i}^{2}-K_{p}^{2}-2x_{{i,p}}x_{q}-2y_{{i,p}}y_{q}-2z_{{i,p}}z_{q},$
(34)

where $x_{{i,p}}=(x_{i}-x_{p})$, $y_{{i,p}}=(y_{i}-y_{p})$ and
$z_{{i,p}}=(z_{i}-z_{q})$. In order to derive the pseudolinear equations
related to TDOA, (34) can be transformed into

$\displaystyle-
x_{{i,p}}x_{q}-y_{{i,p}}y_{q}-z_{{i,p}}z_{q}-R_{{B,i}}R_{{BU}}=-\frac{1}{2}K_{i}^{2}+\frac{1}{2}K_{p}^{2}+\frac{1}{2}R^{2}_{{B,i}}.$
(35)

Additionally, for the sake of analysis, we define

$\displaystyle\mathbf{g}_{{t,i}}=-[x_{{i,p}},y_{{i,p}},z_{{i,p}},R_{{B,i}}]^{T},\quad\mathbf{u}=[x_{q},y_{q},z_{q},R_{BU}]^{T},\quad
h_{i}=-\frac{1}{2}K_{i}^{2}+\frac{1}{2}K_{p}^{2}+\frac{1}{2}R^{2}_{{B,i}}.$
(36)

Accordingly, we can derive the pseudolinear equation with respect to the TDOA
as

$\displaystyle h_{i}-\mathbf{g}^{T}_{{t,i}}\mathbf{u}=0.$ (37)

As $R_{B,i}$ is estimated as $\hat{R}_{B,i}$ by applying the TDOA estimation,
by replacing $R_{B,i}$ with $\hat{R}_{B,i}$, $\mathbf{g}_{{t,i}}$ and
${h}_{i}$ in (37) are re-interpreted as $\hat{\mathbf{g}}_{{t,i}}$ and
$\hat{h}_{i}$, given by

$\displaystyle\hat{\mathbf{g}}_{{t,i}}$
$\displaystyle=-[x_{{i,p}},y_{{i,p}},z_{{i,p}},\hat{R}_{{B,i}}]^{T},\quad\hat{h}_{i}=-\frac{1}{2}K_{i}^{2}+\frac{1}{2}K_{p}^{2}+\frac{1}{2}\hat{R}^{2}_{{B,i}}.$
(38)

Thus, according to [29], (37) can be derived as

$\displaystyle\eta_{t}=\hat{h}_{i}-\hat{\mathbf{g}}_{{t,i}}^{T}\mathbf{u}=R_{{RU,i}}\nu_{i}+\frac{1}{2}\nu_{i}^{2},$
(39)

where $\eta_{t}$ and $\nu_{i}$ denote the residual error and the TDOA
estimation error of ${R}_{B,i}$ in (6), respectively. As $\nu_{i}\ll
R_{{RU,i}}$ is always satisfied in practice, the second order term on the
right hand side of (39) can be neglected, and (39) can be approximated as

$\displaystyle\hat{h}_{i}-\hat{\mathbf{g}}_{{t,i}}^{T}\mathbf{u}\approx
R_{{RU,i}}\nu_{i}.$ (40)

###  III-C Closed-form solution of the position of MU

In this subsection, we derive the compact form of the pseudolinear equations
with respect to the AOA and TDOA estimation. We also derive the closed-form
solution of the MU’s position using the TSWLS algorithm, the details are
summarized in Algorithm 1 and the flowchart of proposed algorithm is shown in
Fig. 3.

Figure 3: Flowchart of TSWLS algorithm. Algorithm 1 Two-Stage Weighted Least
Squares (TSWLS) Algorithm

1: Construct the pseudolinear equations based on the AOA and TDOA estimation;

2: Calculate the preliminary estimation $\breve{\mathbf{u}}$ at _Stage 1_
using Algorithm 2;

3: Construct the pseudolinear equations based on the preliminary estimation at
_Stage 1_ ;

4: Calculate the estimation of ${\bm{\xi}}=\left({\bf q}-{\bf
p}\right)\odot\left({\bf q}-{\bf p}\right)$, denoted by $\breve{\bm{\xi}}$;

5: Derive the final estimation of ${\bf q}$ using $\breve{\bm{\xi}}$ and
$\breve{\mathbf{u}}$.

First, we define

$\displaystyle\hat{\mathbf{g}}_{\theta_{MR,i}}=[\breve{\mathbf{g}}^{T}_{\theta_{MR,i}},0]^{T},\quad\hat{\mathbf{g}}_{\phi_{MR,i}}=[\breve{\mathbf{g}}^{T}_{\phi_{MR,i}},0]^{T}.$
(41)

A better estimation performance can be obtained using (29) together with (40),
both of which contain the unknown location ${\bf q}$. Thus, we can derive the
following compact form of equations as

$\hat{\mathbf{h}}_{\textrm{r}}-\hat{\mathbf{G}}_{\textrm{r}}\mathbf{u}=\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}},$
(42)

where
$\mathbf{z}_{\textrm{r}}=[\mathbf{n}^{T},\bm{\omega}^{T},\bm{\nu}^{T}]^{T}$
and its covariance matrix is written as

$\mathbf{Q}_{\textrm{r}}=\left[\begin{array}[]{ccc}\mathbf{Q}_{n}&\mathbf{O}&\mathbf{O}\\\
\mathbf{O}&\mathbf{Q}_{\omega}&\mathbf{O}\\\
\mathbf{O}&\mathbf{O}&\mathbf{Q}_{\nu}\\\ \end{array}\right].$ (43)

On the left hand side of (42), we have

$\displaystyle\hat{\mathbf{h}}_{\textrm{r}}$
$\displaystyle=[\mathbf{1}^{T}(\breve{\mathbf{G}}_{\theta_{\textrm{r}}}\odot\mathbf{S})^{T},\mathbf{1}^{T}(\breve{\mathbf{G}}_{\phi_{\textrm{r}}}\odot\mathbf{S})^{T},\hat{\mathbf{h}}^{T}_{{t}}]^{T},$
$\displaystyle\hat{\mathbf{G}}_{\textrm{r}}$
$\displaystyle=[\hat{\mathbf{G}}^{T}_{\theta_{\textrm{r}}},\hat{\mathbf{G}}^{T}_{\phi_{\textrm{r}}},\hat{\mathbf{G}}^{T}_{{t}}]^{T},$
(44)

where

$\displaystyle\breve{\mathbf{G}}_{\theta_{\textrm{r}}}$
$\displaystyle=[\breve{\mathbf{g}}_{\theta_{{MR,1}}},\breve{\mathbf{g}}_{\theta_{{MR,2}}},...,\breve{\mathbf{g}}_{\theta_{{MR,M}}}]^{T},\breve{\mathbf{G}}_{\phi_{\textrm{r}}}=[\breve{\mathbf{g}}_{\phi_{{MR,1}}},\breve{\mathbf{g}}_{\phi_{{MR,2}}},...,\breve{\mathbf{g}}_{\phi_{{MR,M}}}]^{T},$
$\displaystyle\hat{\mathbf{G}}_{\theta_{\textrm{r}}}$
$\displaystyle=[\hat{\mathbf{g}}_{\theta_{{MR,1}}},\hat{\mathbf{g}}_{\theta_{{MR,2}}},...,\hat{\mathbf{g}}_{\theta_{{MR,M}}}]^{T},\hat{\mathbf{G}}_{\phi_{\textrm{r}}}=[\hat{\mathbf{g}}_{\phi_{{MR,1}}},\hat{\mathbf{g}}_{\phi_{{MR,2}}},...,\hat{\mathbf{g}}_{\phi_{{MR,M}}}]^{T},$
$\displaystyle\hat{\mathbf{h}}_{{t}}$
$\displaystyle=[\hat{h}_{{1}},\hat{h}_{{2}},...,\hat{h}_{{M}}]^{T},\quad\mathbf{S}=[\mathbf{s}_{{1}},\mathbf{s}_{{2}},...,\mathbf{s}_{{M}}]^{T},\quad\hat{\mathbf{G}}_{{t}}=[\hat{\mathbf{g}}_{{t,1}},\hat{\mathbf{g}}_{{t,2}},...,\hat{\mathbf{g}}_{{t,M}}]^{T}.$
(45)

On the right hand side of (42), we have

$\displaystyle\mathbf{B}_{\textrm{r}}$
$\displaystyle=[\mathbf{B}_{\textrm{r,n}},\mathbf{B}_{\textrm{r,}\omega},\mathbf{B}_{\textrm{r,}\nu}]^{T},$
(46)

where

$\displaystyle\mathbf{B}_{\textrm{r,n}}=[\mathbf{B}_{\theta_{{MR,n}}},\mathbf{O},\mathbf{O}]^{T},\quad\mathbf{B}_{{\textrm{r},}\omega}=[\mathbf{O},\mathbf{B}_{\phi_{{MR,}\omega}},\mathbf{O}]^{T},\quad\mathbf{B}_{{\textrm{r},}\nu}=[\mathbf{O},\mathbf{O},\mathbf{B}_{{t}}]^{T},$
$\displaystyle\mathbf{B}_{\theta_{{MR,n}}}=-\textrm{diag}[R_{{RU,1}}\cos\phi_{{MR,1}},...,R_{{RU,1}}\cos\phi_{{MR,M}}],$
$\displaystyle\mathbf{B}_{\phi_{{MR,}\omega}}=-\textrm{diag}[R_{RU,1},...,R_{RU,M}],\quad\mathbf{B}_{t}=\textrm{diag}[R_{RU,1},...,R_{RU,M}].$
(47)

Using the compact pseudolinear equations (42), we propose the TSWLS algorithm
to estimate the position of the MU as follows.

####  III-C1 Stage 1

In the first stage, we aim to derive the estimation of $\mathbf{u}$ using the
WLS method [23]. Based on (42), the WLS cost function can be formulated as

$f({\mathbf{u}})=(\hat{{\bf
h}}_{\textrm{r}}-\hat{\mathbf{G}}_{\textrm{r}}{\mathbf{u}})^{T}\mathbf{W}_{\textrm{r}}(\hat{\mathbf{h}}_{\textrm{r}}-\hat{\mathbf{G}}_{\textrm{r}}{\mathbf{u}})=\hat{\mathbf{h}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{h}}_{\textrm{r}}-2\hat{\mathbf{h}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}}{\mathbf{u}}+{\mathbf{u}}^{T}\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}}{\mathbf{u}},$
(48)

where $\mathbf{W}_{\textrm{r}}$ denotes the weight matrix of _Stage 1_.

To derive the estimation of ${\bf u}$, the cost function $f({\mathbf{u}})$ in
(48) should be minimized. By setting the first-order derivative of
$f({\mathbf{u}})$ equal to zero, we have $\frac{\partial
f({\mathbf{u}})}{\partial{\mathbf{u}}}=-2(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{h}}_{\textrm{r}})+2(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}}{\mathbf{u}})=0$.
Hence, $\mathbf{u}$ is estimated as

$\breve{\mathbf{u}}=[\breve{x}_{q},\breve{y}_{q},\breve{z}_{q},\breve{R}_{BU}]^{T}=(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1}\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{h}}_{\textrm{r}}.$
(49)

As the estimation error is correlated, the weight matrix
$\mathbf{W}_{\textrm{r}}$ should be equal to the inverse of the covariance
matrix of the estimation error as

$\mathbf{W}_{\textrm{r}}={\bm{\Omega}}_{\textrm{r}}^{-1},{\bm{\Omega}}_{\textrm{r}}=\mathbb{E}[\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}\mathbf{z}_{\textrm{r}}^{T}\mathbf{B}_{\textrm{r}}^{T}]=\mathbf{B}_{\textrm{r}}\mathbf{Q}_{\textrm{r}}\mathbf{B}_{\textrm{r}}^{T},$
(50)

where ${\bm{\Omega}}_{\textrm{r}}$ denotes the covariance matrix of the
estimation error. However, as shown in (46) and (III-C), matrix
$\mathbf{B}_{\textrm{r}}$ contains the true distances
$\\{R_{{RU,1}},R_{{RU,2}},\cdots,R_{{RU,M}}\\}$, which remain unknown.
Therefore, the weight matrix is initialized by taking the estimated distances
as the approximation of true values. Then, the weight matrix can be
iteratively updated by the new estimations. The detailed initialization
algorithm is provided in Algorithm 2.

Algorithm 2 Initial Location Estimation

1: Set $\hat{\mathbf{W}}_{\textrm{r}}=\mathbf{I}$;

2: Derive the initial estimation $\breve{\mathbf{u}}$ in (49) by replacing
$\mathbf{W}_{\textrm{r}}$ with $\hat{\mathbf{W}}_{\textrm{r}}$;

3: Calculate $\\{\hat{R}_{{RU,1}},\hat{R}_{{RU,2}},...,\hat{R}_{{RU,M}}\\}$ to
replace $\\{R_{{RU,1}},R_{{RU,2}},...,R_{{RU,M}}\\}$ using the initial
estimation $\breve{\mathbf{u}}$;

4: Derive $\hat{\mathbf{B}}_{\textrm{r}}$ to replace $\mathbf{B}_{\textrm{r}}$
in (50) using the definition of $\mathbf{B}_{\textrm{r}}$ in (46) and
$\\{\hat{R}_{{RU,1}},\hat{R}_{{RU,2}},...,\hat{R}_{{RU,M}}\\}$;

5: Derive $\hat{\mathbf{W}}_{\textrm{r}}$ to replace $\mathbf{W}_{\textrm{r}}$
in (49) with $\hat{\mathbf{B}}_{\textrm{r}}$ based on (50) ;

6: Repeat.

####  III-C2 Stage 2

As we can see from (3), $R_{BU}$ and ${\bf q}=[x_{q},y_{q},z_{q}]^{T}$ are
related, which is ignored in the first stage. Therefore, we aim to utilize
this relationship to construct a new set of pseudolinear equations in this
stage.

Considering the estimation error, the relationship between the MU’s
preliminary estimated position and the BS’s position is

$\left[{\begin{array}[]{*{20}{c}}(\breve{x}_{q}-x_{p})^{2},(\breve{y}_{q}-y_{p})^{2},(\breve{z}_{q}-z_{p})^{2},\breve{R}_{BU}^{2}\end{array}}\right]^{T}\buildrel\Delta\over{=}\hat{{\bf
h}}_{1}.$ (51)

Without the estimation error, the relationship between the MU’s position and
the BS’s position is represented by

$\left[{\begin{array}[]{*{20}{c}}({x}_{q}-x_{p})^{2},({y}_{q}-y_{p})^{2},({z}_{q}-z_{p})^{2},{R}_{BU}^{2}\end{array}}\right]^{T}\buildrel\Delta\over{=}{\bf
G}_{1}{\bm{\xi}},$ (52)

where

${\bf G}_{1}=\left[{\begin{array}[]{*{20}{c}}{\bf I}_{3\times 3}\\\ {\bf
1}_{1\times
3}\end{array}}\right],\quad{\bm{\xi}}=\left[{\begin{array}[]{*{20}{c}}(x_{q}-x_{p})^{2},(y_{q}-y_{p})^{2},(z_{q}-z_{p})^{2}\end{array}}\right]^{T}=\left({\bf
q}-{\bf p}\right)\odot\left({\bf q}-{\bf p}\right).$ (53)

Using (51) and (52), we can derive the compact form of the pseudolinear
equations in _Stage 2_ as

$\hat{{\bf h}}_{1}-{\bf G}_{1}{\bm{\xi}}={\bf z}_{1},$ (54)

where ${\bf
z}_{1}=[z_{11},z_{12},z_{13},z_{14}]^{T}=[(\breve{x}_{q}-x_{p})^{2},(\breve{y}_{q}-y_{p})^{2},(\breve{z}_{q}-z_{p})^{2},\breve{R}_{BU}^{2}]^{T}-[({x}_{q}-x_{p})^{2},({y}_{q}-y_{p})^{2},({z}_{q}-z_{p})^{2},{R}_{BU}^{2}]^{T}$
denotes the residual vector due to the estimation error. Denote the
positioning estimation error of the first stage as $\tilde{\bf
u}=[e_{1},e_{2},e_{3},e_{4}]^{T}$, i.e ., $\breve{{\bf u}}={{\bf
u}}+\tilde{\bf u}$. Then, the elements of ${\bf z}_{1}$ can be represented as
the functions of the elements of error $\tilde{\bf u}$, which can be expressed
as

$\displaystyle z_{11}$
$\displaystyle=(\breve{x}_{q}-x_{p})^{2}-(x_{q}-x_{p})^{2}=2(x_{q}-x_{p})e_{1}+e_{1}^{2}\approx
2(x_{q}-x_{p})e_{1},$ $\displaystyle z_{12}$
$\displaystyle=(\breve{y}_{q}-y_{p})^{2}-(y_{q}-y_{p})^{2}=2(y_{q}-y_{p})e_{2}+e_{2}^{2}\approx
2(y_{q}-y_{p})e_{2},$ $\displaystyle z_{13}$
$\displaystyle=(\breve{z}_{q}-z_{p})^{2}-(z_{q}-z_{p})^{2}=2(z_{q}-z_{p})e_{3}+e_{3}^{2}\approx
2(z_{q}-z_{p})e_{3},$ $\displaystyle z_{14}$
$\displaystyle=\breve{R}_{BU}^{2}-R_{BU}^{2}=2R_{BU}e_{4}+e_{4}^{2}\approx
2R_{BU}e_{4}.$ (55)

Based on (54), we can again apply the WLS method to estimate ${\bm{\xi}}$,
thus the WLS cost function can be written as

$f({\bm{\xi}})=(\hat{\bf
h}_{1}-\mathbf{G}_{1}{\bm{\xi}})^{T}{\mathbf{W}}_{1}(\hat{\mathbf{h}}_{1}-\mathbf{G}_{1}{\bm{\xi}})=\hat{\mathbf{h}}_{1}^{T}{\mathbf{W}}_{1}\hat{\mathbf{h}}_{1}-2\hat{\mathbf{h}}_{1}^{T}{\mathbf{W}}_{1}\mathbf{G}_{1}{\bm{\xi}}+{\bm{\xi}}^{T}\mathbf{G}_{1}^{T}{\mathbf{W}}_{1}\mathbf{G}_{1}{\bm{\xi}},$
(56)

where ${\mathbf{W}}_{1}$ denotes the weight matrix at _Stage 2_. To derive the
estimation of ${\bm{\xi}}$, $f({\bm{\xi}})$ in (56) should be minimized,
leading to $\frac{\partial
f({\bm{\xi}})}{\partial{\bm{\xi}}}=-2(\mathbf{G}_{1}^{T}{\mathbf{W}}_{1}{\hat{\mathbf{h}}}_{1})+2(\mathbf{G}_{1}^{T}{\mathbf{W}}_{1}\mathbf{G}_{1}{\bm{\xi}})=0$.
Therefore, the estimation of ${\bm{\xi}}$ at _Stage 2_ is given by

$\breve{{\bm{\xi}}}=({\bf G}_{1}^{T}{\mathbf{W}}_{1}{\bf G}_{1})^{-1}{\bf
G}^{T}_{1}{\bf W}_{1}\hat{\bf h}_{1}.$ (57)

By defining ${\bm{\Omega}}_{1}=\mathbb{E}[{\bf z}_{1}{\bf z}_{1}^{T}]$ as the
covariance matrix of ${\bf z}_{1}$, we have
$\mathbf{W}_{1}={\bm{\Omega}}^{-1}_{1}$. Moreover, ${\bm{\Omega}}_{1}$ can be
further derived as

${\bm{\Omega}}_{1}={{\bf B}}_{1}(\hat{{\bf G}}_{\textrm{r}}^{T}{\bf
W}_{\textrm{r}}\hat{{\bf G}}_{\textrm{r}})^{-1}{\mathbf{B}}_{1},\quad{{\bf
B}}_{1}=2\textrm{diag}\left\\{{{\bf u}}-\left[{\begin{array}[]{*{20}{c}}{\bf
p}\\\ 0\end{array}}\right]\right\\},$ (58)

where the proof of (58) is given in Appendix A. Although ${\bf B}_{1}$
contains the true values $\\{x_{q},y_{q},z_{q},R_{BU}\\}$, it can be
approximated as

$\hat{\bf B}_{1}=2\textrm{diag}\left\\{\breve{{\bf
u}}-\left[{\begin{array}[]{*{20}{c}}{\bf p}\\\ 0\end{array}}\right]\right\\}.$
(59)

Then, ${\bm{\Omega}}_{1}$ can be approximated as
$\hat{\bm{\Omega}}_{1}=\hat{\bf B}_{1}(\hat{{\bf G}}_{\textrm{r}}^{T}{{\bf
W}}_{\textrm{r}}\hat{{\bf G}}_{\textrm{r}})^{-1}\hat{\bf B}_{1}$. Hence, ${\bf
W}_{1}$ can be estimated as $\hat{\bf
W}_{1}=\hat{\bm{\Omega}}_{1}^{-1}=\hat{\bf B}_{1}^{-1}(\hat{{\bf
G}}_{\textrm{r}}^{T}{{\bf W}}_{\textrm{r}}\hat{{\bf G}}_{\textrm{r}})\hat{\bf
B}_{1}^{-1}$, and $\breve{{\bm{\xi}}}$ is approximated as

$\breve{{\bm{\xi}}}\approx({\bf G}_{1}^{T}\hat{\mathbf{W}}_{1}{\bf
G}_{1})^{-1}{\bf G}^{T}_{1}\hat{\bf W}_{1}\hat{\bf h}_{1}.$ (60)

####  III-C3 Final Solution

Based on the estimations in _Stages 1_ and _2_ , the final closed-form
expression of the position of MU can be expressed as:

$\hat{{\bf q}}={\bm{\Pi}}\sqrt{\breve{{\bm{\xi}}}}+{\bf
p},\quad{\bm{\Pi}}=\textrm{diag}\\{\textrm{sgn}(\breve{\mathbf{u}}(1:3)-{\bf
p})\\},$ (61)

where $\textrm{sgn}(\textrm{x})$ denotes the signum function111The signum
function $\textrm{sgn}(\textrm{x})$ is a piecewise function which is define
as:
$\textrm{sgn}(\textrm{x})=\left\\{\begin{array}[]{rcl}-1,&&{\textrm{x}<0}\\\
0,&&{\textrm{x}=0}\\\ 1,&&{\textrm{x}>0.}\end{array}\right.$.

##  IV Bias Analysis

To analyze the bias here, we use the second-order statistics of the estimation
error. In addition, it is assumed that the error level is not large, which
enables us to ignore the error terms of the third order or higher.
Furthermore, we need to derive the bias in both stages so that the bias for
the final estimation in the previous section can be determined.

As the estimation error is inevitable, $\hat{\mathbf{G}}_{\textrm{r}}$ defined
in (III-C) consists of two parts: the true matrix ${\bf G}_{\textrm{r}}$ and
the error matrix $\tilde{\bf{G}}_{\textrm{r}}$, i.e.,
$\hat{\mathbf{G}}_{\textrm{r}}=\mathbf{G}_{\textrm{r}}+\tilde{\mathbf{G}}_{\textrm{r}}$.
According to (III-C) and (III-C), to derive the expressions of
$\mathbf{G}_{\textrm{r}}$ and $\tilde{\mathbf{G}}_{\textrm{r}}$, we need to
decompose $\hat{\mathbf{g}}_{\theta_{MR,i}}$, $\hat{\mathbf{g}}_{\phi_{MR,i}}$
and $\hat{\mathbf{g}}_{t,i}$ into the true values and the estimation error at
first. The flowchart of bias analysis is shown in Fig.4.

First, let us decompose $\hat{\mathbf{g}}_{\theta_{MR,i}}$ into the true value
and the estimation error. Using the approximations in (III-A) and the
definition of $\hat{\mathbf{g}}_{\theta_{MR,i}}$ in (41),
$\hat{\mathbf{g}}_{\theta_{MR,i}}$ is approximated as
$\hat{\mathbf{g}}_{\theta_{MR,i}}\approx[-\cos\theta_{MR,i},\sin\theta_{MR,i},0,0]^{T}+n_{i}[\sin\theta_{MR,i},\cos\theta_{MR,i},0,0]^{T}$.
By defining $\bar{\mathbf{g}}_{\theta_{MR,i}}=[-\cos\theta_{MR,i},$
$\sin\theta_{MR,i},0,0]^{T}$ and ${\bf
g}_{1n_{i}}=[\sin\theta_{MR,i},\cos\theta_{MR,i},0,0]^{T}$, we have

$\hat{\mathbf{g}}_{\theta_{MR,i}}\approx\bar{\mathbf{g}}_{\theta_{MR,i}}+n_{i}{\bf
g}_{1n_{i}},$ (62)

where $n_{i}{\bf g}_{1n_{i}}$ denotes the error term.

Similarly, by utilizing the approximations in (III-A), (III-A) and the
definition of $\hat{\mathbf{g}}_{\phi_{MR,i}}$ in (41),
$\hat{\mathbf{g}}_{\phi_{MR,i}}$ is approximated as

$\displaystyle\hat{\mathbf{g}}_{\phi_{MR,i}}\approx$
$\displaystyle[(\sin\phi_{MR,i}+\omega_{i}\cos\phi_{MR,i})(\sin\theta_{MR,i}+n_{i}\cos\theta_{MR,i}),$
$\displaystyle(\sin\phi_{MR,i}+\omega_{i}\cos\phi_{MR,i})(\cos\theta_{MR,i}-n_{i}\sin\theta_{MR,i}),-(\cos\phi_{MR,i}-\omega_{i}\sin\phi_{MR,i}),0]^{T}.$

By ignoring the terms $\omega_{i}n_{i}\cos\phi_{MR,i}\cos\theta_{MR,i}$ and
$-\omega_{i}n_{i}\cos\phi_{MR,i}\sin\theta_{MR,i}$ [31] and defining
$\bar{\mathbf{g}}_{\phi_{MR,i}}=\left[\sin\phi_{MR,i}\sin\theta_{MR,i},\sin\phi_{MR,i}\cos\theta_{MR,i},-\cos\phi_{MR,i},0\right]^{T}$,
$\mathbf{g}_{\omega_{i}}=[\cos\phi_{MR,i}$ $\sin\theta_{MR,i},$
$\cos\phi_{MR,i}\cos\theta_{MR,i},\sin\phi_{MR,i},0]^{T}$ and
$\mathbf{g}_{2n_{i}}=[\sin\phi_{MR,i}\cos\theta_{MR,i},$
$-\sin\phi_{MR,i}\sin\theta_{MR,i},$ $0,0]^{T}$, we have

$\hat{\mathbf{g}}_{\phi_{MR,i}}\approx\bar{\mathbf{g}}_{\phi_{MR,i}}+\omega_{i}\mathbf{g}_{\omega_{i}}+n_{i}\mathbf{g}_{2n_{i}},$
(63)

where $\omega_{i}\mathbf{g}_{\omega_{i}}$ and $n_{i}\mathbf{g}_{2n_{i}}$
denote the error terms, respectively.

Similarly, we can decompose $\hat{\mathbf{g}}_{t,i}$ into the true value and
the estimation error. Based on the definition of $\mathbf{g}_{t,i}$ in (36)
and $\hat{\mathbf{g}}_{t,i}$ in (38), we can rewrite $\hat{\mathbf{g}}_{t,i}$
as
$\hat{\mathbf{g}}_{t,i}=-[x_{{i,p}},y_{{i,p}},z_{{i,p}},R_{{B,i}}]^{T}-[0,0,0,\nu_{i}]^{T}=\mathbf{g}_{t,i}+\nu_{i}[0,0,0,-1]^{T}$.
Letting $\mathbf{g}_{\nu}=[0,0,0,-1]^{T}$, we have

$\displaystyle\hat{\mathbf{g}}_{t,i}=\mathbf{g}_{t,i}+\nu_{i}\mathbf{g}_{\nu},$
(64)

where $\nu_{i}\mathbf{g}_{\nu}$ denotes the error term.

Therefore, using (62), (63) and (64), $\hat{\mathbf{G}}_{\textrm{r}}$ can be
written as

$\hat{\mathbf{G}}_{\textrm{r}}=\mathbf{G}_{\textrm{r}}+\tilde{\mathbf{G}}_{\textrm{r}}.$
(65)

Moreover, using the definition of ${\bf n}$, $\bm{\omega}$ and $\bm{\nu}$ in
(II-C), $\tilde{\mathbf{G}}_{\textrm{r}}$ can be written as

$\displaystyle\tilde{\mathbf{G}}_{\textrm{r}}=$
$\displaystyle\bar{\mathbf{G}}_{n}\odot[(\mathbf{n1}_{1\times
4})^{T},(\mathbf{n1}_{1\times 4})^{T},(\mathbf{n1}_{1\times
4})^{T}]^{T}+\bar{\mathbf{G}}_{\omega}\odot[(\bm{\omega 1}_{1\times
4})^{T},(\bm{\omega 1}_{1\times 4})^{T},(\bm{\omega 1}_{1\times 4})^{T}]^{T}$
$\displaystyle+\bar{\mathbf{G}}_{\nu}\odot[(\bm{\nu 1}_{1\times
4})^{T},(\bm{\nu 1}_{1\times 4})^{T},(\bm{\nu 1}_{1\times 4})^{T}]^{T},$ (66)

where

$\displaystyle\bar{\mathbf{G}}_{n}$
$\displaystyle=[\mathbf{G}_{1n}^{T},\mathbf{G}_{2n}^{T},\mathbf{O}_{4\times
M}]^{T},\bar{\mathbf{G}}_{\omega}=[\mathbf{O}_{4\times
M},\mathbf{G}_{\omega}^{T},\mathbf{O}_{4\times
M}]^{T},\bar{\mathbf{G}}_{\nu}=[\mathbf{O}_{4\times M},\mathbf{O}_{4\times
M},\mathbf{G}_{\nu}^{T}]^{T},$ $\displaystyle\mathbf{G}_{1n}$
$\displaystyle=[\mathbf{g}_{1n_{1}},\mathbf{g}_{1n_{2}},...,\mathbf{g}_{1n_{M}}]^{T},\mathbf{G}_{2n}=[\mathbf{g}_{2n_{1}},\mathbf{g}_{2n_{2}},...,\mathbf{g}_{2n_{M}}]^{T},$
$\displaystyle\mathbf{G}_{\omega}$
$\displaystyle=[\mathbf{g}_{\omega_{1}},\mathbf{g}_{\omega_{2}},...,\mathbf{g}_{\omega_{M}}]^{T},\mathbf{G}_{\nu}=[\mathbf{g}_{\nu},\mathbf{g}_{\nu},...,\mathbf{g}_{\nu}]^{T}.$
(67)

In addition, $\mathbf{G}_{\textrm{r}}$ is given by

$\mathbf{G}_{\textrm{r}}=[\mathbf{G}_{\theta_{\textrm{r}}}^{T},\mathbf{G}_{\phi_{\textrm{r}}}^{T},\mathbf{G}_{{t}}^{T}]^{T},$
(68)

where

$\displaystyle\mathbf{G}_{\theta_{\textrm{r}}}$
$\displaystyle=[\bar{\mathbf{g}}_{\theta_{MR,1}},\bar{\mathbf{g}}_{\theta_{MR,2}},...,\bar{\mathbf{g}}_{\theta_{MR,M}}]^{T},\mathbf{G}_{\phi_{\textrm{r}}}=[\bar{\mathbf{g}}_{\phi_{{MR,1}}},\bar{\mathbf{g}}_{\phi_{{MR,2}}},...,\bar{\mathbf{g}}_{\phi_{{MR,M}}}]^{T},$
$\displaystyle\mathbf{G}_{t}$
$\displaystyle=[\mathbf{g}_{{t,1}},\mathbf{g}_{{t,2}},...,\mathbf{g}_{{t,M}}]^{T}.$
(69)

Figure 4: Flowchart of bias analysis.

###  IV-A Bias in Stage 1

In this subsection, we aim to derive the bias of $\breve{\mathbf{u}}$ obtained
in _Stage 1_. As we can see from the flowchart in Fig. 4, to derive the bias
of $\breve{\bf u}$, the expression of estimation error of _Stage 1_ should be
derived, which is given by $\tilde{\bf u}=\breve{\bf u}-{\bf u}$. Then, the
bias of $\breve{\bf u}$ is given by taking the expectation of $\tilde{\bf u}$,
which is written as $\mathbb{E}(\tilde{\bf u})$. Using (49) and the definition
of $\mathbf{u}$, $\tilde{\bf u}$ can be further derived as

$\tilde{\mathbf{u}}=\breve{\bf u}-{\bf
u}=(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1}\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}(\hat{\mathbf{h}}_{\textrm{r}}-\hat{\mathbf{G}}_{\textrm{r}}\mathbf{u}).$
(70)

Since we have
$\hat{\mathbf{h}}_{\textrm{r}}-\hat{\mathbf{G}}_{\textrm{r}}\mathbf{u}=\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}$
in (42), $\tilde{\bf u}$ can be derived as $\tilde{\bf
u}=(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1}\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}$.
As we have mentioned above, it is necessary to consider the second order term
when analyzing the bias [28]. However, we have used (40) to derive (42) rather
than (39), which ignores the second order term $\frac{1}{2}\nu_{i}^{2}$.
Therefore, to obtain the expression of $\tilde{\bf u}$ and derive the bias of
_Stage 1_ , $\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}$ should be
extended to
$\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}+[\mathbf{0}^{T},\frac{1}{\sqrt{2}}{\bm{\nu}}^{T}]^{T}\odot[\mathbf{0}^{T},\frac{1}{\sqrt{2}}\mathbf{\bm{\nu}}^{T}]^{T}$,
where $\mathbf{0}$ is a $2M\times 1$ column vector. Then, $\tilde{\bf u}$ can
be further derived as
$\tilde{\mathbf{u}}=(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1}\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\bigg{(}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}+[\mathbf{0}^{T},\frac{1}{\sqrt{2}}{\bm{\nu}}^{T}]^{T}\odot[\mathbf{0}^{T},\frac{1}{\sqrt{2}}\mathbf{\bm{\nu}}^{T}]^{T}\bigg{)}$.
Letting
$\mathbf{\bm{\eta}}=[\mathbf{0}^{T},\frac{1}{\sqrt{2}}\bm{\nu}^{T}]^{T}$ and
$\hat{\mathbf{P}}_{\textrm{r}}=\hat{\mathbf{G}}_{\textrm{r}}^{T}{\mathbf{W}}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}}$,
$\tilde{\mathbf{u}}$ can be rewritten as

$\tilde{\mathbf{u}}=\hat{\mathbf{P}}^{-1}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}(\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}+\mathbf{\bm{\eta}}\odot\mathbf{\bm{\eta}}).$
(71)

According to the definition of $\hat{\mathbf{G}}_{\textrm{r}}$ in (65),
$\hat{\mathbf{P}}_{\textrm{r}}$ can be further derived as

$\hat{\mathbf{P}}_{\textrm{r}}=(\mathbf{G}_{\textrm{r}}+\tilde{\mathbf{G}}_{\textrm{r}})^{T}\mathbf{W}_{\textrm{r}}(\mathbf{G}_{\textrm{r}}+\tilde{\mathbf{G}}_{\textrm{r}})=\mathbf{G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}+\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}+\mathbf{G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}+\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}.$
(72)

Here,
$\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}$
can be neglected222If we consider
$\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}$
when deriving the expression of $\tilde{\bf u}$, it will be multiplied by
$\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}$, and the terms including
$\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}$
in the final expression of $\tilde{\bf u}$ is higher than the second order,
thus we neglect this term here., thus $\hat{\mathbf{P}}_{\textrm{r}}$ can be
approximated as

$\hat{\mathbf{P}}_{\textrm{r}}\approx\mathbf{G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}+\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}+\mathbf{G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}=\mathbf{P}_{\textrm{r}}+\tilde{\mathbf{P}}_{\textrm{r}},$
(73)

where
$\mathbf{P}_{\textrm{r}}=\mathbf{G}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}$
and
$\tilde{\mathbf{P}}_{\textrm{r}}=\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}+\mathbf{G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}$,
respectively. However, if we use the definition of
$\hat{\mathbf{P}}_{\textrm{r}}$ in (73) to derive (71), the inverse matrix
$\hat{\mathbf{P}}_{\textrm{r}}^{-1}=(\mathbf{P}_{\textrm{r}}+\tilde{\mathbf{P}}_{\textrm{r}})^{-1}$
is complex and challenging to derive. Fortunately, we can derive the
approximation of $\hat{\mathbf{P}}_{\textrm{r}}^{-1}$ by utilizing the Newmann
expansion [28] when the error level is small, which is written as

$\hat{\mathbf{P}}_{\textrm{r}}^{-1}\approx\left(\mathbf{I}-\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{P}}_{\textrm{r}}\right)\mathbf{P}_{\textrm{r}}^{-1}.$
(74)

By substituting (74) and (65) into (71), $\tilde{\mathbf{u}}$ is derived as

$\tilde{\mathbf{u}}\approx\left(\mathbf{I}-\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{P}}_{\textrm{r}}\right)\mathbf{P}_{\textrm{r}}^{-1}(\mathbf{G}_{\textrm{r}}+\tilde{\mathbf{G}}_{\textrm{r}})^{T}\mathbf{W}_{\textrm{r}}(\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}+\bm{\eta}\odot\bm{\eta}).$
(75)

For the sake of illustration, using $\mathbf{P}_{\textrm{r}}$ in (73), let us
define
$\mathbf{H}_{\textrm{r}}=(\mathbf{G}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}})^{-1}\mathbf{G}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}=\mathbf{P}^{-1}_{\textrm{r}}\mathbf{G}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}$.
Then, $\tilde{\mathbf{u}}$ can be represented as

$\displaystyle\tilde{\mathbf{u}}$
$\displaystyle\approx\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}+\mathbf{H}_{\textrm{r}}(\bm{\eta}\odot\bm{\eta})-\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{P}}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}-\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{P}}_{\textrm{r}}\mathbf{H}_{\textrm{r}}(\bm{\eta}\odot\bm{\eta})+\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}$
$\displaystyle\quad+\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}(\bm{\eta}\odot\bm{\eta})-\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{P}}_{\textrm{r}}\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}-\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{P}}_{\textrm{r}}\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}(\bm{\eta}\odot\bm{\eta}).$
(76)

For simplicity, we can ignore the error terms in (IV-A) which are higher than
the second order. As a result, using $\tilde{\mathbf{P}}_{\textrm{r}}$ in
(73), $\tilde{\mathbf{u}}$ is approximated as

$\tilde{\mathbf{u}}\approx\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}+\mathbf{H}_{\textrm{r}}(\bm{\eta}\odot\bm{\eta})-\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}-\mathbf{P}_{\textrm{r}}^{-1}\mathbf{G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}+\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}.$
(77)

As a result, the bias in _Stage 1_ is written as

$\displaystyle\mathbb{E}[\tilde{\mathbf{u}}]$
$\displaystyle=\mathbb{E}[\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}+\mathbf{H}_{\textrm{r}}(\bm{\eta}\odot\bm{\eta})-\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}-\mathbf{P}_{\textrm{r}}^{-1}\mathbf{G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}+\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}]$
$\displaystyle={\bf E}_{1}+{\bf E}_{2}+{\bf E}_{3},$ (78)

where ${\bf
E}_{1}=\mathbb{E}[\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}]$,
${\bf E}_{2}=\mathbb{E}[\mathbf{H}_{\textrm{r}}(\bm{\eta}\odot\bm{\eta})]$ and
${\bf
E}_{3}=\mathbb{E}[-\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}-\mathbf{P}_{\textrm{r}}^{-1}\mathbf{G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}$
$\mathbf{z}_{\textrm{r}}+\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}]$,
the detailed derivations of which are given in Appendix B.

###  IV-B Bias in Stage 2

In this subsection, we aim to derive the bias of $\breve{\bm{\xi}}$ according
to the solution obtained in _Stage 2_. Similar to the derivations of the bias
in _Stage 1_ , we need to derive the expression of estimation error of
${\bm{\xi}}$, and the bias of ${\bm{\xi}}$ can be derived by taking the
expectation of the estimation error of ${\bm{\xi}}$. First, denote
$\tilde{\bm{\xi}}$ as the estimation error of ${\bm{\xi}}$, e.g.,
$\tilde{\bm{\xi}}=\breve{\bm{\xi}}-{\bm{\xi}}$. Then, using the definition of
$\breve{\bm{\xi}}$ in (60), $\tilde{\bm{\xi}}$ can be further derived as

$\tilde{\bm{\xi}}=({\bf G}_{1}^{T}\hat{\bf W}_{1}{\bf G}_{1})^{-1}{\bf
G}^{T}_{1}\hat{\bf W}_{1}(\hat{\bf h}_{1}-{\bf G}_{1}{\bm{\xi}}).$ (79)

According to (54), $\tilde{\bm{\xi}}$ is derived as

$\tilde{\bm{\xi}}=({\bf G}_{1}^{T}\hat{\bf W}_{1}{\bf G}_{1})^{-1}{\bf
G}^{T}_{1}\hat{\bf W}_{1}{\bf z}_{1}.$ (80)

Using the derivations of (III-C2), the definition of ${\bf B}_{1}$ in (58) and
the definition of $\tilde{\mathbf{u}}$ below (54), we have ${\bf z}_{1}={\bf
B}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}}$. Then,
(80) can be derived as

$\tilde{\bm{\xi}}=({\bf G}_{1}^{T}\hat{\bf W}_{1}{\bf G}_{1})^{-1}{\bf
G}_{1}^{T}\hat{\bf
W}_{1}({\mathbf{B}}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}}).$
(81)

By defining $\hat{{\bf P}}_{1}={\bf G}_{1}^{T}\hat{\bf W}_{1}{\bf G}_{1}$, we
have $\tilde{\bm{\xi}}=\hat{{\bf P}}_{1}^{-1}{\bf G}_{1}^{T}\hat{\bf
W}_{1}({\mathbf{B}}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}})$.
According to the definition of $\hat{\bf W}_{1}$ below (59), $\hat{\bf W}_{1}$
is composed of the true value and the estimation error, which are denoted by
${\mathbf{W}}_{1}$ and $\tilde{\mathbf{W}}_{1}$, respectively. To further
derive the expression of $\tilde{\bm{\xi}}$, the expressions of
${\mathbf{W}}_{1}$ and $\tilde{\mathbf{W}}_{1}$ should be derived, which are
given by

$\hat{\mathbf{W}}_{1}={\mathbf{W}}_{1}+\tilde{\mathbf{W}}_{1},{\mathbf{W}}_{1}=\mathbf{B}_{1}^{-1}{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1},\tilde{\mathbf{W}}_{1}=\mathbf{B}_{1}^{-1}\tilde{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1}-{\mathbf{W}}_{1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}-\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}{\mathbf{W}}_{1},$ (82)

where $\tilde{\bf B}_{1}$ is given in (102). The details of deriving
${\mathbf{W}}_{1}$ and $\tilde{\mathbf{W}}_{1}$ are given in Appendix C.
Similarly, based on the definition of $\hat{{\bf P}}_{1}$ below (81),
$\hat{{\bf P}}_{1}$ can be obtained as a summation of the true value and the
estimation error, which are denoted by ${\mathbf{P}}_{1}$ and
$\tilde{\mathbf{P}}_{1}$, respectively. Hence, we have

$\hat{\bf
P}_{1}={\mathbf{P}}_{1}+\tilde{\mathbf{P}}_{1},\quad{\mathbf{P}}_{1}={\bf
G}_{1}^{T}{\mathbf{W}}_{1}{\bf G}_{1},\quad\tilde{\mathbf{P}}_{1}={\bf
G}_{1}^{T}\tilde{\mathbf{W}}_{1}{\bf G}_{1}.$ (83)

Similar to (74), we can derive the approximation of $\hat{\bf P}_{1}^{-1}$ as

$\hat{\bf P}_{1}^{-1}\approx({\bf
I}-{\mathbf{P}}_{1}^{-1}\tilde{\mathbf{P}}_{1}){\mathbf{P}}_{1}^{-1}={\mathbf{P}}_{1}^{-1}-{\mathbf{P}}_{1}^{-1}\tilde{\mathbf{P}}_{1}{\mathbf{P}}_{1}^{-1}.$
(84)

Therefore, using the approximation of $\hat{\bf P}_{1}^{-1}$ in (84) and the
definition of $\hat{\bf W}_{1}$ in (82), (81) can be derived as

$\tilde{\bm{\xi}}\approx({\mathbf{P}}_{1}^{-1}-{\mathbf{P}}_{1}^{-1}\tilde{\mathbf{P}}_{1}{\mathbf{P}}_{1}^{-1}){\bf
G}_{1}^{T}\hat{\bf
W}_{1}({\mathbf{B}}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}}).$
(85)

For notation simplicity, let ${\bf H}_{1}={\mathbf{P}}_{1}^{-1}{\bf
G}^{T}_{1}{\mathbf{W}}_{1}$, we have

$\displaystyle\tilde{\bm{\xi}}$ $\displaystyle\approx{\bf
H}_{1}({\mathbf{B}}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}})+{\mathbf{P}}_{1}^{-1}{\bf
G}^{T}_{1}\tilde{\mathbf{W}}_{1}({\mathbf{B}}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}})-{\mathbf{P}}_{1}^{-1}\tilde{\mathbf{P}}_{1}{\bf
H}_{1}({\mathbf{B}}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}})$
$\displaystyle\quad-{\mathbf{P}}_{1}^{-1}\tilde{\mathbf{P}}_{1}{\mathbf{P}}_{1}^{-1}{\bf
G}_{1}\tilde{\mathbf{W}}_{1}({\mathbf{B}}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}}).$
(86)

By ignoring the error terms higher than the second order and substituting
$\tilde{\mathbf{P}}_{1}$ in (83) into (IV-B), we have

$\tilde{\bm{\xi}}\approx{\bf
H}_{1}({\mathbf{B}}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}})+{\mathbf{P}}_{1}^{-1}{\bf
G}^{T}_{1}\tilde{\mathbf{W}}_{1}{\mathbf{B}}_{1}\tilde{\mathbf{u}}-{\mathbf{P}}_{1}^{-1}\tilde{\mathbf{P}}_{1}{\bf
H}_{1}{\mathbf{B}}_{1}\tilde{\mathbf{u}}.$ (87)

By assuming that ${\bf P}_{2}={\bf I}_{4\times 4}-{\bf G}_{1}{\bf H}_{1}$, we
have

$\tilde{\bm{\xi}}\approx{\bf
H}_{1}({\mathbf{B}}_{1}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}})+{\mathbf{P}}_{1}^{-1}{\bf
G}^{T}_{1}\tilde{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}.$ (88)

Finally, by taking expectation of $\tilde{\bm{\xi}}$, the bias of _Stage 2_ is
given by

$\mathbb{E}[\tilde{\bm{\xi}}]={\bf
H}_{1}{\mathbf{B}}_{1}\mathbb{E}[\tilde{\mathbf{u}}]+{\bf
H}_{1}\mathbb{E}[\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}}]+{\mathbf{P}}_{1}^{-1}{\bf
G}^{T}_{1}\mathbb{E}[\tilde{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}].$ (89)

To further derive the bias of _Stage 2_ , let us introduce Proposition 1 as
follows.

Proposition 1. _For the vector $\tilde{\bf a}\in\mathbb{C}^{M\times 1}$, the
expectation of $\tilde{\bf a}\odot\tilde{\bf a}$ can be expressed as a column
vector ${\bf c}_{\textrm{a}}$ containing the diagonal elements of
${\bm{\Omega}}_{\textrm{a}}$, which is the expectation of the second order
moment of $\tilde{\bf a}$._

Proof: Please see Appendix D. $\hfill\blacksquare$

Using Proposition 1, $\mathbb{E}[\tilde{\mathbf{u}}\odot\tilde{\mathbf{u}}]$
can be derived as ${\bf c}_{\textrm{u}}$ containing the diagonal elements of
${\bm{\Omega}}_{\textrm{u}}$. Furthermore, the details of deriving
${\bm{\Omega}}_{\textrm{u}}$ are given by (100) in Appendix A. Therefore, the
bias in _Stage 2_ can be further derived as

$\mathbb{E}[\tilde{\bm{\xi}}]={\bf
H}_{1}({\mathbf{B}}_{1}\mathbb{E}[\tilde{\mathbf{u}}]+{\bf
c}_{\textrm{u}})+{\mathbf{P}}_{1}^{-1}{\bf
G}^{T}_{1}\mathbb{E}[\tilde{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}],$ (90)

where $\mathbb{E}[\tilde{\mathbf{u}}]$ is given in (IV-A), and the derivations
of $\mathbb{E}[\tilde{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]$ are given in Appendix E.

###  IV-C Bias of Final Solution

First, let us define the estimation error of MU’s position as $\tilde{\bf
q}=\hat{\bf q}-{\bf q}$. Then using the definition of ${\bm{\xi}}$ in (53),
$\breve{\bm{\xi}}$ can be derived as

$\breve{\bm{\xi}}=\left(\hat{{\bf q}}-{\bf p}\right)\odot\left(\hat{{\bf
q}}-{\bf p}\right)=\tilde{\bf q}\odot\tilde{\bf q}+2\tilde{\bf q}\odot({\bf
q}-{\bf p})+{\bm{\xi}}.$ (91)

By utilizing $\tilde{\bm{\xi}}=\breve{\bm{\xi}}-{\bm{\xi}}$, we have
$\tilde{\bm{\xi}}-\tilde{\bf q}\odot\tilde{\bf q}=2\tilde{\bf q}\odot({\bf
q}-{\bf p})$. Then, by assuming that ${\bf B}_{q}\tilde{\bf q}=2\tilde{\bf
q}\odot({\bf q}-{\bf p})$, where ${\bf B}_{q}=2\textrm{diag}\\{({\bf q}-{\bf
p})\\}$ 333The BSs are generally installed on high places, and the MU’s
location is unknown and impossible to overlap with the BS, hence ${\bf B}_{q}$
is a full-rank matrix., we have the expression of $\tilde{\bf q}$ given by

$\tilde{\bf q}={\bf B}_{q}^{-1}(\tilde{\bm{\xi}}-\tilde{\bf q}\odot\tilde{\bf
q}).$ (92)

Then, we can derive the bias of the final solution as

$\mathbb{E}[\tilde{\bf q}]=\mathbb{E}[{\bf
B}_{q}^{-1}(\tilde{\bm{\xi}}-\tilde{\bf q}\odot\tilde{\bf q})].$ (93)

Using Proposition 1, we have $\mathbb{E}[\tilde{\bf q}\odot\tilde{\bf q}]={\bf
c}_{q}$, where ${\bf c}_{q}$ is a column vector formed by the diagonal
elements of ${\bm{\Omega}}_{q}$, and the details of deriving
${\bm{\Omega}}_{q}$ can be found in Appendix F. Therefore, (93) can be further
derived as

$\mathbb{E}[\tilde{\bf q}]={\bf B}_{q}^{-1}(\mathbb{E}[\tilde{\bm{\xi}}]-{\bf
c}_{q}),$ (94)

where $\mathbb{E}[\tilde{\bm{\xi}}]$ is given in (90).

##  V Simulation Results

This section presents simulation results to evaluate the performance of the
proposed positioning algorithm aided by multiple RISs. Moreover, the MU, the
BS, and the RISs are assumed to be placed in a 3D area. The location of the BS
is ${\bf p}=[10,12,12]^{T}$, while the locations of three RISs are ${\bf
s}_{1}=[2,20,2]^{T}$, ${\bf s}_{2}=[-12,-16,58]^{T}$ and ${\bf
s}_{3}=[-10,-10,50]^{T}$, respectively. The phase shift matrix of the RIS is
set to a unit matrix. The following results are obtained by averaging over
10,000 random estimation error realizations. The positioning accuracy is
assessed in terms of the mean-square error (MSE), and the bias. The TDOA and
AOA estimation error are assumed to follow the Gaussian distribution. In the
figures showing the simulation results, we shall use the log scale for the
error level to indicate the wide range of the levels tested.

Figure 5: MSE of proposed positioning algorithm versus $\sigma_{a}$ in case of
2 RISs with $\sigma_{t}=10^{-1},10^{-2},10^{-3},10^{-4}$.

Figure 6:  MSE of proposed positioning algorithm versus $\sigma_{a}$ in case
of 2 RISs and 3 RISs with $\sigma_{t}=10^{-2},10^{-3}$.

Fig. 6 evaluates the impact of the radio positioning parameters’ estimation
error on the performance of the proposed positioning algorithm by comparing
the MSE. To be specific, in the case of deploying 2 RISs in the considered
system, we consider 4 kinds of standard deviation of the TDOA estimation
error, denoted as ‘2 RIS, $\sigma_{t}=10^{-1}$’, ‘2 RIS,
$\sigma_{t}=10^{-2}$’, ‘2 RIS, $\sigma_{t}=10^{-3}$’, and ‘2 RIS,
$\sigma_{t}=10^{-4}$’, respectively. As shown in Fig. 6, when $\sigma_{a}$
decreases from $1$ rad to $10^{-3}$ rad, the MSE decreases, which means that
the accuracy of the proposed algorithm increases with the accuracy of the AOA
estimation. Besides, when $\sigma_{t}$ becomes larger, the MSE of the proposed
positioning algorithm becomes larger as well. It means that the positioning
accuracy deteriorates as the TDOA estimation error increases. This means that
deploying 3 RISs in the proposed systems can significantly improve the
positioning accuracy.

Fig. 6 illustrates the impact of the number of RISs on the performance of the
proposed positioning algorithm. To be specific, we present the comparison of
the positioning system aided by 2 RISs and 3 RISs with $\sigma_{t}=10^{-2}$
and $\sigma_{t}=10^{-3}$. As it can be seen from Fig. 6, MSE increases with
the number of RISs. This implies that the positioning accuracy is better with
a larger number of RISs. Furthermore, in the case of deploying 2 RISs in the
proposed system, the improvement in TDOA estimation accuracy does not provide
significant performance gain in positioning accuracy. However, in the case of
deploying 3 RISs in the proposed system, the improvement of the TDOA
estimation accuracy brings a great improvement to the localization performance
of the proposed algorithm.

Figure 7: MSE of proposed algorithm versus SNR.

Figure 8:  Comparison of proposed algorithm and geometry algorithm.

Fig. 8 compares the MSE of the proposed positioning algorithm aided by 2 RISs
with that aided by 3 RISs when the signal-noise ratio (SNR) increases from
$-30$ dB to $10$ dB. As shown in the figure, the MSE decreases with SNR as
expected. Furthermore, it is shown that the MSE of 2 RISs is larger than that
of 3 RISs. It implies that increasing the number of RISs can significantly
improve the positioning accuracy of the proposed algorithms. Besides, when
deploying 3 RISs in the proposed system with $\sigma_{a}=10^{-2}$, the MSE of
the proposed algorithm can achieve $10^{-10}$ even if SNR is $-30$ dB. This
means that the proposed positioning systems aided by 3 RISs can still have
good positioning performance with low SNR.

In [32], the position of the MU is derived by using the geometry relationship
between the estimated AOA and the 3D position, which is denoted as the
geometry algorithm in this paper. Fig. 8 presents the positioning performance
comparison of the proposed algorithm with the geometry algorithm. It is seen
from Fig. 8 that the proposed algorithm outperforms the geometry algorithm
even with only 2 RISs.

Figure 9: Bias of the proposed algorithm versus $\sigma_{a}$.

Figure 10:  Bias of the proposed algorithm versus $\sigma_{t}$.

Fig. 10 and Fig. 10 confirm the accuracy of the theoretical analysis on the
bias of the proposed positioning algorithm. As shown in Fig. 10, in the case
of 2 RISs, the theory bias is a little larger than the simulation bias when
SNR$=-20,-15$ dB. In the case of 3 RISs, the theoretical bias matches very
well with the simulation bias. As shown in Fig. 10, in the case of 2 RISs, the
theory bias is a little larger than the simulation bias when SNR is larger
than $-10$ dB. In the case of 3 RISs, the theoretical bias matches very well
with simulation results. Therefore, the proposed positioning algorithm aided
by 3 RISs performs better than that aided by 2 RISs.

##  VI Conclusion

In this paper, we studied an mmWave 3D positioning system aided by multiple
RISs. We proposed a TSWLS algorithm to derive the closed-form expression of
the position of the MU. We investigated the bias analysis of the proposed
algorithm. In the first stage, the preliminary estimation of the position of
the MU was obtained using an approximation algorithm. In the second stage, a
more accurate estimation was derived using geometric relationships. Finally,
we investigated the bias analysis to evaluate the performance of the proposed
algorithm. Simulation results confirmed the accuracy of the derived results
and demonstrated the positioning performance of the proposed algorithm.

## Appendix A

First, using the definition of ${\bf z}_{1}$ in (III-C2), ${\bm{\Omega}}_{1}$
is derived as

$\displaystyle{\bm{\Omega}}_{1}$
$\displaystyle=\textrm{diag}\left[2(x_{q}-x_{p}),2(y_{q}-y_{p}),2(z_{q}-z_{p}),2R_{BU}\right]^{T}\cdot$
$\displaystyle\qquad\mathbb{E}\left[[e_{1},e_{2},e_{3},e_{4}]^{T}\left[e_{1},e_{2},e_{3},e_{4}\right]\right]\textrm{diag}\left[2(x_{q}-x_{p}),2(y_{q}-y_{p}),2(z_{q}-z_{p}),2R_{BU}\right]^{T}$
$\displaystyle={\bf B}_{1}\cdot{\bm{\Omega}}_{\textrm{u}}\cdot{\bf B}_{1},$
(95)

where ${\bm{\Omega}}_{\textrm{u}}$ denotes the covariance matrix of ${\bf u}$.
Here, ${\bf B}_{1}$ is given in (58), and ${\bm{\Omega}}_{\textrm{u}}$ is
defined as follows

$\displaystyle{\bm{\Omega}}_{\textrm{u}}=\mathbb{E}\left[\left[e_{1},e_{2},e_{3},e_{4}\right]^{T}\left[e_{1},e_{2},e_{3},e_{4}\right]\right].$
(96)

Then, we need to derive the expression of ${\bm{\Omega}}_{\textrm{u}}$ as
follows. First, we will derive the error of $\mathbf{u}$, which is denoted as
$\tilde{\mathbf{u}}$. According to (49) and the definition of $\mathbf{u}$,
$\tilde{\mathbf{u}}$ is formulated as

$\tilde{\mathbf{u}}=\breve{\bf u}-{\bf
u}=(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1}\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}(\hat{\mathbf{h}}_{\textrm{r}}-\hat{\mathbf{G}}_{\textrm{r}}\mathbf{u}).$
(97)

According to (42), we have

$\tilde{\mathbf{u}}=(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1}\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}.$
(98)

Using the expression of $\tilde{\mathbf{u}}$ and taking expectation of
$\tilde{\mathbf{u}}\tilde{\mathbf{u}}^{T}$, ${\bm{\Omega}}_{\textrm{u}}$ is
given by

${\bm{\Omega}}_{\textrm{u}}=(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1}\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{Q}_{\textrm{r}}\mathbf{B}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}^{T}\hat{\mathbf{G}}_{\textrm{r}}((\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1})^{T}.$
(99)

Note that
$\mathbf{W}_{\textrm{r}}=(\mathbf{B}_{\textrm{r}}\mathbf{Q}_{\textrm{r}}\mathbf{B}_{\textrm{r}}^{T})^{-1}$
according to (50) and $\mathbf{W}_{\textrm{r}}$ is a diagonal matrix. Thus, we
have $\mathbf{W}_{\textrm{r}}^{T}=\mathbf{W}_{\textrm{r}}$. As a result, (99)
can be derived as

${\bm{\Omega}}_{\textrm{u}}=(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1}.$
(100)

As a result, ${\bm{\Omega}}_{1}$ is derived as

${\bm{\Omega}}_{1}={\bf
B}_{1}(\hat{\mathbf{G}}_{\textrm{r}}^{T}\mathbf{W}_{\textrm{r}}\hat{\mathbf{G}}_{\textrm{r}})^{-1}{\bf
B}_{1}.$ (101)

## Appendix B

By utilizing (46), recall that
$\mathbf{z}_{\textrm{r}}=[\mathbf{n}^{T},\bm{\omega}^{T},\bm{\nu}^{T}]^{T}$,
$\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}$ can be expanded as
$(\mathbf{B}_{{\textrm{r},n}}\mathbf{n}+\mathbf{B}_{{\textrm{r},}\omega}\bm{\omega}+\mathbf{B}_{{\textrm{r},}\nu}\bm{\nu})$.
Then, we have

${\bf
E}_{1}=\mathbb{E}[\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}]=\mathbf{H}_{\textrm{r}}\mathbb{E}[\mathbf{B}_{{\textrm{r},n}}\mathbf{n}+\mathbf{B}_{{\textrm{r},}\omega}\bm{\omega}+\mathbf{B}_{{\textrm{r},}\nu}\bm{\nu}]=\mathbf{0}.$

For ${\bf E}_{2}$, we have

${\bf
E}_{2}=\mathbb{E}[\mathbf{H}_{\textrm{r}}(\bm{\eta}\odot\bm{\eta})]=\frac{1}{2}\mathbf{H}_{\textrm{r}}\mathbb{E}[[\underbrace{0,\cdots,0}_{2M},\nu^{2}_{1},\cdots,\nu^{2}_{M}]^{T}]=\frac{1}{2}\mathbf{H}_{\textrm{r}}\mathbf{q}_{u},$

where

$\mathbf{q}_{u}=\frac{1}{2}[\underbrace{0,\cdots,0}_{2M},\sigma^{2}_{\nu_{1}},\cdots,\sigma^{2}_{\nu_{M}}]^{T}.$

Then, ${\bf E}_{3}$ can be derived as

$\displaystyle{\bf E}_{3}$
$\displaystyle=\mathbb{E}[\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}]-\mathbb{E}[\mathbf{P}_{\textrm{r}}^{-1}\tilde{\mathbf{G}}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}]-\mathbb{E}[\mathbf{P}_{\textrm{r}}^{-1}\mathbf{G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}\tilde{\mathbf{G}}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}]$
$\displaystyle={\bf E}_{31}-{\bf E}_{32}-{\bf E}_{33},$

where the detailed derivations of ${\bf E}_{31}$, ${\bf E}_{32}$ and ${\bf
E}_{33}$ are given as follows.

For ${\bf E}_{31}$, as we have $\tilde{\bf G}_{\textrm{r}}$ in (IV) and
$\mathbf{B}_{\textrm{r}}\mathbf{z}_{\textrm{r}}=(\mathbf{B}_{{\textrm{r},n}}\mathbf{n}+\mathbf{B}_{{\textrm{r},}\omega}\bm{\omega}+\mathbf{B}_{{\textrm{r},}\nu}\bm{\nu})$,
${\bf E}_{31}$ can be derived as

$\displaystyle{\bf E}_{31}$
$\displaystyle=\mathbb{E}[\mathbf{P}_{\textrm{r}}^{-1}(\bar{\mathbf{G}}_{n}^{T}\odot[(\mathbf{n1}_{1\times
4})^{T},(\mathbf{n1}_{1\times 4})^{T},(\mathbf{n1}_{1\times
4})^{T}])\mathbf{W}_{\textrm{r}}\mathbf{B}_{{\textrm{r},n}}\textrm{diag}({\bf
n}){\bf 1}_{M\times 1}]$
$\displaystyle\quad+\mathbb{E}[\mathbf{P}_{\textrm{r}}^{-1}(\bar{\mathbf{G}}^{T}_{\omega}\odot[(\bm{\omega
1}_{1\times 4})^{T},(\bm{\omega 1}_{1\times 4})^{T},(\bm{\omega 1}_{1\times
4})^{T}])\mathbf{W}_{\textrm{r}}\mathbf{B}_{{\textrm{r},\omega}}\textrm{diag}({\bm{\omega}}){\bf
1}_{M\times 1}]$
$\displaystyle\quad+\mathbb{E}[\mathbf{P}_{\textrm{r}}^{-1}(\bar{\mathbf{G}}^{T}_{\nu}\odot[(\bm{\nu
1}_{1\times 4})^{T},(\bm{\nu 1}_{1\times 4})^{T},(\bm{\nu 1}_{1\times
4})^{T}])\mathbf{W}_{\textrm{r}}\mathbf{B}_{{\textrm{r},\nu}}\textrm{diag}({\bm{\nu}}){\bf
1}_{M\times 1}].$

To further derive the expression of ${\bf E}_{31}$, we introduce Proposition 2
and Proposition 3 as follows.

Proposition 2. _For vectors_ ${\bm{\alpha}}\in\mathbb{C}^{M\times 1}$,
${\bm{\beta}}\in\mathbb{C}^{N\times 1}$ _and matrix_ ${\bf
G}\in\mathbb{C}^{M\times N}$, _and corresponding diagonal matrices_
$\textrm{diag}({\bm{\alpha}})$ _and_ $\textrm{diag}({\bm{\beta}})$ _with these
vectors as their main diagonals, we have_
$({\bm{\alpha}}{\bm{\beta}}^{T})\odot{\bf G}=\textrm{diag}({\bm{\alpha}}){\bf
G}\textrm{diag}({\bm{\beta}})$.

Proof: Please see reference [33].$\hfill\blacksquare$

Proposition 3. _For a vector ${\bf x}\in\mathbb{C}^{M\times 1}$ with its
covariance matrix given by ${\bf Q}_{x}$, we have_

$\mathbb{E}[({\bf A}\odot[({\bf x1}_{1\times N})^{T},({\bf x1}_{1\times
N})^{T},({\bf x1}_{1\times N})^{T}]){\bf B}\textrm{diag}({\bf x})]={\bf
A}({\bf B}\odot[{\bf Q}_{x},{\bf Q}_{x},{\bf Q}_{x}]^{T}),$

_where ${\bf A}\in\mathbb{C}^{N\times 3M}$ and ${\bf B}\in\mathbb{C}^{3M\times
M}$._

Proof: First, the expression on the left hand side of the equation can be
derived as

$\displaystyle\mathbb{E}[({\bf A}\odot[({\bf x1}_{1\times N})^{T},({\bf
x1}_{1\times N})^{T},({\bf x1}_{1\times N})^{T}]){\bf B}\textrm{diag}({\bf
x})]$ $\displaystyle=\mathbb{E}[({\bf A}\odot({\bf 1}_{N\times 1}({\bf
1}_{1\times 3}\otimes{\bf x}^{T}))){\bf B}\textrm{diag}({\bf x})]$
$\displaystyle=\mathbb{E}[(({\bf 1}_{N\times 1}({\bf 1}_{1\times 3}\otimes{\bf
x}^{T}))\odot{\bf A}){\bf B}\textrm{diag}({\bf x})]$
$\displaystyle=\mathbb{E}[(({\bf 1}_{N\times 1}({\bf 1}_{3\times 1}\otimes{\bf
x})^{T})\odot{\bf A}){\bf B}\textrm{diag}({\bf x})].$

Then, using Proposition 2, we have

$\displaystyle\mathbb{E}[(({\bf 1}_{N\times 1}({\bf 1}_{3\times 1}\otimes{\bf
x})^{T})\odot{\bf A}){\bf B}\textrm{diag}({\bf x})]$
$\displaystyle=\mathbb{E}[{\bf A}\textrm{diag}({\bf 1}_{3\times 1}\otimes{\bf
x}){\bf B}\textrm{diag}({\bf x})]$ $\displaystyle=\mathbb{E}[{\bf A}([({\bf
1}_{3\times 1}\otimes{\bf x}){\bf x}^{T}]\odot{\bf B})]$
$\displaystyle=\mathbb{E}[{\bf A}(({\bf 1}_{1\times 3}\otimes{\bf
xx}^{T})^{T}\odot{\bf B})].$

Then, by calculating the expectation of ${\bf xx}^{T}$, the expression can be
rewritten as

$\mathbb{E}[{\bf A}(({\bf 1}_{1\times 3}\otimes{\bf xx}^{T})^{T}\odot{\bf
B})]={\bf A}(({\bf 1}_{1\times 3}\otimes{\bf Q}_{x})^{T}\odot{\bf B})={\bf
A}({\bf B}\odot[{\bf Q}_{x},{\bf Q}_{x},{\bf Q}_{x}]^{T}),$

which completes the proof of Proposition 3. $\hfill\blacksquare$

Consequently, using Proposition 3, ${\bf E}_{31}$ can be derived as

$\displaystyle{\bf E}_{31}$
$\displaystyle=\mathbf{P}_{\textrm{r}}^{-1}\bar{\mathbf{G}}_{n}^{T}((\mathbf{W}_{\textrm{r}}\mathbf{B}_{{\textrm{r},n}})\odot[{\bf
Q}_{n},{\bf Q}_{n},{\bf Q}_{n}]^{T}){\bf 1}_{M\times
1}+\mathbf{P}_{\textrm{r}}^{-1}\bar{\mathbf{G}}^{T}_{\omega}((\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r},\omega})\odot[{\bf
Q}_{\omega},{\bf Q}_{\omega},{\bf Q}_{\omega}]^{T}){\bf 1}_{M\times 1}$
$\displaystyle\quad+\mathbf{P}_{\textrm{r}}^{-1}\bar{\mathbf{G}}^{T}_{\nu}((\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r},\nu})\odot[{\bf
Q}_{\nu},{\bf Q}_{\nu},{\bf Q}_{\nu}]^{T}){\bf 1}_{M\times 1}$
$\displaystyle=\mathbf{P}_{\textrm{r}}^{-1}(\bar{\mathbf{G}}_{n}^{T}((\mathbf{W}_{\textrm{r}}\mathbf{B}_{{\textrm{r},n}})\odot\bar{\bf
Q}_{n}){\bf 1}_{M\times
1}+\bar{\mathbf{G}}^{T}_{\omega}((\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r},\omega})\odot\bar{\bf
Q}_{\omega}){\bf 1}_{M\times
1}+\bar{\mathbf{G}}^{T}_{\nu}((\mathbf{W}_{\textrm{r}}\mathbf{B}_{\textrm{r},\nu})\odot\bar{\mathbf{Q}}_{\nu}){\bf
1}_{M\times 1}),$

where

$\displaystyle\bar{\mathbf{Q}}_{n}=[1,1,1]^{T}\otimes\mathbf{Q}_{n},\qquad\bar{\mathbf{Q}}_{\omega}=[1,1,1]^{T}\otimes\mathbf{Q}_{\omega},\qquad\bar{\mathbf{Q}}_{\nu}=[1,1,1]^{T}\otimes\mathbf{Q}_{\nu}.$

Similarly, by applying Proposition 2 and Proposition 3, ${\bf E}_{32}$ can be
derived as

$\displaystyle{\bf E}_{32}$
$\displaystyle=\mathbb{E}[\mathbf{P}_{\textrm{r}}^{-1}(\bar{\mathbf{G}}_{n}^{T}\odot[(\mathbf{n1}_{1\times
4})^{T},(\mathbf{n1}_{1\times 4})^{T},(\mathbf{n1}_{1\times
4})^{T}]+\bar{\mathbf{G}}^{T}_{\omega}\odot[(\bm{\omega 1}_{1\times
4})^{T},(\bm{\omega 1}_{1\times 4})^{T},(\bm{\omega 1}_{1\times 4})^{T}]$
$\displaystyle\quad+\bar{\mathbf{G}}^{T}_{\nu}\odot[(\bm{\nu 1}_{1\times
4})^{T},(\bm{\nu 1}_{1\times 4})^{T},(\bm{\nu 1}_{1\times
4})^{T}])\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}\mathbf{H}_{\textrm{r}}(\mathbf{B}_{{\textrm{r},n}}\mathbf{n}+\mathbf{B}_{{\textrm{r},}\omega}\bm{\omega}+\mathbf{B}_{{\textrm{r},}\nu}\bm{\nu})]$
$\displaystyle=\mathbf{P}_{\textrm{r}}^{-1}(\bar{\mathbf{G}}_{n}^{T}((\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{{\textrm{r},n}})\odot\bar{\bf
Q}_{n}){\bf 1}_{M\times
1}+\bar{\mathbf{G}}^{T}_{\omega}((\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r},\omega})\odot\bar{\bf
Q}_{\omega}){\bf 1}_{M\times 1}$
$\displaystyle\quad+\bar{\mathbf{G}}^{T}_{\nu}((\mathbf{W}_{\textrm{r}}\mathbf{G}_{\textrm{r}}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r},\nu})\odot\bar{\mathbf{Q}}_{\nu}){\bf
1}_{M\times 1}).$

Similarly, using Proposition 2 and Proposition 3, ${\bf E}_{33}$ is formulated
as

$\displaystyle{\bf E}_{33}$
$\displaystyle=\mathbb{E}[\mathbf{P}_{\textrm{r}}^{-1}{\bf
G}^{T}_{\textrm{r}}\mathbf{W}_{\textrm{r}}(\bar{\mathbf{G}}_{n}\odot[(\mathbf{n1}_{1\times
4})^{T},(\mathbf{n1}_{1\times 4})^{T},(\mathbf{n1}_{1\times
4})^{T}]^{T}+\bar{\mathbf{G}}_{\omega}\odot[(\bm{\omega 1}_{1\times
4})^{T},(\bm{\omega 1}_{1\times 4})^{T},(\bm{\omega 1}_{1\times 4})^{T}]^{T}$
$\displaystyle\quad+\bar{\mathbf{G}}_{\nu}\odot[(\bm{\nu 1}_{1\times
4})^{T},(\bm{\nu 1}_{1\times 4})^{T},(\bm{\nu 1}_{1\times
4})^{T}]^{T})\mathbf{H}_{\textrm{r}}(\mathbf{B}_{{\textrm{r},n}}\mathbf{n}+\mathbf{B}_{{\textrm{r},}\omega}\bm{\omega}+\mathbf{B}_{{\textrm{r},}\nu}\bm{\nu})]$
$\displaystyle=\mathbf{P}_{\textrm{r}}^{-1}{\bf
G}^{T}_{\textrm{r}}(((\mathbf{W}_{\textrm{r}}\bar{\mathbf{G}}_{n}\mathbf{H}_{\textrm{r}}\mathbf{B}_{{\textrm{r},n}})\odot\bar{\bf
Q}_{n}){\bf 1}_{M\times
1}+((\mathbf{W}_{\textrm{r}}\bar{\mathbf{G}}_{\omega}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r},\omega})\odot\bar{\bf
Q}_{\omega}){\bf 1}_{M\times 1}$
$\displaystyle\quad+((\mathbf{W}_{\textrm{r}}\bar{\mathbf{G}}_{\nu}\mathbf{H}_{\textrm{r}}\mathbf{B}_{\textrm{r},\nu})\odot\bar{\mathbf{Q}}_{\nu}){\bf
1}_{M\times 1}).$

## Appendix C

First, based on the definition of $\hat{\mathbf{W}}_{1}$ below (59) and the
definition of $\hat{\bf P}_{\textrm{r}}$ above (71), we have
$\hat{\mathbf{W}}_{1}=\hat{\mathbf{B}}_{1}^{-1}\hat{\bf
P}_{\textrm{r}}\hat{\mathbf{B}}_{1}^{-1}$. Then, according to the definition
of ${\mathbf{B}}_{1}$ in (58) and $\hat{\mathbf{B}}_{1}$ in (59), we have

$\hat{\bf B}_{1}=\mathbf{B}_{1}+\tilde{\bf B}_{1},\tilde{\bf
B}_{1}=2\textrm{diag}\\{\tilde{\bf u}\\}.$ (102)

If we use $\hat{\bf B}_{1}=\mathbf{B}_{1}+\tilde{\bf B}_{1}$ to derive
$\hat{\mathbf{W}}_{1}$, the inverse matrix $\hat{\bf
B}_{1}^{-1}=(\mathbf{B}_{1}+\tilde{\bf B}_{1})^{-1}$ would be complex and
challenging. Fortunately, using the Newmann expansion [28], we can derive the
approximation of $\hat{\bf
B}_{1}^{-1}\approx\mathbf{B}_{1}^{-1}-\mathbf{B}_{1}^{-1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}$. As a result, $\hat{\bf W}_{1}$ is given by

$\hat{\bf W}_{1}=(\mathbf{B}_{1}^{-1}-\mathbf{B}_{1}^{-1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1})({\bf P}_{\textrm{r}}+\tilde{\bf
P}_{\textrm{r}})(\mathbf{B}_{1}^{-1}-\mathbf{B}_{1}^{-1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}).$ (103)

By neglecting the error terms that are higher than the first order, $\hat{\bf
W}_{1}$ can be approximated as

$\hat{\bf W}_{1}\approx\mathbf{B}_{1}^{-1}{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1}+\mathbf{B}_{1}^{-1}\tilde{{\bf
P}}_{\textrm{r}}\mathbf{B}_{1}^{-1}-\mathbf{B}_{1}^{-1}{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}-\mathbf{B}_{1}^{-1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1}={\mathbf{W}}_{1}+\tilde{\mathbf{W}}_{1},$
(104)

where ${\bf W}_{1}$ and $\tilde{\mathbf{W}}_{1}$ are defined as

${\mathbf{W}}_{1}=\mathbf{B}_{1}^{-1}{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1},\tilde{\mathbf{W}}_{1}=\mathbf{B}_{1}^{-1}\tilde{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1}-{\mathbf{W}}_{1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}-\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}{\mathbf{W}}_{1}.$ (105)

## Appendix D

Before proving Proposition 1, let us introduce Proposition 4 as follows.

Proposition 4. _The Hadamard product of two vectors_ ${\bf
a}\in\mathbb{C}^{N\times 1}$ _and_ ${\bf b}\in\mathbb{C}^{N\times 1}$, _is the
same as matrix multiplication of one vector by the corresponding diagonal
matrix of the other vector:_

${\bf a}\odot{\bf b}=\textrm{diag}({\bf a}){\bf b}.$

Proof: please see reference [33]. $\hfill\blacksquare$

Let $\mathbb{E}[\tilde{\bf a}\odot\tilde{\bf a}]$ denote the expectation of
$\tilde{\bf a}\odot\tilde{\bf a}$, using Proposition 4, we have
$\mathbb{E}[\tilde{\bf a}\odot\tilde{\bf
a}]=\mathbb{E}[\textrm{diag}(\tilde{\bf a})\tilde{\bf
a}]=\mathbb{E}[\textrm{diag}(\tilde{\bf a})\textrm{diag}(\tilde{\bf a}){\bf
1}_{M\times 1}]$. Using Proposition 2, $\mathbb{E}[\tilde{\bf
a}\odot\tilde{\bf a}]$ can be rewritten as

$\mathbb{E}[\tilde{\bf a}\odot\tilde{\bf a}]=\mathbb{E}[((\tilde{\bf
a}\tilde{\bf a}^{T})\odot{\bf I}_{M\times M}){\bf 1}_{M\times
1}]=({\bm{\Omega}}_{\textrm{a}}\odot{\bf I}_{M\times M}){\bf 1}_{M\times
1}={\bf c}_{\textrm{a}},$

where ${\bf c}_{\textrm{a}}$ denotes the vector containing the diagonal
elements of ${\bm{\Omega}}_{\textrm{a}}$. Here, we complete the proof of
Proposition 1.

## Appendix E

Using the definition of $\tilde{\mathbf{W}}_{1}$ in (105),
$\mathbb{E}[\tilde{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]$ can be derived as

$\displaystyle\mathbb{E}[\tilde{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]$
$\displaystyle=\mathbb{E}[(\mathbf{B}_{1}^{-1}\tilde{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1}-{\mathbf{W}}_{1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}-\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}{\mathbf{W}}_{1}){\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]$
$\displaystyle=\mathbb{E}[\mathbf{B}_{1}^{-1}\tilde{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]-\mathbb{E}[{\mathbf{W}}_{1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]-\mathbb{E}[\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}].$

Let ${\mathbf{a}}_{1}=\mathbb{E}[\mathbf{B}_{1}^{-1}\tilde{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]$,
${\mathbf{a}}_{2}=-\mathbb{E}[{\mathbf{W}}_{1}\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}{\bf P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]$ and
${\mathbf{a}}_{3}=-\mathbb{E}[\tilde{\bf
B}_{1}\mathbf{B}_{1}^{-1}{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]$, we have
$\mathbb{E}[\tilde{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]={\mathbf{a}}_{1}+{\mathbf{a}}_{2}+{\mathbf{a}}_{3}$,
while ${\mathbf{a}}_{1}$, ${\mathbf{a}}_{2}$, ${\mathbf{a}}_{3}$ are given as
follows.

Using the definition of $\tilde{\bf P}_{\textrm{r}}$ below (73),
${\mathbf{a}}_{1}$ can be derived as

$\displaystyle{\mathbf{a}}_{1}$
$\displaystyle=\mathbb{E}[\mathbf{B}_{1}^{-1}\tilde{\bf
P}_{\textrm{r}}\mathbf{B}_{1}^{-1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]=\mathbb{E}[\mathbf{B}_{1}^{-1}(\tilde{\bf
G}^{T}_{\textrm{r}}{\bf W}_{\textrm{r}}{\bf G}_{\textrm{r}}+{\bf
G}^{T}_{\textrm{r}}{\bf W}_{\textrm{r}}\tilde{\bf
G}_{\textrm{r}})\mathbf{B}_{1}^{-1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}].$

Using the definition of $\tilde{\mathbf{u}}$ in (77) and ignoring the error
terms that are higher than the second order, ${\mathbf{a}}_{1}$ can be
approximated as $\mathbb{E}[\mathbf{B}_{1}^{-1}(\tilde{\bf
G}^{T}_{\textrm{r}}{\bf W}_{\textrm{r}}{\bf G}_{\textrm{r}}+{\bf
G}^{T}_{\textrm{r}}{\bf W}_{\textrm{r}}\tilde{\bf
G}_{\textrm{r}})\mathbf{B}_{1}^{-1}{\bf P}_{2}{\mathbf{B}}_{1}{\bf
H}_{\textrm{r}}{\bf B}_{\textrm{r}}{\bf z}_{\textrm{r}}]$. For notation
simplicity, we assume that $\mathbf{B}_{1}^{-1}{\bf P}_{2}{\mathbf{B}}_{1}{\bf
H}_{\textrm{r}}={\bf P}_{3}$. Then, by applying Proposition 2 and Proposition
3, ${\mathbf{a}}_{1}$ can be derived as

$\displaystyle{\mathbf{a}}_{1}$ $\displaystyle=\mathbf{B}_{1}^{-1}(\bar{{\bf
G}}^{T}_{n}(({\bf W}_{\textrm{r}}{\bf G}_{\textrm{r}}{\bf P}_{3}{\bf
B}_{\textrm{r},n})\odot\bar{{\bf Q}}_{n}){\bf 1}_{M\times 1}+\bar{{\bf
G}}^{T}_{\omega}(({\bf W}_{\textrm{r}}{\bf G}_{\textrm{r}}{\bf P}_{3}{\bf
B}_{\textrm{r},\omega})\odot\bar{{\bf Q}}_{\omega}){\bf 1}_{M\times 1}$
$\displaystyle\quad+\bar{{\bf G}}^{T}_{\nu}(({\bf W}_{\textrm{r}}{\bf
G}_{\textrm{r}}{\bf P}_{3}{\bf B}_{\textrm{r},\nu})\odot\bar{{\bf
Q}}_{\nu}){\bf 1}_{M\times 1}+{\bf G}^{T}_{\textrm{r}}(({\bf
W}_{\textrm{r}}\bar{{\bf G}}_{n}{\bf P}_{3}{\bf
B}_{\textrm{r},n})\odot\bar{{\bf Q}}_{n}){\bf 1}_{M\times 1}$
$\displaystyle\quad+{\bf G}^{T}_{\textrm{r}}(({\bf W}_{\textrm{r}}\bar{{\bf
G}}_{\omega}{\bf P}_{3}{\bf B}_{\textrm{r},\omega})\odot\bar{{\bf
Q}}_{\omega}){\bf 1}_{M\times 1}+{\bf G}^{T}_{\textrm{r}}(({\bf
W}_{\textrm{r}}\bar{{\bf G}}_{\nu}{\bf P}_{3}{\bf
B}_{\textrm{r},\nu})\odot\bar{{\bf Q}}_{\nu}){\bf 1}_{M\times 1}).$

Moreover, for ${\mathbf{a}}_{2}$, using the definition of $\tilde{\bf B}_{1}$
in (102), we have

$\displaystyle{\mathbf{a}}_{2}$
$\displaystyle=-\mathbb{E}[{\mathbf{W}}_{1}\mathbf{B}_{1}^{-1}\tilde{\bf
B}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]=-2{\mathbf{W}}_{1}\mathbf{B}_{1}^{-1}\mathbb{E}[\textrm{diag}(\tilde{\bf
u}){\bf P}_{2}{\mathbf{B}}_{1}\textrm{diag}(\tilde{\bf u}){\bf 1}_{4\times
1}].$

Using Proposition 2, we have

$\displaystyle{\mathbf{a}}_{2}=-2{\mathbf{W}}_{1}\mathbf{B}_{1}^{-1}\mathbb{E}[(\tilde{\bf
u}\tilde{\bf u}^{T})\odot({\bf P}_{2}{\bf B}_{1})]{\bf 1}_{4\times
1}=-2{\mathbf{W}}_{1}\mathbf{B}_{1}^{-1}({\bf\Omega}_{\textrm{u}}\odot({\bf
P}_{2}{\bf B}_{1})){\bf 1}_{4\times 1}.$

As ${\bf\Omega_{\textrm{u}}}$ is a diagonal matrix, ${\mathbf{a}}_{2}$ can be
further derived as

$\displaystyle{\mathbf{a}}_{2}$
$\displaystyle=-2{\mathbf{W}}_{1}\mathbf{B}_{1}^{-1}\textrm{diag}({\bf\Omega}_{\textrm{u}}({\bf
P}_{2}{\bf B}_{1})){\bf 1}_{4\times
1}=-2{\mathbf{W}}_{1}\mathbf{B}_{1}^{-1}{\bf p}_{\textrm{4}},$

where ${\bf p}_{\textrm{4}}$ represents the column vector formed by the
diagonal elements of ${\bf P}_{4}={\bf
P}_{2}{\mathbf{B}}_{1}{\bf\Omega}_{\textrm{u}}$.

Similar to the derivations of ${\bf a}_{2}$, ${\mathbf{a}}_{3}$ can be derived
as

$\displaystyle{\mathbf{a}}_{3}=-\mathbb{E}[\mathbf{B}_{1}^{-1}\tilde{\bf
B}_{1}{\mathbf{W}}_{1}{\bf
P}_{2}{\mathbf{B}}_{1}\tilde{\mathbf{u}}]=-2\mathbf{B}_{1}^{-1}{\bf
p}_{\textrm{4,W1}},$

where ${\bf p}_{\textrm{4,W1}}$ represents the column vector formed by the
diagonal elements of ${\mathbf{W}}_{1}{\bf P}_{4}$.

## Appendix F

First, let us derive the covariance matrix of ${\bm{\xi}}$, which is denoted
by ${\bf\Omega}_{\xi}$. Using the definition of $\tilde{\bm{\xi}}$ in (80),
its covariance matrix ${\bf\Omega}_{\xi}$ is given by

$\displaystyle{\bf\Omega}_{\xi}$
$\displaystyle=\mathbb{E}[\tilde{\bm{\xi}}\tilde{\bm{\xi}}^{T}]\approx{\bf
G}_{1}^{T}\hat{\bf W}_{1}{\bf G}_{1})^{-1}{\bf G}^{T}_{1}\hat{\bf
W}_{1}\hat{\bf W}_{1}^{-1}\hat{\bf W}_{1}^{T}{\bf G}_{1}(({\bf
G}_{1}^{T}\hat{\bf W}_{1}{\bf G}_{1})^{-1})^{T}=({\bf G}_{1}^{T}\hat{\bf
W}_{1}{\bf G}_{1})^{-1},$

where $\mathbb{E}({\bf z}_{1}{\bf
z}_{1}^{T})\approx\hat{\bf\Omega}_{1}=\hat{\bf W}_{1}^{-1}$ is given below
(59).

Furthermore, by neglecting the second order term in (92), the $\tilde{\bf q}$
can be approximated as ${\bf B}_{q}^{-1}\tilde{\bm{\xi}}$, therefore, the
covariance matrix of $\hat{\bf q}$ can be derived as

$\displaystyle{\bf\Omega}_{q}$ $\displaystyle=\mathbb{E}[\tilde{\bf
q}\tilde{\bf q}^{T}]={\bf B}_{q}^{-1}({\bf G}_{1}^{T}\hat{\bf W}_{1}{\bf
G}_{1})^{-1}{\bf B}_{q}^{-1}.$

## References

  * [1] C. L. Nguyen, O. Georgiou, G. Gradoni, and M. Di Renzo, “Wireless Fingerprinting Localization in Smart Environments Using Reconfigurable Intelligent Surfaces,” _IEEE Access_ , vol. 9, pp. 135 526–135 541, 2021.
  * [2] H. Ren, K. Wang, and C. Pan, “Intelligent Reflecting Surface-Aided URLLC in a Factory Automation Scenario,” _IEEE Trans. Commun._ , vol. 70, no. 1, pp. 707–723, 2022.
  * [3] S. Han, C.-L. I., T. Xie, S. Wang, Y. Huang, L. Dai, Q. Sun, and C. Cui, “Achieving High Spectrum Efficiency on High Speed Train for 5G New Radio and Beyond,” _IEEE Wireless Commun._ , vol. 26, no. 5, pp. 62–69, 2019.
  * [4] R. Liu, Q. Wu, M. Di Renzo, and Y. Yuan, “A Path to Smart Radio Environments: An Industrial Viewpoint on Reconfigurable Intelligent Surfaces,” _IEEE Wireless Commun._ , vol. 29, no. 1, pp. 202–208, 2022.
  * [5] M. Di Renzo, A. Zappone, M. Debbah, M.-S. Alouini, C. Yuen, J. de Rosny, and S. Tretyakov, “Smart Radio Environments Empowered by Reconfigurable Intelligent Surfaces: How It Works, State of Research, and The Road Ahead,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2450–2525, 2020.
  * [6] C. Pan, G. Zhou, K. Zhi, S. Hong, T. Wu, Y. Pan, H. Ren, M. D. Renzo, A. Lee Swindlehurst, R. Zhang, and A. Y. Zhang, “An Overview of Signal Processing Techniques for RIS/IRS-Aided Wireless Systems,” _IEEE J. Sel. Topics Signal Process._ , vol. 16, no. 5, pp. 883–917, 2022.
  * [7] Y. Pan, K. Wang, C. Pan, H. Zhu, and J. Wang, “Self-Sustainable Reconfigurable Intelligent Surface Aided Simultaneous Terahertz Information and Power Transfer (STIPT),” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2022.
  * [8] H. Ren, K. Wang, and C. Pan, “Intelligent Reflecting Surface-aided URLLC in a Factory Automation Scenario,” _IEEE Trans. Commun._ , pp. 1–1, 2021.
  * [9] F. Jiang and A. L. Swindlehurst, “Optimization of UAV Heading for the Ground-to-Air Uplink,” _IEEE J. Sel. Areas Commun._ , vol. 30, no. 5, pp. 993–1005, 2012.
  * [10] K. Zhi, C. Pan, H. Ren, and K. Wang, “Statistical CSI-Based Design for Reconfigurable Intelligent Surface-Aided Massive MIMO Systems With Direct Links,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 5, pp. 1128–1132, 2021.
  * [11] M. Renzo, M. Debbah, D. T. Phan-Huy, A. Zappone, M. S. Alouini, C. Yuen, V. Sciancalepore, G. C. Alexandropoulos, J. Hoydis, and H. Gacanin, “Smart Radio Environments Empowered by AI Reconfigurable Meta-Surfaces: An Idea Whose Time Has Come,” _EURASIP J Wirel Commun Netw_ , vol. 2019, no. 1, 2019.
  * [12] C. Pan, H. Ren, K. Wang, W. Xu, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Multicell MIMO Communications Relying on Intelligent Reflecting Surfaces,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 8, pp. 5218–5233, 2020.
  * [13] G. Zhou, C. Pan, H. Ren, K. Wang, and Z. Peng, “Secure Wireless Communication in RIS-Aided MISO System With Hardware Impairments,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 6, pp. 1309–1313, 2021.
  * [14] G. Zhou, C. Pan, H. Ren, K. Wang, and A. Nallanathan, “A Framework of Robust Transmission Design for IRS-Aided MISO Communications With Imperfect Cascaded Channels,” _IEEE Trans. Signal Process._ , vol. 68, pp. 5092–5106, 2020.
  * [15] K. Zhi, C. Pan, H. Ren, and K. Wang, “Uplink Achievable Rate of Intelligent Reflecting Surface-Aided Millimeter-Wave Communications With Low-Resolution ADC and Phase Noise,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 3, pp. 654–658, 2021.
  * [16] Z. Zhang and L. Dai, “A Joint Precoding Framework for Wideband Reconfigurable Intelligent Surface-Aided Cell-Free Network,” _IEEE Trans. Signal Process._ , vol. 69, pp. 4085–4101, 2021.
  * [17] D. Fan, F. Gao, Y. Liu, Y. Deng, G. Wang, Z. Zhong, and A. Nallanathan, “Angle Domain Channel Estimation in Hybrid Millimeter Wave Massive MIMO Systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 12, pp. 8165–8179, 2018.
  * [18] G. Zhou, C. Pan, H. Ren, P. Popovski, and A. L. Swindlehurst, “Channel estimation for RIS-aided multiuser millimeter-wave systems,” _IEEE Trans. Signal Process._ , vol. 70, pp. 1478–1492, 2022.
  * [19] J. He, H. Wymeersch, L. Kong, O. Silv$\acute{\mathrm{e}}$n, and M. Juntti, “Large intelligent surface for positioning in millimeter wave MIMO systems,” in _2020 IEEE 91st Veh. Technol. Conf._ , 2020, pp. 1–5.
  * [20] Y. Liu, S. Hong, C. Pan, Y. Wang, Y. Pan, and M. Chen, “Optimization of RIS Configurations for Multiple-RIS-Aided mmWave Positioning Systems based on CRLB Analysis,” 2021. [Online]. Available: https://arxiv.org/abs/2111.14023
  * [21] R. Wang, Z. Xing, and E. Liu, “Joint location and communication study for intelligent reflecting surface aided wireless communication system,” 2021. [Online]. Available: https://arxiv.org/abs/2103.01063
  * [22] Z. Feng, B. Wang, M. Luan, and F. Hu, “Power Optimization for Target Localization with Reconfigurable Intelligent Surfaces,” _Signal Process._ , vol. 189, no. 5, p. 108252, 2021.
  * [23] A. Elzanaty, A. Guerra, F. Guidi, and M.-S. Alouini, “Reconfigurable Intelligent Surfaces for Localization: Position and Orientation Error Bounds,” _IEEE Trans. Signal Process._ , vol. 69, pp. 5386–5402, 2021.
  * [24] H. Zhang, H. Zhang, B. Di, K. Bian, Z. Han, and L. Song, “Towards Ubiquitous Positioning by Leveraging Reconfigurable Intelligent Surface,” _IEEE Commun. Lett._ , vol. 25, no. 1, pp. 284–288, 2021.
  * [25] J. He, H. Wymeersch, T. Sanguanpuak, O. Silven, and M. Juntti, “Adaptive Beamforming Design for mmWave RIS-Aided Joint Localization and Communication,” in _2020 IEEE Wireless Commun. and Net. Conf. Works. (WCNCW)_ , 2020, pp. 1–6.
  * [26] A. Fascista, M. F. Keskin, A. Coluccia, H. Wymeersch, and G. Seco-Granados, “RIS-aided joint localization and synchronization with a single-antenna receiver: Beamforming design and low-complexity estimation,” _IEEE J. Sel. Topics Signal Process._ , vol. 16, no. 5, pp. 1141–1156, 2022.
  * [27] T. Wu, C. Pan, Y. Pan, H. Ren, M. Elkashlan, and C.-X. Wang, “Fingerprint Based mmWave Positioning System Aided by Reconfigurable Intelligent Surface,” 2022. [Online]. Available: https://arxiv.org/abs/2210.11998
  * [28] K. C. Ho, “Bias Reduction for an Explicit Solution of Source Localization Using TDOA,” _IEEE Trans. Signal Process._ , vol. 60, no. 5, pp. 2101–2114, 2012.
  * [29] Y. T. Chan and K. C. Ho, “A simple and efficient estimator for hyperbolic location,” _IEEE Trans. Signal Process._ , vol. 42, no. 8, pp. 1905–1915, 2002.
  * [30] C. Hu, L. Dai, S. Han, and X. Wang, “Two-Timescale Channel Estimation for Reconfigurable Intelligent Surface Aided Wireless Communications,” _IEEE Trans. Commun._ , vol. 69, no. 11, pp. 7736–7747, 2021.
  * [31] Y. Wang and K. C. Ho, “An Asymptotically Efficient Estimator in Closed-Form for 3D AOA Localization Using a Sensor Network,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 12, pp. 6524–6535, 2015.
  * [32] H. Wymeersch and B. Denis, “Beyond 5G Wireless Localization with Reconfigurable Intelligent Surfaces,” in _ICC 2020 - 2020 IEEE International Conf. Commun. (ICC)_ , 2020, pp. 1–6.
  * [33] R. A. Horn and C. R. Johnson, _Matrix Analysis_. Cambridge University Press, 1985.