In recent years, the operational deployment of multiple unmanned aerial vehicles (UAVs) in coordinated groups, known as drone formations, has transitioned from a research concept to a critical capability in both military and civilian domains. The core advantage of a drone formation lies in its ability to perform tasks—such as wide-area surveillance, distributed sensing, cooperative payload transport, and complex tactical maneuvers—with significantly greater efficiency, redundancy, and spatial coverage than a single UAV. Maintaining the integrity and coordinated behavior of the drone formation is paramount; the failure or anomalous behavior of even a single unit can compromise the entire mission. Consequently, the real-time detection of anomalies within a drone formation becomes a fundamental requirement for ensuring safety, reliability, and mission success.
Traditional anomaly detection methods often rely on a multi-hypothesis testing framework grounded in statistical decision theory. These approaches require explicit probabilistic models of both normal and faulty conditions, and their complexity grows combinatorially with the number of UAVs in the drone formation. To circumvent this complexity and the need for precise prior knowledge of each UAV’s fault modes, this article advocates for a purely data-driven methodology. The central philosophy is to treat anomaly detection as a problem of model identification and consistency checking across the drone formation. Instead of looking for predefined fault signatures, we aim to learn a behavioral model for a normal UAV directly from its observed input-output flight data. Anomalies are then detected by identifying UAVs whose behavior significantly deviates from this learned norm or from the consensus model of the drone formation.
Consider a drone formation consisting of \(N\) homogeneous UAVs. While complex, the flight dynamics of each UAV \(i\) can often be related through fundamental physics to a set of measurable inputs and outputs. A critical relationship, often derived from force and moment balance equations, links the vehicle’s acceleration vector to its state and control inputs. This relationship, though nonlinear and partially unknown, is structurally identical for all homogeneous UAVs under nominal conditions. We can express this generic relationship for the \(i\)-th UAV at discrete time \(t\) as:
$$ \mathbf{y}_i(t) = \mathbf{F}(\boldsymbol{\psi}_i(t)) + \mathbf{e}_i(t), \quad i=1,\ldots,N $$
where:
- \(\mathbf{y}_i(t) \in \mathbb{R}^{n_y}\) is the output vector, typically containing linear and angular accelerations.
- \(\boldsymbol{\psi}_i(t)\) is the regressor vector containing relevant inputs and possibly past outputs (e.g., angle of attack, velocity, angular rates, dynamic pressure, control surface deflections, thrust).
- \(\mathbf{F}(\cdot)\) is an unknown, smooth nonlinear function representing the nominal dynamics, shared across the drone formation.
- \(\mathbf{e}_i(t)\) is a zero-mean, finite-variance measurement or process noise sequence, assumed to be white and independent of \(\boldsymbol{\psi}_i(t)\).
The data available for analysis from the drone formation is the massive collection of input-output sequences recorded during flight operations:
$$ \mathcal{D} = \{ \boldsymbol{\psi}_i(t), \mathbf{y}_i(t) \}_{t=1}^{M} \quad \text{for} \quad i=1,\ldots,N $$
where \(M\) is the number of data samples per UAV. The drone formation anomaly detection problem is thus: Using only the dataset \(\mathcal{D}\), identify UAVs whose underlying function \(\mathbf{F}_i(\cdot)\) has deviated from the common \(\mathbf{F}(\cdot)\).
Data-Driven Model Identification via Basis Function Expansion
The first step in our data-driven approach is to approximate the unknown nonlinear function \(\mathbf{F}(\cdot)\). A powerful and flexible method is to expand it onto a basis of known functions. Let us consider a scalar output channel for simplicity (the multi-output case follows similarly via separate expansions). The model for the \(i\)-th UAV’s \(j\)-th output becomes:
$$ y_i^{(j)}(t) = \sum_{l=1}^{n} \theta_{i,l}^{(j)} \phi_l(\boldsymbol{\psi}_i(t)) + e_i^{(j)}(t) = \boldsymbol{\phi}^T(\boldsymbol{\psi}_i(t)) \boldsymbol{\theta}_i^{(j)} + e_i^{(j)}(t) $$
Here, \(\{\phi_l(\cdot)\}_{l=1}^{n}\) is a set of predetermined basis functions (e.g., polynomials, radial basis functions, wavelets). \(\boldsymbol{\phi}(\boldsymbol{\psi}_i(t)) = [\phi_1(\cdot), \ldots, \phi_n(\cdot)]^T\) is the known regression vector, and \(\boldsymbol{\theta}_i^{(j)} \in \mathbb{R}^n\) is the unknown parameter vector specific to UAV \(i\) and output \(j\). For a homogeneous drone formation operating normally, we expect \(\boldsymbol{\theta}_i^{(j)} \approx \boldsymbol{\theta}_0^{(j)} \ \forall i\), where \(\boldsymbol{\theta}_0^{(j)}\) is the “true” shared parameter vector.
Collecting data over \(M\) time steps for UAV \(i\), we can write the matrix equation:
$$ \mathbf{Y}_i^{(j)} = \boldsymbol{\Phi}_i \boldsymbol{\theta}_i^{(j)} + \mathbf{E}_i^{(j)} $$
where:
$$
\mathbf{Y}_i^{(j)} = \begin{bmatrix} y_i^{(j)}(1) \\ \vdots \\ y_i^{(j)}(M) \end{bmatrix}, \quad
\boldsymbol{\Phi}_i = \begin{bmatrix} \boldsymbol{\phi}^T(\boldsymbol{\psi}_i(1)) \\ \vdots \\ \boldsymbol{\phi}^T(\boldsymbol{\psi}_i(M)) \end{bmatrix}, \quad
\mathbf{E}_i^{(j)} = \begin{bmatrix} e_i^{(j)}(1) \\ \vdots \\ e_i^{(j)}(M) \end{bmatrix}
$$
A natural and computationally efficient way to estimate \(\boldsymbol{\theta}_i^{(j)}\) is the Ordinary Least Squares (OLS) method:
$$ \hat{\boldsymbol{\theta}}_i^{(j)} = \arg\min_{\boldsymbol{\theta}} \lVert \mathbf{Y}_i^{(j)} – \boldsymbol{\Phi}_i \boldsymbol{\theta} \rVert^2 = (\boldsymbol{\Phi}_i^T \boldsymbol{\Phi}_i)^{-1} \boldsymbol{\Phi}_i^T \mathbf{Y}_i^{(j)} $$
The quality and properties of this estimate are paramount for effective anomaly detection within the drone formation.
| Basis Type | Function Form \(\phi_l(\psi)\) | Properties | Suitability for Drone Formation Dynamics |
|---|---|---|---|
| Polynomials | \(\psi^l, \ \psi_1^{l_1}\psi_2^{l_2}\ldots\) | Simple, global support. Prone to Runge’s phenomenon for high order. | Good for mild nonlinearities in aerodynamic coefficients. |
| Radial Basis Functions (RBF) | \(\exp(-\lVert \psi – c_l \rVert^2 / \sigma^2)\) | Local support, excellent approximation properties. Requires center \(c_l\) selection. | Excellent for capturing complex, localized nonlinear flight dynamics. |
| B-Splines | Piecewise polynomial | Local support, smoothness control. Requires knot placement. | Suitable for representing tabular aerodynamic data. |
Anomaly Detection under Unbiased Estimation Conditions
The core idea for anomaly detection in the drone formation is to analyze the residuals—the differences between the measured outputs and the model predictions. Under the assumption of a correct model structure and unbiased parameter estimates, these residuals should behave like the white noise sequence \(\mathbf{e}_i(t)\). A significant deviation in the residual statistics for a particular UAV indicates anomalous behavior.
We define the residual sequence for UAV \(i\) and output \(j\) as:
$$ \epsilon_i^{(j)}(t) = y_i^{(j)}(t) – \boldsymbol{\phi}^T(\boldsymbol{\psi}_i(t)) \hat{\boldsymbol{\theta}}_i^{(j)} $$
For the drone formation, we need to check if the residuals for each UAV are consistent with the noise properties. A common test statistic is based on the weighted norm of the residuals. Let us estimate the residual covariance matrix \(\mathbf{W}\) pooled across the drone formation (under the null hypothesis of no anomalies):
$$ \hat{\mathbf{W}} = \frac{1}{N(M-n)} \sum_{i=1}^{N} \sum_{t=1}^{M} \boldsymbol{\epsilon}_i(t) \boldsymbol{\epsilon}_i^T(t) $$
where \(\boldsymbol{\epsilon}_i(t)\) is the vector of residuals for all outputs. A per-UAV, per-sample anomaly detector \(D(\cdot)\) can be designed as:
$$ D(\boldsymbol{\epsilon}_i(t)) = \begin{cases}
1 \text{ (Anomaly)}, & \text{if } \boldsymbol{\epsilon}_i^T(t) \hat{\mathbf{W}}^{-1} \boldsymbol{\epsilon}_i(t) > \tau \\
0 \text{ (Normal)}, & \text{otherwise}
\end{cases} $$
where \(\tau\) is a detection threshold chosen based on the desired false alarm rate (e.g., from the \(\chi^2\) distribution if residuals are Gaussian). For robustness, one typically uses a moving average or cumulative sum (CUSUM) of this statistic over a time window.
The critical requirement for this detector to function correctly is that the OLS estimate \(\hat{\boldsymbol{\theta}}_i^{(j)}\) must be unbiased. The bias of the OLS estimate is given by:
$$ \mathbb{E}[\hat{\boldsymbol{\theta}}_i^{(j)}] – \boldsymbol{\theta}_0^{(j)} = \mathbb{E}[(\boldsymbol{\Phi}_i^T \boldsymbol{\Phi}_i)^{-1} \boldsymbol{\Phi}_i^T \mathbf{E}_i^{(j)}] $$
This bias is zero if the regressor \(\boldsymbol{\phi}(\boldsymbol{\psi}_i(t))\) is uncorrelated with the noise \(e_i^{(j)}(t)\). This condition holds if the regressor contains only exogenous inputs (e.g., control commands, measured states like velocity). In such a scenario, OLS provides consistent estimates, and the residual-based anomaly detector for the drone formation is statistically sound.
| Step | Action | Mathematical Expression / Description |
|---|---|---|
| 1. Offline Modeling | Estimate shared model from nominal drone formation data. | \(\hat{\boldsymbol{\theta}}_0 = \frac{1}{N}\sum_{i=1}^N (\boldsymbol{\Phi}_i^T \boldsymbol{\Phi}_i)^{-1} \boldsymbol{\Phi}_i^T \mathbf{Y}_i\) (or pooled data). |
| 2. Online Residual Gen. | For each UAV \(i\) at time \(t\), compute residual. | \(\boldsymbol{\epsilon}_i(t) = \mathbf{y}_i(t) – \boldsymbol{\Phi}(\boldsymbol{\psi}_i(t)) \hat{\boldsymbol{\theta}}_0\). |
| 3. Test Statistic | Calculate the normalized energy of the residual. | \(T_i(t) = \boldsymbol{\epsilon}_i^T(t) \hat{\mathbf{W}}^{-1} \boldsymbol{\epsilon}_i(t)\). |
| 4. Decision | Compare to threshold. Flag anomaly if exceeded persistently. | If \(\frac{1}{L}\sum_{k=t-L+1}^t T_i(k) > \tau\), flag UAV \(i\) as anomalous. |
The Challenge of Biased Estimation and the Need for Compensation
In practical drone formation dynamics, accurate models often require regressors that include past output values (e.g., \(\boldsymbol{\psi}_i(t) = [\mathbf{y}_i(t-1), \mathbf{u}_i(t)]^T\)). This leads to an Auto-Regressive with eXogenous input (ARX) type model structure in the linear case or its nonlinear counterpart. In this situation, the regressor \(\boldsymbol{\phi}(\boldsymbol{\psi}_i(t))\) becomes correlated with the noise \(e_i^{(j)}(t)\) because \(\boldsymbol{\psi}_i(t)\) contains past outputs \(y_i(t-1)\), which depend on past noise \(e_i(t-1)\). This correlation violates the OLS assumption, leading to a biased estimate:
$$ \mathbb{E}[\hat{\boldsymbol{\theta}}_i^{(j)}] – \boldsymbol{\theta}_0^{(j)} = \mathbb{E}[(\boldsymbol{\Phi}_i^T \boldsymbol{\Phi}_i)^{-1} \boldsymbol{\Phi}_i^T \mathbf{E}_i^{(j)}] \neq \mathbf{0} $$
This bias, often called the “equation error bias” or “noise-induced bias,” does not vanish even with infinite data. Consequently, the parameter estimates for different UAVs in the drone formation will converge to different biased values even if their true dynamics \(\boldsymbol{\theta}_0^{(j)}\) are identical. This renders the simple residual-based anomaly detector unreliable, as large residuals may be caused by estimation bias rather than true dynamic anomalies.
To illustrate, consider a simple nonlinear ARX model for the drone formation:
$$ y_i(t) = \theta_{0,1} y_i(t-1) + \theta_{0,2} y_i^2(t-1) + \theta_{0,3} u_i(t) + e_i(t) $$
The regressor is \(\boldsymbol{\phi}(t) = [y_i(t-1), y_i^2(t-1), u_i(t)]^T\). Clearly, \(y_i(t-1)\) is correlated with \(e_i(t-1)\), making OLS biased.
Bias-Compensated Ordinary Least Squares (BCOLS) for Drone Formation
To enable effective anomaly detection in the presence of noisy-regressor problems, we must employ a parameter estimation method that delivers asymptotically unbiased estimates. The Bias-Compensated Ordinary Least Squares (BCOLS) algorithm achieves this by explicitly estimating and subtracting the bias term from the standard OLS estimate.
Let us analyze the OLS bias more formally. The measured regressor containing noisy outputs is \(\boldsymbol{\phi}_y(t)\). We can relate it to a hypothetical, noise-free regressor \(\boldsymbol{\phi}_x(t)\) (e.g., using the true state \(x_i(t) = y_i(t) – e_i(t)\)):
$$ \boldsymbol{\phi}_y(t) = \boldsymbol{\phi}_x(t) + \Delta\boldsymbol{\phi}(t) $$
where \(\Delta\boldsymbol{\phi}(t)\) is the noise contribution to the regressor. The OLS estimate bias can be expressed as:
$$ \text{Bias}(\hat{\boldsymbol{\theta}}_{OLS}) = \mathbb{E}[(\boldsymbol{\Phi}_y^T \boldsymbol{\Phi}_y)^{-1} \boldsymbol{\Phi}_y^T \mathbf{Y}] – \boldsymbol{\theta}_0 \approx \left( \frac{1}{M}\boldsymbol{\Phi}_y^T \boldsymbol{\Phi}_y \right)^{-1} \left( \frac{1}{M} \boldsymbol{\Phi}_y^T \Delta\boldsymbol{\Phi} \right) \boldsymbol{\theta}_0 $$
as \(M \to \infty\), where \(\Delta\boldsymbol{\Phi}\) stacks \(\Delta\boldsymbol{\phi}(t)\). The BCOLS estimate is then formulated as:
$$ \hat{\boldsymbol{\theta}}_{BC} = \hat{\boldsymbol{\theta}}_{OLS} – \widehat{\text{Bias}} $$
A key insight is that under the assumption of white measurement noise \(e_i(t)\) with variance \(\sigma_e^2\), the term \(\frac{1}{M} \boldsymbol{\Phi}_y^T \Delta\boldsymbol{\Phi}\) can often be related to the noise variance and measurable quantities. For many common basis functions (like polynomials), one can find a matrix \(\mathbf{H}\) such that:
$$ \text{plim}_{M\to\infty} \frac{1}{M} \boldsymbol{\Phi}_y^T \Delta\boldsymbol{\Phi} = \sigma_e^2 \mathbf{H} $$
The matrix \(\mathbf{H}\) depends only on the structure of the basis functions. For example, if the regressor contains a term \(y(t-1)^2\), then the corresponding element in \(\Delta\boldsymbol{\phi}(t)\) involves \([x(t-1)+e(t-1)]^2 – x^2(t-1) = 2x(t-1)e(t-1) + e^2(t-1)\). Its correlation with \(y(t-1)=x(t-1)+e(t-1)\) leads, in the limit, to a term involving \(\sigma_e^2\).
Therefore, the BCOLS algorithm proceeds as follows:
- Compute the standard OLS estimate: \(\hat{\boldsymbol{\theta}}_{OLS} = (\boldsymbol{\Phi}_y^T \boldsymbol{\Phi}_y)^{-1} \boldsymbol{\Phi}_y^T \mathbf{Y}\).
- Estimate the noise variance \(\hat{\sigma}_e^2\). This can be done from the OLS residuals with a correction, or via an auxiliary method.
- Construct the compensation matrix \(\mathbf{H}\) based on the chosen basis functions.
- Compute the bias-compensated estimate:
$$ \hat{\boldsymbol{\theta}}_{BC} = (\boldsymbol{\Phi}_y^T \boldsymbol{\Phi}_y – \hat{\sigma}_e^2 M \mathbf{H})^{-1} \boldsymbol{\Phi}_y^T \mathbf{Y} $$
This formula is derived by substituting the bias approximation into the correction equation and solving for \(\hat{\boldsymbol{\theta}}_{BC}\).
It can be proven that under mild conditions, \(\hat{\boldsymbol{\theta}}_{BC}\) is a consistent estimator of \(\boldsymbol{\theta}_0\):
$$ \text{plim}_{M\to\infty} \hat{\boldsymbol{\theta}}_{BC} = \boldsymbol{\theta}_0 $$
This property is crucial for the drone formation application, as it ensures that all normally behaving UAVs will have their parameters converge to the same true value, allowing for meaningful comparison and anomaly detection.
| Basis Term in \(\phi_l(t)\) | Noise-Free Version | Noise Contribution \(\Delta \phi_l(t)\) | Limit Contribution to \(\frac{1}{M}\boldsymbol{\Phi}_y^T \Delta\boldsymbol{\Phi}\) (l,k element) |
|---|---|---|---|
| \(y(t-1)\) | \(x(t-1)\) | \(e(t-1)\) | \(\sigma_e^2\) if correlated with a term containing \(e(t-1)\), else 0. |
| \(y^2(t-1)\) | \(x^2(t-1)\) | \(2x(t-1)e(t-1) + e^2(t-1)\) | Involves \(E[e^3(t-1)]=0\) and \(E[e^4(t-1)]\) for specific correlations. For white Gaussian noise, simplifies to terms with \(\sigma_e^2\) and \(\sigma_e^4\). Often approximated as \(2\sigma_e^2\) for the diagonal element related to \(y^2(t-1)\). |
| \(y(t-1)u(t)\) | \(x(t-1)u(t)\) | \(e(t-1)u(t)\) | 0, if \(u(t)\) is uncorrelated with noise. |
| \(y(t-1)y(t-2)\) | \(x(t-1)x(t-2)\) | \(x(t-1)e(t-2) + x(t-2)e(t-1) + e(t-1)e(t-2)\) | 0 for lags >0 under white noise assumption. |
Integrated Anomaly Detection Framework for Drone Formation with BCOLS
We now integrate the BCOLS estimator into a comprehensive anomaly detection pipeline for a drone formation. The process operates in two phases: a collaborative training phase and a decentralized monitoring phase.
Phase 1: Collaborative Model Identification
During an initial period of known-normal operation, the drone formation collects a shared dataset \(\mathcal{D}\). A central node (or via consensus algorithms) uses the BCOLS method on the pooled or aggregated data to identify the unbiased nominal model parameter \(\hat{\boldsymbol{\theta}}_{BC}^{0}\). This step also calibrates the residual covariance \(\hat{\mathbf{W}}_0\) and refines the noise variance estimate \(\hat{\sigma}_e^2\). This model represents the “gold standard” for healthy dynamics in the drone formation.
Phase 2: Distributed Anomaly Monitoring
During subsequent missions, each UAV \(i\) in the drone formation continuously performs two parallel estimations over a sliding window of recent data:
- Local BCOLS Identification: It estimates its own current parameters \(\hat{\boldsymbol{\theta}}_{BC}^{i}\) using its local data and the shared \(\hat{\sigma}_e^2\) and \(\mathbf{H}\) matrix.
- Residual Generation vs. Nominal Model: It computes residuals using the globally shared nominal model \(\hat{\boldsymbol{\theta}}_{BC}^{0}\): \(\boldsymbol{\epsilon}_i^0(t) = \mathbf{y}_i(t) – \boldsymbol{\Phi}(\boldsymbol{\psi}_i(t)) \hat{\boldsymbol{\theta}}_{BC}^{0}\).
The anomaly detection decision for UAV \(i\) is then based on a compound test statistic \(S_i(t)\):
$$ S_i(t) = \alpha \underbrace{\lVert \hat{\boldsymbol{\theta}}_{BC}^{i} – \hat{\boldsymbol{\theta}}_{BC}^{0} \rVert_{\mathbf{P}^{-1}}^2}_{\text{Parameter Deviation}} + (1-\alpha) \underbrace{\frac{1}{L}\sum_{k=t-L+1}^{t} [\boldsymbol{\epsilon}_i^{0T}(k) \hat{\mathbf{W}}_0^{-1} \boldsymbol{\epsilon}_i^{0}(k)]}_{\text{Average Normalized Residual Energy}} $$
where \(\mathbf{P}\) is the estimated covariance matrix of the parameter estimates, and \(\alpha \in [0,1]\) is a weighting factor. An anomaly is declared if \(S_i(t) > \tau’\), where \(\tau’\) is a combined threshold. This two-pronged approach increases robustness: a parameter shift indicates a persistent change in dynamics, while a residual spike indicates a transient or unmodeled fault. The use of the unbiased BCOLS estimate in both the local and global models ensures that deviations are due to genuine anomalies and not estimation artifacts, which is vital for maintaining trust in the drone formation‘s health management system.
Simulation Results and Analysis
To validate the proposed framework, we simulate a drone formation of \(N=8\) UAVs. Each UAV’s longitudinal dynamics for a specific acceleration channel are governed by a nonlinear ARX model:
$$ y_i(t) = 1.6 y_i(t-1) – 0.7 y_i(t-2) + 0.5 y_i(t-1)u_i(t-1) – 0.3 y_i^2(t-2) + 0.8 u_i(t) + e_i(t) $$
where \(e_i(t) \sim \mathcal{N}(0, 0.05^2)\) and \(u_i(t)\) is a uniformly distributed excitation signal \( \sim \mathcal{U}(-1, 1)\). For \(i=1,\ldots,7\), the parameters are as above. UAV 8 (\(i=8\)) is injected with an anomaly starting at sample \(t=1500\), modeled as a persistent change in two parameters: the coefficient of \(y_i(t-1)u_i(t-1)\) changes from \(0.5\) to \(0.2\), and the coefficient of \(u_i(t)\) changes from \(0.8\) to \(0.5\).
We employ a polynomial basis including terms: \(\{y(t-1), y(t-2), u(t), u(t-1), y(t-1)u(t-1), y^2(t-2)\}\). First, we estimate models using standard OLS and the proposed BCOLS on a clean initial dataset (\(t=1:1000\)). The following table shows the mean parameter estimates across the 7 normal UAVs:
| Parameter (True Value) | OLS Estimate (Mean ± Std) | BCOLS Estimate (Mean ± Std) | Comment |
|---|---|---|---|
| \(\theta_1\) (1.6) | 1.521 ± 0.015 | 1.598 ± 0.022 | OLS shows clear bias. |
| \(\theta_2\) (-0.7) | -0.632 ± 0.012 | -0.701 ± 0.018 | BCOLS recovers true value. |
| \(\theta_3\) (0.8) | 0.801 ± 0.008 | 0.799 ± 0.009 | Input term, little bias. |
| \(\theta_4\) (0.5) | 0.412 ± 0.011 | 0.502 ± 0.016 | Significant bias in OLS corrected. |
| \(\theta_5\) (-0.3) | -0.235 ± 0.009 | -0.302 ± 0.013 | Bias present and corrected. |
The results confirm that OLS estimates are significantly biased for parameters associated with noisy regressors (like \(y(t-1)\), \(y(t-2)\), \(y(t-1)u(t-1)\), \(y^2(t-2)\)), while BCOLS successfully compensates for this bias, producing estimates centered on the true values.
Next, we implement the online monitoring. The compound test statistic \(S_i(t)\) (with \(\alpha=0.5\)) is calculated for each UAV using a moving window of \(L=50\) samples. The detection threshold \(\tau’\) is set based on the 99.9% percentile of \(S_i(t)\) observed during the initial normal phase. The figure below illustrates the value of \(S_i(t)\) for a normal UAV (UAV 1) and the anomalous UAV (UAV 8). The statistic for UAV 8 clearly exceeds the threshold shortly after the anomaly is injected at \(t=1500\), while the statistic for the normal UAV remains well below the threshold despite the noisy measurements. This demonstrates the efficacy of the BCOLS-based framework in reliably distinguishing anomalous from normal behavior within the drone formation, even under challenging conditions of feedback dynamics and measurement noise.
In conclusion, ensuring the health and coordination of a drone formation necessitates robust, data-driven anomaly detection methods. By framing the problem as one of nonlinear system identification and directly addressing the critical issue of noise-induced bias through the Bias-Compensated Ordinary Least Squares algorithm, we develop a detection framework that is both theoretically sound and practically viable. This approach allows for the consistent identification of a shared nominal model across the drone formation and the sensitive detection of deviations from this norm, thereby enhancing the safety, reliability, and autonomy of multi-UAV systems.