基于反泄漏Fourier變換的地震數據重建方法

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

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

?

基于反泄漏Fourier變換的地震數據重建方法

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

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

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

數據重建； 規則化； Fourier變換； GPU； 地震數據采集

0 引言

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

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

CPU/GPU協同并行加速技術被廣泛應用于地震處理的各個領域。如在多次波預測中，通過應用GPU加速處理矩陣乘法部分[9]，可以大幅提高處理速度。在疊前偏移[10～13]及全波形反演中[14]，GPU加速技術的應用使得數據的處理變得更為方便快捷。基于反泄露Fourier變換方法，筆者應用CPU/GPU并行加速算法提高運算效率，并驗證方法的合理性和有效性。

1 方法原理

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

(1)

(2)

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

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

去除表達式為

(3)

規則化重建表達式為

(4)

2 優化反泄漏算法

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

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

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

(5)

(6)

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

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

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

3 理論模型與實際數據測試

3.1 理論模型測試

為了測試文中算法的有效性，抽取多炮復雜Smaart模型單炮并進行模型試算，去除單炮數據部分道信息并進行抽稀處理，將數據變為非均勻采樣形式(見圖2)。該數據空間方向為361道，道間距為75 m，時間方向采樣率為0.008 s。空間方向上隨機抽出22道，最大道間距為225 m(見圖2(a))，缺失道重建后成圖(見圖2(b))。其中不規則數據f-k譜見圖2(c)，規則化后數據f-k譜見圖2(d)。

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

圖2 Smaart模型理論數據測試結果Fig.2 Smaart model theoretical data test results

3.2 實際數據測試

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

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

圖3 某油田海上疊前實際數據測試結果Fig.3 Test results for real data of oilfield

3.3 數據加速比

為測試基于CUDA的反泄漏地震數據重建算法的性能，在LINUX操作平臺上應用C語言、結合CUDA語言編寫程序，并進行測試。算法測試環境硬件參數見表1。

表1 算法測試環境硬件參數

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

表2 算法CPU和CPU/GPU程序測試運行時間

4 結論

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

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

[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] 孟小紅,郭良輝,張致付,等.基于非均勻快速傅里葉變換的最小二乘反演地震數據重建[J].地球物理學報,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] 劉喜武,劉洪,劉彬.反假頻非均勻地震數據重建方法研究[J].地球物理學報,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] 高建軍,陳小宏,李景葉.三維不規則地震數據重建方法[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] 石穎,王維紅,李瑩,等.基于波動方程三維表面多次波預測方法研究[J].地球物理學報,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數據逆時偏移[J].東北石油大學學報,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].東北石油大學學報,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；編輯：任志平

黑龍江省自然科學基金面上項目(D2015011)

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

王維紅,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