Fusion of Land-Based and Satellite-Based Localization Using Constrained Weighted Least Squares

Combining multiple devices for localization has important applications in the military field. This paper exploits the land-based short-wave platforms and satellites for fusion localization. The ionospheric reflection height error and satellite position errors have a great impact on the short-wave localization and satellite localization accuracy, respectively. In this paper, an iterative constrained weighted least squares (ICWLS) algorithm is proposed for these two kinds of errors. The algorithm converts the nonconvex equation constraints to linear constraints using the results of the previous iteration, thus ensuring convergence to the globally optimal solution. Simulation results show that the localization accuracy of the algorithm can reach the corresponding Constrained Cramér–Rao Lower Bound (CCRLB). Finally, the localization results of the two methods are fused using Kalman filtering. Simulations show that the fused localization accuracy is improved compared to the single-means localization.


Introduction
Source localization has important applications in wireless sensor networks, radar, navigation, and other fields [1][2][3][4].Typical localization parameters include direction of arrival (DOA), time difference of arrival (TDOA), frequency difference of arrival (FDOA), and so on.Over-the-horizon (OTH) localization is very common for emitters such as aircrafts and warships in the military field.For important military targets, it is very necessary to utilize multiple means and multiple platforms to locate and surveillance them, and the positioning accuracy is higher compared to single means.For targets such as ships, the signals they send can usually be received by land-based shortwave stations as well as satellites.Therefore, it is necessary to study accurate and effective multi-platform localization algorithms [5].
For short-wave stations, the frequency range of the received signal is from 3 MHz to 30 MHz, and the signal usually arrives at the observatory station after reflection from the ionosphere.The state of the ionosphere has a great influence on the localization results, so it is necessary to use an appropriate ionosphere model.Recently, TDOA-based localization for short waves as a new method has attracted the interest of many scholars.A mathematical solution for short-wave TDOA-based localization was proposed in the study of Jain A [6]. Huang S [7] proposed a gradient-type algorithm based on the Quasi-Parabolic (QP) ionosphere model.Although the QP model can describe the ionospheric state more accurately and can give analytic equations for a ray path, it involves a number of ionospheric parameters, which brings more errors.Although the ionospheric virtual height (IVH) reflection model [8][9][10] is relatively simple, it contains the ionospheric reflection height as the main parameter, which can accurately reflect the state of the ionosphere and can effectively simplify the localization problem.
However, ionospheric reflection height errors can also have a significant impact on the localization results, so the influence on the localization algorithm needs to be considered.
The time-varying property of the ionosphere leads to a large error in the estimation of the TDOA and thus low accuracy in TDOA-based localization for short-wave sources.In contrast, DOA-based localization for short-wave emitters has gained maturity, and the localization results are more stable.Wang [9] proposed a new algorithm for short-wave sources using orthogonal triangular decomposition (QR decomposition) based on the IVH model.The positioning accuracy of this algorithm can reach the corresponding Cramér-Rao Lower Bound (CRLB).
In addition, satellite-based localization is a high-precision method that typically utilizes the TDOA.A closed-form algorithm for TDOA-AOA hybrid localization was proposed in [11] based on weighted least squares (WLS).Ref. [12] analyzes the effect of earth constraints on multi-satellite localization and obtains the final result by weighting the positioning results under different constraints through weighting coefficients.However, the position of satellites is usually inaccurate due to orbital errors, altitude errors, and other reasons [13].Therefore, satellite position errors need to be considered in order to improve the localization accuracy [14].Ref. [15] corrects for satellite errors using the calibration source and proposed a Lagrange algorithm.Ref. [16] proposed a moving horizon estimation (MHE)-based technique to achieve source tracking using TDOA-FDOA measurements from multiple satellites.
The surveillance of enemy warships is an important application in the military field, which usually utilizes short waves or satellites for communication.Therefore, highprecision positioning is of great significance.However, short-wave positioning is greatly affected by ionospheric interference, while satellite positioning is affected by orbital errors.Combining multiple platforms to localize important sources can effectively improve the localization accuracy, which is a very effective means in practice.By combining landbased short-wave stations with satellites for localization, a hybrid localization algorithm based on a penalty function was proposed in [5].But, Ref. [5] does not exploit the elevation information and does not consider the position errors of the satellite.Moreover, the equation constraints in the localization optimization problem are nonconvex, and thus the method can fall into a local optimal solution, which leads to a degradation of the localization accuracy.
In this paper, an iterative constrained weighted least squares (ICWLS) algorithm is proposed for sources such as warships transmitting multiple types of signals, which are localized using land-based short-wave stations and satellites.By introducing auxiliary variables, the DOA-based localization using short-wave stations and the TDOA-based localization using satellites can be modeled as optimization problems with two quadratic equation constraints, respectively.The quadratic nonconvex equation constraints are linearized using the results of the previous iteration [17,18], which converts the nonconvex constraints to linear constraints and ensures the algorithm can converge to the global optimal solution.Due to the difference in auxiliary variables between the DOA pseudolinear equation and the TDOA pseudo-linear equation, a unified pseudo-linear equation cannot be established.Therefore, it is proposed in this paper that the localization results of the two methods are fused using Kalman filtering to improve the localization accuracy.
The mathematical symbols used in this paper and explanations are shown in Table 1.This paper is organized as follows: Section 2 introduces the localization scenario and the measurement model; Section 3 establishes the corresponding pseudo-linear equations; Section 4 derives the Constrained CRLB (CCRLB); Section 5 gives the proposed localization algorithms and the method for fusion; Section 6 shows the simulation results; and Section 7 is the conclusion.

Symbols Explanations
[•] T Matrix or vector transpose The inverse of the matrix trace(•) The ith element of the vector respectively.There is a ground station for satellites, with a longitude and latitude of ω t,1 and ω t,2 , respectively.Based on the Earth ellipsoid model, (1) provides the transformation relationship between the geodetic coordinate system [ω 1 , ω 2 , H] T (ω 1 and ω 2 denote the longitude and latitude, respectively, and H is the altitude) and the space rectangular coordinate system [19]: α = R e / 1 − (esin(ω 1 )) 2 , where R e = 6378.137km is the radius of the Earth's equator and e = 1 − R 2 p /R 2 e denotes the eccentricity.R p = 6356.752km is the polar radius of the Earth.Therefore, based on (1), the position vector of the land-based station and satellite can be s o j = s o j,x , s o j,y , s o j,z T .The source position vector is represented as u = u x , u y , u z T , and the ground station position vector is t = t x , t y , t z T .Since the source is located on the Earth surface, u satisfies the following equation: where

DOA Measurement Model for Land-Based Station
Figure 1 is the IVH model, which shows the short-wave signal sent by the source reaching the land-based station after being reflected by the ionosphere.θ o i represents the true value of the azimuth and φ o i is the true value of elevation.θ o i is defined as the angle between the projection of the incidence direction on the surface plane and the local north direction of the station.φ o i is defined as the angle between the incidence direction and the surface plane of the station.h o i is the ionosphere reflection height and R o = 6371.393km is the average radius of the Earth.
Combining all position vectors of land-based stations yields: ] T .The ionosphere reflection height is often difficult to measure accurately and is therefore usually modeled as measurements with errors:

The reflection height measurement vector is thus given as
and T denotes the true value of the ionosphere reflection height.
n h = δh 1 , δh 2 , . . ., δh N 1 T denotes the measurement error, and its covariance matrix is denoted as In the local coordinate system of the station, the target coordinate vector u can be transformed:

Ionosphere
The origin of the local coordinate system is the station.Therefore, θ o i can be given as According to (4), the true value vector of azimuth can be given as In practice, there will be measurement noise, so the following measurement vector can be obtained as: where n θ = δθ 1 , δθ 2 , ..., δθ N 1 T denotes the noise vector, which follows a zero-mean Gaussian distribution, and the covariance matrix is

△ABC, based on the sine theorem, yields
where 6) can be converted into Based on ( 7), the true value vector of the elevation angle is expressed as T , and further, the measurement vector is given as ) T denotes the elevation noise vector, and the covariance matrix of Combining ( 5) and ( 8) gives the DOA measurement model for land-based stations as where

TDOA Measurement Model for Land-Based Station
The signals emitted by the source are forwarded by satellites to the ground station, which receives the signals sent by all the satellites and can obtain the TDOAs.Figure 2 shows the schematic diagram of the satellite forwarding the signal to the ground station.Satellites are in high-speed motion and thus the positions are also in error, so the satellite positions are modeled as measurements with errors: T is the random error.Thus, the satellite position vector is obtained as s t = s o t + n s , where From Figure 2, the propagation length can be expressed as Taking the first satellite as the reference satellite, the TDOA can be obtained as where c is the signal propagation velocity and is a known constant.Thus, the TDOA shown in (11) can be converted to the range difference of arrival (RDOA) shown in (12): The measurement vector of RDOA is denoted as where n r = δr 21 , δr 31 , ..., δr N 2 1 T is the measurement noise with the covariance matrix Q r .

DOA Pseudo-Linear Equation for Land-Based Station
The linear equation of azimuth is first established.Based on (4), we have Simplifying (14) gives Combining ( 15) of all land-based stations obtains where Converting (7) gives Substituting Squaring and simplifying both sides of (20) simultaneously yields where is the only positive root of the equation shown in (21).Therefore, we have Squaring both sides of ( 22) simultaneously yields The pseudo-linear equation for elevation is obtained as follows: where . . .

TDOA Pseudo-Linear Equation for Satellite
Let Shifting the terms and squaring both sides simultaneously of (10) yields Expanding and simplifying (30) yields Therefore, the pseudo-linear equation for the satellite can be given as where . . .
Remark 1.Any two equations having the form ( 27) and (32) can be combined into one equation by expanding the matrix dimension.However, this is a simple stacking of matrices and vectors and does not take advantage of the connection between the two.Since the auxiliary variables utilized in ( 27) and (32) are different, the unknown variables in the two equations are also different.Thus, it is not possible to directly combine ( 27) and (32) into one equation.Although it is possible to combine ( 27) and (32) to obtain an equation, the combination of the two equations is a simple matrix combination that does not take advantage of the connection between the two equations.

CCRLB
The CRLB gives the theoretical lower bound that unbiased estimation can achieve, and the CRLB can be given by the inverse of the Fisher information matrix (FIM).

CCRLB for Land-Based Station
Under the Gaussian noise assumption, the log-likelihood function of z o with respect to ξ 1 = u T , (h o ) T T can be expressed as [20,21] ln λ a is a constant.The corresponding FIM is denoted as where The CRLB for DOA-based localization can be given as where gives the CRLB for source localization.Since the source position vector satisfies (2), the CCRLB with the source constraint can be expressed as [9] where Λ = [Λ 1 u] T .

CCRLB for Satellite
In the presence of errors in the satellite position, the TDOA-based log-likelihood function can be expressed as λ t is also a constant and ξ 2 = u T , (s o t ) T T .The FIM can be given as where The CRLB for TDOA-based localization is obtained as where X TDOA 1 ∈ R 3×3 .And, the CCRLB can be given as

CCRLB for DOA-TDOA Hybrid Localization
Define the unknown vector ξ = u T , (h o ) T , (s o t ) T T and the measurement vector is The covariance matrix of n is Q z .The FIM of hybrid localization can be given as where Therefore, the CRLB and CCRLB can be expressed as where Λ = blkdiag Λ, 0 (N 1 +3N 2 )×(N 1 +3N 2 ) .Since in the practical process, only the measurement with errors is available, the following error equation can be established:

Proposed Method
Substituting the measurements with errors into (25) gives y a and G a .Performing first-order Taylor expansion for y a and G a yields where Substituting ( 52) into (51) yields Therefore, the DOA localization problem for land-based stations can be modeled as W a is the weighted matrix, which is defined as As can be seen from ( 55), the cost function is convex and the equation constraints are nonconvex.The equation constraints are written in matrix form as where D a is related to u, so the estimate of u can be substituted into D a to obtain D k a during the iteration process.J a (η a ) in ( 55) can be expressed as y T a W a y a is a constant with respect to η a .Therefore, substituting the results of the kth iteration into (57), (55) can be converted to min where Ga = G T a W a G a and ỹa = G T a W a y a .k is the iteration number.D k a and d k a can be obtained by substituting the kth iteration result u k a into (58).As can be seen in ( 60), the cost function is still convex.However, by using the results of the previous iteration, the nonconvex quadratic equation constraints can be converted to linear constraints, which can ensure that the algorithm converges to the globally optimal solution.Thereby, the optimal solution of ( 60) is [17,22]: where Therefore, the iteration results can be updated by Equation (62): 0 < κ a < 1 denotes the weighting factor.The initial estimation can usually be obtained by making W a = Q −1 z and then obtaining an initial estimation based on η0 The DOA localization algorithm for land-based stations is summarized in Table 2.

Covariance Matrix of DOA Result
The DOA localization error is defined as ηa .Thus, ∆η a = ηa − η o a is the optimal solution to the following problem: where Λ1 = blkdiag(Λ 1 , 0).The constraints in (64) are derived from ( 2) and (29), respectively.∆η a can be represented as [5] ∆η where The covariance matrix of the source estimate can be obtained as

Computational Complexity of Land-Based Localization
Table 3 gives the number of multiplications required for some elements.
Table 3. Computational complexity of land-based localization.

Localization Method for Satellite
Substituting the measurements into (32) gives Performing the first-order Taylor expansion of y t and G t gives where Substituting ( 69) into (68) gives Therefore, the TDOA localization for the satellite can be modeled as W t is defined as Similarly, to ensure that the algorithm converges to the global optimal, a transformation of the nonconvex equality constraint is required.Writing the two equality constraints in matrix form gives where Substituting the results of the dth iteration into (72) gives min where Gt = G T t W t G t and ỹt = G T t W t y t .The optimal solution of (76) is expressed as where The k + 1th iteration result is obtained by the following equation: where 0 ≤ κ t ≤ 1.
In summary, the TDOA localization algorithm for satellites is summarized in Table 4.

Covariance Matrix of TDOA Result
Define the TDOA estimation error as ∆η t = ηt − η o t .Therefore, ∆η t is the optimal solution to the following problem: ∆η t is defined as where T , 0 .
The covariance matrix of ∆η t can be given as Therefore, the covariance matrix of the source estimate is

Computational Complexity of Satellite-Based Localization
Table 5 gives the number of multiplications required for some elements.

Localization Result Fusion
ûa and ût are asymptotically unbiased estimates which, respectively, obey the following distribution: Based on Kalman filtering [23], the DOA-based localization results and TDOA-based localization results can be fused to obtain where K denotes the Kalman gain, which is defined as [24,25] K = trace(cov( ûa )) trace(cov( ûa )) + trace(cov( ût )) The details of the proposed method are given in Figure 3.The computational complexity of the proposed method can be obtained by combining Tables 3 and 5, although the proposed method has an increased computational complexity compared to land-based localization or satellite-based localization.However, since N 1 and N 2 are both small, the increase in computation is acceptable.

Computational Element
Computational Complexity

Simulation
Simulations are set up with a total of 10 observation stations, including 5 land-based stations and 5 satellites.The longitude and latitude of the land-based station and the corresponding ionospheric reflection heights are shown in Table 6.The longitude and latitude of the satellite and its altitude are shown in Table 7.The ground station is located at (110.3 • E, 25.7 • N).The source longitude and latitude are (126 • E, 37 • N), respectively.Remark 2. The main application scenario of this paper is the military field, and the targets are mainly military targets such as warships.In addition, limited to the current level of the authors, no publicly available short-wave dataset has been found.The method proposed in this paper can not be tested on public datasets.Therefore, this paper uses simulation to verify the performance of the proposed method.
Assuming that the noise between each parameter is uncorrelated, the covariance matrix can be set as , where R N 2 −1 denotes the (N 2 − 1) × (N 2 − 1) matrix with a diagonal element of 1 and the rest of the elements are 0.5.The root mean square error (RMSE) is used as the localization accuracy criterion and is defined as follows: where M denotes the number of Monte Carlo experiments.Set δ a = δ t = 10 −6 , κ a = 0.9 k , and κ t = 0.9 k , where k denotes the number of iterations.The maximum number of iterations is set to 50.
As can be seen from Figure 4, the localization results are consistent with the ellipse, proving the effectiveness of the proposed ICWLS algorithm.In addition, the ellipse area of the fusion result is smaller compared to the single-method localization, which demonstrates that the fusion method improves the localization accuracy.

Localization Accuracy
This subsection simulates the localization accuracy of the proposed algorithm under the influence of different errors.Set σ θ = 0.2σ 1 , σ φ = 0.3σ 1 , and σ r = 0.3σ 1 , where σ 1 is the localization parameter error.Let σ h = 0.5σ 2 and σ s = 0.2σ 2 , where σ 2 is the systematic error parameter.Remark 3. The effects of the signal-noise ratio (SNR) and propagation channel will ultimately affect the estimation accuracy of the DOA and TDOA.Therefore, the error parameters σ 1 and σ 2 can ultimately reflect the individual errors.This paper focuses on the localization algorithm and thus mainly simulates the effect of the DOA error and TDOA error on the positioning accuracy.Firstly, set σ h = 3 km and σ s = 2 km and simulate the localization accuracy with σ 1 .Figure 5 is the localization result.Then, let σ θ = 1.5 • , σ φ = 2.25 • , and σ r = 0.5 km, and Figure 6 is the corresponding result.
From Figures 5 and 6, the positioning accuracy of the proposed ICWLS algorithm for land-based DOA localization and satellite TDOA localization can reach the corresponding CCRLB.In addition, it can be seen that the fusion result can significantly improve the localization accuracy, which further proves the effectiveness of the proposed fusion method.Firstly, set σ h = 3 km and σ s = 2 km and simulate the localization accuracy with σ 1 .Figure 5 is the localization result.Then, let σ θ = 1.5 • , σ φ = 2.25 • , and σ r = 0.5 km, and Figure 6 is the corresponding result.
From Figures 5 and 6, the positioning accuracy of the proposed ICWLS algorithm for land-based DOA localization and satellite TDOA localization can reach the corresponding CCRLB.In addition, it can be seen that the fusion result can significantly improve the localization accuracy, which further proves the effectiveness of the proposed fusion method.As can be seen from Figures 7 and 8, both the RMSE and CCRLB increase as the error increases, and the distribution is more dispersed.In addition, the distribution of the RMSE is basically the same as that of the CCRLB, and the localization results do not produce abnormal distribution values, which indicates that the proposed algorithm has good robustness in terms of the source location.

Conclusions
This paper investigated the fusion localization of land-based stations and satellites.And, the corresponding ICWLS algorithm is proposed for DOA-based localization and TDOA-based localization, respectively.Firstly, the pseudo-linear equations are established by using auxiliary variables, and then the localization problem is modeled as an optimization problem with two quadratic equation constraints.To ensure that the algorithm can converge to the global optimal solution, the results of the previous iteration are substituted into the equation constraints so that the nonconvex quadratic constraints can be converted into linear constraints.Simulations show that the localization accuracy of the proposed method can reach the corresponding CCRLB.In addition, this paper proposed to fuse the DOA-based localization and TDOA-based localization results using Kalman filtering.The simulation results show that the accuracy of fusion result is higher than the DOA-based localization and TDOA-based localization.
Although the proposed method can improve the positioning accuracy, the method requires the ability to receive signals from different frequencies.In addition, the method does not take into account the errors caused by the ionospheric multi-path effect and other effects.These will be studied in depth in the next step.In addition, we will further consider the practical application of the proposed method.

Figure 2 .
Figure 2. Schematic diagram of a satellite relaying the source signal.

Figure 3 .
Figure 3. Block diagram for the proposed method.
(a) X-Y plane (b) X-Y plane enlarged (c) X-Z plane (d) X-Z plane enlarged (e) Y-Z plane (f) Y-Z plane enlarged

Figure 4 .
Figure 4. Localization results and uncertainty error ellipses.

Figure 4 .
Figure 4. Localization results and uncertainty error ellipses.

Figure 5 .
Figure 5.The localization accuracy varies with localization parameter error.

Figure 6 .
Figure 6.The localization accuracy varies with systematic error.

6. 3
. Robustness A total of 10 sources are randomly selected within the range of [120 • E∼130 • E] × [20 • N∼30 • N].The localization accuracy of the proposed algorithm and the corresponding CCRLB distribution are tabulated.First, set σ h = 3 km and σ s = 2 km, and the corresponding results are shown in Figure 7.Then, set σ θ = 1.5 • , σ φ = 2.25 • , and σ r = 0.5 km, and the results are shown in Figure 8.
(a) CCRLB for land-based station (b) RMSE for land-based station (c) CCRLB for satellite (d) RMSE for satellite (e) CCRLB for hybrid localization (f) RMSE for hybrid localization

Table 1 .
Mathematical symbols and explanations.

u,1 and ω u,2 ,
We assume a total of N 1 + N 2 stations, including N 1 land-based stations and N 2 satellites.These stations are used to locate the source on the Earth's surface.The longitude and latitude of the source are denoted as ω respectively.The station is represented byO 1 , O 2 , . .., O N 1 , O N 1 +1 , . .., O N 1 +N 2 ,where O 1 , O 2 , . . ., O N 1 represents the land-based stations and O N 1 +1 , . . ., O N 1 +N 2 represents satellites.The longitude and latitude of the jth station are ω j,1 and ω j,2 , 1

Table 2 .
Summary of land-based localization.Set k = 0 and choose appropriate δ a > 0 and κ a .Let W a = Q −1 z to obtain the initial result;Step 2: Set k = k + 1. Substituting ηk−1 a .If not, go to Step 2.

Table 4 .
Summary of satellite-based localization.Set k = 0 and choose appropriate δ t > 0 and κ t .Let W t = Q −1 r to obtain the initial result;Step 2: Set k = k + 1. Substituting ηk−1 t .If not, go to Step 2.

Table 6 .
Longitude and latitude (E: east, N: north) and the corresponding ionospheric reflection height for land-based stations.

Table 7 .
Longitude and latitude (E: east, N: north) and the corresponding altitude for satellites.