戎澤存，胡長青，趙 梅

(1.中國科學(xué)院聲學(xué)研究所東海研究站，上海201815；2.中國科學(xué)院大學(xué)，北京100049)

0 引 言

海洋環(huán)境噪聲作為背景場，在很長一段時(shí)間里是被當(dāng)作干擾項(xiàng)來看待，而事實(shí)上海洋環(huán)境噪聲中包含著豐富的信息，因此它具有重要的研究價(jià)值。海洋環(huán)境噪聲的組成比較復(fù)雜，學(xué)界普遍認(rèn)為航船噪聲和風(fēng)成噪聲是兩大主要組成部分，前者主要影響的是中低頻段噪聲(10～500 Hz)，后者主要影響的是較高頻段噪聲(500～25 000 Hz)。除了以上兩大主要組成外，在低頻段(1～100 Hz)，潮、涌、浪、大尺度湍流以及遠(yuǎn)處的地震、風(fēng)暴均對海洋環(huán)境噪聲產(chǎn)生貢獻(xiàn)[1]。

作為一種隨機(jī)過程，海洋環(huán)境噪聲場不是一成不變的，尤其是隨著航運(yùn)業(yè)的發(fā)展，交通航運(yùn)噪聲對環(huán)境噪聲的影響更加顯著，使近海海域內(nèi)的環(huán)境噪聲較之前發(fā)生較大變化，近海環(huán)境噪聲日益成為人們關(guān)注的重點(diǎn)。近年來的研究發(fā)現(xiàn)，低頻段海洋環(huán)境噪聲級呈上升趨勢，和經(jīng)典的Wenz譜相比，在不同頻段內(nèi)有著小到1～3 dB、大到10～12 dB的增加[2]。

目前水聲設(shè)備向著低頻、遠(yuǎn)程方向發(fā)展，作為低頻段海洋環(huán)境噪聲的主要噪聲源，海上航船對噪聲場的影響不容忽視。目前學(xué)界對風(fēng)成海洋噪聲的研究較多，而對于低頻段的海洋環(huán)境噪聲涉足相對較少。因此本文針對淺海環(huán)境，提出一種計(jì)算航船噪聲對海洋環(huán)境噪聲貢獻(xiàn)的算法，主要關(guān)注頻段為 50～400 Hz，并利用實(shí)驗(yàn)數(shù)據(jù)對該算法進(jìn)行了驗(yàn)證。

1 理論推導(dǎo)

1.1 聲傳播模型

海洋環(huán)境下聲傳播模型主要包括波數(shù)積分模型、簡正波模型、射線模型、拋物方程模型、有限差分模型(有限元模型)[3]。本文選用了拋物方程模型對聲場進(jìn)行計(jì)算。

在研究低頻段海洋環(huán)境噪聲時(shí)，對于低頻、淺海條件，射線模型不再適用，因此本文采用了拋物方程(Parabolic Equations,PE)模型[3]。此方法可以較為精確地計(jì)算聲壓。

拋物方程法是一種近似算法，將聲壓方程近似為拋物線型方程。采用柱坐標(biāo)(r,φ,z)進(jìn)行推導(dǎo)，等密度介質(zhì)的三維亥姆霍茲方程為

假設(shè)方位對稱，可得到：

式中，p(r,z)是聲壓，k0=ω/c0是參考波數(shù)，n(r,z)=c0/ c(r,z)是折射率。

根據(jù)Tappert的研究[4]，可以假設(shè)式(2)的解滿足以下形式：

將式(3)代入式(2)，并作遠(yuǎn)場近似，即 k0r?1，可得到：

引入近軸近似：

最終得到標(biāo)準(zhǔn)的窄角拋物線方程：

為求解式(6)所示的PE方程，推導(dǎo)一系列基于算子形式的拋物波方程。首先定義如下算符：

則式(6)可重寫為

這里，求解方程的關(guān)鍵在于對Q算符的數(shù)學(xué)近似。距離相關(guān)聲學(xué)模型(Range-dependent Acoustic Model,RAM)程序作為一種聲場計(jì)算程序，它基于分步帕德(Padé)算法，對拋物方程組進(jìn)行求解。在數(shù)學(xué)中，Padé近似值是一個(gè)給定階數(shù)的有理函數(shù)對一個(gè)函數(shù)的“最佳”近似值[9-10]。應(yīng)用到實(shí)際場景時(shí)，將聲速梯度、海底聲速、海底底質(zhì)密度和吸收系數(shù)等環(huán)境參數(shù)作為輸入，利用RAM程序?qū)β晥鲞M(jìn)行計(jì)算，得到不同頻率聲波在特定海洋環(huán)境下的聲傳播特性。本文將航船視為點(diǎn)源，利用 RAM程序計(jì)算淺海中航船的傳播損失LT。在海洋環(huán)境中，聲波在傳播過程中的損失主要包括[5]：擴(kuò)展損失、吸收損失、散射、海底反射損失。在這里假定海面為理想界面，不考慮海面的散射。

通常深海環(huán)境下可不考慮海底反射損失。但在淺海環(huán)境下，除了擴(kuò)展損失和海水吸收損失，海底反射損失是一項(xiàng)關(guān)鍵因素。本文關(guān)注的是淺海海洋環(huán)境，因此在利用 RAM進(jìn)行計(jì)算時(shí)，必須將海底損失考慮在內(nèi)，以得到聲波傳播過程中的總傳播損失。

在利用RAM程序計(jì)算聲場時(shí)，還要注意以下兩點(diǎn)[3]：

(1)首先要建立一個(gè)在深度上有限的求解區(qū)域(0≤z≤zmax)，它包括實(shí)際物理域(海水和海底)及“海綿”層，“海綿”層是為了避免在下部計(jì)算邊界上產(chǎn)生附加反射，減少誤差；“海綿”層底部深度zmax=4H/3，其中H是傳播路徑上實(shí)際物理域的最大深度。

(2)其次要選擇合適的網(wǎng)格尺寸(Δr，Δz)，Δr為距離步長，Δz為深度步長。在淺海環(huán)境中，深度步長需滿足：Δz≤λ/4，其中λ是聲波波長；距離步長需滿足：2Δz≤Δr≤5Δz。為了獲得精確的數(shù)值解，唯一的辦法是按一定規(guī)則逐步減小Δr和Δz，對解(傳播損失)的收斂性仔細(xì)的檢驗(yàn)。若Δr和Δz的值設(shè)置不合適，可能會(huì)導(dǎo)致不同頻率下的傳播損失計(jì)算值產(chǎn)生誤差，因此需要針對不同頻率分別進(jìn)行收斂性檢驗(yàn)。

1.2 航船噪聲計(jì)算方法

為了研究海洋環(huán)境噪聲場的水平方向分布，以接收點(diǎn)為中心，將海洋劃分為若干楔形區(qū)域，如圖1(a)所示。圖1中xOy平面代表海面；從海面開始往下直至接收點(diǎn)所在平面，接收點(diǎn)記為R，D為接收點(diǎn)的深度。在z=D處取垂直于z軸的截面，并繪制扇區(qū)分布如圖1(b)所示。航船位置按劃分的扇區(qū)分組，扇區(qū)均勻分布，每個(gè)扇區(qū)作為一個(gè)方向。例如，劃分間隔若為2°，整個(gè)區(qū)域?qū)澐譃?80個(gè)扇區(qū)。扇區(qū)劃分越密，誤差越小。扇區(qū)截面的半徑由選取的計(jì)算范圍決定，選取的范圍越大，半徑越大，黑色散點(diǎn)代表航船在扇區(qū)截面的投影。利用本文算法計(jì)算時(shí)，可提前計(jì)算目標(biāo)海域的聲場，得到傳播損失隨傳播距離的變化，結(jié)合傳播損失和航船聲源級確定扇區(qū)的半徑。

圖1 楔形區(qū)域及扇區(qū)劃分分布圖Fig.1 Division map of wedge-shaped area and sector area

對某一個(gè)扇區(qū)的接收聲壓pR，可按式(9)和式(10)進(jìn)行計(jì)算[6]：

式(10)中，m表示航船序號，m=1,2,3,··,M，M為該扇區(qū)內(nèi)船的數(shù)量，各個(gè)符號含義如下：

(1)rm為該扇區(qū)內(nèi)序號為m的航船到接收點(diǎn)的距離；

(2)rc為該扇區(qū)內(nèi)航船與接收點(diǎn)的最大距離，分別計(jì)算該扇區(qū)內(nèi)的每艘航船到接收點(diǎn)的距離，取最大值即為rc；

(3)Lsm是該扇區(qū)內(nèi)航船 m在某一頻率下的聲源級，計(jì)算方法在1.3節(jié)討論；

(4)ωm指該扇區(qū)內(nèi)航船m對應(yīng)的權(quán)重系數(shù)，由式(11)得到：

按式(9)、(10)計(jì)算得到不同頻率的接收聲壓pR后，進(jìn)而可得到接收點(diǎn)處水聽器的接收聲壓譜密度DRL(f)[7]，計(jì)算公式為

式中：參考聲壓pref=1μ Pa。對其他扇區(qū)可重復(fù)進(jìn)行上述計(jì)算流程。

1.3 航船噪聲源

對航船噪聲源的模型，目前仍缺乏合適的數(shù)學(xué)與物理描述，一般仍采用經(jīng)驗(yàn)公式或半經(jīng)驗(yàn)公式描述。Ross[1]提出典型商船輻射噪聲的回歸公式，后來Hamson[7]更新了Ross的航船輻射噪聲的回歸方程，航船輻射噪聲功率譜密度如式(14)[1,7]所示：

式(14)各參數(shù)代表的意義如下：

(1)cv和cl分別代表航船速度和航船長度的羅斯冪律系數(shù)，cv=6，cl=2。

(2)v，l分別是船舶的速度和長度。v0和l0分別是航船參考速度和參考長度，v0=12 kn(1kn=1.852 km·h-1)、l0=300 ft(1ft=0.3048 m)。

(3)df是隨頻率變化的修正量，由式(15)給出，頻率f的單位是Hz[7]：

(4)dl是隨航船長度變化的修正量，由式(16)給出：

(5)Ls0(f)是僅與頻率相關(guān)的分量，由式(17)給出[7]：

利用式(14)，結(jié)合獲取的航船自動(dòng)識別系統(tǒng)(Automatic Identification System,AIS)數(shù)據(jù)(AIS數(shù)據(jù)將在2.1節(jié)進(jìn)行說明)，可以得到某一頻率的航船輻射噪聲功率譜密度Ls。

2 AIS數(shù)據(jù)介紹及預(yù)處理

2.1 AIS數(shù)據(jù)介紹

利用第1節(jié)構(gòu)建的算法，結(jié)合AIS數(shù)據(jù)，進(jìn)行了仿真計(jì)算。航船的位置、航速、船長等信息來源于船舶自動(dòng)識別系統(tǒng)(AIS)，以上數(shù)據(jù)作為算法的輸入，需要進(jìn)行一定的預(yù)處理，在本節(jié)中將對此進(jìn)行闡述。

圖2所示的是仿真所考察的航船AIS數(shù)據(jù)分布范圍示意圖，其中中心點(diǎn)經(jīng)緯度為[123.16°E,30.02°N]，4個(gè)頂點(diǎn)的坐標(biāo)分別為：[122.16°E,29.02°N]、[124.16°E,29.02°N]、[124.16°E,31.02°N]、[122.16°E,31.02°N]。該區(qū)域經(jīng)度跨度、緯度跨度均為2°。

航船AIS數(shù)據(jù)主要包括船名、海事移動(dòng)服務(wù)識別參數(shù)(Maritime Mobile Service Identity,MMSI)、船舶類型、船舶狀態(tài)、船長、船寬、吃水、航向、速度、經(jīng)度、緯度、時(shí)間戳等信息。文中船名和MMSI號已被隱去，必要時(shí)以手動(dòng)編號代替。

圖2 航船AIS數(shù)據(jù)分布范圍Fig.2 Distribution range of ship AIS data

表1中代表的是某航船在兩個(gè)時(shí)刻對應(yīng)的航船信息，包括船長、船速、吃水深度、經(jīng)緯度和記錄時(shí)間等，在兩個(gè)時(shí)間點(diǎn)的航速和航向不同。表中UTC代表格林威治標(biāo)準(zhǔn)時(shí)間，與北京時(shí)間相差8個(gè)小時(shí)，“UTC+8”即為北京時(shí)間。

表1 某航船AIS數(shù)據(jù)Table 1 AIS data of a ship

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

對于得到的AIS數(shù)據(jù)進(jìn)行預(yù)處理，包括：

(1)剔除存在明顯錯(cuò)誤的數(shù)據(jù)，剔除航船中發(fā)動(dòng)機(jī)已關(guān)閉的船舶；

(2)篩選符合時(shí)間范圍內(nèi)的數(shù)據(jù)，利用 Matlab軟件保存到數(shù)據(jù)庫；

(3)利用航船的船速、航向?qū)?shù)據(jù)進(jìn)行插值處理，作為數(shù)據(jù)的補(bǔ)充；

(4)查看地圖可發(fā)現(xiàn)，島礁是一個(gè)必須考慮的因素，在某一個(gè)時(shí)刻下，如果航船與接收點(diǎn)的連線上出現(xiàn)了島礁，可能會(huì)出現(xiàn)衍射(繞射)。現(xiàn)已證明，聲衍射的強(qiáng)弱和障礙物的大小與聲波波長的比值密切相關(guān)，其關(guān)系式為

式中：k為聲波的波數(shù)；f、c、λ分別為聲波的頻率、速度和波長；a是代表障礙物尺度的量，如為球體a就代表球體的半徑。當(dāng)ka?1時(shí)衍射很弱，聲波幾乎只是沿著直線傳播以至于被球體擋住傳播路徑而在其背面形成明顯的聲影。當(dāng)島礁尺寸滿足ka?1時(shí)，則應(yīng)將此艘航船予以剔除。

經(jīng)預(yù)處理得到航船分布和各艘船的信息后，依據(jù)算法，可獨(dú)立計(jì)算各扇區(qū)內(nèi)船舶輻射噪聲級譜密度DRL(f)。

3 仿真計(jì)算

3.1 仿真水體環(huán)境

水體環(huán)境參數(shù)參考某次海洋環(huán)境噪聲實(shí)驗(yàn)。該次測量點(diǎn)位于圖2所示海域范圍內(nèi)，水深約70 m，聲速梯度由電導(dǎo)率、溫度、深度(Conductivity,Temperature,Depth,CTD)系統(tǒng)測得，如圖3所示。用側(cè)掃聲吶測量海底，發(fā)現(xiàn)該海域海底較為平坦。同時(shí)通過對海底底質(zhì)的采樣，可知其主要成分接近粉砂質(zhì)砂。

圖3 實(shí)驗(yàn)海域聲速梯度剖面圖Fig.3 Acoustic velocity profile in the experimental area

3.2 水平方向分布分析

取某一時(shí)刻海域內(nèi)的航船數(shù)據(jù)作為輸入，包括航船經(jīng)緯度、長度、速度、吃水深度。預(yù)處理之后，運(yùn)用1.3節(jié)的方法進(jìn)行仿真計(jì)算。圖4即為該時(shí)刻預(yù)處理之后的航船位置分布圖，將此位置分布下的航船數(shù)據(jù)作為輸入，進(jìn)行航船噪聲的仿真計(jì)算。

按1/3倍頻程計(jì)，選定中心頻率，分別計(jì)算頻段內(nèi)若干頻點(diǎn)的功率譜密度值，并取平均得到的結(jié)果如圖5所示。

從圖5中可以看出：

(1)在不同頻率下，航船噪聲功率譜密度水平分布圖形狀不同。這是由于不同頻率的聲源在實(shí)驗(yàn)海域內(nèi)的傳播損失曲線不同：給定兩個(gè)不同頻率、相同聲源級的聲源，兩者的傳播損失隨距離的變化不同，導(dǎo)致不同頻率下的水平分布圖不同。

圖4 某一時(shí)刻的航船位置分布圖Fig.4 Distribution of ship positions at a certain time

圖5 不同中心頻率航船噪聲功率譜密度的水平方向分布圖Fig.5 Horizontal distributions of the ship noise power spectral density at different center frequencies

(2)圖5(a)～5(j)中：6°～12°、32°～36°、50°～52°、70°～72°、152°～154°、-152°～-122°、-118°～-30°、-18°～-16°和-12°～-8°方向的水平分布圖為“包絡(luò)”，其余為“線譜”。

進(jìn)一步分析數(shù)據(jù)，計(jì)算方法不變，從圖5(a)～5(h)中提取相應(yīng)方向的數(shù)據(jù)形成的-116°方向、128°方向、-36°方向、34°方向上在這一時(shí)刻的航船信息，結(jié)果如表2所示。表2說明在這一時(shí)間點(diǎn)某一方向上有幾艘船，以及每艘船的船長、船速、吃水深度信息，在-116°方向上有兩艘船：A1，A2；在-36°方向上有四艘船：B1，B2，B3，B4；而在128°方向和34°方向上均只有一艘船，分別是D和E。

4個(gè)方向上航船噪聲功率譜密度隨頻率的變化趨勢如圖6(a)～6(b)所示。

依頻率提取多個(gè)方向上的數(shù)據(jù)，得出如下結(jié)論:

(1)若扇區(qū)內(nèi)存在多艘航船，各艘航船的航行狀態(tài)不同，受到航船噪聲疊加的影響，不同方向上航船噪聲功率譜密度隨頻率的變化趨勢差別較大，如圖6(a)、6(b)所示。

(2)若扇區(qū)內(nèi)僅存在一艘航船，則航船噪聲功率譜密度隨頻率的變化趨勢較為穩(wěn)定，圖形走勢基本一致，如圖6(c)、6(d)所示。

表2 不同方向上的航船信息Table 2 Ship information in different directions

3.3 功率譜密度分析

采用本文算法，航船AIS信息輸入和水體環(huán)境參數(shù)與3.1節(jié)相同，得到了某一時(shí)刻的航船噪聲功率譜密度仿真結(jié)果。該時(shí)刻為2016年7月6日9時(shí)0分0秒，對不同方向的噪聲功率譜密度取平均，仿真結(jié)果如圖7中的星點(diǎn)所示，計(jì)算的頻點(diǎn)包括50、63、80、100、125、160、200、250、315、400 Hz。

圖6 不同方向航船噪聲功率譜密度Fig.6 The ship noise power spectrum densities in different directions

圖7 對某一時(shí)刻不同方向航船噪聲功率譜密度取平均的仿真結(jié)果Fig.7 Simulation result of the ship noise power spectral density averaged in different directions at a certain time

對于海洋環(huán)境噪聲的研究，Wenz噪聲譜級圖是一個(gè)重要的參考。它是由Wenz等[8]總結(jié)出的經(jīng)典海洋環(huán)境噪聲功率譜級圖。經(jīng)典Wenz噪聲譜級圖中包含了航船、風(fēng)、降雨、地震和爆炸等不同因素對環(huán)境噪聲的貢獻(xiàn)，如圖8所示：橫坐標(biāo)代表頻率(以倍頻程記，頻率從1 Hz～100 kHz)，縱坐標(biāo)代表海洋環(huán)境噪聲譜級。Wenz噪聲譜級圖中，不同因素貢獻(xiàn)的噪聲頻段不同，圖中各曲線代表的影響因素如表 3 所示。航船對噪聲的貢獻(xiàn)主要是10 Hz～1 kHz頻段，圖8中描述了強(qiáng)航道、淺海通常航道、深海通常航道等條件下的噪聲譜級。

圖8 Wenz噪聲功率譜級圖[8]Fig.8 Wenz noise power spectrum diagram[8]

為了驗(yàn)證本文算法的可靠性，將Wenz噪聲譜級圖中描述淺海通常航道的曲線取出，與圖7所示的仿真計(jì)算值同時(shí)繪制在圖9中。圖9中Wenz曲線 1、2分別代表淺海通常航道(Usual Traffic Noise-Shallow Water)的上、下限曲線。由圖9可以看出，仿真計(jì)算值與描述淺海通常航道的曲線2較為符合，在淺海通常航道的上、下限曲線之間，這也與水聽器和航船距離較遠(yuǎn)這一事實(shí)符合。實(shí)驗(yàn)所用水聽器為無指向性水聽器，布放位置為圖4中的接收點(diǎn)位置，深度為 41.03 m，水聽器與航船的距離從36.69～147.46 km不等。

表3 Wenz噪聲功率譜級圖說明Table 3 Description of Wenz noise power spectrum diagram

為了進(jìn)一步驗(yàn)證算法的可靠性，對實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了分析。圖10中連續(xù)曲線代表的是在2016年7月6日8時(shí)59分30秒～9時(shí)30秒的時(shí)間段內(nèi)，實(shí)測數(shù)據(jù)的1/3倍頻程功率譜密度。由圖10可見，在100～400 Hz頻段內(nèi)，仿真值與實(shí)驗(yàn)值較為符合。

與圖 10不同，現(xiàn)選取 10個(gè)時(shí)間點(diǎn)(間隔為1 min)，對10個(gè)時(shí)間點(diǎn)重復(fù)上述計(jì)算，并取平均，結(jié)果如圖11中星點(diǎn)所示。圖11中的實(shí)驗(yàn)值是對共計(jì)10 min長度的實(shí)驗(yàn)噪聲數(shù)據(jù)進(jìn)行分析得到的。由圖10、11可以發(fā)現(xiàn)，其結(jié)果與圖10中星點(diǎn)代表的仿真計(jì)算值結(jié)果較為一致，在100～400 Hz頻段內(nèi)兩者較為符合。

圖9 某一時(shí)刻航船噪聲功率譜密度仿真結(jié)果(圖7)與典型淺海航道的Wenz譜曲線對比Fig 9 Comparison between the simulation result of ship noise power spectral density(Fig.7)and the Wenz spectrums in typical shallow sea channels

圖10 某一時(shí)刻航船噪聲功率譜密度仿真結(jié)果(圖7)與實(shí)驗(yàn)值對比Fig.10 Comparison between the simulation result of ship noise power spectral density(Fig.7)and the experimental values

圖11 航船噪聲功率譜密度仿真值(某一時(shí)間段內(nèi)取平均)與實(shí)驗(yàn)值對比Fig 11 Comparison between the simulation result of ship noise power spectral density(averaging over a period of time)and the experimental value

本文中 50～400 Hz頻段內(nèi)的仿真結(jié)果與經(jīng)典Wenz譜總體上較為符合，在Wenz曲線淺海通常航道的上下限曲線之內(nèi)。但50～100 Hz頻段內(nèi)的仿真結(jié)果略低于實(shí)測環(huán)境噪聲數(shù)據(jù)。分析實(shí)驗(yàn)海域附近環(huán)境發(fā)現(xiàn)，由于實(shí)驗(yàn)海域?yàn)榻?，環(huán)境噪聲源種類繁多，噪聲源除了航船以外，還有工業(yè)活動(dòng)、潮、涌、浪、遠(yuǎn)處風(fēng)暴等，并且這些聲源對低頻段的影響較大，這也是誤差產(chǎn)生的原因之一。

4 結(jié) 論

本文提出了一種計(jì)算航船對近海海洋環(huán)境噪聲貢獻(xiàn)的算法。利用本文算法，結(jié)合AIS數(shù)據(jù)和實(shí)際海域水文參數(shù)，可獲取接收點(diǎn)處航船噪聲的水平方向分布特點(diǎn)，并可初步定量分析航船噪聲對海洋環(huán)境噪聲的貢獻(xiàn)。本文主要關(guān)注的頻段為 50～400Hz，理論分析和數(shù)值計(jì)算均反映出了由于航船的非均勻分布，噪聲場在不同頻率下呈現(xiàn)明顯的非均勻分布特征這一特點(diǎn)。對噪聲功率譜密度進(jìn)行了定量分析，并通過與經(jīng)典Wenz譜和實(shí)測噪聲數(shù)據(jù)進(jìn)行比對，在約50～400 Hz時(shí)仿真值和實(shí)驗(yàn)值吻合較好，驗(yàn)證了算法的可靠性。同時(shí)，在應(yīng)用時(shí)需要盡可能多地獲取環(huán)境參數(shù)，以降低數(shù)值計(jì)算時(shí)的誤差。