周 震，賈曉峰

(中國科學(xué)技術(shù)大學(xué) 地球和空間科學(xué)學(xué)院，合肥 230026)

基于OpenMP加速的無單元逆時(shí)偏移成像

周 震，賈曉峰*

(中國科學(xué)技術(shù)大學(xué) 地球和空間科學(xué)學(xué)院，合肥 230026)

有限元法(FEM)和有限差分法(FDM)被用于求解地震波動(dòng)方程，但該兩種方法，在實(shí)際運(yùn)用中，計(jì)算精度和效率會(huì)出現(xiàn)一些不足。事實(shí)上，無單元法(EFM)已經(jīng)被用于地震波模擬和偏移成像，與FDM和FEM相比，EFM具有獨(dú)特的優(yōu)勢(shì)：不需要事先劃分大量網(wǎng)格，僅需要研究區(qū)域中節(jié)點(diǎn)和邊界的信息，具有局部擬合的特點(diǎn)。但當(dāng)前EFM處理地震波模擬和成像仍存在很多問題，如求取系數(shù)矩陣耗費(fèi)過多計(jì)算機(jī)內(nèi)存和計(jì)算量過大的問題。這里對(duì)系數(shù)矩陣采用按行壓縮(CSR)的格式存儲(chǔ)，并采用分步計(jì)算的策略。在時(shí)間遞推過程中，求解線性方程組，為了提高計(jì)算效率，基于已有的FDM數(shù)據(jù)，利用OpenMP多核加速。通過上述方法，有效地解決了內(nèi)存限制和計(jì)算效率的問題。數(shù)值算例表明，這里提出的方法是精確和高效的。

無單元法； 系數(shù)矩陣； 線性方程組； OpenMP； 逆時(shí)偏移

0 引言

在地震勘探領(lǐng)域，利用彈性波方程(近似為聲波方程)描述物理震源在空間傳播地震波場(chǎng)的過程，這是地震波方程正演和反演方法的物理基礎(chǔ)。人們一直在探討各種波動(dòng)方程的數(shù)值解法，其目的是為了得到滿意的地震波正反演結(jié)果。在這個(gè)過程中，不可避免的產(chǎn)生一些其他問題，比如數(shù)值算法的精度、計(jì)算機(jī)資源占用率、算法穩(wěn)定性以及邊界條件等。當(dāng)前應(yīng)用最廣泛的數(shù)值計(jì)算方法為有限元法(FEM)和有限差分法(FDM)[3，8]，這兩種方法發(fā)展較早且相對(duì)比較成熟。近年來，隨著研究問題的深入和廣泛，有限差分法和有限元法在計(jì)算精度和計(jì)算效率上出現(xiàn)了一些不足。

在巖石和工程力學(xué)領(lǐng)域中，發(fā)展起來的無單元法(EFM)已在地震勘探中獲得有效的應(yīng)用[7]。無單元法對(duì)于處理有限元法(FEM)中獨(dú)立變量一次導(dǎo)數(shù)不連續(xù)、預(yù)處理較復(fù)雜等問題，具有良好的效果。同時(shí)利用滑動(dòng)最小二乘法來擬合波場(chǎng)函數(shù)，使得EFM不僅繼承了FEM的一些特點(diǎn)，更具有某些獨(dú)特的優(yōu)勢(shì):①無需事先劃分單元，讓EFM更具有靈活性，這使得在預(yù)處理時(shí)無需再去配置和組集單元，可有效降低工作成本；②EFM降低了構(gòu)造形函數(shù)的難度，只需要考慮一個(gè)最小二乘的精度即可；③在EFM采用最小二乘擬合的思想，可獲得高次連續(xù)解。

自2006年起，EFM首次被用于地震波模擬和成像，這次嘗試為EFM在彈性波動(dòng)方程中的應(yīng)用奠定了基礎(chǔ)。其中包括：①權(quán)函數(shù)的選擇；②影響域的大??；③形函數(shù)的構(gòu)造；④矩陣求逆的運(yùn)算等。同時(shí)因?yàn)橄∈杈仃嚥磺‘?dāng)?shù)拇鎯?chǔ)問題，該方法在2G計(jì)算機(jī)內(nèi)存的存儲(chǔ)情況下，僅能處理81×81網(wǎng)格節(jié)點(diǎn)的模型，極大限制了EFM在較大速度模型中的應(yīng)用，同時(shí)對(duì)于超大稀疏矩陣求逆，也是一項(xiàng)挑戰(zhàn)。

通過對(duì)稀疏矩陣壓縮處理，并將矩陣的求逆過程轉(zhuǎn)化為求解線性方程組，F(xiàn)an[4]解決了因稀疏矩陣不恰當(dāng)存儲(chǔ)所引起的存儲(chǔ)瓶頸，同時(shí)通過‘PARDISO’ (Parallel Direct Sparse Solver)求解器[5]求解線性方程組，在計(jì)算精度得到保證的前提下，降低了計(jì)算成本。該方法對(duì)于節(jié)點(diǎn)和時(shí)間步長的選取，同樣給出相應(yīng)的解決方案。然而在計(jì)算效率方面，因?yàn)樵摲椒ㄔ诶酶咚狗e分計(jì)算剛度矩陣和質(zhì)量矩陣時(shí)，采用了堆棧這一數(shù)據(jù)結(jié)構(gòu)，該過程重復(fù)開辟內(nèi)存存儲(chǔ)空間，嚴(yán)重阻礙了計(jì)算效率。

作者在計(jì)算質(zhì)量矩陣(M)和剛度矩陣(K)時(shí)，摒棄了借助堆棧的解決方案，改用數(shù)組存儲(chǔ)中間稀疏矩陣，并對(duì)最終系數(shù)矩陣按行壓縮存儲(chǔ)(CSR)。在計(jì)算M和K的過程中，分兩步進(jìn)行：①在每一個(gè)高斯節(jié)點(diǎn)各自單獨(dú)影響域內(nèi)對(duì)模型節(jié)點(diǎn)計(jì)算出M和K，同時(shí)在計(jì)算過程中采用模型節(jié)點(diǎn)內(nèi)層循環(huán)局部搜索代替全局搜索的策略，并利用相對(duì)坐標(biāo)這一概念；②對(duì)不同高斯點(diǎn)影響域下包含的相同的模型節(jié)點(diǎn)上的M和K求和。其目的就是在不借助堆棧的情況下，利用有限的計(jì)算機(jī)內(nèi)存，盡可能計(jì)算較大速度模型，從而為提高計(jì)算精度。

當(dāng)求出M和K并壓縮存儲(chǔ)后，可通過隱式差分格式得到稀疏矩陣為系數(shù)的線性方程組。之后通過‘PARDISO’線性求解器求解。為了提高計(jì)算效率，利用已有的有限差分(FDM)程序，參考計(jì)算機(jī)的內(nèi)核(CPU)個(gè)數(shù)，保存相應(yīng)幾個(gè)時(shí)刻的波場(chǎng)值及其對(duì)時(shí)間的一階偏導(dǎo)。然后把已求值設(shè)為EFM中的起始點(diǎn)值，利用OpenMP-Fortran 編寫并行程序，使得EFM在波場(chǎng)模擬和逆時(shí)偏移成像的過程中多核同時(shí)進(jìn)行，進(jìn)而提高計(jì)算效率。

1 無單元法基本理論介紹

考慮下面Ω研究區(qū)域，Γ為邊界的二維標(biāo)量波動(dòng)方程為式(1)。

(1)

式中：u表示位移場(chǎng)；t表示時(shí)間坐標(biāo)；x、y分別表示為空間上水平和豎直方向坐標(biāo)；D表示介質(zhì)中波速的平方。

假設(shè)u(x)為研究區(qū)域Ω上的場(chǎng)函數(shù)，各個(gè)模型節(jié)點(diǎn)上的值可表示為

u(xi)=uii=1,2,…,n

(2)

為了使用更加簡單的形式求解方程，我們?nèi)∑浼?jí)數(shù)形式的近似解uh(x)

(3)

其中：pT(x)為m維基函數(shù)；m為基函數(shù)所含元素的個(gè)數(shù)。

當(dāng)在二維情況下有：

一次基函數(shù) pT(x)=[1,x,y]

(4)

二次基函數(shù) pT(x)=[1,x,y,x2,xy,y2]

(5)

我們通過滑動(dòng)最小二乘原理來確定a(x)。通過定義下面的泛函方程式，并求取泛函方程式的最小值：

(6)

其中：w(x-xI)為點(diǎn)x影響區(qū)域內(nèi)的權(quán)函數(shù)其作用范圍在有限區(qū)域內(nèi)[4]。

為了求得a(x)，對(duì)式(6)的泛函方程式求極小值即

(7)

可得

a(x)=A-1(x)B(x)U

(8)

其中

(9)

B(x)= [w(x-x1)p(x1),w(x-x2)p(x2),

…,w(x-xNinf)p(xNinf)]

(10)

U=[u1,u2,…,uNinf]T

(11)

將a(x)帶入方程(3)可得

uh(x)=p(x)A-1(x)B(x)U≡φ(x)U

(12)

其中：φ(x)稱為形函數(shù)。

從方程(12)可知所得近似場(chǎng)函數(shù)，僅與基函數(shù)、權(quán)函數(shù)以及鄰近點(diǎn)場(chǎng)值有關(guān)。

對(duì)方程(12)使用變分原理，即求解該方程等同于找出下面泛函的極小值：

(13)

其中a是罰因子。將近似解(12)代入(13) 并令求導(dǎo)為零，則可得：

(14)

式(14)中K、M和F分別表示為剛度矩陣、質(zhì)量矩陣和載荷向量，其中K、M均為超大稀疏矩陣可用下面的表達(dá)式表示：

(15)

由公式(14)，可以得到如下方程組：

(16)

把平均加速法的隱式差分格式[7]應(yīng)用到公式(14)，最終得到關(guān)于時(shí)間的遞推公式：

(17)

當(dāng)我們暫時(shí)不考慮邊界條件時(shí)，載荷向量Fq可以忽略不考慮，通過前一時(shí)刻的波場(chǎng)值及其對(duì)時(shí)間的一階偏導(dǎo)，聯(lián)立方程組即可求得下一時(shí)刻的波場(chǎng)值和波場(chǎng)值對(duì)時(shí)間的一階偏導(dǎo)。

2 無單元法中影響域的選擇

在無單元法中，每個(gè)計(jì)算節(jié)點(diǎn)權(quán)函數(shù)的影響區(qū)域稱為影響域，其選擇會(huì)直接影響計(jì)算精度和計(jì)算工作量。影響域的選擇要恰當(dāng)，既不 能太小也不能太大：過小會(huì)使得滑動(dòng)最小二乘法的解值不唯一；過大既偏離了無單元法中強(qiáng)調(diào)的局部擬合的特點(diǎn)，更會(huì)增加計(jì)算的工作量，加大計(jì)算成本。

考慮二維情況，并采用圓形影響域，則影響半徑可定義為：

(18)

這里取β=4、m=6； d為節(jié)點(diǎn)密度，可定義為：

(19)其中：node為整個(gè)研究區(qū)域內(nèi)節(jié)點(diǎn)總數(shù)；a、b分別為區(qū)域的長和寬；β是一經(jīng)驗(yàn)常數(shù)，其大小可根據(jù)速度模型離散點(diǎn)間距大小而定，本文模型點(diǎn)間距為10 m～12.5 m，故β=4；m是基函數(shù)所包含的項(xiàng)數(shù)，因這里采用二次基函數(shù)即滿足需求，故取m=6。

3 分步求取質(zhì)量矩陣和剛度矩陣

無單元算法(EFM)中，剛度矩陣(K)和質(zhì)量矩陣(M)是通過高斯數(shù)值積分計(jì)算得到的。其中K、M的計(jì)算和存儲(chǔ)是限制EFM發(fā)展的瓶頸之一。為了解決內(nèi)存限制并提高計(jì)算效率，采用分步求取K、M的策略：①在各自單獨(dú)高斯節(jié)點(diǎn)影響域內(nèi)對(duì)模型節(jié)點(diǎn)計(jì)算出K和M；②再對(duì)不同高斯點(diǎn)影響域內(nèi)包含的相同的模型節(jié)點(diǎn)的K和M求和。有關(guān)分步計(jì)算K和M的過程可參考圖1。

圖1 分步計(jì)算K和M的流程圖Fig.1 The flow diagram of computing K,M by multi-step

考慮到速度模型的特性(通常水平長度大于豎直深度)，為了便于后期求和，即讓每個(gè)高斯點(diǎn)影響域內(nèi)的模型點(diǎn)坐標(biāo)值相差較小，我們把模型點(diǎn)按照自上而下，然后從左至右依次排序編號(hào)。此時(shí)在程序上，對(duì)于高斯點(diǎn)循環(huán)和速度模型離散點(diǎn)的兩層循環(huán)要滿足水平(X)方向在外，豎直(Y)方向在內(nèi)的原則。

我們采用一種策略：讓內(nèi)層循環(huán)僅出現(xiàn)在相應(yīng)高斯點(diǎn)影響域附近?？紤]到高斯點(diǎn)的影響域是圓形，故可以做圓的外切正方形(圖2(b))，使內(nèi)層模型節(jié)點(diǎn)的循環(huán)僅出現(xiàn)在正方形區(qū)域內(nèi)，讓局部搜索循環(huán)替代之前的全局搜索循環(huán)。這樣就保證了模型節(jié)點(diǎn)上的循環(huán)僅發(fā)生在影響域區(qū)域內(nèi)及周圍。即使隨著模型的增大，內(nèi)部的串行循環(huán)也不會(huì)有所增加，此策略下的串行循環(huán)只和高斯點(diǎn)區(qū)域大小和模型節(jié)點(diǎn)間的間距有關(guān)。

圖2 不同高斯點(diǎn)對(duì)應(yīng)的影響域范圍和局部搜索示意圖Fig.2 The figures of influenle domains of differnet Gauss points and loal search(a)不同高斯點(diǎn)對(duì)應(yīng)的影響域范圍圖；(b)局部搜索示意圖

考慮到B(x)、A(x)、p(x)是有關(guān)模型節(jié)點(diǎn)上的值，因此，它們必須和相應(yīng)模型節(jié)點(diǎn)的坐標(biāo)一一對(duì)應(yīng)才會(huì)有意義。傳統(tǒng)的方式是把上述三個(gè)待求矩陣內(nèi)元素的坐標(biāo)按照位于整個(gè)速度模型上的絕對(duì)位置，予以編號(hào)，然后賦予坐標(biāo)值。這種方式有兩個(gè)缺點(diǎn)；①用數(shù)組來保存時(shí)，消耗較大內(nèi)存，并且極為稀疏；②參與下一步計(jì)算時(shí)，因?yàn)闃O為稀疏，所以耗時(shí)巨大。為了改變這一狀況，把上述待求量在速度模型上的絕對(duì)坐標(biāo)值，轉(zhuǎn)換成每個(gè)影響域內(nèi)的相對(duì)坐標(biāo)。

考慮到二維模型上剛度矩陣和質(zhì)量矩陣的物理含義，單獨(dú)針對(duì)一個(gè)模型節(jié)點(diǎn)而言是沒有意義的，而是應(yīng)該考慮各個(gè)高斯點(diǎn)影響域內(nèi)模型節(jié)點(diǎn)對(duì)模型節(jié)點(diǎn)的關(guān)系(圖2(a))。又因?yàn)镵、M均為稀疏對(duì)稱矩陣，因此我們只需要求取上三角或者下三角矩陣即可。而對(duì)于K、M的坐標(biāo)，則需要用基于整個(gè)速度模型的絕對(duì)坐標(biāo)表示。

當(dāng)求出每個(gè)高斯點(diǎn)影響域內(nèi)的K、M后，因?yàn)镃PU內(nèi)存的限制，無法一次性把整個(gè)模型上不同高斯點(diǎn)影響域內(nèi)包含相同模型節(jié)點(diǎn)處的值相加。從圖2(a)中我們可以發(fā)現(xiàn)，因?yàn)槊總€(gè)高斯點(diǎn)的影響域是有限的，而待求和處的模型節(jié)點(diǎn)(不同高斯點(diǎn)影響域共有的模型節(jié)點(diǎn))位于相應(yīng)高斯點(diǎn)的附近，所以當(dāng)把高斯點(diǎn)分塊時(shí)，相應(yīng)的待求和模型節(jié)點(diǎn)，也會(huì)被相應(yīng)的分開。為了克服CPU存儲(chǔ)空間的限制，分塊后的高斯點(diǎn)應(yīng)該滿足下面的條件：分塊后對(duì)應(yīng)的模型節(jié)點(diǎn)的絕對(duì)坐標(biāo)值(MS)滿足max(MS)- min(MS)≤M(M是和CPU存儲(chǔ)空間大小有關(guān)的整數(shù)，如CPU內(nèi)存為2 GB時(shí)，M=6 561)。從圖2(a)可知，鄰近高斯點(diǎn)影響域總會(huì)有重合區(qū)域，即速度模型中的模型節(jié)點(diǎn)會(huì)屬于鄰近不同的高斯點(diǎn)影響域。為了得到最終的剛度矩陣和質(zhì)量矩陣，需要把不同高斯點(diǎn)的影響域下包含的相同的模型節(jié)點(diǎn)的剛度矩陣和質(zhì)量矩陣求和。

4 稀疏矩陣壓縮存儲(chǔ)

稀疏矩陣可定義為非零元素占全部元素的百分比很小(如5%以下)的矩陣?；蛘哂械木仃嚪橇阍卣既吭氐陌俜直容^大(如近50%)，但它們的分布很有規(guī)律，利用這一特點(diǎn)可以避免存放零元素或避免對(duì)這些零元素進(jìn)行運(yùn)算。由于影響域通常是個(gè)較小區(qū)域(二維情形下包含十幾個(gè)節(jié)點(diǎn)即可滿足需求)，求出的K、M滿足稀疏矩陣的定義。

稀疏矩陣的壓縮方式有很多種，常用的是三元組方式，即用三個(gè)一維數(shù)組來表示稀疏矩陣，同時(shí)數(shù)組僅保存稀疏矩陣中的非零元素。這里采用按行壓縮存儲(chǔ)格式(CSR)。對(duì)于m×n(m行n列)具有q個(gè)非零元素的稀疏矩陣，應(yīng)用CSR 格式存儲(chǔ)時(shí)，使用的三個(gè)數(shù)組分別為：一個(gè)q×1維的值數(shù)組value，它按行序分成了m個(gè)段；一個(gè)q×1維的列下標(biāo)數(shù)組column；一個(gè)(m+1)×l 維的數(shù)組row，該數(shù)組中的元素指向各段中首元素在稀疏矩陣中的順序號(hào)。第i行非零元素的個(gè)數(shù)(row[i+1]-row[i])。

以文中的簡單三層模型為例：對(duì)于(596×300)*(596×300)的雙精度稀疏矩陣，按上三角或下三角保存的結(jié)果計(jì)算，保存在一維數(shù)組的元素個(gè)數(shù)小于10 000 000。如果采用一般方法全存儲(chǔ)所需要的空間為：(596×300)*(596×300)*sizeof(float)*16/(1 024*1 024*1 024)=1 905.52 GB。而采用CSR上三角或下三角存儲(chǔ)所需要的空間最大為10 000 000*sizeof(float)*16/(1 024*1 024*1 024)=0.6 GB。因此壓縮存儲(chǔ)所占全存儲(chǔ)的比例小于0.6/1 905.52=0.032%。

有關(guān)CSR格式的壓縮存儲(chǔ)情況參見表1。

表1 稀疏矩陣C的CSR格式壓縮存儲(chǔ)形式為(按行優(yōu)先)
Tab.1 The CSR format of the sparse matrixC(priority by row)

1-basedindexingvalue1-1-3-25464-4278-5row14691214column1231234513425

其中：value為稀疏矩陣中包含的非零元素值；row為每一行第一個(gè)非零元素在value中的起始偏移位置；column為非零元素值對(duì)應(yīng)的列標(biāo)。

1)在程序中，我們采用1-based設(shè)置 indexing。需要注意的是：最后一行最后一個(gè)非零元素在value的排列位置再加“1”，則為row中最后一個(gè)元素值。

2)實(shí)踐證明，選擇CSR格式較好，理由是CSR存儲(chǔ)格式，不僅相對(duì)來說占用較少內(nèi)存，同時(shí)對(duì)于后面涉及到的矩陣向量相乘，以及解線性方程組，效率更高。主要原因是矩陣向量相乘本身的運(yùn)算規(guī)則，使得CSR存儲(chǔ)格式更占優(yōu)勢(shì)。同時(shí)有關(guān)MKL的庫函數(shù)對(duì)于CSR格式的存儲(chǔ)，在做矩陣向量相乘時(shí)，包含有并行的思想。

5 基于OpenMP 對(duì)時(shí)間遞推的多核加速

疊前偏移可以作用于任意非層狀速度模型，疊前偏移主要是利用共反射點(diǎn)反射波的疊加。有關(guān)疊前偏移成像通常是分別計(jì)算單炮炮點(diǎn)正傳時(shí)的波場(chǎng)函數(shù)和反傳時(shí)檢點(diǎn)數(shù)據(jù)延拓的波場(chǎng)函數(shù)，然后利用成像條件，即可得到地下結(jié)構(gòu)的成像。這里將采取疊前偏移，同時(shí)利用互相關(guān)成條件。零延遲互相關(guān)成像條件的表達(dá)式為：

(20)

其中：tmax為最大記錄時(shí)間；S(x,t)為正傳的震源波場(chǎng)；R(x,tmax-t)為檢點(diǎn)數(shù)據(jù)反傳的接收?qǐng)?；I(x,t)為x處的成像結(jié)果。

將時(shí)間遞推公式(17)轉(zhuǎn)化成線性方程組可得式(21)。

(21)

M、K，待求解線性方程組的系數(shù)矩陣分別為M+((Δt)2/4)K、M。

在做矩陣向量相乘時(shí)，利用 Sparse BLAS 庫函數(shù)，其調(diào)用語句為：

Call mkl_dcsrsymv (uplo,m,a,ia,ja,x,y)

其中各個(gè)參數(shù)含義如表2所示。

在波場(chǎng)模擬和逆時(shí)偏移成像過程中，把原來需要對(duì)稀疏矩陣求逆的方程組轉(zhuǎn)化成線性方程組，并采用‘PARDISO’ 線性求解器求解?？紤]到在每一個(gè)時(shí)間迭代步中需要做三次矩陣向量相乘和求解兩次線性方程，即使K、M是壓縮存儲(chǔ)以后的，該過程仍需要消耗大量的時(shí)間。為了提高計(jì)算效率，我們利用已有的FDM程序，計(jì)算并均勻保存8個(gè)時(shí)間點(diǎn)處的波場(chǎng)值及其對(duì)時(shí)間的一階偏導(dǎo)數(shù)(考慮到服務(wù)器為8核CPU)。

表2 調(diào)用矩陣向量相乘各個(gè)參數(shù)的含義

Tab.2 The meaning of each parameter of calling matrix vector to multiply

輸入?yún)?shù)uplouplo=‘U’或‘u’表示為a只保存對(duì)稱矩陣的上三角;uplo=‘L’或‘l’表示為a只保存對(duì)稱矩陣的下三角;m整型。表示為稀疏方陣A的行或者列a雙精度型。表示為稀疏矩陣A按CSR格式壓縮存儲(chǔ)后的value向量ia整型。表示為稀疏矩陣A按CSR格式壓縮存儲(chǔ)后的row向量ja整型。表示為稀疏矩陣A按CSR格式壓縮存儲(chǔ)后的column向量x雙精度型。表示為已知被乘的向量,且包含m個(gè)元素輸出參數(shù)y雙精度型。表示為待求的向量,且包含m個(gè)元素

假設(shè)波場(chǎng)模擬的總時(shí)間為T，首先利用FDM程序，求出完整的波場(chǎng)正傳和反傳值，并保存波場(chǎng)正傳和反傳時(shí)T/8，2T/8，…，7T/8時(shí)刻的波場(chǎng)及一階偏導(dǎo)值；之后以此為基準(zhǔn)點(diǎn)值，即作為EFM中的起始值，這樣8個(gè)時(shí)間段內(nèi)的時(shí)間遞推即可在OpenMP程序中同時(shí)進(jìn)行。理論上，8個(gè)時(shí)間段內(nèi)的時(shí)間遞推同時(shí)進(jìn)行，可使效率提高到8倍左右，但實(shí)際上并非如此。主要原因是在利用Sparse BLAS 庫函數(shù)執(zhí)行矩陣向量相乘和‘PARDISO’求解線性方程組的過程中，所調(diào)用的函數(shù)本身是多核并行的，以Marmousi模型為例，在調(diào)用這兩個(gè)函數(shù)時(shí)，CPU利用率約為300%，同時(shí)內(nèi)存占用率約為20%(總內(nèi)存為30 G)。綜合考慮，我們一次并行的數(shù)量可選取4個(gè)時(shí)間段。這樣理論上應(yīng)該獲得8/3倍的加速，算例表明通過OpenMP，加速率約為2.5倍，與理論加速比基本吻合。

OpenMP編譯器命令以！$OMP parallel/!$OMP end parallel定義一個(gè)并行域[6]，并利用OpenMP函庫中的OMP_get_num_process( )函數(shù)，獲得CPU核的個(gè)數(shù)，根據(jù)前面討論的4段并行，確定線程數(shù)目(call OMP_set_num_threads(4))。并行部分程序流程見圖3，其中紅綠色箭頭表明計(jì)算過程有先后之分。

在編寫OpenMP并行程序時(shí)，特別需要注意在實(shí)現(xiàn)并行后，循環(huán)內(nèi)部所有變量根據(jù)具體情況，開辟為與循環(huán)維度和次數(shù)相關(guān)的數(shù)組。每個(gè)循環(huán)內(nèi)部的變量必須獨(dú)立于其他循環(huán)，只被該循環(huán)修改，不受其他任何循環(huán)的影響。對(duì)于各循環(huán)共同讀取但不修改的變量，可以不必定義(擴(kuò)充)為數(shù)組。調(diào)試并行程序中對(duì)那些非數(shù)組變量一定要小心謹(jǐn)慎，確保其在循環(huán)中始終為某一常量，否則就會(huì)出現(xiàn)錯(cuò)誤。并行后的程序需滿足以下原則：該循環(huán)內(nèi)部的所有變量，在任何情況下都不會(huì)被任何其他循環(huán)修改。

圖3 并行程序流程圖Fig.3 The flow diagram of parallelism

6 數(shù)值算例

6.1 水平層狀介質(zhì)模型

為了有效說明本文算法及改進(jìn)的合理性，先測(cè)試一個(gè)簡單模型：水平三層介質(zhì)模型(圖4)，該模型研究區(qū)域水平方向長為7.4 375 km，豎直方向深度為2.99 km，模型節(jié)點(diǎn)數(shù)為596×300。三層速度自上而下依次為2.5 km/s、3.5 km/s、4.5 km/s。采用3×3階次的高斯積分，高斯網(wǎng)格數(shù)為300×300。權(quán)函數(shù)為冪權(quán)函數(shù)，震源是主頻為25Hz的高斯子波，總記錄時(shí)間為3 s，時(shí)間間隔為1 ms。從圖5和圖6可以看出，得到利用FDM數(shù)值，基于OpenMP-Fortran編寫的并行EFM程序，對(duì)于三層簡單模型的RTM單炮和多炮疊后成像。圖7 表示基于FDM 20炮RTM疊加后的成像圖。由成像結(jié)果可知，改進(jìn)后EFM算法得到的成像精度與FDM成像精度是吻合的。

圖4 水平三層介質(zhì)速度模型Fig.4 The three-layer velocity model

圖5 單炮利用無單元法RTM成像結(jié)果Fig.5 The single-shot images obtained by the EFM (based on FDM and OpenMp) RTM

有關(guān)圖6、圖7分界面附近的干擾是因?yàn)椴捎秒p程波動(dòng)方程RTM引起的，圖6相比于圖7可發(fā)現(xiàn)無單元RTM結(jié)果的同相軸更細(xì)，主要原因是無單元的淺層噪音相關(guān)性較差，疊加干涉的效應(yīng)在橫向上的規(guī)律不明顯，客觀上避免了因強(qiáng)烈的干涉導(dǎo)致子波展寬的現(xiàn)象發(fā)生；而差分的淺層噪音相關(guān)性較好，疊加干涉后頂點(diǎn)能量加強(qiáng)并對(duì)反射界面主能量有展寬效應(yīng)。

圖6 利用無單元20炮RTM疊加后的成像結(jié)果Fig.6 The stacked image from 20 single-shot images obtained by the EFM (based on FDM and OpenMp) RTM

圖7 利用有限差分法20炮RTM疊加后的成像Fig.7 The stacked image from 20 single-shot images obtained by the finite difference RTM

6.2 Marmousi模型

圖8為Marmousi模型。該模型長為7.425 km，深度為2.99 km。把Marmousi模型離散為595×298個(gè)模型節(jié)點(diǎn)，仍然采用采用3×3階次的高斯積分，300×300高斯網(wǎng)格數(shù)，震源是主頻為25 Hz的高斯子波，總記錄時(shí)間為3 s，時(shí)間間隔為1 ms。利用已有FDM數(shù)值做EFM起始數(shù)值和OpenMP-Fortran編寫的并行程序。圖9 (a)、圖9 (b)、圖9 (c)、圖9 (d)、圖9 (e)分別表示基于OpenMP-Fortran并行的無單元算法在0.6 s、0.8 s、1.0 s、1.3 s、1.7 s的波場(chǎng)快照。圖10表示基于OpenMP-Fortran并行EFM 18炮RTM疊加后的成像圖。從圖10中可以看出，Marmousi模型中三個(gè)高角度逆掩斷層、斷層下面的鹽丘、以及深部兩側(cè)的高速體都得到很好的成像，結(jié)構(gòu)清晰，位置精確。圖11為基于FDM 18炮RTM疊加后的成像圖。從表3、圖10、圖11與圖12可知，和傳統(tǒng)無單元算法相比，在分步求取K、M時(shí)，對(duì)模型節(jié)點(diǎn)的內(nèi)層循環(huán)采用不同的搜索方式，對(duì)于效率的提升有明顯的影響；同時(shí)借助已有的FDM程序計(jì)算結(jié)果，通過OpenMP多核加速，所得成像結(jié)果是可信的，且精度在可接受范圍內(nèi)。

圖8 Marmousi速度模型Fig.8 Marmousi velocity model

圖9 無單位法(基于FDM數(shù)據(jù)做OpenMp多核加速)在炮點(diǎn)位于(100 m，3737.5m)對(duì)應(yīng)的0.6 s,0.8 s,1.0 s,1.3 s,1.7 s的波場(chǎng)快照Fig.9 Snapshots simulatcd by the EFM (based on FDM and OpenMp)in Marmousi model(a)0.6 s時(shí)的波場(chǎng)快照；(b)0.8 s時(shí)的波場(chǎng)快照；(c)1 s時(shí)的波場(chǎng)快照；(d)1.3 s時(shí)的波場(chǎng)快照；(e)1.7 s時(shí)的波場(chǎng)快照；

圖10 利用無單元法18炮RTM疊加后的成像Fig.10 The stacked image from 18 single-shot images obtained by the EFM (based on FDM and OpenMp)

圖11 利用有限差分法18炮RTM疊加后的成像Fig.11 The stacked image from 18 single-shot images obtained by the finite difference RTM

K,M計(jì)算時(shí)間/min 模擬時(shí)間/minRTM成像時(shí)間/min傳統(tǒng)EFM420520620改進(jìn)后EFM(全局搜索求K,M+OpenMp)5595135改進(jìn)后EFM(局部搜索求K,M+OpenMp)94989計(jì)算耗時(shí)比(傳統(tǒng)EFM:全局搜索求K,M+OpenMp)7.635.474.59計(jì)算耗時(shí)比(傳統(tǒng)EFM:局部搜索求K,M+OpenMp)46.6710.616.97

圖12 Marmousi模型RTM計(jì)算效率比較Fig.12 The comparison of computation efficiency of RTM in Marmousi

7 結(jié)論

研究的主要內(nèi)容為分步計(jì)算無單元算法中的剛度矩陣和質(zhì)量矩陣，并借助已有的有限差分程序的計(jì)算結(jié)果，利用OpenMP在時(shí)間遞推過程中多核加速。在地震波模擬和成像中的應(yīng)用，無單元算法一直存在的兩個(gè)問題：①計(jì)算機(jī)內(nèi)存的限制；②計(jì)算效率不高的問題。通過分布計(jì)算K、M解決內(nèi)存的瓶頸，同時(shí)提高了計(jì)算K、M的效率，并通過多核加速提升時(shí)間遞推過程的計(jì)算效率。

在分步計(jì)算K、M時(shí)，采用局部搜索內(nèi)層模型節(jié)點(diǎn)循環(huán)替代全局搜索，并把計(jì)算所得的K、M按行壓縮存儲(chǔ)(CSR)?？紤]到K、M的對(duì)稱性，只保存上三角(或者下三角)。同時(shí)需注意在求取的第二步，即求和的過程，要根據(jù)計(jì)算機(jī)的內(nèi)存，把高斯點(diǎn)相應(yīng)分塊處理去做求和，直到整個(gè)速度模型上所有高斯節(jié)點(diǎn)坐標(biāo)不同但模型節(jié)點(diǎn)坐標(biāo)相同的值求和完畢為止。轉(zhuǎn)化為線性方程組，避免了對(duì)超大稀疏矩陣求逆，通過‘PARDISO’線性求解器，可以有效解決問題。當(dāng)使用OpenMP 做基于CPU多核并行計(jì)算時(shí)，程序中靜態(tài)數(shù)組僅能支持小數(shù)組(比如一維數(shù)組最大為1000000)，當(dāng)數(shù)組較大時(shí)，程序中應(yīng)使用動(dòng)態(tài)數(shù)組。

數(shù)值算例表明，改進(jìn)后的無單元算法在滿足計(jì)算精度的前提下，計(jì)算效率有了很大地提升。這為計(jì)算更大的速度模型和三維速度模型提供了可能。同時(shí)文中的這些改進(jìn)措施同樣可用于有限元(FEM)問題的求解。當(dāng)然在計(jì)算效率方面，無單元法和有限差分法相比，還是存在很大差距，提升空間巨大。

關(guān)于未來的工作，我們可以有以下幾個(gè)方面的展望：

1)可研究類似小波變換等算法，對(duì)質(zhì)量矩陣和剛度矩陣進(jìn)行數(shù)據(jù)壓縮存儲(chǔ)，在可行的壓縮率和壓縮質(zhì)量下，可以有效地減少空間存儲(chǔ)和提高計(jì)算效率。

2)有關(guān)求取剛度矩陣和質(zhì)量矩陣可考慮移植到GPU平臺(tái)計(jì)算。

3)在二維無單元算法的基礎(chǔ)上，把無單元算法推廣到三維空間。

4)把時(shí)間域的無單元算法，推廣到頻率域。因?yàn)閱蝹€(gè)頻率間是相互獨(dú)立，頻率域無單元算法更加適合做并行計(jì)算。

[1]BELL N,GARLAND M.Efficient sparse matrix-vector multiplication on CUDA[R].Nvidia Technical Report NVR-2008-004,Nvidia Corporation,2008.

[2]BELYTSCHKO,T.,LU,Y.Y,et al.Element-free Galerkin methods:Int.J.Numer[J].Methods Eng，1994,37:229-256.

[3]CLAERBOUT,J.Toward a unified theory of reflector mapping[J].Geophysics,1911,36:467-481.

[4]Fan,Z,JIA,X.Element-Free Method and its Efficiency Improvement in Seismic Modeling and Reverse Time Migration[J].Geophys.Eng,2013,10(2):1-12.

[5]GOULD N I M,HU Y,SCOTT J A.A numerical evaluation of sparse direct solvers for the solution of large sparse,symmetric linear systems of equations Technical Report[R].Numerical Analysis Internal Report 2005-1,Ruther ford Appleton Laboratory,2005.

[6]HERMANNS M.Parallel Programming in Fortran 95 Using OpenMP[D].Universidad de Madrid,2002.

[7]JIA,X,HU,T.Element-free precise integration method and its applications in seismic modelling and imaging[J].Geophys,2006,166(1):349-372.

[8]MARFURT,K.J.Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[J].Geophysics,1984,49:533-549.

[9]KRYSL P,BELYTSCHKO T.Element-free Galerkin method:Convergence of the continuous and discontinuous shape[J].Computer methods in applied mechanics and engineering,1997,148(3):257-277.

[10]LU,Y.Y.,BELYTSCHKO,T.Tabbara,M.Element-free Galerkin method for wave propagation and dynamic fracture[J].Comput.Methods Appl.Mech.Engrg,1995,126:131-153.

[11]LU Y Y,BELYTSCHKO T,TABBARA M.Element-free Galerkin method for wave propagation and dynamic fracture Comput[J].Methods Appl.Mech.Eng,1995,126 :131-53.

[12]馮英杰,楊長春,吳萍.地震波有限差分模擬綜述[J].地球物理學(xué)進(jìn)展,2007,22(2):487-491.FENG Y J,YANG C C,WU P.The review of the finite-difference elastic wave motion modeling [J].Progress in Geophysics,2007,22(2):487-491.(In Chinese)

[13]紀(jì)國良,馮仰德.大規(guī)模有限元?jiǎng)偠染仃嚧鎯?chǔ)及其并行求解算法[J].數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用,2012,33(3):230-240.JI G L,FENG Y D.Parallel solving large-scale linear system of FEM based on compressed storage format [J].Journal on Numerical Methods and Computer Applications.2012,33(3):230-240.(In Chinese)

[14]賈曉峰,胡天躍,王潤秋.復(fù)雜介質(zhì)中地震波模擬的無單元法[J].石油地球物理勘探,2006,41(1):43-48.JIA X F,HU T Y,WANG R Q.A mesh free algorithm of seismic wave simulation in complex medium [J].Oil Geophysical Prospecting,2006,41(1):43-48.(In Chinese)

[15]閻樹文.無單元法用于高分辨率地震模擬與成像[D].合肥：中國科學(xué)技術(shù)大學(xué),2010.YAN S W.Element-free method applied in seismic modeling and imaging in high-resolution[D].Hefei：University of Science and Technology of China,2010.(In Chinese)

[16]楊頂輝.雙相各向異性介質(zhì)中彈性波方程的有限元解法及波場(chǎng)模擬[J].地球物理學(xué)報(bào),2002,45(4):575-583.YANG D H.Finite element method of the elastic wave equation and wavefield simulation in two-phase anisotropic media [J].Chinese J.Geophys,2002,45(4):575-583.(In Chinese)

[17]姚松,田紅旗.有限元?jiǎng)偠染仃嚨膲嚎s存貯及組集[J].中南大學(xué):自然科學(xué)版,2006,37(4):826-830.YAO S,TIAN H.Compressed storage scheme and assembly procedure of global stiffness matrix in finite dement analysis [J].Journal of Central South University:Natural Sciences,2006,37(4):826-830.(In Chinese)

[18]張湘?zhèn)?蔡永昌.一種改進(jìn)的無單元[J].計(jì)算力學(xué)學(xué)報(bào),2002,19(1):26-30.ZHANG X W,CAI Y C.An improved element free method [J].Chinese Journal of Computational Mechanics,2002,19(1):26-30.(In Chinese)

[19]周小平,周瑞忠.無單元法研究現(xiàn)狀及展望[J].工程力學(xué),2005,22(1):12-20.ZHOU X P,ZHOU R Z.Current Study and development of Element-Free Galerkin method [J].Engineering Mechanics,2005,22(1):12-20.(In Chinese)

OpenMP-accelerated element-free reverse-time migration

ZHOU Zhen,JIA Xiao-feng*

(School of Earth and Space Science,University of Science and Technology of China,Hefei 230026,China)

The wave-equation-based method was widely used in seismic modeling and imaging in the past years.Many numerical strategies such as the finite element method (FEM) and the finite difference method (FDM) are developed in solving seismic wave equations.However,both methods have some shortcomings of either accuracy or computational cost in practice.In fact,element-free method (EFM) has been applied in seismic modeling and seismic migration imaging.Compared with FEM and FDM,EFM has unique advantages:it doesn't need to be divided a large number of grid in advance and it satisfies local fitting because only the information of the nodes and the boundary of the study area are required in computation.However,there are many problems that EFM is applied in seismic modeling and seismic migration imaging in present,e.g.the problems of storage and computation efficiency about computing the coefficient matrix.The storage is caused by improper storage of sparse coefficient matrix.In this paper we compress the sparse matrices by the compressed sparse row (CSR) format and employ the following strategy:firstly,the value is computed at the corresponding model nodes within the influence domain of each Gauss point; secondly,the summation is performed within the influence domain of different Gauss points that contain the same model node.In the calculation of wave field time recursive,we solve the linear equations with the help of linear sparse solver 'PARDISO'.In order to further improve the computation efficiency,we use OpenMP basing the existing FDM data.Through the above methods,we can effectively solve these problems of memory limit and computational efficiency.The final results indicate that the above methods are accurate and efficient.

element-free method； coefficient matrix； linear system of equations； OpenMP； reverse time migration

2015-10-08 改回日期：2015-11-26

國家自然科學(xué)基金(41374006，41274117)

周震(1986-)，男，碩士，從事地震波偏移成像和并行加速研究，E-mail：zzhen12@mail.ustc.edu.cn。

*通信作者：賈曉峰(1979-)，男，副教授，主要從事地震波場(chǎng)的數(shù)值模似和偏移成像，Email:xjia@ustc.edu.cn。

1001-1749(2016)06-0821-11

P 631.4