海水聲速剖面約束下的速度分析方法

高子傲，孫建國

吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

海水速度是非常重要的參數(shù)，對(duì)反演和成像均有影響。Munk于1947年提出了深海聲道模型的數(shù)學(xué)表達(dá)式。Brandt在1975年首次指出聲波成像研究海洋水體運(yùn)動(dòng)和精細(xì)結(jié)構(gòu)的可能性[1]。Phillips et al.對(duì)淺灘地震資料的研究中發(fā)現(xiàn)利用地震反射波可以反映海水的溫鹽結(jié)構(gòu)[2]。Holbrook et al.在對(duì)芬蘭某地地震剖面資料處理時(shí)發(fā)現(xiàn)了海水層細(xì)結(jié)構(gòu)，推斷可以利用反射地震研究海洋水體運(yùn)動(dòng)[3]。Biescas et al.在重處理Iberian--Atlantic Margin項(xiàng)目的海洋地震反射剖面時(shí)，發(fā)現(xiàn)了得到海洋學(xué)驗(yàn)證的渦旋結(jié)構(gòu)，證實(shí)可利用反射地震方法研究海洋水體運(yùn)動(dòng)[4]。宋海斌等[5]、宋洋等[6]和黃興輝等[7]利用海洋地震數(shù)據(jù)反演海水的溫鹽結(jié)構(gòu)，最終得到的海水速度變化趨勢(shì)與Munk給出的聲速剖面曲線基本一致，為深海聲道的模擬提供理論依據(jù)。魯統(tǒng)祥等[8]提出了海水地震數(shù)據(jù)水層弱信號(hào)保護(hù)技術(shù)，通過有效分離水層反射信號(hào)達(dá)到保真效果。孫建國[9]通過Munk公式和海底反射走時(shí)反演，提出計(jì)算深海水體速度的理論方法。其與常規(guī)走時(shí)反演相比,利用Munk公式進(jìn)行的海底反射走時(shí)反演最多只需要反演4個(gè)參數(shù)，并可利用交替反演將四參數(shù)反演化為單參數(shù)反演，大大地減少了走時(shí)反演的計(jì)算量。應(yīng)用地震學(xué)方法研究深海水體反射有著重要的意義。

目前，將反射地震方法引入海洋學(xué)研究逐步深化。但在利用AVO技術(shù)研究海水物性差異方面，針對(duì)水體反射地震數(shù)據(jù)的處理研究鮮少，尤其是反演水體速度的方法研究鮮少。與海底巖性變化相比，水體的速度變化不大，這導(dǎo)致水層反射信號(hào)弱化。在采集的海洋反射數(shù)據(jù)中，水層反射信號(hào)常與背景噪聲信號(hào)產(chǎn)生混疊，無法計(jì)算獲得準(zhǔn)確的海水速度信息。

本文提出計(jì)算聲速剖面疊加速度，即水層疊加速度的新方法。應(yīng)用Munk公式對(duì)海水聲速剖面進(jìn)行建模，計(jì)算擬合海水層條件下的射線路徑、走時(shí)[10--12]。將聲速剖面正演地震記錄處理后進(jìn)行速度分析[13--14]，即通過識(shí)別并補(bǔ)償閾值化處理后的聲速剖面地震記錄，使間斷的同相軸及其幅值得到恢復(fù)，強(qiáng)化水層反射信號(hào)，實(shí)現(xiàn)對(duì)水體的有效速度分析，得到水層的疊加速度。將此方法應(yīng)用到海洋地震資料中，通過計(jì)算水層反射得到了真實(shí)有效的水體疊加速度值和正確的海水速度分布趨勢(shì)。

1 聲速剖面水體建模

聲速剖面是描述深海聲道現(xiàn)象的模型。深海聲道是在一定水深產(chǎn)生的聲波波導(dǎo)。在溫度、壓力和其他因素影響下，其聲波速度隨深度變化呈現(xiàn)具有極小值點(diǎn)的二次曲線形態(tài)。而此極小值點(diǎn)所處的深度為深海聲道軸深度。

聲速剖面由Munk理論公式[9]約束：

v(z)=v0{1+ε[e-η-(1-η)]}

(1)

通過Munk公式[9]描述縱向水體速度變化，建立初始速度模型，明確水體縱向速度(圖1)。模型橫縱網(wǎng)格間距為5 m，聲道軸深度為1 000 m。在經(jīng)典Munk聲速剖面模型中，聲道軸處的速度對(duì)應(yīng)聲速剖面極小值1 500 m/s。

圖1 Munk初始速度模型(a)和剖面速度曲線(b)Fig.1 Munk initial velocity model(a)and profile velocity curve(b)

二維運(yùn)動(dòng)學(xué)射線追蹤系統(tǒng)[13]表示為：

(2)

式中：xi是位置坐標(biāo)分量；ρi為慢度矢量分量；v是速度矢量。求解此式為求解初值問題的常微分方程。應(yīng)用四階龍格--庫塔(Runge--Kutta)公式，式(2)改寫為：

(3)

式中：λ為x或z；Δτ為時(shí)間步長；Kxi,Kzi,Lxi,

Lzi為求解系數(shù)。通過式(3)可以求得網(wǎng)格內(nèi)的射線路徑。即根據(jù)震源初始值，應(yīng)用龍格庫塔法以固定的時(shí)間步長計(jì)算向外波前的傳播，若射線間距過大，則在鄰近射線間插入新射線。同時(shí)，對(duì)射線網(wǎng)格內(nèi)結(jié)點(diǎn)的參數(shù)進(jìn)行計(jì)算。初始震源為中心坐標(biāo)為(500,0)，半徑為1 m的半圓。取角度間隔為1.5°，即對(duì)0～180°范圍內(nèi)共120條射線以步長2 000進(jìn)行追蹤計(jì)算。震源子波采用主頻為25 Hz的雷克子波。圖2為對(duì)經(jīng)典Munk模型應(yīng)用波前構(gòu)筑法計(jì)算得到的模型走時(shí)及模型內(nèi)射線路徑。其中，模型網(wǎng)格數(shù)為1 000×1 000，橫縱網(wǎng)格間距為5 m。模型最大速度vmax=1 537 m/s,最小速度vmin=1 500 m/s。時(shí)間步長為0.004 s。射線理論要求地震波一個(gè)波長范圍內(nèi)滿足介質(zhì)相對(duì)光滑條件，經(jīng)典Munk模型滿足此條件。

圖2 Munk模型走時(shí)(a)及射線路徑計(jì)算(b)Fig.2 Munk model travel time(a) and ray path(b)

采用八階聲波方程有限差分方法對(duì)經(jīng)典Munk聲速剖面模型進(jìn)行數(shù)值模擬，邊界條件采用PML吸收邊界條件。圖3為聲速剖面有限差分正演。其中，聲速剖面速度模型的橫縱網(wǎng)格數(shù)為800×500，網(wǎng)格間距為5。正演的震源坐標(biāo)為(400,10)，時(shí)間采樣點(diǎn)為1 300，時(shí)間步長為0.002 5 s。震源采用雷克子波，主頻為30 Hz。道間距為5 m。最小炮檢距5 m，最大炮檢距2 000 m。圖3a為正演結(jié)果。圖3b為不同尺度細(xì)節(jié)放大的部分經(jīng)典Munk模型地震記錄。相較于強(qiáng)直達(dá)波與底層反射信號(hào)，水層反射信號(hào)極其微弱。經(jīng)過正演得到聲速剖面多炮反射地震數(shù)據(jù)后，觀察直達(dá)波與底層反射形態(tài)，用函數(shù)切割法切除直達(dá)波等淺層干擾信息及地層強(qiáng)反射信息，進(jìn)行閾值化處理，清除與速度分析計(jì)算無關(guān)的微弱的值，得到僅有水層反射部分的地震記錄(圖4)。圖4a對(duì)應(yīng)圖3a中單邊放炮切除直達(dá)波與底層反射波后的結(jié)果。

圖3 Munk模型地震記錄(a)，Munk模型地震記錄細(xì)節(jié)放大(b)和2.5 s波場(chǎng)快照(c)Fig.3 Munk model seismic record(a), zoom in details of munk model seismic record(b)and 2.5 s wave field snapshot(c)

圖4 原始聲速剖面地震記錄(a)和閾值化處理后的結(jié)果(b)Fig.4 Seismic record of initial sound velocity profile(a)and thresholding result(b)

閾值化原則為，選擇數(shù)據(jù)中明顯低于有效數(shù)值兩個(gè)數(shù)量級(jí)的某一數(shù)值作為閾值，將數(shù)據(jù)中大于該閾值的值保持原狀，將小于該閾值的值變?yōu)?。具體表示為：

(4)

式中：dst(x,y)為閾值化處理后的結(jié)果。src(x,y)為原始數(shù)據(jù)結(jié)果。

由于聲速剖面速度光滑梯度變化，使得水層反射同相軸能量微弱并出現(xiàn)同相軸間斷現(xiàn)象(圖5)。

圖5 水層反射地震記錄同相軸間斷Fig.5 Discontinuity in phase axis of water reflection seismic record

2 聲速剖面速度分析方法

為增強(qiáng)水層反射地震記錄中的有效信號(hào)，對(duì)閾值化處理的聲速剖面地震記錄進(jìn)行同相軸信息補(bǔ)償，消除由速度光滑梯度變化帶來的同相軸弱化影響。具體是選取能直觀判斷正確而有效的同相軸，根據(jù)其趨勢(shì)拾取同相軸坐標(biāo)信息，建立高斯級(jí)數(shù)擬合公式擬合同相軸，對(duì)聲速剖面反射地震記錄水體反射同相軸坐標(biāo)進(jìn)行弱振幅補(bǔ)償，并將原始數(shù)據(jù)中有效同相軸的最大數(shù)值作為擬合同相軸的數(shù)據(jù)數(shù)值。

擬合正確的標(biāo)準(zhǔn)為，擬合后的速度分析結(jié)果與擬合前的趨勢(shì)一致。這是補(bǔ)償方法實(shí)施的唯一原則。高斯級(jí)數(shù)擬合形式為：

(5)

式中：a1,a2,b1,b2,c1,c2為擬合系數(shù)。

高斯函數(shù)表示為：

(6)

高斯函數(shù)描述的對(duì)稱譜圖像可表示為圖6。其中，a1,b1,c1分別為峰高、峰位和區(qū)域?qū)挾?其中c1=2·d2。為了得到擬合系數(shù)，需將高斯函數(shù)進(jìn)行轉(zhuǎn)換以符合多項(xiàng)式結(jié)構(gòu)。將式(6)改寫為：

圖6 對(duì)稱譜圖像Fig.6 Symmetrical spectral image

(7)

(8)

矩陣表示形式為：

(9)

根據(jù)拾取的坐標(biāo)矩陣，通過求解系數(shù)矩陣得到擬合系數(shù)g0,g1,g2，從而得到所需要的a1,b1,c1。擬合精確度以決定系數(shù)(coefficient of determination)判定，定義為：

(10)

由于同相軸擬合是基于原有同相軸坐標(biāo)拾取，對(duì)精度要求嚴(yán)格，必須遵守?cái)M合后的速度分析結(jié)果與擬合前的趨勢(shì)一致的原則。在擬合補(bǔ)償同相軸后的數(shù)據(jù)中添加閾值化消除掉的值，得到還原真實(shí)地震記錄的同相軸補(bǔ)償水體地震記錄。圖7表示將水層反射地震記錄實(shí)施同相軸補(bǔ)償。其中，編號(hào)1～4表示相應(yīng)的水層反射得到補(bǔ)償后的同相軸。

圖7 同相軸補(bǔ)償Fig.7 Events compensation

抽取共中心點(diǎn)道集并進(jìn)行疊加速度譜計(jì)算。在采用的速度分析方法里，歸一化互相關(guān)和方法通過求取每個(gè)地震道內(nèi)的歸一化的離散值，得到歸一化互相關(guān)和K。K達(dá)到最大時(shí)的速度即為所求的疊加速度。歸一化互相關(guān)和[15]由下式定義：

(11)

對(duì)聲速剖面共中心點(diǎn)道集地震記錄和經(jīng)過同相軸補(bǔ)償處理后的聲速剖面地震記錄分別進(jìn)行疊加速度譜計(jì)算(圖8、9)。

a.經(jīng)典聲速剖面速度譜；b.用歸一化選擇性互相關(guān)和方法對(duì)同相軸補(bǔ)償后的速度譜。圖8 速度譜計(jì)算結(jié)果Fig.8 Velocity spectrum calculation results

a.聲速剖面速度分析結(jié)果；b.歸一化選擇性互相關(guān)和方法對(duì)同相軸補(bǔ)償后的速度分析結(jié)果。圖9 速度分析結(jié)果Fig.9 Velocity analysis results

Munk聲速剖面的速度譜中，有數(shù)量較少的有效速度團(tuán)。速度團(tuán)能量分散，不容易分辨。在有效的速度團(tuán)中，可以得到聲速剖面聲道軸處的疊加速度值，對(duì)應(yīng)海水速度分布的最小值。通過計(jì)算可驗(yàn)證速度分析結(jié)果與補(bǔ)償方法的正確性。

地下為均勻介質(zhì)，雙程旅行時(shí)時(shí)距曲線方程[15]為：

(12)

式中：介質(zhì)速度為v；t0為零炮檢距雙程旅行時(shí)間。根據(jù)速度譜計(jì)算原理，走時(shí)接近零炮檢距雙程旅行時(shí)間時(shí)的速度團(tuán)能量最大。

當(dāng)共中心點(diǎn)道集第一道的炮檢距為零炮檢距時(shí)，將聲速剖面進(jìn)行微元化，零炮檢距雙程走時(shí)T的表達(dá)式為：

(13)

引入Munk公式表示為：

(14)

將Munk聲速剖面經(jīng)典參數(shù)代入，得到：

(15)

計(jì)算式(15)積分方程，得到不同深度下雙程走時(shí)T的數(shù)值解。零炮檢距雙程走時(shí)T與深度的具體關(guān)系為：

理論模型速度分析計(jì)算所用的共中心點(diǎn)道集首道炮檢距接近零炮檢距，在大尺度的一定誤差范圍內(nèi)，可根據(jù)表1判斷聲道軸處速度團(tuán)的拾取是否正確。

表1 深度及對(duì)應(yīng)走時(shí)對(duì)照表Table 1 Depth and corresponding travel time comparison table

根據(jù)表1，利用速度分析拾取的聲道軸處的雙程走時(shí)T可以得到該處對(duì)應(yīng)的確切深度。按照聲速剖面速度公式，利用深度信息可以計(jì)算聲速剖面聲道軸處的速度。

Munk聲速剖面模型的正演時(shí)間間隔為2.5 ms，理論模型速度譜拾取的聲道軸雙程走時(shí)--疊加速度坐標(biāo)為(531.453，1 509.026)。應(yīng)用同相軸補(bǔ)償方法后的速度譜拾取坐標(biāo)為(530.246,1 508.926)。計(jì)算得到雙程旅行時(shí)間為1.328 632 5 s及1.325 615 s。根據(jù)表1所示時(shí)深關(guān)系，聲道軸的計(jì)算位置為1 001～1 003 m?？紤]首道共中心點(diǎn)道集的炮檢距并非零炮檢距的計(jì)算誤差，計(jì)算結(jié)果與實(shí)際聲道軸深度1 000 m的位置十分接近。

Munk經(jīng)典聲速剖面模型地震記錄與應(yīng)用同相軸補(bǔ)償方法計(jì)算的地震記錄經(jīng)過速度分析得到了分布一致的疊加速度。聲道軸處的均方根速度值均為1 510±。這證明了同相軸補(bǔ)償方法的有效性。

3 實(shí)際海水?dāng)?shù)據(jù)速度分析

對(duì)于實(shí)際海水反射地震記錄數(shù)據(jù)，以美國東海岸實(shí)測(cè)數(shù)據(jù)為例。1978年，美國地質(zhì)調(diào)查局(USGS)與地球物理服務(wù)公司(GSI)合作采集總長為4 813 km的反射地震數(shù)據(jù)。在美國東海岸共采集38條地震剖面。實(shí)際數(shù)據(jù)采用與海岸線基本平行的第38號(hào)測(cè)線。對(duì)多炮地震數(shù)據(jù)抽取共中心點(diǎn)道集，切除初至波及海底反射，保留水層部分記錄。水層部分地震記錄時(shí)間采樣點(diǎn)為750，時(shí)間采樣間隔為4 ms。通過對(duì)此共中心點(diǎn)道集數(shù)據(jù)應(yīng)用相似系數(shù)法、歸一化互相關(guān)和法及歸一化選擇性互相關(guān)和法進(jìn)行速度分析，得到如圖10b、圖11a、圖11b所示結(jié)果。

圖10 真實(shí)水層反射數(shù)據(jù)(a)和利用相似系數(shù)法計(jì)算水層反射速度譜(b)Fig.10 Real water layer reflection data(a) and calculation of water layer reflection velocity spectrums using similarity coefficient method(b)

圖11 歸一化選擇性互相關(guān)和方法計(jì)算速度譜(a)和歸一化互相關(guān)和方法計(jì)算速度譜(b)Fig.11 Calculation of stacking velocity panel using normalized Cross--Correlation Sum method(a) and calculition of stacking velocity panel using normalized selective Cross--Correlation Sum method(b)

海水速度在大范圍內(nèi)分布在1 500 m/s±。原始數(shù)據(jù)計(jì)算的速度譜中，僅能得出一部分有效速度。由于背景噪聲影響，速度譜中出現(xiàn)大量假能量團(tuán)，無法得到準(zhǔn)確的疊加速度數(shù)值及完整的海水疊加速度分布趨勢(shì)。

應(yīng)用同相軸補(bǔ)償方法，將原始地震記錄同相軸校正并增強(qiáng)幅值。計(jì)算并得到最終的速度分析結(jié)果。閾值化后同相軸補(bǔ)償前后的數(shù)據(jù)結(jié)果如圖12a、圖12b所示。在對(duì)應(yīng)用同相軸補(bǔ)償方法后的數(shù)據(jù)進(jìn)行相似系數(shù)法和歸一化選擇性互相關(guān)和法得到的速度分析結(jié)果如圖13a、圖13b所示。

圖12 閾值化后的真實(shí)水層反射數(shù)據(jù)(a)和同相軸補(bǔ)償后的結(jié)果(b)Fig.12 Real water layer reflection data after thresholding(a)and result after events compensation(b)

圖13 相似系數(shù)法速度分析結(jié)果(a)和歸一化互相關(guān)和法速度分析結(jié)果(b)Fig.13 Velocity analysis results of similarity coefficient method(a)and velocity analysis results of normalized Cross--Correlation Sum method(b)

通過對(duì)初始水層反射地震記錄有效反射同相軸的識(shí)別與加強(qiáng)，削弱了由于背景噪聲所造成的海水層反射同相軸間斷、錯(cuò)位。經(jīng)過速度分析后得到準(zhǔn)確的海水層疊加速度值及海水疊加速度分布趨勢(shì)。這提供了更加真實(shí)、有效的水層疊加速度計(jì)算結(jié)果。

4 結(jié)論

(1)本文在應(yīng)用Munk公式進(jìn)行聲速剖面水體建模的前提下，經(jīng)過正演得到聲速剖面地震記錄后，應(yīng)用同相軸補(bǔ)償方法，使得可以對(duì)速度光滑梯度變化的聲速剖面進(jìn)行有效速度分析。提高了海洋地震資料的分辨率，為水體反射速度分析提供方法。

(2)通過將水體反射同相軸校正及能量補(bǔ)償方法應(yīng)用到實(shí)際海洋地震數(shù)據(jù)，得到了更為真實(shí)準(zhǔn)確的海水疊加速度值及疊加速度分布趨勢(shì)，證實(shí)此方法應(yīng)用于實(shí)際生產(chǎn)中的有效性，為進(jìn)一步獲取水層速度提供了新的手段。