基于隨機采樣的頻率域多路徑波場模擬與偏移成像

李 博

（1．中國石油化工股份有限公司石油物探技術(shù)研究院，江蘇南京211103；2．中國科學(xué)院地質(zhì)與地球物理研究所油氣資源研究重點實驗室，北京100029）

基于隨機采樣的頻率域多路徑波場模擬與偏移成像

李 博1，2

（1．中國石油化工股份有限公司石油物探技術(shù)研究院，江蘇南京211103；2．中國科學(xué)院地質(zhì)與地球物理研究所油氣資源研究重點實驗室，北京100029）

采用射線追蹤方法的常規(guī)Kirchhoff深度偏移不能全面描述焦散區(qū)的波現(xiàn)象，有可能導(dǎo)致偏移成像產(chǎn)生畸變，利用此類方法進行復(fù)雜介質(zhì)條件下的偏移成像，不可避免地存在焦散問題。詳細推導(dǎo)了二次震源的波場模擬理論方法與區(qū)域分解的實施方案，在確保局部空間不存在射線交叉的前提下進行區(qū)域分解；基于惠更斯二次震源理論進行波場外推，逐步完成整個區(qū)域的波場延拓。提出了基于隨機稀疏采樣與低秩分解的波場模擬高效算法，借助于GPU的計算能力實現(xiàn)了頻率域多路徑波場模擬與成像，大幅提高了計算效率。數(shù)模實驗結(jié)果表明，該方法能正確處理焦散問題，實現(xiàn)多路徑情況下反射界面的正確成像。

偏移成像；波場外推；焦散；格林函數(shù)；低秩分解；隨機采樣

Kirchhoff積分法深度偏移通常采用的旅行時是單值走時，即一個炮點針對地下的成像點只能有一條傳播路徑，復(fù)雜介質(zhì)條件下，這種方法很難滿足實際需求。若考慮多路徑問題，積分法偏移的理論會遇到焦散問題。需要引入相空間的概念才能區(qū)分多條路徑的波場。很多學(xué)者對此問題進行了深入研究，發(fā)現(xiàn)在計算高頻近似下的體波傳播算法中，關(guān)鍵是要產(chǎn)生空間局部的方向波場才能從根本上解決焦散問題［13］。HILLN［4］和HALE［5］提出高斯束偏移方法在炮集中對不同檢波點進行加窗局部傾斜疊加，利用炮集的τ－p譜進行偏移成像，可以完全解決焦散問題。后來GAO等［6］和SUN等［7］通過選擇τ－p譜的部分能量進行偏移，實現(xiàn)有選擇性的波場傳播的快速束和控制束偏移。李輝等［8］提出了地震數(shù)據(jù)的特征高斯波包（CGP）分解方法，在Radon域選取出地震數(shù)據(jù)的時空和方向特征，大幅減少了高斯波包框架的數(shù)量，實現(xiàn)了方向波場在Gabor域的稀疏化表達。王雄文等［9］將常規(guī)的線性Radon變換（LRT）轉(zhuǎn)化為壓縮感知問題，在壓縮感知的理論框架下實現(xiàn)了高分辨率的方向波分解。這些方法的核心都是如何正確描述空間局部方向波場疊加效應(yīng)的近似表達，以提高波場方向分辨率、能量分解更準(zhǔn)確、特征更突出為研究目標(biāo)。

本文提出頻率域多路徑疊加的波場模擬及成像方法，無需方向波分解傳播，也能解決焦散問題。首先對震源局部范圍內(nèi)無焦散區(qū)域進行波場模擬，然后將已知波場作為二次震源逐步進行波場外推。實際上，在空間同一個位置完成多次射線的疊加，既避免了射線穿過焦散區(qū)域，又實現(xiàn)了多路徑的波場模擬。實際測試結(jié)果表明本文方法可以有效解決焦散問題，獲得正確的成像結(jié)果。

1 頻率域多路徑波場模擬基本原理

基于惠更斯二次震源原理，波場能夠由空間中的一個閉合曲面上的二次震源進行重構(gòu)，其重構(gòu)的表達式為Kirchhoff積分公式：

式中：S表示二次震源的閉合曲面（如圖1所示）；r′表示曲面S上的點；U（r′；r0）表示曲面S上的波場；n（r′）表示二次震源r′處曲面S 的單位法向量；G（r′；r）表示射線從r到r′的Green函數(shù)。

公式（1）經(jīng)常被應(yīng)用在散射波場描述和成像中。在均勻介質(zhì)中，Green函數(shù)G（r′；r）是柱形或球形貝塞爾函數(shù)組成的顯示解析表達形式［10－12］，并且滿足平移和旋轉(zhuǎn)不變性。這種類型的波場模擬快速算法已經(jīng)成功應(yīng)用于各個領(lǐng)域，如快速多極算法［13－17］等，實際上，這些算法的本質(zhì)是如何提高Kirchhoff積分核的矩陣運算。對于非均勻介質(zhì)而言，Green函數(shù)G（r′；r）目前沒有解析表達式，且不滿足平移和旋轉(zhuǎn)不變性，因此均勻介質(zhì)中的快速算法不能簡單地推廣到非均勻介質(zhì)的情況。

本文提出一種非均勻介質(zhì)中Green函數(shù)與Kirchhoff積分數(shù)值計算方法，能夠快速實現(xiàn)公式（1）的波場計算。

在公式（1）中，U（r′；r0）和是二次震源曲面上的波場函數(shù)，首先我們必須計算G（r′；r）和在曲面S 上的值。根據(jù)幾何光學(xué)的高頻漸進假設(shè)，G（r′；r）和有如下形式：

其中，τ表示地震波旅行時。若僅僅保留頻率ω的高階項，忽略低階項可得：

將公式（4）代入公式（1）得到波場的積分近似表達：

由于τ（r′；r）滿足程函方程，代入公式（5）可得：

式中：θ（r′；r）是從r到r′的射線入射方向與曲面S上r′處的單位法向量的夾角（見圖1）。

圖1 惠更斯二次震源波場構(gòu)建原理

從公式（5）和公式（6）可見，需要首先計算G（r′；r）和cos［θ（r′；r）］才能最終完成波場的計算。分析公式（5）需要r到r′的走時τ、振幅A和入射角θ（r′；r），射線的起點應(yīng)該從二次震源的r′開始，即射線從r′發(fā)出到r。令θ（r′；r）＝π－θ＊（r；r′），那么射線從r′到r的出射角為θ＊（r；r′），根據(jù)震源與接收點的互易原理，有G（r′；r）＝G（r；r′），可以得出：

代入公式（5）也可獲得波場的積分表達式：

可見，公式（8）中，利用二次震源曲面上的點r′處的波場U（r′；r0），Green函數(shù)G（r；r′），法向梯度Δr′U（r′；r0）·n（r′）和出射角cos［θ＊（r；r′）］就可以實現(xiàn)波場積分核函數(shù)。公式（8）給出了通過已知波場進行外推的方法，只要保證局部范圍內(nèi)無“焦散”現(xiàn)象，就可以通過射線追蹤、高頻近似等常規(guī)方法進行計算。

2 基于低秩特征分解的惠更斯掃描積分法

理論上，我們需要在無限大的平面上進行二次震源的積分，而實際上只能在有限范圍內(nèi)的數(shù)據(jù)上實現(xiàn)。根據(jù)局部射線無“焦散”問題，在沒有射線交叉或多次到達的情況下，該物理過程滿足Sommerfeld輻射邊界條件，因此截斷二次震源的無限大平面會對截斷位置的波場產(chǎn)生邊界效應(yīng)，但對中間部分的波場影響很小，因此本文方法仍然適用于有限范圍二次震源的情況。下面介紹如何設(shè)定網(wǎng)格參數(shù)并進行波場計算。

給定速度v（r）和頻率ω，我們可以估計最小波長λmin＝2πvmin／ω，這里vmin表示區(qū)域內(nèi)的最小速度，根據(jù)區(qū)域的范圍和最小波長λmin確定采樣點數(shù)，保證每個波長范圍內(nèi)有8～10個采樣點即可。本文方法基于幾何光學(xué)理論，相關(guān)的部分內(nèi)容都可以在稀疏網(wǎng)格上進行，當(dāng)涉及到與頻率相關(guān)的計算時才需要在全部的網(wǎng)格上進行。

離散化的波場積分公式（8）可以表達為如下形式：

式中：Ms表示二次震源的個數(shù)，用表示二次震源集合，用表示需要計算波場的區(qū)域內(nèi)網(wǎng)格點；h表示網(wǎng)格間距 。如圖2所示的密網(wǎng)格部分是需要輸出波場值的區(qū)域，若直接計算波場值將面臨一個O（N3）或O（N5）計算的難題，計算效率難以滿足實際應(yīng)用。

我們提出一種基于低秩特征分解的Kirchhoff積分算法，可以大幅度提高計算效率并借助GPU的計算能力實現(xiàn)波場的快速計算。算法的實現(xiàn)過程如下。

波場計算的公式（9）實際上可以簡化為大型矩陣的乘法和加法，為此，首先做如下變量定義：

其中，f，g，U 為列向量，ρ和ρ1為矩陣。這樣公式（9）可以簡化為：

從公式（10）可以看出，波場計算轉(zhuǎn)換為兩個矩陣向量積之和。每個矩陣向量積都是大規(guī)模的計算，需要耗費大量的計算機時。因此，我們將問題轉(zhuǎn)化為矩陣向量積的快速算法問題。

在光滑背景速度假設(shè)條件下，局部區(qū)域內(nèi)的Green函數(shù)具有一定的冗余性，即ρ＝［G（ri；sj）］應(yīng)該是一個低秩的稠密矩陣。由于cos［θ＊（ri；sj）］也是一個光滑的可壓縮函數(shù)，因此ρ1應(yīng)該也是低秩稠密矩陣。下面僅討論公式（10）的第1項，第2項用同樣的方法可以實現(xiàn)快速計算。

公式（10）的第1項可表達為：

其中，Green函數(shù)矩陣［G（ri；sj）］（i＝1，2，…，n；j＝1，2，…，m）是由旅行時表、振幅表和出射角余弦表計算得到；［fj］（j＝1，2，…，m）由二次震源位置的已知波場計算得到。我們的主要任務(wù)是尋找一種近似關(guān)系滿足如下的條件：

式中：J和I分別表示矩陣G的行和列。如果公式（12）成立，那么不難得出：

式中：W是一個｜J｜×｜I｜的矩陣，由于G是一個低秩矩陣，我們可以找到部分行和列參與計算，縮?。麶｜和｜I｜的長度，從而能夠減少計算量。如果能夠找到矩陣W，那么可以利用公式（13）在精度范圍內(nèi)近似計算出結(jié)果。

為了得到矩陣W，我們采用隨機采樣的方法，對計算區(qū)域進行劃分，如圖2所示。然后通過隨機采樣的局部空間可以重新組成一個新的矩陣Gs，它是G的一個子矩陣。由于局部區(qū)域內(nèi)的旅行時、振幅、出射角都是低頻光滑函數(shù)，因此局部具有信息冗余，導(dǎo)致矩陣Gs是一個低秩的矩陣，在隨機采樣點數(shù)達到一定數(shù)量后將大于矩陣的秩，此時認為Gs的信息已經(jīng)能夠反映這個局部區(qū)域內(nèi)的特征，利用矩陣Gs來構(gòu)造我們需要的矩陣W是可行的。下面詳細介紹矩陣W的獲取方法。

圖2 二次震源分布與局部區(qū)域的隨機采樣

假設(shè)Gs是一個隨機采樣后ns×ms的矩陣，包含ns行和ms列數(shù)據(jù)，首先對Gs進行可截斷的正交三角（QR）分解：

其中，Q1為正交矩陣，R1是上三角矩陣，對角線元素分別對應(yīng)矩陣Gs的特征值，并且按照大小降序排列，這樣可以剪裁矩陣R1對角元素小于門檻值e的行和列，R1成為｜J｜×｜J｜的矩陣Rc，Q1變?yōu)閚s×｜J｜的矩陣Qc，令Gc＝QcRc，則Gs被裁減為ns×｜J｜的矩陣，可見通過QR分解可以選擇能夠表述矩陣Gs的一個子矩陣Gc＝G（∶，J）。

采用相似的方法處理Gs的共軛矩陣G＊s，同樣利用QR分解處理G＊s：

選取大于門檻值e的行和列得到一個子矩陣G＊r：

那么有：

式中：Gr和Qr為｜I｜×ms的矩陣；Rr為｜I｜×｜I｜的方陣。裁減后的子矩陣可以幫助我們找到一個Gs的近似矩陣，類似于公式（12）：

那么此時W的矩陣維度已經(jīng)變?yōu)椋麶｜×｜I｜。不難得出利用廣義逆矩陣的奇異值分解方法可以獲得W 的表達形式：

上述方法我們稱為低秩特征分解方法，要求Gs矩陣具有低秩的特征，才能通過裁減矩陣降低計算量。回到公式（11），光滑速度模型背景條件下，對于單一頻率ω而言，G（ri；sj）是一個大型的低秩稠密矩陣，可以通過上述低秩特征分解的方法實現(xiàn)矩陣向量積的運算：

實際計算時首先計算公式（20）最右側(cè)的3項矩陣向量乘積，即：

其中，m為二次震源的個數(shù)，其計算量級為｜J｜×｜I｜×m。然后計算：

其中，n為局部區(qū)域內(nèi)的空間離散網(wǎng)格點個數(shù)，其計算量為｜J｜×n。

這樣總體計算量為：

若采用公式（11）直接進行矩陣向量積的計算，則計算量為m×n。當(dāng)矩陣G為低秩矩陣時有｜J｜《m，｜I｜《n，因此可以大幅度降低整體計算量。此外，公式（21）和（22）的矩陣乘法和加法適合于GPU的多核并行計算模式，從而能夠進一步提高計算效率。

3 頻率域多路徑疊前偏移成像

每一種地震波模擬方法都可以開發(fā)研究出一種相應(yīng)的地震偏移成像算法，基本原理就是采用CLEARBOUT提出的波場成像條件［18］：

式中：I（ω，r）表示頻率域成像結(jié)果；PS（ω，r）和PR（ω，r）分別表示震源波場與接收點波場?？紤]炮點震源位置在r0，成像空間位置用r表示的單炮偏移成像問題，將成像空間分成若干水平層狀區(qū)域，如圖3所示。其中，L表示震源與二次震源的層位置。

圖3 成像空間分解方法

依據(jù)前文描述的基于低秩特征分解的二次震源波場構(gòu)建方法，并在第1層內(nèi)采用鏡像邊界條件方法實現(xiàn)震源波場表達形式［19－20］。

當(dāng)在地表炮點位置處，有層位L＝0，第1層內(nèi)的震源波場可以表示為：

在成像空間中二次震源以下的區(qū)域?qū)游籐＞0時：

式中：r′是二次震源位置；nr′表示二次震源個數(shù)；h表示成像網(wǎng)格間距；PS（r0）表示在震源r0處的子波函數(shù)。公式（25）和公式（26）形成逐層遞推的波場延拓公式，第1層和其它層的波場構(gòu)建方法不同。

下面利用同樣的方式構(gòu)建接收點波場表達，不同之處只是把每個接收點當(dāng)作一個二次震源進行波場構(gòu)建。

當(dāng)在地表接收點位置L＝0時存在nR個檢波器接收信號，那么此時的波場有如下表達：

在成像空間中二次震源以下的區(qū)域L＞0時：

將（25）式、（26）式、（27）式和（28）式代入成像條件公式（24）可以得到大步長分層延拓的成像公式：

式中：nS表示單炮數(shù)量：NL表示區(qū)域邊界數(shù)量：nω表示離散采樣的頻率數(shù)量。這樣我們可以通過逐層延拓波場，逐層成像的方式完成整個區(qū)域的偏移成像算法，偽代碼流程如圖4所示。

圖4 多路徑偏移算法偽代碼流程

4 數(shù)模實驗

為了驗證多路徑波場模擬方法，我們設(shè)計了按照正弦函數(shù)變化的速度模型（圖5），其變化規(guī)律描述為：

模型的長度為2．5km，深度為2．5km，成像網(wǎng)格間距為10m，震源位置為（1．0km，0．1km）分別對應(yīng)（x，z）坐標(biāo)。全局射線追蹤可以發(fā)現(xiàn)，速度的非均勻性導(dǎo)致射線存在交叉現(xiàn)象，即“焦散”現(xiàn)象。本文選用20Hz頻率參數(shù)進行波場構(gòu)建，對比分析單路徑和多路徑條件下的波場差異，并統(tǒng)計計算時間。依據(jù)本文的多路徑模擬方法，需要進行區(qū)域劃分，我們將兩個由二次震源組成的平面放在如圖5所示的層位位置，保證在每個區(qū)域內(nèi)不發(fā)生焦散現(xiàn)象。圖5中的紅色線條表示從震源發(fā)出的射線路徑，在通過焦散區(qū)的時候射線發(fā)生了交叉現(xiàn)象。圖5中的藍色線條表示二次震源的位置。

計算結(jié)果如圖6所示，在焦散以外的區(qū)域，單路徑模擬可以得到正確結(jié)果。但在焦散區(qū)，多路徑波場十分復(fù)雜，多條射線干涉現(xiàn)象明顯，常規(guī)的單路徑方法已不能獲得正確波場。多路徑波場模擬的正確性在后文的偏移成像中可以得到驗證。

圖5 數(shù)值速度模型、射線路徑與二次震源的區(qū)域分解方案

圖6 圖5所示模型的單路徑波場（a）與多路徑波場（b）

計算效率對比結(jié)果如表1所示。

表1 計算效率對比結(jié)果

可見，基于低秩分解的多路徑模擬算法采用CPU計算時效率比單路徑方法有所降低，但采用GPU計算后，計算效率已經(jīng)與單路徑方法相當(dāng)，如果直接采用矩陣向量積進行計算，即使采用GPU加速仍然不能獲得理想的效果。因此，本文的低秩分解優(yōu)化算法成功提高了計算效率，達到能夠滿足實際應(yīng)用需求的目的。

為了測試本文成像方法的優(yōu)缺點，我們設(shè)計了2個層次的實驗，分別是脈沖響應(yīng)實驗和簡單模型實驗。

首先，利用脈沖響應(yīng)測試算子的優(yōu)劣，分別與同類偏移算法進行對比。我們?nèi)匀徊捎脠D5所示的速度模型做脈沖響應(yīng)計算，在震源位置設(shè)置一個接收點數(shù)據(jù)，并給予雙程旅行時延遲時間，然后作為地震接收數(shù)據(jù)輸入，完成成像過程。并將之與單程波算法結(jié)果進行比較，效果如圖7所示。

從圖7a，圖7c，圖7d，圖7e和圖7f可見，本文方法在焦散區(qū)的成像有“蝴蝶狀”的脈沖響應(yīng)，與單程波方法相似，而單路徑成像方法并沒有這種現(xiàn)象。由于單程波方程求解能夠自然解決波場的多路徑問題，因此我們認為“蝴蝶狀”脈沖響應(yīng)是正確結(jié)果，同時也驗證了圖6b中多路徑波場模擬的正確性。

第二步，我們設(shè)計的速度模型包含一個水平反射層和一個速度異常體，如圖8a所示。模型寬度為1．5km，深度為2．5km，在深度1．7km處有一個速度反射界面，模型中心有一個速度漸變異常區(qū)，導(dǎo)致射線彎曲和交叉，出現(xiàn)焦散現(xiàn)象。我們利用有限差分時間域正演模擬單炮記錄，檢波器與震源位于同一深度0．02km，震源的水平坐標(biāo)為0．75km，檢波器分布在整個橫向網(wǎng)格點。其中一個單炮記錄如圖8b所示。

圖7 脈沖響應(yīng)測試

圖8 速度模型（a）與單炮記錄（b）

可見，單炮記錄中有明顯的分層信息，并不是一個簡單的反射同相軸，這正是焦散引起的波場畸變，用常規(guī)的單路徑模擬的成像算法無法獲得真實反射界面的形態(tài)特征。我們分別用單路徑和多路徑模擬方法及單路徑和多路徑偏移成像方法進行驗證，各自的波場與成像結(jié)果如圖9和圖10所示。

圖9 圖8所示模型的單路徑波場（a）與多路徑波場（b）

從圖9a可以看出，焦散區(qū)的單路徑波場已經(jīng)發(fā)生畸變，不能正確反映實際波場現(xiàn)象，而本文方法（圖9b）則能夠正確反映焦散區(qū)的波場。從圖10a所示的成像結(jié)果看，單路徑偏移成像方法在異常體下方的成像結(jié)果與實際模型不符。多路徑偏移成像結(jié)果與實際模型位置吻合更好（圖10b），因此本文的多路徑偏移方法對于非均勻介質(zhì)的適應(yīng)性更強。

圖10 單路徑偏移（a）與多路徑偏移（b）成像結(jié)果

5 結(jié)論

本文介紹了一種復(fù)雜介質(zhì)中的基于二次震源多路徑波場模擬方法，該方法采用隨機采樣和低秩分解的Green函數(shù)并利用GPU實現(xiàn)焦散區(qū)波場的快速正確計算，提高了頻率域射線類深度偏移方法的成像精度。數(shù)模實驗結(jié)果顯示：常規(guī)單路徑的Kirchhoff疊前深度偏移成像位置發(fā)生畸變，本文方法能夠獲得正確的成像結(jié)果。該方法在解決復(fù)雜地質(zhì)體的成像問題（如我國西部、南方的復(fù)雜深度成像問題）中有廣泛的應(yīng)用前景。目前，本文方法在推向大規(guī)模實際生產(chǎn)應(yīng)用的過程中還面臨著如地震波的走時、振幅與傳播角度的計算精度問題，以及如何最優(yōu)化設(shè)計二次震源的位置與數(shù)量等理論問題，需要進一步研究。

致謝：感謝美國密歇根州立大學(xué)錢建良教授在本文研究過程中給予的理論支持和寶貴建議，感謝中國科學(xué)院地質(zhì)與地球物理研究所劉洪研究員提供的大力支持與幫助。

［1］ 蔡杰雄，方伍寶，楊勤勇．高斯束深度偏移的實現(xiàn)與應(yīng)用研究［J］．石油物探，2012，51（5）：469－475

CAI J X，F(xiàn)ANG W B，YANG Q Y．Realization and application of Gaussian beam depth migration［J］．Geophysical Prospecting for Petroleum，2012，51（5）：469－475

［2］ 崔興福，劉衛(wèi)東，劉桂寶，等．平面波偏移、分角度成像與AVA道集生成［J］．石油物探，2007，46（6）：615－620

CUI X F，LIU W D，LIU G B，et al．Plan wave migra－tion，angle Imaging and AVA gather production［J］．Geophysical Prospecting for Petroleum，2007，46（6）：615－620

［3］ 袁茂林，黃建平，李振春，等．局部角度域高斯束偏移參數(shù)優(yōu)選研究［J］．石油物探，2015，54（5）：602－612

YUAN M L，HUANG J P，LI Z C，et al．Parameters optimization in local angle－domain Gaussian beam migration［J］．Geophysical Prospecting for Petroleum，2015，54（5）：602－612

［4］ HILLN R．Gaussian beam migration［J］．Geophysics，1990，55（11）：1416－1428

［5］ HALE D．Migration by the Kirchhoff slant stack and Gaussian beam methods［R］．Center for Wave Phenomena，Colorado School of Mines，1992，CWP－126

［6］ GAO F，ZHANG P，WANG B，et al．Fast beam migration a step toward interactive imaging［J］．Expanded Abstracts of 76thAnnual Internat SEG Mtg，2006，2470－2473

［7］ SUN H，SCHUSTER G T．2－D wavepath migration［J］．Geophysics，2001，66（5）：1528－1537

［8］ 李輝，王華忠，馮波，等．特征高斯波包疊前深度偏移方法［J］．地球物理學(xué)報，2014，57（7）：2258－2268

LI H，WANG H Z，F(xiàn)ENG B，et al．Efficient pre－stack depth migration method using characteristic Gaussian packets［J］．Chinese Journal of Geophysics，2014，57（7）：2258－2268

［9］ 王雄文，王華忠．基于壓縮感知的高分辨率平面波分解方法研究［J］．地球物理學(xué)報，2014，57（9）：2946－2960

WANG X W，WANG H Z．A research of high－resolution plane－wave decomposition based on compressed sensing［J］．Chinese Journal of Geophysics，2014，57（9）：2946－2960［10］ CHENG H，CRUTCHFIELD W Y，GIMBUTAS Z，et al．A wideband fast multipole method for the Hemlholtz equation in three dimensions［J］．Journal of Computational Physics，2006，216（1）：300－325

［11］ CHEW W C，LU C C．The use of Huygens’equivalence principle for solving the volume integral equation of scattering［J］．IEEE Transactions on Antennas ＆Propagation，1993，41（7）：897－904

［12］ ROKHLIN V．Diagonal forms of translation operators for the helmholtz equation in three dimensions［J］．Applied ＆Computational Harmonic Analysis，1993，1（1）：82－93

［13］ ENGQUISTB，YING L．Fast directional multilevel algorithms for oscillatory kernels［J］．SIAM Journal on Scientific Computing，2007，29（4）：1710－1737

［14］ GREENGARD L，ROKHLIN V．A fast algorithm for particle simulations［J］．Journal of Computational Physics，1987，73（2）：325－348

［15］ MESSNER M，SCHANZ M，DARVE E．Fast direc－tional multilevel summation for oscillatory kernels based on Chebyshev interpolation［J］．Journal of Computational Physics，2012，231（4）：1175－1196

［16］ MICHIELSSEN E，BOAG A．A multilevel matrix decomposition algorithm for analyzing scattering from large structures［J］．Antennas ＆Propagation IEEE Transactions on，1996，44（8）：1086－1093

［17］ ROKHLIN V．Rapid solution of integral equations of scattering theory in two dimensions［J］．Journal of Computational Physics，1990，86（2）：414－439

［18］ CLAERBOUT J F．Toward a unified theory of reflector mapping［J］．Geophysics，1971，36（3）：467－481

［19］ BERKHOUT A J，WAPENAAR C P A．One－way versions of the Kirchhoff integral［J］．Geophysics，1989，54（4）：460－467

［20］ BLEISTEIN N．On the imaging of reflectors in the earth［J］．Geophysics，1987，52（7）：931－942

（編輯：顧石慶）

Multipath seismic simulation and imaging in frequency domain based on random sampling method

LI Bo1，2
（1．Sinopec Geophysical Research Institute，Nanjing211103，China；2．Key Laboratory of Petroleum Resources Research，Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing100029，China）

Conventional Kirchhoff prestack depth migration cannot handle the wave fields with caustic phenomenon；consequently，the subsurface reflector imaging could be distorted．In the complex media，Kirchhoff migration based on ray tracing must be inevitably affected by the caustic problem．In this paper，a multipath wave field simulation method based on the secondary source Huygens theory in complex media is proposed．Firstly，we ensure that there is no ray－cross in the local scope for domain decomposition．Then we use the secondary source Kirchhoff－Huygens integral method to gradually complete wave field extrapolation for the entire imaging area．In this way，we can correctly handle the caustic phenomena and get the correct imaging reflector．The wave field simulation method and domain decomposition strategy of secondary sources definition we described here have included large－scale matrix calculation that is the main bottleneck of efficiency．We propose a wavefield simulation algorithm based on random sampling and low rank decomposition method for large－scale matrix multiplication．On the other hand，we use GPU as the computing devices to realize the multipath seismic simulation and imaging in frequency domain，which greatly improve the computational efficiency．The numerical simulation test results show caustic could be removed and the correct reflector imaging is realized with this method．

seismic imaging，wave field extrapolation，caustic，Green function，low－rank decomposition，random sampling

P631

A

1000－1441（2017）03－0382－08

10．3969／j．issn．1000－1441．2017．03．008

李博．基于隨機采樣的頻率域多路徑波場模擬與偏移成像［J］．石油物探，2017，56（3）：382－389

LI Bo．Multipath seismic simulation and imaging in frequency domain based on random sampling method［J］．Geophysical Prospecting for Petroleum，2017，56（3）：382－389

2016－10－09；改回日期：2017－03－16。

李博（1981—），男，博士，高級工程師，主要從事復(fù)雜地震波傳播與成像算法、地球物理方法高性能計算方法研究。

國家重點研發(fā)計劃項目（2016YFC0601101）資助。

This research is financially supported by the Nation Key R ＆D Program of China（Grant No．2016YFC0601101）．