基于三階Adams格式的求解聲波方程的多步算法

1、 China University of Mining & Technology-Beijing創(chuàng)新項目論文一種基于三階Adams格式的求解聲波方程的多步算法 摘 要一個準確、高效、低數(shù)值頻散的正演算法能夠提高反演精度、加快反演收斂速度，因此研究地震波場正演模擬技術(shù)具有重要意義。區(qū)別于傳統(tǒng)的空間離散方法，利用空間插值, 用網(wǎng)格點處的函數(shù)值及其梯度共同逼近空間高階偏導數(shù)的方法稱為近似解析離散化方法。聲波方程通過變換，并采用近似解析離散化方法進行空間離散，從而轉(zhuǎn)變成為一個半離散化的常微分方程組，再利用三階顯式Adams格式進行時間推進，求解半離散化的常微分方程組，從而得到了一個新的求解聲波

2、方程的有限差分方法（AD-STEM）。對AD-STEM進行了理論誤差和數(shù)值誤差分析、計算效率比較和數(shù)值波場模擬。研究表明，與傳統(tǒng)方法AD-LWC比較， AD-STEM方法數(shù)值精度更高，數(shù)值頻散更低，更高效，且與解析解匹配更好。AD-STEM方法能夠通過壓制數(shù)值頻散而提高計算效率。在無可見數(shù)值頻散的條件下，AD-STEM的計算速度是AD-LWC的1.88倍，而存儲量只有其72%，更適合在粗網(wǎng)格下進行大規(guī)模地震波場數(shù)值模擬。關(guān)鍵詞：近似解析離散化方法；三階Adams格式；數(shù)值頻散；有限差分目 錄1 緒論11.1選題背景和研究意義 1.2粘彈性介質(zhì)國內(nèi)外研究現(xiàn)狀1.3有限差分國內(nèi)外研究現(xiàn)狀1.4本文

3、主要研究內(nèi)容2 粘彈性介質(zhì)的基本模型33 方法介紹33.1 Stereo-modeling方法簡介33.2 Lax-Wendroff correction方法簡介54 粘彈性介質(zhì)中的波場數(shù)值模擬64.1 波場快照64.2 波形圖64.3 SEG模型的地表地震記錄105 結(jié)論136 參考文獻15中國礦業(yè)大學（北京）2016屆本科畢業(yè)設計（論文）1 緒論1.1 選題背景和研究意義 隨著我國石油和天然氣工業(yè)的迅速發(fā)展，油氣勘探的精度和難度也在不斷地加大。這就要求我們必須深入了解地震波在地下介質(zhì)中的傳播規(guī)律，以便我們能更加準確地進行地下油氣的勘探。地震波數(shù)值模擬技術(shù)是目前地震勘探領(lǐng)域人們研究地震傳播規(guī)

4、律的一個重要手段。傳統(tǒng)地震波數(shù)值模擬技術(shù)，一般假設地層介質(zhì)是均勻和完全彈性的，但是我們發(fā)現(xiàn)：在這種假設情況下，數(shù)值模擬的結(jié)果和我們在野外觀察到的實際記錄有很大差別。造成這種差別的主要原因是因為實際的地層介質(zhì)并不是均勻和完全彈性的，我們應該用粘彈性介質(zhì)來代替，這對于我們研究地震波波場的傳播規(guī)律具有了非常重要的理論和實際意義。1.2 粘彈性介質(zhì)國內(nèi)外研究現(xiàn)狀 為了對實際地震資料的解釋以及地震波的偏移反演等處理提供可靠依據(jù)，我們需要建立不同的介質(zhì)模型來模擬地下介質(zhì)，發(fā)展相應的理論和勘探方法。從以往的模擬結(jié)果來看，以經(jīng)典彈性波理論為基礎(chǔ)得到的理論記錄實際地震波在大地中傳播時所接收到的記錄之間存在很大差

5、別，為了消除這種差別，我們通常要對球面擴散、界面的透射損失等進行補償，但是在補償過后，理論記錄上來自深層的地震波振幅仍然較小、低頻成分較豐富，這說明地下介質(zhì)對地震波有吸收衰減的作用，特別是對高頻成分，吸收更為嚴重。地震波衰減受多種因素影響，事實上現(xiàn)在沒有一種機理可以描述所有環(huán)境條件下產(chǎn)生的損耗，介質(zhì)的粘滯性就是其中的一個主要原因，為了弄清楚地震波在地層中的衰減機制，很多研究人員開始對粘彈性介質(zhì)的理論和應用方法進行研究。（1） 國外研究現(xiàn)狀 早在1845年Stokes就開始研究粘彈性介質(zhì)中的地震波，他提出要考慮具有耗損的固體就必須使這種固體具有類似于粘滯液體的性質(zhì)，因此他給固體的剪切模量附加上一

6、個剪切粘滯，并建立了由于粘滯型內(nèi)摩擦引起的能量耗損的Stokes粘彈性波動方程。但是粘彈介質(zhì)理論的發(fā)展卻是非常緩慢的，直到上個世紀40年代，美國地球物理學家N.K.Ricker首次在粘彈介質(zhì)理論中完整詳細的描述了地震波在地下介質(zhì)傳播中時的衰減問題。其后，人們紛紛開始研究粘彈介質(zhì)中地震波的傳播理論，提出了一系列的理論和研究方法。Aki和Richards(1980)建立了一種新的粘彈性介質(zhì)模型標準線性固體模型，相對于開爾芬佛格特模型和麥克斯韋爾模型，標準線性固體模型更加適合描述地震波衰減的物理機制。粘彈介質(zhì)地震波數(shù)值模擬技術(shù)的研究是從上世紀八十年代開始的，Day和Minster (1984)第一次

7、成功的在粘彈性介質(zhì)中利用Pade近似法進行2-D時間域數(shù)值模擬。EmmeriCh和Korn (1987)發(fā)現(xiàn)這種方法存在明顯的缺陷即質(zhì)量低劣和計算效率低，于是他們提出了一種新的近似稱之為“廣義標準線性固體”(GSLS)，同時對粘彈性模量的有理式進行推導，并發(fā)展二維有限差分算法使之適合標量波的傳播，提高了計算效率。CarCione(1988)等對粘滯聲波在地層中傳播的進行了模擬并推導出了模擬方程。Robertsson(1994)等設計了一種速度一應力有限交錯網(wǎng)法，并利用它研究了二維粘彈性介質(zhì)下縱波和橫波的傳播規(guī)律，隨后他們又將該算法推廣到三維介質(zhì)情況。(二)國內(nèi)研究現(xiàn)狀跟國外關(guān)于粘彈性介質(zhì)理論及

8、其數(shù)值模擬的迅速發(fā)展，國內(nèi)對這方面研究起步較晚。近年來國內(nèi)地球物理學家也開始對地震波在粘彈性介質(zhì)中傳播進行研究。張劍鋒、李幼銘(1994)把地層介質(zhì)假設為水平層狀介質(zhì)，利用混合變量粘彈性波方程，直接逐層求解位移、應力的方法進行數(shù)值模擬四。畢玉英、楊寶俊(1995)給出一種實現(xiàn)二維粘滯介質(zhì)完全波場模型計算的方法。該方法的獨到之處在于將傳播時間分解成了傳播水平時間和傳播垂直時間兩部分，平面源人射改為線源人射，無需繁雜的射線追蹤，只考慮入射角隨偏移距的變化情況便可獲得包括反射縱波、轉(zhuǎn)換波、多次波、折射波及面波等在內(nèi)的多種波場的時間一空間道集的正演模擬記錄，而且還能靈活地模擬井間和VSP地震剖面。宋守

9、根(1996)提出了一種提高地震剖面縱橫向分辨率的粘彈性介質(zhì)波場外推方法，利用該方法進行反演時，無需對積分方程作線性化近似，適應性較強，并通過實例對該理論推斷予以驗正。張建貴(1999)等針對塔里木盆地的沙漠覆蓋面大，地層埋藏深，地質(zhì)構(gòu)造幅度小的特點以大地介質(zhì)為粘彈性介質(zhì)為前提，利用粘彈性介質(zhì)波動方程聚焦成像技術(shù)得到了一套高分辨率、高信噪比和高保真的地震處理剖面。孟凡順(2000)等人根據(jù)粘滯彈性理論，推導出了粘滯彈性波的有限差分運動方程并對任意復雜地質(zhì)體粘滯彈性波進行了正演模擬；程昌鈞(2001)等系統(tǒng)研究了粘彈性介質(zhì)中波的逆散射問題及其求解過程；崔建軍和何繼善 (2001)以二維粘彈性波動

10、方程為基礎(chǔ)，研究了粘彈性波動方程F-K域的正演模擬和偏移方法的思路和理論基礎(chǔ)，并闡明了推導過程；楊午陽(2003)提出了一種利用粘彈性聲波波動方程進行偏移的新方法；奚先和姚姚(2004)通過波動方程的交錯網(wǎng)格有限差分數(shù)值模擬方法，對地震波在二維粘彈性隨機介質(zhì)中的傳播進行了模擬并研究了其波場特征；朱慧卿(2004)利用新定義的各向異性網(wǎng)格一廣義擬一致網(wǎng)格對粘彈性方程各向異性有限元方法的超收斂性進行了分析22；劉財、張智(2005)以積分本構(gòu)方程為基礎(chǔ)，應用對應原理建立波動方程，對線性粘彈體中地震波場進行偽譜法模擬；宋常瑜(2006)等采用二維粘彈性波方程交錯網(wǎng)格高階有限差分數(shù)值模擬方法進行了井間

11、地震粘彈性波場數(shù)值模擬，并對地震波的傳播規(guī)律和衰減特征進行了研究；邵志剛(2007)將以往兩種粘彈介質(zhì)中地震波模擬方法的優(yōu)點結(jié)合起來，以模型理論和積分本構(gòu)方程為基礎(chǔ)，從理論上分析了模型對地震波場的影響并采用交錯網(wǎng)格有限差分法對粘彈介質(zhì)中的地震波進行數(shù)值模擬；孫成禹(2007)針對以往粘彈性模型在描述介質(zhì)品質(zhì)因子對頻率的依賴關(guān)系方面存在不足，與實際觀測結(jié)果不符的問題，提出了一種新的粘彈介質(zhì)模型的數(shù)值方法，該方法較好地描述了品質(zhì)因子對頻率的非依賴性，可以從理論上較準確地對實際介質(zhì)的粘彈性對地震波場的影響進行研究；唐啟軍(2009)將Von Karman型的隨機各向同性背景引入粘彈性單斜各向異性波動

12、方程，并應用交錯網(wǎng)格技術(shù)進行模擬；Yong Yun-dong (2010)等人推導出了交錯網(wǎng)格有限差分的并行算法，并將其用于三維粘彈性隨機溶洞介質(zhì)數(shù)值模擬中，取得了良好的效果。1.3有限差分國內(nèi)外研究現(xiàn)狀 地震數(shù)值模擬作為地震勘探基礎(chǔ)性研究，是我們認識地震波傳播規(guī)律，檢驗各種處理方法正確性的重要工具，地震波的數(shù)值模擬是是研究地震波傳播規(guī)律的最有效工具和手段，在地震勘探和地震學各工作階段都有重要的作用，貫穿于地震資料的采集、處理、解釋的整個過程中，同時也是地震反演的基礎(chǔ)。有限差分法是最常用的一種數(shù)值模擬方法，基本原理就是在解偏微分方程時用有限差分算子代替微分，將微分方程化為相關(guān)的線性代數(shù)方程，通

13、過求解代數(shù)方程，得到偏微分方程的數(shù)值解。有限差分法數(shù)值模擬技術(shù)相對于射線方法具有更高的精度，同時比有限元方法計算量小，因此在實際應用中占很重要的地位。(一)國外研究現(xiàn)狀有限差分法數(shù)值模擬技術(shù)的發(fā)展只有幾十年的歷史，早期對于地震波傳播有限差分方法的研究都是基于二階位移方程。有限差分數(shù)值模擬方法最早由Alterman和Karal(1968)提出，他們首次將有限差分數(shù)值模擬方法應用于層狀介質(zhì)的模擬，并分析彈性波在層狀均勻介質(zhì)中的傳播規(guī)律。他們?yōu)槔貌罘址椒ń鉀Q各種勘探地震學的實際問題奠定了基礎(chǔ)，隨后差分方法被廣泛應用于地震勘探中并在應用中不斷得到發(fā)展。Boore(1972)進行了非均勻介質(zhì)地震波有限

14、差分數(shù)值模擬并研究了地震波非均勻介質(zhì)中的傳播規(guī)律。Alford (1974)對有限差分模擬的影響因素(網(wǎng)格大小和地震波傳播方向)進行分析，網(wǎng)格間距控制地震波數(shù)值模擬的計算精度，相同精度下高階差分與低階差分對網(wǎng)格間距的要求不同型;Kelly(1976)利用有限差分法進行地震記錄的合成，還對合成的地震記錄進行一些常規(guī)的數(shù)據(jù)處理分析，進一步拓寬了數(shù)值模擬在實際地震處理中的應用；前面都是基于笛卡爾坐標系下常規(guī)網(wǎng)格的離散差分，Madariaga(1976)首次提出交錯網(wǎng)格方法，并將此方法應用于彈性介質(zhì)地震波的模擬，其模擬的差分精度為O(t2+x2)。virieux (1984)提出了應力一速度一階方程交

15、錯網(wǎng)格有限差分法，用一階速度一應力方程來代替二階位移方程，由于一階速度一應力方程不需要對彈性參數(shù)進行求導，因此它更加適用于復雜介質(zhì)的模擬。 Levander對Virieux (1988)的方法進行了發(fā)展，采用四階空間有限差分算子模擬彈性波的地震記錄，進一步提高模擬精度。在解決計算精度和運算效率這一矛盾方面，Dablain(1986)利用Taylor級數(shù)給出了聲波方程高階差分格式，成功消除了空間差分頻散對數(shù)值模擬的影響，極大的提高了有限差分數(shù)值模擬的精度，同時還分析了時間和空間離散的頻散特征。Crase (1990)利用該方法推導了二階彈性波方程高階差分形式，并給出了任意階的高階交錯網(wǎng)格差分方法

16、。前面講到的差分方法它們的網(wǎng)格步長都是固定的，Jastram和Tessmer(1994)提出了垂直間距可變網(wǎng)格差分方法，并將此方法應用于彈性介質(zhì)地震波的模擬，Graves (1996)將交錯網(wǎng)格二維模型模擬推廣到三維模型的模擬，該方法模擬能力有了明顯提高。Pitarka（1999)在三維交錯網(wǎng)格基礎(chǔ)上發(fā)展了矩形非規(guī)則交錯網(wǎng)格，并進行了三維彈性波的模擬；在有限差分算法中除了對網(wǎng)格步長的控制外，還可以通過變時間步長來減少計算量提高計算效率，Tessmer(2000)提出了二階運動方程的變時間算法；為了進一步模擬地震波在非完全彈性的實際地層中的傳播，Carcione(1988)等利用交錯網(wǎng)格模擬了粘

17、滯聲波在地層中的傳播；Saenge和Bohlen(2004)提出了旋轉(zhuǎn)交錯網(wǎng)格有限差分方法，進一步提高了交錯網(wǎng)格對復雜介質(zhì)模擬的精確度，并將旋轉(zhuǎn)交錯網(wǎng)格應用于三維各向異性介質(zhì)中的研究，carcione和Helle(1999)提出孔隙介質(zhì)中粘彈性波交錯網(wǎng)格有限差分數(shù)值模擬。(二)國內(nèi)研究現(xiàn)狀 何樵登(1996等對橫向各向同性介質(zhì)中地震波數(shù)值模擬方法進行了系統(tǒng)介紹；董良國等(2000)詳細闡述了一階速度應力方程交錯網(wǎng)格高階有限差分的實現(xiàn)方法，從理論上證明了交錯網(wǎng)格具有較高的計算精度和計算效率；裴正林等(2005)進一步利用任意偶數(shù)階精度交錯網(wǎng)格有限差分實現(xiàn)了對不同復雜介質(zhì)的模擬；殷文、印興耀、吳國

18、忱等(2006)提出了高精度頻率域彈性波方程有限差分方法，從模擬的結(jié)果看，該方法能很好地消除頻散和邊界反射；李勝軍、孫成禹等(2007)等在不做插值處理的前提下實現(xiàn)了垂向任意倍數(shù)網(wǎng)格步長比高階有限差分算法，并解決了高速層中含有低速層在模擬時存在的數(shù)值頻散和穩(wěn)定性不能均衡的問題。1.4本文主要研究內(nèi)容2 粘彈性介質(zhì)的基本模型 1.粘彈性的模型表述 介質(zhì)的線粘彈性質(zhì)，主要表現(xiàn)為粘滯性和彈性兩方面，這兩方面的性質(zhì)可以由離散的彈性元件彈簧與粘性原件阻尼器進行描述。彈簧和阻尼器以不同的方式進行組合可以得到具有不同性質(zhì)的粘彈性模型。彈性元件彈簧服從胡克定律 =E或=G 其中為正應力、為剪應力、為正應變、為

19、剪應變；E為拉壓彈性模量，G為剪切彈性模量。我們可以看出應力和應變的關(guān)系與時間無關(guān)，不隨時間的變化而變化，彈性變形和彈性回復都是瞬間完成的(圖 2-4)。 圖 2-3 彈性元件模型 圖 2-4 彈性元件應力、應變 粘性元件即阻尼器，有時也被稱為粘壺(圖 2-5)，它服從牛頓粘性定律： = 或 圖 2-5 粘性元件模型 粘彈介質(zhì)交錯網(wǎng)格有限差分數(shù)值模擬 圖 2-6 粘性元件應力、應變響應 （左邊應變響應、右邊應力響應） 其中， 或為粘性系數(shù)；為應變率 阻尼器的流變特性，我們可以通過觀察等應力和等應變作用下它們的響應來了解在應力 =H(t) 作用下，應變響應為即呈穩(wěn)態(tài)流動(圖2-6)其中 H(t)

20、 是單位階躍函數(shù)， 關(guān)于 t=0 時的數(shù)值，將在后面具體問題中加以說明。 在(t) =H(t)的條件下，由式(22)得應力 =(t) )，式中 (t) 為單位脈沖函數(shù)，它滿足兩個條件： 因此，阻尼器受階躍應變作用時，應力為無限大而后瞬即為零(圖 2-6)。由于不可能產(chǎn)生一數(shù)值為無限大的力，所以實際上不能瞬時地使粘性元件產(chǎn)生有限應變。 2 麥克斯韋爾介質(zhì)模型 麥克斯韋爾（Mexwell）介質(zhì)模型由一個彈簧和一個阻尼器串聯(lián)而成(圖 2-7)。設在應力 作用下，彈簧和阻尼器的應變分別為和，則總應變?yōu)椋?圖 2-7 麥克斯韋爾（Mexwell）介質(zhì)模型 得： 或 式中，和均表示材料常數(shù)。式(2-9)即

21、為麥克斯韋爾（Mexwell）介質(zhì)模型的應力應變微分關(guān)系式。 3 開爾芬介質(zhì)模型 開爾芬（Kelvin）介質(zhì)模型由一個彈簧和一個阻尼器并聯(lián)而成如圖 2-8 所示。兩個元件的應變都等于模型的總應變，而模型的總應力為兩元件應力之和，即: 代人上式，得: 或這就是開爾芬（Kelvin）介質(zhì)模型的本構(gòu)關(guān)系。式中，圖 2-8 開爾芬（Kelvin）介質(zhì)模型 4 標準線性粘彈性介質(zhì)模型 如圖 2-9 所示，標準線性粘彈性介質(zhì)是由一個彈簧與一個開爾芬粘彈性單元體相互串聯(lián)而成的三元固體。 圖 2-9 標準線性粘彈性介質(zhì)模型 對于兩個組成元件，分別有： 和總應變?yōu)? 進行拉普拉斯變換及其逆變換就可得到標準線性粘

22、彈性介質(zhì)微分形式的本構(gòu)方程，如下式所示： 3 方法介紹3.1 Stereo-modeling方法簡介Stereo-modeling方法與傳統(tǒng)有限差分方法的主要區(qū)別在于空間離散形式，時間推進可以與任一時間格式結(jié)合。傳統(tǒng)方法的空間高階偏導數(shù)近似可用下面的公式表示8： （3.1）其中m是被近似的高階導數(shù)的階數(shù)。當應用n+1個節(jié)點近似函數(shù)f(x)的m階導數(shù)時，是在節(jié)點j處的系數(shù)。m階偏導數(shù)的近似，需要的節(jié)點個數(shù)依賴于m。當近似的導數(shù)階數(shù)m從1增加到M時，使用的節(jié)點個數(shù)n從m增加到N個。將這種空間離散格式從一維推廣到二維和三維的過程中，每個維度只用到沿著坐標軸方向的節(jié)點，是簡單的線性推廣。Stereo-

23、modeling方法對于空間高階導數(shù)的近似，可以用下面的公式表示5-7 (3.2) 從這個公式中可以看出，Stereo-modeling方法對空間高階導數(shù)的近似，不僅應用了函數(shù)值，還應用了函數(shù)的一階導數(shù)。所用的節(jié)點個數(shù)與具體的導數(shù)階數(shù)m并非直接相關(guān)，而是取決于算法的整體空間精度。具體關(guān)系是，應用2L+1個節(jié)點時，可以取得4L階精度。在下文中我們將采用具有4階空間精度Stereo-modeling方法進行方程的空間離散。3.2Lax-Wendroff correction方法簡介 為了說明新方法STEM的有效性和高效性，我們選取經(jīng)典傳統(tǒng)方法Lax-Wendroff correction（LWC）

24、與其作對比。更確切地說，是選取LWC方法的空間離散格式進行波動方程的空間域離散，對于LWC方法，這里只列出參與計算的具有二階精度的二階空間偏導數(shù)的離散格式 ： (3.3) 4 粘彈性介質(zhì)中的波場數(shù)值模擬 接下來應用上述兩種方法求解以下的二維聲波方程，來進行波場數(shù)值模擬： （3.4）在開爾芬粘彈性介質(zhì)中，對應規(guī)則為 因此上述二維波動方程變?yōu)?（3.5） 其中， ，Q為品質(zhì)因子將（3.5）式同除，記，則 設，=，即，則， 在該數(shù)值模擬實驗中,網(wǎng)格點數(shù)為201×201, 空間網(wǎng)格步長為x = z = 50 m，震源中心頻率為f0=18 Hz，介質(zhì)中聲波速度為c0=4 km/s，時間步長為t

25、=0.002 s。在時間時的入射角。圖2.1顯示了由公式（2.19）計算得到的相對誤差曲線，其中紅色虛線為AD-LWC格式所得，黑色實線為AD-STEM格式所得?？梢钥吹紸D-STEM格式的誤差小于AD-LWC格式的誤差，盡管兩種方法的理論誤差相同。 這是一個置于計算區(qū)域中心的爆炸源，其主頻率為=20 Hz，聲波波速是c=4 km/s,計算區(qū)域為,空間步長為，網(wǎng)格點個數(shù)為241×241，時間步長為0.002 s。圖2.2a和圖2.2b分別顯示了由AD-STEM和AD-LWC方法得到的在粗網(wǎng)格下的波場快照。從圖中可以看出，AD-STEM方法沒有明顯的可見數(shù)值頻散（圖2.2a），而AD-

26、LWC方法則產(chǎn)生了嚴重的數(shù)值頻散（圖2.2b）。為了更清楚地觀察波場模擬情況，在計算區(qū)域中的（6.5 km，6.5 km）位置放置一個接收點，觀察接收到的波形記錄圖。圖2.3和圖2.4分別顯示了AD-STEM和AD-LWC方法在相同計算條件下，對應于圖2.2a和圖2.2b在相同接收點處的波形記錄圖。從這兩個圖中，可以更清晰的看到，AD-STEM方法獲得的數(shù)值解與解析解吻合較好（解析解采用De hoop11等提出的方法計算得到），而AD-LWC方法計算得到的波形則與解析解相差甚遠。為了消除AD-LWC方法獲得的波場快照中的數(shù)值頻散，我們采用加細網(wǎng)格的方法。在保持空間步長與時間步長的比值不變的前提

27、下，將空間網(wǎng)格步長加細到15 m時，數(shù)值頻散才徹底消失（圖2.2c），此時AD-LWC方法的計算時間為250 s，網(wǎng)格點數(shù)為801×801，每個網(wǎng)格點上存儲8個數(shù)組。而AD-STEM方法在粗網(wǎng)格50 m的條件下已無可見的數(shù)值頻散，用時為133 s，網(wǎng)格點數(shù)為241×241，每個網(wǎng)格點上存儲24個數(shù)組?？梢娫跓o可見數(shù)值頻散的條件下，AD-STEM方法的計算速度是AD-LWC方法的1.88倍，而存儲量僅為AD-LWC的72%，具體數(shù)據(jù)參見表一。表一：Stems方法與LWC方法計算效率比較數(shù)值方法網(wǎng)格步長（m）CPU TIME (s)存儲量計算效率AD-STEM50133241&

28、#215;241×241.88AD-LWC15250801×801×80.720 距離（km） 120 深度（km） 12(a) AD-STEM（m）0 距離（km） 120 距離（km） 120 深度（km） 120 深度（km） 12 (b) AD-LWC（m） (c) AD-LWC（m）圖2.2 均勻各向同性介質(zhì)中AD-STEM與AD-LWC方法所得到的波場快照圖2.3 AD-STEM方法得到的數(shù)值解與解析解比較圖2.4 AD-LWC方法得到的數(shù)值解與解析解比較2.3.2 波場數(shù)值模擬本節(jié)將在較復雜的均勻橫向各向同性介質(zhì)（兩層模型）中考察比較兩種方法的數(shù)值性

29、質(zhì)。圖2.5為兩層模型的速度分布示意圖。如圖所示，聲波在上層介質(zhì)(km)和下層介質(zhì)（）中的速度分別為2.4 km/s和5.0 km/s。節(jié)點數(shù)為501×501，模型規(guī)模是20 km×20 km，空間步長為40 m，時間步長為0.001 s。震源仍采用上一節(jié)中用到的震源子波，其中主頻率=15 Hz，震源放置在（10 km，11.2 km）處，接收器R放置在（10 km，10 km）處（分別如圖2.5的五角星和R）。圖2.6為在1.8 s時由AD-STEM（圖2.6（a）和Ad-LWC（圖2.6（b）分別得到的波場快照。從圖2.6可以明顯看出，在相同的空間步長和時間步長的條件下

30、，AD-STEM可以有效的壓制數(shù)值頻散，而傳統(tǒng)AD-LWC方法則有明顯的數(shù)值頻散。為進一步觀察數(shù)值頻散情況，繪制接收點R處的波形記錄圖，且與利用王美霞等13提出的計算方法所得到的解析解進行比較。圖2.7（a）和（b）分別為AD-STEM方法和AD-LWC方法對于接收器R所接收到的波形圖。這兩幅波形對比圖顯示，AD-STEM方法在粗網(wǎng)格（50 m）的情形下，獲得的波形記錄的兩個相位與解析解匹配的都較好，誤差很小，而AD-LWC方法與解析解吻合不好，誤差很大。該結(jié)果與圖2.6的波場快照一致。 圖2.5 雙層模型速度分布示意圖0 距離（km） 20 0 距離（km） 20 (a) AD-STEM (

31、b) Ad-LWC0 深度（km） 20 圖2.6 （a）AD-STEM方法和（b）Ad-LWC方法得到的雙層模型波場快照圖2.7 （a）AD-STEM和（b）AD-LWC方法數(shù)值解與解析解比較3 結(jié)論本文提出了一種新的數(shù)值離散方法，AD-STEM方法。區(qū)別于傳統(tǒng)的數(shù)值離散方法，AD-STEM方法在空間上采取Stereo-modeling方法，時間上采用3階Adams格式，從而得到了一個具有高精度、高效率、低數(shù)值頻散的新方法，對于該方法，可獲得以下結(jié)論：（1）數(shù)值誤差分析顯示，雖然理論誤差相同，但AD-STEM方法比AD-LWC方法的數(shù)值精度更高；（2）計算效率比較數(shù)值實驗顯示，AD-STEM

32、在粗空間網(wǎng)格步長（50 m）下，已經(jīng)能較好的壓制數(shù)值頻散，在波形記錄圖中與解析解也匹配較好。而同等條件下的AD-LWC方法數(shù)值頻散十分明顯，波形記錄圖與解析解匹配不好。在無可見數(shù)值頻散的條件下，AD-STEM方法的計算速度約為AD-LWC的1.88倍，而存儲量為其0.72%。（3）較復雜的雙層模型波場模擬顯示，在同等計算條件下，AD-STEM方法可以有效地壓制數(shù)值頻散，獲得清晰準確的波場快照，波形與解析解吻合，而AD-LWC方法產(chǎn)生較嚴重的數(shù)值頻散，波形嚴重震蕩，與解析解不匹配。綜上，本文發(fā)展的AD-STEM方法是一種數(shù)值性質(zhì)優(yōu)良且十分實用的數(shù)值離散方法，尤其可以利用粗網(wǎng)格進行大規(guī)模波場模擬。

33、參考文獻1 杜世通. 地震波動力學M. 石油大學出版社，1996.2 楊桂通，張善元. 彈性動力學M. 中國鐵道出版社，1988. 3 馮英杰，楊長春，吳萍. 地震波有限差分綜述J. 地球物理學進展，2007.4 佘德平. 波場數(shù)值模擬技術(shù)J. 勘探地球物理進展，2004，27(1)：1621.5 Yang D，Lu M，Wu R，et al. An optimal nearly analytic discrete method for 2D acoustic and elastic wave equations J. Bulletin of the Seismological Society

34、 of America，2004，94：1982 1992. 6 Yang D，Peng J，Lu M，et al. Optimal nearly analytic discrete approximation to the scalar wave equation J. Bulletin of the Seismological Society of America，2006，96：11141130. 7 Yang D，Teng J，Zhang Z，et al. A nearly analytic discrete method for acoustic and elastic wave e

35、quations in anisotropic media J. Bulletin of the Seismological Society of America，2003，93：882 890. 8 Fornberg B. Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of computation J. 1988，51：699 706. 9 牛濱華，孫春巖. 半空間介質(zhì)與地震波傳播M. 北京：石油工業(yè)出版社，2002.10 李靜爽. 基于近似解析離散化算子的逆時偏移方法及其應

36、用D. 北京：清華大學數(shù)學科學系，2014.11 Dablain M. A. The application of high-order differencing to scalar wave equation J. Geophysics，1986，51：54 66.12 De Hoop A. T. A modification of Cagniards method for solving seismic pulse problems J. Appl. Sci. Res. B8，1960：349356.13 王美霞，楊頂輝，宋國杰. 二維SH波方程的半解析解及其數(shù)值模擬J. 地球物理學報，2

37、012，55(3)：914 924.致 謝畢業(yè)論文的完成，代表著四年的大學生活劃上了句號。首先，我要感謝我的論文指導老師，李靜爽老師。無論從前期的準備工作，或是論文整體結(jié)構(gòu)的設計，乃至論文具體數(shù)學理論的講解和程序的調(diào)試，李老師無時無刻不在給予我最大的幫助。甚至在我周中由于實習的原因而無法同老師探討論文問題時，李老師犧牲掉周末的休息時間來幫助我解決問題，修改論文。可以說沒有李老師的幫助，就不會有這篇論文的誕生。其次，我要感謝大學四年時間里傳授我知識的每一位老師和輔導員。是你們不辭辛苦、嘔心瀝血的諄諄教誨才讓我在大學的時光里開拓了眼界，增長了學識。我還要感謝大學四年里身邊的每一位同學，每一個朋友。

38、是你們無私的關(guān)懷與幫助讓我在最美好的年華里收獲了最純真的友誼。最后，我要感謝我的父母和家人。在我大學四年的求學生涯中一直給予我物質(zhì)和精神上的幫助，讓我無憂無慮的專心求學。希望通過這篇論文，能作為我大學四年學習生活的一個總結(jié)，能當作對老師，對家人，對自己的一個交代。由于我的學術(shù)水平有限，本論文難免有不足之處，學生懇請各位老師同學不吝批評和指正。 附錄1 設備清單本文所進行的一切數(shù)值實驗均在該電腦上完成，具體參數(shù)為：電腦型號：惠普 HP ENVY 4 NOTEBOOK PC 筆記本電腦操作系統(tǒng)：Windows 7 家庭普通版 64位 sp1處理器：英特爾 第三代酷睿 i5-3317U 1.70Hz

39、 雙核 超低電壓處理器內(nèi)存：4GB主硬盤：日立 HTS545050A7E380 (500GB / 5400轉(zhuǎn)/分)附錄2 程序清單算法實現(xiàn)的部分主要Fortran 程序：1. 空間偏導數(shù)的計算：U2x=(2.0/(h*h)*(U(I+1,J)-2.0*U(I,J)+U(I-1,J)-(Ux(I+1,J)-Ux(I-1,J)/(2.0*h) U2z=(2.0/(z*z)*(U(I,J+1)-2.0*U(I,J)+U(I,J-1)-(Uz(I,J+1)-Uz(I,J-1)/(2.0*z) U3x=(15.0/(2.0*h*h*h)*(U(I+1,J)-U(I-1,J)-(3.0/(2.0*h*h)

40、*(Ux(I+1,J)+8.0*Ux(I,J)+Ux(I-1,J) U2xz=(1.0/(2.0*h*z)*(-Ux(I+1,J+1)-Ux(I-1,J-1)+Ux(I+1,J)+Ux(I-1,J)-2.0*Ux(I,J+1)+4.0*Ux(I,J)-2.0*U x(I,J-1)+(1.0/(h*h)*(Uz(I+1,J)-2.0*Uz(I,J)+Uz(I-1,J)+(1.0/(4.0*h*h*z)*(5.0*U(I+1,J+1) -5.0*U(I-1,J-1)+U(I+1,J-1)-U(I-1,J+1)-6.0*U(I+1,J)+6.0*U(I-1,J)-4.0*U(I,J+1)+4.0*U(I,J-1) Ux2z=(1.0/(2.0*h*z)*(-Uz(I+1,J+1)-Uz(I-1,J-1)+Uz(I,J+1)+Uz(I,J-1)-2.0*Uz(I+1,J)+4.0*Uz(I,J)-2.0*U z(I-1,J)+(1.0/(z*z)*(Ux(I,J+1)-2.0*Ux(I,J)+Ux(I,J-1)+(1.0/(4.0*h*z*z)*(5.0*U(I+1,J+1) -5.0*U(I-1,J-1)-U(I+1,J-1)+U