高精度空間和時(shí)間任意倍數(shù)可變網(wǎng)格有限差分正演方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及地震波正演模擬技術(shù)，特別是涉及到一種高精度空間和時(shí)間任意倍數(shù) 可變網(wǎng)格有限差分正演方法。
【背景技術(shù)】
[0002] 地震波正演模擬技術(shù)是國(guó)內(nèi)外地球物理學(xué)界的一個(gè)重要研究領(lǐng)域，而基于波動(dòng)理 論的計(jì)算方法最具有應(yīng)用遠(yuǎn)景，因?yàn)樗軌虮３值卣鸩ǖ倪\(yùn)動(dòng)學(xué)與動(dòng)力學(xué)特征，可W達(dá)到 精確模擬地震波傳播特性的目的。其中，有限差分法是最常用的地震波正演數(shù)值模擬方法 之一。它具有概念和列出計(jì)算格式簡(jiǎn)單，對(duì)均勻介質(zhì)、非均勻介質(zhì)、各向異性介質(zhì)都適用的 優(yōu)點(diǎn)，因此在地震學(xué)界受到廣泛重視。
[0003] 有限差分法明顯的缺點(diǎn)是難W克服頻散效應(yīng)，要解決頻散效應(yīng)，必須加密計(jì)算網(wǎng) 格，但送會(huì)導(dǎo)致計(jì)算量的增加，效率下降，并且由于只能給出網(wǎng)格點(diǎn)上的解，模擬精度受到 網(wǎng)格大小等的限制，另外當(dāng)?shù)乇砥鸱^大或地質(zhì)結(jié)構(gòu)復(fù)雜時(shí)，由于傳統(tǒng)的有限差分算法采 用了固定的網(wǎng)格步長(zhǎng)，使得該方法在一些實(shí)例中的有效性受到限制。
[0004] 李勝軍等（《工程地球物理學(xué)報(bào)》Vol. 4,No. 3, 2007, 207-212頁(yè)）實(shí)現(xiàn)在一個(gè)方向 上進(jìn)行變網(wǎng)格，此時(shí)變網(wǎng)格不需要進(jìn)行插值計(jì)算，在一定程度上保證了算法精度，但是方法 仍然有很大的局限性，只能在整個(gè)模型范圍進(jìn)行變網(wǎng)格計(jì)算。為了提高方法的靈活性與模 擬精度，趙海波（《科學(xué)通報(bào)》Vol. 52,No. 12, 2007, 1387-1395頁(yè)）采用了 優(yōu)化變網(wǎng) 格優(yōu)化算法，提高了算子精度；朱生旺等（《石油地球物理勘探》Vol. 42,No. 6, 2007, 634-639 頁(yè)）將網(wǎng)格變化方向推廣至二維，并采用了空間高階差分，提高了空間變網(wǎng)格的靈活性與 模擬精度。
[0005] 但上述可變網(wǎng)格的思路僅在空間網(wǎng)格上進(jìn)行可變處理，對(duì)于提高模擬的效率存在 一定的局限性。因?yàn)樵诮橘|(zhì)變化劇烈的區(qū)域采用空間小網(wǎng)格離散方式要求采用較小的時(shí)間 步長(zhǎng)W滿(mǎn)足計(jì)算的穩(wěn)定性。而如果使用傳統(tǒng)的統(tǒng)一時(shí)間步長(zhǎng)，在介質(zhì)緩變區(qū)域就存在時(shí)間 上過(guò)采樣的問(wèn)題，會(huì)大幅度降低模擬效率。時(shí)間變步長(zhǎng)問(wèn)題仍待解決。為此我們發(fā)明了一 種新的高精度空間和時(shí)間任意倍數(shù)可變網(wǎng)格有限差分正演方法，解決了W上技術(shù)問(wèn)題。

【發(fā)明內(nèi)容】

[0006] 本發(fā)明的目的是提供一種將空間和時(shí)間可變網(wǎng)格結(jié)合起來(lái)，實(shí)現(xiàn)聲波的高精度空 間和時(shí)間任意倍數(shù)可變網(wǎng)格有限差分正演方法。
[0007] 本發(fā)明的目的可通過(guò)如下技術(shù)措施來(lái)實(shí)現(xiàn)：高精度空間和時(shí)間任意倍數(shù)可變網(wǎng)格 有限差分正演方法，該高精度空間和時(shí)間任意倍數(shù)可變網(wǎng)格有限差分正演方法包括：步驟 1，建立地下介質(zhì)的正演速度模型，根據(jù)模型區(qū)域劃分大網(wǎng)格和小網(wǎng)格，并根據(jù)網(wǎng)格大小，設(shè) 置不同空間網(wǎng)格區(qū)域的正演模擬時(shí)間采樣步長(zhǎng)；步驟2,對(duì)正演速度模型進(jìn)行二維網(wǎng)格離 散化，對(duì)正演速度模型中的聲波波場(chǎng)進(jìn)行二維網(wǎng)格離散化，該聲波波場(chǎng)位于網(wǎng)格節(jié)點(diǎn)上；步 驟3,對(duì)最佳匹配層邊界條件網(wǎng)格離散化；W及步驟4,通過(guò)聲波波動(dòng)方程進(jìn)行時(shí)域有限差 分正演模擬，其中時(shí)間采樣步長(zhǎng)在不同大小網(wǎng)格為變步長(zhǎng)。
[0008] 本發(fā)明的目的還可通過(guò)如下技術(shù)措施來(lái)實(shí)現(xiàn)：
[0009] 在步驟2中，聲波波場(chǎng)二維網(wǎng)格離散化公式如下所示：
[0010]
[0011] 其中，U是聲波波場(chǎng)；i和k是二維網(wǎng)格序號(hào)；V是速度；Δt是時(shí)間義樣步長(zhǎng)；ΔX 是X方向網(wǎng)格大小；ΔZ是Z方向網(wǎng)格大??；Μ是笠間Μ階精度；Cm是差分系數(shù)；η是時(shí)間切 片。
[0012] 3.根據(jù)權(quán)利要求1所述的高精度空間和時(shí)間任意倍數(shù)可變網(wǎng)格有限差分正演方 法，其特征在于，在步驟2中，對(duì)于空間過(guò)度網(wǎng)格的聲波波場(chǎng)二維網(wǎng)格離散化如下所示：
[0013] 對(duì)大網(wǎng)格內(nèi)過(guò)渡區(qū)各點(diǎn)，差分格式為：
，：
[0017] 其中，U是聲波波場(chǎng)；ΔΖ是Ζ方向網(wǎng)格大??；Μ是笠間Μ階精度；Cm是差分系數(shù)；η 是時(shí)間切片；Ν是空間網(wǎng)格倍數(shù)。
[0018] 在步驟3中，最佳匹配層邊界條件網(wǎng)格離散化公式如下所示：
[0019]
[0020] 其中，U是聲波波場(chǎng)；i和k是二維網(wǎng)格序號(hào)；At是時(shí)間采樣步長(zhǎng)；Λχ是X方向 網(wǎng)格大小；ΔΖ是Ζ方向網(wǎng)格大??；Μ是笠間Μ階精度；Cm是差分系數(shù)；η是時(shí)間切片，d(x) 和d(z)分別為與X和z方向?qū)?shù)有關(guān)的最佳匹配層邊界衰減項(xiàng)。
[0021] 在步驟4中，改變時(shí)間采樣步長(zhǎng)時(shí)在大時(shí)間采樣步長(zhǎng)區(qū)設(shè)置一個(gè)過(guò)渡帶，利用大 時(shí)間采樣步長(zhǎng)求出過(guò)渡帶區(qū)域中下一時(shí)刻的波場(chǎng)值，采用插值方法得到小時(shí)間采樣步長(zhǎng)計(jì) 算所需的各時(shí)刻波場(chǎng)值。
[0022] 在步驟4中，時(shí)間采樣步長(zhǎng)為任意倍數(shù)，時(shí)間采樣步長(zhǎng)插值方法為使用拉格朗日 算子同時(shí)對(duì)多點(diǎn)時(shí)間采樣步長(zhǎng)進(jìn)行高階非線性插值。
[0023] 本發(fā)明中的高精度空間和時(shí)間任意倍數(shù)可變網(wǎng)格有限差分正演方法，是聲波有限 差分正演方法，能在空間上任意倍數(shù)變網(wǎng)格，時(shí)間上任意倍數(shù)變步長(zhǎng)，充分適應(yīng)各類(lèi)復(fù)雜地 下條件的正演模擬計(jì)算。比如，地表低速層模型、裂縫介質(zhì)模型、生物礁模型W及溶洞模型 等。對(duì)正演模型中復(fù)雜構(gòu)造區(qū)域使用小網(wǎng)格，其他區(qū)域使用大網(wǎng)格，從而對(duì)模型分區(qū)域計(jì) 算，同時(shí)在計(jì)算小網(wǎng)格區(qū)域時(shí)加密時(shí)間片循環(huán)，即在大網(wǎng)格區(qū)域采用大的時(shí)間步長(zhǎng)，小網(wǎng)格 區(qū)域采用小時(shí)間步長(zhǎng)。與常規(guī)有限差分方法對(duì)比，能顯著提高計(jì)算效率和計(jì)算精度。本發(fā) 明與其它空間變網(wǎng)格有限差分方法對(duì)比，能在加密網(wǎng)格區(qū)域的正演時(shí)間上變步長(zhǎng)，充分適 應(yīng)各類(lèi)復(fù)雜地下條件的正演模擬計(jì)算。本發(fā)明能有效壓制地表低速層造成的正演頻散，提 高整個(gè)剖面的信噪比。
【附圖說(shuō)明】
[0024] 圖1為高精度空間和時(shí)間任意倍數(shù)可變網(wǎng)格有限差分正演方法的流程圖；
[00巧]圖2為多級(jí)可變網(wǎng)格模型剖分示意圖；
[0026] 圖3為大網(wǎng)格內(nèi)過(guò)渡區(qū)差分格式示意圖；
[0027] 圖4為小網(wǎng)格內(nèi)過(guò)渡區(qū)差分格式示意圖；
[0028] 圖5為完全匹配層吸收邊界加載示意圖；
[0029] 圖6為變時(shí)間步長(zhǎng)過(guò)渡區(qū)處理方式示意圖；
[0030] 圖7為均勻介質(zhì)模型波場(chǎng)快照；
[0031] 圖8為裂縫介質(zhì)模型示意圖；
[0032] 圖9為裂縫介質(zhì)模型單炮記錄；
[0033] 圖10為生物礁模型示意圖；
[0034] 圖11為生物礁模型單炮記錄；
[0035] 圖12為準(zhǔn)北山前帶地區(qū)模型示意圖；
[0036] 圖13為準(zhǔn)北山前帶模型變網(wǎng)格區(qū)域；
[0037] 圖14為準(zhǔn)北山前帶模型常規(guī)網(wǎng)格單炮記錄；
[0038] 圖15為準(zhǔn)北山前帶模型變網(wǎng)格單炮記錄。
【具體實(shí)施方式】
[0039] 為使本發(fā)明的上述和其他目的、特征和優(yōu)點(diǎn)能更明顯易懂，下文特舉出較佳實(shí)施 例，并配合所附圖式，作詳細(xì)說(shuō)明如下。
[0040] 圖1為本發(fā)明的高精度空間和時(shí)間任意倍數(shù)可變網(wǎng)格有限差分正演方法的流程 圖。其具體實(shí)現(xiàn)流程如下：
[0041] 在步驟101，建立地下介質(zhì)的正演速度模型，根據(jù)模型區(qū)域劃分大網(wǎng)格和小網(wǎng)格。 并根據(jù)網(wǎng)格大小，設(shè)置不同空間網(wǎng)格區(qū)域的正演模擬時(shí)間采樣步長(zhǎng)。
[004引如圖2所示，本發(fā)明實(shí)現(xiàn)多級(jí)可變網(wǎng)格模型剖分。將模型中需要利用小網(wǎng)格計(jì)算 的區(qū)域用內(nèi)外兩層小網(wǎng)格進(jìn)行嵌套剖分，通過(guò)兩次變網(wǎng)格達(dá)到最終需要的模擬精度。內(nèi)外 兩層變網(wǎng)格的倍數(shù)均可W為任意常數(shù)。嵌套的多級(jí)可變網(wǎng)格剖分方式可W通過(guò)低倍數(shù)的變 網(wǎng)格算法來(lái)實(shí)現(xiàn)高倍的網(wǎng)格步長(zhǎng)變化。對(duì)于計(jì)算精度要求高，復(fù)雜地質(zhì)體尺度小的情況，送 種方法更為靈活，同時(shí)相對(duì)于常規(guī)變網(wǎng)格算法，送種方法也可W相應(yīng)節(jié)約計(jì)算內(nèi)存。流程進(jìn) 入到步驟102。
[0043] 在步驟102,對(duì)正演模型進(jìn)行二維網(wǎng)格離散化，對(duì)正演模型中的聲波波場(chǎng)進(jìn)行二維 網(wǎng)格離散化，所述聲波波場(chǎng)位于網(wǎng)格節(jié)點(diǎn)上。其中聲波差分離散化公式為：
[0044]
[004引其中，U是聲波波場(chǎng)；i和k是二維網(wǎng)格序號(hào)；V是速度；Λt是時(shí)間采樣步長(zhǎng)；ΛX是X方向網(wǎng)格大小；ΔΖ是Ζ方向網(wǎng)格大小；Μ是笠間Μ階精度；Cm是差分系數(shù)；η是時(shí)間切 片。
[0046] 對(duì)于空間過(guò)度網(wǎng)格的差分格式需要修正，對(duì)大網(wǎng)格內(nèi)過(guò)渡區(qū)各點(diǎn)，差分格式為：
[0047]
[004引其中，Ν是空間網(wǎng)格倍數(shù)。
[0049] 圖3是大網(wǎng)格內(nèi)過(guò)渡區(qū)差分格式示意圖。大網(wǎng)格區(qū)域中的過(guò)渡區(qū)5個(gè)點(diǎn)，其差分 計(jì)算取點(diǎn)在小網(wǎng)格區(qū)域內(nèi)并非順序取點(diǎn)，而是維持著關(guān)于中必點(diǎn)對(duì)稱(chēng)的格式，仍使用大網(wǎng) 格的差分格式和差分系數(shù)。
[0050] 小網(wǎng)格內(nèi)過(guò)渡區(qū)各點(diǎn)差分格式為：
[0051]
[005引其中，Ν是空間網(wǎng)格倍數(shù)。
[0053] 圖4是小網(wǎng)格內(nèi)過(guò)渡區(qū)差分格式示意圖。小網(wǎng)格區(qū)域中的過(guò)渡區(qū)5個(gè)點(diǎn)，當(dāng)差分 計(jì)算需要用到大網(wǎng)格點(diǎn)時(shí)，小網(wǎng)格內(nèi)使用與之對(duì)稱(chēng)的點(diǎn)參與計(jì)算。此時(shí)步長(zhǎng)值發(fā)生了改變， 需要修正差分格式。流程進(jìn)入到步驟103。
[0054] 在步驟103,對(duì)最佳匹配層邊界條件網(wǎng)格離散化，最佳匹配層邊界條件網(wǎng)格離散化 公式如下所示：
[00巧]
[0056] 其中，d(x)和d(z)分別為與X和Z方向?qū)?shù)有關(guān)的最佳匹配層邊界衰減項(xiàng)。
[0057] 最佳匹配層的構(gòu)造如圖5所示。最佳匹配層邊界位于模型計(jì)算區(qū)域外，該邊界被 分成幾個(gè)不同區(qū)域，每個(gè)區(qū)域的對(duì)波場(chǎng)吸收的作用不同，因此衰減因子取值不同。B1和B2 所代表的最佳匹配層吸收區(qū)域中，該區(qū)域只對(duì)沿Z軸方向傳播的波進(jìn)行衰減；B3和B4所代 表的PML吸收區(qū)域中，該區(qū)域?qū)ρ豖軸方向傳播的波進(jìn)行衰減；在C1，C2,C3和C4所代表 的最佳匹配層吸收區(qū)域中，角點(diǎn)對(duì)沿X軸方向和Z軸方向傳播的波均進(jìn)行衰減。流程進(jìn)入 到步驟104。
[0058] 在步驟104,通過(guò)聲波波動(dòng)方程進(jìn)