一種計算航船對近海海洋環(huán)境噪聲貢獻的算法

戎澤存，胡長青，趙 梅

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

0 引 言

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

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

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

1 理論推導

1.1 聲傳播模型

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

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

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

假設方位對稱，可得到：

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

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

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

引入近軸近似：

最終得到標準的窄角拋物線方程：

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

則式(6)可重寫為

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

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

在利用RAM程序計算聲場時，還要注意以下兩點[3]：

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

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

1.2 航船噪聲計算方法

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

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

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

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

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

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

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

(4)ωm指該扇區(qū)內航船m對應的權重系數，由式(11)得到：

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

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

1.3 航船噪聲源

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

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

(1)cv和cl分別代表航船速度和航船長度的羅斯冪律系數，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)是僅與頻率相關的分量，由式(17)給出[7]：

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

2 AIS數據介紹及預處理

2.1 AIS數據介紹

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

圖2所示的是仿真所考察的航船AIS數據分布范圍示意圖，其中中心點經緯度為[123.16°E,30.02°N]，4個頂點的坐標分別為：[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ū)域經度跨度、緯度跨度均為2°。

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

圖2 航船AIS數據分布范圍Fig.2 Distribution range of ship AIS data

表1中代表的是某航船在兩個時刻對應的航船信息，包括船長、船速、吃水深度、經緯度和記錄時間等，在兩個時間點的航速和航向不同。表中UTC代表格林威治標準時間，與北京時間相差8個小時，“UTC+8”即為北京時間。

表1 某航船AIS數據Table 1 AIS data of a ship

2.2 數據預處理

對于得到的AIS數據進行預處理，包括：

(1)剔除存在明顯錯誤的數據，剔除航船中發(fā)動機已關閉的船舶；

(2)篩選符合時間范圍內的數據，利用 Matlab軟件保存到數據庫；

(3)利用航船的船速、航向對數據進行插值處理，作為數據的補充；

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

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

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

3 仿真計算

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

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

圖3 實驗海域聲速梯度剖面圖Fig.3 Acoustic velocity profile in the experimental area

3.2 水平方向分布分析

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

按1/3倍頻程計，選定中心頻率，分別計算頻段內若干頻點的功率譜密度值，并取平均得到的結果如圖5所示。

從圖5中可以看出：

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

圖4 某一時刻的航船位置分布圖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°方向的水平分布圖為“包絡”，其余為“線譜”。

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

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

依頻率提取多個方向上的數據，得出如下結論:

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

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

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

3.3 功率譜密度分析

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

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

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

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

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

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

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

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

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

圖9 某一時刻航船噪聲功率譜密度仿真結果(圖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 某一時刻航船噪聲功率譜密度仿真結果(圖7)與實驗值對比Fig.10 Comparison between the simulation result of ship noise power spectral density(Fig.7)and the experimental values

圖11 航船噪聲功率譜密度仿真值(某一時間段內取平均)與實驗值對比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頻段內的仿真結果與經典Wenz譜總體上較為符合，在Wenz曲線淺海通常航道的上下限曲線之內。但50～100 Hz頻段內的仿真結果略低于實測環(huán)境噪聲數據。分析實驗海域附近環(huán)境發(fā)現，由于實驗海域為近海，環(huán)境噪聲源種類繁多，噪聲源除了航船以外，還有工業(yè)活動、潮、涌、浪、遠處風暴等，并且這些聲源對低頻段的影響較大，這也是誤差產生的原因之一。

4 結 論

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