第5章基于模拟信息转换的阵列信号空频域参数联合估计

随着现代科技与信息技术的发展，各领域中信号的带宽也越来越宽。尤其在阵列中，多个阵元传感器同时接收信号的情况下，多个阵元均向后端数据处理中心发送原始信号，如此大量的数据传输将不可避免地引起诸如传输延迟大和高功耗的问题。因此，尽量降低采样速率又能实现基于阵列的空频域参数的联合估计，具有非常重要的现实应用意义。本章将MWC技术与阵列信号处理相结合，在压缩采样下，利用少量的采样点实现信号源的二维DOA与载频参数的联合估计，确定信源方向的同时，得到信号源准确的时域波形，这有重要的实际应用意义。

5.1相关研究综述
5.1.1阵列信号空频域参数联合估计研究现状

阵列信号处理的理论发展从20世纪60年代开始，Howells于1965年提出了自适应陷波的旁瓣对消器； 其后各位学者开始在该方面进行大量研究并取得了相应的成果，其中Schmidt在1979年提出的多重信号分类算法(Multiple Signal Classification，MUSIC)以及Roy等在1986年提出的旋转不变子空间(Estimating Signal Parameter via Rotational Invariance Techniques，ESPRIT)算法，实现了从测向算法向子空间类超分辨方法的跃进，极大地促进了阵列信号处理理论的发展。21世纪初， Sidiropoulos教授将阵列接收信号模型建模为正则分解/平行因子分析(CANDECOMP/PARAFAC，CP)模型，将CP分解方法引入阵列信号处理领域，为阵列测向问题提供了新的解决思路。
电磁波信号的空域信息参数包括信号的来波方向、距离等，而频域的参数包括信号的载频、带宽、相位、能量、多普勒频率和极化参数等。分析认识信号就必须从这些参数入手，在通信、雷达和认知无线电等领域中，信号的频率和到达角是对电磁波信号进行识别的重要特征，因此研究对频率和到达角的估计是现代阵列信号处理领域的重要内容。
空间谱估计是阵列信号处理的重要内容，对空间中目标的DOA估计是传统的空间谱估计主要内容，它利用空间阵列的不同阵元的接收数据之间的相位差来估计目标信号的DOA参数，一般假设信号源在一个平面内，并且假设所有入射到接收阵列的信号的频率已知且相同。这种假设显然不适合现代实际的雷达、通信等环境要求，因此需要进行阵列多参数估计，以同时得到多个不同目标信号各自的方位信息和频率等特征参数。其中，在雷达、通信等领域的应用中信号频率与DOA联合估计问题出现比较广泛，一般利用线性预测方法、多维MUSIC方法、最大似然方法、ESPRIT方法和CP分解方法等估计方法来进行求解。其中，ESPRIT算法是一种利用均匀阵列接收信号之间的旋转不变特性来对信号参数的数学表达式进行求解的方法，不需要进行谱峰搜索操作，因此计算量相对较小，研究也更为丰富。Lemma提出了基于多维ESPRIT算法的到达角与频率联合估计算法； Viberg等提出了适用于任意阵列的在色噪声情况下的参数估计方法，采用状态空间实现技术实现载频估计，再利用ESPRIT算法进行方位角的估计，并将两个参数配对。国内葛利嘉等人针对窄带信号，利用旋转不变技术，实现了方位角/频率的联合估计方法； 徐友根等人通过时延处理构造伪数据矩阵，通过子空间旋转获得频率和DOA的同时估计。
CP分解的提出是为了解决心理测量学中多维数据分析问题，在分析化学、心理学、食品学等领域均有广泛的应用，直到Sidiropoulos将其引入通信及信号处理领域。基于CP分解的信号处理技术同样利用了均匀阵列结构的旋转不变特性，相较于ESPRIT算法，虽然增加了计算量，但鲁棒性更好，因此近年来备受学者关注。张小飞等人将CP分解方法应用到均匀线阵参数估计问题中，提出了一种角度与频率的盲估计方法。Xu L中提出了一种利用酉变换和正反向技术的基于三线性分解的DOA和载频联合估计方法。
阵列测向理论的基础是阵列结构，为了同时获得信号的频率和到达角信息，需要使用二维阵列对信号进行接收。目前应用最广的二维阵列是平面阵，但平面阵一般来说体积较大，所需阵元多，相应地，需要的设备多、成本高，并且计算量大，不利于阵列信号处理的实时性，因此其他不同的阵列形式例如L形阵列、T形阵列、圆形阵列和十字阵列等结构的研究越来越受到学者的关注。其中，L形阵列的设计最简单，同时阵列方向图主瓣最窄，估计性能较好，因此本章主要基于L形阵列结构进行研究。
但在实际的雷达、通信等系统中，目标并非均位于同一个平面中，因此需要二维DOA，即俯仰角和方位角来联合表示，以准确描述三维空间中目标信号的方向。因此，阵列信号的频率与二维DOA联合估计是阵列信号空频域多参数估计中的一个重要研究内容，使接收侦察机可以同时得到多个信号源的方位和频率信息。一般通过对阵列形式及算法的调整，可以实现从一维DOA载频联合估计到二维DOA载频联合估计的拓展。Haardt等人讨论了移动通信领域中信号频率与二维DOA的联合估计，以解决空分多址技术的问题； Peter Strobach在二维均匀矩形阵列结构上，利用TLSESPRIT算法对频率及二维DOA进行联合估计； Zoltowski等结合非均匀L形阵列研究了欠采样条件下信号频率和二维DOA联合估计的ESPRIT方法。国内学者们也纷纷展开了研究。梁军利计算接收数据的高阶累积量，并用其构造平行因子模型，使用交替三线性最小二乘方法联合估计了信源频率及二维DOA参数。许凌云将平行因子四线性模型分解技术应用到均匀面阵列中，利用该模型低秩分解盲估计信源的二维DOA和频率，实现了参数同时估计与配对。孙晓颖同时进行了时域采样和空域采样，并利用时空的采样值构造了二维时空矩阵，提出了一种基于2DESPRIT方法的频率和二维到达角联合估计算法，适用于任意阵列。Wang L提出了一种基于CP分解模型，采用共形阵列结构的联合频率和二维到达角估计算法，需要4个位置精确的引导阵元。

5.1.2基于模拟信息转换的阵列信号参数联合估计
研究现状
以上研究虽然解决了阵列信号参数估计的问题，但随着现代信息技术越来越发达，如果以传统的香农奈奎斯特采样定理采样，会造成前端ADC设备的采样速率过高，后端设备存储和处理压力大等问题，如雷达探测、脉冲超宽带通信等领域中采样速率都在几GS/s以上。这样高采样速率的ADC，以目前的半导体集成工艺很难做到低成本与低功耗。为了能够在尽可能降低资源消耗的前提下进行信号采集与数据处理，很多学者对欠奈奎斯特采样方法进行了研究。例如多片ADC时间交错采样技术、等效采样技术、多带信号处理方法和频率信道化方法等； 还有一些基于欠采样的解模糊方法，例如利用中国余数定理的解模糊方法、在多个延迟通道结构上利用前后向稀疏线性预测的解模糊方法、利用ProESPRIT算法的解模糊方法等。
以上方法虽然在一定程度上降低了采样速率，基于压缩感知的欠采样方法为欠采样的研究开辟了新的方向。为了降低阵列信号处理的采样、数据处理及存储的压力，有必要将压缩感知理论与阵列信号处理技术相结合，发展基于压缩感知的阵列参数估计方法。目前已有部分学者进行相关研究并取得一定成果，例如2006年，Malioutov等在估计信源角度时将其离散化，利用了稀疏的思想来建立信号重构模型。国内研究方面，2014年，沈志博等人提出一种基于CS理论的信号频率和一维DOA联合估计的CSSVD算法； 陈玉龙提出了一种基于L形阵列的二维估计压缩感知算法； Le X将多时延输出数据构建为CP分解模型，然后进行分区压缩，通过三线性交替最小二乘法实现二维和频率联合估计。Fang J和Wang F利用每个天线的时间延迟，实现了基于CP分解的一维DOA、载波频率和功率谱的联合恢复方法。作者Liu L等提出了基于多陪集采样的窄带远场信号的频率和一维DOA联合估计方法，分别基于子空间分解和CP分解方法，并推导了克拉美罗界(CramérRao Bound，CRB)。接着，他们又进行了扩展，进一步减少了采样结构所需的阵元个数。
2009年，Mishali M和Yonina C. Eldar等针对多频带信号，提出基于压缩感知理论的调制带宽转换器系统，通过周期混合函数建立了连续与离散之间的关系，混频模块将信号的频谱全部搬移到基带内再进行处理，易于硬件实现，结合压缩感知理论压缩与采样同时进行的优点，可以在信号载频未知的情况下通过少量的采样点高效地恢复信号的频谱信息，精确地重构原始信号，这是一个巨大的进步。2015年，Shahar Stein等将MWC技术与DOA估计结合，提出了基于均匀线性阵列MWC结构的信号一维DOA估计方法。紧接着，他们提出了一种基于L形阵列的CP分解方法，能同时实现信号载频和一维DOA的估计。2017年，他们具体讨论了均匀线阵以及L形阵列的MWC系统，并利用旋转不变技术和压缩感知方法成功地恢复了信号二维DOA、载频以及信号时域波形。2018年，Esmaeil R在阵列MWC结构的基础上，提出了一种基于压缩感知方法的二维DOA与载频联合估计方法，然而参数估计精度取决于网格的大小； 还提出了一种不依赖网格的贝叶斯算法，但多维数据的组合矩阵体积较大，计算复杂。国内Tao C提出了一种基于均匀线性MWC阵列的一维DOA和载频联合估计方法，由两阶段估计和相应的参数配对方法组成，然而配对操作十分复杂，精度有限。张小飞提出了一种基于多个MWC结构的一维DOA和载频联合估计方法，提升了鲁棒性，但该方法重构所需的通道数过多，硬件成本高。
综上所述，由于射频多频带信号频域稀疏的特点和阵列信号处理庞大的存储量与运算量，与压缩感知结合有十分重要的应用意义，可以降低对ADC硬件的要求，解决采集数据量巨大难以存储的问题。然而，目前国内外学者对基于阵列MWC结构的多频带信号一维DOA和频率联合估计已经有了充分的研究，但在基于阵列MWC结构的多频带信号二维DOA和频率联合估计方面研究较少，因此，在压缩采样条件下对多频带信号的二维DOA和频率参数联合估计还有待研究。
5.2基于双L型阵列MWC的二维DOA和载频联合估计方法
5.2.1阵列多频带信号模型

阵列多频带信号在雷达探测、超宽带无线通信及认知无线电等领域都有着广泛的应用，多个不同窄带目标信号经各自发射机分别调制后入射到天线接收阵列。此时，接收阵列所接收到的信号的组合即可以看作一个多频带信号。图51为阵列信号示意图。


图51阵列信号示意图


为了简化复杂信号传输模型，便于阵列接收信号的分析，一般对目标信号和接收阵列作如下几个假设。
(1) 窄带假设： 假设信号的带宽远小于信号的载频，即接收阵列中不同阵元之间的接收的信号包络相同，不同传感器阵元接收的信号只有相位不相同。
(2) 假设信号之间互不相关，互相关系数为0。
(3) 平面波假设： 假设目标到接收阵列之间的距离远大于阵列的直径，此时入射到阵列的电磁波可以近似为平面波，即目标信号到达每个阵元的角度均相同。
(4) 忽略阵元之间的互耦，假设接收阵列中的不同传感器阵元之间的工作均处于理想状态，互相之间不存在电磁耦合效应。
(5) 各阵元的增益相等，噪声独立且为加性高斯白噪声。

在以上阵列假设的条件的基础上，设有M个信源发出的幅值为复数的连续时间窄带信号si(t)，其中i∈{1,2,…,M}，调制载频为{fi}Mi=1∈R。假定所有的窄带信号{si(t)}Mi=1均为远场的、广义平稳、均值为零并且互不相关的信号，即对于i≠j，有E［si(t)s-j(t)］=0
，并且σ2i=E［s2i(t)］≠0。假定目标信号满足远场窄带信号模型，存在于三维空间中，类似于球面坐标系内两个角度的定义，定义第i个目标与原点之间的连线在xy平面的投影与x轴正方向的夹角为方位角θi，定义目标与原点之间的连线与其在xy平面的投影之间的夹角为俯仰角φi
，两个角度即可以决定目标的唯一方向。对于i=1,2,…,M，有θi∈(-90°,90°),φi∈(0°,90°)。为了避免阵列模糊，假设对于i≠j，有： 


ficosθicosφi≠fjcosθjcosφj

fisinθicosφi≠fjsinθjcosφj

fisinφi≠fjsinφj(51) 
因此，M个信源发出的信号{si(t)ej2πfit}Mi=1入射到阵元上，阵元接收的信号可以看为M个窄带信号的叠加，可以看作一个多频带信号。假设窄带信号{si(t)}Mi=1的频谱范围有限，其频率范围为B=-12T,12T，调制后的信号{si(t)ej2πfit}Mi=1的频谱互相不重叠。每个窄带信号的带宽不超过B且远小于信号载频，满足前文提到的窄带假设条件。每个阵元接收信号的频率范围为F=-fNyq2,fNyq2，其中fNyq为信号的奈奎斯特频率，定义为信号最高频率的两倍。

5.2.2双L型阵列MWC系统
1. 双L型阵列MWC结构

双L型阵列结构如图52所示，沿x轴、y轴和z轴正方向均有N个均匀分布的传感器阵元{x1，x2,…,xN}、{y1,y2,…,yN}和{z1,z2,…,zN}，


图52双L型阵列结构

3个互相垂直的轴在原点共用一个传感器。为了避免空间角度多值模糊的情况，每个轴的相邻传感器阵元之间的间距应满足d≤cfNyq，其中c为光速。每个传感器阵元后连接一个MWC通道结构，不同阵元接收的信号首先与同一个周期性伪随机混频序列相乘，之后再通过低通滤波器滤除高频成分，只保留基带内的信号，最后进行低速采样。

为了使三轴接收信息相对独立，表示更加简洁，定义与x轴的夹角为αi，与y轴的夹角为βi，与z轴的夹角为γi。由空间几何关系可知角度之间的关系如式(52)和式(53)所示： 
cosαi=cosθicosφi

cosβi=sinθicosφi

cosγi=sinφi(52)

cos2αi+cos2βi+cos2γi=1(53)

因此，可以用信号与x轴正方向的夹角αi以及与y
轴正方向的夹角βi来进一步表示方位角θi和俯仰角φi： 
θi=arctancosβicosαi(54)

φi=arccos(cos2αi+cos2βi)(55)
则x轴的第n个传感器接收的信号xn(t)可以表示为
uxn(t)≈∑Mi=1si(t)ej2πfi(t+τxn(αi))(56)
其中，τxn(αi)=dnccosαi，表示x轴第n个传感器与第一个传感器接收信号之间的相位差。类似的有y轴和z轴传感器接收的信号yn(t)和zn(t)，其中τyn(βi)=dnccosβi，τxn(γi)=dnccosγi。

2.  频域分析
下面对系统进行频域上的分析，首先定义： 

fp=1/Tp,Fp=［－fp/2,fp/2］


fs=1/Ts,Fs=［－fs/2,fs/2］(57)

由于信号满足窄带假设条件，信号的带宽远小于其中心频率，即满足B/fi<0.1，因此阵元接收之间的信号包络变化相对缓慢，信号包络的延迟远小于1/B，其在各阵元上的差异可以忽略，即si(t+τn)≈si(t)。因此信号延迟只反映在相位变化上。则L型阵列的x轴上的第n个传感器xn接收的信号uxn(t)可以表示为式(58)的形式： 
uxn(t)=∑Mi=1si(t+τxn(θi))ej2πfi(t+τxn(θi))≈∑Mi=1si(t)ej2πfi(t+τxn(θi))
(58)
其中： 
τxn(θi)=dnccos(θi)
(59)
表示第n个传感器与第一个传感器接收信号之间的相位差。接收信号xn(t)的傅里叶变换形式为

Uxn(f)=∑Mi=1Si(f－fi)ej2πfiτxn(θi)
(510)



图53阵元接收信号频谱Uxn(f)示意图

其中，Si(f)是si(t)的傅里叶变换形式，当fF时，Si(f)=0。阵元接收信号uxn(t)可以看作一个多频带信号，如图53所示为阵元接收信号频谱Uxn(f)示意图。 
各传感器阵元接收的信号先与同一个周期伪随机序列p(t)混频，p(t)的周期为Tp=1/fp，其傅里叶级数形式为
p(t)=∑∞l=－∞clej2πlfpt
(511)
其中，系数cl可由式(512)求得
cl=1Tp∫Tp0p(t)ej2πlfptdt
(512)
所以混频后的信号x～n(t)=xn(t)p(t)的傅里叶变换为

X～n(f)=∫∞－∞uxn(t)p(t)e－j2πftdt

=∫∞－∞uxn(t)∑∞l=－∞clej2πfplte－j2πftdt

=∑∞l=－∞cl∫∞－∞uxn(t)e－j2π(f－lfp)tdt

=∑∞l=－∞clUxn(f－lfp)
(513)
可见，待测信号与周期伪随机序列混频后，相当于在时域相乘，输出信号的频谱是输入多频带信号的频谱以fp进行整数倍移位后的加权组合，权重是周期序列p(t)的傅里叶级数系数cl，此时输出信号的基带内就包含待测信号的频谱信息。
将式(510)代入式(513)中可得

X～n(f)=∑∞l=－∞cl∑Mi=1Si(f－fi－lfp)ej2πfiτxn(θi)(514)
经过低通滤波器后，混频后的信号仅留下了Fs内的频谱，因此可得
Xn(f)=∑L0l=－L0cl∑Mi=1Si(f－fi－lfp)ej2πfiτxn(θi)

=∑Mi=1ej2πfiτn∑L0l=－L0clSi(f－fi－lfp)

=∑Mi=1S~i(f)ej2πfiτxn(θi)
(515)
其中，L0是使式(515)包含Yn(f)在F内，即可以保证奈奎斯特频率范围内的所有频率被搬移到基带的所有非零项的最小正整数，可以由计算得
L0=fNyq2fp(516)
且有
S～i(f)=∑L0l=－L0clSi(f－fi－lfp)(517)
混频过程和滤波过程的信号频谱示意图如图54所示，混频序列p(t)以周期频率fp为基频的倍频谐波分量，可以将被测信号的奈奎斯特频率分为L=2L0+1
个子频带，每个子频带宽度为fp，以lfp频率为中心，其中l=［－L0,…,L0］。混频步骤将信号不同子频带上的信号进行移位加权，当fs=fp时，截止频率为fs/2滤波器刚好将Fp内的频率保留。从图54中可以看出，Fp内的频谱保留了原始多频带信号的全部频谱信息。在滤波后需要满足奈奎斯特采样定理进行采样才能完全保留基带内的频谱信息，因此采样速率fs需要满足fs≥B，即每通道采样速率fs和混频序列p(t)的周期频率fp都需要大于或等于信号中单个频带的最大带宽B。 


图54混频过程和滤波过程的信号频谱示意图


经低速采样后，第n个通道的信号xn［k］=xn(kTs)的离散傅里叶变换可表示为
Xn(ej2πfTs)=∑Mi=1Wi(ej2πfTs)ej2πfiτxn,f∈Fs(518)
定义wi［k］=S～i(kTs)，则Wi(ej2πfTs)=DTFT{wi［k］}，表示信号被搬移到基带后的频谱。
因此，对于x轴，可以将式(518)写成如下矩阵形式： 
X(f)=AxW(f)，f∈Fs(519)
其中，X(f)是N×1的矩阵，第n个元素为Xn(f)=Xn(ej2πfTs)。未知向量W(f)是M×1的矩阵，第i个元素为Wi(f)=Wi(ej2πfTs) 。矩阵Ax为空间阵列流型矩阵，是一个N×M的矩阵，其中元素由未知载频{fi}Mi=1，到达角{θi}Mi=1决定，其中元素定义为
Ax=ej2πf1τx1(θ1)…ej2πfMτx1(θM)



ej2πf1τxN(θ1)…ej2πfMτxN(θM)N×M(520)

同理，对于z轴传感器阵列接收信号的采样序列也可以写为矩阵形式： 

Z(f)=AzW(f)，f∈Fs
(521)
其中，z轴方向的阵列流型矩阵Az与式(520)中定义的矩阵Ax中的元素类似，不同的是有τzn(θ)=dncsin(θ)。由于矩阵Ax与Az均由未知载频f和到达角θ决定，因此可以写作Ax=Ax(f,θ)，Az=Az(f,θ)。
对式(519)和式(521)中定义的采样序列频率的矩阵形式取离散时间傅里叶逆变换(Discretetime Fourier Transform，DTFT)，在时域有类似的矩阵形式： 

x［k］=Axw［k］，k∈Z(522)

z［k］=Azw［k］，k∈Z(523)
其中，x［k］和z［k］是x轴和z轴的观测值； w［k］是长度为M的向量，第i个元素为wi［k］。根据式(522)和式(523)，利用ESPRIT方法和CP分解方法就可以从低速采样序列x［k］和z［k］中恢复载频{fi}Mi=1和到达角{θi}Mi=1，进而可以求得阵列流型矩阵Ax与Az。代入式(522)和式(523)中获得信号被搬移到基带后的波形{wi(t)}Mi=1，最后经过相应的频带搬移和频率调制获得窄带目标发射信号{si(t)ej2πfit}Mi=1。下面首先介绍得到阵列流型矩阵后如何恢复窄带目标的基带信号{si(t)}Mi=1。
假设对于所有的l∈{－L0,…,L0}，cl≠0。考虑第i个窄带信号，对于f′∈BFp，有

Wi(ej2πfTs)= S～i(f′)=claSi(f′－fi－lafp)
(524)
由于cla≠0，则

Si(f′－fi－lafp)=1claWi(ej2πfTs)(525)
替换变量后可得
Si(f′)=1claWi(ej2π(f′+fi+lafp）Ts)
(526)
其中索引la满足
la=fi+f′+fp/2fp
(527)

利用求得的信号被搬移到基带后的波形{wi(t)}Mi=1，通过式(526)中描述的关系求得{Si(f)}Mi=1后就可以通过傅里叶反变换求得窄带目标信号的时域波形和调制后的发射波形。值得注意的是，由于窄带信号Si(f)
可能会被定义的间隔为fp
的频带划分为两个子频带，式(527)中定义的索引la仅可索引被分割的频带中的一个子频带，也就是对于一个窄带信号si(t)，可能会有1~2个索引la来分别索引被分割的两个子频带。

5.2.3基于ESPRIT的参数估计方法
1.  参数估计方法

将x轴上的传感器阵元接收信号的采样值阵列分为两个子阵列，第一个子阵由传感器{1,…,N-1}接收信号的采样值组成，第二个子阵由传感器{2,…,N}接收信号的采样值组成，对y轴和z轴上的传感器阵列接收信号的采样值进行类似的操作，可得

x1［k］=Ax1w［k］，x2［k］=Ax2w［k］

y1［k］=Ay1w［k］，y2［k］=Ay2w［k］

z1［k］=Az1w［k］，z2［k］=Az2w［k］(528)
其中，向量x1［k］和矩阵Ax1是向量x［k］和矩阵Ax的前N-1行; x2［k］和Ax2是向量x［k］和矩阵Ax的后N-1行，同理，y1［k］、y2［k］、z1［k］
、z2［k］和Ay1、Ay2、Az1、Az2定义类似。
每个轴的两个传感器阵列接收信号的两个子阵列之间仅有一个相位差，其阵列流型矩阵之间有如式(529)所示的旋转不变关系： 


Ax2=Ax1Φx

Ay2=Ay1Φy

Az2=Az1Φz
(529)
其中，Φx、Φy
和Φz为旋转矩阵，是M×M的对角阵，表示各轴传感器接收信号之间的相位差，其中元素如式(530)所示： 

Φx=diagexpj2πf1dccosα1,…,expj2πfMdccosαM

Φy=diagexpj2πf1dccosβ1,…,expj2πfMdccosβM

Φz=diagexpj2πf1dccosγ1,…,expj2πfMdccosγM
(530)
通过观察式(530)可以发现，未知载频{fi}Mi=1
，方位角{θi}Mi=1
和俯仰角{φi}Mi=1
可以从旋转不变矩阵Φx
、Φy
和Φz
中计算。接下来介绍如何从低速采样序列x［k］、y［k］
和z［k］
中获得旋转不变矩阵Φx
、Φy
和Φz
。考虑将三维参数估计问题转化为两个二维参数估计问题，再通过相应的配对操作得到对角元素顺序对应的旋转不变矩阵。
首先利用x
轴和y
轴的传感器样本得到矩阵Φx和Φy，同时实现Φx
和Φy对角元素的配对。计算如式(531)的相关矩阵： 

R1=E{x1yH1}=Ax1RwAHy1

R2=E{x2yH1}=Ax2RwAHy1=Ax1ΦxRwAHy1

R3=E{x1yH2}=Ax1RwAHy2=Ax1ΦHyRwAHy1

R4=E{x2yH2}=Ax2RwAHy2=Ax1ΦxΦHyRwAHy1
(531)
其中，E{·}表示取数学期望，Rw=E{WWH}为信源相关矩阵，由于信源相互独立，Rw可以表示成一个对角阵。构造如式(532)的协方差矩阵R： 

R=［R1; R2; R3; R4］
(532)
相比于一维DOA与载频联合估计方法，增加了一个相关项R4来提升估计的准确性。对利用采样值计算得到的协方差矩阵R执行奇异值分解操作： 

R=［U1U2］Λ0

00VH(533)

协方差矩阵R的左奇异向量可以组合为［U1U2］； 对角阵Λ对应了矩阵R的M个非零的奇异值； 对应地，协方差矩阵R的右奇异向量组成矩阵V。其中M个非零的奇异值对应的左奇异向量U1张成信号子空间，其余的奇异值对应的左奇异向量U2张成噪声子空间。显然两个子空间是正交的，有RHU2=0。而由矩阵Ax1、Ay1和Az1的列向量所张成的子空间也与噪声子空间正交，同时满足列满秩的条件。因此由阵列流型矩阵的列向量张成的子空间和信号子空间是同一个子空间，存在一个M×M维非奇异阵T使得式(534)成立： 

U=U1

U2

U3

U4=Ax1

Ax1Φx

Ax1ΦHy

Ax1ΦxΦHyT(534)
其中，U是一个4（N-1）×M的矩阵； Ui为（N-1）×M的矩阵，i=1,2,3,4。从式(534)可以看出矩阵Ui和对角阵Φx和Φy之间的关系。为了书写简洁，定义如式(535)的矩阵V1和V2： 


V1=U1

U3U2

U4=TΦxT－1(535)

V2=U1

U2
U3

U4=TΦHyT-1(536)
通过对矩阵V1+V2进行特征值分解操作，可以得到对应的特征向量矩阵T^，利用式(537)即可计算得到两个特征值的顺序正确对应的旋转不变矩阵Φ^x和Φ^y： 


Φ^x= T^-1V1T^

Φ^y=(T^-1V2T^)H(537)

类似地，利用x轴和z轴的采样值可得两个特征值的顺序正确对应的特征值矩阵Φ^′x和Φ^′z。

2.  配对方法
目前通过两组特征值分解可以得到两组内部顺序对应的特征值{Φ^x,Φ^y}和{Φ^′x,Φ^′z}，但此时两组特征值之间的顺序并不对应，因此还需要进行配对操作，通过调整对角阵Φ^′x中特征值的顺序来实现对两组特征值{Φ^x，Φ^y}和{Φ^′x，Φ^′z}进行三维参数配对。寻找一个置换矩阵Ξ满足： 
Ξ^=argminΞ‖Φ^x-Φ^′xΞ‖s.t.Ξi,j={0,1}(538)
利用置换矩阵Ξ^根据特征值Φ^′x的顺序相应的对矩阵Φ^′z
中的元素的顺序进行调整，可通过式(539)计算配对后的矩阵Φ^z： 

Φ^z=Φ^′zΞ^(539)

此时就得到了配对成功的旋转不变矩阵{Φ^x,Φ^y,Φ^z}，其中元素顺序完全对应，此时就可以利用估计的矩阵{Φ^x,Φ^y,Φ^z}
来计算目标信号的载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1参数。
令ui、vi和wi分别为Φx、Φy和Φz的第i个对角线元素，i=1,2,…,M，则有


2πfidccosαi=angle(ui)

2πfidccosβi=angle(vi)

2πfidccosγi=angle(wi)

(540)
其中，angle(·)为求复数的相位角。通过将式(53)和式(540)联立，可得

fi=c·angle2(ui)+angle2(vi)+angle2(wi)2πd
(541)
联立式(54)、式(55)和式(540)可得

θi=arctanangle(vi)angle(ui)(542)

φi=arccosangle2(ui)+angle2(vi)angle2(ui)+angle2(vi)+angle2(wi)(543)

因此，通过估计得到的旋转不变矩阵中的元素的值即可通过式(541)、式(542)和式(543)计算载频{fi}Mi=1、方位角{θi}Mi=1
和俯仰角{φi}Mi=1
的值。一旦载频、方位角和俯仰角的值确定，即可利用估计的参数计算矩阵Ax： 


Ax=111

ej2πf1dccosθ1cosφ1…ej2πfMdccosθ1cosφ1



ej2πfMd(N-1)ccosθMcosφM…ej2πfMd(N-1)ccosθMcosφM
(544)

Ay和Az的计算类似。得到阵列流型矩阵后即可通过式(545)恢复信号： 

W(f)=Ax

Ay

AzX(f)

Y(f)

Z(f)(545)
得到向量W(f)后，原信号{Si(f)}Mi=1可以用式(526)求得。
综上所述，将基于双L形阵列MWC的ESPRIT二维DOA与频率联合估计算法总结如下。
输入： x轴、y轴和z轴每个传感器接收信号的Q个采样快拍x［k］、y［k］和z［k］。
输出： 载频{fi}Mi=1、方位角{θi}Mi=1、俯仰角{φi}Mi=1和窄带信号频谱{Si(f)}Mi=1。
步骤1： 利用快拍值x［k］和y［k］根据式(531)和式(532)计算协方差矩阵R； 
步骤2： 对R进行奇异值分解，得到左奇异向量U，并根据式(535)和式(536)计算矩阵V1和V2； 
步骤3： 对V1+V2进行特征值分解求得一对旋转不变矩阵{Φ^x,Φ^y}； 
步骤4： 利用快拍值x［k］和z［k］重复步骤1~步骤3，求得一对旋转不变矩阵{Φ^′x,Φ^′z}； 
步骤5： 根据式(538)和式(539)对旋转不变矩阵中的元素进行配对，得到配对后的矩阵{Φ^x,Φ^y,Φ^z}； 
步骤6： 根据式(541)~式(543)计算载频{fi}Mi=1、方位角{θi}Mi=1、俯仰角{φi}Mi=1； 
步骤7： 利用估计的参数计算阵列流型矩阵Ax、Ay和Az； 
步骤8： 根据式(545)和式(526)计算原窄带目标信号的频谱{Si(f)}Mi=1。

5.2.4基于CP分解的参数估计方法
1.  CP分解模型

张量分解最早在1927年由Hitchcock提出，1944年Cattell在此基础上提出了多维数据分析的模型。直到20世纪60年代，这些概念才被重视起来，并且广泛应用于化学计量学、心理学等领域。直到1998年，Sidiropoulos将其引入通信及信号处理领域后才被信号处理领域的学者关注，并得到大量研究。下面先简单介绍张量分析的基础及CP分解的模型及具体求解方法。
1) 张量基础
张量是多维的数组，其阶数被称为维数或模。一个N阶张量χ∈CI1×I2×…×IN是由N个有自己的坐标系的向量空间的张量积组成的一个元素。例如，零阶张量相当于标量，一阶张量相当于平面中的向量，而二阶张量则相当于常用的二维的矩阵。图55为三阶张量χ∈RI×J×K
示意图。

下面简单介绍几个张量的概念。
N阶张量χ∈CI1×I2×…×IN的范数类似于矩阵的F范数，相当于将它其中所有元素的平方相加后再开方的值： 

‖χ‖=∑I1i1=1∑I2i2=1…∑INiN=1x2i1i2…iN
(546)
其中，xi1i2…iN为张量χ中第(i1,i2,…,iN)位置的元素。
秩1张量是一种特殊的张量类型，一个N阶的秩1张量χ∈CI1×I2×…×IN可以由N个向量的外积来完全表示，如式(547)所示： 

χ=a(1)a(2)…a(N)(547)

其中，表示向量之间的外积； a(n)为In维的向量，n=1,2，…,N。这意味着张量中的每个元素都可以由对应位置的向量中元素的乘积来表示： 

xi1i2…iN=a(1)i1a(2)i2…a(N)iN,1≤in≤In(548)
其中，a(n)in
为向量a(n)中的第in个元素。以一个三阶的秩1张量χ=abc为例，三阶秩1张量结构如图56所示。



图55三阶张量χ∈RI×J×K示意图




图56三阶秩1张量结构



在二维矩阵中，一个大小为I×J矩阵Y的任一元素可以表示为


yi,j=∑Rr=1ai,rbj,r(549)
其中,yi,j为矩阵Y的第(i,j)位置的元素； ai,r为一个二维矩阵A的第(i,r)位置的元素； bj,r为矩阵B的第(j,r)位置的元素。这个二维矩阵Y还可以表示为另外一种形式： 
Y=a1bT1+…+ajbTj+…+aRbTR(550)

其中，aj和bj分别为矩阵A和B的第j列元素组成的向量，其中j=1,2，…,J
。式(550)表示一个二维矩阵Y可以利用R个二元矢量的叉乘相加之和来表示，其中每个二元矢量的叉乘的结构都是一个秩1矩阵。这也是二维矩阵的秩的一种定义方法，即对矩阵Y进行双线性分解需要的最小秩1矩阵的个数R。
类似地，这种秩的定义方法同样可以推广到多维张量上，例如一个三维张量χ∈CI1×I2×…×IN，它其中的任意元素也可以表示为三个二维矩阵A、B和C的元素的乘积，即

xi,j,k=∑Rr=1ai,rbj,rck,r(551)
其中，xi,j,k为三维张量中第（i,j,k）位置的元素； ai,r为一个二维矩阵A的第(i,r)位置的元素； bj,r为矩阵B的第(j,r)位置的元素； ck,r为矩阵C的第(k,r)位置的元素。类似地，分解张量χ所需的最小秩1张量的个数R即为三维张量的秩。
下面首先对两种矩阵的乘积Kronecker积和KhatriRao积做简单的介绍，并给出它们的定义。
两个不同维度的矩阵A∈RI×J和B∈RK×L之间的Kronecker积即为矩阵A中的每一个元素均与矩阵B中的每一个元素相乘的结果，表示为AB。通过Kronecker积可以得到一个大小为(IK)×(JL)的矩阵，如式(552)所示： 

AB=a11Ba12B…a1JB

a21Ba22B…a2JB



aI1BaI2B…aIJB

=［a1b1a1b2a1b3…aJbL-1aJbL］(552)
其中，aj为矩阵A的第j列(j=1,…,J)； bl为矩阵B的第l列(l=1,2，…,L)。
矩阵和矩阵B∈RJ×K的KhatriRao积结果表示为A⊙B，是一个大小为(IJ)×K的矩阵： 

A⊙B=［a1b1a2b2…aKbK］(553)
由式(553)可以看出，KhatriRao积是列向的Kronecker积。如果a和b是向量，那么有ab=a⊙b。

2) CP分解模型
1927年，Hitchcock首次提出了张量分解的概念，直到1970年这个概念以Carroll的正则分解(canonical decomposition，CANDECOMP)和Harshman的平行因子分析 (parallel factors，PARAFAC)的形式引进心理测量学才开始流行起来。将CANDECOMP/PARAFAC分解简称为CP分解。
CP分解把一个高维的张量，分解成多个秩1张量的和，通过这样的分解，可以大大地降低参数的维度。希望可以将指定的三阶张量χ∈RI×J×K写为如下形式： 

χ≈∑Rr=1arbrcr
(554)
其中，R代表张量χ的秩，是一个正整数。对于r=1,2，…,R有ar∈RI，br∈RJ，cr∈RK。因此张量χ中的元素可以写为如下形式，尽可能地使模型残差ei，j，k的数值最小： 

xi,j,k=∑Rr=1ai,rbj,rck,r+ei,j,k,i=1,2，…,I; j=1,2，…,J; k=1,2，…,K(555)

三阶张量CP分解示意图如图57所示。可见，CP分解的形式只是原始张量的一种近似的表示，类似于矩阵的低秩分解，没办法保证完全无失真地复原原来的张量。


图57三阶张量CP分解示意图


因子矩阵指的是秩1分量的向量的组合，即A=［a1a2…aR］，矩阵B和C也类似定义。因此可以将式(554)写为矩阵形式： 

X(1)≈A(C⊙B)T
(556)

X(2)≈B(C⊙A)T
(557)

X(3)≈C(B⊙A)T
(558)

三阶张量χ∈RI×J×K也可以用其各方向的切片来表示： 

Xi≈Bdiag(ai:)CT,i=1,2，…,I

Xj≈Adiag(bj:)CT,j=1,2，…,J

Xk≈Adiag(ck:)BT,k=1,2，…,K(559)
其中，Xi
表示三阶张量水平方向的切片； Xj表示侧向的切片； 相应的Xk表示正向的切片，三阶张量各方向切片示意图如图58所示。


图58三阶张量各方向切片示意图


3) CP分解的唯一性
CP分解的唯一性意味着在不考虑列置换模糊和尺度模糊的条件下，式(554)是唯一可能的秩1张量的和与χ的组合。下面首先要介绍一个重要的概念k秩(Kruskal秩)。
定义： 对于给定的矩阵A∈￡I×R，当且仅当A包含至少r+1个独立的列时，A的秩为rA=Rank(A)=r。如果矩阵A的任意k列独立，则A的k秩kA=k。此时，k=R，或者A包含k+1个独立的列，即kA≤rA≤min(I,R),A。
由定义可以推断出，常用的随机矩阵和范德蒙德矩阵都是不仅满秩，而且满k秩的矩阵。依赖于k秩的概念，针对三阶张量的CP分解，Kruskal给出了满足唯一性的充分条件如下： 
kA+kB+kC≥2R+2
(560)
4) CP分解的计算
Carroll提出的交替最小二乘(Alternating Least Square，ALS)方法是当今CP分解的主流算法。对于一个给定的三阶张量χ∈RI×J×K，希望用最接近χ的R个分量来计算CP分解： 

minχ^‖χ-χ^‖,χ^=∑Rr=1arbrcr
(561)
其中，χ^为待优化的参数，被分解为向量外积的形式。交替最小二乘方法先固定B和C来解A，然后固定A和C来解B。然后固定A和B解出C，继续重复整个过程，直到满足某个收敛准则。
通过固定除一个因子矩阵外的所有矩阵的方式，来将该最小化问题转化为一个简单的线性最小二乘的问题。例如，假设B和C是固定的。根据式(556)可以将上述最小化问题改写为矩阵形式： 

minA^‖X(1)-A^(C⊙B)T‖F
(562)
其中X（1）为含噪矩阵，则矩阵A的最小二乘更新为


A^=X(1)［(C⊙B)T］
(563)
同理，根据式(557)，最小二乘问题为
minB^‖X(2)- B^(C⊙A)T‖F
(564)

矩阵B的最小二乘更新为

B^=X(2)［(C⊙A)T］
(565)

根据式(558)，最小二乘问题为


minC^‖X(3)- C^(B⊙A)T‖F
(566)

矩阵C的最小二乘更新为

C^=X(3)［(B⊙A)T］
(567)

直至算法收敛前不断地重复以上的矩阵更新过程。当满足分解的唯一性的条件时，可以得到关于因子矩阵A、B和C的唯一估计。
2.  参数估计方法
为了使式(531)中的相关矩阵符合式(554)中定义的三线性模型，定义4×M的矩阵Rxy=［R1R2…RM］，其第i列为

RiRwii

ej2πfiτx1(αi)·Rwii

e-j2πfiτy1(βi)·Rwii

ej2πfiτx1(αi)·e-j2πfiτy1(βi)·Rwii
(568)
其中，Rwii为矩阵Rw的第k个对角线元素，Rwii=E|wi|2，i=1,2，…,M。由于窄带目标信号之间互不相关，因此矩阵Rw仅对角线元素不为零。可以定义一个三阶张量χ(N-1)×(N-1)×4，其正向切片Xk=Rk,k={1,2,3,4}，其中元素可以写为如下形式： 

Xk=Rk=Ax1diag［(RTxy)k］AHy1(569)

类似地，相比于一维DOA与载频联合估计CP分解方法，在矩阵Rxy增加了一行自相关项来提升参数估计的准确性。由CP分解的基础理论可知，式(569)中定义的三阶张量χ符合典型的CP分解模型，因此χ中的元素可以表示为

χijk≈∑Mm=1aimbjmckm
(570)
其中，aim表示矩阵Ax1的第i行，第m列的元素； bjm表示矩阵Ay1第j行，第m列的元素； ckm表示定义的矩阵Rxy的第k行，第m列的元素。其中i=j=1,…,N-1； k=1,2,3,4。
由于目标信号的二维DOA与载频参数满足式(51)中的条件，因此阵列流型矩阵Ax1和Ay1为标准的范德蒙德矩阵，满k秩等于M。同时，由式(568)的定义可以看出，矩阵Rxy必有两列以上线性不相关，其k秩大于等于2。因此针对式(569)中的三线性模型，有

kAx1+kAy1+kRxy≥2M+2(571)

该模型满足CP分解的唯一性条件，在不考虑列置换模糊和尺度模糊的条件下，Ax1、Ay1和
Rxy可以从张量χ
中唯一恢复。利用交替最小二乘法对以上三线性模型进行求解，可解得到矩阵Rxy。接下来，通过对矩阵中的元素数值进行计算可得如下两式，为了方便，把它们分别记为u^i和ν^i： 


u^i=ficosθicosφi=c2πd·12∠Rxy(2,i)Rxy(1,i)+∠Rxy(4,i)Rxy(3,i)
(572)

ν^i=fisinθicosφi=c2πd·12∠Rxy(3,i)Rxy(1,i)+∠Rxy(4,i)Rxy(2,i)
(573)

通过平均来提高估计的准确度。同理，对x
轴和z
轴接收信号的采样值重复以上步骤可得矩阵Rxz，从而计算可得u^′i和w^′i： 


u^′i=ficosθicosφi=c2πd·12angleRxz(2,i)Rxz(1,i)+angleRxz(3,i)Rxz(4,i)
(574)

w^′i=fisinφi=c2πd·12angleRxz(3,i)Rxz(1,i)+angleRxz(4,i)Rxz(2,i)
(575)

3.  配对方法
类似地，此时还需要进行配对操作，通过调整u^′i中元素的顺序来实现三维参数配对。令u^=［u^1,u^2,…,u^M］,
u^′=［u^′1,u^′2,…,u^′M］，ν^,ν^′,w^和w^′定义类似。配对寻找一个置换矩阵Ξ满足： 

Ξ^=argminΞ‖u^-u^′Ξ‖s.t.Ξi,j=0,1
(576)

计算配对后的矩阵w^： 

w^= w^′Ξ^(577)
此时就得到了配对成功的特征值矩阵{u^,ν^,w^}，目标信号的载频{fi}Mi=1，方位角{θi}Mi=1和俯仰角{φi}Mi=1参数就可以分别通过式(578)~式(580)计算。

θ^i=arctanν^iu^i
(578)

φ^i=arccosu^2i+ν^2iu^2i+ν^2i+ w^2i
(579)

fi=u^2i+ν^2i+ w^2i(580)

一旦载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1的值确定，即可利用式(544)计算矩阵Ax(矩阵Ay和Az类似)，并通过式(545)和式(526)恢复目标信号。
综上所述，将基于双L形阵列MWC的二维DOA与频率联合估计CP分解算法总结如下。
输入： x轴、y轴和z轴每个传感器接收信号的Q个采样快拍x［k］、y［k］和z［k］。
输出： 载频{fi}Mi=1、方位角{θi}Mi=1、俯仰角{φi}Mi=1和窄带信号频谱{Si(f)}Mi=1。
步骤1： 利用快拍值x［k］和y［k］根据式(531)和式(569)计算三阶张量χ； 
步骤2： 利用交替最小二乘法对三阶张量χ进行求解，可求得矩阵Rxy，并利用Rxy中元素计算一对向量{u^，ν^}； 
步骤3： 利用快拍值x［k］
和y［k］重复步骤1和步骤2，求得矩阵Rxz和向量{u^′,w^′}； 
步骤4： 根据式(576)和式(577)对向量中的元素进行配对，得到{u^,ν^,w^}； 
步骤5： 根据式(578)~式(580)计算载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1； 
步骤6： 利用估计的参数计算阵列流型矩阵Ax、Ay和Az； 
步骤7： 根据式(545)和式(526)计算原窄带目标信号的频谱{Si(f)}Mi=1。

5.2.5仿真实验
本节主要对提出的基于双L型阵列的ESPRIT算法、CP分解算法以及Esmaeil等在双L型阵列结构下提出的基于网格压缩感知的二维DOA和载频联合估计方法进行对比分析。同时探究信噪比、阵元个数、阵元间距等系统参数对方法估计精度的影响。
利用式(581)、式(582)和式(583)中定义的参数作为载频fi、方位角θi和俯仰角φi参数估计准确度的评价指标： 

Ef［%］=1MfNyq∑Mi=1fi-f^i
(581)

Eθ［%］=1M180°∑Mi=1θi-θ^i
(582)

Eφ［%］=1M180°∑Mi=1φi-φ^i
(583)
设置信号个数M=3，奈奎斯特频率fNyq=10GHz，带宽B=140MHz的sinc载波调制信号，如式(584)所示。

un(t)=∑Mi=1EiBsinc(Bt)ej2πfit
(584)
调制载频fi在(0，(fNyq-B)/2)范围内选取，方位角θi在(-90°,90°)范围内选取，俯仰角φi在(0°,90°)范围内选取。设置每轴阵元数N=6，总阵元数即为3N-2=16个。设置阵元间距d=0.03m，每通道快拍数Q=300。设置混频序列每周期有65个点，L=65，混频序列的周期频率fp=1.1B=154MHz。相应的低通滤波器的截止频率为fp/2=77MHz，每通道采样速率为fs=fp。

1.  阵元个数对参数估计精度影响
本节对相同条件下，利用不同的传感器阵元个数来接收信号来进行仿真实验。设置信噪比SNR=10dB
，x轴、y轴和z轴每个轴的正方向均匀分布的传感器阵元个数从临界最小个数N=M+1=4到N=10
递增，步进值为1。实验次数设置为500次，每次实验中对三种方法输入相同的信号，每次实验二维DOA与载频参数在定义范围内随机选取，故噪声随机产生。参数恢复效果利用式(581)~式(583)中定义的参数来评价。图59为不同阵元个数下载频估计效果，图510为不同阵元个数下方位角估计效果，图511为不同阵元个数下俯仰角估计效果。
从图中曲线可以发现，随着传感器阵元个数的增加，ESPRIT方法与CP分解方法的估计效果均越来越好。同时，在信噪比SNR=10dB的条件下，基于CP分解方法的估计效果相比于基于ESPRIT方法的估计效果更好。


图59不同阵元个数下载频估计效果





图510不同阵元个数下方位角估计效果




图511不同阵元个数下俯仰角估计效果


2. 阵元间距对参数估计精度影响
接下来分析相邻传感器阵元之间的间距大小对参数估计精确度的影响。理想状态下阵元间距要满足d≤c/fNyq，将设置的奈奎斯特频率代入，可得d≤3·108/1010=0.03m。设置信噪比SNR=10dB，每通道阵元间距从d=0.01m到d=0.06m递增，步进值为0.01m。实验次数设置为500次，每次实验时，二维DOA与载频参数在定义范围内随机选取，噪声随机产生，每次实验中对两种方法输入相同的信号，参数恢复效果利用式(582)~式(584)中定义的参数来评价。
图512为不同阵元间距下载频估计效果，图513为不同阵元间距下方位角估计效果，图514为不同阵元间距下俯仰角估计效果。由曲线可知，在d=0.03m时，载频、方位角和俯仰角参数估计的效果最好。阵元间距越小，两传感器接收信号之间的相关性越大，阵元之间耦合效果越严重，导致恢复效果越差。而当阵元间距d>0.03m时，会导致阵元接收信号之间的相位差过大，造成空间模糊，参数估计效果也会变差。因此，当相邻传感器阵元间距等于接收信号波长的1/2(即d=0.03m)时，各参数估计的效果最好。


图512不同阵元间距下载频估计效果




图513不同阵元间距下方位角估计效果




图514不同阵元间距下俯仰角估计效果


3.  混频信号频率与信号带宽比值对参数估计精度影响
接下来分析混频序列的周期频率与窄带目标信号的最大带宽的比值fp/B对参数估计精确度的影响。设置fp/B的值从0.5~2.3变化，步进值为0.2。实验次数设置为500次，每次实验中二维DOA与载频参数在定义范围内随机选取，噪声也随机产生。每次实验中对两种方法输入相同的信号，利用式(581)~式(583)中定义的参数作为载频{fi}Mi=1，方位角{θi}Mi=1和俯仰角{φi}Mi=1参数估计准确度的评价指标。
图515为不同fp/B值下载频估计效果，图516为不同fp/B值下方位角估计效果，图517为不同fp/B值下俯仰角估计效果。可以发现，随着fp/B的比值的增大，信号参数的估计效果变得越好。由阵列MWC理论可知，应满足fp/B才能实现阵列MWC的完全重构，此时信号最多会被划分为两个子频带内，搬移到基带后自身频谱不会混叠。当混频序列的周期频率小于信号的最大频带的带宽时，信号所在被划分的频带会分割为两个或以上的子频带。在被分别搬移到基带时，由于带宽大于频域划分的网格，


图515不同fp/B值下载频估计效果





图516不同fp/B值下方位角估计效果





图517不同fp/B值下俯仰角估计效果


因此单个频带的信号自身的频谱会发生混叠，造成估计效果下降。因此，参数估计的效果也取决于fp/B的值，比值越高，估计效果越好。
4.  计算时间对比
本节分析对比两种方法在计算量上的差异，以计算的时间长短为评价指标，设置CP分解算法的循环次数从100~1000次，步进值为100次。实验次数设置为500次，每次实验中二维DOA与载频参数在定义范围内随机选取，噪声也随机产生，每次实验中对两种方法输入相同的信号，不同循环次数下两种方法的计算时间对比如图518所示。可以发现，随着循环次数的增加，CP分解方法所需的计算量越来越大，而ESPRIT方法的计算时间与CP分解方法相比要小得多。


图518不同循环次数下两种方法的计算时间对比


5.  噪声对参数估计精度影响
最后分析系统信噪比大小对各方法重构精度的影响，并将提出的两种方法与基于网格压缩感知的二维DOA和载频联合估计方法进行对比。针对式(584)定义的sinc脉冲调制信号，在信号上叠加加性高斯白噪声，设置信噪比从-10~30dB递增，步进值为5dB，实验次数设置为500次。由于基于网格压缩感知的二维DOA和载频联合估计方法定义估计的方位角范围为θi∈(0°,90°)，为了与该方法的估计效果进行比较，每次实验方位角在(0°,90°)
范围内随机选取，俯仰角在(0°,90°)
范围内随机选取，载频在［0,(fNyq-B)/2］
范围内随机选取，噪声也随机选取，对三种方法输入相同的信号。利用式(581)~式(583)中定义的参数作为载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1参数估计准确度的评价指标。
图519为不同信噪比下载频估计效果，图520为不同信噪比下方位角估计效果，图521为不同信噪比下俯仰角估计效果。随着信号输入信噪比的增加，各方法的参数估计效果均越来越好。在信噪比较低时，本节提出的基于双L型阵列的CP分解方法的参数恢复准确度要高于ESPRIT方法和CS网格搜索算法。相比于CP分解算法，ESPRIT算法仅通过两次奇异值分解即可直接求得参数，而CP分解是一个最优化求解的过程，通过不停地迭代来求解参数，其精度同时也受迭代次数和迭代初值的影响。因此，ESPRIT算法的复杂度和计算量相比于CP分解方法较小，但鲁棒性较差，在信号输入信噪比比较低的情况下表现一般，估计效果欠佳。
基于CS的网格搜索算法也是利用双L型阵列MWC结构进行信号的采集，由于阵列形式相同，因此本节将提出的ESPRIT算法与CP分解方法与其相比较。可以看出，基于CS的网格搜索算法的重构精度取决于网格的大小以及参数是否正好落于网格上，在信噪比比较高的时候，该方法的精确度受到网格大小的限制，因此，相较于本节提出的两种方法估计的误差较大。同时,基于CS的网格搜索算法的方位角只能重构θi∈(0°,90°)的角度，否则会发生角度模糊，提出的两种方法均可以恢复(-90°,90°)之内的角度，因此适用性更好。


图519不同信噪比下载频估计效果





图520不同信噪比下方位角估计效果




图521不同信噪比下俯仰角估计效果



6.  L型阵列与双L型阵列特点比较
最后对基于L型阵列的一维DOA与载频联合估计ESPRIT方法和CP分解方法以及提出的基于双L型阵列的二维DOA与载频联合估计ESPRIT方法和CP分解方法的估计效果以及实用性进行比较与分析。为了保证公平性，设每个轴的阵元个数均为N=6，因此L型阵列的总阵元个数为（2N-1）=11个，双L型阵列的总阵元个数为（3N-2）=16个。针对式(586)定义的sinc脉冲调制信号，在信号上叠加加性高斯白噪声，设置信噪比从-10~30dB递增，步进值为5dB，实验次数设置为500次。每次实验方位角在(-90°,90°)范围内随机选取，俯仰角在(0°,90°)范围内随机选取，载频在［0,(fNyq-B)/2］范围内随机选取，噪声也随机选取，对两种方法输入相同的方位角参数与载频参数。利用式(581)~式(583)中定义的参数作为载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1参数估计准确度的评价指标。
图522为不同信噪比下载频估计效果，图523为不同信噪比下方位角估计效果，图524为不同信噪比下俯仰角估计效果。可以看出，由于双L型阵列需要进行额外的配对操作，因此估计效果比L型阵列结构稍差。但双L型阵列结构阵元较多，在计算时使用的样本较多，因此最终估计效果与L型阵列相差不大。


图522不同信噪比下载频估计效果





图523不同信噪比下方位角估计效果





图524不同信噪比下俯仰角估计效果


但是实际上信号存在于三维空间中，假设所有目标信号均在一个平面内的情况不符合现实中对目标信号到达方向检测任务的要求。显然，二维DOA与载频联合估计的方法比一维DOA与载频联合估计的方法更实用。
5.3基于L型延迟阵列MWC的二维DOA和载频联合估计方法
考虑到基于双L型阵列MWC的二维DOA和频率的联合估计ESPRIT方法和CP分解方法中的配对操作比较复杂，同时重复利用轴阵元接收信号的采样值进行计算会造成冗余。为了省略配对操作，本节进一步提出了基于L型延迟阵列MWC的二维DOA和频率的联合估计CP分解方法和ESPRIT方法。利用延迟通道估计载频，可以直接计算三维参数估计问题，无须额外的配对操作，减小了算法复杂度。下面分别对这两种方法的原理进行推导，并进行仿真验证并对比，探究系统各参数设置对参数估计精度的影响。最后对两种阵列结构进行比较和总结。
5.3.1L型延迟阵列MWC系统
L型延迟阵列结构如图525所示，它由两个相互正交的均匀线性阵列组成，在x轴和y轴的正方向上均有N个以相同间隔均匀分布的传感器阵元{x1,x2,…,xN}和{y1,y2,…,yN}，两个轴的传感器阵列在原点共用同一个传感器。类似地，为了避免空间角度多值模糊的情况，每个轴的相邻阵元之间的间距应满足d≤c/fNyq。


图525L型延迟阵列结构


设有M个互不相关的窄带目标信号si(t)入射到该阵列，每个窄带信号均有一个调制载频为{fi}Mi=1∈R，其方位角θi
和到达角φi与5.2节的定义相同，对于i=1,2，…,M，有θi∈(-90°,90°),φi∈(0°,90°)。为了避免阵列模糊，假设信号满足式(51)。
L型延迟阵列中y轴的每个传感器阵元后连接一个单独的MWC通道结构，由混频器、滤波器和低速采样模块构成。而x轴每个传感器阵元后连接两个MWC通道，其中一个通道在混频模块前增加一个固定的已知的延时模块，x轴延迟通道结构图如图526所示。



图526x轴延迟通道结构图


其中，x1［n］和z1［n］分别表示x轴第一个阵元x1接收信号的非延迟MWC通道与延迟MWC通道的采样值。类似地，定义x［k］、y［k］和z［k］分别为x轴、y轴和x轴延迟通道的观测值； 定义与x轴的夹角为αi
，与y轴的夹角为βi，与z轴的夹角为γi
。可知角度之间的关系如式(52)和式(53)所示。与双L型阵列结构传感器接收的信号不同的是，x轴第n个传感器接收的信号的经延迟之后的信号zn(t)可以表示为

zn(t)≈∑Mi=1si(t)ej2πfi(t+τ)
(585)
其中，τ为设置的已知的延迟时间，满足τ≤1/fNyq。而x轴非延迟通道的接收信号则与式(56)中定义的相同，y轴的接收信号也采用类似的定义方法。
5.3.2基于ESPRIT的参数估计方法
下面首先将ESPRIT方法进行改进，结合L型延迟阵列的特点，省略额外的配对步骤，实现基于ESPRIT算法的二维DOA和载频的联合估计方法。
类似地，将x轴、y轴和z轴延迟通道接收信号的采样值x［k］、y［k］和z［k］分为式(528)中的两个子阵，并利用采样值构造如下的相关矩阵： 

R1=E{x1yH1}=Ax1RwAHy1

R2=E{x2yH1}=Ax2RwAHy1=Ax1ΦxRwAHy1

R3=E{x1yH2}=Ax1RwAHy2=Ax1ΦHyRwAHy1

R4=E{z1yH1}=Az1RwAHy1=Ax1ΦzRwAHy1

R5=E{z2yH1}=Az1RwAHy1=Ax1ΦxΦzRwAHy1

R6=E{z1yH2}=Az1RwAHy1=Ax1ΦzΦHyRwAHy1(586)
其中，旋转矩阵Φx、Φy和Φz中的元素如式(587)所示： 

Φx=diag{exp(j2πf1τx1(α1)),…,exp(j2πfMτx1(αM))}

Φy=diag{exp(j2πf1τy1(β1)),…,exp(j2πfMτy1(βM))}

Φz=diag{exp(j2πf1τ),…,exp(j2πfMτ)}
(587)
其中，τx1(αi)=dccosαi=dccosθicosφi； τy1(βi)=dccosβi=dccosθisinφi。通过观察式(587)可知，由于τ已知，因此从对角阵Φz中可以直接求得载频fi，当三个旋转不变矩阵中的元素对应时，将求出的载频代入矩阵Φx和Φy中即可求得信号的方位角θi和俯仰角φi。由于增加了延迟通道，可以直接利用延迟通道和未延迟通道的采样值估计载频，利用两个轴阵元接收的信号采样值的相关矩阵直接计算三维参数估计问题，无须额外的配对操作，只需计算一次奇异值分解，减小了算法复杂度。
考虑式(584)中构造的相关矩阵，将采样值联合，可以直接构造如下矩阵： 
R=［R1; R2; R3; R4; R5; R6］(588)
对R进行奇异值分解，U为前M个非零奇异值对应的左奇异向量。存在一个M×M的可逆矩阵T使得式(589)成立： 

U=U1

U2

U3

U4

U5

U6=Ax1

Ax1Φx

Ax1ΦHy

Ax1Φz

Ax1ΦxΦz

Ax1ΦzΦHyT
(589)
其中，U是一个6(N-1)×M
的矩阵，Ui为(N-1)×M的矩阵，i=1,2，…,6。从式(587)，可以推导出矩阵Ui和旋转不变阵Φx、Φy和Φz之间的关系。为了书写简洁，定义如式(535)的矩阵V1、V2
和V3： 

V1=U1

U4U2

U5=TΦxT-1
(590)

V2=U1

U4U3

U6=TΦHyT-1
(591)

V3=U1

U2

U3U4

U5

U6=TΦzT-1(592)

为了保证旋转不变矩阵中的元素顺序的对应，对矩阵(V1+V2+V3)进行特征值分解得到对应的特征向量矩阵T^
，此时再利用同一个特征向量T^
对三个矩阵依次进行求解即可使其中对角线元素一一对应。利用式(591)可计算得到三个特征值的顺序一一对应的旋转不变矩阵Φ^x
、Φ^y
和Φ^z： 

Φ^x= T^-1V1T^


Φ^y=(T^-1V2T^)H

Φ^z= T^-1V3T^
(593)
令ui、vi和wi分别为Φx
、Φy
和Φz的第i个对角线元素，i=1,2,…,M，则有

2πfidccosθicosφi=angle(ui)

2πfidcsinθicosφi=angle(vi)

2πfiτ=angle(wi)
(594)
得到顺序一一对应的对角线元素后，目标信号的载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1参数就可以通过式(595)、式(596)和式(597)计算： 


fi=angle(wi)2πτ
(595)

θi=arctanangle(vi)angle(ui)
(596)

φi=arccoscangle2(ui)+angle2(vi)2πd
(597)

一旦载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1的值确定，即可利用估计的值计算矩阵Ax、Ay和Az。其中Az由式(598)计算： 


Az=111

ej2πf1τ…ej2πfMdcτ



ej2πf1d(N-1)cτ…ej2πfMd(N-1)cτ
(598)
得到阵列流型矩阵后，通过式(545)和式(526)即可恢复窄带目标信号。
综上所述，将基于L型延迟阵列MWC的二维DOA与频率联合估计ESPRIT算法总结如下。
输入： x轴、y轴和x轴延迟后的每个通道的Q个采样快拍x［k］、y［k］和z［k］。
输出： 载频{fi}Mi=1、方位角{θi}Mi=1、俯仰角{φi}Mi=1和窄带信号频谱{Si(f)}Mi=1。
步骤1： 利用快拍值x［k］、y［k］和z［k］和根据式(586)和式(588)计算协方差矩阵R； 
步骤2： 对R进行奇异值分解，得到左奇异向量U，并根据式(535)和式(536)计算矩阵V1和V2； 
步骤3： 对V1+V2+V3进行特征值分解，求得旋转不变矩阵{Φ^x,Φ^y,Φ^z}； 
步骤4： 根据式(595)~式(597)计算载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1； 
步骤5： 利用估计的参数计算阵列流型矩阵Ax、Ay和Az； 
步骤6： 根据式(545)和式(526)计算原窄带目标信号的频谱{Si(f)}Mi=1。

5.3.3基于CP分解的参数估计方法
类似地，将CP分解方法也进行了相应的扩展，应用于L型延迟阵列结构，通过一次分解过程，无须配对操作即可实现载频和二维到达角的估计。
为了使式(588)中的相关矩阵满足CP分解的模型，定义了一个大小为6×M的矩阵R=［R1R2…RM］，其第i列元素为

RiΔRwii

ej2πfiτx1(αi)·Rwii

e-j2πfiτy1(βi)·Rwii

ej2πfiτ·Rwii

ej2πfiτx1(αi)·ej2πfiτ·Rwii

ej2πfiτ·e-j2πfiτy1(βi)·Rwii(599)

由于窄带信号互不相关，因此矩阵Rw除对角线元素外的元素均为零，可以定义一个三阶张量χ(N-1)×(N-1)×6
，其正向切片Xk=Rk,k=1,2，…,6，Rk可以表示为阵列流行矩阵和矩阵R的乘积的形式：  

Xk=Rk=Ax1diag［(RT)k］AHy1
(5100)

该模型满足CP分解的唯一性条件，在不考虑列置换模糊和尺度模糊的条件下，矩阵Ax1、Ay1和Rxy可以从张量χ中唯一恢复。利用交替最小二乘法对以上三线性模型进行求解，可解得到矩阵R，通过对矩阵中的数值进行计算可得式(5101)~式(5103)，为了方便，把它们分别记为u^i、v^i和w^i： 

u^i=ficosθicosφi=c2πd·12∠R(2,i)R(1,i)+∠R(5,i)R(4,i)(5101)

v^i=fisinθicosφi=c2πd·12∠R(3,i)R(1,i)+∠R(6,i)R(3,i)(5102)

w^i=fiτ=12π·13∠R(4,i)R(1,i)+∠R(5,i)R(2,i)+∠R(6,i)R(3,i)
(5103)

通过平均来提高估计参数的准确性。载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1就可以通过式(5104)~式(5106)计算。


fi=w^iτ
(5104)

θ^i=arctanv^iu^i
(5105)

φ^i=arccosτu^2i+ v^2iw^i
(5106)
一旦载频、方位角和俯仰角的值确定，即可利用式(544)计算矩阵Ax、矩阵Ay和Az类似，并通过式(545)和式(526)恢复目标信号。
综上所述，将基于L型延迟阵列MWC的二维DOA与频率联合估计CP分解算法总结如下。
输入： x轴、y轴和x轴延迟后的每个通道的Q个采样快拍x［k］、y［k］和z［k］。
输出： 载频{fi}Mi=1、方位角{θi}Mi=1、俯仰角{φi}Mi=1和窄带信号频谱{Si(f)}Mi=1。
步骤1： 利用快拍值x［k］
、y［k］和z［k］根据式(599)和式(5100)计算三阶张量χ； 
步骤2： 利用交替最小二乘法对三阶张量χ进行求解，可求得矩阵R，并利用R中元素计算一组向量{u^，v^,w^}； 
步骤3： 根据式(5104)~式(5106)计算载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1； 
步骤4： 利用估计的参数计算阵列流型矩阵Ax、Ay和Az； 
步骤5： 根据式(545)和式(526)计算原窄带目标信号的频谱{Si(f)}Mi=1。

5.3.4仿真实验
本节主要对基于L型延迟阵列的ESPRIT算法、CP分解算法进行对比分析； 同时探究信噪比、通道数等系统参数对方法估计精度的影响。
接下来对在L型延迟阵列下基于ESPRIT的二维DOA与载频联合估计方法进行仿真实验验证。本节中进行无噪声仿真实验，参数设置与5.2.5节中相同。设置每轴阵元数N=6，总阵元数为2N-1=11个，而总的通道数为3N-1=17个。

1.  阵元个数对参数估计精度影响
接下来分析每坐标轴阵元个数对参数估计精确度的影响。设置信噪比SNR=10dB，x轴、y轴和z轴的每个轴正方向均匀分布的传感器阵元个数从临界最小个数N=M+1=4到N=10递增，步进值为1，其他实验参数设置与5.2.5节中相同。实验次数设置为500次，每次实验中对3种方法输入相同的信号，每次实验二维DOA与载频参数在定义范围内随机选取，噪声随机产生。参数恢复效果利用式(581)~式(583)中定义的参数来评价。
图527为不同阵元个数下载频估计效果，图528为不同阵元个数下方位角估计效果，图529为不同阵元个数下俯仰角估计效果。可以发现，随着传感器阵元个数的增加，ESPRIT方法与CP分解方法的估计效果都越来越好。


图527不同阵元个数下载频估计效果





图528不同阵元个数下方位角估计效果





图529不同阵元个数下俯仰角估计效果


2. 阵元间距对参数估计精度影响
接下来分析阵元之间间距大小对参数估计精确度的影响。设置信噪比SNR=10dB，每通道阵元间距从d=0.01m到d=0.06m递增，步进值为0.01m，本小节其他实验参数设置与5.2.5节相同。实验次数设置为500次，每次实验中二维DOA与载频参数在定义范围内随机选取，噪声随机产生，每次实验中对两种方法输入相同的信号，参数恢复效果利用式(581)~式(583)中定义的参数来评价。
图530为不同阵元间距下载频估计效果，图531为不同阵元间距下方位角估计效果，图532为不同阵元间距下俯仰角估计效果。由曲线可知，在d=0.03m时，载频、方位角和俯仰角参数估计的效果最好。阵元间距越小或越大时，都会导致参数估计效果变差。


图530不同阵元间距下载频估计效果




图531不同阵元间距下方位角估计效果




图532不同阵元间距下俯仰角估计效果



3.  混频信号频率与信号带宽比值对参数估计精度影响
接下来分析混频序列的周期频率与窄带目标信号的最大带宽的相对大小对参数估计精确度的影响。设置fp/B的值从0.5~2.3变化，步进值为0.2。设置信噪比SNR=10dB，实验次数设置为500次，每次实验中二维DOA与载频参数在定义范围内随机选取，利用式(581)~式(583)中定义的参数作为载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1参数估计准确度的评价指标。
图533为不同fp/B值下载频估计效果，图534为不同fp/B值下方位角估计效果，图535为不同fp/B值下俯仰角估计效果。可以看出，随着fp/B值的增大，信号参数的估计效果越好。


图533不同fp/B值下载频估计效果





图534不同fp/B值下方位角估计效果





图535不同fp/B值下俯仰角估计效果



4.  计算时间对比
本节中分析对比两种方法在计算量上的差异。以计算的时间长短为评价指标，设置CP分解算法的循环次数从100~1000次，步进值为100次。实验次数设置为500次，每次实验中二维DOA与载频参数在定义范围内随机选取，噪声也随机产生，每次实验中对两种方法输入相同的信号，不同循环次数下两种方法的计算时间对比如图536所示。可以发现，随着循环次数的增加，CP分解方法所需的计算量越来越大，而ESPRIT方法的计算时间与CP分解方法相比相当小。


图536不同循环次数下两种方法的计算时间对比



5.  噪声对参数估计精度影响
最后分析系统信噪比大小对各方法参数估计精度的影响，并将提出的两种方法与在双L型阵列基础上的基于网格压缩感知的二维DOA和载频联合估计方法进行对比，虽然阵列结构不相同，但有一定的参考价值。设置信噪比从-10~30dB递增，步进值为5dB，实验次数设置为500次。每次实验方位角在(0°,90°)范围内随机选取，俯仰角在(0°,90°)
范围内随机选取，载频在［0,(fNyq-B)/2］
范围内随机选取，噪声也随机选取，对3种方法输入相同的信号。其他实验参数设置与5.2.5节相同。利用式(581)~式(583)中定义的参数作为载频{fi}Mi=1、方位角{θi}Mi=1和俯仰角{φi}Mi=1参数估计准确度的评价指标。
图537为不同信噪比下载频估计效果，图538为不同信噪比下方位角估计效果，图539为不同信噪比下俯仰角估计效果。可以发现，随着信号输入信噪比的增加，各方法的参数估计效果均越来越好。类似于双L型阵列中两种方法噪声对比试验的结果，在有噪声时，本节提出的基于L型延迟阵列的CP分解方法的参数恢复准确度要高于ESPRIT方法，但由于使用的传感器阵元数较少，因此比CS网格搜索算法稍差。而基于CS的网格搜


图537不同信噪比下载频估计效果






图538不同信噪比下方位角估计效果





图539不同信噪比下俯仰角估计效果


索算法在信噪比比较高的时候，该方法的精确度受到网格大小的限制，因此相较于本节提出的两种方法估计的误差较大。同时，基于CS的网格搜索算法的方位角只能重构θi∈(0°,90°)内的角度，提出的两种方法均可以恢复在(-90°,90°)范围的角度，因此适用性更好。

5.3.5双L型阵列与L型延迟阵列结构仿真对比
本节将提出的在双L型阵列结构和L型延迟阵列下的同一种方法的参数估计效果进行对比。实验参数设置与5.2.5节相同，信噪比从-5~30dB递增，步进值为10dB，实验次数设置为500次。每次实验中，二维DOA与载频参数在定义范围内随机选取，参数恢复效果利用式(581)~式(583)中定义的参数来评价。
首先针对ESPRIT方法进行实验。图540为不同信噪比下的基于两种阵列结构的ESPRIT方法载频估计效果，图541为不同信噪比下的基于两种阵列结构的ESPRIT方法方位角估计效果，图542为不同信噪比下的基于两种阵列结构的ESPRIT方法到达角估计效果。



图540不同信噪比下的基于两种阵列结构的ESPRIT方法载频估计效果





图541不同信噪比下的基于两种阵列结构的ESPRIT方法方位角估计效果





图542不同信噪比下的基于两种阵列结构的ESPRIT方法到达角估计效果


对于方位角的估计来说，由于两种阵列结构中均利用x轴和y轴的采样值直接计算方位角，因此两种阵列结构下的估计效果很接近。而对于载频和俯仰角的估计来说，由于双L型阵列中利用三轴的不同传感器接收的采样值共同估计俯仰角和载频，而L型延迟阵列只有x轴和y轴有不同的传感器接收信号，因此基于L型延迟阵列的方法俯仰角和载频参数的估计效果比双L型阵列的稍差。

对于基于双L型阵列结构和L型延迟阵列的CP分解方法，做了类似的仿真验证，结果与ESPRIT方法的规律类似。图543为不同信噪比下的基于两种阵列结构的CP分解方法载频估计效果，图544为不同信噪比下的基于两种阵列结构的CP分解方法方位角估计效果，图545为不同信噪比下的基于两种阵列结构的CP分解方法俯仰角估计效果。由于使用的传感器阵元数相较双L型阵列较少，因此方位角的估计性能较差。


图543不同信噪比下的基于两种阵列结构的CP分解方法载频估计效果





图544不同信噪比下的基于两种阵列结构的CP分解方法方位角估计效果





图545不同信噪比下的基于两种阵列结构的CP分解方法俯仰角估计效果



接着对这几种方法的运算量进行对比，由前面的理论推导可知，双L型阵列下的CP分解方法需要进行两次交替最小二乘法的计算，而L型延迟阵列的CP分解方法只需要进行一次交替最小二乘的计算，因此基于L型延迟阵列的方法计算量较小。为了对比两种阵列下的两种方法在计算量上的差异。以计算的时间长短为评价指标，设置CP分解算法的循环次数从100~1000次，步进值为100次。实验次数设置为500次，每次实验中二维DOA与载频参数在定义范围内随机选取，噪声也随机产生，每次实验中对4种方法输入相同的信号，不同循环次数下4种方法的计算时间对比如图546所示。可以发现，随着循环次数的增加，两种CP分解方法所需的计算量都越来越大。基于L型延迟阵列的CP分解方法计算量要明显小于基于双L型阵列结构的CP分解算法，计算所需时间几乎为双L型延迟阵列的一半。两种ESPRIT方法的计算时间与CP分解方法相比相当小，相互之间也相差不大。


图546不同循环次数下4种方法的计算时间对比


由此可见，虽然L型延迟阵列结构下的信号载频与二维DOA参数估计方法省略了配对步骤，降低了算法复杂度和计算量，但由于使用的传感器阵元数相较双L型阵列较少，因此方位角和载频的估计性能较差。作为总结，表51列出了本章所提出方法的特点。


表51本章所提出方法的特点



方法
算法
传感器个数
通道数
配对
优劣

双LMWC

ESPRIT
CP分解
3N-2

3N-2

需要额外

配对

计算量更小，受噪声影响大； 

鲁棒性更好，效果受迭代次数和初值影响

L延迟MWC
ESPRIT
CP分解
2N-1

3N-1

无须额外

配对

计算量更小，受噪声影响大，载频与俯仰角估计效果较差； 

计算量相对更小，鲁棒性更好，载频与俯仰角估计效果较差

本章小结
本章在阵列MWC理论的基础上进行扩展，研究基于压缩采样的空频域参数联合估计方法。提出了基于双L型阵列和L型延迟阵列的ESPRIT及CP分解方法，可以估计窄带目标信号的载频和二维DOA参数，同时重构目标信号的原始波形。最后通过仿真实验对基于双L型阵列和L型延迟阵列的ESPRIT及CP分解方法的原理进行了验证，并对ESPRIT及CP分解方法进行了对比，探究了系统参数的影响。