In modern defense scenarios, the paradigm of unmanned aerial vehicle (UAV) operations is shifting from single-platform missions to intelligent swarm cooperation, particularly for intelligence, surveillance, and reconnaissance in high-threat environments. To minimize platform vulnerability while ensuring persistent monitoring and precise tracking beyond the defensive perimeter, we must address several interrelated challenges: target state uncertainty, strict standoff distance constraints, and the difficulty of coordinating heterogeneous sensors across multiple UAV drone platforms. Traditional single-UAV drone tracking often requires high-performance sensors, increasing cost, and lacks redundancy. Multi-UAV drone cooperation offers a promising solution by leveraging distributed perception, configuration optimization, and functional complementarity. However, achieving effective standoff cooperative tracking demands rapid convergence from distant approach to stable orbiting, balancing safety with observation quality, and adaptive decision-making under continuous target maneuvers. In this work, we propose a method that integrates receding horizon control (RHC) with an improved heuristic optimizer to maximize the cumulative Fisher information over the prediction horizon, thereby enabling autonomous formation and maintenance of near-optimal observation geometry for a heterogeneous-sensor multi-UAV drone system.
System Modeling and Problem Formulation
We consider a scenario where a team of \(N\) UAV drones cooperatively track a single ground moving target. The target state is denoted as \(\mathbf{o} = [x_t, y_t]^\top\) in the horizontal plane. Each UAV drone \(i\) is modeled as a fixed-wing aircraft with state \(\mathbf{s}_i = [x_i, y_i, \psi_i, v_i, \omega_i]^\top\), where \((x_i, y_i)\) is its planar position, \(\psi_i\) the heading angle, \(v_i\) the speed, and \(\omega_i\) the yaw rate. All UAV drones fly at a constant altitude \(h\), forming an air-ground observation geometry. The standoff constraint requires the horizontal distance \(r_i = \sqrt{(x_i – x_t)^2 + (y_i – y_t)^2}\) to satisfy \(r_i \ge R_s\) for all UAV drones, where \(R_s\) is the defensive perimeter radius.
The UAV drone dynamics are described by a first-order model:
$$
\dot{\mathbf{s}}_i = f(\mathbf{s}_i, \mathbf{u}_i) =
\begin{bmatrix}
v_i \cos\psi_i \\ v_i \sin\psi_i \\ \omega_i \\
-\frac{1}{\tau_v} v_i + \frac{1}{\tau_v} u_{v,i} \\
-\frac{1}{\tau_\omega} \omega_i + \frac{1}{\tau_\omega} u_{\omega,i}
\end{bmatrix}
$$
where \(\tau_v, \tau_\omega\) are time constants and \(\mathbf{u}_i = [u_{v,i}, u_{\omega,i}]^\top\) are control inputs. Constraints include \( |u_{v,i} – v_0| \le v_{\max}\) and \( |u_{\omega,i}| \le \omega_{\max}\). The model is discretized using Euler integration with sampling time \(T_s\).
Each UAV drone carries heterogeneous sensors: a range sensor and a bearing sensor. The measurement vector for UAV drone \(i\) at time \(k\) is
$$
\mathbf{z}_i(k) =
\begin{bmatrix}
z_{r,i}(k) \\ z_{\phi,i}(k)
\end{bmatrix} =
\begin{bmatrix}
r_i(k) \\ \phi_i(k)
\end{bmatrix} +
\begin{bmatrix}
\varepsilon_{r,i}(k) \\ \varepsilon_{\phi,i}(k)
\end{bmatrix}
$$
where \(\varepsilon_{r,i} \sim \mathcal{N}(0, \sigma_{r,i}^2)\) and \(\varepsilon_{\phi,i} \sim \mathcal{N}(0, \sigma_{\phi,i}^2)\). The joint Jacobian matrix for UAV drone \(i\) is:
$$
\mathbf{H}_i =
\begin{bmatrix}
\frac{\Delta x_i}{r_i} & \frac{\Delta y_i}{r_i} \\[4pt]
-\frac{\Delta y_i}{r_i^2} & \frac{\Delta x_i}{r_i^2}
\end{bmatrix}
$$
with \(\Delta x_i = x_t – x_i\), \(\Delta y_i = y_t – y_i\).
We adopt the Extended Information Filter (EIF) for state estimation due to its additive information form. The information matrix \(\mathbf{Y}(k|k)\) is updated by fusing contributions from all UAV drones:
$$
\mathbf{Y}(k|k) = \mathbf{Y}(k|k-1) + \sum_{i=1}^N \mathbf{I}_i(k)
$$
where \(\mathbf{I}_i(k) = \mathbf{H}_i^\top(k) \mathbf{R}_i^{-1}(k) \mathbf{H}_i(k)\) and \(\mathbf{R}_i = \text{diag}(\sigma_{r,i}^2, \sigma_{\phi,i}^2)\). For a single-sensor UAV drone (e.g., only range or only bearing), the missing measurement noise is set to infinity, effectively ignoring that dimension.
Optimal Observation Configuration Analysis
The Fisher Information Matrix (FIM) of the overall system is the sum of individual contributions. For a single UAV drone with both sensors, the FIM is
$$
\mathbf{J}_i = \mathbf{H}_i^\top \mathbf{R}_i^{-1} \mathbf{H}_i =
\begin{bmatrix}
\frac{\sin^2\phi_i}{\sigma_{r,i}^2} + \frac{\cos^2\phi_i}{r_i^2 \sigma_{\phi,i}^2} & \frac{\sin\phi_i \cos\phi_i}{\sigma_{r,i}^2} – \frac{\sin\phi_i \cos\phi_i}{r_i^2 \sigma_{\phi,i}^2} \\
\frac{\sin\phi_i \cos\phi_i}{\sigma_{r,i}^2} – \frac{\sin\phi_i \cos\phi_i}{r_i^2 \sigma_{\phi,i}^2} & \frac{\cos^2\phi_i}{\sigma_{r,i}^2} + \frac{\sin^2\phi_i}{r_i^2 \sigma_{\phi,i}^2}
\end{bmatrix}
$$
The determinant of \(\mathbf{J}_i\) is \(\det(\mathbf{J}_i) = \frac{1}{r_i^2 \sigma_{r,i}^2 \sigma_{\phi,i}^2}\). Hence, for a single UAV drone, the observation quality improves as \(r_i\) decreases, bounded by the standoff constraint \(r_i \ge R_s\), leading to \(\max \det(\mathbf{J}_i) = 1/(R_s^2 \sigma_r^2 \sigma_\phi^2)\).
For a multi-UAV drone system with homogeneous sensors (identical noise variances), the total FIM determinant becomes
$$
\det(\mathbf{J}_{\text{total}}) = \frac{1}{\sigma_\phi^2 \sigma_r^2} \sum_{i} \frac{1}{r_i^2} + \frac{1}{\sigma_\phi^2 \sigma_r^2} \sum_{i<j} $$="" +="" -="" <p="" \(\phi_{ij}="\phi_i" \(\sin^2\phi_{ij}\).="" \cos^2\phi_{ij}="" \frac{1}{\sigma_\phi^4}="" \frac{1}{\sigma_r^4}="" \frac{1}{r_i^2="" \frac{1}{r_i^2}="" \frac{1}{r_j^2}="" \left(="" \phi_j\).="" \right)="" \sin^2\phi_{ij}="" \sum_{iFor heterogeneous sensors, we introduce parameters
$$
a_i = \frac{1}{\sigma_{r,i}^2}, \quad b_i = \frac{1}{r_i^2 \sigma_{\phi,i}^2}, \quad d_i = a_i – b_i
$$
Then the total FIM determinant can be expressed as:
$$
\det(\mathbf{J}_{\text{total}}) = \frac{1}{4} \left[ \left( \sum_i (a_i+b_i) \right)^2 – \sum_i d_i^2 – 2 \sum_{i<j} $$="" (a_i+b_i)="" -="" 2\phi_j)="0" 2\phi_{ij}="" 2\phi_{ij}\).="" <p="" \(\left(="" \(\sum_{iThe parameter \(d_i\) characterizes the sensor’s preference: \(d_i > 0\) indicates a range-favoring sensor, \(d_i < 0\) a bearing-favoring sensor, and \(d_i = 0\) a balanced sensor. Key properties emerge:
- Property 1: When \(d_i\) have the same sign, the optimal angular separation between any two UAV drones tends to 90°, maximizing information orthogonality.
- Property 2: Among same-sign sensors, those with larger \(|d_i|\) tend to separate their angles more aggressively.
- Property 3: When \(d_i\) have opposite signs, the optimal angular separation tends to 0°, enabling complementary range-bearing fusion.
- Property 4: High-precision sensors (larger \(|d_i|\)) dominate the overall configuration.
Proposed Method: RHC with Improved Snake Optimizer
We formulate the cooperative tracking problem as a receding horizon optimization. At each time step, each UAV drone predicts its future trajectory over a horizon \(N_p\) based on candidate control sequences \(\mathbf{U}_i = [\mathbf{u}_i(k), \mathbf{u}_i(k+1), \ldots, \mathbf{u}_i(k+N_p-1)]^\top\). The objective is to maximize the cumulative FIM determinant over the horizon while satisfying dynamics, control bounds, and standoff constraints. The cost function for UAV drone \(i\) is:
$$
J_{\text{pre}}(\mathbf{U}_i) = \sum_{l=k}^{k+N_p-1} \sum_{j=1}^N \det\left( \mathbf{J}_j(l|\mathbf{U}_i) \right)
$$
To convert maximization into minimization, we use the negative logarithm:
$$
C_{\text{total}}(\mathbf{U}_i) = – w_{\text{obs}} \log \det\big( \mathbf{J}_{\text{pre}}(\mathbf{U}_i) \big) + w_{\text{safe}} C_{\text{safe}}(\mathbf{U}_i)
$$
where \(C_{\text{safe}}\) penalizes predicted positions that penetrate obstacle zones modeled as circular regions with collision radius \(R_c\) and avoidance radius \(R_a\) (\(R_a > R_c\)). The safety cost at each step is:
$$
f_{\text{safe}}(\delta_{i,l}) = \begin{cases}
C_h, & \delta_{i,l} \le 0 \\
C_s, & 0 < \delta_{i,l} < R_a \\
0, & \delta_{i,l} \ge R_a
\end{cases}
$$
where \(\delta_{i,l}\) is the distance from UAV drone \(i\) at step \(l\) to the obstacle boundary. The total safety cost accumulates over the horizon and all obstacles.
To ensure smooth transition from approach to orbiting and maintain direction consistency, we introduce a consistency penalty and an orbital band screening mechanism. The directional consistency ratio \(\rho\) is computed over the predicted horizon:
$$
\rho = \frac{1}{N_p-1} \sum_{l=1}^{N_p-1} \mathbb{I}\big[ \text{sgn}(\Delta \phi_k) = 1 \big]
$$
A low \(\rho\) triggers a penalty that biases the optimizer toward trajectories with consistent encircling direction. Furthermore, we apply an annular region screening: candidate trajectories whose predicted positions fall within a narrow band just outside the standoff radius \(R_s\) are evaluated with higher priority; those outside the band are postponed. This mechanism helps the swarm converge to and stay near the optimal observation radius.
Additionally, a smoothness term is added to avoid abrupt heading changes:
$$
C_{\text{orbit}} = \sum_{l=1}^{N_p-1} \| \phi_{l+1} – \phi_l \|
$$
This term is included in the total cost to promote continuous transition.
The resulting optimization problem is high-dimensional, non-convex, and constrained. We employ an Improved Snake Optimizer (ISO) as the solver. ISO simulates snake mating and foraging behavior with gender-based population division and two-phase search: exploration and exploitation. During exploration (control parameter \(E < 0.25\)), male and female groups generate new candidates via multi-strategy chaotic sampling:
$$
X_{m,n+1} = \left| \sin\big( \pi ( 4 X_{m,n}(1-X_{m,n}) + \sin(\pi X_{m,n}) ) \big) \right|
$$
for males, and a symmetric mechanism for females. In the exploitation phase (\(E \ge 0.25\)), the algorithm switches between hot and cold environments based on temperature \(T_{\text{emp}} = \exp(-t/T)\). In hot conditions (\(T_{\text{emp}} > 0.6\)), all individuals converge toward the global best. In cold conditions, with 60% probability, fighting mode is activated where individuals learn from opposite-gender best; otherwise, a bidirectional Pareto evolution (20/80 rule) refines the population. ISO has demonstrated strong global search and local refinement capabilities for complex engineering problems. At each RHC cycle, ISO returns the optimal control sequence, and only the first control input \(u^*_v(0), u^*_\omega(0)\) is applied.
Simulation Results and Discussion
We conduct extensive simulations in a ground vehicle tracking scenario. The simulation parameters are summarized in Table 1.
| Parameter | Value |
|---|---|
| Number of iterations | 40 |
| Population size | 40 |
| Prediction horizon \(N_p\) | 20 |
| Sampling period \(T_s\) | 0.5 s |
| Simulation duration | 1000 s |
| Speed time constant \(\tau_v\) | 6 |
| Yaw rate time constant \(\tau_\omega\) | 6 |
| Speed bounds | (10 m/s, 30 m/s) |
| Yaw rate bounds | (-0.12 rad/s, 0.12 rad/s) |
| Max acceleration | <5 m/s² |
| Max angular acceleration | <0.02 rad/s² |
| UAV start position | (0, 0) |
| Target start position | (5000 m, 0) |
| Target velocity vector | (1.0, 0.2) m/s |
| Base range noise std \(\sigma_r\) | 20 m (\(\sigma_r^2=400\) m²) |
| Base bearing noise std \(\sigma_\phi\) | 5° |
First, we verify the single-UAV drone case: the optimal strategy is to fly as close to the standoff radius as possible, achieving the maximum FIM determinant. Next, for two and three homogeneous UAV drones, the optimal angular separations are 90° and 60°, respectively. Our method successfully drives the formation to these configurations, as shown by the measured angles converging to the theoretical values (Table 2).
| Number of UAV drones | Avg. estimation error (overall) / m | Avg. estimation error (orbiting phase) / m |
|---|---|---|
| 1 | 23.74 | 7.96 |
| 2 | 10.07 | 4.78 |
| 3 | 9.23 | 4.08 |
| 4 | 8.37 | 3.57 |
| 5 | 7.61 | 3.38 |
The results demonstrate that as the number of UAV drones increases, the estimation accuracy improves, but the gain saturates beyond three. For heterogeneous sensor configurations, we examine two scenarios: same-sign \(d_i\) vs. opposite-sign \(d_i\). Table 3 compares the estimation errors between the optimal angular configuration and the equally spaced configuration, highlighting the benefit of phase optimization.
| Configuration | Error (optimal) / m | Error (equal-angle) / m |
|---|---|---|
| Two UAV drones, opposite sign \(\sigma_r^2=[1000,10]\), \(\sigma_\phi^2=[0.1,1]\) |
0.0233 | 0.5076 |
| Three UAV drones, opposite sign \(\sigma_r^2=[1000,100,10]\), \(\sigma_\phi^2=[0.01,0.1,1]\) |
0.4678 | 0.7027 |
| Three UAV drones, same sign \(\sigma_r^2=[100,100,100]\), \(\sigma_\phi^2=[0.11,0.1,0.09]\) |
0.3336 | 0.3952 |
When \(d_i\) have opposite signs, the optimal configuration deviates significantly from uniform spacing, leading to substantial improvement in accuracy. For same-sign configurations, the improvement is more modest because the optimal angles are close to uniform.
We compare our proposed method (ISO-RHC) with the classic Lyapunov Vector Field Guidance (LVFG) method for two homogeneous UAV drones. The convergence time (to within ±3° of optimal phase) and the approaching phase error are listed in Table 4. Our method achieves virtually zero convergence time because the formation is already optimal before reaching the standoff orbit.
| Method | Convergence time / s | Angle error (orbiting) / deg | Approach phase avg. error / m | Time to reach standoff / s |
|---|---|---|---|---|
| ISO-RHC | ≈ 0 | 1.51 | 13.10 | 254 |
| LVFG | 90 | 1.81 | 17.20 | 315 |
Additionally, we test obstacle avoidance scenarios. The annular screening mechanism allows UAV drones to temporarily deviate from the standoff orbit to avoid obstacles while maintaining formation stability. The approach also demonstrates robustness under random target motions (smooth and aggressive).
Conclusion
In this work, we addressed the cooperative observation and tracking problem for a multi-UAV drone system equipped with heterogeneous sensors operating beyond a defensive perimeter. By analytically deriving the Fisher Information Matrix for heterogeneous sensor configurations, we identified the optimal angular geometry and key properties governing formation design. We then proposed a receding horizon control framework coupled with an improved snake optimizer, incorporating an orbital consistency mechanism and smooth transition strategy. The method autonomously guides the UAV drone swarm from distant approach to stable orbiting while maximizing cumulative information gain. Simulation results demonstrate significant improvements in tracking accuracy and robustness compared to the Lyapunov vector field method, especially in heterogeneous sensor scenarios. Future work will extend the approach to full 3D environments, account for heterogeneous UAV drone dynamics, and integrate Doppler velocity sensors for enhanced performance.