基于彈簧網絡模型的正演模擬方法及裝置與流程

本發明涉及石油地震勘探技術領域，特別涉及一種基于彈簧網絡模型的正演模擬方法及裝置。

背景技術：

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

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

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

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

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

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

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

技術實現要素：

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

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

根據所述空間采樣間隔對待模擬的地質模型進行離散化，獲得數字化地質模型；所述數字化地質模型為矩形網格形式或長方體網格形式；

根據所述空間采樣間隔確定所述數字化地質模型的各個方向的彈性參數；

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

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

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

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

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

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

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

附圖說明

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

圖1為本發明實施例提供的一種基于彈簧網絡模型的正演模擬方法的處理流程圖；

圖2為本發明實施例提供的一種基于彈簧網絡模型的正演模擬方法具體處理流程圖；

圖3為本發明實施例提供的二維(D2Q8)彈簧網絡模型示意圖；

圖4為本發明實施例提供的三維(D3Q18)彈簧網絡模型示意圖；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9為本發明實施例提供的一種基于彈簧網絡模型的正演模擬裝置的結構示意圖；

圖10為本發明實施例提供的一種基于彈簧網絡模型的正演模擬裝置的一種具體結構示意圖。

具體實施方式

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

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

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

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

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

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

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

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

當然，本發明方法對正方形網格或者正方體網格同樣適用。

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

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

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

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

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

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

且有，

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

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

且有，

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

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

其中，Δt為模擬采用的時間采樣間隔，F_i為中心節點i所受到的線彈簧的總作用力，m_i為中心節點i的質量，χ為粘滯項系數，用于衰減節點的震動能量，這里取χ＝0。

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

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

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

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

中間變量計算公式如下：

上式中

且有

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

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

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

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

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

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

在實際數值模擬過程中，彈簧網絡模型(LSM)的穩定性條件為：

Δt＜Δd_min/V_max (14)

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

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

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

下面通過比較來說明本方法的優點。

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

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

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

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

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

圖9是本發明實施例的基于彈簧網絡模型的正演模擬裝置的一種結構框圖，如圖9所示，包括：

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

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

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

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

下面對該結構進行說明。

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

將待模擬的地質模型離散成由離散網格節點與線彈簧構成的彈性元，相鄰離散網格節點之間通過線彈簧連接；相鄰離散網格節點之間的距離為所述空間采樣間隔；

矩形網格形式的數字化地質模型由彈性面元構成；

長方體網格形式的數字化地質模型由彈性體元構成。

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

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

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

對于矩形網格形式的數字化地質模型，E為彈性面元的面積，N＝8；

對于長方體網格形式的數字化地質模型，E為彈性體元的體積，N＝18；

對于矩形網格形式的數字化地質模型，各個方向的彈性參數按如下公式確定：

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

對于長方體網格形式的數字化地質模型，各個方向的彈性參數按如下公式確定：

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

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

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

其中，x_i為位移；v_i為速度；a_i為加速度；Δt為時間采樣間隔，F_i為中心節點i所受到的線彈簧的總作用力，m_i為中心節點i的質量，χ為粘滯項系數，χ＝0。

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

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

按如下公式確定離散格式的彈簧網絡模型的相速度頻散曲線按照如下公式確定：

中間變量計算公式如下：

上式中

且有

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

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

具體實施時，所述彈簧網絡模型的穩定性條件為：

Δt＜Δd_min/V_max；

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

綜上所述，本發明實施例具有以下優點：

1.本方法在模擬過程中不依賴于傳統的彈性波方程，不受其相關假設條件限制，能夠適應實際復雜介質模型。

2.依據穩定性條件和相速度頻散曲線優選的時間與空間采樣間隔等模擬參數，保證了模擬的波場的正確性。

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

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

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

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

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

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