TDOA positioning in the presence of atmospheric refraction and observer uncertainty influence analysis and correction

1. Introduction

In the actual battlefield environment, most of the moving radiation sources are flying at equal altitude and parallel to the Earth's surface [1]. By resorting to this constraint, the position of moving sources can be estimated accurately and simultaneously by using TDOA measurements from only three stations. Due to the reduced propagation speed of radio signal in the atmosphere and the influence of uneven atmospheric refraction, atmospheric refraction error will appear in the measured time difference, resulting in position error.

In single-station active radar applications, distance, azimuth and elevation angles are utilized, and the atmospheric refraction errors can be calculated directly according to the explicit atmospheric refraction model. Due to the different target positioning systems, the atmospheric refraction error correction method of single station active radar is not suitable for multi station TDOA positioning system. In order to improve the convergence speed and accuracy, a piecewise refraction model is innovatively proposed. Based on the stability of the error difference of atmospheric refraction TDOA, an iterative correction algorithm of successive approximation is proposed to correct the atmospheric refraction error of TDOA location, and analyze the impact of atmospheric refraction on the location accuracy under different random noises. To a certain extent, it can meet the requirements of modern military high-precision positioning.

2. Problem formulation

As shown in Fig. 1, the moving source radiates signals at the unknown position of $\mathbf{u}=(x,y,z{)}^T$, with its altitude given by $h_o$. A total of three observers are introduced to collect the radio signals emitted from the source. The position of the ith static observer is ${\mathbf{s}}_i=(x_i,y_i,z_i{)}^T$, $i=$1, 2, 3, with the altitude of the ith observer given by $h_i$. The apparent elevational angle, which is essentially the angel between the tangent direction and horizontal direction of the signal propagation path at the receiving station $i$, is denoted by ${\theta }_i$, while the unknown true elevational angle ${\alpha }_i$. The apparent elevational angle at the radiation source is ${\theta }_{Ti}$, while the true elevational angle is ${\alpha }_{Ti}$, with the elevational angle error given by ${\epsilon }_{Ti}$. The center of the earth is denoted by $C$, the radius of the earth by $r_o$, the linear propagation distance from the radiation source to the receiving station by $d_{oi}$, and the actual propagation path (curved distance) by $d_{gi}$. $n_T$ is the refractive index at the radiation source, $c$ is the speed of light, $f_o$ is the carrier frequency of the radiation source, and the wavelength is $\lambda =c/f_o$.

As the atmospheric variations in the vertical direction are of several orders larger in magnitude than those in the horizontal direction, the atmosphere refraction could thus be considered as the spherical stratification, the refractive index can be simplified to a quantity that varies only with elevation $h$, namely $n=n\left(h\right)$.The atmospheric refractive index $n$ is close to 1, the ground value is about 1.00026-1.00046. For instance, the value of $n$ at the elevation of 9 km is about 1.000105 [2]. Another physical quantity of refraction index is introduced, namely:

where $N\left(h\right)$ is the refraction index at altitude $h$.

The model error of piecewise model is proved to be relatively small, thus this model is adopted in this paper, which is given as:

where $h_s$ is the altitude of the ground, $N_o$ is the refraction index of the ground, and $\mathrm{\Delta }{N}_1$ is the gradient of the refraction index when $h_s\le h\le h_s+1\mathrm{}\mathrm{k}\mathrm{m}$. $N_1$ is the refraction index at the altitude $\left(h_s+1\right)\mathrm{}$km, $c_{\alpha 1}$ is the exponential attenuation rate from 1 km to 9 km over the ground, and $N_9$ is the refraction index at the altitude of 9 km. $c_{\alpha 9}$ is the exponential decay rate of altitude between 9 km-60 km. The curved propagation time from the source to station $i$ is $t_{gi}=\frac{1}{c}{\int }_{h_i}^{h_o}g\left(h\right)dh$, the apparent distance from the radiation source to the receiving station $i$ is:

where, $g\left(h\right)=\frac{n\left(h{)}^2\right(r_o+h)}{\sqrt{n\left(h{)}^2\right(r_o+h{)}^2-{n_i}^2(r_o+h_i{)}^2\mathrm{c}\mathrm{o}{\mathrm{s}}^2{\theta }_i}}$, $n_i$ is the index of refraction at altitude $h_i$. ${\theta }_i$ is the apparent elevational angle of the receiving station $i$, the Euclidean distance from the radiation source to the receiving station $i$ is:

with the propagation time given by:

The true elevational angle of the radiation source in the station-center coordinate system at the receiving station $i$ is [3]:

Take ${\alpha }_i$ as the initial value of ${\theta }_i^{\left(k\right)}$ in the following iteration process. Geocentric angle is:

where, $f\left(h\right)=\frac{n_i(r_o+h_i)\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i^{\left(k\right)}}{(r_o+h)\sqrt{n\left(h{)}^2\right(r_o+h{)}^2-{n_i}^2(r_o+h_i{)}^2\mathrm{c}\mathrm{o}{\mathrm{s}}^2{\theta }_i^{\left(k\right)}}}$. The true elevational angle at the receiving station $i$ can be obtained:

The elevation angle error is ${\epsilon }_i^{\left(k\right)}={\theta }_i^{\left(k\right)}-{\stackrel{-}{\alpha }}_i^{\left(k\right)}$. Thus, the apparent elevation angle can be approximated as ${\theta }_i^{(k+1)}={\alpha }_i+{\epsilon }_i^{\left(k\right)}$, combining Eq. (7), and (8), one could obtain ${\theta }_i^{(k+2)}$ from ${\theta }_i^{(k+1)}$. Repeat the above iteration process until the required accuracy is satisfied or the update of parameter is within the given tolerance. Then we have ${\theta }_i={\theta }_i^{\left(m\right)}$. Note that when atmospheric refraction is taken into account, the TDOA curve propagation is $t_{gi1}=t_{gi}-t_{g1}$, when atmospheric refraction is not considered, the TDOA of linear propagation is $t_{oi1}=t_{oi}-t_{o1}$, the model error on TDOA is thus given by $\mathrm{\Delta }{t}_{i1}=t_{gi1}-t_{oi1}$, ($\mathrm{}i=\mathrm{2,3},\cdots N$).

3. Influence of atmospheric refraction

When atmospheric refraction is not considered and is considered, the formulas of positioning measurement are:

The formulas used in actual positioning measurement is:

In actual positioning, the TDOA is measured under atmospheric refraction conditions, while the TDOA formula is used the linear propagation formulas, so the estimation error of position will certainly be caused. It is assumed that the position of the receiving stations are:

The positioning measurement errors of radiation source with elevation of 10 km, 12 km and 14 km are analyzed. The direction of motion is 30° from north to west in the local coordinate system. The simulation results are shown in Fig. 2. In the case of the same elevation, the smaller the elevation angle is, the larger the positioning measurement errors are, the higher the elevation, the larger the positioning measurement error.

Fig. 1Sketch of atmospheric refraction

Fig. 2Influence of atmospheric refraction on position error

4. Position error correction

${\mathbf{d}}_g={\left[d_{g21},d_{g31}\right]}^T$ is the apparent distance difference, ${\mathbf{d}}_g^o={\left[d_{g21}^o,d_{g31}^o\right]}^T$ is the true apparent distance difference, Then the measurement models of TDOA can be expressed as ${\mathbf{d}}_g=c{\mathbf{t}}_g={\mathbf{d}}_g^o+\mathbf{n}$, where, $\mathbf{n}={\left[n_{21},n_{31}\right]}^T$is the corresponding measurement noise, ${\mathbf{t}}_g={\left[t_{g21},t_{g31}\right]}^T$ is the curve arrival time difference. The noise variable is distributed independently and the mean value is zero. Its covariance matric is $\mathrm{E}\left[{\mathbf{n}}^T\mathbf{n}\right]={\mathbf{Q}}_t$, $\mathrm{E}$ is the expectation.

When there is an elevation constraint condition, the source position must be on the elevation constraint surface. In the three-station TDOA position system, there are only two TDOA formulas. There must be one point in the iteration on the elevation constraint surface to stabilize the difference of TDOA error. Therefore, the iterative correction method in the three-station TDOA position system will converge to the optimal point [4].

The specific process of iterative correction algorithm is as follows,

Step 1: According to the measured value of TDOA ${{\mathbf{t}}_g}^{\left(j\right)}$ ($j$ is the iteration number, and ${{\mathbf{t}}_g}^{\left(1\right)}$ is the initial measured value of TDOA), the position ${\mathbf{u}}^{\left(j\right)}$ of the radiation source is calculated by using the positioning Eq. (9);

Step 2: According to the formulas of the Section 2, the TDOA ${{\mathbf{t}}_g}^{\left(j\right)}$ of curve propagation and the TDOA ${{\mathbf{t}}_o}^{\left(j\right)}$ of straight line propagation at position ${\mathit{u}}^{\left(j\right)}$.

Step 3: Calculating the TDOA error at position ${\mathbf{u}}^{\left(j\right)}$, get $\mathrm{\Delta }{\mathbf{t}}^{\left(j\right)}={{\mathbf{t}}_g}^{\left(j\right)}-{{\mathbf{t}}_o}^{\left(j\right)}$.

Step 4: Analyze the calculation results and decide whether to continue iteration the next step.

(1) if $j=1$, the new TDOA is calculated as ${{\mathbf{t}}_g}^{(j+1)}={{\mathbf{t}}_g}^{\left(j\right)}-\mathrm{\Delta }{\mathbf{t}}^{\left(j\right)}$, and return Step1 to continue the calculation.

(2) if $j>1$, calculating the difference between the current TDOA error and the previous TDOA error, get $\mathrm{\Delta }{{\mathbf{t}}_c}^{(j-1)}=\mathrm{\Delta }{\mathbf{t}}^{\left(j\right)}-\mathrm{\Delta }{\mathbf{t}}^{(j-1)}$. Then the decision is made: if $\mathrm{\Delta }{{\mathbf{t}}_c}^{(j-1)}<\mathrm{\Delta }{\mathbf{t}}_{TH}$ ($\mathrm{\Delta }{\mathbf{t}}_{TH}$ is the set threshold), the iteration is stopped and the final positioning result $\mathbf{u}$ is output. Otherwise, the new TDOA is calculated as ${{\mathbf{t}}_g}^{(j+1)}={{\mathbf{t}}_g}^{\left(j\right)}-\mathrm{\Delta }{\mathbf{t}}^{\left(j\right)}$, and go back to Step1 to continue the iteration.

5. CRLB of position estimation

The CRLB has extended the bound to incorporate constraints [5]. Differentiation of Eq. (5), ${\mathbf{H}}_{1}d\mathbf{u}=cd\mathbf{T}$, where:

Taking the position $\mathbf{u}$ derivative of Eqs. (3) and (7). Because the position-dependent parameters are only the integral upper limit $h_o$ and the apparent elevation angle ${\theta }_i$, we can get:

Multiplying the two ends of Eq. (11) by $n_i(r_o+h_i)\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i$, subtract from Eq. (10), get:

From cosine and sine theorem, we have $d_{oi}^2=(r_o+h_o{)}^2+(r_o+h_i{)}^2-2(r_o+h_o)(r_o+h_i)\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_i$ and $(r_o+h_o)\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{Ti}=(r_o+h_i)\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_i$, then we can get:

Substituting Eq. (13) into Eq. (12), and the elevation error is ${\epsilon }_{Ti}={\alpha }_{Ti}-{\theta }_{Ti}$, get:

It can be seen that when there is no atmospheric refraction, taking the position $\mathit{u}$ derivative of Eq. (4), from Fig. 1, $h_o=‖\mathit{u}‖-r_o$, there is:

Substituting Eq. (15) into (14), we can get:

CRLB of radiation source position estimator $\mathbf{u}$ can be expressed as:

6. Simulation results

The position of the receiving station is shown above and Monte Carlo is 500 times. Assuming the position of the radiation source is ${\mathbf{u}}^o=[98.5^{\mathrm{}\circ }\mathrm{E},38.5^{\mathrm{}\circ }\mathrm{N},12000\mathrm{}\mathrm{m}{]}^T$, the elevation constraint condition is $H=$ 12000 m. The results show that the TDOA of linear propagation without considering atmospheric refraction and the TDOA of the curve propagation with considering atmospheric refraction are respectively as follows:

The errors of TDOA caused by atmospheric refraction are:

As shown in Fig. 3, In the case of small time difference noise, the correction effect is very obvious, because the positioning error caused by atmospheric refraction plays a major role. The corrected positioning error is close to CRLB, which also verifies the effectiveness of the modified algorithm.

Fig. 3The curve of position error

7. Conclusions

Under the condition of low elevation and long distance, the time difference location errors caused by atmospheric refraction are large. In this paper, the iterative method is used to correct the positioning error of atmospheric refraction TDOA. The effectiveness of the algorithm is verified through simulation analysis.