海洋渦旋地震成像的時(shí)域有限差分仿真

劉佳琪，郭永超，樸勝春

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

0 引言

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

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

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

1 基本原理

1.1 深海低頻聲場(chǎng)的時(shí)域有限差分理論

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

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

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

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

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

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

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

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

1.2 地震反射成像理論

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

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

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

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

2 實(shí)驗(yàn)與仿真

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

圖3 實(shí)驗(yàn)海域和巡航路線圖 Fig.3 Experimental sea area and shipping route

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

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

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

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

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

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

圖7 單次氣槍信號(hào)激發(fā)的水平陣接收信號(hào)仿真 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

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

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

圖9 第6 000 個(gè)共反射點(diǎn)處的疊加信號(hào) Fig.9 Accumulated signal at the 6 000th co-reflection point location

3 分析與討論

海水中各層的反射系數(shù)隨著海水深度變化，體現(xiàn)在地震剖面圖上條紋的強(qiáng)弱各不相同?？蓪?duì)接收信號(hào)進(jìn)行解卷積得到反射系數(shù)，再進(jìn)一步轉(zhuǎn)換為相對(duì)的聲速擾動(dòng)，合并測(cè)量的背景聲速和相對(duì)聲速擾動(dòng)初步獲得隨深度變化的聲速分布。

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

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

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

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

4 結(jié)論

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

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

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