基于稀疏协方差矩阵迭代的单快拍气流速度估计算法

虞 飞陶建武陈 诚钱立林

①(海军航空工程学院控制工程系 烟台 264001)

②(空军航空大学飞行器控制系 长春 130022)

③(94638部队 南昌 330201)

基于稀疏协方差矩阵迭代的单快拍气流速度估计算法

虞 飞*①②陶建武②陈 诚③钱立林②

①(海军航空工程学院控制工程系 烟台 264001)

②(空军航空大学飞行器控制系 长春 130022)

③(94638部队 南昌 330201)

该文研究基于声传感器阵列的单快拍气流速度估计问题。首先，根据声波在亚音速和超音速气流中的传播特性，针对特定的测量装置，建立了声传感器线性阵列的输出模型。在此基础上，提出一种稀疏协方差矩阵迭代的单快拍(Sparse Covariance Matrix Iteration with a Single Snapshot, SCMISS)气流速度估计算法，与其他稀疏估计方法相比，该文提出的SCMISS算法无需正则化参数选择，计算量更低，具有更强的实时性，且只需单快拍采样数据就可对亚音速和超音速气流速度进行统一估计。最后，为了评价所提算法的估计性能，推导了气流速度估计的克拉美-罗界(Cramér-Rao Bound, CRB)表达式。仿真实验验证了该算法的有效性。

阵列信号处理；气流速度估计；单快拍；稀疏参数估计；声传感器阵列；克拉美-罗界

1 引言

当前，对高速、高机动性飞行器的设计研究引起了世界各主要军事强国的极大关注，而传统的以空速管为代表的空速测量方法无法满足这类先进飞行器的设计要求。现有的空速测量方法大都可以测量飞行器的亚音速空速，但是在超音速飞行状态下，空气流动特性发生质的变化而产生激波，使得目前的测量方法要么难以适用，如空速管，要么需要通过参数校正才能适用，如嵌入式大气数据传感(Flush Air Data Sensing, FADS)系统[1]。因此，有必要建立统一的亚音速和超音速空速测量模型，使模型在各个飞行状态都能输出正确的空速值，同时避免因亚音速和超音速飞行状态的气流差异带来的高成本的参数校正，但是，如何建立统一的跨音速测量模型是当前大气数据测量研究的一个难点。

声矢量传感器作为一种新型传感器，可以同步测量声场中同一点处的声波质点振速和声压[2]。文献[3]采用声矢量传感器构成了一种新型的气流速度测量方法，它可以利用声矢量传感器的测量值计算气流速度，再根据测量装置所在的机体表面形状及飞行姿态进行解算，可由气流速度成功导出飞行器的空速。它的缺点在于仅采用单个声矢量传感器来测量气流速度，使得系统在抗干扰能力和测量精度上难以满足实际要求。文献[4]则建立了声矢量传感器阵列输出数据与气流速度之间的数学模型，并研究了一种数据缺失情形下的气流速度估计算法，此算法可以有效提高气流速度的估计精度。然而，前面提出的气流速度测量模型和算法都只适用于亚音速气流速度的估计，难以满足超音速飞行器的使用要求。针对超音速气流速度测量问题，文献[5]提出了一种声波在超音速气流中传播的分析方法——等效声源法，并由此近似求得马赫锥内一点的质点振速表达式，建立了基于声矢量传感器的信号输出模型，但是该模型只适用于超音速气流速度测量。为此，本文基于声传感器，提出了一种能同时适用于亚音速和超音速气流速度测量的新方法。该方法综合考虑了亚音速和超音速两种状态下气流流动的差异，建立了统一的跨音速测量模型，从原理上，实现了从亚音速到高超音速整个范围内气流速度的统一测量。

前面提到的气流速度估计算法都是基于多快拍采样数据的方法，而在实际应用中，由于飞行器作复杂机动飞行动作，以及气流环境本身复杂多变等因素的影响，使得待测气流速度可能是实时变化的。这时，只有少数快拍甚至单快拍数据可以利用。虽然文献[6]提出的基于鲁棒H∞滤波的气流速度估计算法只需单快拍采样数据，但该算法无法用于超音速气流速度测量。本文在建立统一的跨音速测量模型的基础上，提出了一种基于稀疏协方差矩阵迭代的单快拍气流速度估计算法。此算法只需单快拍采样数据就可以对亚音速和超音速气流进行统一测量，这样不仅节省了数据存储空间，而且能实现对快变气流速度的跟踪测量。

2 阵列接收模型

气流速度测量装置如图1所示，管路内径为D，一般可要求D在1～2 cm范围内，整个管路测量装置可以嵌入式安装在飞机机头或机翼前缘。在管路内剖面建立xOy直角坐标系。在坐标原点O处安装一个声波发射换能器，在管路内壁正上方平行于x轴方向布置L个声波接收换能器，分别编号1,2,…, L，设第l(l=1,2,…,L )个接收换能器相对于发射换能器的方位角为θl，则θl∈(-π/2,π/2)(从y轴测量)。

如图1所示，设声波波阵面由发射换能器到达第l个接收换能器的传播距离为Rl，则Rl=D/ cosθl。假设来流从x轴正向吹来。当气流速度v=0时，声波波阵面由发射换能器发射到达第l个接收换能器的传播时间为

图1 测量装置剖面图

其中，c为声波传播速度。当0vc＜＜时，对应亚音速情形，声波波阵面传播时间变为

当vc≥时，对应超音速情形。这时，波阵面一方面扩展，一方面顺流而下，所有波阵面被挤压而聚集在一个被称为马赫锥的锥面内。马赫锥的半顶角与气流速度之间有如下关系：

式中α称为马赫角。只有当接收换能器位于y轴的左半平面且方位角满足π/2+θl＜α时，换能器才能接收到声波信号，且声波波阵面由发射换能器传播到该接收换能器的时间为

由π/2+θl＜α得cosθl＜c/v，则有v＜c/ cosθl。这个不等式给出了方位角为θl的接收换能器可测量到的超音速气流速度的上界。当待测气流速度v≥c/cosθl时，该接收换能器将接收不到发射的声波。

以下在建立接收换能器基阵的阵列接收模型时，假设在速度为v的气流作用下，L个接收换能器都能接收到发射换能器发射的声波信号。考虑发射换能器发射一个连续单频声波信号()st= Aej(ωt+φ(t))，则第l(l=1,2,…,L)个接收换能器的接收信号可以表示为

式中，A为发射换能器发射的声波信号的幅值，ω=2πf , f为声波信号的频率，φ(t)为在区间[0,2π]之间均匀分布的随机相位。σs(t)为零均值，单位方差复高斯随机过程，“”表示取复数模值。nl(t)为第l个接收换能器上的零均值，方差为σl的加性复高斯白噪声。假设信号与噪声互不相关，以x1(t)为参照，将式(5)变换得

式中，s1(t)表示第1个接收阵元接收的声波信号。Dl1=τl-τ1为声波到达接收换能器1和l的延时差值，其中包含了待估计的气流速度信息v,R1/Rl为相应接收信号的幅度比值。将各个接收换能器的接收信号联立起来，可得到接收换能器基阵的阵列输出模型。

式中，x(t)=[x(t),x(t),…,x(t )]T,a(v)=[1,(R/

12L1为接收换能器基阵的L×1维阵列流形矢量，为加性噪声矢量。

3 稀疏协方差矩阵迭代的单快拍气流速度估计算法

3.1 协方差矩阵的稀疏表示

将待处理的气流速度范围V进行N次均匀网格划分{v1,v2,…,vN}，假设网格划分密度足够大，以至于真实的气流速度v始终可以落在其中某一网格上，则可得如下稀疏表示模型

式中，Φ=[a(v1),a(v2),…,a(vN)]为L×N维稀疏基矩阵，a(vn)表示第n个网格对应的阵列流形矢量。α(t)为N×1维稀疏信号矢量。应注意的是，由于在亚音速气流作用下y轴左、右两端的接收换能器都可以接收到发射的声波信号，而在超音速情形下，只有y轴左半平面的接收换能器才可以接收到声波，因此，两种情形下气流速度范围的网格划分应分开考虑。在亚音速时，只需对整个亚音速域[0,)c进行均匀网格划分；在超音速时，只需对超音速域[c,c/cosθl),(θl＜0)进行均匀网格划分，而无需对整个超音速域进行划分，这里l表示在待测超音速气流作用下可用的接收阵元个数。

一旦求出了α(t)的最稀疏解ˆα(t)，则αˆ()的t归一化稀疏谱的谱峰对齐的速度网格就是当前气流速度的估计值vˆ。目前，求解()tα的最稀疏解使用最广泛的方法是1-范数最小化稀疏恢复方法。该方法将1-范数最小化问题转化为二阶锥规划问题的一般形式，并在二阶锥规划框架下求解。但在实际测量气流速度时，速度网格划分的数量级将达到310，这将明显降低二阶锥规划求解速度。同时，在很多情况下对正则化参数β的准确选择是很困难的。

下面提出的稀疏协方差矩阵迭代的单快拍气流速度估计算法将极大地解决上述这些问题。设n(t)内各个分量互不相关，则E[n(t)nH(t)]=diag[σ1,…,σL]。假设稀疏信号矢量α(t)与n(t)不相关，且α(t)内各个分量也互不相关，则有阵列接收数据的稀疏协方差矩阵

式中，pn=E[αn(t)(t )]表示α(t)中第n个分量的信号功率，IL为L×L维单位阵。

于是

本文提出的稀疏协方差矩阵迭代的气流速度估计算法就是利用单快拍阵列接收数据，通过迭代算法来估计功率矩阵P的对角元素，中最大元素的下标即为真实气流速度v在网格中的位置下标，由此来实现气流速度的估计。

3.2 单快拍气流速度估计算法

对于单快拍阵列接收数据()tx，有如下协方差矩阵拟合准则[7]

其中，tr(·)表示求矩阵的迹，且由式(12)可得

同时考虑到协方差矩阵拟合准则式(13)中目标函数f的自变量是R,ˆR可作为常量处理，则目标函数f的最小化可以等效为式(16)所示的目标函数g的最小化问题：

其证明过程类似于文献[8]，本文不再赘述。为了节省计算量，这里参照文献[8]转而求解一个与优化问题式(17)相关的优化问题来间接求优化问题式(17)的最优解，则归纳出SCMISS算法的更新公式：

本文算法与文献[7]中所述方法的主要区别有两点：(1)在应用场合方面，本文算法可以应用到部分测量数据缺失的场合；而文献[7]算法仅适用于具有完整测量数据的场合。对于不同的待估计气流速度，有效的接收阵元个数会有所变化。因此，有效的测量数据是随着气流速度动态变化的。如果将文献[7]中的算法直接用于本文的气流速度估计，会因有些测量数据的缺失而导致估计性能严重下降甚至算法失效。而本文算法解决了在测量数据缺失的场合下参数估计问题。(2)在网格划分方面，本文采用了动态网格划分方法，而文献[7]采用的是静态网格划分方法。本文算法在气流速度范围进行网格划分时，充分考虑了亚音速和超音速阶段的气流流动差异，采用了动态网格划分方法。这样既保证了算法的估计精度，又有效降低算法的计算复杂度。

4 算法性能分析

4.1 CRB推导

由给出的阵列输出模型式(7)可知，观测信号x(t )符合复高斯随机分布，记为x(t )～CN(μ,Rx)，其中Rx和μ分别表示观测矢量x(t)的L×L维协方差矩阵和L×1维均值矢量，且μ=0, Rx=(A2/定义未知的实值参数矢量，其中，则关于参数矢量η的Fisher信息矩阵的第(i,j)个元素

式中，ηi表示矢量η的第i个元素。在实际应用中，一般只关注气流速度v，其他参数为多余参数，且气流速度v与其他参数不是互耦的，则关于气流速度v的Fisher信息函数为可表示为[9]

故相对于气流速度v的CRB为

4.2 计算量分析与比较

下面分析SCMISS算法的计算量，并与l1-范数最小化方法作比较。l1-范数最小化方法一般是在二阶锥规划的框架下采用内点法来求最优解，其计算复杂度为O(N3)[10]。设MDN表示复数乘和除的次数，由于SCMISS算法计算量主要集中于功率迭代更新计算中，因此通过分析迭代更新的计算量可以得到该算法的计算复杂度。在每次迭代的第1步中：

计算R(i)共需要η1=2L2(N+L)次MDN；计算R-1(i)需要2L3/3次MDN，计算需要L2次MDN，计算需要2L2次MDN，因此计算共需要次 MDN；计算N+L共需要η3=3(N+L)次MDN。因此，SCMISS算法每次迭代的第1步所需的MDN次数近似为

每次迭代的2步所需的MDN次数为

设算法达到要求的估计精度所需迭代次数为cn，则SCMISS算法中迭代更新计算所需的MDN总次数大致为

由于接收换能器基阵阵元数L数量级一般为10，而速度网格划分数N的数量级高达310，同时SCMISS算法一般经过不超过15次迭代即可得到很高精度的全局最优解，因此，SCMISS算法的计算复杂度明显低于1-范数最小化方法，从而具有更好的实时性。

5 仿真实验

考虑到管路测量装置是拟用于机载嵌入式大气数据传感系统，设管路内径D=1.5 cm ，静止声场中的声速c=340 m/s，发射换能器发射的声波信号频率取6.8 kHz。SCMISS算法和L1-Min算法速度网格划分间隔均取1 m/s。在如下各组实验中，L1-Min算法正则化参数的选择依据是：选择正则化参数β，使得以一个较小的概率成立，从而使残留量与噪声范数的期望尽量匹配，具体参考文献[10]。

实验1 亚音速情形下的算法统计性能实验。在均匀稳定气流中发射换能器发射一单频声波信号入射到由10个接收换能器构成的基阵，各阵元方位角分别为θ=[-88°,-80°,-70°,-60°,-50°,50°,60°, 70°,80°,88°]T。亚音速气流速度v=240 m/s，对单快拍阵列观测数据分别采用SCMISS算法、L1-Min算法和鲁棒H∞滤波气流速度估计算法[6]进行100次Monte Carlo仿真实验，得到气流速度估计vˆ的均方根误差(Root-Mean-Square Error, RMSE)随信噪比的变化曲线，如图2所示。由图2可以看出，本文提出的SCMISS算法的估计精度明显高于另外两种算法，且当信噪比高于2.5 dB时，算法对气流速度估计误差的方差接近CRB，说明此时估计误差接近于最小值。

实验2 超音速情形下的算法统计性能。取7个接收换能器构成的基阵，且各阵元的方位角为θ=[-88°,-87°,-86°,-85°,-82°,-80°,-78°]T。设超音速气流速度v=500 m/s，分别采用SCMISS算法和L1-Min算法进行50次Monte Carlo仿真实验，得到气流速度估计vˆ的RMSE随信噪比的变化曲线，如图3所示。由图3可以看出，采用这两种算法估计的气流速度vˆ的RMSE随着信噪比的增大，与CRB越来越接近，说明算法对气流速度的估计精度越来越高。其中，本文提出的SCMISS算法的估计精度要高于L1-Min算法。

实验3 容错性实验。假设接收换能器基阵的第2个阵元失效，则x2(t)=0。亚音速气流速度v= 240 m/s, SNR=7.5 dB，图4(a)给出了气流速度估计的归一化稀疏功率谱图。图4(b)为无阵元失效时的稀疏功率谱图。图4(c)给出了在速度为v= 1435 m/s的超音速气流作用下，气流速度估计的稀疏功率谱图。从图4可以发现，当接收换能器基阵的某个阵元失效时，本文算法仍可以对亚音速气流速度进行有效估计，但无法正确估计出超音速气流速度。

6 结论

本文建立了稳定气流作用下声传感器阵列的输出模型，提出了一种基于稀疏协方差矩阵迭代的单快拍气流速度估计算法(SCMISS)，分析了SCMISS算法的计算量，并与L1-Min算法的计算量进行了比较，推导了气流速度估计的CRB表达式。理论分析和仿真实验表明：本文提出的SCMISS算法可同时适用于亚音速和超音速气流速度的高精度估计；当某个接收阵元失效时，算法仍可以对亚音速气流速度进行有效估计；SCMISS算法比L1-Min算法具有更低的计算量。

图2 气流速度估计的RMSE随信噪比变化曲线

图3 气流速度估计的RMSE随信噪比变化曲线

图4 气流速度估计的稀疏功率谱图

[1] Whitmore S A, Cobleigh B R, and Haering E A. Design and calibration of the X-33 flush airdata sensing (FADS) system [R]. NASA/TM-1998-206540, 1998.

[2] Wu Y I, Wong K T, and Lau S K. The acoustic vectorsensor’s near-field array-manifold[J]. IEEE Transactions on Signal Processing, 2010, 58(7): 3946-3951.

[3] 王立新, 陶建武. 一种新型的大气数据测量方法[J]. 计量学报, 2011, 32(1): 31-35.

Wang Li-xin and Tao Jian-wu. A new measuring airdata algorithm[J]. Acta Metrologica Sinica, 2011, 32(1): 31-35.

[4] 陈诚, 陶建武, 曾宾. 数据缺失情形下的基于声矢量传感器阵列的空气流动速度估计算法[J]. 电子学报, 2014, 42(3): 491-497.

Chen Cheng, Tao Jian-wu, and Zeng Bin. Estimation of airspeed in case of missing data[J]. Acta Electronica Sinica, 2014, 42(3): 491-497.

[5] Zeng Bin, Tao Jian-wu, Yu Fei, et al.. A new estimation method of airspeed based on acoustic vector sensor[C]. 2013 3th Intermational Conference on Signal Processing, Communications and Computing, Kunming, 2013: 307-310.

[6] 陈诚, 陶建武. 基于声矢量传感器阵列的鲁棒H∞空气速度估计算法[J]. 航空学报, 2013, 34(2): 361-370.

Chen Cheng and Tao Jian-wu. Robust H∞estimation of airspeed based on acoustic vector sensor array[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(2): 361-370.

[7] Stoica P, Babu P, and Li J. New method of sparse parameter estimation in separable models and its use for spectral analysis of irregularly sampled data[J]. IEEE Transactions on Signal Processing, 2011, 59(1): 35-47.

[8] Stoica P, Babu P, and Li J. SPICE: a sparse covariance-based estimation method for array processing[J]. IEEE Transactions on Signal Processing, 2011, 59(2): 629-638.

[9] Korso M N El, Boyer R, Renaux A, et al.. Conditional and unconditional Cramer-Rao bounds for near-field source localization[J]. IEEE Transactions on Signal Processing, 2010, 58(5): 2901-2907.

[10] Malioutov D. A sparse signal reconstruction perspective for source localization with sensor arrays[D]. [Master dissertation], Massachusetts Institute of Technology, 2003.

虞 飞： 男，1987年生，博士生，研究方向为传感器与智能测量系统、阵列信号处理.

陶建武： 男，1959年生，博士，教授，博士生导师，从事阵列信号处理及应用、非平稳信号处理等方面的研究工作.

Single Snapshot Airspeed Estimation Based on Sparse Covariance Matrix Iteration

Yu Fei①②Tao Jian-wu②Chen Cheng③Qian Li-lin②①(Department of Control Engineering, Naval Aeronautical and Astronautical University, Yantai 264001, China)
②(Department of Aircraft Control Engineering, Aviation University of Air Force, Changchun 130022, China)
③(94638 Department, Nanchang 330201, China)

The issue of single snapshot airspeed estimation is researched based on acoustic sensor array. According to the propagation property of acoustic waves in subsonic and supersonic air current, the output model of acoustic sensor array is constructed for a given measuring equipment. Then an airspeed estimation algorithm based on Sparse Covariance Matrix Iteration with a Single Snapshot (SCMISS) is presented. SCMISS has several unique features not shared by other sparse estimation methods: it does not require the user to make any difficult selection of regularization parameters, and it has lower computational complexity and better real-time. What is more, the proposed algorithm can be applied to both subsonic and supersonic circumstances with single snapshot measurement. Finally, a compact expression for the Cramér-Rao Bound (CRB) on the estimation error of airspeed is derived to evaluate the performance of the proposed algorithm. Simulations are implemented to show the effectiveness of SCMISS.

Array signal processing; Airspeed estimation; Single snapshot; Sparse parameter estimation; Acoustic sensor array; Cramér-Rao Bound (CRB)

TN911.7; TN06

A

1009-5896(2015)03-0574-06

10.11999/JEIT140668

2014-05-21收到，2014-09-24改回

国家自然科学基金(61172126)和吉林省自然科学基金(20140101073 JC)资助课题

*通信作者：虞飞 yufei19871128@163.com