宋利偉 石穎 陳樹民 柯璇 侯曉慧 劉志奇

1) (東北石油大學(xué)物理與電子工程學(xué)院, 大慶 163318)

2) (東北石油大學(xué)地球科學(xué)學(xué)院, 大慶 163318)

3) (大慶油田有限責任公司勘探開發(fā)研究院, 大慶 163712)

實際介質(zhì)普遍具有黏彈性, 波在傳播過程中常伴有能量的耗散、相位畸變和頻帶變窄等, 對于含有液體和氣體的介質(zhì), 衰減現(xiàn)象尤為突出. 由于經(jīng)典的波動理論未考慮介質(zhì)的黏彈性效應(yīng), 基于完全彈性假設(shè)的模擬波場和實際傳播特征之間的差異明顯, 波動理論在工程技術(shù)中的應(yīng)用效果還有提升的空間. 在巖石物理中,品質(zhì)因子Q是量化地震衰減強度的參數(shù), 為了研究波在地下介質(zhì)中的傳播規(guī)律, 本文從常Q理論出發(fā), 在黏彈性介質(zhì)頻散關(guān)系中, 利用多項式擬合和Taylor展開法將頻率的分數(shù)階轉(zhuǎn)化為整數(shù)階, 進而推導(dǎo)了時間域復(fù)數(shù)形式的地下黏彈性介質(zhì)波動方程. 該近似處理避免了頻散關(guān)系經(jīng)域轉(zhuǎn)換后出現(xiàn)分數(shù)階時間微分項, 能有效地降低計算成本. 最后, 采用有限差分法聯(lián)合偽譜法對均質(zhì)模型實現(xiàn)了波場的數(shù)值模擬, 驗證了方程的有效性.

1 引 言

地下介質(zhì)普遍具有黏彈性, 波在傳播過程中部分能量轉(zhuǎn)變成熱, 波場特征會發(fā)生明顯變化, 如振幅衰減、相位畸變和頻帶變窄等. 當?shù)貙又泻袣怏w或液體時, 黏彈性效應(yīng)更加突出. 忽略黏彈性的波場模擬易導(dǎo)致波場描述失真, 不利于工程技術(shù)中的應(yīng)用[1-4], 因此有必要開展黏彈性介質(zhì)波動理論及波場數(shù)值模擬的研究.

在巖石物理中, 品質(zhì)因子Q用于計算應(yīng)力與應(yīng)變之間的相位差, 是度量地震波衰減特征的參數(shù). Liu等[5]利用廣義標準線性體建立了不依賴頻率的Q模型, 解釋了地震波能量衰減和相位畸變的原因, 理論上Q值與實際介質(zhì)的匹配度取決于線性體的數(shù)目和組合方式[6,7]. Kjartansson[8]在衰減系數(shù)和頻率成線性關(guān)系的假設(shè)下建立了常Q模型, 由于頻散關(guān)系的數(shù)學(xué)形式簡單, 在地震資料處理中得到了廣泛關(guān)注. 從常Q模型的本構(gòu)關(guān)系出發(fā), Carcione等[9,10]推導(dǎo)了黏彈性介質(zhì)波動方程,并采用Grünwald-Letnikov近似實現(xiàn)了地震波場數(shù)值模擬, 由于方程中分數(shù)階時間微分項的計算依賴于歷史時刻波場, 因此需要保存全時程波場, 計算成本較高. 雖然短時截斷法和內(nèi)部變量法能降低內(nèi)存的占用[11,12], 但依然難以應(yīng)對大尺度模型.Yang和Zhu[13]推導(dǎo)了分數(shù)階拉普拉斯算子的黏彈性介質(zhì)波動方程, 并應(yīng)用于地下構(gòu)造成像中.

波動方程是雙曲線偏微分方程[14-17], 常用的數(shù)值模擬方法有偽譜法、有限元法和有限差分法.有限差分法無需數(shù)據(jù)域的轉(zhuǎn)換, 且計算效率高, 是波動方程數(shù)值模擬中最常用的算法[18,19], 為了在相同計算量的情況下提高數(shù)值模擬精度, 先后發(fā)展了交錯網(wǎng)格和旋轉(zhuǎn)交錯網(wǎng)格[20,21]. 偽譜法通過Fourier變換實現(xiàn), 相當于無限階差分算子精度, 計算效率不及有限差分法[22]. 有限元法的優(yōu)勢是網(wǎng)格剖分靈活, 能模擬復(fù)雜介質(zhì)中的波傳播, 但計算成本較高[23].

本文基于Kjartansson提出的常Q模型, 建立了黏彈性介質(zhì)波動方程, 研究了波動方程的數(shù)值模擬方法, 通過數(shù)值實驗分析了品質(zhì)因子Q對波場振幅、相位和頻帶的影響.

2 波動方程建立

在Kjartansson[8]提出的常Q模型中, 模量M(ω)在頻率域的形式為

其中ω0是參考頻率;M0是模量;ω是頻率; i 是 虛數(shù)單位;γ與品質(zhì)因子Q相關(guān), 可通過下式確定[8]:

式中, πγ的物理意義是介質(zhì)應(yīng)力和應(yīng)變之間相位差.

在地下黏彈性介質(zhì)中波速vc與量M(ω) 的關(guān)系為[24]

式中,ρ為介質(zhì)密度. 相速度vp可利用下式計算,

式中,kR是波數(shù)k的實部, 利用等價關(guān)系i=cos(π/2)+isin(π/2)/2 ,k的實部kR和虛部kI根據(jù)(5)式和(6)式確定,

式中, Re和Im分別為取實部和虛部算子;v0=是在參考頻率處的相速度. 將(5)式代入(4)式可得

類似于彈性介質(zhì), 在黏彈性介質(zhì)中波的頻散關(guān)系可表示為

將(1)式和(3)式代入(8)式有

式中,v2=M0/ρ,v是Q趨于無窮大時的相速度,即介質(zhì)趨于彈性介質(zhì)時波的相速度. 若將(9)式從頻率-波數(shù)域變換至?xí)r間-空間域, 可得含有分數(shù)階時間微分的波動方程, 雖然采用Grünwald-Letnikov近似能實現(xiàn)分數(shù)階微分項的計算, 但該方法需保存歷史時間波場. 為此, 將(9)式等號左側(cè)第一項進行多項式擬合有

式中,a1,a2和a3為擬合系數(shù), 表1提供了參考頻率為50 Hz時不同品質(zhì)因子Q對應(yīng)的擬合系數(shù).當Q= 30時, (10)式擬合曲線如圖1所示, 圖中實心點是理論數(shù)值, 實線為擬合曲線, 以上擬合均是將頻率轉(zhuǎn)換至圓頻率所得.

表1 不同Q值擬合系數(shù)Table 1. Fitting coefficients of different Q values.

圖1 擬合曲線 (Q = 30)Fig. 1. Fitting curve (Q = 30).

(9)式中等號左側(cè)第二項中 iω經(jīng)Fourier變換是時間的一階微分算子, 為了使剩余的ω1-2γ易于計算, 先對k2γ-1進行處理,

將(5)式和(6)式代入(11)式得

對于弱衰減介質(zhì)有γ?1 , 利用Taylor展開可將[1-itan(γπ/2)]2γ-1表示為

忽略(13)式的高階無窮小量后代入(12)式有

(4)式的等價形式為

根據(jù)(7)式, 當參考頻率接近于波的頻率時, 有vp≈v0, 將(14)式代入(15)式為

將(10)式和(16)式代入(9)式有

根據(jù)Fourier變換的微分性質(zhì), (18)式給出了算子在頻率-波數(shù)域和時間-空間域的對應(yīng)關(guān)系,F為Fourier變換算子,?是拉普拉斯算子. 將(17)式乘頻率-波數(shù)域波場后, 經(jīng)域轉(zhuǎn)換可得時間-空間域地下黏彈性介質(zhì)波動方程:

(19)式中,u是復(fù)數(shù)波場, 取其實部即物理意義上的波場. 為了在數(shù)值模擬中研究波動方程中各項對波場的主控因素, 可將(19)式改寫為

(20)式中參數(shù)λ和μ的取值為1或0, 兩參數(shù)組合可形成四個方程, 當λ=1,μ=1 時, (20)式和(19)式等價; 當λ=0,μ=0 時, (20)式退化為聲波方程.

3 波場遞推公式

本文采用二階精度有限差分法求解(19)式中的時間微分項, 其離散格式為

式中,n為時間節(jié)點; Δt是時間采樣間隔.

(19)式中的空間微分項采用偽譜法計算, 先通過Fourier變換將時間-空間域波場轉(zhuǎn)換至波數(shù)域, 與波數(shù)k和拉普拉斯算子的階數(shù)作用后, 再經(jīng)反Fourier變換即可實現(xiàn)求解[25,26]. 以上過程可用(23)式表示:

結(jié)合(21)式、(22)式和(23)式, 波動方程(19)式的波場遞推公式為

4 數(shù)值實驗

為研究地下黏彈性介質(zhì)中波的傳播規(guī)律, 設(shè)計水平和垂直方向的距離為3.8 km的均質(zhì)模型, 空間網(wǎng)格間距為10 m, 波在參考頻率50 Hz處的相速度為2216 m/s. 利用主頻是20 Hz的Ricker子波作為激發(fā)震源, 其振幅峰值為1 Pa, 時間采樣間隔1 ms, 子波的函數(shù)形式為

震源位于模型中央位置, 距離震源0.6 km處布置檢波器, 數(shù)值實驗?zāi)Ｐ腿鐖D2所示. 設(shè)置模型品質(zhì)因子Q為20, 編程驗證(19)式的有效性和數(shù)值模擬方案的可行性. 圖3為600 ms的波場, 由于介質(zhì)中不含波阻抗界面, 因此波場特征是外擴的圓.在模型垂直方向1.9 km處, 水平方向0—3.8 km范圍分別提取500和600 ms波場, 如圖4所示,隨著時間推移, 波場的能量由1.9 km處向外傳遞,受到介質(zhì)的吸收衰減, 波場的振幅明顯變小.

圖2 數(shù)值模擬模型Fig. 2. Numerical simulation model.

圖3 600 ms波場快照Fig. 3. Wavefield snapshot at 600 ms.

圖4 不同時刻的1維波場Fig. 4. One dimensional wavefield at different times.

介質(zhì)的黏彈性不僅衰減波場的振幅, 而且影響相位和頻帶. 現(xiàn)采用不同Q值實施波場模擬, 檢波器從震源激發(fā)時開始采集信息, 直至1000 ms結(jié)束, 在模擬區(qū)域外側(cè)添加吸收衰減層以避免邊界產(chǎn)生虛假反射. 圖5是在相同傳播路徑下, 不同Q值所對應(yīng)的記錄, 由于檢波器接收的信息集中在380 ms, 排除無效信息僅對200至600 ms記錄顯示,Q值由30變?yōu)?0后, 波場的振幅明顯減小,相位發(fā)生畸變. 圖5中兩記錄的歸一化頻譜如圖6所示,Q值影響了記錄頻帶的特征, 主頻由初始的20 Hz分別降至15和10 Hz. 在Kjartansson提出的常Q模型中,Q的減小對應(yīng)著應(yīng)力和應(yīng)變之間的滯后變大, 進而介質(zhì)對波場的吸收衰減變強, 此外該模型的假設(shè)是衰減系數(shù)與波場的頻率成正比關(guān)系, 即低頻成分的衰減比高頻成分弱, 通過數(shù)值實驗表明, 建立的波動方程能夠準確地描述波場的衰減規(guī)律, 且數(shù)值模擬方案可行.

圖5 波場記錄Fig. 5. Wavefield record.

圖6 記錄頻譜Fig. 6. Record spectrum.

為了對比分析(20)式中各項對波場的主控因素, 采用控制參數(shù)λ和μ的方式對方程實施數(shù)值模擬, 圖7(a)—(d)分別對應(yīng)λ=0,μ=0 ,λ=0,μ=1 ,λ=1,μ=0 和λ=1,μ=1 的 截 取 波 場,圖7(a)與圖7(b)和圖7(c)的差異分別在于振幅和相位, 圖7(a)和圖7(d)在振幅和相位上均有差異. 基于以上研究, (20)式中等號左側(cè)第二項主控相位畸變, 第三項主控振幅衰減. 圖8為4個波場的局部(圖2中A點至B點)單道顯示, 該圖展示了彈性介質(zhì)波場經(jīng)振幅衰減和相位畸變后演變?yōu)轲椥越橘|(zhì)波場的過程.

圖7 組合波場Fig. 7. Combined wavefield.

圖8 部分波場Fig. 8. Part wavefield.

為了測試數(shù)值方法對于復(fù)雜模型的有效性, 進行了波場數(shù)值模擬. 采用模型在水平和垂直方向的網(wǎng)格點數(shù)分別為401和301, 網(wǎng)格間距為10 m, 即水平距離為4 km, 垂直深度為3 km, 模型速度分布在2500至3500 m/s區(qū)間, 以層狀構(gòu)造為主, 如圖9(a)所示, 品質(zhì)因子Q取25. 震源布置在模型中央, 三個不同時刻的波場如圖9(b),(c)所示. 隨著時間遞增, 波場能量逐步外擴, 由于介質(zhì)的黏滯性效應(yīng), 波場振幅出現(xiàn)了衰減特征. 通過測試表明,構(gòu)建的黏彈性介質(zhì)波動方程能刻畫波傳播至波阻抗界面所產(chǎn)生的透射波和反射波.

圖9 復(fù)雜模型波場模擬 (a) 速度模型; (b) 0.3 s波場; (c) 0.4 s波場; (d) 0.5 s波場Fig. 9. Wavefield simulation for the complex model: (a) Velocity; (b) wavefield at 0.3 s; (c) wavefield at 0.4 s; (a) wavefield at 0.5 s.

5 結(jié) 語

本文在常Q理論基礎(chǔ)上, 建立了黏彈性介質(zhì)波動方程, 該方程能解耦振幅衰減和相位畸變效應(yīng), 有助于分析介質(zhì)中的復(fù)雜波場現(xiàn)象, 可為地震資料處理提供理論基礎(chǔ). 需要注意的是方程中拉普拉斯算子的階數(shù)是Q的函數(shù), 偽譜法僅能處理Q場為常數(shù)的情況, 對于Q場較平滑的介質(zhì), 可將Q場的平均應(yīng)用至全局. 為了應(yīng)對Q場具有強非均質(zhì)性的波場模擬, 如何將Q從拉普拉斯算子的階數(shù)中分離值得深入研究.