Quasi-Real RFI Source Generation Using Orolia Skydel LEO Satellite Simulator for Accurate Geolocation and Tracking: Modeling and Experimental Analysis

Abstract

Accurate geolocation and tracking of Radio-Frequency Interference (RFI) sources, which affect wireless and satellite systems such as Global Navigation Satellite Systems (GNSS) and Satellite Communication (SatCom) systems, are considered to be a significant issue. Several studies connected to civil and military operations on this issue have been investigated recently. The literature review has surveyed many algorithm simulations for optimizing geolocation and target-tracking estimation. Although most of these algorithms have their own advantages, they have weaknesses, such as accuracy, mathematical complexity, difficulties in implementation, and validation in the real environment, etc. This study has been concerned with investigating the accuracy of geolocation and tracking under high speed and powerful rotation using extracted data from the Orolia Skydel simulator, which simulates the space environment involving Low Earth Orbit (LEO) satellites as sensors and Unmanned Aerial Vehicles (UAV) as RFI emitters. Various scenarios modeled using the Orolia Simulator for quasi-real dynamic trajectories of LEO satellites have been created. The assumed approaches have been verified by Cramer–Rao Lower Bound (CRLB) and Posterior CRLB (PCRLB) to determine the increase in Root Mean Square Error (RMSE) value. The simulation scenarios have been performed using the Monte Carlo iteration. Eventually, the overall achieved results of the considered approaches using data acquired from the Orolia Simulator were presented and compared with theoretical simulation.

1. Introduction

Radio-Frequency Interference (RFI) is an unwanted signal. It can reduce the reliability and efficiency of the communication networks in wireless systems [1,2]. This issue has led to the need for location identification and RFI source tracking. From that, it can notify the organization responsible for RF intervention, which can be mitigated with the help of sophisticated techniques. For this reason, researchers and scientists have paid more attention to finding a suitable technique for precision tracking and localization of dynamic RFI, as well as creating different scenarios using a Software-Defined Radio (SDR) system as a quasi-real implementation on proposed algorithms. Even though extensive studies on geolocation and target-tracking techniques have been carried out based on various electromagnetic propagation parameters, which use various statistical methods by MATLAB simulation scenarios, none of the literature covers modeling and experiments using the very advanced and unique kind simulator called “Orolia Skydel simulator”. Various experimental work has been implemented for geolocation and emitter-tracking to emulate navigation and communication satellite signals using other simulators. For example, ref. [3] created a model and analyzed the performance geolocation technique based on Satellite Tool Kit (STK) scenarios. In [4], the authors offered two constellations based on the GPS positioning output by OEM628 of BDStar Navigation Company. The data were processed using the signal software-processing platform, and the fusion positioning was done using opportunity signals (SOPs). In another article, [5], the pirent^TM GSS6700 GNSS signal simulator (Spirent Communications plc, Paignton, UK) was paired with an interference signal generator through a Spirent^TM GSS8366 device to construct quasi-real dynamic trajectories. Signal-power degradation, navigation solution availability, Dilution Of Precision (DOP), and positioning accuracy were all evaluated as part of the performance evaluation.

Table 1. Comparison between the modeling and experiment proposed in [3] with our proposed modelling and experiment.

No	Ref. [3]	Our Proposed Modelling and Experiment
1	They have modeled quasi-real geostationary satellite with three antennas as sensors to receive the signal strengths from stationary Earth station as emitter using STK.	We have modeled quasi-real dynamic trajectories of LEO satellites as sensors and Unmanned Aerial Vehicle (UAV) as an emitter using Orolia Skydel simulator.
2	They worked to geolocate the stationary emitter using the RSS algorithm.	We worked to geolocate and track a dynamic emitter using Cross Ambiguity Function (CAF) and High-degree Nonlinear Tracking Filters based on TDOA/FDOA.
3	They verified the approach performance by increasing latitudes at varying contour widths.	We verified the approach performance based on different measurements of hybrid TDOA/FDOA as well TDOA and FDOA individually.

presents a comparison between the approach STK modeling executed in [3] and those executed by our approach.

In this study, the Orolia Skydel simulator is used elegantly to generate In-phase/ Quadrature (I/Q) data files of RFI signals received by LEO satellites. Figure 1 illustrates the considered model of proposed scenarios.

The simulated scenario is very well adapted for geolocation tracking of maneuver RFI sources. Many vital projects rely on durable and reliable GNSS signals to deliver extremely accurate Positioning, Navigation, and Timing (PNT) data in the next generation of satellite-enabled applications. Orolia GNSS testing and simulation services to assure the system’s performance, resiliency, and accuracy in complex GNSS space-based applications, also generating low earth orbit simulation [6].

Time of Arrival (TOA), Time Difference of Arrival (TDOA), Frequency of Arrival (FOA), Frequency Difference of Arrival (FDOA), Angle of Arrival (AOA), Received Signal Strength (RSS), and Power Difference of Arrival (PDOA), etc. are electromagnetic propagation parameters that are used by various approached of geolocation, and tracking [7,8,9]. The TOA measurement is related to knowing the distance between the source and the sensor and radio propagation velocity, where the radio propagation velocity is well known as 2.998 × 10⁸ m/s. The common lack of synchronization between the source and the receiver is one practical obstacle in TOA measurements. In other meanings, receivers are frequently unaware of the specific start time of source transmission. Because of the uncertainty surrounding the transmission instant start time, all received TOA measurements have a common temporal offset, which can contribute to severe localization accuracy [10]. Subtraction of pairwise TOA measurement results in correlated noise in TDOA. The sensors synchronization is taken into consideration in the TDOA model. In a two-dimensional (2D) scenario, to achieve better precision for TDOA-based geolocation, at least three sensors are needed for receiving similar signals that broadcast from the emitter [11]. The FOA model depends on the carrier frequency, radio propagation velocity, and Doppler shift between emitter and sensor when both or one of them is moving [12]. The FDOA is a difference in pairwise FOA measurement caused by the Doppler shift of the received signal due to differences between sensors and emitter velocities [13].

A signal observation that might be received from a sensor is not enough to compute a better estimation of the emitter source. In the TDOA method, the emitter might estimate over many scattered points of a hyperbola curve. This problem requires increasing the number of sensors or sensing points that mainly target the geolocation measurement of received signals. From that, more than one hyperbola will produce a precise target location estimation [14,15]. In addition, combining two or more geolocation measurement models will achieve the best measurement of electromagnetic propagation parameters and reduce the sensors or sensing points needed [16,17,18,19,20,21].

From the literature review, numerous geolocation and target-tracking algorithms based on electromagnetic propagation parameters have been considered [2,22,23]. In [24], the authors have studied the accuracy of emitter geolocation using the CAF approach. They proposed an approach to eliminate mirrored locations in tri-satellite TDOA localization based on TDOA measurements taken at many times to attain precision. In addition, the authors in [25] addressed the problem of stationary ground-based location using one transmitter (emitter) and two airborne mobile receivers. They have explored the effect of spatially multiplexed Multiple-Input-Multiple-Output (MIMO) signals on the CAF calculation of TDOA and FDOA, and the consequent effect on geolocation; this method did not carry out any modeling and experimenting of the scenario before arriving at the solution.

Moreover, Dexiu Hu et al. [20], proposed an algebraic method using combined TDOA and FDOA with a differential Doppler rate that increases the moving source geolocation accuracy. The method proposed by the author relied on the pseudo-linear set of equations and a two-step weighted least squares (WLS) estimator. In [26], the authors proposed a Linear Least Squares (LLS) approach. This study analyzed the performance of the proposed algorithm, which considered four static sensors to geolocate mobile sources using the LLS method based on the TOA/RSS measurement. The authors of [27] presented a position-estimation method based on a Deep-Learning (DL) algorithm that works directly on TDOA-based locating system channel impulse responses. The authors explained the signal and data possessing in depth and verified the method’s effectiveness in various real-world scenarios. Although their method corresponds to traditional signal-processing methods under line-of-sight conditions, it outperforms earlier methods under significant multi-path propagation.

To estimate the location and tracking of the RF emitter, Kalman filter (KF) and its expansion based on conventional geolocation measurements has been proposed in numerous studies [28,29,30]. In [31], the Extended Kalman filter (EKF) based on TDOA/FDOA was proposed to estimate target position and tracking. The main disadvantage of the EKF is that it relies on pre-estimation to estimate the variance of process and measurement noises, which may be difficult to accurately track the RFI position when predicting with a significant deviation from the actual values (measurement) [32]. In [33], the authors proposed two algorithms to solve the simultaneous mapping and localization (SLAM) problem. The first algorithm is an amortized constant-time coordinate descent algorithm for recovering these state estimates from the information form, and the second is an efficient algorithm for data association in Sparse Extended Information Filters (SEIFs) that takes logarithmic time, assuming that an efficient search tree searches for nearby features. The method was put into practice and compared to the EKF solution. In [34], the authors used EKF to create a unified Factor Graph (FG) architecture in which a practical AOA-based position detector accomplishes the observation process.

Furthermore, the variance of observation errors, which EKF requires, was determined in real time by using both the AOA measurements and the expected target state. Compared to the traditional method, such a dynamic estimating strategy has higher performance robustness, especially when the sensing environment is unstable. Simulations show that the suggested system achieves lower Root Mean Square Error (RMSE) in many assessment settings and has a rapid convergence tendency. The Adaptive Extended Kalman Filter (AEKF) was proposed to reduce error estimation for similar problems [35]. Still, for high non-linearity, the error may increase because of the rotational speed of the emitter. Multiple Quadrature information Kalman filters were proposed in [36]; these can work in high non-linearity with a high-speed mobile sensor. In digital signal-processing, estimation has been considered an important challenge. For that, numerous studies have been deployed for estimations, and KF deemed as a highly powerful tool used for a linear dynamical system to estimate state and deems as a least-squares linear filter that provides an optimal recursive deployment. In [37], a hybrid localization approach based on the Particle Filter (PF) and particle swarm optimization algorithm are presented, focusing on the localization tasks when an a priori environment map is available. In addition, the approach was advocated for Distributed Bayesian filtering (DBF) in [38], where each node operates a distinct particle filter, and the collected consensus is sought on the sensor data alone or in conjunction with intermediate local filter estimates.

However, it is quite challenging to design an optimal filter that functions in a nonlinear system. In addition, a study showed that the cubature Kalman filter (CKF) estimations are better when compared with the estimations of the unscented Kalman filter (UKF) [39]; however, its performances are worse compared with the Gauss–Hermite quadrature filter (GHQF) [40]. Therefore, for better performance, CKFs with random degrees of accuracy were proposed, and therefore it is stated that the high-degree CKFs could estimate better accuracy and stability performance such as GHQF but with less computational cost. High-degree cubature rules are applied to numerical integration and a target-tracking problem to demonstrate the system effectiveness and better performance computations. The measured results of target-tracking simulation illustrated that enough accuracy can be obtained employing fifth-degree CKF rather than EKF, therefore combining the use of third-degree CKF and PF; the performance is much better than GHQF [41].

The main contributions of this paper are: (i) creating various scenarios using Orolia Skydel as a space environment involving LEO satellites to receive emitted RFI signals from UAVs; (ii) applying the assumed scenarios as a quasi-real environment to verify the approaches of optimizing geolocation measurement and high-degree nonlinear tracking filters and compare performance with theoretical simulations.

The paper is structured as follows. The architecture of geolocation and target-tracking using the Orolia Simulator that is considered and presented in Section 2. The mathematical model is detailed in Section 3. Experimental implementation and actual simulation scenarios are explained in Section 4. Experiment results and performance analysis are illustrated in Section 5. Finally, the conclusion and suggestions for future works are provided in Section 6.

2. Proposed Simulation Environment

The multi-LEO satellites (sensors) produce time delay and a frequency Doppler affecting the received signal from the UAV(emitter), which has led us to propose TDOA/FDOA measurement using the Orolia Skydel simulator as a quasi-real environment. TDOA/FDOA measurements can be implemented on a nonlinear tracking filter to track the RFI sources. Figure 2 shows the general diagram of the proposed model and architecture. In this study, we considered the fixed Earth-centered Earth-Fixed (ECEF) coordinates [42,43].

2.1. Geolocation and Target-Tracking Architecture

In the Orolia Skydel simulator, setting up LEO satellite simulation concerning tracking UAV is possible, which spuriously generates the RFI source of interference. With the help of “Orolia”, an advanced software-defined GNSS simulator installed in the LASSENA laboratory, we can generate the simulation for different instances localization algorithms are recorded in raw I/Q files [44]. Moreover, the information of transmitter and receiver can also be downloaded using logging files for each epoch as a reference file, which will be useful to fix the parameters of UAV and LEO satellites coordinates and velocity, etc. in the source localization algorithm implemented in MATLAB. Initially, the number of sources of RFI in the form of UAV and the LEO satellites are set in separate instances by assigning a window to a Radio 0, 1, 2, etc. To the best of our knowledge, this is the first time, this simulator is used to simulate and record RFI signal geolocation and emitter-tracking problems.

2.2. Modeling of Emitter, Sensor and Radio Assign

In the beginning, we need to define the radio and select the type of signal, which could be GNSS Upper L-band, GSSN Lower L-Band, or Interference. After that, starting with Instance-1, we can specify the trajectory of the dynamic emitter from the Interference section. In the study, we consider the UAV as an RF emitter. After that, we can specify the trajectory of the UAV(as a dynamic emitter) from the interference section. Next, we can specify the reference power level (that we considered −90 dBm). From the transmitter section, we can select signal and trajectory. We can observe the power level at the receiver side (LEO satellite) specified at the transmitter side. Additionally, in the transmitter section, we can ignore propagation loss. In the trajectory, we can choose a fixed position of service that we can put as a circle shape or import the trajectory of the dynamic emitter as Keyhole Markup Language (KML) or Comma-Separated Values (CSV) file format with timestamps.

Moreover, the information of transmitter and receiver can also be downloaded using logging files for each epoch as a reference file, which will be useful to fix the parameters of UAV and LEO satellites coordinates and velocity etc. From the Settings-Interference-Transmitter 1, and Signal section, we can select the type of signal we want to transmit. It could be Continuous-Wave (CW), Chirp, Pulse, or Binary Phase-Shift Keying (BPSK) signal, see Figure 3.

As well, the central frequency of the BBSK signal, which was considered in this study as 1575.42 MHz with 2.046 MHz bandwidth as shown in Figure 4.

From the antenna section, we can import the antenna pattern of the emitter, or we can select None to avoid interfering with the power of the signal, see Figure 5.

After configuring the transmitter, we need to configure the receiver. To configure the LEO satellite as a receiver, we can click on the Settings-Vehicle tab and specify the Keplerian Elements. For example, we can change reference time that will change satellite position and adjust the best line of sight (LOS) between the LEO satellite and the emitter; see Figure 6. Regarding antenna pattern, it is the same situation as the transmitter; we import the antenna pattern or select it None to avoid interfering with the signal’s power.

We can repeat the same steps in the next instances, and put the first instance as master and others as slave. From the Settings-Global-Synchronize Simulators section, we can select the instance if it is master or slave, see Figure 7. Eventually, we can save all changes at each instance, and click on start at instance-1 to run the scenario configured. All configurations are saved in the Skydel-SDX/Configurations folder, and the configuration files use the (sdx) filename. After passing the selected time, we can stop the running simulation and save the simulation scenario at I/Q files (CSV format) [44].

3. Mathematical Model of Geolocation and Tracking Scenarios

3.1. Conventional Geolocation Measurements

The signal transmitted from a high-speed dynamic emitter received at multi-sensors has time delay and frequency Doppler. Therefore, we can compute TDOA and FDOA between the emitter and the ith sensor [20].

Let the emitter position at time index k be,

$${\mathit{u}}_k^o={\mathit{u}}_0^o+\zeta k{\dot{\mathit{u}}}_k^o,$$

(1)

where ${\mathit{u}}_0^o={[x_0^o,y_0^o,z_0^o]}^{\prime }$ and ${\dot{\mathit{u}}}_k^o={[{\dot{x}}_k^o,{\dot{y}}_k^o,{\dot{z}}_k^o]}^{\prime }$ are initial location vector and velocity vector of emitter, respectively, $\zeta$ is time between successive measurements, ${(·)}^{\circ }$ represents the true value without noise, and ${[·]}^{\prime }$ denotes the matrix transpose.

3.1.1. TDOA

The TOA measurement at ith sensor is given by

$$\left\{\begin{array}{c}{\tau }_{i,k}^o=\frac{\parallel d_{i,k}^o\parallel }{c},\hfill \\ d_{i,k}^o=({\mathit{u}}_k^o-{\mathit{s}}_i^o),\{i=1,2,\dots ,M\},\hfill \end{array}\right.$$

(2)

where $\parallel ·\parallel$ represents the 2-norm, and c is the propagation speed ≈2.998 × 10⁸ m/s, $d_{i,k}^{\circ }$ is range between emitter and ith sensors in (m), as well ${\mathit{s}}_i^o={[x_i^o,y_i^o,z_i^o]}^{\prime }$ is initial location vector of ith sensor.

The time difference of signal arrival (TDOA) between ith sensors when we use sensor #1 as reference is calculated as

$${\tau }_{i,1,k}={\tau }_{i,k}^o-{\tau }_{1,k}^o+n_{i,1,k},\{i=2,3,\dots ,M;k=1,2,\dots ,K\},$$

(3)

where $n_{i,1,k}\sim N(0,{\sigma }_{\tau }^2)$. The vector form of TDOA can be represented as

$${\mathit{\tau }}_k={\mathit{\tau }}_k^o+{\mathbf{n}}_k,$$

(4)

where

$$\left\{\begin{array}{c}{\mathit{\tau }}_k={[{\tau }_{2,1,k},{\tau }_{3,1,k},\cdots {\tau }_{M,1,k}]}^{\prime },\hfill \\ \\ {\mathit{\tau }}_k^o={[{\tau }_{2,1,k}^o,{\tau }_{3,1,k}^o,\cdots {\tau }_{M,1,k}^o]}^{\prime },\hfill \\ \\ {\mathbf{n}}_k={[n_{2,1,k},n_{3,1,k},\cdots n_{M,1,k}]}^{\prime },\hfill \end{array}\right.$$

(5)

Covariance matrix $({\mathit{R}}_{n_k})$ of TDOA noise is determined as

$${\mathit{R}}_{n_k}=\mathrm{E}[{\mathbf{n}}_k{\mathbf{n}}_k^{\prime }]=\left[\begin{array}{c}\underset{(M-1)(K+1)}{\underbrace{{\sigma }_{\tau }^2\cdots {\sigma }_{\tau }^2}}\end{array}\right],$$

(6)

where $\mathrm{E}[·]$ represents the statistical expectation; and ${\sigma }_{\tau }^2$ is the variance of TDOA measurement noise.

$${\sigma }_{\tau }=\frac{1}{B_S}\frac{1}{\sqrt{B_{N}T}}\frac{1}{\sqrt{S/N}},$$

(7)

where $B_S$ is signal bandwidth, $B_N$ is noise bandwidth, T is integration time, and $S/N$ is interference power [20].

3.1.2. FDOA

The FOA measurement at the ith sensors is given by

$$\left\{\begin{array}{c}{\mathit{f}}_{i,k}^o=f_o-f_o\frac{{\dot{d}}_{i,k}^o}{c},\hfill \\ {\dot{d}}_{i,k}^o=\frac{{\left(d_{i,k}^o\right)}^{\prime }\left({\dot{\mathit{u}}}_k^o-{\dot{\mathit{s}}}_i^o\right)}{\parallel d_{i,k}^o\parallel },\hfill \end{array}\right.\{i=1,2,\dots ,M\},$$

(8)

where $f_o$ is assumed carrier frequency received by the sensor in (Hz), ${\dot{d}}_{i,k}^{\circ }$ is the range rate, also called the Doppler or radial velocity between emitter, and ith sensors in (m/s), and ${\dot{\mathit{s}}}_i^o={[{\dot{x}}_i^o,{\dot{y}}_i^o,{\dot{z}}_i^o]}^{\prime }$ is initial velocity vector of ith sensor (m/s).

The FDOA between ith sensors when use sensor #1 as reference is calculated as

$$f_{i,1,k}=f_{i,k}^o-f_{1,k}^o+{\dot{n}}_{i,1,k},\{i=2,3,\dots ,M;k=1,2,\dots ,K\},$$

(9)

where ${\dot{n}}_{i1,k}\sim N(0,{\sigma }_f^2)$. The vector form of FDOA can be represented as

$${\mathit{f}}_k={\mathit{f}}_k^o+{\dot{\mathbf{n}}}_k,$$

(10)

where

$$\left\{\begin{array}{c}{\mathit{f}}_k={[f_{2,1,k},f_{3,1,k},\cdots f_{M,1,k}]}^{\prime },\hfill \\ \\ {\mathit{f}}_k^o={[f_{2,1,k}^o,f_{3,1,k}^o,\cdots f_{M,1,k}^o]}^{\prime },\hfill \\ \\ {\dot{\mathbf{n}}}_k={[{\dot{n}}_{2,1,k},{\dot{n}}_{3,1,k},\cdots {\dot{n}}_{M,1,k}]}^{\prime },\hfill \end{array}\right.$$

(11)

Covariance matrix $({\mathit{R}}_{{\dot{n}}_k})$ of FDOA noise is determined as

$${\mathit{R}}_{{\dot{n}}_k}=\mathrm{E}[{\dot{\mathbf{n}}}_k{\dot{\mathbf{n}}}_k^{\prime }]=\left[\begin{array}{c}\underset{(M-1)(K+1)}{\underbrace{{\sigma }_f^2\cdots {\sigma }_f^2}}\end{array}\right],$$

(12)

where ${\sigma }_f^2$ is the variance of FDOA measurement noise [20].

$${\sigma }_f=\frac{\sqrt{3}}{\pi T}\frac{1}{\sqrt{B_{N}T}}\frac{1}{\sqrt{S/N}}.$$

(13)

3.1.3. Hybrid TDOA/FDOA

The TDOA and FDOA measurements can be arranged in one vector ${\mathbf{h}}_{i,1}({\mathbf{x}}_k)$ as

$${\mathbf{h}}_{i,1}({\mathbf{x}}_k)=\left(\begin{array}{c}{\mathit{\tau }}_k\\ {\mathit{f}}_k\end{array}\right)=\left(\begin{array}{c}{\tau }_{i,k}-{\tau }_{1,k}\\ f_{i,k}-f_{1,k}\end{array}\right),$$

(14)

$${\mathbf{H}}_k({\mathbf{x}}_k)={[{\mathbf{h}}_{2,1}({\mathbf{x}}_k),{\mathbf{h}}_{3,1}({\mathbf{x}}_k),\cdots {\mathbf{h}}_{M,1}({\mathbf{x}}_k)]}^{\prime },$$

(15)

The corresponding noise vector is

$${\widehat{\mathbf{n}}}_k={\left[\begin{array}{cc}{\mathbf{n}}_k^{\prime }& {\dot{\mathbf{n}}}_k^{\prime }\end{array}\right]}^{\prime },$$

(16)

which assumes zero-mean Gaussian with covariance matrix ${\mathit{R}}_k$ of hybrid (TDOA/FDOA) [45].

$$\left\{\begin{array}{c}\mathrm{E}[{\widehat{\mathbf{n}}}_k]=0,\hfill \\ \mathrm{E}[{\widehat{\mathbf{n}}}_k{\widehat{\mathbf{n}}}_k^{\prime }]={\mathit{R}}_k,\hfill \end{array}\right.$$

(17)

where

$${\mathit{R}}_k=\left[\begin{array}{cc}{\mathit{R}}_{n_k}& \\ & {\mathit{R}}_{{\dot{n}}_k}\end{array}\right]=\left[\begin{array}{cc}\underset{(M-1)(K+1)}{\underbrace{{\sigma }_{\tau }^2\cdots {\sigma }_{\tau }^2}}& \\ & \underset{(M-1)(K+1)}{\underbrace{{\sigma }_f^2\cdots {\sigma }_f^2}}\end{array}\right],$$

(18)

3.2. State Tracking Estimation

The state estimation is used as information in the conventional Kalman filter. The nonlinear system process ${\mathbf{x}}_{k+1}$ and observation process is specified by

$${\mathbf{x}}_k=\mathbf{f}({\mathbf{x}}_{k-1})+{\mathbf{v}}_{k-1},$$

(19)

$${\mathbf{y}}_k={\mathbf{H}}_k({\mathbf{x}}_k)+{\widehat{\mathbf{n}}}_k,$$

(20)

where $\mathbf{f}(·)$ is state transition matrix; ${\mathbf{H}}_k({\mathbf{x}}_k)$ is geolocation measurement matrix; ${\mathbf{v}}_{k-1}$ is system noise and ${\widehat{\mathbf{n}}}_k$ geolocation measurement noise, which are statistically independent Gaussian processes of zero mean and known covariance matrices (${\mathbf{Q}}_{k-1}$, ${\mathit{R}}_k$) [9,35].

From (19), the state–space model for a constant velocity of target-tracking can be defined as

$$\begin{array}{c}\hfill {\mathbf{x}}_k=\left[\begin{array}{ccccc}1& \frac{\sin ({\omega }_{k-1}\Delta t)}{{\omega }_{k-1}}& 0& \frac{\cos ({\omega }_{k-1}\Delta t)-1}{{\omega }_{k-1}}& 0\\ 0& \cos ({\omega }_{k-1}\Delta t)& 0& -\sin ({\omega }_{k-1}\Delta t)& 0\\ 0& \frac{1-\cos ({\omega }_{k-1}\Delta t)}{{\omega }_{k-1}}& 1& \frac{\sin ({\omega }_{k-1}\Delta t)}{{\omega }_{k-1}}& 0\\ 0& \sin ({\omega }_{k-1}\Delta t)& 0& \cos ({\omega }_{k-1}\Delta t)& 0\\ 0& 0& 0& 0& 1\end{array}\right]({\mathbf{x}}_{k-1})+{\mathbf{v}}_{k-1},\end{array}$$

(21)

and system covariance matrix is

$$\left\{\begin{array}{c}\mathrm{E}[{\mathbf{v}}_{k-1}]=0,\hfill \\ \mathrm{E}[{\mathbf{v}}_{k-1}{\mathbf{v}}_{k-1}^{\prime }]={\mathbf{Q}}_{k-1}=\left[\begin{array}{ccccc}\frac{\Delta t^3}{3}& \frac{\Delta t^2}{2}& 0& 0& 0\\ \frac{\Delta t^2}{2}& \Delta t& 0& 0& 0\\ 0& 0& \frac{\Delta t^3}{3}& \frac{\Delta t^2}{2}& 0\\ 0& 0& \frac{\Delta t^2}{2}& \Delta t& 0\\ 0& 0& 0& 0& 1.75\times {10}^{-3}\Delta t\end{array}\right],\hfill \end{array}\right.$$

(22)

where ${\mathbf{x}}_k={[{\mathit{u}}_k,{\dot{\mathit{u}}}_k,{\omega }_k]}^{\prime }$ is the target state that consists of the three-dimensional (3D) position vector (${\mathit{u}}_k^{\prime }$), velocity vector (${\dot{\mathit{u}}}_k^{\prime }$), and rotation rates (${\omega }_k$) at time instant $t_k$; $\Delta t$ is the sampling interval [46].

3.3. Geolocation and Target-Tracking Techniques

3.3.1. Cross Ambiguity Function (CAF) Approach

CAF combines TDOA and FDOA to achieve more accuracy and reduce the requirements for geolocation and mapping measurements. The measurements of TDOA and FDOA are used for cross-correlation equations, as presented in (23). CAF is applied in order to compute the peaks that are assumed to be the sources of the emitter using the CAF function. Then, the evaluation will be performed to do selection of the most accurate position, such as the higher peak, as a precise interference source.

$$\mathrm{CAF}(\mathit{\tau },\mathit{f})={\int }_0^{T}si{g}_1(t)si{g}_2^{*}(t-\mathit{\tau })e^{-j2\pi \mathit{f}t}dt,$$

(23)

where $si{g}_1$ and $si{g}_2$ are continuous-time signals; “$\ast$” means complex conjugate; T is the integration period in seconds; $\mathit{\tau }$ is the time delay in seconds, and $\mathit{f}$ is the frequency offset in Hertz. When the RFI emitter sends its signal towards the sensors, it reaches each sensor at a different time and with shifted frequency (due to their motion). The signals are captured for the selected period T. Binary-phase-shift-keying (BPSK) signals takes two-phase reversals such as 0° and 180° in this model. This is widely used in satellite communications.

A sinusoidal carrier wave is modulated by binary digits “1” and “0”, so the data signal shifts the phase of the carrier waveform is either 0° or 180°. From the trigonometric relationships,

$$\left\{\begin{array}{c}sin(x+\pi )=-sin(x),\hfill \\ cos(x+\pi )=-cos(x),\hfill \end{array}\right.$$

(24)

There are two possible states in a BPSK modulation, where the carrier is multiplied by ±1. Therefore, general analytic expression of BPSK signal can be given as,

$$si{g}_i(t)=Acos\left[2\pi f_{c}t+{\phi }_i(t)\right]\left\{\begin{array}{c}0\le t\le T_{sym},\hfill \\ i=1,2\hfill \end{array}\right.$$

(25)

where A is the amplitude, $f_c$ is the carrier frequency, ${\phi }_i(t)$ is the data symbol period, which takes on the values of “0” or “$\pi$”, and $T_{sym}$ is the data symbol period that equal $(1/R_{sym})$. $R_{sym}$ is a symbol rate for the generated BPSK signals.

By substituting (24) in (25), it can achieve two waveforms transmitted in a BPSK of continuous-time or analog signals

$$\left\{\begin{array}{c}si{g}_1(t)=Acos(2\pi f_{c}t),\hfill \\ si{g}_2(t)=-Acos(2\pi f_{c}t),\hfill \end{array}\right.$$

(26)

The discrete time is represented as

$$\left\{\begin{array}{c}si{g}_1(n)=Acos[2\pi f_c(n{T}_s)],\hfill \\ si{g}_2(n)=-Acos[2\pi f_c(n{T}_s)],\hfill \end{array}\right.$$

(27)

Practically, (27) must be in discrete form to be easy for digital signal analysis and computation. For that, sampling frequency $(f_s)$ must change to sample duration $T_s=1/f_s$, Therefore let

$$\left\{\begin{array}{c}t=n{T}_s,\hfill \\ \mathit{f}=\frac{{\mathit{f}}^D{f}_s}{N},\hfill \end{array}\right.$$

(28)

where ${\mathit{f}}^D$ is the digital FDOA, N is the number of samples for each signal that equal $\frac{T}{T_s}$.

From that, we can produce the discrete CAF.

$${\displaystyle \mathrm{CAF}(\mathit{\tau },{\mathit{f}}^D)=\sum _{n=0}^{N-1}si{g}_1[n]si{g}_2^{*}[n-\frac{\mathit{\tau }}{T_s}]e^{-j2\pi \frac{{\mathit{f}}^{D}n}{N}},}$$

(29)

By substituting (5) and (11) in (29), as well as letting $n_0$ equal the first sample of the receiving signals and $N=f_{s}T$, the CAF surface can modify as [47]

$$\begin{array}{c}\hfill {\displaystyle \mathrm{CAF}(x,y)=\sum _{n=n_0}^{n_0+f_{s}T-1}si{g}_1[n]si{g}_2^{*}[n-\mathrm{Round}(\frac{\mathit{\tau }}{T_s})]e^{-j2\pi \frac{f_{c}n}{f_s}[f_2(n{T}_s)-f_1(n{T}_s)]},}\end{array}$$

(30)

Because the term $\{\frac{\mathit{\tau }}{T_s}\}$ is in the digital domain and corresponds to an offset in the sample index, it must be rounded to the nearest integer. That means fractional offsets are not permitted. Computing the CAF surface in the $xy$ domain is dependent on the number of snapshots taken.

3.3.2. Gauss–Hermite Quadrature Filter (GHQF) Approach

Different versions of derivative-free filters were derived to employ EKF, Unscented Kalman Filter (UKF), and are generally called Sigma Point Kalman Filter (SPKFs), including Central Difference Kalman Filter (CDKF), Divided Difference Filter (DDF), and Particle Filters as well in the literature [41]. The Cubature Information Filter (CIF) was used in [36] to solve the pedestrian integrated navigation problem in a foot-mounted sensor fusion design. The Unscented Information Filter (UIF) was used for target-tracking in [48]. In [36], limitations of variant versions of 3rd- and 5th-degree CKF/CIF were observed in simulation with high initialization error when the superiority of the Gauss–Hermite and the 7th-degree CKF/CIF have been demonstrated and then justified. Based on that, an important question arises: how can one adapt the best cubature rule degree to the state estimation problem parameters (initialization, non-linearity, state dimension) in order to achieve the best estimate with less computational complexity. For that, GHQF 3rd and 5th degrees were considered and verified their performance.

Referring to the various derivations of multiple Gaussian-point filters, and through the following state-of-art of authors in [36], we can conclude a uniform algorithm of different square Kalman filters. The algorithm achieved will be a difference of others in how the points are, and the weights given in (31) are calculated [41]. The Gauss–Hermite rule can be created using a Gaussian weighted integral [49,50].

$${\displaystyle {\int }_{R_n}\mathbf{f}(\mathbf{x})\mathcal{N}(\mathbf{x};\widehat{\mathbf{x}},\Sigma )d\mathbf{x}\approx \sum _{i=1}^N{\mathrm{w}}_i\mathbf{f}(\widehat{\mathbf{x}}+\sqrt{\Sigma }{\xi }_i),}$$

(31)

where $\sqrt{\Sigma }$ is the Cholesky decomposition of $\Sigma$ and satisfies the relation $\Sigma =\left(\sqrt{\Sigma }\right){\left(\sqrt{\Sigma }\right)}^{\prime }$, ${\xi }_i$ is the quadrature point, and ${(·)}^{\prime }$ denotes the matrix transposition. Using the algorithm of higher-degree filtering, we can use square conversions, cubes, and convert various sigma points to create points and weights. The numerical integration methods of the quadrature integral given in (31) are suitable when the measurement noises with covariance matrices ${\mathbf{Q}}_{k-1}$ and ${\mathit{R}}_k$ of process state and localization, respectively. The considered algorithm introduces three steps of the Gauss–Hermite quadrature filtering.

• First step: Initialization
  In the first step,

  (i)

  Initialize the mean (${\widehat{\mathbf{x}}}_{k|k-1}$) and covariance (${\mathbf{P}}_{k-1|k-1}$) of the random variable $\mathbf{X}(0)\le ({\xi }_i,W_i)$ with ${\widehat{\mathbf{x}}}_{0|0}$, and ${\mathbf{P}}_{0|0}$.

  (ii)

  Compute the quadrature points by

  $${\xi }_i=\sqrt{2}{\epsilon }_i,$$

  (32)

  where ${\epsilon }_i$ is the ith eigenvalue of $\mathbf{J}$ that is supposed to be a symmetric tridiagonal matrix with zero diagonal elements given by

  $${\mathbf{J}}_{i,i+1}=\sqrt{i/2}=0,$$

  (33)

  (iii)

  Calculate the respective weights $W_i$ of the quadrature points (${\xi }_i$), where $W_i$ is equal to the square of the first element of the ith normalized eigenvector of $\mathbf{J}$.

• Second step: Time update
  In the second step, it is necessary to evaluate and estimate state

  $${\mathbf{X}}_{k-1|k-1}^j=\sqrt{{\mathbf{P}}_{k-1|k-1}}{\xi }^j+{\widehat{\mathbf{x}}}_{k-1|k-1},$$

  (34)

  $${\mathbf{X}}_{k|k-1}^{*(j)}=\mathbf{f}({\mathbf{X}}_{k-1|k-1}^j,{\mathit{u}}_{k-1},k-1),$$

  (35)

  $${\widehat{\mathbf{x}}}_{k|k-1}=\sum _{j=1}^m{\mathrm{w}}^j{\mathbf{X}}_{k|k-1}^{*(j)},$$

  (36)

  $$\begin{array}{c}\hfill {\mathbf{P}}_{k|k-1}=\sum _{j=1}^m{\mathrm{w}}^j\left({\mathbf{X}}_{k|k-1}^{*(j)}-{\widehat{\mathbf{x}}}_{k|k-1}\right)\times {\left({\mathbf{X}}_{k|k-1}^j-{\widehat{\mathbf{x}}}_{k|k-1}\right)}^{\prime }+{\mathbf{Q}}_{k-1},\end{array}$$

  (37)

  where ${\{{\mathbf{X}}_{k-1|k-1}^j\}}_{j=1}^m$ is the quadrature points; ${\{{\mathbf{X}}_{k|k-1}^{*(j)}\}}_{j=1}^m$ is the propagated quadrature points; ${\widehat{\mathbf{x}}}_{k|k-1}$ is the predicted state mean; and ${\mathbf{P}}_{k|k-1}$ is the predicted error covariance.
• Third step: Measurement update
  In the third step, it is necessary to evaluate and estimate measurement

  $${\mathbf{X}}_{k|k-1}^j=\sqrt{{\mathbf{P}}_{k|k-1}}{\xi }^j+{\widehat{\mathbf{x}}}_{k|k-1},$$

  (38)

  $${\mathbf{Y}}_{k|k-1}^j=\mathbf{h}({\mathbf{X}}_{k|k-1}^j,{\mathit{u}}_k,k),$$

  (39)

  $${\widehat{\mathbf{y}}}_{k|k-1}=\sum _{j=1}^m{\mathrm{w}}^j{\mathbf{Y}}_{k|k-1}^j,$$

  (40)

  $$\begin{array}{c}\hfill {\mathbf{P}}_{\mathbf{y}\mathbf{y},k|k-1}=\sum _{j=1}^m{\mathrm{w}}^j\left({\mathbf{Y}}_{k|k-1}^j-{\widehat{\mathbf{y}}}_{k|k-1}\right)\times {\left({\mathbf{Y}}_{k|k-1}^j-{\widehat{\mathbf{y}}}_{k|k-1}\right)}^{\prime }+{\mathit{R}}_k,\end{array}$$

  (41)

  $$\begin{array}{c}\hfill {\mathbf{P}}_{\mathbf{x}\mathbf{y},k|k-1}=\sum _{j=1}^m{\mathrm{w}}^j\left({\mathbf{X}}_{k|k-1}^j-{\widehat{\mathbf{x}}}_{k|k-1}\right)\times {\left({\mathbf{Y}}_{k|k-1}^j-{\widehat{\mathbf{y}}}_{k|k-1}\right)}^{\prime },\end{array}$$

  (42)

  $${\mathrm{W}}_k={\mathbf{P}}_{\mathbf{x}\mathbf{y},k|k-1}{\mathbf{P}}_{\mathbf{y}\mathbf{y},k|k-1}^{-1},$$

  (43)

  $${\widehat{\mathbf{x}}}_{k|k}={\widehat{\mathbf{x}}}_{k|k-1}+{\mathrm{W}}_k({\mathbf{y}}_{k|k}-{\widehat{\mathbf{y}}}_{k|k-1}),$$

  (44)

  $${\mathbf{P}}_{k|k}={\mathbf{P}}_{k|k-1}-{\mathrm{W}}_k{\mathbf{P}}_{k|k-1}{\mathrm{W}}_k^{\prime },$$

  (45)

  where ${\{{\mathbf{X}}_{k|k-1}^j\}}_{j=1}^m$ is the new quadrature points; ${\{{\mathbf{Y}}_{k|k-1}^j\}}_{j=1}^m$ is the new propagated quadrature points; ${\widehat{\mathbf{y}}}_{k|k-1}$ is the predicted measurement; ${\mathbf{P}}_{\mathbf{y}\mathbf{y},k|k-1}$ is the innovation covariance matrix; ${\mathbf{P}}_{\mathbf{x}\mathbf{y},k|k-1}$ is the cross-covariance matrix; ${\mathrm{W}}_k$ is the Kalman gain; ${\widehat{\mathbf{x}}}_{k|k}$ is the updated state; and ${\mathbf{P}}_{k|k}$ is the corresponding error covariance.
  Eventually, this can obtain the posterior density; therefore, a nonlinear Kalman filter estimates the hidden state using a Gaussian distribution of probability density functions.

  $$\mathbf{P}({\mathbf{x}}_k\mid {\mathbf{y}}_{1:k})=\mathcal{N}({\widehat{\mathbf{x}}}_{k\mid k},{\mathbf{P}}_{k|k}),$$

  (46)

  where ${\mathbf{y}}_{1:k}$ is the set of measurements ${\mathbf{y}}_{1:k}=[{\mathbf{y}}_1,{\mathbf{y}}_2,\cdots {\mathbf{y}}_K]$.

In this work, the visibility and availability of multiple LEO satellite constellations are assumed to be sequential, ensuring a minimum of one LEO satellite signal availability. A distributed sequential high-degree nonlinear filtering algorithm was developed and compared in the simulation with a distributed sequential filter such CKF 3rd degree, CIF 3rd degree, CIF-V1 5th degree, CIF-V2 5th degree, CIF-V3 5th degree, and CIF-MY SOVSKKH 7th degree.

3.4. Algorithm Performance Evaluation

To know the performance of the emitter geolocation algorithm, Cramer–Rao lower bound (CRLB) is the best validator. CRLB is the inverse of the Fisher Information Matrix (FIM) [51]. In addition, an unbiased estimator performance for the emitter position and velocity is the Posterior CRLB (PCRLB). The PCRLB was employed to verify the algorithm performance of optimizing geolocation and tracking estimation by determining the increase in RMSE. The deriving of PCRLB and algorithm performance of nonlinear measurements is presented in Appendix A [52]. RMSE of position, velocity, and rotation rate for target-tracking are defined as [53].

$$\left\{\begin{array}{c}{\mathrm{RMSE}}_k^{pos.}=\sqrt{\frac{1}{N}{\sum }_{n=1}^N{\parallel {\mathit{u}}_k^n-{\mathit{u}}_k^{o^n}\parallel }^2},\hfill \\ \\ {\mathrm{RMSE}}_k^{vel.}=\sqrt{\frac{1}{N}{\sum }_{n=1}^N{\parallel {\dot{\mathit{u}}}_k^n-{\dot{\mathit{u}}}_k^{o^n}\parallel }^2},\hfill \\ \\ {\mathrm{RMSE}}_k^{rot.}=\sqrt{\frac{1}{N}{\sum }_{n=1}^N{\left({\omega }_k^n-{\omega }_k^{o^n}\right)}^2},\hfill \end{array}\right.$$

(47)

where $({\mathit{u}}_k^n,{\dot{\mathit{u}}}_k^n,{\omega }_k^n)$ and $({\mathit{u}}_k^{o^n},{\dot{\mathit{u}}}_k^{o^n},{\omega }_k^{o^n})$ are the true and estimated target positions, velocities, and rotation rate at the N runs of the Monte Carlo iterative [54,55].

The metrics have been used to compare the performance of all proposed algorithms. For a fair comparison, 100 runs of independent Monte Carlo have been made. The performance is verified by the RMSE of the position, velocity, and rotation rate defined in (47). In the considered scenarios, the initial estimate ${\widehat{\mathbf{x}}}_{0|0}$ was created randomly from the Gaussian distribution $\mathcal{N}({\mathbf{x}}_{0|0},{\mathbf{P}}_{0|0}$) in each run.

4. Modeling Setup and Realistic Experimental of Simulation Scenarios

This section introduces measurement results to analyze and demonstrate the accuracy of geolocation and tracking for RFI emitter using extracted data from the Orolia Skydel Simulator. Figure 8 illustrates the experimental setup used in this study.

4.1. Modeling Scenario

Realistic simulation scenarios using the Orolia Simulator considered three LEO satellites to estimate and track a UAV flying in a circle at high speed and powerful rotation. The UAV emits an RFI signal on the central frequency 1575.42 MHz, corresponding to the three GPS receivers onboard the three LEO satellites. The experiment has considered a binary phase-shift keying (BPSK) modulation scheme. Three instants are created and summarized as

• Instance-1: Simulate UAV trajectory, generate RFI signal (common for three scenarios), and create earth-orbiting spacecraft (setup-LEO satellite # 1)
  • Simulation of UAV Circular Trajectory.
  • Create new radio assign # 0.
  • Vehicle- Keplerian orbital elements setup.
  • Save it as Master into the named file (sdx format).
• Instance-2: Creating earth-orbiting spacecraft (setup LEO satellite # 2).
  • Create a new radio and assign # 1.
  • Vehicle-Keplerian orbital elements setup.

    (a)

    Change the reference time.

    (b)

    Other Ephemeris Elements.

  • Save it as Slave into the named file (sdx format).
• Instance-3: Creating earth-orbiting spacecraft (setup-LEO satellite # 3).

  ∗

  Repeating the same steps of the Instance-2.

• Save and Record I/Q files (CSV format), as well export RAW logging files for reference. Figure 9 illustrates snapshots from the three Instances.

I/Q files will be used as input information at the realistic scenarios of geolocation and tracking approaches. Useful quasi-real data extracted from the Orolia Simulator is presented in Appendix B.

4.2. RFI Emitter Geolocation Using CAF

This section covers the accurate RFI emitter geolocation scenario using the CAF approach based on extracted data from the Orolia Skydel Simulator. In this scenario, we have considered two LEO satellites as sensors that receive emitted RFI signals from the motionless UAV as an emitter in this approach. The two sensors $({\mathit{s}}_1$ and ${\mathit{s}}_2)$ and emitter ${\mathit{u}}_k$ are located at the coordinate positions shown in the grid. The symbols $d_1$ and $d_2$ are the relative position vectors between each of the sensors and the emitter, while ${\dot{\mathit{s}}}_1,{\dot{\mathit{s}}}_2$ and ${\dot{\mathit{u}}}_k$ are the respective velocity vectors. dm is the required resolution to achieve the best CAF-Map image. In this derivation, the emitter is considered stationary, so there is no time-variant (k), and the emitter’s velocity equals zero. Figure 10 illustrates the geometry of the geolocation scenario.

For your information, this figure represents one instantaneous snapshot $(Nu{m}_{Snaps})$ of the system geometry. Since the sensors are moving at their respective velocities, the geometry changes with each passing instant of time. That is why the TDOAs and FDOAs in a system are time-varying in nature. Therefore, CAF can be computed based on the TDOA and FDOA measurements. Consequently, it can assume the number of snapshots within a fixed gap in the considered scenario [56]. The steps of emitter geolocation using CAF are detailed as

• Input parameters:$f_c;f_s$; $R_{sym}$; dm; $Nu{m}_{Snaps}$; $S/N$; and create XY grid
• Snapshots Loop for break up the signals into snapshots.
  for$i=1,2,\cdots ,Nu{m}_{Snaps}$do,
• Import I/Q files for UAV and LEO satellites trajectory, as well the delay time and Doppler shift of receiving signal.
• Calculation: Calculating a CAF surface based upon input signals.
  • Calculating relative position vectors for each sensor and the current map location,
  • Start going through each map location in XY grid by creating a loop for XY index
    for$ii=1,2,\cdots ,$ length (indexX) do,
    for$jj=1,2,\cdots ,$ length(indexY) do,
    • Calculate relative position vectors for each sensor and the current map location,
    • Calculate the TOA and FOA between emitter and sensors via (2) and (8),
    • Compute TDOA and FDOA from the relative position vectors via (4) and (10),
    • Compute the CAF value for the current MAP position via (29),
    end loop $jj$
    end loop $ii$
  • Saving CAF mapping result of current snapshot.
  end i, end Snapshots Loop.
• Output result: Save file of CAF mapping; and plotting CAF result with the peak point.

4.3. RFI Emitter-Tracking Using High-Degree Nonlinear Filters

This section covers the scenario of accurate RFI emitter-tracking using the high-degree nonlinear tracking filters approach based on extracted data from the Orolia Skydel Simulator. In this scenario, three LEO satellites are considered to receive RFI emitted a signal from a UAV under high speed and powerful rotation. As mentioned in the previous section, the UAV moves circularly at the initial equal position (1000, 1000) with an altitude of 1000 m and an absolute speed of 100 m/s. The velocity has adjusted vertically to the line from the circle center to the origin position. In addition, the initial position of each LEO satellite is ${\mathit{s}}_1=(798.50,4502.10,5008.70)\mathrm{km},{\mathit{s}}_2=(1236.50,-4547.90,4876.07)\mathrm{km}$, ${\mathit{s}}_3=(1668.80,-4572.70,4721.00)\mathrm{km}$, and their velocity is ${\mathit{v}}_1=(7007.70,-993.50,2020.60)\mathrm{m}/\mathrm{s}$, ${\mathit{v}}_2=(6926.40,-678.30,-2399.70)\mathrm{m}/\mathrm{s}$, ${\mathit{v}}_3=(6813.10,-360.00,-2767.90)\mathrm{m}/\mathrm{s}$. The initial estimate ${\widehat{\mathbf{x}}}_{0|0}$ values with different uncertainty level values were introduced gradually from the Gaussian distribution $\mathcal{N}({\widehat{\mathbf{x}}}_{0|0},{\mathbf{P}}_{0|0})$ in which the true initial value ${\mathbf{x}}_{0|0}$ is ${\mathbf{x}}_{0|0}=(1000\mathrm{m},100\mathrm{m}/\mathrm{s},1000\mathrm{m},$${0\mathrm{m}/\mathrm{s},-3\mathrm{deg}/\mathrm{s})}^{\prime }$ and ${\mathbf{P}}_{0|0}$ being the initial covariance: ${\mathbf{P}}_{0|0}=\mathrm{diag}(100{\mathrm{m}}^2,10{\mathrm{m}}^2/{\mathrm{s}}^2,$$100{\mathrm{m}}^2,10{\mathrm{m}}^2/{\mathrm{s}}^2,100{\mathrm{m}\mathrm{rad}}^2/{\mathrm{s}}^2)$.

The steps of emitter-tracking using GHQF are detailed as

• Input parameters: TDOA and FDOA measurements based on data from the Orolia simulator; ${\mathit{Q}}_{k-1}$; initial position and velocity for sensors and RFI emitter; number of Monte Carlo runs (N); scan times $(K)$; and number of quadrature points $(m)$.

  (a)

  First step:

  • Initialization →${\widehat{\mathbf{x}}}_{0|0}$, and ${\mathbf{P}}_{0|0}$
    for$i=1,2,\cdots ,m$do,
  • Calculating the quadrature points via (31),
  • Calculating first element of the respective weights,
    end loop i,

• Start tracking: Time update and Measurement update.
  for$n=1,2,\cdots ,N$do,
  for$k=1,2,\cdots ,K$do,

  (b)

  Second step:

  • for$j=1,2,\cdots ,m$do,

  • Computing of evaluate state via (33) and (34),
  • Computing of estimate state via (35) and (36),

  • end loop j,

  (c)

  Third step: update.

  • for$j=1,2,\cdots ,m$do,

  • Computing of evaluate state via (37) and (38),
  • Computing of estimate state via (39) to (44),

  • end loop j,

  end loop k,
  end loop n,
• Evaluation: Evaluate the posterior density via (45).
• Output result: The accurate estimates of the emitter position and velocity.

5. Experimental Results

5.1. RFI Emitter Geolocation Scenario

In this scenario, the CAF method was applied to achieve peak of RFI emitter estimation and cross-correlation of TDOA and FDOA measurements based on extracted data from the Orolia Simulator. Figure 11 is illustrated the RMSEs of the position for the RFI emitter using measurements of TDOA and FDOA individually and hybrid, as well as, Figure 12 is illustrated the peak CAF and cross-correlation of TDOA and FDOA, which denoted the estimation position of the RFI emitter.

5.2. RFI Emitter-Tracking Scenario

In this scenario, the high-degree nonlinear tracking filters were applied to estimate the position, velocity, and turn rate of RFI emitter based on extracted data from the Orolia Simulator. Figure 13 shows the emitter-tracking using nonlinear filtering, where Figure 14 presents the zoom of the trajectory and filter tracking. Figure 15, Figure 16 and Figure 17 show the RMSEs of the position, velocity, and turn rate for the RFI emitter using the proposed filters.

• Important notes

It is important to mention that the RFI target-tracking problem in this work is different from other standard problems known in the literature because of the nature of the measurement. TDOA/FDOA measurement is highly nonlinear and different in range and bearing, and is very well analyzed and understood from the RADAR community. Indeed, as nonlinear degree becomes different, observability degree and analysis become different. It is traditionally important to develop particle-filtering methods for target-tracking problems; however, in this paper, we have restricted our analysis to nonlinear derivative-free algorithms such as CKF GHKF of different degrees, as these algorithms present more advantages compared to particle filters and sequential particle filters when speaking about real-time implementation feasibility.

5.3. Discussion and Performance Analysis

From Figure 11, it can be noted that using TDOA to geolocate the UAV as an RFI emitter is problematic since the TDOA is based on the difference in arrival time. In contrast, FDOA estimation is based on frequency Doppler contrast. In this way, it is possible to observe how the LEO satellite high-accelerated linear velocity affects the rate of FDOA calculations when using FDOA alone or in combination with TDOA/FDOA, which is crucial for identifying and following the UAV. As a result, the RMSE convergence speed when using FDOA individually, or a combination of TDOA/FDOA is faster than when using only TDOA in identical situations.

Figure 12 shows the particulars of the peak and cross-correlation of TDOA and FDOA in the CAF-Map. The distance deviation, in this case, was only 650 m, and there were no issues with left–right ambiguity. The CAF-resolution maps were set to 500 m, so the results were extremely encouraging. The implications of slight adjustments in the geometry of the sensing platforms were demonstrated in this scenario.

From Figure 13, Figure 14 and Figure 15, it can be noted that the robust UAV (as an RFI emitter) tracking occurred when applying high-degree nonlinear tracking filtering. From Figure 14, we can note a high diversion of 3rd-degree CKF tracking because of the sensitivity of CKF.

From Figure 16, Figure 17, Figure 18 and Figure 19 the RMSE convergence speed when using FDOA individually or a combination of TDOA/FDOA is faster than when using only TDOA in identical situations. As a result, we can note changes of a high-degree nonlinear filtering tracking based on TDOA/FDOA, where the 3rd-degree CKF displays much larger errors than the 3rd-degree GHQF or 5th-degree GHQF 3rd-degree GHQF maintains a good performance, but it is less accurate than the 5th-degree GHQF, which has the best performance and an indistinguishable accuracy. In Figure 17b and Figure 18c, we can note the evidence that CKF is very sensitive to initialization values compared to GHQF. The GHQF is more accurate and has faster convergence to the real state. CKF requires more time to converge and propagate higher error during the first 5 s. Some solutions exist using adaptive forms of CKF or what we call stochastic CKF sampling, in order to perform better initialization robustness. However, that is another problem to solve and will be considered in future works.

Higher degrees of CKF and higher degrees of GHKF were implemented and it is shown that GHKF 5th degree presented the best results in the three state estimation problems even with loss of observability of the third state.

Table 2. Average RMSE for a position (Pos.), velocity (Vel.), and rotation rate ($\omega$) implemented by experimental and theoretical simulation.

	RMSEs of Algorithm Validation with the Experimental Task			RMSEs of Algorithm Validation with the Theoretical Simulation
	Pos. (m)	Vel. (m/s)	$\mathit{\omega }$ (deg/s)	Pos. (m)	Vel. (m/s)	$\mathit{\omega }$ (deg/s)
3rd Deg CKF	72.2575	43.8975	5.7510	69.1524	40.8675	4.2712
3rd Deg GHQF	60.8740	38.8890	5.2535	58.9539	36.9493	4.2130
5th Deg GHQF	60.0550	38.2520	5.1500	58.6646	36.5513	4.0128

presented the performance validation for the proposed approach implemented by extracted data from the Orolia Simulator and compared it to the theoretical simulation.

Table 3. Comparative of quadrature points and execution time during implementation of the experimental task and theoretical simulations.

	Experimental Task		Theoretical Simulation
	Quadrature Points	Execution Time (ms)	Quadrature Points	Execution Time (ms)
3rd Deg CKF	10	0.95	10	0.8
3rd Deg GHQF	243	17.2	243	16.1
5th Deg GHQF	3125	209.0	3125	206.0

presents quadrature points and average execution time for 100 runs of the Monte Carlo iteration at the experimental task and theoretical simulation of the proposed filters. The main configurations of the computer used in this process are listed as the following: Intel (R) Core (TM) i7-6500 CPU @ 2.59 GHz; 8.00 GB RAM. According to the experimental task and theoretical simulation, the execution time of the filters is approximately proportional to the quadrature points.

It is clear and easy to note that the GHQF approach based on Orolia modeling and TDOA/FDOA measurements for tracking an RFI emitter at high speed and a powerful rotation has achieved perfect tracking compared to conventional filters. It is crucial to compare the performance of the modeling and proposed approaches with the latest and recently developed methods [3]. For that, we have created various scenarios for muti-LEO satellites and a UAV as emitter using the Orolia Skydel simulator and implemented at CAF and GHKF 3rd and 5th degrees based on TDOA/FDOA. In our approach, the results have been presented in terms of position, speed, and rotational force (Omega). Although [3] created their scenario for one Geosynchronous Equatorial Orbit (GEO) satellite with three antennas and stationary Earth station as emitter using STK, they implemented their scenario to localize the emitter position employing the RSS method.

6. Conclusions and Future Works

In this paper, we have investigated the accuracy of geolocation and tracking for RFI emitters at high speed and powerful rotation using extracted data from the Orolia Skydel simulator as a quasi-real environment. Employing the Orolia Skydel simulator, various scenarios for quasi-real dynamic trajectories of LEO satellites as sensors and a UAV as an RFI emitter have been constructed. The CAF and GHQF approaches were applied based on TDOA and FDOA measurements and data acquired from Orolia to geolocate and track a UAV at high speed and strong rotation circumstances. CRLB and PCRLB have verified the proposed approaches to determine the increase in RMSE using Monte Carlo iteration. Finally, it has been concluded that the proposed approaches can work perfectly with the extracted data from quasi-real environments. In future work, we propose to implement more experiments for complex algorithms such as considering the problem of initialization of CKF filter tracking by using stochastic CKF sampling, and considering an error estimated for position and velocity of RFI emitters achieved by the $H_{\infty }$/GHQF filtering with different values of $\gamma$ (performance bound) and we will work to improve emitter-tracking estimation with developed filtering in circumstances of sensor position uncertainty.