基于矩阵信息几何的飞机尾流目标检测方法

刘俊凯; 李健兵; 马梁; 陈忠宽; 蔡益朝

doi:10.12000/JR17058

基于矩阵信息几何的飞机尾流目标检测方法

DOI: 10.12000/JR17058

①.
空军预警学院武汉 430019
②.
国防科技大学电子科学与工程学院长沙 410073

基金项目: 国家自然科学基金（61302193，61401503）

详细信息

作者简介:
刘俊凯(1979–)，男，河北景县，博士，讲师，主要研究方向为新体制雷达探测技术、目标检测、相控阵雷达建模仿真。E-mail: liujkradar@163.com

李健兵(1979–)，男，博士，副研究员，硕士生导师，IEEE Senior Member，中国电子学会高级会员，主要研究方向为大尺度分布式复杂目标的雷达特性与探测。E-mail: jianbingli@nudt.edu.cn

马梁：马梁(1983–)，男，博士，讲师，主要研究方向为极化信息处理、目标识别、相控阵雷达建模仿真

通讯作者:
刘俊凯 liujkradar@163.com

计量
- 文章访问数: 2526
- HTML全文浏览量: 602
- PDF下载量: 524
- 被引次数: 11
出版历程
- 收稿日期: 2017-06-15
- 修回日期: 2017-07-24
- 网络出版日期: 2017-12-28

Radar Target Detection Method of Aircraft Wake Vortices Based on Matrix Information Geometry

①.
Air Force Early Warning Academy, Wuhan 430019, China
②.
School of Electronic Science and Engineering, National University of Denfense Technology, Changsha 410073, China

Funds: The National Natural Science Foundation of China (61302193, 61401503)

摘要

摘要:
矩阵信息几何在雷达信号处理和目标检测中的应用是一个正在引起关注的研究方向。飞机尾流回波经过傅里叶变换后，其功率谱是展宽的，传统动目标检测(MTD)方法未能对展宽的功率谱进行有效积累。针对飞机尾流目标检测问题，基于矩阵信息几何理论，该文提出了一种矩阵恒虚警率(CFAR)检测方法，该方法中观测数据协方差矩阵构成一个矩阵流形，类比CFAR检测的思想，利用检测单元协方差矩阵与参考单元协方差矩阵均值间定义的距离作为检测统计量。最后利用噪声中仿真的尾流回波数据，分析了黎曼均值的迭代估计性能、尾流目标协方差矩阵与噪声协方差矩阵的测地线距离随信噪比的变化，比较了常规MTD检测方法和矩阵CFAR检测方法的检测性能。
- 飞机尾流 /
- 目标检测 /
- 矩阵信息几何 /
- 矩阵流形 /
- 矩阵CFAR
Abstract:
The application of matrix information geometry to radar signal processing and target detection is a new and interesting subject. Wake vortices are Doppler-spread after Fourier transform. The traditional Moving Target Detection (MTD) method cannot adequately accumulate returns power of the whole spectrum. Based on matrix information geometry, a matrix Constant False Alarm Rate (CFAR) detection method is proposed to improve the detection performance of a weak wake target. In this method, covariance matrices of the observed data can be constructed into a matrix manifold; compared with CFAR detection, the geodesic distance between the covariance matrix in the detection cell and the mean of covariance matrices in the secondary cell is regarded as the detection statistics. Using simulated wake vortices, the return data in background noise and the iterative estimation performance of Riemannian mean are analyzed; the geodesic distance of covariance matrices of target return and noise with varying signal-noise rate is analyzed; and the detection performance of the matrix CFAR and the conventional MTD method is compared.
- Aircraft wake vortices /
- Target detection /
- Matrix information geometry /
- Matrix manifold /
- Matrix CFAR

HTML全文

1. 引言

混合源定位在无源雷达中具有广泛的应用前景^[1,2]。按照空间中的辐射源与阵列之间的相对距离，可以将辐射源分为远场源和近场源。对远场源定位，需要对波达方向(Direction Of Arrival, DOA)进行估计；对近场源定位，除了对波达方向进行估计，还需要对距离参数进行估计。阵列的近场区域小于 $2{D^2}/\lambda$ ，其中D为圆阵孔径， $\lambda$ 为波长^[3,4]。虽然远场源可以视为距离为无穷远的近场源，但是当远场源的二维DOA与近场源二维DOA相同时，如果采用近场和远场参数联合估计将无法对混合源进行识别。此外，由于远场源定位场景和近场源定位场景可以看作混合源定位场景特殊模式，因此混合源定位方法同样适用于辐射源全部为远场源定位场景或者辐射源全部为近场源定位场景。

不同于单独对远场源进行定位的方法和单独对近场源进行定位的方法，对混合源进行定位还需要对近场源和远场源进行分类和识别。文献[5]首先利用距离为无穷远时的近场源导向矢量和混合源的噪声子空间，通过一维多重信号分类(One Dimensional MUltiple SIgnal Classification, 1-D MUSIC)方法估计出远场源DOA，接着利用对称阵元的2阶统计量去除近场源的距离参数，估计出混合源DOA并通过斜投影算法对混合源进行分类，最后将近场源DOA代入后导向矢量和混合源噪声子空间，通过1-D MUSIC方法对近场源距离进行估计。文献[6]提出了基于子空间的混合源定位算法，该算法利用线性变换估计出远场源的DOA，通过正交三角分解(QR分解)得到斜投影算子避免了特征值分解，大大减少了计算复杂度，利用1-D MUSIC方法确定近场源DOA，最后将近场源DOA代入导向矢量，确定近场源的距离参数。该方法有效解决了在快拍数有限时出现饱和的情况。文献[7]利用两级对称嵌套线阵，利用距离为无穷远的近场源导向矢量和混合源的噪声子空间，通过1-D MUSIC方法估计出远场源DOA，接着利用斜投影方法去除混合源接收信号中的远场分量，计算近场源在对称阵元下的4阶累积量并进行向量化处理，估计出近场源DOA，最后利用代入近场源DOA的导向矢量和混合源噪声子空间，确定近场源的距离参数。该方法能够在阵元数相同的情况下增大阵列孔径，提高混合源参数估计精度。上述3种方法在得到远场源的波达方向的基础上，利用斜投影算法将远场源和近场源进行分离，但是在远场源参数估计的基础上利用斜投影算法得到的近场源分量会受到噪声分量的影响，引起近场源参数估计误差，导致定位精度下降。文献[8]利用协方差矩阵差分方法消除混合源中远场源和噪声分量，得到近场源差分矩阵，该方法能够有效避免噪声分量对近场源定位精度的影响，文献[9]采用降秩的方法从混合源中分离出远场源分量，文献[10]设计了3个4阶累积量矩阵消除近场源导向矢量中距离参数，进而同时估计出远场源和近场源的波达方向。但是上述方法采用线阵对混合源进行定位只能估计出一维波达方向，并且线阵在0°波达方向的分辨率最高，在60°波达方向的分辨率会降低1/2。

相较于均匀线阵只能对一维波达方向进行估计，均匀圆阵可以实现方位角和俯仰角等二维波达方向估计^[11-13]，并且在不同方位角的分辨率具有各向同性^[14,15]。在现有文献中，由于均匀圆阵的导向矩阵不具有均匀线阵导向矩阵的Vandermonde形式，不便于在数学上进行处理，因此提出的方法相对较少。均匀圆阵下对混合源进行定位需要针对相应的数学模型进行分析，文献[16]采用均匀圆阵，提出了TSMUSIC方法，该方法利用均匀圆阵下的远场源的协方差矩阵具有Hermitian和Toeplitz形式，而近场源的协方差矩阵只具有Hermitian形式，利用传播算子的正交性得到近场源的噪声子空间，进而通过2-D MUSIC方法得到近场源二维DOA的空间谱，进而根据谱峰确定近场源二维DOA；接着，利用将近场源二维DOA代入后的导向矢量和混合源的噪声子空间，通过1-D MUSIC方法确定近场源的距离；最后，利用远场源的导向矢量和混合源的噪声子空间确定远场源的二维DOA。该方法参数估计精度较高，但是计算复杂度较大。文献[17]采用均匀圆阵，利用每个阵元的频谱计算非相干的混合源在每个阵元的相位，利用均匀圆阵中心对称的特点，通过计算对角阵元的相位差来去除近场源的距离参数，得到混合源二维DOA的相位差矩阵，进而运用最小二乘法同时估计出远场源和近场源的方位角与俯仰角；接着，根据1-D MUSIC方法得到距离空间谱的收敛性，对混合源进行识别并估计出近场源的距离参数，当距离空间谱出现峰值，该混合源判断为近场源，峰值所对应的位置即为近场源的距离；若距离空间谱不收敛，则该混合源判断为远场源。基于相位差的方法能够有效减少计算复杂度，但是混合源的定位精度有所降低。

针对上述方法存在的问题，本文根据均匀圆阵中心对称的特点，提出一种基于协方差矩阵差分的混合源定位方法，该方法首先利用二维多重信号分类方法估计出远场源的方位角和俯仰角；接着利用协方差矩阵差分方法提取出近场源差分矩阵，通过改进类旋转不变估计信号参数(Estimation of Signal Parameters via Rotational Invariance Techniques like, ESPRIT-like)方法计算出近场源的方位角和俯仰角；最后，利用一维多重信号分类方法估计出近场源的距离。该方法能够消除混合源协方差矩阵中的远场源的噪声分量，从而提取出近场源差分矩阵，并且相较于基于相位差的方法，本文提出的方法能够提高混合源的定位精度。当远场源的二维DOA与近场源二维DOA相同时，如果采用近场和远场参数联合估计将无法对混合源进行识别，采用本文所提的两步法可以将远场源和近场源进行有效分离，在混合源识别上更具优势。

2. 混合源模型

均匀圆阵下的混合源模型如图1所示，以均匀圆阵的圆心为原点建立三维坐标系，均匀圆阵的半径为 $R$ ，由分布在xy平面上的M个阵元组成，混合源包含K₁个远场源和K₂个近场源，方位角参数 $\phi$ 为混合源与均匀圆阵中心的连线投影到xy平面上，相对于x坐标轴的逆时针方向旋转的角度；俯仰角参数 $\theta$ 为混合源与均匀圆阵中心的连线相对于z坐标轴旋转的角度；当混合源为近场源时，距离参数 $r$ 为近场源与均匀圆阵中心的距离。

图 1 混合源模型

Figure 1. The geometric model of mixed sources

下载: 全尺寸图片幻灯片

均匀圆阵的第 $m$ 个阵元在第 $n$ 个快拍的接收数据可以表示为

${x_m}(n) = \sum\limits_{k = 1}^{{K_1}{\rm{ + }}{K_2}} {{s_k}(n){{\rm{e}}^{{\rm{j}}{\tau _{k,m}}}}} + {w_m}(n)$

(1)

其中， ${s_k}(n)$ 表示第 $k$ 个混合源发射数据， ${w_m}(n)$ 表示零均值复白高斯过程， ${\tau _{k,m}}$ 表示第 $k$ 个混合源到达阵列中心以及到达第 $m$ 个阵元的时差，当混合源为远场源时，时差可以表示为

${\tau _{k,m}}=\frac{{2{\rm{\pi }}}}{\lambda }{\zeta _{k,m}}$

(2)

其中， $k = 1,2, \cdots ,{K_1}$ ， $\lambda$ 表示波长， $2{\rm{\pi /}}\lambda$ 表示波数， ${\zeta _{k,m}} = \cos ({\varphi _m} - {\phi _k})\sin {\theta _k}$ $,\;{\varphi _m} = 2{\rm{\pi }}(m - 1)/ M$ ，当混合源为近场源时，时差可以表示为

${\tau _{k,m}}=\frac{{2{\rm{\pi }}}}{\lambda }\left( {{r_k}{\rm{ - }}{r_{k,m}}} \right)$

(3)

其中， $k = {K_1} + 1,{K_1} + 2, \cdots ,{K_1} + {K_2}$ , ${r_{k,m}}$ 表示第 $k$ 个近场源与第 $m$ 个阵元的距离， ${r_{k,m}}$ 通过2阶泰勒级数^[11]可以展开为

$\begin{split}{r_{k,m}}& = \sqrt {\mathop {{r_k^2 }} + \mathop R\nolimits^2 - 2{r_k}R{\zeta _{k,m}}} \\ & ={r_k} - R{\zeta _{k,m}} + \frac{{{R^2}}}{{2{r_k}}}(1 - \zeta _{k,m}^2) + O\left( {\frac{{{R^2}}}{{r_k^2}}} \right)\end{split}$

(4)

将式(2)和式(4)代入式(1)，混合源模型可以简化为

$\begin{split} {x_m}(n) = &\sum\limits_{k = 1}^{{K_1}} {{s_k}(n){{\rm{e}} ^{{\rm{j}}{\eta _{k,m}}}}} \\ &+ \sum\limits_{k = {K_1} + 1}^{{K_1} + {K_2}} {{s_k}(n){{\rm{e}} ^{{\rm{j}}{\eta _{k,m}}{\rm{ - j}}{\xi _{k,m}}}}} + {w_m}(n) \end{split}$

(5)

其中， ${\eta _{k,m}} = \left( {2{\rm{\pi }}R/\lambda } \right){\zeta _{k,m}}$ , ${\xi _{k,m}} = \left( {{\rm{\pi }}{R^2}/\lambda {r_k}} \right)$ $\cdot \left( {1 - \zeta _{k,m}^2} \right)$ 。

3. 算法描述

在均匀圆阵下的近场和远场混合源数学模型的基础上，本节首先利用2-D MUSIC方法估计出远场源的方位角和俯仰角，实现远场源定位；接着利用协方差矩阵差分方法提取出近场源差分矩阵，进而通过改进的类旋转不变估计信号参数方法计算出近场源的方位角和俯仰角；最后利用1-D MUSIC方法估计出近场源的距离，实现近场源定位。

3.1 远场源二维DOA参数估计

均匀圆阵下混合源的协方差矩阵 ${\boldsymbol{E}}$ 可以由式(6)计算

${\boldsymbol{E}}=\frac{1}{N}{\boldsymbol{X}}{{\boldsymbol{X}}^{\rm{H}}}$

(6)

其中， ${\boldsymbol{X}}=\left[ {{\boldsymbol{x}}(1)\;{\boldsymbol{x}}(2)\; \cdots \;{\boldsymbol{x}}(N)} \right]$ 表示 $M \times N$ 维的混合源接收数据， ${\boldsymbol{x}}(n) = {\left[ {{x_1}(n)\;{x_2}(n)\; \cdots \;{x_M}(n)} \right]^{\rm{T}}}$ 表示第 $n$ 个快拍接收数据， ${( \cdot )^{\rm{H}}}$ 表示共轭转置，N为混合源的快拍数， ${( \cdot )^{\rm{T}}}$ 表示转置运算。

对混合源的协方差矩阵 ${\boldsymbol{E}}$ 进行特征值分解得到

${\boldsymbol{E}}={{\boldsymbol{U}}_{\rm{s}}}{{\boldsymbol{\varSigma}} _{\rm{s}}}{\boldsymbol{U}}_{\rm{s}}^{\rm{H}} + {{\boldsymbol{U}}_{\rm{n}}}{{\boldsymbol{\varSigma}} _{\rm{n}}}{\boldsymbol{U}}_{\rm{n}}^{\rm{H}}$

(7)

其中， ${{\boldsymbol{\varSigma}}_{\rm{s}}} \in {\mathbb{C}^{({K_1} + {K_2}) \times ({K_1} + {K_2})}}$ 表示 ${K_1} + {K_2}$ 个大特征值组成的对角矩阵， ${{\boldsymbol{\varSigma}}_{\rm{n}}} \in {\mathbb{C}^{(M - {K_1} - {K_2}) \times (M - {K_1} - {K_2})}}$ 表示 $M - {K_1} - {K_2}$ 个小特征值组成的对角矩阵， ${{\boldsymbol{U}}_{\rm{s}}}\in$ ${\mathbb{C}^{M \times ({K_1} + {K_2})}}$ 表示 ${K_1} + {K_2}$ 个大特征值对应的特征向量组成的混合源信号子空间， ${{\boldsymbol{U}}_{\rm{n}}} \in {\mathbb{C}^{M \times (M - {K_1} - {K_2})}}$ 表示 $M - {K_1} - {K_2}$ 个小特征值对应的特征向量组成的混合源噪声子空间。

利用3-D MUSIC方法得到混合源的空间谱函数

${\boldsymbol{V}}\left( {\phi ,\theta ,r} \right) = \frac{1}{{{\boldsymbol{a}}_{{\rm{NF}}}^{\rm{H}}\left( {\phi ,\theta ,r} \right){{\boldsymbol{U}}_{\rm{n}}}{\boldsymbol{U}}_{\rm{n}}^{\rm{H}}{{\boldsymbol{a}}_{{\rm{NF}}}}\left( {\phi ,\theta ,r} \right)}}$

(8)

其中， ${{\boldsymbol{a}} _{{\rm{NF}}}}(\phi ,\theta ,r) = {[{a_{{\rm{NF,1}}}}(\phi ,\theta ,r)\;\;{a_{{\rm{NF}},2}}(\phi ,\theta ,r)\; \cdots}$ ${a_{{\rm{NF}},M}} (\phi ,\theta ,r)]^{\rm{T}}$ 为近场源的导向矢量， ${a_{{\rm{NF}},m}} \left( {\phi ,\theta ,r} \right) =$ ${{\rm{e}}^{{\rm{j}}(2{\rm{\pi }}R/\lambda ){\zeta _m}(\phi ,\theta ) - {\rm{j}}\left( {{\rm{\pi }}{R^2}/\lambda r} \right)\left( {1 - \zeta _m^2(\phi ,\theta )} \right)}}$ ， $m =$ $1,2, \cdots ,M$ 。混合源空间谱函数 ${\boldsymbol{V}}\left( {\phi ,\theta ,r} \right)$ 会出现 $K$ 个峰值，其中 ${K_1}$ 个峰值的距离估计 $\hat r \to \infty$ ，表示 ${K_1}$ 个远场源； ${K_2}$ 个峰值的距离估计 $r \in \left[ {0.62{{(4{R^2}/\lambda )}^{0.5}},8{R^2}/\lambda } \right]$ ，表示 ${K_2}$ 个近场源。因此，远场源二维DOA的空间谱函数可以通过2-D MUSIC方法得到

${{\boldsymbol{V}}_{{\rm{FF}}}}\left( {\phi ,\theta } \right) = \frac{1}{{{\boldsymbol{a}}_{{\rm{NF}}}^{\rm{H}}\left( {\phi ,\theta ,\infty } \right){{\boldsymbol{U}}_{\rm{n}}}{\boldsymbol{U}}_{\rm{n}}^{\rm{H}}{{\boldsymbol{a}}_{{\rm{NF}}}}\left( {\phi ,\theta ,\infty } \right)}}$

(9)

其中， ${{\boldsymbol{a}} _{{\rm{NF}}}}(\phi ,\theta ,\infty ) = {[{a_{{\rm{NF,1}}}}(\phi ,\theta ,\infty )\;{a_{{\rm{NF}},2}}(\phi ,\theta ,\infty )\; \cdots}$ ${\;{a_{{\rm{NF}},M}}(\phi ,\theta ,\infty )]^{\rm{T}}}$ 表示距离为无穷远情况下的近场源导向矢量，即远场源导向矢量，此时 ${a_{{\rm{NF}},m}}( {\phi ,\theta ,} {\infty } )$ $= {{\rm{e}}^{{\rm{j}}(2{\rm{\pi }}R/\lambda ){\zeta _m}(\phi ,\theta )}}$ , $m = 1,2, \cdots ,M$ 。只包含方位角和俯仰角等二维DOA。对远场源方位角和俯仰角的空间谱函数 ${{\boldsymbol{V}}_{{\rm{FF}}}}\left( {\phi ,\theta } \right)$ 进行谱峰搜索，峰值所对应的位置即为远场源的方位角估计值和俯仰角估计值，则第 $k$ 个远场源的定位结果为 $({\hat \phi _k},{\hat \theta _k},\infty )$ , $k = 1,2, \cdots ,{K_1}$ 。

3.2 混合源协方差矩阵差分

混合源协方差矩阵 ${\boldsymbol{E}}$ 可以表示为

$\begin{split} {\boldsymbol{E}}=&{{\boldsymbol{E}}_{{\rm{FF}}}}+{{\boldsymbol{E}}_{{\rm{NF}}}}+{\sigma ^2}{{\boldsymbol{I}}_M} \\ =&{{\boldsymbol{A}}_{{\rm{FF}}}}{{\boldsymbol{\varPi}} _{{\rm{FF}}}}{\boldsymbol{A}}_{{\rm{FF}}}^{\rm{H}}+{{\boldsymbol{A}}_{{\rm{NF}}}}{{\boldsymbol{\varPi}} _{{\rm{NF}}}}{\boldsymbol{A}}_{{\rm{NF}}}^{\rm{H}}+{\sigma ^2}{{\boldsymbol{I}}_M}\end{split}$

(10)

其中， ${{\boldsymbol{\varPi}} _{{\rm{FF}}}} \in {\mathbb{C}^{{K_1} \times {K_1}}}$ 表示远场源功率组成的对角矩阵， ${{\boldsymbol{\varPi}} _{{\rm{NF}}}} \in {\mathbb{C}^{{K_2} \times {K_2}}}$ 表示近场源功率组成的对角矩阵， ${\sigma ^2}$ 表示噪声功率， ${{\boldsymbol{I}}_M}$ 表示 $M \times M$ 维的单位矩阵， ${{\boldsymbol{A}}_{{\rm{FF}}}}$ 表示远场源的导向矩阵

${{\boldsymbol{A}}_{{\rm{FF}}}}=\big[{{\boldsymbol{a}}_{{\rm{FF}}}}({\phi _1},{\theta _1})\;{{\boldsymbol{a}}_{{\rm{FF}}}}({\phi _2},{\theta _2})\; \cdots \;{{\boldsymbol{a}}_{{\rm{FF}}}}({\phi _{{K_1}}},{\theta _{{K_1}}})\big]$

(11)

由于均匀圆阵具有中心对称的结构，满足 ${\varphi _{m + M/2}} = {\varphi _m} + {\rm{\pi }}$ 和 ${\zeta _{m + M/2}} = - {\zeta _m}$ ，因此，远场源的导向矢量可以表示为

$\begin{split} {{\boldsymbol{a}}_{{\rm{FF}}}}({\phi _k},{\theta _k}) = &[\underbrace {{{\rm{e}}^{{\rm{j}}{\eta _{k,1}}}}\;{{\rm{e}}^{{\rm{j}}{\eta _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{j}}{\eta _{k,M/2}}}}}_{M/2}\;\\ &\underbrace {{{\rm{e}}^{{\rm{ - j}}{\eta _{k,1}}}}\;{{\rm{e}}^{{\rm{ - j}}{\eta _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{ - j}}{\eta _{k,M/2}}}}}_{M/2}{]^{\rm{T}}} \end{split}$

(12)

同理， ${{\boldsymbol{A}}_{{\rm{NF}}}}$ 表示近场源的导向矩阵为

$\begin{split} {{\boldsymbol{A}}_{{\rm{NF}}}}=& [{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _1},{\theta _1},{r_1})\;{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _2},{\theta _2},{r_2})\; \cdots\\ &\;{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _{{K_2}}},{\theta _{{K_2}}},{r_{{K_2}}})]\end{split}$

(13)

其中，近场源的导向矢量可以表示为

$\begin{split} {{\boldsymbol{a}}_{{\rm{NF}}}}&({\phi _k},{\theta _k},{r_k}) \\ =&[\underbrace {{{\rm{e}}^{{\rm{j}}{\eta _{k,1}}{\rm{ - j}}{\xi _{k,1}}}}\;{{\rm{e}}^{{\rm{j}}{\eta _{k,2}}{\rm{ - j}}{\xi _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{j}}{\eta _{k,M/2}}{\rm{ - j}}{\xi _{k,M/2}}}}}_{M/2}\;\\ &\underbrace {{{\rm{e}}^{{\rm{ - j}}{\eta _{k,1}}{\rm{ - j}}{\xi _{k,1}}}}\;{{\rm{e}}^{{\rm{ - j}}{\eta _{k,2}}{\rm{ - j}}{\xi _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{ - j}}{\eta _{k,M/2}}{\rm{ - j}}{\xi _{k,M/2}}}}}_{M/2}{]^{\rm{T}}}\end{split}$

(14)

由于远场源的导向矢量中第 $m$ 个元素和第 $m +$ $M/2$ 个元素具有共轭的结构，因此远场源的协方差矩阵 ${{\boldsymbol{E}}_{{\rm{FF}}}}$ 可以表示为

${\boldsymbol{E}}_{{\rm{FF}}}^{\rm{T}}={\boldsymbol{J}}{{\boldsymbol{E}}_{{\rm{FF}}}}{\boldsymbol{J}}$

(15)

其中，旋转变换矩阵 ${\boldsymbol{J}}$ 可以表示为

${\boldsymbol{J}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{O}}_{M/2}}}&{{{\boldsymbol{I}}_{M/2}}} \\ {{{\boldsymbol{I}}_{M/2}}}&{{{\boldsymbol{O}}_{M/2}}} \end{array}} \right]$

(16)

${{\boldsymbol{O}}_{M/2}}$ 和 ${{\boldsymbol{I}}_{M/2}}$ 分别为 $M/2 \times M/2$ 维的零矩阵和单位矩阵。

由式(10)和式(15)可以推导近场源差分矩阵为

${{\boldsymbol{E}}_{{\rm{NF,diff}}}}={\boldsymbol{JEJ}}-{{\boldsymbol{E}}^{\rm{T}}}={\boldsymbol{J}}{{\boldsymbol{E}}_{{\rm{NF}}}}{\boldsymbol{J}}-{\boldsymbol{E}}_{{\rm{NF}}}^{\rm{T}}$

(17)

由于近场源协方差矩阵 ${{\boldsymbol{E}}_{{\rm{NF,diff}}}}$ 是Hermitian矩阵，满足 ${\boldsymbol{E}}_{{\rm{NF,diff}}}^{\rm{T}}={\boldsymbol{E}}_{{\rm{NF,diff}}}^{*}$ ，因此，近场源协方差矩阵 ${{\boldsymbol{E}}_{{\rm{NF,diff}}}}$ 可以重构为

$\begin{split} {{\boldsymbol{E}}_{{\rm{NF,diff}}}}=&{{\boldsymbol{JE}}_{{\rm{NF}}}}{\boldsymbol{J}}-{\boldsymbol{E}}_{{\rm{NF}}}^{*} \\ =&{{\boldsymbol{JA}}_{{\rm{NF}}}}{{\boldsymbol{\varPi}} _{{\rm{NF}}}}{\boldsymbol{A}}_{{\rm{NF}}}^{\rm{H}}{\boldsymbol{J}}-{\boldsymbol{A}}_{{\rm{NF}}}^{*}{\boldsymbol{\varPi}} _{{\rm{NF}}}^{*}{\left( {{\boldsymbol{A}}_{{\rm{NF}}}^{\rm{H}}} \right)^*} \\ = &\left[ {\begin{split} {{{\boldsymbol{JA}}_{{\rm{NF}}}}}\qquad{{\boldsymbol{A}}_{{\rm{NF}}}^{*}} \end{split}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\varPi}} _{{\rm{NF}}}}}&{\rm{}}\\ {{\rm{}}}& {-{\boldsymbol{\varPi}} _{{\rm{NF}}}^{*}} \end{array}} \right] \\ & \cdot{\left[ {\begin{split} {{{\boldsymbol{JA}}_{{\rm{NF}}}}}\qquad{{\boldsymbol{A}}_{{\rm{NF}}}^{*}} \end{split}} \right]^{\rm{H}}}\\ =& {{\boldsymbol{A}}_{{\rm{NF,diff}}}}{{\boldsymbol{\varPi}} _{{\rm{NF}}}}{\boldsymbol{A}}_{{\rm{NF,diff}}}^{\rm{H}} \\ \end{split}$

(18)

其中，近场源差分矩阵的导向矩阵可以表示为

$\begin{split} {{\boldsymbol{A}}_{{\rm{NF}},}}{\rm{_{diff}}} =& [\underbrace {{\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _1},{\theta _1},{r_1})\;{\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _2},{\theta _2},{r_2})\; \cdots \;{\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _{{K_2}}},{\theta _{{K_2}}},{r_{{K_2}}})}_{{K_2}}\\ & \underbrace {{\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _1},{\theta _1},{r_1})\;{\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _2},{\theta _2},{r_2})\; \cdots \;{\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _{{K_2}}},{\theta _{{K_2}}},{r_{{K_2}}})}_{{K_2}}] \end{split}$

(19)

其中，第 $k$ 个导向矢量可以表示为

${\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _k},{\theta _k},{r_k}) =[\underbrace {{{\rm{e}}^{{\rm{ - j}}{\eta _{k,1}}{\rm{ - j}}{\xi _{k,1}}}}\;{{\rm{e}}^{{\rm{ - j}}{\eta _{k,2}}{\rm{ - j}}{\xi _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{ - j}}{\eta _{k,M/2}}{\rm{ - j}}{\xi _{k,M/2}}}}}_{M/2}\underbrace {{{\rm{e}}^{{\rm{j}}{\eta _{k,1}}{\rm{ - j}}{\xi _{k,1}}}}\;{{\rm{e}}^{{\rm{j}}{\eta _{k,2}}{\rm{ - j}}{\xi _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{j}}{\eta _{k,M/2}}{\rm{ - j}}{\xi _{k,M/2}}}}}_{M/2}{]^{\rm{T}}}$

(20)

第 $k + {K_2}$ 个导向矢量可以表示为

${\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _k},{\theta _k},{r_k}) =[\underbrace {{{\rm{e}}^{{\rm{ - j}}{\eta _{k,1}}{\rm{ + j}}{\xi _{k,1}}}}\;{{\rm{e}}^{{\rm{ - j}}{\eta _{k,2}}{\rm{ + j}}{\xi _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{ - j}}{\eta _{k,M/2}}{\rm{ + j}}{\xi _{k,M/2}}}}}_{M/2}\underbrace {{{\rm{e}}^{{\rm{j}}{\eta _{k,1}}{\rm{ + j}}{\xi _{k,1}}}}\;{{\rm{e}}^{{\rm{j}}{\eta _{k,2}}{\rm{ + j}}{\xi _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{j}}{\eta _{k,M/2}}{\rm{ + j}}{\xi _{k,M/2}}}}}_{M/2}{]^{\rm{T}}}$

(21)

由于 ${{\boldsymbol{A}}_{{\rm{NF,diff}}}}$ 的第 $k$ 列的指数项第1部分与第 $k + {K_1}$ 列的指数项第1部分相同，各列中第 $m$ 个元素的指数项第2部分与第 $m + M/2$ 个元素的指数项第2部分相同，因此可以得到

$\begin{split} {{\boldsymbol{JA}}_{{\rm{NF,diff}}}}=[\underbrace {{\boldsymbol{D}}({\phi _1},{\theta _1}){\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _1},{\theta _1},{r_1})\;\;{\boldsymbol{D}}({\phi _2},{\theta _2}){\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _2},{\theta _2},{r_2})\;\; \cdots \;\;{\boldsymbol{D}}({\phi _{{K_2}}},{\theta _{{K_2}}}){\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _{{K_2}}},{\theta _{{K_2}}},{r_{{K_2}}})}_{{K_2}} \\ \underbrace { {\boldsymbol{D}}({\phi _1},{\theta _1}){\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _1},{\theta _1},{r_1})\;\;{\boldsymbol{D}}({\phi _2},{\theta _2}){\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _2},{\theta _2},{r_2})\;\; \cdots \;\;{\boldsymbol{D}}({\phi _{{K_2}}},{\theta _{{K_2}}}){\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _{{K_2}}},{\theta _{{K_2}}},{r_{{K_2}}})]}_{{K_2}} \end{split}$

(22)

其中， ${\boldsymbol{D}}({\phi _k},{\theta _k})$ 表示为

$\begin{split} {\boldsymbol{D}}({\phi _k},{\theta _k}) =& {\rm{diag}}\underbrace {[{{\rm{e}}^{{\rm{j}}2{\eta _{k,1}}}}\;{{\rm{e}}^{{\rm{j}}2{\eta _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{j}}2{\eta _{k,M/2}}}}}_{M{\rm{/2}}}\\ &\underbrace {{{\rm{e}}^{{\rm{ - j}}2{\eta _{k,1}}}}\;{{\rm{e}}^{{\rm{ - j}}2{\eta _{k,2}}}}\; \cdots \;{{\rm{e}}^{{\rm{ - j}}2{\eta _{k,M/2}}}}]}_{M{\rm{/2}}} \end{split}$

(23)

${\boldsymbol{D}}({\phi _k},{\theta _k})$ 只包含近场源的二维DOA参数，因此可以利用 ${\boldsymbol{D}}({\phi _k},{\theta _k})$ 将近场源三维位置参数估计分解为二维DOA估计和距离参数估计，能够有效降低计算复杂度。

3.3 近场源二维DOA参数估计

对仅包含近场源差分矩阵 ${{\boldsymbol{E}}_{{\rm{NF,diff}}}}$ 进行特征值分解

$\begin{split} {{\boldsymbol{E}}_{{\rm{NF,diff}}}}=&{{\boldsymbol{U}}_{{\rm{NF,diff,s}}}}{{\boldsymbol{\varSigma}}_{{\rm{NF,diff,s}}}}{\boldsymbol{U}}_{{\rm{NF,diff,s}}}^{\rm{H}} \\ &+{{\boldsymbol{U}}_{{\rm{NF,diff,n}}}}{{\boldsymbol{\varSigma}}_{{\rm{NF,diff,n}}}}{\boldsymbol{U}}_{{\rm{NF,diff,n}}}^{\rm{H}}\end{split}$

(24)

其中， ${{\boldsymbol{\varSigma}}_{{\rm{NF,diff}},{\rm{s}}}} \in {\mathbb{C}^{2{K_2} \times 2{K_2}}}$ 表示 $2{K_2}$ 个大特征值组成的对角矩阵， ${{\boldsymbol{\varSigma}}_{{\rm{NF,diff}},{\rm{n}}}} \in {\mathbb{C}^{(M - 2{K_2}) \times (M - 2{K_2})}}$ 表示 $M - 2{K_2}$ 个小特征值组成的对角矩阵， ${{\boldsymbol{U}}_{{\rm{NF,diff,s}}}} \in$ ${\mathbb{C}^{M \times 2{K_2}}}$ 表示 $2{K_2}$ 个大特征值对应的特征向量组成的近场源信号子空间， ${{\boldsymbol{U}}_{\rm{n}}} \in {\mathbb{C}^{M \times (M - 2{K_2})}}$ 表示M– $2{K_2}$ 个小特征值对应的特征向量组成的近场源噪声子空间。

存在一个满秩矩阵 ${\boldsymbol{G}}$ 使得

${{\boldsymbol{U}}_{{\rm{NF,diff,s}}}}={{\boldsymbol{A}}_{{\rm{NF,diff}}}}{\boldsymbol{G}}$

(25)

由式(22)和式(19)可以得到

$\begin{split} {{\boldsymbol{JU}}}&_{{{\rm{NF,diff,s}}}} - {\boldsymbol{\varPsi}} (\phi ,\theta ){{\boldsymbol{U}}_{{\rm{NF,diff,s}}}}\\ &={{\boldsymbol{JA}}_{{\rm{NF,diff}}}}{\boldsymbol{G}} - {\boldsymbol{\varPsi}} (\phi ,\theta ){{\boldsymbol{A}}_{{\rm{NF,diff}}}}{\boldsymbol{G}}=\left[ {{{\boldsymbol{\varOmega}} _1}\;{{\boldsymbol{\varOmega}} _2}} \right]{\boldsymbol{G}}\end{split}$

(26)

其中，

$\begin{split} {{\boldsymbol{\varOmega}} _1} = &[ {( {{\boldsymbol{D}}({\phi _1},{\theta _1}) - {\boldsymbol{\varPsi}} (\phi ,\theta )} ){\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _1},{\theta _1},{r_1})\;}\\ &\left( {{\boldsymbol{D}}({\phi _2},{\theta _2}) - {\boldsymbol{\varPsi}} (\phi ,\theta )} \right){\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _2},{\theta _2},{r_2})\; \cdots\\ &\left( {{\boldsymbol{D}}({\phi _{{K_2}}},{\theta _{{K_2}}}) - {\boldsymbol{\varPsi}} (\phi ,\theta )} \right){\boldsymbol{J}}{{\boldsymbol{a}}_{{\rm{NF}}}}({\phi _{{K_2}}},{\theta _{{K_2}}},{r_{{K_2}}})]\end{split}$

(27)

$\begin{split}{{\boldsymbol{\varOmega}} _2} =& [ {\left( {{\boldsymbol{D}}({\phi _1},{\theta _1}) - {\boldsymbol{\varPsi}} (\phi ,\theta )} \right){\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _1},{\theta _1},{r_1})\;}\\ &\left( {{\boldsymbol{D}}({\phi _2},{\theta _2}) - {\boldsymbol{\varPsi}} (\phi ,\theta )} \right){\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _2},{\theta _2},{r_2})\; \cdots \\ &\;\left( {{\boldsymbol{D}}({\phi _{{K_2}}},{\theta _{{K_2}}}) - {\boldsymbol{\varPsi}} (\phi ,\theta )} \right){\boldsymbol{a}}_{{\rm{NF}}}^*({\phi _{{K_2}}},{\theta _{{K_2}}},{r_{{K_2}}}) ]\end{split}$

(28)

当 ${\boldsymbol{D}}({\phi _k},{\theta _k}) = {\boldsymbol{\varPsi}} (\phi ,\theta )$ 时， ${{\boldsymbol{JU}}_{{\rm{NF,diff,s}}}} - {\boldsymbol{\varPsi}} (\phi ,\theta )$ ${{\boldsymbol{U}}_{{\rm{NF,diff,s}}}}=0$ 。因此，可以通过改进的ESPRIT-like方法^[15]得到近场源二维DOA的空间谱函数

${{\boldsymbol{V}}_{{\rm{NF}}}}\left( {\phi ,\theta } \right) = \frac{1}{{\det \left( {{{\boldsymbol{Q}}^{\rm{H}}}{{\boldsymbol{JU}}_{{\rm{NF,diff,s}}}} - {{\boldsymbol{Q}}^{\rm{H}}}{\boldsymbol{\varPsi}} (\phi ,\theta ){{\boldsymbol{U}}_{{\rm{NF,diff,s}}}}} \right)}}$

(29)

其中， $\det \left( \cdot \right)$ 表示行列式的值， ${\boldsymbol{Q}}$ 为 $M \times 2{K_2}$ 维随机满秩矩阵， ${\boldsymbol{J}}$ 为3.2节确定的 $M \times M$ 维旋转变换矩阵，对角矩阵 ${\boldsymbol{\varPsi}} (\phi ,\theta )$ 为

$\begin{split} {\boldsymbol{\varPsi}} (\phi ,\theta ) = {\rm{diag([}}\underbrace {{{\rm{e}} ^{{\rm{j}}2{\eta _1}(\phi ,\theta )}}\;{{\rm{e}} ^{{\rm{j}}2{\eta _2}(\phi ,\theta )}}\; \cdots \;{{\rm{e}} ^{{\rm{j}}2{\eta _{M/2}}(\phi ,\theta )}}}_{M/2}\;\\ \underbrace {{{\rm{e}} ^{{\rm{ - j}}2{\eta _1}(\phi ,\theta )}}\;{{\rm{e}} ^{{\rm{ - j}}2{\eta _2}(\phi ,\theta )}}\; \cdots \;{{\rm{e}} ^{{\rm{ - j}}2{\eta _{M/2}}(\phi ,\theta )}}}_{M/2}{\rm{])}}\end{split}$

(30)

${\rm{diag}}\left( \cdot \right)$ 表示将向量转化为对角矩阵运算。对近场源方位角和俯仰角的空间谱函数 ${{\boldsymbol{V}}_{{\rm{NF}}}}\left( {\phi ,\theta } \right)$ 进行谱峰搜索，峰值所对应的位置即为近场源的方位角估计值 ${\hat \phi _k}$ 和俯仰角估计值 ${\hat \theta _k}$ , $k = 1,2, \cdots ,{K_2}$ 。

3.4 近场源距离参数估计

在估计出近场源的方位角和俯仰角的基础上，通过1-D MUSIC方法得到第 $k$ 个近场源的距离空间谱函数

$\begin{split} &{{\boldsymbol{G}}_{{\rm{NF}}}}({\hat \phi _k},{\hat \theta _k},{r_k}) =\\ &\quad\frac{1}{{{{\boldsymbol{a}}_{{\rm{NF}}}}\left({{\hat \phi }_k},{{\hat \theta }_k},{r_k}\right){{\boldsymbol{U}}_{\rm{n}}}{\boldsymbol{U}}_{\rm{n}}^{\rm{H}}{{\boldsymbol{a}}_{{\rm{NF}}}}\left({{\hat \phi }_k},{{\hat \theta }_k},{r_k}\right)}} \end{split}$

(31)

其中， ${{\boldsymbol{U}}_{\rm{n}}}$ 为混合源噪声子空间， ${{\boldsymbol{a}}_{{\rm{NF}}}}({\hat \phi _k},{\hat \theta _k},{r_k})$ 表示第 $k$ 个近场源的方位角估计值 ${\hat \phi _k}$ 和俯仰角估计值 ${\hat \theta _k}$ 代入近场源导向矢量。对第 $k$ 个近场源的距离空间谱函数 ${{\boldsymbol{G}}_{{\rm{NF}}}}({\hat \phi _k},{\hat \theta _k},{r_k})$ 进行谱峰搜索，峰值位置即为第 $k$ 个近场源相对于阵列中心的距离估计值 ${\hat r_k}$ ，则第 $k$ 个近场源的定位结果为 $({\hat \phi _k},{\hat \theta _k},{\hat r_k})$ , $k =$ $1,2, \cdots ,{K_2}$ 。

4. 算法分析

4.1 计算复杂度分析

对于均匀圆阵下的近场源定位，传统的3-D MUSIC方法通过三维谱峰搜索估计出近场源方位角、俯仰角和距离，本文首先通过改进的ESPRIT-like方法得到近场源的二维DOA空间谱，利用二维谱峰搜索估计出近场源的方位角和俯仰角，进而利用1-D MUSIC方法得到近场源的距离空间谱，利用一维谱峰搜索估计出近场源的距离。本文提出的方法需要 ${M^2}N$ 次乘法和 ${M^2}(N - 1)$ 次加法计算混合源协方差矩阵，需要 $2{M^2}N$ 次乘法和 ${M^2}$ 次减法计算近场源差分矩阵，需要 $O({M^2})$ 量级的计算复杂度分别对混合源协方差矩阵以及近场源差分矩阵进行特征值分解，需要 $O({M^2})$ 量级的计算复杂度得到远场源二维DOA空间谱，需要 $4{K_2}{M^2} + 8K_2^2M$ 次乘法、 $(4{K_2}M + 8K_2^2)(M - 1)$ 次加法和 $4K_2^2$ 次减法得到近场源二维DOA空间谱，需要 $O(M)$ 量级的计算复杂度得到近场源距离空间谱。

4.2 混合源能量差异敏感性分析

由于本文所提的两步法通过协方差矩阵差分去除混合源的远场源和噪声分量，进而提取出近场源差分矩阵，因此对近场源进行估计时不受远场源和噪声能量的影响。本文第1步在对远场源进行估计时，当远场源的能量远远小于近场源的能量，此时通过第1步的方法对远场源进行定位的性能急剧恶化。为了解决两个源能量差异敏感的问题，可以首先对噪声能量进行估计，去除混合源协方差矩阵中的噪声分量；利用第2步的近场源估计结果，通过斜投影算法去除无噪混合源协方差矩阵中近场源分量，进而得到远场源分量。具体步骤如下：

混合源的协方差矩阵可以分解为

${\boldsymbol{R}} = \mathop {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{R}}_{11}}}& {{{\boldsymbol{R}}_{12}}}\\ {{{\boldsymbol{R}}_{21}}}& {{{\boldsymbol{R}}_{22}}} \end{array}} \right]}\limits^{K\qquad M - K} {\rm{ }}_{\scriptstyle\hfill\atop \scriptstyle M - K\hfill}^K$

(32)

噪声的能量可以通过式(33)计算

$\hat {\boldsymbol{\sigma}} _{\rm{n}}^2=\frac{{{\rm{tr}}\left\{ {{{\boldsymbol{R}}_{{\rm{22}}}}{\boldsymbol Z}} \right\}}}{{{\rm{tr}}\left\{ {\boldsymbol Z} \right\}}}$

(33)

其中， ${\boldsymbol Z}={{\boldsymbol{I}}_{M - K}}-{{\boldsymbol{R}}_{{\rm{21}}}}{\boldsymbol{R}}_{{\rm{21}}}^\dagger$ , ${\boldsymbol{R}}_{{\rm{21}}}^\dagger ={\left( {{\boldsymbol{R}}_{{\rm{21}}}^{\rm{H}}{{\boldsymbol{R}}_{{\rm{21}}}}} \right)^{ - 1}} {\boldsymbol{R}}_{{\rm{21}}}^{\rm{H}}$ ，(·)^†表示伪逆运算，tr(·)表示矩阵的迹， ${{\boldsymbol{I}}_{M - K}}$ 表示(M–K)维单位矩阵。因此，无噪的混合源协方差矩阵可以表示为

$\bar {\boldsymbol{R}}={{\boldsymbol{R}}_{{\rm{NF}}}} + {{\boldsymbol{R}}_{{\rm{FF}}}} = {\boldsymbol{R}} - \hat {\boldsymbol\sigma} _{\rm{n}}^{\rm{2}}{{\boldsymbol{I}}_M}$

(34)

由于斜投影矩阵 ${{\boldsymbol{E}}_{{{\boldsymbol{A}}_{{\rm{NF}}}}\left| {{{\boldsymbol{A}}_{{\rm{FF}}}}} \right.}}$ 满足

${{\boldsymbol{E}}_{{{\boldsymbol{A}}_{{\rm{NF}}}}\left| {{{\boldsymbol{A}}_{{\rm{FF}}}}} \right.}}{{\boldsymbol{A}}_{{\rm{NF}}}}={{\boldsymbol{A}}_{{\rm{NF}}}},\;{{\boldsymbol{E}}_{{{\boldsymbol{A}}_{{\rm{NF}}}}\left| {{{\boldsymbol{A}}_{{\rm{FF}}}}} \right.}}{{\boldsymbol{A}}_{{\rm{FF}}}}={{\boldsymbol{O}}_{M \times (K - {K_1})}}$

(35)

并可以通过式(36)计算

${\hat {\boldsymbol{E}}_{{{\boldsymbol{A}}_{{\rm{NF}}}}\left| {{{\boldsymbol{A}}_{{\rm{FF}}}}} \right.}}={\hat {\boldsymbol{A}}_{{\rm{NF}}}}{\left( {\hat {\boldsymbol{A}}_{{\rm{NF}}}^{\bf{H}}{{\bar {\boldsymbol{R}}}^\dagger }{{\hat {\boldsymbol{A}}}_{{\rm{NF}}}}} \right)^{ - 1}}\hat {\boldsymbol{A}}_{{\rm{NF}}}^{\bf{H}}{\bar {\boldsymbol{R}}^\dagger }$

(36)

其中， ${{\boldsymbol{O}}_{M \times (K - {K_1})}}$ 表示 $M \times (K - {K_1})$ 维零矩阵， ${\hat {\boldsymbol{A}}_{{\rm{NF}}}}$ 表示将近场源参数估计结果代入后的导向矩阵， ${\bar {\boldsymbol{R}}^\dagger }={\left( {{{\bar {\boldsymbol{R}}}^{\rm{H}}}\bar {\boldsymbol{R}}} \right)^{ - 1}}{\bar {\boldsymbol{R}}^{\rm{H}}}$ 。因此，远场源的接收矩阵可以表示为

$\begin{split} {{\boldsymbol{X}}_{{\rm{FF}}}}={{\boldsymbol{A}}_{{\rm{FF}}}}{{\boldsymbol{S}}_{{\rm{FF}}}}=&\left( {{{\boldsymbol{I}}_M}-{{\hat {\boldsymbol{E}}}_{{{\boldsymbol{A}}_{{\rm{NF}}}}\left| {{{\boldsymbol{A}}_{{\rm{FF}}}}} \right.}}} \right) \\ &\cdot \left( {{{\boldsymbol{A}}_{{\rm{NF}}}}{{\boldsymbol{S}}_{{\rm{NF}}}} + {{\boldsymbol{A}}_{{\rm{FF}}}}{{\boldsymbol{S}}_{{\rm{FF}}}}} \right)\end{split}$

(37)

进而远场源的协方差矩阵可以通过式(38)计算

$\begin{split} {\hat {\boldsymbol{R}}_{{\rm{FF}}}}={{\boldsymbol{X}}_{{\rm{FF}}}}{\boldsymbol{X}}_{{\rm{FF}}}^{\rm{H}}=&\left( {{{\boldsymbol{I}}_M}-{{\hat {\boldsymbol{E}}}_{{{\boldsymbol{A}}_{{\rm{NF}}}}\left| {{{\boldsymbol{A}}_{{\rm{FF}}}}} \right.}}} \right) \\ &\cdot \bar {\boldsymbol{R}}{\left( {{{\boldsymbol{I}}_M}-{{\hat {\boldsymbol{E}}}_{{{\boldsymbol{A}}_{{\rm{NF}}}}\left| {{{\boldsymbol{A}}_{{\rm{FF}}}}} \right.}}} \right)^{\rm{H}}}\end{split}$

(38)

此时，可以减少两个源的能量差异对远场源二维DOA估计的影响。

4.3 算法总结

基于协方差矩阵差分的混合源定位方法流程图如图2所示，对于远场源定位，本文利用均匀圆阵下的混合源接收数据得到混合源协方差矩阵，利用远场源导向矢量以及经过特征值分解得到的混合源噪声子空间，通过2-D MUSIC方法估计出远场源的方位角和俯仰角；对于近场源定位，本文利用混合源协方差矩阵和旋转矩阵得到近场源差分矩阵，利用随机满秩矩阵、对角矩阵以及经过特征值分解的近场源信号子空间，通过改进的ESPRIT-like方法计算出近场源的方位角和俯仰角，利用近场源的导向矢量和混合源的噪声子空间，通过1-D MUSIC方法估计出近场源的距离。本文根据均匀圆阵下的远场源在中心对称阵元的导向参数具有共轭性质，利用协方差矩阵差分方法去除混合源的远场源分量，得到近场源差分矩阵。由于协方差矩阵差分方法没有利用第1步的远场源二维DOA参数估计结果，因此，第2步近场源参数估计性能与远场源DOA估计性能无关。或者说，第1步和第2步的顺序也可以进行互换，可以先通过混合源协方差矩阵差分方法得到近场源差分矩阵，估计出近场源位置参数，再利用2-D MUSIC方法估计远场源DOA参数。另外，在多个混合源能量相同的场景下，由于第1步利用远场源导向矢量与噪声子空间的正交性，因此不会将近场源判断为远场源。当近场源的距离临近 $2{D^2}/\lambda$ ，由于噪声和模型误差的影响，第1步可能会将近场源的方向判断为远场源。当远场源的二维DOA与近场源二维DOA相同时，如果采用近场和远场参数联合估计将无法对混合源进行识别，采用本文所提的两步法可以将远场源和近场源进行有效分离，在混合源识别上更具优势。

图 2 基于协方差矩阵差分的混合源定位方法流程图

Figure 2. Flow chart of the mixed source localization based on the covariance differencing method

下载: 全尺寸图片幻灯片

5. 仿真实验

本节通过MATLAB仿真验证提出方法对混合源的定位效果，仿真中均匀圆阵的阵元个数为8，半径为0.5 m，以均匀圆阵的中心为原点建立三维坐标系，混合源包含一个远场源和一个近场源，远场源的位置为(50°, 30°, $\infty$ )，近场源的位置为(55°, 35°, 10 m)。

5.1 混合源定位有效性

这一小节用于验证本文提出方法得到的空间谱对远场和近场混合源的定位效果，在信噪比为20 dB，快拍数为600时，图3为通过本文提出方法得到的混合源定位的空间谱，图3(a)为通过2-D MUSIC方法得到的远场源二维DOA空间谱，方位角估计范围为0.1°～360°，间隔为0.1°，俯仰角估计范围为0°～90°，间隔0.1°，可以看出二维空间谱中出现的一个峰值对应的位置为(50.0°, 30.2°)，即为远场源方位角和俯仰角估计。图3(b)为通过改进的ESPRIT-like方法得到的近场源二维DOA空间谱，同样地，方位角估计范围为0.1°～360°，间隔为0.1°，俯仰角估计范围为0°～90°，间隔0.1°，可以看出二维空间谱中出现的一个峰值对应的位置为(55.0°, 35.1°)，即为近场源方位角和俯仰角估计。图3(c)为通过1-D MUSIC方法得到的近场源一维距离空间谱，近场源距离估计范围为0.1～30.0 m，间隔为0.1 m，图中红色的线表示近场源真实距离，可以看出距离空间谱中出现的一个峰值对应的位置为9.9 m，即为近场源距离估计。综上可知，本文所提协方差矩阵差分方法能够从混合源中提取出近场源分量，并通过改进的ESPRIT-like方法近场源三维位置参数中分离出二维DOA，最终实现混合源的有效定位。

图 3 混合源空间谱

Figure 3. Spatial spectrum of mixed sources

下载: 全尺寸图片幻灯片

5.2 不同信噪比下混合源定位精度

本节验证所提方法在不同信噪比下对混合源的定位性能，并与TSMUSIC方法^[16]、相位差方法^[17]以及克拉美罗界(Cramer-Rao Lower Bound, CRLB)进行对比，其中TSMUSIC方法是利用2-D MUSIC方法分别估计出远场源和近场源的二维DOA，利用1-D MUSIC方法分别估计出近场源的距离；相位差方法是利用均匀圆阵下对角阵元的相位差，通过最小二乘法同时估计出远场源和近场源的二维DOA，根据1-D MUSIC方法对混合源进行识别并估计出近场源的距离参数；CRLB是参数任意无偏估计量的下限。在仿真中利用均方根误差(Root Mean Square Error，RMSE)对算法的性能进行分析，RMSE可由式(39)计算

${\rm{RMSE}}\left( \beta \right) = \sqrt {\frac{1}{L}\sum\limits_{l = 1}^L {{{\left( {{{\hat \beta }_l} - \beta } \right)}^2}} }$

(39)

其中，L表示独立重复仿真次数， ${\hat \beta _l}$ 表示第 $l$ 次的混合源参数估计结果， $l = 1,2, ··· ,L$ , $\beta$ 表示混合源的真实位置参数。当快拍数为600时，图4(a)—图4(e)分别为通过600次独立重复仿真得到在不同信噪比下的远场源方位角RMSE、远场源俯仰角RMSE、近场源方位角RMSE、近场源俯仰角RMSE以及近场源距离RMSE，从图4可以看出混合源的RMSE随着信噪比的增大而降低，表明混合源参数估计精度不断提高。另外，还可以看出本文所提基于协方差矩阵差分方法的混合源定位精度比基于相位差方法的定位精度高。由于本文方法与TSMUSIC方法都采用2-D MUSIC方法对远场源进行定位，因此两种方法的远场源定位精度相当，本文所提改进ESPRIT-like方法对近场源二维DOA参数估计精度略低于TSMUSIC方法参数估计精度。此外，由于近场源距离估计是在近场源二维DOA估计基础上进行的，相较于近场源二维DOA估计误差，近场源距离估计误差相对较大，从图4中也可以看出，在信噪比为20 dB时，近场源的二维DOA估计误差接近0.01°，而近场源的距离误差接近于0.1 m。

图 4 不同信噪比下混合源RMSE

Figure 4. RMSE of mixed sources versus SNR

下载: 全尺寸图片幻灯片

5.3 不同快拍数下混合源定位精度

本节进一步验证在不同快拍数下混合源的定位性能，本文所提方法同样与TSMUSIC方法、相位差方法以及克拉美罗界进行了对比。当信噪比为20 dB时，图5(a)—图5(e)分别为通过600次独立重复仿真得到在不同快拍数下的远场源方位角RMSE、远场源俯仰角RMSE、近场源方位角RMSE、近场源俯仰角RMSE以及近场源距离RMSE，从图5中可以看出混合源RMSE随着快拍数的增大而降低，表明混合源参数估计精度逐渐提高。另外，从图中可以看出当快拍数小于600时，混合源RMSE快速下降，当快拍数超过600时，混合源RMSE逐渐平缓。此外，从图5可以看出本文提出方法的混合源定位精度比基于相位差方法的定位精度高，对于远场源参数估计精度，本文提出的方法与TSMUSIC方法相当，对于近场源参数估计精度，本文提出的方法略低于TSMUSIC方法。

图 5 不同快拍数下的混合源RMSE

Figure 5. RMSE of mixed sources versus snapshots

下载: 全尺寸图片幻灯片

5.4 混合源能量差异性对比

由于本文所提的两步法通过协方差矩阵差分去除混合源的远场源和噪声分量，进而提取出近场源差分矩阵，因此对近场源进行估计时不受远场源和噪声能量的影响。图6(a)—图6(c)给出了当近场源和远场源能量比(Near-field and Far-field source Power Ratio, NFPR)分别为1.0, 0.2和0.1时的近场源空间谱，可以看出近场源空间谱出现一个峰值，峰值所对应的位置即为近场源二维DOA。仿真结果表明，通过本文提出的算法对近场源进行估计时不受远场源和噪声能量的影响。

图 6 近场源空间谱

Figure 6. Spatial spectrum of near-field source

下载: 全尺寸图片幻灯片

为了解决远场源对能量差异敏感的问题，通过对信噪比和近场源进行估计，去除混合源协方差矩阵中的噪声分量和近场源分量，得到式(38)远场源协方差矩阵，进而估计远场源二维DOA。图7(a)—图7(c)给出了当近场源和远场源能量比分别为1, 5和10时的近场源空间谱，可以看出远场源空间谱出现了一个峰值，峰值所对应的位置即为远场源二维DOA。仿真结果表明，通过本文提出的算法对近场源进行估计时不受远场源和噪声能量的影响。仿真结果表明，所提算法能够有效地解决远场源能量差异敏感的问题。

图 7 远场源空间谱

Figure 7. Spatial spectrum of far-field source

下载: 全尺寸图片幻灯片

6. 结论

本文所提均匀圆阵下利用混合源协方差矩阵差分的方法能够消除混合源中的近场源和噪声分量，提取出近场源差分矩阵，实现远场源和近场源的有效分离，利用2-D MUSIC方法估计出远场源的方位角和俯仰角，通过改进的ESPRIT-like方法提取并估计出近场源的方位角和俯仰角，运用1-D MUSIC方法估计出近场源的距离参数。仿真结果表明，相较于基于相位差的混合源定位算法，本文所提方法能够提高混合源的定位精度。由于协方差矩阵差分方法没有利用第1步的远场源二维DOA参数估计结果，因此，第2步近场源参数估计性能与远场源DOA估计性能无关。或者说，第1步和第2步的顺序也可以进行互换，可以先通过混合源协方差矩阵差分方法得到近场源差分矩阵，估计出近场源位置参数，再利用2-D MUSIC方法估计远场源DOA参数。另外，在多个混合源能量相同的场景下，由于第1步利用远场源导向矢量与噪声子空间的正交性，因此不会将近场源判断为远场源。当近场源的距离临近 $2{D^2}/\lambda$ ，由于噪声和模型误差的影响，第1步可能会将近场源的方向判断为远场源。当远场源的二维DOA与近场源二维DOA相同时，如果采用近场和远场参数联合估计将无法对混合源进行识别，采用本文所提的两步法可以将远场源和近场源进行有效分离，在混合源识别上更具优势。

图 1 仿真的飞机尾流雷达回波

Figure 1. Simulation of the aircraft wake radar echo

下载: 全尺寸图片幻灯片

图 2 马可尼研究中心X波段雷达测量的功率谱^[16]

Figure 2. The Power Spectrum measured by the X-band radar in Marconi Research Center

下载: 全尺寸图片幻灯片

图 3 基于信息几何的矩阵CFAR检测器框图

Figure 3. Block diagram of CA-CFAR detector based on information geometry

下载: 全尺寸图片幻灯片

图 4 脉冲个数为16时的自相关函数与功率谱

Figure 4. Autocorrelation function and power spectrum when the pulse number is 16

下载: 全尺寸图片幻灯片

图 5 协方差矩阵降低维数之前与之后的功率谱

Figure 5. The power spectrum before and after the covariance matrix reduces the dimension

下载: 全尺寸图片幻灯片

图 6 随SNR的变化测地线距离的变化

Figure 6. The variation of geodesic distance with the variation of SNR

下载: 全尺寸图片幻灯片

图 7 矩阵均值迭代估计性能

Figure 7. Iterative estimation performance of matrix mean

下载: 全尺寸图片幻灯片

图 8 有无尾流情况下检测统计量的统计直方图

Figure 8. The statistical histogram of detection statistics with and without the vortex target

下载: 全尺寸图片幻灯片

图 9 基于矩阵CFAR的尾流目标的检测概率

Figure 9. Detection probability of the vortex target based on matrix CFAR

下载: 全尺寸图片幻灯片

参考文献(18)

[1]	孙华飞, 张真宁, 彭林玉, 等. 信息几何导引[M]. 北京: 科学出版社, 2016. Sun Hua-fei, Zhang Zhen-ning, Peng Lin-yu, et al.. An Elementary Introduction to Information Geometry[M]. Beijing: Science Press, 2016.
[2]	Rao C. Information and the accuracy attainable in the estimation of statistical parameters[J]. Bulletin of Calcutta Mathematical Society, 1945, 37: 81–91.
[3]	Chencov N N. Statistical Decision Rules and Optimal Inference[M]. Rhode Island, USA: American Mathematical Society, 1982.
[4]	Amari S I and Nagaoka H. Methods of Information Geometry[M]. Providence, RI: American Mathematical Society, 2000.
[5]	Lenglet C, Rousson M, Deriche R, et al. Statistics on the manifold of multivariate normal distributions: Theory and application to diffusion tensor MRI processing[J]. Journal of Mathematical Imaging and Vision, 2006, 25(3): 423–444. DOI: 10.1007/s10851-006-6897-z
[6]	Moakher M. A differential geometric approach to the geometric mean of symmetric positive-definite matrices[J]. SIAM Journal on Matrix Analysis and Applications, 2005, 26(3): 735–747. DOI: 10.1137/S0895479803436937
[7]	Nielsen F and Bhatia R. Matrix Information Geometry[M]. Berlin Heidelberg: Springer, 2013.
[8]	Barbaresco F. Interactions between symmetric cone and information geometries: Bruhat-Tits and Siegel spaces models for high resolution autoregressive Doppler imagery[C]. Proceedings of ETVC 2008 Conference, Berlin, Heidelberg, 2009: 124–163.
[9]	Barbaresco F. Innovative tools for radar signal processing based on Cartan’s geometry of SPD matrices & information geometry[C]. Proceedings of 2008 IEEE Radar Conference, Rome, 2008: 1–6.
[10]	Pennec X. Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[J]. Journal of Mathematical Imaging and Vision, 2006, 25(1): 127–154. DOI: 10.1007/s10851-006-6228-4
[11]	Barbaresco F. Robust statistical radar processing in Fréchet metric space: OS-HDR-CFAR and OS-STAP processing in Siegel homogeneous bounded domains[C]. Proceedings of 2011 IEEE International Radar Symposium, Leipzig, Germany, 2011: 639–644.
[12]	黎湘, 程永强, 王宏强, 等. 雷达信号处理的信息几何方法[M]. 北京: 科学出版社, 2014. Li Xiang, Cheng Yong-qiang, Wang Hong-qiang, et al.. Methods of Information Geometry of Radar Signal Processing[M]. Beijing: Science Press, 2014.
[13]	Broderick A J, Bevilaqua P, Crouch J, et al.. Wake Turbulence: An Obstacle to Increased Air Traffic Capacity[M]. Washington, DC: The National Academies Press, 2008: 1–59.
[14]	刘俊凯, 李文臣, 王雪松, 等. 基于多普勒特性的飞机尾流回波提取方法[J]. 系统仿真学报, 2011, 23(7): 1323–1328. DOI: 10.16182/j.cnki.joss.2011.07.015 Liu Jun-kai, Li Wen-chen, Wang Xue-song, et al. Extraction of aircraft wake vortices radar returns based on the Doppler characteristics[J]. Journal of System Simulation, 2011, 23(7): 1323–1328. DOI: 10.16182/j.cnki.joss.2011.07.015
[15]	Gerz T, Holzäpfel F, and Darracq D. Commercial aircraft wake vortices[J]. Progress in Aerospace Sciences, 2002, 38(3): 181–208. DOI: 10.1016/S0376-0421(02)00004-0
[16]	Shephard D J, Kyte A P, and Segura C A. Radar wake vortex measurements at F and I band[C]. Proceedings of IEE Colloquium on Radar and Microwave Imaging, London, 1994: 7/1–7/5.
[17]	王首勇, 万洋, 刘俊凯, 等, 著. 现代雷达目标检测理论与方法[M]. 第2版, 北京: 科学出版社, 2015. Wang Shou-yong, Wan Yang, Liu Jun-kai, et al.. Modern Radar Target Detection Theory and Methods[M]. Second Edition, Beijing: Science Press, 2015.
[18]	Calvo M and Oller J M. An explicit solution of information geodesic equations for the multivariate normal model[J]. Statistics&Risk Modeling, 1991, 9(1/2): 119–138.

施引文献

期刊类型引用(5)

1.	马助兴，张立硕，徐红元，郑焕坤，赵智龙，康哲. 面向变电站智能巡检5G大规模MIMO融合定位方法. 全球能源互联网. 2023(03): 308-315 . 百度学术
2.	任笑莹，王英民，张立琛，王奇，廉杰. 均匀双圆环阵虚拟成阵方法研究. 水下无人系统学报. 2023(05): 735-745 . 百度学术
3.	苏晓龙，户盼鹤，刘天鹏，彭勃，程耘，刘振. 基于深度展开ISTA网络的混合源定位方法. 信号处理. 2022(10): 2082-2091 . 百度学术
4.	母采凤，李森，吕梦然. 脉冲噪声环境下基于矩阵差分的远近场混合源定位. 信号处理. 2022(11): 2342-2349 . 百度学术
5.	张国鑫，易伟，孔令讲. 基于1比特量化的大规模MIMO雷达系统直接定位算法. 雷达学报. 2021(06): 970-981 . 本站查看