郭雪豹， 劉洪， 石穎

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

?

基于頻域衰減的時(shí)域全波形反演

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

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

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

時(shí)域全波形反演；頻域衰減；低頻；GPU

1　引言

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

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

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

2　基本原理

2.1時(shí)域全波形反演基本理論

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

(1)

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

(2)

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

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

(3)

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

參數(shù)更新公式：

(4)

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

(5)

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

(6)

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

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)為衰減后的波場(chǎng)數(shù)據(jù)殘差.

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

(8)

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

2.3GPU_CUFFT

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

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

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

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

3　模型測(cè)試

本文采用Marmousi模型測(cè)試所述方法，Marmousi模型如圖1所示，縱、橫向網(wǎng)格點(diǎn)數(shù)分別為251和767，網(wǎng)格間距為12 m，雷克子波主頻為10 Hz，時(shí)間采樣間隔1.5 ms，采樣點(diǎn)數(shù)為2800.設(shè)計(jì)21炮震源激發(fā)，第一炮位于深度36 m，自左向右36 m位置處，向右每隔456 m放一炮.每炮767個(gè)檢波器接收，檢波器深度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)定.對(duì)比速度模型可以看出，其反演結(jié)果并不理想，尤其在淺部存在較強(qiáng)的噪音干擾.由此可以看出，即便采用平滑后的真實(shí)速度模型，時(shí)域的全頻帶反演沒有逐步的反演出不同的波長(zhǎng)分量，依舊會(huì)使結(jié)果陷入局部極小值.

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

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

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

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

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

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

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

4　結(jié)論

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

圖7　4800 m處單道對(duì)比(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處單道對(duì)比(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ù)對(duì)比(a) Q=10, Q=20, Q=40單道對(duì)比; (b) Q=20, Q=40, Q=80單道對(duì)比; (c) Q=80, Q=0, Q=0單道對(duì)比.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處兩種方法單道對(duì)比(a) 原始方法迭代300次與Q=10迭代40次對(duì)比; (b) 原始方法迭代300次與Q=20迭代40次對(duì)比; (c) 原始方法迭代300次與Q=40迭代40次對(duì)比.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處兩種方法單道對(duì)比(a) 原始方法迭代300次與Q=80迭代40次對(duì)比; (b) 原始方法迭代300次與Q=0迭代40次對(duì)比; (c) 原始方法迭代300次與Q=0迭代140次對(duì)比.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單道對(duì)比(a) Q=0與Q=10單道對(duì)比; (b) Q=0與Q=20單道對(duì)比; (c) Q=0與Q=40單道對(duì)比; (d) Q=0與Q=80單道對(duì)比.Fig.12　Contrast of single traces (Distance=2400 meters, 0～1 s)

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

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

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)

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

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

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

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

石穎, 王維紅, 李瑩等. 2013. 基于波動(dòng)方程三維表面多次波預(yù)測(cè)方法研究. 地球物理學(xué)報(bào), 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.

國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)(2012AA061202)，國(guó)家自然科學(xué)基金項(xiàng)目(41574117, 41474118)，黑龍江省杰出青年科學(xué)基金項(xiàng)目(JC2016006)以及大連理工大學(xué)國(guó)家重點(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. 基于頻域衰減的時(shí)域全波形反演. 地球物理學(xué)報(bào)，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.