Laguerre譜方法在地震勘探中的應(yīng)用

閔濤，任菊成，邢星

MIN Tao，REN Jucheng，XING Xing

西安理工大學(xué)，理學(xué)院，陜西 西安 710054

School of Science, Xi’an University of Technology, Xi’an 710054,China

引言

地震勘探[1-2]是利用地下介質(zhì)彈性和密度的差異，通過觀測和分析大地對人工激發(fā)地震波的響應(yīng)，推斷地下巖層的性質(zhì)和形態(tài)的地球物理勘探方法。它是鉆探前勘測石油與天然氣資源的重要手段，在煤田和工程地質(zhì)勘查，區(qū)域地質(zhì)研究和地殼研究等方面得到廣泛的應(yīng)用.地震勘探的深度一般從數(shù)十米到幾千米不等.如果假設(shè)地層是橫向均勻的，且震源力只沿縱向分布，則地震波傳播近似的被描述為如下的一維波動方程定解問題

其中u(x, t)表示質(zhì)點(diǎn)震動的速度，f(x, t)為震源函數(shù)，ρ(x)為介質(zhì)密度，k(x)為彈性系數(shù)，T為地面記錄的最大時間，L為所要研究的最大深度.x軸指向地下方，即x為深度.在地表x=0處，我們可接收到反射長波

地震勘探的任務(wù)就是根據(jù)就是根據(jù)方程（1）和測量條件（2）反演地震參數(shù)ρ(x)和k(x).對此問題已經(jīng)有一些文獻(xiàn)對其進(jìn)行了研究，文獻(xiàn)[3]張大力等人利用正則迭代法對k(x)進(jìn)行了參數(shù)反演，文獻(xiàn)[4]張麗琴利用同倫方法進(jìn)行了波阻抗反演等，它們更多關(guān)注于參數(shù)反演的方法而對正問題的研究關(guān)注較少，然而我們知道，要求反問題必須先解決正問題，對于正問題，不同的解決方法自然影響反演的精度.因?yàn)閮H根據(jù)g(t)便要同時進(jìn)行雙參數(shù)反演是不可能的，因此本文假定介質(zhì)密度ρ(x)和彈性系數(shù)k(x)部分信息已知，且主要研究深度L較大時的情況，比如L=500千米，以米為單位，這時我們可以近似的將問題看作為定義在[0,∞)上的問題.對此用無界區(qū)域上的正交多項(xiàng)式(函數(shù))直接在無界區(qū)域上進(jìn)行逼近求解.關(guān)于這方面的詳細(xì)介紹可見綜述性文獻(xiàn)[5-8].而無界區(qū)域上的譜方法便是一種很好的解決方案，例如Hermite譜方法[9-10]和Laguerre譜方法[11-14]，這是因?yàn)樽V方法具有很高的精度，它有助于我們獲得更加準(zhǔn)確的反演結(jié)果.

本文主要從大深度地震勘探一維波動方程[15-18]出發(fā)，首先采用 Laguerre譜方法對其正問題進(jìn)行研究，給出了求解的離散過程，然后用高斯-牛頓法對其進(jìn)行參數(shù)反演并通過仿真進(jìn)行了數(shù)值模擬.

1 Laguerre譜方法

對于半無界區(qū)域上的問題，這里選擇用Laguerre譜方法，為了獲得差分矩陣，首先利用譜配置法，它是一種基于權(quán)函數(shù)的插值，形式如（3）式：

也就是說pN-1(x)作為f(x)的一種插值，有（4）式成立：

對（1）式，在各節(jié)點(diǎn)xk處求n階導(dǎo)數(shù)可得（5）式:

導(dǎo)數(shù)的求解用一個矩陣(n)D 表示，則可推導(dǎo)出（6）式:

故而數(shù)值差分過程如下：

其中，f為函數(shù)在各節(jié)點(diǎn)xk處的函數(shù)值，f(n)為函數(shù)在各節(jié)點(diǎn)xk處近似得到的導(dǎo)數(shù)值。求解微分方程時，導(dǎo)數(shù)通過（7）的離散求導(dǎo)去近似.對于 Laguerre譜方法，我們?nèi)?x1=0,x2…xN為LN-1(x)的根，其中LN-1(x)為Laguerre 多項(xiàng)式（關(guān)于權(quán)函數(shù)α(x)=e-x2正交），N-1為Laguerre多項(xiàng)式的次數(shù):

權(quán)函數(shù)為α(x)=e-x/2，則插值函數(shù)為：

這樣我們便獲得了n階的Laguerre譜差分矩陣.

其中k'(x)為k(x)的導(dǎo)數(shù)（這里只考慮k(x)可導(dǎo)的情形），對空間變量進(jìn)行Laguerre離散得m階導(dǎo)數(shù)：

其中，D(1),D(2)分別為一階，二階 Laguerre譜差分矩陣，這樣方程（1）變成了（11）和（12）組成的常微分方程組:

2 仿真實(shí)驗(yàn)

圖1 u(x, t)的近似解Fig 1 The approximate solution of u(x, t)

圖2 u(x, t)的真解Fig 2 The true solution of u(x, t)

表1 不同時間點(diǎn)處真解與近似解的誤差Table 1 The error of true solution and approximate solution at different time points

0.6 9.544235569643507e-015 1.555719199907700e-014 0.7 1.089865647698605e-014 1.305177103030653e-014 0.8 1.186760807613160e-014 1.088117471906690e-014 0.9 1.308203817152395e-014 9.477263496488302e-015 1 1.418303078577562e-014 8.322649553556040e-015 CPU時間（秒） 0.119568

從以上我們可以看出，此方法對于求解此類問題具有非常高的精度且計算速度較快.

有了高效的正問題求解，下面將由方程（1）和測量條件（2）反演地震參數(shù)ρ(x)和k(x).顯然僅根據(jù)g(x)便要同時進(jìn)行雙參數(shù)反演是不可能的.因此本文假定密度為 ρ(x)=x+a1，彈性系數(shù)為k(x)=x2+a2x ,其中a1, a2是需要反演的參數(shù)，為了驗(yàn)證方法的有效性，首先給出參數(shù)真值a1=1,a2=0.5，可通過求解正問題得 u(0,t)=g(t)，并把它作為附加條件來反求參數(shù) a1，a2,其中g(shù)(t)=t2.設(shè) u(0,ti)是在地面接收到反射長波,u(0,ti, a1, a2)是以a1, a2為參數(shù)解正問題所得的計算結(jié)果，則參數(shù)反演問題轉(zhuǎn)化為如下非線性優(yōu)化問題

采用高斯-牛頓法求解，具體步驟如下：

第三：為了避免求解的不穩(wěn)定性，利用正則化方法，將方程組轉(zhuǎn)化為(ATA+αI)·Δ=ATG來求解σi，進(jìn)而得到當(dāng)σi值較大時，可令當(dāng)前的ai值代替原來的近似值，重復(fù)上述過程，得到新的σi（進(jìn)而得ai）.這種過程可以反復(fù)迭代，直到指定的迭代次數(shù)為止.

由于在實(shí)際問題中，接收到的反射波會受到各種因素的影響，所以對條件（2）施加干擾，如=(1+δrand(1))g(t)，我們?nèi)〕鯐r猜測(a1, a2)=(0.1,0.1)，參數(shù)α=0.03,σ=0.0001，迭代100次，δ取不同值時計算結(jié)果如下（真解為：(a1, a2)=(1,0.5)）：

表2 不同程度干擾反演結(jié)果Table 2 The inversion results of different degree interference

計算結(jié)果表明：（1）當(dāng)?shù)孛嬗涗沢(t)含一定噪聲時，反演結(jié)果與真實(shí)速度的吻合程度非常好，其相對誤差幾乎可以忽略.從迭代過程來看，雖然要反演的有2個數(shù)，但用該方法只需要很短時間就可以找到全局極小點(diǎn)，當(dāng)然如果迭代次數(shù)更多的話，反演結(jié)果會更接近真解.（2）雖然隨著噪聲水平的增加，反演結(jié)果的精度有所降低，但由上面的計算結(jié)果可以看出，反演結(jié)果仍是令人滿意的，這說明該方法具有一定的抗噪能力.

3 結(jié)論

針對大深度地震勘探問題，本文首先采用Laguerre譜方法對其正問題進(jìn)行研究，給出了求解的離散，然后用高斯-牛頓法對其進(jìn)行參數(shù)反演并通過仿真進(jìn)行了數(shù)值模擬.結(jié)果表明這種方法對于大深度地震勘探具有較好的效果.

[1]戴向峰 .姜太亮 鄭曉英 星全玲.地震勘探方法及應(yīng)用-地震屬性分析[J].青海石油.2011，29（3）：13-17

[2]馬在田等編著，計算地球物理學(xué)概論[M].同濟(jì)大學(xué)出版社，1997

[3]Huang G Y.An algorithm for processing seismic exploration data .Chinese J.Geophys(in Chinese), 1985，28(1):74-83

[4]張大力，吳建成，劉家琦.一維波動方程反問題求解的正則迭代法 [J].計算物理, 2000,17（3）:326-330.

[5]張麗琴，王家映，嚴(yán)德天.一維波動方程波阻抗反演的同倫方法 [J].地球物理學(xué)報, 2004,47（6）:1111-1117.

[6]Jie Shen and Lilian Wang.Some recent advances on spectral methods for unbounded domain [J].Commun.Comput.Phys,2009,5(2-4):195-241

[7]Jie Shen.Stable and efficient spectral methods in unbounded do mains using laguerre functions[ J].SIAM J.NUMER.ANAL, 2000, 38(4) :1113～ 1133.

[8]B.Y.Guo and J.Shen.Laguerre-Galerkin method for nonlinear partial differential equations on a seminfinite interval[ J].Numerische Mathematik, 2000,86(4) :635～ 654.

[9]V.Iranzo and A.Falques.Some spectral approximations for differential equations in unbounded domains[ J].Comp.Meth.in Appl .Mech and Eng, 1992, 98:105-126.

[10]王中慶.無界區(qū)域問題的有理譜方法[ D].上海:上海大學(xué),2002.

[11]張璟，無窮域問題的譜方法研究[D]上海大學(xué):上海大學(xué)博士學(xué)位論文.2003

[12]王沖，無界區(qū)域問題的Laguerre譜方法 [J].山東理工大學(xué)學(xué)報（自然科學(xué)版）2010,24（6）:85-87.

[13]葉小華.四階方程的Legendre-Laguerre復(fù)合譜方法[J].吉林師范大學(xué)學(xué)報（自然科學(xué)版）.2009，5（2）：121-128

[14]徐承龍，郭本瑜.多維區(qū)域中非線性偏微分方程的修正 Laguerre譜與擬譜方法[J].應(yīng)用數(shù)學(xué)和力學(xué).2008，29（3）：281-300

[15]Han B, Yang X J, Liu J Q.Differential continuation-regularization method and the coefficient inverse problem of one-dimensional wave equations.Applied Mathematics-A Journal of Chinese Universities(in Chinese) , 1994, 9(A4) :351～ 360

[16]劉家琦，劉克安等.微分方程反演波阻抗剖面[J].地球物理學(xué)報.1994，37（1）：101-107

[17]Xie GQ，A new iterative method for solving the coefficient inverse problem of the wave equation[J].Comm pure Appl Math，1986，XXXIX,307-332

[18]J.P.Boyd.Chebyshev and Fourier spectral methods [M].Courier Dover Publication,2001