全球卫星导航系统因其导航精度高和不随时间积累误差的优点得到了广泛的应用，然而卫星导航信号到达GNSS接收机时的功率比环境基底噪声还要低20 dB，所以极易受到强干扰信号压制.此时GNSS接收机将无法正常捕获卫星导航信号，必然导致接收机搜星定位失败^[1-2].因此目前普遍在接收机端采用阵列天线抑制强干扰信号.工程中常采用线性约束最小功率准则在强干扰信号来向上自适应生成零陷来抑制干扰^[3]，并采用最小均方算法迭代求解阵列权值.但是在GNSS接收机载体处于高频振动或者高速运动的状态下，干扰信号来向在短时间内会发生快速变化，上述算法可能无法及时收敛，致使干扰信号不能被及时有效地抑制^[4-6]，导致GNSS接收机失锁.

针对这一问题，通过采用基于协方差矩阵锥化的零陷加宽算法可以有效加宽干扰零陷，使得干扰信号变化方向始终保持在零陷之内，从而达到抑制干扰的目的.Mailloux^[7]首先提出了协方差矩阵锥化零陷加宽算法，并假设在干扰源附近均匀分布着若干相同干扰源，求解假设情况下信号协方差矩阵与真实情况下协方差矩阵的关系，再利用真实协方差矩阵得到虚拟协方差矩阵，即可实现干扰零陷拓宽；Zatman^[8]则假设干扰源带宽为真实带宽的若干倍，通过求解虚拟带宽情况下信号协方差矩阵与真实带宽下协方差矩阵的关系来拓宽干扰零陷.Guerci等^[9-10]又系统总结了协方差矩阵锥化算法的一般表达式.在此基础上，李荣锋等^[11]提出了新的零陷加宽抗干扰思路，即假设干扰信号来向变化服从高斯概率分布，实现了干扰零陷拓宽，并且证明了当干扰信号来向变化服从均匀分布时其算法与Zatman算法是等价的；随后，Lu等^[12]提出了干扰信号来向变化服从拉普拉斯分布的零陷加宽方法；武思军等^[13]提出了干扰信号来向变化服从伯努利分布的零陷加宽方法；张柏华等^[14-15]研究了基于空时阵列处理的零陷加宽算法；王妙等^[16]提出了基于降秩共轭梯度法的零陷加宽技术，能够有效降低计算量；柯熙政等^[17]推导了基于均匀圆阵的零陷加宽算法.

需要说明的是，干扰信号来向变化实际上是未知的，并非一定是均匀分布，拉普拉斯分布，抑或是高斯分布.文献[11-13]中借用概率分布统计模型的目的是为了拓展干扰零陷宽度，与Mailloux虚拟干扰源的本质是相同的，只是采用不同概率分布统计模型所形成的零陷形状有所区别，性能也会有差异.正因为干扰来向变化是未知的，干扰零陷深度在零陷宽度范围内更应该保持一致.

本文提出了一种假设干扰信号来向变化服从三角形分布概率统计模型，基于线性约束最小功率准则的协方差矩阵锥化零陷加宽抗干扰算法，并将其与具有代表性的均匀分布、拉普拉斯分布算法进行比较分析，理论分析与仿真实验证明了所提算法拓宽的零陷更加平稳，输出信干噪比也更高，最后结合BD2软件接收机验证了该算法的有效性.

1 信号模型

GNSS接收机阵列天线模型如图 1所示，阵列天线包含M个阵元，阵元之间间距为d，每个阵元对应一个独立的处理通道.远场信号源波面到达阵列天线前端可视为平面波(plane wave)，各阵元接收信号后首先在射频前端(radio frequency front end, RFFE)中对其进行放大去噪，再经模数转换器(analog to digital converter, ADC)抽样提取为数字信号，最后对每一路信号进行加权求和得到输出值.

假设L个卫星导航信号和Q个窄带干扰信号从远场入射，那么阵列天线的接收信号模型可以表示为

$ \mathit{\boldsymbol{x}}\left( t \right) = \sum\limits_{l = 1}^L {\mathit{\boldsymbol{a}}\left( {{\theta _l}} \right){s_l}\left( t \right)} + \sum\limits_{q = 1}^Q {\mathit{\boldsymbol{a}}\left( {{\theta _q}} \right){j_q}\left( t \right)} + \mathit{\boldsymbol{e}}\left( t \right). $

式中：θ_l为第l个卫星导航信号来向; s_l(t)为第l个卫星导航信号包络；θ_q为第q个干扰信号来向; j_q(t)为第q个干扰信号包络；e(t)为噪声向量；a(θ)为信号空域导向矢量，对于半波长均匀线阵，其表达式为[1, e^{-jπdsin θ}, …, e^{-jπ(M-1)dsin θ}].

线性约束最小功率准则是在线性约束条件下，通过使阵列输出功率最小来求解权矢量.工程中常设置约束条件为天线参考阵元无失真，其目标函数为

$ \mathop {\min }\limits_w {\mathit{\boldsymbol{w}}^{\rm{H}}}{\mathit{\boldsymbol{R}}_x}\mathit{\boldsymbol{w}},{\rm{s}}.\;{\rm{t}}.\;{\mathit{\boldsymbol{w}}^{\rm{H}}}{\mathit{\boldsymbol{c}}_M} = 1. $

(1)

式中：上标“H”为Hermit转置；c_M=[1, 0, …, 0]^H为M维约束列向量；R_x为仅包含干扰信号与噪声的理论协方差矩阵.该算法也称功率倒置(power inversion, PI)算法，本文所提算法将基于该准则求解阵列权值.

利用拉格朗日乘数法，设置式(1)的最小化函数为

$ J = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{H}}}{\mathit{\boldsymbol{R}}_x}\mathit{\boldsymbol{w}} + \lambda \left( {{\mathit{\boldsymbol{w}}^{\rm{H}}}{\mathit{\boldsymbol{c}}_M} - 1} \right), $

(2)

式中λ为拉格朗日乘数.为了求解式(2)的最优化问题，求解其关于w的梯度，并令其为0，可得

$ {\nabla _w}J\left( \lambda \right) = {\mathit{\boldsymbol{R}}_x}\mathit{\boldsymbol{w}} + \lambda {\mathit{\boldsymbol{c}}_M} = 0, $

进而可知

$ \mathit{\boldsymbol{w}} = - \lambda \mathit{\boldsymbol{R}}_x^{ - 1}{\mathit{\boldsymbol{c}}_M}. $

(3)

为了求解拉格朗日参数λ，将式(3)代入式(1)中的约束条件，可得

$ {\left( { - \lambda \mathit{\boldsymbol{R}}_x^{ - 1}{\mathit{\boldsymbol{c}}_M}} \right)^{\rm{H}}}{\mathit{\boldsymbol{c}}_M} = 1, $

则拉格朗日参数λ表示如下:

$ \lambda = - \frac{1}{{\mathit{\boldsymbol{c}}_M^{\rm{H}}\mathit{\boldsymbol{R}}_x^{ - 1}{\mathit{\boldsymbol{c}}_M}}}. $

(4)

将式(4)代入式(3)，即可得到权矢量的最优解w_opt为

$ {\mathit{\boldsymbol{w}}_{{\rm{opt}}}} = \frac{{\mathit{\boldsymbol{R}}_x^{ - 1}{\mathit{\boldsymbol{c}}_M}}}{{\mathit{\boldsymbol{c}}_M^{\rm{H}}\mathit{\boldsymbol{R}}_x^{ - 1}{\mathit{\boldsymbol{c}}_M}}}. $

(5)

由于实际中难以得到理论协方差矩阵，一般利用快拍得到的样本协方差矩阵${\mathit{\boldsymbol{\tilde R}}_x}$进行替代，样本协方差矩阵的表达式如下

$ {{\mathit{\boldsymbol{\tilde R}}}_x} = \frac{1}{N}\sum\limits_{n = 1}^N {\mathit{\boldsymbol{x}}\left( n \right){\mathit{\boldsymbol{x}}^{\rm{H}}}\left( n \right)} . $

式中：N为采样得到的快拍数；x(n)为阵列天线实际接收到的包含目标信号、干扰信号与噪声的采样序列.

那么最终得到基于功率倒置算法的阵列权矢量${\mathit{\boldsymbol{\tilde w}}}$为

$ \mathit{\boldsymbol{\tilde w}} = \frac{{\mathit{\boldsymbol{\tilde R}}_x^{ - 1}{\delta _M}}}{{\delta _M^{\rm{H}}\mathit{\boldsymbol{\tilde R}}_x^{ - 1}{\delta _M}}}. $

(6)

2 Triangular-CMT-PI算法

基于信号模型和阵列天线模型，假设干扰信号来向发生变化，即

$ \tilde \theta = {\theta _0} + \Delta \theta . $

式中：θ₀为干扰信号初始来向；Δθ为干扰信号来向变化幅度；${\tilde \theta }$为干扰信号来向变化后的真实值，其规律是未知的.

假设第q个干扰信号来向变化幅度Δθ_q在区间[-β, β]内符合三角形分布，那么其概率分布函数如下

$ f\left( {\Delta {\theta _q}} \right) = \left\{ \begin{array}{l} \frac{{\Delta {\theta _q} + \beta }}{{{\beta ^2}}}, - \beta \le \Delta {\theta _q} \le 0;\\ \frac{{\beta - \Delta {\theta _q}}}{{{\beta ^2}}},0 < \Delta {\theta _q} \le \beta . \end{array} \right. $

阵列天线接收的信号包括卫星导航信号、干扰信号和噪声，由于卫星导航信号的功率很弱，一般比噪声还要低20 dB以上，在计算协方差矩阵时可以将其忽略，所以有

$ \mathit{\boldsymbol{\tilde R}}\left( x \right) \approx {{\mathit{\boldsymbol{\tilde R}}}_j} + {{\mathit{\boldsymbol{\tilde R}}}_{\rm{e}}} = \sum\limits_{q = 1}^Q {\sigma _q^2\mathit{\boldsymbol{a}}\left( {{\theta _q}} \right){\mathit{\boldsymbol{a}}^{\rm{H}}}\left( {{\theta _q}} \right)} + \sigma _{\rm{e}}^2\mathit{\boldsymbol{I}}. $

式中：σ_q²为第q个干扰信号的功率；a(θ_q)为第q个干扰信号的导向矢量；σ_e²为噪声信号的功率；I为单位矩阵.当干扰信号来向发生变化时，式(5)应该变化如下

$ {{\mathit{\boldsymbol{\tilde R}}}_L}\left( x \right) \approx \sum\limits_{q = 1}^Q {\sigma _q^2\int {f\left( {\Delta {\theta _q}} \right)\mathit{\boldsymbol{a}}\left( {{{\tilde \theta }_q}} \right){\mathit{\boldsymbol{a}}^{\rm{H}}}\left( {{{\tilde \theta }_q}} \right){\rm{d}}\Delta \theta } } + \sigma _{\rm{e}}^2\mathit{\boldsymbol{I}}, $

则可得${{\mathit{\boldsymbol{\tilde R}}}_L}\left( x \right)$第(m, n)个元素的表达式为

$ \begin{array}{l} {{\mathit{\boldsymbol{\tilde R}}}_L}\left( {m,n} \right) \approx \\ \sigma _{\rm{e}}^2{\delta _{mn}} + \sum\limits_{q = 1}^Q {\sigma _q^2\int_{ - \beta }^\beta {f\left( {\Delta {\theta _q}} \right)} } \cdot \\ \exp \left( { - j{\rm{ \mathsf{ π} }}\left( {m - n} \right)\sin {{\tilde \theta }_q}} \right){\rm{d}}\Delta {\theta _q} = \\ \sigma _{\rm{e}}^2{\delta _{mn}} + \sum\limits_{q = 1}^Q {\sigma _q^2\int_{ - \beta }^\beta {f\left( {\Delta {\theta _q}} \right)} } \cdot \\ \exp \left( { - j{\rm{ \mathsf{ π} }}\left( {m - n} \right)\sin \left( {{\theta _q} + \Delta {\theta _q}} \right)} \right){\rm{d}}\Delta {\theta _q} \approx \\ \sigma _{\rm{e}}^2{\delta _{mn}} + \sum\limits_{q = 1}^Q {\sigma _q^2\exp \left( { - j{\rm{ \mathsf{ π} }}\left( {m - n} \right)\sin {\theta _q}} \right)} \cdot \\ \int_{ - \beta }^\beta {f\left( {\Delta {\theta _q}} \right)\exp \left( { - j{\rm{ \mathsf{ π} }}\left( {m - n} \right)\Delta {\theta _q}\cos {\theta _q}} \right){\rm{d}}\Delta {\theta _q}} = \\ \sigma _{\rm{e}}^2{\delta _{mn}} + \sum\limits_{q = 1}^Q {\sigma _q^2\exp \left( { - j{\rm{ \mathsf{ π} }}\left( {m - n} \right)\sin {\theta _q}} \right)} \cdot \\ \frac{{2\left[ {1 - \cos \left( {{\rm{ \mathsf{ π} }}\left( {m - n} \right)\beta \cos {\theta _q}} \right)} \right]}}{{{{\left[ {{\rm{ \mathsf{ π} }}\left( {m - n} \right)\beta \cos {\theta _q}} \right]}^2}}}. \end{array} $

(7)

因此，扩张矩阵T_L的第(m, n)个元素可以表示为

$ {\mathit{\boldsymbol{T}}_L}\left( {m,n} \right) = \left\{ \begin{array}{l} \frac{{2\left[ {1 - \cos \left( {{\rm{ \mathsf{ π} }}\left( {m - n} \right)\beta \cos {\theta _q}} \right)} \right]}}{{{{\left[ {{\rm{ \mathsf{ π} }}\left( {m - n} \right)\beta \cos {\theta _q}} \right]}^2}}},m \ne n;\\ 1,m = n. \end{array} \right. $

可以看出当m≠n时，T_L(m, n)与每个具体干扰信号来向都有关系，在满足零陷宽度要求的前提下，为了避免求解干扰信号来向信息，假设每个干扰信号来向变化对应的零陷宽度均为最大值，即用β_max替代β_maxcos θ_q，这样新的扩张矩阵T_L的第(m, n)个元素可以表示为

$ {\mathit{\boldsymbol{T}}_L}\left( {m,n} \right) = \left\{ \begin{array}{l} \frac{{2\left[ {1 - \cos \left( {{\rm{ \mathsf{ π} }}\left( {m - n} \right){\beta _{\max }}} \right)} \right]}}{{{{\left[ {{\rm{ \mathsf{ π} }}\left( {m - n} \right){\beta _{\max }}} \right]}^2}}},m \ne n;\\ 1,m = n. \end{array} \right. $

由于新的扩张矩阵不需要包含干扰信号来向信息，式(7)可以进一步简化如下

$ \begin{array}{l} {{\mathit{\boldsymbol{\tilde R}}}_L}\left( {m,n} \right) \approx \sigma _{\rm{e}}^2{\delta _{mn}} + {\mathit{\boldsymbol{T}}_L}\left( {m,n} \right) \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{q = 1}^Q {\sigma _q^2\exp \left( { - {\rm{j \mathsf{ π} }}\left( {m - n} \right)\sin {\theta _q}} \right)} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{T}}_L}\left( {m,n} \right)\mathit{\boldsymbol{\tilde R}}\left( {m,n} \right). \end{array} $

定义“$ \circ $”表示Hardmard积，则采样协方差矩阵${{\mathit{\boldsymbol{\tilde R}}}}\left( x \right)$与锥化后的采样协方差矩阵${{\mathit{\boldsymbol{\tilde R}}}_L}\left( x \right)$满足下述关系：

$ {{\mathit{\boldsymbol{\tilde R}}}_L}\left( x \right) \approx \mathit{\boldsymbol{\tilde R}}\left( x \right) \circ {\mathit{\boldsymbol{T}}_L}. $

再根据式(6)，可得到协方差矩阵锥化后的阵列权矢量w_T为

$ {\mathit{\boldsymbol{w}}_{\rm{T}}} = \frac{{{{\left( {\mathit{\boldsymbol{\tilde R}}\left( x \right) \circ {\mathit{\boldsymbol{T}}_L}} \right)}^{ - 1}}{\mathit{\boldsymbol{c}}_M}}}{{\mathit{\boldsymbol{c}}_M^{\rm{H}}{{\left( {\mathit{\boldsymbol{\tilde R}}\left( x \right) \circ {\mathit{\boldsymbol{T}}_L}} \right)}^{ - 1}}{\mathit{\boldsymbol{c}}_M}}}. $

以上是基于半波长均匀线阵推导出的关系式，实际上该推导过程适用于任意形状的阵列结构.

3 仿真实验

实验采用七阵元半波长均匀线阵，BD2 B1频点信号来向为0°，信噪比为-20 dB；窄带干扰信号来向为-20°，干噪比为40 dB；BD2数字中频信号频率为2.42 MHz，信号采样率为8.184 MHz，快拍数为1 024.根据文献[15]中关于高动态的阐述，假设权值更新期间干扰来向变化幅度在3°以内；为消除随机误差影响，进行100次蒙特卡洛仿真实验.

仿真1  比较干扰信号来向变化统计模型分别服从三角形分布、均匀分布、拉普拉斯分布时协方差矩阵锥化零陷加宽算法的阵列方向增益特性.为了便于识别图形，将三角形分布、均匀分布、拉普拉斯分布对应的曲线分别标注为Tri、Uni、Lap.

图 2是PI算法、均匀分布算法、拉普拉斯分布算法和三角形分布算法对应的阵列天线方向增益图.首先，相比较于传统PI算法，3种分布都能有效展宽零陷.其次，就零陷深度在零陷宽度范围内的一致性而言，三角形分布和均匀分布算法对应的零陷深度在整个零陷范围内都保持较好，且三角形分布算法对应的零陷深度在整个零陷宽度上基本是优于均匀分布的；另外，虽然三角形分布算法对应的零陷极值没有拉普拉斯分布的零陷极值大，但是拉普拉斯分布的零陷深度也仅在不到1°幅度范围内保持优于三角形分布.

造成上述特点的原因可以从3种分布的概率密度曲线上进行解释，如图 3所示，均匀分布对应的概率密度曲线在整个变化区间上是等概率的，那么产生的零陷深度一致性保持要相对最好；三角形分布和拉普拉斯分布对应的概率密度曲线则在整个变化区间上呈现出中间高、两边低的特点，因此其零陷极值要大于均匀分布，拉普拉斯分布的概率密度曲线更为陡峭，其极值也更大，但是其在整个区间的零陷深度一致性保持上相对较差.由此可以看出，三角形分布在零陷深度以及宽度一致性之间要更加平衡，由于干扰来向变化规律实际未知，在拓展零陷范围内更加平稳的零陷深度将有利于干扰抑制，这是所提算法的优势.

图 4是PI算法、均匀分布算法、拉普拉斯分布算法和三角形分布算法对应的阵列天线输出信干噪比随输入信噪比变化的比较图.首先，3种分布算法对应的阵列输出信干噪比要明显高于PI算法，可知零陷展宽算法不仅展宽了零陷，也提高了信号的输出信干噪比；其次，三角形分布对应的阵列输出信干噪比又要好于均匀分布、拉普拉斯分布.

仿真2  结合BD2软件接收机，验证所提算法的可行性.首先利用PI算法、三角形分布算法、均匀分布算法和拉普拉斯分布算法分别对包含静态干扰的BD2 B1频点数据进行处理，分别得到权矢量；为考察展宽性能，以静态干扰信号来向为初始值，根据以上的仿真条件构造高动态干扰条件下的BD2数据，利用得到的初始权矢量分别处理该数据，最后送入BD2软件接收机对其进行捕获.

图 5(a)显示的是利用PI算法对处理后数据的捕获结果，可以看出PI算法并未有效抑制高动态干扰条件下的干扰信号，因此接收机无法有效捕获BD2信号.图 5(b)、(c)、(d)显示的是分别利用三角形分布、均匀分布和拉普拉斯分布算法得到的捕获结果，可以看出本文所提Triangular-CMT-PI算法能够有效抑制高动态干扰条件下的干扰信号，从而使接收机成功捕获到了BD2信号，验证了算法的可行性.需要注意的是均匀分布算法与拉普拉斯算法也是可行的，这在之前的文献中已经被证明过了.

4 结论

1) 本文提出的一种假设干扰来向变化服从三角形概率分布的协方差矩阵锥化零陷加宽抗干扰算法，能够在干扰来向上形成展宽零陷，在干扰来向快速变化的情况下保证其始终处于零陷之内.

2) 相比较于均匀分布算法和拉普拉斯分布算法，三角形分布算法拓宽后的干扰零陷更加平稳，阵列输出信干噪比要更高.正是由于干扰信号来向变化情况实际未知，平稳的零陷深度更加有利于抑制干扰信号，凸显了该算法的优点.

3) 通过与GNSS软件接收机结合，证明所提算法可以有效抑制干扰来向变化条件下的干扰信号，能够使软件接收机成功捕获卫星导航信号，具有一定的工程应用价值.