雷达信号处理技术与系统

雷达信号处理技术与系统
设计
脉冲多普勒雷达信号处理仿真
一、雷达概述
雷达是Radar（Radio Detection And Ranging）的音译词，意为“无线电检测和测距”，即利用无线电波来检测目标并测定目标的位置，这也是雷达设备在最初阶段的功能。

雷达的任务就是测量目标的距离、方位和仰角，还包括目标的速度，以及从目标回波中获取更多有关目标的信息。

典型的雷达系统如图1，它主要由雷达发射机、天线、雷达接收机、收发转换开关、信号处理机、数据处理机、终端显示等设备组成。

收发转换开关天线发射的电磁波
目标
雷达发射机
接收的电磁波
雷达接收机
信号处理机
数据处理机
终端显示
图1 雷达系统框图
雷达发射机产生符合要求的雷达波形，然后经馈线和收发开关由发射天线辐射出去，遇到目标后，电磁波一部分反射，经接收天线和收发开关由雷达接收机接收，然后对雷达回波信号依次进行信号处理、数据处理，就可以获知目标的相关信息。

二、雷达信号
雷达发射信号可以分为连续信号和脉冲信号，常规雷达信号包括非相参脉冲信号、相参脉冲信号、参差变周期脉冲信号、步进频率脉冲信号、线性调频信号、非线性调频信号、相位编码信号等，这里主要介绍常用的线性调频信号，非线性调频信号，相位编码信号等。

1.线性调频信号
为了实现雷达发射能量与分辨率之间的矛盾，线性调频脉冲压缩体制的发
射信号其载频在脉冲宽度内按线性规律变化即用对载频进行调制(线性调频)的方法展宽发射信号的频谱，使其相位具有色散。

LFM （Linear Frequency Modulation ）信号(也称Chirp 信号)的数学表达式为：
)
2(22)(
)(t K
t f j c e T
t rect t s +=π
式中c f 为载波频率，()t
rect T
为矩形信号，即
11()0,t t rect T
T elsewise
⎧ , ≤⎪=⎨⎪ ⎩
B
K T =
，是调频斜率。

于是，信号的瞬时频率为()22c T T f Kt t + -≤≤，根
据K 的正负可以分为两种典型的chirp 信号，如图2所示。

图2 典型的chirp 信号
（a ）up-chirp(K>0)（b ）down-chirp(K<0)
2. 非线性调频信号
非线性调频脉冲信号是指脉内频率调制函数是非线性函数的一类信号。

可以表示为：
))(exp()()(t j t u t x ϕ=
)(t x 的调频函数：∑+∞
=-+=
=1
1
2sin
)()()(n nt
n K B Bt
f T t f τ
πτ
)(t x 的相位函数：∑
⎰+∞
=∞
-+=
=1
2
2sin )(2)(2)(n t
nt
n n K B t B dv v f t τ
πττ
π
πϕ
上式中，)(f T 为)(t x 的群时延，τ和B 分别为非线性调频信号的时宽和带宽，)(n K 为傅里叶级数的系数，实际应用中只取前几项。

3. 相位编码信号
相位编码信号的调制函数是离散的有限状态，属于离散编码信号。

由于相位编码采用伪随机序列，故亦成为伪随机编码信号。

伪随机相位编码信号按相移取值数目分类。

如果相移只限取0、π两个数
值，称之为二相码信号，如巴克码、M 序列码、L 序列等；如果相移可取两个以上的数值，则称之为多相码信号。

如Taylor 多相码、法兰克多相码、赫夫曼序列等。

图3 相位编码信号
三、目标回波仿真
概述
雷达发射机产生线性调频信号，通过天线辐射出去，如果传播过程中遇到目标，就会反射回一部分电磁波，由雷达接收机接收。

这就是回波信号，回波信号中包含有目标的距离，速度，角度等给方面的信息。

目标信号包括期望目标和非期望目标，如图4所示。

图4 目标回波产生
由于目标和雷达之间的距离和相对速度的影响，回波信号会产生一定的延迟，以及多普勒频移。

使用传播响应函数来描述目标回波产生的过程，相对发射电磁波传播响应函数如式3-1：
()
()()()()exp exp 2()a t c r t P A t j t j f t h t w t φπ⎡⎤⎡⎤=⋅⋅⋅⊗+⎣⎦⎣⎦
(),,1()M
t m m m r m P m h t G t t G L δτσ=⎡⎤=⋅-∆⋅⋅⋅⎣⎦∑ （3-1）
式4-2中，M 表示目标个数，()w t 为噪声信号，,t m G 和,r m G 分别为雷达天线发射
幅度增益和接收幅度增益，p L 为传播衰减，()m t τ∆为各目标的延时时间。

各目标的延时时间()m t τ∆满足式3-2：
()()2 2 m m m m R C t R v t C τ⎧⎪∆=⎨-⎪⎩
对静止目标
对运动目标 （3-2）
其中，m v 表示各个目标相对于雷达的速度。

对于雷达天线发射幅度增益和接收幅度增益，采用低旁瓣天线功率方向图进行仿真，使用sin ()c 函数描述天线方向图，即
()max ()sin /2r dB G G c θθθ= （3-3）
其中，max G 为天线最大幅度增益，r θ为目标目标偏离雷达发射方向的角度，3dB θ表示天线3dB 带宽。

● 传播衰减
回波信号在空气中传播，会发生一定的损耗，称之为传播衰减。

主要包括
两个方面，一是大气损耗，二是功率稀释。

大气损耗A L 是雷达工作频率、目标距离和仰角的函数。

雷达工作频率越高，大气损耗越大。

所以，在频率较低的频段(3GHz 以下)，大气损耗可以不予考虑 ，在频率较高的频段，进行选择性的考虑。

在本次实验中，雷达工作在中重频下，并没有考虑大气损耗的影响。

功率稀释是由于天线是向所有方向均匀发射能量的，也就是说天线具有球
形辐射方向图，所以目标处接收到的电磁波能量的功率密度为2
4t
D P P R π=，其
中t P 为雷达发射功率，R 为雷达和目标之间的距离，即单程的功率稀释为
214one way L R π-=
，类似的，双程的功率稀释为()24
1
4two way L R
π-=。

最后，传播衰减要同时考虑大气损耗和功率稀释两方面的影响，即
p A two way L L L -=。

● 调制到中频
雷达是利用物体反射电磁波的特性来发现并确定目标参数的，雷达发射的
信号应该是一个载波受到调制的大功率射频信号。

雷达工作频率是按照雷达的用途来确定的，为了调高雷达系统的工作性能和抗干扰能力，有时要求它能在几个频率上跳变工作或者同时工作。

调制到中频，就是对发射信号乘上一个载频信号，即
()()exp(2)c s t s t j f t π= （3-4）
● 加入噪声
在雷达接收机中，除了目标回波信号之外的任何其他信号都成为噪声。

它包括雷达系统之外的干扰信号和雷达接收机内部产生的热噪声。

热噪声（电子的热骚动）和散射效应噪声（半导体的载流子密度的变化）是雷达接收机中两种主要的内部噪声源。

在本仿真实验中，使用了简化模型，即假设加入的噪声是服从高斯分布的，模型如式（3-5）：
()()+()r t r t n t =，()()+()r i n t n t jn t = （3-5）
其中，()
2()()~0,2r i n n t n t N δ，，2
0n kBT F δ=
四、信号处理
雷达信号处理的流程如下：
下面具体介绍如下。

● 正交双通道采样
正交双通道处理就是中频回波信号经过两个相似的支路分别处理，其差别仅是其基准的相参电压相位差90°，这两路称为：
同相支路(Inphase Channel)——I 支路 正交支路(Quadrature Channel)——Q 支路
传统方法使用的是模拟正交双通道处理，正交I 、Q 通道处理是将接收机输出的中频回波信号分别与正交的两路相参信号混频（采用模拟乘法器），然后进行低通滤波，从而得到I 、Q 两路基带信号，再通过A/D 变换给出同相分量和正交分量的数字量，如图5所示：
图5 正交双通道采样结构图
低通滤波
低通滤波
A/D
I
Q
中频带通信号
)
2sin(0t f π)
2cos(0t f π)]
(2cos[)(0t t f t A θπ+A/D
接收机 正交双通道
采样
匹配滤波
MTI MTD
测距 测速 测角
正交双通道处理的优点(相对于单通道处理)： ➢ 可区分d f ±，以确定目标相对运动方向。

➢ 能消除盲相（单通道MTI 时目标多普勒信号的相位取样对消导致零输出）。

● 匹配滤波
脉冲压缩的目的是集中单个雷达发射信号的所有能量，获的最大输出信噪比。

方法是进行匹配滤波，在接收机中设置一个与发射信号频率相匹配的压缩网络，使经过调制的宽脉冲的回波信号变成窄脉冲，保持良好的距离分辨力。

脉冲压缩网络实际上就是一个匹配滤波器网络。

匹配滤波器是指输出信噪比最大准则下的最佳线性滤波器。

根据匹配理论，匹配滤波器的传输特性：
0)()(*t j e KS H ωωω-= （4-1）
时域表示（冲激响应）为：
)()(0*t t Ks t h -= （4-2）
其中，K 为幅度归一化常数，S (ω) 是发射信号，x (t )是回波信号。

线性调频信号的匹配滤波有两种方法：时域匹配滤波、频域匹配滤波。

时域匹配滤波)()()(t h t x t y ⊗=：滑动滤波器（FIR ）。

运算量大，难以满足实时处理的要求。

频域匹配滤波：傅立叶变换后频谱相乘。

可采用FFT 算法大幅度降低运算量，满足实时处理的要求。

频域匹配滤波：傅立叶变换后频谱相乘，具有简单的频域解析表达式。

可采用FFT 算法大幅度降低运算量，满足实时处理的要求。

通过加窗能够获得很低的旁瓣，如图6所示。

[]()IFFT ()()y t H f R f =⋅，()
()2()exp H f j f K Win f π=⋅
图6 频域匹配滤波算法
● 动目标显示MTI
动目标显示（MTI ）即Moving Target Indication ，是利用MTI 滤波器滤除相应杂波，从而提高目标检测性能。

固定目标频谱的谱线位于脉冲重复频率的整数倍点处，而运动目标回波信号存在多普勒频移，动目标显示滤波器利用运动
目标回波和杂波在频谱上的区别，有效地抑制杂波而提取信号。

最直接的方法是将相邻重复周期的回波信号相减，则固定目标回波由于振幅不变而互相抵消，运动目标回波相减后剩下相邻重复周期振幅变化的部分。

实验中用到的就是这种传统的非递归型一次对消器，即二脉冲对消。

结构如图7：
图7 二脉冲对消结构图
时域方程为：)1()()(--=n x n x n y ，传输函数为：11)(--=z z H ，它是一个单零点系统，零点的位置在1=z ，频率响应为：
)2
cos 2(sin 2sin 21)(T
j T T e e H t j j ωωωωω+=-=-
频率响应如图8，在脉冲重复频率的整数倍点处有凹口，所以固定目标回波在通过MTI 滤波器后将受到很大的抑制，理想状态下，输出为零。

图8 MTI 滤波器频率响应
动目标检测MTD
MTD 也就是一种相参积累和多普勒滤波的结合，相干积累的目的为：1、集中多个脉冲重复周期/调频周期内雷达发射的所有信号所有能量，获取最大输出信噪比。

2、减小目标RCS 起伏对目标检测的影响。

动目标检测（MTD ）即Moving Target Detection ，根据最佳线性滤波理论，在杂波背景下检测运动目标回波，除了杂波抑制滤波器外，还应串接有对脉冲串信号匹配的滤波器。

MTD 利用了回波脉冲串的相参性进行相参积累。

实际工作中，采用一组相邻且部分重叠的滤波器组覆盖整个多普勒频率范围，这就是窄带多普勒滤波器组。

N 个相邻的多普勒滤波器组的实现是由N 个输出的横向滤波器（N 个脉冲和N-1根迟延线）经过各脉冲不同的加权并求和后形成的。

结构如图9：
1
-Z +-X (z)
Y (z)
+
图9 MTD 滤波器组成框图 图10 MTD 滤波器频率响应N=8
设加在第k 个滤波器的第i 个输出端头的加权值为：
[2(1)]
,0,1,
-1j i k N ik w e i N π--==
k 表示标号从0到N-1的滤波器，每一个k 值对应一组不同的加权值，相应地对应一个不同的多普勒滤波器响应。

图10中所示滤波器响应是N=8时加权所得各标记k 的滤波器频率响应，k 取0~7。

该滤波器的频率覆盖范围为0到 f r 。

在仿真实验中，通常是通过快速傅里叶变换FFT 来实现的。

恒虚警检测CFAR
在对雷达回波信号作了脉冲压缩、MTI 、MTD 滤波后，接着就要对目标的存在进行判决：过门限检测。

这里采用的恒虚警概率检测，即CFAR ，将检测门限计算成使雷达接收机能保持恒定的预定虚警率。

检测的原则如下：
10: ()()()() () :()()() H r t s t c t n t r t H r t c t n t =++>⎧⎨=+<⎩
门限门限 （4-3）
图11 CFAR 检测原理图
恒虚警检测的重点就是确定恒虚警检测的门限。

式4-4给出了门限值V T 和虚警概率P fa 之间的关系：
21
2ln T fa V P ψ⎛⎫
=
⎪ ⎪⎝
⎭
（4-4） Tr Tr Tr
Tr 求和输入
输出
0k w 1k w 2k
w 3k
w (1)N k
w -…
k
y
其中，2ψ为噪声的功率，由于噪声的功率是一直变化的，为了保持恒定的
虚警概率，必须依据噪声方差的估计连续更新门限值。

连续改变门限值以保持恒定虚警概率的过程叫做恒虚警概率（CFAR）。

在仿真实验中，假设()()()
r t c t n t
=+，要计算出门限值，就要先估算出噪声的功率。

()()
()()
()
()(
)
()
2
2222
exp
2
n c n c
n t c t n t c t
p n t c t
σσσσ
⎛⎫
++
⎪
⎡⎤
+=⋅-
⎣⎦ ⎪
++
⎝⎭
()
()
()
()()()
22
fa222222
exp exp
22
T
T
U
n c n c n c
r t r t U
P dr t
σσσσσσ
∞
⎛⎫⎛⎫
⎪ ⎪
=⋅-⋅=-
⎪ ⎪
+++
⎝⎭⎝⎭
⎰
由
fa
P恒定，可以推导出：
()
22
fa
2ln
T n c
U P
σσ
=-+
但是噪声的功率是不知道的，必须要对它进行估计。

因为目标信号的幅度值是比较大的，为了消除目标信号对噪声估计的影响，在被检测目标左右，都有保护单元，不参与对噪声功率的估计，原理图如图12所示，估计的噪声功率如式4-5：
24
13
22
22
2143
=()()()()()
t t
n c t t
r t d r t r t d r t t t t t
σσ⎡⎤
++-+-
⎢⎥
⎣⎦
⎰⎰（4-5）
图12 噪声估计
测距
雷达工作时，发射机经天线向空间发射一串重复周期一定的高频脉冲。

如果在电磁波传播的过程上有目标的存在，那么雷达就可以接受到由目标反射回来的回波。

由于回波信号往返于雷达和目标之间，它将滞后于发射脉冲一段时间r t，如图13所示，电磁波的能量是以光速传播的，设目标的距离为R，则传播的距离为光速乘上时间间隔，即
2
r
R =
（4-6） 式中，R 为目标到雷达站的单程距离，单位为m （米），t ∆的单位为s （秒），因子12是考虑到往返的时间延迟。

图13 雷达测距原理
与往返延迟时间T （脉冲周期）对应的距离称为雷达的非模糊距离u R 。

因为接收到的回波可以看成使离它最近的发射脉冲的回波，也可以是前一个发射脉冲的回波，在这种情况下，
2r
ct R = 或 ()2
r c t T R +=
可以看到，该回波就有了距离模糊。

所以，最大无模糊距离必须对应于脉
冲周期的一半，即
22u r
cT c R f =
= （4-7） 其中，T 为脉冲周期，1r f T =为脉冲重复频率。

测速
有些雷达除了要确定目标的位置外，也需要确定运动目标的相对速度。

目标运动的速度可以从测量确定时间间隔的距离变化量R ∆来定，即v R t =∆∆。

这种办法称为目标距离微分法测速，这种测速方法需要较长的时间，且不能测定其瞬时速度。

一般来说，测量的准确度也差，其数据只能作为粗侧用。

测速常使用的方法是多普勒测频法。

我们知道，当目标与雷达站之间存在相对速度时，接收到回波信号的载频相对于发射信号的载频会产生一个频移，这个频移在物理学上称为多普勒频移，它的数值为
r
d f λ
=
（4-8）
式中，d f 为多普勒频移，单位为Hz ；r v 为雷达和目标之间的径向速度，单位为m/s ；λ为载波波长，单位为m 。

当目标向着雷达站运动时，0r v >，回波载频提高；反之0r v <，回波载频降低，如图14。

雷达只需测量出目标回波的多普勒频移d f ，就可以求得目标相对于雷达的相对速度。

（a ）靠近目标 （b ）后退目标
图14 显示多普勒频移的雷达接收信号频谱
多普勒测频法测量速度可靠性比较高，但是容易出现多普勒模糊，脉冲串线谱的包络为sin x x 形，谱线的间隔为脉冲重复频率r f ，如图15所示。

只要预期的多普勒频移小于各个滤波器带宽的二分之一（即一个FFT 门的宽度的二分之一），则多普勒滤波器组就能够解出目标的多普勒频移。

所以，在无模糊的情况下，目标最大多普勒频移为max 2d r f f =，即脉冲雷达的最大无模糊预期目标的相对径向速度max r v 为：
max max ==24
d r
r f f v λλ （4-15）
（a）多普勒得到了分辨（b）频谱移到了下一个多普勒滤波器中，发生模糊
图15 发射和接收波形的频谱和多普勒组
振幅和差式测角
雷达测角的物理基础是电波在均匀介质中传播的直线性和雷达天线的方向性。

测角的方法可以分为振幅法和相位法两大类。

在本次实验中，用到的是属于振幅法的比幅单脉冲测角法，这里只对这种测角方法进行简单介绍。

比幅单脉冲测角法采用两个相同且彼此部分重叠的波束，其方向图如图16所示。

如果目标处在两波束的交叠轴OA方向，则由两波束收到的信号强度相等，否则一个波束收到的信号强度高于另一个。

所以常常称OA为等信号轴，这种测角方法也叫做等信号测角法。

当两个波束收到的回波信号相等时，等信号轴所指的方向即为目标方向。

如果目标处在OB方向，波束2收到的回波比波束1的强，处在OC方向时，波束2的回波较波束1的弱，因此，比较两个波束的强弱就可以判断目标偏离等信号轴的方向，并可用查表的办法估计出偏离等信号轴角度的大小。

图16 比幅单脉冲测角原理
波束1接收到的回波信号
1()
k t
u KFθθ
=-，波束2接收到的回波信号的电
压值为
2()
k t
u KFθθ
=+，其中，()
Fθ为天线电压方向性函数，
t
θ为目标方向
偏离等信号轴0θ的角度，k θ为0θ与波束最大值方向的偏角。

由1u 和2u 可以求得其差值()t θ∆以及和值()t θ∑，在等信号轴附近可以进行近似，即
()[]0
12()
()()2t k t k t t
dF u u K F F k d θθθθθθθθθθ=∆=-=--+≈
()[]120()()2()t k t k t u u K F F F k θθθθθθ∑=+=-++≈
归一化的和差值为：
0()=()t dF F d θθθθθθ=∆
∑ 因为∆∑正比于目标偏离0θ的角度t θ，故可以用它来判断t θ的方向和大小。

五、仿真结果
● 参数
从PD 雷达系统构成图中，可以看到该雷达接收到回波信号分为和、差信号两路，雷达发射信号是线性调频脉冲信号，和波束和差波束进行信号处理的方式是完全相同的，在进行CFAR 后，进行振幅和差式测角，使用和波束处理后的结果进行测距和测速。

实验中， PD 雷达系统及目标的仿真参数设置如下：
表5-1 雷达系统仿真参数表
雷达系统仿真参数
雷达发射功率Pt 50 W 中心频率Fc 1 GHz 脉冲宽度 Tp 100 us 脉冲重复频率 Fr 2 kHz 带宽B 5 MHz 采样率 Fs 10 MHz 天线孔径D
5 m
天线3dB 波束宽度
0.88*λ/D A 、B 波束与天线轴向的夹角
0.8*B/2 波束主瓣初始指向 20 度 CPI 周期数 64 目标参数
目标距离
50 Km 目标平均后向散射截面积 1 m 2 目标方位角 20度 目标速度
100 m/s
● 结果
发射信号（LFM ）：
0.5
1
1.5
2
2.53
3.5
4
4.5
5
x 10
-3
-1-0.500.51Real part
time [sec]m a g n i t u d e
0.5
1
1.5
2
2.53
3.5
4
4.5
5
x 10
-3
-1-0.500.51
Image part
time [sec]
m a g n i t u d e
1.92
1.94 1.96 1.982
2.02 2.04 2.06 2.08 2.1 2.12x 10
-3
-1-0.500.51Real part
time [sec]m a g n i t u d
e
1.98
2
2.02 2.04
2.06
2.08
2.1x 10
-3
-1
-0.500.51Image part
time [sec]
m a g n i t u d e
回波信号：
00.51 1.52 2.53 3.5x 10
5
-1
-0.5
0.5
1
1.5
0.51 1.52 2.53 3.5x 10
5
-1.5
-1
-0.5
0.5
1
1.5
x 10
-6
差波束回波信号
匹配滤波：
0.51 1.52 2.53 3.5x 10
5
00.2
0.40.60.811.21.41.61.8
00.51 1.52 2.53 3.5
x 10
5
1
2
3
4
5
6
7x 10
-7
差波束回波信号
MTI ：
MTD：
CFAR：
测距：49.9890 m
测速：100.3125 m/s
测角：20.1815 度
六、心得
通过对雷达系统仿真与性能评估这门课的学习，对雷达系统有了一定深度的了解，简单学习了雷达系统各个部分的原理，包括雷达发射机，雷达接收机，天线，目标回波，信号处理机以及数据处理机部分。

对脉压，MTI，MTD，CFAR，测距，测速，测角的原理都有了一定的理解。

在此基础上，学习了雷达系统仿真的方法，包括功能级仿真，信号级仿真，分布式交互仿真和半实物仿真。

课程中，主要是使用Matlab编程实现雷达系统各部分的功能。

七、代码
LFM信号
c=3e8; % 光速
Fc=1e9; % 中心频率【Hz】
Tp=50e-6; % 脉冲宽度【微秒】
Fr=2000; % 脉冲重复频率【Hz】
B=1e6; % 带宽【Hz】
Fs=10e6; % 采样率【Hz】
K=B/Tp; % 调频率【Hz】
Tr=1/Fr; % 脉冲重复周期【秒】
CPI=10*Tr; % 仿真持续时间【秒】
Delta_t=1/Fs; % 时域采样点时间间隔【秒】
s_lfm_I=zeros(length(round(CPI*Fs)),1); % 定义信号实部数组
s_lfm_Q=zeros(length(round(CPI*Fs)),1); % 定义信号虚部数组
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% 数据流仿真方式 %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t_set_Tp=0:Delta_t:Tp; % 一个脉冲内的时间采样点数组
s_lfm_Tp_I=cos(pi*K*(t_set_Tp-Tp/2).^2); % 脉冲信号实部
s_lfm_Tp_Q=sin(pi*K*(t_set_Tp-Tp/2).^2); % 脉冲信号虚部
N=length(s_lfm_Tp_I);
for n=1:10
pulse_start=round((n-1)*Tr*Fs)+1; % 第n个脉冲起始采样点
pulse_stop=round((n-1)*Tr*Fs)+N; % 第n个脉冲结束采样点
s_lfm_I(pulse_start:pulse_stop)=s_lfm_Tp_I;
% 将第n个脉冲信号的实部添加到其在信号数组中的对应位置；
s_lfm_Q(pulse_start:pulse_stop)=s_lfm_Tp_Q;
% 将第n个脉冲信号的虚部添加到其在信号数组中的对应位置；end
figure
plot(0:Delta_t:Delta_t*(length(s_lfm_I)-1),s_lfm_I)
axis([0 CPI -1.25 1.25])
title('Real part')
xlabel('time [sec]')
ylabel('magnitude')
信号处理
clear %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% 雷达系统仿真参数 %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c=3e8; % 光速
k=1.38e-23; % 玻尔兹曼常数
Pt=50; % 发射功率【W】
Fc=1e9; % 中心频率【Hz】
Wavelength=c/Fc; % 工作波长【m】
Tp=100e-6; % 脉冲宽度【微秒】
Fr=2000; % 脉冲重复频率【Hz】
B=5e6; % 带宽【Hz】
Fs=10e6; % 采样率【Hz】
F=10^(6.99/10); % 噪声系数
K=B/Tp; % 调频率【Hz】
Tr=1/Fr; % 脉冲重复周期【秒】
Delta_t=1/Fs; % 时域采样点时间间隔【秒】
D=5; % 天线孔径【m】
Ae=1*pi*(D/2)^2; % 天线有效面积【m^2】
G=4*pi*Ae/Wavelength^2; % 天线增益
BeamWidth=0.88*Wavelength/D; % 天线3dB波束宽度【deg】BeamShift=0.8*BeamWidth/2; % A、B波束与天线轴向的夹角【deg】Theta0=20*pi/180; % 波束主瓣初始指向【度】
Wa=0;2*pi/1; % 天线波束转速【rad/sec】
Num_Tr_CPI=64+1; % CPI周期数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% 目标仿真参数 %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R_set=[50e3]; % 目标距离【m】
RCS=[1]; % 目标平均后向散射截面积【m^2】
Theta_target_set=[20]*pi/180; % 目标方位角【deg】
V_set=[100]; % 目标速度【m/s】
RCS_Ground_0=10^(-30/10); % 地面单位面积后向散射截面积【m^2】
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s_Sigma=zeros(round(Tr*Fs),Num_Tr_CPI); % 定义和通道信号数组
s_Delta=zeros(round(Tr*Fs),Num_Tr_CPI); % 定义差通道信号数组
t_set_Tp=(0:Delta_t:Tp)'; % 一个脉冲内的时间采样点数组
s_lfm=exp(j*pi*K*(t_set_Tp-Tp/2).^2); % 脉冲复信号
N=length(s_lfm); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 仿真目标回波信号 %%%%%
for No_PRI=1:Num_Tr_CPI
Theta_bp=Theta0+Wa*No_PRI*Tr; % 波束主瓣指向【度】
for No_target=1:1
delay_target=2*(R_set(No_target)-V_set(No_target)*No_PRI*Tr)/c;
% 目标时延【sec】
RVP=-2*pi*Fc*delay_target;
% 目标回波视频检波剩余相位
Gt_A=G*(sinc((Theta_target_set(No_target)-(Theta0-BeamShift))/BeamWidth)).^2; Gr_A=G*(sinc((Theta_target_set(No_target)-(Theta0-BeamShift))/BeamWidth)).^2; % 波束A在目标方向上的增益
Gt_B=G*(sinc((Theta_target_set(No_target)-(Theta0+BeamShift))/BeamWidth)).^2; Gr_B=G*(sinc((Theta_target_set(No_target)-(Theta0+BeamShift))/BeamWidth)).^2; % 波束B在目标方向上的增益
Lp=1./((4*pi)^2*R_set(No_target).^4);
% 目标处的电磁波传播损耗
Magnitude_echo_A=sqrt(Pt*RCS(No_target)...
*(Gt_A+Gt_B)*Gr_A*Wavelength^2/(4*pi)*Lp);
% 波束A中目标回波幅度
Magnitude_echo_B=sqrt(Pt*RCS(No_target)...
*(Gt_A+Gt_B)*Gr_B*Wavelength^2/(4*pi)*Lp);
% 波束B中目标回波幅度
Echo_start=round((delay_target)*Fs); % 目标回波起始采样点
Echo_stop=Echo_start+N-1; % 目标回波结束采样点
s_Sigma(Echo_start:Echo_stop,No_PRI)=s_Sigma(Echo_start:Echo_stop,No_PRI)+... (Magnitude_echo_A+Magnitude_echo_B)*exp(j*RVP)*s_lfm;
s_Delta(Echo_start:Echo_stop,No_PRI)=s_Delta(Echo_start:Echo_stop,No_PRI)+... (Magnitude_echo_A-Magnitude_echo_B)*exp(j*RVP)*s_lfm;
% 将第n个目标的复回波的添加到其在信号数组中的对应位置；
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 仿真热噪声信号 %%%%%
n_Sigma=sqrt(k*B*F*290/2)*( randn(size(s_Sigma))+j*randn(size(s_Sigma)) );
n_Delta=sqrt(k*B*F*290/2)*( randn(size(s_Delta))+j*randn(size(s_Delta)) );
s_Sigma=s_Sigma+n_Sigma;
s_Delta=s_Delta+n_Delta;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
as=reshape(s_Sigma,1,size(s_Sigma,1)*size(s_Sigma,2));
bs=reshape(s_Delta,1,size(s_Delta,1)*size(s_Delta,2));
figure;
plot(real(as));
title('和波束回波信号');
figure;
plot(real(bs));
title('差波束回波信号');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 匹配滤波（脉冲压缩) %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Num_sample=round(Tr*Fs);
% 每个脉冲发射回波信号的长度
frange=((-Fs/2+Fs/(2*Num_sample)):Fs/Num_sample:(Fs/2-Fs/(2*Num_sample))).'; % 每个脉冲发射回波信号变换到频域后每个频域采样点对应的频率
Win=hamming(Num_sample);
% 窗函数（降低旁瓣）
H_match=exp(1i*pi*frange.^2/K).*Win;
% 匹配滤波函数
S_Sigma_r=fftshift(fft(s_Sigma,[],1),1);
S_Delta_r=fftshift(fft(s_Delta,[],1),1);
% 将每个脉冲发射回波信号变换到频域
S_Sigma_r=S_Sigma_r.*(H_match*ones(1,Num_Tr_CPI));
S_Delta_r=S_Delta_r.*(H_match*ones(1,Num_Tr_CPI));
% 匹配滤波
s_Sigma_rc=ifft(ifftshift(S_Sigma_r,1),[],1);
s_Delta_rc=ifft(ifftshift(S_Delta_r,1),[],1);
% 将每个脉冲发射回波信号变换回时域
as=reshape(s_Sigma_rc,1,size(s_Sigma_rc,1)*size(s_Sigma_rc,2));
bs=reshape(s_Delta_rc,1,size(s_Delta_rc,1)*size(s_Delta_rc,2));
figure;
plot(abs(real(as)));
title('和波束回波信号');
figure;
plot(abs(real(bs)));
title('差波束回波信号');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 两脉冲MTI %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s_Sigma_mti=zeros(size(s_Sigma_rc));
s_Delta_mti=zeros(size(s_Delta_rc));
for No_Pulse=2:Num_Tr_CPI
for No_tr=1:Num_sample
s_Sigma_mti(No_tr,No_Pulse)=s_Sigma_rc(No_tr,No_Pulse)...
-s_Sigma_rc(No_tr,No_Pulse-1);
s_Delta_mti(No_tr,No_Pulse)=s_Delta_rc(No_tr,No_Pulse)...
-s_Delta_rc(No_tr,No_Pulse-1);
end
end
figure;
mesh(abs(real(s_Sigma_mti)));
title('和波束MTI结果');
figure;
mesh(abs(real(s_Delta_mti)));
title('差波束MTI结果');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 多普勒滤波（脉冲积累） %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
win_doppler=hamming(Num_Tr_CPI-1).';
S_Sigma_a=fftshift(fft(s_Sigma_mti(:,2:Num_Tr_CPI).*(ones(Num_sample,1)*win_dopple r),[],2),2);
S_Delta_a=fftshift(fft(s_Delta_mti(:,2:Num_Tr_CPI).*(ones(Num_sample,1)*win_doppler), [],2),2);
% 将每个脉冲发射回波信号变换到多普勒域
figure;
mesh(abs(real(S_Sigma_a)));
title('和波束MTD结果');
figure;
mesh(abs(real(S_Delta_a)));
title('差波束MTD结果');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% CFAR（恒虚警检测） %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pfa=1e-6;
% 虚警概率
No_dopper_channel_set=1:Num_Tr_CPI-1;
% 进行CFAR检测的多普勒通道序号
Num_ReservedCell=3;
% 待测单元附近的保护单元长度（前或后）
Num_TestWin=50;
% 待测单元附近的统计杂波噪声功率的窗口长度（前或后）
N1=Num_TestWin+Num_ReservedCell;
m=1;
for No_dopper_channel=No_dopper_channel_set
n=1;
for No_tr=1+N1:Num_sample-N1
Tranning_set=[(No_tr-N1:No_tr-Num_ReservedCell),...
(No_tr+N1:No_tr+Num_ReservedCell)];
Power_noise_clutter=mean(abs(S_Sigma_a(Tranning_set,No_dopper_channel)).^2); % 统计杂波噪声功率
Threshold=sqrt(2*Power_noise_clutter*log(1/Pfa));
% 检测门限
if abs(S_Sigma_a(No_tr,No_dopper_channel))>=Threshold
S_output(n,m)=1;
else
S_output(n,m)=0;
end
n=n+1;
end
m=m+1;
end
figure;
mesh(abs(real(S_output)));
title('CFAR结果');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 目标距离、多普勒粗测 %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
No_target=0;
for No_dopper_channel=No_dopper_channel_set
if sum(S_output(:,No_dopper_channel))>0
for No_tr=4:Num_sample-2*N1-3
condition_1=S_output(No_tr,No_dopper_channel)==1;
condition_2=sum(S_output(No_tr-3:No_tr-1,No_dopper_channel))==0;
condition_3=sum(S_output(No_tr-3:No_tr+3,No_dopper_channel-1))==0;
% 判决是否为一个新目标的条件
if condition_1 && condition_2 && condition_3
No_target=No_target+1;
Target_Doppler_No(No_target)=No_dopper_channel;
Target_Range_No(No_target)=No_tr+N1;
%记录每个目标的距离序号和多普勒序号
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 目标距离测量 %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for No_target=1:length(Target_Range_No)
s=ifft(fft(S_Sigma_a(:,Target_Doppler_No(No_target))),10*Num_sample);
% 距离向插值细化
[vmax,pmax]=max(abs(s((10*Target_Range_No(No_target)-
49):(10*Target_Range_No(No_target)+49))));
Target_Range(No_target)=((10*Target_Range_No(No_target)-49)+pmax-
1)/(Fs*10)*c/2-(Tp-Delta_t)/2*c/2;
% 利用细化后的距离向峰值点测量目标距离
end
disp(Target_Range/1e3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 目标速度测量 %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for No_target=1:length(Target_Range)
s=fftshift(fft(ifft(ifftshift(S_Sigma_a(Target_Range_No(No_target),:))),10*(Num_Tr_CPI-1)));
% 多普勒域插值细化
[vmax,pmax]=max(abs(s((10*Target_Doppler_No(No_target)-
100):(10*T arget_Doppler_No(No_target)+100))));
Target_Doppler(No_target)=-Fr/2+(10*Target_Doppler_No(No_target)-100+pmax-
1)*(Fr/(10*(Num_Tr_CPI-1)));
% 利用细化后的多普勒频谱峰值点测量目标距离
end
disp(Target_Doppler*Wavelength/2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 单脉冲测角 %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Theta_set=-BeamWidth/2:BeamWidth/100:BeamWidth/2;
% 查表范围
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%计算不同角度的和差比值 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=1;
for Theta=Theta_set
G_A=G*(sinc((Theta+BeamShift)/BeamWidth)).^2;
G_B=G*(sinc((Theta-BeamShift)/BeamWidth)).^2;
% 计算波束A、B的方向图
Sigma_G=G_A+G_B;
Delta_G=G_A-G_B;
% 计算和、差波束的方向图
Delta_2_Sigma_ratio(n)=Delta_G/Sigma_G;
% 计算和差比
n=n+1;
end
for No_target=1:length(Target_Range_No)
Delta_target=S_Delta_a(T arget_Range_No(No_target),Target_Doppler_No(No_target));
Sigma_target=S_Sigma_a(Target_Range_No(No_target),Target_Doppler_No(No_target)); angle_Sigma=angle(Sigma_target);
angle_Delta =angle(Delta_target);
% 求和差通道输出相角差
if abs(angle_Sigma-angle_Delta)>pi/2
Sign_Delta_Sigma=-1;
else
Sign_Delta_Sigma=1;
end
% 确定目标角度正负
Delta_2_Sigma_ratio_T arget(No_target)=abs(Delta_target/Sigma_target);
% 求和差比
[min_v,p_theta]=min(abs(Sign_Delta_Sigma*Delta_2_Sigma_ratio_Target(No_target)-Delta_2_Sigma_ratio));
Theta_targe(No_target)=Theta_set(p_theta)+Theta0;
% 查表求目标角度
end
disp(Theta_targe*180/pi)。