本發(fā)明屬于地球物理勘探中彈性波波場數(shù)值模擬技術(shù)領(lǐng)域，具體地涉及基于二維各向同性介質(zhì)的彈性波廣角吸收邊界條件的確定方法。

背景技術(shù)：

彈性波波場數(shù)值模擬技術(shù)是地震勘探中的重要技術(shù)，利用彈性波波場數(shù)值模擬技術(shù)中的正演技術(shù)，可以檢驗(yàn)經(jīng)過地震資料處理、地質(zhì)解釋最終得到的地下物性參數(shù)以及反演的速度模型是否能真實(shí)地反映地下的儲(chǔ)油構(gòu)造，如果離開了彈性波波場數(shù)值模擬技術(shù)，我們目前的許多地震資料的處理和解釋工作將寸步難行。目前廣泛應(yīng)用于實(shí)際工作中的彈性波波場數(shù)值模擬技術(shù)主要是彈性波波場有限差分波場數(shù)值模擬技術(shù)，在彈性波有限差分波場數(shù)值模擬中，由于計(jì)算模型的限制，人們所碰到的一個(gè)不可避免的問題是如何解決由計(jì)算網(wǎng)格邊界所引起的人為邊界反射能量，這些來自邊界的人為反射扭曲了波在無限介質(zhì)中傳播的真實(shí)性，因此有必要提出一種波在邊界上傳播的算法能使往外傳播的波的能量透過邊界，否則，我們將不得不擴(kuò)大模型的邊界，大大增加計(jì)算的工作量。

為了解決邊界反射問題，幾十年來，人們進(jìn)行了大量的研究，提出了許多不同的有關(guān)邊界波場的計(jì)算方法，但至今仍沒有一種十分完美的解決方法。而與聲波波場數(shù)值模擬相比，彈性波有限差分波場數(shù)值模擬中的邊界問題更復(fù)雜、也更難解決。目前，彈性波有限差分波場數(shù)值模擬中的邊界問題研究仍是該領(lǐng)域里的一個(gè)重要課題，同時(shí)也相繼誕生了一些針對(duì)不同邊界問題的算法，其中通過波場吸收消除來自邊界反射的方法比較通用且效果也較明顯，國內(nèi)國外都有不少這方面的研究成果出現(xiàn)和文獻(xiàn)發(fā)表，目前勘探中廣泛使用兩類吸收邊界條件：海綿吸收邊界條件和傍軸近似吸收邊界條件。海綿吸收邊界條件利用黏滯邊界或靠近邊界的條帶范圍內(nèi)對(duì)入射波進(jìn)行衰減；傍軸近似吸收邊界條件基于單程波傳播的方法原理，將不同近似條件下的單程波方程作為不同邊界處的吸收邊界條件。在國內(nèi)，董良國在文獻(xiàn)中(董良國.彈性波數(shù)值模擬中的吸收邊界條件.石油地球物理勘探，1999，34(1):45～56)利用特征分析方法將前人提出的一維情況下無邊界反射概念推廣至三維各向異性介質(zhì)彈性波的數(shù)值模擬中，得到ti介質(zhì)中的吸收邊界條件。在國外，clayton和engquist等人在文獻(xiàn)中(robertclayton，bjornengquist.absorbingboundaryconditionsforacousticandelasticwaveequations.bull.seis.soc.am.，1977，67(6):1529～1540)利用彈性波傳播方程傍軸近似理論得到了各向同性介質(zhì)中彈性波數(shù)值模擬的吸收邊界條件。至今，他們所提出的吸收邊界條件仍是在各向同性介質(zhì)中進(jìn)行彈性波有限差分?jǐn)?shù)值模擬用來解決邊界問題最常用的方法。

此后，clayton和engquist在文獻(xiàn)中(robertclayton，bjornengquist，absorbingboundaryconditionsforwaveequationsmigration.geophysics，1980，45(3):895～904)又對(duì)彈性波有限差分波場數(shù)值模擬中的吸收邊界條件進(jìn)行了更深入的研究。在國內(nèi)，也有不少研究人員對(duì)利用有限差分法進(jìn)行彈性波場數(shù)值模擬情況下的吸收邊界條件做過較詳細(xì)的研究，得到了許多有價(jià)值的研究成果，他們的算法對(duì)某些模型得到了一些較好的結(jié)果，本發(fā)明第一發(fā)明人在文獻(xiàn)中(李文杰等.一種新的彈性波數(shù)值模擬吸收邊界條件.石油地球物理勘探，2009，44(4):501～507)利用二維各向同性介質(zhì)的彈性波傳播方程傍軸近似理論得到了另一種的吸收邊界條件，該方法與clayton和engquist等人提出的吸收邊界方法不一樣，經(jīng)過數(shù)值模擬結(jié)果的證明，該吸收邊界條件具有較好的吸收效果。

但遍觀已有的彈性波波場數(shù)值模擬算法，雖然較好地解決了模型四邊的反射能量問題，卻沒有解決好模型角點(diǎn)的反射能量問題，從模擬效果來看，這些邊界條件不能有效消除來自任意角度角點(diǎn)的反射能量，所以，基于各向同性介質(zhì)的彈性波傍軸近似吸收邊界條件有待進(jìn)一步完善。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的是為了改進(jìn)現(xiàn)有的基于二維各向同性介質(zhì)的彈性波有限差分波場數(shù)值模擬技術(shù)，進(jìn)而達(dá)到提高地震勘探技術(shù)精度的目的。目前基于二維各向同性介質(zhì)的彈性波有限差分波場數(shù)值模擬技術(shù)所用的計(jì)算速度模型角點(diǎn)波場的算法存在缺陷，這些算法不能有效消除來自任意角度角點(diǎn)的反射能量，而只能消除來自45度角(即角點(diǎn)的入射波場方向與模型的頂邊之間的夾角為45度角)的角點(diǎn)反射能量，從而導(dǎo)致計(jì)算得到的模型波場值不能正確反映波在無限介質(zhì)中傳播的真實(shí)性。

為了解決這個(gè)問題，在本發(fā)明第一發(fā)明人提出的彈性波吸收邊界條件的基礎(chǔ)上，本發(fā)明提出了彈性波廣角吸收邊界條件的算法，利用本發(fā)明的波場算法，無論激發(fā)源位于模型內(nèi)的什么位置，都可以很好地解決速度模型任意邊界角點(diǎn)的波場反射問題，使得到的模型波場值準(zhǔn)確地反映波在無限介質(zhì)中傳播的真實(shí)性，這樣便可以提高波場數(shù)值模擬的準(zhǔn)確性，從而達(dá)到提高地震勘探精度的目的。

本發(fā)明提供了一種基于二維各向同性介質(zhì)的彈性波波場數(shù)值的模擬方法，該方法包括以下步驟：

根據(jù)給定的速度模型，確定波場激發(fā)源的位置所對(duì)應(yīng)的坐標(biāo)值；

選擇地震子波，確定速度模型的空間采樣間隔和時(shí)間采樣間隔，對(duì)速度模型進(jìn)行有限差分網(wǎng)格離散化；

計(jì)算速度模型內(nèi)部各個(gè)網(wǎng)格點(diǎn)的各個(gè)時(shí)刻的波場值，計(jì)算速度模型四條邊上的網(wǎng)格點(diǎn)所對(duì)應(yīng)的各個(gè)時(shí)刻的波場值，計(jì)算速度模型四個(gè)角點(diǎn)所對(duì)應(yīng)的各個(gè)時(shí)刻的波場值；

由四個(gè)角點(diǎn)的波場值、四條邊上的網(wǎng)格點(diǎn)的波場值以及速度模型內(nèi)部各個(gè)網(wǎng)格點(diǎn)的波場值得到整個(gè)速度模型的各個(gè)網(wǎng)格點(diǎn)在各個(gè)時(shí)刻的波場值。

進(jìn)一步地，利用公式(3)計(jì)算速度模型內(nèi)部各個(gè)網(wǎng)格點(diǎn)的各個(gè)時(shí)刻的波場值：

其中：u、w分別為彈性波場所對(duì)應(yīng)的水平位移場和垂直位移場；

α、β分別為縱、橫波速度，x和z是二維介質(zhì)情況下沿x軸和z軸的坐標(biāo)。

進(jìn)一步地，利用公式(10a-10d)計(jì)算速度模型四個(gè)角點(diǎn)所對(duì)應(yīng)的各個(gè)時(shí)刻的波場值：

右下角角點(diǎn)的波場計(jì)算方程為：

左上角角點(diǎn)的波場計(jì)算方程為：

右上角角點(diǎn)的波場計(jì)算方程為：

左下角角點(diǎn)的波場計(jì)算方程為：

其中：u、w分別為彈性波場所對(duì)應(yīng)的水平位移場和垂直位移場；

α、β分別為縱、橫波速度，x和z是二維介質(zhì)情況下沿x軸和z軸的坐標(biāo)；

其中：

θ為波的入射角。

進(jìn)一步地，根據(jù)激發(fā)源的位置坐標(biāo)以及四個(gè)角點(diǎn)所對(duì)應(yīng)的位置坐標(biāo)，計(jì)算出sinθ和cosθ，由此得到c矩陣向量中的各項(xiàng)值。

進(jìn)一步地，利用公式(10a-10d)計(jì)算速度模型四條邊上的網(wǎng)格點(diǎn)所對(duì)應(yīng)的各個(gè)時(shí)刻的波場值：

速度模型底邊所對(duì)應(yīng)的吸收邊界方程為：

速度模型頂邊所對(duì)應(yīng)的吸收邊界方程為：

速度模型左邊所對(duì)應(yīng)的吸收邊界方程為：

速度模型右邊所對(duì)應(yīng)的吸收邊界方程為：

其中，式(11a)和(11b)所對(duì)應(yīng)的參數(shù)項(xiàng)為：

其中，式(11c)和(11d)所對(duì)應(yīng)的參數(shù)項(xiàng)為：

其中，α、β分別為縱、橫波速度，x和z是二維介質(zhì)情況下沿x軸和z軸的坐標(biāo)。

將本發(fā)明提出的彈性波廣角吸收邊界條件方法應(yīng)用于二維各向同性介質(zhì)的彈性波波場數(shù)值模擬，可以有效地消除來自有限的速度模型邊界及任意角度的邊界角點(diǎn)所產(chǎn)生的反射能量。本發(fā)明所提供的邊界波場以及邊界節(jié)點(diǎn)波場的算法能使往外傳播的波的能量透過邊界和邊界角點(diǎn)，使計(jì)算得到的模型波場值能準(zhǔn)確地反映出彈性波在二維各向同性無限介質(zhì)中傳播的真實(shí)性，從而達(dá)到提高地震勘探技術(shù)精度的目的。

附圖說明

通過結(jié)合附圖對(duì)本公開示例性實(shí)施方式進(jìn)行更詳細(xì)的描述，本公開的上述以及其它目的、特征和優(yōu)勢(shì)將變得更加明顯，其中，在本公開示例性實(shí)施方式中，相同的參考標(biāo)號(hào)通常代表相同部件。

圖1顯示了本發(fā)明基于二維各向同性介質(zhì)的彈性波波場數(shù)值的模擬方法的流程圖。

圖2顯示了根據(jù)一個(gè)實(shí)施例的二維各向同性介質(zhì)常速度模型。

圖3顯示了根據(jù)一個(gè)實(shí)施例的剛性邊界條件下t＝120ms時(shí)的波場傳播快照。

圖4顯示了根據(jù)一個(gè)實(shí)施例的吸收邊界條件下t＝120ms時(shí)的波場傳播快照。

圖5顯示了根據(jù)一個(gè)實(shí)施例的利用常規(guī)吸收邊界條件方法得到的t＝120ms時(shí)的波場傳播快照。

圖6顯示了根據(jù)一個(gè)實(shí)施例的利用本發(fā)明的方法得到的t＝120ms時(shí)的波場傳播快照。

具體實(shí)施方式

下面將參照附圖更詳細(xì)地描述本公開的優(yōu)選實(shí)施方式。雖然附圖中顯示了本公開的優(yōu)選實(shí)施方式，然而應(yīng)該理解，可以以各種形式實(shí)現(xiàn)本公開而不應(yīng)被這里闡述的實(shí)施方式所限制。相反，提供這些實(shí)施方式是為了使本公開更加透徹和完整，并且能夠?qū)⒈竟_的范圍完整地傳達(dá)給本領(lǐng)域的技術(shù)人員。

本發(fā)明提出的彈性波廣角吸收邊界條件算法是對(duì)現(xiàn)有的基于二維各向同性介質(zhì)彈性波傍軸近似理論的吸收邊界條件的一種改進(jìn)，它解決了來自速度模型內(nèi)任意角度的角點(diǎn)的反射能量問題，使通過彈性波數(shù)值模擬得到的模型波場值能正確地反映彈性波在二維各向同性無限介質(zhì)中傳播的真實(shí)性。

參照?qǐng)D1，描述本發(fā)明提供了一種基于二維各向同性介質(zhì)的彈性波波場數(shù)值的模擬方法，該方法包括以下步驟：

根據(jù)給定的速度模型，確定波場激發(fā)源的位置(即炮點(diǎn)位置)所對(duì)應(yīng)的坐標(biāo)值。優(yōu)選地，速度模型可以是各向同性彈性介質(zhì)的速度模型。

選擇地震子波，確定速度模型的空間采樣間隔和時(shí)間采樣間隔，對(duì)速度模型進(jìn)行有限差分網(wǎng)格離散化(子波一般選用雷克子波。優(yōu)選地，時(shí)間采樣間隔δt與空間采樣間隔δx、δz滿足如下的公式：其中：h＝min(δx,δz),α、β分別為縱、橫波速度)。

計(jì)算速度模型內(nèi)部各個(gè)網(wǎng)格點(diǎn)的各個(gè)時(shí)刻的波場值，計(jì)算速度模四條邊上的網(wǎng)格點(diǎn)所對(duì)應(yīng)的各個(gè)時(shí)刻的波場值，計(jì)算速度模型四個(gè)角點(diǎn)所對(duì)應(yīng)的各個(gè)時(shí)刻的波場值；

由四個(gè)角點(diǎn)的波場值、四條邊上的網(wǎng)格點(diǎn)的波場值以及速度模型內(nèi)部各個(gè)網(wǎng)格點(diǎn)的波場值得到整個(gè)速度模型的各個(gè)網(wǎng)格點(diǎn)在各個(gè)時(shí)刻的波場值。

根據(jù)本發(fā)明具體實(shí)施例，從二維各向同性介質(zhì)的彈性波傳播方程出發(fā)，通過系列推導(dǎo)得到了速度模型任意角度的角點(diǎn)的彈性波波場計(jì)算公式，其方法原理如下：

在各向同性介質(zhì)條件下，假設(shè)x和z是二維介質(zhì)情況下沿x軸和z軸的坐標(biāo)，x軸的方向向右，z軸的方向朝下，可以利用兩個(gè)共軛的二階偏微分方程來描述介質(zhì)中p波的運(yùn)動(dòng)以及垂向偏振的sv波的運(yùn)動(dòng)，在這里不考慮沿水平方向偏振的sh波，并假定u、w分別為水平方向與垂直方向的位移，ρ為介質(zhì)的密度，t為時(shí)間，λ和μ為具體介質(zhì)中的拉梅系數(shù)，那么可以得到關(guān)于各向同性非均勻介質(zhì)的彈性波動(dòng)方程如下式表示：

現(xiàn)在假定密度ρ為常數(shù)，這樣就可以把上面的方程組看成是隨空間位置變化的p波和sv波速度的函數(shù)。其中λ和μ與縱、橫波速度α(x,z)和β(x,z)的關(guān)系如下所示：

這樣便可以把(1)式寫成下列形式：

其中：分別為彈性波場所對(duì)應(yīng)的水平位移場和垂直位移場；

α、β分別為縱、橫波速度。

方程(3)的傅氏變換可以寫成下列形式：

為了能利用傍軸近似，將式(4)的左邊項(xiàng)寫成下面的形式：

d2(kz/ω)²+h(kz/ω)(kx/ω)+d1(kx/ω)²＝i(5)

對(duì)(5)式進(jìn)行整理，得：

令：

其中：

很顯然：

令：

其中：

θ為波的入射角；

根據(jù)激發(fā)源的位置坐標(biāo)以及四個(gè)角點(diǎn)所對(duì)應(yīng)的位置坐標(biāo)，計(jì)算出sinθ和cosθ，由此得到c矩陣向量中的各項(xiàng)值。

令：則有：

于是，(6)式可以寫成下面式子：

由(7)式便可以將(4)式改寫成下面的式子：

令d＝i+c(kz/ω)(kx/ω)+a(kz/ω)+b(kx/ω)，在(8)式的兩邊乘上d^-1，則得到下面式子：

將(9)式變換到時(shí)間域，可以得到下面的方程：

式(10a)為右下角角點(diǎn)的波場計(jì)算式，我們得到其他三個(gè)角點(diǎn)的波場計(jì)算式如下：

左上角角點(diǎn)的波場計(jì)算方程為：

右上角角點(diǎn)的波場計(jì)算方程為：

左下角角點(diǎn)的波場計(jì)算方程為：

在(10a-10d)各個(gè)表達(dá)式中，其中：u、w分別為彈性波場所對(duì)應(yīng)的水平位移場和垂直位移場；

α、β分別為縱、橫波速度。

本發(fā)明計(jì)算邊界的波場采用本發(fā)明第一發(fā)明人基于二維各向同性介質(zhì)彈性波傍軸近似理論提出的吸收邊界條件，其對(duì)應(yīng)的四邊方程如下所示：

速度模型底邊所對(duì)應(yīng)的吸收邊界方程為：

速度模型頂邊所對(duì)應(yīng)的吸收邊界方程為：

其中，式(11a)和(11b)所對(duì)應(yīng)的參數(shù)項(xiàng)為：

速度模型左邊所對(duì)應(yīng)的吸收邊界方程為：

速度模型右邊所對(duì)應(yīng)的吸收邊界方程為：

其中，式(11c)和(11d)所對(duì)應(yīng)的參數(shù)項(xiàng)為：

為便于理解本發(fā)明實(shí)施例的方案及其效果，以下給出一個(gè)具體應(yīng)用示例。本領(lǐng)域技術(shù)人員應(yīng)理解，該示例僅為了便于理解本發(fā)明，其任何具體細(xì)節(jié)并非意在以任何方式限制本發(fā)明。下面結(jié)合附圖和實(shí)施例對(duì)本發(fā)明進(jìn)一步說明。

實(shí)現(xiàn)彈性波波場正演，首先需要確定速度模型，圖2為本實(shí)施例所選的二維各向同性介質(zhì)常速度模型，縱波速度為3000米/秒，橫波速度為2000米/秒，激發(fā)源為縱波震源，位于速度模型的中央。

在正演過程中，選取主頻為50hz的雷克子波作為本實(shí)施例的地震子波，空間采樣間隔為1米，時(shí)間采樣間隔為0.25ms，然后根據(jù)所選的空間采樣對(duì)所選的速度模型進(jìn)行網(wǎng)格離散化。

圖3為選用剛性邊界條件(即波場傳播到邊界時(shí)出現(xiàn)全反射的情況)波場傳播到120毫秒時(shí)整個(gè)模型的波場快照；圖4為使用本發(fā)明的彈性波廣角吸收邊界條件，波場傳播到120毫秒時(shí)刻時(shí)整個(gè)模型的波場快照。對(duì)比圖3和圖4，很容易看出，本發(fā)明所使用的吸收邊界條件可以對(duì)傳播到邊界的波場起到很好的吸收作用。

為了展示本發(fā)明在波場正演中的優(yōu)越性，我們將模型1中的激發(fā)源位置進(jìn)行一下變動(dòng)：將激發(fā)源的位置移到距模型頂邊12.5米、距模型的左邊102.5米的位置。圖5為使用常規(guī)吸收邊界條件所得到的波場傳播到120毫秒時(shí)刻時(shí)整個(gè)模型的波場快照；圖6為使用本發(fā)明的廣角吸收邊界條件所得到的波場傳播到120毫秒時(shí)刻時(shí)整個(gè)模型的波場快照。

對(duì)比圖5和圖6，不難發(fā)現(xiàn)常規(guī)吸收邊界條件不能對(duì)來自角點(diǎn)的波場反射起到很好的吸收作用，而使用本發(fā)明的廣角吸收邊界條件則可以對(duì)來自角點(diǎn)的波場反射起到較好的吸收作用。

本發(fā)明提出了彈性波廣角吸收邊界條件波場算法，利用本發(fā)明的波場算法可以有效地消除來自邊界及任意角度的邊界角點(diǎn)的反射能量，使計(jì)算的模型波場值能準(zhǔn)確地反映波在無限介質(zhì)中傳播的真實(shí)性，從而達(dá)到提高勘探精度的目的。

以上已經(jīng)描述了本公開的各實(shí)施例，上述說明是示例性的，并非窮盡性的，并且也不限于所披露的各實(shí)施例。在不偏離所說明的各實(shí)施例的范圍和精神的情況下，對(duì)于本技術(shù)領(lǐng)域的普通技術(shù)人員來說許多修改和變更都是顯而易見的。本文中所用術(shù)語的選擇，旨在最好地解釋各實(shí)施例的原理、實(shí)際應(yīng)用或?qū)κ袌鲋械募夹g(shù)改進(jìn)，或者使本技術(shù)領(lǐng)域的其它普通技術(shù)人員能理解本文披露的各實(shí)施例。