段 心 標

（①同濟大學海洋與地球科學學院波現(xiàn)象與智能反演成像研究組，上海200092；②中國石化石油物探技術研究院，江蘇南京211103）

0 引言

近年來，隨著計算機性能和存儲能力大幅提高，逆時偏移技術迅速發(fā)展成熟并得到廣泛應用［1-7］。目前最常見逆時偏移技術采用顯式有限差分方法進行正反波場傳播，利用互相關成像條件獲取成像值，進而通過濾波方式去除低頻噪聲［8-9］。這類方式在應用中存在兩個局限：一是實際偏移時計算網格較大，顯式有限差分法存在較嚴重的數(shù)值頻散［10-12］，且產生無規(guī)則噪聲，不利于高頻信息的保護與成像；二是成像后濾波在去除低頻噪聲的同時，會損傷有效低頻信息，不利于低頻信息的保護與成像。

空間差分數(shù)值頻散由地震波的傳播方向、空間差分精度和空間采樣密度等三個因素決定。其中：因地震波向四周傳播，故傳播方向無法控制；提高空間差分精度，會使數(shù)值頻散減弱；空間采樣越密，一個波長范圍內離散點數(shù)越多，數(shù)值頻散越弱。因此，要想減弱有限差分算子的數(shù)值頻散，通?？刹捎酶叩牟罘蛛A數(shù)或減小空間網格尺寸［13-14］，但這將導致計算量的顯著增加。

偽譜法［15-20］利用快速傅里葉變換將波場數(shù)據(jù)由空間域變換到波數(shù)域，避免了空間差分運算，對空間微分算子的逼近使其計算精度可達Nyquist譜精度，幾乎沒有空間數(shù)值頻散。但因其時間方向二階導數(shù)是采用有限差分算子近似，故仍存在時間頻散，當時間采樣間隔較大時，會出現(xiàn)不穩(wěn)定現(xiàn)象。Etgen等［21］提出擬解析法，在波數(shù)—空間域直接修正拉普拉斯算子，補償因時間二階差分近似而造成的解析誤差，相對于偽譜法而言可有效消除時間方向數(shù)值頻散。對于速度變化劇烈的介質，可由若干常速擬拉普拉斯算子插值得到變速擬拉普拉斯算子。為了更好處理變速介質情況下的波場傳播問題，Song等［22］提出傅里葉有限差分法，將速度分解為背景速度和擾動速度，分別利用快速傅里葉變換和有限差分進行波場傳播計算，實現(xiàn)了對波場傳播數(shù)值頻散的有效壓制。

逆時偏移互相關成像條件包含炮點與檢波點傳播方向相同波場的成像計算，而這部分成像與Claerbout成像準則不相符，它是逆時偏移傳播路徑上低頻、強振幅偏移噪聲產生的根本原因［4］。常規(guī)疊后濾波法雖簡單快捷，但會損失成像剖面中的低頻有效信息，影響高陡斷面構造成像質量。

為了消除逆時偏移的低頻噪聲，F(xiàn)ei等［23-24］提出在頻率—波數(shù)域將波場分解為不同方向的傳播波場，再做選擇性互相關成像，從而實現(xiàn)對低頻噪聲和偏移假象的壓制。然而，這種頻率—波數(shù)域波場分解方式需保存所有時刻波場，并沿最慢維的時間方向做傅里葉變換，計算量和I／O 量都很大，實際應用中難以實現(xiàn)。對此，Liu等［25］提出僅在z 方向波數(shù)域實現(xiàn)波場分解的隱式分解法，避免了時間域傅里葉變換導致的計算和存儲問題，有利于逆時偏移中消除低頻 噪 聲。胡 江 濤 等［26］、Shen等［27］提 出 基于解析時間波場的顯式波場分解逆時偏移方法，消除低頻噪聲的同時壓制了偏移構造假象，波場外推過程中需同時求解兩次波動方程，利用原波場與構建的解析時間波場進行波場分解，避免了頻率—波數(shù)域波場分解的存儲問題。此類先波場分解再相關成像的逆時偏移方法，盡管成像條件與單程波偏移有些類似，不利于鹽下構造回折波成像，但能滿足一般復雜構造成像需求，因而受到廣泛關注。

本文針對常規(guī)逆時偏移顯式有限差分波場傳播算法和拉普拉斯濾波方法存在的高低頻信息損失問題，采用傅里葉有限差分算法計算波場傳播，波場外推中直接在時間—波數(shù)域進行顯式波場分解和互相關成像，從而實現(xiàn)對高、低頻信息的充分保護與利用；針對傅里葉有限差分計算效率較低問題，研討了依托GPU 的算法，顯著提高了傅里葉有限差分逆時偏移的計算速度。將該方法應用于實際地震資料處理，其成像效果優(yōu)于常規(guī)逆時偏移。

1 方法原理

1.1 傅里葉有限差分波場傳播算法

逆時偏移雙程波方程為

式中：u（x，t）是空間點x 在時刻t 的地震波場；v是地震波在介質中的傳播速度。

據(jù)Soubaras等［28］提出的時間匹配算法和傅里葉變換公式，假設地下介質為常速，則式（1）可寫為式中：k是波數(shù)矢量；︵u（k，t）表示時間—波數(shù)域的地震波場；Δt是時間步長；vc是假設的常速度值。

實際地下介質是變速的，可用v（x）近似代替式（2）中的vc，并引入背景速度v0，則有

上式右端的分式可展開為

式中：n＝1，2，3代表三個不同的方向；Δxn、kn分別是某一方向的網格間距和波數(shù)；a、bn為系數(shù)，可分別表示為

式（2）等號右端可分為兩步實現(xiàn)：首先基于背景速度利用傅里葉正反變換進行波場傳播；然后再利用有限差分算法校正速度擾動影響。因而，雙程波傅里葉有限差分波場傳播算法實現(xiàn)步驟為：

（1）利用FFT 將時間—空間域波場u（x，t）變換成時間—波數(shù)域波場（k，t）；

（4）利用式（5）、式（6）中差分系數(shù)，通過有限差分法由q（x，t）計算速度擾動下的傳播波場q（x，t＋Δt）；

（5）基于時間差分格式，新時刻的時間—空間域波場u（x，t＋Δt）表示為q（x，t＋Δt）＋2u（x，t）－u（x，t－Δt）；

（6）重復步驟（1）～（5），逐個計算各時刻波場。

利用上述步驟開展雙程波波場傳播試驗。首先給定一個三維均勻介質模型，速度為3000m／s，模型尺寸為9000m×9000m×5000m，速度網格單元為20m×20m×20m，震源點位于模型中心（4500m，4500m，2500m），地震子波是主頻為40Hz的雷克子波，時間步長為1ms。

圖1是波場快照剖面，圖2是其水平切片?？梢姵Ｒ?guī)有限差分法波場傳播過程中存在較為嚴重的頻散噪聲，而傅里葉有限差分法波場快照中，頻散得到較徹底壓制，有利于高頻信息的保護。

圖1 均勻速度模型的常規(guī)有限差分法（a）和傅里葉有限差分法（b）波場快照剖面

圖2 均勻速度模型的常規(guī)有限差分法（a）和傅里葉有限差分法（b）波場快照水平切片

進一步利用三維鹽丘模型檢測傅里葉有限差分法針對復雜介質的處理效果。速度網格單元為30m×30m×20m，考慮到海水速度為1500m／s，地震子波采用主頻為20Hz雷克子波，時間步長為2ms。圖3是鹽丘速度模型及對應的T1時刻常規(guī)有限差分法和傅里葉有限差分法波場快照，圖4是應用這兩種方法得到的T2時刻波場快照的水平切片。對比亦能看出，常規(guī)有限差分法波場頻散較嚴重，而傅里葉有限差分法波場信息更保真。

1.2 時間—波數(shù)域解析波場分解方法

解析波場表示為

其實部是地震傳播波場u（x，t），虛部p（x，t）是波場u（x，t）的Hilbert變換。Hilbert變換是一個時間方向的卷積，直接進行Hilbert變換計算需將所有時刻的波場數(shù)據(jù)存盤并做卷積運算，數(shù)據(jù)量和計算量都很大，因而實用性較差。胡江濤等［26］提出對波動方程的源項進行Hilbert變換的方法。考慮到波場傳播算子的線性特征，傳播波場即是解析時間波場的虛部，如此通過兩個波場傳播構建各時刻的解析波場。另外，傅里葉有限差分算法在波數(shù)域和空間域是交替進行波場傳播，本文采用Zhang 等［29］給出的波數(shù)域Hilbert變換方法

式中：F－1表示三維空間傅里葉反變換；Re表示對復數(shù)取實部。其傳播波場為

圖5是圖1b所示數(shù)據(jù)分離后的波場，可見上、下行波得到很好分離。此外，對三維鹽丘模型的復雜波場也進行了波場分離試驗，圖6a、圖6b對應為圖3c波場分離后的上行波和下行波，分離效果整體較好，但也存在一定殘余?？紤]到實際逆時偏移采用的速度模型通常是低頻光滑模型，透射波能量明顯大于反射波能量，且逆時偏移是震源順時傳播的下行波（透射波）與檢波器逆時傳播的下行波（透射波）相關成像，其上行波殘余可忽略不計。

1.3 GPU 加速方法

疊前偏移成像計算量大、效率低，通常采用GPU 并行算法實現(xiàn)加速［30-32］。傅里葉有限差分波場傳播算法需做大量的正反三維傅里葉變換，計算效率低于常規(guī)高階有限差分算法。本文利用GPU運算平臺CUDA（Compute Unified Device Architecture）提供的用于傅里葉變換的函數(shù)庫cu FFT（CUDA Fast Fourier Transform）進行加速。試算結果表明，cu FFT 對信號長度較敏感，當信號長度為奇數(shù)時，運算速度較慢；當信號長度為偶數(shù)時，速度較快；當信號長度為64的整數(shù)倍時，速度最快。因此，執(zhí)行cuFFT 之前需先對三維空間網格進行鑲邊，使x、y、z方向空間網格數(shù)為64的整數(shù)倍，達到實現(xiàn)快速傅里葉變換的目的。

圖3 鹽丘模型及其T1時刻波場快照剖面

圖4 鹽丘模型的常規(guī)有限差分法（a）和傅里葉有限差分法（b）T2時刻波場快照水平切片

圖5 均勻速度模型的下行波（a）和上行波（b）波場分離結果

圖6 鹽丘模型的上行波（a）和下行波（b）波場分離結果

GPU 加速傅里葉有限差分法有兩種計算方案：一是僅將正反三維FFT 置于GPU 中，有限差分計算由CPU 完成，該方案對內存需求小，但每一步波場傳播都需做GPU 數(shù)據(jù)I／O，影響計算效率；二是將整個計算步驟都置于GPU 中，避免每步波場傳播的數(shù)據(jù)I／O，顯然會增加內存需求，特別是FD 有限差分系數(shù)需占用大量內存空間。對此，在每步波場傳播過程中重復計算FD 有限差分系數(shù)，通過增加計算量換取減小內存空間，以降低內存需求。

同樣采用上面的參數(shù)和硬件資源條件測試不同方法的計算效率，得到的對比結果如表1?？梢娚鲜龇桨付侍嵘^大，達到12 倍加速比。GPU內存耗費方面，方案二需存儲速度場、成像值和兩個時刻的復數(shù)域波場，大致是常規(guī)有限差分法的1.5倍，目前常規(guī)GPU 設備就能滿足此內存需求。本文采用方案二進行GPU 加速計算。

表1 不同算法效率對比

2 應用實例

M 探區(qū)構造較簡單，資料品質高，常規(guī)RTM 成像分辨率較低（圖7a）。應用本文方法進行偏移處理得到成像剖面（圖7b），可見本文方法成像結果分辨率更高，細節(jié)信息更豐富。圖8是兩個結果由深度域轉時間域后的頻譜圖，顯然可見本文方法成像頻帶更寬。

N 探區(qū)構造較復雜，資料品質一般。從常規(guī)RTM 法（圖9a）和本文方法（圖9b）成像剖面看出，本文方法分辨率依然較高；從其成像頻譜對比圖（圖10）可見，本文方法具有更寬頻帶。

圖7 M 探區(qū)資料常規(guī)RTM（a）和本文方法（b）成像剖面

圖8 M 探區(qū)資料常規(guī)RTM （藍線）和本文方法（紅線）成像剖面的頻譜

3 結論

圖9 N 探區(qū)資料常規(guī)RTM（a）和本文方法（b）成像剖面

圖10 N 探區(qū)資料常規(guī)RTM （藍線）和本文方法（紅線）成像剖面的頻譜

針對常規(guī)逆時偏移中高階有限差分算法、疊后濾波算法中的不足，本文研究了基于GPU 的傅里葉有限差分波場傳播算法和時間—波數(shù)域顯式波場分解方法，形成了基于GPU 的傅里葉有限差分逆時偏移成像技術，得到以下結論和認識。

（1）逆時偏移的高頻信息損失緣于高階有限差分空間頻散，傅里葉有限差分波場傳播方法將速度分解為背景速度和擾動速度，分別利用快速傅里葉變換和有限差分計算波場傳播，能較好壓制數(shù)值頻散，更好地保護高頻數(shù)據(jù)信息。

（2）基于解析波場的上、下行波分解是逆時偏移壓制低頻噪聲的重要手段，在傅里葉有限差分波場傳播中，直接對時間—波數(shù)域波場做波場分解，實現(xiàn)簡單，額外計算量小。

（3）傅里葉有限差分波場傳播算法中包含大量的FFT 計算，效率較低，通過GPU 平臺cuFFT 函數(shù)加速并選擇合適的I／O 策略，能顯著提高傅里葉有限差分法計算效率，使其實用性增強。

（4）實際資料應用結果表明，該技術有利于高低頻信息的充分保護與利用，提高了地震成像分辨率，成像效果優(yōu)于常規(guī)逆時偏移。