技術領域
[0001] 本
發(fā)明屬于地球物理勘探研究領域，具體提供一種高溫介質
地震波傳播模擬方法。
背景技術
[0002] 地震波波場正演模擬越來越受到地球物理學者們的重視,它的地位也顯得越來越重要，正演模擬技術的日趨成熟也使得反演技術在
地震勘探中應用的準確性提高。而傳統(tǒng)的正演模型往往是常溫介質，并沒有考慮影響彈性介質物理屬性的
溫度因素，特別是高溫狀態(tài)下介質的彈性屬性變化。
[0003] 眾所周知，在不同的溫度下，其
巖石密度、
彈性模量、
波速會有較大的差異，尤其對于有裂縫存在的情況下。隨著溫度的增加，原始的縫洞開裂密度逐漸變大，其必然會引起儲層巖性物性參數(shù)的變化。所以在實際的深層—超深層地震勘探中，必須考慮到溫度的影響，才能更加準確的地反映真實地下構造，故高溫介質中地震波傳播模擬方法至關重要。
[0004] 相應地，本領域需要一種高溫介質地震波傳播模擬方法來解決上述問題。
發(fā)明內容
[0005] 為了解決
現(xiàn)有技術中的上述問題，即為了解決傳統(tǒng)的地震波正演模型為常溫介質導致的模擬結果不準確的問題，本發(fā)明提供了一種高溫介質地震波傳播模擬方法，所述模擬方法包括以下步驟：
[0006] 步驟S10、利用邊界元法將熱彈性
波動方程轉化為邊界元形式方程；
[0007] 步驟S20、利用邊界元形式方程求得地震波在均勻巖石物理模型區(qū)域的邊界的應
力值和位移值；
[0008] 步驟S30、將所述邊界的
應力值和位移值代入邊界元形式方程求得地震波在均勻巖石物理模型內區(qū)域的應力值和位移值；
[0009] 步驟S40、根據(jù)均勻巖石物理模型區(qū)域的邊界和內區(qū)域的全部應力值和位移值確定地震波的
波形記錄。
[0010] 在上述高溫介質地震波傳播模擬方法的優(yōu)選技術方案中，所述步驟S10具體包括以下步驟：
[0011] 步驟S101、將熱彈性波動方程寫成
頻率域形式并轉化為二維的情形如下：
[0012]
[0013] 公式(1)中，
[0014] θ＝T-T0,γ＝(3λ+2μ)αt，
[0015] 其中，Δ代表拉普拉斯算子，i為虛數(shù)單位，λ和μ代表
拉梅常數(shù)，u代表質點的位移，ux代表質點在x方向的位移，uz代表質點在z方向的位移，γ代表應力的熱力系數(shù)，θ代表介質的相對溫度，ρ代表介質的密度，ω代表
角頻率，Xx和Xz分別代表x方向和z方向的體力，i代表第i個離散單元，κ代表熱擴散系數(shù)，Q代表內熱源，k代表導熱系數(shù)，T0代表介質的初始溫度，T代表介質的實時溫度，αt代表線膨脹系數(shù)，η為中間變量。
[0016] 步驟S102、將上述熱彈性波動方程轉化為如下邊界元形式方程：
[0017]
[0018] 公式(2)中，u和p是公式(2)待求的位移值和應力值，i和j為第i個或者第j個離散單元，N為離散的總單元數(shù)，Σ*和V*是熱彈性波動方程的應力基本解和位移的基本解，Γ為離散單元的邊界。
[0019] 在上述高溫介質地震波傳播模擬方法的優(yōu)選技術方案中，所述步驟S20進一步包括：
[0020] 步驟S201、給出均勻巖石物理模型參數(shù)、觀測系統(tǒng)中各檢波點及各炮點的
位置坐標、邊界元法網格單元的各個控制點坐標及區(qū)域的劃分、
震源子波的信息以及均勻巖石物理模型區(qū)域邊界；
[0021] 步驟S202、將上述所有信息代入邊界元形式方程中，求得地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值。
[0022] 在上述高溫介質地震波傳播模擬方法的優(yōu)選技術方案中，在步驟S202之前，所述模擬方法還包括以下步驟：
[0023] 利用傅里葉變換將所述的時間域的熱彈性波動方程轉化為頻率域的熱彈性波動方程，以便在頻率域進行邊界元形式方程的計算。
[0024] 在上述高溫介質地震波傳播模擬方法的優(yōu)選技術方案中，所述步驟S202進一步包括：
[0025] 步驟S2021、利用形函數(shù)對所述均勻巖石物理模型的邊界進行離散，計算離散點坐標，同時離散外邊界吸收層數(shù)；
[0026] 步驟S2022、利用貝塞爾函數(shù)計算所述邊界元形式方程的積分核函數(shù)，得出第一邊界矩陣系數(shù)；
[0027] 步驟S2023、利用ASSEM函數(shù)對內邊界和外邊界處理，即將所述第一邊界矩陣系數(shù)進行組裝；
[0028] 步驟S2024、采用無限元方法進行人工截斷邊界處理得出第二邊界矩陣系數(shù)，并將所述第二邊界矩陣系數(shù)組合到所述第一邊界矩陣系數(shù)中；
[0029] 步驟S2025、利用
塊狀高斯消元法對代入步驟S2024中組合后的第一邊界矩陣系數(shù)的邊界元形式方式進行求解，得出地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值。
[0030] 在上述高溫介質地震波傳播模擬方法的優(yōu)選技術方案中，所述步驟S30進一步包括：
[0031] 步驟S301、將地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值代入邊界元形式方程；
[0032] 步驟S302、利用貝塞爾函數(shù)計算上述邊界元形式方程的積分核函數(shù)，得出第三邊界矩陣系數(shù)；
[0033] 步驟S303、利用塊狀高斯消元法對代入第三邊界矩陣系數(shù)的邊界元形式方式進行求解，得出地震波在均勻巖石物理模型的內區(qū)域的應力值和位移值。
[0034] 在上述高溫介質地震波傳播模擬方法的優(yōu)選技術方案中，所述步驟S40進一步包括：利用傅里葉變換將地震波在均勻巖石物理模型的邊界的應力值和位移值以及內區(qū)域的應力值和位移值從頻率域轉化為時間域；根據(jù)時間域下的地震波在均勻巖石物理模型的邊界的應力值和位移值以及內區(qū)域的應力值和位移值獲得地震波的波形記錄。
[0035] 在上述高溫介質地震波傳播模擬方法的優(yōu)選技術方案中，所述均勻巖石物理模型參數(shù)包括密度ρ、拉梅常數(shù)λ和μ、熱擴散系數(shù)κ、
熱膨脹系數(shù)αt、單位初始體積常應變
比熱cε。
[0036] 在上述高溫介質地震波傳播模擬方法的優(yōu)選技術方案中，所述震源子波為高斯子波或者雷克子波，所述震源子波信息包括震源子波的起始頻率、中心頻率和終止頻率。
[0037] 本領域技術人員能夠理解的是，在本發(fā)明的技術方案中，高溫介質地震波傳播模擬方法包括以下步驟：步驟S10、利用邊界元法將熱彈性波動方程轉化為邊界元形式方程；步驟S20、利用邊界元形式方程求得地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值；步驟S30、將所述邊界的應力值和位移值代入邊界元形式方程求得地震波在均勻巖石物理模型內區(qū)域的應力值和位移值；步驟S40、根據(jù)均勻巖石物理模型區(qū)域的邊界和內區(qū)域的全部應力值和位移值確定地震波的波形記錄。
[0038] 通過這樣的設置，有如下有益效果：1、由于邊界元法的理論
基礎是邊界積分方程理論，其特點是
降維、解析與離散相結合、能自動滿足無窮遠處邊界條件，因此該方法只需在區(qū)域的邊界上劃分單元，不必在無窮遠邊界上劃分單元，且計算
精度高，適于無限域問題。2、引入溫度參量的熱彈性波動方程，計及到溫度的影響，使結果更加貼切真實情況。本發(fā)明實現(xiàn)了基于熱彈性波動方程的高溫介質地震波傳播模擬方法，開展熱彈理論研究，可以準確考察高溫對地震波傳播的影響，對深層及超深層油氣勘探領域拓展提供理論基礎，對解釋地球深部觀測數(shù)據(jù)異?，F(xiàn)象提供合理的解釋。
附圖說明
[0039] 下面參照附圖來描述本發(fā)明的優(yōu)選實施方式，附圖中：
[0040] 圖1是本發(fā)明一種
實施例的高溫介質地震波傳播模擬方法的主要步驟
流程圖；
[0041] 圖2是本發(fā)明一種實施例的高溫介質地震波傳播模擬方法的具體步驟流程圖；
[0042] 圖3是本發(fā)明一種實施例的地震波
水平方向的波形圖；
[0043] 圖4是本發(fā)明一種實施例的地震波垂直方向的波形圖。
具體實施方式
[0044] 為使本發(fā)明實施例的目的、技術方案和優(yōu)點更加清楚，下面將結合本發(fā)明實施例中的附圖，對本發(fā)明實施例中的技術方案進行清楚、完整地描述，顯然，所描述的實施例是本發(fā)明一部分實施例，而不是全部的實施例?；诒景l(fā)明中的實施例，本領域普通技術人員在沒有做出創(chuàng)造性勞動前提下所獲得的所有其他實施例，都屬于本發(fā)明保護的范圍。
[0045] 下面參照附圖來描述本發(fā)明的優(yōu)選實施方式。本領域技術人員應當理解的是，這些實施方式僅僅用于解釋本發(fā)明的技術原理，并非旨在限制本發(fā)明的保護范圍。
[0046] 地震波波場正演模擬方法的種類很多，一般將其分為波動方程法、積分方程法、幾何射線法三大類。我們基于耦合熱彈性理論建立熱彈性波動方程，利用邊界元數(shù)值方法進行高溫下地震波傳播正演模擬。熱彈性耦合理論由
熱力學基本定律、介質的本構理論和Helmholtz自由能表達式導出的熱彈性材料的熱傳導方程中，除了待定的溫度場函數(shù)外，還含有應變率。它表明物體上的溫度場不僅取決于熱源，以及各有關的熱力學物性系數(shù)和換熱邊界條件，而且還受到彈性
變形應變率的影響，或者說，彈性變形的體積應變率將在一定的程度上改變物體上的熱量的傳遞，因此熱傳導方程和熱彈性運動方程必須耦合求解。
[0047] 邊界元法是對邊界積分方程離散求解的現(xiàn)代數(shù)值分析方法，屬于近似解法之一。邊界積分方程從定解問題的控制方程轉換而得，其最大的特點就是降低了求解問題的維數(shù)，將二維問題化為其邊界線上的一維問題，它只以邊界變量為基本變量，域內位置可以在需要時根據(jù)邊界變量求出。基于熱彈性理論及邊界元法，我們發(fā)展一種熱介質中彈性波邊界元法正演模擬。
[0048] 如圖1所示，圖1是本發(fā)明一種實施例的高溫介質地震波傳播模擬方法的主要步驟流程圖。參照圖1，本發(fā)明的高溫介質地震波傳播模擬方法主要包括以下步驟：
[0049] 步驟S10、利用邊界元法將熱彈性波動方程轉化為邊界元形式方程；
[0050] 優(yōu)選地，步驟S10進一步包括：步驟S101、將熱彈性波動方程寫成頻率域形式并轉化為二維的情形如下：
[0051]
[0052] 公式(1)中，
[0053] θ＝T-T0,γ＝(3λ+2μ)αt,
[0054] 其中，Δ代表拉普拉斯算子，i為虛數(shù)單位，λ和μ代表拉梅常數(shù)，u代表質點的位移，ux代表質點在x方向的位移，uz代表質點在z方向的位移，γ代表應力的熱力系數(shù)，θ代表介質的相對溫度，ρ代表介質的密度，ω代表角頻率，Xx和Xz分別代表x方向和z方向的體力，i代表第i個離散單元，κ代表熱擴散系數(shù)，Q代表內熱源，k代表導熱系數(shù)，T0代表介質的初始溫度，T代表介質的實時溫度，αt代表線膨脹系數(shù)，η為中間變量。
[0055] 如下圖所示，在本發(fā)明的實施例中常數(shù)和一些參數(shù)的取值如下：
[0056] ρ(kg·m-3) λ(N·m-2) μ(N·m-2) κ(m2/s) αt(1/K) cε(J/m3/K)2.65 4×109 6×109 0.089 0.33×10-5 1.17×102
[0057] 步驟S102、將上述熱彈性波動方程轉化為如下邊界元形式方程：
[0058]
[0059] 公式(2)中，Σ*和V*是熱彈性波動方程的應力基本解和位移的基本解，u和p是公式(2)待求的位移值和應力值，i和j為第i個或者第j個離散單元，N為離散的總單元數(shù)，Γ為離散單元的邊界。
[0060] 一方面，通過引入溫度參量的熱彈性波動方程，計及到溫度的影響，使模擬結果更加貼切實際情況。另一方面，由于邊界元法的理論基礎是邊界積分方程理論，其特點是降維、解析與離散相結合、能自動滿足無窮遠處邊界條件，因此該方法只需在區(qū)域的邊界上劃分單元，不必在無窮遠邊界上劃分單元，且計算精度高，適于無限域問題。
[0061] 步驟S20、利用邊界元形式方程求得地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值。
[0062] 在本發(fā)明的一種實施例中，步驟S20進一步包括：步驟S201、給出均勻巖石物理模型參數(shù)、觀測系統(tǒng)中各檢波點及各炮點的位置坐標、邊界元法網格單元的各個控制點坐標及區(qū)域的劃分、震源子波的信息以及均勻巖石物理模型區(qū)域邊界。
[0063] 其中，均勻巖石物理模型參數(shù)具體包括：密度ρ、拉梅常數(shù)λ和μ、熱擴散系數(shù)κ、
熱膨脹系數(shù)αt、單位初始體積常應變比熱cε共六個參數(shù)。這六個參數(shù)的信息需要從靶區(qū)的原位巖石獲取，對靶區(qū)的原位巖石進行高溫原位熱巖石物理實驗測量，使求得的參數(shù)更加準確、真實可靠。其中，震源子波的信息包括起止頻率、中心頻率、終止頻率等信息。
[0064] 步驟S202、將上述所有信息代入邊界元形式方程中，求得地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值。
[0065] 需要說明的是，由于在時間域進行邊界元形式方程的計算比較復雜，因此為了方便計算，在步驟S202之前，所述方法還包括：利用傅里葉變換將時間域的熱彈性波動方程轉化為頻率域的熱彈性波動方程，即公式(1)所示，以便在頻率域進行邊界元形式方程的計算。
[0066] 在本發(fā)明的一種實施例中，步驟S202進一步包括：
[0067] 步驟S2021、利用形函數(shù)對均勻巖石物理模型的邊界進行離散，計算離散點坐標，同時離散外邊界吸收層數(shù)；
[0068] 如步驟S201已經給出邊界元法網格單元的各個控制點坐標及區(qū)域的劃分，利用形函數(shù)及區(qū)域邊界的控制點，對區(qū)域邊界中進行插值，從而計算離散點坐標及離散外邊界吸收層數(shù)；
[0069] 步驟S2022、利用貝塞爾函數(shù)計算所述邊界元形式方程的積分核函數(shù)，得出第一邊界矩陣系數(shù)；
[0070] 步驟S2023、利用ASSEM函數(shù)對內邊界和外邊界處理，即將所述第一邊界矩陣系數(shù)進行組裝；
[0071] 步驟S2024、采用無限元方法進行人工截斷邊界處理得出第二邊界矩陣系數(shù)，并將所述第二邊界矩陣系數(shù)組合到所述第一邊界矩陣系數(shù)中；即在得到第一邊界矩陣系數(shù)的基礎上，再增加擴充矩陣邊界系數(shù)。
[0072] 步驟S2025、利用塊狀高斯消元法對代入步驟S2024中組合后的第一邊界矩陣系數(shù)的邊界元形式方式進行求解，得出地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值。
[0073] 通過求出地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值，來進行地震波在均勻巖石物理模型內區(qū)域的應力值和位移值的求解。
[0074] 步驟S30、將邊界的應力值和位移值代入邊界元形式方程求得地震波在均勻巖石物理模型內區(qū)域的應力值和位移值；
[0075] 在本發(fā)明的一種實施例中，步驟S30進一步包括：
[0076] 步驟S301、將地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值代入邊界元形式方程；
[0077] 步驟S302、利用貝塞爾函數(shù)計算上述邊界元形式方程的積分核函數(shù)，得出第三邊界矩陣系數(shù)；
[0078] 步驟S303、利用塊狀高斯消元法對代入第三邊界矩陣系數(shù)的邊界元形式方式進行求解，得出地震波在均勻巖石物理模型的內區(qū)域的應力值和位移值。
[0079] 可以理解的是，要得出多組地震波在均勻巖石物理模型的邊界以及內區(qū)域的應力值和位移值，需要輸入震源子波的某個頻率類型的不同頻率，得出地震波在震源子波的不同頻率下的在均勻巖石物理模型的邊界以及內區(qū)域的應力值和位移值。如改變給定的震源子波的起始頻率：30赫茲、40赫茲、50赫茲等等，計算出地震波在相對應的赫茲下的在均勻巖石物理模型的邊界以及內區(qū)域的應力值和位移值。
[0080] 步驟S40、根據(jù)均勻巖石物理模型區(qū)域的邊界和內區(qū)域的全部應力值和位移值確定地震波的波形記錄。
[0081] 可以理解的是，由于計算出來的應力值和位移值是在頻率下，因此步驟S40進一步包括：利用傅里葉變換將地震波在均勻巖石物理模型的邊界的應力值和位移值以及內區(qū)域的應力值和位移值從頻率域轉化為時間域；根據(jù)時間域下的地震波在均勻巖石物理模型的邊界的應力值和位移值以及內區(qū)域的應力值和位移值獲得地震波的波形記錄。
[0082] 需要說明的是，上述實施例是在模擬的巖石物理模型為均勻介質下進行的，在巖石物理模型為多層不均勻時，需要引入ANGLE函數(shù)對每層地程的傾斜角度進行計算，然后再進行步驟S20、S30和S40。
[0083] 在一種可能的實施例中，如圖2所示，圖2是本發(fā)明一種實施例的高溫介質地震波傳播模擬方法的具體步驟流程圖。參照圖2，高溫介質地震波傳播模擬方法包括：
[0084] 1)給定輸入的均勻巖石物理模型，其模型包括：密度ρ、拉梅常數(shù)λ和μ、熱擴散系數(shù)κ、熱膨脹系數(shù)αt、單位初始體積常應變比熱cε共六個參數(shù)的信息；2)給出觀測系統(tǒng)中各檢波點及各炮點的位置坐標；3)給出邊界元法網格單元的各個控制點坐標及區(qū)域的劃分；4)給出震源子波的信息，包括起止頻率、中心頻率、終止頻率等信息；5)根據(jù)前面4個步驟得到的所有輸入信息，輸入子模塊EDITP中，同時子模塊EDITP包含模擬的均勻巖石物理模型區(qū)域邊界以及模擬震源子波的類型(程序中可
選定高斯子波或者雷克子波)；我們選用形函數(shù)對均勻巖石物理模型的邊界進行離散處理，從而計算邊界上的離散點坐標，同時離散外界吸收
邊界層數(shù)；利用傅里葉變換(FFT)對時間域的熱彈性波動方程轉化為頻率域的熱彈性波動方程，從而在頻率域進行邊界元形式方程的離散；調用ASSEM函數(shù)對第一邊界矩陣系數(shù)進行組裝，其中第一邊界矩陣系數(shù)為通過調用貝塞爾函數(shù)計算邊界元形式方程的積分核函數(shù)得出；本程序中采用無限元方法進行人工截斷邊界處理對第一邊界矩陣系數(shù)進行重新組裝；然后通過塊狀高斯消元求解，得出地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值；最后將地震波在均勻巖石物理模型區(qū)域的邊界的應力值和位移值代入邊界元形式方程中，通過貝塞爾函數(shù)和塊狀高斯消元求解檢波點波場(即地震波在均勻巖石物理模型的區(qū)域內的應力值和位移值)；6)將頻率域的結果轉換為時間域，得出全部
節(jié)點的位移和應力(即地震波在均勻巖石物理模型的區(qū)域內和邊界的應力值和位移值)，為后續(xù)的反演做準備，從而進行后續(xù)的分析。
[0085] 如圖3和圖4所示，圖3是本發(fā)明一種實施例的地震波水平方向的波形圖；圖4是本發(fā)明一種實施例的地震波垂直方向的波形圖。參照圖3和圖4，從特定的輸入參數(shù)中，我們得到地震波水平和垂直方向的記錄，圖中我們可以看出，高溫介質下的地震波傳播記錄中一共存在三種波形，即第一縱波、第二縱波(熱波)及橫波。在不考慮溫度效應中，即傳統(tǒng)的均勻各項同性介質中的僅僅存在兩種記錄即橫波和縱波記錄，我們利用彈性理論模擬地震波在熱介質中的傳播有新的波形出現(xiàn)，為地震勘探提供新思路。
[0086] 通過上述描述可以看出，本發(fā)明的高溫介質地震波傳播模擬方法具有以下優(yōu)點：1、輸入的模型巖石物理參數(shù)準確，其中密度ρ、拉梅常數(shù)λ和μ、熱擴散系數(shù)κ、熱膨脹系數(shù)αt、單位初始體積常應變比熱cε共六個參數(shù)的信息共六個參數(shù)的信息需要靶區(qū)的原位巖石，從而進行熱巖石物理實驗測量，求得的參數(shù)更加真實可靠。2、利用邊界元法，邊界元法的理論基礎是邊界積分方程理論，其特點是降維、解析與離散相結合、能自動滿足無窮遠處邊界條件。該方法只需在區(qū)域的邊界上劃分單元，不必在無窮遠邊界上劃分單元，計算精度高，適于無限域問題。3、引入溫度參量的熱彈性波動方程，計及到溫度的影響，使結果更加貼切真實情況。本發(fā)明實現(xiàn)了基于熱彈性波動方程的高溫介質地震波傳播模擬方法，開展熱彈理論研究，可以準確考察高溫對地震波傳播的影響，對深層及超深層油氣勘探領域拓展提供理論基礎，對解釋地球深部觀測數(shù)據(jù)異?，F(xiàn)象提供合理的解釋。
[0087] 至此，已經結合附圖所示的優(yōu)選實施方式描述了本發(fā)明的技術方案，但是，本領域技術人員容易理解的是，本發(fā)明的保護范圍顯然不局限于這些具體實施方式。在不偏離本發(fā)明的原理的前提下，本領域技術人員可以對相關技術特征作出等同的更改或替換，這些更改或替換之后的技術方案都將落入本發(fā)明的保護范圍之內。
