CN104076338B - 基于数字高程和数字地表覆盖的机载雷达杂波仿真方法 - Google Patents

Abstract

本发明公开了基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，涉及雷达技术领域，其步骤为：步骤1，构建雷达和载机在大地坐标系下的系统参数，以及雷达阵元组成的雷达天线阵面和雷达速度在东北天坐标系下的系统参数；步骤2，求取雷达与杂波单元C_k之间的视线矢量，并计算雷达与杂波单元C_k之间的距离R_k；步骤3，得到杂波单元C_k与雷达之间的入射余角；杂波单元C_k的雷达截面积S_k；杂波单元C_k单位面积的后向散射系数；杂波单元C_k的回波功率；步骤4，根据每一俯仰角判断杂波单元C_k是否被遮挡，得到遮挡标志δ_k；步骤5，构建杂波单元C_k的回波信号；将雷达照射范围内K个杂波单元的回波信号进行累积，得到杂波信号。本发明能够获取逼真的杂波仿真数据。

Description

基于数字高程和数字地表覆盖的机载雷达杂波仿真方法

技术领域

本发明属于雷达技术领域，具体的说是一种基于数字高程模型和数字地表覆盖模型的机载雷达杂波仿真方法，适用于机载相控阵雷达进行真实场景的杂波仿真。

技术背景

雷达是现代战争中不可缺少的装备，对于采用下视工作方式的机载预警雷达来说，地海面杂波对目标检测的影响十分突出，对杂波的抑制能力就成了检验雷达性能的重要指标。为了能够提供有效的杂波抑制方法，提升雷达在杂波中检测微弱信号的能力，必须对雷达工作环境的杂波特性有充分完备的认识。实测的杂波数据不可能在短时间内得到，而且花费极高。但是随着计算机技术的提高，利用计算机进行仿真杂波的方法来研究机载雷达的杂波特性，为雷达系统的设计和信号处理方法提供仿真数据，就显得极为重要了。

J.Ward提出了传统的杂波仿真方法，在极坐标系中将杂波按等距离环和方位角来划分成多个杂波单元，在某一方位角，随着距离的增大，外层的杂波单元面积增大，擦地角也随着距离的增大而减小。该传统仿真方法中假设每个杂波单元所对应的地表类型一致，即杂波服从独立同分布，然而在实际的环境中，尤其当面积逐渐增大之后，该杂波单元所对应的地表类型将可能不止一种，这就背离了杂波单元独立同分布的假设，对于描述远距离杂波单元的回波特性产生了比较大的误差。同时，传统仿真方法也假设地形没有起伏，在一定程度上反应不出真实场景(如高山、丘陵等)中地面起伏和擦地角变化，因此这就极大的限制了获取真实场景仿真的数据。

范国忠等人使用的数字高程模型进行杂波仿真的方法，能有效的描述真实场景中地面的起伏以及擦地角的变化，但是在描述地面起伏的遮挡判断算法计算量大，对于大范围的杂波仿真比较有限。王爱国等人仅仅利用了数字高程模型描述真实的地面起伏状态，对于真实的地貌特征，如城市、河流、沙漠等没有具体建模分析与仿真。

技术内容

针对目前基于独立同分布假设的杂波仿真所存在的缺陷，本发明提出一种基于数字高程模型和数字地表覆盖模型的机载相控阵雷达杂波仿真方法，结合数字高程模型和数字地表覆盖模型，改进了原有特定场景杂波仿真中遮挡判断计算量大的缺点并且加入了数字地表覆盖模型，对地表类型进行归并分类与建模，极大程度的还原了真实场景的各项特性，从而获得逼真的杂波仿真数据。

为达到上述目的，本发明采用以下技术方案予以实现。

一种基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，其特征在于，包括以下步骤：

步骤1，构建雷达在大地坐标系下的系统参数，以及雷达阵元组成的雷达天线阵面和雷达速度在东北天坐标系下的系统参数；

步骤2，将雷达天线阵面所处的东北天坐标系变换到地心坐标系，得到地心坐标系下的雷达阵元间隔矢量和地心坐标系下的雷达速度矢量将雷达所处的大地坐标系转换到地心坐标系，定义数字高程模型在大地坐标系下的杂波单元C_k，求取雷达与杂波单元C_k之间的视线矢量并计算雷达与杂波单元C_k之间的距离R_k；

步骤3，在大地坐标系下构建杂波单元相对于雷达的四边形入射平面，得到四边形的四个顶点坐标，左上顶点与右下顶点形成一条对角线矢量，右上顶点与左下顶点形成另一条对角线矢量，利用两条对角线矢量计算入射平面的法向矢量再计算雷达视线矢量利用入射平面的法向矢量和雷达视线矢量计算杂波单元C_k与雷达之间的入射余角根据入射平面的两条对角线矢量计算杂波单元C_k的雷达截面积S_k；

利用数字地表覆盖模型计算杂波单元C_k单位面积的后向散射系数利用杂波单元后向散射系数杂波单元雷达截面积S_k，杂波单元C_k与雷达之间的距离R_k，计算杂波单元C_k的回波功率；

步骤4，求取雷达与杂波单元C_k之间离散化点集；求取离散化点集中每一离散点相对于雷达的俯仰角；根据每一俯仰角判断杂波单元C_k是否被遮挡，得到遮挡标志δ_k；

步骤5，利用雷达与杂波单元C_k之间的视线矢量和雷达阵元间隔矢量构建回波空域信号导向矢量和回波时域信号导向矢量；利用回波空域信号导向矢量、回波时域信号导向矢量、回波功率、遮挡标志构建杂波单元C_k的回波信号；将雷达照射范围内K个杂波单元的回波信号进行累积，得到杂波信号。

上述技术方案的特点和进一步改进在于：

(1)步骤1中所构建的系统参数包括：

所述雷达为相控阵雷达；雷达与载机一体所处位置在大地坐标系下的坐标为P(l₀,b₀,h₀)，雷达天线阵面在东北天坐标系下以速度矢量飞行，雷达位于载机上，雷达为正侧视安装方式，雷达天线阵面为矩形平面，水平向N_a个阵元，俯仰向N_e个阵元，阵元间间隔矢量为天线主播束方位向与相控阵雷达天线的矩形阵面的法向垂直，俯仰向指向水平视线的无穷远处；其中E、N、V为东北天坐标轴，E轴指向正东，N轴指向正北，V与E、N构成右手坐标系，v_E为雷达在E轴下的速度分量，v_N为雷达在N轴下的速度分量，v_V为雷达V轴下的速度分量；L、B、H为大地坐标轴，L轴为经度轴，B为纬度轴，H为高度轴，l₀为载机经度，b₀为载机纬度，h₀为载机飞行高度。

(2)步骤2具体包括：

将雷达天线阵面所处的东北天坐标系变换到地心坐标系，将雷达所处的大地坐标系变换到地心坐标系，转换公式如下所示：

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{ccc}-\sin L& -\cos L\sin B& \cos L\cos B\\ \cos L& -\sin L\sin B& \sin L\cos B\\ 0& \text{cosB}& \sin B\end{array}\right]\×\left[\begin{array}{c}E\\ N\\ V\end{array}\right]---\left(1\right)$

其中X、Y、Z为地心坐标系，原点在地球中心，X轴指向本初子午线与赤道的交点，Z轴与地球的自转轴重合，指向北极，Y轴位于赤道平面与X轴垂直，形成一个右手坐标系；L代表大地坐标系的经度轴，B代表大地坐标系的纬度轴；

地心坐标系下，雷达阵元间隔矢量为：

$\stackrel{\→}{d}=d_x\stackrel{\→}{x}+d_y\stackrel{\→}{y}+d_z\stackrel{\→}{z}---\left(2\right)$

其中，分别代表地心坐标轴的3个单位矢量，d_x、d_y、d_z分别代表阵元间隔在地心坐标轴上的3个分量；

地心坐标系下，雷达的速度矢量为：

$\stackrel{\→}{v}=v_x\stackrel{\→}{x}+v_y\stackrel{\→}{y}+v_z\stackrel{\→}{z}---\left(3\right)$

其中，v_x、v_y、v_z分别代表速度在地心坐标轴上的3个分量；

大地坐标系下，数字高程模型中的每一条经度轴L和每一条纬度轴B相交形成矩形网格点，每一个网格点则对应一个高度，从数字高程模型中通过定位经纬度直接读取该经纬度所对应的高度，而每一个网格点就是一个杂波单元，设定杂波单元C_k在大地坐标系下坐标(l_k,b_k,h_k),k＝1,2,...,K，K为雷达照射范围内所有杂波单元的个数；

由大地坐标系转换到地心坐标系的转换公式如式(4)所示：

X＝(N+H)cosBcosL

Y＝(N+H)cosBsinL (4)

Z＝[N(1-e²)+H)]sinB

将雷达的大地坐标P(l₀,b₀,h₀)代入公式(4)计算雷达在地心坐标系下的坐标为(x₀,y₀,z₀)，将数字高程模型中的杂波单元C_k的大地坐标(l_k,b_k,h_k)代入公式(4)计算得到杂波单元在地心坐标系下的坐标为(x_k,y_k,z_k)；雷达与杂波单元C_k之间的视线矢量为：

${\stackrel{\→}{\mathrm{\α}}}_k=(x_k-x_0)\stackrel{\→}{x}+(y_k-y_0)\stackrel{\→}{y}+(z_k-z_0)\stackrel{\→}{z}---\left(5\right)$

计算杂波单元C_k与雷达之间的距离R_k：

$R_k=\sqrt{{(x_0-x_k)}^2+{(y_0-y_k)}^2+{(z_0-z_k)}^2}.---\left(6\right)$

(3)步骤3包括以下子步骤：

3a)在大地坐标系下，杂波单元C_k作为四边形入射平面的左上顶点与杂波单元C_k右相邻的杂波单元作为入射平面的右上顶点与杂波单元C_k下相邻的杂波单元作为入射平面的左下顶点与C_k右下相邻的杂波单元作为入射平面的右下顶点将雷达入射平面中的四个顶点对应的杂波单元按照公式(4)转换到地心坐标系下后，入射平面的左上顶点与右下顶点形成一条对角线矢量右上顶点与左下顶点形成另一条对角线矢量两条对角线矢量和叉乘得到入射平面的法向矢量杂波单元C_k与雷达P构成的矢量为雷达视线矢量；

3b)利用入射平面的法向矢量和雷达视线矢量计算杂波单元C_k与雷达之间的入射余角

其中，·表示点乘；

3c)根据入射平面的两条对角线矢量和计算杂波单元C_k的雷达截面积S_k：

$S_k=\frac{1}{2}\·\left|\stackrel{\→}{C_k^{\mathrm{ru}}{C}_k^{\mathrm{ld}}}\right|\·\left|\stackrel{\→}{C_k^{\mathrm{lu}}{C}_k^{\mathrm{rd}}}\right|\·\sin \left({\cos }^{-1}\left(\frac{\stackrel{\→}{C_k^{\mathrm{ru}}{C}_k^{\mathrm{ld}}}\·\stackrel{\→}{C_k^{\mathrm{lu}}{C}_k^{\mathrm{rd}}}}{\left|\stackrel{\→}{C_k^{\mathrm{ru}}{C}_k^{\mathrm{ld}}}\right|\·\left|\stackrel{\→}{C_k^{\mathrm{lu}}{C}_k^{\mathrm{rd}}}\right|}\right)\right)---\left(8\right)$

3d)利用数字地表覆盖模型计算杂波单元C_k单位面积的后向散射系数

3e)利用杂波单元后向散射系数杂波单元雷达截面积S_k，杂波单元C_k与雷达之间的距离R_k，计算杂波单元C_k的回波功率：

${\mathrm{\ξ}}_k=\frac{P_t{G}_t{\mathrm{\λ}}^2{\mathrm{\σ}}_k^0{S}_k}{{\left(4\mathrm{\π}\right)}^2{L}_s{R}_k^4}---\left(10\right)$

式(10)中，P_t为雷达发射峰值功率，G_t为杂波单元C_k的方向图增益，λ为载频波长，为杂波单元后向散射系数，S_k为杂波单元雷达截面积，R_k为杂波单元C_k与雷达之间的距离，L_s为雷达损耗。

(4)子步骤3d)具体包括：

根据杂波单元C_k的大地经纬度坐标l_k,b_k从数字地表覆盖模型中获取杂波单元C_k对应的地表标签；根据该地表标签计算杂波单元C_k的单位面积散射系数通过下式(9-a)和(9-b)表示：

非水体 (9-a)

水体 (9-b)

在非水体的式(9-a)中：为入射余角，

θ_c＝sin^-1(λ/4πh_e)，f₀为雷达工作频率，单位GHz，λ为雷达工作波长，W≈1；A为幅度参数，B为相位参数，β₀为镜面反射参数，为角度参数；

在水体的式(9-b)中：其中，为入射余角，SS是水情级数，

θ_c＝sin^-1(λ/4πh_e)，h_e＝0.025+0.046SS^1.72，β＝[2.44(SS+1)^1.08]/57.29为镜面反射参数，W＝1.9。

(5)步骤4包括以下子步骤：

4a)在大地坐标系下，雷达位置坐标为P(l₀,b₀,h₀)，当杂波单元C_k在第一卦限时，则经度轴L的坐标大于纬度轴B的坐标，连接雷达位置P与杂波单元C_k对应的点，则与经度轴L和纬度轴B相交于J个点，根据公式(11)计算雷达P到杂波单元C_k投影的离散化点集{l_j,b_j,h_j},j＝0,1,...,J：

当杂波单元C_k在第二卦限时，则纬度轴B的坐标大于经度轴L的坐标，连接雷达位置P与杂波单元C_k，则与经度轴L和纬度轴B相交于J个点，根据公式(12)计算雷达P到杂波单元C_k投影的离散化点集{l_j,b_j,h_j},j＝0,1,...,J：

其中表示向下取整，j＝0时即为雷达坐标，j＝J即为杂波单元C_k的坐标；

当杂波单元C_k在其他三至八任一卦限中时，在经度轴L的坐标大于纬度轴B的坐标情况下，用式(11)计算离散化点集，在纬度轴B的坐标大于经度轴L的坐标情况下，用式(12)计算离散化点集；

4b)通过公式(4)先将雷达与离散化点集中的离散点转换到地心坐标系下，然后通过公式(13)计算得到雷达与离散化点集中的离散点之间的距离R_j：

$R_j=\sqrt{{(x_0-x_j)}^2+{(y_0-y_j)}^2+{(z_0-z_j)}^2}---\left(13\right)$

离散化点集中每一离散点相对于雷达的俯仰θ_j表示为下式，j＝0,1,...,J，J表示离散化点集中离散点的数目；

${\mathrm{\θ}}_j={\sin }^{-1}\left(\frac{h_0-h_j}{R_j}\right),j=\mathrm{0,1},...,J---\left(14\right)$

4c)将离散化点集中的每一个俯仰角与杂波单元C_k的俯仰角进行比较，只要在离散化点集中存在一个俯仰角小于杂波单元C_k的俯仰角，则杂波单元C_k被遮挡，否则就不被遮挡，遮挡标志为δ_k：

(6)步骤5包括以下子步骤：

回波空域信号导向矢量为

$s_{\mathrm{sk}}\left(f_{\mathrm{sk}}\right)={[1,e^{j2\mathrm{\π}\·f_{\mathrm{sk}}},...,e^{j2\mathrm{\π}\·(N_a{N}_e-1)f_{\mathrm{sk}}}]}^T---\left(16\right)$

回波时域信号导向矢量为

$s_{\mathrm{tk}}\left(f_t\right)={[1,e^{j2\mathrm{\π}\·f_{\mathrm{tk}}},...,e^{j2\mathrm{\π}\·(M-1)f_{\mathrm{tk}}}]}^T---\left(17\right)$

其中，为空间频率，为归一化多普勒频率，N_a为水平向阵元个数，N_e为俯仰向阵元个数，M为一个脉冲重复间隔内发射脉冲的个数，λ为载频波长，f_r为脉冲重复频率，为雷达与杂波单元C_k之间的视线矢量，为雷达阵元间隔矢量，为雷达速度矢量；

杂波单元C_k的回波信号为：

$c_k={\mathrm{\δ}}_k\sqrt{{\mathrm{\ξ}}_k}{s}_{\mathrm{sk}}\⊗s_{\mathrm{tk}}---\left(18\right)$

其中，δ_k为遮挡标志，ξ_k为回波功率，表示Kronecker积；

将雷达照射范围内K个杂波单元的回波信号进行累积，得到杂波信号为

$c=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{c}_k.---\left(19\right)$

本发明主要针对上述已有方法的缺点进行改进与提升，包括利用数字高程模型描述真实场景中地面起伏的遮挡判断算法的优化，以及利用数字地表覆盖模型描述真实场景中的地貌特征，获得较逼真的杂波仿真数据。

附图说明

下面结合附图和具体实施方式对本发明做进一步说明。

图1为本发明的流程图；

图2为离散点集遮挡判断示意图；横坐标表示离散点集，纵坐标表示高度；

图3为数字高程模型遮挡判断示意图；横坐标表示经度轴，纵坐标表示纬度轴；

图4为入射平面示意图；横坐标表示经度轴，纵坐标表示纬度轴；

图5为遮挡判断后的数字高程模型；横坐标表示经度轴，纵坐标表示纬度轴；

图6为数字地表覆盖图；横坐标表示经度轴，纵坐标表示纬度轴；

图7为本发明仿真方法的单通道杂波距离-多普勒图；横坐标表示多普勒，纵坐标表示距离门；

图8为Ward方法的单通道距离-多普勒图；横坐标表示多普勒，纵坐标表示距离门；

图9为实测MACARM数据的单通道距离-多普勒图；横坐标表示多普勒，纵坐标表示距离门。

具体实施方式

参照图1，说明本发明的一种基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，用于机载雷达杂波仿真，包括以下步骤：

步骤1，构建雷达和载机在大地坐标系下的系统参数，以及雷达阵元组成的雷达天线阵面和雷达速度在东北天坐标系下的系统参数。

雷达与载机一体所处位置在大地坐标系下的坐标为P(l₀,b₀,h₀)，雷达天线阵面在东北天坐标系下以速度矢量飞行，载机上的相控阵雷达为正侧视安装方式，雷达天线阵面为矩形平面，水平向N_a个阵元，俯仰向N_e个阵元，阵元间间隔矢量为天线主播束方位向与相控阵雷达天线的矩形阵面的法向垂直，俯仰向指向水平视线的无穷远处。其中E、N、V为东北天坐标轴，E轴指向正东，N轴指向正北，V与E、N构成右手坐标系，v_E为雷达在E轴下的速度分量，v_N为雷达在N轴下的速度分量，v_V为雷达V轴下的速度分量；L、B、H为大地坐标轴，L轴为经度轴，B为纬度轴，H为高度轴，l₀为载机经度，b₀为载机纬度，h₀为载机飞行高度。

步骤2，将雷达阵面所处的东北天坐标系变换到地心坐标系，得到地心坐标系下的雷达阵元间隔矢量和地心坐标系下的雷达的速度矢量将雷达所处的大地坐标系转换到地心坐标系，定义数字高程模型在大地坐标系下的杂波单元C_k，求取雷达与杂波单元C_k之间的视线矢量并计算雷达与杂波单元C_k之间的距离R_k。

雷达处于东北天坐标系，将雷达天线阵面以及雷达所处的东北天坐标系变换到地心坐标系，转换公式如下所示：

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{ccc}-\sin L& -\cos L\sin B& \cos L\cos B\\ \cos L& -\sin L\sin B& \sin L\cos B\\ 0& \text{cosB}& \sin B\end{array}\right]\×\left[\begin{array}{c}E\\ N\\ V\end{array}\right]---\left(1\right)$

其中X、Y、Z为地心坐标系，原点在地球中心，X轴指向本初子午线与赤道的交点，Z轴与地球的自转轴重合，指向北极，Y轴位于赤道平面与X轴垂直，形成一个右手坐标系；L代表大地坐标系的经度轴，B代表大地坐标系的纬度轴。

地心坐标系下，雷达阵元间隔矢量为：

$\stackrel{\→}{d}=d_x\stackrel{\→}{x}+d_y\stackrel{\→}{y}+d_z\stackrel{\→}{z}---\left(2\right)$

其中，分别代表地心坐标轴的3个单位矢量，d_x、d_y、d_z分别代表阵元间隔在地心坐标轴上的3个分量。

地心坐标系下，雷达的速度矢量为：

$\stackrel{\→}{v}=v_x\stackrel{\→}{x}+v_y\stackrel{\→}{y}+v_z\stackrel{\→}{z}---\left(3\right)$

其中，v_x、v_y、v_z分别代表速度在地心坐标轴上的3个分量。

大地坐标系下，数字高程模型中的每一条经度轴L和每一条纬度轴B相交形成矩形网格点，每一个网格点则对应一个高度，从数字高程模型中通过定位经纬度直接读取该经纬度所对应的高度，而每一个网格点就是一个杂波单元，设定杂波单元C_k在大地坐标系下坐标(l_k,b_k,h_k),k＝1,2,...,K，K为雷达照射范围内所有杂波单元的个数。

由大地坐标系转换到地心坐标系的转换公式如式(4)所示：

X＝(N+H)cosBcosL

Y＝(N+H)cosBsinL (4)

Z＝[N(1-e²)+H)]sinB

将雷达的大地坐标P(l₀,b₀,h₀)代入公式(4)计算雷达在地心坐标系下的坐标为(x₀,y₀,z₀)，将数字高程模型中的杂波单元C_k的大地坐标(l_k,b_k,h_k)代入公式(4)计算得到杂波单元在地心坐标系下的坐标为(x_k,y_k,z_k)。那么雷达与杂波单元C_k之间的视线矢量为：

${\stackrel{\→}{\mathrm{\α}}}_k=(x_k-x_0)\stackrel{\→}{x}+(y_k-y_0)\stackrel{\→}{y}+(z_k-z_0)\stackrel{\→}{z}---\left(5\right)$

计算杂波单元C_k与雷达之间的距离R_k：

$R_k=\sqrt{{(x_0-x_k)}^2+{(y_0-y_k)}^2+{(z_0-z_k)}^2}---\left(6\right)$

步骤3，构建杂波单元C_k相对于雷达的四边形入射平面，得到入射平面的法向矢量雷达视线矢量和入射平面的两条对角线矢量；利用入射平面的法向矢量和雷达视线矢量计算杂波单元C_k与雷达之间的入射余角根据入射平面的两条对角线矢量计算杂波单元C_k的雷达截面积S_k；利用数字地表覆盖模型计算杂波单元C_k单位面积的后向散射系数利用杂波单元后向散射系数杂波单元雷达截面积S_k，杂波单元C_k与雷达之间的距离R_k，计算杂波单元C_k的回波功率。

3a)在大地坐标系下，构建如图4所示的四边形入射平面的几何模型，其中杂波单元C_k作为入射平面的左上顶点与杂波单元C_k右相邻的杂波单元作为入射平面的右上顶点与杂波单元C_k下相邻的杂波单元作为入射平面的左下顶点与C_k右下相邻的杂波单元作为入射平面的右下顶点将雷达平面中的四个杂波单元转换到地心坐标系下后，入射平面的左上顶点与右下顶点形成一条对角线矢量右上顶点与左下顶点形成另一条对角线矢量两条对角线矢量和叉乘得到入射平面的法向矢量杂波单元C_k与雷达位置P构成的矢量为雷达视线矢量。与构成入射角。

3b)利用入射平面的法向矢量和雷达视线矢量计算杂波单元C_k与雷达之间的入射余角

其中，·表示点乘；

3c)根据入射平面的两条对角线矢量和计算杂波单元C_k的雷达截面积S_k：

$S_k=\frac{1}{2}\·\left|\stackrel{\→}{C_k^{\mathrm{ru}}{C}_k^{\mathrm{ld}}}\right|\·\left|\stackrel{\→}{C_k^{\mathrm{lu}}{C}_k^{\mathrm{rd}}}\right|\·\sin \left({\cos }^{-1}\left(\frac{\stackrel{\→}{C_k^{\mathrm{ru}}{C}_k^{\mathrm{ld}}}\·\stackrel{\→}{C_k^{\mathrm{lu}}{C}_k^{\mathrm{rd}}}}{\left|\stackrel{\→}{C_k^{\mathrm{ru}}{C}_k^{\mathrm{ld}}}\right|\·\left|\stackrel{\→}{C_k^{\mathrm{lu}}{C}_k^{\mathrm{rd}}}\right|}\right)\right)---\left(8\right)$

3d)利用数字地表覆盖模型计算杂波单元C_k单位面积的后向散射系数

本发明所使用的数字地表覆盖模型为欧洲宇航局产品，由联合国地表分类系统将地表覆盖类型细分为的22种，即有22种地表标签。在大地坐标系下，与数字高程模型类似，每一个网格点(杂波单元)对应一个经纬度坐标以及一种地表标签，根据杂波单元C_k的经纬度坐标l_k,b_k从数字地表覆盖模型中获取杂波单元C_k对应的地表标签，本发明中将这些地表标签粗分为5种地表覆盖类型，确定杂波单元C_k的地表标签所属的地表覆盖类型，并与所建立后向散射系数模型一一对应。分类后的地表标签、地表覆盖类型与后向散射系数模型对应的地类及参数如表1所示。

根据杂波单元C_k的经纬度坐标l_k,b_k从数字地表覆盖模型中获取杂波单元C_k对应的地表标签；根据该地表标签计算杂波单元C_k的单位面积散射系数为：

非水体 (9-a)

水体 (9-b)

在非水体的式(9-a)中：为入射余角，

θ_c＝sin^-1(λ/4πh_e)， f₀为雷达工作频率，单位GHz，λ为雷达工作波长，W≈1。A为幅度参数，B为相位参数，β₀为镜面反射参数，为角度参数，参数值见表1。

在水体的式(9-b)中：其中，为入射余角，SS是水情级数，

θ_c＝sin^-1(λ/4πh_e)，h_e＝0.025+0.046SS^1.72，β＝[2.44(SS+1)^1.08]/57.29为镜面反射参数，W＝1.9。

表1

3e)利用杂波单元后向散射系数杂波单元雷达截面积S_k，杂波单元C_k与雷达之间的距离R_k，计算杂波单元C_k的回波功率：

${\mathrm{\ξ}}_k=\frac{P_t{G}_t{\mathrm{\λ}}^2{\mathrm{\σ}}_k^0{S}_k}{{\left(4\mathrm{\π}\right)}^2{L}_s{R}_k^4}---\left(10\right)$

式(10)中，P_t为雷达发射峰值功率，G_t为杂波单元C_k的方向图增益，λ为载频波长，为杂波单元后向散射系数，S_k为杂波单元雷达截面积，R_k为杂波单元C_k与雷达之间的距离，L_s为雷达损耗。

步骤4，求取雷达与杂波单元C_k之间离散化点集；求取离散化点集中每一离散点相对于雷达的俯仰角；根据每一俯仰角判断杂波单元C_k是否被遮挡，得到遮挡标志δ_k。

4a)在大地坐标系下，雷达位置坐标为P(l₀,b₀,h₀)，若杂波单元C_k在第一卦限，则经度轴L的坐标大于纬度轴B的坐标，连接雷达位置P与杂波单元C_k对应的点，则与经度轴L和纬度轴B相交于J个点，根据公式(11)计算雷达P到杂波单元C_k投影的离散化点集{l_j,b_j,h_j},j＝0,1,...,J：

若杂波单元C_k在第二卦限，则纬度轴B的坐标大于经度轴L的坐标，连接雷达位置P与杂波单元C_k，则与经度轴L和纬度轴B相交于J个点，根据公式(12)计算雷达P到杂波单元C_k投影的离散化点集{l_j,b_j,h_j},j＝0,1,...,J：

其中表示向下取整，j＝0时即为雷达坐标，j＝J即为杂波单元C_k的坐标。

若杂波单元C_k在其他三至八任一卦限中时，在经度轴L的坐标大于纬度轴B的坐标情况下，就用式(11)计算离散化点集，在纬度轴B的坐标大于经度轴L的坐标情况下，就用式(12)计算离散化点集。

求得雷达到杂波单元C_k投影的离散化点集后，判断杂波单元C_k是否被遮挡。

4b)通过公式(4)先将雷达与离散化点集中的离散点C_j转换到地心坐标系下，然后通过公式(13)计算得到雷达与离散化点集中的离散点C_j之间的距离R_j：

$R_j=\sqrt{{(x_0-x_j)}^2+{(y_0-y_j)}^2+{(z_0-z_j)}^2}---\left(13\right)$

如图2所示，再计算离散化点集中每一离散点相对于雷达的俯仰θ_j,j＝0,1,...,J，J表示离散化点集中杂波单元的数目。

${\mathrm{\θ}}_j={\sin }^{-1}\left(\frac{h_0-h_j}{R_j}\right),j=\mathrm{0,1},...,J---\left(14\right)$

4c)将离散化点集中的每一个俯仰角与杂波单元C_k的俯仰角进行比较，只要在离散化点集中存在一个俯仰角小于杂波单元C_k的俯仰角，则杂波单元C_k被遮挡，否则就不被遮挡，那么遮挡标志为δ_k。

在图2中，杂波单元C_k的仰角为θ_J，可以发现θ_J-2<θ_J，那么该杂波单元被遮挡，遮挡标志δ_k＝0。

步骤5，利用雷达与杂波单元C_k之间的视线矢量和雷达阵元间隔矢量构建回波空域信号导向矢量和回波时域信号导向矢量；利用回波空域信号导向矢量、回波时域信号导向矢量、回波功率、遮挡标志构建杂波单元C_k的回波信号；将雷达照射范围内K个杂波单元的回波信号进行累积，得到杂波信号。

回波空域信号导向矢量为

$s_{\mathrm{sk}}\left(f_{\mathrm{sk}}\right)={[1,e^{j2\mathrm{\&pi;}\&CenterDot;f_{\mathrm{sk}}},...,e^{j2\mathrm{\&pi;}\&CenterDot;(N_a{N}_e-1)f_{\mathrm{sk}}}]}^T---\left(16\right)$

回波时域信号导向矢量为

$s_{\mathrm{tk}}\left(f_t\right)={[1,e^{j2\mathrm{\&pi;}\&CenterDot;f_{\mathrm{tk}}},...,e^{j2\mathrm{\&pi;}\&CenterDot;(M-1)f_{\mathrm{tk}}}]}^T---\left(17\right)$

其中为空间频率，为归一化多普勒频率，N_a为水平向阵元个数，N_e为俯仰向阵元个数，M为一个脉冲重复间隔内发射脉冲的个数，λ为载频波长，f_r为脉冲重复频率，为雷达与杂波单元C_k之间的视线矢量，为雷达阵元间隔矢量，为雷达速度矢量。

本发明中，回波空域信号导向矢量和回波时域信号导向矢量就是杂波单元C_k的回波信号相位。

杂波单元C_k的回波信号为

$c_k={\mathrm{\&delta;}}_k\sqrt{{\mathrm{\&xi;}}_k}{s}_{\mathrm{sk}}\&CircleTimes;s_{\mathrm{tk}}---\left(18\right)$

其中，δ_k为遮挡标志，ξ_k为回波功率，表示Kronecker积。

将雷达照射范围内K个杂波单元的回波信号进行累积，得到杂波信号为

$c=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{c}_k---\left(19\right)$

下面结合仿真实验对本发明的效果做进一步说明。

1、仿真参数

在本实验中，载机和雷达在大地坐标系下的坐标为(l₀,b₀,h₀)＝(-76.7167°,38.9645°,3589.8m)，载频为1.24GHz，相控阵方位向和俯仰向阵元数为M＝11,N＝2，阵元间隔分别为0.1407m和0.1092m。脉冲个数为P＝128，脉冲重复频率f_r＝1984Hz，带宽B＝0.8MHz，采样频率为f_s＝1.25MHz，脉宽为50.4us，雷达峰值功率为1500W，雷达损耗为L_s＝10dB，飞机以(v_E,v_N,v_V)＝(120.7,-37.4,2.4)m/s的速度飞行，主波束方位向指向与飞机轴向垂直，俯仰向指向无穷远处。上述所述参数基本与MCARM(Multi-Channel AirborneRadar Measurements)数据的参数基本一致。

2、仿真数据处理结果及分析

A.为了说明本发明的优越性，首先根据上述雷达系统参数使用Ward方法进行杂波仿真，仿真的场景没有地形的起伏并且为单一地貌。然后采用本发明方法仿真，该方法采用数字高程模型能体现真实的地形起伏，并结合数字地表覆盖模型能体现真实的地表类型特征。

图5为本发明仿真方法得到的高程遮挡判断图，图右边的矩形条代表高度，颜色越亮处高度越高，颜色最深处代表该处的杂波单元被遮挡的或者在雷达照射范围之外，从图中可以看出遮挡效果正确。图6为本发明仿真方法得到的数字地表覆盖图，图右边的矩形条代表地表标签，能够真实体现雷达所处位置周围的地表类型。

图7为本发明仿真方法仿真得到的杂波数据距离多普勒图，图8为Ward方法仿真得到的距离多普勒图，对比图7和图8可以看出Ward方法仿真的杂波距离多普勒图较均匀而不能体现出雷达照射环境的地表类型的变化。而本发明方法则明显可以体现出地形的起伏以及地表类型的变化，如图7中，距离门200至500，多普勒80至120区域为深色区域，这片区域表示的则是水体。

B.为了进一步说明本发明的优势，图9为MCARM实测数据获得的距离多普勒图。对比图8和图9，Ward仿真方法得到的距离多普勒图较均匀，没有特别深或者特别浅的地方，而MCARM实测数据中距离门300至500，多普勒80至120区域颜色较深，代表地表类型水体，故Ward仿真方法完全不能体现实测数据的杂波特性。对比图7和图9，本发明方法仿真得到的距离多普勒图与MCARM实测数据地距离多普勒图在大致形状上相似度很高，图9中距离门300至500，多普勒80至120区域颜色较深，代表地表类型水体，图7中，距离门300至500，多普勒80至120区域颜色较深，也代表地表类型水体。故发明方法能体现实测数据的杂波特性，相较于Ward仿真方法，本发明方法在体现地表类型特征上具有较大的优势。

通过上述分析可以得出结论：本发明所述的杂波仿真方法较Ward杂波仿真方法有较大优势，可以使用本方法获得较逼真的杂波仿真数据。

Claims (7)

1.一种基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，其特征在于，包括以下步骤：

步骤1，构建雷达在大地坐标系下的系统参数，以及雷达阵元组成的雷达天线阵面和雷达速度在东北天坐标系下的系统参数；

步骤2，将雷达天线阵面所处的东北天坐标系变换到地心坐标系，得到地心坐标系下的雷达阵元间隔矢量和地心坐标系下的雷达速度矢量将雷达所处的大地坐标系转换到地心坐标系，定义数字高程模型在大地坐标系下的杂波单元C_k，求取雷达与杂波单元C_k之间的视线矢量并计算雷达与杂波单元C_k之间的距离R_k；

步骤3，在大地坐标系下构建杂波单元相对于雷达的四边形入射平面，得到四边形的四个顶点坐标，左上顶点与右下顶点形成一条对角线矢量，右上顶点与左下顶点形成另一条对角线矢量，利用两条对角线矢量计算入射平面的法向矢量再计算雷达视线矢量利用入射平面的法向矢量和雷达视线矢量计算杂波单元C_k与雷达之间的入射余角根据入射平面的两条对角线矢量计算杂波单元C_k的雷达截面积S_k；

利用数字地表覆盖模型计算杂波单元C_k单位面积的后向散射系数利用杂波单元后向散射系数杂波单元雷达截面积S_k，杂波单元C_k与雷达之间的距离R_k，计算杂波单元C_k的回波功率；

步骤4，求取雷达与杂波单元C_k之间离散化点集；求取离散化点集中每一离散点相对于雷达的俯仰角；根据每一俯仰角判断杂波单元C_k是否被遮挡，得到遮挡标志δ_k；

步骤5，利用雷达与杂波单元C_k之间的视线矢量和雷达阵元间隔矢量构建回波空域信号导向矢量和回波时域信号导向矢量；利用回波空域信号导向矢量、回波时域信号导向矢量、回波功率、遮挡标志构建杂波单元C_k的回波信号；将雷达照射范围内K个杂波单元的回波信号进行累积，得到杂波信号。

2.根据权利要求1所述的一种基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，其特征在于，步骤1中所构建的系统参数包括：

所述雷达为相控阵雷达；雷达与载机一体所处位置在大地坐标系下的坐标为P(l₀，b₀，h₀)，雷达天线阵面在东北天坐标系下以速度矢量飞行，雷达位于载机上，雷达为正侧视安装方式，雷达天线阵面为矩形平面，水平向N_a个阵元，俯仰向N_e个阵元，阵元间隔矢量为天线主波束方位向与相控阵雷达天线的矩形阵面的法向垂直，俯仰向指向水平视线的无穷远处；其中E、N、V为东北天坐标轴，E轴指向正东，N轴指向正北，V与E、N构成右手坐标系，v_E为雷达在E轴下的速度分量，v_N为雷达在N轴下的速度分量，v_V为雷达V轴下的速度分量；L、B、H为大地坐标轴，L轴为经度轴，B为纬度轴，H为高度轴，l₀为载机经度，b₀为载机纬度，h₀为载机飞行高度。

3.根据权利要求1所述的一种基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，其特征在于，步骤2具体包括：

将雷达天线阵面所处的东北天坐标系变换到地心坐标系，转换公式如下所示：

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{ccc}-\sin L& -\cos L\sin B& \cos L\cos B\\ \cos L& -\sin L\sin B& \sin L\cos B\\ 0& \cos B& \sin B\end{array}\right]\&times;\left[\begin{array}{c}E\\ N\\ V\end{array}\right]---\left(1\right)$

其中X、Y、Z为地心坐标系，原点在地球中心，X轴指向本初子午线与赤道的交点，Z轴与地球的自转轴重合，指向北极，Y轴位于赤道平面与X轴垂直，形成一个右手坐标系；L代表大地坐标系的经度轴，B代表大地坐标系的纬度轴；

地心坐标系下，雷达阵元间隔矢量为：

$\stackrel{\&OverBar;}{d}=d_x\stackrel{\&RightArrow;}{x}+d_y\stackrel{\&RightArrow;}{y}+d_z\stackrel{\&RightArrow;}{z}---\left(2\right)$

其中，分别代表地心坐标轴的3个单位矢量，d_x、d_y、d_z分别代表阵元间隔在地心坐标轴上的3个分量；

地心坐标系下，雷达的速度矢量为：

$\stackrel{\&RightArrow;}{v}=v_x\stackrel{\&RightArrow;}{x}+v_y\stackrel{\&RightArrow;}{y}+v_z\stackrel{\&RightArrow;}{z}---\left(3\right)$

其中，v_x、v_y、v_z分别代表速度在地心坐标轴上的3个分量；

大地坐标系下，数字高程模型中的每一条经度轴L和每一条纬度轴B相交形成矩形网格点，每一个网格点则对应一个高度，从数字高程模型中通过定位经纬度直接读取该经纬度所对应的高度，而每一个网格点就是一个杂波单元，设定杂波单元C_k在大地坐标系下坐标(l_k，b_k，h_k)，k＝1，2，...，K，K为雷达照射范围内所有杂波单元的个数；

将雷达所处的大地坐标系变换到地心坐标系，由大地坐标系转换到地心坐标系的转换公式如式(4)所示：

$\begin{array}{c}X=(N+H)\cos B\cos L\\ Y=(N+H)\cos B\sin L\\ Z=\&lsqb;N(1-e^2)+H)\&rsqb;\sin B\end{array}---\left(4\right)$

将雷达的大地坐标P(l₀，b₀，h₀)代入公式(4)计算雷达在地心坐标系下的坐标为(x₀，y₀，z₀)，将数字高程模型中的杂波单元C_k的大地坐标(l_k，b_k，h_k)代入公式(4)计算得到杂波单元在地心坐标系下的坐标为(x_k，y_k，z_k)；雷达与杂波单元C_k之间的视线矢量为：

${\stackrel{\&RightArrow;}{\mathrm{\&alpha;}}}_k=(x_k-x_0)\stackrel{\&RightArrow;}{x}+(y_k-y_0)\stackrel{\&RightArrow;}{y}+(z_k-z_0)\stackrel{\&RightArrow;}{z}---\left(5\right)$

计算杂波单元C_k与雷达之间的距离R_k：

$R_k=\sqrt{{(x_0-x_k)}^2+{(y_0-y_k)}^2+{(z_0-z_k)}^2}.---\left(6\right)$

4.根据权利要求3所述的一种基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，其特征在于，步骤3包括以下子步骤：

3a)在大地坐标系下，杂波单元C_k作为四边形入射平面的左上顶点与杂波单元C_k右相邻的杂波单元作为入射平面的右上顶点与杂波单元C_k下相邻的杂波单元作为入射平面的左下顶点与C_k右下相邻的杂波单元作为入射平面的右下顶点将雷达入射平面中的四个顶点对应的杂波单元按照公式(4)转换到地心坐标系下后，入射平面的左上顶点与右下顶点形成一条对角线矢量右上顶点与左下顶点形成另一条对角线矢量两条对角线矢量和叉乘得到入射平面的法向矢量杂波单元C_k与雷达P构成的矢量为雷达视线矢量；

3b)利用入射平面的法向矢量和雷达视线矢量计算杂波单元C_k与雷达之间的入射余角

其中，·表示点乘；

3c)根据入射平面的两条对角线矢量和计算杂波单元C_k的雷达截面积S_k：

$S_k=\frac{1}{2}\&CenterDot;\left|\stackrel{\&RightArrow;}{C_k^{\mathrm{ru}}{C}_k^{\mathrm{ld}}}\right|\&CenterDot;\left|\stackrel{\&RightArrow;}{C_k^{\mathrm{lu}}{C}_k^{\mathrm{rd}}}\right|\&CenterDot;\sin \left({\cos }^{-1}\left(\frac{\stackrel{\&RightArrow;}{C_k^{\mathrm{ru}}{C}_k^{\mathrm{ld}}}\&CenterDot;\stackrel{\&RightArrow;}{C_k^{\mathrm{lu}}{C}_k^{\mathrm{rd}}}}{\left|\stackrel{\&RightArrow;}{C_k^{\mathrm{ru}}{C}_k^{\mathrm{ld}}}\right|\&CenterDot;\stackrel{\&RightArrow;}{C_k^{\mathrm{lu}}{C}_k^{\mathrm{rd}}}|}\right)\right)---\left(8\right)$

3d)利用数字地表覆盖模型计算杂波单元C_k单位面积的后向散射系数

3e)利用杂波单元后向散射系数杂波单元雷达截面积S_k，杂波单元C_k与雷达之间的距离R_k，计算杂波单元C_k的回波功率：

${\mathrm{\&xi;}}_k=\frac{P_t{G}_t{\mathrm{\&lambda;}}^2{\mathrm{\&sigma;}}_k^0{S}_k}{{\left(4\mathrm{\&pi;}\right)}^3{L}_s{R}_k^4}---\left(10\right)$

式(10)中，P_t为雷达发射峰值功率，G_t为杂波单元C_k的方向图增益，λ为载频波长，为杂波单元后向散射系数，S_k为杂波单元雷达截面积，R_k为杂波单元C_k与雷达之间的距离，L_s为雷达损耗。

5.根据权利要求4所述的一种基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，其特征在于，子步骤3d)具体包括：

根据杂波单元C_k的大地经纬度坐标l_k，b_k从数字地表覆盖模型中获取杂波单元C_k对应的地表标签；根据该地表标签计算杂波单元C_k的后向散射系数通过下式(9-a)和(9-b)表示：

在非水体的式(9-a)中：为入射余角，

θ_c＝sin^-1(λ/4πh_e)，f₀为雷达工作频率，单位GHz，λ为载频波长，W≈1；A为幅度参数，B为相位参数，β₀为镜面反射参数，为角度参数；

在水体的式(9-b)中：其中，为入射余角，SS是水情级数，

θ_c＝sin^-1(λ/4πh_e)，h_e＝0.025+0.046SS^1.72，β＝[2.44(SS+1)^1.08]/57.29为镜面反射参数，W＝1.9。

6.根据权利要求3所述的一种基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，其特征在于，步骤4包括以下子步骤：

4a)在大地坐标系下，雷达位置坐标为P(l₀，b₀，h₀)，当杂波单元C_k在第一卦限时，则经度轴L的坐标大于纬度轴B的坐标，连接雷达位置P与杂波单元C_k对应的点，则与经度轴L和纬度轴B相交于J个点，根据公式(11)计算雷达位置P到杂波单元C_k投影的离散化点集{l_j，b_j，h_j}，j＝0，1，...，J：

当杂波单元C_k在第二卦限时，则纬度轴B的坐标大于经度轴L的坐标，连接雷达位置P与杂波单元C_k，则与经度轴L和纬度轴B相交于J个点，根据公式(12)计算雷达位置P到杂波单元C_k投影的离散化点集{l_j，b_j，h_j}，j＝0，1，...，J：

其中表示向下取整，j＝0时即为雷达坐标，j＝J即为杂波单元C_k的坐标；

当杂波单元C_k在其他三至八任一卦限中时，在经度轴L的坐标大于纬度轴B的坐标情况下，用式(11)计算离散化点集，在纬度轴B的坐标大于经度轴L的坐标情况下，用式(12)计算离散化点集；

4b)通过公式(4)先将雷达与离散化点集中的离散点转换到地心坐标系下，然后通过公式(13)计算得到雷达与离散化点集中的离散点之间的距离R_j：

$R_j=\sqrt{{(x_0-x_j)}^2+{(y_0-y_j)}^2+{(z_0-z_j)}^2}---\left(13\right)$

离散化点集中每一离散点相对于雷达的俯仰角θ_j表示为下式，j＝0，1，...，J，J表示离散化点集中离散点的数目；

${\mathrm{\&theta;}}_j={\sin }^{-1}\left(\frac{h_0-h_j}{R_j}\right),j=0,1,...,J---\left(14\right)$

4c)将离散化点集中的每一个俯仰角与杂波单元C_k的俯仰角进行比较，只要在离散化点集中存在一个俯仰角小于杂波单元C_k的俯仰角，则杂波单元C_k被遮挡，否则就不被遮挡，遮挡标志为δ_k：

7.根据权利要求1所述的一种基于数字高程和数字地表覆盖的机载雷达杂波仿真方法，其特征在于，步骤5包括以下子步骤：

回波空域信号导向矢量为

$s_{sk}\left(f_{sk}\right)={\&lsqb;1,e^{j2\mathrm{\&pi;}\&CenterDot;f_{sk}},...,e^{j2\mathrm{\&pi;}\&CenterDot;(N_a{N}_e-1)f_{sk}}\&rsqb;}^T---\left(16\right)$

回波时域信号导向矢量为

$s_{tk}\left(f_t\right)={\&lsqb;1,e^{j2\mathrm{\&pi;}\&CenterDot;f_{tk}},...,e^{j2\mathrm{\&pi;}\&CenterDot;(M-1)f_{tk}}\&rsqb;}^T---\left(17\right)$

其中，为空间频率，为归一化多普勒频率，N_a为水平向阵元个数，N_e为俯仰向阵元个数，M为一个脉冲重复间隔内发射脉冲的个数，λ为载频波长，f_r为脉冲重复频率，为雷达与杂波单元C_k之间的视线矢量，为雷达阵元间隔矢量，为雷达速度矢量；

杂波单元C_k的回波信号为：

$c_k={\mathrm{\&delta;}}_k\sqrt{{\mathrm{\&xi;}}_k}{s}_{sk}\&CircleTimes;s_{tk}---\left(18\right)$

其中，δ_k为遮挡标志，ξ_k为回波功率，表示Kronecker积；

将雷达照射范围内K个杂波单元的回波信号进行累积，得到杂波信号为：

$c=\underset{k=1}{\overset{K}{\&Sigma;}}{c}_k.---\left(19\right)$