Cooperative ISAC-empowered Low-Altitude Economy

摘要—本文提出了一种用于低空感知场景的协作集成感知与通信（ISAC）方案，旨在通过协作估计无人机（UAVs）的参数并提升感知性能。所提出的方案包括两个阶段。第一阶段，我们通过使用张量分解模型来制定单点参数估计问题。通过利用因子矩阵的范德蒙德结构，引入了一个空间平滑张量分解方案来估计UAV的参数。为了进一步减少计算复杂度，我们设计了一个基于广义Rayleigh商（GRQ）的降维（RD）到达角（AoA）估计算法。第二阶段，通过多个基站（BS）之间的数据融合确定UAV的位置和真实速度。具体而言，我们首先开发了一种基于最小生成树（minimum spanning tree，MST）删除假数据的关联方法，以准确匹配基站的参数估计到同一个UAV。然后，提出了一种帕累托最优方法和残差加权方案，分别用于促进位置和速度估计。我们进一步扩展了我们的方法到双极化系统。仿真结果验证了所提出方案在与传统技术比较时的有效性。

关键词—集成感知与通信（ISAC），协作感知，张量分解，参数估计，数据融合。

I. INTRODUCTION

在低空感知场景中，根据发射机和接收机的空间部署，ISAC配置可以分为三种主要类别：monostatic [3], [4], bistatic [5], [6] and multistatic [7], [8]。具体而言，在双点和多点系统中，一个或多个基站（BS）充当发射机，其他基站充当接收机接收感知符号。然而，这些配置对多个基站的部署位置和同步性提出了较高的要求。不同的是，在单点配置中，基站同时发射感知符号，同时接收无人机（UAV）反射的回波信号，这简化了基站的部署，并提高了系统对不同步的容忍度。然而，传统的单点配置仍然面临一些挑战。首先，由于路径损耗，远距离UAV的感知性能较差[9]。此外，如果UAV与基站之间存在障碍物（如高楼或树木），检测概率将大幅下降。此外，UAV的综合状态（如真实速度）不能完全捕获，因为只能从基站的回波信号中获得径向速度（radial velocities）[7]。

为了克服传统单点配置的上述局限性，提出了协作ISAC框架[10]–[14]。在协作ISAC网络中，每个基站（BS）首先从接收的回波信号中估计目标的参数。然后，估计的参数在云端融合，以提高感知性能。在这种情况下，远离某些基站的目标可能相对靠近其他基站，从而增强ISAC网络的检测概率和覆盖范围。此外，基站的分布式部署进一步增强了感知性能，从而减少了估计误差，并增强了恢复更全面信息的能力，如UAV在三维（3-D）空间中的真实速度。最近在协作ISAC领域的研究旨在设计感知方案，以促进不同场景下的感知能力。在单天线配置中，[12]、[13]设计了在正交频分复用（OFDM）网络中的双阶段目标定位方案，其中目标的距离通过压缩感知（CS）和二维（2-D）快速傅里叶变换（FFT）算法分别得出，然后基于距离估计目标的位置。同样，[14]首先推导了目标的距离和径向速度，然后提出了一种符号级别的融合方案，通过多个基站协作促进位置和速度的估计。在多天线配置中，除了距离和速度外，基站还需要捕捉目标的到达角（AoA），以进一步增强感知性能，多个经典的AoA估计算法，如多重信号分类（MUSIC）算法，已被采用来推导目标的AoA[15]。

超越上述传统的估计方法，张量技术因其在提取多维参数方面的优越性而越来越受到关注，并已被应用于信道和目标参数估计。其原因在于，这些技术可以实现更高的估计精度，并且推导出的多维参数能够自动配对[16]–[19]。作者[19]对基于张量的ISAC进行了初步探索，其中信道和目标参数的估计是通过统一的张量方法实现的。然而，[19]并未揭示通过多基站合作提高感知精度的好处。

II. SYSTEM MODEL

A. Signal Model

B. Sensing Channel Model

III. PRELIMINARY STEPS FOR THE COOPERATIVE SENSING SCHEME

A. Beam Scanning and UAV Detection

需要注意的是，在UAV靠近基站（例如，少于一百米）的场景中，回波信号的信噪比（SNR）足以满足UAV检测和参数估计的要求。在这种情况下，可以直接通过几个参数估计方案实现UAV的精确定位，而无需额外的波束扫描或对准过程。然而，考虑到UAV可能远离基站的场景，以下的波束扫描和对准就显得至关重要，以抵消严重的路径损耗，从而提高SNR，促进后续的参数估计。

如图3(a)所示，我们首先通过波束扫描确定UAV的大致位置。为了进一步提高接收信号的SNR，我们将每个射频链上接收到的信号符号进行合并。通过这种方式，接收信号可以表示为：

为了检测UAV是否位于波束范围内，我们将二元假设检验（binary hypothesis testing，BHT）问题公式化为如下形式[23]：

$${\tilde{y}}_{m,n}=\{\begin{array}{l}{H}_0:{\tilde{m}}_{m,n},\\ H_1:f_{RX}^H{H}_{m,n}{f}_{TX}+{\tilde{m}}_{m,n},\end{array}$$

其中，零假设$H_0$假设基站（BS）接收到的只有噪声，而备择假设$H_1$则表明基站接收到了反射的回波信号和噪声。

为了求解上述BHT问题，我们首先需要构造一个检测器$T(\cdot )$，将${\tilde{y}}_{m,n}$映射为一个实数，然后将其与预定义的阈值$\gamma$进行比较，以确定是否接受$H_0$或$H_1$ [23]，即：

$$T({\tilde{y}}_{m,n})\underset{H_0}{\stackrel{H_1}{\ge }}\gamma .$$

通过考虑特定的场景和可用的先验知识，可以利用各种假设检验方法，例如似然比检验（LRT）来开发检测器[24]。为了进一步处理，我们引入了检测标志$d_{i,j}$，具体来说，如果UAVs在索引$(i,j)$范围内被检测到，则$d_{i,j}=1$；否则，$d_{i,j}=0$。

B. Estimation of the Number of UAVs

对于检测到UAV存在的波束范围（即，$d_{i,j}=1$），需要进一步估计波束范围内UAV的数量，这可以通过考虑信息理论（information-theoretic criteria）标准[4]来实现。具体来说，我们首先估计接收信号${\tilde{y}}_{m,n}$的协方差，表示为：

$$\hat{R}=\frac{1}{MN}\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}{\tilde{y}}_{m,n}{\tilde{y}}_{m,n}^H\in {\mathbb{C}}^{R\times R}.$$

然后，我们对$\hat{R}$执行奇异值分解（SVD），即，$\hat{R}=U\Lambda U^H$，其中$\Lambda =D({\lambda }_1,\dots ,{\lambda }_R)$是一个对角矩阵，特征值按降序排列，即${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_R$。随后，通过采用最小描述长度（MDL）准则[25]，波束范围内UAV的数量估计为：

Ki,j=arg⁡min⁡k∈{1,…,R−1}MDL(k),K_{i,j} = \arg \min_{k \in \{1, \dots, R-1\}} \text{MDL}(k), Ki,j​=argk∈{1,…,R−1}min​MDL(k),

${}^1$ 需要注意的是（Noting），我们假设当前波束范围外的无人机（UAV）将不会被当前波束的旁瓣检测到。这是（This）因为旁瓣功率远小于主瓣功率。此外（Additionally），我们还可以通过[4]中提出的剪枝过程进一步消除相邻波束内同一无人机的冗余检测。

在上述步骤之后，可以通过记录在{di,j}i=0,j=0NH−1,NV−1\{d_{i,j}\}_{i=0,j=0}^{N_H-1,N_V-1}{di,j​}i=0,j=0NH​−1,NV​−1​中等于1的元素的索引来确定UAV在三维空间中的大致位置，并通过对每个波束范围内UAV的估计数量求和来推导出UAV的总数，公式如下：

$$K=\sum _{i=0,j=0}^{N_H-1,N_V-1}{d}_{i,j}{K}_{i,j},$$

这些步骤为后续的单点参数估计提供了有价值的先验信息和指导。如图3(b)所示，我们可以将整个天线阵列分为$N^{\prime }$组，其中每组$di,j$与方向$({\theta }_i,{\varphi }_j)$对齐。举例来说，我们可以假设数字预编码矩阵为单位矩阵，而模拟预编码矩阵中的非零元素填充为对应于方向$({\theta }_i,{\varphi }_j)$的归一化引导向量。然而，组合矩阵的非零元素则从归一化单位圆中选择，以缓解HBF结构中由于维度减少而导致的到达角（AoA）估计歧义。

需要注意的是，前述的与波束相关的步骤已经被广泛研究，且并非本文的重点，在接下来的部分中，我们假设所需的先验信息，即波束检测标志{di,j}i=1,j=1NH,NV\{d_{i,j}\}_{i=1,j=1}^{N_H,N_V}{di,j​}i=1,j=1NH​,NV​​和UAV的总数$K$，已经通过上述步骤准确推导出来。

IV. TENSOR DECOMPOSITION APPROACH FOR PARAMETER ESTIMATION

A. Tensor Preliminaries

为了增强本文的可读性，提供了一些关于张量的基本理论和关键定义。我们建议读者参考[26]以获取更多详细信息。

1. 展开：张量$\mathcal{X}\in {\mathbb{C}}^{I_1\times I_2\times \cdots \times I_N}$的模式-$n$展开表示将张量重新排列为一个矩阵。具体来说，张量元素$(i_1,i_2,\dots ,i_N)$映射到矩阵元素$(i_n,j)$，其中
   $$j=1+\sum _{k=1,k\ne n}^N(i_k-1)J_k,\text{且}{J}_k=\prod _{m=1,m\ne n}^{k-1}{I}_m.$$ (19)

2. 秩-1张量：一个$N$阶张量$\mathcal{X}$是一个秩-1张量，如果它可以表示为$N$个向量的外积，即：
   $$\mathcal{X}=a^{(1)}\circ a^{(2)}\circ \cdots \circ a^{(N)}.$$ (20)

3. CP分解：CANDECOMP/PARAFAC（CP）分解表示将张量分解为一组秩-1张量的和，即：
   $$\mathcal{X}=\sum _{r=1}^R{\lambda }_r{a}_r^{(1)}\circ a_r^{(2)}\circ \cdots \circ a_r^{(N)},$$ (21)
   其中$R$表示$\mathcal{X}$的秩。对应的$n$-阶模式的因子矩阵定义为
   $$A^{(n)}=[a_1^{(n)},\dots ,a_R^{(n)}],n=1,\dots ,N.$$
   $\mathcal{X}$的模式-$n$展开形式为：
   $$X^{(n)}=A^{(n)}\Lambda {\left(A^{(N)}\circ \cdots \circ A^{(n+1)}\circ A^{(n-1)}\circ \cdots \circ A^{(1)}\right)}^T,$$ (22)
   其中$\Lambda =D(\left[{\lambda }_1,\dots ,{\lambda }_R\right])$。

B. Tensor Formulation

B. Position Estimation

在精确的数据关联之后，我们接着在多个基站 (BS) 之间执行数据融合。通常，有两种主要的数据融合策略：硬融合和软融合 [7]。具体来说，前者指的是每个 BS 估计并将最终的定位结果，即 {p^k,j}k=1,j=1K,J\{\hat{\mathbf{p}}_{k,j}\}_{k=1,j=1}^{K,J}{p^​k,j​}k=1,j=1K,J​，发送到云端。然后，云端通过一个 (加权) 平均融合方案来整合这些位置估计。例如，依赖于平均融合的位置估计由下式给出

$${\hat{\mathbf{p}}}_k=\frac{1}{J}\sum _{j=1}^J{\hat{\mathbf{p}}}_{k,j},k=1,\dots ,K,\text{\hspace{1em}(72)}$$

不同的是，对于软融合策略，所有的 BS 将原始的参数估计，包括到达角 (AoAs) 和距离，发送到云端。然后，云端通过整合上传的参数来执行对无人机 (UAV) 位置的估计。

${}^2$假设 BS${}_j$ 成功检测到第 $k$ 个无人机，那么总是存在至少另一个 BS，其对同一无人机的定位结果与 BS${}_j$ 的估计相似，这保证了 $W_{\min }(a)\le \varsigma$，其中 $a$ 表示对应于 BS${}_j$ 估计的顶点。相反， $W_{\min }(a)>\varsigma$ 表示 BS${}_j$ 未能检测到第 $k$ 个无人机。因此，为防止错误检测影响后续的数据融合，顶点 $a$ 应在更新的图 $\tilde{G}$ 中被识别为一个孤立顶点。

由于软融合为数据融合提供了更细粒度的信息，我们在本文的后续内容中主要关注这种融合策略。特别地，注意到位置估计的精度主要取决于距离和径向方向 (由 AoAs 确定) 的估计，我们分别定义了以下的距离损失函数和方向损失函数

$$f_r({\mathbf{p}}_k)=\frac{{\sum }_{j=1}^J{\alpha }_{k,j}|d_{k,j}({\mathbf{p}}_k)-{\hat{d}}_{k,j}|}{{\sum }_{j=1}^J{\alpha }_{k,j}},\text{\hspace{1em}(73a)}$$$$f_d({\mathbf{p}}_k)=\frac{{\sum }_{j=1}^J{\alpha }_{k,j}\Vert {\mathbf{r}}_{k,j}({\mathbf{p}}_k)-{\hat{\mathbf{r}}}_{k,j}{\Vert }_2}{{\sum }_{j=1}^J{\alpha }_{k,j}},\text{\hspace{1em}(73b)}$$

其中 $d_{k,j}({\mathbf{p}}_k)=\Vert {\mathbf{p}}_k-{\mathbf{p}}_j^B{\Vert }_2$ 表示无人机和 BS 之间的真实距离，而 ${\mathbf{r}}_{k,j}({\mathbf{p}}_k)$ 表示真实径向方向向量。符号 ${\alpha }_{k,j}$ 表示分配给 BS${}_j$ 的权重。特别地，由于路径损耗，回波信号的信噪比 (SNR) 随着无人机和 BS 之间距离的增加而减小。因此，我们定义加权系数为

$${\alpha }_{k,j}=\frac{1}{(d_{k,j}({\mathbf{p}}_k){)}^{{\beta }_1}},\text{\hspace{1em}(74)}$$

其中 ${\beta }_1>0$ 是一个为控制加权强度而引入的因子。然后，我们解决以下问题来估计第 $k$ 个无人机的位置，即即

min⁡pk{fr(pk),fd(pk)},(75)\min_{\mathbf{p}_k} \{f_r(\mathbf{p}_k), f_d(\mathbf{p}_k)\}, \tag{75} pk​min​{fr​(pk​),fd​(pk​)},(75)

这是一个高度非线性的多目标优化问题，其目标函数 (OFs) 具有不同的维度。为了解决问题 (75)，我们提出一种帕累托最优策略（Pareto optimality strategy）来确定无人机的位置。具体来说，如图 6 所示，我们首先将在 (72) 中的平均融合结果作为粗略的位置估计。然后，我们围绕这个粗略估计构建 $L$ 个格点（lattice points），并计算对应于这些格点的两个 OF 值。为了同时最小化这两个 OF 值，我们首先识别出格点中的被支配解（dominated solutions），即

U~k={pkl1∣fr(pkl1)>fr(pkl2),fd(pkl1)>fd(pkl2),pkl1,pkl2∈Uk},(76)\tilde{U}_k=\{\mathbf{p}_k^{l_1}|f_r(\mathbf{p}_k^{l_1})>f_r(\mathbf{p}_k^{l_2}), f_d(\mathbf{p}_k^{l_1})>f_d(\mathbf{p}_k^{l_2}), \mathbf{p}_k^{l_1}, \mathbf{p}_k^{l_2} \in U_k \}, \tag{76} U~k​={pkl1​​∣fr​(pkl1​​)>fr​(pkl2​​),fd​(pkl1​​)>fd​(pkl2​​),pkl1​​,pkl2​​∈Uk​},(76)

其中 Uk={pk1,…,pkL}U_k = \{\mathbf{p}_k^1, \dots, \mathbf{p}_k^L\}Uk​={pk1​,…,pkL​} 表示由 $L$ 个格点形成的集合。然后，我们移除被支配解并保留解的帕累托集为

$${\bar{U}}_k=U_k\setminus {\tilde{U}}_k.\text{\hspace{1em}(77)}$$

最后，我们根据系统配置从帕累托集中确定解作为最终的位置估计。具体来说，在每个 BS 被分配了大量子载波但只配备了相对少量天线的场景中，距离估计总是比 AoA 估计表现出更高的精度。因此，距离估计在定位中应该更具主导性。在这种情况下，我们可以从帕累托集中选择具有最小距离损失函数值的解，即

$${\hat{\mathbf{p}}}_k=\arg \underset{{\mathbf{p}}_k^l\in {\bar{U}}_k}{\min }{f}_r({\mathbf{p}}_k^l).\text{\hspace{1em}(78)}$$

另一方面（On the flip side），当 BS 配备了大量天线但子载波数较少时，AoA 估计成为关键因素，引导我们从帕累托集中选择具有最小方向损失函数值的解，即

$${\hat{\mathbf{p}}}_k=\arg \underset{{\mathbf{p}}_k^l\in {\bar{U}}_k}{\min }{f}_d({\mathbf{p}}_k^l).\text{\hspace{1em}(79)}$$

此外，我们也可以在设计的单基地参数估计算法中考虑 AoA 和距离估计的精度，以实现更好的定位性能。注意到计算对应于格点的 OF 值只需要几次简单的计算，因此复杂度是可以接受的。