高階彈性波數(shù)值模擬及逆時偏移成像研究

薛 凡，張 智，李學(xué)來，李 奇，李喆祥

(桂林理工大學(xué) 地球科學(xué)學(xué)院，廣西 桂林 541006)

0 引言

地震數(shù)據(jù)正演、逆時偏移及全波形反演等方法均需要以高精度的地震波場數(shù)值模擬算法為基礎(chǔ)。有限差分法因其計算效率高、實現(xiàn)簡單而被廣泛應(yīng)用。Alterman和Karal[1]在對均勻介質(zhì)進(jìn)行正演模擬時首先使用了有限差分法；Boore[2]用二階彈性波有限差分法實現(xiàn)了地震波在非均勻介質(zhì)條件下的正演模擬；Kelly[3]等使用有限差分法制作合成人工地震記錄；隨著研究的不斷發(fā)展，Madariaga[4]提出了一種較為先進(jìn)的交錯網(wǎng)格有限差分法；Virieux[5]首先使用交錯網(wǎng)格差分法求解一階速度-應(yīng)力方程，相對于常規(guī)網(wǎng)格解法，精度提高了4倍，成功模擬了轉(zhuǎn)換橫波在非均勻介質(zhì)中的傳播；Dablain[6]為解決計算精度與運(yùn)算效率這一矛盾，提出了聲波方程高階差分解法，其優(yōu)點是可以采用較大的空間步長來提高計算效率，同時計算精度也得到保障；董良國[7]將時間三階導(dǎo)數(shù)通過一階速度-應(yīng)力方程轉(zhuǎn)化為空間導(dǎo)數(shù)，解決了在不增加儲存空間的前提下提高計算精度的問題；印興耀[8]等提出了旋轉(zhuǎn)交錯網(wǎng)格與緊致有限差分相結(jié)合的方法，并基于模擬退火算法進(jìn)行全局優(yōu)化，有效地壓制數(shù)值頻散，提高了模擬精度。

逆時偏移的思想是1983年Whitmore[9]在SEG年會上提出的，利用逆時偏移結(jié)果拾取層位更新速度模型，通過迭代求取速度；同年McMechan[10]采用時間和空間域二階精度有限差分法求解全聲波方程，實現(xiàn)了逆時偏移；Chang[11-12]等先后采用有限差分實現(xiàn)了聲波、彈性波的疊前逆時偏移；張會星[13]與杜啟振[14]利用高階有限差分法實現(xiàn)標(biāo)量和矢量波場的疊前逆時偏移；劉伊克[15]提出利用地震多次波對鹽丘下部成像的新方法，為地震成像研究開拓了新方向；楊江蜂[16]等對碳酸鹽巖縫-洞儲層進(jìn)行疊前逆時偏移成像研究。逆時偏移是基于全波動方程，可以處理多次波現(xiàn)象，且不存在傾角限制，是對復(fù)雜介質(zhì)成像非常理想的偏移方法。

本文在前人研究成果的基礎(chǔ)上，選取不同地質(zhì)構(gòu)造模型，對正演模擬中的各模型進(jìn)行逆時偏移成像，模型試算結(jié)果表明，交錯網(wǎng)格有限差分法可以準(zhǔn)確的模擬地震波在介質(zhì)中的傳播，分析偏移成像結(jié)果，驗證了逆時偏移成像的準(zhǔn)確性和有效性。

1 正演模擬

1.1 基本方程

本文研究的介質(zhì)為各向同性介質(zhì)，選用彈性波本構(gòu)方程：

(1)

將柯西方程帶入公式(1)中，與運(yùn)動平衡微分方程聯(lián)立，再將Vx=?ux/?x、Vz=?uz/?z代入，可以得到二維情況下的各向同性介質(zhì)彈性波一階速度-應(yīng)力方程[公式(2)]：

式中：Vx、Vz表示質(zhì)點速度，σxx、σzz、σxz表示應(yīng)力向量，ρ為密度，λ、μ為拉梅常數(shù)。

(2)

同縱、橫波速度之間的關(guān)系：

1.2 離散公式

Vx、Vz、σxx、σzz、σxz的離散值分別以U、V、R、T、H代表，i、j分別表示x、z方向的網(wǎng)格節(jié)點，k表示時間節(jié)點。通過對時間導(dǎo)數(shù)采用二階近似，得到彈性波一階速度-應(yīng)力的交錯網(wǎng)格差分格式(即離散形式)[公式(3)]：

(3)

在設(shè)置交錯網(wǎng)格時，將介質(zhì)的波場條件和相應(yīng)的波場參數(shù)定義在圖1的主體網(wǎng)格(圖1a)和交錯網(wǎng)格(圖1b)上。表1網(wǎng)格點的數(shù)量表示波場的彈性分量和彈性參數(shù)在網(wǎng)格中的空間位置。

1.3 PML邊界條件

大量研究者將Berenger[17]提出的完全匹配層(PML)邊界條件應(yīng)用到地震學(xué)研究中[18]。在二維情況下，PML吸收邊界條件的主要思想是將邊界中的未知量分解成垂直方向和平行界面方向兩個部分，在不同的邊界匹配層按照不同的方向進(jìn)行衰減。用方程表示：

圖1 主體網(wǎng)格(a)與交錯網(wǎng)格設(shè)計圖(b)Fig.1 Main grid design map(a)andstaggered grid design map(b)

表1 彈性波場分量及彈性參數(shù)空間位置Table 1 Spatial location of elastic wave fieldcomponent and elastic parameter

(4)

其中：衰減因子為

(5)

式中：Vp為最大縱波速度，δ為匹配層寬度，R為理想的邊界反射系數(shù)(一般取10-6)，當(dāng)dx、dz不等于0時表示衰減，等于0時表示不衰減。創(chuàng)建完全匹配層吸收邊界條件基本的做法是在計算區(qū)域四周加入匹配層，圖2b中左右兩側(cè)匹配層4、6區(qū)域中在x方向衰減，在z方向不衰減，表示為dx≠0，dz=0；在上下匹配層的2、8區(qū)域中在z方向衰減，x方向不衰減，表示為dx=0，dz≠0；在四周1、3、7、9區(qū)域內(nèi)在x方向和z方向都衰減，表示為dx≠0，dz≠0。

圖2 衰減邊界(a)和PML邊界設(shè)計方法(b)Fig.2 Attenuation boundary design method(a)andPML boundary design method(b)

為了驗證PML邊界條件的吸收效果和適用性，建立大小為300 m×300 m的均勻介質(zhì)模型，模型參數(shù)設(shè)置為Vp=3000 m/s，Vs=2500 m/s，ρ=2000 kg/m2。網(wǎng)格間距Δx=Δz=1 m，時間采樣間隔為Δt=0.1 ms，震源采用雷克子波，震源主頻為35 Hz，震源位置在(150 m，150 m)處。設(shè)定PML邊界吸收寬度為20 m。

由圖3可見，采用PML邊界條件對均勻介質(zhì)的數(shù)值模擬效果良好，且邊界處無明顯反射，適用于本文數(shù)值模擬方法。

2 疊前逆時偏移成像

2.1 基本原理

假設(shè)地層內(nèi)存在一個反射界面，在界面位置下行波到達(dá)時間和上行波到達(dá)時間是一致的，在偏移時，將檢波器接收到的疊前或疊后地震資料，以地下真實速度沿著逆時間軸方向，從記錄道的末尾時刻延拓到零時刻，在空間波場的每一時刻的波場快照中，選擇出適用于成像條件的信息，再將這些波場信息疊加，得到反射界面的偏移成像結(jié)果。疊前偏移處理是將震源正演模擬和檢波器逆時外推同時進(jìn)行，當(dāng)兩者處在同一時刻(TS=TR)得到成像結(jié)果(圖4)。

圖3 均勻介質(zhì)模型在t=0.08 s時刻的波場快照

圖4 時間一致性原理Fig.4 Time consistency principle

2.2 接收波場的逆時延拓

根據(jù)Whitmore提出的逆時偏移基本物理思想：若函數(shù)g(x，z，t)是波動方程的一個解，那么函數(shù)g(x，z，t-Δt)也是波動方程的一個解。意味著波動現(xiàn)象可以按照時間前進(jìn)的方向來觀察事件發(fā)生的過程和最終結(jié)果。同理反之，沿著時間倒退的方向來觀察這一事件發(fā)生的過程?；谶@種思想，彈性波逆時偏移可以視為彈性波正演的逆過程。接收波場的逆時延拓是地震波逆時傳播的過程，本質(zhì)上是正演模擬的逆時計算，求解波動方程邊值問題?；诙S各向同性介質(zhì)彈性波方程進(jìn)行全波動方程逆時偏移，表示為

(6)

T代表的是接收點記錄的最大時間，Rx、Rz分別代表接收點記錄的x分量和z分量，xR、zR分別代表接收點的位置。實現(xiàn)過程可以解釋如下：逆時延拓是從檢波器所接收到的最終時間t=T開始直至t=0，逆時延拓就是隨著設(shè)置的時間步長進(jìn)行運(yùn)算，每一次計算就代表著將該時刻接收點的波場值插入到接收點的位置上，將該值作為接收點位置的波場值，即完成這一時刻的計算，當(dāng)計算完整個時間時，也就完成了逆時延拓過程。

2.3 互相關(guān)成像條件

成像結(jié)果的質(zhì)量對后續(xù)解釋有很大的影響，而成像結(jié)果的質(zhì)量取決于成像條件的好壞，因此對成像條件的研究是逆時偏移方法中最重要的內(nèi)容。自Claerbout[19]提出時間一致性原理后，經(jīng)過不斷地發(fā)展，目前疊前偏移成像所應(yīng)用的成像條件主要分為三種類型：互相關(guān)成像條件、激發(fā)時間成像條件和上下行波振幅比成像條件。本文主要使用互相關(guān)成像條件。

震源正演模擬波場和接收點逆時延拓波場零延遲互相關(guān)是實現(xiàn)時間一致性成像原理基本的方法，用方程可以表示為

(7)

式中：S為震源波場，R為接收波場，T為地震記錄的最大時間。對彈性波場而言，就是震源波場和接收波場的x分量或z分量。互相關(guān)成像振幅的單位是震源波場或接收波場振幅的平方，其強(qiáng)度取決于震源強(qiáng)度，沒有物理上的意義。

3 復(fù)雜模型

本文基于Matlab語言編程建立幾種典型的介質(zhì)模型，對地震波場進(jìn)行數(shù)值模擬計算和逆時偏移成像處理。模型大小取300 m×300 m，網(wǎng)格間距為Δx=Δz=1 m，采樣間隔為Δt=0.1 ms，采樣時間為3.0 s。采用雷克子波作為震源子波，震源主頻為35 Hz，設(shè)定的PML邊界吸收寬度為20 m，檢波器均勻分布在z=125 m處。所有成像只取x方向(水平分量)結(jié)果。

3.1 單一裂縫模型

建立單一裂縫模型(圖5a)，模型參數(shù)見表2。介質(zhì)分為三層，第一層深度為160 m，第二層介質(zhì)厚40 m，在第二層介質(zhì)中間存在一個裂縫，裂縫大小為10 m×40 m。參數(shù)與第三層介質(zhì)相同。圖5b為單一裂縫模型的地震記錄圖，震源位置在(150 m，20 m)處，圖中可見直達(dá)波，分界面反射P、PS波，透射P、PS波，由裂縫邊緣的地形尖滅點產(chǎn)生的繞射P、PS波及多次反射波等；圖5c為含單一裂縫的地質(zhì)模型在地震記錄中去除直達(dá)波和分界面反射波后的地震記錄圖；圖5d為多炮逆時偏移成像結(jié)果；圖5e為散射波成像結(jié)果。

分析可知，偏移結(jié)果對地層和裂縫有較好的還原，但是可以看到下層部位有低頻干擾。分析散射波成像結(jié)果，在去除反射波之后，裂縫形態(tài)更加明顯，與模型基本吻合，并且逆時偏移成像的位置清晰準(zhǔn)確，偏移取得了很好的結(jié)果。

圖5 單一裂縫模型正演記錄及其偏移結(jié)果

表2 單一裂縫模型參數(shù)設(shè)置Table 2 Parameter setting of single fracture model

3.2 凹陷地質(zhì)模型

建立凹陷模型(圖6a)，參數(shù)見表3。第一層深度為100 m，第二層厚度為75 m，為了研究地震波在凹陷地質(zhì)模型中的傳播，第三層介質(zhì)中嵌入一個凹陷模型，凹陷區(qū)長度從x=100 m到x=200 m，深度從z=175 m到z=225 m。介質(zhì)參數(shù)與第二層介質(zhì)相同。圖6b為凹陷地質(zhì)模型震源位置在(150 m，20 m)處時的地震記錄圖，根據(jù)第二分界面的上行反射波和凹陷區(qū)底層界面的反射波，可以判斷出凹陷區(qū)的位置和大體范圍；圖6c為凹陷地質(zhì)模型去除直達(dá)波和分界面反射波后的地震記錄圖；圖6d為多炮逆時偏移成像結(jié)果；圖6e為散射波成像結(jié)果。

偏移結(jié)果能夠準(zhǔn)確的反映出凹陷地形的分界面，與給定的凹陷地質(zhì)模型完全一致，多炮逆時偏移成像對地層分界面有較好的成像，凹陷區(qū)分界線明顯，但也存在低頻噪聲的干擾；散射波成像效果較好，不受反射波的干擾，對凹陷區(qū)域能夠準(zhǔn)確還原，位置與模型吻合。

3.3 斷層地質(zhì)模型

建立斷層地質(zhì)模型(圖7a)，模型參數(shù)見表4。斷層起始位置為(100 m，100 m)，終止位置為(200 m，200 m)，斷層面傾角為45°。斷層上、下盤中填充物質(zhì)與最下方介質(zhì)參數(shù)相同，寬度為5 m。圖7b為斷層地質(zhì)模型震源位置在(150 m，20 m)處時的地震記錄，圖中包含直達(dá)波、斷層上、下界面反射P、PS，透射P、PS波和尖滅點產(chǎn)生的繞射P、PS波等。分析各波形分布，可以看出上盤上界面反射波和繞射波混合且能量較小，不易分辨。圖7c為去除直達(dá)波和壓制反射波后的地震記錄圖；圖7d為多炮逆時偏移成像結(jié)果；圖7e為散射波成像結(jié)果。

各偏移結(jié)果對斷層上、下盤都有很好的還原，對斷層模型繞射點(地形尖滅點)也可以還原，但對下盤下界面的偏移效果不明顯，存在低頻噪聲的干擾；散射波成像結(jié)果對繞射點還原得更好，對斷層界面的還原也有一定的效果。與給定的斷層模型位置基本吻合。

4 結(jié)語

高階彈性波動方程正演數(shù)值模擬和偏移成像一直是地球物理學(xué)界的一個重要研究領(lǐng)域。既可以作為一種認(rèn)知手段來提高研究者對各種復(fù)雜介質(zhì)中地震波傳播規(guī)律的認(rèn)識，也可以作為一項檢測工具，檢驗各種正演和偏移算法的應(yīng)用效果。本文從幾種速度模型入手，通過對各模型成像結(jié)果分析并結(jié)合前人的研究成果，從中得到以下認(rèn)識：

圖6 凹陷模型正演記錄及偏移結(jié)果

圖7 斷層地質(zhì)模型正演記錄及其偏移結(jié)果Fig.7 Forward recording and migration results of geological fault model

表3 凹陷地質(zhì)模型參數(shù)設(shè)置Table 3 Parameter setting of geological depression model

表4 斷層地質(zhì)模型參數(shù)設(shè)置Table 4 Parameter setting of geological fault model

1)交錯網(wǎng)格高階有限差分法是種計算精度較高的數(shù)值模擬方法，采用這種方法得到的地震波波場信息豐富，對研究地震波在地下介質(zhì)中的傳播規(guī)律和波場特征十分有利。其缺點是難以克服頻散效應(yīng)，要解決這一問題，必須加密計算網(wǎng)格，但這會導(dǎo)致計算量的增加，使效率下降；另外，當(dāng)?shù)乇砥鸱^大或地質(zhì)構(gòu)造復(fù)雜時，有限差分法難以確保數(shù)值模擬精度，適應(yīng)性和靈活性較差。

2)通過對幾種繞射點模型進(jìn)行疊前逆時偏移成像處理，驗證了本文的疊前逆時偏移方法是切實有效的，該方法可以對理論模型的反射點、繞射點進(jìn)行準(zhǔn)確成像。對不同的模型進(jìn)行正演模擬，再進(jìn)行疊前逆時偏移成像，可以對整個正反演的過程有一個系統(tǒng)的了解。利用疊前逆時偏移方法結(jié)合速度分析以及其他的反演方法，可以解決對更復(fù)雜的介質(zhì)模型成像，進(jìn)而對實際工作中遇到的地質(zhì)問題提供參考依據(jù)，為解決地球物理反演多解性的問題提供有效的手段。