王官超,杜啟振

(中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580)

?

基于包絡(luò)目標(biāo)函數(shù)的彈性波波形反演

王官超,杜啟振

(中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580)

摘要:彈性波場(chǎng)的復(fù)雜性致使反演問(wèn)題非線(xiàn)性性增強(qiáng),尋優(yōu)過(guò)程極易陷入局部極值。多尺度反演能在一定程度上提高反演的穩(wěn)定性及分辨率,但多尺度反演策略無(wú)法解決實(shí)際數(shù)據(jù)的低頻缺失問(wèn)題。針對(duì)該問(wèn)題,提出了一種新的基于包絡(luò)目標(biāo)函數(shù)的彈性波多尺度反演策略。首先將包絡(luò)反演方法應(yīng)用到彈性波,恢復(fù)地下彈性介質(zhì)參數(shù)中的低波數(shù)信息；然后將包絡(luò)反演結(jié)果作為常規(guī)多尺度反演的輸入,進(jìn)一步提高反演的穩(wěn)定性。Marmousi2模型試算結(jié)果表明該策略能夠更加準(zhǔn)確地反演縱、橫波速度以及密度信息。

關(guān)鍵詞：包絡(luò)反演；低頻缺失；多尺度；彈性波全波形反演

全波形反演技術(shù)(Full Waveform Inversion,FWI)是一種從全波場(chǎng)地震記錄中定量提取地下介質(zhì)參數(shù)的數(shù)據(jù)匹配技術(shù),該技術(shù)利用波場(chǎng)的走時(shí)、振幅以及相位等多種信息,通過(guò)不同優(yōu)化方法獲得地下彈性介質(zhì)參數(shù),其分辨率明顯高于傳統(tǒng)基于射線(xiàn)理論的反演方法,反演結(jié)果不僅能為疊前深度偏移提供更為可靠的偏移速度,而且能直接提供精度較高的地下介質(zhì)彈性參數(shù)信息。但是,該技術(shù)存在計(jì)算量大、非線(xiàn)性性強(qiáng)、初始速度場(chǎng)構(gòu)建困難以及波場(chǎng)之間的不確定性等缺點(diǎn),導(dǎo)致波形反演容易陷入循環(huán)跳動(dòng)和局部極小值等問(wèn)題,當(dāng)涉及到彈性多參數(shù)反演時(shí)這些問(wèn)題表現(xiàn)得尤為嚴(yán)重,其中初始模型的準(zhǔn)確與否直接影響反演結(jié)果的正確性。

通常情況下,背景宏觀速度模型可以利用初至走時(shí)層析[1]、立體層析[2-3]、反射/折射層析、疊加速度以及偏移速度分析等方法獲得。當(dāng)初始模型遠(yuǎn)離真實(shí)模型時(shí),全波形反演目標(biāo)函數(shù)存在多極值,往往導(dǎo)致最優(yōu)化方法收斂到局部極小[4-5]。多尺度反演策略[6-8]能在一定程度上解決該問(wèn)題。但是,由于采集技術(shù)的制約，實(shí)際地震勘探中很難記錄到5Hz以下的頻率成分,準(zhǔn)確的含低波數(shù)信息的初始速度模型的構(gòu)建仍然是FWI中亟需解決的關(guān)鍵問(wèn)題。針對(duì)初始模型構(gòu)建困難這一瓶頸,許多學(xué)者進(jìn)行了大量研究,SHIN等[9-10]提出了利用Laplace域阻尼波場(chǎng)零頻分量反演長(zhǎng)波長(zhǎng)宏觀速度模型作為頻率域全波形反演的初始速度模型。事實(shí)上,將Laplace域和Fourier域的波形反演集合成一個(gè)過(guò)程,既可以提高計(jì)算效率,同時(shí)也可以提高算法穩(wěn)定性以及降低對(duì)初始速度模型的依賴(lài)性。XU[11]在反演過(guò)程中引入反偏移技術(shù)，進(jìn)行了基于反射波的波形反演研究,并在實(shí)際資料應(yīng)用中取得較好效果。胡光輝等[12]采用基于早至波的特征波波形反演建模方法解決收斂過(guò)程中由于初始模型不準(zhǔn)確導(dǎo)致的局部極小問(wèn)題,使基于特征波的全波形反演逐漸受到人們的關(guān)注。地震信號(hào)經(jīng)過(guò)包絡(luò)變換后的能量主要集中在極低頻區(qū)域,這種含有極低頻信息的包絡(luò)信號(hào)可以有效反映地下介質(zhì)的大尺度信息。因此,WU等[13-14]在聲波假設(shè)下提出基于包絡(luò)目標(biāo)函數(shù)的全波形反演方法建立低頻速度模型。通過(guò)對(duì)低頻信息進(jìn)行補(bǔ)償、頻率外推等方式,減少低頻信息缺失帶來(lái)的不利影響,同樣是一種可行策略。

基于前人的研究成果,本文提出了一種基于彈性波包絡(luò)反演方法的多尺度反演策略。首先將聲波包絡(luò)反演方法推廣應(yīng)用到彈性波,恢復(fù)地下彈性介質(zhì)參數(shù)中的大尺度信息,降低彈性多參數(shù)反演的非線(xiàn)性性；其次將彈性波包絡(luò)反演獲得的含有低頻信息的反演結(jié)果作為常規(guī)多尺度反演的初始模型,進(jìn)一步解決反演過(guò)程中存在的不穩(wěn)定性問(wèn)題。Marmousi2模型試算結(jié)果表明,結(jié)合彈性波包絡(luò)反演和多尺度反演的反演策略能夠更加準(zhǔn)確地反演縱、橫波速度以及密度信息,實(shí)現(xiàn)彈性波多參數(shù)反演。

1基于包絡(luò)目標(biāo)函數(shù)的低頻反演方法

常規(guī)全波形反演通過(guò)匹配觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù),使其數(shù)據(jù)殘差δu(t)=u(t)-d(t)達(dá)到最小,定量提取地下介質(zhì)彈性參數(shù)m,目標(biāo)函數(shù)定義為:

(1)

在彈性介質(zhì)中,多分量地震數(shù)據(jù)可通過(guò)速度應(yīng)力方程進(jìn)行數(shù)值模擬合成,即:

式中:ρ(x)是密度；u(x,t)和T(x,t)分別為速度分量和應(yīng)力分量；F(x,t)為震源函數(shù)；C(x)為彈性常數(shù)。利用高階交錯(cuò)網(wǎng)格有限差分方法進(jìn)行波動(dòng)方程求解[15],采用PML吸收邊界條件。

包絡(luò)反演與常規(guī)時(shí)間域反演的區(qū)別在于其匹配的是地震數(shù)據(jù)的包絡(luò)而非地震數(shù)據(jù)。根據(jù)聲波包絡(luò)反演的目標(biāo)函數(shù)[12]可以直接給出彈性波包絡(luò)反演目標(biāo)函數(shù):

(3)

式中:d(t)=[dx(t),dz(t)]和u(t)=[ux(t),uz(t)]分別為觀測(cè)和合成的多分量地震記錄；eobs(t)和esyn(t)分別為多分量地震記錄d(t)和u(t)的包絡(luò)；dh(t)和uh(t)為相應(yīng)波形的希爾伯特變換；E為包絡(luò)殘差；參數(shù)p控制包絡(luò)信號(hào)的能量,可以是任意正數(shù),p值越大包絡(luò)信號(hào)中深層反射的能量損失越嚴(yán)重,導(dǎo)致包絡(luò)反演無(wú)法刻畫(huà)深層結(jié)構(gòu),在本文的數(shù)值測(cè)試中p=2。

通過(guò)推導(dǎo)彈性波包絡(luò)目標(biāo)函數(shù)((3)式)相對(duì)于模型參數(shù)m的一階偏導(dǎo)可得:

(4)

(4)式可表示為雅克比矩陣J以及有效殘差向量η的乘積:

(5)

(6a)

(6b)

(6c)

式中:v=(vx,vz),w=(wx,wz,wxx,wzz,wzx)分別代表正傳波場(chǎng)和反傳殘差波場(chǎng)。因?yàn)椴煌瑓?shù)之間存在著相互關(guān)系,目標(biāo)函數(shù)對(duì)于其它相關(guān)參數(shù)m=[vP,vS,ρ]的梯度可以通過(guò)鏈?zhǔn)椒▌t進(jìn)行求解。另外,迭代過(guò)程采用拋物線(xiàn)插值方法估計(jì)更新步長(zhǎng)。

包絡(luò)反演假設(shè)地震波形為調(diào)幅信號(hào),包絡(luò)算子是一種非線(xiàn)性的解調(diào)算子,可反映隱藏在調(diào)幅信號(hào)中的低頻調(diào)制信號(hào)的波形變化。圖1a和圖1b分別展示了兩個(gè)余弦衰減信號(hào)疊加的調(diào)幅信號(hào)和單道地震記錄的波形及其包絡(luò),圖1c和圖1d為對(duì)應(yīng)的頻譜特征。從頻譜上可看出,包絡(luò)變換后信號(hào)的能量主要集中在低頻區(qū)域,橫坐標(biāo)越接近0，振幅譜能量就越強(qiáng)(圖1c和圖1d),包絡(luò)反演正是利用這些極低頻信息來(lái)恢復(fù)地下介質(zhì)參數(shù)的長(zhǎng)波長(zhǎng)分量。

圖1　調(diào)幅信號(hào)及其包絡(luò)(a),單道地震記錄及其包絡(luò)(b),對(duì)應(yīng)圖1a的頻譜特征(c)以及對(duì)應(yīng)圖1b的頻譜特征(d)

2基于包絡(luò)目標(biāo)函數(shù)的多尺度反演策略

BUNKS等[6]采用Hamming窗低通濾波器對(duì)地震子波和數(shù)據(jù)進(jìn)行時(shí)間尺度分解,首次實(shí)現(xiàn)了時(shí)間域多尺度全波形反演策略。BOONYASIRIWAT等[19]指出利用Hamming窗進(jìn)行濾波會(huì)有高頻泄漏。為解決Hamming窗濾波的高頻泄漏問(wèn)題,可以選用其它無(wú)泄漏或少泄漏的濾波窗函數(shù)。本文中采用維納濾波器進(jìn)行低通濾波,實(shí)現(xiàn)多尺度反演策略,維納濾波器表達(dá)式為:

(7)

式中:fWiener是維納濾波器；WOriginal是原始信號(hào)；WTarget為期望得到的信號(hào)；ω是角頻率；ε是為了避免不穩(wěn)定現(xiàn)象發(fā)生而設(shè)定的一個(gè)很小的值。利用維納濾波器將波場(chǎng)數(shù)據(jù)分解成不同頻帶,由低頻至高頻帶進(jìn)行逐頻帶反演，從而實(shí)現(xiàn)多尺度反演策略,降低反演過(guò)程陷入局部極小值的風(fēng)險(xiǎn)。圖2a為利用20Hz Ricker子波生成的單炮Z分量記錄,圖2b 為利用5Hz Ricker子波對(duì)圖2a的記錄進(jìn)行維納濾波后的地震記錄。由于實(shí)際地震數(shù)據(jù)往往缺失低頻信息,基于低通濾波的常規(guī)多尺度反演策略無(wú)法反演模型參數(shù)中的大尺度信息,導(dǎo)致反演陷入局部極值或循環(huán)跳動(dòng)。

圖2　利用20Hz Ricker子波生成的單炮記錄(a)以及利用5Hz Ricker子波對(duì)圖2a進(jìn)行維納濾波后的單炮記錄(b)

由于包絡(luò)算子具有解調(diào)地震數(shù)據(jù)中隱含低頻信息的作用，因此，利用包絡(luò)算子解調(diào)出的低頻信息可以有效恢復(fù)地下介質(zhì)的長(zhǎng)波長(zhǎng)成分,進(jìn)而解決常規(guī)多尺度反演無(wú)法克服的低頻缺失或不足問(wèn)題。同時(shí),包絡(luò)變換后的信號(hào)在富含低頻信息的同時(shí)也存在相應(yīng)的高頻信息。圖3是主頻為7Hz(實(shí)線(xiàn)和虛線(xiàn))、20Hz(點(diǎn)線(xiàn)和點(diǎn)劃線(xiàn))Ricker子波及其包絡(luò)的頻譜響應(yīng),可以看出包絡(luò)信號(hào)攜帶的高頻信息與原始信號(hào)主頻帶的高低是相對(duì)的。如果原信號(hào)主頻較高，變換后的信號(hào)頻率也相對(duì)較高，反之亦然。

圖3　20Hz Ricker子波及其包絡(luò)的頻譜響應(yīng)以及7Hz Ricker子波及其包絡(luò)的頻譜響應(yīng)

我們通過(guò)數(shù)值算例來(lái)說(shuō)明包絡(luò)中高頻信息對(duì)包絡(luò)反演結(jié)果的影響。圖4是以線(xiàn)性模型作為初始模型分別利用圖3中兩種震源子波進(jìn)行彈性包絡(luò)反演獲得的縱波速度,可以看出,當(dāng)高頻信息較多時(shí)，包絡(luò)反演結(jié)果不穩(wěn)定(圖4b黑色矩形框)。為提高包絡(luò)反演的穩(wěn)定性,在進(jìn)行包絡(luò)反演之前利用類(lèi)似于時(shí)間域多尺度反演對(duì)原始地震數(shù)據(jù)進(jìn)行低通濾波處理,這樣就把彈性包絡(luò)反演與常規(guī)多尺度反演有機(jī)結(jié)合在一起,自適應(yīng)的實(shí)現(xiàn)了一種從極低頻開(kāi)始的多尺度反演策略。

因此,基于包絡(luò)目標(biāo)函數(shù)的多尺度反演策略的具體實(shí)現(xiàn)過(guò)程為:首先,對(duì)地震數(shù)據(jù)進(jìn)行時(shí)間尺度分解；其次,采用低通濾波后的地震數(shù)據(jù)實(shí)施包絡(luò)反演,恢復(fù)地下彈性介質(zhì)參數(shù)中的大尺度信息；最后,將彈性波包絡(luò)反演結(jié)果作為初始模型進(jìn)行常規(guī)多尺度反演以獲取精確的最終反演結(jié)果。

圖4　彈性波包絡(luò)反演得到的縱波速度a 7Hz震源子波; b 20Hz震源子波

3Marmousi2模型測(cè)試

在反演測(cè)試過(guò)程中采用兩種不同類(lèi)型的震源子波合成多分量地震數(shù)據(jù)(主頻為20Hz的全頻帶Ricker子波和低頻缺失(<7Hz)Ricker子波)來(lái)驗(yàn)證彈性波包絡(luò)反演恢復(fù)地下模型參數(shù)大尺度信息的能力及其在地震數(shù)據(jù)缺失低頻情況下的適用性。為保證包絡(luò)反演的穩(wěn)定性,首先利用7Hz Ricker子波對(duì)兩種類(lèi)型震源子波生成的地震記錄進(jìn)行維納低通濾波,然后進(jìn)行包絡(luò)反演。圖6和圖7 為以線(xiàn)性模型作為初始模型進(jìn)行包絡(luò)反演迭代30次后得到的平滑背景場(chǎng)?？梢钥吹?針對(duì)兩種不同震源子波包絡(luò)反演獲得的背景模型構(gòu)造非常相近,這說(shuō)明了包絡(luò)反演對(duì)低頻信息的依賴(lài)性很弱,即使在低頻缺失的情況下同樣可以反演得到模型中的低頻信息,為常規(guī)多尺度波形反演提供含有低頻信息的初始速度模型,降低全波形反演的非線(xiàn)性性。

圖5　真實(shí)Marmousi2縱波速度模型(a)及其線(xiàn)性模型(b)

圖6　線(xiàn)性模型作為初始模型的彈性波包絡(luò)反演30次迭代后的結(jié)果(全頻帶震源子波)a 縱波速度; b 橫波速度; c 密度

圖7　線(xiàn)性模型作為初始模型的彈性波包絡(luò)反演30次迭代后的結(jié)果(低頻缺失震源子波)a 縱波速度; b 橫波速度; c 密度

圖8和圖9分別為以圖6和圖7中包絡(luò)反演結(jié)果作為初始模型,選取5,10,20Hz的Ricker子波對(duì)地震記錄進(jìn)行維納濾波進(jìn)而實(shí)施多尺度反演得到的最終結(jié)果。相比于直接利用線(xiàn)性模型作為初始模型進(jìn)行多尺度反演的結(jié)果(圖10和圖11),利用包絡(luò)反演加多尺度反演得到的最終結(jié)果在構(gòu)造(黑色橢圓處)上得到了明顯提高,構(gòu)造假象明顯減少。對(duì)比圖10和圖11可以看出,利用低頻缺失震源子波的反演結(jié)果相比于全頻帶震源子波反演結(jié)果的構(gòu)造假象更為嚴(yán)重,這說(shuō)明傳統(tǒng)全波形反演對(duì)震源的頻譜特性更為敏感。

圖8　以圖6所示包絡(luò)反演結(jié)果作為初始模型得到的多尺度反演結(jié)果(全頻帶震源子波)a 縱波速度; b 橫波速度; c 密度

圖9　以圖7所示包絡(luò)反演結(jié)果作為初始模型得到的多尺度反演結(jié)果(低頻缺失震源子波)a 縱波速度; b 橫波速度; c 密度

圖10　以線(xiàn)性模型作為初始模型得到的多尺度反演結(jié)果(全頻帶震源子波)a 縱波速度; b 橫波速度; c 密度

圖11　以線(xiàn)性模型作為初始模型得到的多尺度反演結(jié)果(低頻缺失震源子波)a 縱波速度; b 橫波速度; c 密度

4結(jié)論

通過(guò)理論研究與模型試算得到以下認(rèn)識(shí):①包絡(luò)算子可以解調(diào)出地震數(shù)據(jù)中隱含的低頻信息。包絡(luò)反演利用解調(diào)出的低頻信息能有效恢復(fù)模型中的大尺度信息,在地震數(shù)據(jù)缺失低頻的情況下同樣適用。②包絡(luò)信號(hào)在富含低頻信息的同時(shí)也存在相應(yīng)的高頻信息。高頻信息較多時(shí)會(huì)導(dǎo)致包絡(luò)反演不穩(wěn)定,通過(guò)將包絡(luò)反演與基于低通濾波的多尺度反演相結(jié)合,自適應(yīng)的實(shí)現(xiàn)了一種從極低頻開(kāi)始的多尺度反演策略,最大限度地降低了反演對(duì)初始模型的依賴(lài)性。

實(shí)際應(yīng)用是全波形反演的挑戰(zhàn),野外數(shù)據(jù)低頻缺失或不足是影響波形反演獲得全波數(shù)帶參數(shù)估計(jì)結(jié)果的主要因素。包絡(luò)反演技術(shù)可恢復(fù)隱藏在地震資料中的低頻信息,建立低頻初始模型,具有很好的發(fā)展前景。

參考文獻(xiàn)

[1]NOLET G.Seismic tomography:with applications in global seismology and exploration geophysics [M].Bonston:Kluwer Academic Publishers,1987,1-24

[2]WOODWARD M J.Wave-equation tomography[J].Geophysics,1992,57(1):15-26

[3]PRIEUX V,OPERTO S,BROSSIER R,et al.Application of acoustic full waveform inversion to the synthetic Valhall model[J].Expanded Abstracts of 77thAnnual Internat SEG,Mtg,2007:2268-2272

[4]VIRIEUX J,OPERTO S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

[5]楊勤勇,胡光輝,王立歆.全波形反演研究現(xiàn)狀及發(fā)展趨勢(shì)[J].石油物探,2014,53(1):77-83

YANG Q Y,HU G H,WANG L X.Research status and development trend of full waveform inversion [J].Geophysical Prospecting for Petroleum,2014,53 (1):77-83

[6]BUNKS C,SALECK F M,ZALESKI S,et al.Multiscale seismic waveform inversion[J].Geophysics,1995,60(5):1457-1473

[7]SIRGUE L,PRATT R G.Efficient waveform inversion and imaging:a strategy for selecting temporal frequencies[J].Geophysics,2004,69(1):231-248

[8]李媛媛,李振春,張凱.頻率域多尺度彈性波全波形反演[J].石油物探,2015,54(3):317-323

LI Y Y,LI Z C,ZHANG K.Multi-scale elastic full waveform inversion in frequency domain [J].Geophysical Prospecting for Petroleum,2015,54 (3):317-323

[9]SHIN C,HA W.A comparison between the behavior of objective functions for waveform inversion in the frequency and Laplace domains[J].Geophysics,2008,73(5):VE119-VE133

[10]SHIN C,CHA Y H.Waveform inversion in the Laplace domain [J].Geophysical Journal International,2008,173(3):922-931

[11]XU S,WANG D,CHEN F,et al.Full waveform inversion for reflected seismic data[J].Expanded Abstracts of 74thEAGE Annual Conference,2012:114-117

[12]胡光輝,王立歆,王杰,等.基于早至波的特征波波形反演建模方法[J].石油物探,2015,54 (1):71-76

HU G H,WANG L X,WANG J,et al.Characteristics waveform inversion based on early arrival waves [J].Geophysical Prospecting for Petroleum,2015,54 (1):71-76

[13]WU R S,LUO J R,WU B.Seismic envelope inversion and modulation signal model[J].Geophysics,2014,79(3):WA13-WA24

[14]LUO J R,WU R S.Seismic envelope inversion:reductionoflocal minima and noise resistance [J].Geophysical Prospecting,2015,63(3):597-614

[15]董良國(guó),馬在田,曹景忠.一階彈性波交錯(cuò)網(wǎng)格高階有限差分解法[J].地球物理學(xué)報(bào),2000,43(3):411-419

DONG L G,MA Z T,CAO J Z.Higher order difference method for first order elastic wave equation with staggered grid [J].Chinese Journal of Geophysics,2000,43(3):411-419

[16]TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[17]PLESSIX R E.A review of the adjoint-state method for computing the gradient of a functional with geophysical applications [J].Geophysical Journal International,2006,167(2):495-503

[18]WANG J,ZHOU H,TIAN Y,et al.A new scheme for elastic full waveform inversion based on velocity-stress wave equations in time domain[J].Expanded Abstracts of 82ndAnnual Internat SEG,Mtg,2012:1-5

[19]BOONYASIRIWAT C,VALASEK P,ROUTH P,et al.An efficient multiscale method for time-domain waveform tomography[J].Geophysics,2009,74(6):WCC59-WCC68

(編輯：朱文杰)

Elastic full waveform inversion based on envelope objective function

WANG Guanchao,DU Qizhen

(SchoolofGeosciences,ChinaUniversityofPetroleum,Qingdao266580,China)

Abstract:The complexity of elastic wavefield will increase the nonlinearity of inversion and make the inverse problem fall into local minima.Multi-scale inversion algorithms can improve the stability and resolution of inversion to a certain extent,but it cannot solve the problem of low frequency missing in real data.To alleviate this problem,we propose a new multi-scale scheme based on envelope inversion to obtain the elastic parameters accurately.The main procedures are as follows:Firstly,we extend the envelope inversion method to elastic medium to estimate the low wavenumber components of elastic parameters; then,the envelope inversion results are used as the initial model of conventional multi-scale inversion to further improve the stability of the inversion.The Marmousi2 model test results show that the new multi-scale scheme can estimate P- and S- wave velocities more accurately,as well as the density information.

Keywords:envelope inversion,low frequency missing,multi-scale,elastic full waveform inversion

文章編號(hào)：1000-1441(2016)01-0133-09

DOI:10.3969/j.issn.1000-1441.2016.01.017

中圖分類(lèi)號(hào)：P631

文獻(xiàn)標(biāo)識(shí)碼：A

基金項(xiàng)目：國(guó)家自然科學(xué)基金(41174100和41074087)和中國(guó)石化地球物理重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金聯(lián)合資助。

作者簡(jiǎn)介：王官超(1990—),男,碩士在讀,主要從事全波形反演研究工作。通訊作者：杜啟振(1969—),男,教授,博士生導(dǎo)師,主要從事地震彈性波傳播理論與成像等研究工作。

收稿日期：2015-07-07;改回日期：2015-08-24。

王官超,杜啟振.基于包絡(luò)目標(biāo)函數(shù)的彈性波波形反演[J].石油物探,2016,55(1):-141

WANG Guanchao,DU Qizhen.Elastic full waveform inversion based on envelope objective function[J].Geophysical Prospecting for Petroleum,2016,55(1):-141

This research is financially supported by the National Science Foundation of China (Grant Nos.41174100,41074087) and the SINOPEC Key Laboratory of Geophysics Project.