在阵列信号处理中，传统的信源参数估计算法以点目标模型为基础。当目标的空间散射特性可以忽略、目标与接收阵列之间为直线路径, 点目标模型可以简化计算，有效进行参数估计。在无线通信领域，目标与接收阵列周围存在较多障碍物，目标到达阵列通过多条路径传播；在水下声呐工作环境中, 信号传播存在海底与海面的反射路径，随着目标与接收阵列之间的距离缩短，目标反射信号的各个部位都可以看作散射体，此时信源的参数估计就必须要考虑目标的空间分布特征。当信源的空间散射特性不能被忽略时，点目标模型不能有效描述信号源特征，基于点目标模型的信源参数估计也会存在较大误差。

为了克服点目标模型的局限，学者们提出了分布源模型^[1]。分布源信号可以看作点源信号在一定空间分布内的聚集。分布源依据其散射单元的时间相关性主要分为非相干分布源和相干分布源。非相干分布源和相干分布源的空间分布特征分别用角功率密度函数和角信号分布函数来描述。根据分布源的空间分布特点，角功率密度函数或者角信号分布函数可以表示为多种分布函数，通常有高斯分布、均匀分布和非对称分布^[2]。针对分布源模型，学者们提出多种参数估计算法，主要有子空间类算法如DSPE^[1]和DSPARE^[3]、协方差匹配算法^[4-5]、极大似然估计算法^[6-7]和波束形成算法^[8]。以上参数估计算法假设分布源与接收阵列处于同一平面即采用一维分布源模型，分布源的角功率密度函数或者角信号分布函数用波达方向或称为中心角以及角度扩展2个参数来描述。分布源与接收阵列处于不同平面即二维分布源模型更具有一般性，此时分布源用波达方向包括中心方位角和中心俯仰角、角度扩展包括方位角扩展和俯仰角扩展来描述，复杂程度显著增加。

对于二维相干分布源模型，代表性的成果为国内学者利用广义方向矢量的泰勒近似得出信号子空间的旋转不变关系，针对不同的阵列结构提出多种参数估计算法^[9-12]。对于二维非相干分布源模型国内外研究较少。文献[13]提出二维协方差匹配算法，该方法采用四维非线性优化求解二维中心角和角度扩展，计算量巨大。文献[14]利用均匀平面阵提出一种类ESPRIT算法，该算法针对二维非相干分布源求解中心角而没有涉及角度扩展参数。文献[15]利用双平行线阵、三平行线阵和L型阵对分布源的低复杂度参数估计做了系统的研究，并针对双平行阵列提出一种非相干分布源中心角估计算法，首先利用旋转不变关系计算出中心俯仰角，利用一维Capon谱峰搜索来计算中心方位角，再通过二维Capon谱峰搜索来得到扩展参数，计算复杂度也较高。

本文研究二维非相干分布源模型的参数估计。首先针对双L型阵列结构特点推导出阵列广义方向矢量关于中心方位角和中心俯仰角的旋转不变关系，在此基础上利用接收矢量协方差矩阵的信号子空间解算出旋转算子，依据Capon估计准则提出一种参数配对方法，从而得出中心俯仰角和中心方位角。本文二维波达方向的估计不需要谱峰搜索，适用于角功率函数未知和多个分布源的情景。在中心角已知的基础上，通过Capon谱峰搜索确定扩展参数。

1 二维非相干分布源模型

双L型阵列结构如图 1所示。阵列由位于xoz平面的两对平行线阵X1和X2、Z1和Z2组成, 其中X1和Z1分别位于x轴和z轴, 共有位于原点的第一个阵元。每对平行线阵间的间隔为δ, 每条线阵有间距为d的M个阵元。

假设远场有k个窄带二维分布源以中心角为(θ_i, ϕ_i) (i=1, 2, …, k)入射到阵列。其中θ_i为第i个分布源的中心方位角, ϕ_i为第i个分布源中心俯仰角, θ_i∈[0, π], ϕ_i∈[0, π]。则阵列的接收矢量可以表示为

(1)

式中, n_xm(t)与n_zm(t)为与信号不相关的加性高斯白噪声, 噪声功率为σ_n²。m=1, 2。x₁(t)和x₂(t)分别表示阵列X1和X2的接收矢量; z₁(t)和z₂(t)分别为阵列Z1和Z2的接收矢量; α₁(θ, ϕ)和α₂(θ, ϕ)分别表示点源信号在X1和X2的方向矢量; β₁(θ, ϕ)和β₂(θ, ϕ)分别表示点源信号在Z1和Z2的方向矢量。

(2)

(3)

s_i(θ, ϕ, t)是第i个分布源的角信号密度函数, s_i(θ, ϕ, t)为一随机过程, 反应第i个分布源在t时刻不同角度上散射分量的贡献。非相干分布源模型假设同一分布源不同角度的散射完全不相关, 则存在以下关系

(4)

式中，σ_i²是第i个分布源的功率, p_i(θ, ϕ, t)是归一化的角功率密度函数满足以下关系

(5)

2 算法原理 2.1 广义方向矢量的旋转不变关系

将阵列X1和X2的方向矢量在(θ_i, ϕ_i)处做一阶泰勒级数展开如下

(6)

(7)

式中, [·]′_θ和[·]′_ϕ分别表示函数在(θ_i, ϕ_i)处关于θ和ϕ的一阶偏导数。

(8)

如果δ/λ≪1, [α₂(θ_i, ϕ_i]′_ϕ的第二项可以忽略不计则存在以下近似关系

(9)

定义广义信源矢量s=[s₁, s₂, s₃]^H

(10)

式中

(11)

可以证明, 存在以下关系^[15]

(12)

式中

(13)

(14)

由于不同的分布源是不相关的, 则

(15)

(16)

定义X1阵列的广义方向矢量[A₁₁, A₁₂, A₁₃]为

(17)

则X1阵列的接收矢量可以表示为

(18)

定义X2阵列的广义方向矢量[A₂₁, A₂₂, A₂₃]为

(19)

则X2阵列的接收矢量可以表示为

(20)

根据(8)~(9)式可以得到如下关系

(21)

式中

(22)

将阵列Z1和Z2的方向矢量在(θ_i, ϕ_i)做一阶泰勒级数展开如下

(23)

(24)

式中

(25)

(26)

如果d/λ≪1, 则(25)式右端及(26)式第二项可以忽略不计则有以下关系

(27)

定义Z1阵列方向矢量[B₁₁, B₁₂]

(28)

则Z1阵列的接收矢量可以表示为

(29)

定义Z2阵列方向矢量[B₁₁, B₁₂]为

(30)

Z2阵列的接收矢量可以表示为

(31)

根据(27)式得出以下关系

(32)

(33)

2.2 中心角参数估计

将X1和X2阵列的接收矢量合并如下

(34)

式中

(35)

(36)

X1和X2阵列的接收矢量协方差矩阵可以表示为

(37)

在小角度扩展的条件下, M_θi和M_ϕi都是小于1的数, 对协方差矩阵R_x¹²进行特征分解, 得到协方差矩阵最大的k个特征值对应的特征向量构成的子空间E_x。E_x子空间与A_x1张成的子空间相同, 则存在一个k×k非奇异矩阵T满足以下关系

(38)

令E_x1和E_x2分别为E_x的上k行和下k行则有

(39)

(40)

可以得到以下关系

(41)

定义Ω_x=T^-1ΦT则

(42)

其中(·)⁺表示伪逆算子, 从而可以求解出分布源的中心俯仰角

(43)

式中, η_i为Ω_x第i个特征值。angle(·)表示取复变量的辐角。

将Z1和Z2阵列的接收矢量合并如下

(44)

式中

(45)

(46)

Z1和Z2阵列的接收矢量协方差矩阵可以表示为

(47)

对协方差矩阵R_z¹²进行特征分解, 得到与协方差矩阵最大的k个特征值对应的特征向量构成的子空间E_z。E_z子空间与B_z1张成的子空间相同, 则存在一个k×k个的非奇异矩阵Q满足以下关系

(48)

令E_z1和E_z2分别为E_z的上k行和下k行则有以下关系

(49)

(50)

进一步可以得到

(51)

定义Ω_z=Q^-1ΨQ则

(52)

进而得出

(53)

式中，μ_i为Ω_z第i个特征值。

2.3 角度配对

考虑平行阵列X1和X2的广义方向矢量[A_x1, A_x2, A_x3], Capon估计准则可以表示为

(54)

(54) 式可以通过拉格朗日函数进行求解, 得出Capon代价函数为

(55)

式中，为平行阵列X1和X2的协方差矩阵估计值

(56)

角度配对的步骤可以归纳如下

1) 从集合{ϕ₁, ϕ₂, …, ϕ_k}任选一, 遍历集合{angle(μ₁), angle(μ₂), …, angle(μ_k)}依据(53)式求解出k个与的角度配对(j=1, 2, …, k)。

2)将(j=1, 2, …, k)带入(55)式选择使L(θ, ϕ)最大的θ_j为的配对, 记为。

3) 将从集合{ϕ₁, ϕ₂, …, ϕ_k}剔除, 同时将从集合{angle(μ₁), angle(μ₂), …, angle(μ_k)}剔除, 重复步骤1~步骤2。

经过(k+2)(k-1)/2次计算就可以配对出所有的中心方位角和中心俯仰角。

2.4 扩展参数估计

考虑阵列X1, 信号协方差矩阵可以表示为

(57)

对于高斯分布

(58)

式中，u_i=[θ_i, ϕ_i, σ_θi, σ_ϕi], 假设d/λ=0.5且在小角度扩展的前提下有

(59)

式中，[·]_lh表示矩阵第l行h列的元素。

(60)

对于均匀分布

(61)

(62)

(63)

方位角扩展和俯仰角扩展可以采用Capon估计准则通过二维谱峰搜索得出

(64)

式中，表示矩阵的最大特征值。

2.5 算法步骤和复杂度分析

根据上述分析, 本文算法步骤归纳如下

1) 计算协方差矩阵的估计值。

2) 对进行特征分解得出子空间E_x和E_z, 根据(42)、(52)式计算得出Ω_x和Ω_z。

3) 对Ω_x和Ω_z进行特征分解得出特征值η_i和μ_i。

4) 依据参数配对方法确定中心角度。

5) 依据(59)式、(60)式或者(62)式、(63)式确定R(θ, ϕ; u_i), 根据(64)式得出扩展参数和。

本文算法求解中心角复杂度主要包括三部分:①计算协方差矩阵, 复杂度为O(NM²); ②特征分解, 复杂度为O(M³); ③求解及特征分解Ω_x和Ω_z, 复杂度为O(k³); 文献[15]提出的双平行线的估计方法复杂度包括三部分:①计算协方差矩阵复杂度为O(NM²); ②特征分解复杂度为O(M³)。③一维搜索复杂度为O(kL_θM²), 其中L_θ为方位角的搜索点数。显然, 本文提出的算法复杂度低于文献[15]。在得出中心角的基础上求解扩展参数的计算复杂度主要包括两部分:①协方差矩阵求逆, 复杂度为o(M³)；②通过二维搜索得出扩展参数, 复杂度为o(kM³L_σθL_σϕ), L_σθ和L_σϕ为方位角扩展和俯仰角扩展搜索的点数。

3 仿真实验与分析

本文通过3个仿真实验来验证算法的有效性, 考虑阵列结构如图 1所示, 每条线阵阵元数M=8, 阵元间距d=λ/2, 平行阵列间隔δ=λ/10。

中心角的均方根误差的G_a定义为

(65)

M_c为蒙特卡罗仿真次数, 表示第ζ次蒙特卡罗仿真实验估计的中心方位角和中心俯仰角。角度扩展的均方根误差G_s定义为

(66)

表示第ζ次蒙特卡罗仿真实验估计的方位角扩展和俯仰角扩展。

实验1  考察信噪比和快拍数对实验的影响。假设目标参数为[30°, 45°, 2°, 3°]的二维非相干高斯分布源信号。实验采用100次独立的蒙特卡罗仿真验证完成。图 2反应快拍数等于200时中心角均方根误差和角度扩展均方根误差随信噪比从0~30 dB的变化情况。图 3反应信噪比等于20 dB时中心角均方根误差和角度扩展均方根误差随快拍数从20~200的变化情况。为了验证算法的性能, 图中还给出采用文献[15]的方法利用双平行线X1和X2得出的结果。可以看到随着快拍数和信噪比的提高算法估计精度不断提高, 本文算法在相同实验条件下的估计精度优于文献[15]提出的算法。

实验2  考察角度扩展的变化对中心角估计的影响。针对中心角度为(60°, 35°)一个非相干二维分布源, 分别考虑2种角功率密度函数:高斯分布和均匀分布。考察方位角扩展和俯仰角扩展同时从0°变化到10°时中心角均方根误差的变化。实验采用100次独立的蒙特卡罗仿真验证完成, 信噪比为20 dB, 快拍数为200。图 4表明对于角功率密度函数为高斯分布的分布源, 随着角度扩展的增加估计精度会降低, 角度扩展为5°时, 中心角均方根误差为0.03;角度扩展达为10°时误差为0.4。对于均匀分布的分布源, 随着角度扩展的增加估计精度同样降低, 角度扩展为5°时, 中心角均方误差为0.04, 10°时误差为0.7。实验表明随着角度扩展的增加估计性能会下降, 但是在10°以内仍具有满意的估计效果。本文提出的算法中心角度的估计在小角度扩展的情况下具有鲁棒性。

实验3  考察算法对多个分布源的估计能力。考虑参数为[25°, 70°, 2°, 3°], [35°, 40°, 3.5°, 4.5°]和[55°, 80°, 1.5°, 2°]的3个均匀非相干二维分布源，信噪比为20 dB, 快拍数为200。实验采用100次独立的蒙特卡罗仿真完成。图 5和图 6分别给出3个均匀分布源中心角和角度扩展估的结果。实验结果表明本文算法可以有效估计多个分布源的参数。

需要说明的是本文旋转不变关系的得出是建立在假设平行线阵间隔很小, 与信号波长之比远小于1即δ/λ≪1的条件之上。在工程实际中如果阵元间距很小, 一方面阵元间的互耦效应将不能被忽略; 另一方面, 阵元的安装误差也会影响估计效果, 所以在一些高频信号领域如移动通信很难适用。但是在低频水下探测领域, 声纳的探测频率可以降到100 Hz, 这意味着信号波长可以达到14.5 m; 在这种情况下, δ/λ实际上可以达到很小的值, 完全可以满足假设条件。

4 结论

本文利用双L型阵列，提出一种二维非相干分布源参数估计方法。首先推导出广义方向矢量关于中心角的在双L型阵列下的旋转不变关系，在协方差矩阵特征分解的基础上计算得出旋转算子，参数配对后解出中心角。在中心角已知的基础上通过二维搜索得出角度扩展参数。仿真结果表明，在平行线阵间隔很小的条件下，本文提出的参数估计方法具有良好的估计性能。