白洪濤，李 昂，歐陽(yáng)丹彤，邢書(shū)豪，劉雪飛

（1.吉林大學(xué) 公共計(jì)算機(jī)教學(xué)與研究中心，長(zhǎng)春 130012；2.吉林大學(xué) 符號(hào)計(jì)算與知識(shí)工程教育部重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)春 130012；3.吉林大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130012；4.吉林大學(xué) 軟件學(xué)院，長(zhǎng)春 130012）

可控源音頻大地電磁法（controlled source audio－frequency magnetotelluric，CSAMT）是一種應(yīng)用廣泛、適用性強(qiáng)的大地電磁勘探法，目前已有許多研究成果.Coggon［1］提出了利用有限元法解決2.5維電磁問(wèn)題的思想；Stoyer等［2］利用有限差分法計(jì)算了二維介質(zhì)中垂直磁偶極子頻率域的2.5維電磁響應(yīng)，并成功地用于地下水資源勘查；Unsworth等［3］提出了模擬二維地球介質(zhì)且考慮源是三維性質(zhì)的2.5維有限元數(shù)值模擬系統(tǒng)，并討論了有限長(zhǎng)度電性源激發(fā)時(shí)場(chǎng)的特征；Torres－Verdín等［4］提出了2.5維數(shù)值模擬的快速算法；文獻(xiàn)［5］利用有限元法實(shí)現(xiàn)了CSAMT法的二維正演，并導(dǎo)出了波數(shù)域中電磁場(chǎng)第三類(lèi)邊值條件系數(shù)的公式；毛先進(jìn)等［6］利用格林函數(shù)導(dǎo)出了2.5維問(wèn)題的公式，以研究區(qū)域網(wǎng)格化中任意電性分布條件下空間電位所滿(mǎn)足的邊界積分方程，提高了成像質(zhì)量，擴(kuò)大了電阻率成像技術(shù)的實(shí)用范圍；底青云等［7－8］進(jìn)行了CSAMT2.5維有限元數(shù)值模擬.除傳統(tǒng)的數(shù)值解法外，文獻(xiàn)［9－11］探討了模擬退火、蟻群和粒子群等非線(xiàn)性方法在地球物理問(wèn)題反演中的應(yīng)用.隨著勘探范圍和精度的增加，傳統(tǒng)串行算法已不能滿(mǎn)足越來(lái)越復(fù)雜的勘探要求.Newman等［12］將并行計(jì)算技術(shù)運(yùn)用到三維電磁場(chǎng)反演中；劉羽等［13］提出了一種PC集群上的二維Occam反演方法，取得了一定性能的加速；Tan等［14－15］也將并行計(jì)算技術(shù)運(yùn)用到可控源音頻大地電磁和大地電磁正、反演中.

在多種并行方法中，OpenMP是可移植多線(xiàn)程應(yīng)用程序開(kāi)發(fā)的行業(yè)標(biāo)準(zhǔn)，在細(xì)粒度（循環(huán)級(jí)別）與粗粒度（函數(shù)級(jí)別）技術(shù)上都具有較高的效率.本文通過(guò)分析2.5維CSAMT串行正演算法的性能瓶頸，在串行算法的頻率域計(jì)算部分進(jìn)行多線(xiàn)程粗粒度的并行分解.OpenMP主線(xiàn)程將頻率各次循環(huán)分配給從線(xiàn)程和超線(xiàn)程，并進(jìn)行結(jié)果匯總處理.實(shí)驗(yàn)結(jié)果表明，通過(guò)合理的規(guī)劃并行粒度和任務(wù)分配，在保證并行計(jì)算精度的前提下，算法獲得了接近于CPU核心數(shù)目的線(xiàn)性加速比.

1 2.5維CSAMT正演

本文針對(duì)2.5維復(fù)電阻率電磁場(chǎng)正演算法，原理如下：首先通過(guò)有限元法，將待求解區(qū)域采用四邊形網(wǎng)格法剖分成多個(gè)小單元，在對(duì)起伏地表模擬上，采用地表節(jié)點(diǎn)高程值的不同實(shí)現(xiàn).單元剖分結(jié)束后，在其上建立電磁場(chǎng)的Maxwells方程組，以此推導(dǎo)出相應(yīng)的電磁場(chǎng)耦合頻率域微分方程.基于多個(gè)頻率，對(duì)推導(dǎo)出的方程采用雙二次插值法進(jìn)行離散化，得到有關(guān)各節(jié)點(diǎn)函數(shù)的線(xiàn)性代數(shù)方程組，再將復(fù)電阻率代入方程組中，利用基于不完全Cholesky分解的雙共軛梯度法求解電磁場(chǎng)響應(yīng)，即可得出正演的最終結(jié)果.算法描述如下：

1）讀入模型參數(shù)，包括發(fā)射線(xiàn)源和接收偶極、網(wǎng)絡(luò)剖分、觀測(cè)點(diǎn)、頻率和波數(shù)等；

2）離散化并生成節(jié)點(diǎn)間的關(guān)系；

3）提取當(dāng)前頻率；

4）在某波數(shù)下，生成有限元方程系數(shù)矩陣和右端項(xiàng)；

5）代入并生成邊界條件；

6）進(jìn)行大型稀疏線(xiàn)性有限元方程求解；

7）計(jì)算輔助場(chǎng)；若所有波數(shù)計(jì)算完成，則轉(zhuǎn)8）；否則轉(zhuǎn)4）；

8）反Fourier變換；

9）若所有頻率計(jì)算完成，則轉(zhuǎn)10）；否則轉(zhuǎn)3）；

10）輸出正演結(jié)果.

實(shí)際測(cè)試結(jié)果表明，串行算法的計(jì)算負(fù)載集中在算法步驟4）～8）部分.

2 并行算法

在OpenMP并行模式中，較常見(jiàn)的并行分解模式有兩種，一種是SPMD（single program multiple data），由于這種分解模式對(duì)于循環(huán)的并行較繁瑣，且對(duì)部分實(shí)例支持較差，因此本文采取另一種更適合的并行分解模式模式.

如圖1所示，F(xiàn)ork－Join模式中由編譯指導(dǎo)語(yǔ)句規(guī)定哪些代碼為并行區(qū)域，直至遇到并行區(qū)域，否則程序保持串行執(zhí)行.程序開(kāi)始時(shí)，僅有一個(gè)主線(xiàn)程執(zhí)行程序的串行部分，如變量初始化、文件的讀寫(xiě)操作和有限單元的分解等.在遇到編譯指導(dǎo)語(yǔ)句后，自動(dòng)進(jìn)入并行模塊，由Fork主線(xiàn)程派生出多個(gè)子線(xiàn)程共同完成并行模塊中的計(jì)算任務(wù).本文將頻率循環(huán)設(shè)置為并行模塊，將頻率循環(huán)任務(wù)分配給并行執(zhí)行的多個(gè)線(xiàn)程.計(jì)算完成后，這些派生線(xiàn)程會(huì)被主線(xiàn)程以Join形式注銷(xiāo)并收回資源，然后主線(xiàn)程繼續(xù)完成余下的程序，直至程序結(jié)束.

圖1 OpenMP并行的Fork－Join模式Fig.1 Fork－Join parallel－model of OpenMP

3 并行算法優(yōu)化

并行算法實(shí)現(xiàn)后，本文對(duì)其進(jìn)行進(jìn)一步優(yōu)化，以挖掘算法的并行潛能.

3.1 負(fù)載均衡

在多核多線(xiàn)程并行程序執(zhí)行過(guò)程中，會(huì)有同步點(diǎn)的存在，只有在每個(gè)線(xiàn)程都運(yùn)行到同步點(diǎn)時(shí)，程序才會(huì)向下執(zhí)行.而在實(shí)際計(jì)算中，很難保證各線(xiàn)程都在同一時(shí)間到達(dá)同步點(diǎn)，從而導(dǎo)致已經(jīng)到達(dá)同步點(diǎn)的線(xiàn)程等待尚未到達(dá)同步點(diǎn)的線(xiàn)程，產(chǎn)生空等時(shí)間，影響整個(gè)程序的運(yùn)行效率.若線(xiàn)程分配到的任務(wù)不均衡，這種現(xiàn)象就會(huì)頻繁出現(xiàn)，一些線(xiàn)程會(huì)長(zhǎng)期處于等待狀態(tài)，降低了CPU的工作效率.因此，在編寫(xiě)OpenMP代碼時(shí)，應(yīng)保證負(fù)載平衡.表1列出了120×28模型串行算法各頻率的計(jì)算時(shí)間.

表1 不同頻率的計(jì)算時(shí)間Table 1 Calculating time for different frequencies

由表1可見(jiàn)，頻率值不同對(duì)頻率循環(huán)的計(jì)算時(shí)間略有影響.因此，在任務(wù)分配上，本文采用OpenMP中的動(dòng)態(tài)調(diào)度策略，將循環(huán)迭代按給定CHUNKSIZE（默認(rèn)為1）劃分為若干個(gè)迭代塊.在塊內(nèi)部使用任務(wù)隊(duì)列按先來(lái)先服務(wù)的原則進(jìn)行計(jì)算.若一個(gè)線(xiàn)程完成了當(dāng)前計(jì)算的任務(wù)塊，系統(tǒng)會(huì)自動(dòng)分配給其一個(gè)新的任務(wù)塊，直到所有的任務(wù)塊都被執(zhí)行結(jié)束.最后一個(gè)迭代塊的大小會(huì)按實(shí)際情況平均分配，可能不等于指定的CHUNKSIZE.從而保證了各線(xiàn)程不會(huì)因同步等待而產(chǎn)生大量空等時(shí)間.

3.2 超線(xiàn)程

超線(xiàn)程技術(shù)（hyper－threading，HT）使得一個(gè)實(shí)體核上同時(shí)運(yùn)行兩個(gè)邏輯線(xiàn)程.這樣，在一個(gè)物理核心上，通過(guò)模擬兩個(gè)邏輯內(nèi)核而實(shí)現(xiàn)了多核技術(shù).超線(xiàn)程的工作原理：通常CPU在工作時(shí)執(zhí)行單元并未被充分利用，導(dǎo)致CPU的很多資源在某些時(shí)間閑置，為了充分調(diào)動(dòng)這些資源，HT技術(shù)模擬出另一個(gè)邏輯核心，在CPU執(zhí)行一個(gè)線(xiàn)程的空閑時(shí)間里，邏輯CPU則會(huì)利用這些空閑資源去執(zhí)行其他線(xiàn)程，從而極大提高了CPU的執(zhí)行效率.由外部可見(jiàn)，兩個(gè)邏輯線(xiàn)程是并發(fā)的.實(shí)際上，超線(xiàn)程并非真正的多線(xiàn)程，其性能也不能完全達(dá)到多線(xiàn)程的性能.

OpenMP本身對(duì)于超線(xiàn)程有很好的支持，為了盡量提高多核多線(xiàn)程的并行性能，本文采用超線(xiàn)程技術(shù)，在實(shí)驗(yàn)評(píng)測(cè)部分驗(yàn)證了超線(xiàn)程技術(shù)的效果.

4 實(shí)驗(yàn)分析

本文選用多個(gè)包括含有異常體和不含異常體的模型進(jìn)行測(cè)試，以更全面地驗(yàn)證算法的有效性.實(shí)驗(yàn)環(huán)境為CPU：I72670QM 4核8線(xiàn)程（超線(xiàn)程）；內(nèi)存：8Gb RAM；操作系統(tǒng)：64位 Windows7 Professional Version；編譯環(huán)境：Vistual Studio 2010＋I(xiàn)ntel Vistual Fortran 11.0.

4.1 計(jì)算精度

為了驗(yàn)證算法的數(shù)值準(zhǔn)確性，本文構(gòu)造一個(gè)含異常體的地電模型，模型參數(shù)為：網(wǎng)絡(luò)剖分120（x軸）×28（z軸），空氣層10層，大地18層，波數(shù)16個(gè)，頻率數(shù)16個(gè).初始條件為：空氣的電阻率值為107Ω·m，大地背景電阻率為100Ω·m.源中心位于地表x向剖分的第20塊單元上.異常體包含兩層：外層x軸上為－100～950m，z軸上為0～540m；內(nèi)層x軸上為300～600m，z軸上為50～350m.為了更直觀地對(duì)比結(jié)果，本文將并行和串行計(jì)算得出的結(jié)果分別畫(huà)出了等值線(xiàn)，如圖2所示.

圖2 并行（A），（C）和串行（B），（D）的正演計(jì)算對(duì)比結(jié)果Fig.2 Comparison of parallel（A），（C）and serial（B），（D）calculating results

由圖2可見(jiàn)，異常體分布周?chē)目醽喴曄辔缓鸵曤娮杪是€(xiàn)分布有較明顯的波動(dòng)，異常體所在位置及形狀較明顯.并行計(jì)算獲得了與串行計(jì)算一致的解釋結(jié)果.

4.2 性能分析

下面進(jìn)行多組模型測(cè)試以更準(zhǔn)確地測(cè)試并行算法的效率，并從加速比和并行效率兩方面討論并行算法的效率，不含異常體和含異常體模型的測(cè)試結(jié)果分別列于表2和表3.

表2 不含異常體模型的測(cè)試結(jié)果Table 2 Results based on models with no anomalies

表3 含異常體模型的測(cè)試結(jié)果Table 3 Results based on models with anomalies

由表2和表3可見(jiàn)：隨著實(shí)例規(guī)模的增長(zhǎng)，相應(yīng)的并行加速效果越來(lái)越明顯，同時(shí)并行效率也越高；模型的規(guī)模越大，OpenMP的加速效果越明顯.隨著核心數(shù)目的增加，由于線(xiàn)程調(diào)度開(kāi)銷(xiāo)和同步等原因，并行效率不斷降低，尤其在核心數(shù)目為7時(shí)，甚至?xí)a(chǎn)生性能倒退的現(xiàn)象，這是因?yàn)槿蝿?wù)數(shù)目不是7的整數(shù)倍，導(dǎo)致難以實(shí)現(xiàn)負(fù)載平衡.OpenMP中超線(xiàn)程的并行能力相當(dāng)突出，開(kāi)啟超線(xiàn)程最高能達(dá)到35%的并行效率提升，但其性能與實(shí)際物理核心性能還存在差距.

綜上所述，本文基于多核多線(xiàn)程技術(shù)，實(shí)現(xiàn)了基于OpenMP的2.5維電磁場(chǎng)頻率域并行正演算法.在保證計(jì)算結(jié)果正確可靠的前提下，通過(guò)合理的算法設(shè)計(jì)和任務(wù)調(diào)度獲得了較高的并行計(jì)算效率.

［1］Coggon J H.Electromagnetic and Electrical Modeling by the Finite－Element Method ［J］.Geophysics，1971，36（1）：132－155.

［2］Stoyer C H，Greenfield R J.Numerical Solutions of the Response of a Two－Dimensional Earth to an Oscillating Magnetic Dipole Source［J］.Geophysics，1976，41（3）：519－530.

［3］Unsworth M J，Travis B J，Chave A D.Electromagnetic Induction by a Finite Electric Dipole Source over a 2－D Earth［J］.Geophysics，1993，58（2）：198－214.

［4］Torres－Verdín C，Habashy T M.Rapid 2.5－Dimensional Forward Modeling and Inversion via a New Nonlinear Scattering Approximation［J］.Radio Science，1994，29（4）：1051－1079.

［5］LI Yuguo，Key K.2DMarine Controlled－Source Electromagnetic Modeling：Part 1.An Adaptive Finite－Element Algorithm ［J］.Geophysics，2007，72（2）：51－62.

［6］毛先進(jìn)，鮑光淑.2.5 維問(wèn)題電阻率正演的新方法 ［J］.中南工業(yè)大學(xué)學(xué)報(bào)，1997，28（4）：307－310.（MAO Xianjin，BAO Guangshu.A New Method for 2.5－Dimensional Resistivity Forward Modelling［J］.Journal of Central South University of Technology，1997，28（4）：307－310.）

［7］底青云，Unsworth M，王妙月.復(fù)雜介質(zhì)有限元法2.5維可控源音頻大地電磁法數(shù)值模擬 ［J］.地球物理學(xué)報(bào)，2004，47（4）：723－730.（DI Qingyun，Unsworth M，WANG Miaoyue.2.5－D CSAMT Modeling with the Finite Element Method over 2－D Complex Earth Media［J］.Chinese Journal of Geophysics，2004，47（4）：723－730.）

［8］底青云，Unsworth M，王妙月.有限元法2.5維 CSAMT數(shù)值模擬 ［J］.地球物理學(xué)進(jìn)展，2004，19（2）：317－324.（DI Qingyun，Unsworth M，WANG Miaoyue.2.5DCSAMT Modeling with Finite Element Method［J］.Progress in Geophysics，2004，19（2）：317－324.）

［9］師學(xué)明，王家映.地球物理資料非線(xiàn)性反演方法講座（三）：模擬退火法 ［J］.工程地球物理學(xué)報(bào)，2007，4（3）：165－174.（SHI Xueming，WANG Jiaying.Lecture on Non－linear Inverse Methods in Geophysics（3）：Simulated Annealing Method［J］.Chinese Journal of Engineering Geophysics，2007，4（3）：165－174.）

［10］王書(shū)明，劉玉蘭，王家映.地球物理資料非線(xiàn)性反演方法講座（九）：蟻群算法 ［J］.工程地球物理學(xué)報(bào)，2009，6（2）：131－136.（WANG Shuming，LIU Yulan，WANG Jiaying.Lecture on Non－linear Inverse Methods in Geophysics（9）：Ant Colony Optimization ［J］.Chinese Journal of Engineering Geophysics，2009，6（2）：131－136.）

［11］易遠(yuǎn)元，王家映.地球物理資料非線(xiàn)性反演方法講座（十）：粒子群反演方法 ［J］.工程地球物理學(xué)報(bào)，2009，6（4）：385－389.（YI Yuanyuan，WANG Jiaying.Lecture on Non－linear Inverse Methods in Geophysics（10）：Particle Swarm Optimization Inversion Method ［J］.Chinese Journal of Engineering Geophysics，2009，6（4）：385－389.）

［12］Newman G A，Alumbaugh D L.Three－Dimensional Massively Parallel Electromagnetic Inversion.Ⅰ.Theory［J］.Geophysical Journal International，2007，128（2）：345－354.

［13］劉羽，王家映，孟永良.基于PC機(jī)群的大地電磁Occam反演并行計(jì)算研究 ［J］.石油物探，2006，45（3）：311－315.（LIU Yu，WANG Jiaying，MENG Yongliang.PC Cluster Based Magnetotelluric 2－D Occam’s Inversion Parallel Calculation［J］.Geophysical Prospecting for Petroleum，2006，45（3）：311－315.）

［14］TAN Handong，TONG Tuo，LIN Changhong.The Parallel 3DMagnetotelluric Forward Modeling Algorithm［J］.Applied Geophysics，2006，3（4）：197－202.

［15］LIN Changhong，TAN Handong，TONG Tuo.Parallel Rapid Relaxation Inversion of 3DMagnetotelluric Data［J］.Applied Geophysics，2009，6（1）：77－83.