地震波傳播數值模擬方法及其應用

                                                   裴正林

摘要 本文概述了地震數值模擬及其理論基礎，闡述了地震數值模擬方法及特點。綜述了地震數值模擬方法、三維建模和計算機硬件平臺等方面的現狀和進展。最后，給出了地震波數值模擬在地震勘探方法研究、地震觀測系統優化設計、地震數據處理、地震資料解釋、開發地震等方面的應用。

關鍵詞 地震數值模擬，研究進展，應用領域，三維建模, 集群計算

1 地震數值模擬及其理論基礎

1.1 地震數值模擬

地震數值模擬是地震勘探和地震學的重要基礎。所謂地震數值模擬就是在假定地下介質結構模型和相應物理參數已知的情況下，模擬研究地震波在地下各種介質中的傳播規律，并計算在地面或地下各觀測點所應觀測到的數值地震記錄的一種地震模擬方法。這種地震數值模擬方法已在地震勘探和天然地震領域中得到廣泛的應用。它不但在石油、天然氣、煤、金屬和非金屬等礦產資源及工程和環境地球物理中得到普遍的應用，而且在地震災害預測、地震區帶劃分以及地殼構造和地球內部結構研究中，也得到相當廣泛的應用。地震數值模擬在地震勘探和地震學各工作階段中都有重要的作用。在地震數據采集設計中，地震數值模擬可用于野外地震觀測系統的設計和評估，并進行地震觀測系統的優化。在地震數據處理中，地震數值模擬可以檢驗各種反演方法的正確性。在地震數據處理結果的解釋中，地震數值模擬又可以對地震解釋結果的正確性進行檢驗。

地震數值模擬與地震物理模擬同屬于地震正演模擬，即在已知地下介質結構和物理參數模型情況下，預測在地面或井中各觀測點所應觀測到的地震記錄。而地震勘探的目的則是根據地面或井中各觀測點所觀測的地震記錄來刻畫地下介質結構模型，并描述其狀態或巖性。這是一個反演過程，但是，這個反演過程是建立在地震正演模型的基礎上的。因此地震數值模擬不僅可以進行地震正演模擬研究，同時也是地震反演的基礎。

地震數值模擬的發展非常迅速，現在已經研制出了各種各樣的地震數值模擬方法，并均在地震勘探和地震學中得到廣泛而有效的應用，這些地震數值模擬方法可以歸納為地震波方程數值解法、積分方程法和射線追蹤法三大類（Carcione,2002）。

1.2 地震數值模擬的理論基礎

地震數值模擬是在已知地下介質結構情況下，研究地震波在地下各種介質中傳播規律的一種地震模擬方法，其理論基礎就是表征地震波在地下各種介質中傳播的地震波傳播理論。

上述三類地震數值模擬方法相應的地震波傳播理論的數學物理表達方式也不盡相同。地震波方程數值解法是建立在以彈性或黏彈性理論和Newton力學為基礎的雙曲型偏微分方程—地震波傳播方程的理論基礎上的。由于地下介質性質不同，其相應的地震波傳播方程也不同。如聲學介質中的聲波波動方程；彈性介質中的彈性波波動方程；黏彈性介質中的黏彈性波波動方程；孔隙彈性介質（雙相或多相介質）中的雙相（或多相）介質彈性波方程；各向異性介質中的各向異性彈性波波動方程等。積分方程法是建立在以Huygens原理為基礎的波疊加原理基礎上的，其數學表達形式為波動方程的Green函數域積分方程式和邊界積分方程式。射線追蹤法是建立在以射線理論為基礎的波動方程高頻近似理論基礎上的，其數學表達形式為程函方程和傳輸方程。

由于地震波傳播方程在三維復雜介質中地震波傳播的廣泛適應性及地震波方程數值解法在地震波數值模擬中應用的廣泛性和有效性，下面，我們將著重研究地震波方程數值解法，同時，對積分方程和射線追蹤法也進行適當的討論。

2 地震數值模擬方法及其特點

2.1 地震數值模擬方法

地震數值模擬是以地震波傳播理論為基礎的。描述地震波在各種介質中傳播的波動方程屬于雙曲型變系數偏微分方程。地震波波動方程的定解問題（或正演方程）包括：微分算子；微分算子的系數（可變或常系數）；震源項；邊界條件；初始條件等。

地震波正演過程是由“因”推“果”的確定過程，求解微分方程可以算出系統中表示狀態的參量隨時間的演化。地震波的正演過程數學上可描寫為：

其中，—合成地震數據向量；—正演算子；—模型向量。的精度受的離散化精度和正演算子計算精度的影響。

地震波傳播方程正演問題的研究內容:

1. 地震波數值模擬原理

2. 地震波數值模擬算法

3. 計算程序的質量。它受三個主要因素影響，即計算精度，計算效率和計算的穩定性。

另外，現代計算機的飛速發展極大地促進了地震數值模擬研究的進展。

現代地震波數值模擬是利用先進的計算機技術將野外一定規模的復雜地質構造、地質體和復雜巖性通過可視化建模技術建立二維、三維地質-地球物理模型，并用先進的數值模擬方法對野外地震勘探方法進行模擬（如圖1所示）。

2.2 地質—地球物理模型離散化

不同的研究目的，構成地質—地球物理模型的地球物理量也不同。對于均勻各向同性介質聲波方程，地質—地球物理模型可表示為而對均勻各向同性介質彈性波方程而言，其地質—地球物理模型則表示為：

即模型離散化為包含有限個模型物理參數的向量。

地質—地球物理模型的離散化是通過對模型的空間剖分來實現的。目前地質—地球物理模型的空間剖分方法可分為兩種，即正交網格剖分和非正交網格剖分。正交網格在平面上是矩形網格，而非正交網格剖分在平面上是三角形網格和不規則四邊形網格。對于地下介質進行非正交網格剖分可以充分考慮到地下介質分布的幾何形狀，并且不受邊界幾何形狀的限制。從這一點來說，基于非正交網格的數值模擬方法要優于基于正交網格的數值模擬方法。

為了準確地刻畫地下介質物理性質空間變化，模型網格剖分必須要求足夠精細，只有這樣才能以較小的誤差逼近地下介質的實際物性分布。模型網格剖分的越細，則空間網格點的數目越多，這必然會占用大量的內存量，加大計算的工作量，降低計算的效率，增加計算的成本。因此，模型離散化時既要考慮數值模擬的分辨率（或網格大?。┯忠紤]計算成本。

3 地震數值模擬的發展

地震數值模擬是模擬地震波在介質中傳播的一種數值模擬技術。隨著地震波理論在天然地震和地震勘探中的應用，地震模擬技術便應運而生，并隨著地震波理論和計算機技術的發展，地震數值模擬技術自20世紀60年代以來也得到了飛速發展，形成了目前具有有限差分法、有限元法、虛譜法和積分方程法等各種數值模擬方法的現代地震數值模擬技術。

有限差分法是偏微分方程的主要數值解法之一。在各種地震數值模擬方法中，最早出現的數值模擬方法是有限差分法。Alterman和Karal(1968)首先將有限差分法應用于層狀介質彈性波傳播的數值模擬中。此后，Boore（1972）又將有限差分法用于非均勻介質地震波傳播的模擬。Alford等（1974）研究了聲波方程有限差分法模擬的精確性。Kelly等（1976）研究了用有限差分法制作人工合成地震記錄的方法。Virieux（1986）提出了應用速度—應力一階方程交錯網格有限差分法模擬P—SV波在非均勻介質中的傳播。交錯網格方法提高了地震模擬的精度和穩定性，并消除了部分假象。Dablain（1986）提出了應用高階有限差分算子進行標量波模擬的方法，此后，Bayliss等（1986）、Levander（1988）采用四階空間有限差分算子計算彈性波的地震記錄。裴正林等（2002b，2003）研究了非均勻介質中地震波傳播的交錯網格高階有限差分法。Graves（1996）給出了三維速度—應力方程交錯網格有限差分法彈性波傳播的模擬方法。三維速度—應力方程交錯網格有限差分彈性波數值模擬方法使得地震數值模擬對地震波在彈性介質中傳播模擬能力有了明顯提高。為了進一步模擬地震波在非完全彈性的實際地層中的傳播，Carcione等（1988）提出了黏滯聲波在地層中傳播的模擬方法。Tal-Ezer等（1990）進行了線性黏彈性介質中地震波傳播的方法研究。Robertsson等（1994）給出了黏彈性波有限差分模擬方法。而為了模擬地震波在實際地層中的各向異性性質，Mora（1989）研究了三維各向異性介質中地震波的模擬方法。Igel等（1995）給出了各向異性介質中地震波傳播的有限差分模擬方法。此后，為了模擬地震波在雙相或多相孔隙彈性介質中的傳播，Dai等（1995）給出了非均勻孔隙介質中速度—應力方程有限差分地震模擬方法。Carcione和Quiroga-Goode（1996）研究了孔隙彈性介質中Biot縱波的數值模擬問題。OZdenvar和McMechan（1997）給出了孔隙彈性介質中地震波的交錯網格有限差分模擬方法。Carcione和Helle（1999）提出了孔隙黏彈性介質中地震波傳播的交錯網格有限差分模擬方法。何樵登等（1996）研究了橫向各向同性介質中地震波數值模擬方法。董良國等（2000）給出了一階彈性波方程交錯網格高階差分解法。裴正林等（2004）給出了三維任意各向異性介質中彈性波傳播的的交錯網格高階有限差分模擬方法。

上述的黏彈性、各向異性和孔隙彈性介質中地震波傳播的模擬方法進一步增強了地震數值模擬在實際復雜巖性介質中的地震波傳播模擬的能力。但是一般的有限差分地震模擬方法基于笛卡爾坐標系中的規則網格，在模擬復雜地質構造和復雜地質體的復雜界面時，必然會出現階梯狀的邊界，在這種邊界上必然引起人為的虛假繞射波，為了減弱這種虛假繞射波，則必須采用精細網格，而這將導致數值模擬計算儲存量的增加和計算量的增大。為此，則必須發展基于可變網格和不規則網格的地震數值模擬方法。

Jastram 和Tessmer（1994）提出了垂直間距可變網格的彈性波模擬方法。Falk等（1996）應用可變網格有限差分法模擬井筒波的傳播。Hayashi（1999）采用可變網格對起伏地表下彈性波場進行了模擬研究。張劍鋒（1998）、Oprsal和Zahradnik（1999）提出了非均勻介質中彈性波的矩形不規則網格有限差分法模擬方法。Pitarka（1999）給出了三維各向同性介質中彈性波的矩形非規則交錯網格有限差分模擬方法。Nordstrom和Carpenter（2001）提出了曲線坐標下變形網格高階有限差分法地震數值模擬方法。裴正林等（2002a,2002b）提出了地震波傳播小波多尺度自適應網格數值解法。

隨著山地、沙漠等地表復雜區地震勘探的發展，起伏地表地震數值模擬技術受到了地球物理勘探學家的廣泛關注和重視。Hestholm和Ruud（1994）對起伏地表下彈性波傳播有限差分法數值模擬進行了研究。Tessmer和Kosloff（1994）利用Chebychev譜法對起伏地表下彈性波傳播進行了模擬。Robertsson（1996）給出了含起伏地表彈性和黏彈性介質中地震波傳播的有限差分模擬方法。裴正林等（2002b，2004）對任意起伏地表下2D彈性波傳播交錯網格高階有限差分法數值模擬進行了研究。

有限元法也是偏微分方程的數值解法之一。Lysmer和Drake（1972）最早將有限元法應用于地震數值模擬。Marfurt（1984）研究對比了模擬彈性波傳播的有限差分法和有限元法的精度。Seron等（1990，1996）給出了彈性波傳播有限元模擬方法。Padovani等（1994）研究了地震波模擬的低階和高階有限元法。Sarma等（1998）給出了彈性波傳播有限元模擬中的無反射邊界條件。王尚旭（1990）利用有限元法研究了雙向各向同性介質中彈性波傳播規律。周輝等（1997）給出了各向異性介質中波動方程的有限元數值解法。楊頂輝（1998）對孔隙各向異性介質中基于BISQ模型的彈性波傳播進行了有限元方法模擬。

虛譜法是偏微分方程的另一種數值解法。Gazdag（1981）提出了聲波傳播的虛譜法地震模擬方法。Kosloff等（1982）給出了虛譜法地震波模擬方法。Fornberg（1987）對彈性波傳播模擬的虛譜法和有限差分法進行了對比研究。Reshef等（1988）給出了三維聲波模擬的虛譜法。

積分方程法是建立在波動方程的積分表達式的基礎上的，其理論基礎是Huygens原理。積分方程法也是有限元法之后發展起來的一種地震數值模擬方法。Pao和Varatharajulu（1976）提出了彈性波散射的積分表達式。Bouchon（1987）給出了裂隙或孔洞彈性波繞射的離散波數法模擬方法。Bakamjian（1992）給出了三維地震波傳播模擬的邊界積分方程法。符力耘等（1994）提出了彈性波正演模擬的邊界元法。

射線追蹤法是建立在波動方程的高頻近似基礎上的一種地震數值模擬方法（Cerveny等, 1977）。這種方法實際只計算了最奇異部分的解，即旅行時和振幅函數的特征曲線，它們分別是程函方程和傳輸方程的解。這種方法計算效率高。但是，一些復雜的本構方程由于積分方程法和射線追蹤法不滿足假設條件而限制了這些方法的應用。

上述這些地震數值模擬方法各有優缺點。對于三維復雜構造、復雜地質體和復雜巖性地震模擬而言，交錯網格高階有限差分法其綜合性能（占內存大小、模擬精度、計算效率和并行算法實現）最好，是實用性最好的方法。

4 地震數值模擬三維建模

實際含油氣構造及地質體一般是很復雜的。為了實現對復雜油氣藏儲集體的地震數值模，就必須解決三維地質—地球物理建模問題。三維地質—地球物理建模技術是三維地震數值模擬的重要組成部分。美國SEG/EAEG和法國IFP（Institut Francais du Petrole）等地球物理組織，在20世紀90年代，先后推出了二維Marmousi、三維鹽丘和三維逆掩推覆體等模型，業已成為勘探地球物理界的標準模型。

這些模型主要針對基于聲波方程的縱波地震勘探方法而建立的。為了研究多波多分量地震勘探方法，SEG/EAEG于20世紀90年代末對原有的二維聲介質Marmousi模型進行了修改，給出了二維彈性介質Marmousi模型。上述模型主要反映的是海相沉積特點。我國“陸相油儲地球物理”學術委員會于20世紀90年代末提出了反映我國陸相沉積特點的二維陸相斷陷模型（“油儲地球物理”項目辦公室，2000）。

在進行地震數值模擬時，必須首先建立地質—地球物理模型。為了準確地建立地震數值模擬所需要的地質—地球物理模型，則需要根據工作地區的地質、鉆井、地震、測井和巖石物理資料，進行綜合分析，確定地質—地球物理模型。

5 地震數值模擬的計算機硬件平臺

地震波數值模擬是通過計算機來實現的。因此，地震波數值模擬是介于地震波傳播理論和計算機科學之間的一門交叉學科。它的發展直接受到計算機科學水平的影響。另外，現代數值分析的方法和技術對于地震波數值模擬算法的研究起主要作用。地震波數值模擬算法性能的評價主要考慮三個因素：精度、穩定性和計算效率。

高性能計算機（HPC）和大型并行計算機（MPP）已經成為衡量一個國家綜合國力和科學技術進步的重要標志。并行是實現計算能力突破的根本手段，現在的超級計算機和高性能計算機必定是并行機（如IBM SP2，我國的曙光2000系列）。然而，隨著計算機技術的發展，具有高性能計算能力不再是昂貴的巨型機的專項。Beowulf微機集群（PC-Cluster）技術是近年來涌現出來的利用廉價個人微機，針對一些特定應用來產生與巨型計算具有相同能力的技術。

地震波數值模擬，尤其是大規模三維模型的數值模擬是一個非常費機時的過程，但它在地震勘探、天然地震中確實又是一個不可或缺、非常有效的手段，因此，需要高性能計算機硬件平臺的支持。隨著個人計算機硬件性能的提高（特別是64位CPU的面市），PC-Cluster將成為三維地震數值模擬的一個可供選擇的高效而經濟的硬件平臺。

過去由于計算速度存在的瓶頸，三維地震波數值模擬方法一直未能得到實際生產應用。但近年來隨著計算機硬件的快速發展，三維地震波數值模擬方法及其應用在國際上又形成了研究的熱點。由此，可以預見，今后幾年內三維地震波數值模擬方法即將得到廣泛的實際生產應用，也必將為三維地震勘探帶來新的生機和活力。

6 地震數值模擬的應用

現代地震數值模擬在地震勘探和天然地震領域都有著廣泛的應用。在地震勘探領域中，主要應用于石油天然氣和煤等礦產資源的勘探及工程地震勘探。在天然地震中，主要用來進行地震災害預報及地殼深部構造和地球內部結構研究。

在石油和天然氣勘探方面，隨著我國石油天然氣勘探工作的不斷發展，我們面臨的勘探對象和開發條件越來越復雜，面臨著復雜構造油氣藏、巖性油氣藏和裂縫油氣藏等復雜油氣藏的勘探及尋找和開發“剩余油”的艱巨而復雜的任務。為了解決這些復雜油氣藏勘探開發問題，必須將復雜介質地震波傳播理論與地震物理模擬試驗研究結合起來，必須將復雜介質地震波數值模擬研究與地震物理模擬研究緊密結合起來。通過三維復雜介質地震模擬，包括地震數值模擬和地震物理模擬對三維復雜介質地震勘探方法、三維地震觀測系統的優化設計、三維地震數據處理方法和三維地震資料解釋方法及復雜油氣藏開發地震方法進行研究。

目前，三維復雜介質地震數值模擬的應用主要有以下五個方面。

6.1 三維復雜介質地震勘探方法研究中的應用

通過三維復雜介質——黏彈性介質、孔隙彈性（雙相）介質或各向異性介質中地震波傳播的地震波數值模擬研究，可以發現三維復雜介質中形成的一些特殊波及一些波傳播的特殊規律，因而為三維復雜介質地震勘探方法提供依據。具體應用包括以下6個方面。

1. 三維地震勘探方法研究；

2. 多波地震勘探包括轉換波法勘探方法研究；

3. 橫波勘探包括橫波分裂及橫波再分裂應用研究；

4. 砂巖油氣儲層研究；

5. 碳酸鹽巖裂縫儲層研究；

6. 消除高速層屏蔽（包括陸地和海上玄武巖層、膏巖層和凍土層）影響研究。

6.2 三維地震觀測系統優化設計中的應用

野外地震數據采集是整個地震勘探的基礎，其中激發、接收和觀測系統是構成野外采集系統的三個不可缺少的組成部分。激發和接收通過觀測系統將它們與地下構造、地質目標聯系起來。觀測系統設計的優劣直接關系到勘探任務的成敗。成功的觀測系統設計就要通過野外采集經濟、有效地獲得勘探目標的地震信息，經過后續的數據處理和資料解釋，完成勘探任務。由于三維復雜地區勘探工作量大、技術難度大、投資高，因而三維復雜地區的地震觀測系統設計就需要通過三維復雜介質地震數值模擬和地震物理模擬進行優化設計，以期達到最佳野外采集效果。具體應用如下：

1. 復雜斷塊構造三維地震觀測系統優化設計；

2. 復雜逆掩斷裂帶三維地震觀測系統優化設計；

3. 復雜潛山構造三維地震觀測系統優化設計；

4. 復雜鹽丘構造三維地震觀測系統優化設計；

5. 起伏地表檢波器組合接收優化設計。

6.3 三維地震數據處理中的應用

地震數值模擬在三維地震數據處理中的應用為各種處理方法包括靜校正、動校正、速度分析、疊加、偏移和反演方法提供各種典型三維復雜地質構造和地質體的標準數據體，用以檢驗各種地震數據處理方法的正確性。具體應用如下：

1. 三維復雜介質標準數據體的建立；

2. 對三維復雜介質地震數據處理方法進行檢驗，包括①地面地震勘探方法；②VSP方法；③井間地震層析成像方法；④各向異性介質偏移成像方法；⑤復雜鹽丘及鹽下偏移成像方法；⑥消除多次波方法。

6.4 三維地震資料解釋中的應用

地震數值模擬在三維復雜介質地震資料解釋中的應用主要為檢驗三維復雜介質地震資料解釋的正確性。首先按地震資料解釋的結果獲得的地質—地球物理模型進行地震數值模擬，如果所得到的數值模擬結果與實際地震資料相一致，則可以認為地震資料解釋是正確的（或者是多解性中的一個解）。而如果所得到的地震數值模擬結果與實際地震資料不一致，則需要修正地質—地球物理模型重新進行地震數值模擬，如此反復迭代，直到最后所得到的地震數值模擬結果與實際地震資料相一致（或其誤差在某一允許范圍內），則其最終地質—地球物理模型即可以認為是正確的地震解釋結果。由于地震解釋存在多解性問題，為了避免多解性必須將地震資料解釋與鉆井和測井資料結合起來，才能獲得正確的單一地震資料解釋結果。具體應用如下：

1. 對三維復雜介質地震資料解釋結果進行檢驗；

2. 對三維復雜介質地震資料解釋結果進行迭代修改，直到獲得正確的地震解釋結果；

3. 三維復雜構造解釋，包括①隱蔽斷塊圈閉解釋；②逆掩斷裂帶解釋；③鹽丘和鹽內幕解釋。

4. AVO、AVA屬性解釋，包括①巖性解釋；②孔隙流體性質解釋。

5. 地震層序解釋；

6. 儲層預測，包括①砂巖儲層預測；②砂泥巖薄互層儲層預測；③碳酸鹽巖裂縫儲層預測

6.5 開發地震中的應用

地震數值模擬在開發地震中也有廣泛的應用，主要包括時移地震油氣藏動態監測和油藏管理。具體應用如下：

1. 油藏精細描述；

2. 應用地震屬性進行油藏動態監測；

3. 應用時移地震進行油藏動態監測；

4. 應用三維地震正演模擬進行堵水和其它油藏管理。