技术领域
[0001] 本
发明属于雷达技术领域，具体的说是一种用于雷达对低仰角目标的波达方向估计方法。
背景技术
[0002] 目标
定位是雷达的一项基本功能。当探测低仰角目标时，雷达波束打地，使得目标的直达波和地(海)面反射的多径反射波在天线波束主瓣内
叠加，且它们的多普勒
频率及回波时延近似相等，这组强相关的
信号同时被雷达天线接收，使得雷达角度估计系统不能正确的估计出目标的波达方向，给目标定位带来了困难。为了解决这个问题，
现有技术有以下三种方法：
[0003] 1.吴向东等人在“一种基于线性预处理的米波雷达低仰角处理
算法，
电子学报，2006，34(9)：1668-1671”文章中，提出先对接收数据进行差分预处理，再通过多重信号分类MUSIC算法估计目标的波达方向。
[0004] 2.赵光辉等人在“基于交替投影的DOA估计方法及其在米波雷达中的应用，电子与信息学报，2008，30(1)：224-227”文章中，提出利用广义MUSIC算法估计低仰角目标的波达方向。
[0005] 3.赵永波等人在“雷达低角
跟踪环境下的最大似然波达方向估计方法，电子学报，2004，32(9)：1520-1523”的文章中，提出了一种
时空级联最大似然算法，即先进行多普勒频率估计和滤波，然后再利用最大似然算法估计目标的波达方向。
[0006] 上述方法1和方法2都是基于数据协方差矩阵的方法，由于没有充分利用雷达可获取的信息，因而对
信噪比要求比较高，且不能用于单样本数情况下的波达方向估计。
[0007] 方法3虽然具有低信噪比
门限和所需快拍数少的优点，但是当阵元较多和搜索角度间隔较小时，角度搜索的运算量将会很大，这给波达方向估计系统的工程实现带来了困难。
发明内容
[0008] 本发明的主要目的在于克服已有方法的缺点，提供一种无需进行角度搜索的雷达低仰角目标的波达方向估计方法，以在保证波达方向估计
精度的情况下，减小运算量，便于波达方向估计系统的工程实现。
[0009] 为实现上述目的，本发明采用的技术方案是：充分利用脉冲雷达回波信号时域和频域的信息对雷达回
波数据进行积累，然后利用均匀线阵的结构特点构造矩阵束，最后直接得到目标波达方向的解析解，不需要角度搜索，具体实现步骤包括如下：
[0010] 1.雷达发射一组带宽不超过
载波频率1％的窄带相参脉冲信号；
[0011] 2.通过均匀线阵天线接收雷达回波数据，并将其混频到基带进行离散
采样；
[0012] 3.对离散采样后的雷达回波数据进行波束形成，波束指向0°；
[0013] 4.根据波束形成后的数据，测量出目标所在的距离采样单元和多普勒通道；
[0014] 5.根据测量得到的目标所在的距离采样单元和多普勒通道，对离散采样后的雷达回波数据重新进行积累，使每个接收阵元只获得一个样本数；
[0015] 6.利用重新积累后的数据，构造矩阵束：
[0016] 其 中，γ 为 一 个 任 意 复 数，* T
() 表示矩阵共轭， () 表示
矩阵转置，qk为第k个阵元对离散采样后的雷达回波数据重新进行积累的输出值，k＝1，
2，...，L+1，L为自由参数，L取值为M/3～M/2间的整数，M是接收雷达回波数据所采用的均匀线阵天线的阵元个数；
[0017] 7.采用广义特征值分解的总体最小二乘法对构造的矩阵束求解目标的波达方向：
[0018] (7a)将 奇异值分解为 其中，()H表示矩阵共轭转置，∑2s为 中两个最大的奇异值所组成的对角阵，∑2n为 中其他奇异值组成的对角阵，U2s由 中两个最大的大奇异值所对应的左奇异向量组成，V2s由 中两个最大的奇异值所对应的右奇异向量组成，U2n由 中其他奇异值所对应的左奇异向量组成，V2n由 中其他奇异值所对应的右奇异向量组成；
[0019] (7b)将矩阵 先左乘 再右乘V2s，得到新的矩阵束：
[0020] (7c)令 为广义特征值矩阵，()-1表示矩阵求逆，对该广义特征值矩阵Φ进行特征值分解，得到特征值γi(i＝1，2)；
[0021] (7d)根据特征值γi，计算目标的波达方向和多径反射波的波达方向的角度值λ为发射信号载波
波长，d为接收雷达回波数据所采用的均匀线阵的阵元间距，arg()表示取
相位，arcsin()表示取反正弦；
[0022] (7e)将计算出的θi(i＝1，2)中的正角度值作为目标的波达方向。
[0023] 本发明与现有技术相比具有以下优点：
[0024] 1.本发明充分利用脉冲雷达可获取的目标距离信息和多普勒信息，通过从雷达回波数据中提取目标回波信号，进行
脉冲压缩和多普勒滤波，以提高信噪比。
[0025] 2.本发明不通过构造空间谱函数进行角度搜索获得目标的波达方向，而是针对均匀线阵的阵元等间隔分布的结构特点，通过构造矩阵束和求解矩阵束的广义特征值，进而得到目标波达方向的解析解，免去了角度搜索带来的巨大运算量。
[0026] 3.本发明将2(M-L)×L维的矩阵束 变换为2×2维的新的矩阵束使一个大维数矩阵的特征值求解的问题转化为一个小维数矩阵的特征值求
解的问题，然后通过广义特征值得到目标的波达方向，大大降低了求解目标波达方向解析解所需的运算量。
[0027] 理论分析和仿真结果表明，本发明与现有技术相比，在保证波达方向估计精度的情况下运算量小。
附图说明
[0028] 图1是本发明使用的雷达低仰角目标回波的多径
几何模型；
[0030] 图3是用本发明方法与现有时空级联最大似然算法仿真得到的目标波达方向估计精度随信噪比变化图。
具体实施方式
[0031] 参照图1，本发明使用的雷达低仰角目标回波的多径几何模型，包括一个垂直放置的均匀线列天线和一个高度为ht的目标，其中，均匀线阵天线作为雷达的接收天线，天线的阵元数为M，阵元的间距为d，天线的中心高度为ha，目标与雷达的直线距离为Rd，目标回波经过地(海)面反射后到达雷达的距离为Rs，目标的波达方向和多径反射波的波达方向分别为θd和θs，假设
水平方向以上的波达方向为角度的正方向。
[0032] 参照图2，本发明结合图1的多径几何模型进行雷达低仰角目标的波达方向估计，具体步骤如下：
[0033] 步骤1，雷达发射一组带宽不超过载波频率1％的窄带相参脉冲信号。
[0034] 雷达以Tr为脉冲重复周期，连续发射Nd个带宽不超过载波频率1％的相参脉冲信号：s(t)＝rect(t)·g(t)·exp(j2πfot)，式中忽略了相参脉冲信号的初相，t表示时间，To为脉冲宽度，·表示数乘，g(t)为发射信号复包络， fo为载波频率。
[0035] 步骤2，通过均匀线阵天线接收雷达回波数据，并将其混频到基带进行离散采样。
[0036] (2a)采用均匀线阵天线接收雷达回波数据；
[0037] (2b)将雷达回波信号混频到基带；
[0038] (2c)对混频到基带后的数据进行离散采样，使在第m(m＝1，2，…，M)个阵元处理得到的雷达回波数据为：
[0039] xm((i-1)Tr+nTs)＝rect((i-1)Tr+nTs-noTs)·g((i-1)Tr+nTs-τdm)·[0040] exp(-j2πfoτdm)·exp(j2πfd((i-1)Tr+nTs))
[0041] (1)
[0042] +ρ·rect((i-1)Tr+nTs-noTs)·g((i-1)Tr+nTs-τsm)·
[0043] exp(-j2πfoτsm)·exp(j2πfd((i-1)Tr+nTs))
[0044] 式中，i＝1，2，…，Nd表示第i个脉冲重复周期，n为距离采样单元，Ts为离散采样周期，fd为目标运动引起的多普勒频移，ρ为表面反射系数，τdm＝τo+τm(θd)，τsm＝τo+τm(θs)+Δτ，τ0为阵列参考点接收目标直达波产生的距离时延，Δτ为目标直达波和多径反射波的波程差产生的时延， 为第m个阵元接收到的目标直达波相对于阵列参考点的波程差产生的时延， 为第m个阵元接收到
的目标多径反射波相对于阵列参考点的波程差产生的时延，c为光速，no为目标所在的距离采样单元。
[0045] 考虑到发射的脉冲信号是带宽不超过载波频率1％的窄带信号，因此有g((i-1)Tr+nTs-τdm)≈g((i-1)Tr+nTs-τsm)≈g((i-1)Tr+nTs-τo)，故将式(1)重写如下：
[0046]
[0047] 式中，am(θd)＝exp(-j2πfoτm(θd))，am(θs)＝exp(-j2πfoτm(θs))，[0048]
[0049] 步骤3，对离散采样后的回波信号进行波束形成，波束方向指向0°。
[0050] 对每个阵元得到的离散采样后的雷达回波数据xm((i-1)Tr+nTs)进行加权求和，得到波束形成的结果如下：
[0051]
[0052] 式中，wm为对第m个阵元所加的权值，由于要求波束指向0°，所以wm＝1。步骤4，根据波束形成后的数据，测量出目标所在的距离采样单元和多普勒通道。
* *
[0053] (4a)以h(nTs)＝g(-nTs)为脉冲压缩的
滤波器系数，() 表示共轭，对波束形成后的数据进行脉冲压缩，得到脉冲压缩后的结果为：η((i-1)Tr+nTs)；
[0054] (4b)将距离采样单元n相同的所有脉冲压缩结果η((i-1)Tr+nTs)组成一个向量[η(nTs)，η(Tr+nTs)，…，η((Nd-1)Tr+nTs)]，采用快速
傅立叶变换方法，对其进行相参积累；
[0055] (4c)对相参积累后的数据采用恒虚警检测方法，得到目标所在的距离采样单元no和多普勒通道Nx。
[0056] 步骤5，根据测量得到的目标所在的距离采样单元no和多普勒通道Nx，对采样后的雷达回波数据重新进行积累。
[0057] (5a)根据己测得的目标所在距离采样单元no，从每个阵元离散采样后的雷达回波数据中提取目标回波信号，即提取xm((i-1)Tr+nTs)中所有满足的距离采样单元的数据；
[0058] (5b)对提取的目标回波信号进行脉冲压缩；
[0059] (5c)在每个脉冲重复周期内，从脉冲压缩后的数据中取一个幅度最大的值，使第m个阵元得到一个1×Nd维信号矢量：
[0060] Ym＝[ym(noTs)，ym(Tr+noTs)，…，ym((Nd-1)Tr+noTs)](4)
[0061] 式中，
[0062] Q＝B·exp(-j2πfoτo)·exp(j2πfdnoTs)，
[0063] B为目标回波信号脉冲压缩后得到的幅度值；
[0064] (5d)根据己测得的目标所在的多普勒通道Nx，以[1，exp(-j2πfxTr)，…，exp(-j2π(Nd-1)fxTr)]为滤波器系数， 对Ym进行多普勒滤波，使第m个阵元的多普勒滤波输出值为
[0065]
[0066]
[0067] 式中， 因此，M个阵元可以获得M个输出值。
[0068] 步骤6，利用重新积累后的数据构造矩阵束。
[0069] 经多普勒滤波后，每个阵元将只获得一个样本数，此时大量的基于协方差矩阵的方法己不再适用，而最大似然算法的运算量又过大。本发明通过构造如下两个矩阵，利用矩阵束的方法来进行波达方向估计。
[0070]
[0071]
[0072] 式中， k＝1，2，...，L+1，L为自由参数，()T表示矩阵转置，根据Sarkar T K等人在“Using the matrix pencil method to estimate the parameters of a sum of complex exponentials，IEEE Antennas and Propagation Magazine，1995，37(1)：48-54”中的分析，L取值为M/3～M/2间的整数。于是，得到如下矩阵束：
[0073]
[0074]
[0075]
[0076]
[0077]
[0078] 式中，以第一个阵元为参考点，则a1(θ)＝1， 为载波波长，diag[]表示矩阵对角化，I2为2×2维的单位矩阵，γ为一个任意复数。
[0079] 步骤7，采用广义特征值分解的总体最小二乘法对构造的矩阵束求解目标的波达方向。显然，当 或 时，式(9)的秩将为1，因此可利用基于广义特征值分解的方法来估计目标的波达方向。本发明采用张贤达所著的“矩阵论分析与应用，北京：清华大学出版社，2004，9”中给出的广义特征值分解的总体最小二乘法来减小噪声的影响和运算量。
[0080] (7a)将 奇异值分解为
[0081]
[0082] 其中，()H表示矩阵共轭转置，∑2s为 中两个最大的奇异值所组成的对角阵，∑2n为 中其他奇异值组成的对角阵，U2s由 中两个最大的大奇异值所对应的左奇异向量组成，V2s由 中两个最大的奇异值所对应的右奇异向量组成，U2n由 中其他奇异值所对应的左奇异向量组成，V2n由 中其他奇异值所对应的右奇异向量组成；
[0083] (7b)将矩阵 先左乘 再右乘V2s，在不改变广义特征值的条件下，得到新的矩阵束如下：
[0084]
[0085] 这样，将式(8)中2(M-L)×L维矩阵束的广义特征值分解的问题转
化成了2×2维矩阵束的广义特征值分解的问题，且变换中只需要一次奇异值分解；
[0086] (7c)令 为广义特征值矩阵，()-1表示矩阵求逆，该广义特征值矩阵Φ的特征值即为式(14)中新的矩阵束的广义特征值，对该广义特征值矩阵Φ进行特征值分解，得到特征值γi(i＝1，2)；
[0087] (7d)根据特征值γi，计算目标的波达方向和多径反射波的波达方向的角度值arg()表示取相位，arcsin()表示取反正弦；
[0088] (7e)由于目标的波达方向角度为正，而多径反射波的波达方向角度为负，所以将计算出的θi(i＝1，2)中的正角度值作为目标的波达方向。
[0089] 本发明的效果可通过以下理论分析和仿真实验进一步说明：
[0090] 一.运算量分析
[0091] 本发明的运算量远远小于时空级联最大似然算法的运算量。由于本发明的步骤1、2、3、4和5与时空级联最大似然算法基本类似，所以下面只考虑步骤6和7的运算量差别。由上面分析可知，本发明的运算量主要集中在2(M-L)×L维矩阵 的奇异值分解上，不
3 3 2 2
需要额外的搜索过程，运算量(次复乘)约为8(M-L)+L+4(M-L)L+2(M-L)L，当自由参数L
3
＝M/2时，运算量(次复乘)小于2M。而时空级联最大似然算法需要一维搜索，在每个搜索角度下需要一次投影矩阵和最大似然函数值的计算，搜索次数与精度需求有关，假设搜
2
索次数为Ω，则时空级联最大似然算法的运算量约为Ω·(3M+9M+8)，而一般Ω远远大于M，所以本发明的运算量要明显低于时空级联最大似然算法。
[0092] 二.仿真结果
[0093] 1.仿真条件
[0094] 仿真中，接收天线采用M＝16个阵元的均匀线性阵列，阵元间距d＝λ/2，雷达架高ha＝100m，载波频率为fo＝150MHz，发射信号复包络g(t)为线性调频信号，脉冲宽度To＝500μs，调频带宽为500kHz，脉冲重复周期Tr＝4ms，相干脉冲积累数为Nd＝32，
采样频率fs＝1MHz，反射系数ρ＝0.9exp(jπ)，自由参数L＝8。假设目标距天线阵中心的距离Rd＝100km，目标多普勒频率fd＝100Hz，目标高度ht＝3580m，则目标的波达方向θd≈1.99°，多径反射波的波达方向θs≈-2.11°。假设每个阵元的接收噪声均为独立同分布的高斯白噪声。目标波达方向估计的精度定义为 为第n次实验得到的估计值，MC为Monte-Carlo实验的总次数，以下仿真实例中，MC＝200。
[0095] 2.仿真内容
[0096] 仿真中，将本发明方法与现有时空级联最大似然算法得到的目标波达方向估计精度进行比较，以说明本发明基本不损失现有技术可达到的波达方向估计精度。本发明方法与现有时空级联最大似然算法得到的目标波达方向估计精度随信噪比变化图见图3。从图3中看出，在低信噪比情况下，本发明与时空级联最大似然算法获得的波达方向估计精度的差别小于0.007°；随着信噪比的增加，本发明获得的波达方向估计精度越来越接近时空级联最大似然算法的精度；在15dB以上时，本发明获得的波达方向估计精度超过了时空级联最大似然算法的精度，这是因为时空级联最大似然算法需要角度搜索，搜索步长限制了其估计精度，仿真中采用的搜索步长为0.01°，而本发明是直接求解波达方向的解析解，因此获得的估计精度更高。
