標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值解的計算方法及系統(tǒng)的制作方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及地震勘探基礎(chǔ)應(yīng)用技術(shù)領(lǐng)域，尤其涉及一種標(biāo)準(zhǔn)線性固體模型的有限 差分解的穩(wěn)定性條件數(shù)值解的計算方法及系統(tǒng)。
【背景技術(shù)】
[0002] 地震數(shù)值模擬是地震勘探和地震學(xué)的重要基礎(chǔ)，同時也是了解復(fù)雜介質(zhì)中地震波 傳播規(guī)律的重要工具，其作用貫穿于整個地震采集、處理和解釋中。隨著地震勘探開發(fā)的深 入，常規(guī)的彈性介質(zhì)理論難W滿足實際介質(zhì)需求。地震波在實際地層中傳播時，能量和相位 都發(fā)生改變，直接影響地震資料的分辨率，實際地層表現(xiàn)出一定的黏彈特性，因此通過對黏 彈介質(zhì)中的地震波進行數(shù)值模擬，來研究和分析地震波傳播過程中的衰減特征，對實際地 震資料分辨率的提高非常有意義。但是確保黏彈介質(zhì)數(shù)值模擬穩(wěn)定（即時間步長的確定） 是一個急需解決的問題，送關(guān)系到模擬算法的成敗，也關(guān)系到模擬效率的高低。
[0003] 目前主要運用黏彈固體模型來表征實際介質(zhì)的粘滯特性，應(yīng)用最廣的黏彈固體模 型主要包括；Kelvin-Voigt固體模型、Maxwell固體模型和標(biāo)準(zhǔn)線性固體模型（也稱為標(biāo)準(zhǔn) 線性黏彈固體模型）。其中，Kelvin-Voigt固體模型不能考慮應(yīng)力作用下應(yīng)變的突然變化， 也不能表示應(yīng)力消失后的剩余應(yīng)變，Maxwell固體模型不具備蠕變特征，兩者都不足W描述 大多數(shù)黏彈介質(zhì)的特征，但是標(biāo)準(zhǔn)線性固體模型可W同時考慮具備應(yīng)變突然變化和剩余應(yīng) 變及蠕動特征，因此，標(biāo)準(zhǔn)線性固體模型能更加全面表征實際地層的黏彈特性，相比較而言 更符合實際情況，所W進行標(biāo)準(zhǔn)線性固體模型的黏彈介質(zhì)數(shù)值模擬是非常有必要的。
[0004] 現(xiàn)有技術(shù)中標(biāo)準(zhǔn)線性固體模型的黏彈介質(zhì)數(shù)值模擬主要是全Η維、全波場的模 擬，此種方法的計算量非常大，從而產(chǎn)生的費用（主要體現(xiàn)在電力、存儲設(shè)備、人力等資源 的消耗）是巨大的，數(shù)值模擬的效率低，因此有必要提供一種計算方法定量化地給出標(biāo)準(zhǔn) 線性固體模型數(shù)值模擬時間步長（即標(biāo)準(zhǔn)線性固體模型的有限差分解的穩(wěn)定性條件數(shù)值 解），從而更高效、更成功地指導(dǎo)數(shù)值模擬的進行。

【發(fā)明內(nèi)容】

[0005] 本發(fā)明所要解決的技術(shù)問題是現(xiàn)有技術(shù)中采用全Η維、全波場的模擬方法產(chǎn)生的 費用高、數(shù)值模擬效率低的缺陷。
[0006] 為了解決上述技術(shù)問題，本發(fā)明提供了一種標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值 解的計算方法及系統(tǒng)。
[0007] 本發(fā)明的技術(shù)方案為：
[0008] -種標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值解的計算方法，包括：
[0009] 構(gòu)建標(biāo)準(zhǔn)線性固體模型，并使所述標(biāo)準(zhǔn)線性固體模型包括彼此串聯(lián)的第一彈性體 和第二彈性體、W及與所述第一彈性體并聯(lián)的阻尼器；
[0010] 確定所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條件的狀態(tài)傳遞矩陣；
[0011] 獲取多組參數(shù)集，并使每組參數(shù)集中包括所述第一彈性體的彈性系數(shù)、所述第二 彈性體的彈性系數(shù)、所述阻尼器的黏滯系數(shù)、頻率和介質(zhì)密度；
[0012] 在設(shè)定的空間差分精度和空間網(wǎng)格步長下，依次利用各組所述參數(shù)集并通過計算 機計算使得所述狀態(tài)傳遞矩陣的特征值的模小于1的時間步長；
[0013] 確定計算得出的時間步長為標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值解。
[0014] 優(yōu)選的是，所述確定所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條件的狀態(tài)傳遞 矩陣包括：
[0015] 確定所述標(biāo)準(zhǔn)線性固體模型的本構(gòu)方程
其中P為總應(yīng)力，ε為總應(yīng)變，彈性系數(shù)為Ml的 所述第一彈性體的應(yīng)變與黏滯系數(shù)為M2的所述阻尼器的應(yīng)變相等，Ms為所述第二彈性體的 彈性系數(shù)；
[0016] 根據(jù)總應(yīng)變ε與質(zhì)點位移（u，v，w)間的關(guān)系方程
W及所述本 構(gòu)方程，得到第一方程：
[0017]
[0018] 對所述第一方程的左右兩邊分別對時間求二次偏導(dǎo)數(shù)，得到第二方程：
[0019]
[0020] 利用所述第二方程和聲波的納維爾方程
得到第Η方程：
?
[0021] 其中Ρ為所述介質(zhì)密 度；
[0022] 將所述第S方程中的應(yīng)力取空間傅里葉變換，得到第四方程：
其中#為總應(yīng)力Ρ的空間傅里葉變換， k為波數(shù)；
[0023] 對所述第四方程中的時間偏導(dǎo)數(shù)用差分近似，得到第五方程：
[0024]
其 中護聲"、；ri分別為第n-2, n-1，η, n+1時刻的夢值，并且所述波數(shù)k滿足：
在所述空間差分精度為2N的情況下，X，y，Z Η個方向上的空 間網(wǎng)格步長分別為Δχ，Ay, Δζ，曰1為對應(yīng)所述空間差分精度2Ν的空間差分系數(shù)，At為 所述時間步長；
[0025] 根據(jù)所述第五方程，得到所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條件的狀態(tài) 傳遞矩陣
其中：
[0029] 優(yōu)選的是，所述方法還包括；在計算得出的所述時間步長后，根據(jù)所述參數(shù)集計算 所述標(biāo)準(zhǔn)線性固體模型的品質(zhì)因子，并根據(jù)所述參數(shù)集和所述品質(zhì)因子計算所述標(biāo)準(zhǔn)線性 固體模型的介質(zhì)速度。
[0030] 優(yōu)選的是，在設(shè)定的空間差分精度和空間網(wǎng)格步長下，依次利用各組所述參數(shù)集， 運用安裝在所述計算機上的Matlab仿真軟件編程計算使得所述狀態(tài)傳遞矩陣的特征值的 模小于1的時間步長。
[0031] 一種標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值解的計算系統(tǒng)，包括：
[0032] 模型構(gòu)建單元，用于構(gòu)建標(biāo)準(zhǔn)線性固體模型，并使所述標(biāo)準(zhǔn)線性固體模型包括彼 此串聯(lián)的第一彈性體和第二彈性體、W及與所述第一彈性體并聯(lián)的阻尼器；
[0033] 狀態(tài)傳遞矩陣確定單元，用于確定所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條 件的狀態(tài)傳遞矩陣；
[0034] 參數(shù)集獲取單元，用于獲取多組參數(shù)集，并使每組參數(shù)集中包括所述第一彈性體 的彈性系數(shù)、所述第二彈性體的彈性系數(shù)、所述阻尼器的黏滯系數(shù)、頻率和介質(zhì)密度；
[0035] 時間步長計算單元，用于在設(shè)定的空間差分精度和空間網(wǎng)格步長下，依次利用各 組所述參數(shù)集并通過計算機計算使得所述狀態(tài)傳遞矩陣的特征值的模小于1的時間步長；
[0036] 穩(wěn)定性條件數(shù)值解確定單元，用于確定計算得出的時間步長為標(biāo)準(zhǔn)線性固體模型 的穩(wěn)定性條件數(shù)值解。
[0037] 優(yōu)選的是，所述狀態(tài)傳遞矩陣確定單元包括：
[003引本構(gòu)方程確定單元，用于確定所述標(biāo)準(zhǔn)線性固體模型的本構(gòu)方程
其中P為總應(yīng)力，ε為總應(yīng)變，彈性系數(shù)為Ml的 所述第一彈性體的應(yīng)變與黏滯系數(shù)為M2的所述阻尼器的應(yīng)變相等，Ms為所述第二彈性體的 彈性系數(shù)；
[0039] 第一方程確定單元，用于根據(jù)總應(yīng)變ε與質(zhì)點位移（u，v，w)間的關(guān)系方程
W及所述本構(gòu)方程確定單元確定的本構(gòu)方程，得到第一方程：
[0040]
[0041] 第二方程確定單元，用于對所述第一方程確定單元得到的第一方程的左右兩邊分 別對時間求二次偏導(dǎo)數(shù)，得到第二方程：
[0042]
[0043] 第Η方程確定單元，用于利用所述第二方程確定單元得到的第二方程和聲波的納 維爾方程
得到第Η方程：
>
[0044] 其中Ρ為所述介質(zhì)密 度；
[0045] 第四方程確定單元，用于將所述第Η方程確定單元得到的第Η方程中的應(yīng)力取空 間傅里葉變換，得到第四方程：
其中# 為總應(yīng)力Ρ的空間傅里葉變換，k為波數(shù)；
[0046] 第五方程確定單元，用于對所述第四方程確定單元得到的第四方程中的時間偏導(dǎo) 數(shù)用差分近似，得到第五方程：
[0047]
其中r-氣護1.護擴'' 分別為第η-化-lnn+1時刻的#值，并且所述波數(shù)k滿足：
在所述空間差分精度為2N的情況下，X，y，Z Η個方向上的空 間網(wǎng)格步長分別為Δχ，Ay, Δζ，曰1為對應(yīng)所述空間差分精度2Ν的空間差分系數(shù)，At為 所述時間步長；
[0048] 狀態(tài)傳遞矩陣子確定單元，用于根據(jù)所述第五方程確定單元得到的第五方程，得 到所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條件的狀態(tài)傳遞矩陣
癢中：
[0052] 優(yōu)選的是，所述系統(tǒng)還包括：
[0053] 品質(zhì)因子計算單元，用于在所述時