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

周 震，賈曉峰

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

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

周 震，賈曉峰*

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

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

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

0 引言

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

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

自2006年起，EFM首次被用于地震波模擬和成像，這次嘗試為EFM在彈性波動方程中的應用奠定了基礎。其中包括：①權函數的選擇；②影響域的大小；③形函數的構造；④矩陣求逆的運算等。同時因為稀疏矩陣不恰當的存儲問題，該方法在2G計算機內存的存儲情況下，僅能處理81×81網格節點的模型，極大限制了EFM在較大速度模型中的應用，同時對于超大稀疏矩陣求逆，也是一項挑戰。

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

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

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

1 無單元法基本理論介紹

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

(1)

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

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

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

(2)

為了使用更加簡單的形式求解方程，我們取其級數形式的近似解uh(x)

(3)

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

當在二維情況下有：

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

(4)

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

(5)

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

(6)

其中：w(x-xI)為點x影響區域內的權函數其作用范圍在有限區域內[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)稱為形函數。

從方程(12)可知所得近似場函數，僅與基函數、權函數以及鄰近點場值有關。

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

(13)

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

(14)

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

(15)

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

(16)

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

(17)

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

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

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

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

(18)

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

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

3 分步求取質量矩陣和剛度矩陣

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

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

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

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

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

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

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

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

4 稀疏矩陣壓縮存儲

稀疏矩陣可定義為非零元素占全部元素的百分比很小(如5%以下)的矩陣。或者有的矩陣非零元素占全部元素的百分比較大(如近50%)，但它們的分布很有規律，利用這一特點可以避免存放零元素或避免對這些零元素進行運算。由于影響域通常是個較小區域(二維情形下包含十幾個節點即可滿足需求)，求出的K、M滿足稀疏矩陣的定義。

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

以文中的簡單三層模型為例：對于(596×300)*(596×300)的雙精度稀疏矩陣，按上三角或下三角保存的結果計算，保存在一維數組的元素個數小于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%。

有關CSR格式的壓縮存儲情況參見表1。

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

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

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

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

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

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

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

(20)

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

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

(21)

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

在做矩陣向量相乘時，利用 Sparse BLAS 庫函數，其調用語句為：

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

其中各個參數含義如表2所示。

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

表2 調用矩陣向量相乘各個參數的含義

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

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

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

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

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

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

6 數值算例

6.1 水平層狀介質模型

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

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

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

有關圖6、圖7分界面附近的干擾是因為采用雙程波動方程RTM引起的，圖6相比于圖7可發現無單元RTM結果的同相軸更細，主要原因是無單元的淺層噪音相關性較差，疊加干涉的效應在橫向上的規律不明顯，客觀上避免了因強烈的干涉導致子波展寬的現象發生；而差分的淺層噪音相關性較好，疊加干涉后頂點能量加強并對反射界面主能量有展寬效應。

圖6 利用無單元20炮RTM疊加后的成像結果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個模型節點，仍然采用采用3×3階次的高斯積分，300×300高斯網格數，震源是主頻為25 Hz的高斯子波，總記錄時間為3 s，時間間隔為1 ms。利用已有FDM數值做EFM起始數值和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模型中三個高角度逆掩斷層、斷層下面的鹽丘、以及深部兩側的高速體都得到很好的成像，結構清晰，位置精確。圖11為基于FDM 18炮RTM疊加后的成像圖。從表3、圖10、圖11與圖12可知，和傳統無單元算法相比，在分步求取K、M時，對模型節點的內層循環采用不同的搜索方式，對于效率的提升有明顯的影響；同時借助已有的FDM程序計算結果，通過OpenMP多核加速，所得成像結果是可信的，且精度在可接受范圍內。

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

圖9 無單位法(基于FDM數據做OpenMp多核加速)在炮點位于(100 m，3737.5m)對應的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傳統EFM420520620改進后EFM(全局搜索求K,M+OpenMp)5595135改進后EFM(局部搜索求K,M+OpenMp)94989計算耗時比(傳統EFM:全局搜索求K,M+OpenMp)7.635.474.59計算耗時比(傳統EFM:局部搜索求K,M+OpenMp)46.6710.616.97

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

7 結論

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

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

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

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

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

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

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

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].地球物理學進展,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].數值計算與計算機應用,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]賈曉峰,胡天躍,王潤秋.復雜介質中地震波模擬的無單元法[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].合肥：中國科學技術大學,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]楊頂輝.雙相各向異性介質中彈性波方程的有限元解法及波場模擬[J].地球物理學報,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].中南大學:自然科學版,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]張湘偉,蔡永昌.一種改進的無單元[J].計算力學學報,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]周小平,周瑞忠.無單元法研究現狀及展望[J].工程力學,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

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

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

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

1001-1749(2016)06-0821-11

P 631.4