CN103018741A - 一种基于后向投影的InSAR成像去平地一体化方法 - Google Patents

Abstract

本发明公开了一种基于反向投影的InSAR成像去平地一体化方法，它是通过将后向投影算法引入InSAR成像，以场景地平面为成像空间，主、副天线分别单独进行后向投影成像，然后提取主、副复图像的干涉相位，不仅提高了干涉相位的保相精度，而且不需要进行去平地效应处理，简化了InSAR数据处理流程，提高了地形高程反演精度。通过残差点统计和主、副复图像相干性测试，可以验证本发明所提供方法得到的干涉相位质量明显好于通过传统InSAR成像方法得到的干涉相位质量。

Description

一种基于后向投影的InSAR成像去平地一体化方法

技术领域

本发明属于雷达技术领域，它特别涉及合成孔径雷达（SAR）成像技术领域。

背景技术

合成孔径雷达（SAR）是一种高分辨率的微波成像***。合成孔径雷达利用大时宽带宽信号实现距离向高分辨率，依靠雷达和目标之间的相对运动来合成虚拟阵列获取方位向高分辨率，可以对照射场景进行二维成像。合成孔径雷达干涉测量（InSAR）是一般SAR功能的延伸和扩展，是利用两个或者多个位置不同的天线观测同一个目标场景，根据目标到不同天线的斜距差获得测量数据的干涉相位，再通过平台与地面观测场景的几何关系反演出地面场景的数字高程信息的技术。由于具有全天时、全天候的特点，InSAR已经成为当前提取大面积地表三维图像和地形高程变化信息的一项重要遥感技术，在地形测绘、自然灾害监测和自然资源调查等领域发挥越来越大的作用。

高保相的干涉相位是InSAR获取高精度地形数字高程模型（简称DEM）的基础，随着高程测量精度的不断提高，InSAR数据处理对干涉相位保相精度的要求也越来越高。传统的InSAR单视成像处理通常采用距离多普勒（简称RD）和变尺度（简称CS）等频域算法获取单视复图像，这些成像算法由于基于参考点成像，同时对平台运动轨迹作了近似直线处理，降低了InSAR成像精度和干涉相位保持精度。另一方面，传统的InSAR数据处理需要在进行单视成像处理以后，进行去平地效应操作，通常是采用基于轨道参数等先验信息或基于离散傅里叶变换频移的去平地效应方法，这些方法均会引入相位误差，同时增加了InSAR数据处理的复杂度。

反向投影（简称BP）算法是一种基于时域相干处理的成像算法，其基本思想是通过计算成像区域内每一采样点到合成孔径长度内雷达天线相位中心之间的双程时延，然后将对应的时域回波信号进行相干累加，从而恢复出每个采样点的散射系数信息。后向投影算法不基于参考点成像，能够实现每个采样点的精确聚焦，同时后向投影算法的雷达天线相位中心物理意义清晰，能够精确计算每个方位时刻的回波延时相位，补偿平台抖动引入的相位误差，将其应用于InSAR单视成像处理，能够显著地提高干涉相位的保相精度，提高干涉相位质量。此外，通过在后向投影算法方位聚焦时，主天线补偿主天线的多普勒相位因子，副天线补偿副天线的多普勒相位因子，即主、副天线分别单独进行后向投影成像，得到主、副复图像，然后提取的干涉相位不需要进行去平地效应，即主、副复图像提取的干涉相位只包含地形高程信息，没有平地相位。

发明内容

为了得到高保相的去除平地效应的InSAR干涉相位，本发明提出了一种基于后向投影的InSAR成像方法，该方法实现在反向投影（BP）成像的同时实现去除平地效应，提高了InSAR成像精度和干涉相位保持精度，最后提取的干涉相位只包含地形高程信息，简化了InSAR数据处理流程。

为了方便描述本发明的内容，首先作以下术语定义：

定义1、合成孔径雷达干涉测量（简称InSAR）

合成孔径雷达干涉测量（InSAR）指利用两个或者两个以上的SAR数据中的相位信息进行相干处理，结合雷达参数和雷达几何位置信息反演地表三维及其变化信息的遥感技术，详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。

定义2、数字高程模型（简称DEM）

数字高程模型（Digital Elevation Model，DEM）是指利用一组有序数值阵列形式表示地表或地面高程的一种实体地面模型。本发明中DEM表示成一系列地面点的平面坐标X、Y和高程坐标Z组成的数据阵列。对于一个地面区域D，地形DEM表示为

DEM={D_i |(x_i,y_i,z_i),i∈D}

其中(x_i,y_i)是第i个地面像素点对应的平面坐标，z_i是对应的高程坐标。

定义3、雷达成像空间

雷达成像空间是指将场景空间中的散射点投影到距离向－方位向－高度向的三维空间坐标系，该空间由合成孔径雷达成像空间中的三个相互正交的坐标基确定。目前典型的合成孔径雷达成像空间包括距离向－方位向－高度向投影空间。本发明中用以下数学关系表示成像空间M：

$M=\left\{\stackrel{\‾}{P}(v,u)\right|\stackrel{\‾}{P}(v,u)=v\·{\stackrel{\‾}{\mathrm{\ζ}}}_v+u\·{\stackrel{\‾}{\mathrm{\ζ}}}_u+H(v,u)\·{\stackrel{\‾}{\mathrm{\ζ}}}_h,u,v\∈R\}$

其中和

表示构成成像空间M的三个相互正交的坐标基，分别表示距离向、方位向和高度向。为成像空间中的采样点向量，u,v分别表示该点的距离向坐标和方位向坐标，H(v,u)表示该点的高度向坐标，H(v,u)与距离向坐标u和方位向坐标v存在一一对应关系，R表示实数。

定义4、合成孔径雷达后向投影算法

后向投影算法是基于匹配滤波原理的合成孔径雷达成像算法，其主要通过相干累加实现合成孔径雷达数据的聚焦成像。详细内容可参考文献：“Research on Anovel fast back projection algorithm for strip map bistatic SAR imaging”，Huang Yulin等。

定义5、天线相位中心

天线相位中心是指雷达天线向外辐射信号的中心，本发明中天线相位中心指雷达平台天线的轨迹位置。

定义6、合成孔径雷达标准距离压缩方法

合成孔径雷达标准距离压缩方法是指利用合成孔径雷达发射参数，采用以下公式生成参考信号，并采用匹配滤波技术对合成孔径雷达的距离向信号进行滤波的过程。

$f\left(t\right)=\exp (j\·\mathrm{\π}\·\frac{B}{T_r}\·t^2),t\∈[-\frac{T_r}{2},\frac{T_r}{2}]$

其中，j为虚数单位（即-1开根），f(t)为距离压缩参考函数，B为雷达发射基带信号的信号带宽，T_r为雷达发射信号脉冲宽度，t为时间变量，取值范围从

到

详见文献“雷达成像技术”，保铮等编著，电子工业出版社出版。

定义7、辛格插值（sinc插值）方法

辛格插值方法是指对于一个带限信号，在满足采样定理的情况下，采用卷积核为sinc的函数h(x)，h(x)的长度即窗长为W。

$h\left(x\right)=\sin c\left(x\right)=\frac{\sin \left(\mathrm{\πx}\right)}{\mathrm{\πx}}$

进行对已离散的信号g_d(i)插值，得到插值后所要的信号

$g\left(x\right)=\underset{i}{\mathrm{\Σ}}{g}_d\left(i\right)\sin c(x-i)$

详见文献“合成孔径雷达成像----算法与实现”，Frank H.Wong等编著，电子工业出版社出版。

定义8、合成孔径与慢时刻

合成孔径是指对于测绘场景中的一个散射点从进入雷达波束照射范围至离开雷达波束照射范围的这段时间内，雷达波束中心所走过的长度。

雷达平台飞过一个合成孔径所需要的时间称为慢时间，雷达***以一定的重复周期T_r发射接收脉冲，因此慢时间可以表示为一个以重复周期T_r为步长的离散化时间变量，其中每一个离散时间变量值为一个慢时刻。

详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。

定义9、距离史与距离门

距离史是指不同天线相位中心到场景中散射点的距离组成的序列。

距离门是指对应距离史的回波数据在整个回波数据中的位置，详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。

定义10、合成孔径雷达成像场景参考点

合成孔径雷达成像场景参考点是指合成孔径雷达成像空间中的某个散射点，作为分析和处理场景中其他散射点的参照。

定义11、共轭相乘

共轭相乘是指将两个复数或复数矩阵对应数据的幅度相乘、相位相减的过程。

本发明中的共轭相乘指的是对两幅复图像对应像素点的数据，分别进行幅度相乘、相位相减的过程。

定义12、范数

设X是数域C上线性空间，称||·||为X上的范数(norm)，若它满足：1.正定性：||X||≥0，且||X||=0<=>X=0；2.齐次性：||aX||=a||X||；3.次可加性(三角不等式)：||X+Y||≤||X||+||Y||。若X=[x₁,x₂,…,x_N]^T为N×1维离散信号，向量X LP范数为 ${\left|\right|X\left|\right|}_P={\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\right|x^P\left|\right)}^{1/P},$ L1范数为 ${\left|\right|X\left|\right|}_1=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|x\right|,$ L2范数为 ${\left|\right|X\left|\right|}_2={\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left|x\right|}^2\right)}^{\frac{1}{2}}.$

定义13、近似取整函数round(x)

近似取整函数round(x)是指按四舍五入原则将实数x映射到相近的整数的函数。

定义14、合成孔径雷达原始回波仿真方法

合成孔径雷达原始回波仿真方法是指基于合成孔径雷达成像原理仿真出一定***参数条件下具有SAR回波信号特性的原始信号的方法，详细内容可参考文献：“InSAR回波信号与***仿真研究”，张剑琦，哈尔滨工业大学硕士论文。

本发明提供的一种基于后向投影的InSAR成像去平地一体化方法，它包括以下步骤：

步骤1、初始化InSAR成像***参数

InSAR成像空间由InSAR成像空间中的三个相互正交的坐标基确定，定义与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基，记做

该坐标基方向为方位向；定义在地平面内，并与InSAR成像空间的第一个坐标基垂直的单位向量作为InSAR成像空间的第二个坐标基，记做

该坐标基方向为距离向；定义垂直于地平面向上的单位向量作为InSAR成像空间的第三个坐标基，记做

该坐标基方位为高度向。

InSAR雷达平台包含两组天线，即主天线和副天线，两组天线之间的距离为基线长，记做B_l，主天线发射脉冲信号，经过T_d时间的延迟，主天线和副天线同时接收回波延迟信号。雷达平台主天线接收的回波数据，记做雷达平台副天线接收的回波数据，记做

其中

和

均为二维矩阵，二维矩阵的第一维均对应方位向，第二维均对应距离向，即二维矩阵

和

的行存储的是方位向数据，二维矩阵

和

的列存储的是距离向数据。

初始化InSAR成像***参数包括：雷达***工作的信号波长，记做λ，雷达平台主天线发射信号带宽，记做B，雷达平台主天线发射脉冲时宽，记做T_r，雷达平台接收***采样频率，记做F_s，雷达***脉冲重复频率，记做PRF，雷达平台一个合成孔径长度内的慢时刻个数，记做N_l，雷达平台速度矢量，记做

雷达平台主天线初始位置矢量，记做雷达平台副天线初始位置矢量，记做

场景参考点位置矢量，记做

雷达***距离向采样点数，记做N_r，雷达***方位向采样点数，记做N_a，场景距离向散射点间隔，记做d_r，场景方位向散射点间隔，记做d_a，场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离，记做R_mc，场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离，记做R_sc，场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门相同，距离门位置记做I_c。对于场景散射点P(a,r)，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数，场景散射点P(a,r)的低精度高程，记做H(a,r)。上述参数中，雷达***工作的信号波长λ，雷达平台主天线发射的信号带宽B，雷达平台主天线发射的脉冲时宽T_r，雷达平台接收***的采样频率F_s，雷达***的脉冲重复频率PRF，两组天线之间的基线长B_l以及接收***接收波门相对于发射信号发射波门的延迟T_d在InSAR雷达***设计过程中已经确定；雷达平台一个合成孔径长度内的慢时刻个数N_l，雷达平台速度矢量

雷达平台主天线初始位置矢量雷达平台副天线初始位置矢量

场景参考点位置矢量

雷达***距离向采样点数N_r，雷达***方位向采样点数N_a，场景距离向散射点间隔d_r，场景方位向散射点间隔d_a，场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离R_mc，场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离R_sc，场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门位置I_c以及场景散射点P(a,r)的低精度高程H(a,r)，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数，在InSAR雷达成像观测方案设计中已经确定。根据InSAR雷达***方案和InSAR雷达成像观测方案，以上基于后向投影的InSAR成像去平地一体化方法需要的初始化成像***参数均为已知。

步骤2：InSAR原始回波数据进行距离压缩

采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据

进行压缩，得到平台主天线距离压缩后数据，记做

采用传统的合成孔径雷达标准距离压缩方法对雷达平台副天线距离向回波数据进行压缩，得到平台副天线距离压缩后数据，记做

步骤3、计算InSAR成像空间中散射点的距离史

采用公式

计算得到雷达平台主天线第n个慢时刻的天线相位中心矢量

采用公式

计算得到雷达平台副天线第n个慢时刻的天线相位中心矢量

其中n表示第n个慢时刻，n=1,...,N_a，N_a为步骤1初始化得到的雷达***方位向采样点数，PRF为步骤1初始化得到的雷达***脉冲重复频率，为步骤1初始化得到的雷达平台速度矢量，为步骤1初始化得到的雷达平台主天线初始位置矢量，

为步骤1初始化得到的雷达平台副天线初始位置矢量。采用公式 $\stackrel{\&OverBar;}{P}(a,r)={\stackrel{\&OverBar;}{P}}_c+(r-1)\&CenterDot;d_r\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u+(a-1)\&CenterDot;d_a\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v+H(a,r)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_h$ 计算得到场景采样点P(a,r)的位置矢量，其中a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数，d_r为步骤1初始化得到的场景距离向散射点间隔，d_a为步骤1初始化得到的场景方位向散射点间隔，

为步骤1初始化得到的场景参考点位置矢量，H(a,r)为步骤1初始化得到的场景散射点P(a,r)的低精度高程，

为步骤1中定义的InSAR成像空间的第一个坐标基，

为步骤1中定义的InSAR成像空间的第二个坐标基，

为步骤1中定义的InSAR成像空间的第三个坐标基。i为雷达平台距离散射点P(a,r)前后半个合成孔径内的慢时刻，i的取值满足条件：i＝a-round(N₁/2),…,a+round(N_l/2)，当a-round(N_l/2)<1时，则我们取a-round(N_l/2)为1，当a+round(N_l/2)>N_a时，则我们取a+round(N_l/2)为N_a，其中round(·)为近似取整函数，N_l为步骤1初始化得到的雷达平台一个合成孔径长度内的慢时刻个数，N_a为步骤1初始化得到的雷达***方位向采样点数，即i只取雷达***采样范围内的雷达平台距离散射点P(a,r)前后半个合成孔径内的慢时刻。采用公式

计算得到雷达平台主天线i时刻处散射点P(a,r)的距离史

采用公式

计算得到雷达平台副天线i时刻处散射点P(a,r)的距离史其中||·||₂为L2范数。

步骤4、距离压缩后数据插值重采样

为雷达平台主天线i时刻处散射点P(a,r)的距离史

对应的距离门，计算公式为 ${\stackrel{\&OverBar;}{I}}_m(i;a,r)=I_c+\mathrm{round}(({\stackrel{\&OverBar;}{R}}_m(i;a,r)-R_{\mathrm{mc}})/d_r),$ 为雷达平台副天线i时刻处散射点P(a,r)的距离史

对应的距离门，计算公式为

其中round(·)为近似取整函数，

为步骤3中得到的雷达平台主天线i时刻处散射点P(a,r)的距离史，

为步骤3中得到的雷达平台副天线i时刻处散射点P(a,r)的距离史，R_mc为步骤1初始化得到的场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离，R_sc为步骤1初始化得到的场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离，I_c为步骤1初始化得到的场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门位置，d_r为步骤1初始化得到的场景距离向散射点间隔，i为步骤3中得到的雷达平台距离采样点P(a,r)前后半个合成孔径内的慢时刻，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数。

从步骤2中得到的平台主天线距离压缩后数据

的第i行中取第

个数据的前W₀个数据和后W₀个数据，采用标准辛格插值方法对这组数据进行插值，得到主天线插值重采样后的数据从步骤2中得到的平台副天线距离压缩后数据

的第i行中取第

个数据的前W₀个数据和后W₀个数据，采用标准辛格插值方法对这组数据进行插值，得到副天线插值重采样后的数据

其中W₀为标准辛格插值的半个窗长。

步骤5、插值重采样后数据相干求和

根据主天线补偿相位因子的计算公式

得到散射点P(a,r)在慢时刻i应补偿的主天线相位因子

根据副天线补偿相位因子的计算公式 ${\stackrel{\&OverBar;}{F}}_s(i;a,r)=\exp \{j\&CenterDot;4\&CenterDot;\mathrm{\&pi;}\&CenterDot;{\stackrel{\&OverBar;}{R}}_s(i;a,r)/\mathrm{\&lambda;}\},$ 得到散射点P(a,r)在慢时刻i应补偿的副天线相位因子

其中j为虚数单位（即-1开根），

为步骤3中得到的雷达平台主天线i时刻处散射点P(a,r)的距离史，

为步骤3中得到的雷达平台副天线i时刻处散射点P(a,r)的距离史，i为步骤3中得到的雷达平台距离采样点P(a,r)前后半个合成孔径内的慢时刻，λ为步骤1初始化得到的雷达***工作的信号波长，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数。

将步骤4中得到的主天线插值重采样后的数据

与散射点P(a,r)在慢时刻i应补偿的主天线相位因子

相乘，得到主天线相位补偿后的数据将步骤4中得到的副天线插值重采样后的数据

与散射点P(a,r)在慢时刻i应补偿的副天线相位因子

相乘，得到副天线相位补偿后的数据

把雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的所有慢时刻的主天线相位补偿后的数据

相加，得到雷达平台主天线场景散射点P(a,r)成像后的数据

即

把雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的所有慢时刻的副天线相位补偿后的数据相加，得到雷达平台副天线场景散射点P(a,r)成像后的数据

即 ${\stackrel{\&OverBar;}{B}}_s(a,r)=\underset{i}{\mathrm{\&Sigma;}}{\stackrel{\&OverBar;}{T}}_s(i;a,r).$

步骤6、干涉相位提取

将步骤5中得到的雷达平台主天线场景散射点P(a,r)成像后的数据

和雷达平台副天线场景散射点P(a,r)成像后的数据

进行共轭相乘，得到InSAR成像的干涉相位

完成基于后向投影的InSAR成像去平地一体化处理，其中a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数。

本发明的创新点在于提出了基于后向投影的InSAR成像去平地一体化方法，该方法首先具有高精度的干涉相位保相特性，同时其提取的干涉相位不需要进行去平地效应操作，简化了InSAR数据处理流程，提高了高程反演精度。

本发明的优点在于针对高精度InSAR成像***，通过在场景地平面上进行后向投影成像，实现在高精度InSAR成像的同时完成去平地效应操作，提取的干涉相位，不仅具有高保相精度，而且不需要进行去平地效应操作，简化了InSAR数据处理流程。

附图说明

图1为发明所提供方法的流程示意框图；

图2为本发明具体实施方式采用的干涉合成孔径雷达平台飞行几何关系及仿真圆锥场景图；

其中，横坐标（X轴）为方位向，纵坐标（Y轴）为距离向，垂直坐标（Z轴）为高度向，B_l为雷达平台主天线和副天线之间的基线，基线的两个端点分别为雷达平台主天线和副天线，P(a,r)为场景散射点，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为仿真场景方位向总的散射点数，s_a=512，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为仿真场景距离向总的散射点数，s_r=512，d_a为仿真场景方位向散射点间隔，d_a=0.3m，d_r为仿真场景距离向散射点间隔，d_r=0.5m，

为雷达平台速度矢量，

为雷达平台主天线第i个慢时刻的天线相位中心矢量，为雷达平台副天线第i个慢时刻的天线相位中心矢量，

为雷达平台主天线i时刻处散射点P(a,r)的距离史，为雷达平台副天线i时刻处散射点P(a,r)的距离史，i为雷达平台距离采样点P(a,r)前后半个合成孔径内的慢时刻。

具体实施方式

本发明主要采用仿真实验的方法进行验证，所有步骤、结论都在MATLAB7.0上验证正确。具体实施步骤如下：

步骤1、初始化InSAR成像***参数

选择与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基 ${\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v=\left[\begin{array}{ccc}0& 1& 0\end{array}\right],$ 该坐标基方向为方位向；选择InSAR成像空间的第二个坐标基 ${\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u=\left[\begin{array}{ccc}1& 0& 0\end{array}\right],$ 该坐标基方向为距离向；选择InSAR成像空间的第三个坐标基 ${\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_h=\left[\begin{array}{ccc}0& 0& 1\end{array}\right],$ 该坐标基方向为高度向。

初始化InSAR成像***参数包括：雷达***工作的信号波长λ＝0.03m，雷达平台主天线发射信号带宽B=150MHz，雷达平台主天线发射脉冲时宽T_r=1μs，雷达平台接收***采样频率F_s=300MHz，雷达***脉冲重复频率PRF=500Hz，雷达平台一个合成孔径长度内的慢时刻个数N_l=500，雷达平台速度矢量 $\stackrel{\&OverBar;}{V}=\left[\begin{array}{ccc}0& 150& 0\end{array}\right],$ 速度的单位为m/s，雷达平台主天线初始位置矢量 ${\stackrel{\&OverBar;}{P}}_m\left(0\right)=\left[\begin{array}{ccc}0& 0& 6000\end{array}\right],$ 雷达平台副天线初始位置矢量 ${\stackrel{\&OverBar;}{P}}_s\left(0\right)=\left[\begin{array}{ccc}10& 0& 6000\end{array}\right],$ 场景参考点位置矢量 ${\stackrel{\&OverBar;}{P}}_c=\left[\begin{array}{ccc}8000& 0& 0\end{array}\right],$ 位置的单位为m，雷达***距离向采样点数N_r=1024，雷达***方位向采样点数N_a=2048，场景距离向散射点间隔d_r=0.5m，场景方位向散射点间隔d_a=0.3m，场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离R_mc=10000m，场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离R_sc=9992，场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门I_c=1。

仿真场景为一立体圆锥，仿真场景方位向总的散射点数s_a=512，仿真场景距离向总的散射点数s_r=512，斜坡高度由0m到50m，对于场景散射点P(a,r)，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，低精度高程H(a,r)为真实仿真高度的0.8倍，仿真坡面场景如图3所示。

采用传统的合成孔径雷达原始回波仿真方法，生成雷达平台主天线仿真回波数据

采用传统的合成孔径雷达原始回波仿真方法，生成雷达平台副天线仿真回波数据

步骤2：InSAR原始回波数据进行距离压缩

采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据

进行压缩，得到平台主天线距离压缩后数据，记做

采用传统的合成孔径雷达标准距离压缩方法对雷达平台副天线距离向回波数据

进行压缩，得到平台副天线距离压缩后数据，记做

步骤3、计算InSAR成像空间中散射点的距离史

采用公式计算得到雷达平台主天线第n个慢时刻的天线相位中心矢量

采用公式计算得到雷达平台副天线第n个慢时刻的天线相位中心矢量

其中n表示第n个慢时刻，n=1,...,2048，

为步骤1初始化得到的雷达平台速度矢量，

为步骤1初始化得到的雷达平台主天线初始位置矢量，

为步骤1初始化得到的雷达平台副天线初始位置矢量。采用公式 $\stackrel{\&OverBar;}{P}(a,r)={\stackrel{\&OverBar;}{P}}_c+(r-1)\&CenterDot;0.5\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u+(a-1)\&CenterDot;0.3\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v+H(a,r)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_h$ 计算得到场景采样点P(a,r)的位置矢量，其中a表示散射点位于场景方位向的第a个位置，a=1，...,512，r表示散射点位于场景距离向的第r个位置，r=1，...,512，

为步骤1初始化得到的场景参考点位置矢量，H(a,r)为步骤1初始化得到的场景散射点P(a,r)的低精度高程，

为步骤1中定义的InSAR成像空间的第一个坐标基，

为步骤1中定义的InSAR成像空间的第二个坐标基，

为步骤1中定义的InSAR成像空间的第三个坐标基。i为雷达平台距离散射点P(a,r)前后半个合成孔径内的慢时刻，i的取值满足条件：i＝a-250,…,a+250，若a-250<1，则我们取a为251，若a+250>2048，则我们取a为1798。采用公式

计算得到雷达平台主天线i时刻处散射点P(a,r)的距离史

采用公式

计算得到雷达平台副天线i时刻处散射点P(a,r)的距离史

其中||·||₂为L2范数。

步骤4、距离压缩后数据插值重采样

为雷达平台主天线i时刻处散射点P(a,r)的距离史对应的距离门，计算公式为 ${\stackrel{\&OverBar;}{I}}_m(i;a,r)=1+\mathrm{round}(({\stackrel{\&OverBar;}{R}}_m(i;a,r)-10000)/0.5),$

为雷达平台副天线i时刻处散射点P(a,r)的距离史

对应的距离门，计算公式为

其中round(·)为近似取整函数，

为步骤3中得到的雷达平台主天线i时刻处散射点P(a,r)的距离史，

为步骤3中得到的雷达平台副天线i时刻处散射点P(a,r)的距离史，i为步骤3中得到的雷达平台距离采样点P(a,r)前后半个合成孔径内的慢时刻，a表示散射点位于场景方位向的第a个位置，a=1,...,512，r表示散射点位于场景距离向的第r个位置，r=1,...,512。

从步骤2中得到的平台主天线距离压缩后数据

的第i行中取第

个数据的前4个数据和后4个数据，采用标准辛格插值方法对这组数据进行插值，得到主天线插值重采样后的数据

从步骤2中得到的平台副天线距离压缩后数据

的第i行中取第

个数据的前4个数据和后4个数据，采用标准辛格插值方法对这组数据进行插值，得到副天线插值重采样后的数据

步骤5、插值重采样后数据相干求和

根据主天线补偿相位因子的计算公式

得到散射点P(a,r)在慢时刻i应补偿的主天线相位因子

根据副天线补偿相位因子的计算公式 ${\stackrel{\&OverBar;}{F}}_s(i;a,r)=\exp \{j\&CenterDot;4\&CenterDot;\mathrm{\&pi;}\&CenterDot;{\stackrel{\&OverBar;}{R}}_s(i;a,r)/0.03\},$ 得到散射点P(a,r)在慢时刻i应补偿的副天线相位因子

其中j为虚数单位（即-1开根），

为步骤3中得到的雷达平台主天线i时刻处散射点P(a,r)的距离史，为步骤3中得到的雷达平台副天线i时刻处散射点P(a,r)的距离史，i为步骤3中得到的雷达平台距离采样点P(a,r)前后半个合成孔径内的慢时刻，a表示散射点位于场景方位向的第a个位置，a=1,...,512，r表示散射点位于场景距离向的第r个位置，r=1,...,512。

将步骤4中得到的主天线插值重采样后的数据

与散射点P(a,r)在慢时刻i应补偿的主天线相位因子相乘，得到主天线相位补偿后的数据

将步骤4中得到的副天线插值重采样后的数据

与散射点P(a,r)在慢时刻i应补偿的副天线相位因子

相乘，得到副天线相位补偿后的数据

把雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的所有慢时刻的主天线相位补偿后的数据

相加，得到雷达平台主天线场景散射点P(a,r)成像后的数据即

把雷达平台副天线距离散射点P(a，r)前后半个合成孔径内的所有慢时刻的副天线相位补偿后的数据相加，得到雷达平台副天线场景散射点P(a,r)成像后的数据

即 ${\stackrel{\&OverBar;}{B}}_s(a,r)=\underset{i}{\mathrm{\&Sigma;}}{\stackrel{\&OverBar;}{T}}_s(i;a,r).$

步骤6、干涉相位提取

将步骤5中得到的雷达平台主天线场景散射点P(a,r)成像后的数据

和雷达平台副天线场景散射点P(a,r)成像后的数据

进行共轭相乘，得到InSAR成像的干涉相位完成基于后向投影的InSAR成像去平地一体化处理，其中a表示散射点位于场景方位向的第a个位置，a=1,...,512，r表示散射点位于场景距离向的第r个位置，r=1,...,512。

通过本发明具体实施方式的仿真及测试，本发明所提供的基于反向投影的InSAR成像去平地一体化方法，得到的干涉相位没有平地效应，通过残差点统计和主、副复图像相干性测试，验证本发明所提供方法得到的干涉相位质量明显好于通过传统InSAR成像方法得到的干涉相位质量。本发明通过将后向投影算法引入InSAR成像，以场景地平面为成像空间，主、副天线分别单独进行后向投影成像，然后提取主、副复图像的干涉相位，不仅提高了干涉相位的保相精度，而且不需要进行去平地效应处理，简化了InSAR数据处理流程，提高了地形高程反演精度。

Claims (1)

1.一种基于后向投影的InSAR成像去平地一体化方法，其特征是它包括以下步骤：

步骤1、初始化InSAR成像***参数

InSAR成像空间由InSAR成像空间中的三个相互正交的坐标基确定，定义与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基，记做该坐标基方向为方位向；定义在地平面内，并与InSAR成像空间的第一个坐标基

垂直的单位向量作为InSAR成像空间的第二个坐标基，记做

该坐标基方向为距离向；定义垂直于地平面向上的单位向量作为InSAR成像空间的第三个坐标基，记做

该坐标基方位为高度向；

InSAR雷达平台包含两组天线，即主天线和副天线，两组天线之间的距离为基线长，记做B_l，主天线发射脉冲信号，经过T_d时间的延迟，主天线和副天线同时接收回波延迟信号；雷达平台主天线接收的回波数据，记做

雷达平台副天线接收的回波数据，记做

其中

和

均为二维矩阵，二维矩阵的第一维均对应方位向，第二维均对应距离向，即二维矩阵

和

的行存储的是方位向数据，二维矩阵

和

的列存储的是距离向数据；

初始化InSAR成像***参数包括：雷达***工作的信号波长，记做λ，雷达平台主天线发射信号带宽，记做B，雷达平台主天线发射脉冲时宽，记做T_r，雷达平台接收***采样频率，记做F_s，雷达***脉冲重复频率，记做PRF，雷达平台一个合成孔径长度内的慢时刻个数，记做N_l，雷达平台速度矢量，记做

雷达平台主天线初始位置矢量，记做

雷达平台副天线初始位置矢量，记做

场景参考点位置矢量，记做

雷达***距离向采样点数，记做N_r，雷达***方位向采样点数，记做N_a，场景距离向散射点间隔，记做d_r，场景方位向散射点间隔，记做d_a，场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离，记做R_mc，场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离，记做R_sc，场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门相同，距离门位置记做I_c；对于场景散射点P(a,r)，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数，场景散射点P(a,r)的低精度高程，记做H(a,r)；上述参数中，雷达***工作的信号波长λ，雷达平台主天线发射的信号带宽B，雷达平台主天线发射的脉冲时宽T_r，雷达平台接收***的采样频率F_s，雷达***的脉冲重复频率PRF，两组天线之间的基线长B_l以及接收***接收波门相对于发射信号发射波门的延迟T_d在InSAR雷达***设计过程中已经确定；雷达平台一个合成孔径长度内的慢时刻个数N_l，雷达平台速度矢量

雷达平台主天线初始位置矢量

雷达平台副天线初始位置矢量

场景参考点位置矢量

雷达***距离向采样点数N_r，雷达***方位向采样点数N_a，场景距离向散射点间隔d_r，场景方位向散射点间隔d_a，场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离R_mc，场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离R_sc，场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门位置Ic以及场景散射点P(a,r)的低精度高程H(a,r)，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数，在InSAR雷达成像观测方案设计中已经确定；根据InSAR雷达***方案和InSAR雷达成像观测方案，以上基于后向投影的InSAR成像去平地一体化方法需要的初始化成像***参数均为已知；

步骤2：InSAR原始回波数据进行距离压缩

采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据

进行压缩，得到平台主天线距离压缩后数据，记做

采用传统的合成孔径雷达标准距离压缩方法对雷达平台副天线距离向回波数据

进行压缩，得到平台副天线距离压缩后数据，记做

步骤3、计算InSAR成像空间中散射点的距离史

采用公式计算得到雷达平台主天线第n个慢时刻的天线相位中心矢量

采用公式计算得到雷达平台副天线第n个慢时刻的天线相位中心矢量

其中n表示第n个慢时刻，n=1,...,Na，N_a为步骤1初始化得到的雷达***方位向采样点数，PRF为步骤1初始化得到的雷达***脉冲重复频率，

为步骤1初始化得到的雷达平台速度矢量，为步骤1初始化得到的雷达平台主天线初始位置矢量，

为步骤1初始化得到的雷达平台副天线初始位置矢量；采用公式 $\stackrel{\&OverBar;}{P}(a,r)={\stackrel{\&OverBar;}{P}}_c+(r-1)\&CenterDot;d_r\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_u+(a-1)\&CenterDot;d_a\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_v+H(a,r)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_h$ 计算得到场景采样点P(a,r)的位置矢量，其中a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数，d_r为步骤1初始化得到的场景距离向散射点间隔，d_a为步骤1初始化得到的场景方位向散射点间隔，

为步骤1初始化得到的场景参考点位置矢量，H(a,r)为步骤1初始化得到的场景散射点P(a,r)的低精度高程，为步骤1中定义的InSAR成像空间的第一个坐标基，为步骤1中定义的InSAR成像空间的第二个坐标基，为步骤1中定义的InSAR成像空间的第三个坐标基；i为雷达平台距离散射点P(a,r)前后半个合成孔径内的慢时刻，i的取值满足条件：i＝a-round(N_l/2),...,a+round(N_l/2)，当a-round(N_l/2)<1时，则我们取a-round(N_l/2)为1，当a+round(N_l/2)>N_a时，则我们取a+round(N_l/2)为N_a，其中round(·)为近似取整函数，N_l为步骤1初始化得到的雷达平台一个合成孔径长度内的慢时刻个数，N_a为步骤1初始化得到的雷达***方位向采样点数，即i只取雷达***采样范围内的雷达平台距离散射点P(a,r)前后半个合成孔径内的慢时刻；采用公式计算得到雷达平台主天线i时刻处散射点P(a,r)的距离史

采用公式

计算得到雷达平台副天线i时刻处散射点P(a,r)的距离史

其中||·||₂为L2范数；步骤4、距离压缩后数据插值重采样

为雷达平台主天线i时刻处散射点P(a,r)的距离史

对应的距离门，计算公式为 ${\stackrel{\&OverBar;}{I}}_m(i;a,r)=I_c+\mathrm{round}(({\stackrel{\&OverBar;}{R}}_m(i;a,r)-R_{\mathrm{mc}})/d_r),$

为雷达平台副天线i时刻处散射点P(a,r)的距离史

对应的距离门，计算公式为 ${\stackrel{\&OverBar;}{I}}_s(i;a,r)=I_c+\mathrm{round}(({\stackrel{\&OverBar;}{R}}_s(i;a,r)-R_{\mathrm{sc}})/d_r),$ 其中round(·)为近似取整函数，

为步骤3中得到的雷达平台主天线i时刻处散射点P(a,r)的距离史，

为步骤3中得到的雷达平台副天线i时刻处散射点P(a,r)的距离史，R_mc为步骤1初始化得到的场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离，R_sc为步骤1初始化得到的场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离，I_c为步骤1初始化得到的场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门位置，d_r为步骤1初始化得到的场景距离向散射点间隔，i为步骤3中得到的雷达平台距离采样点P(a,r)前后半个合成孔径内的慢时刻，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数；

从步骤2中得到的平台主天线距离压缩后数据

的第i行中取第

个数据的前W₀个数据和后W₀个数据，采用标准辛格插值方法对这组数据进行插值，得到主天线插值重采样后的数据

从步骤2中得到的平台副天线距离压缩后数据

的第i行中取第

个数据的前W₀个数据和后W₀个数据，采用标准辛格插值方法对这组数据进行插值，得到副天线插值重采样后的数据

其中W₀为标准辛格插值的半个窗长；

步骤5、插值重采样后数据相干求和

根据主天线补偿相位因子的计算公式 ${\stackrel{\&OverBar;}{F}}_m(i;a,r)=\exp \{j\&CenterDot;4\&CenterDot;\mathrm{\&pi;}\&CenterDot;{\stackrel{\&OverBar;}{R}}_m(i;a,r)/\mathrm{\&lambda;}\},$ 得到散射点P(a,r)在慢时刻i应补偿的主天线相位因子

根据副天线补偿相位因子的计算公式 ${\stackrel{\&OverBar;}{F}}_s(i;a,r)=\exp \{j\&CenterDot;4\&CenterDot;\mathrm{\&pi;}\&CenterDot;{\stackrel{\&OverBar;}{R}}_s(i;a,r)/\mathrm{\&lambda;}\},$ 得到散射点P(a,r)在慢时刻i应补偿的副天线相位因子

其中j为虚数单位（即-1开根），

为步骤3中得到的雷达平台主天线i时刻处散射点P(a,r)的距离史，

为步骤3中得到的雷达平台副天线i时刻处散射点P(a,r)的距离史，i为步骤3中得到的雷达平台距离采样点P(a,r)前后半个合成孔径内的慢时刻，λ为步骤1初始化得到的雷达***工作的信号波长，a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数；

将步骤4中得到的主天线插值重采样后的数据

与散射点P(a,r)在慢时刻i应补偿的主天线相位因子

相乘，得到主天线相位补偿后的数据

将步骤4中得到的副天线插值重采样后的数据与散射点P(a,r)在慢时刻i应补偿的副天线相位因子

相乘，得到副天线相位补偿后的数据

把雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的所有慢时刻的主天线相位补偿后的数据

相加，得到雷达平台主天线场景散射点P(a,r)成像后的数据即

把雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的所有慢时刻的副天线相位补偿后的数据相加，得到雷达平台副天线场景散射点P(a,r)成像后的数据

即 ${\stackrel{\&OverBar;}{B}}_s(a,r)=\underset{i}{\mathrm{\&Sigma;}}{\stackrel{\&OverBar;}{T}}_s(i;a,r);$

步骤6、干涉相位提取将步骤5中得到的雷达平台主天线场景散射点P(a,r)成像后的数据

和雷达平台副天线场景散射点P(a,r)成像后的数据

进行共轭相乘，得到InSAR成像的干涉相位

完成基于后向投影的InSAR成像去平地一体化处理，其中a表示散射点位于场景方位向的第a个位置，a=1,...,s_a，s_a为场景方位向总的散射点数，r表示散射点位于场景距离向的第r个位置，r=1,...,s_r，s_r为场景距离向总的散射点数。

Images