劉　慶，徐鳳姣，Osborne Kachaje

(長(zhǎng)江大學(xué) 油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢　430100)

20世紀(jì)50年代初Tikhonov(1950)[1]和Cagnird(1953)[2]分別提出利用天然電磁場(chǎng)進(jìn)行勘探的方法。

有限單元法由于精度高，能夠模擬復(fù)雜地形，因此在國(guó)內(nèi)外得到了廣泛應(yīng)用。徐世浙[3-4]對(duì)有限單元法進(jìn)行了深入研究，為國(guó)內(nèi)學(xué)者使用有限單元法進(jìn)行MT的正演奠定了理論基礎(chǔ)。

劉鵬茂[5]、顧關(guān)文[6]等基于CPU實(shí)現(xiàn)了大地電磁的二維與三維反演，取得了一定的加速效果。本文采用CPU(OpenMP)和GPU(CUDA)異構(gòu)并行的方式進(jìn)行MT的二維正演模擬，并采用全稀疏存儲(chǔ)的方式來(lái)存儲(chǔ)有限元中的剛度矩陣，有效地減少了單次MT二維正演所需的時(shí)間和對(duì)內(nèi)存的占用。

1　MT二維正演

根據(jù)麥克斯韋方程，角頻率為ω(時(shí)諧因子為e-iωt)的定態(tài)電磁場(chǎng)的方程是

▽×E=iwμH

(1)

▽×H=(σ-iωε)E

(2)

假定電性結(jié)構(gòu)是二維的，取走向?yàn)閥軸，x軸與y軸垂直，保持水平，z軸在垂直方向上?？芍?u/?y=0，將(1)、(2)展開(kāi)可得到兩個(gè)獨(dú)立的方程組，即E型波和H型波，方程組的具體形式如下。

E型：

(3)

H型：

(4)

本文選用H型波的方程，并基于有限單元法進(jìn)行大地電磁二維正演。其邊值問(wèn)題可以歸納為：

(5)

邊值問(wèn)題(5)與下列變分問(wèn)題等價(jià)：

(6)

2　有限單元法

利用4節(jié)點(diǎn)的等參單元來(lái)進(jìn)行數(shù)值模擬，采用雙線(xiàn)性插值，將式(6)中的區(qū)域積分分解為各單元積分之和：

(7)

單元積分：

(8)

(9)

單元積分：

(10)

(11)

單元積分：

(12)

(13)

(14)

δF(u)=δuTKu=0

(15)

由δu得任意性，可得Ku=0。解線(xiàn)性方程組前，將AB線(xiàn)上得邊界值代入。解得線(xiàn)性代數(shù)方程組后，得各節(jié)點(diǎn)得u。

3　稀疏存儲(chǔ)

本文采用三元組的存儲(chǔ)方式來(lái)儲(chǔ)存總體剛度矩陣K。并分別采用全稀疏存儲(chǔ)和不做稀疏處理的存儲(chǔ)對(duì)不同網(wǎng)格大小的二維區(qū)域進(jìn)行MT二維正演，所得的剛度矩陣所占用的內(nèi)存對(duì)比如圖1所示。

圖1　完全存儲(chǔ)與稀疏存儲(chǔ)對(duì)比

由圖1可知，采用完全存儲(chǔ)的方式時(shí)，剛度矩陣中的元素將隨著網(wǎng)格節(jié)點(diǎn)的個(gè)數(shù)以二次方的量級(jí)增長(zhǎng)，而全稀疏存儲(chǔ)的方式則對(duì)內(nèi)存的需求較小。

4　基于CPU和GPU的并行計(jì)算

本文采用CPU(OpenMP)和GPU(CUDA)異構(gòu)并行的方式來(lái)進(jìn)行MT二維正演計(jì)算，加速效果如表1所示。

表1　MT二維正演并行效率對(duì)比

5　MT二維并行正演結(jié)果分析

在二維情況下，設(shè)計(jì)了如圖2所示的二維正演模型，圖3為其二維正演視電阻率斷面圖。

圖2　二維正演模型

圖3　正演結(jié)果

從圖3中可以看出正演的視電阻率斷面圖可以很好的對(duì)正演模型做出響應(yīng)，進(jìn)一步地驗(yàn)證了MT二維并行正演程序的正確性。

6　結(jié)　語(yǔ)

1)采用全稀疏存儲(chǔ)的方式避免了由于存儲(chǔ)剛度矩陣中的零元素而造成的內(nèi)存的浪費(fèi)，為模擬范圍更大、網(wǎng)格更多、精度更高的MT二維正演提供了便利。

2)采用CPU和GPU異構(gòu)并行的方式較明顯地提高了單次大地電磁二維正演的速度。