梁文全,王彥飛,楊長春

(1.龍巖學(xué)院資源工程學(xué)院,福建龍巖364000；2.中國科學(xué)院地質(zhì)與地球物理研究所,中國科學(xué)院油氣資源研究重點實驗室,北京100029)

基于優(yōu)化方法的時間-空間域隱格式有限差分算子確定方法

梁文全1,王彥飛2,楊長春2

(1.龍巖學(xué)院資源工程學(xué)院,福建龍巖364000；2.中國科學(xué)院地質(zhì)與地球物理研究所,中國科學(xué)院油氣資源研究重點實驗室,北京100029)

聲波方程數(shù)值模擬構(gòu)成了地震逆時偏移成像技術(shù)和全波形反演的基礎(chǔ)。對于有限差分法而言,在滿足一定的穩(wěn)定性條件時,普遍存在著因網(wǎng)格化而形成的數(shù)值頻散效應(yīng)。如何有效地壓制數(shù)值頻散是有限差分方法研究的關(guān)鍵所在。為了進(jìn)一步抑制數(shù)值頻散,利用隱式有限差分比顯式有限差分更能壓制數(shù)值頻散的特點,采用前人提出的新的有限差分模板(在保持相同精度的情況下增大了時間步長),應(yīng)用信賴域優(yōu)化方法在時間-空間域確定隱格式有限差分系數(shù)。頻散分析和數(shù)值模擬試算的結(jié)果表明,這種新模板隱格式有限差分優(yōu)化方法既提高了聲波數(shù)值模擬精度又提高了計算效率。

聲波數(shù)值模擬；數(shù)值頻散；隱格式有限差分；時間-空間域；信賴域優(yōu)化方法

地震波場的數(shù)值模擬技術(shù)是在已知地下介質(zhì)結(jié)構(gòu)和參數(shù)的情況下,利用波動方程研究地震波在地下介質(zhì)中的傳播規(guī)律并合成地震記錄的一種技術(shù)。有限差分方法因為計算效率高、所需內(nèi)存小、實現(xiàn)簡單而廣泛應(yīng)用于地震正演研究[1-6],同時構(gòu)成了逆時偏移成像和全波形反演的基礎(chǔ)[7-11]。

有限差分?jǐn)?shù)值模擬的關(guān)鍵問題是數(shù)值頻散和吸收邊界。數(shù)值頻散是對波動方程的時間和空間偏導(dǎo)數(shù)的離散化所造成的,使得相速度變成了網(wǎng)格間距的函數(shù),這會導(dǎo)致地震波的數(shù)值相速度不等于地球介質(zhì)的真實相速度,使得波場模擬的精度降低。一般來說,如果存在時間頻散,則高頻的相速度增大；如果存在空間頻散,則高頻的相速度減小[1]。邊界條件大致有兩種處理方法,一種是衰減邊界條件方法(如完全匹配層等),另一種是基于單程波近似理論的吸收邊界條件(ABC)方法。

偏微分方程的有限差分格式分為顯式有限差分和隱式有限差分。顯式有限差分格式中,某點導(dǎo)數(shù)的計算需要該點和它周圍點的值。隱式有限差分格式中,某點導(dǎo)數(shù)的計算除了需要該點和它周圍點的值之外,還需要它周圍點的導(dǎo)數(shù)值[5]。眾所周知,隱式有限差分格式中空間偏導(dǎo)數(shù)的確定需要求解三對角陣。傳統(tǒng)的算法中,三對角陣的求逆需要較大的計算量。但是,若要達(dá)到相同的精度,使用隱式有限差分算子比使用顯式有限差分算子一般需要更少的網(wǎng)格點。

起初,波動方程空間偏導(dǎo)數(shù)的有限差分算子在空間域確定。但是,地震波在時間域和空間域同時傳播,因此在時間-空間域確定空間偏導(dǎo)數(shù)的有限差分算子能同時抑制時間頻散和空間頻散。Liu等[12]通過平面波理論和泰勒展開方法在時間-空間域確定有限差分算子系數(shù)。梁文全等[13]提出根據(jù)震源頻率范圍、波速和網(wǎng)格間距確定一個合適的波數(shù)范圍,在此波數(shù)范圍內(nèi)使用線性方法確定有限差分系數(shù),并取得良好效果。為了進(jìn)一步壓制頻散和增加時間步長,Liu等[14]提出了使用菱形有限差分模板并確定了對應(yīng)的有限差分系數(shù)。Liu等[4]也提出了新的有限差分模板,在保持相同精度的情況下增大了時間步長,并用最小二乘方法確定了相應(yīng)的有限差分系數(shù)。本文應(yīng)用信賴域優(yōu)化方法在時間-空間域確定隱格式有限差分新模板對應(yīng)的系數(shù),以期在保證有限差分算子精度的同時提高其計算效率。

1 方法原理

二維的聲波方程可以寫為：

(1)

式中：p是波場；v是波速。時間二階導(dǎo)數(shù)的顯式有限差分算子表達(dá)式為：

(2)

二階空間偏導(dǎo)數(shù)的隱式有限差分(與以前格式略有不同,在每個方向上多了4個點)[4-5,14-16]為：

(3)

(4)

其中，

(5)

(6)

式中：cm,b和d是需要確定的交錯網(wǎng)格有限差分系數(shù)(當(dāng)d=0時為傳統(tǒng)有限差分模板,如圖1a所示；當(dāng)d≠0時為新的有限差分模板,如圖1b所示)。將(2)式、(3)式和(4)式代入(1)式,得到：

(7)

整理公式(7)得到：

(8)

對公式(8)進(jìn)行傅里葉變換,得到：

圖1 傳統(tǒng)有限差分模板(a)和新的有限差分模板(b)

[cos(khcosθ)-1]}[2cos(?τ)-2]

(9)

式中：r=vτ/h。為了求得有限差分系數(shù),建立目標(biāo)函數(shù)使得速度誤差最小,即：

(10)

式中：q表示需要優(yōu)化的范圍。其中,

(11)

(12)

g2=2bcos(khcosθ)-2+2bcos(khsinθ)-2+

(13)

為求解φ(c)→min,我們使用信賴域方法可以得到優(yōu)化的解。信賴域方法指的是[17-18]：對函數(shù)φ(c)做二次逼近,即通過第k次迭代求解(14)式所示的球約束問題,獲得一個試探步長ξ,并根據(jù)公式ck+1=ck+ξk來更新優(yōu)化系數(shù)。

(14)

subjectto‖ξ‖l2≤Δk

(15)

并且問題的解的表達(dá)式可以寫成：

(16)

2 頻散分析

有限差分方法在一個空間網(wǎng)格上傳播的時間誤差ε為[2]：

(17)

式中：δ為頻散速度vFD與真實速度v的比值,即：

(18)

圖2比較了傳統(tǒng)模板的時間-空間域線性方法和新模板的時間-空間域隱格式優(yōu)化方法的頻散誤差,其中算子長度M=7,空間間距h=20m,時間步長τ=2ms。由圖2可見,時間-空間域新模板隱格式優(yōu)化方法的頻散誤差為傳統(tǒng)模板線性方法頻散誤差的1/10(注意兩圖的縱坐標(biāo)數(shù)值相差一個數(shù)量級)?？梢姴捎眯履０咫[格式優(yōu)化方法可以抑制數(shù)值頻散,提高地震波數(shù)值模擬的精度。

圖2 傳統(tǒng)模板線性方法(a,c,e)和新模板隱格式優(yōu)化方法(b,d,f)的頻散誤差分析(M=7,h=20m,τ=2ms；其中,a,b的波速為1500m/s,c,d的波速為2000m/s,e,f的波速為2500m/s)

3 數(shù)值模擬試算

3.1 均勻介質(zhì)模型

首先考慮均勻介質(zhì)模型,對于所有的有限差分算子,v=2500m/s,τ=2ms。顯式有限差分算子M=10,隱式有限差分算子M=7。計算區(qū)域劃分為530×530網(wǎng)格,網(wǎng)格間距h=10m。震源在模型的中心位置,震源是高斯函數(shù)的二階導(dǎo)數(shù):

(19)

式中：f0=130Hz,f0/π是主頻頻率；t0=4/f0；c是常數(shù)。

圖3顯示了1000ms時不同有限差分算法對應(yīng)的波場快照。由于使用了稀疏存儲的方法和較小的算子長度,隱式有限差分方法的數(shù)值模擬耗時(141s)比顯式有限差分?jǐn)?shù)值模擬耗時(145s)略少,提高了計算效率。由圖3可見,使用傳統(tǒng)模板線性方法的時間頻散和空間頻散較大(圖3a),而使用新模板隱格式優(yōu)化方法同時減小了時間頻散和空間頻散(圖3b)。

圖3 均勻介質(zhì)模型時間-空間域傳統(tǒng)模板線性方法(a)和時間-空間域新模板隱格式優(yōu)化方法(b)的波場快照

3.2 鹽丘模型

圖4顯示了美國勘探地球物理學(xué)會的BP鹽丘模型,其速度從1500m/s變化到4800m/s。震源位置用紅色菱形表示。震源函數(shù)采用公式(12),其中震源頻率f0=45Hz。對于所有的有限差分算子,空間步長h=20m。圖5顯示了3種不同有限差分算子在3000ms時的波形快照(z=400m),其中時間-空間域泰勒展開方法[5](藍(lán)色曲線)的時間步長τ=1ms；時間-空間域傳統(tǒng)模板線性方法(黑色曲線)的時間步長τ=1ms；時間-空間域新模板隱格式優(yōu)化方法(紅色曲線)的時間步長τ=2ms。從圖5可以看到(圖中紅圈所示位置),新模板隱格式優(yōu)化方法在時間步長τ=2ms時的頻散明顯小于時間-空間域泰勒展開方法在時間步長τ=1ms時的數(shù)值頻散,也小于時間-空間域傳統(tǒng)模板線性方法在時間步長τ=1ms時的數(shù)值頻散。由此可見,新模板隱格式有限差分優(yōu)化方法既保證了精度又提高了計算效率。

圖4 BP鹽丘速度模型

圖5 鹽丘模型3種不同有限差分算子在3000ms時的波形快照對比

4 結(jié)論

本文使用信賴域優(yōu)化方法確定了新模板隱格式有限差分系數(shù),在其它條件相同的情況下,時間步長可以從1ms增加到2ms,同時采用稀疏存儲的方法,使得隱格式有限差分的計算效率得到提高。數(shù)值模擬試算結(jié)果表明,采用新模板隱格式有限差分優(yōu)化方法既保證了模擬精度又提高了計算效率。

鑒于新模板隱格式有限差分優(yōu)化方法能夠有效抑制時間頻散和空間頻散,在保持相同精度的情況下增加了時間步長,因此可以代替?zhèn)鹘y(tǒng)的有限差分算子應(yīng)用于聲波數(shù)值模擬中。

[1] Dablain M A.The application of high-order differencing to the scalar wave equation[J].Geophysics,1986,51(1):54-66

[2] Liu Y.Globally optimal finite-difference schemes based on least squares[J].Geophysics,2013,78(4):T113-T132

[3] Liu Y.Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling[J].Geophysical Journal International,2014,197(2):1033-1047

[4] Liu H,Dai N,Niu F,et al.An explicit time evolution method for acoustic wave propagation[J].Geophysics,2014,79(3):T117-T124

[5] Liu Y,Sen M K.A practical implicit finite difference method:examples from seismic modeling[J].Journal of Geophysics Engineering,2009,6(3):231-249

[6] 趙建國,史瑞其,陳競一,等.輔助微分方程完全匹配層在聲波方程數(shù)值模擬中的應(yīng)用[J].吉林大學(xué)學(xué)報(地球科學(xué)版),2014,44(2):675-682 Zhao J G,Shi R Q,Chen J Y,et al.Application of auxiliary differential equation perfectly matched layer for numerical modeling of acoustic wave equations[J].Journal of Jilin University(Earth Science Edition),2014,44(2):675-682

[7] Du X,Bancroft J C,Lines L R.Anisotropic reverse-time migration for tilted TI media[J].Geophysical Prospecting,2007,55(6):853-869

[8] Liu H W,Li B,Liu H,et al.The issues of prestack reverse time migration and solutions with Graphic Processing Unit implementation[J].Geophysical Prospecting,2012,60:906-918

[9] Xie W,Yang D,Liu F,et al.Reverse-time migration in acoustic VTI media using a high-order stereo operator[J].Geophysics,2014,79(3):WA3-WA11

[10] Yan H,Liu Y.Acoustic VTI modeling and pre-stack reverse-time migration based on the time-space domain staggered-grid finite-difference method[J].Journal of Applied Geophysics,2013,90(1):41-52

[11] Yan H,Liu Y.Visco-acoustic prestack reverse-time migration based on the time-space domain adaptive high-order finite-difference method[J].Geophysical Prospecting,2013,61(5):941-954

[12] Liu Y,Sen M K.Scalar wave equation modeling with time-space domain dispersion-relation-based staggered-grid finite-difference schemes[J].Bulletin of the Seismological Society of America,2011,101(1):141-159

[13] 梁文全,楊長春,王彥飛,等.用于聲波方程數(shù)值模擬的時間-空間域有限差分系數(shù)確定新方法[J].地球物理學(xué)報,2013,56(10):3497-3506 Liang W Q,Yang C C,Wang Y F,et al.Acoustic wave equation modeling with new time-space domain finite difference operators[J].Chinese Journal of Geophysics,2013,56(6):840-850

[14] Liu Y,Sen M K.Time-space domain dispersion-relation-based finite-difference method with arbitrary even-order accuracy for the 2D acoustic wave equation[J].Journal of Computational Physics,2013,232(1):327-345

[15] Tan S,Huang L.An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation[J].Geophysical Journal International,2014,197(2):1250-1267

[16] Tan S,Huang L.A staggered-grid finite-difference scheme optimized in the time-space domain for modeling scalar-wave propagation in geophysical problems[J].Journal of Computational Physics,2014,276:613-634

[17] Wang Y F,Yuan Y X.Convergence and regularity of trust region methods for nonlinear ill-posed inverse problems[J].Inverse Problems,2005,21(3):821-838

[18] 王彥飛.反演問題的計算方法及其應(yīng)用[M].北京：高等教育出版社,2007:1-370 Wang Y F.Computational methods for Inverse problems and applications[M].Beijing:Higher Education Press,2007:1-370

(編輯：顧石慶)

Determination on the implicit finite-difference operator based on optimization method in time-space domain

Liang Wenquan1,Wang Yanfei2,Yang Changchun2

(1.CollegeofResourceEngineering,LongyanUniversity,Longyan364000,China；2.KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China)

Numerical simulation of acoustic wave equation is the basis of the reverse time migration and full waveform inversion. With some stability conditions,grid dispersion often exists because of the discretization of the time and the spatial derivatives in the wave equation. How to suppress the grid dispersion is therefore a key problem for finite difference approaches. Finite-difference scheme of partial differential equations are either explicit or implicit. Implicit finite difference scheme can suppress the grid dispersion further. Meanwhile,a new finite difference stencil was proposed to increase the time step while preserving the accuracy by previous researchers. Therefore,we propose to determine the optimal implicit finite difference operators in time-space domain with a new finite difference stencil. Dispersion analysis and seismic numerical simulation demonstrate that the accuracy of acoustic wave numerical simulation and the computation efficiency can be improved by the method.

acoustic wave numerical simulation,grid dispersion,implicit finite difference scheme,time-space domain,trust region optimization method

2014-10-17;改回日期：2014-12-29。

梁文全(1983—)，男，講師，現(xiàn)主要研究方向為地震波數(shù)值模擬及偏移成像。

國家自然科學(xué)基金(41325016,11271349)、國家重大科研裝備研制項目(ZDYZ2012-1-02-04)、龍巖學(xué)院博士科研啟動基金(LB2014010)和福建省教育廳JK項目(2014050)聯(lián)合資助。

P631

A

1000-1441(2015)03-0254-06

10.3969/j.issn.1000-1441.2015.03.002