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

背景技術(shù)：

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

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

彈簧網(wǎng)絡(luò)模型(Lattice Spring Model，簡(jiǎn)稱LSM)是一種從微觀角度研究介質(zhì)彈塑性的方法，其相關(guān)原理可追溯到20世紀(jì)初Born等學(xué)者提出的晶體動(dòng)力學(xué)理論。早期的Born彈簧模型，相鄰質(zhì)點(diǎn)之間只考慮中心力(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)聚能和晶格動(dòng)力學(xué)等打下了堅(jiān)實(shí)基礎(chǔ)。之后，Kittel、Lax和Keating等人對(duì)Born彈簧模型作了進(jìn)一步地發(fā)展與完善，保證了形變后晶體的應(yīng)變能量計(jì)算公式具有旋轉(zhuǎn)不變性。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

附圖說明

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖8a)為本發(fā)明實(shí)施例提供的Marmousi模型的FDM模擬得到的波場(chǎng)(v_x)快照?qǐng)D；

圖8b)為本發(fā)明實(shí)施例提供的Marmousi模型的FDM模擬得到的波場(chǎng)(v_z)快照?qǐng)D；

圖8c)為本發(fā)明實(shí)施例提供的Marmousi模型的LSM模擬得到的波場(chǎng)(v_x)快照?qǐng)D；

圖8d)為本發(fā)明實(shí)施例提供的Marmousi模型的LSM模擬得到的波場(chǎng)(v_z)快照?qǐng)D；

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

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

具體實(shí)施方式

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

且有，

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

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

且有，

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

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

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

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

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

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

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

中間變量計(jì)算公式如下：

上式中

且有

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

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

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

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

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

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

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

Δt＜Δd_min/V_max (14)

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

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

最后，輸出并保存計(jì)算獲得的各個(gè)時(shí)刻的波場(chǎng)值。

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

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

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

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

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

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

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

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

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

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

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

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

具體實(shí)施時(shí)，離散模塊902具體用于：

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

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

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

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

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

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

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

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

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

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

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

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

具體實(shí)施時(shí)，所述波場(chǎng)值確定模塊904具體用于：

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

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

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

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

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

中間變量計(jì)算公式如下：

上式中

且有

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

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

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

Δt＜Δd_min/V_max；

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

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

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

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

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

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

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

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

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

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