屈科 樸勝春 朱鳳芹
1)(哈爾濱工程大學,水聲技術重點實驗室,哈爾濱 150001)
2)(海洋信息獲取與安全工信部重點實驗室(哈爾濱工程大學),工業和信息化部,哈爾濱 150001)
3)(哈爾濱工程大學水聲工程學院,哈爾濱 150001)
4)(廣東海洋大學,廣東省近海海洋變化與災害預警重點實驗室,湛江 524088)
聲速剖面是海洋聲場分析必要的先驗信息,海洋內波、潮汐、鋒面等動力活動的聲學監測往往也是通過聲速剖面或者其擾動的反演來實施的[1-3].從理論上講,聲速剖面可以表示為時間或空間隨深度變化的矩陣形式,矩陣越大精度越高.然而,參數的增加會極大地增加逆問題的求解難度,必須通過一定的降維技術來構建聲速剖面才能保證其快捷有效地獲取.因此,聲速剖面構建方法的準確性和適用性,以及其能否有效反映海洋動力過程的物理特征,成為了海洋動力活動聲監測中的一個關鍵問題.
目前最廣泛使用的聲速剖面構建方法[4-6]是正交經驗函數(empirical orthogonal functions,EOF).通過若干組聲速剖面正交基以及對應的投影系數,極大地減少了描述聲速剖面垂直結構所需參數,可有效地應用于海洋聲層析[7-10]、聲速剖面時空變化分析[11,12]、海洋活動監測[13,14]、聲源定位[15-17]等領域.大量的理論研究及實際應用證實了EOF方法的有效性和可行性,同時也暴露了其在海洋環境聲監測領域存在著以下兩方面的不足:1)方法需要一定量的樣本進行特征提取,當樣本過少或者沒有完整覆蓋海洋活動周期時,難以有效地對聲速剖面進行構建;2)EOF實質上是樣本數據矩陣的特征,是“數據”的展開基函數,它并沒有明確的物理意義,從反演獲得的若干個投影系數上很難直觀地獲得海洋動力活動的信息.在實際海上應用中,最亟待監測動力活動的海域往往是缺乏現場測量數據的.缺乏樣本的情況下EOF方法獲取基函數較為困難,對反演結果的進一步分析處理也會影響監測的實時性.除了EOF方法,聲速剖面也可以通過一定的解析函數模型進行構建.解析函數模型一般通過一系列的數學表達式及參數描述聲速隨深度的變化關系,例如Munk模型[18]、連續介質力學的離散元方法模型[19]、分層聲速剖面模型[20]等.然而,解析模型通常有一定的適用區域,為保證精度參數往往較多,且涉及的參數也與海洋動力特征沒有直接聯系,因此很少應用于海洋動力活動的監測中.
本文提出一種全新的基于內潮動力特征的聲速剖面構建方法.根據流體動力學原理,結合聲速擾動與水質子運動的關系,利用內潮簡正模式建立聲速剖面的水動力模式基函數(hydrodynamic mode bases,HMB)進行聲速剖面構建.由于潮汐運動具有較強的周期性和季節性特征,HMB可以從WOA13氣候態數據[21]等歷史資料中直接計算獲得,實現了不借助實時現場測量樣本的聲速剖面展開基獲取.同時,這種基函數及其對應投影系數與海洋動力過程有直接的聯系,較之傳統的方法具有更明確的物理意義.結合2001年東中國海中美聯合實驗數據的聲速剖面重構證實了方法的有效性,對比EOF方法對HMB的準確性和物理意義進行了討論.最后,將HMB應用于匹配場聲層析,初步探討了這種聲速剖面基函數應用于海洋動力活動監測的可行性.與EOF等傳統方法相比,基于內潮動力特征的HMB可以僅依靠數據庫歷史數據獲得,其參數與海洋動力活動的物理特性直接相關,這些特點在海洋動力活動的聲學監測中是非常有意義的.
按照時間尺度的不同,海水聲速剖面的變化可以分為大尺度的背景變化以及小尺度動力活動引起的變化[22].大尺度的變化主要影響較長時間段內聲速剖面的穩態背景特性,例如溫躍層的深度與厚度等,它帶有明顯的季節特性;而小尺度的變化通常指短時間內的變化,例如線性和非線性內波等引起的聲速剖面擾動,它對穩態背景剖面沒有影響.以EOF為代表的傳統聲速剖面構建方法主要采用的是數據特征提取的原理,在構建某時刻聲速剖面時必須對臨近時間和空間內的聲速剖面進行特征提取,這就決定了其樣本依賴性.從海洋物理動力過程的角度出發,季節性的海洋背景層化特征是較為穩定的,而短時間的內波等垂向動力過程特性是受穩態的背景層化特性所控制的,所以已知穩態背景層化特性的情況下可以對某時刻聲速剖面的瞬態擾動特性進行推斷.與此同時,隨著WOA13氣候態數據等數據庫的建立,全球絕大多數海域的季節性層化背景特征是已知的,這就為無現場實時測量的情況下僅憑數據庫歷史數據來獲得聲速剖面展開基函提供了可能.
取直角坐標原點于靜止海平面,不考慮背景流場影響的水質子運動滿足方程[23]:

其中水質子的振速包含三個分量 (x,y,z),x軸沿著緯度向東,y軸沿著經度向北,z軸垂直于海面向上;W是水質子垂向運動速度w的振幅;σ是頻率;k和l對應為水平方向 (x,y)的波數且χ2=k2+l2;g是重力加速度;f是科氏力;N浮力頻率.引入海洋研究中常用的剛蓋近似等假定,方程可進一步簡化為

中尺度的海洋動力活動會引起聲速剖面的變化,其變化過程是滿足上述流體力學控制方程的.同時,海水質子受到海底和海面兩個邊界的制約,其垂向的運動具有簡正模式,質子的垂向運動可以表示為模式的疊加.考慮到內潮幾乎存在于所有淺海海域,且具有穩定的半日或全日周期的內潮波往往是造成等聲速線波動的主要因素[24].基于歷史數據的季節性層化信息,考慮內潮波的運動,引入長波近似并忽略地轉效應,本征函數ψn可以利用更為簡單的Sturm-Li-ouville方程計算:

其中cn是相速度.這里本征函數隨頻率的變化可以被忽略,質子的垂向運動可以表示為[25]

其中Wj是第j號簡正模式ψj幅值.
由于溫度對聲速的影響遠大于鹽度,且通常情況海洋的鹽度變化區間相對較小,下面的分析只考慮了溫度的影響(如果要考慮鹽度的影響,分析方法類似).
海水物態滿足熱力學公式:

其中ρ是相速度,T是溫度,cv是等體積比熱容,kr是熱導率,QT代表熱源.將水質子活動視作沒有熱源的絕熱過程,kr和QT均為零,水質子運動引起的溫度隨時間t的變化可以簡化為[26]

將(4)式代入(6)式得到:

與聲速的變化類似,溫度可以表示為季節性的穩態背景T0以及短時間內動力活動造成的擾動T′,那么 dT0/dz>>dT′/dz,對(7)式兩邊進行積分可以得到溫度剖面的表達式:

這里溫度就被表示成與EOF方法類似的展開形式:

展開基函數φn為背景剖面梯度與對應階數本征函數的乘積:

投影系數為

在淺海環境中,海水聲速與溫度有非常強的相關性,溫度每升高1 ℃聲速也隨之增加約4.0 m/s,因此也可以近似地用φn來作為聲速剖面的展開基.聲速剖面C可以展開為如下的形式:

實際應用中C0為季節性的穩態背景聲速剖面(也可以通過一個海洋動力活動的周期平均獲得),ηn為第n階投影系數.因為是從海水動力方程中推導獲得,這里將這種聲速剖面展開基稱作HMB.在文獻[25]中,基于數據庫歷史數據計算的水動力模式被直接用來構建聲速剖面.由于聲速不是一個動力學量,所以其有效性難以保證.文獻[26]解決了動力模式轉化為聲速擾動的問題,但是由于動力方程適用性的限制,仍然必須依靠現場實時測量的樣本數據才能求得基函數.HMB方法結合了兩種方法的特點,實現了依靠數據庫歷史數據獲得有效的聲速剖面展開基函數.
較之EOF聲速剖面基函數,HMB最大的特點在于它有明確的物理表達式,而不再是單純的樣本數據特征提取.從季節性的穩態分層特征可以計算浮力頻率N,這就決定了水質子運動的各階模式ψ.考慮溫度為聲速擾動主要因素的情況下,結合背景溫度剖面的梯度就可以計算水動力基函數φ.季節性的穩態背景剖面可以從WOA13等數據庫的多年統計平均或者同化數據中獲得,這就使不依靠現場實時的測量而獲得聲速剖面展開基函數成為了可能.此外,較之EOF方法,HMB有更明確的物理含義.一方面,季節性層化特性決定了水動力模式,HMB中包含了水質子運動的模態信息.另一方面,季節性層化特征對內波等海洋動力活動有控制作用,HMB對應的投影系數也與海洋動力活動存在著一定的內在聯系.下面將通過海上實測數據的處理,對HMB進行分析.
實驗數據來自2001年夏季中美兩國在東中國海進行的聯合科學考察實驗[27],本節將利用CTD(conductance-temperature-depth)測量數據進行聲速剖面重構分析,并基于聲傳播實驗數據進行匹配場層析.
從2001年6月2日19點35至6月6日23點23分,實驗3號在拋錨點(29°40.47′N,126°49.21′E,海深105 m)按照約1 h的間隔(輔機停電時,暫停測量)進行了52次CTD測量,圖1為換算所得的聲速剖面數據.可以看到聲速剖面隨時間的變化較大,20—70 m的深度范圍內有明顯的躍層,海深80 m到海底變化相對較小.實驗觀測中發現30 m和60 m深度附近可能存在冷水團,這可能是導致部分數據相應深度出現較低聲速的原因.

圖1 CTD測得的聲速剖面(細線)和平均聲速(粗線)Fig.1.Sound speed profiles (thin line)and average sound speed profile (rough line)measured by CTD.
為了分析不依靠現場實時樣本情況下獲得的聲速剖面展開基,這里從WOA13夏季氣候態數據中提取了實驗點溫度和鹽度剖面,并計算了海水的浮力頻率.WOA13是美國國家海洋大氣局發布的海洋氣候態數據集產品(https://www.nodc.noaa.gov/OC5/woa13),包含了全球海洋溫度、鹽度、密度、溶氧量、磷酸鹽等多種海洋要素.它是一種整合了多種數據集和實測數據的平均格點化數據,分為年平均數據、季節平均數據以及月平均數據,空間分辨率有5°,1°,0.25°三種.深度上,數據集采用內插法,其中0—100 m的分辨率為5 m[21].本文選用了多年(1955—2012)空間分辨率0.25°的季節平均數據,圖2為實驗點夏季的平均溫度、鹽度以及浮力頻率剖面.從海水的層結特性來看,浮力頻率在20—70 m的深度有較大的值,這與實測數據中溫度躍層位置對應,體現了數據的季節性穩態背景特征.
通過(10)式可以計算得到HMB,這里將其與EOF方法的聲速基函數進行比較,前三階的計算結果如圖3所示.可以看出,由于水動力聲速基源自海水層結的特征函數,隨著模態數的增加其拐點相應地遞增.其中第一階聲速展開基的極大值在58 m的深度,這與EOF方法獲得正交基有一些不同,但兩者的垂直分布上有一定的相似之處,這表明季節性穩態背景剖面雖然能體現一定的穩態特性,但是并不能非常精確地描述剖面在某一特定時刻的瞬態變化.其他階聲速剖面基函數的分布也體現了相同的特點,接下來將通過聲速剖面的重構分析對這種穩態特性獲得的聲速剖面展開基的有效性進行分析.

圖2 溫度、鹽度以及浮力頻率剖面Fig.2.Temperature,salinity and buoyancy frequency profile.

圖3 不同方法獲得的聲速剖面展開基比較(黑線,HMB方法;紅線,EOF方法)Fig.3.Comparison of the sound speed profile bases obtained by different methods (black line,HMB;red line,EOF).
為了檢驗聲速剖面構建效果,本節分別基于HMB和EOF方法對聲速剖面進行重構分析.考慮到CTD測量中不可能完全在等距離點處進行采樣,因此采用樣條插值法得到了從海面到海底采樣距離為1m的106個等距離點的聲速值,并基于這些采樣點與聲速剖面展開基函數進行了重構.為分析重構剖面的準確性,定義重構的均方根誤差為

其中ci中分別為測量和重構剖面對應深度的聲速值,N為深度上的采樣點個數.
表1為HMB及EOF方法的重構誤差分析.可以看出除了只用一階重構以外,兩種方法的重構誤差相差不多,EOF方法效果稍好于HMB方法.對比其他文獻對同組數據的重構[28],兩種方法使用前六階基函數均能較好地構建出聲速剖面.由于HMB方法在計算基函數時只利用了穩態的季節性分層特性,包括年際變化、海氣相互作用、水團等因素的影響勢必會導致HMB方法重構效果不及現場樣本的特征提取.但是,從重構結果來看,從數據庫歷史數據獲得的HMB仍然能保證較高的重構精度,基于穩態層結信息的聲速剖面構建是可行的.表2為不同方法中各階所占聲速剖面變化比重,可以看到兩種方法中前兩階基函數所占百分比遠大于其他階.結合表1和表2的數據,EOF方法由于提取了包含重構剖面在內的多個樣本特征,重構效果稍好于HMB方法.但是,HMB方法基于歷史數據獲得的基函數是可以獲得較好的重構效果的.

表1 不同方法重構效果的誤差分析(單位:m/s)Table 1.Error analysis of different reconstruction methods (m/s).

表2 不同方法各階所占聲速剖面變化比重Table 2.Proportion of sound speed variety for each order in different methods.
參考逆問題應用中的通常做法,這里選擇了兩種方法的前三階基函數重構結果進行分析.圖4給出了均方根誤差隨樣本的變化,總體來看重構效果與表1結果一致,絕大多數的重構均方誤差都在0.5 m/s以下.HMB和EOF重構均出現了一些誤差稍大的情況,下面將對HMB誤差最大(1.14 m/s)的第8號樣本進行分析.

圖4 不同樣本重構效果的誤差分析Fig.4.Error analysis of different samples.
圖5為第8號樣本兩種方法重構結果與真實聲速剖面的對比.由于湍流混合、水團等因素的作用,聲速剖面有時會在很小的深度范圍內出現劇烈的拐點或者大梯度變化.因為這種聲速的變化并不受季節性層結特征控制,用若干階的HMB重構難以表現這種細微的突變結構;與之對應,由于這種聲速剖面的突變具有隨機性,并非大量樣本具有的共有特征,所以若干階的EOF重構常常也難以避免這種誤差.圖5中EOF重構結果在一定程度上體現了這種聲速剖面的突變,主要是由于圖1所示的樣本數據中很多都在30 m和60 m附近出現了低溫水團的影響,因此這一特征被保留進了EOF基函數并反映在聲速剖面的重構上.值得一提的是,這種聲速的突變是具有隨機性的,所以圖4中某些樣本HMB方法的重構效果要好于EOF方法.總體來講,由于EOF方法是數據特征的提取,在突變的細結構描述上它可能是好于HMB方法的.但是,當突變結構并非樣本普遍特征時,包括HMB方法在內的幾乎所有聲速剖面重構方法都是難以描述.

圖5 第8號聲速剖面重構Fig.5.The 8th sound speed profile reconstruction.
從聲速剖面的重構來看,HMB方法可以保證較高的精度,但是無法避免在細結構上的一些誤差.由于聲場的計算是重構聲速剖面最核心的目標,重構剖面能夠準確計算聲場就成為了評價方法有效性的核心問題.本節采用重構誤差最大的第8號樣本的聲場計算進行分析.重構仍然采用前三階基函數,波導參數與實驗保持一致,聲源深度為50 m,頻率200 Hz,采用水平不變波導,海底密度1.86 g/cm3、聲速1610 m/s、衰減系數0.15 dB/λ,海深105 m.
圖6為利用所得聲速剖面計算的傳播損失.可知兩種方法計算的聲場與它們對應的CTD數據計算聲場基本一致.在50 km的范圍內,HMB方法計算非相干傳播損失的誤差不超過1.6 dB,EOF方法的誤差不超過0.8 dB.聲場計算結果表明雖然在重構過程中一些細結構存在誤差,但仍然能保證較高的聲場計算精度,兩種方法均能有效地表示聲速剖面結構.從聲學應用的角度來講,HMB方法在依靠數據庫歷史數據獲得基函數的情況下重構的聲速剖面仍然能保證較高的聲場計算精度.

圖6 聲速剖面計算傳播分析(單位:dB)(a)CTD測量;(b)HMB重構;(c)EOF重構Fig.6.Analysis of transmission loss calculated by sound speed profiles (dB):(a)CTD measurement;(b)HMB reconstruction;(c)EOF reconstruction.
根據(10)式的推導,可以看出HMB方法的聲速基函數與水質子的運動模態相關,聲速基函數包含了水質子運動各階模態的規律.在實際的聲學逆問題應用中,各階投影系數通常是反演獲得的直接結果,建立它們與海洋動力活動的直接聯系將極大地便利海洋動力活動的監測.
這里對前兩階垂直模態及其投影系數進行分析.圖7為前兩階模態的歸一化幅值,可以看到隨著序號的增加,幅值的拐點增加.第一階模態的幅值在整個水體中的符號是一致的,它代表了不同深度水體以不同振速的同向運動.因為模態間的正交關系,在第一階模態取最大值的深度第二階模態幅值為零.由于前兩階基函數在聲速剖面構建中所占權重非常高,如果忽略其他階的貢獻,可以利用第一階投影系數來近似描述特定深度等聲速線的振幅.圖8為重構聲速剖面中1525 m/s等聲速線深度與第一階投影系數隨樣本的變化.兩者的變化趨勢非常相近,兩條曲線的Pearson相關系數為0.96,說明第一階模態的系數與等聲速線的變化具有非常強的相關性.這里1525 m/s是第一階模態最大值50 m深度附近平均聲速值,仿真表明距離一階幅值最大深度越近,兩者的相關性越高.第二階模態具有兩個極值,兩個極值的符號相反且分別位于一階極值深度的上方和下方.它代表了等聲速線在不同深度相反方向的運動,這意味著躍層厚度的變化,同時也代表了躍層內聲速梯度的變化.據此原理,可以利用第二階投影系數來描述溫躍層內的聲速梯度變化.圖9為躍層(30—70 m)內聲速梯度與第二階投影系數隨樣本的變化.兩者的變化趨勢基本一致,兩條曲線的Pearson相關系數為0.99,說明第二階模態的系數可以用于推導溫躍層內聲速梯度的變化.基于HMB方法獲得的反演結果本身與海洋動力活動直接相關,結合模態的分布可以對一些基本海洋動力參數進行求解,這在海洋動力活動監測方面是非常實用的.

圖7 前兩階模態振幅Fig.7.Amplitude of the first two modes.

圖8 1525 m/s等聲速線深度和第一階投影系數隨樣本的變化Fig.8.Variations of the depth at a sound speed of 1525 m/s and the first-order projection coefficient with samples.

圖9 聲速梯度和第二階投影系數隨樣本的變化Fig.9.Variations of the sound speed gradient and the second-order projection coefficient with samples.
上面的討論中得到HMB方法可以有效構建聲速剖面的結論.然而,聲速剖面展開基函數最重要的應用場合是在水聲學逆問題的求解中.反演過程中存在的多種復雜因素與基函數的相互作用影響著反演的效果,因此有必要通過實際的反演對HMB方法在逆問題應用中的效果進行分析.
在前述東中國海實驗中也進行了聲傳播實驗,聲源為38 g定深50 m的爆炸聲源,接收采用32元垂直陣.水聽器覆蓋的深度為4.5—90.5 m,32.5 m以上水聽器間隔約2 m,32.5 m以下間隔約4 m,其中36.5 m及52.25 m的兩個水聽器損壞,所以反演使用了余下的30個水聽器的聲壓數據.傳播實驗過程中,海深可近似為恒定105 m.為檢驗HMB方法的普遍適用性,采用常規的匹配場聲層析方法.反演使用寬帶非相關Bartlett處理器,匹配處理使用頻段為99—201 Hz,共35個頻率點數,整體流程如圖10所示.使用遺傳算法在尋優空間中快速尋找一組參數向量m使得代價函數

達到最小值,其中L為頻率點數,N為水聽器陣元數,Pnel(m)及Pncl(m)為取參數向量m時第n號水聽器在第l個頻率點的實測聲壓和拷貝場聲壓,*號表示共軛.反演中采用三階聲速剖面基函數反演聲速剖面,投影系數的搜索區間均為—5—5,海底采用前述海底參數,選擇概率為0.5,變異概率為0.05,交叉概率為0.8,種群數為100,遺傳5000代,采用多組并行計算的辦法保證結果的收斂性并為后驗概率分析采集樣本.
圖11為聲源在不同距離時的反演結果,爆炸聲源分別為10.2 km (2個)及12 km處.從圖11可以看出,兩種方法反演的結果基本一致,在不同距離HMB方法獲得的聲速剖面也基本一致.

圖10 反演流程Fig.10.Inversion process.
為了考察結果的可靠性,對反演結果的后驗概率分布進行分析.將遺傳算法產生的采樣值按照對應代價函數的大小進行排序,在形成概率分布時通過Boltzmann函數進行加權,可以得到第k組的參量的概率分布為[29]

其中Nobs為樣本數量,T為類似模擬退火算法中的溫度控制參數,可以取優化過程中最佳目標函數的平均值.樣本向量中第i個參數取值κ的邊緣概率分布為


圖11 不同距離的反演結果 (a)10.2 km;(b)10.2 km;(c)12 kmFig.11.Inversion results at different ranges:(a)10.2 km;(b)10.2 km;(c)12 km.
其中δ為狄拉克函數.圖12給出了對10.2 km處(圖11(a))寬帶爆炸聲源信號進行匹配場聲層析所得前三階投影系數的歸一化一維邊緣概率密度分布.可以看出:1)3個參數均以最高概率收斂于全局最優值,反演結果可信度高;2)隨著參數遠離最優結果,概率密度迅速下降且基本無旁瓣,說明參數敏感和結果唯一性較好.根據反演的分析,基于HMB方法聲速逆問題求解結果有效且唯一性好,這有利于基于反演結果和聲速基函數對海洋動力參數進行直接的求解分析.

圖12 3個反演參數的邊緣概率密度分布 (紅線為代價函數最優值對應參數)Fig.12.Probability distribution for the three inversion parameters (the red line is the optimum value for the objective function).
利用HMB方法反演所得的聲速剖面對64.5 m深度的200 Hz單頻信號傳播損失進行了預報,預報值與水聽器真實測量值的對比如圖13所示.預報值與實測值符合得較好,除個別異常的點,預報曲線與實測值比較貼近.實驗中由于爆炸聲源源級和爆炸深度有一定的波動,這可能會導致預報的誤差.同時,由于走航路線海深是緩慢增加的,這導致了用水平不變模型計算的傳播損失預報值比實測值略微偏大的情況.但是總體來看,預報值能較好地反映傳播損失的變化趨勢.傳播損失預報結果和圖11中反演結果的一致性可以證實HMB方法在水聲學逆問題應用中的有效性,基于數據庫歷史數據獲得的聲速剖面展開基能有效反演出準確的聲速剖面.

圖13 傳播損失預報與觀測值Fig.13.Transmission loss predicted and measured.
針對海洋監測中對聲速剖面表示的降維需求,本文提出了一種基于水質子運動模態的聲速剖面展開基函數.該基函數僅需歷史數據中的季節平均剖面即可獲得.從聲速剖面的構建效果來看,HMB能夠較為準確地進行聲速剖面的構建.較之經典的EOF方法,HMB方法具有更明確的物理含義:基函數中包含了水體運動的各階模態特征,而前兩階投影系數可以有效描述等聲速線及聲速梯度的變化.最后利用匹配場聲層析實驗對HMB方法在逆問題求解中的有效性進行了分析,結果表明反演的投影系數有較好的敏感性,結果可信度高且唯一性好.
HMB方法對現場實時測量的要求相對低,且直接聯系海洋動力活動特性,適合用于現場實時樣本不太充分海區的動力活動監測.此外,結合反演所獲投影系數以及聲速剖面基函數,可以進一步推導海洋動力學參數用于水下動力活動監測以及其他一些用途,相關方法是下一步研究的重點.