DEM是实际地形高程信息的离散模型，主要描述地表地貌形态的空间分布。有别于传统地形图，DEM将表征实际地形的地理信息以离散的数字形式按照一定的结构组织在一起。对DEM数据进行分析和计算可以提取出海拔、坡度、地面起伏和地表粗糙度等描述地貌类型的地形因子。根据这些地形因子，可以对地表的地貌类型进行划分，并根据划分结果计算散射单元的后向散射系数，最终实现杂波的建模。考虑到计算的便捷性，选择以规则栅格数据形式存储的DEM数据来提取海拔、坡度、地面起伏和地表粗糙度等描述地貌类型的地形因子。这类数据网格间距相等，可以直接根据原始坐标和网格间距确定每一栅格的平面坐标。

DEM是由一系列的离散数据点组成，每一个数据点表示该地点的高程值，即海拔。因此，可以直接利用DEM数据获得每一个地面点的海拔数据。

地面坡度是一个矢量，由地面各点的坡度组成，可用过该点的切平面与水平地面的夹角表示，反映了地面在该点的倾斜程度。地面每一个点都有坡度，它是点上的概念。若用$ H = f(x, y) $表示地表曲面函数，H为地面某点的高程值，x为东西方向，y为南北方向，则坡度为：

式中，$ p = \partial H/\partial x $，$ q = \partial H/\partial y $分别是高程H在x方向上和y方向上的变化率。

对于每一个地面栅格点，首先确定其分析窗口的大小，本文以栅格点H_{i, j}为中心，以包含周围相邻8个栅格点的3×3窗口作为坡度的分析窗口，如表 1所示。采用文献[12]提出的三阶反距离平方权差分算法来计算x和y方向上的高程变化率p、q，再代入式(1)中，就可以得到该点的坡度。

由于采用的DEM数据是以矩阵形式存储的规则栅格数据，因此，将分析窗口进行上下左右平移就可方便地获得整块区域的坡度数据。

地形起伏描述的是一定区域内地形的起伏变化情况，是一个区域地形特征的一个宏观指标。地形起伏是用一定区域内高程值的最大值与最小值之间的差值来描述的，有：

式中，TR_i代表该区域的地形起伏；A_i为分析区域；H_max、H_min分别是该区域的最大和最小高程值。

由式(2)可以看出，分析区域A_i对地形起伏因子的值有直接的影响，因此，选择一个合适的分析区域，对地形起伏信息的有效性起着关键性的作用，直接影响到对地面起伏特征的描述是否准确。文献[13]表明分析区域的形状对地形起伏的影响较小，而分析区域的大小则对地形起伏有明显的影响，甚至可能影响到地形的分类。地形起伏最大值与分析区域的面积之间存在着某种依赖关系，且存在一个明显的分界点，使分析区域超过分界点后，地形起伏不再有明显的变化。因此，以分界点界定的区域作为分析区域是一种最佳选择。文献[14]表明全中国区域内的最佳分析统计区域单元面积为21 km²。考虑到分析区域形状对地形起伏的影响不大和DEM数据的特点，将分析区域设置为矩形且单元面积取为21 km²，来提取地形起伏因子TR_i。

地表粗糙度主要用来描述一定区域内的地表侵蚀程度，也是一个区域地形特征的一个宏观指标。它被定义为地形单元中4个点(H_{i-1, j-1}、H_{i-1, j+1}、H_{i+1, j-1}、H_{i+1, j+1})组成的矩形的两条对角线的两个中点之间的距离，本文用D表示地表粗糙度。

根据D的定义可得：

类似地，可以取一个分析窗口从DEM数据中提取地表粗糙度因子，窗口的大小根据DEM数据的水平分辨率来决定。

本文主要根据地貌对杂波功率谱的实质性影响对地形进行分类，与常规根据地理地貌学分类地形的方法不同。如按地理地貌学、平原和草原应该分类为不同的地貌，但这两种地貌对杂波功率谱(后向散射系数)产生的影响几乎没有差别，因此把两者划分为同一类。根据提取出的海拔高度、地面坡度、地形起伏和地表粗糙度等地形因子数据，再结合中国陆地基本地貌形态类型的分类类别^{[13, 14, 15, 16, 17, 18]}，把地貌分为5种类型，分别为海面、平原、森林、丘陵以及高山，并分别量化编码为0、1、2、3、4。地貌类型编码及其对应的地形因子数值如表 2所示。

根据表 2可以看出，在地貌类型分类中，平原、森林和丘陵之间存在着交叉，这是由于实际场景中很少存在单一的地形，丘陵地区一般都覆盖着树木等植被。因此，在根据地形因子进行地貌分类时，首先依据海拔高度进行粗划分，然后根据坡度、地形起伏和地表粗糙度等其他地形因子进行细划分。

一个地点的方位主要由经度和纬度决定。由于每一份DEM数据是以该份数据的经纬度起始点为名保存的，本身不包含经纬度信息，又考虑到每份数据的经纬度跨度是确定的，因此，根据起始点经纬度、横向和纵向的数据点数以及终点经纬度信息，可以计算得到每一份DEM数据中每个高程数据对应的经纬度信息，即：

式中，N_i为具体的经度(或纬度)信息；N₀为该份DEM数据的起始点的经度(或纬度)；L_d是该份DEM数据的经度(或纬度)终点与起始点的差值；N_k为该份DEM数据的经度(或纬度)方向上的数据点数。

由于地貌在较大的时间跨度内可认为是稳定不变的，并且计算不同地形因子的数据在来源、类型和格式等方面存在差异性，再加上一些地形因子的计算十分耗时，因此需要把构建杂波模型的地貌信息事先存贮起来，建立不同信息之间的内在联系，提高先验信息的利用率和后续杂波建模的实时性。本文主要采用开源数据库MySQL平台来存贮和管理地貌信息。

由于需要存储的信息是未处理的DEM数据和经过处理之后获得的地形分类信息，所以库中一条数据记录被设计为由编号、经度、纬度、高程值和地形类型码(ID)等字段构成。其中，经度和纬度主要用于确定一个地点的方位，高程值和地形用于后续杂波建模，如擦地角的计算。为了方便存储和数据管理，数据表设计为二维数组，每行代表一条记录，每列由一个字段组成。由于知识辅助的机载认知雷达自适应处理过程具有读多写少操作的特点，所以选择MyISAM为存储引擎对信息进行更新、删除、添加和提取等操作。

在实际应用中，本文仅从该数据库中提取与载机当前位置和雷达扫描范围相关的记录，并重组为DEM数据和地形分类编码信息后用于杂波功率的计算，以降低运算量和提高处理的实时性。

假设目标信号和干扰信号均为零，收发共用同一天线，雷达发射功率为P_t，雷达收发天线增益为G_t=G_t=G，收发天线的有效接收面积为A_r=A_t，信号波长为λ，散射单元截面积为σ，则根据雷达方程^[19]，天线接收的杂波功率计算模型为如下形式：

由式(5)可知，基于DEM数据进行杂波建模的关键是如何利用DEM提供的真实地貌信息计算杂波散射单元的截面积σ。σ主要与杂波散射单元的散射系数和杂波散射单元的几何面积有关，有：

式中，σ⁰为该杂波单元的后向散射系数；S_m为该杂波单元的几何面积。

影响杂波后向散射系数的因素有很多，包括雷达工作频率、波束擦地角以及地面环境等。而地面环境因素又包括波束照射区域的地貌类型和天气情况等。后向散射系数模型就是后向散射系数与这些影响因素之间的关系模型。

与地形相关的散射系数模型有F.T.Ulaby模型、Morchin模型、修正的Morchin模型^[20]以及Barton模型^[21]等。这几种散射模型都包含了与地面场景有关的参数，得到的散射系数将随地形的变化而变化。考虑载机在实际环境中的擦地角为0°~90°，采用修正的Morchin模型来估算地杂波散射系数。

修正的Morchin模型为^[20]：

式中，$ {\theta _g} $为擦地角；当地面类型为沙漠并且$ {\theta _g} < {\theta _c} $时，$ \sigma _c^0 = {\theta _g}/{\theta _c} $；对于其他地面类型或地面类型为沙漠，而 ${\theta _g} > {\theta _c} $时，$ \sigma _c^0 = 1 $； ${\theta _c} = \arcsin (\lambda /4{\text{\pi }}{h_e}) $；$ {h_e} \approx 9.3\beta _0^{2.2} $；$\mu = \sqrt {{f_0}/4.7} $；$ {f_0} $为雷达工作频率，单位为GHz；A、B和β₀是与地貌类型密切相关的参数，其数值可以利用第2.2节提供的地形分类信息确定。

由此可知，只要获得波束擦地角，可计算出杂波的后向散射系数。然而，实际雷达波束照射的区域存在起伏，地面某一散射单元的擦地角并不一定等于雷达下视角，擦地角的计算需要考虑到杂波散射单元与载机的位置关系。为了获得准确的地面场景散射系数，需要计算场景中每个散射单元(DEM数据点)相对于载机位置的擦地角。为此，本文首先从数据库中提取出散射单元的高程和坐标信息，再根据散射单元到载机位置的向量与散射单元平面的法线向量之间的关系，确定散射单元的擦地角，具体方法如下：从DEM数据中取任意相邻的3个点，假设分别为A(x_a, y_a, z_a)，B(x_b, y_b, z_b)，C(x_c, y_c, z_c)，这3个点组成一个三角形平面。令平面法向量u=(x, y, z)。由于该法向量u分别与三角形平面的AB和AC边垂直，则有：

通过变换可以得到平面法向量为：

式中，d_s为DEM数据的分辨率。

根据载机位置(x_p, y_p, z_p)与散射单元坐标(x_t, y_t, z_t)，可得它们之间的视线向量ww。与u之间的夹角即为雷达波束入射角，其计算公式如下：

一个散射单元的u与w之间的关系为：

式中，ε为雷达波束入射角。因此，该散射单元的擦地角为：

在获得了各点的后向散射系数之后，只要计算出散射单元面积，就可求得杂波散射单元的几何面积S_m。为此，设计了一种网格映像法来计算S_m，具体步骤如下：

1)对雷达照射区域进行网格划分。采用等距离环划分法，在雷达照射区域内将地表划分为多个ΔR ×Δϕ的栅格单元，其中ΔR为最小距离环宽度，Δϕ为最小方位角间隔。如图 1中的扇形代表一个栅格单元。

图  1  杂波散射单元几何面积计算过程示意图

2)统计一个栅格单元内点散射体的个数。图 1中黑色点代表来自于DEM的杂波数据，只有落入一个栅格单元内的杂波数据，才对杂波散射单元的几何面积有贡献。由于距离环按弧形划分，使一个栅格单元呈环形，难以用计算方法统计其中的数据点数，因此以一个栅格单元为基础划定一个最小矩形，把一个栅格单元内的数据点映像到最小矩形中，将最小矩形中的数据点个数视为一个栅格单元内的数据点数，如图 1所示。

3)将一个栅格单元的点散射体相加，形成一个新的散射体。因此，杂波散射单元的几何面积就转化为栅格单元内所有点散射体面积的相干叠加。根据DEM数据的特点，一个栅格单元对应的散射体面积为：

式中，d_s是DEM数据分辨率；N表示散射单元内点散射体的个数。

最后，将计算得到的杂波后向散射系数和单元几何面积代入式(6)，得到杂波单元的雷达散射截面积σ，再将结果代入式(5)便可计算出杂波功率。

综上所述，本文提出的杂波建模方法流程如图 2所示。

图  2  杂波建模方法流程图

本文对上述杂波建模的主要过程进行了仿真实验，以说明本文方法的可行性和有效性。首先，为了验证地貌类型划分的有效性，本文对地形因子算法以及地形量化分类算法进行了仿真。为了验证数据库设计的合理性，对数据库数据的存贮以及读取重组算法进行了仿真。为了验证选择的后向散射系数模型是否符合实际，本文仿真了擦地角对地形后向散射系数的影响。最后对选定范围内的杂波功率进行了仿真，用以验证整个算法流程的可行性以及杂波建模方法的准确性。

本文根据某地的DEM数据形成的原始地形图如图 3所示。该地形图范围为：北纬26°～27°，东经119°～120°。所选区域为沿海某地，在该区域内既存在海洋环境，同时也有山地等地形的存在，基本覆盖了大部分地形类型，具有较好的普适性。

图  3  根据某地DEM数据形成的原始地形图

图 4是图 3分别基于坡度、地表粗糙度和地形起伏因子显示的地形图。从图 4a的坡度图可以看出，对于山地地形，由于地形复杂，坡度的变化也较大，根据坡度值的大小可以大致区分出陡峭区域和平缓区域。坡度值较小的为平缓区域(图中蓝绿色所示)，较大的为陡峭区域(图中橙色所示)。海平面高程值为0，因此，其坡度值也显示为0(图中深蓝色所示)。从图 4b的地表粗糙度图可以看出，地势平缓区域的地表粗糙度变化比较平缓，相对值较小(图中蓝色所示)，而山区丘陵地带的地表粗糙度变化明显，相对值较大(图中蓝绿色所示)。海岛地区由于区域面积小，地表粗糙度变化不明显，相对值也较小(图中右下角蓝色所示)。1 km²和21 km²的地形起伏图分别如图 4c和图 4d所示。对比分析可知，当分析单元面积为21 km²时，地形起伏因子图(图 4d)可以较明显地体现图 3中的几个山体，但是对于其右下角部分的海岛，显示效果却很一般，这是由于此时的分析区域面积相对于海岛面积过大而导致的。然而，当分析单元面积为1 km²时，图 3中的海岛信息则可以在地形起伏图(图 4c中被明显的显示出来。因此，分析区域的大小对地形起伏影响比较大。

图  4  基于地形因子显示的地形图

根据表 2的分类依据，本文对图 3所示的DEM地形图进行地貌类型划分，结果如图 5所示。图 5中，深蓝色的区域为海面，淡蓝色的区域为平原，浅绿色的区域为森林，橙黄色的区域为丘陵，红棕色区域为高山。从地貌类型划分结果看，图 5中大部分地区为丘陵地带，森林基本伴随着丘陵边缘分布，平原以一定的块状区域分布在丘陵中。通过与该地区的卫星地图进行比较，发现基于地貌类型量化分类地貌的结果与实际区域的地貌类型基本吻合，只是在平原的分类中存在微小差别。其主要原因是把农田、草地以及水泥地等地貌均划分到平原的类型中，没有细分。事实上，把平原划分为大类是有科学依据的，后面的仿真表明，农田、草地和水泥地等地形随擦地角的变化和数值大小基本一致，因此从杂波建模的角度，它们可以被划分到平原地形中。

图  5  基于地貌量化分类的地形图

在载机坐标为[26.623°N, 119.210°E]和雷达扫描范围为20 km的条件下，本文从数据库中提取出DEM地形图和量化分类的地形图，并与完整的原图进行比较，发现提取的地形图特征和经纬度范围与原地形的一致。因此，本文数据库存储的信息能够准确地反映载机所处的真实地貌，可以有效地提高杂波建模的准确性。

本文进一步测试了MyISAM存储引擎存取一份DEM数据样本(数据点数为1 601×1 601)的时间，并在相同运行环境下与InnoDB存储引擎进行了比较，结果如表 3所示。由此可知，本文采用MyISAM的读取效率明显优于InnoDB。

不同地形散射系数随擦地角的变化曲线如图 6所示。除了高山地形以外，大部分地形的散射系数与擦地角之间的关系曲线走势相同。当雷达波束擦地角较小时，高山的反射系数最大，其他几种地形的差异比较小；当擦地角较大时，大部分地形的反射系数增大很快，而高山增大不多，甚至略有降低。这主要是由于擦地角较小时，高山对雷达波速的镜面反射起主要作用；当擦地角较大时，大部分地形的镜面反射成分占了主导地位，因此散射系数迅速增加。图 6结果符合实际情况，说明采用修正的Morchin模型计算地杂波散射系数具有合理性。

图  6  不同地形散射系数随擦地角的变化曲线

假设载机坐标为[26.623°N, 119.210°E]，载机高度为5 000 m，最大扫描距离为20 000 m，最小扫描距离为5 500 m，发射频率为1 GHz，距离分辨率为90 m，方位角划分为100份。基于DEM数据提取的地形图及其对应的空间杂波功率谱如图 7所示。

图  7  基于DEM数据计算的杂波功率谱

根据载机的坐标数据可以看出，载机位于图 7a的中点(图中箭头所示)，且由西向东进行扫描(箭头指向方向)。由于设定了扫描距离的范围，使得俯仰角范围缩小为(15°～16°)，因此该仿真实验只模拟了载机俯仰角范围内的杂波功率谱，如图 7b所示。对比图 7a和图 7b可知，在接近载机处，该区域地形相对平坦，镜面反射占主导地位，因此在图 7b中出现了强杂波。当俯仰角在25°～55°区间时，在方位角为30°和-30°左右，存在两条较强杂波区域，对应到图 7a中，刚好此时在载机方位角30°线上存在一条山脊，而在载机方位角-30°线存在一条峡谷，这两种地貌的存在产生了较强的杂波。进一步从图 7b中可知，在俯仰角15°～20°范围内有比较强烈的杂波，此时波束扫描距离接近最大值，而雷达波束以半圆形向外发射，对应到图 7a中，该区域的地形起伏剧烈，上半部分(方位角0°～90°)基本由山地构成(红色区域)，下半部分(方位角0°～90°)海拔落差达到了300m左右，上下部分结合，近似形成了一圈山脊，而这圈山脊起到了类似于镜面反射的作用，因此会有较剧烈的杂波存在。综上所述，本文构建的杂波谱与真实地貌之间存在一一对应关系，能够反映真实地貌的变化。

本文通过分析和挖掘DEM数据，提取了海拔、坡度、地面起伏和地表粗糙度等4种地形因子用于划分地貌类型，并在MySQL平台上建立了地形信息数据库，存储DEM数据和相关地貌信息。仿真结果表明，本文提取的地形因子和划分的地貌类型符合实际场景地形，数据存取方便快捷。在此基础上，本文建立了一种基于DEM数据的地杂波建模方法。仿真结果表明，采用该方法构建的地杂波模型符合机载认知雷达所处的实际场景，具有真实、稳定和可靠的优势，当用作机载认知雷达的先验信息时可有效辅助机载认知雷达提高工作性能和生存力。