王 婧，劉紀(jì)元 ，焦學(xué)峰，張 皓

（1.沈陽航空航天大學(xué) 電子信息工程學(xué)院，遼寧 沈陽 110136；2.中國科學(xué)院聲學(xué)研究所綜合聲納實(shí)驗(yàn)室 北京 100080；3.北京航瑞博泰科技有限公司 北京 100102）

在注重節(jié)省能源，提高經(jīng)濟(jì)效益和產(chǎn)品質(zhì)量的今天，流量計量與測試的重要更加突出并為越來越多的人所認(rèn)識。特別是隨著現(xiàn)代工業(yè)生產(chǎn)的飛速發(fā)展，人們對流量測量的要求越來越高，對流量測量技術(shù)和儀表的研究和開發(fā)不斷深入，流量測量方法儀表的種類也越來越多。在一些大型引水場合，如：城市供水引水渠，火電廠引水和排水渠，污水治理流入和排放渠，工礦企業(yè)水排放以及水利工程和農(nóng)業(yè)灌溉用渠道，明渠流量計發(fā)揮極其重要的作用。而超聲波明渠流量計中影響流量的一個重要因素就是水速的測量及其精度，因此合理獲取精確的水流速顯得至關(guān)重要。

圖1 多普勒原理圖Fig.1 Doppler schematic

1 水聲測速的基本原理

水聲測速的原理圖如圖1所示。

在換能器收發(fā)共置的情況下，聲波經(jīng)發(fā)射后，測量聲波被水中移動的懸浮物質(zhì)反射，此時回波信號與發(fā)射信號就會產(chǎn)生頻率差，即多普勒頻移，利用頻移就可以計算流速。

在一個測量周期內(nèi)，首先換能器以發(fā)射頻率fc發(fā)射一束超聲波，這束波被一個懸浮質(zhì)以與發(fā)射頻率不同的頻率接收到（或反射），這一過程就產(chǎn)生了多普勒頻移，然后波束又一次反射回?fù)Q能器，再一次產(chǎn)生一個多普勒頻移。兩個多普勒頻移存在于這個交換過程中，前一個是信號通過換能器傳到懸浮質(zhì)，第二個是信號從懸浮質(zhì)反射回?fù)Q能器。

應(yīng)用多普勒原理，同時考慮發(fā)射角α可以推導(dǎo)以下結(jié)論：

式中，c為聲波在水中的傳播速度，常溫下大約1 500 m/s；v水為水流流速；fc為聲源發(fā)射信號頻率；α為發(fā)射信號與水平面的夾角。

在測量的過程中，由于v水＜＜c，公式簡化為：

因此，要計算多普勒頻移，只要將所接收到的回波做FFT后與水靜止時的FFT做差，就可以得到頻移f頻移，從而得到水的流速V水。

2 三種方法概述

2.1 疊加法

顧名思義，疊加法就是取多次回波頻譜將其相加，這樣做的好處是，經(jīng)過多次疊加，一些隨機(jī)噪聲的幅值隨著疊加次數(shù)的增多增長很緩慢，但信號的頻譜卻增長明顯，因此信號越來越突出，而噪聲信號則越來越不明顯了[1]，證明如下。

設(shè)每次測得頻率為fi，每次測量值是互相獨(dú)立的隨機(jī)量，設(shè)其方δ2差均為，進(jìn)行N次獨(dú)立測量后的樣本均值為

方差為

可見疊加的次數(shù)越多，效果越好，但是為能夠反映當(dāng)前水速值，測量時間不宜過長（次數(shù)不宜太多）。

2.2 譜減法

與語音信號的“譜減法”[2]不同的是，“回波譜減法”不是減去“寂靜段”的噪聲功率譜，而是帶信號去噪。步驟如下：

1）首先發(fā)射換能器不發(fā)射信號，觀察回波頻譜，確定干擾的大概范圍R。

2）調(diào)整水流速，使其對應(yīng)的頻譜在R之外，保存R內(nèi)的數(shù)值。本步驟執(zhí)行多次。

3）根據(jù)上一步驟整理數(shù)值，看干擾R內(nèi)的數(shù)值與主頻值的比例關(guān)系及主頻與真實(shí)多普勒頻率值的關(guān)系。

4）在接收的回波信號中依據(jù)比例減去干擾。

2.3 全相位頻譜校正法

通過式（4）可以看出，水速的獲得是通過回波的頻率獲得的，要想提高測速的精確度，就要知道回波的精確頻率。如果真正的頻率在兩根譜線之間，則誤差最大，為半個頻率分辨率[3]。因此頻譜校正對于提高測量精度是很必要的。傳統(tǒng)FFT結(jié)果可以通過一些算法實(shí)現(xiàn)頻譜校正，如能量重心法、比值法等。但以上的頻率估計都是在傳統(tǒng)的FFT架構(gòu)下進(jìn)行的，因而傳統(tǒng)FFT固有的頻譜泄露效應(yīng)無疑會在很大程度上影響這些校正法的精度。全相位FFT算法（apFFT）具有初始相位不變和有效防止頻譜泄露的特性[4]。

在圖 2 中，ω*代表信號的數(shù)字角頻率，k*=[ω*/Δω]（“[]”代表四舍五入），Δω 是角頻率分辨率ω^*，A^，及θ^0為校正后的角頻率、幅值及相位。

圖2 FFT/apFFT綜合相位差法的頻譜校正流程Fig.2 The spectra calibrating process of FFT/apFFT method

由圖2可以看出，F(xiàn)FT/apFFT綜合相位差法需要對序列進(jìn)行2次譜分析，一次傳統(tǒng)FFT分析，一次全相位FFT譜分析。雙窗全相位FFT過程，即對數(shù)據(jù)進(jìn)行全相位預(yù)處理，然后進(jìn)行傳統(tǒng)的FFT運(yùn)算。過程如下：

1）構(gòu)成一個N點(diǎn)的漢寧窗；

2）漢寧窗對自己求卷積，得到2N－1點(diǎn)的卷積窗；

3）求2N－1點(diǎn)的卷積窗的和；

4）將卷積窗的每一項(xiàng)除以卷積窗的和，得到2N－1點(diǎn)的歸一化卷積窗；

5）將數(shù)據(jù)的1到2N－1項(xiàng)和歸一化卷積窗相乘，得到加窗的2N－1項(xiàng)；

6）將第 1項(xiàng)和 N+1項(xiàng)，第 2項(xiàng)和 N+2項(xiàng) ...第 N－1項(xiàng)和第2N－1項(xiàng)相加，得到經(jīng)過全相預(yù)處理的N點(diǎn)序列；

7）對6）中的經(jīng)過全相位處理的N點(diǎn)序列做FFT運(yùn)算；

以下簡單介紹FFT/apFFT。

已知單頻復(fù)指數(shù)信號為

信號的數(shù)字頻率ω*表示為β倍頻率間隔2π/N的形式（β可以是小數(shù)）。

其歸一化的FFTapFFT譜：

其歸一化的apFFT譜：

FFT/apFFT校正法需要2N－1個樣點(diǎn)。對后N點(diǎn)作FFT，對全部的2N－1點(diǎn)作apFFT。對頻率值β取整或從振幅譜可確定譜峰位置位于 k0=[β]處，[]表示取整，ang（.）表示取相操作。

由式（9）可知，在k0的apFFT的初始相位值φ0，則相位校正值：

由式（8）和（9）可知，在k0的FFT的相位值減去apFFT的相位值可得校正后的頻率值β：

由式（7）和（8）可知，在的FFT的振幅值平方除以apFFT振幅即得校正后的振幅A:

3 結(jié)果分析

3.1 未加任何處理的回波波形

圖3為未加任何處理的靜水中的回波信號頻譜圖。

圖3 靜水中的回波信號頻譜圖Fig.3 Echo signal spectrum of still water

經(jīng)過采樣后的回波信號的余頻所在的譜線數(shù)在409.6≈410根處。

圖3以及下面的圖4、圖5、圖6的橫軸代表的是譜線根數(shù)，單位為頻率分辨率，即fs/N Hz，N為FFT點(diǎn)數(shù)。縱軸是對采集到的數(shù)字信號做FFT后取絕對值。單位為10/216V。

在上圖的回波波形中可以看出，回波有噪聲干擾。

3.2 疊加后的頻譜圖

疊加后的靜水及水流動的回波頻譜圖如圖4所示。

圖4 疊加后的靜水及水流動的回波頻譜圖Fig.4 Superposed echo spectrum of still water and flowing water

由上圖可以看出，疊加后的頻譜圖輪廓清晰，峰值明顯，噪聲與信號的頻譜比較，微乎其微。但是此時的主頻處410根譜線的峰值仍然大于真實(shí)譜線處的峰值 （410根右側(cè)的那個峰值）。

3.3 疊加+譜減法后的頻譜圖

圖5為靜水中疊加及譜減法后的回波頻譜圖。

圖5 靜水中疊加及譜減法后的回波頻譜圖Fig.5 Echo spectrum after superposition and spectral subtraction of still water

從圖6中可以看到，頻譜圖稍亂一些，這是由于系數(shù)選取的緣故，因?yàn)樗x取的系數(shù)是多次測量取的平均值。但是回波的波形不是固定不變的，所以經(jīng)過處理后的頻譜圖為當(dāng)前的狀態(tài)，但是這不會影響判斷結(jié)果。

圖6 水流動疊加及譜減法后的回波頻譜圖Fig.6 Echo spectrum after superposition and spectral subtraction of flowing water

現(xiàn)在已經(jīng)可以看出主頻處410根譜線的峰值不再大于真實(shí)值譜線處的峰值了。當(dāng)測量真實(shí)結(jié)果在R譜線范圍內(nèi)，譜減法后的頻譜也不會影響真實(shí)值的判斷。因此此種辦法有一定效果[6-7]。

3.4 FFT/apFFT校正前后對比

發(fā)射的信號頻率為：1 MHz；信號幅值：0.8 V；初相：180；采樣頻率：25.51 kHz；FFT/apFFT分析校正后的頻譜如圖7所示，表1列出了校正前后數(shù)據(jù)（靜水）。

在圖7中，橫坐標(biāo)k表示的是譜線根數(shù)，單位為頻率分辨率，即fs/N Hz，N為FFT點(diǎn)數(shù)?？v坐標(biāo)y1表示的是加窗后的FFT振幅譜，單位為V；y2表示的是振幅譜，單位為V；y3代表的是apFFT的相位譜，單位為（°）；y4表示的是歸一化的相位偏離校正值，單位為Δf Hz，Δf代表頻率分辨率；y5表示的是振幅校正值，單位為V[8-9]。

注意這里是欠采樣，所以回波頻率指的都是余頻。為了考核此種方法的可行性及穩(wěn)定性，在常溫下做了多組實(shí)驗(yàn)來驗(yàn)證。調(diào)節(jié)水速從而使回波產(chǎn)生不同的頻率，通過綜合應(yīng)用方法對回波頻率進(jìn)行校正（噪聲的幅值大約是信號幅值的0.1倍）。整理后的結(jié)果如表2所示。其中公式（4）中 α 的取 30°。

圖7 FFT/apFFT頻譜校正Fig.7 FFT/apFFTspectrum correction

表1 FFT/apFFT頻譜校正結(jié)果前后數(shù)值Tab.1 Results before and after FFT/apFFT correction

由于測速關(guān)心的主要問題是回波的頻率，故表2里只對頻率的校正前后做了對比。此處，相對誤差=（測量值-理論理）/理論值。

由式（4）可知，多普勒頻移的相對誤差就是水速的相對誤差。因此在計算水速的相對誤差時，可以通過計算多普勒頻移的相對誤差來獲得水速的相對誤差，減少計算量。

以第1組數(shù)據(jù)為例：

f頻移=5 159.735-5 102=57.735 Hz

由表中的數(shù)據(jù)可以看出，在有噪聲的情況下，校正后的相對誤差比校正前的相對誤差減低了很多，相對誤差基本降低了一個數(shù)量級。此方法在水速測量上對提高測度精度是起到了一定的作用。

表2 FFT/apFFT頻譜校正方法測試數(shù)據(jù)Tab.2 Test data of FFT/apFFT correction

4 結(jié) 論

本文在設(shè)計中將“疊加法”、“譜減法”及“全相位頻譜校正法”結(jié)合應(yīng)用于水速測量中。這些處理方法的結(jié)合，一方面可以降低噪聲干擾，突出信號，另一方面降低了主頻對其他頻率的干擾。特別是全相位頻譜校正方法的引用，可以避免

離散頻譜引起的誤差，能夠在有噪聲的情況下，將誤差降低一個數(shù)量級。

[1] 高晉占.微弱信號檢測[M].北京:清華大學(xué)出版社，2002.

[2] 楊行峻，遲惠生.語音信號數(shù)字處理[M].北京:電子工業(yè)出版社，1995.

[3] 丁康，謝明，楊志堅(jiān).頻譜分析校正理論與技術(shù)[M].北京:科學(xué)出版社，2008.

[4]王兆華，黃翔東.數(shù)字信號全相位頻譜分析與濾波技術(shù)[M].北京:電子工業(yè)出版社，2009.

[5] 黃翔東，王兆華.基于全相位頻譜分析的相位差頻譜校正法[J].電子與信息學(xué)報，2008，30（2）：293-297.HANG Xiang-dong，WANG Zhao-hua.Phasedifference correcting spectrum method based on all-phase spectrum analysis[J].Journal of Electronics&Information Technology，2008，30（2）:293-297.

[6] 李啟虎.聲納信號處理引論[M].北京:海洋出版社，1985.

[7] 俞小鼎.多普勒天氣雷達(dá)原理與業(yè)務(wù)應(yīng)用[M].北京:氣象出版社，2006.

[8] 張賢達(dá).現(xiàn)代信號處理[M].北京：清華大學(xué)出版社，2002.

[9] 劉波，文忠，等.MATLAB信號處理[M].北京：電子工業(yè)出版社，2006.