雷达、声呐和无线通信等应用对自适应波束形成的抗干扰能力和实时性提出了更高要求。传统基于最速迭代的自适应波束形成算法存在“过拟合”特性，导致在相干干扰条件下的干扰抑制性能急剧下降;另外，当干扰存在扰动且导向向量失配时，也无法有效抑制干扰。针对上述问题，本文提出一种基于共轭梯度(Conjugate Gradient, CG)加速的二次约束宽零陷干扰抑制自适应波束形成方法。该方法首先利用CG算法的快速收敛特性，完成采样协方差矩阵与导向向量间线性方程组的求解;其次将CG算法输出的权矢量作为迭代最速波束形成方法的初始权值，利用该方法的“过拟合”特性，确保对期望信号的强锁定;最后提出一种强化干扰特征的波达方向(Direction of Arrival, DOA)估计方法，实现宽带相干干扰下的干扰来波方向估计，并将该方法与二次约束零陷展宽方法结合，用于捕获干扰特征，形成自适应零陷。仿真实验验证了所提方法在单快拍、宽带相干干扰条件下，能够自适应抑制干扰且稳健性较好。
关键词：单快拍;宽带相干信号;自适应波束形成;零陷展宽
论文《基于共轭梯度加速的自适应二次约束零陷展宽方法研究》发表在《信号处理》，版权归《信号处理》所有。本文来自网络平台，仅供参考。
1 引言
自适应波束形成是阵列信号处理领域中的一项核心技术。自适应波束形成旨在确保接收信号在特定方向上实现最大增益，从而提高信号的强度和质量;此外，该技术还具备出色的干扰抑制能力，能够有效地减少或消除来自其他方向的干扰信号，进而提升系统的整体性能和可靠性。与传统的固定波束形成方法不同，自适应波束形成能够根据实时接收的信号特征和环境变化，自主调整波束的方向和形状，这种自适应性使其在复杂和多变的信号环境中，能更有效地应对干扰和噪声，提高信号质量和系统性能[1-3]。其中，最小均方(Least Mean Square, LMS)算法作为一种简单且计算负担较小的自适应算法，虽然收敛速度不如其他算法迅速，但在低计算资源需求的应用中具有显著优势[4]。
当前，随着设备和外部信号环境的复杂化，自适应波束形成的稳健性也越来越受到关注。针对实际需求，研究者提出了多种方法提高自适应波束形成的实时性和稳健性。变步长最小均方(Variable Step Size Least Mean Square, VSSLMS)算法通过设计不同的函数模型动态调整步长，相较于传统的LMS算法，可以更灵活地适应外部环境，快速求取最优权值[5-6]。然而，在信号剧烈变化的条件下，变步长算法可能会因为步长调整过大而导致系统不稳定;当误差接近零点时，系统可能会出现过度响应;且对初始步长较为敏感等。因此，需要设计有效的步长控制策略，使得系统在动态环境中能够平衡收敛速度和稳定性。
除了LMS类方法外，对角加载(Diagonal Loading, DL)也是一种广泛使用的自适应波束形成技术，可以处理协方差矩阵病态的问题[7-8]。近年来，研究者提出了自动计算加载量的算法，能够在低快拍和信号导向向量失配环境下，对干扰进行有效抑制，但计算量会显著增加[9]。基于不确定集的波束形成算法是一种通过建模期望信号方向向量失配误差问题，来增强波束形成鲁棒性的技术[10]，文献[11]提出一种算法，通过引入单快拍更新和迭代梯度方法，巧妙地减少了协方差矩阵运算的复杂度。此外，特征子空间算法在阵列信号处理中有着重要应用，文献[12]提出一种提高波束形成技术稳定性的方法，通过重构协方差矩阵，可以在协方差矩阵估计不准确的情况下，改善算法的性能。
上述各种方法均基于窄带信号模型，在低快拍、宽带强相干干扰条件下均存在局限性，无法有效抑制干扰。文献[13]提出的基于对数函数更新步长的LMS算法在宽带相干干扰条件下，由于信号与干扰的相似性，算法的步长更新以及权值更新均受到影响，导致算法无法有效抑制干扰。文献[14]引入一种自适应波束形成算法，利用最速下降法，将多变量二次凸优化问题简化为单变量二次问题，从而简化了问题的求解过程，并提高了计算效率。通过反复迭代来寻找极值，并逐步逼近高维凸优化问题的全局最优点，直至最终收敛到最优解。该方法在低快拍、宽带非相干干扰条件下，干扰抑制效果较好，但存在显著缺陷：算法的核心是采用一元二次迭代逼近最优权矢量，损失了数据的潜在联系，导致处理后的数据趋向于模型特征，即为了实现对多方向干扰的抑制，以损失干扰的信号特征为代价，导致无法对干扰实现自适应抑制，即使在非相干干扰条件下，其波束图也不能在干扰信号方向自适应地生成零陷，这种与期望信号过度拟合而导致的干扰特征缺失，可称为“过拟合”;其次，当期望信号与干扰信号相干时，算法的收敛速度急剧下降且无法对相干干扰进行有效抑制。
为解决这些问题，本文提出一种利用共轭梯度(Conjugate Gradient, CG)加速的二次约束宽零陷干扰抑制自适应波束形成算法。首先，利用CG算法的快速收敛特性，完成采样协方差矩阵与导向向量间线性方程组的求解;其次，将CG算法输出的权矢量作为迭代最速波束形成方法的初始权值，利用该方法的“过拟合”特性，确保对期望信号的强锁定;最后，提出一种强化干扰特征的波达方向(Direction of Arrival, DOA)估计方法，实现宽带相干干扰下的干扰来波方向估计，并将该方法与二次约束零陷展宽方法结合，用于捕获干扰特征，形成自适应零陷。通过上述过程，有效解决了迭代自适应波束形成算法的缺陷，并且在单快拍条件下能够有效抑制宽带相干干扰，性能显著提升。
2 信号模型与迭代最速下降算法
2.1 相干干扰环境下的信号模型
在阵列信号处理领域，相干干扰导致的信号抵消现象是一个常见且重要的问题。相干干扰指干扰信号与期望信号具有相同或相似的频率和相位关系，从而在接收端产生叠加效应，导致信号抵消或畸变，这种现象严重影响了信号接收的质量和系统性能。为深入理解和分析相干干扰导致的信号抵消现象，本文采用最佳波束形成的理论模型，该方法的原理可为分析相干干扰信号抵消现象提供理论支持[15-17]。
假设在一个M元阵列中，有一个宽带期望信号和n个宽带干扰信号分别从角度( heta_{0})与([ heta_{1}, heta_{2}, cdots, heta_{n}])入射。该阵列由各向同性阵元组成，即每个阵元在所有方向上具有相同的接收性能，期望信号与干扰信号之间是相干的。因此可得信源相干模型为：
[
X(t)=sum_{k=0}^{n} r_{k} aleft( heta_{k} ight) s_{0}(t)+sum_{k=n+1}^{K-1} aleft( heta_{k} ight) s_{k}(t)+N(t)
]
其中，期望信号的导向向量为(a( heta_{0}))，非相干干扰信号的导向向量分别为(a( heta_{n+1}), a( heta_{n+2}), cdots, a( heta_{K-1}))，而相干干扰信号的导向向量分别为(a( heta_{1}), a( heta_{2}), cdots, a( heta_{n}));相干干扰信号与期望信号之间的关系由复常数(r_{k}=s_{k}(t)/s_{0}(t))表示，其中(k=0,1,cdots,n)。
为简化信号处理模型并便于进一步的理论分析，通常假设期望信号、干扰信号(无论是相干还是非相干干扰)与噪声在统计上彼此独立。具体来说，期望信号的平均能量用(sigma_{s}^{2})表示，该值也代表了期望信号的功率;对于非相干干扰信号，其功率被定义为(sigma_{k}^{2})(其中(k=n+1, n+2, cdots, K-1))，用于描述干扰信号的平均能量特性。这种独立性假设不仅有助于简化计算和推导过程，还为复杂环境下的波束形成算法提供了基础支持。为使得噪声在不同阵元之间没有相关性，并且在时间上是平稳的，假设各阵元接收的噪声是相互独立的白噪声，噪声功率均为(sigma_{n}^{2})。在此设定下，阵列协方差矩阵可以定义为：
[
egin{gathered}
R_{X}=Eleft[X(t) X^{H} ight]= \
sigma_{s}^{2} a_{cb} a_{cb}^{H}+sum_{k=n+1}^{K-1} sigma_{k}^{2} aleft( heta_{k} ight) a^{H}left( heta_{k} ight)+sigma_{n}^{2} I= \
R_{s}+R_{i}+R_{n}
end{gathered}
]
其中，(a_{cb}=a( heta_{0})+sum_{k=1}^{n} r_{k} a( heta_{k}))，([·]^{H})表示共轭转置矩阵。对(R_{X})(协方差矩阵)进行特征值分解得到如下结果：
[
R_{X}=sum_{k=1}^{M} lambda_{k} u_{k} u_{k}^{H}
]
[
sum_{k=K-n+1}^{M} u_{k} u_{k}^{H}=I-sum_{k=1}^{K-n} u_{k} u_{k}^{H}
]
其中，特征值的排列顺序满足：(lambda_{1} geq lambda_{2} geq cdots geq lambda_{K-n} geq lambda_{K-n+1}=cdots=lambda_{M}=sigma_{n}^{2})。假设干扰信号强度远大于噪声，由式(5)可知：
[
egin{gathered}
R_{X}^{-1}=sum_{k=1}^{K-n} frac{1}{lambda_{k}} u_{k} u_{k}^{H}+frac{1}{sigma_{n}^{2}} sum_{k=K-n+1}^{M} u_{k} u_{k}^{H}= \
sum_{k=1}^{K-n} frac{1}{lambda_{k}} u_{k} u_{k}^{H}+frac{1}{sigma_{n}^{2}}left(I-sum_{k=1}^{K-n} u_{k} u_{k}^{H} ight)= \
frac{1}{sigma_{n}^{2}}left(I-sum_{k=1}^{K-n} frac{lambda_{k}-sigma_{n}^{2}}{lambda_{k}} u_{k} u_{k}^{H} ight) approx \
frac{1}{sigma_{n}^{2}}left(I-sum_{k=1}^{K-n} u_{k} u_{k}^{H} ight)=frac{1}{sigma_{n}^{2}} sum_{k=K-n+1}^{M} u_{k} u_{k}^{H}
end{gathered}
]
如果在相干干扰环境下选择应用最小方差无失真响应(Minimum Variance Distortless Response, MVDR)波束形成器，最优权矢量(w_{opt})在合成导向向量(a_{cb})方向的阵列响应则为：
[
F_{opt}left(a_{cb} ight)=w_{opt}^{H} a_{cb}=alpha^{*} a^{H}left( heta_{0} ight) R_{X}^{-1} a_{cb} approx 0
]
可知，最优权矢量(w_{opt})在合成方向矢量(a_{cb})方向上的阵列响应接近为0。
在非相干干扰方向上(w_{opt})的阵列响应表示为：
[
egin{gathered}
F_{opt}left( heta_{l} ight)=w_{opt}^{H} aleft( heta_{l} ight)=alpha^{*} a^{H}left( heta_{0} ight) R_{X}^{-1} aleft( heta_{l} ight) approx 0, \
l=n+1, n+2, cdots, K-1
end{gathered}
]
可以看出，在无关干扰方向上，应用该权重(w_{opt})后，阵列的响应几乎为0，能够有效抑制这些干扰，(w_{opt})在非相干干扰方向上生成零陷，进一步降低干扰的影响。因此，波束形成器的输出为：
[
egin{gathered}
y(t)=w_{opt}^{H} X(t)=w_{opt}^{H} a_{cb} s_{0}(t)+ \
sum_{k=n+1}^{K-1} w_{opt}^{H} aleft( heta_{k} ight) s_{k}(t)+w_{opt}^{H} N(t)= \
frac{alpha^{*}}{sigma_{n}^{2}} sum_{k=K-n+1}^{M} a^{H}left( heta_{0} ight) u_{k} u_{k}^{H} N(t)= \
frac{alpha^{*}}{sigma_{n}^{2}} sum_{k=K-n+1}^{M} ho_{k} u_{k}^{H} N(t)
end{gathered}
]
其中( ho_{k}=a^{H}( heta_{0}) u_{k})。在这种情况下，阵列输出的特性表现为只包含噪声的线性组合，不包含期望信号。
在相干干扰环境下，期望信号和干扰信号由于相干性而发生叠加，可能导致信号的相消。具体来说，这种相消现象表现为期望信号与干扰信号在某些阵列方向上的振幅相互抵消，导致接收到的期望信号能量大幅减少，甚至消失。
2.2 基于迭代最速方法的自适应波束形成
基于迭代最速下降的自适应波束形成方法的核心在于将多变量二次凸优化问题转换为单变量问题[14]，下面简要介绍该算法的推导过程。
假设多变量二次函数的一般形式是：
[
Fleft(x_{1}, x_{2}, cdots, x_{m} ight)=f_{1}^{2}+f_{2}^{2}+cdots+f_{n}^{2}
]
其中：
[
f_{1}=a_{11} x_{1}+a_{21} x_{2}+cdots+a_{m1} x_{m}+c_{1}
]
在多元函数优化过程中，最速下降算法的一个关键步骤是简化多元函数，使其成为单变量二次函数。这种简化过程通过选择一个主要变量，同时将其余变量固定为常数，实现对复杂多元函数的逐步优化。具体来说，将变量(x_{1})作为主要变量，同时将其余变量固定为常数值(x_{2}^{(0)}, x_{3}^{(0)}, cdots, x_{m}^{(0)})，此时，原函数简化为仅关于(x_{1})的二次函数，二次函数的极值点可以通过求解其一阶导数为零的条件来确定。
将(x_{2})作为自变量，所有其他变量被设定为初始值(x_{3}^{(0)}, x_{4}^{(0)}, cdots, x_{m}^{(0)})，求解一元二次函数一阶导数为零的条件确定极值点(x_{2}^{(1)})。
按照相同的步骤重复操作，确定每个变量的极值点(x_{1}^{(1)}, x_{2}^{(1)}, cdots, x_{m}^{(1)})。对上述过程进行循环迭代求取极值点，随着迭代次数的增加，极值点越来越接近多元函数的极值点。迭代公式为：
[
left{
egin{aligned}
x_{1}^{(k+1)}&=x_{1}^{(k)}-frac{left.sum_{l=1}^{n} f_{l} a_{1l} ight|_{x_{1}=x_{1}^{(k)}, x_{2}=x_{2}^{(k)}, cdots, x_{m}=x_{m}^{(k)}}}{sum_{l=1}^{n} a_{1l}^{2}} \
x_{2}^{(k+1)}&=x_{2}^{(k)}-frac{left.sum_{l=1}^{n} f_{l} a_{2l} ight|_{x_{1}=x_{1}^{(k+1)}, x_{2}=x_{2}^{(k)}, cdots, x_{m}=x_{m}^{(k)}}}{sum_{l=1}^{n} a_{2l}^{2}} \
&vdots \
x_{M}^{(k+1)}&=x_{M}^{(k)}-frac{left.sum_{l=1}^{n} f_{l} a_{Ml} ight|_{x_{1}=x_{1}^{(k+1)}, cdots, x_{M-1}=x_{M-1}^{(k+1)}, x_{M}=x_{M}^{(k)}}}{sum_{l=1}^{n} a_{Ml}^{2}}
end{aligned}
ight.
]
LMS自适应波束形成中阵列输出信号与期望信号的均方误差函数为：
[
J(w)=Eleft[left|w^{H}( heta) x(n)-d(n) ight|^{2} ight]=frac{1}{N} sum_{n=1}^{N}left[left|sum_{i=1}^{M} w_{i}^{*}( heta) x_{i}(n)-d(n) ight|^{2} ight]
]
均方误差函数展现了多元二次函数的典型结构，可以采用上述迭代方式求取极值点。进一步分析可得：
[
frac{partial}{partial w} J(w)=frac{2}{N} sum_{n=1}^{N} x(n)left[x^{H}(n) w-d^{*}(n) ight]
]
[
frac{partial^{2}}{partial w^{2}} J(w)=frac{2}{N} sum_{n=1}^{N} x(n) x^{H}(n)
]
令：
[
e(n)=x^{H}(n) w-d^{*}(n)
]
权系数迭代公式如下：
[
left{
egin{aligned}
w_{1}^{(k+1)}&=w_{1}^{(k)}-frac{sum_{n=1}^{N} e(n) x_{1}(n)ig|_{w_{1}=w_{1}^{(k)}, w_{2}=w_{2}^{(k)}, cdots, w_{M}=w_{M}^{(k)}}}{sum_{l=1}^{N} x_{1}(l) x_{1}^{*}(l)} \
w_{2}^{(k+1)}&=w_{2}^{(k)}-frac{sum_{n=1}^{N} e(n) x_{2}(n)ig|_{w_{1}=w_{1}^{(k+1)}, w_{2}=w_{2}^{(k)}, cdots, w_{M}=w_{M}^{(k)}}}{sum_{l=1}^{N} x_{2}(l) x_{2}^{*}(l)} \
&vdots \
w_{M}^{(k+1)}&=w_{M}^{(k)}-frac{sum_{n=1}^{N} e(n) x_{M}(n)ig|_{w_{1}=w_{1}^{(k+1)}, cdots, w_{M-1}=w_{M-1}^{(k+1)}, w_{M}=w_{M}^{(k)}}}{sum_{l=1}^{N} x_{M}(l) x_{M}^{*}(l)}
end{aligned}
ight.
]
为满足波束形成的性能要求，这里设定N个约束方程。接下来，为权系数w赋予随机初始值(w^{(0)})。在迭代过程中，通过权值迭代公式逐步计算每个权系数变量。根据约束方程的数量更新迭代方程中的系数，随后使用更新后的系数重新进行下一次迭代运算。重复上述过程，在处理完所有快拍数据后，记录最终计算得到的权系数w。
上述基于迭代最速下降的LMS自适应波束形成算法收敛速度较快，能够在低快拍条件下对强干扰形成有效抑制。然而，当出现相干干扰时，其性能急剧下降，算法无法在小快拍条件下捕捉到足够的信息特征，不能对干扰进行有效抑制。
2.3 共轭梯度(Conjugate Gradient, CG)算法
由于具有结构简单且收敛速度快的优点，共轭梯度(CG)算法在信号处理领域得到了广泛应用。CG算法最早被引入波束形成研究时，主要依托于广义旁瓣相消器(GSC)结构。然而，GSC的实现需要构造阻塞矩阵，这在一定程度上增加了计算复杂度，影响了其应用效率[18-20]。下面将对这种改进及其应用作简要概述。
在Capon波束形成器中，最优权值的计算是波束形成过程中至关重要的一步。为降低计算复杂度，常通过引入辅助矢量x，将原本需要进行的矩阵求逆操作简化为一个线性方程组的求解问题。
Capon波束形成器的最优权值表达式为：
[
w=frac{R^{-1} aleft( heta_{0} ight)}{a^{H}left( heta_{0} ight) R^{-1} aleft( heta_{0} ight)}
]
引入辅助矢量：
[
x=R^{-1} aleft( heta_{0} ight)
]
则：
[
w=frac{x}{a^{H}left( heta_{0} ight) x}
]
对式(20)两边同乘以R，可得：
[
R x=aleft( heta_{0} ight)
]
式(22)表示一个复数线性方程组，通过对这一方程组的求解，能够得到辅助矢量x的具体值。接着，将x的值代入式(21)中，可以计算出相应的波束形成器的权值矢量w。值得注意的是，式(22)对应的线性方程组求解过程也可以重新表述为一个最优化问题，这一转换不仅简化了计算过程，还使得问题的求解更加高效。通过这种优化方法，能够更快地获得所需的波束形成器权值，从而提升系统的整体性能。
[
F(x)=frac{1}{2} x^{T} R x-x^{T} a
]
设定集合({x_{1}, x_{2}, x_{3}, cdots})表示通过共轭梯度(CG)算法得到的最优化问题解的一系列迭代估计值。同时，定义集合({p_{1}, p_{2}, p_{3}, cdots})为约束问题式(23)所定义的二次表面上的一组共轭矢量，重要的是约束条件要求这些共轭矢量(p_{k})关于R正交，这一特性使得在求解过程中能够保持一定的几何特征，从而影响最终的解的性质。
[
p_{k}^{H} R p_{j}=0, quad k eq j
]
当选定(p_{k})作为搜索方向时，可以沿着这一方向进行搜索，从而得到迭代更新的形式(x_{k+1})。具体而言，这一过程涉及在(p_{k})的指导下，对当前估计值进行调整，以逐步逼近最优解。通过这样的迭代方式，能够有效地优化目标函数，并逐步改善解的精度。
[
x_{k+1}=x_{k}+alpha_{k} p_{k}, quad k=1,2,3, cdots, M
]
步长因子(alpha_{k})是本迭代方法中的重要参数，其值通过对代价函数(F(x))相对于(alpha_{k})的偏导数求解而得。
[
alpha_{k}=frac{r_{k}^{H} r_{k}}{p_{k}^{H} R p_{k}}
]
剩余矢量被定义为当前估计值与实际解之间的差距，通常用以衡量算法在每次迭代中产生的误差。这一矢量的计算不仅有助于评估当前解的质量，还能为后续迭代提供重要的反馈信息，进而指导优化过程的调整。定义剩余矢量为：
[
r_{k+1}=-left(R x_{k+1}-a ight)=-left(Rleft(x_{k}+alpha_{k} p_{k} ight)-a ight)=r_{k}-alpha_{k} R p_{k}
]
共轭矢量(p_{k})的迭代关系表示如下：
[
p_{k+1}=r_{k+1}+eta_{k} p_{k}
]
其中，(eta_{k})为步长系数，表示如下：
[
eta_{k}=frac{r_{k+1}^{H} r_{k+1}}{r_{k}^{H} r_{k}}
]
利用共轭梯度算法求解式(23)对应的最优化问题后，将得到的解代入式(19)，这一过程将直接计算出基于CG算法的Capon波束形成器的加权系数。通过这种方法，能够高效地获取所需的加权系数，从而优化波束形成器的性能。
具体计算步骤如下：
1. 计算协方差矩阵R，初始化：(x_{0}=0)，(r_{0}=a=p_{1})，(delta_{0}=left|r_{0} ight|^{2})，设定迭代终止条件(delta_{N})。
2. 对于(k=1,2,cdots)，有：
[
v_{k}=R p_{k}
]
[
alpha_{k}=frac{delta_{k-1}}{p_{k}^{H} v_{k}}
]
3. 更新(x_{k+1}=x_{k}+alpha_{k} p_{k})，计算剩余矢量(r_{k+1}=r_{k}-alpha_{k} v_{k})。
4. 计算(delta_{k}=left|r_{k+1} ight|^{2})。
5. 在算法执行过程中，判断(delta_{k} leq delta_{N})或者(k geq M)，如果满足任一条件，则立即终止运算，并跳转到步骤6进行后续处理;如果不满足这些条件，则计算(eta_{k})，并依据此更新共轭矢量(p_{k+1}=r_{k+1}+eta_{k} p_{k})，完成更新后，再次返回步骤2继续迭代。
6. 输出方程的解(x_{k+1})。
7. 计算加权系数：
[
w_{CG}=frac{x_{k+1}}{a^{H}left( heta_{0} ight) x_{k+1}}
]
采用共轭梯度(CG)算法求解Capon波束形成器的加权矢量，展现出良好的收敛速度和计算复杂性，同时保持了优秀的数值稳定性。该算法的显著优势在于其高数据一致性：在迭代的初始步骤中进行一次矩阵与向量的乘法运算，而在随后的步骤中仅需执行向量间的乘法操作，便可完成向量r、x和p的更新。这种高效的计算流程使得通过CG算法加速收敛并获取初始权值成为可能，为后续处理打下了坚实基础。
3 基于二次约束的零陷展宽方法
由上述论述可知，迭代自适应波束形成在宽带相干干扰条件下收敛速度急剧下降，且无法有效抑制干扰;同时，算法采用一元二次逼近最优解的方法，割裂了数据的潜在联系，使得算法本身即使在非相干干扰下依旧无法自适应生成零陷。因此，本文提出一种强化干扰特征的DOA估计方法，可以有效捕捉宽带相干干扰条件下的干扰来波方向;同时，为避免可能存在的DOA估计误差，将该方法与二次约束零陷展宽方法相结合，以抵消DOA估计误差对干扰抑制的影响，实现宽带相干干扰下的零陷展宽;将该方法与迭代自适应波束形成结合在一起，有效解决了迭代自适应方法的自适应缺陷问题，实现了干扰的自适应抑制。
下面对该算法的原理和过程进行论述：
1. 对数据矩阵进行解相干，这里使用常规的前后向空间平滑方法得到解相干的数据矩阵(R^{fb})。
2. 对前后空间平滑的数据协方差矩阵进行奇异值分解，提取信号特征：
[
R^{fb}=U Sigma V^{T}
]
去除U中信源列向量，得到(U_{n})，(U_{n})包含的信息仅与干扰、噪声相关。通过(U_{n})求取投影矩阵(G_{n})为：
[
G_{n}=U_{n} U_{n}^{H}
]
3. 将投影矩阵(G_{n})作为MUSIC算法的数据协方差矩阵，对其进行特征分解，即：
[
G_{n}=U_{s} Lambda_{s} U_{s}^{H}+U_{no} Lambda_{no} U_{no}^{H}
]
其中，(U_{s})是信号子空间的特征向量矩阵，(Lambda_{s})是信号特征值矩阵，(U_{no})是噪声子空间的特征向量矩阵，(Lambda_{no})是噪声特征值矩阵。
4. 由于噪声子空间与信号子空间正交，信号的导向向量(a( heta))与噪声子空间的特征向量正交，因此有：
[
U_{no}^{H} a( heta)=0
]
MUSIC算法利用噪声子空间构造谱函数为：
[
P_{MUSIC}( heta)=frac{1}{a^{H}( heta) U_{no} U_{no}^{H} a( heta)}
]
通过扫描所有可能的方向( heta)，当(a( heta))接近真实信号来波方向时，分母趋于零，谱函数(P_{MUSIC})出现峰值，其对应的角度即为来波方向(DOA)。
下面简述二次约束零陷展宽方法[15]：
设干扰信号方向为( heta_{i})，则理论上式(left|W_{0}^{H} a( heta_{i}) ight|=0)成立，式中(|cdot|)表示欧式范数，(W_{0})为前面波束形成算法求取的权矢量。这里希望在一定的角度范围( heta_{a} in[ heta_{i}-Delta heta / 2, heta_{i}+Delta heta / 2])内满足(left|W^{H} a( heta_{a}) ight|=0)，(Delta heta)为预设的零陷宽度。
考虑如下优化问题：
[
left{
egin{aligned}
min f(W)&=left|W-W_{0} ight|^{2} \
ext{s.t.} quad W^{H} Q W&=varepsilon
end{aligned}
ight.
]
[
Q=sum_{j=1}^{J} int_{ heta_{j}-frac{Delta heta_{j}}{2}}^{ heta_{j}+frac{Delta heta_{j}}{2}} a( heta) a^{H}( heta) d heta
]
其中，Q为干扰信号导向向量的积分约束矩阵，J表示干扰源的数量。根据具体的应用需求，可以调整不同干扰源的零陷宽度(Delta heta_{j})，以实现零陷的加宽效果，这一调整能够有效抑制特定方向上的干扰信号。为解决式(38)的优化问题，采用拉格朗日乘子法，通过引入拉格朗日乘子，将约束条件纳入优化过程，以便更好地获得所需的解。
[
egin{gathered}
f(W)=left(W-W_{0} ight)^{H}left(W-W_{0} ight)+gammaleft(W^{H} Q W-varepsilon ight)= \
W^{H} W-W^{H} W_{0}-W_{0}^{H} W+W_{0}^{H} W_{0}+gammaleft(W^{H} Q W-varepsilon ight)
end{gathered}
]
其中，(gamma)为拉格朗日乘子，(varepsilon)为零陷深度系数。
对(f(W))关于(W^{H})求导，并令结果为0，可得：
[
W_{opt}=(I+gamma Q)^{-1} W_{0}
]
其中，I是(M imes M)的单位阵。
从式(41)可以看出，阵列输出的权值与拉格朗日乘子(gamma)之间存在密切关系，因此，如何选择(gamma)将对波束形成的性能产生重要影响。输出权值需要满足约束条件(W^{H} Q W=varepsilon)。通过将式(41)代入这一约束条件，能够推导出求解(gamma)的具体公式，这一过程确保了在优化波束形成算法时，能够合理地调整拉格朗日乘子的值，以获得最佳的性能。
[
g(gamma)=W_{0}^{H}(I+gamma Q)^{-1} Q(I+gamma Q)^{-1} W_{0}=varepsilon
]
Q是对称矩阵，通过矩阵分解得到：
[
Q=U Lambda U^{H}=sum_{m=1}^{M} v_{m} u_{m} u_{m}^{H}
]
其中，(Lambda= ext{diag}(v_{1}, v_{2}, cdots, v_{M}))是对角阵，(U=(u_{1}, u_{2}, cdots, u_{M}))是Hermitian矩阵，(v_{m}(m=1,2,cdots,M))((v_{1} geq v_{2} geq cdots geq v_{M}))和(u_{m}(m=1,2,cdots,M))分别为Q的特征值和所对应的特征向量。
式(42)可写为：
[
egin{gathered}
g(gamma)=W_{0}^{H} U(I+gamma Lambda)^{-1} U^{H} U Lambda U^{H} U(I+gamma Lambda)^{-1} U^{H} W_{0}= \
sum_{m=1}^{M} frac{v_{m}left|W_{0}^{H} u_{m} ight|^{2}}{left(1+gamma v_{m} ight)^{2}}
end{gathered}
]
令(z_{m}=left|W_{0}^{H} u_{m} ight|)，则式(44)可写为：
[
g(gamma)=sum_{m=1}^{M} frac{v_{m}left|z_{m} ight|^{2}}{left(1+gamma v_{m} ight)^{2}}
]
综上所述，针对文献[14]提出算法在宽带相干干扰条件下干扰抑制能力急剧下降的问题，以及其自适应缺陷的问题，本文提出一种利用CG加速的二次约束宽零陷干扰抑制自适应波束形成算法。首先，利用CG算法的快速收敛特性，完成采样协方差矩阵与导向向量间线性方程组的求解;其次，将CG算法输出的权矢量作为迭代最速波束形成方法的初始权值，利用该方法的“过拟合”特性，确保对期望信号的强锁定;最后，提出一种强化干扰特征的DOA估计方法，实现宽带相干干扰下的干扰来波方向估计，并将该方法与二次约束零陷展宽方法结合，用于捕获干扰特征，形成自适应零陷。所提方法的复杂度为(O(3 M^{3}+(3+G+N z) M^{2}))，其中M为阵元个数，G为DOA估计中空域采样点数，(N z)为二次约束零陷展宽时的迭代求和次数。
本文所提方法的步骤见表1所示：
表1 算法步骤 Tab. 1 Algorithm steps
|步骤|详细内容|
|Step1|计算协方差矩阵R，初始化：(x_{0}=0)，(r_{0}=a=p_{1})，(delta_{0}=left|r_{0} ight|^{2})，迭代终止条件(delta_{N})。|
|Step2|根据式(30)、(31)、(25)更新(x_{k+1}=x_{k}+alpha_{k} p_{k})，计算剩余矢量(r_{k+1}=r_{k}-alpha_{k} v_{k})。|
|Step3|计算(delta_{k}=left|r_{k+1} ight|^{2})。|
|Step4|判断是否满足终止条件，即(delta_{k} leq delta_{N})或者(k geq M)，满足条件直接根据式(25)输出结果，计算权矢量(w_{CG});如果不满足条件，根据式(27)和(28)更新约束(p_{k+1})，并跳转Step2继续执行。|
|Step5|设定约束方程个数为N，使用式(18)进行权值迭代，一次迭代完成后根据约束方程数对迭代方程进行权系数更新，得到最佳权值(W_{0})。|
|Step6|根据式(33)和(34)得到解相干、强化干扰特征的投影矩阵(G_{n})。|
|Step7|将投影矩阵(G_{n})作为MUSIC算法的数据协方差矩阵，根据式(35)对其进行特征分解。|
|Step8|根据式(37)构造谱函数进行DOA估计。|
|Step9|依据DOA估计结果，根据式(39)求取干扰信号导向向量的积分约束矩阵Q。|
|Step10|依据式(42)求取拉格朗日乘子(gamma)。|
|Step11|根据式(41)求取最后的最优权矢量(W_{opt})。|
|Step12|根据(W_{opt})生成波束。|
4 仿真分析
本节对上述波束形成算法进行了实验分析和验证。仿真实验场景设定为：16个阵元组成均匀线性天线阵列;加性噪声为复高斯零均值的白噪声，并且每个阵元噪声方差相同;干扰设定为带宽50 MHz且与目标信号存在频谱重叠的线性调频宽带干扰;目标信号的中心频率设定为50 MHz，带宽设定为30 MHz，脉宽设定为10 µs，信号类型为线性调频脉冲信号;算法迭代步长为使算法收敛的最大允许步长。在该场景下，对经典LMS算法[1]、变步长LMS算法[6]、迭代LMS算法[14]、本文所提方法进行仿真分析。信噪比(Signal to Noise Ratio, SNR)、干扰噪声比(Interference to Noise Ratio, INR)等仿真参数设置见表2~表6所示。所有实验结果均为500次蒙特卡罗仿真后取平均的结果。
表2 系统及信号模型仿真参数 Tab. 2 Model simulation parameters
|参数|详细信息|
|阵元个数|16|
|信号带宽|30 MHz|
|信号脉宽|10 µs|
|信号方向|0°|
|信号中心频率|50 MHz|
|信号类型|线性调频信号|
|干扰带宽|50 MHz|
|干扰中心频率|40 MHz|
|干扰类型|压制性宽带干扰|
|噪声|加性复高斯零均值空时白噪声|
|快拍数|1|
表3 抗干扰性能仿真参数 Tab. 3 Anti-interference performance simulation parameters
|参数|详细信息|
|干扰方向|50°|
|信噪比(SNR)|0 dB|
|干扰噪声比(INR)|20 dB、50 dB、70 dB、100 dB、150 dB|
表4 自适应性能仿真参数 Tab. 4 Adaptive performance simulation parameters
|参数|详细信息|
|干扰方向|−60°、30°、50°|
|信噪比(SNR)|0 dB|
|干扰噪声比(INR)|50 dB|
表5 DOA估计性能与临近角度干扰抑制仿真参数 Tab. 5 Simulation parameters for DOA estimation performance and suppression of close-angle interference
|参数|详细信息|
|干扰方向|48°、50°|
|信噪比(SNR)|0 dB|
|干扰噪声比(INR)|50 dB|
表6 收敛速度仿真参数 Tab. 6 Convergence speed simulation parameters
|参数|详细信息|
|干扰方向|50°|
|信噪比(SNR)|0 dB|
|干扰噪声比(INR)|20 dB|
4.1 抗干扰能力对比
第一个实验分析算法的抗干扰能力。仿真实验中，将LMS算法的INR分别设置为20 dB、30 dB、40 dB、50 dB、70 dB;将变步长LMS算法的INR分别设置为20 dB、30 dB、40 dB、50 dB、60 dB;迭代LMS算法与本文所提方法的INR设置相同，分别为20 dB、50 dB、70 dB、100 dB、150 dB。其他参数设置参见表2和表3。各算法在不同INR下形成的波束如图1(a)~(d)所示。
由图1(a)、图1(b)可知，LMS算法、变步长LMS算法均无法对干扰进行有效抑制。图1(c)所示结果表明，在INR为20 dB时，迭代LMS算法能够在期望信号方向形成清晰的主瓣，同时旁瓣压制效果较好;然而，当INR上升到50 dB及以上时，迭代LMS算法的性能显著下降，主瓣位置发生偏移，旁瓣水平升高，无法正常工作。相比之下，如图1(d)所示，不论INR如何变化，本文所提方法始终能够保持其波束主瓣稳定地指向期望信号的方向，同时，在干扰方向上形成的零陷不仅稳定且深，显示出卓越的干扰抑制能力。综上所述，本文所提方法在抗干扰性能上具有明显的优势，抗干扰性能较强。
4.2 自适应能力对比
第二个实验分析算法的自适应能力，参数设置见表2和表4所示。仿真实验验证各个算法是否能根据干扰方位变化自适应生成对应零陷，结果见图2所示。
由图2(a)、图2(b)可以看出，LMS算法和变步长LMS算法在单快拍宽带相干干扰下无法收敛，在干扰方向不能形成零陷。由图2(c)可以看出，迭代LMS算法均无法根据干扰信号方向变化生成对应的零陷，自适应能力存在缺陷。由图2(d)可以看出，本文所提方法能够根据干扰方向变化分别在−60°、30°、50°方向生成零陷，在单快拍宽带相干干扰条件下工作正常。
图1 抗干扰性能分析 Fig. 1 Anti-jamming performance analysis
图2 自适应能力分析 Fig. 2 Adaptive capability analysis
4.3 输出信干噪比(Signal to Interference and Noise Ratio, SINR)随SNR变化对比
第三个实验分析各个算法输出SINR随SNR变化关系，此时干扰方向为50°，INR为20 dB，其他参数见表2。实验结果见图3所示。
可以看出，当SNR小于0 dB时，本文所提方法输出SINR较高，有效保证了低SNR下的期望信号的质量;同时，本文所提方法在低SNR下输出SINR较为稳定，能够有效抑制干扰与噪声。其他算法在低SNR下输出SINR较小且不稳定。比较而言，本文所提方法输出SINR较高且稳定，能够有效保证低SNR下的期望信号质量。
图3 输出SINR随SNR变化对比 Fig. 3 Comparison of output SINR and SNR
4.4 DOA估计性能与临近角度干扰抑制分析
第四个实验分析所提方法DOA估计性能和近角度多干扰的抑制效果，此时(INR=50 ext{dB})，(SNR=0 ext{dB})，其他参数见表2和表5。实验结果见图4所示。由图4(b)可以看出，改进后的MUSIC算法能够有效区分48°方向和50°方向的干扰，角度分辨率优于传统MUSIC方法。由图4(c)可以看出，本文所提方法能够对近角度多干扰形成宽零陷，零陷深度超过−70 dB，可以有效抑制干扰。因此，本文所提方法有效提升了传统MUSIC算法的角度估计分辨率，同时能够有效抑制临近角度的多个相干干扰。
图4 DOA估计性能与临近角度干扰抑制分析 Fig. 4 DOA estimation performance and suppression of close-angle interference
4.5 收敛速度对比
第五个实验分析阵列输出均方误差与快拍数之间的关系，详细参数见表2和表6。实验结果见图5所示，可以看出，本文所提方法在单快拍条件下输出均方误差已达到(10^{-2})以下，而变步长LMS算法、LMS算法的输出均方误差均在(10^{-1})以上。进一步比较可得，本文所提方法收敛速度快于LMS算法、迭代LMS算法，与变步长LMS算法的收敛速度相当。
图5 相关算法收敛速度对比图 Fig. 5 Comparison of convergence speed among algorithms
4.6 展宽零陷前后的主波束增益和宽度的性能变化
第六个实验分析展宽零陷前后的主波束增益和宽度的性能变化，采用3 dB凹口宽度与对应的零陷深度的乘积表示整体的干扰抑制效果。其中(SNR=0 ext{dB})，(INR=20 ext{dB})，干扰方向为50°，其他参数见表2。实验结果见图6所示。由图6(a)可知，展宽前3 dB凹口处的零陷深度为−44.2 dB，3 dB凹口宽度为0.3°，因此3 dB凹口宽度与对应的零陷深度的乘积为13.26。由图6(b)可知，展宽后3 dB凹口处的零陷深度为−94.8 dB，3 dB凹口宽度为0.2°，因此3 dB凹口宽度与对应的零陷深度的乘积为18.96。所以，展宽后的干扰抑制效果优于展宽前。由图6(c)、图6(d)可知，展宽前后主瓣的3 dB宽度无明显变化，且主副瓣比无明显变化，可见零陷展宽对主瓣增益并无影响。综上，零陷展宽后，干扰抑制性能提高，且对主瓣增益无明显影响。
图6 展宽零陷前后的主波束增益和宽度的性能变化对比图 Fig. 6 Comparison of main beam gain and width performance before and after null broadening
4.7 零陷展宽算法性能对比
第七个实验为不同零陷展宽算法性能对比分析，其中(SNR=0 ext{dB})，(INR=20 ext{dB})，干扰方向分别为60°、30°、50°，其他参数见表2。实验结果见图7所示。由图7(a)可以看出，现有二次约束零陷展宽算法形成的方向图性能较差，二次约束虽然形成零陷，但对主瓣和副瓣产生较大影响;同时，该算法在宽带相干干扰条件下无法有效捕捉干扰特征，只有在假定干扰方向已知的情况下才能形成图7(a)所示的效果，如果不假定干扰方向已知，则该算法在相干干扰条件下无法正常工作。由图7(b)可以看出，本文所提方法通过CG加速收敛，之后采用改进的DOA估计方法获取干扰方向，并使用二次约束算法强化干扰方向零陷，能够形成较好的干扰抑制效果，且主瓣和副瓣波形不受影响。
图7 零陷展宽算法对比图 Fig. 7 Comparison of zero-null broadening algorithms
5 结论
传统的最速迭代自适应波束形成算法在相干干扰的环境中，无法有效维持其通常表现出的干扰抑制能力，从而影响系统的性能。为了解决这一问题，本文提出一种利用CG加速的二次约束宽零陷干扰抑制自适应波束形成算法。首先，利用CG算法的快速收敛特性，完成采样协方差矩阵与导向向量间线性方程组的求解;其次，将CG算法输出的权矢量作为迭代最速波束形成方法的初始权值，利用该方法的“过拟合”特性，确保对期望信号的强锁定;最后，提出一种强化干扰特征的DOA估计方法，实现宽带相干干扰下的干扰来波方向估计，并将该方法与二次约束零陷展宽方法结合，用于捕获干扰特征，形成自适应零陷。通过上述方法，有效解决了迭代自适应波束形成算法的缺陷，并且在单快拍条件下能够有效抑制宽带相干干扰，性能较好。通过仿真实验可知，相比于已有算法，所提方法能够在单快拍、低SNR条件下对宽带相干干扰进行有效抑制且稳健性较好。
参考文献
[1] 王永良, 丁前军, 李荣锋. 自适应阵列处理[M]. 北京: 清华大学出版社, 2009. WANG Yongliang, DING Qianjun, LI Rongfeng. Adaptive Array Processing[M]. Beijing: Tsinghua University Press, 2009. (in Chinese)
[2] 孟昊宇, 杨小鹏, 高升, 等. 基于特征值斜投影的主瓣干扰抑制方法[J]. 信号处理, 2022, 38(2): 439-444. MENG Haoyu, YANG Xiaopeng, GAO Sheng, et al. Main lobe Interference suppression method based on eigenvalue oblique projection[J]. Journal of Signal Processing, 2022, 38(2): 439-444. (in Chinese)
[3] DAKULAGI V. Improved adaptive beamforming algorithms for wireless systems[J]. Wireless Personal Communications, 2023, 130(1): 625-633.
[4] 段彦亮. 稳健自适应波束形成算法研究[D]. 桂林: 桂林电子科技大学, 2024. DUAN Yanliang. Research on robust adaptive beamforming[D]. Guilin: Guilin University of Electronic Technology, 2024. (in Chinese)
[5] 陈东伟, 李明洁, 方庆园, 等. 一种改进变步长的自适应电磁干扰对消算法[J/OL]. 微波学报, 2023: 1-7 (2023-09-15)[2024-07-11]. https://kns.cnki.net/kcms/detail/32.1493.tn.20230913.2102.004.html. CHEN Dongwei, LI Mingjie, FANG Qingyuan, et al. An improved variable step size adaptive electromagnetic interference cancellation algorithm[J/OL]. Journal of Microwaves, 2023: 1-7 (2023-09-15)[2024-07-11]. (in Chinese)
[6] 肖玮, 刘思蔚, 唐宇龙, 等. 基于对称非线性函数的变步长LMS自适应滤波算法[J]. 火力与指挥控制, 2023, 48(2): 66-72, 85. XIAO Wei, LIU Siwei, TANG Yulong, et al. Variable step size LMS adaptive filtering algorithm based on the symmetric nonlinear function[J]. Fire Control & Command Control, 2023, 48(2): 66-72, 85. (in Chinese)
[7] JIAN Li, STOICA P, WANG Zhisong. On robust Capon beamforming and diagonal loading[J]. IEEE Transactions on Signal Processing, 2003, 51(7): 1702-1715.
[8] YANG Yue, XU Xu, YANG Huichao, et al. Robust adaptive beamforming via covariance matrix reconstruction with diagonal loading on interference sources covariance matrix[J]. Digital Signal Processing, 2023, 136: 103977.
[9] 殷冰洁. 基于多普勒特征恢复的声矢量阵鲁棒自适应波束形成方法[J]. 声学技术, 2018, 37(6): 596-600. YIN Bingjie. Doppler feature restoration based robust adaptive beamforming for acoustic vector-sensor array[J]. Technical Acoustics, 2018, 37(6): 596-600. (in Chinese)
[10] LI J, STOICA P, WANG Z. Doubly constrained robust Capon beamformer[C]//The Thirty-Seventh Asilomar Conference on Signals, Systems & Computers. Pacific Grove, CA, USA. IEEE, 2003: 1335-1339.
[11] 郭翔宇, 鄢社锋, 王文侠. 基于迭代梯度方法的线性约束稳健Capon波束形成快速算法[J]. 信号处理, 2021, 37(5): 712-723. GUO Xiangyu, YAN Shefeng, WANG Wenxia. A fast algorithm for linear constrained robust capon beamforming based on iterative gradient method[J]. Journal of Signal Processing, 2021, 37(5): 712-723. (in Chinese)
[12] 孔繁, 冯祺. 基于协方差矩阵重构的特征子空间投影稳健波束形成算法[J]. 中国电子科学研究院学报, 2016, 11(4): 371-375. KONG Fan, FENG Qi. Robust adaptive beamforming based on covariance matrix reconstruction and feature subspace projection[J]. Journal of China Academy of Electronics and Information Technology, 2016, 11(4): 371-375. (in Chinese)
[13] 王海旭. 基于对数函数的自适应波束形成算法研究[D]. 哈尔滨: 哈尔滨工程大学, 2022. WANG Haixu. Investigation on adaptive beamforming algorithm based on logarithmic function[D]. Harbin: Harbin Engineering University, 2022. (in Chinese)
[14] 王帅, 向建军, 彭芳, 等. 基于新最速下降算法的自适应波束形成[J]. 系统工程与电子技术, 2022, 44(7): 2104-2111. WANG Shuai, XIANG Jianjun, PENG Fang, et al. Adaptive beamforming based on new steepest descent algorithm[J]. Systems Engineering and Electronics, 2022, 44(7): 2104-2111. (in Chinese)
[15] 姚俊. 自适应波束形成的零陷加宽技术研究[D]. 南京: 南京信息工程大学, 2016. YAO Jun. Research on nulling widening technique based on adaptive beamforming[D]. Nanjing: Nanjing University of Information Science & Technology, 2016. (in Chinese)
[16] CHOI Y H. Adaptive beamforming for maximum SINR in the presence of correlated interferences[J]. Digital Signal Processing, 2024, 146: 104370.
[17] DU Yuxi, CUI Weijia, MEI Fengtong, et al. Robust adaptive beamforming algorithm for coherent signals based on virtual array[J]. Annals of Telecommunications, 2023, 78(9): 641-651.
[18] 汪倩. 雷达信号处理中自适应波束形成算法的研究与实现[D]. 西安: 西安电子科技大学, 2018. WANG Qian. The research and implementation of adaptive beamforming algorithm in radar signal processing[D]. Xi’an: Xidian University, 2018. (in Chinese)
[19] APOLINARIO J A, DE CAMPOS M L R, BERNAL O C P. The constrained conjugate gradient algorithm[J]. IEEE Signal Processing Letters, 2000, 7(12): 351-354.
[20] LANDAU L, DE LAMARE R C, WANG Lei, et al. Low-complexity robust beamforming based on conjugate gradient techniques[C]//2011 International ITG Workshop on Smart Antennas. Aachen, Germany. IEEE, 2011: 1-5.
