地震正演與高斯束偏移

地震成像與速度反演李振春中國石油大學（華東）地球科學與技術(shù)學院提綱一、雙變網(wǎ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ì)波動問題的有效方法！雙變網(wǎng)格正演模擬原理

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

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

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

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

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

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

已知時刻的波場值已知時刻的波場值

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

已知時刻的波場值已知時刻的波場值

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

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

波場快照對比圖

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

炮記錄對比圖

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

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

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

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

加窗局部傾斜疊加并計算由震源和接收點出射的高斯束對于每個震源射線參數(shù){

對于每個接收點射線參數(shù){

對于局部孔徑內(nèi)的每個成像點{

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

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

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

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

加窗局部傾斜疊加并計算由震源和接收點出射的高斯束對于每個中心點射線參數(shù){

對于每個粗網(wǎng)格點{

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

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

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

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

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

偏移過程中消除地震波幾何擴散造成的能量損失，使得成像振幅正比于界面的反射系數(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）復雜地表高斯束偏移4）彈性波多分量高斯束偏移5）3D高斯束偏移復雜地表條件下地震成像方法波動方程基準面校正靜校正法波場沿拓法射線類偏移方法準確性差計算量大，靈活性差，對數(shù)據(jù)要求高靈活，高效復雜地表高斯束偏移高斯束偏移成像高斯束偏移成像局部靜校正法基本原理（Gray，Geophysics，2005）復雜地表高斯束偏移RayleighII積分（岳玉波等，地球物理學報，2011第三期）復雜地表保幅成像公式起伏地表平面波分解

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

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

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