閔濤，任菊成，邢星

MIN Tao，REN Jucheng，XING Xing

西安理工大學，理學院，陜西 西安 710054

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

引言

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

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

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

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

1 Laguerre譜方法

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

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

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

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

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

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

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

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

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

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

2 仿真實驗

圖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 不同時間點處真解與近似解的誤差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)便要同時進行雙參數(shù)反演是不可能的.因此本文假定密度為 ρ(x)=x+a1，彈性系數(shù)為k(x)=x2+a2x ,其中a1, a2是需要反演的參數(shù)，為了驗證方法的有效性，首先給出參數(shù)真值a1=1,a2=0.5，可通過求解正問題得 u(0,t)=g(t)，并把它作為附加條件來反求參數(shù) a1，a2,其中g(t)=t2.設 u(0,ti)是在地面接收到反射長波,u(0,ti, a1, a2)是以a1, a2為參數(shù)解正問題所得的計算結果，則參數(shù)反演問題轉化為如下非線性優(yōu)化問題

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

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

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

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

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

3 結論

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

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

[2]馬在田等編著，計算地球物理學概論[M].同濟大學出版社，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]張麗琴，王家映，嚴德天.一維波動方程波阻抗反演的同倫方法 [J].地球物理學報, 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].上海:上海大學,2002.

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

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

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

[14]徐承龍，郭本瑜.多維區(qū)域中非線性偏微分方程的修正 Laguerre譜與擬譜方法[J].應用數(shù)學和力學.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].地球物理學報.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