周 震，賈曉峰

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

基于OpenMP加速的無單元逆時偏移成像

周 震，賈曉峰*

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

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

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

0 引言

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

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

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

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

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

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

1 無單元法基本理論介紹

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

(1)

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

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

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

(2)

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

(3)

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

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

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

(4)

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

(5)

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

(6)

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

為了求得a(x)，對式(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)可知所得近似場函數(shù)，僅與基函數(shù)、權(quán)函數(shù)以及鄰近點場值有關(guān)。

對方程(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)于時間的遞推公式：

(17)

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

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

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

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

(18)

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

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

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

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

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

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

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

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

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

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

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

4 稀疏矩陣壓縮存儲

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

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

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

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

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

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

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

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

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

5 基于OpenMP 對時間遞推的多核加速

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

(20)

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

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

(21)

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

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

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

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

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

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

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

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

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

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

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

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

6 數(shù)值算例

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

為了有效說明本文算法及改進(jìn)的合理性，先測試一個簡單模型：水平三層介質(zhì)模型(圖4)，該模型研究區(qū)域水平方向長為7.4 375 km，豎直方向深度為2.99 km，模型節(jié)點數(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的高斯子波，總記錄時間為3 s，時間間隔為1 ms。從圖5和圖6可以看出，得到利用FDM數(shù)值，基于OpenMP-Fortran編寫的并行EFM程序，對于三層簡單模型的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分界面附近的干擾是因為采用雙程波動方程RTM引起的，圖6相比于圖7可發(fā)現(xiàn)無單元RTM結(jié)果的同相軸更細(xì)，主要原因是無單元的淺層噪音相關(guān)性較差，疊加干涉的效應(yīng)在橫向上的規(guī)律不明顯，客觀上避免了因強烈的干涉導(dǎo)致子波展寬的現(xiàn)象發(fā)生；而差分的淺層噪音相關(guān)性較好，疊加干涉后頂點能量加強并對反射界面主能量有展寬效應(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個模型節(jié)點，仍然采用采用3×3階次的高斯積分，300×300高斯網(wǎng)格數(shù)，震源是主頻為25 Hz的高斯子波，總記錄時間為3 s，時間間隔為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的波場快照。圖10表示基于OpenMP-Fortran并行EFM 18炮RTM疊加后的成像圖。從圖10中可以看出，Marmousi模型中三個高角度逆掩斷層、斷層下面的鹽丘、以及深部兩側(cè)的高速體都得到很好的成像，結(jié)構(gòu)清晰，位置精確。圖11為基于FDM 18炮RTM疊加后的成像圖。從表3、圖10、圖11與圖12可知，和傳統(tǒng)無單元算法相比，在分步求取K、M時，對模型節(jié)點的內(nèi)層循環(huán)采用不同的搜索方式，對于效率的提升有明顯的影響；同時借助已有的FDM程序計算結(jié)果，通過OpenMP多核加速，所得成像結(jié)果是可信的，且精度在可接受范圍內(nèi)。

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

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

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

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

7 結(jié)論

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

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

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

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

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

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

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

4)把時間域的無單元算法，推廣到頻率域。因為單個頻率間是相互獨立，頻率域無單元算法更加適合做并行計算。

[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ī)模有限元剛度矩陣存儲及其并行求解算法[J].數(shù)值計算與計算機應(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ì)中彈性波方程的有限元解法及波場模擬[J].地球物理學(xué)報,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]姚松,田紅旗.有限元剛度矩陣的壓縮存貯及組集[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].計算力學(xué)學(xué)報,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-)，男，副教授，主要從事地震波場的數(shù)值模似和偏移成像，Email:xjia@ustc.edu.cn。

1001-1749(2016)06-0821-11

P 631.4