星载合成孔径雷达(synthetic aperture radar, SAR)具有全天时、全天候、高分辨对地观测等优越性能，为提高SAR图像的信息量，高分辨率和宽测绘带成像是星载SAR的重要发展趋势^[1-4]。在传统单通道SAR中，由于天线孔径尺寸和波束宽度为反比关系，为了实现宽测绘带所要求的宽发射波束，需要选取较小的天线俯仰向尺寸，导致发射增益降低使得雷达成像性能急剧恶化。同时，较宽的俯仰向波束宽度也会引入较明显的距离模糊干扰，降低成像质量。为此，高分辨率宽测绘带SAR成像实现的一个关键问题是如何在实现宽测绘带的同时保持较高信号功率和减少距离模糊。

数字波束形成(digital beam forming, DBF)技术^[5]是实现高分宽幅成像的一个重要手段，由数字处理技术和阵列天线理论结合而来，被应用于通信、声呐、雷达等多个领域^[5-7]。DBF技术可以将阵列天线的信息充分保留并利用复杂的数字处理方法在数字域对阵列信号进行处理从而实现所期望的功能。为解决宽测绘带实现所面临的问题，在距离向上，发射脉冲时使用较低增益宽发射波束，接收时将大尺寸接收天线分为多个子孔径接收，对各孔径接收的回波信号进行加权形成高接收增益窄波束，并使窄接收波束随着回波信号波达方向沿测绘带从近到远跟踪回波，这种DBF技术又称作扫描接收(scan-on-receive, SCORE)技术^[8]。该技术为DBF技术在DBF-SAR系统中所使用的一个主要技术，后续的很多改进也是在SCORE技术基础上进行研究以满足各种工作需求^[9-11]。

DBF处理中各个子孔径都独立产生数据，经过处理后合为一路数据。如果将DBF处理放在地面进行，则需要将各路数据都发送至地面，数据量会成倍增长，给星上存储和数据传输带来巨大压力。因此需要将DBF处理放在星上实时进行，而传统的SCORE算法和零陷算法都难以实时实现。为此文献[12]改进DBF处理结构, 并提出SCORE扫描实时处理方案；文献[13]提出基于线性约束最小方差波束形成算法(linearly constrained minimum variance, LCMV)的多零点DBF处理结构；文献[14]设计了基于LCMV的多零点DBF硬件处理方案，但此方案并没有解决LCMV方法中线速矩阵求逆对硬件资源的消耗问题，无法在卫星平台上实现。

为此，本文提出一套改进的工程上可实现的俯仰向多零点DBF实时处理方法，分析计算多零点DBF-SAR系统与传统SAR系统波位选择限制的区别；论证LDLT分解可用于波达矩阵求逆，并基于LDLT分解迭代进行线速矩阵求逆，解决了现有多点零陷DBF实时处理方案中的线速矩阵求逆问题。相较于已有方案，本方案改进了数据处理结构并大幅减少资源消耗，具备工程可实现性。

1 多零点DBF处理技术 1.1 SCORE扫描处理

对于俯仰向配置了N个通道的DBF-SAR系统，发射信号时使用宽波束覆盖整个测绘带，接收信号时使用高增益窄波束沿着时变的回波波达方向跟踪扫描接收，在快时间域上对整个测绘带进行扫描。DBF-SAR的几何关系如图 1所示。

图中O为星下点、P为目标点，R_e为地球半径，h为轨道高度，H_r=h+R_e为卫星到地心的距离，θ为下视角，β为天线安装角，阵列间隔为d，R_c为斜距。根据几何原理可以求得θ与R_c的关系为

$ \theta=\arccos \left(\frac{H_{\mathrm{r}}^{2}+R_{\mathrm{c}}^{2}-R_{\mathrm{e}}^{2}}{2 H_{\mathrm{r}} R_{\mathrm{c}}}\right) . $

(1)

在DBF-SAR系统中最常使用线性调频信号作为原始信号，DBF-SAR所发射的脉冲信号为

$ s_{\mathrm{T}}(t)=\operatorname{rect}\left(\frac{t}{T_{\mathrm{r}}}\right) \cdot \cos \left(2 {\rm{\mathsf{π}}} f_{\mathrm{c}} t+{\rm{\mathsf{π}}} K t^{2}\right)，$

(2)

其中：f_c为发射信号载频，K为发射信号调频率，T_r为脉冲宽度，t为快时间。发射信号传播至地面点目标发生后向散射形成回波信号被天线接收，第k个通道所接收到的射频回波信号为

$\begin{array}{c} s_{\mathrm{RF}, k}(t)=\operatorname{rect}\left(\frac{t-t_{k}}{T_{\mathrm{r}}}\right) \cdot \cos \left(2 {\rm{\mathsf{π}}} f_{\mathrm{c}}\left(t-t_{k}\right)+\right. \\ \left.{\rm{\mathsf{π}}} K\left(t-t_{k}\right)^{2}\right), \end{array} $

(3)

$ t_{k}=t_{0}+\Delta t_{k}=\frac{R_{0}+\Delta R_{k}}{c} , $

(4)

$ \Delta R_{k}=(k-1) \cdot d \cdot \sin \left(\theta_{\mathrm{t}}-\beta\right). $

(5)

其中：t_k表示回波被第k个通道接收的收发双程时延，t₀为参考通道的收发双程时延，R₀为参考通道的收发双程距离，ΔR_k为第k个通道与参考通道之间的收发双程距离差，θ_t为点目标实际的下视角。参考通道一般取第1个通道或者最中间的通道，本文取第1通道为参考通道，即最靠近星下点的接收通道。取中间通道为参考通道的情况可以由取第1通道为参考通道的情况加上一个相位差推导得到。射频回波信号经过下变频后得到中频信号为

$ \begin{gather*} s_{\mathrm{IF}, k}(t)=\operatorname{rect}\left(\frac{t-t_{k}}{T_{\mathrm{r}}}\right). \\ \cos \left(2 {\rm{\mathsf{π}}} f_{\mathrm{IF}} t+{\rm{\mathsf{π}}} K\left(t-t_{k}\right)^{2}-2 {\rm{\mathsf{π}}} f_{\mathrm{c}} t_{k}\right), \end{gather*} $

(6)

其中f_IF为中频信号的中心频率。中频信号经过数字IQ解调得到一对正交的I、Q基带信号，分别作为复信号的实部和虚部，可以得到基带信号的复信号形式为

$ s_{k}(t)=\operatorname{rect}\left(\frac{t-t_{k}}{T_{\mathrm{r}}}\right) \cdot \exp \left(-\mathrm{j} 2 {\rm{\mathsf{π}}} f_{\mathrm{c}} t_{k}+\mathrm{j} {\rm{\mathsf{π}}} K\left(t-t_{k}\right)^{2}\right), $

(7)

参考通道接收到的基带复信号s₀(t)为

$ s_{0}(t)=\operatorname{rect}\left(\frac{t-t_{0}}{T_{\mathrm{r}}}\right) \cdot \exp \left(-\mathrm{j} 2 {\rm{\mathsf{π}}} f_{\mathrm{c}} t_{0}+\mathrm{j} {\rm{\mathsf{π}}} K\left(t-t_{0}\right)^{2}\right). $

(8)

由此，第k个通道接收到的回波信号则可以表示为

$ s_{k}(t)=s_{0}\left(t+\Delta t_{k}\right) \cdot v_{k}(t), $

(9)

其中

$ v_{k}(t)=\exp \left(\mathrm{j} 2 {\rm{\mathsf{π}}}(k-1) \frac{d}{\lambda} \sin (\theta(t)-\beta)\right) . $

(10)

v_k(t)为该通道与参考通道之间接收回波的时间差所带来的相位调制，θ(t)为t时刻的波束指向角，λ为载波波长。在SCORE扫描中为实现波束扫描，需要使各通道回波信号相位对齐以进行相干叠加。所以SCORE扫描处理中所采用的DBF处理权重w_k(t)为波达相位v_k(t)的共轭，即

$ w_{k}(t)=v_{k}^{*}(t) \text {. } $

(11)

经过DBF处理后的回波信号为

$ \begin{gather*} s_{k}^{\prime}(t)=\operatorname{rect}\left(\frac{t-t_{k}}{T_{\mathrm{r}}}\right) \cdot \exp (-\mathrm{j} 2 {\rm{\mathsf{π}}} f_{\mathrm{c}} t_{0}+ \\ \mathrm{j} {\rm{\mathsf{π}}} K\left(t-t_{k}\right)^{2}) \cdot \exp (-\mathrm{j} 2 {\rm{\mathsf{π}}}(k-1) \frac{d}{\lambda} .\\ \left.\left[\sin (\theta(t)-\beta)-\sin \left(\theta_{\mathrm{t}}-\beta\right)\right]\right) . \end{gather*} $

(12)

文献[14]指出，为补偿脉冲延展损失(pulse extension loss, PEL)带来的影响，需要给各通道添加一个时延组件。T₀为场景中心回波时延，第k个通道的时延为

$ \Delta t_{k, \mathrm{PEL}}=-\left.\frac{(k-1) d \partial \theta(t)}{\lambda K}\right|_{t=T_{0}} . $

(13)

补偿后求和得到的SCORE扫描模式的输出为

$ \begin{gather*} s(t)=\sum\limits_{k=1}^{N}\left\{\operatorname{rect}\left(\frac{t-t_{k}-\Delta t_{k, \mathrm{PEL}}}{T_{\mathrm{r}}}\right) \cdot\right. \\ \exp \left(-\mathrm{j} 2 {\rm{\mathsf{π}}} f_{\mathrm{c}} t_{k}+\mathrm{j} {\rm{\mathsf{π}}} K\left(t-t_{k}-\Delta t_{k, \mathrm{PEL}}\right)^{2}\right). \\ \left.\exp \left(-\mathrm{j} 2 {\rm{\mathsf{π}}}(k-1) \frac{d}{\lambda} \sin \left(\theta\left(t-\Delta t_{k, \mathrm{PEL}}\right)-\beta\right)\right)\right\} . \end{gather*} $

(14)

1.2 多零陷点DBF处理

星载SAR实际工作中会收到干扰信号影响成像性能，比如距离模糊带来的干扰、星下点的干扰和外部射频信号干扰等。其中，以星下点干扰为代表，一般使用斑马图来分析这种干扰，如图 2所示。DBF系统往往要配合方位多通道系统实现高分宽幅，因此，脉冲重复频率(pulse repetition frequency, PRF)的选择相比方位单通道受到限制。以图 2(a)为例，在PRF为1 400~1 800 Hz的范围内，为了观测下视角25°~30°的区域，波位选择无法避开星下点，所以需要对星下点回波进行抑制。

如图 3，在接收斜距为$R_{\mathrm{c}}$的目标点$P$的回波信号时，同时会接收到来自斜距为$R_{m_{\mathrm{a}}}=R+m_{\mathrm{a}}$. $c /(2 \cdot \mathrm{PRF}), \left(m_{\mathrm{a}}= \pm 1, \pm 2, \cdots, \pm n\right)$的距离模糊点$P_{m_{\mathrm{a}}}$的回波信号。在传统单天线SAR中，距离模糊会对图像产生负面影响，所以要合理选取系统参数以避免距离模糊。但随着测绘带宽度的增加，距离模糊的影响会越来越大，所以需要对距离模糊干扰进行抑制。

由式(1)的星载SAR几何关系可以得知，由于斜距$R_{m_{\mathrm{a}}}$不同，所以距离模糊点$P_{m_{\mathrm{a}}}$的下视角$\theta_{m_{\mathrm{a}}}$与目标点$P$的下视角$\theta$不同。同时，星下点的下视角为$0^{\circ}$也与$\theta$不同，所以可以通过DBF对各通道信号加权，使接收波束在干扰信号下视角处形成零陷从而抑制这些干扰信号。

若一个$N$通道的DBF-SAR同时接收到了$M$个方向的回波信号，则该DBF-SAR系统接收到的原始回波信号$\boldsymbol{s}(t) \in \mathbb{C}^{N \times 1}$为

$ \boldsymbol{s}(t)=\sum\nolimits_{m=1}^{M} s_{0, m}(t) \cdot \boldsymbol{v}_{m}(t), $

(15)

其中

$ \begin{align*} \boldsymbol{v}_{m}(t)= & {\left[1, \exp \left(\mathrm{j} 2 {\rm{\mathsf{π}}} \frac{d}{\lambda} \sin \left(\theta_{m}(t)-\beta\right)\right), \cdots, \right.} \\ & \left.\exp \left(\mathrm{j} 2 {\rm{\mathsf{π}}}(N-1) \frac{d}{\lambda} \sin \left(\theta_{m}(t)-\beta\right)\right)\right]^{\mathrm{T}} \\ = & {\left[v_{1 m}(t), v_{2 m}(t), \cdots, v_{N m}(t)\right]^{\mathrm{T}} . } \end{align*} $

(16)

式中：$\theta_{m}$表示第$m$个波达方向的角度，$\boldsymbol{v}_{m} \in \mathbb{C}^{N \times 1}$为该波达方向所对应的天线导向矢量，$s_{0, m}(t)$为第$m$个点目标的参考回波，取

$ \boldsymbol{V}(t)=\left[\boldsymbol{v}_{1}(t), \boldsymbol{v}_{2}(t), \cdots, \boldsymbol{v}_{m}(t)\right], $

(17)

$ \boldsymbol{s}_{0}(t)=\left[s_{0, 1}(t), s_{0, 2}(t), \cdots, s_{0, m}(t)\right]^{\mathrm{T}} . $

(18)

$\boldsymbol{V}(t) \in \mathbb{C}^{N \times M}$为该DBF-SAR的波达矩阵，反映该系统接收到所有回波的空间信息。$\boldsymbol{s}_{0}(t) \in \mathbb{C}^{N \times 1}$为参考信号矢量，则回波信号可以写成

$ \boldsymbol{s}(t)=\boldsymbol{V}(t) \cdot \boldsymbol{s}_{0}(t) . $

(19)

DBF加权后的信号可以写成

$ s^{\prime}(t)=\boldsymbol{w}(t) \cdot \boldsymbol{s}(t)=\boldsymbol{w}(t) \cdot \boldsymbol{V}(t) \cdot \boldsymbol{s}_{0}(t), $

(20)

其中$\boldsymbol{w}(t) \in \mathbb{C}^{1 \times N}$为加权矢量。为使DBF加权后的信号在第$m$个波达方向形成最大增益，同时在其他波达方向形成零陷且不损失能量，即使得$s^{\prime}(t)=N \cdot \boldsymbol{e}_{m} \cdot \boldsymbol{s}_{0}(t)$，其中$\boldsymbol{e}_{m}$为单位矩阵的第$m$列，$\boldsymbol{w}_{m}(t)$应当满足

$ \boldsymbol{w}_{m}(t) \cdot \boldsymbol{V}(t)=N \cdot \boldsymbol{e}_{m} . $

(21)

根据LCMV方法，可以得到加权系数$\boldsymbol{w}_{m}(t)$的约束方程为

$ \left\{\begin{array}{l} \min \boldsymbol{w}_{m}(t) \cdot \boldsymbol{R} \cdot \boldsymbol{w}_{m}^{\mathrm{H}}(t) \\ \text { s. t. } \boldsymbol{w}_{m}(t) \cdot \boldsymbol{V}(t)=N \cdot \boldsymbol{e}_{m}, \\ M<N \end{array}\right. $

(22)

其中$\boldsymbol{R}=\boldsymbol{\sigma}_{0}^{2} \cdot \boldsymbol{I}$为接收端热噪声的相关阵^[15]。根据LCMV算法可以求得$\boldsymbol{w}_{m}(t)$的解为

$ \begin{gather*} \boldsymbol{w}_{m}(t)=N \cdot \boldsymbol{e}_{m}^{\mathrm{T}}\left(\boldsymbol{V}^{\mathrm{H}} \boldsymbol{V}\right)^{-1} \boldsymbol{V}^{\mathrm{H}}(t), \\ \boldsymbol{w}_{m}(t)=\left[w_{1, m}(t), w_{2, m}(t), \cdots, w_{N, m}(t)\right] . \end{gather*} $

(23)

其中$w_{n, m}(t)$为指向第$m$个波达角时第$n$个通道所用的加权系数。基于该加权矢量设计出来的方向图为

$ P(\theta)=P_{\mathrm{sub}}(\theta) \cdot G(\theta), $

(24)

其中

$ P_{\text {sub }}(\theta)=\left|\operatorname{sinc}\left(\frac{d \cdot \sin (\theta-\beta)}{\lambda}\right)\right|, $

(25)

$\begin{array}{c} G(\theta)=\sum\limits_{k=1}^{N} \exp \left(\mathrm{j} 2 {\rm{\mathsf{π}}}(k-1) \frac{d}{\lambda}(\sin (\theta-\beta))\right) .\\ w_{k, m}(t).\end{array} $

(26)

P_sub(θ)为天线的子阵图，G(θ)为阵列因子，该方向图的阵列因子具备周期性，所以当该方向图在θ₀处实现某一功能时，会存在θ_g处可以实现相似的功能。如θ₀处出现扫描的主瓣，则θ_g处也会出现对应的栅瓣；θ₀处实现零陷，θ_g处也会出现零陷的“栅瓣”，本文称其为“零陷栅瓣”。θ₀与θ_g的关系为

$ \begin{align*} & \exp \left(\mathrm{j} 2 {\rm{\mathsf{π}}}(k-1) \frac{d}{\lambda}\left(\sin \left(\theta_{0}-\beta\right)\right)\right) \cdot w_{k, m}(t)= \\ & \exp \left(\mathrm{j} 2 {\rm{\mathsf{π}}}(k-1) \frac{d}{\lambda}\left(\sin \left(\theta_{\mathrm{g}}-\beta\right)\right)\right) \cdot w_{k, m}(t) . \end{align*} $

(27)

仅考虑w_{k, m}(t)≠0, k≠1的情况，由周期性，该式可以化简为

$ 2 {\rm{\mathsf{π}}} \frac{d}{\lambda} \sin \left(\theta_{g}-\beta\right)=2 {\rm{\mathsf{π}}} \frac{d}{\lambda} \sin \left(\theta_{0}-\beta\right) \pm 2 n {\rm{\mathsf{π}}}, $

(28)

即

$ \theta_{\mathrm{g}}=\beta+\arcsin \left(\sin \left(\theta_{0}-\beta\right) \pm n \frac{\lambda}{d}\right), $

(29)

其中n为整数。据式(9)，当对星下点θ₀=0°进行零陷抑制时在别的角度也会出现“零陷栅瓣”，如果“零陷栅瓣”与主瓣重合则主瓣也会被抑制从而干扰成像，将由星下点抑制所产生的“零陷栅瓣”加入斑马图中可得图 2(b)。

在可以对星下点进行抑制的情况下，星下点回波可以不再是限制波位选择的一个因素，但仍然有一个特殊情况，即图 2(b)中星下点回波线与“零陷栅瓣”线相交的部分，这部分会受到星下点回波的干扰。但对星下点回波抑制则会受到“零陷栅瓣”的干扰，所以波位选择需要避开这种区域，如图 2(c)实心区域所示。

同理，对距离模糊进行零陷抑制也会产生“零陷栅瓣”，应有

$ \left\{\begin{array}{l}\theta_{\mathrm{g}}=\beta+\arcsin \left(\sin \left(\theta_{0}-\beta\right) \pm n \frac{\lambda}{d}\right) \\ R_{\mathrm{c}} \pm m_{\mathrm{a}} \cdot \frac{c}{2 \cdot \mathrm{PRF}} \neq H_{\mathrm{r}} \cos \theta_{\mathrm{g}}-\sqrt{R_{\mathrm{e}}^{2}-\left(H_{\mathrm{r}} \sin \theta_{\mathrm{g}}\right)^{2}}\end{array}\right., $

(30)

其中：m_a, n为整数，由此可以求得扫描角θ₀与PRF之间的限制关系，将其加入斑马图中得到考虑“零陷栅瓣”后用于波位选择的图 2(d)。竖直实线为第3节仿真实验所选择的波位。

2 多零点DBF实时处理方案

为实现星载DBF实时处理，DBF处理过程需要以较低的硬件资源消耗实时实现。分析前文加权矢量$\boldsymbol{w}_{m}(t)$的求解过程，可以将该过程分为2步来讨论，第1步求出波达矩阵$\boldsymbol{V}(t)$，第2步根据$\boldsymbol{V}(t)$求出$\boldsymbol{w}_{m}(t)$。为实现实时处理，这2步均需能进行线速处理。

计算波达矩阵即为求出各天线导向矢量，第m个方位角对应的天线导向矢量v_m(t)的第k项v_{k, m}(t)为

$ v_{k, m}(t)=\exp \left(\mathrm{j} \varphi_{k, m}(t)\right). $

(31)

φ_{k, m}(t)为第k通道接收到第m个目标点回波的波达相位

$ \begin{align*} \varphi_{k, m}(t) & =2 {\rm{\mathsf{π}}}(k-1) \frac{d}{\lambda} \sin \left(\theta_{m}(t)-\beta\right) \\ & =A F_{k} \cdot f_{m}(t) , \end{align*} $

(32)

其中

$ A F_{k}=(k-1) \cdot\left(2 {\rm{\mathsf{π}}} \frac{d}{\lambda}\right), $

(33)

$ f_{m}(t)=\sin \left(\theta_{m}(t)-\beta\right) . $

(34)

可以看出，$A F_{k}$是对应第$k$通道的常系数，所以对天线导向矢量$\boldsymbol{v}_{m}(\mathrm{t})$的实时计算便化为对系数$f_{m}(t)$的实时计算。

将SCORE扫描的波达角度设为$\theta_{1}$，干扰信号的波达角度则设为$\left\{\theta_{2}, \theta_{3}, \cdots, \theta_{M}\right\}$，根据星载几何关系，$\theta_{1}$所对应的斜距$R_{1}(t)=c\left(t+T_{0}\right) / 2$，其中$T_{0}$为场景中心所对应的双向时延，$f_{1}(t)$的表达式为

$ f_{1}(t)=\sin \left(\arccos \left(\frac{H_{\mathrm{r}}^{2}+\left(\frac{c\left(t+T_{0}\right)}{2}\right)^{2}-R_{\mathrm{e}}^{2}}{2 H_{\mathrm{r}}\left(\frac{c\left(t+T_{0}\right)}{2}\right)}\right)-\beta\right). $

(35)

此时f₁(t)尽管本身很复杂, 但函数图像为接近直线的单调递增的弧形，且t的取值范围也很小一般不会超过ms级, 所以泰勒展开近似具备很好的拟合效果。根据泰勒展开式，展开到3阶近似得

$ g_{1}(t)=A_{0}+A_{1} t+A_{2} t^{2}+A_{3} t^{3} \approx f(t). $

(36)

其中{A₀, A₁, A₂, A₃}分别为f(t)在t=0处展开的0至3阶泰勒级数系数。文献[12]对该系数进行了详细计算，这些系数极为复杂，但是在一次DBF-SAR工作过程中这些系数不会变动，不需要在星上实时计算，可以由辅助计算机求得后输入硬件处理模块即可。

基于上述理论，SCORE波束扫描权重生成所用到的v_{k, 1}(t)实时生成的结构如图 4所示。

为生成波达矩阵，需要生成其余波达方向所对应的波达相位，观察星下点和距离模糊点与目标点之间的关系可以发现，星下点的波达方向固定为θ_sub(t)=0, 即星下点所对应的f_sub(t)=sin(-β)为常数，不需要额外计算。而m_a阶距离模糊点与目标点之间相差固定时间，$\Delta t_{m}= \pm m_{\mathrm{a}}$ • $\mathrm{PRI} / 2, \mathrm{PRI}=1 / \mathrm{PRF}$，因此$f_{\mathrm{amb}, m_{\mathrm{a}}}(t)=f_{1}\left(t+\Delta t_{m}\right)$，对于低阶距离模糊点，可以求得距离模糊点所用的泰勒近似$g_{\mathrm{amb}, m_{\mathrm{a}}}(t)$为

$ \begin{align*} & g_{\mathrm{amb}, m_{\mathrm{a}}}(t)=g_{1}\left(t+\Delta t_{m}\right) \\ = & A_{0}+A_{1}\left(t+\Delta t_{m}\right)+A_{2}\left(t+\Delta t_{m}\right)^{2}+A_{3}\left(t+\Delta t_{m}\right)^{3} \\ = & A_{0}^{\prime}+A_{1}^{\prime} t+A_{2}^{\prime} t^{2}+A_{3}^{\prime} t^{3} \\ \approx & f\left(t+\Delta t_{m}\right)=f_{\mathrm{amb}, m_{\mathrm{a}}}(t) . \end{align*} $

(37)

其中系数$\left\{A_{0}^{\prime}, A_{1}^{\prime}, A_{2}^{\prime}, A_{3}^{\prime}\right\}$由$\left\{A_{0}, A_{1}, A_{2}, A_{3}\right\}$与时间差$\Delta t_{m}$简单运算得到

$ \left\{\begin{array}{l} A_{0}^{\prime}=A_{0}+A_{1} \Delta t_{m}+A_{2} \Delta t_{m}^{2}+A_{3} \Delta t_{m}^{3} \\ A_{1}^{\prime}=A_{1}+A_{2} \Delta t_{m}+A_{3} \Delta t_{m}^{2} \\ A_{2}^{\prime}=A_{2}+A_{3} \Delta t_{m} \\ A_{3}^{\prime}=A_{3} \end{array}\right. . $

(38)

据此，可以通过$\left\{A_{0}, A_{1}, A_{2}, A_{3}\right\}$求得由扫描角、星下点、距离模糊构成的波达矩阵$\boldsymbol{V}(t)$。为实现多点零陷，需要根据$\boldsymbol{V}(t)$求得加权矢量$\boldsymbol{w}(t)$使得$\boldsymbol{w}(t) \cdot \boldsymbol{s}(t)=N \cdot \boldsymbol{e}_{1} \cdot \boldsymbol{s}_{0}(t)$。为便于计算，文献[13]将多点零陷DBF处理分为2步，将加权矢量$\boldsymbol{w}(t)$变为

$ \boldsymbol{w}(t)=N \cdot \boldsymbol{e}_{1}\left(\boldsymbol{V}^{\mathrm{H}} \boldsymbol{V}\right)^{-1} \boldsymbol{V}^{\mathrm{H}}(t) \\ =N \cdot \boldsymbol{e}_{1} \cdot \boldsymbol{Y}(t) \cdot \boldsymbol{X}(t), $

(39)

$ \boldsymbol{X}(t)=\boldsymbol{V}^{\mathrm{H}}(t)=\left(\begin{array}{ccc} x_{11}(t) & \cdots & x_{1 N}(t) \\ \vdots & & \vdots \\ x_{M 1}(t) & \cdots & x_{M N}(t) \end{array}\right), $

(40)

$ \boldsymbol{Y}(t)=\left(\boldsymbol{V}^{\mathrm{H}} \boldsymbol{V}\right)^{-1}(t)=\left(\begin{array}{ccc} y_{11}(t) & \cdots & y_{1 M}(t) \\ \vdots & & \vdots \\ y_{M 1}(t) & \cdots & y_{M M}(t) \end{array}\right), $

(41)

所以有

$ \begin{align*} \boldsymbol{w}(t) & =N \cdot \boldsymbol{e}_{1} \cdot \boldsymbol{Y}(t) \cdot \boldsymbol{X}(t) \\ & =N \cdot\left(y_{11}, y_{12}, \cdots, y_{1 M}\right) \boldsymbol{X}(t) . \end{align*} $

(42)

因$\boldsymbol{X}(t)$为波达矩阵的共轭转置，在硬件处理中可以直接通过对波达矩阵重新排线和虚部转为负数完成，核心问题就是求$\left(y_{11}, y_{12}, \cdots, y_{1 M}\right)$可硬件实现的简单表达式，即求矩阵$\left(\boldsymbol{V}^{\mathrm{H}} \boldsymbol{V}\right)^{-1}$的第1行，本文采取LDLT分解法来实现线速矩阵求逆，LDLT分解法是一种用于对称复正定矩阵的分解方法，计算量较小且数值稳定，该方法将矩阵分解为单位下三角阵$\boldsymbol{L}$与对角阵$\boldsymbol{D}$的乘积，而这2种矩阵的求逆较为简单，从而将对复杂矩阵的求逆变为对多个简单矩阵的求逆与相乘。

用矩阵$\boldsymbol{Z}$表示$\boldsymbol{Z}=\boldsymbol{V}^{\mathrm{H}} \boldsymbol{V}, ~$因为$\boldsymbol{Z}=\boldsymbol{Z}^{\mathrm{H}}, ~$所以$\boldsymbol{Z}$为埃尔米特矩阵，且由于$\boldsymbol{Z}$对角线元素为实数，所以矩阵$\boldsymbol{Z}$的特征值均为实数，取任意非零复向量$\boldsymbol{x} \in \mathbb{C}^{M \times 1}$，考察$\boldsymbol{x}^{\mathrm{H}} \boldsymbol{Z} \boldsymbol{x}$可得

$ \boldsymbol{x}^{\mathrm{H}} \boldsymbol{Z} \boldsymbol{x}=\boldsymbol{x}^{\mathrm{H}} \boldsymbol{V}^{\mathrm{H}} \boldsymbol{V} \boldsymbol{x}=(\boldsymbol{V} \boldsymbol{x})^{\mathrm{H}} \boldsymbol{V} \boldsymbol{x} \geqslant 0 . $

(43)

由此可以证明Z为半正定矩阵，设

$ \boldsymbol{V}=\left[\begin{array}{c} \boldsymbol{V}_{M} \\ \boldsymbol{V}_{N-M} \end{array}\right], $

(44)

其中：$\boldsymbol{V}_{M} \in \mathbb{C}^{M \times M}, \boldsymbol{V}_{N-M} \in \mathbb{C}^{N-M \times M}$，分别为矩阵$\boldsymbol{V}$的前$M$行和后$N-M$行，则

$ \boldsymbol{V} \boldsymbol{x}=\left[\begin{array}{c} \boldsymbol{V}_{M} \boldsymbol{x} \\ \boldsymbol{V}_{N-M} \boldsymbol{x} \end{array}\right] . $

(45)

且由于矩阵$\boldsymbol{V}_{M}$为范德蒙矩阵且列不相关，所以$r\left(\boldsymbol{V}_{M}\right)=\boldsymbol{M}, \boldsymbol{V}_{M}$为满秩矩阵，对于$M$维向量$\boldsymbol{x}$，一定使得$\boldsymbol{V}_{M} \boldsymbol{x} \neq \boldsymbol{0}$，所以也可得$\boldsymbol{V} \boldsymbol{x} \neq \boldsymbol{0}, \boldsymbol{V} \boldsymbol{x}$不为零向量，因此$(\boldsymbol{V} \boldsymbol{x})^{\mathrm{H}} \boldsymbol{V} \boldsymbol{x} \neq \boldsymbol{0}$，所以$\boldsymbol{x}^{\mathrm{H}} \boldsymbol{Z} \boldsymbol{x}>\boldsymbol{0}$，所以矩阵$\boldsymbol{Z}$为复正定矩阵，$\boldsymbol{Z}$可唯一分解为$\boldsymbol{Z}= \boldsymbol{L} \boldsymbol{D} \boldsymbol{L}^{\mathrm{H}}$，其中

$ \boldsymbol{L}=\left(\begin{array}{cccc} 1 & 0 & \cdots & 0 \\ l_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{M 1} & l_{M 2} & \cdots & 1 \end{array}\right), $

(46)

$ \boldsymbol{D}=\operatorname{diag}\left(d_{1}, d_{2}, \cdots, d_{M}\right) . $

(47)

$\boldsymbol{L} 、\boldsymbol{D}$矩阵中的元素和$\boldsymbol{Z}$矩阵中第$i$行第$j$列元素$z_{i j}$有以下关系

$ d_{j}=z_{j j}-\sum\nolimits_{m=1}^{j-1} l_{j m} l_{j m}^{*} d_{m} , $

(48)

$ l_{i j}=\frac{z_{i j}-\sum\nolimits_{m=1}^{j-1} l_{i m} l_{j m}^{*} d_{m}}{d_{j}}. $

(49)

其中$1 \leqslant j \leqslant M, j+1 \leqslant i \leqslant M$。除$d_{j}$涉及倒数外，其余运算均为复数乘法和加法，为减少乘法次数，引入辅助量$u_{i j}=l_{i j} d_{j}$，将上式化简为

$ d_{j}=z_{j j}-\sum\nolimits_{m=1}^{j-1} u_{j m} l_{j m}^{*}, $

(50)

$ u_{i j}=z_{i j}-\sum\nolimits_{m=1}^{j-1} u_{i m} l_{j m}^{*}, $

(51)

$ l_{i j}=\frac{u_{i j}}{d_{j}}. $

(52)

由此便得到了矩阵$\boldsymbol{L}$和$\boldsymbol{D}$使得$\left(y_{11}, y_{12}, \cdots\right.$，$\left.y_{1 M}\right)=\boldsymbol{e}_{1}{ }^{\mathrm{H}}\left(\boldsymbol{L} \boldsymbol{D} \boldsymbol{L}^{\mathrm{H}}\right)^{-1}$，设$\boldsymbol{L}$的逆矩阵为$\boldsymbol{K}=\boldsymbol{L}^{-1}$，便有

$ \left(y_{11}, y_{12}, \cdots, y_{1 M}\right)=\boldsymbol{e}_{1}{ }^{\mathrm{H}} \boldsymbol{K}^{\mathrm{H}} \boldsymbol{D}^{-1} \boldsymbol{K} . $

(53)

其中，$\boldsymbol{K}$也为下三角矩阵，$\boldsymbol{K}$的第$i$行第$j$列元素$k_{i j}, ~(1 \leqslant j \leqslant i \leqslant M)$为

$ k_{i j}=-\sum\nolimits_{m=j}^{i-1} l_{i m} k_{m j} . $

(54)

D的逆矩阵为

$ \boldsymbol{D}^{-1}=\operatorname{diag}\left(1 / d_{1}, 1 / d_{2}, \cdots, 1 / d_{M}\right), $

(55)

由此可以求得

$ \left\{\begin{array}{l} y_{11}=\frac{1}{d_{1}}+\frac{k_{21}^{*} k_{21}}{d_{2}}+\cdots+\frac{k_{M 1}^{*} k_{M 1}}{d_{M}}, \\ y_{12}=\frac{k_{21}^{*}}{d_{2}}+\frac{k_{31}^{*} k_{32}}{d_{3}}+\cdots+\frac{k_{M 1}^{*} k_{M 2}}{d_{M}}, \\ \quad \vdots \\ y_{1 m}=\frac{k_{m 1}^{*}}{d_{m}}+\frac{k_{m+1, 1}^{*} k_{m+1, 2}}{d_{m+1}}+\cdots+\frac{k_{M 1}^{*} k_{M m}}{d_{M}}, \\ \quad \vdots \\ y_{1 M}=\frac{k_{M 1}^{*}}{d_{M}} . \end{array}\right. $

(56)

由此即可求得对应的$\left(y_{11}, y_{12}, \cdots, y_{1 M}\right)$。上述求逆过程是一个迭代过程，其中LDLT分解的复杂度为$O\left(\frac{1}{3} N^{3}\right)$，求逆矩阵$\boldsymbol{K}$的复杂度为$O\left(\frac{1}{3} N^{3}\right)$，求$y$的复杂度为$O\left(\frac{3}{2} N^{2}\right)$，整个求逆过程的综合复杂度为$O\left(\frac{2}{3} N^{3}\right)$。具体迭代过程如图 5所示，迭代过程从左至右，标注数字为在基于现场可编程逻辑门阵列(field programmable gate array，FPGA)进行浮点数运算的情况下所需要的FPGA时钟周期数，完成一次迭代共需要108个FPGA时钟周期。

图中$\boldsymbol{D}_{n} 、\boldsymbol{L}_{n} 、\boldsymbol{U}_{n} 、\boldsymbol{K}_{n} 、\boldsymbol{Y}_{n}^{-1}$分别为第$n$次迭代时对应矩阵或向量有进行迭代计算的元素，以$\boldsymbol{D}_{n}$为例，$\boldsymbol{D}_{n}$为

$ \boldsymbol{D}_{n}=\left\{d_{n}, d_{n, n+1}^{\prime}, \cdots, d_{n, M}^{\prime}\right\}, $

(57)

其中$d_{n, j}^{\prime}$为$d_{j}$在第$n$次迭代时的数值

$ d_{n, j}^{\prime}=a_{j j}-\sum\nolimits_{m=1}^{n-1} u_{j m} l_{j m}^{*}=d_{n-1, j}^{\prime}-u_{j, n-1} l_{j, n-1}^{*}, $

(58)

其中n+1≤j≤M，所以有

$ d_{j}=d_{n, j}^{\prime}-\sum\limits_{m=n}^{j-1} u_{j m} l_{j m}^{*} . $

(59)

因此，第$n$次迭代将求出$\boldsymbol{D}$矩阵的第$n$个元素，$\boldsymbol{U}$矩阵和$\boldsymbol{L}$矩阵的第$n$列，$\boldsymbol{K}$矩阵的第$n$行；$\left(y_{11}, y_{12}, \cdots, y_{1 M}\right)$的前$n$个元素开始迭代，在第$M$次迭代时$\left(y_{11}, y_{12}, \cdots, y_{1 M}\right)$会同时算出并输出。此外，$D_{n}^{-1}$为求矩阵$\boldsymbol{D}$第$n$项的倒数即

$ D_{n}^{-1}=1 / d_{n} \text {. } $

(60)

D_n^-1求倒数过程为本方法中占用资源较高的部分，但其运算次数与方向角个数M相等，不会占用大量硬件资源。

由于LDLT分解求逆是一个迭代过程，所以无法即时实现对波达矩阵的求逆。但该迭代过程仅包含加法、乘法和少量取倒数这种FPGA可以线速实现的运算，且迭代过程中仅向前迭代。整个迭代过程中每个元素只会用于该元素被求出以后的那一次迭代，而不需要在多次迭代中重复使用。所以该迭代过程可以在不额外消耗大量硬件资源的情况下设计为线速迭代过程，实现波达矩阵以线速输入，经过迭代后以线速输出(y₁₁, y₁₂, …, y_1M)，达成吞吐量平衡，满足多零陷DBF实时处理工作需求。同时，LDLT分解迭代过程是完全解析的，没有使用近似引入额外误差。本方案的工程误差来源为产生波达矩阵时的近似处理和FPGA处理本身的硬件环境限制。产生波达矩阵时的近似处理所带来的误差与文献[12]的误差一致，处于可以接受的范围。

3 仿真实验

为验证本文实时处理方法的可行性和成像效果，现设计仿真实验，仿真参数如表 1所示，波位根据图 2(d)选择，且由图 2(a)可知，在PRF为1 600 Hz的情况下，24.1°~28.1°的下视角范围会受到星下点回波干扰。设置多点目标如表 2所示。

本文硬件部分采用FPGA实现，使用vivado作为仿真软件，基于xilinx virtex-7 xc7vx690t进行仿真开发，将点目标仿真数据输入到vivado进行FPGA仿真DBF处理，然后对处理后的数据进行SAR成像以判断可行性以及分析成像效果。DBF-SAR的FPGA处理整体结构如图 6所示。

FPGA处理时序如图 7所示，输入信号为16通道8位模数转换信号，输出信号为单通道32位定点IQ信号，从输入原始回波到输出处理后的信号需要2 910 ns，FPGA时钟主频为200 MHz，即输入输出间隔582个时钟。在FPGA处理中降低权重更新率使其与时钟主频对齐，即6个时钟的原始回波对应1个权重，以雷达灵敏度降低不超过0.1 dB为代价使得FPGA硬件资源消耗将降为约1/6，且权重更新率可以进一步降低以降低FPGA中的时序限制，根据文献[14]中的分析，权重更新率和时钟主频的比值Φ小于400都可以使灵敏度降低不超过0.1 dB。

4波达角零陷DBF实时处理时FPGA资源消耗如表 3所示，本文采用的LDLT迭代求逆方案所需资源为第1行，I/O资源消耗受原始回波信号影响，即取决于通道数和模数转换采样数，不受后续处理方案的影响，其余种类的资源中最大资源消耗占比为DSP资源占比15.22 %，主要用于矩阵相乘时的快速计算，若DSP资源有压力可以将部分DSP资源由逻辑组合代替。根据资源消耗量可以认为该处理方案足以满足星载平台硬件资源限制，且具备可拓展性。在FPGA工程设计中，资源使用量一般不超过80 %，根据表 3的资源消耗量，该处理方案可以在xc7vx690 t单板上实现16通道7波达角零陷DBF实时处理，或单板复用5组16通道4波达角零陷DBF实时处理。同时表 3列出了不改进矩阵求逆方案，直接使用高斯消元法进行矩阵求逆以生成DBF处理权重时所需的FPGA资源量，第2行为单个高斯消元矩阵求逆模块所需的资源，除LUTRAM外，其余的资源消耗都明显大于LDLT迭代所需的资源。4波达角的情况下，高斯消元矩阵求逆需要344个FPGA时钟，而该过程无法线速进行，这意味着如果要达到与本文LDLT迭代一致的权重更新率的话，需要复用344个高斯消元矩阵求逆模块，所需资源量过大，完全没有实现意义。根据文献[14]中的分析，Φ小于400可以使灵敏度降低不超过0.1 dB，满足此要求需要将高斯消元求逆模块复用4次，此时Φ=336，所需的资源量如第4行所示。16通道4波达角零陷DBF实时处理需要2块xc7vx690t才能实现，性能远差于LDLT迭代。同时，如果时钟主频降低为100 MHz，高斯消元矩阵求逆所需的资源会翻倍，而LDLT迭代的Φ会变为12但消耗资源量不变。

扫描角为27.1°时，多点零陷DBF-SAR方向图如图 8(a)实线所示，虚线方向图为子阵方向图。零陷目标星下点和距离模糊点已在图上表示出来，除零陷目标外产生的多个零陷点为“零陷栅瓣”，在当前仿真参数下主瓣可以避开“零陷栅瓣”，如果未避开会导致主瓣被抑制从而降低成像质量。对回波信号做距离向成像，P₁、P₂、P₃点目标回波强度为0 dB，星下点回波强度为28 dB，图 8(b)为经过SCORE扫描处理后的距离向成像结果，图 8(c)为经过星下点抑制的距离向成像结果，这里的成像结果以星下点抑制后的点目标P₃为基准进行了归一化。可以看出SCORE扫描中星下点回波由旁瓣接收，强度降到-7.5 dB，仍然会影响成像质量；而经过星下点抑制后，星下点回波强度被抑制到-48 dB，相较抑制前降低了40.5 dB。可以看出本实时处理方案可以实现星下点回波抑制，满足高分辨率宽测绘带SAR成像需要。图 8(d)为距离模糊比(range ambiguity to signal ratio, RASR)曲线，用以表征距离向上DBF-SAR系统的模糊性，值为距离模糊信号和主信号强度之比。图中虚线为SCORE扫描时的RASR曲线，实线为距离模糊抑制时的RASR曲线，可以看出距离模糊抑制使RASR降低10~30 dB，距离模糊抑制前后的回波信号强度如图 8(e)和8(f)所示，成像结果以星下点抑制后的点目标P₃为基准进行了归一化。图中为点P₁、P₂、P₃的远端1阶距离模糊信号，可以看出距离模糊信号强度被抑制了24 dB，降低到-50 dB左右，可以看出本实时处理方案可以实现距离模糊抑制，降低距离模糊信号对雷达成像质量的影响。

4 总结

本文论证了多零点DBF处理中的波达矩阵可以进行LDLT分解，提出一套改进的工程上可实现的俯仰向多点零陷DBF实时处理方案，用以在实现SCORE扫描中对星下点回波和距离模糊点回波进行抑制。该方案使用2次DBF处理，基于LCMV方法求2次DBF处理所需的权重，使用LDLT分解迭代来解决求第2次DBF权重时面对的线速矩阵求逆问题。相较于已有方案改进了数据处理结构并大幅减少资源消耗，具备工程可实现性。该方案在实时处理的条件下保持SCORE扫描方向最大增益的同时，实现对星下点和距离模糊点的抑制。本文将该实时处理方案在FPGA上进行仿真实现，验证数字资源使用量和处理时延，确认其满足星载平台资源约束，并通过一个受星下点回波干扰的仿真实例验证了本处理方案可以消除星下点回波干扰，降低距离模糊，为高分辨率宽测绘带星载DBF-SAR系统的工程实现提供了依据。