CN111999712B - 一种三维gtd散射中心模型的散射中心点的参数估计方法 - Google Patents

Abstract

本发明公开了一种三维GTD散射中心模型的散射中心点的参数估计方法，包括以下步骤，S1：基于雷达目标后向电磁数据构建Hankel矩阵，得到原目标回波数据矩阵X^x；S2：定义大小为PQL×PQL的置换矩阵J，得到目标回波数据共轭数据信息的矩阵E_conj；S3：将步骤S2中得到的矩阵E_conj及原目标回波数据矩阵X^x的自相关矩阵进行相加，并作取平均处理，得到一个过程协方差矩阵R，S4：令最终总协方差矩阵R₁＝RR^H＝R²；S5：利用经典3D‑ESPRIT算法对矩阵R₁进行处理，即可得到三维GTD散射中心模型中的各个参数。本发明中的方法可提高对目标原始回波数据的利用率并增大信号特征值与噪声特征值之间差距，等效为增大信噪比，从而提高了散射中心的分辨率。

Description

一种三维GTD散射中心模型的散射中心点的参数估计方法

技术领域

本发明涉及雷达成像技术领域，尤其涉及一种三维GTD散射中心模型的散射中心点的参数估计方法。

背景技术

随着现代雷达技术的不断发展以及对于全空间、全极化目标电磁散射特性数据需求的不断增长，存储、利用这些海量目标电磁散射特性数据的压力也日益加大。在高频区，雷达目标的后向电磁散射回波可视为有限个强散射点的相干叠加，而这些强散射点通常可称为雷达目标的散射中心。散射中心作为雷达目标的电磁散射特性之一，在雷达目标识别，RCS频率内插与外推，雷达目标三维重构等军事领域应用十分广泛。

通过构建合理的散射中心模型，可更有效的描述目标的散射特性，对深入研究其散射机理起到重要的作用。目前，基于几何绕射理论的(GTD)散射中心模型对目标在高频区散射特性的刻画最为准确。GTD模型可细分为一维、二维、三维GTD模型，随着维数的增大，其对目标电磁散射特性刻画的越来越精准，但运算复杂度也相应的变大，参数估计的难度也随之增大。因此，大多数研究者主要利用各类算法，如MUSIC算法、ESPRIT算法、MEMP算法等对一维GTD、二维GTD模型进行参数估计提取。

此外，文献《一种基于3D-ESPRIT的散射中心参数估计算法[J]》([14].温晓杨,石志广,赵宏钟,等；雷达科学与技术,2007,5(2):119-123.)中利用3D-ESPRIT算法对散射中心模型参数进行估计，但存在参数失配与低信噪比时参数估计精度不高的问题；文献《光学区雷达目标散射中心提取及其应用研究[D]》(王菁；南京：南京航空航天大学，2010.)提出一种改进的3D-ESPRIT算法，有效解决了参数失配问题，但在信噪比情况下，该文献中改进算法的参数估计也较低。

根据文献《光学区雷达目标散射中心提取及其应用研究[D]》中提出的3D-GTD散射中心模型，目标后向电磁散射数据可等效为I个强散射中心叠加合成的，其数学表述为

式中：

表示目标后向电磁散射回波，I代表散射中心个数，{A_i，α_i，x_iy_i，z_i}分别表示第i个散射中心的散射强度、散射类型、横向距离、纵向距离及垂直距离。f_m＝f₀+mΔf，m＝1，2，…，M，f₀为起始频率，Δf为步进频率，m代表频率下标，M代表总频率个数；

θ_n＝θ₀+nΔθ，n＝1，2，...，N，其中θ₀为起始方位角，Δθ为步进方位角，n为方位角下标，N代表总方位角个数；

其中/>

为起始俯仰角，/>

为步进俯仰角，k为俯仰角下标，K为总俯仰角个数；nΔθ、/>

分别为方位方向上的小转角、俯仰方向上的小转角。

c＝3×10⁸m/s为电磁波传播速度，

为复高斯白噪声。α_i为0.5的整数倍，根据不同散射体可分为5种，二面角、三面角、平面法向反射，其α_i取值为1；单曲面反射、圆柱面反射，其α_i取值为0.5；双曲面反射、球面反射，其α_i取值为0；边缘绕射，其α_i取值为-0.5；尖顶绕射，其α_i取值为-1；

由于选择的雷达工作频率满足Δf/f₀＝1，因此可作如下近似

将式(2)中的近似结果代入式(1)中，将得到的公式变换到笛卡尔坐标并进行插值规范化处理、利用重采样技术，则目标的电磁回波数据可用下式(3)表示。

式中，m＝0，...，M-1，n＝0，...，N-1，k＝0，...，K-1；

f_x0、f_y0、f_z0分别为x、y、z方向的雷达信号初始频率；

P_yi＝exp(-4πjΔf_yy_i/c) (6)

P_zi＝exp(-4πjΔf_zz_i/c) (7)

其中B_i、P_xi、P_yi、P_zi均为过程参数，Δf_x、Δf_y、Δf_z分别代表雷达坐标系下x、y、z方向的步进频率，表达式如式(8)所示：

设f_c为雷达的中心频率，B为雷达工作带宽，则：

式(5)-(7)包含散射中心的类型参数及三类位置参数，因此可由下式(12)-(15)求解：

α_i＝(|P_xi|-1)f₀/Δf (12)

发明内容

针对上述存在的问题，本发明旨在提供一种三维GTD散射中心模型的散射中心点的参数估计方法，通过构建原始回波数据的共轭矩阵、对协方差矩阵叠加、取平均、平方处理，可提高对目标原始回波数据的利用率并增大信号特征值与噪声特征值之间差距，等效为增大了信噪比，从而提高了参数估计精度。

为了实现上述目的，本发明所采用的技术方案如下：

一种三维GTD散射中心模型的散射中心点的参数估计方法，其特征在于，包括以下步骤，

S1：基于雷达目标后向电磁数据构建Hankel矩阵，得到原目标回波数据矩阵X^x；

S2：定义大小为PQL×PQL的置换矩阵J，其反对角线位置上的元素为1，其他位置上的元素为0，即

则可得到目标回波数据共轭数据信息的矩阵E_conj，E_conj＝J·X^x；

S3：将步骤S2中得到的目标回波数据共轭数据信息矩阵E_conj及原目标回波数据矩阵X^x的自相关矩阵进行相加，并作取平均处理，得到一个过程协方差矩阵R

其中，R＝R^H；

S4：令最终总协方差矩阵R₁＝RR^H＝R²，则矩阵R₁、矩阵R两者特征值与特征向量具有以下关系式：

式中，λ₁、λ分别代表矩阵R₁与R的特征值，Λ₁、Λ分别代表矩阵R₁与R的特征向量；

S5：利用经典3D-ESPRIT算法对最终总协方差矩阵R₁进行处理，即可得到三维GTD散射中心模型中的散射中心点的各个参数。

进一步的，步骤S1中所述的原目标回波数据矩阵X^x大小为

PQL×(M-P+1)(N-Q+1)(K-L+1)，

且

式中，

其中，M/2≤P≤2M/3，N/2≤Q≤2N/3，K/2≤L≤2K/3；

表示目标后向电磁散射回波，f_m＝f₀+mΔf，m＝1，2，...，M，f₀为起始频率，Δf为步进频率，m代表频率下标，M代表总频率个数；θ_n＝θ₀+nΔθ，n＝1，2，...，N，其中θ₀为起始方位角，Δθ为步进方位角，n为方位角下标，N代表总方位角个数；

其中/>

为起始俯仰角，/>

为步进俯仰角，k为俯仰角下标，K为总俯仰角个数；nΔθ、/>

分别为方位方向上的小转角、俯仰方向上的小转角。

进一步的，步骤S5的具体操作包括，

S51：对最终总协方差矩阵R₁进行奇异值分解，则得到

其中，U_xS，V_xS代表R₁的信号子空间，分别由R₁的前I个主左特征向量与前I个主右特征向量构成；其中U_xN，V_xN代表R₁的噪声子空间，分别由R₁的非主左特征向量与非主右特征向量构成；

S52：构造矩阵F_x为

式中，U _xS、

为矩阵U_xS去除后Q×L行、去除前Q×L行得到的矩阵，/>

代表U _xS的广义逆；

S53：根据三维条件下置换矩阵之间的关系，得到y方向、z方向对应的信号子空间

U_ys＝E_xyU_xS

U_zs＝E_yzU_yS；

式中，E_xy、E_yz分别为三维条件下的置换矩阵；

进而得到y方向和z方向的矩阵F_y、F_z

式中，U _vS，

为矩阵U_vS去除后Q×L行、去除前Q×L行得到的矩阵；U _zS，/>

为矩阵U_zS去除后Q×L行、去除前Q×L行得到的矩阵；/>

分别代表U _vS和U _zS的广义逆；

S54：计算矩阵F_x、F_y、F_z前I个元素的主特征值向量Ψ_x、Ψ_y、Ψ_z

Ψ_x＝T_xF_xT_x ^-1

Ψ_y＝T_yF_yT_y ^-1

Ψ_z＝T_zF_zT_z ^-1；

S55：利用求得的矩阵Ψ_x、Ψ_y、Ψ_z求解其对应主对角线上的元素P_xi、P_yi和P_zi

P_xi＝diag(Ψ_x)，i＝1，K，I

P_yi＝diag(Ψ_y)，i＝1，K，I

P_zi＝diag(Ψ_z)，i＝1，K，I；

S56：根据求得的P_xi、P_yi和P_zi，求解类型参数α_i，横向距离参数x_i、纵向距离参数y_i与垂直距离参数z_i；

S57：根据类型参数α_i，横向距离参数x_i、纵向距离参数y_i与垂直距离参数z_i，利用最小二乘法，求得散射中心模型中的强度参数

式中，G＝[a₁，…，a_I]，

其中，G^H代表G的转置，(G^HG)^-1代表矩阵G^HG的共轭转置。

进一步的，步骤S53中三维条件下置换矩阵E_xy、E_yz分别为

式中，

表示Kronecker乘积，/>

代表(q，l)位置上元素为1，其他位置上元素为0的Q×L矩阵，/>

代表(l，p)位置上元素为1，其他位置上元素为0的LP矩阵，/>

代表(p，q)位置上元素为1，其他位置上元素为0的P×Q矩阵。

本发明的有益效果是：

本发明对经典的3D-ESPRIT算法进行了改进，利用改进后的3D-ESPRIT算法对三维GTD散射中心模型的参数进行估计，通过构建原始回波数据的共轭矩阵、对协方差矩阵叠加、取平均、平方处理，可提高对目标原始回波数据的利用率并增大信号特征值与噪声特征值之间差距，等效为增大了信噪比，从而提高了参数估计精度。将本发明中的参数研究方法应用到目标HRRP特征提取中，相比于传统逆傅立叶变换方法，其对目标散射中心的分辨率更高。

附图说明

图1为本发明仿真试验中x1-x4的均方差对比图；

图2为本发明仿真试验中y1-y4的均方差对比图；

图3为本发明仿真试验中z1-z4的均方差对比图；

图4为本发明仿真试验中α₁-α₄的均方差对比图；

图5为本发明仿真试验中A1-A4的均方差对比图；

图6为本发明仿真试验中两散射中心距离0.1m时的不同算法参数估计精度对比图；

图7为本发明仿真试验中两散射中心距离0.2m时的不同算法参数估计精度对比图。

具体实施方式

为了使本领域的普通技术人员能更好的理解本发明的技术方案，下面结合附图和实施例对本发明的技术方案做进一步的描述。

一种三维GTD散射中心模型的散射中心点的参数估计方法，包括以下步骤，

S1：基于雷达目标后向电磁数据构建Hankel矩阵，得到原目标回波数据矩阵X^x；

具体的，为更为准确地估计散射中心模型参数，首先要对目标的后向电磁散射数据进行空间平滑处理以去相关，而Hankel矩阵可替代空间平滑起到去相关的作用。因此，首先基于目标后向电磁数据构建Hankel矩阵。先沿x方向进行平滑，构建一个PQL×(M-P+1)(N-Q+1)(K-L+1)的原目标回波数据矩阵X^x，如下式(16)所示。

式中，

其中，M/2≤P≤2M/3，N/2≤Q≤2N/3，K/2≤L≤2K/3；

表示目标后向电磁散射回波，f_m＝f₀+mΔf，m＝1，2，…，M，f₀为起始频率，Δf为步进频率，m代表频率下标，M代表总频率个数；θ_n＝θ₀+nΔθ，n＝1，2，...，N，其中θ₀为起始方位角，Δθ为步进方位角，n为方位角下标，N代表总方位角个数；

其中/>

为起始俯仰角，/>

为步进俯仰角，k为俯仰角下标，K为总俯仰角个数；nΔθ、/>

分别为方位方向上的小转角、俯仰方向上的小转角。

进一步的，步骤S2：定义大小为PQL×PQL的置换矩阵J，其反对角线位置上的元素为1，其他位置上的元素为0，即

基于式(19)，可得到目标回波数据共轭数据信息的矩阵E_conj

E_conj＝J·X^x (20)

进一步的，步骤S3：将步骤S2中得到的目标回波数据共轭数据信息矩阵E_conj及原目标回波数据矩阵X^x的自相关矩阵进行相加，并作取平均处理，得到一个过程协方差矩阵R

式(21)中得到的过程协方差矩阵R为Hermittan矩阵，因此其满足R＝R^H。

进一步的，步骤S4：令最终总协方差矩阵R₁＝RR^H＝R² (22)

则矩阵R₁、矩阵R两者特征值与特征向量具有以下关系式：

式中，λ₁、λ分别代表矩阵R₁与R的特征值，Λ₁、Λ分别代表矩阵R₁与R的特征向量；

用R₁代替R，可增大信号特征值与噪声特征值之间的差距，并且不会改变原有的特征向量，因此在信噪比较低时，可更容易区分信号特征值与噪声特征值，可等效为增大了信噪比，有效提高了参数的估计精度。

进一步的，步骤S5：利用经典3D-ESPRIT算法对最终总协方差矩阵R₁进行处理，即可得到三维GTD散射中心模型中的各个参数。

具体的，S51：对矩阵R₁进行奇异值分解，则得到

其中，U_xS，V_xS代表R₁的信号子空间，分别由R₁的前I个主左特征向量与前I个主右特征向量构成；其中U_xN，V_xN代表R₁的噪声子空间，分别由R₁的非主左特征向量与非主右特征向量构成；

S52：构造矩阵F_x，如下式(25)所示，

式中，U _xS、

为矩阵U_xS去除后Q×L行、去除前Q×L行得到的矩阵，/>

代表U _xS的广义逆；

S53：根据三维条件下置换矩阵之间的关系求得y方向、z方向对应的信号子空间；由于三维条件下的置换矩阵E_xy、E_yz分别为

式中，

表示Kronecker乘积，/>

代表(q，l)位置上元素为1，其他位置上元素为0的Q×L矩阵，/>

代表(l，p)位置上元素为1，其他位置上元素为0的L×P矩阵，/>

代表(p，q)位置上元素为1，其他位置上元素为0的P×Q矩阵；

根据根据三维条件下置换矩阵之间的关系，得到y方向、z方向对应的信号子空间

U_ys＝E_xyU_xS (28)

U_zs＝E_yzU_yS (29)

式中，E_xy、E_yz分别为三维条件下的置换矩阵；

进而得到y方向和z方向的矩阵F_y、F_z

式中，U _vS，

为矩阵U_vS去除后Q×L行、去除前Q×L行得到的矩阵；U _zx，/>

为矩阵U_zS去除后Q×L行、去除前Q×L行得到的矩阵；/>

分别代表U _vS和U _zS的广义逆；

S54：计算矩阵F_x、F_y、F_z前I个元素的主特征值向量Ψ_x、Ψ_y、Ψ_z

Ψ_x＝T_xF_xT_x ^-1 (32)

Ψ_y＝T_yF_yT_y ^-1 (33)

Ψ_z＝T_zF_zT_z ^-1 (34)

S55：利用求得的矩阵Ψ_x、Ψ_y、Ψ_z求解其对应主对角线上的元素P_xi、P_yi和P_zi

P_xi＝diag(Ψ_x)，i＝1，K，I (35)

P_yi＝diag(Ψ_y)，i＝1，K，I (36)

P_zi＝diag(Ψ_z)，i＝1，K，I (37)

S56：根据求得的P_xi、P_yi和P_zi，将其代入式(12)-(15)，求解散射中心模型中的散射中心点的类型参数α_i，横向距离参数x_i、纵向距离参数y_i与垂直距离参数z_i；

S57：根据类型参数α_i，横向距离参数x_i、纵向距离参数y_i与垂直距离参数z_i，利用最小二乘法，求得散射中心模型中的散射中心点的强度参数

式中：

G＝[a₁，…，a_I] (39)

其中，G^H代表G的转置，(G^HG)^-1代表矩阵GHG的共轭转置。

仿真试验：

利用四个散射中心合成目标的电磁散射回波数据，来验证本发明中提出的方法对三维GTD散射中心模型的参数的估计性能。首先产生目标的后向电磁散射数据：设定雷达的初始工作频率f₀为10GHz，频率步进为16MHz，总频率步进数M＝11；初始方位角θ₀、初始俯仰角

均为90度，角度间隔均为0.01度，总方位角度数N＝11，总俯仰角度数K＝11。四个散射中心的各模型参数如下表1所示，且在回波数据中加入-10dB-20dB的信噪比SNR，其中SNR定义如下：

式中，

表1散射中心参数

为比较在不同信噪比下，本发明中提出的改进算法与其他算法对参数进行估计后的参数估计性能，分别在-10dB-20dB进行200次蒙特卡洛实验并比较不同算法估计得到的参数的均方差(root-mean-errors，RMSE)，仿真结果如附图1-5所示。其中，定义均方差RMSE如下：

式中：z_i表示第i次仿真实验估计得到的散射中心参数，z表示真实设定的散射中心参数，D表示每个信噪比对应的蒙特卡洛实验次数。

由附图1-5可知，本发明中改进的3D-ESPRIT算法估计得到的各散射中心模型参数的均方差均小于经典3D-ESPRIT算法及文献《光学区雷达目标散射中心提取及其应用研究[D]》中的改进的免配对3D-ESPRIT算法。并且可以发现，在信噪比为-10dB-0dB的仿真条件下，本发明中提出的改进算法优势更为明显，因此可以验证通过增大原始回波数据的利用率及对总协方差矩阵取平方处理方法，可有效提高算法的噪声鲁棒性与参数估计性能。

进一步的，为探究本发明中改进算法的参数估计分辨率性能，设置散射中心位置参数为：1.0，1.1，2.5，2.9和1.0，1.2，2.5，2.9两种情况，散射中心类型与强度参数设置同表1。

首先，将式(1)作简化处理得到一维GTD散射中心模型；进而可得目标的后向电场回波数据，从而利用传统逆傅立叶变换方法变换得到目标一维距离像；最后比较目标的一维距离像与改进算法得到的各散射中心位置、幅度，结果如附图6和附图7所示。

从附图6中可以看出，当4个散射中心其中的2个距离为0.1m时，基于逆傅立叶变换方法只可分辨出3个散射中心，无法分辨出全部的4个散射中心，而本发明的改进算法则可较好的分辨出全部的四个散射中心。

从附图7中可以看出，当4个散射中心其中的2个距离为0.2m时，逆傅立叶变换方法和本发明改进算法均可分辨出4个散射中心。

由此可见，基于本发明中三维GTD散射中心模型的中心参数研究方法得出的参数分辨率要优于传统的逆傅立叶变换方法。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (2)

1.一种三维GTD散射中心模型的散射中心点的参数估计方法，其特征在于，包括以下步骤，

S1：基于雷达目标后向电磁数据构建Hankel矩阵，得到原目标回波数据矩阵X^x；

S2：定义大小为PQL×PQL的置换矩阵J，其反对角线位置上的元素为1，其他位置上的元素为0，即

则可得到目标回波数据共轭数据信息的矩阵E_conj，E_conj＝J·X^x；

S3：将步骤S2中得到的目标回波数据共轭数据信息矩阵E_conj及原目标回波数据矩阵X^x的自相关矩阵进行相加，并作取平均处理，得到一个过程协方差矩阵R

其中，R＝R^H；

S4：令最终总协方差矩阵R₁＝RR^H＝R²，则矩阵R₁、矩阵R两者特征值与特征向量具有以下关系式：

式中，λ₁、λ分别代表矩阵R₁与R的特征值，Λ₁、Λ分别代表矩阵R₁与R的特征向量；

S5：利用经典3D-ESPRIT算法对最终总协方差矩阵R₁进行处理，即可得到三维GTD散射中心模型中的散射中心点的各个参数；

其中，步骤S5的具体操作包括，

S51：对最终总协方差矩阵R₁进行奇异值分解，则得到

其中，U_sS，V_xS代表R₁的信号子空间，分别由R₁的前I个主左特征向量与前I个主右特征向量构成；其中U_xN，U_xN代表R₁的噪声子空间，分别由R₁的非主左特征向量与非主右特征向量构成；

S52：构造矩阵F_x为

式中，U _xS、

为矩阵U_xS去除后Q×L行、去除前Q×L行得到的矩阵，/>

代表U _xS的广义逆；

S53：根据三维条件下置换矩阵之间的关系，得到y方向、z方向对应的信号子空间

U_ys＝E_xyU_xS

U_zs＝E_yzU_yS；

式中，E_xy、E_yz分别为三维条件下的置换矩阵；

进而得到y方向和z方向的矩阵F_y、F_z

式中，U _vS，

为矩阵U_vS去除后Q×L行、去除前Q×L行得到的矩阵；U _zS，/>

为矩阵U_zS去除后Q×L行、去除前Q×L行得到的矩阵；/>

分别代表U _vS和U _zS的广义逆；

S54：计算矩阵F_x、F_y、F_z前I个元素的主特征值向量Ψ_x、Ψ_y、Ψ_z

Ψ_x＝T_xF_xT_x ^-1

Ψ_y＝T_yF_yT_y ^-1

Ψ_z＝T_zF_zT_z ^-1；

S55：利用求得的矩阵Ψ_x、Ψ_y、Ψ_z求解其对应主对角线上的元素P_xi、P_yi和P_zi

P_xi＝diag(Ψ_x)，i＝1，K，I

P_yi＝diag(Ψ_y)，i＝1，K，I

P_zi＝diag(Ψ_z)，i＝1，K，I；

S56：根据求得的P_xi、P_yi和P_zi，求解类型参数α_i，横向距离参数x_i、纵向距离参数y_i与垂直距离参数z_i；

S57：根据类型参数α_i，横向距离参数x_i、纵向距离参数y_i与垂直距离参数z_i，利用最小二乘法，求得散射中心模型中的强度参数

式中，G＝[a₁，...，a_I]，

其中，G^H代表G的转置，(G^HG)^-1代表矩阵G^HG的共轭转置；

步骤S53中三维条件下置换矩阵E_xy、E_yz分别为

式中，

表示Kronecker乘积，/>

代表(q，l)位置上元素为1，其他位置上元素为0的Q×L矩阵，/>

代表(l,p)位置上元素为1，其他位置上元素为0的L×P矩阵，/>

代表(p,q)位置上元素为1，其他位置上元素为0的P×Q矩阵。

2.根据权利要求1所述的一种三维GTD散射中心模型的散射中心点的参数估计方法，其特征在于：步骤S1中所述的原目标回波数据矩阵X^x大小为

PQL×(M-P+1)(N-Q+1)(K-L+1)，

且

式中，

其中，M/2≤P≤2M/3，N/2≤Q≤2N/3，K/2≤L≤2K/3；

表示目标后向电磁散射回波，f_m＝f₀+mΔf，m＝1，2，...，M，f₀为起始频率，Δf为步进频率，m代表频率下标，M代表总频率个数；θ_n＝θ₀+nΔθ，n＝1，2，...，N，其中θ₀为起始方位角，Δθ为步进方位角，n为方位角下标，N代表总方位角个数；

其中/>

为起始俯仰角，/>

为步进俯仰角，k为俯仰角下标，K为总俯仰角个数；nΔθ、/>

分别为方位方向上的小转角、俯仰方向上的小转角。

Images