劉佳琪，郭永超，樸勝春

（1.哈爾濱工程大學(xué)水聲技術(shù)重點實驗室，黑龍江哈爾濱 150001；2.海洋信息獲取與安全工信部重點實驗室（哈爾濱工程大學(xué)），工業(yè)和信息化部，黑龍江哈爾濱 150001；3.哈爾濱工程大學(xué)水聲工程學(xué)院，黑龍江哈爾濱 150001）

0 引言

2003 年，Holbrook 等[1]從海上地震成像剖面中提取出了水平分層的溫度分布，自此海洋地震反射技術(shù)作為重要的海洋現(xiàn)象觀測手段被引入到海洋動力學(xué)研究中，用于研究海洋動力學(xué)現(xiàn)象及其過程，國內(nèi)有學(xué)者將其稱為地震海洋學(xué)[2]。相比于利用拋棄式溫度測量儀（Expendable Bathy Thermo-graph,XBT）獲取水文數(shù)據(jù)，該方法具有更好的時空分辨率，為熱鹽侵入鋒、內(nèi)波和渦旋研究提供了新的手段。

海水地震成像剖面可以反映水團的內(nèi)部結(jié)構(gòu)，呈現(xiàn)出與地質(zhì)成像類似的分層結(jié)構(gòu)。海洋中廣泛存在的動力學(xué)現(xiàn)象導(dǎo)致了空間水平方向的不連續(xù)聲反射。為了仿真中尺度渦存在海域聲信號的傳播過程，可選擇一種直接求解波動方程的數(shù)值方法進行反射聲場的研究。時域有限差分（Time-Domain Finite Difference,TDFD）作為一種數(shù)值求解方法，首先被 Yee[3]于1966 年開發(fā)并應(yīng)用于電磁波相關(guān)領(lǐng)域，在接下來的十幾年里用于求解麥克斯韋（Maxwell）方程電磁散射問題，到80 年代后期才被引入空氣聲學(xué)[4]和水下散射聲場[5]，90 年代水聲學(xué)家提出此方法同樣適用于求解海洋的標(biāo)準(zhǔn)楔形海底的聲傳播過程[6]。時域有限差分方法[7]作為計算海洋聲場的新方法，更適合刻畫有傳播損失的海底的聲傳播問題。時域有限差分方法能用于直接求解瞬時聲壓，反映復(fù)雜海域的聲波傳播特征，具有計算速度快，占用內(nèi)存小等優(yōu)點。

本文首先建立了水中寬帶脈沖聲傳播的時域有限差分計算模型，該模型不同于通常的水聲傳播模型，考慮了聲波在水團結(jié)構(gòu)之間的反射；其次針對GO project 的海上地震勘探實驗的水文環(huán)境，仿真計算了雷克子波氣槍聲源發(fā)射脈沖信號在水團結(jié)構(gòu)之間的反射過程，獲取實驗海域不同水層間反射信號的多途到達結(jié)構(gòu)；最后利用時間結(jié)構(gòu)進行疊前偏移的地震成像。仿真獲得的成像結(jié)果與海上實驗數(shù)據(jù)的成像結(jié)果呈現(xiàn)的水團結(jié)構(gòu)相符合，初步檢驗了本算法在海洋地震成像仿真計算方面應(yīng)用的可行性。

1 基本原理

1.1 深海低頻聲場的時域有限差分理論

在過去的幾十年里，學(xué)者們提出了許多用于分析水聲中的傳播特性的計算方法，例如簡正波法[8]和拋物方程方法[9]。為了解決時域的聲傳播問題也提出了一系列的數(shù)值算法：有限差分方法和有限元算法。有限差分方法求解聲波在分層海洋環(huán)境的傳播基于兩個基本方程：運動方程和形變方程。運動方程和形變方程的表達式分別如式（1）和式（2）：

其中：P為聲壓，Vx、Vz分別為x和z方向上的質(zhì)點振速分量，κ為體積模量，ρ為密度，t為時間。在方程式（2）右側(cè)的第2 部分，考慮了介質(zhì)吸收系數(shù)α引起的聲衰減。將式（2）代入式（1），得到：

采用泰勒（Taylor）級數(shù)進行時間和空間網(wǎng)格離散：

得到中心四階的Crank-Nicolson 離散空間網(wǎng)格為

二階中心差分格式離散時間網(wǎng)格為

在空間上采用有限差分的交錯網(wǎng)格進行數(shù)值離散[10]，將ρx、Vx、ρz和Vz以及1/κ=ρc2放在同一個計算網(wǎng)格上，聲壓場在計算過程中采取與質(zhì)點振速不同的時間差分步長。時域有限差分方法空間上的網(wǎng)格劃分如圖1 所示。聲壓的第一個計算時間為n?t其中n=0，那么由聲壓場推算得到的質(zhì)點振速場的時間為（n+1/2）?t，下一個時刻的聲壓（n+1）?t不僅與（n?t）的聲壓有關(guān)也與（n+1/2）?t的質(zhì)點振速有關(guān)，此方法考慮了水平、垂直方向的聲速和密度變化，能夠更好地求解復(fù)雜海域環(huán)境的聲學(xué)問題。

圖1 時域有限差分方法空間上的網(wǎng)格劃分 Fig.1 Spatial mesh division of TDFD method

此外還應(yīng)關(guān)注邊界條件的設(shè)置，選擇在海底和兩側(cè)邊界位置添加完美匹配層[11]，吸收反射到水中的能量和無窮遠場的輻射近似。完美匹配層位置的兩個基本方程為

1.2 地震反射成像理論

地震海洋學(xué)中常用拖曳氣槍陣列的發(fā)射信號作為聲源，可激發(fā)氣泡脈動被很好壓制的低頻聲能，并被數(shù)千千米的水聽器陣列接收。由于分層海水介質(zhì)中的聲阻抗變化，各層界面處也會發(fā)生聲反射現(xiàn)象，聲源和接收陣列之間的時間延遲取決于到達不同陣元的聲線路徑。聲源和各陣元在海底的反射點近似地分布在氣槍的海底虛源與陣元的連線和海底交界的位置。如圖2（a）中，當(dāng)僅關(guān)注某一反射位置R，存在若干聲源（O）和密集排列的陣元（O′）都在此位置（R）出現(xiàn)反射，此位置即被稱為共反射點（Common Midpoint,CMP）。由于一次反射波不存在剩余時差，利用共反射點疊加后，能量增強，而多次反射波由于存在剩余時差，疊加后能量相對減弱，累加在共反射點位置的所有能量，可獲得較高信噪比的一次反射波，圖2（b）展示了共反射點反射波形成和一次反射波的時距曲線。

在深度z反射的信號到達接收位置的時延tx和隨遠離聲源的距離x的關(guān)系表示為

其中：cs為假定的恒定疊加聲速，約為1 500 m·s-1，tx被稱為信號雙程時間（Two Way Travel Time,TWTT）；反射深度為z=cstx/2，雙程時間隨著聲源到接收器的距離的增加而增加，隨著共反射點位置的改變，反射信號時距曲線形成雙曲形狀。首先，在每個共反射位置進行正常時差校正，則圖2（b）中的時距圖中反射信號雙曲線的曲率將被校正，修正的時間為?t=tx? 2z/cs，時差校正為獲取反射位置和聲速分布奠定基礎(chǔ)。此外，需要去除的還有淺層水體中從氣槍激發(fā)直接被水聽器陣列接收到的氣泡脈動。最后，疊加過程每隔10 m 選取一個共反射點進行抽道集，合成地震成像剖面圖。

圖2 一個共反射點位置的一組源和接收器路徑。Fig.2 A set of source to receiver paths corresponding to a co-reflection point location

2 實驗與仿真

2008 年在大西洋開展了多道地震勘探的實驗，在實驗過程中每間隔37.5 m 發(fā)射一次氣槍信號，接收水平陣的陣元數(shù)為192，陣元間距為12.5 m，采樣率為10 kHz。接收信號需要5 Hz 的高通濾波器去除海面波浪帶來的影響，地震圖像采用標(biāo)準(zhǔn)處理算法，利用Seismic Unix 進行速度分析、動校正和共反射點疊前時間偏移處理。巡航路線及實驗海域如圖3 所示，紅色直線為實驗海域裝載氣槍和拖曳陣的實驗船行進軌跡。XBT 測量的聲速分布如圖4 所示，獲得的地震成像結(jié)果如圖5 所示。

圖3 實驗海域和巡航路線圖 Fig.3 Experimental sea area and shipping route

圖4 XBT 測量的聲速剖面 Fig.4 Sound speed profile measured by XBT

圖5 零偏移地震成像結(jié)果 Fig.5 Seismic imaging profile with zero offset

采用地震反射成像技術(shù)，濾除氣槍信號的氣泡脈動波、消除海底地層反射和噪聲，保留真實的反射系數(shù)，校準(zhǔn)幾何擴散和偏移時間，得到的成像圖5 中能夠觀測到與XBT 測量一致的中尺度渦結(jié)構(gòu)。

利用有限差分方法根據(jù)巡航實驗過程中的XBT 測量數(shù)據(jù)，對整個研究區(qū)域的接收信號進行仿真建模。圖6 為距離聲源84 m 位置處陣元接收信號，在1.3 s 處的信號（右上方為前1.3 s 接收信號的 放大圖）可觀測到除直達波外的層間反射結(jié)構(gòu)，在1 s 位置處，渦旋的上邊緣存在較大的溫度梯度，垂直方向上阻抗差最大，因此存在高反射信號成分。

圖6 聲源附近接收到的時域波形 Fig.6 Time-domain waveform received near the sound source

共炮點表征的是一次氣槍爆炸被水聽器接收的、聲源位置相同的一系列反射路徑。仿真的共炮點接收信號如圖7 所示，氣槍位于水下8 m，接收水聽器陣位于水下10 m 處。水聽器陣分布在距離聲源84 m～2.47 km。該范圍內(nèi)海深超過4 km，為消除海底反射信號的影響，僅截取TWTT 為3 s 的部分。能夠從圖7 中明顯看到最主要能量成分為時距圖中的直達波，其次為1 s 位置由于反射形成的多道明暗間隔的雙曲線結(jié)構(gòu)。對所有位置氣槍激發(fā)的聲信號進行仿真，采用同樣的地震成像流程得到仿真接收信號的成像結(jié)果，如圖8 所示。

圖7 單次氣槍信號激發(fā)的水平陣接收信號仿真 Fig.7 Simulation of horizontal array receiving signal corresponding to a single air-gun excitation

圖8 利用TDFD 仿真獲取的地震成像結(jié)果 Fig.8 Seismic imaging profile obtained from TDFD simulation data

地震成像時壓制了多次反射波，海水的層間一次反射波得到相對加強。仿真結(jié)果也證明了，分層的海水中存在著不可忽視的反射聲場。

圖9 給出了渦旋中心位置共反射點信號疊加后的結(jié)果，在渦旋的上下邊緣處疊加信號的能量最強，阻抗梯度最大。由于投放的XBT 水平方向的測量限制，聲學(xué)仿真得到了渦旋結(jié)構(gòu)的成像剖面，但沒有實驗數(shù)據(jù)得到結(jié)果的反射結(jié)構(gòu)清晰，遺漏掉了水平方向的細節(jié)成分。后續(xù)我們將討論如何利用接收到的地震信號獲取聲速的精細分布，進一步構(gòu)建高分辨的成像結(jié)果。

圖9 第6 000 個共反射點處的疊加信號 Fig.9 Accumulated signal at the 6 000th co-reflection point location

3 分析與討論

海水中各層的反射系數(shù)隨著海水深度變化，體現(xiàn)在地震剖面圖上條紋的強弱各不相同?？蓪邮招盘栠M行解卷積得到反射系數(shù)，再進一步轉(zhuǎn)換為相對的聲速擾動，合并測量的背景聲速和相對聲速擾動初步獲得隨深度變化的聲速分布。

以其中一個接收器為例，理想的雷克子波信號形式如圖10（a）所示，利用該位置的XBT 溫度數(shù)據(jù)可計算得到近似的反射系數(shù)，在該位置除了聲傳播距離導(dǎo)致的衰減，還存在水層間的反射，水平分層的反射系數(shù)隨時間的變化如圖10（b）所示。理想的 接收信號為發(fā)射信號和反射系數(shù)的卷積形式，如圖10（c）所示。對接收信號的最小平方反卷積運算可獲得水平分層的反射系數(shù)，如圖10（d）所示，其中紅色表示正反射系數(shù)、綠色代表負的反射系數(shù)。

圖10 利用雷克子波和層反射系數(shù)獲取接收信號，對接收信號 反卷積計算反射系數(shù) Fig.10 Received signals obtained from Ricker wavelets and layered reflection coefficients and the reflection coefficients obtained by deconvolution of the received signals

最小平方反卷積基于無噪聲干擾的最小相位接收信號，因此在應(yīng)用于實際信號時應(yīng)進行低通濾波處理。對比圖10（b）和10（d）反射系數(shù)存在較好的一致性。在海水中密度變化較小的情況下，利用雙程時間和深度的關(guān)系最終可反演得到的隨深度分布的聲速結(jié)構(gòu)如圖11 所示。

圖11 利用仿真接收信號反演獲取的聲速剖面 Fig.11 Sound velocity profile from simulated receiving signal inversion

4 結(jié)論

本文建立了基于時域有限差分的波動方程數(shù)值求解方法，采用該方法仿真計算了2008 年在大西洋進行的多道地震勘探實驗中接收到的低頻寬帶脈沖信號。仿真計算過程結(jié)合了溫鹽空間分布的測量數(shù)據(jù)和實驗測線的地形起伏，計算了氣槍發(fā)射信號到達水平接收陣列時各陣元接收信號的時域波形。在此基礎(chǔ)上，參照地震勘探技術(shù)中的地震反射波法，將獲取的時域信號進行疊前時間偏移，抽取并累加共反射位置的接收信號時域波形，按共反射點位置排序后濾波，最終獲得海水中水體分層結(jié)構(gòu)的地震成像結(jié)果。仿真計算的海上地震成像剖面圖呈現(xiàn)出與多道地震勘探實驗接收陣接收信號成像類似的中尺度渦分層結(jié)構(gòu)。

海水中水平分層的結(jié)構(gòu)間存在聲阻抗差，而傳統(tǒng)水聲傳播模型忽略了水層間的微弱聲反射現(xiàn)象。本文采用的聲傳播時域建模預(yù)報方法，不僅包含水層間的折射現(xiàn)象，也考慮了水層帶來的聲反射影響，是一種能夠模擬水層分層結(jié)構(gòu)多次反射影響的聲傳播預(yù)報數(shù)值方法。

本文驗證了層間聲反射對基于反射波法地震成像技術(shù)的影響，也可為后續(xù)水體結(jié)構(gòu)地震波反射成像結(jié)果的解釋和成像新方法的研究提供有效的手段。