李 博

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

基于隨機(jī)采樣的頻率域多路徑波場(chǎng)模擬與偏移成像

李 博1，2

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

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

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

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

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

1 頻率域多路徑波場(chǎng)模擬基本原理

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

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

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

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

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

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

將公式（4）代入公式（1）得到波場(chǎng)的積分近似表達(dá)：

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

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

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

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

代入公式（5）也可獲得波場(chǎng)的積分表達(dá)式：

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

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

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

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

離散化的波場(chǎng)積分公式（8）可以表達(dá)為如下形式：

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

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

波場(chǎng)計(jì)算的公式（9）實(shí)際上可以簡(jiǎn)化為大型矩陣的乘法和加法，為此，首先做如下變量定義：

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

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

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

公式（10）的第1項(xiàng)可表達(dá)為：

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

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

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

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

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

假設(shè)Gs是一個(gè)隨機(jī)采樣后ns×ms的矩陣，包含ns行和ms列數(shù)據(jù)，首先對(duì)Gs進(jìn)行可截?cái)嗟恼蝗牵≦R）分解：

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

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

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

那么有：

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

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

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

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

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

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

這樣總體計(jì)算量為：

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

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

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

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

圖3 成像空間分解方法

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

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

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

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

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

當(dāng)在地表接收點(diǎn)位置L＝0時(shí)存在nR個(gè)檢波器接收信號(hào)，那么此時(shí)的波場(chǎng)有如下表達(dá)：

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

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

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

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

4 數(shù)模實(shí)驗(yàn)

為了驗(yàn)證多路徑波場(chǎng)模擬方法，我們?cè)O(shè)計(jì)了按照正弦函數(shù)變化的速度模型（圖5），其變化規(guī)律描述為：

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

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

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

圖6 圖5所示模型的單路徑波場(chǎng)（a）與多路徑波場(chǎng)（b）

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

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

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

為了測(cè)試本文成像方法的優(yōu)缺點(diǎn)，我們?cè)O(shè)計(jì)了2個(gè)層次的實(shí)驗(yàn)，分別是脈沖響應(yīng)實(shí)驗(yàn)和簡(jiǎn)單模型實(shí)驗(yàn)。

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

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

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

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

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

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

圖9 圖8所示模型的單路徑波場(chǎng)（a）與多路徑波場(chǎng)（b）

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

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

5 結(jié)論

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

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

［1］ 蔡杰雄，方伍寶，楊勤勇．高斯束深度偏移的實(shí)現(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é)報(bào)，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é)報(bào)，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ī)采樣的頻率域多路徑波場(chǎng)模擬與偏移成像［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—），男，博士，高級(jí)工程師，主要從事復(fù)雜地震波傳播與成像算法、地球物理方法高性能計(jì)算方法研究。

國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（2016YFC0601101）資助。

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