CN103941243A - 一种基于sar三维成像的自旋式飞行器测高方法 - Google Patents

Abstract

本发明提出了一种基于SAR三维成像的自旋式飞行器测高方法，它是将自旋式飞行器天线相位中心轨迹等效为天线面阵，通过等效面阵天线对目标场景三维成像，将飞行器测高问题转化为对三维成像场景中心点高度的求解，通过对目标场景高精度三维成像，完成对自旋式飞行器的测高。在成像过程中，本发明提供了三维BP算法，该算法能适用于等效天线面阵非均匀分布的情况。利用自旋式飞行器运动过程中形成的等效天线阵列，结合自旋式飞行器三维SAR的系统参数和观测场景目标参数，通过在时域对有散射点的场景单元格补偿时延相位，完成了对目标场景的三维成像。本发明具有对于变加速、非直线运动的自旋式飞行器能成像优点。

Description

一种基于SAR三维成像的自旋式飞行器测高方法

技术领域

本发明属于雷达系统中引导测高雷达的技术领域，它特别涉及合成孔径雷达（SAR）三维成像技术领域。

背景技术

自旋式飞行器平台具有运动速度快、非匀速直线运动等特点，尤其是由于飞行器自旋等复杂运动，造成雷达回波相位补偿困难，回波不相参，能量难以聚集，使得传统的测高手段失效。针对自旋式飞行器平台，配合三维数字高程图（DEM）对目标场景的三维成像，通过计算成像场景中心点高度以得出飞行器的测量高度，是一种测高新体制、新方法。

合成孔径雷达（Synthetic Aperture Radar,SAR）成像技术是一种常用的雷达成像技术。合成孔径雷达是一种主动式微波成像传感器，具有全天候、高精度的特点，可在光学能见度极差的气象条件下得到高分辨率雷达图像。合成孔径雷达利用大时宽带宽信号实现距离向高分辨率，依靠雷达和目标之间的相对运动来合成虚拟阵列获取方位向高分辨率，可以对照射场景进行二维成像。

阵列三维合成孔径雷达（3D-ASAR）是一般SAR功能的延伸和扩展，其基本原理是通过沿运动平台轨迹和切航迹向放置天线阵列合成二维虚拟面阵以获得沿-切航迹向二维分辨率，然后高度向使用脉冲相参积累的方法获得高度向分辨率。与传统的二维SAR成像技术相比，它可以克服传统二维SAR技术存在的地面阴影遮挡效应和空间模糊问题。

自旋式飞行器平台由于自身的复杂运动特性，雷达回波较难相参积累，不能用传统的成像算法进行成像。后向投影（简称BP）算法是一种基于时域相干处理的成像算法，其基本思想是通过计算成像区域内每一采样点到合成孔径长度内雷达天线相位中心之间的双程时延，然后将对应的时域回波信号进行相干累加，从而恢复出每个采样点的散射系数信息。后向投影算法不基于参考点成像，能够实现每个采样点的精确聚焦，同时后向投影算法的雷达天线相位中心物理意义清晰，能够精确计算每个方位时刻的回波延时相位，补偿平台抖动引入的相位误差。同时，在天线相位中心形成过程中，由于自旋式飞行器在变加速直线飞行时伴随飞行器自旋运动，阵列天线是非等间隔分布的，将三维BP算法应用于自旋式飞行器平台三维成像处理，克服了传统的三维RD等频域算法要求阵列天线等间隔采样的缺点，具有聚焦效果好、适用于非均匀或欠采样条件、对于变加速、非直线、自旋的运动轨迹能成像优点。

发明内容

为了得到对自旋式飞行器的高度测量值，本发明提出了一种基于SAR三维成像的自旋式飞行器测高方法，它是针对自旋式飞行器成像平台，将自旋式飞行器的天线相位中心运动轨迹看做非均匀分布的等效天线阵列，利用等效天线阵列对场景三维成像，成像场景中心点高度即为雷达天底点高度，结合飞行器的水平高度数据，得出飞行器飞行的相对高度。同时，针对等效阵列非均匀分布的特点，本发明采用了三维BP（后像投影）成像算法。与传统方法相比，本发明克服了传统的三维RD等频域算法要求阵列天线等间隔采样的缺点，具有聚焦效果好、适用于非均匀或欠采样条件、对于变加速、非直线、自旋的运动轨迹能成像优点。

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

定义1、阵列三维合成孔径雷达（简称ASAR）

阵列合成孔径雷达（Array Synthetic Aperture Radar）是指通过控制阵列天线在空间中运动轨迹形成虚拟二维面阵，获得观测目标二维分辨率，并结合脉冲压缩技术得到观测目标的第三维分辨率，完成对目标成像的技术。详见文献“A new airborne radar for3D imaging-imageformation using the ARTINO principle，J.Klare,A.R.Brenner,J.Ender。”

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

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

DEM={D_i|x_iO,y_iOz_iO),iO∈D}

其中(x_i0,y_i0,z_i0)是第i0个地面像素点对应的三维空间坐标。

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

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

定义4、一帧回波数据

在本发明中，一帧回波数据是指在一个雷达扫描周期内，雷达接收机所接收、采样并存储的在这一个雷达扫描周期内所有发射脉冲的回波数据。

定义5、飞行器天底点

飞行器上雷达正下方垂直于地面的交点。

定义6、合成孔径雷达相参积累与非相参积累

在雷达信号处理中脉冲的能量积累包括非相参积累和相参积累，前者指的是仅仅对数据的幅度（也可能是幅度的平方或者幅度对数）进行积累，而后者是指对复数据（即包含幅度和相位的数据）进行积累。

定义7、天线相位中心

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

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

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

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

其中，j为虚数单位（即-1开平方根），f(t)为距离压缩参考函数，B为雷达发射基带信号的信号带宽，T_r为雷达发射信号脉冲宽度，t₀为时间变量，取值范围从-T_r/2到T_r/2，详见文献“雷达成像技术，保铮等编著，电子工业出版社出版”

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

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

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

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

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

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

定义10、三维SAR慢时刻与快时刻

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

定义11、三维SAR观测场景目标空间

三维SAR观测场景目标空间是指现实空间中所有待观测场景目标散射点的位置集合。观测场景目标空间在不同空间坐标系下有不同的表示，但一旦坐标系确立以后其表示是唯一的。一般情况下为了方便成像，阵列SAR观测场景目标空间取为地面坐标系。本发明中用以下数学关系表示场景目标空间Ω：

$\mathrm{\Ω}=\left\{{\stackrel{\_}{P}}_s\right|{\stackrel{\_}{P}}_s=x\·{\stackrel{\_}{e}}_x+y\·{\stackrel{\_}{e}}_y+z\·{\stackrel{\_}{e}}_z;x,y,z\∈R\}$

其中表示构成观测场景目标空间Ω的地表正交坐标基，分别表示水平横向、水平纵向和垂直地表的高度向，为场景目标空间中一个分辨单元位置向量，x、y和z分别表示该分布单元的水平横向、水平纵向和高度向坐标，R表示实数域。

定义12、后像投影算法

BP算法起源于计算机断层扫描(Computer-Aided Tomography,CAT)技术,是一种精确的时域成像方法。详见文献“LARS M.H.ULANDER,HANS HELLSTEN.GUNNARSTENSTRO¨M.Synthetic-ApertureRadarProcessing Using FastFactorized Back-Projection.IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMSVOL.39,NO.3JULY2003”。

定义13、取向量或矩阵的指定元素

a、b、c、d为四个正整数，其中a<b，c<d。对于向量表示向量的第a个元素，表示取向量的第a到b个元素所构成的向量。对于二维矩阵A，A(a,b)表示矩阵A第a行b列的元素，A(a,：)表示取矩阵A第a行所有元素所构成的向量，A(：,c)表示取矩阵的第c列所有元素所构成的向量，A(a：b,：)表示取矩阵的a到b行所有元素所构成的矩阵，A(：,c：d)表示取矩阵的c到d列所有元素所构成的矩阵。详见文献“MATLAB宝典（第3版）”，陈杰等编著，电子工业出版社。”

定义14、脉冲压缩

脉冲压缩是一种现代雷达信号处理技术，简单来说就是雷达发射宽脉冲，然后再接收端“压缩”为窄脉冲，从而改善雷达的两种性能：作用距离和高度向分辨率。详见“皮亦鸣,杨建宇,付毓生,杨晓波.合成孔径雷达成像原理.第一版.电子科技大学出版社.2007.3”

定义15、范数

设X是数域上线性空间，表示复数域，若它满足如下性质：||X||≥0，且||X||=0仅有X=0；||aX||=|a|||X||，a为任意常数；||X₁+X₂||≤||X₁||+||X₂||，则称||X||为X空间上的范数(norm)，其中X₁和X₂为X空间上的任意两个值。对于N×1维离散信号向量X=[x₁,x₂...x_N]^T，向量X的LP范数表达式为其中x_i为向量X的第i个元素，Σ||表示绝对值求和符号，向量X的L1范数表达式为向量X的L2范数表达式为向量X的L0范数表达式为且x_i≠0。详见文献“矩阵理论，黄廷祝等编著，高等教育出版社出版”。

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

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

定义17、三维SAR传统理论成像分辨率

三维SAR传统理论成像分辨率是指利用经典匹配滤波理论成像算法得到阵列SAR系统在高度向、沿航迹向和切航迹向的成像分辨率。对于收发共用天线，阵列SAR高度向的分辨率记为ρ_r，近似表达式为其中C为光在空气中的传播速度，B_r为阵列SAR发射信号的带宽；沿航迹向的分辨率记为ρ_a，近似表达式为其中D_a为天线在沿航迹向的真实孔径；切航迹向的分辨率记为ρ_c，近似表达式为其中λ为阵列SAR雷达载频波长，R₀为阵列SAR平台到观测场景中心的参考斜距，L为阵列天线长度。详见参考文献“Shi,J.,et al.,APC trajectory design for One-Active linear-array three-dimensionalimaging SAR,IEEE Transactions on Geoscience and Remote Sensing,Vol.48,No.3,pp:1470-1486,2010”。

本发明提供的一种基于SAR三维成像的自旋式飞行器测高方法，它包括以下步骤：

步骤1、初始化自旋式飞行器SAR成像系统各参数

初始化自旋式飞行器SAR成像系统参数包括：飞行器水平高度，记做H，飞行器半径，记做r，飞行器飞行速度，记做v，最大加速度，记做a，自旋角频率，记做ω，由v、a、ω决定的飞行器速度矢量，记做雷达系统工作的信号波长，记做λ，雷达发射线性调频信号的载频，记做fc，线性调频信号的调频斜率，记做fdr，雷达平台天线发射信号带宽，记做B，雷达平台天线发射脉冲时宽，记做Tr，雷达平台接收系统采样频率，记做Fs，雷达发射脉冲重复频率，记做prf，光在空气中的传播速度，记做C，第t个慢时刻，记做t，t=1,2...T，T为慢时刻总数，雷达平台天线相位中心位置矩阵，记做APC，沿航迹向分辨率，记做ρ_a，切航迹向分辨率，记做ρ_b，高度向分辨率，记做ρ_r，接收波门相对发射信号延时Td，距离门位置记做Ic；上述参数均为阵列SAR系统标准参数，其中雷达发射线性调频信号的载频fc，雷达载频波长λ，雷达发射基带信号的信号带宽B，雷达发射信号脉冲宽度Tp，雷达发射信号调频斜率f_dr，雷达接收系统的采样频率Fs，雷达系统的脉冲重复频率prf，光在空气中的传播速度C，雷达系统的脉冲重复时间PRI和雷达接收系统接收波门相对于发射信号发散波门的延迟TD，均在三维SAR系统设计过程中已经确定；飞行器飞行高度H，飞行器半径r，飞行器速度v，加速度a，自旋角频率ω，平台速度矢量沿航迹向分辨率ρ_a，切航迹向分辨率ρ_b，高度向分辨率ρ_r，雷达平台天线相位中心位置矩阵APC，均在自旋式飞行器三维SAR观测方案设计中已经确定。根据三维SAR雷达系统方案和三维SAR雷达成像观测方案，以上基于三维成像的自旋式飞行器测高系统参数均为已知。

步骤2：初始化阵列三维SAR的观测场景目标空间

初始化阵列三维SAR的观测场景目标空间参数，包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标作为阵列三维SAR的观测场景目标空间Ω；将三维SAR观测场景目标空间Ω均匀划分成大小相等的立体单元网格(亦称为分辨单元)，立体单元网格在水平横向、水平纵向和高度向边长分别记为dx、dy和dz，立体单元网格的每一对水平横向、水平纵向所对应高度单元格中，只有一个有点目标散射体，立体单元网格大小选择为阵列SAR系统传统理论成像分辨率的二分之一，初始化数字高程模型（DEM）数据，其中，场景沿航迹总的散射点数，记做sx，场景切航迹向总的散射点数，记做sy，场景沿航迹向上的第x0个散射点序号，记做x0，x₀=1,2,...s_x，场景沿航迹向上的第y0个散射点序号，记做y₀，y₀=1,2,...s_y，沿航迹向散射点位置x₀与切航迹向散射点位置y₀所唯一确定的目标散射点的高度，记做z(x₀,y₀)。观测场景目标空间Ω中第m个散射点的坐标矢量，记做m表示第m个目标散射点，m=1,2,...M，M为观测场景目标空间Ω中的目标散射点总数，构成观测场景目标矩阵P_n，矩阵P_n的单个元素为P_n(l₀,m)，(l₀=1,2,3，l₀代表空间三维坐标；m=1，2....M,代表第m个目标散射点)。观测场景目标空间Ω在阵列SAR成像方案设计中已经确定。利用传统的合成孔径雷达原始回波仿真方法产生阵列SAR的原始回波数据，经过相参积累得到的雷达回波数据帧数为K，K帧回波数据中第ü帧回波数据记做Sⁱⁱ(ü=1,2...K)，初始化T行M列测量矩阵ψ，用于存放经残余相位补偿后的回波数据。

步骤3：初始化天线相位中心轨迹

自旋式飞行器雷达的天线收发模式如图1所示，初始化天线相位中心轨迹的具体步骤如下：

步骤3.1初始化步骤1中的雷达平台天线相位中心位置矩阵APC，矩阵APC的第一行代表自旋式飞行器沿航迹向坐标，第二行代表自旋式飞行器切航迹向坐标，第三行代表自旋式飞行器高度向坐标。矩阵APC的每个元素为APC(l,t)，（l=1,2,3,代表空间坐标的3个维度；t=1,2....T,代表第t个慢时刻，T为步骤1中定义的为慢时刻总数），APC(1:3,t)代表第t个慢时刻所对应的天线相位中心坐标。初始化慢时刻t，令t=1，代表第1个慢时刻。

步骤3.2令 $\mathrm{APC}(1,t)=\frac{v\&times;t}{\mathrm{prf}}+0.5\&times;a\&times;{\left(\frac{t}{\mathrm{prf}}\right)}^2,\mathrm{APC}(2,t)=r\&times;\cos \left(\frac{\mathrm{\&omega;}\&times;t}{\mathrm{prf}}\right),\mathrm{APC}(3,t)=r\&times;\sin \left(\frac{\mathrm{\&omega;}\&times;t}{\mathrm{prf}}\right)+H,$ 若APC(3,t)＞H，令 $(2,t)=r\&times;\cos (\frac{\mathrm{\&omega;}\&times;t}{\mathrm{prf}}+\mathrm{\&pi;})(3,t)=r\&times;\sin (\frac{\mathrm{\&omega;}\&times;t}{\mathrm{prf}}+\mathrm{\&pi;})$ 其中，v、ω、a、r、H、prf分别为步骤1中定义的飞行器飞行速度、飞行器自旋角速度、飞行器飞行加速度、飞行器自旋半径、飞行器水平高度、雷达系统的脉冲重复频率。

步骤3.3令t=t+1，若t<T+1，转步骤3.2，否则，存储雷达平台天线相位中心位置矩阵APC，转步骤4。

步骤4：观测场景目标空间序列化

观测场景目标空间序列化处理步骤如下：

步骤4.1初始化步骤2中的观测场景目标矩阵P_n，P_n第一行代表目标沿航迹向坐标，第二行代表目标切航迹向坐标，第三行为由沿航迹与切航迹向坐标所决定的高度坐标。P_n的单个元素为P_n(l₀,m)，(l₀=1,2,3，l₀代表空间三个维度的坐标；m=1,2....M,代表第m个目标散射点，M代表观测场景目标空间Ω中的目标散射点总数)，P_n(1:3,m)代表第m个目标空间Ω中的目标散射点。

步骤4.2令x=1，y=1，m=1,观测场景目标空间Ω中，x代表场景沿航迹向散射点序号，y代表场景切航迹向散射点序号，m表示第m个目标散射点，z(x,y)表示步骤2中由沿航迹向散射点序号与切航迹向散射点序号共同确定的目标散射点的高度。转步骤4.3。

步骤4.3观测场景目标矩阵P_n，矩阵第一行 $P_n(1,m)=(x-\frac{5x}{2})\&times;{\mathrm{\&rho;}}_a,$ 第二行 $P_n(2,m)=(y-\frac{5y}{2})\&times;{\mathrm{\&rho;}}_b,$ 第三行P_n (3,m)=round(z(x,y))×ρ_r，其中，s_x为场景沿航迹散射点总数，s_y为场景切航迹散射点总数，ρ_a为沿航迹向分辨率，ρ_b为切航迹向分辨率，ρ_r为高度向分辨率，转步骤4.4。

步骤4.4令x=x+1，m=m+1，若x<s_x+1，转步骤4.3，否则，转步骤4.5。

步骤4.5令x=1，y=y+1，若m<M+1，转步骤4.3，否则，转步骤5。

步骤5：用后像投影算法进行相参积累

后像投影方法步骤如下：

步骤5.1取出步骤2中的K帧原始回波数据；

首先，对K帧回波数据s¹,s²...s^K的高度向采用传统的匹配滤波压缩算法进行脉冲压缩，得到脉冲压缩后的K帧回波数据矩阵x¹,x²...x^K，匹配滤波的参考函数为：

H(t₀)=exp(-jπf_drt₀ ²)

其中，exp(·)表示e指数运算符号，f_dr为步骤1种初始化的线性调频信号的调频斜率，t₀为时间变量，T_r为步骤1中初始化的雷达平台天线发射脉冲时宽；

然后，对脉冲压缩后的K帧回波数据矩阵X¹,X²...X^k进行辛格插值，得到扩展为对应到T个慢时刻所对应的M列（M代表观测场景目标空间Ω中的目标散射点总数）场景目标向量其中，为第1个慢时刻点所对应的目标散射点回波数据，为第2个慢时刻点所对应的目标散射点回波数据，为第T个慢时刻点所对应的目标散射点回波数据，(t=1,2...T，t为第t个慢时刻，T为步骤1中定义的慢时刻总数)的每一个元素对应着场景目标散射点的回波；

最后，将T个慢时刻所对应的场景目标向量存储到计算机中，转步骤5.2。

步骤5.2令t=1，t代表第t个慢时刻，1≤t≤T，T为步骤1中定义的慢时刻总数，初始化迭代变量m=0，m代表第m个目标散射点，1≤m≤M，M为步骤1中定义的目标散射点总数，转步骤5.3。

步骤5.3令m=m+1，计算第m个目标散射点P_n(1:3,m)到的第t个慢时刻的天线相位中心APC(1:3,t)的距离，记为R_m：

$\begin{array}{c}{R}_m={\left|\right|\mathrm{APC}(1:3,t)-P_n(1:3,m)\left|\right|}_2\\ =\sqrt{{[\mathrm{APC}(1,t)-P_n(1,m)]}^2+{[\mathrm{APC}(2,t)-P_n(2,m)]}^2+{[\mathrm{APC}(3,t)-P_n(3,m)]}^2}\end{array}$

所以，第t个慢时刻到第m个目标散射点的回波延时：

$\mathrm{\&tau;}(t,m)=\frac{{2R}_m}{C}$

APC为步骤3中得到的雷达平台天线相位中心位置矩阵，雷达平台天线相位中心位置矩阵APC的每列（APC(1:3,t)）代表一个慢时刻点，P_n为步骤4中的序列化后的观测场景目标矩阵，矩阵P_n的每1列代表1个目标散射点的坐标，转步骤5.4。

步骤5.4计算得到观测场景目标空间Ω中第m个目标散射点与第t个慢时刻对应的时延函数，记为公式为：

(m=1,2...M)

其中，exp(·)表示e指数运算符号，f_c为步骤1中初始化得到的雷达工作中心频率，τ(t,m)为步骤5.3得到的第t个慢时刻到第m个目标散射点的回波延时，f_dr为步骤1中初始化得到的发射信号调频斜率，j为虚数单位(即-1的开平方根)，若迭代变量m<M+1，转步骤5.3，否则，转步骤5.5。

步骤5.5取出步骤5.1得到的第t个慢时刻所对应的场景目标向量其中，的第m个元素为a_t(m)，所以，将a_t(m)与步骤5.4得到的进行共轭相乘后得到经残余相位补偿后的回波数据存入步骤2中的测量矩阵ψ的第t行m列。当m=M时，令t=t+1，m=0。若t=T+1，转步骤5.6，否则，转步骤5.3。

步骤5.6输出测量矩阵ψ，ψ为经残余相位补偿后的测量矩阵，矩阵的每个元素对应一个散射点的回波数据φ_t(m)，记测量矩阵ψ的第t行所代表的第t个慢时刻所对应的场景目标向量为β_t，ψ化简为：

其中，φ₁(1)为观测场景目标空间Ω中第1个散射点在第1个慢时刻点所测得的回波函数，φ₁(2)为观测场景目标空间Ω中第2个散射点在第1个慢时刻点所测得的回波函数，φ₁(M)为观测场景目标空间Ω中第M个散射点在第1个慢时刻点所测得的回波函数；φ₂(1)为观测场景目标空间Ω中第1个散射点在第2个慢时刻点所测得的回波函数，φ₂(2)为观测场景目标空间Ω中第2个散射点在第2个慢时刻点所测得的回波函数；φ₂(M)为观测场景目标空间Ω中第M个散射点在第2个慢时刻点所测得的回波函数；φ_T(1)为观测场景目标空间Ω中第1个散射点在第T个慢时刻点所测得的回波函数，φ_T(2)为观测场景目标空间Ω中第2个散射点在第T个慢时刻点所测得的回波函数，φ_T(M)为观测场景目标空间Ω中第M个散射点在第T个慢时刻点所测得的回波函数；β₁为与第1个慢时刻点对应的场景目标回波向量，β₂为与第2个慢时刻点对应的场景目标回波向量，β_T为与第T个慢时刻点对应的场景目标回波向量，三维SAR的测量矩阵ψ为T行M列的二维矩阵。

步骤5.6取出步骤5.5中的场景目标向量β₁，β₂......β_T，将向量β₁，β₂......β_T相加，得到T个慢时刻相参积累后的三维空间成像场景向量：（t为第t个慢时刻），σ中每个元素的绝对值对应一个目标散射点的高度。

步骤6、输出飞行器测量高度

取出步骤5中对T个慢时刻相参积累后的三维空间成像场景向量σ，由σ的生成过程可知，的绝对值为飞行器天底点高度，即为测得的飞行器对地相对高度高度，其中round为定义的取整运算符，H为步骤1中初始化的飞行器飞行高度，M为观测场景目标空间Ω中散射点总数。

本发明的创新点在于针对自旋式飞行器提出了基于对目标场景三维成像的测高新模型，通过对目标场景三维成像，提取所需目标散射点的高度数据，完成对飞行器的高度测量。它是通过计算成像区域内每一采样点到合成孔径长度内雷达天线相位中心之间的双程时延，然后将对应的时域回波信号进行相干累加，从而恢复出每个采样点的散射系数信息。后向投影算法不基于参考点成像，能够实现每个采样点的精确聚焦，同时后向投影算法的雷达天线相位中心物理意义清晰，能够精确计算每个方位时刻的回波延时相位，补偿平台抖动引入的相位误差。同时，在天线相位中心形成过程中，由于自旋式飞行器在变加速直线飞行时伴随飞行器自旋运动，阵列天线是非等间隔分布的，将三维BP算法应用于自旋式飞行器平台三维成像处理，克服了传统的三维RD等频域算法要求阵列天线等间隔采样的缺点，具有聚焦效果好、适用于非均匀或欠采样条件、对于变加速、非直线、自旋的运动轨迹能成像优点。

本发明的优点在于针对自旋式飞行器成像平台，由于自旋运动导致一般的脉冲积累方法无法实现回波能量聚集，传统的测高方法不再适用的问题，将自旋式飞行器的天线相位中心运动轨迹看做非均匀分布的等效天线阵列，利用等效天线阵列对场景三维成像，成像场景中心点高度即为雷达天底点高度，结合飞行器的水平高度数据，得出飞行器飞行的相对高度。同时，针对等效阵列非均匀分布的特点，本发明采用了三维BP（后像投影）成像算法。与传统方法相比，本发明克服了传统的三维RD等频域算法要求阵列天线等间隔采样的缺点，具有聚焦效果好、适用于非均匀或欠采样条件、对于变加速、非直线、自旋的运动轨迹能成像优点。本发明可以应用于合成孔径雷达成像和引导测高雷达等领域。

附图说明

图1为自旋式飞行器雷达收发天线结构

其中，自旋式飞行器SAR平台采用对称分布双天线架构。飞行器在沿箭头所指方向直线飞行的同时伴随着自旋，指向朝下的天线(天线A)处于工作状态，发射并接收回波，指向朝上的天线(天线B)不工作。当飞行器转动到双天线同时位于水平面时，原来工作的天线(天线A)停止发射或接受回波，原来不工作的天线（天线B）开始工作。即飞行器每旋转180度，天线在工作/停止状态间交替一次。

图2为三维成像回波延时模型

其中，任取场景中的3个散射点（点目标散射点1、点目标散射点2、点目标散射点3），在任一慢时刻，雷达到这3个散射点的距离分别记为（R₁、R₂、R₃），则回波延时分别为 $\mathrm{\&tau;}\left(1\right)=\frac{{2R}_1}{C},\mathrm{\&tau;}\left(2\right)=\frac{{2R}_2}{C},\mathrm{\&tau;}\left(3\right)=\frac{{2R}_3}{C},$ C为光速。

图3为天线相位中心存储矩阵APC

其中，矩阵APC每一列代表一个慢时刻点的自旋式飞行器轨迹坐标，如果有T个慢时刻点，那么矩阵APC就有T列。目标轨迹动态存储矩阵APC每一列有3行：矩阵APC第一行代表自旋式飞行器沿航迹向坐标，第二行为自旋式飞行器切航迹向坐标，第三行为自旋式飞行器高度向坐标。

图4为三维BP算法流程图

具体步骤为对先原始回波数据进行匹配滤波得到经过处理的回波数据，再计算每个天线相位中心到所有目标散射点的回波延时得到相位补偿因子，然后对经过匹配滤波处理的回波数据进行插值重采样得到目标散射点对应的回波数据，最后将目标散射点对应的回波数据与相位补偿因子相乘得到相位补偿后的回波数据，相干叠加后取绝对值输出三维场景成像结果。

具体实施方式

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

步骤1、初始化自旋式飞行器SAR成像系统各参数

该方法需要初始化的系统参数包括：光在空气中的传播速度C=3000Km/s，飞行器水平高度H=3000m，飞行器自旋半径r=2m，飞行器飞行速度v=2000m，最大加速度a=200m/s，自旋角频率ω=6πrad/s，雷达系统工作的信号波长λ=0.0554m，雷达发射线性调频信号的载频f_c=5.415×10⁹Hz，线性调频信号的调频斜率f_dr=3×10¹³Hz，雷达平台天线发射信号带宽B=3×10⁷Hz，雷达平台天线发射脉冲时宽T_r=10^-6s，雷达平台接收系统采样频率F_s=33×10⁶Hz，雷达发射脉冲重复频率prf=1000Hz，第t个慢时刻t=1,2...T，慢时刻总数T=2000，雷达平台天线相位中心位置矩阵APC的初始位置矢量APC(1:3,1)=[0,0,H]，沿航迹向分辨率ρ_a=50m，切航迹向分辨率ρ_b=100m，高度向分辨率ρ_r=5m。

步骤2、初始化阵列三维SAR的观测场景目标空间

该方法需要初始化的阵列三维SAR观测场景目标空间参数包括：初始化观测场景目标空间Ω的大小为256×256×256像素，观测场景目标空间Ω的中心坐标位置位于[0,0,0]，每一个单元网格在水平横向、水平纵向和高度向边长为d_x=50m、d_y=100m和d_z=5m，在观测场景目标空间Ω的每个单元格中加入仿真点目标散射体，单元网格的每一对水平横向、水平纵向所对应的64个高度单元格中，只有一个有点目标散射体，故观测场景目标空间Ω中的目标散射点总数为M=65536，它们散射系数值均为1，观测场景目标空间Ω中没有包含点目标单元格的散射系数设置为0。初始化DEM数据，观测场景目标空间Ω中，第m个目标散射点的位置矢量为 ${\stackrel{\&RightArrow;}{P}}_m=[(x\&prime;-128)\&CenterDot;{\mathrm{\&rho;}}_a,(y\&prime;-128)\&CenterDot;{\mathrm{\&rho;}}_b,z(x\&prime;,y\&prime;)\&CenterDot;{\mathrm{\&rho;}}_r],$ 其中x′为场景沿航迹向散射点序号，x′=1,2...256，y′为场景切航迹向散射点序号，y′=1,2...256，z(x′,y′)为由x′、y′所对应的散射点高度向坐标，z(x′,y′)取值为1,2...64，第m个散射点与x′、y′的关系为m=(x′-1)·256+y′，m的取值为1,2,...M，M为观测场景目标空间Ω中的目标散射点总数，ρ_a、ρ_b、ρ_r分别为步骤1中初始化的沿航迹向分辨率、切航迹向分辨率、高度向分辨率，观测场景目标空间Ω中，向量构成观测场景目标矩阵P_n，矩阵P_n的单个元素为P_n(l₀,m)，(l₀=1,2,3，l₀代表空间三维坐标；m=1,2....M,代表第m个目标散射点)。利用传统的合成孔径雷达原始回波仿真方法产生阵列SAR的原始回波数据，相参处理的雷达回波数据帧数K=5，K帧原始回波数据中第ü帧原始回波数据记做Sⁱⁱ(ü=1,2...K),初始化T行M列测量矩阵ψ，用于存放经残余相位补偿后的回波数据。

步骤3：初始化天线相位中心轨迹

初始化天线相位中心轨迹分个3小步骤：

步骤3.1初始化步骤1中的雷达平台天线相位中心位置矩阵APC。矩阵APC的第一行代表自旋式飞行器沿航迹向坐标，第二行代表自旋式飞行器切航迹向坐标，第三行代表自旋式飞行器高度向坐标。APC的单个元素为APC(l,t)，（l=1,2,3,代表空间坐标的3个维度；t=1,2...T,代表第t个慢时刻），T为步骤1中定义的为慢时刻总数，APC(1:3,t)代表第t个慢时刻所对应的天线相位中心坐标。初始化迭代变量t=1，代表第1个慢时刻。

步骤3.2令 $\mathrm{APC}(1,t)=\frac{v\&times;t}{\mathrm{prf}}+0.5\&times;a\&times;{\left(\frac{t}{\mathrm{prf}}\right)}^2,\mathrm{APC}(2,t)=r\&times;\cos \left(\frac{\mathrm{\&omega;}\&times;t}{\mathrm{prf}}\right),\mathrm{APC}(3,t)=r\&times;\sin \left(\frac{\mathrm{\&omega;}\&times;t}{\mathrm{prf}}\right)+H,$ 若APC（3，t＞H，令 $(2,t)=r\&times;\cos (\frac{\mathrm{\&omega;}\&times;t}{\mathrm{prf}}+\mathrm{\&pi;})(3,t)=r\&times;\sin (\frac{\mathrm{\&omega;}\&times;t}{\mathrm{prf}}+\mathrm{\&pi;})$ 其中，v、ω、a、r、H、prf分别为步骤1中初始化的飞行器速度、飞行器自旋角速度、飞行器加速度、飞行器自旋半径、飞行器水平高度、雷达系统的脉冲重复频率，转步骤3.3。

步骤3.3令t=t+1，若t<T+1，转步骤3.2，否则，存储矩阵APC，转步骤4。

步骤4：观测场景目标空间序列化

观测场景目标空间序列化处理分5个小步骤：

步骤4.1取出步骤2中的观测场景目标矩阵P_n，P_n第一行代表目标沿航迹向坐标，第二行代表目标切航迹向坐标，第三行为由沿航迹与切航迹向坐标所决定的高度坐标。P_n的单个元素为P_n(l₀,m)，(l₀=1,2,3，l₀代表空间三维坐标；m=1,2....M,代表第m个目标散射点，M为观测场景目标空间Ω中的目标散射点总数)。P_n(1:3,m)代表第m个目标空间Ω中的目标散射点。

步骤4.2初始化迭代变量x=1，y=1，m=1,x代表场景沿航迹向散射点序号，y代表场景切航迹向散射点序号，m代表第m个目标散射点，z(x,y)代表步骤2中由沿航迹向散射点序号和切航迹向散射点序号共同确定的目标散射点的高度，转步骤4.3。

步骤4.3目标矩阵P_m第一行 $P_n(1,m)=(x-\frac{5x}{2})\&times;{\mathrm{\&rho;}}_a,$ 第二行 $P_n(2,m)=(y-\frac{5y}{2})\&times;{\mathrm{\&rho;}}_b,$ 第三行P_m(3,m)=round(z(x,y))×ρ_r。其中，s_x为场景沿航迹散射点总数，s_y为场景切航迹散射点总数，ρ_a为沿航迹向分辨率，ρ_b为切航迹向分辨率，ρ_r为高度向分辨率，转步骤4.4。

步骤4.4令x=x+1，m=m+1，若x<s_x+1，转步骤4.3，否则，转步骤4.5。

步骤4.5令x=1，y=y+1，若m小等于场景总的散射点数M，转步骤4.3，否则，转步骤5。

步骤5：用后像投影算法进行相参积累

步骤5.1取出步骤2中的K帧原始回波数据。对K帧回波数据s¹,s²...s^K的高度向采用传统的匹配滤波压缩算法进行脉冲压缩，得到脉冲压缩后的K帧回波数据矩阵x¹,x²...x^K，匹配滤波的参考函数为：

H(t₀)=exp（-jπf_dr t₀ ²）

其中，exp(·)表示e指数运算符号，f_dr为步骤1种初始化的线性调频信号的调频斜率，t₀为时间，T_r为步骤1中初始化的雷达平台天线发射脉冲时宽。然后将脉冲压缩后的K帧回波数据矩阵X¹,X²...X^k进行辛格插值，扩展为对应到T个慢时刻的M维场景目标向量其中，为第1个慢时刻点所对应的场景目标回波数据，为第2个慢时刻点所对应的场景目标回波数据，为第T个慢时刻点所对应的场景目标回波数据，(t=1,2...T，t为第t个慢时刻)的每一个元素对应着场景目标一个散射点的回波。

最后将这T个慢时刻所对应的场景目标向量存储到计算机中，转步骤5.2。

步骤5.2。初始化迭代变量t=1，t代表第t个慢时刻，1≤t≤T，T为步骤1中定义的慢时刻总数，初始化迭代变量m=0，m代表第m个目标散射点，1≤m≤M，M为步骤1中定义的观测场景目标空间Ω中的目标散射点总数，转步骤5.3。

步骤5.3令m=m+1，计算第m个目标散射点P_n(1:3,m)到第t个慢时刻的天线相位中心APC(1:3,t)的距离，记为R_m：

$\begin{array}{c}{R}_m={\left|\right|\mathrm{APC}(1:3,t)-P_n(1:3,m)\left|\right|}_2\\ =\sqrt{{[\mathrm{APC}(1,t)-P_n(1,m)]}^2+{[\mathrm{APC}(2,t)-P_n(2,m)]}^2+{[\mathrm{APC}(3,t)-P_n(3,m)]}^2}\end{array}$

所以，第t个慢时刻到第m个目标散射点的回波延时：

$\mathrm{\&tau;}(t,m)=\frac{{2R}_m}{C}$

APC为步骤3中生成的天线相位中心矩阵，矩阵APC的每列代表一个慢时刻点，P_n为步骤4中的序列化后的观测场景目标矩阵，矩阵P_n的每1列代表1个场景点目标的坐标，转步骤5.4。

步骤5.4计算得到观测场景目标空间Ω中第m个目标散射点与第t个慢时刻对应的时延函数，记为公式为：

(m=1,2...M)

其中，exp(·)表示e指数运算符号，f_c为步骤1中初始化得到的雷达工作中心频率，τ(t,m)为第t个慢时刻到第m个目标散射点的回波延时，f_dr为步骤1中初始化得到的发射信号调频斜率，j为虚数单位(即-1的开平方根)，若迭代变量m<M+1，转步骤5.3，否则，转步骤5.5。

步骤5.5取出步骤5.1得到的第t个慢时刻所对应的场景目标向量其中，的第m个元素为a_t(m)，所以，将步骤5.4的与a_t(m)进行共轭相乘后得到经残余相位补偿后的回波数据存入步骤2中的测量矩阵ψ的第t行m列。当m=M时，令t=t+1，m=0。若t=T+1，转步骤5.6，否则，转步骤5.3。

步骤5.6输出测量矩阵ψ，矩阵ψ的每个元素对应一个慢时刻的一个散射点，测量矩阵ψ的第t行表示第t个慢时刻所对应的场景目标向量为β_t，ψ化简为：

其中，φ₁(1)为观测场景目标空间Ω中第1个散射点在第1个慢时刻点所测得的回波函数，φ₁(2)为观测场景目标空间Ω中第2个散射点在第1个慢时刻点所测得的回波函数，φ₁(M)为观测场景目标空间Ω中第M个散射点在第1个慢时刻点所测得的回波函数，φ₂(1)为观测场景目标空间Ω中第1个散射点在第2个慢时刻点所测得的回波函数，φ₂(2)为观测场景目标空间Ω中第2个散射点在第2个慢时刻点所测得的回波函数，φ₂(M)为观测场景目标空间Ω中第M个散射点在第2个慢时刻点所测得的回波函数，φ_T(1)为观测场景目标空间Ω中第1个散射点在第T个慢时刻点所测得的回波函数，φ_T(2)为观测场景目标空间Ω中第2个散射点在第T个慢时刻点所测得的回波函数，φ_T(M)为观测场景目标空间Ω中第M个散射点在第T个慢时刻点所测得的回波函数。β₁为与第1个慢时刻点对应的场景目标回波向量，β₂为与第2个慢时刻点对应的场景目标回波向量，β_T为与第T个慢时刻点对应的场景目标回波向量，三维SAR的测量矩阵ψ为T行M列的二维矩阵。

步骤5.6取出步骤5.5中的场景目标向量β₁，β₂......β_T，将向量β₁，β₂......β_T相加，得到T个慢时刻相参积累后的三维空间成像场景向量：三维空间成像场景向量σ每个元素的绝对值对应一个目标散射点的高度。

步骤6、输出飞行器测量高度

取出步骤5中T个慢时刻相参积累后的三维空间成像场景向量σ，即为飞行器天底点高度，，即为测得的飞行器对地相对高度高度，其中round为定义的取整运算符，H为步骤1中初始化的飞行器飞行高度，M为观测场景目标空间Ω中散射点总数。

通过本发明具体实施方式可以看出，本发明首先将自旋式飞行器天线相位中心轨迹等效为天线面阵，通过等效面阵天线对目标场景三维成像，提出了基于三维成像的自旋式飞行器天底点测高方法，将飞行器测高问题转化为对三维成像场景中心点高度的求解，通过对目标场景高精度三维成像，完成对自旋式飞行器的测高。在成像过程中，本发明提供了三维BP算法，该算法能适用于等效天线面阵非均匀分布的情况。利用自旋式飞行器运动过程中形成的等效天线阵列，结合自旋式飞行器三维SAR的系统参数和观测场景目标参数，通过在时域对有散射点的场景单元格补偿时延相位，完成了对目标场景的三维成像。本发明具有对于变加速、非直线运动的自旋式飞行器能成像优点。

Claims (1)

1.一种基于SAR三维成像的自旋式飞行器测高方法，其特征是它包括如下步骤： 

步骤1、初始化自旋式飞行器SAR成像系统各参数 

初始化自旋式飞行器SAR成像系统参数包括：飞行器水平高度，记做H，飞行器半径，记做r，飞行器飞行速度，记做v，最大加速度，记做a，自旋角频率，记做ω，由v、a、ω决定的飞行器速度矢量，记做雷达系统工作的信号波长，记做λ，雷达发射线性调频信号的载频，记做fc，线性调频信号的调频斜率，记做fdr，雷达平台天线发射信号带宽，记做B，雷达平台天线发射脉冲时宽，记做Tr，雷达平台接收系统采样频率，记做Fs，雷达发射脉冲重复频率，记做prf，光在空气中的传播速度，记做C，第t个慢时刻，记做t，t=1，2...T,T为慢时刻总数，雷达平台天线相位中心位置矩阵，记做APC，沿航迹向分辨率，记做ρa，切航迹向分辨率，记做ρ_b，高度向分辨率，记做ρ_r，接收波门相对发射信号延时T_d，距离门位置记做Ic；上述参数均为阵列SAR系统标准参数，其中雷达发射线性调频信号的载频fc，雷达载频波长λ，雷达发射基带信号的信号带宽B，雷达发射信号脉冲宽度Tp，雷达发射信号调频斜率f_dr，雷达接收系统的采样频率Fs，雷达系统的脉冲重复频率prf，光在空气中的传播速度C，雷达系统的脉冲重复时间PRI和雷达接收系统接收波门相对于发射信号发散波门的延迟TD，均在三维SAR系统设计过程中已经确定；飞行器飞行高度H，飞行器半径r，飞行器速度v，加速度a，自旋角频率ω，平台速度矢量沿航迹向分辨率ρ_a，切航迹向分辨率ρ_b，高度向分辨率ρ_r，雷达平台天线相位中心位置矩阵APC，均在自旋式飞行器三维SAR观测方案设计中已经确定；根据三维SAR雷达系统方案和三维SAR雷达成像观测方案，以上基于三维成像的自旋式飞行器测高系统参数均为已知； 

步骤2：初始化阵列三维SAR的观测场景目标空间 

初始化阵列三维SAR的观测场景目标空间参数，包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标作为阵列三维SAR的观测场景目标空间Ω；将三维SAR观测场景目标空间Ω均匀划分成大小相等的立体单元网格(亦称为分辨单元)，立体单元网格在水平横向、水平纵向和高度向边长分别记为dx、dy和dz，立体单元网格的每一对水平横向、水平纵向所对应高度单元格中，只有一个有点目标散射体，立体单元网格大小选择为阵列SAR系统传统理论成像分辨率的二分之一，初始化数字高程模型（DEM）数据，其中，场景沿航迹总的散射点数，记做sx，场景切航迹向总的散射点数，记做sy，场景沿航迹向上的第x0个散射点序号，记做x0，x₀=1,2,...s_x，场景沿航迹向上的第y0个散射点序号，记做y₀，y₀=1,2,...s_y，沿航迹向散射点位置x₀与切航迹向散射点位置y₀所唯一确定的目标散射点的高度，记做z(x₀,y₀)；观测场景目标空间Ω中第m个散射点的坐标矢量，记做m表示第m个目标散射点，m=1,2,...M，M为观测场景目标空间Ω中的目标散射点总数，构成观测场景目标矩阵P_n，矩阵P_n的单个元素为P_n(l₀,m)，l₀=1,2,3，l₀代表空间三维坐标；m=1,2,...M,代表第m个目标散射点；观测场景目标空间Ω在阵列SAR成像方案设计中已经确定；利用传统的合成孔径雷达原始回波仿真方 法产生阵列SAR的原始回波数据，经过相参积累得到的雷达回波数据帧数为K，K帧回波数据中第ü帧回波数据记做sⁱⁱ(ü=1,2...K)，初始化T行M列测量矩阵ψ，用于存放经残余相位补偿后的回波数据； 

步骤3：初始化天线相位中心轨迹 

自旋式飞行器雷达的天线收发模式，初始化天线相位中心轨迹的具体步骤如下： 

步骤3.1初始化步骤1中的雷达平台天线相位中心位置矩阵APC，矩阵APC的第一行代表自旋式飞行器沿航迹向坐标，第二行代表自旋式飞行器切航迹向坐标，第三行代表自旋式飞行器高度向坐标；矩阵APC的每个元素为APC(l,t)，（l=1,2,3,代表空间坐标的3个维度；t=1,2....T,代表第t个慢时刻，T为步骤1中定义的为慢时刻总数），APC(1:3,t)代表第t个慢时刻所对应的天线相位中心坐标；初始化慢时刻t，令t=1，代表第1个慢时刻； 

步骤3.2令若APC(3,t)＞H,令其中，v、ω、a、r、H、prf分别为步骤1中定义的飞行器飞行速度、飞行器自旋角速度、飞行器飞行加速度、飞行器自旋半径、飞行器水平高度、雷达系统的脉冲重复频率； 

步骤3.3令t=t+1，若t<T+1，转步骤3.2，否则，存储雷达平台天线相位中心位置矩阵APC，转步骤4； 

步骤4：观测场景目标空间序列化 

观测场景目标空间序列化处理步骤如下： 

步骤4.1初始化步骤2中的观测场景目标矩阵P_n，P_n第一行代表目标沿航迹向坐标，第二行代表目标切航迹向坐标，第三行为由沿航迹与切航迹向坐标所决定的高度坐标；P_n的单个元素为P_n(l₀,m)，l₀=1,2,3，l₀代表空间三个维度的坐标；m=1,2....M,代表第m个目标散射点，M代表观测场景目标空间Ω中的目标散射点总数，P_n(1:3,m)代表第m个目标空间Ω中的目标散射点； 

步骤4.2令x=1，y=1，m=1,观测场景目标空间Ω中，x代表场景沿航迹向散射点序号，y代表场景切航迹向散射点序号，m表示第m个目标散射点，z(x,y)表示步骤2中由沿航迹向散射点序号与切航迹向散射点序号共同确定的目标散射点的高度；转步骤4.3； 

步骤4.3观测场景目标矩阵P_n，矩阵第一行第二行第三行P_n (3,m)=round(z(x,y))×ρ_r，其中，s_x为场景沿航迹散射点总数，s_y为场景切航迹散射点总数，ρ_a为沿航迹向分辨率，ρ_b为切航迹向分辨率，ρ_r为高度向分辨率，转步骤4.4； 

步骤4.4令x=x+1，m=m+1，若x<s_x+1，转步骤4.3，否则，转步骤4.5； 

步骤4.5令x=1，y=y+1，若m<M+1，转步骤4.3，否则，转步骤5； 

步骤5：用后像投影算法进行相参积累 

后像投影方法步骤如下： 

步骤5.1取出步骤2中的K帧原始回波数据； 

首先，对K帧回波数据s¹,s²...s^K的高度向采用传统的匹配滤波压缩算法进行脉冲压缩，得到脉冲压缩后的K帧回波数据矩阵x¹,x²...x^K，匹配滤波的参考函数为： 

H(t₀)=exp(-jπf_drt₀ ²) 

其中，exp(·)表示e指数运算符号，f_dr为步骤1种初始化的线性调频信号的调频斜率，t₀为时间变量， T_r为步骤1中初始化的雷达平台天线发射脉冲时宽； 

然后，对脉冲压缩后的K帧回波数据矩阵x¹,x²...x^K进行辛格插值，得到扩展为对应到T个慢时刻所对应的M列（M代表观测场景目标空间Ω中的目标散射点总数）场景目标向量其中，为第1个慢时刻点所对应的目标散射点回波数据，为第2个慢时刻点所对应的目标散射点回波数据，为第T个慢时刻点所对应的目标散射点回波数据，(t=1,2...T，t为第t个慢时刻，T为步骤1中定义的慢时刻总数)的每一个元素对应着场景目标散射点的回波； 

最后，将T个慢时刻所对应的场景目标向量存储到计算机中，转步骤5.2； 

步骤5.2令t=1，t代表第t个慢时刻，1≤t≤T，T为步骤1中定义的慢时刻总数，初始化迭代变量m=0，m代表第m个目标散射点，1≤m≤M，M为步骤1中定义的目标散射点总数，转步骤5.3； 

步骤5.3令m=m+1，计算第m个目标散射点P_n(1:3,m)到的第t个慢时刻的天线相位中心APC(1:3,t)的距离，记为R_m： 

所以，第t个慢时刻到第m个目标散射点的回波延时： 

APC为步骤3中得到的雷达平台天线相位中心位置矩阵，雷达平台天线相位中心位置矩阵APC的每列（APC(1:3,t)）代表一个慢时刻点，P_n为步骤4中的序列化后的观测场景目标矩阵，矩阵P_n的每1列代表1个目标散射点的坐标，转步骤5.4； 

步骤5.4计算得到观测场景目标空间Ω中第m个目标散射点与第t个慢时刻对应的时延函数，记为公式为： 

(m=1,2...M) 

其中，exp(·)表示e指数运算符号，f_c为步骤1中初始化得到的雷达工作中心频率，τ(t,m)为步骤5.3得到的第t个慢时刻到第m个目标散射点的回波延时，f_dr为步骤1中初始化得到的发射信号调频斜率，j为虚数单位(即-1的开平方根)，若迭代变量m<M+1，转步骤5.3，否则，转步骤5.5； 

步骤5.5取出步骤5.1得到的第t个慢时刻所对应的场景目标向量其中，的第m个元素为a_t(m)，所以，将a_t(m)与步骤5.4得到的进行共轭相乘后得到经残余相位补偿后的回波数据 存入步骤2中的测量矩阵ψ的第t行m列；当m=M时，令t=t+1，m=0；若t=T+1，转步骤5.6，否则，转步骤5.3； 

步骤5.6输出测量矩阵ψ，ψ为经残余相位补偿后的测量矩阵，矩阵的每个元素对应一个散射点的回波数据φ_t(m)，记测量矩阵ψ的第t行所代表的第t个慢时刻所对应的场景目标向量为β_t，ψ化简为： 

其中，φ₁(1)为观测场景目标空间Ω中第1个散射点在第1个慢时刻点所测得的回波函数，φ₁(2)为观测场景目标空间Ω中第2个散射点在第1个慢时刻点所测得的回波函数，φ₁(M)为观测场景目标空间Ω中第M个散射点在第1个慢时刻点所测得的回波函数；φ₂(1)为观测场景目标空间Ω中第1个散射点在第2个慢时刻点所测得的回波函数，φ₂(2)为观测场景目标空间Ω中第2个散射点在第2个慢时刻点所测得的回波函数；φ₂(M)为观测场景目标空间Ω中第M个散射点在第2个慢时刻点所测得的回波函数；φ_T(1)为观测场景目标空间Ω中第1个散射点在第T个慢时刻点所测得的回波函数，φ_T(2)为观测场景目标空间Ω中第2个散射点在第T个慢时刻点所测得的回波函数，φ_T(M)为观测场景目标空间Ω中第M个散射点在第T个慢时刻点所测得的回波函数；β₁为与第1个慢时刻点对应的场景目标回波向量，β₂为与第2个慢时刻点对应的场景 目标回波向量，β_T为与第T个慢时刻点对应的场景目标回波向量，三维SAR的测量矩阵ψ为T行M列的二维矩阵； 

步骤5.6取出步骤5.5中的场景目标向量β₁，β₂......β_T，将向量β₁，β₂......β_T相加，得到T个慢时刻相参积累后的三维空间成像场景向量：（t为第t个慢时刻），σ中每个元素的绝对值对应一个目标散射点的高度； 

步骤6、输出飞行器测量高度 

取出步骤5中对T个慢时刻相参积累后的三维空间成像场景向量σ，由σ的生成过程可知， 的绝对值为飞行器天底点高度，即为测得的飞行器对地相对高度高度，其中round为定义的取整运算符，H为步骤1中初始化的飞行器飞行高度，M为观测场景目标空间Ω中散射点总数。