吳偉文，任群言，鹿力成，馬 力

(1.中國科學院聲學研究所中國科學院水聲環(huán)境特性重點實驗室，北京 100190；2.中國科學院大學，北京 100049)

0 引言

聲信號在水中具有良好的傳播性能，因此常用來進行水下目標探測或信息傳遞。淺海環(huán)境中，聲傳播受海底的影響較大，準確的地聲參數(shù)有助于分析淺海環(huán)境中聲傳播特性，因此地聲參數(shù)估計是海洋聲學中的一個基礎研究問題。

地聲參數(shù)可以通過直接測量法[1]獲得，但此方法時間成本和經(jīng)濟代價都很高。反演算法如匹配場處理[2]是比直接測量法更有效率的方式。匹配場處理將地聲參數(shù)反演問題作為優(yōu)化問題，其目的是找到一個模型，此模型包含一組地聲參數(shù)(假設模型的真實地聲參數(shù)向量包含在參數(shù)集合中)，使前向模型產(chǎn)生的數(shù)據(jù)與觀測數(shù)據(jù)達到最優(yōu)匹配。

匹配場處理可分為兩步。(1)選擇目標函數(shù)，對觀測數(shù)據(jù)與前向模型生成數(shù)據(jù)之間的誤差進行衡量；(2)基于最小化目標函數(shù)的準則搜索最優(yōu)的參數(shù)向量。目標函數(shù)中待匹配的數(shù)據(jù)，一般可選擇為對地聲參數(shù)變化較敏感的物理量。例如，海底反射系數(shù)隨沉積層分層特性、聲速、密度的變化而明顯變化，因此海底反射系數(shù)[3]可以作為目標函數(shù)中待匹配的物理量。本文將海底反射系數(shù)作為待匹配的物理量，利用優(yōu)化算法搜索地聲參數(shù)向量。由于地聲反演問題具有維度高，局部最優(yōu)多，前向模型非線性等特點，所以參數(shù)搜索方式一般選遺傳算法[4]、模擬退火[4]、MCMC[5]采樣等全局算法。上述優(yōu)化算法一般稱之為批處理反演算法，不僅搜索過程費時而且一批數(shù)據(jù)只能生成一個反演結果，反演結果之間互相獨立，無法高效處理環(huán)境隨距離變化時的參數(shù)估計問題。

為了利用相鄰兩段距離之間地聲參數(shù)的聯(lián)系，減少運算量，可以將序貫尋優(yōu)算法如粒子濾波應用到地聲參數(shù)的反演中。在地聲參數(shù)隨距離變化緩慢的情況下粒子濾波算法已經(jīng)取得了良好的效果[6]。但是在地聲參數(shù)隨距離變化劇烈，比如沉積層分層個數(shù)改變的情況下，傳統(tǒng)粒子濾波可能會失效。為了解決此問題，Dettmer等[7]將帶橋接的粒子濾波應用到地聲參數(shù)反演中，通過橋接過渡劇烈變化的目標分布，以實現(xiàn)地聲參數(shù)隨距離變化劇烈時的反演，并通過仿真數(shù)據(jù)證明了此方法的有效性。本文在粒子濾波中增加橋接的同時改進了粒子的采樣方式，由傳統(tǒng)的固定方差采樣改進為隨機方差采樣，同時對沉積層厚度的采樣采用高維與低維混合采樣的方式。將此方法應用到地聲參數(shù)反演中，通過仿真分析得出，此方法在地聲參數(shù)隨距離變化劇烈的情況下仍可得到較穩(wěn)定的反演結果。

本文利用反射系數(shù)構造似然函數(shù)，分析其對沉積層厚度、聲速、密度的敏感性；給出了改進粒子濾波的算法流程；然后進行仿真數(shù)據(jù)處理；得到改進粒子濾波在地聲參數(shù)隨距離變化劇烈的情況下仍可得到較穩(wěn)定的反演結果的結論。

1 海底反射系數(shù)

1.1 反射系數(shù)計算

海底反射系數(shù)對沉積層厚度以及聲速都較為敏感，因此常作為地聲參數(shù)反演中待匹配的物理量。反射系數(shù)R為反射到海水的所有反射波與入射波的比值，基于射線聲學理論，平面波反射系數(shù)計算原理如圖1所示[8]，R的計劃公式為

圖1 沉積層數(shù)為1時的海底反射圖Fig.1 Bottom reflection diagram when the number of sediment layer is 1

其中：ρ1、ρ2、ρ3分別為海水、沉積層和基底的密度；c1、c2、c3分別為海水、沉積層和基底的聲速(當考慮海底聲衰減系數(shù)α2、α3時，海底聲速為復數(shù))；h2為沉積層厚度；R12為聲波由海水入射到沉積層的反射系數(shù)，T12為聲波由海水入射到沉積層的透射系數(shù)，R21為聲波由沉積層入射到海水的反射系數(shù)，T21為聲波由沉積層入射到海水的透射系數(shù)，R23為聲波由沉積層入射到基底的反射系數(shù)，T23為聲波由沉積層入射到基底的透射系數(shù)，R32為聲波由基底入射到沉積層的反射系數(shù)，T23為聲波由基底入射到沉積層的透射系數(shù)。相位延遲為

其中：k2=2πf/c2，f為聲源頻率,θ2為折射波方向(已知入射波方向θ1，聲速c1和c2，可根據(jù)斯涅爾定律計算θ2)。根據(jù)R12=-R21和T12T21=1-R212，反射系數(shù)R進一步可表示為

多層沉積層的反射系數(shù)可以依公式(3)遞推獲得，也可通過OASR工具箱[9]計算獲得。

1.2 敏感性分析

反射系數(shù)隨沉積層層厚，聲速，密度的變化而變化，本文用反射系數(shù)構造如式(4)所示的似然函數(shù)：

其中：xt為第t幀的地聲參數(shù)向量，包含沉積層深度、聲速、密度、衰減系數(shù)；yt為第t時刻的觀測值，對應于反演問題中的第t幀的反射系數(shù)觀測向量(包含N個掠射角)；g為測量函數(shù)。研究似然函數(shù)對沉積層聲速、密度、衰減系數(shù)、厚度的敏感性，對反演的可行性以及反演結果分析具有重要作用。下面在四層沉積層的情況下，對似然函數(shù)進行敏感性分析。

令聲源頻率分別為500 Hz和1 000 Hz，仿真參數(shù)如表1所示。

表1 仿真信道參數(shù)Table 1 Parameters of the simulation channel

圖2 Φ隨ci、ρi、αi、hi各參數(shù)變化的曲線圖Fig.2Variation curves of Φ with parameters ci、ρi、αi、hi

可以清楚地看到，Φ對聲速和沉積層深度敏感，對密度、衰減系數(shù)不敏感。同時可以看到Φ隨h變化時，局部最優(yōu)明顯。為了防止算法在對沉積層深度的搜索時陷入局部最優(yōu)，本文提出了沉積層厚度混合搜索方式，此方法將在第2.2節(jié)介紹。

2 改進粒子濾波

2.1 帶橋接的粒子濾波

在隨距離變化的地聲反演問題中，可以將地聲參數(shù)分為很多幀。本文的幀是指一小段范圍，在此范圍內(nèi)地聲參數(shù)可以近似認為是不變的。由于每一幀之間地聲參數(shù)具有相關性，所以可以將粒子濾波等序貫尋優(yōu)方式應用到地聲反演中。根據(jù)粒子濾波原理建立狀態(tài)方程和觀測方程分別為

式中：xt是第t時刻的狀態(tài)，在地聲反演問題中，即為第t幀的地聲參數(shù)向量，包含沉積層厚度、聲速、密度、衰減系數(shù)；xt-1為上一時刻的狀態(tài)，對應于反演問題中上一幀的地聲參數(shù)向量；yt為第t時刻的觀測值，對應于反演問題中的第t幀的反射系數(shù)觀測向量(包含N個掠射角)；f為狀態(tài)轉(zhuǎn)移函數(shù)，g為測量函數(shù)；g(xt)為以xt為地聲參數(shù)向量，代入反射系數(shù)計算公式，所求得的反射系數(shù)向量；nt和vt分別是狀態(tài)方程與觀測方程的噪聲，vt一般為高斯噪聲，則yt即為以g(xt)為均值的高斯分布。假設各掠射角的高斯噪聲獨立，似然函數(shù)L可以表示為

粒子濾波可能會由于似然變化較大而出現(xiàn)粒子枯竭的情況。為了應對粒子枯竭，可以采用重采樣的方式，重采樣是指維持總粒子數(shù)量不變的情況下，復制權重大的粒子，同時舍去權重非常小的粒子，是一種有效防止粒子枯竭的方法。

將粒子濾波應用到地聲參數(shù)反演中，在地聲參數(shù)隨距離緩慢變化時，傳統(tǒng)粒子濾波具有良好的效果。但當?shù)芈晠?shù)隨距離劇烈變化時，似然函數(shù)變化劇烈，即使使用了重采樣，傳統(tǒng)粒子濾波仍可能會出現(xiàn)粒子枯竭[11]的情況，無法正確收斂到目標分布。為了使粒子濾波可以在地聲參數(shù)劇烈變化的情況下收斂到目標分布，有文獻使用了帶橋接的粒子濾波的方法。使用橋接可以使粒子濾波結果緩慢收斂到目標分布，實現(xiàn)目標分布之間的平滑過渡。帶橋接的粒子濾波原理為[12]

其中，1＞α1＞α2…＞αM＞0，M為橋接層數(shù)。橋接的作用是使目標分布從逐漸收斂到，即為重要性采樣概率密度，當m=1時，αM+1=1，。當m=0時，，即為最終的目標分布。通過橋接，目標分布從q(x1:t)逐漸收斂到，即最終收斂到參數(shù)后驗概率密度。橋接原理如圖3所示，每一步橋接中都使用重采樣從q(x1:t)平緩過渡到，以防止粒子枯竭。

圖3 橋接原理圖Fig.3 Bridging priciple diagram

2.2 方差隨機、單多層層厚混合采樣的粒子濾波

粒子濾波解決地聲反演問題時，狀態(tài)方程式(5)可以表示為式(15)形式：

其中：nt一般選擇為高斯分布，對于地聲反演這種局部最優(yōu)多的優(yōu)化問題，選擇高斯分布可能會因隨機性的不足而導致反演結果易陷入局部最優(yōu)中。尤其是對于地聲參數(shù)隨距離劇烈變化的情況，即xt與xt-1相距甚遠的情況，這時候采用高斯分布噪聲，更可能由于采樣隨機性的不足而導致反演陷入局部最優(yōu)。容易聯(lián)想到增大方差來增加粒子隨機性的方式，但如此一來由于搜索精度不足，也容易使反演問題陷入局部最優(yōu)。有文獻使用在粒子濾波后增加馬爾可夫鏈蒙特卡洛采樣[9]的方法，利用馬爾科夫鏈蒙特卡洛方法可以在目標分布附近采樣，從而增加粒子多樣性。但是馬爾可夫鏈蒙特卡洛方法計算量巨大，代價高昂。

為解決此問題，在每一步中令nt的方差隨機，使用隨機方差采樣。令z=0+nt，用方差固定和方差高斯這兩種高斯采樣方式，分別進行100 000個點的采樣。方差固定采樣是在方差固定為1的零均值高斯分布中采樣，隨機方差采樣相對于固定方差采樣，其方差由固定的1變?yōu)殡S機。隨機的方差可以在方差為1的零均值高斯分布中采樣獲得。如圖4所示為兩種采樣時，變量z的樣本分布情況。

圖4 兩種采樣方式(方差固定和方差高斯)時，變量z的樣本分布Fig.4 Sample distributions of two sampling methods(fixed variance and Gaussian variance)

從圖4中可以看到，在方差為高斯分布的情況下，樣本在中心更集中，同時，其樣本又有較好的兩端擴展能力，即在距離中心較遠的地方仍然具有相當數(shù)量的樣本點。在距離相關的地聲參數(shù)反演問題中，每幀之間的地聲參數(shù)既可以緩慢變化，也可以劇烈變化，所以方差隨機的高斯采樣方式十分適合距離相關的地聲參數(shù)反演問題。本文中的方差隨機采樣方法，同樣達到了增加粒子多樣性的效果，進一步防止了粒子枯竭，有助于目標分布收斂。相較于固定方差的采樣，3.3節(jié)中的模擬結果驗證了方差隨機采樣可以使反演結果具有更好的穩(wěn)定性。

對于采樣的維度，傳統(tǒng)粒子濾波中一般使用多維高斯分布采樣，采樣維度與待反演地聲參數(shù)的維度相同。但是根據(jù)敏感度曲線，沉積層厚度有局部最優(yōu)，使用高維高斯分布采樣雖然保證了全局性能，但由于搜索精度不夠，不容易找到最優(yōu)解。如果在搜索中混合低維采樣，相當于提高了搜索精度，但同時全局性能不足，容易陷入局部最優(yōu)。本文提出了單層層厚搜索與多層層厚搜索混合的方式，層厚混合搜索的粒子濾波采樣時使用了沉積層厚度的單層多層混合改變的搜索方式，單層搜索指所有沉積層的層厚只改變?nèi)我庖粋€，多層搜索指所有沉積層的層厚全部改變，采用50%概率所有層層厚都變化，50%概率只改變?nèi)我粚訉雍竦姆绞?。此方法在增加了全局搜索能力的同時，增加了搜索結果的精度。3.2節(jié)的仿真結果驗證了此方法的有效性。

基于地聲參數(shù)反演的特點，本文使用了帶橋接、方差隨機、單層層厚與多層層厚混合搜索的粒子濾波方法，在提高搜索精度的同時保證了解的全局性，更好地處理了地聲參數(shù)隨環(huán)境劇烈變化的問題。

3 仿真結果分析

3.1 仿真環(huán)境及傳統(tǒng)粒子濾波反演結果

采用的仿真環(huán)境設置如下：地聲參數(shù)隨距離變化的100幀海底環(huán)境，在第10幀和第11幀之間沉積層分層數(shù)目改變，同時沉積層聲速和密度也有較大變化。在第20幀和第21幀之間以及第80幀和第81幀之間沉積層分層數(shù)目也發(fā)生改變。第21幀沉積層由4層變?yōu)?層，第81幀沉積層變?yōu)?層。第1幀共4層沉積層，地聲參數(shù)參見表1中。在所有幀中，沉積層總厚度設置為固定的10 m。粒子數(shù)量設置為10 000個。水聽器數(shù)量設置為13個，收集在36°～79°之間共13個掠射角的反射系數(shù)。發(fā)射聲源發(fā)射兩個單頻信號，頻率分別為1 000 Hz和500 Hz。實驗裝置如圖5所示。

圖5 實驗布設示意圖Fig.5 Layout diagram of experiments

模擬的真實地聲參數(shù)分布如圖6～8所示。利用傳統(tǒng)粒子濾波算法反演結果如圖9～11所示。圖9中，每幀聲速用10 000個粒子的均值表示。圖10中，每幀密度用10 000個粒子的均值表示。圖11中，每幀衰減系數(shù)用10 000個粒子的均值表示。

圖6 模擬真實環(huán)境聲速Fig.6 Simulated sound speed in real environment

圖7 模擬真實環(huán)境密度Fig.7 Simulated density in real environment

圖8 模擬真實沉積層衰減系數(shù)Fig.8 Simulated attenuation in real environment

圖9 傳統(tǒng)粒子濾波沉積層聲速c反演結果Fig.9 Acoustic velocity c inversion results of sedimentary layer obtained by inversion with traditional particle filtering

圖10 傳統(tǒng)粒子濾波沉積層密度ρ反演結果Fig.10 Density ρ of sediment layer obtained by inversion with traditional particle filtering

圖11 傳統(tǒng)粒子濾波沉積層衰減系數(shù)α反演結果Fig.11 Attenuation coefficient α of sediment layer obtained by inversion with traditional particle filtering

可見傳統(tǒng)粒子濾波在地聲參數(shù)劇烈變化時，無法正確估計地聲參數(shù)，獲得的地聲參數(shù)有較大誤差。

3.2 層數(shù)混合采樣

使用帶橋接的粒子濾波，一定程度上能夠處理地聲參數(shù)隨距離劇烈變化的反演問題。但是僅僅使用橋接并不夠穩(wěn)定，粒子濾波會出現(xiàn)無法收斂或收斂到局部最優(yōu)的情況?；诘芈暦囱莸奶攸c，本節(jié)將2.2節(jié)提出的混合層厚搜索的方式，與只有橋接但沒有混合層厚搜索的粒子濾波做比較。

由于仿真環(huán)境設置為第11幀變化最劇烈，第11幀將是粒子濾波最具挑戰(zhàn)性的1幀。接下來對比這三種采樣方式在第11幀時的性能差異。隨機生成30個不同的仿真環(huán)境分別用單層層厚改變、多層層厚均改變、層厚混合改變的采樣方式進行帶橋接的粒子濾波，記錄通過這三種粒子濾波后，第11幀的反射系數(shù)最小均方誤差(反射系數(shù)最小均方誤差即指通過多層橋接后得到的最小均方誤差)。為了使均方誤差不至于太小，這里的均方誤差沒有除以樣本數(shù)量(26)，而是直接用|R-R?|2表示，其中R為觀測反射系數(shù)，R?為用所有粒子均值求得的反射系數(shù)。各掠射角的反射系數(shù)的標準差設置為0.01。將這三種方式所得的30次實驗的最小均方誤差按照從大到小排列，得到如圖12所示曲線。

圖12 三種采樣方式得出的第11幀的最小均方誤差Fig.12 The minimum mean square errors of the 11th frame obtained by the three sampling methods

從圖12可以看出，層厚的混合改變采樣方式相較于層厚均改變的采樣方式有更小的均方誤差，前者更容易獲得全局最小的反演結果；相較于層厚單層改變的采樣方式反演效果也有小幅度提升。粒子濾波采用層厚混合采樣的方式時，第11幀地聲參數(shù)有更大的概率收斂到全局最優(yōu)。

3.3 方差隨機搜索

狀態(tài)方程中的噪聲一般設為方差較大的高斯擾動，設置此方差較大的原因是為了增加粒子濾波的搜索范圍，使地聲參數(shù)至少有從某一幀采樣到下一幀真值的可能性。為了對比隨機方差的粒子濾波與固定方差的粒子濾波對反演效果的影響，進行了17次實驗，分別使用隨機方差(均值為0，方差為1的高斯分布)采樣、固定方差0.6(狀態(tài)方程噪聲方差的0.6倍，下同)采樣、固定方差0.8采樣、固定方差1采樣、固定方差1.2采樣、固定方差1.4采樣。沉積層厚度采用層厚均改變的采樣方式。得到的最小均方誤差結果如圖13所示。

由圖13可以看出，隨機方差的采樣方式大大提高了第11幀地聲參數(shù)收斂到全局最優(yōu)的可能性。注意到方差過大(1.4)或過小(0.6)將會出現(xiàn)第11幀最小均方誤差大于1的情況，說明此時粒子濾波效果已十分不佳。同時還可以看出，所有固定方差的采樣方式，最小均方誤差曲線都不足17個點，這是由于某些實驗中，固定方差的采樣方式已經(jīng)不能使粒子濾波穩(wěn)定運行，即沒有搜索到一個最優(yōu)或局部最優(yōu)點，這將導致粒子濾波失效。而隨機方差采樣的方式則可大幅提高粒子濾波的穩(wěn)定性，幾乎每次粒子濾波都能正常收斂。

圖13 不同方差采樣得出的第11幀的最小均方誤差Fig.13 The minimum mean square errors of the 11th frame obtained by different variance sampling methods

為了減少運算量，可以令橋接數(shù)量可變，橋接數(shù)量最大為9層，令反射系數(shù)均方誤差基本收斂時停止橋接。實驗采用3.2中的沉積層層厚混合改變的采樣方式，統(tǒng)計所用橋接數(shù)量。圖14是方差固定和方差隨機時橋接數(shù)量隨幀數(shù)的變化，圖15是將這些橋接數(shù)量按照從大到小排列的結果。

圖14 方差固定和方差隨機時橋接數(shù)量隨幀數(shù)的變化Fig.14 Variation of the number of bridges with the number of frames when variance is fixed or random

圖15 當幀數(shù)增加到100時橋的數(shù)量按降序排列圖Fig.15 Diagram of the number of bridges arranged in descending order when the number of frames is increased to 100

從圖15可以看出，方差隨機時平均橋接數(shù)量小于方差為固定值時的平均橋接數(shù)量。隨機方差采樣平均橋接數(shù)量為4.65層，固定方差采樣平均橋接數(shù)量為5.43層，方差隨機時算法收斂速度快?？梢钥吹皆诘芈晠?shù)隨距離變化劇烈時需要的橋接數(shù)量較多，變化緩慢時需要的橋接數(shù)量較少，這是可以預見的。

3.4 改進粒子濾波反演結果

在粒子濾波中使用橋接的同時，同時使用狀態(tài)方程方差隨機以及沉積層層厚混合采樣的方式，得到如圖16～18的反演結果。圖16～18中，每一幀的聲速、粒子、密度均用粒子均值代替。

圖16 改進粒子濾波聲速c反演結果Fig.16 Acoustic velocity c of sediment layer obtained by inversion with improved particle filtering

圖17 改進粒子濾波密度ρ反演結果Fig.17 Density ρ of sediment layer obtained by inversion with improved particle filtering

圖18 改進粒子濾波聲衰減系數(shù)α反演結果Fig.18 Attenuation coefficient α of sediment layer obtained by inversion with improved particle filtering

由反演結果可見，改進粒子濾波在沉積層參數(shù)隨距離劇烈變化的情況下，仍可較準確地反演地聲參數(shù)。真實環(huán)境中，沉積層層數(shù)會發(fā)生改變。但在改進粒子濾波反演結果中，沉積層層數(shù)一直使用的都是4層?？梢钥吹皆?0～80幀，以及80～100幀中，會出現(xiàn)某些層極薄且和上層或下層聲速、密度相近的情況，此時反演結果中這些層將不含明顯的分層顯示，所以一直使用4層沉積層仍能取得較好的反演結果。改進粒子濾波在第11幀、第21幀、第81幀聲速和密度都有較好的反演結果。改進粒子濾波在第11幀時反射系數(shù)擬合效果如圖19所示，反射系數(shù)曲線擬合效果較好。

圖19 第11幀反射系數(shù)R擬合效果Fig.19 Fitting results of the 11th frame reflection coefficient R

3.5 反演結果的邊緣概率密度

為了獲得某個地聲參數(shù)的分布情況，分析了反演結果的邊緣概率密度。圖20～22為第11幀聲速、密度和衰減系數(shù)反演結果的邊緣概率密度圖。圖20～22中粉色虛線為真值。

從圖20～22中可見聲速密度和衰減系數(shù)均在真值附近，其中聲速的不確定性較小，密度和衰減系數(shù)的不確定性大，這是由于似然函數(shù)對聲速的敏感性高，對密度和衰減系數(shù)的敏感性低。不確定性分析對應于第1節(jié)敏感性的分析，證實了反演結果的可靠性。

圖20 聲速c邊緣概率密度圖Fig.20 Marginal probability density diagram of sound velocityc

圖21 密度ρ邊緣概率密度圖Fig.21 Marginal probability density diagram of density ρ

圖22 衰減系數(shù)α邊緣概率密度圖Fig.22 Marginal probability density diagram of attenuation coefficient α

4 結論

針對海底地聲參數(shù)隨距離變化的地聲參數(shù)反演問題，本文提出了一種改進的粒子濾波序貫估計方法。在帶橋接的粒子濾波基礎上，使用了沉積層厚度混合變化采樣的方式，增加全局搜索能力的同時使搜索結果更加精確。在狀態(tài)方程中使用噪聲方差隨機采樣，增加了粒子的多樣性，同時減少了所需的橋接系數(shù)數(shù)量。相較于傳統(tǒng)粒子濾波，該方法在地聲參數(shù)隨距離劇烈變化的情況下，依舊能較穩(wěn)定地反演地聲參數(shù)。相較于只有橋接的粒子濾波，本文方法的反演結果更加穩(wěn)定，計算效率更高。

利用不同掠射角和不同頻率的反射系數(shù)作為待擬合的觀測數(shù)據(jù)，采用改進粒子濾波的方法進行地聲參數(shù)反演。仿真結果證明了在沉積層參數(shù)劇烈變化的環(huán)境中，本文方法仍然具有良好的效果。