王　穎,周　輝,盛善波

(1.中國石油大學(北京)油氣資源與探測國家重點實驗室,中國石油大學(北京)CNPC物探重點實驗室,中國石油大學(北京),北京102249；2.中國石油天然氣勘探開發(fā)公司,北京100034)

?

貼體坐標系二維聲波方程SBP有限差分法的穩(wěn)定性分析

王穎1,周輝1,盛善波2

(1.中國石油大學(北京)油氣資源與探測國家重點實驗室,中國石油大學(北京)CNPC物探重點實驗室,中國石油大學(北京),北京102249；2.中國石油天然氣勘探開發(fā)公司,北京100034)

摘要:有限差分方法(Finite Difference Method,FDM)是波動方程正演數(shù)值模擬領(lǐng)域應用最為廣泛的方法之一,然而,當模擬區(qū)域不規(guī)則或者地表起伏不平時,規(guī)則網(wǎng)格有限差分法求解波動方程會產(chǎn)生階梯狀近似,影響模擬的精度。借助貼體網(wǎng)格技術(shù),將不規(guī)則的物理區(qū)域轉(zhuǎn)換為規(guī)則的計算域,給出了貼體坐標系下的二維聲波方程及其二階精度的分部求和(Summation by Parts,SBP)有限差分離散格式，采用Fourier譜分析方法分析了該離散格式的穩(wěn)定性,得到了貼體網(wǎng)格二維聲波方程SBP有限差分方法的穩(wěn)定性條件。數(shù)值實驗結(jié)果表明：①當時間采樣間隔的選取滿足穩(wěn)定性條件時,貼體網(wǎng)格SBP有限差分的數(shù)值計算過程是穩(wěn)定的；②與貼體網(wǎng)格中心差分方法相比，貼體網(wǎng)格SBP有限差分方法的穩(wěn)定性更好。

關(guān)鍵詞:貼體坐標；聲波方程；分部求和有限差分；中心差分；穩(wěn)定性條件

自從1968年ALTERMAN等[1]首次將有限差分方法(FDM)應用于層狀介質(zhì)的彈性波動方程求解問題以來,FDM已經(jīng)成功地發(fā)展成為應用最為廣泛的地震波場數(shù)值模擬方法之一,并在復雜不均勻介質(zhì)波場模擬問題上發(fā)揮著越來越重要的作用。FDM將求解域劃分為規(guī)則差分網(wǎng)格,大多用來處理水平地表的模擬問題。在地表起伏不平或者模擬區(qū)域不規(guī)則的情況下,使用FDM需對求解區(qū)域不規(guī)則的邊界進行階梯狀近似,會在階梯狀網(wǎng)格的角點處產(chǎn)生虛假的散射和繞射,從而影響波場模擬的精度[2-3]。

有限元、邊界元等方法適用于不規(guī)則區(qū)域的地震波場數(shù)值模擬。其中有限元法處理不規(guī)則地表問題比較方便,能有效地模擬起伏地表情況下的地震波傳播,但模擬效率低,精度不高,對機器內(nèi)存需求也較大。與有限差分方法相比,當二者精度相同時,有限元法更加費時[4-7]。邊界積分或邊界元法適用于特定形狀下的地震波傳播,但不適用于地表速度變化較大的情況[8]。相比而言,有限差分法計算效率較高,易于編程和并行運算,并且能夠很好地處理各種強不均勻介質(zhì)問題,因此仍然是解決地震波場數(shù)值模擬問題最有效的方法。

貼體坐標技術(shù)的提出為有限差分法處理不規(guī)則區(qū)域問題開辟了新的局面[9-13]。20世紀80年代初期,貼體網(wǎng)格開始用于求解不規(guī)則區(qū)域的物理問題,并逐漸成為最主要的數(shù)值方法之一。TESSMER等[2]基于坐標變換的方法將起伏地表下的曲線網(wǎng)格映射到水平地表下的矩形網(wǎng)格,并使用偽譜法實現(xiàn)了不規(guī)則區(qū)域中二維彈性波場的數(shù)值模擬。HESTHOLM等[4]基于縱向坐標變換,使用高階交錯網(wǎng)格有限差分算子實現(xiàn)了二維起伏地表下的彈性波場數(shù)值模擬。1998年,HESTHOLM等[14]將此方法擴展到三維情形,給出了所能模擬的起伏地表最大起伏程度。2003年,HESTHOLM[15]針對不同的dz/dx(z方向和x方向空間采樣間隔比值)和縱橫波速度比值進行了基于映射的彈性波場數(shù)值模擬,并通過實驗給出了達到長時間穩(wěn)定的dz/dx取值范圍,但未進行理論上的推導,也未討論時間采樣間隔的影響。董良國[16]在橫波品質(zhì)因子Q值的基礎上,借助縱向坐標變換,利用交錯網(wǎng)格高階差分方法對起伏地表情況下的彈性波傳播進行了數(shù)值模擬。針對規(guī)則網(wǎng)格中二階形式的彈性波方程和自由邊界條件,NILSSON等[17]給出了二階精度的分部求和(SBP)差分格式,討論了非均勻介質(zhì)中該差分算法的穩(wěn)定性條件。APPEL?等[18]將NILSSON等[17]研究的方法擴展到曲線坐標系,通過構(gòu)造曲線網(wǎng)格中的二階SBP離散格式,實現(xiàn)了不規(guī)則區(qū)域的彈性波模擬,但未涉及曲線坐標系中相應的穩(wěn)定性條件。TARRASS等[19]借鑒空氣聲學中的非中心差分算子,使用高階的空間和時間優(yōu)化算子模擬起伏地表下的波場傳播,給出了差分格式的穩(wěn)定性條件,但是該方法的網(wǎng)格剖分基于坐標映射,導致該方法的穩(wěn)定性在一定程度上依賴于地形的起伏,當?shù)匦屋^陡時,該方法穩(wěn)定性變差。LAN等[8]將APPEL?等[18]研究的方法擴展到三維非均質(zhì)各向異性介質(zhì),采用曲線坐標離散并求解垂直對稱軸橫向各向同性(VTI)介質(zhì)中的彈性波方程和自由邊界條件,實現(xiàn)了起伏地表下VTI介質(zhì)中地震波場的顯式有限差分模擬。此后,LAN等[20]采用類似的方法實現(xiàn)了起伏地表下水平對稱軸橫向各向同性(HTI)介質(zhì)中地震波場的顯式有限差分模擬。劉一峰等[21]借助貼體網(wǎng)格技術(shù),研究了曲線坐標系中程函方程的求解方法。唐文等[22]通過坐標變換,實現(xiàn)了起伏地表二階彈性波方程的地震波場數(shù)值模擬,研究了起伏地表下的自由邊界條件,定性地比較了采用不同自由邊界條件時模擬算法的穩(wěn)定性,但未定量給出差分格式的穩(wěn)定性條件。

本文基于貼體坐標的概念,首先給出貼體網(wǎng)格中的二維聲波方程以及完全匹配層(PML)吸收邊界條件；然后采用二階精度的SBP有限差分對貼體網(wǎng)格的聲波方程進行離散,推導了該離散格式的穩(wěn)定性條件；最后通過數(shù)值模擬驗證該穩(wěn)定性條件的正確性,并通過與貼體網(wǎng)格中心差分方法對比,顯示貼體網(wǎng)格SBP有限差分方法在穩(wěn)定性方面的優(yōu)勢。

1技術(shù)方法

1.1貼體坐標系中二維聲波方程、PML及其離散格式

1.1.1貼體坐標系的生成

貼體坐標系是由物理區(qū)域的邊界點坐標經(jīng)過數(shù)值計算生成的坐標系,能夠很好地與邊界形狀相適應。貼體網(wǎng)格的生成可以看作是一種坐標變換,它把物理平面上的不規(guī)則區(qū)域變換為計算平面上的規(guī)則區(qū)域。一般采用偏微分方程法實現(xiàn)貼體坐標變換,變換關(guān)系為[9-10]:

(1)

(2)

其中,

(3)

1.1.2貼體坐標系中的聲波方程、PML和中心有限差分算子

經(jīng)過坐標變換,貼體坐標系中的二維聲波方程可以表示為:

(4)

其中,u為聲壓場,v為聲波的傳播速度,t為時間變量。

為了消除正演模擬過程中的人工邊界反射,在貼體坐標系中引入了PML邊界條件,經(jīng)過一系列的推導,得到方程(4)對應的時間域分裂式PML匹配方程:

(5)

其中,u=u1+u2+u3+u4+u5,u1,u2,…,u5為分裂的聲壓場；dq,dr代表PML區(qū)域中q,r方向的衰減函數(shù)；“*”表示時間褶積。(4)式和(5)式中的時間導數(shù)和空間導數(shù)均采用二階中心差分,即:

(6)

將(6)式代入方程(4)和方程(5)中即可得到貼體坐標系中二維聲波方程和PML匹配方程的二階中心差分離散格式。

1.1.3貼體坐標系中的SBP有限差分算子

地下介質(zhì)的強非均質(zhì)性使得介質(zhì)參數(shù)在遠小于地震波波長的空間范圍內(nèi)變化巨大,這種變化在計算區(qū)域的網(wǎng)格中更加劇烈[17]。為了確保計算的穩(wěn)定性,必須研究合適的數(shù)值方法。1974年,KREISS等[24]提出了一種穩(wěn)定的有限差分方法,可以求解連續(xù)時間系統(tǒng)中能量守恒的偏微分方程。由于該方法滿足SBP性質(zhì),因而被稱為SBP有限差分方法。

為了方便地使用SBP有限差分方法,我們將貼體坐標系中的二維聲波方程寫為:

(7)

其中,

(8)

可以證明,方程(7)和方程(4)是等價的。該方程對應的時間域分裂式PML匹配方程可以表示為:

(9)

其中,u=u1+u2+u3。

將(7)式中的空間變量導數(shù)項替換為SBP有限差分,時間變量導數(shù)項替換為二階中心差分,可以得到貼體坐標系中二維聲波方程的SBP有限差分離散格式:

(10)

其中,D+,D-為標準的向前、向后差分算子,以r方向為例,

(11)

(12)

類似地,可以得到PML匹配方程(9)的SBP有限差分離散格式:

(13a)

(13b)

(13c)

其中,

(14)

通過能量估計可以證明,該SBP差分格式滿足離散能量守恒[17],與傳統(tǒng)的中心差分方法相比,具有更好的穩(wěn)定性,形式也更簡潔[25]。

1.2穩(wěn)定性分析

本節(jié)討論貼體坐標系中二維聲波方程的SBP有限差分格式(10)的穩(wěn)定性條件。常用的有限差分穩(wěn)定性分析方法是VONNEUMANN等[26]和SJ?GREEN[27]提出的Fourier譜分析方法,該方法適用于常系數(shù)初值和帶周期邊值條件的混合問題,因此假定公式(10)中的介質(zhì)參數(shù)Mqq,Mrr和J為常量。由于

(15)

而在曲線坐標系中坐標線正交的條件是β=0,因此Mqr=0,公式(10)簡化為：

(16)

公式(16)中的空間差分算子在Fourier域?qū)獮?

(17)

(18)

經(jīng)過Fourier變換,公式(16)對應為:

(19)

其中,

(20)

(21)

公式(21)右端的矩陣稱為增長矩陣,簡記為G。當增長矩陣的譜半徑ρ(G)≤1,即G的所有特征值的模小于等于1時,差分格式穩(wěn)定。矩陣G的特征值λ滿足的特征方程為:

(22)

特征方程的解為:

(23)

如果0≤δ≤2,判別式的值為負值或0,特征方程有一對復共軛的根:

(24)

根的模滿足

(25)

此時穩(wěn)定的必要條件得到滿足。進一步分析0≤δ≤2時對應的時間采樣間隔:

(26)

(27)

穩(wěn)定性條件為:

(28)

值得注意的是,受Fourier譜分析方法本身的限制,針對貼體網(wǎng)格聲波方程SBP有限差分離散格式的穩(wěn)定性分析較為粗略,并帶有一定的近似。能量法可以精確地分析穩(wěn)定性問題,但是能量法中譜半徑的求解是個難題,即使在笛卡爾坐標系中也無法得到其解析解[27]。因此,Fourier譜分析方法仍不失為貼體網(wǎng)格中聲波方程SBP有限差分格式穩(wěn)定性分析的一種有效手段。

2數(shù)值模擬

為了驗證貼體坐標系中二維聲波方程SBP有限差分格式的穩(wěn)定性條件,即公式(28)的正確性,選取不同的時間采樣間隔進行正演模擬,通過與二階中心差分格式的貼體坐標聲波方程正演模擬結(jié)果進行對比,顯示出貼體網(wǎng)格SBP有限差分方法在穩(wěn)定性方面的優(yōu)勢。鑒于本文主要是為了驗證差分格式的穩(wěn)定性,因此在正演模擬過程中,模型的四周均采用了PML吸收邊界條件。

為了簡單起見,設計帶有正弦函數(shù)起伏地表的二維均勻模型,如圖1所示。模型大小為4000m×3000m,介質(zhì)的速度為2500m/s。該模型的地表可用正弦函數(shù)y=Asin(2πx/T)表示,其中,A=-60m,T=520m,即該模型地表的起伏周期為520m,相對最大高差為120m。震源采用主頻為20Hz的Ricker子波,震源位置坐標為(2080m,36m)。正演模擬采用的空間采樣間隔為Δq=Δr=4m。

圖1　正弦函數(shù)起伏地表模型

首先利用(28)式計算出該模型SBP有限差分格式的穩(wěn)定性條件為Δt≤0.93ms。為了對比驗證該穩(wěn)定性條件,分別選取時間采樣間隔Δt=0.90ms以及Δt=1.00ms進行正演模擬,記錄長度為6s。圖2a和圖2b分別為時間采樣間隔為0.90ms和1.00ms時的模擬地震記錄,可以看出,在滿足穩(wěn)定性條件的情況下,利用貼體網(wǎng)格SBP有限差分方法可以得到清晰的直達波記錄(圖2a),而在時間采樣間隔不滿足穩(wěn)定性條件時,出現(xiàn)數(shù)值結(jié)果的溢出(圖2b)。圖2c和圖2d分別是從圖2a和圖2b兩個地震記錄中抽取的5個地震道,抽取的位置如圖2a和圖2b中倒三角所示。

從抽取的單道記錄中可以得出相同的結(jié)論:當時間采樣間隔為0.90ms時,滿足穩(wěn)定性條件,此時貼體網(wǎng)格SBP有限差分方法模擬的整個記錄是穩(wěn)定的;而當時間采樣間隔為1.00ms時,超出了穩(wěn)定性條件的范圍,此時貼體網(wǎng)格SBP有限差分方法的模擬結(jié)果不穩(wěn)定,造成數(shù)值結(jié)果溢出,波動方程的數(shù)值求解失敗。

圖3a和圖3b分別為貼體網(wǎng)格的SBP有限差分方法和中心有限差分方法正演模擬得到的地震記錄對比,兩種方法的時間采樣間隔均為Δt=0.93ms。正演模擬記錄時間長度為30s,為了美觀,只顯示到23s左右。圖3c和圖3d分別是從圖3a和圖3b兩個地震記錄中抽取的3個地震道,抽取的位置如圖3a和圖3b中倒三角所示。從圖3 可以看出,在相同時間采樣間隔的情況下,貼體網(wǎng)格SBP有限差分方法和中心有限差分方法均可以得到清晰的直達波記錄,但中心有限差分方法在18s左右開始出現(xiàn)不穩(wěn)定的數(shù)值,并且這種不穩(wěn)定性隨著時間的推移越來越嚴重。從抽取的單道記錄中可以得出相同的結(jié)論:貼體網(wǎng)格SBP有限差分方法模擬的整個記錄都是穩(wěn)定的,而中心有限差分方法模擬的記錄在后期出現(xiàn)不穩(wěn)定,數(shù)值結(jié)果溢出。

圖2　時間采樣間隔為0.90ms(a)和1.00ms(b)正演模擬的地震記錄及其5個地震道(c)和(d)

圖3　貼體網(wǎng)格SBP有限差分法(a)和中心差分法(b)正演模擬的地震記錄及其3個地震道(c)和(d)

3結(jié)束語

本文針對貼體坐標系中聲波方程的SBP有限差分離散格式,采用Fourier譜分析方法討論了該離散格式的穩(wěn)定性條件，通過對起伏地表均勻介質(zhì)模型選取不同的時間采樣間隔進行正演模擬對比,驗證了該穩(wěn)定性條件的正確性。與貼體網(wǎng)格中心差分方法相比,貼體網(wǎng)格SBP有限差分方法更加穩(wěn)定,適用于求解不規(guī)則區(qū)域的波動方程。

本文的SBP有限差分算子僅限于二階精度,貼體坐標系中聲波方程的高精度SBP有限差分算法及對應的穩(wěn)定性條件將是我們今后研究的重點。

參考文獻

[1]ALTERMAN Z,KARAL F C.Propagation of elastic waves in layered media by finite difference methods[J].Bulletin of the Seismological Society of America,1968,58(1):367-398

[2]TESSMER E,KOSLOFF D,BEHLE A.Elastic wave propagation simulation in the presence of surface topography[J].Geophysical Journal International,1992,108(2):621-632

[3]ZHANG W,CHEN X F.Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[J].Geophysical Journal International,2006,167(1):337-353

[4]HESTHOLM S,RUUD B.2D finite-difference elastic wave modeling including surface topography[J].Geophysical Prospecting,1994,42(5):371-390

[5]KOMATITSCH D,TROMP J.Introduction to the spectral element method for three-dimensional seismic wave propagation[J].Geophysical Journal International,1999,139(3):806-822

[6]BROSSIER R,VIRIEUX J,OPERTO S.Parsimonious finite-volume frequency-domain method for 2D P-SV-wave modeling[J].Geophysical Journal International,2008,175(2):541-559

[7]LIU T,HU T Y,SEN M K,et al.A hybrid scheme for seismic modelling based on Galerkin method[J].Geophysical Journal International,2011,186(3):1165-1178

[8]LAN H Q,ZHANG Z J.Three-dimensional wave-field simulation in heterogeneous transversely isotropic medium with irregular free surface[J].Bulletin of the Seismological Society of America,2011,101(3):1354-1370

[9]THOMPSON J F,THAMES F C,MASTIN C W.Automatic numerical generation of body-fitted curve linear coordinate system for field containing any number of arbitrary two-dimensional bodies[J].Journal of Computational Physics,1974,15(3):299-319

[10]THOMPSON J F,THAMES F C,MASTIN C W.TOMCAT-A code for numerical generation of boundary-fitted curvilinear coordinate systems on fields containing any number of arbitrary two-dimensional bodies[J].Journal of Computational Physics,1977,24(3):274-302

[11]THOMPSON J F,WARSI Z U A,MASTIN C W.Numerical grid generation:foundations and applications[M].New York:Elsevier North-Holland,1985:1-3,95

[12]LAKNER M,PLAZL I.The finite differences method for solving systems on irregular shapes[J].Computers and Chemical Engineering,2008,32(12):2891-2896

[13]KAUL U K.Three-dimensional elliptic grid generation with fully automatic boundary constraints[J].Journal of Computational Physics,2010,229(17):5966-5979

[14]HESTHOLM S,RUUD B.3-D finite-difference elastic wave modeling including surface topography[J].Geophysics,1998,63(2):613-622

[15]HESTHOLM S.Elastic wave modeling with free surfaces:stability of long simulations[J].Geophysics,2003,68(1):314-321

[16]董良國.復雜地表條件下地震波傳播數(shù)值模擬[J].勘探地球物理進展,2005,28(3):187-194

DONG L G.Numerical simulation of seismic wave propagation under complex near surface conditions[J].Progress in Exploration Geophysics,2005,28(3):187-194

[17]NILSSON S,PETERSSON N A,SJ?GREEN B,et al.Stable difference approximations for the elastic wave equation in second order formulation[J].SIAM Journal on Numerical Analysis,2007,45(5):1902-1936

[18]APPEL? D,PETERSSON N A.A stable finite difference method for the elastic wave equation on complex geometries with free surfaces[J].Communications in Computational Physics,2009,5(1):84-107

[19]TARRASS I,GIRAUD L,THORE P.New curvilinear scheme for elastic wave propagation in presence of curved topography[J].Geophysical Prospecting,2011,59(5):889-906

[20]LAN H Q,ZHANG Z J.Seismic wavefield modeling in media with fluid-filled fractures and surface topography[J].Applied Geophysics,2012,9(3):301-312

[21]劉一峰,蘭海強.曲線坐標系程函方程的求解方法研究[J].地球物理學報,2012,55(6):2014-2026

LIU Y F,LAN H Q.Study on the numerical solutions of the eikonal equation in curvilinear coordinate system[J].Chinese Journal of Geophysics,2012,55(6):2014-2026

[22]唐文,王尚旭,袁三一.起伏地表二階彈性波方程差分策略穩(wěn)定性分析[J].石油物探,2013,52(5):457-463

TANG W,WANG S X,YUAN S Y.Stability analysis of differential strategy for rugged topography by second-order elastic wave equation based on coordinate transformation[J].Geophysical Prospecting for Petroleum,2013,52(5):457-463

[23]WEI W L,JIN Z Q.Numerical solution for unsteady 2-D flow using the transformed shallow water equations[J].Journal of Hydrodynamics,1995,7(3):65-71

[24]KREISS H O,SCHERER G.Finite element and finite difference methods for hyperbolic partial differential equations[M].New York:Academic Press,1974:195-212

[25]DIENER P,DORBAND E N,SCHNETTER E,et al.Optimized high-order derivative and dissipation operators satisfying summation by parts,and applications in three-dimensional multi-block evolutions[J].Journal of Scientific Computing,2007,32(1):109-145

[26]VON NEUMANN J,RICHTMYER R D.A method for the numerical calculation of hydrodynamic shocks[J].Journal of Applied Physics,1950,21(3):232-237

[27]SJ?GREEN B,PETERSSON N A.A fourth order accurate finite difference scheme for the elastic wave equation in second order formulation[J].Journal of Scientific Computing,2012,52(1):17-48

(編輯:戴春秋)

Stability analysis of SBP finite difference scheme for two-dimensional acoustic wave equation in boundary-conforming grids

WANG Ying1,ZHOU Hui1,SHENG Shanbo2

(1.StateKeyLaboratoryofPetroleumResourceandProspecting,CNPCKeyLabofGeophysicalProspecting,ChinaUniversityofPetroleum,Beijing102249,China;2.ChinaNationalOil&GasExplorationandDevelopmentCorporation,Beijing100034,China)

Abstract:Traditionally,finite difference method is widely used as a fast and accurate method for numerical simulation and migration of wave equation.However,finite difference method faces obstacles when there are surface topography and irregular interfaces.The staircase approximations caused by regular grids influence the accuracy of finite difference method.Boundary-conforming grids by elliptic method offer an optimal alternative for finite difference wavefield simulation in irregular regions.By such grids,the calculations of spatial derivatives are transformed by a chain rule into those in the regular computational space,where traditional finite difference schemes are still applicable.The two-dimensional acoustic wave equation is reformulated in boundary-conforming grids.The Summation-by-Parts (SBP) finite difference scheme with second order accuracy is used for discretization.The stability of the discretization formula is discussed by Fourier spectral method to obtain the stability condition of the SBP finite difference scheme for two-dimensional acoustic wave equation in boundary-conforming grids.A numerical test shows that the accuracy and effectiveness can be guaranteed by reasonably choosing the parameters in numerical simulation,such as the temporal and spatial intervals.Comparisons between the numerical simulations of SBP and central finite difference method in boundary-conforming grids are performed,which shows that SBP finite difference scheme is more stable.

Keywords:boundary-conforming grids,acoustic wave equation,SBP finite difference,central finite difference,stability condition

文章編號:1000-1441(2016)01-0033-08

DOI:10.3969/j.issn.1000-1441.2016.01.005

中圖分類號:P631

文獻標識碼:A

基金項目:國家重點基礎研究發(fā)展計劃(973計劃)(2013CB228603)、國家自然科學基金(41174119)和中石油物探新方法新技術(shù)研究(2014A-3609)項目共同資助。

作者簡介:王穎(1987—),女,博士在讀,主要從事地震數(shù)據(jù)處理新方法研究。通訊作者:周輝(1966—),男,教授,博士生導師,主要從事地震正演模擬、波形反演、偏移成像研究工作。

收稿日期:2015-06-03;改回日期:2015-08-10。

王穎,周輝,盛善波.貼體坐標系二維聲波方程SBP有限差分法的穩(wěn)定性分析[J].石油物探,2016,55(1):-40

WANG Ying,ZHOU Hui,SHENG Shanbo.Stability analysis of SBP finite difference scheme for two-dimensional acoustic wave equation in boundary-conforming grids[J].Geophysical Prospecting for Petroleum,2016,55(1):-40

This research is financially supported by the National Key Basic Research Program of China (973 Program)(Grant No.2013CB228603),the National Science Foundation of China (Grant No.41174119) and CNPC Geophysical Exploration New Technology Research Program (Grant No.2014A-3609).