薛 凡，張 智，李學來，李 奇，李喆祥

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

0 引言

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

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

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

1 正演模擬

1.1 基本方程

本文研究的介質為各向同性介質，選用彈性波本構方程：

(1)

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

式中：Vx、Vz表示質點速度，σxx、σzz、σxz表示應力向量，ρ為密度，λ、μ為拉梅常數。

(2)

同縱、橫波速度之間的關系：

1.2 離散公式

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

(3)

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

1.3 PML邊界條件

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

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

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

(4)

其中：衰減因子為

(5)

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

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

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

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

2 疊前逆時偏移成像

2.1 基本原理

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

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

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

2.2 接收波場的逆時延拓

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

(6)

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

2.3 互相關成像條件

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

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

(7)

式中：S為震源波場，R為接收波場，T為地震記錄的最大時間。對彈性波場而言，就是震源波場和接收波場的x分量或z分量?；ハ嚓P成像振幅的單位是震源波場或接收波場振幅的平方，其強度取決于震源強度，沒有物理上的意義。

3 復雜模型

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

3.1 單一裂縫模型

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

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

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

表2 單一裂縫模型參數設置Table 2 Parameter setting of single fracture model

3.2 凹陷地質模型

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

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

3.3 斷層地質模型

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

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

4 結語

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

圖6 凹陷模型正演記錄及偏移結果

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

表3 凹陷地質模型參數設置Table 3 Parameter setting of geological depression model

表4 斷層地質模型參數設置Table 4 Parameter setting of geological fault model

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

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