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

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

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