CN101876704A - 干涉合成孔径雷达三维陆地场景回波仿真方法 - Google Patents

Abstract

本发明提供一种干涉合成孔径雷达三维陆地场景的回波仿真方法，技术方案是：首先，对三维陆地场景进行几何建模，把三维陆地场景划分为彼此相邻的小面单元。然后，以小面单元为散射元进行电磁建模，计算散射元的后向散射系数。最后，根据干涉合成孔径雷达系统信号模型，通过计算单条回波，并行生成回波。本发明计算量小，并可采用并行计算方式，能够在现有计算能力下在较短的时间内完成；同时本方法的计算精度高，能够满足干涉合成孔径雷达系统仿真与研究的要求。

Description

干涉合成孔径雷达三维陆地场景回波仿真方法

技术领域

本发明属于微波遥感和信号处理的交叉技术领域，特别涉及一种干涉合成孔径雷达针对三维陆地场景的回波仿真方法。

背景技术

干涉合成孔径雷达能够全天候、全天时地获取具有较高空间分辨率和高程精度的地球表面三维地形图。利用三维地形图信息，可监测地球陆地表面和冰雪表面的变化情况，进行地震、火山爆发、滑坡和大洪水等灾害预测，为农业、林业、渔业等生产提供帮助，为军事活动提供信息支持。总之，利用干涉合成孔径雷达获得的三维地形图具有重要的应用价值和广泛的应用领域。

早期的干涉合成孔径雷达主要用于机载系统，在应用范围和应用空间有较大的局限性。20世纪90年代随着技术的进步，干涉合成孔径雷达成功应用于星载系统，为全球提供了丰富的干涉合成孔径雷达的星载数据。2000年2月11日美国国家航空航天局和国家影像与测绘局联合进行了航天飞机雷达地形测绘任务，这次基于干涉合成孔径雷达技术的试验首次实现了真正意义上的全球地形三维测绘，试验的部分数据和结果已经公布以供各科研机构处理和评估。2007年，德国发射了分辨率达1米的X波段合成孔径雷达卫星TerraSAR-X，这是迄今为止分辨率最高的民用合成孔径雷达卫星。德国计划于2010年发射第二颗几乎一样的卫星，两颗合成孔径雷达卫星将组成双星编队的TanDEM-X系统进行全球高精度地形测绘，公布的预期指标将高于美国的航天飞机雷达地形测绘任务。目前德国已经将部分的TerraSAR-X的重复轨道干涉合成孔径雷达数据在互联网公布供全球研究机构使用。

虽然有国外的多种实际数据资源可以利用，但它们的主要用途在于促进干涉合成孔径雷达的信号处理研究。目前，我国处于干涉合成孔径雷达系统的研发阶段，更关心的是系统层面的问题，如：总体参数选取、系统性能评估、系统误差分析与补偿等。具有高逼真度、高效率的仿真系统将为系统总体设计提供有力工具，仿真的灵活性和可控性能够为系统不同方案的对比试验提供便利，能够为修正系统的理论分析模型提供依据，能够为信号处理研究提供丰富的、针对性的数据源。对于干涉合成孔径雷达系统，三维陆地场景原始回波仿真是干涉合成孔径雷达系统仿真中最耗时、最基础的环节，是构建整个高精度、高效率仿真系统的核心。

干涉合成孔径雷达系统回波仿真的耗时与仿真场景中目标的数目成正比。在进行系统的分辨率等性能验证和成像算法研究中，点数在10²量级的点目标或点目标阵列的回波仿真可满足研究要求。这时即使采用传统的无近似误差的时域叠加方法，仿真的耗时也仅需要数十秒。但是，在对干涉测高任务功能的仿真分析中，研究对象主要是大范围三维陆地场景，点目标或点目标阵列的点数可达10⁸量级，此时应用传统的时域叠加方法需要数百年时间，仿真已经失去实际意义。因此研究基于星载干涉合成孔径雷达系统的三维陆地场景高精度、高效雷达回波仿真方法具有重要意义。

在干涉合成孔径雷达系统回波仿真方法的研究中，意大利学者Franceschetti1992年提出了一种基于二维傅立叶变换的快速回波模拟方法；2002年，提出了一种基于二维傅立叶变换的聚束式合成孔径雷达的快速回波模拟算法；在2005年，他对非匀速直线运动情况下的快速回波模拟方法进行了研究。上述这些方法针对回波信号模拟中某一个具体的问题，有着很大的局限性：难以体现卫星、地球、目标之间的相对运动，不易注入多种系统误差，模拟的精度达不到实际运用要求，不能得到普适性的应用。Alessandro于2004年采用直接时域叠加的方法进行干涉合成孔径雷达的回波模拟，该方法同Franceschetti的方法相比，存在计算效率低下的问题。国内在这方面做的工作与上述方法类似，都在精度和效率之间进行折衷处理，难以取得较好的平衡。并且上述研究大多针对回波仿真中的某个具体问题，未能形成完整、实用、适用于系统仿真的回波仿真方法和流程。

发明内容

本发明要解决的技术问题是提供干涉合成孔径雷达三维陆地场景的回波仿真方法，本方法计算量小，并可采用并行计算方式，能够在现有计算能力下在较短的时间内完成；同时本方法的计算精度高，能够满足干涉合成孔径雷达系统仿真与研究的要求。

本发明技术方案的思路是：首先，对三维陆地场景进行几何建模，把三维陆地场景划分为彼此相邻的小面单元。然后，以小面单元为散射元进行电磁建模，计算散射元的后向散射系数。最后，根据干涉合成孔径雷达系统信号模型，通过计算单条回波，并行生成回波。

本发明技术方案是：

利用已知的三维陆地场景的数字高程模型数据，完成以下步骤：

第一步：几何建模

本步骤对数字高程模型数据进行多次插值、小面单元划分和平面拟合，建立三维陆地场景的几何模型。三维陆地场景的几何模型由彼此相邻的小面单元组成，每个小面单元由其法矢量和高程描述。

第(一)步，数字高程模型数据插值

若数字高程模型数据的采样间距为Δx_origin，干涉合成孔径雷达的地距分辨率为Δρ，需要进行的插值次数N_interp为

$N_{\mathrm{interp}}=\mathrm{ceil}\left({\log }_2\left(\frac{{9\mathrm{\Δx}}_{\mathrm{origin}}}{\mathrm{\Δ\ρ}}\right)\right)$

其中，ceil(x)为：

ceil(x)＝k，k-1＜x≤k，k为整数

对数字高程模型数据进行N_interp次的插值。每次插值的算法如下：

设本次插值前数字高程模型数据的点数为I×J，且任一点的位置序号为(i₀，j₀)，1≤i₀≤I，1≤j₀≤J，高程为g(i₀，j₀)；本次插值后数字高程模型数据点数为(2I-1)×(2J-1)，且任一点位置序号为(i′，j′1≤i′≤2I-1，1≤j′≤2J-1，高程为f(i′，j′)；用下述步骤计算f(i′，j′)：

对i′和j′都为奇数，用下式进行赋值：

$f(i^{\′},j^{\′})=g(\frac{i^{\′}+1}{2},\frac{j^{\′}+1}{2})$

对i′和j′都为偶数，用下式进行赋值：

$f(i^{\′},j^{\′})=\frac{1}{4}[f(i^{\′}-1,j^{\′}-1)+f(i^{\′}+1,j^{\′}-1)+f(i^{\′}-1,j^{\′}+1)+f(i^{\′}+1,j^{\′}+1)]+$

$2^{-H/2}\sqrt{1-2^{2H-2}}{\left|\right|\mathrm{\Δx}\left|\right|}^H\mathrm{\σ}{G}_N$

对i′和j′只有一个为偶数，用下式进行赋值：

$f(i^{\′},j^{\′})=\frac{1}{4}[f(i^{\′},j^{\′}-1)+f(i^{\′}-1,j^{\′})+f(i^{\′}+1,j^{\′})+f(i^{\′},j^{\′}+1)]+$

$2^{-H/2}\sqrt{{1-2}^{2H-2}}{\left|\right|\mathrm{\Δx}\left|\right|}^H{\mathrm{\σG}}_N$

上述两个公式中的σ为本次插值前数字高程模型数据的标准差，H是三维陆地场景的频谱指数，G_N为N(0，1)分布的随机变量，Δx为本次插值前数字高程模型数据的采样间距。

第(二)步，平面拟合

对经过N_interp次插值后数字高程模型数据，以3行×3列的点阵为基本单位进行划分，形成若干个小面单元。

对任意一个小面单元，设该小面单元内任一点的位置序号为(i_s，j_s)，1≤i_s≤3，1≤j_s≤3，高程为f_s(i_s，j_s)，记：

$X=\left[\begin{array}{ccc}1& 0& 0\\ 1& 0& 1\\ 1& 0& 2\\ 1& 1& 0\\ 1& 1& 1\\ 1& 1& 2\\ 1& 2& 0\\ 1& 2& 1\\ 1& 2& 2\end{array}\right],$ $Y=\left[\begin{array}{c}{f}_s\left(\mathrm{1,1}\right)\\ f_s\left(\mathrm{1,2}\right)\\ f_s\left(\mathrm{1,3}\right)\\ f_s\left(\mathrm{2,1}\right)\\ f_s\left(\mathrm{2,2}\right)\\ f_s\left(\mathrm{2,3}\right)\\ f_s\left(\mathrm{3,1}\right)\\ f_s\left(\mathrm{3,2}\right)\\ f_s\left(\mathrm{3,3}\right)\end{array}\right],$

用下式计算：

$\mathrm{\β}=\left[\begin{array}{c}{\mathrm{\β}}_1\\ {\mathrm{\β}}_2\\ {\mathrm{\β}}_3\end{array}\right]={\left(X^{\′}X\right)}^{-1}{X}^{\′}Y$

可求得该小面单元的法矢量

高程z_s＝β₁。

第二步，电磁建模

本步骤通过计算三维陆地场景的几何模型中每个小面单元的雷达散射截面积，建立三维陆地场景的电磁模型。三维陆地场景的电磁模型中，每个小面单元由其雷达散射截面积描述。

设三维陆地场景的中心在地心固连坐标系下的坐标为(x_c，y_c，z_c)，经过第一步得到小面单元的数量为M×N，小面单元距离向间隔为ρ_x，方位向间隔为ρ_y。

对任一小面单元，设其法矢量为n₀(m，n)，位置序号为(m，n)，1≤m≤M，1≤n≤N，则该小面单元在地心固连坐标系下的坐标(x_t(m，n)，y_t(m，n)，z_t(m，n))为(x_c+(m-M/2)ρ_x，y_c+(n-N/2)ρ_y，z_c+z_s(m，n))，进行如下步骤的计算：

首先，计算入射角θ_n(m，n)，计算公式为：

${\mathrm{\θ}}_n(m,n)=\frac{n_0(m,n)\·n_i(m,n)}{\left|n_0(m,n)\right|\left|n_i(m,n)\right|}$

其中，n_i(m，n)为波束方向矢量，设本次仿真的中心时刻干涉合成孔径雷达天线相位中心在地心固连坐标系下的坐标为(x_a，y_a，z_a)，则

$n_i(m,n)=\left[\begin{array}{c}{x}_t(m,n)-x_a\\ y_t(m,n)-y_a\\ z_t(m,n)-z_a\end{array}\right]$

然后，采用下式计算后向散射系数σ⁰(m，n)：

σ⁰(m，n)＝P₁+P₂exp[-P₃θ_n(m，n)]+P₄cos(P₅θ_n(m，n)+P₆)

其中，P₁～P₆为与地形有关的参数。

最后，采用下式计算雷达散射截面积σ_F(m，n)：

σ_F(m，n)＝σ⁰(m，n)ρ_xρ_y。

第三步：回波生成

本步骤先生成单条回波，再并行计算多条回波，从而生成三维陆地场景的仿真回波。

设干涉合成孔径雷达的发射信号为p(τ)，τ为快时间。

第(一)步，生成单条回波

对每个慢时间l，进行如下步骤的计算：

第(1)步，计算系统传输函数

首先，计算每个小面单元的斜距，采用下述方法：

即利用下式计算位置序号为(m，n)的小面单元的斜距：

$R(l,m,n)=\sqrt{{[x_l-x_t(m,n)]}^2+{[y_l-y_t(m,n)]}^2+{[z_l-z_t(m,n)]}^2}$

上式中，(x_l，y_l，z_l)为在慢时间l时，天线相位中心在地心固连坐标系下的坐标。

然后，确定每个小面单元的增益总和，采用下述方法：

位置序号为(m，n)的小面单元的增益总和G(l，m，n)为

$G(l,m,n)=\frac{\sqrt{P_t{G}_t{G}_r}}{2\sqrt{\mathrm{\π}}\mathrm{\λR}{(l,m,n)}^2}$

其中，P_t为干涉合成孔径雷达的发射功率，G_t为发射天线增益，G_r为接收天线增益，λ为发射信号波长。

接着，确定每个小面单元的冲激响应时，采用下述方法：

对位置序号为(m，n)的小面单元，计算冲激信号的幅度A(l，m，n)和相位为：

A(l，m，n)＝σ_F(m，n)G(l，m，n)

上式中，f₀为干涉合成孔径雷达的中心频率，c为光速；

采用延迟近似处理计算冲激响应的延迟Δτ(l，m，n)为：

$\mathrm{\Δ\τ}(l,m,n)=\mathrm{round}\left(\frac{2R(l,m,n)}{\mathrm{c\ΔT}}\right)\mathrm{\ΔT},$

round(x)＝k，k-0.5≤x＜k+0.5，k为整数

上述ΔT为干涉合成孔径雷达接收信号的采样时间间隔，

可得到冲激响应为

其中，δ(τ)为单位冲激函数。

最后可以得到在慢时间l的系统传输函数：

$h(l,\mathrm{\τ})=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{h}_s(l,m,n,\mathrm{\τ}).$

第(2)步，计算回波信号

首先，对发射信号p(τ)进行快速傅利叶变换，得到发射信号的频谱P(f)；对系统传输函数h(l，τ)进行快速傅利叶变换，得到系统传输函数的频谱H(l，f)。

然后，把发射信号频谱P(f)和系统传输函数的频谱H(l，f)相乘，得到回波的频谱(l，f)。

最后，对(l，f)进行快速傅利叶反变换，得到在慢时间l的单条回波s(l，τ)。

第(二)步，并行计算多条回波

因为每个慢时间生成的单条回波是相对独立的，因此本步骤可在不同的计算机上并行计算。

首先计算在所有慢时间的单条回波，再把每个慢时间生成的单条回波按时间顺序合成，得到多条回波。

作为对本发明的进一步改进，对高程为常数z₀时的三维陆地场景，可以采用斜距递推方法计算斜距，减少计算量：

将相邻的M_r行×N_r列个小面单元划分为一个斜距计算单元，其中：

$M_r=\mathrm{floor}\left(\frac{15000\mathrm{\λ}}{{\mathrm{\ρ}}_x}\right)$

$N_r=\mathrm{floor}\left(\frac{15000\mathrm{\λ}}{{\mathrm{\ρ}}_y}\right)$

其中，floor(x)为：

floor(x)＝k，k≤x＜k+1，k为整数。

在每个斜距计算单元内，递推计算斜距，采用下述方法：

首先，设当前斜距计算单元的中心小面单元位置序号为(m_c，n_c)，计算其慢时间l时的斜距R(l，m_c，n_c)：

$R(l,m_c,n_c)=\sqrt{{[x_l-x_t(m_c,n_c)]}^2+{[y_l-y_t(m_c,n_c)]}^2+{[z_l-z_t(m_c,n_c)]}^2}$

然后，用下列式子计算递推系数A至E：

$A=-\frac{x_l-x_t(m_c,n_c)}{R(l,m_c,n_c)}{\mathrm{\ρ}}_x,$ $B=-\frac{y_l-y_t(m_c,n_c)}{R(l,m_c,n_c)}{\mathrm{\ρ}}_y$

$C=(-\frac{{(x_l-x_t(m_c,n_c))}^2}{R^3(l,m_c,n_c)}+\frac{1}{R({l,m}_c,n_c)}){{\mathrm{\ρ}}_x}^2,$ $D=(-\frac{{(y_l-y_t(m_c,n_c))}^2}{R^3(l,m_c,n_c)}+\frac{1}{R(l,m_c,n_c)}){{\mathrm{\ρ}}_y}^2$

$E=-\frac{\mathrm{AB}}{R(l,m_c,n_c)}$

最后，递推计算当前斜距计算单元内其他小面单元的斜距：

依次令q＝1，2…(N_r-n_c-1)，用下式计算位置序号为(m_c，n_c+q+1)的小面单元在慢时间l时的斜距：

R(l，m_c，n_c+q+1)＝R(l，m_c，n_c+q)+qD+(D/2+B)

依次令q＝-1，-2…1-n_c，用下式计算位置序号为(m_c，n_c+q-1)的小面单元在慢时间l时的斜距：

R(l，m_c，n_c+q-1)＝R(l，m_c，n_c+q)-qD+(D/2+B)

依次令q＝1-n_c，2-n_c…(N_r-n_c-1)，p＝1，2…(M_r-m_c-1)，用下式计算位置序号为(m_c+p+1，n_c+q)的小面单元在慢时间l时的斜距

R(l，m_c+p+1，n_c+q)＝R(l，m_c+p，n_c+q)+pC+(qE+C/2+A)

依次令q＝1-n_c，2-n_c…(N_r-n_c-1)，p＝-1，-2…1-m_c，用下式计算位置序号为(m_c+p-1，n_c+q)的小面单元在慢时间l时的斜距

R(l，m_c+p-1，n_c+q)＝R(l，m_c+p，n_c+q)-pC+(qE+C/2-A)。

采用本发明可取得以下技术效果：

本发明的回波仿真方法中，回波生成步骤中的每个慢时间生成单条回波与其它慢时间生成单条回波是独立的，可以分别完成，因此本发明的回波仿真方法十分有利于进行并行计算，可以充分发挥集群计算机的优势，设计并行程序大幅提高程序运行效率。

本发明的回波仿真方法在回波生成这一步中采用了延迟近似处理，通过延迟近似处理，使得本发明可以利用快速傅利叶变换实现卷积，避免了传统时域叠加法中逐点计算累加的繁琐步骤，有效降低了计算量。特别是对高程为常数的场景，如果采用斜距递推方法，可以进一步简化干涉合成孔径雷达回波仿真中较耗时的斜距计算过程，降低回波仿真的计算量。因此，采用本发明可以高效地完成回波仿真，实现在现有计算能力下，短时间内完成回波仿真；同时精度能够满足干涉合成孔径雷达系统仿真任务的要求。

附图说明

图1为本发明提供的合成孔径雷达三维陆地场景回波仿真方法流程示意图；

图2为斜距递推方法的递推顺序图；

图3为仿真实验中使用的数字高程模型数据灰度图；

图4为仿真实验得到的回波信号处理后的干涉相位图；

图5为常见陆地场景的地形参数P1～P6的取值；

图6为利用本发明一具体实施方式进行点目标仿真成像结果的性能分析；

图7为利用本发明一具体实施方式进行面目标仿真成像结果的性能分析。

图1为本发明提供的合成孔径雷达三维陆地场景回波仿真方法流程示意图。整个流程分为三步。第一步，几何建模；这一步包括数字高程模型数据插值和平面拟合两个步骤。经过几何建模，三维陆地场景被抽象为彼此相邻的小面单元，形成三维陆地场景的几何模型。第二步，电磁建模；在三维陆地场景的几何模型基础上，利用几何关系和电磁散射特性，计算每个小面单元的雷达散射截面积，形成三维陆地场景的电磁模型。第三步，回波生成，这一步包括单条回波生成和多条回波并行计算两步，通过这一步，可以得到三维陆地场景的回波。

具体实施方式

图2为斜距递推方法的递推顺序图。在每个斜距计算单元内，首先计算中心小面单元(m_c，n_c)的斜距。然后，分四步计算其它小面单元的斜距：(一)计算与中心小面同行，位于中心小面单元右边的小面单元的斜距；(二)计算与中心小面单元同行，位于中心小面单元左边的小面单元的斜距；(三)按列计算中心小面单元所在行上方的小面单元的斜距；(四)按列计算中心小面单元所在行下方的小面单元的斜距。

图3～图4是在实验室进行的回波仿真实验输入的数字高程模型数据和得到的干涉相位灰度图。仿真使用的集群计算机由1个主节点DELL PE2650及15个运算节点PE1750组成，每个节点配有2个Xeon3.0GHz的CPU、4GB内存。使用的操作系统为红帽3.0系统，仿真程序用C/C++开发。用本发明的仿真方法进行仿真，仿真时间为672.426秒。如果用传统时域叠加法计算，理论计算如果用1个CPU需要595.720天，即使也用32个CPU进行并行计算，仍需要23.270天。实验说明本算法有效提高了干涉合成孔径雷达三维陆地场景的回波仿真速度。

图3为仿真实验中使用的数字高程模型数据灰度图。该数字高程模型数据的点数为1105×1055，描述的三维陆地场景范围为2.5公里×2.0公里，高程最小值为0.11米，高程最大值为224.6米，高程平均值为99.5米，高程标准差为56.6米，坡度最大值为0.802696，坡度平均值为0.162504。

图4为仿真实验得到的回波信号处理后的干涉相位图。基于该干涉相位图进行数字高程模型数据重建结果的高程误差最小值为-3.7788米，高程误差最大值为4.1136米，高程误差均值为0.0556米，高程误差标准差为0.5795米，高程误差绝对值均值为0.4386米，高程误差绝对值标准差0.3829米。从上述结果可以看出，本发明的算法能够完成干涉合成孔径雷达系统仿真任务。

图5为常见陆地场景下的地形参数P1～P6的取值。采用堪萨斯大学的F.T.Ulaby和密执安大学的M.C.Dobson建立的经验模型。图中列出了在土壤与岩石、树林、草地、灌木、矮植被、路面、城区和湖泊7种陆地场景下的地形参数P1～P6的取值。

图6为利用本发明一具体实施方式进行点目标仿真成像结果的性能分析。使用传统的时域叠加法和本发明提供的回波仿真方法对点目标进行仿真。对比了点目标的相位误差，方位向和距离向的分辨率、积分旁瓣比、峰值旁瓣比。图中的数据表明，利用本发明一具体实施方式进行的仿真结果与传统时域叠加法的仿真结果几乎没有差别，说明本发明的回波仿真方法具有较高的精度。

图7为利用本发明一具体实施方式进行面目标仿真成像结果的性能分析。利用坡度为-10°、0°、-10°的场景进行成像结果，统计相干系数估计值与理论值的对比结果。一般要求仿真的相干系数估计值与理论值相差在10％以内，本次仿真的相干系数估计值与理论值误差相差都在1.3％以内，说明本发明算法精度能够满足干涉合成孔径雷达系统仿真要求。

Claims (3)

1.一种干涉合成孔径雷达三维陆地场景回波仿真方法，利用已知的三维陆地场景的数字高程模型数据，完成以下步骤：

第一步：几何建模

第(一)步，数字高程模型数据插值

若数字高程模型数据的采样间距为Δx_origin，干涉合成孔径雷达的地距分辨率为Δρ，需要进行的插值次数N_interp为

$N_{\mathrm{int\; erp}}=\mathrm{ceil}\left({\log }_2\left(\frac{9\mathrm{\Δ}{x}_{\mathrm{origin}}}{\mathrm{\Δ\ρ}}\right)\right)$

其中，ceil(x)为：

ceil(x)＝k，k-1＜x≤k，k为整数

对数字高程模型数据进行N_interp次的插值，每次插值的算法如下：

设本次插值前数字高程模型数据的点数为I×J，且任一点的位置序号为(i₀，j₀)，1≤i₀≤I，1≤j₀≤J，高程为g(i₀，j₀)；本次插值后数字高程模型数据点数为(2I-1)×(2J-1)，且任一点位置序号为(i′，j′)，1≤i′≤2I-1，1≤j′≤2J-1，高程为f(i′，j′)；用下述步骤计算f(i′，j′)：

对i′和j′都为奇数，用下式进行赋值：

$f(i^{\′},j^{\′})=g(\frac{i^{\′}+1}{2},\frac{j^{\′}+1}{2})$

对i′和j′都为偶数，用下式进行赋值：

$f(i^{\′},j^{\′})=\frac{1}{4}[f(i^{\′}-1,j^{\′}-1)+f(i^{\′}+1,j^{\′}-1)+f(i^{\′}-1,j^{\′}+1)+f(i^{\′}+1,j^{\′}+1)]+$

$2^{-H/2}\sqrt{1-2^{2H-2}}{\left|\right|\mathrm{\Δx}\left|\right|}^H\mathrm{\σ}{G}_N$

对i′和j′只有一个为偶数，用下式进行赋值：

$f(i^{\′},j^{\′})=\frac{1}{4}[f(i^{\′},j^{\′}-1)+f(i^{\′}-1,j^{\′})+f(i^{\′}+1,j^{\′})+f(i^{\′},j^{\′}+1)]+$

$2^{-H/2}\sqrt{1-2^{2H-2}}{\left|\right|\mathrm{\Δx}\left|\right|}^H\mathrm{\σ}{G}_N$

上述两个公式中的σ为本次插值前数字高程模型数据的标准差，H是三维陆地场景的频谱指数，C_N为N(0，1)分布的随机变量，Δx为本次插值前数字高程模型数据的采样间距；

第(二)步，平面拟合

对经过N_interp次插值后数字高程模型数据，以3×3的点阵为基本单位进行划分，形成若干个小面单元；

对任意一个小面单元，设该小面单元内任一点的位置序号为(i_s，j_s)，1≤i_s≤3，1≤j_s≤3，高程为f_s(i_s，j_s)，记：

$X=\left[\begin{array}{ccc}1& 0& 0\\ 1& 0& 1\\ 1& 0& 2\\ 1& 1& 0\\ 1& 1& 1\\ 1& 1& 2\\ 1& 2& 0\\ 1& 2& 1\\ 1& 2& 2\end{array}\right],Y=\left[\begin{array}{c}{f}_s\left(\mathrm{1,1}\right)\\ f_s\left(\mathrm{1,2}\right)\\ f_s\left(\mathrm{1,3}\right)\\ f_s\left(\mathrm{2,1}\right)\\ f_s\left(\mathrm{2,2}\right)\\ f_s\left(\mathrm{2,3}\right)\\ f_s\left(\mathrm{3,1}\right)\\ f_s\left(\mathrm{3,2}\right)\\ f_s\left(\mathrm{3,3}\right)\end{array}\right],$

用下式计算：

$\mathrm{\β}=\left[\begin{array}{c}{\mathrm{\β}}_1\\ {\mathrm{\β}}_2\\ {\mathrm{\β}}_3\end{array}\right]={\left(X^{\′}X\right)}^{-1}{X}^{\′}Y$

可求得该小面单元的法矢量

高程z_s＝β₁；

第二步，电磁建模

设三维陆地场景的中心在地心固连坐标系下的坐标为(x_c，y_c，z_c)，经过第一步得到小面单元的数量为M×N，小面单元距离向间隔为ρ_x，方位向间隔为ρ_y；

对任一小面单元，设其法矢量为n₀(m，n)，位置序号为(m，n)，1≤m≤M，1≤n≤N，则该小面单元在地心固连坐标系下的坐标(x_t(m，n)，y_t(m，n)，z_t(m，n))为(x_c+(m-M/2)ρ_x，y_c+(n-N/2)ρ_y，z_c+z_s(m，n))，进行如下步骤的计算：

首先，计算入射角θ_n(m，n)，计算公式为：

${\mathrm{\θ}}_n(m,n)=\frac{n_0(m,n)\·n_i(m,n)}{\left|n_0(m,n)\right|\left|n_i(m,n)\right|}$

其中，n_i(m，n)为波束方向矢量，设本次仿真的中心时刻干涉合成孔径雷达天线相位中心在地心固连坐标系下的坐标为(x_a，y_a，z_a)，则

$n_i(m,n)=\left[\begin{array}{c}{x}_t(m,n)-x_a\\ y_t(m,n)-y_a\\ z_t(m,n)-z_a\end{array}\right]$

然后，采用下式计算后向散射系数σ⁰(m，n)：

σ⁰(m，n)＝P₁+P₂ exp[-P₃θ_n(m，n)]+P₄cos(P₅θ_n(m，n)+P₆)

其中，P₁～P₆为与地形有关的参数；

最后，采用下式计算雷达散射截面积σ_F(m，n)：

σ_F(m，n)＝σ⁰(m，n)ρ_xρ_y；

第三步：回波生成

设干涉合成孔径雷达的发射信号为p(τ)，τ为快时间；

第(一)步，生成单条回波

对每个慢时间l，进行如下步骤的计算：

第(1)步，计算系统传输函数

首先，计算每个小面单元的斜距；

然后，确定每个小面单元在慢时间l时的增益总和，采用下述方法：

对位置序号为(m，n)，斜距为R(l，m，n)小面单元在慢时间l时的增益总和G(l，m，n)为

$G(l,m,n)=\frac{\sqrt{P_t{G}_t{G}_r}}{2\sqrt{\mathrm{\π}}\mathrm{\λR}{(l,m,n)}^2}$

其中，P_t为干涉合成孔径雷达的发射功率，G_t为发射天线增益，G_r为接收天线增益，λ为发射信号波长；

接着，确定每个小面单元的冲激响应，采用下述方法：

对位置序号为(m，n)的小面单元，计算冲激信号的幅度A(l，m，n)和相位

为：

A(l，m，n)＝σ_F(m，n)G(l，m，n)

上式中，f₀为干涉合成孔径雷达的中心频率，c为光速；

采用延迟近似处理计算冲激响应的延迟Δτ(l，m，n)为：

$\mathrm{\Δ\τ}(l,m,n)=\mathrm{round}\left(\frac{2R(l,m,n)}{\mathrm{c\ΔT}}\right)\mathrm{\ΔT},$

round(x)＝k，k-0.5≤x＜k+0.5，k为整数

上述ΔT为干涉合成孔径雷达接收信号的采样时间间隔，

可得到冲激响应为

其中，δ(τ)为单位冲激函数；

最后，可以得到在慢时间l的系统传输函数：

$h(l,\mathrm{\τ})=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{h}_s(l,m,n,\mathrm{\τ})$

第(2)步，计算回波信号

首先，对发射信号p(τ)进行快速傅利叶变换，得到发射信号的频谱P(f)；对系统传输函数h(l，τ)进行快速傅利叶变换，得到系统传输函数的频谱H(l，f)；

然后，把发射信号频谱P(f)和系统传输函数的频谱H(l，f)相乘，得到回波的频谱S(l，f)；

最后，对S(l，f)进行快速傅利叶反变换，得到在慢时间l的单条回波s(l，τ)；

第(二)步，并行计算多条回波

首先计算在所有慢时间的单条回波，再把每个慢时间生成的单条回波按时间顺序合成，得到多条回波。

2.根据权利要求1所述的干涉合成孔径雷达三维陆地场景回波仿真方法，其特征在于计算计算每个小面单元的斜距，采用下述方法：

即利用下式计算位置序号为(m，n)的小面单元的斜距：

$R(l,m,n)=\sqrt{{[x_l-x_t(m,n)]}^2+{[y_l-y_t(m,n)]}^2+{[z_l-z_t(m,n)]}^2}$

上式中，(x_l，y_l，z_l)为在慢时间l时，天线相位中心在地心固连坐标系下的坐标。

3.根据权利要求1所述的干涉合成孔径雷达三维陆地场景回波仿真方法，其特征在于，对高程为常数z₀时的三维陆地场景，采用下述方法计算斜距：

将相邻的M_r行×N_r列个小面单元划分为一个斜距计算单元，其中：

$M_r=\mathrm{floor}\left(\frac{15000\mathrm{\λ}}{{\mathrm{\ρ}}_x}\right)$

$N_r=\mathrm{floor}\left(\frac{15000\mathrm{\λ}}{{\mathrm{\ρ}}_y}\right);$

对每个斜距计算单元，进行下述步骤：

首先，设当前斜距计算单元的中心小面单元位置序号为(m_c，n_c)，计算其慢时间l时的斜距R(l，m_c，n_c)：

$R(l,m_c,n_c)=\sqrt{{[x_l-x_t(m_c,n_c)]}^2+{[y_l-y_t(m_c,n_c)]}^2+{[z_l-z_t(m_c,n_c)]}^2};$

然后，用下列式子计算递推系数A至E：

$A=-\frac{x_l-x_t(m_c,n_c)}{R(l,m_c,n_c)}{\mathrm{\ρ}}_x,B=-\frac{y_l-y_t(m_c,n_c)}{R(l,m_c,n_c)}{\mathrm{\ρ}}_y$

$C=(-\frac{{(x_l-x_t(m_c,n_c))}^2}{R^3(l,m_c,n_c)}+\frac{1}{R(l,m_c,n_c)}){{\mathrm{\ρ}}_x}^2,D=(-\frac{{(y_l-y_t(m_c,n_c))}^2}{R^3(l,m_c,n_c)}+\frac{1}{R(l,m_c,n_c)}){{\mathrm{\ρ}}_y}^2$

$E=-\frac{\mathrm{AB}}{R(l,m_c,n_c)}$

最后，递推计算当前斜距计算单元内其他小面单元的斜距：

依次令q＝1，2…(N_r-n_c-1)，用下式计算位置序号为(m_c，n_c+q+1)的小面单元在慢时间l时的斜距：

R(l，m_c，n_c+q+1)＝R(l，m_c，n_c+q)+qD+(D/2+B)

依次令q＝-1，-2…1-n_c，用下式计算位置序号为(m_c，n_c+q-1)的小面单元在慢时间l时的斜距：

R(l，m_c，n_c+q-1)＝R(l，m_c，n_c+q)-qD+(D/2+B)

依次令q＝1-n_c，2-n_c…(N_r-n_c-1)，p＝1，2…(M_r-m_c-1)，用下式计算位置序号为(m_c+p+1，n_c+q)的小面单元在慢时间l时的斜距：

R(l，m_c+p+1，n_c+q)＝R(l，m_c+p，n_c+q)+pC+(qE+C/2+A)

依次令q＝1-n_c，2-n_c…(N_r-n_c-1)，p＝-1，-2…1-m_c，用下式计算位置序号为(m_c+p-1，n_c+q)的小面单元在慢时间l时的斜距：

R(l，m_c+p-1，n_c+q)＝R(l，m_c+p，n_c+q)-pC+(qE+C/2-A)。

Images