王天野， 白軍輝， 王維紅， 賀旎妮， 宋元東， 井洪亮

( 1. 東北石油大學(xué) 地球科學(xué)學(xué)院，黑龍江 大慶 163318； 2. 大慶油田有限責(zé)任公司 第三采油廠，黑龍江 大慶 163113； 3. 北京大學(xué) 地球與空間科學(xué)學(xué)院，北京 100871； 4. 東方地球物理公司 吐哈物探處，新疆 哈密 839009 )

?

基于反泄漏Fourier變換的地震數(shù)據(jù)重建方法

王天野1， 白軍輝2， 王維紅1， 賀旎妮3， 宋元東1， 井洪亮4

( 1. 東北石油大學(xué) 地球科學(xué)學(xué)院，黑龍江 大慶 163318； 2. 大慶油田有限責(zé)任公司 第三采油廠，黑龍江 大慶 163113； 3. 北京大學(xué) 地球與空間科學(xué)學(xué)院，北京 100871； 4. 東方地球物理公司 吐哈物探處，新疆 哈密 839009 )

應(yīng)用反泄漏Fourier重建方法，對空間方向采樣不規(guī)則的地震數(shù)據(jù)進(jìn)行規(guī)則化重建，引入GPU加速算法，將非均勻Fourier變換部分轉(zhuǎn)化為矩陣乘法形式，將矩陣乘法部分應(yīng)用于GPU加速算法加速處理。理論模型與實(shí)際數(shù)據(jù)測試結(jié)果驗(yàn)證該方法的有效性及合理性。

數(shù)據(jù)重建； 規(guī)則化； Fourier變換； GPU； 地震數(shù)據(jù)采集

0 引言

在地震數(shù)據(jù)采集過程中，通常受障礙物、禁采區(qū)、海洋拖纜羽狀漂移影響及經(jīng)濟(jì)條件等制約。地震數(shù)據(jù)在空間方向上的采集通常是稀疏的或不規(guī)則的，對后續(xù)偏移及多次波處理等產(chǎn)生影響。因此，對地震數(shù)據(jù)進(jìn)行合理化重建、規(guī)則地震數(shù)據(jù)及合理恢復(fù)稀疏道尤為重要。

目前，基于變換域重建的地震數(shù)據(jù)重建方法分為Radon變換重建、Fourier變換重建及Curvelet變換重建。這類方法一般將數(shù)據(jù)進(jìn)行正變換，在變換域中根據(jù)信號特點(diǎn)進(jìn)行相應(yīng)處理；再通過反變換將數(shù)據(jù)變換到原始數(shù)據(jù)域，通過正反變換實(shí)現(xiàn)地震數(shù)據(jù)重建。第一步，將數(shù)據(jù)以恰當(dāng)?shù)男问秸儞Q估算得到變換域的系數(shù)。假使數(shù)據(jù)是不規(guī)則的，那么正變換計算得到的系數(shù)可能是不精確的。第二步，通過反變換合理重建理想規(guī)則數(shù)據(jù)。基于變換域的重建方法的優(yōu)勢是計算速度快，對輸入數(shù)據(jù)可以在一定程度上稀疏采樣，既可以較好地處理規(guī)則采樣數(shù)據(jù)，又可以對非規(guī)則采樣數(shù)據(jù)進(jìn)行合理重建。當(dāng)?shù)卣鹗怯邢迬挄r重建效果更好。Fourier變換重建技術(shù)已經(jīng)成功應(yīng)用于一維、二維非規(guī)則采樣數(shù)據(jù)重建。Duijndam A J W等[1]采用分頻譜反演方法估計空間Fourier系數(shù)，并成功將Fourier重建方法應(yīng)用于非規(guī)則采樣地震數(shù)據(jù)的規(guī)則化問題和三維數(shù)據(jù)重建[2]，可以解決參數(shù)選擇等問題，但存在空間假頻的地震數(shù)據(jù)而使重建效果變差。Liu B等[3]提出最小加權(quán)范數(shù)插值的Fourier重建方法，雖然可以對所得解進(jìn)行優(yōu)化，但沒有解決數(shù)據(jù)存在假頻的情況。孟小紅[4]等將帶限地震數(shù)據(jù)的重建歸結(jié)為最小二乘反演問題，并將非均勻離散Fourier變換變?yōu)榉蔷鶆蚩焖貴ourier變換[5]，有效提高計算效率，但抗假頻能力不足。劉喜武等將非均勻Fourier變換與預(yù)測濾波方法結(jié)合，以解決抗假頻重建問題[6]。高建軍[7]等應(yīng)用抗泄漏Fourier變換方法，有效解決不規(guī)則地震道假頻問題，并結(jié)合貝葉斯反演對三維地震數(shù)據(jù)進(jìn)行重建。Xu S[8]等應(yīng)用過采樣技術(shù)及權(quán)值求取方法，消除吉布斯現(xiàn)象及假頻影響，并給出二維切片形式，以重建三維數(shù)據(jù)體。

CPU/GPU協(xié)同并行加速技術(shù)被廣泛應(yīng)用于地震處理的各個領(lǐng)域。如在多次波預(yù)測中，通過應(yīng)用GPU加速處理矩陣乘法部分[9]，可以大幅提高處理速度。在疊前偏移[10～13]及全波形反演中[14]，GPU加速技術(shù)的應(yīng)用使得數(shù)據(jù)的處理變得更為方便快捷?；诜葱孤禙ourier變換方法，筆者應(yīng)用CPU/GPU并行加速算法提高運(yùn)算效率，并驗(yàn)證方法的合理性和有效性。

1 方法原理

反泄漏Fourier變換的目的是對不規(guī)則地震數(shù)據(jù)進(jìn)行重建，并能夠壓制波數(shù)域譜泄漏。由于時間方向是規(guī)則采樣，因此可以應(yīng)用FFT將時間域數(shù)據(jù)變換到頻率域數(shù)據(jù)，以提高計算效率。在不規(guī)則網(wǎng)格中，單個頻率的反泄漏Fourier變換公式[11]可以表示為

(1)

(2)

式(1-2)中：F(k)為波數(shù)k的Fourier系數(shù)；xn為采樣點(diǎn)，ΔX=∑Δxn；Np為輸入點(diǎn)的個數(shù)；fk(xl)為k的Fourier系數(shù)對采樣點(diǎn)xn的函數(shù)值貢獻(xiàn)。對于不規(guī)則數(shù)據(jù)，由于正交性被破壞，fk(xl)能量泄漏到其他頻率成分中。

為了消除譜泄漏現(xiàn)象，首先對所有Fourier系數(shù)進(jìn)行估計，確定最大能量及其對應(yīng)位置；然后反變換到輸入數(shù)據(jù)網(wǎng)格中，從原始數(shù)據(jù)中去除最大能量對應(yīng)的貢獻(xiàn)值，得到更新后的輸入數(shù)據(jù)并反復(fù)迭代；最終得到消除譜泄漏的規(guī)則化數(shù)據(jù)。

去除表達(dá)式為

(3)

規(guī)則化重建表達(dá)式為

(4)

2 優(yōu)化反泄漏算法

將反泄漏Fourier變換算子及輸入數(shù)據(jù)拆分成矩陣相乘形式，并從主機(jī)端輸入到GPU端，根據(jù)分塊線程原則和矩陣乘法在GPU中的計算方法，合理分配內(nèi)存及賦值處理。該算法原理見圖1。

圖1 GPU棋盤劃分矩陣乘法流程示意Fig.1 GPU board division matrix multiplication process

將反泄漏變換公式拆分，正反變換矩陣相乘形式為

(5)

(6)

式(5-6)中：offset[i]為空間方向不規(guī)則采樣間隔；offsetinter[i]為規(guī)則化后空間方向的采樣間隔；dk為波數(shù)方向的采樣間隔。

在科學(xué)計算方面，GPU具有高速度、高性能及較高的并行計算效率。利用GPU加速計算NDFT運(yùn)算的矩陣乘法、GPU提供的CUFFT(基于CUDA語言的快速Fourier變換)函數(shù)庫，以及乘除開方等快速計算函數(shù)優(yōu)化算法指令，提高地震數(shù)據(jù)重建算法效率。

由圖1可以看出，該方法對全局存儲器的數(shù)據(jù)依次按照線程ID號，將A塊數(shù)據(jù)依次放入BLOCK對應(yīng)的共享存儲器；然后BLOCK中每個Thread按照索引計算一個Csub，進(jìn)行多次循環(huán)，每次循環(huán)后將數(shù)據(jù)加到上次循環(huán)得到的Csub上，直到所有數(shù)據(jù)在共享存儲器中得到計算并將獲得結(jié)果相加；最后得到的求和即為(0,0)塊結(jié)果。

3 理論模型與實(shí)際數(shù)據(jù)測試

3.1 理論模型測試

為了測試文中算法的有效性，抽取多炮復(fù)雜Smaart模型單炮并進(jìn)行模型試算，去除單炮數(shù)據(jù)部分道信息并進(jìn)行抽稀處理，將數(shù)據(jù)變?yōu)榉蔷鶆虿蓸有问?見圖2)。該數(shù)據(jù)空間方向?yàn)?61道，道間距為75 m，時間方向采樣率為0.008 s。空間方向上隨機(jī)抽出22道，最大道間距為225 m(見圖2(a))，缺失道重建后成圖(見圖2(b))。其中不規(guī)則數(shù)據(jù)f-k譜見圖2(c)，規(guī)則化后數(shù)據(jù)f-k譜見圖2(d)。

由圖2(a)、(b)可以看出，經(jīng)過文中算法處理后，圖像缺失道數(shù)據(jù)的缺失信息得到有效補(bǔ)償，重建后圖像清晰、誤差小，補(bǔ)償后圖像光滑連續(xù)，符合實(shí)際數(shù)據(jù)圖形走向。由圖2(c)、(d)可以看出，不規(guī)則缺失數(shù)據(jù)造成的能量泄露現(xiàn)象得到有效壓制，能量收斂，波數(shù)方向0.3～0.4 m-1部分能量回歸到有效能量中。規(guī)則化后f-k譜有效能量集中在小視速度范圍內(nèi)。因此，該算法對于復(fù)雜數(shù)據(jù)具有很好的恢復(fù)效果。

圖2 Smaart模型理論數(shù)據(jù)測試結(jié)果Fig.2 Smaart model theoretical data test results

3.2 實(shí)際數(shù)據(jù)測試

為了驗(yàn)證該算法對復(fù)雜數(shù)據(jù)的適用性，對某油田海上疊前實(shí)際數(shù)據(jù)進(jìn)行處理(見圖3)。該數(shù)據(jù)道間距為20 m，共計800道。由于數(shù)據(jù)道數(shù)較大，對數(shù)據(jù)進(jìn)行去除和前后統(tǒng)一處理，保留前200道。在規(guī)則的200道數(shù)據(jù)中隨機(jī)抽出40道數(shù)據(jù)，其中最大道間距為60 m，抽稀后不規(guī)則數(shù)據(jù)見圖3(a)。對去除后數(shù)據(jù)進(jìn)行規(guī)則化處理，缺失道數(shù)據(jù)得到有效補(bǔ)償，重建后規(guī)則數(shù)據(jù)見圖3(b)。重建前f-k譜見圖3(c)，重建后f-k譜見圖3(d)。

由圖3可以看出，應(yīng)用文中算法對復(fù)雜疊前不規(guī)則數(shù)據(jù)進(jìn)行規(guī)則化重建，重建前、后時域數(shù)據(jù)中缺失道信息得到有效補(bǔ)償，恢復(fù)后圖像同相軸連續(xù)，與原始圖像對比誤差較小。f-k譜泄漏能量收斂，視速度較大部分的泄漏能量回歸到有效能量中，有效部分能量集中在視速度較小位置。

圖3 某油田海上疊前實(shí)際數(shù)據(jù)測試結(jié)果Fig.3 Test results for real data of oilfield

3.3 數(shù)據(jù)加速比

為測試基于CUDA的反泄漏地震數(shù)據(jù)重建算法的性能，在LINUX操作平臺上應(yīng)用C語言、結(jié)合CUDA語言編寫程序，并進(jìn)行測試。算法測試環(huán)境硬件參數(shù)見表1。

表1 算法測試環(huán)境硬件參數(shù)

采用文中算法分別對Smaart模型、疊前實(shí)際數(shù)據(jù)進(jìn)行測試，得到程序運(yùn)行加速比，將運(yùn)行平均時間作為最后的運(yùn)算時間，將最后的運(yùn)算時間作為性能分析的依據(jù)。反泄漏離散Fourier變換部分分別應(yīng)用CPU及CPU/GPU并行處理，兩種數(shù)據(jù)在CPU上實(shí)現(xiàn)的總計算時間與在GPU上實(shí)現(xiàn)的總計算時間見表2，其中CPU運(yùn)行時間與GPU運(yùn)行時間之比為加速比。由表2可以看出，基于CUDA的反泄漏并行算法的有效性，程序在原基礎(chǔ)上得到20～40倍的性能提升，數(shù)據(jù)處理的加速效果明顯，在共享存儲器中計算還可以將算法效率提升1.5倍。

表2 算法CPU和CPU/GPU程序測試運(yùn)行時間

4 結(jié)論

(1)基于反泄漏Fourier變換的疊前地震數(shù)據(jù)規(guī)則化算法，可以對非規(guī)則網(wǎng)格采樣數(shù)據(jù)進(jìn)行規(guī)則化處理，對采樣缺失數(shù)據(jù)進(jìn)行插值重建，從而滿足后續(xù)算法對規(guī)則網(wǎng)格數(shù)據(jù)的要求。應(yīng)用CPU/GPU協(xié)同并行加速算法可以對不規(guī)則數(shù)據(jù)進(jìn)行加速處理，對不規(guī)則及稀疏數(shù)據(jù)進(jìn)行有效重建。該算法準(zhǔn)確、易于實(shí)現(xiàn)，重建效果誤差小、精度高，能夠滿足工業(yè)化處理需求。

(2)通過引用CPU/GPU協(xié)同并行加速算法，該算法將CPU中反復(fù)處理的乘法加法過程拆分成多線程運(yùn)算形式并進(jìn)行處理，在運(yùn)算過程中大幅減少算法的時間復(fù)雜度，在大數(shù)據(jù)及三維數(shù)據(jù)中存在實(shí)現(xiàn)的可能性。

[1] Duijndam A J W, Schonewille M A, Hindriks C O H. Reconstruction of band-limited signals, irregularly sampled along one spatial direction [J]. Geophysics, 1999,64:524-538.

[2] Hindriks K, Ajw. D. Reconstruction of 3D seismic signals irregularly sampled along two spatial coordinates [J]. Geophysics, 2000,65(1):253-263.

[3] Liu B, Sacchi M D. Minimum weighted norm interpolation of seismic records [J]. Geophysics, 2004,69(10):1560-1568.

[4] 孟小紅,郭良輝,張致付,等.基于非均勻快速傅里葉變換的最小二乘反演地震數(shù)據(jù)重建[J].地球物理學(xué)報,2008,51(1):235-241.

Meng Xiaohong, Guo Lianghui, Zhang Zhifu, et al. The nonuniform fast Fourier transform least square inversion of seismic data reconstruction based on [J]. Chinese Journal of Geophysics, 2008,51(1):235-241.

[5] Duijndam A J W, Schonewille M A. Nonuniform fast Fourier transform [J]. Geophysics, 1997,64(2):551.

[6] 劉喜武,劉洪,劉彬.反假頻非均勻地震數(shù)據(jù)重建方法研究[J].地球物理學(xué)報,2004,47(2):299-305.

Liu Xiwu, Liu Hong, Liu Bin. Study on the reconstruction method of the non uniform seismic data in the anti false frequency [J]. Journal of Geophysics, 2004,47(2):299-305.

[7] 高建軍,陳小宏,李景葉.三維不規(guī)則地震數(shù)據(jù)重建方法[J].石油地球物理勘探,2011,46(1):40-47.

Gao Jianjun, Chen Xiaohong, Li Jingye. Seismic data reconstruction method for 3D seismic data [J]. Oil Geophysical Prospecting, 2011,46(1):40-47.

[8] Xu S, Zhang Y, Lambaré G. Antileakage Fourier transform for seismic data regularization in higher dimensions [J]. Geophysics, 2010,75(4):87.

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

Shi Ying, Wang Weihong, Li Ying, et al. Study on multi wave prediction method for 3D surface based on wave equation [J]. Journal of Geophysics, 2013,56(6):2023-2032.

[10] Clapp R G, Fu H, Lindtjorn O. Selecting the right hardware for reverse time migration [J]. Leading Edge, 2010,29(29):48-58.

[11] Knibbe H, Mulder W A, Oosterlee C W, et al. Closing the performance gap between an iterative frequency-domain solver and an explicit time-domain scheme for 3D migration on parallel architectures [J]. Geophysics, 2014,79(2):47-61.

[12] 郭雪豹,王建民,王維紅,等.基于GPU并行加速的VSP數(shù)據(jù)逆時偏移[J].東北石油大學(xué)學(xué)報,2014,38(2):58-62.

Guo Xuebao, Wang Jianmin, Wang Weihong, et al. Reverse time migration of VSP data based on GPU parallel acceleration [J]. Journal of Northeast Petroleum University, 2014,38 (2):58-62.

[13] 田東升,王云專,李義鵬,等.單程和雙程波動方程疊前深度偏移方法[J].東北石油大學(xué)學(xué)報,2014,38(4):39-44.

Tian Dongsheng, Wang Yunzhuan, Li Yipeng, et al. Journal of one-way and two-way wave equation prestack depth migration method [J]. Northeast Petroleum University, 2014,38(4):39-44.

[14] Yang P, Gao J, Wang B. A graphics processing unit implementation of time-domain full-waveform inversion [J]. Geophysics, 2015,80(3):31-39.

2016-07-07；編輯：任志平

黑龍江省自然科學(xué)基金面上項(xiàng)目(D2015011)

王天野(1992-)，男，碩士研究生，主要從事地震資料數(shù)字處理方面的研究。

王維紅,E-mail: wwhsy@sina.com

P631.4+14

A

2095-4107(2016)05-0102-06

DOI 10.3969/j.issn.2095-4107.2016.05.012