張東玉，邢傳璽，吳耀文

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

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

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

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

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

1 基本原理

1.1 傳播模型

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

(1)

(2)

Kraken方法使用有限差分來對簡正波方程進行求解，波動方程的解可以寫為Φ(r)和Ψ(z)的乘積：

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

(3)

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

(4)

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

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

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

(11)

對上式進行傅里葉變換得到：

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

(12)

(13)

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

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

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

c(x)=c0+Qx,

(14)

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

pobs=g(x)+n,

(15)

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

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

(16)

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

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

(17)

(18)

(19)

(20)

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

2 方法實現(xiàn)

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

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

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

圖1 實驗環(huán)境圖 圖2 前3個EOF

圖3 SSP反演流程圖

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

1) 采用Kraken模型對水深 30 m 的海洋聲場進行了數(shù)值模擬，得到水聽器垂直陣列(VLA)中的復(fù)聲壓，并使用一階泰勒展開對其進行線性化.

pobs=g(x)+n,

(21)

(22)

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

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

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

(23)

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

(24)

(25)

3 結(jié)果分析

圖4 發(fā)射信號 圖5 水聽器接收的聲壓

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

圖6 8個水聽器的接收信號

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

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

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

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

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

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

圖10 波形預(yù)報對比圖

4 結(jié)語

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