三維超聲陣列風速風向測量方法

李新波,朱閣彥,李厚禹,賈云龍

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

風是一種由許多在時空上隨機變化的小尺度脈動疊加而成的自然現(xiàn)象,也是在大尺度規(guī)則氣流上的一種三維矢量,主要包括風速、風向角和俯仰角這3個參數(shù)。作為一種很常見的自然現(xiàn)象,風的準確測量在工業(yè)、氣象和航運等領域中,正發(fā)揮著越來越重要的作用[1-3]。近年來,基于超聲檢測的風參數(shù)測量技術越來越成熟,且與基于機械式、熱敏式和激光多普勒式等的測風技術相比,具有無磨損、測量范圍大和維護成本低等優(yōu)勢,因此被廣泛關注[4-8]。在測風技術中,最常用的測風方法是時差法[9-11],即通過測量在順風和逆風情況下,超聲波信號在傳感器發(fā)射端和接收端之間傳播的不同時間來進行風參數(shù)的測量,因此基于時差法的測風儀器的測量精度完全取決于超聲波傳播時間的測量精度,但由于環(huán)境噪聲的影響,想要通過儀器精準測量超聲波傳播時間并不容易。相對于平面風參數(shù)測量,國內外對空間三維風參數(shù)測量的研究較少。鄧云逸等利用正四面體的超聲波傳感器陣列結構,并依據(jù)時差法的測量原理,分別測出水平面和垂直面的風速風向,然后通過矢量合成得到了空間三維風參數(shù)[12];Frank等通過聯(lián)合安裝在水平面和垂直面上的測風儀的測量結果,利用空間三維矢量合成計算出了三維風參數(shù)[13];Lopes等討論分析了超聲測風儀結構對測風精度的影響,并提出了一種幾何結構非正交的超聲測風儀,顯著提高了基于時差法進行風參數(shù)測量的準確性[14]。但是,這些研究還是在基于時差法進行平面風參數(shù)測量的基礎上,通過空間矢量的分解與合成計算出三維風參數(shù),并沒有克服基于時差法進行風參數(shù)測量的弊端。Li等在2015年提出利用陣列信號處理算法中的多重信號分類(multiple signal classification,MUSIC)的思想來計算風參數(shù),該方法通過提取特征子空間和掃描功率譜來估計平面風參數(shù),結果表明,所提方法相對于傳統(tǒng)的時差法在抑制噪聲和提高測量精度方面都有很大的優(yōu)越性[15],但MUSIC算法需要對陣列接收信號進行奇異值分解[16],較大的計算量在一定程度上限制了實際應用中對系統(tǒng)實時性的要求。

本文將陣列信號處理中典型的波束形成算法應用到三維風參數(shù)測量領域,在計算量相對較小的情況下,實現(xiàn)了三維風參數(shù)的測量。首先,設計了一種空間超聲波傳感器陣列結構,并基于該結構分析了利用波束形成算法進行三維風參數(shù)測量的原理。然后,在Matlab軟件中進行了所提測風方法的可行性實驗及測量性能實驗。最后,搭建了三維超聲測風實驗平臺,對三維空間風進行了實物測量實驗。

1 超聲陣列結構設計及測風原理

1.1 超聲陣列結構設計

三維超聲測風裝置結構如圖1所示,主要包括超聲波發(fā)射單元、超聲波接收單元以及支撐單元。超聲波發(fā)射單元包括1個超聲波發(fā)射傳感器,也稱發(fā)射陣元,能在激勵信號的作用下產(chǎn)生振動,發(fā)射出特定頻率的超聲波信號。超聲波接收單元由8個超聲波接收傳感器(a1～a4,b1～b4)組成,每4個組成一組,分別安裝在相互垂直的弧形支撐架上,用來接收發(fā)射陣元所發(fā)射的超聲波信號,超聲波接收傳感器也稱接收陣元,發(fā)射陣元和接收陣元是型號相同、且都具有收發(fā)一體化功能的超聲波傳感器。支撐單元是支撐架,用來連接、支撐和固定各個超聲波傳感器之間的相對位置。

圖1 三維超聲測風裝置結構圖

由圖1可知:處于垂直面yoz上的4個接收陣元a1、a2、a3、a4組成子陣列Aa,處于水平面xoy上的4個接收陣元b1、b2、b3、b4組成子陣列Ab,子陣列Aa和子陣列Ab共同組成該測風裝置的傳感器陣列A;在同一子陣列中,任意兩個相鄰的接收陣元與發(fā)射陣元之間連線的夾角為20°;發(fā)射陣元到各接收陣元之間的距離相等,即8個接收陣元在同一個以發(fā)射陣元o為圓心、r為半徑的球面上。根據(jù)該測風裝置的各陣元之間的相對結構可知:當無風時,從發(fā)射陣元發(fā)射的超聲波信號同時傳播到8個接收陣元;當有風時,超聲波信號受風的影響,傳播速度發(fā)生變化,從而到達各接收陣元的時間不一致。

1.2 基于波束形成算法的測風原理

超聲波在空氣中傳播時,風速會在其聲速上疊加,使超聲波傳播速度變快或變慢。因此,利用超聲波傳播速度和超聲波聲速之間的關系,便能測量出風參數(shù)[17]。

由于超聲波聲速c會受到傳播介質的溫度、濕度等因素的影響,因此有必要其進行修正來提高測風裝置的測量精度[18]。目前,關于超聲波聲速與溫度之間關系的研究較為成熟,工程上常用的表述溫度T與超聲波聲速c之間的經(jīng)驗公式[10]為

(1)

本文用式(1)來修正超聲波聲速,下文中所使用的c即為修正之后的超聲波聲速。

空間風的來向可用圖2表示,圖中:V為風速;θ為風向角;φ為俯仰角,0°≤φ≤180°;α為在同一子陣列中的任意兩個相鄰的接收陣元與發(fā)射陣元之間連線的夾角,α=20°;Va為空間來風在垂直面yoz上的風速分量;Vb為空間來風在水平面xoy上的風速分量。根據(jù)空間幾何分解的相關知識可知,Vb在水平面xoy上的分量如圖3所示,Vb1～Vb4分別表示4個接收陣元方向的風速分量。

圖2 空間風來向圖

圖3 Vb在水平面xoy上的分量示意圖

空間風在xoy平面上各接收陣元方向上的風速分量可表示為

(2)

由于發(fā)射陣元所發(fā)射的超聲波信號s(t)為單頻窄帶信號,因此可表示為

s(t)=u(t)ej[ω t+φ(t)]

(3)

式中:t為時間;u(t)為信號的幅度;φ(t)為信號的相位;角頻率ω=2πf,f=40 kHz為信號的頻率。

經(jīng)過時間延遲τ之后,超聲波信號可表示為

s(t-τ)=s(t)e-jωτ

(4)

子陣列Ab所接收的陣列信號xb(t)可表示為

(5)

式中:xbi為子陣列Ab中陣元bi所接收的超聲波信號;nbi(t)為子陣列Ab中陣元bi上的空域采樣噪聲,且各陣元接收的噪聲彼此獨立;τi為超聲波信號傳播到陣元bi相對于傳播到參考陣元的時間延遲,簡稱時延。

將式(5)以矢量形式表示,可以得到

xb(t)=s(t)Db+nb(t)

(6)

(7)

(8)

式中:Db為子陣列Ab的陣列流型矢量;nb(t)為子陣列Ab所接收的采樣噪聲的矢量形式。

超聲波信號從發(fā)射陣元傳播到子陣列Ab中的接收陣元bi(i=1,2,3,4)所需的時間ti可表示為

(9)

若將子陣列Ab中的接收陣元b4作為參考陣元,則時延可表示為

(10)

(11)

(12)

τ4=t4-t4=0

(13)

由以上分析可知,子陣列Ab的陣列流型矢量Db僅與時延τi有關,而τi是由空間風在平面xoy上的風矢量分量信息風速Vb和風向θ確定的,所以在超聲波傳感器陣列結構確定的前提下,子陣列Ab的陣列流型矢量Db直接由風速Vb和風向θ確定,即Db=Db(Vb,θ)。反過來說,確定了子陣列Ab的陣列流型矢量Db,也就確定了空間風在平面xoy上的風矢量分量信息風速Vb和風向θ。

根據(jù)波束形成算法中加權系數(shù)的常規(guī)確定方法可知,在空間陣列方向辨識應用中,若空間中只有一個來自方向為θk的信號s(t),其方向矢量為a(θk),則當加權矢量W取作a(θk)時,陣列輸出為X(t)=WHs(t)=aH(θk)s(t),對其進行功率求解后得到的功率最大[19-20]。類比這種思想,可認為Db是風速風向矢量,當取加權矢量W=Db時,陣列輸出為y(t)=WHxb(t),此時進行功率求解后功率最大。

根據(jù)常規(guī)波束形成算法,結合待估計參數(shù)為xoy平面上具有2參數(shù)的風速風向可知,子陣列Ab中各接收陣元所接收信號的加權輸出為

(14)

將式(14)以矢量表示,可以得到

(15)

此時,子陣列Ab的陣列輸出功率為

P(θ,Vb)=E[y2(t)]=WHRW=

(16)

式中矩陣R為子陣列Ab的陣列接收信號xb(t)的協(xié)方差矩陣,即

(17)

根據(jù)式(16),在Vb、θ分別處于0～60 m/s和0～359°的范圍內,以確定步長遍歷風速和風向,即可得到子陣列Ab的陣列接收信號的輸出功率,其中最大功率所對應的風速和風向即為空間風在水平面xoy上的風矢量分量Vb和θ。

同理,對子陣列Aa接收的陣列信號做同樣的處理,便能確定出空間風在垂直面yoz上的風矢量分量信息風速Va和風向φyoz,本節(jié)不再贅述。

圖4 空間三維風矢量的分解與合成示意圖

圖4是空間三維風矢量的分解與合成示意圖,可以看出,在風參數(shù)矢量信息(Vb,θ)和(Va,φyoz)已知的情況下,球坐標系下的空間矢量(V,θ,φ)的終點在直角坐標系下x、y、z軸上的投影坐標可分別表示為

Vx=Vbsinθ

(18)

Vy=Vbcosθ=Vacosφyoz

(19)

Vz=Vasinφyoz

(20)

所以,球坐標系下的空間矢量(V,θ,φ)也可以在空間直角坐標系中表示為(Vbsinθ,Vacosφyoz,Vasinφyoz),該矢量的模為

(21)

空間風矢量的來向與z軸正方向之間所成夾角φ的余弦值可表示為

(22)

夾角φ可表示為

(23)

通過上述分析可知,在已知空間三維風矢量在水平面xoy和垂直面yoz上的風矢量分量信息(Vb,θ)和(Va,φyoz)時,便能確定出空間中的三維風參數(shù),即風速、風向角和俯仰角。

2 仿真實驗結果及分析

2.1 可行性實驗

為了驗證本文測風方法的可行性,在Matlab軟件上進行了可行性仿真驗證。首先,隨機給出確定的空間風矢量在水平面xoy和垂直面yoz上的風矢量分量信息(Vb,θ)和(Va,φyoz),并計算出三維風參數(shù)(V,θ,φ)作為真實空間中三維風參數(shù)。然后,將子陣列Ab和Aa的陣列接收信號分別進行波束形成算法處理,得到2個相互垂直平面上的風矢量分量信息,并利用這2個風矢量分量信息計算出空間三維風參數(shù),以此作為仿真結果。最后,將仿真結果與真實結果(V,θ,φ)進行比較,以此驗證本文測風方法的可行性。

實驗條件如下:模擬發(fā)射的超聲波信號頻率為40 kHz;發(fā)射陣元到各接收陣元之間的距離為r=10 cm;各接收陣元上的噪聲為加性高斯白噪聲;速度的掃描范圍為0～60 m/s,步長為0.1 m/s;角度的掃描范圍為0～359°,步長為1°。當快拍數(shù)為5 000、信噪比rsn=5 dB時,仿真結果如表1所示。

表1 三維風參數(shù)的真實值和仿真值

由仿真結果可知,當信噪比rsn=5 dB時,空間風的風速、風向角和俯仰角都能無差地估計出來。因此,可認為本文提出的方法是可行的。

2.2 成功率實驗

由于存在環(huán)境噪聲,估計結果相對于真實結果會有一定的偏差,為了驗證本文方法的估計性能,設計了在不同信噪比下三維風矢量估計結果成功率的實驗。由于速度掃描步長為0.1 m/s、角度掃描步長為1°,所以在該實驗中,當速度估計的誤差小于0.1 m/s、角度估計誤差小于1°時,都認為估計成功。

實驗條件如下:信噪比rsn的取值范圍為-6～20 dB,步長為1 dB。在每個信噪比下做400次蒙特卡羅實驗,分別估計可行性實驗中的3組隨機參數(shù),結果如圖5～7所示。

圖5 風速估計成功率

圖6 風向角估計成功率

圖7 俯仰角估計成功率

由仿真結果可知,在可接受誤差存在的前提下,三維風參數(shù)的估計成功率隨著信噪比的增加而增加,且在信噪比rsn≥5 dB時,成功率幾乎達到100%。

2.3 均方根誤差實驗

為了討論估計結果與真實結果的偏差度,設計了在不同信噪比的條件下,三維風矢量信息估計結果均方根誤差RMSE的仿真實驗,實驗條件與2.2小節(jié)相同。均方根誤差的公式

(24)

圖8 風速均方根誤差

圖9 風向角均方根誤差

圖10 俯仰角均方根誤差

由仿真結果可知:在低信噪比條件下,速度的估計效果好于角度的估計效果;隨著信噪比的增加,三維風參數(shù)的估計結果越來越趨向于真值,在信噪比大于5 dB時,速度的均方根誤差小于0.1 m/s,角度的均方根誤差值小于1°。

通過以上3個仿真實驗及仿真結果可知,本文所提的測風方法在理論上具有較好的三維風參數(shù)測量效果。

3 實驗與結果

3.1 實驗平臺

為了驗證本文提出的方法是否具有工程可實現(xiàn)性,搭建了基于超聲波傳感器陣列的三維風參數(shù)測量實驗平臺,如圖11所示。直流可調穩(wěn)壓電源(HY3005ET)為風機提供穩(wěn)定的電壓,可調穩(wěn)壓電源(HSPY-200-02)作為控制單元(ATMEGA328P)的電源,使控制單元發(fā)送頻率為40 kHz的脈寬調制,并通過驅動電路驅動發(fā)射陣元發(fā)射頻率為40 kHz的超聲波信號。風洞的整體框架由透明的亞克力板制作而成,如圖12所示。當測風裝置置于風洞中時,發(fā)射陣元發(fā)射的超聲波信號在空氣中傳播會受風的影響,從而傳播速度發(fā)生變化,所以超聲波信號在發(fā)射陣元和接收陣元之間的傳播時間就會發(fā)生變化,具體表現(xiàn)為接收信號之間相位差的變化。本實驗采用4通道同步采樣的示波器(MSO7054B),同步采樣頻率為20 MHz,每路信號最多可同時采集1 000個離散的數(shù)據(jù)點,保存成.csv格式的文件。

圖11 三維風參數(shù)測量實驗平臺

圖12 實驗風洞

3.2 風洞中的風參數(shù)測量實驗

在采集信號之前,需要先將測風裝置校準,使同一平面上4路通道所接收的信號在無風時保持同相。但是,由于沒有精確的安裝技術,在人為固定各超聲波傳感器的過程中,會存在較大的安裝誤差。經(jīng)過一系列的結構調整和校準,無風時水平面4路接收信號在示波器上顯示的波形如圖13所示,可以看出,最終只能達到接收陣元b2、b3所接收的信號保持同相的水平。同時,在垂直面上只能保證接收陣元a2、a3所接收的信號保持同相。此時,接收陣元b2、b3、a2、a3接收的信號如圖14所示。因此,在后續(xù)的實驗過程中,只利用接收陣元b2、b3、a2、a3所接收的數(shù)據(jù)進行運算。

(a)陣元b1

(b)陣元b2

(c)陣元b3

(d)陣元b4圖13 無風時水平面4路陣元接收信號的波形

(a)陣元b1

(b)陣元b2

(c)陣元a2

(d)陣元a3圖14 無風時接收陣元b2、b3、a2、a3接收信號的波形

(a)陣元b1

(b)陣元b2

(c)陣元a2

(d)陣元a3圖15 有風時接收陣元b2、b3、a2、a3接收信號的波形

有風時接收陣元b2、b3、a2、a3接收的信號如圖15所示,可以看出,陣元b2、b3接收的信號和陣元a2、a3接收的信號之間出現(xiàn)明顯的相位差。利用Matlab軟件調用有風時的采樣數(shù)據(jù).csv文件,便能計算出風洞中的三維風參數(shù)。測量環(huán)境的溫度為14.9 ℃,根據(jù)式(1)可得此時的超聲波聲速c=340.323 6 m/s,利用該修正后的超聲波聲速及采樣數(shù)據(jù),計算可得測量的風洞中的三維風參數(shù),結果如表2所示。考慮到超聲測風裝置的安裝誤差以及實物實驗測量中使用的接收陣元數(shù)量減半等因素的影響,超聲測風裝置測量的三維風參數(shù)相對于真值會不可避免地存在誤差。由于風洞中風場的不穩(wěn)定性和時變性,很難獲得某一時刻風參數(shù)的真值,因此可通過借助市場上的測風儀來測量風洞中風速的最大值,進行本文方法可行性的驗證。由于采用機械式測風儀(SMART SENSOR AS8556)測量的風洞中最大風速為6.443 m/s,所以當本文實驗平臺測量的風速小于6.443 m/s且在其附近時,便認為是合理的。本文實驗平臺所測風速為6.363 7 m/s,與機械式測風儀測量結果僅相差0.079 3 m/s,因此本文所提測風方法是具有工程實現(xiàn)性的。風向角和俯仰角會因測風裝置所處的位置和風洞中氣流的變化而變化,所以本節(jié)不討論。本節(jié)僅說明本文方法具有工程實現(xiàn)性,更精確的測量需要更專業(yè)的儀器在更高規(guī)格的實驗室或機構中進行。

表2 本文實驗平臺測量結果

通過以上分析可知,在允許一定結構誤差和設備誤差的前提下,本文測風方法具有工程實現(xiàn)性。同時,如何克服結構誤差和減小設備誤差也是后續(xù)工作中需要努力改進的地方。

4 結 論

(1)本文設計了一種新的超聲波傳感器陣列結構,并將波束形成算法應用于三維風參數(shù)測量中。

(2)充分利用了超聲波傳感器陣列接收信號的冗余信息和波束形成算法的空域濾波、干擾抑制、期望信號增強等優(yōu)點,在提高風參數(shù)的測量精度方面具有一定的理論優(yōu)勢。

(3)通過Matlab軟件對本文提出的三維風參數(shù)測量方法進行了仿真實驗,結果表明,在信噪比大于5 dB時,三維風參數(shù)測量信息幾乎能達到100%的成功率,風速的均方根誤差小于0.2 m/s,風向角的均方根誤差小于1°。

(4)搭建了三維超聲測風實驗平臺并對風洞中的風參數(shù)進行了實際測量,結果表明,本文系統(tǒng)與市場上成型的機械式測風儀的風速測量結果僅相差0.078 3 m/s,說明本文方法具有工程可行性。