李新波， 王曉玉， 朱閣彥

(吉林大學(xué) 通信工程學(xué)院，吉林 長(zhǎng)春 130022)

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

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

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

1 信號(hào)模型

1.1 構(gòu)建陣列結(jié)構(gòu)

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

圖1 多發(fā)多收陣列結(jié)構(gòu)Fig.1 Multiple send and receive array structure

發(fā)射陣元發(fā)射的超聲波波束經(jīng)過(guò)反射物后，波束方向發(fā)生改變，使發(fā)射波束被接收傳感器陣列所接收，實(shí)現(xiàn)多發(fā)多收的超聲波傳感器陣列結(jié)構(gòu)。所有陣元分布在半徑為R的圓周上。發(fā)射信號(hào)Mi發(fā)生全反射的位置恰好為Ni。發(fā)射陣元位置關(guān)于水平線對(duì)稱(chēng)，則通過(guò)全反射得到4個(gè)接收陣元位置關(guān)于豎直線對(duì)稱(chēng)。由于發(fā)射陣元同時(shí)發(fā)射超聲波信號(hào)，在無(wú)風(fēng)時(shí)刻，接收陣列同時(shí)接收信號(hào)，不存在時(shí)間延遲；在有風(fēng)時(shí)刻，風(fēng)信號(hào)的存在會(huì)影響發(fā)射信號(hào)到達(dá)反射物的時(shí)間，以及反射信號(hào)到達(dá)接收陣元的時(shí)間，產(chǎn)生時(shí)間延遲。根據(jù)超聲波信號(hào)在傳播過(guò)程中產(chǎn)生的時(shí)間延遲構(gòu)建發(fā)射陣列流型矩陣以及接收陣列流型矩陣，通過(guò)DOA估計(jì)算法進(jìn)行風(fēng)參數(shù)的測(cè)量。雖然超聲波在傳播過(guò)程中存在衰減和遮擋的問(wèn)題，但對(duì)于時(shí)間延遲的計(jì)算影響不大。因此不需要考慮這類(lèi)問(wèn)題。

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

1.2 陣列流型推導(dǎo)

圖1中V為待測(cè)風(fēng)速的大小，V與豎直順時(shí)針?lè)较虻膴A角θ為待測(cè)風(fēng)向角，定義風(fēng)速與水平逆時(shí)針?lè)较虻膴A角θ′，則θ+θ′=π/2。其發(fā)射陣列的陣列流型和接收陣列的陣列流型都包含風(fēng)速、風(fēng)向信息。由圖1中幾何關(guān)系，根據(jù)矢量分解定理，可以得到風(fēng)速在反射物與發(fā)射傳感器Mi連線上的分量：

根據(jù)τMi得到發(fā)射陣列的陣列流型為：

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

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

式中f為超聲波信號(hào)頻率。

1.3 接收信號(hào)的匹配濾波處理

發(fā)射信號(hào)S=[s1,s2,s3,s4]T，其中s1、s2、s3、s4之間相互正交，即相鄰兩發(fā)射信號(hào)之間相位差均為90°。目標(biāo)回波數(shù)學(xué)模型可表示為：

x=ηa(θ,V)TS

(1)

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

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

(2)

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

對(duì)于接收陣元為N的系統(tǒng)，若同時(shí)發(fā)射M個(gè)線性獨(dú)立正交波形，在接收端通過(guò)匹配濾波將來(lái)自不同發(fā)射天線的信號(hào)分離開(kāi)，因此可以形成MN個(gè)相互獨(dú)立的通道，相當(dāng)于合成了虛擬陣元，增大了陣列孔徑。利用發(fā)射信號(hào)與陣列接收數(shù)據(jù)進(jìn)行匹配濾波，第i個(gè)快拍數(shù)下得到的匹配濾波輸出為：

(3)

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

vec(Xi)=ηiw+v

(4)

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

(5)

2 測(cè)風(fēng)算法及性能分析

2.1 二級(jí)空間MUSIC算法

多輸入多輸出系統(tǒng)根據(jù)式(5)得到經(jīng)過(guò)匹配濾波處理的陣列輸出快拍數(shù)據(jù)Y。將Y的數(shù)據(jù)協(xié)方差矩陣進(jìn)行特征值分解，得到信號(hào)子空間和噪聲子空間，利用2個(gè)子空間的正交性估計(jì)風(fēng)速風(fēng)向參數(shù)。由于實(shí)際接收數(shù)據(jù)矩陣Y是有限長(zhǎng)的，陣列接收數(shù)據(jù)的協(xié)方差矩陣采用最大似然估計(jì)，即：

(6)

對(duì)R進(jìn)行特征值分解得到噪聲子空間V2，進(jìn)而得到風(fēng)參數(shù)估計(jì)的MUSIC譜估計(jì)公式為：

(7)

在理想條件下，發(fā)射陣列與接收陣列的聯(lián)合導(dǎo)向矢量w與噪聲子空間V2正交，P(θ,V)趨于無(wú)窮。由于在實(shí)際情況下，w與V2不能完全正交，可以通過(guò)尋找空間譜函數(shù)P(θ,V)的最大值對(duì)應(yīng)的參數(shù)值得到待估計(jì)的風(fēng)速、風(fēng)向值。

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

2.2 算法計(jì)算量分析

若發(fā)射陣元數(shù)為M,接收陣元數(shù)為N，每計(jì)算一次空間譜函數(shù)時(shí)，各步驟需要的計(jì)算量如下：

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

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

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

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

2)Kronecker積的計(jì)算：

乘除法次數(shù)：NM。

3)空間譜函數(shù)的計(jì)算：

乘除法次數(shù)：

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

加減法次數(shù)：

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

所以計(jì)算一次空間譜函數(shù)需要進(jìn)行的計(jì)算量如下：

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

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

指數(shù)運(yùn)算次數(shù)：N+M-2。

風(fēng)向估計(jì)范圍為0°～360°，步長(zhǎng)為1°，風(fēng)速估計(jì)范圍為0～60 m/s，步長(zhǎng)為0.1 m/s。在進(jìn)行計(jì)算量分析時(shí)，選取一級(jí)搜索時(shí)風(fēng)速步長(zhǎng)為0.2 m/s，風(fēng)向步長(zhǎng)為2°，二級(jí)搜索時(shí)σ=10°，ε=1 m/s。MUSIC算法空間譜構(gòu)建計(jì)算量分析如表1所示。

表1 算法計(jì)算量分析Table 1 Algorithm calculation analysis

從表1可以看出，二級(jí)空間譜的構(gòu)建極大減少了算法計(jì)算量，并且隨著陣元數(shù)的增加、一級(jí)搜索步長(zhǎng)的增加以及二級(jí)搜索參數(shù)的減小，該方法對(duì)系統(tǒng)實(shí)時(shí)性的提高效果越明顯。

2.3 估計(jì)方差與克拉美羅界

在大快拍數(shù)的情況下，MUSIC算法的估計(jì)誤差是均值為零的漸進(jìn)聯(lián)合高斯分布,文獻(xiàn)[10]推導(dǎo)了超聲波傳感器陣列結(jié)合MUSIC算法風(fēng)參數(shù)估計(jì)方差簡(jiǎn)化公式，并且證明了估計(jì)方差只與陣列流型，噪聲功率以及快拍數(shù)有關(guān)。因此，通過(guò)推導(dǎo)得出當(dāng)風(fēng)向一定時(shí)，風(fēng)速估計(jì)方差為：

(8)

(9)

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

(10)

(11)

(12)

對(duì)于參數(shù)估計(jì)問(wèn)題，方差是針對(duì)某種特定的估計(jì)量而言的，而克拉美羅界與具體估計(jì)方式無(wú)關(guān)，它反映的是已有信息所能估計(jì)參數(shù)的最好效果，為任何無(wú)偏估計(jì)量的方差確定一個(gè)下限，估計(jì)方差越接近克拉美羅界，估計(jì)效果越好。因此，可以用克拉美羅界作為估計(jì)性能好壞的標(biāo)準(zhǔn)。

根據(jù)文獻(xiàn)[20]給出的克拉美羅界定義式，得到風(fēng)速的克拉美羅界公式以及風(fēng)向的克拉美羅界公式分別為：

(13)

(14)

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

3 仿真結(jié)果分析

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

3.1 可行性實(shí)驗(yàn)

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

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

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

由仿真得到的譜峰圖可知，對(duì)于隨機(jī)選取的風(fēng)參數(shù)，在低信噪比情況下，多發(fā)多收陣列可以無(wú)差估計(jì)，而一發(fā)多收的陣列無(wú)法進(jìn)行風(fēng)參數(shù)估計(jì)。

3.2 估計(jì)方差與克拉美羅界實(shí)驗(yàn)

在信噪比為10 dB，快拍數(shù)為5 000的條件下，通過(guò)式(8)～(14)計(jì)算出測(cè)量范圍內(nèi)每個(gè)風(fēng)速和風(fēng)向值所對(duì)應(yīng)的估計(jì)方差與克拉美羅界，并計(jì)算二者的差值，如圖4和圖5所示。

圖4 風(fēng)向角估計(jì)方差與CRB之差Fig.4 Difference between wind direction angle estimation variance and CRB

圖5 風(fēng)速估計(jì)方差與CRB之差Fig.5 Difference between wind speed estimation variance and CRB

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

3.3 不同信噪比下風(fēng)參數(shù)估計(jì)均方根誤差實(shí)驗(yàn)

采用蒙特卡洛仿真實(shí)驗(yàn)來(lái)評(píng)估本文所提方法的風(fēng)速、風(fēng)向估計(jì)性能。均方根誤差公式為：

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

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

3.4 不同信噪比下風(fēng)參數(shù)估計(jì)成功率實(shí)驗(yàn)

實(shí)驗(yàn)條件與3.3節(jié)相同，風(fēng)速估計(jì)誤差小于其步長(zhǎng)0.1 m/s時(shí)，認(rèn)為風(fēng)速估計(jì)成功；風(fēng)向估計(jì)誤差小于其步長(zhǎng)1°時(shí)，認(rèn)為風(fēng)向估計(jì)成功。估計(jì)成功率隨信噪比變化曲線如圖7所示，圖中虛線為采用二級(jí)空間MUSIC算法對(duì)風(fēng)參數(shù)的估計(jì)結(jié)果；實(shí)線為直接采用符合參數(shù)精度的MUSIC算法得到的估計(jì)結(jié)果。

圖7 風(fēng)速、風(fēng)向估計(jì)成功率對(duì)比Fig.7 Wind speed and wind direction estimation success rate

由仿真結(jié)果可以看出，隨著信噪比的增加，風(fēng)參數(shù)估計(jì)成功率增大。雖然在低信噪比下，二級(jí)空間MUSIC算法對(duì)參數(shù)估計(jì)成功率略低于MUSIC算法，當(dāng)信噪比分別增大到-25 dB和-20 dB時(shí)，風(fēng)向和風(fēng)速估計(jì)成功率分別達(dá)到100%。說(shuō)明二級(jí)空間譜的構(gòu)建方式不僅可以減小算法計(jì)算量，而且具有較高的估計(jì)精度。通過(guò)以上實(shí)驗(yàn)和理論分析，驗(yàn)證了本文所提方法具有較好的風(fēng)參數(shù)測(cè)量效果。

4 結(jié)論

1)本文針對(duì)風(fēng)速風(fēng)向測(cè)量精度低的問(wèn)題，基于陣列信號(hào)處理理論，提出多輸入多輸出的超聲波傳感器陣列結(jié)構(gòu)的風(fēng)速風(fēng)向測(cè)量方法，提高了低信噪比下的風(fēng)參數(shù)的測(cè)量精度。

2)本文基于MUSIC算法，提出二級(jí)空間譜搜索方法，來(lái)估計(jì)風(fēng)速風(fēng)向值，避免了全局譜峰搜索，減小了算法計(jì)算量，提高了系統(tǒng)實(shí)時(shí)性。

3)為了驗(yàn)證本文所提陣列結(jié)構(gòu)和算法的風(fēng)速風(fēng)向測(cè)量性能，進(jìn)行了風(fēng)參數(shù)估計(jì)方差和克拉美羅界的理論分析和計(jì)算公式的推導(dǎo)，并通過(guò)仿真對(duì)比試驗(yàn)，驗(yàn)證了所提方法在測(cè)量精度以及系統(tǒng)實(shí)時(shí)性方面的優(yōu)勢(shì)。