一種適用于交錯(cuò)網(wǎng)格有限差分的全吸收pml方法
【專利摘要】本發(fā)明涉及一種適用于交錯(cuò)網(wǎng)格有限差分的全吸收PML方法，綜合PML中的d⊥，C-PML中的α，k，M-PML中的m因素，得到與復(fù)雜的大斜度井聲波測井吻合的數(shù)值模擬結(jié)果；將C-PML當(dāng)中的X、Y、Z向吸收邊界的拉伸因子和M-PML當(dāng)中的X、Y、Z向吸收邊界的拉伸因子分別融合為混合X向、Y向、Z向吸收邊界的拉伸因子：本發(fā)明的有益效果是：通過參數(shù)優(yōu)化，吸收C-PML和M-PML優(yōu)勢，實(shí)現(xiàn)了PML吸收邊界在求解TTI介質(zhì)時(shí)的高效性。在極端入射角入射的情況下計(jì)算精確，對(duì)極端各向異性介質(zhì)也能保持良好的穩(wěn)定性。
【專利說明】-種適用于交錯(cuò)網(wǎng)格有限差分的全吸收PML方法

【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于地球物理（測井）勘探方法，尤其涉及彈性波數(shù)值模擬領(lǐng)域。

【背景技術(shù)】
[0002] 石油是一種重要的戰(zhàn)略資源。從經(jīng)濟(jì)和技術(shù)的角度，無論是對(duì)外合作、通過國際投 標(biāo)購買國外油田，還是在國內(nèi)進(jìn)行勘探尋找新的油氣區(qū)以及對(duì)老油田進(jìn)行重新評(píng)價(jià)，數(shù)值 模擬方法都是一個(gè)重要的關(guān)鍵技術(shù)，由于它可以為實(shí)際數(shù)據(jù)提供理論依據(jù)，所以無論是對(duì) 國民經(jīng)濟(jì)的發(fā)展，還是石油儲(chǔ)備量的合理確定，都有十分重要的意義。
[0003] 聲反射成像測井已經(jīng)成為一種具有潛力的遠(yuǎn)探測方法，無論是對(duì)地質(zhì)導(dǎo)向，還是 對(duì)于大斜度井和水平井環(huán)境下的頁巖氣開采，都起到了重要的作用。在各種常用的聲反射 成像測井技術(shù)的數(shù)值模擬方法中，交錯(cuò)網(wǎng)格有限差分方法是一種相對(duì)高效的方法，并且可 以對(duì)復(fù)雜的三維模型進(jìn)行模擬，更接近真實(shí)情況。我們在研究大斜度井穿越頁巖層（往往 是VTI(垂向橫向各向同性）介質(zhì)）所產(chǎn)生的井孔聲場時(shí)，所遇到的一個(gè)較大問題是，有限 差分網(wǎng)格為正四邊形，在大斜度井中的建模會(huì)產(chǎn)生鋸齒狀網(wǎng)格，導(dǎo)致無法避免的離散誤差。 通過將參考坐標(biāo)系進(jìn)行坐標(biāo)變換，可以將VTI介質(zhì)中的模擬問題轉(zhuǎn)換為TTI(水平橫向各向 同性）介質(zhì)中。


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

[0004] 本發(fā)明的目的是提出一種適用于交錯(cuò)網(wǎng)格有限差分的全吸收PML方法的技術(shù)方 案。同時(shí)具備穩(wěn)定性和準(zhǔn)確性，使PML吸收邊界在求解TTI介質(zhì)時(shí)的實(shí)現(xiàn)高效性。
[0005] 為了實(shí)現(xiàn)上述目的，本發(fā)明的技術(shù)方案是：一種適用于交錯(cuò)網(wǎng)格有限差分的全吸 收PML方法，綜合PML中的d丄，C-PML中的a，k，M-PML中的m因素，得到與復(fù)雜的大斜度 井聲波測井吻合的數(shù)值模擬結(jié)果；
[0006] 將C-PML當(dāng)中的X向吸收邊界的拉伸因子和M-PML當(dāng)中的X向吸收邊界的拉伸因 子融合為一個(gè)混合X向吸收邊界的拉伸因子： ,d, rd, +m, _d_
[0007] ? + ---~ +ιω
[0008] 將C-PML當(dāng)中的Y向吸收邊界的拉伸因子和M-PML當(dāng)中的Y向吸收邊界的拉伸因 子融合為一個(gè)混合Y向吸收邊界的拉伸因子： rn ,d+mvλ?λ +m,Jz ar + no
[0010] 將C-PML當(dāng)中的Z向吸收邊界的拉伸因子和M-PML當(dāng)中的X向吸收邊界的拉伸因 子融合為一個(gè)混合X向吸收邊界的拉伸因子： d, + m_ cl + m_ d y
[0011] Sz =kz+----， a. + ico
[0012] 上述式中：dxSpml在X向的衰減函數(shù)，dySpml在Y向的衰減函數(shù)，dzSpml在 z向的衰減函數(shù)，
[0013] CtxSX向復(fù)頻移函數(shù)，
[0014] Cty為Y向復(fù)頻移函數(shù)，
[0015] azSZ向復(fù)頻移函數(shù)，
[0016] kx為X向抑制倏逝波的函數(shù)，
[0017] ky為Y向抑制倏逝波的函數(shù)，
[0018] kz為Z向抑制倏逝波的函數(shù)，
[0019] mx/y為X向與Y向間的切向衰減控制因子，
[0020] mx/z為X向與Z向間的切向衰減控制因子，
[0021] my/x為Y向與X向間的切向衰減控制因子，
[0022] my/z為Y向與Z向間的切向衰減控制因子，
[0023] mz/x為Z向與X向間的切向衰減控制因子，
[0024] mz/y為Z向與Y向間的切向衰減控制因子，
[0025] ω為角頻率，為頻率的2π倍，
[0026] i為虛數(shù)單位。
[0027] 更進(jìn)一步，所述切向衰減控制因子mx/y、mx/z、my/x、my/z、mz/x、mz/y的取值范圍為 0. 005 ?0. 020。
[0028] 本發(fā)明的有益效果是：通過參數(shù)優(yōu)化，兼具C-PML和M-PML的優(yōu)勢，實(shí)現(xiàn)了PML吸 收邊界在求解TTI介質(zhì)時(shí)的高效性。在極端入射角入射的情況下計(jì)算精確，對(duì)極端各向異 性介質(zhì)也能保持良好的穩(wěn)定性；可準(zhǔn)確地對(duì)TTI介質(zhì)進(jìn)行數(shù)值模擬，獲得復(fù)雜的大斜度井 模型。
[0029] 下面結(jié)合附圖和實(shí)施例對(duì)本發(fā)明作一詳細(xì)描述。

【專利附圖】

【附圖說明】
[0030] 圖1是本發(fā)明一個(gè)各向同性均一模型分別采用C-PML和M-PML的效果圖，上側(cè)為 C-PML的吸收效果，下側(cè)為M-PML的吸收效果；
[0031] 圖2是本發(fā)明一個(gè)極端各向異性均一模型分別采用C-PML和M-PML的效果圖，上 側(cè)為C-PML的吸收效果，下側(cè)為M-PML的吸收效果；
[0032] 圖3是一個(gè)極端TTI各向異性井眼模型圖；
[0033] 圖4是本發(fā)明當(dāng)m= 0. 000 (即只采用C-PML時(shí)）的波場傳播情況；
[0034] 圖5是本發(fā)明當(dāng)m= 0. 005 (即采用全吸收融合PML時(shí)）的波場傳播情況；
[0035] 圖6是當(dāng)模型較大時(shí)的參考解的波形圖；
[0036] 圖7是本發(fā)明當(dāng)m= 0. 000 (即只采用C-PML時(shí)）的波形圖；
[0037] 圖8是本發(fā)明當(dāng)m= 0. 005 (即采用全吸收融合PML時(shí)），吸收參數(shù)恰當(dāng)時(shí)的波形 圖；
[0038] 圖9是本發(fā)明當(dāng)m= 0. 010 (即采用全吸收融合PML時(shí)），吸收參數(shù)不恰當(dāng)時(shí)的波 形圖。

【具體實(shí)施方式】
[0039] 聲反射成像測井已經(jīng)成為一種具有潛力的遠(yuǎn)探測方法，在各種常用的聲反射成像 測井技術(shù)的數(shù)值模擬方法中，交錯(cuò)網(wǎng)格有限差分方法是一種相對(duì)高效的方法，并且很容易 應(yīng)用到研究復(fù)雜的三維情況中。大斜度井穿越頁巖層（往往是VTI介質(zhì)）所產(chǎn)生的井孔聲 場時(shí)，由于有限差分網(wǎng)格為正四邊形，因此大斜度井的建模會(huì)產(chǎn)生鋸齒狀，并產(chǎn)生無法避免 的離散誤差。
[0040] 在本發(fā)明中，我們采用一種方法來避免這種問題，它是通過旋轉(zhuǎn)參考坐標(biāo)系，使得 井軸總是與直角坐標(biāo)系中的其中一個(gè)坐標(biāo)軸保持一致，這樣VTI介質(zhì)將通過Bond變換矩陣 轉(zhuǎn)化為TTI介質(zhì)。于是我們的這種方法將可以廣泛適用于（1)考查VTI對(duì)稱軸與井軸之間 不同角度的影響；（2)考查井軸與井外地層界面不同角度的影響。
[0041] 本發(fā)明主要致力于解決PML(完全匹配吸收層）吸收邊界在求解TTI介質(zhì)時(shí)的高 效性。C-PML(卷積完全匹配吸收層）在極端入射角入射的情況下依舊計(jì)算精確，然而它在 某些極端各向異性介質(zhì)中會(huì)出現(xiàn)數(shù)值不穩(wěn)定。另一方面，M-PML(多軸完全匹配吸收層）即 使對(duì)極端各向異性介質(zhì)也能保持卓越的穩(wěn)定性，但在PML層厚度不夠厚，以及修正因子不 夠優(yōu)化的情況下，卻犧牲了準(zhǔn)確性。而過厚的PML層厚度又大幅降低了計(jì)算效率。為了達(dá) 到穩(wěn)定性和準(zhǔn)確性的共贏，我們引入一種混合型吸收邊界，通過參數(shù)優(yōu)化，使兩種PML的優(yōu) 勢能夠共存。
[0042] 本發(fā)明的特點(diǎn)是同時(shí)反映了PML中的d丄，C-PML中的a，k，M_PML中的d//等因 素的影響，可以最大限度克服穩(wěn)定性和準(zhǔn)確性難以兼顧的難題。同時(shí)這種新型吸收邊界在 吸收層較薄的情況下（10個(gè)網(wǎng)格左右的厚度，其厚度遠(yuǎn)遠(yuǎn)小于波長）依然可以獲得較好的 吸收效果，故計(jì)算效率也優(yōu)于PMUC-PML和M-PML的任何一個(gè)。在這種方法中，傳統(tǒng)的PML 是在此新型PML下的特例，即本方法涵括了PMUC-PML和M-PML，同時(shí)在上述PML之一不能 適用的情況下，本方法依然適用。
[0043] 三維交錯(cuò)網(wǎng)格有限差分方法中的速度-應(yīng)力方程為： IdvtIiv,Idv_
[0044] -脈'τ=--- + --- + --^ ox sy oy OZ
[0045] 其中以sx為例，它為PML吸收邊界的拉伸因子。在最為傳統(tǒng)的PML當(dāng) 中，&二1 + 4，而在C-PML和M-PML當(dāng)中，它分別被修正為A+ 和 ?ω α, +ιω dr +/?7dv +m. _d. sx=\ + ~~~^。同理可以得到Sy和Sz。ιω
[0046] 本發(fā)明的技術(shù)方案是：一種適用于交錯(cuò)網(wǎng)格有限差分的全吸收PML方法。綜合PML 中的d丄，C-PML中的a，k，M-PML中的m因素，得到與復(fù)雜的大斜度井聲波測井吻合的數(shù) 值模擬結(jié)果； ,d,
[0047] 將C-PML當(dāng)中的X向吸收邊界的拉伸因子A=屹+ -:―和M-PML當(dāng)中的X向吸 CXx +10) 收邊界的拉伸因子心=1+< 融合為一個(gè)混合X向吸收邊界的拉伸因子: I(I) d+m,'d'+m!d
[0048] sx=kx+---~; αΛ+ιω d"
[0049] 將C-PML當(dāng)中的Y向吸收邊界的拉伸因子& = & + -=-和M-PML當(dāng)中的Y向吸 I(〇 收邊界的拉伸因子\ =1 +< j-_融合為一個(gè)混合Y向吸收邊界的拉伸因子: 1(0 d+m,d+m,d
[0050] 八=人',+ ---~； CXy +1(0 (L
[0051] 將C-PML當(dāng)中的Z向吸收邊界的拉伸因子\ +-"-和M-PML當(dāng)中的X向吸 α_ +ιω 收邊界的拉伸因子& =1+4+W_ 融合為一個(gè)混合X向吸收邊界的拉伸因子: ico 「n , dz+m:A、+m:'d'
[0052] Sz =kz+-〇 a. +?ω
[0053] 上述式中：dxSPML在X向的衰減函數(shù)，dySPML在Y向的衰減函數(shù)，dzSPML在 Z向的衰減函數(shù)；
[0054]αχ為X向復(fù)頻移函數(shù)，對(duì)應(yīng)于C-PML中的X向復(fù)頻移函數(shù)；
[0055]ay為Y向復(fù)頻移函數(shù)，對(duì)應(yīng)于C-PML中的Y向復(fù)頻移函數(shù)；
[0056]αz為Z向復(fù)頻移函數(shù)，對(duì)應(yīng)于C-PML中的Z向復(fù)頻移函數(shù)；
[0057]kx為X向抑制倏逝波的函數(shù)，對(duì)應(yīng)于C-PML中的X向抑制倏逝波的函數(shù)；
[0058]ky為Y向抑制倏逝波的函數(shù)，對(duì)應(yīng)于C-PML中的Y向抑制倏逝波的函數(shù)；
[0059]kz為Z向抑制倏逝波的函數(shù)，對(duì)應(yīng)于C-PML中的Z向抑制倏逝波的函數(shù)；
[0060]mx/y為X向與Y向間的切向衰減控制因子；對(duì)應(yīng)于M-PML中X向與Y向間的切向衰 減控制因子；
[0061]mx/z為X向與Z向間的切向衰減控制因子；對(duì)應(yīng)于M-PML中X向與Z向間的切向衰 減控制因子；
[0062]my/x為Y向與X向間的切向衰減控制因子；對(duì)應(yīng)于M-PML中Y向與X向間的切向衰 減控制因子；
[0063]my/z為Y向與Z向間的切向衰減控制因子；對(duì)應(yīng)于M-PML中Y向與Z向間的切向衰 減控制因子；
[0064]mz/x為Z向與X向間的切向衰減控制因子；對(duì)應(yīng)于M-PML中Z向與X向間的切向衰 減控制因子；
[0065]mz/y為Z向與Y向間的切向衰減控制因子；對(duì)應(yīng)于M-PML中Z向與Y向間的切向衰 減控制因子；
[0066]ω為角頻率，為頻率的2π倍，
[0067]i為虛數(shù)單位。
[0068] 所述切向衰減控制因子mx/y、mx/z、my/x、my/z、mz/x、mz/y的取值范圍為0· 005?0· 020。
[0069] 本發(fā)明在實(shí)現(xiàn)C-PML的基礎(chǔ)上，將其中的dx修正為dx+mx/ydy+mx/zdz，并同理修改dy 和dz為相應(yīng)的形式，即實(shí)現(xiàn)了兩種方法的有機(jī)融合。
[0070] 大量實(shí)驗(yàn)表明，相比于傳統(tǒng)PML和C-PML，盡管優(yōu)化的α和k能夠在無精度損失 的前提下大幅提高吸收效果，但在TTI介質(zhì)的模擬中仍不足夠精確。于是需要M-PML中的 m(包括mx/y、mx/z、my/x、my/z、mz/x、mz/y)取一個(gè)足夠小的值,經(jīng)過研究與實(shí)驗(yàn)，m值取0. 005到 0. 020,就可以在最小的精度損失中達(dá)到最大的穩(wěn)定性。
[0071] 實(shí)施例一：
[0072] 我們考慮兩個(gè)較為典型的均一模型，如圖1和圖2所示。其中圖1為各向同性介 '4 _ 質(zhì)（其廣義Hooke矩陣為C= 4 的波場快照，此時(shí)C-PML和M-PML均是穩(wěn)定 2 的，C-PML吸收效果較好，而M-PML產(chǎn)生了部分虛假反射。圖2為極端各向異性介質(zhì)（其廣 _ 4 7.5 " 義Hooke矩陣為C= 7.5 20 (/"〃）的波場快照，此時(shí)C-PML是不穩(wěn)定的，M-PML是穩(wěn)定 2 的，C-PML吸收效果較好，而M-PML產(chǎn)生了部分虛假反射。此例子印證了兩種PML吸收邊界 融合的必要性。
[0073] 我們進(jìn)一步考慮一個(gè)典型的聲波測井模型，如圖3所示。模型大小為0. 6mX6.Om， 井內(nèi)流體內(nèi)有一實(shí)心剛性儀器存在，聲源采用的是中心頻率為8kHz的偶極子聲源。井外地 '17.47 11.93 2.14' 層為將一個(gè)VTI地層旋轉(zhuǎn)45°所得，其廣義Hooke矩陣為C= 11.93 17.47 2.14 ￠7/?。另 2.14 2.14 4.91 夕卜，我們在這個(gè)問題中，吸收邊界層只占10個(gè)網(wǎng)格的厚度，使得其厚度遠(yuǎn)遠(yuǎn)小于波長，這對(duì) 我們的數(shù)值模擬結(jié)果是一個(gè)非常大的挑戰(zhàn)。
[0074] 圖4和圖5為每Ims間隔所得到的波場快照，可見當(dāng)m= 0. 000時(shí)（即C-PML情 形），模型中出現(xiàn)了數(shù)值不穩(wěn)定，而當(dāng)m= 0. 005時(shí)，模擬效果較好。
[0075] 我們進(jìn)一步考查模擬波形的結(jié)果。由于當(dāng)TI介質(zhì)的對(duì)稱軸與井軸不一致，難以求 取解析解，因此本專利與參考解進(jìn)行相互驗(yàn)證。參考解是通過將模型設(shè)置的較大，從而使得 PML吸收邊界不影響模擬結(jié)果而求取的。如圖6-9所示。其中當(dāng)m= 0. 000時(shí)，數(shù)值不穩(wěn)定 導(dǎo)致T= 3.Oms時(shí)出現(xiàn)了較大的誤差。當(dāng)m= 0.005時(shí)，兩者吻合較好。而m= 0.010時(shí)， 兩者吻合度略差于m= 0. 005的情況，說明該混合吸收邊界的參數(shù)優(yōu)化對(duì)于數(shù)值模擬的準(zhǔn) 確性是至關(guān)重要的。
[0076] 本發(fā)明的方法已被應(yīng)用于一個(gè)實(shí)際的三維大斜度井問題中，并取得了良好的效 果。
【權(quán)利要求】
1. 一種適用于交錯(cuò)網(wǎng)格有限差分的全吸收PML方法，綜合PML中的di，C-PML中的α， k，M-PML中的m因素，得到與復(fù)雜的大斜度井聲波測井吻合的數(shù)值模擬結(jié)果；其特征在于： 將C-PML當(dāng)中的X向吸收邊界的拉伸因子和M-PML當(dāng)中的X向吸收邊界的拉伸因子融 合為一個(gè)混合X向吸收邊界的拉伸因子： ,+/?Λ Λ7, +/?. β, ? +- ， αχ +ιω 將C-PML當(dāng)中的Υ向吸收邊界的拉伸因子和M-PML當(dāng)中的Υ向吸收邊界的拉伸因子融 合為一個(gè)混合Υ向吸收邊界的拉伸因子： dν + 川、、ci、十 m、Μ. Sy=ky+-~---；~， a κ +ιω 將C-PML當(dāng)中的Z向吸收邊界的拉伸因子和M-PML當(dāng)中的X向吸收邊界的拉伸因子融 合為一個(gè)混合X向吸收邊界的拉伸因子： d_ +m_ +m:l d sz=K+-:-， (X.. +1(0 上述式中：dxSPML在X向的衰減函數(shù)，dySPML在Y向的衰減函數(shù)，dzSPML在Z向 的衰減函數(shù)， α x為X向復(fù)頻移函數(shù)， a y為Y向復(fù)頻移函數(shù)， α z為Z向復(fù)頻移函數(shù)， kx為X向抑制倏逝波的函數(shù)， ky為Y向抑制倏逝波的函數(shù)， kz為Z向抑制倏逝波的函數(shù)， mx/y為X向與Y向間的切向衰減控制因子， mx/z為X向與Z向間的切向衰減控制因子， my/x為Y向與X向間的切向衰減控制因子， my/z為Y向與Z向間的切向衰減控制因子， mz/x為Z向與X向間的切向衰減控制因子， mz/y為Z向與Y向間的切向衰減控制因子， ω為角頻率，為頻率的2 π倍， i為虛數(shù)單位。
2. 根據(jù)權(quán)利要求1所述的一種適用于交錯(cuò)網(wǎng)格有限差分的全吸收PML方法，其特征在 于，所述切向衰減控制因子mx/ y、mx/z、my/x、m y/z、mz/x、mz/y的取值范圍為0· 005?0· 020。
【文檔編號(hào)】G01V1/40GK104237944SQ201410526208
【公開日】2014年12月24日 申請(qǐng)日期:2014年10月9日 優(yōu)先權(quán)日:2014年10月9日
【發(fā)明者】王兵, 張闊 申請(qǐng)人:王兵