地震正演與高斯束偏移

地震成像與速度反演李振春中國(guó)石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院提綱一、雙變網(wǎng)格正演模擬二、高斯束偏移三、波動(dòng)方程偏移四、繞射波成像五、多次波成像六、層析速度分析與全波形反演雙變網(wǎng)格正演模擬算法主要包括：基于聲波、彈性波、各向異性介質(zhì)的2D/3D有限差分算法，基于粘聲、粘彈介質(zhì)的2D正演模擬算法。適用于起伏地表模擬的變網(wǎng)格算法、2D譜元素法及偽譜法；適用于裂縫型介質(zhì)的雙變網(wǎng)格正演模擬算法；以及較為完善的邊界處理方法包括：吸收衰減邊界條件、PML以及CPML邊界條件等。2D3D聲波粘彈性/各向異性彈性介質(zhì)模型傳統(tǒng)有限差分地震波模擬方法基于規(guī)則的矩形網(wǎng)格，是處理非均勻介質(zhì)，特別是連續(xù)非均勻介質(zhì)波動(dòng)問題的有效方法！雙變網(wǎng)格正演模擬原理

難題2難題1為保證計(jì)算精度和計(jì)算穩(wěn)定性，需要精細(xì)網(wǎng)格剖分，導(dǎo)致計(jì)算量增加和局部過采樣需要發(fā)展等效模型，將非均質(zhì)體的物性參數(shù)與波動(dòng)方程聯(lián)系起來速度場(chǎng)離散方式1速度場(chǎng)離散方式2雙變網(wǎng)格算法即根據(jù)地震地質(zhì)模型將速度場(chǎng)進(jìn)行不同尺度的網(wǎng)格剖分，即變網(wǎng)格，同時(shí)又結(jié)合局部變時(shí)間算法，提高精度和效率。模型網(wǎng)格劃分的優(yōu)劣與波場(chǎng)正演的效率密切相關(guān)。

速度場(chǎng)多尺度離散策略Virieux（1986）和Graves（1996b）提出的一階速度應(yīng)力方程，各向同性非均勻介質(zhì)方程可寫為：變網(wǎng)格彈性波方程離散化原理在傳統(tǒng)交錯(cuò)網(wǎng)格技術(shù)中，變量關(guān)于的2n階空間偏導(dǎo)數(shù)差分近似式為交錯(cuò)網(wǎng)格高階差分格式

差分系數(shù)是常數(shù)差分算子是固定步長(zhǎng)彈性波方程離散化在空間變網(wǎng)格技術(shù)中，變量關(guān)于的2n階空間偏導(dǎo)數(shù)差分近似式為變交錯(cuò)網(wǎng)格高階差分格式

差分系數(shù)是網(wǎng)格步長(zhǎng)的函數(shù)差分算子是空間變化的算法優(yōu)點(diǎn)：

在全局時(shí)間計(jì)算時(shí)刻粗細(xì)網(wǎng)格過渡區(qū)域不需要重復(fù)設(shè)置計(jì)算區(qū)域空間的差分精度不需要有任何改變只需要前期求解出差分系數(shù)，便可以在每一計(jì)算時(shí)刻引用彈性波方程離散化

空間網(wǎng)格的變化需要時(shí)間采樣滿足精細(xì)網(wǎng)格采樣對(duì)應(yīng)的穩(wěn)定性條件！全局時(shí)間采樣不需要滿足最短波長(zhǎng)原則，同傳統(tǒng)交錯(cuò)網(wǎng)格算法的穩(wěn)定性條件是一致的。但是在全局時(shí)間層內(nèi)，局部精細(xì)網(wǎng)格的時(shí)間計(jì)算步長(zhǎng)需要滿足最短波長(zhǎng)原則。雙變網(wǎng)格算法局部可變的時(shí)間采樣步長(zhǎng)首先給出一些假設(shè)：已知時(shí)刻的波場(chǎng)值已知時(shí)刻的波場(chǎng)值

已知時(shí)刻的波場(chǎng)值已知時(shí)刻的波場(chǎng)值

具體實(shí)現(xiàn)步驟如下：1：粗網(wǎng)格平面的求解時(shí)刻的波場(chǎng)值時(shí)刻的波場(chǎng)值2：利用過渡區(qū)域的波場(chǎng)值判斷是否進(jìn)入精細(xì)時(shí)間采樣局部時(shí)間采樣變化（LVTS）首先給出一些假設(shè)：已知時(shí)刻的波場(chǎng)值已知時(shí)刻的波場(chǎng)值

已知時(shí)刻的波場(chǎng)值已知時(shí)刻的波場(chǎng)值

具體實(shí)現(xiàn)步驟如下：3：精細(xì)時(shí)間采樣平面內(nèi)的精細(xì)網(wǎng)格點(diǎn)波場(chǎng)值求解時(shí)刻的波場(chǎng)值利用雙線性公式：獲得精細(xì)時(shí)間采樣平面邊界處的波場(chǎng)值4：波場(chǎng)傳遞局部時(shí)間采樣變化（LVTS）速度場(chǎng)網(wǎng)格分為粗糙網(wǎng)格、精細(xì)網(wǎng)格，中間存在過渡區(qū)域時(shí)間層波場(chǎng)分為全局時(shí)間層波場(chǎng)和精細(xì)時(shí)間層波場(chǎng)總波場(chǎng)滿足彈性波波動(dòng)方程時(shí)間層局部波場(chǎng)求解空間精細(xì)網(wǎng)格剖分變交錯(cuò)網(wǎng)格求解時(shí)間層全局波場(chǎng)求解空間粗糙網(wǎng)格剖分傳統(tǒng)交錯(cuò)網(wǎng)格求解波場(chǎng)傳遞總波場(chǎng)算法流程模型參數(shù)：水平距離：401*6m縱向深度：400*6m精細(xì)網(wǎng)格采樣：2m全局時(shí)間采樣：0.3ms局部時(shí)間采樣：0.1ms炮點(diǎn)位置：（201*6m，5*6m）模擬時(shí)間長(zhǎng)度：2.5s傳統(tǒng)算法x分量雙變算法x分量傳統(tǒng)算法z分量雙變算法z分量典型模型一：低降速帶模型2、雙變網(wǎng)格模型試算

對(duì)比可知：基于相同的常網(wǎng)格大小，而對(duì)目標(biāo)低速區(qū)域采用精細(xì)網(wǎng)格剖分的雙變網(wǎng)格算法，能夠較好的壓制頻散，得到振幅連續(xù)性更好的同相軸分布。速度場(chǎng)模型典型模型二：局部低速體綠色框內(nèi)表示雙變區(qū)域模型參數(shù)：水平距離：930*6m縱向深度：550*6m精細(xì)網(wǎng)格采樣：2m全局時(shí)間采樣：0.3ms局部時(shí)間采樣：0.06ms炮點(diǎn)位置：（101*6m，3*6m）模擬時(shí)間長(zhǎng)度：2.4s左圖雙變算法波場(chǎng)快照中圖常網(wǎng)格波場(chǎng)快照右圖二者波場(chǎng)快照之差

波場(chǎng)快照對(duì)比圖

由于采用了局部精細(xì)化的雙變網(wǎng)格技術(shù)，在局部目標(biāo)區(qū)域網(wǎng)格可以從米級(jí)變化到毫米級(jí)（而在全局采用毫米級(jí)的常網(wǎng)格算法將導(dǎo)致巨大的計(jì)算量，因而是不可行的），從波場(chǎng)快照的對(duì)比圖可以看出：雙變網(wǎng)格算法能夠較好模擬局部小尺度低速構(gòu)造。Z分量X分量左圖雙變算法波場(chǎng)快照中圖常網(wǎng)格波場(chǎng)快照右圖二者波場(chǎng)快照之差

炮記錄對(duì)比圖

對(duì)比可知：兩種網(wǎng)格剖分方式的炮記錄差異主要來自于局部小尺度的低速異常構(gòu)造，為實(shí)現(xiàn)高精度的正演模擬及偏移成像，基于變網(wǎng)格剖分的雙變網(wǎng)格算法，具有較好的應(yīng)用前景。Z分量X分量雙變網(wǎng)格比精細(xì)網(wǎng)格、精細(xì)采樣節(jié)省80.1%的時(shí)間比局部精細(xì)網(wǎng)格，全局采樣節(jié)省60.3%的時(shí)間從計(jì)算耗時(shí)對(duì)比可知：在保證計(jì)算精度的前提下，模型網(wǎng)格非均勻剖分的雙變網(wǎng)格算法，能夠極大的提高正演模擬的效率。雙變區(qū)域模型參數(shù)：水平距離：601*8m縱向深度：700*6m精細(xì)網(wǎng)格采樣：2m全局時(shí)間采樣：0.5ms局部時(shí)間采樣：0.16ms炮點(diǎn)位置：（301*8m，10*6m）模擬時(shí)間長(zhǎng)度：3.5s典型模型三：陡坡帶砂礫巖體薄層圖示速度場(chǎng)最薄層厚30-60m左圖雙變算法波場(chǎng)快照中圖常網(wǎng)格波場(chǎng)快照右圖二者波場(chǎng)快照之差

波場(chǎng)快照對(duì)比圖Z分量X分量雙變網(wǎng)格算法能夠較好的模擬出來自陡坡帶砂礫巖體的地震波響應(yīng)，振幅更加均勻、同相軸一致性更好。左圖雙變炮記錄中圖常網(wǎng)格炮記錄右圖二者炮記錄之差Z分量X分量雙變網(wǎng)格算法能夠展示更多的來自于地下細(xì)微構(gòu)造的波場(chǎng)信息，局部分辨率較常網(wǎng)格大大提高。左圖全局速度場(chǎng)右圖含裂縫局部速度場(chǎng)圖示含裂縫速度場(chǎng)裂縫參數(shù)：裂縫位于3000m深度，密度10m/條，共200個(gè)裂縫，縫寬2mm典型模型五：裂縫發(fā)育層左圖無縫炮記錄中圖裂縫延伸150m炮記錄右圖裂縫對(duì)應(yīng)某時(shí)刻波場(chǎng)快照裂縫頂部裂縫底部上界面反射P波上界面反射S波裂縫的波場(chǎng)響應(yīng)從波場(chǎng)快照及單炮記錄可以看出：雙變網(wǎng)格高階有限差分算法，能夠較好的模擬毫米級(jí)（2mm）的裂縫型儲(chǔ)層，在裂縫密度達(dá)到一定程度時(shí)，裂縫地震響應(yīng)特征較為清晰?；陔p變網(wǎng)格算法的小結(jié)1、雙變網(wǎng)格有限差分方法，依據(jù)模型特征采用非均勻網(wǎng)格剖分，可以在不顯著增加計(jì)算成本的基礎(chǔ)上，極大提高目標(biāo)研究區(qū)域的正演模擬精度并抑制計(jì)算頻散；2、雙變網(wǎng)格有限差分方法，能夠較好的對(duì)毫米級(jí)裂縫進(jìn)行正演模擬，裂縫波場(chǎng)響應(yīng)特征清晰；3、雙變網(wǎng)格有限差分方法，較容易與物理模擬方法結(jié)合，并為物理模擬提供技術(shù)支持；4、雙變網(wǎng)格有限差分方法，可以進(jìn)一步發(fā)展為雙變網(wǎng)格的RTM，極大提高RTM的成像效率，為RTM的生產(chǎn)實(shí)用化提供幫助。提綱一、雙變網(wǎng)格正演模擬二、高斯束偏移三、波動(dòng)方程偏移四、繞射波成像五、多次波成像六、層析速度分析與全波形反演地震偏移成像方法高斯束偏移成像波動(dòng)方程的數(shù)值解法單程波動(dòng)方程偏移雙程波動(dòng)方程偏移幾何射線類方法Kirchhoff偏移（復(fù)雜構(gòu)造難以成像）高斯束偏移（高效靈活精確）缺陷：計(jì)算量巨大主要內(nèi)容射線類偏移高斯束偏移★高斯束的基本性質(zhì)★高斯束偏移成像技術(shù)★高斯束偏移共成像點(diǎn)道集

高斯束逆時(shí)偏移★聲波高斯束逆時(shí)偏移★彈性波高斯束逆時(shí)偏移高斯束偏移成像垂向傳播的高斯束為虛數(shù)代表空間局部化、方向局部化的地震波場(chǎng)1、高斯束的基本性質(zhì)高斯束偏移成像高斯束偏移成像1、高斯束的基本性質(zhì)高斯束表示球面波高斯束表示平面波高斯束偏移成像高斯束偏移發(fā)展歷程GaussianbeammigrationRossHill,1990,GeophysicsPrestackGaussian-beamdepthmigrationRossHill,2001,GeophysicsGaussianbeammigrationofcommonshotrecordingSamuelGray,2005,GeophysicsTrueamplitudeGaussianbeammigrationSamuelGray,NormanBleistein,2009,Geophysics高斯束偏移成像2、高斯束偏移基本過程1）不同道集成像公式2）保真振幅高斯束偏移3）復(fù)雜地表高斯束偏移4）彈性波多分量高斯束偏移5）3D高斯束偏移高斯束偏移成像高斯束疊后偏移加窗局部?jī)A斜疊加疊后偏移高斯束偏移成像高斯束疊后偏移洼陷模型疊后記錄

30個(gè)波束最終成像結(jié)果高斯束偏移成像常規(guī)高斯束偏移計(jì)算流程震源成像點(diǎn)高斯束疊前偏移高斯束偏移成像不同道集成像過程共炮道集成像公式共接收點(diǎn)道集成像公式高斯束偏移成像不同道集成像過程共偏移距道集成像公式共偏移距道集傾斜疊加高斯束偏移成像不同道集成像過程共中心點(diǎn)道集成像公式共中心點(diǎn)道集傾斜疊加高斯束偏移成像共偏移距道集高斯束偏移Marmousimodel共炮集高斯束偏移Kirchhoff偏移不同道集成像過程高斯束偏移成像對(duì)于每個(gè)束中心位置{

加窗局部?jī)A斜疊加并計(jì)算由震源和接收點(diǎn)出射的高斯束對(duì)于每個(gè)震源射線參數(shù){

對(duì)于每個(gè)接收點(diǎn)射線參數(shù){

對(duì)于局部孔徑內(nèi)的每個(gè)成像點(diǎn){

根據(jù)高斯束的走時(shí)以及振幅利用所對(duì)應(yīng)的傾斜疊加道進(jìn)行成像

}成像點(diǎn)循環(huán)結(jié)束

}接收點(diǎn)高斯束循環(huán)結(jié)束

}震源高斯束循環(huán)結(jié)束}束中心循環(huán)結(jié)束共炮集偏移局部孔徑的成像運(yùn)算偏移計(jì)算量為：束中心為nb，震源高斯束為nps,接收點(diǎn)高斯束為npr全波至算法計(jì)算流程全波至算法高效算法高斯束偏移成像對(duì)于每個(gè)束中心位置{

加窗局部?jī)A斜疊加并計(jì)算由震源和接收點(diǎn)出射的高斯束對(duì)于每個(gè)中心點(diǎn)射線參數(shù){

對(duì)于每個(gè)粗網(wǎng)格點(diǎn){

掃描并確定使虛值走時(shí)最小的偏移距射線參數(shù)

}粗網(wǎng)格循環(huán)結(jié)束對(duì)于孔徑內(nèi)的每個(gè)成像點(diǎn){

根據(jù)高斯束的走時(shí)振幅信息利用所對(duì)應(yīng)的傾斜疊加道進(jìn)行成像

}成像點(diǎn)循環(huán)結(jié)束

}中心點(diǎn)射線參數(shù)循環(huán)結(jié)束}束中心循環(huán)結(jié)束局部孔徑的成像運(yùn)算偏移計(jì)算量為：高頻復(fù)值振蕩積分最速下降法高效算法計(jì)算流程(Hill,Geophysics,2001)算法比較高斯束偏移成像Sigsbee2amodel全波至算法（182.73s）高效算法（51.03s）算法比較高斯束偏移成像全波至算法高效算法常規(guī)Kirchhoff偏移單程波偏移高斯束偏移成像2、高斯束偏移基本過程1）不同道集成像公式2）保真振幅高斯束偏移3）復(fù)雜地表高斯束偏移4）彈性波多分量高斯束偏移5）3D高斯束偏移保幅高斯束偏移高斯束偏移成像保幅偏移定義

偏移過程中消除地震波幾何擴(kuò)散造成的能量損失，使得成像振幅正比于界面的反射系數(shù)。保幅偏移應(yīng)用條件

保幅成像條件+保幅的延拓算子。反褶積成像條件(Claerbout,Geophysics1971)互相關(guān)成像條件(Claerbout,Geophysics1971)高斯束偏移成像炮域保幅角度域保幅（SamuelGray,NormanBleistein,Geophysics,2009)保幅高斯束偏移高斯束偏移成像單炮記錄提取的反射系數(shù)四層模型（常速變密度）保幅高斯束偏移保幅高斯束偏移高斯束偏移成像保幅高斯束偏移單炮成像結(jié)果（互相關(guān)）角度域共成像道集（互相關(guān)）歸一化反射系數(shù)歸一化反射系數(shù)高斯束偏移成像保幅高斯束偏移2D鹽丘模型常規(guī)高斯束偏移保幅高斯束偏移（反褶積）高斯束偏移成像保幅高斯束偏移常規(guī)高斯束偏移保幅高斯束偏移（互相關(guān)）CDP=350CDP=450CDP=530CDP=640高斯束偏移成像30~40°20~30°0~10°10~20°不同角度范圍的疊加成像結(jié)果保幅高斯束偏移高斯束偏移成像2、高斯束偏移基本過程1）不同道集成像公式2）保真振幅高斯束偏移3）復(fù)雜地表高斯束偏移4）彈性波多分量高斯束偏移5）3D高斯束偏移復(fù)雜地表?xiàng)l件下地震成像方法波動(dòng)方程基準(zhǔn)面校正靜校正法波場(chǎng)沿拓法射線類偏移方法準(zhǔn)確性差計(jì)算量大，靈活性差，對(duì)數(shù)據(jù)要求高靈活，高效復(fù)雜地表高斯束偏移高斯束偏移成像高斯束偏移成像局部靜校正法基本原理（Gray，Geophysics，2005）復(fù)雜地表高斯束偏移RayleighII積分（岳玉波等，地球物理學(xué)報(bào)，2011第三期）復(fù)雜地表保幅成像公式起伏地表平面波分解

起伏地表模型單炮記錄復(fù)雜地表高斯束偏移高斯束偏移成像沿反射層提取的歸一化振幅曲線保幅高斯束偏移常規(guī)高斯束偏移復(fù)雜地表高斯束偏移高斯束偏移成像加拿大逆掩斷層模型單炮記錄地表高程高斯束偏移成像復(fù)雜地表高斯束偏移Kirchhoff偏移局部靜校正法保幅延拓法“直接下延”法波動(dòng)方程偏移高斯束偏移成像復(fù)雜地表高斯束偏移速度模型地表高程保幅延拓法“直接下延”法波動(dòng)方程偏移復(fù)雜地表高斯束偏移高斯束偏移成像高斯束偏移成像2、高斯束偏移基本過程1）不同道集成像公式2）保真振幅高斯束偏移3）復(fù)雜地表高斯束偏移4）彈性波多分量高斯束偏移5）3D高斯束偏移（1）彈性波Kirchhoff偏移具有極高的計(jì)算效率；成像噪聲壓制；容易解決轉(zhuǎn)換波極性問題。精度不高，不能對(duì)多次波至進(jìn)行成像。（2）彈性波逆時(shí)偏移具有極高的成像精度；計(jì)算量巨大；難以壓制不同波型的噪聲串?dāng)_；難以解決轉(zhuǎn)換波極性問題。（3）彈性波高斯束偏移兼具較高成像精度與計(jì)算效率。彈性波高斯束偏移高斯束偏移成像公式推導(dǎo)彈性波Kirchhoff-Helmohltz積分彈性動(dòng)力學(xué)高斯束表示格林張量P波位移S波位移權(quán)系數(shù)高斯束偏移成像彈性波高斯束偏移P-P成像公式P-S成像公式沿不同波型極性方向投影互相關(guān)成像條件求取標(biāo)量成像值高斯束偏移成像彈性波高斯束偏移

速度模型（縱橫波速度比:1.5/1）垂直分量水平分量高斯束偏移成像彈性波高斯束偏移單分量高斯束偏移彈性波高斯束偏移P-SP-SP-PP-PZ分量X分量高斯束偏移成像彈性波高斯束偏移(縱橫波速度比為：1.7/1.0)斜層模型高斯束偏移成像彈性波高斯束偏移彈性波高斯束偏移P垂直分量水平分量S高斯束偏移成像彈性波高斯束偏移彈性波高斯束偏移P-SP-P高斯束偏移成像彈性波高斯束偏移高斯束偏移成像2、高斯束偏移基本過程1）不同道集成像公式2）保真振幅高斯束偏移3）復(fù)雜地表高斯束偏移4）彈性波多分量高斯束偏移5）3D高斯束偏移SEG/EAGEC3-NA模型速度模型觀測(cè)系統(tǒng)第1300炮記錄模型大小：676*676*210截?。?51*481*210數(shù)據(jù)大小：6.25G偏移孔徑：Inline2000m,Crossline1000m計(jì)算時(shí)間：100個(gè)計(jì)算節(jié)點(diǎn)，炮域并行算法，30h高斯束偏移成像3D高斯束偏移高斯束偏移成像3D高斯束偏移高斯束偏移成像3D高斯束偏移Inline337常規(guī)Kirchhoff偏移高斯束偏移高斯束偏移成像3D高斯束偏移Depth=2080m常規(guī)Kirchhoff偏移高斯束偏移高斯束偏移成像3D高斯束偏移高斯束偏移成像3、高斯束偏移共成像道集1）偏移距域共成像點(diǎn)道集2）角度域共成像點(diǎn)道集二維反射角域共成像點(diǎn)道集三維方位角/反射角域共成像點(diǎn)道集多分量反射角域共成像點(diǎn)道集3）平面波域共成像點(diǎn)道集

映射轉(zhuǎn)換法空移成像條件時(shí)移成像條件L-ODCIGsTSCIGs高斯束偏移成像角道集的提取直接提取法走時(shí)拋物展開水平及垂直慢度傳播角度高斯束偏移成像角道集的提取傳播方向矢量3D反射張角3D反射方位角高斯束偏移成像角道集的提取P、S波傳播角度Snell定律P波反射張角S波反射張角高斯束偏移成像角道集的提取洼陷模型L-ODCIGsTSCIGs成像結(jié)果高斯束偏移成像角道集的提取空移成像條件時(shí)移成像條件直接提取高斯束偏移成像角道集的提取速度模型成像結(jié)果CDP=500CDP=850高斯束偏移成像角道集的提取28°21°49°15°22°38°Azimuth=0°Azimuth=135°ADCIGs高斯束偏移成像角道集的提取Inline337成像結(jié)果ADCIGs高斯束偏移成像角道集的提取勝利模型P-P成像結(jié)果P-S成像結(jié)果PSPSPSPSCDP=300CDP=500CDP=700高斯束偏移成像角道集的提取高斯束偏移成像3、高斯束偏移共成像道集1）偏移距域共成像點(diǎn)道集2）角度域共成像點(diǎn)道集二維反射角域共成像點(diǎn)道集三維方位角/反射角域共成像點(diǎn)道集多分量反射角域共成像點(diǎn)道集3）平面波域共成像點(diǎn)道集高斯束偏移平面波偏移高斯束偏移成像平面波共成像道集CDP=300CDP=400CDP=500CDP=600Marmousi成像結(jié)果高斯束偏移成像平面波共成像道集不同平面波成像結(jié)果-15~0°0~15°15~30°30~45°高斯束偏移成像平面波共成像道集主要內(nèi)容高斯束偏移成像1、聲波高斯束逆時(shí)偏移2、彈性波高斯束逆時(shí)偏移聲波高斯束逆時(shí)偏移高斯束偏移成像對(duì)于邊界，取地面水平面，邊界的側(cè)面忽略(如下圖)，Kirchhoff方程中的格林函數(shù)由高斯束漸近形式代替，其滿足Kirchhoff近似中地震邊界條件。結(jié)果有以下近似方程，即偏移域內(nèi)一點(diǎn)時(shí)間間隔為的反向延拓波場(chǎng)，地下一點(diǎn)，時(shí)間時(shí)的波場(chǎng)為：上行波場(chǎng)聲波高斯束逆時(shí)偏移高斯束偏移成像下行波場(chǎng)由地面一點(diǎn)源產(chǎn)生的直達(dá)波場(chǎng)