CN101900805B - 一种球面波成像数学模型及近场效应的补偿方法 - Google Patents

Abstract

本发明涉及一种球面波成像数学模型及近场效应的补偿方法，属于微波成像技术领域。本发明的目的是为了解决微波成像在不满足远场条件的情况下，现有基于转台目标成像的经典数学模型难以提供目标成像的精度，测量结果近场效应明显，距离越近，位置偏移量越大，幅度估计误差越大等问题，提出了一种球面波成像的数学模型，并根据此模型对成像的近场效应进行修正补偿的方法。其原理为在原有数学模型的基础上，通过改进测量距离Rθ的精度，构造新的球面波成像的数学模型，并用此模型对成像的近场效应采用二维ESPRIT超分辨算法来进行修正补偿。本方法能够有效地补偿微波成像的近场效应，且具有运算量小，误差修正效果明显等特点。

Description

一种球面波成像数学模型及近场效应的补偿方法

技术领域

本发明涉及一种球面波成像数学模型及近场效应的补偿方法，属于微波成像技术领域。

背景技术

微波是频率在300MHz～300GHz，相应波长为1m～1mm的电磁波。微波具有频率高、频带宽、信息容量大、波长短、能穿透电离层和方向性好等特点，微波成像是指以微波作为信息载体的一种成像手段，其原理是用微波照射被测物体，然后通过物体外部散射场的测量值来重构物体的形状或(复)介电常数分布。微波成像是依赖电磁波与目标的相互作用，从散射回波信号中挖掘、提取目标信息，重构目标特征。

按原理来说，微波成像通常应在远场条件下实施，照射目标的电磁波为平面波；但是，在实际情况中，由于发射功率的限制或者场地的限制，微波成像或微波诊断测量只能在一定的距离条件下进行，难以满足经典远场条件，照射目标的电磁波为球面波，在目标处产生波前弯曲，从而对目标回波成像时发生位置偏移和幅度估计不准确的现象，产生近场效应，造成误差。

为了补偿微波成像的近场效应，科研工作者做了很多相关研究。文献“NearField Imaging For Conducting Objects”，Li Hsueh-Jyh，Lin Feng-Li，IEEE TransAntenna and Propagation，1991，39(5)：600-605.通过补偿球面波的波前相位实现了目标的近场二维成像；文献：“引入双程传播损耗的近场大型目标球面逆投影成像算法研究”，孙厚军，李世勇，胡伟东，吕昕，电子学报，2009，37(3)：449-453.又考虑了近场大尺寸目标情况下对球面波双程传播损耗的补偿问题。

但这些算法都没有给出不同距离条件下的误差及近场效应并加以修正，包括修正位置偏移和进行精确的幅度估计。

超分辨类算法主要是基于空间谱估计的算法，空间谱估计技术是近三十年来发展起来的一门新兴的空域信号处理技术，其主要目标是研究提高在信号带宽内空间信号(包括独立、部分相关和相干)的估计精度、角度分辨力和提高运算速度的各种算法。超分辨算法最典型的有MUSIC、ESPRIT及WSF等算法。

多信号分类(MUSIC)成像算法具有运算量大、分辨门限高，且散射点数目难以确定等局限性；旋转不变参数估计技术(Estimation of Signal Parameter viaRotational Invariance Techniaques，缩写ESPRIT)已被应用于雷达成像中，但是研究者往往只是考虑远场条件下的微波成像，如文献“基于一种二维ESPRIT算法的ISAR超分辨成像”，焦芸，田野，宿富林，哈尔滨师范大学自然科学学报，第21卷，第04期，2005年，没有考虑近场成像中的应用问题。

基于矩阵束的ESPRIT算法运算量比较小，而且能够实现二维坐标及散射强度的自动配对，无须额外的人为配对，因此可以将该方法直接应用于近场二维成像的修正，补偿近场效应。

为了更好地说明本发明的内容，下面对近场转台目标成像的经典数学模型和超分辨ESPRIT算法进行简要介绍。

图1为近场转台目标二维成像的示意图，目标置于转台之上，可以随着转台绕中心旋转，雷达固定不动，采用单天线系统发射和接收电磁波。

假设g(x，y)为目标二维散射密度函数，目标与x-y坐标系绕O点顺时钟旋转，θ为u轴相对于x轴的旋转角度；x-y和u-v坐标系的原点均为O，R₀为雷达到目标旋转中心O的距离，R_θ为目标在相对于u轴转过角度θ时，目标上任一点到达雷达天线的距离；φ为R_θ与R₀之间的夹角；

为目标x-y各点相对于坐标系的极坐标；弧线S代表到达雷达天线等距离的散射点的连线。

x-y和u-v坐标系的关系为：

$\left\{\begin{array}{c}u=x\cos \mathrm{\θ}+y\sin \mathrm{\θ}\\ v=-x\sin \mathrm{\θ}+y\cos \mathrm{\θ}\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{c}x=u\cos \mathrm{\θ}+v\sin \mathrm{\θ}\\ y=u\sin \mathrm{\θ}+v\cos \mathrm{\θ}\end{array}\right.---\left(2\right)$

当目标位于远场时，雷达回波信号经混频及坐标变换以后可以表示为：

$z(k,\mathrm{\θ})=\underset{D}{\text{\∫\∫}}g(x,y)\exp [-j2\mathrm{\πk}(x\cos \mathrm{\θ}+y\sin \mathrm{\θ})]$

$\·\mathrm{dxdy}\·\exp (-j2\mathrm{\πk}{R}_0)/R_0^2---\left(3\right)$

其中k＝2f/c，f为电磁波的频率，c为光速。

但是在近场情况下时，由于照射目标的电磁波近似为球面波，一般难以满足远场测量条件，球面波前弯曲难以消除。根据合成孔径或者逆合成孔径的原理，接收到的目标回波可以被合成出二维像，那么(3)式就应该表示为：

$z(k,\mathrm{\θ})=\underset{D}{\text{\∫\∫}}g(x,y)\frac{\exp (-j2\mathrm{\π}{R}_{\mathrm{\θ}})}{R_{\mathrm{\θ}}^2}\mathrm{dxdy}---\left(4\right)$

上式为目前通用的近场转台目标成像的经典数学模型。但是该模型中R_θ的近似值不精确，造成测量结果的近场效应尤为明显，目标多散射中心的位置和强度和实际值相比都会发生偏差，而且测量距离越近，这种偏差越大。

超分辨ESPRIT算法最基本的假设就是存在两个完全相同的子阵，且两个子阵的间距是已知的。由于两个子阵的结构完全相同，且子阵的阵元数为m，对以同一信号而言，两个子阵的输出只有一个相位差φ_i，i＝1，2，...，n。假设子阵1的接收数据为X₁，子阵2的的接收数据为X₂，则可将两个子阵表示为：

X₁＝[a(θ₁)…a(θ_n)]S+N₁＝AS+N₁

$X_2=[a\left({\mathrm{\θ}}_1\right)e^{j{\mathrm{\φ}}_1}\·\·\·a\left({\mathrm{\θ}}_n\right)e^{j{\mathrm{\φ}}_n}]S+N_2=\mathrm{A\ΦS}+N_2$

其中，N₁、N₂分别为子阵1和子阵2的噪声数据。上式中，子阵1的阵列流型A₁＝A，子阵2的阵列流型A₂＝AΦ，且式中

$\mathrm{\Φ}=\mathrm{diag}[e^{j{\mathrm{\φ}}_1}\·\·\·e^{j{\mathrm{\φ}}_2}]$

由于需要求解的是信号的方向，而信号的方向信息包含在A和Φ中，只要得到两个子阵间的旋转不变关系Φ，就可以方便地得到关于信号的方向信息。

先将两个子阵的模型进行合并，即：

$X=\begin{array}{c}\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right]\end{array}=\left[\begin{array}{c}A\\ \mathrm{A\Φ}\end{array}\right]S+N=\stackrel{\‾}{A}S+N$

在理想条件下，可得上式的协方差矩阵：

$R=E\left[{\mathrm{XX}}^H\right]=\stackrel{\‾}{A}{R}_S{\stackrel{\‾}{A}}^H+R_N$

对上式进行特征值分解可得：

$R=\underset{i=1}{\overset{2m}{\∑}}{\mathrm{\λ}}_i{e}_i{e}_i^H=U_S{\∑}_S{U}_S^H+U_N{\∑}_N{U}_N^H$

上式中得到的特征值有如下关系：λ₁≥…≥λ_N＞λ_N+1＝…＝λ_2m，U_S为大特征值对应的特征向量张成的信号子空间，U_N为小特征值对应矢量张成的噪声子空间。对于实际的快拍数据，上式应修正如下：

$\stackrel{\mathrm{\Λ}}{R}={\stackrel{\mathrm{\Λ}}{U}}_S{\stackrel{\mathrm{\Λ}}{\∑}}_S{\stackrel{\mathrm{\Λ}}{U}}_S^H+{\stackrel{\mathrm{\Λ}}{U}}_N{\stackrel{\mathrm{\Λ}}{\∑}}_N{\stackrel{\mathrm{\Λ}}{U}}_N^H$

上述特征分解中大特征矢量张成的信号子空间与阵列流型张成的信号子空间是相等的，即：

$\mathrm{span}\left\{U_S\right\}=\mathrm{span}\left\{\stackrel{\‾}{A}\left(\mathrm{\θ}\right)\right\}$

此时，存在一个惟一的非奇异矩阵T，使得

$U_S=\stackrel{\‾}{A}\left(\mathrm{\θ}\right)T$

显然上述的结构对两个子阵都成立，所以有：

$U_S=\left[\begin{array}{c}{U}_{S1}\\ U_{S2}\end{array}\right]=\left[\begin{array}{c}\mathrm{AT}\\ \mathrm{A\ΦT}\end{array}\right]$

由两个子阵在阵列流型上的关系可知：

A₂＝A₁Φ

于是有：U_S2＝U_S1T^-1ΦT＝U_S1Ψ

上式反映了两个子阵的阵列接收数据的信号子空间的旋转不变性。

如果阵列流型A是满秩矩阵，则可以得到：

Φ＝TΨT^-1

所以上式中Ψ的特征值组成的对角阵一定等于Φ，而矩阵T的各列就是Ψ的特征矢量。所以一旦得到上述的旋转不变关系矩阵Ψ，就可以直接得到信号的方向信息。

该方法无需经过多次矩阵变换就可以达到参数自动配对的目的，在保证二维ESPRIT算法参数估计性能的前提下，简化了参数的配对过程，经实践验证了该方法的有效性。

发明内容

本发明的目的是为了解决微波成像在不满足远场条件的情况下，现有基于转台目标成像的经典数学模型难以提供目标成像的精度，测量结果近场效应明显，距离越近，位置偏移量越大，幅度估计误差越大等问题，提出了一种球面波成像的数学模型，并根据此模型对成像的近场效应进行修正补偿的方法。

该方法的原理为在原有数学模型的基础上，通过改进测量距离R_θ的精度，构造新的球面波成像的数学模型，并用此模型对成像的近场效应采用二维ESPRIT超分辨算法来进行修正补偿。

一种球面波成像数学模型及近场效应的补偿方法，其具体实现步骤如下：

步骤1，对现有转台目标二维成像数学模型

作改进如下：令

其中，g(x，y)为目标二维散射密度函数，目标与x-y坐标系绕O点顺时钟旋转，θ为u轴相对于x轴的旋转角度；x-y和u-v坐标系的原点均为O，R₀为雷达到目标旋转中心O的距离，R_θ为目标在相对于u轴转过角度θ时，目标上任一点到达雷达天线的距离；φ为R_θ与R₀之间的夹角；

为目标x-y各点相对于坐标系的极坐标；弧线S代表到达雷达天线等距离的散射点的连线。

将R_θ做二阶泰勒级数做近似可得：

步骤2，由于实际情况中一般满足小角度成像的条件，θ非常小，于是步骤1得到的R_θ可以化简为：

$R_{\mathrm{\θ}}=R_0+x+\mathrm{y\θ}+\frac{y^2-2\mathrm{xy\θ}}{2(R_0+x)}$

步骤3，在步骤2的基础上，对

项进一步做近似：

$R_{\mathrm{\θ}}^2\≈{(R_0+x+\frac{y^2}{{2R}_0+2x})}^2$

这样做的近似能使幅度估计误差非常小。

步骤4，在步骤3的基础上，对于频率步进体制测量雷达，满足k＝k₀+nΔk(Δk远小于k₀)，(其中k＝2f/c，f为电磁波的频率，c为光速，n＝0，1，…N-1)；θ＝mΔθ，m＝0，1，…M-1。设P为目标区域独立的散射点个数(P＜N，M)。对步骤1所述的现有转台目标二维成像数学模型进行离散化，并进一步忽略二阶小量，得：

$z(\mathrm{n\Δk},\mathrm{m\Δ\θ})=\underset{p=0}{\overset{P-1}{\∑}}g(x_p,y_p)\mathrm{\Δx\Δy}$

$/{(R_0+x_p+\frac{y_p^2}{{2R}_0+2x_p})}^2$

·exp[-j2π(k₀+nΔk)R₀]

$\·\exp [-j2\mathrm{\π}(k_0+\mathrm{n\Δk})\·(x_p+\frac{y_p^2}{{2R}_0+{2x}_p})]$

$\·\exp [-j2\mathrm{\π}{k}_0\mathrm{m\Δ\θ}(y_p-\frac{x_p{y}_p}{R_0+x_p})]$

其中，x_p，y_p为第P个目标散射点的横纵坐标。上式即为改进后的近场目标二维成像的数学模型。

针对这个模型可以采用MUSIC，ESPRIT，WSF等多种超分辨算法进行参数的估计，经分析基于矩阵束的ESPRIT算法运算量比较小，而且能够实现二维坐标及散射强度的自动配对，无须额外的人为配对，因此可以将该方法直接应用于近场二维成像的修正，补偿近场效应。

步骤5，进一步变形步骤4得到的数学模型。令：

$\left\{\begin{array}{c}{z}^{\text{\′}}(n,m)=z(\mathrm{n\Δk},\mathrm{m\Δ\θ})\exp \left[j2\mathrm{\π}(k_0+\mathrm{n\Δk})R_0\right]\\ x_p^{\′}=x_p+\frac{y_p^2}{{2R}_0+{2x}_p}\\ y_p^{\′}=y_p-\frac{x_p{y}_p}{R_0+x_p}\end{array}\right.$

由于ΔxΔy为常数量，在数学变换中将被合并为一系数，因此在模型中可以省略，则步骤4所述的改进后的近场目标二维成像的数学模型变为：

$z^{\′}(n,m)=\underset{p=0}{\overset{P-1}{\∑}}g(x_p,y_p)/{(R_0+x_p^{\′})}^2$

·exp(-j2πk₀x′_p)

·exp(-j2πnΔkx′_p)

·exp(-j2πk₀mΔθy′_p)

步骤6，在步骤5的基础上，进行子空间矩阵分裂。

由于z′(n，m)与z′(n+1，m)及z′(n，m+1)分别相差一个固定相位，则将z′(n，m)按附图3所示方式排列，并形成三个子空间矩阵X，Y，Z。

根据步骤5所得z′(n，m)的表示式，可以将X，Y，Z表示成如下形式：

$\left\{\begin{array}{c}X=\mathrm{AS}+N_x\\ Y=\mathrm{A\ΦS}+N_y\\ Z=\mathrm{A\ΘS}+N_z\end{array}\right.$

X，Y，Z均为(N-1)×(M-1)维矩阵，A为(N-1)×(M-1)×P维矩阵，N_i(i＝x，y，z)为测量噪声矩阵。而且有：

$\left\{\begin{array}{c}A(n,m,p)=\exp [-j2\mathrm{\π}(k_0+\mathrm{n\Δk})x_p^{\′}]\\ \·\exp (-j2\mathrm{\π}{k}_0\mathrm{m\Δ\θ}{y}_p^{\′})\\ S=[g(x_0,y_p)/{(R_0+x_p^{\′})}^2,g(x_1,y_1)/{(R_0+x_1^{\′})}^2,\\ \·\·\·,g(x_P,y_P)/{(R_0+x_P^{\′})}^2]\\ \mathrm{\Φ}=\mathrm{diag}(\exp (-j2\mathrm{\π\Δk}{x}_0^{\′}),\exp (-j2\mathrm{\π\Δk}{x}_1^{\′}),\\ \·\·\·,\exp (-j2\mathrm{\π\Δ}{\mathrm{kx}}_P^{\′})\\ \mathrm{\Θ}=\mathrm{diag}(\exp (-j2\mathrm{\π}{k}_0\mathrm{\Δ\θ}{y}_0^{\′}),(\exp (-j2\mathrm{\π}{k}_0\mathrm{\Δ\θ}{y}_1^{\′}),\\ \·\·\·,(\exp (-j2\mathrm{\π}{k}_0\mathrm{\Δ\θ}{y}_P^{\′})]\end{array}\right.$

步骤7，在步骤6的基础上进行奇异值分解，得到目标散射点的个数。

将X，Y，Z构成如下两个矩阵：

$\left[\begin{array}{ccc}X& Y& Z\end{array}\right]=U_1{D}_1{V}_1^H$ $\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=U_2{D}_2{V}_2^H$

其中，U₁，V₁为奇异值矩阵。根据上式估计出目标散射点个数P。

步骤8，通过步骤7中奇异值矩阵U₁，V₁计算的特征值，得到步骤6中Φ和Θ的对角元素的值。

取步骤7得到的奇异值矩阵U₁，V₁中与P个散射点对应的奇异值的列矢量构成

从而可以得到三个新矩阵：

$\left\{\begin{array}{c}{E}_x={\stackrel{\mathrm{\Λ}}{U}}_1^{H}X{\stackrel{\mathrm{\Λ}}{V}}_2\\ E_y={\stackrel{\mathrm{\Λ}}{U}}_1^{H}Y{\stackrel{\mathrm{\Λ}}{V}}_2\\ E_z={\stackrel{\mathrm{\Λ}}{U}}_1^{H}Z{\stackrel{\mathrm{\Λ}}{V}}_2\end{array}\right.$

并定义

$E_{\mathrm{\θ}}=E_x^{-1}{E}_z$

计算

、E_θ的特征值就可以得到Φ和Θ的对角元素的值。

步骤9，在步骤8的基础上，确定变换矩阵Q使得对角化，并将这种变换应用于E_θ，即：

Q^HE_θQ＝T_θ

这里T_θ是一个近似上三角矩阵，和T_θ的主对角元素分别等于Φ和Θ的主对角元素。由于采用了相同的变换矩阵Q，

和T_θ的主对角元素是一一对应的，故形成了参数的自动配对。

实际上T_θ已经非常接近上三角矩阵，再做进一步的变换反而可能造成参数失配情况，这里直接取T_θ的对角元素为特征值。

步骤10，在步骤9的基础上，得到目标散射点的位置估计坐标x_p和y_p。

设

θ_p分别为

和T_θ对应的主对角元素的相位，那么可以得到如下关系：

当求解完x′_p和y′_p后，再解一个二次方程组即可以得到对应的x_p和y_p。

步骤11，完成步骤10后，计算得到目标散射点的散射强度。

在小角度成像条件下，可以认为任意一个采样角度散射点的散射强度是不变的，那么令m＝0，A可以降为一个(N-1)×P维矩阵，然后可以得到S的最小二乘估计值：

S＝(A^HA)^-1A^HX

由步骤6中第二个方程组便可求得散射强度g(x_p，y_p)。

有益效果

本发明针对近场情况下转台目标二维成像问题，给出了一种新的近场二维成像数学模型，该模型引入了目标散射点的横纵坐标信息，在此基础上，利用超分辨ESPRIT算法，直接代入目标散射点的横纵坐标信息完成近场参数估计。该方法给出了近场目标的散射中心的个数、位置及强度的精确估计，能够有效地补偿微波成像的近场效应，且具有运算量小，误差修正效果明显等特点。

本发明的方法可以应用到近场微波成像、SAR成像、ISAR成像等诸多工程领域中。

附图说明

图1为现有技术的近场转台目标二维成像图；

图2为本发明一种球面波成像数学模型及近场效应补偿方法的流程图；

图3为本发明近场效应补偿方法中步骤2所述的z′(n，m)排列方式图；其中，“…”点画线所示范围表示子空间矩阵X，“-”长画线所示范围表示子空间矩阵Y，“---”短画线所示范围表示子空间矩阵Z。

图4为具体实施方式中R₀为33600m时未经近场修正情况下二维坐标估计值与标准值图；

图5为具体实施方式中R₀为33600m时经过近场修正情况下二维坐标估计值与标准值图；

图6为具体实施方式中R₀为100m时未经近场修正情况下二维坐标估计值与标准值图；

图7为具体实施方式中R₀为100m时经过近场修正情况下二维坐标估计值与标准值图；

图8为具体实施方式中R₀为50m时未经近场修正情况下二维坐标估计值与标准值图；

图9为具体实施方式中R₀为50m时经过近场修正情况下二维坐标估计值与标准值图；

图10为具体实施方式中近场修正前后散射强度对比图。

具体实施方式

为了更好的说明本发明的目的和优点，下面结合附图和实施例对本发明作进一步的详细说明：

本发明一种球面波成像数学模型及近场效应补偿方法的实现步骤如图2所示。其中，步骤2所述的z′(n，m)排列方式如图3所示。“…”点画线所示范围表示子空间矩阵X，“-”长画线所示范围表示子空间矩阵Y，“---”短画线所示范围表示子空间矩阵Z。

本实施例中，一个近场情况下由三个独立散射点组成的模型，这三个点按照横向坐标从小到大的顺序依次命名为点A、B、C。对应的坐标为(-6-6)、(00)、(66)，散射强度都为0dBsm。起始频率f₀设为35GHz，跳频步进Δf设为2MHz，采样角度间隔Δθ设为0.004°。取M和N都为128。由于目标直径为12米，那么对应远场条件为R₀≥33600m。为了对比，分别取R₀为33600m，100m，50m，通过在这三种距离条件下近场修正前后的参数估计来说明该方法的有效性。

R₀为33600m距离条件下：

1.将该距离条件直接代入步骤5中的变形后的数学模型，并根据步骤6的方法形成128×128维数据矩阵z′(n，m)；进行子空间矩阵分裂，将z′(n，m)按附图3所示方式排列，并形成三个127×127维子空间矩阵X，Y，Z，表示成如下形式：

$\left\{\begin{array}{c}X=\mathrm{AS}+N_x\\ Y=\mathrm{A\ΦS}+N_y\\ Z=\mathrm{A\ΘS}+N_z\end{array}\right.$

A为127×127×3维矩阵，N_i(i＝x，y，z)为测量噪声矩阵。

2.根据步骤7进行奇异值分解，得到特征目标散射点的个数为3个。

3.根据步骤7得到的奇异值矩阵U₁，V₁中与3个散射点对应的奇异值的列矢量构成

从而可以得到三个新矩阵：

$\left\{\begin{array}{c}{E}_x={\stackrel{\mathrm{\Λ}}{U}}_1^{H}X{\stackrel{\mathrm{\Λ}}{V}}_2\\ E_y={\stackrel{\mathrm{\Λ}}{U}}_1^{H}Y{\stackrel{\mathrm{\Λ}}{V}}_2\\ E_z={\stackrel{\mathrm{\Λ}}{U}}_1^{H}Z{\stackrel{\mathrm{\Λ}}{V}}_2\end{array}\right.$

定义 $E_{\mathrm{\θ}}=E_x^{-1}{E}_z$

计算

E_θ的特征值就可以得到Φ和Θ的对角元素的值；

4.根据步骤9所述方法确定变换矩阵Q使得

对角化，并将这种变换应用于E_θ，由于采用了相同的变换矩阵Q，

和T_θ的主对角元素是一一对应的，形成了参数的自动配对；

5.解步骤10所述二次方程组，得到对应的x_p和y_p，即目标的位置估计；

6.根据步骤11所述方法，由步骤6中第二个方程组便可求得散射强度g(x_p，y_n)。

图4为R₀33600m距离条件下未经近场修正情况下二维坐标估计值与标准值图，图5为经近场修正情况下二维坐标估计值与标准值图，其中″+″代表估计值，″o″代表标准值。

R₀为100m距离条件下，按照R₀为33600m距离条件时所述方法一样的步骤进行处理，得到目标的位置估计以及散射强度g(x_p，y_p)。

图6为R₀为100m距离条件下未经近场修正情况下二维坐标估计值与标准值图，图7为经近场修正情况下二维坐标估计值与标准值图，其中″+″代表估计值，″o″代表标准值。

R₀为50m距离条件下，按照R₀为33600m距离条件时所述方法一样的步骤进行处理，得到目标的位置估计以及散射强度g(x_p，y_p)。

图8为R₀为50m距离条件下未经近场修正情况下二维坐标估计值与标准值图，图9为经近场修正情况下二维坐标估计值与标准值图，其中″+″代表估计值，″o″代表标准值。

图10为在A、B、C三个点近场修正前后散射强度对比图。

通过上述三种距离条件下近场修正前后的结果对比看出：当R₀为33600m时，目标位于远场，这时近场修正前后参数估计值相差不大。当R₀为100m时，目标处于近场，经过近场修正后的估计值很接近标准值，明显好于未修正的情况。当距离更近，R₀取为50m时，图9及图10所示的修正效果更为明显。由此可见，本文提出的近场修正方法取得了很好的效果。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种球面波成像数学模型及近场效应的补偿方法，其特征在于：包括如下步骤：

步骤1，对现有转台目标二维成像数学模型 $z(k,\mathrm{\θ})=\underset{D}{\∫\∫}g(x,y)\frac{\exp (-j2\mathrm{\πk}{R}_{\mathrm{\θ}})}{R_{\mathrm{\θ}}^2}\mathrm{dxdy}$ 作改进如下：令

$R_{\mathrm{\θ}}=\sqrt{{(R_0+u)}^2+v^2}$

其中，g(x，y)为目标二维散射密度函数，目标与x-y坐标系绕O点顺时钟旋转，θ为u轴相对于x轴的旋转角度；x-y和u-v坐标系的原点均为O，R₀为雷达到目标旋转中心O的距离，R_θ为目标在相对于u轴转过角度θ时，目标上任一点到达雷达天线的距离；φ为R_θ与R₀之间的夹角；

为目标各点相对于x-y坐标系的极坐标；弧线S代表到达雷达天线等距离的散射点的连线；

将R_θ做二阶泰勒级数近似可得：

$=R_0+x\cos \mathrm{\θ}+y\sin \mathrm{\θ}+\frac{{(y\cos \mathrm{\θ}-x\sin \mathrm{\θ})}^2}{2(R_0+x\cos \mathrm{\θ}+y\sin \mathrm{\θ})}$

步骤2，由于实际情况一般满足小角度成像的条件，θ非常小，于是步骤1得到的R_θ可以化简为：

$R_{\mathrm{\θ}}=R_0+x+\mathrm{y\θ}+\frac{y^2-2\mathrm{xy\θ}}{2(R_0+x)}$

步骤3，在步骤2的基础上，对

项进一步做近似：

$R_{\mathrm{\θ}}^2\≈{(R_0+x+\frac{y^2}{2R_0+2x})}^2$

步骤4，在步骤3的基础上，对于频率步进体制测量雷达，满足k＝k₀+nΔk，其中k＝2f/c，f为电磁波的频率，c为光速，n＝0，1，…N-1，Δk远小于k₀；θ＝mΔθ，其中Δθ为采样角度间隔，m＝0，1，…M-1；设P为目标区域独立的散射点个数，且P＜N，P＜M；对步骤1所述的现有转台目标二维成像数学模型进行离散化，并进一步忽略二阶小量，得：

$z(\mathrm{n\Δk},\mathrm{m\Δ\θ})=\underset{p=0}{\overset{P-1}{\mathrm{\Σ}}}g(x_p,y_p)\mathrm{\Δx\Δy}$

${/(R_0+x_p+\frac{y_p^2}{2R_0+2x_p})}^2$

$\·\exp [-j2\mathrm{\π}(k_0+\mathrm{n\Δk})R_0]$

$\·\exp [-j2\mathrm{\π}(k_0+\mathrm{n\Δk})\·(x_p+\frac{y_p^2}{2R_0+2x_p})]$

$\·\exp [-j2\mathrm{\π}{k}_0\mathrm{m\Δ\θ}(y_p-\frac{x_p{y}_p}{R_0+x_p})]$

其中，x_p，y_p为第P个目标散射点的横纵坐标；上式即为改进后的近场目标二维成像的数学模型；

步骤5，进一步变形步骤4得到的数学模型；令：

$\left\{\begin{array}{c}{z}^{\′}(n,m)=z(\mathrm{n\Δk},\mathrm{m\Δ\θ})\exp \left[j2\mathrm{\π}(k_0+\mathrm{n\Δk})R_0\right]\\ x_p^{\′}=x_p+\frac{y_p^2}{2R_0+2x_p}\\ y_p^{\′}=y_p-\frac{x_p{y}_p}{R_0+x_p}\end{array}\right.$

由于ΔxΔy为常数量，在数学变换中将被合并为一系数，因此在模型中可以省略，则步骤4所述的改进后的近场目标二维成像的数学模型变为：

$z^{\′}(n,m)=\underset{p=0}{\overset{P-1}{\mathrm{\Σ}}}g(x_p,y_p)/{(R_0+x_p^{\′})}^2$

$\·\exp (-j2\mathrm{\π}{k}_0{x}_p^{\′})$

$\·\exp (-j2\mathrm{\πn\Δk}{x}_p^{\′})$

$\·\exp (-j2\mathrm{\π}{k}_0\mathrm{m\Δ\θ}{y}_p^{\′})$

步骤6，在步骤5的基础上，进行子空间矩阵分裂；

由于z′(n，m)与z′(n+1，m)及z′(n，m+1)分别相差一个固定相位，则将z′(n，m)按一定方式排列，并形成三个子空间矩阵X，Y，Z；

根据步骤5所得z′(n，m)的表示式，可以将X，Y，Z表示成如下形式：

$\left\{\begin{array}{c}X=\mathrm{AS}+N_x\\ Y=\mathrm{A\ΦS}+N_y\\ Z=\mathrm{A\ΘS}+N_z\end{array}\right.$

X，Y，Z均为(N-1)×(M-1)维矩阵，A为(N-1)×(M-1)×P维矩阵，N_i(i＝x，y，z)为测量噪声矩阵；而且有：

$\left\{\begin{array}{c}A(n,m,p)=\exp [-j2\mathrm{\π}(k_0+\mathrm{n\Δk})x_p^{\′}]\\ \·\exp (-j2\mathrm{\π}{k}_0\mathrm{m\Δ\θ}{y}_p^{\′})\\ S=[g(x_0,y_p)/{(R_0+x_p^{\′})}^2,g(x_1,y_1)/{(R_0+x_1^{\′})}^2,\\ ...,g(x_P,y_P)/{(R_0+x_P^{\′})}^2]\\ \mathrm{\Φ}=\mathrm{diag}(\exp (-j2\mathrm{\π\Δk}{x}_0^{\′}),\exp (-j2\mathrm{\π\Δk}{x}_1^{\′})),\\ ...,\exp (-j2\mathrm{\π\Δk}{x}_P^{\′})]\end{array}\right.$

步骤7，在步骤6的基础上进行奇异值分解，得到目标散射点的个数；

得到目标散射点个数的具体方法为：将步骤6中得到的X，Y，Z构成如下两个矩阵：

[X Y Z]＝U₁D₁V₁ ^H $\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=U_2{D}_2{V}_2^H$

其中，U₁，V₁为奇异值矩阵；根据上式估计出目标散射点个数P；

步骤8，通过步骤7中奇异值矩阵U₁，V₁计算的特征值，得到步骤6中Φ和Θ的对角元素的值；

取步骤7得到的奇异值矩阵U₁，V₁中与P个散射点对应的奇异值的列矢量构成从而可以得到三个新矩阵：

$\left\{\begin{array}{c}{E}_x={\stackrel{\mathrm{\Λ}}{U}}_1^{H}X{\stackrel{\mathrm{\Λ}}{V}}_2\\ E_y={\stackrel{\mathrm{\Λ}}{U}}_1^{H}Y{\stackrel{\mathrm{\Λ}}{V}}_2\\ E_z={\stackrel{\mathrm{\Λ}}{U}}_1^{H}Z{\stackrel{\mathrm{\Λ}}{V}}_2\end{array}\right.$

并定义

计算

Eθ的特征值就可以得到Φ和Θ的对角元素的值；

步骤9，在步骤8的基础上，确定变换矩阵Q使得

对角化，并将这种变换应用于E_θ，即：

Q^HE_θQ＝T_θ

这里T_θ是一个近似上三角矩阵，和T_θ的主对角元素分别等于Φ和Θ的主对角元素；由于采用了相同的变换矩阵Q，

和T_θ的主对角元素是一一对应的，故形成了参数的自动配对；

实际上T_θ已经非常接近上三角矩阵，再做进一步的变换反而可能造成参数失配情况，这里直接取T_θ的对角元素为特征值；

步骤10，在步骤9的基础上，得到目标散射点的位置估计坐标x_p和y_p；

步骤11，完成步骤10后，计算得到目标散射点的散射强度。

2.根据权利要求1所述的一种球面波成像数学模型及近场效应的补偿方法，其特征在于：所述的步骤6中z (n，m)按一定方式排列，形成三个子空间矩阵X，Y，Z的排列方式为：

3.根据权利要求1所述的一种球面波成像数学模型及近场效应的补偿方法，其特征在于：所述的步骤10中得到目标散射点的位置估计坐标的方法为：设

θ_p分别为

和T_θ对应的主对角元素的相位，那么可以得到如下关系：

当求解完x′_p和y′_p后，再解一个二次方程组即可以得到对应的x_p和y_p。

4.根据权利要求1所述的一种球面波成像数学模型及近场效应的补偿方法，其特征在于：所述的步骤11中得到目标散射强度的方法为：

在小角度成像条件下，可以认为任意一个采样角度散射点的散射强度是不变的，那么令m＝0，A可以降为一个(N-1)×P维矩阵，然后可以得到S的最小二乘估计值：

S＝(A^HA)^-1A^HX

由所述步骤6中第二个方程组便可求得散射强度g(x_p，y_p)。

Images