陸上地震資料全波形反演策略研究

王 杰,胡光輝,劉定進,邵文潮,王鵬燕

(中國石油化工股份有限公司石油物探技術研究院,江蘇南京211103)

陸上地震資料全波形反演策略研究

王 杰,胡光輝,劉定進,邵文潮,王鵬燕

(中國石油化工股份有限公司石油物探技術研究院,江蘇南京211103)

全波形反演已成功應用于海上地震資料，而陸上地震資料全波形反演應用還面臨諸多難點。研究給出了一種針對陸上地震資料的全波形反演策略：首先采用相位擬合互相關目標泛函增強全波形反演的適用性；其次構建動態波數域濾波梯度預處理方法壓制梯度噪聲；最后采用構造傾角信息約束的自適應高斯平滑算子改善反演結果質量。二維陸上地震資料的反演結果表明，該策略可有效實現陸上地震資料高波數全波形反演速度建模。

全波形反演;互相關目標泛函;分區波數域濾波;自適應各向異性高斯平滑

全波形反演(FWI)具有求取高分辨率、高精度地下物性參數的能力。TARANTOLA[1]于20世紀80年代首次提出了基于廣義最小二乘的時間域全波形反演技術。近10年來,隨著計算機技術的發展,全波形反演成為國內外地震勘探領域研究的熱點之一。目前全波形反演理論研究日臻成熟,但相較于線性射線層析,全波形反演由于匹配走時、振幅和相位等信息,非線性程度高,極易陷入局部極值產生周期跳問題,此外,所需計算資源巨大,因此,全波形反演的工業化應用面臨諸多挑戰。

傳統全波形反演是一高度非線性反問題,數據缺失低頻信息或初始模型不準時,模型在迭代更新過程中容易陷入局部極小值[2]。為降低反演的非線性性,BUNKS等[3]、SIRGUE等[4]分別從時間空間域和頻率空間域給出從低頻到高頻的多尺度反演策略。BROSSIER等[5]采用凸性更好、性態更優的互相關目標泛函完成全波形反演,相較常規L2范數而言,其局部極小值更少。WU等[6]從原始地震數據中提取包含豐富低頻信息的信號包絡反演宏觀背景速度以解決低頻缺失造成的局部極小值問題。LI等[7]采用相位追蹤和頻率外推方法從地震波高頻數據中人工模擬出低頻信息并完成了低波數背景速度重構,削弱了傳統全波形反演的非凸性。

巨大的計算成本是全波形反演在實際應用中的重要制約因素。為提高全波形反演計算效率,KREBS等[8]和BEN-HADJ-ALI等[9]分別在時間域和頻率域采用超級震源技術顯著減少每次迭代正演模擬的計算時間,并采用隨機相位編碼有效壓制反演結果中的串擾噪聲,在保證反演質量的前提下大大提高了全波形反演計算效率。

目前國際上海洋地震資料全波形反演成功應用實例越來越多。在SIRGUE等[10]采用三維寬方位海底電纜數據對北海油田成功構建全波形反演高精度速度模型之后,大偏移距、寬方位、寬頻海上拖纜數據、海底電纜數據、常規拖纜采集數據及寬頻數據、窄方位海底電纜數據等淺水環境和深水環境全波形反演成功應用案例越來越多。

相較海上資料,陸上地震資料全波形反演的成功應用還存在不小的挑戰[11]。在陸上地震數據采集過程中,受震源強度、震源激發響應、表層吸收衰減、各種干擾噪聲、復雜近地表傳播機制及檢波器耦合效應等影響,不同震源、不同檢波點地震信號振幅能量嚴重不均衡,局部地區存在震源稀疏,數據不完備現象,且受采集儀器限制,地震信號往往缺失低頻及大偏移距數據信息。PLESSIX等[12-13]采用專門為全波形反演采集的低頻、大偏移距、高密度陸地可控震源數據完成了內蒙古地區二維高精度速度模型建模,并取得較好應用效果。HE等[14]指出,當陸上地震資料數據不完備時,反演結果中存在較強噪聲,迭代收斂性和穩健性都有一定程度影響,為此他們給出了一種有效的噪聲算子濾波方法,然而當地質條件較為復雜、構造起伏劇烈時,該方法仍存在一定的局限性。我國以陸上探區為主,研發適合于陸上地震資料的全波形反演技術迫在眉睫。

本文給出適用于不完備陸上地震資料的穩健的全波形反演策略:采用相位匹配的非零延遲歸一化互相關目標泛函提高反演問題的凸性和穩健性,降低不同震源和檢波器振幅能量不均衡對反演造成的影響;基于已知工區構造先驗信息將反演工區分塊并對不同區塊在波數域完成反演噪聲壓制;對整體模型更新量采用自適應高斯平滑技術減弱不同區塊模型更新量的耦合效應,從而使該反演策略適用于任意地質構造條件下陸上地震資料的高精度速度建模。最后采用東北某探區實際二維地震資料進行了應用測試。

1 方法原理

1.1 互相關目標泛函反演框架

L2范數全波形反演在迭代過程中需要與地震數據的振幅、波形、相位信息匹配,其對數據振幅信息非常敏感。由于陸上地震資料激發和接收條件等因素的限制,原始采集的地震記錄存在炮間、道間能量不均衡,在局部區域非常稀疏。當采用此類數據進行凸性相對較差的L2范數全波形反演時,不均衡的殘差能量會引起模型更新量在不同地區存在較大差異,從而導致反演不收斂或不穩定,使反演結果陷入局部極值。為提高反演的穩健性,在進行陸上地震資料反演時需要給旅行時更多的權重,以降低振幅的影響?；谙辔粩M合的互相關目標泛函更強調匹配數據的相位信息,凸性更優,且相對于L2范數具有更廣的波谷,因而反演結果更加穩健[5,15]。本文采用BROSSIER等[5]給出的非零延遲歸一化互相關目標泛函理論框架實現陸上地震資料的全波形反演。非零延遲歸一化互相關目標泛函為:

(1)

式中:m為要反演的模型參數;τ為延遲時間;xr為檢波器坐標;xs為震源坐標;W(τ)為懲罰函數;dcal(t,xr;xs)為人工合成地震記錄;dobs(t+τ,xr;xs)為時延觀測地震記錄。

即使觀測地震記錄和模擬地震記錄相位差大于半個地震波波長,該目標泛函亦可穩定收斂,而且歸一化的引入也減小了不同震源和檢波器振幅能量不均衡引起的噪聲干擾問題。梯度表達式可用伴隨狀態法求解,但與L2范數不同,非零延遲歸一化互相關目標泛函需修改反傳震源項[5]。采用目標泛函對模型參數的導數可求解梯度表達式:

(2)

1.2 分區波數域濾波算子

NANGOO[16]采用稀疏震源子集對常規缺乏低頻和大偏移距的窄方位海上拖纜數據進行了全波形反演研究,采用水平光滑預處理技術壓制梯度大傾角噪聲,保留反演速度中的地質構造信息。當工區地質背景較為復雜時,水平光滑預處理技術不再適用。相對海上拖纜數據,陸上資料的稀疏性更強,而且信噪比低,數據質量相對更差。當工區震源較為稀疏,且覆蓋次數較低時,全波形反演會在梯度中引入強烈的反演噪聲[14,16]。為得到較為準確的、符合地質構造背景的速度模型,需要壓制迭代過程中產生的反演噪聲。HE等[14]提出了在波數域壓制噪聲的方法,但該方法的不足之處在于若選擇的構造傾角較小,則高陡構造也被壓制;若選擇的構造傾角過大,雖然在濾波后的反演結果中保留了高陡構造特征,但會帶來保留較強噪聲干擾的結果。本文為突破上述單角度波數域濾波的限制,基于已有地質構造信息,對地下介質分區,采用分而治之思想,對不同區塊選取相應角度的波數域濾波算子壓制噪聲,在有效壓制噪聲的前提下,盡可能保留地質構造特征信息。繼而采用整體自適應高斯平滑濾波技術減小不同區塊梯度波數域濾波后的不一致性。

分區波數域濾波需要在每次迭代過程中按照成像剖面沿地質構造傾角分區。局部構造傾角場可采用局部平面波解構濾波法[17]或圖像梯度結構張量算法[18]求解,本文采用局部平面波解構法求解構造傾角場。平面波解構濾波器利用局部平面波的疊加表征地震數據,可以較好地估計平滑連續同相軸的局部傾角信息。根據局部平面波的物理模型,用局部有限差分方程定義平面波解構濾波器為:

(3)

式中:P(x,t)為局部平面波波場;σ(x,t)為平面波局部斜率場。(3)式可化簡為基于最小二乘的傾角估計

問題:

(4)

式中:σ0為斜率初始值;Δσ為斜率增量;C為局部傾角的函數;d為原始數據。利用(4)式可構建局部傾角場。根據構造的局部傾角場將整個地質工區分為不同區塊,每個區塊對應不同的主要構造傾角。利用不同區塊的主要構造傾角和各自相應的波數域濾波算子對不同區塊梯度壓制噪聲。各區塊采用的噪聲壓制算子公式為:

(5)

1.3 自適應各向異性高斯平滑算子

分區濾波梯度場在整合為整體梯度時不同區塊交界處存在邊界效應,這將在反演結果中引入明顯的區塊邊界。為削弱這種邊界效應,本文對整體梯度場采用自適應各向異性高斯平滑濾波算子[19]沿傾角方向濾波。各向異性高斯濾波算子在x,y平面上的投影為橢圓(如圖1a所示),該濾波算子表達式為:

(6)

針對不同的構造傾角θ,采用旋轉的方式得到如圖1b所示的投影,由坐標變換可得旋轉后的坐標u-v和x-y坐標的關系為:

(7)

旋轉后的高斯濾波算子為:

(8)

圖1 各向異性高斯平滑算子無旋轉(a)和旋轉θ(b)在x-y平面上的投影

采用構建的自適應各向異性高斯濾波算子可對整合后的梯度場沿著構造傾角平滑:

(9)

式中:sn為整合梯度場;dn表示高斯濾波梯度場;G為自適應各向異性高斯濾波算子。

在構建出構造較為連續、噪聲污染較少的梯度場后,采用預條件擬牛頓(L-BFGS)迭代最優化算法對速度模型迭代更新,迭代步長采用非精確線性搜索技術求取。

2 實際資料應用

為驗證本文方法的有效性,選取東北某探區一條二維線進行全波形反演。沿探區東西方向抽取44炮單炮地震記錄,每炮180道左右。震源在x方向的分布范圍為13000m,炮點距350m左右,相對稀疏。圖2為該二維線覆蓋次數圖。其中,最大覆蓋次數為16次。圖3a為去除面波、線性噪聲、異常振幅,并進行地表一致性振幅補償、球面擴散補償后的原始單炮地震記錄。圖3b為6Hz以上頻率段觀測地震記錄。該記錄中面波得到了有效壓制,信噪比得以提高。圖3c 為圖3b的低頻段信息,頻率范圍為6～20Hz。本文采用6～20Hz頻率段數據進行反演測試。圖4為二維深度域克?；舴虔B前偏移剖面。對偏移剖面采用局部平面波解構濾波可獲取地質構造傾角信息(圖5)。從構造傾角場中可以看出,淺層的主要反射層構造傾角較小,大部分集中在10°之內,中、深層構造傾角較大,為-25°～35°。

地震子波是傳統基于L2范數實際資料全波形反演成敗的一個關鍵因素[11]。由于實際資料的地震子波是時變和空變的,而且在相同條件下,不同子波正演模擬出的人工合成地震記錄在振幅、相位和波形上均存在較大差別,因而如何做好地震子波估計是全波形反演的重要環節。本文采用基于相位擬合的非零延遲歸一化互相關目標反演理論框架完成速度場的重構,在反演迭代過程中削弱了原始地震記錄和觀測地震記錄振幅和波形的影響,更加強調相位的作用,在一定程度上減小了對地震子波估計的要求,因此本文采用合適主頻的雷克子波完成正演模擬及速度重構,非零延遲的引入在一定程度上允許地震子波具有一定的相位延遲。本文采用的初始速度模型為反射波旅行時層析速度反演模型(圖6),該模型具有較為充分的低波數信息及宏觀地質構造背景信息。

圖2 研究區某二維線覆蓋次數分布

圖3 研究區某二維線單炮地震記錄a 去噪和振幅補償后的觀測地震記錄; b 消除面波后的觀測地震記錄; c 6～20Hz觀測地震記錄

采用原始稀疏觀測地震記錄進行互相關全波形反演時,由于覆蓋次數較少,反演的非線性程度較高,反演結果存在較強的噪聲。圖7顯示了6～20Hz數據迭代20次的全波形反演結果。由圖7可見,反演結果信噪比較低,有效地質構造信息被反演噪聲湮沒。圖8給出了在每次迭代過程中采用單一角度(10°,15°,25°,35°)波數域濾波迭代20次得到的反演結果。相較于原始資料直接重構的速度模型,采用單一角度波數域濾波后構建的反演速度有效壓制了串擾噪聲,增強了有效地質構造的連續性。然而單一角度波數域濾波本身存在缺陷,當采用的濾波算子角度過小時,雖然反演結果信噪比顯著提高,但濾波算子消除了大角度的地質構造信息(圖8a,圖8b)。由于該工區的構造傾角范圍為-25°～35°,當采用大角度波數域濾波算子時,雖然大傾角的地質構造信息得以保留,但是反演結果中也相應引入了較為嚴重的大角度線性噪聲(圖8c,圖8d)。

圖4 深度域克?；舴虔B前偏移剖面

圖5 疊前偏移剖面構造傾角場

圖6 反射波旅行時層析速度模型

圖7 互相關目標泛函全波形反演結果

該工區構造傾角在淺層集中在10°以下,深層主要集中在25°左右,因此本文采用淺、深層分區波數域濾波算子完成全波形反演建模。圖9給出了分區波數域濾波建模結果。由圖9可見,淺層和深層都保留了有效地質構造信息,與單一大角度波數域濾波反演結果(圖8c,圖8d)相比,分辨率有一定改善,但是在紅色箭頭處,出現低速異常體,反演結果陷入局部極值,且白色箭頭的位置速度出現不耦合現象。為此,每次迭代在分區波數域濾波的基礎上引入各向異性自適應高斯平滑算子。圖10顯示了最終的重構模型,該重構結果不僅信噪比提高而且構造更加連續,局部極值問題得到削弱。

圖8 采用單一角度波數域濾波后互相關目標泛函全波形反演結果a 濾波角度10°; b 濾波角度15°; c 濾波角度25°; d 濾波角度35°

圖11給出了初始模型、反演速度模型正演模擬歸一化地震記錄和歸一化觀測地震記錄的對比結果。

圖9 分區波數域濾波互相關全波形反演結果(淺層10°,深層15°)

圖10 分區波數域濾波后采用各向異性自適應高斯平滑處理得到的互相關全波形反演結果

圖11 觀測地震記錄與模擬地震記錄結果對比a 初始模型正演模擬地震記錄(左)與觀測地震記錄(右); b 反演模型正演模擬地震記錄(左)與觀測地震記錄(右)

與初始模型正演模擬地震記錄相比,反演模型的模擬地震記錄主要反射層的相位與觀測地震記錄匹配得更好,從相位擬合的角度驗證了反演的準確性與穩健性。

3 結論

本文給出陸上地震資料不完備條件下采用全波形反演進行高精度速度建模的方案。采用基于相位擬合的非零延遲歸一化互相關目標泛函取代傳統的基于數據波形、振幅、相位匹配的L2范數目標泛函,提高了反演的穩健性,且降低了對地震子波精度和原始地震資料振幅能量均衡的要求。采用分區波數域濾波算子和各向異性自適應高斯平滑技術有效減少了原始地震資料稀疏時反演結果中較強的反演噪聲并保留了有效地質構造信息,反演結果更加合理、準確,而且反演過程較為穩健。二維陸上地震資料反演測試結果表明本文給出的陸上地震資料反演策略不僅有效消除了反演結果中強烈的噪聲,而且構建了具有地質意義的高分辨率速度模型,波形的匹配表明本文反演方法相對層析反演結果更優。

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

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

[3] BUNKS C,SALECK F.Multi-scale seismic waveform inversion[J].Geophysics,1995,60(5):1457-1473

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

[5] BROSSIER R,OPERTO S.Velocity model building from seismic reflection data by full-waveform inversion[J].Geophysical Prospecting,2015,63(2):354-367

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

[7] LI Y Y,DEMANET L.Full-waveform inversion with extrapolated low-frequency data[J].Geophysics,2016,81(6):R339-R348

[8] KREBS J R,ANDERSON J E,HINKLEY D,et al.Fast full-wavefield seismic inversion using encoded sources[J].Geophysics,2009,74(6):WCC177-WCC188

[9] BEN-HADJ-ALI H,OPERTO S,VIRIEUX J.An efficient frequency-domain full waveform inversion method using simultaneous encoded sources[J].Geophysics,2011,76(4):R109-R124

[10] SIRGUE L,BARKVEL O.3D waveform inversion in Valhall wide-azimuth OBC[J].Expanded Abstracts of 71stAnnual Internat EAGE Mtg,2009:U038

[11] 楊勤勇,胡光輝,王立歆.全波形反演研究現狀及發展趨勢[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

[12] PLESSIX R E,BAETEN G.Application of acoustic full waveform inversion to a low-frequency large-offset land dataset[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:930-934

[13] PLESSIX R E,BAETEN G.Full waveform inversion and distance separated simultaneous sweeping:a study with a land seismic data set[J].Geophysical Prospecting,2012,60(4):733-747

[14] HE B H,WU G C.Robust full-waveform inversion with incomplete refraction[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:150-153

[15] MASONI I,BROSSIER R.Alternative misfit functions for FWI applied to surface wave[J].Expanded Abstracts of 75thAnnual Internat EAGE Mtg,2013:Th P10 13

[16] NANGOO T P.Seismic full-waveform inversion of 3D field data——from the near surface to the reservoir [D].London:University of Imperial College,2013

[17] FOMEL S.Applications of plane-wave destruction filters[J].Geophysics,2002,67(10):1946-1960

[18] HALE D.Local dip filtering with directional Laplacians[R].CWP report,2007

[19] 王懷野,張科,李言俊.一種自適應各向異性高斯濾波方法[J].計算機工程與應用,2004,40(10):18-19 WANG H Y,ZHANG K,LI Y J.An adaptive anisotropic Gaussian filter for noise reduction[J].Computer Engineering and Applications,2004,40(10):18-19

(編輯：陳 杰)

Strategy study on full waveform inversion for the land seismic data

WANG Jie,HU Guanghui,LIU Dingjin,SHAO Wenchao,WANG Pengyan

(SinopecGeophysicalResearchInstitute,Nanjing211103,China)

Full waveform inversion of marine seismic data has been successfully applied,however,the application of this technology for onshore seismic data faces many difficulties.In this study,a robust strategy of full waveform inversion for onshore seismic data is presented.Firstly,the phase-matched cross-correlation cost function is adopted to improve the applicability of full waveform inversion.Secondly,dynamic wavenumber domain gradient preprocessing method is used to suppress gradient noise.Finally,an adaptive Gaussian smoothing operator is applied to improve the quality of inversion results.The inversion case of two-dimensional onshore seismic data shows that this strategy presented in this study can effectively realize high wavenumber onshore seismic data velocity modeling.

fullwaveform inversion,cross-correlation cost function,sub-regional wavenumber domain filtering,adaptive anisotropic Gaussian smoothing

2016-11-10;改回日期：2016-12-01。

王杰(1988—),男,碩士,工程師,主要從事全波形反演速度建模研究。

國家科技重大專項(2016ZX05014-001-002)資助。

P631

A

1000-1441(2017)01-0081-08

10.3969/j.issn.1000-1441.2017.01.010

This research is financially supported by National Science and Technology Major Project of China (Grant No.2016ZX05014-001-002).