高子傲，孫建國

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

0 引言

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

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

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

1 聲速剖面水體建模

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

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

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

(1)

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

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

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

(2)

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

(3)

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

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

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

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

圖3 Munk模型地震記錄(a)，Munk模型地震記錄細(xì)節(jié)放大(b)和2.5 s波場快照(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ù)值兩個數(shù)量級的某一數(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)水層反射地震記錄中的有效信號，對閾值化處理的聲速剖面地震記錄進(jìn)行同相軸信息補(bǔ)償，消除由速度光滑梯度變化帶來的同相軸弱化影響。具體是選取能直觀判斷正確而有效的同相軸，根據(jù)其趨勢拾取同相軸坐標(biāo)信息，建立高斯級數(shù)擬合公式擬合同相軸，對聲速剖面反射地震記錄水體反射同相軸坐標(biāo)進(jìn)行弱振幅補(bǔ)償，并將原始數(shù)據(jù)中有效同相軸的最大數(shù)值作為擬合同相軸的數(shù)據(jù)數(shù)值。

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

(5)

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

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

(6)

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

圖6 對稱譜圖像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)拾取，對精度要求嚴(yán)格，必須遵守擬合后的速度分析結(jié)果與擬合前的趨勢一致的原則。在擬合補(bǔ)償同相軸后的數(shù)據(jù)中添加閾值化消除掉的值，得到還原真實地震記錄的同相軸補(bǔ)償水體地震記錄。圖7表示將水層反射地震記錄實施同相軸補(bǔ)償。其中，編號1～4表示相應(yīng)的水層反射得到補(bǔ)償后的同相軸。

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

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

(11)

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

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

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

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

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

(12)

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

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

(13)

引入Munk公式表示為：

(14)

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

(15)

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

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

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

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

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

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

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

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

圖10 真實水層反射數(shù)據(jù)(a)和利用相似系數(shù)法計算水層反射速度譜(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)和方法計算速度譜(a)和歸一化互相關(guān)和方法計算速度譜(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ù)計算的速度譜中，僅能得出一部分有效速度。由于背景噪聲影響，速度譜中出現(xiàn)大量假能量團(tuán)，無法得到準(zhǔn)確的疊加速度數(shù)值及完整的海水疊加速度分布趨勢。

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

圖12 閾值化后的真實水層反射數(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)

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

4 結(jié)論

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

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