1.   引言

无源雷达是指通过第三方非合作辐射源信号进行探测的被动雷达系统，自身不发射电磁波，相对传统主动(有源)雷达具有如下优势：(1)无需构建发射机，成本低，尺寸小；(2)无需特定频率分配，无电磁污染；(3)具有良好的四抗能力；(4)辐射源覆盖广泛，易于信息融合、雷达组网以及多站定位。

近年来，基于民用通信信号的无源雷达吸引了国内外的广泛研究^[1]。因其辐射源几乎无死角布设全国，具有覆盖范围广、运行成本低、低空检测性能好、反隐身性能优等优势，如基于第2/3/4代移动电话通信信号(Global System for Mobile communications, GSM/Code Division Multiple Access, CMDA/Long Term Evolution, LTE)^[2-4]、数字音频广播(Digital Audio Broadcasting, DAB)^[5]、数字电视信号(Digital Video Broadcasting-Terrestrial, DVB-T)^[5]、无线局域网信号(Wireless Fidelity, Wi-Fi)^[6]的无源雷达等。绝大多数通信信号的调制方式决定其映射到物理层后呈高斯分布，如，DAB信号、DVB-T信号、LTE信号下行链路、Wi-Fi(IEEE 802.11a/g/n)信号使用正交频分复用(Orthogonal Frequency Division Multiplexing, OFDM)调制、CDMA信号使用码分多址(Code Division Multiple Access, CDMA)调制，根据大数定律，OFDM调制信号在子载波数较多时、CDMA调制信号在扩频码数较多时均呈高斯分布^[7,8]。

无源雷达通过主通道与参考通道分别接收信号。通常情况下，主通道除接收目标回波外还不可避免地接收辐射源直达波及其多径，参考通道天线对准辐射源接收直达波，主通道接收信号须经杂波对消^[9-11]后与参考信号进行互模糊运算获取目标信息。然而，对基于民用通信信号的无源雷达来说，辐射源分布密集，且同频辐射源信号无法在频域区分，如：采用蜂窝布站方式且多个基站共用同一频段的CDMA、LTE信号基站、密集分布且仅有2.4 GHz与5 GHz两个频段的Wi-Fi信号源、工作于同一频段且空间相近但广播内容不同的地面数字多媒体广播(Digital Terrestrial Multimedia Broadcast, DTMB)电台等。因此，主通道还可能接收同频辐射源信号^[12]，参考通道也会不同程度地受同频干扰污染。在本文中，规定目标回波来自主辐射源，其余辐射源均为干扰辐射源。对主通道而言，干扰信号能量通常比目标回波高几个数量级，需对其进行抑制才能有效检测目标，否则将会抬高底噪，掩盖目标，造成漏警。但通过已有方法如最小均方(Least Mean Square, LMS)、递推最小二乘(Recursive Least Square, RLS)、拓展相消(Extensive Cancellation Algorithm, ECA)类算法^[9-11]只能抑制主辐射源直达波及其多径，对同频干扰无能为力；对参考通道而言，虽干扰信号能量比直达波低1个数量级以上，但参考通道信号是互模糊运算的基准信号，若混入同频干扰不仅会造成底噪无法通过对消有效下降，还可能形成相关峰造成虚警，需对参考通道信号进行信号提纯，但现有提纯方法^[13,14]只能针对参考通道未被同频干扰污染的情况。

现有无源雷达同频干扰抑制方法一般仅针对主通道接收信号，未考虑参考通道混入同频干扰的情况。文献[12]对各种干扰峰的形成机理及特点进行了分析，提出了一种信号域帧头预处理的去同频干扰峰方法，但该方法仅能消除相关峰虚警，不能有效降低底噪减小漏警，还会导致目标信噪比降低。文献[15]提出了一种基于自适应波束形成的干扰抑制方法，使用多天线接收信号，在感兴趣的方位角范围内选取相邻天线接收信号的多个子集进行空间平滑波束形成，通过在干扰来波方向上形成零陷抑制干扰。但一方面，当干扰来波与目标回波方向在空域相差不大时，便无法区分；另一方面，天线自由度限制了可形成零陷的数目，不能抑制所有方向的干扰。文献[16]提出了一类基于独立分量分析(Independent Component Analysis, ICA)的同频干扰抑制方法，将各辐射源直达波及其多径、目标回波分别视作一个源信号，使用基于瞬时混合模型的快速独立分量分析 (FastICA) 算法依据源信号的统计独立性与非高斯性实现时域分离进而执行后续处理。但一方面，算法需求天线自由度不少于源信号数，实际环境中多径数目很大，天线自由度与运算能力均不满足分离条件，基于瞬时混合模型的算法将会造成欠定情况严重影响性能；另一方面，民用通信信号几乎都不满足FastICA算法的非高斯性分离条件。

综上所述，现有无源雷达信号处理流程未考虑实际情况下参考通道接收的同频干扰信号，且存在天线硬件要求高、无法应对目标干扰源空域混叠、需求信号非高斯统计特性的问题。本文针对上述问题，在不改变现有硬件条件的情况下提出了一种符合系统实际接收情况、加入同频干扰抑制的信号处理流程。改进流程首先对所有(主、参考)通道信号联合处理，使用多通道盲反卷积算法估计各个辐射源直达波，再利用参考通道信号能量绝大部分来自主辐射源的事实，使用最大互相关法识别主辐射源直达波作为参考信号，其余为各干扰辐射源直达波，然后对主通道信号中各辐射源杂波进行逐个对消，最后用参考信号(主辐射源直达波)与剩余信号进行互模糊运算，进行目标检测。改进流程将所有通道联合处理，更符合系统实际接收情况，最大化地利用了硬件自由度，简化了直达波提纯步骤；需求通道总数目大于等于辐射源数目，保证了改进流程可以在不改变现有硬件条件下使用；最小互信息准则仅需求辐射源信号统计独立；加入的惩罚正则项与带动量项的批量梯度下降优化算法在减少计算量的同时保证了算法的收敛。仿真分析与实测数据验证说明了改进流程可以有效提升对消比、降低底噪、减少漏警，为基于民用通信信号的无源雷达同频干扰抑制问题提供了一种处理框架。

2.   接收信号模型

假定主通道数目为${N_\rm s}$，每个主通道接收各辐射源直达波、多径及目标回波信号。约定目标回波来自主辐射源，其余辐射源均为干扰辐射源。由于目标反射干扰辐射源信号能量微弱，对系统模型影响小，视作系统噪声的一部分。第$k$个$\left( {k = 1,2,·\!·\!· ,{N_{\rm{s}}}} \right)$主通道接收信号可表示为

其中，${a_{i,1,k}}$,${\tau _{i,1,k}}$为第$k$个通道接收主辐射源直达波与多径的衰减系数与离散化时延，${b_{i,q,k}}$,${\tau _{i,q,k}}$为第$k$个通道接收第$q$个干扰辐射源直达波与多径的衰减系数与离散化时延，${c_{i,k}}$, ${\tau '\!\!_{i,k}}$, ${f_{i,k}}$为第$k$个通道接收目标回波的衰减系数、离散化时延与多普勒频率，${s_1}\left( {n - {\tau _{i,1,k}}} \right)$为第$k$个通道接收主辐射源直达波与多径，${s_q}\left( {n - {\tau _{i,q,k}}} \right)$为第$k$个通道接收第$q$个干扰辐射源直达波与多径，${s_1}\left( {n - {{\tau '}\!\!_{i,k}}} \right){{\rm{e}}^{{{{\rm{j2{\text{π}} }}{f_{i,k}}n} / {{f_{\rm{s}}}}}}}$为第$k$个通道接收的目标回波信号，${f_{\rm{s}}}$为采样率，${N_{\rm{a}}}\left( k \right)$, ${N_{\rm{b}}}\left( k \right)$为第$k$个通道接收主、干扰辐射源多径数目，${N_{\rm{t}}}\left( k \right)$为第$k$个通道接收目标回波数目，$Q$为辐射源数目，$e\left( {k,n} \right)$为第$k$个通道接收噪声。

混入同频干扰的参考通道同主通道类似，接收信号可表示为

其中，${a_{i,1}}$,${\tau _{i,1}}$为参考通道接收主辐射源直达波与多径的衰减系数与离散化时延，${b_{i,q}}$,${\tau _{i,q}}$为接收第$q$个干扰源直达波与多径的衰减系数与离散化时延，${s_1}\left( {n - {\tau _{i,1}}} \right)$为参考通道接收主辐射源直达波与多径，${s_q}\left( {n - {\tau _{i,q}}} \right)$为参考通道接收第$q$个干扰辐射源直达波与多径，${N_{\rm{a}}}$,${N_{\rm{b}}}$为接收到主、干扰辐射源的多径数目，$Q$为辐射源数目，$e\left( n \right)$为接收噪声。

3.   加入同频干扰抑制的无源雷达信号处理流程

传统无源雷达信号处理流程如图1(a)。

图  1  无源雷达信号处理流程

Figure  1.  Signal processing flows of passive radar

下载: 全尺寸图片 幻灯片

由式(1)、式(2)分析可知，除噪声外所有通道接收信号均可看作各辐射源直达波的卷积混合，但参考通道中主辐射源信号能量占比最大。因而将所有通道接收信号联合处理估计各辐射源直达波，再利用通道间能量占比差异识别主辐射源信号作为参考信号，进行后续处理。经改进后加入同频干扰抑制的信号处理流程如图1(b)。

改进流程如下：首先将所有通道接收信号联合处理，使用多通道盲反卷积算法估计各辐射源直达波；由于盲反卷积输出信号具有顺序不确定性^[17]，利用参考通道中主辐射源信号能量占比最大，判定与参考通道信号互相关最大的信号为主辐射源直达波，其余为各干扰辐射源直达波；然后利用已有杂波对消算法^[9-11]抑制主通道中各辐射源直达波与多径，剩余目标回波信号；最后将剩余信号与提取出的主辐射源直达波进行互模糊运算，检测目标。

现有杂波对消算法^[9-11]已较为成熟，非本文核心，不作赘述。本节重点介绍基于卷积混合模型的多通道盲反卷积算法。

3.1   基于卷积混合模型的多通道盲反卷积算法

式(1)、式(2)可用卷积混合模型统一表达。卷积混合模型是一个离散线性时不变(Linear Time-Invariant, LTI)多输入多输出(Multiple-Input Multiple-Output, MIMO)有限脉冲响应(Finite Impulse Response, FIR)滤波系统，第$p$个通道接收信号可表示为

其中，$P = {N_{\rm{s}}} + 1$为通道数目，$K - 1$为信道滤波器最大阶数，$Q$为辐射源数目，${e_p}\left( n \right)$为第$p$个通道接收噪声。式(3)可写作

其中，${\text{H}}\left( z \right)$为$P \times Q$滤波系统传递函数矩阵，${\text{H}}\left( z \right)$第$\left( {i,j} \right)$个元素${h_{ij}}\left( z \right)$为第$j$个辐射源信号到第$i$个接收通道信号的传递函数，${\text{x}}\left( n \right) = \left[ {{x_1}\left( n \right)}\right.$$\left.{{x_2}\left( n \right)}\ ·\!·\!· \ {{x_P}\left( n \right)} \right]^{\rm{T}}$为通道接收信号向量，${\text{s}}\left( n \right) =$${\left[ {{s_1}\left( n \right)}\ {{s_2}\left( n \right)} ·\!·\!· \ {{s_Q}\left( n \right)}\right]^{\rm{T}}}$为辐射源信号向量，${\text{e}}\left( n \right) = {\left[ {{e_1}\left( n \right)}\ {{e_2}\left( n \right)}\ ·\!·\!· \ {{e_P}\left( n \right)} \right]^{\rm{T}}}$为噪声信号向量。

多通道盲反卷积的任务可表述为：在源信号与信道参数完全未知的情况下，通过使代价函数达到最优求解一个离散线性时不变多输入多输出有限脉冲响应滤波系统，使得接收信号通过系统后为源信号的估计。即，求解传递函数矩阵${\text{W}}\left( z \right)$，使得

为${\text{s}}\left( n \right)$的估计。其中，$L$为分离滤波器长度。多通道盲反卷积算法要求接收信号维度大于等于源信号维度即$P \ge Q$，其算法示意图如图2，其中$\Sigma$表示求和。

图  2  多通道盲反卷积算法示意图

Figure  2.  Schematic diagram of multi-channel blind deconvolution algorithm

下载: 全尺寸图片 幻灯片

多通道盲反卷积一般很难得到源信号的精确估计，其完成标准为估计信号${\text{y}}\left( n \right)$为源信号${\text{s}}\left( n \right)$经重新排列、尺度缩放、时延后的版本^[17]，即

其中，${\text{W}}\left( z \right)$为分离系统传递函数矩阵，${\text{G}}\left( z \right)$为全局系统传递函数矩阵，${\text{P}}$为一置换矩阵，${\text{Λ}}$为一非奇异对角尺度变换矩阵，${\text{D}}\left( z \right) = {\rm diag}\left( {{z^{ - {\tau _1}}},{z^{ - {\tau _2}}}, ·\!·\!· ,{z^{ - {\tau _P}}}} \right)$为一时延对角矩阵。根据自适应滤波理论，源信号的重新排列、尺度缩放和时延不会对对消产生影响^[18]，保证了后续处理的可行性。

考虑到源信号的高斯性，选取互信息作为代价函数^[19]，该代价函数仅依赖于源信号间统计独立，分离准则为输出信号间互信息最小化。互信息定义为

其中，${p_{\text{y}}}\left( {\text{y}} \right)$为${\text{y}}$的联合概率密度函数，${p_{{y_i}}}\left( {{y_i}} \right)$为${\text{y}}$各分量的边缘概率密度函数。互信息$I\left( {\text{y}} \right) \ge {\rm{0}}$具有非负性，取等号时当且仅当${p_{\text{y}}}\left( {\text{y}} \right) = \displaystyle\prod\nolimits_i {{p_{{y_i}}}\left( {{y_i}} \right)}$，此时${\text{y}}$的联合概率密度函数等于各分量边缘概率密度函数的乘积，各分量统计独立。实际环境中接收信号为随机过程，其独立性判据为

对任意${t_1},{t_2}, ·\!·\!· ,{t_P} \in {\text{T}}$都成立，其中${\text{T}}$为一实参数集。即，统计独立条件为分量在每个时刻独立。由于信号源本身相互统计独立，导致接收信号分量间相关的原因是滤波系统${\text{H}}\left( z \right)$引起的卷积混合。所以要使得信号分量间统计独立，只需令其在可能导致相关的时延集合$\left\{ { - T, - T + 1, ·\!·\!· ,T - 1,T} \;\right\}$内消除相关，这里取$T = 2L$。于是得到代价函数

其中，${{\text{y}}^{\left( {{{\text{h}}_i}} \right)}}\left( n \right) = \left[ {{y_1}\left( {n - {h_{i,1}}} \right)}\ {{y_{\rm{2}}}\left( {n - {h_{i,{\rm{2}}}}} \right)}\ ·\!· · \right.$$\left.{{y_P}\left( {n - {h_{i,P}}} \right)} \right]^{\rm{T}}$, ${{\text{h}}_i}$为任一可能导致信号分量相关的时延序列即${{\text{h}}_i} = \left[ {{h_{i,1}}}\ {{h_{i,{\rm{2}}}}}\ ·\!·\!· \ {{h_{i,P}}} \right]$, ${h_{i,j}} \in \left\{ - T, \right.$$\left. {- T + 1,·\!·\!· ,T - 1,T} \,\right\}$, ${\text{κ}}$为所有不重复的${{\text{h}}_i}$组成的集合，${{\text{y}}^{\left( {{{\text{h}}_i}} \right)}}$为${{\text{h}}_i}$对应的信号向量时延版本。

然而，只保证互信息最小化非保证收敛到最优解的充分条件，还可能收敛到两类非最优解：${\text{y}}$中分量为互相独立的高斯噪声、${\text{y}} = {\text{0}}$。为规避这一问题，引入含惩罚正则项的代价函数

其中，${\left\| \cdot \right\|_2}$为范2-数，$\lambda > 0$和$\gamma > 0$为实常数。式(10)中第2项借鉴了最小失真原理^[20]的相关思想：在所有可行解中规定一个最优解的结构，使得输出${\text{y}}\left( n \right)$尽可能为输入${\text{x}}\left( n \right)$的主要能量部分，这里是使得输出分量尽可能为占接收信号大部分能量的辐射源直达波，从而避免收敛至第1类非最优解。第3项即归一化项保证输出不会收敛到第2类非最优解，也避免了算法发散。此外，式(10)中最后一项中含多个期望运算，不易求解。为简化运算，对${\text{x}}\left( n \right)$进行去均值、归一化的预处理，式(10)代价函数便可去掉第3项中${\rm{E}}\left( {{y_i}\left( n \right)} \right)$项。

综上所述，含惩罚正则项的多通道盲反卷积算法可表达为

3.2   互信息梯度及其估计方法

算法使用梯度下降迭代求解式(5)所示矩阵形式滤波系统参数，需研究互信息梯度及其估计方法。

直接计算时延版本信号互信息梯度较为困难，首先计算无时延版本信号互信息梯度。假设接收信号为均方连续随机过程(实际中容易取得)。用${w_{k,gh}}$表示${{\text{W}}_k}$的第$\left( {g,h} \right)$个元素，易得${y_1}\left( n \right),{y_2}\left( n \right),$$·\!·\!· ,{y_P}\left( n \right)$中仅有${y_g}\left( n \right)$依赖于${w_{k,gh}}$。根据均方连续随机过程期望与微分运算次序的可交换性，有

所以，

其中，$P \times {\rm{1}}$向量${{\text{K}}_{\text{y}}}\left( n \right) = \left[ {{k_{{\text{y}},1}}\left( n \right)}\ {{k_{{\text{y}},{\rm{2}}}}\left( n \right)}\ ·\!·\!· \ \right.$$\left.{{k_{{\text{y}},P}}\left( n \right)} \right]^{\rm{T}}$，它的第g个元素为

下面对时延版本信号向量的互信息梯度进行计算。类似式(12)，有

式(15)最后一步利用了${\text{x}}\left( n \right)$与${\text{y}}\left( n \right)$的短时平稳性。因此

其中，$P \times {\rm{1}}$向量${\text{K}}_{{{\text{y}}^{\left( {{{\text{h}}_i}} \right)}}}^{\left( { - {{\text{h}}_i}} \right)}\left( n \right) = \left[ {{k_{{{\text{y}}^{\left( {{{\text{h}}_i}} \right)}},1}}\left( {n + {h_{i,1}}} \right)}\ \ {{k_{{{\text{y}}^{\left( {{{\text{h}}_i}} \right)}},{\rm{2}}}}\left( {n + {h_{i,{\rm{2}}}}} \right)}\ ·\!·\!· \ {{k_{{{\text{y}}^{\left( {{{\text{h}}_i}} \right)}},P}}\left( {n + {h_{i,P}}} \right)}\right]^{\rm{T}}$。

下面介绍互信息梯度估计方法。考虑到${\text{K}}_{{{\text{y}}^{\left( {{{\text{h}}_i}} \right)}}}^{\left( { - {{\text{h}}_i}} \right)}\left( n \right)$为${{\text{K}}_{{{\text{y}}^{\left( {{{\text{h}}_i}} \right)}}}}\left( n \right)$的时延版本，${{\text{y}}^{\left( {{{\text{h}}_i}} \right)}}\left( n \right)$为${\text{y}}\left( n \right)$的时延版本，因此只需研究${{\text{K}}_{\text{y}}}\left( n \right)$的估计方法。${{\text{K}}_{\text{y}}}\left( n \right)$的第$g$个元素为

于是问题可简化为估计$p\left( {{{\text{y}}_g}\left( n \right)|{y_g}\left( n \right)} \right)$与${p^{\left( {\rm{1}} \right)}}\left( {{{\text{y}}_g}\left( n \right)|{y_g}\left( n \right)} \right)$，其中，${{\text{y}}_g}\left( n \right)$为$P - 1$维变量

这里假定$p\left( {{{\text{y}}_g}\left( n \right)|{y_g}\left( n \right)} \right)$对任意${{\text{y}}_g}\left( n \right)$与${y_g}\left( n \right)$都是连续的且${p^{\left( {\rm{1}} \right)}}\left( {{{\text{y}}_g}\left( n \right)|{y_g}\left( n \right)} \right)$存在。容易得到

其中，$\delta \left( {\text{y}} \right)$是多维狄拉克$\delta$函数(Dirac delta function)。将${\rm E}\left\{ {\delta \left( {{{\text{y}}_g}\left( n \right) - {{\text{y}}_g}\left( {{n_{\rm{0}}}} \right)} \right)|{y_g}\left( n \right)} \right\}$视作${y_g}\left( n \right)$的函数，将${\rm E}\left\{ {\delta \left( {{{\text{y}}_g}\left( n \right) - {{\text{y}}_g}\left( {{n_{\rm{0}}}} \right)} \right)|{y_g}\left( n \right)} \right\}$在${y_g}\left( n \right) = {y_g}\left( {{n_0}} \right)$泰勒展开

因此，可利用式(20)同时估计$p\left( {{{\text{y}}_g}\left( n \right)|{y_g}\left( n \right)} \right)$与${p^{\left( {\rm{1}} \right)}}\left( {{{\text{y}}_g}\left( n \right)|{y_g}\left( n \right)} \right)$。但$\delta \left( {\text{y}} \right)$非严格意义上的函数而是广义函数。为了使用式(20)一般采用序列函数的极限对$\delta \left( {\text{y}} \right)$进行建构。假设${\text{y}}$为$P - 1$维变量，$\eta \left( {\text{y}} \right)$为一个在${\text{R}}$上总积分为1的绝对可积函数，并定义

这里，$\varepsilon$定义为${\eta _\varepsilon }\left( {\text{y}} \right)$的带宽。$\delta \left( {\text{y}} \right)$为函数序列${\eta _\varepsilon }\left( {\text{y}} \right)$的极限

于是，

将式(23)运用到式(20)中有

在工程中通常使用带宽$\varepsilon$足够小的核函数${\eta _\varepsilon }\left( {\text{y}} \right)$及有限阶泰勒展开使用式(24)进行估计，问题转化为如下最优化问题

其中，$N$为估计样本数，${K_b}\left( {\text{y}} \right) = {b^{1 - P}}{\kappa _1}\left( {\dfrac{{\text{y}}}{b}} \right)$, ${W_h}\left( x \right) = {h^{ - 1}}{\kappa _2}\left( {\dfrac{x}{h}} \right)$为核函数，$b$,$h$分别其带宽，本文采用高斯核函数${\kappa _1}\left( {\text{y}} \right) = \dfrac{1}{{{{\left( {2{\text{π}} } \right)}^{{{\left( {P - 1} \right)} / 2}}}}}\exp \left( { - \dfrac{{{{\text{y}}^{\rm T}}{\text{y}}}}{2}} \right)$与${\kappa _2}\left( x \right) = \dfrac{1}{{\sqrt {2{\text{π}} } }}\exp \left( { - \dfrac{{{x^2}}}{2}} \right)$，估计参数${\text{θ}} = {\left[ {{\theta _0}}\ {{\theta _{\rm{1}}}}\ ·\!·\!· \ {{\theta _r}} \right]^{\rm{T}}} = \left[ {\hat p\left( {{{\text{y}}_g}\left( {{n_0}} \right)|{y_g}\left( {{n_0}} \right)} \right)}\right.\ \ {{\hat p}^{\left( {\rm{1}} \right)}} \left( {{{\text{y}}_g}\left( {{n_0}} \right)|}\right.$$\left.\left.{y_g}\left( {{n_0}} \right) \right)\ ·\!·\!· \ {{{\hat p}^{\left( r \right)}}\left( {{{\text{y}}_g}\left( {{n_0}} \right)|{y_g}\left( {{n_0}} \right)} \right)} \right]^{\rm{T}}$。使用最小二乘法求解式(24)，参数${\text{θ}}$的最优估计值为

其中，${\text{W}} \! =\! {\rm diag}\left[ {{W_h}\left( {{y_g}\left( {{n_{\rm{1}}}} \right) \!-\! {y_g}\left( {{n_0}} \right)} \right)}\ {W_h}\left( {{y_g}\left( {{n_{\rm{2}}}} \right) -} \right.\right.$$\left. \left.{{y_g}\left( {{n_0}} \right)} \right) ·\!·\!· {{W_h}\left( {{y_g}\left( {{n_N}} \right) - {y_g}\left( {{n_0}} \right)} \right)} \right]$, ${\text{Y}} \!=\! \left[ {{K_b}} \right.\left( {{{\text{y}}_g}} \right.\left( {{n_1}} \right) -$${\left. { {{\text{y}}_g}\left( {{n_{\rm{0}}}} \right)} \right)}\;\;\, {{K_b}\left( {{{\text{y}}_g}\left( {{n_{\rm{2}}}} \right) - {{\text{y}}_g}\left( {{n_{\rm{0}}}} \right)} \right)}\;\;\; ·\!·\!· \;\;\; {K_b}\left( {{\text{y}}_g}\left( {{n_N}} \right) \right.$–$\left.\left.{{\text{y}}_g}\left( {{n_{\rm{0}}}} \right) \right) \right]^{\rm{T}}$, ${\text{X}}$为$N \times \left( {r + 1} \right)$的矩阵且其第$i$行为$\left[ {\rm{1}}\ {\left( {{y_g}\left( {{n_i}} \right) - {y_g}\left( {{n_0}} \right)} \right)}\ ·\!·\!· \ {{{{{\left( {{y_g}\left( {{n_i}} \right) - {y_g}\left( {{n_0}} \right)} \right)}^r}}/ {r!}}} \right]$。

综上所述，${\text{K}}_{{{\text{y}}^{\left( {{{\text{h}}_i}} \right)}}}^{\left( { - {{\text{h}}_i}} \right)}\left( n \right)$的估计流程可表示为式(27)。

3.3   惩罚正则项的梯度

正则项第1项的梯度为

第2项展开计算梯度，最后将$E\left( {y_i^2\left( n \right)} \right)$并入期望运算有

其中，$\delta \left( {g - i} \right) = \left\{ \begin{gathered} 1,g = i \\ 0,g \ne i \\ \end{gathered} \right.$为克罗内克$\delta$函数。因此

其中，${{\text{Γ}}_i}\left( n \right) \;= \left[ { {4\, \rm E}}\left\{\, \left[\, {{\rm E}\left( {y_i^2\left( n \right)} \right) \,-\, 1} \right] {y_i}\left( n \right){x_1}\right.\right. ( n$$\left. - k ) \right\}\, {{\rm{4E}}\left\{ {\left[ {{\rm E}\left( {y_i^2\left( n \right)} \right) - 1} \right]{y_i}\left( n \right){x_2}\left( {n - k} \right)} \right\}\; ·\!· · {\rm{4E}}}$$\left.\left\{ \left[ {{\rm E}\left( {y_i^2\left( n \right)} \right) } - 1\right]{y_i}\left( n \right){x_P}\left( {n - k} \right) \right\} \right]$为梯度矩阵第$i$行。最后得到梯度矩阵

其中，${\text{Φ}}\left( n \right)\! =\! \left[ {{\rm{4}}\left[ {{\rm{E}}\left( {y_1^2\left( n \right)}\! \right) - 1} \right]{y_1}\left( n \right)}\,{\rm{4}}\left[ {\rm{E}}\left( {y_2^2}\left( n \right) \right)\right.\right.$$\left.\left.{- 1} \right]{y_2}\left( n \right)\ ·\!·\!·\ {4\left[ {{\rm{E}}\left( {y_P^2\left( n \right)} \right) - 1} \right]{y_P}\left( n \right)} \right]^{\rm{T}}$。

3.4   加入动量项的批量梯度下降优化算法

综合3.1—3.3节相关分析，得到代价函数的梯度

直接使用式(32)迭代需要昂贵的计算成本。为解决此问题并保证收敛速度，采用随机批量梯度下降(Batch Gradient Descent, BGD)与动量梯度下降(gradient descent with momentum)相结合的方式进行迭代求解。改进方法每次迭代对${{\text{W}}_k}\left( {k = 0,1, }\right.$$\left.{·\!·\!·,L - 1} \right)$进行更新，计算梯度时随机选取${\text{κ}}$中个数为${\rm{batchsize}}$的向量元素生成多个时延版本信号向量，将其互信息之和替代$J$的第1项，即

此外，在更新量中额外加入1个动量项${{\text{υ}}_k}$, ${{\text{υ}}_k}$的初始值为${{\text{υ}}_k} = {\text{0}}$，每一步迭代公式为

加入动量项的作用是避免因使用部分梯度代替原代价函数而可能带来的局部震荡。用变量的上标表示迭代的步数，实际上将${\text{υ}}_k^{[t]}$展开有

式(35)表明，梯度的计算使用了指数加权平均方法，每一步计算梯度与所有历史信息都有关联，对信息的保留量随时间指数变化，越近时间的历史信息保留地越多。不仅如此，真正的梯度下降方向也能很好地保持。

综上所述，多通道盲反卷积算法步骤如表1。

表  1  多通道盲反卷积算法步骤

Table  1.  Algorithm procedure of multi-channel blind deconvolution

参数：${\text{κ}},{\rm{batchsize}},L,\alpha ,\beta ,b,h,r$ 　输入：${\text{x}}\left( n \right)$ 　初始化：${{\text{W}}_{\rm{0}}} = {\text{I}},{{\text{W}}_k} = {\text{0}}\left( {k = 1,2, ·\!·\!· ,L - 1} \right),{\text{y}}\left( n \right) = {\text{x}}\left( n \right),{{\text{υ}}_k} = {\text{0}}$ 　循环： 　(1) 随机选取${{\text{q}}_1},{{\text{q}}_2}, ·\!·\!· ,{{\text{q}}_{{\rm{batchsize}}}} \in {\text{κ}}$； 　(2) 对每个${{\text{h}}_i} \in \left\{ {{{\text{q}}_1},{{\text{q}}_2}, ·\!·\!· ,{{\text{q}}_{{\rm{batchsize}}}}} \right\}$，估计${\text{K}}_{{{\text{y}}^{\left( {{{\text{h}}_i}} \right)}}}^{\left( { - {{\text{h}}_i}} \right)}\left( n \right)$； 　(3) 对$k = 0{\rm{ ,1,}} ·\!·\!· {\rm{,}}L - 1$，更新${{\text{υ}}_k}$,${{\text{W}}_k}$: 　　 ${{\text{υ}}_k} \leftarrow \beta {{\text{υ}}_k} + \left( {1 - \beta } \right){\rm{E}}\left\{ {\left[ {\sum\nolimits_{{{\text{h}}_i} \in \left\{ {{{\text{q}}_1},{{\text{q}}_2}, ··· ,{{\text{q}}_{{\rm{batchsize}}}}} \right\}} {{\text{K}}_{{{\text{y}}^{\left( {{{\text{h}}_i}} \right)}}}^{\left( { - {{\text{h}}_i}} \right)}\left( n \right) + 2\left( {{\text{y}}\left( n \right) - {\text{x}}\left( n \right)} \right) + {\text{Φ}}\left( n \right)} } \right]{\text{x}}{{\left( {n - k} \right)}^{\rm{T}}}} \right\}$ 　　 ${{\text{W}}_k} \leftarrow {{\text{W}}_k} - \alpha {{\text{υ}}_k}$; 　(4) 更新${\text{y}}\left( n \right)$:${\text{y}}\left( n \right) \leftarrow \sum\limits_{k = {\rm{0}}}^{L - 1} {{{\text{W}}_k}{\text{x}}\left( {n - k} \right)}$； 　(5) 判断：收敛或达到最大迭代次数后结束循环。

下载: 导出CSV 

| 显示表格

4.   仿真分析

对加入同频干扰抑制的无源雷达信号处理流程进行仿真分析。仿真设置3个同频高斯信号源(信号源1视作主辐射源直达波，信号源2和信号源3视作同频干扰辐射源1, 2直达波)，3个通道(主通道1，主通道2与参考通道)接收信号，满足本文算法对系统自由度的要求。仿真信号源带宽4 MHz，采样率10 MHz，积累时间0.02 s。3通道接收信号成分见表2—表4，各天线接收信噪比为30 dB。

表  2  主通道1接收信号成分

Table  2.  Signal component of primary channel 1

信号源	信号成分
	参数	直达波	多径1	多径2	多径3	强目标1	弱目标2
主辐射源	时延(μs)	0.1	0.3	1.6	2.1	46.1	59.1
	幅度(dB)	–0.41	–19.53	–28.29	–35.09	–33.15	–43.10
干扰辐射源1	时延(μs)	0	0.7	1.2	2.4	19.9
	幅度(dB)	–3.10	–19.27	–20.57	–29.79	–33.98
干扰辐射源2	时延(μs)	0	0.5	1.7	2.6
	幅度(dB)	–2.55	–18.63	–21.75	–33.31

下载: 导出CSV 

| 显示表格

表  3  主通道2接收信号成分

Table  3.  Signal component of primary channel 2

信号源	信号成分
	参数	直达波	多径1	多径2	多径3	强目标1	弱目标2
主辐射源	时延(μs)	0.1	0.8	1.6	2.2	46.1	59.1
	幅度(dB)	–0.16	–20.17	–28.62	–32.52	–30.74	–40.56
干扰辐射源1	时延(μs)	0	0.3	1.5	2.4
	幅度(dB)	–3.28	–19.99	–23.13	–36.34
干扰辐射源2	时延(μs)	0	0.5	1.8	2.3
	幅度(dB)	–2.56	–19.58	–25.93	–30.22

下载: 导出CSV 

| 显示表格

表  4  参考通道接收信号成分

Table  4.  Signal component of reference channel

信号源	信号成分
	参数	直达波	多径1	多径2	多径3
主辐射源	时延(μs)	0	1.1	1.7	2.9
	幅度(dB)	0	–22.05	–28.64	–34.42
干扰辐射源1	时延(μs)	0.2	0.7	1.3	2.8
	幅度(dB)	–10.17	–22.82	–23.99	–33.00
干扰辐射源2	时延(μs)	0.1	0.5	1.5	2.3
	幅度(dB)	–9.37	–24.40	–26.03	–35.21

下载: 导出CSV 

| 显示表格

分别使用传统无源雷达信号处理流程以及本文流程进行目标检测，其中杂波对消均使用归一化最小均方(Normalized LMS, NLMS)算法，对消参数均为：步长0.01，滤波器长度350，归一化分母中正则化参数为1。

图3(a)，图3(b)为盲反卷积前、后3通道信号与3信号源间的散点图。盲反卷积前3通道信号与3信号源散点图均为椭圆形，椭圆越扁长说明相关度越大；由于3通道信号均为3信号源的卷积混合，所以与3信号源间均有一定相关性，且主辐射源直达波与参考通道信号相关度最大，符合仿真设置。经盲反卷积后，输出信号1与主辐射源直达波、输出信号2、输出信号3与干扰辐射源1、干扰辐射源2直达波的散点图呈细长条状(相关性强)，而其余散点图基本呈正圆形(相关性弱)，说明分离信号1～3经过盲反卷积基本已被提纯为主辐射源、干扰辐射源1、干扰辐射源2直达波，消除了因卷积混合引起的相关。

下载: 全尺寸图片 幻灯片

图4(a)、图5(a)为使用传统处理流程的处理结果，图4(a)为互模糊函数距离-多普勒平面，图5(a)为互模糊函数距离剖面。使用传统流程处理，对消比2.04 dB，互模糊函数底噪62.77 dB，强目标幅度73.12 dB，从底噪中显露出来，弱目标被底噪掩盖，造成漏警。

图  4  处理结果(互模糊函数距离-多普勒平面)

Figure  4.  Processing results (range-Doppler plane of the cross-ambiguity function)

下载: 全尺寸图片 幻灯片

图  5  处理结果(互模糊函数距离剖面)

Figure  5.  Processing result (range profile of the cross-ambiguity function)

下载: 全尺寸图片 幻灯片

图4(b)、图5(b)为使用本文处理流程的处理结果，图4(b)为互模糊函数距离-多普勒平面，图5(b)为互模糊函数距离剖面。使用本文流程处理，3次对消总对消比20.60 dB，互模糊函数底噪44.30 dB，强目标幅度72.90 dB，弱目标幅度60.72 dB，从底噪中显露出来。对比传统流程处理结果，对消比提升了18.56 dB，互模糊函数底噪下降了18.47 dB，强目标幅度基本不变，弱目标因底噪下降而显露出来，避免了漏警。

5.   实测数据验证

5.1   实验场景

为验证本文所提信号处理流程的有效性，开展了目标检测外场实验。实验选用中国电信频分双工长期演进(Frequency Division Duplexing-Long Term Evolution, FDD-LTE)信号作为第3方照射源，中心频率为1867.5 MHz，带宽为15 MHz，采样频率50 MHz。实验场景如图6(a)—图6(c)所示，使用双天线接收信号数据，其中左侧天线位于LTE信号基站对面接收参考通道信号，右侧天线指向低空监测区域接收主通道信号，实验设置1个合作目标(汽车)，可以根据实验需要按照设定路线与速度行驶，通过理论结果与两种流程处理结果对比的方式验证本文的正确性与有效性。

图  6  实验场景

Figure  6.  Experiment scenes

下载: 全尺寸图片 幻灯片

5.2   实验结果及分析

使用0.2 s的接收数据进行对比实验，分别使用传统无源雷达信号处理流程以及本文流程进行处理，其中杂波对消均使用NLMS算法，对消参数均为：步长0.05，滤波器长度256，归一化分母中正则化参数为1。

图7(a)、图8(a)为使用传统处理流程的处理结果，图7(a)为互模糊函数三维图，图8(a)为互模糊函数距离剖面。使用传统流程，对消比1.17 dB，几乎没有对消效果。由于通道中同频干扰信号的影响，互模糊函数底噪73.75 dB，目标幅度72.75 dB，目标被高底噪淹没，无法检测，造成漏警。

图  7  处理结果(互模糊函数三维图)

Figure  7.  Processing result (3-D graph of the cross-ambiguity function)

下载: 全尺寸图片 幻灯片

图  8  处理结果(互模糊函数距离剖面图)

Figure  8.  Processing result (range profile of the cross-ambiguity function)

下载: 全尺寸图片 幻灯片

图7(b)、图8(b)为使用本文处理流程的处理结果，图7(b)为互模糊函数三维图，图(8)为互模糊函数距离剖面。使用本文流程，对消比13.35 dB，互模糊函数底噪63.70 dB，目标幅度72.44 dB，从底噪中显露出来，避免了漏警。对比传统流程处理结果，对消比提升了12.18 dB，底噪降低了10.05 dB，大部分杂波被抑制，目标幅度基本不变，最终检测得到目标距离546 m，速度13.03 m/s，与合作汽车目标的实验设置一致，验证了本文所提方法的正确性与有效性。

6.   结束语

本文针对基于民用通信信号的无源雷达主通道与参考通道易同时受同频干扰污染的实际现状，提出一种加入同频干扰抑制的无源雷达信号处理流程。改进流程首先建立针对所有通道接收信号的卷积混合模型，对所有通道联合处理，考虑多数民用通信信号的高斯统计特性，以互信息为代价函数使用盲反卷积算法估计各辐射源直达波，再利用信号能量差异识别主、干扰辐射源直达波，然后对各辐射源杂波进行逐个对消，最后对剩余信号与识别出的主辐射源直达波进行互模糊检测。改进流程将所有通道联合处理，更加符合系统实际接收情况，弥补了硬件自由度，简化了直达波提纯步骤；仅需求天线自由度大于等于辐射源数目，可以在不改变现有硬件接收条件的情况下使用；加入的惩罚正则项与带动量项的批量梯度下降优化算法减少了计算量的同时避免了迭代震荡，保证了算法的收敛。仿真分析及实测数据验证说明算法可有效抑制干扰，提升对消比，降低底噪，减少漏警，为基于民用通信信号的无源雷达同频干扰抑制问题提供了一种有效解决方案，也为无源雷达信号处理流程提供了一种新的思路。