地震波場數(shù)值模擬方法

1、第 42卷第 2期 石 油 物 探 V ol. 42,N o. 2文章編號 :100021441(2003 0220143206地震波場數(shù)值模擬方法張永剛(中國石油化工股份有限公司科技發(fā)展部 , 北京 100029摘要 :簡要總結(jié)了地震波場數(shù)值模擬的各種方法的基本原理及其主要特點(diǎn) , 對最近在該領(lǐng)域出現(xiàn)的一些方法和 研究結(jié)果做了簡要的闡述 , 并對比了各種方法的優(yōu)缺點(diǎn) 。 在此基礎(chǔ)上提出了運(yùn)用波動方程數(shù)值模擬作為基礎(chǔ) , 結(jié)合射線方法輔助識別波場類型 , 用于分析異常波的產(chǎn)生機(jī)理和出現(xiàn)特點(diǎn)的基本思想 , 這對復(fù)雜條件下的地震 勘探具有指導(dǎo)和借鑒意義 。關(guān)鍵詞 :地震波場 ; 數(shù)值模擬 ; 射線

2、追蹤 ; 有限元 ; 偽譜法 ; 正演模擬中圖分類號 :P63114+1 文獻(xiàn)標(biāo)識碼 :AOn numerical simulations of seismic w avefieldZhang Y onggang(Department of Science and T echnology Development , , ,Abstract :This paper reviews the principles of of seismic wavefield , and com 2 pares the merits and defects the and results are briefly d

3、iscussed. The author pro 2 poses to study of abnormal waves based on wave equation numerical simulation supplemented byK ey w ords :seismic wavefield ;numerical simulation ;ray tracing ;finite element ;pseudo 2spectrum ;forward m odeling 地震波場數(shù)值模擬是研究復(fù)雜地區(qū)地震資料 采集 、 處理和解釋的有效輔助手段 , 地震波場數(shù)值 模擬的主要方法包括 2大類 ,

4、 即波動方程法和幾何 射線法 。 波動方程數(shù)值模擬方法實(shí)質(zhì)上是求解地 震波波動方程 , 因此模擬的地震波場包含了地震波 傳播的所有信息 , 但其計(jì)算速度相對于幾何射線法 要慢 。 幾何射線法也就是射線追蹤法 , 屬于幾何地 震學(xué)方法 , 由于它將地震波波動理論簡化為射線理 論 , 主要考慮的是地震波傳播的運(yùn)動學(xué)特征 , 缺少 地震波的動力學(xué)信息 , 因此該方法計(jì)算速度快 。 因 為波動方程模擬包含了豐富的波動信息 , 為研究地 震波的傳播機(jī)理和復(fù)雜地層的解釋提供了更多的 佐證 , 所以波動方程數(shù)值模擬方法一直在地震模擬 中占有重要地位 。1 幾何射線法幾何射線法又稱為射線追蹤法 , 其主要理論

5、基 礎(chǔ)是 , 在高頻近似條件下 , 地震波的主能量沿射線 軌跡傳播 。 基于這種認(rèn)識 , 運(yùn)用惠更斯原理和費(fèi)馬 原理來重建射線路徑 , 并利用程函方程來計(jì)算射線 的旅行時(shí) 。 在旅行時(shí)計(jì)算中可以應(yīng)用有限差分等 方法 , 以期獲得快速的解 。 隨著地震勘探技術(shù)的發(fā) 展 , 新的射線追蹤技術(shù)也不斷涌現(xiàn) , 以滿足大的數(shù) 據(jù)處理 (如三維數(shù)據(jù) 和較高精度要求下對復(fù)雜地 質(zhì)體研究的需要 。這些技術(shù)主要研究焦點(diǎn)是如何 精確地劃分地質(zhì)體 , 如何實(shí)現(xiàn)旅行時(shí)場的快速準(zhǔn)確 的計(jì)算以及對已有方法的改良 。射線法的主要優(yōu)點(diǎn)是概念明確 , 顯示直觀 , 運(yùn) 算方便 , 適應(yīng)性強(qiáng) ; 其缺陷是應(yīng)用有一定限制條件 ,

6、計(jì)算結(jié)果在一定程度上是近似的 , 對于復(fù)雜構(gòu)造進(jìn) 行兩點(diǎn)三維射線追蹤往往比較麻煩 。為了計(jì)算波 沿射線的旅行時(shí)和波沿射線的振幅變化 , 首先都必 須知道波的傳播路徑 。 所謂射線追蹤 , 狹義來說指 的就是根據(jù)地震波的傳播規(guī)律確定地震波在實(shí)際 地層中傳播的射線路徑 。在地震學(xué)中 , 有 2類地震射線追蹤問題 :一類 是 1點(diǎn)射線追蹤 , 即已知射線初始點(diǎn)位置和初始出 射方向求地震波的傳播路徑問題 ; 另一類是 2點(diǎn)射 線追蹤 , 即已知射線初始點(diǎn)和另一觀察點(diǎn) (接收點(diǎn) 的位置 , 不知射線初始出射方向 , 求 2點(diǎn)之間的射 線路徑問題 。用射線理論制作 VSP 模型 , 不論是收稿日期 :20

7、030306。作者簡介 :張永剛 (1956 , 男 , 高級工程師 , 現(xiàn)從事科技發(fā)展規(guī)劃 和管理工作。零偏移距或非零偏移距 , 因?yàn)檎鹪春徒邮拯c(diǎn)不是同 一點(diǎn) (除零偏移距井口記錄道外 , 遇到的都是 2點(diǎn) 射線追蹤問題 。1 有關(guān)初值問題的試射法 (Shooting method 。 這種方法根據(jù)由震源出發(fā)的一束射線到達(dá)接收點(diǎn) 的情況對射線出射角及其密度進(jìn)行調(diào)整 , 最后由最 靠近接收點(diǎn)的 2條射線走時(shí)內(nèi)插求出接收點(diǎn)走時(shí) 。 2 有關(guān)邊值問題的彎曲法 (Bending method 。 這種方法從震源與接收點(diǎn)之間的一條假想初始路 徑開始 , 根據(jù)最小時(shí)準(zhǔn)則對路徑進(jìn)行擾動 , 從而求 出接收

8、點(diǎn)處的走時(shí)及射線路徑 。費(fèi)馬時(shí)間穩(wěn)定原理是彎曲法依據(jù)的基本原理 , 即在固定的 2個(gè)端點(diǎn)之間 , 波實(shí)際傳播的路徑是在 真實(shí)射線附近變動的路徑中 , 能使波的旅行時(shí)穩(wěn)定 的路徑 , 也就是說 , 時(shí)間取極值 (通常是極小值 的 路徑 。以上 2的速度變化 ,時(shí) , 計(jì)算效率低 , 。僅 , 最近幾年在該方面 的研究主要關(guān)注多值走時(shí)計(jì)算方面 。20世紀(jì) 80年代末以來 , 隨著 K irchhoff 積分疊 前深度偏移在解決復(fù)雜構(gòu)造成像中獲得一系列成 功 , 作為其算法基礎(chǔ)之一的射線追蹤方法也得到了 很大的促進(jìn)和發(fā)展 , 出現(xiàn)了大量不同于傳統(tǒng)算法的 新型算法 。 這些算法的主要特點(diǎn)在于不再局限于

9、 地震波的路徑描述 , 而是直接從惠更斯原理或費(fèi)馬 原理出發(fā) , 采用等價(jià)的波前描述地震波場的特征 。 近年來隨著三維勘探的發(fā)展 , 為了適應(yīng)三維數(shù)據(jù)處 理和復(fù)雜地質(zhì)體研究 , 新的射線追蹤算法不斷涌 現(xiàn) 。 這些改進(jìn)主要有 :在傳統(tǒng)的試射法和彎曲法的 基礎(chǔ)之上產(chǎn)生了很多改進(jìn)的射線追蹤方法 , 如波前 重建法 ; 對最小走時(shí)算法的改進(jìn) , 使之可適應(yīng)多值 走時(shí)計(jì)算 , 如慢度匹配法 。 下面對這些方法作一簡 要的介紹 。1 Vidale 方法 。 與傳統(tǒng)方法不同 ,Vidale 方法 計(jì)算的是波陣面而不是射線路徑 , 以二維情況為 例 , 用正方形的網(wǎng)格對速度模型進(jìn)行離散化 。 2 改進(jìn)的 V

10、idale 方法 。 當(dāng)介質(zhì)中存在較大的 速度間斷時(shí) ,Vidale 方法會出現(xiàn)不穩(wěn)定 , 因此 , 在 Vidale 之后 , 相當(dāng)一部分關(guān)于程函方程有限差分法 的研究主要是針對上述問題的 。其中 ,Qin (1992 在 Vidale 擴(kuò)展方陣的基礎(chǔ)上實(shí)現(xiàn)了擴(kuò)展波前的遞 推方法 , 但計(jì)算量增加很大 ;P odvin (1991 也是按擴(kuò) 展方陣的方式求取走時(shí) , 對每一個(gè)網(wǎng)格點(diǎn) , 系統(tǒng)地 比較來自各個(gè)方向的透射波 、 衍射波和首波 。 P od 2 vin 采用的擴(kuò)展方陣方式與 Vidale 方法相同 , 即從 上一方陣面的走時(shí)相對極小值點(diǎn)開始到走時(shí)相對 極大點(diǎn)結(jié)束 , 考慮到這種情況

11、 , 則只需比較 2個(gè)可 能的首波與 1個(gè)衍射波 , 與 Vidale 方法相比 ,P odvin 的方法同樣增加了很多計(jì)算量 , 但其穩(wěn)定性相當(dāng) 好 。3 Van Trier 法 。 Van Trier (1991 首先將程函 方程化為守恒型程函方程 , 然后用有限差分 (上風(fēng) 法 直接求解變換后的方程 , 進(jìn)而求出地震波場的 最小走時(shí) 。 該算法的主要局限性在于守恒通量的 計(jì)算 , , 守恒通量可能變?yōu)?, , 在實(shí)際應(yīng)用中 , 常。 這種坐標(biāo)系中實(shí)現(xiàn) 。這樣雖不需要進(jìn)行不同坐標(biāo)系下 走時(shí)的轉(zhuǎn)換 , 但為保證差分格式的穩(wěn)定性 、 相容性 等又需對源點(diǎn)及計(jì)算網(wǎng)格作特殊處理 , 所以并沒有 太

12、大的優(yōu)越性 。4 WFRT 法 。 WFRT 方法的基本出發(fā)點(diǎn)為惠 更斯原理 。 根據(jù)介質(zhì)的非均勻程度 , 將所要研究的 介質(zhì)分割成大小相等的矩形網(wǎng)格 , 每個(gè)矩形網(wǎng)格單 元內(nèi)的速度可視為均勻的 , 稱為第一次分割 ; 然后 再根據(jù)計(jì)算精度的要求將每一個(gè)矩形網(wǎng)格進(jìn)一步 分成等份的小矩形網(wǎng)格 , 稱為第二次分割 。 5 最短路徑法 。 最短路徑法的基礎(chǔ)是費(fèi)馬原 理及圖論中的最短路徑理論 。利用最短路徑的思 路求解程函方程 , 很多人都作了大量工作 , 這些方 法的基本思路相同 , 只不過具體實(shí)現(xiàn)步驟上存在差 異 。6 HWT (Huygens Wavefront Tracing 法 。 Sava

13、 和 F omel (1998 認(rèn)為 , 傳統(tǒng)求解射線方程的射線追 蹤方 法 雖 然 能 求 出 多 值 走 時(shí) (Vinje , 1992; Sun , 1992;Lambare ,et al ,1996 , 但缺乏穩(wěn)健性 , 程函方 程的有限差分法雖具有穩(wěn)健性 , 但只能計(jì)算最小走 時(shí) 。 因此 ,Sava 和 F omel 將兩者結(jié)合起來 , 提出了 一種新方法 。7 慢度匹配 (Slowness Matching 法 。 慢度匹配 法的目的仍是求解多值走時(shí)問題 ,Symes (1998 假 設(shè)地震波的射線都是下行的 , 該算法的計(jì)算量很 441石 油 物 探 第 42卷大 , 因此在目

14、前還不是一種可實(shí)用化的方法 。 8 基于圖形結(jié)構(gòu)的三維射線追蹤方法 。節(jié)點(diǎn) 走時(shí)利用的是 Bresenham 畫線算法 。 最小走時(shí)點(diǎn)查 尋中運(yùn)用了快速排序算法 , 節(jié)點(diǎn)設(shè)置方式可引入速 度界面 , 可實(shí)現(xiàn)反射波射線追蹤 。9 三維多值走時(shí)地震波場重建及格林函數(shù)計(jì) 算 。 這種方法根據(jù)三維復(fù)雜介質(zhì)地震波傳播數(shù)值 計(jì)算中出現(xiàn)的多值走時(shí)情況 , 以及地震波場重建和 格林函數(shù)計(jì)算中的困難 , 提出在相位空間拉格朗日 流形上的波場重建及格林函數(shù)計(jì)算方法 。該方法 基于波前重建思路 , 通過在相位空間中將位形空間 中具有多值且不可微的射線參數(shù)曲面展為單值可 微曲面 , 以克服地震波場射線解的復(fù)雜性態(tài) 。

15、 10 辛幾何算法 13。 在非均勻介質(zhì)中 , 常常 需要用 Maslov 漸進(jìn)理論解決求地震波路徑時(shí)出現(xiàn) 的焦散問題 。 Maslov 漸進(jìn)理論需要對空間坐標(biāo)進(jìn) 行傅氏變換 , 。 此 , 在數(shù)值計(jì)算中 ,間關(guān)系的算法 。用辛格式計(jì)算出的數(shù)值解滿足一擾動的哈密 頓系統(tǒng) , 因此總能量的誤差不會有很大的變化 。 由 于地震勘探中涉及的介質(zhì)通常是具有不規(guī)則界面 和復(fù)雜結(jié)構(gòu)的地球介質(zhì) , 而這些介質(zhì)中的地震波傳 播過程非常復(fù)雜 , 對其解析求解幾乎不可能 , 一般 只能通過數(shù)值計(jì)算來實(shí)現(xiàn) 。波場計(jì)算方法是整個(gè) 地震成像技術(shù)的基石 , 也是人們歷年來研究的焦 點(diǎn) 。 辛幾何算法是在哈密頓體系下表述動

16、力學(xué)體 系的一種方法 。辛算法是基于如下理念的算法 , 即動力學(xué)系統(tǒng) 在哈密頓表述下 , 能量守恒 (也即哈密頓能量守恒 對應(yīng)的數(shù)學(xué)表述為系統(tǒng)狀態(tài)的演化過程是辛變換 , 因而在計(jì)算系統(tǒng)狀態(tài)時(shí)的離散化方程也應(yīng)該是辛 變換 , 滿足這一條件的數(shù)值算法統(tǒng)稱為辛算法 。 比較辛格式與有限差分格式可以看出 , 有限差 分格式變換算子與辛格式的變換算子一樣 , 都滿足 辛性質(zhì) 。但是 , 辛格式保持的是相點(diǎn)的辛性質(zhì) , 即 t 時(shí)刻的波場和波場的廣義動量構(gòu)成的向量 , 而有 限差分格式保持的不是相點(diǎn)的辛性質(zhì) , 而是 t 時(shí)刻 的波場和 t -d t/2時(shí)刻波場的廣義動量構(gòu)成的向 量 。 在計(jì)算步長 d

17、t 較大時(shí) , 有限差分格式帶來的 偏差是不可忽視的 。 但是在大多數(shù)情況下 , 步長很 小 , t 時(shí)刻和 t -d t/2時(shí)刻的廣義動量差別不大 , 因 而有限差分格式是辛格式的一個(gè)很好的近似 。 2 波動方程法介質(zhì)中聲波或彈性波場的數(shù)值模擬 , 對于人們 理解波動傳播規(guī)律 , 解釋實(shí)際地震資料 , 表征地下 介質(zhì)結(jié)構(gòu)與巖性以及地球資源開發(fā)等 , 均具有重要 的理論和實(shí)際意義 。地震波動方程的數(shù)值模擬方 法主要有 :有限差分法 、 有限元法 、 反射率法 、 傅里 葉偽譜法等 , 但這些方法都各具優(yōu)缺點(diǎn) 。 傅里葉變 換法由于利用空間的全部信息對波場函數(shù)進(jìn)行三 角函數(shù)插值 , 所以能更加精

18、確地模擬地震波的傳播 規(guī)律 。 同時(shí) , 利用快速傅里葉變換 (FFT 進(jìn)行計(jì) 算 , 還可以提高運(yùn)算效率 。 該方法的主要優(yōu)點(diǎn)是精 度高 , ; 。波動方, 保證復(fù)雜地 ; 缺點(diǎn)是占用內(nèi)存和運(yùn)算量均 較大 。 有限差分法的主要優(yōu)點(diǎn)是計(jì)算速度快 , 占用 內(nèi)存小 ; 缺點(diǎn)是精度低 , 僅適合于相對較簡單的地 質(zhì)模型 。 下面對幾種方法的特點(diǎn)及最近的發(fā)展作 一簡要的介紹 。2. 1 有限差分方法有限差分?jǐn)?shù)值仿真技術(shù)是聲波或彈性波場數(shù) 值模擬中最為流行的方法之一 , 然而傳統(tǒng)的有限差 分方法在求解波動方程時(shí) , 會產(chǎn)生不期望的數(shù)值頻 散或稱網(wǎng)格頻散 , 導(dǎo)致了數(shù)值模擬結(jié)果分辨率的降 低 。 究其根

19、源 , 是因?yàn)榛诓▌臃匠痰挠邢薏罘智?解過程 , 通常是利用一離散化的有限差分方程去逼 近波動方程 , 從而使得相速度變成了離散空間間隔 的函數(shù) 。 因此 , 當(dāng)每一波長內(nèi)空間采樣太少 (即空 間網(wǎng)格太粗 時(shí) , 就會產(chǎn)生數(shù)值頻散 。所謂數(shù)值頻散實(shí)質(zhì)上是一種因離散化求解波 動方程而產(chǎn)生的偽波動 , 這種頻散既不同于波動方 程本身引起的頻散 , 也不同于因波傳播的速度 、 頻 率和角度而引起的頻散 , 它是有限差分方法求解波 動方程時(shí)所固有的本質(zhì)特征 , 無法避免 。 為了消除 這種數(shù)值頻散 , 人們進(jìn)行了大量的研究 , 他們的結(jié) 論是基本一致的 , 即為了消除數(shù)值頻散 , 在使用二 階有限差

20、分方法時(shí) , 每個(gè) 功率對應(yīng)的波長至少 必須使用 11個(gè)網(wǎng)格點(diǎn) , 而四階有限差分則可用二 階差分網(wǎng)格點(diǎn)數(shù)的一半 。王才經(jīng) 4(1990 也得到 了基本一致的結(jié)論 。541第 2期 張永剛 1地震波場數(shù)值模擬方法二階精度有限差分解實(shí)際上是對波動方程中 的波函數(shù) (位移函數(shù) , 位函數(shù)或其他函數(shù) 在其 T ay 2 lor 展開式中略去了二階以上的高階項(xiàng)后得到的一 種近似表達(dá)式 。 為了提高有限差分的精度 , 減少數(shù) 值模擬中頻散的發(fā)生 , 需要進(jìn)一步保留一部分高階 項(xiàng) , 也就是說 , 要保留三階 、 四階甚至更高階的微 量 , 使差分方程的近似程度更加接近實(shí)際的偏微分 方程 。同時(shí) , 波動

21、方程里還含有關(guān)于時(shí)間的二次導(dǎo)數(shù) 項(xiàng) , 但時(shí)間導(dǎo)數(shù)項(xiàng)的高精度差分卻不能象空間差分 那樣通過增加參與計(jì)算的樣點(diǎn)數(shù)來得到 , 因?yàn)樵跁r(shí) 間軸上 , 函數(shù)值是永遠(yuǎn)向前迭代的 , 而不能利用下 一個(gè)時(shí)刻或以后幾個(gè)時(shí)刻的函數(shù)值 。 因此 , 高精度 的時(shí)間差分難以在時(shí)間軸上直接實(shí)現(xiàn) 。但可以通 過適當(dāng)?shù)淖儞Q , 借助各空間坐標(biāo)點(diǎn)上的波場值及導(dǎo) 數(shù)值的組合 , 來達(dá)到提高差分精度的目的 。 通常情 況下 , 只取等號右端的前三項(xiàng) ,精度 ,高精度 ,間差分運(yùn)算將時(shí)間差分的精度也提高到四階或更 高階 。從上述研究結(jié)論中 , 我們似乎可以得出 , 通過 提高有限差分格式的精度 , 可充分減少每一個(gè)波長 內(nèi)的采

22、樣點(diǎn)而不產(chǎn)生數(shù)值頻散 , 因此問題可以輕而 易舉地得到解決 。 但事實(shí)上并非如此 , 因?yàn)樵偬岣?差分格式的精度 , 將會大大增加計(jì)算量 。因此 , 不 可能無限制地通過提高有限差分格式的精度來消 除在粗網(wǎng)格計(jì)算條件下的網(wǎng)格頻散 。在流體動力學(xué)連續(xù)方程的求解 中 , Boris 和 Book 5(1973 ,Book 6等 (1975 發(fā)展了一種通量校 正傳輸方法 (FCT , 并將他們的 FCT 方法應(yīng)用于聲 波方程的求解中 , 有效地壓制了在粗網(wǎng)格情況下差 分計(jì)算產(chǎn)生的數(shù)值頻散 。楊頂輝 7等 (1997 將這 種 FCT 技術(shù)與求解各向異性波動方程組的有限差 分方法結(jié)合 , 獲得了一種適

23、用于求解各向異性介質(zhì) 中二階聲波和彈性波方程的 FCT 有限差分算法 。 其過程主要包括 3步 :有限差分計(jì)算 、 平滑方程的 解和反漫射處理 。正如上面所述 , 在使用有限差分求解波動方程 時(shí) , 如果每一波長的采樣點(diǎn)太少 , 會產(chǎn)生數(shù)值頻散 。 因此 , 為了消除差分計(jì)算引起的頻散 , 必須作漫射 通量的計(jì)算 , 通過漫射通量來修改差分方程的解 (即平滑方程的解 。由于這種平滑處理是對計(jì)算 區(qū)域內(nèi)的每一個(gè)網(wǎng)格點(diǎn)實(shí)施的 , 從而一定程度地壓 制了真實(shí)的波動 , 造成振幅精度的部分損失 。 但這 種精度的損失并不影響 FCT 技術(shù)的推廣 , 因?yàn)閾p 失是不明顯的 。 漫射是以守恒方式實(shí)現(xiàn)的 ,

24、 即任何 時(shí)候在任一點(diǎn)去掉一部分波場時(shí) , 就必須有相同數(shù) 量加回到某一地方 。因此 , 一旦解發(fā)生漫射 , 就必 須在任何不需要漫射的地方引入反漫射以抵消漫 射 , 使真實(shí)波場得到恢復(fù) 。因此 , 對解做反漫射修 正 , 比漫射計(jì)算更復(fù)雜 , 這通常是通過一非線性局 部搜索來實(shí)現(xiàn)的 。有限差分法是一種最常用的正演模擬方法 , 現(xiàn) 已比較成熟 , 正向提高精度的方向發(fā)展 。 1970年 , Alterman 8等首先將有限差分法應(yīng)用于地震波動方 , 9度 , (空間和時(shí)間 彈 , , 將 Virieux 10,11的方法推廣到空間四階 、 時(shí) 間二階的情況 。 Crase 12則發(fā)展了精度可達(dá)

25、任意階 的高階交錯(cuò)網(wǎng)格法 , 但其計(jì)算量和內(nèi)存要求比低階 有限差分法大幅度增加 。 Magnier 13等提出了最小 網(wǎng)格有限差分法 , 它能壓制非最小網(wǎng)格的人為現(xiàn) 象 。 周家紀(jì)和賀振華 14用大網(wǎng)格快速差分算法模 擬地震波的傳播 , 其空間網(wǎng)格可以取得很大 , 達(dá)到 每最短波長只需 3個(gè)網(wǎng)格點(diǎn) , 大大縮短了計(jì)算時(shí) 間 , 等等 。2. 2 偽譜法有限差分法的算法簡單快速 , 但難以克服頻散 效應(yīng) , 而要解決頻散問題 , 須加密數(shù)值計(jì)算的網(wǎng)格 , 這勢必會導(dǎo)致計(jì)算量增加 , 效率下降 。因此 , 選擇 一種既能精確計(jì)算 , 又有較高計(jì)算效率的方法就顯 得非常必要 , 而偽譜法正好符合這種

26、要求 。 偽譜法是一種有效的數(shù)值模擬方法 , 權(quán)衡精度 和效率 , 有其他方法不可替代的優(yōu)點(diǎn) 。 在二維介質(zhì) 中 , 用偽譜法做波傳播的正演數(shù)值模擬由來已久 , 由于條件限制 , 以前的研究者僅限于二維算法程序 的開發(fā)和研究 。 20世紀(jì) 80年代末 , K osloff 1518等 人用三維聲波方程和三維彈性波方程做均勻各向 同性介質(zhì)中波傳播的模擬 , 并與解析計(jì)算結(jié)果和超 聲物 理 模 擬 進(jìn) 行 比 較 , 證 明 了 方 法 的 正 確 性 。 Reshef 18等用偽譜法作了三維聲波模擬 , 所設(shè)計(jì)的 模型在水平和垂直方向有任意的密度和波速 , 并應(yīng) 用了吸收邊界條件 。 Nie1s

27、en 等將偽譜法用到了曲 641石 油 物 探 第 42卷線網(wǎng)格的二維聲波方程的模擬中 。 T essmen 等則用 偽譜法在存在表面地形的情況下模擬了彈性波的 傳播 。 1987年 , T al 2E zer 等提出了新的快速展開 法 , 該方法可看作譜法類中的另一種 。 與偽譜法不 同的是 , 其時(shí)間導(dǎo)數(shù)用契比雪夫展開法來計(jì)算 , 可 取較大的時(shí)間步長 , 當(dāng)其與高階的空間微商結(jié)合 時(shí) , 有較高的精度 。 該方法實(shí)質(zhì)上避免了偽譜法計(jì) 算中時(shí)間精度是二階而空間精度可達(dá)無限階的不 平衡性 。 但是該法處理吸收邊界和自由表面邊界 條件比較困難 , 計(jì)算量和占用內(nèi)存也較大 。張文生等 19(19

28、98 用偽譜法進(jìn)行了二維橫向 各向同性介質(zhì)波動方程的正演模擬 , 特別是對邊界 吸收問題作了有效的處理 , 證明了偽譜法是有限差 分法近似階數(shù)趨于無限時(shí)的極限 , 它用快速傅氏變 換來計(jì)算空間導(dǎo)數(shù) , 計(jì)算精度要高于有限差分法 。 為滿足快速傅氏變換的條件 , 偽譜法的空間網(wǎng)格點(diǎn) 數(shù)必須限制取 2的正整數(shù)冪對偽譜法提出來的 ,性的邊界條件 ,界反射或一些假像 , 2倍的計(jì)算量 來消除這些效應(yīng) , 但這種處理不能用于其它的正演 模擬方法中 。 當(dāng)有必要考慮底面邊界反射時(shí) , 可采 用衰減型吸收邊界條件和反周期擴(kuò)展相結(jié)合的方 法來消除底邊界的效應(yīng) 。同時(shí) , 在偽譜法計(jì)算中 , 震源可更集中于一點(diǎn)

29、 , 這對提高精度有好處 ; 還可 取相對較大的空間網(wǎng)格步長 , 以便用來減少計(jì)算 量 。 另外 , 對波動方程系數(shù)間斷的情況 , 偽譜法比 有限差分法更有效 , 這也是偽譜法比差分法精度高 的另一原因 。2. 3 有限元及譜元法有限元法也是正演模擬的有效手段 。由于剖 分的任意性及它所依據(jù)的變分原理 , 對含有多種介 質(zhì)和自然邊界條件的處理 , 非常方便有效 , 已成為 解決地震波傳播數(shù)值模擬的一種重要方法 。它是 目前為止最精確的一種正演模擬方法 , 但它對計(jì)算 機(jī)內(nèi)存要求很高 , 計(jì)算量大 , 讓人難以承受 。 同時(shí) , 為滿足不同研究的需要 , 應(yīng)當(dāng)設(shè)置不同的邊界條 件 。 比如要研究

30、波在界面處的反射強(qiáng)度時(shí) , 就應(yīng)當(dāng) 在界面處加力 , 加力的時(shí)間為炮點(diǎn)到加力點(diǎn)的走 時(shí) , 該走時(shí)可以由射線追蹤來獲得 。有限差分法和有限元法的主要缺點(diǎn)在于對高 頻分辨的限制 , 它們對地震勘探中典型的速度和頻 率 , 計(jì)算中需要大量的網(wǎng)格點(diǎn) , 而偽譜法則相對更 有效 。有限元法的主要優(yōu)點(diǎn)是適宜于模擬任意地質(zhì) 體形態(tài) , 可以任意三角形逼近地層界面 , 保證復(fù)雜 地層形態(tài)模擬的逼真性 。 近年來 , 國內(nèi)外許多著名 學(xué)者在復(fù)雜介質(zhì)中地震波傳播問題的研究方面已 作了很多杰出的工作 。以吳如山等為代表發(fā)展了 廣義屏方法 , 基于 DeW olf 近似推導(dǎo)了雙域方程的 表達(dá)式 , 將復(fù)雜介質(zhì)用背景

31、速度和擾動速度近似 , 相應(yīng)的波場用背景場和擾動場近似 , 背景場在波數(shù) 域中延拓 , 擾動場在空間域延拓 , 同時(shí)處理邊界等 。 這種方法對于地下復(fù)雜的介質(zhì)己取得明顯的效果 。 而對于地表復(fù)雜的介質(zhì)也已開始注意 。 Stig Het 2 holm 等 (1998,2001 用交錯(cuò)網(wǎng)格高階有限差分法解 速度 應(yīng)力方程 , 研究遠(yuǎn)場 P,的 K omatisch 等 (1998,2000 用譜元法實(shí)現(xiàn) 了起伏地表二維和三維彈性波數(shù)值模擬 , 討論了復(fù) 雜地表引起的各種干擾 , 該方法具有費(fèi)時(shí)少 、 精度 高等優(yōu)點(diǎn) , 并利用該方法模擬了按正弦規(guī)律起伏的 海底反射 、 透射特點(diǎn) 。 流體區(qū)域和固體

32、區(qū)域分別用 速度場的標(biāo)量方程和位移場的矢量方程來表征 , 兩 區(qū)域的交界面處則通過一種邊界積分來得到邊界 條件 。 模擬記錄上波場信息豐富 , 包含了反射波 、 透射波 、 轉(zhuǎn)換波 、 首波和多次波等多種成分 。作者 認(rèn)為該方法將來可能用于研究其他非均勻流體的 波場特征 。 柯本喜等 (2001 用基于混合單元 (三角 形單元和四邊形單元 的有限元法實(shí)現(xiàn)了二維起伏 地表聲波數(shù)值模擬 , 并給出了較理想的人工邊界條 件 , 理論和實(shí)際模型的模擬結(jié)果令人鼓舞 , 說明有 限元法也可用于研究起伏山區(qū)波傳播問題 。他們 正在研究 , 當(dāng)表層結(jié)構(gòu)復(fù)雜 , 存在大速度反差時(shí) , 出 現(xiàn)的網(wǎng)格頻散和計(jì)算不收

33、斂問題以及地表引起的 強(qiáng)繞射干擾淹沒有效反射問題 。 可以看出 , 復(fù)雜介 質(zhì)中波的傳播問題是當(dāng)前的一個(gè)研究熱點(diǎn) , 其中地 下復(fù)雜問題己取得重大進(jìn)展 , 而地表復(fù)雜或地表 、 地下都復(fù)雜的問題雖已引起人們的強(qiáng)烈注意 , 但研 究程度較低 , 地形變化 、 表層不均勻引起的強(qiáng)繞射 干擾 、 發(fā)散等問題尚侍解決 , 復(fù)雜表層問題己成為 人們主攻的難點(diǎn) 。741第 2期 張永剛 1地震波場數(shù)值模擬方法148 石 42 卷 油 物 探 第 283 3 結(jié)束語 射線方法雖然是一種經(jīng)典的方法 ,但不能誤認(rèn) 為它是低級的和行將被淘汰的方法 。近幾年射線 法更是取得了很多新的發(fā)展 ,例如最近提出的高斯 射線

34、束法 、 傍軸射線法 、 馬斯羅夫漸進(jìn)理論等就是 明顯的例證 。另外射線法和波動方程理論相結(jié)合 , 解決地震學(xué)和地震勘探中的各種正 、 反演問題 , 也 是目前發(fā)展的一種趨勢 。波動方程法是研究地震 波場特點(diǎn)的最根本方法 。雖然費(fèi)時(shí) , 易引進(jìn)干擾 波 ,但其波場齊全 ,信息豐富 ,因此對于研究復(fù)雜條 件下的各種波場最為有效 , 具有廣闊的發(fā)展前景 。 隨著計(jì)算機(jī)運(yùn)算能力的提高和各種新方法的誕生 , 波動方程法會越來越顯示出它的優(yōu)勢 。若能 2 種 方法同時(shí)運(yùn)用 , 便可以結(jié)合射線方法簡單易行 , 波 場清晰的特點(diǎn) ,用于輔助識別波動方程法模擬出的 波的類型 ,分析異常波的產(chǎn)生機(jī)理和出現(xiàn)特點(diǎn)

35、。這 對于推動地震勘探技術(shù)的進(jìn)一步發(fā)展將具有積極 的作用 。 參 考 文 獻(xiàn) 1 陳景波 , 秦孟兆 . 辛幾何算法在射線追蹤中的應(yīng)用 J . 7 楊頂輝 , 滕吉文 . 各向異性介質(zhì)中三分量地震記錄的 FCT 有限差分模擬 J . 石油地球物理勘探 ,1997 ,32 ( 2 : 181190 8 Altermen Z S ,Loewenthal D. Seismic wave in a quarter and three quarter p1ane J . Geophysics J Roy Astr Soc ,1970 ,20 (1 :101126 9 Alford R M , Kelly

36、 K R ,Boore D M. Accuracy of finite differ2 1974 ,39 (6 :834842 10 Virieux J . P2SV wave propagation in heterogeneous media : ing of elastic wave equation J . Expanded Abstacts of 60th SEG Annual Meeting ,1990. 987991 ma1 gridsJ . Geophysics ,1994 ,59 (9 :1 4351 443 14 周家紀(jì) ,賀振華 . 模擬地震波傳播的大網(wǎng)格快速差分算 12

37、 Crase E. High2order ( space and time finite2difference model2 13 Magnier S A ,Mora P ,Tarantola A. Finite difference on mini2 ence modeling of the acoustic wave equation J . Geophysics , 數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用 ,2000 ,43 (4 : 254265 2 羅明秋 , 劉洪 , 李幼銘 . 地震波傳播的哈密頓表述及辛 幾何算法 J . 地球物理學(xué)報(bào) ,2001 ,44 (1 :120127 3 秦孟兆 , 陳景波 . Maslov 漸近理論與辛