楊永壽 王忠康
(杭州瑞利海洋裝備有限公司,杭州,310023)
多普勒聲吶是一種測量自然和人工水道中水流速度及剖面流量的聲學(xué)設(shè)備,由于具備不干擾水流場、測量精度高和測量成本低等優(yōu)勢,已經(jīng)廣泛應(yīng)用于水文測量、海洋研究、水資源勘探和水下導(dǎo)航等諸多領(lǐng)域[1-3]。
寬帶聲學(xué)多普勒測流信號是一隨機過程,通常具有有限的平均功率,可以用功率譜密度(Power Spectral Density,PSD)相關(guān)方法來描述、分析和處理?;诠β首V矩估計的測流信號處理方法因具有精度高、應(yīng)用廣等優(yōu)勢而逐漸獲得研究人員的重視。文獻[4]基于準平穩(wěn)和高斯隨機信號模型,利用單延遲自相關(guān)和全譜估計等譜矩估計方法研究了大氣探測雷達中的參數(shù)估計和干擾抑制問題。文獻[5]利用協(xié)方差法推導(dǎo)了偽噪聲序列一階譜矩和二階譜矩的均值和方差表達式,指出一階譜矩的均值是多普勒頻移的無偏估計,方差與子序列的相關(guān)系數(shù)、信噪比和時寬帶寬積有關(guān)。文獻[6]提出了一種基于譜密度函數(shù)理論的譜矩函數(shù)分析方法,旨在同時計算電荷密度和譜函數(shù),并提出了一種利用前 4個譜矩計算譜函數(shù)的有效算法。文獻[7]針對具有加性復(fù)高斯噪聲的隨機信號,基于譜矩估計原理推導(dǎo)了其頻率和譜峰寬度的均值和方差的漸進估計表達式,其方差估計的主導(dǎo)項與Cramer-Rao下界相同。文獻[8]提出了一種利用自相關(guān)函數(shù)計算高階譜矩的算法,用于分析非高斯功率譜天氣雷達觀測信號。與離散傅里葉變換(Discrete Fourier Transform,DFT)和脈沖對方法相比,該算法具有計算量小、準確度高等優(yōu)勢。文獻[9]提出了一種利用隨機信號的協(xié)方差樣本估計任意階譜矩的算法,并利用該算法計算了極地?zé)o線電信號經(jīng)電離層衰減所得接收信號的譜矩,驗證了算法的有效性。文獻[10]利用四階譜矩估計和深度學(xué)習(xí)研究了利用腦電波信號進行情感檢測的方法,譜矩估計用于檢測腦電信號的非線性和非高斯性。
現(xiàn)有的底速估計方法(如匹配濾波法)通過優(yōu)選能量最強或者與發(fā)射波形相關(guān)性最大的底回波分段得到底速。然而,能量最強、相關(guān)性最大回波分段的速度估計誤差卻不一定是最小的。受旁瓣干擾的影響,底回波前段的信噪比有所下降。在底速估計時,如能避開受污染的回波段則可提高底速估計的精度。本文所述寬帶回波是指寬帶編碼發(fā)射信號的散射回波經(jīng)數(shù)字解調(diào)后所得信號。
本文從譜矩估計的原理出發(fā),推導(dǎo)了寬帶回波的歸一化一階及二階譜矩表達式,據(jù)此設(shè)計了一種基于互相關(guān)系數(shù)的底速估計方法和一種基于信噪比的測流數(shù)據(jù)統(tǒng)計方法,兩種方法可顯著減小底速和流速的估計誤差。
在多普勒測流中,主波束通常傾斜向下,波束軸線與垂向成 20°~30°夾角。現(xiàn)有換能器存在與主波束成 30°~40°夾角的旁瓣。在主波束沒有達到底部之前,旁瓣已經(jīng)到達底部。在一定深度條件下,較強的旁瓣底回波會淹沒主波束回波的接收,從而污染底回波信號。這種污染可能導(dǎo)致3種影響:(1)受污染底回波分段的能量增大;(2)局部底回波的信噪比下降;(3)底回波功率譜的譜峰展寬,因為旁瓣底回波與主瓣底回波中包含了不同的多普勒信息。
旁瓣干擾底回波的示意圖如圖1所示,圖中水底以下部分僅為示意,并非實際聲波傳播路徑。由圖1中的紅色弧線可知,旁瓣底回波先于主瓣底回波到達換能器,由圖1中的黑色弧線可知,旁瓣底回波污染了主瓣底回波的前段,進而導(dǎo)致現(xiàn)有匹配濾波方法底速估計誤差的增大。
圖1 旁瓣對底回波的干擾示意圖
假設(shè)寬帶聲學(xué)多普勒測流信號具有時間平穩(wěn)性和空間均勻性,對于任意非負可積多普勒測流回波的PSDS(ω),定義其k階矩為
寬帶測流中的多普勒頻移均值可由其 PSD的歸一化一階矩來估計[7]:
由于 PSD的計算復(fù)雜度較高且其噪聲抑制能力沒有自相關(guān)方法強,所以使用自相關(guān)方法進行參數(shù)估計更加實用。根據(jù) Wiener-Khinchin定理,自相關(guān)函數(shù)R(τ)與PSD的關(guān)系如下:
式中,τ為自相關(guān)延遲。分別求式(4)的一階和二階導(dǎo)數(shù),并令結(jié)果中的τ=0,得到PSD零階矩、一階矩和二階矩與自相關(guān)函數(shù)的關(guān)系式:
綜合式(2)、(3)和式(5),得到
接著推導(dǎo)自相關(guān)函數(shù)一階和二階導(dǎo)數(shù)的計算式。寬帶回波的復(fù)自相關(guān)函數(shù)可表示為極坐標形式:
式中,A(τ)為幅度函數(shù),φ(τ)為相位函數(shù)。對式(8)兩邊分別求一階和二階導(dǎo)數(shù),并利用A(τ)是τ的實偶函數(shù),φ(τ)是τ的實奇函數(shù)的特性可得
綜合式(6)、(7)和式(9),可得
對φ(τ)在τ=0處進行泰勒展開并取前兩項,得到φ'(0)≈[φ(τ)-φ(0)]/τ,且φ(0)=0;對A(τ)在τ=0 處進行泰勒展開并取前三項,得到A"(0)≈2[A(τ)-A(0)]/τ2;由式(8)可知A(τ)=|R(τ)|,A(0)=R(0)。綜合可得利用自相關(guān)函數(shù)計算PSD譜矩的表達式為
式(12)~(13)分別為多普勒頻移均值和方差的估計表達式。式(13)中,|R(τ)|表示時延為τ的兩段回波互相關(guān)的模值,R(0)等于回波中噪聲與信號的總功率,|R(τ)|/R(0)表征以上兩段回波的相關(guān)性。由式(13)可知,|R(τ)|/R(0)的值與頻率估計方差成反比。
根據(jù)以上原理,設(shè)計了一種基于互相關(guān)系數(shù)的底速估計方法和一種基于信噪比的測流數(shù)據(jù)統(tǒng)計方法,可減小底速和流速的估計誤差。
由第1節(jié)的分析可知,受旁瓣污染的影響,主瓣底回波前段的信噪比下降、功率譜峰展寬,其與后段未被污染底回波的相關(guān)性將顯著下降。|R(τ)|/R(0)可以用來表征兩段回波的相關(guān)性,本文在式(13)的基礎(chǔ)上設(shè)計了如式(14)所示的底回波優(yōu)選判據(jù)ηb,用于優(yōu)選受旁瓣干擾最小的底回波分段。
式中,R1(0)和R2(0)分別為參與復(fù)自相關(guān)運算的兩段回波信號樣本的自相關(guān)。底回波信號的ηb判據(jù)值越大,其前后兩段的相關(guān)性越強,表示受旁瓣污染的程度越輕。利用ηb判據(jù)可以優(yōu)選受旁瓣干擾最小的底回波段,進而提高底跟蹤速度估計的精度。
由于實際測流回波的信號功率一般較小,式(13)中等號右側(cè)的表達式對于回波信噪比的敏感度不高,不適合作為數(shù)據(jù)統(tǒng)計分析的判斷依據(jù)。本文在式(13)的基礎(chǔ)上設(shè)計了一種測流數(shù)據(jù)統(tǒng)計分析判據(jù)ηv。
式中,|R(τ)|/R(0)表征分層回波前后兩段信號的相關(guān)性,其值越大即ηv越大則分層回波的信噪比越高,流速估計的誤差也就越小。測流儀器實時計算水體分層回波的ηv,對ηv數(shù)值小于設(shè)定閾值的測流數(shù)據(jù)予以剔除,或替換為合格數(shù)據(jù)的均值。
利用外場試驗結(jié)果驗證基于互相關(guān)系數(shù)的底速估計方法和基于信噪比的測流數(shù)據(jù)統(tǒng)計方法的性能。外場試驗河段長約1 km,河面寬度約60 m,測船靜止試驗處的水深約為 3.3 m,測船運動試驗河段的平均深度約為 5.5 m。所用試驗儀器為一臺600 kHz寬帶聲學(xué)多普勒測流試驗系統(tǒng)。
試驗樣機波束1在測船靜止條件下的單次回波樣本如圖2(a)所示,此時底速的理論值為零。逐點滑動計算與發(fā)射脈沖等長的回波段的判據(jù)值ηb與底速估計誤差,并與匹配濾波法進行了對比,結(jié)果如圖2(b)所示。圖2(b)中紅色粗實線表示ηb判據(jù),灰色實線表示匹配濾波判據(jù),ηb判據(jù)的值基本與底速估計誤差呈反比,這與理論預(yù)測一致。紅色和灰色虛線分別表示本文底速估計方法和匹配濾波法根據(jù)各自判據(jù)最大值確定的底回波起點3.995 ms和3.92 ms,其對應(yīng)的底速估計誤差分別為-0.076 m/s和-0.222 m/s。本文方法的底速估計誤差相對減小了66%。
圖2 單個樣本的底回波優(yōu)選方法
分別利用本文方法和匹配濾波方法對測船靜止試驗中的100個樣本進行底回波優(yōu)選,結(jié)果如圖3所示。
圖3 多個樣本的底回波優(yōu)選結(jié)果分析
圖3(a)中紅色方形和黑色圓圈分別表示本文方法和匹配濾波法優(yōu)選的底回波時間起點,紅色和黑色實線分別表示對應(yīng)的最小二乘擬合結(jié)果。本文方法所選底回波起點總體上比匹配濾波方法延遲約0.065 ms。將圖3(a)中匹配濾波法所得底回波起始時刻減去互相關(guān)系數(shù)法所得底回波起始時刻,結(jié)果如圖 3(b)中黑色三角所示。圖 3(b)中大部分的差值為負,說明本文方法確定的底回波起點大部分處于相對靠后的位置。由第1節(jié)的分析可知,底回波的前段受到接收旁瓣的影響,雖然其與發(fā)射信號的相關(guān)性最大,但信噪比明顯下降,最終導(dǎo)致底速估計的精度不如底回波的后段。圖3(a)中底回波起點隨樣本延后的趨勢是因為試驗時測船隨風(fēng)浪搖晃。
分別使用本文底速估計方法、匹配濾波法估計測船靜止和運動條件下單波束、100個樣本點的底速,結(jié)果如圖4所示。靜止條件下的試驗結(jié)果對比如圖4(a)所示,兩種方法估計的速度均值都接近零,但匹配濾波方法的波動幅度明顯較大。運動條件下的試驗結(jié)果對比如圖4(b)所示。
圖4 兩種底速估計方法實測底速對比
圖4中兩種方法估計的速度均值約等于理論值,但本文方法的波動幅度明顯較小。測船運動條件下,底速曲線前后部分的速度均值變化明顯,是因為在第50個樣本點附近測船轉(zhuǎn)向,船速發(fā)生了改變。
兩種方法的底速估計標準差對比見表 1。相比匹配濾波法,測船靜止和運動條件下本文方法的速度標準差分別減小了60%和50%。
表1 兩種方法底速估計標準差對比 單位:m/s
利用第2.3節(jié)所示方法對測流試驗系統(tǒng)波束3的第2層徑向流速估計結(jié)果進行統(tǒng)計分析,總樣本數(shù)為100,徑向流速的理論值是1 m/s。計算得到的統(tǒng)計分析判據(jù)曲線如圖5(a)所示。判據(jù)大小介于0~1之間,根據(jù)經(jīng)驗確定判據(jù)小于 0.3的樣本不合格。不合格樣本用統(tǒng)計分析后合格樣本的均值替換。統(tǒng)計分析前后的徑向流速標準差對比如圖5(b)所示。圖5(b)中灰色點劃線表示統(tǒng)計分析前的徑向流速,其標準差為 0.779 m/s;圖 5(b)中紅色實線表示統(tǒng)計分析后的徑向流速,其標準差為 0.352 m/s。經(jīng)本文方法統(tǒng)計分析后,徑向流速標準差減小了約55%。
圖5 測流數(shù)據(jù)的統(tǒng)計分析
為了直觀展示統(tǒng)計分析判據(jù)與速度估計標準差之間的關(guān)系,將圖 5(a)中判據(jù)ηv按大小進行排序,并按排序后的樣本順序?qū)D5(b)中統(tǒng)計分析前的數(shù)據(jù)進行重排,再計算重排后每個樣本點前向15個樣本的速度估計標準差,第1~14個樣本的前向標準差均用第15個樣本的結(jié)果代替,處理結(jié)果如圖6所示。結(jié)果表明判據(jù)越大,回波的信噪比越高,對應(yīng)的流速估計誤差越小,與理論預(yù)測一致。
圖6 流速標準差與統(tǒng)計分析判據(jù)的關(guān)系
為提高寬帶多普勒測流方法的估計性能,本文基于隨機信號功率譜矩估計原理推導(dǎo)并設(shè)計了一種基于互相關(guān)系數(shù)的底速估計方法和一種基于信噪比的測流數(shù)據(jù)統(tǒng)計方法。兩種方法均可利用自相關(guān)運算實現(xiàn),具有計算復(fù)雜度低、實時性好、估計精度高等優(yōu)點。本文獲得的結(jié)論如下:
(1)由于寬帶測流中采用傾斜波束布局,底回波前段受到接收旁瓣的影響而信噪比明顯下降,底速估計時避開受旁瓣污染的回波段可顯著降低底速估計的誤差。
(2)|R(τ)|/R(0)表示延遲為τ的兩段回波分段的相關(guān)性,寬帶多普勒測流中,|R(τ)|/R(0)的值越大說明回波信噪比越高,流速估計的誤差也就越小。
(3)本文設(shè)計的基于互相關(guān)系數(shù)的底速估計方法可顯著提高底跟蹤速度的估計精度。實測數(shù)據(jù)分析結(jié)果表明,在測船靜止和測船運動兩種條件下,底跟蹤速度標準差分別減小了60%和50%。
(4)本文設(shè)計的基于信噪比的測流數(shù)據(jù)統(tǒng)計方法可明顯提高水流速度的估計精度。實測數(shù)據(jù)分析結(jié)果表明,在徑向流速理論值為1 m/s條件下,數(shù)據(jù)篩選后徑向流速標準差減小了約55%。