CN113608218B - 一种基于后向投影原理的频域干涉相位稀疏重构方法 - Google Patents

Abstract

本发明公开了一种基于后向投影原理的频域干涉相位稀疏重构方法，它是通过对接收到原始回波数据进行距离向上的脉冲压缩，然后在选定的成像网格中，计算网格上每个像素点到在合成孔径长度内平台天线的斜距，进而计算并精确补偿每个像素点的回波延迟相位，然后再进行相干累积操作，提高信噪比，最终实现观测区域的SAR成像。本发明通过利用BP算法精确计算每个像素点的补偿相位，采用稀疏重构算法极大地提高算法的运算效率，同时有效提高了InSAR复图像干涉相位的精度，并且可以自由地选择成像空间。

Description

一种基于后向投影原理的频域干涉相位稀疏重构方法

技术领域

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

背景技术

作为一种工作在微波波段的有源雷达，合成孔径雷达(Synthetic ApertureRadar,SAR)具有全天时、全天候的成像能力，即无论是白天或黑夜、晴天还是雷雨风雪天气，都可以随时随地成像，克服了光学和红外系统不能在晚上和复杂天气条件进行成像的缺点。传统的SAR成像一般只具有二维成像分辨率，在一些起伏比较大的地方比如陡峭的山峰、峡谷以及城市中矗立挺拔的高楼时，传统SAR成像存在的失真(阴影遮挡效应、空间模糊、顶底倒置等)导致空间的一些重要信息(比如高度)丢失，所以能对目标进行三维成像是非常有必要的，为了适应这种需求，目前常见的三维成像技术有圆周SAR(Circular SAR)三维成像、层析SAR(Tomography SAR)三维成像、线阵SAR(Array SAR，ASAR)三维成像。

BP算法是一种在时域处理的雷达成像算法，其优势在于成像上高聚焦深度，得到的SAR复图像相位精度更好和能够自由地选择成像空间。但是其缺点在于时域计算时，需要精确计算每个像素点的补偿相位，成像网格划得越细，消耗时间越长。随着计算机硬件和软件系统发展，计算机的运算能力逐步提高，从而使得BP算法逐渐应用到SAR的多种模式中，此处将BP算法应用于宽测绘带InSAR干涉相位稀疏重构。

基于RD原理的传统InSAR成像方法虽然整体的运算复杂度小且成像算法的流程和基本原理简单，但是在距离徙动严重时，RD原理的校正误差较大，特别是在宽测绘场景下，这个误差更加明显，从而使得干涉相位图精度不高，而且对带限信号的采样频率需要大于等于带限信号最高频率的两倍，即奈奎斯特采样定理。然而在实际SAR成像中，若仍然采用奈奎斯特采样定理，海量原始回波数据的采样、存储和传输等操作对硬件设备和信号传输速率要求极高。如果一个信号是稀疏的或者是可压缩的，那么这个信号就能以低于Nyquist采样定理要求的采样率精确的重构出该信号，这就是稀疏重构的基本思想。针对稀疏重构理论用于SAR成像，目前的重构算法大概可以分为以下几类：贪婪追踪算法、凸松弛算法、贝叶斯框架算法、组合算法。

在信号处理领域中，稀疏信号处理是近年来很好地利用欠采样数据重构出原始信号的方式之一。稀疏信号处理理论其中关键是重构信号具有稀疏性，将稀疏信号处理理论引入到扫描式SAR模式中，来解决回波数据在方位向上的欠采样问题，目前的针对SAR复图像稀疏重构方法，可以较好地重构出SAR复图像幅度信息，但是重构出的SAR复图像相位信息损失严重，无法保证干涉相位的准确提取，进而影响干涉相位的精度，因此该思路不适用于宽测绘InSAR中。本发明提出了一种基于后向投影原理的频域干涉相位稀疏重构方法。

发明内容

为了解决扫描式SAR模式下回波数据欠采样问题和提高InSAR成像的干涉相位精度，本发明提出的基于后向投影原理的频域干涉相位稀疏重构方法，该方法先对接收到原始回波数据进行距离向上的脉冲压缩，然后在选定的成像网格中，计算网格上每个像素点到在合成孔径长度内平台天线的斜距，从而计算并精确补偿每个像素点的回波延迟相位，然后再进行相干累积操作，提高信噪比，最终实现观测区域的SAR成像。该算法通过利用BP算法精确计算每个像素点的补偿相位，有效提高了InSAR复图像干涉相位的精度，并且可以自由地选择成像空间。

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

定义1、合成孔径雷达(SAR)

合成孔径雷达是将雷达固定于载荷运动平台上，结合平台的运动以合成等效阵列以实现阵列向的分辨率，再利用雷达波束向回波延时实现距离一维成像，从而实现对观测目标二维成像的一种合成孔径雷达技术，详见文献“合成孔径雷达成像原理，杨建宇等编著，电子科技大学出版社”。

定义2、标准合成孔径雷达回波数据距离向脉冲压缩方法

标准合成孔径雷达回波数据距离向脉冲压缩是指利用合成孔径雷达发射信号参数，采用匹配滤波技术对合成孔径雷达的距离向信号进行信号聚焦成像的过程。详见文献“雷达成像技术”，保铮，邢孟道，王彤，电子工业出版社，2005。

定义3、范数

设X是复数域

上线性空间，其中/>

表示复数域，若它满足如下性质：||X||≥0，且当X||＝0仅有X＝0；||aX||＝|a|||X||，其中a为任意常数；||X₁+X₂||≤||X₁||+||X₂||，则称||X||为X空间上的范数(norm)，其中X₁和X₂为X空间上的任意两个值。对于定义1中的N×1维离散信号向量X＝[x₁,x₂,…,x_N]^T，向量X的LP范数表达式为/>

其中x_i为向量X的第i个元素，∑|·|表示绝对值求和运算符号，向量X的L1范数表达式为

向量X的L2范数表达式为/>

向量X的L0范数表达式为

且x_i≠0。详见文献“矩阵理论”，黄廷祝等编著，高等教育出版社出版。

定义4、方位向、距离向

将雷达平台运动的方向叫做方位向，将垂直于方位向的方向叫做距离向。详见文献“合成孔径雷达成像原理，杨建宇等编著，电子科技大学出版社”。

定义5、压缩感知稀疏重构理论

如果一个信号是稀疏的或可压缩的，那么该信号就可以用远低于奈奎斯特采样定理所要求的采样率来无失真的重构出该信号。如果信号稀疏，并且测量矩阵满足不相干和RIP属性，使用压缩感知恢复的信号稀疏重建可以通过解决以下最优化问题来实现：

其中，α是估计信号，y是测量信号，Θ是测量矩阵，ε是噪声门限。详见文献“阵列三维合成孔径雷达稀疏成像技术研究”韦顺军，2013。

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

合成孔径雷达原始回波仿真方法是指基于合成孔径雷达成像原理仿真出一定系统参数条件下具有合成孔径雷达回波信号特性的原始信号的方法，详见文献“张朋，合成孔径雷达回波信号仿真研究，西北工业大学博士论文，2004”。

定义8、线阵SAR的快时刻和慢时刻

线阵SAR运动平台飞过一个方位向合成孔径长度所需要的时间称为慢时间，雷达系统以一定时间长度的重复周期发射接收脉冲，因此慢时间可以表示为一个以脉冲重复周期为步长的离散化时间变量，其中每一个脉冲重复周期离散时间变量值为一个慢时刻。快时刻是指在一个脉冲重复周期内，距离向采样回波信号的时间间隔变量。详见文献“合成孔径雷达成像原理”，皮一鸣等编著，电子科技大学出版社出版。

定义9、经典软阈值迭代方法((Iterative Soft Thresholding Algorithm,ISTA)

经典软阈值迭代方法(ISTA)是求解范数稀疏性正则化的反问题目标函数，其要求所要求解的原始反问题是信号的稀疏表示。详见文献“K.Bredies,D.A.Lorenz.Linearconvergence of iterative soft-thresholding,2008”。

定义10、传统理论成像分辨率

线阵SAR成像传统理论分辨率是指利用经典匹配滤波理论成像算法得到线阵SAR系统在距离向、方位向和切航迹向的成像分辨率。对于收发共用天线，线阵SAR距离向的分辨率记为ρ_r，近似表达式为

其中C为电磁波在空气中传播的速度，B_r为线阵SAR发射信号带宽；方位向分辨率记为ρ_a，近似表达为/>

其中D_a为天线在方位向的真实孔径。详见文献“合成孔径雷达成像原理，杨建宇等编著，电子科技大学出版社”。

定义11、距离向重采样

重采样分为上采样和下采样，本发明中使用的是上采样，需要对距离向脉冲压缩后的信号进行插值，使计算时延能准确对应距离压缩回波数据，常采用的重采样方式为sinc插值。详见网站“https://blog.csdn.net/weixin_33918114/article/details/89697178”。

定义12、距离徙动校正

距离徙动校正是补偿InSAR回波数据的距离徙动量，使得同一目标点的回波信号位于同一个距离门。通常距离徙动校正分为距离走动校正和距离弯曲校正。详见文献“合成孔径雷达成像原理，杨建宇等编著，电子科技大学出版社”。

定义13、线性插值方法

线性插值方法是指插值函数为一次多项式的插值方式，其在插值节点上的插值误差为零。线性插值的几何意义即为概述图中利用过A点和B点的直线来近似表示原函数。线性插值可以用来近似代替原函数，也可以用来计算得到查表过程中表中没有的数值。详见文献“数值计算方法”，蔡锁章等编著，国防工业出版社。

本发明提出的一种基于后向投影原理的频域干涉相位稀疏重构方法，它包括如下步骤：

步骤1、初始化SAR系统参数：

初始化SAR系统参数包括：平台速度矢量记为

v_r为方位向速度；平台发射天线和接受天线初始位置矢量，记为/>

和/>

两天线的基线长度，记为L；雷达发射信号载频，记为f_c；脉冲重复时间记为PRI；雷达系统的脉冲重复频率为PRF；雷达发射信号带宽记为B_r；电磁波在空气中的传播速度记为c；距离向快时间记为t，t＝t₁,…,t_T，T为距离向快时刻总数，方位向慢时刻记为τ，τ＝τ₁,…,τ_K，K为方位向慢时刻总数，τ＝[τ₁,τ₂,…τ_K]^T为方位向慢时刻矢量；

上述参数均为SAR系统标准参数，其中雷达信号载频f_c，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽B_r；平台速度矢量

在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；

步骤2、划分SAR的成像场景空间：

以雷达波束照射场区域地平面的单位向量所构成的平面直角坐标系作为InSAR的成像场景目标空间Ω；初始化距离向成像场空间长度为L_x，方位向成像场空间长度为L_y；将成像场景目标空间Ω进行划分为大小相等的平面单元格，成像场景平面在距离向、方位向单元格数分别为M_x，M_y；

采用公式

计算得到距离向、方位向的单元格大小，分别记为d_x和d_y；

采用公式P_uv＝(u,v)＝(m_xd_x,m_yd_y)，计算得到划分后的成像场景空间Ω中的成像平面中第m_x个距离向单元格第m_y个方位单元格所对应的元素的位置(u,v)，记为P_uv，其中m_x＝1，…,M_x，m_y＝1，…,M_y；

步骤3、InSAR的距离向脉冲压缩：

在实际InSAR成像中，原始回波数据由数据接收机提供，在第τ个方位向慢时刻和第t个距离向快时刻中InSAR的原始回波数据记为s_r(τ,t)；

采用定义2中定义的标准合成孔径雷达距离向脉冲压缩方法对s_r(τ,t)进行距离向脉冲压缩，得到距离向压缩后的InSAR数据，记为s_AC(τ,t)；

步骤4、网格中所有目标点的斜距计算：

步骤4.1、采用公式

计算得到平台发射天线在τ时刻的位置，记为P_T(τ)，其中τ＝τ₁,…,τ_K为步骤1中方位向慢时刻，K为步骤1中方位向慢时刻总数，/>

为步骤1中平台发射天线的初始位置，/>

为步骤1中平台速度矢量；

采用公式

计算得到平台接收天线在τ时刻的位置，记为P_R(τ)，其中/>

为步骤1中平台接收天线的初始位置。

步骤4.2、采用公式R_uv(τ)＝||P_T(τ)-P_uv||₂+||P_R(τ)-P_uv||₂，计算得到每个目标点在τ时刻与平台天线的斜距，记为R_uv(τ)，其中τ＝τ₁,…,τ_K为步骤1中方位向慢时刻，P_uv＝(u,v)为步骤2中划分后的成像场景空间Ω中任意一像素点的位置坐标，u＝m_xd_x，m_x＝1，…,M_x，M_x为步骤2中距离向单元格总数，d_x为距离向单元格大小，v＝m_yd_y，m_y＝1，…,M_y，M_y为步骤2中方位向单元格总数，d_y为方位向单元格大小，K为步骤1中方位向慢时刻总数，P_T(τ)为步骤4中平台发射天线在τ时刻的位置，P_R(τ)为步骤4中平台接收天线在τ时刻的位置，||·||₂为定义3中定义的范数；

步骤5、构造观测矩阵；

步骤5.1、采用公式

计算得到每个目标点与平台天线的最近斜距，记为R_uv(i)，其中R_uv(τ)为步骤4.2中每个目标点在τ时刻与平台天线的斜距，τ＝τ₁,…,τ_K为步骤1中方位向慢时刻，min()为求取目标函数的最小值运算符号；

采用公式

计算得到h_i，其中w_a(·)表示原始回波信号在方位向上的包络，τ为步骤1中的方位向慢时刻矢量，x_i＝(i-1-K/2)d_y为平台的辐射区间内的第i个方位时刻位置，i＝1,2,…,K，K为步骤1中方位向慢时刻总数，d_y为步骤2中成像时的方位向间隔，v_r为步骤1中的平台方位向速度；

采用h_i构建矩阵Φ＝[h₁,…h_i,…,h_K]；

步骤5.2、采用公式

计算得到载波波数，记为K₀，其中/>

为载波波长，c为步骤1中电磁波在空气中的传播速度，f_c为步骤1中载波频率；

采用公式φ_i＝K₀R_uv(τ_i)，计算得到每个目标点相位的第i个元素，记为φ_i，其中i＝1,…,K，K为步骤1中方位向慢时刻总数，R_uv(τ_i)为步骤5中τ_i时刻网格中点目标与平台天线的斜距；

采用公式P＝diag{exp(jφ_i)}，计算得到天线的相位对角化形式矩阵，记为P；

步骤6、生成干涉相位图；

采用定义9中经典的软阈值迭代方法求解稀疏重构问题

估计得到InSAR干涉相位图，记为/>

其中/>

为使目标函数f(x)最小值时的x的值，s_AC为步骤3中距离向压缩后的回波信号，Φ为步骤5.1中构建的矩阵，P为步骤5.2中的相位对角化矩阵，Ψ为傅里叶变换基，λ为正则化约束参数，||·||₂和||·||_p为定义3中定义的范数；

至此，我们得到全场景InSAR的干涉相位结果，整个重构方法结束。

本发明的创新点是：针对RD算法在距离徙动严重时，校正误差较大，特别是在宽测绘场景下，这个误差更加明显，从而使得干涉相位图精度不高。本发明在定义9中的ISTA的基础上，结合BP算法，对成像区域进行网格划分，计算网格上每个像素点到在合成孔径长度内平台天线的斜距，从而计算并精确补偿每个像素点的回波延迟相位，然后再进行相干累积操作，提高信噪比，最终实现观测区域的SAR成像，有效提高了InSAR成像干涉图的精度。

本发明优点在于针对RD算法在距离徙动严重时，校正误差较大，特别是在宽测绘场景下，这个误差更加明显，从而使得干涉相位图精度不高。本发明在定义9中的ISTA的基础上，结合BP算法，对成像区域进行网格划分，计算网格上每个像素点到在合成孔径长度内平台天线的斜距，从而计算并精确补偿每个像素点的回波延迟相位，然后再进行相干累积操作，提高信噪比，最终实现观测区域的SAR成像。该算法采用稀疏重构算法极大地提高算法的运算效率，同时有效提高了InSAR成像干涉图的精度。

附图说明

图1为本发明流程图；

图2为系统参数表；

具体实施方式

本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-2017b上验证正确。具体实施步骤如下：

步骤1、初始化SAR系统参数：

初始化SAR系统参数包括：平台速度矢量记为

v_r＝50m/s为方位向速度；平台发射天线和接受天线初始位置矢量，记为P_T0＝[0,3000,0]m和P_R0＝[0.707,0,3000.707]m；两天线的基线长度，记为L＝1m；雷达发射信号载频记为f_c＝35GHz；脉冲重复时间记为PRI＝5ms；雷达系统的脉冲重复频率为PRF＝200Hz；雷达发射信号带宽记为B_r＝0.4GHz；电磁波在空气中的传播速度记为c＝299792458m/s；距离向快时间记为t，t＝t₁,…,t_T，T＝800为距离向快时刻总数，方位向慢时刻记为τ，τ＝τ₁,…,τ_K＝-3.7575，-3.7525，…，3.7475，3.7525，K＝1503为方位向慢时刻总数，τ＝[τ₁,τ₂,…τ_K]^T＝[-3.7575，-3.7525，…，3.7475，3.7525]^T为方位向慢时刻矢量；

上述参数均为SAR系统标准参数，其中雷达信号载频f_c，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽B_r；平台速度矢量

在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；

步骤2、划分SAR的成像场景空间：

以雷达波束照射场区域地平面的单位向量所构成的平面直角坐标系作为InSAR的成像场景目标空间Ω；初始化距离向成像场空间长度为L_x＝300m，方位向成像场空间长度为L_y＝300m；将成像场景目标空间Ω进行划分为大小相等的平面单元格，成像场景平面在距离向、方位向单元格数分别为M_x＝801，M_y＝801；

采用公式

计算得到距离向、方位向的单元格大小，分别记为d_x和d_y；

采用公式P_uv＝(u,v)＝(m_xd_x,m_yd_y)＝(0.3747m_x,0.3747m_y)得到划分后的成像场景空间Ω中的成像平面中第m_x个距离向单元格第m_y个方位向单元格所对应的元素的位置(u,v)，记为P_uv，其中m_x＝1，…,M_x，m_y＝1，…,M_y；

步骤3、InSAR的距离向脉冲压缩：

在实际InSAR成像中，原始回波数据由数据接收机提供，在第τ个方位向慢时刻和第t个距离向快时刻中InSAR的原始回波数据记为s_r(τ,t)；

采用定义2中定义的标准合成孔径雷达距离向脉冲压缩方法对s_r(τ,t)进行距离向脉冲压缩，得到距离向压缩后的InSAR数据，记为s_AC(τ,t)；

步骤4、网格中所有目标点的斜距计算：

步骤4.1、采用公式

计算得到平台发射天线在τ时刻的位置，记为P_T(τ)，其中τ＝τ₁,…,τ_K＝-3.7575，-3.7525，…，3.7475，3.7525为步骤1中方位向慢时刻，K＝1503为步骤1中方位向慢时刻总数，/>

为步骤1中平台发射天线的初始位置，

为步骤1中平台速度矢量；

采用公式

计算得到平台接收天线在τ时刻的位置，记为P_R(τ)，其中/>

为步骤1中平台接收天线的初始位置。

步骤4.2、采用公式R_uv(τ)＝||P_T(τ)-P_uv||₂+||P_R(τ)-P_uv||₂，计算得到每个目标点在τ时刻与平台天线的斜距，记为R_uv(τ)，其中τ＝τ₁,…,τ_K＝-3.7575，-3.7525，…，3.7475，3.7525为步骤1中方位向慢时刻，P_uv＝(u,v)为步骤2中划分后的成像场景空间Ω中任意一像素点的位置坐标，u＝m_xd_x，m_x＝1，…,M_x，M_x＝801为步骤2中距离向单元格总数，d_x＝0.3747m为距离向单元格大小，v＝m_yd_y，m_y＝1，…,M_y，M_y＝801为步骤2中方位向单元格总数，d_y＝0.3747m为方位向单元格大小，K＝1503为步骤1中方位向慢时刻总数，P_T(τ)为步骤4中平台发射天线在τ时刻的位置，P_R(τ)为步骤4中平台接收天线在τ时刻的位置，||·||₂为定义3中定义的范数；

步骤5、构造观测矩阵；

步骤5.1、采用公式

计算得到每个目标点与平台天线的最近斜距，记为R_uv(i)，其中R_uv(τ)为步骤4.2中每个目标点在τ时刻与平台天线的斜距，τ＝τ₁,…,τ_K＝-3.7575，-3.7525，…，3.7475，3.7525为步骤1中方位向慢时刻，min(·)为求取目标函数的最小值运算符号；

采用公式

计算得到h_i，其中w_a(·)表示原始回波信号在方位向上的包络，τ为步骤1中的方位向慢时刻矢量，x_i＝(i-1-K/2)d_y为平台的辐射区间内的第i个方位时刻位置，i＝1,2,…,K，K＝1503为步骤1中方位向慢时刻总数，d_y＝0.3747m为步骤2中成像时的方位向间隔，v_r＝50m/s为步骤1中的平台方位向速度；

采用h_i构建矩阵Φ＝[h₁,…h_i,…,h_K]；

步骤5.2、采用公式

计算得到载波波数，记为K₀，其中/>

为载波波长，c＝299792458m/s为步骤1中电磁波在空气中的传播速度，f_c＝35GHz为步骤1中载波频率；

采用公式φ_i＝K₀R_uv(τ_i)，计算得到每个目标点相位的第i个元素，记为φ_i，其中i＝1,…,K，K＝1503为步骤1中方位向慢时刻总数，R_uv(τ_i)为步骤5中τ_i时刻网格中点目标与平台天线的斜距；

采用公式P＝diag{exp(jφ_i)}，计算得到天线的相位对角化形式矩阵，记为P；

步骤6、生成干涉相位图；

采用定义9中的软阈值迭代算法求解稀疏重构问题

估计得到InSAR干涉相位图，记为/>

其中/>

为使目标函数f(x)最小值时的x的值，s_AC为步骤3中距离向压缩后的回波信号，Φ为步骤5.1中构建的矩阵，P为步骤5.2中的相位对角化矩阵，Ψ为傅里叶变换基，λ为正则化约束参数，||·||₂和||·||_p为定义3中定义的范数；

至此，我们得到全场景InSAR的干涉相位结果，整个重构方法结束。

经过计算机仿真及实测数据结果证明，本发明首先距离向脉冲压缩后的回波信号，并利用BP算法成像得到天线的相位图，并利用迭代软阈值算法ISTA进行成像。和RD算法比较，本发明通过利用BP算法的精确补偿相位和稀疏重构有效提高了InSAR成像干涉图的精度并降低了运算量，极大地提升了InSAR成像的运算效率。

Claims (1)

1.一种基于后向投影原理的频域干涉相位稀疏重构方法，其特征是它包括如下步骤：

步骤1、初始化SAR系统参数：

初始化SAR系统参数包括：平台速度矢量记为

v_r为方位向速度；平台发射天线和接受天线初始位置矢量，记为/>

和/>

两天线的基线长度，记为L；雷达发射信号载频，记为f_c；脉冲重复时间记为PRI；雷达系统的脉冲重复频率为PRF；雷达发射信号带宽记为B_r；电磁波在空气中的传播速度记为c；距离向快时间记为t，t＝t₁,…,t_T，T为距离向快时刻总数，方位向慢时刻记为τ，τ＝τ₁,…,τ_K，K为方位向慢时刻总数，τ＝[τ₁,τ₂,…τ_K]^T为方位向慢时刻矢量；

上述参数均为SAR系统标准参数，其中雷达信号载频f_c，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽B_r；平台速度矢量

在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；

步骤2、划分SAR的成像场景空间：

以雷达波束照射场区域地平面的单位向量所构成的平面直角坐标系作为InSAR的成像场景目标空间Ω；初始化距离向成像场空间长度为L_x，方位向成像场空间长度为L_y；将成像场景目标空间Ω进行划分为大小相等的平面单元格，成像场景平面在距离向、方位向单元格数分别为M_x，M_y；

采用公式

计算得到距离向、方位向的单元格大小，分别记为d_x和d_y；

采用公式P_uv＝(u,v)＝(m_xd_x,m_yd_y)，计算得到划分后的成像场景空间Ω中的成像平面中第m_x个距离向单元格第m_y个方位单元格所对应的元素的位置(u,v)，记为P_uv，其中m_x＝1，…,M_x，m_y＝1，…,M_y；

步骤3、InSAR的距离向脉冲压缩：

在实际InSAR成像中，原始回波数据由数据接收机提供，在第τ个方位向慢时刻和第t个距离向快时刻中InSAR的原始回波数据记为s_r(τ,t)；

采用标准合成孔径雷达距离向脉冲压缩方法对s_r(τ,t)进行距离向脉冲压缩，得到距离向压缩后的InSAR数据，记为s_AC(τ,t)；

步骤4、网格中所有目标点的斜距计算：

步骤4.1、采用公式

计算得到平台发射天线在τ时刻的位置，记为P_T(τ)，其中τ＝τ₁,…,τ_K为步骤1中方位向慢时刻，K为步骤1中方位向慢时刻总数，/>

为步骤1中平台发射天线的初始位置，/>

为步骤1中平台速度矢量；

采用公式

计算得到平台接收天线在τ时刻的位置，记为P_R(τ)，其中

为步骤1中平台接收天线的初始位置

步骤4.2、采用公式R_uv(τ)＝||P_T(τ)-P_uv||₂+||P_R(τ)-P_uv||₂，计算得到每个目标点在τ时刻与平台天线的斜距，记为R_uv(τ)，其中τ＝τ₁,…,τ_K为步骤1中方位向慢时刻，P_uv＝(u,v)为步骤2中划分后的成像场景空间Ω中任意一像素点的位置坐标，u＝m_xd_x，m_x＝1，…,M_x，M_x为步骤2中距离向单元格总数，d_x为距离向单元格大小，v＝m_yd_y，m_y＝1，…,M_y，M_y为步骤2中方位向单元格总数，d_y为方位向单元格大小，K为步骤1中方位向慢时刻总数，P_T(τ)为步骤4中平台发射天线在τ时刻的位置，P_R(τ)为步骤4中平台接收天线在τ时刻的位置，||·||₂为范数；

步骤5、构造观测矩阵；

步骤5.1、采用公式

计算得到每个目标点与平台天线的最近斜距，记为R_uv(i)，其中R_uv(τ)为步骤4.2中每个目标点在τ时刻与平台天线的斜距，τ＝τ₁,…,τ_K为步骤1中方位向慢时刻，min(·)为求取目标函数的最小值运算符号；

采用公式

计算得到h_i，其中w_a(·)表示原始回波信号在方位向上的包络，τ为步骤1中的方位向慢时刻矢量，x_i＝(i-1-K/2)d_y为平台的辐射区间内的第i个方位时刻位置，i＝1,2,…,K，K为步骤1中方位向慢时刻总数，d_y为步骤2中成像时的方位向间隔，v_r为步骤1中的平台方位向速度；

采用h_i构建矩阵Φ＝[h₁,…h_i,…,h_K]；

步骤5.2、采用公式

计算得到载波波数，记为K₀，其中/>

为载波波长，c为步骤1中电磁波在空气中的传播速度，f_c为步骤1中载波频率；

采用公式φ_i＝K₀R_uv(τ_i)，计算得到每个目标点相位的第i个元素，记为φ_i，其中i＝1,…,K，K为步骤1中方位向慢时刻总数，R_uv(τ_i)为步骤5中τ_i时刻网格中点目标与平台天线的斜距；

采用公式P＝diag{exp(jφ_i)}，计算得到天线的相位对角化形式矩阵，记为P；

步骤6、生成干涉相位图；

采用经典的软阈值迭代方法求解稀疏重构问题

估计得到InSAR干涉相位图，记为/>

其中/>

为使目标函数f(x)最小值时的x的值，s_AC为步骤3中距离向压缩后的回波信号，Φ为步骤5.1中构建的矩阵，P为步骤5.2中的相位对角化矩阵，Ψ为傅里叶变换基，λ为正则化约束参数，||·||₂和||·||_p为范数；

至此，得到全场景InSAR的干涉相位结果，整个重构方法结束。

Images