doi:10.12000/JR17070

Circulation Retrieval of Wake Vortices under Rainy Conditions with an X Band Radar

DOI: 10.12000/JR17070

BU Advanced Radar Concept, Surface Radar, Thales Air Systems, Limours, France

More Information

Author Bio:
Jean-Yves Schneider, Engineer specialist in Advanced Studies

Dr Gilles Beauquet, Engineer specialist in Electromagnetic modelling

Frédéric Barbaresco, Thales senior expert in the processing & computing domain for Thales Land & Air Systems

Corresponding author: Jean-Yves Schneider. E-mail: jean-yves-j.schneider@thalesgroup.com

摘要

Abstract:
At airports, runway operation is the limiting factor for the overall throughput; specifically the fixed and overly conservative ICAO wake turbulence separation minima. The wake turbulence hazardous flows can dissipate quicker because of decay due to air turbulence or be transported out of the way on oncoming traffic by cross-wind, yet wake turbulence separation minima do not take into account wind conditions. Indeed, for safety reasons, most airports assume a worst-case scenario and use conservative separations; the interval between aircraft taking off or landing therefore often amounts to several minutes. However, with the aid of accurate wind data and precise measurements of wake vortex by radar sensors, more efficient intervals can be set, particularly when weather conditions are stable. Depending on traffic volume, these adjustments can generate capacity gains, which have major commercial benefits. This paper presents the use of Electronic scanning radar for detecting wake vortices. In this method, the raindrops Doppler spectrogram is used to retrieve the strength of the wake vortex. Numerical simulation are performed to establish an empirical model used during the retrieval method. This paper presents also the results obtained during the trials of the PARIS-CDG data set recorded from October 2014 to November 2015 with an X-band RADAR developed and deployed by THALES.
- Wake-vortex hazard /
- Airport capacity /
- Airport safety /
- X-band radar /
- Wake-vortex circulation /
- Eddy dissipation rate

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相同时，如果采用近场和远场参数联合估计将无法对混合源进行识别，采用本文所提的两步法可以将远场源和近场源进行有效分离，在混合源识别上更具优势。

Figure 1. THALES X-band RADAR

下载: 全尺寸图片幻灯片

Figure 2. Location of THALES X-band RADAR

下载: 全尺寸图片幻灯片

Figure 3. RADAR scanning domain (horizontal)

下载: 全尺寸图片幻灯片

Figure 4. RADAR scanning domain (vertical)

下载: 全尺寸图片幻灯片

Figure 5. Global architecture of RADAR system deployed in PARIS-CDG

下载: 全尺寸图片幻灯片

Figure 6. X-band RADAR modes available

下载: 全尺寸图片幻灯片

Figure 7. Notations used to define a scenario

下载: 全尺寸图片幻灯片

Figure 8. Mean Doppler velocity as a function of range-Case of A319-Beam elevation 5°

下载: 全尺寸图片幻灯片

Figure 9. Comparison between T_L(max) and its fitted function-All beams

下载: 全尺寸图片幻灯片

Figure 10. Comparison between T_R(min) and its fitted function- All beams

下载: 全尺寸图片幻灯片

Figure 11. Block diagram of wake vortex processing

下载: 全尺寸图片幻灯片

Figure 12. Range evaluation-A380-RADAR beam elevation 5°

下载: 全尺寸图片幻灯片

Figure 13. Elevation evaluation-A380-RADAR beam elevation 5°

下载: 全尺寸图片幻灯片

Figure 14. Circulation evaluation-A380-RADAR beam elevation 5°

下载: 全尺寸图片幻灯片

Figure 15. Range/Doppler map-Impact of ground clutter

下载: 全尺寸图片幻灯片

Figure 16. Pre-processing-A342-2015/01/28 05:22-All beams-scan 2

下载: 全尺寸图片幻灯片

Figure 17. V_min & V_max-A342-2015/01/28 05:22-All beams

下载: 全尺寸图片幻灯片

Figure 18. D_L & D_R-A342-2015/01/28 05:22-All beams

下载: 全尺寸图片幻灯片

Figure 19. Data availability-Threshold 27R landing recordings

下载: 全尺寸图片幻灯片

Figure 20. Histogram of instantaneous wind direction (MET data)

下载: 全尺寸图片幻灯片

Figure 21. Histogram of instantaneous wind speed (mean=5.4 m/s, MET data)

下载: 全尺寸图片幻灯片

Figure 22. Durations of vortices for all landing aircraft per category-Threshold 27R

下载: 全尺寸图片幻灯片

Figure 23. Coefficient a versus b₀ (all categories)-Threshold 27R

下载: 全尺寸图片幻灯片

Figure 24. Histogram of measured V_min & V_max (all categories)-Threshold 27R

下载: 全尺寸图片幻灯片

Table 1. Statistics about observed landing aircrafts

Number of recordings (Any SNR)	category M:	category H:	category S:		Number of recordings in presence of rain (SNRe≥5 dB)	category M:	category H:	category S:
26924	18411	7380	1131		7936	5403	2193	340

下载: 导出CSV

Table 2. Calculation of detection probability 1/2

WTC (Wake Turbulence Category)	SNR > 5 dB	SNR > 10 dB
WTC (Wake Turbulence Category)	detected/total	detected/total
Cat. M	1734/2332 (74.4%)	1362/1803 (75.5%)
Cat. H	740/820 (90.2%)	573/632 (90.7%)
Cat. S	154/169 (91.1%)	103/113 (91.2%)

下载: 导出CSV

Table 3. Calculation of detection probability 2/2

WTC	SNR > 15 dB	SNR > 20 dB
WTC	detected/total	detected/total
Cat. M	1001/1282 (78.1%)	472/632 (74.7%)
Cat. H	394/436 (90.4%)	211/226 (93.4%)
Cat. S	64/70 (91.4%)	27/30 (90.0%)

下载: 导出CSV

Table 4. Statistical results

Category	N_points	$\nu$ _avg (m/s)	$\sigma$ _avg (m/s)	$\nu$ _mode (m/s)
All	2639	2,10	0,89	1,75
M	1879	1,94	0,79	1,75
H	638	2,49	0,95	1,75
S	2639	2,10	0,89	1,75

下载: 导出CSV

参考文献(28)

[1]	12.2.2, Preliminary System Requirements of Runway Wake Vortex Detection, Prediction and decision support tools, D1,00.01.00, 11/07/2010.
[2]	12.2.2, Preliminary System Architecture of Runway Wake Vortex Detection, Prediction and Decision Support Tools, D2,00.01.00, 27/10/2010.
[3]	12.2.2, Recommendations on sensor technologies, D4,00.01.01, 14/10/2011.
[4]	Shariff K and Wray A. Analysis of the Radar reflectivity of aircraft vortex wakes[J]. Journal of Fluid Mechanics, 2002, 463(1): 121–161.
[5]	Liu Z X. Modélisation des signatures radar des tourbillons de sillage par temps de pluie[D]. [Ph.D. dissertation], 2013.
[6]	Li J L, Wang X S, Wang T, et al. Circulation retrieval of wake vortex under rainy condition with a vertically pointing radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2017, 53(4): 1893–1906. DOI: 10.1109/TAES.2017.2675198
[7]	Barbaresco F, Brion V, and Jeannin N. Radar wake-vortices cross-section/Doppler signature characterisation based on simulation and field tests trials[J]. IET Radar,Sonar&Navigation, 2016, 10(1): 82–96.
[8]	Oude Nijhuis A C P, Thobois L P, Barbaresco F, et al.. Wind hazard and turbulence monitoring at airports with lidar and radar sensors and Mode-S downlinks: The UFO Project[J]. Bulletin of the American Meteorological Society. (in press)
[9]	Kovalev D, Vanhoenacker-Janvier D, Wilson R, et al.. Electromagnetic wind radar simulator validation using meteorological data and a zenith X-band radar[C]. 2016 European Radar Conference (EuRAD 2016), London, 2016.
[10]	Barbaresco F, Thobois L, Dolfi-Bouteyre A, et al.. Monitoring wind, turbulence and aircraft wake vortices by high resolution RADAR and LIDAR remote sensors in all weather conditions[Z]. Journées scientifiques URSI, 2015.
[11]	Barbaresco F, Jeantet A, and Meier U. Wake vortex detection & monitoring by X-band doppler radar: Paris orly radar campaign results[C]. Proceedings of IET International Conference on Radar Systems, Edinburgh, UK, 2007.
[12]	Barbaresco F, Wasselin J P, Jeantet A, et al.. Wake vortex monitoring & profiling by Doppler X-band radar in all weather conditions[C]. Eurocontrol Innovative Research Workshop, Bretigny, 2007.
[13]	Barbaresco F, Wasselin J P, Jeantet A, et al.. Wake vortex profiling by Doppler X-band radar: Orly trials at initial take-off & ILS interception critical areas[C]. Proceedings of IEEE Radar Conference, Rome, 2008.
[14]	Barbaresco F, Jeantet A, and Meier U. Wake vortex X-band radar monitoring: Paris-CDG airport 2008 campaign results & prospectives[C]. Proceedings of International Radar Conference-Surveillance for a Safer World, Bordeaux, France, 2009.
[15]	Barbaresco F and Meier U. Radar monitoring of a wake vortex: Electromagnetic reflection of wake turbulence in clear air[J]. Comptes Rendus Physique, 2010, 11(1): 54–67. DOI: 10.1016/j.crhy.2010.01.001
[16]	Barbaresco F. Airport radar monitoring of wake vortex in all weather conditions[C]. Proceedings of 2010 European Radar Conference (EuRAD), Paris, 2010: 85–88.
[17]	Barbaresco F, Juge P, Klein M, et al.. Optimising runway throughput through wake vortex detection, prediction and decision support tools[C]. Proceedings of 2011 Tyrrhenian International Workshop on Digital Communications-Enhanced Surveillance of Aircraft and Vehicles (TIWDC/ESAV), Capri, Italy, 2011.
[18]	Barbaresco F, Juge P, Klein M, et al.. Wake vortex detection, prediction and decision support tools[C]. Proceedings of the 2013 IEEE/AIAA 32nd Digital Avionics Systems Conference (DASC), East Syracuse, NY, 2013: 6B1-1–6B1-15.
[19]	Barbaresco F, Juge P, Klein M, et al.. Boom of airport capacity based on wake vortex hasards mitigation sensors and systems[C]. AUN2014: Airports in Urban Networks, Paris, 2014.
[20]	Barbaresco F, Juge P, Bruchec P, et al.. Eddy Dissipation Rate (EDR) retrieval with ultra-fast high range resolution electronic-scanning X-band airport radar: Results of European FP7 UFO toulouse airport trials[C]. Proceedings of 2015 European Radar Conference (EuRAD), Paris, 2015.
[21]	Speiker L J P, Barbaresco F, Frech M, et al. ATC-Wake: Integrated wake vortex safety & capacity system[J]. Journal of Air Traffic Control, 2007, 49(1): 17–32.
[22]	Winckelmans G, Desenfans O, Barbaresco F, et al.. The ATC-Wake Predictor system and its potential use to increase the capacity at airports[C]. International Symposium on Sensors and Systems for Airport Surveillance, Paris, France, 2005.
[23]	Barbaresco F. Compréhension et maîtrise des tourbillons de sillage: Cinq siècles d’aventures de Léonard de Vinci à Jean Leray[Z]. Revue REE, 2013(3): 84–88.
[24]	Barbaresco F. Airport radar monitoring of wake vortex in all weather conditions[C]. Proceedings of 2010 European Radar Conference (EuRAD), Paris, 2010.
[25]	Vanhoenacker-Janvier D, Djafri K, della Faille de Leverghem R, et al.. Simulation of the radar cross-section of wake vortices in clear air[C]. Proceedings of the Seventh European Conference on Radar in Meteorology and Hydrology, Toulouse, 2012.
[26]	Barbaresco F, Juge P, Klein M, et al.. Wake vortex detection, prediction and decision support tools in SESAR program[C]. Proceedings of the 32th Digital Avionics Systems Conference (DASC), East Syracuse, NY, 2013.
[27]	Barbaresco F. Radar monitoring of wake vortex: Electromagnetic reflection of wake turbulence in clear air[Z]. CR Physique Académie des Sciences, Elsevier, 2010.
[28]	Liu Z, Jeannin N, Vincent F, et al.. Development of a radar simulator for monitoring wake vortices in rainy weather[C]. Proceedings of 2011 IEEE CIE International Conference on Radar (Radar), Chengdu, China, 2011.

施引文献

期刊类型引用(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 . 本站查看