岳曉鵬，付亞運(yùn)，陳康靖

（許昌學(xué)院數(shù)理學(xué)院，河南 許昌 461000）

在地球內(nèi)部廣泛存在各向異性，如西非洲海岸、墨西哥灣及中國(guó)南海部分區(qū)域都存在各向異性地層，橫向各向同性(TI)最為多見的，這種介質(zhì)具有柱對(duì)稱軸，對(duì)稱軸在空間垂直的稱為VTI 介質(zhì)；對(duì)稱軸在空間水平的稱為和HTI 介質(zhì)。正演是偏移和反演的基礎(chǔ)，有限差分正演可以較準(zhǔn)確地描述波傳播規(guī)律，是波場(chǎng)正演時(shí)應(yīng)用廣泛的方法。頻率域波動(dòng)方程正演早期由Lysmer 等提出并利用該方法研究了波在各種介質(zhì)中的特征；Jo 等人提出了減小數(shù)值頻散的最優(yōu)化9點(diǎn)格式[1]；Min 等提出了減小了空間采樣點(diǎn)數(shù)的25 點(diǎn)差分格式；吳國(guó)忱教授用25 點(diǎn)差分格式正演了VTI和TTI 介質(zhì)[2]，還有國(guó)內(nèi)的諸多學(xué)者對(duì)于VTI 介質(zhì)的數(shù)值模擬進(jìn)行了研究[3-5]。

頻率域正演優(yōu)勢(shì)在于沒有時(shí)間累積誤差[6]，每次用獨(dú)立的頻率正演，更適用于并行計(jì)算[7]。在并行計(jì)算研究方面，出現(xiàn)過很多的并行框架，GPU 并行、CPU 并行、集群并行等。很多框架設(shè)計(jì)結(jié)構(gòu)繁雜，對(duì)開發(fā)人員要求很高[8，9]。

本研究在詳細(xì)分析各種常見并行模式的基礎(chǔ)上，設(shè)計(jì)了VTI 介質(zhì)多核并行計(jì)算方法。與CPU 串行算法相比，具有較好的加速比。

1 VTI 介質(zhì)頻率域25 點(diǎn)格式

對(duì)于二維VTI 介質(zhì)質(zhì)頻率域波方程如下

在式(1)中所有項(xiàng)，利用25 點(diǎn)格式離散時(shí)所需的點(diǎn)分布見圖1。

ρω2U 系數(shù)見圖1(a)，表達(dá)式為

圖1 25 點(diǎn)加權(quán)離散示意圖

同理，可以展開其它幾項(xiàng)，將所有算子帶入公式(1)并進(jìn)行整理后，可以得到矩陣方程

如果網(wǎng)格為N2×Nx，N=Nz×Nx，A 為2N 階矩陣，Y是未知2N×1 矩陣，F(xiàn) 是2N×1 矩陣。

2 數(shù)值算例

地質(zhì)模型為:縱波為4 000 m/s；橫波2 000 m/s；介 質(zhì) 2 500 kg/m3；NX=160，NZ=160； 空 間 步 長(zhǎng)Δx=Δz=5 m，ε=0.5，δ=0.1。邊界是PML 邊界。記錄時(shí)間0.3 s，采樣間隔0.001 s，震源是Richer 子波，主頻25 Hz，檢波器在地表按網(wǎng)格排列，震源坐標(biāo)為(80，80)，檢波點(diǎn)坐標(biāo)(100，100)，圖2 和圖3 為數(shù)值模擬的波場(chǎng)快照和地震記錄。

圖2 VTI 介質(zhì)的波場(chǎng)快照（125 ms 時(shí)刻）

圖3 VTI 介質(zhì)的地震記錄

3 多核并行正演對(duì)比分析

這里采用的多核并行方式是使用Matlab 的多核并行方法。正演模擬時(shí)使用的平臺(tái)為雙路Xeon E5-2650 v4，32 GB 內(nèi)存的服務(wù)器，共有24 個(gè)計(jì)算核心。程序運(yùn)行的主要耗時(shí)在矩陣方程(6)的求解上，例如當(dāng)網(wǎng)格數(shù)為160*160 時(shí)，矩陣A就是一個(gè)80 000 階（200*200*2）的大型稀疏矩陣，求解這樣系數(shù)的矩陣方程相當(dāng)耗費(fèi)時(shí)間。對(duì)于VTI 介質(zhì)頻率域方程，可以采用并行的方式在每個(gè)頻點(diǎn)獨(dú)立求解。

假設(shè)計(jì)算網(wǎng)格數(shù)分別為100*100、130*130、160*160，190*190、220*220，四周加上PML 邊界后，實(shí)際計(jì)算網(wǎng)格數(shù)為140*140、170*170、200*200，230*230、260*260，對(duì)比順序執(zhí)行和并行求解方程(6)平均耗時(shí)及程序運(yùn)行時(shí)間，見表1。

表1 兩種方法耗時(shí)對(duì)比

為了更直觀的說明問題，將以上數(shù)據(jù)作圖表示，見圖4。

圖4 兩種方法耗時(shí)對(duì)比

4 結(jié)論

目前，CPU 正在朝著多核心方向快速發(fā)展，最新的i9-13900k 已經(jīng)是24 核32 線程，然而，現(xiàn)有的波動(dòng)方程正演及偏移成像代碼大多是采用MPI 技術(shù)實(shí)現(xiàn)并行，這種進(jìn)程級(jí)粒度的計(jì)算方式在集群之類的分布式系統(tǒng)上效果很好，在單個(gè)計(jì)算節(jié)點(diǎn)又要受到內(nèi)存等限制，只能利用少數(shù)的計(jì)算核心，不能發(fā)揮多核處理器的最大效能。本研究將占據(jù)波動(dòng)方程頻率正演模擬主要計(jì)算時(shí)間的大型矩陣方程求解過程分配給多個(gè)計(jì)算核心，利用多核心并行計(jì)算方法，多個(gè)頻點(diǎn)的求解過程可以同步計(jì)算完畢，即在進(jìn)程中實(shí)現(xiàn)多核心并行，有效提高了多核處理器的利用效率和代碼執(zhí)行速度。