偽深度域交錯網格逆時偏移成像方法及并行優化

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

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

0 引言

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

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

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

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

1 方法原理

1.1 偽深度域變換

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

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

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

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

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

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

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

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

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

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

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

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

式中：U、W分別為偽深度域中ξ、η方向對應速度分量；P為偽深度域的壓力波場；vτ為偽深度域速度場。

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

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

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

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

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

圖2 偽深度域交錯網格剖分示意圖

表1 交錯網格法差分系數表［27］

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

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

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

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

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

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

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

1.3 C＋＋AMP并行加速

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

表2 GPU端不同并行架構的優勢及主要限制對比

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

圖3 并行優化算法示意圖

圖4 不同總網格數時C＋＋AMP GPU并行架構和與CPU常規正演計算耗時對比

2 模型試算

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

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

2.1 鹽丘模型

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

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

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

圖5 逆掩構造模型不同速度場反變換結果及其相對誤差

圖6 鹽丘模型

2.2 黔北鳳崗斷塊模型

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

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

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

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

圖8 鹽丘模型逆時偏移結果

圖9 鳳崗斷塊模型速度場變換結果

2.3 復雜逆掩斷層模型

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

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

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

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

圖12 逆掩構造模型速度場變換結果

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

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

圖14 逆掩構造模型逆時偏移結果

表3 逆時偏移算法網格數目及耗時對比統計

3 結論

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

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

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

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

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