CN112578352B - 基于fpga的箔条云雷达回波实时模拟方法 - Google Patents

Abstract

本发明公开了一种基于FPGA的动态箔条云雷达回波模拟方法，主要解决现有技术模拟箔条云雷达回波可移植性差，且只模拟单发箔条云回波的问题。其通过模拟运行在FPGA平台上，并对计算过程进行硬件优化来实现：1)计算箔条的初始姿态、每根箔条的初始直角坐标和每根箔条的雷达散射截面值；2)根据箔条高度计算得到大气密度；3)计算出箔条在一个采样时间点内的位移增量，进而计算出每根箔条的直角坐标；4)根据箔条位置的直角坐标与雷达坐标计算两者距离，根据该距离和箔条雷达散射截面值计算箔条云的雷达信号回波，完成一个时间点的计算；重复2)～4)得到不同时刻箔条云的回波。本发明可移植性好，可用于模拟多发箔条云的雷达回波。

Description

基于FPGA的箔条云雷达回波实时模拟方法

技术领域

本发明属于雷达技术领域，尤其涉及一种雷达对箔条云的回波实时模拟方法，可用于获得不同时刻箔条云中各个箔条的动态回波。

背景技术

箔条云是一种使用最为广泛的无源干扰，由于它具有滞空时间长、作用面积大、回波强的特点，在防空、反舰等领域被广泛应用。箔条云回波计算主要分为三个部分，分别是箔条运动扩散计算、箔条雷达散射截面计算和箔条回波计算。其运动扩散与箔条释放平台的箔条抛洒速度、抛洒方向、箔条所在环境的风力、空气密度和粘度相关，是影响箔条的雷达回波特性的主要因素之一。箔条云内的箔条雷达散射截面统计值受箔条云内箔条总数、入射波频率因素影响。其回波特性与雷达的工作体制以及工作波形有关。

在已有的技术中，大多数箔条云回波模拟都是基于PC平台进行的。中国专利CN108020819A，公开了一种弹载雷达跟踪到运动目标后的箔条云干扰回波模拟方法，该方法首先计算每个时刻的箔条云质心坐标，再求出每个距离门上的箔条数；然后得到单根箔条的天线增益、幅度、频率，求得单根箔条的干扰调制值；最后通过叠加所有箔条的回波值得到整个箔条云的回波。此专利实施是基于电脑平台，资源消耗大，系统结构复杂，此外只能对单发箔条云的回波进行仿真，不能对多发箔条的情况进行仿真。中国专利CN110850376A，公开了一种箔条干扰模拟方法及箔条干扰模拟器，该专利首先基于箔条电磁特性进行离线仿真，产生箔条干扰滤波器系数；之后在线产生具有箔条干扰特性的回波系数；最后产生箔条干扰模拟信号。此专利基于箔条的回波系数生成箔条回波，未考虑箔条的实际运动情况，不能准确反映箔条在空中扩散时雷达回波的变化情况；此外，基于电脑平台的仿真资源消耗大、移动性较差。

发明内容

本发明的目的在于针对上述现有技术的不足，提供一种基于FPGA的箔条云雷达回波实时模拟方法，以在准确反映箔条云在空中运动扩散的条件下，实现对多发箔条的雷达回波进行仿真，并减小硬件资源消耗。

为实现上述目的，本发明包括如下步骤：

(1)计算箔条的初始姿态，包括原点到箔条的连线与正z轴之间的天顶角θ、箔条在xOy平面的投影线与正x轴之间的方位角

箔条转动角度ζ、水平扩散速度v，并将这些计算结果保存到RAM中；

(2)计算每个箔条的初始直角坐标Lx₀(i)、Ly₀(i)、Lz₀(i)，分别保存到RAM中，其中i为箔条的序号；

(3)计算每个箔条的雷达散射截面RCS值，把结果保存到RAM中；

(4)计算每条箔条所在高度z(i)处的大气密度ρ(i)；

(5)计算每条箔条受力引起的三个方向上的位移增量dx₁(i)、dy₁(i)、dz₁(i)；

(6)计算箔条螺旋转动引起的水平位移增量dx₂(i)、dy₂(i)；

(7)根据(5)得到的箔条位移增量和(6)得到的位移增量，计算本时间采样点的箔条直角坐标x_t(i)、y_t(i)、z_t(i)；

(8)根据雷达的空间坐标和每个箔条的空间坐标计算雷达到每个箔条的距离R_t(i)；

(9)根据步骤(3)得到的每根箔条的雷达散射截面RCS值和步骤(8)得到的雷达到箔条的距离R_t(i)，计算箔条的窄带信号回波S(t,i)，并把结果存储到FPGA外部的DDR存储器；

(10)循环步骤(4)到步骤(10)，计算不同时间采样点的箔条回波。

本发明与现有技术相比，具有如下有优点：

1)本发明由于基于FPGA平台进行模拟，因而资源消耗小，可移植性好。

2)本发明由于在计算过程中采用多种硬件计算单元，效率相比于电脑更高。

3)本发明从箔条的空间扩散模型出发，首先对箔条云在不同时刻的空间分布进行计算，然后根据空间分布计算每个时间采样点的箔条回波，能准确仿真空间中多发箔条云在不同时刻的回波。

附图说明

图1是本发明的实现流程图；

图2是本发明中第1000个时间采样点时的箔条在空间中的分布；

图3是本发明中在第100个时间采样点的箔条云回波信号；

图4是图3箔条云回波幅度的直方图；

图5是现有箔条云雷达回波仿真的箔条云回波幅度直方图。

具体实施方式

以下结合附图对本发明的实施例和效果做进一步详细描述。

参照图1，本发明的具体实施步骤如下:

步骤1.计算箔条的初始姿态、每个箔条的初始直角坐标和每个箔条的雷达散射截面值。

这三个计算过程互不依赖、相互独立，在FPGA内采用并行方式同时计算，具体实现如下：

1.1)计算箔条的初始姿态：在计算初始姿态参数时，需要生成均匀分布随机数和高斯分布随机数；本实例使用xorshiro128算法生成范围是(0,1)均匀分布随机数，模块名称为rand()，使用Box-Muller算法在FPGA内计算标准正态分布随机数z₀～N(0,1)，并根据公式z＝z₀σ+μ得到均值为μ、标准差为σ的高斯分布随机数z，模块名称为normrnd(μ,σ)；

调用随机数模块，计算各初始姿态的如下参数：

原点到第i个箔条的连线与正z轴之间的天顶角：

箔条在xOy平面的投影线与正x轴之间的方位角：

箔条转动角度：ζ(i)＝2π·rand(i)；

水平扩散速度：v(i)＝normrnd(0,5)；

把这些计算结果保存到RAM中；

1.2)计算每个箔条的初始直角坐标Lx₀(i)、Ly₀(i)、Lz₀(i)：

首先，调用rand()在球坐标中计算箔条云中各个箔条的随机初始位置(r(i),s₁(i),s₂(i))：

r(i)＝2rand(i)，s₁(i)＝2πrand(i)，s₂(i)＝2πrand(i)；

然后，根据预先设定的抛洒箔条处的直角坐标(hx,hy,hz)和箔条的随机初始位置(r(i),s₁(i),s₂(i))，计算每个箔条的初始直角坐标：

其中三角函数调用CORDIC IP核计算；

1.3)计算箔条云中各个箔条的雷达散射截面值：

首先，计算单根箔条平均雷达散射截面值：σ_{single_mean}＝0.153λ²，其中λ为入射波长；

其次，计算箔条云平均雷达散射截面值：σ_{cloud_mean}＝0.65Nσ_{single_mean}，其中N为单发箔条云中的箔条总数；

然后，调用对数IP核计算箔条云中的各个箔条雷达散射截面值：σ(i)＝-σ_{cloud_ mean}log(1-rand(i))，并将计算结果存入RAM中，供后续步骤调用。

步骤2.计算每条箔条所在高度z_t(i)处的大气密度ρ(i)。

2.1)计算每根箔条所在高度z_t(i)对应的位势高度H(i)：

其中，有效地球半径R′_m＝8500km，z_t(i)是每根箔条在时间采样点t时的高度；但t＝0时，由于箔条还未下落，所以z₀(i)＝Lz₀(i)；当t≥1时，由于箔条所在高度的大气密度会影响箔条的受力情况，进而改变箔条的瞬时速度，最终影响箔条在两个时间采样点中的位移增量，所以时间采样点t时的箔条位势高度采用t-1时的箔条高度z_t-1(i)计算得出；

2.2)根据每根箔条位势高度H(i)计算的大气密度ρ(i)：

ρ(i)＝ρ₀exp(aH(i))

其中，海面标准大气压ρ₀＝1.2250，拟合系数a＝-1.1553×10^-4，指数运算使用HLSDSP库中的exp函数实现。

步骤3.计算箔条云中每根箔条受力引起的三个方向上的位移增量dx₁(i)、dy₁(i)、dz₁(i)和箔条自身螺旋转动产生的水平位移增量dx₂(i)、dy₂(i)。

3.1)计算每根箔条受力引起的三个方向上的位移增量dx₁(i)、dy₁(i)、dz₁(i)：

3.1.1)根据上一采样时间点的箔条速度计算本采样时间点的箔条受力引起的三个方向上的瞬时速度vx_t(i)、vy_t(i)、vz_t(i)：

其中，t＝0时，vx₀(i)、vy₀(i)、vz₀(i)是箔条初速度；k_t-1(i)是采用上一个时间采样点中计算得到结果，首次计算时t＝1，k₀(i)使用箔条初速度和初始直角坐标计算得到；

3.1.2)按照如下步骤计算箔条在γ(i)平面上的投影面积S(i)：

γ(i)平面在空间中的位置按照如下向量确定：

计算在每个时间采样点处，每根箔条在γ(i)平面上的投影面积S(i)：

其中，d为箔条的直径，l为箔条的长度；上述公式在FPGA实现过程中，平方根计算使用浮点平方根功能核FSqrt实现，三角函数调用CORDIC IP核计算，计算过程中采用并行计算的优化方式，具体实现方式如下：

通过设立临时变量temp1＝dl，temp_x＝vx_t(i)-U_x，temp_y＝vy_t(i)-U_y，temp_z＝vz_t(i)-U_z，

来实现计算的并行，其中temp1、temp_x、temp_y、temp_z、temp_sc和temp_ss采用双端口RAM实现；

3.1.3)计算箔条受力引起的三个方向上的位移增量dx₁(i)、dy₁(i)、dz₁(i)：

把3.1.1)公式带入使用箔条速度增量计算箔条受力引起的三个方向上的如下位移增量的公式中：

得到计算箔条受力引起的三个方向上的最终位移增量的公式：

其中，U_x、U_y、U_z是大气风速在三个方向上的速度分量，ρ(i)是箔条所在高度的大气密度，g是重力加速度，m是箔条质量，dt是采样时间间隔；

对该公式的计算使用多阶流水线乘法器功能核Mulns进行优化。

3.2)计算箔条旋转引起的水平位移增量：

其中，dt是采样时间间隔，三角函数使用FPGA内部CORDIC IP核计算。

步骤4.计算本时间采样点的箔条直角坐标x_t(i)、y_t(i)、z_t(i)：

在初始状态t＝0时，计算初始箔条直角坐标x₀(i)＝Lx₀(i)，y₀(i)＝Ly₀(i)，z₀(i)＝Lz₀(i)；

在t≥1时，根据上一时间采样点的箔条直角坐标x_t-1(i)、y_t-1(i)、z_t-1(i)和步骤3.1)计算得到的箔条受力引起的三个方向上的位移增量，以及3.2)计算得到的箔条旋转引起的水平位移增量，计算本时间采样点的箔条直角坐标x_t(i)、y_t(i)、z_t(i)：

步骤5.箔条直角坐标x_t(i)、y_t(i)、z_t(i)，计算雷达到每根箔条的距离R_t(i)：

其中，(R_x0，R_y0，R_z0)为雷达所在位置的直角坐标，该公式中的平方运算采用FPGA内的DSP48乘法器单元实现，平方根运算使用浮点平方根功能核FSqrt实现。

步骤6.根据每根箔条到雷达的距离R_t(i)和箔条云中各个箔条的雷达散射截面值σ(i)，计算每个时间采样点的箔条的窄带信号回波s(t,i)：

其中，j是虚数单位，f_c是入射波频率，c代表光速。

本发明的效果可通过以下仿真进一步说明：

一.仿真参数，如表1、表2、表3所示：

表1系统仿真参数

系统仿真参数	参数值
		入射波频率	5GHz
发射功率	1w
		时间采样间隔	0.1s
仿真总时长	1200s
		箔条云数量	3
风速	(0,-5,0)m/s
		雷达位置	(0,11700,26000)

表2箔条参数

箔条参数	参数值
		单根箔条长度(米)	0.05
圆柱形箔条直径(米)	2.5×10-5
		单根箔条质量(千克)	1.45×10-8

表3箔条云运动参数

箔条云仿真参数	参数值
		单发箔条云箔条数量	1500
抛洒高度	12km
		转动速度起伏标准差	5
箔条云在xz平面投影与z轴夹角	π/3
		箔条云在xz平面投影与发射方向夹角	π/6
箔条初始速度	(0,0,0)m/s

二.仿真内容：

仿真1，根据表1中的时间采样间隔、箔条云数量、仿真总时长和风速，以及表2和表3中的参数，对箔条在空间中扩散的运动状态进行仿真，选取仿真过程中的第1000个时间采样点时的箔条在空间中的分布，在直角坐标系中绘图，结果如图2，从图2中可以看出，箔条云中的箔条在下落的过程中向四周扩散，呈现上方密、下部稀疏的状态。

仿真2，根据表1、表2、和表3中的参数进行箔条云回波仿真，选取仿真结果中的第100个时间采样点时的箔条云回波进行绘图，结果如图3；根据图3计算直方图，结果如图4。

从图4与图5所示的现有技术根据箔条云的雷达回波模型在电脑平台仿真的直方图进行比较，可见本发明与现有箔条云回波幅度的直方图结果相符，证明了本发明基于FPGA进行箔条云回波模拟的可行性和有效性。

Claims (10)

1.一种基于FPGA的箔条云雷达回波实时模拟方法，其特征在于，包括如下：

(1)计算箔条的初始姿态，包括原点到箔条的连线与正z轴之间的天顶角θ、箔条在xOy平面的投影线与正x轴之间的方位角

箔条转动角度ζ、水平扩散速度v，并将这些计算结果保存到RAM中；

(2)计算每个箔条的初始直角坐标Lx₀(i)、Ly₀(i)、Lz₀(i)，分别保存到RAM中，其中i为箔条的序号；

(3)计算每个箔条的雷达散射截面RCS值，把结果保存到RAM中；

(4)计算每条箔条所在高度z(i)处的大气密度ρ(i)；

(5)计算每条箔条受力引起的三个方向上的位移增量dx₁(i)、dy₁(i)、dz₁(i)；

(6)计算箔条螺旋转动引起的水平位移增量dx₂(i)、dy₂(i)；

(7)根据(5)得到的箔条位移增量和(6)得到的位移增量，计算本时间采样点的箔条直角坐标x_t(i)、y_t(i)、z_t(i)；

(8)根据雷达的空间坐标和每个箔条的空间坐标计算雷达到每个箔条的距离R_t(i)；

(9)根据步骤(3)得到的每根箔条的雷达散射截面RCS值和步骤(8)得到的雷达到箔条的距离R_t(i)，计算箔条的窄带信号回波S(t,i)，并把结果存储到FPGA外部的DDR存储器；

(10)循环步骤(4)到步骤(10)，计算不同时间采样点的箔条回波。

2.根据权利要求1所述的方法，其特征在于，所述(1)中计算箔条的初始状态，实现如下：

计算原点到第i个箔条的连线与正z轴之间的天顶角：

其中，rand()为生成取值(0,1)范围内的均匀分布随机数；

计算箔条在xOy平面的投影线与正x轴之间的方位角：

计算箔条转动角度：ζ(i)＝2π·rand(i)；

计算水平扩散速度：v(i)＝normrnd(0,5)，其中normrnd(0,5)为生成服从均值为0、标准差为5的正态分布的随机数。

3.根据权利要求1所述的方法，其特征在于，所述(2)中计算每个箔条的直角坐标Lx₀(i)、Ly₀(i)、Lz₀(i)，实现如下：

(2a)在球坐标中计算箔条云中各个箔条的随机初始位置(r(i),s₁(i),s₂(i))：

r(i)＝2rand(i)，s₁(i)＝2πrand(i)，s₂(i)＝2πrand(i)；

(2b)根据预先设定的抛洒箔条处的直角坐标(hx,hy,hz)和箔条的随机初始位置(r(i),s₁(i),s₂(i))

计算每个箔条的初始直角坐标：

4.根据权利要求1所述的方法，其特征在于，所述(3)根据统计方法计算箔条云中各个箔条雷达散射截面值，按照如下步骤进行：

(3a)计算单根箔条的平均雷达散射截面值：

σ_{single_mean}＝0.153λ²，其中λ为入射波长；

(3b)根据(3a)中的结果计算箔条云的平均雷达散射截面值：

σ_{cloud_mean}＝0.65Nσ_{single_mean}，其中N为单发箔条云中的箔条总数；

(3c)根据(3b)中的计算结果计算箔条云中的每根箔条的雷达散射截面值：

σ(i)＝-σ_{cloud_mean}log(1-rand(i))。

5.根据权利要求1所述的方法，其特征在于，所述(4)计算每条箔条所在高度z(i)处的大气密度ρ(i)，公式如下：

ρ(i)＝ρ₀exp(aH(i))，

其中，海面标准大气压ρ₀＝1.2250，拟合系数a＝-1.1553×10^-4，H(i)为与每根箔条所在高度z_t(i)对应的位势高度，

有效地球半径R′_m＝8500km。

6.根据权利要求1所述的方法，其特征在于，所述(5)中三个方向上的受力产生的位移增量,按照如下公式步骤计算：

其中，U_x、U_y、U_z是大气风速在三个方向上的速度分量，ρ(i)是箔条所在高度的大气密度，g是重力加速度，m是箔条质量，dt是采样时间间隔，vx_t(i)、vy_t(i)、vz_t(i)是箔条受力引起的三个方向上的瞬时速度，S(i)是箔条在γ(i)平面上的投影面积，

每根箔条的投影面积S(i)按照如下公式计算：

其中，d为箔条的直径，l为箔条的长度。

7.根据权利要求1所述的方法，其特征在于，所述(6)中由箔条旋转扩散引起的水平位移增量，按照如下公式计算：

其中，dt是时间采样间隔。

8.根据权利要求1所述的方法，其特征在于，所述(7)中计算本时间采样点的箔条直角坐标x_t(i)、y_t(i)、z_t(i)，实现如下：

在初始状态t＝0时，计算初始箔条直角坐标x₀(i)＝Lx₀(i)，y₀(i)＝Ly₀(i)，z₀(i)＝Lz₀(i)；

在t≥1时，根据上一时间采样点的箔条直角坐标x_t-1(i)、y_t-1(i)、z_t-1(i)和(5)、(6)中计算得到的位移增量，计算本时间采样点的箔条直角坐标x_t(i)、y_t(i)、z_t(i)：

9.根据权利要求1所述的方法，其特征在于，所述(8)中计算雷达到每个箔条的距离R_t(i)按下式进行：

其中，(R_x0，R_y0，R_z0)为雷达所在位置的直角坐标。

10.根据权利要求1所述的方法，其特征在于，所述(9)中计算箔条的窄带信号回波窄带信号回波S(t,i)，通过下式进行：

其中，f_c是入射波频率，c代表光速，σ(i)是箔条云中每根箔条的雷达散射截面值，R(i)是雷达与箔条云中各个箔条的距离。

Images