基于稀疏協方差矩陣迭代的單快拍氣流速度估計算法

虞 飛陶建武陳 誠錢立林

①(海軍航空工程學院控制工程系 煙臺 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