利用Lissajous图形的雷达辐射源DOA与极化参数联合估计方法

Joint Estimation Method of Radar Emitter DOA and Polarization Parameters Based on Lissajous Figure

摘要

稀疏压缩感知类算法无法对相干信号进行波达方向（DOA）估计，极化参数精度估计受幅相误差影响严重。传统算法拘束研究在一维信号层面，对上述问题解决效果不佳，所以尝试利用Lissajous图形获取更高维度的信息。首先，分析信号在空间的电场矢量形成的Lissajous图形与DOA空间关系，结合离网压缩感知（OGJS）算法进行DOA参数估计，改善原算法因灵活度拓展带来的精度损失，信噪比为5 dB时DOA估计绝对平均误差由3.83°降至0.27°；然后，利用调参变分自编码器（β-VAE）对灰度图形式的Lissajous图形仿真数据集进行学习，实现对扰动数据的重构，消除幅相误差。对比卷积自编码器网络处理一维信号，实现信噪比为10 dB、极化相差分辨率在1.5°时，小样本100轮学习训练下，匹配率由80.4%提升至99.2%。

Abstract

Sparse compressed sensing algorithms cannot estimate the direction of arrival (DOA) of coherent signals, and the accuracy estimation of polarization parameter is greatly affected by amplitude and phase errors. Traditional algorithms, which are mostly restricted to the research of one-dimensional signals, fail to effectively solve the above problems. Lissajous figure is used to obtain higher dimensional information. Firstly, the spatial relationship between the Lissajous figure formed by the electric field vector of the signal in space and the DOA is analyzed, and the off-grid joint-sparse (OGJS) algorithm is used for DOA parameter estimation, which reduces the accuracy loss caused by the flexibility expansion of the original algorithm. When the signal-to-noise ratio (SNR) is 5 dB, the absolute average error of DOA estimation is reduced from 3.83° to 0.27°. Then, the β-variational auto-encoder (β-VAE) is employed to learn the simulation dataset of Lissajous figure in the form of grayscale images, so as to reconstruct the perturbed data and eliminate the amplitude and phase errors. Compared with the processing of one-dimensional signals by convolutional autoencoder networks, when the SNR is 10 dB and the polarization difference resolution is 1.5°, after 100 rounds of small sample learning training, the matching rate is improved from 80.4% to 99.2%.

Graphical abstract

关键词

Lissajous图形 / 离网压缩感知 / 变分自编码器 / 联合估计 / 雷达辐射源 / 波达方向

Key words

Lissajous figure / OGJS / VAE / joint estimation / radar emitter / DOA

现代战场电磁环境日益复杂，天线扫描方式和信号参数的变化给雷达信号的分选和识别带来很多困难^[1]。极化参数作为相对稳健的雷达信号特征，在一定程度上反映雷达的类别和工作模式，可结合信号波达方向（Direction of Arrival, DOA）参数进行同向、同频雷达信号的分选和识别^[2]。

基于极化敏感阵列的雷达辐射源DOA与极化参数联合估计方法近年来发展迅速：肖雄等^[3]提出一种基于三正交电偶极子天线多信号分类估计算法，从而引入了极化信息，但是存在计算复杂度高、耗时长的问题。逯岩斌等^[4]提出一种基于互质极化敏感阵列的模值约束降维求根MUSIC参数联合估计算法，利用互质阵列解DOA参数模糊，但实际机器计算中，由于约束矩阵接近奇异，极化参数估计误差较大。郭玉华等^[5]将极化域信息和空域信息相结合，利用目标信号和干扰信号的极化域和空域特征差别提高了参数估计精度，但无法处理相干干扰。苑清扬^[6]提出空间平滑算法恢复存在相干干扰下信号协方差矩阵的秩，从而提高参数估计精度，郑秩松等^[7]从空域滤波角度将DOA参数估计精度进一步提高，但两者都是通过阵列灵活度损失换取相干源的分离。Chalise等^[8]将压缩稀疏感知算法与极化敏感阵列相结合，拓展了阵列灵活度，但参数估计精度受网格限制。Wagner等^[9]通过将基于原子范数最小化算法应用到稀疏阵列中，相比压缩感知类算法达到了更高的性能。肖思棋等^[10]建立离网接收信号模型，使参数估计不再受到网格限制，提高参数估计精度，但由于有限等距性质（Restricted Isometry Property, RIP）准则^[11]的稀疏要求，精度依然不高。另外，针对极化敏感阵列由于阵元空间位置、阵元抖动造成的通道幅相误差近年来不断有专家学者利用深度学习对其进行研究消除。在前端处理中，最新成果有利用长短期记忆网络（Long Short-Term Memory, LSTM）捕捉极化敏感阵列的幅相误差随时间变化的动态特性^[12]和设计深度自编码器网络重构校准后的信号^[13]，其中深度自编码器网络通过设计损失函数驱动模型学习误差特征，无需依赖人工标注的误差参数，更便于训练和模型更新。

基于三正交电偶极子接收合成的空间信号Lissajous图形，直接完整地表现雷达信号的空间矢量信息和极化特征，可利用其估计雷达信号的DOA和极化参数。在雷达辐射源信号参数估计中，模型的稳定性和时效性是两个重要指标，相较传统自编码器网络结构，变分自编码器（Variational Auto-Encoder, VAE）具有收敛速度快、模型稳定性强等优势。本文尝试利用Lissajous图形对来波信号在空间的电磁波矢量状态进行表征提取，从新角度研究雷达辐射源信号极化状态。首先，通过分析三正交电偶极子形成的三维Lissajous图形确定来波信号在空间的传播方向，从理论和实验上论证利用Lissajous图形信息提取来波信号DOA信息的可行性，并将其用于改进离网压缩感知算法（Off-Grid Joint-Sparse, OGJS），相当于在算法中加入空域约束，突破RIP准则限制，提高参数估计精度。最后，利用调参变分自编码器（β-VAE）络对Lissajous图形进行降噪，通过匹配的方式估计信号极化参数。

1 信号Lissajous图形特征分析与建模

1.1 基于Lissajous图获取信号极化特征

数学上，Lissajous图形^[14]是两个沿着互相垂直方向的正弦振动的合成的轨迹。Lissajous图形由以下参数方程定义：

\left\{\begin{array}{l}x(\theta )=a\mathrm{s}\mathrm{i}\mathrm{n}(\theta );\\ y(\theta )=b\mathrm{s}\mathrm{i}\mathrm{n}(n\theta +\varphi ).\end{array}\right.（1）

式中：\theta表示角度变量；\varphi表示相位差，\mathrm{0}\le \varphi \le \frac{\mathrm{\pi }}{\mathrm{2}}；a、b表示幅度；n表示图形的参数，是两个正弦振动的频率比，n\ge \mathrm{1}。

若n为有理数，则参数方程可以写为

\left\{\begin{array}{l}x(\theta )=a\mathrm{s}\mathrm{i}\mathrm{n}(p\theta );\\ y(\theta )=b\mathrm{s}\mathrm{i}\mathrm{n}(q\theta +\varphi ).\end{array}\right.（2）

式中：\mathrm{0}\le \theta \le \mathrm{2}\mathrm{\pi }；\mathrm{0}\le \varphi \le \frac{\mathrm{\pi }}{\mathrm{2}p}。

对于式（1），若a=b、n=\mathrm{1}，则这两个正弦振动所合成的Lissajous图形是椭圆。由Lissajous图形可以求出这两个正弦振动的相位差^[15]，其原理如下：当n\theta =\mathrm{0}时，x\left(\theta \right)=a\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{0}=\mathrm{0}，y\left(\theta \right)=b\mathrm{s}\mathrm{i}\mathrm{n}\varphi。由此可求出x\left(\theta \right)在x轴上的截距A=b\mathrm{s}\mathrm{i}\mathrm{n}\varphi，A表示归一化的截距。设质点在水平方向上的最大偏移量为B，则有B=b。综上，有：

\frac{A}{B}=\frac{b\mathrm{s}\mathrm{i}\mathrm{n}\varphi }{b}=\mathrm{s}\mathrm{i}\mathrm{n}\varphi（3）

\varphi =\mathrm{s}\mathrm{i}{\mathrm{n}}^{-\mathrm{1}}\frac{A}{B}（4）

假设雷达辐射源发射的电磁波为单色完全极化电磁波，沿z轴正方向传播，其矢量形式为

\mathit{E}=E_x{\mathit{e}}_x+E_y{\mathit{e}}_y（5）

式中，{\mathit{e}}_x、{\mathit{e}}_y分别为x轴和y轴正方向的单位基矢量。E_x、E_y分别为：

E_x=E_{\mathrm{0}x}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\mathrm{\pi }{f}_{\mathrm{0}}t-Kz+{\phi }_x\right)（6）

E_y=E_{\mathrm{0}y}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\mathrm{\pi }{f}_{\mathrm{0}}t-Kz+{\phi }_y\right)（7）

式中：f_{\mathrm{0}}表示电磁波频率；K表示空间传播常数；z表示在传播方向上的相位延迟；E_{\mathrm{0}x}和{\phi }_x分别表示x轴正方向上的电场振幅和相位；E_{\mathrm{0}y}和{\phi }_y分别表示y轴正方向上的电场振幅和相位。

将式（5）写成复矢量形式

\mathit{e}\left(t,z\right)=\mathit{e}\cdot {\mathrm{e}}^{\mathrm{j}\left(\mathrm{2}\mathrm{\pi }{f}_{\mathrm{0}}t-Kz\right)}（8）

\mathit{e}=\left[\begin{array}{c}{\mathit{e}}_x\\ {\mathit{e}}_y\end{array}\right]=\left[\begin{array}{c}{e}_{\mathrm{0}x}{\mathrm{e}}^{\mathrm{j}{\phi }_x}\\ e_{\mathrm{0}y}{\mathrm{e}}^{\mathrm{j}{\phi }_y}\end{array}\right]（9）

式中：\mathit{e}表示复电场矢量；e_{\mathrm{0}x}和{\phi }_x分别表示x轴正方向上的电场振幅和相位；e_{\mathrm{0}y}和{\phi }_y分别表示y轴正方向上的电场振幅和相位。x轴和y轴正方向的单位基矢量的比值为

\eta =\frac{{\mathit{e}}_x}{{\mathit{e}}_y}=\frac{e_{\mathrm{0}y}{\mathrm{e}}^{\mathrm{j}{\phi }_y}}{e_{\mathrm{0}x}{\mathrm{e}}^{\mathrm{j}{\phi }_x}}=\mathrm{t}\mathrm{a}\mathrm{n}\gamma \cdot {\mathrm{e}}^{\mathrm{j}\phi }（10）

式中：\mathrm{t}\mathrm{a}\mathrm{n}\gamma ={\mathit{e}}_y/{\mathit{e}}_x，表示极化比，\gamma \in [\mathrm{0,2}\mathrm{\pi })；\phi ={\phi }_y-{\phi }_x，表示极化相差，\phi \in \left[\mathrm{0,2}\mathrm{\pi }\right]。

\left(\gamma ,\phi \right)即为电磁波的琼斯矢量。雷达辐射源来波信号电场矢量端点随时间随机变化所形成的轨迹为一参数恒定的椭圆。调整接收方向，可以使得极化比为1，Lissajous图形就是在这种情况下接收天线处的电场矢量端点轨迹图，信号极化相差即为电场在x、y方向正交分解后两电场矢量的相位差。

1.2 三正交电偶极子天线雷达辐射源信号接收模型

电磁矢量传感器（Electro Magnetie Vector Senor, EMVS）包含3个电偶极子和3个磁偶极子，可以同时接收6维电磁场矢量。假设EMVS构成的阵元位于坐标中心，3个电偶极子方向分别与x、y、z三坐标轴平行，利用琼斯矢量描述电磁波的极化信息，极化信息描述子为\left(\gamma ,\eta \right)，假定信源数为k，来波方位、俯仰角分别为{\phi }_k、{\theta }_k，沿单位方向矢量{\mathit{a}}_k的反方向入射：

{\mathit{a}}_k=\left(\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_k\right)\mathit{x}+\left(\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_k\right)\mathit{y}+\left(\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\right)\mathit{z}（11）

传播方向上两正交单位电场矢量分别为：

{\mathit{a}}_{{\theta }_k}=\left(\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_k\right)\mathit{x}+\left(\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_k\right)\mathit{y}-\left(\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\right)\mathit{z}（12）

{\mathit{a}}_{{\phi }_k}=\left(-\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_k\right)\mathit{x}+\left(\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_k\right)\mathit{y}（13）

式中，\mathit{x}\mathrm{、}\mathit{y}\mathrm{、}\mathit{z}表示三坐标轴的单位矢量。6维电场矢量\mathit{e}\left(k\right)和磁场矢量\mathit{h}\left(k\right)构成的极化—角度导向矢量{\mathit{\alpha }}_k可表示为

{\mathit{\alpha }}_k=\left[\begin{array}{c}\mathit{e}\left(k\right)\\ \mathit{h}\left(k\right)\end{array}\right]=\left[\begin{array}{c}{E}_x\\ E_y\\ E_z\\ H_x\\ H_y\\ H_z\end{array}\right]=\left[\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_k& -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_k\\ \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_k& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_k\\ -\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_k& -\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_k\\ \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_k& -\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_k\\ \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\end{array}\right]\left[\begin{array}{c}{E}_{{\theta }_k}\\ E_{{\phi }_k}\end{array}\right]（14）

式中：E_x、E_y、E_z表示平行于x、y、z轴的电场分量；H_x、H_y、H_z表示平行于x、y、z轴的磁场分量。

琼斯矢量可表示为

\mathit{J}={\left[E_{{\theta }_k},E_{{\phi }_k}\right]}^{\mathrm{T}},E_{{\theta }_k}={\mathrm{e}}^{\mathrm{j}\eta }\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_k,E_{{\phi }_k}=\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_k（15）

由于阵元所对应的通道并不完全平行，会出现各通道的不一致性，从而产生通道幅相误差。该误差与通道增益不一致造成的误差形式相同，因此可以在极化—角度导向矢量{\mathit{\alpha }}_k前乘以与方位无关的时变幅相误差矢量{\mathit{C}}_t^{\text{'}}，来表示含有时变幅相误差的极化—角度导向矢量{\tilde{\mathit{\alpha }}}_k。{\mathit{C}}_t^{\text{'}}和{\tilde{\mathit{\alpha }}}_k的表达式为

{\mathit{C}}_t^{\text{'}}={\left[\mathrm{1},{\mathit{C}}_{\mathrm{2}},\cdots ,{\mathit{C}}_t\right]}^{\mathrm{T}}=\left[\mathrm{1},{\rho }_{\mathrm{2}}{\mathrm{e}}^{\mathrm{j}{\phi }_{\mathrm{2}}},\cdots ,{\rho }_t{\mathrm{e}}^{\mathrm{j}{\phi }_t}\right]（16）

{\tilde{\mathit{\alpha }}}_k={\mathit{\alpha }}_k\cdot {\mathit{C}}_t^{\text{'}}（17）

式中：{\rho }_t表示幅度误差；{\mathrm{e}}^{\mathrm{j}{\phi }_t}表示相位误差。两者均随时间变化。

在实际情况中，通道间不平行常常是天线载体以一定规律振动造成的，故还可参考天线振动模型^[16]进行设置：

C\left(x,t\right)=Z_{\mathrm{0}}\left(x\right)+\alpha \left(t\right)+q_{\mathrm{1}}\cdot Z_{\mathrm{1}}\left(t\right)（18）

式中：Z_{\mathrm{0}}\left(x\right)表示低频的弯曲模式振动；\alpha \left(t\right)={\omega }_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t\right)表示位于坐标原点处的强制振动；q_{\mathrm{1}}表示振动模式系数；Z_{\mathrm{1}}\left(t\right)表示高频第一振动模式。由于通道时变幅相误差对接收信号的本质影响是乘性噪声，为简化模型，只考虑强制振动部分：

C\left(t\right)=\alpha \left(t\right)={\omega }_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t\right)（19）

因此，含有时变幅相误差的信号数据表示为

\tilde{\mathit{z}}\left(t\right)={\tilde{\mathit{\alpha }}}_k\cdot C\left(t\right)\cdot s_k\left(t\right)+n\left(t\right)（20）

式中：C\left(t\right)表示通道时变幅相误差；s_k\left(t\right)表示原始信号；n\left(t\right)表示加性噪声。

1.3 利用Lissajous图形的雷达辐射源DOA与极化参数联合估计流程

利用Lissajous图形的雷达辐射源DOA与极化参数联合估计主要分3步。首先，依据改进的OGJS算法进行DOA估计；其次，对信号进行空域滤波；最后，对空域滤波后提取出的单一信号的空间Lissajous图形在来波方向投影，运用VAE算法进行通道时变幅相误差校正，匹配对应参数的标准Lissajous图形得到极化相差。具体流程如下。

步骤1：利用信号空间Lissajous点云图得到法向单位矢量；

步骤2：利用法向单位矢量构造修正项，改进OGJS算法估计来波信号DOA；

步骤3：利用空域滤波得到该DOA的信号空间Lissajous点云图；

步骤4：利用β-VAE对该DOA的信号空间Lissajous点云图进行幅相误差消除和重构，并匹配估计其极化相差。

2 通道时变幅相误差校正

2.1 误差校正网络结构

经实验对比取优，通道时变幅相误差校正最终设计了7层网络结构：输入层、两层全连接层、数据特征提取层、两层全连接层和输出层，如图1所示。

数据处理模块将投影在来波方向上的Lissajous点云图作网格划分，再对每个网格内的散点作数量统计，再将散点数量映射为相应灰度：

g=\mathrm{225}\times \frac{n_i}{n_{\mathrm{m}\mathrm{a}\mathrm{x}}}（21）

式中：g表示灰度值；n_i表示第i个格子中散点数；n_{\mathrm{m}\mathrm{a}\mathrm{x}}表示所有计数网格中最多散点网格的散点数。

将所有计数网格中的散点数量作归一化处理，以及将散点图化为灰度图。考虑到数据为二维图形数据，噪声点集中在真实值周围，经过实验验证，将网络参数学习模块的编码过程设计为两层全连接层结构，以便更好地消除噪声影响。将M\times N^{\mathrm{2}}维数据转换为M\times N^{\mathrm{2}}\times \mathrm{1}维，利用全连接层进行特征提取工作。通过一维特征提取，以same模式对N^{\mathrm{2}}\times \mathrm{1}维的原始数据进行处理，保证充分提取特征的前提下，前后的样本尺度不变。全连接层通道数分别为4 000、80，分别提取出M\times \mathrm{4}\mathrm{000}维、M\times \mathrm{80}维特征，之后采用ReLU函数进行分别激活，由于层数较多，在全连接层与激活函数之间设置批正则化（Batch Normalization）层，降低神经网络训练过程中的梯度弥散问题，让训练结果更加稳定可靠。

ReLU激活函数采用具有分段线性的整流线性单元^[17]，是比较通用的激活函数，具有较强的鲁棒性和泛化能力；tanh激活函数具有更强的稳定性和饱和性，利于数据重构。在于前两层的激活函数为ReLU，最后一层的激活函数为tanh，这样设计兼顾了收敛速度和收敛效果，也是比较流行的设计方法。解码时经过3层全连接层输出的重构数据尺度为M\times N^{\mathrm{2}}\times \mathrm{1}维，再经过数据尺度变换得到M\times N^{\mathrm{2}}维重构数据。

2.2 β-VAE算法优势及理论推导

作为深度生成模型，变分自编码器（Variational Auto-Encoder, VAE）^[18]是基于变分贝叶斯推断的生成式网络结构。与传统的自编码器网络结构通过数值向量的方式表征潜在空间不同，这种通过表征潜在空间重构数据的方式具有较好的稳定性和收敛效果。β-VAE算法是对VAE算法损失函数进行优化得到的。如果将整个系统产生某组数据的概率定义为P\left(x\right)，用P\left(z\right)表示从标准正态分布中随机采样得到编码z的概率。由最大似然理论，系统目标函数如下：

\begin{array}{l}\mathrm{m}\mathrm{a}\mathrm{x}L=\sum \mathrm{l}\mathrm{g}P\left(x\right)\\ \mathrm{s}.\mathrm{t}.P\left(x\right)={\int }_z^{x}P\left(z\right)P\left(x\left|z\right.\right)\mathrm{d}z(\mathrm{22})\end{array}

用q表示编码器中先验、后验概率，用P表示解码器中先验、后验概率，则有：

\begin{array}{l}\mathrm{l}\mathrm{g}P\left(x\right)={\int }_{z}q\left(z\left|x\right.\right)lg\left(\frac{P\left(z,x\right)}{q\left(z\left|x\right.\right)}\right)\mathrm{d}z+\\ {\int }_{z}q\left(z\left|x\right.\right)lg\left(\frac{q\left(z\left|x\right.\right)}{P\left(z\left|x\right.\right)}\right)\mathrm{d}z=\\ L_{\mathrm{b}}+\mathrm{K}\mathrm{L}\left(q\left(z\left|x\right.\right)‖P\left(z\left|x\right.\right)\right)(\mathrm{23})\end{array}

\begin{array}{l}{L}_{\mathrm{b}}=\int q\left(z\left|x\right.\right)\mathrm{l}\mathrm{g}\left(\frac{P\left(x\left|z\right.\right)P\left(z\right)}{q\left(z\left|x\right.\right)}\right)\mathrm{d}z=\\ -\mathrm{K}\mathrm{L}\left(q\left(z\left|x\right.\right)\Vert P\left(z\right)\right)+\int q\left(z\left|x\right.\right)\mathrm{l}\mathrm{g}P\left(x\left|z\right.\right)\mathrm{d}z(\mathrm{24})\end{array}

对于式（24）右第1项，调整q\left(z\left|x\right.\right)使L_{\mathrm{b}}增大，由KL散度的定义，也就是使q\left(z\left|x\right.\right)接近于P\left(z\right)，将P\left(z\right)设置为标准正态分布，损失函数为

l_{\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{1}}}={\displaystyle \sum _{i=\mathrm{1}}^l}\left(\mathrm{e}\mathrm{x}\mathrm{p}\left({\sigma }_i\right)-\left(\mathrm{1}+{\sigma }_i\right)+m_i^{\mathrm{2}}\right)（25）

式中，{\sigma }_i和m_i可以看作是现实分布和正态分布间的方差、均值偏差，显然二者均为零时损失函数达到最小。对于式（24）右第2项，调整P\left(x|z\right)使L_{\mathrm{b}}增大，由蒙特卡罗方法可以得到损失函数为

\begin{array}{l}{l}_{\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{2}}}=\int q\left(z\left|x\right.\right)lgP\left(x\left|z\right.\right)dz=\\ E_{q\left(z\left|x\right.\right)}\left[\mathrm{l}\mathrm{g}P\left(x\left|z\right.\right)\right]\simeq \frac{\mathrm{1}}{L}{\displaystyle \sum _{l=\mathrm{1}}^L}\mathrm{l}\mathrm{g}P\left(x^{\left(i\right)}\left|z^{\left(i,l\right)}\right.\right)(\mathrm{26})\end{array}

式中，L为数据长度。式（26）在实际训练中用重构数据与原数据之间的交叉熵代替，也就是传统自编码器的损失函数。显然，传统自编码器没有建立表征潜在空间的损失函数来约束迭代过程（即式（25）），这也就从机理上解释了变分自编码器收敛速度快、稳定性强的原因。

综上，β-VAE算法的损失函数可表示为

\begin{array}{l}{l}_{\beta -\mathrm{V}\mathrm{A}\mathrm{E}}={\displaystyle \sum _{i=\mathrm{1}}^l}\left(\mathrm{e}\mathrm{x}\mathrm{p}\left({\sigma }_i\right)-\left(\mathrm{1}+{\sigma }_i\right)+m_i^{\mathrm{2}}\right)-\\ \beta {\displaystyle \sum _{i=\mathrm{1}}^N}\left[x_i\mathrm{l}\mathrm{g}\left({\stackrel{̑}{x}}_i\right)+(\mathrm{1}-x_i)\mathrm{l}\mathrm{g}\left(\mathrm{1}-{\stackrel{̑}{x}}_i\right)\right](\mathrm{27})\end{array}

式中，\beta表示超参数。相较于传统VAE算法损失函数中网络和数据两部分似然函数比重固定，β-VAE算法通过调整\beta的值平衡特征提取和图形重构效果两者权重，\beta越大说明该网络更加注重图像重构效果，但网络的泛化能力和稳定性会随之下降。实验获取最佳\beta值，使β-VAE算法在Lissajous图形重构效果达到最优。

3 OGJS算法的优化

压缩感知算法应用于稀疏阵列中，解决传统子空间类算法损失自由度的问题。其设计一个感知矩阵来压缩表示域中稀疏的信号再还原，在对压缩后信号采样时可以低于奈奎斯特采样频率，这是压缩感知算法的两大优势。OGJS算法对网格压缩感知理论中的网格失配问题进行改进。本文利用空间Lissajous图为原本凸优化算法添加空域约束项，突破RIP准则的限制，实现了更为精确DOA估计。

3.1 OGJS的DOA估计算法

将空域网格化后，由稀疏阵列的DOA估计理论，先将接收信号的协方差矩阵写成向量化形式再进行分解：

\begin{array}{l}{\mathit{r}}_{xx}=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{R}}_{xx}\right)=\left({\mathit{\varphi }}^{\mathrm{*}}\otimes \mathit{\varphi }\right)\left({\left(\overline{\mathit{A}}\right)}^{\mathrm{*}}\circ \overline{\mathit{A}}\right)\mathit{s}+\\ p_n\left({\mathit{\varphi }}^{\mathrm{*}}\otimes \mathit{\varphi }\right)\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{I}}_M\right) \mathit{\Psi }\tilde{\mathit{A}}\mathit{s}+p_n\mathit{\Psi }\mathit{l}(\mathrm{28})\end{array}

式中：\mathit{\varphi }表示压缩矩阵；\mathit{\Psi }={\mathit{\varphi }}^{\mathrm{*}}\otimes \mathit{\varphi }；\tilde{\mathit{A}}={\overline{\mathit{A}}}^{\mathrm{*}}\circ \overline{\mathit{A}}；\mathit{s}={\left[p_{\mathrm{1}},p_{\mathrm{2}},\cdots ,p_k\right]}^{\mathrm{T}}表示由信源功率组成的向量；\mathit{l}=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{I}}_M\right)表示噪声向量。将空域网格化，即\overline{\mathit{\theta }}=\left[{\overline{\theta }}_{\mathrm{1}},{\overline{\theta }}_{\mathrm{2}},\cdots ,{\overline{\theta }}_N\right]，其中N\gg K。设网格间距为\mathrm{2}r，即{\overline{\theta }}_{i+\mathrm{1}}-{\overline{\theta }}_i=\mathrm{2}r，\mathrm{1}\le i\le N-\mathrm{1}，如果信号真实入射方向为\mathit{\theta }=\left[{\theta }_{\mathrm{1}},{\theta }_{\mathrm{2}},\cdots ,{\theta }_k\right]，真实阵列流型矩阵记作\mathit{A}，\overline{\mathit{A}}为网格化流形矩阵，偏差矢量记作\mathit{D}=\mathit{\theta }-\overline{\mathit{\theta }}。

利用一阶泰勒展开原理：

\mathit{A}\left(\mathit{\theta }\right)=\mathit{A}\left(\overline{\mathit{\theta }}\right)+\mathit{B}\left(\overline{\mathit{\theta }}\right)\mathit{D}（29）

即

{\mathit{r}}_{xx}=\mathit{\Psi }\left(\tilde{\mathit{A}}+\tilde{\mathit{B}}\mathit{D}\right)\mathit{s}+p_n\mathit{\Psi }\mathit{l}（30）

式中，

\tilde{\mathit{B}}=\left\{\frac{\partial \left[{\mathit{a}}^{\mathrm{*}}\left({\overline{\theta }}_{\mathrm{1}}\right)\otimes \mathit{a}\left({\overline{\theta }}_{\mathrm{1}}\right)\right]}{\partial {\overline{\theta }}_{\mathrm{1}}},\cdots ,\frac{\partial \left[{\mathit{a}}^{\mathrm{*}}\left({\overline{\theta }}_N\right)\otimes \mathit{a}\left({\overline{\theta }}_N\right)\right]}{\partial {\overline{\theta }}_N}\right\}（31）

至此，根据压缩感知信号重构理论^[19]可以得到信号源\mathit{s}的优化模型：

\begin{array}{l}\underset{\mathit{s},\mathit{D}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{\Psi }\left(\tilde{\mathit{A}}+\tilde{\mathit{B}}\mathit{D}\right)\mathit{s}+p_n\mathit{\Psi }\mathit{l}-{\mathit{r}}_{xx}‖}_{\mathrm{2}}^{\mathrm{2}}+\mu {‖\mathit{s}‖}_{\mathrm{1}}\\ \mathrm{s}.\mathrm{t}.-r\le \mathit{s}\le r,-r\le \mathit{D}\le r(\mathrm{32})\end{array}

上式存在双线性变量，属于非凸问题。OGJS为解决这一问题定义了一个修正向量\mathit{p}=\mathit{D}\odot \mathit{s}，令\mathit{x}={\left[{\mathit{s}}^{\mathrm{T}},{\mathit{p}}^{\mathrm{T}}\right]}^{\mathrm{T}}，其中\mathit{p}与\mathit{s}具有相同的稀疏结构，并且两者中非零值位置相同。定义{‖\mathit{x}‖}_{\mathrm{2,1}}为\mathit{x}的l_{\mathrm{2}}/l_{\mathrm{1}}混合范数^[20]:

{‖\mathit{x}‖}_{\mathrm{2,1}}={\displaystyle \sum _{i=\mathrm{1}}^N}\sqrt[]{x_i^{\mathrm{2}}+x_{N+\mathrm{1}}^{\mathrm{2}}}（33）

将式（32）转换为一个凸优化问题：

\begin{array}{l}\underset{\mathit{s},\mathit{D}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{\Psi }\left(\tilde{\mathit{A}}\mathit{s}+\tilde{\mathit{B}}\mathit{p}\right)+p_n\mathit{\Psi }\mathit{l}-{\mathit{r}}_{xx}‖}_{\mathrm{2}}^{\mathrm{2}}+\mu {‖\mathit{x}‖}_{\mathrm{2,1}}\\ \mathrm{s}.\mathrm{t}.-r\left|\mathit{s}\right|\le \mathit{p}\le r\left|\mathit{s}\right|,\mathit{s}\ge \mathrm{0},\mathit{x}={\left[{\mathit{s}}^{\mathrm{T}},{\mathit{p}}^{\mathrm{T}}\right]}^{\mathrm{T}}(\mathrm{34})\end{array}

OGJS算法如算法1所示。

算法1 OGJS算法

输入：压缩测量向量\mathit{x}\left[t\right]

输出：DOA估计值

1.向量化{\mathit{R}}_{xx}得到{\mathit{r}}_{xx}

2.定义修正矢量\mathit{p}

3.求解凸优化问题最优解\widehat{\mathit{s}}、\widehat{\mathit{p}}

4.由\widehat{\mathit{s}}、\widehat{\mathit{p}}得到\mathit{D}，搜寻\widehat{\mathit{s}}中前k个最大值，利用\mathit{\theta }=\overline{\mathit{\theta }}+\mathit{D}得到离网DOA估计

3.2 基于空间Lissajous图形优化OGJS的DOA估计算法

仿真来波方向为\theta =\mathrm{45}°、\phi =\mathrm{45}°的两个信号，极化信息描述子分别为\left(\gamma =\mathrm{45}°,\eta =\mathrm{0}°\right)（线极化）、\left(\eta =\mathrm{90}°\right)（圆极化），接收时使用20倍采样，考虑通道时变幅相误差和噪声。

在纯净的Lissajous点云图中，线极化波的Lissajous图是一维直线，无法通过找点云所在平面法线的方式进行来波方向估计。但由于实际中存在噪声，实际接收中线极化波的Lissajous点云图也可以确定平面进行DOA估计。存在噪声的情况下，Lissajous点云图总是严格的在一个平面上，下面进行理论阐述。实际接收中，三正交电偶极子三通道接收到的是电场矢量信息，由式（14）可知，对于第k个来波信号有：

{\mathit{e}}_k=\left[\begin{array}{c}{E}_x\\ E_y\\ E_z\end{array}\right]=\left[\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_k& -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_k\\ \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_k& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_k\\ -\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k& \mathrm{0}\end{array}\right]\left[\begin{array}{c}{E}_{{\theta }_k}\\ E_{{\phi }_k}\end{array}\right]（35）

Lissajous点云矢量即所接收电场矢量的集合为\left\{{\mathit{z}}_k\left(t\right)\left|{\mathit{z}}_k\left(t\right)\right.\in p_n{\mathit{e}}_k\cdot \left[s_k\left(t\right)+n\left(t\right)\right]\right\}。

单位方向矢量{\mathit{a}}_k未知，p_n为信源功率，易得{\mathit{a}}_k\cdot {\mathit{e}}_k=\mathrm{0}，这解释了存在噪声的情况下Lissajous点云图总是严格在一个平面上的原因。也就是说，求DOA问题可以转化为一个凸优化问题：

\begin{array}{l}\underset{\mathit{\theta },\mathit{\phi }}{\mathrm{m}\mathrm{i}\mathrm{n}}\sum \left|{\mathit{a}}_k\cdot \left\{{\mathit{z}}_k\left(t\right)\left|{\mathit{z}}_k\left(t\right)\right.\in {\mathit{e}}_k\cdot \left[s_k\left(t\right)+n\left(t\right)\right]\right\}\right|\\ \mathrm{s}.\mathrm{t}.\hspace{1em}\left|{\mathit{a}}_k\right|=\mathrm{1}(\mathrm{36})\end{array}

上述讨论中，单一信号理论上最少只需要3个采样点确定点云平面即可确定DOA；多信号时，结合压缩感知理论，向量化协方差矩阵，如式（37）所示，这里省略压缩矩阵。式（38）是电磁场的能量传递的能量守恒式。

{\mathit{r}}_{xx,yy,zz}=\left({\left(\overline{\mathit{A}}\right)}^{\mathrm{*}}\circ \overline{\mathit{A}}\right){\mathit{s}}_{x,y,z}+p_n\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{I}}_M\right)（37）

\overline{\mathit{A}}\cdot {\mathit{s}}_{x,y,z}={\overline{\mathit{A}}}_{\stackrel{̑}{\mathit{a}}}\frac{\mathrm{1}}{T}\sum {\left|{\mathit{z}}_{k_{x,y,z}}\left(t\right)\right|}^{\mathrm{2}}（38）

式中：{\mathit{s}}_{x,y,z}表示三通道的稀疏信号功率向量；\stackrel{̑}{\mathit{a}}表示磁场方向，可以看作反映信号的DOA，即有

\sum \left|\stackrel{̑}{\mathit{a}}\cdot \left\{{\stackrel{̑}{z}}_k\left(t\right)\left|{\stackrel{̑}{z}}_k\left(t\right)\right.\in \sum {\mathit{e}}_k\cdot \left[s_k\left(t\right)+n\left(t\right)\right]\right\}\right|=\mathrm{0}（40）

综上可以在传统压缩感知形成的凸优化目标函数上加上修正项，可表示为

\underset{\mathit{s}}{\mathrm{m}\mathrm{i}\mathrm{n}}\lambda \left|\overline{\mathit{A}}\cdot \mathit{s}-{\overline{\mathit{A}}}_{\stackrel{̑}{\mathit{a}}}\frac{\mathrm{1}}{T}\sum {\left|s_k\left(t\right)+n\left(t\right)\right|}^{\mathrm{2}}\right|（41）

4 实验分析

4.1 通道时变误差校正效果分析

仿真极化相差分别为0°（线极化）、1.5°（椭圆极化），极化比均为1的两目标辐射源信号数据集和测试集各100份。将数据集在信噪比（Signal-to-Noise Ratio, SNR）为0、5、10、15、20 dB下利用CAE和β-VAE网络分别做预训练，迭代100次。再利用测试集分别在信噪比为0、5、10、15、20 dB下做交叉测试实验。

当SNR为10 dB时，两种算法收敛效果对比图如图2所示。本文算法收敛速度和收敛结果均优于CAE算法，可以快速完成对信号通道幅相误差的校正。在不同信噪比下训练的网络重构结果的均方误差进行SVM分类时，发现极化相差0°（线极化）、1.5°（椭圆极化）的分类均方误差值稳定，故分类时采用对重构均方误差预设硬阈值的方式。

固定SNR为10 dB，记录两种极化相差相近信号随\beta值变化的匹配率，如表1所示。\beta=0.5时分类效果最好，明显优于经典VAE算法（\beta=1时）。

图3为\beta分别取0.5和1时对于重构均方误差的分类散点图。从图3可以看出，\beta=0.5比\beta=1时，样本远离预设分类界线\mathrm{0.89}\times {\mathrm{10}}^{-\mathrm{5}}。

在SNR为0、5、10、15、20 dB这5种信噪比情况下，β-VAE算法与CAE算法消除通道时变幅相误差后的匹配率对比如图4所示。由图4可知，尤其在低SNR的情况下，β-VAE算法消除幅相误差效果明显优于CAE算法。对比卷积自编码器网络处理一维信号，β-VAE算法实现了SNR为10 dB、极化相差分辨率为1.5°时，小样本100轮学习训练下，匹配率由80.4%提升至99.2%的高精度参数估计。

4.2 OGJS算法优化效果分析

图5为SNR为10 dB下仿真到达角分别为86.4°、94.7°、102.5°这3个DOA来波信号下改进前后的OGJS算法空间谱图。由图5可知，改进后的OGJS算法峰值更为明显。

经过100次蒙特卡洛仿真，改进前后的OGJS算法DOA估计均方误差随SNR的变化如图6所示。由图6可知，改进后的OGJS算法DOA估计均方误差明显减小，当信噪比在0~10 dB之间时尤其明显。当SNR为5 dB时，DOA估计平均绝对误差由3.83°减少至0.27°，从而体现了改进后算法的优越性。

5 结束语

由实验验证可知，利用信号空间电场矢量的变化信息——Lissajous图形，可以对传统来波信号DOA测量算法和极化参数联合估计算法进行改进，将深度学习用于信号射频端处理、误差校正。相对大多数传统算法将关注点放在信号的一维特征，利用Lissajous图形获取信号更高维度的信息，对于估计DOA和极化参数这种高维参数提供新的手段和思路。经过对比近年提出的先进估计算法，体现出所提方法的实用性和优越性。但由于本文VAE网络中采取规定阈值判决，泛化性能上还有待提升，是下一步工作的改进方向。

参考文献

[1] 杨永杰.辐射源极化信息的检测和识别技术研究［D］.西安：西安电子科技大学，2012：57-68.

[2] 李永祯，苏翔，张燕琴，等.雷达辐射源的极化特征及全极化侦收方法［J］.航天电子对抗，2013，29（1）：45-48.

[3] 肖雄，刘桐辛，姜春华，等.基于三正交电偶极子天线的MUSIC估计算法［J］.现代电子技术，2022，45（23）：1-7.

[4] 逯岩斌，陈文东，杨赟秀，等.基于互质极化敏感阵列的参数降维估计算法［J］.电子测量技术，2024，47（5）：173-180.

[5] 郭玉华，常青美，余道杰，等.一种改进的极化域-空域联合的自适应波束形成算法［J］.电子学报，2012，40（6）：1279-1283.

[6] 苑清扬.相干辐射源极化-DOA联合估计算法研究［D］.哈尔滨：哈尔滨工业大学，2021：17-20.

[7] 郑轶松，陈伯孝，杨明磊.一种基于空域滤波的空间临近相干源角度估计方法［J］.电子与信息学报，2016，38（12）：3100-3106.

[8] CHALISE B K， ZHANG Y D， HIMED B. Compressed sensing based joint DOA and polarization angle estimation for sparse arrays with dual-polarized antennas［C］∥Proceedings of the 2018 IEEE Global Conference on Signal and Information Processing. Piscataway， USA： IEEE， 2018：251-255.

[9] WAGNER M， PARK Y， GERSTOFT P. Gridless DOA estimation and root-MUSIC for non-uniform linear arrays［J］. IEEE Transactions on Signal Processing， 2021，69：2144-2157.

[10] 肖思棋.压缩稀疏阵列下基于离网压缩感知的DOA估计算法研究［D］.哈尔滨：哈尔滨工程大学，2022：23-40.

[11] 王强，李佳，沈毅.压缩感知中确定性测量矩阵构造算法综述［J］.电子学报，2013，41（10）：2041-2050.

[12] LI T， LYU J L， WANG R， et al. Bi-directional long short-term memory-based amplitude and phase prediction of antenna array［J］. IEEE Antennas and Wireless Propagation Letters， 2025，24（7）：2034-2038.

[13] 张梓轩，齐子森，许华，等.基于自编码器的阵列时变幅相误差校正算法［J］.西北工业大学学报，2023，41（6）：1134-1145.

[14] 徐定藩.李萨如图形的相位参量［J］.物理与工程，2001，11（5）：27-29.

[15] 郑俊达.基于李萨如图形的相位测量及实现［J］.电子世界，2012（3）：80-81.

[16] 侯青松，郭英，王布宏，等.共形阵列天线振动条件下稳健的DOA估计及位置误差校正［J］.信号处理，2010，26（11）：1756-1760.

[17] XU L， CAO M Y， SONG B Y， et al. Open-circuit fault diagnosis of power rectifier using sparse autoencoder based deep neural network［J］. Neurocomputing， 2018，311：1-10.

[18] FABIUS O， VAN AMERSFOORT J R. Variational recurrent auto-encoders［DB/OL］. （2015-06-15）［2025-03-03］.

[19] WANG Y， LEUS G， PANDHARIPANDE A. Direction estimation using compressive sampling array processing［C］∥Proceedings of the 2009 IEEE/SP 15th Workshop on Statistical Signal Processing. Piscataway， USA： IEEE， 2009：626-629.

[20] TAN Z， YANG P， NEHORAI A. Joint sparse recovery method for compressed sensing with structured dictionary mismatches［J］. IEEE Transactions on Signal Processing， 2014，62（19）：4997-5008.