張東玉，邢傳璽，吳耀文

(云南民族大學(xué) 電氣信息工程學(xué)院，云南 昆明 650500)

海洋擁有著豐富的資源，如礦物資源、生物資源、水資源，海洋資源的開(kāi)發(fā)、利用和保護(hù)是水聲學(xué)研究的最終目的，水下目標(biāo)探測(cè)、導(dǎo)航、定位以及聲吶系統(tǒng)性能的預(yù)測(cè)等研究成為近幾年的熱點(diǎn)，這些研究的開(kāi)展都必須首先獲得海洋環(huán)境參數(shù)[1].聲速作為水聲研究中重要的環(huán)境參數(shù)之一，是了解和研究水聲信道和水聲傳播的基礎(chǔ)，是所有水聲學(xué)研究領(lǐng)域都不可缺少的重要環(huán)境參數(shù)，海洋聲速剖面在聲場(chǎng)建模和水聲設(shè)備性能的評(píng)估方面占據(jù)著十分重要的地位[2].

采用直接測(cè)量的方法測(cè)量水聲環(huán)境中的聲速，需要多次實(shí)地測(cè)量，需要耗費(fèi)大量的物力人力，相較而言，反演的方法只需一次測(cè)量就可多次重復(fù)使用，能降低成本，節(jié)省資源[3].隨著水聲理論日漸成熟，水聲設(shè)備技術(shù)不斷完善，為了更好的研究、利用和保護(hù)海洋，大量學(xué)者提出利用己有的歷史聲速剖面，通過(guò)反演的方法獲得海洋的垂直聲速剖面，最常用的方法是使用經(jīng)驗(yàn)正交函數(shù)模型反演聲速剖面，但經(jīng)驗(yàn)正交函數(shù)模型反演聲速剖面的方法容易遺漏SSP中微小擾動(dòng)的特征點(diǎn)，同時(shí)計(jì)算效率低下，反演誤差較大[5].SSP反演問(wèn)題傳統(tǒng)上是通過(guò)最小化最小二乘代價(jià)函數(shù)的解的能量來(lái)正則化，需要對(duì)復(fù)雜的SSP結(jié)構(gòu)進(jìn)行低采樣或使用少量基函數(shù)來(lái)解釋SSP的變化[6].這種分辨率的降低會(huì)導(dǎo)致SSP的變化不確定，當(dāng)內(nèi)部波或水流產(chǎn)生強(qiáng)烈的、隨時(shí)間變化時(shí)，SSP的變化具有不確定性，這種不確定性會(huì)嚴(yán)重影響其他參數(shù)反演的準(zhǔn)確性.

近年來(lái)，有學(xué)者開(kāi)始將壓縮感知理論應(yīng)用于海洋聲學(xué)中，如應(yīng)用于波束形成，地聲參數(shù)反演.壓縮感知理論能夠?qū)⒖蓧嚎s信號(hào)在降采樣觀測(cè)后以很高的精度重建，打破了傳統(tǒng)采樣定理的局限.目前，壓縮感知理論已經(jīng)廣泛應(yīng)用于醫(yī)學(xué)、導(dǎo)航、遙感等領(lǐng)域.文獻(xiàn)[5]證明了CS在海洋聲學(xué)中的應(yīng)用，可以將重構(gòu)問(wèn)題可以轉(zhuǎn)化為求解最小l0范數(shù)問(wèn)題.在與距離無(wú)關(guān)的淺海環(huán)境中，SSP反演可以用CS來(lái)解決.在稀疏假設(shè)下，海洋中聲學(xué)SSP反演可以表述為一個(gè)不確定的線(xiàn)性問(wèn)題，在稀疏域中的SSP擾動(dòng)可以用形狀基函數(shù)來(lái)表示.

本文采用Kraken模型對(duì)淺海水聽(tīng)器的垂直線(xiàn)陣壓力觀測(cè)進(jìn)行了正演模擬.用一階泰勒展開(kāi)方法對(duì)正演模型的非線(xiàn)性響應(yīng)進(jìn)行線(xiàn)性化[7].在稀疏域中使用形狀基函數(shù)表示SSP的擾動(dòng).并用帶有l(wèi)1稀疏約束的最小二乘代價(jià)函數(shù)的凸優(yōu)化來(lái)求解字典系數(shù)向量，最后使用CS方法反演得到SSP.實(shí)驗(yàn)結(jié)果表明，使用CS方法可以很好地估計(jì)SSP，并驗(yàn)證了CS方法對(duì)于聲速剖面反演的優(yōu)點(diǎn).

1 基本原理

1.1 傳播模型

在反演海洋環(huán)境參數(shù)的研究中，選取適當(dāng)?shù)膫鞑ツＰ蛯?duì)水聲傳播實(shí)驗(yàn)中的聲場(chǎng)進(jìn)行建模是反演研究的重要前提，反演結(jié)果的精度由傳播模型的計(jì)算精度決定.從水聽(tīng)器的垂直線(xiàn)陣列(vertical line array，VLA)觀察到海洋環(huán)境中的聲壓，然后利用簡(jiǎn)正波傳播模型對(duì)海洋聲場(chǎng)進(jìn)行正向建模.簡(jiǎn)正波傳播模型的基礎(chǔ)是求解式(1)的波動(dòng)方程[3].

(1)

(2)

Kraken方法使用有限差分來(lái)對(duì)簡(jiǎn)正波方程進(jìn)行求解，波動(dòng)方程的解可以寫(xiě)為Φ(r)和Ψ(z)的乘積：

P(r,z)=Φ(r)Ψ(z).

(3)

將(3)式代入聲壓方程，可得：

(4)

(5)

把式(5)代入式(4)可得到深度方向的模式方程如式(6)：

(6)

假設(shè)海面是絕對(duì)軟邊界，即，ψ(0)=0，海底為剛性邊界，水深為D，滿(mǎn)足條件：

(7)

式(7)的解有無(wú)窮多個(gè)，每一階簡(jiǎn)正波都有對(duì)應(yīng)的模態(tài)函數(shù)Ψm(z)和相應(yīng)的一個(gè)水平波數(shù)krm.Ψm(z)也可叫做為本征函數(shù)，krm也可稱(chēng)做本征值.各階簡(jiǎn)正波的波數(shù)都是正交的，因此有：

(8)

聲壓函數(shù)可由所有簡(jiǎn)正波疊加表示：

(9)

利用貝賽爾函數(shù)推導(dǎo)出簡(jiǎn)正波的時(shí)域聲壓表達(dá)式為:

(10)

水聲傳播信道用沖擊響應(yīng)h(t)來(lái)描述，信號(hào)x(t)=e-iω0t作為激勵(lì)信號(hào)，用來(lái)激勵(lì)信道h(t)，則相應(yīng)的輸出為：

y(t)=x(t)*h(t)=p(r,z)e-iω0t.

(11)

對(duì)上式進(jìn)行傅里葉變換得到：

Y(ω)=X(ω)*H(ω)=2πδ(ω+ω0)p(r,z).

(12)

(13)

即求得信道沖擊響應(yīng)也即聲壓值P.

1.2 參數(shù)稀疏化和壓縮感知

由于CS是一種尋找非定值線(xiàn)性系統(tǒng)稀疏解的技術(shù)，因此利用一階泰勒展開(kāi)對(duì)SSP的非線(xiàn)性響應(yīng)進(jìn)行線(xiàn)性化[7].對(duì)于k點(diǎn)離散化海洋SSP，c(x)∈RK建模為：

c(x)=c0+Qx,

(14)

其中，c0∈RK是使用CTD探測(cè)儀測(cè)量經(jīng)轉(zhuǎn)換得出的離散參考聲速.Q=[q1,…,qn]∈RK×N是形狀基函數(shù)qn的字典，x∈RN是字典系數(shù)向量.SSP由形狀基函數(shù)建模，形狀基函數(shù)描述了在x中具有少數(shù)的非零系數(shù)的SSP擾動(dòng)[9].在M元VLA處接收的聲壓為pobs.

pobs=g(x)+n,

(15)

其中，pobs∈CM，g(x)∈CM是模型傳播函數(shù)，n∈CM是高斯白噪聲，聲源和其他海洋環(huán)境已知.

模型傳播函數(shù)的一階泰勒展開(kāi)的線(xiàn)性化響應(yīng)為[8]：

(16)

矩陣D=[d1,…,dn]∈CK×N包含了M個(gè)壓力觀測(cè)值相對(duì)于字典Q中的N個(gè)形狀基函數(shù)的導(dǎo)數(shù).dn通過(guò)對(duì)參考剖面進(jìn)行qn的一部分?jǐn)_動(dòng),用雙邊有限差分計(jì)算.

假設(shè)dn充分非相干，利用l1范數(shù)凸優(yōu)化方法找到x的稀疏估計(jì).式(15)的稀疏解可以表示為

(17)

(18)

(19)

(20)

其中，λ控制解能量和度量擬合的相對(duì)重要性.

2 方法實(shí)現(xiàn)

使用Kraken模型對(duì)簡(jiǎn)正波聲場(chǎng)進(jìn)行仿真計(jì)算，首先編寫(xiě)程序提前設(shè)置好式(10)中的密度、聲源距離、聲源深度、接收陣元深度、衰減系數(shù)等參數(shù)，調(diào)用軟件計(jì)算聲場(chǎng)可以得到模態(tài)函數(shù)和水平波數(shù)，再把得到的相關(guān)值和參數(shù)代入求聲壓的式中即可得到某一頻率對(duì)應(yīng)的聲壓值.

采用Kraken模型對(duì)水深為 30 m 的海洋聲場(chǎng)進(jìn)行了數(shù)值模擬.聲場(chǎng)由 100 Hz 的聲源產(chǎn)生，且聲源位于水下深度 10 m 處，在 1 km 范圍內(nèi)采樣，采用M列8個(gè)均勻間隔的VLA單元，跨度6～20 m.實(shí)驗(yàn)參數(shù)如下：海底聲速 1 700 m/s，密度 1.5 g/cm3，衰減 0.1 dB/λ.環(huán)境圖如圖1所示.

使用2013年在中國(guó)某海域?qū)崪y(cè)的聲速剖面數(shù)據(jù)作為SSP訓(xùn)練數(shù)據(jù)集，通過(guò)對(duì)SSP訓(xùn)練數(shù)據(jù)集進(jìn)行主成分分析，計(jì)算出一組EOF，前3個(gè)EOF如圖2所示.通過(guò)傳播模型計(jì)算獲得聲壓場(chǎng)的聲壓值，使用CS方法結(jié)合壓縮估算得到的稀疏系數(shù)x對(duì)SSP進(jìn)行反演.

圖1 實(shí)驗(yàn)環(huán)境圖 圖2 前3個(gè)EOF

圖3 SSP反演流程圖

聲速剖面反演步驟的過(guò)程流程圖如圖3所示，主要用到的聲速剖面反演具體算法步驟如下：

1) 采用Kraken模型對(duì)水深 30 m 的海洋聲場(chǎng)進(jìn)行了數(shù)值模擬，得到水聽(tīng)器垂直陣列(VLA)中的復(fù)聲壓，并使用一階泰勒展開(kāi)對(duì)其進(jìn)行線(xiàn)性化.

pobs=g(x)+n,

(21)

(22)

2) 從SSP訓(xùn)練集中使用數(shù)據(jù)壓縮算法主成分分析(principal component analysis,PCA)方法計(jì)算得到訓(xùn)練集的字典系數(shù)向量X和特征向量V.

3) 離散的K個(gè)點(diǎn)從CTD轉(zhuǎn)換得出的線(xiàn)性壓力值離散的K個(gè)點(diǎn)的參考聲速結(jié)合傳播模型得N個(gè)離散為K個(gè)點(diǎn)的形狀基函數(shù)qn的字典Q和M個(gè)壓力觀測(cè)值相對(duì)于Q中的N個(gè)形狀基函數(shù)的字典D.

Q=[q1,…,qn]∈RK×N,

(23)

D=[d1,...,dn]∈CK×N.

(24)

(25)

3 結(jié)果分析

圖4 發(fā)射信號(hào) 圖5 水聽(tīng)器接收的聲壓

發(fā)射信號(hào)為頻率 120 Hz，持續(xù)時(shí)間 2 s 的加Blackman窗進(jìn)行幅度調(diào)制后的線(xiàn)性調(diào)頻脈沖，時(shí)域信號(hào)及頻譜如圖4所示.圖5是水聽(tīng)器陣列接收到的聲壓值P，圖6是8個(gè)陣元接收到的信號(hào)，陣元接收到的信號(hào)來(lái)自同一信號(hào)，大致輪廓和信號(hào)形式是統(tǒng)一的，但由于陣元擺放位置不同，信號(hào)傳播的路徑不同，最終接收到的信號(hào)有誤差.由接收信號(hào)圖可以看出接收信號(hào)的時(shí)間約為 6.8 s，結(jié)合水下的平均聲速 1 500 m/s，可計(jì)算得到聲源水平距離在 10 000 m 附近，與 10 000 m 的聲源位置一致，計(jì)算結(jié)果證明，仿真得到的水聽(tīng)器陣列接收到的信號(hào)的聲壓值可進(jìn)行下一步的反演實(shí)驗(yàn).

圖6 8個(gè)水聽(tīng)器的接收信號(hào)

圖7 聲速剖面反演結(jié)果圖

為了驗(yàn)證實(shí)驗(yàn)結(jié)果的可靠性和有效性，進(jìn)一步驗(yàn)證CS方法對(duì)于聲速剖面反演的優(yōu)點(diǎn)，分析CS的性能.對(duì)圖7中的實(shí)測(cè)聲速剖面進(jìn)行30次聲速剖面計(jì)算反演過(guò)程，對(duì)最優(yōu)估計(jì)值統(tǒng)計(jì)分析，求解獲得反演聲速剖面的置信區(qū)間，置信區(qū)間定義為均值±標(biāo)準(zhǔn)差的結(jié)果，近似代表置信度為70%[10-11].反演得到的聲速剖面均值與測(cè)量聲速剖面、置信區(qū)間與測(cè)量聲速剖面對(duì)比如圖8所示.

從圖8看出，反演得到的聲速剖面均值與實(shí)際測(cè)量的聲速剖面非常接近，并且實(shí)際測(cè)量的反演聲速剖面置信區(qū)間內(nèi).在30組反演計(jì)算中，得到聲速剖面最大誤差為 0.642 8 m/s，最小誤差為 0.442 2 m/s，平均誤差 0.509 5 m/s.

圖8 反演結(jié)果對(duì)比分析

使用2013年在黃海實(shí)驗(yàn)的真實(shí)海試數(shù)據(jù)作為待反演聲速，分別使用傳統(tǒng)的EOF方法和本文方法進(jìn)行聲速剖面反演實(shí)驗(yàn)，反演結(jié)果圖如圖9所示，(a)圖是使用EOF方法選前3階經(jīng)驗(yàn)正交函數(shù)作為基函數(shù)得到的反演剖面.(b)圖使用EOF結(jié)合壓縮感知方法得到的反演剖面.從圖給出的2種方法的仿真結(jié)果來(lái)看，壓縮感知方法的重構(gòu)效果明顯優(yōu)于傳統(tǒng)的EOF方法.使用壓縮感知方法反演得到的聲速剖面和實(shí)測(cè)的聲速剖面，對(duì) 30 m 深度的 120 Hz 的單頻信號(hào)的傳播進(jìn)行波形預(yù)報(bào)，波形預(yù)報(bào)對(duì)比圖如圖10所示，反演聲速剖面得到的預(yù)報(bào)波形與實(shí)測(cè)聲速剖面得到的預(yù)報(bào)波形基本相同，研究結(jié)果表明使用文中方法對(duì)聲速進(jìn)行反演是一種行之有效的方法，在水聲學(xué)研究領(lǐng)域具有重要作用.

圖9 反演結(jié)果對(duì)比分析

圖10 波形預(yù)報(bào)對(duì)比圖

4 結(jié)語(yǔ)

本文使用一階泰勒展開(kāi)對(duì)SSP的非線(xiàn)性響應(yīng)進(jìn)行線(xiàn)性化，使用海洋聲速統(tǒng)計(jì)數(shù)據(jù)的先驗(yàn)知識(shí)以稀疏表示SSP的形狀基函數(shù)字典，將PCA方法應(yīng)用于海洋SSP數(shù)據(jù)以計(jì)算EOF，結(jié)合海洋聲傳播模型通過(guò)適當(dāng)更新的EOF訓(xùn)練字典，使用l1范數(shù)凸優(yōu)化方法找到稀疏估計(jì)值，結(jié)合訓(xùn)練字典使用CS算法對(duì)聲速剖面進(jìn)行反演.計(jì)算了SSP估計(jì)誤差，誤差的70%置信區(qū)間在 ±0.6 m/s 之內(nèi).實(shí)驗(yàn)結(jié)果表明，本文算法對(duì)SSP的反演是一種行之有效的方法.