專利名稱：一種地震波動(dòng)方程生成方法及系統(tǒng)的制作方法
技術(shù)領(lǐng)域：
本發(fā)明有關(guān)于地質(zhì)勘探技術(shù)領(lǐng)域，尤其是有關(guān)于地震勘探技術(shù)領(lǐng)域，具體地講是一種地震波動(dòng)方程生成方法及系統(tǒng)。

背景技術(shù)：
在油氣勘探中，能更好地、更準(zhǔn)確地知道地質(zhì)結(jié)構(gòu)尤為重要，而知道了地質(zhì)結(jié)構(gòu)以后就可以判斷出氣田、油田的位置，以及地震傳播的地震波動(dòng)傳播過程。
而目前，研究地震波動(dòng)傳播過程是通過求出地震波動(dòng)方程來(lái)進(jìn)行，即通過地震波場(chǎng)數(shù)值模擬技術(shù)來(lái)進(jìn)行模擬。
上述的地震波場(chǎng)數(shù)值模擬技術(shù)是在地下介質(zhì)結(jié)構(gòu)和參數(shù)已知的情況下，利用理論計(jì)算的方法研究地震波在地下介質(zhì)中的傳播規(guī)律，并獲得人工合成地震記錄的一種技術(shù)。國(guó)內(nèi)外在地震勘探生產(chǎn)實(shí)踐中所采用的地震模型數(shù)值模擬方法有射線方法(射線追蹤法)和波動(dòng)方程數(shù)值解法(克?；舴蚍e分法、有限差分法、有限單元法、邊界單元法和虛譜法)。
射線追蹤的理論基礎(chǔ)是在高頻近似條件下，地震波場(chǎng)的主能量沿射線軌跡傳播。主要問題在于①難以處理介質(zhì)中較強(qiáng)的速度變化；②難以求出多值走時(shí)中的全局最小走時(shí)；③復(fù)雜波場(chǎng)的振幅信息不準(zhǔn)確；④陰影區(qū)內(nèi)射線覆蓋密度不足。
波動(dòng)方程數(shù)值解法除了獲得走時(shí)信息，還可得到波場(chǎng)的振幅信息，能描述在復(fù)雜介質(zhì)中波的傳播過程。
有限差分法是一種最常用的正演模擬方法，現(xiàn)已比較成熟，正向提高精度的方向發(fā)展。但差分法是從整體上考慮用一組差分方程來(lái)代替微分方程，對(duì)于復(fù)雜的結(jié)構(gòu)形式及邊界條件難以處理，對(duì)于高頻信號(hào)分辨能力受到限制。
偽譜法做波動(dòng)方程傳播的正演已有一段時(shí)間，早在70年代初期就有人介紹該方法，它在氣象預(yù)報(bào)、非線性波動(dòng)和渦流模擬等領(lǐng)域得到廣泛應(yīng)用(最初由Kreiss和Oliger(1972)建議)。大約十年后，由Gazdag(1981)、Kosloff和Baysal(1982)最先應(yīng)用于地震波動(dòng)問題，虛譜法開始進(jìn)入地震模擬方面。Kosloff等用虛譜法對(duì)聲波方程進(jìn)行了正演模擬，并與解析計(jì)算結(jié)果和超聲物理模擬進(jìn)行了比較，證明了方法的正確性。但是該方法處理吸收邊界和自由表面邊界條件比較困難，計(jì)算量和占用內(nèi)存也較大。
有限元法在70年代將有限元法引入地震勘探領(lǐng)域，對(duì)其理論和應(yīng)用做了大量研究和探索。有限元法可以從標(biāo)量、矢量運(yùn)動(dòng)方程出發(fā)，不限于特定的坐標(biāo)系，對(duì)整個(gè)區(qū)域進(jìn)行靈活剖分，適用于均勻、非均勻地層所組成的復(fù)雜構(gòu)造形態(tài)，但是當(dāng)計(jì)算域較大時(shí)，用有限元法求解，因其空間坐標(biāo)離散后未知數(shù)多，勢(shì)必導(dǎo)致計(jì)算機(jī)耗時(shí)量大，計(jì)算成本高。目前對(duì)三維地震正演問題硬件條件尚不夠。
目前三維正演模擬由于計(jì)算效率低，精度也不高，相對(duì)模型也比較簡(jiǎn)單，無(wú)法滿足生產(chǎn)的需要。因此迫切需要有快速的、高精度的三維正演方法。


發(fā)明內(nèi)容
本發(fā)明正是基于上述現(xiàn)有的問題而提出，其目的在于提供一種地震波動(dòng)方程生成方法及系統(tǒng)，以快速、準(zhǔn)確的求出地震波動(dòng)方程。
為了實(shí)現(xiàn)上述目的，本發(fā)明的實(shí)施例提供一種地震波動(dòng)方程生成方法，包括以下步驟獲取聲波方程數(shù)據(jù)；獲取地質(zhì)參數(shù)信息；根據(jù)所述的聲波方程數(shù)據(jù)和所述的地質(zhì)參數(shù)信息進(jìn)行任意差分精細(xì)積分，得出地震傳播方程；求出穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件；由所述的地震傳播方程以及所述的穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件生成地震波動(dòng)方程。
為了實(shí)現(xiàn)上述目的，本發(fā)明的實(shí)施例還提供一種地震波動(dòng)方程生成系統(tǒng)，該地震波動(dòng)方程生成系統(tǒng)包括聲波方程獲取單元，用于獲取聲波方程數(shù)據(jù)；地質(zhì)參數(shù)獲取單元，用于獲取地質(zhì)參數(shù)信息；傳播方程生成單元，用于根據(jù)所述的聲波方程數(shù)據(jù)和所述的地質(zhì)參數(shù)信息進(jìn)行任意差分精細(xì)積分，得出地震傳播方程；條件生成單元，用于求出穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件；波動(dòng)方程生成單元，用于由所述的地震傳播方程以及所述的穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件生成地震波動(dòng)方程。
本發(fā)明的任意差分精細(xì)積分方法有較高的精度和較好的數(shù)值穩(wěn)定性，可以在幾乎不增加計(jì)算量的情況下較大地提高了計(jì)算精度和穩(wěn)定性。完全匹配邊界條件可以用于波場(chǎng)邊界上進(jìn)行振幅衰減。用任意差分精細(xì)積分法導(dǎo)出的完全匹配層邊界條件對(duì)邊界反射有很好的衰減作用。本發(fā)明可以為復(fù)雜區(qū)地震波傳播規(guī)律研究提供了一個(gè)高精度、穩(wěn)定性能好的數(shù)據(jù)體。



此處所說明的附圖用來(lái)提供對(duì)本發(fā)明的進(jìn)一步理解，構(gòu)成本申請(qǐng)的一部分，并不構(gòu)成對(duì)本發(fā)明的限定。在附圖中 圖1A所示的是本發(fā)明的實(shí)施方式的地震波動(dòng)方程生成系統(tǒng)的結(jié)構(gòu)框圖； 圖1B所示的是本發(fā)明的實(shí)施方式的地質(zhì)參數(shù)獲取單元的結(jié)構(gòu)框圖； 圖1C所示的是本發(fā)明的實(shí)施方式的傳播方程生成單元的結(jié)構(gòu)框圖； 圖2所示的是本發(fā)明的實(shí)施方式的地震波動(dòng)方程生成方法的流程圖； 圖3所示的是本發(fā)明的實(shí)施方式的進(jìn)行任意差分精細(xì)積分的結(jié)果示意圖； 圖4所示的是本發(fā)明的實(shí)施方式的進(jìn)行任意差分精細(xì)積分的結(jié)果示意圖； 圖5所示的是本發(fā)明的實(shí)施方式的一維波動(dòng)方程完全匹配邊界條件應(yīng)用效果對(duì)比圖； 圖6所示的是本發(fā)明的實(shí)施方式的完全匹配層吸收邊界示意圖。

具體實(shí)施例方式 為使本發(fā)明的目的、技術(shù)方案和優(yōu)點(diǎn)更加清楚明白，下面結(jié)合實(shí)施方式和附圖，對(duì)本發(fā)明做進(jìn)一步詳細(xì)說明。在此，本發(fā)明的示意性實(shí)施方式及其說明用于解釋本發(fā)明，但并不作為對(duì)本發(fā)明的限定。
本發(fā)明實(shí)施例提供一種地震波動(dòng)方程生成方法及系統(tǒng)。以下結(jié)合附圖對(duì)本發(fā)明進(jìn)行詳細(xì)說明。
圖1A所示的是本發(fā)明的實(shí)施方式的地震波動(dòng)方程生成系統(tǒng)的結(jié)構(gòu)框圖、圖1B所示的是本發(fā)明的實(shí)施方式的地質(zhì)參數(shù)獲取單元的結(jié)構(gòu)框圖、圖1C所示的是本發(fā)明的實(shí)施方式的傳播方程生成單元的結(jié)構(gòu)框圖。
其中，本發(fā)明的實(shí)施方式的地震波動(dòng)方程生成系統(tǒng)包括聲波方程獲取單元101，用于獲取聲波方程數(shù)據(jù)；地質(zhì)參數(shù)獲取單元102，用于獲取地質(zhì)參數(shù)信息；傳播方程生成單元103，用于根據(jù)所述的聲波方程數(shù)據(jù)和所述的地質(zhì)參數(shù)信息進(jìn)行任意差分精細(xì)積分，得出地震傳播方程；條件生成單元104，用于求出穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件；波動(dòng)方程生成單元105，用于由所述的地震傳播方程以及所述的穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件生成地震波動(dòng)方程。
上述地質(zhì)參數(shù)獲取單元包括地震勘探模塊106，采集地震數(shù)據(jù)；數(shù)據(jù)處理模塊107，利用promax軟件處理所述地震數(shù)據(jù)得出地震剖面；分析模塊108，對(duì)所述地震剖面進(jìn)行分析得出對(duì)應(yīng)的地質(zhì)參數(shù)信息。
而上述傳播方程生成單元包括有限差分模塊109，用于進(jìn)行時(shí)間有限差分；精細(xì)積分模塊110，用于進(jìn)行空間精細(xì)積分，所述傳播方程生成單元進(jìn)行時(shí)間有限差分和空間精細(xì)積分的交替來(lái)進(jìn)行所述的任意差分精細(xì)積分。
另外，上述邊界條件是指完全匹配層吸收邊界條件，以在完全匹配邊界條件下進(jìn)行任意差分精細(xì)積分，并且求出的地震波動(dòng)方程為一半解析解。
圖2所示的是本發(fā)明的實(shí)施方式的地震波動(dòng)方程生成方法的流程圖，一下結(jié)合圖2詳細(xì)說明本發(fā)明的地震波動(dòng)方程生成方法。
S201獲取聲波方程數(shù)據(jù)； S202采集地震數(shù)據(jù)； S203利用promax軟件處理所述地震數(shù)據(jù)得出地震剖面； S204對(duì)所述地震剖面進(jìn)行分析得出對(duì)應(yīng)的地質(zhì)參數(shù)信息； S205進(jìn)行時(shí)間有限差分； S206交替進(jìn)行空間精細(xì)積分得出地震傳播方程； S207求出穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件； S208由所述的地震傳播方程以及所述的穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件生成地震波動(dòng)方程。
上述邊界條件是指完全匹配層吸收邊界條件，以在完全匹配邊界條件下進(jìn)行任意差分精細(xì)積分，且上述地震波動(dòng)方程為一半解析解。
以下詳細(xì)說明本發(fā)明解決其技術(shù)問題所采用的技術(shù)手段 一維波動(dòng)方程任意差分精細(xì)積分原理 本節(jié)以一維波動(dòng)方程的初邊值問題為例來(lái)闡述任意差分精細(xì)積分法的基本原理。
設(shè)有如下初邊值問題 式中，u代表位移，t是時(shí)間變量，x是空間坐標(biāo)，v是波的傳播速度，L是一維區(qū)域的長(zhǎng)度。
對(duì)(1-1)中的第一式進(jìn)行空間坐標(biāo)離散，在任一點(diǎn)j(點(diǎn)序號(hào))處，有 這里，J是總離散點(diǎn)數(shù)。將函數(shù)u(x，t)在j點(diǎn)的鄰域內(nèi)作泰勒展開，得到 其中，Δxi＝xi-xj，i是j點(diǎn)鄰域內(nèi)的任意點(diǎn)，對(duì)(1-3)式移項(xiàng)并取n個(gè)加權(quán)參數(shù)后再求和得 其中αi是加權(quán)系數(shù)，取m＝n并省略掉(1-4)式右端項(xiàng)的高階項(xiàng)O[(Δxi)m+1]，再令其等于(1-2)式的右端，得到 上式是一個(gè)關(guān)于αi的m維線性方程組，解(1-5)可求出加權(quán)系數(shù)αi。
(對(duì)于均勻網(wǎng)格，m＝2的情形，有)。
(1-2)式可近似為 在tn時(shí)刻附近的小領(lǐng)域內(nèi)，(1-6)可表示為下式 (上式u的第一個(gè)下標(biāo)代表空間變量，第二個(gè)下標(biāo)代表時(shí)間變量，以下同) (1-7)式右端的uj移到方程的左端得tn時(shí)刻小領(lǐng)域內(nèi)常微分方程 在tn時(shí)刻小領(lǐng)域內(nèi)，將上式右端看成是不隨時(shí)間變化的常數(shù)項(xiàng)，其解析解為 式中， 可以看出(1-9)式有兩個(gè)未知數(shù)c1和c2，選擇在子域[tn-1，tn+1]積分，假設(shè)tn-1，tn時(shí)刻的值已知，求出tn+1時(shí)刻的遞推公式 為了進(jìn)一步提高差分格式的精度，對(duì)(1-6)移項(xiàng) 將(1-11)式的右端項(xiàng)在tn的鄰域內(nèi)做泰勒展開 根據(jù)(1-11)式的右端項(xiàng)不同的截?cái)?，可以得到不同的?jì)算格式，展開到t-tn的二階導(dǎo)數(shù)，得 解常微分方程 得解析解(此即精細(xì)積分的意義)為 上式中，對(duì)(1-14)式選擇積分子域?yàn)閇tn-1，tn+1]，利用tn-1，tn時(shí)刻的值，確定兩個(gè)未知數(shù)c1和c2，tn+1時(shí)刻的遞推計(jì)算公式為 上式中a，bn同(1-9)式，而dn待求，下面求dn 由于

的計(jì)算可由(1-6)式求得即求得(1-16) 同理可證，當(dāng)(1-11)式的右端項(xiàng)進(jìn)行t-tn的零階展開和t-tn的一階展開時(shí)，得到的遞推公式相同，都為 由于(1-15)是時(shí)間二階展開，其截?cái)嗾`差較小，因此理論上，(1-15)式較(1-10)具有更高的時(shí)間精度。
二維及三維波動(dòng)方程的精細(xì)積分方法的導(dǎo)出 精細(xì)積分法是采用差分的方法計(jì)算空間偏導(dǎo)數(shù)，采用子域精細(xì)積分的方法計(jì)算時(shí)間偏導(dǎo)數(shù)。對(duì)于二維波動(dòng)方程 首先對(duì)上式中的第一式進(jìn)行二維空間坐標(biāo)離散(按均勻網(wǎng)格劃分)，在任一離散點(diǎn)j(點(diǎn)序號(hào))處，有 上式中，(0≤j≤J)，J是總離散點(diǎn)數(shù)。
對(duì)(2-2)式u在j點(diǎn)的x鄰域作泰勒展開 上式中，Δxi＝xi-xj，i是j的x鄰域內(nèi)的任意點(diǎn)，對(duì)上式移項(xiàng)并加權(quán)得 其中αi是加權(quán)系數(shù)，m原則上可任意選取，但m>n使下式的方程組成為矛盾方程組，m<n使下式的方程組為欠定的，在此一般取m＝n(以m＝n為例來(lái)展開討論)， 省略掉(2-3)式右端項(xiàng)的高階項(xiàng)O[(Δxi)m+1]，再令其等于(2-2)式的右端項(xiàng)

則有 上式是一個(gè)關(guān)于αi的m維線性方程組，解之可求出加權(quán)系數(shù)αi。
令 同理，

也可以這樣處理，將u在j點(diǎn)的z鄰域作泰勒展開 上式中，Δzp＝zp-zj，p是j的z鄰域內(nèi)的任意點(diǎn)，對(duì)上式移項(xiàng)并加權(quán)得 其中αp是加權(quán)系數(shù)，m原則上可任意選取，但m>n使下式的方程組成為矛盾方程組，m<n使下式的方程組為欠定的，在此一般取m＝n(以m＝n為例來(lái)展開討論)，省略掉(2-4)式右端項(xiàng)的高階項(xiàng)O[(Δzp)m+1])，再令其等于(2-2)式的右端項(xiàng)

則有 上式是一個(gè)關(guān)于αp的m維線性方程組，解之可求出加權(quán)系數(shù)αp。對(duì)αp，αi重新排列，并統(tǒng)一用αi表示。(2-2)右端的二階導(dǎo)數(shù)可用相鄰點(diǎn)的差分來(lái)表示。暫不考慮震源項(xiàng)δ(x-x0，z-z0)f(t)，那么(2-2)式可表示為 該式的右端可近似為tn時(shí)刻，關(guān)于uj微分方程，因此上式可寫為 在tn時(shí)刻的小鄰域內(nèi)(2-5)的右端項(xiàng)可展開到t-tn的二階項(xiàng)來(lái)近似，即 將上式代入方程(2-5)，求得的解析解為 式中， 選擇積分子域[tn-1，tn+1]，假定tn-1，tn的波場(chǎng)已知，可確定c1和c2，進(jìn)而可以求出tn+1時(shí)刻的值 對(duì)(2-5)式的右端項(xiàng)進(jìn)行t-tn的零階展開和t-tn的一階展開時(shí)，得到的遞推公式相同，都為 從(2-15)和(2-6)式可以看出，一維與二維波動(dòng)方程的ADPI(任意差分精細(xì)積分法)方法計(jì)算公式形式上完全相同，但實(shí)際上公式中的bn和dn是與空間維數(shù)密切相關(guān)的。
對(duì)于三維波動(dòng)方程 具體做法為 首先對(duì)上式中的第一式進(jìn)行三維空間坐標(biāo)離散(即按均勻網(wǎng)格劃分)，在任一離散點(diǎn)j(點(diǎn)序號(hào))處，有 上式中，(0≤j≤J)，J是總離散點(diǎn)數(shù) (2-9)式的右端項(xiàng)

可用任意差分的方法求得，將u在j點(diǎn)的x鄰域作泰勒展開 上式中，Δxi＝xi-xj，i是j的x鄰域內(nèi)的任意點(diǎn)，對(duì)上式移項(xiàng)并加權(quán)得 其中αi是加權(quán)系數(shù)，一般取m＝n，省略掉(2-3)式右端項(xiàng)的高階項(xiàng)O[(Δxi)m+1]，再令其等于(2-2)式的右端項(xiàng)

則有 上式是一個(gè)關(guān)于αi的m維線性方程組，解之可求出加權(quán)系數(shù)αi。
令 同理，

也可以這樣處理，將u在j點(diǎn)的z鄰域作Taylor展開 上式中，Δzp＝zp-zj，p是j的z鄰域內(nèi)的任意點(diǎn)，對(duì)上式移項(xiàng)并加權(quán)得 取m＝n，省略掉(2-11)式右端項(xiàng)的高階項(xiàng)O[(Δzp)m+1])，再令其等于(2-9)式的右端項(xiàng)

則有 上式是一個(gè)關(guān)于αp的m維線性方程組，解之可求出加權(quán)系數(shù)αp。
令 對(duì)

做與

相同的處理，將u在j點(diǎn)的y鄰域作泰勒展開 上式中，Δyq＝y(tǒng)q-zj，q是j的y鄰域內(nèi)的任意點(diǎn)，對(duì)上式移項(xiàng)并加權(quán)得 取m＝n，省略掉(2-12)式右端項(xiàng)的高階項(xiàng)O[(Δyq)m+1]，再令其等于(2-9)式的右端項(xiàng)

則有 上式是一個(gè)關(guān)于αq的m維線性方程組，解之可求出加權(quán)系數(shù)αq。
令 對(duì)αp，αq，αi重新排列，并統(tǒng)一用αi表示。則αi是大小為3*m的系數(shù)。
(2-2)右端的二階導(dǎo)數(shù)可用相鄰點(diǎn)的任意差分來(lái)表示，暫不考慮震源，那么可表示為 該式的右端可近似為tn時(shí)刻，關(guān)于uj微分方程，因此上式可寫為 將(2-13)的右端項(xiàng)展開到t-tn的二階，即 這時(shí)的解析解為 上式中， 選擇積分子域[tn-1，tn+1]，假定tn-1，tn的波場(chǎng)已知，可確定c1和c2，進(jìn)而可以求出tn+1時(shí)刻的值 上式中， 當(dāng)對(duì)(2-6)式的右端項(xiàng)進(jìn)行t-tn的零階展開和t-tn的一階展開時(shí)，得到的遞推公式相同，都為 運(yùn)用公式(2-14)進(jìn)行數(shù)值求解時(shí)以下幾個(gè)問題必須注意 ①對(duì)

項(xiàng)，公式(2-14)在鄰域[tn-1，tn+1]內(nèi)二階展開，隨著展開的階數(shù)增高，數(shù)值求解的精度也會(huì)相應(yīng)提高。當(dāng)然計(jì)算公式的復(fù)雜性和計(jì)算花費(fèi)的時(shí)間也會(huì)增加。
②dn的計(jì)算涉及到求

用預(yù)報(bào)--校正方法可提高計(jì)算精度，首先用(2-15)式預(yù)測(cè)uj，n+1，用下列中心差分計(jì)算公式預(yù)測(cè)

和
(2-7)的第一個(gè)公式簡(jiǎn)化為校正公式 預(yù)報(bào)一校正方法比較實(shí)用，精度較高，且節(jié)約了計(jì)算量，本文使用了這種方法。
③要使得ADPI方法成為實(shí)用的地震正演技術(shù)，一方面，需要用理論解析解作為參照，證明方法的正確性；另一方面，必須考慮如何確定將實(shí)際半無(wú)窮空間波場(chǎng)轉(zhuǎn)化為有限區(qū)域波場(chǎng)的邊界條件。同時(shí)要研究方法的數(shù)值計(jì)算穩(wěn)定性，并與實(shí)際科研生產(chǎn)中正在使用的其它地震正演技術(shù)對(duì)比，說明ADPI方法的精度。
圍繞問題③展開ADPI方法的進(jìn)一步研究，下面的內(nèi)容將逐一回答上述的計(jì)算公式的正確性、邊界條件、穩(wěn)定性、與數(shù)值計(jì)算的精確度問題。
ADPI方法穩(wěn)定性分析 采用Von Neumann方法來(lái)分析(1-15)式的數(shù)值穩(wěn)定性。首先介紹下面的引理。
引理實(shí)系數(shù)二次方程μ2-bμ-c＝0的根，按其模小于或等于1的充分且必要條件是|b|≤1-c，|c|≤1。
為了便于與有限差分法比較，這里考慮均勻網(wǎng)格、當(dāng)m＝2時(shí)，可以得到 令uj，n＝VneikjΔx，代入(1-15)式，得到下式 根據(jù)歐拉公式eikΔx＝cos(kΔx)+i sin(kΔx)，消去公因子eikjΔx得
其中λ＝cos(aΔt)，β＝cos(kΔx)。由于(4-1)式是三層格式，Vn+1與Vn、Vn-1有關(guān)，故令 Vn＝1×Vn+0×Vn-1(4-2) 將(4-1)式和(4-2)式合并寫成 式(4-3)中過渡矩陣為 G的特征值μ滿足 μ2-bμ-c＝0 其中b＝2λ+2(2β-β2)(1-λ)+2αβ(β-1)(Δt)2，c＝-1 當(dāng)時(shí)，(aΔt)2≈2(1-cos(aΔt)) b≈2λ+2β(1-λ) 所以|b|≤2＝1-c，|c|＝1 由引理知，當(dāng)時(shí)，G的特征值小于等于1。
定理1由Von Neumann條件知(1-15)式穩(wěn)定。即當(dāng)時(shí)，四點(diǎn)二階精細(xì)積分波動(dòng)方程解法穩(wěn)定。
對(duì)零階展開和一階展開得到的計(jì)算格式(1-17)，穩(wěn)定性條件仍為 對(duì)于(1-1)式的有限差分格式，由參考文獻(xiàn)可知其穩(wěn)定性條件為對(duì)比可知精細(xì)積分法比有限差分法更為穩(wěn)定。
下面用數(shù)值求解的例子直觀地說明了ADPI方法的穩(wěn)定性比有限差分格式更好。圖3表明，當(dāng)dt＝0.00025s，時(shí)，有限差分法不穩(wěn)定，結(jié)果發(fā)散。而ADPI方法穩(wěn)定，結(jié)果正確；圖4表明，當(dāng)dt＝0.00075s，有限差分法不穩(wěn)定，結(jié)果發(fā)散。而ADPI方法穩(wěn)定，結(jié)果正確。
例1dt＝0.00025s， 其中，圖3的左圖有限差分法的結(jié)果，不穩(wěn)定，解發(fā)散； 圖3的右圖二階ADPI法的結(jié)果，穩(wěn)定，解正確。
例2dt＝0.00075s， 其中，圖4的左圖有限差分法的結(jié)果，不穩(wěn)定，解發(fā)散 圖4的右圖二階ADPI法的結(jié)果，穩(wěn)定，解正確 邊界條件分析 利用波動(dòng)方程模擬地震記錄時(shí)，吸收邊界條件是非常重要的。在一個(gè)有限的區(qū)域中，利用數(shù)值方法進(jìn)行計(jì)算時(shí)，如果不對(duì)邊界進(jìn)行處理，就會(huì)產(chǎn)生邊界反射?，F(xiàn)有許多方法用以消除邊界反射。Cerjan等人通過在計(jì)算邊界附近引入損耗介質(zhì)來(lái)衰減向外傳播的波。由于在兩個(gè)具有不同吸收系數(shù)的損耗介質(zhì)的界面上會(huì)產(chǎn)生反射，這就要求有較多的吸收層，才能有較好的吸收效果，而較多的吸收層意味著較大的計(jì)算量。另一種比較著名的吸收邊界條件是Clayton吸收邊界條件，該吸收邊界條件在旁軸近似理論的基礎(chǔ)上導(dǎo)出，在特定的入射角和頻率范圍內(nèi)，具有較好的吸收效果。Berenger針對(duì)電磁波傳播情況，給出了一種高效的完全匹配層吸收邊界條件，并在理論上證明該方法可以完全吸收來(lái)自各個(gè)方向、各種頻率的電磁波，而不產(chǎn)生任何反射。這里將完全匹配層的思想引入到三維波動(dòng)方程ADPI解法中，并推導(dǎo)了三維波動(dòng)方程ADPI解法的完全匹配層吸收邊界條件。
①一維波動(dòng)方程完全匹配層原理 在一維聲波介質(zhì)中，波的傳播可以用如下方程來(lái)描述 式中u(x，t)表示位移；V(x)表示速度。通過引入中間變量A(x，t)此方程也可寫為兩個(gè)一階偏微分方程組的形式 容易證明，式(5-1)和式(5-2)是等價(jià)的。由式(5-2)構(gòu)建新的方程組得 式中p*(x，t)表示新方程的解；d(x)為一個(gè)不隨時(shí)間變化的邊界衰減函數(shù)。
下面研究式(5-2)的解p(x，t)和式(5-3)的解p*(x，t)的關(guān)系。首先對(duì)式(5-3)做關(guān)于時(shí)間的Fourier變換得 其中，

分別代表p*(x，ω)、A*(x，ω)的Fourier變換。
對(duì)x做坐標(biāo)變換 則有 將式(5-7)、式(5-8)分別代入式(5-4)并整理得 對(duì)式(5-9)做關(guān)于時(shí)間的Fourier逆變換，可得 現(xiàn)在來(lái)比較式(5-2)和式(5-10)，兩對(duì)式子(方程組)在形式上完全相同，所以他們有相同形式的解，但是在不同的空間坐標(biāo)中，假設(shè)式(5-2)的解是p(x，t)，則式(5-10)的解就是

即式(5-3)的解為 也就是說，可以先求解式(5-2)，然后再對(duì)它們的解做坐標(biāo)變換，就可得到式(5-3)的解?，F(xiàn)在考察式(5-3)解的性質(zhì)。式(5-2)在均勻介質(zhì)中有如下形式的特解 p(x，t)＝p0exp[-i(kxx-ωt)] (5-11) 由上述討論可知，對(duì)應(yīng)的式(5-3)的解為 將式(5-5)代入上式，并簡(jiǎn)化可得 上述解的振幅比為 此式說明，式(5-3)的解相對(duì)式(5-2)的解是衰減的，而d(x)起到了一個(gè)衰減系數(shù)的作用，波場(chǎng)按傳播距離的呈指數(shù)衰減，衰減速度很快。因而不會(huì)產(chǎn)生邊界反射。
下面用一個(gè)簡(jiǎn)單的示例說明一維波動(dòng)方程完全匹配的效果。
圖5是一維波動(dòng)方程完全匹配邊界條件應(yīng)用效果對(duì)比圖。
其中，(a)未加衰減量的一維波動(dòng)方程正演；(b)加入衰減量后的一維波動(dòng)方程正演；(c)進(jìn)一步加大衰落減量的一維波動(dòng)方程正演。
從圖5可以看出加入吸收系數(shù)后振幅開始衰減，衰減系數(shù)越大，振幅衰減越快。示例表明完全匹配層邊界條件是可以用于波場(chǎng)邊界上進(jìn)行振幅衰減的，而衰減函數(shù)d(x)選擇可以靈活調(diào)整完全匹配層邊界條件。
②二維波動(dòng)方程完全匹配層原理 對(duì)下列的二維波動(dòng)方程 其中p(x，z，t)是位移函數(shù)；V(x，z)是介質(zhì)的速度，此方程可寫為等價(jià)的一階偏微分方程組的形式 其中u1、u2、A1、A2是引入的中間變量。
相應(yīng)的完全匹配層控制方程為 和一維情況相同，式(5-17)的解也是衰減的。d1(x)和d2(z)分別為x方向和z方向的衰減系數(shù)，也就是說d1(x)起到衰減x方向傳播的波、d2(z)起到衰減z方向傳播的波的作用。對(duì)于任意方向傳播的波，可以通過矢量分解，分解成x方向和z方向傳播的波，分別進(jìn)行衰減。波場(chǎng)是按傳播距離的指數(shù)規(guī)律衰減，衰減速度很快。當(dāng)衰減系數(shù)d1(x)、d2(z)隨空間位置變化時(shí)，不會(huì)在介質(zhì)中產(chǎn)生任何反射。這些性質(zhì)使得這種介質(zhì)特別適合用作波動(dòng)方程的邊界吸收介質(zhì)。
圖6是完全匹配層吸收邊界示意圖。利用完全匹配層作為吸收邊界的基本做法是在所研究區(qū)域的四周引入完全匹配層。如圖6所示，區(qū)域ABCD為所要研究的區(qū)域，即我們要在此區(qū)域中研究波的傳播問題。在區(qū)域的周圍加上完全匹配層，在區(qū)域1中，令d1(x)≠0；d2(z)≠0，速度V都等于角點(diǎn)的速度。在區(qū)域2中，令d1(x)＝0；d2(z)≠0，速度V在z方向?yàn)槌?shù)，在x方向和邊界的速度相等。在區(qū)域3中，令d1(x)≠0；d2(z)＝0，速度V在x方向?yàn)槌?shù)，在z方向和邊界的速度相等。這樣在計(jì)算邊界的周圍都有完全匹配層吸收介質(zhì)，波由區(qū)域內(nèi)通過邊界傳播到完全匹配層時(shí)，不會(huì)產(chǎn)生任何反射。波在完全匹配層中傳播時(shí)，不會(huì)產(chǎn)生反射，并且按傳播距離的指數(shù)規(guī)律衰減。當(dāng)波傳播到完全匹配層的邊界時(shí)，波場(chǎng)近似為零，也不會(huì)產(chǎn)生反射。
同理，可以推出到三維波動(dòng)方程完全匹配層的邊界，其形式與二維類同。
③二維波動(dòng)方程任意差分精細(xì)積分法求解完全匹配層計(jì)算公式 用任意差分精細(xì)積分進(jìn)行完全匹配邊界計(jì)算，需要對(duì)(5-17)式方程組在計(jì)算區(qū)域內(nèi)和邊界用同樣方法求解，若僅在邊界上的用完全匹配層解(5-17)式，就將出現(xiàn)波場(chǎng)的不相容性，從而出現(xiàn)數(shù)值計(jì)算引起的波場(chǎng)反射。為了消除這種反射，必須推導(dǎo)方程組(5-17)式的任意差分精細(xì)積分完全匹配邊界計(jì)算公式。
由(5-17)式，令v＝u*1，w＝u*2，px＝A*1，pz＝A*2整理后得到 上式的右端項(xiàng)可在點(diǎn)的x鄰域作泰勒展開 移項(xiàng)加權(quán)得 令(5-20)右端展開并等于(5-18)關(guān)于x的偏導(dǎo)數(shù)得 解該方程組即可求出加權(quán)系數(shù)αi即 同理在j點(diǎn)的z點(diǎn)鄰域作泰勒展開 移項(xiàng)加權(quán)得 令上式右端展開并等于關(guān)于z的偏導(dǎo)數(shù)得 解該方程組即可求出加權(quán)系數(shù)αl即 整理得 因?yàn)樗杂? 設(shè) d(s)＝0.5×pml_max×[1-cos(2×π×s)/pml_thick)] (5-27) 參數(shù)s的選擇根據(jù)波的傳播方向而定，若是向x方向傳播，則s為x。若為z方向則s選為z。pml_max為衰減系數(shù)，其大小的選定可根據(jù)邊界大小而定，但不可過大，會(huì)造成計(jì)算的不穩(wěn)定，pml_thick為完全匹配層的寬度。
即 解此方程得 整理得 選擇積分子域[tn，tn+1]求出tn+1時(shí)刻的值
進(jìn)一步提高精度可將展開到t-tn的二階導(dǎo)數(shù)項(xiàng) 同理將展開為t-tn的二階導(dǎo)數(shù)項(xiàng)得 代入并整理(5-18)得 令 即 解此方程得 代換C得任意差分精細(xì)積分進(jìn)行完全匹配邊界遞推求解計(jì)算公式 ④三維波動(dòng)方程任意差分精細(xì)積分法求解完全匹配層計(jì)算公式 根據(jù)以上二維方程推導(dǎo)同理可得三維方程完全匹配邊界條件的求解公式。
三維波動(dòng)方程等價(jià)的一階偏微分方程組可由下式表示 令 作泰勒展開并整理得 令 解得 選擇積分子域[tn，tn+1]確定C得
進(jìn)一步提高精度后得 令 即 公式(5-44)是三維波動(dòng)方程完全匹配邊界條件的求解遞推公式，在波場(chǎng)內(nèi)部其解滿足方程(2-1)。波場(chǎng)在邊界傳送過程中由于d(s)的衰減作用使得波場(chǎng)振幅很快衰減從而達(dá)到去除邊界反射波的目的。
本發(fā)明的任意差分精細(xì)積分方法有較高的精度和較好的數(shù)值穩(wěn)定性，可以在幾乎不增加計(jì)算量的情況下較大地提高了計(jì)算精度和穩(wěn)定性。完全匹配邊界條件可以用于波場(chǎng)邊界上進(jìn)行振幅衰減。用任意差分精細(xì)積分法導(dǎo)出的完全匹配層邊界條件對(duì)邊界反射有很好的衰減作用。本發(fā)明可以為復(fù)雜區(qū)地震波傳播規(guī)律研究提供了一個(gè)高精度、穩(wěn)定性能好的數(shù)據(jù)體。
以上所述的具體實(shí)施方式
，對(duì)本發(fā)明的目的、技術(shù)方案和有益效果進(jìn)行了進(jìn)一步詳細(xì)說明，所應(yīng)理解的是，以上所述僅為本發(fā)明的具體實(shí)施方式
而已，并不用于限定本發(fā)明的保護(hù)范圍，凡在本發(fā)明的精神和原則之內(nèi)，所做的任何修改、等同替換、改進(jìn)等，均應(yīng)包含在本發(fā)明的保護(hù)范圍之內(nèi)。
權(quán)利要求
1.一種地震波動(dòng)方程生成方法，其特征在于，所述方法包括以下步驟
獲取聲波方程數(shù)據(jù)；
獲取地質(zhì)參數(shù)信息；
根據(jù)所述的聲波方程數(shù)據(jù)和所述的地質(zhì)參數(shù)信息進(jìn)行任意差分精細(xì)積分，得出地震傳播方程；
求出穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件；
由所述的地震傳播方程以及所述的穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件生成地震波動(dòng)方程。
2.根據(jù)權(quán)利要求1所述的地震波動(dòng)方程生成方法，其特征在于，所述獲取地質(zhì)參數(shù)信息的步驟包括
采集地震數(shù)據(jù)；
利用promax軟件處理所述地震數(shù)據(jù)得出地震剖面；
對(duì)所述地震剖面進(jìn)行分析得出對(duì)應(yīng)的地質(zhì)參數(shù)信息。
3.根據(jù)權(quán)利要求1所述的地震波動(dòng)方程生成方法，其特征在于，所述進(jìn)行任意差分精細(xì)積分的步驟包括
進(jìn)行時(shí)間有限差分和空間精細(xì)積分的交替來(lái)進(jìn)行所述的任意差分精細(xì)積分。
4.根據(jù)權(quán)利要求1所述的地震波動(dòng)方程生成方法，其特征在于，所述邊界條件是指完全匹配層吸收邊界條件，以在完全匹配邊界條件下進(jìn)行任意差分精細(xì)積分。
5.根據(jù)權(quán)利要求1所述的地震波動(dòng)方程生成方法，其特征在于，所述地震波動(dòng)方程為一半解析解。
6.一種地震波動(dòng)方程生成系統(tǒng)，其特征在于，該地震波動(dòng)方程生成系統(tǒng)包括
聲波方程獲取單元，用于獲取聲波方程數(shù)據(jù)；
地質(zhì)參數(shù)獲取單元，用于獲取地質(zhì)參數(shù)信息；
傳播方程生成單元，用于根據(jù)所述的聲波方程數(shù)據(jù)和所述的地質(zhì)參數(shù)信息進(jìn)行任意差分精細(xì)積分，得出地震傳播方程；
條件生成單元，用于求出穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件；
波動(dòng)方程生成單元，用于由所述的地震傳播方程以及所述的穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件生成地震波動(dòng)方程。
7.根據(jù)權(quán)利要求6所述的地震波動(dòng)方程生成系統(tǒng)，其特征在于，所述地質(zhì)參數(shù)獲取單元包括
地震勘探模塊，采集地震數(shù)據(jù)；
數(shù)據(jù)處理模塊，利用promax軟件處理所述地震數(shù)據(jù)得出地震剖面；
分析模塊，對(duì)所述地震剖面進(jìn)行分析得出對(duì)應(yīng)的地質(zhì)參數(shù)信息。
8.根據(jù)權(quán)利要求6所述的地震波動(dòng)方程生成系統(tǒng)，其特征在于，所述傳播方程生成單元包括
有限差分模塊，用于進(jìn)行時(shí)間有限差分；
精細(xì)積分模塊，用于進(jìn)行空間精細(xì)積分，
所述傳播方程生成單元進(jìn)行時(shí)間有限差分和空間精細(xì)積分的交替來(lái)進(jìn)行所述的任意差分精細(xì)積分。
9.根據(jù)權(quán)利要求6所述的地震波動(dòng)方程生成系統(tǒng)，其特征在于，所述邊界條件是指完全匹配層吸收邊界條件，以在完全匹配邊界條件下進(jìn)行任意差分精細(xì)積分。
10.根據(jù)權(quán)利要求6所述的地震波動(dòng)方程生成系統(tǒng)，其特征在于，所述地震波動(dòng)方程為一半解析解。
全文摘要
本發(fā)明提供一種地震波動(dòng)方程生成方法及系統(tǒng)，該地震波動(dòng)方程生成方法包括獲取聲波方程數(shù)據(jù)；獲取地質(zhì)參數(shù)信息；根據(jù)所述的聲波方程數(shù)據(jù)和所述的地質(zhì)參數(shù)信息進(jìn)行任意差分精細(xì)積分，得出地震傳播方程；求出穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件；由所述的地震傳播方程以及所述的穩(wěn)定性條件、初始輸入條件、邊界處理?xiàng)l件生成地震波動(dòng)方程。本發(fā)明有較高的精度和較好的數(shù)值穩(wěn)定性，可以在幾乎不增加計(jì)算量的情況下較大地提高了計(jì)算精度和穩(wěn)定性。用任意差分精細(xì)積分法導(dǎo)出的完全匹配層邊界條件對(duì)邊界反射有很好的衰減作用。本發(fā)明可以為復(fù)雜區(qū)地震波傳播規(guī)律研究提供了一個(gè)高精度、穩(wěn)定性能好的數(shù)據(jù)體。
文檔編號(hào)G01V1/28GK101369024SQ20081011976
公開日2009年2月18日 申請(qǐng)日期2008年9月9日 優(yōu)先權(quán)日2008年9月9日
發(fā)明者王潤(rùn)秋, 胡天躍 申請(qǐng)人:中國(guó)石油天然氣集團(tuán)公司, 中國(guó)石油大學(xué)(北京), 北京大學(xué)