Precise-Point-Positioning Estimations for Recreational Drones Using Optimized Cubature-Extended Kalman Filtering

In modern-day multi-dimensional recreational drones (UAVs), the global navigation satellite system <italic>(GNSS)</italic> units in- use are commonly fraught with precise-point-positioning <italic>(PPP)</italic> data errors or inaccuracies, hence, necessitating this work. These data inaccuracies, occasioned by the system’s drawbacks such as sudden <italic>GPS</italic> lock or jamming, embedded device misalignment, drone’s limited communication coverage, signaling and detection, all contributes to the system’s <italic>PPP</italic> computation complexity. To mitigate PPP complexity, an intelligent and robust accurate continuous-discrete <italic>(ACD)</italic> based hybrid cubature-extended Kalman filter <italic>(C-EKF)</italic> computation model for an integrated <italic>GNSS</italic> unit is corroborated in this article. More precisely, time updates to the state and parameter sub-vectors for the <italic>GNSS</italic> unit is accomplished using the third-degree spherical-radial cubature rule. The system’s testbed simulation is then conducted using tightly-coupled units of (i) ring laser gyroscope (RLG) and (ii) micro-electro-mechanical system (MEMS) variants of the inertial measurement unit (IMU) to ascertain the PPP cooperative tendencies. Optimized performance comparisons of the proposed hybrid <italic>C-EKF</italic> over the existing cubature Kalman filter <italic>(CKF)</italic> and extended Kalman filter <italic>(EKF)</italic> models with-respect-to (w.r.t) its probabilistic outages, <italic>Yaw</italic> error-differences and ergodic capacities are demonstrated and presented.


I. INTRODUCTION
The frequent use of the inertial measurement unit (IMU) for industrial and space-based vehicular estimations have gained traction [1]- [6]. The IMU hardware, also called as the inertial navigation system (INS), are usually embedded in tightlycoupled GNSS devices such as the GPS unit. The Kalman filtering (KF) derivative is a direct differential approach that requires gradient-based methodologies in estimating position accuracies using the corresponding filter sensitivity computations. With KF, the GPS-IMU thus assigns the carrier phase as its observation phase. Additional phases connecting the state-vectors are further included thereby improving the The associate editor coordinating the review of this manuscript and approving it for publication was Min Wang . unit's PPP to its crucial sub-meter level. However, inherent nonlinear drawbacks downplay the GPS-IMU performance. Areas affected are in outage probability ratios and device misalignment scenario [4] leading to miscalculated cross-track errors. The error may seem little but can collapse systemcritical industries such as the avionics, pervasive systems and AI-powered systems [5], [6].
Recent introduction of up-scaled versions of CKF and EKF [7]- [10] for drone and radar-tracking systems have been recorded. For instance, [7], [8] focused at predicting the state mean and covariance matrix with measurement updates of radar tracking limitations showing when an aircraft executes a coordinated turn. They however failed to progress further with the integration of their scheme into GPS-IMU devices in order to ascertain its reliability. In [9], [10], attempt at VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ reviving the original EKF scheme with enhanced derivative computations was made. Modification was made to enhance time updates and measurement constraints for robust radartracking capacities, hence minimizing system's along-track errors. But no sufficient detail at mitigating existing CKF and EKF computational complexities were provided. The need to achieve unified GPS-IMU synchronization and error analysis with proportionate deployment resolution is important [2]. At present, there are three suitable GPS-IMU techniques for hoisting this work. They are the loosely-coupled, the tightly-coupled integrated and the ultratightly coupled model techniques [2]. The loosely-coupled and the tightly-coupled GPS-IMU models exhibits nonlinear constraints with inconsistent system errors as prevalent in all nonlinear systems. These are usually resolved using the EKF [11]. However, prior to implementing EKF, the system's position and velocity computations must attain certain PPP accuracy threshold for loosely-coupled, unlike in tightlycoupled. The research limitations on GPS-based [2] timeupdates at sub-vector [7]- [9] levels necessitates our work.
The main contributions in this paper as follows: • To formulate a hybrid and scalable cubature-extended Kalman filter (C-EKF) model by concurrently estimating the state and the parameter sub-vectors with recorded thresholds using the third-degree spherical-radial cubature rule. This is a further improvement to the works of [7]- [9] where only the differential features of EKF computation was investigated and introduced.
• The proposed C-EKF model is then simulated separately into each one of a tightly-coupled i) ring laser gyroscope (RLG) and (ii) micro-electro-mechanical system (MEMS) IMU devices to ascertain its PPP's cooperative tendencies. In [7], [8], tightly-coupled IMU devices were not used for tracking and detection, making them susceptible to inaccurate measurements.
• By deploying an integrated C-EKF into a GPS-IMU , further significant system accuracies such as the drone's motion dynamics, along-track and cross-track errors were analyzed and recorded. Unlike [2], our testbed was carried-out in a confined zone where frequent GPS-lock losses due to long GPS baselines are common occurrence, and factored in the cause of estimations.
• We then investigated the system's probabilistic capacity outages, Yaw error-differences and ergodic capacities in dissimilar situations of PPP caused by GNSS distortions and signal-to-noise ratio (SNR) thresholds. The proposed model is able to cross-interact with other orientation and position-estimation devices, such as the Light Detection And Ranging (LiDAR) and the Interferometric Synthetic Aperture Radar (InSAR), where direct sensor orientation measurements are applied on the exterior orientation parameters (three orientation angles and an added positioncoordinates of embedded navigation multi-sensor) of the system. The imaging sensors resolution is proportional to the hybrid model's navigation accuracy. Effective extraction of useful numerical information about the system's states from available measurements is successfully conducted, hence guaranteeing the unit's PPP.
The rest of this article is structured as follows: Section II discusses the problem formulation with a presentation of the current GPS-IMU integrated system and a review of time update mechanism on existing state-space filters. In Section III the proposed hybrid C-EKF model was derived and articulated. Section IV proposes novel random trajectory and guidance mechanism for the drones. The paper's performance evaluation with focuses on the testbed deployment, computations and numerical results were exhausted in Section V. In Section VI, broad look at this work's open research issues was investigated, while Section VII concludes the work with an expansive future works.

II. GPS-IMU SYSTEM AND EXISTING FILTERS: A REVIEW A. THE INTEGRATED GPS-IMU SYSTEM
The need to achieve unified GPS error analysis with proportionate resolution necessitated the GPS-IMU [2] deployment. The model is advanced and beneficial. The three (3) known techniques for achieving an integrated GPS + IMU are summarized as: the loosely-coupled, the tightly-coupled and the ultra-tightly coupled approaches.

• The Loosely-coupled Technique
This is when raw data (codes, sensed acceleration and angular rates for the inertial sensors, and carrier-phase pseudo ranges for GPS sensors) are processed to obtain the PPP's velocity and orientation estimations. The data is collated and pre-processed prior to being merged.
To obtain a merged data, KF are applied. This approach is common as it require at-least four (4) GPS satellites with the embedded sensors for data computation.
• The Tightly-coupled Technique The raw data from both GPS and inertial sensors are simultaneously processed to obtain desired goals. Obtained goals are more accurate estimations of velocity, orientation with the sensor's inertial error, located in the single filter process. It is more advanced since partial GPS raw data are available for pre-processing [2].
• The Ultra-tightly coupled Technique It is more complex approach of attempting to access GPS hardware information loop, in estimating the raw data (codes, sensed acceleration and angular rates for inertial sensors, and carrier-phase pseudo ranges for GPS sensors) for a more dynamic (or insecure) location. The general public cannot gain access to the GPS hardware information, making the approach impracticable. The loosely-coupled and the tightly-coupled integrated GPS-IMU models exhibit nonlinear constraints with inconsistent system errors as prevalent in all nonlinear systems. These are usually resolved with the EKF methodology [11]. While the system's position and velocity computations are required to achieve accurate PPP in the looselycoupled, they are however not required in the tightly-coupled approach.

B. THE EXTENDED KALMAN FILTER (EKF) SCHEME
The EKF scheme is deployed to sub-optimally provide nonlinear state estimations for nonlinear dynamic systems by using linearized Taylor approximation technique. For instance, it is used to resolve position accuracies in the integrated GPS-IMU units. However, the evolving dronebased GPS-IMU with its added surrounding environmental dynamism requires enhanced EKF capacity. To enhance its capacity, its non-linearity continuous-discrete feature from the Taylor linearized models are obtained as follows: where subscripts denoting the epochs; x,v,z and n are state sub-vector, process description noise vector, measurement vector and the measurement noise vector respectively. xk represents the state vector at epoch k,x − x, v and n are mutually independent of themselves. F and H denotes the nonlinear system and the measurement functions, process and measurement noises are assumed to possess a zero-mean value, x k represents the optimal state estimates at epoch k while thê x − k is the predicted state estimate at epoch k, and G are the matrices forming the linearized system model, while those of F and H are the matrices forming the linearized measurement model. The time update of a typical continuous-discrete filter is split into major computing groups. The group-I solves the Fokker-Plank's theorem by computing its conditional density using available numerical methods [12]- [14]. The group-II merely computes finite values of the conditional moments after discretizing the process equation using either higher order Runge-Kutta methods or the Eular approach [15], [16]. The underlying computational complexity constraints of the group-I filters increases with-respect-to (w.r.t) its state vector dimension, while that of group-II filters varies polynomially [17]. Similar to the group-II filter family, our proposed hybrid C + CKF model are computationally less-complex since only finite values are computed implying that its time updates are similar to those of continuous-discrete extended Kalman filter (CD-EKF) [11] and the continuous-discrete unscented Kalman filter (CD-UKF) [18].

• (CD)-Extended Kalman Filter:
The first known nonlinear derivation of CD-EKF in [11], [chapter 6 -9] was obtained from the EKF fundamentals, where its predicted state estimate and its covariance collectively form the time-update. The state estimate is: whereX k|k are expanded Taylor series derivatives [11]. The Taylor series is also harnessed to obtain the state error covariance matrix as: • (CD)-Unscented Kalman Filter: Similar to the CD-EKF scheme, the CD-UKF scheme obtains its predicted state-error and covariance matrix values from [11] as follows: where var[χ(t + δ)] remains Taylor derivatives.

C. THE CUBATURE KALMAN FILTER (CKF) SCHEME
The recently introduced (CKF) is a nonlinear discrete-time based state-space model that specifically targets the thirddegree spherical-radial cubature rule for computing all Gaussian-weighted integrals. The rule leads to an even number of equal cubature points 2n, (where n represents the state vector's dimensionality). The points are even and are distributed at the center of the ellipsoid. Some of the cubature rule's desirable properties are listed as follows: • The cubature rule is non-derivative, allowing it to seamlessly apply to areas where Jacobian and Hessians rules are commonly computed.
• The rule require 2n cubature points, where n is the state vector dimensions, making 2n to be imperatively evaluated at each update cycle. A third-degree cubature rules has a theoretical lower bound of 2n cubature points. These properties further justify why the CKF is new and optimal numerical scheme for nonlinear filtering.

5)
Estimate the predicted error covariance as:

D. BENEFITS OF HYBRID C-EKF DEPLOYMENT
The importance and benefits of deploying the proposed hybrid C-EKF model over the more conventional models is discussed. Key parameters and factors such as the model's interoperability with the introduced GPS-IMU devices, the signal's beamforming models, the relative computational complexity and the model's deployment purposes are all considered. As obtained in Table 1, two EKF-based schemes of square root (SR) and singular value decomposition (SVD), one upper diagonal (UD) factorization CKF-based scheme, and the original differential KF model are all compared using the proposed hybrid C-EKF approach as case model. The major drawbacks for all but the SVD-based EKF and the proposed scheme is in the deployment purposes, where the schemes can be implemented at the relaying (control towers), source (control center) and destination (drones) segments of the system model. Other key factor to observe are the system's interoperability with other navigation devices. While other schemes can either conduct one of a transmit or a receive operation at a time, the proposed C-EKF model simultaneously conducts both operations at the same time, making it robust. Aside from the CKF scheme which recorded a higher computational complexity, others schemes recorded a low or moderate complexities which are based on the conducted computations per minute.
To summarize, the previously differentiated-KF, the differentiated-SVD, and the differentiated-EKF based computations all provides quality estimations for position accuracy models by harnessing the and solving the gradientbased parameter estimation constraints, using adaptive filter techniques. The differentiated-CKF on the other attempts to harness and solve a combination of factor constraints such as time and vehicle trajectory. To make the system more robust, intelligent-driven and precisely accurate, combination of more advanced sets of the accurate continuous discreet (ACD) based CKF with EKF technique was deployed.

III. PROPOSED HYBRID C-EKF MODEL
The model's operational flowchart as captured in Fig. 1 presents the methodology for achieving desired drone accuracy. First, the drones are deployed, while trajectory and directions are navigated and controlled using embedded GPS-IMU devices. As previously stated in section 1, the IMU devices are either the tightly-coupled ring laser gyroscope (RLG) or the micro-electro-mechanical system (MEMS) units. A computation of the GPS-IMU signal beamforming with the proposed hybrid scheme is then deployed. On the deployment, PPP accuracy for the drones is achieved. To ascertain/verify these accuracies, goal posts are set for the drones to cross. Once they're able to cross the goal post line, target PPP accuracy is met. However, in the event goals are not met, the drone's GPS-IMU with C-EKF computations are re-calibrated by increasing the cubature node points by +1. The updated value is then sent to be re-computed by the GPS-IMU with C-EKF decision block.
The system architecture in Fig. 2 depicts the concurrent estimation of state-space and parameter sub-vectors predictions. The state's previous (priori) and aftermath (posteriori) estimations, with the covariance matrix of CKF's state, are all fed back-and-forth into EKF loop [1], [2], [6] as the system commences its execution. Estimation of the mixed state/parameter mobile drone sink is achieved with its linear observation as: where x k is state function with covariance matrix P 1,k , while f ( .) and θ k are the parametrized covariance matrix P 3,k . Also note that w and v are the Gaussian white noise, while the covariance matrices of Q and R refers to unity transition matrix. Taylor's 1 st -order term for non-linear function f are used for the cross covariance matrix (P 2,k ) predictions.

A. TO ESTIMATE STATE SUB-VECTOR
To guarantee optimal PPP using drones, state and parameter sub-vector positions for deployed drone are analyzed. Calculation of state sub-vector (a non-linear system part of the drone), embedded in the tightly-fixed GPS-IMU is achieved first. To update CKF time-line, the 2n points known as the cubature points are then deduced as: where n is the non-linear state vector's dimensionality. SVD is a single value method for matrix decomposition, while X are the cubature points. Note that, subsets ξ = m 2 (1) i , m = 2n, where n equals state vector's dimensionality and (1) i is system's generator. However, when n = 2, sets of points to denote the system generator (1) i are presented as: Applying non-linear state transition function in order to implement cubature points, latest predicted cubature becomes; To calculate both the predicted state and its related or relevant covariance matrix, mathematical models are expressed thus; where both P 1,k|k−1 and x k|k−1 remains the predicted covariance matrix and the predicted state matrix respectively. Since unwanted noise plays a pivotal role at not achieving reliable systems during cross-correlation with other drones and which occurs in the state transition matrix during estimation, enlargement of the non-linear state's noise covariance VOLUME 9, 2021 matrix (Q x ) system will provide some level of cushion in the overall system noise reduction and ultimately guarantee a reliable system. Eqn. (10) is only for generic systems devoid of tightly-coupled GPS-IMU unit. For our proposed tightly-coupled GPS-IMU enabled drones, further upgrade on eqn. (10) is made as: However, to update the measurement part, cubature points is required to be continuously regenerated using the predicted z k|k−1 and P 1,k|k−1 as follows: (17) To now calculate the predicted observation using newly defined cubature points, the following steps are introduced: In an assumed parameter, then states are only estimated. To have these measurement methodology revised by also estimating units for parameters, it is expressed as: On successfully modifying the measurement approach, a noise occurs in the form of estimation error θ emanating from θ k|k−1 . To convert noise into a good one, otherwise termed as 'white noise', the following expressions are introduced as: The expression implements further changes to the measurement approach in (11) as: The latest measurement approach in (13) is essential in the proposed model as it influences the system's estimation accuracy for better when changes to observation occurs. The state's white noise in (13) is accurately vilified as it possesses similar expressions with other existing (non-hybrid) approaches, when their respective cross-correlation (P 2 ) is zero. Hence, the need for state covariance matrix (P 1 ). Hence, picking our approach from the CKF, the state's Kalman gain is now expressed as: To finally express estimations for both a-posteriori state and covariance matrix, these follow-up steps are introduced:

B. TO ESTIMATE PARAMETER SUB-VECTOR
To guarantee collaboration and consistency of operations between EKF and the already finalized CKF, the parameter sub-vector of EKF must be analyzed with a pre-existing assumption of a known state's sub-vector. As predicted by EKF, parameters and covariance matrix are generated as: As a follow-up on the analyzed CKF, the generated outputs shown as predicted states and covariance matrix, acts as inputs for the system, together with state vector x k|k−1 and covariance vector P 3,k|k−1 . However, since the states' estimation errors are hindrances to prediction accuracy of parameters and covariance matrix, then a revision of the measurement equation is required and expressed thus: while its covariance matrix measurement is deduced as: Furthermore, expressions for Kalman gain, a-posteriori parameters and covariance measurement are updated as:

C. TO ESTIMATE CROSS-COVARIANCE MATRIX
The predictions and estimations for the cross-covariance matrix is achieved in this section. To begin, a computation of the CKF for the cross covariance is initialized as: Resulting EKF is implemented in our proposed model to further optimize the filtering method's reliability, often  degraded due to additional computational load. Mathematical step to obtain the predictions for cross covariance matrix is: Finally, a computation formula for a-posteriori cross covariance matrix is presented as: The features of (22) and (23) show resemblance to parameter covariance matrix. In generic systems, when cross covariance matrix (P 2 ) is computed, one major performancelimiting factor is parameter covariance matrix (P 3 ), which in-turn affects the estimation and outputs of Kalman gain. Extension of these limitation is observed in the proposed tightly-fixed GPS-IMU -enabled drones, where the negligibility of cross covariance matrix causes smaller Kalman gain thereby, hence limiting the convergence speed. However, in our approach, the flexibility of adopting different state's estimation for the two introduced sub-vectors meant is capable at multi-updating, guaranteeing estimation accuracy and further ensuring that computational loads are kept at barest minimum.

IV. TRAJECTORY AND DEPLOYMENT OF 4-D DRONES A. TESTBED DEPLOYMENT AND INVESTIGATIONS
This work was deployed into a region-of-interest of 6 km square-size (3500 m above sea level), mirroring a typically hazardous/conflict terrain of limited LOS and more signal interferences (e.g., noise, poor visibility, mountainous and in foggy/wintery condition). Detailed numerical configurations of EKF, CKF and hybrid C-EKF was compiled with the GPS-IMU processor codes to analyze variance in performance. Novatel GPS receiver and MEMS based IMU, called as imar-VRU, was then embedded into the drones. Data samplings from the GPS was kept at 1 Hz, while that of IMU was kept at 100 Hz. The lever-arm top of the drone, closer to the IMU was kept at nearly 1 m 2 apart. While the imar-RQH (a version of imar-VRU ) is deployed to observe the altitude reference and phase differential-GPS (DGPS). Inbuilt inertial DGPS features was used in obtaining results from imar-VRU /imar-RQH devices.
The integrated GPS-IMU based imar-VRU is depicted as relay device for estimating base stations to target coordinates orientations. Estimation accuracy hinges upon scanerror angle, footprints sizes, scan-range angle, imar-VRU range accuracy measurements, sensor-to-sensor proximity and the system's navigation trajectories (both vertical and horizontal). The scheme is implemented to retrieve the 4-D coordinates of the deployed aircraft, where base station is able to coordinate the parameters and actions of imar-VRU and have them transmitted to ground imar-VRU for analysis. The coordinates for geo-referencing is formulated as: where r M ,k is 3-D vector coordinate with k mapping. The drone's spatial resolution as determined by imar-VRU calibration tool, justifies higher target accuracy. The high accuracy is determined by the estimations of the GPS-IMU . By using Monte Carlo simulation techniques, sensitive orientation estimations to base station was demonstrated in Fig. 3. The simulation was done with 200 particle parameters to demonstrate both horizontal and vertical accuracies of GPS-IMU on the base station. The features were flight height (1000 m) to ground station, translational (±1 cm) and rotational (±10 arcsec) sensor misalignment accuracies, the imar-VRU range (±1 cm) and scan angle measurements (±5 arcsec). Others are the horizontal features (roll and pitch), the navigation selection methodology and mapping coordinates. To achieve higher accuracy readings, drone outputs with GPS-IMU is required. Where unit of centimeter-to-decimeter is replaced by centimeter-to-arcsec (roughly 0.0028 deg). Since the units are currently available in the more expensive technologies (eg., IMU devices), achieving maximal PPP accuracy for 4-D drone becomes imperative.

B. RANDOM TRAJECTORY MODEL
This section proposes random trajectory model of Fig. 4   GPS-IMU antennas are embedded to at-least two or more drones under consideration. During deployment, antennas may encounter obstacles [21] in the LOS, hindering optimal data transmission. To estimate the recalled drone's capability at data transmission with newly deployed ones, received signal analysis of the recalled drone must be deduced: where P U denotes drone's transmit power, x(t) is the transmit symbol, and n(t) is the drone's zero-mean additive white Gaussian noise (AWGN), consisting E[|n(t)| 2 ] = σ 2 . Complex coefficient between drones is h 1 , while LOS links are presumed to be contained as the atmospheric open space inside the spherical radius. The Rician factor (K ), normally accounting for the influencing reflections and scattered signals within the confining environment is corroborated, where the environment was assumed to have experienced Rician fading. Simplifying K in Rician fading, which is the ratio of LOS signal-power to power of scattered beam path, an assumption of the receiver-side CSI is made, leading to steps for probability density function (PDF) presented by Beard and McLain [22] is deduced as follows: whereθ = P U σ 2 l α , representing average SNR where path loss exponent is α.

1) INCLINATION AND DECLINATION STRATEGY
All drones designed for surveillance or combat purposes (e.g., unmanned combat aerial vehicles (UCAVs)), possess dissimilar pre-designed mechanism for inclination (climbing) and declination (descending). But their operational strategies are relatively similar in that their flight path angel (γ ) is always positive during take-off or when in-motion, and gives  negative value when about to descend. Fig. 5 gives the factored drone parameters as: overall weight of the craft (W), drag constant (D), craft's lift (L), and craft's thrust (T), which is proportional to its center of gravitational force. The craft's path-angle is obtained assuming constant velocity (V) and further re-construction of craft's repelling forces (x-axis) as: where inclination angle is sin(γ ). Further steps to obtain path angle for incline is: where m and g are flight angle variables by base station. Path angle (γ ) for declination using Y-axis is: where cos(γ ) is the declination. System's path angle for decline is obtained as: To illustrate true motion or air speed for drones, features for (V) are obtained. These features are either based on the horizontal velocity plane (V x ) or vertical velocity plane (V y ).
(V x ) is obtained either with the craft's rate of incline (RoI) and rate of decline (RoD). This is expressed as: 2) ACHIEVING MINIMUM TURNING RADIUS A competing drone system achieves a turn when it changes its course in a circular twist and over a pre-determined airspeed vector ratio, known as (V). When this is done, it ends up defining the drone's velocity. However, to accomplish this level twist in a circular path, the forces acting on the sides of the craft equals to 0, implying that the craft ends up tilting the lift (L) vector ratio in an angle (φ) as suggested by Crofton [23]. To generate the radial force (F), products of the titled lift (L) with the drone weight (W) force is obtained as n = L W . The Fig. 6a expresses the use of load factor (n) for the mathematical calculation of the radial acceleration, equaling the proportionality ratio of lift (L) and weight force (W). To illustrate the craft's flyable path measurement, a mathematical estimation of the craft's turn-rate (u), the bank angle (φ) and craft's turning radius (r) are induced. The craft's envelop or commencement strategy as consistent with the estimations are depicted in Fig. 6b. Parameters such as craft's limits and maximum drive velocity are illustrated as: where load peak factor (n) and its airspeed value (V) are placed at 2.45 m/s and 28.0 m/s. Assuming gravitational acceleration (g) equals 9.8 m/s 2 , then peak parameters for minimum turning radius, maximum turn rate and maximum bank-angle becomes: 37.0 m, 46.15 deg/sec and 62.75 deg/sec respectively.

3) RANDOM TRAJECTORY-AWARE COMMUNICATION LINK
A tightly-fixed illustration of the ergodic and the outage capacities for both the existing and the fully-recharged drones can be deduced. First, their respective distances using the Lemma 1 theorem which embeds the investigations made using Crofton Fixed Point Theorem is articulated. Lemma 1: If N vertices ζ equaling i = 1, 2, . . . , N , achieved a remote-distribution and are deployed at random in zone A, with volume of |A|, where height (H) is dependent on ζ 1 , . . . , ζ N . And that A ⊂ A, where δA denotes infinite boundary value of A, but not of A . On exploiting both Jensen's inequality distribution in [24] and Crofton's theorem of [11], the underlining postulation is made.
Postulation: The 4-D multi-drone using randomized motion trajectory and with the Rician fading propagation is determined thus: On careful observations, model's resulting outage capacity is: given that l th = ( τ 3 −1 τ 1 τ 2 τ 3 τ 4 ), and where; Proof: In proving (43), Euclidean distance among drones is denoted by L 1 , while fL 1 (l) denotes the PDF of L 1 value, while, P is probability between two drones spaced by distance l. Whereas, P 1 is known as probability of having atleast one drone around network area. Referencing Crofton's in Lemma 1, the numerical model is: Given |V| = volumetric 4-D sphere, i.e., |V| = 4 3 πr 3 s . d|V| = 4πr 2 s dr s . To evaluate P, we assume that at-least one drone is located within network zone, while at-least another drone is located just on the tip of network area and within radius l. The rest drones can be within or just outside network area. Exact location of strongest network point is dl, where volumetric network area for drones is expressed as 2πl 2 1 − l 2r s dl, implying that P 1 is numerically obtained as: (47) VOLUME 9, 2021 Average distance between existing and newly deployed drones are derived by numerical expansion of substituting (46) for (47) and further applying integral functions as: However, the existing Jensen's inequality in [25] is classified as being below-par, while its low bounds for the Ergodic capacity is now deduced as: The Ergodic derivatives is concluded since (43) is proved.

C. DRONE-TO-DRONE GUIDANCE MECHANISM
By utilizing the MATLAB simulink tool, investigations of 4-D drone guidance mechanism was accomplished in this section. All existing PPP models based on the already discussed randomized trajectory were consistent with our model. The along-track and cross-track errors, shown in Fig. 7 are the commonest of all errors that are prevalent in all drones and ways to permanently resolve these or to bring them to the barest minimum are solved using the embedded random trajectory tracking device. At any drone's point or location, vector error for cross-track is computed as the orthogonal projection depicted with positions P e (X e , Y e , Z e ) and is proportional to the winged drone's velocity ratio V r situated at a position P r (X r , Y r , Z r ) within a time-frame t, and in an earthbased coordinates of North-East-Down (NED). In contrast to cross-track error, the along-track error obtains the value difference between the winged tangential drone's velocity V r (V rx , V ry , V rz ) and its original velocity V e (V ex , V ey , V ez ) within the same time-frame t. A numerical estimation of the drone's relative velocity in NED coordinates is calculated assuming the drones's proximity as d = P r − P e , then the drone's original positions in NED will be expressed as dV = V r − V e , implying that NED conversion to frame R r from frame-size (u, v, n) is achieved. The frame R r with its frame size (u, v, n) are deduced using the drone's turning rotation and shown using matrix R as: where R is further enhancement of vector units with frame size (u, v, n), while u represents a parallel vector to V r , and n indicates a downward pointing vector. To mitigate errors in cross-track, a closed-loop of proportionate integral derivative (PID) controller is used. A closed-loop of proportionate (P) controller is used for along-track errors. Prior to this, operator ensures controllers are considerably raised to their respective optimal levels and allowed to implement velocity limits just to avoid excess speed for cross-track error, and to also implement angular and anti-wind features by avoiding steep/sharp turns for along-track error.

V. PERFORMANCE EVALUATION
This section discusses the deployment scenario and the performance evaluation of our proposed scheme. Details such as the system's operational procedure, investigations and graphical illustration/outputs were presented. Upon investigating the proposed hybrid (C-EKF), and as consistent with the detailed low-cost IMU devices of Table 2 and drone parameters of Table 3, overall system performance was then analyzed accordingly. Preliminary and final results presented optimal performance features at implementing the hybrid C-EKF over EKF and CKF schemes, in that order.

A. TESTBED DEPLOYMENT 1) DRONE SYSTEM FEATURES AND DEPLOYMENT a: MULTIPLE FLIGHT SPEEDS AND MODES
The E58 Pro Drone includes the new 3-D rolling special effects. One-key return function allows the drone to return automatically without losing the drone. The Headless mode option makes the drone's front side the same as the remote control, making it easier to fly to a preferred destination. With the Trimming function, the drones are adjusted according to tilt direction thus, making the drone stabilize during flight.

b: DEPLOYMENT, APP CONTROL, AND EASY-TO-USE MECHANICS
APPs are downloaded to control the drone. A click on the gyroscope icon flies the drone according to the gravity of the mobile phone. The drone can flip the screen 180 degrees when the REV icon is clicked. Click on  the VR icon to turn on mobile phone split screen-mode. VR glasses are essential to experience real-time transmission in 3-D visualization [19], [20]. Click on the Trajectory flight icon to draw flight trajectory on the phone. The drone will then follow the trajectory. On altitude hold-mode, operator can accurately lock the height and location, stable hover and capture video/images from any angle, making the experience quite easy and convenient. Non-experts can easily operate the drones. The drone automatically takes-off and lands with one handy click. There is an emergency landing button to prevent collision with obstacles.

2) TESTBED EVALUATION
A recapitulation of our work and its optimal performances were buttressed in the Fig. 8a -Fig. 8f, depicting the deployment, analysis and recorded outputs of the proposed hybrid C-EKF scheme. In Fig. 8a & Fig. 8b, the recreational drones were deployed and controlled using a console for its marked target goalpost. In Fig. 8c & Fig. 8d, a drone is seen to perfectly mirror the C-EKF computations and learn the target goalpost. The Fig. 8e depicts the drone's accuracy to it target goalpost. In Fig. 8f, perfection is achieved when a recreational-drone crossed the target goalpost's opening. The kinematic baselines and abilities of PPP accuracy of the deployed drones to conduct and complete a task achieved significant 45% improvement with the proposed hybrid C-EKF model. It is also observed that the system's speed-ratios improved while the heading convergence rate, using the lowcost INS/IMU was higher. Around 18% of the system's relative improvement was found in the average position error for 1-second GPS gaps using the hybrid C-EKF model versus 11% that was recorded using other filters. During stationary periods and without sufficient kinematic calibrations, a nearly 36% improvement in the heading component's drift rate was recorded.
B. NUMERICAL RESULTS AND DISCUSSION 1) YAW DIFFERENCE COMPARISON The Fig. 9 and Fig. 10 depicts the yaw differences in both stationary and in constant velocity, averaging a total of 20 filter runs. A system should maintain normal yaw to attain credibility. In Fig. 9, The output suggests that only the hybrid C-EKF maintained a normal yaw , hence, outperforming both EKF and CKF. When comparing the system in a constant velocity and in a stationary motion, end results suggests the optimality of the system in stationary motion over velocity mode.

2) OUTAGE PROBABILITY VS. SNR THRESHOLD ϑ
The system's outage capacity in (bit/s/Hz) is plotted against the SNR threshold ϑ ratios (dB) in Fig. 11, with an applicable Rician factor of K = 5 dB, P U = 0.1 W , and where r s = 500 m. An increased outage capacity of all the approach was recorded with a proportionate increase in threshold ϑ ratios. At the lowest and highest values of threshold ϑ, all the approaches maintained a steady increase until the value of ϑ = 6, when they all started displaying some variance in outage capacity. These values as recorded negates the 2-D approach with maintained a totally dissimilar ϑ ratevalue. It is seen that the 2-D approach maintained the lowest outage of 0.5 when ϑ = 0, and also maintained the highest outage capacity of nearly 2.5 when ϑ = 10 dB.
3) ERGODIC CAPACITY VS. SNR THRESHOLD ϑ Fig. 12 depicts a confluence comparison of the ergodic capacity against the SNR threshold ϑ, with an application of the Rician factor where K = 5 dB and P U = 0.01 W , and ϑ th = 0 dB. The varying densities to the network link is evaluated by implementing a multi-hop (relay-based) and direct Ad-Hoc communication linking network. It is observed that, with a multi-hop (relay-based) ergodic network, not much rapid changes occur. However, with direct communication link of both C erg6 and C erg1 , under a pre-determined radius, rapid changes occurred. A sharp decline of values was recorded with an increased radius r s , when in C erg1 mode for all the approaches. The decreased values, as recorded were consistent and was never broken for all the approaches. The peak point of zero (0) value was reached for the 4-D approach, when the coverage radius was at its widest radius of nearly 600, as against other 3-D and 2-D approaches. However for direct C erg6 ergodic capacity, little changes occurred as no sharp decline was recorded, even its radius at its widest. As seen, the values maintained a nearly-linear value for all the approaches and only started to considerably decline when the radius r s was nearing the mark of 300. It can also be observed that a confluence both the ergodic values of C erg6 and C erg1 , occurred when at a radius of 350 m, 400 m and 500 m, for analytical at N = 50, analytical at N = 20 and simulation at N = 20, respectively. There was no confluence of system for when simulation at N = 50.

C. ILL-CONDITIONING COMPUTATIONS
The numerical testbed computation procedure is organized as follows. The unknown system parameter θ is estimated from all available experimental data Z N 1 = {z 1 , z 2 , . . . , z N }, through gradient-based adaptive KF technique and as consistent with the works of [9], [26], [27]. For easier analysis, all differentiated algorithms observes the same data process Z N 1 while the initial value of the optimization method, θ (0) = 1. Next, the obtained optimal estimateθ ( * ) is then computed, while comparison is drawn against the ''real'' value ofθ ( * ) = 5 for ascertaining each approach's estimation quality. The experiment is repeated for M = 100 times and then numerically investigate the estimate's a posterior mean, the root-mean-squared error RMSE and the mean absolute percentage error (MAPE%) over 100 Monte Carlo runs.
The results obtained in Table 4 depicts careful experimental observations as consistent with the available data and MATLAB computations. To this effect, the following conclusions are presented. Firstly, when not ill-conditioned i.e δ at 10 −1 and 10 −2 , all differentiated variants of KF work equal and normal, ascertaining that the techniques are algebraically similar. Secondly, the differential KF is the least performing approach, implying its fastest degradation when δ = ∈ roundoff . A computer roundoff for floating-point arithmetic is characterized by a single parameter ∈ roundoff , that is often defined in different sources as the largest number such that 1 + ∈ roundoff = 1 or 1 + ∈ roundoff /2 = 1 in machine precision terms. The unaccounted lines for the differential KF and differential CKF approaches suggests MATLAB's inability to compute the algorithms. Thirdly, on investigation, it is observed that the UD-based CKF approach marginally outperformed the conventional differential KF approach. It posted a bit more robust computations with range of ill-conditioning put at 10 −10 in comparison to the KF with range of just 10 −6 . However, our proposed CKF-EKF model not only outperformed UD-based EKF and KF approaches, it also exhibited robust and quite competitive MAPE%, RMSE and mean values against the others.

VI. OPEN RESEARCH ISSUES A. INTERFERENCE LIMITATION FROM BASE STATION
In [28]- [31], authors confirm that adhoc interference vulnerability of drones with remotely located ground station (GS) transmitters are major research concern. The effects of RF interferences in our proposed 4-D based competing recreational drones are investigated starting with a step to illustrate the drone's strongest received adhoc signal as: y = √ P U h 1 g 1 x(t) + √ P I h I g I x(t) + n(t), where ground station transmitter power is P I , channel fading coefficient from the transmitter to drones is h I , path loss equaling g 1 = L −α 1 is denoted, where L 1 is spherical proximity between independent drones. Similarly, path loss equaling g 1 = L −αI 3 is the proximity between a drone and its designated GS, where αI is the path loss exponent.
To estimate the system's signal-to-interference-plus-noise (SINR) denoted as (Z), an expression is formulated as: where a denotation of (X ) as (P U L −α 1 |h 1 | 2 ), while that of the (Y ) as (P I L −α I 3 |h I | 2 + σ 2 ) are formulated. But with channel model h I that VOLUME 9, 2021 observes Rician distribution, where parameter K denotes K I , an expression for PDF is then deduced as: , where x is denoted by (P U L −α 1 |h 1 | 2 . In a similar manner, the PDF of Y can be formulated with respective replacements of K, x and X with K I , y and Y. On further mathematical expansions, a probabilistic expectation for X is generated as:X = E[X ] = where Jensen's inequality ratio expressions [32] are core inequality generators.

B. ERROR-PRONE CSI AT THE PPP'S RECEIVER-END
Error-prone channel estimation are prevalent as a result of the stipulated requirements of the minimum mean-square error when they are being projected. Recently, authors in [32]- [36]. have modeled and projected these error-prone CSI to be h 1 =ĥ 1 + h e , where the receiver side's estimated CSI isĥ 1 , the zero-mean value is h e , the Gaussian estimation complex error is h 1 , and the CSI's variance is . To now calculate the SNR for the received average, these expressions are presented as follows:θ = P U L −α 1 | h 1 | 2 P U L −α 1 +σ 2 = P U | h 1 | 2 P U +σ 2 L α 1 , However, recall that the ergodic capacity C erg6 of the system is essential and deduced as: C erg6 = C erg1 (θ ) = 1 ln (2)  ].
The lower bounds of the ergodic capacity are further analyzed using Jensen's inequality formulation method expressed as: The two discussed open issues above were resolved using Jensen's inequality model, together with the enhanced computation C-EKF model as stated in the previous sections. They were also instrumental in shaping this work. It must be mentioned in clear terms that other open-issues are very much prevalent and they are subjective to the operations and research investigations of all other authors.

VII. CONCLUSION
This work's major objectives and accomplishments were: • To design and deploy hybrid C-EKF model into a GPS-IMU tightly-coupled four-dimensional (4-D) motion recreational drones to further enhance their ability to cooperate and execute specified tasks.
• To analyze network communication link capacity between drones and base stations (BSs) using randomized 4-D trajectory model, as consistent with simulation.
• To further improve GPS-IMU reliability, which is utilized for mapping applications in confined environment where frequent losses of GPS lock and PPP inaccuracy due to long GPS baselines are common occurrence. Furthermore, this article recapitulates the derivatives for proposed C-EKF for modern GNSS-enabled drone systems, while demonstrating its inherent real-time and computational load capacities. The PPP's kinematic baselines and abilities for the mobility system to conduct and complete a task achieved significant 45% improvement with the proposed model. It is also observed that significant improvement of the system's speed ratios was attained while its heading convergence rate remained higher when IMUs are introduced. Up to 18% of the system's relative improvement was achieved in the average position error for every 1 Sec GPS gap versus the 11% that was recorded using other filters. During stationary periods and without sufficient kinematic calibrations, up to 36% improvement in the drift rate of the heading component was also recorded. For future research considerations, concerted effort at resolving imperfect channel state information (CSI) estimations and system signal de-noising from ground station terminals are issues to be investigated.