CN114594438A - 一种双基合成孔径雷达回波模拟方法 - Google Patents

Abstract

本发明公开了一种双基合成孔径雷达回波模拟方法，本发明的方法通过分别判断阴影，利用真实场景的DEM数据进行模拟，能够反映出双基SAR成像过程中的阴影几何畸变特点；针对在平坦地形区域回波模拟的问题，通过加入地物分类数据，利用不同地表的散射特性不同，实现了平坦地形区域的回波模拟。本发明的方法综合利用DEM数据和地物分类数据，提高了对SAR回波模拟的准确度，能够对任意地区进行SAR回波模拟，整个模拟过程还原了合成孔径雷达实际工作过程，对验证双基SAR系统参数指标、最优化设计和验证成像算法有重要的意义。

Description

一种双基合成孔径雷达回波模拟方法

技术领域

本发明属于合成孔径雷达(Synthetic Aperture Radar，SAR)技术领域，具体涉及SAR回波模拟的技术领域。

背景技术

合成孔径雷达具有穿透性、全天时和全天候工作能力的独特优点，目前已被广泛应用在地球遥感、资源勘探、侦察、测绘、灾情预报等领域。在进行SAR雷达系统方案设计和成像算法验证时，都需要符合特定要求的大量SAR回波数据。然而由于飞行成本等因素的限制，不可能通过直接的飞行试验来获取回波数据，所以通过模拟来获得所需的SAR回波数据是一种重要的手段。

目前SAR回波模拟的方法主要有两种，一种是利用光学图像或真实的雷达图像，将其灰度图作为场景后向散射系数的估计值，进行原始信号的仿真。如文献《扩展场景的SAR回波信号快速仿真算法》(遥测遥控，2005,26(6)：33-38)中，将光学照片应用于SAR回波信号的仿真，并以此验证时域和频域的回波生成算法。此方法可以快速模拟回波信号，具有效率高和处理简洁的特点，然而无法反映目标场景的真实地形特征，以及无法模拟雷达图像特有的几何畸变现象。

另一种方法是先根据真实场景的三维地物结构，计算目标后向散射系数，然后模拟雷达回波信号。如文献“具有地面起伏特征的三维场景回波模拟研究》(雷达科学与技术，2007,5(5)：342-248)”中，给出了一种针对有起伏地形的三维场景回波模拟方法。该方法基于数字高程数据(DEM)拟合小面单元，并计算小面单元的后向散射系数，然后逐点生成SAR回波信号。该方法只适用于单基SAR的回波模拟，然而对于双基SAR，由于发射和接收方向均可能产生阴影，阴影产生规律与单基SAR不同，因此该方法不适用于双基SAR；同时该方法只利用了DEM数据来获取场景的后向散射系数，因此无法反映均质平坦地形区域的模拟效果。

发明内容

针对现有SAR回波模拟技术中存在的上述问题，本发明提出了一种基于DEM和地物分类数据的双基SAR回波信号模拟方法。

本发明的具体技术方案为：一种双基合成孔径雷达回波模拟方法，具体包括如下步骤：

A.获取待模拟区域的DEM数据，并对DEM数据进行分形插值，使其数据间隔小于或等于SAR模拟系统的分辨率；

B.获取待模拟区域的地物分类数据，并对地物分类数据进行最近邻插值，得到与步骤A中区域大小一致的地物分类数据；

C.利用步骤A中得到的DEM插值数据拟合小面单元；

D.根据步骤C中得到的拟合小面单元系数，计算各小面单元的表面积、局部入射角和局部散射角；

E.根据步骤B中插值后的地物分类数据和雷达工作频段，确定步骤C中得到的各个拟合小面单元的单位面积初始散射系数；再根据步骤D得到的拟合小面单元局部入射角和局部散射角，确定小面单元的单位面积散射系数；

F.根据SAR平台和场景的空间几何关系，计算阴影函数；

G.根据步骤D、E、F中得到的各拟合小面单元的表面积、单位面积散射系数和阴影函数，计算待模拟区域的散射系数；

H.根据步骤G中得到的待模拟区域散射系数模拟回波信号。

本发明的有益效果：本发明的双基SAR回波信号模拟方法通过分别判断阴影，利用真实场景的DEM数据进行模拟，能够反映出双基SAR成像过程中的阴影几何畸变特点；针对在平坦地形区域回波模拟的问题，通过加入地物分类数据，利用不同地表的散射特性不同，实现了平坦地形区域的回波模拟。本发明的方法综合利用DEM数据和地物分类数据，提高了对SAR回波模拟的准确度，能够对任意地区进行SAR回波模拟，整个模拟过程还原了合成孔径雷达实际工作过程，对验证双基SAR系统参数指标、最优化设计和验证成像算法有重要的意义。

附图说明

图1为本发明实施例中的双基SAR几何配置示意图。

图2为本发明实施例中双基SAR回波信号模拟的流程示意图。

图3为本发明实施例中待模拟区域的DEM数据图。

图4为本发明实施例中同一区域的地物分类数据图。

图5是本发明实施例中DEM插值数据拟合小面单元示意图。

图6是本发明实施例步骤D中计算局部入/散射角示意图。

图7是本发明实施例阴影判断示意图。

图8为具体实施例中待模拟区域的阴影函数图。

图9为具体实施例中计算得到的该区域散射系数仿真图。

图10为具体实施例中计算得到的该区域回波数据幅度仿真图。

图11为具体实施例中计算得到的该区域的模拟SAR图像。

具体实施方式

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

本实施例中双基SAR的几何配置如图1所示，发射机以正侧视照射待模拟区域，接收机下前视接收回波信号，C波段(f_c＝5.3GHz)机载雷达仿真参数如表1所示。

表1

参数	符号	数值	单位
				发射机高度	H<sub>t</sub>	10000	m
发射机速度	v<sub>t</sub>	150	m/s
				接收机高度	H<sub>r</sub>	8000	m
接收机速度	v<sub>r</sub>	150	m/s
				载波频率	f<sub>c</sub>	5.3	GHz
脉冲时宽	T<sub>r</sub>	5	μs
				信号带宽	B	10	MHz
距离采样率	F<sub>r</sub>	20	MHz
				脉冲重复调频率	PRF	1000	Hz
距离采样点数	Nfast	4096
				方位采样点数	Nsolw	61440

在一个xyz空间坐标系中，O表示坐标系的原点，目标场景中心位置为X_c＝0,Y_c＝30500m，Z_c＝0；雷达发射机以速度v_t＝150m/s，高度H_t＝10000m沿X轴方向匀速直线飞行，波束指向为正侧视。接收机以速度v_r＝150m/s，高度H_r＝8000m与发射机同向匀速直线飞行，其位置在发射机右后方，其中，Y_r＝Y_c＝30500m，X_r＝-30000m，波束指向为下前视；雷达发射信号波长λ＝0.0566m，脉冲时宽T_r＝5μs，带宽B＝10MHz。

本发明实施例的基于DEM和地物分类数据的双基SAR回波信号模拟方法，流程示意图如图2所示，包括如下步骤：

A.获取待模拟区域的DEM数据，并对DEM数据进行分形插值，使其数据间隔小于或等于SAR模拟系统的分辨率；

B.获取待模拟区域的地物分类数据，并对地物分类数据进行最近邻插值，得到与步骤A中区域大小一致的地物分类数据；

C.利用步骤A中得到的DEM插值数据拟合小面单元；

D.根据步骤C中得到的拟合小面单元系数，计算各小面单元的表面积、局部入射角和局部散射角；

E.根据步骤B中插值后的地物分类数据和雷达工作频段，确定步骤C中得到的各个拟合小面单元的单位面积初始散射系数；再根据步骤D得到的拟合小面单元局部入射角和局部散射角，确定小面单元的单位面积散射系数；

F.根据SAR平台和场景的空间几何关系，计算阴影函数；

G.根据步骤D、E、F中得到的各拟合小面单元的表面积、单位面积散射系数和阴影函数，计算待模拟区域的散射系数；

H.根据步骤G中得到的待模拟区域散射系数模拟回波信号。

步骤A中，待模拟区域的DEM如图3(a)所示，对应的三维地形如图3(b)所示。这里从地理空间数据云(www.gscloud.cn)上获取待模拟区域30m分辨率的DEM(数字高程)数据，然后对DEM数据进行分形插值。

本实施例中采用领域点递归中值置换法来实现分形插值，这里需要先计算出待模拟区域的分形参数，然后确定插入点的高程值。设DEM数据中某一点坐标为(x,y)，z(x，y)为该点的数字高程，根据分形布朗模型：

其中，F(ζ)是随机变量ζ的概率分布函数，服从N(0，F_σ ²)的高斯分布；F_H为自相似参数，用来描述地面的粗糙程度；||(Δx，Δy)||表示样本的间隔，本实施例中Δx＝Δy。

根据分形布朗场的随机统计特性，式(1)可等价为：

其中，E(·)表示数学期望，常数C表示随机变量|ζ|的均值，有

对式(2)两边同时取对数，可以得到：

ln E[|z(x+Δx，y+Δy)-z(x，y)|]-F_Hln||(Δx，Δy)||＝ln C (3)

式(3)可以看作是以X＝ln||(Δx，Δy)||与Y＝lnE[|z(x+Δx，y+Δy)-z(x，y)|]为变量，以F_H为斜率的直线方程，ln C即为直线在Y轴上的截距。因此，对于给定的DEM数据，以不同的间隔i(i＝1,2,…,10)算出不同的X_i及对应的Y_i，再用线性回归模型拟合点对(X_i，Y_i)，求出拟合直线的斜率k和Y轴截距b；根据分形布朗模型可知：F_H＝k，ln C＝b，则

本实施例中原始的DEM数据矩阵为300×300，首先建立一个600×600的空矩阵，将原来的DEM数据放置在新矩阵行列坐标都为奇数的位置，然后通过分形内插公式得到剩余点的高程值。具体地，当行列坐标都为偶数时，其高程值可通过下式计算得到：

当行列坐标只有一个为偶数时，其高程值可通过下式计算得到：

其中，G表示服从标准正态分布的随机变量；F_σ和F_H为待模拟区域的分形参数；||Δx||为插值后的样本间距。

在本实施例中，为减少计算量，只对DEM数据进行一次分形插值，即插值后的数据间隔为15m。插值后的待模拟区域DEM如图3(C)所示,对应的三维地形如图3(d)所示。

步骤B中，与步骤A中相同区域的地物分类数据如图4(a)所示。这里从清华大学(http://data.ess.tsinghua.edu.cn/)提供的全球30m分辨率地物覆盖数据集中，获取待模拟区域的地物分类数据，然后对地物分类数据进行最近邻插值。

本实施例中插值倍数为2，假设原图像上的像素点坐标为(srcX,srcY)，插值后目标图像上的像素点坐标为(dstX,dstY)，则目标图像与原图像上像素点的对应关系为：

srcX＝round(dstX/2) (5)

srcY＝round(dstY/2) (6)

式中round表示将值四舍五入为整数。即插值后的各像素点灰度值用原图像中最近邻像素点的灰度值来填充，整体实现了对原图像2倍的放大处理。进行一次最近邻插值后，待模拟区域的地物分类数据如图4(b)所示。

本实施例中，步骤C的具体过程如下：

如图5所示，设DEM中某一区域的小面单元由相邻的4个点A、B、C、D组成，相应的高程分别为z₁、z₂、z₃、z₄，DEM的分辨率为ρ。假设四个点A、B、C、D的坐标分别为(x₁,y₁,z₁)、(x₂,y₂,z₂)、(x₃,y₃,z₃)、(x₄,y₄,z₄)，其中x₁＝x₂，x₃＝x₄，y₁＝y₄，y₂＝y₃，设这四个点拟合得到的平面为z＝ax+by+c，利用最小二乘原理，建立目标函数

令

可计算得到拟合小面单元的系数a、b、c。

本实施例中，步骤D的具体过程如下：

步骤C中拟合小面单元的法向量为

与z轴的夹角(即小面的坡度角γ)的余弦

由投影定理可知小面单元的表面积为

如图6所示，雷达入射波对拟合小面单元的下视角

其中，z₀为拟合小面单元中心点的高程，z₀＝a·x₀+b·y₀+c，

H_t为发射机的飞行高度，d_t为发射机到小面中心点的地距。雷达散射波对拟合小面单元的下视角

其中，H_r为接收机的飞行高度，d_r为接收机到小面中心点的地距。因此雷达入射波向量为(0，sinβ_t，-cosβ_t)，雷达散射波向量为(-sinβ_r，0，cosβ_r)。将雷达入射波向量和拟合小面单元的法向量做点积运算，可得到拟合小面单元的局部入射角：

将雷达散射波向量和拟合小面单元的法向量做点积运算，可得到拟合小面单元的局部散射角：

本实施例中，步骤E的具体过程如下：

本实施例中采用经验散射模型，由经验散射模型可知，各拟合小面单元的单位面积的散射系数σ⁰随局部入射角θ_t和局部散射角θ_r变化的公式为：

σ⁰(θ_t，θ_r)＝σ⁰(0)·cosθ_t·cosθ_r (7)

式中的σ⁰(0)为入(散)射角为0度(即垂直小面照射)时，小面元的单位面积初始散射系数，是与地表材质有关的量。

地物分类数据的值(范围为1-10)代表不同的地物类型，可根据表2确定各个拟合小面单元所对应的地物类型，由此确定各拟合小面单元对应的单位面积初始散射系数σ⁰(0)。

表2

本实施例中，根据待模拟区域包含的地物类型，确定其在雷达波长为C波段(本例中波长为λ＝0.0566m)的条件下，对应的单位面积初始散射系数如表3所示。

表3

本实施例中，步骤F的具体过程如下：

由于双基SAR的发射装置和接收装置不在同一平台，因此成像时可能出现双影。如图7所示，发射机和接收机的阴影判断都分为两步进行：

本实施例中采用局部入/散射角判断和下视角比较法来计算阴影。对于发射机产生的阴影，首先判断发射机波束能否照射到小面单元。根据第(i，j)个小面单元的局部入射角θ^t _i，j(即步骤D中计算的θ_t)，如果

则认为小面单元能够被发射机波束照射，若

则认为小面单元不能被发射机波束照射，即该小面单元为阴影区。

若小面单元能够被发射机波束照射，则需进一步判断。设小面单元中心点P(i，j)，计算雷达入射波对第(i，j)个小面单元的下视角β^t _i，j(即步骤D中计算的β_t)，然后按照多普勒中心线方向计算比P(i，j)点地距小的所有小面中心点的下视角，找到其中最大值β^t _max。如果β^t _max＞β^t _i，j，则认为P(i，j)点没有被遮挡，若β^t _max＞β^t _i，j，则认为P(i，j)点被遮挡，即该小面单元为阴影区。此处定义一个发射机阴影函数：

对于接收机产生的阴影，首先判断小面单元能否散射雷达波束给接收机。根据第(i，j)个小面单元局部散射角θ^r _i，j(即步骤D中计算的θ_r)，如果

则认为小面单元能够散射雷达波束给接收机，若

则认为小面单元不能散射雷达波束给接收机，即该小面单元为阴影区。

若小面单元能够散射雷达波束给接收机，则需进一步判断。对同一拟合小面单元中心点P(i，j)，计算雷达散射波对第(i，j)个小面单元的下视角β^r _i，j(即步骤D中计算的β_r)，然后按照散射波束中心线方向计算比P(i，j)点地距小的所有小面中心点的下视角，找到其中最大值β^r _max。如果β^r _max≤β^r _i，j，则认为P(i，j)点没有被遮挡，若β^r _max＞β^r _i，j，则认为P(i，j)点被遮挡，即该小面单元为阴影区。此处定义一个接收机阴影函数：

因此双基SAR的整体阴影函数为：

本实施例中计算得到待模拟区域的阴影函数如图8所示。

本实施例中，步骤G的具体过程如下：

根据步骤E得到的拟合小面单元的单位面积散射系数，则第(i，j)个小面单元的单位面积散射系数σ⁰ _i，j(θ_t，θ_r)，再考虑每个小面的表面积大小A_i，j，可以计算得到单个拟合小面的散射系数：

σ_i，j＝A_i，j·σ⁰ _i，j(θ_t，θ_r) (8)

考虑阴影之后的散射系数为：

ε_i，j＝σ_i，j·δ_i，j (9)

由此方法计算所有拟合小面单元的散射系数，即可得到整个待模拟区域的散射系数。

将得到的散射系数图显示成灰度图像，进行归一化处理后，图像强度值变换到0-255之间，最终得到待模拟区域的散射系数灰度图如图9所示。

本实施例中，步骤H的具体过程如下：

本实施例中SAR发射的信号为线性调频(LFM)脉冲信号：

其中，f_c表示信号的载波频率，rect[·]为矩形窗函数，K_r为调频率，τ为距离向快时间变量；T_r为发射信号的脉冲宽度。

实际中雷达发射的是一串线性调频脉冲信号：

其中，T_p为发射脉冲串的间隔周期；n为发射信号的脉冲序号；Nslow为发射信号脉冲数(也即方位采样点数)。

假设雷达发射机发射的信号为s_pul(τ)，发射信号遇到目标散射产生回波，雷达接收机收到的单个目标的回波信号为：

其中，ε为散射系数，η为方位向慢时间变量，R(η)为目标与雷达收发平台的距离之和；W_t[·]为发射天线单程方向加权函数，W_r[·]为接收天线单程方向加权函数。

假设地面场景有M×N个散射点，地面场景的回波信号是M×N个散射点回波信号的矢量和，则整个场景的回波信号为：

其中，ε_i，j表示第(i，j)位置目标的散射系数；R_i，j(η)为第(i，j)位置目标与雷达收发平台的距离之和。

本实施例中，通过回波信号模拟方法得到的待模拟区域回波信号幅度图如图10所示。为验证模拟回波数据的正确性，采用BP(后向投影)成像算法对生成的目标区域回波数据进行成像，得到的成像结果如图11所示。

采用本实施例中的双基SAR回波模拟方法，得到的模拟图像能够反映出双基SAR的成像特点，有助于对双基SAR系统几何构型和工作模式进行验证。而在相同场景下采用本实施方式的方法，与只采用DEM数据进行回波模拟的传统方法相比，本实施方式的方法对平坦地形区域的模拟效果更好，可以区分出平坦区域上不同的地物类型，以获得更多的有用信息，使得模拟的结果更加接近真实的SAR图像。

Claims (6)

1.一种双基合成孔径雷达回波模拟方法，具体包括如下步骤：

A.获取待模拟区域的DEM数据，并对DEM数据进行分形插值，使其数据间隔小于或等于SAR模拟系统的分辨率；

B.获取待模拟区域的地物分类数据，并对地物分类数据进行最近邻插值，得到与步骤A中区域大小一致的地物分类数据；

C.利用步骤A中得到的DEM插值数据拟合小面单元；

D.根据步骤C中得到的拟合小面单元系数，计算各小面单元的表面积、局部入射角和局部散射角；

E.根据步骤B中插值后的地物分类数据和雷达工作频段，确定步骤C中得到的各个拟合小面单元的单位面积初始散射系数；再根据步骤D得到的拟合小面单元局部入射角和局部散射角，确定小面的单位面积散射系数；

F.根据SAR平台和场景的空间几何关系，计算阴影函数；

G.根据步骤D、E、F中得到的各拟合小面单元的表面积、单位面积散射系数和阴影函数，计算待模拟区域的散射系数；

H.根据步骤G中得到的待模拟区域散射系数模拟回波信号。

2.根据权利要求1所述的一种双基合成孔径雷达回波模拟方法，其特征在于，步骤C的具体过程如下：

设DEM中某一区域的小面单元由相邻的4个点A、B、C、D组成，相应的高程分别为z₁、z₂、z₃、z₄，假设四个点A、B、C、D的坐标分别为(x₁,y₁,z₁)、(x₂,y₂,z₂)、(x₃,y₃,z₃)、(x₄,y₄,z₄)，其中x₁＝x₂，x₃＝x₄，y₁＝y₄，y₂＝y₃，设这四个点拟合得到的平面为z＝ax+by+c，利用最小二乘原理，建立目标函数

令

计算得到拟合小面单元的系数a、b、c。

3.根据权利要求2所述的一种双基合成孔径雷达回波模拟方法，其特征在于，步骤D的具体过程如下：

步骤C中拟合小面单元的法向量为

与z轴的夹角(即小面的坡度角γ)的余弦

由投影定理可知小面单元的表面积为

雷达入射波对拟合小面单元的下视角

其中，z₀为拟合小面单元中心点的高程，z₀＝a·x₀+b·y₀+c，

H_t为发射机的飞行高度，d_t为发射机到小面中心点的地距，雷达散射波对拟合小面单元的下视角

其中，H_r为接收机的飞行高度，d_r为接收机到小面中心点的地距，因此雷达入射波向量为(0，sinβ_t，-cosβ_t)，雷达散射波向量为(-sinβ_r，0，cosβ_r)，将雷达入射波向量和拟合小面单元的法向量做点积运算，得到拟合小面单元的局部入射角：

将雷达散射波向量和拟合小面单元的法向量做点积运算，可得到拟合小面单元的局部散射角：

4.根据权利要求3所述的一种双基合成孔径雷达回波模拟方法，其特征在于，步骤F的具体过程如下：

对于发射机产生的阴影，首先判断发射机波束能否照射到小面单元，根据第(i，j)个小面单元的局部入射角θ^t _i，j，如果

则认为小面单元能够被发射机波束照射，若

则认为小面单元不能被发射机波束照射，即该小面单元为阴影区；

若小面单元能够被发射机波束照射，则需进一步判断：设小面单元中心点P(i，j)，计算雷达入射波对第(i，j)个小面单元的下视角β^t _i，j，然后按照多普勒中心线方向计算比P(i，j)点地距小的所有小面中心点的下视角，找到其中最大值β^t _max，如果β^t _max≤β^t _i，j，则认为P(i，j)点没有被遮挡，若β^t _max＞β^t _i，j，则认为P(i，j)点被遮挡，即该小面单元为阴影区；此处定义一个发射机阴影函数：

对于接收机产生的阴影，首先判断小面单元能否散射雷达波束给接收机，根据第(i，j)个小面单元局部散射角θ^r _i，j，如果

则认为小面单元能够散射雷达波束给接收机，若

则认为小面单元不能散射雷达波束给接收机，即该小面单元为阴影区；

若小面单元能够散射雷达波束给接收机，则需进一步判断：对同一拟合小面单元中心点P(i，j)，计算雷达散射波对第(i，j)个小面单元的下视角β^r _i，j，然后按照散射波束中心线方向计算比P(i，j)点地距小的所有小面中心点的下视角，找到其中最大值β^r _max，如果β^r _max≤β^r _i，j，则认为P(i，j)点没有被遮挡，若β^r _max＞β^r _i，j，则认为P(i，j)点被遮挡，即该小面单元为阴影区；此处定义一个接收机阴影函数：

因此双基SAR的整体阴影函数为：

5.根据权利要求4所述的一种双基合成孔径雷达回波模拟方法，其特征在于，步骤G所述的计算待模拟区域的散射系数具体过程如下：

根据步骤E得到的拟合小面单元的单位面积散射系数，则第(i，j)个小面单元的单位面积散射系数σ⁰ _i，j(θ_t，θ_r)，考虑每个小面的表面积大小A_i，j，计算得到单个拟合小面的散射系数：

σ_i，j＝A_i，j·σ⁰ _i，j(θ_t，θ_r)

考虑阴影之后的散射系数为：

ε_i，j＝σ_i，j·δ_i，j

由此计算所有拟合小面单元的散射系数，即可得到整个待模拟区域的散射系数。

6.根据权利要求5所述的一种双基合成孔径雷达回波模拟方法，其特征在于，步骤H的具体过程如下：

SAR发射的信号为线性调频脉冲信号：

其中，f_c表示信号的载波频率，rect[·]为矩形窗函数，K_r为调频率，τ为距离向快时间变量；T_r为发射信号的脉冲宽度；

雷达发射的是一串线性调频脉冲信号：

其中，T_p为发射脉冲串的间隔周期；n为发射信号的脉冲序号；Nslow为发射信号脉冲数；

假设雷达发射机发射的信号为s_pul(τ)，发射信号遇到目标散射产生回波，雷达接收机收到的单个目标的回波信号为：

其中，ε为散射系数，η为方位向慢时间变量，R(η)为目标与雷达收发平台的距离之和；W_t[·]为发射天线单程方向加权函数，W_r[·]为接收天线单程方向加权函数；

假设地面场景有M×N个散射点，地面场景的回波信号是M×N个散射点回波信号的矢量和，则整个场景的回波信号为：

其中，ε_i，j表示第(i，j)位置目标的散射系数；R_i，j(η)为第(i，j)位置目标与雷达收发平台的距离之和。

Images