疊前地震數(shù)據(jù)特征波場(chǎng)成像與層析反演

疊前地震數(shù)據(jù)特征波場(chǎng)成像與層析反演

1局部平面波的成像與成像地震資料包含豐富有效的信號(hào)，但時(shí)間間隔內(nèi)的信號(hào)表達(dá)方法無(wú)法直接描述地震資料的特征。因此，可以通過(guò)特定的數(shù)學(xué)變換來(lái)提取地震信號(hào)的主要特征，進(jìn)行地震數(shù)據(jù)的稀疏表示，這有助于地震數(shù)據(jù)處理、波場(chǎng)分析和反演圖像。事實(shí)上，地震數(shù)據(jù)可以在轉(zhuǎn)換區(qū)域中稀疏表示。常用的轉(zhuǎn)換區(qū)域包括foyer、小波場(chǎng)、curvelet、rando、gaor等。一般來(lái)說(shuō)，轉(zhuǎn)換域中的地震數(shù)據(jù)可以是稀疏的，數(shù)據(jù)在轉(zhuǎn)換區(qū)域中的稀疏性可以是稀疏的。利用數(shù)據(jù)在轉(zhuǎn)換區(qū)域中的稀疏性，可以壓縮和重建地震數(shù)據(jù)，并具有插值的能力。xu等人（2005）反演不規(guī)則地震數(shù)據(jù)的頻譜，并進(jìn)行數(shù)據(jù)重建。將震源和卡盤(pán)上的地震數(shù)據(jù)稀疏度（2008）和白蘭淑（2014）用于地震數(shù)據(jù)重建。然而,Curvelet變換無(wú)法直觀描述地震波的傳播特征,因而沒(méi)有明確的地球物理含義.本文中,我們把疊前地震數(shù)據(jù)抽象為不同出射角度(即出射慢度矢量)以及不同到達(dá)時(shí)的局部平面波的線性疊加,而接收到的地震記錄可以看作局部平面波波前在空間離散的檢波器處的采樣.在立體層析理論(Billette和Lambaré,1998)中,地震記錄被簡(jiǎn)化為炮點(diǎn)、檢波點(diǎn)坐標(biāo),炮點(diǎn)、檢波點(diǎn)處的射線參數(shù)以及射線的雙程走時(shí).上述參數(shù)構(gòu)成了立體層析中的數(shù)據(jù)空間,但忽略了地震子波的波形信息.若用Gabor框架作為基函數(shù),則可以地震數(shù)據(jù)精確重構(gòu).ˇZáˇcek(2003)用相干態(tài)變換(coherentstatetransform)方法構(gòu)造一組完備的Gabor框架函數(shù),實(shí)現(xiàn)了地震數(shù)據(jù)的高斯波包分解.由于Gabor框架函數(shù)的冗余度高(Sondergaard,2007),該方法無(wú)法實(shí)現(xiàn)數(shù)據(jù)的稀疏表達(dá).Geng(2011)忽略了高斯波包中一些不重要的參數(shù),壓縮了重構(gòu)觀測(cè)波場(chǎng)的高斯波包泛函空間,從而減少框架函數(shù)的數(shù)目.李輝等(2014)提出了地震數(shù)據(jù)的特征高斯波包(CGP)分解方法,在Radon域選取出地震數(shù)據(jù)的時(shí)空和方向特征,大幅減少了高斯波包框架的數(shù)量,實(shí)現(xiàn)了Gabor域的稀疏化表達(dá).基于線性Radon變換(LRT)的局部平面波分解的方法是主要的地震信號(hào)的線性表達(dá)方法,以其明確的物理含義(可以描述局部平面波的傳播方向)而被廣泛應(yīng)用.非線性的Radon變換,如拋物Radon變換(Sacchi和Ulrych,1995),雙曲Radon變換(Liu和Sacchi,2002),多項(xiàng)式Radon變換(牛濱華等,2001)等,則主要用于多次波壓制、去噪、數(shù)據(jù)重構(gòu)等方面.不同于Fourier變換和正交小波變換,LRT的基函數(shù)是非完備正交的,因而導(dǎo)致LRT結(jié)果難以實(shí)現(xiàn)稀疏表達(dá).王雄文和王華忠(2014)將常規(guī)的LRT轉(zhuǎn)化為壓縮感知問(wèn)題,在壓縮感知的理論框架下實(shí)現(xiàn)了高分辨率的LRT.在Radon域中,可以測(cè)量局部平面波在炮點(diǎn)處的入射射線參數(shù)P本文中,首先考慮地震數(shù)據(jù)在局部平面波域中的稀疏表達(dá),并引入地震數(shù)據(jù)的特征分解方法,隨后給出基于特征波數(shù)據(jù)的偏移成像與反演成像方法.2局部平面波的成像在天然地震學(xué)中,通常用檢波器排列的某種組合來(lái)分離相干信號(hào)和噪聲.最基本的方法就是波束合成方法(BeamForming,Rost和Thomas,2002;Krüger等,1993),它利用局部平面波波前到達(dá)不同檢波器的時(shí)差,對(duì)檢波器接收到的地震記錄進(jìn)行時(shí)移疊加,使局部平面波相干疊加而其他信號(hào)相互抵消,從而提取某種特定的波型(震相).波束合成的思想也可應(yīng)用于勘探地震學(xué)中,如高斯束偏移方法(Hill,1990;Hale,1992),在炮集中對(duì)不同檢波器中心進(jìn)行加窗局部?jī)A斜疊加(即波束合成),利用炮集的τ-p譜進(jìn)行偏移.在射線束偏移(Sun等,2001)和控制束偏移(Gao,2006)中,通過(guò)選擇τ-p譜的能量,實(shí)現(xiàn)有選擇性的波場(chǎng)傳播與成像.波束合成算法需要估計(jì)局部平面波的方向信息.本文中,我們考慮Radon譜約束下的局部平面波反演方法(胡江濤等,2013):式中,X(ω)為局部空間窗內(nèi)的地震信號(hào),X(ω)=[x上述Radon譜約束下的局部平面波反演方法可以提高Radon譜的精度并壓制泄漏的噪聲.更進(jìn)一步地,為了得到更稀疏的局部平面波解(王雄文和王華忠,2014),可以把式(1)改寫(xiě)為壓縮感知問(wèn)題,在L0范數(shù)意義下求解如下反問(wèn)題:利用匹配追蹤方法求解(2)式,可以得到局部空間窗內(nèi)的稀疏平面波解.通過(guò)傳統(tǒng)的線性Radon變換或者高分辨率的平面波反演方法,可以得到局部平面波及其慢度矢量(或射線參數(shù)p).此時(shí),我們可以考慮對(duì)炮集的波束合成.在共炮點(diǎn)集中,對(duì)中心檢波點(diǎn)x其中,d(x我們把在中心檢波點(diǎn)處進(jìn)行一次波束合成之后的地震數(shù)據(jù)記為d若將地震數(shù)據(jù)分選為共檢波點(diǎn)道集,此時(shí)可以再次在中心炮點(diǎn)x其中,Φ(x同公式(3),經(jīng)過(guò)中心炮點(diǎn)處的再次波束合成后,實(shí)現(xiàn)了入射方向?yàn)閜將(3)式中的波束合成代入(4)式中,可得同時(shí)在中心炮點(diǎn)和檢波點(diǎn)處的波束合成算法:在三維時(shí),上述公式中的坐標(biāo)與射線參數(shù)均為矢量,表示為x從公式(5)可以看出,特征波場(chǎng)合成與局部平面波的入射與出射射線參數(shù)有關(guān),其本質(zhì)是利用局部平面波的入射和出射慢度計(jì)算道間時(shí)差(Δt=p·Δx),并將波束合成孔徑內(nèi)的地震道按照道間時(shí)差進(jìn)行時(shí)移疊加.在時(shí)移疊加的過(guò)程中,只有給定慢度矢量(局部平面波入射與出射射線參數(shù))的局部平面波波前實(shí)現(xiàn)了同相疊加,而其他的信號(hào)相互抵消,從而實(shí)現(xiàn)不同震相(波型)的分離,即特征波場(chǎng)合成.事實(shí)上,上述方法也可以推廣到任意線性排列的觀測(cè)系統(tǒng)中(如VSP或井間觀測(cè)系統(tǒng)等),只要滿(mǎn)足局部的炮點(diǎn)/檢波點(diǎn)排列具有線性特征.相應(yīng)的,局部平面波的慢度矢量(射線參數(shù))定義為地震波走時(shí)對(duì)地震排列(炮點(diǎn)或檢波點(diǎn)坐標(biāo))的方向?qū)?shù).此外,在頻率域特征波場(chǎng)合成公式(5)對(duì)應(yīng)的相移矩陣為A3特征波成像優(yōu)勢(shì)地震數(shù)據(jù)偏移成像的過(guò)程可以描述為地表波場(chǎng)的反傳播加上成像條件.以運(yùn)動(dòng)學(xué)射線追蹤為傳播算子的Kirchhoff疊前深度偏移(PSDM)技術(shù)在橫向變速劇烈情況下存在多路徑和振幅焦散等問(wèn)題,基于射線束的成像方法(劉少勇等,2012)在一定程度上解決了這些問(wèn)題.高斯束偏移技術(shù)(Hill,1990;Hale,1992)已經(jīng)成為工業(yè)界一個(gè)重要的疊前深度偏移成像工具.Sun和Schuster(2001)提出菲涅爾帶控制的波路徑偏移方法,將積分法偏移在成像空間中沿等時(shí)面的投影轉(zhuǎn)化為沿胖射線路徑的投影,在計(jì)算效率上有一定優(yōu)勢(shì).Liu和Palacharla(2011)對(duì)高斯束偏移進(jìn)行簡(jiǎn)化,只保留其運(yùn)動(dòng)學(xué)特征,推導(dǎo)出簡(jiǎn)單的解析表達(dá)式來(lái)控制射線束的寬度和振幅,理論上說(shuō)該實(shí)現(xiàn)方式效率要優(yōu)于高斯束偏移,成像精度基本有保障.快速射線束(Fast-beam)偏移(Gao等,2006)便是該思想的一種實(shí)現(xiàn).特征波偏移成像利用射線束波傳播算子來(lái)反傳局部特征波場(chǎng),既保留了射線類(lèi)成像的靈活性,又保證了成像精度.特征波場(chǎng)是地震波場(chǎng)在高維空間中的稀疏表達(dá),且局部平面波的入射與出射慢度矢量已知.利用這個(gè)性質(zhì),射線束類(lèi)的偏移方法可以實(shí)現(xiàn)從“畫(huà)弧”向“搬家”的轉(zhuǎn)變,進(jìn)一步提高成像效率.與上節(jié)的特征波場(chǎng)合成方法結(jié)合,匹配合適的局部射線束傳播算子,可以實(shí)現(xiàn)特征波疊前深度偏移成像(characteristicwaveimaging,簡(jiǎn)記為CWI).更進(jìn)一步的,“搬家”或一對(duì)一映射的偏移方法與基于角度道集的層析速度反演可以自然的融合成一體,共享同樣的數(shù)據(jù)空間和成像空間的數(shù)據(jù)體.成像與層析速度反演一體化也是特征波成像的一大優(yōu)勢(shì).傳統(tǒng)基于炮道集的射線束(包括高斯束)成像條件可以描述為其中:x表示成像點(diǎn)坐標(biāo),W其中,W傳統(tǒng)成像條件和特征波成像條件的幾何解釋可由圖1表示.圖中紅色反射層為地層真實(shí)位置,傳統(tǒng)射線束成像方法是將局部平面波投影到整個(gè)橢圓(三維情況為橢球)上.由于特征波場(chǎng)在地表的入射與出射慢度已知,因而特征波成像只對(duì)特征波場(chǎng)進(jìn)行“能量搬家”.紅色小橢圓的位置代表一對(duì)特征波射線束在成像域的菲涅爾帶范圍.從圖1可以看出,基于特征波場(chǎng)表達(dá),實(shí)現(xiàn)了成像方式從“畫(huà)弧”到“搬家”的轉(zhuǎn)變.傳統(tǒng)積分類(lèi)成像的畫(huà)弧成像過(guò)程必然會(huì)帶來(lái)畫(huà)弧噪聲,而搬家的方式可最大限度減少此類(lèi)噪聲同時(shí)能兼顧菲涅爾帶的成像分辨率,其成像效率也會(huì)得到很大的提高.4特征波場(chǎng)走時(shí)反演傳統(tǒng)的波動(dòng)方程走時(shí)反演方法(Luo和Schuster,1991)利用互相關(guān)方法估計(jì)觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的時(shí)差,需要對(duì)兩者進(jìn)行加窗處理,從而提取初至波或者折射波.該方法需要拾取特定的波型并在時(shí)間域加窗,繁重的拾取工作限制了該方法的應(yīng)用.利用特征波場(chǎng)合成,單道地震記錄被分離為一系列單獨(dú)的震相.此時(shí),可直接利用獨(dú)立的震相進(jìn)行反演,無(wú)需對(duì)地震記錄加窗處理.更重要的是,特征數(shù)據(jù)使得多尺度、逐層反演成為可能,如先利用折射波反演淺層速度,再利用由淺到深的反射波反演中深層速度.由于特征波場(chǎng)只含有單獨(dú)的波型(震相),可用互相關(guān)函數(shù)(Luo和Schuster,1991)測(cè)量觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的時(shí)差,無(wú)需對(duì)地震數(shù)據(jù)加窗.本文中,我們將特征波場(chǎng)合成算法應(yīng)用初至波場(chǎng)分解,并應(yīng)用于透射波走時(shí)反演中.至于特征波場(chǎng)分解在反射波走時(shí)反演中的應(yīng)用,可以參考馮波(2015),本文不再贅述.若觀測(cè)數(shù)據(jù)d與模擬數(shù)據(jù)u經(jīng)過(guò)特征波場(chǎng)合成后,得到的初至震相分別可以表示為d觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的初至走時(shí)殘差Δt定義為基于走時(shí)殘差L2范數(shù)的目標(biāo)函數(shù)可以表示為經(jīng)過(guò)推導(dǎo),目標(biāo)泛函梯度(馮波,2015)為其中,g對(duì)于初至走時(shí)反演,d梯度類(lèi)的迭代算法為:其中,α5單獨(dú)震相的成像反演根據(jù)上文的分析,我們可以利用特征波場(chǎng)合成實(shí)現(xiàn)不同波型(震相)的分離,而分離后的單獨(dú)震相可以用于后續(xù)的成像及反演.在數(shù)值實(shí)驗(yàn)中,我們首先分析不同信噪比的數(shù)據(jù)對(duì)射線參數(shù)估計(jì)以及特征波場(chǎng)合成的影響,然后,用一個(gè)簡(jiǎn)單模型(凹陷模型)展示特征波場(chǎng)成像的效果及其優(yōu)勢(shì).最后,用特征波場(chǎng)合成方法提取初至波場(chǎng),并用于波動(dòng)方程初至走時(shí)反演中.5.1特征波場(chǎng)成算法我們用Marmousi模型正演觀測(cè)數(shù)據(jù).速度模型為:X方向網(wǎng)格數(shù)目nx=737,Z方向網(wǎng)格數(shù)目nz=250;網(wǎng)格間隔均為12.5m.觀測(cè)系統(tǒng)設(shè)置為:第一炮的坐標(biāo)為625m,炮間隔25m,總共301炮.檢波器為固定排列,第一道的坐標(biāo)為625m,檢波器間隔12.5m,總共601道.用10Hz主頻的Ricker子波作為震源.首先,選取中心地震道的炮點(diǎn)和檢波點(diǎn)坐標(biāo)為(1850,3000).圖2a中,從左到右分別是x為了測(cè)量局部平面波的入射射線參數(shù),需要抽取共檢波點(diǎn)道集.圖2b中,從左到右分別是x利用圖2中拾取的入射和出射射線參數(shù),根據(jù)公式(5),可進(jìn)行特征波場(chǎng)合成.圖3a中,第1—4道為合成的特征波場(chǎng).其中,第1道代表對(duì)直達(dá)波的特征波場(chǎng)合成,它將直達(dá)波從地震記錄中提取.第2—4道代表對(duì)反射波的特征波場(chǎng)合成,它將反射波從地震記錄中提取.第5道為重構(gòu)的中心道(第1—4道疊加),與原始中心地震道(第6道)相比,只保留了我們提取的震相,提高了特征波的信噪比.同時(shí),考慮到特征波場(chǎng)(第1—4道)的時(shí)間局部特征,可以(在變換域)進(jìn)一步壓縮存儲(chǔ).接著,為了測(cè)試噪聲對(duì)特征波場(chǎng)合成的影響,我們對(duì)原始數(shù)據(jù)加入了高斯白噪聲(S/N=100).圖4中展示了相應(yīng)的道集及其τ-p譜(所有參數(shù)含義同圖2).可以看出,高斯隨機(jī)噪聲(S/N=100)降低了τ-p譜的分辨率,因而可能會(huì)對(duì)τ-p譜的自動(dòng)拾取有一定的影響,但對(duì)人工拾取的影響較小.圖3b中,展示了對(duì)圖4中的數(shù)據(jù)進(jìn)行特征波場(chǎng)合成的結(jié)果(所有參數(shù)同圖2a).顯然,特征波場(chǎng)合成之后,仍然能夠得到高信噪比的單獨(dú)震相.這是由于特征波場(chǎng)合成的過(guò)程中,同相的信號(hào)相互增強(qiáng),而其他信號(hào)相互干涉.若進(jìn)一步增強(qiáng)噪聲強(qiáng)度(S/N=20),此時(shí),隨機(jī)噪聲完全壓制了有效信號(hào),即使在τ-p譜也無(wú)法拾取射線參數(shù).若我們?nèi)匀挥谜鎸?shí)的射線參數(shù)(圖2中拾取的結(jié)果),特征波場(chǎng)合成結(jié)果如圖4c所示.其中,直達(dá)波的波形仍然具有較為保真的形態(tài),而后續(xù)的幾個(gè)反射波則受到了一定程度的噪聲干擾,略有震蕩.上述試驗(yàn)說(shuō)明,本文中的特征波場(chǎng)合成算法十分穩(wěn)健,隨機(jī)噪聲對(duì)其影響較小.當(dāng)噪聲能量逐漸增強(qiáng)時(shí),τ-p譜的分辨率會(huì)逐漸降低,直至無(wú)法拾取準(zhǔn)確的射線參數(shù),進(jìn)而導(dǎo)致特征波場(chǎng)合成算法失效.5.2特征波場(chǎng)成像優(yōu)勢(shì)本文用洼陷模型測(cè)試特征波場(chǎng)成像方法的有效性,其速度模型如圖5a所示.模擬數(shù)據(jù)共201炮,炮間距20m,每一炮301道接收,道間距10m.本節(jié)分別對(duì)該模型數(shù)據(jù)進(jìn)行傳統(tǒng)的KirchhoffPSDM和特征波成像來(lái)進(jìn)行效率效果的對(duì)比.圖5b是KirchhoffPSDM的成像剖面,由于其成像過(guò)程是通過(guò)等時(shí)面疊加收斂繞射波,成像剖面上留下了畫(huà)弧噪聲.圖5c是特征波場(chǎng)成像得到的成像剖面,由于其通過(guò)對(duì)特征域數(shù)據(jù)進(jìn)行“搬家”,其畫(huà)弧范圍控制在局部范圍內(nèi),其畫(huà)弧噪聲明顯弱于傳統(tǒng)KirchhoffPSDM.由于特征波的成像過(guò)程中可以得到射線方向,張角道集可以方便的輸出.傳統(tǒng)KirchhoffPSDM一般只能輸出地表偏移距道集,在對(duì)后續(xù)的速度分析或?qū)游龀上竦倪m應(yīng)性方面上來(lái)說(shuō),張角道集質(zhì)量要優(yōu)于地表偏移距道集.圖5(d,e)分別為洼陷模型3500m處的KirchhoffPSDM地表偏移距道集和特征波成像的張角道集,對(duì)比兩圖可見(jiàn),張角道集的成像質(zhì)量要高于地表偏移距道集.特征波成像由于只對(duì)特征數(shù)據(jù)進(jìn)行反投影(特征數(shù)據(jù)是對(duì)原地震數(shù)據(jù)的一個(gè)壓縮處理),其成像效率要遠(yuǎn)遠(yuǎn)高于傳統(tǒng)的KirchhoffPSDM或者傳統(tǒng)的射線束成像.本文對(duì)比了KirchhoffPSDM和特征波成像在洼陷模型上的成像效率,具體數(shù)據(jù)見(jiàn)表1,數(shù)據(jù)得到了有效的壓縮,同時(shí)成像效率也有近一個(gè)數(shù)量級(jí)的提高.由于特征波場(chǎng)包含了數(shù)據(jù)方向信息,特征波場(chǎng)成像方法可以方便的推廣到目標(biāo)體成像,圖5f展示了僅僅對(duì)最后一個(gè)反射層的面向目標(biāo)成像結(jié)果.對(duì)特定目標(biāo)成像將會(huì)進(jìn)一步提高效率,縮短處理周期.5.3初至走時(shí)反演速度模型在本節(jié)中,我們將特征波場(chǎng)合成算法分離初至波型,并用于波動(dòng)方程初至走時(shí)層析中.速度模型(BP2004速度模型的右端部分,橫向約26km,深度約12km)如圖6a所示.觀測(cè)系統(tǒng)為:炮點(diǎn)深度12.5m,炮點(diǎn)間隔為200m,共107炮.檢波器深度12.5m,間隔400m,每炮54道.觀測(cè)系統(tǒng)總共5778道.初始速度模型采用一維模型,其中,水體速度以及水底深度已知.總共進(jìn)行18次反演迭代,反演得到的速度模型見(jiàn)圖6b.反演得到了整體背景速度的變化趨勢(shì),其中模型右側(cè)的高速隆起反演得較好,但沒(méi)有反演出淺層小尺度的低速異常.圖7中為初至波走時(shí)對(duì)比曲線.其中,紅色、綠色和藍(lán)色曲線分別是真實(shí)模型、初始模型和第18次反演的速度模型中的初至走時(shí).可以看出,真實(shí)的初至走時(shí)與反演得到的速度模型中的初至走時(shí)基本吻合,實(shí)現(xiàn)了初至走時(shí)反演的目標(biāo).6特征波場(chǎng)輔助偏移成像本文提出了一種地震數(shù)據(jù)的稀疏表達(dá)方法.不同于傳統(tǒng)的變換域中的數(shù)據(jù)稀疏表達(dá)理論,該方法利用局部平面波的傳播方向(慢度矢量),在中心炮檢點(diǎn)處同時(shí)進(jìn)行波束合成,從而將地震數(shù)據(jù)投影到局部平面波域(高維空間,二維觀測(cè)對(duì)應(yīng)5D數(shù)據(jù);三維觀測(cè)對(duì)應(yīng)9D數(shù)據(jù))中.由于波束合成后的地震數(shù)據(jù)描述了局部平面波的方向特征,因此我們稱(chēng)之為特征波場(chǎng).相應(yīng)的,我們稱(chēng)這種波場(chǎng)合成(分離)算法為特征波場(chǎng)合成方法.特征