石屹然，齊金偉，曲思凝，趙 洋

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

1 引 言

風(fēng)速、風(fēng)向作為氣象要素的重要組成部分，在軍事、航空、航海、交通等國(guó)計(jì)民生的各個(gè)方面均有著重大的影響［1］。相比于傳統(tǒng)的如熱敏測(cè)量［2］、多普勒測(cè)量［3］、皮托管測(cè)量［4］以及機(jī)械測(cè)量［5］等測(cè)量方法，超聲波風(fēng)矢量測(cè)量方法具有的測(cè)量精度高、范圍寬、適用于各種野外復(fù)雜情況等優(yōu)點(diǎn)，使之成為國(guó)內(nèi)外學(xué)者的熱點(diǎn)研究問(wèn)題。

針對(duì)超聲波風(fēng)矢量測(cè)量問(wèn)題，Ruttgen［6］提出了一種基于相對(duì)時(shí)差法的超聲波風(fēng)速風(fēng)向測(cè)量方法，該方法利用超聲波在順、逆風(fēng)情況下渡越時(shí)間的不同來(lái)測(cè)量風(fēng)速，從而消除了環(huán)境溫濕度等參數(shù)對(duì)測(cè)量精度的影響，有效地提高了測(cè)量精度。但是該方法忽略了環(huán)境噪聲對(duì)超聲波渡越時(shí)間測(cè)量的影響，在較低信噪比條件下，該方法測(cè)量性能急劇下降。對(duì)此，Castagnede 等［7］將超聲波渡越時(shí)間測(cè)量問(wèn)題凝練為含噪信號(hào)時(shí)延估計(jì)問(wèn)題，提出了一種基于自相關(guān)時(shí)延參數(shù)估計(jì)的風(fēng)矢量測(cè)量方法，有效抑制了高斯噪聲對(duì)時(shí)延參數(shù)估計(jì)精度的影響。然而，該方法在一定程度上銳化了信號(hào)的峰值，存在峰值誤判的問(wèn)題［8］。Brassier 等［9］提出了一種基于互相關(guān)時(shí)延參數(shù)估計(jì)的超聲波風(fēng)矢量測(cè)量方法，該方法將發(fā)射信號(hào)與接收信號(hào)進(jìn)行互相關(guān)運(yùn)算，直接估計(jì)渡越時(shí)間，減小了計(jì)算結(jié)果的方差。曹長(zhǎng)宏等［10］在互相關(guān)時(shí)延估計(jì)算法的基礎(chǔ)上提出了一種二次互相關(guān)時(shí)延估計(jì)算法。該方法通過(guò)對(duì)互相關(guān)和自相關(guān)運(yùn)算結(jié)果再次進(jìn)行互相關(guān)運(yùn)算以估計(jì)超聲波渡越時(shí)間，進(jìn)一步降低了時(shí)延估計(jì)算法的信噪比閾值。但是，基于時(shí)間延遲估計(jì)的超聲波風(fēng)矢量測(cè)量方法需在時(shí)間延遲為微秒級(jí)精度下進(jìn)行，因而難以實(shí)現(xiàn)高精度、寬范圍的測(cè)量［11］。

同時(shí)，在實(shí)際應(yīng)用過(guò)程中發(fā)現(xiàn)，環(huán)境中不僅存在著高斯噪聲，還存在著大量的具有脈沖沖擊特性的非高斯噪聲。研究結(jié)果表明［12］［13］，這種噪聲在時(shí)域上表現(xiàn)出大量的脈沖尖峰特征，在頻域上顯示出更為厚重的拖尾現(xiàn)象。Nikias 指出［12］，該類噪聲可利用α 穩(wěn)定分布過(guò)程進(jìn)行表征。特別是當(dāng)特征指數(shù) α=2 時(shí)，α 穩(wěn)定分布噪聲（以下簡(jiǎn)稱α 噪聲）退變?yōu)閭鹘y(tǒng)高斯噪聲。因此，高斯噪聲是α 噪聲的一個(gè)特例，α 噪聲具有更為寬泛的適用范圍。但是，α 穩(wěn)定分布不存在有限的二階矩，甚至當(dāng) α<1 時(shí)，連均值的運(yùn)算都不能進(jìn)行［13］。這直接導(dǎo)致傳統(tǒng)的那些基于二階矩或高階累積量的超聲波測(cè)量方法估計(jì)性能退化甚至失效。

針對(duì)上述問(wèn)題，本文提出了一種基于分?jǐn)?shù)低階矩（Fractional Lower Order Moment，F(xiàn)LOM）的雙相測(cè)量方法。首先，本文利用FLOM 算子對(duì)α 與高斯混合噪聲進(jìn)行抑制，提升了抑制噪聲的能力；然后，將風(fēng)矢量測(cè)量問(wèn)題轉(zhuǎn)化為相位估計(jì)問(wèn)題，進(jìn)而有效地?cái)U(kuò)大了算法的測(cè)量范圍，利用參考信號(hào)的正交性提出了一種基于FLOM 的雙相估計(jì)方法（以下簡(jiǎn)稱分?jǐn)?shù)低階雙相估計(jì)方法）。該方法可消除超聲波信號(hào)傳播過(guò)程中的幅度變化對(duì)相位估計(jì)的影響，同時(shí)實(shí)現(xiàn)了高精度、寬范圍的風(fēng)矢量測(cè)量。

2 模型建立

2.1 超聲波測(cè)量模型

如圖1 所示，分別在南北方向以及東西方向各布置一組超聲波換能器，即S1，S2，W1，W2，其中S1 與S2，W1 與W2 之間距離相等，且同時(shí)具備收發(fā)功能，并采用相對(duì)時(shí)差法對(duì)風(fēng)矢量進(jìn)行測(cè)量。

圖1 風(fēng)矢量測(cè)量模型Fig.1 Model of wind vector measurement

以東西方向?yàn)槔?，設(shè)由東到西（W2→W1）為順風(fēng)。當(dāng)W2 發(fā)射超聲波，W1 接收時(shí)，超聲波在W2→W1 方向上傳播的真實(shí)速度為：

其中：C為聲速，VW為風(fēng)速在東西方向上的分量。

設(shè)W2 與W1 之間的距離為L(zhǎng)，則超聲波在W2→W1 方向上傳播的時(shí)間為：

當(dāng)W1 發(fā)射超聲波，W2 接收時(shí)，超聲波在W1→W2 方向上傳播的真實(shí)速度為：

則超聲波在W1→W2 方向上傳播的時(shí)間為：

由式（2）與式（4）可求得風(fēng)速在東西方向上分量為：

由式（5）可知，采用相對(duì)時(shí)差法可消除環(huán)境因素（溫度、濕度等）對(duì)聲速的影響。

同理，設(shè)由南到北（S2→S1）為順風(fēng)。則風(fēng)速在南北方向上的速度分量為：

其中：TS2S1為超聲波沿S2→S1 方向傳播的時(shí)間，TS1S2為超聲波沿S1→S2 方向傳播的時(shí)間。

由于W1-W2 與S1-S2 之間正交，由式（5）與式（6）可得風(fēng)速V與風(fēng)向角θ分別為：

2.2 含噪信號(hào)模型

設(shè)超聲波的發(fā)射信號(hào)為：

其中：A為參考信號(hào)的幅值，ω=2πf為參考信號(hào)的頻率，f≥20 KHz 為中心頻率。

設(shè)接收到的含噪超聲信號(hào)y（t）為：

其中：s(t)為經(jīng)過(guò)延遲的超聲波信號(hào)，B為超聲波傳播后經(jīng)過(guò)變化的幅值，Δt為超聲波傳播的延遲時(shí)間，n(t)為加性噪聲。

不失一般性，假設(shè)噪聲n(t)=nα(t)+ng(t)，其中nα(t)服從α穩(wěn)定分布，為 SαS 噪聲，ng(t)為高斯白噪聲，且ng(t)，nα(t)與x(t)，s(t)之間統(tǒng)計(jì)獨(dú)立。

傳統(tǒng)的基于時(shí)間延遲估計(jì)的超聲波風(fēng)矢量測(cè)量方法是通過(guò)直接估計(jì)超聲波傳播的延遲時(shí)間Δt來(lái)計(jì)算出風(fēng)速與風(fēng)向。

注意到，由于

其中，φ=ωΔt為接收和發(fā)射超聲信號(hào)的相位差。

因此，與傳統(tǒng)的時(shí)間延遲估計(jì)方法相比，采用相位估計(jì)方法來(lái)估計(jì)φ等價(jià)于將Δt放大了ω倍，在相同的誤差范圍下，能夠獲得更高的估計(jì)精度，進(jìn)而使得風(fēng)矢量具有更高的測(cè)量精度和更寬的測(cè)量范圍。

2.3 α 穩(wěn)定分布

α穩(wěn)定分布沒(méi)有閉式的概率密度解析式，一般通過(guò)特征函數(shù)進(jìn)行表征［14］：

式（12）可記為：Sα（γ，β，a）。其中0<α≤2為特征指數(shù)，-1≤β≤ 1 為偏斜參數(shù)，γ≥ 0 為分散系數(shù)，-∞

3 基于分?jǐn)?shù)低階雙相估計(jì)方法的超聲波風(fēng)矢量測(cè)量

根據(jù) FLOM 算子，由式（9）與式（11）可得超聲波含噪信號(hào)的分?jǐn)?shù)低階互相關(guān)函數(shù)為［15］：

其中 0

將式（9）與式（11）代入式（15）可得：

由廣義二項(xiàng)展開(kāi)式可知［16］，式（16）可轉(zhuǎn)化為：

由前文假設(shè)，式（17）可以化簡(jiǎn)為：

由式（18）可知，利用 FLOM 算子，可將附加噪聲n(t)充分抑制。為了獲取超聲波信號(hào)的相位差，將式（18）進(jìn)一步化簡(jiǎn)可得：

在傳統(tǒng)的相位估計(jì)方法中，均未考慮信號(hào)的幅值變化，即幅值A(chǔ)=B。然而，在超聲波測(cè)量過(guò)程中，接收信號(hào)的幅值會(huì)受到傳播路徑、超聲波換能器的性能等因素而發(fā)生變化，即幅值A(chǔ)≠B。如果仍認(rèn)為A=B，在對(duì)φ進(jìn)行估計(jì)時(shí)則會(huì)產(chǎn)生誤差，進(jìn)而對(duì)風(fēng)速風(fēng)向測(cè)量的精確性產(chǎn)生較大影響。因此，為了消除信號(hào)幅值變化帶來(lái)的影響，本文在分?jǐn)?shù)低階互相關(guān)函數(shù)的基礎(chǔ)上提出了分?jǐn)?shù)低階雙相估計(jì)方法。

將式（9）的發(fā)射參考信號(hào)進(jìn)行90°相移，可得

將式（20）與式（11）代入式（15）可得：

利用前文所述式（16）到式（19）的推導(dǎo)過(guò)程可得：

由式（22）與式（19）聯(lián)立可得：

對(duì)比式（23）與式（19）可知，本文所提的分?jǐn)?shù)低階雙相估計(jì)方法僅與相位差φ有關(guān)，有效的消除了超聲波傳播過(guò)程中幅值變化對(duì)相位估計(jì)的影響，極大提高了相位估計(jì)的精度。

為了確定式（24）中的n，以零風(fēng)速條件為參考，由式（11）可得：

其中，φ0為零風(fēng)速條件下的相位差。

設(shè)W1~W2 與S1~S2 方向上所測(cè)的最大風(fēng)速為v，分別取順風(fēng)與逆風(fēng)方向時(shí)，有：

其中，φv為最大風(fēng)速條件下的相位差。

則式（25）與式（26）需滿足如下條件：

化簡(jiǎn)式（27）可得：

由式（28）可知，在0~v風(fēng)速段內(nèi)滿足條件的f，L可使所測(cè)得的相位差φ僅在φ0的 ±2π內(nèi)波動(dòng)，此時(shí)，相位差φ在式（24）中的n值與零風(fēng)速條件下的n值相同，避免了因選取的超聲波中心頻率過(guò)高、換能器之間距離過(guò)長(zhǎng)等因素引起的n值不確定的情況。其中，零風(fēng)速條件下的n值可由式（25）確定，由式（24）與零風(fēng)速下的n值即可確定φ。

由渡越時(shí)間與相位差關(guān)系式，風(fēng)速在W1-W2 與S1-S2 方向上的分量分別為：

其中：φW2W1為W2→W1 方向上的相位差，φW1W2為W1→W2 方向上的相位差；

其中：φS2S1為S2→S1 方向上的相位差，φS1S2為S1→S2 方向上的相位差。

將式（29）與式（30）代入式（7）與式（8）即可計(jì)算出風(fēng)速V與風(fēng)向角θ。

4 仿真分析與實(shí)際測(cè)試

4.1 仿真分析

實(shí)驗(yàn)1 驗(yàn)證分?jǐn)?shù)低階雙相估計(jì)方法對(duì)α 與高斯混合噪聲的抑制能力。

實(shí)驗(yàn)條件：超聲波的中心頻率f=25.0 KHz，幅度系數(shù)A=1，幅度系數(shù)B=0.6，超聲波傳播的距離L=0.07 m，聲速C=340 m/s；附加噪聲nα(t)為 SαS 噪聲，它可按文獻(xiàn)［17］產(chǎn)生，特征指數(shù)α=1.2，其信噪比采用如下廣義信噪比定義：

取GSNR=2.5 dB。附加噪聲ng(t)為高斯噪聲，可根據(jù)文獻(xiàn)［17］產(chǎn)生，其信噪比采用如下定義：

其中：PS為信號(hào)的有效功率，PN為噪聲的有效功率。取SNg R=2.5 dB。FLOM 算子的階次p=0.5，采樣頻率取10 MHz，數(shù)據(jù)長(zhǎng)度取10 000。測(cè)量風(fēng)速V=24.16 m s，風(fēng)向角θ=65.56°。其中f，L，V，θ滿足式（28）。

實(shí)驗(yàn)方法：使用本文方法與傳統(tǒng)基于二階矩的時(shí)延估計(jì)方法分別對(duì)風(fēng)速以及風(fēng)向角進(jìn)行估計(jì)，進(jìn)行100 次蒙特卡洛實(shí)驗(yàn)。其中基于二階矩的時(shí)延估計(jì)方法可由文獻(xiàn)［18］定義。

實(shí)驗(yàn)結(jié)果：圖2 為基于二階矩的時(shí)延估計(jì)方法和本文分?jǐn)?shù)低階雙相估計(jì)方法在α與高斯混合噪聲下超聲波測(cè)風(fēng)性能的仿真結(jié)果。

由圖2 可知，由于α與高斯混合噪聲的影響，傳統(tǒng)的二階矩時(shí)延估計(jì)方法將會(huì)在超聲波進(jìn)行風(fēng)矢量測(cè)量時(shí)發(fā)生嚴(yán)重退化。而本文的分?jǐn)?shù)低階雙相估計(jì)方法仍然能夠較為準(zhǔn)確地利用超聲波進(jìn)行風(fēng)矢量的測(cè)量，所測(cè)得的風(fēng)速的均方根誤差為0.009 1，風(fēng)向角的均方根誤差為0.13。

圖2 實(shí)驗(yàn)1 估計(jì)結(jié)果Fig.2 The estimation results of experiment 1

實(shí)驗(yàn)2 驗(yàn)證分?jǐn)?shù)低階雙相估計(jì)方法的高精度測(cè)量性能。

實(shí)驗(yàn)條件：同實(shí)驗(yàn)1。

實(shí)驗(yàn)方法：使用本文方法與基于FLOM 的時(shí)延估計(jì)方法分別對(duì)風(fēng)速以及風(fēng)向角進(jìn)行估計(jì)，進(jìn)行100 次蒙特卡洛實(shí)驗(yàn)。其中基于FLOM 的時(shí)延估計(jì)方法可由文獻(xiàn)［18］定義。

實(shí)驗(yàn)結(jié)果：圖3 為基于FLOM 的時(shí)延估計(jì)方法和本文分?jǐn)?shù)低階雙相估計(jì)方法在α與高斯混合噪聲下超聲波測(cè)風(fēng)性能的仿真結(jié)果。

由圖3 可知，在α與高斯混合噪聲的影響下，基于分?jǐn)?shù)低階時(shí)延估計(jì)的超聲波風(fēng)矢量測(cè)量方法所測(cè)得的風(fēng)速的均方根誤差為0.027，風(fēng)向角的均方根誤差為0.361，本文方法所測(cè)得的風(fēng)速的均方根誤差為0.009 1，風(fēng)向角的均方根誤差為0.13。兩種方法均具有較好的測(cè)量性能，但是本文分?jǐn)?shù)低階雙相測(cè)量方法精確度更高。

圖3 實(shí)驗(yàn)2 估計(jì)結(jié)果Fig.3 The estimation results of experiment 2

實(shí)驗(yàn)3 驗(yàn)證分?jǐn)?shù)低階雙相估計(jì)方法的寬范圍測(cè)量性能。

實(shí)驗(yàn)條件：在實(shí)驗(yàn)1 的基礎(chǔ)上將測(cè)量風(fēng)速提升到V=70.71 m s，風(fēng)向角θ=45°。

實(shí)驗(yàn)方法：使用本文方法與基于FLOM 的時(shí)延估計(jì)方法分別對(duì)風(fēng)速以及風(fēng)向角進(jìn)行估計(jì)，進(jìn)行100 次蒙特卡洛實(shí)驗(yàn)。

實(shí)驗(yàn)結(jié)果：圖4 為基于FLOM 的時(shí)延估計(jì)方法和本文分?jǐn)?shù)低階雙相估計(jì)方法在α與高斯混合噪聲下超聲波測(cè)風(fēng)性能的仿真結(jié)果。

圖4 實(shí)驗(yàn)3 估計(jì)結(jié)果Fig.4 The estimation results of experiment 3

由圖4 可知，在α 與高斯混合噪聲的影響下，當(dāng)測(cè)量風(fēng)速提高到V=70.71 m s，風(fēng)向角為θ=45°時(shí)，基于分?jǐn)?shù)低階時(shí)延估計(jì)的超聲波風(fēng)矢量測(cè)量方法已經(jīng)失效，而本文的分?jǐn)?shù)低階雙相測(cè)量方法仍然能夠較為準(zhǔn)確地利用超聲波進(jìn)行風(fēng)矢量的測(cè)量，所測(cè)得的風(fēng)速的均方根誤差為0.083，風(fēng)向角的均方根誤差為0.363。

對(duì)比實(shí)驗(yàn)2 與實(shí)驗(yàn)3 可知，本文分?jǐn)?shù)低階雙相估計(jì)方法較基于分?jǐn)?shù)低階時(shí)延估計(jì)的超聲波風(fēng)矢量測(cè)量方法的測(cè)量范圍寬。

實(shí)驗(yàn)4 驗(yàn)證本文分?jǐn)?shù)低階雙相估計(jì)方法在不同信噪比條件下的超聲波測(cè)風(fēng)性能。

實(shí)驗(yàn)條件：在實(shí)驗(yàn)1 條件下，使附加噪聲的信噪比SNR=GSNR+SNg R以5 dB 間隔從-10 dB 到 10 dB 增加。

實(shí)驗(yàn)方法：在每一信噪比下，采用本文分?jǐn)?shù)低階雙相估計(jì)方法進(jìn)行100 次蒙特卡洛實(shí)驗(yàn)。將實(shí)驗(yàn)的結(jié)果進(jìn)行均方根誤差（RMSE）［19］運(yùn)算：

其中：Vr為真實(shí)風(fēng)速，θr為真實(shí)風(fēng)向角。

實(shí)驗(yàn)結(jié)果：圖5 為在不同信噪比下，采用本文分?jǐn)?shù)低階雙相估計(jì)方法進(jìn)行超聲波風(fēng)矢量得到的風(fēng)速與風(fēng)向角的均方根誤差曲線。

圖5 風(fēng)速與風(fēng)向角的均方根誤差曲線Fig.5 RMSE curve of wind speed and wind direction angle

由圖5 可知，隨著α 與高斯混合噪聲信噪比（SNR）的逐漸增大，測(cè)得的風(fēng)速與風(fēng)向角的均方根誤差（RMSE）逐漸減小。即使在SNR=-10 dB 的條件下，本文方法在進(jìn)行超聲波風(fēng)矢量測(cè)量時(shí)仍有很高的精度，充分體現(xiàn)了本文方法對(duì)高斯與α 混合噪聲的強(qiáng)抑制能力。

4.2 實(shí)際測(cè)試

將本文方法應(yīng)用在超聲波風(fēng)矢量測(cè)量實(shí)驗(yàn)臺(tái)上進(jìn)行測(cè)試，實(shí)驗(yàn)平臺(tái)如圖6 所示。

圖6 超聲波測(cè)風(fēng)平臺(tái)Fig.6 Ultrasonic wind measuring platform

實(shí)驗(yàn)方法：在長(zhǎng)春氣象所EDE1-5 型風(fēng)洞中，在風(fēng)速V=71.06 m s，風(fēng)向角θ=50.71°的環(huán)境下進(jìn)行測(cè)試。其中L=0.07 m，f=25.0 KHz。采用在測(cè)量環(huán)境周圍放置多個(gè)交流接觸器，使其進(jìn)行反復(fù)開(kāi)關(guān)動(dòng)作產(chǎn)生高斯與α混合噪聲。

實(shí)驗(yàn)結(jié)果：圖7 為在超聲波風(fēng)矢量實(shí)測(cè)平臺(tái)所得的實(shí)驗(yàn)結(jié)果。

圖7 風(fēng)矢量實(shí)測(cè)圖Fig.7 Actual wind vector test map

由圖7 可知，在α與高斯混合噪聲的干擾下，采用本文分?jǐn)?shù)低階雙相估計(jì)方法所測(cè)得的均方根誤差為0.104，風(fēng)向角的均方根誤差為0.54?？梢?jiàn)本文提出的分?jǐn)?shù)低階雙相估計(jì)方法在實(shí)物平臺(tái)實(shí)測(cè)風(fēng)矢量時(shí)具有較強(qiáng)的泛化能力。

5 結(jié) 論

本文針對(duì)α 與高斯混合噪聲下超聲波風(fēng)矢量測(cè)量問(wèn)題，首先利用FLOM 算子對(duì)α與高斯混合噪聲進(jìn)行抑制，彌補(bǔ)了傳統(tǒng)的二階矩和高階累積量方法在超聲波風(fēng)矢量測(cè)量中退化失效的不足。然后在此基礎(chǔ)上，將時(shí)延估計(jì)方法轉(zhuǎn)化為相位估計(jì)方法，并利用參考信號(hào)的正交性，提出了一種分?jǐn)?shù)低階雙相估計(jì)方法，解決了超聲波傳播過(guò)程中的幅值變化降低風(fēng)矢量測(cè)量精度的問(wèn)題，有效地?cái)U(kuò)大了風(fēng)矢量的測(cè)量范圍。仿真與實(shí)測(cè)結(jié)果表明：

（1）與傳統(tǒng)的二階矩時(shí)延估計(jì)方法相比，分?jǐn)?shù)低階雙相測(cè)量方法對(duì)高斯與α混合噪聲具有強(qiáng)抑制能力。與基于FLOM 的時(shí)延估計(jì)方法相比，分?jǐn)?shù)低階雙相測(cè)量方法在中低風(fēng)速段具有更高的精度，而在高風(fēng)速段，基于FLOM 的時(shí)延估計(jì)方法已經(jīng)失效，但是分?jǐn)?shù)低階雙相測(cè)量方法仍然具有良好的測(cè)風(fēng)能力。

（2）在信噪比為-10 dB 的條件下，風(fēng)速的均方根誤差小于1.5 m/s，風(fēng)向角的均方根誤差小于2°，仍然具有較高的精度。

（3）在實(shí)測(cè)平臺(tái)上，本文方法在風(fēng)速V=70.71 m s，風(fēng)向角θ=45°的環(huán)境下所測(cè)得的風(fēng)速的均方根誤差為0.104，風(fēng)向角的均方根誤差為0.54，具有較強(qiáng)的泛化能力。