中国电波传播研究所，山东 青岛

收稿日期：2021年5月17日；录用日期：2021年6月6日；发布日期：2021年6月17日

摘要

针对短波沃森瓦特测向系统，提出了一种基于波束形成的测向算法，该算法将沃森瓦特接收机的三个通道输出进行线性组合，形成阵列测量矢量，采用最小方差无失真波束形成器进行空域滤波测向。与传统沃森瓦特算法相比，该方法充分利用了爱德考克天线的阵列孔径。计算机仿真表明，所提算法的测向精度优于传统算法。

关键词

无线电测向，沃森瓦特，谱估计，短波

Direction Finding of Shortwave Emitters Based on Spectrum Estimation

Longtao Li

China Research Institute of Radiowave Propagation, Qingdao Shandong

Received: May 17^th, 2021; accepted: Jun. 6^th, 2021; published: Jun. 17^th, 2021

ABSTRACT

An algorithm based on beamforming is proposed for shortwave direction finding with a Watson-Watt (WW) system. It constructs an array measurement vector by the linear combination of the outputs of the three receiving channels of a WW system, and then utilizes the Minimum Variance Distortionless Response (MVDR) beamforming to perform spatial filtering. Compared to the traditional WW algorithm, the proposed method takes advantage of the full aperture of the Adcock antenna array, which leads to its better accuracy of direction finding, and it has been demonstrated by computer simulations.

Keywords:Radio Direction Finding, Watson-Watt, Spectrum Estimation, Shortwave

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

短波的频率范围为1.5 MHz到30 MHz，该频段的无线电波可以被电离层反射或者散射，可以用于长距离通信。短波可以用于大范围音频广播、天波超视距雷达、业余无线电通信、航空通信、海上通信等等 [1] [2] [3] [4]。

短波测向系统主要有两类形态，一种是以大型圆阵为基础的固定站，采用干涉仪和空间谱体制，占地面积较大，空间谱体制可对多个同频信号同时测向；另一种是以爱德考克(Adcock)天线阵列为基础的小型可搬移站或者车载/舰载站，主要采用沃森瓦特(Watson-Watt)体制，使用机动灵活，主要缺点是精度相比大型固定站有所下降 [5]。

沃森瓦特体制的发展迄今已有百年历史。1919年Adcock F提出利用四个天线组成的阵列的机械旋转来进行幅度测向 [6]。1925至1926年，Watson-Watt R改进了Adcock F的方法，通过观察示波器上的利萨如(Lissajous)图形的椭圆倾斜角来指示信号方向而不再需要机械旋转，并提出增加第五个天线来实现去模糊 [7]。随着现代数字信号处理技术的发展，现代沃森瓦特测向设备采用全数字的方式进行处理，其测向稳定度、精度、灵敏度都获得了显著提升，其主要部件有爱德考克天线、三通道测向处理接收机等 [8]。

本文着力于在不改变现代沃森瓦特测向设备软硬件技术体制的情况下，从信号处理的角度去提升其测向精度，将其角度估计问题转换成为波束形成问题，并采用多项式求根实现闭式求解。

2. 沃森瓦特测向设备原理

2.1. 爱德考克天线阵列

爱德考克天线阵列由五个单极子天线组成，其中一个天线(称为sensor天线)放置在中心，用于解180度模糊，剩下四个天线每对沿着参考坐标轴放置，不失一般性，假设Y轴沿着北–南(NS)方向，X轴沿着东–西(EW)方向(图1)。

图1. Adcock天线示意图

另外也可以采用一副交叉环天线来代替N-S天线和E-W天线，从而形成一种三阵元天线，被称为沃森瓦特天线，其中：交叉环天线为磁加载天线，两天线正交，采用高频低噪声、低失真半导体管组成的输入放大加驱动电路的形式，具有大信号输入时线性度高，非线性失真小的优点。沃森瓦特天线具有体积小、准确度和灵敏度高，性能稳定等特点，适用于短波车载无线电信号的监测和测向，因其后续接收处理与Adcock天线相似，本文不再赘述。

2.2. 三通道测向处理接收机

三通道测向处理接收机由三路接收通道模块、校准源模块、本振模块、数字处理板卡、单板计算机等模块组成，其中：三路接收通道分别用于接收N-S差信号、E-W差信号以及经过90度移相的sensor天线信号并将其变换到中频(一般为21.4 MHz)，校准源可以产生标准单音信号，用于三路接收通道幅相校准，本振模块用于提供变频所需要的本振信号，数字信号处理板卡用于模数转换、信号处理等，单板计算机用于控制、通信与数据存储等，如图2所示。

图2. 三通道测向处理接收机原理框图

具体原理如下：

1) 爱德考克天线接收外界环境中的短波信号，首先产生5路射频信号 $x_N\left(t\right)$ 、 $x_S\left(t\right)$ 、 $x_E\left(t\right)$ 、 $x_W\left(t\right)$ 和 $x_{sensor}\left(t\right)$，然后， $x_S\left(t\right)$ 经过180˚移相器后与 $x_N\left(t\right)$ 合路形成N-S差信号 $x_N\left(t\right)-x_S\left(t\right)$，$x_W\left(t\right)$ 经过180˚移相器后与 $x_E\left(t\right)$ 合路形成N-S差信号 $x_E\left(t\right)-x_W\left(t\right)$，$x_{sensor}\left(t\right)$ 经过90˚移相器后形成信号 $j{x}_{sensor}\left(t\right)$，上述三路射频信号的频率范围均为1.5~30 MHz；

2) 根据设定的接收机频率，本振模块产生22.9~51.4 MHz的本振信号，在短波信号1/2/3中分别与射频信号 $x_N\left(t\right)-x_S\left(t\right)$ 、 $x_E\left(t\right)-x_W\left(t\right)$ 、 $j{x}_{sensor}\left(t\right)$ 进行混频，混频输出为21.4 MHz的正常中频以及24.4~81.4 MHz的高中频，再经过中心频率为21.4 MHz、带宽为2 MHz的模拟滤波器后，得到中频输出信号并放大到合适的电平；

3) 三路中频信号进入数字信号处理板卡，进行3通道AD采集，采样率为5.12MHz，并经过4.12 MHz数字下变频(Digital Down Converter, DDC)，形成零中频IQ信号，然后经过抽取、级联积分梳状(Cascaded Integrator-Comb, CIC)滤波器、低通滤波器等，形成带宽为10KHz的数字信号；

4) 三路数字信号送入单板计算机，经过数字信号算法，形成测向结果，发送至上位机软件系统进行显示。

2.3. 基本测向算法原理

设短波信号相对于参考坐标系的入射方向为 ，那么在理想情况下，四个测向单极子天线和sensor天线的电压值分别为：

$x_N\left(t\right)=s\left(t\right){\text{e}}^{jkR\sin \theta }$ (1)

$x_S\left(t\right)=s\left(t\right){\text{e}}^{-jkR\sin \theta }$ (2)

$x_E\left(t\right)=s\left(t\right){\text{e}}^{jkR\cos \theta }$ (3)

$x_W\left(t\right)=s\left(t\right){\text{e}}^{-jkR\cos \theta }$ (4)

$x_{sensor}\left(t\right)=s\left(t\right)$ (5)

其中：R为天线到阵列中心的距离， $s\left(t\right)$ 为信号时域波形，k为波数(wavenumber)， $k=2\pi /\lambda$，$\lambda$ 为信号波长，j为虚部单位。

首先计算N-S差信号：

$x_N\left(t\right)-x_S\left(t\right)=s\left(t\right)\left({\text{e}}^{jkR\sin \theta }-{\text{e}}^{-jkR\sin \theta }\right)$ (6)

利用欧拉公式可得，

$x_N\left(t\right)-x_S\left(t\right)=2js\left(t\right)\sin \left(kR\sin \theta \right)$ (7)

在短波频段，波长很长，波数k很小。因为对于一个很小的角度 $\beta$，有 $\sin \beta \approx \beta$，所以：

$x_N\left(t\right)-x_S\left(t\right)\approx 2js\left(t\right)kR\sin \theta$ (8)

类似的，可以得到E-W差信号：

$x_E\left(t\right)-x_W\left(t\right)\approx 2js\left(t\right)kR\cos \theta$ (9)

由两路差信号可以得出：

$\theta ={\tan }^{-1}\left(\tan \theta \right)={\tan }^{-1}\left(\frac{2js\left(t\right)kR\sin \theta }{2js\left(t\right)kR\cos \theta }\right)\approx {\tan }^{-1}\left(\frac{x_N\left(t\right)-x_S\left(t\right)}{x_E\left(t\right)-x_W\left(t\right)}\right)$ (10)

考虑到噪声的存在以及信号本身存在的非平稳性，实际进行信号处理时，采用两路都进行FFT计算的方式，读取关心频点的谱线的幅度，计算两路幅度之间的比值，然后进行反正切运算。

由于反正切函数存在180度模糊，所以需要利用第三路(中心天线)进行判断。如果N-S差信号 $x_N\left(t\right)-x_S\left(t\right)$ 与经过90度移相的sensor天线信号 $j{x}_{Sensor}\left(t\right)=js\left(t\right)$ 同相，则方向位于NW/NE象限，反之则位于SW/SE象限。由于 $kR\sin \theta$ 很小，这种相位的超前/滞后的绝对值都不大，一般不会出现模糊。

3. 基于谱估计的短波无线电测向算法

首先做出如下假设，采集的IQ信号已经过通带足够窄的频域数字滤波，使得输入到测向算法模块的IQ信号中只包含有单个频率的短波信号。我们从阵列信号处理的角度重新对三路数字信号进行建模，其阵列输出矢量 $x\left(t\right)$ 可以表示为：

$x\left(t\right)=\left[\begin{array}{c}-\frac{x_N\left(t\right)-x_S\left(t\right)+j\left[x_E\left(t\right)-x_W\left(t\right)\right]}{2kR}\\ \frac{x_N\left(t\right)-x_S\left(t\right)-j\left[x_E\left(t\right)-x_W\left(t\right)\right]}{2kR}\\ -j\cdot j{x}_{Sensor}\left(t\right)\end{array}\right]\approx \underset{=a\left(\theta \right)}{\left[\begin{array}{c}{\text{e}}^{-j\theta }\\ {\text{e}}^{j\theta }\\ 1\end{array}\right]}s\left(t\right)+n\left(t\right)=a\left(\theta \right)s\left(t\right)+n\left(t\right)$ (10)

其中： $a\left(\theta \right)$ 为方向矢量， $n\left(t\right)$ 为阵列噪声矢量。

由式(10)可以看出，该模型有效利用了爱德考克天线的三维阵列孔径，而不是像传统沃森瓦特体制那样，仅仅将sensor天线用来解模糊。

对式(10)中表示的阵列输出矢量模型通常可以采用子空间类(Subspace)方法进行测向处理，常用的方法主要为多重信号分类(Multiple Signal Classification, MUSIC)方法 [9]，其基本原理是：阵列输出矢量的协方差矩阵的特征矢量(Eigenvectors)是完备的线性空间，其中，信号方向矢量的张成空间(Spanned Space)，与噪声对应的特征矢量的张成空间是正交的，也就是信号子空间与噪声子空间是正交子空间。考虑到 $x\left(t\right)$ 的三个分量的噪声水平并不相同，属于色噪声背景，采用子空间类方法在信噪比较低时容易出现子空间的误判，因此本文将利用波束形成的方法进行测向。设波束形成器的加权矢量为 $w$，则波束形成后滤出的信号为：

$y\left(t\right)=w^{\text{H}}x\left(t\right)$ (11)

其中： ${(\cdot )}^{\text{H}}$ 表示共轭转置。可以采用标准波束形成器进行波束扫描测向，设计准则为最小响应无失真(MVDR) [10]，即：

$\min \left\{{\left|y\left(t\right)\right|}^2\right\}=w^{\text{H}}Cw\text{s}\text{.t}.w^{\text{H}}a\left(\theta \right)=1$ (12)

其中： $C$ 为阵列输出协方差矩阵，

$C=E\left\{x\left(t\right)x^{\text{H}}\left(t\right)\right\}$ (13)

利用拉格朗日方法，易知

$w=\frac{C^{-1}a\left(\theta \right)}{a^{\text{H}}\left(\theta \right)C^{-1}a\left(\theta \right)}$ (14)

所以其扫描谱为：

$P\left(\theta \right)=w^{\text{H}}Cw=\frac{a^{\text{H}}\left(\theta \right)a\left(\theta \right)}{a^{\text{H}}\left(\theta \right)C^{-1}a\left(\theta \right)}$ (15)

通过寻找 $P\left(\theta \right)$ 的谱峰即可估计到达方向。由于

$a^{\text{H}}\left(\theta \right)a\left(\theta \right)=3$ (16)

所以，上述问题等价于寻找下述 $Q\left(\theta \right)$ 的最小值：

$Q\left(\theta \right)=a^{\text{H}}\left(\theta \right)C^{-1}a\left(\theta \right)$ (17)

令

$C^{-1}=\left[\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{21}& c_{22}& c_{23}\\ c_{31}& c_{32}& c_{33}\end{array}\right]$ (18)

于是有

$Q\left(\theta \right)=c_{12}{\text{e}}^{j2\theta }+c_{21}{\text{e}}^{-j2\theta }+\left(c_{13}+c_{32}\right){\text{e}}^{j\theta }+\left(c_{23}+c_{31}\right){\text{e}}^{-j\theta }+c_{11}+c_{22}+c_{33}$ (19)

令

$z={\text{e}}^{j\theta }$ (20)

于是有

$Q\left(z\right)=c_{12}{z}^2+c_{21}{z}^{-2}+\left(c_{13}+c_{32}\right)z+\left(c_{23}+c_{31}\right)z^{-1}+c_{11}+c_{22}+c_{33}$ (21)

由于谱谷对应的导数为零，所以可以利用下式求解：

$\begin{array}{c}\frac{\text{d}Q\left(z\right)}{\text{d}z}=\text{2}{c}_{12}z+\left(c_{13}+c_{32}\right)-\left(c_{23}+c_{31}\right)z^{-2}-\text{2}{c}_{21}{z}^{-3}\\ =z^{-3}\left[\text{2}{c}_{12}{z}^4+\left(c_{13}+c_{32}\right)z^3-\left(c_{23}+c_{31}\right)z-\text{2}{c}_{21}\right]\\ =0\end{array}$ (22)

由于 $z={\text{e}}^{j\theta }\ne 0$，所以，

$\text{2}{c}_{12}{z}^4+\left(c_{13}+c_{32}\right)z^3-\left(c_{23}+c_{31}\right)z-\text{2}{c}_{21}=0$ (23)

这是一个关于z的一元四次方程，存在四个解，对上式进行多项式求根，然后找到单位圆内(或外)最贴近单位元的根，求取幅角即可获得角度值估计。本文方法的原理框图如图3所示。

图3. 本文方法原理框图

4. 计算机仿真

假设爱德考克天线阵列半径为0.6 m；接收通道中频输出为21.4 MHz，采样率为5.12 MHz，采样点数为8192；入射信号的带宽10 kHz、调制类型为模拟调频。

图4、图5、图6给出了传统方法(即2.3节中介绍的沃森瓦特短波测向方法)和本文方法(即谱估计短波测向方法)的测向误差随着信噪比的变化曲线，其中信号中心频率分别为1.5 MHz、10 MHz以及30 MHz。由于噪声等因素的影响，测向结果为随机变量，通常认为测向误差 $\Delta \theta$ 近似服从零均值高斯分布 $N\left(0,{\sigma }^2\right)$，所以通常采用两种方式描述误差的大小：

1) 均方根误差(Root Mean Squared Error, RMSE)：也就是高斯分布 $N\left(0,{\sigma }^2\right)$ 的标准差 $\sigma$，即：测向误差以68%的概率小于等于RMSE误差；

2) p%误差：也就是测向误差以p%的概率小于等于该误差，通常p可以取85、90、95等，当p = 95时，误差相当于两倍RMSE。

本文选取RMSE误差作为测向误差性能指标进行仿真。仿真表明，本文方法与传统方法相比，精度取得明显提升，且收敛速度较快。

图4. 传统方法和本文方法的精度对比(1.5 MHz)

图5. 传统方法和本文方法的精度对比(10 MHz)

图6. 传统方法和本文方法的精度对比(30 MHz)

5. 结论

本文为沃森瓦特短波测向体制设计了一种数字信号处理算法，运算逻辑简单，便于算法在FPGA、DSP或者单板计算机等软硬件平台中实施。与传统算法相比，本文算法测向精度和随着信噪比的收敛速度都有显著提升。

文章引用

李龙涛. 基于谱估计的短波信号测向算法
Direction Finding of Shortwave Emitters Based on Spectrum Estimation[J]. 无线通信, 2021, 11(03): 43-50. https://doi.org/10.12677/HJWC.2021.113006

参考文献