MIMO超聲波傳感器陣列風速風向測量方法

李新波， 王曉玉， 朱閣彥

(吉林大學 通信工程學院，吉林 長春 130022)

風是相對于大氣表面的空氣運動，在航空領域、航海領域、工業領域、農業領域、環境監控領域，風參數的測量都發揮著重要作用[1]。常用的測風方法有機械式、熱式、超聲波式等。機械式測風儀以風杯式、螺旋槳式為主[2]，主要基于機械旋轉部件來感知風參數，但其機械磨損損耗大，影響測量精度和使用壽命。熱式風速儀[3-4]不存在機械磨損，適用于高分辨率以及低風速測量場合，但是隨著風速的增加，靈敏度下降，而且其中的熱量元件較脆弱容易損壞[5]。超聲波測量技術與機械式、熱式測量相比，無機械磨損，測量速度快，維護成本低，但測量難點在于對超聲波飛行時間(time of flight ，TOF)的檢測精度會直接影響風參數信息的測量精度[6-7]。

為了發掘多維空域信息，提高風參數的測量精度，基于陣列信號處理思想，文獻[8]設計了一發多收的弧形超聲波傳感器陣列結構。在此基礎上，文獻[9-10]分別采用MVDR(minimum variance directionless response)算法和MUSIC(multiple signal classification)算法實現了二維風矢量的測量。但目前所提出的陣列結構大多只存在一個發射陣元，在低信噪比或環境因素較差時，風參數估計精度會有所下降。而且空間譜的構建及譜峰搜索占用了MUSIC算法的絕大多數時間，是影響波達方向估計實時性的主要因素[11]。文獻[12]通過計算陣元接收數據之間的延時相關函數得到新的陣列輸出矩陣，利用變尺度混沌優化算法進行譜峰搜索，但是以上算法的實現也非常復雜。文獻[13]提出基于傳播算子的空間平滑技術避免特征值分解，但仍需要全局譜峰搜索。本文通過二級譜峰搜索的方式提高MUSIC算法的實時性。

由于超聲波具有指向性強、遇到障礙物會被反射[14]的特點，文獻[15]提出一種基于反射的超聲波風向風速測量裝置，通過時差法原理進行風參數的測量；文獻[16]證明了縱波在0～45°的入射角度范圍內的第1個反射面可以發生全反射。而且超聲波在波束角內遇到障礙物，波束角內的超聲波都會發生全反射[17]。因此，本文提出基于反射的多發多收超聲波傳感器陣列測風方法，形成虛擬陣元，增大陣列孔徑，通過CRB理論分析和仿真實驗驗證了本文所提方法的有效性以及測量精度的優越性。

1 信號模型

1.1 構建陣列結構

陣列結構如圖1所示，其中M1～M4為發射陣元，N1～N4為接收陣元，O為反射物，其中反射物的尺寸與各接收陣元和發射陣元相同，其端面與豎直方向逆時針夾角為45°。圖中陣元以及反射物只代表其位置，并不代表其幾何形狀。由于無線信道中信號的傳輸情況是極其復雜的，為了得到一個比較有用的參數化模型，必需建立關于接收天線陣的假設，即陣元由位于空間已知坐標處的無源陣元按一定的形式排列而成，其接收特性僅與其位置有關而與其尺寸無關(認為其是一個點)，并且陣元都是全向陣元，增益均相等，相互之間的互耦忽略不計，這是陣列信號處理理論的基礎[18]。文獻[10，19]在進行風信號檢測時也忽略了陣元結構對風參數估計的影響。因此，本文可以近似認為接收陣元、發射陣元以及反射物為一個點，接收特性只與其位置有關，與尺寸無關。

圖1 多發多收陣列結構Fig.1 Multiple send and receive array structure

發射陣元發射的超聲波波束經過反射物后，波束方向發生改變，使發射波束被接收傳感器陣列所接收，實現多發多收的超聲波傳感器陣列結構。所有陣元分布在半徑為R的圓周上。發射信號Mi發生全反射的位置恰好為Ni。發射陣元位置關于水平線對稱，則通過全反射得到4個接收陣元位置關于豎直線對稱。由于發射陣元同時發射超聲波信號，在無風時刻，接收陣列同時接收信號，不存在時間延遲；在有風時刻，風信號的存在會影響發射信號到達反射物的時間，以及反射信號到達接收陣元的時間，產生時間延遲。根據超聲波信號在傳播過程中產生的時間延遲構建發射陣列流型矩陣以及接收陣列流型矩陣，通過DOA估計算法進行風參數的測量。雖然超聲波在傳播過程中存在衰減和遮擋的問題，但對于時間延遲的計算影響不大。因此不需要考慮這類問題。

為了使發射信號Mi被所有接收陣元所接收，以反射物與Ni為水平線其接收陣元Ni的±45°范圍內都應存在超聲波信號。根據選取的接收陣元個數N以及發射信號角度范圍，計算得到接收陣列中相鄰兩陣元與反射點連線的夾角α=45°/(N-1)。α隨著接收陣元個數的增加而減小，考慮到硬件成本和陣列孔徑問題，選取四陣元的接收陣列，此時α=15°，Mi發射的信號可以被所有接收陣元Ni所接收。

1.2 陣列流型推導

圖1中V為待測風速的大小，V與豎直順時針方向的夾角θ為待測風向角，定義風速與水平逆時針方向的夾角θ′，則θ+θ′=π/2。其發射陣列的陣列流型和接收陣列的陣列流型都包含風速、風向信息。由圖1中幾何關系，根據矢量分解定理，可以得到風速在反射物與發射傳感器Mi連線上的分量：

根據τMi得到發射陣列的陣列流型為：

a(θ,V)=[e-j2πfτM1e-j2πfτM2e-j2πfτM3e-j2πfτM4]T根據τNi得到接收陣列的陣列流型為：

b(θ,V)=[e-j2πfτN1e-j2πfτN2e-j2πfτN3e-j2πfτN4]T

式中f為超聲波信號頻率。

1.3 接收信號的匹配濾波處理

發射信號S=[s1,s2,s3,s4]T，其中s1、s2、s3、s4之間相互正交，即相鄰兩發射信號之間相位差均為90°。目標回波數學模型可表示為：

x=ηa(θ,V)TS

(1)

式中：η?C表示反射信號的復幅度；a(θ,V)為發射信號的導向矢量。接收陣列信號可表示為：

y=b(θ,V)ηa(θ,V)TS+E

(2)

式中：b(θ,V)為接收信號的導向矢量；y∈CN×L是接收天線的輸出;E∈CN×L為零均值，空間白和時間白的復高斯噪聲。

對于接收陣元為N的系統，若同時發射M個線性獨立正交波形，在接收端通過匹配濾波將來自不同發射天線的信號分離開，因此可以形成MN個相互獨立的通道，相當于合成了虛擬陣元，增大了陣列孔徑。利用發射信號與陣列接收數據進行匹配濾波，第i個快拍數下得到的匹配濾波輸出為：

(3)

式中：L為快拍數；ηi為反射系數;Res=Ei(n)Si(n)H;Xi∈CN×M，把Xi的列依次排成一列得到：

vec(Xi)=ηiw+v

(4)

Y=[vec(X1)vec(X2)…vec(XL)]

(5)

2 測風算法及性能分析

2.1 二級空間MUSIC算法

多輸入多輸出系統根據式(5)得到經過匹配濾波處理的陣列輸出快拍數據Y。將Y的數據協方差矩陣進行特征值分解，得到信號子空間和噪聲子空間，利用2個子空間的正交性估計風速風向參數。由于實際接收數據矩陣Y是有限長的，陣列接收數據的協方差矩陣采用最大似然估計，即：

(6)

對R進行特征值分解得到噪聲子空間V2，進而得到風參數估計的MUSIC譜估計公式為：

(7)

在理想條件下，發射陣列與接收陣列的聯合導向矢量w與噪聲子空間V2正交，P(θ,V)趨于無窮。由于在實際情況下，w與V2不能完全正交，可以通過尋找空間譜函數P(θ,V)的最大值對應的參數值得到待估計的風速、風向值。

在進行二級譜峰搜索時，先在待估計參數范圍內進行一級搜索，即采用大于其估計精度的搜索步長找到空間譜極大值出現的大致位置，粗略估計出風速值θ和風向值V；然后確定θ的鄰域[θ-σ,θ+σ]為二級風向角搜索范圍，V的鄰域[V-ε,V+ε]為二級風速搜索范圍，其中，σ和ε均為正數。最后，采用符合其估計精度的搜索步長在二級搜索范圍內精確估計出風速風向值。

2.2 算法計算量分析

若發射陣元數為M,接收陣元數為N，每計算一次空間譜函數時，各步驟需要的計算量如下：

1)a(θ,V)和b(θ,V)的計算：

乘除法次數：9(N-1)+9(M-1)=9(N+M-2)；

加減法次數：3(N-1)+3(M-1)=3(N+M-2)；

指數運算次數：(N-1)+(M-1)=(N+M-2)。

2)Kronecker積的計算：

乘除法次數：NM。

3)空間譜函數的計算：

乘除法次數：

2NM(NM-1)+NM+1=MN(2NM-1)+1

加減法次數：

(NM-1)2+(NM-2)NM+(NM-1)=NM(2NM-3)

所以計算一次空間譜函數需要進行的計算量如下：

乘除法次數：D=9(M+N-2)+2(MN)2+1；

加減法次數：A=3(M+N-2)+MN(2MN-3)；

指數運算次數：N+M-2。

風向估計范圍為0°～360°，步長為1°，風速估計范圍為0～60 m/s，步長為0.1 m/s。在進行計算量分析時，選取一級搜索時風速步長為0.2 m/s，風向步長為2°，二級搜索時σ=10°，ε=1 m/s。MUSIC算法空間譜構建計算量分析如表1所示。

表1 算法計算量分析Table 1 Algorithm calculation analysis

從表1可以看出，二級空間譜的構建極大減少了算法計算量，并且隨著陣元數的增加、一級搜索步長的增加以及二級搜索參數的減小，該方法對系統實時性的提高效果越明顯。

2.3 估計方差與克拉美羅界

在大快拍數的情況下，MUSIC算法的估計誤差是均值為零的漸進聯合高斯分布,文獻[10]推導了超聲波傳感器陣列結合MUSIC算法風參數估計方差簡化公式，并且證明了估計方差只與陣列流型，噪聲功率以及快拍數有關。因此，通過推導得出當風向一定時，風速估計方差為：

(8)

(9)

式中：(·)*表示矩陣的共軛轉置運算；⊙為Hadamard積；σ2為噪聲功率；L為快拍數；w*Tw、HV和Hθ為：

(10)

(11)

(12)

對于參數估計問題，方差是針對某種特定的估計量而言的，而克拉美羅界與具體估計方式無關，它反映的是已有信息所能估計參數的最好效果，為任何無偏估計量的方差確定一個下限，估計方差越接近克拉美羅界，估計效果越好。因此，可以用克拉美羅界作為估計性能好壞的標準。

根據文獻[20]給出的克拉美羅界定義式，得到風速的克拉美羅界公式以及風向的克拉美羅界公式分別為：

(13)

(14)

式中信號源S(i)=diag{s1(i),s2(i),s3(i),s4(i)}。

3 仿真結果分析

實驗條件：超聲波傳感器陣列結構如圖1所示。選取發射與接收陣列的弧形半徑0.1 m，發射陣列中各發射陣元發射的相互正交的超聲波頻率40 kHz，聲速340 m/s；陣元噪聲為零均值，復高斯隨機噪聲；反射系數在采樣時刻恒定不變，不同快拍數下獨立變化；風向角掃描范圍0°～360°，風速掃描范圍0～60 m/s，快拍數為5 000。

3.1 可行性實驗

為了驗證本文所提方法的可行性，設計仿真實驗。選取一級搜索風速步長為0.2 m/s，風向步長為2°，二級搜索參數σ=10°，ε=1 m/s。在信噪比為-10 dB時隨機選取風速風向參數(45°,5.4 m/s)，得到如圖2所示的多發多收二級譜峰圖和如圖3所示的一發多收全局譜峰圖。

圖2 θ=45°，V=5.4 m/s，二級譜峰圖Fig.2 θ=45°，V=5.4 m/s， second-level peak

圖3 θ=45°，V=5.4 m/s，一發多收譜峰圖Fig.3 θ=45°，V=5.4 m/s， one-transmit and multiple-receive second-level peak

由仿真得到的譜峰圖可知，對于隨機選取的風參數，在低信噪比情況下，多發多收陣列可以無差估計，而一發多收的陣列無法進行風參數估計。

3.2 估計方差與克拉美羅界實驗

在信噪比為10 dB，快拍數為5 000的條件下，通過式(8)～(14)計算出測量范圍內每個風速和風向值所對應的估計方差與克拉美羅界，并計算二者的差值，如圖4和圖5所示。

圖4 風向角估計方差與CRB之差Fig.4 Difference between wind direction angle estimation variance and CRB

圖5 風速估計方差與CRB之差Fig.5 Difference between wind speed estimation variance and CRB

由仿真結果可知，所有的風向角估計方差與CRB之差均小于0.04；風速估計方差與CRB之差均小于1×10-3。風速和風向的估計方差與克拉美羅界較為接近，由此可知，本文所提方法對測量范圍內所有風速、風向值的估計具有較小的偏差。

3.3 不同信噪比下風參數估計均方根誤差實驗

采用蒙特卡洛仿真實驗來評估本文所提方法的風速、風向估計性能。均方根誤差公式為：

圖6 風速、風向均方根誤差Fig.6 Root mean square error of wind speed and direction

由仿真結果可知，隨著信噪比的增大，風速和風向角的均方根誤差都逐漸減小。本文所提方法當信噪比分別增大到-25 dB和-20 dB時，風向角和風速均方根誤差分別減小為零；而一發四收的弧形陣列當信噪比增大到11 dB時，風向角和風速均方根誤差才減小為零。因此，在低信噪比情況下本文所提方法有效提高了風參數的估計精度，獲得了良好的估計性能。雖然風速、風向參數是隨機選取的，但從3.2節估計方差仿真實驗中可以看出，在同一信噪比條件下，對整個測量范圍內參數的估計誤差相差不大。因此，隨機選取的風速、風向參數其估計誤差可以代表整個測量范圍內參數估計誤差的趨勢。

3.4 不同信噪比下風參數估計成功率實驗

實驗條件與3.3節相同，風速估計誤差小于其步長0.1 m/s時，認為風速估計成功；風向估計誤差小于其步長1°時，認為風向估計成功。估計成功率隨信噪比變化曲線如圖7所示，圖中虛線為采用二級空間MUSIC算法對風參數的估計結果；實線為直接采用符合參數精度的MUSIC算法得到的估計結果。

圖7 風速、風向估計成功率對比Fig.7 Wind speed and wind direction estimation success rate

由仿真結果可以看出，隨著信噪比的增加，風參數估計成功率增大。雖然在低信噪比下，二級空間MUSIC算法對參數估計成功率略低于MUSIC算法，當信噪比分別增大到-25 dB和-20 dB時，風向和風速估計成功率分別達到100%。說明二級空間譜的構建方式不僅可以減小算法計算量，而且具有較高的估計精度。通過以上實驗和理論分析，驗證了本文所提方法具有較好的風參數測量效果。

4 結論

1)本文針對風速風向測量精度低的問題，基于陣列信號處理理論，提出多輸入多輸出的超聲波傳感器陣列結構的風速風向測量方法，提高了低信噪比下的風參數的測量精度。

2)本文基于MUSIC算法，提出二級空間譜搜索方法，來估計風速風向值，避免了全局譜峰搜索，減小了算法計算量，提高了系統實時性。

3)為了驗證本文所提陣列結構和算法的風速風向測量性能，進行了風參數估計方差和克拉美羅界的理論分析和計算公式的推導，并通過仿真對比試驗，驗證了所提方法在測量精度以及系統實時性方面的優勢。