廖志宇, 朱廣平, 殷敬偉, 王成, 趙宿辰, 秦振林
(1.哈爾濱工程大學(xué) 水聲重點(diǎn)實(shí)驗(yàn)室, 黑龍江 哈爾濱 150001; 2.海洋信息獲取與安全工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室 (哈爾濱工程大學(xué)), 黑龍江 哈爾濱 150001; 3.哈爾濱工程大學(xué) 水聲工程學(xué)院, 黑龍江 哈爾濱 150001)
近年來由于地球變暖的趨勢,北極冰區(qū)面積持續(xù)減少,北極航道有望開通。又由于北極高緯度地區(qū)豐富的自然資源,引起了各海洋大國的高度關(guān)注。水聲技術(shù)是開發(fā)極地的有力手段,北極水聲學(xué)已成為新的研究熱點(diǎn)。北極地區(qū)地理位置獨(dú)特、氣候寒冷,北冰洋的大部分區(qū)域終年被海冰覆蓋,具有特殊的海洋環(huán)境,很多新問題和新現(xiàn)象漸漸地被從事該領(lǐng)域研究工作者所注意,研究發(fā)現(xiàn)了很多新的研究成果[1-4]。由于冰蓋的作用,北極也形成了獨(dú)特的聲場環(huán)境,造成冰下噪聲劇烈起伏以及強(qiáng)混響效應(yīng), 并形成了北冰洋獨(dú)有的半波導(dǎo)聲道。近些年來我國很多科研工作者在北極開展科學(xué)研究,劉洪寧等[5]研究發(fā)現(xiàn)加拿大海盆海冰邊緣區(qū)是聲體積后向散射強(qiáng)度的明顯過渡區(qū);劉崇磊等[6]開展了基于聲場模型的擴(kuò)頻通信技術(shù)在北極冰下聲信息傳輸中獲得了良好的通信性能的成果;黃海寧等[7]分析北極冰下聲信道多途結(jié)構(gòu),提出了基于OASES-Bellhop耦合模型的冰下聲信道多途結(jié)構(gòu)快速分析方法;朱廣平等[8]將Burke-Twersky(BT)散射模型與射線理論相結(jié)合研究了北極典型冰下的水聲信道特性,對認(rèn)知極地冰下水聲信道特性以及開展極地水聲系統(tǒng)性能預(yù)報(bào)具有一定意義。隨著水聲技術(shù)的不斷發(fā)展,冰下噪聲對聲吶檢測性能的影響也受到了越來越廣泛的關(guān)注[9]。在傳統(tǒng)的聲吶信號檢測性能分析中,常用高斯分布對其進(jìn)行建模和理論研究。然而,冰下噪聲信號脈沖值出現(xiàn)較頻繁,具有很強(qiáng)的脈沖特性,統(tǒng)計(jì)上表現(xiàn)出較厚的拖尾,具有很強(qiáng)的非高斯性[10]并且常包含很多的瞬變信號[11-12]特性。采用高斯模型對冰下噪聲進(jìn)行建模研究分析會導(dǎo)致檢測性能不準(zhǔn)確,在嚴(yán)重時(shí)甚至?xí)13]。因此有必要對冰下脈沖噪聲特性及其對檢測性能的影響進(jìn)行深入分析。
Alpha穩(wěn)定分布由Levy于1925年在研究廣義中心極限定理時(shí)提出,它是滿足廣義中心極限定理的唯一分布。田亞男等[14]分析了北極脈沖噪聲的特性并對其進(jìn)行了聲信道估計(jì),描述了北極冰噪聲具有厚重拖尾的特性。宋國麗等[15]表明了Alpha穩(wěn)定分布對尖峰態(tài)噪聲建模的可行性和適應(yīng)性。劉曄等[16]分析了淺海環(huán)境下混響形成過程及其非高斯特性,采用對稱 Alpha 穩(wěn)定分布理論對其進(jìn)行了建模,驗(yàn)證了模型的有效性。所以,采用Alpha穩(wěn)定分布對有關(guān)冰下噪聲方面的建模具有可行性。
本文分析了第8次北極科考聲學(xué)實(shí)驗(yàn)采集的冰下噪聲及松花江冰下噪聲統(tǒng)計(jì)特性,采用Alpha穩(wěn)定分布對噪聲進(jìn)行建模,采用蒙特卡羅方法在實(shí)際噪聲及仿真噪聲下分析匹配濾波器的檢測性能,重點(diǎn)分析了不同參數(shù)的Alpha穩(wěn)定分布噪聲、不同的發(fā)射信號及信噪比對檢測性能的影響。
Alpha穩(wěn)定分布除了幾個特例外,沒有統(tǒng)一、封閉的解析表達(dá)式,通常使用特征函數(shù)來描述其分布特性,其特征函數(shù)表達(dá)式形式為[17-18]:
φ(x)=exp{iδx-|γx|α[1+iβsgn(x)ω(x,α)]}
(1)
式中:
(2)
(3)
式中:α、β、γ、δ確定了Alpha穩(wěn)定分布的特征函數(shù)。
假設(shè)變量X服從Alpha穩(wěn)定分布,則有X~Sα(γ,β,δ),各個參數(shù)的性質(zhì)如下:
1)α是Alpha穩(wěn)定分布的特征指數(shù),決定了噪聲脈沖特性的強(qiáng)弱即概率密度函數(shù)的拖尾厚度。其數(shù)值范圍是0<α≤2,α趨近于0的時(shí)候,噪聲的脈沖特性越來越明顯,概率密度函數(shù)的拖尾程度加重。當(dāng)α=2時(shí),Alpha穩(wěn)定分布退化為高斯分布即S2(γ,0,δ)=N(δ,2γ2);當(dāng)α=1時(shí),Alpha穩(wěn)定分布退化為柯西分布S1(γ,0,δ);當(dāng)α=0.5,β=1的時(shí)候Alpha穩(wěn)定分布退化為Levy分布S0.5(γ,1,δ)。
2)β是Alpha穩(wěn)定分布的偏斜參數(shù),表示分布的傾斜程度。其數(shù)值范圍是:-1≤β≤1,當(dāng)β<0時(shí),Alpha穩(wěn)定分布的概率密度函數(shù)的峰向右傾斜,分布向左傾斜;當(dāng)β>0時(shí),概率密度函數(shù)的峰向左傾斜,分布向右傾斜。當(dāng)β=0時(shí),概率密度函數(shù)為對稱的,此時(shí)分布稱為對稱Alpha穩(wěn)定分布。
3)γ是Alpha穩(wěn)定分布的尺度參數(shù),描述了Alpha穩(wěn)定分布的離散程度,其數(shù)值范圍是0≤γ<∞,類似于高斯分布的方差。
4)δ是Alpha穩(wěn)定分布的位置參數(shù),決定了概率密度函數(shù)峰值所處位置,其范圍是-∞<δ<+∞。當(dāng)β=0,γ=1,δ=0,Alpha穩(wěn)定分布對應(yīng)標(biāo)準(zhǔn)SαS分布。
可以通過S1參數(shù)系、S2參數(shù)系、標(biāo)準(zhǔn)參數(shù)系生成服從Alpha穩(wěn)定分布的隨機(jī)變量。S1參數(shù)系消除了標(biāo)準(zhǔn)參數(shù)系特征函數(shù)存在的不連續(xù)性[19],S2參數(shù)系下隨機(jī)變量產(chǎn)生的方法便于理論證明[20],采用S2參數(shù)系以及Alpha穩(wěn)定分布的基本性質(zhì)生成服從Alpha穩(wěn)定分布的隨機(jī)變量。
產(chǎn)生2個互相獨(dú)立的隨機(jī)變量V和W,其中V是服從(-π/2,π/2)的均勻分布,W是服從均值為1的指數(shù)分布。
當(dāng)α≠1時(shí),
(4)
(5)
(6)
當(dāng)α=1時(shí),β2=β:
(7)
(8)
通過式(8)得到X即是服從Alpha穩(wěn)定分布的隨機(jī)變量。
滿足位置參數(shù)δ=0的實(shí)對稱Alpha穩(wěn)定分布隨機(jī)變量具有有限的分?jǐn)?shù)低階矩,本文采用分?jǐn)?shù)低階矩法對Alpha穩(wěn)定分布的α、γ進(jìn)行參數(shù)估計(jì)[21],測試過程中發(fā)現(xiàn)α、γ過大時(shí)發(fā)生對參數(shù)γ的估計(jì)嚴(yán)重不準(zhǔn)確,故對公式加以稍加改動,在不影響對α估計(jì)性能的前提下改善了對參數(shù)γ的估計(jì)性能,具體過程如下:
令隨機(jī)變量的正負(fù)階數(shù)為p、q,則正負(fù)階矩分別為:
(9)
(10)
(11)
令p=-q,以保證正負(fù)階矩均有限,則有0
0
(12)
進(jìn)而根據(jù)式(13) 估計(jì)出參數(shù)γ的值:
(13)
p值越小,對于α、γ的估計(jì)越準(zhǔn)確[22]。在估計(jì)的過程中,樣本點(diǎn)數(shù)為5 000,p=0.1,當(dāng)估計(jì)α的數(shù)值時(shí)令γ=1,估計(jì)γ的數(shù)值時(shí)令α=1對α、γ的估計(jì)情況如表1。
表1 參數(shù)α、γ估計(jì)情況統(tǒng)計(jì)表Table 1 Parameter α,γ estimation statistics table
分析冰下噪聲的高階統(tǒng)計(jì)量如偏度值和峰度值,也有利于充分認(rèn)知其非高斯特性,峰度和偏度均無量綱。偏度(Skewness)反映了隨機(jī)變量分布對稱情況,定義為樣本的三階矩,其計(jì)算公式為:
(14)
峰度(Kurtosis) 反映了隨機(jī)變量分布的尖銳程度,峰度越大,表現(xiàn)在分布上面是中心點(diǎn)越尖銳。在相同方差的情況下,峰度越大,中間一大部分的值方差都很小,必須有一些值離中心點(diǎn)越遠(yuǎn),所以這就是所說的“厚尾”,定義為四階標(biāo)準(zhǔn)矩,其計(jì)算公式為:
(15)
式中:μ是均值;δ為標(biāo)準(zhǔn)差;E是均值操作。
當(dāng)輸入為確定信號加上平穩(wěn)噪聲時(shí),能夠使輸出信噪比達(dá)到最大值的線性系統(tǒng)稱為匹配濾波器。
(16)
式中:s0(T)為輸出信號;n0(T)為輸出噪聲;h0(t)為沖激響應(yīng)。使SNR最大等效于在輸出信號功率E[s02(T)]為常數(shù)的約束條件下,使輸出噪聲功率E[n02(T)]最小。這是一個有約束條件的最優(yōu)化問題,應(yīng)用拉格朗日(Lagrange)乘數(shù)法構(gòu)建目標(biāo)函數(shù)對其求導(dǎo),當(dāng)施瓦茲不等式中等號成立時(shí):
(17)
式(17)就是匹配濾波器傳輸函數(shù)的普遍形式。匹配濾波器的作用主要有2個方面:1)使濾波器輸出有用信號成分盡可能強(qiáng);2)抑制噪聲,使濾波器輸出噪聲成分盡可能小,減小噪聲對信號處理的影響。在聲吶雷達(dá)信號處理中這2個作用都是以輸出信噪比最大為準(zhǔn)則來提高信噪比。
首先將噪聲樣本(北極噪聲樣本或仿真噪聲樣本)與拷貝信號(發(fā)射信號)進(jìn)行疊加形成目標(biāo)回波樣本集(以下統(tǒng)稱樣本集),然后將拷貝信號分別與噪聲樣本以及目標(biāo)回波樣本集做匹配濾波進(jìn)行信號檢測。采用蒙特卡羅方法做1 000次重復(fù)匹配濾波實(shí)驗(yàn)。根據(jù)信號樣本和噪聲樣本的檢測結(jié)果,設(shè)置一定的檢測閾值。在這一過程中,每次計(jì)算出2個重要量的值,ROC曲線的縱軸是“檢測概率”(true positive rate, TPR),橫軸是“虛警概率”(false positive rate,FPR)分別以它們?yōu)闄M、縱坐標(biāo)作圖,就得到了“ROC曲線”,該曲線描述了在不同虛警情況下對信號的檢測情況。檢測性能分析流程如圖1所示。
圖1 檢測性能分析流程Fig.1 Flow chart of test performance analysis
研究過程中所采用的數(shù)據(jù)為松花江實(shí)驗(yàn)采集的冰下噪聲數(shù)據(jù)和北極冰下噪聲數(shù)據(jù),松花江噪聲數(shù)據(jù)為2020年1月份在松花江冰面破洞,由水聽器采集獲得,溫度為-18 ℃左右,冰層厚度在46 cm左右,北極噪聲數(shù)據(jù)為第8次北極科考采集的數(shù)據(jù)。對噪聲數(shù)據(jù)進(jìn)行脈沖特性分析,展示部分冰下噪聲樣本圖,繪制了其相應(yīng)的統(tǒng)計(jì)分布直方圖并采用Alpha穩(wěn)定分布和高斯分布對噪聲數(shù)據(jù)進(jìn)行了擬合并對比了二者擬合效果,噪聲數(shù)據(jù)的采樣頻率fs=50 kHz。
圖2描述了松花江冰下噪聲數(shù)據(jù)及擬合情況。圖2(a)、(b)是在松花江冰下噪聲數(shù)據(jù)時(shí)域圖,圖2(c)、(d) 描述了樣本方差隨松花江噪聲樣本數(shù)的增加的變化情況,圖2(e)、(f)繪制了2段數(shù)據(jù)其相應(yīng)的統(tǒng)計(jì)分布直方圖并分別用Alpha分布和高斯分布對冰下噪聲數(shù)據(jù)進(jìn)行了擬合。圖2(e)中高斯分布的均值μ=0、方差δ=0.211 1,采用分?jǐn)?shù)階矩法估計(jì)噪聲數(shù)據(jù)獲得的Alpha穩(wěn)定分布的參數(shù)α=1.661 7、γ=0.245 8。圖2(f)中高斯分布的均值μ=0、方差δ=0.274 4,采用分?jǐn)?shù)階矩法估計(jì)噪聲數(shù)據(jù)獲得的Alpha穩(wěn)定分布的參數(shù)α=1.600 5、γ=0.268 5。圖2(g)、(h)分別對應(yīng)松花江兩段噪聲的頻譜情況,可知能量主要集中在低頻部分。
圖3描述了北極冰下噪聲數(shù)據(jù)及擬合情況。圖3中,(a)、(b)是在北極冰下噪聲數(shù)據(jù)時(shí)域圖,圖3(c)、(d) 描述了樣本方差隨北極噪聲樣本數(shù)的增加的變化情況,圖3(e)、(f)繪制了2段數(shù)據(jù)其相應(yīng)的統(tǒng)計(jì)分布直方圖并分別用Alpha分布和高斯分布對冰下噪聲數(shù)據(jù)進(jìn)行了擬合。圖3(e)中高斯分布的均值μ=0、方差δ=0.126 8,采用分?jǐn)?shù)階矩法估計(jì)噪聲數(shù)據(jù)獲得的Alpha穩(wěn)定分布的參數(shù)α=1.721 7、γ=0.188 2。圖3(f)中高斯分布的均值μ=0、方差δ=0.093 3,采用分?jǐn)?shù)階矩法估計(jì)噪聲數(shù)據(jù)獲得的Alpha穩(wěn)定分布的參數(shù)α=1.825 1、γ=0.186 0,圖3(g)、(h)分別對應(yīng)北極2段噪聲的頻譜情況,能量主要集中在低頻部分,相比于松花江噪聲,北極冰下噪聲在低頻部分的能量更加集中。
根據(jù)圖2、3(c)、(d)及估計(jì)出的分布參數(shù)可以看出:噪聲數(shù)據(jù)的樣本方差隨噪聲樣本數(shù)的增加不收斂于常數(shù),根據(jù)文獻(xiàn)[15]可知松花江冰下噪聲和北極冰下噪聲均為非高斯分布;表2描述了噪聲數(shù)據(jù)對應(yīng)的峰度值和偏度值,從冰下噪聲的峰度值亦可知冰下噪聲為非高斯分布且具有沖擊特性,從偏度值可知冰下噪聲數(shù)據(jù)分布大致符合對稱分布;對冰下噪聲數(shù)據(jù)的擬合時(shí),Alpha穩(wěn)定分布擬合效果明顯優(yōu)于高斯分布,Alpha穩(wěn)定分布更加貼近于冰下噪聲實(shí)際情況。
表2 冰下噪聲的偏度和峰度Table 2 Skewness and kurtosis of under ice noise
圖4、5描述松花江和北極2種冰下噪聲在不同頻段條件下的統(tǒng)計(jì)和擬合情況,從圖4、5(a)、(b)中看出,1~10 kHz的噪聲能量比較大,對比圖4、5(c)、(d),可知隨著頻率的提高,噪聲的能量逐漸減小,即噪聲能量主要集中在低頻段范圍內(nèi)。
圖2 松花江冰下噪聲數(shù)據(jù)統(tǒng)計(jì)及擬合Fig.2 Statistics and fitting of under ice noise data of Songhua River
根據(jù)冰下噪聲的脈沖特性以及Alpha穩(wěn)定分布的擬合情況,根據(jù)2.3節(jié)可基于不同參數(shù)的Alpha穩(wěn)定分布仿真出的冰下噪聲用來研究對匹配濾波器檢測性能的影響因素。圖6是令參數(shù)β=0,δ=0,改變參數(shù)α、γ仿真得到的服從Alpha穩(wěn)定分布的冰下噪聲。
采用第8次北極科考采集的冰下噪聲數(shù)據(jù)和松花江冰下噪聲數(shù)據(jù)(采樣頻率fs=50 kHz)為噪聲背景,選取噪聲數(shù)據(jù)時(shí)寬為20 ms,采用CW和線性調(diào)頻脈沖信號與實(shí)際噪聲疊加后作為接收信號進(jìn)行匹配濾波檢測,繪制相應(yīng)的ROC曲線,分析檢測性能。檢測性能分析過程如2.2節(jié)所示。實(shí)驗(yàn)中設(shè)置CW信號的頻率為1 kHz,脈寬為20 ms,LFM信號的起始頻率為1 kHz,帶寬為4 kHz,脈寬為20 ms。
圖3 北極冰下噪聲數(shù)據(jù)統(tǒng)計(jì)及擬合Fig.3 Statistics and fitting of under ice noise data of Arctic
圖7(a)、(c)描述了不同SNR下CW信號在松花江噪聲背景和北極噪聲背景下的檢測情況,圖7(b)、(d)描述了不同SNR下LFM信號在松花江噪聲背景和北極噪聲背景下的檢測情況。從4幅圖中都可以看出隨著SNR的增加,對信號的檢測性能增強(qiáng)。表3中,虛警概率固定為5%,在松花江和北極冰下噪聲下的LFM信號和CW信號的檢測概率可以看出,2種背景噪聲下對LFM信號的檢測能力比對CW信號的檢測情況要好得多,CW信號的檢測對SNR的要求比LFM信號的檢測高,在低信噪比的情況下, CW的檢測能力很差, LFM的檢測情況較好。
圖4 不同頻段下松花江冰下噪聲1的統(tǒng)計(jì)及擬合Fig.4 Statistics and fitting of under ice noise data 1 of Songhua River at different frequency bands
圖5 不同頻段下北極冰下噪聲數(shù)據(jù)1的統(tǒng)計(jì)及擬合Fig.5 Statistics and fitting of under ice noise data 1 of Arctic at different frequency bands
圖6 服從Alpha穩(wěn)定分布的冰下噪聲Fig.6 Under ice noise with Alpha stable distribution
圖7 不同情形下匹配濾波檢測性能Fig.7 Matched filter detection performance in different situations
通過仿真生成采樣頻率fs=50 kHz,符合Alpha穩(wěn)定分布的隨機(jī)變量作為冰下噪聲數(shù)據(jù),根據(jù)Alpha穩(wěn)定分布統(tǒng)計(jì)參數(shù)生成不同噪聲背景,研究參數(shù)對檢測性能的影響??刂菩旁氡萐NR=-16 dB,改變Alpha穩(wěn)定分布的參數(shù)α生成對應(yīng)的冰下噪聲數(shù)據(jù),與LFM信號疊加后,采用匹配濾波進(jìn)行檢測并作ROC曲線圖分析其檢測性能。LFM信號的起始頻率為1 kHz,帶寬為4 kHz,脈寬為20 ms。設(shè)定Alpha穩(wěn)定分布的參數(shù)β=0,γ=1,δ=0,α的數(shù)值依次為0.4、0.6、1.0、1.4、1.6、2.0,生成冰下噪聲數(shù)據(jù),得到不同特征指數(shù)α對應(yīng)的ROC曲線;設(shè)定Alpha穩(wěn)定分布的參數(shù)α=1.0,β=0,δ=0,γ的數(shù)值依次為0.2、0.4、0.6、0.8、1.0,生成冰下噪聲數(shù)據(jù),得到不同尺度參數(shù)γ對應(yīng)的ROC曲線;根據(jù)1.2節(jié),當(dāng)α的數(shù)值固定為2.0時(shí),此時(shí)的Alpha穩(wěn)定分布即為高斯分布。以下對應(yīng)不同噪聲背景下對LFM信號檢測的ROC曲線。
圖8(a)描述了不同特征指數(shù)α對應(yīng)的LFM信號檢測情況,圖8(b)描述了虛警概率為5%時(shí)各參數(shù)所對應(yīng)的檢測概率,可以看出隨著α的逐漸增加,檢測概率隨之降低,檢測能力下降,也說明了采用高斯分布生成的變量作為冰下噪聲數(shù)據(jù)會對檢測性能產(chǎn)生一定的影響;
圖8(c)描述了不同尺度參數(shù)γ對應(yīng)的LFM信號檢測情況,圖8(d)描述了虛警概率為5%時(shí)各參數(shù)所對應(yīng)的檢測概率,可以得到隨著γ的逐漸增加,檢測概率浮動較小,對檢測性能的影響較小。
表3虛警概率5%時(shí)檢測情況統(tǒng)計(jì)表
Table3Statisticaltableofdetectionwhenfalsealarmprobabilityis5%
接收信號信噪比/dB檢測概率/%接收信號信噪比/dB檢測概率/%松花江噪聲+CW信號北極噪聲+CW信號-899.1-1097.1-1289.2-1462.5-1628.9-898.8-1093.3-1273.0-1441.1-1616.1松花江噪聲+LFM信號北極噪聲+LFM信號-12100.0-1499.4-1691.8-1849.3-2013.1-18100-20100-2297.1-2477.5-2553.5
1) 北極噪聲和松花江噪聲數(shù)據(jù)的統(tǒng)計(jì)特性都具有非高斯特性,且偏度在0附近,基本上滿足對稱分布;本次收集的松花江噪聲數(shù)據(jù)的峰度值略大于北極噪聲數(shù)據(jù)的峰度值,北極噪聲數(shù)據(jù)估計(jì)得的α數(shù)值大于松花江噪聲數(shù)據(jù)估計(jì)的數(shù)值。
2) 在統(tǒng)計(jì)特性方面,采用Alpha穩(wěn)定分布建模的北極噪聲和松花江噪聲數(shù)據(jù)相比于高斯噪聲,更加貼近于冰下噪聲實(shí)際情況。
3) 在本次收集的北極和松花江冰下噪聲數(shù)據(jù)下對不同的信號形式及參數(shù)的檢測性能分析發(fā)現(xiàn),匹配濾波器對LFM信號的檢測能力優(yōu)于CW信號;信噪比相同的情形下,隨著Alpha穩(wěn)定分布參數(shù)α的逐漸增加,對信號的檢測能力逐漸減弱。