任意對數(shù)階精度交錯網(wǎng)格差分格式的討論

任意對數(shù)階精度交錯網(wǎng)格差分格式的討論

地震波場值模擬是地震學(xué)和地震勘探的重要基礎(chǔ)，在地震學(xué)和地震勘探的所有過程中起著重要作用（pedrin等人，2004）。通過正演模擬技術(shù),可以模擬地震波在各種地質(zhì)結(jié)構(gòu)中的傳播規(guī)律,可為地震資料處理解釋中各種方法的應(yīng)用提供標(biāo)準(zhǔn)的先驗(yàn)數(shù)據(jù)體,檢測方法的實(shí)用性(趙海波等,2011;張華等,2012)。有限差分方法因其計(jì)算效率高、占用計(jì)算內(nèi)存小等特點(diǎn),成為當(dāng)前最常用的數(shù)值模擬方法(張英杰等,2007;呂玉增等,2011)。交錯網(wǎng)格技術(shù)相對于規(guī)則網(wǎng)格在復(fù)雜模型模擬和模擬精度等方面都有較大的優(yōu)勢(李紅星等,2008;李國平等2010),交錯網(wǎng)格引入的頻散明顯低于規(guī)則網(wǎng)格,因此地震波場數(shù)值模擬中交錯網(wǎng)格有限差分法得到廣泛應(yīng)用。為得到高精度的正演模擬結(jié)果,本文基于一階應(yīng)力-速度聲波方程組,運(yùn)用任意偶數(shù)階精度的交錯網(wǎng)格有限差分格式進(jìn)行了各向同性介質(zhì)模型中的正演模擬。1波動方程的有限差分格式為1.1電壓速度的初級聲波方程根據(jù)彈性力學(xué)的分析,可得到動力學(xué)方程組和應(yīng)變應(yīng)力方程組。二維情況下一階應(yīng)力-速度聲波方程組(張伯軍,2010)如下:式中K=ρv1.2交錯網(wǎng)格設(shè)置一階應(yīng)力-速度聲波方程組各個物性參量在交錯網(wǎng)格中的位置設(shè)置如下:空間網(wǎng)格中,參數(shù)P放置在規(guī)則的i,j節(jié)點(diǎn)上。而參數(shù)v時間網(wǎng)格中,參數(shù)P在規(guī)則的k節(jié)點(diǎn)上,v1.3)陳可洋,陳可洋本文從Taylor級數(shù)展開式出發(fā),實(shí)現(xiàn)了任意偶數(shù)階精度交錯網(wǎng)格下的差分格式算法(劉洋等,1998;陳可洋,2009)。設(shè)函數(shù)U(x)為存在2L+1階導(dǎo)數(shù)的函數(shù),其在空間i+1/2半節(jié)點(diǎn)處,分別對i與i+1點(diǎn)進(jìn)行泰勒展開處理,得到U將式(2)與式(3)相減得:類似上述推導(dǎo)過程,可得任意U將上式展開成線性方程式(7):解方程式得:式中a1.4時間二階空間任意級數(shù)階交錯網(wǎng)格差分算法實(shí)現(xiàn)由式(1)給出的二維情況下一階應(yīng)力-速度方程組,可知需對v由式(1)結(jié)合式(9)及式(10)給出的一階導(dǎo)數(shù)差分計(jì)算公式,便可實(shí)現(xiàn)時間二階空間任意偶數(shù)階精度交錯網(wǎng)格下的聲波方程差分格式(11)。2在數(shù)值模擬中，一些重要問題2.1震源函數(shù)的選取震源函數(shù)有多種取法,如高斯子波、雷克子波等。在波場數(shù)值模擬中,震源函數(shù)的選取對模擬結(jié)果有重要的影響(張海燕等,2007;董良國等,2004)。本文采用雷克子波作為震源函數(shù)。2.2空間域網(wǎng)格步長與穩(wěn)定性條件對文中一階應(yīng)力-速度聲波方程組任意偶數(shù)階精度交錯網(wǎng)格有限差分方程,通過平面諧潑的分析(裴正林等,2003),得出其數(shù)值解的穩(wěn)定性條件:其中,Δt為時間域網(wǎng)格步長,Δx、Δz表示空間域網(wǎng)格步長,v2.3衰減因子的選擇基于在模擬區(qū)外圍附加吸收層的思想,Cerjan等(1985)提出了衰減邊界條件。在模擬區(qū)域內(nèi),波按照正常的波動傳播;當(dāng)波到達(dá)外圍衰減區(qū)時,波場乘以一個衰減因子G,使其按指數(shù)級逐漸衰減。其中:式中N為衰減區(qū)總網(wǎng)格數(shù),i為衰減區(qū)的網(wǎng)格編號(1ue025iue025N),a為衰減系數(shù)。圖3為衰減邊界示意圖。李信富等(2007)研究表明,衰減因子參數(shù)的選擇應(yīng)采取折中方法,即在取某一較小衰減系數(shù)的前提下,設(shè)置較少的衰減帶網(wǎng)格數(shù)使邊界的能量反射最小。經(jīng)對模擬速度和精度的綜合考慮,本文a值選用0.15,N值為50。那么在衰減區(qū)波場計(jì)算公式為:3選擇測試和限制效果的比較3.1差分精度選取測試?yán)碚撋险f,差分算子的階數(shù)越高,波場模擬結(jié)果也越精確。但階數(shù)過高,會造成一定的計(jì)算困難(左瑩,2007)。那么在效果與效率平衡的基礎(chǔ)上,選取合適的差分精度是很有必要的。為確定合適的差分精度,本文進(jìn)行了差分精度選取測試。測試模型為均勻各向同性介質(zhì)模型,模型介質(zhì)參數(shù):v=1000m/s,ρ=1g/m通過對比分析,易得出隨著空間差分階數(shù)的提高,模擬精度也逐漸增高?？臻g二階精度時,存在嚴(yán)重的頻散現(xiàn)象;空間四階精度時,頻散現(xiàn)象明顯減弱?？臻g八階精度時,幾乎無頻散現(xiàn)象出現(xiàn);階數(shù)再增加,模擬精度變化不明顯。因此,為得到高精確的正演模擬結(jié)果,空間差分精度不得低于八階。本文的波場模擬都是采用空間八階精度進(jìn)行的。測試模型為均勻各向同性介質(zhì)模型。震源放置在模型中心點(diǎn)處,震源函數(shù)頻率為25Hz。選取640ms時未使用衰減邊界條件和使用衰減邊界條件的波場快照進(jìn)行邊界效果對比(圖5)。4模型建模分別對不同類型的模型進(jìn)行了正演模擬。模擬過程均采用時間二階空間八階差分精度。4.1震源模型正演模擬如圖6a所示水平層狀介質(zhì)模型,模型大小為1000×1000m,水平分界面在埋深500m處。上層介質(zhì)參數(shù):v震源位于水平方向500m,埋深10m處,震源函數(shù)頻率為25Hz。對模型進(jìn)行正演模擬,得到不同時刻的波場快照。圖6b和圖6c分別為水平層狀介質(zhì)模型560ms和600ms時刻的波場快照。波場快照中,界面反射波和透射波清晰可見,邊界反射被有效吸收,未出現(xiàn)頻散現(xiàn)象。圖6d為水平層狀介質(zhì)模型的共炮點(diǎn)記錄,直達(dá)波及水平分界面z=500m處的反射波同相軸清晰。4.2震源表面組合介質(zhì)模型依舊采用201×201的網(wǎng)格剖分,空間網(wǎng)格步長為5m,時間網(wǎng)格步長0.4ms,第一層介質(zhì)模型參數(shù):v傾斜介質(zhì)模型中傾斜分界面從左邊界埋深300m處延伸到右邊界埋深550m處(圖7a);水平分界面在埋深z=550m處。震源位于水平方向600m,埋深10m處,震源函數(shù)頻率為25Hz。同樣采用空間八階時間二階差分精度進(jìn)行波場模擬。圖7b和圖7c分別為傾斜介質(zhì)模型420ms和600ms時刻的波場快照,可清晰反映出波在模型中的傳播規(guī)律。邊界反射被有效吸收,但存在輕微頻散現(xiàn)象。圖7d為傾斜介質(zhì)模型的共炮點(diǎn)記錄,圖中直達(dá)波同相軸和界面反射波同相軸清晰明顯。邊界反射被有效吸收,整體效果良好。4.3介質(zhì)模型的建立如圖8a所示斷層介質(zhì)模型,斷層分界面中,垂直段在x=600m處,從埋深380m延伸到520m。斷層水平段一部分在埋深380m處,水平坐標(biāo)x=0至x=600之間;水平段另一部分在埋深520m處,水平坐標(biāo)x=600至x=1000m之間。第二、三層介質(zhì)的水平分界面在z=600m處。第一層模型的介質(zhì)參數(shù)為:v震源位于水平方向500m,埋深10m處,震源函數(shù)頻率為25Hz。采用同樣的差分精度對斷層介質(zhì)模型進(jìn)行模擬試算。圖8b和圖8c分別為420ms和600ms時刻的模型波場快照,充分反映出斷層介質(zhì)模型中波的傳播規(guī)律。圖中可清晰分辨出界面反射波和透射波,以及斷層的斷點(diǎn)上產(chǎn)生的繞射波;邊界反射吸收效果良好,未出現(xiàn)頻散現(xiàn)象。圖8d為斷層介質(zhì)模型的共炮點(diǎn)記錄。圖中直達(dá)波、斷層界面的反射波及水平界面反射波同相軸清晰。斷層上下水平段界面對應(yīng)的反射波同相軸末端,都伴隨著明顯的繞射波同相軸,與實(shí)際情況相符合。5數(shù)值模擬效果高偶數(shù)階交錯網(wǎng)格有限差分?jǐn)?shù)值模擬方法可以清晰模擬出波在地下介質(zhì)中的傳播規(guī)律。本文方法模擬精度高,計(jì)算效率快,能有效壓制頻散噪音,邊界反射吸收效果良好,充分證實(shí)了該數(shù)值模擬方法的準(zhǔn)確性和可行性,具有較高的應(yīng)用價值?？梢酝茝V應(yīng)用于三維聲波高精度正演模擬中,幫助我們深入學(xué)習(xí)和研究地震波