基于多層水平變化淺海海底模型的低頻反演方法*

康娟 彭朝暉 何利 李晟昊 于小濤

1) (中國科學院聲學研究所,聲場聲信息國家重點實驗室,北京 100190)

2) (中國科學院大學物理學院,北京 101407)

準確獲取海底聲學參數對聲場分析、聲納應用等具有重要意義.反演是獲取海底聲學參數的重要手段之一.為了避免反演中的多值問題,地聲反演往往盡量減少待反演參數,多采用單層或雙層水平不變地聲模型,反演后的聲學參數僅能在有限頻段應用,難以同時應用于較寬頻段.本文以實驗中淺地層剖面儀測量的海底沉積層分層結構為基礎建立地聲模型,采用多物理量分步聯合手段反演各沉積層聲學參數,并通過匹配場定位、后驗概率分析、采樣等手段對反演結果進行驗證.將反演結果與單層等效地聲模型反演結果進行對比分析,結果說明,本文建立的多層水平變化海底聲學參數更接近真實情況,可以應用在較寬頻段范圍,解釋不同聲學現象,且精度更高.

1 引言

在淺海,海底作為海洋信道的下邊界,其地聲參數會對聲波傳播產生顯著影響.反演是一種有效獲取淺海地聲參數的方法,可以快速和低成本地估計出一定范圍內的海底沉積物聲速、密度和衰減系數等聲學參數,受到水聲學界的廣泛關注.國內外已經發展了多種反演方法,包括匹配場反演法[1,2]、海底反射系數反演法[3]、簡正波到達結構反演法[4]、傳播損失反演法[5,6]、脈沖時間到達結構反演法[7]、聲場頻率-距離干涉結構反演法[8]等等.采用單個物理量構成代價函數進行多個物理量反演時,往往要面對反演結果的多值性問題,因此有學者利用不同代價函數對海底參數敏感度不同,提出了多個代價函數分步反演地聲參數的方法.2004 年,Li 和Zhang[9]提出的多種反演方法聯合反演方案,經多個海區進行實驗驗證,一定程度上解決了反演中的多值性問題.2019 年,李夢竹等[10]先利用海底反射損失得到層厚以減少待反演參數維度,再利用匹配場反演其他參數,有效地降低了反演結果的不確定性.

在對淺海海底參數反演時,一般將海底等效為如圖1 所示三種水平不變海底模型[11].對于較窄頻率范圍的聲學應用,采用這三種模型通?？梢赃_到比較好的效果.但是實際應用時,往往需要更接近實際的多層地聲模型.崔智強等[12]利用傳播損失作為代價函數,仿真表明在假設三層沉積層的海底環境下,能夠提高地聲參數的反演精度.鹿力成等[13]通過研究淺海海底模型對于聲傳播的影響,發現較高頻率聲場主要受海底表層影響,而影響低頻聲場的海底深度能夠達到幾十米.此時,要想同時滿足較寬頻段的聲傳播,需要采用多層海底模型.

圖1 三種較為常見海底模型Fig.1.Three typical bottom models.

實際地聲反演中,除了要考慮海底分層結構外,往往還需要考慮海底水平變化情況.祝捍皓等[14]針對淺海水平變化波導下的低頻聲能量傳輸問題,基于有限元方法,討論了三種不同海底變化下帶來的聲能量差異.邱海賓和楊坤德[15]在利用拖線陣對海底參數進行反演時,建立了水平變化海底模型,在聲源頻率增多、拖線陣變長的情況下,依舊得到了收斂性好的匹配場結果.Kim等[16]考慮到海底地形在空間上的差異,將反演區域分為幾段進行反演,反演結果很好地解釋了實驗場地隨水平變化的聲學現象.Bonnel等[17]在利用啁啾聲反射手段得到多層結構的基礎上,提出一種在模態號數未知情況下利用頻散曲線反演隨距離變化多層海底參數的方法,成功得到沉積層的地聲參數.

2022 年8 月,中國科學院聲學研究所聲場聲信息實驗室在我國南海北部某大陸架海區進行了一次海底聲學反演實驗.實驗中利用淺地層剖面儀對海底進行了測量,如圖2 所示,可見海底具有明顯的分層結構及水平變化.本文利用本次實驗測量的淺地層剖面儀、海洋環境及聲學數據,開展水平變化分層海底地聲參數反演方法研究,建立起更接近真實海底并適用于較寬頻帶(50—500 Hz)聲場計算的多層水平變化海底地聲模型.

圖2 淺地層剖面儀獲取的海底分層結構Fig.2.Seabed layering obtained from sub-bottom profiler.

2 海上實驗及海底地聲模型

海上實驗選用沒有多次氣泡脈動干擾的氣槍作為聲源,其標稱深度為7.5 m.實驗期間獲取的氣槍信號聲源譜級曲線如圖3 所示.使用潛標垂直接收陣接收記錄聲信號,垂直陣共有15 個陣元,深度分布在5—80 m 之間.本文反演過程中將選用10 km 范圍內的氣槍信號用于計算海底聲速、密度和厚度,20 km 范圍內的氣槍信號用于計算衰減系數.

圖3 氣槍信號聲源譜級Fig.3.The frequency response curve of source level of air gun.

垂直陣不同深度處共綁扎了16 個自容式溫度壓力計(TD)組成溫度鏈,實驗期間獲取的聲速剖面如圖4 所示.可見是典型的淺海夏季溫躍層聲速剖面,存在線性內波,但聲速剖面起伏不大.因此,本文不考慮水文環境的水平變化,聲速剖面取信號發射期間的平均聲速剖面,如圖4 中黑線所示.

圖4 聲速剖面Fig.4.Sound-speed profile.

實驗期間通過淺地層剖面儀獲取了海底沉積層結構,如圖2 所示.對該海底分層結構進行簡化,得到地聲模型如圖5 所示.圖5 中,接收點位處為0 km,傳播路徑上的海底主要由3 層沉積層和1 層基底層組成.假設各沉積層內的海底聲速、密度及衰減系數不變,則待確定的參數主要包括各層的厚度、海底聲速、密度及衰減系數,共15 個參數.

反演前,可根據先驗信息減少待反演參數數量.首先,利用淺地層剖面儀得到的各沉積層聲波傳播時間,將各沉積層的厚度變換為海底聲速的函數:

式中,tm+1-tm為m層兩個邊界的往返時間差;cbm為m層的海底聲速.

其次,考慮到海底密度敏感性遠小于海底聲速,海底密度可以通過Hamilton 經驗公式由海底聲速表示[18]:

第三,利用較近距離聲信號反演海底聲學參數,一般代價函數對海底衰減系數不敏感.反演時,可以先將海底衰減系數設定為定值,先對圖5 中各層聲速cb1—cb4進行反演,然后將反演后的各層聲速作為已知條件,利用不同距離和頻率的聲傳播損失來反演海底衰減系數.

反演中,需要選擇合適的聲場計算模型.在對淺海低頻聲傳播進行研究時通常采用簡正波理論,由于所考慮的海底水平變化較弱,且將各層海底參數視為隨距離不變,因此在計算聲場時不考慮簡正波模態之間的耦合,采用絕熱簡正波進行表示[19]:

其中krm為簡正波特征值;φm為特征函數.

3 反演方法及過程

3.1 反演流程

基于上述分析,對分層海底地聲模型的聲學參數采取分步反演的方法,反演流程如圖6 所示.首先,利用淺地層剖面儀獲得海底分層結構,選取簡正波頻散結構對各層聲速進行反演.再利用匹配場反演方法對表層聲速進一步反演.最后將利用實測的傳播損失反演海底衰減系數.在反演過程中采用多種方式對反演結果進行驗證.

圖6 反演流程圖Fig.6.Flowchart of inversion.

為避免待反演參數陷入局部最優,采用并行多種群遺傳算法[20]進行優化搜索.具體操作為:編碼采取20 位二進制數;種群數目選擇為3;每個種群有24 個個體;多個種群交叉概率在0.7—0.9 間隨機產生;多個種群變異概率在0.001—0.05 間隨機產生;精英種群代價函數連續20 代不再變化時停止搜索.

3.2 簡正波頻散結構反演各層聲速

在淺海,接收的聲壓信號p可以用一系列簡正波的疊加來表示[21]:

式中,An(t),Φn(t) 分別代表在t時刻、第n號簡正波的幅度和相位.

利用簡正波頻散結構的反演海底聲速,其基本原理是:由于簡正波存在頻散現象,不同階數不同頻率的簡正波群速度不同,對一定距離外接收的脈沖聲信號進行時頻分析,提取出來頻散曲線與簡正波群速度密切相關,可用于海底聲速反演,采用如下代價函數[22]:

其中M為待反演的參數向量; ΔTn(f) 為第n號簡正波不同頻率信號到達時間差; ΔTnm(f) 為頻率f下第n號簡正波和第m號簡正波信號到達時間差.上標 e 代表實驗提取數據,c 代表拷貝場數據,這里通過簡正波模型KEAKENC[23]計算拷貝場.

在距離為r的情況下,不同頻率同號簡正波到達時間差 ΔTn(f)為

同一頻率不同號簡正波信號到達時間差ΔTnm(f)為

其中水平變化海底群速度近似為

利用簡正波頻散結構進行反演,關鍵在于如何精確地提取出各號簡正波的到達時間差(f)及(f) .本文采用warping 變換結合Wigner-Ville 分布的方法提取簡正波頻散結構并計算各號簡正波到達時間差.

Warping 變換是一種可逆的變換,可以通過逆變換算子對warping 變換后的信號采取時頻域濾波再恢復的方式分離出各號簡正波.分離得到的各號簡正波可以利用Wigner-Ville分布做時頻分析,提取出其頻散結構,Wigner-Ville分布的表達式為[25]

式中,*代表復數共軛,它是一種對信號瞬時自相關函數的傅里葉變換.

在實際反演過程中,首先選取距接收點右側約6 km 處的氣槍信號如圖7(a),提取接收深度處于75.93 m 的陣元接收信號的頻散曲線.如圖7(b)和圖7(c)為warping 變換后的頻域分布圖和時頻分布圖,其中第一階簡正波受噪聲干擾程度較大且與第二階簡正波形成干涉現象,故選取三到八階簡正波頻散結構用于反演,圖7(d)為提取出的三到八階簡正波頻散曲線.

圖7 采用warping 變換提取頻散曲線 (a) 接收信號;(b) warping 變化后的頻域信號(圖中紅色虛線為各號簡正波截止頻率);(c) warping 變換后信號時頻圖;(d) 接收信號時頻圖(圖中白色虛線為提取出的頻散曲線)Fig.7.Estimating dispersion curves using warping:(a) Received signal;(b) spectrum of the warped signal,and the cutoff frequency;(c) spectrogram of the warped signal;(d) spectrogram of the signal after warping and Wigner-Ville,and estimated dispersion curves.

采用后驗概率分析方法可以評估反演方法的有效性[26].在遺傳算法反演進化過程中,保存所有的搜索值Mk以及對應的目標函數 cost(Mk),對反演結果在搜索空間Ω內進行概率統計,將向量M按照其代價函數的大小進行排列,按照Boltzmann 分布進行加權,得到反演結果的第k組向量后驗概率分布:

其中Nobs是保存下來的模式向量的個數;T是溫度控制參數,一般取反演過程中的50 個最佳代價函數的均值較為合適;向量M中第n個參量值為m的邊緣概率密度為

針對頻散曲線圖7(d)進行反演時,保留樣本的邊緣概率密度如圖8 所示,其收斂程度依次為cb2,cb3,r,cb1,cb4.其中cb4的收斂 性最 差,且收斂于最高值2000 m/s 處,說明它對于代價函數敏感性較弱,后續將cb4固定為2000 m/s,不再對其進行反演.而相比于收斂性較好、收斂結果近似為單值的cb2和cb3,cb1存在一定的離散性.

圖8 反演參數的邊緣概率密度Fig.8.The marginal probability distribution of inversion parameters.

為了進一步檢驗反演方法的有效性,將反演得到的各海底參數作為輸入,計算距離約為6.2 km、深度處于75.93 m 的陣元接收到的信號對應的簡正波頻散結構,并與實驗數據提取的頻散曲線進行比較,如圖9 所示,紅線為模型計算的頻散曲線,藍線為實驗提取的頻散曲線,可以看出兩者符合較好.再將反演的海底參數用于匹配場定位,結果如圖10 所示,圖10 中用紅色圓圈標出匹配場定位的距離為5.9 km、模糊深度中值為7.75 m,與實際氣槍聲源的距離5.98 km、標定深度7.5 m 相比,誤差較小.

圖9 頻散曲線擬合圖Fig.9.Dispersion curves as estimated from the data with warping (blue) and predicted after inversion (red).

圖10 定位結果Fig.10.Localization result of Bartlett processor.

進一步對不同距離多個氣槍信號進行反演,得到的結果如表1 所列.統計反演結果發現,cb1產生了40 m/s 范圍的浮動,cb2產生了12 m/s 范圍的浮動,cb3產生了33 m/s 范圍的浮動,結合前面的收斂性分析,也表明了利用簡正波頻散結構反演二三層聲速較第一層更為準確,因此需要利用其他方式完成對于第一層聲速的進一步探究反演.

表1 不同距離下海底參數反演結果Table 1.Inversion results at different distances.

3.3 匹配場反演表層聲速

以寬帶相關Bartlett 匹配場處理器作為代價函數:

M為待反演參數向量;N為接收通道個數;pa(fi)為第a個垂直陣元接收到的復聲壓,e 代表實際數據; c 代表拷貝場數據; * 代表復數共軛.該代價函數φ位于0—1 之間,越靠近零越接近于真實場.

匹配場反演之前,考慮到不同頻率聲波在傳播過程中作用于不同海底深度,首先通過仿真分析匹配場在不同頻率時對海底參數的敏感性,以獲得更適用于反演表層聲速cb1的方案.

假設真實的海底結構如圖5 所示,共分四層,各層聲速cb1—cb4分別為1580,1630,1700,1800 m/s.與實驗情況相同,聲源深度為7.5 m,接收器位于5.9 km 處,深度位于5—80 m,共15 個水聽器,仿真計算垂直陣接收的信號用于模擬實驗獲取的聲信號,計算頻率為50—500 Hz.

首先,采用半無限大的單層海底,在不同頻段范圍內(50 Hz 帶寬)進行等效聲速cbe的反演,聲速搜索范圍為1500—2000 m/s,步長為5 m/s,結果如表2 所列.可以發現,50 Hz 頻率下反演得到的cbe為1725 m/s,位于真實海底結構中下層聲速cb3和cb4之間,隨著頻率增加,cbe逐漸減小,當頻率大于250 Hz 時,cbe接近于真實海底的表層聲速1580 m/s,說明較高頻率下聲波主要受表層作用,此時匹配場反演對底層聲速較不敏感,因此在利用匹配場進行多層海底地聲反演時,應選取較低頻率進行.而在較低頻率上下層海底共同產生影響時,需要對各層聲速進行敏感性分析,分析時選取100—150 Hz 頻段.

表2 仿真中不同頻率下反演獲取的單層等效聲速Table 2.Single-layer equivalent sound speeds at different frequencies.

下面,利用100—150 Hz 較低頻段對各層聲速進行敏感性分析,通過改變其中一個聲速,保持其他聲速不變的情況下,來測試φ(M) 的敏感性,結果如圖11 所示,圓圈代表其真實值.由圖11 可以看出,待反演的參數中第一層聲速cb1敏感性遠遠大于cb2和cb3.因此,可以利用高低頻匹配場結合的方式對表層聲速進行反演.

圖11 各層聲速敏感性分析Fig.11.One-dimension cost function value of sound speeds in each layer.

接下來,再利用實驗數據進行分析.首先將海底視為一個半無限大的單層海底,利用匹配場方法反演不同頻率下的單層等效聲速cbe,結果如圖12所示.圖12 中,紅點代表各頻段中代價函數最小值對應的cbe.圖12 中可以看出,與上面的仿真結果類似,頻率低于250 Hz 時反演得到的cbe應該受較深沉積層聲速影響較大,而頻率大于250 Hz 時反演的cbe應該只能反映表層聲速cb1的變化情況.盡管圖12 中不同距離下得到的cbe隨頻率的變化關系存在一定的差異,但仍然可以總結出如下規律:250—300 Hz的cbe在1600—1630 m/s之間,與3.2 節利用簡正波頻散結構反演得到的cb1較為接近;在350 Hz處cbe存在一個明顯的聲速躍變,此頻率下的cbe在1560—1580 m/s 之間;大于350 Hz后隨著頻率增加,cbe逐漸減小.因此可以根據隨著頻率增加,聲線可穿透的海底厚度減小,猜想表層可能存在一個明顯的分層,該分層以上聲速呈正梯度分布,平均聲速在1560—1580 m/s 之間;而該分層以下聲速較大且較為穩定.

圖12 不同頻率下的單層等效聲速 (a) r=3.5 km ;(b) r=3.7 km ;(c) r=3.9 km ;(d)r=4.1 kmFig.12.Single-layer equivalent sound speeds at different frequencies:(a) r=3.5 km ; (b) r=3.7 km ; (c) r=3.9 km ;(d) r=4.1 km .

為此,對淺地層剖面儀數據的表層分層結果進一步分析,建立更細化的沉積層分層結構,如圖13所示.設上表層存在cb11—cb12的正聲速梯度,平均厚度為h11,下表層聲速為cb13,平均厚度為h12.在反演過程中,將簡正波頻散結構反演得到的二三層平均聲速(cb2=1572 m/s,cb3=1620 m/s)作為已知,利用100—150 Hz 頻帶下的信號通過匹配場方式反演表層聲速變化情 況(即cb11,cb12,cb13),得到的結果如表3,其中cba為cb11和cb12的平均聲速.和圖12 進行對比,可以看出表層情況同單層等效聲速變化較為一致,上表層聲速呈正梯度分布,平均聲速在1562.5—1585 m/s 之間,與350 Hz下的cbe較為一致,下表層聲速較為穩定在1670 m/s附近,且h11和h12的和近似于表1中h1的結果.

表3 在100—150 Hz 頻段表層參數的反演結果Table 3.Inversion results of the surface layer at 100—150 Hz.

圖13 細化的海底分層結構Fig.13.Finer-grained seabed layering along the consider track.

實驗期間進行了海底重力取樣,如表4 為5 km 位點處采樣樣品分析結果(100 kHz 換能器測量).由表4 可以看出,0—4 m 深度內整體采樣趨勢為隨深度增加,聲速逐漸增大,其中的平均聲速為1560 m/s,另外,0 km 處海底取樣樣品分析0—2.35 m 深度內的平均聲速為1576 m/s,10 km處海底取樣樣品分析0—4 m 深度內的平均聲速為1574 m/s.均與本文反演得到的cba較為接近.

表4 5 km 站位采樣樣品分析結果Table 4.Sample results at 5 km site.

3.4 傳播損失反演衰減系數

一定距離內的傳播損失曲線是反演低頻海底衰減系數的有效途徑.將實驗獲得的聲傳播損失用于海底衰減系數反演,構造代價函數如下:

式中,N為聲傳播距離上的采樣點數;M為陣元個 數; e 代表實 際數據; c 代表拷 貝場數 據; TL 表示以f0為中心頻率的1/3 倍頻程帶寬內的平均聲傳播損失.

將簡正波頻散結構反演得到的二三層聲速及其密度和匹配場反演得到的表層聲速及其密度的均值作為已知,將各層衰減系數視作相等,利用100—500 Hz 各陣元在20 km 內共計39 個距離處的傳播損失反演衰減系數,結果如表5,平均誤差為2 dB.圖14 所示為深度8 和54 m 處陣元在225和375 Hz 處傳播損失的對比結果,其中圓圈代表實驗數據,直線代表利用反演結果計算的聲傳播損失曲線,可以看出兩者符合較好.

表5 不同頻率下反演的衰減系數及反演結果預報聲傳播損失的平均誤差Table 5.The inversion result of attenuation coefficient and their errors.

圖14 不同頻率的傳播損失 (a) 225 Hz;(b) 375 HzFig.14.Comparison of theoretical and experimental transmission-loss at different frequencies:(a) 225 Hz;(b) 375 Hz.

將100—500 Hz 間衰減系數隨頻率的變化關系進行擬合,如圖15 所示,為α=0.36f1.92dB/m,其中,f的單位是kHz,α與文獻[27]總結中的砂質海底條件下衰減系數與頻率的非線性關系α=(0.38±0.08)f(1.87±0.13)dB/m 基本一致.

圖15 不同頻率的衰減系數Fig.15.Inverted attenuation coefficient at different frequencies.

為了進一步對比本文反演結果的效果,建立傳統單層海底模型,根據圖12 中不同頻段下等效聲速的反演結果,對于聲速較為穩定的100—300 Hz頻段范圍,利用傳播損失反演衰減系數,反演結果及對應的聲傳播損失預報誤差如表6 所示.比較表5 和表6,可以看出相同頻率下多層海底模型反演結果預報傳播損失的誤差小于單層等效海底反演結果預報傳播損失的誤差,且多層海底聲速可以同時適用于較寬頻段.

表6 單層海底模型下反演結果Table 6.Inversion results for single-layer bottom model.

4 結論

2022 年8 月在南海北部大陸架斜坡海區開展的海上實驗中,淺地層剖面儀獲取的海底呈現出明顯的水平變化分層結構.本文據此建立的多層水平變化海底地聲模型,提出一種低頻聯合反演方法,即先利用簡正波頻散結構反演各層海底聲速、再利用匹配場方法反演對簡正波頻散結構較不敏感的表層海底聲速,最后將得到的結果用到傳播損失中反演衰減系數.形成結論如下.

1) 在復雜海底環境中,利用淺地層剖面儀可以預先獲取海底分層結構,在此基礎上進行反演,可以有效地減少待反演參數,獲得更接近實際海底情況的反演結果.

2) 根據敏感性不同,利用簡正波頻散結構和匹配場聯合反演的方式,可以提高各層聲速反演結果的可信性,反演結果可以驗證多種聲學現象.

3) 相比于等效單層海底模型,多層水平變化海底地聲模型的反演結果在較低頻傳播損失上的精度得到明顯提升,且多層海底聲速可以同時適用于較寬頻段.

感謝自然資源部第一海洋研究所提供的淺地層剖面儀數據,感謝中山大學海洋與技術學院提供的采樣數據,感謝南方海洋科學與工程廣東省實驗室(廣州)組織實施的“關鍵海區底質聲學特性調查航次”的全體海上工作人員,是他們的辛勤勞動為本論文提供了可靠的數據支撐.