在高频区，雷达目标的后向电磁散射回波可视为有限个强散射点的相干叠加，而这些强散射点通常可称为雷达目标的散射中心^[1].而通过构建合理近精准的散射中心模型，可更有效地描述目标的散射特性，对深入研究其散射机理起到重要的作用.目前，基于几何绕射理论的(GTD，geometric theory of diffraction)散射中心模型^[1]对目标在高频区散射特性的刻画最为准确.而如何准确地从雷达回波数据提取出目标的散射中心模型参数成为研究热点.目前存在的主流算法有多重信号分类(MUSIC，multiple signal classification)算法^[2-7]、旋转不变子空间(ESPRIT，estimating signal parameter via rotational invariance techniques)算法^[8-15]、增广矩阵束(MEMP，matrix enhancement matrix pencil)^[16]算法等，通过上述算法即可对GTD散射中心模型进行参数估计提取.

利用三维GTD模型对雷达目标的电磁散射数据进行描述.针对经典3D-ESPRIT算法存在的参数失配问题，王菁^[15]提出一种改进的免配对3D-ESPRIT算法，有效解决了参数失配问题，但在信噪比情况下，改进免配对3D-ESPRIT算法的参数估计较低.为解决免配对3D-ESPRIT算法在低信噪比时，参数估计精度降低、算法鲁棒性较差这一问题，在免配对3D-ESPRIT算法的基础上，笔者提出了一种改进的3D-ESPRIT算法.改进算法通过构建原始回波数据的共轭矩阵、对协方差矩阵叠加、取平均、平方处理，提高了雷达目标原始回波数据的利用率，从而有效地提高了算法的参数估计性能和噪声鲁棒性.

1 信号模型

GTD模型作为一种经典的散射中心模型，可有效描述雷达目标在高频区的电磁后向散射特性，以步进频率雷达信号为例，目标的三维GTD模型可表示为

$ \begin{array}{*{20}{c}} {E({f_m},{\theta _n},{\varphi _k}) = }\\ {\sum\limits_{i = 1}^I {{A_i}} {{\left( {{\rm{j}}\frac{{{f_m}}}{f}} \right)}^{{\alpha _i}}}{\rm{exp}}[ - 4\pi {\rm{j}}{f_m}({x_i}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\kern 1pt} {\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \varphi + }\\ {{y_i}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta {\kern 1pt} {\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \varphi + {z_i}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \varphi )/c] + \omega ({f_m},{\theta _n},{\varphi _k})} \end{array} $

(1)

其中：E(f_m, θ_n, φ_k)表示目标的后向散射回波，I代表散射中心个数，{A_i, α_i, x_i, y_i, z_i}分别表示第i个散射中心的散射强度、散射类型、横向距离、纵向距离及垂直距离. f_m=f₀+mΔf，其中f₀为起始频率，Δf为步进频率，m代表频率下标；θ_n=θ₀+nΔθ，n=0, 1, …, N，其中θ₀为起始方位角，Δθ为步进方位角，n为方位角下标；θ_n=θ₀+nΔθn=0, 1, …, N，其中φ₀为起始俯仰角，Δφ为步进俯仰角，k为俯仰角下标，φ_k=φ₀+kΔφ，k=0, 1, …, K；nΔθ、kΔφ分别为方位方向上的小转角、俯仰方向上的小转角. c=3×10⁸ m/s为电磁波传播速度，ω(f_m, θ_n, φ_k)为复高斯白噪声. α_i为0.5的整数倍，共有5种^[1].式(1)中的{A_i, α_i, x_i, y_i, z_i}即为待估计散射中心参数.

2 改进算法原理

针对经典3D-ESPRIT算法在信噪比较低时参数估计性能下降这一问题，提出了一种改进的3D-ESPRIT算法，有效提高了算法的噪声鲁棒性与参数估计性能.改进算法主要对经典3D-ESPRIT算法前2个步骤进行改进，具体改进步骤如下.

步骤1  定义大小为PQL×PQL的置换矩阵J：

$ \mathit{\boldsymbol{J}} = {\left[ {\begin{array}{*{20}{c}} 0& \cdots &0&1\\ \vdots &0&1&0\\ 0& {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & \ddots & \vdots \\ 1&0& \cdots &0 \end{array}} \right]_{PQL \times PQL}} $

(2)

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

基于式(2)，可得到包含目标回波数据共轭数据信息的矩阵为

$ {\mathit{\boldsymbol{E}}_{{\rm{conj}}}}{\rm{ }} = \mathit{\boldsymbol{J}}{\mathit{\boldsymbol{X}}^x} $

(3)

其中

$ {\mathit{\boldsymbol{X}}^x} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{X}}_0^x}&{\mathit{\boldsymbol{X}}_1^x}& \cdots &{\mathit{\boldsymbol{X}}_{M - P}^x}\\ {\mathit{\boldsymbol{X}}_1^x}&{\mathit{\boldsymbol{X}}_2^x}& \cdots &{\mathit{\boldsymbol{X}}_{M - P + 1}^x}\\ \vdots & \vdots & \cdots & \vdots \\ {\mathit{\boldsymbol{X}}_{P - 1}^x}&{\mathit{\boldsymbol{X}}_P^x}& \cdots &{\mathit{\boldsymbol{X}}_{M - 1}^x} \end{array}} \right] $

(4)

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{X}}_m^x = \\ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{x}}(m,0)}&{\mathit{\boldsymbol{x}}(m,1)}& \cdots &{\mathit{\boldsymbol{x}}(m,K - L)}\\ {\mathit{\boldsymbol{x}}(m,1)}&{\mathit{\boldsymbol{x}}(m,2)}& \cdots &{\mathit{\boldsymbol{x}}(m,K - L + 1)}\\ \vdots & \vdots & \cdots & \vdots \\ {\mathit{\boldsymbol{x}}(m,L - 1)}&{\mathit{\boldsymbol{x}}(m,L)}& \cdots &{\mathit{\boldsymbol{x}}(m,K - 1)} \end{array}} \right] \end{array} $

(5)

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{x}}(m,k) = \\ \left[ {\begin{array}{*{20}{c}} {E(m,0,k)}&{E(m,1,k)}& \cdots &{E(m,N - Q,k)}\\ {E(m,1,k)}&{E(m,2,k)}& \cdots &{E(m,N - Q + 1,k)}\\ \vdots & \vdots & \cdots & \vdots \\ {E(m,Q - 1,k)}&{E(m,Q,k)}& \cdots &{E(m,N - 1,k)} \end{array}} \right] \end{array} $

(6)

通过对新得到的矩阵E_conj及原目标回波数据矩阵X^x的自相关矩阵相加，并作取平均处理，可得到一个新的总协方差矩阵R：

$ \mathit{\boldsymbol{R}} = \frac{{{\mathit{\boldsymbol{X}}^x}{{({\mathit{\boldsymbol{X}}^x})}^{\rm{H}}} + {E_{{\rm{conj}}}}{{({E_{{\rm{conj}}}})}^{\rm{H}}}}}{2} $

(7)

步骤2  式(7)得到的新总协方差矩阵为Hermittan矩阵，因此其满足R=R^H，即R₁=RR^H=R²，则矩阵R₁、R两者的特征值与特征向量具有以下关系：

$ \left. {\begin{array}{*{20}{l}} {{\lambda _1} = {\lambda ^2}}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_1} = \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \end{array}} \right\} $

(8)

其中：λ₁、λ分别代表矩阵R₁和R的特征值，Λ₁、Λ分别代表矩阵R₁和R的特征向量.

用R₁代替R，可增大信号特征值与噪声特征值之间的差距，并且不会改变原有的特征向量，因此在信噪比较低时，更容易区分信号特征值和噪声特征值.构造下式：

$ {\mathit{\boldsymbol{R}}_1} = \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{H}}} = {\mathit{\boldsymbol{R}}^2} $

(9)

对新的总协方差矩阵R₁ 3D-ESPRIT算法的三维GTD模型参数估计提取步骤，即可得到三维GTD模型中的各参数值.

提出的算法主要针对经典3D-ESPRIT算法中的前2个步骤，即协方差矩阵的构造作改进.因此，仅对2种算法的协方差矩阵构造这一步骤的运算量作比较分析，经典3D-ESPRIT算法在协方差矩阵构造这一步的运算量为P²Q²L²(M-P+1)(N-Q+1)(K-L+1)，而提出的改进3D-ESPRIT算法在此步骤的运算量为P²Q²L²[3(M-P+1)(N-Q+1)(K-L+1)+2PQL+1].由此可见，改进算法在构造协方差矩阵这一步的运算量要稍大于经典3D-ESPRIT算法.

3 仿真实验

假设目标由4个散射中心构成，具体的参数设置见表 1.仿真实验设定雷达的初始工作频率f₀为10 GHz，频率步进为16 MHz，总频率步进数M=11；初始方位角θ₀、初始俯仰角φ₀均为90°，角度间隔均为0.01°，总方位角度数N=11，总俯仰角度数K=11. 4个时散射中心的各模型参数见表 1.在回波数据中加入-10~20 dB的信噪比(SNR，signal-to -noise ratio)，SNR定义为

$ S = 10{\rm{lg}}\left( {\frac{{\sigma _n^2}}{{\sigma _s^2}}} \right) $

(10)

其中$\sigma_{s}^{2}=\sum\limits_{i=1}^{I} A_{i}^{2}$.

分别在-10~20 dB信噪比条件下进行200次蒙特卡洛实验，并比较不同算法估计得到参数的均方差(RMSE，root-mean-errors)，仿真结果如图 1~图 5所示.定义均方差RMSE为

$ E(z) = \sqrt {\frac{{\sum\limits_{i = 1}^D {{{({z_i} - z)}^2}} }}{D}} $

(11)

其中：z_i为估计参数，z为设定参数，D为蒙特卡洛次数.

由图 1~图 5可知，改进的3D-ESPRIT算法估计得到的各散射中心模型参数的均方差均小于经典3D-ESPRIT算法及文献[15]中的改进算法.可以发现，在信噪比为-10~0 dB的仿真条件下，改进的算法优势更为明显，因此可以验证，通过增大原始回波数据的利用率及对总协方差矩阵取平方处理的方法，能有效提高算法的噪声鲁棒性和参数估计性能.

4 结束语

笔者提出了一种改进3D-ESPRIT算法，该算法不仅解决了参数失配问题，而且在低信噪比的条件下参数估计性能更优，噪声鲁棒性更加稳定.所提算法在运算复杂度增加不大的条件下，解决了传统旋转不变子空间算法低信噪比时参数估计精度较低这一问题.仿真实验验证了改进算法在估计三维GTD模型参数过程中的先进性和有效性.相比于经典3D-ESPRIT算法，提出的改进算法的运算量稍有增加，因此，如何提出快速的散射中心模型参数估计算法将是此领域的重点研究方向.