地震波模擬方法及系統(tǒng)
技術(shù)領(lǐng)域
[0001] 本
發(fā)明涉及油汽地球物理技術(shù)領(lǐng)域，更具體地，涉及一種
地震波模擬方法及系統(tǒng)。
背景技術(shù)
[0002] 偽譜法與有限差分法、有限元法并列為地震
波數(shù)值模擬的三種主要方法。它通過數(shù)學(xué)變換的方法把物理變量對空間的微分轉(zhuǎn)化為轉(zhuǎn)換域中的代數(shù)運(yùn)算，然后再把結(jié)果通過逆變換轉(zhuǎn)換到物理空間，從而求得對應(yīng)量的空間微分。相對于有限差分法和有限元法，偽譜法具有
精度高的特點，運(yùn)算速度比較快的特點。因此，偽譜法在
地震勘探中有著比較廣泛的應(yīng)用。
[0003] 目前地震領(lǐng)域中的偽譜法主要指傅里葉偽譜法。傅里葉偽譜法中的基函數(shù)為三
角函數(shù)，具有周期特性，因此，傅里葉偽譜法求得的空間導(dǎo)數(shù)值也是周期的，這種偽譜法一般用于周期邊界條件的情況，而不適合有固定條件的情況。地球介質(zhì)是有邊界的，如地表、固定底面，因此邊界條件往往在波場模擬中需要考慮。用傳統(tǒng)的傅里葉偽譜法不能很好處理這類含邊界的問題。因此，有必要開發(fā)一種地震波模擬方法及系統(tǒng)。
[0004] 公開于本發(fā)明背景技術(shù)部分的信息僅僅旨在加深對本發(fā)明的一般背景技術(shù)的理解，而不應(yīng)當(dāng)被視為承認(rèn)或以任何形式暗示該信息構(gòu)成已為本領(lǐng)域技術(shù)人員所公知的
現(xiàn)有技術(shù)。
發(fā)明內(nèi)容
[0005] 本發(fā)明提出了一種地震波模擬方法及系統(tǒng)，其能夠通過傅里葉偽譜法計算
水平導(dǎo)數(shù)，用切比
雪夫偽譜法計算垂直方向?qū)?shù)，在解決自由面模擬的同時，保證了模擬的高精度。
[0006] 根據(jù)本發(fā)明的一方面，提出了一種地震波模擬方法。所述方法可以包括：步驟1：建立含自由表面的
地層介質(zhì)中的地震波
波動方程，以g表示其中任意一個變量；步驟2：對三維空間進(jìn)行網(wǎng)格剖分，其中x方向采用等距網(wǎng)格，z方向采用高斯-切比雪夫-洛巴托點獲得空間離散點；步驟3：針對某一時刻點的內(nèi)部網(wǎng)格點，x方向采用傅里葉方程求解下一時刻點的內(nèi)部網(wǎng)格點的變量g，z方向采用采用切比雪夫方程求解下一時刻點的內(nèi)部網(wǎng)格點的變量g；步驟4：施加物理邊界條件和數(shù)值邊界條件，針對某一時刻點進(jìn)行邊界方程求解，獲得下一時刻點的邊界上的變量g；步驟5：隨時間步進(jìn)，重復(fù)步驟3-4，直至獲得每個時刻點的內(nèi)部網(wǎng)格點和邊界上的變量g。
[0007] 優(yōu)選地，x方向采用傅里葉方程求解變量g的一階微分包括：
[0008]
[0009] 其中，k代表波數(shù)，Δk為波數(shù)間隔，Δk＝2π/(NΔx)，i是虛數(shù)單位，G(lΔk)表示函數(shù)g(nΔx)的傅里葉變換，g(nΔx)為g在x方向的離散表達(dá)式。
[0010] 優(yōu)選地，g(nΔx)的表達(dá)式為：
[0011]
[0012] 式中，n,l為
節(jié)點索引，N為傅里葉多項式的最高階數(shù)。
[0013] 優(yōu)選地，z方向的高斯-切比雪夫-洛巴托離散點的表達(dá)式為：
[0014]
[0015] 其中，計算區(qū)間z被限制在[-1,1]之間，端點z0＝1，zM＝-1，M是運(yùn)算中所取切比雪夫多項式的階數(shù)。
[0016] 優(yōu)選地，z方向采用切比雪夫方程求解變量g的一階微分包括：
[0017]
[0018] 其中， 的離散表達(dá)式為
[0019]
[0020] 其中，Dc為(M+1)×(M+1)的矩陣，
[0021] 根據(jù)本發(fā)明的另一方面，提出了一種地震波模擬系統(tǒng)，其上存儲有
計算機(jī)程序，其特征在于，所述程序被處理器執(zhí)行時實現(xiàn)以下步驟：步驟1：建立含自由表面的地層介質(zhì)中的地震波波動方程，以g表示其中任意一個變量；步驟2：對三維空間進(jìn)行網(wǎng)格剖分，其中x方向采用等距網(wǎng)格，z方向采用高斯-切比雪夫-洛巴托點獲得空間離散點；步驟3：針對某一時刻點的內(nèi)部網(wǎng)格點，x方向采用傅里葉方程求解下一時刻點的內(nèi)部網(wǎng)格點的變量g，z方向采用采用切比雪夫方程求解下一時刻點的內(nèi)部網(wǎng)格點的變量g；步驟4：施加物理邊界條件和數(shù)值邊界條件，針對某一時刻點進(jìn)行邊界方程求解，獲得下一時刻點的邊界上的變量g；步驟5：隨時間步進(jìn)，重復(fù)步驟3-4，直至獲得每個時刻點的內(nèi)部網(wǎng)格點和邊界上的變量g。
[0022] 優(yōu)選地，x方向采用傅里葉方程求解變量g的一階微分包括：
[0023]
[0024] 其中，k代表波數(shù)，Δk為波數(shù)間隔，Δk＝2π/(NΔx)，i是虛數(shù)單位，G(lΔk)表示函數(shù)g(nΔx)的傅里葉變換，g(nΔx)為g在x方向的離散表達(dá)式。
[0025] 優(yōu)選地，g(nΔx)的表達(dá)式為：
[0026]
[0027] 式中，n,l為節(jié)點索引，N為傅里葉多項式的最高階數(shù)。
[0028] 優(yōu)選地，z方向的高斯-切比雪夫-洛巴托離散點的表達(dá)式為：
[0029]
[0030] 其中，計算區(qū)間z被限制在[-1,1]之間，端點z0＝1，zM＝-1，M是運(yùn)算中所取切比雪夫多項式的階數(shù)。
[0031] 優(yōu)選地，z方向采用切比雪夫方程求解變量g的一階微分包括：
[0032]
[0033] 其中， 的離散表達(dá)式為
[0034]
[0035] 其中，Dc為(M+1)×(M+1)的矩陣，
[0036] 本發(fā)明具有其它的特性和優(yōu)點，這些特性和優(yōu)點從并入本文中的
附圖和隨后的具體實施方式中將是顯而易見的，或者將在并入本文中的附圖和隨后的具體實施方式中進(jìn)行詳細(xì)陳述，這些附圖和具體實施方式共同用于解釋本發(fā)明的特定原理。
附圖說明
[0037] 通過結(jié)合附圖對本發(fā)明示例性
實施例進(jìn)行更詳細(xì)的描述，本發(fā)明的上述以及其它目的、特征和優(yōu)勢將變得更加明顯，其中，在本發(fā)明示例性實施例中，相同的參考標(biāo)號通常代表相同部件。
[0038] 圖1示出了根據(jù)本發(fā)明的地震波模擬方法的步驟的
流程圖。
[0039] 圖2示出了根據(jù)本發(fā)明的一個實施例的網(wǎng)格剖分的示意圖。
[0040] 圖3示出了根據(jù)本發(fā)明的一個實施例的二維x-z空間模型的示意圖。
[0041] 圖4示出了根據(jù)本發(fā)明的一個實施例的
震源子波在0.35秒時刻的波場快照的示意圖。
[0042] 圖5a、5b、5c分別示出了根據(jù)本發(fā)明的一個實施例的模擬得到的
檢波器R1、R2、R3的位移記錄和理論結(jié)果的對比圖。
具體實施方式
[0043] 下面將參照附圖更詳細(xì)地描述本發(fā)明。雖然附圖中顯示了本發(fā)明的優(yōu)選實施例，然而應(yīng)該理解，可以以各種形式實現(xiàn)本發(fā)明而不應(yīng)被這里闡述的實施例所限制。相反，提供這些實施例是為了使本發(fā)明更加透徹和完整，并且能夠?qū)⒈景l(fā)明的范圍完整地傳達(dá)給本領(lǐng)域的技術(shù)人員。
[0044] 圖1示出了根據(jù)本發(fā)明的地震波模擬方法的步驟的流程圖。
[0045] 在該實施例中，根據(jù)本發(fā)明的地震波模擬方法可以包括：步驟1：建立含自由表面的地層介質(zhì)中的地震波波動方程，以g表示其中任意一個變量；步驟2：對三維空間進(jìn)行網(wǎng)格剖分，其中x方向采用等距網(wǎng)格，z方向采用高斯-切比雪夫-洛巴托點獲得空間離散點；步驟3：針對某一時刻點的內(nèi)部網(wǎng)格點，x方向采用傅里葉方程求解下一時刻點的內(nèi)部網(wǎng)格點的變量g，z方向采用采用切比雪夫方程求解下一時刻點的內(nèi)部網(wǎng)格點的變量g；步驟4：施加物理邊界條件和數(shù)值邊界條件，針對某一時刻點進(jìn)行邊界方程求解，獲得下一時刻點的邊界上的變量g；步驟5：隨時間步進(jìn)，重復(fù)步驟3-4，直至獲得每個時刻點的內(nèi)部網(wǎng)格點和邊界上的變量g。
[0046] 在一個示例中，x方向采用傅里葉方程求解變量g的一階微分包括：
[0047]
[0048] 其中，k代表波數(shù)，Δk為波數(shù)間隔，Δk＝2π/(NΔx)，i是虛數(shù)單位，G(lΔk)表示函數(shù)g(nΔx)的傅里葉變換，g(nΔx)為g在x方向的離散表達(dá)式。
[0049] 在一個示例中，g(nΔx)的表達(dá)式為：
[0050]
[0051] 式中，n,l為節(jié)點索引，N為傅里葉多項式的最高階數(shù)。
[0052] 在一個示例中，z方向的高斯-切比雪夫-洛巴托離散點的表達(dá)式為：
[0053]
[0054] 其中，計算區(qū)間z被限制在[-1,1]之間，端點z0＝1，zM＝-1，M是運(yùn)算中所取切比雪夫多項式的階數(shù)。
[0055] 在一個示例中，z方向采用切比雪夫方程求解變量g的一階微分包括：
[0056]
[0057] 其中， 的離散表達(dá)式為
[0058]
[0059] 其中，Dc為(M+1)×(M+1)的矩陣，
[0060] 具體地，以一階SH波方程為例，根據(jù)本發(fā)明的地震波模擬方法可以包括：
[0061] 對三維空間進(jìn)行網(wǎng)格剖分，其中x方向采用等距網(wǎng)格，z方向采用高斯-切比雪夫-洛巴托點離散。根據(jù)彈性
力學(xué)原理，針對二維x-z空間建立一階SH波方程為公式(6)：
[0062]
[0063] 其中，V表示介質(zhì)沿垂直x-z平面的y方向的速度，τxy、τyz是
應(yīng)力分量，fy是單位
質(zhì)量上沿y方向的外力，ρ代表介質(zhì)
密度，μ是介質(zhì)的
剪切模量。
[0064] 針對某一時刻點，對二維x-z空間的內(nèi)部網(wǎng)格點進(jìn)行一階SH波方程求解，獲得下一時刻點的內(nèi)部網(wǎng)格點的速度、應(yīng)力，其中，x方向采用傅里葉方程求解，z方向采用切比雪夫方程求解。
[0065] 對于x方向，采用等距網(wǎng)格，得到對應(yīng)的函數(shù)離散值τxy、V，然后將離散后的函數(shù)變換至波數(shù)域分別為為公式(7)、(8)：
[0066]
[0067]
[0068] 其中， 表示函數(shù)τxy、V的傅里葉變換。根據(jù)函數(shù)微分的傅里葉變換和函數(shù)的傅里葉變換的關(guān)系，將 和ilΔk相乘，并將結(jié)果逆傅里葉變換回空間域，得到函數(shù)的1階空間微分結(jié)果，即通過公式(9)求解x方向的應(yīng)力分量的一階微分：
[0069]
[0070] 其中，τxy是應(yīng)力分量，k代表波數(shù)，波數(shù)間隔Δk＝2π/(NΔx)，i是虛數(shù)單位，通過公式(10)求解速度的一階微分：
[0071]
[0072] 其中，V表示介質(zhì)沿垂直x-z平面的y方向的速度，k代表波數(shù)，波數(shù)間隔Δk＝2π/(NΔx)。
[0073] 對于z方向，z方向采用高斯-切比雪夫-洛巴托點離散，其表達(dá)形式為公式(3)，其中計算區(qū)間z被限制在[-1,1]之間，端點z0＝1，zM＝-1，[a,b]為z方向的物理區(qū)間(a、b可以為有界的任意實數(shù))，通過坐標(biāo)轉(zhuǎn)換 映射至計算坐標(biāo)z’，其中z’屬于[-1,1]，然后通過公式(11)求解z方向的應(yīng)力分量的一階微分：
[0074]
[0075] 其中 的離散表達(dá)式為
[0076]
[0077] 其中，公式(11)中 [a,b]為z方向的物
理區(qū)間。
[0078] 通過公式(13)求解速度在z方向的一階微分：
[0079]
[0080] 其中 的離散表達(dá)式為
[0081]
[0082] 其 中 ，公式(14)中 [a,b]為z方向的物
理區(qū)間。
[0083] 求得 和 后，公式(6)右邊項的值就可以得到了(fy，ρ和μ是已知的)，于是公式(6)由偏微分方程組變?yōu)榱顺Ｎ⒎址匠探M，再利用數(shù)值積分技術(shù)(如歐拉法，4階龍格庫塔法等)求得下一時刻的V、τxy和τyz。
[0084] 施加物理邊界條件和數(shù)值邊界條件，針對某一時刻點，對二維x-z空間的邊界進(jìn)行邊界方程(15)求解：
[0085] τyz＝0,當(dāng)z＝0???(15)
[0086] 獲得下一時刻點的邊界上的位移和應(yīng)力。
[0087] 隨時間步進(jìn)，獲得每個時刻點的內(nèi)部網(wǎng)格點的速度、應(yīng)力和邊界上的位移、應(yīng)力。根據(jù)每個時刻點的內(nèi)部網(wǎng)格點的速度、應(yīng)力和邊界上的位移、應(yīng)力，獲得SH波的波場記錄，得到空間各場量移隨時間的變化關(guān)系。
[0088] 本方法通過傅里葉偽譜法計算水平導(dǎo)數(shù)，用切比雪夫偽譜法計算垂直方向?qū)?shù)，在解決自由面模擬的同時，保證了模擬的高精度。
[0089] 應(yīng)用示例
[0090] 為便于理解本發(fā)明實施例的方案及其效果，以下給出一個具體應(yīng)用示例。本領(lǐng)域技術(shù)人員應(yīng)理解，該示例僅為了便于理解本發(fā)明，其任何具體細(xì)節(jié)并非意在以任何方式限制本發(fā)明。
[0091] 圖2示出了根據(jù)本發(fā)明的一個實施例的網(wǎng)格剖分的示意圖。
[0092] 圖3示出了根據(jù)本發(fā)明的一個實施例的二維x-z空間模型的示意圖。
[0093] 圖4示出了根據(jù)本發(fā)明的一個實施例的震源子波在0.35秒時刻的波場快照的示意圖。
[0094] 根據(jù)本發(fā)明的地震波模擬方法包括：
[0095] 對三維空間進(jìn)行網(wǎng)格剖分，其中x方向采用等距網(wǎng)格，z方向采用高斯-切比雪夫-洛巴托點離散，如圖2所示，對如圖2所示的二維x-z空間進(jìn)行SH波數(shù)值模擬，安置了3個虛擬檢波器R1、R2、R3，來記錄空間點的位移，如圖3所示，模擬中震源采用50Hz的Ricker子波，如圖4所示。
[0096] 根據(jù)彈性力學(xué)原理，針對二維x-z空間建立一階SH波方程為公式(6)。
[0097] 針對某一時刻點，對二維x-z空間的內(nèi)部網(wǎng)格點進(jìn)行一階SH波方程求解，獲得下一時刻點的內(nèi)部網(wǎng)格點的速度、應(yīng)力，其中，x方向采用傅里葉方程求解，z方向采用切比雪夫方程求解。
[0098] 對于x方向，采用等距網(wǎng)格，得到對應(yīng)的函數(shù)離散值τxy、V，然后將離散后的函數(shù)變換至波數(shù)域為公式(7)、(8)。根據(jù)函數(shù)微分的傅里葉變換和函數(shù)的傅里葉變換的關(guān)系，將和ilΔk相乘，并將結(jié)果逆傅里葉變換回空間域，得到函數(shù)的1階空間微分結(jié)果，即通過公式(9)求解x方向的應(yīng)力分量的一階微分，通過公式(10)求解速度的一階微分。
[0099] 對于z方向，z方向采用高斯-切比雪夫-洛巴托點離散，其中計算區(qū)間z被限制在[-1,1]之間，端點z0＝1，zM＝-1，[a,b]為z方向的物理區(qū)間(a、b可以為有界的任意實數(shù))，通過坐標(biāo)轉(zhuǎn)換 映射至計算坐標(biāo)z’，其中z’屬于[-1,1]，然后通過公式(11)求解z方向的應(yīng)力分量的一階微分，通過公式(13)求解速度在z方向的一階微分，其中，公式(12)、(14)表示 的離散表達(dá)式。
[0100] 求得 和 后，公式(6)右邊項的值就可以得到了(fy，ρ和μ是已知的)，于是公式(6)由偏微分方程組變?yōu)榱顺Ｎ⒎址匠探M，再利用數(shù)值積分技術(shù)(如歐拉法，4階龍格庫塔法等)求得下一時刻的V、τxy和τyz。
[0101] 施加物理邊界條件和數(shù)值邊界條件，針對某一時刻點，對二維x-z空間的邊界進(jìn)行邊界方程(15)求解，獲得下一時刻點的邊界上的位移和應(yīng)力。
[0102] 隨時間步進(jìn)，獲得每個時刻點的內(nèi)部網(wǎng)格點的速度、應(yīng)力和邊界上的位移、應(yīng)力。根據(jù)每個時刻點的內(nèi)部網(wǎng)格點的速度、應(yīng)力和邊界上的位移、應(yīng)力，獲得SH波的波場記錄，得到空間各場量移隨時間的變化關(guān)系。
[0103] 圖5a、5b、5c分別示出了根據(jù)本發(fā)明的一個實施例的模擬得到的檢波器R1、R2、R3的位移記錄和理論結(jié)果的對比圖，由圖可知，數(shù)值解和理論解符合得很好，證明本方法具有很高的模擬精度。
[0104] 綜上所述，本發(fā)明通過傅里葉偽譜法計算水平導(dǎo)數(shù)，用切比雪夫偽譜法計算垂直方向?qū)?shù)，在解決自由面模擬的同時，保證了模擬的高精度。
[0105] 本領(lǐng)域技術(shù)人員應(yīng)理解，上面對本發(fā)明的實施例的描述的目的僅為了示例性地說明本發(fā)明的實施例的有益效果，并不意在將本發(fā)明的實施例限制于所給出的任何示例。
[0106] 根據(jù)本發(fā)明的一種地震波模擬系統(tǒng)，其上存儲有計算機(jī)程序，所述程序被處理器執(zhí)行時實現(xiàn)以下步驟：步驟1：建立含自由表面的地層介質(zhì)中的地震波波動方程，以g表示其中任意一個變量；步驟2：對三維空間進(jìn)行網(wǎng)格剖分，其中x方向采用等距網(wǎng)格，z方向采用高斯-切比雪夫-洛巴托點獲得空間離散點；步驟3：針對某一時刻點的內(nèi)部網(wǎng)格點，x方向采用傅里葉方程求解下一時刻點的內(nèi)部網(wǎng)格點的變量g，z方向采用采用切比雪夫方程求解下一時刻點的內(nèi)部網(wǎng)格點的變量g；步驟4：施加物理邊界條件和數(shù)值邊界條件，針對某一時刻點進(jìn)行邊界方程求解，獲得下一時刻點的邊界上的變量g；步驟5：隨時間步進(jìn)，重復(fù)步驟3-4，直至獲得每個時刻點的內(nèi)部網(wǎng)格點和邊界上的變量g。
[0107] 在一個示例中，x方向采用傅里葉方程求解變量g的一階微分包括：
[0108]
[0109] 其中，k代表波數(shù)，Δk為波數(shù)間隔，Δk＝2π/(NΔx)，i是虛數(shù)單位，G(lΔk)表示函數(shù)g(nΔx)的傅里葉變換，g(nΔx)為g在x方向的離散表達(dá)式。
[0110] 在一個示例中，g(nΔx)的表達(dá)式為：
[0111]
[0112] 式中，n,l為節(jié)點索引，N為傅里葉多項式的最高階數(shù)。
[0113] 在一個示例中，z方向的高斯-切比雪夫-洛巴托離散點的表達(dá)式為：
[0114]
[0115] 其中，計算區(qū)間z被限制在[-1,1]之間，端點z0＝1，zM＝-1，M是運(yùn)算中所取切比雪夫多項式的階數(shù)。
[0116] 在一個示例中，z方向采用切比雪夫方程求解變量g的一階微分包括：
[0117]
[0118] 其中， 的離散表達(dá)式為
[0119]
[0120] 其中，Dc為(M+1)×(M+1)的矩陣，
[0121] 本系統(tǒng)通過傅里葉偽譜法計算水平導(dǎo)數(shù)，用切比雪夫偽譜法計算垂直方向?qū)?shù)，在解決自由面模擬的同時，保證了模擬的高精度。
[0122] 以上已經(jīng)描述了本發(fā)明的各實施例，上述說明是示例性的，并非窮盡性的，并且也不限于所披露的各實施例。在不偏離所說明的各實施例的范圍和精神的情況下，對于本技術(shù)領(lǐng)域的普通技術(shù)人員來說許多
修改和變更都是顯而易見的。
