本發(fā)明涉及石油地震勘探技術(shù)領(lǐng)域，特別涉及一種基于彈簧網(wǎng)絡(luò)模型的正演模擬方法及裝置。

背景技術(shù)：

地震正演模擬是地震資料采集、處理和地質(zhì)解釋的基礎(chǔ)，是降低地震多解性的重要手段。目前，地震波數(shù)值模擬的主流方法都是通過近似求解波動方程來實現(xiàn)的，比如常用的有限差分法、有限元法等，其固有缺點是必須假設(shè)地質(zhì)體是相對均勻而連續(xù)的單元，以確保描述地震波傳播的偏微分方程具有可解性。然而，隨著地質(zhì)體非均勻性的增強，當介質(zhì)中存在強間斷面或者含有多相流體時，這一問題變得更加突出。因此，怎樣將實際復(fù)雜的地下孔隙介質(zhì)，用簡單理想化且不失一般性的模型來描述，一直是地震波場正演模擬研究的重要課題與方向。

彈簧網(wǎng)絡(luò)模型作為一種微觀力學(xué)模型，通過與速度Verlet結(jié)合可用于模擬彈性波的傳播過程。由于該方法是從最基本的牛頓運動定律和胡克定律出發(fā)來進行波場模擬，不會受到波動方程的近似假設(shè)條件限制，因而能夠適應(yīng)任意復(fù)雜的地下介質(zhì)模型。

彈簧網(wǎng)絡(luò)模型(Lattice Spring Model，簡稱LSM)是一種從微觀角度研究介質(zhì)彈塑性的方法，其相關(guān)原理可追溯到20世紀初Born等學(xué)者提出的晶體動力學(xué)理論。早期的Born彈簧模型，相鄰質(zhì)點之間只考慮中心力(Central-Force)作用，這種力也稱為線彈簧作用力。這種模型模擬的介質(zhì)泊松比固定為0.25，其應(yīng)用受到限制。1939年，Kirkwood在Born彈簧模型的基礎(chǔ)上引入了一種鍵彎曲作用力(Bond-Bending Force)，即角彈簧作用力，使得模擬材料的泊松比能夠在一定范圍內(nèi)變化，該模型因此被稱為鍵彎曲模型(Bond-Bending Model)。1954年，Born和Huang合著的圖書《Dynamical Theory of Crystal Lattices》比較全面地介紹了這樣一種微觀模型，為該模型用于研究晶體的彈性、內(nèi)聚能和晶格動力學(xué)等打下了堅實基礎(chǔ)。之后，Kittel、Lax和Keating等人對Born彈簧模型作了進一步地發(fā)展與完善，保證了形變后晶體的應(yīng)變能量計算公式具有旋轉(zhuǎn)不變性。

Born彈簧模型在早期主要用于研究介質(zhì)的彈性性能與斷裂行為。1984年，美國學(xué)者Grest和Webman提出了將固體離散成一系列由彈簧與節(jié)點連接成的立方體元的想法，并將其用于研究滲透集群的震動模式。Kantor和Webman(1984)、Arbabi和Sahimi(1988)先后通過將滲流系統(tǒng)分解成微觀的鍵彎曲模型，研究了其彈性性能。Hassold和Srolovitz(1989)、Murat等(1992)通過Born彈簧模型模擬得到了介質(zhì)的破裂統(tǒng)計特性。Starzewski等(1996)采用Spring Network Model研究了復(fù)合材料和晶體的伸縮性和斷裂性。Ladd等(1997)進一步研究了線彈簧與角彈簧的綜合模型，并最先使用了“Lattice Spring Model”一詞。Buxton等(2001)，Chung等(2002)采用LSM描述了非均勻介質(zhì)的彈塑性形變和脆性破裂效應(yīng)。

2000年，愛爾蘭的兩位學(xué)者Toomey和Bean采用一種離散粒子模型，通過將介質(zhì)離散成微觀球體緊密排列的正六邊形構(gòu)造單元，借助虎克定律和速度Verlet算法，實現(xiàn)了泊松體中彈性波的模擬，該方法被稱為離散粒子方法(Discrete Particle Scheme，簡稱DPS)，也有學(xué)者稱其為彈性網(wǎng)格方法(Elastic Lattice Method，簡稱ELM)。其后，O’Brien和Bean對相關(guān)理論作了進一步推廣，解除了泊松比的限制，并在火山震源機制、彈性波模擬等方面得到了應(yīng)用。

此外，Yim和Sohn(2000)構(gòu)造了一種新的質(zhì)點彈簧模型(Mass Spring Lattice Model，簡稱MSLM)，將傳統(tǒng)的彈性波方程中的彈性常數(shù)用彈簧的彈性常數(shù)替換，采用有限差分方法實現(xiàn)了復(fù)雜介質(zhì)中的超聲波模擬。2008年，O’Brien通過在LSM中加入粘彈性彈簧模擬了地震波在粘彈性介質(zhì)中的傳播。Pazdniakou和Adler(2012)對LSM模型的相關(guān)理論作了比較全面的介紹，并嘗試將其用于模擬干燥孔隙介質(zhì)中的彈性波。Zhao等(2011)在原始LSM的基礎(chǔ)上增加了一種獨特的剪切彈簧并應(yīng)用Distinct Lattice Spring Model(DLSM)研究了介質(zhì)的彈性和動態(tài)破裂性質(zhì)。Zhu等(2011)借助DLSM模擬了P波在不同巖石交界處的傳播過程。2012年，Zhao等對DLSM模型的演化作了簡要的概述，并提出了一種采用高階DLSM模擬介質(zhì)彈性的方法。同年，Zhao和Khalili(2012)對DLSM的GPU并行加速算法進行了研究。

綜上所述，目前基于彈簧網(wǎng)絡(luò)模型或者類似方法的正演模擬手段，主要集中于正方形網(wǎng)格或者正三角形網(wǎng)格情形，而在實際生產(chǎn)中，為了適應(yīng)特定的觀測系統(tǒng)，常常需要將網(wǎng)格剖分成矩形或長方體，以節(jié)省工作量以滿足特定方向的計算精度，但現(xiàn)有技術(shù)中尚欠缺適應(yīng)矩形網(wǎng)格或長方體網(wǎng)格的問題研究。

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

本發(fā)明實施例提供了一種基于彈簧網(wǎng)絡(luò)模型的正演模擬方法，將常規(guī)的僅僅適用于正方形網(wǎng)格或者正三角形網(wǎng)格的彈簧網(wǎng)絡(luò)模型，推廣到了適用于矩形網(wǎng)格或長方體網(wǎng)格，能夠更好地適應(yīng)實際生產(chǎn)需要，該方法包括：

設(shè)定正演模擬采用的時間采樣間隔和空間采樣間隔；

根據(jù)所述空間采樣間隔對待模擬的地質(zhì)模型進行離散化，獲得數(shù)字化地質(zhì)模型；所述數(shù)字化地質(zhì)模型為矩形網(wǎng)格形式或長方體網(wǎng)格形式；

根據(jù)所述空間采樣間隔確定所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)；

根據(jù)震源子波函數(shù)、時間采樣間隔及所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)，更新所述數(shù)字化地質(zhì)模型的各個時刻的波場值。

本發(fā)明實施例提供了一種基于彈簧網(wǎng)絡(luò)模型的正演模擬裝置，將常規(guī)的僅僅適用于正方形網(wǎng)格或者正三角形網(wǎng)格的彈簧網(wǎng)絡(luò)模型，推廣到了適用于矩形網(wǎng)格或長方體網(wǎng)格，能夠更好地適應(yīng)實際生產(chǎn)需要，該裝置包括：

參數(shù)設(shè)定模塊，用于設(shè)定正演模擬采用的時間采樣間隔和空間采樣間隔；

離散模塊，用于根據(jù)所述空間采樣間隔對待模擬的地質(zhì)模型進行離散化，獲得數(shù)字化地質(zhì)模型；所述數(shù)字化地質(zhì)模型為矩形網(wǎng)格形式或長方體網(wǎng)格形式；

彈性參數(shù)確定模塊，用于根據(jù)所述空間采樣間隔確定所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)；

波場值確定模塊，用于根據(jù)震源子波函數(shù)、時間采樣間隔及所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)，更新所述數(shù)字化地質(zhì)模型的各個時刻的波場值。

在本發(fā)明實施例中，設(shè)定正演模擬采用的時間采樣間隔和空間采樣間隔；根據(jù)所述空間采樣間隔對待模擬的地質(zhì)模型進行離散化，獲得數(shù)字化地質(zhì)模型；所述數(shù)字化地質(zhì)模型為矩形網(wǎng)格形式或長方體網(wǎng)格形式；根據(jù)所述空間采樣間隔確定所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)；根據(jù)震源子波函數(shù)、時間采樣間隔及所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)，更新所述數(shù)字化地質(zhì)模型的各個時刻的波場值。采用本發(fā)明方法可以將常規(guī)的僅僅適用于正方形網(wǎng)格或者正三角形網(wǎng)格的彈簧網(wǎng)絡(luò)模型，推廣到了適用于更一般性的矩形網(wǎng)格或長方體網(wǎng)格，這樣能夠更好地與實際野外地震資料采集觀測系統(tǒng)采集的數(shù)據(jù)進行匹配，滿足了實際生產(chǎn)中對于不同方向的計算精度以及計算效率的需求。

附圖說明

為了更清楚地說明本發(fā)明實施例或現(xiàn)有技術(shù)中的技術(shù)方案，下面將對實施例或現(xiàn)有技術(shù)描述中所需要使用的附圖作簡單地介紹，顯而易見地，下面描述中的附圖僅僅是本發(fā)明的一些實施例，對于本領(lǐng)域普通技術(shù)人員來講，在不付出創(chuàng)造性勞動的前提下，還可以根據(jù)這些附圖獲得其他的附圖。

圖1為本發(fā)明實施例提供的一種基于彈簧網(wǎng)絡(luò)模型的正演模擬方法的處理流程圖；

圖2為本發(fā)明實施例提供的一種基于彈簧網(wǎng)絡(luò)模型的正演模擬方法具體處理流程圖；

圖3為本發(fā)明實施例提供的二維(D2Q8)彈簧網(wǎng)絡(luò)模型示意圖；

圖4為本發(fā)明實施例提供的三維(D3Q18)彈簧網(wǎng)絡(luò)模型示意圖；

圖5a)為本發(fā)明實施例提供的LSM與FDM在網(wǎng)格比(dz/dx＝5/8)下的P波的歸一化相速度頻散曲線對比圖；

圖5b)為本發(fā)明實施例提供的LSM與FDM在網(wǎng)格比(dz/dx＝5/8)下的S波的歸一化相速度頻散曲線對比圖；

圖5c)為本發(fā)明實施例提供的LSM與FDM在網(wǎng)格比(dz/dx＝1)下的P波的歸一化相速度頻散曲線對比圖；

圖5d)為本發(fā)明實施例提供的LSM與FDM在網(wǎng)格比(dz/dx＝1)下的S波的歸一化相速度頻散曲線對比圖；

圖5e)為本發(fā)明實施例提供的LSM與FDM在網(wǎng)格比(dz/dx＝8/5)下的P波的歸一化相速度頻散曲線對比圖；

圖5f)為本發(fā)明實施例提供的LSM與FDM在網(wǎng)格比(dz/dx＝8/5)下的S波的歸一化相速度頻散曲線對比圖；

圖6a)為本發(fā)明實施例提供的二階空間精度的FDM模擬得到的波場(v_x)切片圖；

圖6b)為本發(fā)明實施例提供的二階空間精度的FDM模擬得到的波場(v_z)切片圖；

圖6c)為本發(fā)明實施例提供的四階空間精度的FDM模擬得到的波場(v_x)切片圖；

圖6d)為本發(fā)明實施例提供的四階空間精度的FDM模擬得到的波場(v_z)切片圖；

圖6e)為本發(fā)明實施例提供的LSM模擬得到的波場(v_x)切片圖；

圖6f)為本發(fā)明實施例提供的LSM模擬得到的波場(v_z)切片圖；

圖7a)為本發(fā)明實施例提供的二維層狀介質(zhì)情形下的界面之上的接收點處的LSM與二階及四階精度FDM的單道地震記錄對比圖(v_x)；

圖7b)為本發(fā)明實施例提供的二維層狀介質(zhì)情形下的界面之下的接收點處的LSM與二階及四階精度FDM的單道地震記錄對比圖(v_x)；

圖7c)為本發(fā)明實施例提供的二維層狀介質(zhì)情形下的界面之上的接收點處的LSM與二階及四階精度FDM的單道地震記錄對比圖(v_z)；

圖7d)為本發(fā)明實施例提供的二維層狀介質(zhì)情形下的界面之下的接收點處的LSM與二階及四階精度FDM的單道地震記錄對比圖(v_z)；

圖8a)為本發(fā)明實施例提供的Marmousi模型的FDM模擬得到的波場(v_x)快照圖；

圖8b)為本發(fā)明實施例提供的Marmousi模型的FDM模擬得到的波場(v_z)快照圖；

圖8c)為本發(fā)明實施例提供的Marmousi模型的LSM模擬得到的波場(v_x)快照圖；

圖8d)為本發(fā)明實施例提供的Marmousi模型的LSM模擬得到的波場(v_z)快照圖；

圖9為本發(fā)明實施例提供的一種基于彈簧網(wǎng)絡(luò)模型的正演模擬裝置的結(jié)構(gòu)示意圖；

圖10為本發(fā)明實施例提供的一種基于彈簧網(wǎng)絡(luò)模型的正演模擬裝置的一種具體結(jié)構(gòu)示意圖。

具體實施方式

下面將結(jié)合本發(fā)明實施例中的附圖，對本發(fā)明實施例中的技術(shù)方案進行清楚、完整地描述，顯然，所描述的實施例僅是本發(fā)明一部分實施例，而不是全部的實施例。基于本發(fā)明中的實施例，本領(lǐng)域普通技術(shù)人員在沒有做出創(chuàng)造性勞動前提下所獲得的所有其他實施例，都屬于本發(fā)明保護的范圍。

現(xiàn)有的基于彈簧網(wǎng)絡(luò)模型或者類似原理的地震波場正演模擬方法，幾乎全都是采用正方形或者正三角形網(wǎng)格，這樣的方法在某種意義上不能滿足實際野外地震資料采集的需求?；诖?，本發(fā)明提出一種基于彈簧網(wǎng)絡(luò)模型的正演模擬方法，以彌補前人提出方法的不足之處，從而更好地與實際地震資料進行匹配，并達到提高數(shù)值計算效率的目的。

圖1是本發(fā)明實施例提供的一種基于彈簧網(wǎng)絡(luò)模型的正演模擬方法的處理流程圖，如圖1所示，該方法包括：

步驟101、設(shè)定正演模擬采用的時間采樣間隔和空間采樣間隔。

步驟102、根據(jù)所述空間采樣間隔對待模擬的地質(zhì)模型進行離散化，獲得數(shù)字化地質(zhì)模型；所述數(shù)字化地質(zhì)模型為矩形網(wǎng)格形式或長方體網(wǎng)格形式。

步驟103、根據(jù)所述空間采樣間隔確定所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)。

步驟104、根據(jù)震源子波函數(shù)、時間采樣間隔及所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)，更新所述數(shù)字化地質(zhì)模型的各個時刻的波場值。

具體實施時，在得到一個地質(zhì)模型對其進行正演模擬時，首先需要設(shè)定正演模擬所要采用的時間采樣間隔和空間采樣間隔(步驟101)，然后需要根據(jù)預(yù)設(shè)的空間采樣間隔對待模擬的地質(zhì)模型(模擬地質(zhì)區(qū)域)進行離散化，獲得數(shù)字化地質(zhì)模型(步驟102)。進行離散化的方法是：對于二維(2D)數(shù)值模擬問題，將待模擬的地質(zhì)模型剖分為矩形網(wǎng)格(D2Q8模型，見圖3)，或者，對于三維(3D)數(shù)值模擬問題，將待模擬的地質(zhì)模型剖分為長方體網(wǎng)格(D3Q18模型，見圖4)，二維及三維模型的網(wǎng)格邊長可以相等，也可以不相等。然后將待模擬的地質(zhì)模型離散成由離散網(wǎng)格節(jié)點與線彈簧構(gòu)成的彈性元，相鄰離散網(wǎng)格節(jié)點之間通過線彈簧連接；相鄰離散網(wǎng)格節(jié)點之間的距離為預(yù)設(shè)的空間采樣間隔。矩形網(wǎng)格(D2Q8模型)形式的數(shù)字化地質(zhì)模型由彈性面元構(gòu)成；長方體網(wǎng)格(D3Q18模型)形式的數(shù)字化地質(zhì)模型由彈性體元構(gòu)成。

當然，本發(fā)明方法對正方形網(wǎng)格或者正方體網(wǎng)格同樣適用。

具體實施時，在獲得數(shù)字化地質(zhì)模型后，首先需要根據(jù)預(yù)設(shè)的空間采樣間隔求取數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)。假設(shè)，對數(shù)字化地質(zhì)模型(離散模型)施加外力或者擾動后，相應(yīng)位置的線彈簧會發(fā)生形變，彈性面元(2D)或者彈性體元(3D)的中心節(jié)點i受到的線彈簧的總作用力F_i可以按照如下公式計算：

其中，φ_i為第i個彈性元的能量密度，k_j為中心節(jié)點i與相鄰節(jié)點j之間的線彈簧的彈性常數(shù)，u_i為節(jié)點的位移，x_ij為節(jié)點i指向節(jié)點j的向量，為歸一化的方向向量，∑代表對中心節(jié)點的所有鄰點求和，符號“·”代表向量內(nèi)積運算；

對于二維D2Q8模型，上式中，E為彈性面元的面積，N＝8；

對于三維D3Q18模型，上式中，E為彈性體元的體積，N＝18。

線彈簧的彈性參數(shù)是一個與彈簧長度有關(guān)的常數(shù)，在離散模型中具體表現(xiàn)為與鄰點位置或者方向有關(guān)的常數(shù)，在計算彈性體元或者彈性面元時，用到了鄰近的8個節(jié)點(D2Q8模型)或者18個節(jié)點(D3Q18模型)上的波場值，且不同鄰點與中心節(jié)點之間的線彈簧的彈性常數(shù)不相等。假設(shè)相等長度的線彈簧具有相同的彈性參數(shù)時，其計算公式如下：

①對于二維D2Q8模型，各個方向彈性參數(shù)的計算公式如下：

且有，

其中，Δx和Δz分別為沿著x軸和z軸方向的網(wǎng)格長度，ρ為介質(zhì)的質(zhì)量密度，V_P為縱波傳播速度，k₁₀、k₂₀和k₃₀分別為以中心節(jié)點為頂點，沿著x軸正負方向、z軸正負方向以及對角正負方向的線彈簧彈性參數(shù)，如圖3所示。

②對于三維D3Q18模型，各個方向彈性常數(shù)的計算公式如下：

且有，

其中，Δx、Δy和Δz分別為沿著x軸、y軸和z軸方向的網(wǎng)格長度，ρ為介質(zhì)的質(zhì)量密度，V_P為縱波傳播速度，k₁₀₀、k₂₀₀和k₃₀₀分別為以中心節(jié)點為頂點，沿著x軸、y軸和z軸方向的線彈簧彈性參數(shù)，k₄₀₀、k₅₀₀和k₆₀₀分別代表沿著xoy、xoz和yoz坐標平面對角方向的線彈簧彈性參數(shù)，如圖4所示。

具體實施時，計算完模擬區(qū)域內(nèi)各個節(jié)點的總的作用力后，根據(jù)震源子波函數(shù)、時間采樣間隔及數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)，利用速度Verlet算法更新數(shù)字化地質(zhì)模型的各個時刻的波場值(步驟104)。該方法不是基于有限差分方法(FDM)或者有限元法(FEM)中常用的彈性波動方程，而是依據(jù)最基本的牛頓運動定律和胡克定律來模擬彈性波在固體介質(zhì)中的傳播，具體實施過程中，引用流體力學(xué)中經(jīng)常用到的速度Verlet算法，用于更新模擬區(qū)域內(nèi)各個時刻的波場值，波場值包括位移x_i、速度v_i和加速度a_i：

其中，Δt為模擬采用的時間采樣間隔，F(xiàn)_i為中心節(jié)點i所受到的線彈簧的總作用力，m_i為中心節(jié)點i的質(zhì)量，χ為粘滯項系數(shù)，用于衰減節(jié)點的震動能量，這里取χ＝0。

上述速度Verlet算法是一個時間二階計算精度的數(shù)值算法，而彈簧網(wǎng)絡(luò)模型在數(shù)值模擬時的空間計算精度取決于中心節(jié)點i所受作用力的計算精度。

具體實施時，如圖2所示，該方法還包括：

步驟105：根據(jù)所述時間采樣間隔和空間采樣間隔確定所述數(shù)字化地質(zhì)模型的相速度頻散曲線。

具體的，在二維數(shù)值模擬情況下，針對一般性的矩形網(wǎng)格(Δx≠Δz)，數(shù)字化地質(zhì)模型的P波和S波的相速度頻散曲線按照如下公式確定：

中間變量計算公式如下：

上式中

且有

其中，q_p為LSM模擬的波場對應(yīng)的P波的相速度頻散曲線，q_s為LSM模擬的波場對應(yīng)的S波的相速度頻散曲線，r＝Δz/Δx為z軸方向與x軸方向上的網(wǎng)格邊長之比，Δt為時間采樣間隔，V_P和V_S分別為縱波和橫波速度，θ為平面波與x軸正方向的夾角，λ為波數(shù)。

需要指出的是，上述針對一般性矩形網(wǎng)格的相速度頻散曲線，對于正方形網(wǎng)格(Δx＝Δz)情形同樣適用。

倘若僅僅考慮線彈簧作用力時，縱橫波速度之比(V_P/V_S)為也就是模擬介質(zhì)的泊松比在三維情況下為0.25。

具體實施時，如圖2所示，該方法還包括：

步驟106：對所述時間采樣間隔和空間采樣間隔進行驗證，確定時間采樣間隔和空間采樣間隔的選擇是否合理。具體的，根據(jù)彈簧網(wǎng)絡(luò)模型的穩(wěn)定性條件和相速度頻散曲線，驗證時間采樣間隔和空間采樣間隔是否合理。當所述時間采樣間隔不滿足彈簧網(wǎng)絡(luò)模型的穩(wěn)定性條件，且所述空間采樣間隔不滿足相速度頻散精度要求時，重新設(shè)定時間采樣間隔和空間采樣間隔，對重新設(shè)定的時間采樣間隔和空間采樣間隔進行再次驗證，直到重新設(shè)定的時間采樣間隔滿足彈簧網(wǎng)絡(luò)模型的穩(wěn)定性條件，且重新設(shè)定的空間采樣間隔滿足相速度頻散精度要求為止。

比如，驗證時間采樣間隔和空間采樣間隔選擇過大或過小，則重新設(shè)定時將時間采樣間隔和空間采樣間隔減小或增大，然后再次驗證，直到時間采樣間隔和空間采樣間隔滿足數(shù)值計算穩(wěn)定性條件和數(shù)值頻散精度要求。

在實際數(shù)值模擬過程中，彈簧網(wǎng)絡(luò)模型(LSM)的穩(wěn)定性條件為：

Δt＜Δd_min/V_max (14)

其中，Δt為時間采樣間隔，Δd_min為最小的空間采樣間隔，V_max為最大傳播速度。

數(shù)值模擬中，該方法(LSM)的穩(wěn)定性條件(Δt＜Δd_min/V_max)比FDM(對于2D情形：對于3D情形：更寬松，與FEM的穩(wěn)定性條件(Δt＜Δd_min/V_max)在同一個水平。

最后，輸出并保存計算獲得的各個時刻的波場值。

下面通過比較來說明本方法的優(yōu)點。

圖5a)至5f)顯示的是LSM與時間二階、空間二階精度的FDM在不同網(wǎng)格比下的P波或S波的歸一化相速度頻散曲線對比圖。其中，不同網(wǎng)格比分別為：dz/dx＝5/8、dz/dx＝1、dz/dx＝8/5。從圖5c)和5d)中可看出，當采用正方形網(wǎng)格(Δx＝Δz)時，兩種方法的頻散曲線很接近；從圖5a)和5b)、5e)和5f)中可看出，而當采用矩形網(wǎng)格(Δx≠Δz即dz/dx＝5/8或dz/dx＝8/5)時，兩種方法的數(shù)值頻散特性各有優(yōu)劣，總體差別不大。

接下來通過一個二維雙層介質(zhì)模型來驗證本發(fā)明方法的有效性。模型上層與下層縱波速度分別為4000m/s與5000m/s，而橫波速度取時間采樣間隔Δt＝0.5ms；x軸方向與z軸方向的空間采樣間隔取值為：Δx＝5m,Δz＝8m；模型大小為：5600m×5000m；震源采用主頻為15Hz雷克子波。針對相同的模型參數(shù)，采用LSM、二階(FDM-2nd)與四階(FDM-4th)空間精度的FDM模擬得到的波場切片對比圖見圖6a)至6f)，相應(yīng)的單道地震記錄見圖7a)至7d)。

根據(jù)圖6a)至6f)和圖7a)至7d)，不難看出，對于二維層狀介質(zhì)模型，LSM模擬得到的地震波場切片和單道地震記錄與二階及四階精度的FDM吻合很好，兩者的誤差很小，驗證了該方法的正確性。

為了更進一步說明本方法對實際復(fù)雜介質(zhì)模型的適用性，我們采用Marmousi模型進行實驗測試。模擬采用的時間采樣間隔為0.5ms；x軸方向與z軸方向的空間采樣間隔取值為：Δx＝5m,Δz＝7m；15Hz主頻的雷子子波震源加載在模型中間，其坐標為(1915m，700m)。為了作對比，采用相同參數(shù)的FDM也進行了相應(yīng)的波場模擬。其中，為了減弱模擬區(qū)域的邊界反射波，F(xiàn)DM采用的是CPML吸收邊界，而LSM采用的是傳統(tǒng)的指數(shù)衰減吸收邊界(ABC)，模擬得到的波場快照見圖8a)至8d)。通過對比，可發(fā)現(xiàn)兩種方法模擬得到的復(fù)雜介質(zhì)模型的波場快照吻合，這說明了LSM能夠適應(yīng)實際復(fù)雜地質(zhì)模型的彈性波正演工作。

基于同一發(fā)明構(gòu)思，本發(fā)明實施例中還提供了一種基于彈簧網(wǎng)絡(luò)模型的正演模擬裝置，如下面的實施例所述。由于基于彈簧網(wǎng)絡(luò)模型的正演模擬裝置解決問題的原理與基于彈簧網(wǎng)絡(luò)模型的正演模擬方法相似，因此基于彈簧網(wǎng)絡(luò)模型的正演模擬裝置的實施可以參見基于彈簧網(wǎng)絡(luò)模型的正演模擬方法的實施，重復(fù)之處不再贅述。以下所使用的，術(shù)語“單元”或者“模塊”可以實現(xiàn)預(yù)定功能的軟件和/或硬件的組合。盡管以下實施例所描述的裝置較佳地以軟件來實現(xiàn)，但是硬件，或者軟件和硬件的組合的實現(xiàn)也是可能并被構(gòu)想的。

圖9是本發(fā)明實施例的基于彈簧網(wǎng)絡(luò)模型的正演模擬裝置的一種結(jié)構(gòu)框圖，如圖9所示，包括：

參數(shù)設(shè)定模塊901，用于設(shè)定正演模擬采用的時間采樣間隔和空間采樣間隔；

離散模塊902，用于根據(jù)所述空間采樣間隔對待模擬的地質(zhì)模型進行離散化，獲得數(shù)字化地質(zhì)模型；所述數(shù)字化地質(zhì)模型為矩形網(wǎng)格形式或長方體網(wǎng)格形式；

彈性參數(shù)確定模塊903，用于根據(jù)所述空間采樣間隔確定所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)；

波場值確定模塊904，用于根據(jù)震源子波函數(shù)、時間采樣間隔及所述數(shù)字化地質(zhì)模型的各個方向的彈性參數(shù)，更新所述數(shù)字化地質(zhì)模型的各個時刻的波場值。

下面對該結(jié)構(gòu)進行說明。

具體實施時，離散模塊902具體用于：

將待模擬的地質(zhì)模型離散成由離散網(wǎng)格節(jié)點與線彈簧構(gòu)成的彈性元，相鄰離散網(wǎng)格節(jié)點之間通過線彈簧連接；相鄰離散網(wǎng)格節(jié)點之間的距離為所述空間采樣間隔；

矩形網(wǎng)格形式的數(shù)字化地質(zhì)模型由彈性面元構(gòu)成；

長方體網(wǎng)格形式的數(shù)字化地質(zhì)模型由彈性體元構(gòu)成。

具體實施時，所述彈性參數(shù)確定模塊903具體用于：

對所述數(shù)字化地質(zhì)模型施加外力或者擾動，相應(yīng)位置的線彈簧發(fā)生形變，彈性元的中心節(jié)點i受到的線彈簧的總作用力F_i按照如下公式確定：

其中，φ_i為第i個彈性元的能量密度，k_j為中心節(jié)點i與相鄰節(jié)點j之間的線彈簧的彈性常數(shù)，u_i為中心節(jié)點i的位移，u_j為相鄰節(jié)點j的位移，x_ij為中心節(jié)點i指向相鄰節(jié)點j的向量，為中心節(jié)點i指向相鄰節(jié)點j的歸一化的方向向量，∑代表對中心節(jié)點i的所有相鄰節(jié)點j求和，符號“·”代表向量內(nèi)積運算，N為相鄰節(jié)點j的個數(shù)；中心節(jié)點的個數(shù)和彈性元的個數(shù)相同；

對于矩形網(wǎng)格形式的數(shù)字化地質(zhì)模型，E為彈性面元的面積，N＝8；

對于長方體網(wǎng)格形式的數(shù)字化地質(zhì)模型，E為彈性體元的體積，N＝18；

對于矩形網(wǎng)格形式的數(shù)字化地質(zhì)模型，各個方向的彈性參數(shù)按如下公式確定：

其中，Δx和Δz分別為沿著x軸和z軸方向的網(wǎng)格長度，ρ為介質(zhì)的質(zhì)量密度，V_P為縱波傳播速度，k₁₀、k₂₀和k₃₀分別為以中心節(jié)點為頂點，沿著x軸方向、z軸方向以及對角方向的線彈簧彈性常數(shù)；

對于長方體網(wǎng)格形式的數(shù)字化地質(zhì)模型，各個方向的彈性參數(shù)按如下公式確定：

其中，Δx、Δy和Δz分別為沿著x軸、y軸和z軸方向的網(wǎng)格長度，ρ為介質(zhì)的質(zhì)量密度，V_P為縱波傳播速度，k₁₀₀、k₂₀₀和k₃₀₀分別為以中心節(jié)點為頂點，沿著x軸、y軸和z軸方向的線彈簧彈性常數(shù)，k₄₀₀、k₅₀₀和k₆₀₀分別代表沿著xoy、xoz和yoz坐標平面對角方向的線彈簧彈性常數(shù)。

具體實施時，所述波場值確定模塊904具體用于：

利用速度Verlet算法按如下公式更新所述數(shù)字化地質(zhì)模型的各個時刻的波場值：

其中，x_i為位移；v_i為速度；a_i為加速度；Δt為時間采樣間隔，F(xiàn)_i為中心節(jié)點i所受到的線彈簧的總作用力，m_i為中心節(jié)點i的質(zhì)量，χ為粘滯項系數(shù)，χ＝0。

具體實施時，如圖10所示，該裝置還包括：相速度數(shù)值頻散分析模塊905，用于根據(jù)所述時間采樣間隔和空間采樣間隔確定所述數(shù)字化地質(zhì)模型的相速度頻散曲線；

所述相速度數(shù)值頻散分析模塊905具體用于：

按如下公式確定離散格式的彈簧網(wǎng)絡(luò)模型的相速度頻散曲線按照如下公式確定：

中間變量計算公式如下：

上式中

且有

其中，q_p為LSM模擬的波場對應(yīng)的P波的相速度頻散曲線，q_s為LSM模擬的波場對應(yīng)的S波的相速度頻散曲線，r＝Δz/Δx為z軸方向與x軸方向上的網(wǎng)格邊長之比，Δt為時間采樣間隔，V_P和V_S分別為縱波和橫波速度，θ為平面波與x軸正方向的夾角，λ為波數(shù)。

具體實施時，如圖10所示，該裝置還包括：驗證模塊906，用于對所述時間采樣間隔和空間采樣間隔進行驗證，當所述時間采樣間隔不滿足彈簧網(wǎng)絡(luò)模型的穩(wěn)定性條件，且所述空間采樣間隔不滿足相速度頻散精度要求時，重新設(shè)定時間采樣間隔和空間采樣間隔，對重新設(shè)定的時間采樣間隔和空間采樣間隔進行再次驗證，直到重新設(shè)定的時間采樣間隔滿足彈簧網(wǎng)絡(luò)模型的穩(wěn)定性條件，且重新設(shè)定的空間采樣間隔滿足相速度頻散精度要求為止。

具體實施時，所述彈簧網(wǎng)絡(luò)模型的穩(wěn)定性條件為：

Δt＜Δd_min/V_max；

其中，Δt為時間采樣間隔，Δd_min為空間采樣間隔的最小值，V_max為傳播速度的最大值。

綜上所述，本發(fā)明實施例具有以下優(yōu)點：

1.本方法在模擬過程中不依賴于傳統(tǒng)的彈性波方程，不受其相關(guān)假設(shè)條件限制，能夠適應(yīng)實際復(fù)雜介質(zhì)模型。

2.依據(jù)穩(wěn)定性條件和相速度頻散曲線優(yōu)選的時間與空間采樣間隔等模擬參數(shù)，保證了模擬的波場的正確性。

3.鑒于本發(fā)明提出的正演方法可以適用于矩形網(wǎng)格，具有一定的網(wǎng)格選擇靈活性，因而能夠滿足實際生產(chǎn)中對于不同方向的計算精度以及計算效率的要求。

本領(lǐng)域內(nèi)的技術(shù)人員應(yīng)明白，本發(fā)明的實施例可提供為方法、系統(tǒng)、或計算機程序產(chǎn)品。因此，本發(fā)明可采用完全硬件實施例、完全軟件實施例、或結(jié)合軟件和硬件方面的實施例的形式。而且，本發(fā)明可采用在一個或多個其中包含有計算機可用程序代碼的計算機可用存儲介質(zhì)(包括但不限于磁盤存儲器、CD-ROM、光學(xué)存儲器等)上實施的計算機程序產(chǎn)品的形式。

本發(fā)明是參照根據(jù)本發(fā)明實施例的方法、設(shè)備(系統(tǒng))、和計算機程序產(chǎn)品的流程圖和/或方框圖來描述的。應(yīng)理解可由計算機程序指令實現(xiàn)流程圖和/或方框圖中的每一流程和/或方框、以及流程圖和/或方框圖中的流程和/或方框的結(jié)合。可提供這些計算機程序指令到通用計算機、專用計算機、嵌入式處理機或其他可編程數(shù)據(jù)處理設(shè)備的處理器以產(chǎn)生一個機器，使得通過計算機或其他可編程數(shù)據(jù)處理設(shè)備的處理器執(zhí)行的指令產(chǎn)生用于實現(xiàn)在流程圖一個流程或多個流程和/或方框圖一個方框或多個方框中指定的功能的裝置。

這些計算機程序指令也可存儲在能引導(dǎo)計算機或其他可編程數(shù)據(jù)處理設(shè)備以特定方式工作的計算機可讀存儲器中，使得存儲在該計算機可讀存儲器中的指令產(chǎn)生包括指令裝置的制造品，該指令裝置實現(xiàn)在流程圖一個流程或多個流程和/或方框圖一個方框或多個方框中指定的功能。

這些計算機程序指令也可裝載到計算機或其他可編程數(shù)據(jù)處理設(shè)備上，使得在計算機或其他可編程設(shè)備上執(zhí)行一系列操作步驟以產(chǎn)生計算機實現(xiàn)的處理，從而在計算機或其他可編程設(shè)備上執(zhí)行的指令提供用于實現(xiàn)在流程圖一個流程或多個流程和/或方框圖一個方框或多個方框中指定的功能的步驟。

以上所述僅為本發(fā)明的優(yōu)選實施例而已，并不用于限制本發(fā)明，對于本領(lǐng)域的技術(shù)人員來說，本發(fā)明實施例可以有各種更改和變化。凡在本發(fā)明的精神和原則之內(nèi)，所作的任何修改、等同替換、改進等，均應(yīng)包含在本發(fā)明的保護范圍之內(nèi)。