本發(fā)明涉及一種波動方程地震響應(yīng)模擬方法，尤其涉及一種時(shí)間域聲波方程顯式有限差分地震響應(yīng)模擬方法。

背景技術(shù)：

波動方程正演模擬是勘探地球物理領(lǐng)域重要的研究方向，對波動現(xiàn)象的模擬、逆時(shí)偏移和全波形反演的研究和實(shí)施具有重要的意義。目前，波動方程正演模擬方法主要有有限差分方法、有限元方法、邊界元法和偽譜法等。在諸多正演方法中，有限差分方法相比于其他方法，具有更小的計(jì)算復(fù)雜度和內(nèi)存需求，能夠滿足逆時(shí)偏移和全波形反演對計(jì)算速度和內(nèi)存占用的要求。因此，有限差分方法是當(dāng)前逆時(shí)偏移和全波形反演中應(yīng)用最廣泛的方法。

有限差分方法通過差分近似代替偏導(dǎo)數(shù)，從而存在精度誤差。時(shí)間上的精度誤差稱為時(shí)間頻散，通常表現(xiàn)為相位超前；空間上的精度誤差稱為空間頻散，通常表現(xiàn)為相位滯后。時(shí)間二階精度和空間任意偶數(shù)階精度差分形式廣泛使用于波動方程模擬，要求時(shí)間步長和網(wǎng)格空間間隔足夠小，從而避免時(shí)間頻散和空間頻散。

目前有限差分系數(shù)的求取主要有以下幾類：一是基于泰勒展開的有限差分系數(shù)求取，該方法獨(dú)立求取時(shí)間和空間方向的差分系數(shù)，精度較低。二是基于時(shí)-空域頻散關(guān)系的有限差分系數(shù)求取，采用常規(guī)的差分網(wǎng)格策略，該方法可以實(shí)現(xiàn)部分方向的時(shí)間高階精度，相比較泰勒展開方法，該方法沒有增加額外的計(jì)算量。三是基于最優(yōu)化方法的有限差分系數(shù)求取，該方法對高頻成分仍能實(shí)現(xiàn)較高的精度，但是計(jì)算量較大，某些情況下容易陷入局部最優(yōu)，無法取得全局最優(yōu)值。

波動方程正演是逆時(shí)偏移和全波形反演的主要環(huán)節(jié)，實(shí)現(xiàn)波動方程的高效率、高精度的數(shù)值模擬，對逆時(shí)偏移和全波形反演具有重要的意義。由于逆時(shí)偏移和全波形反演需要大量的正演計(jì)算，當(dāng)前有限差分方法計(jì)算量和內(nèi)存需求仍然巨大，計(jì)算速度和計(jì)算效率低，限制了逆時(shí)偏移和全波形反演的廣泛應(yīng)用，因此高精度和高效率的有限差分正演是逆時(shí)偏移和全波形反演廣泛實(shí)施的關(guān)鍵。

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

針對上述問題，本發(fā)明的目的是提供一種時(shí)間域聲波方程顯式有限差分地震響應(yīng)模擬方法，基于k空間算子補(bǔ)償?shù)目臻g-波數(shù)域函數(shù)擬合方法求取差分系數(shù)，同時(shí)針對最大擬合誤差，相對應(yīng)地提出了自適應(yīng)空間差分算子長度，能夠?qū)崿F(xiàn)高精度、高效率的時(shí)間域有限差分?jǐn)?shù)值模擬。

為實(shí)現(xiàn)上述目的，本發(fā)明采取以下技術(shù)方案：一種時(shí)間域聲波方程顯式有限差分地震響應(yīng)模擬方法，包括以下步驟：

1)設(shè)定空間差分算子長度2m+1的初始值和允許的最大近似誤差ξ；

2)根據(jù)目標(biāo)區(qū)域的聲波速度，基于空間-波數(shù)域函數(shù)近似的空間有限差分系數(shù)求取方法，求取基于k空間算子補(bǔ)償數(shù)值模擬誤差建立的顯式有限差分公式的差分系數(shù)和波數(shù)域誤差函數(shù)值；

3)如果波數(shù)域誤差函數(shù)值大于允許的最大近似誤差ξ，則將m的值加1，得到新的空間差分算子長度2m+1，返回步驟2)；如果波數(shù)域誤差函數(shù)值小于允許的最大近似誤差ξ，且m>1，則將m的值減1，得到新的空間差分算子長度2m+1，返回步驟2)；否則，確定滿足允許的最大近似誤差ξ的最小空間差分算子長度2m+1及其對應(yīng)的差分系數(shù)，繼續(xù)下一步；

4)根據(jù)滿足允許的最大近似誤差ξ的最小空間差分算子長度2m+1的值及其對應(yīng)的差分系數(shù)，對目標(biāo)區(qū)域進(jìn)行地震響應(yīng)數(shù)值模擬，得到目標(biāo)區(qū)域的模擬地震記錄或地震波場。

所述步驟2)中基于k空間算子補(bǔ)償數(shù)值模擬誤差建立的顯式有限差分公式為：

式中，p(x,t)為壓力波場；x表示笛卡爾空間坐標(biāo)系統(tǒng)；t為時(shí)間；δt為時(shí)間間隔；c(x,δxi)為加權(quán)系數(shù)；δxi為位置x和第i個(gè)差分點(diǎn)之間的距離，δx0＝0；m為空間差分算子長度的一半，為整數(shù)；n是模型維度；k為波數(shù)矢量，k＝|k|；e(x,k)為x處波數(shù)域誤差函數(shù)；為壓力波場p(x,t)的傅里葉響應(yīng)；e^ik·x為空間反傅里葉變換的一部分。

所述步驟2)中基于空間-波數(shù)域函數(shù)近似的空間有限差分系數(shù)求取方法求取差分系數(shù)的方程為：

a^twac＝a^tws；

其中，

式中，t表示矩陣轉(zhuǎn)置；n為離散波數(shù)的數(shù)量；ν為聲波速度；sinc是一個(gè)算子，在函數(shù)數(shù)學(xué)上表示為sinc(x)＝sin(x)/x；ki為第i個(gè)波數(shù)分量；w為相對權(quán)重系數(shù)；a(ki)為波數(shù)分量ki的振幅值；num(ki)是波數(shù)分量ki的出現(xiàn)次數(shù)；kmax為每個(gè)維度方向的最大波數(shù)；是波數(shù)的最大值。

本發(fā)明由于采取以上技術(shù)方案，其具有以下優(yōu)點(diǎn)：1、本發(fā)明的一種時(shí)間域聲波方程顯式有限差分地震響應(yīng)模擬方法，通過采用自適應(yīng)空間差分算子長度和基于k空間算子的空間-波數(shù)域函數(shù)擬合方法求取差分系數(shù)，可以高效率、高精度地模擬實(shí)際地質(zhì)模型中的地震響應(yīng)，相比較常規(guī)的顯式有限差分方法，具有更高的時(shí)間精度，模擬中可以采用更大的時(shí)間步長，計(jì)算速度更快。2、本發(fā)明的一種時(shí)間域聲波方程顯式有限差分地震響應(yīng)模擬方法，可以廣泛應(yīng)用在均勻模型、層狀模型和實(shí)際地質(zhì)模型中。

附圖說明

圖1(a)是均勻介質(zhì)模型中采用本發(fā)明方法且δt＝2ms時(shí)得到的波場快照圖；

圖1(b)、(c)、(d)是均勻介質(zhì)模型中采用基于頻散關(guān)系的有限差分方法且時(shí)間間隔分別為δt＝2ms、δt＝1ms、δt＝0.25ms時(shí)得到的波場快照圖；

圖2(a)是2d層狀介質(zhì)中的速度模型示意圖；

圖2(b)是2d層狀介質(zhì)中采用基于頻散關(guān)系的有限差分方法得到的波場快照圖；

圖2(c)、(d)是2d層狀介質(zhì)中采用本發(fā)明方法以及分別使用定長度算子和變長度算子得到的波場快照圖；

圖2(e)、(f)分別是對比道取自不同位置的圖2(b)、(c)、(d)的波場快照單道對比圖；

圖3是2dsigsbee2模型示意圖；

圖4(a)是2dsigsbee2模型中參考波場示意圖；

圖4(b)、(c)是2dsigsbee2模型中采用基于頻散關(guān)系的有限差分方法和分別使用2ms和0.5ms時(shí)間步長得到的波場快照圖；

圖4(d)、(e)是2dsigsbee2模型中采用本發(fā)明方法、時(shí)間間隔為2ms以及分別采用定長度算子和變長度算子得到的波場快照圖；

圖5(a)、(b)、(c)和(d)分別是與圖4(b)、(c)、(d)和(e)相對應(yīng)的地震記錄對比圖。

具體實(shí)施方式

下面結(jié)合附圖和實(shí)施例對本發(fā)明進(jìn)行詳細(xì)的描述。

本發(fā)明方法包括基于k空間算子補(bǔ)償?shù)挠邢薏罘窒禂?shù)的求取和新的自適應(yīng)空間差分算子長度策略。

首先，基于k空間算子補(bǔ)償數(shù)值模擬誤差建立顯式有限差分公式。

均勻無吸收衰減介質(zhì)中，時(shí)間-空間域二階常密度聲波方程為：

式中，p(x,t)為壓力波場；x表示笛卡爾空間坐標(biāo)系統(tǒng)；t為時(shí)間；ν為介質(zhì)速度；為空間拉普拉斯算子。

通過空間傅里葉變換，可以將時(shí)間-空間域二階常密度聲波方程轉(zhuǎn)換為時(shí)間-波數(shù)域二階常密度聲波方程：

式中，為壓力波場p(x,t)的傅里葉響應(yīng)；k為波數(shù)矢量，k＝|k|。

采用二階時(shí)間離散表示時(shí)間導(dǎo)數(shù)，可以得到時(shí)間二階和空間準(zhǔn)確的半差分方程：

式中，δt為時(shí)間間隔。

常規(guī)的聲波方程采用二階精度中心差分離散時(shí)間導(dǎo)數(shù)，對該離散采用二階k空間算子進(jìn)行補(bǔ)償，可以有效提高時(shí)間精度。其中，二階k空間算子為sinc²(νkδt/2)。因此，在均勻介質(zhì)中，可以采用二階k空間算子補(bǔ)償時(shí)間精度，從而得到時(shí)間和空間準(zhǔn)確的方程為：

式中，sinc是一個(gè)算子，在函數(shù)數(shù)學(xué)上表示為sinc(x)＝sin(x)/x。

采用空間傅里葉變換，可以得到如下的時(shí)間域波場更新公式：

式中，e^ik·x為空間反傅里葉變換的一部分。

與有限差分方法相比，基于傅里葉變換的數(shù)值模擬方法計(jì)算量大。為避免使用空間傅里葉變換，本發(fā)明從時(shí)間域波場更新公式出發(fā)，采用中心對稱差分樣板，可以將位置x處的空間2m^th精度差分方程表示為：

式中，m為正整數(shù)，2m+1為空間算子長度，精度為2m；n是模型維度；δxi為位置x和第i個(gè)差分點(diǎn)之間的距離，δx0＝0；e(x,t)為時(shí)空域的誤差函數(shù)；c(x,δxi)為加權(quán)系數(shù)。

利用空間傅里葉變換，位置x處的空間2m^th精度差分方程可以轉(zhuǎn)換為顯式有限差分公式：

式中，e(x,k)為x處波數(shù)域的誤差函數(shù)。

對比準(zhǔn)確的時(shí)間域波場更新公式和轉(zhuǎn)換后位置x處的空間2m^th精度差分方程，可以得到：

在空間2m^th精度差分方程中，在每個(gè)時(shí)間點(diǎn)，波場更新為逐點(diǎn)更新，局部速度條件即可以滿足空間2m^th精度差分方程的要求，因此本方法也適用于復(fù)雜介質(zhì)模型。

然后，空間-波數(shù)域函數(shù)近似方法求取差分系數(shù)的一般公式為：

對于給定的網(wǎng)格剖分和時(shí)間步長，空間差分系數(shù)可以通過最小化網(wǎng)格點(diǎn)處的近似誤差得到。在本實(shí)施例中可以采用二范數(shù)最小化近似誤差：

則加權(quán)系數(shù)c(x,δx0)可以通過以下矩陣-向量方程求得：

式中，n為離散波數(shù)的數(shù)量。

上述矩陣-向量方程可以簡寫為：

an×nmcnm×1＝sn×1

通常情況下，矩陣-向量方程是一個(gè)超定方程，因?yàn)殡x散波數(shù)的數(shù)量遠(yuǎn)大于差分系數(shù)的個(gè)數(shù)，即n>nm。為了提高求解精度，本發(fā)明引入相對權(quán)重系數(shù)：

式中，ki為第i個(gè)波數(shù)分量；a(ki)為波數(shù)分量ki的振幅值；num(ki)是波數(shù)分量ki的出現(xiàn)次數(shù)；kmax為每個(gè)維度方向的最大波數(shù)；因?yàn)槊總€(gè)方向的最大波數(shù)為kmax，所以波數(shù)的最大值是

通過引入權(quán)重系數(shù)，差分系數(shù)可以通過以下方程求取：

a^twac＝a^tws

式中，t表示矩陣轉(zhuǎn)置。

最后，基于近似誤差的自適應(yīng)空間差分算子長度策略：

在有限差分?jǐn)?shù)值中，波場外推和差分系數(shù)計(jì)算均和空間差分算子長度成正比，空間算子長度越大，計(jì)算量越大。為了減小計(jì)算量，本發(fā)明采用自適應(yīng)差分算子長度策略?？臻g差分算子長度可以通過下列公式求得：

max(|e(x,k,m)|)<ξwhenk≤kmax

式中，ξ為允許的最大近似誤差。

通過實(shí)際試驗(yàn)得到，最大擬合誤差通常在最大波數(shù)處取得，因此，差分系數(shù)可以通過以下公式求得，從而進(jìn)一步減小計(jì)算量：

max(|e(x,kmax,m)|)<ξ

基于上述原理，本發(fā)明提供的一種時(shí)間域聲波方程顯式有限差分地震響應(yīng)模擬方法，具體包括以下步驟：

1)設(shè)定空間差分算子長度2m+1的初始值和允許的最大近似誤差ξ；

2)根據(jù)目標(biāo)區(qū)域的聲波速度，基于空間-波數(shù)域函數(shù)近似的空間有限差分系數(shù)求取方法，求取基于k空間算子補(bǔ)償數(shù)值模擬誤差建立的顯式有限差分公式的差分系數(shù)和波數(shù)域誤差函數(shù)值(即擬合誤差值)。

3)如果波數(shù)域誤差函數(shù)值大于允許的最大近似誤差ξ，則將m的值加1，得到新的空間差分算子長度2m+1，返回步驟2)；如果波數(shù)域誤差函數(shù)值小于允許的最大近似誤差ξ，且m>1，則將m的值減1，得到新的空間差分算子長度2m+1，返回步驟2)；否則，確定滿足允許的最大近似誤差ξ的最小m的值及其對應(yīng)的空間差分算子長度2m+1和差分系數(shù)，繼續(xù)下一步。

4)根據(jù)滿足允許的最大近似誤差ξ的最小的空間差分算子長度2m+1的值及其對應(yīng)的差分系數(shù)，對目標(biāo)區(qū)域進(jìn)行地震響應(yīng)數(shù)值模擬，得到目標(biāo)區(qū)域的模擬地震記錄或地震波場。

本發(fā)明的一種時(shí)間域聲波方程顯式有限差分地震響應(yīng)模擬方法，可以廣泛應(yīng)用在均勻模型、層狀模型和實(shí)際地質(zhì)模型中。

在均勻介質(zhì)模型中，分別采用本發(fā)明的方法和基于頻散關(guān)系的有限差分方法進(jìn)行地震響應(yīng)模擬，模擬參數(shù)為：v＝1500m/s，h＝17.5m，m＝4；時(shí)間間隔分別為：δt＝2ms、δt＝1ms、δt＝0.25ms；震源位于模型中心，模擬過程中沒有設(shè)置吸收邊界。得到的波場快照如圖1(a)、(b)、(c)和(d)所示。從圖中可以看出，相比較基于時(shí)-空域頻散關(guān)系方法，本發(fā)明方法的精度更高，當(dāng)時(shí)間步長較大時(shí)，仍能取得很好的模擬效果。

在2d層狀介質(zhì)中，分別采用基于頻散關(guān)系的有限差分方法，以及本發(fā)明的方法分別采用不同的空間差分算子長度進(jìn)行地震響應(yīng)模擬，采用的速度模型如圖2(a)所示?？臻g差分算子長度和相對應(yīng)的cpu運(yùn)算時(shí)間如下表1所示，得到的波場快照如圖2(b)、(c)和(d)所示；波場快照單道對比如圖2(e)、(f)所示，其中，道1、2、3分別來自圖2(b)、(c)、(d)，道4為道1和2的差值，道5為道1和2的差值。從圖中可以看出，相比較基于時(shí)-空域頻散關(guān)系方法，本發(fā)明方法的數(shù)值精度更高，采用空間自適應(yīng)差分算子，仍能取得較好地?cái)?shù)值精度。

表1:2d水平層狀介質(zhì)模擬采用的空間算子長度和cpu運(yùn)算時(shí)間

如圖3所示，為2d的sigsbee2模型，圖中大黑點(diǎn)表示炮點(diǎn)位置，小黑點(diǎn)表示檢波點(diǎn)位置。在2dsigsbee2模型中，采用基于頻散關(guān)系的有限差分方法分別采用不同的時(shí)間步長，以及采用本發(fā)明的方法分別采用不同長度的空間差分算子進(jìn)行地震響應(yīng)模擬；基于頻散關(guān)系的有限差分方法采用的時(shí)間步長為2m/s和0.5m/s，本發(fā)明的方法采用的時(shí)間間隔為2m/s；采用的參考波場如圖4(a)所示?？臻g差分算子長度和相對應(yīng)的cpu運(yùn)算時(shí)間如下表2所示，得到的波場快照如圖4(b)、(c)、(d)和(e)所示；相對應(yīng)地得到的地震記錄對比如圖(b)、(c)、(d)和(e)所示。此處圖4(a)、(b)、(c)、(d)和(e)中類似馬賽克的部分為真實(shí)的波場模擬結(jié)果，通過數(shù)值模擬得到的結(jié)果如此，沒有經(jīng)過后續(xù)處理。從圖中可以看出，相比較基于時(shí)-空域頻散關(guān)系方法，本方法的數(shù)值精度更高。采用空間自適應(yīng)差分算子，仍能取得較好地?cái)?shù)值精度。自適應(yīng)差分算子沒有明顯降低數(shù)值精度。

表2:sigsbee2模型數(shù)值模擬采用的空間算子長度和cpu運(yùn)算時(shí)間

上述各實(shí)施例僅用于說明本發(fā)明，其中各部件的結(jié)構(gòu)、設(shè)置位置及其連接方式等都是可以有所變化的，凡是在本發(fā)明技術(shù)方案的基礎(chǔ)上進(jìn)行的等同變換和改進(jìn)，均不應(yīng)排除在本發(fā)明的保護(hù)范圍之外。