崔永福， 彭更新， 吳國(guó)忱， 尚帥， 郭念民， 趙銳銳

1 中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 山東 青島　266580　2 中國(guó)石油塔里木油田分公司勘探開(kāi)發(fā)研究院， 新疆 庫(kù)爾勒　841000

?

全波形反演在縫洞型儲(chǔ)層速度建模中的應(yīng)用

崔永福1,2， 彭更新2， 吳國(guó)忱1， 尚帥2， 郭念民2， 趙銳銳2

1 中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 山東 青島2665802 中國(guó)石油塔里木油田分公司勘探開(kāi)發(fā)研究院， 新疆 庫(kù)爾勒841000

摘要速度是地震偏移成像準(zhǔn)確與否的關(guān)鍵所在.全波形反演綜合利用地震波場(chǎng)運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)信息，能夠得到相比傳統(tǒng)速度建模方法更高頻的成分.全波形反演的理論比較成熟，但實(shí)際應(yīng)用成功的例子相對(duì)較少，特別是對(duì)于陸上地震資料.塔里木盆地地震地質(zhì)條件復(fù)雜，為了實(shí)現(xiàn)縫洞型儲(chǔ)層的準(zhǔn)確成像，本文開(kāi)展了針對(duì)目標(biāo)靶區(qū)的全波形反演精細(xì)速度建場(chǎng)研究.采用一種時(shí)間域分層多尺度全波形反演流程：首先通過(guò)層析成像建立初始速度模型；其次利用折射波反演淺層速度模型；最后利用反射波反演中深層速度模型.偏移成像結(jié)果表明基于全波形反演的速度建模技術(shù)能有效改善火成巖下伏構(gòu)造的成像精度，顯示了全波形反演在常規(guī)陸上采集資料的應(yīng)用潛力.

關(guān)鍵詞全波形反演； 陸上數(shù)據(jù)； 速度建場(chǎng)； 縫洞型儲(chǔ)層； 火成巖

CUI Yong-Fu1,2， PENG Geng-Xin2， WU Guo-Chen1， SHANG Shuai2，

GUO Nian-Min2， ZHAO Rui-Rui2

1SchoolofGeosciences，ChinaUniversityofPetroleum(Huadong)，QingdaoShandong266580,China2ResearchInstituteofExplorationandDevelopment,TarimOilfieldCompany,PetroChina,KorlaXinjiang841000,China

By detailed analysis of FWI theory and its application challenges for land seismic data, taking into account of the actual situation of seismic data in YM area, we propose special strategies for successful application of FWI in this area: (1) Joint denoising technology based on surface wave modeling and curvelet transform is used in preprocessing stage to retain low frequency effective signal as far as possible, the five-dimensional regularization method is also used; (2) we utilize the high frequency component of static correction to eliminate the high frequency of near-surface velocity model and then use first arrival and reflected wave tomography to build a more accurate initial velocity model; (3) With the initial velocity model provided by traveltime tomography, the hierarchical multi-scale FWI method in time domain is applied. We utilize refraction waveform information to invert for shallow velocity while reflection waveform for deep velocity model.

The actual data processing result shows that the joint denoising technology based on surface wave modeling and curvelet transform can protect low frequency information more effectively, maintaining the dynamics and kinematics information and meeting the basic requirement of FWI method. The initial velocity model provided by tomography can satisfy the accuracy of FWI. The hierarchical multi-scale FWI in time domain can characterize the velocity of igneous rock with high accuracy. The result of pre-stack depth migration indicates that the imaging of target formation under the igneous rock is improved obviously, eliminating the “fault” phenomenon caused by rough velocity model, and the imaging of fractured-cavernous reservoir is better.

Although FWI is a high-precision velocity modeling technology with perfect theory, its application for land seismic data is still a challenge. The accuracy of initial velocity model and useful low frequency seismic information are keys to affect its result. We believe that the productive application of FWI for land seismic data will be developed gradually with the development of acquisition and processing technology.

1引言

速度是描述地下介質(zhì)構(gòu)造和儲(chǔ)層特性的重要物性參數(shù)，準(zhǔn)確的速度模型對(duì)于偏移至關(guān)重要.目前常用的速度建模方法主要包括疊加速度分析、偏移速度分析和層析成像等(符力耘等，2013).而隨著勘探目標(biāo)的日益復(fù)雜及對(duì)成像精度要求越來(lái)越高，傳統(tǒng)的速度建模方法已經(jīng)難以滿足地震資料處理與解釋的需求，而全波形反演(Full-waveform inversion，F(xiàn)WI)可以有效地彌補(bǔ)其不足.全波形反演綜合利用疊前地震波場(chǎng)的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)信息，通過(guò)不斷匹配模擬記錄與實(shí)測(cè)資料來(lái)更新速度模型，最終獲取可以準(zhǔn)確描述地下介質(zhì)速度分布的深度域速度場(chǎng)模型，具有較高的反演精度，開(kāi)始越來(lái)越多被關(guān)注并應(yīng)用在實(shí)際資料中.全波形反演利用地震記錄中的全部有效信息，通過(guò)局部尋優(yōu)來(lái)進(jìn)行初始模型的迭代，逐步逼近真實(shí)模型，因此本質(zhì)上是一個(gè)極度“病態(tài)”的非線性優(yōu)化問(wèn)題.地震波場(chǎng)模擬精度、計(jì)算效率和地震資料低頻成分的保留情況是影響全波形反演方法實(shí)際資料應(yīng)用的重要因素.

20世紀(jì)80年代，Tarantola(1984)首先提出了基于廣義最小二乘的時(shí)間域全波形反演思想，通過(guò)對(duì)正演波場(chǎng)和實(shí)測(cè)記錄偏差的最小化約束來(lái)更新速度模型.此后國(guó)內(nèi)外研究學(xué)者分別從正演和反演兩個(gè)方向?qū)ζ溥M(jìn)行了發(fā)展.在正演模擬方面，Virieux(1984)將二階交錯(cuò)網(wǎng)格有限差分方法運(yùn)用于地震波場(chǎng)模擬；之后，Levander(1988)將該方法推廣至四階.現(xiàn)在高階交錯(cuò)網(wǎng)格有限差分法廣泛應(yīng)用于全波形反演中，此外還有偽譜法(Danecek and Seriani, 2008)、有限元法(Marfurt，1984)等數(shù)值模擬方法.在反演方面，考慮時(shí)間域全波形反演可以對(duì)數(shù)據(jù)進(jìn)行靈活處理，Bunks等(1995)提出了時(shí)間域的多尺度反演方法，通過(guò)將問(wèn)題分解成不同的頻率尺度，降低其非線性程度；為了進(jìn)一步提高計(jì)算效率，Pratt和Goulty(1991)將Tarantola的理論從時(shí)間域發(fā)展到了頻率域，開(kāi)啟了頻率域全波形反演的研究熱潮；Sirgue和Pratt(2004)提出基于頻率迭代的多尺度全波形反演方法并且給出了反演頻率的選擇標(biāo)準(zhǔn).頻率域波形反演直接在頻率域進(jìn)行求解，容易實(shí)現(xiàn)從低頻到高頻的多尺度反演，僅需要幾個(gè)離散的頻率便可以完成速度模型的重建，但Operto等(2007)指出，當(dāng)處理三維較大模型時(shí)，由于三維頻率域正演對(duì)內(nèi)存的超大需求，使得三維頻率域波形反演有較大的局限性.國(guó)內(nèi)的一些學(xué)者(胡光輝等，2013；董烈乾等，2011；王薇等，2013；魏哲楓等，2014；曹書紅和陳景波，2014；劉玉柱等，2014；楊勤勇等，2014)針對(duì)全波形反演理論研究同樣取得了一些成果.

與理論研究相比，全波形反演的實(shí)際應(yīng)用起步較晚，主要原因是全波形反演方法本身的一些局限以及對(duì)數(shù)據(jù)的要求過(guò)高.Gauthier等(1986)和Mora(1987)于20世紀(jì)80年代率先實(shí)現(xiàn)了二維地震資料的全波形反演，證明了全波形反演的高精度建模的潛力，但同時(shí)也指出全波形反演的“病態(tài)”性使其在地震資料缺少低頻情況下很難成功；20世紀(jì)90年代，全波形反演主要被用于井間地震數(shù)據(jù)，利用寬角透射波信息重建速度模型(Pratt and Goulty，1991；Song et al.，1995)；Sirgue等(2010)率先對(duì)挪威北海油田OBC數(shù)據(jù)實(shí)現(xiàn)了三維全波形反演，對(duì)該地區(qū)淺層氣云和周邊充氣的斷裂構(gòu)造進(jìn)行了精細(xì)刻畫，極大地鼓舞了全波形反演的研究熱潮，之后全波形反演在海上地震資料的應(yīng)用越來(lái)越多(Prieux et al.，2010；Ratcliffe et al.，2011；Takougang and Calvert，2011；韓淼，2014).相比海上地震資料全波形反演的應(yīng)用，陸地資料由于缺少足夠多的低頻成分和初始模型的限制，使得陸上全波形反演方法的實(shí)現(xiàn)較難.2012年，殼牌公司和東方物探合作實(shí)現(xiàn)了二維陸上資料的全波形反演，其數(shù)據(jù)采集基于低頻大偏移距的特殊觀測(cè)系統(tǒng)，但在實(shí)際地震采集中很難推廣，因此不具有普遍性，但至少驗(yàn)證了陸地全波形反演應(yīng)用的可行性(Plessix et al.，2012).

塔里木盆地哈拉哈塘地區(qū)碳酸鹽巖儲(chǔ)層普遍埋藏深(超過(guò)5500 m)、非均質(zhì)性強(qiáng)，目的層上覆地層發(fā)育厚度不均、速度不同的多期火成巖高速層，該火成巖對(duì)速度建模帶來(lái)極大挑戰(zhàn)，火成巖速度的準(zhǔn)確與否直接影響縫洞型儲(chǔ)層偏移成像位置的準(zhǔn)度.常規(guī)基于沿層層析和網(wǎng)格層析建模技術(shù)刻畫火成巖速度精度不高，難以消除由于火成巖速度不準(zhǔn)確引起的成像假象.因此，針對(duì)該區(qū)地質(zhì)條件，地震資料噪音類型、近地表、能量不均衡等問(wèn)題，建立了面向全波形反演的疊前保幅去噪、靜校正等預(yù)處理技術(shù)流程，提出了適合該區(qū)的全波形反演速度建場(chǎng)技術(shù)流程.實(shí)際資料全波形反演結(jié)果證實(shí)了基于波形反演的速度建模方法相比常規(guī)層析成像速度建場(chǎng)精度更具優(yōu)勢(shì)，其精細(xì)刻畫了火成巖的構(gòu)造形態(tài)和速度分布，有效改善了火成巖下伏地層的成像精度.

2陸上資料全波形反演應(yīng)用挑戰(zhàn)及對(duì)策

全波形反演理論上很完美，但它對(duì)數(shù)據(jù)質(zhì)量、初始速度模型、正演子波等是有理論假設(shè)的.目前，F(xiàn)WI反演很難收斂到正確的結(jié)果上，其核心問(wèn)題在于誤差泛函存在非常多的局部極值點(diǎn)，因此需要很好的初始模型來(lái)降低誤差泛函的非線性.

與海上資料相比，陸上地震資料應(yīng)用波形反演的主要挑戰(zhàn)在于：

(1) 要求寬方位、大偏移距的觀測(cè)系統(tǒng).常規(guī)陸上三維多為滾動(dòng)采集，縱橫比和偏移距較小，這在很大程度上限制了全波場(chǎng)信息的有效獲??； (2) 近地表結(jié)構(gòu)和地下構(gòu)造的雙重復(fù)雜導(dǎo)致信噪比低，簡(jiǎn)單的正演模擬難以描述實(shí)際資料中的復(fù)雜波；同時(shí)地震子波空間不一致性加重了地震波場(chǎng)與反演參數(shù)之間的非線性； (3) 缺失低頻.海上資料一般最低有效頻帶為3或者3.5 Hz，常規(guī)陸上地震資料的頻帶有效范圍在6 Hz以上，因此，沒(méi)有低頻數(shù)據(jù)與初始模型的很好耦合，容易陷入局部極小解； (4) 數(shù)據(jù)預(yù)處理難度大.全波形反演要運(yùn)用地震波的走時(shí)信息和動(dòng)力學(xué)信息，在數(shù)據(jù)預(yù)處理時(shí)，既要消除噪聲干擾又不能破壞波的動(dòng)力學(xué)特征; (5) 陸上炸藥震源穩(wěn)定性較差，隨機(jī)干擾嚴(yán)重，地面檢波器與地面的耦合性遠(yuǎn)遠(yuǎn)不如水上檢波器與海水的一致性耦合；同時(shí)，陸上近地表介質(zhì)極其復(fù)雜，許多不確定因素導(dǎo)致炮、檢點(diǎn)分布不規(guī)則或者有空洞，炮、檢點(diǎn)的高程、坐標(biāo)不準(zhǔn)等誤差在構(gòu)建誤差函數(shù)時(shí)累積，影響其收斂.

在深入分析FWI理論及其陸上資料應(yīng)用挑戰(zhàn)后，提出了實(shí)現(xiàn)陸上波形反演實(shí)用化的針對(duì)性對(duì)策.

(1) 優(yōu)選盡量滿足其基本假設(shè)的地震數(shù)據(jù).

本文研究的YM工區(qū)位于塔北隆起南斜坡帶，地面海拔在900～1100 m，地形較為平坦，近地表簡(jiǎn)單；震源子波空間一致性好；低頻信息較豐富、低頻信息低至3 Hz；縱向最大炮檢距為6375 m、最大非縱距4775 m、方位較寬(橫縱比達(dá)到0.75)，最大炮檢距為7965；本區(qū)二疊系除了發(fā)育多期火成巖外，地下構(gòu)造簡(jiǎn)單，相對(duì)于目的層埋深超過(guò)6000 m，偏移距偏小.

(2) 設(shè)計(jì)保真的數(shù)據(jù)預(yù)處理流程.

陸上資料普遍存在面波、散射面波、環(huán)境干擾等，能量高、多在低頻段，常規(guī)去噪技術(shù)往往不同程度地傷害低頻有效信號(hào)，設(shè)計(jì)保持動(dòng)力學(xué)和運(yùn)動(dòng)學(xué)特征的數(shù)據(jù)預(yù)處理流程是關(guān)鍵.

(3) 建立高精度的初始速度模型.

由于該區(qū)偏移距偏小，地質(zhì)體的照明度不足，導(dǎo)致全局收斂慢，需要建立高精度的初始速度模型減少迭代次數(shù)和避免局部最優(yōu)解.

通常陸上資料近地表速度低且變化劇烈，炮、檢點(diǎn)的位置和高程也可能不準(zhǔn)，現(xiàn)有方法難以描述速度的高頻變化，處理中采用靜校正的高頻分量來(lái)消除近地表速度高頻，近地表速度的中、長(zhǎng)波長(zhǎng)分量放在速度建模中解決；在該區(qū)采用基于層析成像的速度建模方法建立的初始速度模型，已經(jīng)基本滿足該區(qū)的疊前深度偏移成像對(duì)速度精度要求.

隨著陸上地震采集技術(shù)水平的不斷提高，全波場(chǎng)、寬頻帶等觀測(cè)手段的實(shí)現(xiàn)以及保真預(yù)處理技術(shù)的發(fā)展，陸上全波形反演的應(yīng)用將逐步走向?qū)嵱没?

3時(shí)間域全波形反演原理及流程實(shí)現(xiàn)

全波形反演方法的基本原理是通過(guò)給定一個(gè)初始模型，利用正演模擬得到其傳播波場(chǎng)，將正演波場(chǎng)與實(shí)測(cè)波場(chǎng)進(jìn)行對(duì)比分析，通過(guò)不斷修正模型使得二者誤差達(dá)到可允許范圍，從而得到反演結(jié)果.該方法主要包括三部分(韓淼，2014)：首先，基于已知的初始速度模型，通過(guò)波動(dòng)方程正演模擬獲得理論地震波場(chǎng)；其次，基于理論地震波場(chǎng)與觀測(cè)地震波場(chǎng)之間的波場(chǎng)殘差，構(gòu)建目標(biāo)函數(shù)，并將殘差逆時(shí)傳播得到反向波場(chǎng)，并通過(guò)正向波場(chǎng)和反向波場(chǎng)計(jì)算模型參數(shù)梯度；最后，選擇一種優(yōu)化算法對(duì)該反問(wèn)題進(jìn)行求解.全波形反演方法可以分為時(shí)間域全波形反演和頻率域全波形反演.為更好適應(yīng)三維，本文采用時(shí)間域全波形反演方法，為了增加問(wèn)題求解的穩(wěn)定性降低其非線性程度，保證求解不會(huì)陷入局部收斂應(yīng)用中采取多尺度反演(Bunks et al.，1995).

3.1時(shí)間域全波形反演基本原理

時(shí)間域正演通過(guò)在時(shí)間方向的迭代來(lái)求各個(gè)時(shí)間的波場(chǎng)值，內(nèi)存需求相對(duì)較小，能夠更好地適應(yīng)三維問(wèn)題.

時(shí)間域的三維聲波方程可表示為(Vigh et al.，2011)

(1)

方程中P(x,y,z,t)是地震波場(chǎng)，V(x,y,z)是速度項(xiàng).假設(shè)P(xr,yr,zr,t)是xr處記錄的地震波場(chǎng)，構(gòu)建如下所示目標(biāo)函數(shù)：

(2)

公式中Pobs是野外觀測(cè)到的地震數(shù)據(jù)，Pcal是通過(guò)聲波波動(dòng)方程正演的記錄.全波形反演通過(guò)尋找目標(biāo)函數(shù)E的極小值，使得滿足目標(biāo)函數(shù)極小值的速度即為最終的反演速度.

3.2波形反演速度建模實(shí)現(xiàn)流程

陸上資料的全波形反演應(yīng)用主要受限于子波空變、信噪比低、低頻信息缺失、數(shù)據(jù)不規(guī)則、低頻數(shù)據(jù)與初始模型的耦合.針對(duì)FWI在陸上資料應(yīng)用受限原因及哈拉哈塘地震資料特點(diǎn)，提出地震保幅預(yù)處理流程和實(shí)現(xiàn)FWI速度建模流程，流程如圖1所示.

圖1　保幅預(yù)處理流程(a)與實(shí)現(xiàn)FWI速度建模流程(b)Fig.1　Fidelity preprocessing strategy (a) and FWI velocity modeling strategy (b)

4預(yù)處理關(guān)鍵技術(shù)及效果

YM地區(qū)地表高程變化不大，北部為農(nóng)田、中部為浮土、南部為小沙，低降速帶厚度變化不大，靜校正難度不大；工區(qū)內(nèi)疏松的近地表低降速帶對(duì)地震波吸收衰減嚴(yán)重，導(dǎo)致深層有效波能量弱、主頻較低；主要干擾波為強(qiáng)能量的面波、折射波及異常干擾；奧陶系縫洞型碳酸鹽巖是主要目的層，埋藏深、非均質(zhì)強(qiáng)，因此噪聲衰減、恢復(fù)弱信號(hào)和保護(hù)低頻信號(hào)是資料預(yù)處理的關(guān)鍵技術(shù).

4.1保幅去噪技術(shù)

針對(duì)該區(qū)噪音類型，本文首先采用面波模擬方法(Strobbia et al.，2009)去除高能量頻散面波；然后利用穩(wěn)健的地表一致性反褶積(Tnaer and Koehler，1981)壓制淺層多次折射波；再利用曲波變換(董烈乾等，2011)去除剩余面波；最后進(jìn)行異常振幅消除處理，預(yù)處理流程見(jiàn)圖1a所示.

圖2a為靜校正后單炮數(shù)據(jù)，圖2b為面波模擬法壓制模擬的基階、一階和二階面波噪聲后的結(jié)果.衰減剩余更高階面波和頻散面波之前，先對(duì)數(shù)據(jù)壓制淺層多次折射波，如圖2c為淺層多次折射波壓制后的炮集.圖2d為曲波變換法衰減剩余面波之后的炮集，可以看出高階面波和頻散面波衰減效果好，大大提高了資料信噪比.圖3、圖4為本文方法與常規(guī)三維FKK方法去噪后低頻信號(hào)(0～10 Hz)和頻譜圖對(duì)比，可以看出本文方法可以保留更多的有效低頻信息，更有利于后續(xù)疊前反演、流體檢測(cè)、波形反演等技術(shù)應(yīng)用.

圖2　噪聲壓制效果圖 (a) 靜校正后單炮； (b) 面波模擬法去除噪聲后單炮； (c) 壓制淺層多次折射后單炮； (d) 壓制剩余面波后單炮.Fig.2　(a) Data with field static correction； (b) Noise reduction with surface wave modeling method； (c) Data after surface consistent deconvolution； (d) Residual surface wave reduction

圖3　常規(guī)FKK去噪(a)與本文方法去噪(b)PSTM低頻部分對(duì)比(0～10 Hz)Fig.3　Low-frequency component (0～10 Hz) of seismic data after FKK denoising (a) and our method (b)

圖4　常規(guī)三維FKK去噪(藍(lán))與本文方法去噪后(紅)頻譜對(duì)比Fig.4　Spectrum of seismic data after FKK denoising (blue) and our method (right)

4.2五維數(shù)據(jù)規(guī)則化與偏移面校正

YM地區(qū)地表簡(jiǎn)單，高程變化不大、地表障礙極少，不存在嚴(yán)重的數(shù)據(jù)不規(guī)則或較大面積的空洞，僅存在人為施工誤差造成原始地震資料的炮點(diǎn)、檢波點(diǎn)分布不規(guī)則或壞道造成部分?jǐn)?shù)據(jù)缺失，給后續(xù)正演模擬與實(shí)際地震資料匹配帶來(lái)一定困難.本文采用基于OMP(正交匹配追蹤)算法的五維數(shù)據(jù)規(guī)則化技術(shù)(Hollander et al.，2012)極大地提高了重建質(zhì)量，解決了正演模擬與實(shí)際資料的匹配問(wèn)題.

前期靜校正處理已經(jīng)對(duì)地表高程變化和低降速帶異常時(shí)差進(jìn)行了校正，這對(duì)地震波場(chǎng)時(shí)間信息的破壞非常大，此時(shí)以CMP面或者固定基準(zhǔn)面作為偏移的基準(zhǔn)面，往往不能正確求取近地表速度，甚至誤差很大.理論上講，應(yīng)選取真實(shí)地表面，但是現(xiàn)有資料無(wú)法實(shí)現(xiàn)，目前最現(xiàn)實(shí)的情況是選取平滑地表作為偏移的基準(zhǔn)面，把數(shù)據(jù)從CMP面校正到近地表平滑面，同時(shí)數(shù)據(jù)中僅僅保留地表平滑面附近因高程變化和存在微小異常體產(chǎn)生的高頻靜校正時(shí)差，而去掉因低降速帶異常體引起的中長(zhǎng)波長(zhǎng)靜校正，此部分的速度異常在深度偏移淺層速度建模中考慮，這是陸上地震資料預(yù)處理的關(guān)鍵環(huán)節(jié)之一.

5波形反演精細(xì)速度建模

全波形反演理論體系完美，但它對(duì)數(shù)據(jù)限制和應(yīng)用瓶頸很多.陸上資料的全波形反演應(yīng)用主要受限于子波空變、信噪比低、低頻觀測(cè)數(shù)據(jù)和滿足全波形反演要求的初始速度模型，其中低頻數(shù)據(jù)與初始模型的耦合是全波形反演在實(shí)際資料應(yīng)用中遇到的最大瓶頸.

YM工區(qū)地表?xiàng)l件相對(duì)簡(jiǎn)單，子波空間一致性較好，信噪比較高，野外采集沒(méi)有低截處理；在0～3 Hz低頻端地震信號(hào)中折射波的能量較強(qiáng)，以折射波信號(hào)為主，反射信號(hào)能量弱，在在0～5 Hz低頻端地震信號(hào)中反射波能量強(qiáng)，信號(hào)低頻較豐富(如圖6所示).上述數(shù)據(jù)預(yù)處理最大程度地保護(hù)了數(shù)據(jù)低頻成分，使該區(qū)地震資料滿足波形反演基本要求.

圖5　近地表及炮、檢點(diǎn)布設(shè)示意圖Fig.5　Diagram of near-surface and seismic geometry

圖6　原始數(shù)據(jù)低頻濾波結(jié)果對(duì)比分析 (a) 原始單炮記錄； (b) 0～3 Hz低通濾波后單炮記錄； (c) 0～5 Hz低通濾波后單炮記錄； (d) 0～7 Hz低通濾波后單炮記錄.Fig.6　Low frequency pass filter of raw shot data (a) Raw shot record； (b) Low pass filter of 0～3 Hz； (c) Low pass filter of 0～5 Hz； (d) Low pass filter of 0～7 Hz.

在實(shí)際資料處理中，高精度正演模擬和全波形反演的計(jì)算量非常大，所以盡可能通過(guò)建立高精度的初始速度模型，減少迭代次數(shù)，增加求解的穩(wěn)定性，提高收斂速度，從而減少計(jì)算量.

5.1初始速度模型的建立

初始速度模型的建立主要利用走時(shí)層析成像技術(shù)，具體流程如圖1b.先采用回轉(zhuǎn)波層析成像技術(shù)建立淺層速度模型；然后利用反射波層析成像建立中深層各向同性速度模型；再利用井?dāng)?shù)據(jù)標(biāo)定各向異性參數(shù)，建立中深層各向異性速度模型(如圖7a所示)；最后，對(duì)該各向異性速度模型做平滑處理，去掉速度高頻成分的背景速度作為波形反演的最終初始速度模型(如圖7b所示).

圖7　層析成像建立的各向異性速度模型(a)與FWI初始速度模型(b)Fig.7　Anisotropy velocity model generated by tomography (a) and initial velocity of FWI (b)

初始速度模型建立以后，采用時(shí)間域分層多尺度全波形反演精細(xì)建模流程，主要分為以下兩步：(1)淺層折射波FWI，(2)中深層反射波FWI.最后利用疊前深度偏移對(duì)地下構(gòu)造進(jìn)行偏移成像.

5.2關(guān)鍵技術(shù)細(xì)節(jié)

波形反演成功與否，與波形反演的很多技術(shù)細(xì)節(jié)密切相關(guān).結(jié)合本工區(qū)實(shí)際資料情況，具體技術(shù)細(xì)節(jié)如下：

① 地震正演子波，指導(dǎo)原則是盡可能接近實(shí)際資料地震子波.本文從實(shí)際地震資料獲得，主要考慮地震資料有效頻帶范圍； ② 地震數(shù)據(jù)匹配過(guò)程中的規(guī)則化.由于YM工區(qū)野外近地表?xiàng)l件簡(jiǎn)單，資料信噪比也高，五維數(shù)據(jù)規(guī)則化技術(shù)已很成熟，數(shù)據(jù)規(guī)則化難度不大，在數(shù)據(jù)預(yù)處理階段，采用正交匹配追蹤的五維數(shù)據(jù)規(guī)則化技術(shù)完成規(guī)則化； ③ 計(jì)算效率/迭代次數(shù)問(wèn)題，波形反演由于包含波動(dòng)方程地震正演計(jì)算，計(jì)算量非常大.提高計(jì)算效率主要通過(guò)兩方面：一是提高初始速度模型的精度；二是采用GPU/CPU聯(lián)合運(yùn)算.迭代次數(shù)主要根據(jù)對(duì)反演結(jié)果的質(zhì)控和反演效果來(lái)決定，否則容易陷入局部極小值; ④ 局部尋優(yōu)策略選用共軛梯度法.

5.3淺層折射波波形反演

淺層速度模型的反演精度是影響地下地質(zhì)結(jié)構(gòu)成像效果的主要因素之一.本文利用折射波，由低頻到高頻反演4000 m以內(nèi)較低頻的淺層速度，由于該區(qū)地震信號(hào)中折射波的能量較強(qiáng)，較豐富的低頻折射波有利于反演淺層速度信息.該方法只拾取初至后一段時(shí)間的波形信息反演，降低了全波形反演的非線性.為了進(jìn)一步提高反演穩(wěn)定性和收斂效果，采用時(shí)間域多尺度全波形反演算法，反演主頻依次為3 Hz(0～6 Hz)，5 Hz(0～10 Hz)，7 Hz(0～14 Hz).全波形反演前后折射波信號(hào)的匹配情況如圖8所示，圖中箭頭左側(cè)為實(shí)際數(shù)據(jù)、右側(cè)為正演模擬記錄，可以看出初始速度模型的匹配精度已經(jīng)很高，這也驗(yàn)證了我們初始模型的有效性，此外，模型更新后正演模擬記錄與實(shí)際數(shù)據(jù)更加匹配，說(shuō)明了淺層速度模型經(jīng)FWI更新后更加精準(zhǔn)；每個(gè)頻段經(jīng)10次左右迭代后，誤差目標(biāo)函數(shù)趨于收斂(如圖9所示).

圖8　FWI前(a)、后(b)不同主頻折射波匹配情況(箭頭左側(cè)為實(shí)際數(shù)據(jù)，右側(cè)為正演模擬記錄)Fig.8　Comparison of real refraction waves and modeling refraction waves before (a) and after (b) FWI (data left of arrow represents real waves, data right of arrow represents modeling waves)

圖9　不同主頻折射波FWI誤差目標(biāo)函數(shù)(左為主頻3 Hz、中為主頻5 Hz、右為主頻7 Hz)Fig.9　Misfit of FWI with refraction waves of different domain frequency (left:3 Hz, middle:5 Hz, right:7 Hz)

5.4中深層反射波FWI

針對(duì)中深層豐富的反射信號(hào)，采用折射波FWI的最終結(jié)果作為反射波FWI的初始模型，同樣進(jìn)行時(shí)間域多尺度反演，利用主頻依次為3.5 Hz(0～7 Hz)、4 Hz(0～10 Hz)、5 Hz(0～12 Hz)，反演較高頻的速度信息.全波形反演前后反射波信號(hào)的匹配情況如圖10所示，圖中箭頭左側(cè)為實(shí)際數(shù)據(jù)、右側(cè)為正演記錄，同樣可以看到模型更新后正演模擬記錄與實(shí)際數(shù)據(jù)匹配更好，說(shuō)明經(jīng)FWI更新后，中深層速度模型更加精準(zhǔn)；除了主頻為4 Hz，其余兩個(gè)主頻經(jīng)10次左右迭代后，誤差目標(biāo)函數(shù)趨于收斂(如圖11所示).

圖10　FWI前(a)、后(b)不同頻段反射波匹配情況(箭頭左側(cè)為實(shí)際數(shù)據(jù)，右側(cè)為正演模擬記錄)Fig.10　Comparison of real reflection waves and modeling reflection waves before (a) and after (b) FWI (data left of arrow represents real waves, data right of arrow represents modeling waves)

圖11　不同頻段反射波FWI誤差目標(biāo)函數(shù)(左為主頻3.5 Hz、中為主頻4 Hz、右為主頻5 Hz)Fig.11　Misfit of FWI with reflection waves of different domain frequency (left:3.5 Hz, middle:4 Hz, right:5 Hz)

經(jīng)過(guò)反射波FWI以后得到的速度模型即為全波形反演的最終速度模型.圖12為過(guò)某測(cè)線層析成像和全波形反演得到的速度模型對(duì)比圖，采用克?；舴虔B前深度偏移進(jìn)行偏移成像處理，結(jié)果如圖13和圖14所示.可以看出：基于FWI精細(xì)速度模型的疊前深度偏移有效改善了火成巖以下目的層的成像，消除了由于速度不準(zhǔn)造成的火成巖邊緣及以下斷層和破碎帶(圖13橢圓所示)的構(gòu)造假象以及奧陶系以下假低幅度隆起構(gòu)造(圖13箭頭所示)，同時(shí)對(duì)于深層縫洞型碳酸鹽巖儲(chǔ)層的成像也更加清晰(如圖14所示).

圖12　層析成像速度模型(a)和波形反演速度模型(b)對(duì)比Fig.12　Comparison of tomography result (a) and FWI result (b)

圖14　不同速度模型的PSDM 7600 m深度切片對(duì)比 (a) 層析成像速度PSDM深度切片； (b) 波形反演速度的PSDM深度切片.Fig.14　Depth slice (7600 m) of PSDM results with tomography velocity (a) and FWI velocity (b)

6火成巖對(duì)成像影響的正演論證

塔里木盆地YM地區(qū)二疊系大規(guī)模發(fā)育火成巖，巖性空間變化大、厚度變化也大.奧陶系目的層埋藏深，上覆地層速度的精度嚴(yán)重影響偏移成像的效果，不同巖性的火成巖速度差異大、厚度空間變化大，火成巖空間分布范圍和速度精度描述不準(zhǔn)會(huì)嚴(yán)重影響低幅構(gòu)造的落實(shí)和縫洞型碳酸鹽巖儲(chǔ)層空間成像位置.

為了驗(yàn)證圖13、圖14成像結(jié)果的可靠性和進(jìn)一步從地球物理角度證明，選擇橫跨火成巖的一條測(cè)線，根據(jù)初步地質(zhì)解釋成果、層析成像結(jié)果及測(cè)井資料，構(gòu)建簡(jiǎn)化2.5維的速度模型，如圖15a所示.采用25Hz雷克子波進(jìn)行有限差分正演模擬，然后對(duì)正演結(jié)果進(jìn)行克?；舴虔B前深度偏移，偏移結(jié)果如圖15b所示，可以看到箭頭所示處沒(méi)有斷層.

為了展示火成巖的空間邊界和速度精度對(duì)偏移成像結(jié)果的影響，分別對(duì)火成巖速度區(qū)域進(jìn)行1000 m局部平滑(圖16a)和3000 m局部平滑(圖16c)，采用平滑的速度模型分別對(duì)正演模擬的數(shù)據(jù)進(jìn)行疊前深度偏移處理，對(duì)應(yīng)的成像結(jié)果如圖16b和16d所示.可以看出，當(dāng)火成巖空間邊界和速度刻畫不準(zhǔn)時(shí)，容易導(dǎo)致火成巖邊緣以下對(duì)應(yīng)位置出現(xiàn)“假斷層”(如箭頭所示)，因此精細(xì)刻畫火成巖的空間邊界和速度對(duì)火成巖周緣及下伏地層的成像至關(guān)重要，火成巖刻畫不準(zhǔn)也會(huì)導(dǎo)致深部的縫洞型儲(chǔ)層成像誤差，同時(shí)證明了本文中全波形反演得到的速度模型的合理性和高精度，實(shí)際鉆井也已證實(shí)該結(jié)果.

7結(jié)論與認(rèn)識(shí)

全波形反演是一種理論完美，高精度的速度建模技術(shù)，初始速度模型的精度和足夠低頻的觀測(cè)地震數(shù)據(jù)是影響波形反演效果的兩個(gè)關(guān)鍵因素.本文針對(duì)塔里木盆地YM工區(qū)縫洞型儲(chǔ)層成像和地震資料特點(diǎn)，從影響波形反演關(guān)鍵因素角度出發(fā)，探索實(shí)現(xiàn)了面向縫洞型儲(chǔ)層的波形反演技術(shù)流程，并得到如下結(jié)論：

(1) 初始速度模型的精度是全波形反演成功與否的最大挑戰(zhàn)，淺層采用回折波層析成像、中深層采用反射波層析成像技術(shù)聯(lián)合建立各向異性速度模型精度高；

(2) 陸上資料低頻信息不足是全波形反演的瓶頸之一，與常規(guī)FKK方法相比，基于面波模擬和曲波變換聯(lián)合去噪技術(shù)能夠較好地保護(hù)低頻信息和波的動(dòng)力學(xué)特征，滿足全波形反演對(duì)數(shù)據(jù)的要求；

(3) 實(shí)現(xiàn)陸上資料波形反演的關(guān)鍵是建立合理流程，在建立較高精度初始速度模型的基礎(chǔ)上，采用折射波波形反演速度低頻信息，再用反射波反演速度高頻信息的波形反演速度迭代流程是一種行之有效的方法；

(4) 雖然目前大部分常規(guī)陸上地震采集的數(shù)據(jù)不能完全滿足經(jīng)典全波形反演的要求，但隨著采集新技術(shù)的應(yīng)用和處理技術(shù)的發(fā)展，尤其是波形反演在塔里木盆地YM地區(qū)的成功應(yīng)用，有望為陸上資料波形反演的工業(yè)化發(fā)展起到一定的推動(dòng)和借鑒作用.

References

Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion. Geophysics, 60(5): 1457-1473.

Cao S H, Chen J B. 2014. Studies on complex frequencies in frequency domain full waveform inversion. Chinese Journal of Geophysics (in Chinese), 57(7): 2302-2313, doi: 10.6038/cjg20140724.

Danecek P, Seriani G. 2008. An efficient parallel Chebyshev pseudo-spectral method for large scale 3D seismic forward modelling.∥70th EAGE Conference & Exhibition. SPE, EAGE.Dong L Q, Li Z C, Wang D Y, et al. 2011. Ground-roll suppression based on the second generation Curvelet transform. Oil Geophysical Prospecting (in Chinese), 46(6): 897-904.

圖13　不同速度模型的PSDM成像對(duì)比 (a) 層析成像速度PSDM結(jié)果； (b) 波形反演速度的PSDM結(jié)果.Fig.13　PSDM results with tomography velocity (a) and FWI velocity (b)

圖15　簡(jiǎn)化的速度模型(a)及對(duì)應(yīng)的疊前深度偏移結(jié)果(b)

圖16　平滑速度模型及其對(duì)應(yīng)的PSDM剖面

Fu L Y, Xiao Y J, Sun J W, et al. 2013. Seismic imaging studies of complex high and steep structures in Kuqa depression. Chinese Journal of Geophysics (in Chinese), 56(6): 1985-2001, doi: 10.6038/cjg20130620. Gauthier O, Virieux J, Tarantola A. 1986. Two-dimensional nonlinear inversion of seismic waveform: numerical results. Geophysics, 51(7): 1387-1403.

Han M. 2014. Methods and application of full waveform inversion for abyssal seismic data[Ph. D. thesis] (in Chinese). Changchun: Jilin University.Hollander Y, Kosloff D, Koren Z, et al. 2012. Seismic data interpolation by orthogonal matching pursuit.∥74th EAGE Conference & Exhibition. SPE, EAGE.

Hu G H, Jia C M, Xia H R, et al. 2013. Implementation and validation of 3D acoustic full waveform inversion. Geophysical Prospecting for Petroleum (in Chinese), 52(4): 417-425.

Levander A R. 1988. Fourth-order finite-difference P-SV seismogram. Geophysics, 53(11): 1425-1436.

Liu Y Z, Xie C, Yang J Z. 2014. Gaussian beam first-arrival waveform inversion based on Born wavepath. Chinese Journal of Geophysics (in Chinese), 57(9): 2900-2909, doi: 10.6038/cjg20140915.

Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics, 49(5): 533-549.

Mora P. 1987. Nonlinear two-dimensional elastic inversion of multi-offset seismic data. Geophysics, 52(9): 1211-1228.

Operto S, Virieux J, Amestoy P, et al. 2007. 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: a feasibility study. Geophysics, 72(5): SM195-SM211.

Plessix R é, Baeten G, de Maag J, et al. 2012. Full waveform inversion and distance separated simultaneous sweeping: a study with a land seismic data set. Geophysical Prospecting, 60(4): 733-747.

Pratt R G, Goulty N R. 1991. Combining wave-equation imaging with traveltime tomography to form high-resolution images from crosshole data. Geophysics, 56(2): 208-224.

Prieux V, Operto S, Brossier R, et al. 2010. Application of 2D acoustic frequency-domain full-waveform inversion to OBC wide-aperture data from the Valhall field.∥SEG Technical Program Expanded Abstracts, Denver, Colorado: SEG, 920-924.

Ratcliffe A, Win C, Vinje V, et al. 2011. Full waveform inversion: a North Sea OBS case study.∥Proceedings of the 81st Annual International Meeting. SEG, 2384-2388.

Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies. Geophysics, 69(1): 231-248. Sirgue L, Barkved O I, Dellinger J, et al. 2010. Full waveform inversion: the next leap forward in imaging at Valhall. First Break, 28(4): 65-70.

Song Z M, Williamson P R, Pratt R G. 1995. Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part II-Inversion method, synthetic experiments and real-data results. Geophysics, 60(3): 796-809.

Strobbia C L, Laake A, Vermeer P L, et al. 2009. Surface waves-use them then lose them.∥71th EAGE Conference & Exhibition. Amsterdam, The Newtherlands.

Takougang E M T, Calvert A J. 2011. Application of waveform tomography to marine seismic reflection data from the Queen Charlotte Basin of western Canada. Geophysics, 76(2): B55-B70. Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266.

Tnaer M T, Koehler F. 1981. Surface consistent corrections. Geophysics, 46(1): 12-22. Vigh D, Kapoor J, Li H Y. 2011. Full-waveform inversion application in different geological settings.∥SEG Technical Program Expanded Abstracts. SEG, 2374-2378.

Virieux J. 1984. SH-wave propagation in heterogeneous media, velocity-stress finite-difference method. Geophysics, 49(11): 1933-1942.

Wang W, Han B, Tang J P. 2013. Regularization method with sparsity constraints for seismic waveform inversion. Chinese Journal of Geophysics (in Chinese), 56(1): 289-297, doi: 10.6038/cjg20130130.

Wei Z F, Gao H W, Zhang J F. 2014. Time-domain full waveform inversion based on an irregular-grid acoustic modeling method. Chinese Journal of Geophysics (in Chinese), 57(2): 586-594, doi: 10.6038/cjg20140222.

Wu R S. 2003. Wave propagation, scattering and imaging using dual-domain one-way and one-return propagators. Pure and Applied Geophysics, 160(3-4): 509-539.

Yang Q Y, Hu G H, Wang L X. 2014. Research status and development trend of full waveform inversion. Geophysical Prospecting for Petroleum (in Chinese), 53(1): 77-83.

附中文參考文獻(xiàn)

曹書紅, 陳景波. 2014. 頻率域全波形反演中關(guān)于復(fù)頻率的研究. 地球物理學(xué)報(bào), 57(7): 2302-2313, doi: 10.6038/cjg20140724.

董烈乾, 李振春, 王德?tīng)I(yíng)等. 2011. 第二代Curvelet變換壓制面波方法. 石油地球物理勘探, 46(6): 897-904.

符力耘, 肖又軍, 孫偉家等. 2013. 庫(kù)車坳陷復(fù)雜高陡構(gòu)造地震成像研究. 地球物理學(xué)報(bào), 56(6): 1985-2001, doi: 10.6038/cjg20130620.

韓淼. 2014. 深水區(qū)地震全波形反演策略與應(yīng)用[博士論文]. 長(zhǎng)春: 吉林大學(xué).

胡光輝, 賈春梅, 夏洪瑞等. 2013. 三維聲波全波形反演的實(shí)現(xiàn)與驗(yàn)證. 石油物探, 52(4): 417-425. 劉玉柱, 謝春, 楊積忠. 2014. 基于Born波路徑的高斯束初至波波形反演. 地球物理學(xué)報(bào), 57(9): 2900-2909, doi: 10.6038/cjg20140915. 王薇, 韓波, 唐錦萍. 2013. 地震波形反演的稀疏約束正則化方法. 地球物理學(xué)報(bào), 56(1): 289-297, doi: 10.6038/cjg20130130.

魏哲楓, 高紅偉, 張劍鋒. 2014. 基于非規(guī)則網(wǎng)格聲波正演的時(shí)間域全波形反演. 地球物理學(xué)報(bào), 57(2): 586-594, doi: 10.6038/cjg20140222.

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

(本文編輯胡素芳)

基金項(xiàng)目國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(2013CB228604)，國(guó)家油氣重大專項(xiàng)(2011ZX05030-004-002, 2011ZX05019-003, 2011ZX05009-003-004)聯(lián)合資助.

作者簡(jiǎn)介崔永福， 1978年生，男，內(nèi)蒙古寧城縣人，工程師，中國(guó)石油大學(xué)(華東)在讀博士生，現(xiàn)主要從事石油物探方面的研究工作. E-mail：cuiyongfu-tlm@petrochina.com.cn

doi:10.6038/cjg20160734 中圖分類號(hào)P631

收稿日期2015-06-23，2016-01-12收修定稿

Application of full waveform inversion velocity model-building technology for the fractured-vuggy reservoir

AbstractVelocity is the key to pre-stack depth migration. Carbonate reservoir in Halahatang region of Tarim basin is buried deeply (deeper than 5500m) and has strong heterogeneity. Multi-period igneous rock with uneven thickness and different velocity in Permian had developed overlying target formation. Velocity modeling based on tomography is difficult to obtain accurate velocity field. In this paper, we carry out full waveform inversion (FWI) method to improve imaging accuracy of the fractured-vuggy reservoir.

KeywordsFull-waveform inversion; Land data; Velocity model-building; Fractured-vuggy reservoir; Igneous rock

崔永福， 彭更新， 吳國(guó)忱等. 2016. 全波形反演在縫洞型儲(chǔ)層速度建模中的應(yīng)用. 地球物理學(xué)報(bào)，59(7):2713-2725,doi:10.6038/cjg20160734.

Cui Y F, Peng G X, Wu G C, et al. 2016. Application of full waveform inversion velocity model-building technology for the fractured-vuggy reservoir. Chinese J. Geophys. (in Chinese),59(7):2713-2725,doi:10.6038/cjg20160734.