微型无人机阵列位置扰动下的DoA估计混合Cramér-Rao界

Hybrid Cramér-Rao Bound for DoA Estimation with Micro-UAV Array Under Positional Errors

摘要

随着无人机技术的快速发展，基于微型无人机集群的分布式天线阵列因其灵活部署、低成本、可重构等优势，在无线通信、目标探测、环境监测等领域展现出广阔的应用前景。然而，在实际部署过程中，受外界环境干扰和无人机自身控制精度限制，无人机集群难以精确到达预设位置，从而产生随机的位置扰动误差。这些位置扰动不仅会恶化阵列的空间频谱特性，还会显著降低对目标信号的探测、定位和参数估计性能，成为制约无人机集群阵列实用化的关键问题之一。针对上述挑战，基于位置扰动误差在相邻无人机间传播的物理特性，建立了考虑位置扰动的无人机集群阵列信号模型，重点研究了位置扰动对目标信号到达角（DoA）估计性能的影响。通过引入位置扰动相关矩阵，推导了DoA估计在有限采样情况下的混合克拉美罗下界（HCRB），构建了评估位置扰动影响的理论基础。为了验证理论分析的有效性，本文通过数值仿真，从位置扰动方差、信噪比、无人机数量，以及DoA角度等多个角度，对比分析了有无位置扰动情况下的DoA估计性能。仿真结果表明：位置扰动显著恶化了DoA估计精度，且随着信噪比的提高，DoA估计对位置误差更加敏感。此外，增加无人机数量可在一定程度上降低HCRB，但难以消除位置扰动带来的估计误差。

Abstract

Objective In recent years, applications of unmanned aerial vehicles (UAVs) have experienced rapid growth in both military and civilian domains, including aerial surveillance, remote sensing, and search and rescue. Beyond individual use, UAV swarms are widely employed as micro-UAV array radars. High-precision detection is achieved through beamforming by treating a swarm as a large movable antenna aperture. Therefore, high-precision flight control becomes essential. However, outdoor formations often rely on relative position-based strategies to avoid GPS inaccuracies. Although this approach remains theoretically feasible, micro-UAVs can deviate because of turbulence, and positional errors can be inherited by subsequent UAVs, leading to significant degradation in detection performance. Therefore, this study aims to evaluate the impact of such inherent positional perturbations on the direction of arrival (DoA) estimation performance. Methods This study focused on the lower bound of estimation accuracy to assess the efficiency and precision of estimators in the presence of calibration errors. Because the analytical calculation of the Cramér‒Rao bound (CRB) with random nuisance parameters was challenging, the hybrid Cramér‒Rao bound (HCRB) was adopted, as it represented an effective balance between computability and tightness. Based on the assumption of the inheritability of UAV array perturbations, the HCRB was derived for the DoA of a single unknown source. A swarm of N UAVs was considered as a uniform linear array, where each UAV carried an omnidirectional antenna. The model jointly analyzed perturbations in the parallel and perpendicular directions to the aperture under a far-field narrowband signal scenario. In addition, numerical simulations and Monte Carlo experiments based on the matched filter (MF) estimator were conducted to validate the theoretical derivation. Results and Discussions Numerical simulations analyzed the HCRB performance under various influencing factors. The results demonstrated that as the positional perturbation variance increases, the HCRB for both low and high SNR cases increases from the CRB level without interference, confirming that the error-free CRB represents the lower bound of the HCRB. Monte Carlo experiments exhibited trends consistent with these theoretical results. In addition, as the SNR increases, the HCRB of DoA becomes more sensitive to positional errors. Although the HCRB decreases with higher SNR, smaller positional variance consistently produces lower HCRB values. Regarding the number of UAVs, although the HCRB monotonically decreases as N increases from 1 to 128, it did not converge to the CRB without perturbations. Finally, the impact of different DoA angles indicated that the HCRB is non-uniform and axially symmetric, achieving its lowest value at 0° and increasing as the absolute value of the DoA increases. Conclusions The inherent nature of positional errors arising from relative position-based formations, which has been overlooked in previous literature, significantly affects UAV array performance. The derivation and quantitative analysis of the HCRB in this study confirm that the inheritability of perturbations leads to a performance floor that cannot be eliminated simply by increasing the SNR or the number of UAVs. These findings provide an important theoretical basis for assessing efficiency and predicting the performance of mobile array radars in practical flight environments, particularly for applications such as interferometric synthetic aperture radar (InSAR).

Graphical abstract

关键词

无人机 / 阵列 / 位置扰动 / 到达角 / Cramér-Rao界

Key words

UAV / array / positional errors / DoA / Cramér-Rao bound

本刊网刊

近年来，随着无人机技术的飞速发展，无人机集群在军、民用领域受到了广泛关注，应用包括空中监视、遥感和搜救行动等^[1‒2]。相比单个无人机，无人机集群能够克服单个无人机的局限性，实现更复杂的任务。无人机集群中的多个无人机通过协同工作，可以组成一个虚拟的大规模天线阵列，从而提高系统的性能。编队无人机阵列已广泛用于民用领域，如现场控制或监督、测绘等^[3]。然而，由于无人机的随机分布和动态移动，无人机集群面临着诸多挑战，其中，到达角（DoA）估计是一个重要的研究问题。本文对无人机集群到达角估计的最新研究进展进行综述，重点讨论无人机位置扰动、信道估计、多目标跟踪、重构智能表面等热点问题。

Zeng等^[1]设计了一种集成感知和通信的分布式协同跟踪框架，通过无人机间资源协调，在降低通信开销的同时提高了感知和跟踪性能。Rihan等^[2]分析了位置扰动对DoA估计性能的影响，理论和仿真结果表明位置扰动导致空间谱性能恶化，存在峰值功率和主瓣宽度的折中，且随无人机数量的增加，性能恶化加剧。不同DoA估计方法的归一化空间谱渐近收敛到相同形状。Bahiki等^[3]先利用辅助用户估计无人机取向角，再通过求根法估计DoA，该方法与无人机间隔无关。Zhou^[4]、Du^[5]等通过实验室和外场试验，验证了无人机集群DoA估计的有效性，重点讨论了传感器误差和位置扰动等因素的影响。Fan等^[6]将信道分解为DoA、位置、增益等多个参数，并使用秩约简估计器估计DoA。Chen等^[7]提出一种多阵型跟踪方法，通过阵型变换实现全向估计、高增益跟踪和包围。Chen等^[8]综述了射频检测定位技术，分析了面临的诸多挑战，提出了一些潜在的解决方案，如融合射频与其他传感器，并引入人工智能等。Fan等^[9]提出一种重构智能表面辅助的无人机集群到达角估计新方法，该系统仅需一个接收通道，定义了含位置扰动参数的原子范数，大大降低了复杂度，改善了估计性能。

微型无人机阵列雷达是研究中非常重要的一项技术。具体来说，就是将微型无人机阵列视为空中一个大的可移动天线阵列孔径，通过微型无人机阵列的波束成形和信号处理技术来实现窄波束、高精度的目标探测^[4‒5]。然而，无人机集群到达角估计与传统的阵列信号处理有所不同，主要面临的挑战之一是阵元位置扰动。由于无人机的随机分布和动态移动，阵元位置很难准确知晓，存在随机扰动，位置扰动会恶化DoA估计性能。

由此可见，对无人机编队进行高精度的飞行控制是非常重要的环节。到目前为止，室外无人机阵列的编队控制主要有两种技术，即基于GPS位置的部署和基于相对位置的部署^[6‒9]。由于GPS位置数据不准确，误差在2～5 m之间，可能会引发无人机间的碰撞。因此，大多数设计师和研究人员采用基于相对位置的控制策略。

尽管基于相对位置的编队在理论上是可行的，但在实践中仍然存在一些问题。由于气流扰动，飞行中的微型无人机可能会偏离其预定位置。此外，在基于相对位置编队的情况下，一架无人机的位置误差会被其后续无人机继承，这可能会大大降低移动阵列雷达的探测性能。

在许多情况下，天线阵列可以有效地通过联合时间和空间维度来检索表征信号源的参数。尽管已有许多研究分析了位置扰动对无人机或天线阵列性能的影响^[10‒12]，但对于相对位置编队引起的位置误差仍具有继承特性这一特点，大多数文献都没有给予充分考虑。针对此类含有随机扰动的复杂估计问题，文献[13‒15]对混合Cramér‒Rao下界（HCRB）的计算方法进行了深入讨论，为后续的性能评估提供理论计算基础。

研究表明，传感器微小的位置误差都会严重降低传感器阵列的性能^[16‒19]。因此，在存在校准误差的情况下，估计精度的下限有助于判断估计器的效率、预测性能及精度，特别是在干涉合成孔径雷达（InSAR）领域^[20]。众所周知，Cramér‒Rao下界对于一大类估计器来说是严格的。

然而，在有随机干扰参数的情况下，很难通过数值积分分析计算得到下界。为此，Jiang等^[21]修正了Cramér‒Rao下界（CRB）被应用于雷达问题中的参数估计。在修正方法中，HCRB被证明代表了其所指估计器类的易计算性、紧密性和宽度（即估计无偏的条件）之间的良好权衡。不过，研究也观察到，除非满足某些特定条件，否则HCRB通常不是可达界，甚至不是渐近可达的^[22]。

在无人机编队的实际应用中，这种继承性位置误差的建模与分析至关重要。Gini^[23]探讨了修正Cramér-Rao下界在非高斯杂波参数估计中的雷达应用，为处理复杂背景下的扰动提供了思路。针对基于相对位置编队的微型无人机阵列，Feng等^[24]深入分析了位置扰动对阵列构型及通信性能的具体影响。此外，Ferreol等^[25]考虑到实际探测环境中，位置误差会直接导致经典子空间算法（如MUSIC算法）在建模误差存在时的渐近性能显著退化。

为了进一步提升系统的鲁棒性，针对机载传感器阵列在流形扰动下的被动源定位问题^[26]、卫星相控阵天线机械失真的在线校准方法^[27]均被提出。在MIMO雷达成像等领域，利用最大后验估计方法进行稀疏自校准也成为了应对位置不确定性的重要手段^[28]。在性能评估的理论框架下，Gini等^[29]详细阐述了存在随机干扰参数时类Cramér-Rao界限的使用准则；Gini^[23]进一步通过修正下界量化了非高斯环境对估计精度的限制。

因此，通过整合这些理论，本文基于无人机阵列扰动的可继承性这一假设，推导出单个未知源DoA的HCRB，从而得出位置误差对DoA估计性能影响的定量分析。为简单起见，假设无人机阵列的编队为均匀线性阵列，在扰动继承的假设下，对孔径的平行和垂直方向的扰动进行联合分析。此外，讨论了位置方差、信噪比（SNR）、无人机数量和到达角（DoA）等因素对DoA的HCRB性能的影响；深入剖析继承性扰动如何通过改变阵列几何形态，最终影响无人机阵列的DoA估计性能。

1 DoA的HCRB推导

考虑无人机阵列有N个无人机，每个无人机都携带一根全向天线。假设第n个无人机可能会随机偏离其预定的位置，在平行于孔径的方向上偏移d_n，在垂直于孔径的方向上偏移r_n，考虑远场情况下的窄带信号到达角为\theta，微型无人机阵列扰动示意图如图1所示。

阵列天线接收到的信号\mathit{x}在时间t时刻可以表示为：

\begin{array}{rr}\mathit{x}(t)& =\overline{\mathit{a}}(\theta )s(t)+\mathit{n}(t)\end{array}（1）

式中：信号波形s(t)和加性白噪声\mathit{n}(t)都服从复高斯分布，即s(t)~\mathit{𝒞}\mathit{𝒩}(\mathrm{0},{\sigma }_s^{\mathrm{2}})，\mathit{n}(t)~\mathit{𝒞}\mathit{𝒩}(\mathrm{0},{\sigma }_n^{\mathrm{2}}{\mathit{I}}_N)，其中，{\mathit{I}}_N为N\times N单位矩阵，{\sigma }_s^{\mathrm{2}}和{\sigma }_n^{\mathrm{2}}分别为输入信号和每个天线接收噪声的方差；\overline{\mathit{a}}为扰动的导向向量。

由于微型无人机阵列的编队是基于相对位置关系进行部署的，这意味着位置误差是可继承的，即d_n=d_{n-\mathrm{1}}+u_n，r_n=r_{n-\mathrm{1}}+v_n，其中，u_n和v_n分别为第n架无人机的平行和垂直位置噪声的扰动，假设它们都服从均值为0且方差为{\sigma }_d^{\mathrm{2}}的高斯分布，即u_n,v_n ~\mathit{𝒩}(\mathrm{0},{\sigma }_d^{\mathrm{2}})，则第n个孔径的平行和垂直位置误差满足：

\left\{\begin{array}{l}{d}_n~\mathit{𝒩}(\mathrm{0},n{\sigma }_d^{\mathrm{2}}),\\ r_n~\mathit{𝒩}(\mathrm{0},n{\sigma }_d^{\mathrm{2}})\end{array}\right.（2）

接收信号的协方差矩阵可以表示为：

\mathit{C}(\mathit{\xi })={\mathbb{E}}_{\mathit{x}}[\mathit{x}(t){\mathit{x}}^{\mathrm{H}}(t)]={\sigma }_s^{\mathrm{2}}\overline{\mathit{a}}(\theta ){\overline{\mathit{a}}}^{\mathrm{H}}(\theta )+{\sigma }_n^{\mathrm{2}}{\mathit{I}}_N（3）

式中：{\mathbb{E}}_{\mathit{x}}(\cdot )代表对\mathit{x}的期望；\mathit{\xi }为未知参数向量，\mathit{\xi }=[\theta ,d_{\mathrm{0}},d_{\mathrm{1}},\cdots ,d_{N-\mathrm{1}},r_{\mathrm{0}},r_{\mathrm{1}},\cdots ,r_{N-\mathrm{1}}{]}^{\mathrm{T}}；上标H表示共轭转置。关于参数的“无偏性”已经在文献[15]中进行了讨论，本文重点讨论参数的均方误差的性能。

关于未知参数向量\mathit{\xi }的条件似然函数为：

\begin{array}{c}p(\mathit{x}|\mathit{\xi })=\frac{\mathrm{e}\mathrm{x}\mathrm{p}(-{\mathit{x}}^{\mathrm{H}}{\mathit{C}}^{-\mathrm{1}}(\mathit{\xi })\mathit{x})}{{\mathrm{\pi }}^N|\mathit{C}(\mathit{\xi })|}\end{array}（4）

需要注意的是，DoA没有先验信息，位置误差假设为独立的高斯分布。因此，先通过构造\mathit{\rho }=[d_{\mathrm{0}},d_{\mathrm{1}},\cdots ,d_{N-\mathrm{1}},r_{\mathrm{0}},r_{\mathrm{1}},\cdots ,r_{N-\mathrm{1}}{]}^{\mathrm{T}}，得出\mathit{\rho }的联合先验概率密度函数为：

p(\mathit{\rho })=\frac{\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\rho }}^{\mathrm{H}}{\mathit{\Sigma }}^{-\mathrm{1}}\mathit{\rho })}{{(\mathrm{2}\mathrm{\pi })}^N|\mathit{\Sigma }{|}^{\mathrm{1}/\mathrm{2}}}（5）

式中，\mathit{\Sigma }为整体协方差矩阵，具体计算式如下：

\mathit{\Sigma }=\left[\begin{array}{cc}{\mathit{\Sigma }}_d& \\ & {\mathit{\Sigma }}_r\end{array}\right]（6）

\left\{\begin{array}{l}{\mathit{\Sigma }}_d={\sigma }_d^{\mathrm{2}}{\displaystyle \sum _{n=\mathrm{0}}^{N-\mathrm{1}}}{\mathit{c}}_n{\mathit{c}}_n^{\mathrm{T}},\\ {\mathit{\Sigma }}_r={\mathit{\Sigma }}_d\end{array}\right.（7）

式（6）、（7）中：{\mathit{\Sigma }}_d为轴向扰动协方差；{\mathit{\Sigma }}_r为垂直轴向的协方差；{\mathit{c}}_n=[{\mathit{o}}_n;{\mathit{l}}_{N-n}]，其中，{\mathit{o}}_n为长度为n的全0向量，{\mathit{l}}_{N-n}为长度为N-n的全1向量，分号代表列的分隔符。

然后，参数向量\mathit{\xi }的混合信息矩阵可以表示为^[13]：

\mathit{J}={\mathit{J}}_{\mathit{x}|\mathit{\xi }}+{\mathit{J}}_{\mathit{\xi }}（8）

式中，{\mathit{J}}_{\mathit{x}|\mathit{\xi }}和{\mathit{J}}_{\mathit{\xi }}分别为观测和先验统计的信息矩阵。矩阵{\mathit{J}}_{\mathit{x}|\mathit{\xi }}的第(n,m)个元素，可以表示为：

{J_{\mathit{x}|\mathit{\xi }}}_{,n,m}=-{\mathbb{E}}_{\mathit{x},\mathit{\xi }}[\frac{{\partial }^{\mathrm{2}}\mathrm{l}\mathrm{n}p(\mathit{x}|\mathit{\xi })}{\partial {\mathit{\xi }}_n\partial {\mathit{\xi }}_m}]（9）

式中，n,m=\mathrm{0,1},\cdots ,\mathrm{2}N，{\mathbb{E}}_{\mathit{x},\mathit{\xi }}为关于\mathit{x}和\mathit{\xi }的联合期望，{\mathit{\xi }}_n、{\mathit{\xi }}_m分布为\mathit{\xi }的第n、m项。

对于\mathit{x}|\mathit{\xi }，其观测信息矩阵的第(n,m)个元素可以进一步表示为^[14]：

J_{\mathit{x}|\mathit{\xi },n,m}=-{\mathbb{E}}_{\mathit{\xi }}[\frac{{\partial }^{\mathrm{2}}}{\partial {\xi }_n\partial {\xi }_m}{\mathbb{E}}_{\mathit{x}}\mathrm{l}\mathrm{n}p(\mathit{\xi }|\mathit{x})]={\mathbb{E}}_{\mathit{\xi }}\{\mathrm{T}\mathrm{r}[{\mathit{C}}^{-\mathrm{1}}(\mathit{\xi })\frac{\partial \mathit{C}(\mathit{\xi })}{\partial {\xi }_n}{\mathit{C}}^{-\mathrm{1}}(\mathit{\xi })\frac{\partial \mathit{C}(\mathit{\xi })}{\partial {\xi }_m}]\}（10）

式中，Tr表示矩阵的迹。

为了方便推导，在下文中用\mathit{C}表示\mathit{C}(\mathit{\xi })。则\mathit{C}关于\theta、d_n和r_n的导数分别为：

{\dot{\mathit{C}}}_{\theta }={\sigma }_s^{\mathrm{2}}(\dot{\overline{\mathit{a}}}{\overline{\mathit{a}}}^{\mathrm{H}}+\overline{\mathit{a}}{\dot{\overline{\mathit{a}}}}^{\mathrm{H}})（11）

{\dot{\mathit{C}}}_{dn}=\mathrm{j}{k}_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\theta {\sigma }_s^{\mathrm{2}}({\overline{a}}_n{\mathit{e}}_n{\overline{\mathit{a}}}^{\mathrm{H}}-\overline{\mathit{a}}{\mathit{e}}_n^{\mathrm{H}}{\overline{a}}_n^{*})（12）

{\dot{\mathit{C}}}_{rn}=\mathrm{j}{k}_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\theta {\sigma }_s^{\mathrm{2}}({\overline{a}}_n{\mathit{e}}_n{\overline{\mathit{a}}}^{\mathrm{H}}-\overline{\mathit{a}}{\mathit{e}}_n^{\mathrm{H}}{\overline{a}}_n^{*})（13）

式（11）～（13）中：\dot{\overline{\mathit{a}}}为\overline{\mathit{a}}相对于\theta的导数；{\overline{a}}_n为\overline{\mathit{a}}的第n个元素，{\overline{a}}_n=\mathrm{e}\mathrm{x}\mathrm{p}\{\mathrm{j}{k}_{\mathrm{0}}((nd+d_n)\mathrm{s}\mathrm{i}\mathrm{n}\theta +r_n\mathrm{c}\mathrm{o}\mathrm{s}\theta )\}，k_{\mathrm{0}}=\mathrm{2}\mathrm{\pi }/{\lambda }_{\mathrm{0}}，其中，{\lambda }_{\mathrm{0}}为信号波长，d为无人机间距{\mathit{e}}_n为第n列N维单位矩阵；上标*表示共轭。

注意到

{\dot{\overline{a}}}_n=\mathrm{j}{k}_{\mathrm{0}}[(nd+d_n)\mathrm{c}\mathrm{o}\mathrm{s}\theta -r_n\mathrm{s}\mathrm{i}\mathrm{n}\theta ]\cdot {\overline{a}}_n=b_n{\overline{a}}_n（14）

式中，b_n\triangleq \mathrm{j}{k}_{\mathrm{0}}[(nd+d_n)\mathrm{c}\mathrm{o}\mathrm{s}\theta -r_n\mathrm{s}\mathrm{i}\mathrm{n}\theta ]。

为便于表示，令

A\triangleq {\dot{\overline{\mathit{a}}}}^{\mathrm{H}}\dot{\overline{\mathit{a}}}={\displaystyle \sum _{n=\mathrm{0}}^{N-\mathrm{1}}}|b_n{|}^{\mathrm{2}}（15）

B\triangleq {\overline{\mathit{a}}}^{\mathrm{H}}\dot{\overline{\mathit{a}}}={\displaystyle \sum _{n=\mathrm{0}}^{N-\mathrm{1}}}{b}_n（16）

式（15）、（16）中，{\dot{\overline{\mathit{a}}}}^{\mathrm{H}}\overline{a}=-{\overline{\mathit{a}}}^{\mathrm{H}}\dot{\overline{\mathit{a}}}。

通过一些简单的推导，可以得到A和B^{\mathrm{2}}对\mathit{\xi }的期望如下：

\begin{array}{l}{\mathbb{E}}_{\mathit{\xi }}[A]=\frac{\mathrm{1}}{\mathrm{2}}N(N-\mathrm{1})k_{\mathrm{0}}^{\mathrm{2}}{\sigma }_d^{\mathrm{2}}+\\ \frac{\mathrm{1}}{\mathrm{6}}N(N-\mathrm{1})(\mathrm{2}N-\mathrm{1})k_{\mathrm{0}}^{\mathrm{2}}{d}^{\mathrm{2}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\theta \end{array}（17）

{\mathbb{E}}_{\mathit{\xi }}[B^{\mathrm{2}}]=-\frac{\mathrm{1}}{\mathrm{2}}N(N-\mathrm{1})k_{\mathrm{0}}^{\mathrm{2}}{\sigma }_d^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{4}}{N}^{\mathrm{2}}{(N-\mathrm{1})}^{\mathrm{2}}{k}_{\mathrm{0}}^{\mathrm{2}}{d}^{\mathrm{2}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\theta（18）

根据Sherman-Morrison公式，\mathit{C}的逆矩阵可以计算为：

{\mathit{C}}^{-\mathrm{1}}={\sigma }_n^{-\mathrm{2}}({\mathit{I}}_N-\frac{\mathrm{1}}{N_{\mathrm{0}}}\overline{\mathit{a}}{\overline{\mathit{a}}}^{\mathrm{H}})（19）

式中，N_{\mathrm{0}}=\frac{\mathrm{1}}{e_{\mathrm{0}}}+N，e_{\mathrm{0}}\triangleq {\sigma }_s^{\mathrm{2}}/{\sigma }_n^{\mathrm{2}}，则有：

{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{\theta }=e_{\mathrm{0}}(\frac{\overline{\mathit{a}}{\dot{\overline{\mathit{a}}}}^{\mathrm{H}}}{N_{\mathrm{0}}{e}_{\mathrm{0}}}+\dot{\overline{\mathit{a}}}{\overline{\mathit{a}}}^{\mathrm{H}}-\frac{\overline{\mathit{a}}{\overline{\mathit{a}}}^{\mathrm{H}}\dot{\overline{\mathit{a}}}{\overline{\mathit{a}}}^{\mathrm{H}}}{N_{\mathrm{0}}})（20）

{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dn}=\mathrm{j}{k}_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\theta e_{\mathrm{0}}({\overline{a}}_n{\mathit{e}}_n{\overline{\mathit{a}}}^{\mathrm{H}}-\frac{\overline{\mathit{a}}{\overline{\mathit{a}}}^{\mathrm{H}}}{N_{\mathrm{0}}}-\frac{{\overline{a}}_n^{*}\overline{\mathit{a}}{\mathit{e}}_n^{\mathrm{H}}}{N_{\mathrm{0}}{e}_{\mathrm{0}}})（21）

由上，{\mathit{J}}_{\mathit{x}|\mathit{\xi }}的前部分元素可以推导如下：

\begin{array}{l}{\mathit{J}}_{\theta \theta }\triangleq {\mathbb{E}}_{\mathit{\xi }}\mathrm{T}\mathrm{r}[{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{\theta }{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{\theta }]=\\ \frac{e_{\mathrm{0}}{k}_{\mathrm{0}}^{\mathrm{2}}}{N_{\mathrm{0}}}(N{(N-\mathrm{1})}^{\mathrm{2}}{\sigma }_d^{\mathrm{2}}+\frac{N^{\mathrm{2}}(N^{\mathrm{2}}-\mathrm{1})}{\mathrm{6}}{d}^{\mathrm{2}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\theta )\end{array}（22）

{\mathit{J}}_{\theta dn}\triangleq {\mathbb{E}}_{\mathit{\xi }}\mathrm{T}\mathrm{r}[{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{\theta }{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dn}]=\frac{N(\mathrm{2}n+\mathrm{1}-N)}{\mathrm{2}{N}_{\mathrm{0}}}{e}_{\mathrm{0}}{k}_{\mathrm{0}}^{\mathrm{2}}d\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}\theta（23）

\begin{array}{l}{\mathit{J}}_{dd}\triangleq [{\mathbb{E}}_{\mathit{\xi }}\mathrm{T}\mathrm{r}({\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dn}{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dm}{)]}_{N,N}=\\ \frac{\mathrm{2}}{N_{\mathrm{0}}}{e}_{\mathrm{0}}{k}_{\mathrm{0}}^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta (N{\mathit{I}}_N-\mathit{l}{\mathit{l}}^{\mathrm{T}})\end{array}（24）

式中，\mathit{l}为全1向量。

注意到{\dot{\mathit{C}}}_{dn}和{\dot{\mathit{C}}}_{rn}除了标量系数外其余几乎相同，因此，可类比推导出{\mathit{J}}_{\mathit{x}|\mathit{\xi }}。为此，先定义C_s\triangleq \mathrm{c}\mathrm{o}\mathrm{s}\theta /\mathrm{s}\mathrm{i}\mathrm{n}\theta，然后有：

{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{rn}=C_s{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dn}（25）

因此，{\mathit{J}}_{\mathit{x}|\mathit{\xi }}后部分元素可以分别表示为：

{\mathit{J}}_{\theta rn}=C_s{\mathbb{E}}_{\mathit{\xi }}\mathrm{T}\mathrm{r}[{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{\theta }{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dn}]（26）

{\mathit{J}}_{dnrm}=C_s{\mathbb{E}}_{\mathit{\xi }}\mathrm{T}\mathrm{r}[{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dn}{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dm}]（27）

{\mathit{J}}_{rnrm}=C_s^{\mathrm{2}}{\mathbb{E}}_{\mathit{\xi }}\mathrm{T}\mathrm{r}[{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dn}{\mathit{C}}^{-\mathrm{1}}{\dot{\mathit{C}}}_{dm}]（28）

由于{\mathit{J}}_{\mathit{x}|\mathit{\xi }}是Hermitian矩阵，完整的{\mathit{J}}_{\mathit{x}|\mathit{\xi }}可以表示为：

{\mathit{J}}_{\mathit{x}|\mathit{\xi }}=\left[\begin{array}{ccc}{\mathit{J}}_{\theta \theta }& {\mathit{J}}_{\theta d}& {\mathit{J}}_{\theta r}\\ {\mathit{J}}_{\theta d}^{\mathrm{H}}& {\mathit{J}}_{dd}& {\mathit{J}}_{dr}\\ {\mathit{J}}_{\theta r}^{\mathrm{H}}& {\mathit{J}}_{dr}^{\mathrm{H}}& {\mathit{J}}_{rr}\end{array}\right]（29）

对于{\mathit{J}}_{\mathit{\xi }}的先验信息矩阵，第(n,m)个元素可以表示为：

{J_{\xi }}_{,nm}=-{\mathbb{E}}_{\mathit{\xi }}[\frac{{\partial }^{\mathrm{2}}\mathrm{l}\mathrm{n}p(\mathit{\rho })}{\partial {\mathit{\xi }}_n\partial {\mathit{\xi }}_m}],n,m=\mathrm{0,1},\cdots ,\mathrm{2}N（30）

根据式（2）、（5），矩阵{\mathit{J}}_{\mathit{\xi }}可以简化为：

{\mathit{J}}_{\mathit{\xi }}=\left[\begin{array}{cc}\mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathit{\Sigma }}^{-\mathrm{1}}\end{array}\right]（31）

最后，\mathit{\xi }的第i个元素的估计量{\widehat{\mathit{\xi }}}_i的HCRB满足：

\mathrm{v}\mathrm{a}\mathrm{r}({\widehat{\mathit{\xi }}}_i)\ge [({\mathit{J}}_{\mathit{x}|\mathit{\xi }}+{\mathit{J}}_{\mathit{\xi }}{)}^{-\mathrm{1}}{]}_{ii},i=\mathrm{0,1},\cdots ,\mathrm{2}N（32）

式中，\mathrm{v}\mathrm{a}\mathrm{r}({\widehat{\mathit{\xi }}}_i)为{\mathit{\xi }}_i的无偏估计量的方差。具体来说，对于DoA的HCRB，有：

\mathrm{v}\mathrm{a}\mathrm{r}(\widehat{\theta })\ge \mathrm{H}\mathrm{C}\mathrm{R}\mathrm{B}(\theta )=\frac{\mathrm{1}}{J_{\theta \theta }-{\mathit{q}}^{\mathrm{T}}(\overline{\mathit{J}}+{\mathit{\Sigma }}^{-\mathrm{1}}{)}^{-\mathrm{1}}\mathit{q}}（33）

式中，

\mathit{q}=[{\mathit{J}}_{\theta d},{\mathit{J}}_{\theta r}{]}^{\mathrm{T}}（34）

扰动信息矩阵为：

\overline{\mathit{J}}=\left[\begin{array}{cc}{\mathit{J}}_{dd}& {\mathit{J}}_{dr}\\ {\mathit{J}}_{dr}^{\mathrm{H}}& {\mathit{J}}_{rr}\end{array}\right]（35）

注意：通过设置{\sigma }_d^{\mathrm{2}}=\mathrm{0}，DoA的HCRB退化为与没有扰动的DoA的CRB相同的表达，可以表示为：

\mathrm{C}\mathrm{R}\mathrm{B}(\theta )=\frac{\mathrm{6}{N}_{\mathrm{0}}}{N^{\mathrm{2}}(N^{\mathrm{2}}-\mathrm{1})e_{\mathrm{0}}{k}_{\mathrm{0}}^{\mathrm{2}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\theta d^{\mathrm{2}}}（36）

2 仿真结果与分析

通过数值仿真，分析位置扰动对无人机阵列DoA估计性能的影响，包括位置方差{\sigma }_d^{\mathrm{2}}、SNR、无人机数量和DoA等参数。仿真中，考虑一个由10架无人机组成的均匀线性阵列，相邻无人机之间的距离设置为电磁信号波长的一半，即半波长距离。信号波长选取为1 m，对应频率为300 MHz。这种阵列构型和参数设置能够在保证无人机编队紧凑性的同时，避免阵元间的电磁耦合效应。

在仿真中，假设位置扰动服从0均值高斯分布，标准差在0到1倍波长之间变化。这一设置综合考虑了实际无人机集群中的定位误差水平和模型适用范围。此外，信噪比（SNR）分别设置为-10、0 dB，以分析不同噪声环境下位置扰动的影响。

实验主要分析位置扰动对无人机阵列中DoA估计性能的影响。实验内容包括：1）位置扰动方差对HCRB的影响，以无位置误差的CRB为基准，调节位置扰动方差，分别对SNR为-10 dB和0 dB的情况进行数值仿真，使用匹配滤波（MF）估计器作为角度估计器的参照，每个参数点进行10 000次蒙特卡洛仿真。2）比较不同SNR情况下的HCRB，随着SNR的增加，DoA的HCRB对位置误差会变得更加敏感。以无位置误差的CRB为基准，在两种位置扰动方差下，对SNR从-10 dB到10 dB进行数值分析。3）考察无人机数量对HCRB的影响，在不同位置扰动方差下，分析无人机数量对HCRB的影响，假设SNR为0 dB，无人机数量从1架增加到128架，观察HCRB的变化趋势。4）分析不同DoA对HCRB的影响，在默认配置下，对不同的DoA进行分析，同时考虑不同位置扰动方差的情况。

首先，考虑位置扰动方差{\sigma }_d^{\mathrm{2}}对HCRB的影响。以没有位置误差的CRB为基准，从\mathrm{0}到\mathrm{1}对扰动方差进行调节，分别在SNR为-10、0 dB两种情况，以最为经典的匹配滤波估计器（MF）为对角度估计器协方差下限的参照，每个点进行10 000次的蒙特卡洛仿真，具体结果如图2所示。由图2可知：随着{\sigma }_d^{\mathrm{2}}的增加，\mathrm{S}\mathrm{N}\mathrm{R}\mathrm{为}-\mathrm{10}dB和\mathrm{S}\mathrm{N}\mathrm{R}\mathrm{为}\mathrm{0}dB两种情况下的HCRB都从无干扰的CRB水平往上增加，这印证了没有位置误差的CRB是有位置错误的HCRB的下限。基于MF的蒙特卡洛实验在走势上与理论值一致，进一步印证了结论推导的正确性。同时，基于MF的角度估计的蒙特卡洛实验，其协方差在数值上与理论下限有近一个数量级的差距（其中，部分误差来自角度-功率谱网格的划分误差），说明MF估计器不是最佳估计器。

然后，比较不同SNR对HCRB的影响。随着SNR的增加，DoA的HCRB对位置误差会变得更加敏感。为了进一步阐明不同SNR对DoA的HCRB的影响，以没有位置误差时的CRB为基准，在{\sigma }_d^{\mathrm{2}}分别为0.1、1.0两种情况下，对SNR从-\mathrm{10} dB到\mathrm{10} dB进行了数值分析，具体结果如图3所示。

从图3可以看出：DoA的HCRB随着SNR的增加而减小；{\sigma }_d^{\mathrm{2}}越小，DoA的HCRB越小，即{\sigma }_d^{\mathrm{2}}\mathrm{为}\mathrm{0.1}的HCRB在不同的信噪比下均小于{\sigma }_d^{\mathrm{2}}\mathrm{为}\mathrm{1}.0的值，符合预期。

此外，还考察了HCRB在不同方差值下，即{\sigma }_d^{\mathrm{2}}\mathrm{分别}\mathrm{为}\mathrm{0.1}和\mathrm{1.0}时，无人机数量N对HCRB的影响，为了便于观察，假设SNR为0 dB，具体结果如图4所示。由图4可知：随着无人机数量从1架增加到128架，HCRB呈单调下降，但有位置误差的HCRB无法收敛到没有扰动情况下的CRB；另外，位置误差的方差越小，其HCRB的值也越小，与预期相符。

最后，仿真了不同DoA对HCRB的影响，也对不同的位置方差的情况进行了分析。延续上述默认的配置，假设SNR为0 dB，分析不同的DoA对HCRB的影响，具体结果如图5所示。

从图5中可以看到：不同DoA所对应的HCRB是不均匀的，即\theta =\mathrm{0}°的信号在相同{\sigma }_d^{\mathrm{2}}下达到最低的HCRB；随着DoA绝对值的增加，HCRB也增加。此外，HCRB关于\theta =\mathrm{0}°呈轴对称，这与没有位置误差时的CRB的趋势一致。

3 结 论

本文研究了基于相对位置部署的无人机阵列DoA估计，并重点分析了无人机阵元相对于孔径的水平和垂直方向的位置扰动误差所带来的影响。在此基础上，推导出了无人机阵列DoA估计的HCRB，得到了位置误差对DoA估计影响的定量分析。通过数值仿真讨论了位置误差的方差、SNR、无人机数量和DoA等因素对HCRB的影响。理论和仿真结果均表明，位置误差的方差对HCRB影响较大，通过增加无人机数量可以一定程度地对冲该影响，但无法收敛到无干扰的情况。本研究仅考虑单源一维均匀线阵，且扰动假设为独立同分布且忽略气动耦合，未来研究可扩展至三维非均匀阵与多源场景，引入实时扰动估计与补偿闭环。

参考文献

[1] Zeng Yong, Zhang Rui, Lim T J.Wireless communications with unmanned aerial vehicles:Opportunities and challenges[J].IEEE Communications Magazine,2016,54(5):36-42. doi:10.1109/mcom.2016.7470933

[2] Rihan M, Selim M M, Xu Chen,et al.D2D communication underlaying UAV on multiple bands in disaster area:Stochastic geometry analysis[J].IEEE Access,2019,7:156646-156658. doi:10.1109/access.2019.2949813

[3] Bahiki M R, Talib N A, Azrad S.Relative positioning-based system with tau control for collision avoidance in swarming application[J].Materials Science and Engineering,2016,152:012025. doi:10.1088/1757-899x/152/1/012025

[4] Zhou Longyu, Leng Supeng, Wang Qing,et al.Integrated sensing and communication in UAV swarms for cooperative multiple targets tracking[J].IEEE Transactions on Mobile Computing,2023,22(11):6526-6542. doi:10.1109/TMC.2022.3193499

[5] Du Jiaxin, Guo Wenbo, Yan Mu,et al.Location error analysis for collaborative beamforming in UAVs random array[J].IEEE Wireless Communications Letters,2024,13(3):904-907. doi:10.1109/lwc.2024.3350163

[6] Fan Dian, Guo Gan, Song Jiaming,et al.DOA estimation for arbitrarily distributed subarrays in UAV swarm[C]//Proceedings of the 2020 IEEE/CIC International Conference on Communications in China(ICCC).Chongqing:IEEE,2020:794-799. doi:10.1109/ICCC49849.2020.9238865

[7] Chen Zhong, Liu Wenbo, Huff G H.Experimental campaign to evaluate the fundamental capabilities and limitations of synthetic DOA using swarming UAVs[C]//Proceedings of the 2021 15th European Conference on Antennas and Propagation(EuCAP).Dusseldorf:IEEE,2021:1-4. doi:10.23919/eucap51087.2021.9411152

[8] Chen Zhong, Yeh S, Chamberland J F,et al.Impact of position errors on synthetic aperture DOA convergence based on swarming UAVs[C]//Proceedings of the 2020 IEEE USNC-CNC-URSI North American Radio Science Meeting(Joint with AP-S Symposium).Montreal:IEEE,2021:83-84. doi:10.23919/usnc/ursi49741.2020.9321600

[9] Fan Dian, Gao Feifei, Ai Bo,et al.Channel estimation and self-positioning for UAV swarm[J].IEEE Transactions on Communications,2019,67(11):7994-8007. doi:10.1109/tcomm.2019.2933203

[10] Keskin F, Filik T.Isotropic and directional DOA estimation of the target by UAV swarm-based 3-D antenna array[C]//Proceedings of the 2020 4th International Symposium on Multidisciplinary Studies and Innovative Technologies (ISMSIT).Istanbul:IEEE,2020:1-7. doi:10.1109/ismsit50672.2020.9254341

[11] Fratu O, Nor A M, Halunga S,et al.RF propagation and interferences challenges for UAVs swarm detection and positioning[C]//Proceedings of the 2022 25th International Symposium on Wireless Personal Multimedia Communications(WPMC).Herning:IEEE,2023:501-504. doi:10.1109/wpmc55625.2022.10014750

[12] Chen Peng, Chen Zhimin, Zheng Beixiong,et al.Efficient DOA estimation method for reconfigurable intelligent surfaces aided UAV swarm[J].IEEE Transactions on Signal Processing,2022,70:743-755. doi:10.1109/tsp.2022.3146791

[13] Vincent J A, Arnold E J.Beamforming sensitivity of airborne distributed arrays to flight tracking and vehicle dynamics[C]//Proceedings of the 2017 IEEE Aerospace Conference.Big Sky:IEEE,2017:1-14. doi:10.1109/aero.2017.7943588

[14] Leuschen C, Hale R, Keshmiri S,et al.UAS-based radar sounding of the polar ice sheets[J].IEEE Geoscience and Remote Sensing Magazine,2014,2(1):8-17. doi:10.1109/mgrs.2014.2306353

[15] Hauert S, Leven S, Varga M,et al.Reynolds flocking in reality with fixed-wing robots:Communication range vs.maximum turning rate[C]//Proceedings of the 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems.San Francisco:IEEE,2011:5015-5020. doi:10.1109/iros.2011.6095129

[16] Wilson D B, Göktoğan A H, Sukkarieh S.A vision based relative navigation framework for formation flight[C]//Proceedings of the 2014 IEEE International Conference on Robotics and Automation(ICRA).Hong Kong:IEEE,2014:4988-4995. doi:10.1109/icra.2014.6907590

[17] Meng Wei, He Zhirong, Su Rong,et al.Decentralized multi-UAV flight autonomy for moving convoys search and track[J].IEEE Transactions on Control Systems Technology,2017,25(4):1480-1487. doi:10.1109/tcst.2016.2601287

[18] Strader J, Gu Yu, Gross J N,et al.Cooperative relative localization for moving UAVs with single link range measurements[C]//Proceedings of the 2016 IEEE/ION Position,Location and Navigation Symposium(PLANS).Savannah:IEEE,2016:336-343. doi:10.1109/plans.2016.7479718

[19] Namin F, Petko J S, Werner D H.Analysis and design optimization of robust aperiodic micro-UAV swarm-based antenna arrays[J].IEEE Transactions on Antennas and Propagation,2012,60(5):2295-2308. doi:10.1109/tap.2012.2189715

[20] Petko J S, Werner D H.Positional tolerance analysis and error correction of micro-UAV swarm based antenna arrays[C]//Proceedings of the 2009 IEEE Antennas and Propagation Society International Symposium.North Ch-arleston:IEEE,2009:1-4. doi:10.1109/aps.2009.5172135

[21] Jiang Wei, Haimovich A M.Cramér-Rao bound for noncoherent direction of arrival estimation in the presence of sensor location errors[J].IEEE Signal Processing Letters,2017,24(9):1303-1307. doi:10.1109/lsp.2017.2716944

[22] Van Trees H L, Bell K L.Bayesian bounds for parameter estimation and nonlinear filtering/tracking[M].New York:Wiley,2007. doi:10.1109/9780470544198

[23] Gini F.A radar application of a modified Cramér-Rao bound:Parameter estimation in non-Gaussian clutter[J].IEEE Tr-ansactions on Signal Processing,1998,46(7):1945-1953. doi:10.1109/78.700966

[24] Feng F Y P, Rihan M, Huang Lei.Positional perturbations analysis for micro-UAV array with relative position-based formation[J].IEEE Communications Letters,2021,25(9):2918-2922. doi:10.1109/lcomm.2021.3086676

[25] Ferreol A, Larzabal P, Viberg M.On the asymptotic performance analysis of subspace DOA estimation in the presence of modeling errors:Case of MUSIC[J].IEEE Transactions on Signal Processing,2006,54(3):907-920. doi:10.1109/tsp.2005.861798

[26] Mir H S, Sahr J D, Hatke G F,et al.Passive source localization using an airborne sensor array in the presence of manifold perturbations[J].IEEE Transactions on Signal Processing,2007,55(6):2486-2496. doi:10.1109/tsp.2007.893936

[27] Takahashi T, Nakamoto N, Ohtsuka M,et al.On-board calibration methods for mechanical distortions of satellite phased array antennas[J].IEEE Transactions on Antennas and Propagation,2012,60(3):1362-1372. doi:10.1109/tap.2011.2180303

[28] Liu Changchang, Yan Jin, Chen Weidong.Sparse self-calibration by map method for MIMO radar imaging[C]//Proceedings of the 2012 IEEE International Conference on Acoustics,Speech and Signal Processing(ICASSP).Kyoto:IEEE,2012:2469-2472. doi:10.1109/icassp.2012.6288416

[29] Gini F, Reggiannini R.On the use of Cramér-Rao-like bounds in the presence of random nuisance parameters[J].IEEE Transactions on Communications,2000,48(12):2120-2126. doi:10.1109/26.891222

基金资助

深圳市高等院校稳定支持计划(20220818232202001)

深圳市科技计划项目(GJHZ20220913143013024)