作为世界上最大的单口径射电望远镜，500 m口径球面射电望远镜(five-hundred-meter aperture spherical radio telescope, FAST)具有极高的灵敏度，在接收天体信号的同时，不可避免地会受到自然或人工的其他无线电信号干扰。得益于FAST电磁宁静区的建立和运行，以及台址内部电子设备的电磁兼容措施，望远镜周围的电磁环境维持在稳定良好的水平^[1-2]。但随着周边经济发展和基础建设的不断推进，无线电信号干扰也随之增加，需要对FAST周边的无线电环境进行长期稳定的监测，为电磁环境保护和无线电频谱管理工作提供依据。新一代的干扰监测系统需要对干扰信号具备高灵敏度的探测和识别能力，并能对干扰源进行高精度定位。

由于FAST台址地处喀斯特地貌山区，地形复杂，山体对信号存在遮挡、多径等效应，单个监测站对干扰信号的探测能力不足^[3]，而且周边环境中信号种类繁多，不能单纯地检测信号是否存在。

可以根据干扰程度对识别出的信号进行排序，并结合定位系统对主要的干扰信号源进行定位。目前已有不少基于深度学习进行信号调制方式识别的研究，基本思路是将时域信号转换为星座图，利用图像识别的方式进行干扰信息源定位^[4]。但这种方法忽略了信号的大量有用信息，如时间或频域上的相关性等。传统的信号源定位采用到达时间(time of arrival, TOA)法，基于信号到达不同接收机所用的时间，得到信号发射点与各接收机之间的距离，然后以各接收机为圆心，以测得的距离为半径作圆，各圆的交点即为被测信号源所在的位置。由于在复杂环境定位中，非视距传播、电磁波的反射及绕射等引起的多径效应会使监测站接收到延迟的信号，故而大大降低定位精度^[5]。由此可知，TOA定位算法需要精确的时间检测装置，对硬件设备要求较高。因此，本文采用多站点的协同谱感知技术，实现在低信噪比下的干扰信号探测^[6-7]；采用到达时间差(time difference of arrival, TDOA)法实现高精度的干扰源定位；并将信号的瀑布图作为网络的输入，基于深度神经网络强大的学习能力，实现不同制式信号的识别。

1 干扰信号探测识别技术 1.1 多站点协同频谱感知

基于多站点的频谱感知首先要考虑如何快速传输数据，本文采用基于压缩感知理论的频谱感知技术，利用信号的稀疏性，通过找到合适的基矩阵，实现对信号的采样和压缩，以不低于Nyquist采样率的速率获得更精确的信号。

多站点数据融合的处理方法可以分为非相干融合和相干融合(coherent fusion)2类。非相干融合的实现较为简单，适用于信噪比较高的场景，其中常见的一种处理方法是基于能量检测的频谱感知法^{[3, 8]}，即通过计算信号在特定时间内的能量，与设定的门限值进行比较，作出判决；相干融合是利用多站点信号之间的相关性，分析出信号间的相位和时延差别，基于此对信号校准后通过相关函数的尖锐程度来判断是否有信号存在^[9-12]，实现起来相对复杂，但更适应于低信噪比的场景。本文采用相干融合的方法判断信号是否存在。

假设站点i的接收信号序列为

$x_i(n)=h_i s_i(n)+n_i(n) .$

(1)

其中：h_i为站点i接收到的信号的衰减，n为时间序列，s_i(n)为信号，n_i(n)为噪声。对x_i(n)进行离散Fourier变换(discrete Fourier transform, DFT)，得到

$X_i(\omega)=\operatorname{DFT}\left(x_i(n)\right).$

(2)

其中：X_i(ω)为站点i的信号频谱，ω为对应的角频率。通过谱相关和离散Fourier反变换(inverse discrete Fourier transform, IDFT)可以求解出站点i和站点j信号间的相关函数R_ij(τ)，可表示为

$R_{i j}(\tau)=\operatorname{IDFT}\left(X_i(\omega) X_j^*(\omega)\right) .$

(3)

其中：τ为时间量，X_j(ω)为站点j的信号频谱，*为其共轭。设计了一种相关函数叠加公式来实现站点i和站点j接收到的信号的相关叠加, 具体表示为：

$R(\tau)=\sum\limits_i \sum\limits_j R_{i j}\left(\tau-\Delta \tau_{j i}\right) \rho_{i j}^*,$

(4)

$\rho_{i j}=\underset{R_{i j}(\tau)}{\arg \max }\left(\left|R_{i j}(\tau)\right|\right).$

(5)

其中：R(τ)为相关函数叠加的结果，Δτ_ji为站点j和站点i的TDOA，ρ_ij^*为幅度加权系数ρ_ij的共轭。

Δτ_ji实现了时延校准，ρ_ij实现了相位校准，ρ_ij的模值表示融合的功率比，得到的融合互相关函数R(τ)可以作为信号存在的依据。如果信号存在，则R(τ)存在明显的峰值，图 1给出了有信号和无信号2种条件下有无相关系数的仿真结果。可以看出，相比于图 1a和1b中信号强度的差异，图 1c和1d中具有相关系数的R(τ)之间的差异更为突出，故利用后者作为判别信号是否存在的依据更为可靠。因此，本文利用相关函数的峰均比R_PAP作为信号是否存在的依据，可表示为

$R_{\mathrm{PAP}}=\frac{\max \left(|R(\tau)|^2\right)}{|\overline{R(\tau)}|^2} .$

(6)

R_PAP值越大，包含峰值的概率越高，最终可以得到融合判决公式，表示为

$U= \begin{cases}1, & R_{\mathrm{PAP}}>T_{\mathrm{cf}} ; \\ 0, & R_{\mathrm{PAP}}<T_{\mathrm{cf}} .\end{cases}$

(7)

其中：T_cf为相干融合判决门限，当R_PAP＞T_cf，U=1时，表示探测到信号；当R_PAP＜T_cf，U=0时，表示没有探测到信号。

1.2 基于深度神经网络的信号识别

为了充分利用信号在时间和频域上的相关性，本文采用信号的瀑布图作为网络的输入，具体做法为令S_t={S_{t, 1}, S_{t, 2}, …, S_{t, N}}表示t时刻待测频段所接收到的频谱序列，包含相位和幅度响应，因此S_{t, n}为复数，总的频率通道数目n=1, 2, …, N。对应的瀑布图表示如下：

$\boldsymbol{S}=\left[\begin{array}{c} \boldsymbol{s}_1 \\ \boldsymbol{s}_2 \\ \vdots \\ \boldsymbol{s}_T \end{array}\right]=\left[\begin{array}{cccc} s_{1, 1} & s_{1, 2} & \cdots & s_{1, N} \\ s_{2, 1} & s_{2, 2} & \cdots & s_{2, N} \\ \vdots & \vdots & & \vdots \\ s_{T, 1} & s_{T, 2} & \cdots & s_{T, N} \end{array}\right].$

(8)

其中：行向量代表了信号频域上的相关性，列向量代表了时间上的相关性。由于多个站点的接收机均可以计算出瀑布图，因此需要先将多个瀑布图有效融合，可采用概率最大比合并方式，也可采用重构源信号的方式融合多个瀑布图。本文采用的是重构源信号的方式，将在2.1节中详细介绍此方式的算法。

由于卷积神经网络(convolutional neural network，CNN)无法处理复数的输入，故S不能直接作为网络输入；如果将复数的实部和虚部分开，又忽略了信号时间和频域间的相关性。所以，需要对传统的CNN进行复数化改造，主要是对前端的卷积层进行改造，具体方式是使用复数作为卷积核中的系数，乘法运算改为复数乘法运算，同时激活函数的输入是复数的模值，网络的输出可以是实数，也可以是激活前同相位的复数，这主要取决于后续的模块采用的是常规还是复数模块。基于复数值的特征表达有以下几个特点：更易于优化、更具泛化特性、更易于学习，以及更利于提升噪声鲁棒的记忆机制。

噪声干扰会降低信号的感知质量和可识别性，这正是本文所需要解决的信号感知的关键问题。为解决此问题，本文采用了包含复数卷积层的CNN(C-CNN)的网络，结构如图 2所示。其中复数卷积层可对信号进行二维滤波处理，全连接层则可对提取后的特征进行分类。

信号识别神经网络，即C-CNN网络的具体设计如表 1所示。网络输入大小为256×4，使用4层复数卷积层从输入的频谱信息中提取高级特征。其中复数卷积网络使用全0填充，以更好地提取边缘信息。神经网络的步长为2，使每层输出的特征图的长度减半，但是每层卷积核的个数加倍，保证每层复数卷积网络的特征容纳量保持不变。为了提高模型的泛化能力，并减少过拟合，使用L2正则化和Dropout。Dropout是指在网络训练过程中，对于神经网络单元，按照一定的概率将其暂时从网络中丢弃，从而有效防止过拟合。在C-CNN模型中根据预测结果调整L2正则化的权重为0.001，Dropout的权重为0.5。激活函数使用ReLU，最后一层使用Sigmoid函数进行归一化。将监测站收集到的频谱数据传入神经网络，包括内时刻、中心频率、带宽和功率，网络的输出为干扰信号的类型。

2 干扰源定位技术

TDOA是一种双曲线定位方法，首先计算目标辐射源发射的信号到达各基站的时间差，之后以各基站为焦点，以时间差对应的距离差为长轴作双曲线，不同双曲线相交，交点即目标辐射源所在的位置^[13]。TDOA的定位精度主要受时差的计算精度、信号源和各接收系统的相对几何位置(布站方式)等因素的影响^[14-15]。为了获得高精度的时差估计，本文提出了一种基于广义互相关^[16]和重构源信号的迭代算法，提高时差的计算精度。

2.1 高精度时差估计

为了准确估计到达时间差，最重要的是确保源信号的准确性。因此首先从多个接收机的信号中将抵达波的时延和衰减估计出来，再重构源信号。在源信号满足误差要求之后，再用其估计各接收机信号的到达时间差，算法的实现流程如图 3所示。

(1) 选取功率或信噪比最大的一路信号作为源信号，利用广义相关法依次对各站点信号进行时差求解。设s(t)为目标源发出的随时间t变化的信号，2个具有一定距离的监测站(监测站1和2)收到的目标源发出的信号s₁(t)和s₂(t)以及信号x₁(t)和x₂(t)分别为

$\begin{aligned} &s_1(t)=A_1 s\left(t-\tau_1\right), \\ &s_2(t)=A_2 s\left(t-\tau_2\right), \\ &x_1(t)=A_1 s\left(t-\tau_1\right)+n_1(t), \\ &x_2(t)=A_2 s\left(t-\tau_2\right)+n_2(t) . \end{aligned}$

(9)

其中：A₁和A₂分别为监测站1和2对干扰信号的增益，τ₁和τ₂为监测站1和2接收到的信号的传输时延，n₁(t)和n₂(t)分别为监测站1和2的噪声。假定s₁(t)、s₂(t)、n₁(t)和n₂(t)是实信号，是联合平稳的、零均值的随机过程，并且相互之间是统计独立的，则可以得到监测站1和2接收到的信号的相关函数为

$R_{x_1 x_2}(\tau)=R_{s_1 s_2}(\tau)+R_{s_1 n_2}(\tau)+R_{s_2 n_1}(\tau).$

(10)

其中：R_x1x2(τ)为接收到总信号的相关函数，R_s1s2(τ)为2个监测站接收到的目标源信号的相关函数；R_s2n1(τ)和R_s1n2(τ)为目标源信号与监测站噪声之间的相关函数。由于噪声和信号、噪声和噪声之间互不相关，式(10)可简化为

$R_{x_1 x_2}(\tau)=R_{s_1 s_2}(\tau)=A_1 A_2 R_s(\tau-D) .$

(11)

其中：D=τ₂－τ₁为所求时间差；当τ=D时，相关函数取得最大值，D即为所求时间差。

(2) 求解误差函数。当误差足够小并满足要求时，将所得高精度时差输出；当误差较大时，对信号的相位、时差进行纠正，重构源信号并重复之前的操作。

具体求解过程为假定同一时刻站点i收到的信号为x_i(t)，估计的源信号为$ \hat x\left( t \right)$，对应接收机噪声为n_i(t)，源信号估计误差为n(t)，则

$\begin{gathered} \hat{x}(t)=s(t)+n(t), \\ x_i(t)=d_i s(t-\Delta \tau)+n_i(t) . \end{gathered}$

(12)

其中：Δτ为通过广义相关法求解的信号到达监测站的时间差；s(t-Δτ)为监测站点接收到的目标源信号；d_i为2个信号的相对衰减，包含幅度和相位的相对关系。假设n_i(t)和n(t)互不相关，则误差函数为

$L_i\left(d_i\right)=\int\left[\frac{x_i(t+\Delta \tau)}{d_i}-\hat{x}(t)\right]^2 \mathrm{~d} t .$

(13)

其中：L_i为信号估计的误差函数，x_i(t+Δτ)为监测站i接收到的总信号。当d_i的估计最准时，误差函数取最小值，即

$\hat{d}_i=\underset{d_i}{\operatorname{argmin}} L_i\left(d_i\right) .$

(14)

之后，根据估计的$ {\hat d_i}$作为相位和衰减的校准依据，对源信号进行重构，

$\hat{x}(t)=\sum\limits_i \frac{x_i(t+\Delta \tau)}{\hat{d}_i} .$

(15)

经反复迭代，直到误差函数的最小值满足条件为止，最终根据重构的源信号，得到各监测站信号到达时间差的高精度估计。

2.2 基于Taylor级数展开的信号定位技术

在获得信号到达时间差的高精度估计之后，就可以换算为干扰源到不同监测站之间的距离差，进而建立干扰源位置的定位方程，求解该方程即可得到干扰源的估算位置。

假定干扰源在Cartesian坐标系下的坐标为(x, y, z)，站点i的坐标为(x_i, y_i, z_i)，参考监测站的坐标为(x₀, y₀, z₀)。监测站i和参考监测站之间的距离差可以写为

$\begin{gathered} f_i(x, y, z)=\sqrt{\left(x-x_i\right)^2+\left(y-y_i\right)^2+\left(z-z_i\right)^2}- \\ \sqrt{\left(x-x_0\right)^2+\left(y-y_0\right)^2+\left(z-z_0\right)^2} . \end{gathered}$

(16)

令r_i=f_i(x, y, z)，并进行Taylor展开，Taylor级数可以展开为无穷阶，应根据需要确定。X_k－1表示进行一阶展开基准点的三维坐标，k-1表示展开的阶数为1，所得方程即为定位方程，表示如下：

$\boldsymbol{F}(\boldsymbol{X})=\boldsymbol{F}\left(\boldsymbol{X}_{k-1}\right)+\boldsymbol{F}^{\prime}\left(\boldsymbol{X}_{k-1}\right)\left(\boldsymbol{X}-\boldsymbol{X}_{k-1}\right).$

(17)

可得

$\boldsymbol{F}^{\prime}\left(\boldsymbol{X}_{k-1}\right)\left(\boldsymbol{X}-\boldsymbol{X}_{k-1}\right)=\boldsymbol{F}(\boldsymbol{X})-\boldsymbol{F}\left(\boldsymbol{X}_{k-1}\right).$

(18)

如果进行二阶展开，选择基准点为X_k－2，则有

$\begin{aligned} \boldsymbol{F}(\boldsymbol{X})=& \boldsymbol{F}\left(\boldsymbol{X}_{k-2}\right)+\boldsymbol{F}^{\prime}\left(\boldsymbol{X}_{k-2}\right)\left(\boldsymbol{X}-\boldsymbol{X}_{k-2}\right)+\\ & \frac{1}{2} \boldsymbol{F}^{\prime \prime}\left(\boldsymbol{X}_{k-2}\right)\left(\boldsymbol{X}-\boldsymbol{X}_{k-2}\right)^2. \end{aligned}$

(19)

以一阶展开为例，其中：

$\boldsymbol{F}^{\prime}\left(\boldsymbol{X}_{k-1}\right)=\left[\begin{array}{ccc} \frac{\partial f_1}{\partial x} & \frac{\partial f_1}{\partial y} & \frac{\partial f_1}{\partial z} \\ \frac{\partial f_2}{\partial x} & \frac{\partial f_2}{\partial y} & \frac{\partial f_2}{\partial z} \\ \vdots & \vdots & \vdots \\ \frac{\partial f_N}{\partial x} & \frac{\partial f_N}{\partial y} & \frac{\partial f_N}{\partial z} \end{array}\right]_{\substack{x=x_{k-1} \\ y=y_{k-1} \\ z=k_{k-1}}},$

(20)

$\boldsymbol{X}-\boldsymbol{X}_{k-1}=\left[\begin{array}{c} x-x_{k-1} \\ y-y_{k-1} \\ z-z_{k-1} \end{array}\right].$

(21)

之后利用非线性最小二乘法解算上述定位方程，得到干扰源坐标的估计值。

3 仿真实验验证

目前已组建了4个无线电干扰智能监测站，并完成了对每个监测站的硬件、驱动以及软件的性能测试。监测站的硬件部分包括玻璃钢天线(70~900 MHz)、宽带全向监测天线(430 MHz~5.5 GHz)、授时天线和监测服务器。硬件的性能测试包括天线的接收频段、采样时间以及响应时间；驱动的性能测试包括驱动的软硬通信以及控制硬件启动；软件的性能测试包括软件的时间和空间复杂度、稳定性等。

本文基于已知类型的干扰信号数据和算法模型，对C-CNN模型的性能以及布站方式对定位精度的影响进行了仿真验证。

3.1 神经网络性能验证实验

本文对C-CNN的性能进行了测试，一方面，得到神经网络的输入大小、准确率以及收敛所需训练次数的关系；另一方面，将C-CNN与其他传统用于图像分类的神经网络模型进行对比。

图 4展示了输入长度为64、128和256时模型的测试准确率表现。输入长度为128的神经网络收敛时平均准确率为84.4%，输入长度为256的神经网络收敛时平均准确率为84.5%，二者差别不大。当输入长度为64时，准确率最低，收敛时平均准确率为79.3%。尽管增加输入长度会降低收敛速度，但在输入长度低于一定尺寸时，神经网络无法有效提取信息特征，使得识别准确率无法得到有效保证。故虽然在实验中输入长度为128的神经网络相比于输入长度为256的网络收敛速度快，识别准确率也相当，但在实际应用中为了保证能充分提取信息特征，最终选取网络的输入长度为256。

表 2和图 5展示了本文提出的C-CNN与其他神经网络模型如实数卷积神经网络(CNN)、全卷积网络(fully convolutional network, FCN)以及长短时记忆网络(long short-term memory networks, LSMT)等的性能对比结果。其中，表 2给出了这几种网络单个样本的识别准确率和平均运行时间，可以看出，基于C-CNN算法的神经网络模型的识别准确率最高，并且略高于其他神经网络模型的运行效率，这是因为其神经元较少。由于缺乏在复杂领域中提取信息特征的能力，CNN的识别准确率相较于C-CNN降低了7%。LSMT常用于处理语音等一维数据，在二维图像数据的处理上表现一般。综上所述，C-CNN不仅具有更高的识别准确率，而且运行的时间更短。

图 5给出了不同模型预测信号类型的混淆矩阵，包含背景噪声及其他3类干扰信号。C-CNN能够在背景噪声及其他3类干扰的信号中更准确地区分出信号，具有更好的识别能力和鲁棒性。

3.2 布站方式对定位精度的仿真实验

图 6展示了分别部署3、4、5及6个监测站对空间定位点精度的影响，实验模型采用的都是正多边形的布站形状。横纵坐标表示在Cartesian坐标系下的空间坐标(X, Y)，尖点表示站点的部署位置，热值图代表定位精度，等高线为几何精度因子(geometric dilution of precision, GDOP)，其值越小说明定位精度越高^[17]。由图 6可以直观地看出：3站点布站的定位区域小于其他3种布站方式，且定位误差远高于其他3种布站方式；6站点布站的定位区域最大，精度最高。但现实中在考虑了经济因素及计算复杂度的基础上，决定采用4站点布站方式。此外，目标源越靠近几何中心，布站方式对定位精度的影响越小，因此在实际应用时，布站应尽量保证信号源在站点的包围圈内，以获得较好的定位精度。

图 7为在Cartesian坐标系下4个监测站不同布局形状的GDOP等高线仿真结果。可以看出：相同大小的区域，不同布站方式之间GDOP值的范围差距不大，但所形成的定位区域的形状差别较大，监测站的具体分布还要综合考虑干扰源的主要分布范围和现场的实际地理情况。

4 结论

本文针对FAST电磁环境保护和无线电频谱管理的需求，对无线电干扰智能监测和定位技术展开了研究。

在干扰信号的探测和制式识别方面，研究了多站点融合的协同频谱感知技术，采用相干融合的方法实现了在低信噪比条件下的信号探测；在信号识别方面，基于现有的各种制式数据，结合深度神经网络强大的学习能力，提出了以瀑布图作为输入的C-CNN模型，实现了信号的特征提取和分类；并通过实验得到输入长度对模型收敛速度以及精度的影响。此外，还将本文提出的C-CNN模型与其他传统网络模型进行性能上的对比，并证明在根据瀑布图识别信号类别的任务中，C-CNN模型在准确性和运行速度方面均优于其他模型。

在干扰源高精度定位技术的研究中，从定位精度和布站方式2方面展开。本文通过广义互相关法和重构源信号的迭代算法估计信号到达时间差，基于Taylor级数展开的定位方程，得到了高精度的信号源位置估计；对不同布站方式的定位范围和精度等性能进行仿真实验，综合经济及计算复杂度方面的因素决定采用4站点的布站方式。

下一步，本团队将在FAST周边进行监测站的布站选址，并对多站之间进行联合监测实验，对干扰检测识别和定位精度进行实测，最终将其应用于FAST电磁环境的日常监测中，为FAST干扰信号数据库的建设积累数据，指导电磁环境保护和频谱管理工作。