基于正向法的仿三维场景SAR雷达回波的快速实现方法与流程

本发明属于雷达信号处理技术领域，特别涉及一种基于正向法的仿三维场景SAR雷达回波的快速实现方法，适用于大场景三维地形SAR雷达回波的模拟。

背景技术：

合成孔径雷达(SAR)能够全天时、全天候地对目标或场景进行观测，在灾害监测、资源勘查，尤其是军事和民用领域都有着巨大的应用空间。随着合成孔径雷达(SAR)技术日益发展，在对新型SAR设计方案和成像算法验证和评估时，需要满足特定参数的SAR回波信号，而这无法直接通过雷达录取回波得到，因此需要进行SAR回波仿真工作；SAR回波仿真分为逆向法和正向法，其中，逆向法是基于二维SAR图像进行回波录取，其简化了后向散射系数的计算，但丢失了场景的高程信息，最终获取的信息量只是基准SAR图像的子集，成像效果有限；正向法是基于原始地形结构的三维数据进行回波录取，可以实现场景的重构，在复杂地形探测、景象匹配制导中有重大作用。相比逆向法，正向法虽然复杂，但却拥有逆向法无法实现的优势，因此研究三维地面场景的正向法回波仿真是非常必要的。

然而对三维地面大场景目标进行正向法回波仿真时，面对的第一个难题就是如何有效地模拟后向散射的过程。关于此电磁散射的计算有很多方法，如物理光学法、几何绕射法等，都涉及到具体的散射机理，求解过程相对困难；考虑到合成孔径雷达常以地面场景为观测对象，因此可以采用实测数据结合雷达工作频率、波束射线侧偏角等因素，估算出后向散射系数的参考数值，这其中的计算量也是巨大的，此外还存在着阴影遮挡如何判断，大场景的回波仿真如何快速实现等问题。

技术实现要素：

针对以上现有技术存在的不足，本发明的目的在于提出一种基于正向法的仿三维场景SAR雷达回波的快速实现方法，使用该种基于正向法的仿三维场景SAR雷达回波的快速实现方法能够快速得到高信噪比的三维场景SAR雷达回波。

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

一种基于正向法的仿三维场景SAR雷达回波的快速实现方法，包括以下步骤：

步骤1，获取SAR雷达的载机飞行速度V、SAR雷达合成孔径方位时间T、SAR雷达的载机飞行方向和SAR雷达波束范围内的数字高程模型DEM数据；所述SAR雷达波束范围内的数字高程模型DEM数据包含K个点目标；K为自然数；

步骤2，对SAR雷达波束范围内的数字高程模型DEM数据进行插值处理，得到插值处理后的数字高程模型DEM数据；所述插值处理后的数字高程模型DEM数据包含K个点目标，然后将所述K个点目标随机分布在地理坐标系中的二维地形地形上，并任意选取其中一个点目标作为基准点目标，将该基准点目标坐标记为(x,y)，z(x,y)为该基准点目标对应的数字高程，计算得到时刻t时所述基准点目标在SAR雷达的载机上的累积概率密度函数F(t)，进而计算得到所述基准点目标在SAR雷达的载机上的三维地形的粗糙度方差σ；t为时间变量，K为自然数；

步骤3，根据基准点目标在SAR雷达的载机上的三维地形的粗糙度方差σ，对插值处理后的数字高程模型DEM数据进行精细化处理，得到精细化处理后SAR雷达波束范围内包含的DEM数据；

步骤4，将精细化处理后SAR雷达波束范围内的DEM数据划分为Q个小面元，其中第α个小面元的波束侧偏角为θ_α，然后分别计算波束侧偏角为θ₁时第1个小面元的后向散射系数σ₁，到波束侧偏角为θ_Q时第Q个小面元的后向散射系数σ_Q，并将所述后向散射系数σ₁到所述后向散射系数σ_Q的和，作为精细化处理后SAR雷达波束范围内的DEM数据的后向散射特性分布；其中，α∈{1,2,…,Q}，Q为自然数；

步骤5，得到精细化处理后SAR雷达波束范围内的DEM数据的后向散射特性分布后，利用改进的下视角比较法对Q个小面元进行阴影区域判定，进而计算得到经过阴影区域判定后Q个小面元的最终实际后向散射系数；

步骤6，根据经过阴影区域判定后Q个小面元的最终实际后向散射系数，分别计算Q个小面元各自的SAR雷达回波，然后使用线程外推技术对Q个小面元各自的SAR雷达回波进行累加，得到Q个小面元的SAR雷达回波，完成整个三维场景的回波仿真。

本发明的有益效果：本发明方法利用分形布朗运动模型FBM对原始DEM数据进行分形插值，保证了仿真精度，又利用小面元模型来计算后向散射系数；同时考虑到三维地形，本发明还进行了阴影区域的判断，最后本发明方法引入图形处理单元(Graphics Processing Unit,GPU)架构进行深度优化，极大地提高了仿真效率，同时满足三维场景下SAR雷达回波的实时处理需求。

附图说明

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

图1是本发明的一种基于正向法的仿三维场景SAR雷达回波的快速实现方法流程图；

图2(a)是原始DEM数据示意图；

图2(b)是分形插值结果示意图；

图3(a)是SAR雷达波束范围内的DEM数据在2^gM×2^gN维零矩阵中的分布图；

图3(b)是奇数位置坐标点插值结果示意图；

图3(c)是DEM数据分形插值三维结果示意图；

图4是本发明中剖分的小面元的空间几何关系示意图；

图5是使用基于下视角的阴影遮挡示意图；

图6(a)三维仿真场景设置示意图；横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元；

图6(b)轨迹1对应的SAR成像结果示意图；横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元；

图6(c)轨迹2对应的SAR成像结果示意图；横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元；

图7是本发明中实测数据仿真的某地区的实测DEM数据示意图；；其中，横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元，高度轴表示高程向，单位是采样单元；

图8(a)轨迹1对应的SAR成像结果示意图；横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元；

图8(b)轨迹2对应的SAR成像结果示意图；横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元；

图8(c)轨迹3对应的SAR成像结果示意图；横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元；

图8(d)轨迹4对应的SAR成像结果示意图；横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元。

具体实施方式

参照图1，为本发明的一种基于正向法的仿三维场景SAR雷达回波的快速实现方法流程图；该种基于正向法的仿三维场景SAR雷达回波的快速实现方法，包括以下步骤：

步骤1，确定图形处理单元，并从图形处理单元中获取SAR雷达参数，所述SAR雷达参数包括SAR雷达的载机飞行速度V、SAR雷达合成孔径方位时间T、SAR雷达的载机飞行方向和SAR雷达波束范围内的数字高程模型(DEM)数据；所述SAR雷达波束范围内的数字高程模型(DEM)数据包含K个点目标；K为自然数。

具体地，本发明的一种基于正向法的仿三维场景SAR雷达回波快速实现方法，首先在CPU处理器中设定SAR雷达参数，所述SAR雷达参数包括SAR雷达的载机飞行速度V、SAR雷达合成孔径方位时间T、SAR雷达的载机飞行方向和SAR雷达波束范围内包含的数字高程模型(Digital Elevation Model，DEM)数据，并据此分配图形处理单元(GPU)内的物理内存空间。

所述CPU处理器即常规计算机里的中央处理器，CPU处理器作为主机端，图形处理单元(GPU)是一个并行数据计算器；CPU处理器和图形处理单元(GPU)构成协同处理的工作模式，CPU处理器负责处理逻辑性强的事务处理和串行计算，图形处理单元(GPU)负责处理并行化的任务和数据计算。CPU处理器和图形处理单元(GPU)通过新一代总线或接口标准(PCI-Express，PCI-E)总线连接。其中，线程(thread)是图形处理单元(GPU)的最小单位，是一个相对独立可调度的执行单元；GPU中的线程块是由若干线程组成，同一个线程块中的所有线程享用一块内存，彼此通信或快速同步。

然后将所述SAR雷达的载机飞行速度V、SAR雷达合成孔径方位时间T、SAR雷达的载机飞行方向和SAR雷达波束范围内包含的DEM数据存入图形处理单元(GPU)中；参照图2(a)，为原始DEM数据示意图；所述SAR雷达波束范围内包含的DEM数据包含K个点目标；当SAR雷达波束的距离向单元数小于4096个点数时使用图形处理单元(GPU)中的共享内存；当SAR雷达波束的距离向单元数大于4096时调用图形处理单元(GPU)中的全局内存，所述全局内存和所述共享内存为GPU的高速存储器。

步骤2，对SAR雷达波束范围内的数字高程模型(DEM)数据进行插值处理，得到插值处理后的数字高程模型(DEM)数据；所述插值处理后的数字高程模型(DEM)数据包含K个点目标，然后将所述K个点目标随机分布在地理坐标系中的二维地形地形上，并任意选取其中一个点目标作为基准点目标，将该基准点目标坐标记为(x,y)，z(x,y)为该基准点目标对应的数字高程，计算得到时刻t时所述基准点目标在SAR雷达的载机上的累积概率密度函数F(t)，进而计算得到所述基准点目标在SAR雷达的载机上的三维地形的粗糙度方差σ；t为时间变量，K为自然数。

具体地，由于数字高程模型(DEM)数据具有统计自相似性(即随机分形特性)，因此数字高程模型(DEM)数据满足分形布朗运动(Fractional Brownian Motion，FBM)模型的统计特性；利用所述分形布朗运动(FBM)模型的统计自相似性对SAR雷达波束范围内包含的数字高程模型(DEM)数据进行插值处理，得到插值处理后的数字高程模型(DEM)数据；参照图2(b)，为分形插值结果示意图；所述插值处理后的数字高程模型(DEM)数据包含K个点目标，同时也是所述数字高程模型(DEM)数据的精细化处理结果；然后所述K个点目标随机分布在地理坐标系中的二维地形地形上，任意选取其中一个点目标作为基准点目标，并将该基准点目标坐标记为(x,y)，z(x,y)为该基准点目标对应的数字高程，计算得到时刻t时所述基准点目标在SAR雷达的载机上的累积概率密度函数F(t)，其表达式为：

$F\left(t\right)=P\[\frac{z(x+\mathrm{\Δ}x,y+\mathrm{\Δ}y)-z(x,y)}{\left|\right|(\mathrm{\Δ}x,\mathrm{\Δ}y)|{|}^H}\<t\]$

其中，P[·]表示概率分布函数，z(x,y)表示基准点目标对应的数字高程，(x,y)表示基准点目标在二维地形地形上的坐标，所述时刻t时所述基准点目标在SAR雷达的载机上的累积概率密度函数F(t)服从N(0,σ²)的高斯分布，σ表示三维地形的粗糙度方差，平坦区域的方差σ较小，陡峭区域方差的σ值较大；H表示基准点目标在SAR雷达的载机上的自相似参数，0<H<1，基准点目标在SAR雷达的载机上的自相似参数H说明地形的自相似程度，并能够用来描述地面的粗糙程度，平滑区域的自相似参数H值较大，粗糙区域的自相似参数H值较小；表示基准点目标在二维地形上的采样，△x表示基准点目标在二维地形上的横坐标间隔长度，△y表示基准点目标在二维地形上的纵坐标间隔长度；此外D为基准点目标在二维地形上的地形分形维数，且满足D＝3-H。

然后再利用所述分形布朗运动(FBM)模型的统计自相似性，并根据时刻t时所述基准点目标在SAR雷达的载机上的累积概率密度函数F(t)，计算得到所述基准点目标在SAR雷达的载机上的分形布朗模型等价式：

E[|z(x+△x,y+△y)-z(x,y)|]·||(△x,△y)||^-H＝C

其中，E(·)表示数学期望，z(x,y)表示基准点目标对应的数字高程，△x表示基准点目标在二维地形上的横坐标间隔长度，△y表示基准点目标在二维地形上的纵坐标间隔长度；C表示随机变量的均值，H表示基准点目标在SAR雷达的载机上的自相似参数，0<H<1。

对所述基准点目标在SAR雷达的载机上的分形布朗模型等价式两边同时取对数，得到所述分形布朗模型等价式的对数形式：

lnE[|z(x+△x,y+△y)-z(x,y)|]-Hln||(△x,△y)||＝lnC，进而得到以X＝ln||(△x,△y)||和Y＝lnE[|z(x+△x,y+△y)-z(x,y)|]为变量、以基准点目标在SAR雷达的载机上的自相似参数H为斜率的回归直线方程，△x表示基准点目标在二维地形上的横坐标间隔长度，△y表示基准点目标在二维地形上的纵坐标间隔长度，lnC表示基准点目标在SAR雷达的载机上的自相似参数H为斜率的回归直线方程在Y轴上的截距。

然后对所述分形布朗模型等价式的对数形式

lnE[|z(x+△x,y+△y)-z(x,y)|]-Hln||(△x,△y)||＝lnC进行线性回归模型拟合，计算得到基准点目标在SAR雷达的载机上的三维地形的粗糙度方差σ，

步骤3，根据基准点目标在SAR雷达的载机上的三维地形的粗糙度方差σ，对插值处理后的数字高程模型(DEM)数据进行精细化处理，得到精细化处理后SAR雷达波束范围内包含的DEM数据。

具体地，在确定了基准点目标在SAR雷达的载机上的自相似参数H和基准点目标在SAR雷达的载机上的三维地形的粗糙度方差σ后，对插值处理后的数字高程模型(DEM)数据进行精细化处理，所述精细化处理为分形插值处理，其过程为：

3.1：创建一个2^gM×2^gN维零矩阵，然后将SAR雷达波束范围内的DEM数据放置在所述2^gM×2^gN维零矩阵的行列坐标均为偶数的位置上，到包含K个点目标的2^gM×2^gN维矩阵。

其中，g为分形插值的迭代次数索引，g的初值为1，g∈{1,2,…,MAX}，MAX为设定的最大插值迭代次数；本实施例中MAX＝10；如图3(a)中所示，图3(a)为SAR雷达波束范围内的DEM数据在2^gM×2^gN维零矩阵中的分布图。

3.2：通过内插法计算得到所述包含K个点目标的2^gM×2^gN矩阵中所有行列坐标均为奇数处的高度值，然后将2^gM×2^gN矩阵中所有行列坐标均为奇数处的高度值和所述包含K个点目标的2^gM×2^gN维矩阵，作为第g次迭代后的2^gM×2^gN维插值矩阵。

具体地，所述包含K个点目标的2^gM×2^gN矩阵中所有行列坐标均为奇数处的高度值，其中将2^gM×2^gN矩阵中行列坐标(i,j)均为奇数处的高度值z(i,j)，其表达式为：

$z(i,j)=\frac{1}{4}\&lsqb;z(i-1,j-1)+z(i+1,j-1)+z(i-1,j+1)+z(i+1,j+1)\&rsqb;+\sqrt{1-2^{2H-2}}\left|\right|\mathrm{\&Delta;}x|{|}^H\mathrm{\&sigma;}\&CenterDot;G$

进而计算得到所述包含K个点目标的2^gM×2^gN维矩阵中行列坐标(i,j)均为奇数的位置上的数据，所述行列坐标(i,j)均为奇数的位置上的数据为所述包含K个点目标的2^gM×2^gN维矩阵中行列坐标(i,j)均为奇数处的高度值，如图3(b)中所示，图3(b)是奇数位置坐标点插值结果示意图。

其中，G代表服从标准正态分布的随机变量；||△x||表示插值后的2M×2N维矩阵的样本间隔，(i,j)表示所述包含K个点目标的2^gM×2^gN维矩阵中行列坐标均为奇数的点目标坐标，z(i-1,j-1)表示2^gM×2^gN维矩阵中行列坐标为(i-1,j-1)处的高度值，z(i+1,j-1)表示2^gM×2^gN维矩阵中行列坐标为(i+1,j-1)处的高度值，z(i-1,j+1)表示2^gM×2^gN维矩阵中行列坐标为(i-1,j+1)处的高度值，z(i+1,j+1)表示2^gM×2^gN维矩阵中行列坐标为(i+1,j+1)处的高度值，i为小于2M的正整数和j小于2N的正整数。

3.3：通过内插法计算得到所述包含K个点目标的2^gM×2^gN维矩阵中所有行列坐标中有一个为奇数处的高度值，然后将第1次迭代后的2^gM×2^gN维插值矩阵和所述2^gM×2^gN维矩阵中所有行列坐标中有一个为奇数处的高度值，作为第g次迭代后的最终2^gM×2^gN维插值矩阵，进而得到第g次迭代后的最终2_gM×2_gN维插值矩阵的样本间隔

具体地，所述包含K个点目标的2^gM×2^gN维矩阵中所有行列坐标中有一个为奇数处的高度值，其中将2^gM×2^gN矩阵中行列坐标(i',j')中有一个为奇数处的高度值记为z(i',j')，其表达式为：

$\hat{z}(i^{\&prime;},j^{\&prime;})=\frac{1}{4}\&lsqb;z(i^{\&prime;},j^{\&prime;}-1)+z(i^{\&prime;}-1,j^{\&prime;})+z(i^{\&prime;}+1,j^{\&prime;})+z(i^{\&prime;},j^{\&prime;}+1)\&rsqb;+2^{-H/2}\sqrt{1-2^{2H-2}}\left|\right|\mathrm{\&Delta;}\hat{x}|{|}^H\mathrm{\&sigma;}\&times;G$

其中，G代表服从标准正态分布的随机变量；表示第g次迭代后的最终2^gM×2^gN维插值矩阵的样本间隔，(i',j')表示所述包含K个点目标的2^gM×2^gN维矩阵中行列坐标有一个为奇数的点目标坐标，z(i',j'-1)表示2^gM×2^gN维矩阵中行列坐标为(i',j'-1)处的高度值，z(i'-1,j')表示2^gM×2^gN维矩阵中行列坐标为(i'-1,j')处的高度值，z(i'+1,j')表示2^gM×2^gN维矩阵中行列坐标为(i'+1,j')处的高度值，z(i',j'+1)表示2^gM×2^gN维矩阵中行列坐标为(i',j'+1)处的高度值，i'为小于2M的正整数，j'为小于2N的正整数。

3.4：设定仿真精度||△x_T||，如果第g次迭代后的最终2^gM×2^gN维插值矩阵的样本间隔大于或等于设定的仿真精度||△x_T||，则令g加1，返回子步骤3.2.

如果第g次迭代后的最终2^gM×2^gN维插值矩阵的样本间隔小于设定的仿真精度||△x_T||，则迭代停止，并将此时第g次迭代后的最终2^gM×2^gN维插值矩阵中包含的K个点目标和所有高度值，作为精细化处理后SAR雷达波束范围内包含的DEM数据，如图3(c)所示的DEM数据分形插值三维结果示意图。

步骤4，将精细化处理后SAR雷达波束范围内的DEM数据划分为Q个小面元，其中第α个小面元的波束侧偏角为θ_α，然后分别计算波束侧偏角为θ₁时第1个小面元的后向散射系数σ₁，到波束侧偏角为θ_Q时第Q个小面元的后向散射系数σ_Q，并将所述后向散射系数σ₁到所述后向散射系数σ_Q的和，作为精细化处理后SAR雷达波束范围内的DEM数据的后向散射特性分布；其中，α∈{1,2,…,Q}，Q为自然数。

具体地，结合实验数据与反映真实场景三维结构的DEM数据得到电波照射场景的后向散射系数。通常，三维地面可以看作是由大量平面小面片组成的，地面场景的电磁散射特性正是所有小面元的后向散射场相干叠加的结果，这种小面元剖分的思想是电磁计算方法的基础。本发明以四边形面元剖分为例。考虑到SAR雷达照射的目标是大范围的三维地面场景，所以在小面元剖分的时候，小面元的尺寸往往大于信号波长却远远小于SAR系统的分辨单元，这些要求均已在DEM数据的精细化处理中解决。每个小面元的空间几何参数由它的中心位置和法向矢量确定，它的后向散射系数可由下述方法近似得到。

4.1对精细化处理后SAR雷达波束范围内的DEM数据进行小面元剖分，得到Q个小面元，将Q个小面元随机分布在三维空间中，α∈{1,2,…,Q}，α表示第α个小面元，α的初值为1，Q表示精细化处理后SAR雷达波束范围内的DEM数据划分的小面元个数，每个小面元包含m×n个散射点；参照图4，为本发明中剖分的小面元的空间几何关系示意图，其中，坐标系为第α个小面元的局部坐标系，为坐标处第α个小面元的中心，X轴、Y轴、Z轴的方向与地理坐标系相同；S点为雷达载机位置；并设坐标处第α个小面元的散射点高度为

4.2计算得到坐标处小面元α的平面方程为

$q_{i_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}},j_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}}}=a\&times;i_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}}{\mathrm{\&Delta;x}}_{\mathrm{\&alpha;}}+b\&times;j_{i^{\&prime;\&prime;}}^{\mathrm{\&alpha;}}{\mathrm{\&Delta;y}}_{\mathrm{\&alpha;}}+c$

其中，a表示第α个小面元在X轴方向的斜率，b表示第α个小面元在Y轴方向的斜率，c表示第α个小面元的中心高度，△x_α表示第α个小面元沿X轴方向的采样间隔，△y_α表示第α个小面元沿Y轴方向的采样间隔，表示第α个小面元在坐标系O_i”-XYZ中的坐标，表示坐标处第α个小面元的中心。

然后利用最小二乘法计算小面元α的中心最佳方程minF(a,b,c)，其表达式为：

$\min F(a,b,c)=\underset{i^{\&prime;\&prime;}=0}{\overset{m-1}{\&Sigma;}}\underset{j^{\&prime;\&prime;}=0}{\overset{n-1}{\&Sigma;}}{(q_{i^{\&prime;\&prime;},j^{\&prime;\&prime;}}-h_{i^{\&prime;\&prime;},j^{\&prime;\&prime;}})}^2=\underset{i^{\&prime;\&prime;}=0}{\overset{m-1}{\&Sigma;}}\underset{j^{\&prime;\&prime;}=0}{\overset{n-1}{\&Sigma;}}{(a\&times;i^{\&prime;\&prime;}\mathrm{\&Delta;}x+b\&times;j^{\&prime;\&prime;}\mathrm{\&Delta;}y+c-h_{i^{\&prime;\&prime;},j^{\&prime;\&prime;}})}^2$

其中，h_i”,j”表示坐标(i”,j”)处第α个小面元的散射点真实高度，a表示第α个小面元在X轴方向的斜率，b表示第α个小面元在Y轴方向的斜率，c表示第α个小面元的中心高度，△x表示第α个小面元沿X轴方向的采样间隔，△y表示第α个小面元沿Y轴方向的采样间隔，(i”,j”)表示第α个小面元在坐标系O_i”-XYZ中的坐标，O_i”表示坐标(i”,j”)处第α个小面元的中心。

进而计算得到第α个小面元所在平面的单位法向矢量为

由雷达相位中心发出射向第α个小面元中心的电波的入射场矢量为则将电波在第α个小面元上的局部入射角，记为第α个小面元的波束侧偏角θ_α，SO_α表示雷达的载机位置到第α个小面元的中心长度，→表示矢量。

设第α个小面元在地面上的实际投影面积为S_α，则坐标处第α个小面元的面积为假定相同性质的电波垂直照射单位面积的该植被地面的后向散射系数的实验数据值为σ(0)，则计算得到波束侧偏角为θ_α时第α个小面元的后向散射系数为σ_α，

${\mathrm{\&sigma;}}_{\mathrm{\&alpha;}}\&ap;\mathrm{\&sigma;}\left(0\right)\&times;{\hat{S}}_{\mathrm{\&alpha;}}\&times;{\cos }^2\left({\mathrm{\&theta;}}_{\mathrm{\&alpha;}}\right).$

4.3令α加1，重复子步骤4.2，直到得到波束侧偏角为θ_Q时第Q个小面元Q的后向散射系数为σ_Q，并将此时得到的波束侧偏角为θ₁时第1个小面元的后向散射系数为σ₁，到波束侧偏角为θ_Q时第Q个小面元的后向散射系数为σ_Q的和，作为精细化处理后SAR雷达波束范围内的DEM数据的后向散射特性分布。

步骤5，得到精细化处理后SAR雷达波束范围内的DEM数据的后向散射特性分布后，利用改进的下视角比较法对Q个小面元进行阴影区域判定，进而计算得到经过阴影区域判定后Q个小面元的最终实际后向散射系数。

具体地，雷达发射机发出的电波在空间中沿着直线传播，传播过程中遇到山峰、高大建筑物等地形结构阻挡时，这些地形结构的背面就会接收不到电磁波，无法获得背面的后向散射系数，在SAR图像上出现阴影现象。对地形结构遮挡区域的判断是一个相当复杂的过程，它将导致回波模拟过程中计算量的急剧增加，特别是在随场景复杂度加深导致多重遮挡区域的判断中更为复杂，因此本发明利用了下视角比较法来判断三维场景阴影区域。

三维地形的阴影模拟示意图如图5所示，图5是使用基于下视角的阴影遮挡示意图；SAR雷达波束划分后某一方位角度间隔内的阴影模拟示意图。

5.1初始化：m'∈{1,2,…,Q}，m'表示第m'个小面元，也表示小面元的索引，m'的初值为1，Q表示精细化处理后SAR雷达波束范围内的DEM数据划分的小面元个数。

5.2任意选取第m'个小面元，并设定第m'个小面元地距为P_m'，计算得到第m'个小面元的下视角为β_m'，进而计算其余Q-1个小面元对应的Q-1个地距中地距小于第m'个小面元地距P_m'的Q'个小面元，Q'<＝Q-1，然后计算所述Q'个小面元各自的下视角，得到Q'个下视角，将所述Q'个下视角中的最大下视角记为β_max，若β_m'≤β_max，则第m'个小面元被遮挡；否则，第m'个小面元没有被遮挡；进而计算得到第m'个小面元的遮挡判定函数：

其中，为第m'个小面元的方位角度，r为第m'个小面元的地距。步骤4中通过剖分小面元的方法计算出第m'个小面元的后向散射系数为σ_m'，则计算得到经过阴影区域判定后第m'个小面元最终的后向散射系数为

5.3；令m'加1，重复子步骤5.2，直到得到经过阴影区域判定后第Q个小面元最终的后向散射系数为并将此时得到的经过阴影区域判定后第1个小面元最终的后向散射系数为到经过阴影区域判定后第Q个小面元最终的后向散射系数为作为经过阴影区域判定后Q个小面元的最终实际后向散射系数。

步骤6，根据经过阴影区域判定后Q个小面元的最终实际后向散射系数，并利用常规的同心圆算法分别计算Q个小面元各自的SAR雷达回波，然后使用线程外推技术对Q个小面元各自的SAR雷达回波进行累加，得到Q个小面元的SAR雷达回波，此时在图形处理单元GPU中完成了整个三维场景的回波仿真，最后将GPU中的SAR雷达回波信号发送至CPU处理器进行输出。

通过对以下实测数据的仿真实验来进一步验证本发明方法的有效性。

(一)理想三维模型的仿真实验

在阴影区域判断的仿真实验中，仿真平台和环境参数雷达基本参数如表一所示；选用如图6(a)所示的圆锥山丘场景作为基准图像，共有1024×1024个像素点；图6(a)三维仿真场景设置示意图；横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元；仿真圆锥山丘场景参数如表二所示。

表一仿真平台和环境参数

表二仿真圆锥山丘场景参数

利用本发明的基于正向法的仿三维场景SAR回波的快速实现方法，仿真过程中按轨迹1和轨迹2扫描场景。仿真过程中都扫描1024×1024个点目标。然后利用常规的线频调变标算法对产生的回波数据分别进行成像处理，得到如图6(b)-图6(c)所示的两幅成像结果，两幅图像均已做了地距校正；其中，图6(b)表示按轨迹1扫描得到的SAR成像结果；图6(c)表示按轨迹2扫描得到的SAR成像结果；横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元；

从仿真结果上来看，SAR雷达的载机沿着预设轨迹1飞行后的SAR成像结果如图6(b)所示，图中的阴影区域位于图像正下方，且阴影的朝向也是正下方，这与轨迹1的照射相对应一致；SAR雷达的载机沿着预设轨迹2飞行后的SAR成像结果如图6(c)所示，图中的阴影区域位于右下方，朝向右下角，这与轨迹2的照射相对应的。

(二)实际DEM数据的仿真实验

利用某地区的实测DEM数据进行仿真验证及相关分析。仿真仍然选用表一所示的仿真平台及环境参数，SAR回波模拟时选用的插值后的DEM数据的基本信息见下表三所示，为了对比阴影效果，预设的SAR雷达载机飞行轨迹的基本信息如表四所示，表中的H代表载机飞行高度，β代表波束中心线下视角，θ代表飞行轨迹与X轴的夹角；实测DEM数据的插值后的三维图如图7所示，图7是本发明中实测数据仿真的某地区的实测DEM数据示意图；其中，横轴表示方位向，单位是采样单元，纵轴表示距离向，单位是采样单元，高度轴表示高程向，单位是采样单元。

表三实测DEM数据的场景参数

表四载机飞行轨迹的基本信息

为了对比区分不同照射角度下的SAR雷达回波的差异，将SAR雷达的载机沿着预设飞行轨迹1～轨迹4的SAR回波的成像结果显示于图8(a)-图8(d)中；图8(a)轨迹1对应的SAR成像结果示意图；图8(b)轨迹2对应的SAR成像结果示意图；图8(c)轨迹3对应的SAR成像结果示意图；横轴表示方位向，图8(d)轨迹4对应的SAR成像结果示意图；图8(a)-图8(d)中的横轴均表示方位向，单位是采样单元；纵轴均表示距离向，单位是采样单元。

图8(a)～图8(d)分别是载机在飞行轨迹1～轨迹4下模拟的SAR回波的成像结果，四幅图均已经做了地距校正。对比图8(a)与图8(b)、图8(c)与图8(d)可以发现，虽然SAR雷达高度、波束中心的下视角相同，但由于飞行方向不同，导致阴影效果不同，所以成像结果中有着细节差异；对比图8(a)与图8(c)、图8(b)与图8(d)可以发现，由实测DEM数据的高程值相比载机高度不是很大，上下两幅图的纹理基本相同，但阴影效果有稍微差异，下视角为60°的阴影效果明显一些，这些事实充分说明了本发明方法的准确性。

综上所述，仿真实验验证了本发明的正确性，有效性和可靠性。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围；这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。