利用海洋環(huán)境噪聲估計(jì)海流流速?

汪 愷 李風(fēng)華 楊習(xí)山

(1 中國(guó)科學(xué)院聲學(xué)研究所 北京 100190)

(2 中國(guó)科學(xué)院大學(xué) 北京 100049)

0 引言

海洋聲層析方法自提出至今得到了長(zhǎng)久的發(fā)展，主要被用來(lái)估計(jì)海洋聲速結(jié)構(gòu)和海流[1?2]。過(guò)去數(shù)十年，諸多遠(yuǎn)場(chǎng)聲傳播實(shí)驗(yàn)已經(jīng)驗(yàn)證論述了海洋聲層析方法的有效性[2]。由于深海聲道軸的存在，大部分海流測(cè)量實(shí)驗(yàn)都在深海進(jìn)行[3]。而淺海環(huán)境下，由于聲速?gòu)暮１砻娴胶５鬃兓容^小、聲線分辨性較差、噪聲干擾強(qiáng)烈，因此淺海聲層析在海流估計(jì)上的應(yīng)用較少。淺海海域的海流結(jié)構(gòu)以水平的潮汐余流為主。在聲傳播過(guò)程中，海流的流速相對(duì)于聲速來(lái)說(shuō)是一個(gè)極小量，對(duì)于聲傳播時(shí)間并不敏感，而且海流方向多變，因此難以將海流流速與聲速剖面一起反演。從理論計(jì)算上，聲層析方法可以應(yīng)用在海流估計(jì)中[4]。Kaneko 等[5]在日本瀨戶內(nèi)海進(jìn)行聲層析實(shí)驗(yàn)并進(jìn)行了10 km的二維水平流速場(chǎng)重建，Zhu 等[6]完成了對(duì)中國(guó)舟山群島海區(qū)的二維流場(chǎng)重建。Taniguchi 等采用聲層析方法成功測(cè)量了Kuroshio海域的流速垂直分布結(jié)構(gòu)[7]，并用同樣的方法實(shí)現(xiàn)了對(duì)呂宋海峽的流速估計(jì)[8]。但上述海流測(cè)量基于的聲層析實(shí)驗(yàn)或是需要主動(dòng)發(fā)射聲源，或者需要長(zhǎng)時(shí)間接收信號(hào)才能對(duì)流速進(jìn)行估計(jì)，難以觀測(cè)海流在短時(shí)間尺度上的變化。

通過(guò)長(zhǎng)時(shí)間累積平均的環(huán)境噪聲互相關(guān)函數(shù)(Ambient noise cross-correlation function)，其時(shí)間微分形式與兩個(gè)接收器間經(jīng)驗(yàn)格林函數(shù)(Empirical Green’s function)的時(shí)間到達(dá)結(jié)構(gòu)已經(jīng)被論證過(guò)是一致的[9?10]。在此基礎(chǔ)上，Godin 等從公路環(huán)境噪聲成功反演了風(fēng)速[11]，適合進(jìn)行長(zhǎng)期聲速觀測(cè)的被動(dòng)海洋聲層析也得到了迅速發(fā)展[12?13]。并且，被動(dòng)聲層析方法也更適合人為因素干擾較多的淺海海域。Li等通過(guò)波束形成加速能量累積，在2 h時(shí)間內(nèi)提取出兩個(gè)平行水平陣列間的經(jīng)驗(yàn)格林函數(shù)，準(zhǔn)確地反演了2 h 平均的聲速剖面[14]。

本文利用環(huán)境噪聲互相關(guān)提取格林函數(shù)在2 h時(shí)間尺度上估計(jì)海流流速。首先介紹波束形成增強(qiáng)水平陣列間環(huán)境噪聲互相關(guān)函數(shù)，計(jì)算得到經(jīng)驗(yàn)格林函數(shù)的理論方法。然后依據(jù)南中國(guó)海海域聲層析實(shí)驗(yàn)數(shù)據(jù)，獲取經(jīng)驗(yàn)格林函數(shù)的時(shí)間到達(dá)結(jié)構(gòu)，估計(jì)反演2 h 平均的海流流速。再給出該海域海流流速仿真計(jì)算，驗(yàn)證實(shí)驗(yàn)數(shù)據(jù)處理所獲海流流速和理論仿真估計(jì)結(jié)果的一致性、有效性，并分析海深失配、距離失配對(duì)流速估計(jì)的影響。在此基礎(chǔ)上，歸納總結(jié)主要結(jié)論：基于聲層析原理，利用海洋環(huán)境噪聲可以有效反演估計(jì)海流流速。

1 理論

兩點(diǎn)間經(jīng)驗(yàn)格林函數(shù)的提取主要是通過(guò)累積端射方向上的噪聲源能量，因此兩個(gè)陣列之間的信號(hào)可以在端射方向上做波束形成來(lái)加快能量累積[14]。如圖1(a)的設(shè)計(jì)方案所示，當(dāng)滿足遠(yuǎn)場(chǎng)假設(shè)時(shí)，陣列A 和陣列B 在θ0方向上波束的頻域噪聲互相關(guān)函數(shù)可以用式(1)表達(dá)：

圖1 被動(dòng)聲層析設(shè)計(jì)示意圖與海洋分層示例Fig.1 Schematic design of passive acoustic tomography and an example of ocean stratification

這里，(ω,θ0)是波束形成的權(quán)重向量，PA(B)(ω)是陣列的頻域接收信號(hào)，(·)H是共軛轉(zhuǎn)置，ω是角頻率，θ0是陣列A 的第m個(gè)陣元指向陣列B的第n個(gè)陣元的端射方向。

那么，θ0方向上的格林函數(shù)就可以從時(shí)域噪聲互相關(guān)函數(shù)的時(shí)域微分提取出來(lái)：

由于波束形成會(huì)過(guò)濾掉其他方向上的能量，為了得到反向的格林函數(shù)，通過(guò)將波束形成指向θ0的反方向，陣列B 到陣列A 的格林函數(shù)也可以同樣的形式被提取出來(lái)：

其中，GBA(t)相當(dāng)于未做波束形成時(shí)GAB(t)的時(shí)間負(fù)半軸，它們間的關(guān)系是GAB(t)=?GBA(?t)。

聲線到達(dá)時(shí)間可以依據(jù)正反兩個(gè)方向的格林函數(shù)得到。在對(duì)海洋進(jìn)行分層后，可以將聲線軌跡劃分為N層，圖1(b)給出了將海洋劃分為3 層的例子。此時(shí)正反兩個(gè)方向上第i條本征聲線的到達(dá)時(shí)間可以通過(guò)下式計(jì)算[6]:

其中，cj和uj分別是第j層上距離平均聲速和海流流速，ri,j代表第i條本征聲線在第j段路徑上的傳播路程，terr是兩個(gè)接收器之間的時(shí)鐘誤差。由于cj ?uj，忽略小量uj后，兩個(gè)陣列間的聲線傳播時(shí)間差可以表示為

其中，terr可以通過(guò)從所有聲線中減去特定的聲線傳播時(shí)間差來(lái)消除，通過(guò)求解式(5)，即可得到相關(guān)海域的各層海流流速。如果對(duì)海洋不進(jìn)行分層的話，這里求解式(5)將退化為

2 實(shí)驗(yàn)

2.1 實(shí)驗(yàn)描述

在海南島東南近海進(jìn)行了一次聲層析實(shí)驗(yàn)，實(shí)驗(yàn)環(huán)境如圖1(a)所示，在海深97 m的平坦海底布放兩個(gè)連接的水平線陣，兩者大致平行，相距3529 m。并且每個(gè)陣列有15個(gè)等間距間隔6 m的陣元，實(shí)際陣列長(zhǎng)度大約為84 m。水聽(tīng)器以5000 Hz的采樣率同步連續(xù)不斷地記錄寬帶的海洋環(huán)境噪聲信號(hào)。同時(shí)在實(shí)驗(yàn)期間，一條包含38個(gè)溫深儀的溫度鏈被豎直布放在兩個(gè)水平陣列的中央，用以記錄水體大約2.5～95.0 m 深度范圍的水文數(shù)據(jù)。分別使用了4 段海洋噪聲接收信號(hào)，具體測(cè)量時(shí)間如表1 所示，對(duì)應(yīng)的部分水文數(shù)據(jù)如圖2所示。

表1 噪聲信號(hào)接收時(shí)間Table 1 Noise signal receiving time

圖2 溫度鏈實(shí)際測(cè)量的聲速剖面Fig.2 Sound speed profiles measured by temperature chain

2.2 數(shù)據(jù)處理

將每個(gè)陣元接收的海洋噪聲信號(hào)按每2 h 分為一個(gè)時(shí)間段后，按10 s 一拍對(duì)各段時(shí)域信號(hào)進(jìn)一步劃分。在快速傅里葉變換后，對(duì)每個(gè)快拍的頻域信號(hào)做40～350 Hz的帶寬濾波。對(duì)濾波后的頻域信號(hào)做波束形成來(lái)加快能量累積，提升信噪比，波束形成的參考聲速設(shè)置為1530 m/s。對(duì)各快拍的波束輸出累積求和后，可以在2 h 時(shí)間上平均得到累積波束輸出。因此，2 h時(shí)間累積的噪聲互相關(guān)函數(shù)可以通過(guò)式(1)計(jì)算得到。接著通過(guò)式(2)～(3)端射方向上的經(jīng)驗(yàn)格林函數(shù)即可被提取出來(lái)。

由于兩個(gè)水平線陣在物理上相互連接，因此兩個(gè)陣列的水聽(tīng)器記錄的噪聲信號(hào)都是同步且連續(xù)的，因此式(5)～(6)中terr可以忽略不計(jì)。由于聲線傳播的多途效應(yīng)，沒(méi)有經(jīng)過(guò)海面反射的本征聲線以及前兩次海面反射波都混在了第一個(gè)信號(hào)包絡(luò)中[14]。第一個(gè)信號(hào)包絡(luò)中各組本征聲線難以辨別區(qū)分，而后續(xù)包絡(luò)的本征聲線都穿過(guò)了水體所有深度，因此難以對(duì)海洋進(jìn)行分層。第二個(gè)信號(hào)包絡(luò)對(duì)應(yīng)了海面三次反射波，物理意義清晰。所以，通過(guò)式(6)，在利用幾何關(guān)系計(jì)算聲線路徑長(zhǎng)度后，不對(duì)海洋進(jìn)行分層，第二個(gè)信號(hào)包絡(luò)的時(shí)延差被用來(lái)計(jì)算估計(jì)海域深度上平均的海流流速。

2.3 實(shí)驗(yàn)結(jié)果

定義從陣列A 到陣列B 為正向，提取出端射方向的經(jīng)驗(yàn)格林函數(shù)。圖3 展示了自2 h 累積噪聲互相關(guān)函數(shù)提取出的經(jīng)驗(yàn)格林函數(shù)的正半軸部分，共有48 h。可以從圖中看到，在2 h 累積時(shí)間下，到達(dá)波的信號(hào)包絡(luò)都比較清晰穩(wěn)定。

圖3 從2 h 噪聲互相關(guān)函數(shù)中提取的正向格林函數(shù)，共持續(xù)48 hFig.3 Positive Green’s functions extracted from the 2-hour noise cross-correlation functions for a total of 48 hours

同時(shí)，圖4 展示了端射方向上提取出的一組經(jīng)驗(yàn)格林函數(shù)，其中圖4(a)將兩個(gè)經(jīng)驗(yàn)格林函數(shù)以正負(fù)時(shí)間軸的形式統(tǒng)一歸一化，圖4(b)為正反向分別歸一化后的格林函數(shù)。可以從圖4(a)看出，正反兩個(gè)方向上信號(hào)在統(tǒng)一歸一化后反方向的格林函數(shù)信噪比更高。

圖4 自2 h 噪聲互相關(guān)函數(shù)提取出的一組格林函數(shù)Fig.4 A Green’s function extracted from the 2-hour noise cross-correlation function

圖5 展示了兩個(gè)水平陣列之間端射方向上海流流速的反演結(jié)果。通過(guò)計(jì)算如圖4(b)所示的第二個(gè)信號(hào)包絡(luò)的時(shí)延差，在2 h 時(shí)間尺度內(nèi)計(jì)算估計(jì)出了海流流速。

圖5 從2 h 經(jīng)驗(yàn)格林函數(shù)計(jì)算得到的海流流速Fig.5 Current velocities calculated from 2-hour empirical Green’s functions

2.4 實(shí)驗(yàn)分析

圖3 展示了連續(xù)48 h 的經(jīng)驗(yàn)格林函數(shù)，可以發(fā)現(xiàn)在2 h 的時(shí)間累積下，大部分時(shí)間到達(dá)結(jié)構(gòu)都比較清晰，但也存在個(gè)別到達(dá)波的包絡(luò)不清晰。這是因?yàn)樵谠摾鄯e時(shí)間段內(nèi)，端射方向上噪聲源較弱，沒(méi)有足夠的能量積累，可以通過(guò)延長(zhǎng)累積時(shí)間來(lái)提高格林函數(shù)的信噪比。

從圖4(a)中可以看到在正反兩個(gè)方向上所提取出的經(jīng)驗(yàn)格林函數(shù)信噪比并不一致，這是由于在實(shí)際海洋環(huán)境中噪聲源分布是不均勻的。同時(shí)從圖4(b)中可以看到盡管信噪比存在差異，但正反兩個(gè)方向上得到的經(jīng)驗(yàn)格林函數(shù)時(shí)間到達(dá)結(jié)構(gòu)都清晰明確，已經(jīng)足夠提取得到到達(dá)時(shí)延來(lái)進(jìn)行后續(xù)的流速計(jì)算。

圖5 展示了從第二個(gè)信號(hào)包絡(luò)時(shí)延差計(jì)算估計(jì)的海流流速，其數(shù)值范圍在50 cm/s以內(nèi)。附近海域歷史測(cè)量平均流速為52 cm/s，最大流速109 cm/s，最小流速為18 cm/s[15]，考慮到海流的方向與實(shí)驗(yàn)波束形成指向的不同，反演結(jié)果落在歷史測(cè)量的合理范圍內(nèi)。從圖5 中可以看到海流流速在2 h 時(shí)間尺度上隨時(shí)間的變化趨勢(shì)，同時(shí)隨著季節(jié)改變，海流方向也發(fā)生了改變。已有的研究結(jié)果表明，海南島東南海域受南中國(guó)海大尺度表層環(huán)流的影響，海流在夏季和冬季分別以東北向和西南向流為主[15]，海流方向發(fā)生改變符合歷史實(shí)際觀測(cè)結(jié)果。

3 射線理論仿真

3.1 仿真與結(jié)果

通過(guò)射線追蹤方法可以對(duì)兩點(diǎn)之間的聲傳播過(guò)程進(jìn)行仿真[16?17]。仿真使用的是由溫度鏈實(shí)際測(cè)量的一組聲速剖面，海流流速則采用恒定流速剖面，分別設(shè)置為1 m/s、0.5 m/s 和0.2 m/s。兩點(diǎn)之間的距離設(shè)置為3500 m，海底聲速由頭波反演得到設(shè)置為1603 m/s，海底密度為1.9 g/cm3，吸收系數(shù)為0.5 dB/λ，其余參數(shù)與實(shí)驗(yàn)參數(shù)一致。

圖6 展示了仿真環(huán)境下的本征聲線軌跡?？梢詮膱D6 中發(fā)現(xiàn)，雖然本征聲線可以直觀分為直達(dá)波、海底一次反射以及海面反射3 組聲線，通過(guò)3 組聲線穿過(guò)水體的深度不同來(lái)對(duì)水體進(jìn)行分層來(lái)估計(jì)流速，但在與仿真對(duì)應(yīng)的實(shí)際實(shí)驗(yàn)環(huán)境中，直達(dá)波、海底一次反射波以及海面二次反射都對(duì)應(yīng)著圖3 和圖4(b)中第一個(gè)到達(dá)的信號(hào)包絡(luò)，難以區(qū)分，因此在實(shí)際的淺海環(huán)境中很難通過(guò)將射線分為3 組來(lái)對(duì)海流流速進(jìn)行分層估計(jì)。

圖6 基于射線理論仿真得到的本征聲線軌跡Fig.6 The simulated eigen rays trajectory based on the ray theory

表2 列出了3 組聲線的到達(dá)時(shí)間和利用3 次海面反射波的到達(dá)時(shí)間估計(jì)出的海流流速。從表2 中可以看出反演的海流流速與仿真的實(shí)際環(huán)境十分接近。

表2 仿真反演的結(jié)果Table 2 Results of simulation inversions

3.2 誤差分析

格林函數(shù)的到達(dá)時(shí)間主要受聲線在水體中傳播的距離影響，對(duì)應(yīng)的主要不確定性在于海水深度。在聲線傳播過(guò)程中，如果海底深度發(fā)生變化，聲線軌跡將會(huì)發(fā)生變化，極大地影響聲線傳播時(shí)間。而實(shí)際實(shí)驗(yàn)時(shí)，海底深度失配也是通常無(wú)法回避的問(wèn)題。

表3 比較了不同海底深度條件下對(duì)3 次海面反射波傳播時(shí)間和反演結(jié)果的影響，?z為正表示計(jì)算時(shí)海底深度大于實(shí)際海底深度，仿真所用參數(shù)與3.1 節(jié)參數(shù)一致，仿真流速設(shè)定為恒定1 m/s?？梢詮谋? 中發(fā)現(xiàn)本實(shí)驗(yàn)環(huán)境下深度失配對(duì)正反向聲傳播時(shí)間影響較小。雖然通過(guò)計(jì)算時(shí)延差，可以消除部分誤差，但如果時(shí)延差未能完全消除誤差時(shí)，0.1 ms 的時(shí)延誤差將帶來(lái)約0.033 m/s 的流速估計(jì)誤差。仿真表明在本實(shí)驗(yàn)環(huán)境下，5 m左右的深度失配可能會(huì)輕微影響流速反演的結(jié)果，因此實(shí)際實(shí)驗(yàn)時(shí)，應(yīng)當(dāng)選擇海底平坦、深度穩(wěn)定的海域，盡量精確地測(cè)量海水深度，以此來(lái)更好地保證流速反演結(jié)果的準(zhǔn)確性。

表3 不同海深仿真反演結(jié)果Table 3 Results of simulation inversions with different sea depths

同樣在實(shí)驗(yàn)環(huán)境中，兩個(gè)陣列間的距離也往往受實(shí)際布放原因而產(chǎn)生誤差，通過(guò)改變仿真中兩點(diǎn)之間的距離，比較了距離擾動(dòng)下流速的估計(jì)結(jié)果仿真的結(jié)果。表4 列出了在仿真條件正負(fù)100 m 距離擾動(dòng)下對(duì)反演結(jié)果的影響，其中?d為正表示計(jì)算時(shí)距離大于實(shí)際距離?？梢钥吹郊词咕嚯x有著100 m的擾動(dòng)，正反向聲傳播時(shí)間有著較大誤差，但通過(guò)計(jì)算時(shí)延差，兩者的誤差可以消除。在實(shí)驗(yàn)中通常距離不確定性遠(yuǎn)小于100 m，因此距離上的失配對(duì)流速估計(jì)精度影響很小，一般情況下可以忽略。

表4 不同距離仿真反演結(jié)果Table 4 Results of simulation inversions with different distances

4 結(jié)論

本文提出了一種利用淺海海洋環(huán)境噪聲互相關(guān)函數(shù)估計(jì)海流流速的方法。本方法基于聲層析原理，采用波束形成增加能量積累，可以在2 h 時(shí)間尺度內(nèi)提取海洋環(huán)境噪聲互相關(guān)函數(shù)；計(jì)算獲得正反兩個(gè)方向上的格林函數(shù)及其時(shí)間到達(dá)結(jié)構(gòu)，反演估計(jì)出每2 h 的海流流速，從而可以觀測(cè)流速較長(zhǎng)期的變化趨勢(shì)；流速反演估計(jì)值和變化趨勢(shì)符合歷史觀測(cè)研究的結(jié)果，且和理論仿真計(jì)算結(jié)果一致性好；海水深度和距離失配誤差對(duì)流速估計(jì)精度影響小，可忽略?；诼晫游鲈?，利用海洋環(huán)境噪聲可以有效反演估計(jì)海流流速。

本文沒(méi)有使用聲學(xué)多普勒聲速儀[18]對(duì)海流流速進(jìn)行實(shí)時(shí)同步的測(cè)量，后續(xù)的實(shí)驗(yàn)應(yīng)當(dāng)實(shí)測(cè)海流流速，用以對(duì)比論證實(shí)驗(yàn)反演結(jié)果的精度。此外，實(shí)驗(yàn)使用的水聽(tīng)器采樣率為5000 Hz，也就是時(shí)間采樣間隔為0.2 ms，因此目前只能對(duì)海流流速的變化趨勢(shì)進(jìn)行大致觀測(cè)，后續(xù)實(shí)驗(yàn)將通過(guò)使用更高采樣率的設(shè)備，來(lái)獲取更精確的聲傳播時(shí)間，以此實(shí)現(xiàn)對(duì)海流流速的高精度估計(jì)。