CN103018740A - 一种基于曲面投影的InSAR成像方法 - Google Patents

Abstract

本发明公开了一种基于曲面投影的InSAR成像方法，该方法通过在低精度的高程曲面上进行反向投影成像，将不同雷达天线的测量数据投影到以低精度高程曲面建立的InSAR成像空间中，提高了InSAR成像精度和干涉相位保持精度，同时提取干涉相位以后，陡变地形区域的干涉条纹稀疏，后续相位解缠的难度降低，相位解缠的精度提高，使得高程反演的精度提高，为InSAR高质量干涉相位生成提供了一种新方法。

Description

一种基于曲面投影的InSAR成像方法

技术领域

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

背景技术

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

高质量的干涉相位是InSAR获取高精度地形数字高程模型（DEM）的基础，随着高程测量精度的不断提高，InSAR数据处理对干涉相位质量的要求也越来越高，一方面需要提高InSAR成像的干涉相位保持精度，另一方面需要稀疏化干涉条纹，以降低后续InSAR数据处理中相位解缠的难度，提高高程反演精度。传统的InSAR成像处理通常采用距离多普勒（RD）算法、变尺度（CS）算法、距离徙动（RM）算法等获取SAR图像，这些成像算法由于基于参考点成像，同时对平台运动轨迹作了近似直线处理，降低了InSAR成像精度和干涉相位保持精度。在提取干涉相位以后，干涉相位的表现形式是干涉条纹，通常陡变地形区域的干涉条纹密集，使得后续相位解缠的难度增加，同时相位解缠的精度也大大降低，导致最后的高程反演精度也明显降低，因此，如何稀疏化干涉条纹并提高相位解缠精度一直是InSAR技术研究的重点。

反向投影算法是一种基于时域相干处理的成像算法，其基本思想是通过计算成像区域内每一散射点到合成孔径长度内雷达平台天线相位中心之间的双程时延，然后将对应的时域回波信号进行相干累加，从而恢复出每个散射点的散射系数信息。曲面投影算法是反向投影算法的改进算法，曲面投影算法通过在低精度的高程曲面上进行反向投影成像，精确计算投影曲面上每个散射点与一个合成孔径长度内雷达平台天线相位中心的斜距信息，实现每个散射点的精确聚焦。将曲面投影算法应用于InSAR数据成像，不仅可以减小由于参考点近似和平台运动轨迹误差对成像的影响，提高InSAR成像精度和干涉相位保持精度，而且可以稀疏化陡变地形区域的干涉条纹，使得提取的干涉相位易于相位解缠，同时提高相位解缠精度，从而提高高程反演精度。

发明内容

为了得到高精度的InSAR干涉相位，同时稀疏化陡变地形区域的干涉条纹，本发明提出了一种基于曲面投影的InSAR成像方法，该方法通过在低精度的高程曲面上进行反向投影成像，将不同雷达天线的测量数据投影到以低精度高程曲面建立的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 onA novel 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、图像配准

图像配准是指将不同时间、不同传感器(成像设备)或不同条件下(天候、照度、摄像位置和角度等)对同一个地区获取的两幅或多幅图像进行地理坐标的匹配，使不同图像中的相同目标或者特征点位于图像中同一个位置的过程。

本发明中的图像配准特指将不同天线SAR图像空间中对应的相同目标进行匹配的过程，详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。

定义12、天线相位中心

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

定义13、范数

设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维离散信号，向量XLP范数为 ${\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回波信号与系统仿真研究”，张剑琦，哈尔滨工业大学硕士论文。

定义15、相位解缠

相位解缠是指将干涉相位由主值或相位差值恢复到真实值的过程，详见文献“星载合成孔径雷达干涉成像”，王超等编著，科学出版社出版。

本发明提供的一种基于曲面投影的InSAR成像方法，它包括以下几个步骤：

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

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

该坐标基方向为方位向；定义在地平面内，并与InSAR成像空间的第一个坐标基

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

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

该坐标基方位为高度向。

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

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

其中和

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

和

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

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

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

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

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

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

雷达系统距离向采样点数，记做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_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}(i;a,r)=\exp \{j\&CenterDot;4\&CenterDot;\mathrm{\&pi;}\&CenterDot;{\stackrel{\&OverBar;}{R}}_m(i;a,r)/\mathrm{\&lambda;}\},$ 得到散射点P(a,r)在慢时刻i应补偿的相位因子

其中j为虚数单位（即-1开根），为步骤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成像方法，该方法不仅具有良好的SAR图像干涉相位保持特性，还具有稀疏化陡变地形区域干涉条纹的优点，解决了高精度InSAR成像中，传统InSAR成像方法干涉相位提取精度差和陡变地形区域干涉条纹密集的缺点，提高了高程反演精度。

本发明的优点在于针对高精度InSAR成像系统，通过在低精度的高程曲面上进行反向投影成像，将不同雷达天线的测量数据投影到以低精度高程曲面建立的InSAR成像空间中，提高了InSAR成像精度和干涉相位保持精度，同时提取干涉相位以后，陡变地形区域的干涉条纹稀疏，后续相位解缠的难度降低，相位解缠的精度提高，使得高程反演的精度提高。

附图说明

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

图2为InSAR成像空间中实际测绘地形与低精度高程在距离向－高度向的二维坐标系中的关系图；

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

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

为雷达平台速度矢量，

为雷达平台主天线第i个慢时刻的天线相位中心矢量，为雷达平台副天线第i个慢时刻的天线相位中心矢量，为雷达平台主天线i时刻处散射点P(a,r)的距离史，

为雷达平台副天线i时刻处散射点P(a,r)的距离史，i为雷达平台距离采样点P(a,r)前后半个合成孔径内的慢时刻，平面ABCD为实际仿真坡面，平面ABC′D′为InSAR成像空间中低精度高程组成的坡面。

具体实施方式

本发明主要采用仿真实验的方法进行验证，所有步骤、结论都在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，雷达平台一个合成孔径长度内的慢时刻个数B_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-250=1，若a+250>2048，则a+250＝2048。采用公式计算得到雷达平台主天线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)的距离史

对应的距离门，计算公式为 ${\stackrel{\&OverBar;}{I}}_s(i;a,r)=1+\mathrm{round}(({\stackrel{\&OverBar;}{R}}_s(i;a,r)-9992)/0.5),$ 其中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、插值重采样后数据相干求和

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

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

为步骤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高质量干涉相位生成提供了一种新方法。

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，场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门位置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_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(·)为近似取整函数， ${\stackrel{\&OverBar;}{R}}_m(i;a,r)$ 为步骤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}(i;a,r)=\exp \{j\&CenterDot;4\&CenterDot;\mathrm{\&pi;}\&CenterDot;{\stackrel{\&OverBar;}{R}}_m(i;a,r)/\mathrm{\&lambda;}\},$ 得到散射点P(a,r)在慢时刻i应补偿的相位因子其中j为虚数单位（即-1开根），

为步骤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