鄭廣學(xué)，樸勝春1,2,，祝捍皓，張海剛1,2,，李楠松1,2,

(1.哈爾濱工程大學(xué) 水聲技術(shù)重點(diǎn)實(shí)驗(yàn)室，哈爾濱 150001；2.海洋信息獲取與安全工信部重點(diǎn)實(shí)驗(yàn)室(哈爾濱工程大學(xué))，工業(yè)和信息化部，哈爾濱 150001；3.哈爾濱工程大學(xué) 水聲工程學(xué)院，哈爾濱 150001；4.浙江海洋大學(xué) 海洋科學(xué)與技術(shù)學(xué)院,浙江 舟山 316022；5.浙江大學(xué) 浙江省海洋觀測(cè)-成像試驗(yàn)區(qū)重點(diǎn)實(shí)驗(yàn)室，浙江 舟山 316021；6.中國(guó)科學(xué)院聲學(xué)研究所 聲場(chǎng)聲信息國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100190)

海洋水體環(huán)境參數(shù)與地聲參數(shù)共同構(gòu)成了水下聲傳播計(jì)算最重要的環(huán)境參數(shù)，相比水體環(huán)境參數(shù)，地聲參數(shù)更難直接測(cè)量，因而對(duì)其的有效獲取研究一直是水聲學(xué)中的熱點(diǎn)課題[1]。由于地聲參數(shù)的變化會(huì)對(duì)聲壓場(chǎng)的分布產(chǎn)生重要影響，因此可以考慮以水聽器接收到的復(fù)聲壓作為研究對(duì)象來反演地聲參數(shù)。相比于直接測(cè)量的方法，利用聲學(xué)方法反演可以快速、高效的獲取大面積海區(qū)的底質(zhì)參數(shù)，避免了人力、物力的浪費(fèi)，因此受到了廣泛的關(guān)注[2-5]。近年來，利用聲學(xué)方法反演地聲參數(shù)已取得了長(zhǎng)足的發(fā)展。各種反演方法不斷涌現(xiàn)，如：利用傳播損失的地聲參數(shù)反演方法[6]；利用艦船輻射噪聲等各類海洋噪聲信息的地聲參數(shù)反演方法[7-10]；利用波導(dǎo)頻散特性的反演方法等。但上述地聲參數(shù)反演方法均著重研究反演問題中正演模型的選擇，多仍采用傳統(tǒng)尋優(yōu)算法，如下降單純形法(downhill simplex，DS)、遺傳算法(genetic algorithm，GA)等對(duì)目標(biāo)函數(shù)進(jìn)行求解。但此類算法一般只能給出待反演參數(shù)的點(diǎn)估計(jì)，無法對(duì)參數(shù)的不確定性進(jìn)行進(jìn)一步描述，且在尋優(yōu)過程中對(duì)算法初始模型依賴很強(qiáng)，容易陷入局部最優(yōu)解。貝葉斯反演方法有效地估計(jì)最大后驗(yàn)概率(maximum a posteriori，MAP)模型參數(shù)，進(jìn)而從統(tǒng)計(jì)角度定性和定量分析參數(shù)反演結(jié)果的不確定性。

基于上述原因，本文以單水聽器接收到各距離點(diǎn)上的聲壓場(chǎng)為研究對(duì)象，基于貝葉斯理論，建立了針對(duì)淺?？紤]彈性海底下的地聲參數(shù)反演方法，反演對(duì)象包括聲速、聲速衰減和密度3類海底聲學(xué)參數(shù)。

1 反演方法

根據(jù)聲場(chǎng)信息估計(jì)未知海底參數(shù)的聲學(xué)反演方法一般由4部分組成：1)構(gòu)建聲傳播模型；2)建立合適的目標(biāo)函數(shù)；3)全局優(yōu)化算法；4)反演結(jié)果的不確定性分析。在本文的反演方法研究中，同樣針對(duì)上述4部分展開工作。

1.1 淺海聲場(chǎng)正演模型

淺海波導(dǎo)中，海底橫波聲速對(duì)低頻聲信號(hào)的傳播影響不可忽略[11]，本文研究中選擇了如圖1所示的淺海參數(shù)化模型。其中，z=0表征海面，z軸表深度，向下為正，r正軸表示聲場(chǎng)向外傳播方向。設(shè)水深為H、聲源位于水深zs處、密度ρ1、聲速c1；海底為半無限彈性介質(zhì)，各向同性?？v波聲速、橫波聲速、密度、縱波聲速衰減和橫波聲速衰減分別用cp、cs、ρb、αp和αs表示，該5項(xiàng)海底參數(shù)也是本文反演時(shí)的研究對(duì)象。

圖1 單層彈性海底的淺海參數(shù)化模型Fig.1 Parametric model in shallow sea with elastic bottom

圖1模型下，流體層中的各位置聲壓p(r,z)為：

(1)

(2)

式中：ξ為水平波數(shù)；r為水平距離；J0為貝塞爾函數(shù)；文獻(xiàn)[12]給出了A、B、C、β等參數(shù)的具體推導(dǎo)過程及表示式。對(duì)式(1)中聲壓場(chǎng)p(r,z)的數(shù)值計(jì)算，通常采用簡(jiǎn)正波方法(normal mode method，NMM)及快速場(chǎng)方法(fast field method，F(xiàn)FM)[13]?？紤]此模型為淺海環(huán)境，F(xiàn)FM更適合用于聲場(chǎng)的快速計(jì)算，故本文選擇FFM完成對(duì)上述參數(shù)化模型下聲壓場(chǎng)p(r,z)的數(shù)值計(jì)算。

在完成對(duì)聲壓場(chǎng)p(r,z)數(shù)值計(jì)算的基礎(chǔ)上，需建立恰當(dāng)?shù)哪繕?biāo)函數(shù)來衡量實(shí)測(cè)聲壓p(r,z)與理論預(yù)報(bào)聲壓p′(r,z)間的適配程度，并通過對(duì)其尋優(yōu)求解便可實(shí)現(xiàn)該模型下各海底參數(shù)的獲取。

1.2 目標(biāo)函數(shù)與優(yōu)化算法

貝葉斯方法用于反演問題時(shí)，觀測(cè)數(shù)據(jù)向量d和模型參數(shù)向量m被視為隨機(jī)的。根據(jù)貝葉斯理論有:

P(m|d)=P(d|m)P(m)/P(d)

(3)

其中，P(m|d)為后驗(yàn)概率密度(posterior probability density，PPD)。條件概率密度P(d|m)定義為似然函數(shù)L(m)，為某一(固定)測(cè)量數(shù)據(jù)d下m的函數(shù)，代表數(shù)據(jù)提供的信息，用于衡量在模型參數(shù)m下模型預(yù)報(bào)結(jié)果和數(shù)據(jù)的適配程度。P(m)為先驗(yàn)概率密度，表示模型參數(shù)先驗(yàn)信息，P(d)為測(cè)量數(shù)據(jù)的概率密度函數(shù)[14-15]。因此，后驗(yàn)概率密度既包含數(shù)據(jù)信息，也包含模型參數(shù)的先驗(yàn)信息。且P(d)與m無關(guān)，可以看作1個(gè)常數(shù)，所以式(3)可化為：

P(m|d)∝P(d|m)P(m)

(4)

根據(jù)貝葉斯理論，目標(biāo)函數(shù)是在高斯數(shù)據(jù)誤差的假設(shè)下基于似然函數(shù)建立的，似然函數(shù)作為某一指定測(cè)量數(shù)據(jù)下模型參數(shù)的函數(shù)，代表了數(shù)據(jù)不確定性分布[14-15]。本文使用的目標(biāo)函數(shù)為：

E(m)=ln[B(m)|p|2]

(5)

式中B(m)表示歸一化的Bartlett失配器：

(6)

式中：p表示單個(gè)水聽器接收到的聲壓場(chǎng)(實(shí)測(cè)場(chǎng))，p′(m)為待反演參數(shù)向量m的預(yù)報(bào)聲壓(拷貝場(chǎng))，上標(biāo)“+”代表共軛轉(zhuǎn)置。

從式(6)可以看出，B(m)取得最優(yōu)匹配時(shí)的條件為p(r,z)=p′(r,z)，同時(shí)，E(m)取得極小值。目標(biāo)函數(shù)的建立，將對(duì)m的反演問題轉(zhuǎn)化為在尋優(yōu)范圍內(nèi)對(duì)式(5)的極小值求解問題。根據(jù)貝葉斯理論，上述待反演參數(shù)的PPD將依據(jù)Metropolis準(zhǔn)則采樣[16]給出，并從統(tǒng)計(jì)角度對(duì)其進(jìn)行不確定性分析。

2 聲場(chǎng)正演模型仿真分析

本文在完成了聲場(chǎng)建模、確定了目標(biāo)函數(shù)及尋優(yōu)算法的基礎(chǔ)上，利用數(shù)值模擬，對(duì)比圖1模型下各海底參數(shù)與其反演結(jié)果，對(duì)研究方法進(jìn)行驗(yàn)證。

2.1 聲傳播特性分析

圖2給出了利用FFM計(jì)算所得圖1所示淺海環(huán)境下，1～1 500 m 100 Hz聲信號(hào)聲壓場(chǎng)的仿真結(jié)果，以傳播損失(transmission loss，TL)曲線的形式給出。

圖2 5項(xiàng)參數(shù)對(duì)TL的影響程度分析Fig.2 The impact of five parameters on transmission loss

仿真中，設(shè)聲源深度為zs=20 m，接收點(diǎn)zr位于水下10 m處，表1給出了仿真中各海洋環(huán)境參數(shù)。圖2(a)～(e)中的實(shí)線對(duì)應(yīng)利用表1中各參數(shù)真值的計(jì)算結(jié)果，5項(xiàng)地聲參數(shù)取討論值時(shí)的計(jì)算結(jié)果分別用點(diǎn)線(較小值)和虛線(較大值)表示。圖2(f)中應(yīng)用“距平”的定義[17]，視設(shè)定的環(huán)境參數(shù)真值計(jì)算所得傳播損失為均值，視各參數(shù)的討論值計(jì)算所得傳播損失為離散值來計(jì)算距平，用5項(xiàng)地聲參數(shù)討論值的距平大小來衡量各參數(shù)變化對(duì)傳播損失的影響程度。

表1中各參數(shù)的討論值均設(shè)定為真值偏移±10%，在此基礎(chǔ)上討論各參數(shù)變化對(duì)聲壓場(chǎng)的影響程度。從圖2(a)～(e)的中可以看出，5項(xiàng)參數(shù)的變化對(duì)TL曲線影響程度各不相同，從圖2(f)中5項(xiàng)參數(shù)討論值的距平大小可以看出，在各參數(shù)均偏移±10%的情況下，改變某一參數(shù)對(duì)TL曲線的影響程度從大到小依次為cp、cs、ρb、αp、αs。因此，5項(xiàng)海底參數(shù)對(duì)p(r,z)的影響程度可初步定性為：cp、cs>ρb>αp、αs。

2.2 目標(biāo)函數(shù)對(duì)各參數(shù)的敏感度分析

在完成了正演建模、聲壓場(chǎng)數(shù)值計(jì)算的基礎(chǔ)上，本文通過數(shù)值模擬仿真分析本文所給目標(biāo)函數(shù)對(duì)5項(xiàng)地聲參數(shù)的敏感度，并討論該目標(biāo)函數(shù)在反演中的可靠性。仿真中，海洋環(huán)境參數(shù)真值與各參數(shù)尋優(yōu)區(qū)間仍取表1中數(shù)值，取各參數(shù)設(shè)定真值下模擬計(jì)算得到p(r,z)為參考真值。本文采用固定變量法來分析目標(biāo)函數(shù)對(duì)單個(gè)地聲參數(shù)的敏感度時(shí)，即其他參數(shù)固定，只在討論區(qū)間內(nèi)改變某一參數(shù)計(jì)算聲壓場(chǎng)理論預(yù)報(bào)值p′(r,z)，由式(5)計(jì)算p(r,z)與p′(r,z)間的目標(biāo)函數(shù)值E(m) 作為敏感度評(píng)價(jià)標(biāo)準(zhǔn)。

表1 海洋環(huán)境仿真參數(shù)Table 1 The simulation parameters of ocean environment

圖3(a)～(e)分別給出了目標(biāo)函數(shù)E(m)隨表1中各單項(xiàng)參數(shù)變化時(shí)的數(shù)值變化規(guī)律。從圖3可以看出，在5項(xiàng)地聲參數(shù)討論范圍內(nèi)，目標(biāo)函數(shù)E(m)均只在其真值處取得極小值，且無偽峰或第2極小值，避免了局部最優(yōu)解的陷阱，證明了該目標(biāo)函數(shù)的可行性。但從圖3(f)可知目標(biāo)函數(shù)隨各參數(shù)變化時(shí)的波動(dòng)程度不一致，在5項(xiàng)參數(shù)的討論范圍內(nèi)，目標(biāo)函數(shù)值波動(dòng)范圍從大到小依次為cp、cs、ρb、αp、αs。因此，5項(xiàng)海底參數(shù)對(duì)p(r,z)的影響程度可定性為：cp>cs>ρb>αp>αs，進(jìn)一步證明了圖2的討論結(jié)果。

圖3 單參數(shù)敏感度Fig.3 The sensitivity of single parameter

2.3 仿真反演結(jié)果分析

在完成目標(biāo)函數(shù)對(duì)5項(xiàng)海底聲學(xué)參數(shù)敏感度的研究基礎(chǔ)上，本節(jié)通過數(shù)值仿真，分析本文所研究的反演方法對(duì)目標(biāo)函數(shù)尋優(yōu)結(jié)果的可靠性。海洋環(huán)境參數(shù)如表1所設(shè)，當(dāng)尋優(yōu)結(jié)束時(shí)，各參數(shù)尋優(yōu)結(jié)果與仿真真值的對(duì)比如圖4給出。圖中水平線段表示了各參數(shù)PPD的均值位置及其方差大小，垂直于橫坐標(biāo)直線段標(biāo)注仿真真值所在位置，各參數(shù)反演結(jié)果在表2中給出。

圖4 5項(xiàng)參數(shù)及其概率密度分布Fig.4 Five parameters and their probability density distribution

表2 尋優(yōu)結(jié)果Table 2 The results of optimization

從圖4和表2中的結(jié)果可以看出，各參數(shù)的仿真真值與反演結(jié)果吻合程度較高，均處于其尋優(yōu)結(jié)果的均值與方差的和或差范圍內(nèi)。且各參數(shù)的概率密度峰值的尖銳程度也反映了目標(biāo)函數(shù)對(duì)各參數(shù)的敏感程度，如圖4所示，目標(biāo)函數(shù)對(duì)各參數(shù)的敏感程度為：cp、cs>ρb>αp、αs，這也與2.1節(jié)、2.2節(jié)的討論結(jié)果相吻合。

為驗(yàn)證數(shù)值仿真中反演結(jié)果的可靠性，圖5給出了仿真環(huán)境參數(shù)下TL曲線與仿真反演結(jié)果計(jì)算所得TL曲線的對(duì)比圖。等間距取12處繪制誤差棒，并標(biāo)注反演結(jié)果計(jì)算所得TL曲線的誤差范圍。從圖中可以看出，仿真中TL曲線均處于反推TL曲線誤差范圍之內(nèi)，證明了該反演方法在數(shù)值仿真應(yīng)用中的可靠性。

圖5 TL曲線對(duì)比Fig.5 The comparison of TL

3 消聲水池縮比實(shí)驗(yàn)驗(yàn)證

3.1 水池實(shí)驗(yàn)

本文結(jié)合消聲水池縮比實(shí)驗(yàn)[18]，利用實(shí)驗(yàn)數(shù)據(jù)對(duì)本文所研究反演方法進(jìn)行討論。實(shí)驗(yàn)中選取質(zhì)地均勻、高硬度的塑料板模擬“半無限彈性海底”。

表3給出了實(shí)驗(yàn)中各設(shè)備的布放參數(shù)，其中聲源深度zs為87 mm，接收深度zr為84 mm，水深z為182 mm，水中聲速c1為1 450.212 m/s由測(cè)量水中溫度(11.15 ℃)后通過經(jīng)驗(yàn)公式求得。測(cè)量過程中，聲源固定不動(dòng)，發(fā)射信號(hào)為135 kHz的脈沖信號(hào)；接收水聽器固定在可移動(dòng)走架上，采集卡采樣率fs為20 MHz，單個(gè)測(cè)量點(diǎn)分別測(cè)量10次以減少隨機(jī)擾動(dòng)帶來的誤差；完成一點(diǎn)測(cè)量后，走架帶動(dòng)水聽器向遠(yuǎn)離聲源方向移動(dòng)2 mm(誤差不超過20 μm)至下1個(gè)測(cè)量點(diǎn)，實(shí)驗(yàn)中測(cè)量點(diǎn)共720個(gè)。測(cè)量中TL曲線如圖6(b)所示，圖6(a) 給出了測(cè)量中所接收到的第1～100道接收信號(hào)時(shí)域圖。

圖6 實(shí)驗(yàn)中接收到的信號(hào)Fig.6 The signals received in experiment

3.2 反演結(jié)果驗(yàn)證與分析

表3列出了各地聲參數(shù)的尋優(yōu)范圍，以及利用本文所研究反演方法對(duì)上述實(shí)驗(yàn)數(shù)據(jù)反演所得到的“海底”地聲參數(shù)值，同時(shí)還給出了利用GA對(duì)該實(shí)驗(yàn)數(shù)據(jù)尋優(yōu)的結(jié)果作為對(duì)比。

表3 尋優(yōu)結(jié)果Table 3 The results of optimization

圖7給出了算法終止時(shí)，各參數(shù)在尋優(yōu)區(qū)間上的分布及其概率密度分布，與圖4相仿，圖7中垂直于橫坐標(biāo)直線段標(biāo)注GA尋優(yōu)結(jié)果所在位置，水平線段表示了各參數(shù)取值PPD的均值位置及其方差大小。圖8定量給出了參數(shù)間的相關(guān)系數(shù)矩陣圖。

從表4、圖7和圖8的結(jié)果中可以看到：首先，各待反演參數(shù)的GA尋優(yōu)結(jié)果均處于均值與方差的和或差范圍內(nèi)，證明了利用貝葉斯反演方法獲得的結(jié)果與傳統(tǒng)反演方法得到結(jié)果相吻合。另外，方差的相對(duì)大小從一方面也反映出該參數(shù)反演結(jié)果不確定性的大小。其次敏感度較小的參數(shù)對(duì)目標(biāo)函數(shù)貢獻(xiàn)也較小，在圖8中表現(xiàn)為該參數(shù)與其他參數(shù)的相關(guān)性偏小，如αp和αs；換言之，較敏感參數(shù)間則表現(xiàn)出較明顯的相關(guān)性，如cp、cs、ρb，綜上，可以看出5項(xiàng)待反演參數(shù)的敏感度大小分別為cp、cs>ρb>αp、αs，這也與第2節(jié)結(jié)論一致。

圖7 算法終止時(shí)5項(xiàng)參數(shù)及其概率密度分布Fig.7 The five parameters and their probability density distribution at the termination of the algorithm

圖8 參數(shù)間相關(guān)矩陣Fig.8 The correlation matrix of parameters

為進(jìn)一步驗(yàn)證表4反演結(jié)果的可信度，圖9對(duì)比了利用反演結(jié)果計(jì)算得到的TL曲線和實(shí)測(cè)TL曲線；圖10對(duì)比了同一測(cè)量點(diǎn)(第470號(hào)接收點(diǎn))，反演所得時(shí)域波形與該點(diǎn)實(shí)際接收信號(hào)的時(shí)域波形圖(均進(jìn)行歸一化處理)。圖9中，實(shí)線為實(shí)測(cè)TL曲線，點(diǎn)線為應(yīng)用表4反演結(jié)果計(jì)算所得的TL曲線并等間距取12處繪制誤差棒，從圖中對(duì)比結(jié)果可以看到實(shí)驗(yàn)中TL曲線基本處于反推TL曲線誤差范圍之內(nèi)，證明了該反演方法在實(shí)際應(yīng)用中的可行性；圖10中，實(shí)線為第470號(hào)接收點(diǎn)實(shí)測(cè)時(shí)域波形，點(diǎn)線為應(yīng)用表4反演結(jié)果“預(yù)報(bào)”的第470號(hào)接收點(diǎn)時(shí)域波形，從圖中對(duì)比結(jié)果可以看到時(shí)域波形計(jì)算預(yù)報(bào)結(jié)果與實(shí)測(cè)信號(hào)波形匹配較好，若不考慮信號(hào)“拖尾”，相關(guān)系數(shù)可達(dá)0.89，證明了本反演方法在實(shí)際應(yīng)用中的可靠性。

圖9 實(shí)測(cè)TL與反演所得TL對(duì)比Fig.9 The comparison of TL between measured and inversed

圖10 時(shí)域波形對(duì)比Fig.10 The comparison of time domain waveform

4 結(jié)論

1)貝葉斯反演能夠從統(tǒng)計(jì)角度給出反映地聲參數(shù)向量綜合信息的后驗(yàn)概率來表征反演結(jié)果，同時(shí)還能夠分析各參數(shù)反演結(jié)果的不確定性，相比傳統(tǒng)只能給出待反演參數(shù)點(diǎn)估計(jì)的反演方法，反演信息更加全面、科學(xué)。

2)根據(jù)貝葉斯反演理論，利用本文給出的目標(biāo)函數(shù)，運(yùn)用Metropolis準(zhǔn)則采樣，可以較為準(zhǔn)確得到待反演參數(shù)的PPD。數(shù)值模擬和水池實(shí)驗(yàn)的結(jié)果均表明：PPD分布特性與GA尋優(yōu)結(jié)果吻合程度較好，利用反演參數(shù)計(jì)算得到的聲場(chǎng)特性與測(cè)量場(chǎng)基本一致，時(shí)域波形相關(guān)系數(shù)達(dá)0.89。

3)從5項(xiàng)地聲參數(shù)反演結(jié)果的不確定性大小來看，在仿真模擬及水池實(shí)驗(yàn)的環(huán)境下，縱波聲速cp、橫波聲速cs以及密度ρb的不確定性相對(duì)較小，對(duì)聲壓場(chǎng)的變化更敏感，5項(xiàng)參數(shù)依據(jù)敏感程度大小依次為cp>cs>ρb>αp>αs，聲場(chǎng)研究結(jié)果與反演結(jié)果均證實(shí)了這一研究結(jié)果。