CN104155653A - 一种基于特征距离子空间的sar后向投影成像方法 - Google Patents

Abstract

本发明公开了一种基于特征距离子空间的SAR后向投影成像方法，它是通过以距离历史向量为基础，从传统BP算法能否聚焦的角度提出了距离历史向量匹配度的概念；然后基于最小均方误差准则和主成分分析准则，在SAR成像领域，本发明提出了特征距离子空间的概念；因为在特征距离子空间中最优地修正了测量距离历史向量，所以很好地抑制了测量噪声，提高了距离历史向量的匹配度，从而对补偿的相位进行了准确地估计，最终获得了良好的聚焦效果。由于本发明适用于任意体制的SAR二维成像和三维成像，所以它具有很强的可推广性。最后，由于在特征距离子空间中修正测量距离历史向量是一种线性运算，所以本发明具有很强的可操作性。

Description

一种基于特征距离子空间的SAR后向投影成像方法

技术领域：

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

背景技术：

在众多的合成孔径雷达(synthetic aperture radar，SAR)成像算法中，后向投影(back projection，BP)算法以其适应能力强和聚焦性能卓越而著称(详见参考文献“P.Zhang,X.Zhang,and G.Fang.Comparison of the imaging resolutionsof time reversal and back-projection algorithms in EM inverse scattering.IEEEGeoscience and Remote Sensing Letters,2013,10(6):357-361”)。然而，在现实中，由于用于相位补偿的距离历史只能通过实际的测量来获得，所以这种距离历史就不可避免地要受到测量噪声的干扰。因为由测量噪声所引起的相位误差会直接影响BP算法的聚焦效果，于是只有对相位误差进行准确地估计，BP算法才能获得良好地聚焦。

近年来，发展起来的自聚焦算法可以根据回波数据对相位误差进行准确地估计，其在SAR成像处理中占据着非常重要的地位。现有的自聚焦算法主要分为两类。第一类自聚焦算法是基于图像特显点的方法，其代表算法为相位梯度自聚焦(详见参考文献“D.E.Wahl,P.H.Eichel,D.C.Ghiglia,and C.V.Jakowatz.Phasegradient autofocus—A robust tool for high resolution SAR phase correction.IEEETransactions on Aerospace and Electronic Systems,1994,30(3):827–835”)，该算法利用多个强点目标的相位历程来对相位误差进行估计。第二类自聚焦算法是基于图像质量的估计方法，该类算法首先将图像的整体信息作为评价准则，目前常用的评价准则包括：熵最小化(详见参考文献“M.Liu,C.Li,and X.Shi.Aback-projection fast autofocus algorithm based on minimum entropy for SAR imaging.3rd International Asia-Pacific Conference on Synthetic Aperture Radar(APSAR),2011:1-4”)、对比度最大化(详见参考文献“J.Kolman.PACE:An autofocusalgorithm for SAR.IEEE International Radar Conference,2005:310-314”)、图像锐度最大化(详见参考文献“N.A.Joshua.An autofocus method for back projectionimagery in synthetic aperture radar.IEEE Geoscience and Remote Sensing Letters,2012,9(1):104-108”)，然后该类算法通过调节相位的估计值使得图像的质量达到最优，以此来确定补偿的相位。然而，第一类自聚焦算法过度依赖场景中的强点目标，第二类自聚焦算法其本质属于图像域算法，它们都没有从直接抑制测量噪声的角度出发去估计补偿的相位。

如何通过直接抑制测量噪声来估计补偿的相位，从而达到BP算法自聚焦的目的，据发明人所知这个方向目前还没有公开出版的文献进行过报道。

发明内容：

本发明提出了一种基于特征距离子空间的SAR后向投影成像方法，它是基于传统的BP算法，首先以距离历史向量为基础，从BP算法能否聚焦的角度引入了距离历史向量匹配度的概念；然后根据传统的最小均方误差准则与主成分分析准则，在SAR成像领域，提出了特征距离子空间的概念，进而提出了一种新的BP算法，即基于特征距离子空间的BP成像方法。由于本发明较好地抑制了测量噪声，所以提高了距离历史向量的匹配度，从而对补偿的相位进行了准确地估计，最终获得了良好的聚焦效果。

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

定义1、回波距离历史

通常，SAR的发射信号为线性调频(linear frequency modulation，LFM)信号。设观测场景为Ω，点目标P∈Ω。如果不考虑天线方向图对回波信号的加权影响，那么点目标P的基带回波信号S(t,u；P)为：

$\begin{array}{c}S(t,u;P)=\mathrm{\σ}\left(P\right)\·\mathrm{rect}\left(\frac{t-R(u;P)/c}{T_P}\right)\·\mathrm{rect}\left(\frac{u}{T}\right)\\ \×\exp [\mathrm{j\π}\frac{B}{T_P}{(t-\frac{R(u;P)}{c})}^2-{j2\mathrm{\πf}}_c\frac{R(u;P)}{c}]\end{array}$

其中，t表示距离向快时间；u表示方位向慢时间；σ(P)表示点目标P的后向散射系数；B表示雷达发射信号的距离带宽；T_P表示雷达发射信号的脉冲宽度；T表示合成孔径时间；f_c表示雷达工作中心频率；c表示电磁波的传播速度；rect(u/T)表示慢时间u的窗函数。在上式中，R(u；P)表示点目标P的距离历史，由于R(u；P)包含在回波信号中，所以本发明将R(u；P)定义为回波距离历史。

定义2、回波距离历史向量

针对回波距离历史R(u；P)，通过在[-T/2,T/2]区间内对慢时间u进行均匀采样，可以得到向量R(P)＝[R(u₁；P),R(u₂；P),…,R(u_N；P)]^T，其中是u的采样点，N表示方位向采样点数，采样间隔为雷达的脉冲重复时间。由于向量R(P)是由回波距离历史R(u；P)所产生的，所以本发明把R(P)定义为回波距离历史向量。

定义3、传统BP算法模型

在这里，首先对基带回波信号S(t,u；P)进行距离脉冲压缩，经过距离脉冲压缩以后的信号S_c(t,u；P)为：

$\begin{array}{c}{S}_c(t,u;P)=\mathrm{\σ}\left(P\right)\·\mathrm{rect}\left(\frac{t-R(u;P)/c}{T_P}\right)\·\mathrm{rect}\left(\frac{u}{T}\right)\\ \×\sin c\left[\mathrm{\πB}(t-\frac{R(u;P)}{c})\right]\·\exp (-j2\mathrm{\π}\frac{R(u;P)}{\mathrm{\λ}})\end{array}$

其中，t表示距离向快时间；u表示方位向慢时间；σ(P)表示点目标P的后向散射系数；B表示雷达发射信号的距离带宽；T_P表示雷达发射信号的脉冲宽度；T表示合成孔径时间；c表示电磁波的传播速度；λ表示雷达载波波长；rect(u/T)表示慢时间u的窗函数；R(u；P)表示回波距离历史。然后，通过求解下式，便可获得点目标P的SAR图像I(P)，即：

$I\left(P\right)=\underset{-T/2}{\overset{T/2}{\∫}}{S}_c(\frac{R_m(u;P)}{c},u;P)\·\exp \left(j2\mathrm{\π}\frac{R_m(u;P)}{\mathrm{\λ}}\right)\mathrm{du}$

其中，R_m(u；P)为测量距离历史，对于R_m(u；P)的描述由下面的定义4给出。本发明把上式定义为传统BP算法模型。

定义4、测量距离历史

在传统BP算法模型中，R_m(u；P)表示通过实际测量所获得的点目标P的距离历史，本发明把R_m(u；P)定义为测量距离历史。在现实中，由于回波距离历史R(u；P)是未知的，而且是无法获取的，所以在传统BP算法中用于相位补偿的距离历史只能是已知的测量距离历史R_m(u；P)。

定义5、测量距离历史向量

针对测量距离历史R_m(u；P)，通过在[-T/2,T/2]区间内对慢时间u进行均匀采样，可以得到向量R_m(P)＝[R_m(u₁；P),R_m(u₂；P),…,R_m(u_N；P)]^T，其中是u的采样点，N表示方位向采样点数，采样间隔为雷达的脉冲重复时间。由于向量R_m(P)是由测量距离历史R_m(u；P)所产生的，所以本发明把R_m(P)定义为测量距离历史向量。

定义6、测量噪声

现实中，因为测量距离历史R_m(u；P)是通过实际测量而获得的，所以由于一些客观因素的影响，R_m(u；P)中不可避免地要包含噪声，本发明把该噪声定义为测量噪声，同时令该测量噪声为e_m(u；P)。需要指出的是，测量距离历史R_m(u；P)、回波距离历史R(u；P)、测量噪声e_m(u；P)三者之间满足如下关系：

R_m(u；P)＝R(u；P)+e_m(u；P)

在本发明中，e_m(u；P)为一个随机过程。

定义7、测量噪声向量

针对测量噪声e_m(u；P)，通过在[-T/2,T/2]区间内对慢时间u进行均匀采样，可以得到向量e_m(P)＝[e_m(u₁；P),e_m(u₂；P),…,e_m(u_N；P)]^T，其中是u的采样点，N表示方位向采样点数，采样间隔为雷达的脉冲重复时间。由于向量e_m(P)是由测量噪声e_m(u；P)所产生的，所以本发明把e_m(P)定义为测量噪声向量。需要特别指出的是，e_m(P)为一个随机向量。

由定义6可以得出测量距离历史向量R_m(P)、回波距离历史向量R(P)、测量噪声向量e_m(P)三者之间满足如下关系：

R_m(P)＝R(P)+e_m(P)

定义8、匹配度

令相位φ(u；P)＝2πR(u；P)/λ和相位φ_m(u；P)＝2πR_m(u；P)/λ，其中R(u；P)为回波距离历史，R_m(u；P)为测量距离历史，λ为雷达载波波长。利用BP算法在对回波数据进行成像的过程中，在相位补偿阶段由测量噪声e_m(u；P)所引起的相位误差Δφ_m(u；P)为：

${\mathrm{\Δ\φ}}_m(u;P)={\mathrm{\φ}}_m(u;P)-\mathrm{\φ}(u;P)=\frac{2\mathrm{\π}}{\mathrm{\λ}}(R_m(u;P)-R(u;P))=\frac{2\mathrm{\π}}{\mathrm{\λ}}{e}_m(u;P)$

假定测量噪声e_m(u；P)的构成如下：

e_m(u；P)＝γ+e′_m(u；P)

其中，u∈[-T/2,T/2]；γ为一个常数；e′_m(u；P)为一个零均值的随机过程。由此相位误差Δφ_m(u；P)可以表示为：

${\mathrm{\Δ\φ}}_m(u;P)=\frac{2\mathrm{\π}}{\mathrm{\λ}}\mathrm{\γ}+\frac{2\mathrm{\π}}{\mathrm{\λ}}{e^{\′}}_m(u;P)$

上式中，相位2πγ/λ为一个常数，该相位不会影响BP算法的聚焦效果，而真正影响BP算法聚焦效果的是相位2πe′_m(u；P)/λ。因此，为了方便起见，在本发明中，假定测量噪声e_m(u；P)为一个零均值的随机过程。

由于测量噪声e_m(u；P)是零均值的随机过程，所以测量噪声向量e_m(P)是一个零均值的随机向量，于是可以建立如下条件(详见参考文献“I.G.Cumming,F.H.Wong.Digital signal processing of synthetic aperture radar data:algorithms andimplementation.in Artech House,USA:Boston,2004:567-570”)：

$\left|{\mathrm{\Δ\φ}}_m(u_j;P)\right|=\frac{2\mathrm{\π}}{\mathrm{\λ}}\left|e_m(u_j;P)\right|\≤\frac{\mathrm{\π}}{2}$

其中，u_j(j＝1,2,…,N)为慢时间u的采样点(详见定义7)；对上式中的运算符号进行说明，如果假定x可以代表任意一个标量，|x|表示标量x的绝对值。由上式进一步可以得到如下条件：

$\left|e_m(u_j;P)\right|=|R_m(u_j;P)-R(u_j;P)|\≤\frac{\mathrm{\λ}}{4}$

在这里，令向量R_m(P)-R(P)中满足上式的分量的数目为N_α，其中R_m(P)为测量距离历史向量，R(P)为回波距离历史向量。显然，R_m(P)-R(P)的总的分量的数目为N，那么在本发明中，向量R(P)和R_m(P)之间的匹配度α被定义为：α＝N_α/N，其中α的取值范围为：α∈[0,1]。

定义9、匹配度门限值

由匹配度的定义可知，BP算法的聚焦效果随着匹配度α的减小而逐渐变差，并且当α＝1时，BP算法的聚焦效果最好。因此，根据BP算法能否聚焦，可以定义一个匹配度门限值η。当α≥η时，BP算法能够聚焦，此时本发明认为回波距离历史向量R(P)和测量距离历史向量R_m(P)是匹配的；而当α＜η时，BP算法不能聚焦，此时本发明认为R(P)和R_m(P)是不匹配的，即失配的。需要说明的是，R(P)和R_m(P)之间的匹配度的大小决定于测量噪声，如果测量噪声过大，将很有可能引起R(P)和R_m(P)的失配，从而导致BP算法不能正常聚焦。

定义10、特征距离

由定义7可知，测量距离历史向量是一个随机向量，对其去均值处理以后的距离历史向量R′_m(P)可以表示为R′_m(P)＝R_m(P)-R_m，其中R_m＝E[R_m(P)]表示R_m(P)的均值，显然也是一个随机向量。接着，建立R′_m(P)的协方差矩阵Ψ：Ψ＝E[R′_m(P)R′_m(P)^T]。令矩阵Ψ的N个特征值为并且满足条件d₁≥d₂≥…≥d_k≥…≥d_N，取出Ψ的k个最大的特征值其中参数k为整数，且1≤k＜N。根据传统的最小均方误差准则(minimummean square error，MMSE)和主成分分析准则(principal component analysis，PCA)，只要参数k取值合理，那么特征值所对应的特征向量就能够包含距离历史向量R′_m(P)中的绝大部分信息，所以本发明把特征向量定义为特征距离。详细内容可参考文献“R.Bustin,M.Payaro,D.P.Palomar,and S.Shamai.OnMMSE crossing properties and implications in parallel vector Gaussian channels.IEEE Transactions on Information Theory,2013,59(2):818-844”和“R.P.Good,D.Kost,and G.A.Cherry.Introducing a unified PCA algorithm for model size reduction.IEEE Transactions on Semiconductor Manufacturing,2010,23(2):201-209”。

定义11、特征距离子空间

令Τ＝span(ξ₁,ξ₂…,ξ_k)表示由特征距离所张成的子空间，本发明把该子空间Τ定义为特征距离子空间。在这里需要强调的是，本发明所提出的特征距离子空间的概念适用于任意体制的SAR系统。详细内容可参考文献“Y.Liao,X.Lin.Blind image restoration with eigen-face subspace.IEEE Transactions on ImageProcessing,2005,14(11):1766-1772”。

定义12、修正距离历史向量

在特征距离子空间Τ＝span(ξ₁,ξ₂…,ξ_k)中，可以最优地修正测量距离历史向量R_m(P)，修正以后的结果为：

${\hat{R}}_m\left(P\right)=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{c}_j\left(P\right)\·{\mathrm{\ξ}}_j+R_m$

其中， ${\hat{R}}_m\left(P\right)={[{\hat{R}}_m(u_1;P),{\hat{R}}_m(u_2;P),...,{\hat{R}}_m(u_N;P)]}^T$ 为向量R_m(P)的修正值，所以本发明定义其为修正距离历史向量；R_m为测量距离历史向量的均值；参数c_j(P)＝(ξ_j)^T(R_m(P)-R_m)(j＝1,2,…,k)为修正系数。

定义13、修正距离历史

修正距离历史向量 ${\hat{R}}_m\left(P\right)={[{\hat{R}}_m(u_1;P),{\hat{R}}_m(u_2;P),...,{\hat{R}}_m(u_N;P)]}^T$ 可以恢复成关于慢时间u的连续函数在本发明中，将定义为修正距离历史。

定义14、基于特征距离子空间的BP算法模型

利用修正距离历史来取代传统BP算法模型中的测量距离历史R_m(u；P)，可以得到如下新的BP算法模型，即：

$\hat{I}\left(P\right)=\underset{-T/2}{\overset{T/2}{\∫}}{S}_c(\frac{{\hat{R}}_m(u;P)}{c},u;P)\·\exp \left(j2\mathrm{\π}\frac{{\hat{R}}_m(u;P)}{\mathrm{\λ}}\right)\mathrm{du}$

其中，T表示合成孔径时间；u表示方位向慢时间；c表示电磁波的传播速度；λ表示雷达载波波长；S_c(t,u；P)表示经过距离脉冲压缩后的回波信号(详见定义3)；表示利用新的BP算法所获得的点目标P的SAR图像。由于对于R_m(u；P)的修正是在特征距离子空间中进行的，所以本发明把这种新的BP算法模型定义为基于特征距离子空间的BP算法模型。

本发明提出了一种基于特征距离子空间的SAR后向投影成像方法，它包括以下步骤：

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

为了突出本发明的普遍适用性和有效性，在本发明中，实验用的SAR系统选择为更具有一般意义的双基地SAR系统。初始化双基地SAR系统参数主要包括：平台发射机速度，记作V_T；平台接收机速度，记作V_R；平台发射机初始位置，记作P_T；平台接收机初始位置，记作P_R；雷达工作中心频率，记作f_c；雷达发射信号的距离带宽，记作B；雷达发射信号的脉冲宽度，记作T_P；雷达载波波长，记作λ；电磁波的传播速度，记作c；雷达的脉冲重复频率，记作PRF；雷达的脉冲重复时间，记作PRT；SAR的合成孔径时间，记作T；方位向采样点数，记作N_X；距离向采样点数，记作N_Y。上述参数均为双基地SAR系统的标准参数，在双基地SAR的系统设计和观测过程当中已经确定。根据双基地SAR的成像系统方案和观测方案，BP算法所需要的系统初始化参数均为已知。

步骤2、初始化观测场景目标空间参数以及产生回波距离历史：

在二维场景成像实验中，观测场景Ω中的像素点的数目为K＝X×Y，其中，X为方位向的像素点数目，Y为距离向的像素点数目，方位向的分辨率为Δ_X，距离向的分辨率为Δ_Y。观测场景Ω在SAR的成像方案设计中已经确定。由于在BP算法成像过程中，只要相位补偿精确，那么BP算法对观测场景中的每一个点目标的聚焦能力则完全相同。因此，在匹配度仿真实验中，只需要研究观测场景Ω中的任意一个点目标即可，本发明任意选取一个点目标P∈Ω作为匹配度仿真实验的研究对象，其坐标位置为P，后向散射系数为σ(P)。

对于双基地SAR，令平台发射机的位置为P_T(u)＝(x_T(u),y_T(u),z_T(u))，其中u∈[-T/2,T/2]表示方位向慢时间，由于平台发射机的位置P_T(u)的每个分坐标与平台发射机的初始位置P_T和速度V_T相对应的分坐标为线性关系，所以平台发射机的位置P_T(u)由简单计算得到。同理，可以计算得到平台接收机的位置P_R(u)＝(x_R(u),y_R(u),z_R(u))。令点目标P的位置为P＝(x,y,z)，点目标P的坐标位置可以根据需要事先进行设定；在二维场景成像实验中，根据参考坐标系的坐标系统、观测场景Ω的中心位置以及方位向的分辨率Δ_X与距离向的分辨率Δ_Y，得到观测场景Ω中的所有像素点的坐标位置。

点目标P的回波距离历史R(u；P)可以通过求解下式得到，即：

$\begin{array}{c}R(u;P)={\left|\right|P_T\left(u\right)-P\left|\right|}_2+{\left|\right|P_R\left(u\right)-P\left|\right|}_2\\ =\sqrt{{(x_T\left(u\right)-x)}^2+{(y_T\left(u\right)-y)}^2+{(z_T\left(u\right)-z)}^2+}\\ \sqrt{{(x_R\left(u\right)-x)}^2+{(y_R\left(u\right)-y)}^2+{(z_R\left(u\right)-z)}^2}\end{array}$

对上式中的运算符号进行说明，如果假定x可以代表任意一个向量，||x||₂表示向量x的l₂范数。利用上式，获得观测场景Ω中的其它点目标的回波距离历史。

步骤3、产生回波数据：

通常，SAR的发射信号为LFM信号，在不考虑天线方向图对回波信号的加权影响的条件下，点目标P的基带回波数据S(t,u；P)为：

$\begin{array}{c}S(t,u;P)=\mathrm{\σ}\left(P\right)\·\mathrm{rect}\left(\frac{t-R(u;P)/c}{T_P}\right)\·\mathrm{rect}\left(\frac{u}{T}\right)\\ \×\exp [\mathrm{j\π}\frac{B}{T_P}{(t-\frac{R(u;P)}{c})}^2-{j2\mathrm{\πf}}_c\frac{R(u;P)}{c}]\end{array}$

其中，t表示距离向快时间；u表示方位向慢时间；σ(P)表示点目标P的后向散射系数；B表示雷达发射信号的距离带宽；T_P表示雷达发射信号的脉冲宽度；T表示合成孔径时间；f_c表示雷达工作中心频率；c表示电磁波的传播速度；rect(u/T)表示慢时间u的窗函数；R(u；P)表示回波距离历史。基带回波数据S(t,u；P)采用传统的合成孔径雷达回波仿真方法产生得到。观测场景Ω中的其它点目标的基带回波数据的产生方法同点目标P。

步骤4、对SAR基带回波数据进行距离脉冲压缩：

采用传统的合成孔径雷达距离脉冲压缩方法，对步骤3中所得到的SAR基带回波数据S(t,u；P)进行距离脉冲压缩处理。得到经过距离脉冲压缩以后的数据S_c(t,u；P)为：

$\begin{array}{c}{S}_c(t,u;P)=\mathrm{\σ}\left(P\right)\·\mathrm{rect}\left(\frac{t-R(u;P)/c}{T_P}\right)\·\mathrm{rect}\left(\frac{u}{T}\right)\\ \×\sin c\left[\mathrm{\πB}(t-\frac{R(u;P)}{c})\right]\·\exp (-j2\mathrm{\π}\frac{R(u;P)}{\mathrm{\λ}})\end{array}$

其中，t表示距离向快时间；u表示方位向慢时间；σ(P)表示点目标P的后向散射系数；B表示雷达发射信号的距离带宽；T_P表示雷达发射信号的脉冲宽度；T表示合成孔径时间；c表示电磁波的传播速度；λ表示雷达载波波长；rect(u/T)表示慢时间u的窗函数；R(u；P)表示回波距离历史。同理，利用上面求解经过距离脉冲压缩以后的数据S_c(t,u；P)的公式，得到观测场景Ω中的其它点目标经过距离脉冲压缩后的回波数据。

步骤5、测量噪声：

测量噪声e_m(u；P)包括两种类型：第一种类型服从高斯分布N(μ,δ²)，其中μ表示测量噪声的均值，δ²表示测量噪声的方差；第二种类型服从均匀分布U(a,b)，其中a表示测量噪声取值的下限，b表示测量噪声取值的上限。

步骤6、获取测量距离历史：

由步骤2获得回波距离历史R(u；P)，同时由步骤5获得测量噪声e_m(u；P)，通过求解下式，获得点目标P的测量距离历史R_m(u；P)，即：

R_m(u；P)＝R(u；P)+e_m(u；P)

利用上式，获得观测场景Ω中的其它点目标的测量距离历史。

步骤7、利用传统BP算法获取SAR图像：

利用传统BP算法，对距离脉冲压缩后的回波数据通过相关积累，获得点目标P的SAR图像I(P)，即：

$I\left(P\right)=\underset{-T/2}{\overset{T/2}{\∫}}{S}_c(\frac{R_m(u;P)}{c},u;P)\·\exp \left(j2\mathrm{\π}\frac{R_m(u;P)}{\mathrm{\λ}}\right)\mathrm{du}$

其中，T表示合成孔径时间；u表示方位向慢时间；c表示电磁波的传播速度；λ表示雷达载波波长；R_m(u；P)表示测量距离历史；S_c(t,u；P)表示经过距离脉冲压缩后的回波数据(详见步骤4)。利用上式，获得观测场景Ω中的其它点目标的SAR图像。

步骤8、获取回波距离历史向量和测量噪声向量：

由步骤2获得回波距离历史R(u；P)，针对回波距离历史R(u；P)，通过在[-T/2,T/2]区间内对慢时间u进行均匀采样，得到回波距离历史向量R(P)＝[R(u₁；P),R(u₂；P),…,R(u_N；P)]^T，其中是u的采样点，N表示方位向采样点数(N＝N_X)，采样间隔为雷达的脉冲重复时间PRT。

由步骤5获得测量噪声e_m(u；P)，针对测量噪声e_m(u；P)，通过在[-T/2,T/2]区间内对慢时间u进行均匀采样，得到测量噪声向量e_m(P)＝[e_m(u₁；P),e_m(u₂；P),…,e_m(u_N；P)]^T，其中是u的采样点，N表示方位向采样点数(N＝N_X)，采样间隔为雷达的脉冲重复时间PRT。

步骤9、获取测量距离历史向量：

由步骤8获得回波距离历史向量R(P)以及测量噪声向量e_m(P)，令点目标P的测量距离历史向量为R_m(P)＝[R_m(u₁；P),R_m(u₂；P),…,R_m(u_N；P)]^T，通过求解下式，获得R_m(P)，即：

R_m(P)＝R(P)+e_m(P)

观测场景Ω中的其它点目标的测量距离历史向量的获取方法同点目标P。

步骤10、测量距离历史向量去均值处理：

对测量距离历史向量R_m(P)进行去均值处理，即：

R′_m(P)＝R_m(P)-R_m

上式中，R_m＝E[R_m(P)]表示R_m(P)的均值；R′_m(P)表示去均值以后的测量距离历史向量。

步骤11、求解特征距离子空间：

建立去均值以后的测量距离历史向量R′_m(P)的协方差矩阵Ψ，用下式所示：

Ψ＝E[R′_m(P)R′_m(P)^T]

首先求出矩阵Ψ的N个特征值同时把这N个特征值按照由大到小的次序进行排列，即d₁≥d₂≥…≥d_k≥…≥d_N。然后取出矩阵Ψ的k个最大的特征值其中参数k为整数，且1≤k＜N。特征值所对应的特征向量就是特征距离，从而得到特征距离子空间Τ＝span(ξ₁,ξ₂,…,ξ_k)。

步骤12、求解修正距离历史向量和修正距离历史：

在特征距离子空间Τ＝span(ξ₁,ξ₂,…,ξ_k)中，采用下面的公式修正测量距离历史向量R_m(P)，即：

${\hat{R}}_m\left(P\right)=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{c}_j\left(P\right)\·{\mathrm{\ξ}}_j+R_m$

其中， ${\hat{R}}_m\left(P\right)={[{\hat{R}}_m(u_1;P),{\hat{R}}_m(u_2;P),...,{\hat{R}}_m(u_N;P)]}^T$ 为向量R_m(P)的修正值，即就是所求的修正距离历史向量；R_m为测量距离历史向量的均值；参数c_j(P)＝(ξ_j)^T(R_m(P)-R_m)(j＝1,2,…,k)为修正系数。

利用传统的插值运算把向量恢复成关于慢时间u的连续函数就是所求的修正距离历史。观测场景Ω中的其它点目标的修正距离历史向量和修正距离历史的获取方法同点目标P。

步骤13、求解补偿的相位：

利用修正距离历史计算出本发明所提出的新的BP算法在成像过程中需要补偿的相位 ${\widehat{\mathrm{\φ}}}_m(u;P):{\widehat{\mathrm{\φ}}}_m(u;P)=2\mathrm{\π}{\hat{R}}_m(u;P)/\mathrm{\λ},$ 其中λ为雷达载波波长。

步骤14、利用新的BP算法获取SAR图像：

求出需要补偿的相位以后，通过采用传统的脉冲相关积累方法对距离脉冲压缩后的回波数据进行相关积累，获得点目标P的SAR图像即：

$\hat{I}\left(P\right)=\underset{-T/2}{\overset{T/2}{\∫}}{S}_c(\frac{{\hat{R}}_m(u;P)}{c},u;P)\·\exp \left(j{\widehat{\mathrm{\φ}}}_m(u;P)\right)\mathrm{du}$

其中，T表示合成孔径时间；u表示方位向的慢时间；c表示电磁波的传播速度；表示修正距离历史；S_c(t,u；P)表示经过距离脉冲压缩后的回波数据(详见步骤4)；表示利用新的BP算法所获得的点目标P的SAR图像。利用求取点目标P的SAR图像的上述公式，获得观测场景Ω中的其它点目标的SAR图像。

本发明的创新点在于针对如果测量噪声过大，导致传统BP算法不能正常聚焦的情况，提供了一种新的BP算法，即基于特征距离子空间的BP成像方法。新的BP算法首先以距离历史向量为基础，从BP算法能否聚焦的角度提出了距离历史向量匹配度的概念。然后基于传统的最小均方误差准则和主成分分析准则，在SAR成像领域，新的BP算法提出了特征距离子空间的概念。由于新的BP算法在特征距离子空间中最优地修正了测量距离历史向量，所以很好地抑制了测量噪声，提高了距离历史向量的匹配度，从而对补偿的相位进行了准确地估计，最终获得了良好的聚焦效果。

本发明的优点在于首先新的BP算法理论完善，实验结果稳定。其次，由于新的BP算法适用于任意体制的SAR二维成像和三维成像，所以新的BP算法具有很强的可推广性。最后，由于在特征距离子空间中修正测量距离历史向量是一种线性运算，所以新的BP算法具有很强的可操作性。

附图说明：

图1为本发明对双基地SAR回波数据进行成像处理的流程图。

图2为本发明在具体实施实验过程中所采用的双基地SAR系统的仿真参数列表

其中，m表示米；Km表示千米；s表示秒；m/s表示米每秒；MHz表示兆赫兹；GHz表示吉赫兹。

图3为本发明在测量噪声服从高斯分布时，利用蒙特-卡罗实验所计算出的匹配度

其中，“times of experiments”表示实验次数，实验次数为100次；“matchingratio”表示匹配度；α表示测量距离历史向量R_m(P)和回波距离历史向量R(P)之间的匹配度；β表示修正距离历史向量和回波距离历史向量R(P)之间的匹配度；测量噪声服从高斯分布N(0m,3×10^-4m²)。

图4为本发明在测量噪声服从均匀分布时，利用蒙特-卡罗实验所计算出的匹配度

其中，“times of experiments”表示实验次数，实验次数为100次；“matchingratio”表示匹配度；α表示测量距离历史向量R_m(P)和回波距离历史向量R(P)之间的匹配度；β表示修正距离历史向量和回波距离历史向量R(P)之间的匹配度；测量噪声服从均匀分布U(-0.02m,0.02m)。

图5为本发明在测量噪声服从零均值高斯分布时，利用传统BP算法对二维观测场景所成的SAR图像

其中，“Range”表示距离向，“Azimuth”表示方位向；测量噪声服从零均值高斯分布N(0m,3×10^-4m²)。

图6为本发明在测量噪声服从零均值高斯分布时，利用本发明所提出的新的BP算法对二维观测场景所成的SAR图像

其中，“Range”表示距离向，“Azimuth”表示方位向；测量噪声服从零均值高斯分布N(0m,3×10^-4m²)。

图7为本发明在测量噪声服从零均值均匀分布时，利用传统BP算法对二维观测场景所成的SAR图像

其中，“Range”表示距离向，“Azimuth”表示方位向；测量噪声服从零均值均匀分布U(-0.02m,0.02m)。

图8为本发明在测量噪声服从零均值均匀分布时，利用本发明所提出的新的BP算法对二维观测场景所成的SAR图像

其中，“Range”表示距离向，“Azimuth”表示方位向；测量噪声服从零均值均匀分布U(-0.02m,0.02m)。

图9为本发明在测量噪声服从非零均值高斯分布时，利用传统BP算法对二维观测场景所成的SAR图像

其中，“Range”表示距离向，“Azimuth”表示方位向；测量噪声服从非零均值高斯分布N(0.02m,3×10^-4m²)。

图10为本发明在测量噪声服从非零均值高斯分布时，利用本发明所提出的新的BP算法对二维观测场景所成的SAR图像

其中，“Range”表示距离向，“Azimuth”表示方位向；测量噪声服从非零均值高斯分布N(0.02m,3×10^-4m²)。

图11为本发明在测量噪声服从非零均值均匀分布时，利用传统BP算法对二维观测场景所成的SAR图像

其中，“Range”表示距离向，“Azimuth”表示方位向；测量噪声服从非零均值均匀分布U(0m,0.04m)。

图12为本发明在测量噪声服从非零均值均匀分布时，利用本发明所提出的新的BP算法对二维观测场景所成的SAR图像

其中，“Range”表示距离向，“Azimuth”表示方位向；测量噪声服从非零均值均匀分布U(0m,0.04m)。

具体实施方式：

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

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

在本发明中，初始化双基地SAR系统参数主要包括：平台发射机速度V_T＝(200m/s,0m/s,0m/s)；平台接收机速度V_R＝(80m/s,100m/s,0m/s)；平台发射机初始位置P_T＝(-3Km,-20Km,20Km)；平台接收机初始位置P_R＝(0Km,-5Km,8Km)；雷达工作中心频率f_c＝10GHz；雷达发射信号的距离带宽B＝150MHz；雷达发射信号的脉冲宽度T_P＝1×10^-6s；雷达载波波长λ＝0.03m；电磁波的传播速度c＝3×10⁸m/s；雷达的脉冲重复频率PRF＝400Hz；雷达的脉冲重复时间PRT＝2.5×10^-3s；SAR的合成孔径时间T＝1.28s；方位向采样点数N_X＝512；距离向采样点数N_Y＝1024。上述参数均为双基地SAR系统的标准参数，在双基地SAR的系统设计和观测过程当中已经确定。根据双基地SAR的成像系统方案和观测方案，BP算法所需要的系统初始化参数均为已知。

步骤2、初始化观测场景目标空间参数以及产生回波距离历史：

在匹配度仿真实验中，任意选取一个点目标P∈Ω作为研究对象，其坐标位置为P＝(50m,50m,0m)，后向散射系数为σ(P)＝1。需要说明的是，点目标P的后向散射系数的大小与匹配度无关，所以为了方便起见，在匹配度仿真实验中，点目标P的后向散射系数设定为σ(P)＝1。在二维场景成像实验中，观测场景Ω中的像素点的数目为K＝X×Y＝256×256，其中，方位向的像素点数目为X＝256，距离向的像素点数目为Y＝256，方位向的分辨率为Δ_X＝3m，距离向的分辨率为Δ_Y＝3m。

对于双基地SAR，令平台发射机的位置为P_T(u)＝(x_T(u),y_T(u),z_T(u))，其中u∈[-T/2,T/2]表示方位向慢时间，由于平台发射机的位置P_T(u)的每个分坐标与平台发射机的初始位置P_T和速度V_T相对应的分坐标为线性关系，所以平台发射机的位置P_T(u)可以通过简单的计算得到。同理，可以计算出平台接收机的位置P_R(u)＝(x_R(u),y_R(u),z_R(u))。接下来，令点目标P的位置为P＝(x,y,z)，需要说明的是：在本发明中，对于匹配度仿真实验，点目标P的坐标位置分别为x＝50m、y＝50m、z＝0m；对于二维场景成像实验，根据参考坐标系的坐标系统、观测场景Ω的中心位置以及方位向的分辨率Δ_X与距离向的分辨率Δ_Y，可以得到观测场景Ω中的所有像素点的坐标位置。

点目标P的回波距离历史R(u；P)可以通过求解下式来得到，即：

$\begin{array}{c}R(u;P)={\left|\right|P_T\left(u\right)-P\left|\right|}_2+{\left|\right|P_R\left(u\right)-P\left|\right|}_2\\ =\sqrt{{(x_T\left(u\right)-x)}^2+{(y_T\left(u\right)-y)}^2+{(z_T\left(u\right)-z)}^2+}\\ \sqrt{{(x_R\left(u\right)-x)}^2+{(y_R\left(u\right)-y)}^2+{(z_R\left(u\right)-z)}^2}\end{array}$

对上式中的运算符号进行说明，如果假定x可以代表任意一个向量，||x||₂表示向量x的l₂范数。利用上式，可以获得观测场景Ω中的其它点目标的回波距离历史。

在这里需要特别说明的是，在现实中，回波距离历史R(u；P)是未知的，也是无法获取的。由于本发明主要采用仿真实验的方法进行计算和验证，因此在该步骤中，之所以要事先获得回波距离历史R(u；P)，是为了在下面的步骤中产生测量距离历史。然而，在实际中，利用BP算法对回波数据真正进行聚焦时，在相位补偿阶段所用到的距离历史是测量距离历史而不是回波距离历史。

步骤3、产生回波数据：

采用传统的合成孔径雷达回波仿真方法产生SAR的回波数据。点目标P的基带回波数据S(t,u；P)为：

$\begin{array}{c}S(t,u;P)=\mathrm{\σ}\left(P\right)\·\mathrm{rect}\left(\frac{t-R(u;P)/c}{T_P}\right)\·\mathrm{rect}\left(\frac{u}{T}\right)\\ \×\exp [\mathrm{j\π}\frac{B}{T_P}{(t-\frac{R(u;P)}{c})}^2-{j2\mathrm{\πf}}_c\frac{R(u;P)}{c}]\end{array}$

其中，t表示距离向快时间；u表示方位向慢时间；σ(P)表示点目标P的后向散射系数；B表示雷达发射信号的距离带宽；T_P表示雷达发射信号的脉冲宽度；T表示合成孔径时间；f_c表示雷达工作中心频率；c表示电磁波的传播速度；rect(u/T)表示慢时间u的窗函数；R(u；P)表示回波距离历史。利用上式，可以产生观测场景Ω中的其它点目标的基带回波数据。

步骤4、对SAR基带回波数据进行距离脉冲压缩：

采用传统的合成孔径雷达距离脉冲压缩方法，对步骤3中所得到的SAR基带回波数据S(t,u；P)进行距离脉冲压缩处理。经过距离脉冲压缩以后的数据S_c(t,u；P)为：

$\begin{array}{c}{S}_c(t,u;P)=\mathrm{\σ}\left(P\right)\·\mathrm{rect}\left(\frac{t-R(u;P)/c}{T_P}\right)\·\mathrm{rect}\left(\frac{u}{T}\right)\\ \×\sin c\left[\mathrm{\πB}(t-\frac{R(u;P)}{c})\right]\·\exp (-j2\mathrm{\π}\frac{R(u;P)}{\mathrm{\λ}})\end{array}$

其中，t表示距离向快时间；u表示方位向慢时间；σ(P)表示点目标P的后向散射系数；B表示雷达发射信号的距离带宽；T_P表示雷达发射信号的脉冲宽度；T表示合成孔径时间；c表示电磁波的传播速度；λ表示雷达载波波长；rect(u/T)表示慢时间u的窗函数；R(u；P)表示回波距离历史。利用上式，可以产生观测场景Ω中的其它点目标经过距离脉冲压缩后的回波数据。

步骤5、测量噪声：

在本发明的仿真试验中，以测量噪声的均值是否为零，把测量噪声分为两组。第一组测量噪声e_m(u；P)分别服从零均值的高斯分布N(0m,3×10^-4m²)和零均值的均匀分布U(-0.02m,0.02m)；第二组测量噪声e_m(u；P)分别服从非零均值的高斯分布N(0.02m,3×10^-4m²)和非零均值的均匀分布U(0m,0.04m)。

步骤6、获取测量距离历史：

由步骤2可以获得回波距离历史R(u；P)，同时由步骤5可以获得测量噪声e_m(u；P)，那么通过求解下式，便可获得点目标P的测量距离历史R_m(u；P)，即：

R_m(u；P)＝R(u；P)+e_m(u；P)

利用上式，可以获得观测场景Ω中的其它点目标的测量距离历史。

步骤7、利用传统BP算法获取SAR图像：

利用传统BP算法，通过对距离脉冲压缩后的回波数据进行相关积累，便可获得点目标P的SAR图像I(P)，即：

$I\left(P\right)=\underset{-T/2}{\overset{T/2}{\∫}}{S}_c(\frac{R_m(u;P)}{c},u;P)\·\exp \left(j2\mathrm{\π}\frac{R_m(u;P)}{\mathrm{\λ}}\right)\mathrm{du}$

其中，T表示合成孔径时间；u表示方位向慢时间；c表示电磁波的传播速度；λ表示雷达载波波长；R_m(u；P)表示测量距离历史；S_c(t,u；P)表示经过距离脉冲压缩后的回波数据(详见步骤4)。利用上式，可以获得观测场景Ω中的其它点目标的SAR图像。

步骤8、获取回波距离历史向量和测量噪声向量：

由步骤2可以获得回波距离历史R(u；P)。针对回波距离历史R(u；P)，通过在[-T/2,T/2]区间内对慢时间u进行均匀采样，可以得到回波距离历史向量R(P)＝[R(u₁；P),R(u₂；P),…,R(u_N；P)]^T，其中是u的采样点，N表示方位向采样点数(N＝N_X)，采样间隔为雷达的脉冲重复时间PRT。

由步骤5可以获得测量噪声e_m(u；P)。针对测量噪声e_m(u；P)，通过在[-T/2,T/2]区间内对慢时间u进行均匀采样，可以得到测量噪声向量e_m(P)＝[e_m(u₁；P),e_m(u₂；P),…,e_m(u_N；P)]^T，其中是u的采样点，N表示方位向采样点数(N＝N_X)，采样间隔为雷达的脉冲重复时间PRT。

步骤9、获取测量距离历史向量：

由步骤8可以获得回波距离历史向量R(P)以及测量噪声向量e_m(P)，令点目标P的测量距离历史向量为R_m(P)＝[R_m(u₁；P),R_m(u₂；P),…,R_m(u_N；P)]^T，那么通过求解下式，便可获得R_m(P)，即：

R_m(P)＝R(P)+e_m(P)

观测场景Ω中的其它点目标的测量距离历史向量的获取方法同点目标P。

步骤10、求解测量距离历史向量均值的估计值：

理论上，测量距离历史向量R_m(P)的均值为R_m＝E[R_m(P)]。实际中，理论值R_m很难得到，所以可以利用相关统计学知识对其进行有效地估计。由步骤2可知，观测场景Ω中共有K个像素点，即那么第j(j＝1,2,…,K)个像素点P_j的测量距离历史向量可以表示为R_m(P_j)。因此，可以获得测量距离历史向量均值R_m的估计值即：

${\hat{R}}_m=\frac{1}{K}\underset{j=1}{\overset{K}{\mathrm{\Σ}}}{R}_m\left(P_j\right)$

步骤11、求解去均值后的测量距离历史向量协方差矩阵的估计值：

理论上，去均值以后的测量距离历史向量为R′_m(P)＝R_m(P)-R_m，其中R_m(P)为测量距离历史向量，R_m＝E[R_m(P)]为R_m(P)的均值，由R′_m(P)所产生的协方差矩阵为Ψ＝E[R′_m(P)R′_m(P)^T]。实际中，理论值Ψ很难得到，所以可以利用相关统计学知识对其进行有效地估计。利用步骤10的结果，结合相关统计学知识，可以获得协方差矩阵Ψ的估计值即：

$\widehat{\mathrm{\Ψ}}=\frac{1}{K}\underset{j=1}{\overset{K}{\mathrm{\Σ}}}(R_m\left(P_j\right)-{\hat{R}}_m){(R_m\left(P_j\right)-{\hat{R}}_m)}^T$

其中，K为观测场景Ω中的像素点数目；R_m(P_j)为Ω中的第j(j＝1,2,…,K)个像素点P_j的测量距离历史向量，为测量距离历史向量均值的估计值。

步骤12、求解特征距离子空间：

首先求出协方差矩阵的估计值的N个特征值同时把这N个特征值按照由大到小的次序进行排列，即然后取出矩阵的k个最大的特征值其中参数k为整数，且1≤k＜N。那么特征值所对应的特征向量就是特征距离，从而得到特征距离子空间由于本发明是利用BP算法对二维观测场景进行成像的，那么参数k的最优取值为k＝2，所以本发明中所用到的特征距离子空间为 $\hat{T}=\mathrm{span}({\widehat{\mathrm{\ξ}}}_1,{\widehat{\mathrm{\ξ}}}_2).$

步骤13、求解修正距离历史向量和修正距离历史：

在特征距离子空间中，本发明可以最优地修正测量距离历史向量R_m(P)，修正以后的结果为：

${\hat{R}}_m\left(P\right)=\underset{j=1}{\overset{2}{\mathrm{\Σ}}}{\hat{c}}_j\left(P\right)\·{\widehat{\mathrm{\ξ}}}_j+{\hat{R}}_m$

其中， ${\hat{R}}_m\left(P\right)={[{\hat{R}}_m(u_1;P),{\hat{R}}_m(u_2;P),...,{\hat{R}}_m(u_N;P)]}^T$ 为向量R_m(P)的修正值，即就是所求的修正距离历史向量；参数 ${\hat{c}}_j\left(P\right)={\left({\widehat{\mathrm{\ξ}}}_j\right)}^T(R_m\left(P\right)-{\hat{R}}_m)(j=\mathrm{1,2})$ 为修正系数；为测量距离历史向量均值的估计值。可以利用插值运算把向量恢复成关于慢时间u的连续函数就是所求的修正距离历史。观测场景Ω中的其它点目标的修正距离历史向量和修正距离历史的获取方法同点目标P。

步骤14、匹配度计算和聚焦效果分析：

根据定义8，计算出测量距离历史向量R_m(P)和回波距离历史向量R(P)之间的匹配度α，同时计算出修正距离历史向量和回波距离历史向量R(P)之间的匹配度β。如果由于测量噪声过大，从而导致匹配度α小于匹配度门限值η，那么传统BP算法就不能正常聚焦了，在这种情况下，如果匹配度β大于匹配度门限值η，那么就可以采用本发明所提出的新的BP算法进行聚焦。

由于测量距离历史向量R_m(P)是一个随机向量，所以R_m(P)的每一次具体的观测值都不相同，所以为了充分验证本发明所提出的新的BP算法的正确性和有效性，附图中的图3和图4都通过100次蒙特-卡罗实验，分别绘制出在不同类型测量噪声的情况下匹配度α和β的曲线。

在这里需要特别说明的是，本发明通过大量的仿真实验，发现匹配度门限值η的取值为η∈[0.6,0.8]。

观察附图中的图3和图4可以看出，当测量噪声过大的时候，匹配度α的每次实验的计算结果均远远小于匹配度门限值η∈[0.6,0.8]，此时BP算法不能聚焦。而匹配度β的每次实验的计算结果均等于1，此时BP算法能够进行精确地聚焦。由此可以得出，本发明所提出的新的BP算法确实可以有效地抑制测量噪声。

步骤15、求解补偿的相位：

利用修正距离历史可以求出本发明所提出的新的BP算法在成像过程中需要补偿的相位其中λ为雷达载波波长。

步骤16、利用新的BP算法获取SAR图像和成像结果分析：

求出需要补偿的相位以后，利用本发明所提出的新的BP算法，通过对距离脉冲压缩后的回波数据进行相关积累，便可获得点目标P的SAR图像即：

$\hat{I}\left(P\right)=\underset{-T/2}{\overset{T/2}{\∫}}{S}_c(\frac{{\hat{R}}_m(u;P)}{c},u;P)\·\exp \left(j{\widehat{\mathrm{\φ}}}_m(u;P)\right)\mathrm{du}$

其中，T表示合成孔径时间；u表示方位向的慢时间；c表示电磁波的传播速度；表示修正距离历史；S_c(t,u；P)表示经过距离脉冲压缩后的回波数据(详见步骤4)；表示利用新的BP算法所获得的点目标P的SAR图像。利用上式，可以获得观测场景Ω中的其它点目标的SAR图像。

对比附图中的图5和图6以及图7和图8，很显然，当测量噪声过大的时候，传统BP算法已经不能聚焦了，然而本发明所提出的新的BP算法却依然能够进行精确地聚焦。由图6和图8还可以看出，对于不同类型的测量噪声，新的BP算法的成像效果基本相当，其原因阐述如下：由附图中的图3和图4可以很明显地看出，无论测量噪声是服从高斯分布N(0m,3×10^-4m²)还是服从均匀分布U(-0.02m,0.02m)，修正距离历史向量和回波距离历史向量R(P)之间的匹配度β都等于1，因此在这种情况下，BP算法能够进行精确地聚焦，且基本无差别。

观察附图中的图10和图12可以看出，当测量噪声的均值为一个非零的常数时，本发明所提出的新的BP算法依然能够进行精确地聚焦，其原因阐述如下：在图10和图12中，测量噪声的均值等于0.02m，其为一个常数；由定义8可知，此时参数γ＝0.02m，γ所对应的相位2πγ/λ＝4π/3为一个常数相位，该相位对BP算法的聚焦效果基本上无影响。

综上所述，本发明所提出的新的BP算法，确实具有优良地抑制测量噪声的能力。

Claims (1)

1.一种基于特征距离子空间的SAR后向投影成像方法，其特征是它包括以下步骤：

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

初始化双基地SAR系统参数主要包括：平台发射机速度，记作V_T；平台接收机速度，记作V_R；平台发射机初始位置，记作P_T；平台接收机初始位置，记作P_R；雷达工作中心频率，记作f_c；雷达发射信号的距离带宽，记作B；雷达发射信号的脉冲宽度，记作T_P；雷达载波波长，记作λ；电磁波的传播速度，记作c；雷达的脉冲重复频率，记作PRF；雷达的脉冲重复时间，记作PRT；SAR的合成孔径时间，记作T；方位向采样点数，记作N_X；距离向采样点数，记作N_Y；上述参数均为双基地SAR系统的标准参数，在双基地SAR的系统设计和观测过程当中已经确定；根据双基地SAR的成像系统方案和观测方案，BP算法所需要的系统初始化参数均为已知；

步骤2、初始化观测场景目标空间参数以及产生回波距离历史：

在二维场景成像中，观测场景Ω中的像素点的数目为K＝X×Y，其中，X为方位向的像素点数目，Y为距离向的像素点数目，方位向的分辨率为Δ_X，距离向的分辨率为Δ_Y；观测场景Ω在SAR的成像方案设计中已经确定；任意选取一个点目标P∈Ω作为研究对象，其坐标位置为P，后向散射系数为σ(P)；

对于双基地SAR，令平台发射机的位置为P_T(u)＝(x_T(u),y_T(u),z_T(u))，其中u∈[-T/2,T/2]表示方位向慢时间，由于平台发射机的位置P_T(u)的每个分坐标与平台发射机的初始位置P_T和速度V_T相对应的分坐标为线性关系，所以平台发射机的位置P_T(u)由简单计算得到；同理，得到平台接收机的位置P_R(u)＝(x_R(u),y_R(u),z_R(u))；令点目标P的位置为P＝(x,y,z)，点目标P的坐标位置可以根据需要事先进行设定；在二维场景成像实验中，根据参考坐标系的坐标系统、观测场景Ω的中心位置以及方位向的分辨率Δ_X与距离向的分辨率Δ_Y，得到观测场景Ω中的所有像素点的坐标位置；

点目标P的回波距离历史R(u；P)可以通过求解下式得到，即：

$\begin{array}{c}R(u;P)={\left|\right|P_T\left(u\right)-P\left|\right|}_2+{\left|\right|P_R\left(u\right)-P\left|\right|}_2\\ =\sqrt{{(x_T\left(u\right)-x)}^2+{(y_T\left(u\right)-y)}^2+{(z_T\left(u\right)-z)}^2+}\\ \sqrt{{(x_R\left(u\right)-x)}^2+{(y_R\left(u\right)-y)}^2+{(z_R\left(u\right)-z)}^2}\end{array}$

对上式中的运算符号进行说明，如果假定x可以代表任意一个向量，||x||₂表示向量x的l₂范数；利用上式，获得观测场景Ω中的其它点目标的回波距离历史；

步骤3、产生回波数据：

通常，SAR的发射信号为LFM信号，在不考虑天线方向图对回波信号的加权影响的条件下，点目标P的基带回波数据S(t,u；P)为：

$\begin{array}{c}S(t,u;P)=\mathrm{\σ}\left(P\right)\·\mathrm{rect}\left(\frac{t-R(u;P)/c}{T_P}\right)\·\mathrm{rect}\left(\frac{u}{T}\right)\\ \×\exp [\mathrm{j\π}\frac{B}{T_P}{(t-\frac{R(u;P)}{c})}^2-{j2\mathrm{\πf}}_c\frac{R(u;P)}{c}]\end{array}$

其中，t表示距离向快时间；u表示方位向慢时间；σ(P)表示点目标P的后向散射系数；B表示雷达发射信号的距离带宽；T_P表示雷达发射信号的脉冲宽度；T表示合成孔径时间；f_c表示雷达工作中心频率；c表示电磁波的传播速度；rect(u/T)表示慢时间u的窗函数；R(u；P)表示回波距离历史；基带回波数据S(t,u；P)采用传统的合成孔径雷达回波仿真方法产生得到；观测场景Ω中的其它点目标的基带回波数据的产生方法同点目标P；

步骤4、对SAR基带回波数据进行距离脉冲压缩：

采用传统的合成孔径雷达距离脉冲压缩方法，对步骤3中所得到的SAR基带回波数据S(t,u；P)进行距离脉冲压缩处理；得到经过距离脉冲压缩以后的数据S_c(t,u；P)为：

$\begin{array}{c}{S}_c(t,u;P)=\mathrm{\σ}\left(P\right)\·\mathrm{rect}\left(\frac{t-R(u;P)/c}{T_P}\right)\·\mathrm{rect}\left(\frac{u}{T}\right)\\ \×\sin c\left[\mathrm{\πB}(t-\frac{R(u;P)}{c})\right]\·\exp (-j2\mathrm{\π}\frac{R(u;P)}{\mathrm{\λ}})\end{array}$

其中，t表示距离向快时间；u表示方位向慢时间；σ(P)表示点目标P的后向散射系数；B表示雷达发射信号的距离带宽；T_P表示雷达发射信号的脉冲宽度；T表示合成孔径时间；c表示电磁波的传播速度；λ表示雷达载波波长；rect(u/T)表示慢时间u的窗函数；R(u；P)表示回波距离历史；同理，利用上面求解经过距离脉冲压缩以后的数据S_c(t,u；P)的公式，得到观测场景Ω中的其它点目标经过距离脉冲压缩后的回波数据；

步骤5、测量噪声：

测量噪声e_m(u；P)包括两种类型：第一种类型服从高斯分布N(μ,δ²)，其中μ表示测量噪声的均值，δ²表示测量噪声的方差；第二种类型服从均匀分布U(a,b)，其中a表示测量噪声取值的下限，b表示测量噪声取值的上限；

步骤6、获取测量距离历史：

由步骤2获得回波距离历史R(u；P)，同时由步骤5获得测量噪声e_m(u；P)，通过求解下式，获得点目标P的测量距离历史R_m(u；P)，即：

R_m(u；P)＝R(u；P)+e_m(u；P)

利用上式，获得观测场景Ω中的其它点目标的测量距离历史；

步骤7、利用传统BP算法获取SAR图像：

利用传统BP算法，对距离脉冲压缩后的回波数据通过相关积累，获得点目标P的SAR图像I(P)，即：

$I\left(P\right)=\underset{-T/2}{\overset{T/2}{\∫}}{S}_c(\frac{R_m(u;P)}{c},u;P)\·\exp \left(j2\mathrm{\π}\frac{R_m(u;P)}{\mathrm{\λ}}\right)\mathrm{du}$

其中，T表示合成孔径时间；u表示方位向慢时间；c表示电磁波的传播速度；λ表示雷达载波波长；R_m(u；P)表示测量距离历史；S_c(t,u；P)表示由步骤4得到的经过距离脉冲压缩后的回波数据；利用上式，获得观测场景Ω中的其它点目标的SAR图像；

步骤8、获取回波距离历史向量和测量噪声向量：

由步骤2获得回波距离历史R(u；P)，针对回波距离历史R(u；P)，通过在[-T/2,T/2]区间内对慢时间u进行均匀采样，得到回波距离历史向量R(P)＝[R(u₁；P),R(u₂；P),…,R(u_N；P)]^T，其中是u的采样点，N表示方位向采样点数(N＝N_X)，采样间隔为雷达的脉冲重复时间PRT；

由步骤5获得测量噪声e_m(u；P)，针对测量噪声e_m(u；P)，通过在[-T/2,T/2]区间内对慢时间u进行均匀采样，得到测量噪声向量e_m(P)＝[e_m(u₁；P),e_m(u₂；P),…,e_m(u_N；P)]^T，其中是u的采样点，N表示方位向采样点数(N＝N_X)，采样间隔为雷达的脉冲重复时间PRT；

步骤9、获取测量距离历史向量：

由步骤8获得回波距离历史向量R(P)以及测量噪声向量e_m(P)，令点目标P的测量距离历史向量为R_m(P)＝[R_m(u₁；P),R_m(u₂；P),…,R_m(u_N；P)]^T，通过求解下式，获得R_m(P)，即：

R_m(P)＝R(P)+e_m(P)

观测场景Ω中的其它点目标的测量距离历史向量的获取方法同点目标P；

步骤10、测量距离历史向量去均值处理：

对测量距离历史向量R_m(P)进行去均值处理，即：

R′_m(P)＝R_m(P)-R_m

上式中，R_m＝E[R_m(P)]表示R_m(P)的均值；R′_m(P)表示去均值以后的测量距离历史向量；

步骤11、求解特征距离子空间：

建立去均值以后的测量距离历史向量R′_m(P)的协方差矩阵Ψ，用下式所示：

Ψ＝E[R′_m(P)R′_m(P)T]

首先求出矩阵Ψ的N个特征值同时把这N个特征值按照由大到小的次序进行排列，即d₁≥d₂≥…≥d_k≥…≥d_N；然后取出矩阵Ψ的k个最大的特征值其中参数k为整数，且1≤k＜N；特征值所对应的特征向量就是特征距离，从而得到特征距离子空间Τ＝span(ξ₁,ξ₂,…,ξ_k)；

步骤12、求解修正距离历史向量和修正距离历史：

在特征距离子空间Τ＝span(ξ₁,ξ₂,…,ξ_k)中，采用下面的公式修正测量距离历史向量R_m(P)，即：

${\hat{R}}_m\left(P\right)=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{c}_j\left(P\right)\·{\mathrm{\ξ}}_j+R_m$

其中， ${\hat{R}}_m\left(P\right)={[{\hat{R}}_m(u_1;P),{\hat{R}}_m(u_2;P),...,{\hat{R}}_m(u_N;P)]}^T$ 为向量R_m(P)的修正值，即就是所求的修正距离历史向量；R_m为测量距离历史向量的均值；参数c_j(P)＝(ξ_j)^T(R_m(P)-R_m)(j＝1,2,…,k)为修正系数；

利用传统的插值运算把向量恢复成关于慢时间u的连续函数就是所求的修正距离历史；观测场景Ω中的其它点目标的修正距离历史向量和修正距离历史的获取方法同点目标P；

步骤13、求解补偿的相位：

利用修正距离历史可以求出本发明所提出的新的BP算法在成像过程中需要补偿的相位 ${\widehat{\mathrm{\φ}}}_m(u;P):{\widehat{\mathrm{\φ}}}_m(u;P)=2\mathrm{\π}{\hat{R}}_m(u;P)/\mathrm{\λ},$ 其中λ为雷达载波波长；

步骤14、利用新的BP算法获取SAR图像：

求出需要补偿的相位以后，通过采用传统的脉冲相关积累方法对距离脉冲压缩后的回波数据进行相关积累，获得点目标P的SAR图像即：

$\hat{I}\left(P\right)=\underset{-T/2}{\overset{T/2}{\∫}}{S}_c(\frac{{\hat{R}}_m(u;P)}{c},u;P)\·\exp \left(j{\widehat{\mathrm{\φ}}}_m(u;P)\right)\mathrm{du}$

其中，T表示合成孔径时间；u表示方位向慢时间；c表示电磁波的传播速度；表示修正距离历史；S_c(t,u；P)表示由步骤4得到的经过距离脉冲压缩后的回波数据；表示利用新的BP算法所获得的点目标P的SAR图像；利用求取点目标P的SAR图像的上述公式，获得观测场景Ω中的其它点目标的SAR图像。