偽深度域交錯網(wǎng)格逆時偏移成像方法及并行優(yōu)化

金宗瑋 黃金強 王甘露* 夏 鵬 牟雨亮

（①貴州大學資源與環(huán)境工程學院，貴州貴陽550025；②貴州省石油學會，貴州貴陽550000）

0 引言

逆時偏移法自二十世紀七、八十年代提出以來，已經(jīng)持續(xù)發(fā)展了近四十余年［1-2］，作為當今對復雜構(gòu)造成像中最為常用手段之一，因其較好的成像效果及實現(xiàn)容易等特點受到了青睞［3-5］。逆時偏移可根據(jù)目標的不同及模型復雜程度分為聲介質(zhì)偏移、黏彈性介質(zhì)偏移以及各向異性介質(zhì)偏移等［6］，黏彈性介質(zhì)及各向異性介質(zhì)的偏移目前尚處于探索實驗階段，基于聲波的逆時偏移技術(shù)是當前的研究熱點，但聲波逆時偏移技術(shù)的推廣和應用仍受較高的存儲占用和較低的計算效率制約［7-8］。

與Kirchhoff偏移相比，基于雙程波動方程的逆時偏移技術(shù)無成像傾角限制，且能夠同時對回轉(zhuǎn)波、多次波及棱柱波精確成像［9-10］，可以在橫向速度變化劇烈時取得不錯的成像結(jié)果［11］。在正演及逆時偏移中，為保證薄層、縫洞及高傾角層位的有效信息精確獲取，往往需取較小的模型網(wǎng)格，由此便導致了網(wǎng)格劃分數(shù)目和計算量的增加。針對該問題，學者們提出了不同的網(wǎng)格劃分優(yōu)化方法：不規(guī)則網(wǎng)格法通過引入非規(guī)則網(wǎng)格差分算子，將規(guī)則四邊形網(wǎng)格映射至任意形狀四邊形網(wǎng)格，對崎嶇地表有較強適應能力［12-14］；變網(wǎng)格法，通過對模型速度不同區(qū)域采用變網(wǎng)格步長的方式，在保證目標區(qū)域信息準確性的情況下，實現(xiàn)了網(wǎng)格數(shù)目的整體減少，能夠有效降低存儲占用和計算量［15-16］；自適應網(wǎng)格法，通過閾值控制的方式，自動將原始模型均勻網(wǎng)格細化改造，使計算中所有網(wǎng)格誤差值均處在設定閾值之內(nèi)，有效地提高了計算效率及地層分辨能力［17-18］。但上述各種網(wǎng)格改進方法都是在笛卡爾坐標系下進行的，并不能徹底解決波場在不同速度層中采樣不均勻的問題。

偽深度域法不同于傳統(tǒng)的算法改進，通過引入曲坐標變換，將原始笛卡爾坐標系下模型速度場變換為縱向等間隔的單程旅行時速度場（又稱偽深度域速度場或垂直時間域速度場），既避免了波場在高速度層中的過采樣現(xiàn)象，又減少了模型縱向采樣點數(shù)，從而達到降低計算量及提高計算效率的目的。偽深度域最早由Alkhailfah等［19］提出，Ma等［20］推導并分析了偽深度域下速度—應力方程，解決了縱向過采樣問題。李慶洋等［21］為解決橫向過采樣問題，引入自適應差分算子計算橫向空間微分，綜合考慮了橫、縱向空間采樣問題，進一步提升了計算效率。但是，使用CPU（MPI并行）進行算法實現(xiàn)仍需消耗大量時間，使該方法仍然難以投入實際應用。因此，運用GPU并行就顯得格外重要。GPU并行多采用CUDA和OpenCL兩種架構(gòu)。二者均有較強的數(shù)據(jù)并行處理能力，但都有編程環(huán)境配置繁復、編程語言不夠精煉、初學者門檻較高等缺點。C＋＋AMP是基于Direct X技術(shù)實現(xiàn)的并行計算庫，集成于微軟VisualStudio軟件開發(fā)工具中，用戶無需進行復雜的環(huán)境配置，只需簡單的引用相應頭文件即可實現(xiàn)對應并行函數(shù)的調(diào)用，相比而言，上手難度低且加速效果并不差［22-23］。

本文將曲坐標變換引入到聲波正演及逆時偏移，同時結(jié)合C＋＋AMP并行編程架構(gòu)，實現(xiàn)了偽深度域下的聲波方程正演及逆時偏移算法，通過標準鹽丘模型、黔北頁巖勘探區(qū)鳳崗斷塊模型及復雜逆掩構(gòu)造模型的試算，證實了本文方法在存儲空間、I/O及計算成本上都具有顯著的優(yōu)勢。

1 方法原理

1.1 偽深度域變換

二維時空域一階速度—應力聲波方程可表示為

式中：p為壓力波場；u、w分別為橫、縱向振動速度分量場；ρ為介質(zhì)密度；v為聲波速度；S為震源函數(shù)。

偽深度域法需先利用積分公式

將笛卡爾坐標系下平滑速度場模型進行縱向積分［19］，以此得到單程旅行時間場τ（x，t），然后再將原速度模型用樣條插值法插值成縱向等間隔（Δτ）的偽深度域速度模型。式中vsm為深度域平滑速度場。相應地，可再次利用樣條插值法將偽深度域速度場反變換回深度域速度場（圖1）。

圖1 偽深度域與深度域正反變換示意圖

在得到偽深度域速度模型后，還需將深度域聲波方程改寫為偽深度域聲波方程。

將深度域波動方程轉(zhuǎn)化到偽深度域時，對應的坐標轉(zhuǎn)換為

為方便起見，將ηx記為α，則有

假設笛卡爾坐標系x、z方向單位矢量為i、k；偽深度域坐標系中ξ、η方向的單位矢量為e1、e2，則其對應關系可表示為［24］

由散度定理可知，任意矢量場M在偽深度域坐標系下的散度和梯度的表達式分別為

式中M1和M2分別為M的ξ和η方向分量。

將式（7）和式（8）代入式（1），在考慮密度為常值的情況下可得到偽深度域二維時空域一階速度—應力方程

式中：U、W分別為偽深度域中ξ、η方向?qū)俣确至浚籔為偽深度域的壓力波場；vτ為偽深度域速度場。

在不考慮參數(shù)α（α＝0）的影響時，式（10）可簡化為

相比于聲波方程（式（1）），簡化后的偽深度域聲波方程（式（11））只在縱向速度場更新中多乘了一個系數(shù)。這樣只需更新更少的網(wǎng)格即可實現(xiàn)相同的正、反演效果，而計算量保持不變，從而減少存儲空間占用和提高計算效率。

1.2 方程離散及吸收邊界條件

本文偽深度域聲波方程離散和網(wǎng)格剖分與常規(guī)深度域聲波方程所運用的有限差分交錯網(wǎng)格法［25-26］相類似。但值得注意的是，在考慮式（10）中α參數(shù)不為零的情況時，為保證計算的合理性，在對偽深度域中U、W更新時需額外利用到壓力場網(wǎng)格中間數(shù)據(jù)P′。因此，以四周壓力場P的均值作為中間壓力場P′，具體網(wǎng)格剖分方式如圖2所示。

對于時間和空間上的差分近似，采用與傳統(tǒng)深度域交錯網(wǎng)格法相同的離散規(guī)則，即壓力和速度場的計算分別對應于t和t＋Δt/2時刻進行。而空間差分近似則以高階半網(wǎng)格方式實現(xiàn)，相應的空間差分系數(shù)如表1所示。

圖2 偽深度域交錯網(wǎng)格剖分示意圖

表1 交錯網(wǎng)格法差分系數(shù)表［27］

將偽深度域二維一階速度—應力方程離散后可得

式中：Δx、Δτ分別為偽深度域縱橫向采樣間隔；Cl為有限差分系數(shù)；N為差分階次的一半。

為解決人工引起的邊界反射，本文采用分裂式PML吸收邊界條件［28-29］，將U、W和P分別沿ξ、η方向分解為兩個部分

同時引入邊界衰減因子σ1、σ2，對應帶衰減因子的分裂式偽深度域聲波方程可寫為

其中衰減因子σ1和σ2分別為

式中：vmax為速度場中最大速度值；δ為PML吸收層厚度；R為理想邊界反射系數(shù)（通常取值10－6）；ξ′和η′分別表示橫、縱向計算點到正常計算區(qū)域之間的網(wǎng)格距離。

最終，得到帶PML吸收邊界條件的偽深度域二維時空域一階速度—應力離散方程

1.3 C＋＋AMP并行加速

現(xiàn)階段并行計算架構(gòu)的選擇有：MPI、OpenCL、CUDA和近幾年出現(xiàn)的C＋＋AMP。其中MPI可以利用多個CPU并行協(xié)同處理同任務，而另外三種則屬于GPU并行計算架構(gòu)。與CPU并行架構(gòu)相比，GPU并行架構(gòu)計算效率更高，是當前的主要發(fā)展方向。三種GPU并行架構(gòu)對比如表2所示。C＋＋AMP并行計算架構(gòu)相對容易實現(xiàn)且無過多限制，同時C＋＋AMP架構(gòu)還提供了一定的顯存及線程優(yōu)化控制方法，其良好的硬件兼容性使代碼能夠在AMD及NVIDIA等不同硬件設備上運行，大大降低了并行成本，結(jié)合其簡單易學及編程環(huán)境搭建簡單等優(yōu)勢，C＋＋AMP并行架構(gòu)成為Windows系統(tǒng)下并行計算的一個不錯選擇。本文在C＋＋AMP架構(gòu)下實現(xiàn)并行加速優(yōu)化，具體實現(xiàn)過程可大致分為CPU端算法與GPU端算法兩部分組成（圖3）：①CPU端負責簡單數(shù)據(jù)的處理及傳遞，如文件讀取、寫入及數(shù)據(jù)交換等；②GPU端負責大量數(shù)據(jù)的并行計算，如波場P、U、W的迭代更新等。

表2 GPU端不同并行架構(gòu)的優(yōu)勢及主要限制對比

圖4對比了相同正演參數(shù)下CPU常規(guī)計算與C＋＋AMP架構(gòu)下并行計算耗時，可以看出，C＋＋AMP架構(gòu)下的并行計算效率約為并行前的10倍，且隨著模型網(wǎng)格總數(shù)的增加而增高。當然，加速效果同樣與所使用的GPU性能有關。同時需要指出的是，該測試僅是在完成算法且未進行編程優(yōu)化的情況下得到的，在完成編程優(yōu)化后加速效果有望進一步提高。

圖3 并行優(yōu)化算法示意圖

圖4 不同總網(wǎng)格數(shù)時C＋＋AMP GPU并行架構(gòu)和與CPU常規(guī)正演計算耗時對比

2 模型試算

為了驗證并行架構(gòu)下偽深度域算法的可靠性，分別采用二維鹽丘模型、黔北鳳崗斷塊模型和復雜逆掩構(gòu)造模型進行正演及逆時偏移成像試算，采用帶照明補償?shù)某上駰l件。逆向波場重構(gòu)采用邊界存儲法［30］，相比直接儲存每一時刻波場，可降低約90%～95%的存儲空間，大大減少了數(shù)據(jù)讀寫量。

需要指出的是，若采用未平滑模型作為初始速度場，則會在變換過程中因初始速度模型縱向速度值變化不均勻?qū)е伦儞Q結(jié)果存在一定插值誤差，且正反變換兩次的樣條插值會導致一定的誤差積累，在縱向速度變換劇烈處尤為嚴重。通過實驗發(fā)現(xiàn)，采用低平滑速度模型作為正演、偏移速度場，高平滑速度模型作為積分速度場，能夠有效減少變換帶來的插值誤差，同時也能減小參數(shù)α的橫向變化，提高算法的穩(wěn)定性。以復雜逆掩構(gòu)造模型為例，對比低平滑初始速度場與非平滑初始速度場反變換結(jié)果（圖5）可看出，采用低平滑初始速度場較非平滑初始速度場反變換后結(jié)果誤差較小、精度較高，能相對準確地將偽深度域成像結(jié)果反變換回深度域。

2.1 鹽丘模型

模型中鹽丘與圍巖速度差較大，能檢驗正演和偏移算法對高速異常體的適應性。該模型網(wǎng)格數(shù)為676×200，橫、縱向間距均為3m。正演采樣間隔為1ms、長度為1.2s，震源選擇主頻為25Hz的雷克子波；采用雙邊接收方式，最小炮檢距為0，最大炮檢距為1014m，道間距為3m。偽深度采樣間隔Δτ＝2ms，由

計算的偽深度網(wǎng)格點數(shù)nτ＝142，可見偽深度變換可減少約29.0%的存儲需求。式中：τmax為最大單程旅行時；INT（·）為向上取整函數(shù)。原始速度模型和偽深度域正、反變換模型如圖6所示，偽深度變換導致高速區(qū)抽稀采樣、高速體下方呈現(xiàn)“界面上凸”現(xiàn)象（圖6b橢圓所示），做反變換后（圖6c）與原始速度場（圖6a）的相對誤差（圖6d）最大約為4.2%。采用相同模擬參數(shù)，圖6a速度場的深度域正演記錄如圖7a所示，偽深度域正演記錄如圖7b所示，二者相對誤差不超過1.5%（圖7c）。

逆時偏移實驗中，共計模擬169炮，起始炮點位于0，炮點距為12m，其他參數(shù)與正演保持一致，深度域和偽深度域成像結(jié)果如圖8所示。對比常規(guī)深度域與偽深度域反變換偏移結(jié)果可以看出偽深度域逆時偏移對高速異常體及其下方層位（紅色箭頭所示）均有較精確的成像效果。

圖5 逆掩構(gòu)造模型不同速度場反變換結(jié)果及其相對誤差

圖6 鹽丘模型

2.2 黔北鳳崗斷塊模型

根據(jù)黔北鳳崗頁巖勘探區(qū)的前期鉆孔及地震解釋成果，設計了鳳崗斷塊模型（圖9a）。該模型主要展示了研究區(qū)斷層發(fā)育特征，其中上部多為高速灰?guī)r地層，下部為低速泥頁巖地層。模型網(wǎng)格數(shù)為737×275，橫、縱向網(wǎng)格間距均為10m。時間采樣間隔置為0.5ms，選擇主頻為25Hz的雷克子波為震源，模擬時長為1.8s；采用雙邊接收方式，最小炮檢距為0，最大炮檢距為3690m，道間距為10m。偽深度域采樣間隔為3.76ms，偽深度網(wǎng)格點數(shù)為214（圖9b），可見偽深度域變換能節(jié)省22.18%的存儲空間。

同樣，由圖9b能很明顯看出偽深度域速度場的變形現(xiàn)象，尤其是淺色低速地層（泥頁巖層）導致下方地層呈現(xiàn)“下拽”現(xiàn)象（圖中紅圈所示），且其程度隨著斷層傾角的變大而增加。深度域和偽深度域模擬記錄及差值如圖10所示，二者的相對誤差約為1.5%。

起始炮點位于0，炮點距為40m，共模擬184炮，其逆時偏移結(jié)果如圖11所示。鳳崗斷塊模型深部低速頁巖層的成像效果同樣很好（圖11紅色箭頭所示），細節(jié)部分得以充分體現(xiàn)。

圖7 鹽丘模型常規(guī)深度域（a）與偽深度域（b）正演記錄及其差值（c）

圖8 鹽丘模型逆時偏移結(jié)果

圖9 鳳崗斷塊模型速度場變換結(jié)果

2.3 復雜逆掩斷層模型

復雜逆掩斷層模型（圖12a）淺、表層速度變化劇烈，褶皺構(gòu)造、高陡構(gòu)造及逆掩斷層發(fā)育［31］，能夠檢測不同模擬、偏移算法的適用性。該模型網(wǎng)格數(shù)為834×500，橫縱向網(wǎng)格間距均為10m。時間采樣間隔置為0.8ms，選擇主頻為30 Hz的雷克子波為震源，模擬時長為4s；采用雙邊接收方式，最小炮檢距為0m，最大炮檢距為4170m，道間距為10m。偽深度域采樣間隔為2.5ms，偽深度網(wǎng)格點數(shù)為426（圖12b），可見偽深度域變換能夠減少14.8%的存儲需求。同時，為減少速度變換過程中可能出現(xiàn)的抽樣誤差，正演及偏移成像均采用低平滑模型作為初始模型，高平滑模型作為變換速度場。

逆掩構(gòu)造模型深度域和偽深度域模擬記錄及差值如圖13所示，二者的相對誤差約為2.2%，可忽略不計。417炮模擬記錄（間隔20m）的逆時偏移結(jié)果如圖14所示，對比兩種方法成像結(jié)果可知，對于復雜逆掩構(gòu)造模型而言，常規(guī)深度域與偽深度域反變換成像結(jié)果均能很好對應。

圖10 鳳崗斷塊模型常規(guī)深度域（a）與偽深度域（b）正演記錄及其差值（c）

圖11 鳳崗斷塊模型逆時偏移成像結(jié)果

圖12 逆掩構(gòu)造模型速度場變換結(jié)果

由以上三個模型試算結(jié)果可以看出，偽深度域法對斷層構(gòu)造、深部薄互層、高速體邊界及復雜高陡逆掩構(gòu)造都有很好的成像效果，反變換后成像層位依然清晰準確，驗證了算法的正確性和適用性。對比兩種偏移方法（表3）可以看出偽深度域節(jié)省約20%～35%的整體耗時，減少了15%～30%存儲占用。

圖13 逆掩構(gòu)造模型常規(guī)深度域（a）與偽深度域（b）正演記錄及其差值（c）

圖14 逆掩構(gòu)造模型逆時偏移結(jié)果

表3 逆時偏移算法網(wǎng)格數(shù)目及耗時對比統(tǒng)計

3 結(jié)論

本文通過引入曲坐標變換，推導了偽深度域聲波方程及其離散形式，實現(xiàn)了偽深度域正、反向波場延拓。

（1）相比傳統(tǒng)深度域正演，偽深度域法可以起到波場“均衡采樣”作用，變換后波場在偽深度域計算中任何位置都呈現(xiàn)出縱向上均勻采樣的狀態(tài)，能夠避免深度域下高速層位波場的過采樣，減少深度方向的樣點數(shù)。在保證計算精度的情況下，降低存儲空間需求及提高計算效率。標準鹽丘模型以及黔北鳳崗斷塊模型的逆時偏移試算表明，與深度域方法相比，偽深度域法能夠節(jié)省15%～30%的數(shù)據(jù)存儲量，效率提高了20%～35%。

（2）將C＋＋AMP并行架構(gòu)與偽深度域法正演、偏移算法相結(jié)合，大幅提高了計算效率，方便了算法的推廣與應用。

（3）采用低平滑速度場作為正演、偏移速度場能夠有效降低正反變換中抽樣插值帶來的誤差，增加成像結(jié)果反變換精度。

在考慮到現(xiàn)階段實際工程中越發(fā)復雜的地下構(gòu)造及更加精細化的成像需求，將該方法進一步推廣到最小二乘法逆時偏移及全波形反演具有重要意義。