基于马尔科夫判别谱聚类的极化SAR图像分类方法

张向荣; 于心源; 唐旭; 侯彪; 焦李成

doi:10.12000/JR19059

基于马尔科夫判别谱聚类的极化SAR图像分类方法

DOI: 10.12000/JR19059

西安电子科技大学人工智能学院，智能感知与图像理解教育部重点实验室，国际智能感知与计算联合研究中心西安 710071

基金项目: 国家自然科学基金(61772400)，陕西省重点研发计划(2019ZDLGY03-08)

详细信息

作者简介:
张向荣(1978–)，女，西安电子科技大学人工智能学院，教授/博导，IEEE高级会员，IEEE GRSS西安分会副主席。主要从事遥感影像智能解译、机器学习、模式识别相关方向研究。E-mail: xrzhang@mail.xidian.edu.cn

于心源(1993–)，女，西安电子科技大学研究生。研究领域为极化SAR图像处理、机器学习等

唐　旭(1985–)，男，西安电子科技大学人工智能学院讲师，IEEE会员。主要从事遥感影像内容解译方向的研究

侯　彪(1974–)，男，西安电子科技大学人工智能学院，教授/博导，IEEE会员。从事人工智能、类脑计算、遥感脑、图像和视频分析、智能教育等研究

焦李成(1959–)，男，西安电子科技大学人工智能学院，教授/博导，IEEE Fellow, IEEE GRSS西安分会主席，IEEE Transactions on Geoscience and Remote Sensing副主编。主要从事人工智能相关领域的研究

通讯作者:
张向荣　xrzhang@mail.xidian.edu.cn

中图分类号: TN958
计量
- 文章访问数: 2586
- HTML全文浏览量: 868
- PDF下载量: 250
- 被引次数: 11
出版历程
- 收稿日期: 2019-06-01
- 修回日期: 2019-07-22
- 网络出版日期: 2019-08-28

PolSAR Image Classification Method Based on Markov Discriminant Spectral Clustering

Artificial Intelligence Institute, Xidian University, Key Laboratory of Intelligent Perception and Image Understanding of Ministry of Education, International Research Center for Intelligent Perception and Computation, Xi’an 710071, China

Funds: The National Natural Science Foundation of China (61772400), The Key Research and Development Plans of Shaanxi Province (2019ZDLGY03-08)

More Information

Corresponding author: ZHANG Xiangrong, xrzhang@mail.xidian.edu.cn

摘要

摘要: 该文针对现有的谱聚类方法用于极化SAR图像分类时精度较低的问题，提出一种基于马尔科夫的判别谱聚类方法(MDSC)，具有低秩和稀疏分解的特点。该方法首先恢复一个真实的低秩概率转移矩阵，将其作为标准马尔科夫谱聚类方法的输入，以减少噪声对分类结果的影响；然后在目标函数中引入判别信息，使极化SAR图像的数据信息能够得到更加充分地利用；最后采用增广拉格朗日乘子法来解决低秩和概率单纯形约束下的目标函数优化问题。在荷兰小农田、德国、西安和荷兰大农田4个不同数据集上的实验证明，该方法具有较好的准确率，且参数敏感性较低，表现出了良好的分类性能。
- 极化SAR /
- 谱聚类 /
- 判别谱聚类 /
- 多视角谱聚类
Abstract: Due to the existing spectral clustering methods have low accuracy for PolSAR image classification, a Markov-based Discriminative Spectral Clustering(MDSC) method is proposed, which has the characteristics of low-rank and sparse decomposition. Firstly, a real low-rank probability transfer matrix is restored as an input to the standard Markov spectral clustering method to reduce the influence of noise on the classification result. Then the discriminative information is introduced into the objective function to make the polarimetric SAR image data can be more fully used. Finally, the augmented Lagrangian multiplier method is used to solve the objective function optimization problem under low-rank and probability simplex constraints. Experiments on three different data sets of Flevoland, Oberpfaffenhofen, and Xi’an show that our method has good accuracy and low sensitivity, which having a good classification performance.
- PolSAR /
- Spectral clustering /
- Discriminant spectral clustering /
- Multi-view spectral clustering

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. Real probability transfer matrix construction profile

下载: 全尺寸图片幻灯片

图 2 本文算法框架图

Figure 2. Algorithm frame diagram

下载: 全尺寸图片幻灯片

图 3 荷兰Flevoland地区农田小图的伪彩图、类标图以及不同算法的分类结果图

Figure 3. Pseudo-color map, class diagram and classification results of different algorithms for farmland maps in the Flevoland region of the Netherlands

下载: 全尺寸图片幻灯片

图 4 德国Oberpfaffenhofen地区数据的伪彩图、类标图以及不同算法的分类结果图

Figure 4. Pseudo-color map, class diagram and data classification results of different algorithms in the Oberpfaffenhofen region of Germany

下载: 全尺寸图片幻灯片

图 5 西安地区数据的伪彩图、类标图以及不同算法的分类结果图

Figure 5. Pseudo-color map, class diagram and data classification results of different algorithms in Xi’an area

下载: 全尺寸图片幻灯片

图 6 荷兰 Flevoland 地区大农田数据的伪彩图、类标图以及不同算法的分类结果图

Figure 6. Pseudo-color map, class diagram and classification results of different algorithms for large farmland data in the Flevoland region of the Netherlands

下载: 全尺寸图片幻灯片

图 7 荷兰小农田中不同 $\lambda$ 和 $\beta$ 下的分类结果图

Figure 7. Classification results of different $\lambda$ and $\beta$ below in small Dutch farmland

下载: 全尺寸图片幻灯片

图 10 不同正则项参数 $\xi$ 的分类结果图

Figure 10. Classification results of different regular item parameter $\xi$

下载: 全尺寸图片幻灯片

图 8 德国地区中不同 $\lambda$ 和 $\beta$ 下的分类结果

Figure 8. Classification results for different $\lambda$ and $\beta$ below in the German region

下载: 全尺寸图片幻灯片

图 9 西安地区中不同 $\lambda$ 和 $\beta$ 下的分类结果图

Figure 9. Classification results of different $\lambda$ and $\beta$ subordinates in Xi’an area

下载: 全尺寸图片幻灯片

表 1 4种算法对Flevoland小农田图的分类结果

Table 1. Classification results of four algorithms for Flevoland small farmland map

	裸土	土豆	甜菜	大麦	豌豆	小麦	OA	AA	Kappa
Co-Reg	0.8860	0.9452	0.7551	0.7906	0.8262	0.8960	0.8400	0.8498	0.8775
MMC	0.9180	0.9580	0.7242	0.9623	0.8756	0.6994	0.8708	0.8396	0.9005
SR-MO	0.9034	0.9049	0.8845	0.9561	0.8362	0.9554	0.9130	0.9067	0.9331
本文算法	0.9048	0.9088	0.8834	0.9604	0.8920	0.9382	0.9243	0.9146	0.9418

下载: 导出CSV

表 2 4种算法对德国Oberpfaffenhofen地区的分类结果

Table 2. Classification results of four algorithms for the Oberpfaffenhofen region of Germany

	Co-Reg	MMC	SR-MO	本文算法
农田	0.6111	0.6018	0.6859	0.7016
居民区	0.6072	0.6521	0.7336	0.7389
林地	0.8162	0.7993	0.9055	0.9108
道路	0.5311	0.5681	0.6049	0.6418
其他	0.8791	0.8789	0.8673	0.8814
OA	0.7363	0.7471	0.7822	0.7974
AA	0.6889	0.7000	0.7594	0.7749
Kappa	0.6170	0.6348	0.6920	0.7205

下载: 导出CSV

表 3 4种算法对西安地区的分类结果

Table 3. Classification results of four algorithms for Xi’an area

	Co-Reg	MMC	SR-MO	本文算法
河流	0.9372	0.8876	0.9189	0.8890
城区	0.7300	0.6622	0.8128	0.8550
植被	0.6876	0.8052	0.8006	0.8555
OA	0.7400	0.7670	0.8227	0.8503
AA	0.7259	0.7584	0.8136	0.8436
Kappa	0.7341	0.7710	0.8071	0.8471

下载: 导出CSV

表 4 4种算法对荷兰 Flevoland 地区大农田图的分类结果

Table 4. Classification results of four algorithms for large farmland maps in the Flevoland region of the Netherlands

	Co-Reg	MMC	SR-MO	本文算法
蚕豆	0.7459	0.8942	0.9614	0.9584
油菜籽	0.1393	0.7194	0.7094	0.8337
裸地	0.2056	0.9616	0.9541	0.9583
土豆	0.2479	0.8912	0.8796	0.9086
甜菜	0.1079	0.9481	0.9656	0.9515
小麦2	0.2427	0.6251	0.8525	0.7941
豌豆	0.7932	0.9517	0.8887	0.9571
小麦3	0.5412	0.9231	0.9180	0.9300
苜蓿	0.9541	0.8940	0.8391	0.9284
大麦	0.9226	0.6311	0.9660	0.8524
小麦	0.1158	0.8458	0.8660	0.8796
草地	0.3944	0.6459	0.7470	0.8773
森林	0.4041	0.8833	0.8329	0.9122
水域	0.5403	0.9757	0.9035	0.9620
建筑物	0.5954	0.7761	0.5865	0.7912
OA	0.4204	0.8441	0.8501	0.8923
AA	0.4165	0.8398	0.8697	0.9043
Kappa	0.4269	0.8409	0.8441	0.9136

下载: 导出CSV

参考文献(13)

[1]	邹焕新, 罗天成, 张月, 等. 基于组合条件随机场的极化SAR图像监督地物分类[J]. 雷达学报, 2017, 6(5): 541–553. doi: 10.12000/JR16109 ZOU Huanxin, LUO Tiancheng, ZHANG Yue, et al. Combined conditional random fields model for supervised PolSAR images classification[J]. Journal of Radars, 2017, 6(5): 541–553. doi: 10.12000/JR16109
[2]	FREEMAN A and DURDEN S L. A three-component scattering model for polarimetric SAR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 1998, 36(3): 963–973. doi: 10.1109/36.673687
[3]	BUONO A, NUNZIATA F, MIGLIACCIO M, et al. Classification of the Yellow River delta area using fully polarimetric SAR measurements[J]. International Journal of Remote Sensing, 2017, 38(23): 6714–6734. doi: 10.1080/01431161.2017.1363437
[4]	JOULIN A, BACH F, and PONCE J. Discriminative clustering for image co-segmentation[C]. Proceedings of 2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, San Francisco, USA, 2010: 1943–1950.
[5]	HARALICK R M, SHANMUGAM K, and DINSTEIN I. Textural features for image classification[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1973, SMC-3(6): 610–621. doi: 10.1109/TSMC.1973.4309314
[6]	RATHA D, BHATTACHARYA A, and FRERY A C. Unsupervised classification of PolSAR data using a scattering similarity measure derived from a geodesic distance[J]. IEEE Geoscience and Remote Sensing Letters, 2018, 15(1): 151–155. doi: 10.1109/LGRS.2017.2778749
[7]	AHMED N and CAMPBELL M. Variational Bayesian learning of probabilistic discriminative models with latent softmax variables[J]. IEEE Transactions on Signal Processing, 2011, 59(7): 3143–3153. doi: 10.1109/TSP.2011.2144587
[8]	LIN Zhouchen, CHEN Minming, and MA Ya. The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices[J]. Eprint Arxiv, 2010(9): 26.
[9]	XIA Rongkai, PAN Yan, DU Lei, et al. Robust multi-view spectral clustering via low-rank and sparse decomposition[C]. Proceedings of the Twenty-Eighth AAAI Conference on Artificial Intelligence, Québec City, 2014: 2149–2155.
[10]	CAI Jianfeng, CANDÈS E J, and SHEN Zuowei. A singular value thresholding algorithm for matrix completion[J]. SIAM Journal on Optimization, 2010, 20(4): 1956–1982. doi: 10.1137/080738970
[11]	ZHU Ciyou, BYRD R H, LU Peihuang, et al. Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization[J]. ACM Transactions on Mathematical Software, 1997, 23(4): 550–560. doi: 10.1145/279232.279236
[12]	KUMAR Y, MULLER U, BEN J, et al. spectral clustering MIT Press[C]. 24th Annual Conference on Neural Information Processing Systems, Vancouver, Canada, 2010.
[13]	MATSUGU M, MORI K, ISHII M, et al. Convolutional spiking neural network model for robust face detection[C]. Proceedings of the 9th International Conference on Neural Information Processing, Singapore, Singapore, 2002: 660–664.

施引文献

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