Refining Inaccurate Transmitter and Receiver Positions Using Calibration Targets for Target Localization in Multi-Static Passive Radar

Abstract

Transmitter and receiver position errors have been known to significantly deteriorate target localization accuracy in a multi-static passive radar (MPR) system. This paper explores the use of calibration targets, whose positions are known to the MPR system, to counter the loss in target localization accuracy arising from transmitter/receiver position errors. This paper firstly evaluates the Cramér–Rao lower bound (CRLB) for bistatic range (BR)-based target localization with calibration targets, which analytically indicates the potential of calibration targets in enhancing localization accuracy. After that, this paper proposes a novel closed-form solution, which includes two steps: calibration step and localization step. Firstly, the calibration step is devoted to refine the inaccurate transmitter and receiver locations using the BR measurements from the calibration targets, and then in the calibration step, the target localization can be accurately achieved by using the refined transmitter/receiver positions and the BR measurements from the unknown target. Theoretical analysis and simulation results indicate that the proposed method can attain the CRLB at moderate measurement noise level, and exhibits the superiority of localization accuracy over existing algorithms.

1. Introduction

Passive radar technology, which allows operators to detect and localize potential targets using already existing transmitters such as commercial frequency modulation (FM) broadcast/digital audio broadcast (DAB)/terrestrial digital video broadcast (DVB-T) [1,2] and non-cooperative radar transmission [3,4], has been interesting to both civilian and military fields in the last few decades [5]. This sort of radar, compared to active radar technology, offers numerous advantages including lower cost, lower power usage and more covert surveillance capability, which suggests the possibility of employing passive radars on a wide range of concerned applications such as homeland security, costal surveillance and early warning system for vehicles detection, etc.

One of the remarkable characters of passive radar is the deployment of two receiving channels, with one for capturing the direct path signal from the transmitter and the other for collecting the potential target echoes [6]. By performing a cross-correlation operation between the direct path signal and the target echoes, the time delay (TD) could be measured, which holds information about the unknown target position. By multiplying with the signal propagation speed, the TD can be directly converted into the bistatic range (BR) [7]. Each BR measurement traces out an ellipsoid equation, with its foci located at the transmitter and the receiver positions. Theoretically, for multi-static passive radar (MPR), if enough BR measurements with respect to multiple transmitter-receiver pairs are available, the target position can be determined by solving the set of nonlinear ellipsoid equations. However, due to the high nonlinearity implied in the BR measurement equations, determining the target position from the BR measurements obtained at a single time instant is not a trivial topic.

In recent years, numerous algorithms have been developed to address this challenging topic, mainly including iterative methods [8,9] and closed-form solution methods [10,11,12,13,14]. Iterative methods [8,9] rely on an initial position guess close to the true solution but such a good guess may not always be available in reality. By contrast, closed-form solution methods [10,11,12,13,14] have always been more attractive to researchers due to their computational efficiency, independence on initial guess and absence of divergence problem. Illuminated by Ho and Xu’s two-step weighted least squares (2WLS) idea [15], these closed-form solution methods [10,11,12,13,14] generally follow the basic two-step framework below: in the first step, the non-linear BR equations are linearized into pseudo-linear ones by introducing some proper nuisance parameters and a coarse estimate of target position can be obtained from the pseudo-linear equation set via weighted least squares (WLS) minimization; then in the second step, the function relation between the target position and the introduced nuisance parameters is explored to refine the initial estimate.

Nevertheless, all the aforementioned methods are designed based on the assumption that the positions of the transmitters and receivers are exactly known, but such exact priori knowledge may not be available in reality. In fact, the positions of the transmitters and receivers are inevitably perturbed by errors to some extent, and these errors (also referred to as position uncertainties) are often non-negligible, especially when the antennas are mounted on moving platforms [16,17] or the exploited transmitters are highly non-cooperative (such as the hostile radar radiation whose position could usually only roughly determined by electronic reconnaissance technique [18,19]). On the other hand, Rui and Ho [20] quantitatively analyze the influence of the transmitter and receiver position error on the localization accuracy, indicating that the target localization accuracy can be very sensitive to the transmitter/receiver position error and a slight error in transmitter/receiver position could remarkably deteriorate the localization accuracy. More recently, some novel methods [21,22] that take the statistical distributions of transmitter/receiver position error into consideration are developed to reduce the target localization error, and they are shown to attain the Cramér–Rao lower bound (CRLB) under small measurement noise and transmitter/receiver position error assumption. Nevertheless, these methods [21,22] only present the solutions when the transmitter/receiver position errors exist, but cannot fundamentally compensate the localization performance loss arising from transmitter/receiver position error at the CRLB level.

The use of calibration sensors has been a common technique in wireless sensor network self-localization, where each sensor broadcasts signals and receives signals from other sensors so as to determine their positions collaboratively [23,24]. Syldatk [25,26] considers the calibration of ground sensor networks where an accurate calibration of sensor positions and orientations is required for target tracking. For source localization problem, Hasan first considered in [27] utilizing calibration sensors to improve the angle-based source localization performance. Ho [28], Yang [29] as well as Li [30] et al. further expanded and applied the calibration technique to range difference (RD)-based source localization problem, where additional RD measurements from the calibration sensors were incorporated to reduce the receivers’ position error and thus improve the source localization accuracy. The successful use of calibration techniques in these fields inspires us the possibility of employing calibration technique in the target localization for multi-static passive radar. When it comes to BR-based target localization in multi-static passive radar, using a ‘calibration target’ with known position may also be able to mitigate the target localization performance loss arising from the transmitter/receiver position error. In theory, any target appearing in the radar coverage area and meanwhile broadcasting its position can be taken as a calibration target. Typically, for example, to avoid potential accidents and collisions, the commercial aircrafts will report their positions and other information to the ground stations and other aircraft by using the automatic dependent surveillance broadcast (ADS-B) system [31]. Hence, the commercial aircraft broadcasting ADS-B signal can be regarded as an off-the-shelf calibration target. If no such off-the-shelf calibration targets are available in the radar coverage area, we can manually launch some strong scatterers with known positions as calibration targets. However, despite that, up to now there exists no publication in the open literature that addresses refining inaccurate transmitter and receiver positions using calibration targets for target localization in multi-static passive radar.

Motivated by these facts, in this paper, taking the transmitter and receiver position error into consideration, we explore using calibration targets to counter the loss in BR-based target localization accuracy arising from the transmitter/receiver position error. We begin our work by deriving the CRLB for BR-based target position estimation when the BR measurements from the calibration targets are available. The interpretation on the CRLB demonstrates that the use of calibration targets can significantly mitigate the influence of the transmitter/receiver position error and dramatically enhance the localization accuracy, at least in the sense of CRLB. We then proceed to develop a novel localization method to alleviate the transmitter and receiver position error and enhance the localization accuracy using calibration targets. It mainly includes two processing stages, referred to as calibration stage and localization stage respectively. In the calibration stage, the BR measurements from the calibration targets are exploited to refine the inaccurate transmitter and receiver positions; in the localization stage, the refined transmitter and receiver positions and BR measurements corresponding to the unknown target are exploited to determine the target position. Both processing stages are closed-form, which brings the proposed method computational efficiency and high robustness. Furthermore, the accuracy of the proposed solution is shown analytically to reach the CRLB when the transmitters’/receivers’/calibration targets’ position errors and the BR measurement noises are sufficiently small. Simulations will be conducted to verify the effectiveness and superiority of the proposed solution over existing methods.

Notations: There will be a lot of notations throughout this paper. Without exception, vectors (matrices) are denoted by bold lower (upper) case letters, respectively. Also, notations ${(·)}^{°}$, ${(·)}^{\mathrm{T}}$, $\Vert ·\Vert$, ${(·)}^{-1}$, $\mathrm{E}(·)$, $O_{p\times q}$, $I_{p\times p}$, $0_{p\times 1}$, $\mathrm{diag}(·)$ and $\mathrm{tr}(·)$, represent the true value of a noisy or an estimated variable, transpose operation, Euclidean norm, inverse of matrix, statistical expectation, a p-by-q zero matrix, an identity matrix of size p, a p-by-1 zero vector, diagonal matrix and the trace of a matrix, respectively.

The remainder of this paper is organized as follows. Section 2 is about the localization scenario in the presence of transmitter/receiver position error and calibration targets. In Section 3, the corresponding CRLB is evaluated, indicating the potential of calibration targets in improve localization accuracy. In Section 4, a closed-form solution is developed for target localization in the presence of transmitter/receiver position error and calibration targets, and theoretical performance analysis is also given. Section 5 describes the results of Monte Carlo simulations that compare the proposed solution with existing methods. Section 6 is the conclusion.

2. Problem Formulation

Address a typical multi-static passive radar localization scenario as presented in Figure 1, where M non-cooperative transmitters located at $s_{\mathrm{t},m}^{\mathrm{o}}={[x_{\mathrm{t},m}^{\mathrm{o}},y_{\mathrm{t},m}^{\mathrm{o}},z_{\mathrm{t},m}^{\mathrm{o}}]}^{\mathrm{T}}$ ($m=1,2,\dots ,M$) are employed to illuminate the surveillance area, and N receivers located at $s_{\mathrm{r},n}^{\mathrm{o}}={[x_{\mathrm{r},n}^{\mathrm{o}},y_{\mathrm{r},n}^{\mathrm{o}},z_{\mathrm{r},n}^{\mathrm{o}}]}^{\mathrm{T}}$ ($n=1,2,\dots ,N$) are deployed to determine a single target’s position denoted by $u^{\mathrm{o}}={[x^{\mathrm{o}},y^{\mathrm{o}},z^{\mathrm{o}}]}^{\mathrm{T}}$. In fact, the exact positions of the transmitters and receivers might not be known, and only the inaccurate measured versions, i.e., $s_{\mathrm{t},m}={[x_{\mathrm{t},m},y_{\mathrm{t},m},z_{\mathrm{t},m}]}^{\mathrm{T}}$ and $s_{\mathrm{r},n}={[x_{\mathrm{r},n},y_{\mathrm{r},n},z_{\mathrm{r},n}]}^{\mathrm{T}}$, are available for processing. Formulaically, we arrive at

$$s_{\mathrm{t},m}=s_{\mathrm{t},m}^{\mathrm{o}}+\mathbf{\Delta }{s}_{\mathrm{t},m}$$

(1)

$$s_{\mathrm{r},n}=s_{\mathrm{r},n}^{\mathrm{o}}+\mathbf{\Delta }{s}_{\mathrm{r},n}$$

(2)

where $\mathbf{\Delta }{s}_{\mathrm{t},m}$ and $\mathbf{\Delta }{s}_{\mathrm{r},n}$ are the position error of the mth transmitter and the nth receiver respectively and also referred to as position uncertainty. Stacking (1) and (2) with respect to all the transmitters and receivers, yields a 3(M + N)-by-1 transmitter and receiver position vector as

$$s=s^{\mathrm{o}}+\mathbf{\Delta }s$$

(3)

where $s={[s_{\mathrm{t}}^{\mathrm{T}},s_{\mathrm{r}}^{\mathrm{T}}]}^{\mathrm{T}}$ with $s_{\mathrm{t}}={[s_{\mathrm{t},1}^{\mathrm{T}},s_{\mathrm{t},2}^{\mathrm{T}},\dots ,s_{\mathrm{t},M}^{\mathrm{T}}]}^{\mathrm{T}}$ and $s_{\mathrm{r}}={[s_{\mathrm{r},1}^{\mathrm{T}},s_{\mathrm{r},2}^{\mathrm{T}},\dots ,s_{\mathrm{r},N}^{\mathrm{T}}]}^{\mathrm{T}}$ is the noisy transmitter and receiver position vector, $s^{\mathrm{o}}={[{(s_{\mathrm{t}}^{\mathrm{o}})}^{\mathrm{T}},{(s_{\mathrm{r}}^{\mathrm{o}})}^{\mathrm{T}}]}^{\mathrm{T}}$ with $s_{\mathrm{t}}^{\mathrm{o}}={[{(s_{\mathrm{t},1}^{\mathrm{o}})}^{\mathrm{T}},{(s_{\mathrm{t},2}^{\mathrm{o}})}^{\mathrm{T}},\dots ,{(s_{\mathrm{t},M}^{\mathrm{o}})}^{\mathrm{T}}]}^{\mathrm{T}}$ and $s_{\mathrm{r}}^{\mathrm{o}}={[{(s_{\mathrm{r},1}^{\mathrm{o}})}^{\mathrm{T}},{(s_{\mathrm{r},2}^{\mathrm{o}})}^{\mathrm{T}},\dots ,{(s_{\mathrm{r},N}^{\mathrm{o}})}^{\mathrm{T}}]}^{\mathrm{T}}$ is the true transmitter and receiver position vector, and $\mathbf{\Delta }s={[\mathbf{\Delta }{s}_{\mathrm{t}}^{\mathrm{T}},\mathbf{\Delta }{s}_{\mathrm{r}}^{\mathrm{T}}]}^{\mathrm{T}}$ with $\mathbf{\Delta }{s}_{\mathrm{t}}={[\mathbf{\Delta }{s}_{\mathrm{t},1}^{\mathrm{T}},\mathbf{\Delta }{s}_{\mathrm{t},2}^{\mathrm{T}},\dots ,\mathbf{\Delta }{s}_{\mathrm{t},M}^{\mathrm{T}}]}^{\mathrm{T}}$ and $\mathbf{\Delta }{s}_{\mathrm{r}}={[\mathbf{\Delta }{s}_{\mathrm{r},1}^{\mathrm{T}},\mathbf{\Delta }{s}_{\mathrm{r},2}^{\mathrm{T}},\dots ,\mathbf{\Delta }{s}_{\mathrm{r},N}^{\mathrm{T}}]}^{\mathrm{T}}$ is the transmitter and receiver position error vector that can be assumed zero-mean Gaussian with covariance $Q_s$ without loss of generality.

Using the above notations, the distance from the mth transmitter to the target of interest is equal to

$$R_{\mathrm{t},m}^{\mathrm{o}}=\Vert u^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}}\Vert$$

(4)

the distance from the target of interest to the nth receiver is equal to

$$R_{\mathrm{r},n}^{\mathrm{o}}=\Vert u^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$$

(5)

and the baseline distance with respect to the mth transmitter and nth receiver is

$$R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}}=\Vert s_{\mathrm{t},m}^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$$

(6)

According to this, the BR measurement with respect to the mth transmitter and nth receiver, i.e., the sum of the distances from the mth transmitter to the target and the target to the nth receiver, can be formulized as

$$\begin{array}{l}{r}_{m,n}=r_{m,n}^{\mathrm{o}}+\mathbf{\Delta }{r}_{m,n}\\ =R_{\mathrm{t},m}^{\mathrm{o}}+R_{\mathrm{r},n}^{\mathrm{o}}-R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}}+\mathbf{\Delta }{r}_{m,n}\\ =\Vert u^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}}\Vert +\Vert u^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert -\Vert s_{\mathrm{t},m}^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert +\mathbf{\Delta }{r}_{m,n}\end{array}$$

(7)

where $r_{m,n}^{\mathrm{o}}=R_{\mathrm{t},m}^{\mathrm{o}}+R_{\mathrm{r},n}^{\mathrm{o}}-R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}}$ represents the true BR with respect to the mth transmitter and nth receiver, $\mathbf{\Delta }{r}_{m,n}$ is the BR measurement noise. Herein, it is important to emphasize that the time delay measurement comes from the cross correlation operation between the target signal and the direct path reference signal [7], and its error characteristics are not affected by the transmitter/receiver position error since the transmitter/receiver position error is not involved in the estimation of time delay. Thus, the BR measurement noise is only related to the time delay measurement noise, and not related to the transmitter/receiver position error. Obviously, there will be MN BR measurements to be produced with respect to the M transmitters and N receivers, which can be recast into a MN-by-1 vector as

$$r=r^{\mathrm{o}}+\mathbf{\Delta }r$$

(8)

where $r={[r_1^{\mathrm{T}},r_2^{\mathrm{T}},\dots ,r_M^{\mathrm{T}}]}^{\mathrm{T}}$ with $r_m={[r_{m,1},r_{m,2},\dots ,r_{m,N}]}^{\mathrm{T}}$ is the BR measurement vector, $r^{\mathrm{o}}={[{(r_1^{\mathrm{o}})}^{\mathrm{T}},{(r_2^{\mathrm{o}})}^{\mathrm{T}},\dots ,{(r_M^{\mathrm{o}})}^{\mathrm{T}}]}^{\mathrm{T}}$ with $r_m^{\mathrm{o}}={[r_{m,1}^{\mathrm{o}},r_{m,2}^{\mathrm{o}},\dots ,r_{m,N}^{\mathrm{o}}]}^{\mathrm{T}}$ is the true BR vector, and $\mathbf{\Delta }r={[\mathbf{\Delta }{r}_1^{\mathrm{T}},\mathbf{\Delta }{r}_2^{\mathrm{T}},\dots ,\mathbf{\Delta }{r}_M^{\mathrm{T}}]}^{\mathrm{T}}$ with $\mathbf{\Delta }{r}_m={[\mathbf{\Delta }{r}_{m,1},\mathbf{\Delta }{r}_{m,2},\dots ,\mathbf{\Delta }{r}_{m,N}]}^{\mathrm{T}}$ is the BR measurement noise vector, which is usually assumed follow a Gaussian distribution with zero-mean and covariance $Q_r$.

As presented in Figure 1, to alleviate the transmitter/receiver position error and enhance localization accuracy, K calibration targets located at $c_k^{\mathrm{o}}={[x_{\mathrm{c},k}^{\mathrm{o}},y_{\mathrm{c},k}^{\mathrm{o}},z_{\mathrm{c},k}^{\mathrm{o}}]}^{\mathrm{T}}$ ($k=1,2,\dots ,K$) are employed, and the BRs among the calibration targets and the transmitter-receiver pairs are also measured. Similarly, the exact positions of the calibration targets are not known to us, and the nominal versions denoted by $c_k={[x_{\mathrm{c},k},y_{\mathrm{c},k},z_{\mathrm{c},k}]}^{\mathrm{T}}$ ($k=1,2,\dots ,K$) are given as

$$c_k=c_k^{\mathrm{o}}+\mathbf{\Delta }{c}_k$$

(9)

where $\mathbf{\Delta }{c}_k$ is the position error of the kth calibration target. Collecting (9) for all the K calibration targets forms a 3K-by-1 calibration target position vector as

$$c=c^{\mathrm{o}}+\mathbf{\Delta }c$$

(10)

where $c={[c_1^{\mathrm{T}},c_2^{\mathrm{T}},\dots ,c_K^{\mathrm{T}}]}^{\mathrm{T}}$ is the nominal calibration target position vector, $c^{\mathrm{o}}={[{(c_1^{\mathrm{o}})}^{\mathrm{T}},{(c_2^{\mathrm{o}})}^{\mathrm{T}},\dots ,{(c_K^{\mathrm{o}})}^{\mathrm{T}}]}^{\mathrm{T}}$ is the true calibration target position vector, and $\mathbf{\Delta }c={[\mathbf{\Delta }{c}_1^{\mathrm{T}},\mathbf{\Delta }{c}_2^{\mathrm{T}},\dots ,\mathbf{\Delta }{c}_K^{\mathrm{T}}]}^{\mathrm{T}}$ is the calibration target position error vector that is usually supposed to obey Gaussian distribution with zero-mean and covariance $Q_c$. Herein, it should be pointed out that, the positions of calibration targets are generally considered to be more precise compared with those of the transmitters and receivers, although they are also contaminated by errors.

Then, the distance from the mth transmitter to the kth calibration target is given by

$$R_{\mathrm{c},k,\mathrm{t},m}^{\mathrm{o}}=\Vert c_k^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}}\Vert$$

(11)

and the distance from the kth calibration target to the nth receiver is given by

$$R_{\mathrm{c},k,\mathrm{r},n}^{\mathrm{o}}=\Vert c_k^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$$

(12)

Based on this, the BR measurement corresponding to the kth calibration target, mth transmitter and nth receiver can be modeled as

$$\begin{array}{l}{r}_{\mathrm{c},k,m,n}=r_{\mathrm{c},k,m,n}^{\mathrm{o}}+\mathbf{\Delta }{r}_{\mathrm{c},k,m,n}\\ =R_{\mathrm{c},k,\mathrm{t},m}^{\mathrm{o}}+R_{\mathrm{c},k,\mathrm{r},n}^{\mathrm{o}}-R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}}+\mathbf{\Delta }{r}_{\mathrm{c},k,m,n}\\ =\Vert c_k^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}}\Vert +\Vert c_k^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert -\Vert s_{\mathrm{t},m}^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert +\mathbf{\Delta }{r}_{\mathrm{c},k,m,n}\end{array}$$

(13)

where $\mathbf{\Delta }{r}_{\mathrm{c},k,m,n}$ represents measurement noise in $r_{\mathrm{c},k,m,n}$, $r_{\mathrm{c},k,m,n}^{\mathrm{o}}=R_{\mathrm{c},k,\mathrm{t},m}^{\mathrm{o}}+R_{\mathrm{c},k,\mathrm{r},n}^{\mathrm{o}}-R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}}$ represents the true BR with respect to the kth calibration target, mth transmitter and nth receiver. Collecting (13) for the set of K calibration targets, M transmitters and N receivers, results in a KMN-by-1 vector as

$$r_{\mathrm{c}}=r_{\mathrm{c}}^{\mathrm{o}}+\mathbf{\Delta }{r}_{\mathrm{c}}$$

(14)

where $r_{\mathrm{c}}={[r_{\mathrm{c},1}^{\mathrm{T}},r_{\mathrm{c},2}^{\mathrm{T}},\dots ,r_{\mathrm{c},K}^{\mathrm{T}}]}^{\mathrm{T}}$ with $r_{\mathrm{c},k}={[r_{\mathrm{c},k,1}^{\mathrm{T}},r_{\mathrm{c},k,2}^{\mathrm{T}},\dots ,r_{\mathrm{c},k,M}^{\mathrm{T}}]}^{\mathrm{T}}$ and $r_{\mathrm{c},k,m}={[r_{\mathrm{c},k,m,1},r_{\mathrm{c},k,m,2},\dots ,r_{\mathrm{c},k,m,N}]}^{\mathrm{T}}$ denotes the BR measurement vector from the calibration targets, $r_{\mathrm{c}}^{\mathrm{o}}={[{(r_{\mathrm{c},1}^{\mathrm{o}})}^{\mathrm{T}},{(r_{\mathrm{c},2}^{\mathrm{o}})}^{\mathrm{T}},\dots ,{(r_{\mathrm{c},K}^{\mathrm{o}})}^{\mathrm{T}}]}^{\mathrm{T}}$ with $r_{\mathrm{c},k}^{\mathrm{o}}={[{(r_{\mathrm{c},k,1}^{\mathrm{o}})}^{\mathrm{T}},{(r_{\mathrm{c},k,2}^{\mathrm{o}})}^{\mathrm{T}},\dots ,{(r_{\mathrm{c},k,M}^{\mathrm{o}})}^{\mathrm{T}}]}^{\mathrm{T}}$ and $r_{\mathrm{c},k,m}^{\mathrm{o}}={[r_{\mathrm{c},k,m,1}^{\mathrm{o}},r_{\mathrm{c},k,m,2}^{\mathrm{o}},\dots ,r_{\mathrm{c},k,m,N}^{\mathrm{o}}]}^{\mathrm{T}}$ denotes the corresponding true value vector, $\mathbf{\Delta }{r}_{\mathrm{c}}={[\mathbf{\Delta }{r}_{\mathrm{c},1}^{\mathrm{T}},\mathbf{\Delta }{r}_{\mathrm{c},2}^{\mathrm{T}},\dots ,\mathbf{\Delta }{r}_{\mathrm{c},K}^{\mathrm{T}}]}^{\mathrm{T}}$ with $\mathbf{\Delta }{r}_{\mathrm{c},k}={[\mathbf{\Delta }{r}_{\mathrm{c},k,1}^{\mathrm{T}},\mathbf{\Delta }{r}_{\mathrm{c},k,2}^{\mathrm{T}},\dots ,\mathbf{\Delta }{r}_{\mathrm{c},k,M}^{\mathrm{T}}]}^{\mathrm{T}}$ and $\mathbf{\Delta }{r}_{\mathrm{c},k,m}={[\mathbf{\Delta }{r}_{\mathrm{c},k,m,1},\mathbf{\Delta }{r}_{\mathrm{c},k,m,2},\dots ,\mathbf{\Delta }{r}_{\mathrm{c},k,m,N}]}^{\mathrm{T}}$ denotes the corresponding error vector, which is presumed to be a Gaussian random vector with zero mean and covariance $Q_{r\mathrm{c}}$.

Now, the purpose of this work is to determine the target position from the noisy BR measurements and the inaccurate transmitter/receiver positions. In particular, the calibration targets with known position and the corresponding BR measurements are also available to reduce the transmitter/receiver position error and improve localization accuracy.

3. Evaluation of the CRLB with Calibration Targets

The CRLB does not address the specific estimators employed, but simply reflects minimum possible variance that an unbiased estimator can achieve with existing observations. In this section, in order to justify the necessity of refining the inaccurate transmitter and receiver positions using calibration targets, we shall first set up the CRLB for the target localization problem described above. Besides the BR measurement noise, the position errors of transmitter, receivers and calibration targets are also included. From the localization scenario presented in Section 2, the deterministic but unknown parameters for the CRLB evaluation, collected into a 3(M + N + K + 1)-by-1 vector $\mathsf{\phi }={[{(u^{\mathrm{o}})}^{\mathrm{T}},{(s^{\mathrm{o}})}^{\mathrm{T}},{(c^{\mathrm{o}})}^{\mathrm{T}}]}^{\mathrm{T}}$, include the target position vector $u^{\mathrm{o}}$, the transmitter and receiver position vector $s^{\mathrm{o}}$, and the calibration target position vector $c^{\mathrm{o}}$; the observations, collected into a (MN + KMN + 3M + 3N + 3K)-by-1 vector $z={[r^{\mathrm{T}},r_{\mathrm{c}}^{\mathrm{T}},s^{\mathrm{T}},c^{\mathrm{T}}]}^{\mathrm{T}}$, include the BR measurement vector r from the unknown target, the BR measurement vector $r_{\mathrm{c}}$ from the calibration targets, the inaccurate measured transmitter and receiver position vector s, and the nominal calibration target position vector c, which are Gaussian distributed and independent with one another. Based on this, the joint probability density function (pdf) of the observations parameterized by the unknown parameter vector is readily shown to be

$$\begin{array}{l}p(z|\mathsf{\phi })=p(r|u^{\mathrm{o}},s^{\mathrm{o}})\cdot p(r_{\mathrm{c}}|s^{\mathrm{o}},c^{\mathrm{o}})\cdot p(s|s^{\mathrm{o}})\cdot p(c|c^{\mathrm{o}})\\ =\kappa \cdot \exp [-\frac{1}{2}{(r-r^{\mathrm{o}})}^{\mathrm{T}}{Q}_r^{-1}(r-r^{\mathrm{o}})-\frac{1}{2}{(r_{\mathrm{c}}-r_{\mathrm{c}}^{\mathrm{o}})}^{\mathrm{T}}{Q}_{r\mathrm{c}}^{-1}(r_{\mathrm{c}}-r_{\mathrm{c}}^{\mathrm{o}})\\ -\frac{1}{2}{(s-s^{\mathrm{o}})}^{\mathrm{T}}{Q}_s^{-1}(s-s^{\mathrm{o}})-\frac{1}{2}{(c-c^{\mathrm{o}})}^{\mathrm{T}}{Q}_c^{-1}(c-c^{\mathrm{o}})]\end{array}$$

(15)

where $\kappa$ is a constant with respect to the unknown parameters. By taking the logarithm of (15), partial derivatives with respect to the unknown parameters twice, and then expectation, the Fisher information matrix (FIM) can be calculated as

$$\begin{array}{l}FIM(\mathsf{\phi })=\mathrm{E}\left[\frac{\partial \ln p(z|\mathsf{\phi })}{\partial \mathsf{\phi }}\left(\frac{\partial \ln p(z|\mathsf{\phi })}{\partial \mathsf{\phi }}\right)\right]\\ =\left[\begin{array}{ccc}X& Y& O_{3\times 3}\\ Y^{\mathrm{T}}& Z& R^{\mathrm{T}}\\ O_{3\times 3}& R& P\end{array}\right]\end{array}$$

(16)

where the blocks X, Y, Z, R and P are respectively given by

$$X={\left(\frac{\partial r^{\mathrm{o}}}{\partial u^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_r^{-1}\left(\frac{\partial r^{\mathrm{o}}}{\partial u^{\mathrm{o}}}\right)$$

(17)

$$Y={\left(\frac{\partial r^{\mathrm{o}}}{\partial u^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_r^{-1}\left(\frac{\partial r^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)$$

(18)

$$Z=Q_s^{-1}+{\left(\frac{\partial r^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_r^{-1}\left(\frac{\partial r^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)+{\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_{r\mathrm{c}}^{-1}\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)$$

(19)

$$R={\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial c^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_{r\mathrm{c}}^{-1}\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)$$

(20)

$$P=Q_c^{-1}+{\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial c^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_{r\mathrm{c}}^{-1}\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial c^{\mathrm{o}}}\right)$$

(21)

Denote $i_{m,n}=(m-1)N+n$, and $i_{k,m,n}=(k-1)MN+(m-1)N+n$. From the formulations of (7) and (13), the elements of the partial derivatives $\partial r^{\mathrm{o}}/\partial u^{\mathrm{o}}$, $\partial r^{\mathrm{o}}/\partial s^{\mathrm{o}}$, $\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial c^{\mathrm{o}}$ and $\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial s^{\mathrm{o}}$ in (17)–(21), can be determined as

$$\frac{\partial r^{\mathrm{o}}}{\partial u^{\mathrm{o}}}(i_{m,n},1:3)=\frac{{(u^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{t},m}^{\mathrm{o}}}+\frac{{(u^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{r},n}^{\mathrm{o}}}$$

(22)

$$\frac{\partial r^{\mathrm{o}}}{\partial s^{\mathrm{o}}}=\left[\begin{array}{cc}\frac{\partial r^{\mathrm{o}}}{\partial s_{\mathrm{t}}^{\mathrm{o}}}& \frac{\partial r^{\mathrm{o}}}{\partial s_{\mathrm{r}}^{\mathrm{o}}}\end{array}\right]$$

(23)

$$\frac{\partial r^{\mathrm{o}}}{\partial s_{\mathrm{t}}^{\mathrm{o}}}(i_{m,n},3m-2:3m)=\frac{{(s_{\mathrm{t},m}^{\mathrm{o}}-u^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{t},m}^{\mathrm{o}}}-\frac{{(s_{\mathrm{t},m}^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}}}$$

(24)

$$\frac{\partial r^{\mathrm{o}}}{\partial s_{\mathrm{r}}^{\mathrm{o}}}(i_{m,n},3n-2:3n)=\frac{{(s_{\mathrm{r},n}^{\mathrm{o}}-u^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{r},n}^{\mathrm{o}}}-\frac{{(s_{\mathrm{r},n}^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}}}$$

(25)

$$\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial c^{\mathrm{o}}}(i_{k,m,n},3k-2:3k)=\frac{{(c_k^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{c},k,\mathrm{t},m}^{\mathrm{o}}}+\frac{{(c_k^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{c},k,\mathrm{r},n}^{\mathrm{o}}}$$

(26)

$$\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s^{\mathrm{o}}}=\left[\begin{array}{cc}\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s_{\mathrm{t}}^{\mathrm{o}}}& \frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s_{\mathrm{r}}^{\mathrm{o}}}\end{array}\right]$$

(27)

$$\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s_{\mathrm{t}}^{\mathrm{o}}}(i_{k,m,n},3m-2:3m)=\frac{{(s_{\mathrm{t},m}^{\mathrm{o}}-c_k^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{c},k,\mathrm{t},m}^{\mathrm{o}}}-\frac{{(s_{\mathrm{t},m}^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}}}$$

(28)

$$\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s_{\mathrm{r}}^{\mathrm{o}}}(i_{k,m,n},3n-2:3n)=\frac{{(s_{\mathrm{r},n}^{\mathrm{o}}-c_k^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{c},k,\mathrm{r},n}^{\mathrm{o}}}-\frac{{(s_{\mathrm{r},n}^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}})}^{\mathrm{T}}}{R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}}}$$

(29)

for $k=1,2,\dots ,K$, $m=1,2,\dots ,M$ and $n=1,2,\dots ,N$, and zeros elsewhere.

By definition, the CRLB of $\mathsf{\phi }$, denoted by $CRL{B}_{\mathrm{c}}(\mathsf{\phi })$, is given as $\mathrm{FIM}{(\mathsf{\phi })}^{-1}$, where only the upper left 3-by-3 block is for the target position $u^{\mathrm{o}}$. Invoking the partitioned matrix inversion formula as well as the matrix inversion lemma [32] twice on (16), leads to the CRLB of $u^{\mathrm{o}}$ as

$$CRL{B}_{\mathrm{c}}(u^{\mathrm{o}})=X^{-1}+X^{-1}Y{(Z-Y^{\mathrm{T}}{X}^{-1}Y-R^{\mathrm{T}}{P}^{-1}R)}^{-1}{Y}^{\mathrm{T}}{X}^{-1}$$

(30)

For comparison purposes, the CRLB of $u^{\mathrm{o}}$ with transmitter/receiver position error but without calibration derived in [22], denoted by $CRL{B}_{\mathrm{s}}(u^{\mathrm{o}})$, is also given below

$$CRL{B}_{\mathrm{s}}(u^{\mathrm{o}})=X^{-1}+X^{-1}Y{(\stackrel{⌢}{Z}-Y^{\mathrm{T}}{X}^{-1}Y)}^{-1}{Y}^{\mathrm{T}}{X}^{-1}$$

(31)

where $\stackrel{⌢}{Z}=Q_s^{-1}+{(\partial r^{\mathrm{o}}/\partial s^{\mathrm{o}})}^{\mathrm{T}}{Q}_r^{-1}(\partial r^{\mathrm{o}}/\partial s^{\mathrm{o}})$. For the sake of comparison, we proceed to construct an equivalent form of $CRL{B}_{\mathrm{c}}(u^{\mathrm{o}})$ by denoting $\stackrel{⌣}{Z}$ as $\stackrel{⌣}{Z}=Z-R^{\mathrm{T}}{P}^{-1}R$. After invoking the matrix inversion lemma [32] to ${(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial s^{\mathrm{o}})}^{\mathrm{T}}{Q}_{r\mathrm{c}}^{-1}(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial s^{\mathrm{o}})-R^{\mathrm{T}}{P}^{-1}R$ and some algebraic manipulations, we further represent $\stackrel{⌣}{Z}$ as

$$\stackrel{⌣}{Z}=Q_s^{-1}+{\left(\frac{\partial r^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_r^{-1}\left(\frac{\partial r^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)+{\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)}^{\mathrm{T}}{\left(Q_{r\mathrm{c}}+\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial c^{\mathrm{o}}}\right)Q_c{\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial c^{\mathrm{o}}}\right)}^{\mathrm{T}}\right)}^{-1}\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)$$

(32)

Using (32), we obtain an equivalent expression of $CRL{B}_{\mathrm{c}}(u^{\mathrm{o}})$ as

$$CRL{B}_{\mathrm{c}}(u^{\mathrm{o}})=X^{-1}+X^{-1}Y{(\stackrel{⌣}{Z}-Y^{\mathrm{T}}{X}^{-1}Y)}^{-1}{Y}^{\mathrm{T}}{X}^{-1}$$

(33)

Through the comparison of (31) and (33), it is readily to observe that the two CRLBs are identical in structure, except that $\stackrel{⌣}{Z}$ is substituted by $\stackrel{⌢}{Z}$. More specifically, the use of calibration targets introduces an additional component into the bracketed matrix expression to be inverted as

$$\begin{array}{l}\tilde{Z}=\stackrel{⌣}{Z}-\stackrel{⌢}{Z}\\ ={(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial s^{\mathrm{o}})}^{\mathrm{T}}{Q}_{r\mathrm{c}}^{-1}(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial s^{\mathrm{o}})-R^{\mathrm{T}}{P}^{-1}R\\ ={(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial s^{\mathrm{o}})}^{\mathrm{T}}{\left(Q_{r\mathrm{c}}+(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial c^{\mathrm{o}})Q_c{(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial c^{\mathrm{o}})}^{\mathrm{T}}\right)}^{-1}(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial s^{\mathrm{o}})\end{array}$$

(34)

Using (34), we can rewrite ${(\stackrel{⌣}{Z}-Y^{\mathrm{T}}{X}^{-1}Y)}^{-1}$ in (33) as ${((\stackrel{⌢}{Z}-Y^{\mathrm{T}}{X}^{-1}Y)+\tilde{Z})}^{-1}$. Invoking the matrix inversion lemma [32] to the term ${((\stackrel{⌢}{Z}-Y^{\mathrm{T}}{X}^{-1}Y)+\tilde{Z})}^{-1}$ in (33), we obtain after some algebraic manipulations,

$$CRL{B}_{\mathrm{c}}(u^{\mathrm{o}})-CRL{B}_{\mathrm{s}}(u^{\mathrm{o}})=X^{-1}Y\Gamma Y^{\mathrm{T}}{X}^{-1}$$

(35)

where

$$\Gamma ={{\rm H}}^{-1}\mathsf{\Upsilon }{(I+{\mathsf{\Upsilon }}^{\mathrm{T}}{{\rm H}}^{-1}\mathsf{\Upsilon })}^{-1}{\mathsf{\Upsilon }}^{\mathrm{T}}{{\rm H}}^{-1}$$

(36)

$${\rm H}=(\tilde{Z}-Y^{\mathrm{T}}{X}^{-1}Y)$$

(37)

$$\mathsf{\Upsilon }={\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)}^{\mathrm{T}}{L}_{r\mathrm{c}}$$

(38)

and $L_{r\mathrm{c}}$ is the Cholesky decomposition of ${(Q_{r\mathrm{c}}+(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial c^{\mathrm{o}})Q_c{(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial c^{\mathrm{o}})}^{\mathrm{T}})}^{-1}$, i.e., $L_{r\mathrm{c}}{L}_{r\mathrm{c}}^{\mathrm{T}}={(Q_{r\mathrm{c}}+(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial c^{\mathrm{o}})Q_c{(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial c^{\mathrm{o}})}^{\mathrm{T}})}^{-1}$. In form, the right side of (35) is just the performance enhancement because of the use of calibration targets. It is positive semi-definite (PSD) since it has a symmetric structure and ${\mathsf{\Upsilon }}^{\mathrm{T}}$ is not full column rank. Even if the nominal positions of calibration targets and the corresponding BR measurements are very noisy, (35) can still remain PSD. In theory, only in the edge case when ${(Q_{r\mathrm{c}}+(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial c^{\mathrm{o}})Q_c^{-1}{(\partial r_{\mathrm{c}}^{\mathrm{o}}/\partial c^{\mathrm{o}})}^{\mathrm{T}})}^{-1}$ tends to zero and then $L_{r\mathrm{c}}\to O$ and $\mathsf{\Upsilon }\to O$, the performance enhancement in (34) would tend to zero. However, this edge case hardly exists in reality. Thus, mathematically, we can arrive at

$$CRL{B}_{\mathrm{s}}(u^{\mathrm{o}})\ge CRL{B}_{\mathrm{c}}(u^{\mathrm{o}})$$

(39)

The matrix inequality $A\ge B$ means that $A-B$ is PSD. It can be further deduced from (39) that $\mathrm{tr}(CRL{B}_{\mathrm{s}}(u^{\mathrm{o}}))\ge \mathrm{tr}(CRL{B}_{\mathrm{c}}(u^{\mathrm{o}}))$. The trace of $CRL{B}_{\mathrm{c}}(u^{\mathrm{o}})$ and $CRL{B}_{\mathrm{s}}(u^{\mathrm{o}})$ respectively represents minimum possible variance of target position estimation with and without using calibration targets. Therefore, we can conclude that using calibration targets brings potential enhancement to the target localization accuracy, at least at the CRLB level.

Example 1. To substantiate the evaluation on the CRLB presented above, a numerical example using a typical multi-static passive radar localization scenario was conducted, as presented in Figure 2. There are M = 3 transmitters, N = 4 receivers and K = 3 calibration targets in the scenario, and their true positions are listed in

Table 1. Positions (in meters) of the transmitters, receivers and calibration targets.

TX	${\mathit{x}}_{\mathbf{t},\mathit{m}}^{\mathbf{o}}$	${\mathit{y}}_{\mathbf{t},\mathit{m}}^{\mathbf{o}}$	${\mathit{z}}_{\mathbf{t},\mathit{m}}^{\mathbf{o}}$	RX	${\mathit{x}}_{\mathbf{r},\mathit{n}}^{\mathbf{o}}$	${\mathit{y}}_{\mathbf{r},\mathit{n}}^{\mathbf{o}}$	${\mathit{z}}_{\mathbf{r},\mathit{n}}^{\mathbf{o}}$	Calibration Targets	${\mathit{x}}_{\mathbf{c},\mathit{k}}^{\mathbf{o}}$	${\mathit{y}}_{\mathbf{c},\mathit{k}}^{\mathbf{o}}$	${\mathit{z}}_{\mathbf{c},\mathit{k}}^{\mathbf{o}}$
1	20,000	0	100	1	2000	2000	0	1	10,000	10,000	2500
2	15,000	5000	1000	2	2000	−2000	500	2	15,000	30,000	3000
3	15,000	−5000	2000	3	5000	5000	1000	3	20,000	−10,000	3500
--	--	--	--	4	5000	−5000	1500	--	--	--	--

. The noise covariance matrix of the BR measurements from the unknown target are given by $Q_r={\sigma }_{\mathrm{r}}^2{V}_r$, where ${\sigma }_{\mathrm{r}}$ reflects BR measurement noise level and $V_r$ is set to 1 in the diagonal elements and 0.5 elsewhere. The covariance matrix of the transmitter/receiver position error is given as $Q_s={\sigma }_{\mathrm{s}}^2{V}_s$ where ${\sigma }_s$ reflects the transmitter/receiver position error level and $V_{\mathrm{s}}=\mathrm{diag}(5I_{3M\times 3M},I_{3N\times 3N})$. The covariance matrix of calibration target position error is $Q_c={\sigma }_{\mathrm{c}}^2{V}_c$ where ${\sigma }_{\mathrm{c}}$ reflects the calibration target position error level and $V_c=I_{3K\times 3K}$, and the covariance matrix of the corresponding calibration BR measurement noise is $Q_{r\mathrm{c}}={\sigma }_{\mathrm{rc}}^2{V}_{r\mathrm{c}}$ where ${\sigma }_{\mathrm{rc}}={\sigma }_{\mathrm{r}}$ reflects calibration BR measurement noise level and $V_{r\mathrm{c}}$ is set to 1 in the diagonal elements and 0.5 elsewhere. The target of interest is located at position $u^{\mathrm{o}}={[50000,15000,5000]}^{\mathrm{T}}\mathrm{m}$. The effect of the calibration targets on the target localization accuracy, in the sense of CRLB, is presented in Figure 3.

Figure 3a compares the CRLB curves with and without using calibration targets when the BR measurement noise level ${\sigma }_{\mathrm{r}}$ is varied from ${10}^{-2}\mathrm{m}$ to ${10}^3\mathrm{m}$ while the transmitter/receiver position error level and calibration target position error level are fixed at ${\sigma }_{\mathrm{s}}=20\mathrm{m}$ and ${\sigma }_{\mathrm{c}}=10\mathrm{m}$ respectively. It can be observed from Figure 3a that the CRLB with calibration targets is generally below the one without, this coincides with the analytical conclusion given in (39). However, in the edge case where the BR measurement noise is very large, two CRLBs would tend to be the same. This is because in this case, the BR measurement noise dominates and effect of transmitter/receiver position error on the localization accuracy is relatively small. The CRLB curves versus the transmitter/receiver position error level ${\sigma }_{\mathrm{s}}$ are plotted in Figure 3b where the BR measurement noise level and calibration target position error level are fixed at ${\sigma }_{\mathrm{r}}=10\mathrm{m}$ and ${\sigma }_{\mathrm{c}}=10\mathrm{m}$ respectively. A similar trend, i.e., two CRLBs would tend to be the same, appears in Figure 3b, when the transmitter/receiver position error is sufficiently small. A reasonable explanation is that, in this case, the transmitter/receiver positions are known very accurately and their influence on the localization accuracy can be ignored compared to the BR measurement noise. The CRLB comparison versus calibration target position error level ${\sigma }_{\mathrm{c}}$ is provided in Figure 3c where ${\sigma }_{\mathrm{r}}=10\mathrm{m}$ and ${\sigma }_{\mathrm{s}}=20\mathrm{m}$. Interestingly, the trend of CRLB curves implies that, even when the calibration target position error is extremely large, the CRLB with utilization of calibration targets are still remarkably below the one without. This justifies again the analysis under (35), and similar results have also been presented in previous studies [23,24,25,26,27,28,29,30] on source localization and sensor network localization issues. Generally, from Figure 3, the use of calibration targets brings a significant improvement in the localization accuracy in the normal case, at least at the CRLB level.

4. Proposed Localization Method

The evaluation of the CRLB in Section 3 has demonstrated the potential of calibration targets in improving localization accuracy. In what follows, we will proceed to develop a novel closed-form solution for the aforementioned practical localization scenario where the positions of transmitters and receivers are inaccurate but calibration targets are used to refine the transmitter/receiver position and enhance the localization accuracy. After that, a theoretical analysis will be performed to show that the proposed solution achieves the CRLB when satisfying some mild conditions.

4.1. Algorithm Development

The proposed solution mainly includes two processing stages, referred to as calibration stage and localization stage, respectively. The calibration stage is devoted to refining the inaccurate transmitter and receiver positions, and then the localization stage is devoted to determining the target position on the basis of the refined transmitter and receiver positions.

4.1.1. Calibration Stage

To make use of the BR measurements from the calibration targets, the calibration stage begins by reorganizing (13) as

$$r_{\mathrm{c},k,m,n}-\Vert c_k^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}}\Vert -\Vert c_k^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert +\Vert s_{\mathrm{t},m}^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert =\mathbf{\Delta }{r}_{\mathrm{c},k,m,n}$$

(40)

Since only the erroneous versions of $c_k^{\mathrm{o}}$, $s_{\mathrm{t},m}^{\mathrm{o}}$ and $s_{\mathrm{r},n}^{\mathrm{o}}$ are available, we put $c_k^{\mathrm{o}}=c_k-\mathbf{\Delta }{c}_k$, $s_{\mathrm{t},m}^{\mathrm{o}}=s_{\mathrm{t},m}-\mathbf{\Delta }{s}_{\mathrm{t},m}$ and $s_{\mathrm{r},n}^{\mathrm{o}}=s_{\mathrm{r},n}-\mathbf{\Delta }{s}_{\mathrm{r},n}$ into (40), and then expand it around erroneous values $c_k$, $s_{\mathrm{t},m}$ and $s_{\mathrm{r},n}$ to the linear error terms as

$$\begin{array}{l}{r}_{\mathrm{c},k,m,n}-\Vert c_k-s_{\mathrm{t},m}\Vert -\Vert c_k-s_{\mathrm{r},n}\Vert +\Vert s_{\mathrm{t},m}-s_{\mathrm{r},n}\Vert -({\mathsf{\rho }}_{\mathrm{c},k,\mathrm{t},m}^{\mathrm{T}}+{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{T}})\mathbf{\Delta }{s}_{\mathrm{t},m}-({\mathsf{\rho }}_{\mathrm{c},k,\mathrm{r},n}^{\mathrm{T}}-{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{T}})\mathbf{\Delta }{s}_{\mathrm{r},n}\\ =-{({\mathsf{\rho }}_{\mathrm{c},k,\mathrm{t},m}+{\mathsf{\rho }}_{\mathrm{c},k,\mathrm{r},n})}^{\mathrm{T}}\mathbf{\Delta }{c}_k+\mathbf{\Delta }{r}_{\mathrm{c},k,m,n}\end{array}$$

(41)

where

$${\mathsf{\rho }}_{\mathrm{c},k,\mathrm{t},m}=\frac{c_k-s_{\mathrm{t},m}}{\Vert c_k-s_{\mathrm{t},m}\Vert }$$

(42)

$${\mathsf{\rho }}_{\mathrm{c},k,\mathrm{r},n}=\frac{c_k-s_{\mathrm{r},n}}{\Vert c_k-s_{\mathrm{r},n}\Vert }$$

(43)

$${\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},n}=\frac{s_{\mathrm{t},m}-s_{\mathrm{r},n}}{\Vert s_{\mathrm{t},m}-s_{\mathrm{r},n}\Vert }$$

(44)

Stacking (41) for all the k, m and n, we can formulate them in matrix form as

$$h_0-G_0\mathbf{\Delta }s=\mathbf{\Delta }{h}_0$$

(45)

The elements of $h_0$, $G_0$ and $\mathbf{\Delta }{h}_0$ are given by

$$h_0(i_{k,m,n},1)=r_{\mathrm{c},k,m,n}-\Vert c_k-s_{\mathrm{t},m}\Vert -\Vert c_k-s_{\mathrm{r},n}\Vert +\Vert s_{\mathrm{t},m}-s_{\mathrm{r},n}\Vert$$

(46)

$$G_0=\left[\begin{array}{cc}{G}_{0,\mathrm{t}}& G_{0,\mathrm{r}}\end{array}\right],G_{0,\mathrm{t}}(i_{k,m,n},3m-2:3m)={\mathsf{\rho }}_{\mathrm{c},k,\mathrm{t},m}^{\mathrm{T}}+{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{T}},G_{0,\mathrm{r}}(i_{k,m,n},3n-2:3n)={\mathsf{\rho }}_{\mathrm{c},k,\mathrm{r},n}^{\mathrm{T}}-{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{T}}$$

(47)

$$G_{\mathrm{c}}(i_{k,m,n},3k-2:3k)=-{({\mathsf{\rho }}_{\mathrm{c},k,\mathrm{t},m}+{\mathsf{\rho }}_{\mathrm{c},k,\mathrm{r},n})}^{\mathrm{T}}$$

(48)

$$\mathbf{\Delta }{h}_0(i_{k,m,n},1)=-{({\mathsf{\rho }}_{\mathrm{c},k,\mathrm{t},m}+{\mathsf{\rho }}_{\mathrm{c},k,\mathrm{r},n})}^{\mathrm{T}}\mathbf{\Delta }{c}_k+\mathbf{\Delta }{r}_{\mathrm{c},k,m,n}$$

(49)

for $i_{k,m,n}=(k-1)MN+(m-1)N+n$, $k=0,1,\dots ,K-1$, $m=0,1,\dots ,M-1$, $n=0,1,\dots ,N-1$, and zeros elsewhere. Furthermore, the error vector $\mathbf{\Delta }{h}_0$ can be recast using a compact representation as follows

$$\mathbf{\Delta }{h}_0=G_{\mathrm{c}}\mathbf{\Delta }c+\mathbf{\Delta }{r}_{\mathrm{c}}$$

(50)

from which we have the mean $\mathrm{E}(\mathbf{\Delta }{h}_0)=0_{KMN\times 1}$ and the covariance $\mathrm{cov}(\mathbf{\Delta }{h}_0)=G_{\mathrm{c}}{Q}_c{G}_{\mathrm{c}}^{\mathrm{T}}+Q_{r\mathrm{c}}$. In (45), $\mathbf{\Delta }s$ represents the difference between the true and the nominal transmitter/receiver positions.

In order to refine the transmitter and receiver positions, $\mathbf{\Delta }s$ shall be estimated as accurately as possible. Recall that $\mathbf{\Delta }s$ is a Gaussian distributed random vector with mean $\mathrm{E}(\mathbf{\Delta }s)=0_{3(M+N)\times 1}$ and covariance matrix $\mathrm{cov}(\mathbf{\Delta }s)=Q_s$, and it is independent of the error vector $\mathbf{\Delta }{h}_0$. Thus according to the Bayesian Gauss–Markov theorem [33], the linear minimum mean square error (LMMSE) estimate of $\mathbf{\Delta }s$ can be obtained from (45) as

$$\begin{array}{l}\mathbf{\Delta }\widehat{s}=\mathrm{E}(\mathbf{\Delta }s)+{\left(\mathrm{cov}{(\mathbf{\Delta }s)}^{-1}+G_0^{\mathrm{T}}\mathrm{cov}{(\mathbf{\Delta }{h}_0)}^{-1}{G}_{\mathrm{s}}\right)}^{-1}{G}_0^{\mathrm{T}}\mathrm{cov}{(\mathbf{\Delta }{h}_0)}^{-1}\left(h_0-G_0\mathrm{E}(\mathbf{\Delta }s)\right)\\ ={\left(Q_s^{-1}+G_0^{\mathrm{T}}{(G_{\mathrm{c}}{Q}_c{G}_{\mathrm{c}}^{\mathrm{T}}+Q_{r\mathrm{c}})}^{-1}{G}_0\right)}^{-1}{G}_0^{\mathrm{T}}{(G_{\mathrm{c}}{Q}_c{G}_{\mathrm{c}}^{\mathrm{T}}+Q_{r\mathrm{c}})}^{-1}{h}_0\end{array}$$

(51)

Under the assumption that the noise in $G_{\mathrm{c}}$ and $G_0$ is sufficiently small to be ignored, the covariance matrix of $\mathbf{\Delta }\widehat{s}$ can be given as

$$\mathrm{cov}(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})={\left(Q_s^{-1}+G_0^{\mathrm{T}}{(G_{\mathrm{c}}{Q}_c{G}_{\mathrm{c}}^{\mathrm{T}}+Q_{r\mathrm{c}})}^{-1}{G}_0\right)}^{-1}$$

(52)

Using the estimate of transmitter and receiver position error in (51), we can refine the transmitter and receiver positions as

$$\widehat{s}=s-\mathbf{\Delta }\widehat{s}$$

(53)

Utilizing the fact $s=s^{\mathrm{o}}+\mathbf{\Delta }s$, we can rewrite $\widehat{s}$ in (53) as $\widehat{s}=s^{\mathrm{o}}+\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s}$. Hence, the refined estimate of transmitter/receiver positions $\widehat{s}$ has a covariance matrix identical with (52). Forming the inverse of $\mathrm{cov}(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})$ and then comparing it to $Q_s^{-1}$ results in $\mathrm{cov}{(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})}^{-1}-Q_s^{-1}=G_0^{\mathrm{T}}{(G_{\mathrm{c}}{Q}_c{G}_{\mathrm{c}}^{\mathrm{T}}+Q_{r\mathrm{c}})}^{-1}{G}_0$. It is natural to deduce that $\mathrm{cov}{(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})}^{-1}\ge Q_s^{-1}$ is PSD since $G_0^{\mathrm{T}}{(G_{\mathrm{c}}{Q}_c{G}_{\mathrm{c}}^{\mathrm{T}}+Q_{r\mathrm{c}})}^{-1}{G}_0$ has a symmetric structure and $G_{\mathrm{c}}$ is not full column rank. According to the PSD matrix property [34], $\mathrm{cov}{(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})}^{-1}\ge Q_s^{-1}$ is equivalent to $Q_s\ge \mathrm{cov}(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})$. That is to say, the refined positions of transmitters and receivers performs leastwise as well as, if not better than, the original ones, in terms of target localization accuracy.

4.1.2. Localization Stage

The localization stage starts by linearizing the BR equations from the unknown target. Firstly, reorganize (7) as

$$(r_{m,n}+R_{\mathrm{t},m,\mathrm{r},n}^{\mathrm{o}})-R_{\mathrm{t},m}^{\mathrm{o}}=R_{\mathrm{r},n}^{\mathrm{o}}+\mathbf{\Delta }{r}_{m,n}$$

(54)

Since we have obtained the refined estimate of $s_{\mathrm{t},m}^{\mathrm{o}}$ and $s_{\mathrm{r},n}^{\mathrm{o}}$ from calibration stage, we plug $s_{\mathrm{t},m}^{\mathrm{o}}={\widehat{s}}_{\mathrm{t},m}-(\mathbf{\Delta }{s}_{\mathrm{t},m}-\mathbf{\Delta }{\widehat{s}}_{\mathrm{t},m})$ and $s_{\mathrm{r},n}^{\mathrm{o}}={\widehat{s}}_{\mathrm{r},n}-(\mathbf{\Delta }{s}_{\mathrm{r},n}-\mathbf{\Delta }{\widehat{s}}_{\mathrm{r},n})$ into (54) and ignoring the second and higher error terms as

$$\begin{array}{l}2{({\widehat{s}}_{\mathrm{t},m}-{\widehat{s}}_{\mathrm{r},n})}^{\mathrm{T}}{u}^{\mathrm{o}}+2(r_{m,n}+R_{\mathrm{t},m,\mathrm{r},n})R_{\mathrm{t},m}^{\mathrm{o}}={(r_{m,n}+R_{\mathrm{t},m,\mathrm{r},n})}^2+{\widehat{s}}_{\mathrm{t},m}^{\mathrm{T}}{\widehat{s}}_{\mathrm{t},m}-{\widehat{s}}_{\mathrm{r},n}^{\mathrm{T}}{\widehat{s}}_{\mathrm{r},n}-2R_{\mathrm{r},n}^{\mathrm{o}}\mathbf{\Delta }{r}_{m,n}\\ +2{(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{t},m}-R_{\mathrm{r},n}^{\mathrm{o}}{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},n})}^{\mathrm{T}}(\mathbf{\Delta }{s}_{\mathrm{t},m}-\mathbf{\Delta }{\widehat{s}}_{\mathrm{t},m})-2{(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{r},n}-R_{\mathrm{r},n}^{\mathrm{o}}{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},n})}^{\mathrm{T}}(\mathbf{\Delta }{s}_{\mathrm{r},n}-\mathbf{\Delta }{\widehat{s}}_{\mathrm{r},n})\end{array}$$

(55)

By forming an auxiliary vector as ${\mathsf{\theta }}^{\mathrm{o}}={[{(u^{\mathrm{o}})}^{\mathrm{T}},R_{\mathrm{t},1}^{\mathrm{o}},R_{\mathrm{t},2}^{\mathrm{o}},\dots ,R_{\mathrm{t},M}^{\mathrm{o}}]}^{\mathrm{T}}$, we can collect (55) for all the m and n into a matrix form as

$$G_1{\mathsf{\theta }}^{\mathrm{o}}=h_1+\mathbf{\Delta }{h}_1$$

(56)

where

$$G_1=\left[\begin{array}{cc}{G}_{1,\mathrm{s}}& G_{1,\mathrm{r}}\end{array}\right]$$

(57)

$$G_{1,\mathrm{s}}=2\left[\begin{array}{c}{\widehat{s}}_1\\ {\widehat{s}}_2\\ ⋮\\ {\widehat{s}}_M\end{array}\right],s_m=\left[\begin{array}{c}{({\widehat{s}}_{\mathrm{t},m}-{\widehat{s}}_{\mathrm{r},1})}^{\mathrm{T}}\\ {({\widehat{s}}_{\mathrm{t},m}-{\widehat{s}}_{\mathrm{r},2})}^{\mathrm{T}}\\ ⋮\\ {({\widehat{s}}_{\mathrm{t},m}-{\widehat{s}}_{\mathrm{r},N})}^{\mathrm{T}}\end{array}\right],$$

(58)

$$G_{1,\mathrm{r}}=2\left[\begin{array}{cccc}{r}_1& 0_{N\times 1}& \cdots & 0_{N\times 1}\\ 0_{N\times 1}& r_2& \cdots & 0_{N\times 1}\\ ⋮& ⋮& \ddots & ⋮\\ 0_{N\times 1}& 0_{N\times 1}& \cdots & r_M\end{array}\right],r_m=\left[\begin{array}{c}{r}_{m,1}\\ r_{m,2}\\ ⋮\\ r_{m,N}\end{array}\right]$$

(59)

$$h_1=\left[\begin{array}{c}{h}_{1,1}\\ h_{1,2}\\ ⋮\\ h_{1,M}\end{array}\right],h_{1,m}=\left[\begin{array}{c}{(r_{m,1}+R_{\mathrm{t},m,\mathrm{r},1})}^2+{\widehat{s}}_{\mathrm{t},m}^{\mathrm{T}}{\widehat{s}}_{\mathrm{t},m}-{\widehat{s}}_{\mathrm{r},1}^{\mathrm{T}}{\widehat{s}}_{\mathrm{r},1}\\ {(r_{m,2}+R_{\mathrm{t},m,\mathrm{r},2})}^2+{\widehat{s}}_{\mathrm{t},m}^{\mathrm{T}}{\widehat{s}}_{\mathrm{t},m}-{\widehat{s}}_{\mathrm{r},2}^{\mathrm{T}}{\widehat{s}}_{\mathrm{r},2}\\ ⋮\\ {(r_{m,N}+R_{\mathrm{t},m,\mathrm{r},N})}^2+{\widehat{s}}_{\mathrm{t},m}^{\mathrm{T}}{\widehat{s}}_{\mathrm{t},m}-{\widehat{s}}_{\mathrm{r},N}^{\mathrm{T}}{\widehat{s}}_{\mathrm{r},N}\end{array}\right]$$

(60)

and the error vector $\mathbf{\Delta }{h}_1$ is related to the target position as

$$\mathbf{\Delta }{h}_1=B_1\mathbf{\Delta }r+D_1(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})$$

(61)

where

$$B_1=\left[\begin{array}{cccc}{B}_{1,1}& O_{N\times N}& \cdots & O_{N\times N}\\ O_{N\times N}& B_{1,2}& \cdots & O_{N\times N}\\ ⋮& ⋮& \ddots & ⋮\\ O_{N\times N}& O_{N\times N}& \cdots & B_{1,M}\end{array}\right]$$

(62)

$$D_1=\left[\begin{array}{ccccc}{D}_{1,\mathrm{t},1}& O_{N\times 3}& \cdots & O_{N\times 3}& D_{1,\mathrm{r},1}\\ O_{N\times 3}& D_{1,\mathrm{t},2}& \cdots & O_{N\times 3}& D_{1,\mathrm{r},2}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ O_{N\times 3}& O_{N\times 3}& \cdots & D_{1,\mathrm{t},M}& D_{1,\mathrm{r},M}\end{array}\right]$$

(63)

with

$$B_{1,m}=-2\mathrm{diag}\left(R_{\mathrm{r},1}^{\mathrm{o}},R_{\mathrm{r},2}^{\mathrm{o}},\dots ,R_{\mathrm{r},N}^{\mathrm{o}}\right)$$

(64)

$$D_{1,\mathrm{t},m}=2\left[\begin{array}{c}{(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{t},m}-R_{\mathrm{r},1}^{\mathrm{o}}{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},1})}^{\mathrm{T}}\\ {(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{t},m}-R_{\mathrm{r},2}^{\mathrm{o}}{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},2})}^{\mathrm{T}}\\ ⋮\\ {(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{t},m}-R_{\mathrm{r},N}^{\mathrm{o}}{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},N})}^{\mathrm{T}}\end{array}\right]$$

(65)

$$D_{1,\mathrm{r},m}=-2\left[\begin{array}{cccc}{(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{r},1}-R_{\mathrm{r},1}^{\mathrm{o}}{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},1})}^{\mathrm{T}}& 0_{3\times 1}^{\mathrm{T}}& \cdots & 0_{3\times 1}^{\mathrm{T}}\\ 0_{3\times 1}^{\mathrm{T}}& {(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{r},2}-R_{\mathrm{r},2}^{\mathrm{o}}{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},2})}^{\mathrm{T}}& \cdots & 0_{3\times 1}^{\mathrm{T}}\\ ⋮& ⋮& \ddots & ⋮\\ 0_{3\times 1}^{\mathrm{T}}& 0_{3\times 1}^{\mathrm{T}}& \cdots & {(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{r},N}-R_{\mathrm{r},N}^{\mathrm{o}}{\mathsf{\rho }}_{\mathrm{t},m,\mathrm{r},N})}^{\mathrm{T}}\end{array}\right]$$

(66)

From the set of linear equations in (56), the WLS estimate of ${\mathsf{\theta }}^{\mathrm{o}}$, denoted by $\mathsf{\theta }$, which minimizes $\mathbf{\Delta }{h}_1^{\mathrm{T}}{W}_1\mathbf{\Delta }{h}_1$ can be produced as

$$\mathsf{\theta }={(G_1^{\mathrm{T}}{W}_1{G}_1)}^{-1}{G}_1^{\mathrm{T}}{W}_1{h}_1$$

(67)

where $W_1$ represents the weighting matrix and it can be computed by

$$\begin{array}{l}{W}_1={\left[\mathrm{E}(\mathbf{\Delta }{h}_1\mathbf{\Delta }{h}_1^{\mathrm{T}})\right]}^{-1}\\ ={\left[B_1{Q}_{\alpha }{B}_1^{\mathrm{T}}+D_1\mathrm{cov}(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})D_1^{\mathrm{T}}\right]}^{-1}\end{array}$$

(68)

However, to compute $W_1$, the unknown target position has to be acquired in advance. To resolve this contradiction, we preliminarily let $W_1=I_{MN\times MN}$ and use (67) to compute a least squares estimate of ${\mathsf{\theta }}^{\mathrm{o}}$, and then use the estimated ${\mathsf{\theta }}^{\mathrm{o}}$ to update $W_1$ for another repetition.

Based on the WLS theorem, it can be deduced that the estimate $\mathsf{\theta }$ is approximately unbiased and the corresponding covariance matrix can be obtained, given sufficiently small BR measurement noise and transmitter/receiver position error, as

$$\mathrm{cov}(\mathsf{\theta })={(G_1^{\mathrm{T}}{W}_1{G}_1)}^{-1}$$

(69)

Next, the functional relation between the target position $u^{\mathrm{o}}$ and the introduced nuisance parameters $R_{\mathrm{t},1}^{\mathrm{o}},R_{\mathrm{t},2}^{\mathrm{o}},\dots ,R_{\mathrm{t},M}^{\mathrm{o}}$, is explored to compute the final estimate of target position. To this end, reorganize the functional relation in (4) as

$$2{(s_{\mathrm{t},m}^{\mathrm{o}})}^{\mathrm{T}}{u}^{\mathrm{o}}={(u^{\mathrm{o}})}^{\mathrm{T}}{u}^{\mathrm{o}}+{(s_{\mathrm{t},m}^{\mathrm{o}})}^{\mathrm{T}}{s}_{\mathrm{t},m}^{\mathrm{o}}-{(R_{\mathrm{t},m}^{\mathrm{o}})}^2$$

(70)

Denoting the estimation error of $\mathsf{\theta }$ by $\mathbf{\Delta }\mathsf{\theta }$, mathematically we arrive at

$$u^{\mathrm{o}}=\mathsf{\theta }(1:3)-\mathbf{\Delta }\mathsf{\theta }(1:3)$$

(71)

$$R_{\mathrm{t},m}^{\mathrm{o}}=\mathsf{\theta }(3+m)-\mathbf{\Delta }\mathsf{\theta }(3+m)$$

(72)

Putting (71), (72) into the right side of (70) and $s_{\mathrm{t},m}^{\mathrm{o}}={\widehat{s}}_{\mathrm{t},m}-(\mathbf{\Delta }{s}_{\mathrm{t},m}-\mathbf{\Delta }{\widehat{s}}_{\mathrm{t},m})$ into the both sides, we have after ignoring second-order error terms,

$$\begin{array}{l}2{\widehat{s}}_{\mathrm{t},m}^{\mathrm{T}}{u}^{\mathrm{o}}=\mathsf{\theta }{(1:3)}^{\mathrm{T}}\mathsf{\theta }(1:3)-\mathsf{\theta }{(m+3)}^2+{\widehat{s}}_{\mathrm{t},m}^{\mathrm{T}}{\widehat{s}}_{\mathrm{t},m}\\ -2\mathsf{\theta }(1:3)\mathbf{\Delta }\mathsf{\theta }(1:3)+2\mathsf{\theta }(m+3)\mathbf{\Delta }\mathsf{\theta }(m+3)+2{(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{t},m})}^{\mathrm{T}}(\mathbf{\Delta }{s}_{\mathrm{t},m}-\mathbf{\Delta }{\widehat{s}}_{\mathrm{t},m})\end{array}$$

(73)

The final estimate of target position should satisfy (73) and meanwhile retain as close as possible to the estimated values of target position in $\mathsf{\theta }$. In line with this principle, one has the following set of equations

$$G_2{u}^{\mathrm{o}}=h_2+\mathbf{\Delta }{h}_2$$

(74)

where

$$G_2=\left[\begin{array}{c}{I}_{3\times 3}\\ 2s_{t,1}^{\mathrm{T}}\\ 2s_{t,2}^{\mathrm{T}}\\ ⋮\\ 2s_{t,M}^{\mathrm{T}}\end{array}\right]$$

(75)

$$h_2=\left[\begin{array}{c}\mathsf{\theta }(1:3)\\ \mathsf{\theta }{(1:3)}^{\mathrm{T}}\mathsf{\theta }(1:3)-\mathsf{\theta }{(1+3)}^2+{\widehat{s}}_{\mathrm{t},1}^{\mathrm{T}}{\widehat{s}}_{\mathrm{t},1}\\ \mathsf{\theta }{(1:3)}^{\mathrm{T}}\mathsf{\theta }(1:3)-\mathsf{\theta }{(2+3)}^2+{\widehat{s}}_{\mathrm{t},2}^{\mathrm{T}}{\widehat{s}}_{\mathrm{t},2}\\ ⋮\\ \mathsf{\theta }{(1:3)}^{\mathrm{T}}\mathsf{\theta }(1:3)-\mathsf{\theta }{(M+3)}^2+{\widehat{s}}_{\mathrm{t},M}^{\mathrm{T}}{\widehat{s}}_{\mathrm{t},M}\end{array}\right]$$

(76)

$$\mathbf{\Delta }{h}_2=B_2\mathbf{\Delta }\mathsf{\theta }+D_2(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})$$

(77)

$$B_2=\left[\begin{array}{ccccc}-I_{3\times 3}& 0_{3\times 1}& 0_{3\times 1}& \cdots & 0_{3\times 1}\\ -2\mathsf{\theta }(1:3)& 2\mathsf{\theta }(1+3)& 0& \cdots & 0\\ -2\mathsf{\theta }(1:3)& 0& 2\mathsf{\theta }(2+3)& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -2\mathsf{\theta }(1:3)& 0& 0& \cdots & 2\mathsf{\theta }(M+3)\end{array}\right]$$

(78)

$$D_2=\left[\begin{array}{cc}{O}_{3\times 3M}& O_{3\times 3N}\\ 2\mathrm{diag}\left\{{(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{t},1})}^{\mathrm{T}},{(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{t},2})}^{\mathrm{T}},\dots ,{(u^{\mathrm{o}}-{\widehat{s}}_{\mathrm{t},M})}^{\mathrm{T}}\right\}& O_{3M\times 3N}\end{array}\right]$$

(79)

Invoking the WLS theorem again, one has the solution of target position, denoted by u, from (74) as

$$u={(G_2^{\mathrm{T}}{W}_2{G}_2)}^{-1}{G}_2^{\mathrm{T}}{W}_2{h}_2$$

(80)

where $W_2$ is the weighting matrix and it is determined by

$$\begin{array}{l}{W}_2={\left[E(\mathbf{\Delta }{h}_2\mathbf{\Delta }{h}_2^{\mathrm{T}})\right]}^{-1}\\ =[B_2\mathrm{cov}(\mathsf{\theta })B_2^{\mathrm{T}}+D_2\mathrm{cov}(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})D_2^{\mathrm{T}}+B_2{(G_1^{\mathrm{T}}{W}_1{G}_1)}^{-1}{G}_1^{\mathrm{T}}{W}_1{D}_1\mathrm{cov}(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})D_2^{\mathrm{T}}\\ {+D_2\mathrm{cov}(\mathbf{\Delta }s-\mathbf{\Delta }\widehat{s})D_1^{\mathrm{T}}{W}_1{G}_1{(G_1^{\mathrm{T}}{W}_1{G}_1)}^{-1}{B}_2^{\mathrm{T}}]}^{-1}\end{array}$$

(81)

But as presented in (81), the unknown target position is required in the computation of $W_2$. Herein, to circumvent this dilemma, we preliminarily exploit the target position estimate contained in $\mathsf{\theta }$ to form $W_2$ and use (80) to estimate target position. After that we can utilize the estimated target position to update $W_2$ for another repetition.

From the WLS theorem, the covariance matrix of u can be approximated, given sufficiently small BR measurement noise and transmitter/receiver position error, as

$$\mathrm{cov}(u)={(G_2^{\mathrm{T}}{W}_2{G}_2)}^{-1}$$

(82)

4.2. Performance Analysis

As mentioned above, the CRLB traces out a lower bound for minimum possible variance that an unbiased estimator can achieve. Next, we will analyze the efficiency of the proposed solution by comparing its covariance matrix with the benchmark, i.e., CRLB. For derivation simplicity, we would compare their inverse, rather than directly compare the two separately. The CRLB has been presented in (33). By invoking the matrix inversion lemma [32] to (33) and using the definitions of X and Y, we have after mathematical simplifications,

$$CRL{B}_{\mathrm{c}}{(u^{\mathrm{o}})}^{-1}={\left(\frac{\partial r^{\mathrm{o}}}{\partial u^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_r^{-1}\left(\frac{\partial r^{\mathrm{o}}}{\partial u^{\mathrm{o}}}\right)-{\left(\frac{\partial r^{\mathrm{o}}}{\partial u^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_r^{-1}\left(\frac{\partial r^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right){\stackrel{⌣}{Z}}^{-1}{\left(\frac{\partial r^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right)}^{\mathrm{T}}{Q}_r^{-1}\left(\frac{\partial r^{\mathrm{o}}}{\partial u^{\mathrm{o}}}\right)$$

(83)

where the expression of $\stackrel{⌣}{Z}$ has been given in (32).

On the other hand, using (82), (81), (69), (68) and (52) successively, we can reformulate the inverse of $\mathrm{cov}(u)$ as

$$\mathrm{cov}{(u)}^{-1}=G_3^{\mathrm{T}}{Q}_r^{-1}{G}_3-G_3^{\mathrm{T}}{Q}_r^{-1}{G}_4{\overline{Z}}^{-1}{G}_4^{\mathrm{T}}{Q}_r^{-1}{G}_3$$

(84)

where $G_3=B_1^{-1}{G}_1{B}_2^{-1}{G}_2$, $G_4=B_1^{-1}{D}_1$, and $\overline{Z}=Q_s^{-1}+G_4^{\mathrm{T}}{Q}_r^{-1}{G}_4+G_0^{\mathrm{T}}{(G_c{Q}_c{G}_c^{\mathrm{T}}+Q_{r\mathrm{c}})}^{-1}{G}_0$.

Comparing (83) with (84), we observe that $CRL{B}_{\mathrm{c}}{(u^{\mathrm{o}})}^{-1}$ and $\mathrm{cov}{(u)}^{-1}$ are identical in structure. Next, we proceed to prove their equivalency under the following conditions:

(C1) $\Vert \mathbf{\Delta }{s}_{\mathrm{t},m}\Vert \ll \Vert c_k^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}}\Vert$, $\Vert \mathbf{\Delta }{s}_{\mathrm{t},m}\Vert \ll \Vert s_{\mathrm{t},m}^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$, $\Vert \mathbf{\Delta }{s}_{\mathrm{r},n}\Vert \ll \Vert c_k^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$, $\Vert \mathbf{\Delta }{s}_{\mathrm{r},n}\Vert \ll \Vert s_{\mathrm{t},m}^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$, and $\Vert \mathbf{\Delta }{c}_k\Vert \ll \Vert c_k^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}}\Vert$, $\Vert \mathbf{\Delta }{c}_k\Vert \ll \Vert c_k^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$, for $k=1,2,\dots ,K$, $m=1,2,\dots ,M$ and $n=1,2,\dots ,N$;

(C2) $\Vert \mathbf{\Delta }{r}_{m,n}\Vert \ll \Vert u^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}}\Vert$, $\Vert \mathbf{\Delta }{r}_{m,n}\Vert \ll \Vert u^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$, $\Vert \mathbf{\Delta }{r}_{m,n}\Vert \ll \Vert s_{\mathrm{t},m}^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$, and $\Vert \mathbf{\Delta }{s}_{\mathrm{t},m}-\mathbf{\Delta }{\widehat{s}}_{\mathrm{t},m}\Vert \ll$ $\Vert u^{\mathrm{o}}-s_{\mathrm{t},m}^{\mathrm{o}}\Vert$, $\Vert \mathbf{\Delta }{s}_{\mathrm{r},n}-\mathbf{\Delta }{\widehat{s}}_{\mathrm{r},n}\Vert \ll \Vert u^{\mathrm{o}}-s_{\mathrm{r},n}^{\mathrm{o}}\Vert$ for $m=1,2,\dots ,M$ and $n=1,2,\dots ,N$;

The condition C1 implies the transmitter/receiver position error and the calibration target position error are negligibly small compared with the range between the calibration target and the transmitter/receiver. The condition C2 implies the BR measurement noise and the error in the refined transmitter/receiver position are negligibly small compared to the range between the calibration target and the transmitter/receiver. Using the conditions C1 and C2, we obtain, after some involved algebraic manipulations, that

$$G_3=\frac{\partial r^{\mathrm{o}}}{\partial u^{\mathrm{o}}},G_4=-\frac{\partial r^{\mathrm{o}}}{\partial s^{\mathrm{o}}},G_0=-\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial s^{\mathrm{o}}}\right),G_c=-\left(\frac{\partial r_{\mathrm{c}}^{\mathrm{o}}}{\partial c^{\mathrm{o}}}\right)$$

(85)

By this point, we can draw the conclusion that

$$\mathrm{cov}{(u)}^{-1}=CRL{B}_{\mathrm{c}}{(u^{\mathrm{o}})}^{-1}$$

(86)

That is, the proposed solution accomplishes the CRLB accuracy if the two conditions C1 and C2 are satisfied. In reality, localization scenarios, which satisfy the conditions C1 and C2, are not rare. These two conditions can be satisfied if the unknown target and the calibration targets are far from the transmitters and receivers, if not these conditions can still be satisfied if the BR measurement noise and the transmitter/receiver/calibration target position errors are sufficiently small.

5. Simulation Results

In this section, the efficiency and superiority of the proposed solution will be corroborated through Monte Carlo simulations. Amiri’s method presented in [14], which does not consider the transmitter and receiver position error and Zhao’s method proposed in [22], which considers the statistical distributions of transmitter/receiver position error but does not use any calibration targets, are chosen as references for comparison. The exact positions of transmitters/receivers/calibration targets are the same as those in Table 1. Localization accuracy is quantitatively evaluated using root mean squares error (RMSE), which comes from 1000 independent Monte Carlo runs. In each run, the zero-mean Gaussian random errors with covariance matrices $Q_r={\sigma }_{\mathrm{r}}^2{V}_r$, $Q_{r\mathrm{c}}={\sigma }_{\mathrm{rc}}^2{V}_{r\mathrm{c}}$, $Q_s={\sigma }_{\mathrm{s}}^2{V}_s$ and $Q_c={\sigma }_{\mathrm{c}}^2{V}_c$ are added to the BRs from unknown target, BRs from calibration targets, actual transmitter and receiver positions, and actual calibration target positions, respectively, in order to simulate a real localization scenario. The setting of $V_r$, $V_{r\mathrm{c}}$, $V_s$ and $V_c$ are also the same as that in Example 1.

First of all, in order to intuitively show the difference between target localization with and without the use of calibration targets, we plot in Figure 4 the estimated target positions from each Monte Carlo run, which forms a scatter plot for target position estimation. For comparison, the scatterplots of Amiri’s method and Zhao’s method are also plotted. The transmitter/receiver position error level is set to ${\sigma }_{\mathrm{s}}=20\mathrm{m}$, the noise level of BR measurements from unknown target and calibration targets is set as ${\sigma }_{\mathrm{r}}={\sigma }_{\mathrm{rc}}=10\mathrm{m}$, and the calibration target position error level is set to be ${\sigma }_{\mathrm{c}}=10\mathrm{m}$. The true position of the unknown target is $u^{\mathrm{o}}={[50000,15000,5000]}^{\mathrm{T}}\mathrm{m}$, which is marked with red pentagram in Figure 4 for comparison. By comparing the scatterplots of the methods, we find that with the use of calibration targets, the scattered dots of target position estimation are more closely around the target’s true position, which intuitively illustrates the performance gain from the use of calibration targets. Without the use of calibration targets, considering the statistical distributions of transmitter/receiver position error can also reduce dispersion of estimated target position dots to some extent, but compared to using calibration targets, this degree of reduction in dispersion is not sufficiently impressive.

Now, in order to quantitatively evaluate the localization accuracy of the methods, we calculate the RMSE of the proposed solution under different error or noise conditions, and compare it with Amiri’s method, Zhao’s method, as well as the CRLB. As mentioned in Section 4.2, the localization accuracy of the proposed solution is related to the distance between the target and MPR system. Hence, in order to achieve a more comprehensive insight on the performance of the proposed solution, we consider two cases, i.e., the near-field case where the target is close to the MPR system, and the far-field case where the target is far away from the MPR system. The exact positions of transmitters/receivers/calibration targets remain the same as before. We first address the far-field target, whose position is set to $u^{\mathrm{o}}={[120000,120000,12000]}^{\mathrm{T}}\mathrm{m}$. The results are presented in Figure 5.

Figure 5a plots the RMSE curves of the methods versus the BR measurement noise level. It shows that the localization RMSE of the proposed solution matches the CRLB very well and is about an order of magnitude lower than that of Amiri’s method and Zhao’s method at a low-to-moderate BR measurement noise level. Although it deviates from the CRLB when the BR measurement noise level is large, it is still much smaller than that of other two methods. The deviation from the CRLB, known as the thresholding phenomenon, is due to the ignored second order error terms in the design of the solution, which is invalid for large error levels. Owing to considering the statistical distributions of the transmitter/receiver position error, the RMSEs produced by Zhao’s method is generally lower than that by Amiri’s method. But compared with the use of the calibration targets in the proposed solution, the localization accuracy improvement brought by the consideration of transmitter/receiver position error in Zhao’s method is not so significant. Figure 5b gives the RMSE curves of the methods versus the transmitter/receiver position error level. It can be seen that, the superiority of the proposed solution in localization accuracy is mainly reflected at moderate to high transmitter/receiver position error level. When the transmitter/receiver position error is small, the localization accuracy of the proposed solution and the other two methods is comparable. This again agrees very well with the theoretical performance in Section 3. Figure 5c compares the RMSEs from the methods with respect to different calibration target position error levels. As is illustrated in Figure 5c, the proposed solution always offers a remarkable advantage over the other two methods at different calibration target position error level, even when the calibration target position error is extremely large. This is in agreement with the previous simulation results for the CRLB in Section 3.

Next, the same set of simulations was repeated for a near-field target, whose position is set to be $u^{\mathrm{o}}={[12000,1200,1200]}^{\mathrm{T}}\mathrm{m}$. The results are provided in Figure 6, from which we observe that the proposed solution still performs much better than the other methods. However, comparing with the corresponding results in Figure 5, we find the localization accuracy for near-field target is generally better than a far-field target, given the same noise and error levels. One reason may be that, when the target is close to the MPR system, the transmitters/receivers are far apart relative to the distance between the target and the MPR system. Thus, the localization geometry would become more regular and the corresponding geometric dilution of precision (GDOP) value would be smaller compared to the far-field case. However, on the other hand, comparing the thresholding values in Figure 5 and Figure 6 indicates that the RMSE curves for the near-field target deviate from the CRLB at smaller values than those for the far-field target. This phenomenon is consistent with the analysis under (83) that the equivalency between the estimate variance and the CRLB is more affected by the BR measurement noises when the target is close to the MPR system.

At an intuitive level, the more calibration targets are used, the better the localization accuracy is. In what follows, we will quantitatively analyze the effect of number of calibration targets on the localization accuracy by varying the number of calibration targets from 1 to 10. The positions of the transmitters and receivers remain the same as before. The positions of calibration targets and unknown target are chosen randomly from the 50 km × 40 km × 5 km volume as presented in Figure 2. The simulation results are depicted in Figure 7.

Figure 7 shows the RMSE, as well as the CRLB, versus the number of calibration targets. As expected, when the number of calibration targets is small, the localization accuracy improves significantly as the number of calibration targets increases. However, it is seen that there is no obvious dependence on the number of calibration targets as soon as the number of calibration targets is larger than 3. This indicates that when the number of calibration targets reaches 3, the use of more calibration targets would only increase the computational expense and not remarkably enhance the localization accuracy. Therefore, in the absence of any other consideration, it is reasonable to set the number of calibration targets as 3.

6. Conclusions

This paper explores the use of calibration targets with known positions to refine the inaccurate transmitter/receiver positions and thus enhance target localization accuracy in MPR systems. We start our research by evaluating target localization CRLB in the presence of calibration targets, which justifies the potential of calibration targets in enhancing localization accuracy. Then, in order to fulfill this potential, a novel closed-form solution was designed for target localization using BR measurements from the unknown target as well as the calibration targets. The proposed solution was shown both analytically and numerically to attain the CRLB under some mild conditions, and verified to outperform existing methods in terms of localization accuracy. Furthermore, from the view of engineering practice, if the employed calibration targets are off-the-shelf, such as the commercial aircrafts broadcasting an ADS-B signal, the use of calibration targets would bring little added cost or complexity to the MPR system, but could bring a significant enhancement to target localization accuracy.