基于頻域衰減的時域全波形反演

郭雪豹， 劉洪， 石穎

1 中國科學(xué)院地質(zhì)與地球物理研究所中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京　100029 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京　100029 3 東北石油大學(xué)， 大慶　163318

?

基于頻域衰減的時域全波形反演

郭雪豹1,2， 劉洪1,2， 石穎3

1 中國科學(xué)院地質(zhì)與地球物理研究所中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京100029 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京100029 3 東北石油大學(xué)， 大慶163318

時域全波形反演由于采用了全頻段信息，因此在迭代過程中不同波長的信息不能由低到高的逐步重建，極易陷入局部極小值.本文通過分頻段的方式，對地震數(shù)據(jù)做正反傅里葉變換，利用頻域指數(shù)衰減的方法逐級分離出地震數(shù)據(jù)中的高頻成分，在時域上實(shí)現(xiàn)由低頻向高頻的波形反演，從而降低了反演的非線性，使不同波長的信息得到穩(wěn)步恢復(fù).同時，在高頻成分衰減的過程中，后至波的能量也被削弱，由此也降低了深層反射在初始反演過程中的干擾.整個反演僅增加對數(shù)據(jù)做正反傅里葉變換過程，相較于混合域反演，無需提取全部波場的相應(yīng)頻率成分.在計算效率方面，利用GPU進(jìn)行加速，并采用CUDA自帶函數(shù)庫中cufft來提高計算效率.通過對Marmousi模型測試，驗(yàn)證了所述方法的有效性.

時域全波形反演；頻域衰減；低頻；GPU

1　引言

全波形反演是獲取地下參數(shù)的重要手段之一，其主要通過最小化觀測數(shù)據(jù)與預(yù)測數(shù)據(jù)的殘差來得到高精度的地下參數(shù).Tarantola(1984,1986)起初是在時域上實(shí)現(xiàn)了聲波的全波形反演，隨后被其應(yīng)用到彈性波中.全波形反演是一個高度非線性問題，限于計算壓力，無法在實(shí)際應(yīng)用中通過全局最優(yōu)化的方法得到最優(yōu)解，而普遍采用局部優(yōu)化的方法.因此，反演的效果很大程度上取決于初始模型的選擇.許多地球物理學(xué)家對此問題進(jìn)行研究，以降低反演對初始速度模型的依賴.Biondi和Almomin(2012)借助層析來獲得一個較好的初始速度.Almomin和Biondi (2012)采用波動方程偏移速度分析的方法.Bunks(1995)提出了多尺度的反演策略以降低反演的非線性.Pratt(1999a,1999b)、Brenders和Pratt(2007)將全波形反演應(yīng)用到頻域，并發(fā)展了由低頻向高頻的反演策略.Shin和Cha(2008,2009)提出了拉普拉斯域全波形反演，即便采用線性初始模型也可以得到一個好的反演結(jié)果.在這些方法中，多尺度反演大都需要低頻信息去恢復(fù)模型中的長波長分量.而針對常規(guī)地震數(shù)據(jù)中的低頻缺失，Wu等(2013,2014)和Luo等(2013)也提出包絡(luò)反演，并將其應(yīng)用到全波形反演中以彌補(bǔ)低頻的缺失.

最小化殘差本質(zhì)上是一個迭代的過程，在每次迭代過程中，都需要計算對應(yīng)模型參數(shù)的正傳波場和殘差的反傳波場，相關(guān)來計算梯度，因此波場模擬的效率決定了全波形反演的計算效率.對于二維情況來說，頻域可以采用直接法進(jìn)行求解，并且可以多炮并行處理，較時域方法具有明顯的優(yōu)勢.但對于三維情況，頻域求解變得十分困難，因此，通常在時域完成三維反演.而對于時域波場模擬而言，GPU的特點(diǎn)使得其十分適用于這種大規(guī)模的并行計算，Moghaddam(2013)將GPU應(yīng)用到全波形反演(Yang et al.,2015)中，Boonyasiriwat等人(2010)將其應(yīng)用于逆時偏移(Yang et al.,2014; Shin et al.,2014; 李博等,2010; 劉紅偉等,2011; 劉守偉等,2013)中，石穎等人(2010,2013)則將其用于多次波消除的加速計算中.

本文采用了一種基于頻域衰減的時域全波形反演方法.其主要通過對觀測數(shù)據(jù)進(jìn)行高頻衰減來達(dá)到在時域上實(shí)現(xiàn)由低頻向高頻的反演.在頻域衰減的同時，后至波的能量也在一定程度上被衰減，由此也降低了深層反射在初始反演過程中的干擾.與傳統(tǒng)時域全波形反演相比，實(shí)現(xiàn)了由低頻到高頻、由淺層到深層的反演，并同時具有時域正演可高度并行的優(yōu)勢，可直接應(yīng)用于三維.與混合域全波形反演相比，無需提取全部波場的相應(yīng)頻率成分，僅需對地震記錄做傅里葉變換即可.本文算法基于GPU進(jìn)行加速，并采用CUDA自帶函數(shù)庫的cufft，從而保證算法的穩(wěn)定高效.通過對Marmousi模型進(jìn)行測試，驗(yàn)證了本文方法的有效性.

2　基本原理

2.1時域全波形反演基本理論

常密度聲波方程可表示為如下形式：

(1)

因此，接收點(diǎn)處殘差為

(2)

其中，ucal和uobs分別代表計算得到的地震記錄和實(shí)際觀測地震記錄，xs和xr分別為震源坐標(biāo)和檢波點(diǎn)坐標(biāo).

(2)式在最小二乘意義下的目標(biāo)函數(shù)為

(3)

其中，m為待反演的參數(shù).

參數(shù)更新公式：

(4)

其中，在第k+1次迭代中步長αk和方向dk可以分別表示如下：

(5)

為補(bǔ)償照明，梯度歸一化公式(Gauthieretal.,1986;Baietal.,2014;Shinetal.,2001)為

(6)

其中，γ為穩(wěn)定性系數(shù).為了減少噪音干擾，本文在梯度計算后，采用高斯平滑方法壓制噪音.

2.2頻域衰減

頻域衰減公式為

Δu′(xs,xr,t)=

IFFT[FFT[ucal(xs,xr,t)-uobs(xs,xr,t)]G(f)]，

(7)

其中，F(xiàn)FT和IFFT分別為正快速傅里葉變換和反快速傅里葉變換，G(f)為衰減函數(shù)，f為頻率，Δu′(xs,xr,t)為衰減后的波場數(shù)據(jù)殘差.

本文中采取的衰減函數(shù)為

(8)

其中，Q為頻域衰減因子，其隨著頻率的增加而增加.本文采用的是e指數(shù)衰減，也可采用其他形式的衰減函數(shù).

2.3GPU_CUFFT

GPU的引入能夠在極大程度上提高并行算法的計算效率(Yang et al.,2014; Shin et al.,2014; 李博等,2010; 劉紅偉等,2011; 劉守偉等,2013; 石穎等, 2010, 2013; Yang et al.,2015; Boonyasiriwat et al., 2010).本文算法本身搭建在GPU平臺下，并通過存儲器優(yōu)化，進(jìn)一步提高計算過程中的數(shù)據(jù)讀取速度，使得計算效率得到明顯提升.

文中所述的頻域衰減，需要對多炮地震數(shù)據(jù)做傅里葉變換，這一過程的計算效率對整個算法的計算效率有較大影響.由于整個算法利用GPU進(jìn)行加速，因此最好將這一變換過程定義在GPU上進(jìn)行，這樣不僅能夠減少數(shù)據(jù)的讀寫，也可充分利用GPU并行計算的能力.在這里對GPU的離散傅里葉變換(以下簡稱CUDA_DFT)與利用GPU庫的cufft進(jìn)行對比，可以看出利用CUDA庫函數(shù)的CUFFT相較于CUDA_DFT，計算優(yōu)勢明顯.

計算速度對比見表1(以一炮767道，每道2800采樣點(diǎn)為例)，完成40炮計算時，CUFFT較CUDA_DFT計算效率提高約49倍.

表1　CUDA_DFT與CUFFT計算機(jī)時對比Table 1　Contrast of cost of CUDA_DFT and CUFFT

3　模型測試

本文采用Marmousi模型測試所述方法，Marmousi模型如圖1所示，縱、橫向網(wǎng)格點(diǎn)數(shù)分別為251和767，網(wǎng)格間距為12 m，雷克子波主頻為10 Hz，時間采樣間隔1.5 ms，采樣點(diǎn)數(shù)為2800.設(shè)計21炮震源激發(fā)，第一炮位于深度36 m，自左向右36 m位置處，向右每隔456 m放一炮.每炮767個檢波器接收，檢波器深度36 m，自左向右沿橫向依次排列，檢波器間距12 m.反演初始速度為真實(shí)速度平滑后的結(jié)果，如圖2所示.地震正演得到的第5炮和第10炮的炮集記錄如圖3所示.

圖1　Marmousi模型Fig.1　Marmousi model

圖2　初始速度模型Fig.2　Initial model

首先不考慮頻域衰減，即在無任何衰減情況下反演，迭代100次，200次，300次的結(jié)果如圖4所示.觀察易知，隨著迭代次數(shù)的增加，其反演結(jié)果逐漸趨于穩(wěn)定.對比速度模型可以看出，其反演結(jié)果并不理想，尤其在淺部存在較強(qiáng)的噪音干擾.由此可以看出，即便采用平滑后的真實(shí)速度模型，時域的全頻帶反演沒有逐步的反演出不同的波長分量，依舊會使結(jié)果陷入局部極小值.

圖3　炮集記錄(a) 第5炮地震記錄; (b) 第10炮地震記錄.Fig.3　Shot gathers(a) The fifth shot gather; (b) The tenth shot gather.

圖4　傳統(tǒng)時域反演(a) 迭代100次; (b) 迭代200次; (c) 迭代300次.Fig.4　Results of conventional FWI(a) 100 iterations; (b) 200 iterations; (c) 300 iterations.

圖5為不同迭代次數(shù)的單道對比圖，當(dāng)?shù)螖?shù)超過100次后，雖然反演依舊向真實(shí)值趨近，但隨著迭代次數(shù)的增加，這種改善越來越小，收斂速度變慢，并趨近于局部極小值，沒有收斂到全局最優(yōu)解.在700 m以上的淺部，當(dāng)?shù)螖?shù)達(dá)到300時，已經(jīng)基本與真實(shí)速度吻合或逐漸逼近.但隨著深度的增加，其收斂程度明顯降低，并與真實(shí)速度存在較大的偏差，說明時域的全頻段反演對深層影響更大.

而后采取頻率衰減來進(jìn)行反演，Q分別取10,20,40,80,0(0代表不做任何衰減)，Q為10時，迭代40次；Q為20時，迭代40次；Q為40時，迭代40次；Q為80時，迭代40次；不衰減時，迭代140次，同樣迭代300次，每次迭代均以上一次迭代后的結(jié)果為初始速度，結(jié)果分別如圖6所示，對比可以看出，由于逐級分離出地震記錄的高頻成分，將其數(shù)據(jù)限制在低頻端，降低了高頻成分對反演的影響，使得不同的波長分量得以逐級反演出來.從圖6a中可以看出，僅40次迭代，淺部形態(tài)就已經(jīng)基本反演出來.并且隨著衰減程度的降低，在淺層精度提高的同時，深部形態(tài)也逐漸趨于真實(shí).

在頻域衰減過程中由于采用CUDA自帶的庫函數(shù)，保證了對地震記錄做正反傅里葉變換的計算效率，每次迭代僅增加約2 s.不同迭代次數(shù)單道對比如圖7和圖8所示，圖9為不同迭代次數(shù)對比圖.頻域衰減幅度較大時，少量迭代即可快速收斂，與真實(shí)速度吻合較好，隨著迭代次數(shù)的增加，深層較淺層改善明顯.當(dāng)無衰減40次迭代時，收斂速度已趨于穩(wěn)定.繼續(xù)迭代140次后，反演結(jié)果變化較小.并且對比衰減下的深淺層可以看出，衰減參數(shù)直接作用于記錄也消弱了后至波對淺層反演干擾，實(shí)現(xiàn)了由淺層到深層的反演，提高了深層反演的精度.

圖5　4800 m處單道對比圖(a) 迭代100次; (b) 迭代200次; (c) 迭代300次; (d) 不同迭代次數(shù)對比.Fig.5　Contrast of single traces between different iterations (Distance=4800 meters)(a) 100 iterations; (b) 200 iterations; (c) 300 iterations; (d) Contrast of different iterations.

圖6　頻域衰減反演結(jié)果(a) Q=10迭代40次; (b) Q=20迭代40次; (c) Q=40迭代40次; (d) Q=80迭代40次; (e) Q=0迭代40次; (f) Q=0迭代140次.Fig.6　Results of FWI using frequency damping(a) 40 iterations (Q=10); (b) 40 iterations (Q=20); (c) 40 iterations (Q=40);(d) 40 iterations (Q=80); (e) 40 iterations (Q=0); (f) 140 iterations (Q=0).

為了進(jìn)一步對比迭代效果，將無衰減300次迭代結(jié)果與逐級衰減得到的結(jié)果相對比，單道對比如圖10和圖11所示.可以看出，Q為10時，曲線大致趨勢已與真實(shí)速度開始吻合，在淺部一些突變構(gòu)造上幅度雖然有所缺失，但位置卻較無衰減300次迭代的結(jié)果更加準(zhǔn)確.并且在深部吻合的更好.隨著衰減程度的減弱，當(dāng)Q為40時，已經(jīng)與真實(shí)速度的位置貼近.當(dāng)?shù)螖?shù)繼續(xù)增加時，深層開始慢慢改善，與真實(shí)速度的差異慢慢減小.而無衰減300次迭代的結(jié)果在深度依舊保留了一定的誤差.

頻域衰減后，不同衰減因子的單道信號匹配圖，如圖12和13所示，分別為偏移距2400 m處的單道衰減0～1 s、3～4 s部分.對比圖12a、12b、12c、12d可以看出，當(dāng)Q逐漸增大時，地震記錄中的高頻成分被逐級的剝離，僅保留需要頻段成分的地震記錄，類似于頻域的由低頻向高頻反演，大大減少了高頻成分對反演的影響，降低了波形錯位的可能.對比圖13a、13b、13c、13d，在頻域衰減過程中，后至波被一定程度地消弱，也降低了后至波在反演初始階段對淺層反演的影響，等同于由淺層至深層反演.

4　結(jié)論

本文將頻域衰減應(yīng)用于時域全波形反演，以此來實(shí)現(xiàn)由低頻到高頻的逐步反演，并在頻域衰減的同時，后至波的能量也在一定程度上被削弱，即頻域衰減也同時壓制了深層反射，實(shí)現(xiàn)由淺層到深層的反演.當(dāng)衰減程度較大時，僅需少量迭代即可得到一個較好的結(jié)果，很大程度上減少了波形錯位現(xiàn)象，提高了收斂速度.在計算效率方面，采用GPU進(jìn)行加速，利用CUFFT實(shí)現(xiàn)頻域衰減，因此計算效率方面得到有效保證.本文方法屬于時域反演，但其避免了傳統(tǒng)時域方法全頻帶反演帶來的局部極小值問題.與混合域、Laplace-Fourier域全波形反演相比，不存在其他域的波場計算，完全建立在時域的基礎(chǔ)上，并且由于采用了一個頻段的信息，反演更為穩(wěn)健，在選取參數(shù)方面也更為靈活.同時，本文方法也保留了時域方法在三維方面的優(yōu)勢.

圖7　4800 m處單道對比(a) Q=10迭代40次; (b) Q=20迭代40次; (c) Q=40迭代40次.Fig.7　Contrast of single traces between real and Q (Distance=4800 meters).(a) 40 iterations (Q=10); (b) 40 iterations (Q=20); (c) 40 iterations (Q=40).

圖8　4800 m處單道對比(a) Q=80迭代40次; (b) Q=0迭代40次; (c) Q=0迭代140次.Fig.8　Contrast of single traces between real and different Q(Distance=4800 meters).(a) 40 iterations (Q=80); (b) 40 iterations (Q=0); (c) 140 iterations (Q=0).

圖9　4800 m處單道不同迭代次數(shù)對比(a) Q=10, Q=20, Q=40單道對比; (b) Q=20, Q=40, Q=80單道對比; (c) Q=80, Q=0, Q=0單道對比.Fig.9　Contrast of single traces between different Q values (Distance=4800 meters)(a) Contrast of single traces (Q=10, Q=20, Q=40); (b) Contrast of single traces (Q=20, Q=40, Q=80); (c) Contrast of single traces (Q=80, Q=0, Q=0).

圖10　4800 m處兩種方法單道對比(a) 原始方法迭代300次與Q=10迭代40次對比; (b) 原始方法迭代300次與Q=20迭代40次對比; (c) 原始方法迭代300次與Q=40迭代40次對比.Fig.10　Contrast of single traces between two methods (Distance=4800 meters).(a) Contrast of singe traces (300 iterations of conventional method and 40 iterations of Q=10); (b) Contrast of singe traces (300 iterations of conventional method and 40 iterations of Q=20); (c) Contrast of singe traces (300 iterations of conventional method and 40 iterations of Q=40).

圖11　4800 m處兩種方法單道對比(a) 原始方法迭代300次與Q=80迭代40次對比; (b) 原始方法迭代300次與Q=0迭代40次對比; (c) 原始方法迭代300次與Q=0迭代140次對比.Fig.11　Contrast of single traces between two methods (Distance=4800 meters).(a) Contrast of singe traces (300 iterations of conventional method and 40 iterations of Q=80); (b) Contrast of singe traces (300 iterations of conventional method and 40 iterations of Q=0); (c) Contrast of singe traces (300 iterations of conventional method and 140 iterations of Q=0).

圖12　2400 m處0～1 s單道對比(a) Q=0與Q=10單道對比; (b) Q=0與Q=20單道對比; (c) Q=0與Q=40單道對比; (d) Q=0與Q=80單道對比.Fig.12　Contrast of single traces (Distance=2400 meters, 0～1 s)

圖13　2400 m處3-4s單道對比(a) Q=0與Q=10單道對比; (b) Q=0與Q=20單道對比; (c) Q=0與Q=40單道對比; (d) Q=0與Q=80單道對比.Fig.13　Contrast of single traces (Distance=2400 meters, 3～4 s)

對于實(shí)際數(shù)據(jù)應(yīng)用而言，本文方法尚需進(jìn)一步探索.若采用二維算法，則需要對實(shí)際數(shù)據(jù)的振幅和相位進(jìn)行校正，因此有必要在此基礎(chǔ)上開展三維算法.同時，三維的大數(shù)據(jù)處理對算法適應(yīng)性提出更高的要求，將GPU和MPI相結(jié)合是一種有效的方式，顯然時域全波形反演方法更加適合.

Almomin A, Biondi B. 2012. Tomographic full waveform inversion: Practical and computationally feasible approach. ∥ 82nd Annual International Meeting. SEG, 1-5, doi: 10.1190/segam2012-0976.1.

Bai J Y, Yingst D, Bloor R, et al. 2014. Viscoacoustic waveform inversion of velocity structures in the time domain.Geophysics, 79(3): R103-R119, doi: 10.1190/geo2013-0030.1. Biondi B, Almomin A. 2012. Tomographic full waveform inversion (TFWI) by combining full waveform inversion with wave-equation migration velocity analysis. ∥ 82nd Annual International Meeting. SEG, 1-5, doi: 10.1190/segam2012-0275.1. Boonyasiriwat C, Zhan G, Hadwiger M, et al. 2010. Multisource reverse-time migration and full-waveform inversion on a GPGPU. ∥ 72nd EAGE Conference and Exhibition. SPE, EAGE. Brenders A J, Pratt R G. 2007. Full-waveform tomography for lithospheric imaging: Results from a blind test in a realistic crustal model.GeophysicalJournalInternational, 168(1): 133-151, doi: 10.1111/j.1365-246X.2006.03156.x.

Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion.Geophysics, 60(5): 1457-1473, doi: 10.1190/1.1443880. Gauthier O, Virieux J, Tarantola A. 1986. Two-dimensional nonlinear inversion of seismic waveforms: Numerical results.Geophysics, 51(7): 1387-1403, doi: 10.1190/1.1442188.

Hager W W, Zhang H C. 2006. A survey of nonlinear conjugate gradient methods.PacificJournalofOptimization, 2(1): 35-58.

Li B, Liu H W, Liu G F, et al. 2010. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU.ChineseJ.Geophys. (in Chinese), 53(12): 2938-2943, doi: 10.3969/j.issn.0001-5733.2010.12.017.

Liu H W, Liu H, Li B, et al. 2011. Pre-stack reverse time migration for rugged topography and GPU acceleration technology.ChineseJ.Geophys. (in Chinese), 54(7): 1883-1892, doi: 10.3969/j.issn.0001-5733.2011.07.022. Liu S W, Wang H Z, Chen S C, et al. 2013. Implementation strategy of 3D reverse time migration on GPU/CPU clusters.ChineseJ.Geophys. (in Chinese), 56(10): 3487-3496, doi: 10.6038/cjg20131023.

Luo J R, Wu R S. 2013. Envelope inversion-Some application issues. ∥ 83rd Annual International Meeting. SEG, 1042-1047, doi: 10.1190/segam2013-0718.1.

Moghaddam P P, Keers H, Herrmann F J, et al. 2013. A new optimization approach for source-encoding full-waveform inversion.Geophysics, 78(3): R125-R132, doi: 10.1190/geo2012-0090.1.

Pratt R G. 1999a. Seismic waveform inversion in the frequency domain, Part I: Theory and verification in a physical scale model.Geophysics, 64(3): 888-901, doi: 10.1190/1.1444597. Pratt R G, Shipp R M. 1999b. Seismic waveform inversion in the frequency domain, Part 2: Fault delineation in sediments using crosshole data.Geophysics, 64(3): 902-914, doi: 10.1190/1.1444598.

Shi Y, Liu H, Zou Z. 2010. Surface-related multiples prediction based on wave equation and adaptive subtraction investigation.ChineseJ.Geophys. (in Chinese), 53(7): 1716-1724, doi: 10.3969/j.issn.0001-5733.2010.07.023.

Shi Y, Wang W H, Li Y, et al. 2013. 3D surface-related multiple prediction approach investigation based on wave equation.ChineseJ.Geophys. (in Chinese), 56(6): 2023-2032, doi: 10.6038/cjg20130623.

Shin C, Jang S, Min D J. 2001. Improved amplitude preservation for prestack depth migration by inverse scattering theory.GeophysicalProspecting, 49 (5): 592-606, doi: 10.1046/j.1365-2478.2001.00279.x.Shin C, Cha Y H. 2008. Waveform inversion in the Laplace domain.GeophysicalJournalInternational, 173(3): 922-931, doi: 10.1111/j.1365-246X.2008.03768.x.

Shin C, Cha Y H. 2009. Waveform inversion in the Laplace-Fourier domain.GeophysicalJournalInternational, 177(3): 1067-1079, doi: 10.1111/j.1365-246X.2009.04102.x. Shin J, Ha W, Jun H, et al. 2014. 3D Laplace-domain full waveform inversion using a single GPU card.Computers&Geosciences, 67: 1-13, doi: 10.1016/j.cageo.2014.02.006.Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation.Geophysics, 49(8): 1259-1266, doi: 10.1190/1.1441754. Tarantola A. 1986. A strategy for nonlinear elastic inversion of seismic reflection data.Geophysics, 51(10): 1893-1903, doi: 10.1190/1.1442046.

Wu R S, Luo J R, Wu B Y. 2013. Ultra-low-frequency information in seismic data and envelop inversion. ∥ 83nd Annual International Meeting. SEG, 3078-3082, doi: 10.1190/segam2013-0825.1.Wu R S, Luo J R, Wu B Y. 2014. Seismic envelope inversion and modulation signal model.Geophysics, 79(3): WA13-WA24, doi: 10.1190/geo2013-0294.1. Yang P L, Gao J H, Wang B L. 2014. RTM using effective boundary saving: A staggered grid GPU implementation.Computers& Geosciences, 68: 64-72,doi: 10.1016/j.cageo.2014.04.004.YangPL,GaoJH,WangBL. 2015.Agraphicsprocessingunitimplementationoftime-domainfull-waveforminversion. Geophysics, 80(3):F31-F39,doi: 10.1190/GEO2014-0283.1.

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

李博, 劉紅偉, 劉國峰等. 2010. 地震疊前逆時偏移算法的CPU/GPU實(shí)施對策. 地球物理學(xué)報, 53(12): 2938-2943, doi: 10.3969/j.issn.0001-5733.2010.12.017.

劉紅偉, 劉洪, 李博等. 2011. 起伏地表疊前逆時偏移理論及GPU加速技術(shù). 地球物理學(xué)報, 54(7): 1883-1892, doi: 10.3969/j.issn.0001-5733.2011.07.022.

劉守偉, 王華忠, 陳生昌等. 2013. 三維逆時偏移GPU/CPU機(jī)群實(shí)現(xiàn)方案研究. 地球物理學(xué)報, 56(10): 3487-3496, doi: 10.6038/cjg20131023.

石穎, 劉洪, 鄒振. 2010. 基于波動方程表面多次波預(yù)測與自適應(yīng)相減方法研究. 地球物理學(xué)報, 53(7): 1716-1724, doi: 10.3969/j.issn.0001-5733.2010.07.023.

石穎, 王維紅, 李瑩等. 2013. 基于波動方程三維表面多次波預(yù)測方法研究. 地球物理學(xué)報, 56(6): 2023-2032, doi: 10.6038/cjg20130623.

(本文編輯胡素芳)

Time domain full waveform inversion based on frequency attenuation

GUO Xue-Bao1,2， LIU Hong1,2， SHI Ying3

1KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029，China2InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029，China3NortheastPetroleumUniversity,Daqing163318，China

Because the full frequency range information is used in time domain full waveform inversion, different wavelengths can not be reconstructed gradually. So it is easy to fall into local minimum value. In this paper, we extract different frequency bands from seismic data by Fourier transform. Different level of high frequency components can be separated by attenuating corresponding frequency ranges, and time domain full waveform inversion is implemented from low to high frequency. At the same time, latter arrivals are attenuated and this can reduce the interference in the process of shallow part inversion. Fourier transform is only required for the shot gathers, compared with hybrid domain inversion, and don′t need to extract the corresponding frequency components of all wavefields. In the aspect of computing efficiency, the algorithm is based on GPU and uses cufft to further improve the efficiency of calculation. Through the Marmouse model test, the validity of this method is verified.

Time domain full waveform inversion; Frequency attenuation; Low frequency; GPU

10.6038/cjg20161022.

國家高技術(shù)研究發(fā)展計劃(863計劃)(2012AA061202)，國家自然科學(xué)基金項(xiàng)目(41574117, 41474118)，黑龍江省杰出青年科學(xué)基金項(xiàng)目(JC2016006)以及大連理工大學(xué)國家重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目(LP1509)聯(lián)合資助.

郭雪豹，1990年生，在讀博士，主要從事地震數(shù)據(jù)成像相關(guān)研究工作.E-mail：guoxuebao1108@163.com

10.6038/cjg20161022

P631

2015-11-16，2016-06-28收修定稿

郭雪豹， 劉洪， 石穎. 2016. 基于頻域衰減的時域全波形反演. 地球物理學(xué)報，59(10):3777-3787,

Guo X B， Liu H， Shi Y. 2016. Time domain full waveform inversion based on frequency attenuation.ChineseJ.Geophys. (in Chinese),59(10):3777-3787,doi:10.6038/cjg20161022.